基于改進NAD方法的頻率域聲波逆時偏移

郎超，劉少林，楊曉婷，徐錫偉

1 北京信息科技大學理學院，北京 100192 2 應急管理部國家自然災害防治研究院，北京 100085

0 引言

隨著我國社會經濟的快速發展，對石油和天然氣資源的需求也與日俱增，大力提升油氣勘探水平已成為當務之急，對國家資源能源儲備與安全戰略具有重要意義(李幼銘和劉洪，2010).疊前逆時偏移是油氣勘探的有效途徑之一，它以模擬地震波的傳播規律為基礎，并結合適當的偏移條件來獲得地下介質的成像模型(底青云和王妙月,1997).偏移成像方法的研究可追溯至20世紀70年代(Claerbout,1971)，受當時計算條件的限制，早期的地震偏移成像主要基于射線理論(Schneider,1978).由于射線模型只能粗略刻畫實際地震波的物理屬性，這類方法較難對復雜地質結構進行高分辨率成像(Liu et al.,2011).為進一步提高成像結果的質量，基于雙程波方程的逆時偏移成像方法引起研究者們的關注(Baysal et al.,1983；陳生昌和周華敏,2018;Lu et al.,2020)，然而求解波動方程需要消耗較多的計算資源，相應成像過程的計算效率有待提升(黃金強和李振春,2019).因此，研究高效的波動方程逆時偏移成像方法對于提高油氣勘探水平至關重要(杜啟振和秦童,2009).

相比于傳統時間域偏移成像方法，頻率域逆時偏移的優勢主要體現在以下方面：(1)當時間變量經過Fourier變換至頻率域時，波動方程關于頻率取值相互解耦，此時具有良好的并行特性(Chen and Cao,2016)；(2)不同頻率取值下的地震波場可獨立計算，避免了數值誤差的累積與傳播(Pratt et al.,1998)；(3)借助復速度或者復頻率，可更加方便地刻畫頻散或衰減效應(Song and Williamson,1995)；(4)只需選取適量的頻率點進行成像計算，效率較高(Sirgue and Pratt,2004).由此可見，對于頻率域逆時偏移成像方法的研究尤為必要(Shin et al.,2001).

在逆時偏移成像中，絕大部分計算操作集中于模擬地下介質中地震波的傳播情況(正演過程).因此，偏移效率以及成像結果的質量取決正演算法的有效性 (嚴紅勇和劉洋,2013;Qiu et al.,2020).目前地球物理領域常用的正演數值方法包括有限差分方法(Alford et al.,1974;Moczo et al.,2007;劉少林等,2013)、有限元方法(Marfurt,1984;Liu et al.,2014;劉少林等，2014)、偽譜法(Fornberg,1975;Kosloff and Baysal,1982;黃繼偉和劉洪，2020)、譜元法(Komatitsch et al.,2005;Liu et al.,2017a)等.區別于時間域波動方程的數值離散，選取頻率域數值方法須在保證離散格式準確性的同時，還要兼顧到求解離散后所得大型線性代數方程組的計算代價可控(Pratt and Worthington,1990)，這是大尺度區域模型頻率域計算所面臨的主要挑戰之一(Ernst and Gander,2012).求解該線性方程組的計算代價主要依賴于其系數矩陣(也稱為“阻抗矩陣”(Pratt et al.,1998))的稀疏性與條件數，這通常由離散算子的結構與性質所決定，從而頻率域波動方程的數值離散應當盡可能避免復雜格式以及線性方程組的病態程度過高(Lang and Yang,2016).有限差分方法具有原理直觀、格式構造簡單、易于實現、方便并行等優點，有利于頻率域波動方程的數值離散.這類方法已廣泛應用于頻率域地震波場模擬、逆時偏移成像以及全波形反演等方面(Yang et al.,2006;Xu et al.,2010;Lang et al.,2020b).然而普通差分格式在粗網格下較難對復雜地質結構(例如：強間斷面、孔縫等)中地震波的傳播情況進行準確模擬(Liu et al.,2017b)，容易引起嚴重的數值頻散假象(Moczo et al.,2000)，進而導致偏移成像結果的分辨率受損，此時若加密離散網格，問題規模也隨之增加，因此計算效率偏低.

針對普通有限差分方法所面臨的難點，一些學者提出利用極小化數值速度(或者數值波數)與理論速度(或理論波數)之間相對誤差的頻散關系式，在不顯著增加差分算子結構復雜性的前提下調整網格差分模板，可提升相應差分格式壓制數值頻散的能力(Zhang and Yao,2013;Jing et al.,2017;Chen and Cao,2018)；另一方面，楊頂輝等(Yang et al.,2003)提出采用波場及其梯度值的線性組合來逼近高階偏導數項，構造近似解析離散化(Nearly Analytic Discrete，NAD)方法，該方法不僅繼承了普通有限差分方法格式構造的簡易性，相應的離散算子還具有更強的緊致性與更高的數值精度(Yang et al.,2004)；通過這兩種途徑所構造的新型差分方法，均可在粗網格下壓制數值頻散，明顯提高正演模擬的計算效率.

本文根據頻率域計算的特點將以上兩種思路相結合，發展一種有助于大尺度區域模型頻率域波場模擬與偏移成像的高效算法.具體地，在固定NAD網格差分模板的前提下通過優化模板系數得到改進NAD方法，相比于標準NAD格式，改進NAD方法可緩解離散后所形成線性方程組的病態程度(降低阻抗矩陣的條件數)，頻率域正演模擬效率得到進一步提升；離散后的大型稀疏線性方程組將統一采用“不精確旋轉分塊三角預處理子”(Inexact Rotated Block Triangular Preconditioner，IRBTP)結合Krylov子空間迭代方法進行快速求解(Lang and Ren,2015)；再通過推導適當的頻率域逆時偏移成像條件，并結合基于改進NAD方法的正演計算過程，分別對凹陷模型、Sigsbee2B模型以及Marmousi模型進行偏移成像；同時也將相應結果與四階標準NAD方法以及四階普通有限差分(Ordinary Finite Difference，OFD)方法(Jo et al.,1996)進行對比，以此來驗證改進NAD方法在提高偏移效率以及成像質量方面的優勢.

1 正演算法

1.1 NAD離散過程

若將常密度介質中的二維頻率域聲波方程記為

(1)

其中Δ表示二維Laplace算子，ω=2πf為角頻率，c=c(x,z)表示波速，s為震源項，Ω表示二維計算區域.在NAD類方法的離散過程中，需要對(1)式分別關于x和z求偏導數(Yang et al.,2004)，則有

(2)

為消除人工截斷計算區域所產生的虛假反射波，需對Ω的邊界進行處理.本文采用“完美匹配層(Perfectly Matched Layer)”吸收邊界條件(Berenger,1994).應該指出的是，在頻率域施加吸收邊界條件來刻畫波的衰減特性具有天然的便捷性，可避免因反卷積運算所產生的誤差(Lang and Yang,2016).首先引入復坐標(以水平方向為例)

(3)

(4a)

(4b)

(4c)

在本文的討論中，設置x與z方向上的空間離散步長相等，統一記作h.若將(4)式中x方向上高階偏導數在(i,j)點處進行離散的NAD網格差分模板表示為

(5)

z方向討論類似，則(4)式中的高階混合偏導數可借助(5)式以及方向導數的性質進行推導(Tong et al.,2011)

(6)

(7)

以及

(8)

需要指出的是，當網格模板系數取作

(9)

此時對應于四階標準NAD方法(Yang et al.,2006).

采用所有高階偏導數的NAD網格差分模板(5)—(8)對(4)式進行離散(具體過程可見附錄A)，再將二維區域Ω中的網格節點按照從上至下逐行依次排列(每行從左至右)的規則且同一節點處的波場及其梯度項相鄰放置，即可形成線性方程組

Cu=b,(10)

矩陣C的稀疏分塊結構如下：

(11)

1.2 改進NAD方法

(12)

觀察(12)式，阻抗矩陣的主對角元素取值僅依賴于網格差分模板系數a0與d0(忽略高階小量)且二者均未出現在矩陣C的非主對角元素表達式中，從而可通過主對角占比目標函數

(13)

來表征阻抗矩陣的主對角元素取值集中程度，其中自變量coeff2=[a-1,a0,a1,b-1,b0,b1]T與coeff3=[c-1,c0,c1,d-1,d0,d1]T分別表示離散二階與三階偏導數的NAD網格模板系數所組成的向量.通過極小化(13)式可降低阻抗矩陣的條件數.若考慮理想情形，當函數dr為零時C將退化為對角矩陣從而條件數降為1，此時線性方程組的解立刻得出；另一方面，在期望函數dr減小的同時，還需保持相應數值格式離散的準確性，這里引入數值波數與理論波數的偏差來衡量差分格式壓制數值頻散的能力并以此作為約束條件(Zhan and Yao,2013;Jing et al.,2017).若對(5)式關于x作Fourier變換至波數域，則有波數表達式

(14)

其中(kx)num與kx分別表示x方向上數值波數與理論波數.將相位補償θx=kxh∈(0,π)代入(14)式中，則

(15)

根據(15)式可定義波數偏差函數

(16)

(17)

可利用MATLAB優化工具箱即可對(16)與(17)式進行有效求解，所得最優網格模板系數為

(18)

將(18)式代入(5)—(8)式，即可得改進NAD網格差分模板.事實上，上述改進思想不僅適用于四階NAD方法，對任意階標準NAD格式(例如：六階、八階)均可在其基礎上進行來提升頻率域正演模擬的效率.

對于線性方程組(10)的求解下文均采用不精確旋轉分塊下三角(IRBLT)預處理子加速GMRES方法進行，關于該預處理迭代方法的具體實施步驟以及理論分析結果可參見(Lang et al.,2020a).

1.3 改進NAD方法的正演效率評估

本節通過對均勻介質模型進行波場模擬來考查改進NAD(簡記為NAD*)方法的正演計算效率并與四階標準NAD(記作NAD4)方法進行對比.所考慮二維計算區域規模為5 km×5 km，網格參數201×201，背景速度c=3.5 km·s-1.

在下文所有數值試驗中，(1)式中的震源右端項均取作Ricker子波，其頻率表達式為(Lang et al.,2020b)

(19)

其中A為振幅，f0表示震源主頻.對該震源子波的頻譜分析結果顯示其大部分能量集中于0～2.5f0頻段內，因此可采用這一頻段的單頻波波場通過離散Fourier逆變換(IDFT)合成時間域波場來觀察數值頻散效應.為量化離散網格的疏密程度，需定義“網格頻率”(記作Gf)為最短步長所包含的網格點數(Liu,1997)

(20)

其中λmin表示最短波長，cmin為最小波速.若設置A=1,f0=20 Hz，則Gf=2.8為NAD4方法壓制數值頻散時網格疏密度的臨界值(Lang and Yang,2016).與此同時，還將數值格式的頻散誤差定義為不同波的傳播角度β下數值速度關于理論波速的最大相對值(Jing et al.,2017)

(21)

對于給定的數值離散格式，該頻散誤差隨網格頻率Gf增長而增加，這意味著離散網格越粗則頻散誤差越大，反之亦然.經過反復測試，當NAD4方法的網格頻率Gf(NAD4)=2.8以及NAD*方法的網格頻率Gf(NAD*)=2.92時，這兩種方法的頻散誤差相等

ε(NAD*)=ε(NAD4)=0.0258,(22)

此時在頻率取值分別為15,25,35,45 Hz下兩種方法所對應的阻抗矩陣條件數如表1所列.表中顯示各種情形下由NAD*方法離散所生成阻抗矩陣的條件數均小于NAD4方法；當f=8,35 Hz時由NAD*方法計算得到的單頻波波場快照如圖1所示；利用震源能量集中頻段(1～50 Hz)的單頻波波場快照分別合成對應于NAD*與NAD4方法無明顯數值頻散的時間域波場快照，如圖2所示.此時NAD4方法所需計算時間為3.19 min，相比之下，NAD*方法只需2.93 min，可節省約8.15%.

圖1 均勻介質中由NAD*方法計算得到的單頻波波場快照Fig.1 Mono-frequency wave-field snapshots by NAD* method in homogeneous medium

圖2 均勻介質中由兩種方法得到在t=0.5 s時刻的時間域波場快照Fig.2 Time-domain wave-field snapshots at t=0.5 s point by two methods in homogeneous medium

表1 不同頻率取值下兩種方法離散所得阻抗矩陣的條件數Table 1 Condition number of impedance matrix by two methods for various values of frequency

2 頻率域逆時偏移成像原理

經典的時間域逆時偏移互相關成像條件為(Claerbout,1971;Whitmore and Lines,1986)

(23)

=us(x,z,ω)·ur(x,z,ω),(24)

I(x,z)=F-1[us(x,z,ω)·ur(x,z,ω)]

(25)

為消除不同位置處震源強度對成像結果的影響，需要將(25)式進行歸一化(也稱為照明補償(Zhou et al.,2018)).根據Parseval等式，歸一化因子在時間域與頻率域保持能量不變性，即

(26)

則帶有照明補償的頻率域逆時偏移成像條件為

(27)

其中σ為取值很小的正數，引入該參數可避免(27)式分母為0.下文討論均采用帶有照明補償的互相關偏移成像條件.實際計算中，成像結果可通過選取(27)式的被積函數在所考慮頻段范圍0～ωmax中的不同頻率組分進行疊加而獲得.

在采用互相關逆時偏移條件對復雜地質結構進行成像時，由于首波、潛水波以及繞射波的存在，偏移過程將會產生嚴重的低頻噪聲，導致成像結果的分辨率受損(Qiu et al.,2020).為抑制低頻噪聲，通?？刹捎肔aplace濾波對波場關于傳播角度進行衰減(Zhang and Sun,2009).具體地，若將離散Laplace算子記為

(28)

則濾波后成像結果的矩陣形式為

IL(x,z)=L·I(x,z)+I(x,z)·L.

(29)

由(29)式可見，Laplace濾波操作簡單，所需計算量很少.

3 頻率域逆時偏移算例

為檢驗NAD*方法頻率域逆時偏移的數值效果，本節運用該方法分別對凹陷模型和Sigsbee2B模型進行偏移成像并將相應結果與NAD4方法和四階普通有限差分(簡記為OFD)方法進行對比.此外，我們還運用NAD*方法對Marmousi模型的標準數據進行計算.以下偏移計算中，接收器位于計算區域頂部與PML分界面處并且逐個節點密集放置，震源在接收器下一層并根據模型的復雜程度進行選取，PML區域厚度為δ=10 h，對于各種模型其他計算參數的取值情況見表2.

表2 對各種介質模型進行頻率域逆時偏移的基本參數表Table 2 Parameter table of frequency-domain reverse time migration for various media models

3.1 凹陷模型

首先考慮凹陷模型(如圖3所示)，其上半區域背景速度取值為2.5 km·s-1,下半區域為3.0 km·s-1.由于模型相對簡單，在偏移成像中每隔25個節點設有1個震源并且沿水平方向中心對稱放置，共計7個震源.在表2所列的基本參數設置下，網格頻率Gf=2.8.

圖3 凹陷模型Fig.3 Sunken model

本小節分別采用NAD*、NAD4與OFD方法進行偏移計算，所得成像結果如圖4所示.從圖中可見兩種NAD方法的圖像界面清晰，分辨率較高；而OFD方法的成像結果中在界面上方存在明顯的數值頻散并且區域底部出現少許假象，這是由于此時離散網格對于OFD方法而言過于粗糙而導致波場計算的精度偏低，正演模擬的數值誤差體現在最終偏移成像結果中.

根據以上分析，可通過加密離散網格消除圖4中OFD方法成像結果的數值頻散.經過測試，當步長減少至27.0 m時頻散開始消失，此時的網格規模為259×259,相應成像結果如圖5所示.我們統計三種方法在得到準確偏移結果時所需的最少計算時間(如表3所列)，發現NAD*方法用時最少，相比于NAD4與OFD方法，大致分別節省10.05%和19.18%.

圖4 凹陷模型偏移成像結果(a)改進NAD方法;(b)四階NAD方法;(c)四階OFD方法.Fig.4 Migration imaging results of sunken model(a)Improved NAD method;(b)4th-order NAD method;(c)4th-order OFD method.

圖5 細網格(h=27.0 m)下由OFD方法計算凹陷模型的偏移成像結果Fig.5 Migration imaging results of sunken model on fine grid (h=27.0 m)by OFD method

表3 不同方法計算準確偏移結果所需最少計算時間(單位：min)Table 3 Least time of computing accurate migration imaging for various methods (unit:min)

3.2 Sigsbee2B模型

Sigsbee2B模型是偏移方法基準測試的國際通用模型(Paffenholz et al.,2002)，接下來對其進行成像計算.首先將原始模型進行重采樣(現規模為801×301)并將速度值換算至國際標準單位(km·s-1).如圖6所示，計算區域中部出現一個高速異常的鹽丘體，對鹽丘底部層進行準確成像是該模型成像的難點.為此，本次試驗放置更多震源，間隔10個節點，共計79個.

圖6 Sigsbee2B模型速度結構Fig.6 Velocity of Sigsbee2B model

圖7中分別顯示了由NAD*、NAD4以及OFD方法計算的偏移成像結果，可以看到由NAD類方法所得到的結果相對清晰，并對鹽丘底部區域進行了準確刻畫，而OFD方法的圖像中在鹽丘頂部與底部(例如：箭頭所指處)均出現虛假的層狀結構，此為數值頻散，需加密離散網格來消除這些假象.經過試驗，當網格步長降至12.2 m時頻散開始消失，相應偏移結果可見圖8.此時三種方法得到準確結果所需最少計算時間如表3所列，同樣NAD*方法用時最少，相比于NAD4與OFD方法可分別加速14.63%和25.86%.

圖7 Sigsbee2B模型偏移成像結果(a)NAD*方法;(b)NAD4方法;(c)OFD方法.Fig.7 Migration imaging results of Sigsbee2B model(a)NAD* method;(b)NAD4 method;(c)OFD method.

圖8 細網格(h=12.2 m)下由四階OFD方法計算的Sigsbee2B模型偏移成像結果Fig.8 Migration imaging results of Sigsbee2B model on fine grids (h=12.2 m)computed by 4th-order OFD method

3.3 Marmousi模型

最后運用NAD*方法對Marmousi模型(Versteeg,1994)(如圖9所示)進行偏移成像，本小節采用標準數據進行計算(某單炮記錄如圖10所示)，其他參數設置見表2.由NAD*方法計算得到的偏移成像結果如圖11所示，盡管區域頂部受噪聲干擾明顯，但模型的分層界面基本顯現，并且未見數值頻散，這也進一步檢驗了文中所述頻率域偏移成像方法的正確性.

圖9 Marmousi速度模型Fig.9 Marmousi velocity model

圖10 標準數據的某單炮記錄Fig.10 Single-shot record of standard data

圖11 由NAD*方法計算得到Marmousi模型偏移成像結果Fig.11 Migration imaging results of Marmousi model computed by NAD* method

4 結論

逆時偏移成像的關鍵在于正演過程中如何快速有效地求解波動方程，作為一種高精度數值離散方法，NAD方法尤其適用于頻率域正演模擬與偏移成像計算.本文通過優化NAD格式的網格差分模板系數，構造改進NAD方法.相比于標準NAD差分格式，改進NAD在數值離散精度不受損失的前提下降低了離散后所得阻抗矩陣的條件數，從而節省相應線性方程組求解所需的計算代價并使得迭代求解過程更加穩定.因此，采用改進NAD方法進行正演模擬可提高頻率域計算效率,有助于增進大尺度區域頻率域波場模擬與成像的實用性.隨后，本文將改進NAD方法與頻率域逆時偏移原理相結合，分別對凹陷模型、Sigsbee2B模型以及Marmousi模型實施偏移成像，并與標準四階NAD方法和四階OFD方法進行對比.數值結果顯示了改進NAD方法在壓制數值頻散與提高偏移成像效率方面的顯著優勢；特別是對于Sigsbee2B模型的偏移試驗，當得到準確成像結果時改進NAD方法相比于四階NAD與OFD方法可分別節省約14.63%和25.86%的計算時間；Marmousi模型的可靠偏移成像結果進一步表明改進NAD方法對于復雜介質模型成像的有效性.

附錄A 阻抗矩陣的稀疏結構推導

當采用NAD網格差分模板(5)—(8)式離散頻率域聲波方程時，(4a)式可離散為

(A1)

(4b)式離散為

(A2)

(4c)式離散為

(A3)

若將二維區域Ω中的網格節點按照逐行次序排列(每行從左至右)，則有網格節點序號對應關系k=nx·(j-1)+i，并且同一節點處的波場及其梯度項如下放置

(A4)

可形成阻抗矩陣的稀疏分塊結構(11)式.

附錄B 阻抗矩陣的非零子塊具體元素表達式

(B1)

(B2)

(B3)

(B4)

(B5)

(B6)

(B7)

(B8)

(B9)

其中j為分塊行號，跟隨區域內節點的位置而變化，上述表達式中變量dx,dz,tx,tz的取值依賴于該指標.