韋 成，林毓華，韋 寧，江 杰，蔡杏珍，陳銘熙

（1.廣西工業(yè)設(shè)計集團(tuán)有限公司，廣西南寧 530020； 2.廣西大學(xué)土木建筑工程學(xué)院，廣西南寧 530004）

0 引言

滑坡是我國主要的地質(zhì)災(zāi)害形式，每年造成巨大的生命財產(chǎn)損失，因此，如何準(zhǔn)確地預(yù)測滑坡時間成為滑坡災(zāi)害防治的重要任務(wù)之一［1-3］。自1961年SAITO 和UEZAWA 提出第一個實(shí)用的滑坡時間預(yù)測模型以來，國內(nèi)外學(xué)者已經(jīng)提出了大量經(jīng)驗(yàn)或半經(jīng)驗(yàn)預(yù)測模型用于滑坡時間預(yù)測［4］。FUKUZONO基于室內(nèi)降雨滑坡模型試驗(yàn)研究，發(fā)現(xiàn)滑坡在加速變形階段的加速度與速度存在普適性經(jīng)驗(yàn)關(guān)系：v=Avα，其中A 和α 為常數(shù)，并對該式積分提出了著名的逆速度法（INV）［5］。VOIGHT 對FUKUZONO 的理論進(jìn)行數(shù)學(xué)推廣，建立了表征材料破壞出現(xiàn)自加速行為的廣義定律［6］。李天斌和陳明東認(rèn)為典型的滑坡位移-時間曲線與Verhulst 函數(shù)的反函數(shù)曲線形態(tài)類似，提出采用Verhulst 反函數(shù)模型對滑坡時間進(jìn)行預(yù)測并驗(yàn)證了模型的有效性［7］。MUFUNDIRWA 等在對FUKUI 預(yù)測模型研究的基礎(chǔ)上，提出了可以不考慮滑坡破壞機(jī)制、巖性和任何其他參數(shù)的斜率法（SLO）［8-9］。雖然提出的預(yù)測模型眾多，但由于受滑坡演化非線性特征、監(jiān)測誤差、環(huán)境噪音以及預(yù)測模型固有局限性等多種因素的影響，準(zhǔn)確預(yù)測滑坡的發(fā)生時間仍是世界難題［10-12］。近年來，一些學(xué)者逐漸關(guān)注到滑坡預(yù)測存在的不確定性問題，采用區(qū)間預(yù)測、概率預(yù)測等方法對滑坡時間進(jìn)行不確定性估計［13-14］。例如，采用移動平均法對監(jiān)測數(shù)據(jù)進(jìn)行處理，提出了“破壞時間窗口”的區(qū)間預(yù)測方法對滑坡時間可能發(fā)生的時間段進(jìn)行估計［13］；ZHANG 等同時考慮觀測不確性和模型不確定性對預(yù)測的影響，采用逆速度法結(jié)合極大似然估計對滑坡時間進(jìn)行概率預(yù)測，獲取滑坡最可能發(fā)生時刻和時間段［14］。貝葉斯理論作為一種概率反分析的重要手段，廣泛應(yīng)用于巖土體的參數(shù)估計、變形預(yù)測及災(zāi)害預(yù)警中［15-17］，但應(yīng)用于滑坡時間概率預(yù)測的研究較少。ZHANG 等采用貝葉斯機(jī)器學(xué)習(xí)方法對滑坡時間進(jìn)行概率預(yù)測，但未量化實(shí)時預(yù)測的性能。當(dāng)前，在量化滑坡預(yù)測不確定性方面取得了一定的進(jìn)展，但是考慮滑坡實(shí)時預(yù)測過程不確定性方面的研究依然有限，更多的是利用最后一次觀測獲取的數(shù)據(jù)對滑坡時間進(jìn)行事后分析［18］。另外，采用貝葉斯理論進(jìn)行滑坡時間預(yù)測的研究正處于起步階段。

本文的目的是提出一種基于貝葉斯更新并結(jié)合多種預(yù)測模型的滑坡時間概率預(yù)測方法，對實(shí)時預(yù)測存在的不確定性進(jìn)行量化并進(jìn)行預(yù)測決策。首先，選用逆速度法、FUKUI 模型和斜率法三種模型用于滑坡時間預(yù)測。其次，基于貝葉斯理論建立了三種預(yù)測模型的不確定性預(yù)測框架，并采用兩個滑坡案例驗(yàn)證提出方法的可行性。最后，提出多模型預(yù)測決策的建議。

1 滑坡時間預(yù)測模型

1.1 逆速度法(INV)

逆速度法是由福囿輝旗基于室內(nèi)降雨滑坡模型試驗(yàn)研究的基礎(chǔ)上提出的，主要用于滑坡加速變形階段的滑坡時間預(yù)測［5］。該假定滑坡破壞時速度無限，則速度倒數(shù)1/v趨于0 對應(yīng)的時刻即為滑坡時間tf，預(yù)測模型表達(dá)式：

式中：t為觀測時間。模型預(yù)測曲線1/v-t曲線見圖1a。研究表明，α值大多為1.5～2.2，預(yù)測曲線近似為直線［5］，因此實(shí)際預(yù)測中常假定α=2，模型變換為直線型逆速度法（LINV）：

圖1 （a）逆速度法和（b）斜率法預(yù)測曲線［5，8］Figure 1 （a）Reverse velocity method and（b）slope method for predicting curves［5，8］

1.2 FUKUI預(yù)測模型(FM)

FUKUI 等在研究多胡砂巖、水泥砂漿、秋芳大理石等材料單軸壓縮試驗(yàn)的蠕變特性時發(fā)現(xiàn)這些材料在臨近破壞時應(yīng)變ε 與剩余壽命(tf-t) 的對數(shù)成線性關(guān)系［9］：

式中：a和b為常數(shù)。模型在滑坡時間預(yù)測的有效性已得到證明［14］，滑坡預(yù)測中ε替換成滑坡位移s。同樣，F(xiàn)UKUI 預(yù)測模型更適用于加速變形階段的滑坡時間預(yù)測。

1.3 斜率法(SLO)

MUFUNDIRWA 等對FUKUI 預(yù)測模型（式（3））進(jìn)行求導(dǎo)并變換可得［8］：

該模型預(yù)測曲線是一條直線（圖1b），其中速度v與時間t的乘積（vt）為因變量，v為自變量，滑坡時間為斜率，因此該方法被稱為斜率法。由于該模型由Fukui 預(yù)測模型得出，同樣適用于加速變形階段的滑坡時間預(yù)測。

2 貝葉斯更新預(yù)測框架與集成預(yù)測

2.1 貝葉斯更新方法

貝葉斯理論用于滑坡預(yù)測的基本思想為：根據(jù)先驗(yàn)知識對預(yù)測模型的未知參數(shù)假定一個先驗(yàn)分布，而后基于新加入監(jiān)測數(shù)據(jù)持續(xù)地對模型參數(shù)進(jìn)行更新，更新后參數(shù)的分布為后驗(yàn)分布［15］。假設(shè)M（θ）為預(yù)測模型，θ為模型參數(shù)向量，D=［d1，d2，…，dk］為監(jiān)測數(shù)據(jù)的向量。預(yù)測模型和監(jiān)測數(shù)據(jù)二者關(guān)系如下：

式中：εr為預(yù)測誤差?；谪惾~斯理論，隨機(jī)變量θ的后驗(yàn)分布p（θ|D）可更新為

式中：p（D|θ）為似然函數(shù)；p（θ）為先驗(yàn)分布；p（D）為歸一化因子。p（D）可表示為

假定預(yù)測誤差相互獨(dú)立且服從正態(tài)分布，則似然函數(shù)可寫成

式中：σε為預(yù)測誤差的標(biāo)準(zhǔn)差。

2.2 預(yù)測模型先驗(yàn)分布確定與后驗(yàn)分布求解方法

模型LINV 的參數(shù)向量由A、tf和σε三個參數(shù)組成，即θ=｛A，tf，σε｝。在加速變形階段，模型中的A大于0，且據(jù)統(tǒng)計顯示A值大小大多不超過80［10，19］。A的具體分布未知，可以采用均勻分布作為先驗(yàn)分布，均勻分布可取大于80 的數(shù)值A(chǔ)∞，則A～U（0，A∞）。tf的先驗(yàn)分布同樣可以采用均勻分布表示。滑坡進(jìn)入加速變形階段，歷時一般不會太長，大多數(shù)情況小于300 d，因此tf的先驗(yàn)分布可表示為tf～U（tp，tf∞），其中tp，為當(dāng)前觀測時間，tf∞為大于300 d的數(shù)值。σε采用正態(tài)分布作為先驗(yàn)分布，均值和標(biāo)準(zhǔn)差分別取0 和1，即σε～N（0，1）。對于模型FM 和SLO，tf和σε的先驗(yàn)分布與LINV 的一致。模型FM的a和b在加速變形階段均為大于0，二者分布未知，因此它們的先驗(yàn)分布同樣假定為均勻分布，即a～U（0，a∞）和b～U（0，b∞），其中a∞和b∞為相對較大的數(shù)值，與監(jiān)測位移的量級相關(guān)，a∞和b∞一般可取位移值的2～3 倍。后驗(yàn)分布采用No-U-Turn Sampler（NUTS）算法進(jìn)行求解。NUTS 是由HOFFMAN 和GELMAN 提出的一種改進(jìn)漢密爾頓蒙特卡羅算法。NUTS 引入“Build binary tree”技術(shù)，可動態(tài)自適應(yīng)設(shè)置迭代次數(shù)和步長，而不需要人工干預(yù)，大大提升了收斂速度［20］。NUTS 可通過Python 語言的軟件庫Pymc3進(jìn)行調(diào)用。

2.3 集成預(yù)測與預(yù)測決策步驟

滑坡時間預(yù)測步驟如下：

1）步驟1。預(yù)測模型確定和初始監(jiān)測數(shù)據(jù)準(zhǔn)備。預(yù)測模型采用LINV、SLO 和FM。監(jiān)測數(shù)據(jù)采用三個案例加速變形階段的監(jiān)測數(shù)據(jù)，初始預(yù)測序列為加速階段前4 個位移數(shù)據(jù)點(diǎn)或前三個速度數(shù)據(jù)點(diǎn)。

2）步驟2。假定模型參數(shù)為不確定變量，并確定模型參數(shù)的先驗(yàn)分布。

3）步驟3。采用式（8）定義各預(yù)測模型的似然函數(shù)。

4）步驟4。使用NUTS 算法對各模型參數(shù)的后驗(yàn)分布進(jìn)行采樣。NUTS 進(jìn)行單鏈MCMC 仿真，樣本數(shù)為2 000。依據(jù)每次加入的最新監(jiān)測數(shù)據(jù)獲取新的后驗(yàn)分布，并保留每次獲取后驗(yàn)分布的樣本。預(yù)測值均值與真實(shí)值的誤差采用殘差絕對值E表示：

5）步驟5。三個預(yù)測模型每次更新后驗(yàn)分布的樣本疊加，形成新的后驗(yàn)分布。

6）步驟6。對預(yù)測結(jié)果進(jìn)行決策。

3 實(shí)例分析

3.1 天荒坪抽水蓄能電站開關(guān)站滑坡

該滑坡地處浙江省湖州市安吉縣，水電站下水庫左岸進(jìn)水口之上、上下庫連接公路之下［21］。坡體由花崗斑巖和含礫流紋質(zhì)熔凝灰?guī)r組成。1996年3月10 日，受開挖、爆破、降雨等因素的作用，坡體發(fā)生整體失穩(wěn)滑動，體積約為6.2 萬m3。這里選取裂縫計J4 的位移時間序列進(jìn)行分析（圖2）。本案例以1996 年2 月1 作為加速起點(diǎn)，并令加速起點(diǎn)時間t0=1 d，則真實(shí)滑坡時間Tf=28 d。

圖2 天荒坪抽水蓄能電站開關(guān)站滑坡的累積位移時間序列Figure 2 Accumulated displacement time series of landslide in Tianhuangping pumped storage power station

表1 和圖3 為模型LINV、SLO、FM 以及三模型集成最后一次觀測的預(yù)測結(jié)果和概率分布圖。在單一模型預(yù)測中，模型LINV 和FM 的采樣區(qū)間和95%置信區(qū)間（CI）將真實(shí)滑坡時間Tf包含其中，而SLO預(yù)測滑坡時間的采樣區(qū)間和95%置信區(qū)間均早于Tf。LINV 的預(yù)測均值（tf= 29.26 d）晚于Tf，F(xiàn)M（tf= 27.38 d）和SLO（tf=26.27d）的預(yù)測均值略早于Tf。LINV 的預(yù)測標(biāo)準(zhǔn)差（5.78 d）遠(yuǎn)大于另外兩個模型（0.39 d），表明LINV 受模型本身的誤差和監(jiān)測數(shù)噪音的影響較大，預(yù)測具有明顯不確定性。各模型預(yù)測最可能的滑坡時間（即最大后驗(yàn)估計值）均早于Tf，但總體上與預(yù)測均值的差異不大。三個預(yù)測模型集成的預(yù)測結(jié)果顯示其整合的采樣區(qū)間和95%置信區(qū)間均將Tf包含其中。最可能發(fā)生滑坡的預(yù)測時間為26.85 d，早于Tf。集成模型的預(yù)測均值（tf= 27.64 d）比三個單一模型的預(yù)測均值更接近Tf。

表1 最后一次觀測的滑坡時間預(yù)測結(jié)果Table 1 Prediction results of the last observation landslide time

圖3 最后一次觀測滑坡時間預(yù)測結(jié)果：（a）LINV；（b）FM；（c）SLO；（d）模型集成Figure 3 Prediction results of the last observation landslide time：（a）LINV；（b）FM；（c）SLO；（d）model integration

由圖4、圖5 可看出，無論是單一預(yù)測模型還是集成預(yù)測，在初始加速階段，95%置信區(qū)間相對較寬，表明預(yù)測具有明顯的不確定性，滑坡的破壞趨勢仍很模糊，不可能獲得可靠的預(yù)測；在臨滑階段，置信區(qū)間變窄或?qū)挾茸兓呌诜€(wěn)定，表明滑坡破壞趨勢變得清晰，可獲得相對可靠的預(yù)測結(jié)果。與另外兩個單一預(yù)測模型相比，LINV預(yù)測的置信區(qū)寬度較大，代表預(yù)測不確定性非常強(qiáng)，但在預(yù)測的時間尺度內(nèi)均將Tf容納其中。三個模型集成的預(yù)測結(jié)果，降低了LINV 預(yù)測的不確定性，并且在t=14 d之后將Tf包含其中。各模型預(yù)測的殘差絕對值也基本隨時間的增長而呈減小趨勢，其中集合預(yù)測的殘差絕對值的大小介于其它三個模型殘差絕對值之間，且演化趨勢較為平穩(wěn)。

圖4 實(shí)時預(yù)測結(jié)果：（a）LINV；（b）FM；（c）SLO；（d）模型集成Figure 4 Real-time prediction results：（a）LINV；（b）FM；（c）SLO；（d）model integration

圖5 開關(guān)站滑坡的滑坡時間預(yù)測值與觀測值的殘差絕對值Figure 5 Absolute residual value between predicted and observed landslide time of switch station

單一模型和模型集成的預(yù)測在前期加速變形階段同樣表現(xiàn)出明顯的不確定性，但在t= 80.01 d之后預(yù)測的置信區(qū)間變得非常窄，并且Tf基本位于模型集成預(yù)測的置信區(qū)間之內(nèi)，代表破壞趨勢已經(jīng)非常明顯，可對滑坡時間進(jìn)行可靠的預(yù)測。

3.2 準(zhǔn)格爾旗某露天煤礦滑坡

該滑坡案例來源于內(nèi)蒙古自治區(qū)鄂爾多斯市準(zhǔn)格爾旗的某露天煤礦［22］，于2020 年1 月4 日早上發(fā)生滑坡事件?；伦冃纹茐牡睦鄯e位移時間序列見圖6。根據(jù)原文獻(xiàn)提供的數(shù)據(jù)，初始位移的相對時間為1.08 h，滑坡從t=60.09 h 開始加速，真實(shí)滑坡時間Tf為94.23 h。

圖6 準(zhǔn)格爾旗某露天煤礦滑坡的累積位移時間序列Figure 6 Accumulated displacement time series of landslide in an opencast coal mine in Zhungeer

表2和圖7為模型LINV、SLO、FM以及三模型集成最后一次觀測的預(yù)測結(jié)果和概率分布圖。在單一模型預(yù)測中，Tf位于模型LINV 和SLO 的采樣區(qū)間之內(nèi)，但只有SLO 的95%置信區(qū)間能包含Tf。LINV 預(yù)測的置信區(qū)間晚于Tf，而FM 預(yù)測的置信區(qū)早于Tf。然而，三個模型最后一次觀測的預(yù)測均值均非常接近于Tf，誤差在2 h以內(nèi)。與上一案例相同，LINV 的預(yù)測標(biāo)準(zhǔn)差是三個模型最大的。最大后驗(yàn)估計值與Tf接近。模型集成的預(yù)測置信區(qū)間將容納其中，最大后驗(yàn)估計值略早于Tf。集成的預(yù)測結(jié)果在一定程度上降低了單個預(yù)測模型帶來的不可靠性。

表2 最后一次觀測的滑坡時間預(yù)測結(jié)果Table 2 Prediction results of the last observation landslide time

圖7 最后一次觀測滑坡時間預(yù)測結(jié)果：（a）LINV；（b）FM；（c）SLO；（d）模型集成Figure 7 Prediction results of the last observation landslide time：（a）LINV；（b）FM；（c）SLO；（d）model integration

圖8為實(shí)時預(yù)測結(jié)果，圖9為預(yù)測均值與真實(shí)滑坡時間的殘差值。結(jié)果顯示，單一模型和模型集成的預(yù)測在前期加速變形階段同樣表現(xiàn)出明顯的不確定性，但在t= 80.01 d 之后預(yù)測的置信區(qū)間變得非常窄，并且Tf基本位于模型集成預(yù)測的置信區(qū)間之內(nèi)，代表破壞趨勢已經(jīng)非常明顯，可對滑坡時間進(jìn)行可靠的預(yù)測。與前一個案例類似，該案例預(yù)測的殘差絕對值同樣表現(xiàn)為E隨時間的增長而呈減小趨勢，其中集合預(yù)測的殘差絕對值的大小介于其它三個模型殘差絕對值之間，且演化趨勢較為平穩(wěn)。

圖8 實(shí)時預(yù)測結(jié)果：（a）LINV；（b）FM；（c）SLO；（d）模型集成Figure 8 Real-time prediction results：（a）LINV；（b）FM；（c）SLO；（d）model integration

圖9 某露天煤礦滑坡的滑坡時間預(yù)測值與觀測值的殘差絕對值Figure 9 Absolute residual value between predicted and observed landslide time in an opencast coal mine

4 討論

由上述預(yù)測分析可知，利用滑坡加速變形階段的觀測數(shù)據(jù)預(yù)測滑坡時間具有如下規(guī)律：在加速變形階段前期，各模型的預(yù)測均值和最大后驗(yàn)估計值明顯偏離Tf，預(yù)測的置信區(qū)間跨度相對較寬，預(yù)測具有明顯的不確定性；在滑坡進(jìn)入臨滑階段后，各模型的預(yù)測均值和最大后驗(yàn)估計值向Tf收斂并在其附近波動，預(yù)測的置信區(qū)間寬度收縮或趨于穩(wěn)定。因此，在滑坡時間預(yù)測時，需判識滑坡所處的演化階段。許強(qiáng)等將滑坡演化的加速變形階段劃分為初加速、中加速和臨滑加速三個階段［23］。兩個案例的預(yù)測結(jié)果表明，在中加速變形階段，滑坡時間預(yù)測值開始向真實(shí)值收斂；進(jìn)入臨滑階段后，各模型預(yù)測出滑坡時間相對可靠。

在三個模型的預(yù)測中，根據(jù)預(yù)測的殘差絕對值E的動態(tài)演化規(guī)律，在兩個滑坡中FM 模型的E的值均相對較小，而LINV模型和SLO模型變化較大?？傮w而言，LINV模型在天荒坪抽水蓄能電站開關(guān)站滑坡的預(yù)測分析中E是三個模型中最小的，而在準(zhǔn)格爾旗某露天煤礦滑坡中的E是三個模型中最大的；而SLO 模型的情況恰好相反。這些模型的預(yù)測能力需要更多的滑坡案例進(jìn)行分析，以確定滑坡時間預(yù)測的最佳模型。雖然采用貝葉斯更新方法量化了預(yù)測不確定性，但是模型本身的預(yù)測性能仍對預(yù)測結(jié)果占主控作用。另外，集成預(yù)測的E基本位于三個單一模型預(yù)測的E之間，其演化相對平穩(wěn)，且其值大小處于一個相對較低的水平。因此，采用各模型預(yù)測結(jié)果的集成并進(jìn)行綜合評價對滑坡時間進(jìn)行分析是極為關(guān)鍵的。

針對上述分析，本研究對滑坡預(yù)測決策提出如下建議：

1）對比多模型預(yù)測的滑坡時間均值或最大后驗(yàn)估計值。在兩個案例中，二者之間的差距相對較小，可取二者之一或平均值為確定性預(yù)測結(jié)果。當(dāng)多時刻、多模型的預(yù)測均值或最大后驗(yàn)估計值在某一固定值附近波動，可認(rèn)為獲得較為可靠預(yù)測結(jié)果。

2）預(yù)測的95%置信區(qū)間作為滑坡時間窗口。窗口越寬或者隨時間波動性越大，代表預(yù)測不確定越高，邊坡破壞的可能性仍相對較低。窗口越窄或者隨時間波動性降低，可認(rèn)為當(dāng)前觀測數(shù)據(jù)已經(jīng)能較為清晰地反映滑坡的破壞趨勢，預(yù)測結(jié)果較為可靠。不同案例的時間窗口的寬度與變形特征、觀測噪音等因素相關(guān)，很難確定一個固定的閾值窗口，因此這窗口的寬窄是相對不同預(yù)測時刻而言的。此外，集成預(yù)測模型95%置信區(qū)間的下限值可為一個臨界預(yù)警閾值。

5 結(jié)論

本文針對滑坡時間實(shí)時預(yù)測存在不確定性的問題展開研究，采用貝葉斯理論對滑坡時間進(jìn)行概率預(yù)測，并對比三種模型的預(yù)測能力，最后集成三個的模型計算結(jié)果進(jìn)行預(yù)測決策，主要結(jié)論如下：

1）加速變形階段前期，各預(yù)測模型的預(yù)測值明顯偏離真實(shí)滑坡時間；進(jìn)入臨滑階段，預(yù)測值向真實(shí)滑坡時間收斂并在其附近波動，預(yù)測不確定性降低。

2）貝葉斯更新提供了預(yù)測均值、最大后驗(yàn)估計值、95%置信區(qū)間下限及置信區(qū)間寬度等多種預(yù)測指標(biāo)，能為滑坡預(yù)測決策和預(yù)警提供豐富的信息。

3）雖然本研究采用貝葉斯更新量化了預(yù)測不確定性，但受限于模型本身預(yù)測性能及其它難以考慮的影響因素，單一預(yù)測模型預(yù)測能力相對較弱，而多模型概率預(yù)測決策可有效地提高預(yù)測可靠性。