張明媚,薛永安
(1.太原理工大學(xué) 礦業(yè)工程學(xué)院,太原 030024;2.山西能源學(xué)院 地質(zhì)測繪工程系,山西 晉中 030600)
斜坡地質(zhì)災(zāi)害敏感性是指在特殊地形或某些因素作用下發(fā)生滑坡、崩塌災(zāi)害的可能性[1],地質(zhì)災(zāi)害敏感性評(píng)價(jià)的本質(zhì)就是利用數(shù)學(xué)語言來表達(dá)在給定的地質(zhì)環(huán)境條件下地質(zhì)災(zāi)害在空間上的發(fā)生概率[2-3]。目前,對(duì)于地質(zhì)災(zāi)害敏感性評(píng)價(jià)研究已有許多成果。ACHOUR et al[4]基于層次分析法和信息量模型開展了滑坡敏感性研究。YALCIN et al[5]的對(duì)比研究表明層次分析法對(duì)滑坡敏感性的分區(qū)與實(shí)際情況更加符合。NICU et al[6]采用坡度、高程、曲率、巖性、降水、土地利用、地形潤濕性指數(shù)、地形、坡向、平面曲率和距離河流距離等11個(gè)條件因子對(duì)滑坡的敏感性程度進(jìn)行了統(tǒng)計(jì)分析。HONG et al[7]選取坡度、坡向、海拔、平面曲率、剖面曲率、巖性、距離斷層距離、距離河流距離、距離道路距離、土地利用、歸一化植被指數(shù)等15個(gè)因子作為滑坡災(zāi)害敏感性評(píng)價(jià)因子集合。THAI et al[8]選用坡度、曲率、海拔等15個(gè)影響因子作為滑坡災(zāi)害敏感性評(píng)價(jià)因子集合。而KALANTAR et al[9]對(duì)Dodangeh流域進(jìn)行滑坡災(zāi)害敏感性評(píng)價(jià)研究中采用的因子為曲率、平面曲率、剖面曲率、高程、坡度、坡向、距離斷層距離、距離河流距離、地形濕度指數(shù)、河流水力指數(shù)、地形粗糙度指數(shù)、沉積物搬運(yùn)指數(shù)、巖性及土地利用。在上述這些研究中,以高程、坡度、坡向、曲率等為代表的數(shù)字地形因子是斜坡地質(zhì)災(zāi)害敏感性評(píng)價(jià)的重要因子。
然而,作為描述一個(gè)區(qū)域地形特征的宏觀性指標(biāo),地勢起伏度指在一個(gè)特定的區(qū)域內(nèi),最高點(diǎn)海拔高度與最低點(diǎn)海拔高度的差值,它能反映地表起伏變化,與斜坡地質(zhì)災(zāi)害的發(fā)生存在很大的相關(guān)性[10]。按照地貌發(fā)育的基本理論,一種地貌類型存在一個(gè)使最大高差達(dá)到相對(duì)穩(wěn)定的最佳分析區(qū)域,傳統(tǒng)方法是利用柵格窗口的遞增方法來尋找最佳的分析窗口[11]。這是地貌分析中地勢起伏度計(jì)算的最佳統(tǒng)計(jì)單元思路,即隨著統(tǒng)計(jì)單元半徑的增大,地勢起伏度隨之變化,到達(dá)某一臨界點(diǎn)后趨于穩(wěn)定,臨界點(diǎn)即為地勢起伏度的最佳統(tǒng)計(jì)單元,也是真實(shí)反映地形起伏的DEM地勢起伏度提取單位面積。針對(duì)地勢起伏度最佳提取單元的研究[12-13],眾多學(xué)者以不同尺度的DEM在不同實(shí)驗(yàn)區(qū)開展了研究工作,結(jié)果各異。涂漢明等[14]通過對(duì)全國600個(gè)樣點(diǎn)和2個(gè)小區(qū)的詳細(xì)研究,運(yùn)用模糊數(shù)學(xué)方法,確定了全國地勢起伏度最佳統(tǒng)計(jì)單元為21 km2.張軍等[15]則得出基于1∶250 000 DEM的新疆地勢起伏度計(jì)算的最佳統(tǒng)計(jì)單元為2.56 km2.趙斌濱等[16]基于90 m分辨率SRTM-DEM數(shù)據(jù)提取得到中國地勢起伏度的最佳統(tǒng)計(jì)單元面積為2.25 km2.王讓虎等[17]基于ASTER GDEM數(shù)據(jù)提取得到中國東北地區(qū)地形起伏度的最佳統(tǒng)計(jì)單元面積為2.62 km2.陳學(xué)兄等[18]以30 m分辨率ASTER GDEM數(shù)據(jù)為基礎(chǔ),運(yùn)用均值變點(diǎn)分析法計(jì)算得到0.260 1 km2為山西地形起伏度提取的最佳統(tǒng)計(jì)單元。楊艷林等[19]基于ASTER GDEM數(shù)據(jù)分析得出長江中游的地形起伏度最佳分析窗口為19×19像元,面積約0.324 9 km2.而郭芳芳等[20]從地勢起伏度分析在區(qū)域斜坡災(zāi)害評(píng)價(jià)中的應(yīng)用出發(fā),尋找反映斜坡點(diǎn)對(duì)應(yīng)的真實(shí)地勢起伏度單元,對(duì)不同統(tǒng)計(jì)單元地形起伏度值與斜坡個(gè)數(shù)統(tǒng)計(jì)曲線峰值分布特征進(jìn)行分析,確定了研究區(qū)地形起伏度最佳統(tǒng)計(jì)單元面積為4 km2.
基于上述研究,本文以山西省太原市西山地質(zhì)塊體為研究區(qū),以ASTER GDEM V2為地勢起伏度提取數(shù)據(jù)源,以2012年斜坡地質(zhì)災(zāi)害分布信息為斜坡災(zāi)害發(fā)育敏感性評(píng)價(jià)數(shù)據(jù)源,以均值變點(diǎn)法開展斜坡災(zāi)害敏感性評(píng)價(jià)中地勢起伏度因子最佳提取單元研究。
研究區(qū)位于山西省中部、太原盆地的北端,北接太原市陽曲縣,西鄰太原市古交市、清徐縣,東鄰太原盆地,南靠太原市清徐縣,包括太原市尖草坪區(qū)、萬柏林區(qū)和晉源區(qū),總面積441.063 km2,見圖1.研究區(qū)的地理坐標(biāo)為東經(jīng)112°15′03″-112°29′25″,北緯37°38′54″-38°01′14″.
圖1 研究區(qū)地理位置圖Fig.1 Geographical map of the studied area
1.1.1DEM數(shù)據(jù)
ASTER GDEM(先進(jìn)星載熱輻射和反射儀全球數(shù)字高程模型)是根據(jù)NASA的新一代對(duì)地觀測衛(wèi)星Terra的詳盡觀測結(jié)果制作完成的,其全球空間分辨率約為30 m,垂直分辨率為20 m,空間參考為WGS84/EGM96[21].ASTER GDEM目前共有兩版數(shù)據(jù),本文研究采用第2版,于2011年10月發(fā)布,數(shù)據(jù)來源于中國科學(xué)院計(jì)算機(jī)網(wǎng)絡(luò)信息中心地理空間數(shù)據(jù)云平臺(tái)。以研究區(qū)ASTER GDEM V2重采樣后30 m分辨率DEM為基準(zhǔn)數(shù)據(jù)(數(shù)據(jù)參數(shù)及精度見表1).
表1 研究區(qū)主要地形參數(shù)及信息源精度Table 1 Main terrain parameters and information source accuracy of the studied area
1.1.2斜坡地質(zhì)災(zāi)害數(shù)據(jù)
從太原市規(guī)劃和自然資源局收集到研究區(qū)2012年斜坡地質(zhì)災(zāi)害分布信息,見圖2.
圖2 研究區(qū)斜坡地質(zhì)災(zāi)害空間分布圖Fig.2 Spatial distribution of slope geological hazards in the studied area
目前,最佳地勢起伏度統(tǒng)計(jì)單元的尋找主要依賴人工判斷,以平均地勢起伏度隨統(tǒng)計(jì)單元面積變化的轉(zhuǎn)折點(diǎn)作為最佳統(tǒng)計(jì)單元面積,也就是人工判斷曲線上由陡變緩的位置;這種判斷方法顯然主觀性太強(qiáng),不夠準(zhǔn)確。隨著最佳地勢起伏度統(tǒng)計(jì)單元研究工作的深入[22],統(tǒng)計(jì)學(xué)領(lǐng)域的均值變點(diǎn)法逐漸被用作計(jì)算曲線上由陡變緩點(diǎn)位的客觀方法。
均值變點(diǎn)法的原理及步驟如下[23]:
1) 檢驗(yàn)是否存在變點(diǎn),即檢驗(yàn)原假設(shè)H0是否存在變點(diǎn)。若H0被接受,則該數(shù)據(jù)序列沒有變點(diǎn);若H0被否定,則假設(shè)該序列中至多存在q個(gè)變點(diǎn),對(duì)m1,m2,…,mq變點(diǎn)進(jìn)行估計(jì)。
2) 估計(jì)變點(diǎn)個(gè)數(shù)。
3) 估計(jì)變點(diǎn)處均值的跳躍度。
在實(shí)際估計(jì)中,如果先驗(yàn)知識(shí)足以肯定序列中變點(diǎn)的存在,則可跳過步驟1)而直接進(jìn)入后幾步。
均值變點(diǎn)問題的離散數(shù)據(jù)模型如下。
設(shè)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) 此處特別強(qiáng)調(diào)本檢驗(yàn)與多樣本檢驗(yàn)的差別,即此檢驗(yàn)中m1,m2,…,mq為未知。 原假設(shè)H0的檢驗(yàn)步驟如下。 1) 令i=2,…,N,每個(gè)i將樣本分為兩段: x1,x2,…,xi-1和xi,xi+1,…,xN. (3) (4) 2) 計(jì)算統(tǒng)計(jì)量: (5) 3) 計(jì)算期望值: E(S-Si),i=2,3,…,N. (6) 4) 求極大值: (7) 在平均意義下認(rèn)為: S*=min(S2,…,SN) . (8) 5) 取檢驗(yàn)顯著性水平為α,計(jì)算C值: C=σ2(2lg lgN+lg lg lgN-lg π-2lg(-0.5lg(1-α))) . (9) 如果σ2未知,則用下面估計(jì)來替代: (10) 6) 如果S-S*>C,則否定H0,認(rèn)為無變點(diǎn),否則接受H0,該檢驗(yàn)有漸近水平α. 平均來說,如果序列中存在變點(diǎn),則S和Si的差距會(huì)因變點(diǎn)的存在而增大,所以均值變點(diǎn)法原理明確,操作簡單。同時(shí),均值變點(diǎn)法對(duì)只有一個(gè)變點(diǎn)的檢驗(yàn)最為有效,存在多個(gè)變點(diǎn)時(shí)有可能因?yàn)榫档亩啻紊刀窒薙和Si之間的差距。 基于ArcGIS平臺(tái),利用鄰域分析法提取研究區(qū)不同面積單元的地勢起伏度,提取單元分別為30 m分辨率DEM數(shù)據(jù)網(wǎng)格的2×2,3×3,4×4,5×5,…,25×25,對(duì)應(yīng)的提取單元面積(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個(gè)面積單元。經(jīng)統(tǒng)計(jì),研究區(qū)不同提取單元大小與平均地勢起伏度之間關(guān)系如表2所示。 由表2可以看出,隨著提取的地勢起伏度統(tǒng)計(jì)單元的增大,研究區(qū)平均地勢起伏度呈現(xiàn)上升趨勢,從最小的13.53 m上升到173.96 m,這有助于分析研究區(qū)斜坡地質(zhì)災(zāi)害發(fā)育與地形起伏之間的關(guān)系。至于哪個(gè)尺度的分析窗口是最佳窗口,可以通過區(qū)域地勢起伏度與斜坡地質(zhì)災(zāi)害發(fā)育之間的關(guān)系進(jìn)行客觀的分析評(píng)價(jià)。 表2 研究區(qū)不同統(tǒng)計(jì)單元與平均地勢起伏度之間的關(guān)系Table 2 Relationship between statistical unit and average relief amplitude of the area 以統(tǒng)計(jì)分析單元面積為橫坐標(biāo),相應(yīng)的平均地勢起伏度計(jì)算值為縱坐標(biāo),對(duì)地勢起伏度隨統(tǒng)計(jì)單元面積的變化情況進(jìn)行擬合,結(jié)果如圖3所示。 對(duì)于不同統(tǒng)計(jì)單元面積與所提取地勢起伏度之間的關(guān)系,通過對(duì)數(shù)方程進(jìn)行擬合,得到擬合曲線方程。通過統(tǒng)計(jì)學(xué)檢驗(yàn),擬合效果良好,擬合方程為: y=34.097lnx-52.50,R2=0.972 9 . (11) 式中:y為平均地勢起伏度,m;x為統(tǒng)計(jì)單元面積,km2;R2為決定系數(shù)。 從圖3的曲線變化趨勢來看,應(yīng)有1個(gè)變點(diǎn)存在,過了變點(diǎn)后地勢起伏度的增速變緩,而這個(gè)變點(diǎn)所對(duì)應(yīng)的面積就是研究區(qū)提取地勢起伏度的最佳統(tǒng)計(jì)單元面積。 圖3 地勢起伏度與統(tǒng)計(jì)單元面積對(duì)應(yīng)關(guān)系擬合曲線Fig.3 Fitting curve of the corresponding relation between relief amplitude and area of statistical unit 對(duì)不同統(tǒng)計(jì)單元的平均地勢起伏度、不同統(tǒng)計(jì)單元地勢起伏度分級(jí)中的斜坡地質(zhì)災(zāi)害峰值、地勢起伏度峰值進(jìn)行統(tǒng)計(jì),結(jié)果見圖4. 圖4 平均地勢起伏度、地勢起伏度峰值、斜坡地質(zhì)災(zāi)害峰值與統(tǒng)計(jì)單元關(guān)系圖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 根據(jù)圖4統(tǒng)計(jì)結(jié)果分別構(gòu)建樣本序列,利用均值變點(diǎn)法的公式編制Excel程序,計(jì)算樣本序列的統(tǒng)計(jì)量S和Si,求得S-Si(見表3),構(gòu)建S與Si差值的變化擬合曲線。 表3 均值變點(diǎn)法分析統(tǒng)計(jì)結(jié)果Table 3 Statistical results of mean change point 平均地勢起伏度樣本序列S值為16.994 5,S與Si差值的變化曲線見圖5.地勢起伏度峰值樣本序列S的值為20.130 4,S與Si差值的變化曲線見圖6.斜坡地質(zhì)災(zāi)害峰值樣本序列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 基于斜坡地質(zhì)災(zāi)害峰值的S-Si變化曲線Fig.7 Variation curve based on peak value of slope geological hazards S-Si 1) 圖4顯示,每一種統(tǒng)計(jì)單元都存在一個(gè)斜坡地質(zhì)災(zāi)害發(fā)育的峰值區(qū)間,隨統(tǒng)計(jì)單元的增大,峰值逐漸減小,到達(dá)某一個(gè)臨界值后峰值衰減趨于穩(wěn)定。 2) 圖5顯示,i為11時(shí),S與Si差值達(dá)最大,i=11處即為變點(diǎn),對(duì)應(yīng)統(tǒng)計(jì)單元為12×12網(wǎng)格。因此,由平均地勢起伏度以均值變點(diǎn)法來分析得到研究區(qū)的最佳地勢起伏度提取單元為12×12網(wǎng)格,最佳統(tǒng)計(jì)面積為0.129 6 km2. 3) 圖6顯示,i為11時(shí),S與Si差值達(dá)最大,i=11處即為變點(diǎn),對(duì)應(yīng)統(tǒng)計(jì)單元為12×12網(wǎng)格。因此,由地勢起伏度峰值以均值變點(diǎn)法來分析得到研究區(qū)的最佳地勢起伏度提取單元為12×12網(wǎng)格,最佳統(tǒng)計(jì)面積為0.129 6 km2. 4) 圖7顯示,i為8時(shí),S與Si差值達(dá)最大,i=8處即為變點(diǎn),對(duì)應(yīng)統(tǒng)計(jì)單元為9×9網(wǎng)格。因此,通過斜坡地質(zhì)災(zāi)害峰值以均值變點(diǎn)法來分析得到研究區(qū)的最佳地勢起伏度提取單元為9×9網(wǎng)格,最佳統(tǒng)計(jì)面積為0.072 9 km2. 5) 對(duì)比三種樣本序列所計(jì)算得到的最佳地勢起伏度提取單元,使用平均地勢起伏度和地勢起伏度峰值樣本序列計(jì)算結(jié)果一致,均為12×12網(wǎng)格,而斜坡地質(zhì)災(zāi)害峰值樣本序列計(jì)算結(jié)果則為9×9網(wǎng)格。作為客觀存在的區(qū)域地貌特征因子,斜坡地質(zhì)災(zāi)害演變頻繁,而平均地勢起伏度演變較小。因此,以研究區(qū)數(shù)字地貌角度分析,基于平均地勢起伏度所計(jì)算的最佳提取單元更合理;而以地質(zhì)災(zāi)害發(fā)育角度分析,基于斜坡地質(zhì)災(zāi)害峰值所計(jì)算的最佳提取單元?jiǎng)t更合理。 6) 針對(duì)斜坡地質(zhì)災(zāi)害敏感性評(píng)價(jià)中地勢起伏度最佳統(tǒng)計(jì)單元開展研究,重點(diǎn)關(guān)注的是研究區(qū)地勢起伏度與斜坡地質(zhì)災(zāi)害發(fā)育之間的關(guān)系,因此所尋找的是能反映地勢起伏度與斜坡地質(zhì)災(zāi)害發(fā)育之間關(guān)系的最佳統(tǒng)計(jì)單元,而不是數(shù)字地貌研究中區(qū)域地勢起伏度計(jì)算的最佳統(tǒng)計(jì)單元。因此,選擇9×9網(wǎng)格作為研究區(qū)斜坡地質(zhì)災(zāi)害敏感性評(píng)價(jià)中地勢起伏度因子最佳提取單元。 1) 地勢起伏度與統(tǒng)計(jì)單元面積之間存在顯著的對(duì)數(shù)變化關(guān)系,研究區(qū)內(nèi)地勢起伏度隨統(tǒng)計(jì)單元面積的增大而迅速增大,當(dāng)統(tǒng)計(jì)單元面積達(dá)到一定數(shù)值后,增大速度逐漸變緩,最后近似趨于平穩(wěn),這為開展地勢起伏度最佳提取單元提供了基礎(chǔ)。 2) 地勢起伏度提取最佳統(tǒng)計(jì)單元與地質(zhì)災(zāi)害分布峰值、平均地勢起伏度和地勢起伏度峰值之間存在明顯的相關(guān)性,目前的研究主要基于平均地勢起伏度,利用均值變點(diǎn)法或主觀方式分析確定最佳提取單元,這給斜坡地質(zhì)災(zāi)害敏感性評(píng)價(jià)帶來一定的不確定性,所以應(yīng)綜合三種情況開展統(tǒng)計(jì)并綜合分析后確定研究區(qū)地勢起伏度最佳統(tǒng)計(jì)尺度。 3) 以研究區(qū)ASTER GDEM重采樣30 m分辨率DEM和2012年斜坡類地質(zhì)災(zāi)害空間分布信息為數(shù)據(jù)源,以2×2,3×3,…,25×25共24個(gè)分析窗口提取了研究區(qū)的地勢起伏度,以均值變點(diǎn)法對(duì)不同統(tǒng)計(jì)單元的平均地勢起伏度、不同統(tǒng)計(jì)單元地勢起伏度分級(jí)中的斜坡地質(zhì)災(zāi)害峰值、地勢起伏度峰值進(jìn)行統(tǒng)計(jì)分析。綜合三種樣本序列計(jì)算結(jié)果,本文選擇9×9網(wǎng)格作為研究區(qū)斜坡類地質(zhì)災(zāi)害敏感性評(píng)價(jià)時(shí)的最佳地勢起伏度提取單元,這為研究區(qū)斜坡類地質(zhì)災(zāi)害敏感性評(píng)價(jià)研究提供了可靠的地勢起伏度因子,同時(shí)也為同類研究提供了數(shù)字地形因子較為準(zhǔn)確的提取方法。2 結(jié)果與討論
2.1 平均地勢起伏度與統(tǒng)計(jì)單元之間的關(guān)系
2.2 地勢起伏度提取最佳統(tǒng)計(jì)單元分析
2.3 討論
3 結(jié)論