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

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

0 引言

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

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

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

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

1 非均質聲波傳播理論

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

(1)

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

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

(2)

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

(3)

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

(4)

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

(5)

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

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

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

(6)

為加入PML吸收邊界條件，可在頻率域對坐標進行復延伸，以x方向為例

(7)

衰減函數取

(8)

式中：DPML為層寬度；x表示PML域內質點水平位置； 系數cPML為實驗選取經驗值。定義衰減系數

(9)

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

(10)

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

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

2 兩類不同網格剖分的正演算法

2.1 交錯網格法

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

(11)

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

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

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

(12)

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

將該離散公式代入式(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)

其中各系數分別如下

(15a)

(15b)

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

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

圖3 聲波變密度四階交錯網格示意圖

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

(16)

其余一階導數差分格式依此類推。將離散式(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)

其中各系數對應如下

(19)

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

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

(20)

式中：ah和am是差分系數；N表示差分階數。表1給出十階以內對應的差分系數。

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

表1 十階以內差分系數

2.2 混合網格法

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

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

(21)

式中a為加權系數。

旋轉網格框架與傳統(tǒng)網格框架下偏導數關系為

(22)

圖4 聲波變密度混合網格示意圖菱形表示45°方向質點速度

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

(23)

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

(24)

同時對質量加速度項按9點格式加權分配，可得

(25)

式中c、e為加權平均系數。結合以上兩式，按圖4所示旋轉混合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)

其中各系數對應如下

(27)

式中系數a=0.5461，c=0.6248，e=0.09381。

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

(28)

3 計算精度及效率分析

3.1 計算精度分析

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

(29)

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

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

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

3.2 數值頻散分析

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

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

(31)

式中G=λ/d=2π/(k/d)，為每個波長內的網格數。

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

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

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

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

3.3 矩陣結構分析

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

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

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

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

4 數值模擬結果及分析

4.1 層狀模型

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

圖11 層狀模型

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

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

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

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

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

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

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

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

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

4.2 Marmousi模型

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

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

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

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

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

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

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

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

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

5 結論

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

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

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

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