■于華良 潘軍 祝銘 劉興春 張魁珂

（吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院 吉林長(zhǎng)春130026）

基于DEM的火山口提取的研究

■于華良 潘軍 祝銘 劉興春 張魁珂

（吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院 吉林長(zhǎng)春130026）

數(shù)字高程模型即DEM，是地表形態(tài)高程屬性的數(shù)字化表達(dá)。本文基于二次高程異常的方法對(duì)DEM數(shù)據(jù)進(jìn)行火山地形提取的實(shí)驗(yàn)。DEM中火山和平原區(qū)域的處理以及DEM分辨率的大小都直接影響著火山口提取的質(zhì)量和效率。本文利用目視解譯、坡度分析重分類以及Laplacian算子運(yùn)算，并結(jié)合熔巖臺(tái)地地貌分類提取出火山口。對(duì)長(zhǎng)白山地區(qū)進(jìn)行試驗(yàn)，分析結(jié)果證明了此方法的可行性。

DEM火山口個(gè)數(shù)Laplacian模板分析

1　引言

DTM即數(shù)字地形模型，是基于數(shù)學(xué)模型建立起來(lái)的地表地形起伏模型。DTM是以高等數(shù)學(xué)和計(jì)算機(jī)等學(xué)科為依托，隨著這些學(xué)科和計(jì)算機(jī)成圖的發(fā)展成長(zhǎng)起來(lái)的。數(shù)字地形模型是將地形的真實(shí)的空間展布情況利用數(shù)字的形式表現(xiàn)出來(lái)。

數(shù)字高程模型DEM是構(gòu)成DTM的基礎(chǔ)，是區(qū)域地形表面簡(jiǎn)單的、離散的數(shù)字表示模型,它由規(guī)則水平間隔處地面點(diǎn)的抽樣高程矩陣組成，一般以柵格數(shù)據(jù)的形式表示。DEM作為數(shù)字化的地形圖，蘊(yùn)含著大量的、各種各樣的地形結(jié)構(gòu)和特征信息，是定量描述地貌結(jié)構(gòu)、水文過(guò)程、生物空間分布等空間變化的基礎(chǔ)數(shù)據(jù)。對(duì)于地貌的相關(guān)研究工作，極為重要的一種手段就是地形分析。利用遙感圖像配合DEM進(jìn)行地形地貌的解譯的方法，可以用來(lái)整合數(shù)據(jù)，提高遙感圖像和DEM數(shù)據(jù)所隱含的地貌信息的提取。本文研究目的是提取長(zhǎng)白山地區(qū)的火山口，將采用DEM數(shù)據(jù)的地形來(lái)提取。

在歸納遙感與DEM的特點(diǎn)運(yùn)用于地形地貌方面，程維明等學(xué)者利用SRTMDEM與遙感數(shù)據(jù)探究提取基本地貌類型的方法，基本地貌類型被由SRTMDEM派生的各種模型自動(dòng)提取出來(lái)，而且具有很大的實(shí)效性。利用高程模型及起伏度模型相結(jié)合的方式，初步實(shí)現(xiàn)了基本地貌類型的自動(dòng)提取，對(duì)區(qū)域地貌宏觀規(guī)律有較好的表達(dá)。2003年EndreD等利用遙感和DEM綜合解譯制作了區(qū)域土壤圖，2004年劉學(xué)軍等人對(duì)提取坡度、坡向的算法在DEM數(shù)據(jù)上進(jìn)行的試驗(yàn)、研究并對(duì)得到的精度進(jìn)行了分析,2009年周增波，程維明，周成虎等在長(zhǎng)白山靖宇縣地區(qū)運(yùn)用DEM數(shù)據(jù)進(jìn)行火山口地貌識(shí)別，并對(duì)提取和識(shí)別結(jié)果進(jìn)行探。2010年王博，羅微進(jìn)行了基于Aster G-DEM的海南島地形地貌提，龍恩，張洪巖，程維明等進(jìn)行了基于遙感與DEM的長(zhǎng)白山區(qū)數(shù)字地貌制圖方法的研。不同比例尺的DEM對(duì)提取地面坡度精度有較大影響，湯國(guó)安就這一問(wèn)題進(jìn)行研究并取得了一定結(jié)果。2003年劉澤慧提出利用DEM數(shù)據(jù)輔助山脊線、山谷線提取的方法。王培法等人利用不同比例尺和不同空間分辨率提取出了流域面積、河網(wǎng)密度、河道坡度、河道長(zhǎng)度、平均高程、平均坡度和河道長(zhǎng)度并做出了詳細(xì)的分。

本研究以DEM地形因子的劃分為主要技術(shù)手段對(duì)長(zhǎng)白山地區(qū)進(jìn)行研究，使用的是Landsat ETM+數(shù)據(jù)和30米分辨率DEM數(shù)據(jù)。采用基于二次微分的高程異常的方法，最后利用遙感圖像對(duì)長(zhǎng)白山地區(qū)的火山口個(gè)數(shù)進(jìn)行研究和探討。

2　提取方法

二次微分高程異常提取方法主要包括目視解譯、Laplacian模版分析運(yùn)算以及地形起伏度的計(jì)算。目視解譯需要統(tǒng)計(jì)出火山口的大致半徑。模板的如何選取是由火山口的形狀大小確定的，需要將研究區(qū)域范圍內(nèi)火山口的相關(guān)信息進(jìn)行統(tǒng)計(jì)。作為邊緣增強(qiáng)濾波，Laplacian濾波在運(yùn)用過(guò)程中可以忽略邊緣的方向。經(jīng)典的邊緣檢測(cè)算子的基本原理是考察圖像上每個(gè)像素的相應(yīng)鄰域內(nèi)的灰度值的變化情況，可以利用一階或二階方向?qū)?shù)變化規(guī)律，以此獲得邊緣信息。經(jīng)典的邊緣算子是通過(guò)各向同性的拉普拉斯算子演變出來(lái)的，最簡(jiǎn)單的各向同性微分算子是拉普拉斯算子。一個(gè)二維圖像函數(shù)f(x,y)的拉普拉斯算子定義為：

因?yàn)閷?duì)任意階微分都是進(jìn)行線性的操作，所以拉普拉斯變換也是一個(gè)線性算子。為了離散形式描述這一公式，可以用差分的形式來(lái)定義。在x方向上，

同理在y方向上，

由上述三個(gè)公式，可以推導(dǎo)出含兩個(gè)變量的離散Laplacian算子的形式如下

Laplacian濾波強(qiáng)調(diào)圖像中的最大值，它通過(guò)運(yùn)用一個(gè)具有最高中心值的變換核來(lái)完成，所有的Laplacian濾波卷積核的維數(shù)都必須是奇數(shù)。下面是3*3的模版卷積核及5*5的卷積核。由這種規(guī)律性可以推導(dǎo)出更高階的卷積核。

DEM地形起伏度可以在一定程度上反映出地貌的特征，地形起伏度是區(qū)域內(nèi)最大高程與最小高程之差。尋找這個(gè)固定區(qū)域空間單位是關(guān)鍵。在研究區(qū)內(nèi)以3*3的模版大小為起始依次增加到59*59來(lái)計(jì)算模版窗口的最大高差，再擬合成數(shù)學(xué)圖像找到最大高差所對(duì)的模版大小，確定空間單位的大小。

其數(shù)學(xué)模型為：

熔巖臺(tái)地在DEM上可以根據(jù)其地形起伏度來(lái)定義，長(zhǎng)白山地區(qū)熔巖臺(tái)地的地形起伏度在30m和70m之。

3　火山口提取實(shí)驗(yàn)及步驟

3.1研究區(qū)概況

本研究區(qū)位于吉林省東部長(zhǎng)白山地區(qū)內(nèi),其地理坐標(biāo)為東經(jīng)127°～129°、北緯41°～43°。研究區(qū)面積4008.00km2。以長(zhǎng)白山火山主峰為中心。

研究區(qū)主要地貌類型為火山地貌、構(gòu)造剝蝕地貌和流水地貌。以長(zhǎng)白山天池主峰為中心大致呈同心環(huán)狀劃分為巨型復(fù)式火山錐、山麓傾斜熔巖高原和熔巖臺(tái)地三大地貌單元。

在火山活動(dòng)時(shí)期，巖漿由火山口噴出。火山口的位置主要處于火山錐頂端，還有一部分位于其側(cè)面，當(dāng)巖漿不再噴發(fā)時(shí)，其流經(jīng)的通道中的巖漿冷卻，由此造成的環(huán)狀凹陷的區(qū)域即為火山口。

在DEM上火山口通常是類似于山頂點(diǎn)但是又不同于山頂點(diǎn)這主要取決于其形狀。長(zhǎng)白山地區(qū)的火山口形態(tài)各異，一部分火山口為四周高中間低的形態(tài)特征，另一部分則表現(xiàn)出中間高四周低的形態(tài)特征，但是在數(shù)學(xué)模擬圖像上其二階導(dǎo)數(shù)不等于零，也就是說(shuō)火山口的實(shí)際形態(tài)是上下起伏的，本文正是基于這一原理進(jìn)行提取。

具體實(shí)驗(yàn)步驟如下：

①對(duì)DEM進(jìn)行預(yù)處理后與遙感圖像匹配融合進(jìn)行目視解譯；

②利用坡度重分類后與原圖像相乘進(jìn)行加強(qiáng)處理；

③進(jìn)行Laplacian變換得到圖像；

④統(tǒng)計(jì)分析出火山口的值的范圍；

⑤通過(guò)波段運(yùn)算進(jìn)行進(jìn)一步提?。?/p>

⑥對(duì)DEM進(jìn)行地形起伏度的分析提取出熔巖臺(tái)地；

圖3-1　長(zhǎng)白山地區(qū)的ETM+數(shù)據(jù)7、4、2波段合成影像

⑦提取后圖像進(jìn)行疊加從而精確結(jié)果，最后得出結(jié)論。

3.2數(shù)據(jù)源及預(yù)處理

研究采用的遙感數(shù)據(jù)為1999年9月2日的長(zhǎng)白山地區(qū)ETM+遙感圖像，ETM+遙感數(shù)據(jù)有7個(gè)波段，其中第6波段是熱紅外波段，故將其去掉，保留其余的6個(gè)波段。圖3-1長(zhǎng)白山地區(qū)的ETM+數(shù)據(jù)7、4、2波段合成影像。

圖3-2　長(zhǎng)白山地區(qū)DEM預(yù)處理圖像

本研究采用的DEM數(shù)據(jù)為ASTER DEM數(shù)據(jù)分別是：41°N 127°E、41°N 128°E、42°N 127°E、42°N 128°E這四景圖像，然后進(jìn)行融合，鑲嵌邊界忽略零值。在運(yùn)算中，如果單一考慮坡度因素，那么會(huì)出現(xiàn)許多噪聲即零星的像元。為消除零星像元帶來(lái)的干擾，可以在生成坡度之前，要對(duì)DEM進(jìn)行預(yù)處理。如圖3-2為長(zhǎng)白山地區(qū)DEM預(yù)處理圖像。

圖3-3　目視解譯效果圖圖像

3.2.2目視解譯

在Arcgis中把landsat ETM+遙感圖像與處理后的DEM圖形進(jìn)行校正和配準(zhǔn)然后進(jìn)行裁剪然后進(jìn)行目視解譯，在Arcgis中新建一個(gè)矢量文件以表示火山口。統(tǒng)計(jì)出火山口的個(gè)數(shù)，統(tǒng)計(jì)出大致的大小。

統(tǒng)計(jì)柵格的大小從而來(lái)確定模版的選擇，把矢量數(shù)據(jù)轉(zhuǎn)化成柵格文件，改為30米的分辨率和圖像進(jìn)行配準(zhǔn)比較得到正確的統(tǒng)計(jì)值，經(jīng)統(tǒng)計(jì)分析后得出最大的涵蓋范圍從而確定選擇11*11的窗口進(jìn)行試驗(yàn)。

3.3坡度分析加強(qiáng)處理

圖3-4　坡度分析加強(qiáng)處理后圖像

由于要突出火山口的地形坡度，所以要先分析出區(qū)域中DEM的坡度然后進(jìn)行重分類，根據(jù)相關(guān)文獻(xiàn)可知最佳坡度為3.5°[2]。重分類后的圖像與原DEM相乘以突出平原和山地的差別，從而增強(qiáng)火山口的形狀使其更容易識(shí)別。

3.4運(yùn)用模版來(lái)提取

在envi下進(jìn)行l(wèi)aplacian變換處理圖像，通過(guò)試驗(yàn)后選擇11*11的窗口。處理后進(jìn)行訓(xùn)練樣本的統(tǒng)計(jì)分析，從中分析出處理后火山口的最大值、最小值、均值。進(jìn)行波段運(yùn)算（B1 GT 10352 AND B1 LT 15923603）*B1+(B1 GT 15923603 OR B1 LT 10352)*0后設(shè)立閾值提取火山口。興趣區(qū)統(tǒng)計(jì)量如下：

3.5用DEM地表起伏度因子來(lái)對(duì)的熔巖臺(tái)地進(jìn)行提取

在研究區(qū)內(nèi)根據(jù)相關(guān)文獻(xiàn)可知最佳的地形起伏度窗口選擇為19*19，從而確定空間單位的大小。對(duì)DEM圖像進(jìn)行重分類，進(jìn)而提取出熔巖臺(tái)地來(lái)。

3.6疊加分析出火山口個(gè)數(shù)

在Arcgis中對(duì)上面兩個(gè)結(jié)果圖層進(jìn)行疊加分析，然后分析提取出火山口，火山口的個(gè)數(shù)為184個(gè)。與目視解譯的結(jié)果進(jìn)行拓?fù)錂z查，確定提取后的點(diǎn)與原來(lái)目視解譯的火山口重疊個(gè)數(shù)為106。

圖3-5　處理以后的圖像

圖3-6　為熔巖臺(tái)地的提取圖

4　結(jié)果分析

本文在以長(zhǎng)白山地區(qū)為背景研究提取火山口的方法，通過(guò)目視解譯1999年的ETM+圖像可以大體知道火山口個(gè)數(shù)為132個(gè)。二次微分的高程異常的方法經(jīng)過(guò)Laplacian算子提取后與地貌分類相疊加提取火山口結(jié)果為184個(gè)，提取后與目視解譯的火山口進(jìn)行拓?fù)浞治?，得到重合的火山口?06個(gè)，其精度是在可以接受的范圍內(nèi)。某些分布于熔巖臺(tái)地邊沿或者分布在其他復(fù)雜的地理位置的火山口，是造成誤差的主要原因，另一方面，有些火山口的規(guī)模較小，無(wú)法被提取是造成誤差的另一原因。以上這些都給提取增加了難度，使得結(jié)果存在一定誤差，但是總體上來(lái)說(shuō)，此方法提取火山口的精度仍是較為準(zhǔn)確的。

圖3-7　疊加后提取火山口圖像

5　結(jié)論

對(duì)于某種地貌的提取，如火山口，根據(jù)其不同于其他地質(zhì)結(jié)構(gòu)的特點(diǎn)，火山口的實(shí)際形態(tài)是上下起伏的，在數(shù)學(xué)模擬圖像上其二階導(dǎo)數(shù)不為零；以及火山口在遙感和DEM圖像上的形態(tài)，最終確定相關(guān)參數(shù)的選取。利用模板提取和DEM地表起伏度因子提取熔巖臺(tái)地的結(jié)果相疊加，最終確定火山口的數(shù)目和位置。本文采用二次高程異常方法進(jìn)行提取，主要包括目視解譯、坡度分析重分類以及Laplacian算子運(yùn)算。在長(zhǎng)白山地區(qū)這種方法的應(yīng)用效果較好。證明二次高程異常方法具有較高的實(shí)用價(jià)值，在今后的工作中要繼續(xù)深入研究，將此方法的準(zhǔn)確性進(jìn)一步提高。

[1]劉學(xué)軍，王葉飛，等，基于DEM的坡度坡向誤差空間分布特征研究 [B]，測(cè)繪通報(bào)，2004.12:22-28.

[2]周增波，程維明，周成虎等，基于DEM的火山口地貌識(shí)別方法探討 [J]，地球信息科學(xué)報(bào)，2009.11(6):773-776.

[3]王博，羅微，基于Aster G-DEM的海南島地形地貌提取。碩士學(xué)位論文，海南大學(xué)，2010.

[4]龍恩，張洪巖，程維明，基于遙感與DEM的長(zhǎng)白山區(qū)數(shù)字地貌制圖方法研究，碩士學(xué)位論文，東北師范大學(xué)，2005.

[5]湯國(guó)安，楊勤科，不同比例尺DEM提取地面坡度的精度研究-以在黃土丘陵溝壑區(qū)的試驗(yàn)為例 [J]。水土保持通報(bào)，2001.01:54-55.

[6]劉澤慧，黃培之，DEM數(shù)據(jù)輔助的山脊線和山谷線提取方法的研究測(cè)繪科學(xué) [A]，山地學(xué)報(bào)，2003.12:34-36.

[7]湯國(guó)安，楊瑋瑩.一種基于DEM的明暗等高線制圖方法 [B]，測(cè)繪通報(bào)，2001.07: 39-40.

[8]丁峰.高志海.魏懷東.DEM在旱耕地荒漠化評(píng)價(jià)中的應(yīng)用-以武威、古浪二市 (縣)為例 [J].水土保持學(xué)報(bào)，2000.14(3):31-35.

[9]王培法，柵格DEM的尺度與水平分辨率對(duì)流域特征提取的分析-以黃土嶺流域?yàn)槔?[J]，江西師范大學(xué)學(xué)報(bào) (自然科學(xué)版)2004.06.

[10]張會(huì)平，張恒安，楊農(nóng)等，基于GIS的岷江上游地貌形態(tài)初步分析 [J]，中國(guó)地質(zhì)災(zāi)害與防治學(xué)報(bào)，2004.9:116-119.

[11]高楊，呂寧，薛重生，不同比例尺數(shù)字高程模型對(duì)土壤侵蝕強(qiáng)度分級(jí)的影響 [J]。中國(guó)水土保持，2007,(10):26-28.

[12]周啟鳴，劉學(xué)軍，數(shù)字地形分析，科學(xué)出版社，2008.

[13]基于遙感和DEM數(shù)據(jù)的柬埔寨洞里薩湖地區(qū)數(shù)字地貌制圖研究，碩士學(xué)位論文，南京師范大學(xué)，2008.

P208[文獻(xiàn)碼]B

1000-405X（2016）-2-405-3

于華良（1990～），男，碩士研究生，研究方向?yàn)檫b感解譯與制圖。

潘軍（1971～），男，副教授，研究方向?yàn)檫b感與地理信息系統(tǒng)教學(xué)和科研。

吉林大學(xué)研究生創(chuàng)新基金資助項(xiàng)目，項(xiàng)目號(hào)：2015137