基于改進(jìn)灰色系統(tǒng)模型的庫岸邊坡變形分析

劉冠豪++楊品杰++葉險(xiǎn)峰

摘 要：邊坡變形監(jiān)測與預(yù)測預(yù)報(bào)是邊坡災(zāi)害防治的重要技術(shù)措施，也是目前研究的熱點(diǎn)和難點(diǎn)。該文針對邊坡變形監(jiān)測時間間隔不定、數(shù)據(jù)樣本少、預(yù)測精度要求高等問題，以某庫岸邊坡上的D02監(jiān)測點(diǎn)為例，采用拉格朗日插值算法進(jìn)行非等時距序列轉(zhuǎn)換，結(jié)合灰色系統(tǒng)原理，建立非等時距灰色預(yù)測模型進(jìn)行變形分析與預(yù)測。研究結(jié)果表明：（1）庫岸邊坡在該點(diǎn)的垂直位移變化經(jīng)歷了加速變形-減速變形-緩慢變形三個階段，其變形速率分別為-1.599mm/月、-0.103mm/月和-0.039mm/月，總體呈暫時趨穩(wěn)態(tài)勢；（2）建立的非等時距灰色預(yù)測模型有效解決了因天氣、施工和人為等因素影響，導(dǎo)致觀測周期難以保持一致，出現(xiàn)原始數(shù)據(jù)時距不等的問題，且預(yù)測結(jié)果具有很高的精度，有一定的推廣應(yīng)用價值。

關(guān)鍵詞：庫岸邊坡 拉格朗日插值 灰色系統(tǒng)模型 變形分析與預(yù)測

中圖分類號：G64 文獻(xiàn)標(biāo)識碼：A 文章編號：1674-098X（2015）01（a）-0051-04

隨著社會經(jīng)濟(jì)發(fā)展和能源需求增大，我國先后開工或建成了許多大型水利水電工程，如三峽大壩、小浪底水庫、南水北調(diào)工程、金沙江梯級電站等。大型水利水電工程庫區(qū)的水位與水文變化，導(dǎo)致庫岸山體穩(wěn)定性降低，誘發(fā)大量不穩(wěn)定邊坡、崩塌、滑坡、泥石流等地質(zhì)災(zāi)害，常造成山體變形破壞、河道堵塞、泥沙淤積、工程施工建設(shè)與運(yùn)營使用受損，甚至引起工程失效和潰壩等災(zāi)難，已成為社會關(guān)注的熱點(diǎn)。1963年意大利262 m高的瓦依昂拱壩因左岸山體滑坡而潰壩，造成1925人死亡。1976年美國93m高的提堂土壩潰決造成11人死亡、25000人無家可歸[1]。2008年汶川地震造成震區(qū)山體破壞，誘發(fā)大量次生地質(zhì)災(zāi)害，導(dǎo)致近2000座水庫受損[2]。因此，對庫岸邊坡進(jìn)行變形監(jiān)測，并根據(jù)監(jiān)測數(shù)據(jù)進(jìn)行邊坡變形分析與預(yù)測預(yù)報(bào)是一項(xiàng)重要的技術(shù)工作，對保障水利水電工程施工與運(yùn)營安全具有重要意義。

庫岸邊坡變形破壞受眾多不確定性因素影響，是一個復(fù)雜的非線性過程。其變形分析涉及變形數(shù)據(jù)處理與分析、變形物理解釋和變形預(yù)測預(yù)報(bào)等方面[3]。目前，國內(nèi)外學(xué)者對邊坡變形監(jiān)測數(shù)據(jù)處理與分析做了深入研究，取得了一定的理論研究成果，并發(fā)揮了實(shí)用效益[4]。早期基于巖土力學(xué)原理，結(jié)合工程經(jīng)驗(yàn)，采用有限元（1967）、離散元[5]、邊界元[6]和非連續(xù)變形分析[7]等方法，建立邊坡變形值與相關(guān)環(huán)境參數(shù)值的函數(shù)關(guān)系式進(jìn)行邊坡變形分析。但因受其力學(xué)機(jī)理粗淺或假設(shè)不合理、以及參數(shù)獲取等限制，致使計(jì)算結(jié)果與實(shí)際相差較大，預(yù)測精度受到很大影響。隨著概率與數(shù)理統(tǒng)計(jì)、時間序列分析、非線性理論和人工智能等現(xiàn)代數(shù)學(xué)與計(jì)算理論的誕生和發(fā)展，國內(nèi)外學(xué)者在此基礎(chǔ)上建立了眾多邊坡變形分析模型與方法[8-11]。每種變形分析模型都有一定的適用條件和局限性，單一的模型或方法已難以準(zhǔn)確預(yù)測出庫岸邊坡的變形情況，需要改進(jìn)或建立已有的模型方法，使其更好地分析邊坡實(shí)際變形狀態(tài)與趨勢[1]。庫岸邊坡變形的影響因素眾多，往往不能用單一變量表示，同時監(jiān)測時間常為不等間隔，因此灰色預(yù)測理論需在多變量非等步長研究上取得突破。該文針對邊坡變形監(jiān)測時間間隔不等，且影響變形因素眾多、數(shù)據(jù)樣本少、預(yù)測精度要求高等問題，引入拉格朗日插值算法進(jìn)行非等時距數(shù)據(jù)序列轉(zhuǎn)換，結(jié)合灰色系統(tǒng)原理，建立非等時距灰色預(yù)測模型，進(jìn)行邊坡變形分析與趨勢預(yù)測，并采用后驗(yàn)差方法檢驗(yàn)預(yù)測精度。

1 研究方法

1.1 拉格朗日插值法

1.1.1 插值算法原理

在庫岸邊坡變形研究中，可以用函數(shù)y=f（x）來表示其變形值在某點(diǎn)隨時間變化的數(shù)量關(guān)系。但這種函數(shù)關(guān)系沒有明顯解析式，通常通過實(shí)驗(yàn)觀測得到y(tǒng)=f（x）在某區(qū)間[a，b]上的一系列點(diǎn)的函數(shù)值yi=f（xi），即一張表（xi，yi），（i=0，1，2，…，n）。根據(jù)該函數(shù)表“插出”（構(gòu)造）一個既能反映函數(shù)f（x）特性、又計(jì)算簡單的函數(shù)（x）來近似f（x），并使得（xi）=f（xi），（i=0，1，2，…，n）。這種求函數(shù)近似式的方法稱為插值法，其數(shù)學(xué)表達(dá)式為：

設(shè)函數(shù)y=f（x）定義在區(qū)間[a，b]上，x0，x1，…，xn是[a， b]上取定的n+1個互異節(jié)點(diǎn)，且在這些點(diǎn)處的函數(shù)值f（x0），f（x1）， …，f（xn）為已知，即yi=f（xi），（i=0，1，2，…， n）。若存在一個f（x）的近似函數(shù)（x），滿足

（1）

則稱（x）為f（x）的一個插值函數(shù)。f（x）為被插值函數(shù)，點(diǎn)xi為插值節(jié)點(diǎn)，式（1）為插值條件，而誤差函數(shù)R（x）=f（x）-（x）稱為插值余項(xiàng)，區(qū)間[a，b]為插值區(qū)間，插值點(diǎn)在插值區(qū)間內(nèi)的稱為內(nèi)插，否則稱外插。

1.1.2 拉格朗日插值多項(xiàng)式

通過n+1個不同的節(jié)點(diǎn)（xi，yi），構(gòu)造一個特殊n次多項(xiàng)式li（x），使其在各節(jié)點(diǎn)上滿足：

（2）

即 （3）

由條件lk（xi）=0 （i≠k）可知：x0，x1，…，xk-1，xk+1，…，xn都是n次多項(xiàng)式lk（x）的零點(diǎn)，故可設(shè)

（4）

其中Ak為待定常數(shù)。由lk（xk）=1，可求得

（5）

將公式（5）代入公式（4），則

（6）

因此，以n+1個n次基本插值多項(xiàng)式（公式（6））為基礎(chǔ)的n次代數(shù)插值多項(xiàng)式為：

（7）

公式（7）稱為n次拉格朗日插值多項(xiàng)式，用于解決該文灰色系統(tǒng)模型中的非等時距問題。

1.2 灰色模型

1.2.1 灰色原理與方法

灰色系統(tǒng)理論以“部分信息已知、部分信息未知”的“小樣本、貧信息”不確定型系統(tǒng)為研究對象，用灰色數(shù)學(xué)模型將觀測數(shù)據(jù)序列看做隨時間變化的灰色量（過程）。通過累加生成和累減生成逐步使灰色量白化，從而建立相應(yīng)于微分方程解的模型并作出預(yù)測預(yù)報(bào)。其主要模型方法與計(jì)算步驟如下。endprint

1）生成原始數(shù)列：

（8）

2）累加生成序列：

（9）

并且X（0）與X（1）滿足：

（k=a，…，n） （10）

若生成的數(shù)列仍不滿足條件，則可進(jìn)行多次累加：

（11）

3）建立微分方程：

（12）

4）平滑差分：

（13）

5）最小二乘法求解：

（14）

公式（11）和公式（12）即為灰色預(yù)測的兩個基本模型。當(dāng)kn時，稱為模型預(yù)測值。

1.2.2 精度評定

采用后驗(yàn)差方法檢驗(yàn)灰色模型的精度，它由后驗(yàn)差比值C和小誤差概率P共同描述。

設(shè)由灰色模型得到：

（15）

計(jì)算殘差：

（16）

記原始數(shù)列及殘差數(shù)列的方差分別為和，則

（17）

式中

（18）

計(jì)算后驗(yàn)差比值C和小誤差概率P：

（19）

模型精度等級=max{所在的級別，所在的級別}，其判別式如下（表1）。

1.3 改進(jìn)的灰色模型

1.3.1 等時距變換

灰色模型要求原始數(shù)據(jù)必須是等時距，但由于沉降觀測持續(xù)時間長，易受到天氣、施工及人為因素影響，常導(dǎo)致觀測周期難以保持一致，從而出現(xiàn)原始數(shù)據(jù)時距不等。因此，在按等時距灰色理論建模型之前，需要對實(shí)測數(shù)據(jù)進(jìn)行等時距變換，其具體過程如下。

取庫岸邊坡沉降觀測斷面在觀測時刻的累計(jì)沉降作為原始數(shù)據(jù)序列，記為

（20）

各時段時間間隔為ti=ti+1-ti，⊿tj=tj+1-tj，i，j{1，2，…n}。當(dāng)i≠j時，⊿ti≠tj，即表示各時段的時間間隔不相等。

平均時間間隔：

（21）

計(jì)算等時間間隔點(diǎn)的累計(jì)灰沉降值XL（1）（k），（k=1，2，…，n）：

（22）

式中，i的取值應(yīng)滿足，從而得到等時距數(shù)據(jù)序列：

（23）

1.3.2 非等時距建模

將擬合成一階線性微分方程，即

（24）

式中a，u為待定常數(shù)。a為發(fā)展系數(shù)，其大小反映了序列的增長速度，u為灰色作用量。按照最小二乘法原理求解可得：

（k為時間間隔點(diǎn)） （25）

取時所對應(yīng)的時間t=0，則k時間間隔點(diǎn)所對應(yīng)的時間t=（k-1）。故t時刻所對應(yīng)的預(yù)測沉降值為：

（26）

聯(lián)合公式25和公式26，得到非等時距預(yù)測模型：

（27）

由于非等時距灰色模型是經(jīng)過拉格朗日插值與等時距變換而成，故其精度評定方法與等時距灰色模型相同。

2 實(shí)例分析

2.1 形變分析

以湖南省張家界市慈利縣境內(nèi)的某水庫壩址右岸山體一邊坡為研究對象，該邊坡存在變形失穩(wěn)風(fēng)險(xiǎn)，直接影響水庫大壩及其下游5km處的城鎮(zhèn)安全。自2000年以來，在該邊坡體上布設(shè)了沉降監(jiān)測點(diǎn)，并采用GPS、水準(zhǔn)測量和變形監(jiān)測方法進(jìn)行各監(jiān)測點(diǎn)的形變監(jiān)測，以分析邊坡的變形規(guī)律與趨勢?，F(xiàn)以該邊坡上的D02號監(jiān)測點(diǎn)的垂直位移變化值為例進(jìn)行邊坡變形分析與預(yù)測（表2）。

利用公式（2）～（7）和公式（20）～ （23），進(jìn)行等時距變化處理（表3）。

為了便于分析計(jì)算結(jié)果，則將D02監(jiān)測點(diǎn)垂直位移變形的等時距變換值以時間為橫坐標(biāo)、變形值為縱坐標(biāo)，繪制成變形過程線圖（圖1），以直觀反映出觀測點(diǎn)垂直位移變形的規(guī)律、趨勢和幅度。

從圖1和表3可知：D02監(jiān)測點(diǎn)垂直位移在前5個月內(nèi)沉降顯著，呈加速變形趨勢，變形速率為-1.599 mm/月；在第5～49月期間，D02監(jiān)測點(diǎn)沉降速度減少，呈減速變形趨勢，變形速率為-0.103 mm/月；在第49～62月期間，D02監(jiān)測點(diǎn)沉降速度減少，呈緩慢變形趨勢，變形速率為-0.039 mm/月。總體而言，邊坡體上的D02監(jiān)測點(diǎn)在監(jiān)測后的5年內(nèi)經(jīng)歷了加速變形、減速變形和緩慢變形等階段，并趨于穩(wěn)定。

2.2 模型預(yù)測

由于D02監(jiān)測時間間隔各不相等，而等時距變換后的時間為非整月或半月的月份，即非規(guī)則月份，不便于模型預(yù)測和應(yīng)用分析?，F(xiàn)以0.5個月為單位，經(jīng)拉格朗日插值，計(jì)算各規(guī)則監(jiān)測時間間隔，則可以建立插值后的時間序列：{t}={0， 0.5， 1， 4， 5， 7.5， 9， 11.5， 14.5， 17， 18.5， 19.5， 22.5， 28， 32.5， 35.5，37.5， 40.5， 44， 47.5， 50.5，53，56，60.5， 62}。

由于該時間序列觀測值仍不等時距，則采用灰色非等時距建模方法，進(jìn)行累加生成和累減生成變換處理，且保證處理用于灰色非等時距預(yù)測模型的原始序列值滿足非負(fù)性和遞增性要求。由表3可知，經(jīng)拉格朗日插值后的觀測數(shù)據(jù)仍不能滿足非負(fù)性要求，因此，需要在建模前進(jìn)行如下處理。

注：平均時間間隔=62/24=2.58（月）

①移軸。

尋找數(shù)列中絕對值最大的負(fù)數(shù)（-13.3），然后取該負(fù)數(shù)的絕對值，并將該絕對值設(shè)為min，再對數(shù)列中每個數(shù)值加上min，得到數(shù)列{x（0）（t）} ={11.00， 6.56， 5.03， 5.27， 4.26， 4.31， 4.08， 3.47， 2.71， 2.40， 2.40， 2.38， 2.27， 2.06， 1.83， 1.69， 1.50， 1.79， 1.16， 0.51， 0.65， 1.72， 0.40， 0.40， 0.00}。

由于x（0）（2）< x（0）（1），表明移軸后數(shù)列仍然不滿足建模所需的遞增條件，不能建立一階灰色非等距模型，因此，需要對移軸后的數(shù)列進(jìn)行累加生成。

②累加生成。

一次累加生成，結(jié)果為：{x（1）（t）} ={11，17.56， 22.59， 27.86， 32.12， 36.42， 40.50， 43.97， 46.68， 49.08， 51.48， 53.86， 56.12， 58.18， 60.01， 61.70， 63.20， 64.99， 66.15， 66.66， 67.30， 69.02， 69.42， 69.82， 69.82}。該數(shù)列x（1）（23）=x（1）（24），仍不滿足建模所需的遞增條件，則需對其進(jìn)行二次累加生成。

二次累加生成，結(jié)果為：{x（2）（t）} ={11， 28.56， 51.15， 79.01， 111.12， 147.55， 188.05， 232.02， 278.70， 327.77， 379.25， 433.11， 489.23， 547.41， 607.42， 669.13， 732.33， 797.32， 863.46， 930.12， 997.42， 1066.44， 1135.87， 1205.69， 1275.51}。該數(shù)列已滿足非負(fù)、遞增的建模要求。

③二階建模。

利用二次累加生成的數(shù)列{x（2）（t）}，進(jìn)行二階灰色非等距建模，可得：

最小二乘原理求解，可得：a=-0.0206，u=-8.3530，則對應(yīng)的非等距灰色預(yù)測模型為：

用公式（28）計(jì)算D02監(jiān)測點(diǎn)各非等時距灰色預(yù)測值（表4、圖2）。

注：表中“實(shí)際值”為監(jiān)測點(diǎn)實(shí)際監(jiān)測時間經(jīng)等時距插值處理后得到的垂直位移形變值；“預(yù)測值”為按監(jiān)測點(diǎn)時間序列進(jìn)行非等時距灰色模型計(jì)算的垂直位移形變值；“-”表示沒有實(shí)際觀測值；殘差值=實(shí)際值-觀測值；相對誤差=|差值/實(shí)際值|。

由表4和圖2知：D02監(jiān)測點(diǎn)的實(shí)際值和預(yù)測值比較吻合，殘差值圍繞y=0上下波動，初步預(yù)測模型預(yù)測良好。

2.3 精度檢驗(yàn)

由公式（16）～（18），計(jì)算可得：

=-0.064，=-10.66，=5.14，=0.24

計(jì)算后驗(yàn)差比值：，即C<0.35。

計(jì)算小誤差概率：1，即0.95

模型預(yù)測精度等級=max{P所在的級別，C所在的級別}=1級（好）。因此，非等距灰色模型在邊坡變形預(yù)測中具有很高精度，可以滿足實(shí)際工作需要。例如，該模型預(yù)測D02監(jiān)測點(diǎn)在2007年7月9日的變形值為-15.4mm，而該點(diǎn)在2007年8月16日采集的實(shí)際數(shù)據(jù)計(jì)算出的變形值為-14.6mm，殘差值為0.8mm、相對誤差為5.48%，滿足邊坡變形分析的精度要求。

3 結(jié)論與討論

1）D02監(jiān)測點(diǎn)近5年的觀測結(jié)果表明，庫岸邊坡在該點(diǎn)的垂直位移變化經(jīng)歷了加速變形-減速變形-緩慢變形三個階段，在監(jiān)測至第5個月、第5～49個月和第49～62個月的變形速率分別為-1.599mm/月、-0.103mm/月和-0.039mm/月，故該邊坡總體呈暫時趨穩(wěn)態(tài)勢。

2）采用拉格朗日插值算法進(jìn)行非等時距序列轉(zhuǎn)換，結(jié)合灰色系統(tǒng)原理，建立非等時距灰色預(yù)測模型進(jìn)行邊坡變形分析與預(yù)測，有效解決了因天氣、施工和人為等影響因素而導(dǎo)致觀測周期難以保持一致出現(xiàn)原始數(shù)據(jù)時距不等的問題，且預(yù)測結(jié)果具有很高的精度，能夠滿足實(shí)際需要。

3）隨著預(yù)測時間的延長，非等時距灰色模型的預(yù)測精度會有所降低，需要在實(shí)踐應(yīng)用中對模型反復(fù)進(jìn)行試驗(yàn)分析，厘定出滿足預(yù)測精度要求的數(shù)據(jù)序列的有效時間長度，以提高模型的精度與實(shí)用性。

參考文獻(xiàn)

[1] 高玉云.基于小波理論的變形監(jiān)測數(shù)據(jù)處理模型的研究[D].西南石油大學(xué)，2011.

[2] 孔亞平，韓用順，朱穎彥.震區(qū)公路走廊次生地質(zhì)災(zāi)害監(jiān)測與評估地圖集[M]. 人民交通出版社，2013.

[3] 王孟穹.模糊分析方法及時序分析模型在變形監(jiān)測中應(yīng)用[D].西南交通大學(xué)，2009.

[4] 王忠學(xué)，梁東偉.基于小波神經(jīng)網(wǎng)絡(luò)在滑坡變形監(jiān)測預(yù)報(bào)中的應(yīng)用[J].測繪與空間地理信息，2012，35（11）：190-192.

[5] Cundall P. A. A computer model for simulating progressive large scale movements in blocky rock systems[C]//ISRM Symp.， Nancy， France， Proc.2.1971：129-136.

[6] Crouch S.L. Solution of plane elasticity problems by the displacement discontinuity method. I. Infinite body solution[J]. International Journal for Numerical Methods in Engineering， 1976， 10（2）：301-343.

[7] Shi G.H. & Goodman R.E. Two dimensional discontinuous deformation analysis[J]. International Journal for Numerical and Analysis Methods in Geomechanics， 1985，9（6）：541-556.

[8] 石雙忠，岳東杰，梅紅.時間序列分析在變形監(jiān)測數(shù)據(jù)處理中的應(yīng)用[J].工程勘察，2004（3）：59-61.

[9] 孫樹海，曹蘭柱，張立新.露天礦邊坡穩(wěn)定性的模糊綜合評判[J].遼寧工程技術(shù)大學(xué)學(xué)報(bào)，2007，26（2）：177-179.

[10] Prakash N. & Manikandan S.A.， et al. Prediction of biosorption efficiency for the removal of copper using artificial neural networks[J].Hazard Mater，2008：1268-1275.

[11] 付建軍，邱山鳴，趙海斌.基于灰色關(guān)聯(lián)度的邊坡穩(wěn)定影響因素分析[J].長江科學(xué)院院報(bào)，2011，28（1）：53-58.