基于DEM數(shù)據(jù)的地貌類型快速劃分系統(tǒng)研究

仲偉敬, 邢立新, 潘 軍, 王 婷, 王 凱, 張文哲

(1. 吉林大學(xué) 地球探測科學(xué)與技術(shù)學(xué)院, 長春 130026; 2. 西安衛(wèi)星測控中心 第一活動站, 陜西 渭南 714000)

0 引 言

地貌劃分方法很多, 但通常都很復(fù)雜, 且人為主觀影響較大。DEM(Digital Elevation Model)為地貌學(xué)的定量描述提供了豐富的數(shù)據(jù)支持, 通過DEM數(shù)據(jù)不僅可以獲得高程信息, 還可以派生出許多反應(yīng)地形地貌特征的地形因子。為了提高地貌提取方法的普適性, 近年來, 不少學(xué)者在以DEM數(shù)據(jù)為基礎(chǔ), 進(jìn)行地貌分類的方法上取得了一定的進(jìn)展。此方法以DEM數(shù)據(jù)為基礎(chǔ), 通過相關(guān)算法模型進(jìn)行一系列操作后, 選取最佳地形因子組合分別實現(xiàn)了對西南、 三峽等地區(qū)地貌類型的自動劃分, 分類結(jié)果較為理想[1-3]。該方法雖然技術(shù)相對成熟, 而且可以流程化操作, 但是在實際提取過程中： 部分地形因子不可直接提取, 而且中間步驟過多也使得操作顯得異常繁瑣。筆者在充分理解相關(guān)地形因子模型和相關(guān)算法的基礎(chǔ)上, 以ENVI(The Environment for Visualizing Images)二次開發(fā)平臺, 采用IDL(Interactive Data Language)語言, 將流程開發(fā)集成, 實現(xiàn)一鍵化尋找某一地區(qū)的最佳地形因子組合, 全(半)自動對地貌信息快速提取。

1 基于DEM的地貌類型快速提取原理

提取流程如圖1所示。微觀地形因子選取3×3窗口、 采用均值變點分析法分別獲得各宏觀地形因子的最佳窗口, 將各地形因子采用極差法進(jìn)行標(biāo)準(zhǔn)化處理, 然后將統(tǒng)一量綱后的地形因子獲得相關(guān)系數(shù)矩陣、 采用雪式熵值法剔除相關(guān)系數(shù)過大的地形因子, 進(jìn)而選取最佳地形因子組合。最后采用ENVI中的非監(jiān)督分類實現(xiàn)地貌信息快速提取。上述過程大致可以概括為最佳窗口下的地形因子獲取、 選取最佳地形因子組合和地貌類型自動劃分。

圖1 地貌類型快速提取流程圖Fig.1 Rapid extraction flow chart of landform types

1.1 相關(guān)地形因子概念以及算法

地形因子種類繁多, 筆者結(jié)合前人地貌分類經(jīng)驗, 選取高程、 坡度、 坡度變率、 全累積曲率、 地形起伏度、 高程變異系數(shù)、 地表切割深度和地表粗糙度8種地形因子[1,2]。其中高程、 坡度、 坡度變率為微觀地形因子, 剩余5個為宏觀地形因子。地形因子概念以及算法[4]如表1所示。

表1 地形因子概念以及算法

1.2 均值變點分析法確定宏觀地形因子最佳窗口

對于微觀地形因子, 普適的分析窗口為3×3窗口； 而對于宏觀地形因子, 選取最佳分析窗口可以有效提升分類效果, 針對不同區(qū)域, 如何確定最佳窗口是宏觀地形因子算法模型的關(guān)鍵。筆者選擇均值變點分析法確定宏觀地形因子的最佳窗口, 并從3×3窗口一直計算到33×33窗口, 步長為2, 以求取地形起伏度的最佳窗口為例, 計算步驟如下[5-7]。

1) 計算序列X

Xi=ln(ti/si)(1)

其中序列X為{Xi,i=3,5,7,…,33}；ti為各窗口下的地形起伏度均值(m)；si為各窗口面積大小(m2)；i為窗口大小。

(2)

(4)

4) 計算S與Si的差值ΔSi

ΔSi=S-Si(5)

找出ΔSi的最大值, 即為拐點, 對應(yīng)的i值即為最佳窗口大小。

1.3 極差法統(tǒng)一量綱

由于各地形因子的值域不同, 在實際地貌分類時往往會造成權(quán)重不同, 為了消除這種影響, 需要對其統(tǒng)一量綱。筆者采用極差法進(jìn)行標(biāo)準(zhǔn)化處理, 算法如下[6]

=255(xij-xmin)/(xmax-xmin)(6)

1.4 雪式熵值法選取最佳地形因子組合

在實際地貌分類過程中, 地形因子維度過多可能會造成分類混淆或者分類不理想, 往往只選取幾個彼此隔離程度好、 又最具有代表性的地形因子組合參與分類。在求取各地形因子相關(guān)系數(shù)矩陣后, 某幾個地形因子彼此間相關(guān)系數(shù)較大, 這時采取雪式熵值法確定各地形因子, 算法如下

S=n/2+n/2ln(2n)+1/2ln|Ms|(7)

其中n為選取的因子數(shù)；Ms為n×n的方差-協(xié)方差矩陣；S為n維數(shù)據(jù)子集的熵： 行列式值|Ms|越大,S值越大,n維數(shù)據(jù)子集含的信息越多, 這時n維因子即為最佳地形因子組合。

2 基于DEM的地貌類型快速提取設(shè)計與實現(xiàn)

圖2 系統(tǒng)流程以及關(guān)鍵步驟圖Fig.2 System processes and key step diagrams

2.1 系統(tǒng)設(shè)計

按照圖1中地貌類型快速提取流程, 系統(tǒng)流程以及關(guān)鍵算法設(shè)計如下(見圖2)： 1)輸入DEM影像后, 通過均值變點分析法確定宏觀地形因子的最佳窗口； 2) 根據(jù)各地形因子的概念和算法按照最佳窗口大小分別提取最佳窗口下的微觀和宏觀地形因子； 3) 通過極差法對提取結(jié)果進(jìn)行統(tǒng)一量綱； 4) 取得各地形因子間相關(guān)系數(shù)矩陣后, 根據(jù)雪式熵值法確定最佳地形因子組合； 5) 最后將最佳窗口下而且統(tǒng)一量綱后的8個地形因子、 相關(guān)系數(shù)矩陣和最佳地形因子組合輸出。

2.2 系統(tǒng)實現(xiàn)

系統(tǒng)以ENVI二次開發(fā)平臺為基礎(chǔ), 采用IDL語言進(jìn)行開發(fā), 調(diào)用IDL函數(shù)包、 ENVI二次開發(fā)函數(shù)接口, 在函數(shù)里設(shè)置相應(yīng)的參數(shù), 根據(jù)相關(guān)算法實現(xiàn)[8-12], 關(guān)鍵技術(shù)如下。

1) DEM影像讀取。通過讀取文件地址來讀取影像和投影信息, 部分代碼如下。

dem_image=Read_tiff(dem_address,GeoTiff=geoInfor);讀取影像存儲到dem_image

dem_geo=ENVI_GET_PROJECTION(FID=dem_fid,PIXEL_SIZE=pixel_size);dem_geo存儲影像投影信息

2) 均值變點分析法確定宏觀地形因子最佳窗口, 根據(jù)算法公式, 以地形起伏度為例, 部分代碼如下。

tr_Xi=fltarr(16);tr_Xi用來存放序列Xi

FORi=1,16,1 DO BEGIN

tr_Xi[i-1]=ALOG(tr_mean_0[i-1]/(3+(i-1)*2)/(3+(i-1)*2)/8100);先計算各窗口下地表切割深度均值與各窗口面積的比值、 然后將其結(jié)果取對數(shù)

ENDFOR

tr_X_average=mean(tr_Xi);tr_X_average用來存放序列Xi的算術(shù)平均值

tr_X1_average=mean(tr_Xi[0:7]);tr_X1_average用來存放第1段樣本的算術(shù)平均值

tr_X2_average=mean(tr_Xi[8:15]);tr_X2_average用來存放第2段樣本的算術(shù)平均值

tr_Si_DS=fltarr(16);tr_Si_DS用來存放序列Xi與總體平均值的離差平方,DS是離差平方的縮寫

tr_X1_DS=fltarr(16);tr_X1_DS用來存放序列Xi與第1段平均值的離差平方

tr_X2_DS=fltarr(16);tr_X2_DS用來存放序列Xi與第2段平均值的離差平方

FOR i=1, 16 DO BEGIN

tr_Si_DS[i-1]=(tr_Xi[i-1]-tr_X_average)*(tr_Xi[i-1]-tr_X_average)

tr_X1_DS[i-1]=(tr_Xi[i-1]-tr_X1_average)*(tr_Xi[i-1]-tr_X1_average)

tr_X2_DS[i-1]=(tr_Xi[i-1]-tr_X2_average)*(tr_Xi[i-1]-tr_X2_average)

ENDFOR

tr_S=total(tr_Si_DS);tr_S用來存放序列Xi的離差平方和

tr_Si_DV=fltarr(15);tr_Si_DV用來存儲tr_S與各窗口的離差平方和的差值

FOR i=1,15 DO BEGIN

tr_Si_DV[i-1]=tr_S-total(tr_X1_DS[0:i-1])-total(tr_X2_DS[i:15])

ENDFOR

tr_Si_DV_max=max(tr_Si_DV,index);在S與Si的差值中找最大值點并且找到它的索引值

tr_bestwindow_size=(index+1)*2+1;通過索引值來計算最佳窗口大小

3) 地形因子提取。微觀地形因子提?。?對于坡度、 坡度變率以及全累積曲率, 根據(jù)算法其中間結(jié)果可以在ENVI軟件的地形模型模塊中直接提取。本系統(tǒng)在實現(xiàn)時需要將該功能抽調(diào)出使用, 方法是直接調(diào)用ENVI二次開發(fā)提供的函數(shù)包ENVI_DOIT, 此函數(shù)包已經(jīng)把涉及到的算法集成到一個函數(shù)里, 只需根據(jù)需要設(shè)置幾個參數(shù)即可, 以坡度為例, 提取代碼如下。

ENVI_DOIT,‘TOPO_DOIT’,BPTR=[0],DIMS=dims,FID=dem_fid,IN_MEMORY=1,pos=[0],R_FID=slope_fid,PIXEL_SIZE=pixel_size,KERNEL_SIZE=3

slope_image=ENVI_GET_DATA(FID=slope_fid,DIMS=dims,pos=[0])

宏觀地形因子提取。對于地形起伏度、 地表切割深度以及高程變異系數(shù), 根據(jù)算法首先需要計算得到中間結(jié)果： 窗口下高程的最大值、 最小值、 均值和標(biāo)準(zhǔn)差。然后根據(jù)相關(guān)公式再一次計算, 將計算結(jié)果填充到窗口中間的格網(wǎng), 然后在整幅影像下進(jìn)行移動、 循環(huán)分別填充到窗口中心。以3×3窗口為例, 如圖3所示。

以提取地形起伏度為例, 部分代碼如下。

tr_image=fltarr(length,width)

n1=tr_bestwindow_size;n1用來存儲地形起伏度的最佳窗口

FORi=1,width-n1+1,1 DO BEGIN;外側(cè)行循環(huán)

FORj=1,length-n1+1,1 DO BEGIN;里側(cè)列循環(huán), 從第1列開始循環(huán), 步長為1, 設(shè)置窗口寬度為n, 當(dāng)窗口右側(cè)頂?shù)秸跋褡钔鈧?cè)停止, 停止位置為length-n1+1

window_max=max(dem_image[j-1:j+n1-2,i-1:i+n1-2]); 算窗口內(nèi)的最大值

window_min=min(dem_image[j-1:j+n1-2,i-1:i+n1-2]); 算窗口內(nèi)的最小值

window_tr=window_max-window_min;算窗口內(nèi)地形起伏度： 在一個特定區(qū)域內(nèi), 最高點與最低點高程的差值

tr_image[j+(n1-1)/2-1,i+(n1-1)/2-1]=window_tr;窗口中心點坐標(biāo)存進(jìn)地形起伏度數(shù)組

ENDFOR

ENDFOR

圖3 窗口移動示意圖Fig.3 Window move diagram

宏觀地形因子： 地表粗糙度提取, 先求取不同窗口下的坡度, 然后以坡度為基礎(chǔ)根據(jù)公式再求取相應(yīng)窗口下的地表粗糙度, 部分代碼如下。

ENVI_DOIT,‘TOPO_DOIT’,BPTR=[0],DIMS=dims,FID=dem_fid,IN_MEMORY=1,pos=[0],R_FID=slope_fid,PIXEL_SIZE=pixel_size,KERNEL_SIZE=33

slope_33=ENVI_GET_DATA(FID=slope_fid,DIMS=dims,pos=[0])

surface_roughness_33=1/cos(slope_33*!DTOR)

surface_roughness_mean_0[(33-3)/2]=mean(surface_roughness_33)

4) 極差法統(tǒng)一量綱, 以高程因子為例, 部分代碼如下。

dem_min=min(dem)

dem_max=max(dem)

dem_gui=fltarr(after_n);dem_gui存儲歸一化后的數(shù)據(jù)

FORi=1,after_n,1 DO BEGIN

dem_gui[i-1]=(dem[i-1]-dem_min)*255/(dem_max-dem_min)

ENDFOR

5) 相關(guān)系數(shù)矩陣, 雪式熵值法確定最佳地形因子組合直接調(diào)用IDL函數(shù)即可, 求矩陣相關(guān)系數(shù)函數(shù)“CORRELATE”、 求方差-協(xié)方差矩陣函數(shù)“IMSL_COVARIANCES”、 求矩陣行列式值函數(shù)“DETERM”, 部分代碼如下。

r[i-1,j-1]=CORRELATE(p,q);r為系數(shù)

i_c=IMSL_COVARIANCES(all_n1);算方差-協(xié)方差矩陣

zhi=DETERM(i_c);算行列式值

6) 結(jié)果輸出, 以輸出坡度變率和不同地形因子組合以及其行列式值為例, 部分代碼如下。

write_tiff,‘e: empwork坡度變率3x3窗口.tif’,sv_gui_zhuan,GeoTiff=geoInfor,/FLOAT

openw,lun,‘E: empwork從大到小排列后的地形因子組合以及對應(yīng)行列式值.txt’,/append,/Get_lun

printf,lun,lian

free_lun,lun

3 應(yīng)用實例

3.1 研究區(qū)以及數(shù)據(jù)源

長白山位于吉林省東南部, 地處中朝兩國的邊境線上, 最高峰白頭峰2 700 m左右, 研究區(qū)以長白山天池景區(qū)、 望天鵝景區(qū)為中心, 地貌整體形態(tài)呈起伏狀、 高度逐漸向四周遞減, 不同位置相對高程跨度大, 是提取基本地貌類型典型的試驗區(qū)。筆者數(shù)據(jù)源選取SRTMDEM數(shù)字高程數(shù)據(jù)產(chǎn)品, 分辨率為90 m×90 m, 精度相當(dāng)于1 ∶25萬DEM, 數(shù)據(jù)下載自地理空間數(shù)據(jù)云, 其已經(jīng)經(jīng)過一定的數(shù)據(jù)加工, 根據(jù)研究區(qū)覆蓋范圍以及國家線為界限, 將影像裁剪, 地理范圍為41°～43°N、 127°～129°E。

3.2 系統(tǒng)快速提取研究區(qū)內(nèi)最佳地形因子組合

將裁剪的DEM影像輸入系統(tǒng)后, 不用像以往經(jīng)過很多步驟、 繁瑣操作后得出相應(yīng)結(jié)果, 本系統(tǒng)可一鍵化先后輸出最佳窗口大小下而且歸一化后的8個地形因子、 相關(guān)系數(shù)矩陣以及根據(jù)行列式值從大到小排列后的最佳地形因子組合。

1) 宏觀地形因子最佳窗口確定：以提取地形起伏度的最佳窗口為例, 如表2所示。在窗口為15×15時ΔSi值最大7.81, 此點即為拐點, 所以地形起伏度的最佳窗口選取為15×15, 面積相當(dāng)于1.82 km2。其他宏觀地形因子最佳窗口提取結(jié)果為： 高程變異系數(shù)15×15, 地表切割深度15×15, 地表粗糙度13×13。

表2 平均地形起伏度與窗口面積對應(yīng)關(guān)系

2) 相關(guān)系數(shù)矩陣提取結(jié)果如表3所示。可看到該研究區(qū)內(nèi)坡度、 地形起伏度、 高程變異系數(shù)、 地表切割深度和地表粗糙度這5個地形因子彼此間相關(guān)性強, 相關(guān)系數(shù)均大于0.6。

表3 地形因子間的相關(guān)系數(shù)矩陣

3) 將歸一化后最佳窗口下提取的8個地形因子都輸入系統(tǒng)后, 地形因子組合數(shù)選取為5, 按照從大到小排列后的方差-協(xié)方差行列式值輸出結(jié)果如表4所示。

表4 不同因子組合的行列式值(已從大到小排列)

其中12 357組合, 行列式值最大, 為最佳地形因子組合； 從序號[1]、[2]、[3]、[4]可看到組合結(jié)果分別為12 357、12 367、12 356、12 358： 除1、2、3一樣外, 在5、6、7、8中選出兩個即可, 側(cè)面印證了2)的提取結(jié)果： 即該研究區(qū)內(nèi)坡度、 地形起伏度、 高程變異系數(shù)、 地表切割深度和地表粗糙度這5個地形因子彼此間相關(guān)性較強； 從序號[54]、[55]、[56]可以看到： 行列式值較小的結(jié)果中都有4這個地形因子, 也就是全累積曲率對于該研究區(qū)地貌劃分最不適用。

綜上所述, 選取高程、 坡度、 坡度變率、 地形起伏度和地表切割深度這5個地形因子用于本研究區(qū)地貌劃分。

3.3 根據(jù)最佳地形因子組合對地貌類型全(半)自動劃分

起伏度與海拔是區(qū)分基本地貌類型的兩個重要參數(shù), 長白山最高峰不超過3 500 m, 所以研究區(qū)內(nèi)涉及到的基本地貌類型如表5所示[13,14]。

表5 研究區(qū)內(nèi)基本地貌類型

將系統(tǒng)提取的最佳地形因子組合在ENVI軟件下先后進(jìn)行波段疊合、 非監(jiān)督分類(ISODAT: Iterative Self-Organizing Data Analysis Technique)和分類后處理即可得到研究區(qū)內(nèi)地貌類型劃分圖： 1) 波段疊合后如圖4所示, 分別以高程、 坡度和地形起伏度作為RGB 3波段假彩色顯示； 2) 在非監(jiān)督分類過程中, 根據(jù)長白山區(qū)域地貌特征, 將初始分類數(shù)目設(shè)置為16, 迭代次數(shù)設(shè)置為15, 循環(huán)收斂閾值設(shè)置為0.95； 3) 初始分類后逐步進(jìn)行類別合并/定義、 更改分類顏色、 聚類、 小斑塊去除： 在處理過程中會存在圖斑(孤島)現(xiàn)象, 由于生成結(jié)果圖比例尺為1 ∶20萬, 為了統(tǒng)計分析的需要, 在這里圖斑過濾窗口取50×50, 得到分類結(jié)果如圖5所示。

圖4 假彩色合成地形因子疊加圖 圖5 長白山地區(qū)地貌類型分布圖 Fig.4 False color composite topographic Fig.5 Distribution map of landform types in factor overlay changbaimountain area

3.4 結(jié)果分析

前人雖然對長白山地區(qū)研究較多, 但尚缺少一個公認(rèn)的地貌劃分結(jié)果, 只有1 ∶50萬吉林省地貌成圖僅供參考。1) 從圖4、 圖5中可見, 研究區(qū)覆蓋范圍內(nèi)由左上角逐漸到影像中心海拔逐漸上升； 2) 從圖5中可見, 海拔臺地占據(jù)了絕大部分, 此區(qū)域海拔相對較高、 起伏度變化較小, 這與長白山海拔絕大部分處于1 000～1 500 m之間、 坡度一般在10°左右基本吻合； 3) 從圖5中也可看出, 低海拔臺地分布較少, 這與研究區(qū)內(nèi)平原幾乎沒有也基本吻合； 4) 從圖4中可以看到, A點(天池景區(qū))、 B點(望天鵝景區(qū))變化梯度較大, 這是海拔迅速升高的表現(xiàn)： A點海拔高度在2 200 m, B點海拔在2 050 m左右, 區(qū)域海拔呈現(xiàn)快速上升狀態(tài), 這在圖5中分類結(jié)果也有體現(xiàn)。

綜上所述, 分類結(jié)果相對可靠。

3.5 微觀地貌提取

在實際地貌提取過程中, 地貌分類的結(jié)果跟數(shù)據(jù)分辨率有很大關(guān)系。筆者又分別選取天池景區(qū)和望天鵝景區(qū)30 m×30 m DEM數(shù)據(jù)： 將兩個相對較小的研究區(qū)DEM數(shù)據(jù)分別輸入系統(tǒng)后, 得到的最佳地形因子組合均為： 高程3×3、 坡度3×3、 坡度變率3×3、 高程變異系數(shù)15×15和地表切割深度15×15； 然后將5個地形因子疊合后將研究區(qū)分成4類, 得出分類結(jié)果, 如圖6所示。

a 天池景區(qū)實際地物 b 天池景區(qū)地貌分類效果圖 c 望天鵝景區(qū)實際地物 d 望天鵝景區(qū)地貌分類效果圖圖6 實際地物與微觀地貌劃分對比圖Fig.6 Comparison of actual and microscopic features

圖6a、 圖6c分別是在Google Earth查看的實際地貌： 天池和望天鵝景區(qū)； 分類結(jié)果圖6b、 圖6d均未進(jìn)行任何分類后處理, 圖中存在許多斑點, 這里如若對斑點進(jìn)行聚類、 去除, 其分類結(jié)果也是一定程度上的假結(jié)果, 所以為便于對結(jié)果更好的進(jìn)行分析, 這里不對分類結(jié)果進(jìn)行任何后續(xù)處理。

圖6b、 圖6d中沒有對分類結(jié)果進(jìn)行嚴(yán)格的定義類別, 以圖6a、 圖6c地圖為參考, 無論是山脈、 溝谷還是相對平坦的苔地在圖6b、 圖6d中都可以清晰的分辨； 但從圖6a中也可看出, 天池中心湖水與天池景區(qū)四周的苔地化為一類, 所以該方法還需要加入一定的人工影像判別。

4 結(jié) 語

利用已有的算法和模型, 借用ENVI二次開發(fā)平臺、 采用強大的IDL語言編程實現(xiàn), 可將流程集成化, 使研究人員從繁瑣步驟中解脫。筆者結(jié)合前人研究成果, 開發(fā)了基于DEM數(shù)據(jù)的最佳地形因子組合地貌類型快速劃分系統(tǒng), 經(jīng)不斷測試, 算法和模型均編寫正確, 系統(tǒng)無BUG, 達(dá)到了初始的設(shè)計要求； 并以長白山作為試驗區(qū), 以提取的研究區(qū)內(nèi)最佳地形因子組合進(jìn)行地貌分類, 宏觀地貌和微觀地貌提取結(jié)果均較為理想。

由微觀地貌提取結(jié)果, 可以看出該方法還存在一定的誤判, 需加入適當(dāng)?shù)娜斯び跋衽凶x；筆者在對均值變點分析法編程實現(xiàn)過程中, 窗口只計算到33×33, 如果按照1 ∶100萬比例尺大范圍內(nèi)地貌制圖, 此時宏觀地形因子的最佳窗口選取還有待考究, 本論文編程步驟僅供參考。