VTI介質修正聲學近似qP波波動方程與模擬

梁 鍇,陳浩然,孫上饒

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

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

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

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

1 VTI介質修正聲學近似qP波波動方程

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

(1)

(2)

(3)

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

針對精確頻散關系表達式復雜、難以應用的問題,Alkhalifah[3]提出了著名的聲學假設近似,令VS0=0(qSV波沿VTI介質對稱軸相速度為0),得到qP波聲學近似頻散關系方程:

(4)

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

(5)

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

(6)

將VS0進行如下近似:

(7)

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

(8)

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

(9)

其中,ωPma為修正聲學假設近似的qP波圓頻率。

(10)

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

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

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

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

根據有限差分基本原理,將空間導數和時間導數用差分近似表征為:

(11)

對VTI介質修正聲學近似qP波波動方程采用混合有限差分/偽譜法求解,最終的差分格式為:

(12)

其中,G表示用偽譜法計算的非橢圓項,即

(13)

3 數值示例

3.1 頻散關系分析

為了驗證VTI介質修正聲學近似qP波波動方程的正確性和適用性,進行近似qP波頻散關系的精度分析。設計了兩個均勻VTI介質模型A和B,其相同參數為VP0=3 000 m/s,VS0=1 500 m/s,其各向異性參數ε、δ不同,如表1所示。其中,模型A中ε>δ,模型B中ε<δ。

表1 模型各向異性參數

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

圖1 模型A精確和近似頻散關系及其相對誤差

3.2 正演模擬

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

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

圖3 模型A的波場快照

圖4 模型B的波場快照

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

圖5 網格點(130,120)處地震記錄

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

圖6 BP模型

圖7 波場快照

圖8 單炮記錄

4 結論

本研究對聲學近似進行了修正,推導了VTI介質修正聲學近似qP波頻散關系和波動方程,采用混合有限差分/偽譜算法實現了VTI介質純qP波正演模擬。通過理論分析和數值示例得到如下結論。

1) 基于修正聲學近似的qP波波動方程是關于時間二階導數的偏微分方程,其一組共軛解對應發散和匯聚的qP波,因此該方程的解不包含退化qSV波,而是純qP波方程。

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

3) 因為退化qSV波在ε<δ的VTI介質中存在數值溢出,導致模擬結果不穩定,而基于修正聲學近似的qP波波動方程不包含退化qSV波,因此該方程在ε≥δ和ε<δ的VTI介質中均是穩定的,能夠適用于復雜VTI介質。

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