基于有效降雨強(qiáng)度和邏輯回歸的降雨型滑坡預(yù)測模型

盛逸凡，李遠(yuǎn)耀，徐 勇，吳吉明，林 巍

(1.中國地質(zhì)大學(xué)(武漢)地質(zhì)調(diào)查研究院，湖北 武漢 430074；2.中國地質(zhì)調(diào)查局武漢地質(zhì)調(diào)查中心， 湖北 武漢 430205)

滑坡是我國山地丘陵區(qū)主要的地質(zhì)災(zāi)害類型之一[1]，突發(fā)性強(qiáng)、分布廣泛，常常造成大量的人員傷亡和財(cái)產(chǎn)損失。大量調(diào)查研究發(fā)現(xiàn)，90%左右的滑坡發(fā)生在降雨期間或降雨之后，這種由降雨誘發(fā)的滑坡稱為降雨型滑坡[2]?；诮涤晔录c滑坡災(zāi)害間的關(guān)聯(lián)性，通過區(qū)域降雨特征來開展大范圍滑坡災(zāi)害的預(yù)測預(yù)報，是有效預(yù)防降雨型滑坡災(zāi)害的重要途徑。

圍繞降雨誘發(fā)滑坡的機(jī)理及預(yù)測預(yù)報方法，國內(nèi)外學(xué)者開展了大量研究。如：Thomas Glade等[3]采用統(tǒng)計(jì)方法研究得出降雨誘發(fā)、非降雨誘發(fā)這兩種情況下滑坡發(fā)生的降雨量閾值。F Guzzetti等[4]應(yīng)用貝葉斯統(tǒng)計(jì)方法從大量數(shù)據(jù)中分析得出閾值曲線，發(fā)現(xiàn)山區(qū)的有效降雨強(qiáng)度要低于以地中海氣候?yàn)樘卣鞯牡貐^(qū)。國內(nèi)方面，殷坤龍等[5]以浙江省為研究區(qū)，采用統(tǒng)計(jì)學(xué)方法分析了臺風(fēng)、非臺風(fēng)兩種情況下滑坡與降雨強(qiáng)度、降雨量的相關(guān)性，得到了誘發(fā)滑坡發(fā)生的降雨閾值曲線及臨界值。張珍等[6]、林孝松等[7]從降雨量、降雨周期、暴雨頻次、雨型、以及降雨持續(xù)時間等方面進(jìn)行了研究，指出了降雨與滑坡事件間的耦合。謝守益等[8]、李曉等[9]、謝劍明等[10]研究確定了區(qū)域性滑坡發(fā)生的臨界降雨量值及降雨誘發(fā)滑坡發(fā)生的概率等，為地方氣象部門滑坡預(yù)警提供了重要依據(jù)。

綜合而言，針對降雨型滑坡的區(qū)域預(yù)測預(yù)報，最關(guān)鍵是確定誘發(fā)滑坡的臨界降雨閾值。然而，在前人研究中，大多以統(tǒng)計(jì)方法和經(jīng)驗(yàn)?zāi)Ｐ蜑橹?，對于降雨型滑坡本身的?guī)模、形態(tài)等地形特征缺乏深入分析，常造成統(tǒng)計(jì)模型得出的宏觀臨界降雨閾值存在一定缺陷；此外，經(jīng)驗(yàn)性的定性閾值模型與定量的預(yù)測模型也結(jié)合較少，從而造成統(tǒng)計(jì)閾值模型的適用性不足。

本文以湖南桑植縣近30年的長時段歷史滑坡災(zāi)害為研究對象，分析了降雨因子與滑坡事件概率及滑坡數(shù)量間的相關(guān)性和偏相關(guān)性；系統(tǒng)考慮滑坡地形條件對降雨閾值結(jié)果的重要影響，對降雨臨界閾值統(tǒng)計(jì)結(jié)果按滑坡規(guī)模、滑體厚度、滑坡坡度進(jìn)行分類，得出了不同類型滑坡的分類分級臨界降雨閾值；并建立了基于邏輯回歸方法的滑坡定量預(yù)測預(yù)報模型。

1 研究區(qū)自然環(huán)境與滑坡災(zāi)害

1.1 研究區(qū)概況

桑植縣隸屬湖南省張家界市，地處武陵山脈北麓，鄂西山地南端，全縣總面積3 474 km2。區(qū)域?qū)賮啛釒Ъ撅L(fēng)氣候區(qū)，光熱充足，春暖夏熱多雨，秋涼冬冷雨少;年平均氣溫16.8 ℃.降雨量充沛，是湖南省的四大暴雨中心之一，多年降雨天數(shù)平均為140 d，多年平均降雨量1 417 mm，降雨量多集中于5—8月份，占全年降雨量的50%～60%。

區(qū)域地勢總趨勢為西北高、東南低，絕大部分地區(qū)地勢高峻，群山起伏，溝谷深切，構(gòu)造位于新華夏系一級構(gòu)造第三隆起帶的南段，湘西北弧形構(gòu)造的東北端，屬武陵山褶皺帶的東南端，主體為桑植復(fù)向斜。境內(nèi)地層除缺失白堊系、侏羅系、泥盆系下統(tǒng)及石炭系外，從第四系到寒武系地層其他均有出露，其中三疊系、志留系、寒武系分布范圍較廣，區(qū)內(nèi)地質(zhì)災(zāi)害主要發(fā)育在三疊系至志留系地層之中。

1.2 滑坡災(zāi)害發(fā)育特征

圖1 1987—2015年月平均降雨量與滑坡數(shù)量統(tǒng)計(jì)Fig.1 Monthly average rainfall and the number of landslides from 1987 to 2015

根據(jù)湖南省地質(zhì)調(diào)查院提供的桑植縣地質(zhì)災(zāi)害調(diào)查與區(qū)劃報告得知，桑植縣為地質(zhì)災(zāi)害高發(fā)區(qū)，地質(zhì)災(zāi)害類型包括滑坡、崩塌、泥石流和地面塌陷、地面沉降、地裂縫等，其中滑坡是最主要地質(zhì)災(zāi)害，占災(zāi)害總數(shù)的94.1%?；骂愋鸵詼\層土質(zhì)滑坡為主，規(guī)模為中小型，一般發(fā)生在暴雨期間。根據(jù)張家界市氣象局和桑植縣雨量站觀測數(shù)據(jù)，結(jié)合中國地質(zhì)調(diào)查局武漢地調(diào)中心提供的滑坡資料，統(tǒng)計(jì)分析了研究區(qū)內(nèi)近30年的月平均降雨量與滑坡發(fā)生數(shù)量(圖1)，可發(fā)現(xiàn)滑坡事件頻數(shù)與降雨量的大小密切相關(guān)。現(xiàn)場調(diào)查表明，地形地貌、地層巖性等地質(zhì)條件是控制區(qū)內(nèi)滑坡形成和發(fā)展的內(nèi)因，降雨是滑坡形成的主要外部條件之一?；驴臻g分布受地形地貌控制，分布于河溪兩岸，發(fā)育于向斜、背斜及斷裂附近，65%的滑坡的坡度分布在15°～30°之間。

2 滑坡災(zāi)害閾值判別模型

2.1 降雨量與滑坡的關(guān)系

在張家界市氣象局和桑植縣雨量站提供的近30年降雨數(shù)據(jù)和滑坡災(zāi)害資料的基礎(chǔ)上，為了排除非降雨因素的影響，本文遴選出137例降雨誘發(fā)滑坡事件，表1列出了5例典型的降雨誘發(fā)滑坡事件和5例未發(fā)生滑坡的降雨事件。表中P0為當(dāng)日降雨量，P1為滑坡前第1日降雨量，P2為滑坡前第2日降雨量，……，依此類推。由于桑植縣降雨誘發(fā)滑坡的統(tǒng)計(jì)數(shù)據(jù)只精確到日，本文假設(shè)滑坡是在某日0-24時內(nèi)發(fā)生，因此降雨量數(shù)據(jù)也取同時段內(nèi)的數(shù)據(jù)。

表1 桑植縣部分降雨誘發(fā)滑坡數(shù)據(jù)Table 1 Partial data of rainfall-induced landslides in Sangzhi district

表1表明前期降雨對滑坡發(fā)生起主要作用，由于前期降雨的持續(xù)作用，土體逐漸處于飽水狀態(tài)，巖土體抗剪強(qiáng)度下降，即使在當(dāng)日雨量很小的情況下，一旦有效降雨量高于臨界閾值線，就可能誘發(fā)滑坡。所以，對于當(dāng)日降雨量小或暴雨類型的滑坡，前期有效降雨量具有顯著作用。

2.2 有效降雨量模型

有效降雨量是指對誘發(fā)滑坡發(fā)生作用的那部分降雨，這一概念最早是由Crozier在1980年引入地質(zhì)災(zāi)害降雨閾值研究中的，近年來被廣泛用于降雨誘發(fā)滑坡的研究中。降雨導(dǎo)致滑坡發(fā)生的過程可分兩個階段：第一階段是前期降雨期，是指滑坡發(fā)生前一定時間內(nèi)的降雨，降雨通常會使坡體土顆粒之間的空隙增大，隨之孔隙水壓力增加，有效應(yīng)力降低，剪切強(qiáng)度降低，繼而誘發(fā)滑坡災(zāi)害的發(fā)生[13-14]；第二階段為持續(xù)降雨期，這部分降雨在一定程度上直接激發(fā)了滑坡的發(fā)生。由于每次降雨事件只有部分降雨量對滑坡的發(fā)生起作用，所以本文引入了有效降雨量，有效降雨量的計(jì)算公式為：

Pa=R0+KR1+K2R2+K3R3+…+KnRn

(1)

式中：Pa——有效降雨量；

K——降雨衰減系數(shù)；

R0——滑坡發(fā)生當(dāng)日降雨量；

R1——次日降雨量，依此類推，n為滑坡發(fā)生前的天數(shù)。

自Crozier以來，計(jì)算前期有效降雨量時大多采用式(1)，通常需要根據(jù)相關(guān)地區(qū)災(zāi)害特點(diǎn)和地質(zhì)條件的不同合理選取K值，或者對式(1)中降雨衰減系數(shù)的形式稍作改變。

桑植縣由于山區(qū)面積較大，河流沖刷強(qiáng)，地下水動態(tài)變化復(fù)雜，雨季主要集中在5—8月份，特別是7月氣溫較高，降雨量充沛，蒸發(fā)作用明顯，具有較為明顯的地域特征，合理地選擇前期降雨衰減系數(shù)可使滑坡預(yù)測更精確。因此本文用SPSS22.0軟件做降雨因子與滑坡是否發(fā)生及滑坡數(shù)量的相關(guān)分析和偏相關(guān)分析[15-16]。

2.3 有效降雨衰減系數(shù)

在構(gòu)建降雨型滑坡預(yù)測模型時，合理選取降雨因子是首要的問題。本文中的降雨因子考慮了當(dāng)日降雨量R0、滑坡發(fā)生前第1日的降雨量R1、……、第6日的降雨量R6；為進(jìn)一步分析滑坡是否發(fā)生及滑坡發(fā)生數(shù)量與降雨因子的相關(guān)關(guān)系，引入累積降雨因子和前期降雨因子兩個概念[11]。其中，累積降雨因子是指滑坡發(fā)生當(dāng)日及前一段時間內(nèi)的降雨總量，前期降雨因子是指累積降雨量減去滑坡當(dāng)日的降雨量。利用SPSS22.0軟件分析了降雨因子與滑坡數(shù)量及滑坡是否發(fā)生間的相關(guān)性和偏相關(guān)性[1]，結(jié)果見表2、表3。

偏相關(guān)分析是指不考慮其它變量的影響下，對兩個所關(guān)心的變量進(jìn)行相關(guān)性分析，所得結(jié)果能夠更好地反映真實(shí)情況。從表2可以看出，滑坡是否發(fā)生及滑坡發(fā)生數(shù)量與當(dāng)日降雨量呈顯著相關(guān)和顯著偏相關(guān)(在0.01顯著性水平下)，這說明誘發(fā)滑坡最為顯著的因子是當(dāng)日降雨量R0。由表2可知，從R0到R6，相關(guān)系數(shù)及偏相關(guān)系數(shù)呈現(xiàn)逐漸下降的趨勢；從R4開始，降雨因子對滑坡發(fā)生的影響已變得很小，可忽略不計(jì)。同時發(fā)現(xiàn)，當(dāng)日降雨量及前一天降雨量與滑坡發(fā)生數(shù)量之間呈顯著偏相關(guān)，這說明大部分滑坡發(fā)生在暴雨當(dāng)日或暴雨的后一天。

表2 降雨因子與滑坡相關(guān)性及偏相關(guān)性Table 2 Correlation and partial correlation between rainfall factors and landslides

注：** 表示在0.01顯著性水平下，相關(guān)性及偏相關(guān)性顯著。

表3 累積降雨因子與滑坡的相關(guān)性Table 3 Correlation between effective cumulative rainfall factors and landslides

注：** 表示在0.01顯著性水平下，相關(guān)性顯著。

由表3可知，滑坡是否發(fā)生與累積降雨因子呈顯著相關(guān)，表現(xiàn)出先升后降的趨勢，在累積4天達(dá)到最大值；前期5天降雨與滑坡發(fā)生數(shù)量及滑坡是否發(fā)生的相關(guān)系數(shù)很小，這說明前期降雨對誘發(fā)滑坡的作用有限，滑坡當(dāng)日的降雨R0不能忽視。綜上，選用累積4天降雨量(即滑坡前3日及當(dāng)日降雨量)來建立滑坡的預(yù)測預(yù)報模型較為合理，即式(1)中的n確定為3。

式(1)中前期降雨衰減系數(shù)K的選擇會對滑坡預(yù)報的準(zhǔn)確性產(chǎn)生重要影響，因此前期降雨衰減系數(shù)的合理選擇尤為關(guān)鍵。分別設(shè)置不同的前期有效降雨衰減系數(shù)K值(從1到0.5依次減小)，并計(jì)算出4 d有效累積降雨量，然后分別分析不同衰減系數(shù)下，滑坡發(fā)生數(shù)量及滑坡是否發(fā)生與累積4 d降雨量的相關(guān)性。從表4可以發(fā)現(xiàn)前期有效降雨衰系數(shù)K= 0.8時與滑坡是否發(fā)生及滑坡發(fā)生數(shù)量的相關(guān)性最高，綜合分析統(tǒng)計(jì)結(jié)果，本地區(qū)前期有效降雨衰減系數(shù)K設(shè)置為0.8最為合適。

表4 不同衰減系數(shù)下累積4天降雨量與滑坡間的相關(guān)性Table 4 Correlation of 4 days cumulative rainfall and landslide under different attenuation coefficients

注：**表示在0.01顯著性水平下，相關(guān)性顯著。

2.4 計(jì)算降雨閾值的具體步驟

本文中降雨閾值指標(biāo)為有效降雨強(qiáng)度(有效累計(jì)降雨量與降雨歷時的比值)和降雨歷時，計(jì)算降雨閾值的步驟如下：

(1)統(tǒng)計(jì)桑植縣已有滑坡，結(jié)合滑坡災(zāi)害的規(guī)模、坡度和厚度，按體積將研究區(qū)滑坡分為大、中、小規(guī)模三類，按坡度將滑坡分為緩、中、陡坡度三類，按厚度將滑坡分為較小、中等、較大厚度三類。規(guī)模方面，將小于30 000 m3的滑坡定為小規(guī)模，30 000～100 000 m3之間的定為中規(guī)模，大于100 000 m3的定為大規(guī)模；坡度方面，將0～30°的滑坡定為緩坡度，30°～45°的定為中等坡度，45°～82°定為陡坡度；厚度方面，將0～3 m的滑坡定為厚度較小，3～5 m定為中等厚度，5～15 m定為厚度較大。

(2)以降雨持續(xù)時間的對數(shù)為橫軸，以有效降雨強(qiáng)度的對數(shù)為縱軸，按規(guī)模、坡度、厚度分別統(tǒng)計(jì)滑坡有效降雨量(I)和降雨歷時(D)的關(guān)系，得到不同規(guī)模、坡度、厚度下有效降雨強(qiáng)度和持續(xù)時間的滑坡事件的I-D圖(圖3～5)。

(3)根據(jù)I-D圖，分別求取10%、50%、90%概率下誘發(fā)滑坡的有效降雨強(qiáng)度閾值線。其方程形式為I=aD-b，最后采用頻數(shù)法進(jìn)行降雨閾值的擬合。不同條件下誘發(fā)滑坡發(fā)生的有效降雨強(qiáng)度閾值方程見表5。

圖2 滑坡規(guī)模劃分曲線及小、中、大規(guī)模體積滑坡I-D圖Fig.2 Intensity-duration(I-D) thresholds of landslides with large, medium and small volume

圖3 滑坡坡度劃分曲線及緩、中、陡坡度滑坡I-D圖Fig.3 Intensity-duration (I-D) thresholds of landslides with moderate, medium and steep slopes

圖4 滑坡厚度劃分曲線及小中大厚度滑坡I-D圖Fig.4 Intensity-duration (I-D) of landslides with large, medium and small thickness

(4)根據(jù)有效降雨強(qiáng)度閾值對滑坡災(zāi)害的危險性進(jìn)行分級，可分為高、中、低3個等級[1，17]。當(dāng)預(yù)警點(diǎn)落在90%閾值線的上方區(qū)域時，則說明有效降雨強(qiáng)度已經(jīng)超過滑坡發(fā)生概率為90%時的臨界值，此預(yù)警區(qū)為高危險區(qū)；當(dāng)預(yù)警點(diǎn)落在50%閾值線的上方時，則說明有效降雨強(qiáng)度超過了滑坡發(fā)生概率為50%時的降雨強(qiáng)度，此預(yù)警區(qū)為中危險區(qū)；10%到50%之間的為低危險區(qū)；10%閾值線以下的區(qū)域即無危險區(qū)(表6)。

表5 不同規(guī)模、坡度和厚度滑坡的有效降雨強(qiáng)度閾值線Table 5 Effective rainfall intensity threshold line for landslides of different sizes, slopes and thicknesses

表6 滑坡危險性分級Table 6 Risk grades distribution for landslide

3 基于統(tǒng)計(jì)學(xué)方法的定量預(yù)測模型

3.1 基于邏輯回歸的滑坡預(yù)測模型

邏輯回歸是一種基于統(tǒng)計(jì)學(xué)和確定性模型相結(jié)合的分析方法[18]?？紤]到降雨因子為非連續(xù)變量，本文使用邏輯回歸分析中的二元邏輯回歸方法[1]。

假設(shè)某研究區(qū)發(fā)生滑坡的概率為p，p的取值范圍為0到1，則不發(fā)生滑坡的概率為1-p，p/(1-p)為滑坡邏輯回歸值。對其取自然對數(shù)，自變量為X1，X2，…，Xk，其中X1為當(dāng)日降雨，X2為前一天降雨，……，依次類推；因變量為滑坡發(fā)生概率p，則邏輯回歸方程可表示為[20]：

(2)

式中：β1，β2，…，βk——回歸系數(shù)。

變換式(2)可得：

(3)

利用式(3)可以定量預(yù)測降雨型滑坡災(zāi)害發(fā)生的概率。

3.2 滑坡預(yù)測模型的建立

將滑坡發(fā)生定義為邏輯值1，未發(fā)生定義為邏輯值0，把概率p= 0.5作為預(yù)測降雨型滑坡是否發(fā)生的臨界值點(diǎn)，p大于0.5時被認(rèn)為會發(fā)生滑坡，p小于0.5時則認(rèn)為不會發(fā)生滑坡。

本文首先從桑植縣137例樣本中遴選出120例發(fā)生滑坡事件。根據(jù)前文的分析，發(fā)現(xiàn)選用累積4天(即滑坡當(dāng)日及前3日)降雨量、有效降雨衰減系數(shù)k=0.8時，建立滑坡的預(yù)測預(yù)報模型最為合理的結(jié)論。這兩個關(guān)鍵因子的確定為基于二元邏輯回歸的滑坡定量預(yù)測模型的準(zhǔn)確度提供了基本保證。

將當(dāng)日降雨量(R0)和前3日有效累計(jì)降雨量(Ra3)作為自變量，將滑坡是否發(fā)生(0或1，即二元邏輯)作為因變量，根據(jù)選取的120組樣本值建立回歸方程，另選10組降雨事件作為檢驗(yàn)樣本。經(jīng)回歸分析，得到滑坡發(fā)生的概率預(yù)測模型[13]為：

(4)

為了驗(yàn)證式(4)的準(zhǔn)確性和實(shí)用性，選取桑植縣誘發(fā)滑坡的5次降雨事件及未發(fā)生滑坡的5次降雨事件進(jìn)行驗(yàn)證。將降雨因子代入式(4)，結(jié)果見表7，概率預(yù)測值與實(shí)際情況較為吻合。以概率值0.5作為滑坡是否發(fā)生的分界值，結(jié)果表明本文建立的滑坡預(yù)測模型準(zhǔn)確性與實(shí)際情況吻合，見表8。

表7 10組樣本預(yù)測值Table 7 Predicted of 10 sets of samples

表8 模型預(yù)測分類表Table 8 Classification of model prediction

3.3 臨界閾值的確定

通常所說的臨界值是指當(dāng)研究區(qū)的降雨量到達(dá)某一個值時滑坡會發(fā)生，它實(shí)際上就是一個界限閾值，具體指滑坡事件發(fā)生在量上所需要的一個最大值或者一個最小值。因此，本文定義降雨量臨界值為最小降雨量值。當(dāng)降雨量大于或等于這個界限值時，會發(fā)生滑坡。桑植縣符合要求的降雨型滑坡事件共120例，因此選用這120個滑坡發(fā)生當(dāng)日降雨量和前三日有效累積降雨量這兩個評價因子，繼而得到臨界值[20]表達(dá)式為：

R0=-0.566Ra3+70.61

(5)

式(5)中R0為當(dāng)日降雨量，Ra3為前三日有效累積降雨量。當(dāng)Ra3= 0時，即無前期降雨量，誘發(fā)滑坡的日降雨量臨界值為70.61 mm；當(dāng)R0= 0時，即無當(dāng)日降雨，誘發(fā)滑坡前期降雨值為124.75 mm。

4 結(jié)論

(1)以湖南省張家界市桑植縣為研究區(qū)，分別分析了滑坡發(fā)生數(shù)量及滑坡是否發(fā)生與降雨因子的相關(guān)性和偏相關(guān)性，結(jié)果表明，距滑坡發(fā)生后天數(shù)的增加，降雨因子與滑坡的相關(guān)性呈現(xiàn)出下降的趨勢；當(dāng)前期有效降雨衰減系數(shù)K=0.8 時，4天有效累積降雨量(即滑坡當(dāng)日及前3日降雨量)與滑坡是否發(fā)生的相關(guān)系數(shù)最高。

(2)按滑坡的規(guī)模、坡度和厚度大小統(tǒng)計(jì)了桑植縣滑坡災(zāi)害點(diǎn)的有效降雨強(qiáng)度和降雨持續(xù)時間數(shù)據(jù)，分別給出了不同類型滑坡的I-D圖，定性得到了不同概率下誘發(fā)滑坡發(fā)生的有效降雨強(qiáng)度經(jīng)驗(yàn)閾值，并對滑坡災(zāi)害進(jìn)行了危險性等級劃分。

(3)從定量角度建立了基于當(dāng)日降雨量和前三日有效累積降雨量為自變量的二元邏輯回歸方法的滑坡預(yù)報模型，經(jīng)驗(yàn)證表明其預(yù)測結(jié)果較為準(zhǔn)確。在此基礎(chǔ)上，給出了臨界值的表達(dá)式，為桑植縣滑坡預(yù)測預(yù)報提供了依據(jù)。