基于降維分析的大壩安全特征提取方法及其應(yīng)用研究

方衛(wèi)華，王潤英，孫一清

(1.水利部南京水利水文自動(dòng)化研究所，江蘇 南京 210008 ;2.河海大學(xué)水利水電學(xué)院，江蘇 南京 210098)

近年的潰壩失事事件表明，大壩安全風(fēng)險(xiǎn)依然存在。對(duì)于運(yùn)行期的大壩安全問題，目前的安全特征提取方法主要包括基于實(shí)測資料的分析方法和數(shù)值模擬方法。前者主要包括特征值統(tǒng)計(jì)、回歸分析和監(jiān)控指標(biāo)擬定等，后者主要采用數(shù)值模擬法、剛體極限平衡法[1]和強(qiáng)度折減法等獲取強(qiáng)度和穩(wěn)定參數(shù)。盡管上述研究已取得很多研究成果，但前者依然受到小樣本、多重共線性、非正態(tài)分布、測值誤差以及分析過程和結(jié)果不直觀等問題的困擾，而后者由于受計(jì)算復(fù)雜度加卸載過程等因素的限制，選用的計(jì)算模型及參數(shù)往往難以完全反映客觀實(shí)際[2]，特別是難以對(duì)材料的本構(gòu)模型、結(jié)構(gòu)材料分區(qū)和多場耦合問題進(jìn)行準(zhǔn)確的考慮。隨著大壩監(jiān)測自動(dòng)化技術(shù)的應(yīng)用和數(shù)據(jù)采集時(shí)間長度的積累，監(jiān)測數(shù)據(jù)多、維數(shù)高、計(jì)算分析不直觀等問題更加突出，而管理現(xiàn)代化對(duì)分析的實(shí)時(shí)性卻提出更高的要求。隨著多種降維分析方法的不斷提出和應(yīng)用，為解決工程實(shí)際問題，特別是信息的傳輸和可視化問題，對(duì)現(xiàn)有的大壩安全高維實(shí)測數(shù)據(jù)或無窮維大壩基于降維的大壩安全特征提取不僅是必要的，而且是可行的。本文在特征正交分解(POD)、等概率粗?；?、主成分/核主成分分析(PCA/KPCA)分析的基礎(chǔ)上，將這些降維分析方法與大壩安全特征提取方法結(jié)合，提出改進(jìn)的新方法，并結(jié)合典型工程實(shí)例檢驗(yàn)新方法的有效性。

1 研究現(xiàn)狀

1.1 降維分析方法

降維分析是指通過數(shù)學(xué)、數(shù)值或統(tǒng)計(jì)方法對(duì)高維信息進(jìn)行低維變換或通過剔出冗余、誤差或無用信息，從而獲得描述結(jié)構(gòu)整體特征的低維(階)模型、特征向量或閾值的過程。常用的降維方法主要包括矩陣非負(fù)分解、自組織映射、等距映射、局部線性嵌入算法、主成分分析和核主成分分析等，本文涉及的主要降維方法包括特征正交分解(POD)、等概率粗?；⒅鞒煞?核主成分分析(PCA/KPCA)、矩陣范數(shù)(如F-范數(shù)和譜范數(shù))[3-4]。上述方法介紹的文獻(xiàn)較多[5-6]，為節(jié)約篇幅，不再贅述。

1.2 大壩安全特征提取方法

特征提取是指基于實(shí)測數(shù)據(jù)或者計(jì)算分析，提取表征大壩安全狀態(tài)的模型、參數(shù)向量或閾值的過程。特征提取方法主要包括基于實(shí)測數(shù)據(jù)的回歸分析、統(tǒng)計(jì)分析、矩陣變換和監(jiān)控(預(yù)警)指標(biāo)擬定，以及針對(duì)大壩結(jié)構(gòu)本身的剛體極限平衡法、數(shù)值仿真模擬和強(qiáng)度折減法等。

1.2.1 回歸分析

回歸分析通過擬合樣本來確定自變量和因變量之間的統(tǒng)計(jì)關(guān)系，通過回歸系數(shù)向量或方程獲得了隱含在高維樣本內(nèi)的大壩安全信息。然而，當(dāng)前常用的基于最小二乘法的逐步回歸分析提取的是一段時(shí)間序列的平均信息。實(shí)際上，大壩失事是一個(gè)典型極端情況下的小概率事件，尾部特征更加重要，因此必須提出考慮最不利工況組合及極端情況下的大壩安全特征提取方法[5-6]。

1.2.2 預(yù)警指標(biāo)擬定

監(jiān)控指標(biāo)的擬定實(shí)際上也是通過統(tǒng)計(jì)分析或數(shù)值模擬計(jì)算提取實(shí)測資料的極值分布、低維統(tǒng)計(jì)邊界[7]或屈服破壞臨界值[8-9]?；趯?shí)測資料的監(jiān)控指標(biāo)分析目前均利用單一監(jiān)測變量長時(shí)間多維向量，沒有考慮到多個(gè)變量之間信息的互補(bǔ)性，而基于數(shù)值模擬的方法基本都是單調(diào)加荷和理想彈塑性模型，與實(shí)際情況的多場耦合相差較大。

1.2.3 復(fù)雜網(wǎng)絡(luò)分析

時(shí)間序列復(fù)雜網(wǎng)絡(luò)分析是一種直觀的特征提取方法，將實(shí)測資料轉(zhuǎn)化成復(fù)雜網(wǎng)絡(luò)后，通過最短路徑、平均聚集系數(shù)、介數(shù)等復(fù)雜網(wǎng)絡(luò)有限低維參數(shù)提取時(shí)間序列隱含的特征[10]。

1.2.4 安全系數(shù)獲取

剛體極限平衡法和強(qiáng)度折減法都是獲取大壩抗滑穩(wěn)定安全特征提取的常用方法，由于極限平衡法沒有考慮材料的特性和應(yīng)力～應(yīng)變關(guān)系以及邊界條件的變化，不能反映應(yīng)力分布變化及屈服過程，因此逐漸被強(qiáng)度折減法取代。但當(dāng)前強(qiáng)度折減法主要是針對(duì)單一荷載[11]，沒考慮多場耦合。實(shí)際上大壩都在多場耦合作用下工作，而考慮多場耦合，尤其是全耦合方法計(jì)算工作量比較大，因此必須解決數(shù)據(jù)計(jì)算工作量問題。

2 基于降維分析的特征提取方法

2.1 基于POD的分位數(shù)回歸分析法

為考慮極端水位或氣溫等自變量對(duì)大壩安全的影響，深入挖掘自變量尖峰厚尾等非正態(tài)分布以及不同分位點(diǎn)所包含的安全信息，決定選用分位數(shù)回歸模型對(duì)實(shí)測數(shù)據(jù)進(jìn)行分析，分位數(shù)回歸具有如下優(yōu)點(diǎn)：①對(duì)模型中的隨機(jī)擾動(dòng)項(xiàng)無需做任何分布的假定，整個(gè)回歸模型就具有很強(qiáng)的穩(wěn)健性；②本身沒有使用一個(gè)連接函數(shù)來描述因變量的均值和方差的相互關(guān)系，因此具有比較好的彈性性質(zhì)；③由于對(duì)所有分位數(shù)進(jìn)行回歸，因此對(duì)于數(shù)據(jù)中出現(xiàn)的異常點(diǎn)具有耐抗性；④不同于普通的最小二乘回歸，分位數(shù)回歸對(duì)于因變量具有單調(diào)變換性；⑤分位數(shù)回歸估計(jì)出來的參數(shù)具有在大樣本理論下的漸進(jìn)優(yōu)良性。

為得到回歸系數(shù)，均值回歸可采用最小平方法，中位數(shù)回歸可采用最小一乘法(最小絕對(duì)偏差法)，對(duì)于分位數(shù)回歸可采用線性規(guī)劃法等。為進(jìn)一步減少自變量的個(gè)數(shù)和克服自變量之間的多重共線性，將POD的小樣本逼近、特征提取與分位數(shù)回歸的優(yōu)勢(shì)結(jié)合起來，首先對(duì)實(shí)測數(shù)據(jù)求取特征正交基，再利用特征正交基進(jìn)行分位數(shù)回歸分析，從而獲得不同分位點(diǎn)上各綜合自變量對(duì)因變量的影響程度。

2.2 基于實(shí)測資料的復(fù)雜網(wǎng)絡(luò)分析法

該方法可以有效克服離散誤差對(duì)數(shù)據(jù)分析結(jié)論的影響，實(shí)現(xiàn)數(shù)據(jù)分析的可視化，具體步驟如下。

a) 選取典型和重要的大壩測點(diǎn)及其測值。

b) 運(yùn)用等概率粗粒化的方法，將測值時(shí)間序列轉(zhuǎn)化成符號(hào)序列[10]。①設(shè)xmax和xmin分別是測值時(shí)間序列的最大值和最小值，轉(zhuǎn)化后的符號(hào)共有N種，記為s1，s2，…，sN，定義:

Si=Sj,xmin+(j-1)d≤xi

(1)

②記大壩測值時(shí)間序列為{x(l),l=1,…N}，通過式(2)計(jì)算測值時(shí)間序列波動(dòng)k(l)：

(2)

式中 Δl——時(shí)間間隔；

③通過式(3)計(jì)算不同波動(dòng)值可能出現(xiàn)的概率p(k):

(3)

式中Num(x)——測值序列的波動(dòng)模態(tài)x發(fā)生的次數(shù)；

④把測值時(shí)間序列的波動(dòng)k(l)分為5個(gè)區(qū)間，定義5個(gè)特征字符為式：

(4)

式中t——測值數(shù)值快速增加；r——測值數(shù)值緩慢增加；e——測值數(shù)值不變；d——測值數(shù)值緩慢降低；f——測值數(shù)值快速降低。從而將測值時(shí)間序列轉(zhuǎn)換成符號(hào)序列為：

S=(S1,S2,S3,…)，Si∈(t，r，e，d，f)

(5)

c) 將符號(hào)序列編碼通過excel2pajek軟件轉(zhuǎn)換成pajek軟件可處理的格式構(gòu)造復(fù)雜網(wǎng)路。

d) 分析復(fù)雜網(wǎng)絡(luò)的特征參數(shù)，獲取該大壩監(jiān)測項(xiàng)目的敏感測點(diǎn)及其性態(tài)演化規(guī)律，判別大壩的工程安全性態(tài)。

2.3 基于PCA/KPCA的預(yù)警指標(biāo)提取法

置信區(qū)間法、典型監(jiān)測效應(yīng)量小概率法、極限狀態(tài)法和極值理論等都是針對(duì)單個(gè)測點(diǎn)提出的，沒有考慮多個(gè)測點(diǎn)之間的信息冗余或互補(bǔ)，更沒有考慮土石體在失穩(wěn)破壞中的應(yīng)變局部化等效應(yīng)。本文將多個(gè)監(jiān)測項(xiàng)目和測點(diǎn)數(shù)據(jù)按時(shí)間次序形成數(shù)據(jù)帶，選擇一定時(shí)間長度的窗口，根據(jù)主成分和核主成分算法，在85%貢獻(xiàn)率的基礎(chǔ)上得到PCA/KPCA主成分個(gè)數(shù)和對(duì)應(yīng)的矩陣范數(shù)2個(gè)指標(biāo)作為預(yù)警指標(biāo)。

核函數(shù)常用選取方法有Cross-Validation法、多次實(shí)驗(yàn)法和混合核函數(shù)方法，通過對(duì)比分析，本文選擇高斯徑向基核函數(shù)進(jìn)行核主成分分析[5]。

2.4 基于POD的多場耦合的強(qiáng)度折減法

由于溫度荷載對(duì)大壩抗滑穩(wěn)定安全系數(shù)影響不大，暫不考慮溫度荷載，考慮飽和-非飽和流固耦合方程：

(6)

式中uw,j——孔隙壓力；Sw——飽和度；εv——土骨架體應(yīng)變；kwij——滲透系數(shù)；ρw——水的密度；bj——單位體積力；n——孔隙率。壩體及壩基采用非關(guān)聯(lián)準(zhǔn)則的彈塑性本構(gòu)模型、Mohr-Coulomb屈服準(zhǔn)則：

(7)

(8)

基于POD的強(qiáng)度折減法首先利用大壩安全監(jiān)測資料構(gòu)造瞬像集合，采用POD方法提取特征正交基，結(jié)合Galerkin投影法建立有限元方程[12-13]，再結(jié)合強(qiáng)度折減法即可獲得大壩抗滑穩(wěn)定安全系數(shù)，從而可以大幅度節(jié)約計(jì)算時(shí)間，提高計(jì)算效率。

3 工程實(shí)例

本文選擇壩型分別為常態(tài)混凝土重力拱壩、均質(zhì)土壩和深厚覆蓋層上碾壓混凝土重力壩3個(gè)典型工程實(shí)例檢驗(yàn)上述方法的有效性。

3.1 安徽某重力拱壩

3.1.1 工程存在問題

3.1.2 研究內(nèi)容和結(jié)果

3.1.2.1 基于POD的分位數(shù)回歸分析

a) 擬合精度分析：通過對(duì)發(fā)現(xiàn)裂縫過程中間時(shí)段拱冠梁和左右1/4拱壩頂測點(diǎn)徑向分析，基于POD的分位數(shù)回歸分析方法可以很好地適應(yīng)小樣本性質(zhì)和非平穩(wěn)性，擬合精度較高，參見圖1—3。

b) 大壩安全特征變化分析：將發(fā)現(xiàn)裂縫的時(shí)段分成前后2個(gè)時(shí)段，分別進(jìn)行回歸分析，通過模型參數(shù)前后變化分析裂縫對(duì)大壩安全性態(tài)的整體影響。前一時(shí)段典型測點(diǎn)回歸方程為：

(9)

后一時(shí)段回歸方程為：

(10)

通過上述系數(shù)矩陣可見，采用相同的基Ψ1、Ψ2，裂縫存在前后2個(gè)時(shí)段大壩典型測點(diǎn)的徑向變形并未受到明顯影響，說明大壩整體安全性態(tài)穩(wěn)定，這一結(jié)論與大壩仿真計(jì)算和專家綜合評(píng)價(jià)結(jié)果一致。

3.1.2.2 基于實(shí)測資料的復(fù)雜網(wǎng)絡(luò)分析

同樣利用8、18、26號(hào)壩塊壩頂徑向?qū)崪y位移數(shù)據(jù)，分成前后兩段，分別采用復(fù)雜網(wǎng)絡(luò)進(jìn)行分析。

通過圖4可知，裂縫存在前后2個(gè)時(shí)段大壩典型測點(diǎn)的復(fù)雜網(wǎng)絡(luò)未出現(xiàn)顯著變化，說明大壩整體安全性態(tài)穩(wěn)定，這一結(jié)論與上文分析方法得到的結(jié)論一致[14]。

3.2 某均質(zhì)試驗(yàn)土壩

3.2.1 研究目的

為研究潰壩過程中大壩安全性態(tài)的變化，建造了一座水庫并進(jìn)行了潰壩試驗(yàn)。該試驗(yàn)水庫最大庫容10萬m3，大壩為黏性均質(zhì)土壩，壩高9.7 m，壩頂長120 m、寬3 m。壩體分層碾壓填筑。試驗(yàn)過程包括大壩填筑(儀器埋設(shè))、首次蓄水、水位穩(wěn)定、第2次蓄水和潰決等過程，共歷時(shí)46 d[14]。

3.2.2 研究內(nèi)容和結(jié)果

3.2.2.1 基于POD的分位數(shù)回歸分析

假設(shè)實(shí)測數(shù)據(jù)分析前，并不知曉潰壩開始時(shí)間，為了檢驗(yàn)本方法正確性，將得到的變形數(shù)據(jù)分成前80組進(jìn)行基于POD的回歸分析。前80組位移進(jìn)行模擬(圖5)。

由圖5可知，模型的擬合值與實(shí)測值擬合較好，表明所建模型能夠較好地反映前80組位移監(jiān)測值的變化。為了檢驗(yàn)?zāi)Ｐ偷念A(yù)測能力，預(yù)測了后15 d的位移值，見圖6。

由圖6可知，前29組預(yù)測值與實(shí)測值相差較小，說明所建立的模型精度較高，能夠較好的預(yù)測后期數(shù)據(jù)。而前29組數(shù)據(jù)之后預(yù)測值與實(shí)測值相差較大，表明此時(shí)大壩開始處于潰壩前期，瀕臨潰壩。第59組數(shù)據(jù)的模型預(yù)測值突然跳動(dòng)很大，此時(shí)已經(jīng)潰壩。

將上述第59組數(shù)據(jù)開始將數(shù)據(jù)分成兩段，得到位移與POD正交基的回歸方程分別為：

y=-6.841Ψ1-9.151Ψ2-168.16

(11)

y=-3.286Ψ1-5.473Ψ2-288.52

(12)

式(11)與式(12)存在明顯差異，說明通過降維分析獲得的方程系數(shù)反映了土壩安全性態(tài)的變化。

3.2.2.2 基于實(shí)測資料的復(fù)雜網(wǎng)絡(luò)分析

根據(jù)上文時(shí)間結(jié)點(diǎn)，將土壩滲流實(shí)測資料分成潰壩前后兩段，對(duì)應(yīng)的復(fù)雜網(wǎng)絡(luò)見圖7、8。

通過圖7、8不難發(fā)現(xiàn)，兩者之間存在明顯差別，也就是說基于復(fù)雜網(wǎng)絡(luò)的降維分析方法能發(fā)現(xiàn)的潰壩時(shí)間與實(shí)際一致。

3.2.2.3 潰壩預(yù)警指標(biāo)擬定

土壩潰壩實(shí)驗(yàn)過程見文獻(xiàn)[15]，通過在潰壩過程中埋設(shè)安裝在壩體上的63臺(tái)儀器進(jìn)行主成分分析。將63個(gè)測點(diǎn)每連續(xù)兩天(每天48次)的測值構(gòu)成一個(gè)窗口大小為96×63的數(shù)據(jù)窗口，其中63為測點(diǎn)數(shù)，對(duì)應(yīng)列向量個(gè)數(shù)；96為測次，對(duì)應(yīng)行向量個(gè)數(shù)。從2012年10月16日—2012年11月17日共32個(gè)窗口。

應(yīng)用PCA對(duì)監(jiān)測數(shù)據(jù)進(jìn)行主成分分析，結(jié)果見圖9，其中橫坐標(biāo)為時(shí)間窗口，縱坐標(biāo)為主成分累計(jì)貢獻(xiàn)率(下同)。通過圖9可知，10月16—31日滿足一定累計(jì)貢獻(xiàn)率的主成分個(gè)數(shù)一直增加，10月31日至11月14日逐漸下降，在11月17日又突然增加。同時(shí)第一主成分貢獻(xiàn)率呈現(xiàn)先緩慢減少、再緩慢增加又突然減少的變化過程。可見主成分個(gè)數(shù)反映了大壩安全性態(tài)變化，隨著大壩安全性態(tài)的變化，主成分個(gè)數(shù)先增加、后減少，而在潰壩時(shí)突然增加。根據(jù)分形理論，在潰壩前主成分?jǐn)?shù)量減少對(duì)應(yīng)了大壩整體降維和應(yīng)變局部化過程，而17日主成分個(gè)數(shù)突然增加則反映了大壩結(jié)構(gòu)的離散化和增維過程[15]。

一定累計(jì)貢獻(xiàn)率下，96×63窗口核主成分個(gè)數(shù)隨時(shí)間變化過程線見圖10。圖10與圖9結(jié)果基本一致，10月16—30日核主成分個(gè)數(shù)一直增加，11月3—14日逐漸下降，在11月17日又突然增加。同時(shí)第一主成分貢獻(xiàn)率呈現(xiàn)先緩慢減少、再增加又突然減少的趨勢(shì)。即隨著大壩安全性態(tài)的變化，主成分個(gè)數(shù)先增加、后減少，而在潰壩時(shí)突然增加。

由圖9、10可知，總體上10月17—31日主成分個(gè)數(shù)一直增加，10月31—11日逐漸下降，在11月17日又突然增加。同時(shí)主成分范數(shù)也呈現(xiàn)先緩慢減少、再增加又突然減少的過程，上述過程和潰壩過程吻合。在統(tǒng)計(jì)主成分個(gè)數(shù)的同時(shí)還計(jì)算相應(yīng)的矩陣。通過圖11可知，盡管采用不同的范數(shù)計(jì)算，但計(jì)算結(jié)果基本一致，即范數(shù)的變化對(duì)具體范數(shù)計(jì)算方法不敏感，因此主成分范數(shù)可以作為潰壩分級(jí)預(yù)警指標(biāo)。

3.3 山西某碾壓混凝土重力壩

3.3.1 工程存在的問題

大壩為深厚覆蓋層上碾壓混凝土壩，由于在大壩主體工程竣工時(shí)壩體壩肩固結(jié)和帷幕灌漿未完成，導(dǎo)致廊道內(nèi)漏水嚴(yán)重。2014年開始進(jìn)行以壩肩帷幕灌漿、壩基固結(jié)灌漿和壩后連續(xù)墻施工等應(yīng)急除險(xiǎn)加固措施。為了解除險(xiǎn)加固后大壩的安全性態(tài)，采用數(shù)值方法進(jìn)行分析研究。

3.3.2 研究內(nèi)容和結(jié)果

大壩埋設(shè)安裝了大量的位移、滲流、溫度和應(yīng)變傳感器，為基于POD的多場耦合方法提供了條件。分別用除險(xiǎn)加固前后的多場實(shí)測數(shù)據(jù)建立瞬像集合，采用基于特征正交分解方法提取正交基，并將正交基作為試函數(shù)，利用有限元法獲得大壩除險(xiǎn)加固后的變形、滲流分布和抗滑穩(wěn)定安全系數(shù)。表1給出了不同計(jì)算工況下采用POD的多場耦合方法與直接有限元法的計(jì)算機(jī)機(jī)時(shí)比統(tǒng)計(jì)。

表1 不同工況計(jì)算時(shí)間花費(fèi)比統(tǒng)計(jì)

3.3.3 研究結(jié)論

a) 通過采用基于POD的多場耦合分析方法，實(shí)際6個(gè)工況計(jì)算結(jié)果與未采用POD的基于COMSOL的計(jì)算結(jié)果一致，說明基于POD的多場耦合方法是可行的。

b) 采用POD的多場耦合分析方法可以大幅度降低計(jì)算工作量，使得計(jì)算時(shí)間縮小為原來的1/80～1/116，對(duì)計(jì)算結(jié)果幾乎無影響。

c) 除險(xiǎn)加固后整體安全系數(shù)滿足要求，整體應(yīng)力分布基本合理，但也有少部分壩體區(qū)域存在拉應(yīng)力，后期應(yīng)加強(qiáng)監(jiān)測。

4 結(jié)語

基于大壩安全監(jiān)測數(shù)據(jù)結(jié)合相關(guān)理論分析和數(shù)值模擬，分別采用特征正交分解、分位數(shù)回歸、主成分/核主成分等多種方法，對(duì)大壩安全信息降維的理論、方法和應(yīng)用進(jìn)行了研究，主要結(jié)論如下。

a) 基于特征正交分解的分位數(shù)回歸建模分析方法能解決大壩實(shí)測數(shù)據(jù)分析中的小樣本、多重共線性和變量尖峰厚尾分布等問題。

b) 采用滑動(dòng)窗口建立大壩實(shí)測數(shù)據(jù)的復(fù)雜網(wǎng)絡(luò)模型，可以通過復(fù)雜網(wǎng)絡(luò)參數(shù)分析直觀顯示大壩安全性態(tài)變化。

c) 采用主成分/核主成分分析方法，對(duì)潰壩試驗(yàn)多監(jiān)測項(xiàng)目、多測點(diǎn)數(shù)據(jù)進(jìn)行整體降維分析，經(jīng)模型試驗(yàn)驗(yàn)證，主成分個(gè)數(shù)和特征向量矩陣范數(shù)可作為均質(zhì)土壩潰壩預(yù)警指標(biāo)。

d) 基于POD的多場耦合強(qiáng)度折減法計(jì)算結(jié)果與COMSOL Multiphysics軟件計(jì)算結(jié)果一致，但新方法顯著節(jié)省了計(jì)算機(jī)時(shí)。