高拱壩參數(shù)反演的Jaya-高斯過程回歸模型

馬建婷，康 飛，姜成磊，向正林，王一帆

(1.大連理工大學(xué)水利工程學(xué)院，遼寧 大連 116024； 2.吉林省水利水電勘測設(shè)計(jì)研究院，吉林 長春 130021；3.南方電網(wǎng)調(diào)峰調(diào)頻發(fā)電有限公司，廣東 廣州 510630)

21世紀(jì)以來，我國高拱壩建設(shè)取得了飛速發(fā)展，先后建成了小灣、錦屏一級、溪洛渡等多個(gè)特大型工程[1-2]。正確估算壩體混凝土和壩基巖石的材料參數(shù)是進(jìn)行拱壩安全評價(jià)與建立大壩安全監(jiān)控模型的關(guān)鍵[3]。拱壩物理力學(xué)參數(shù)反演分析是一個(gè)多參數(shù)組合的大空間搜索問題，求解算法十分重要[4-5]。在拱壩的智能優(yōu)化反演分析方面，劉耀儒等[6]將改進(jìn)的并行遺傳算法用于溪洛渡拱壩位移反演分析；李火坤等[7]采用粒子群優(yōu)化算法反演得到二灘拱壩的實(shí)際動彈性模量；Kang等[8]提出了一種求解反分析問題的混合單純形人工蜂群算法；練繼建等[9]建立了反映大壩彈性模量與位移之間非線性關(guān)系的BP神經(jīng)網(wǎng)絡(luò)模型，并將其用于李家峽拱壩材料參數(shù)反演；Liu等[10]在反演過程中采用無約束拉格朗日支持向量回歸來代替大型有限元計(jì)算，以獲得彈性模量與變形的關(guān)系。

目前的反演分析方法存在兩方面不足：①基于有限元模型求解大壩靜動力響應(yīng)往往需要耗費(fèi)大量的計(jì)算時(shí)間，對于某些大型結(jié)構(gòu)甚至無法求解；②隨著新型群體智能算法的提出，反演分析的精度還有待進(jìn)一步提高。為此，本文將Jaya智能優(yōu)化方法以及高斯過程回歸(Gaussian process regression，GPR)理論引入大壩安全監(jiān)控領(lǐng)域，提出了基于Jaya-高斯過程回歸代理模型的拱壩物理力學(xué)參數(shù)反演分析方法。

1 混凝土壩反演分析基本理論

對于大壩工程，建立大壩與壩基位移監(jiān)控模型以后，依據(jù)正分析的結(jié)果，通過相應(yīng)的理論反求大壩和地基的材料物理力學(xué)參數(shù)，即為大壩位移反演分析[11]。若計(jì)算分析方法采用有限元法，則反演分析即為大壩的有限元模型修正。有限元模型修正是指利用靜動力測試結(jié)果修改有限元模型的剛度、質(zhì)量、邊界條件等參數(shù)，修正后的模型不僅能夠反映結(jié)構(gòu)的實(shí)際狀態(tài)，還可以用于損傷識別、安全評價(jià)等方面。

優(yōu)化反演分析法是一種直接的反演分析方法，其適用性較強(qiáng)，通過建立目標(biāo)函數(shù)，將反演問題轉(zhuǎn)化為參數(shù)最優(yōu)化問題，通過不斷迭代逐次修正試算值，直至獲得最佳值[12]。靜力計(jì)算采用有限元模型，基于位移監(jiān)測數(shù)據(jù)分離水壓分量反演壩體參數(shù)的數(shù)學(xué)模型可表達(dá)為

(1)

式中：K為剛度矩陣；u為位移向量；R為荷載向量；ωi為權(quán)重，通常取值為1；Δuie為實(shí)測位移的變化值；Δuin為有限元計(jì)算得到的位移變化值；m為測點(diǎn)數(shù)目。

2 基于Jaya-高斯過程回歸代理模型構(gòu)建

2.1 Jaya算法

Jaya算法是一種新型智能優(yōu)化算法，由Rao[13]于2016年提出。Jaya算法設(shè)計(jì)的基本思想是群體應(yīng)趨向最優(yōu)解移動而遠(yuǎn)離最差解。該算法，除了常規(guī)的群體規(guī)模和搜索終止條件參數(shù)，無需輸入其他特定參數(shù)；同時(shí)，更新方程只有一個(gè)，算法簡單、易實(shí)現(xiàn)。

Jaya算法的數(shù)學(xué)描述如下:令g(z)為待優(yōu)化的目標(biāo)函數(shù)，在任意迭代次數(shù)a中，Zj=(zb1,zb2,…,zbD)(b=1,2,…,N)為第b個(gè)候選解(D為設(shè)計(jì)變量的個(gè)數(shù)，N為群體規(guī)模)，Pb=(pb1,pb2,…,pbD)為第b個(gè)候選解尋找的最優(yōu)解。通過式(2)生成新的候選解，不斷更新，逐漸收斂到最優(yōu)解。

Za+1,b=Pa,b+r1(Za,best-|Pa,b|)-r2(Za,worst-|Pa,b|)

(2)

式中：Za,best、Za,worst為第a次迭代整個(gè)群體中的最優(yōu)解和最差解；r1、r2為在(0,1)內(nèi)的2個(gè)隨機(jī)數(shù)。

2.2 GPR

GPR最早由Rasmussen等[14]提出，是一種機(jī)器學(xué)習(xí)新回歸方法。其最大的優(yōu)點(diǎn)是能夠從訓(xùn)練數(shù)據(jù)中學(xué)習(xí)平滑參數(shù)和噪聲以提供不確定性估計(jì)[15]。

GPR預(yù)測均值和協(xié)方差為

(3)

cov(f*)=K(X*,X*)-

(4)

由(3)式可知，協(xié)方差函數(shù)對GPR預(yù)報(bào)結(jié)果有一定的影響。因此，協(xié)方差函數(shù)的選擇對于GPR建模的成功具有重要意義[15]。常用的協(xié)方差函數(shù)為

(5)

2.3 反演分析步驟

反演分析基本步驟如下：

步驟1生成樣本點(diǎn)。確定待反演區(qū)域彈性模量的取值范圍，采用拉丁超立方抽樣生成樣本點(diǎn)。

步驟2構(gòu)建樣本集。樣本集由樣本點(diǎn)及其各監(jiān)測點(diǎn)的相對位移組成。采用有限元進(jìn)行數(shù)值模擬，計(jì)算出各監(jiān)測點(diǎn)的相對位移。

步驟3建立GPR代理模型。利用前2步得到的訓(xùn)練樣本集建立基于GPR的代理模型。

步驟4檢查GPR代理模型的準(zhǔn)確性。利用步驟1和步驟2得到的測試樣本，采用均方根誤差(RMSE)作為樣本精度的評價(jià)指標(biāo)[16]，檢驗(yàn)所建立的GPR代理模型的準(zhǔn)確性。如果精度不滿足，則返回到步驟1，生成更多的訓(xùn)練樣本點(diǎn)，得到更準(zhǔn)確的代理模型。

步驟5參數(shù)優(yōu)化反演。利用可代替有限元計(jì)算的GPR代理模型以及引入的智能優(yōu)化算法進(jìn)行參數(shù)優(yōu)化反演，得到最優(yōu)解，即為反演分析結(jié)果。

3 算例驗(yàn)證

某特高拱壩最大壩高277 m，為混凝土雙曲拱壩，有限元模型如圖1所示，壩體和基礎(chǔ)共分為11 013個(gè)單元，51 917個(gè)結(jié)點(diǎn)。采用有限元正分析的計(jì)算結(jié)果來模擬大壩位移觀測值的水壓分量[3]。選取的大壩變形監(jiān)測點(diǎn)分別為P1、P2和P3。選取一組彈模E=(E1，E2，E3，E4)=(21，20，13，15)GPa，分別為壩體、底部基巖、左岸基巖和右岸基巖的彈性模量。求得上游水位為820 m和790 m的2種工況下P1、P2、P3測點(diǎn)的相對水平位移分別為16.60 mm、3.22 mm、-2.64 mm，相對垂直位移分別為42.79 mm、45.23 mm、24.81 mm。反演時(shí)，選取壩體和基巖的彈模范圍為5～40 GPa[16]。

圖1 某特高壩體有限元模型

表1 算例3種算法反演結(jié)果

為建立Jaya-高斯過程回歸代理模型，研究了建模樣本數(shù)量對精度的影響，將訓(xùn)練樣本和測試樣本(取為訓(xùn)練樣本的25%)的RMSE繪成曲線，如圖2所示。根據(jù)圖2可知，均方根誤差呈下降趨勢，樣本的精度逐漸增加。當(dāng)樣本數(shù)量為100D時(shí)，均方根誤差變化較小。因此，選取樣本數(shù)量為100D的GPR代理模型代替有限元模型。

圖2 樣本均方根誤差趨勢

在該高拱壩反演分析中引入Jaya算法，并將它與粒子群算法[17](particle swarm optimization，PSO)及灰狼算法[18](grey wolf optimizer，GWO)做比較。3種算法群體規(guī)模N=20，最大迭代次數(shù)為200，粒子群算法的其他參數(shù)參考之前的研究[19]設(shè)置。為了避免算法的隨機(jī)性，將優(yōu)化算法運(yùn)行30次[16]。3種算法的反演結(jié)果見表1，最優(yōu)適應(yīng)值變化曲線和壩體混凝土彈性模量識別過程曲線分別如圖3、圖4所示。

圖3 算例拱壩最優(yōu)適應(yīng)值變化曲線

圖4 算例拱壩彈性模量E1識別過程曲線

由表1可知，本文提出的反分析策略在該拱壩反演分析中是可行的。且Jaya算法比PSO算法、GWO算法的精度更高。由圖3及圖4可知，達(dá)到收斂時(shí)，Jaya算法的迭代次數(shù)約為80次，而PSO算法的約為130次，GWO算法卻始終沒有達(dá)到收斂狀態(tài)。此外，在30次的運(yùn)行結(jié)果中，Jaya算法和PSO算法的反演結(jié)果較為穩(wěn)定，相對誤差變化較小；而GWO算法反演結(jié)果的相對誤差變化較大，最大誤差為9.314 6%，最小誤差為0.467 1%。

為了研究3種優(yōu)化算法的抗噪聲能力，在監(jiān)測點(diǎn)的位移分別加3%的噪聲，同樣運(yùn)行30次，結(jié)果見表2。由表2可知，在加入3%噪聲時(shí)，Jaya、PSO及GWO算法的反演結(jié)果的平均相對誤差分別為2.651 0%、2.651 1%和3.285 6%，均較小。因此，在加入3%噪聲時(shí)3種算法對該混凝土拱壩的反演結(jié)果是可靠的。

表2 加入3%噪聲3種算法反演結(jié)果

在反演用時(shí)方面，由表3可知， GPR代理模型比有限元模型節(jié)省約87.48%的時(shí)間。由此可以看出，基于Jaya-GPR代理模型的拱壩物理力學(xué)參數(shù)反演分析方法在保持計(jì)算精度的同時(shí)，計(jì)算效率也具有較大的優(yōu)勢。

表3 算例反演用時(shí)計(jì)算結(jié)果

4 實(shí)例應(yīng)用

以西南某高拱壩為例，該高拱壩壩頂高程610.0 m，最大壩高285.5 m，頂拱中心線弧長681.5 m，厚高比為0.216，弧高比為2.451。依據(jù)壩體幾何參數(shù)及設(shè)計(jì)相關(guān)資料，建立有限元分析模型。其中壩基范圍為：向上下游、左右岸和河床深部分分別延伸1.5倍壩高，整個(gè)模型劃分為233 538個(gè)單元，261 734個(gè)結(jié)點(diǎn)，有限元網(wǎng)格如圖5、圖6所示。邊界條件定義為：壩基底部為三向固定約束，4個(gè)側(cè)面為法向固定約束。其中壩體混凝土和壩基的材料參數(shù)除待反演的彈性模量，其余值均采取設(shè)計(jì)值。

圖5 拱壩整體有限元模型

圖6 壩體有限元模型

于拱壩左岸壩肩、5號壩段、10號壩段、15號壩段、22號壩段、27號壩段、右壩肩拱布置7組29個(gè)正垂線測點(diǎn)用以監(jiān)測壩體水平位移(圖7)。根據(jù)監(jiān)測數(shù)據(jù)顯示7組測線中位于壩體中心線的15號壩段垂線監(jiān)測到的壩體水平位移最大，且由于切向位移較小，實(shí)際觀測值的誤差相對較大，因此本文將以15號壩段徑向位移監(jiān)測數(shù)據(jù)作為研究對象開展分析。

圖7 拱壩正垂線及變形測點(diǎn)布置

對于15號壩段的溫度監(jiān)測數(shù)據(jù)，取與它相鄰的16號壩段布置在壩體表面的32個(gè)溫度傳感器和布置在壩體內(nèi)部的19個(gè)應(yīng)力應(yīng)變傳感器，溫度傳感器分布圖如圖8所示，從圖8可見拱壩溫度測點(diǎn)較多，在應(yīng)用回歸模型時(shí)難以考慮全部溫度測點(diǎn)，經(jīng)過分析削減以后，仍有15個(gè)，故考慮提取主成分的方式，采用SPSS軟件，對15個(gè)已有的溫度因子進(jìn)行降維處理，提取前4個(gè)主要成分(保留溫度測量信息的90.599%)。按照傳統(tǒng)多元回歸模型的因子選擇方法，最終自變量因子設(shè)計(jì)為

V=(h,h2,h3,h4,TPC1,TPC2,TPC3,TPC4,

-exp(-0.5t))

(6)

圖8 拱壩溫度傳感器分布

式中：h為水壓因子，取上游水深；TPC為主成分分析得到的溫度因子；t為時(shí)間，取從監(jiān)測日到始測日的天數(shù)除以365。選取測點(diǎn)PL15-1及測點(diǎn)PL15-2在2014年7月6日至2017年11月29日的實(shí)測徑向位移建立統(tǒng)計(jì)回歸模型。其中實(shí)測值、擬合值及各位移分量結(jié)果如圖9所示。為計(jì)算簡便，荷載工況選取上游水位為599.76 m和559.44 m情況下的水壓差值作為位移反分析的荷載條件，其計(jì)算結(jié)果是在水壓荷載作用下拱壩的位移差值。

圖9 拱壩PL15-1測點(diǎn)各位移分量分離結(jié)果

徑向位移對應(yīng)拱壩的極坐標(biāo)系，需將其進(jìn)行坐標(biāo)轉(zhuǎn)換，得到有限元模型坐標(biāo)系下的位移值，即Δu15-1x= -25.371 8 mm，Δu15-2x= -23.047 2 mm。

根據(jù)相關(guān)試驗(yàn)資料，混凝土和巖體彈性模量的取值范圍設(shè)置為：壩體E1=35～60 GPa，地基E2=15～40 GPa。訓(xùn)練樣本數(shù)量為100D，建立GPR代理模型。3種算法群體規(guī)模N=20，最大迭代次數(shù)為200，每種算法隨機(jī)運(yùn)行30次[16]。3種算法彈性模量的反演結(jié)果見表4，最優(yōu)適應(yīng)值變化曲線及彈性模量識別過程曲線分別如圖10和圖11所示。

圖10 實(shí)例拱壩最優(yōu)適應(yīng)值變化曲線

圖11 實(shí)例拱壩彈性模量識別過程曲線

為了進(jìn)一步驗(yàn)證反演結(jié)果，將上述反演結(jié)果，用有限元計(jì)算壩體測點(diǎn)在上游水位為599.76 m和559.44 m狀態(tài)下的位移差值，并與實(shí)測位移水壓分量差值進(jìn)行比較。從表4中看出，兩者接近，總體誤差較小，說明反演的成果是合理的。

表4 實(shí)例3種算法反演結(jié)果

在反演用時(shí)方面，由表5可知，采用GPR代理模型相比有限元模型節(jié)省時(shí)間約93.75%。進(jìn)一步證明了所提出的基于Jaya-GPR代理模型能大幅度提高體系反演分析的效率。

表5 實(shí)例反演用時(shí)計(jì)算結(jié)果 單位:s

5 結(jié) 論

a.高斯過程回歸用于混凝土拱壩反演分析的代理模型是有效的，與有限元模型計(jì)算結(jié)果一致。

b.與其他算法相比，Jaya算法具有控制參數(shù)少、反演精度高、收斂速度快、穩(wěn)定性好等優(yōu)點(diǎn)。

c.提出的基于Jaya-GPR代理模型的拱壩物理力學(xué)參數(shù)反演分析方法能夠大幅度的提高反演分析的效率，相比于基于傳統(tǒng)有限元的反演方法，用時(shí)節(jié)省了80%以上。