基于敏感性分析與粒子群算法的拱壩原型動彈性模量反演方法

李火坤，王 剛，魏博文，黃 偉，陳良捷

（南昌大學 建筑工程學院，江西 南昌 330031）

1 研究背景

我國的高拱壩多位于西南高山峽谷地區，隨著一批新建高拱壩相繼投入運行，復雜的地質條件及自然環境給長期服役的高拱壩帶來不同程度安全隱患[1-2]。高拱壩的物理力學參數對精準評估大壩結構整體運行性態至關重要，目前獲取大壩物理力學參數的主要方法有現場取樣測試法和基于監測資料的反演法[3-4]?，F場取樣測試法工作量大，損壞原結構，且只能了解大壩局部結構參數特性，不經濟也不全面。而原型監測資料蘊含了大壩真實的工作性態，因此利用大壩監測資料反演大壩及地基物理力學參數是一種有效的方法。

利用監測資料進行大壩物理力學參數反演的方法主要分為基于靜態監測資料的反演方法和基于振動測試的參數反演方法兩種。其中，基于靜態監測資料（如靜力位移）的大壩物理力學參數反演在近幾十年來一直是大壩安全監控研究領域的熱點，并取得了豐碩的研究成果[5-10]，為保障我國大壩的安全運行做出了重要貢獻。近年來，隨著結構動力參數辨識技術的發展，基于振動測試的結構物理力學參數反演研究發展迅速[11]，由于振動實測響應所包含的振動模態信息更全面地反映了結構整體力學行為特征，因此根據振動模態參數進行大壩物理力學參數反演所得到的結果更符合實際。馮新等[12]探討了一種基于不完全模態觀測數據的混凝土壩分區彈性模量反演方法，建立了基于不完全模態數據的大壩彈性模量優化反演模型；康飛等[13]基于不完全模態測試數據提出了一種混合單純形人工蜂群算法進行動力材料參數識別方法；程琳等[14-15]根據大壩強震監測數據識別了某混凝土壩的振動模態參數，并建立了反映大壩動彈性模量和各階模態參數之間非線性關系的多輸出支持向量機代理模型（M-SVM），結合遺傳算法對大壩彈性模量進行了反演；王登剛等[16]考慮模態觀測數據的不確定性，引入不確定問題求解的區間分析思想，建立了利用先驗約束條件的混凝土重力壩動態參數的區間反演模型，采用約束變尺度法求解獲得了壩體混凝土和基巖動彈性模量參數的區間范圍。文獻[17]提出了基于遺傳算法的拱壩材料參數反演方法，但結果驗證中仍存在少部分測點精度不高等問題，仍需進一步改進。

本文以二灘拱壩為研究對象，提出一種基于敏感性分析與粒子群算法的拱壩原型動彈性模量反演方法，建立反映拱壩及地基各分區動彈性模量與拱壩模態參數非線性映射關系的三階多項式響應面數學模型，基于正交試驗法分析拱壩及地基的動彈性模量區域對拱壩模態參數的敏感性，確定待反演的動彈性模量區域；基于拱壩原型振動模態參數，提出基于響應面數學模型計算模態參數和原型振動模態參數的拱壩及地基分區動彈性模量反演最優化數學模型，建立參數反演的目標函數并采用自適應慣性權重的粒子群優化（AWPSO）算法進行尋優求解，最終反演出拱壩及地基各分區實際動彈性模量。

2 拱壩及地基分區動彈性模量反演方法

2.1 響應面數學模型響應面法是以數學理論和統計方法為基礎，通過近似構造一個具有顯式的函數表達式來擬合變量之間的非線性數學關系。本文采用響應面法來表述拱壩及地基分區動彈性模量與拱壩模態參數之間的非線性映射關系，以拱壩固有頻率和振型作為模態參數變量，在準確表達拱壩動彈性模量和模態參數之間的非線性函數關系并充分考慮計算效率的前提下，采用三階多項式響應面方程來表達動彈性模量E與頻率f、振型φ之間的關系，建立的響應面模型如下：

本文采用拉丁超立方抽樣法生成動彈性模量樣本，通過建立拱壩三維有限元模型并計算不同動彈性模量樣本組合下的拱壩模態參數，根據式（1）和式（2），采用多元回歸方法擬合動彈性模量與頻率、振型之間的響應面模型，并以此替代有限元模型。

為確保拉丁超立方抽樣法抽取的動彈性模量樣本均勻覆蓋取樣區間，在動彈性模量的取值范圍內均勻生成U個動彈性模量樣本點時，將歸一至并分成U 等份（即抽樣區間為在每個抽樣區間內以不等且隨機的抽樣概率生成樣本，如圖1所示。

圖1 拉丁超立方抽樣方法

在擬合響應面方程時，響應面的精度是評價響應面模型是否合理的重要指標，本文以有限元模型模態參數（即固有頻率fp和節點歸一化振型φpq）計算值與響應面模型模態參數輸出值之間的相對誤差來衡量響應面數學模型的精度，其表達式如下：

式中： Re為響應面模型的精度，本文響應面模型精度控制5‰以內；yRS為響應面模型模態參數的輸出值；yFEM為有限元模型模態參數的計算值。

2.2 敏感性分析式（1）和式（2）中動彈性模量的個數m 直接影響到動彈性模量反演的計算效率和精度，為合理選取對拱壩模態參數敏感性高的待反演動彈性模量區域，本文采用正交試驗法[18]分析不同動彈性模量區域對拱壩模態參數的敏感性，以拱壩第一階固有頻率作為評價指標評價各動彈性模量區域對拱壩模態參數的敏感性，采用方差分析法對評價指標進行分析。

首先，計算所有動彈性模量區域組合下拱壩第一階固有頻率總離差的平方和Sr、各動彈性模量區域對拱壩第一階固有頻率離差的平方和SA及試驗相對誤差離差的平方和SE，其表達式如下：

其次，計算試驗總自由度ZT=K-1、各動彈性模量區域數自由度ZA=u-1，并確定試驗相對誤差自由度以此計算第A個動彈性模量區域的方差試驗相對誤差的方差

最后，構造動彈性模量第A個區域對拱壩第一階固有頻率敏感性的統計量TA：

由上述過程可判別拱壩及地基不同動彈性模量區域對拱壩模態參數的敏感性。根據文獻[19]關于方差顯著性分析的研究成果，本文選取顯著性水平λ分別為0.025、0.05、0.1 來確定不同動彈性模量區域對拱壩第一階固有頻率的影響并將其劃分為4個等級，即：（1）視為影響高度顯著；（2）視為影響顯著；（3）視為有一定影響；（4）視為無影響。由此，可確定拱壩及地基待反演的動彈性模量區域。

2.3 最優化目標函數以拱壩原型振動模態參數與響應面數學模型計算模態參數的相對偏差作為目標函數J，將拱壩固有頻率和各測試節點歸一化振型的響應面數學模型計算值與拱壩原型識別值相比較，從而將拱壩及地基分區動彈性模量的反演問題轉換為目標函數的最優化求解問題，構造的目標函數如下：

本文采用基于經驗模態分解法（EMD）和奇異熵定階的隨機子空間法開展原型拱壩振動模態參數識別，具體過程可參考文獻[20-21]。

2.4 自適應慣性權重的粒子群算法采用粒子群優化（PSO）算法對目標函數J進行最優化求解。為實現最優解的收斂和避免陷入局部最優，利用自適應慣性權重對PSO 算法進行優化，使其在搜索過程中能夠依據粒子的位置動態變化，從而搜尋到最優解，其基本過程如下。

假定待反演的動彈性模量區域個數為m，由n個粒子組成一個初始群落，則在n×m 空間中，其初始粒子群E為：

式中： c1、 c2為學習因子；r1、 r2為[0，1]范圍內的均勻隨機數；分別為粒子i 第k 次迭代中第d 維的速度與位置；為粒子i 第d 維的個體適應值最小處位置；為粒子群體第d 維的種群適應值最小處位置；ω為慣性權重；ωmin、 ωmax分別為慣性權重的最小值和最大值Jt表示粒子實時的適應值；Javg、 Jmin分別為當前所有粒子適應值函數的平均值和最小值。

當達到迭代次數或者停止條件時，則停止搜索并輸出最優值，否則將繼續迭代并更新粒子，直至搜索到滿足收斂準則的值為止。

2.5 反演流程拱壩及地基分區動彈性模量反演流程如圖2所示?；静襟E說明如下：（1）步驟1?；贓MD對拱壩原型振動響應信號進行降噪處理，計算降噪后信號奇異熵增量并確定拱壩振動模態階次，采用隨機子空間法識別拱壩原型振動模態參數；（2）步驟2。根據拱壩及地基彈性模量分區設計情況，建立拱壩流-固耦合三維有限元模型，計算拱壩及地基不同動彈性模量分區情況下的拱壩模態參數；采用正交試驗法確定拱壩模態參數敏感性高的待反演動彈性模量區域；（3）步驟3。采用拉丁超立方抽樣法生成動彈性模量樣本集，以此作為有限元模型參數輸入，計算各樣本集下所對應的拱壩模態參數（頻率和振型）；（4）步驟4。構建反映拱壩模態參數與各分區動彈性模量之間非線性映射關系的響應面數學模型，以步驟3的計算結果為依據，采用多元回歸法確定響應面數學模型；（5）步驟5?；诠皦卧驼駝幽B參數與響應面數學模型計算模態參數的相對偏差最小，構建拱壩及地基分區動彈性模量反演的目標函數；利用AWPSO 算法進行優化求解，最終反演出拱壩及地基實際動彈性模量。

圖2 拱壩及地基分區動彈性模量反演流程

3 二灘拱壩原型工程實例

3.1 原型振動響應測試情況本文以二灘拱壩為研究對象，開展了泄流激勵下的拱壩原型振動測試，將徑向動位移傳感器布置于壩頂及拱冠梁處；其中，壩頂共布置7個測點，編號為B1 至B7，如圖3所示；拱冠梁共布置4個測點，編號為B8 至B11，如圖4所示。傳感器為DP 型低頻動位移傳感器，其參數見文獻[11]。測試工況如表1所示，采樣頻率為200 Hz，拱壩典型測點（B5）原型振動響應時程線如圖5所示。

表1 原型振動測試工況

圖3 壩頂測點布置

圖4 拱冠梁測點布置

3.2 拱壩振動模態參數辨識基于EMD法對拱壩原型振動響應信號進行降噪處理，并引入信噪比（SNR）作為評價降噪效果的量化指標，計算結果如表2所示。從降噪效果來看，降噪后的信號具有更高的信噪比，典型測點降噪前后時程線（局部放大）如圖6所示。以降噪后的信號為基礎，采用奇異熵定階的隨機子空間法對拱壩振動模態參數進行辨識，識別結果如表3所示。以工況1為例，當信號階次達到6時，奇異熵增量隨階次變化出現穩定（如圖7所示），剔除特征值中非模態項和共軛項[22]，得到拱壩實際振動模態階次為3階；根據所識別的拱壩固有頻率穩定圖（如圖8所示），可得該工況下拱壩前三階固有頻率分別為1.39、1.51和2.13 Hz。

圖5 工況1 下B5 測點原型振動響應時程

圖6 工況1 下降噪前后B5 測點時程（局部放大）

表2 降噪前后信號指標對比 （單位：dB）

表3 拱壩振動模態參數識別結果

3.3 拱壩計算模型本文中響應面數學模型的確定是以拱壩有限元模型為基礎。根據二灘拱壩混凝土的分區設計及大壩地質資料[23]，二灘拱壩及地基彈性模量的分區分別由3個壩體區域和5個地基區域組成，其中壩體區域編號分別為A、B、C，地基區域編號分別為D、E、F、G、H，其設計值如表4所示。本文運用ANSYS 有限元軟件建立拱壩-地基-水體有限元模型，如圖9所示，采用Fluid30 水體單元來模擬流固耦合效應，地基單元按照無質量彈性地基進行模擬，地基模擬范圍為：深度取200 m，上游取424 m，下游取534 m，左、右壩肩取100 m；地基四周采用法向約束，底部采用固定約束。拱壩有限元模型共劃分有282 049個單元和264 720個節點。

圖7 工況1 下奇異熵增量隨階次變化

圖8 工況1 下固有頻率

圖9 拱壩-地基-水體有限元模型

表4 拱壩及地基彈性模量區域設計值

3.4 動彈性模量區域的敏感性分析以二灘拱壩及地基的8個彈性模量分區作為正交試驗因素，根據各區域彈性模量設計值及文獻[24-25]，確定拱壩及地基各動彈性模量區域的敏感性分析取值范圍，并以其上、下限和設計值作為3個試驗水平，如表5所示。根據正交試驗因素和試驗水平，形成了27 組動彈性模量區域敏感性分析的試驗組合并進行有限元模態計算，得到各試驗組合下的拱壩第一階固有頻率。

表5 動彈性模量試驗水平 （單位：GPa）

表6 敏感性分析結果

圖10 動彈性模量區域顯著性分析

基于方差分析法分析各試驗組合下的拱壩第一階固有頻率，得出二灘拱壩及地基各動彈性模量區域對拱壩模態參數的敏感性分析結果，如表6所示。按統計量由大到小進行排序，可得到對拱壩模態參數敏感性由高到低的動彈性模量區域依次為A→B→G→C→E→D→F→H。根據統計量分布表[19]，得到T1-0.025=5.4564，T1-0.005=4.1028，T1-0.1=2.9245，如圖10所示，可見對模態參數影響高度顯著的區域為A、B、C、E、G，有一定影響的區域為D，無影響的區域為F和H。由此，本文確定待反演的動彈性模量區域為A、B、C、D、E和G 等6個區域。

3.5 響應面數學模型的構建根據所確定的待反演動彈性模量區域，采用拉丁超立方抽樣法隨機生成500 組動彈性模量樣本集，參照各動彈性模量區域的敏感性分析取值范圍，擬定A、B、C、D、E和G 共6個區域動彈性模量樣本抽樣區間分別為[28.34，42.5]、[27.37，41.05]、[26.4，39.6]、[28，42]、[20，30]、[16，24]（單位：GPa）。將500 組動彈性模量樣本集輸入拱壩有限元模型，計算動彈性模量樣本集下所對應的500 組拱壩前三階固有頻率和振型，根據式（1）和式（2）擬合動彈性模量E與拱壩固有頻率 fp及各測試節點歸一化振型值φpq之間的響應面數學模型，并通過回歸分析計算得到數學模型的待定系數和從而確定響應面數學模型，并以此替代有限元計算模型。

圖11 頻率響應面精度

圖12 B5 測點歸一化振型系數響應面精度

根據式（3）計算得到各階頻率和典型測點（B5）歸一化振型系數響應面精度均在5‰之內，如圖11、圖12所示。

圖13 粒子群優化算法對比

3.6 優化求解采用AWPSO 算法對二灘拱壩及地基動彈性模量反演的最優化數學模型進行迭代求解，結果如表7所示。同時與標準的PSO 算法適應度值的計算結果進行對比，如圖13所示。從圖13可知，兩種算法都在第12 次迭代時收斂，但采用AWPSO 算法所得的適應值更小，具有更高的精度。

表7 動彈性模量反演結果

3.7 反演結果驗證為驗證反演結果的準確性，以工況2 作為驗證工況，將動彈性模量反演結果輸入拱壩有限元模型進行模態計算，并將模態計算結果與該工況下拱壩原型振動模態參數（前三階固有頻率和各測點歸一化振型系數）識別結果相對比，如表8、表9所示。從對比結果來看，對于固有頻率而言，計算結果與模態識別結果在數值上比較吻合，第一階固有頻率相對誤差為0.29%，第二階固有頻率相對誤差為0.76%，第三階固有頻率相對誤差為2.05%。對于拱壩測點歸一化振型系數而言，大部分測點識別值和計算值之間的相對誤差都在10%以內；少部分測點的相對誤差稍大，如第一階B2 測點的相對誤差為10.75%，第二階B6 測點的相對誤差為-10.38%，主要是由于位于壩肩的測點振幅較小，受噪聲干擾大，導致振型識別存在一定誤差?？傮w上，基于動彈性模量反演結果的拱壩計算模態參數與原型振動模態參數吻合較好，該方法的反演結果是合理可靠的。

表8 拱壩前三階固有頻率對比

表9 拱壩測點歸一化振型系數對比

4 結論

本文提出了基于敏感性分析與粒子群算法的拱壩原型動彈性模量反演方法，主要結論如下：（1）構建了描述拱壩及地基分區動彈性模量與拱壩模態參數之間非線性映射關系的響應面數學模型，基于正交試驗法深度挖掘了各動彈性模量區域對拱壩模態參數的敏感性，提高了動彈性模量反演的計算效率。（2）提出了基于響應面數學模型計算模態參數和拱壩原型振動模態參數之間相對誤差的目標函數，建立了拱壩及地基分區動彈性模量反演的最優化數學模型，將動彈性模量反演問題轉換成目標函數最優化求解問題。（3）利用AWPSO 算法對最優化數學模型進行尋優，克服了PSO 搜索精度不高的不足，具有更好的適用性。工程實例驗證結果表明，該方法是合理可靠的，可為拱壩動彈性模量反演提供一條新思路。