基于最大Lyapunov指數(shù)的高邊坡安全監(jiān)控優(yōu)化模型

牛景太

（1.南昌工程學(xué)院 水利與生態(tài)工程學(xué)院，南昌 330099；2.河海大學(xué) 水利水電學(xué)院，南京 210098）

高邊坡是工程建設(shè)的主要地質(zhì)環(huán)境和工程承載體，高邊坡發(fā)生滑坡將造成嚴(yán)重的后果，因此，高邊坡安全監(jiān)控和預(yù)警已成為當(dāng)前工程建設(shè)極其關(guān)注的問題［1].影響高邊坡穩(wěn)定的因素十分復(fù)雜，高邊坡系統(tǒng)的演化過程表現(xiàn)出復(fù)雜的非線性動(dòng)力學(xué)特性［2].研究表明［3]高邊坡穩(wěn)定具有混沌特性，難以用確定性分析方法對(duì)其進(jìn)行模擬分析.而依據(jù)長期積累的安全監(jiān)測資料，對(duì)高邊坡在多因素復(fù)雜環(huán)境作用下的穩(wěn)定性進(jìn)行監(jiān)控和預(yù)測越來越受到重視，并進(jìn)行了廣泛研究.然而，已有成果［4－6]多建立在整個(gè)監(jiān)測資料時(shí)間序列之上，沒有考慮高邊坡動(dòng)力學(xué)系統(tǒng)在外界或其自身影響下是否發(fā)生改變.

實(shí)際工程中，高邊坡經(jīng)常受到例如爆破、地震等外界因素的擾動(dòng)而使其演變過程產(chǎn)生突變，致使整個(gè)監(jiān)測資料時(shí)間序列具有以突變點(diǎn)為分界的跳躍性，如果不考慮這些突變，而把監(jiān)控預(yù)測模型建立在整個(gè)監(jiān)測資料時(shí)間序列之上，將會(huì)顯著影響監(jiān)控模型的精準(zhǔn)度.因此，有效辨識(shí)突變測值位置，消除或削弱位移突變對(duì)測值序列整體數(shù)值特征的影響，是提高高邊坡位移監(jiān)控模型擬合和預(yù)測精度的關(guān)鍵問題之一.

針對(duì)以上問題，對(duì)高邊坡監(jiān)控模型進(jìn)行研究，基本思路如下：基于高邊坡系統(tǒng)演化過程中的非線性動(dòng)力學(xué)特性，對(duì)監(jiān)測數(shù)據(jù)序列進(jìn)行相空間重構(gòu)，采用動(dòng)力學(xué)結(jié)構(gòu)突變診斷方法與云模型理論，診斷出整個(gè)原始監(jiān)測時(shí)間序列中的突變點(diǎn)；由于突變的產(chǎn)生，高邊坡系統(tǒng)也由一個(gè)動(dòng)力學(xué)屬性轉(zhuǎn)變?yōu)榱硪粋€(gè)新的動(dòng)力學(xué)屬性，以最近一次突變之后相對(duì)穩(wěn)定的動(dòng)力系統(tǒng)結(jié)構(gòu)，考慮降雨和時(shí)效等因素的影響建立高安全監(jiān)控模型，以期提高預(yù)測結(jié)果的準(zhǔn)確性；通過算例證明了該方法具有更高的預(yù)測精度，對(duì)于高邊坡位移的實(shí)時(shí)監(jiān)測和預(yù)警都有重要意義.

1 高邊坡位移突變辨識(shí)方法

高邊坡系統(tǒng)是一個(gè)開放的、耗散的以及復(fù)雜的非線性動(dòng)力系統(tǒng)［2，7]，具有混沌特性.因此，利用動(dòng)力學(xué)理論對(duì)高邊坡在多因素復(fù)雜環(huán)境作用下的動(dòng)力學(xué)結(jié)構(gòu)進(jìn)行分析，深入了解高邊坡的位移演變過程.

1.1 位移監(jiān)測序列的相空間重構(gòu)

令｛ε′（t），t＝1，2，…，N｝為一維位移監(jiān)測時(shí)間序列，混沌時(shí)間序列的預(yù)測建立在高維相空間重構(gòu)的基礎(chǔ)上.為此，需對(duì)ε′（t）進(jìn)行相空間重構(gòu)，將原有序列延拓成m維相空間的一個(gè)相形分布，其表達(dá)式為

式中，τ＝kΔt為時(shí)間延遲，k（k＝1，2，…，n）為延遲參數(shù)，Δt為采樣間隔時(shí)間；m為嵌入維數(shù)；Xi為m維相空間中的相點(diǎn)；M＝N－（m－1）k為相點(diǎn)個(gè)數(shù).

其中，嵌入維數(shù)m的確定與關(guān)聯(lián)維數(shù)D2有關(guān)，可根據(jù) Grassberger－Procaccia方法［8]進(jìn)行選?。欢鴷r(shí)間延遲τ可由互信息［9]確定.據(jù)此，當(dāng)嵌入維數(shù)m和時(shí)間延遲τ確定后，虛殘差序列即可重構(gòu)為如式（1）所示的〈N－k（m－1）〉×m維的向量矩陣.

1.2 高邊坡位移測值序列突變點(diǎn)辨識(shí)

在動(dòng)力學(xué)上，高邊坡內(nèi)部動(dòng)力結(jié)構(gòu)發(fā)生演化或外界的擾動(dòng)過大導(dǎo)致位移突變，使高邊坡的狀態(tài)在相空間中不再趨向于原來的吸引子而是趨向于新的吸引子，從而整個(gè)系統(tǒng)也由一種動(dòng)力學(xué)屬性轉(zhuǎn)變?yōu)榱硪环N動(dòng)力學(xué)屬性.因此，基于混沌動(dòng)力學(xué)中關(guān)于Lyapunov指數(shù)的理論，提出診斷高邊坡位移時(shí)間序列突變的動(dòng)力學(xué)最大Lyapunov指數(shù)法.最大Lyapunov指數(shù)能定量地描述復(fù)雜系統(tǒng)相空間相鄰軌道呈指數(shù)發(fā)散或收斂的性質(zhì)，是描述系統(tǒng)混沌動(dòng)力學(xué)特性的重要參數(shù)之一，與Kolmogorov熵、Hausdorff維、信息維、關(guān)聯(lián)維、高階信息維等相比，最大Lyapunov指數(shù)能夠提供混沌系統(tǒng)更為有用的動(dòng)力學(xué)診斷［10].本文采用小數(shù)據(jù)量法計(jì)算最大Lyapunov指數(shù)，其具體計(jì)算步驟如下：

1）對(duì)位移監(jiān)測序列的相空間重構(gòu)，重構(gòu)后的相空間記為：｛Xj，j＝1，2，…，M｝，M為相點(diǎn)的數(shù)量.

2）找相空間中每個(gè)點(diǎn)Xj的最近臨近點(diǎn)X＾j，p為時(shí)間軌道平均周期，計(jì)算dj（0）的值：

式中，‖·‖表示向量2范數(shù)，inf（·）表示自變量在某特定域中取下界值，即最小值.

3）對(duì)相空間中的每個(gè)點(diǎn)Xj，計(jì)算出該鄰域點(diǎn)對(duì)的第i個(gè)離散時(shí)間步后的距離dj（i）為

4）由公式（2）和（3），前進(jìn)距離dj（i）與dj（0）有以下近似關(guān)系：

式中，Δt為觀測序列的采樣間隔或步長；i為相點(diǎn)沿時(shí)間軌道的滑動(dòng)步長序數(shù)；λ為最大Lyapunov指數(shù).對(duì)式（4）兩邊取對(duì)數(shù)：

考慮到局部計(jì)算的影響，最大Lyapunov指數(shù)的最后經(jīng)驗(yàn)公式為

式中，〈·〉表示按相空間的點(diǎn)數(shù)求平均.

根據(jù)Lyapunov穩(wěn)定判據(jù)［11]，當(dāng)λ≤0時(shí)，表明系統(tǒng)處于穩(wěn)定狀態(tài)；當(dāng)0＜λ＜∞時(shí)，表明監(jiān)測點(diǎn)變形演化系統(tǒng)處于不穩(wěn)定狀態(tài)，且其值越大不穩(wěn)定度或失穩(wěn)度越高.

如果某段測值序列最大Lyapunov指數(shù)為正，但其值較小，屬于微弱不穩(wěn)定，那么本段測值序列不能認(rèn)為發(fā)生突變，因此，通過以上方法求出的最大Lyapunov指數(shù)雖然具有一定的區(qū)分潛在動(dòng)力學(xué)的能力，但是還不能作為辨識(shí)高邊坡測值序列突變的重要標(biāo)志，必須引入不穩(wěn)定度［12]的概念.在重構(gòu)的相空間｛Xj，j＝1，2，…，M｝里，選初始時(shí)刻t0為相點(diǎn)的基準(zhǔn)點(diǎn)，在余下的相點(diǎn)中按與基準(zhǔn)點(diǎn)距離最小的原則找出變動(dòng)點(diǎn)X′（t0＋iΔt），令L（ti）為基準(zhǔn)點(diǎn)X（t0＋iΔt）與變動(dòng)點(diǎn)X′（t0＋iΔt）之間的長度，λi為第i步的 Lyapunov指數(shù)，則：

對(duì)上式進(jìn)行變形，并令Δt＝1，則：

在突變點(diǎn)的辨析過程中，常常關(guān)心數(shù)值較大的Lyapunov指數(shù)，因此，為簡化計(jì)算，舍棄數(shù)值小的Lyapunov指數(shù)，用最大Lyapunov指數(shù)λ代替λi，從而，可得不穩(wěn)定度的表達(dá)式為

式（8）和（9）中，R為不穩(wěn)定度，ΔL為擾動(dòng)效應(yīng)量經(jīng)單位時(shí)間后的變化量，即擾動(dòng)效應(yīng)量變化的平均速率.

式（9）中的R只表明了動(dòng)力學(xué)時(shí)間序列上不穩(wěn)定度，但不能表征高邊坡動(dòng)力系統(tǒng)的突變轉(zhuǎn)異程度，本文借助于云模型理論［13]來實(shí)現(xiàn)這一目的.云模型理論認(rèn)為：云是由若干云滴組成，云滴是某個(gè)定性概念的一次隨機(jī)實(shí)現(xiàn)，多次產(chǎn)生的云滴可以綜合反映這個(gè)定性概念的整體特征.對(duì)于不穩(wěn)定度序列｛Ri｝，設(shè)Yi為Ri的隸屬函數(shù)，則（Ri，Yi）為一個(gè)云滴.從而，以N個(gè)云滴在數(shù)域空間的定量位置及每個(gè)云滴代表概念的確定度為輸入量，采用逆向云發(fā)生器算法［14]，求解云模型的3個(gè)特征值ER、En、He，即可得到云期望曲線方程.據(jù)此，在給定高邊坡變形轉(zhuǎn)異顯著性水平α下，可獲得不穩(wěn)定度R的臨界值Rα，從而基于動(dòng)力學(xué)結(jié)構(gòu)突變的高邊坡變形突變診斷判據(jù)為

1.3 高邊坡測值序列突變點(diǎn)辨識(shí)步驟

綜合以上分析，高邊坡測值序列突變點(diǎn)的辨析步驟歸納如下：1）在位移監(jiān)測序列中取一寬度為n的滑動(dòng)窗口W，從點(diǎn)n至N（從左至右）依次滑動(dòng)；2）對(duì)每個(gè)寬度為W的窗口進(jìn)行嵌入空間上的動(dòng)力學(xué)相空間重構(gòu)，計(jì)算各個(gè)滑動(dòng)窗口的最大Lyapunov指數(shù)｛λi｝；3）根據(jù)每個(gè)窗口的最大Lyapunov指數(shù)｛λi｝，計(jì)算每個(gè)滑動(dòng)窗口的不穩(wěn)定度｛Ri｝；4）利用云模型對(duì)｛Ri｝采用逆向云發(fā)生器算法，求出云期望曲線方程，然后在給定高邊坡變形突變顯著性水平α后，求得不穩(wěn)定度臨界值｛Rα｝；5）利用判據(jù)式（10）、（11）判斷高邊坡變形突變點(diǎn)位置.

2 構(gòu)建考慮突變影響的高邊坡位移安全監(jiān)控模型

通過前文建立突變診斷方法，將高邊坡監(jiān)測序列分解為多個(gè)動(dòng)力結(jié)構(gòu)穩(wěn)定的子序列，以最近一次突變之后相對(duì)穩(wěn)定的動(dòng)力系統(tǒng)結(jié)構(gòu)，考慮降雨和時(shí)效等因素的影響建立高邊坡位移安全監(jiān)控模型，以期提高預(yù)測結(jié)果的準(zhǔn)確性.以高邊坡的變形效應(yīng)量為例，其位移監(jiān)控模型可表示為

式中，δ為高邊坡變形；δθ為高邊坡變形的時(shí)效分量；δP為高邊坡變形的降雨分量.

粘彈塑性力學(xué)［15]認(rèn)為，高邊坡的變形是巖體流變的結(jié)果，其失穩(wěn)主要表現(xiàn)為隨時(shí)間持續(xù)發(fā)展的時(shí)效變形，當(dāng)變形累積到一定程度時(shí)邊坡將進(jìn)入加速變形階段，直至發(fā)生破壞或滑坡，如圖1所示.

圖1 巖體應(yīng)變（位移）－時(shí)間曲線圖

從圖1可以看出，巖體變形從開始直至破壞共經(jīng)歷了3個(gè)階段，第Ⅰ階段（AB段）為減速蠕變階段；第Ⅱ階段（BC段）為穩(wěn)定蠕變階段；第Ⅲ階段（CD段）為加速蠕變階段.其中，AB和BC段可采用對(duì)數(shù)巖石經(jīng)驗(yàn)流變公式［16]進(jìn)行擬合，即：

而CD段可采用冪數(shù)型的巖石經(jīng)驗(yàn)流變公式［17]進(jìn)行擬合，即：

王燕茹：今年3月份，黃宇因?yàn)榇蛭冶恍芯?日，那時(shí)候想著先休息一段時(shí)間。 6月份，黃宇的親戚威脅說，黃宇不會(huì)有事，黃道龍也不會(huì)受到什么影響，到時(shí)候會(huì)讓我在揚(yáng)州待不下去。我怕被打擊報(bào)復(fù)，就跑到北京，選擇在北京中紀(jì)委和國家信訪局繼續(xù)反映情況。

式（13）～（14）中：c0、c1、c2、c3為回歸系數(shù)；δθ1、δθ2為時(shí)變位移；t為自初始監(jiān)測開始的累計(jì)天數(shù).

而高邊坡變形的降雨分量可采用有效降雨量方法［18]計(jì)算，具體做法是選取觀測之日前一個(gè)月內(nèi)日降雨量，并與觀測之日間隔天數(shù)乘以0.8進(jìn)行依次折減并累計(jì)求和，作為降雨因子P的值，即

式中，T為觀測當(dāng)月的天數(shù)；Ri為距觀測當(dāng)日第i天的降雨量.

綜上分析，高邊坡位移監(jiān)控模型可表示為

式中，c4為降雨分量的回歸系數(shù)，其余符號(hào)含義同前.

當(dāng)最近一次突變發(fā)生較晚時(shí)，將會(huì)導(dǎo)致突變后的子序列測值數(shù)據(jù)較少，從而影響模型監(jiān)控和預(yù)測精度，此時(shí)，應(yīng)增大突變后期觀測頻率，增加觀測次數(shù)，這既是特殊情況的應(yīng)對(duì)措施，又為提高監(jiān)控模型的預(yù)測精度奠定基礎(chǔ).

3 工程實(shí)例

本節(jié)以文獻(xiàn)［5]中的實(shí)際工程為例進(jìn)行分析.某水電站工程區(qū)為典型的深切V型峽谷，相對(duì)高差1 500～1 700m.左岸為1km以上的高陡邊坡，基巖裸露，坡度為55～70°.拱壩左岸開挖邊坡高度達(dá)到500m.左岸拱壩開挖邊坡、左岸導(dǎo)流洞出口開挖邊坡都存在變形拉裂巖體，工程邊坡的安全穩(wěn)定性問題十分突出.因此，除在左岸壩肩、導(dǎo)流洞出水口等邊坡上布置了大量的內(nèi)觀監(jiān)測儀器外，在重點(diǎn)區(qū)域也布置了大量外觀監(jiān)測點(diǎn).選取電站左岸邊坡一個(gè)典型外觀測點(diǎn)TP12－1的觀測數(shù)據(jù)（表1）進(jìn)行分析.

表1 TP12－1位移觀測數(shù)據(jù)

續(xù)表1 TP12－1位移觀測數(shù)據(jù)

首先取寬度W＝5的滑動(dòng)窗口，并對(duì)其進(jìn)行相空間重構(gòu)，計(jì)算中取嵌入維數(shù)m＝3，延遲時(shí)間τ＝2；然后按照動(dòng)力學(xué)結(jié)構(gòu)突變的辨識(shí)方法并結(jié)合云模型進(jìn)行突變診斷.得到的不穩(wěn)定度R，其結(jié)果如圖2所示.

圖2 高邊坡監(jiān)測序列不穩(wěn)定度R值曲線圖

取突變顯著水平α＝0.05，從而得到高邊坡位移監(jiān)測數(shù)據(jù)序列中的突變點(diǎn).診斷結(jié)果表明，整個(gè)監(jiān)測序列共被檢測出2個(gè)動(dòng)力結(jié)構(gòu)突變點(diǎn)，突變時(shí)間為2006年1月和2006年6月.

考慮降雨和時(shí)效因素等因素的影響，利用高邊坡位移監(jiān)控模型（式5）分別對(duì)整個(gè)實(shí)測邊坡時(shí)間序列和2006年6月以后的監(jiān)測數(shù)據(jù)子序列進(jìn)行擬合，從而得到該模型表達(dá)式為

利用以上兩安全監(jiān)控模型分別對(duì)2007年9月到2008年2月的邊坡位移進(jìn)行預(yù)測，兩種情況下預(yù)測結(jié)果及殘差計(jì)算結(jié)果見表2.由表2可以看出，兩者的相對(duì)誤差絕對(duì)值的平均值分別為2.34%、1.20%，考慮突變影響以最近一次位移突變后的監(jiān)測資料建立的監(jiān)控預(yù)測模型預(yù)測精度顯著提高.

表2 兩種情況下高邊坡變形預(yù)測結(jié)果表

4 結(jié) 論

1）研究并提出了基于最大Lyapunov指數(shù)的高邊坡突變診斷方法，并利用該方法對(duì)某高邊坡的監(jiān)測時(shí)間序列進(jìn)行了結(jié)構(gòu)突變動(dòng)力學(xué)分析和診斷，結(jié)果表明，利用該方法能夠?qū)崿F(xiàn)高邊坡的動(dòng)態(tài)診斷.

2）通過動(dòng)力學(xué)突變診斷，把邊坡動(dòng)力系統(tǒng)以突變點(diǎn)為界分割為多個(gè)動(dòng)力系統(tǒng)屬性，基于最近的相對(duì)穩(wěn)定的動(dòng)力系統(tǒng)屬性，提出了利用最近一次位移突變后的監(jiān)測資料建立監(jiān)控預(yù)測模型.

3）以最近一次突變之后相對(duì)穩(wěn)定的動(dòng)力系統(tǒng)結(jié)構(gòu)，基于實(shí)測資料，考慮降雨和時(shí)效等因素的影響建立高邊坡位移安全監(jiān)控模型，據(jù)此定量評(píng)價(jià)高邊坡的安全狀態(tài)、監(jiān)控高邊坡的位移性狀、發(fā)現(xiàn)高邊坡的異常征兆，對(duì)于高邊坡位移的實(shí)時(shí)監(jiān)測和預(yù)警都有重要意義.

［1]黃潤秋.中國西南巖石高邊坡的主要特征及其演化［J].地球科學(xué)進(jìn)展，2005（3）：292－297.

［2]付義祥，劉志強(qiáng).邊坡位移的混沌時(shí)間序列分析方法及應(yīng)用研究［J].武漢理工大學(xué)學(xué)報(bào)：交通科學(xué)與工程版，2003，27（4）：473－476.

［3]吳中如，潘衛(wèi)平.應(yīng)用Lyapunov指數(shù)研究巖土邊坡的穩(wěn)定判據(jù)［J].巖石力學(xué)與工程學(xué)報(bào)，1997，16（3）：217－223.

［4]Moody J，Darken C.Fast learning in Networks of Locally－tuned Processing Units［J].Neural Computation，1989，1（2）：281－294.

［5]談小龍，徐衛(wèi)亞.高邊坡變形非線性時(shí)變統(tǒng)計(jì)模型研究［J].巖土力學(xué)，2010，31（5）：1633－1650.

［6]鄭東健，顧沖時(shí)，吳中如.邊坡變形的多因素時(shí)變預(yù)測模型［J].巖石力學(xué)與工程學(xué)報(bào)，2005，24（17）：3180－3184.

［7]周創(chuàng)兵，陳益峰.基于相空間重構(gòu)的邊坡位移預(yù)測［J].巖土力學(xué)，2000，21（3）：205－208.

［8]Albano A M，Muench J，Schwartz C.Singular－value Decomposition and the Grassberger－Procaccia Algorithm［J].Phys Rev A，1998，（38）：3017－3026.

［9]Fraser A M.Information and Entropy in Strange Attractors［J].IEEE Trans Inform Theory，1989，35（2）：245－262.

［10]Wolf A，Swift J B，Swinney H L，et al.Determining Lyapunov Exponents from a Time Series［J].Physica D，1985，16（3）：285－317.

［11]吳中如，潘衛(wèi)平.應(yīng)用Lyapunov指數(shù)研究巖土邊坡的穩(wěn)定判據(jù)［J].巖石力學(xué)與工程學(xué)報(bào)，1997，16（3）：217－223.

［12]顧沖時(shí).大壩與壩基安全監(jiān)控理論和方法及其應(yīng)用［M].南京：河海大學(xué)出版社，2006：177－181.

［13]李德毅，杜 鵝.不確定性人工智能［M].北京：國防工業(yè)出版社，2005.

［14]呂輝軍，王 曄.逆向云在定性評(píng)價(jià)中的應(yīng)用［J].計(jì)算機(jī)學(xué)報(bào)，2003，26（8）：1009－1014.

［15]賈乃文.粘塑性力學(xué)及工程應(yīng)用［M].北京：地震出版社，2000.

［16]VOIGHT B.A Relation to Describe Rate－dependent Material failure［J].Science，1989，243：200－203.

［17]SAITO M.Forecasting the Time of Occurrence of Slope Failure［C].Proceedings of 6th International Congress of Soil Mechanics and Foundation Engineering.Montreal：［s.n.]，1965：537－541.

［18]Zêzere J L，Trigo R M，Trigo I F.Shallow and Deep Landslides Induced by Rainfall in the Lisbon Region（Portugal）：Assessment of Relationships with the North Atlantic Oscillation［J].Natural Hazards and Earth System Sciences，2005，5：331－344.