基于變點分析法提取廢棄采石場地形起伏度的方法

劉育成, 趙廷寧

(1.大連工業(yè)大學(xué), 遼寧 大連 116034; 2.北京林業(yè)大學(xué) 水土保持與荒漠化防治教育部重點實驗室, 北京 100083)

?

基于變點分析法提取廢棄采石場地形起伏度的方法

劉育成1, 趙廷寧2

(1.大連工業(yè)大學(xué), 遼寧 大連 116034; 2.北京林業(yè)大學(xué) 水土保持與荒漠化防治教育部重點實驗室, 北京 100083)

地形起伏度是廢棄采石場水土保持治理和生態(tài)重建的重要地形指標(biāo)之一。以0.5 m，0.7 m，1 m，1.5 m等12種不同水平分辨率DEM為數(shù)據(jù)源，采用GIS的窗口遞增分析法和均值變點分析法對12種不同分辨率DEM的地形起伏度的最佳分析窗口進(jìn)行了分析。結(jié)果表明：12種不同分辨率的最佳分析窗口都是為9×9網(wǎng)格，相應(yīng)地提取地形起伏度的最佳統(tǒng)計面積與其相應(yīng)的水平分辨率關(guān)系比較密切，呈冪函數(shù)的關(guān)系。以0.5 m水平分辨率的DEM數(shù)據(jù)提取地形起伏度的最佳分析窗口面積為20.25 m2，并計算出大于5 m地形起伏度的區(qū)域占總面積的13.12%。均值變點分析法確定地形起伏度也同樣適用于廢棄采石場的地形分析。

廢棄采石場; DEM; 地形起伏度; 最佳統(tǒng)計單元

地形起伏度是宏觀描述地形變化的整體趨勢，也是定量描述地貌形態(tài)的重要指標(biāo)之一[1]。近年來，利用數(shù)字高程模型(DEM)進(jìn)行地形起伏度研究在各個學(xué)術(shù)應(yīng)用流域越來越受到重視，側(cè)重于地形起伏度對區(qū)域地質(zhì)災(zāi)害評價和水土流失定量評價[2-3]、地形起伏度對地區(qū)人口密度與人口分布影響[4-5]、地形起伏度的不同算法對比[6-7]、地形起伏度最佳分析區(qū)域及適宜計算尺度研究[8-10]、地形起伏度與景觀空間格局特征間關(guān)系[11]等。而涉及廢棄礦區(qū)這種特殊土地利用類型的地形起伏度研究相對較少。目前，地形破碎、地勢高差變化大的廢棄礦區(qū)生態(tài)恢復(fù)中有關(guān)水土流失評價中的地形因子比較單一，側(cè)重于坡度、坡長、坡向、海拔高度、溝道密度、流域面積等[12-14]，但國內(nèi)外還沒有一套完整的地形指標(biāo)體系來對廢棄礦區(qū)的水土流失評價。由于DEM提取的坡度已失去原來地貌學(xué)的意義[15]，區(qū)域性的地形起伏度在研究中能更好地分析地形地貌和小流域形態(tài)特征，地形起伏度分析與區(qū)域性滑坡發(fā)生和分布存在良好的相關(guān)性，并且與坡度分析是互相補(bǔ)充的[2]。地形起伏度作為廢棄礦區(qū)水土流失評價的地形因子，能夠反映其土壤侵蝕特征。由于不同的研究目的和應(yīng)用領(lǐng)域，地形起伏度在提取方法和表達(dá)方式等方面有一定的差異性。為此，本研究對地質(zhì)災(zāi)害頻繁發(fā)生和大量水土流失嚴(yán)重的廢棄采石場地形起伏度進(jìn)行提取和計算，并對其不同分辨率DEM的起伏度進(jìn)行對比分析，研究結(jié)果不僅可為廢棄礦區(qū)的水土流失評價提供準(zhǔn)確依據(jù)，還對廢棄礦區(qū)的生態(tài)恢復(fù)和地質(zhì)災(zāi)害評價中具有重要的指導(dǎo)意義。

1　研究區(qū)概況

本研究區(qū)是太行山北段東麓的典型廢棄采石場，位于北京市房山區(qū)(以北緯39°40′23′′、東經(jīng)115°42′35′′為中心)，面積約為39.6 hm2。屬于暖溫帶半濕潤、半干旱大陸性季風(fēng)氣候，四季分明，年平均氣溫4～11.7℃；年均降水量分布不均，主要集中在6—8月份，多年平均降水量為655 mm。地貌類型是以華北低山丘陵區(qū)的石灰?guī)r為主，土壤是多礫石砂壤土，大多是采石場廢棄渣土。天然植被稀疏，植被覆蓋率低，主要有狗尾草(SetariaviridisL．Beauv.)、高羊茅(Festucaarundinacea)、荊條(HibiscussyriacusLinn.)、委陵菜(PotentillaaiscolorBunge)等。

2　研究方法

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

研究區(qū)數(shù)據(jù)是采用Topcon影像全站儀測量所得到的點數(shù)據(jù)，通過ArcGIS 10將這些點數(shù)據(jù)導(dǎo)入，在Beijing1954投影坐標(biāo)系下將其轉(zhuǎn)換為不規(guī)則三角網(wǎng)TIN格式，通過與現(xiàn)場的比對和校正后再將TIN格式內(nèi)插生成規(guī)則柵格DEM。為了對比廢棄采石場不同分辨率的地形起伏度，因此，再將TIN格式轉(zhuǎn)換成0.5 m，0.7 m，1.0 m，1.5 m，2.0 m，2.5 m，3.0 m，3.5 m，4.0 m，4.5 m，5.0 m和7.0 m共12種水平分辨率的柵格DEM。

2.2數(shù)據(jù)方法

地形起伏度算法有好多種[16]，常用的、計算方便的和比較流行的算法是局部高差法。利用ArcGIS 10的空間分析模塊的鄰域分析(Forcalrange函數(shù))來計算地形起伏度，首先采用柵格窗口遞增分析方法，對提取不同分析窗口的地形起伏度進(jìn)行統(tǒng)計分析確定最佳分析窗口，最后根據(jù)最佳分析窗口面積計算出相應(yīng)起伏度。其中，最佳分析窗口面積大小是研究地形起伏度的關(guān)鍵因子[8]。因此，本文采用統(tǒng)計學(xué)上的均值變點分析法來確定地形起伏度隨面積變化的Logarithmic曲線上由陡變緩的點位。均值變點分析法的計算可以按照以下步驟來進(jìn)行[17]：

(1)

(2)

(3)

(2) 根據(jù)式(1)，(2)和(3)計算期望值：E(S-Si)，i=1，2，3，…，N。原始樣本的統(tǒng)計量S與樣本分段后的統(tǒng)計量Si之間最大差值對應(yīng)的點稱為變點。本文所用樣本是指不同分析窗口面積下對應(yīng)的地形平均起伏度值。

3　結(jié)果與分析

3.1不同分辨率下小流域形態(tài)的起伏特征對比

從小流域地貌形態(tài)特征來看，高程和坡度是影響其形態(tài)要素起伏特征的關(guān)鍵因子[18]，也是基本反映了小流域內(nèi)地貌形態(tài)在垂直方向的高差和徑流的過程及變化。因此，對0.5 m，0.7 m，1 m，1.5 m等12種不同分辨率下的小流域面積、高程平均值、高程標(biāo)準(zhǔn)差和平均坡度進(jìn)行計算統(tǒng)計分析，其結(jié)果見表1。從表1可看出，隨著柵格尺寸的逐漸增大，流域面積呈現(xiàn)出降低的趨勢，其變化范圍為0.57%～2.66%；不同柵格尺寸提取的小流域面積平均值為36.49 hm2，標(biāo)準(zhǔn)差1.76。從表1中的高程平均值來看，不同分辨率下的高程平均值變化不顯著；隨著柵格尺寸增大，高程標(biāo)準(zhǔn)差變化總體也是不明顯，變化范圍為0%～0.06%。對于不同分辨率下的坡度提取都是在3×3窗口中進(jìn)行的，隨著柵格尺寸的變大，從表1得出平均坡度呈現(xiàn)遞減趨勢，其算術(shù)平均值為23.84°，標(biāo)準(zhǔn)差為1.13，最大偏差為0.93°，變化范圍為0.28%～4.10%。分析以上產(chǎn)生原因，主要是廢棄采石場地形復(fù)雜多變，隨著柵格尺寸的變大而分辨率的降低，小流域的一些地形特征被平坦化，導(dǎo)致提取的小流域面積、高程平均值和平均坡度隨著柵格尺寸變化而變化。

3.2不同分辨率下地形平均起伏度與網(wǎng)格單元面積關(guān)系

本文采用ArcGIS柵格窗口遞增分析方法，對0.5 m，0.7 m，1 m，1.5 m等12種不同分辨率下研究區(qū)的DEM，采用矩形窗口分析并依次提取2×2，3×3，4×4，…，22×22網(wǎng)格下地形起伏度，并統(tǒng)計出其最佳分析窗口面積。以分辨率為0.5 m的DEM數(shù)據(jù)為例，提取不同窗口下的地形平均起伏度，見表2。

從表2可以得出，隨著網(wǎng)格尺寸增大，地形平均起伏度值也逐漸增大。對地形平均起伏度和分析窗口面積進(jìn)行分析，發(fā)現(xiàn)二者之間存在對數(shù)函數(shù)關(guān)系，得出擬合曲線為：

y=1.4819ln(x)-0.7759(R2=0.9351)

(4)

從式(4)來看，0.5 m分辨率下提取地形起伏度與分析窗口面積之間有良好的擬合效果，通過統(tǒng)計學(xué)的檢驗。

同理，對其他11種不同分辨率下DEM數(shù)據(jù)，分別提取不同分析窗口下的地形起伏度，并對地形平均起伏度與分析窗口面積進(jìn)行分析，也都存在著對數(shù)函數(shù)關(guān)系，并都通過統(tǒng)計學(xué)檢驗，其相應(yīng)曲線函數(shù)擬合效果良好，如圖1所示。

表1　不同分辨率提取DEM的小流域面積、高程和平均坡度

表2　分辨率為0.5 m的地形平均起伏度與柵格單元對應(yīng)關(guān)系

圖1　不同分辨率下網(wǎng)格單元面積與地形平均起伏度對應(yīng)關(guān)系擬合曲線

3.3均值變點分析法下不同分辨率的最佳分析面積

根據(jù)式(1)，(2) 和(3) 分別計算出12種不同分辨率下的分段統(tǒng)計量S和Si，通過S-Si的值確定出陡緩變化的拐點，即為最佳分析面積所對應(yīng)的網(wǎng)格窗口的點。以分辨率為0.5 m的DEM數(shù)據(jù)為例，對表2數(shù)據(jù)計算出單位面積的地形平均起伏度(表3)，并取相應(yīng)對數(shù)構(gòu)建新樣本序列X(X={xi，i=2，3，4，…，22}。按照均值變點分析方法，利用式(1)，(2)和(3)統(tǒng)計出S的值7.675，Si的值見表3。從表3中i和S-Si差值的變化得出曲線，如圖2(左上面第一個)所示。從圖2得出9×9網(wǎng)格為最佳分析窗口，因此本研究區(qū)的分辨率為0.5 m時，地形起伏度的最佳分析面積為20.25 m2。

同理，根據(jù)圖1中數(shù)據(jù)對其他11種分辨率下的最佳分析面積進(jìn)行統(tǒng)計分析，S與Si差值的變化擬合曲線見圖2。從圖2明顯看出，0.5 m，0.7 m，1 m，1.5 m等12種不同分辨率下的S和Si差值達(dá)到最大的點，都是在第8個點，其對應(yīng)的分析窗口為9×9網(wǎng)格。因此廢棄研究區(qū)提取地形起伏度的最佳統(tǒng)計面積與其相應(yīng)的分辨率關(guān)系比較密切，呈冪函數(shù)的關(guān)系。

表3　分辨率為0.5 m的均值變點分析的統(tǒng)計結(jié)果

注：i為變點數(shù)；Ri單位面積上的地形平均起伏度；Si樣本分段后統(tǒng)計量。

圖2　不同分辨率下S和Si差值的變化曲線

根據(jù)表1和表3的結(jié)果，研究區(qū)地形起伏度研究采用分辨率為0.5 m的DEM數(shù)據(jù)、最佳分析面積為20.25 m2，其具體分級見圖3。從圖3計算出地形起伏度，大于5 m的所占百分比為13.12%,是水土流失最嚴(yán)重的區(qū)域，也是后期廢棄采石場土地復(fù)墾治理的重點區(qū)域；在3～5 m之間所占百分比為16.52%，在1～3 m之間所占百分比為51.01%，在0～1 m范圍內(nèi)所占總面積的19.35%。根據(jù)不同地形起伏度和廢棄采石場復(fù)綠規(guī)劃，采取相應(yīng)的水土保持措施。

圖3　分辨率為0.5 m的地形起伏度分級

4　結(jié) 論

(1) 對0.5 m，0.7 m，1 m，1.5 m等12種不同水平分辨率下小流域形態(tài)要素的起伏特征如流域面積、平均高程值、程標(biāo)準(zhǔn)差和平均坡度等進(jìn)行計算和比較。小流域面積隨著柵格尺寸增大而呈現(xiàn)出降低的趨勢，變化范圍為0.57%～2.66%。平均高程對12種不同分辨率沒有顯著影響。平均坡度隨著柵格尺寸增大而逐漸降低，變化范圍為0.28%～4.10%。

(2) 0.5 m，0.7 m，1 m，1.5 m等12種不同分辨率下的最佳分析窗口是9×9網(wǎng)格，相應(yīng)地提取地形起伏度的最佳統(tǒng)計面積與其相應(yīng)的分辨率關(guān)系比較密切，呈冪函數(shù)的關(guān)系。

(3) 對廢棄采石場的分辨率為0.5 m的DEM數(shù)據(jù)進(jìn)行提取地形起伏度，其統(tǒng)計最佳分析面積為20.25 m2，統(tǒng)計出地形起伏度大于5 m的區(qū)域占總面積的13.12%，也是水土流失治理重點區(qū)域；在3～5 m

之間所占百分比為16.52%，在1～3 m之間所占百分比為51.01%。

(4) 地形起伏是導(dǎo)致水土流失的最直接原因，因此針對廢棄礦山生態(tài)修復(fù)或工程綠化規(guī)劃和設(shè)計常用的水平分別率，提出比較或最適宜分析地形起伏度的方法。

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

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

[3]郭芳芳,楊農(nóng),孟暉,等.地形起伏度和坡度分析在區(qū)域滑坡災(zāi)害評價中的應(yīng)用[J].中國地質(zhì),2008,35(1):131-143.

[4]封志明,唐焰,楊艷昭,等.中國地形起伏度及其與人口分布的相關(guān)性[J].地理學(xué)報,2007,62(10):1073-1082.

[5]周自翔,李晶,任志遠(yuǎn).基于GIS的關(guān)中—天水經(jīng)濟(jì)區(qū)地形起伏度與人口分布研究[J].地理科學(xué),2012,32(8):951-957.

[6]朗玲玲,程維明,朱啟疆,等.多尺度DEM提取地勢起伏度的對比分析:以福建低山丘陵區(qū)為例[J].地球信息科學(xué),2007,9(6):1-6.

[7]蔣好忱,楊勤科.基于DEM的地形起伏度算法的比較研究[J].水土保持通報,2014,34(6):162-166.

[8]涂漢明,劉振東.中國地勢起伏度最佳統(tǒng)計單元的求證[J].湖北大學(xué)學(xué)報:自然科學(xué)版,1990,12(3):266-271.

[9]張錦明,游雄.地形起伏度最佳分析區(qū)域研究[J].測繪科學(xué)技術(shù)學(xué)報,2011,28(5):369-373.

[10]陳學(xué)兄,常慶瑞,郭碧云,等.基于STRM DEM數(shù)據(jù)的中國地形起伏度分析研究[J].測繪科學(xué)技術(shù)學(xué)報,2013,21(4):670-678.

[11]趙衛(wèi)權(quán),蘇維詞,袁俊.基于地形起伏度的貴州省景觀空間格局分異特征[J].水土保持研究,2010,17(2):105-110.

[12]楊翠霞.露天開采礦區(qū)廢棄地近自然地形重塑研究[D].北京:北京林業(yè)大學(xué),2014.

[13]呂春娟,白中科,陳衛(wèi)國,等.黃土區(qū)露天礦排土場水分調(diào)控技術(shù)研究[J].水土保持通報,2011,31(1):160-164.

[14]楊翠霞,趙廷寧,謝寶元,等.基于流域自然形態(tài)的廢棄礦山地形恢復(fù)模擬[J].農(nóng)業(yè)工程學(xué)報,2014,30(1):236-244.

[15]Tomislav H, Hannes I, Reuter. Geomorphometry Concepts, Software, Applications[M]. Amsterdam: Elsevier,2009.

[16]Tang Guoan. A Reseach on the Accuracy of Digital Elevation Models [M]. Beijing: Science Press,2000.

[17]項靜恬,史久恩.非線性系統(tǒng)中數(shù)據(jù)處理的統(tǒng)計方法[M]:北京:科學(xué)出版社,2000.

[18]承繼成,江美球.流域地貌數(shù)學(xué)模型[M].北京:科學(xué)出版社,1986.

Method for Extraction of Relief Amplitude of Abandoned Quarry Based on Change Point Method

LIU Yucheng1, ZHAO Tingning2

(1.DalianPolytechnicUniversity,Dalian,Liaoning116034,China; 2.KeyLabofSoilandWaterConservationandDesertificationCombatingofMinistryofEducation,CollegeofSoilandWaterConservation,BeijingForestryUniversity,Beijing100083,China)

The relief amplitude of abandoned quarry is one of the important terrain indexes of soil and water conservation management and ecological restoration in mining area. Twelve different levels of DEMs such as 0.5 m, 0.7 m, 1.0 m, 1.5 m and so on, were taken as the data sources. GIS′s grid window increment analysis method and change point method were used for the analysis of the twelve different levels of DEMs to extract relief amplitude. The results showed that the optimum statistical unit was 9×9 grid. The relationship between the optimum statistical unit area of the relief amplitude and the corresponding resolution was relatively close, which could be described by a power function. The optimum statistical unit area of 0.5 m DEM relief amplitude was 20.25 m2. The relief amplitude greater than 5 m accounted for 13.12% of the total area. The change point method was used to determine the optimum statistical unit of the relief amplitude, which is also applicable to the terrain analysis of the abandoned quarry.

abandoned quarry; DEM; relief amplitude; optimum statistical unit

2015-10-30

2015-11-18

林業(yè)公益性行業(yè)科研專項“建設(shè)工程損毀林地植被修復(fù)關(guān)鍵技術(shù)研究與示范”(200904030);大連工業(yè)大學(xué)青年基金(QNJJ-201423)

劉育成(1975—),男,陜西省扶風(fēng)人,碩士,講師,研究方向為環(huán)境設(shè)計及數(shù)字化景觀研究。E-mail:dllyc@qq.com

趙廷寧(1964—),男,河北省陽原人,博士生導(dǎo)師,教授,主要從事水土保持與工程綠化研究。E-mail:zhtning@bjfu.edu.cn

TD88

A

1005-3409(2016)03-0269-05