加速遺傳算法在反演計(jì)算土層的剪切模量、阻尼比模型參數(shù)值中的應(yīng)用研究

李大虎， 何 強(qiáng)， 顧勤平

(1.四川省地震局工程地震研究院，四川成都 610041;2.江蘇省地震局地震工程研究院，江蘇南京 210014)

加速遺傳算法在反演計(jì)算土層的剪切模量、阻尼比模型參數(shù)值中的應(yīng)用研究

李大虎1， 何 強(qiáng)1， 顧勤平2

(1.四川省地震局工程地震研究院，四川成都 610041;2.江蘇省地震局地震工程研究院，江蘇南京 210014)

針對成都地區(qū)土層特有的沉積環(huán)境及物理力學(xué)性質(zhì)的差異性，采用Martin-Davidenkov模型和阻尼比經(jīng)驗(yàn)公式，使用基于實(shí)碼、并引入了更新概率(Pu)的加速遺傳算法(AGA)，編制了相應(yīng)的C++程序，反演計(jì)算出成都地區(qū)各類沉積土的動剪切模量比G/Gmax和阻尼比λ與剪應(yīng)變幅值γ的關(guān)系曲線模型的參數(shù)值。結(jié)果表明，RAGA收斂速度快，計(jì)算精度高，是一種既可以較大概率搜索全局最優(yōu)解，又能進(jìn)行局部細(xì)致搜索的非線性優(yōu)化方法。

成都地區(qū);加速遺傳算法;反演;動剪切模量比;阻尼比

工程場地地震安全性評價(jià)，最重要的一項(xiàng)工作是確定設(shè)計(jì)地震動參數(shù)，而這些參數(shù)是在考慮了局部場地影響的條件下，通過土層地震反應(yīng)分析計(jì)算出來的。在進(jìn)行土層地震反應(yīng)分析的過程中，除了要用到土層的剪切波速值以外，還需要用到表示土的動力學(xué)特性曲線，即G/Gmax-γ曲線和λ-γ曲線。而位于龍門山造山帶前緣的成都新生代斷陷盆地，第四紀(jì)以來一直處于穩(wěn)定的沉降狀態(tài)，沉積了厚約數(shù)百米沖洪積砂卵礫石層，因此造就了成都地區(qū)特有的土層結(jié)構(gòu)特征?？紤]到這種特殊沉積環(huán)境，在對該地區(qū)有代表性的大量土樣進(jìn)行系統(tǒng)而深入的試驗(yàn)研究基礎(chǔ)上，對其關(guān)系曲線模型參數(shù)進(jìn)行了反演計(jì)算;然而，在對大量試驗(yàn)數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析的時(shí)候，針對以往二進(jìn)制編碼的標(biāo)準(zhǔn)遺傳算法(SGA)尋優(yōu)效果明顯依賴于模型參數(shù)的初始變化區(qū)間的大小，并且局部搜索能力較差、可能會出現(xiàn)過早收斂等問題，作者提出了一種改進(jìn)的算法，即基于實(shí)碼、并引入了更新概率(Pu)的加速遺傳算法(AGA)，該算法是在不斷收縮搜索范圍，直到目標(biāo)函數(shù)值小于預(yù)先設(shè)定值或加速達(dá)預(yù)定次數(shù)的新方法研究，并且通過對土的剪切模量、阻尼比隨剪應(yīng)變幅值變化的關(guān)系曲線進(jìn)行反演計(jì)算，最終給出了成都地區(qū)各類沉積土G/Gmax-γ和λ-γ平均擬合曲線及其參數(shù)的推薦值。

1 基于實(shí)碼的加速遺傳算法原理和實(shí)現(xiàn)流程

1.1 基于實(shí)數(shù)編碼的加速遺傳算法原理

設(shè)一般非線形模型參數(shù)優(yōu)化問題(即目標(biāo)函數(shù))為:

式中，xl為p個(gè)待優(yōu)化變量;Xi，Yi為m對模型輸入、輸出的觀測數(shù)據(jù)，q為任意實(shí)常數(shù)。

本文以計(jì)算值和觀測值之差的平方和最小為標(biāo)準(zhǔn)來求取模型的最優(yōu)參數(shù)，使得目標(biāo)函數(shù)值最小，如式(1)所示，并根據(jù)目標(biāo)函數(shù)的最小值保存最好的模型(周明等，1999)。

1.2 實(shí)現(xiàn)的流程

為提高剪切模量和阻尼比模型參數(shù)的計(jì)算速度，本文應(yīng)用改進(jìn)的遺傳算法，即參數(shù)編碼采用了實(shí)數(shù)(浮點(diǎn)數(shù))編碼，并在運(yùn)算過程中引入了更新概率和加速算子，形成了加速遺傳算法，因此AGA的循環(huán)可逐步調(diào)整、縮小優(yōu)化變量的尋優(yōu)區(qū)間，解的精度隨著循環(huán)次數(shù)的增加可望逐步提高(王福剛等，2002;周雙喜等，1996)。

隨機(jī)產(chǎn)生一個(gè)初始種群，首先執(zhí)行遺傳算法的基本操作(選擇、交叉、變異)生成新的群體，同時(shí)選記錄適應(yīng)度較高的個(gè)體;然后根據(jù)設(shè)計(jì)好的局部搜索概率在對中間群體中選取部分個(gè)體進(jìn)行局部搜索計(jì)算，并將所得的解替換原來的適應(yīng)度，跟未進(jìn)入局部搜索的個(gè)體和記錄的個(gè)體共同組成含有新的種群再次進(jìn)行更新計(jì)算;有限次運(yùn)算后進(jìn)行加速遺傳，縮小優(yōu)秀個(gè)體選擇的區(qū)間(分別將一定次數(shù)演化迭代的優(yōu)秀個(gè)體變化區(qū)間作為下一次加速遺傳的變量區(qū)間)，這樣演化迭代與加速遺傳的反復(fù)交替進(jìn)行可實(shí)現(xiàn)遺傳進(jìn)化逐步向最優(yōu)個(gè)體逼近，并且隨著接近優(yōu)秀個(gè)體，個(gè)體的密度加大，這樣可在一定程度上減少早熟收斂的機(jī)率。具體的流程詳見圖1(王錦國等，2002;楊曉華等，1998)。

圖1 加速遺傳算法流程圖Fig.1 The flow chart of AGA

1.3 兩個(gè)關(guān)鍵問題

為了計(jì)算得到模型最優(yōu)參數(shù)值，本文在遺傳算法的基礎(chǔ)上做了改進(jìn)，不但在計(jì)算的速度上有了提高，而且引入了更新概率，并通過對比研究設(shè)置參數(shù)。

1.3.1 加速循環(huán)

由于SGA受初始區(qū)間影響很大，如果能隨迭代演化的進(jìn)行而逐步縮小搜索范圍，將大大提高搜索效率，所以在每次加速循環(huán)中，只進(jìn)行一次演化迭代，根據(jù)目標(biāo)函數(shù)值，選出最優(yōu)秀的10%的個(gè)體，并將這些優(yōu)秀個(gè)體的最小值和最大值賦給新搜索區(qū)間兩端，由此進(jìn)行搜索。需要注意的是，兩端的值會先賦給二個(gè)個(gè)體，這樣避免漏掉在區(qū)間中很偏的最優(yōu)解。如果迭代次數(shù)已經(jīng)滿足設(shè)定值，或目標(biāo)函數(shù)小于某一特定值，則轉(zhuǎn)入執(zhí)行接下來的步驟。

1.3.2 更新概率

本文引入更新概率，就是一個(gè)再次篩選的過程，如果子代模型的適應(yīng)度值比較大，那么子代模型總是會被保留，如果父代的適應(yīng)度值比較大，那么就根據(jù)更新概率來替代子代個(gè)體，即控制著父代個(gè)體在子代個(gè)體中的影響，那么如何設(shè)置更新概率呢?本文通過反復(fù)試驗(yàn)，研究成果如(圖2)。

在圖2(A)中，當(dāng)Pu=0時(shí)，迭代100次的時(shí)候目標(biāo)函數(shù)平均值才有了比較明顯的降低，在100次迭代以后目標(biāo)函數(shù)平均值仍起伏不定，時(shí)高時(shí)低很不穩(wěn)定，到迭代200次完成了仍然在增大;在圖2(B)中，當(dāng)Pu=0.5時(shí)，在迭代50次的時(shí)候目標(biāo)函數(shù)平均值就迅速減小為0.3，并且在以后的迭代次數(shù)中有逐漸降低的、比較穩(wěn)定的趨勢，到迭代200次是已經(jīng)變得很小了;在圖2(C)中，當(dāng)Pu=0.9時(shí)，一開始目標(biāo)函數(shù)平均值就迅速下降，在迭代不到50次時(shí)收斂至0.1，在以后的迭代中基本都能保持較穩(wěn)定?？梢姡罡赂怕试O(shè)置較高值能加速收斂，提高反演的效率，故在反復(fù)試算的基礎(chǔ)上，本文將Pu設(shè)置為0.9。

2 實(shí)例計(jì)算

在對現(xiàn)有每個(gè)采樣土層的試驗(yàn)數(shù)據(jù)資料統(tǒng)計(jì)分析的基礎(chǔ)上，采用式(3)的Martin-Davidenkov模型和式(4)的阻尼比經(jīng)驗(yàn)公式(陳國興等，2004)，使用加速遺傳算法程序，分別對G/Gmax-γ和λ-γ模型參數(shù)進(jìn)行了反演計(jì)算，其中，A，B，γ0，λmin，λ0，n為模型參數(shù)。本次計(jì)算中的加速遺傳算法的部分參數(shù)定義如下:初始群體大小300個(gè)，迭代次數(shù)80次，交叉概率0.6，變異概率0.01，更新概率0.9，具體的參數(shù)計(jì)算結(jié)果詳見表1。

圖2 目標(biāo)函數(shù)最優(yōu)值隨迭代次數(shù)的變化曲線Fig.2 Optimal value of the objective function curves with the number of iterations

表1 成都地區(qū)各類沉積土G/Gmax-γ和λ-γ曲線的模型參數(shù)值及目標(biāo)函數(shù)值Tab.1 The model parameter values and objective function value of G/Gmax-γand λ-γ of deposited soil types in Chengdu region

現(xiàn)以中密卵石為例，演示遺傳算法在反演計(jì)算模型參數(shù)過程中應(yīng)用，其中圖3和圖4分別顯示的是反演計(jì)算G/Gmax-γ和λ-γ模型參數(shù)時(shí)的目標(biāo)函數(shù)在迭代過程中的變化情況，目標(biāo)函數(shù)值是根據(jù)剪應(yīng)變和利用遺傳算法迭代出的擬合參數(shù)值求得的新動剪切模量比與原始的動剪切模量比的差的平方，用遺傳算法的目的就是使目標(biāo)函數(shù)值最小，通過圖3和圖4可以看出，在整體上，目標(biāo)函數(shù)曲線是一個(gè)減小的趨勢，即RAGA在每次迭代后目標(biāo)函數(shù)值減小，這是由于加速循環(huán)中選出了部份最佳個(gè)體，并根據(jù)這些個(gè)體的最佳值收縮搜索范圍，使得每次迭代的最佳值總是向真值靠近，故目標(biāo)函數(shù)值一般會減小;圖3為計(jì)算動剪切模量比與剪應(yīng)變關(guān)系模型參數(shù)時(shí)，80次迭代過程中目標(biāo)函數(shù)值的變化，可見整體趨勢上是逐步下降的，最后的目標(biāo)函數(shù)值能達(dá)到0.012 145;圖4則代表計(jì)算阻尼比與動剪切模量比關(guān)系模型參數(shù)時(shí)，80次迭代過程中目標(biāo)函數(shù)值的變化，最后的目標(biāo)函數(shù)值達(dá)到0.004 322。

從圖5和圖6中可以看出，通過加速遺傳算法計(jì)算出來的動剪切模量比、阻尼比與剪應(yīng)變的變化關(guān)系模型曲線與試驗(yàn)統(tǒng)計(jì)值的變化趨勢是一致的(袁曉銘等，2000;陳國興等，1995)。

通過改進(jìn)的加速遺傳算法反演計(jì)算出的各類土層的模型參數(shù)值結(jié)合大量試驗(yàn)數(shù)據(jù)統(tǒng)計(jì)結(jié)果的平均化處理①四川大學(xué)水電學(xué)院.2002.成都地鐵一期工程地基土動力特性試驗(yàn)報(bào)告.②中國科學(xué)院水利部成都山地災(zāi)害與環(huán)境研究所.2008.成都雙流國際機(jī)場T2航站樓地基土動力特性試驗(yàn)報(bào)告.，得到了成都地區(qū)各類沉積土的動剪切模量比和阻尼在剪應(yīng)變1×10－6～1×10－2范圍內(nèi)變化的參考值(表2)。

表2 成都地區(qū)各類沉積土G/Gmax-γ和λ-γ關(guān)系曲線的參考值Tab.2 The value of reference curve of G/Gmax-γand λ-γof deposited soil types in Chengdu region

4 結(jié)論

(1)反演計(jì)算的第一步是正演，只有確定了正演模型，才能與算法結(jié)合形成遺傳反演系統(tǒng)，本文則采用Martin-Davidenkov模型和阻尼比經(jīng)驗(yàn)公式，使用了基于實(shí)碼、并引入了更新概率(Pu)的加速遺傳算法(AGA)，反演計(jì)算出成都地區(qū)各類沉積土的動剪切模量比G/Gmax和阻尼比λ與剪應(yīng)變幅值γ的關(guān)系曲線模型的參數(shù)值，為了方便工程的使用，也給出了各類土動剪切模量比和阻尼比隨剪應(yīng)變幅值變化關(guān)系的參考值。期望通過本文的研究結(jié)果不僅能為加深對成都地區(qū)各類第四系沉積土動力特性的認(rèn)識提供依據(jù)，而且為成都地區(qū)一般工程建設(shè)場地的地震安全性評價(jià)工作提供基礎(chǔ)資料。

(2)針對標(biāo)準(zhǔn)的遺傳算法(SGA)早熟現(xiàn)象(即很快收斂到局部最優(yōu)解，而不是全局最優(yōu)解)，以及在快要接近最優(yōu)解的時(shí)候，在最優(yōu)解附近左右擺動，收斂較慢的問題，在消化和吸收已有成果的基礎(chǔ)上，做了算法上的改進(jìn)研究并順利實(shí)現(xiàn)了加速遺傳計(jì)算。當(dāng)然，反演計(jì)算仍存在著非唯一性，本文中只采取了設(shè)置固定迭代次數(shù)的方法，是為了更好的觀察目標(biāo)函數(shù)值在每一次迭代過程中的變化和趨勢，當(dāng)然也可以設(shè)置一個(gè)固定值，當(dāng)目標(biāo)函數(shù)值滿足了這一條件就停止搜索。但如何評價(jià)反演出來的最后結(jié)果，是需要更多的約束條件的，具體的內(nèi)容筆者今后還會另文敘述。

在論文的撰寫過程中，得到了成都理工大學(xué)地球探測與信息技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室李才明教授的悉心指導(dǎo)和幫助，感謝審稿專家對本文提出的寶貴建議，在此向他們表示衷心的感謝!

陳國興，劉雪珠.2004.南京及其鄰近地區(qū)新近沉積土的動剪切模量和阻尼比試驗(yàn)研究［J］.巖石力學(xué)與工程學(xué)報(bào)，23(8):1403-1410.

陳國興，謝君斐，張克緒，等.1995.土的動模量和阻尼比的經(jīng)驗(yàn)估計(jì)［J］.地震工程與工程振動，15(1):73-84.

陳鈞，劉愛華.2010.基于Zlib的圖像差異壓縮算法研究與實(shí)現(xiàn)［J］.東華理工大學(xué)學(xué)報(bào):自然科學(xué)版，33(1):92-96.

阮周生，孫海.2010.同倫攝動法在一類線性積微分方程初值問題中的應(yīng)用［J］.東華理工大學(xué)學(xué)報(bào):自然科學(xué)版，33(3):298-300.

王福剛，曹劍鋒.2002.改進(jìn)的遺傳算法在地下水?dāng)?shù)值模擬中的應(yīng)用［J］.吉林大學(xué)學(xué)報(bào):地球科學(xué)版，21(1):64-68.

王錦國，周志芳.2002.基于實(shí)碼的遺傳算法的地下水污染物運(yùn)移參數(shù)反演［J］.水文，22(5):9-11.

楊曉華，金菊良，張國桃，等.1998.加速遺傳算法及其在暴雨強(qiáng)度公式參數(shù)優(yōu)化中的應(yīng)用［J］.自然災(zāi)害學(xué)報(bào)，7(3):71-76.

袁曉銘，孫銳，孫靜，等.2000.常規(guī)土類動剪切模量和阻尼比的試驗(yàn)研究［J］.地震工程與工程振動，20(4):133-139.

周明，孫樹棟.1999.遺傳算法原理及應(yīng)用［M］.北京:國防工業(yè)出版社.

周雙喜，楊彬.1996.影響遺傳算法性能的因素及改進(jìn)措施［J］.電力系統(tǒng)自動化，20(7):24-31.

Applied Research on Accelerating Genetic Algorithm Inversion Layer in the Shear Modulus，Damping Ratio Values of Model Parameters

LI Da-hu1， HE Qiang1， GU Qin-ping2
(1.Institute of Engineering Seismology of Earthquake Administration of Sichuan Province，Chengdu，SC 610041，China;2.Institute of Engineering Seismology of Earthquake Administration of Jiangsu Province，Nanjing，JS 210014，China)

In light of unique sedimentary environment of Chengdu region soil and the differences of the physical and mechanical properties，Martin-Davidenkov model and the empirical formula of damping ratio are used，the corresponding C++program is done through fact-based code and the introduction of the updated probability of accelerating genetic algorithm，by inversion，we can calculate dynamic shear modulus ratio G/Gmaxand damping ratio λ and the shear strain amplitude γ of the curve model parameters of the various types of deposits in Chengdu.The results show that:because of its fast convergence and high accuracy，RAGA is a nonlinear optimization method that not only can search the global optimal solution in greater probability，but also can search local details.

Chengdu region;accelerating genetic algorithm;inversion;dynamic shear modulus ratio;damping ratio

TU441+.3

A

1674-3504(2011)03-282-06

10.3969/j.issn.1674-3504.2011.03.013

2011-03-07

中國地震局地震科技星火計(jì)劃項(xiàng)目(XH1021Y)

李大虎(1982—)，男，理學(xué)碩士，從事工程場地地震安全性評價(jià)工作。E-mail:lixiang2006@sina.com