頻率—空間域非均質(zhì)聲波有限差分模擬

吳 悠 吳國忱* 李青陽 楊凌云 賈宗鋒

(①中國石油大學(華東)地球科學與技術(shù)學院，山東青島 266580；②海洋國家實驗室海洋礦產(chǎn)資源評價與探測技術(shù)功能實驗室，山東青島 266071)

0 引言

模擬地震波的傳播對理解地下介質(zhì)中的復雜波現(xiàn)象具有重要意義。有限差分法的模擬結(jié)果能反映完整的地震波場信息，且模擬過程相對簡單、高效，因而得到廣泛應用。頻率域波動方程有限差分數(shù)值模擬相對于時間域具有獨特優(yōu)勢，如可靈活地選擇計算頻段、各單頻分量相互獨立、適于多震源模擬、適于頻率相關(guān)介質(zhì)的地震波模擬等[1]。密度是地下巖石的一個重要彈性參數(shù)，對油氣藏描述具有重要意義，因此進行包括密度參數(shù)的正演模擬是必不可少的[2]?，F(xiàn)有關(guān)于聲波方程頻率域正演研究的文獻很少考慮密度的空間變化，因此本文考慮空間變密度對波傳播的影響，構(gòu)建了基于交錯網(wǎng)格的高階有限差分模擬算法的一般框架，并應用于頻率域非均質(zhì)聲波正演。

在頻率域有限差分模擬中，大規(guī)模稀疏矩陣的結(jié)構(gòu)直接影響矩陣方程的求解[3]，而矩陣結(jié)構(gòu)的復雜度又依賴于空間導數(shù)的近似方法。Luo等[4]提出頻率域聲波方程有限差分的二階交錯網(wǎng)格簡化式，該式只涉及二階波動方程中的壓力場。Pratt等[5-6]針對聲波和彈性波方程，采用中心有限差分得到聲波方程的五點格式和彈性波方程的九點格式。Jo等[7]提出的最優(yōu)化九點方法由經(jīng)典笛卡爾坐標系和45°旋轉(zhuǎn)坐標系上二階導數(shù)算子的兩種離散方法線性組合而成，這兩種離散方法的結(jié)合降低了網(wǎng)格的各向異性，進一步結(jié)合質(zhì)量加速度項的加權(quán)平均，可得到一個具有最佳各向異性和頻散特性的緊湊模板，然后利用相速度頻散最小優(yōu)化方法，得到最優(yōu)權(quán)重系數(shù)(下文將該方法稱為混合網(wǎng)格法)。Stekl等[8]將混合網(wǎng)格法推廣到彈性波方程，在彈性情況下須考慮矢量場，因而該方法更復雜。為使混合網(wǎng)格適應流體情況，只能使用旋轉(zhuǎn)網(wǎng)格框架下的離散方法。Shin等[9]對混合網(wǎng)格法做第二次擴展，即在對聲波方程模擬時運用一個25點的模板，每個波長需要的最少網(wǎng)格點數(shù)減至2.5，阻抗矩陣的帶寬保持不變，所需核心存儲單元的數(shù)量僅為原來的四分之一。

針對頻率域波動方程正演模擬，人們也進行了諸多研究。吳國忱等[10]將25點優(yōu)化差分算子應用于頻率—空間域正演，模擬彈性波在TTI介質(zhì)中的傳播。殷文等[11]做了頻率—空間域二維彈性波方程高精度25點有限差分正演模擬研究。梁鍇等[12]基于TTI介質(zhì)做了頻率—空間域qP波方程的波場傳播特性研究。吳建魯?shù)萚13]研究了上覆液相頻率域聲—彈耦合正演并推廣到全波形反演。李青陽等[14]提出準規(guī)則網(wǎng)格高階有限差分法非均值彈性波波場模擬方法。張衡等[15]發(fā)展了一種基于平均導數(shù)方法(Average-derivative method，簡稱ADM)的25點有限差分格式以實現(xiàn)聲波方程頻率域高精度正演。劉璐等[16]綜合利用加權(quán)平均算子、平均加速度項和優(yōu)化系數(shù)三種方法，提出了優(yōu)化15點差分格式。馬曉娜等[17]從頻散關(guān)系、計算效率和存儲量三方面對比分析了四種差分方法，為高精度數(shù)值模擬和聲波頻率域全波形反演提供方法選擇上的參考。范娜等[18]提出一種通用差分系數(shù)優(yōu)化方法，只要給定差分模板形式，即可直接構(gòu)造頻散方程，求取優(yōu)化系數(shù)。

基于以上頻率—空間域波動方程正演模擬研究相關(guān)方法和思想，本文構(gòu)建了基于稀疏交錯網(wǎng)格法離散頻域二階聲波變密度方程的一般框架，且提供了一種適用于該差分框架的完全匹配層(Perfectly matched layers，PML)吸收邊界條件[19]。通過經(jīng)典諧波分析比較了混合網(wǎng)格和四階交錯網(wǎng)格的頻散特性，結(jié)果表明四階交錯網(wǎng)格模板的精度略高于混合網(wǎng)格模板。比較幾種模板阻抗矩陣的結(jié)構(gòu)，混合網(wǎng)格法的計算效率顯著優(yōu)于四階交錯網(wǎng)格法，這是由于使用四階交錯網(wǎng)格法時擴大了矩陣帶寬。另外，基于其緊湊性，混合網(wǎng)格模板在CPU和RAM性能方面是最佳的。

1 非均質(zhì)聲波傳播理論

在各向同性介質(zhì)假設下，時間域聲波方程可表示為如下二階雙曲偏微分方程

(1)

式中：ρ(x，z)是介質(zhì)密度的空間分布函數(shù)；u(x，z，t)=[u1(x，z，t)，u3(x，z，t)]為質(zhì)點位移，且u1(x，z，t)和u3(x，z，t)分別表示u(x，z，t)在x、z方向的位移分量；K(x，z)為彈性介質(zhì)體積模量。

式(1)可分解為以下兩個標量方程

(2)

為降低方程階數(shù)并簡化后續(xù)離散過程，定義

(3)

式中：P(x，z，t)為聲壓場；Q(x，z，t)和S(x，z，t)分別為x、z方向質(zhì)點速度。式(3-1)對時間求一階導數(shù)，并結(jié)合式(2)和式(3)，可得時間域二維聲波方程的一階形式

(4)

對該式做傅里葉變換，得到頻域的二維聲波變密度方程

(5)

式中： i為虛數(shù)單位；ω為角頻率。針對該式即可利用交錯網(wǎng)格進行后續(xù)離散。

地震波數(shù)值模擬只能處理有限區(qū)域，為減弱人工邊界造成的虛假反射，需在模擬過程中添加邊界條件，現(xiàn)應用最廣泛的是PML吸收邊界條件，該邊界理論上對任意角度、任意頻率的波場都能完全吸收。

針對頻率域聲波方程，首先將P(x，z，ω)分解為Px(x，z，ω)和Pz(x，z，ω)，即式(5)分裂為

(6)

為加入PML吸收邊界條件，可在頻率域?qū)ψ鴺诉M行復延伸，以x方向為例

(7)

衰減函數(shù)取

(8)

式中：DPML為層寬度；x表示PML域內(nèi)質(zhì)點水平位置； 系數(shù)cPML為實驗選取經(jīng)驗值。定義衰減系數(shù)

(9)

得到具有PML吸收邊界條件的頻率—空間域二維聲波變密度方程

(10)

為說明該PML邊界對人工邊界反射的吸收效果，在簡單二維半空間介質(zhì)中進行模擬實驗(圖1)，可知該PML邊界吸收效果良好。

圖1 最佳匹配層(PML)邊界吸收效果示意圖(a)和(b)分別表示添加邊界前、后的頻率切片(50Hz)； (c)和(d)分別表示添加邊界前、后的波場快照(250ms)

2 兩類不同網(wǎng)格剖分的正演算法

2.1 交錯網(wǎng)格法

由前文易得具有PML吸收邊界條件的頻率—空間域二維聲波變密度方程

(11)

利用二階交錯網(wǎng)格做模擬，其差分格式如圖2所示。

圖2 聲波變密度二階交錯網(wǎng)格示意圖黑點為聲壓場，三角、長方形對應x、z方向質(zhì)點速度。圖3、圖4同

采用二階中心差分格式將式(11)離散化

(12)

式中：d為步長；i、j分別表示差分中心點在x、z方向的坐標。其余一階導數(shù)差分格式依此類推。

將該離散公式代入式(11)，并按Luo等[4]給出的化簡技巧(Parsimonious technique)將式中的Q、S消掉，再合并，得到二階聲壓差分方程

(13)

對該式按圖2所示5點差分格式整理，可得

C1Pi,j+C2Pi-1,j+C3Pi+1,j+

C4Pi,j-1+C5Pi,j+1=-Si,j

(14)

其中各系數(shù)分別如下

(15a)

(15b)

由以上系數(shù)點陣構(gòu)造阻抗矩陣，并求解線性方程組，即可進行頻率域聲波變密度波場模擬。

為降低數(shù)值頻散，提高模擬精度，在利用二階交錯網(wǎng)格做模擬的基礎(chǔ)上，自然想到通過提高差分階數(shù)的方法達到此目的，四階交錯網(wǎng)格差分格式如圖3所示。

圖3 聲波變密度四階交錯網(wǎng)格示意圖

采用四階交錯網(wǎng)格差分格式將式(11)離散化

(16)

其余一階導數(shù)差分格式依此類推。將離散式(16)代入式(11)并按Luo等[4]的方法將式中的Q、S消掉，合并得到四階聲壓差分方程

(17)

對該式按圖3所示13點差分格式整理可得

C1Pi,j+C2Pi-1,j+C3Pi+1,j+C4Pi,j-1+C5Pi,j+1+C6Pi-2,j+C7Pi+2,j+C8Pi,j-2+

C9Pi,j+2+C10Pi-3,j+C11Pi+3,j+C12Pi,j-3+C13Pi,j+3=-Si,j

(18)

其中各系數(shù)對應如下

(19)

由以上系數(shù)點陣構(gòu)造阻抗矩陣，并求解線性方程組，即可進行頻率域聲波變密度波場模擬。

進而，由前述二階(式(13))、四階(式(17))聲壓差分方程，可推廣至交錯網(wǎng)格任意階聲壓差分方程

(20)

式中：ah和am是差分系數(shù)；N表示差分階數(shù)。表1給出十階以內(nèi)對應的差分系數(shù)。

不同差分階數(shù)的方程對應不同差分系數(shù)，將差分系數(shù)代入式(20)并按相應差分格式構(gòu)造系數(shù)點陣，進一步構(gòu)造阻抗矩陣并求解，即可進行頻率—空間域聲波變密度方程交錯網(wǎng)格高階有限差分模擬。

表1 十階以內(nèi)差分系數(shù)

2.2 混合網(wǎng)格法

Jo等[7]提出一種最優(yōu)化9點有限差分格式，該方法將每個波長需要的最小網(wǎng)格點數(shù)減至約4個，它由經(jīng)典笛卡爾坐標系和45°旋轉(zhuǎn)坐標系上二階導數(shù)算子的兩個離散方法線性組合而成。據(jù)此，本文總結(jié)一種混合網(wǎng)格有限差分法，利用兩種二階交錯網(wǎng)格的混合，實現(xiàn)頻率域非均勻介質(zhì)聲波場的高精度數(shù)值模擬。混合網(wǎng)格差分格式如圖4所示。

將傳統(tǒng)二階網(wǎng)格與旋轉(zhuǎn)45°網(wǎng)格結(jié)合，得到具有PML吸收邊界條件的混合網(wǎng)格頻率域二維聲波變密度方程

(21)

式中a為加權(quán)系數(shù)。

旋轉(zhuǎn)網(wǎng)格框架與傳統(tǒng)網(wǎng)格框架下偏導數(shù)關(guān)系為

(22)

圖4 聲波變密度混合網(wǎng)格示意圖菱形表示45°方向質(zhì)點速度

采用二階中心差分格式將式(21)離散化

(23)

其余一階導數(shù)差分格式依此類推。將該離散式代入式(21)，并按Luo等[4]的方法將式中的Q、S消掉，合并得到混合四階聲壓差分方程

(24)

同時對質(zhì)量加速度項按9點格式加權(quán)分配，可得

(25)

式中c、e為加權(quán)平均系數(shù)。結(jié)合以上兩式，按圖4所示旋轉(zhuǎn)混合9點差分格式整理，可得

C1Pi,j+C2Pi-1,j+C3Pi+1,j+C4Pi,j-1+C5Pi,j+1+

R1Pi-1,j-1+R2Pi+1,j-1+R3Pi-1,j+1+R4Pi+1,j+1=-Si,j

(26)

其中各系數(shù)對應如下

(27)

式中系數(shù)a=0.5461，c=0.6248，e=0.09381。

優(yōu)化求取方法參照Jo等[7]方法，由以上系數(shù)點陣構(gòu)造阻抗矩陣，并求解線性方程組，即可做頻率域聲波變密度波場模擬。進一步地，由二階混合差分方程(式(27))格式，可推廣至混合格式任意階聲壓差分方程

(28)

3 計算精度及效率分析

3.1 計算精度分析

頻率—空間域二維聲波方程解析解為

(29)

將二階交錯網(wǎng)格法和混合網(wǎng)格法的數(shù)值模擬結(jié)果與解析解結(jié)果進行對比(圖5和圖6)，可知二階交錯網(wǎng)格法和混合網(wǎng)格法的模擬精度與解析解結(jié)果總體相差不大，具有較高模擬精度，且混合網(wǎng)格法的模擬精度明顯好于二階交錯網(wǎng)格法。

圖5 二階交錯(a)、四階(b)、混合(c)三種網(wǎng)格及解析法(d)得到的頻率切片(f=30Hz)對比

圖6 二階交錯(a)、四階交錯(b)、混合(c)三種網(wǎng)格及其解析解的頻率切片抽道(第81道)對比

3.2 數(shù)值頻散分析

為有效壓制數(shù)值模擬過程中的頻散現(xiàn)象，同時提高模擬精度和效率，須通過頻散分析確定合適的模擬參數(shù)。以四階交錯網(wǎng)格差分法為例，將平面簡諧波解P=P0e-i(kx+kz)(e為自然常數(shù))代入式(17)所示四階交錯網(wǎng)格離散方程，并假設dx=dz=d，其中dx、dz分別為x和z方向的空間采樣間隔，θ為傳播角度，k為波數(shù)，α1=9/8、α2=-1/24為差分系數(shù)，離散情況下有Pi，j=P0e-idk(isinθ+jcosθ)，得到頻散方程

4α1α2cos(kdsinθ)-4α1α2cos(2kdsinθ)-

2α2α2cos(3kdsinθ)+4α2α2-2α1α1cos(kdcosθ)+

4α1α2cos(kdcosθ)-4α1α2cos(2kdcosθ)-

2α2α2cos(3kdcosθ)]

(30)

定義數(shù)值相速度和群速度分別為vph=ω/k和vgr=?ω/?k，可得

(31)

式中G=λ/d=2π/(k/d)，為每個波長內(nèi)的網(wǎng)格數(shù)。

同理，可對二階交錯網(wǎng)格法及混合網(wǎng)格法進行分析，得到如圖7～圖9所示歸一化相速度和群速度頻散曲線。由圖可知相同采樣間隔情況下四階交錯網(wǎng)格法和混合網(wǎng)格法的數(shù)值頻散情況遠優(yōu)于二階交錯網(wǎng)格。

圖7 二階交錯網(wǎng)格法求取的歸一化相速度(a)和歸一化群速度(b)頻散曲線

圖8 四階交錯網(wǎng)格法求取的歸一化相速度(a)和歸一化群速度(b)頻散曲線

圖9 混合網(wǎng)格法求取的歸一化相速度(a)和歸一化群速度(b)頻散曲線

3.3 矩陣結(jié)構(gòu)分析

在頻率域有限差分模擬過程中，大規(guī)模稀疏矩陣的結(jié)構(gòu)直接影響矩陣方程的求解，而矩陣結(jié)構(gòu)的復雜度與空間導數(shù)的近似方法密切相關(guān)。表2對比了三種差分格式對應的阻抗矩陣特征，其中nx為水平方向的網(wǎng)格點數(shù)。從非零條帶的分布可以看出相較于二階交錯網(wǎng)格，混合網(wǎng)格非零條帶分布復雜程度稍有增加，而四階交錯網(wǎng)格非零條帶分布復雜度大大增加，這也相應的反映在了計算時間上。

圖10對比了三種差分格式對應的阻抗矩陣示意圖，其中模型網(wǎng)格尺寸為15×15，即對應于一個225階矩陣。圖中可見，相較于二階格式矩陣，四階格式的帶寬增加了三倍，而混合網(wǎng)格格式的矩陣帶寬僅稍有增加，這極大降低求解相應線性方程組的計算量和內(nèi)存耗用量。

表2 三種差分格式對應的阻抗矩陣特征

圖10 二階(a)、四階(b)、混合(c)三種差分格式對應的阻抗矩陣示意圖

4 數(shù)值模擬結(jié)果及分析

4.1 層狀模型

構(gòu)建一個200×200網(wǎng)格單元的簡單層狀模型(圖11)，震源在(1010m，210m)處，空間步長統(tǒng)一為10m，震源采用主頻為20Hz的雷克子波(圖11)。

圖11 層狀模型

圖12為層狀介質(zhì)模型得到的30Hz頻率切片，可觀察到速度界面和密度界面反射現(xiàn)象。

圖13為頻率切片經(jīng)傅里葉逆變換得到的時間域t=500ms波場快照，速度層界面和密度層界面反射波明顯，從中可見二階交錯網(wǎng)格法得到的時間域波場存在明顯頻散現(xiàn)象，模擬精度低，相比之下四階交錯網(wǎng)格法和混合網(wǎng)格法得到的波場快照顯示數(shù)值頻散得到有效壓制，模擬精度明顯提高。

圖14為層狀介質(zhì)模型得到的共中心點道集，其中直達波、速度層反射波和密度層反射波同相軸明顯，二階交錯網(wǎng)格法得到的道集中數(shù)值頻散干擾明顯，而四階交錯網(wǎng)格和混合網(wǎng)格法則有效壓制了數(shù)值頻散。

圖15是對三個共中心點道集抽取第81道記錄對比，結(jié)果顯示四階交錯網(wǎng)格法和混合網(wǎng)格法模擬精度接近且遠高于二階交錯網(wǎng)格法，二階交錯網(wǎng)格法在反射波處存在明顯數(shù)值頻散現(xiàn)象。

因此，對層狀模型的正演模擬測試說明了方法的精度和準確性。

圖12 二階交錯(a)、四階交錯(b)、混合(c)三種網(wǎng)格對應的層狀模型f=30Hz頻率切片

圖13 二階交錯(a)、四階交錯(b)、混合(c)三種網(wǎng)格對應的層狀模型t=500ms波場快照

圖14 二階交錯(a)、四階交錯(b)、混合(c)三種網(wǎng)格對應的層狀模型地震記錄

圖15 二階交錯(紅)、四階交錯(黑)、混合(綠)三種網(wǎng)格對應的層狀模型地震記錄抽道對比(第81道)

4.2 Marmousi模型

為進一步驗證和說明本文方法對復雜模型的適用性和穩(wěn)定性，采用復雜Marmousi模型(圖16)進行測試。模型尺寸為500×250，網(wǎng)格間距為10m，震源采用主頻為20Hz的雷克子波，震源激發(fā)點位于坐標(2510m，251m)處，檢波器置于模型表面。

圖17為由該模型得到的30Hz頻率切片，從中可觀察到界面反射現(xiàn)象。

圖18為頻率片經(jīng)過傅里葉逆變換得到的時間域t=500ms波場快照，其界面反射波明顯，且可見二階交錯網(wǎng)格法得到的時間域波場存在明顯的頻散現(xiàn)象，模擬精度低，相比之下四階交錯網(wǎng)格法和混合網(wǎng)格法得到的波場快照顯示數(shù)值頻散得到有效壓制，模擬精度明顯提高。

圖19為由模型得到的共中心點道集，其中直達波、反射波同相軸明顯，二階交錯網(wǎng)格法所得道集中數(shù)值頻散干擾明顯，而四階交錯網(wǎng)格法和混合網(wǎng)格法則有效壓制了數(shù)值頻散。

總之，對Marmousi模型進行正演模擬測試驗證了方法的穩(wěn)定性和可靠性。

圖16 縱波速度(a)和密度(b)表征的Marmousi模型

圖17 二階交錯(a)、四階交錯(b)、混合(c)三種網(wǎng)格對應的Marmousi模型f=30Hz頻率切片

圖18 二階交錯(a)、四階交錯(b)、混合(c)三種網(wǎng)格對應的Marmousi模型t=500ms波場快照

圖19 二階交錯(a)、四階交錯(b)、混合(c)三種網(wǎng)格對應的Marmousi模型共中心點道集記錄

5 結(jié)論

針對頻率—空間域非均勻聲介質(zhì)中地震波有限差分正演模擬問題，本文運用交錯網(wǎng)格思想，分別用笛卡爾坐標系提高差分階數(shù)和結(jié)合45°旋轉(zhuǎn)坐標系提高差分階數(shù)兩種方式提高模擬精度，導出兩類方法的有限差分公式，對比分析了兩類方法的計算精度和效率，并通過模型測試驗證了方法的精度和穩(wěn)定性，得到如下認識。

(1)考慮密度參數(shù)在油氣地球物理勘探中的重要性，本文研究頻率—空間域非均勻聲介質(zhì)地震波有限差分正演模擬，并構(gòu)建了分別應用兩種方式提高正演模擬精度的一般框架，得到了任意階數(shù)下兩種方法的一般差分公式，對后續(xù)正演和反演研究都具有現(xiàn)實意義。

(2)對交錯網(wǎng)格和混合網(wǎng)格模型的數(shù)值頻散分析表明，混合網(wǎng)格的相速度和群速度的數(shù)值離散度都小于二階交錯網(wǎng)格，因為混合網(wǎng)格法將笛卡爾坐標系和旋轉(zhuǎn)坐標系上的兩種二階交錯網(wǎng)格模板線性組合，同時對質(zhì)量加速度項進行了加權(quán)平均，加強了兩種交錯網(wǎng)格模板之間的耦合。

(3)幾種差分方法的計算效率由矩陣結(jié)構(gòu)復雜度決定，也就是由離散模板幾何的緊湊性控制，構(gòu)建網(wǎng)格時涉及的周圍點越少越好。二階的混合網(wǎng)格模板需9個網(wǎng)格點，四階近似的交錯網(wǎng)格模板需13個網(wǎng)格點，因此用于矩陣分解的內(nèi)存需求和浮點運算顯著增加。在內(nèi)存占用和計算效率兩方面，混合網(wǎng)格策略對于二維頻率域有限差分建模均展現(xiàn)出明顯的優(yōu)越性。