基于SRTM的地形因子提取方法研究

張宏鳴 常 毅 楊勤科 張 泉 董 良 許伊昆

(1.西北農(nóng)林科技大學(xué)信息工程學(xué)院， 陜西楊凌 712100； 2.西北大學(xué)城市與環(huán)境學(xué)院， 西安 710069)

0 引言

地形是土壤侵蝕和水土保持措施布設(shè)的主要影響因素。在通用土壤流失方程及其修訂版(Universal soil loss equation/revised universal soil loss equation，USLE/REUSLE)[1-2]和中國土壤流失方程(Chinese soil loss equation，CSLE)[3]中，地形(LS)因子是定量計算土壤流失的重要指標(biāo)。區(qū)域尺度上，通常使用數(shù)字高程模型(Digital elevation model，DEM)進(jìn)行地形因子提取[4]。隨著無人機(jī)技術(shù)、航空攝影測量技術(shù)、激光雷達(dá)技術(shù)的發(fā)展，高精度、大范圍DEM獲取逐漸容易[5]，如何在大區(qū)域上使用高精度DEM進(jìn)行地形因子提取成為研究重點。

SRTM(Shuttle radar topography mission)是美國太空總署和國防部國家測繪局以及德國和意大利航天機(jī)構(gòu)共同合作完成的地球表面遙感測量結(jié)果[6]。該數(shù)據(jù)憑借在全球尺度上精度高、范圍大等優(yōu)勢[7]，可為大區(qū)域甚至全球范圍地形因子提取提供數(shù)據(jù)基礎(chǔ)[8]。SRTM處于地理坐標(biāo)系，而已有的LS因子提取方法通常基于投影坐標(biāo)系，通過提取坡度、坡長，獲取LS因子[9]?，F(xiàn)階段，坡度提取算法較多，如最陡坡降、三階反距離平方權(quán)、三階不帶權(quán)差分等，獲取較容易；坡長由于計算原理復(fù)雜，需要考慮截斷等多種復(fù)雜因素，獲取較困難[10]。文獻(xiàn)[11-16]通過匯水面積法替換或徑流線換算的方法代替坡長求得LS因子，但這些方法無法考慮截斷問題，也不能獲得坡長。文獻(xiàn)[17-19]從坡長基本定義出發(fā)，考慮坡長的起點和終點，獲取精確的坡長，根據(jù)坡長和坡度計算LS因子。文獻(xiàn)[20-21]研究了坡度截斷和溝道截斷對侵蝕坡長提取的影響，可進(jìn)一步提高坡長精度。文獻(xiàn)[22-23]表明，使用分段坡長法提取到的地形因子準(zhǔn)確性及空間分布更加合理。前人對地形因子提取的相關(guān)機(jī)理作了大量研究，投影坐標(biāo)系下地形因子提取方法較為成熟。但在實際應(yīng)用時，要求柵格單元以米為單位，且邊長一致，無法直接對SRTM進(jìn)行計算。通常需對SRTM進(jìn)行坐標(biāo)轉(zhuǎn)換后方可應(yīng)用，計算效率低，不利于大區(qū)域地形因子快速提取。文獻(xiàn)[24]通過分析地理坐標(biāo)系特點，構(gòu)建地球球體模型，提出SRTM的坡度和坡向算法，但未對坡長與LS提取方法進(jìn)行研究與應(yīng)用。

基于此，本文提出一種基于SRTM的LS因子提取算法(LSA-SRTM)。通過模擬地球球體模型，構(gòu)建換算關(guān)系及解算方法，直接在地理坐標(biāo)系下提取SRTM柵格單元邊長、獲取坡度與流向、提取坡長，進(jìn)而提取LS因子。

1 研究數(shù)據(jù)和方法

1.1 實驗數(shù)據(jù)

為了驗證算法的精度，需使用不同地貌類型的高程數(shù)據(jù)進(jìn)行評估[16]。但實際地形的微觀特征過多，數(shù)據(jù)本身存在一定的隨機(jī)噪聲[25]，文獻(xiàn)[26]嘗試構(gòu)建數(shù)學(xué)理論曲面代替實際地形，并在曲面上考慮實際地形可能存在的基本地貌類型。由于數(shù)學(xué)曲面噪聲、誤差等影響較小，文獻(xiàn)[16，27]基于數(shù)學(xué)曲面進(jìn)行了流向算法分析以及改進(jìn)方法研究。Himmelblau-Orlandini數(shù)學(xué)曲面是文獻(xiàn)[28]利用Himmelblau函數(shù)，通過仿射變換后生成的離散曲面，是一種復(fù)雜數(shù)學(xué)曲面，兼具凹凸面、發(fā)散匯集等特征，曲面上同時有山頂、鞍部、谷底等地形特征，與實際地形條件吻合，且通過其生成模型，可計算理論真值，方便進(jìn)行誤差、精度檢驗，同時可以避免投影轉(zhuǎn)換的誤差，減少誤差的累積效應(yīng)。

Himmelblau-Orlandini曲面的原始函數(shù)為

Z(X,Y)=(X2+Y-4)2+(Y2+X-7)2
(-5≤X≤5；-5≤Y≤5)

(1)

為了方便該曲面進(jìn)行DEM離散表達(dá)，使用

x=50X+250

(2)

y=50Y+250

(3)

z=-0.75Z+450

(4)

式中x——x軸坐標(biāo)值，取0～500 m

y——y軸坐標(biāo)值，取0～500 m

z——(x,y)的高程

對曲面進(jìn)行仿射變換[28],生成0～500 m尺度的SRTM如圖1所示。

圖1 數(shù)學(xué)曲面Fig.1 Mathematical surface

1″(分辨率約30 m)SRTM(來源：https:∥e4ftl01.cr.usgs.gov/)為HGT格式，采用WGS 1984坐標(biāo)系。整套數(shù)據(jù)覆蓋60°N～56°S，總面積超過1.19×108km2，單個數(shù)據(jù)文件覆蓋1°(經(jīng)度)×1°(緯度)的地區(qū)。目前，基于該數(shù)據(jù)的地形因子結(jié)果分析已有大量研究。文獻(xiàn)[29]的研究表明，DEM格網(wǎng)擴(kuò)張會降低LS的提取精度，1″ SRTM的LS因子計算精度高于3″ SRTM。文獻(xiàn)[30]的研究表明，在特定的海拔范圍，1″ SRTM提取的地形因子甚至優(yōu)于10 m分辨率的DEM。文獻(xiàn)[31]的研究表明1″ SRTM稍差于25 m的Hc-DEM(水文地貌關(guān)系正確的DEM)，但比相同分辨率的ASTER GDEM(先進(jìn)星載熱發(fā)射和反射幅射儀全球數(shù)字高程模型)地形表達(dá)能力更好。文獻(xiàn)[32]對1″ SRTM的侵蝕表達(dá)能力作了分析，認(rèn)為研究樣區(qū)較大、缺少連續(xù)的高精度DEM時，1″ SRTM可以代替使用。因此，目前全球尺度的土壤侵蝕多基于此數(shù)據(jù)進(jìn)行[8]。

選取全國土壤侵蝕強(qiáng)烈、侵蝕地貌發(fā)育的5個典型樣區(qū)為實驗區(qū)(圖2)。包括：東北漫崗丘陵區(qū)(黑龍江省齊齊哈爾市拜泉縣)，代表波狀起伏的緩坡丘陵；黃土丘陵區(qū)(陜西省延安市安塞縣)，代表半表深厚黃土覆蓋、強(qiáng)烈水蝕作用的陡坡丘陵；西南紫色土丘陵區(qū)(四川省綿陽市鹽亭縣)與南方紅壤丘陵區(qū)(浙江省金華市東陽縣)，代表亞熱帶濕潤多雨環(huán)境下的丘陵；北方土石山區(qū)(山東省臨沂市蒙陰縣)，代表暖溫帶半濕潤、半干旱環(huán)境下的山地地形景觀[33]。各個樣區(qū)尺寸均為1°×1°，樣區(qū)高程數(shù)據(jù)如圖3所示。

圖2 樣區(qū)分布圖Fig.2 Location of samples

圖3 典型樣區(qū)數(shù)字高程模型Fig.3 DEM samples

1.2 實驗方法

LSA-SRTM的核心過程主要分為5個步驟，算法流程如圖4所示：利用經(jīng)緯度信息計算柵格的單元邊長及單元坡長；利用Deterministic 8(D8)算法計算流向和坡度；依據(jù)坡度和截斷因子設(shè)置坡度截斷點，根據(jù)匯水面積和指定的閾值設(shè)置溝道截斷點；參考截斷位置，計算累積坡長；根據(jù)CSLE，利用坡度和累積坡長計算LS因子。輸入的SRTM數(shù)據(jù)是ASCII數(shù)據(jù)，在計算前對輸入數(shù)據(jù)的合法性進(jìn)行判斷。

圖4 基于SRTM的地形因子提取流程圖Fig.4 Flow chart of topographic factor extraction based on SRTM

1.2.1SRTM單元邊長計算

SRTM基于地理坐標(biāo)系，柵格以弧度作為基本單位存儲，與投影坐標(biāo)系統(tǒng)下的DEM柵格不同，每個柵格單元近似看作梯形，相鄰單元之間的距離只能通過地球半徑和經(jīng)緯度進(jìn)行計算。

圖5 地球球體模型Fig.5 Earth sphere model

圖6 地球經(jīng)線面和緯線面示意圖Fig.6 Schematics of longitude and latitude planes of earth

(5)

(6)

由r=Rcosφ有

(7)

由于SRTM單元在經(jīng)緯度上的跨度相同，以1″SRTM為例，其跨度為0.000 277 777 777 778°。取地球半徑R=6 371 000 m，ω=1(°)/3 600，SRTM柵格單元寬度在南北方向恒定為30.887 479 1 m，在東西方向隨緯度變化。

1.2.2坡度和流向提取

坡度和流向采用D8算法[34]計算。D8算法基于最陡坡降思想，在3×3的DEM網(wǎng)格上，計算中心單元格與各相鄰單元格間的坡度，取坡度結(jié)果最大值作為中心單元格的坡度，并且將最大值所在方向確定為中心單元格的流向[35]。如圖7所示，C為當(dāng)前單元所在的位置，其流出方向為周圍8個單元中的一個，分別標(biāo)記為1、2、4、8、16、32、64、128[11]。為確保每個柵格單元都與河網(wǎng)連接，當(dāng)單元格坡度為0時，取0.1[10]。

NW32N64NE128W16CE1SW8S4SE2

1.2.3坡長提取

坡長是從起點沿著垂直等高線方向到達(dá)坡長截斷終點的積分，基于柵格數(shù)據(jù)進(jìn)行計算時，坡長的解算主要由每一個柵格的單元坡長累積來計算完成[10]。其公式可表示為

(8)

式中λi,j——坐標(biāo)點(i,j)的坡長

λc——每個柵格的單元坡長

m——坡長指數(shù)

k——匯入坐標(biāo)點(x,y)的柵格數(shù)

本研究考慮截斷對坡長的影響，坡度截斷通過判斷坡度變化率和中斷因子的關(guān)系來確定，當(dāng)坡度變化率大于中斷因子時，該點標(biāo)記為截斷點[21]。溝道截斷通過設(shè)置閾值的方法確定，當(dāng)匯水面積大于閾值時，標(biāo)記該點為截斷點[20]。本文將匯水面積閾值設(shè)置為105m2。

累積坡長的計算方法是以起點柵格為基礎(chǔ)，沿周圍8個不同方向的最大坡降方向累加坡長。對于SRTM數(shù)據(jù)，無法確定某條流路方向的最大坡長，因此，需對柵格點正向、反向遍歷，從柵格數(shù)據(jù)起點逐點計算來完成累計坡長計算[19]。

1.2.4地形因子提取

在CSLE中，坡長和坡度共同決定了侵蝕地形因子的計算[29]。為避免只考慮均勻坡長帶來的誤差[36]，使用分段坡長因子公式來計算坡長因子。LS因子的計算公式為

(9)

(10)

其中

式中θ——坡度S——坡度因子

λin——入口坡長，m

λout——出口坡長，m

L——坡長因子

1.2.5實驗對比方法

(1)精度對比

在格網(wǎng)DEM上，完整獲取各個位置的地形因子真值非常困難，精度對比需對檢驗樣區(qū)選定樣點測量完成。參照文獻(xiàn)[37]，在實驗樣區(qū)中選取較為規(guī)整的坡面樣點，樣本空間設(shè)置為150，需避免將樣點設(shè)在明確的分水線和溝谷中[38]。根據(jù)DEM繪制(等間距小于cellsize)等高線，在高程變化平緩的區(qū)域減小等高線間距，采用ArcGIS 10.2的量尺工具，量算采樣?xùn)鸥裰行狞c沿與等高線垂直方向到明確分水線或山頂(水流起始點)的投影距離作為該柵格的實測坡長，取多次測量的均值作為最終結(jié)果，結(jié)果保留小數(shù)點后兩位，不考慮人為誤差時，有效測量精度為0.1 m?；趯崪y坡長計算得到的LS因子為實測LS因子。以投影坐標(biāo)系下的LS算法[19](簡稱為LSA-DEM)作為參照方法，將兩種方法的計算結(jié)果與實測結(jié)果進(jìn)行對比。

確定計算精度時，使用回歸分析方法以及決定系數(shù)R2，衡量計算結(jié)果與測量結(jié)果的相關(guān)性；使用標(biāo)準(zhǔn)差(SD)以及絕對差(AD)確定計算誤差。

(2)效率對比

在64位Dell Precision工作站(系統(tǒng)：Windows 10；RAM：64 GB；CPU：Intel Xeon Gold 5115；主頻：2.4 GHz)對兩種方法的運算時間進(jìn)行統(tǒng)計。樣本大小分別為：500 m×500 m，1°×1°、2°×2°、5°×5°、14°×14°。為避免不可控因素對運算時間測量的干擾，對每一測試樣本均進(jìn)行5次運算，取其運算時間的平均值作為實際執(zhí)行時間。

2 結(jié)果與討論

2.1 數(shù)學(xué)曲面提取結(jié)果分析

數(shù)學(xué)曲面上，LSA-SRTM計算結(jié)果見圖8a～8c，坡度最大值為84.97°，最小值為0.1°，平均值為50.91°(表1)。坡度高值均分布在4個局部高點外側(cè)陡坡區(qū)域，局部高點內(nèi)側(cè)坡度變化不明顯(圖8a)。僅在坡度截斷作用下，坡長最大值為407.65 m，最小值為0.48 m，坡長均值64.96 m。坡長從局部高點向坡度變化最陡的方向累加，在匯集坡面可累計至流域邊界，能反映地表起伏狀態(tài)(圖8b)。LSA-DEM提取結(jié)果見圖8d～8f，兩種方法的紋理特征無明顯差異，3種地形指標(biāo)的均值與標(biāo)準(zhǔn)差十分接近。

表1 數(shù)學(xué)曲面上LSA-SRTM和LSA-DEM結(jié)果Tab.1 Results of LSA-SRTM and LSA-DEM in mathematical surface

圖8 數(shù)學(xué)曲面上地形因子提取結(jié)果Fig.8 Topographic factors extraction results in mathematical surface

兩種方法與測量值對比如表2所示。由表2可見，兩種方法的結(jié)果與測量值均呈現(xiàn)較強(qiáng)的相關(guān)性。LSA-SRTM與測量坡長的R2為0.855 2，與測量LS因子的R2為0.890 7，由此可見，計算結(jié)果與測量結(jié)果有較高的相關(guān)性，可以較好地表達(dá)真實地形。計算結(jié)果與實測坡長的SD與AD分別為26.34、18.98 m，與實測LS因子的SD與AD分別為6.41、3.94。

表2 數(shù)學(xué)曲面上LSA-SRTM和LSA-DEM結(jié)果對比Tab.2 Differences between LSA-SRTM and LSA-DEM results in mathematical surface

計算結(jié)果與實測結(jié)果有較高的一致性，誤差范圍較小，從回歸分析的結(jié)果來看，兩種方法具有高度相似的誤差來源。誤差的主要原因在于單流向D8方法應(yīng)用的局限性，D8方法的流向只能選擇周圍8個柵格中的一個，這使得發(fā)散坡面上易出現(xiàn)柵格缺失上游水流入情況[39]，影響坡長的累積結(jié)果。

2.2 典型樣區(qū)提取結(jié)果分析

典型樣區(qū)上3種地形指標(biāo)統(tǒng)計結(jié)果如表3～5所示。

從表3可看出，漫崗丘陵區(qū)坡度最為平緩，平均坡度和坡度標(biāo)準(zhǔn)差均最小，坡度基本在6°以下；其次為北方土石山區(qū)，平均坡度較小，小于5°的面積占23%；整體地形由緩到陡依次為紫色土丘陵區(qū)、黃土丘陵區(qū)和紅壤丘陵區(qū)；平均坡度和標(biāo)準(zhǔn)差最大的是紅壤丘陵區(qū)。

表3 典型樣區(qū)LSA-SRTM與LSA-DEM坡度Tab.3 Slope gradient results of LSA-SRTM and LSA-DEM in typical samples (°)

從表4可知，平均坡長最大的是漫崗丘陵區(qū)，其次為紫色土丘陵區(qū)、紅壤丘陵區(qū)、北方土石山區(qū)，黃土丘陵區(qū)的平均坡長最小，與其支離破碎的地形相適應(yīng)，5個樣區(qū)的坡長大多在200 m以內(nèi)。

表4 典型樣區(qū)LSA-SRTM與LSA-DEM的坡長Tab.4 Slope length results of LSA-SRTM and LSA-DEM in typical samples m

從表5可知，LS因子由漫崗丘陵區(qū)、北方土石山區(qū)、紫色土丘陵區(qū)、黃土丘陵區(qū)和南方紅壤丘陵區(qū)依次增加。漫崗丘陵區(qū)LS因子在0.5以下的面積約75%；黃土丘陵區(qū)LS因子主要分布在5以上，其中大于20的部分占35%左右；西南紫色土丘陵區(qū)LS因子在0.5以下的部分約占5%，大于20的部分占18%左右；南方紅壤丘陵區(qū)LS因子在0.5以下的面積占19%，大于20的部分占23%左右；北方土石山區(qū)LS 因子主要分布在0.5～20之間，0.5以下占11%左右。該結(jié)果與文獻(xiàn)[31]的研究結(jié)果基本一致。

表5 典型樣區(qū)LSA-SRTM與LSA-DEM的LS因素Tab.5 LS factor results between LSA-SRTM and LSA-DEM in typical samples

由于篇幅原因，本文給出黃土丘陵區(qū)的結(jié)果示意圖(圖9)。從宏觀角度觀測，兩種算法得到的3

圖9 黃土丘陵區(qū)地形因子提取結(jié)果Fig.9 Extraction results of topographic factors in loess hilly

種地形指標(biāo)的空間格局和紋理分布十分接近。

典型樣區(qū)上兩種方法與測量值的對比結(jié)果見表6。由表6可見，5個樣區(qū)上LSA-SRTM與坡長測量值的R2分別為：0.778 8、0.726 9、0.702 4、0.690 9、0.725 5；對應(yīng)的SD分別為：46.12、57.38、43.66、44.63、39.41 m；AD分別為：32.16、38.20、31.19、26.69、19.75 m。與LS因子測量值的R2分別為：0.820 9、0.821 3、0.714 2、0.714 5、0.821 2；對應(yīng)的SD分別為：7.97、8.83、8.54、7.69、6.52；AD分別為：5.14、7.63、6.13、4.79、3.56。兩種方法與測量值的相關(guān)性接近，從計算誤差上看，LSA-DEM的SD與AD指標(biāo)均高于LSA-SRTM。主要原因是投影轉(zhuǎn)換導(dǎo)致高程在轉(zhuǎn)換前后發(fā)生改變、柵格點位偏移，對后續(xù)計算引起連鎖反應(yīng)。相較于數(shù)學(xué)曲面的結(jié)果，典型樣區(qū)上，坡長與LS的精度都有所下降，主要與實驗數(shù)據(jù)的分辨率下降有關(guān)。

表6 典型樣區(qū)上LSA-SRTM和LSA-DEM結(jié)果對比Tab.6 Differences between LSA-SRTM and LSA-DEM results in typical samples

2.3 效率分析

LSA-SRTM與LSA-DEM的運行時間如表7所示。在500 m×500 m(數(shù)學(xué)曲面)的樣本上，由于兩種方法計算的柵格數(shù)一樣，且不需要投影轉(zhuǎn)換時，LSA-DEM與LSA-SRTM的運行時間持平。在1″SRTM上，LSA-DEM的運算時間更長，主要與投影變換的時間開銷，以及變換引起的輸入計算的柵格數(shù)增多有關(guān)。LSA-SRTM可以有效節(jié)省坐標(biāo)轉(zhuǎn)換時間，計算效率更高。

表7 LSA-SRTM與LSA-DEM運行時間Tab.7 Running time of LSA-SRTM and LSA-DEM

3 結(jié)束語

針對傳統(tǒng)LS方法計算SRTM需投影轉(zhuǎn)換、效率低的問題，本文提出了一種基于SRTM的地形因子提取算法(LSA-SRTM)。并對該算法的精度和效率與投影坐標(biāo)系下的LS提取算法(LSA-DEM)進(jìn)行對比分析。通過實驗驗證得出，在5個典型樣區(qū)上，LSA-SRTM的LS結(jié)果與實測LS的R2分別為：0.820 9、0.821 3、0.714 2、0.714 5、0.821 2，精度高于LS-ADEM；SD分別為：7.97、8.83、8.54、7.69、6.52，誤差小于LSA-DEM。同時，LSA-SRTM保持計算前后柵格的規(guī)模、節(jié)省投影轉(zhuǎn)換的時間開銷，當(dāng)柵格尺寸為50 401×50 401時，計算時間比LSA-DEM少約1 166.443 s，提高了地形因子提取效率。本文方法在保證計算精度的同時，提高了地形因子提取效率，可為大區(qū)域上地形因子研究提供支撐。