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

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

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

1 研究背景

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

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

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

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

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

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

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

圖1 拉丁超立方抽樣方法

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

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

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

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

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

最后，構(gòu)造動彈性模量第A個區(qū)域?qū)皦蔚谝浑A固有頻率敏感性的統(tǒng)計量TA：

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

2.3 最優(yōu)化目標函數(shù)以拱壩原型振動模態(tài)參數(shù)與響應面數(shù)學模型計算模態(tài)參數(shù)的相對偏差作為目標函數(shù)J，將拱壩固有頻率和各測試節(jié)點歸一化振型的響應面數(shù)學模型計算值與拱壩原型識別值相比較，從而將拱壩及地基分區(qū)動彈性模量的反演問題轉(zhuǎn)換為目標函數(shù)的最優(yōu)化求解問題，構(gòu)造的目標函數(shù)如下：

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

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

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

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

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

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

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

3 二灘拱壩原型工程實例

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

表1 原型振動測試工況

圖3 壩頂測點布置

圖4 拱冠梁測點布置

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

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

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

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

表3 拱壩振動模態(tài)參數(shù)識別結(jié)果

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

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

圖8 工況1 下固有頻率

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

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

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

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

表6 敏感性分析結(jié)果

圖10 動彈性模量區(qū)域顯著性分析

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

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

圖11 頻率響應面精度

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

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

圖13 粒子群優(yōu)化算法對比

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

表7 動彈性模量反演結(jié)果

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

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

表9 拱壩測點歸一化振型系數(shù)對比

4 結(jié)論

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