多目標約束下的土石壩抗液化措施數值建模與優化分析方法

李明超，楊 琳，任秋兵，趙 宇

（1.水利工程仿真與安全國家重點實驗室，天津大學，天津 300350；2.天津水務投資集團有限公司，天津 300200）

1 研究背景

地震是一種影響極大的自然災害，同時會伴隨著嚴重的次生災害，其中地基液化問題較為突出。土石壩作為當今應用最廣泛的壩型，其壩基液化問題不容忽視。一旦壩基發生液化，土體中的孔隙水壓力短時間內迅速升高，有效應力急劇降低，土體的抗剪強度嚴重削弱，導致沉降不均勻、地基失效、壩體滑坡、土壩潰決等一系列災害，對地基、岸坡和壩體等造成極大的破壞［1］。

在目前的土石壩工程中，對于可能發生液化的土層，常常采用挖除、置換、壓重等處理措施［2］。在已有工程中，長河壩礫石土心墻堆石壩挖除壩基幾十米埋深范圍內的可液化砂層；黃金坪瀝青混凝土心墻堆石壩的基礎液化砂層基本挖除，并在下游設置壓重。設置壓重增加了坡腳下的豎向應力，提高了壩基土的抗液化能力；同時也增大了壩體和壩基的整體抗滑穩定性，防止了壩基土體從坡腳擠壓而出。該方式施工簡單，且可以利用開挖的土石料，造價低，經濟有效，因此被廣泛用于我國的中小型土石壩的坡腳加固［3-4］。

近年來，針對土石壩抗液化措施的分析國內外進行了很多研究。Mulunda等［5］主要通過改善液化土的性質和初始應力來控制液化，提出了地表開挖回填、設置壓力平臺、防滲平臺和加強排水設計等液化控制措施；付磊等［6］采用擬靜力法和動力反應分析法相結合確定土石壩深水條件下拋石壓坡的抗震加固斷面；Ye等［7］利用大壩自重固結來提高土石壩抗液化性能，并使用有效應力方法來分析了液化修復的效果；高大水等［8］采用動力有限元分析對花涼亭工程上游壩坡幫坡壓重抗震加固方案進行技術論證，有效解決了水庫蓄水條件下的大壩抗震加固問題；岑威鈞等［9］根據地基液化深度確定液化水平影響長度，再進行不同填土厚度下抗震效果分析，最終確定了安全經濟的壓重平臺厚度和長度。Zhang等［10］探討了土石壩壓重平臺的長度、高度和材料內摩擦角對土石壩靜動力特性的影響，并給出了建議的取值范圍；Amin 等［11］進行穩定滲流條件下Siahoo 大壩的地震動力分析，識別地基的液化區并采用應力重分布方法估計震后位移；Kazunari 等［12］提出了一種根據地質統計參數及地震參數的土石壩液化概率計算方法，并應用于實際設計問題中；王富強等［13］分析了孔隙水轉移對液化土層強度和變形的影響，提出了壩基覆蓋層地震液化的挖除置換、加密、壓重等工程措施；Hadi 等［14］通過振動臺試驗研究了橡膠排水柱和碎石排水柱對飽和砂土液化勢的影響，表明密實度是控制液化的最重要因素；Liam Finn［15］綜合薩迪斯大壩等多個工程，提出了一套包括性能標準選擇、分析程序及非線性動力本構的液化場地土石壩設計方法。

然而，已有研究主要是根據工程經驗確定土石壩抗液化措施的方案，目前在抗液化措施效果優化方面的研究較少?，F有工程中，主要根據設計人員經驗選取主要設計參數，沒有具體的理論設計指導，若參數選擇過于保守，則增加工程投資，延長工期，增加了施工難度；若參數選擇過于經濟，則仍然存在土石壩安全隱患。因此，若能提出一套可以綜合考慮安全性和經濟性的土石壩的抗液化措施優化設計方法，實現對土石壩抗液化措施的定量化、精細化設計，具有重要意義和實踐價值。

針對復雜的工程優化問題，作為常用的機器學習算法之一，遺傳算法（Genetic Algorithm，簡稱GA）通過建立對應具體問題的適應度函數并利用遺傳算子聯合搜索的方法求解最優解，具有原理簡單、適應性廣的特點，在工程優化領域逐漸得到廣泛應用。由于工程中優化問題常常需要考慮成本和效果的平衡，一般都屬于多目標優化問題。多目標遺傳優化算法作為一種群體優化算法，對多目標優化問題的求解與傳統優化方法相比有明顯優勢［16］。例如，諶文武等［17］提出一種以性價比為評價指標，運用快速拉格朗日差分分析和遺傳算法對斜坡加固方案進行優化的新方法；官夏菲等［18］采用SVR 回歸模型建立優化變量（跌坎高度d、突擴比β和尾坎坡度θ）與優化目標（消能率ΔE/E1、臨底流速v）之間的近似模型，并采用遺傳算法求解該近似模型，對消力池體型進行優化；Trung-Thanh［19］等通過存檔微遺傳算法（Archive-based Micro-Genetic Algorithm，簡稱AMGA）與克里金模型的結合對加工參數進行優化，改進了銑削過程。

總的來說，目前工程抗液化措施主要根據經驗設計，尚未建立定量的土石壩抗液化措施優化分析方法。因此，本文有針對性地開展了多目標、多約束條件下綜合安全性與經濟性的土石壩抗液化措施數值優化研究。以壩體坡面水平位移梯度、壩頂震陷和抗液化措施成本最小化為優化目標，以靜動力抗滑穩定安全系數、地基最大液化范圍、液化深度、典型點地震過程中平均液化度和震后永久水平位移為約束條件，建立土石壩抗液化措施多目標優化數學模型；進而通過VC++集成有限差分法軟件FLAC3D進行開發，運用完全非線性動力算法和多目標遺傳算法對該數學模型進行求解，并將其應用于工程實例，為改善土石壩壩基液化問題提供新的手段。

2 研究方法框架

綜合考慮經濟性和安全性，根據壓重平臺和地基置換的特點建立相應的多目標優化數學模型，采用多目標遺傳優化算法AMGA進行求解，對壓重平臺和地基置換的設計參數進行優化，研究框架如圖1所示，主要包括以下內容：

（1）根據以往工程經驗，確定需要優化的設計參數（壓重平臺的長度、高度和坡度以及地基置換的深度和位置等），基于工程資料確定參數變量空間，即參數的取值范圍。優化過程中，算法在該范圍中生成設計變量的方案。

（2）基于設計變量方案，進行參數化建模并采用FLAC3D進行土石壩靜動力計算及壩基液化的分析。首先實現了土石壩填筑及蓄水的靜力模擬，在此基礎上分別進行了安全系數及動力部分的計算。安全系數包括靜力抗滑安全系數及擬靜力法動力抗滑安全系數，在動力計算中編寫了壩基液化監控程序，計算得到壩基的液化程度、范圍等情況。

（3）綜合考慮經濟性與安全性，選取壩體坡面水平位移梯度、壩頂震陷與抗液化措施方案的造價作為優化目標，從地基液化程度、液化范圍及大小、靜動力抗滑穩定安全系數、壩體永久變形的方面給定約束條件，建立土石壩抗液化措施的多目標優化數學模型。

（4）運用多目標遺傳算法AMGA 控制求解上述流程。通過VC++集成FLAC3D 進行開發，實現算法自動搜索設計變量空間，控制FLAC3D進行土石壩靜動力及壩基液化的計算，并在后處理程序中計算目標函數及約束條件，多次迭代直至算法收斂。該算法將求出一組Pareto 解集，根據工程需要從求出的Pareto 解集中選取最優解。將該優化方法應用于工程實例，實現對某土石壩抗液化措施的多目標優化，并通過對比分析驗證了該方法的有效性。

圖1 研究方法總體框架

3 土石壩抗液化措施多目標優化數學模型

（1）確定設計變量及變量空間。以土石壩工程中采用壓重平臺及地基開挖置換措施為例，選取9個控制性設計變量，見表1和圖2。具體的參數空間，即每個變量的變化范圍及初始值，需根據特定的工程確定。

（2）給定約束條件，選取優化目標。常見的單目標壩體優化主要是以滿足安全性指標為約束條件，以降低工程造價為優化目標，這種方法會傾向于方案的經濟性，土石壩的抗液化方案應該綜合考慮安全性和經濟性兩方面的目標，優化方案應在兩者之間尋找平衡。地震導致的土石壩坡面裂縫可能會導致滲漏等現象，從而嚴重危害壩體安全，常用壩體坡面位移梯度作為土石壩安全評價指標。液化過程中，土體由于抗剪強度的降低，常常出現側向大位移，而液化結束后，可液化土體由于震動變得密實，往往出現大沉降［20］，壩頂震陷也常常作為評價土石壩抗震性能的重要指標［21］。因此安全性目標與約束條件將從震中與震后兩方面選取。

表1 設計變量

圖2 設計變量

通過前期試算已知，地震過程中，由于上游庫水壓力的作用，上游坡面水平位移的最大值通常出現在上游壓重平臺與坡面接觸點處，而下游面上水平位移分布從坡腳到坡頂呈現逐漸減小的趨勢。這里將安全性指標選取為：地震過程中上游面平臺接觸點與壩體坡腳及坡頂水平位移差值之和的最大值G1；地震過程中下游面壩體坡腳與坡頂水平位移最大差值G2；壩頂中心震后永久沉降Z。

此外，出于減少工程成本的考慮，設置壓重平臺與地基開挖置換的斷面總面積A為經濟性指標，且考慮到筑壩的成本高于地基處理的成本，將其體現在該指標中：

式中：A、B、C、D、E點分別為上游坡腳、上游坡頂、下游坡頂、下游坡腳和上游面平臺與坡面接觸點；XA、XB、XC、XD、XE分別為A—E點的地震過程中水平位移；s1、s2、s3分別為壩前壓重平臺、壩后壓重平臺和地基置換斷面的面積；λ為壓重平臺與地基處理的每平方米成本單價之比，需根據工程實際確定。

優化目標即為使4個指標G1、G2、Z、A最小化，在優化過程中，根據不同目標的重要性，對每個目標定義權重或比例因子。在設置約束條件時，除根據工程需要設置指標G1、G2、Z和A的約束條件，首先在靜力分析的基礎上，設置壩體靜動力抗滑穩定安全系數的約束條件；其次在動力分析中為了同時實現對地基液化的控制，除計算整體地基液化面積及深度以外，在壩體下方設置兩個監測點，監測并計算整場地震過程中該點液化程度并取平均值作為該點平均液化度，對此設置約束條件。此外，為了控制震后壩體永久變形，對上下游坡腳坡頂的水平絕對位移均設置約束條件。

（3）將上述優化設計的數學模型表示為：

式中：［Xi］，i=a，b，c，d分別為各點水平震后永久位移限制條件；［Ks］［Kd］分別為壩體靜力抗滑穩定安全系數限制條件及擬靜力法計算的壩體動力抗滑穩定安全系數；［tk］，k=1，2為兩個壩基監測點的地震過程中平均液化度限制條件，分別位于上下游坡腳下方壩基深2 m處；［Al］［Hl］分別為地震過程中壩基最大液化面積及液化深度。

4 模型求解方法

4.1 存檔微遺傳算法由于遺傳算法具有自適應性，是一種采用群體搜索策略的全局概率優化算法，它可以不了解問題的內在性質，而直接操作結構對象，目前在多目標優化問題中被廣泛應用［22］。作為一種非歸一化的多目標優化遺傳算法，存檔微遺傳算法（AMGA）非常適合高度非線性、非連續、非凸及高度約束的搜索空間，適用于有多個局部最優的函數。由于本文中模型設計變量眾多，變量與目標間關系復雜，非線性程度較高，因此選用該算法進行優化。

在多目標優化問題中，采用支配關系評估不同方案結果的優劣。在AMGA算法中，對每個目標參數分別進行變異和交叉標準遺傳操作，選擇過程使用三層機制。第一層適應度是根據解在種群中的支配水平來分配的，主要通過設置存檔來進行。第二層適應度的根據為目標方案對算法搜索歷史的貢獻程度，第三層適應度以種群的多樣性為主要考慮因素。因此，AMGA 算法在進化過程之外還需要設立一個存檔，使得檔案中的個體在進化過程中始終保持非支配的地位［23］。

4.2 集成FLAC3D的土石壩抗液化措施優化流程在整個優化流程中，需將FLAC3D作為一個子程序調用，故AMGA優化算法與FLAC3D的串接是關鍵問題。整個流程通過VC++控制實現，并編譯成exe文件，具體思路如下：

（1）在優化控件中設定AMGA算法參數、設計變量及變量空間、約束條件、權重與比例因子、目標函數等；該控件在每一代的進化中輸出設計變量數值，存入txt文件；

（2）調用前處理exe讀取存放設計變量的txt文件，生成參數化建模的命令流txt文件；

（3）利用創建子進程的方式在VC++中調用FLAC3D.exe，并控制其依次讀取命令流文件、自編的液化指示器dll文件以及地震波文件；

（4）FLAC3D進行靜動力計算。在靜力計算的基礎上，利用強度折減法進行安全系數的計算；動力計算中，在壩體上設置a—e五個壩體監測點，監測地震過程中的水平位移值，并在每一時步計算壩基液化范圍、深度及壩基液化監測點的超孔壓比；

（5）計算結束后，利用hist命令提取上述監測值，并讀取安全系數、上下游坡頂坡腳的震后永久水平變形、壩頂震陷量，存入相應結果txt文件；

（6）結束FLAC3D子進程，在VC++主進程中，調用后處理exe，讀取（5）中結果txt 文件，計算平均液化度，提取各約束指標值，同時計算目標指標G1、G2、Z和A的數值，存入存放目標函數及約束指標的txt文件；

（7）將（6）中目標函數及約束指標輸入到優化控件中，判斷是否滿足各約束條件，進行下一代計算。重復步驟（2）—（7），直至算法收斂。上述流程如圖3所示。

圖3 集成優化流程

5 工程實例分析

5.1 數值模型建立某水庫工程是天津市南水北調中線配套工程的重要組成部分［24］。該工程的圍壩型式為復合土工膜防滲斜墻土石壩，總壩長6570 m。選取某典型斷面進行計算，斷面壩高8.6 m，壩頂寬度8.0 m，上游臨水面壩坡坡度為1∶3.0，壩高6.5 m處設置了2.0 m寬的馬道，下游背水側壩坡坡度為1∶3.0。由于圍壩長度方向的尺寸遠大于橫斷面尺寸，把該模型按平面應變問題處理，整體厚度方向為1 m。計算中采用正常運行水位水深為6.0 m。土石壩模型地基截取長度一般取0.3 ～0.5 倍的壩體順河向底部長度［25］，地基深度取30.0 m，模型總長度為173.5 m。上游坡腳開挖錨固槽進行壩基防滲處理，槽深約2.5 m，槽底寬度4.0 m，上游臨水側壩坡設置約0.5 m的鋼筋混凝土板護砌、墊層及復合土工膜，建模過程中將其簡化為整體護砌。

該工程壩基15 m范圍內存在兩層大范圍的液化土層，表層粉土分布連續，第二層液化土層埋藏較淺且上下游貫通，壩基液化問題較為突出。為減輕地震液化的影響，該工程采用挖除置換表層可液土體結合壓重平臺的處理方案，壩后壓重平臺坡面坡度為1∶1，高3.0 m，寬12.0 m，壩前壓重平臺坡面坡度為1∶1.17，高2.50 m，寬15.0 m，挖除表層2.1 m 的可液化粉土。所建立的模型如圖4 所示，模型中材料參數見表2。

圖4 工程實例模型示意圖

表2 材料參數

工程抗震設防烈度為Ⅶ度，庫區地震峰值加速度為0.15g。由于該工程的抗震設防烈度低于Ⅷ度，因此數值模擬計算中僅考慮水平地震作用。文中輸入的地震波（見圖5）為濾波及基線調整后的某南北水平向的天津波（1976年）。

圖5 輸入地震波

5.2 靜動力數值計算計算中編寫了Fish 函數進行參數化建模。讀取優化設計變量后，判斷壩前壓重平臺高度和壩后壓重平臺高度相關關系，以此決定不同的壩體分塊劃分方法。該方法可使相鄰塊體接觸面長度、網格數量和比率相等，以滿足接觸面連續性的要求，同時控制整體模型網格尺寸相似并滿足動力計算的要求。該參數化建模流程僅需輸入設計參數即可實現，便于在優化過程中調用。

土石壩模型的數值計算包括靜力和動力計算兩部分。靜力計算包括壩體填筑及蓄水的模擬，其目的主要是為確定動力計算初始的應力場和滲流場。模擬得到大壩施工填筑完成后的最終沉降為57.96 cm。在圍壩壩段的實測數據中，18個壩體斷面施工期和沉降期的總沉降量平均值為52.83 cm，仿真結果與實測結果是一致的，表明了所建立模型及邊界約束條件的合理性。

動力計算流程包括在震前初始狀態的基礎上，設置自由場邊界及阻尼［26］，輸入地震波，并完成流固耦合計算。其中，采用Finn 模型模擬在動力作用下可液化土體的孔壓積累直至土體液化的過程，該模型假定動力作用中孔壓的上升與塑性體積應變增量相關，增加了動孔壓的上升計算，孔壓上升到一定程度后，土體發生液化。采用有效應力法進行一般應力條件下可液化土體的液化判定。為了便于衡量土體液化的程度，在數值計算中引入超孔壓比的概念。將超孔壓比定義為［27］：

根據上述原理，利用Fish語言編寫了超靜孔隙水壓力和超孔壓比的監測程序。利用自定義單元額外變量在每個計算步中計算各點的超靜孔隙水壓力和超孔壓比，進而通過數據處理得到每個點的液化時長。利用C++語言編寫液化指示器。設置臨界超孔壓比為0.6，對于可液化土層，計算過程中當單元超孔壓比超過臨界超孔壓比，單元狀態將顯示為“液化中”，記錄每一時刻的液化中單元總面積及最大深度。最后將該C++項目文件編譯成dll文件，以便計算時在命令流中直接調用。

采用初步抗液化方案后，地震過程中上游坡腳殘余水平變形為29.88 cm，下游坡腳殘余水平變形為43.79 cm，液化土層在地震強度較大時出現大范圍液化，引起較大的壩體變形。因此，需要開展土石壩抗液化措施優化研究。

5.3 抗液化措施多目標優化分析結合初步抗液化方案設置設計變量空間，見表3，表中坐標值均為相對于模型零點的數值；地面高程為30 m，λ取1.5。優化中需對目標定義合適的權重與比例因子，將各目標數值統一在相似水平［28］。

式中：Xi和為各目標調整前后的數值；Wi為各目標對應的權重因子；SFi為各目標對應的比例因子。通過前期試算，確定權重與比例因子見表4。根據初始方案計算結果給定目標、液化、安全系數與殘余變形的約束條件。

表3 設計變量約束及優化結果

表4 目標函數參數及優化結果

采用AMGA 算法對多目標優化數學模型進行求解。種群規模設置為40，進化到第129代算法達到收斂。經過三次計算，AMGA算法平均尋優率為40.63%，相較于常用的第二代非劣排序遺傳算法（NSGA-Ⅱ）的三次平均尋優率85.49%，其收斂性相對較差，但其探索能力強，搜索并不集中局部最優解，得到的優化方案效果也更好，較為適合本文的設計變量較多及多峰情況的數學模型搜索。由算法搜索得到的解集的目標函數分布如圖6所示，并根據約束條件及目標函數得出Pareto前沿。由圖可知，AMGA 算法前進性較強，在Pareto 前沿處探索較為密集，且Pareto 前沿呈現出兩兩制約關系，難以同時達到最優，需要根據工程實際在Pareto前沿中選擇優化方案。

舍去不滿足各約束及取值較為極端的方案，得到的最終方案見表4。從選取的優化指標來看，壩頂震陷有了較大改善，震陷量減少了13.22 cm，優化幅度達到21.03%；壩體上下游面的水平變形梯度較初始方案均有改善，分別減少了16.59%和24.33%，壩體安全性得到提高；抗液化措施成本優化明顯，降低了18.53%，經濟性得到了改善。優化前后的G1、G2及Z的時程如圖7所示，在地震過程中，G1、G2及Z的數值均是不斷發展的，下游坡面水平位移梯度較大，震后為壩體坡面變形梯度最大的時刻。采取優化方案后，上下游坡面壩體位移梯度在整個地震過程中均得到了控制。從圖8中可以看出，壩體深層滑動趨勢向后移動，并未貫穿至上游坡面，且下游壩體下方水平位移明顯減少，未形成完整滑動面，壩體變形及整體抗滑穩定得到改善。震后沉降主要集中在上游壩體，整體沉降量明顯減少。

各點殘余變形及平均液化度見表5。其中，采取初始方案后震后上游坡腳出現了79.85 cm的較大沉降，下游坡腳出現了43.79 cm的較大水平變形，結合圖8，壩體出現向下游的深層滑動。采取優化方案后，壩體下方深層滑動趨勢有所降低，壩體下方偏下游的位置位移控制較為明顯。優化方案滿足了殘余變形、平均液化度、最大液化面積、液化深度及安全系數等各項約束條件，并有了一定程度的改善，其中對水平殘余變形的改善最為明顯。

圖6 算法搜索結果及Pareto前沿分布

圖7 優化前后G1、G2、Z時程圖

圖8 優化前后殘余位移云圖（單位：m）

表5 各約束條件優化結果

圖9顯示下游液化監測點的超孔壓比明顯高于上游監測點，優化前后下游監測點的最大超孔壓比為0.983和0.724，震后超孔壓比分別為0.484和0.431，存在液化的時刻，上游監測點最大超孔壓比分別為0.340和0.273，震后超孔壓比分別為0.086和0.015，地震過程中并未發生液化。下游的液化程度更大，與圖8中下游水平位移較大相對應。采用優化方案后，兩個監測點在地震過程中超孔壓比的發展均有了一定程度的改善，平均超孔壓比分別由0.115和0.578降至0.095和0.430，壩基最大液化面積也從1618.59 m2降低到了1560.32 m2，表明了優化方案對于壩基液化問題的有效改善。

圖9 優化前后超孔壓比時程對比

6 結論

針對目前工程中主要依靠經驗進行土石壩抗液化措施設計的現狀，提出了一種考慮多目標、多約束條件下的土石壩抗液化措施數值建模與優化分析方法。該方法以壩體坡面水平位移梯度、壩頂震陷和抗液化措施成本最小化為優化目標，以靜動力抗滑穩定安全系數、地基最大液化范圍、深度、典型點平均液化度和壩體永久水平變形為約束條件，建立了相應的數學模型，運用完全非線性動力算法和多目標遺傳算法對土石壩的壓重平臺和壩基置換參數進行優化設計，實現了綜合考慮土石壩抗液化措施的安全性能與經濟效益。該方法應用于工程實例，對工程中壓重平臺及地基置換的9個控制性設計參數進行優化，最終方案壩體水平位移指標G1較初始方案降低了16.59%，水平位移指標G2較初始方案降低了24.33%，壩頂沉降指標Z較初始方案降低了21.03%，節約成本18.53%，壩基液化與震后殘余變形均得到了一定改善。

傳統的土石壩抗液化方案設計方法存在主要根據經驗設計、不能定量化分析、易造成浪費等問題，所提出的綜合考慮安全性和經濟性的土石壩的抗液化措施優化分析方法，初步實現了對土石壩抗液化措施的定量化、精細化設計與優化，并可進一步推廣應用于地基中存在可液化土的各類建筑工程抗液化措施設計與優化中，具有重要的工程實踐價值。

需要指出的是，所建立的優化模型中目標函數及約束條件的權重可根據實際工程需求來進行選取和調整；現有數值模擬中一般采用超孔壓比或超孔隙水壓力作為液化的標準進行抗液化措施的效果評價，但在實際工程中液化的影響因素眾多，不同的液化位置及液化程度對工程安全的影響各不相同，因此該指標與工程安全的具體關系仍有待深入探討；由于優化模型的目標函數與約束條件眾多，導致優化中對算法需要提出更高的要求。因此，后續研究中需要完善地基液化狀況的評估方法并提出有針對性的抗液化措施，同時在算法的探索性及并行計算方面需要進一步研究。