王迎

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

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

1 參數(shù)識別模型

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

固體和流體物質(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ì)波動方程的正問題模型，即由已知孔隙率β(x,z)可以求解出u(x,z,t).若孔隙率β(x,z)未知，給出附加條件

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

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

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

其中

其中

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

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

令

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

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

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

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

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

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

Jα(B)的梯度為

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

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

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

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

步驟1:給定初始溫度T0，初始點(diǎn)B0為對應(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)生擾動的系數(shù)，Δm為目標(biāo)函數(shù)之差（評價(jià)函數(shù)的增量），σ為指數(shù)降溫的系數(shù)，N為B的維數(shù).我們可以用圖1簡要描述此算法流程：

圖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í)，取水平差分步長為hx=10m，豎直差分步長為hz=10m，接收器在地表的每個(gè)差分點(diǎn)上.正則化參數(shù)α=0.000 95，相對誤差ε=1%，T0=1 050，c=1.5，σ=0.85，τ=0.002 s，參數(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.對此模型使用正則化-模擬退火算法進(jìn)行反演得到的結(jié)果如圖2所示.

4 結(jié)論

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

圖1 孔隙率模型

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

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

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

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

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