斜坡地質災害敏感性評價中地勢起伏度提取最佳尺度研究

張明媚，薛永安

(1.太原理工大學 礦業工程學院，太原 030024；2.山西能源學院 地質測繪工程系，山西 晉中 030600)

斜坡地質災害敏感性是指在特殊地形或某些因素作用下發生滑坡、崩塌災害的可能性[1]，地質災害敏感性評價的本質就是利用數學語言來表達在給定的地質環境條件下地質災害在空間上的發生概率[2-3]。目前，對于地質災害敏感性評價研究已有許多成果。ACHOUR et al[4]基于層次分析法和信息量模型開展了滑坡敏感性研究。YALCIN et al[5]的對比研究表明層次分析法對滑坡敏感性的分區與實際情況更加符合。NICU et al[6]采用坡度、高程、曲率、巖性、降水、土地利用、地形潤濕性指數、地形、坡向、平面曲率和距離河流距離等11個條件因子對滑坡的敏感性程度進行了統計分析。HONG et al[7]選取坡度、坡向、海拔、平面曲率、剖面曲率、巖性、距離斷層距離、距離河流距離、距離道路距離、土地利用、歸一化植被指數等15個因子作為滑坡災害敏感性評價因子集合。THAI et al[8]選用坡度、曲率、海拔等15個影響因子作為滑坡災害敏感性評價因子集合。而KALANTAR et al[9]對Dodangeh流域進行滑坡災害敏感性評價研究中采用的因子為曲率、平面曲率、剖面曲率、高程、坡度、坡向、距離斷層距離、距離河流距離、地形濕度指數、河流水力指數、地形粗糙度指數、沉積物搬運指數、巖性及土地利用。在上述這些研究中，以高程、坡度、坡向、曲率等為代表的數字地形因子是斜坡地質災害敏感性評價的重要因子。

然而，作為描述一個區域地形特征的宏觀性指標，地勢起伏度指在一個特定的區域內，最高點海拔高度與最低點海拔高度的差值，它能反映地表起伏變化，與斜坡地質災害的發生存在很大的相關性[10]。按照地貌發育的基本理論，一種地貌類型存在一個使最大高差達到相對穩定的最佳分析區域，傳統方法是利用柵格窗口的遞增方法來尋找最佳的分析窗口[11]。這是地貌分析中地勢起伏度計算的最佳統計單元思路，即隨著統計單元半徑的增大，地勢起伏度隨之變化，到達某一臨界點后趨于穩定，臨界點即為地勢起伏度的最佳統計單元，也是真實反映地形起伏的DEM地勢起伏度提取單位面積。針對地勢起伏度最佳提取單元的研究[12-13]，眾多學者以不同尺度的DEM在不同實驗區開展了研究工作，結果各異。涂漢明等[14]通過對全國600個樣點和2個小區的詳細研究，運用模糊數學方法，確定了全國地勢起伏度最佳統計單元為21 km2.張軍等[15]則得出基于1∶250 000 DEM的新疆地勢起伏度計算的最佳統計單元為2.56 km2.趙斌濱等[16]基于90 m分辨率SRTM-DEM數據提取得到中國地勢起伏度的最佳統計單元面積為2.25 km2.王讓虎等[17]基于ASTER GDEM數據提取得到中國東北地區地形起伏度的最佳統計單元面積為2.62 km2.陳學兄等[18]以30 m分辨率ASTER GDEM數據為基礎，運用均值變點分析法計算得到0.260 1 km2為山西地形起伏度提取的最佳統計單元。楊艷林等[19]基于ASTER GDEM數據分析得出長江中游的地形起伏度最佳分析窗口為19×19像元，面積約0.324 9 km2.而郭芳芳等[20]從地勢起伏度分析在區域斜坡災害評價中的應用出發，尋找反映斜坡點對應的真實地勢起伏度單元，對不同統計單元地形起伏度值與斜坡個數統計曲線峰值分布特征進行分析，確定了研究區地形起伏度最佳統計單元面積為4 km2.

基于上述研究，本文以山西省太原市西山地質塊體為研究區，以ASTER GDEM V2為地勢起伏度提取數據源，以2012年斜坡地質災害分布信息為斜坡災害發育敏感性評價數據源，以均值變點法開展斜坡災害敏感性評價中地勢起伏度因子最佳提取單元研究。

1 數據與方法

1.1 研究區與數據源

研究區位于山西省中部、太原盆地的北端，北接太原市陽曲縣，西鄰太原市古交市、清徐縣，東鄰太原盆地，南靠太原市清徐縣，包括太原市尖草坪區、萬柏林區和晉源區，總面積441.063 km2，見圖1.研究區的地理坐標為東經112°15′03″-112°29′25″，北緯37°38′54″-38°01′14″.

圖1 研究區地理位置圖Fig.1 Geographical map of the studied area

1.1.1DEM數據

ASTER GDEM(先進星載熱輻射和反射儀全球數字高程模型)是根據NASA的新一代對地觀測衛星Terra的詳盡觀測結果制作完成的，其全球空間分辨率約為30 m，垂直分辨率為20 m，空間參考為WGS84/EGM96[21].ASTER GDEM目前共有兩版數據，本文研究采用第2版，于2011年10月發布，數據來源于中國科學院計算機網絡信息中心地理空間數據云平臺。以研究區ASTER GDEM V2重采樣后30 m分辨率DEM為基準數據(數據參數及精度見表1).

表1 研究區主要地形參數及信息源精度Table 1 Main terrain parameters and information source accuracy of the studied area

1.1.2斜坡地質災害數據

從太原市規劃和自然資源局收集到研究區2012年斜坡地質災害分布信息，見圖2.

圖2 研究區斜坡地質災害空間分布圖Fig.2 Spatial distribution of slope geological hazards in the studied area

1.2 均值變點法

目前，最佳地勢起伏度統計單元的尋找主要依賴人工判斷，以平均地勢起伏度隨統計單元面積變化的轉折點作為最佳統計單元面積，也就是人工判斷曲線上由陡變緩的位置；這種判斷方法顯然主觀性太強，不夠準確。隨著最佳地勢起伏度統計單元研究工作的深入[22]，統計學領域的均值變點法逐漸被用作計算曲線上由陡變緩點位的客觀方法。

均值變點法的原理及步驟如下[23]：

1) 檢驗是否存在變點，即檢驗原假設H0是否存在變點。若H0被接受，則該數據序列沒有變點；若H0被否定，則假設該序列中至多存在q個變點，對m1，m2，…，mq變點進行估計。

2) 估計變點個數。

3) 估計變點處均值的跳躍度。

在實際估計中，如果先驗知識足以肯定序列中變點的存在，則可跳過步驟1)而直接進入后幾步。

均值變點問題的離散數據模型如下。

設xt=at+et,t=1,2,…,N，則

a1=…=am1-1=b1，am1=…=am2-1=b2，…，amq=…=aN=bq+1.

(1)

式中：1

H0∶b1=b2=…=bq+1.

(2)

此處特別強調本檢驗與多樣本檢驗的差別，即此檢驗中m1,m2,…,mq為未知。

原假設H0的檢驗步驟如下。

1) 令i=2,…,N，每個i將樣本分為兩段：

x1,x2,…,xi-1和xi,xi+1,…,xN.

(3)

(4)

2) 計算統計量：

(5)

3) 計算期望值：

E(S-Si),i=2,3,…,N.

(6)

4) 求極大值：

(7)

在平均意義下認為：

S*=min(S2,…,SN) .

(8)

5) 取檢驗顯著性水平為α，計算C值：

C=σ2(2lg lgN+lg lg lgN-lg π-2lg(-0.5lg(1-α))) .

(9)

如果σ2未知，則用下面估計來替代：

(10)

6) 如果S-S*>C，則否定H0，認為無變點，否則接受H0，該檢驗有漸近水平α.

平均來說，如果序列中存在變點，則S和Si的差距會因變點的存在而增大，所以均值變點法原理明確，操作簡單。同時，均值變點法對只有一個變點的檢驗最為有效，存在多個變點時有可能因為均值的多次升降而抵消了S和Si之間的差距。

2 結果與討論

2.1 平均地勢起伏度與統計單元之間的關系

基于ArcGIS平臺，利用鄰域分析法提取研究區不同面積單元的地勢起伏度，提取單元分別為30 m分辨率DEM數據網格的2×2，3×3，4×4，5×5，…，25×25，對應的提取單元面積(km2)分別為：0.003 6,0.008 1,0.014 4,0.022 5,0.032 4,0.044 1,0.057 6,0.072 9,0.090 0,0.108 9,0.129 6,0.152 1,0.176 4,0.202 5,0.230 4,0.260 1,0.291 6,0.324 9,0.360 0,0.396 9,0.435 6,0.476 1,0.518 4,0.562 5，共24個面積單元。經統計，研究區不同提取單元大小與平均地勢起伏度之間關系如表2所示。

由表2可以看出，隨著提取的地勢起伏度統計單元的增大，研究區平均地勢起伏度呈現上升趨勢，從最小的13.53 m上升到173.96 m，這有助于分析研究區斜坡地質災害發育與地形起伏之間的關系。至于哪個尺度的分析窗口是最佳窗口，可以通過區域地勢起伏度與斜坡地質災害發育之間的關系進行客觀的分析評價。

表2 研究區不同統計單元與平均地勢起伏度之間的關系Table 2 Relationship between statistical unit and average relief amplitude of the area

以統計分析單元面積為橫坐標，相應的平均地勢起伏度計算值為縱坐標，對地勢起伏度隨統計單元面積的變化情況進行擬合，結果如圖3所示。

對于不同統計單元面積與所提取地勢起伏度之間的關系，通過對數方程進行擬合，得到擬合曲線方程。通過統計學檢驗，擬合效果良好，擬合方程為：

y=34.097lnx-52.50,R2=0.972 9 .

(11)

式中：y為平均地勢起伏度，m；x為統計單元面積，km2；R2為決定系數。

從圖3的曲線變化趨勢來看，應有1個變點存在，過了變點后地勢起伏度的增速變緩，而這個變點所對應的面積就是研究區提取地勢起伏度的最佳統計單元面積。

圖3 地勢起伏度與統計單元面積對應關系擬合曲線Fig.3 Fitting curve of the corresponding relation between relief amplitude and area of statistical unit

2.2 地勢起伏度提取最佳統計單元分析

對不同統計單元的平均地勢起伏度、不同統計單元地勢起伏度分級中的斜坡地質災害峰值、地勢起伏度峰值進行統計，結果見圖4.

圖4 平均地勢起伏度、地勢起伏度峰值、斜坡地質災害峰值與統計單元關系圖Fig.4 Relationship between average value of relief amplitude, peak value of relief amplitude, and peak value of slope geological hazards as one side and statistical unit as another side

根據圖4統計結果分別構建樣本序列，利用均值變點法的公式編制Excel程序，計算樣本序列的統計量S和Si，求得S-Si(見表3)，構建S與Si差值的變化擬合曲線。

表3 均值變點法分析統計結果Table 3 Statistical results of mean change point

平均地勢起伏度樣本序列S值為16.994 5，S與Si差值的變化曲線見圖5.地勢起伏度峰值樣本序列S的值為20.130 4，S與Si差值的變化曲線見圖6.斜坡地質災害峰值樣本序列S的值為89.359 3，S與Si差值的變化曲線見圖7.

圖5 基于平均地勢起伏度的S-Si變化曲線Fig.5 Variation curve based on average value of relief amplitude S-Si

圖6 基于地勢起伏度峰值的S-Si變化曲線Fig.6 Variation curve based on peak value of relief amplitude S-Si

圖7 基于斜坡地質災害峰值的S-Si變化曲線Fig.7 Variation curve based on peak value of slope geological hazards S-Si

2.3 討論

1) 圖4顯示，每一種統計單元都存在一個斜坡地質災害發育的峰值區間，隨統計單元的增大，峰值逐漸減小，到達某一個臨界值后峰值衰減趨于穩定。

2) 圖5顯示，i為11時，S與Si差值達最大，i=11處即為變點，對應統計單元為12×12網格。因此，由平均地勢起伏度以均值變點法來分析得到研究區的最佳地勢起伏度提取單元為12×12網格，最佳統計面積為0.129 6 km2.

3) 圖6顯示，i為11時，S與Si差值達最大，i=11處即為變點，對應統計單元為12×12網格。因此，由地勢起伏度峰值以均值變點法來分析得到研究區的最佳地勢起伏度提取單元為12×12網格，最佳統計面積為0.129 6 km2.

4) 圖7顯示，i為8時，S與Si差值達最大，i=8處即為變點，對應統計單元為9×9網格。因此，通過斜坡地質災害峰值以均值變點法來分析得到研究區的最佳地勢起伏度提取單元為9×9網格，最佳統計面積為0.072 9 km2.

5) 對比三種樣本序列所計算得到的最佳地勢起伏度提取單元，使用平均地勢起伏度和地勢起伏度峰值樣本序列計算結果一致，均為12×12網格，而斜坡地質災害峰值樣本序列計算結果則為9×9網格。作為客觀存在的區域地貌特征因子，斜坡地質災害演變頻繁，而平均地勢起伏度演變較小。因此，以研究區數字地貌角度分析，基于平均地勢起伏度所計算的最佳提取單元更合理；而以地質災害發育角度分析，基于斜坡地質災害峰值所計算的最佳提取單元則更合理。

6) 針對斜坡地質災害敏感性評價中地勢起伏度最佳統計單元開展研究，重點關注的是研究區地勢起伏度與斜坡地質災害發育之間的關系，因此所尋找的是能反映地勢起伏度與斜坡地質災害發育之間關系的最佳統計單元，而不是數字地貌研究中區域地勢起伏度計算的最佳統計單元。因此，選擇9×9網格作為研究區斜坡地質災害敏感性評價中地勢起伏度因子最佳提取單元。

3 結論

1) 地勢起伏度與統計單元面積之間存在顯著的對數變化關系，研究區內地勢起伏度隨統計單元面積的增大而迅速增大，當統計單元面積達到一定數值后，增大速度逐漸變緩，最后近似趨于平穩，這為開展地勢起伏度最佳提取單元提供了基礎。

2) 地勢起伏度提取最佳統計單元與地質災害分布峰值、平均地勢起伏度和地勢起伏度峰值之間存在明顯的相關性，目前的研究主要基于平均地勢起伏度，利用均值變點法或主觀方式分析確定最佳提取單元，這給斜坡地質災害敏感性評價帶來一定的不確定性，所以應綜合三種情況開展統計并綜合分析后確定研究區地勢起伏度最佳統計尺度。

3) 以研究區ASTER GDEM重采樣30 m分辨率DEM和2012年斜坡類地質災害空間分布信息為數據源，以2×2，3×3，…，25×25共24個分析窗口提取了研究區的地勢起伏度，以均值變點法對不同統計單元的平均地勢起伏度、不同統計單元地勢起伏度分級中的斜坡地質災害峰值、地勢起伏度峰值進行統計分析。綜合三種樣本序列計算結果，本文選擇9×9網格作為研究區斜坡類地質災害敏感性評價時的最佳地勢起伏度提取單元，這為研究區斜坡類地質災害敏感性評價研究提供了可靠的地勢起伏度因子，同時也為同類研究提供了數字地形因子較為準確的提取方法。