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

劉育成, 趙廷寧

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

?

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

劉育成1, 趙廷寧2

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

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

廢棄采石場; DEM; 地形起伏度; 最佳統計單元

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

1　研究區概況

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

2　研究方法

2.1數據來源

研究區數據是采用Topcon影像全站儀測量所得到的點數據，通過ArcGIS 10將這些點數據導入，在Beijing1954投影坐標系下將其轉換為不規則三角網TIN格式，通過與現場的比對和校正后再將TIN格式內插生成規則柵格DEM。為了對比廢棄采石場不同分辨率的地形起伏度，因此，再將TIN格式轉換成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數據方法

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

(1)

(2)

(3)

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

3　結果與分析

3.1不同分辨率下小流域形態的起伏特征對比

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

3.2不同分辨率下地形平均起伏度與網格單元面積關系

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

從表2可以得出，隨著網格尺寸增大，地形平均起伏度值也逐漸增大。對地形平均起伏度和分析窗口面積進行分析，發現二者之間存在對數函數關系，得出擬合曲線為：

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

(4)

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

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

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

表2　分辨率為0.5 m的地形平均起伏度與柵格單元對應關系

圖1　不同分辨率下網格單元面積與地形平均起伏度對應關系擬合曲線

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

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

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

表3　分辨率為0.5 m的均值變點分析的統計結果

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

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

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

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

4　結 論

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

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

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

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

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

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

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

[3]郭芳芳,楊農,孟暉,等.地形起伏度和坡度分析在區域滑坡災害評價中的應用[J].中國地質,2008,35(1):131-143.

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

[5]周自翔,李晶,任志遠.基于GIS的關中—天水經濟區地形起伏度與人口分布研究[J].地理科學,2012,32(8):951-957.

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

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

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

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

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

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

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

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

[14]楊翠霞,趙廷寧,謝寶元,等.基于流域自然形態的廢棄礦山地形恢復模擬[J].農業工程學報,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]項靜恬,史久恩.非線性系統中數據處理的統計方法[M]:北京:科學出版社,2000.

[18]承繼成,江美球.流域地貌數學模型[M].北京:科學出版社,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

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

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

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

TD88

A

1005-3409(2016)03-0269-05