地形起伏度最佳統(tǒng)計單元算法的比較研究

陳學(xué)兄, 常慶瑞, 畢如田, 劉正春, 張小軍

(1.山西農(nóng)業(yè)大學(xué) 資源環(huán)境學(xué)院, 山西 太谷 030801; 2.西北農(nóng)林科技大學(xué)資源環(huán)境學(xué)院, 陜西 楊凌 712100; 3.山西農(nóng)業(yè)大學(xué) 園藝學(xué)院, 山西 太谷 030801)

地形起伏度是定量描述地貌形態(tài)，宏觀描述地形變化整體趨勢的主要指標(biāo)[1-2]，也是導(dǎo)致水土流失的最直接因素，它不僅具有數(shù)學(xué)意義且具備土壤侵蝕和地貌學(xué)的意義，已被廣泛應(yīng)用于各個研究領(lǐng)域，如水土流失定量評價[3]、潛在水土流失評價[4]、地質(zhì)環(huán)境質(zhì)量定量評價[5]、區(qū)域滑坡災(zāi)害評價[6]、生態(tài)敏感性評價[7]、水土流失敏感性評價[8-9]等。在此基礎(chǔ)上，國內(nèi)學(xué)者對地形起伏度最佳統(tǒng)計單元的計算方法有了一定的研究[10-12]，最佳統(tǒng)計單元的確定屬于典型的變點(diǎn)分析問題，其計算方法有：人工作圖法、最大高差法、模糊數(shù)學(xué)法、均值變點(diǎn)分析法和累積和CUSUM分析算法，除人工作圖法外，其他方法都以定量的形式尋找突變點(diǎn)，哪種方法更適合山西省地形起伏度最佳統(tǒng)計單元的計算還需研究。鑒于此，本研究使用SRTM DEM數(shù)據(jù)，分析并探討最佳統(tǒng)計單元各計算方法的準(zhǔn)確性及適用性，確定山西省地形起伏度提取的最佳統(tǒng)計單元大小，以期為生態(tài)敏感性評價、水土流失定量評價等相關(guān)研究提供數(shù)據(jù)資料，也為山西省地形起伏度的進(jìn)一步應(yīng)用研究提供參考。

1 研究區(qū)概括

山西省，位于我國華北西部的黃土高原地帶，東隔太行山與河北為鄰，西隔黃河與陜西相望，南與河南省接壤，北與內(nèi)蒙古自治區(qū)相連，疆域輪廓呈由東北斜向西南的平行四邊形，全省總面積約16萬km2，地貌類型復(fù)雜多樣，有平原、臺地、山地、丘陵、高原等，其中山地丘陵約占全省總面積的80%之余。山西屬溫帶暖溫帶、半濕潤半干旱大陸性季風(fēng)氣候，降水偏少，年降水量為400～600 mm，自西北向東南遞增。境內(nèi)礦產(chǎn)資源豐富多樣，且分布廣泛，尤以煤炭資源儲量大，質(zhì)量優(yōu)，品種齊全，易開采而著稱于世。

2 數(shù)據(jù)與研究方法

2.1 基礎(chǔ)數(shù)據(jù)

本研究使用分辨率為90 m的SRTM DEM數(shù)據(jù)，數(shù)據(jù)采樣格網(wǎng)大小為90 m，格式為ArcGIS的Grid格式，地理坐標(biāo)系為：WGS_1984_Albers。

2.2 研究方法

2.2.1 地形起伏度的提取 使用鄰域分析法提取山西省地形起伏度，所使用矩形窗口大小為n×n像元(n=2，3，…，30)。地形起伏度提取的步驟如下：(1) 統(tǒng)計n×n(n=2，3，…，30)窗口內(nèi)像元的最大和最小值；(2) 計算各窗口內(nèi)最大值和最小值的差值；(3) 統(tǒng)計不同窗口的平均起伏度值(表1)。

2.2.2 最佳統(tǒng)計單元的計算 運(yùn)用人工作圖法[10-11]、最大高差法[10]、模糊數(shù)學(xué)法[10]、均值變點(diǎn)分析法[1,11,13]和累積和CUSUM分析算法[12]分別計算最佳統(tǒng)計單元，研究各計算方法在確定最佳統(tǒng)計單元上的區(qū)別與適用性。

(1) 人工作圖法。利用統(tǒng)計軟件對平均起伏度進(jìn)行回歸分析，提取擬合度最好的擬合曲線，從擬合曲線上人工判斷拐點(diǎn)。

(2) 最大高差法。在(i，j)樣本集中，樣點(diǎn)i在鄰域面積j處的高差變化值(dH)為：dHj=Hj-H(j-1)，比較i點(diǎn)所有的高差變化值dHj，取最大的dHj所對應(yīng)的j點(diǎn)為拐點(diǎn)，則拐點(diǎn)所對應(yīng)的面積為最佳統(tǒng)計單元的面積。

(3) 模糊數(shù)學(xué)法。樣本集(i，j)中，樣點(diǎn)i在j鄰域面積處，有：

(1)

(2)

(3)

式中：模糊系數(shù)(ak)均勻分布于(0,1]區(qū)間，其中2≥k≤N，N為鄰域面積個數(shù)(N=29)，a2=1；定義拐點(diǎn)為最大Rj值對應(yīng)的那一點(diǎn)。

(4) 均值變點(diǎn)分析法。均值變點(diǎn)分析法對恰有一個變點(diǎn)的檢驗(yàn)最有效，該方法已被廣泛應(yīng)用于計算地形起伏度最佳統(tǒng)計單元，根據(jù)均值變點(diǎn)分析法的原理構(gòu)建樣本序列，具體過程為[11]：首先，計算各統(tǒng)計單元下的單位地勢度(T)，單位地勢度(Ti)=平均起伏度(ti)/鄰域面積(si)，(i=2,3,…,30)；其次，對單位地勢度(T)取對數(shù)，得到序列X，X為{xk,k=1，2，3，…，29}；最后，根據(jù)下面公式對表1數(shù)據(jù)進(jìn)行計算處理，計算出統(tǒng)計量S和Si的值。

(4)

(5)

(6)

3 結(jié)果與分析

3.1 地形起伏度與鄰域面積關(guān)系分析

計算不同統(tǒng)計單元的鄰域面積(表1)。由表1可知，平均起伏度值隨鄰域面積的增大逐漸增大，但當(dāng)鄰域面積達(dá)到一定的閾值時，其變化趨勢逐漸平緩。

3.2 山西地形起伏度最佳統(tǒng)計單元的探討

3.2.1 最佳統(tǒng)計單元的確定

(1) 人工作圖法。利用統(tǒng)計軟件對表1中的平均起伏度值進(jìn)行回歸分析，提取擬合度最好的2條擬合曲線，見圖1。

表1 統(tǒng)計單元與地形起伏度對應(yīng)關(guān)系

圖1 鄰域面積與平均地形起伏度對應(yīng)關(guān)系擬合曲線

由圖1可知，使用SRTM DEM數(shù)據(jù)分析研究山西省地形起伏度時，采用冪模型和對數(shù)模型對平均起伏度值和鄰域面積進(jìn)行擬合，冪模型的擬合效果最好，而對數(shù)模型次之，擬合方程分別為y=22.777x0.3853，y=46.95ln(x)-60.119；平均起伏度值隨鄰域面積的增大逐漸增大，但其增幅一直在降低并漸漸趨于平緩，由于從擬合曲線上人工判斷拐點(diǎn)容易受主觀因素的影響，所以此方法不參與比對分析。

(2) 最大高差法。根據(jù)最大高差法的含義，由表1數(shù)據(jù)得知平均起伏度的最大變化值在3×3處，即3×3(0.072 9 km2)為使用此方法計算得出的最佳統(tǒng)計單元。

(3) 模糊數(shù)學(xué)法。根據(jù)模糊數(shù)學(xué)法的概念，各變量的取值為ak=1，27/28，…，1/28(k=2，3，4，…，N)，N=29。計算結(jié)果見表2。

表2 模糊數(shù)學(xué)法統(tǒng)計結(jié)果

由表2可知，在j=3時，Rj的值為最大，即使用該方法計算得到的最佳統(tǒng)計單元為3×3(0.072 9 km2)。

(4) 均值變點(diǎn)分析法。根據(jù)均值變點(diǎn)分析法的計算過程，計算得出統(tǒng)計量S=22.23，Si的值見表3。

由表3可以得知，在第10個點(diǎn)時S與Si的差值達(dá)到最大，這一點(diǎn)就是所要求的點(diǎn)，即利用該方法計算得到的最佳統(tǒng)計單元為11×11(0.980 1 km2)。

(5) 累積和(CUSUM)分析算法。 為確定其突變情況，依據(jù)累積和算法分析突變點(diǎn)的原理對表1中的數(shù)據(jù)序列進(jìn)行了分析，分析原始數(shù)據(jù)序列中累加和Si與統(tǒng)計單元窗口的關(guān)系，見圖2。由圖2可知，隨著鄰域面積的增大，CUSUM值逐漸減小，最小值(-827.34)出現(xiàn)在14×14處，然后隨著鄰域面積的增大，累積和(CUSUM)的值也逐漸增大，因此可確定14×14(1.587 6 km2)為原始數(shù)據(jù)序列發(fā)生突變的位置。

表3 均值變點(diǎn)分析法的統(tǒng)計結(jié)果

圖2 平均起伏度的累積和CUSUM圖

3.2.2 計算結(jié)果對比分析 不同計算方法計算所得的最佳統(tǒng)計單元面積見表4。與已有研究成果[13]進(jìn)行比對分析，以驗(yàn)證不同計算方法其結(jié)果的準(zhǔn)確性。在ArcGIS支持下，首先將不同計算方法所得的最佳統(tǒng)計單元下的山西省地形起伏度分為5個等級：平坦(0～30 m)、微起伏(30～70 m)、小起伏(70～200 m)、中起伏(200～500 m)和大起伏(>500 m)，統(tǒng)計各級所對應(yīng)的像元數(shù)，并計算各級所占的面積百分比(表5)，然后依據(jù)《中國1∶100萬地貌制圖規(guī)范(試行)中的地貌基本形態(tài)劃分體系》及其他相關(guān)研究成果[14-15]，對基本地貌類型進(jìn)行了簡單的劃分，并計算各類地貌所占的面積百分比，結(jié)果見表5。

表4 不同方法得到的統(tǒng)計單元面積結(jié)果對比

表5 不同方法計算結(jié)果分類

為比對不同方法其計算結(jié)果的準(zhǔn)確性，使用SPSS統(tǒng)計軟件分別分析最大高差法(或模糊數(shù)學(xué)法)、均值變點(diǎn)分析法和累積和算法與文獻(xiàn)[13]數(shù)據(jù)的相關(guān)性，分析得出均值變點(diǎn)分析法與文獻(xiàn)[13]中的起伏度最接近，相關(guān)系數(shù)為0.89，累積和分析算法次之，相關(guān)系數(shù)為0.757，最大高差法與模糊數(shù)學(xué)法相關(guān)性最差，相關(guān)系數(shù)為0.39，說明最大高差法與模糊數(shù)學(xué)法已不適用，這與前人研究結(jié)果相一致[10]，均值變點(diǎn)分析法和累積和分析算法是相對比較有效的方法，而且均值變點(diǎn)分析法更適合山西省地形起伏度最佳統(tǒng)計單元的確定。

通過上述分析可以得出：不同方法計算所得的最佳統(tǒng)計單元面積有所不同(表4)，人工作圖法易受主觀因素的影響，很難確定最佳統(tǒng)計單元，最大高差法和模糊數(shù)學(xué)法所得的最佳統(tǒng)計單元均為3×3(0.072 9 km2)，均值變點(diǎn)分析法所得的最佳統(tǒng)計單元為11×11(0.980 1 km2)，而累積和分析算法的最佳統(tǒng)計單元為14×14(1.587 6 km2)；由表5得知，均值變點(diǎn)分析法和累積和分析算法所得的最佳統(tǒng)計單元下的起伏度與文獻(xiàn)[13]中的起伏度較接近，以小起伏為主，其面積占研究區(qū)總面積的1/2左右，且這兩種方法劃分的基本地貌類型所占面積比例與實(shí)際較相符，山地丘陵占80%之余；均值變點(diǎn)分析法是計算山西省地形起伏度最佳統(tǒng)計單元的一種最為理想的方法。

4 結(jié)論與討論

(1) 基于90 m分辨率SRTM DEM數(shù)據(jù)，利用鄰域分析法在不同大小窗口下(2×2，3×3，4×4，…，30×30)對山西省地形起伏度進(jìn)行了提取(表1)，平均起伏度值起初隨鄰域面積的增大逐漸增大，但鄰域面積達(dá)到一定的閾值時，其變化趨勢逐漸平緩，采用冪模型對平均起伏度值和鄰域面積的擬合效果最好，對數(shù)模型次之；

(2) 運(yùn)用人工作圖法、最大高差法、模糊數(shù)學(xué)法、均值變點(diǎn)分析法和累積和(CUSUM)分析算法分別計算最佳統(tǒng)計單元。結(jié)果表明：人工作圖法易受主觀因素的影響，最大高差法和模糊數(shù)學(xué)法計算所得最佳統(tǒng)計單元下的地形起伏度分級結(jié)果與實(shí)際很不相符，運(yùn)用均值變點(diǎn)分析法與累積和分析算法計算所得的最佳統(tǒng)計單元分別為11×11(0.980 1 km2)和14×14(1.587 6 km2)，最佳統(tǒng)計單元下的地形起伏度分級結(jié)果與實(shí)際較相符，且使用均值變點(diǎn)分析法和累積和分析算法劃分的基本地貌類型與實(shí)際較相符，均值變點(diǎn)分析法和累積和分析算法是相對比較有效的方法，其他方法已不適用，這與前人研究結(jié)果相一致；

(3) 均值變點(diǎn)分析法的結(jié)果與已有研究結(jié)果[13]最接近，研究區(qū)地形起伏度以小起伏為主，說明該方法是計算最佳統(tǒng)計單元的一種最為理想的方法；

(4) 已有研究[13]所用數(shù)據(jù)為30 m分辨率的ASTERG DEM數(shù)據(jù)，與本文所用數(shù)據(jù)類型及數(shù)據(jù)分辨率有所不同，而文獻(xiàn)[13]使用均值變點(diǎn)分析法計算所得的地形起伏度最佳統(tǒng)計單元為17×17(0.260 1 km2)，可見對于同一地區(qū)而言，不同DEM數(shù)據(jù)類型、不同分辨率DEM，所使用的最佳統(tǒng)計單元大小是不一樣的。

本研究分析研究出了適合山西省地形起伏度最佳統(tǒng)計單元的最為理想的計算方法，對均值變點(diǎn)分析法和累積和分析算法計算所得最佳統(tǒng)計單元大小存在差異還需進(jìn)一步的研究。

參考文獻(xiàn):

[1] 劉育成,趙廷寧.基于變點(diǎn)分析法提取廢棄采石場地形起伏度的方法[J].水土保持研究,2016,23(3):269-273.

[2] 涂漢明,劉振東.中國地勢起伏度研究[J].測繪學(xué)報,1991,20(4):311-319.

[3] 劉新華,楊勤科,湯國安.中國地形起伏度的提取及在水土流失定量評價中的應(yīng)用[J].水土保持通報,2001,21(1):57-59,62.

[4] 馬曉微,楊勤科,劉寶元.基于GIS的中國潛在水土流失評價研究[J].水土保持學(xué)報,2002,16(4):49-53.

[5] 閆滿存,李華梅,王光謙.廣東沿海陸地地質(zhì)環(huán)境質(zhì)量定量評價研究[J].工程地質(zhì)學(xué)報,2000,8(4):416-425.

[6] 祁元,劉勇,楊正華,等.基于GIS的蘭州滑坡與泥石流災(zāi)害危險性分析[J].冰川凍土,2012,34(1):96-104.

[7] 劉康,歐陽志云,王效科,等.甘肅省生態(tài)環(huán)境敏感性評價及其空間分布[J].生態(tài)學(xué)報,2003,23(12):2711-2718.

[8] 凡非得,王克林,熊鷹,等.西南喀斯特區(qū)域水土流失敏感性評價及其空間分異特征[J].生態(tài)學(xué)報,2011,31(21):6353-6362.

[9] 徐天獻(xiàn),王玉寬,傅斌.汶川地震重災(zāi)區(qū)土壤侵蝕敏感性評價[J].中國水土保持,2011,18(1):39-42.

[10] 趙斌濱,程永鋒,丁士君,等.基于SRTM-DEM的我國地勢起伏度統(tǒng)計單元研究[J].水利學(xué)報,2015,46(S1):284-290.

[11] 陳學(xué)兄,常慶瑞,郭碧云,等.基于SRTM DEM數(shù)據(jù)的中國地形起伏度分析研究[J].應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報,2013,21(4):670-678.

[12] 王志恒,胡卓瑋,趙文吉,等.應(yīng)用累積和分析算法的地形起伏度最佳統(tǒng)計單元確定[J].測繪科學(xué),2014,39(6):59-64.

[13] 陳學(xué)兄,畢如田,劉正春,等.基于ASTER GDEM數(shù)據(jù)的山西地形起伏度分析研究[J].山西農(nóng)業(yè)大學(xué)學(xué)報:自然科學(xué)版,2016,36(6):417-421.

[14] 中國科學(xué)院地理研究所.中國1:1000000地貌圖制圖規(guī)范(試行)[M].北京:科學(xué)出版社,1987:33-34.

[15] 李炳元,潘保田,韓嘉福.中國陸地基本地貌類型及其劃分指標(biāo)探討[J].第四紀(jì)研究,2008,28(4):535-543.