基于UAV高分辨率DEM的復(fù)雜微地貌形態(tài)特征分析
——以恐龍谷南緣山區(qū)為例

羅為東,甘 淑?,袁希平,高 莎,胡 琳,袁新悅

(1.昆明理工大學(xué) 國(guó)土資源工程學(xué)院，650000,昆明；2.云南省高校高原山地空間信息測(cè)繪技術(shù)應(yīng)用工程研究中心，650000,昆明；3.滇西應(yīng)用技術(shù)大學(xué)地球科學(xué)與工程技術(shù)學(xué)院，671000,云南大理)

近年來(lái)，許多國(guó)內(nèi)外學(xué)者結(jié)合不同分辨率的DEM數(shù)據(jù)提取地形特征要素，并且根據(jù)組合不同地形因子實(shí)現(xiàn)對(duì)地形特征綜合分析，促進(jìn)高分辨率DEM數(shù)據(jù)在土地適宜性分析和水土流失分析領(lǐng)域的發(fā)展。地形因子是土壤侵蝕模型中的重要變量：我國(guó)的黃土高原地區(qū)，當(dāng)陡坡長(zhǎng)度在10～60 m范圍內(nèi)，坡長(zhǎng)與土壤流失之間的關(guān)系與USLE土壤侵蝕經(jīng)驗(yàn)?zāi)Ｐ途哂泻芎玫亟菩?，而與RUSLE土壤侵蝕經(jīng)驗(yàn)?zāi)Ｐ徒Y(jié)果則不太理想[1]。采用數(shù)字地形分析方法提取山區(qū)坡譜、坡長(zhǎng)和坡度因子，探討3種因子的關(guān)系：坡譜信息熵變化范圍反映地形起伏由平緩到強(qiáng)烈的變化趨勢(shì)；坡長(zhǎng)和坡度因子平均值介于2.72～18.61之間，在一定程度上表征土壤流失量的大小[2]。

此外，DEM分辨率大小影響著土壤侵蝕模型應(yīng)用中的地形因子參數(shù)精度。DENG等[3]通過(guò)對(duì)DEM重采樣的方式，研究地形因子對(duì)DEM分辨率依賴(lài)性，得出不同地形因子受DEM分辨率變化的影響不一致，各地形因子間的相關(guān)性與DEM分辨率呈線性關(guān)系，DEM分辨率降低即地形因子間的相關(guān)性降低。李蒙蒙等[4]為研究DEM分辨率對(duì)地形因子(LS)提取精度的影響，將30 m分辨率DEM數(shù)據(jù)進(jìn)行重采樣為7個(gè)不同分辨率的DEM，結(jié)果表明隨著DEM分辨率的降低，坡長(zhǎng)和坡度的相關(guān)系數(shù)呈現(xiàn)降低的趨勢(shì)。對(duì)于構(gòu)建土壤流失模型，DEM分辨率對(duì)LS因子值有重要的影響：在云貴高原區(qū)不同等級(jí)侵蝕面積隨DEM分辨率發(fā)生變化，DEM分辨率由5 m降低至50 m，所計(jì)算微度侵蝕面積增加49.2%，中度侵蝕面積減小92.5%，輕度增加13.7%[5]。該地區(qū)中度土壤侵蝕主要發(fā)生在海拔較高區(qū)域，同時(shí)海拔較高區(qū)域也是地形起伏度較大區(qū)域，土壤侵蝕分布與地形起伏度分布具有一致性。綜合國(guó)內(nèi)外研究現(xiàn)狀不難發(fā)現(xiàn)，土壤侵蝕與地形有重要關(guān)系，土壤侵蝕模數(shù)會(huì)隨著DEM分辨率降低而減小，高精度DEM更為真實(shí)描述微地貌形態(tài)特征，從而更精確反映土壤侵蝕模數(shù)。

對(duì)于TCI模型，SHAN等[6]用坡度因子和曲率因子來(lái)變化定義地形復(fù)雜性，結(jié)果表明地形復(fù)雜度依賴(lài)于DEM精度。王雷等[7]提取的地形TCI模型，有效通過(guò)DEM描述黃土丘陵溝壑區(qū)的地表形態(tài)變化。盧華興等[8]利用多因子評(píng)價(jià)方法選取 4 種局部地形因子并分析融合4種因子最終得到每個(gè)格網(wǎng)的地形TCI指標(biāo)。

筆者以云南祿豐恐龍國(guó)家地質(zhì)公園南緣山區(qū)為研究區(qū)，基于0.5 m高分辨率DEM，遵循地形因子有效性、易計(jì)算性和因子相互獨(dú)立性原則，選取坡度、地表粗糙度、標(biāo)準(zhǔn)曲率、地表切割度等單地形因子進(jìn)行加權(quán)綜合分析，對(duì)地形因子通過(guò)組合地形復(fù)雜度(Terrain Complexity Index,TCI)模型來(lái)描述和表達(dá)地表形態(tài)，以地形表面褶皺和破碎程度的指標(biāo)判斷土壤侵蝕區(qū)域。

1 研究區(qū)概況

研究區(qū)域位于云南祿豐恐龍國(guó)家地質(zhì)公園南緣山區(qū)，隸屬云南楚雄彝族自治州,地理坐標(biāo)為E 101°38′06″～102°24′34″，N 24°51′33″～25°30′45″。測(cè)區(qū)范圍內(nèi)地貌類(lèi)型復(fù)雜多樣，以構(gòu)造侵蝕地貌、方山地貌為主，地勢(shì)東北高，西南低，最高海拔2 200 m，最低海拔1 302 m。測(cè)區(qū)內(nèi)存在小型中生代紅色沉積盆地，成土母巖由中生代紫色砂頁(yè)巖和元古代的碳酸鹽巖交錯(cuò)分布。受母巖影響，土壤類(lèi)型呈帶分布，紫色土、紅壤分布廣泛。測(cè)區(qū)位置如圖1示所示，其中圖1(a)來(lái)源于百度地圖，圖1(b)為實(shí)地踏勘手機(jī)拍攝。該區(qū)域?qū)賮啛釒Ъ撅L(fēng)性氣候，干濕分明，且土壤沙化嚴(yán)重，裸土、裸巖分布廣泛。

圖1 實(shí)驗(yàn)區(qū)概況Fig.1 Overview of the experimental area

2 材料與方法

2.1 基于UAV的數(shù)據(jù)獲取及DEM構(gòu)建精度分析

試驗(yàn)區(qū)影像數(shù)據(jù)通過(guò)DJI Phantom 4 RTK進(jìn)行影像數(shù)據(jù)采集。首先根據(jù)測(cè)區(qū)實(shí)際周邊地理環(huán)境，交通狀況，結(jié)合Google earth平臺(tái)對(duì)測(cè)區(qū)進(jìn)行航線規(guī)劃，其次通過(guò)前期航線規(guī)劃以及后期數(shù)據(jù)質(zhì)量要求，選擇合適的參數(shù)設(shè)置。鑒于本次數(shù)據(jù)獲取主要面向地形特征分析應(yīng)用，試驗(yàn)區(qū)選取恐龍谷南緣環(huán)狀山區(qū)東北方向的山包，面積0.28 km2。本次飛行參數(shù)設(shè)置為：航向和旁向重疊度均設(shè)80%，平均飛行高度150 m，航測(cè)天氣條件良好，共采集392張影像，影像平均分辨率為0.07 m。

使用Context Capture Master軟件通過(guò)集群方式對(duì)實(shí)測(cè)數(shù)據(jù)進(jìn)行影像匹配、自動(dòng)空三等步驟構(gòu)建密集影像匹配點(diǎn)云。點(diǎn)云需通過(guò)濾波處理，將地物點(diǎn)(植被、建筑物、信號(hào)塔、高大獨(dú)立地物等)剔除構(gòu)建DEM。具體步驟如下：1)應(yīng)用TerraSolid 軟件的Terrascan模塊中漸進(jìn)加密不規(guī)則三角網(wǎng)(Progressive TIN Densification: PTD)濾波算法，設(shè)迭代角6°，迭代距離1.4 m的參數(shù)進(jìn)行濾波處理，得到地面點(diǎn)數(shù)據(jù)集；2)對(duì)于DEM的構(gòu)建，首先由地面點(diǎn)數(shù)據(jù)構(gòu)建不規(guī)則三角網(wǎng)，其次選擇反距離權(quán)重法的像元分配類(lèi)型和線性函數(shù)的填充方法，通過(guò)自然領(lǐng)域插值法計(jì)算像元值，采樣距離選定為0.5、1.0、1.5和2.0，以此構(gòu)建分辨率為0.5、1.0、1.5及2.0 m的DEM。

為驗(yàn)證DEM精度，試驗(yàn)首先在野外用RTK采集30個(gè)點(diǎn)數(shù)據(jù)，點(diǎn)位分布疊加到正射影像顯示，如圖2示，再分別將0.5、1.0、1.5和2.0 m分辨率DEM高程值賦30個(gè)同名點(diǎn)。

圖2 實(shí)測(cè)點(diǎn)位Fig.2 Location of actual measuring points

依據(jù)重復(fù)測(cè)量原理，取觀測(cè)值的算數(shù)平均值作為真值[9-10]，即RTK采集野外30個(gè)實(shí)測(cè)點(diǎn)為真值，0.5、1.0、1.5和2.0 m的DEM中選30個(gè)同名點(diǎn)的高程值為觀測(cè)值，分別按式(1)、(2)、(3)和(4)計(jì)算相關(guān)的誤差值，結(jié)果保留4位小數(shù)(表1)。

di=(Hi-Hture)；

(1)

(2)

(3)

(4)

式中：d為高程誤差，m；M為平均誤差，m，可平均體現(xiàn)數(shù)據(jù)集中的誤差；SD為標(biāo)準(zhǔn)偏差，m，表示高程誤差數(shù)據(jù)離散程度；R為均方根誤差，m，以突顯對(duì)比值與真值之間離散程度[11]；i為4種分辨率序號(hào)(i=0.5,1.0,1.5,2.0)。

表1 分析DEM精度指標(biāo)計(jì)算結(jié)果

通過(guò)表1計(jì)算結(jié)果結(jié)合表2比例尺精度綜合分析表明，除2 m分辨率DEM超出7 mm外，其余三者均滿足國(guó)家1∶500的DEM精度標(biāo)準(zhǔn)?？紤]高程值離散程度和與真值誤差離散程度，1.5 m分辨率DEM的SD=0.04 m和R=0.040 8 m，與其相比1.0 m分辨率DEM的SD和R略低，分別為0.036 6和0.037 1 m。但0.5 m分辨率DEM的SD和R分別僅為0.016 7和0.017 3，所以后續(xù)實(shí)驗(yàn)均以0.5 m分辨率DEM為基礎(chǔ)進(jìn)行。

表2 比例尺精度

2.2 TCI模型研究方法

TCI模型用于描述和表達(dá)地表形態(tài)的復(fù)雜程度，是評(píng)價(jià)地形表面褶皺和破碎程度的指標(biāo)參數(shù)，描述地形復(fù)雜度基于統(tǒng)計(jì)方法有：高程標(biāo)準(zhǔn)差、地形起伏度、等高線密度、溝壑密度等地形因子，基于幾何方法描述因子有：形態(tài)復(fù)雜指數(shù)、地表粗糙度、曲率等因子。

試驗(yàn)區(qū)以構(gòu)造侵蝕地貌為主，構(gòu)建TCI模型主要選取以下地形因子[12-13]：1)其土壤侵蝕量的大小受地面坡度制約，坡度(Slope)和土壤侵蝕度成正比，受重力侵蝕地區(qū)一般坡度值高，所以坡度因子是反映地形復(fù)雜度的典型因子；2)地表粗糙度(Surface Roughness)在微地形地貌中把地形表面崎嶇程度定為粗糙度，粗糙度客觀反映地形表面抗風(fēng)蝕的能力，具體計(jì)算如式(5)；3)地表切割深度(Surface Cut Depth)[14]用合適地面點(diǎn)領(lǐng)域范圍的高程均值與相同范圍內(nèi)高程最小值相減，按式(6)計(jì)算；4)曲率主要包括剖面曲率、平面曲率、全曲率、標(biāo)準(zhǔn)曲率等。曲率反映地形扭曲變化的程度，其中剖面曲率和平面曲率分別沿最大坡度的方向、垂直于最大坡度的方向描述地形的曲面形態(tài)，而標(biāo)準(zhǔn)曲率(Standard Curvature)[15-16]是結(jié)合剖面和平面曲率，更準(zhǔn)確詮釋地形變化程度，如圖3所示，來(lái)源于ArcGIS官方網(wǎng)址對(duì)曲率函數(shù)的解釋。

SR=1/cos(Slopeπ/180)；

(5)

SCD=Dmean-Dmin。

(6)

式中：SR為地表粗糙度；Slope為坡度；SCD為地表切割深度；Dmean和Dmin分別為DEM的均值和最小值；量綱均為1。

源于ArcGIS官方網(wǎng)站https:∥www.esri.com/en-us/arcgis/products/arcgis-pro/resources From the official ArcGIS website https:∥www.esri.com/en-us/arcgis/products/arcgis-pro/resources圖3 標(biāo)準(zhǔn)曲率圖Fig.3 Standard curvature diagram

以0.5 m的高分辨率DEM為基礎(chǔ)，基于ArcGIS平臺(tái)，采用因子分析方法對(duì)試驗(yàn)區(qū)進(jìn)行地形復(fù)雜度分析。利用統(tǒng)計(jì)和幾何綜合的方法描述地形復(fù)雜度指標(biāo)，針對(duì)試驗(yàn)區(qū)的特性，選取反映地形本質(zhì)特征的坡度、地表粗糙度、標(biāo)準(zhǔn)曲率、地表切割深度等地形因子作為評(píng)價(jià)因子，其中坡度用來(lái)描述試驗(yàn)區(qū)關(guān)于曲面的傾斜程度，地形粗糙度用來(lái)表示每個(gè)地形表面單元中地表崎嶇的復(fù)雜程度。通過(guò)變異系數(shù)法，如式(7)用于確定各項(xiàng)因子權(quán)重，通過(guò)各地形因子的均值和標(biāo)準(zhǔn)差計(jì)算變異系數(shù)，單因子變異系數(shù)比所有因子變異系數(shù)總和，以此確定單因子權(quán)重。最后，應(yīng)用ArcGIS求取TCI指數(shù)。

(7)

3 單因子指標(biāo)結(jié)果分析

3.1 坡度形態(tài)特征分析

基于ArcGIS平臺(tái)的坡度分析，如表3所示。平緩坡、緩坡和斜坡面積所占比例較為均勻，在20%左右。平緩坡面積比例最大，為25.10%，急陡坡面積比例最小，僅為3.76%。

表3 坡度數(shù)據(jù)分析

3.2 地表粗糙度形態(tài)特征分析

針對(duì)試驗(yàn)區(qū)裸土、裸巖分布廣泛的特性，在ArcGIS中使用空間分析，基于坡度提取粗糙度，通過(guò)高程值按低、中、高進(jìn)行分級(jí)[17]，如表4所示。低粗糙度比例最大，為74.03%。

表4 地表粗糙度分析

3.3 標(biāo)準(zhǔn)曲率形態(tài)特征分析

通過(guò)ArcGIS輸出每個(gè)像元的表面曲率，曲率是表面的二階偏導(dǎo)數(shù)，該值是通過(guò)將選中像元與8個(gè)相鄰像元擬合而得。依據(jù)表5分析結(jié)果，標(biāo)準(zhǔn)曲率>0面積比例為37.10%，<0面積比例達(dá)到67.18%，而表面平滑地區(qū)面積比例僅有4.82%。

表5 標(biāo)準(zhǔn)曲率分析

3.4 地表切割深度形態(tài)特征分析

對(duì)于地表切割深度，分析最佳窗口是關(guān)鍵步驟[18]。應(yīng)用 ArcGIS的鄰域分析，對(duì)實(shí)驗(yàn)區(qū)進(jìn)行從3×3窗口到41×41窗口的20組數(shù)據(jù)分析，獲得試驗(yàn)區(qū)地表切割深度和分析窗口的關(guān)系如圖4所示。最大高差法輔助人工作圖判斷，地表切割度值隨窗口大小在拐點(diǎn)5×5—19×19呈現(xiàn)上升趨勢(shì)，除2個(gè)端點(diǎn)區(qū)間內(nèi)整體數(shù)值高于擬合函數(shù)。在19×19—35×35窗口區(qū)間整體數(shù)值低于擬合函數(shù)，在35×35拐點(diǎn)后又迅速遞增。

圖4 領(lǐng)域窗口和地表切割深度分析Fig.4 Domain window and surface cut depth analysis

根據(jù)表6的統(tǒng)計(jì)數(shù)據(jù)，結(jié)合統(tǒng)計(jì)軟件可知地表切割深度與網(wǎng)格單元大小與對(duì)數(shù)方程擬合得到

Y=A0lnx+A1。

(8)

式中：Y為每個(gè)網(wǎng)格的平均地勢(shì)起伏度；x為每個(gè)網(wǎng)格的面積；A0和A1為參數(shù)。

表6 地表切割深度分析

擬合曲線為Y=6.015 2lnx+11.238，R2=0.93。在5×5、19×19窗口和擬合曲線出現(xiàn)交點(diǎn)，綜合考慮試驗(yàn)區(qū)面積較小和DEM精度高2個(gè)因素，較大窗口會(huì)造成DEM模糊，故選取5×5拐點(diǎn)作為最佳窗口并按照Ⅰ、Ⅱ、Ⅲ、Ⅳ等級(jí)對(duì)該地區(qū)地表切割深度進(jìn)行劃分(表6)。地表切割深度Ⅰ、Ⅱ級(jí)所占比例相對(duì)更多，分別占總面積的32.54%和31.79%。

4 TCI運(yùn)用結(jié)果分析與驗(yàn)證討論

4.1 TCI構(gòu)建結(jié)果分析

TCI模型進(jìn)行定量化存在模糊性[19]，對(duì)此，筆者對(duì)4項(xiàng)地形因子先按自然間斷點(diǎn)法分級(jí)，再進(jìn)行手動(dòng)調(diào)整，最后以2為步長(zhǎng)按低、較低、中、較高和高復(fù)雜度分為5級(jí)[12,20]，并對(duì)其進(jìn)行賦值。對(duì)4項(xiàng)地形因子統(tǒng)一色帶，添加注記，分級(jí)結(jié)果如圖5和表7所示。

圖5 地形復(fù)雜度Fig.5 Terrain complexity

表7 地形復(fù)雜度因子分級(jí)

通過(guò)變異系數(shù)法計(jì)算坡度、地表粗糙度、標(biāo)準(zhǔn)曲率和地表切割深度4項(xiàng)地形因子的權(quán)重。經(jīng)計(jì)算，在TCI模型中的權(quán)重分別為0.282、0.171、0.259和0.288。構(gòu)建TCI模型如圖6所示，圖中F1、F2和F3標(biāo)出部分復(fù)雜度高的區(qū)域，F(xiàn)1、F2為山脊和山谷主要分布區(qū)域。該區(qū)域復(fù)雜度較高不僅因脊- 谷區(qū)表面崎嶇，經(jīng)實(shí)地野外踏勘，試驗(yàn)區(qū)裸土、裸巖分布廣泛，且紅壤和紫色土黏性低。 F1、F2區(qū)域經(jīng)雨水沖刷后，腳踩凹陷明顯，紅壤和紫色土沙化嚴(yán)重。土壤侵蝕是土壤在水力、風(fēng)力等外力作用下，被剝蝕、破壞和分離的過(guò)程。TCI模型中可知，試驗(yàn)區(qū)較高和高復(fù)雜度區(qū)域主要分布于山體兩側(cè)，山體兩側(cè)凹凸不平說(shuō)明山體受亞熱帶季風(fēng)性氣候影響，經(jīng)風(fēng)化作用導(dǎo)致土壤受到侵蝕，從而反映出以F1、F2為主的較高和高復(fù)雜度區(qū)域內(nèi)水土流失嚴(yán)重。F3為試驗(yàn)區(qū)位于試驗(yàn)區(qū)東北角，實(shí)地踏勘時(shí)發(fā)現(xiàn)，該區(qū)域不僅有沙化嚴(yán)重的紅壤和紫色土，且分布著白色碎化的沉積巖。據(jù)表8統(tǒng)計(jì)分析表明較高復(fù)雜度和高復(fù)雜度總體面積比例相對(duì)集中，二者總占測(cè)區(qū)面積的22.95%。

F1、F2和F3是復(fù)雜度高的區(qū)域. F1, F2 and F3 are areas of high complexity.圖6 TCI分級(jí)圖Fig.6 TCI grading chart

表8 TCI分級(jí)表

4.2 TCI模型運(yùn)用驗(yàn)證分析討論

采用可視化對(duì)比分析方法[8]驗(yàn)證TCI模型的有效性。依據(jù)TCI模型，通過(guò)試錯(cuò)法由0～100%調(diào)整模型透明度，發(fā)現(xiàn)60%透明度模型既能清晰展示模型，在山體陰影和等高線疊加又能較好突出山體和等高線變化趨勢(shì)，如圖7示。經(jīng)對(duì)比分析可知TCI模型符合等高線變化程度，對(duì)于F1和F2區(qū)域在TCI模型中為高復(fù)雜度區(qū)域，等高線走向同樣貼合模型。因此，能表明TCI可較好反映地表形態(tài)變化。

F1和F2是復(fù)雜度高的區(qū)域. F1 and F2 are areas of high complexity. 圖7 TCI模型驗(yàn)證Fig.7 TCI model validation

5 結(jié)論

無(wú)人機(jī)具有攜帶輕便、人工成本低的特點(diǎn)，可快速獲取高精度DEM。以此為基礎(chǔ)，依據(jù)對(duì)復(fù)雜地形地貌進(jìn)行完整性和適用性的分級(jí)處理，結(jié)合曲面傾斜度和地表崎嶇狀況考慮，針對(duì)4項(xiàng)地形單因子指標(biāo)，分別賦予不同權(quán)重構(gòu)建TCI模型。試驗(yàn)研究結(jié)果表明：1)依據(jù)變異系數(shù)計(jì)算的各項(xiàng)單因子權(quán)重分為0.282、0.171、0.259和0.288；2)TCI模型中復(fù)雜度較高的F1、F2區(qū)域脊- 谷的位置表面崎嶇，經(jīng)實(shí)地野外踏勘，試驗(yàn)區(qū)裸土、裸巖分布廣泛，且紅壤和紫色土黏性低；3)結(jié)合恐龍谷南緣山區(qū)的氣候、土壤條件和實(shí)地踏勘情況，發(fā)現(xiàn)TCI模型中較高復(fù)雜、高復(fù)雜度區(qū)域的土壤受一定程度侵蝕導(dǎo)致水土流失。大面積侵蝕地貌的情況下，以小面積TCI模型快速判斷土壤受侵蝕區(qū)域主要集中于脊- 谷位置。無(wú)人機(jī)影像構(gòu)建DEM，相對(duì)于其余技術(shù)手段，不僅省時(shí)省力，且精度相對(duì)較高。在大面積侵蝕地貌的情況下，小面積、高精度的DEM構(gòu)建TCI模型，以地表地形變化結(jié)合氣候、土壤條件從地表直觀角度分析土壤受侵蝕區(qū)域，依此為基礎(chǔ)快速判斷大面積土壤侵蝕區(qū)域?；诟叻直媛蔇EM在突出微地貌地形特征的同時(shí)相對(duì)更符合實(shí)地情況反映土壤侵蝕狀況。對(duì)于TCI模型，可依據(jù)不同的地形條件，選取適宜的地形因子。