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

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

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

1 研究現(xiàn)狀

1.1 降維分析方法

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

1.2 大壩安全特征提取方法

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

1.2.1 回歸分析

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

1.2.2 預警指標擬定

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

1.2.3 復雜網(wǎng)絡分析

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

1.2.4 安全系數(shù)獲取

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

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

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

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

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

2.2 基于實測資料的復雜網(wǎng)絡分析法

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

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

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

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

(1)

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

(2)

式中 Δl——時間間隔；

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

(3)

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

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

(4)

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

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

(5)

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

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

2.3 基于PCA/KPCA的預警指標提取法

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

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

2.4 基于POD的多場耦合的強度折減法

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

(6)

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

(7)

(8)

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

3 工程實例

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

3.1 安徽某重力拱壩

3.1.1 工程存在問題

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

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

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

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

(9)

后一時段回歸方程為：

(10)

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

3.1.2.2 基于實測資料的復雜網(wǎng)絡分析

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

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

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

3.2.1 研究目的

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

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

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

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

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

由圖6可知，前29組預測值與實測值相差較小，說明所建立的模型精度較高，能夠較好的預測后期數(shù)據(jù)。而前29組數(shù)據(jù)之后預測值與實測值相差較大，表明此時大壩開始處于潰壩前期，瀕臨潰壩。第59組數(shù)據(jù)的模型預測值突然跳動很大，此時已經(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 基于實測資料的復雜網(wǎng)絡分析

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

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

3.2.2.3 潰壩預警指標擬定

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

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

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

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

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

3.3.1 工程存在的問題

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

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

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

表1 不同工況計算時間花費比統(tǒng)計

3.3.3 研究結(jié)論

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

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

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

4 結(jié)語

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

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

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

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

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