VTI介質(zhì)修正聲學(xué)近似qP波波動(dòng)方程與模擬

梁 鍇,陳浩然,孫上饒

(中國石油大學(xué)(華東) 地球科學(xué)與技術(shù)學(xué)院,山東 青島 266580)

各向異性介質(zhì)中縱橫波通常是耦合在一起傳播的,因此各向異性彈性波波場解耦在地震成像和全波形反演中十分重要[1-2]。對(duì)應(yīng)具有垂直對(duì)稱軸的橫向各向同性(transverse isotropy with a vertical axis of symmetry,VTI)介質(zhì)是實(shí)際地質(zhì)條件的有效簡化,得到了廣泛的關(guān)注。Alkhalifah[3]提出VTI介質(zhì)聲學(xué)近似方法,即通過設(shè)置垂向準(zhǔn)橫波(qSV波)速度(VS0)為0來解耦qP波場。然而,在聲學(xué)假設(shè)下,準(zhǔn)縱波(qP波)方程有兩組共軛解,分別對(duì)應(yīng)qP波和退化qSV波。后者是不需要的干擾波,并且各向異性參數(shù)ε<δ時(shí),該方法隨著迭代次數(shù)的增加會(huì)出現(xiàn)嚴(yán)重的數(shù)值溢出現(xiàn)象。

為了克服退化qSV波的干擾以及ε<δ時(shí)的不穩(wěn)定現(xiàn)象,許多學(xué)者進(jìn)行了各向異性介質(zhì)改進(jìn)或優(yōu)化的純qP波方程的研究[4-11]。然而,這些方程大多含有分?jǐn)?shù)階非橢圓項(xiàng),導(dǎo)致模擬的計(jì)算成本較高。利用有限差分(finite difference,FD)方法求解這類純波波動(dòng)方程具有高效、適應(yīng)復(fù)雜介質(zhì)等特點(diǎn)[12],相比之下利用偽譜(pseudo-spectral,PS)法求解這類方程空間導(dǎo)數(shù)的精度更高[13],但對(duì)于復(fù)雜介質(zhì)則會(huì)出現(xiàn)嚴(yán)重的數(shù)值發(fā)散。為了降低計(jì)算成本、提高計(jì)算精度,有限差分和偽譜混合的正演方法[14]以及低秩近似方法[15]等多種數(shù)值模擬方法相繼被提出。

本研究將垂向qSV波速度VS0考慮成波數(shù)kx、kz和各向異性參數(shù)ε、δ的函數(shù),推導(dǎo)了VTI介質(zhì)修正聲學(xué)近似qP波頻散關(guān)系,通過傅里葉反變換得到修正聲學(xué)近似的qP波波動(dòng)方程,然后采用混合有限差分/偽譜算法實(shí)現(xiàn)了數(shù)值模擬,最后進(jìn)行了精度和近似的頻散關(guān)系分析,并給出了均勻介質(zhì)模型和復(fù)雜介質(zhì)模型的正演模擬數(shù)值示例。

1 VTI介質(zhì)修正聲學(xué)近似qP波波動(dòng)方程

各向異性介質(zhì)的Christoffel方程通常用于各向異性介質(zhì)彈性波相速度或頻散關(guān)系表征。根據(jù)VTI介質(zhì)Christoffel方程可以確定VTI介質(zhì)qP波和qSV波精確頻散關(guān)系為:

(1)

(2)

(3)

其中:ωP、ωSV分別為qP波、qSV波的圓頻率;kx=ksinθ和kz=kcosθ分別為沿x和z方向的波數(shù),k為波數(shù);θ為傳播方向與z軸的夾角;VP0和VS0分別為qP波和qSV波沿VTI介質(zhì)對(duì)稱軸方向(z軸方向)的相速度;ε和δ為Thomsen參數(shù)[16]。

針對(duì)精確頻散關(guān)系表達(dá)式復(fù)雜、難以應(yīng)用的問題,Alkhalifah[3]提出了著名的聲學(xué)假設(shè)近似,令VS0=0(qSV波沿VTI介質(zhì)對(duì)稱軸相速度為0),得到qP波聲學(xué)近似頻散關(guān)系方程:

(4)

其中,ωPa為聲學(xué)假設(shè)近似qP波圓頻率。傳統(tǒng)聲學(xué)近似頻散關(guān)系形式簡單、精度較高,但是對(duì)應(yīng)的波動(dòng)方程存在退化qSV波等問題。本研究對(duì)該方法進(jìn)行了修正,不是令VS0直接為0,而是將VS0考慮成波數(shù)kx、kz和各向異性參數(shù)ε、δ的函數(shù)。首先根據(jù)式(2)得到VS0和ωSV的關(guān)系為:

(5)

將式(5)代入式(1)可得:

(6)

將VS0進(jìn)行如下近似:

(7)

將式(7)代入式(5),此時(shí)ωSV可表示為:

(8)

將式(8)代入式(6),可以得到VTI修正聲學(xué)近似qP波頻散關(guān)系:

(9)

其中,ωPma為修正聲學(xué)假設(shè)近似的qP波圓頻率。

(10)

其中:P=P(x,z,t)為時(shí)空域的qP波波場;F、F-1分別表示傅里葉正變換和反變換。

VTI介質(zhì)修正聲學(xué)近似qP波波動(dòng)方程(10)包含橢圓項(xiàng)和非橢圓項(xiàng)兩部分。當(dāng)ε=δ,即橢圓各向異性介質(zhì)時(shí),非橢圓項(xiàng)為0,該方程退化為橢圓波動(dòng)方程,說明該方程在橢圓各向異性介質(zhì)中是完全精確的。

2 混合有限差分/偽譜法正演模擬算法

對(duì)VTI介質(zhì)修正聲學(xué)近似qP波波動(dòng)方程(10)進(jìn)行數(shù)值模擬,由于該方程包含橢圓項(xiàng)和非橢圓項(xiàng)兩部分,因此采用混合有限差分/偽譜法求解,即橢圓項(xiàng)采用有限差分算法求解,而非橢圓項(xiàng)包含波數(shù)kx、kz的高階項(xiàng),不適合直接采用有限差分算法求解,則采用偽譜法進(jìn)行求解。

根據(jù)有限差分基本原理,將空間導(dǎo)數(shù)和時(shí)間導(dǎo)數(shù)用差分近似表征為:

(11)

對(duì)VTI介質(zhì)修正聲學(xué)近似qP波波動(dòng)方程采用混合有限差分/偽譜法求解,最終的差分格式為:

(12)

其中,G表示用偽譜法計(jì)算的非橢圓項(xiàng),即

(13)

3 數(shù)值示例

3.1 頻散關(guān)系分析

為了驗(yàn)證VTI介質(zhì)修正聲學(xué)近似qP波波動(dòng)方程的正確性和適用性,進(jìn)行近似qP波頻散關(guān)系的精度分析。設(shè)計(jì)了兩個(gè)均勻VTI介質(zhì)模型A和B,其相同參數(shù)為VP0=3 000 m/s,VS0=1 500 m/s,其各向異性參數(shù)ε、δ不同,如表1所示。其中,模型A中ε>δ,模型B中ε<δ。

表1 模型各向異性參數(shù)

利用式(1)、式(4)和式(9)分別計(jì)算精確qP波頻散關(guān)系ωP、聲學(xué)近似qP波頻散關(guān)系ωPa和修正聲學(xué)近似頻散關(guān)系ωPma及其相對(duì)誤差,結(jié)果如圖1和圖2所示。其中,圖1(a)和圖2(a)為頻散關(guān)系曲線,圖1(b)和圖2(b)為相對(duì)誤差曲線,藍(lán)色實(shí)線為精確qP波頻散關(guān)系,紅色虛線為聲學(xué)近似qP波頻散關(guān)系及其誤差,黑色虛線為修正聲學(xué)近似頻散關(guān)系及其誤差。由圖1和圖2發(fā)現(xiàn),聲學(xué)近似和修正聲學(xué)近似的頻散關(guān)系與精確頻散關(guān)系數(shù)值十分接近,本例中相對(duì)誤差均小于0.25%,說明兩種近似均具有很高的精度;并且對(duì)于ε>δ的VTI介質(zhì)和ε<δ的VTI介質(zhì),修正聲學(xué)近似頻散關(guān)系的最大相對(duì)誤差數(shù)值以及平均相對(duì)誤差數(shù)值均小于聲學(xué)近似的誤差,說明修正聲學(xué)近似的整體精度高于聲學(xué)近似的精度。

圖1 模型A精確和近似頻散關(guān)系及其相對(duì)誤差

3.2 正演模擬

為了進(jìn)行波場對(duì)比,利用有限差分法求解VTI介質(zhì)彈性波波動(dòng)方程和聲學(xué)近似qP波波動(dòng)方程,利用有限差分和偽譜混合法求解修正聲學(xué)近似qP波波動(dòng)方程式(式(12)和式(13))。首先考慮表1中的兩種均勻VTI介質(zhì)模型,網(wǎng)格大小為301×301,橫向和縱向網(wǎng)格間距均為10 m,采樣間隔為1 ms,震源位于模型中心,采用主頻20 Hz的雷克子波。

波場快照如圖3和圖4所示,其中圖3(a)和圖4(a)為彈性波波場,圖3(b)和圖4(b)為基于聲學(xué)近似模擬的波場,圖3(c)和圖4(c)為基于修正聲學(xué)近似模擬的qP波波場。在圖3中,聲學(xué)近似和修正聲學(xué)近似模擬的qP波波場與彈性波模擬的qP波波場均吻合較好,說明兩種方法均適用于ε>δ的模型。但是聲學(xué)近似模擬的波場中存在退化的qSV波波場,說明聲學(xué)近似模擬的波場不是純qP波;修正聲學(xué)近似模擬的波場中不存在退化的qSV波波場,說明修正聲學(xué)近似模擬的波場是純qP波。在圖4中,聲學(xué)近似模擬的波場存在數(shù)值溢出,使模擬的qP波波場淹沒在數(shù)值溢出的干擾中,即聲學(xué)近似無法適用于ε<δ的模型。而修正聲學(xué)近似模擬的qP波波場與彈性波模擬的qP波波場也吻合較好,且不存在數(shù)值溢出現(xiàn)象。該數(shù)值示例表明修正聲學(xué)近似qP波波動(dòng)方程刻畫了純qP波波場,且既適用于ε≥δ的VTI介質(zhì),也適用于ε<δ的VTI介質(zhì)。

圖3 模型A的波場快照

圖4 模型B的波場快照

網(wǎng)格點(diǎn)(130,120)處的地震記錄如圖5所示,其中圖5(a)為模型A模擬結(jié)果,圖5(b)為模型B模擬結(jié)果,藍(lán)色實(shí)線為彈性波記錄,紅色點(diǎn)虛線為聲學(xué)近似記錄,黑色虛線為修正聲學(xué)近似記錄。由圖5可見,彈性波記錄中包含qP波和qSV波,聲學(xué)近似記錄中包含qP波和退化qSV波,而修正聲學(xué)近似記錄中只包含qP波,并且與彈性波模擬中的qP波記錄吻合較好。

圖5 網(wǎng)格點(diǎn)(130,120)處地震記錄

將本研究方法應(yīng)用于復(fù)雜介質(zhì)模型(BP模型)。網(wǎng)格大小為601×601,橫向和縱向網(wǎng)格間距均為10 m,采樣間隔為1 ms。震源位于圖6星號(hào)所示處,采用主頻20 Hz的雷克子波。BP模型如圖6所示,圖6(a)為VP0模型,圖6(b)為ε模型,圖6(c)為δ模型。波場快照如圖7所示,其中圖7(a)為彈性波波場,圖7(b)為聲學(xué)近似模擬的波場,圖7(c)為基于修正聲學(xué)近似模擬的qP波波場。圖7表明,復(fù)雜介質(zhì)中聲學(xué)近似模擬的qP波波場與彈性波模擬的qP波波場運(yùn)動(dòng)學(xué)規(guī)律一致,波前面吻合較好,但是聲學(xué)近似模擬波場中可能存在退化qSV波(圖7(b)中白色箭頭所示),會(huì)增加純qP波處理和解釋的難度。而修正聲學(xué)近似模擬的qP波波場的運(yùn)動(dòng)學(xué)規(guī)律與彈性波模擬的qP波波場一致,二者波前面吻合較好,且該方法模擬結(jié)果中不存在退化qSV波,說明該方法在復(fù)雜VTI介質(zhì)中也能很好地實(shí)現(xiàn)純qP波正演模擬。圖8為震源深度處的單炮記錄,其中圖8(a)為彈性波記錄,圖8(b)為聲學(xué)近似記錄,圖8(c)為修正聲學(xué)近似記錄。由圖8可以看出,彈性波記錄中qSV波比較強(qiáng);聲學(xué)近似記錄中qP波比較明顯,但是仍然存在退化qSV波的干擾;修正聲學(xué)近似記錄中只有qP波,并且與彈性波模擬中的qP波記錄基本一致。

圖6 BP模型

圖7 波場快照

圖8 單炮記錄

4 結(jié)論

本研究對(duì)聲學(xué)近似進(jìn)行了修正,推導(dǎo)了VTI介質(zhì)修正聲學(xué)近似qP波頻散關(guān)系和波動(dòng)方程,采用混合有限差分/偽譜算法實(shí)現(xiàn)了VTI介質(zhì)純qP波正演模擬。通過理論分析和數(shù)值示例得到如下結(jié)論。

1) 基于修正聲學(xué)近似的qP波波動(dòng)方程是關(guān)于時(shí)間二階導(dǎo)數(shù)的偏微分方程,其一組共軛解對(duì)應(yīng)發(fā)散和匯聚的qP波,因此該方程的解不包含退化qSV波,而是純qP波方程。

2) 基于修正聲學(xué)近似的qP波波動(dòng)方程包含橢圓項(xiàng)和非橢圓項(xiàng)兩部分,非橢圓項(xiàng)在橢圓各向異性介質(zhì)中等于零,因此該方程在橢圓各向異性介質(zhì)中是完全精確的。

3) 因?yàn)橥嘶痲SV波在ε<δ的VTI介質(zhì)中存在數(shù)值溢出,導(dǎo)致模擬結(jié)果不穩(wěn)定,而基于修正聲學(xué)近似的qP波波動(dòng)方程不包含退化qSV波,因此該方程在ε≥δ和ε<δ的VTI介質(zhì)中均是穩(wěn)定的,能夠適用于復(fù)雜VTI介質(zhì)。

4) 直接用有限差分法求解波動(dòng)方程中高階導(dǎo)數(shù)是比較困難的,而混合有限差分/偽譜算法能夠較好地解決波數(shù)kx、kz高階項(xiàng)的計(jì)算問題。