基于改進鯨魚優化算法的重力壩變形參數反演

高林鋼，李同春,2，林潮寧，朱致遠，鄭 斌

(1.河海大學 水利水電學院，江蘇 南京 210098; 2.水安全與水科學協同創新中心，江蘇 南京 210098)

1 研究背景

隨著水利資源的持續開發，越來越多的大壩投入建設與運行，建壩的選址條件也愈發復雜，而大壩的規模也越來越高。大壩的建設投資十分巨大，但創造的效益亦非常可觀，在實現其工程價值的同時，安全管理工作也很重要。在混凝土壩的安全評價中，常采用數值模擬方法研究大壩的強度和穩定性[1]，而壩體和基巖力學參數選取的合理程度[2]，是影響數值模擬準確性的關鍵因素[3]。由于環境、施工、老化、徐變等復雜因素的影響，壩體及基巖材料的力學參數設計值往往與實際值存在一定差異[4]。如果僅利用原位試驗或實驗室測得的參數值，由于受到采樣范圍和試驗條件等因素的限制[5-6]，有失代表性，難以滿足數值分析的精度需求[7]。此外，現場取樣試驗成本高昂，工作復雜，并會對原有結構造成破壞[8]。因此，基于位移實測資料的參數反演分析方法成為識別混凝土壩壩體及基巖力學參數的主要方法之一[9]。在早期的研究中，Bonaldi等[10]使用確定性模型反演混凝土彈性模量、溫度線膨脹系數；吳中如等[11]利用臨界荷載法和小概率試件法，反演混凝土的斷裂韌度；蘇懷智等[12]采用遺傳模擬退火算法反演大壩及壩基的物理力學參數。目前，反演分析方法主要有傳統比值法和優化反演法。

由于難以建立待求解力學參數與觀測數據之間的顯性關系[13]，因此在反演過程中需要重復進行數值模擬過程[14]，其計算量大，耗時長，收斂速度慢，且容易陷入局部極值[15]。近年來，隨著計算機技術的快速發展，遺傳算法(GA)、粒子群算法(PSO)、BP神經網絡等智能算法被應用于各類參數反演分析中[16]。龔曉雯等[13]利用改進的粒子群算法對混凝土重力壩參數進行反演；魏博文等[14]將ANSYS有限元程序作為子模塊嵌入到自適應遺傳粒子群算法(GA-PSO)中對混凝土壩參數進行反演；Zhuang等[17]提出了基于支持向量機回歸的位移反分析模型，并將多群策人工魚群算法(MAFSA)應用于隧洞參數識別；劉健[18]將BP神經網絡進行改進之后，應用于壩體位移測值的處理，從而反演相應參數。在上述這些算法中，遺傳算法經常遇到計算效率低、耗時長的情況，而粒子群算法在求解復雜函數的優化問題時，也往往出現收斂過早、精度低、魯棒性差等問題[9]。

鯨魚優化算法(WOA)是一種新型的群智能優化算法，具有參數容忍度強、魯棒性好、收斂速度快等優點[19]。但是智能算法在求解尋優過程中都會表現出一定的趨同性，易使搜索過程陷入局部最優。一般情況下，增加種群數量，可以降低搜索過程的趨同性，提高求解質量。然而在串行計算中，增加種群會增加計算時間，使計算效率低下，尤其是在反演基巖和壩體力學參數問題中，優化算法需調用有限元程序進行正分析，當種群增大時，有限元計算的次數也隨之增加，計算時間增加更加明顯。利用計算機多核處理器進行并行計算，線程控制簡單，計算環境穩定，正逐漸被應用于科學計算中。

本文將多核并行計算的方法引入標準的鯨魚優化算法中，并加入非線性收斂因子，提出了基于多核處理器的并行式改進鯨魚優化算法。在此基礎上建立了重力壩壩體、基巖力學參數與監測點位移之間的非線性映射模型。最后，將所建立的模型應用于某混凝土重力壩力學參數反演中，并與用粒子群算法反演的結果進行比較，以驗證該模型的可行性與合理性。

2 重力壩參數反演基本原理

在水壓力、泥沙壓力、揚壓力、溫度荷載等外部荷載的共同作用下，混凝土重力壩任一點的變形δ按照形成原因可分為水壓分量δH、溫度分量δT和時效分量δθ3部分[20]，即：

δ=δH+δT+δθ

(1)

將位移監測點實測變形數據進行回歸分析，分離出水壓分量、溫度分量和時效分量，可建立重力壩位移的HST統計模型。在總位移δ中扣除溫度分量δT和時效分量δθ之后，得到的水壓分量δH為：

δH=δ-δT-δθ

(2)

(3)

壩體變形是最能直接反映大壩安全狀態的指標，一般采用倒垂線作為混凝土壩變形監測的基準線，倒垂線從壩體廊道一直貫穿到基巖，不易受外界干擾，較為可靠，其測值可以反映壩基的變形。在外部荷載相似的情況下，材料的彈性模量是影響大壩變形的主要因素，當數值計算中采用的彈性模量與真實的彈性模量較為接近時，計算得到的大壩位移與實測值也比較接近；當數值計算中采用的彈性模量與真實彈性模量相差較大時，計算得到的測點位移值也會與實測值相差較大。

(4)

X=(Ec,Er)T

(5)

式中：P為荷載組合，有限元計算結果中的水壓分量利用每年高水位和低水位情況下的測點位移差值計算得到。

對有限元計算得到相應測點的位移與從實測值中分離得到的水壓分量求殘差加權平方和，作為優化算法的目標函數。

(6)

式中：n為目標函數中包含的位移監測點的數量；wi為第i觀測點測量值的權重。

壩體材料及基巖的彈性模量反演的步驟如下：

(1)對待反演壩段的倒垂線測值進行回歸分析，建立重力壩位移的HST統計模型，分離得到δH、δT和δθ。

3 改進鯨魚優化算法

鯨魚優化算法(whale optimization algorithm, WOA)是Mirjalili等[19]在2016年提出的一種新型元啟發式算法，該算法模擬了自然界中座頭鯨捕食近水面魚蝦的行為，采用了29個測試函數對該算法進行測試，結果表明，無論是收斂速度還是尋優能力，與PSO、GA、GSA等算法相比，WOA算法均具有一定優勢。

鯨是世界上最大的哺乳動物，座頭鯨是鯨類中極具代表性的一種。座頭鯨捕食時具有一種特殊的的狩獵方法，座頭鯨在海洋中水深12 m左右的位置游動并搜尋食物[21]，當它發現目標后，會開始沿著螺旋軌跡上升，并不停吐出氣泡，形成一張螺旋狀的氣泡網，將目標食物驅趕至近水面處的氣泡網中心，最終捕獲食物。圖1為鯨魚螺旋上升進行捕獵的示意圖。WOA算法模擬了座頭鯨捕獵的行為，通過迭代更新位置來搜尋獵物，最終獲得最優解。

圖1 座頭鯨螺旋上升捕食方法示意圖

在WOA算法中，每一頭座頭鯨都作為種群的搜索代理，可視為一個候選解。座頭鯨捕獵過程分為3個步驟，分別為：獵物包圍、氣泡網狩獵策略、獵物搜索。

3.1 獵物包圍

座頭鯨可以識別獵物的位置，并對其進行包圍。本算法以當前迭代循環最優解作為目標獵物模擬位置，其余搜索代理根據最優解進行更新迭代。該搜索行為的數學表達如下[22]：

D=|C·X*(t)-X(t)|

(7)

X(t+1)=X*(t)-A·D

(8)

式中：t為當前迭代次數;D為搜索代理與當前最優解的差值;X*(t)為當前最優解位置向量;X(t)為搜索代理當前位置向量;A、C為系數向量。

當某次迭代過程中，最優解發生更新時，X*(t)更新為新的最優解。

系數向量A由公式(9)、(10)得到，當|A|≥1時，算法會擴大搜索范圍，進行全局搜索，尋找更優的獵物位置，由隨機產生的搜索代理充當關鍵點，當|A|<1時，算法會縮小搜索范圍，進行局部搜索，以當前最優解位置作為樞紐點。系數向量C由公式(11)得到。

A=2a·r1-a

(9)

(10)

C=2·r2

(11)

式中：r1,r2為[0,1]上的隨機向量;Tmax為總迭代次數，在迭代過程中，a從2到0呈線性遞減。

3.2 氣泡網狩獵策略

座頭鯨在搜索獵物過程中盤旋上升，在上升時吐出大量氣泡形成陷阱，這個階段被稱為氣泡網狩獵策略。公式(9)中a的減小即為收縮環繞機制的模擬，而A的波動范圍也隨著a的減小而縮小。環繞收縮過程可表達為：

X(t+1)=D′·ebl·cos(2πl)+X*(t)

(12)

式中：D′=|X*(t)-X(t)|，b為定義對數螺旋線形狀的常數;l為[-1,1]之間的隨機值。

座頭鯨搜索獵物環繞收縮過程可建立如下更新模型：

(13)

式中：變量p為[0,1]之間的隨機數，p隨機在0和1直接切換。

3.3 獵物搜索

在初始階段，座頭鯨處于隨機搜索獵物階段，在數學表達中，利用隨機數原理進行初始化。

D=|C·Xrand(t)-X(t)|

(14)

X(t+1)=Xrand(t)-A·D

(15)

式中：Xrand(t)為從種群中隨機選擇的搜索代理位置向量。

3.4 權重因子策略

在群智能優化問題設計中，如何平衡全局尋優能力和局部尋優能力一直是研究的熱點。當算法的全局尋優能力較強時，搜索種群能盡可能廣地分布在每個區域，從而能快速找出全局最優解所在的位置范圍，而局部尋優能力較強時，能夠從最優解的范圍內對最優解位置進行快速定位，提高算法的收斂速度和結果精確程度。為了同時保證鯨魚優化算法的全局尋優能力和局部尋優能力，引入權重因子w。當權重因子較大時，算法的全局尋優能力較強，有利于搜索代理進行全局搜索；當權重因子較小時，算法的局部尋優能力較強，能快速收斂，并保證結果精度。

(16)

利用權重因子進行對鯨魚優化算法進行改進之后，座頭鯨向獵物螺旋環繞收縮過程由公式(11)更新為公式(15)。

X(t+1)=

(17)

權重因子w隨著迭代的進行不斷減小，使搜尋過程從前期的有利于全局尋優過渡到后期的有利于局部尋優，鯨魚優化算法改進前后權重因子變化過程對比如圖2所示，由圖2可見，改進后的權重因子w在搜尋初期下降較快，這樣的趨勢更有利于算法在確定最優解范圍之后進行局部尋優，提高算法的收斂速度和結果精確程度。

圖2 鯨魚優化算法改進前后權重因子變化過程對比

3.5 算法的并行改進

隨著計算機軟硬件技術的發展，近年來，并行計算技術逐漸被應用于科學計算中。并行計算與串行計算的區別在于并行計算通過提高空間復雜度的方式來降低時間復雜度，提高算法效率，而串行計算保持算法順序邏輯的簡單，但需耗費更多的時間。并行計算中，計算線程數與CPU內核數相同，通過將任務分解為小部分，以并發方式，分配給CPU的各個內核同時進行計算，并最終將結果合并。標準鯨魚優化算法反演力學參數的過程是串行計算，在一個迭代步中，一個搜索代理搜尋結束，下一個搜索代理才能啟動搜尋，因此同一時間計算機只能進行一次有限元正分析，反演過程耗時極長。為了提高該優化算法反演混凝土壩力學參數的效率，利用并行計算的原理，提出基于并行式改進鯨魚算法反演混凝土壩力學參數的方法，將搜索代理種群分為若干個子種群，將子種群分配到CPU不同內核上執行有限元正分析計算，待各子種群計算結束，將之合并形成新的種群，然后進入下一迭代步。

并行式改進鯨魚優化算法流程如圖3所示。

圖3 并行式改進鯨魚優化算法流程圖

4 實例分析

某水利樞紐主要由碾壓混凝土重力壩、泄水閘、引水建筑物、廠房和開關站等建筑物組成。水庫校核洪水位151.88 m，設計洪水位150.00 m，正常蓄水位130.00 m，死水位130.00 m，總庫容7.173×108m3。碾壓混凝土重力壩壩底高程41.00 m，壩頂高程153.00 m，最大壩高112.00 m，壩頂寬度6.00 m，壩頂長568.00 m。選用該重力壩6#非溢流壩段作為實例進行驗證，并選用該壩段60.00 m高程處倒垂線測點IP8的觀測資料進行計算分析。

該實例計算所建立的有限元模型如圖4所示，地基分別向上、下游延伸100 m，地基深度取為100 m。上游最高水位為150.00 m，最低水位為132.08 m。有限元模型采用四節點四邊形單元，節點1 220個，單元1 135個，x向為順河向，y向為豎直向。

圖4 實例計算所建立的有限元模型

壩體材料為C20碾壓混凝土，密度為2 400 kg/m3，泊松比為0.167，基巖材料參數如表1所示。荷載考慮上游壩面水荷載和庫水重。地基側面邊界為法向約束，底面邊界為固定約束。由于倒垂線IP8深入建基面以下40 m，故數值計算中以倒垂測點所在60 m高程上的點A相對于地基深度40 m處點B的位移作為IP8測點處的計算位移值。

表1 基巖材料參數

先以數值算例驗證本算法的可行性。假設壩體混凝土材料彈性模量為30.0 GPa，基巖巖體彈性模量為15.0 GPa，其余參數均使用默認值。設定水位變化范圍為130.00～150.00 m，水位變化間隔為1.00 m，使用有限元程序計算各水位工況下測點位移，由于外荷載只存在水壓作用，因此有限元計算值全部為水壓分量，可直接作為基準值。使用改進鯨魚優化算法進行參數反演，圖5為反演迭代過程。

圖5 改進鯨魚優化算法反演迭代過程

由圖5可以看出，在迭代次數達到30時，適應度已收斂于0。反演結果為壩體混凝土彈性模量30.0 GPa，基巖巖體彈性模量為15.0 GPa，與預設模量相同，表明算法可應用于混凝土壩參數反演。

將算法應用于本工程實例中，利用實測數據進行反演。倒垂線IP8于2014年7月開始有測值，故分別選取2015、2016、2017、2018 共4個年份所對應的最低水位和最高水位值進行計算。表2為2015-2018年最低和最高水位對應測點的測值，以及最低和最高水位工況的水位差和測值變化。

對60 m高程的6#壩段倒垂線IP8測點的實測數據進行回歸分析，分離得到δH、δT和δθ。在有限元計算中考慮上游壩面水荷載和庫底水壓力，分別計算每年對應最高水位和最低水位情況下，IP8測點對應位置的位移，計算其差值。

將有限元計算得到的水壓分量差值與實測水壓分量差值的離差平方和作為目標函數，分別用鯨魚優化算法和粒子群算法求解最低適應度時的壩體和基巖彈性模量。為提高計算效率，對粒子群算法同樣進行并行化處理，不影響其計算精度。圖6為兩種算法的反演過程對比，表3為2015-2018年6#壩段彈模反演計算結果。

由圖6可以看出，粒子群算法較早進入局部極值，一直在局部范圍內尋找最優解，而改進鯨魚優化算法尋找的適應度比粒子群算法尋找的適應度小。從表3可得，由改進鯨魚優化算法對2015-2018年倒垂線IP8監測數據進行反演得到結果的適應度分別為0.013，0.012，0.018，0.015，由粒子群算法反演得到結果的適應度分別為0.036，0.086，0.031，0.075，表明改進鯨魚優化算法收斂速度比粒子群算法快，其反演結果的適應度也明顯優于粒子群算法，且改進鯨魚算法得到得反演結果更接近于該大壩的實際參數。

圖6 改進鯨魚優化算法和粒子群算法的反演過程對比

表3 改進鯨魚優化算法和粒子群算法的2015-2018年6#壩段彈性模量反演計算結果

5 結 論

(1)通過對鯨魚優化算法進行改進，結合有限元正分析，提出了基于改進鯨魚優化算法反演混凝土壩變形參數的新方法，通過對某碾壓混凝土重力壩實例進行驗證，計算得到該重力壩壩體和基巖彈性模量，其反演結果優于粒子群算法，表明基于改進鯨魚優化算法反演大壩力學參數的方法具有搜索精度高，反演速度快等特點，是合理可行的。

(2)基于多核CPU實現算法的并行計算，以子種群的形式對搜索代理種群進行細分，使復雜的計算能在CPU的多個內核上同時進行，極大提高了利用該算法進行復雜問題計算的求解速率，當程序在計算集群上進行計算時，可在極短時間內得到力學參數反演結果，可應用于大壩監測系統中以實時反饋各部位力學參數。

(3)多個力學參數同時反演時，解具有高維度的特點，如何繼續改進該優化反演方法，使該算法在同時反演多個力學參數時，仍具有較高的搜索精度和收斂速度，是進一步研究的重點。