多孔隙介質(zhì)參數(shù)識(shí)別的正則化模擬退火法①

王迎

(黑龍江工業(yè)學(xué)院機(jī)電工程學(xué)院，黑龍江雞西158100)

隨著科學(xué)技術(shù)的飛速發(fā)展和人們生活水平的不斷提高，世界范圍內(nèi)面臨著油氣資源需耗量大幅增加的趨勢(shì)，而現(xiàn)有的油氣資源是有限的，如果不去進(jìn)一步勘察探尋油氣資源的分布，那么勢(shì)必會(huì)經(jīng)受資源急劇短缺的嚴(yán)峻問(wèn)題.可見(jiàn)，資源勘探技術(shù)尤為重要，其中涉及到在介質(zhì)的波動(dòng)方程中識(shí)別參數(shù)的問(wèn)題.研究人員在勘探資源之前需要預(yù)先知道一些參數(shù)，諸如孔隙率、飽和度等來(lái)還原對(duì)應(yīng)的介質(zhì)波動(dòng)方程，以此來(lái)預(yù)測(cè)所研究區(qū)域的地質(zhì)層面結(jié)構(gòu)，而要預(yù)知這些參數(shù)，就要通過(guò)參數(shù)識(shí)別方法來(lái)完成[1].另外，所勘探的地層通常是多孔介質(zhì)，更接近真實(shí)的地層介質(zhì)，因此其研究更具有實(shí)際意義[2].本文中有效結(jié)合了具有全局收斂特性的模擬退火法以及求解非適定問(wèn)題的正則化方法，形成了二維多孔隙介質(zhì)波動(dòng)方程的參數(shù)孔隙率識(shí)別問(wèn)題的正則化-模擬退火算法.

1 參數(shù)識(shí)別模型

描述二維多孔隙介質(zhì)波動(dòng)方程的模型為

固體和流體物質(zhì)的物理參數(shù)及各個(gè)參數(shù)間的關(guān)系如下所示，這些參數(shù)均為關(guān)于x和z的變量.

Kb：骨架體變模量；ρs：固體密度；Gb：骨架切變模量；m：有效流體密度；Ks：固體體變模量；ρf：液體密度；λb：固體骨架的Láme系數(shù)；λ?：雙相介質(zhì)的Láme系數(shù)；μ：孔隙介質(zhì)的切變模量；Kf：液體體變模量；β：孔隙率；α?、M：飽和多孔隙介質(zhì)壓縮性系數(shù)，且系數(shù)滿足

邊界條件：

初始條件：

(1)～(10)式構(gòu)成了二維雙相介質(zhì)波動(dòng)方程的正問(wèn)題模型，即由已知孔隙率β(x,z)可以求解出u(x,z,t).若孔隙率β(x,z)未知，給出附加條件

若對(duì)二維雙相介質(zhì)波動(dòng)方程(1)～(4)式考慮平行波入射時(shí)，f=g(t)δ(x)，其中g(shù)(t)為震源函數(shù)，δ(x)為δ函數(shù)，且滿足g(t)=0,t＜0，則(1)～(11)式就構(gòu)成了二維多孔隙介質(zhì)波動(dòng)方程一般意義下的參數(shù)識(shí)別問(wèn)題.

2 模型的離散化與參數(shù)識(shí)別算法

2.1 離散化模型設(shè)即分別為沿x,z方向的步長(zhǎng)，τ為時(shí)間步長(zhǎng)，記對(duì)方程(1)～(11)進(jìn)行差分離散.為了簡(jiǎn)單方便，將離散結(jié)果統(tǒng)一記作(12)式.

其中

其中

2.2 參數(shù)識(shí)別的算法記

對(duì)β(i,j)及分別按適當(dāng)順序排列，可排列成向量形式，分別記為B,G.

令

由(12)式定義了一個(gè)非線性向量值函數(shù)A:B→G.假設(shè)已知也進(jìn)行對(duì)應(yīng)排列，組成向

于是反演βi,j的問(wèn)題轉(zhuǎn)化成求解泛函

的非線性優(yōu)化問(wèn)題.

為了克服反問(wèn)題的不適定性，進(jìn)而達(dá)到數(shù)值求解的穩(wěn)定性，利用正則化方法在(13)式中引入光滑泛函，獲得泛函

其中，α為正則化參數(shù)，M1、M2分別為x方向、z方向的二階光滑矩陣[4].

通過(guò)求解滿足下式的最優(yōu)解B*，將要求解的二維流體飽和多孔隙介質(zhì)波動(dòng)方程反問(wèn)題再次轉(zhuǎn)化為求解非線性優(yōu)化問(wèn)題

Jα(B)的梯度為

其中，A′(B)是A(B)關(guān)于B的Jacobi矩陣.

模擬退火法作為隨機(jī)化算法[4]的通用方法之一，主要思想是通過(guò)設(shè)置隨機(jī)迭代函數(shù)來(lái)產(chǎn)生隨機(jī)解，最終在有限時(shí)間內(nèi)搜尋到符合終止條件的最優(yōu)解.根據(jù)所研究的問(wèn)題給出一個(gè)初始猜測(cè)值B0，通過(guò)(16)式迭代產(chǎn)生B1，B2，….我們期望存在一個(gè)Bk是(15)式的解或者存在某一收斂于解的點(diǎn)列.

其中c為產(chǎn)生擾動(dòng)的系數(shù)，rand(q+1,p+1)為隨機(jī)數(shù)矩陣.

在上述模型的基礎(chǔ)上，結(jié)合正則化方法，得到多孔隙介質(zhì)波動(dòng)方程反演的正則化模擬退火算法，其算法流程為：

步驟1:給定初始溫度T0，初始點(diǎn)B0為對(duì)應(yīng)階初值矩陣，及ε＞0，正則化參數(shù)α＞0,k?0.

否則，如果k+1=N，則Tn+1=σTn，B0?Bk+1,k?0，轉(zhuǎn)步驟2；

否則，轉(zhuǎn)步驟4.

步驟 4：計(jì)算 Δm=Jα(Bk+1)-Jα(Bk)，

若Δm＜0，則Bk=Bk+1，k?k+1，轉(zhuǎn)步驟2；

否則，轉(zhuǎn)步驟2.

上述算法中的kB=1.380 7×10-23J·K-1為Boltzmann常數(shù)，rand( )q+1,p+1為隨機(jī)數(shù)矩陣，c為產(chǎn)生擾動(dòng)的系數(shù)，Δm為目標(biāo)函數(shù)之差（評(píng)價(jià)函數(shù)的增量），σ為指數(shù)降溫的系數(shù)，N為B的維數(shù).我們可以用圖1簡(jiǎn)要描述此算法流程：

圖1 正則化模擬退火法流程圖

3 數(shù)值模擬

數(shù)值模擬使用形式為g(t)=sin(2πvt)exp(-2πvt)的雷克子波，其中入射波的頻率v取為75Hz.研究區(qū)域取為水平方向距離為260m、豎直方向?yàn)?60m的區(qū)域，即a=260，b=260.在進(jìn)行有限差分離散時(shí)，取水平差分步長(zhǎng)為hx=10m，豎直差分步長(zhǎng)為hz=10m，接收器在地表的每個(gè)差分點(diǎn)上.正則化參數(shù)α=0.000 95，相對(duì)誤差ε=1%，T0=1 050，c=1.5，σ=0.85，τ=0.002 s，參數(shù)識(shí)別時(shí)取孔隙率的初始值B0為所有元素都是0.2的相應(yīng)階矩陣.

所選模型如圖1所示，它是一種異常體，其中兩塊具有尖點(diǎn)，其相應(yīng)的各介質(zhì)參數(shù)如表1所示.

表1 多孔隙介質(zhì)模型的參數(shù)表

假設(shè)背景介質(zhì)的孔隙率的精確值為0.4，位于上方的異常體內(nèi)的孔隙率的精確值為0.25，位于下方的異常體內(nèi)的孔隙率的精確值為0.3.對(duì)此模型使用正則化-模擬退火算法進(jìn)行反演得到的結(jié)果如圖2所示.

4 結(jié)論

基于全局收斂的模擬退火方法與求解不適定性問(wèn)題的正則化方法的思想，將兩者有機(jī)結(jié)合構(gòu)造了正則化-模擬退火算法.為增加參數(shù)識(shí)別的收斂性和穩(wěn)定性，利用此算法識(shí)別二維多孔隙介質(zhì)波動(dòng)方程的參數(shù)孔隙率.數(shù)值模擬的結(jié)果顯示出此算法能克服一般的參數(shù)識(shí)別方法的局部收斂性，能有效地解決具有非適定性、非線性性質(zhì)的參數(shù)識(shí)別問(wèn)題.

圖1 孔隙率模型

圖2 正則化模擬退火算法對(duì)孔隙率模型的參數(shù)識(shí)別結(jié)果

[1]王立榮.反問(wèn)題的數(shù)值解法[D].深圳:哈爾濱工業(yè)大學(xué)深圳研究生院，2011.

[2]田迎春，章梓茂，趙成剛，等.黏彈性流體飽和孔隙介質(zhì)動(dòng)力反應(yīng)分析的顯式有限元法[J].科學(xué)技術(shù)與工程，2010，10(5):1112-1112.

[3]張新明，劉克安，劉家琦.流體飽和多孔隙介質(zhì)波動(dòng)方程多尺度反演[J].應(yīng)用力學(xué)學(xué)報(bào)，2007，24(1):87-90.

[4]黃鋒振，袁艷，張修寶，等.基于模擬退火的傅里葉變換成像光譜儀干涉圖相位修正[J].光譜學(xué)與光譜分析，2011，31(7):2012-2012.