基于優(yōu)化時(shí)空域頻散關(guān)系的聲波方程有限差分最小二乘逆時(shí)偏移

薛 浩 劉 洋 楊宗青

(①中國(guó)石油大學(xué)(北京)油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 102249； ②中國(guó)石油大學(xué)(北京)克拉瑪依校區(qū)石油學(xué)院，新疆克拉瑪依 834000； ③中國(guó)石油大學(xué)(北京)CNPC物探重點(diǎn)實(shí)驗(yàn)室，北京 102249；④中國(guó)石油集團(tuán)東方地球物理公司西南物探分公司，四川成都 610213)

1 引言

隨著油氣勘探程度的不斷深入，石油勘探的重點(diǎn)已經(jīng)轉(zhuǎn)向復(fù)雜地質(zhì)構(gòu)造目標(biāo)和巖性勘探，并且對(duì)勘探精度和計(jì)算效率的要求也逐漸提高。偏移方法是對(duì)復(fù)雜構(gòu)造成像的一種有效手段，從基于射線的Kirchhoff偏移[1]和Beam偏移[2，3]，發(fā)展到基于波動(dòng)方程的單程波偏移[4-6]和逆時(shí)偏移(RTM)[7-9]。由于常規(guī)偏移算子是正演算子的共軛轉(zhuǎn)置而不是其逆[10]，受復(fù)雜構(gòu)造及采集條件的影響，常規(guī)偏移方法難以滿足巖性油氣藏勘探開(kāi)發(fā)的需要。為了解決這一問(wèn)題，國(guó)外學(xué)者基于Tarantola[11]提出的最小二乘反演理論，提出了最小二乘偏移(LSM)，通過(guò)模擬數(shù)據(jù)與實(shí)際數(shù)據(jù)的匹配求取最佳成像剖面[10]。經(jīng)過(guò)不斷的發(fā)展，最小二乘偏移從早期的Kirchhoff偏移[12]和單程波波動(dòng)方程偏移[13]逐漸推廣應(yīng)用到逆時(shí)偏移[14-21]中。相比常規(guī)逆時(shí)偏移方法，基于最小二乘的逆時(shí)偏移方法能有效壓制偏移噪聲，生成分辨率更高、保幅性更好的成像剖面，在成像效果方面更具優(yōu)勢(shì)[14-21]。

最小二乘逆時(shí)偏移(LSRTM)的關(guān)鍵技術(shù)之一是波動(dòng)方程的求解，因此波動(dòng)方程的求解精度和效率必然對(duì)LSRTM的成像精度和計(jì)算效率產(chǎn)生重要影響。波動(dòng)方程求解方法主要分為三類(lèi)：有限差分法、偽譜法和有限元法[22]。有限差分法由于計(jì)算量較小且實(shí)現(xiàn)簡(jiǎn)單，在地震數(shù)值模擬、偏移和反演中已得到廣泛應(yīng)用[22]。在求解波動(dòng)方程時(shí)，通常利用高階有限差分法計(jì)算空間導(dǎo)數(shù)[23]。但地震波的傳播不僅與空間有關(guān)，還與時(shí)間有關(guān)，采用傳統(tǒng)的高階有限差分法求解差分系數(shù)時(shí)，一般是基于空間域頻散關(guān)系，通過(guò)對(duì)空間導(dǎo)數(shù)項(xiàng)進(jìn)行Taylor級(jí)數(shù)展開(kāi)得到，這樣會(huì)導(dǎo)致較大的頻散和誤差。為了提高有限差分的模擬精度，Liu等[24，25]基于時(shí)空域頻散關(guān)系，利用平面波理論和Taylor級(jí)數(shù)展開(kāi)法提出了一種時(shí)空域有限差分法，并將其用于求解聲波方程。與傳統(tǒng)有限差分法相比，時(shí)空域有限差分法的精度和穩(wěn)定性都得到了提高。為了進(jìn)一步提高求解精度，許多學(xué)者通過(guò)采用優(yōu)化等方式，在求取有限差分系數(shù)時(shí)，使差分算子最大程度地逼近空間偏導(dǎo)數(shù)算子[26-30]。Liu[31]用最小二乘法優(yōu)化時(shí)空域頻散關(guān)系求解聲波方程優(yōu)化差分系數(shù)。與Taylor級(jí)數(shù)展開(kāi)方法相比，利用優(yōu)化時(shí)空域頻散關(guān)系求解差分系數(shù)可進(jìn)一步減小頻散，提高模擬精度。

本文在上述研究基礎(chǔ)上，將優(yōu)化的時(shí)空域頻散關(guān)系推廣到LSRTM中，在波動(dòng)方程求解過(guò)程中利用優(yōu)化時(shí)空域頻散關(guān)系求解差分系數(shù)。在進(jìn)行波場(chǎng)外推時(shí)，利用混合吸收邊界條件處理邊界反射。通過(guò)Clearbout共軛梯度算法[32]對(duì)LSRTM的計(jì)算流程進(jìn)行改進(jìn)，減少波動(dòng)方程的求解次數(shù)。數(shù)值模型的計(jì)算結(jié)果表明了本文方法的有效性和優(yōu)越性。

2 基于優(yōu)化時(shí)空域的聲波波動(dòng)方程有限差分法

2.1 基于優(yōu)化時(shí)空域頻散關(guān)系的聲波波動(dòng)方程差分系數(shù)求解

波動(dòng)方程的求解是最小二乘逆時(shí)偏移的關(guān)鍵技術(shù)，波動(dòng)方程求解的精度和效率在一定程度上影響了偏移成像、反演的精度和計(jì)算效率。在二維各向同性介質(zhì)中，聲波方程可以寫(xiě)為

(1)

式中：v為速度；p為標(biāo)量波場(chǎng)；t為時(shí)間；x、z為空間坐標(biāo)。

有限差分算法由于其較高的計(jì)算效率和簡(jiǎn)單的實(shí)現(xiàn)方式而被廣泛應(yīng)用于波動(dòng)方程求解[22]。采用二階中心差分計(jì)算時(shí)間二階導(dǎo)數(shù)，可得

(2)

采用高階差分計(jì)算空間導(dǎo)數(shù)，可得x和z方向的空間導(dǎo)數(shù)為

(3)

(4)

式中cm為差分系數(shù)。

將式(2)～式(4)代入式(1)，可得

(5)

≈r-2[-1+cos(ωτ)]

(6)

式中

令kx=kcosθ和kz=ksinθ，其中θ為平面波的傳播角度。化簡(jiǎn)式(6)得到

≈r-2[-1+cos(ωτ)]

(7)

利用泰勒(Taylor)級(jí)數(shù)展開(kāi)法展開(kāi)余弦函數(shù)

(8)

比較k2j的系數(shù)，可得

(9)

j=1,2,…,M

(10)

當(dāng)r=0時(shí)，相當(dāng)于只逼近了空間域的頻散關(guān)系，即只利用了空間導(dǎo)數(shù)來(lái)求解差分系數(shù)。此時(shí)求得的是傳統(tǒng)有限差分法的差分系數(shù)，即由Taylor展開(kāi)空間域有限差分法求解的。

當(dāng)r≠0時(shí)，相當(dāng)于逼近了時(shí)間域和空間域的頻散關(guān)系，即利用了空間導(dǎo)數(shù)和時(shí)間導(dǎo)數(shù)同時(shí)求解差分系數(shù)，此時(shí)求得的差分系數(shù)為T(mén)aylor展開(kāi)時(shí)空域有限差分而得到的。

有限差分的差分系數(shù)cm在一定程度上決定了有限差分法的精度和效率，傳統(tǒng)的有限差分法是有計(jì)算簡(jiǎn)單、實(shí)現(xiàn)方便的特點(diǎn)，但在求解差分系數(shù)時(shí)只利用泰勒級(jí)數(shù)展開(kāi)逼近空間域的頻散關(guān)系，在中高頻區(qū)域數(shù)值頻散迅速增大。為了提高求解波動(dòng)方程的精度，Liu[31]發(fā)展了一種最小二乘優(yōu)化方法，通過(guò)極小化式(7)左端與右端的相對(duì)誤差這一目標(biāo)函數(shù)逼近時(shí)空域的頻散關(guān)系，從而求得基于優(yōu)化時(shí)空域頻散關(guān)系的差分系數(shù)。優(yōu)化時(shí)空域有限差分法的目標(biāo)函數(shù)為

(11)

令目標(biāo)函數(shù)等于0，可得

(12)

通過(guò)求解上述方程組可得優(yōu)化差分系數(shù)cm。

2.2 數(shù)值頻散分析

相速度的相對(duì)誤差ξ可用來(lái)衡量有限差分法的精度，公式如下[31]

(13)

式中vFD代表數(shù)值模擬速度，且有A=cos(mβsinθ)，B=cos(mβcosθ)。

采用一個(gè)數(shù)值算例對(duì)其頻散誤差進(jìn)行分析。圖1為固定M時(shí)，傳統(tǒng)有限差分法、Taylor展開(kāi)時(shí)空域有限差分法和優(yōu)化時(shí)空域有限差分法隨角度θ變化的頻散誤差曲線。圖2為角度θ固定時(shí)，傳統(tǒng)有限差分法、Taylor展開(kāi)時(shí)空域有限差分法和優(yōu)化時(shí)空域有限差分法隨M變化的頻散誤差曲線。從圖1和圖2可看出，隨著M的增大，三種方法都能在更大的波數(shù)范圍內(nèi)達(dá)到較小的誤差，Taylor級(jí)數(shù)展開(kāi)時(shí)空域的有限差分法相比傳統(tǒng)有限差分法數(shù)值頻散要小，優(yōu)化時(shí)空域有限差分法的數(shù)值頻散更小，精度更高。在誤差相差不大時(shí)，優(yōu)化時(shí)空域方法的算子長(zhǎng)度相比傳統(tǒng)的有限差分法的算子長(zhǎng)度更短，因此優(yōu)化時(shí)空域方法能夠利用較短的算子長(zhǎng)度達(dá)到傳統(tǒng)方法更高的精度，進(jìn)而提高波動(dòng)方程的求解效率。

圖1 M固定時(shí)傳統(tǒng)有限差分法(a)、Taylor級(jí)數(shù)展開(kāi)時(shí)空域有限差分法(b)和優(yōu)化時(shí)空域有限差分法(c)隨θ變化的頻散誤差曲線

圖2 θ固定時(shí)傳統(tǒng)有限差分法(a)、Taylor級(jí)數(shù)展開(kāi)時(shí)空域有限差分法(b)和優(yōu)化時(shí)空域有限差分法(c)隨M變化的頻散誤差曲線

3 混合吸收邊界條件

對(duì)不采用吸收邊界、采用PML吸收邊界[34]和采用混合吸收邊界的情況進(jìn)行數(shù)值模擬對(duì)比，說(shuō)明混合吸收邊界的吸收效果。為了方便對(duì)比，PML吸收邊界采用二階聲波方程PML邊界條件[34]，邊界部分采用指數(shù)衰減函數(shù)。 均勻介質(zhì)模型尺寸為

2000m×2000m，速度為3000m/s，空間采樣間隔為10m，時(shí)間采樣間隔為1ms，震源位于模型正中央，震源為一個(gè)主頻為20Hz的Ricker子波。圖4a、圖4b和圖4c分別為不采用吸收邊界、采用10層PML吸收邊界和采用10層混合吸收邊界情況下三個(gè)時(shí)刻的波場(chǎng)快照，圖4d為在(z=500m，x=1000m)位置處記錄的歸一化后的波形記錄。從圖4中可看出，在沒(méi)有吸收邊界的情況下，會(huì)有明顯的邊界反射，采用10層PML吸收邊界后，仍然存在一定的邊界反射，采用10層混合吸收邊界條件后，邊界反射幾乎完全被吸收。

圖3 混合吸收邊界條件示意圖[33]

圖4 吸收邊界對(duì)比圖

4 最小二乘法逆時(shí)偏移原理及計(jì)算流程優(yōu)化

散射波場(chǎng)的Born正演可寫(xiě)為[16]

(14)

式中：G(g|x)表示震源在x位置、檢波點(diǎn)在g位置的格林函數(shù)；m為反射系數(shù)模型；d為Born正演數(shù)據(jù)。

方程可寫(xiě)為矩陣形式

d=Lm

(15)

式中L是正演算子。LSRTM成像的目標(biāo)是用正演數(shù)據(jù)d估計(jì)反射系數(shù)

(16)

寫(xiě)成矩陣形式為

m=LTd

(17)

式中偏移算子LT是正演算子L的伴隨陣。

為了獲得更好的成像結(jié)果，成像問(wèn)題可以當(dāng)作最小二乘反演問(wèn)題來(lái)處理，目標(biāo)函數(shù)定義為

(18)

迭代最速下降解為

m(k+1)=m(k)-αLT[Lm(k)-d]

(19)

式中：α是步長(zhǎng)；k是迭代次數(shù)。該式是LSRTM在數(shù)據(jù)域的迭代。反射系數(shù)模型用來(lái)生成Born正演數(shù)據(jù)，通過(guò)對(duì)觀測(cè)數(shù)據(jù)與Born正演數(shù)據(jù)求差，得到數(shù)據(jù)的殘差； 用殘差結(jié)果作為輸入地震記錄進(jìn)行逆時(shí)偏移計(jì)算求得梯度，梯度結(jié)果用來(lái)修正反射系數(shù)模型。這個(gè)過(guò)程一直重復(fù)，直至找到一個(gè)確定的反射系數(shù)模型。

最小二乘逆時(shí)偏移通常采用共軛梯度法(CG)求取反射系數(shù)模型[18]，計(jì)算流程可用圖5表示[35]。梯度的計(jì)算相當(dāng)于一次逆時(shí)偏移，用邊界重建方法實(shí)現(xiàn)需要求解正傳、重建和反傳一共3次波動(dòng)方程，Born正演需要求解2次波動(dòng)方程，用線性搜索方法計(jì)算步長(zhǎng)也需要進(jìn)行2次Born正演，即求解4次波動(dòng)方程，每次迭代總共需要求解9次波動(dòng)方程。Clearbout對(duì)共軛梯度方法進(jìn)行改進(jìn)，通過(guò)求解線性方程組來(lái)直接求取搜索步長(zhǎng)[32]。

設(shè)si為第i次迭代時(shí)反射系數(shù)模型的更新量，gi為梯度，R為殘差，則有

Si=Lsi

(20)

Gi=Lgi

(21)

R=d-Lm

(22)

定義α和β為搜索步長(zhǎng)，則αG+βS為共軛梯度法搜索平面，求解目標(biāo)是找到合適步長(zhǎng)，使Q(α,β)=(R-αG-βS)·(R-αG-βS)的值最小。

(23)

求解該方程組可得到系數(shù)α和β，然后更新殘差和反射系數(shù)模型。利用優(yōu)化CG方法進(jìn)行LSRTM計(jì)算的流程如圖6所示。優(yōu)化CG方法避免了線性搜索中的多次正演，每次迭代減少了4次波動(dòng)方程的求解，提高了最小二乘逆時(shí)偏移的計(jì)算效率。

圖5 常規(guī)LSRTM計(jì)算流程[35]圖6 優(yōu)化CG方法的LSRTM計(jì)算流程

5 模型試算

正演模擬的精度和效率直接影響反演的精度和效率，利用Taylor級(jí)數(shù)展開(kāi)時(shí)空域有限差分法和優(yōu)化時(shí)空域有限差分法，分別對(duì)均勻模型和Marmousi模型進(jìn)行數(shù)值模擬，分析對(duì)比兩種有限差分法對(duì)正演模擬和LSRTM成像結(jié)果的影響。計(jì)算過(guò)程中，為了提高計(jì)算效率，采用了GPU進(jìn)行加速。

5.1 均勻模型

均勻介質(zhì)模型的尺寸為4000m×4000m，模型速度為1800m/s，空間采樣間隔為20m，時(shí)間采樣間隔為1ms，震源為50Hz的Ricker子波，震源位于模型正中央。圖7a～圖7c分別為利用Taylor級(jí)數(shù)展開(kāi)時(shí)空域有限差分法(M=8、M=30)和優(yōu)化時(shí)空域有限差分法(M=8)得到的2000ms時(shí)刻的波場(chǎng)快照； 圖7d為從圖7a～圖7c中x=2000m位置處抽取的單道波形對(duì)比。從圖7d的波形對(duì)比可以看出，當(dāng)算子長(zhǎng)度相同時(shí)(M=8)，Taylor級(jí)數(shù)展開(kāi)時(shí)空域有限差分法(圖7a)的數(shù)值頻散明顯高于優(yōu)化時(shí)空域有限差分法(圖7c)。算子長(zhǎng)度較短時(shí)的優(yōu)化時(shí)空域有限差分法(圖7c，M=8)與算子長(zhǎng)度較長(zhǎng)的Taylor級(jí)數(shù)展開(kāi)時(shí)空域限差分法(圖7b，M=30)數(shù)值頻散相當(dāng)。在數(shù)值頻散相差不大的情況下，優(yōu)化時(shí)空域有限差分法可以采用較短的算子長(zhǎng)度提高計(jì)算效率，在算子長(zhǎng)度相同的情況下，優(yōu)化時(shí)空域有限差分法有更高的精度。

圖7 均勻介質(zhì)波場(chǎng)快照對(duì)比圖

5.2 Marmousi模型

Marmousi速度模型如圖8所示，模型尺寸為5000m×3520m，空間采樣間隔為10m，震源為50Hz的Ricker子波，震源位于(x=2500m，z=0)。圖9a和圖9b分別為利用Taylor級(jí)數(shù)展開(kāi)時(shí)空域有限差分法(M=4)和優(yōu)化時(shí)空域有限差分法(M=4)得到的1800ms時(shí)刻的波場(chǎng)快照。圖10為得到的單炮地震記錄對(duì)比圖。從圖9和圖10中可看出，在復(fù)雜速度模型情況下，算子長(zhǎng)度相同時(shí)(M=4)，Taylor級(jí)數(shù)展開(kāi)時(shí)空域有限差分法(圖9a、圖10a)的數(shù)值頻散明顯高于優(yōu)化時(shí)空域有限差分法(圖9b、圖10b)，采用相同的算子長(zhǎng)度，對(duì)于復(fù)雜速度模型，優(yōu)化時(shí)空域有限差分法仍然能獲得更高的模擬精度。

用同一模型進(jìn)行逆時(shí)偏移和最小二乘逆時(shí)偏移成像，炮點(diǎn)和檢波點(diǎn)數(shù)目分別為50和353，炮點(diǎn)和檢波點(diǎn)均勻分布于地表，炮點(diǎn)間隔為100m，檢波點(diǎn)間隔為10m，最小二乘逆時(shí)偏移的迭代次數(shù)為20次。

圖11為反射系數(shù)模型，圖12為常規(guī)逆時(shí)偏移經(jīng)Laplacian濾波后的成像結(jié)果。圖13a和圖13b分別為利用Taylor級(jí)數(shù)展開(kāi)時(shí)空域有限差分法(M=4)和優(yōu)化時(shí)空域有限差分法(M=4)得到的最小二乘逆時(shí)偏移成像結(jié)果； 圖13c和圖13d分別為圖13a和圖13b的局部放大圖。圖14為x= 3300m處的反射系數(shù)模型(REF)、常規(guī)逆時(shí)偏移求得的反射系數(shù)(RTM)和優(yōu)化時(shí)空域最小二乘逆時(shí)偏移求得的反射系數(shù)(LSM)對(duì)比圖。從圖11～圖14的對(duì)比可看出，相比逆時(shí)偏移成像，最小二乘逆時(shí)偏移成像精度更高、保幅性更好，最小二乘逆時(shí)偏移的結(jié)果更接近于真實(shí)的反射系數(shù)模型。

圖8 Marmousi速度模型

圖9 波場(chǎng)快照對(duì)比圖

圖10 單炮記錄對(duì)比圖

圖11 反射系數(shù)模型 圖12 逆時(shí)偏移經(jīng)Laplacian濾波后成像結(jié)果

圖13 Marmousi模型的LSRTM成像結(jié)果

圖15為x= 3300m處的反射系數(shù)模型(REF)、Taylor級(jí)數(shù)展開(kāi)時(shí)空域最小二乘逆時(shí)偏移求得的反射系數(shù)(TEM)和優(yōu)化時(shí)空域最小二乘逆時(shí)偏移求得的反射系數(shù)(LSM)對(duì)比圖。從圖11、圖13和圖15的對(duì)比可看出，在差分算子長(zhǎng)度相同時(shí)，Taylor級(jí)數(shù)展開(kāi)時(shí)空域有限差分最小二乘逆時(shí)偏移和優(yōu)化時(shí)空域有限差分最小二乘逆時(shí)偏移在強(qiáng)反射系數(shù)界面處都有著較好的成像結(jié)果，但在弱反射界面和圖13c、圖13d紅框部分中的高速夾層處，優(yōu)化時(shí)空域有限差分最小二乘逆時(shí)偏移具有更高精度的成像結(jié)果。圖16為T(mén)aylor級(jí)數(shù)展開(kāi)時(shí)空域有限差分最小二乘逆時(shí)偏移與優(yōu)化時(shí)空域有限差分最小二乘逆時(shí)偏移的收斂曲線對(duì)比，可看出優(yōu)化時(shí)空域方法具有更快的收斂速度和更小的殘差。

圖14 常規(guī)逆時(shí)偏移與優(yōu)化時(shí)空域最小二乘逆時(shí)偏移的單道成像結(jié)果(x=3300m處反射系數(shù))對(duì)比

圖15 Taylor級(jí)數(shù)展開(kāi)時(shí)空域最小二乘逆時(shí)偏移與優(yōu)化時(shí)空域最小二乘逆時(shí)偏移的

圖16 收斂曲線

6 結(jié)論

本文將基于優(yōu)化時(shí)空域頻散關(guān)系的有限差分法推廣應(yīng)用于聲波最小二乘逆時(shí)偏移中。頻散分析和正演模擬表明了優(yōu)化時(shí)空域有限差分法相對(duì)于傳統(tǒng)的有限差分法和Taylor級(jí)數(shù)展開(kāi)時(shí)空域有限差分法，能夠更有效地壓制數(shù)值頻散，提高數(shù)值模擬精度。混合吸收邊界條件能較好地壓制邊界反射。優(yōu)化的最小二乘逆時(shí)偏移計(jì)算流程能夠減少波動(dòng)方程的求解次數(shù)，提高計(jì)算效率。通過(guò)對(duì)逆時(shí)偏移、Taylor時(shí)空域最小二乘逆時(shí)偏移和優(yōu)化時(shí)空域最小二乘逆時(shí)偏移的對(duì)比分析，表明最小二乘逆時(shí)偏移成像精度高、保幅性好，基于優(yōu)化時(shí)空域頻散關(guān)系的最小二乘逆時(shí)偏移能夠?qū)θ醴瓷浣缑婧透咚賷A層更好成像，能夠獲得更高精度的成像結(jié)果和提高收斂速度。