彭汝發(fā)，許小健

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

彭汝發(fā)，許小健

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

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

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

1 概 述

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

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

2 微進化算法基本原理

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

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

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

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

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

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

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

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

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

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

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

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

對于標準線性體三元件模型，蠕變試驗時，其本構(gòu)方程為

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

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

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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

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

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

［1］ 王宏貴，巍麗敏，赫曉光.根據(jù)長期單向壓縮試驗結(jié)果確定三維流變模型參數(shù)［J］.巖土工程學(xué)報，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ā)，陳新萬.多項式回歸間接求解巖石流變力學(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］ 李云鵬，王芝銀，丁秀麗.巖體原位流變荷載試驗的力學(xué)參數(shù)與模型反演［J］.實驗力學(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ù)識別［J］.巖石力學(xué)與工程學(xué)報，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é)報，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］.大壩觀測與土工測試，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］ 楊成祥，馮夏庭，陳炳瑞.基于擴展卡爾曼濾波的巖石流變模型參數(shù)識別［J］.巖石力學(xué)與工程學(xué)報，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ù)辨識［D］.濟南：山東科技大學(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é)報（自然科學(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-），男，安徽桐城人，高級工程師，主要從事巖土工程勘察設(shè)計方面的工作，（電話）0553-3827968（電子信箱）907137382＠qq.com

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