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

1 北京信息科技大學(xué)理學(xué)院，北京 100192 2 應(yīng)急管理部國家自然災(zāi)害防治研究院，北京 100085

0 引言

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

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

在逆時(shí)偏移成像中，絕大部分計(jì)算操作集中于模擬地下介質(zhì)中地震波的傳播情況(正演過程).因此，偏移效率以及成像結(jié)果的質(zhì)量取決正演算法的有效性 (嚴(yán)紅勇和劉洋,2013;Qiu et al.,2020).目前地球物理領(lǐng)域常用的正演數(shù)值方法包括有限差分方法(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)等.區(qū)別于時(shí)間域波動(dòng)方程的數(shù)值離散，選取頻率域數(shù)值方法須在保證離散格式準(zhǔn)確性的同時(shí)，還要兼顧到求解離散后所得大型線性代數(shù)方程組的計(jì)算代價(jià)可控(Pratt and Worthington,1990)，這是大尺度區(qū)域模型頻率域計(jì)算所面臨的主要挑戰(zhàn)之一(Ernst and Gander,2012).求解該線性方程組的計(jì)算代價(jià)主要依賴于其系數(shù)矩陣(也稱為“阻抗矩陣”(Pratt et al.,1998))的稀疏性與條件數(shù)，這通常由離散算子的結(jié)構(gòu)與性質(zhì)所決定，從而頻率域波動(dòng)方程的數(shù)值離散應(yīng)當(dāng)盡可能避免復(fù)雜格式以及線性方程組的病態(tài)程度過高(Lang and Yang,2016).有限差分方法具有原理直觀、格式構(gòu)造簡(jiǎn)單、易于實(shí)現(xiàn)、方便并行等優(yōu)點(diǎn)，有利于頻率域波動(dòng)方程的數(shù)值離散.這類方法已廣泛應(yīng)用于頻率域地震波場(chǎng)模擬、逆時(shí)偏移成像以及全波形反演等方面(Yang et al.,2006;Xu et al.,2010;Lang et al.,2020b).然而普通差分格式在粗網(wǎng)格下較難對(duì)復(fù)雜地質(zhì)結(jié)構(gòu)(例如：強(qiáng)間斷面、孔縫等)中地震波的傳播情況進(jìn)行準(zhǔn)確模擬(Liu et al.,2017b)，容易引起嚴(yán)重的數(shù)值頻散假象(Moczo et al.,2000)，進(jìn)而導(dǎo)致偏移成像結(jié)果的分辨率受損，此時(shí)若加密離散網(wǎng)格，問題規(guī)模也隨之增加，因此計(jì)算效率偏低.

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

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

1 正演算法

1.1 NAD離散過程

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

(1)

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

(2)

為消除人工截?cái)嘤?jì)算區(qū)域所產(chǎn)生的虛假反射波，需對(duì)Ω的邊界進(jìn)行處理.本文采用“完美匹配層(Perfectly Matched Layer)”吸收邊界條件(Berenger,1994).應(yīng)該指出的是，在頻率域施加吸收邊界條件來刻畫波的衰減特性具有天然的便捷性，可避免因反卷積運(yùn)算所產(chǎn)生的誤差(Lang and Yang,2016).首先引入復(fù)坐標(biāo)(以水平方向?yàn)槔?

(3)

(4a)

(4b)

(4c)

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

(5)

z方向討論類似，則(4)式中的高階混合偏導(dǎo)數(shù)可借助(5)式以及方向?qū)?shù)的性質(zhì)進(jìn)行推導(dǎo)(Tong et al.,2011)

(6)

(7)

以及

(8)

需要指出的是，當(dāng)網(wǎng)格模板系數(shù)取作

(9)

此時(shí)對(duì)應(yīng)于四階標(biāo)準(zhǔn)NAD方法(Yang et al.,2006).

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

Cu=b,(10)

矩陣C的稀疏分塊結(jié)構(gòu)如下：

(11)

1.2 改進(jìn)NAD方法

(12)

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

(13)

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

(14)

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

(15)

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

(16)

(17)

可利用MATLAB優(yōu)化工具箱即可對(duì)(16)與(17)式進(jìn)行有效求解，所得最優(yōu)網(wǎng)格模板系數(shù)為

(18)

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

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

1.3 改進(jìn)NAD方法的正演效率評(píng)估

本節(jié)通過對(duì)均勻介質(zhì)模型進(jìn)行波場(chǎng)模擬來考查改進(jìn)NAD(簡(jiǎn)記為NAD*)方法的正演計(jì)算效率并與四階標(biāo)準(zhǔn)NAD(記作NAD4)方法進(jìn)行對(duì)比.所考慮二維計(jì)算區(qū)域規(guī)模為5 km×5 km，網(wǎng)格參數(shù)201×201，背景速度c=3.5 km·s-1.

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

(19)

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

(20)

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

(21)

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

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

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

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

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

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

2 頻率域逆時(shí)偏移成像原理

經(jīng)典的時(shí)間域逆時(shí)偏移互相關(guān)成像條件為(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)

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

(26)

則帶有照明補(bǔ)償?shù)念l率域逆時(shí)偏移成像條件為

(27)

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

在采用互相關(guān)逆時(shí)偏移條件對(duì)復(fù)雜地質(zhì)結(jié)構(gòu)進(jìn)行成像時(shí)，由于首波、潛水波以及繞射波的存在，偏移過程將會(huì)產(chǎn)生嚴(yán)重的低頻噪聲，導(dǎo)致成像結(jié)果的分辨率受損(Qiu et al.,2020).為抑制低頻噪聲，通常可采用Laplace濾波對(duì)波場(chǎng)關(guān)于傳播角度進(jìn)行衰減(Zhang and Sun,2009).具體地，若將離散Laplace算子記為

(28)

則濾波后成像結(jié)果的矩陣形式為

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

(29)

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

3 頻率域逆時(shí)偏移算例

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

表2 對(duì)各種介質(zhì)模型進(jìn)行頻率域逆時(shí)偏移的基本參數(shù)表Table 2 Parameter table of frequency-domain reverse time migration for various media models

3.1 凹陷模型

首先考慮凹陷模型(如圖3所示)，其上半?yún)^(qū)域背景速度取值為2.5 km·s-1,下半?yún)^(qū)域?yàn)?.0 km·s-1.由于模型相對(duì)簡(jiǎn)單，在偏移成像中每隔25個(gè)節(jié)點(diǎn)設(shè)有1個(gè)震源并且沿水平方向中心對(duì)稱放置，共計(jì)7個(gè)震源.在表2所列的基本參數(shù)設(shè)置下，網(wǎng)格頻率Gf=2.8.

圖3 凹陷模型Fig.3 Sunken model

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

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

圖4 凹陷模型偏移成像結(jié)果(a)改進(jìn)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 細(xì)網(wǎng)格(h=27.0 m)下由OFD方法計(jì)算凹陷模型的偏移成像結(jié)果Fig.5 Migration imaging results of sunken model on fine grid (h=27.0 m)by OFD method

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

3.2 Sigsbee2B模型

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

圖6 Sigsbee2B模型速度結(jié)構(gòu)Fig.6 Velocity of Sigsbee2B model

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

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

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

3.3 Marmousi模型

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

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

圖10 標(biāo)準(zhǔn)數(shù)據(jù)的某單炮記錄Fig.10 Single-shot record of standard data

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

4 結(jié)論

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

附錄A 阻抗矩陣的稀疏結(jié)構(gòu)推導(dǎo)

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

(A1)

(4b)式離散為

(A2)

(4c)式離散為

(A3)

若將二維區(qū)域Ω中的網(wǎng)格節(jié)點(diǎn)按照逐行次序排列(每行從左至右)，則有網(wǎng)格節(jié)點(diǎn)序號(hào)對(duì)應(yīng)關(guān)系k=nx·(j-1)+i，并且同一節(jié)點(diǎn)處的波場(chǎng)及其梯度項(xiàng)如下放置

(A4)

可形成阻抗矩陣的稀疏分塊結(jié)構(gòu)(11)式.

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

(B1)

(B2)

(B3)

(B4)

(B5)

(B6)

(B7)

(B8)

(B9)

其中j為分塊行號(hào)，跟隨區(qū)域內(nèi)節(jié)點(diǎn)的位置而變化，上述表達(dá)式中變量dx,dz,tx,tz的取值依賴于該指標(biāo).