彭汝發(fā)，許小健

用微進(jìn)化算法反演巖石蠕變模型非定常參數(shù)

彭汝發(fā)，許小健

（蕪湖市勘察測(cè)繪設(shè)計(jì)研究院，安徽蕪湖 241000）

巖石蠕變模型的參數(shù)較多，為得到參數(shù)的全局最優(yōu)解，應(yīng)用微進(jìn)化算法（Microevolution Algorithm，MA）對(duì)巖石蠕變模型非定常參數(shù)進(jìn)行了反演分析。算法以實(shí)測(cè)蠕變值與理論計(jì)算值之間的最小二乘誤差為優(yōu)化準(zhǔn)則函數(shù)，直接反演計(jì)算蠕變模型參數(shù)。計(jì)算結(jié)果表明，微進(jìn)化算法可最大限度地利用所有試驗(yàn)數(shù)據(jù)，避免傳統(tǒng)優(yōu)化算法初始參數(shù)選取的困難，且算法簡(jiǎn)單有效，計(jì)算精度高于混沌粒子群優(yōu)化算法。該方法也可推廣應(yīng)用于其它蠕變模型的參數(shù)反演，具有較高的工程應(yīng)用價(jià)值。

巖石蠕變模型；非定常參數(shù)反演；微進(jìn)化算法

1 概 述

巖石的蠕變研究對(duì)于合理評(píng)價(jià)巖體的長(zhǎng)期穩(wěn)定性，即巖石力學(xué)行為的時(shí)間效應(yīng)具有重要意義。巖石蠕變研究一般是從試驗(yàn)出發(fā)獲取有關(guān)數(shù)據(jù)，然后構(gòu)建巖石蠕變的本構(gòu)模型并確定相應(yīng)的模型參數(shù)。因此，根據(jù)蠕變?cè)囼?yàn)資料，如何正確確定巖石蠕變模型參數(shù)是巖石蠕變研究領(lǐng)域的關(guān)鍵課題之一。目前，確定巖石蠕變模型參數(shù)的方法主要有最小二乘法［1］、多項(xiàng)式回歸法［2］、流變曲線分解法［3］、復(fù)合形優(yōu)化法［4］、模式搜索算法［5］等。但巖石蠕變模型的參數(shù)眾多，為得到問(wèn)題的全局最優(yōu)解，須嘗試應(yīng)用新的優(yōu)化算法。近些年來(lái)，隨著優(yōu)化計(jì)算技術(shù)的發(fā)展，人工神經(jīng)網(wǎng)絡(luò)［6］、遺傳算法（Genetic Algorithm，GA）［7］、擴(kuò)展卡爾曼濾波技術(shù)［8］、粒子群優(yōu)化（Particle Swarm Optimization，PSO）算法［9］等逐漸在巖石蠕變模型參數(shù)確定中得到了應(yīng)用，為解決問(wèn)題提供了一種新的思路。

基于算法簡(jiǎn)單性和尋優(yōu)精度方面的考慮，本文給出了巖石蠕變模型非定常參數(shù)智能反演的一種新算法——微進(jìn)化算法 （Microevolution Algorithm，MA），并利用頁(yè)巖的蠕變?cè)囼?yàn)數(shù)據(jù)進(jìn)行了參數(shù)反演分析。計(jì)算結(jié)果表明，微進(jìn)化算法作為蠕變模型參數(shù)反演的新方法是可行有效的，且算法計(jì)算精度高于混沌粒子群優(yōu)化算法，具有較高的工程應(yīng)用價(jià)值。

2 微進(jìn)化算法基本原理

微進(jìn)化算法，是受人類(lèi)社會(huì)發(fā)展進(jìn)程中人類(lèi)趨同和趨異學(xué)習(xí)行為過(guò)程的啟發(fā)構(gòu)造出來(lái)的。它所采用的基本搜索策略是：群體中的每一個(gè)個(gè)體向群體中的優(yōu)秀個(gè)體學(xué)習(xí)，利用優(yōu)秀個(gè)體積累的經(jīng)驗(yàn)知識(shí)來(lái)改變自身，使自身得以處于不斷運(yùn)動(dòng)之中。具體來(lái)說(shuō)，MA基于實(shí)數(shù)編碼，采用群體社會(huì)學(xué)習(xí)機(jī)制，對(duì)種群中的每個(gè)個(gè)體i，以其自身所處歷史最優(yōu)位置為基礎(chǔ)，以群體最佳位置與當(dāng)前個(gè)體i的矢量差異為信息參照，進(jìn)行動(dòng)態(tài)隨機(jī)搜索，以實(shí)現(xiàn)種群的進(jìn)化。

為便于描述，設(shè) Xi＝（xi1，xi2，…，xiD）為個(gè)體 i＝1，2，…，NP的 D維矢量，每個(gè) Xi代表一個(gè)潛在的解；設(shè)Xpbesti為個(gè)體i迄今為止搜索到的最優(yōu)狀態(tài)、Xgbest為整個(gè)群體中的所有個(gè)體迄今為止搜索到的最優(yōu)位置狀態(tài)。這樣，在每一次的迭代搜索中，以式（1）更新各個(gè)體所處狀態(tài)：

式中：iter為迭代次數(shù)；r為0～1之間的正實(shí)數(shù)，一般為了使用方便，可令 r∈U（0，1）；N（0，σ）為正態(tài)分布隨機(jī)數(shù)，其中 σ∈（0.5，2.5）。

分析式（1）可知，當(dāng) Ni（0，σ）＞0時(shí)，群體中的各個(gè)體 i以矢量差（Xgbest－Xi，iter）為方向趨同；當(dāng)Ni（0，σ）＜0時(shí)，各個(gè)體以矢量差（Xgbest－Xi，iter）的反方向趨異。可見(jiàn)，式（1）將趨同與趨異有機(jī)地結(jié)合了起來(lái)，以實(shí)現(xiàn)局部的開(kāi)采與全局的探索進(jìn)化功能。算法實(shí)現(xiàn)步驟如下。

步驟1（初始化）：令 iter＝0，采用實(shí)數(shù)編碼，在可行域空間隨機(jī)初始化種群 Xiter＝（X1，iter，X2，iter，…，Xi，iter，…，XNP，iter）T，若搜索空間為 D維，則每個(gè)個(gè)體Xi，iter中包含了 D個(gè)變量，即 Xi，iter＝（xi1，iter，xi2，iter，…，xij，iter，…，xiD，iter），j＝1，2，…，D。

步驟2（種群評(píng)價(jià)）：評(píng)價(jià)種群中每一個(gè)個(gè)體，求出每一個(gè)個(gè)體到目前為止找到的最優(yōu)解，記為Xpbesti，i＝1，2，…，NP；將所有 Xpbesti中目標(biāo)函數(shù)值最 優(yōu) 的 個(gè) 體 記 為 Xgbest，即 Xgbest＝min｛Xpbest1，Xpbest2，…，Xgbesti，…XpbestNP｝。

步驟3（種群演化）：執(zhí)行如下操作，產(chǎn)生第iter＋1代種群 Xiter＋1。

步驟4（終止檢驗(yàn)）：①判斷算法循環(huán)執(zhí)行次數(shù)是否達(dá)到最大進(jìn)化代數(shù)ITERMAX；②判斷最優(yōu)目標(biāo)函數(shù)值f（Xgbest）是否達(dá)到預(yù)設(shè)精度VTR。若滿足終止條件其中之一，則結(jié)束算法，輸出Xgbest及其目標(biāo)函數(shù)值，否則，轉(zhuǎn)入步驟2直至滿足終止條件。

由算法的實(shí)現(xiàn)步驟可見(jiàn)，微進(jìn)化算法采用的基本策略是：①采用實(shí)數(shù)編碼；②使用固定規(guī)模的種群NP；③采用貪婪的種群更新策略，即用NP個(gè)新個(gè)體Xi，iter＋1完全替換 NP個(gè)舊個(gè)體 Xi，iter；④采用精英保留策略，即更新個(gè)體自身所處歷史最優(yōu)狀態(tài)Xpbesti和群體最佳狀態(tài)Xgbest。

MA需要設(shè)置的算法控制參數(shù)較少，除算法的運(yùn)行控制參數(shù)（ITERMAX、VTR）外，僅需設(shè)置群體數(shù)目NP和σ。過(guò)大的NP會(huì)影響算法的運(yùn)算速度，經(jīng)試驗(yàn)，對(duì)于一般問(wèn)題，建議 NP取20～40左右；而σ∈（0.5，2.5），一般取 1即可獲得較滿意解?？梢?jiàn)，與GA和PSO算法相比，MA控制參數(shù)的選擇更為簡(jiǎn)單，這使得算法易于編程實(shí)現(xiàn)和便于用戶使用。

3 微進(jìn)化算法在巖石蠕變模型非定常參數(shù)反演中的應(yīng)用

為探討本文的微進(jìn)化算法優(yōu)化反演巖石蠕變模型參數(shù)的可行性，利用文獻(xiàn)［10］中的頁(yè)巖蠕變?cè)囼?yàn)數(shù)據(jù)，采用標(biāo)準(zhǔn)線性體三元件模型，并引入非定常參數(shù)后進(jìn)行反演分析。

對(duì)于標(biāo)準(zhǔn)線性體三元件模型，蠕變?cè)囼?yàn)時(shí)，其本構(gòu)方程為

在蠕變本構(gòu)方程中引入非定常參數(shù)EK（t）＝p1＋p2ep3t，并令 ηK＝p4，EH＝p5，式（2）變?yōu)椋?］

設(shè)式（3）等式左邊和右邊項(xiàng)分別記為fL和fR，則根據(jù)式（3），由 m對(duì)實(shí)測(cè)蠕變數(shù)據(jù)｛tl，εl｝（l＝1，2，…，m）可建立反演問(wèn)題的優(yōu)化準(zhǔn)則函數(shù)為

式中，X＝（p1，p2，p3，p4，p5）為待反演參數(shù)。由于式（4）具有高度的非線性特征，傳統(tǒng)優(yōu)化方法如高斯-牛頓或麥夸特法等，對(duì)輸入的模型參數(shù)初始值有一定的要求，若輸入的參數(shù)值與參數(shù)的真實(shí)值偏離較大，可能造成迭代不收斂，而且也很難得到全局優(yōu)化解。而本文微進(jìn)化算法對(duì)模型是否線性、連續(xù)等沒(méi)有限制，也不受優(yōu)化變量數(shù)目的束縛，直接在優(yōu)化準(zhǔn)則函數(shù)的引導(dǎo)下進(jìn)行全局自適應(yīng)尋優(yōu)，該方法簡(jiǎn)便、通用、適應(yīng)性強(qiáng)。因此，現(xiàn)利用MA進(jìn)行反演計(jì)算，步驟如下。

步驟1：初始化NP組待反演的蠕變參數(shù)Xi∈（p1，p2，p3，p4，p5），i＝1，…，NP，作為初始種群。

步驟2：將m組蠕變?cè)囼?yàn)數(shù)據(jù)和初始化得NP組待反演的蠕變參數(shù)代入已建立的反演問(wèn)題的優(yōu)化準(zhǔn)則函數(shù)式（4）中，計(jì)算得到相應(yīng)應(yīng)力水平下的最小二乘誤差。

步驟3：評(píng)價(jià)各組待反演蠕變參數(shù)Xi，更新Xgbest。

步驟4：利用式（1）策略進(jìn)行微進(jìn)化尋優(yōu)操作，并重復(fù)執(zhí)行步驟3，直至滿足微進(jìn)化算法的終止條件，然后輸出待反演蠕變參數(shù)的結(jié)果以達(dá)到最優(yōu)化的目的。

根據(jù)上述蠕變參數(shù)的優(yōu)化反演步驟，可反演計(jì)算得到頁(yè)巖的蠕變模型非定常參數(shù)，如表1。為便于與文獻(xiàn)［9］提出的混沌粒子群優(yōu)化（Chaos Particle Swarm Optimization，CPSO）算法進(jìn)行比較，表1中也列出了CPSO算法反演的參數(shù)結(jié)果。比較可見(jiàn)，在各應(yīng)力水平作用下，MA與CPSO算法反演的參數(shù)p1，p5結(jié)果差別不大，而參數(shù) p2，p3，p4差別較大，相差達(dá)一個(gè)數(shù)量級(jí)之多。由優(yōu)化準(zhǔn)則函數(shù)值min J結(jié)果可見(jiàn)，MA計(jì)算所得非定常參數(shù)的最小二乘誤差遠(yuǎn)小于CPSO計(jì)算所得的最小二乘誤差，MA可將優(yōu)化準(zhǔn)則函數(shù)值min J較少50%左右，這說(shuō)明MA較CPSO算法在處理式（4）所示這種較為復(fù)雜的非線性模型時(shí)，具有更高的尋優(yōu)精度。

為便于直觀反映MA與CPSO算法反演不同應(yīng)力水平各時(shí)刻的應(yīng)變值，圖1繪出了兩種算法所得蠕變計(jì)算值與試驗(yàn)值結(jié)果。從圖1可以看出，各應(yīng)力水平下，MA計(jì)算所得理論反演變形曲線與試驗(yàn)變形曲線較為吻合，誤差較小，說(shuō)明本文提出的MA是有效可靠的；而在應(yīng)力水平σc＝90 MPa，σc＝300 MPa時(shí)，CPSO算法計(jì)算所得的后期理論變形曲線偏離試驗(yàn)曲線相對(duì)較大。

表1 微進(jìn)化算法與混沌粒子群算法反演結(jié)果的比較Table 1 Comparison between MA and CPSO inversion results

圖1 試驗(yàn)數(shù)據(jù)及擬合曲線Fig.1 The test data and the fitting curves

4 結(jié) 論

（1）以實(shí)測(cè)蠕變值與理論計(jì)算值之間的最小二乘誤差為優(yōu)化準(zhǔn)則函數(shù)，用微進(jìn)化算法對(duì)巖石蠕變模型非定常參數(shù)進(jìn)行優(yōu)化反演，能最大限度地利用所有試驗(yàn)數(shù)據(jù)，且可避免傳統(tǒng)優(yōu)化算法初始參數(shù)選取的困難。算例計(jì)算結(jié)果表明，微進(jìn)化算法不受初始選值的影響，計(jì)算精度高于混沌粒子群優(yōu)化算法，具有較高的工程應(yīng)用價(jià)值。

（2）本文雖以標(biāo)準(zhǔn)線性體三元件模型為例，進(jìn)行了非定常參數(shù)的反演分析，但方法同樣適用于其他蠕變模型，只需采用與模型對(duì)應(yīng)的蠕變本構(gòu)方程作為優(yōu)化準(zhǔn)則即可。

（3）微進(jìn)化算法簡(jiǎn)單有效，且易于編程實(shí)現(xiàn)和便于用戶使用。

［1］ 王宏貴，巍麗敏，赫曉光.根據(jù)長(zhǎng)期單向壓縮試驗(yàn)結(jié)果確定三維流變模型參數(shù)［J］.巖土工程學(xué)報(bào)，2006，28（5）：669－673.（WANG Hong-gui，WEI Li-min，HE Xiao-guang. Determination of Three-Dimensional Rheologic Parameters Based on Long Term Oedometer Test［J］.Chinese Journal of Geotechnical Engineering，2006，28（5）：669－673.（in Chinese））

［2］ 許宏發(fā)，陳新萬(wàn).多項(xiàng)式回歸間接求解巖石流變力學(xué)參數(shù)的方法［J］.有色金屬，1994，46（4）：19－22.（XU Hong-fa，CHEN Xin-wan.Method to Find Indirect Solution of Rock Rheology’s Parameters with Polynomial Regression［J］.Nonferrous Metals，1994，46（4）：19－22.（in Chinese））

［3］ BOUKHAROV G N，CHANDA M W，BOUKHAROV N G.The Three Processes of Brittle Crystalline Rock Creep［J］.International Journal of Rock Mechanics，Mining Science and Geomechanics Abstract，1995，32（4）：325－335.

［4］ 李云鵬，王芝銀，丁秀麗.巖體原位流變荷載試驗(yàn)的力學(xué)參數(shù)與模型反演［J］.實(shí)驗(yàn)力學(xué)，2005，20（2）：297－303.（LI Yun-peng，WANG Zhi-yin，DING Xiuli.Back Analysis of Mechanical Parameters and Model of Rheological In-situ Loading Test for Rock Mass［J］.Journal of Experimental Mechanics，2005，20（2）：297－303.（in Chinese））

［5］ 陳炳瑞，馮夏庭，丁秀麗，等.基于模式搜索的巖石流變模型參數(shù)識(shí)別［J］.巖石力學(xué)與工程學(xué)報(bào)，2005，24（2）：207－211.（CHEN Bing-rui，F(xiàn)ENG Xia-ting，DING Xiu-li，et al.Rheological Model and Parameters I-dentification of Rock Based on Pattern Search and Least-Square Techniques［J］.Chinese Journal of Rock Mechanics and Engineering，2005，24（2）：207－211.（in Chinese））

［6］ 陳炳瑞，馮夏庭，丁秀麗，等.基于模式-遺傳-神經(jīng)網(wǎng)絡(luò)的流變參數(shù)反演［J］.巖石力學(xué)與工程學(xué)報(bào)，2005，24（4）：553－558.（CHEN Bing-rui，F(xiàn)ENG Xia-ting，DINGXiu-li，et al.Back Analysis on Rheological Parameters Based on Pattern-Genetic-Neural Network［J］.Chinese Journal of Rock Mechanics and Engineering，2005，24（4）：553－558.（in Chinese））

［7］ 劉 杰，王 媛.巖體流變參數(shù)反演的加速遺傳算法［J］.大壩觀測(cè)與土工測(cè)試，2001，25（6）：24－27.（LIU Jie，WANG Yuan.Back Analysis of Rock Mass Rheologic Parameters Based on Genetic Algorithm［J］.Dam Observation and Geotechnical Tests，2001，25（6）：24－27.（in Chinese））

［8］ 楊成祥，馮夏庭，陳炳瑞.基于擴(kuò)展卡爾曼濾波的巖石流變模型參數(shù)識(shí)別［J］.巖石力學(xué)與工程學(xué)報(bào)，2007，26（4）：754－761.（YANG Cheng-xiang，F(xiàn)ENG Xia-ting，CHEN Bing-rui. Parameter Identification of Rock Rheological Model Based on Extended Kalman Filter［J］.Chinese Journal of Rock Mechanics and Engineering，2007，26（4）：754－761.（in Chinese））

［9］ 丁志坤.巖石粘彈性性非定常蠕變方程的參數(shù)辨識(shí)［D］.濟(jì)南：山東科技大學(xué)，2003.（DING Zhi-kun.Parameter Identification of Non-stationary Creep Equation of Rock.Jinan：Shandong University of Science and Technology，2003.（in Chinese））

［10］李志敬，朱珍德，周偉華.基于CPSO算法的巖石蠕變模型非定常參數(shù)反演分析［J］.河海大學(xué)學(xué)報(bào)（自然科學(xué)版），2008，36（3）：346－349.（LI Zhi-jing，ZHU Zhen-de，ZHOU Wei-hua.Back Analysis of Nonstationary Parameters of Rock Creep Model Based on Chaos Particle Swarm Optimization Algorithm［J］.Journal of Hohai University（Natural Sciences），2008，36（3）：346－349.（in Chinese） ）

Microevolution Algorithm for Inversion of Non-stationary Parameters in Rock Creep Model

PENG Ru-fa，XU Xiao-jian
（Wuhu Geotechnical and Survey Design Institute，Wuhu 241000，China）

Rock creep model generally contains several parameters.To obtain the global optimal solution of the parameters，Microevolution Algorithm（MA）was employed for the inversion of non-stationary parameters.In this paper，parameters of the creep model are directly inversed with the least square error between measured creep values and calculated creep values as the optimization criterion function.The computation results demonstrate that microevolution algorithm can maximize the use of all test data and avoid the difficulty of selecting initial parameters in traditional optimization algorithm.Moreover，microevolution algorithm is simple and effective，and offers higher accuracy than Chaos Particle Swarm Optimization（CPSO）.In this sense，it can be applied to the inversion of parameters in other creep models and has high application value for engineering.

rock creep model；inversion of non-stationary parameter；microevolution algorithm

TU425

A

1001－5485（2011）06－0050－05

2010-07-21

彭汝發(fā)（1964-），男，安徽桐城人，高級(jí)工程師，主要從事巖土工程勘察設(shè)計(jì)方面的工作，（電話）0553-3827968（電子信箱）907137382＠qq.com

（編輯：趙衛(wèi)兵）