基于DEM的小班坡度自動提取算法及其驗(yàn)證

陳 晨 陳永剛 徐文兵 梁 丹

( 1. 浙江農(nóng)林大學(xué)浙江省森林生態(tài)系統(tǒng)碳循環(huán)與固碳減排重點(diǎn)實(shí)驗(yàn)室，浙江 杭州 311300；2. 同濟(jì)大學(xué)測繪與地理信息學(xué)院，上海 200092)

地面坡度（簡稱坡度）是對地面傾斜程度的定量描述，是一個基本的地貌形態(tài)指標(biāo)[1]，其大小直接影響著地表物質(zhì)流動和能量轉(zhuǎn)換的規(guī)模和強(qiáng)度，是工程設(shè)計(jì)和生產(chǎn)力空間布局的重要因子[2]，也是研究水文模型[3]、土壤侵蝕[4]、監(jiān)測滑坡[5]和規(guī)劃土地利用[6]的基礎(chǔ)。小班坡度是森林資源規(guī)劃設(shè)計(jì)調(diào)查（二類調(diào)查）中的一項(xiàng)重要因子，它是反映小班地表形態(tài)、小班立地質(zhì)量、生產(chǎn)作業(yè)條件及水土流失的重要指標(biāo)之一，也是林木采伐首先要考慮的因素之一[7]。在森林資源連續(xù)清查的小班調(diào)查中，小班坡度是重要的調(diào)查因子。因此，快速準(zhǔn)確的提取小班坡度顯得至關(guān)重要。

在以往的二類調(diào)查中，提取小班坡度的傳統(tǒng)方法主要有兩種：一是外業(yè)實(shí)地測量小班坡度。該方法主要是在實(shí)地進(jìn)行小班區(qū)劃后，利用坡度儀或者手持羅盤儀實(shí)地測量小班的平均坡度；二是利用地形圖進(jìn)行室內(nèi)估算，即在小班內(nèi)兩處等高線的水平距離及垂直距離，根據(jù)地形圖提供的坡度尺，計(jì)算小班的坡度值。利用傳統(tǒng)方法提取的二類小班坡度雖然很直觀，但主要依靠手工操作完成，導(dǎo)致工作效率較低，工作量大，時間周期長等問題，并且在野外實(shí)施也有一定的難度；另外，由于依靠調(diào)查人員的經(jīng)驗(yàn)判斷，導(dǎo)致結(jié)果具有一定的主觀性及隨意性，受人為因素影響大，數(shù)據(jù)不準(zhǔn)確，難以滿足二類調(diào)查的要求。研究人員也提出了一些自動化提取小班坡度的方法。陳建義等[7]針對傳統(tǒng)方法的不足，探討了一種基于ArcGIS軟件的小班坡度自動求算方法。熊安華[8]、葉江霞等[9]均利用ArcGIS進(jìn)行小班坡度的計(jì)算。這些基于數(shù)字高程模型（DEM）并結(jié)合GIS分析工具的方法能夠高效方便地計(jì)算出小班的坡度，但是這些方法只計(jì)算了小班內(nèi)所有柵格點(diǎn)坡度的平均數(shù)或者眾數(shù)。這種以簡單數(shù)學(xué)統(tǒng)計(jì)結(jié)果作為整個圖斑坡度表達(dá)的方法忽略了地形因素，并不能作為小班內(nèi)實(shí)際坡度的真實(shí)反映，這種缺陷在地形復(fù)雜度較高的地區(qū)尤為明顯[10]。

針對林業(yè)小班坡度計(jì)算的現(xiàn)狀，本研究提出了一種基于DEM的小班坡度自動提取算法。該算法先將小班內(nèi)真實(shí)地形中存在的異常高低起伏點(diǎn)進(jìn)行剔除，然后將小班整體作為計(jì)算對象，對剔除異常點(diǎn)的DEM進(jìn)行回歸，得到小班所在的斜平面模型，最后根據(jù)坡度公式計(jì)算該小班的坡度值。該算法能夠利用DEM快速準(zhǔn)確地提取小班的坡度，并且能夠適應(yīng)地形復(fù)雜度較高的區(qū)域，實(shí)現(xiàn)代替人工實(shí)測提取小班坡度。

1 基于DEM的小班坡度自動提取的算法

1.1 算法思路

在計(jì)算小班坡度值時，理想的小班內(nèi)部地形表面是平整光滑的，然而在真實(shí)地形中，小班內(nèi)的地形表面總是存在著一些高低起伏，見圖1～2。在計(jì)算小班的坡度時，這些異常的高低起伏在一定程度上影響小班坡度計(jì)算的準(zhǔn)確度。為了更加準(zhǔn)確地計(jì)算小班的坡度，需要將這些高低起伏的地方進(jìn)行剔除。因此，本研究提出了一種剔除異常高低起伏點(diǎn)的小班坡度自動提取算法（FSS），其算法流程見圖3。采用線性回歸獲得小班所在的斜平面模型，從而能夠綜合考慮小班整體情況。

圖 1 理想小班內(nèi)部地形表面Fig. 1 Ideal subcompartment internal terrain surface

圖 2 山坡上異常高低起伏點(diǎn)Fig. 2 Abnormal high and low undulations on the slope

圖 3 算法流程Fig. 3 The flow of experiment

1.2 識別并剔除異常高低起伏點(diǎn)

首先，借助ArcGIS的綜合工具，使用3×3的鄰域窗口，將DEM進(jìn)行均值綜合。然后，使用雙線性插值的方法將綜合后的DEM重采樣至原始DEM的分辨率。雙線性插值方法基于4個最近輸入單元中心的加權(quán)距離平均來確定單元新值。最后，將重采樣后的DEM與原始DEM作差，由此得到原始DEM與重采樣DEM之間的高程差。對DEM高程差進(jìn)行分析統(tǒng)計(jì)，將超過2倍方差的點(diǎn)視為異常的高低起伏點(diǎn)。在本研究中，小班內(nèi)的DEM柵格數(shù)量均大于30，假設(shè)服從正態(tài)分布。剔除異常高低起伏點(diǎn)的公式如下：

式中：X為小班內(nèi)DEM的高程差，為小班內(nèi)高程差的算術(shù)平均數(shù)，為小班下高程差的標(biāo)準(zhǔn)差。

1.3 建立小班斜平面模型

剔除異常高低起伏點(diǎn)后，利用線性回歸來建立小班所在斜平面的回歸模型，回歸模型的方程如下：

式中：Z是斜平面模型的因變量，此處為小班內(nèi)所有點(diǎn)的高程值集合，X1和X2是斜平面模型中的自變量，此處分別為小班內(nèi)所有點(diǎn)的地理坐標(biāo)X值和Y值集合，n是單個小班內(nèi)DEM柵格點(diǎn)的個數(shù)，zi是第i個點(diǎn)的高程值，xi和yi則分別表示第i個點(diǎn)的地理坐標(biāo)的X值和Y值。b1和b2是X1和X2的系數(shù)，ε是誤差值，a 是常數(shù)項(xiàng)。

在算得每個小班的斜平面模型后，根據(jù)坡度計(jì)算公式得到小班的坡度。坡度計(jì)算公式如下：

式中：fx為南北方向高程變化率，fy為東西方向高程變化率。此處的fx和fy分別等于式（2）中的 b1和 b2。

2 檢驗(yàn)數(shù)據(jù)來源及方法

2.1 數(shù)據(jù)來源

采用寧?？h2007年二類調(diào)查成果數(shù)據(jù)，該數(shù)據(jù)為shapefile格式，坐標(biāo)系為1980年國家大地坐標(biāo)系，該數(shù)據(jù)作為本研究的實(shí)測數(shù)據(jù)。選取68個小班作為樣本，這68個小班均未出現(xiàn)跨山脊和山谷的現(xiàn)象，具體小班分布見圖4。DEM數(shù)據(jù)是通過對1∶50 000的等高線數(shù)據(jù)進(jìn)行加工而來，同樣為1980年國家大地坐標(biāo)系。DEM由國家測繪部門提供。用等高線數(shù)據(jù)生成不規(guī)則三角網(wǎng)（TIN）；然后將TIN進(jìn)行柵格化，采用自然鄰域法，從而生成像元大小為10 m×10 m的DEM數(shù)據(jù)；最后，對DEM數(shù)據(jù)進(jìn)行降噪，填補(bǔ)洼地，以提高DEM數(shù)據(jù)的精度。

圖 4 實(shí)驗(yàn)區(qū)Fig. 4 Study area

2.2 檢驗(yàn)方法

檢驗(yàn)主要側(cè)重于算法準(zhǔn)確性和實(shí)用性的檢驗(yàn)。準(zhǔn)確性檢驗(yàn)就是將FSS算法與陳建義等[7]的均值法以及常用的中心點(diǎn)法（將小班中心點(diǎn)的坡度值作為小班的坡度值）計(jì)算小班坡度并進(jìn)行對比驗(yàn)證。由于無法獲得小班的真實(shí)坡度值，本研究利用模擬斜平面進(jìn)行驗(yàn)證，通過向斜平面上添加隨機(jī)噪聲來模擬現(xiàn)實(shí)世界中地表的起伏。若計(jì)算所得結(jié)果更加接近預(yù)設(shè)值，則認(rèn)為該方法得到的小班坡度更加準(zhǔn)確，能夠適應(yīng)復(fù)雜度較高的區(qū)域。為測試FSS算法的實(shí)用性，采用森林資源二類調(diào)查數(shù)據(jù)與本算法的結(jié)果進(jìn)行對比。驗(yàn)證方法均為t檢驗(yàn)。

為了進(jìn)行算法檢驗(yàn)，本研究首先構(gòu)造了一個30°的斜平面和一組隨機(jī)噪聲。其次，為模擬真實(shí)地形，向斜平面中分別加入隨機(jī)噪聲，噪聲分別為指數(shù)噪聲、伽瑪噪聲、整數(shù)噪聲、正常噪聲、泊松噪聲和均勻噪聲（圖5）。模擬斜平面和噪聲的像元大小與實(shí)驗(yàn)區(qū)DEM數(shù)據(jù)相同。最后，在模擬斜平面上繪制1塊圖斑，用于模擬小班，保存成shapefile格式。

圖 5 增加噪聲的模擬地形Fig. 5 Simulated terrain with increased noise

3 檢驗(yàn)結(jié)果與分析

3.1 準(zhǔn)確性檢驗(yàn)

使用均值法、中心點(diǎn)法和FSS算法分別計(jì)算模擬斜面上小班的坡度。在95%置信區(qū)間內(nèi)，通過配對t檢驗(yàn)對這3組數(shù)據(jù)進(jìn)行檢驗(yàn)。其中，原假設(shè)和備擇假設(shè)是：H0: abs_A≥abs_B，H1: abs_A＜abs_B。其中，abs_A和abs_B分別表示數(shù)據(jù)A和數(shù)據(jù)B與真值差的絕對值。因?yàn)楸狙芯啃枰容^FSS算法和其他2種方法計(jì)算結(jié)果的準(zhǔn)確度，所以此處的abs_A和abs_B有以下2種情況：|FSS-真值|和|均值法-真值|、|FSS-真值|和|中心點(diǎn)法-真值|。

配對t檢驗(yàn)的FSS算法與均值法和中心的檢驗(yàn)結(jié)果見表1。由表1可知，P值均小于0.05，即拒絕原假設(shè)，應(yīng)用FSS方法所得坡度值的均值更接近于真值。由表2可知，F(xiàn)SS方法提取的小班坡度與真值偏差的均值、標(biāo)準(zhǔn)差、最大值和最小值均比其他2種方法的結(jié)果小。FSS方法能夠準(zhǔn)確地提取小班坡度，同時面對實(shí)驗(yàn)預(yù)設(shè)的復(fù)雜地形，其相較于其他方法有更好的適應(yīng)性。

表 1 配對t檢驗(yàn)的結(jié)果比較Table 1 Comparison of the results of paired t-test

表 2 各方法提取小班坡度結(jié)果與真值的偏差Table 2 Each method extracts the deviation of the subcompartment slope result from the true value

3.2 實(shí)用性檢驗(yàn)

使用FSS算法對實(shí)驗(yàn)區(qū)數(shù)據(jù)進(jìn)行試驗(yàn)，得到了對應(yīng)小班的坡度值。用FSS算法計(jì)算得到的值與人工實(shí)測值進(jìn)行配對t檢驗(yàn)，其中，原假設(shè)和備擇假設(shè)是：H0: data_A=data_B，H1: data_A≠data_B。其中，data_A和data_B分別表示FSS算法結(jié)果和實(shí)測值數(shù)據(jù)。由表3可知，在95%的置信區(qū)間內(nèi)，這2組數(shù)據(jù)的P值為0.281，P值大于0.05，接收原假設(shè)，即FSS算法結(jié)果與人工實(shí)測數(shù)據(jù)無明顯差異。

表 3 實(shí)測數(shù)據(jù)與FSS算法結(jié)果配對t檢驗(yàn)的結(jié)果比較Table 3 Comparison of the results of paired t-test between measured data and FSS algorithm results

4 結(jié)論與討論

針對現(xiàn)有計(jì)算小班坡度算法的不足，本研究提出了一種基于DEM的小班坡度提取算法。該算法能夠準(zhǔn)確的提取小班坡度值，并且其結(jié)果與人工實(shí)測值無區(qū)別，能夠替代人工野外實(shí)測。FSS算法比常用的算法更加合理，原因可以從以下三方面進(jìn)行解釋：

1）小班內(nèi)的地形并不是光滑平坦的，坡面上存在著異常的高低起伏點(diǎn)。這些高低起伏點(diǎn)會影響小班坡度計(jì)算的準(zhǔn)確性。本研究的算法能夠剔除這些異常高低起伏點(diǎn)。

2）FSS方法利用線性回歸獲取了小班所在斜平面模型，綜合地考慮小班整體情況，而非通過簡單的統(tǒng)計(jì)計(jì)算得到小班的坡度值。因此，本研究所提出的算法結(jié)果比常用的算法更加準(zhǔn)確。

3）本研究所提出的算法是一種自動化的提取方法，對于野外人工調(diào)查具有一定的參考價值?？梢怨?jié)省人力、物力和時間，給森林資源調(diào)查提供了便利。

實(shí)際上，小班坡度算法準(zhǔn)確性檢驗(yàn)是一項(xiàng)很具有挑戰(zhàn)性的研究。由于無法獲得小班坡度的真實(shí)值，現(xiàn)有條件得到的小班坡度值均為測量值或其他方法所得，這些類型的數(shù)據(jù)并不準(zhǔn)確。在相關(guān)的基于DEM的坡度算法研究中，格網(wǎng)的坡度值通常使用模擬曲面來檢驗(yàn)，以此來評估算法的合理性，例如Chen等[11]，F(xiàn)lorinsky[12]，Zhou等[13]為了檢驗(yàn)DEM提取坡度值的精度，在模擬曲面上進(jìn)行算法比較。借鑒以上研究的思路，本研究構(gòu)建了一種適合于小班坡度算法的模擬地形曲面，與之前不同的是該曲面是一個添加了噪聲的斜平面，雖然它無法完全模擬真實(shí)地形，但是其坡度值可以作為一個理想真值。在異常高起伏點(diǎn)的剔除過程中，本研究選用了最為常用的3×3鄰域窗口，并且該窗口在實(shí)驗(yàn)中達(dá)到了預(yù)期效果。在實(shí)際應(yīng)用研究中，窗口大小的選擇受地形和DEM分辨率等因素影響[14-15]。因此，在今后研究中，將對鄰域窗口的選擇與DEM分辨率、異常區(qū)域大小等因素進(jìn)行深入探討。