薛 浩 劉 洋 楊宗青

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

1 引言

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

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

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

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

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

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

(1)

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

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

(2)

采用高階差分計算空間導(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θ，其中θ為平面波的傳播角度。化簡式(6)得到

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

(7)

利用泰勒(Taylor)級數(shù)展開法展開余弦函數(shù)

(8)

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

(9)

j=1,2,…,M

(10)

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

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

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

(11)

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

(12)

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

2.2 數(shù)值頻散分析

相速度的相對誤差ξ可用來衡量有限差分法的精度，公式如下[31]

(13)

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

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

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

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

3 混合吸收邊界條件

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

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

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

圖4 吸收邊界對比圖

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

散射波場的Born正演可寫為[16]

(14)

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

方程可寫為矩陣形式

d=Lm

(15)

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

(16)

寫成矩陣形式為

m=LTd

(17)

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

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

(18)

迭代最速下降解為

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

(19)

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

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

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

Si=Lsi

(20)

Gi=Lgi

(21)

R=d-Lm

(22)

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

(23)

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

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

5 模型試算

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

5.1 均勻模型

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

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

5.2 Marmousi模型

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

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

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

圖8 Marmousi速度模型

圖9 波場快照對比圖

圖10 單炮記錄對比圖

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

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

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

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

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

圖16 收斂曲線

6 結(jié)論

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