陳家琦, 岑威鈞, 李鄧軍, 潘正陽

(河海大學水利水電學院，江蘇 南京 210098)

為了評估土石壩的安全性，預測大壩的應力變形，需要獲得土石料真實力學參數，為此常用現(xiàn)場試驗和室內三軸試驗[1-2]進行測試。但是，基于這類方法得到的參數計算得到的大壩應力變形與原位監(jiān)測信息常有出入。為了準確計算和預測土石壩的應力變形，基于原位監(jiān)測信息的參數反演成為獲得壩體材料真實力學參數的有效方法[3]。

在巖土材料參數反演方面，不少學者引入了各種智能優(yōu)化算法。張美英[4]采用免疫遺傳算法對三峽茅坪溪瀝青心墻壩進行了反分析，預測了大壩實際的運行性態(tài)；張丙印等[5]構造了基于人工神經網絡和演化算法的位移反演分析方法，并采用該方法對三峽茅坪溪防護土石壩的變形進行了反演分析，驗證了反演分析方法的有效性；賈陸鋒等[6]將單純形算法和模擬退火算法結合為一種混合優(yōu)化算法，為反演地下工程圍巖力學參數提供了有效途徑；杜好[7]采用改進的單群體自適應微粒群算法對某面板堆石壩的4個鄧肯E-B材料參數進行了反演分析，得到了比較合理的參數結果。此外，蟻群算法等其他優(yōu)化方法在參數反演中也有不少應用[8-9]。

高維復雜巖土工程的材料參數反演問題通常存在最優(yōu)解不唯一且計算量龐大的問題，常用優(yōu)化算法容易陷入局部最優(yōu)且計算效率十分低下。針對這類問題，本文對自適應混沌變異粒子群算法(adaptive chaotic mutation particle swarm optimization, ACMPSO)進行了改進，并采用改進的ACMPSO并行算法對某實際工程的面板堆石壩的4個材料分區(qū)共20個鄧肯E-B材料參數進行了反演分析。

1 ACMPSO算法及其改進

1.1 ACMPSO算法

Kennedy等[10]根據鳥群覓食的行為提出了粒子群基本算法(particle swarm optimization,PSO)，該算法通過粒子種群間的相互指導和學習的方式來迭代尋找優(yōu)化問題的全局最優(yōu)解。ACMPSO是在PSO算法的基礎上，引入了混沌序列和全局極值自適應變異策略，以改善PSO算法遍歷性差、易于陷入局部最優(yōu)等問題。

1.1.1混沌序列

混沌序列是指由某一初值和映射方程得到的一系列隨機數，因其具有隨機性、邏輯確定性、強敏感性、遍歷性和不可預測性等特點，已被廣泛運用于密碼學當中。常用的混沌映射方程有Double-Bottom映射、Henon映射和Lorenz映射等[11]。

Yang等[12]對Double-Bottom映射進行了改進，構造出對初值敏感性更強的混沌映射，其表達式為

(1)

式中：rp為混沌序列第p個數。該映射在0

本文采用式(1)生成多條隨機序列，用于PSO算法中粒子位置、速度的初始化，并代替速度更新公式中的隨機數，保證種群在可行域內均勻分布，提高算法的搜索能力。

1.1.2 全局極值自適應變異策略

為了改善PSO算法易于陷入局部最優(yōu)的問題，呂振肅等[13]引入了全局極值自適應變異策略，采用粒子群體適應度方差σ2反映粒子群體的收斂程度，當σ2小于某一閾值時，按某一概率對全局極值進行變異，全局極值第j維度的變異公式為

gbestj=gbestj(1+δω)

(2)

式中：gbestj為第j維度的全局極值；δ為變異系數，一般取0.5；ω為服從Gauss(0,1)分布的隨機數。

1.2 改進的AMCPSO算法

1.2.1 全局極值自適應變異的改進

ACMPSO算法中，全局極值執(zhí)行條件是粒子群體適應度方差σ2小于某一閾值，但是對于一些維度高、計算量大的材料參數反演問題，閾值敏感性很強，可能粒子群體還沒有趨于收斂就會陷入局部最優(yōu)，而這類問題的敏感性分析會耗費大量時間，這是不可取的。另一方面，變異公式(式(2))采用的是靜態(tài)變異系數0.5，因此粒子所有維度具有統(tǒng)一的變異幅度和方向，對于某些粒子維度較高的復雜問題，統(tǒng)一的變異幅度和方向不一定適合于全局極值的所有維度，執(zhí)行變異策略后全局極值的某些維度仍有可能處于邊界附近，導致算法不能跳出局部最優(yōu)或進入另一個局部最優(yōu)。

針對上述問題，本文采用全局極值變異率Pm作為策略執(zhí)行條件執(zhí)行變異操作，雖然增加了算法的計算時間，但可以全程提高算法的全局搜索能力。在計算過程中，由于粒子迭代的位置是隨機的，很可能一部分粒子會停留在可行域的邊界附近，一部分粒子停留在可行域的中心，那么對于靠近可行域邊界的粒子，變異幅度應該較大；靠近可行域中心的粒子，變異幅度應該較小。因此，采用動態(tài)變異系數要比靜態(tài)變異系數更具有普適性，使粒子的各維度都得到最合適的變異幅度。全局極值第j維度的動態(tài)變異系數δbestj表達式為

(3)

式中：uj、lj分別為可行域第j維度的上界和下界；τ為最小擾動系數，0<τ<1；μ為控制系數，1<μ<2。

為了引導全局極值朝著正確的方向變異，采用如下改進的變異公式：

gbestj=gbestj[1+ψ(gbestj)δbestj|ω|]

(4)

1.2.2粒子位置自適應變異

為了提高PSO算法前期的局部搜索能力，加快收斂速度，借鑒遺傳算法的變異思想[14]，引入動態(tài)變異函數控制粒子位置變異，在算法前期，變異率較高，之后隨著迭代次數的增加逐步降低。動態(tài)變異函數V(t)為

V(t)=V0[1-(t/tmax)α]β

(5)

式中：V0為初始變異率；t為當前迭代次數；tmax為最大迭代次數；α、β為控制參數，一般取α=1.7，β=10。

在動態(tài)變異函數控制下，算法初期會篩選較多的粒子進行變異，擴大粒子在可行域內搜索范圍，加快了算法的收斂速度；在算法中后期，變異的粒子數量逐漸降低，防止粒子異常波動影響算法的收斂能力。

變異粒子i的第j維度的動態(tài)變異系數δij為

(6)

式中：zij為變異粒子i的第j維度的位置。

根據變異控制函數選擇出需要變異的粒子后，對變異粒子的位置進行變異，變異粒子i的第j維度的粒子位置變異公式為

zij=zij[1+ψ(zij)δij|ω|]

(7)

1.2.3并行計算

采用Fortran語言編寫PSO、ACMPSO和改進的ACMPSO算法，將生成的各粒子參數組寫為材料文件，利用批處理文件間接調用ADINA軟件讀取材料文件，進行有限元計算和后處理[15]。由于對每個粒子進行的有限元計算是相互獨立的，因此可以在反演程序中插入并行指令進行并行計算。

共享內存編程應用程序接口的工業(yè)標準協(xié)議(OpenMP)為一種多線程的程序設計指令，可以在已編寫的Fortran程序段內以標識符的方式添加并行指令，用于指示Fortran程序開始和結束多線程運行的時間和地址[16-17]。為了防止多線程運行可能引起的文件讀寫沖突，反演算法選用了OpenMP Fortran API中具有獨立計算區(qū)域的Sections并行指令。

1.3 算法驗證

為了驗證改進的ACMPSO并行算法的合理性，以某面板堆石壩為例，布置35個模擬測點，有限元網格模型及測點布置如圖1所示，堆石料本構模型采用鄧肯E-B模型[18-19]，對壩體主堆石和次堆石的切線彈性模量參數K、n以及切線體積模量參數Kb進行反演計算。算法驗證中，參數浮動范圍為準確值的±50%，模擬實測值為有限元正算得到的竣工期壩體沉降，目標函數公式為

(8)

圖1 面板壩算例有限元網格模型及測點布置

分別采用PSO算法、ACMPSO算法和改進的ACMPSO并行算法對面板壩算例進行參數反演，各算法參數均設置為：種群數量M=4，初始變異率V0=0.5，全局極值變異率Pm=0.5，慣性因子w最大值、最小值分別為0.9和0.1，最小擾動系數τ=0.1，控制系數μ=1.5，固定初始粒子位置和參數浮動范圍，3種算法的目標函數迭代過程曲線如圖2所示。

圖2 不同算法迭代過程曲線對比

由圖2可以看出，PSO算法的收斂速度很慢，在迭代15次之后陷入了局部最優(yōu)，難以跳出；ACMPSO算法的收斂速度與PSO算法大致相同，但是在第6次迭代時跳出了局部最優(yōu)，目標函數值大幅下降，隨后算法再次陷入局部最優(yōu)，說明ACMPSO算法不能有效避免算法陷入局部最優(yōu)，不適用于高維度的材料參數反演；改進的ACMPSO并行算法收斂速度很快，迭代10次后，目標函數值降幅已達總降幅的90%以上，中后期算法沒有陷入局部最優(yōu)，目標函數值呈下降趨勢，最終收斂至0.75%，說明改進的ACMPSO并行算法前期局部搜索能力強、后期全局搜索能力強，在高維度的材料參數反演中能夠有效地避免算法陷入局部最優(yōu)。

2 實例驗證

2.1 工程概況

某面板堆石壩壩頂高程757.0 m，最大壩高143.0 m，壩頂長度460.0 m。壩體上游壩坡1∶1.5，下游壩坡1∶1.3。壩0+261.0 m斷面處布置有3層水管式沉降儀，高程分別為650.0 m、687.8 m和721.7 m，共18個測點。通過對測點實測值的分析整理，發(fā)現(xiàn)ES3-3和ES3-7沉降儀損壞，ES1-3、ES1-4、ES3-6和ES3-8沉降儀測值異常，反演分析時剔除，其余12個測點作為本文反演分析參考點。沉降儀布置及實測結果如圖3所示。

圖3 壩0+261.00 m斷面沉降儀布置及實測結果(高程單位：m；沉降單位：cm)

根據大壩實際填筑方式，劃分大壩三維有限元計算網格，大壩網格共有7 026個節(jié)點，6 610個單元，壩體實際測點分布見圖4。

圖4 壩體實際測點布置

2.2 反演參數選取與設置

壩體填筑料主要包括墊層料、過渡層、主堆石和次堆石。E-B模型共有9個參數，其中初始內摩擦角φ、內摩擦角增量Δφ和土體黏聚力c一般可由試驗確定，卸載彈性模量參數Kur的敏感性很低，可不用反演，因此將切線彈性模量參數K、n，破壞比Rf和切線體積模量參數Kb、m作為待反演參數，參考現(xiàn)場試驗數據和類似的工程經驗，擬定了壩體材料待反演參數的取值范圍以及其余參數的取值如表1所示。

表1 筑壩料E-B模型參數

由于反演分析共有20個待反演參數，且部分待反演參數的敏感性較低，需要適當擴大迭代次數和粒子種群數，算法參數設置為：粒子種群數M=8，初始變異率V0=0.8，全局極值變異率Pm=0.7，最小擾動系數τ=0.4，最大并行線程數為4，其余參數均與1.3節(jié)一致。為了確保反演結果的可靠性和穩(wěn)定性，一般需要在相同條件下進行多次(不少于3次)獨立的反演計算，計算的次數視反演計算量的大小而定，最終反演結果應取為多次獨立反演計算結果的平均值。

2.3 反演結果分析

在相同條件下進行3次獨立的參數反演計算，得到的迭代過程曲線如圖5所示。由圖5可知，3次反演計算的迭代過程比較相近，迭代計算5次后，目標函數值降幅可以達到7%，占總降幅的一半以上，再次說明改進的ACMPSO并行算法收斂速度很快，最終經過50次迭代計算，目標函數的平均值低于5%，滿足工程誤差限。

圖5 3次參數反演計算的迭代過程曲線

表2給出了各材料參數反演結果的平均值，可以看出反演計算得到的參數均比較合理，符合一般規(guī)律，且3次反演的參數結果沒有出現(xiàn)異常波動，說明改進的ACMPSO并行算法穩(wěn)定性好，能夠得到合理的反演結果。

表2 3次參數反演結果的平均值

圖6給出了壩0+261.00 m監(jiān)測斷面上有效測點的實測沉降和反演沉降(3次反演計算沉降的平均值)的分布對比，可以看出，壩體實測沉降與反演計算的沉降吻合度很高，計算結果符合面板堆石壩的沉降分布規(guī)律，測點沉降的平均相對誤差為4.6%，最大相對誤差為10.7%，位于ES3-1測點，說明改進的ACMPSO并行算法具有較高的搜索精度，能夠處理高維度、計算量龐大的參數反演問題。

圖6 壩0+261.00 m斷面竣工期實測沉降與反演沉降分布對比(高程單位：m；沉降單位：cm)

3次反演計算的平均時長為9.26 h，而單線程的反演計算則需要28.43 h，并行計算的平均加速比為3.07，并行效率為77%。若提高最大并行線程數和計算機性能，加速比可以繼續(xù)提高，但并行效率會降低。

3 結 論

a.改進的ACMPSO并行算法融合了新的混沌映射和全局極值變異策略，引入了粒子位置變異策略和并行算法。與PSO算法和ACMPSO算法相比，改進的ACMPSO并行算法具有收斂速度快，能有效避免算法陷入局部最優(yōu)等特點。

b.將改進的ACMPSO并行算法應用于某面板堆石壩筑壩材料的參數反演中，3次獨立反演計算的迭代過程比較相近，監(jiān)測點沉降的平均相對誤差為4.6%；利用反演值得到的壩體變形與原位監(jiān)測量匹配度高，合理反映了壩體實際的變形情況；并行計算的平均加速比為3.07，說明改進算法穩(wěn)定性好，搜索精度高，能夠大幅提升計算效率，為處理高維度、計算量龐大的復雜參數反演問題提供了一個可行的方法。