結(jié)合泊松算法和波場(chǎng)分解成像條件的各向異性優(yōu)化純聲波方程逆時(shí)偏移

徐世剛 包乾宗 任志明 劉 洋

(①長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院地球物理系，陜西西安 710054； ②中國(guó)石油大學(xué)(北京)油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 102249)

0 引言

地下介質(zhì)普遍呈現(xiàn)各向異性，橫向各向同性(TI)介質(zhì)是地震勘探學(xué)界研究最廣的一類各向異性模型，根據(jù)其對(duì)稱軸方向可以分為VTI、HTI和TTI模型。逆時(shí)偏移技術(shù)受傾角限制小，能夠適應(yīng)復(fù)雜速度模型，獲得高精度成像剖面，因此各向異性逆時(shí)偏移研究具有重要意義。

相比于各向異性彈性波方程，各向異性聲波方程形式簡(jiǎn)單、計(jì)算量較小、波場(chǎng)特征明確。目前在各向異性波場(chǎng)模擬和偏移中主要采用偽聲波方程，該類方程的一種推導(dǎo)思路是將耦合的各向異性P-SV波頻散關(guān)系中的橫波速度設(shè)置為0簡(jiǎn)化而來(lái)[1]。在此基礎(chǔ)上，發(fā)展了多種形式的偽聲波方程[2-8]。另一種偽聲波方程的推導(dǎo)思路是基于聲學(xué)近似的胡克定律[9-10]。但是當(dāng)給定的各向異性參數(shù)不滿足偽聲波假設(shè)時(shí)，現(xiàn)有偽聲波方程在進(jìn)行P波模擬時(shí)，容易出現(xiàn)偽橫波干擾及數(shù)值不穩(wěn)定。為了徹底解決偽聲波存在的問(wèn)題，Liu等[11]解耦了VTI介質(zhì)P-SV波頻散關(guān)系，獲得單獨(dú)的P波和SV波方程。相比于顯式表達(dá)的偽聲波方程，各向異性純聲波方程包含復(fù)雜的擬微分算子，常規(guī)差分法難以直接求解。迄今為止，已經(jīng)發(fā)展了多種數(shù)值算法求解純聲波方程，如偽譜法/偽解析法[12-13]、低秩分解法[14]、迭代求解法[15]、混合求解法[16]和快速泊松算法[17-18]等。目前關(guān)于純聲波方程的研究主要集中于平衡計(jì)算精度和效率方面，以及在偏移成像中的進(jìn)一步應(yīng)用等[19-23]。

傳統(tǒng)逆時(shí)偏移技術(shù)對(duì)震源波場(chǎng)和檢波點(diǎn)波場(chǎng)做互相關(guān)，結(jié)果中容易出現(xiàn)較強(qiáng)的低頻噪聲。Zhang等[24]采用拉普拉斯濾波壓制低頻噪聲； Yoon等[25]、康瑋等[26]利用Poynting矢量切除大角度成像結(jié)果，從而壓制低頻噪聲； Liu等[27]提出波場(chǎng)分解構(gòu)建逆時(shí)偏移成像條件，選擇不同波場(chǎng)做互相關(guān)，能夠較好地抑制成像噪聲。在Liu等[27]的研究基礎(chǔ)上，基于Hilbert變換構(gòu)建復(fù)數(shù)波場(chǎng)的策略被用于提高計(jì)算效率[28-32]。波場(chǎng)分解成像條件已應(yīng)用于各向同性聲波和彈性波逆時(shí)偏移，考慮到該方法的優(yōu)勢(shì)，本文將其推廣到各向異性介質(zhì)逆時(shí)偏移。

作為對(duì)比，本文首先回顧了兩種常用的TTI偽聲波方程，分析了其存在的問(wèn)題； 然后利用最小二乘方法求取了TTI優(yōu)化純聲波頻散關(guān)系，并采用泊松算法和有限差分相結(jié)合的數(shù)值方案進(jìn)行求解； 最后結(jié)合介質(zhì)各向異性特征對(duì)波場(chǎng)分解成像條件進(jìn)行改進(jìn)，將其推廣至純聲波逆時(shí)偏移。相關(guān)理論和算例表明，將優(yōu)化純聲波方程和波場(chǎng)分解成像條件聯(lián)合應(yīng)用于各向異性純聲波逆時(shí)偏移，能有效壓制偽橫波及低頻噪聲，獲得高質(zhì)量成像結(jié)果。

1 理論方法

1.1 基于限制性橫波速度的TTI偽聲波方程

通過(guò)將原始P-SV波頻散關(guān)系中的橫波速度設(shè)置為0，Alkhalifah[1]推導(dǎo)了各向異性偽聲波方程，該方程能夠有效簡(jiǎn)化計(jì)算，其多種改進(jìn)版本被廣泛應(yīng)用于各向異性波場(chǎng)模擬和偏移成像[2-8]。但直接將垂向橫波速度設(shè)置為0并不能完全滿足各向異性介質(zhì)中的物理學(xué)規(guī)律，當(dāng)各向異性參數(shù)不滿足近似條件時(shí)(即橢圓近似條件)，偽聲波方程模擬容易出現(xiàn)SV波干擾及不穩(wěn)定。為此，學(xué)者們提出了不同的解決方案，其中一種常用的基于限制性橫波速度的TTI偽聲波方程形式[4]為

(1)

(2)

其中θ為對(duì)稱軸傾角。特別地，通過(guò)引入限制性橫波速度能夠有效提高穩(wěn)定性。當(dāng)vSz=0時(shí)，式(1)退化為傳統(tǒng)偽聲波方程，形式為

(3)

采用有限差分即可求解式(1)和式(3)。圖1為兩個(gè)方程計(jì)算的雙層TTI模型(表1)的波場(chǎng)切片，其中上、下層介質(zhì)厚度均為1500m。

表1 雙層TTI模型參數(shù)

由圖1可以看出：當(dāng)介質(zhì)具有橢圓各向異性特征(ε=δ)時(shí)，兩種方程均可產(chǎn)生高精度P波響應(yīng)，如圖1左所示； 當(dāng)各向異性參數(shù)不滿足橢圓假設(shè)時(shí)(ε≠δ)，兩種方程會(huì)產(chǎn)生較強(qiáng)的偽橫波污染，如圖1中和圖1a右所示； 當(dāng)ε<δ時(shí)，橫波速度為0的偽聲波方程(式(3))模擬出現(xiàn)不穩(wěn)定，如圖1b右所示。

1.2 TTI優(yōu)化純聲波方程及其結(jié)合泊松算法和有限 差分的解法

多種改進(jìn)版本的各向異性偽聲波方程僅能在一定程度上壓制偽橫波干擾，提高穩(wěn)定性，但不能從根本上解決問(wèn)題。眾多研究表明，發(fā)展純聲波方程能夠徹底消除偽橫波污染及不穩(wěn)定。相比于顯式表達(dá)的偽聲波方程，純聲波方程包含復(fù)雜的偏微分算子，常規(guī)差分法難以直接求解?？紤]到計(jì)算精度和效率，本文給出一種改進(jìn)的、結(jié)合泊松算法和有限差分的各向異性優(yōu)化純聲波方程求解方法[17-18]。

解耦的TTI純聲波頻散關(guān)系[11]為

(4)

(5)

將式(5)代入式(4)，可得到目前最常用的各向異性純聲波頻散關(guān)系?？紤]到一階泰勒展開的精度較低，高階泰勒展開的計(jì)算量較大，為了平衡純聲波模擬時(shí)的精度和效率，本文借鑒優(yōu)化的思想[17-18]獲得一種高精度的純聲波頻散關(guān)系。對(duì)式(4)中的根號(hào)項(xiàng)做如下近似[17-18]

圖1 雙層TTI模型式(1)(a)和式(3)方程模擬的波場(chǎng)切片(b)左：模型1； 中：模型2； 右：模型3

(6)

式(6)相對(duì)于待求近似系數(shù)(f1，f2，f3)為線性方程，因此可以通過(guò)極小化兩端的誤差獲得近似系數(shù)值。構(gòu)建的目標(biāo)函數(shù)為

(7)

通過(guò)最小二乘方法求解上式，即可獲得優(yōu)化系數(shù)(f1，f2，f3)。

將式(6)代入式(4)，整理可得優(yōu)化的純聲波頻散關(guān)系

(8)

式中

(9)

特別地，當(dāng)w1=1+2ε、w2=1、w3=-2(ε-δ)時(shí)，式(8)退化為常用的一階泰勒展開純聲波頻散式。式(8)在時(shí)間—空間域可寫為[17]

(10)

引入輔助波場(chǎng)Q，將式(10)分解為[15，17]

(11)

(12)

式(11)為顯式方程，采用常規(guī)差分即可求解。式(12)為隱式方程，其計(jì)算過(guò)程涉及線性方程組求解。Chu等[15]采用線性迭代求解隱式方程，但計(jì)算量較大。為了提高效率，本文給出一種基于泊松算法的實(shí)現(xiàn)策略[17]。

將式(11)和式(12)進(jìn)一步整理為

(13)

?2Q=P

(14)

式(13)和式(14)組成了新的TTI優(yōu)化純聲波方程，其中式(13)可采用常規(guī)差分法求解。式(14)為泊松方程，基于前人研究成果，可通過(guò)泊松算法實(shí)現(xiàn)，計(jì)算步驟[17]如下。

(1)采用五點(diǎn)差分離散式(14)

Qi-1,j+Qi+1,j+Qi,j-1+Qi,j+1-4Qi,j=h2Pi,j

(15)

式中：h為網(wǎng)格間距；i、j分別為x、z方向網(wǎng)格點(diǎn)序號(hào)。

(16)

式中：Lx和Lz分別為x、z方向網(wǎng)格點(diǎn)數(shù)；n、l分別為波數(shù)域兩個(gè)方向的網(wǎng)格點(diǎn)序號(hào)。

(3)計(jì)算輔助波場(chǎng)Q的傅里葉響應(yīng)

(17)

(4)求解輔助波場(chǎng)Q

(18)

通過(guò)顯式差分離散式(13)，通過(guò)式(15)～式(18)可實(shí)現(xiàn)式(14)的求解。這樣能夠保證純聲波方程式(13)和式(14)的高效求解。

圖2對(duì)比了一階泰勒展開方程(式(5))和優(yōu)化展開方程(式(6))對(duì)根號(hào)部分的近似精度。由圖可見，與精確值(圖2a)相比，一階泰勒展開存在較大誤差(圖2b)，而最小二乘優(yōu)化展開能夠獲得更高的近似精度(圖2c)。圖3為優(yōu)化純聲波方程模擬的雙層TTI模型波場(chǎng)切片，與圖1偽聲波方程的計(jì)算結(jié)果相比，純聲波方程在不同參數(shù)條件下均能產(chǎn)生精確有效的P波響應(yīng)，同時(shí)能夠完全消除偽橫波和數(shù)值不穩(wěn)定。

1.3 改進(jìn)波場(chǎng)分解成像條件

在逆時(shí)偏移成像中，互相關(guān)成像條件具有實(shí)現(xiàn)簡(jiǎn)單、計(jì)算高效等優(yōu)點(diǎn)，得到廣泛應(yīng)用。常規(guī)互相關(guān)成像條件[34]為

(19)

式中：x為空間坐標(biāo)；I(x)為成像結(jié)果；Ps(x,t)為震源波場(chǎng)；Pr(x,t)為檢波點(diǎn)波場(chǎng)；tmax為記錄長(zhǎng)度。

常規(guī)互相關(guān)成像條件易引起較強(qiáng)低頻噪聲，原因是互相關(guān)過(guò)程中沿所有方向傳播的震源波場(chǎng)和檢波點(diǎn)波場(chǎng)均參與成像，同向傳播的波場(chǎng)之間互相關(guān)疊加，易形成較強(qiáng)的低頻噪聲。為了壓制傳統(tǒng)互相關(guān)成像過(guò)程中的低頻噪聲，除了對(duì)低頻噪聲進(jìn)行濾波處理外[24]，另一種有效策略是對(duì)參與成像的全波場(chǎng)沿不同方向進(jìn)行分解，選取傳播方向不同的分量波場(chǎng)參與成像，能夠有效減弱低頻噪聲干擾[27-32]。

以二維情況為例，對(duì)式(19)中的震源波場(chǎng)Ps(x,t)和檢波點(diǎn)波場(chǎng)Pr(x,t)進(jìn)行分解

(20)

(21)

式中上標(biāo)“↑”、“↓”、“←”、“→”分別表示上行、下行、左行和右行波場(chǎng)分量。

由式(19)～式(21)可知，震源波場(chǎng)分量與檢波點(diǎn)波場(chǎng)分量之間兩兩互相關(guān)可得16種成像剖面。其中傳播方向相反的波場(chǎng)分量互相關(guān)能夠產(chǎn)生有效的成像結(jié)果，而同向傳播的波場(chǎng)分量之間互相關(guān)主要產(chǎn)生低頻噪聲[27-32]。在前人關(guān)于波場(chǎng)分解成像條件的研究基礎(chǔ)上，本文采用一種基于Hilbert變換構(gòu)建復(fù)數(shù)波場(chǎng)的分解算法。該分解方法能夠進(jìn)行全波場(chǎng)的顯式分離，目前被廣泛用于各向同性介質(zhì)的

圖2 均勻各向異性模型的精確微分算子(a)及一階泰勒展開誤差(b)與最小二乘優(yōu)化展開誤差(c)的對(duì)比模型參數(shù)為ε=0.2，δ=0.1，θ=60°

圖3 雙層TTI模型優(yōu)化純聲波方程模擬的波場(chǎng)切片(a)模型1； (b)模型2； (c)模型3

逆時(shí)偏移[27-32]。本文將其推廣到各向異性介質(zhì)的逆時(shí)偏移，實(shí)現(xiàn)步驟如下。

(1)基于式(13)和式(14)，震源波場(chǎng)的復(fù)數(shù)外推形式可以表示為

(22)

(2)基于步驟(1)構(gòu)建的復(fù)數(shù)波場(chǎng)，以震源波場(chǎng)為例，上行波場(chǎng)可以表示為

(23)

(24)

圖4 雙層TTI模型4的波場(chǎng)切片(a)及分離的上(b)、下(c)、左(d)、右行波場(chǎng)(e)

2 模型算例

2.1 地塹模型

各向異性地塹模型的尺寸為4000m×3000m，相關(guān)參數(shù)如圖5所示。網(wǎng)格間距為10m，時(shí)間步長(zhǎng)為1ms，震源選用主頻為28Hz的Ricker子波。

考慮到橫波速度為0的偽聲波方程存在穩(wěn)定性問(wèn)題，圖6為基于限制性橫波速度的偽聲波和優(yōu)化純聲波方程模擬的波場(chǎng)切片，可以看出，相比于偽聲波方程的模擬結(jié)果(圖6a)，純聲波方程能夠有效壓

圖5 各向異性地塹模型

制偽橫波干擾，產(chǎn)生高精度P波響應(yīng)(圖6b)。

圖7為圖6b的純聲波波場(chǎng)沿不同方向分解結(jié)果，可以發(fā)現(xiàn)，上、下、左、右行波得到了有效分離。傳統(tǒng)互相關(guān)成像條件和波場(chǎng)分解成像條件的偏移結(jié)果如圖8所示。由圖8a可知，基于傳統(tǒng)互相關(guān)條件的成像結(jié)果容易產(chǎn)生較強(qiáng)低頻噪聲干擾。對(duì)比部分分量波場(chǎng)的成像結(jié)果可見同向傳播的波場(chǎng)之間互相關(guān)主要產(chǎn)生低頻噪聲，波場(chǎng)分解成像條件能夠較好地將這部分干擾剔除，反向傳播波場(chǎng)之間互相關(guān)能夠產(chǎn)生有效的成像剖面。對(duì)圖8a中的成像剖面進(jìn)行濾波處理，結(jié)果如圖8b所示。由圖可知，不同剖面的低頻噪音得到較好地去除。反向傳播波場(chǎng)成像

圖6 地塹模型不同方程模擬的波場(chǎng)切片(a)偽聲波方程(式(1))； (b)優(yōu)化純聲波方程(式(13)、式(14))

圖7 地塹模型分離的上(a)、下(b)、左(c)、右行波場(chǎng)(d)

圖8 地塹模型優(yōu)化純聲波方程不同條件成像結(jié)果(a)及其拉普拉斯濾波結(jié)果(b)左：全波場(chǎng)互相關(guān)； 中：同向傳播波場(chǎng)互相關(guān)之和； 右：反向傳播波場(chǎng)互相關(guān)之和

結(jié)果之和略優(yōu)于傳統(tǒng)全波場(chǎng)成像結(jié)果。相關(guān)結(jié)果表明，不同方向行波之間的互相關(guān)結(jié)果可以作為傳統(tǒng)方法的補(bǔ)充，成像結(jié)果精度更高。

2.2 修改的Marmousi TTI模型

修改的Marmousi TTI模型，如圖9所示。將計(jì)算區(qū)域剖分為500×350個(gè)網(wǎng)格單元。用于逆時(shí)偏移的速度場(chǎng)和各向異性參數(shù)場(chǎng)均基于五點(diǎn)平滑獲得。圖10為各向同性聲波、各向異性偽聲波和優(yōu)化純聲波的成像結(jié)果，相比之下，基于優(yōu)化純聲波方程的逆時(shí)偏移能夠獲得更高精度的成像。圖11為基于優(yōu)化純聲波方程的不同條件成像結(jié)果的對(duì)比，可見，傳統(tǒng)全波場(chǎng)參與的互相關(guān)成像，淺層存在嚴(yán)重的低頻噪聲污染。同向傳播的波場(chǎng)之間互相關(guān)主要產(chǎn)生低頻成像噪聲，波場(chǎng)分解成像條件能夠較好地將低頻噪聲去除，反向傳播波場(chǎng)之間互相關(guān)能夠產(chǎn)生高精度的成像結(jié)果。相比于經(jīng)過(guò)拉普拉斯濾波的成像剖面而言，基于波場(chǎng)分解成像條件的偏移結(jié)果不需要濾波處理，能夠較好地保證振幅信息?？傮w而言，基于波場(chǎng)分解的成像方法可以作為傳統(tǒng)方法的有效補(bǔ)充，以獲得不同展布方向地下構(gòu)造的高精度成像。

圖9 修改的Marmousi TTI模型(a)vPz； (b)ε； (c)δ； (d)θ

圖10 修改的Marmousi TTI模型不同方程的成像結(jié)果(拉普拉斯濾波后)(a)各向同性聲波方程； (b)偽聲波方程(式(1))；(c)優(yōu)化純聲波方程(式(13)、式(14))

圖11 修改的Marmousi TTI模型的優(yōu)化純聲波方程不同條件成像結(jié)果

3 結(jié)束語(yǔ)

本文首先回顧了TTI介質(zhì)中兩種常用的偽聲波方程，在各向異性參數(shù)不滿足聲學(xué)近似時(shí)，利用該方程模擬容易產(chǎn)生偽橫波干擾及不穩(wěn)定現(xiàn)象。為了解決現(xiàn)有偽聲波方程存在的問(wèn)題，應(yīng)用最小二乘方法獲得了優(yōu)化的純聲波頻散關(guān)系，在此基礎(chǔ)上采用一種泊松算法和有限差分相結(jié)合的高效精確各向異性純聲波求解算法，能夠有效避免偽橫波干擾及數(shù)值不穩(wěn)定。將各向同性介質(zhì)中基于Hilbert變換的復(fù)數(shù)波場(chǎng)分解成像條件推廣到TTI介質(zhì)優(yōu)化純聲波逆時(shí)偏移成像，選擇傳播方向相反的分量波場(chǎng)參與最終成像。理論和數(shù)值算例驗(yàn)證了本文方法能夠有效壓制各向異性逆時(shí)偏移的偽橫波干擾和低頻噪聲，獲得高精度的成像結(jié)果。