徐云貴， 廖建平， 周林， 劉和秀， 張青， 謝敬濤, 王立歆

1 西南石油大學(xué)地球科學(xué)與技術(shù)學(xué)院, 成都 610500 2 西南石油大學(xué)油氣藏地質(zhì)及開發(fā)工程國(guó)家重點(diǎn)實(shí)驗(yàn)室, 成都 610500 3 湖南科技大學(xué)地球科學(xué)與空間信息工程學(xué)院, 湖南湘潭 411201 4 中石化地球物理重點(diǎn)實(shí)驗(yàn)室, 中石化南京物探研究院, 南京 211103 5 重慶科技學(xué)院非常規(guī)油氣開發(fā)研究院， 重慶 401331

0 引言

地下巖石中的天然裂縫可形成流體的有效通道或儲(chǔ)集空間，而含裂縫巖石的物理性質(zhì)因流體的存在和流體的流動(dòng)發(fā)生變化.為了解含流體裂縫巖石的物理性質(zhì)以及其中的流體流動(dòng)特性，首先必須研究裂縫本身屬性(Li, 1999).在地震波動(dòng)力學(xué)中，巖石裂縫的存在影響地震波傳播方式，一種有效的研究裂縫方法是通過(guò)地震正演模擬來(lái)研究裂縫系統(tǒng)的地震響應(yīng)，并與現(xiàn)場(chǎng)采集的實(shí)際地震資料進(jìn)行比較，以確定裂縫空間發(fā)育規(guī)律(Tsvankin et al., 2010).裂縫介質(zhì)地震正演模擬在裂縫研究中具有兩個(gè)重要的作用：一是揭示各種裂縫介質(zhì)對(duì)應(yīng)的地震響應(yīng)規(guī)律，二是驗(yàn)證現(xiàn)有裂縫反演方法的有效性(Liu et al., 2000; Mavko et al., 2009).但到目前為止，很多裂縫反演和正演模擬研究都局限于二維情況(Nihei et al., 2002; Vlastos et al., 2003; Vlastos, 2005; Rao and Wang, 2009).三維裂縫介質(zhì)中地震正演模擬因存在計(jì)算量大等問(wèn)題限制了其應(yīng)用規(guī)模.然而，裂縫介質(zhì)的三維地震響應(yīng)研究意義重大，三維地震數(shù)據(jù)更能體現(xiàn)地震波在真實(shí)意義的三維各向異性介質(zhì)中的傳播機(jī)理，從而能依據(jù)各向異性理論和方法來(lái)表征裂縫空間發(fā)育規(guī)律.一個(gè)典型的例子是，在地震勘探中基于三維實(shí)際疊前地震縱波數(shù)據(jù)，通過(guò)分析地震波振幅或速度的方位各向異性來(lái)反演高角度裂縫密度，提高對(duì)高角度裂縫空間發(fā)育規(guī)律的認(rèn)識(shí).

裂縫介質(zhì)的地震正演模擬大多采用有限差分方法(Finite Difference，簡(jiǎn)稱FD).Willis等(2006)和Grandi-Karam (2008)進(jìn)行了三維離散裂縫模型(Discrete Fracture Model，簡(jiǎn)稱DFM)的正演模擬研究.他們討論了DFM模擬數(shù)據(jù)如何利用散射波確定裂縫方向和裂縫間距，其研究?jī)?nèi)容使用現(xiàn)有有限差分工具實(shí)現(xiàn)三維FD正演模擬，并未討論模擬方法的具體實(shí)現(xiàn)過(guò)程.本文將討論裂縫介質(zhì)中FD正演模擬的更多細(xì)節(jié)，譬如正演模擬實(shí)現(xiàn)、優(yōu)化、模擬方法的選擇等.

裂縫介質(zhì)的FD正演模擬常采用三種交錯(cuò)網(wǎng)格方法.第一種是Virieux(1984,1986)提出的標(biāo)準(zhǔn)交錯(cuò)網(wǎng)格方法(SSG)(Levander, 1988; Dong and McMechan, 1995; Graves, 1996; Liu, 2014; Ekanem and Xu, 2018).在SSG方法中，將波場(chǎng)分量離散到不同的數(shù)值網(wǎng)格中，計(jì)算所需網(wǎng)格位置的波場(chǎng)空間導(dǎo)數(shù)，迭代求解波動(dòng)方程，即得到模擬結(jié)果.SSG方法簡(jiǎn)潔且容易編寫程序?qū)崿F(xiàn)，能有效計(jì)算正交各向異性介質(zhì)模型.第二種是Saenger等(2000)提出了旋轉(zhuǎn)交錯(cuò)網(wǎng)格方法(RSG)，該方法通過(guò)組合網(wǎng)格對(duì)角線上的波場(chǎng)分量求解空間導(dǎo)數(shù)，有利于處理任意復(fù)雜各向異性介質(zhì)模型，但相對(duì)SSG方法，RSG方法計(jì)算資源耗費(fèi)更高(Saenger和Bohlen，2004).第三種是Lisitsa(2007)將Lebedev格式(LS，Lebedev, 1964)引入地震正演模擬，波場(chǎng)分量呈菱形分布在FD網(wǎng)格中，可沿坐標(biāo)軸求解空間導(dǎo)數(shù)，Lebedev格式也是一種交錯(cuò)網(wǎng)格方法(Lisitsa et al., 2009；Lisitsa and Vishnevsky，2011；Lisitsa et al., 2012a；Lisitsa et al., 2012b；Vishnevsky et al., 2014；Quintanilla and Leckey，2018；Koene et al., 2021a; Koene et al., 2021b).Lisitsa 和Vishnevskiy(2010)定量比較了任意各向異性介質(zhì)中的RSG和LS方法.然而，三維地震正演模擬需要厘清不同的FD交錯(cuò)網(wǎng)格方法在模擬各種各向異性介質(zhì)(如TTI和正交介質(zhì))中的適用性，因?yàn)椴煌母飨虍愋越橘|(zhì)剛度矩陣形式不同，矩陣形式?jīng)Q定交錯(cuò)網(wǎng)格方法的選取和計(jì)算資源的耗費(fèi)大小.這一點(diǎn)非常重要，特別是在三維模擬的情況下，大型三維各向異性模型的正演模擬通常需要幾天、幾周甚至長(zhǎng)達(dá)數(shù)月.這種對(duì)比研究有利于方法的選擇和計(jì)算過(guò)程的優(yōu)化，在基于各向異性正演模擬的地震資料處理和解釋時(shí)尤為重要(Liao et al., 2022; Wang et al., 2022)，如逆時(shí)偏移 (RTM, Reverse Time Migration, Shi et al., 2019) 和全波形反演(FWI, Full Waveform Inversion, Guo et al., 2019).

國(guó)內(nèi)的學(xué)者們對(duì)Lebedev方法進(jìn)行了深入的研究，并成功應(yīng)用于地震波數(shù)值模擬中.李娜等(2014a)實(shí)現(xiàn)了Lebedev網(wǎng)格與標(biāo)準(zhǔn)交錯(cuò)網(wǎng)格耦合機(jī)制下的復(fù)雜各向異性正演模擬的研究.李娜等(2014b)進(jìn)行了Lebedev網(wǎng)格改進(jìn)差分系數(shù)TTI介質(zhì)正演模擬方法研究.黃建平等(2016)開展了基于M-PML邊界的Lebedev網(wǎng)格起伏地表正演模擬方法及穩(wěn)定性分析的研究.楊宇等(2016)實(shí)現(xiàn)了Lebedev網(wǎng)格黏彈性介質(zhì)起伏地表正演模擬.黃金強(qiáng)等(2017)實(shí)現(xiàn)了TTI介質(zhì)Lebedev網(wǎng)格高階有限差分正演模擬及波型分離.劉東洋等(2018)基于Lebedev網(wǎng)格的TTI介質(zhì)二維三分量正演模擬，引入多軸完全匹配層(M-PML)吸收邊界條件后，在不影響模擬效果的情況下邊界反射現(xiàn)象被有效地壓制.胡書華等(2018)進(jìn)行了TTI介質(zhì)穩(wěn)定的純qP波波場(chǎng)模擬方法研究，他們?cè)贚ebedev交錯(cuò)網(wǎng)格框架下推導(dǎo)了高階有限差分方程，并給出2次計(jì)算波場(chǎng)梯度的數(shù)值算法，以保證波場(chǎng)模擬精度.

三種交錯(cuò)網(wǎng)格求解波場(chǎng)空間導(dǎo)數(shù)的不同求解方式?jīng)Q定了它們對(duì)特定類型各向異性介質(zhì)的模擬的能力.SSG很容易模擬簡(jiǎn)單各向異性介質(zhì)(VTI，HTI，正交各向異性)，而其他兩種方法可以處理任意各向異性介質(zhì)，但計(jì)算成本較高(Xu and Wang, 2017).本文比較了三種主要有限差分交錯(cuò)網(wǎng)格方法，厘清了每種方法的優(yōu)勢(shì)和局限，提出了一種優(yōu)化的正演模擬計(jì)算流程，并改進(jìn)了LS方法的數(shù)值實(shí)現(xiàn).通過(guò)二維與三維合成實(shí)例驗(yàn)證了該方法的正確性與穩(wěn)健性.

1 理論

1.1 三維彈性各向異性波動(dòng)方程的速度-應(yīng)力公式

地震各向異性模擬基于各向異性彈性波動(dòng)方程，波動(dòng)方程可用一階速度-應(yīng)力關(guān)系描述，包括胡克定律(方程(1))和牛頓第二定律(方程(2))：

(1)

(2)

1.2 各向異性介質(zhì)中的地震正演模擬的三種交錯(cuò)網(wǎng)格

有限差分求解上述波動(dòng)方程主要有三種不同的FD交錯(cuò)網(wǎng)格方法(SSG，RSG和LS)，使用交錯(cuò)網(wǎng)格是為了實(shí)現(xiàn)基于中心差分求解空間導(dǎo)數(shù)，三種交錯(cuò)網(wǎng)格方法的共同特征是：(1)正演模擬中至少有兩個(gè)FD交錯(cuò)網(wǎng)格；(2)一個(gè)網(wǎng)格位于另一個(gè)網(wǎng)格的中點(diǎn)處，確?？梢曰谥行牟罘智蠼饪臻g導(dǎo)數(shù)；(3)有限差分求解空間導(dǎo)數(shù)時(shí)，只在導(dǎo)數(shù)位置處使用相鄰的幾個(gè)點(diǎn)的波場(chǎng)分量或模型變量，計(jì)算效率高.

以下詳細(xì)討論三種方法在模擬各種各向異性介質(zhì)模型的優(yōu)勢(shì)和局限.

1.2.1 SSG方法

Virieux(1984, 1986)做了較早的SSG研究工作，用于求解二維各向同性介質(zhì)中的地震波傳播問(wèn)題，SSG中有4種不同的相互交錯(cuò)的網(wǎng)格參數(shù)(如圖1a).二維SSG方法由Dong和McMechan(1995)擴(kuò)展到三維黏彈性各向異性正演模擬，涉及7種不同的網(wǎng)格(圖1b)，這里和后面提到的分量或參數(shù)見方程(1)和(2).不同的波場(chǎng)分量或模型參數(shù)都分布在交錯(cuò)的網(wǎng)格上，以便使用中心差分來(lái)求解速度或應(yīng)力的空間導(dǎo)數(shù)，三維模型的7種網(wǎng)格交錯(cuò)分布的特定位置是SSG方法的特點(diǎn).SSG方法容易實(shí)現(xiàn)、計(jì)算速度快、計(jì)算成本低，因此大多數(shù)正演模擬都是基于SSG方法.SSG方法在處理垂直對(duì)稱軸的橫向各向同性(VTI)介質(zhì)、水平對(duì)稱軸的橫向各向同性(HTI)介質(zhì)和正交介質(zhì)時(shí)都表現(xiàn)良好，因?yàn)檫@些介質(zhì)的剛度矩陣c的非零常數(shù)分量可以放在圖1b的1位置，便于求解中心差分導(dǎo)數(shù).三維SSG的四個(gè)特點(diǎn)是：(1)SSG方法包括7種交錯(cuò)的網(wǎng)格(見圖1b中的矩形和圓形編號(hào)為1—7的網(wǎng)格)；(2)計(jì)算速度快，內(nèi)存成本??；(3)能夠很好地處理各向同性、VTI、HTI和正交介質(zhì)；(4)比正交各向異性更復(fù)雜的介質(zhì)(TTI、單斜和三斜介質(zhì))的正演模擬，需要對(duì)波場(chǎng)進(jìn)行插值(牟永光, 裴正林, 2005).例如，TTI、單斜和三斜介質(zhì)中的正演模擬，速度分量需要插值，從而導(dǎo)致產(chǎn)生誤差.因此，盡量采用其他的交錯(cuò)網(wǎng)格方法去避免插值.

圖1 三種不同的交錯(cuò)網(wǎng)格(a) 二維標(biāo)準(zhǔn)交錯(cuò)網(wǎng)格； (b) 三維標(biāo)準(zhǔn)交錯(cuò)網(wǎng)格； (c) 二維旋轉(zhuǎn)交錯(cuò)網(wǎng)格； (d) 三維旋轉(zhuǎn)交錯(cuò)網(wǎng)格；(e) 二維菱形交錯(cuò)網(wǎng)格； (f) 三維菱形交錯(cuò)網(wǎng)格.Fig.1 Three staggered-grids The SSG in 2D case (a) and 3D case (b), the RSG in 2D case (c) and 3D case (d), and the LS in 2D case (e) and 3D case (f).

1.2.2 RSG方法

為了能夠適應(yīng)更復(fù)雜的各向異性介質(zhì)的模擬，Saenger等(2000)引入了一個(gè)新的有限差分交錯(cuò)網(wǎng)格實(shí)現(xiàn)方法，即RSG，在二維或三維時(shí)，只有2種交錯(cuò)的網(wǎng)格(圖1c和圖1d).網(wǎng)格越少，方法實(shí)現(xiàn)越簡(jiǎn)單，所有的速度分量在一個(gè)交錯(cuò)網(wǎng)格點(diǎn)上，所有的應(yīng)力分量在另一個(gè)網(wǎng)格上.不同的是，在求解應(yīng)力分量的空間導(dǎo)數(shù)時(shí)(例如在圖1d中圓形位置的應(yīng)力空間導(dǎo)數(shù))，需要輸入該導(dǎo)數(shù)位置周圍的4個(gè)對(duì)角線網(wǎng)格點(diǎn)上的8個(gè)應(yīng)力分量(8個(gè)矩形)的線性組合；在求解速度分量的空間導(dǎo)數(shù)時(shí)，則需要這些網(wǎng)格點(diǎn)上的速度分量組合.此特點(diǎn)可以實(shí)現(xiàn)任意各向異性介質(zhì)中地震正演模擬(包括TTI、單斜和三斜介質(zhì))，而不需要對(duì)波場(chǎng)分量進(jìn)行插值.該方法的缺點(diǎn)是求解空間導(dǎo)數(shù)時(shí)，需要對(duì)角線上的4個(gè)分量參與計(jì)算，導(dǎo)致計(jì)算量比SSG大很多.此外，由于使用對(duì)角線分量計(jì)算導(dǎo)數(shù)，對(duì)角線的長(zhǎng)度作為求導(dǎo)的單位長(zhǎng)度，為了保持與SSG方法相同的精度，RSG模型的單位體積內(nèi)需要更多的網(wǎng)格點(diǎn)和計(jì)算量.RSG的四個(gè)特點(diǎn)是：(1)模型中只有矩形和圓形兩種網(wǎng)格，簡(jiǎn)化了模型；(2)比SSG方法需要更多的計(jì)算和內(nèi)存需求；(3)可以很好地處理任意各向異性介質(zhì)；(4)需要線性組合對(duì)角線的波場(chǎng)分量求解空間導(dǎo)數(shù).

1.2.3 LS方法

方程(1)和(2)的重要特點(diǎn)是：方程右邊波場(chǎng)分量的中心差分的空間導(dǎo)數(shù)乘以密度或剛度常量，如果落在方程左邊分量的位置，可以利用交錯(cuò)網(wǎng)格來(lái)解方程(1)和(2)，否則需要插值.在三維LS中，取一個(gè)分量，將另外一個(gè)分量布置在該分量的上下、前后和左右的6個(gè)位置(組成菱形形狀)，則可以完全滿足交錯(cuò)網(wǎng)格的要求，而且對(duì)于任意介質(zhì)的剛度矩陣都適用.因此，LS方法適用于任意各向異性介質(zhì)的波場(chǎng)正演模擬.

在LS方法中(圖1e和圖1f)，相同的場(chǎng)分量組成菱形網(wǎng)格(圖1e和圖1f的虛線的對(duì)象)，這些網(wǎng)格的特點(diǎn)是：可以在不對(duì)波場(chǎng)分量插值的前提下，利用中心差分方法在期望的網(wǎng)格位置處求取所有的導(dǎo)數(shù)，因而可以模擬任意各向異性介質(zhì)，它不同于SSG中的元素需要插值，也不同于RSG中的導(dǎo)數(shù)組合.例如，在圖1f中虛的菱形網(wǎng)格，三個(gè)速度空間導(dǎo)數(shù)沿三個(gè)軸，在三維菱形的6個(gè)頂點(diǎn)處，可以用6個(gè)速度分量在菱形的中心處求解三個(gè)速度導(dǎo)數(shù).應(yīng)力空間導(dǎo)數(shù)求取與此相同.因而，LS在數(shù)值實(shí)現(xiàn)上更加直接而且更容易理解.

從圖1e和1f可以看出菱形網(wǎng)格中，每個(gè)單元網(wǎng)格有更多的分量網(wǎng)格節(jié)點(diǎn)，因此需要更多的計(jì)算內(nèi)存，同時(shí)更多的分量網(wǎng)格會(huì)導(dǎo)致更多的計(jì)算時(shí)間.由于計(jì)算機(jī)適合處理使用二維/三維規(guī)則矩形網(wǎng)格，這種不尋常的菱形特性會(huì)給正演模擬數(shù)值實(shí)現(xiàn)帶來(lái)困難.LS方法的優(yōu)點(diǎn)是網(wǎng)格少，只有兩個(gè)網(wǎng)格，處理簡(jiǎn)單，在任意各向異性介質(zhì)中正演模擬時(shí)，不需要SSG方法中的插值，也不需要RSG方法中的波場(chǎng)分量的線性組合.LS的四個(gè)特點(diǎn)是：(1)只涉及兩個(gè)網(wǎng)格(矩形和圓形兩種網(wǎng)格)；(2)與SSG相比，計(jì)算量更多，存儲(chǔ)成本更高，但均比RSG低；(3)能夠更好地處理任意各向異性介質(zhì)；(4)網(wǎng)格是菱形.

在任意各向異性介質(zhì)中，方程(1)和方程(2)的第一個(gè)表達(dá)式的離散LS形式為

參數(shù)與方程(1)、(2)相同.

當(dāng)i+j+k為奇數(shù)時(shí)，

(3)

當(dāng)i+j+k為偶數(shù)時(shí)，

(4)

值得注意的是，方程(3)和方程(4)的假設(shè)是奇數(shù)與偶數(shù)條件，這里符號(hào)參數(shù)與方程(1)和(2)相同.

總之，對(duì)于三種交錯(cuò)網(wǎng)格，在定義應(yīng)力分量的位置上用中心差分求解速度分量的導(dǎo)數(shù)，反之亦然.SSG網(wǎng)格數(shù)比RSG和LS多.RSG在網(wǎng)格的對(duì)角線求解波場(chǎng)的導(dǎo)數(shù).LS沿坐標(biāo)軸求解波場(chǎng)導(dǎo)數(shù)，但導(dǎo)致非矩形網(wǎng)格，即菱形網(wǎng)格.

1.2.4 三種方法的定量比較

以上定性討論了三種交錯(cuò)網(wǎng)格方法的適用性問(wèn)題，為定量確定各種方法在模擬時(shí)的真實(shí)計(jì)算量和計(jì)算機(jī)內(nèi)存消耗，以下按圖2所示的三維單元網(wǎng)格為基礎(chǔ)展開討論.

圖2 三種方法中的單元格(a) SSG； (b) RSG； (c) LS.由于交錯(cuò)網(wǎng)格的原因，這里單元格的長(zhǎng)度是網(wǎng)格間距的一半.Fig.2 Unit cells in the three schemes: SSG(a), RSG(b) and LS(c) Note the length of the unit cells here are half the cell spacing due to the staggered-grid.

為了有效對(duì)比三種方法，以同樣幾何尺度的正交各向異性模型為例，并且以同樣的精度為基準(zhǔn).通過(guò)比較最終的內(nèi)存消耗和計(jì)算量，就可以確定每種方法的計(jì)算適應(yīng)性特征.以下進(jìn)行詳細(xì)討論與比較(參考表1)：

表1 相同大小模型、相同精度條件下三種方法的計(jì)算機(jī)內(nèi)存消耗和計(jì)算量比較Table 1 Memory cost and computation for the three scheme when the model with the same size and the same accuracy are considered

(1)首先討論三維SSG方法的計(jì)算要求.假設(shè)模型大小為n×n×n(n為每個(gè)邊的樣點(diǎn)數(shù)，單元網(wǎng)格邊長(zhǎng)相等)，參考圖2a，共有7種可能的交錯(cuò)網(wǎng)格，對(duì)應(yīng)11個(gè)波場(chǎng)和模型變量，內(nèi)存消耗量為n×n×n×11.在模擬過(guò)程中，計(jì)算量主要取決于空間導(dǎo)數(shù)的計(jì)算，因此可以利用求解導(dǎo)數(shù)的次數(shù)確定計(jì)算量的大小，依據(jù)方程(1)和(2)，SSG方法中一個(gè)單元網(wǎng)格(如圖2a)需要求取18次空間導(dǎo)數(shù)，即總計(jì)算量為n×n×n×18.

(3)若要達(dá)到相同的精度，LS方法在每個(gè)單元網(wǎng)格上所需的變量網(wǎng)格數(shù)是4(如圖2c，網(wǎng)格點(diǎn)為4)，故三維模型共需要n×n×n×4個(gè)網(wǎng)格點(diǎn)，涉及11個(gè)波場(chǎng)變量和模型變量，內(nèi)存消耗量為n×n×n×4×11，共須求取18次空間導(dǎo)數(shù)，總計(jì)算量是n×n×n×18×4.

綜上所述，對(duì)于SSG、RSG和LS而言，內(nèi)存消耗比例為1∶5.2∶4，計(jì)算量大小比例為1∶20.8∶4，由此可見，SSG計(jì)算效率最高，LS其次，RSG最低.SSG只能有效模擬簡(jiǎn)單各向異性介質(zhì)(各向同性、VTI、HTI和正交各向異性)，RSG和LS則可以模擬任意復(fù)雜的各向異性介質(zhì)(TTI，單斜和三斜介質(zhì)).對(duì)于復(fù)雜的各向異性介質(zhì)，LS消耗的內(nèi)存和計(jì)算量相對(duì)RSG較小，因此LS方法更適合模擬復(fù)雜各向異性介質(zhì)模型.

2 改進(jìn)后的LS實(shí)施方法

Lisitsa和Vishnevskiy(2010)給出了任意各向異性介質(zhì)中使用LS方法進(jìn)行正演模擬的詳細(xì)描述.然而，他們并沒(méi)有直接描述將其應(yīng)用于實(shí)際正演模擬時(shí)如何處理菱形網(wǎng)格.

圖1都是FD交錯(cuò)網(wǎng)格的空間結(jié)構(gòu)，每個(gè)圖中顯示了所有波場(chǎng)速度分量和介質(zhì)參數(shù)的相對(duì)空間位置；在計(jì)算實(shí)現(xiàn)時(shí)，每個(gè)分量或介質(zhì)參數(shù)由2D或3D數(shù)組表示；對(duì)LS方法，如圖1(e、f)，每個(gè)分量或參數(shù)都是菱形，但在這種LS方法中，由于網(wǎng)格為菱形，所以不便于保存這些矩形數(shù)組中的元素或參數(shù)，導(dǎo)致菱形網(wǎng)格有一半的網(wǎng)格點(diǎn)沒(méi)有存儲(chǔ)任何變量，即有一半內(nèi)存占用但是并未使用.本文提出合并速度分量和應(yīng)力分量，形成新的網(wǎng)格.例如，一個(gè)速度分量和一個(gè)應(yīng)力分量合并形成一個(gè)矩形的二維或三維數(shù)組，以保持?jǐn)?shù)組中的每個(gè)元素被充分利用，不額外增加內(nèi)存消耗.根據(jù)我們的研究，菱形網(wǎng)格在實(shí)施過(guò)程中，可以特殊處理：兩個(gè)不同的菱形網(wǎng)格可以合并形成一個(gè)規(guī)則的矩形，以處理所有的波場(chǎng)分量.

這里詳細(xì)說(shuō)明LS的簡(jiǎn)化實(shí)現(xiàn)過(guò)程.基于三維空間中的LS(圖1f)，將所有的σ網(wǎng)格和網(wǎng)格v結(jié)合形成新的網(wǎng)格T.因此在三維模型每個(gè)網(wǎng)格位置都有一個(gè)矩形形狀的T(如圖3所示).這樣一個(gè)組合的條件是，如果i+j+k的和為奇數(shù)，T(i,j,k)代表應(yīng)力分量；如果i+j+k的和為偶數(shù)，則T(i,j,k)代表速度分量.反之亦然，如果i+j+k的和為偶數(shù)代表應(yīng)力分量，i+j+k的和是奇數(shù)代表速度分量.上述表述的優(yōu)點(diǎn)在于，如果方程(5)的三個(gè)表達(dá)式中的T表示速度分量，則其空間導(dǎo)數(shù)即三個(gè)表達(dá)式的結(jié)果正好位于應(yīng)力分量的位置.其優(yōu)點(diǎn)是使用奇數(shù)與偶數(shù)條件，方程(5)可以用來(lái)近似三個(gè)坐標(biāo)軸上的應(yīng)力分量或速度分量的空間導(dǎo)數(shù).

圖3 改進(jìn)的LS數(shù)值實(shí)施方法兩個(gè)網(wǎng)格合并成一個(gè)新網(wǎng)格.極大地簡(jiǎn)化了使用新網(wǎng)格的下標(biāo)而計(jì)算導(dǎo)數(shù).Fig.3 The modified implementation of the LS scheme, where two grids are merged into one This greatly simplifies the calculation of the derivatives with the subscripts of the new grids.

(5)

另外，SSG的實(shí)現(xiàn)涉及到多個(gè)波場(chǎng)分量，而RSG的實(shí)現(xiàn)需要組合多個(gè)波場(chǎng)分量.為了實(shí)現(xiàn)SSG和RSG，須清晰地確定每個(gè)波場(chǎng)分量的空間確切位置；而對(duì)于LS，使用方程(5)和奇數(shù)與偶數(shù)條件，無(wú)需參考空間相對(duì)位置，比如，當(dāng)需要求取(i,j,k)位置的三個(gè)空間導(dǎo)數(shù)，總能通過(guò)所在菱形的6個(gè)頂點(diǎn)元素來(lái)實(shí)現(xiàn)中心差分導(dǎo)數(shù)的求解.基于這種方式，LS方法是這三個(gè)方法中實(shí)現(xiàn)各向異性介質(zhì)模擬的最簡(jiǎn)單有效的方法.

3 模擬流程、震源設(shè)置和吸收邊界條件

3.1 優(yōu)化的模擬流程

鑒于LS法在任意各向異性介質(zhì)模擬方面的優(yōu)越性能，本研究以LS法為基礎(chǔ)，開展三維FD數(shù)值模擬研究.

在三維各向異性介質(zhì)正演模擬時(shí)，計(jì)算量和內(nèi)存開銷是三維正演模擬時(shí)必須關(guān)注的兩個(gè)重要方面.三維各向異性正演模擬多數(shù)應(yīng)用情況下(如基于AVO分析和巖石物理的各向異性研究中)，模型由幾個(gè)層或塊組成.基于層和塊的組成條件，三維模型可使用不同的索引來(lái)區(qū)分不同的層或塊(這里的索引可稱為介質(zhì)構(gòu)造參數(shù)).在模擬計(jì)算過(guò)程中，依據(jù)索引取得介質(zhì)的彈性參數(shù).因此，可將9個(gè)彈性常數(shù)和1個(gè)介質(zhì)密度合并成1個(gè)介質(zhì)構(gòu)造參數(shù)，加上9個(gè)波場(chǎng)變量，共需10個(gè)三維變量，約為原始19個(gè)變量的一半，并且復(fù)雜度顯著降低.

本文提出了一種新的適用于三種交錯(cuò)網(wǎng)格有限差分格式的優(yōu)化正演模擬計(jì)算流程，如圖4所示：首先，定義三種參數(shù)類型：介質(zhì)構(gòu)造參數(shù)(索引)、介質(zhì)彈性參數(shù)和正演模擬參數(shù)；其次，這些參數(shù)被輸入到方程(1)的右邊，以求解速度導(dǎo)數(shù)；然后，用這些速度導(dǎo)數(shù)和前一個(gè)時(shí)間步的應(yīng)力分量求解方程(1)左邊的應(yīng)力分量，用求解的應(yīng)力分量求解方程(2)右邊的導(dǎo)數(shù)，代入方程(2)的左邊，計(jì)算速度分量；最后，在地震檢波器位置輸出速度分量作為地震記錄.在每次循環(huán)過(guò)程中，考慮邊界條件和震源設(shè)置.在圖4的工作流程中，將介質(zhì)參數(shù)分為介質(zhì)構(gòu)造參數(shù)和介質(zhì)彈性參數(shù)的原因是，在建立比較大的三維模型時(shí)，能有效節(jié)省計(jì)算機(jī)內(nèi)存.

圖4 優(yōu)化的三維正演模擬工作流程Fig.4 The optimized 3D modelling workflow

3.2 震源設(shè)置和吸收邊界條件

三個(gè)交錯(cuò)網(wǎng)格方法存在差異，但對(duì)于震源設(shè)置來(lái)說(shuō)，均可在鄰近震源網(wǎng)格處添加Ricker子波的爆炸震源來(lái)實(shí)現(xiàn)(方程(12)，Sheriff，2002)，將子波的振幅添加在正應(yīng)力之上(二維情況正應(yīng)力為σxx和σzz，三維情況正應(yīng)力為σxx,σyy和σzz)，表2給出了具體每種方法震源施加的位置.Ricker子波形式為

r(t)=[1-2(πft)2]e-(πft)2,

(6)

其中，f表示子波主頻，t是時(shí)間，[1-2(πft)2]是子波振幅值.

一旦觸發(fā)震源，地震波將從震源位置開始向四周傳播，如果傳播時(shí)間足夠長(zhǎng)，波場(chǎng)將遇到模型邊界并發(fā)生反射.通常而言，來(lái)自模型邊界的反射是無(wú)用信號(hào)，故須對(duì)邊界能量進(jìn)行處理.我們選擇吸收函數(shù)衰減在二維/三維模型邊界的能量.聲波或各向同性介質(zhì)模擬通常使用PML邊界條件，在本研究中三維各向異性介質(zhì)的模擬涉及到參數(shù)眾多，計(jì)算量巨大，如果使用PML邊界條件，計(jì)算量會(huì)增加數(shù)倍，因此從實(shí)用性的角度出發(fā)，最終選取了常用的“海綿”吸收邊界條件(Cerjan et al.，1985).本研究應(yīng)用“海綿”吸收邊界條件，使用厚度為20個(gè)網(wǎng)格，吸收邊界使用的函數(shù)方程如下：

表2 不同F(xiàn)D格式下應(yīng)用法向應(yīng)力震源子波的網(wǎng)格位置Table 2 Grid locations to be applied to the source wavelet on the normal stresses in different FD schemes

G=e-[0.015(20-i)]2,

(7)

其中，i是離邊界網(wǎng)格點(diǎn)標(biāo)號(hào).

4 實(shí)例比較和三維算法驗(yàn)證

首先在二維情況下實(shí)現(xiàn)這三種方法，以便了解它們?cè)诘卣鹫菽M中的優(yōu)缺點(diǎn).用不同的FD方法在不同的各向異性介質(zhì)中進(jìn)行了一系列二維試驗(yàn)，并且將3D結(jié)果與反射率方法進(jìn)行對(duì)比以驗(yàn)證新方法的精度.

4.1 使用SSG、RSG和LS的二維地震正演模擬

4.1.1 每種方法的潛力

圖5a和5b分別展示了兩個(gè)具有SSG和RSG的各向同性正演模擬示例.所有模型參數(shù)相同，只有方法(SSG或RSG)不同.介質(zhì)參數(shù)為：縱波速度3500 m·s-1，橫波速度2000 m·s-1，密度2.2 g·cm-3.爆炸震源位于模型中心；時(shí)間間隔為1 ms；Ricker子波的主頻為40 Hz；記錄時(shí)間為0.35 s；模型尺寸為600×600，網(wǎng)格尺寸為5 m×5 m.所有的模型參數(shù)相同，但方法不同.在圖5b波前面內(nèi)，出現(xiàn)了許多類似于散射波，這是由數(shù)值頻散引起.正如預(yù)期的那樣，在二維各向同性介質(zhì)的波場(chǎng)快照中縱波的波前是個(gè)完美的圓.圖5b中波前內(nèi)觀察到更多的數(shù)值頻散.其原因是在RSG中，為了保持與SSG或LS中相同的精度，頻散關(guān)系需要更小的網(wǎng)格間距.圖5c顯示了從圖5a和5b中兩個(gè)紅線位置提取的兩條模擬曲線，其中的RSG方法對(duì)應(yīng)的曲線存在數(shù)值頻散現(xiàn)象.

圖5 各向同性介質(zhì)中速度分量Z的波場(chǎng)快照SSG (a)和RSG (b)；(c)提取兩個(gè)波場(chǎng)快照在第400垂直道位置的振幅曲線，繪制形成兩條曲線(黑色曲線表示SSG的振幅，灰色曲線表示RSG的振幅)比較SSG和RSG的地震道，波動(dòng)非常明顯，但SSG準(zhǔn)確性良好.Fig.5 Velocity Z-component snapshots taken with the SSG (a) and RSG (b) in an isotropic medium, and two traces extracted from the locations Trace 400 in (a) and (b) are plotted in (c) The black curve in (c) represents the SSG component amplitude and the gray one in (c) represents the RSG component amplitude.

圖6為兩個(gè)使用SSG和RSG方法進(jìn)行各向異性正演模擬的算例.圖6a和圖6b分別為SSG正演模擬VTI介質(zhì)的X分量和Z分量波場(chǎng)快照.其中模擬參數(shù)是，網(wǎng)格尺寸500×500，爆炸源在每個(gè)模型中心，時(shí)間間隔為1 ms，Ricker 子波的主頻為40 Hz，網(wǎng)格間距為5 m×5 m，記錄時(shí)間為0.3 s，軸的單位為網(wǎng)格數(shù).該VTI模型的Thomsen參數(shù)為α0=2.4495 km·s-1,β0=1.4142 km·s-1,ε=0.3333,δ=0.0885,γ=0.2500和ρ=1 g·cm-3.圖6c和圖6d分別為RSG正演模擬TTI介質(zhì)的X分量和Z分量波場(chǎng)快照.TTI介質(zhì)是由VTI介質(zhì)通過(guò)逆時(shí)針繞Y軸旋轉(zhuǎn)30°而產(chǎn)生.所有其他參數(shù)都相同.各向異性參數(shù)作為Thomsen參數(shù)給出.從圖6c和圖6d，可以在波場(chǎng)快照中清楚地看到旋轉(zhuǎn)角度.在兩個(gè)例子中，VTI和TTI介質(zhì)，無(wú)論是縱波和橫波的觀察，縱波波前長(zhǎng)軸指向系統(tǒng)平面方向.因?yàn)镾SG和RSG的觀測(cè)系統(tǒng)和參數(shù)都相同，而RSG的差分單元長(zhǎng)度是單個(gè)網(wǎng)格對(duì)角線，因而更容易導(dǎo)致數(shù)值頻散，如圖6c和圖6d所示.這兩個(gè)例子表明，SSG能夠有效處理各向同性、VTI、HTI和正交各向異性介質(zhì)，而RSG能夠處理任意復(fù)雜的各向異性介質(zhì)，包括TTI、單斜和三斜介質(zhì).

圖6 RSG方法的速度分量波場(chǎng)快照(a)和(b)為使用SSG得到的VTI介質(zhì)速度的x和z分量； (c)和(d)分別是(a)和(b)的VTI介質(zhì)繞y軸逆時(shí)針旋轉(zhuǎn)30°，即TTI介質(zhì)，使用RSG得到的波場(chǎng)快照.Fig.6 Velocity component snapshots in RSG (a) and (b) are the snapshots of velocity component x and z respectively with the SSG in a VTI medium. (c) and (d) are the snapshots of velocity component x and z respectively with the RSG in the same medium as in the SSG but with a 30 degree anticlockwise rotation applied around the y-axis, namely, a TTI medium.

圖7a和7b分別顯示了VTI介質(zhì)中的X分量和Z分量的波場(chǎng)快照.模型對(duì)應(yīng)的Thomsen參數(shù)為，α0=3.368 km·s-1,β0=1.829 km·s-1,ε=0.11,δ=-0.035,γ=0.255和ρ=2.5 g·cm-3.圖7c和7d分別顯示了在TTI介質(zhì)通過(guò)順時(shí)針45°旋轉(zhuǎn)圖7a和7b的VTI介質(zhì)產(chǎn)生的相應(yīng)的波場(chǎng)快照，模型的尺寸是1000×1000；爆炸震源在每個(gè)模型中心；時(shí)間間隔為1 ms；網(wǎng)格間距為7 m×7 m；Ricker子波的主頻為20 Hz；記錄時(shí)間為0.5 s.圖7中VTI和TTI模型模擬結(jié)果除角度差別外，振幅和形狀一致可以證實(shí)二維LS方法的數(shù)值試驗(yàn)正確．兩個(gè)例子中的頻散關(guān)系相同，所以它們的精度相同.兩個(gè)例子均表明LS能夠模擬任意各向異性介質(zhì).在此，只演示了每種正演模擬方法的優(yōu)勢(shì)與缺點(diǎn)，尚未對(duì)其結(jié)果進(jìn)行定量比較或驗(yàn)證.后續(xù)的三維模擬中，我們將通過(guò)與基于反射率法模擬結(jié)果的比較，來(lái)驗(yàn)證本文方法的有效性.

圖7 使用LS方法的速度分量的波場(chǎng)快照(a)和(b)分別為VTI介質(zhì)中得到的速度x分量和z分量； (c)和(d)為(a)和(b)代表的VTI介質(zhì)繞y軸順時(shí)針旋轉(zhuǎn) 45°得到的TTI介質(zhì).Fig.7 Velocity component snapshots with the LS scheme (a) and (b) are the snapshots of velocity component x and z respectively in a VTI medium. (c) and (d) are the snapshots of velocity component x and z respectively with a 45 degree clockwise rotation applied around the y-axis, namely, a TTI medium.

4.1.2 矩形網(wǎng)格的正演模擬

在某些情況下，數(shù)值模型的網(wǎng)格邊長(zhǎng)不一定相等.更具體地說(shuō)，網(wǎng)格沿x坐標(biāo)的邊長(zhǎng)不等于沿y坐標(biāo)的長(zhǎng)度，也不等于沿z坐標(biāo)的長(zhǎng)度.用四個(gè)基于二維SSG和二維LS的例子演示了正演模擬，但網(wǎng)格邊長(zhǎng)不相等.對(duì)應(yīng)的VTI介質(zhì)模型參數(shù)VTI介質(zhì)參數(shù)用Thomsen參數(shù)表示，α0=2449.5 m·s-1,β0=1414.2 m·s-1,ε=0.3333,δ=0.0885,γ=0.25，ρ=1 g·cm-3.模擬參數(shù)為，采樣間隔0.001 s，記錄長(zhǎng)度0.27 s，主頻率20 Hz的Ricker子波.模擬結(jié)果見圖8.對(duì)于SSG方法，圖8a采用10 m×10 m的網(wǎng)格.在10 m×5 m的網(wǎng)格下，圖8b所示對(duì)應(yīng)的地方被拉伸，正好是垂直長(zhǎng)度的2倍.除拉伸外，圖8a和圖8b中的波場(chǎng)特征都是一樣的，如振幅和相位.這證實(shí)了考慮矩形網(wǎng)格的SSG方法的有效性.LS方法生成圖8c和圖8d，除了模擬方法不一樣之外，模擬和模型參數(shù)與圖8a和圖8b相同，其結(jié)果和圖8a和b具有相同的特征和模式.

圖8 四個(gè)矩形網(wǎng)格的SSG和LS方法數(shù)值試驗(yàn)SSG方法生成(a)和(b)兩個(gè)樣例， LS方法生成(c)和(d)另外兩個(gè)樣例.Fig.8 Four modelling examples with the SSG and LS methods (a) and (b) are generated with the SSG method, and (c) and (d) generated with the LS method.

4.2 三維地震LS正演模擬及算法驗(yàn)證

本文給出了一個(gè)三層模型的數(shù)值實(shí)現(xiàn)來(lái)展示新方法在3D條件下的優(yōu)越性能，同時(shí)與反射率法(Kennett, 1983；Koketsu et al., 1991；Sen and Pal, 2009)進(jìn)行對(duì)比，以驗(yàn)證本方法的正確性.

4.2.1 三維模型和模擬參數(shù)設(shè)置

圖9顯示了3層模型，中間層為HTI(或裂縫)層，其他兩層為各向同性.表3列出了三層參數(shù)的詳細(xì)信息.模型大小是1000 m×1000 m×2000 m，對(duì)應(yīng)的網(wǎng)格數(shù)為101×101×201.震源(圖9a中的加黑點(diǎn))位于網(wǎng)格(0 m, 0 m, 0 m)，位于模型表面的角落.圖9a地震檢波器分布在模型表面的每個(gè)網(wǎng)格位置.其他模擬參數(shù)為，爆炸震源設(shè)置在網(wǎng)格點(diǎn)(0,0,0)，震源主頻為15 Hz的Ricker子波，所有檢波器與震源在同一深度，采樣時(shí)間間隔為1 ms，記錄長(zhǎng)度為2 s.

表3 三維三層模型的參數(shù)Table 3 The parameters of the 3D 3-layer model

正如前文討論過(guò)，層狀或塊狀模型有助于在建立大型3D模型時(shí)節(jié)省計(jì)算機(jī)內(nèi)存和提高計(jì)算效率.在實(shí)際工作中，存在一個(gè)包圍整個(gè)三維模型的吸收邊界層，用來(lái)將能量衰減到邊界以避免邊界的反射.在參數(shù)中沒(méi)有考慮衰減層.震源為主頻為15 Hz的Ricker子波.3D炮道集記錄的時(shí)間長(zhǎng)度為2 s.

4.2.2 結(jié)果與分析

從圖9模型的三維炮道集抽取兩個(gè)z分量剖面進(jìn)行分析，道間距為20 m，圖10表示X-Z剖面，圖11表示Y-Z剖面，兩者都通過(guò)震源位置.X-Z剖面平行于裂縫法向，Y-Z平面平行于裂縫走向.這里，為了便于對(duì)比分析，我們使用反射率法進(jìn)行正演，提取對(duì)應(yīng)位置的剖面.反射率法正演精度高，常用于驗(yàn)證其他正演方法的模擬結(jié)果，但是反射率法通常只適用于水平層狀模型.圖中黑色的波形(圖10和圖11)是我們使用LS正演模擬產(chǎn)生的剖面，紅色的剖面是反射率法所產(chǎn)生的剖面.從圖10和圖11重疊的波形，我們可以看到無(wú)論是振幅還是相位，二者一致性非常好，這驗(yàn)證了我們的3D正演數(shù)值模擬準(zhǔn)確.

圖9 (a) 模型網(wǎng)格數(shù)為101×101×201，網(wǎng)格間距為10 m×10 m×10 m的三維三層模型； (b)為X-Z平面顯示層的厚度在二維模型剖面Fig.9 The 3D three-layer model (a) with the model size 101×101×201 and the grid size 10 m×10 m×10 m. The 2D section (b) in the X-Z plane shows the layers′ thickness

圖10 X-Z平面的z分量合成剖面(a); (b) a圖中藍(lán)色矩形區(qū)域的放大視圖Fig.10 The synthetic section (a) contains the z-component traces in the X-Z plane from both the LS scheme and the reflectivity method； (b) is an zoom-in view of the blue rectangular zone in (a)

圖11 Y-Z平面的z分量合成剖面(a)； (b) a圖中藍(lán)色矩形區(qū)域的放大視圖Fig.11 The synthetic section (a) contains the z-component traces in the Y-Z plane from both the LS scheme and the reflectivity method； (b) is an zoom-in view of the blue rectangular zone in (a)

圖10和圖11中的黑色、藍(lán)色和紅色箭頭分別指向第一界面反射的縱波、第一界面的PS轉(zhuǎn)換波、第二界面的縱波.比較X-Z剖面和Y-Z剖面發(fā)現(xiàn)，在HTI介質(zhì)層的頂部，反射縱波有一個(gè)顯著的幅度差異，特別是在遠(yuǎn)偏移距處(大橢圓處)，揭示了在這兩個(gè)不同的方位，縱波AVO響應(yīng)的差異.為了看清HTI介質(zhì)層頂部的全方位AVO(AVOA)響應(yīng)，在所有檢波器位置處，提取第一個(gè)縱波的最大振幅，然后繪制在圖12.圖中顯示在400～600 m偏移距(入射角范圍大約在25°～30°)范圍內(nèi)AVOA的橢圓變化，響應(yīng)在水平面上，顏色表示振幅從0°到90°，顏色顯示一個(gè)橢圓形的變化，Y軸表示裂縫走向.這表明，橢圓的長(zhǎng)軸變化與裂縫走向一致.在圖10和圖11中，第一個(gè)縱波的傳播時(shí)間相同，而第二個(gè)反射縱波在X-Z剖面(對(duì)應(yīng)裂縫法向)的遠(yuǎn)距處(較小的藍(lán)色橢圓區(qū))傳播時(shí)間稍大一些，這與裂縫法向縱波速度較小的事實(shí)相吻合.

LS正演模擬結(jié)果與反射率法模擬結(jié)果的定性比較如圖10和圖11所示.為了定量地檢查結(jié)果，從圖11a中的z分量中提取了20道(LS和反射率法結(jié)果的第11道到第30道).我們計(jì)算了這兩種情況之間的相關(guān)系數(shù)，圖13顯示了對(duì)比結(jié)果，與0.97以上的相關(guān)系數(shù)有很高的相關(guān)性.兩種方法依然存在差別，除了算法的差別外，部分原因是因?yàn)樵谡鹪吹膶?shí)現(xiàn)時(shí)，LS中震源的添加是在震源點(diǎn)附近的若干個(gè)網(wǎng)格點(diǎn)都施加了應(yīng)力，和反射率法在一點(diǎn)施加震源有區(qū)別.

圖12 HTI層頂部的縱波方位角AVO響應(yīng)Fig.12 P-wave azimuthal AVO response at the top of the HTI layer

圖13 從圖11a中提取的20道地震記錄中觀察到兩種方法的高相關(guān)系數(shù)Fig.13 High correlation coefficients are observed between the 20 extracted traces from figure 11a

5 結(jié)論

地震正演模擬是研究裂縫介質(zhì)地震響應(yīng)的重要方式.本文比較了三種不同類型的交錯(cuò)網(wǎng)格在復(fù)雜裂縫介質(zhì)中模擬時(shí)的優(yōu)點(diǎn)和局限，數(shù)值結(jié)果表明：SSG方法適用于簡(jiǎn)單各向異性裂縫介質(zhì)(各向同性、VTI、HTI和正交各向異性)，對(duì)于復(fù)雜各向異性介質(zhì)(單斜、TTI和三斜介質(zhì))，LS方法優(yōu)于RSG方法，雖然菱形網(wǎng)格特性在數(shù)值實(shí)現(xiàn)中帶來(lái)了不便，但可以通過(guò)本文提出的改進(jìn)方法加以克服.改進(jìn)后的方法極大簡(jiǎn)化了LS的實(shí)現(xiàn)，為開發(fā)基于正演模擬的復(fù)雜各向異性介質(zhì)的方法提供了新的思路，有利于TTI、單斜和三斜介質(zhì)中的RTM和FWI等基于交錯(cuò)網(wǎng)格的數(shù)據(jù)處理方法的應(yīng)用.同時(shí)，對(duì)于大規(guī)模的復(fù)雜三維模型，本文提出了一種優(yōu)化的模擬計(jì)算流程，將介質(zhì)參數(shù)分解為介質(zhì)構(gòu)造參數(shù)和介質(zhì)彈性參數(shù)，可以節(jié)約內(nèi)存消耗并提高計(jì)算效率.三維LS正演模擬實(shí)例與反射率法的結(jié)果具有很好的一致性，驗(yàn)證了提出的三維LS方法的正確性.三維LS模擬工具為研究三維裂縫介質(zhì)中地震波傳播的各種波現(xiàn)象提供了一種新的實(shí)用手段.

致謝感謝審稿專家們和編輯們的寶貴意見.