基于改進(jìn)ACMPSO并行算法的土石壩本構(gòu)參數(shù)反演

陳家琦, 岑威鈞, 李鄧軍, 潘正陽(yáng)

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

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

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

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

1 ACMPSO算法及其改進(jìn)

1.1 ACMPSO算法

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

1.1.1混沌序列

混沌序列是指由某一初值和映射方程得到的一系列隨機(jī)數(shù)，因其具有隨機(jī)性、邏輯確定性、強(qiáng)敏感性、遍歷性和不可預(yù)測(cè)性等特點(diǎn)，已被廣泛運(yùn)用于密碼學(xué)當(dāng)中。常用的混沌映射方程有Double-Bottom映射、Henon映射和Lorenz映射等[11]。

Yang等[12]對(duì)Double-Bottom映射進(jìn)行了改進(jìn)，構(gòu)造出對(duì)初值敏感性更強(qiáng)的混沌映射，其表達(dá)式為

(1)

式中：rp為混沌序列第p個(gè)數(shù)。該映射在0

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

1.1.2 全局極值自適應(yīng)變異策略

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

gbestj=gbestj(1+δω)

(2)

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

1.2 改進(jìn)的AMCPSO算法

1.2.1 全局極值自適應(yīng)變異的改進(jìn)

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

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

(3)

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

為了引導(dǎo)全局極值朝著正確的方向變異，采用如下改進(jìn)的變異公式：

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

(4)

1.2.2粒子位置自適應(yīng)變異

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

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

(5)

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

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

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

(6)

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

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

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

(7)

1.2.3并行計(jì)算

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

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

1.3 算法驗(yàn)證

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

(8)

圖1 面板壩算例有限元網(wǎng)格模型及測(cè)點(diǎn)布置

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

圖2 不同算法迭代過(guò)程曲線對(duì)比

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

2 實(shí)例驗(yàn)證

2.1 工程概況

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

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

根據(jù)大壩實(shí)際填筑方式，劃分大壩三維有限元計(jì)算網(wǎng)格，大壩網(wǎng)格共有7 026個(gè)節(jié)點(diǎn)，6 610個(gè)單元，壩體實(shí)際測(cè)點(diǎn)分布見(jiàn)圖4。

圖4 壩體實(shí)際測(cè)點(diǎn)布置

2.2 反演參數(shù)選取與設(shè)置

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

表1 筑壩料E-B模型參數(shù)

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

2.3 反演結(jié)果分析

在相同條件下進(jìn)行3次獨(dú)立的參數(shù)反演計(jì)算，得到的迭代過(guò)程曲線如圖5所示。由圖5可知，3次反演計(jì)算的迭代過(guò)程比較相近，迭代計(jì)算5次后，目標(biāo)函數(shù)值降幅可以達(dá)到7%，占總降幅的一半以上，再次說(shuō)明改進(jìn)的ACMPSO并行算法收斂速度很快，最終經(jīng)過(guò)50次迭代計(jì)算，目標(biāo)函數(shù)的平均值低于5%，滿(mǎn)足工程誤差限。

圖5 3次參數(shù)反演計(jì)算的迭代過(guò)程曲線

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

表2 3次參數(shù)反演結(jié)果的平均值

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

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

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

3 結(jié) 論

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

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