基于GRU-WPA法的重力壩地基力學參數反演分析方法

徐喆成，劉曉青，季 威，王雪紅

(1.河海大學水利水電工程學院，江蘇 南京 210024；2.南京市長江河道管理處，江蘇 南京 210011)

0 引 言

大壩在建設期和運行期的監測數據是評價大壩安全性能的重要指標。利用監測數據對大壩及其地基的力學參數進行反演對大壩安全運行具有重要意義[1]。反演分析的目標函數往往是一個非線性多峰函數，監測數據與力學參數之間的顯式關系是難以識別的，常采用正演優化反分析方式實現參數識別過程，而傳統優化算法易產生局部極值問題[2]。對于許多需要在龐大的搜索空間中搜尋最優解的多極值實際工程問題，合理設計的智能優化算法可以實現搜尋全局最優。同時，地基力學參數反演分析過程中需要調用的正演分析模型往往是有限元模型，計算耗時較長，效率較低，在實際工程中的應用較為局限，而基于機器學習算法的代理模型來代替有限元模型，可以極大地提高計算效率，近年來在反演分析領域得到了越來越多的應用。代理模型的基本原理是通過尋找一組輸入變量和響應變量之間的響應關系得出近似模型，從而替代真實系統快速給出求解，常用的代理模型有神經網絡模型、SVR模型、MARS模型等[3]。

本文針對某實際工程大壩建立有限元模型，隨機在一定范圍內均勻抽取部分地基彈性模量和水壓狀態組合，利用自研程序進行有限元計算，得出各個狀態下模型節點的位移，并將結果作為訓練和驗證樣本對GRU神經網絡模型進行訓練，構建代理模型以替代有限元模型作為調用正分析模型，并結合狼群算法提高彈性模量的空間搜索能力，實現對大壩地基彈性模量的反演分析。

1 研究理論

1.1 參數反演理論

綜合考慮多種外部荷載的作用，壩體順河向位移disp可看作水壓分量δH、時效分量δH以及溫度分量δT三部分累加[4]，即

(1)

(2)

(3)

(4)

(5)

式中，βi、β1i、β2i、C1、C2、C3、C4均為回歸系數；Hi為壩前水深的i次方；t為測日期數值；t0為起始日期；τ為隨機誤差項。

隨著上游水位的升高，壩體承受的順河向水體壓力同步增長，相應的順河向位移加大。因此水壓分量δH的大小與水位的變動、結構以及地基材料彈性模量和目標位置的坐標存在著強關聯性。混凝土壩本構模型數學表達為

uc=F(E，H，x，y)

(6)

E=[E1，E2，…，En]

(7)

式中，F(*)為有限元模型節點在不同的材料參數和環境荷載作用下與節點位移之間的映射關系；E為有限元模型各個區域不同彈性模量組成的向量；H為壩前水位高度(不考慮下游水位高度)；(x，y)為有限元模型目標節點的坐標；uc為目標節點的順河向位移值。

反演的最終目標是找出合適的力學參數，使得通過有限元計算的目標節點位移值序列同實際測點測量所得位移值所經過分離得來的荷載分量的誤差fe達到最小。誤差fe為

(8)

式中，ucal為計算位移值；utrue為實測位移值。

1.2 GRU神經網絡模型理論

圖1 GRU神經網絡重復模塊

GRU神經網絡模型的核心數學表達[5]為

rt=σ(Wr·[ht-1，xt])

(9)

zt=σ(Wz·[ht-1，xt])

(10)

(11)

(12)

式中，Wr、Wz、W為網絡參數；σ為igmoid激活函數。

1.3 狼群算法理論

狼群算法結構流程如圖2所示。

圖2 狼群算法結構流程示意

狼群算法通過5個步驟來實現求解最優化問題[6]，分別是初始化狼群、游走行為、頭狼召喚、圍攻獵物以及狼群更新，其算法步驟為：

(1)初始化狼群。以D作為目標變量維數，N作為狼的數量，構建一個D×N的解空間。在解空間中初始化第i頭人工狼位置Xi=(xi1，xi2，xi3，…，xiD)，i=(1，2，3，…，N)，若初始目標函數值為Yi=f(X)，則設使目標函數值最大的人工狼為頭狼Ylead。

(2)游走行為。將一部分目標函數值較大的人工狼設為探狼，每個探狼將分別以活動步長stepa在解空間中的h個方向進行游走，每到達一個新位置都將重新計算該位置的目標函數值并返回起點重新出發，當發現某個位置的Yi>Ylead時或游走次數達到上限時，更新頭狼位置。該過程的數學表達為

(13)

(3)頭狼召喚。頭狼會對猛狼進行召喚，猛狼會以步長stepb向頭狼逼近，若在逼近過程中猛狼的目標函數值大于頭狼的目標函數值，當頭狼與猛狼的距離小于設定閾值時，進入算法下一步。猛狼的逼近過程計算式為

(14)

(4)圍攻獵物。將頭狼的位置標記為獵物位置，猛狼聯合探狼對獵物進行步長為stepc攻，圍攻過程為

(15)

式中，λ為[-1，1]之間的隨機數。stepa、stepb、stepc的關系為

(16)

式中，S為步長因子；[mind，maxd]為第d個變量的取值范圍。

通過以上5個步驟不斷地往復循環，算法在迭代次數達到預設的最大值或算法的尋優個體達到精度要求時結束。

1.4 代理模型構建與彈性模量反演分析

利用系統抽樣的方法抽取部分1.1節所述的計算分量，確定最佳的GRU網絡參數，從而構建GRU代理模型，并結合狼群算法對地基彈性模量進行反演。

1.4.1 GRU代理模型構建

選取目標工程的典型斷面，通過Hypermesh軟件，建立該斷面的有限元模型F。在合理的材料力學參數范圍內，選取m組不同的參數組合，同時再從斷面的壩底到壩頂的范圍內隨機選取n個水位高度H。

在有限元計算步驟中，通過依次選取不同的混凝土結構材料的彈性模量和不同高度的壩前水位，相應修改有限元模型F中的單元材料參數和水壓荷載參數。經過有限元方法進行力學運算之后，根據目標節點的坐標(x，y)選取有限元模型對應位置的節點的計算位移值uc，依次選取所有目標節點的位移值之后，將計算結果和輸入變量以[E，H，x，y?uc]的格式儲存為樣本，即

uc=F(E，H，x，y)

(17)

根據上述樣本構造流程可知，樣本[E，H，x，y?uc]中，[E，H，x，y]為自變量，確定了網絡的輸入節點的數量，[uc]為因變量，確定了網絡的輸出節點的數量。因為不同的網絡結構訓練結果影響較大，且該結構的確定不具備明顯的推理性，需要人為選取不同的的網絡結構等超參數，即中間各層的節點數量，以及各層的激活函數，根據損失函數的種類來試算不同的迭代步長α，最大迭代次數N，以及迭代閾值，經比較后選擇最佳的網絡結構狀態。因為一批訓練中，較大的樣本數量可以促使網絡快速收斂至較優參數狀態，但是對于服務器的計算內存要求較高，而同一批次中較小的樣本數量會使得模型訓練時間加長，且網絡收斂狀態有偏差，所以根據實際GPU的內存大小來確定每批合適的訓練數量m。在0～1范圍內對網絡的權值和閾值參數矩陣W進行隨機初始化，為了防止網絡失效，每個參數的初始化均需獨立計算。將樣本中自變量組成的向量[E，H，x，y]導入輸入節點，進入后向傳播階段，通過每層節點之間的網絡參數矩陣W對自變量進行組合計算，并從網絡輸入端傳至輸出端，得到每個樣本對應的網絡計算值O，即

O=f(E，H，x，y)

(18)

將O結合因變量uc計算本批次樣本的損失值loss，即

(19)

最后，進入網絡前向傳播步驟，對各個網絡參數W進行偏導計算，結合步長α對網絡權重及閾值參數進行調整。進而繼續重復后向傳播和前向傳播，直至N次迭代終止。當網絡參數W達到當前最佳狀態，保存網絡結構及參數，即GRU代理模型構造完畢。

1.4.2 地基彈性模量反演流程

構建GRU代理模型后，根據狼群算法理論，先對彈性模量進行初始化設置，將初始彈性模量通過GRU代理模型計算得到位移值并與實際位移進行對比，若誤差不滿足要求且迭代次數未達到設定的上限，則更新作為頭狼的彈性模量，在多次迭代尋優后將輸出表現最優的頭狼的位置定為最佳彈性模量，該計算流程如圖3所示，該計算流程簡稱為GRU-WPA法。

圖3 GRU-WPA法反演分析計算流程示意

表1 地基材料巖體力學參數

2 GRU-WPA法反演實例分析

選取某碾壓混凝土壩為模型驗證研究對象，該壩壩底高程41.00 m，壩頂高程153.00 m，最大壩高112 m，壩頂寬度6.0 m，上游折坡頂點高程為84.0 m，坡度為1∶0.3，下游面坡度1∶0.75，折坡點高程為145.00 m。該工程共計10個壩段。首次蓄水日期為2011年4月1日。地基材料力學參數如表1所示。在壩體壓力以及地下水長期作用下，該工程的順河向位移在監測期內呈現緩慢上升狀態，因此對該處壩段地基材料參數應當重點關注，本文模型驗證目標為地基彈性模量。

2.1 研究思路

本文計算選取該建筑物的某一個壩段，該位置地基高程為45.5 m，壩高107.5 m，壩基長88 m，選取地基長488 m，地基深300 m。針對該壩段的倒垂線進行分析，該儀器于2014年布置于靠近壩體上游面位置。選取該倒垂線從2014年7月25日到2019年10月31日監測的共計221個順河向位移測值。選定地基主要由石英砂巖構成，材料密度為2 500 kg/m3，泊松比為0.167；壩體材料密度為2 400 kg/m3，泊松比為0.167。

本文抽取部分監測的位移值通過處理輸入有限元軟件進行狀態組合計算，將計算結果結合假定的彈性模量作為神經網絡的訓練和驗證樣本，構造初始參數的GRU代理模型，如圖3的步驟，利用狼群算法繼續對彈性模量進行搜索尋優，從而得到最終的反演彈性模量。

為驗證模型反演結果的合理性，本文利用工程常用的粒子群算法(PSO法)對該地基彈性模量進行反演分析，對比反演結果的差異。

2.2 GRU代理模型樣本構造

因為地基的實際彈性模量未知，訓練樣本應包含可能的目標數值，所以根據表1所示地基可能的總體彈性模量，從3～10 GPa范圍內隨機均勻抽取200個不同數值。針對共計44 200種材料狀態與水壓狀態組合，采用有限元軟件GeHoMadrid對上述模型進行所有的狀態組合計算，得出各個狀態組合下模型節點的位移值disp，儲存格式為[E，H，x，y，disp]，作為網絡模型的訓練和驗證樣本。

2.3 彈性模量反演

模型采用的多元線性回歸模型見式(5)。根據該式提取水壓分量構建目標樣本格式為[E，H，x，y，δH]，其中，E為隨機設定初始值；H為實測水位高度；(x，y)為點D有限元模型中的坐標；δH為上述模型提取的水壓分量。

將取得的樣本排列順序進行隨機打亂以后，所有數據按特征分別歸一化到[0， 1]范圍，抽取樣本對網絡進行訓練驗證，經過訓練后保留網絡參數以便后續調用。

經過如1.4節所示的過程得出訓練完成后的GRU代理模型后，代替有限元模型進行基于狼群算法的彈性模量反演分析計算。利用系統抽樣的方法選取圖4中有限元模型測點A的221個不同水位的樣本序列作為計算樣本，其水位變化范圍為125.33～145.96 m，迭代彈性模量預測位移與實測位移的平均絕對誤差變化過程見圖5。

圖4 有限元模型

圖5 WPA-GRU代理模型反演計算適應度迭代過程

由圖5可知，迭代彈性模量預測位移值與實測值的平均絕對誤差值從0.151 mm下降到0.137 mm，平均絕對誤差在計算全程都不大，說明計算精度較高，且初定50次迭代的計算在迭代30次時數據已趨于穩定，主要是因為狼群算法具有較強的全局搜索能力。

為了進行對比，將同樣的計算樣本利用PSO法進行反演計算，其中設定粒子數為50個，個體學習因子為0.8，群體學習因子為0.2，慣性權重為0.5。由于每需要調用有限元模型50次迭代計算用時約4 min，而GRU-WPA反演分析的方法在30 s內即可完成50次迭代計算。

GRU-WPA代理模型反演計算得到的地基彈性模量為5.27 GPa，PSO法反演計算得到的彈性模量為5.16 GPa，結合不同的水位與兩種不同反演彈性模量進行有限元計算，得到的計算位移結果與實測位移進行對比，位移值和相對誤差結果如圖6所示。由圖6a可知，實測位移值與GRU-WPA反演彈性模量有限元計算位移值整體貼合度較好，僅在少數樣本點處出現了較為明顯的分離情況，說明這兩組數據序列在相同水位狀態下的數值較為接近；而PSO算法的反演彈性模量有限元計算位移值與實測數據整體變化趨勢相近但有一定誤差。由圖6b可知，GRU-WPA反演彈性模量的計算位移值誤差大部分都處于較低水平，大部分都遠低于10%，部分峰值也不會超過40%，而PSO法反演彈性模量計算位移相對誤差大部分略超10%，最大相對誤差達到了68.9%，顯然GRU-WPA反演分析計算結果具有更高的精度。

圖6 實測位移與不同反演方法計算位移值的結果對比

3 結 論

(1)與PSO法相對比，為提高重力壩地基彈性模量反演計算的效率而提出的GRU-WPA代理模型的參數空間搜索能力會更高，由于不需要調用有限元模型，計算用時縮短，提高了計算效率。

(2)通過GRU-WPA反演分析計算得出的反演彈性模量與PSO法反演分析計算得出的反演彈性模量相比相差不大，說明其結果是合理的，且反演彈性模量的有限元計算位移值的相對誤差整體更小，因此，GRU-WPA代理模型的計算精度是可以保證的。

(3)GRU-WPA代理模型在單參數反演中進行了嘗試，后續可繼續開展對于多參數高維度的反演計算研究，探究如何在高維度的情況下保證計算精度和求解效率。