基于Lebedev網(wǎng)格的TTI介質(zhì)二維三分量正演模擬

劉東洋　彭蘇萍　師素珍　趙太郎

(中國礦業(yè)大學(xué)(北京)煤炭資源與安全開采國家重點(diǎn)實(shí)驗(yàn)室，北京 100083)

1　引言

隨著地震勘探面臨的勘探目標(biāo)越來越復(fù)雜，地震波在各向異性介質(zhì)中的傳播規(guī)律成為研究熱點(diǎn)。沉積巖地層由于受礦物顆粒、裂隙、晶體等的排列和周期性的薄層等因素的影響而普遍呈現(xiàn)各向異性特征。在理論研究中通常使用VTI和HTI介質(zhì)模型分別描述由周期性薄互層和由定向排列的垂直裂隙引起的各向異性介質(zhì)。但是當(dāng)介質(zhì)的本構(gòu)坐標(biāo)系與觀測坐標(biāo)系存在一定角度時就會引起極大的觀測誤差，因此需要根據(jù)Tsvankin等[1]的研究，對VTI介質(zhì)進(jìn)行坐標(biāo)變換，以獲得TTI介質(zhì)的彈性系數(shù)矩陣。

常采用數(shù)值模擬方法研究地震波傳播規(guī)律，有限差分法是其中最常用的方法。Alterman等[2]最早給出了均勻各向同性介質(zhì)的二維二階波動方程有限差分算法。Virieux[3，4]基于一階速度—應(yīng)力方程提出了交錯網(wǎng)格方法，將應(yīng)力和速度分別定義在不同的網(wǎng)格體系上，并對時間進(jìn)行交錯，與常規(guī)網(wǎng)格法相比差分精度提高了4倍。Dalain[5]給出了求解聲波方程的高階差分方法，該方法降低了單位波長內(nèi)的網(wǎng)格點(diǎn)數(shù)，提高了計(jì)算效率。Graves[6]提出了三維交錯網(wǎng)格有限差分模擬方法。董良國等[7]提出了一階彈性波波動方程交錯網(wǎng)格高階有限差分解法，同時給出了相應(yīng)的穩(wěn)定性條件。裴正林等[8-10]采用高階交錯網(wǎng)格有限差分算法對三維各向異性介質(zhì)進(jìn)行數(shù)值模擬，并推導(dǎo)了三維各向異性介質(zhì)完全匹配層(PML)邊界條件公式和相應(yīng)的交錯網(wǎng)格差分格式。李東等[11]利用交錯網(wǎng)格有限差分法對煤儲層雙相各向異性介質(zhì)進(jìn)行數(shù)值模擬，并用于煤層氣富集區(qū)。綜上所述，對于各向異性介質(zhì)的研究多采用交錯網(wǎng)格方法，但是在交錯網(wǎng)格中每個變量的不同分量都是交錯定義的，對于沒有定義的點(diǎn)需要進(jìn)行變量插值，從而降低了模擬精度。為此，筆者在前人的研究基礎(chǔ)上，給出了TTI介質(zhì)的二維三分量應(yīng)力—速度彈性波方程的Lebedev網(wǎng)格有限差分解法，并對地震波在單層各向異性介質(zhì)以及含透鏡體的復(fù)雜介質(zhì)中的傳播特征進(jìn)行了數(shù)值模擬，驗(yàn)證了Lebedev網(wǎng)格有限差分算法的有效性。

2　方法原理

2.1　TTI介質(zhì)二維三分量應(yīng)力—速度方程

(1)

式中：v、σ分別代表速度和應(yīng)力；Cij為TTI介質(zhì)彈性張量系數(shù)矩陣的元素。

2.2　TTI介質(zhì)Lebedev網(wǎng)格有限差分格式

在使用常規(guī)交錯網(wǎng)格進(jìn)行數(shù)值模擬時，由于需要對未定義在交錯位置上的變量進(jìn)行復(fù)雜的插值計(jì)算，極大地降低了運(yùn)行效率，同時引入了數(shù)值誤差，降低了運(yùn)算精度[12，13]。為此，人們基于Lebedev的思想提出了Lebedev網(wǎng)格[14-16]，該網(wǎng)格的特點(diǎn)是將應(yīng)力、速度分量定義在同一網(wǎng)格上，但是應(yīng)力、速度分量在網(wǎng)格點(diǎn)上交錯分布，相當(dāng)于把需要插值的分量直接在網(wǎng)格點(diǎn)上定義，可以準(zhǔn)確計(jì)算每個變量在各個方向的數(shù)值，而無需進(jìn)行插值。根據(jù)Lebedev網(wǎng)格變量示意圖(圖1)，可對式(1)進(jìn)行網(wǎng)格離散，由此得到Lebedev網(wǎng)格機(jī)制下TTI介質(zhì)波動方程的有限差分格式

(2a)

(2b)

圖1　Lebedev網(wǎng)格變量分布示意圖

2.3　穩(wěn)定條件

在數(shù)值模擬過程中，如果選取的離散參數(shù)不合適，那么數(shù)值計(jì)算結(jié)果會出現(xiàn)巨大偏差，嚴(yán)重時會造成數(shù)值溢出而使計(jì)算無法進(jìn)行，因此穩(wěn)定性問題是數(shù)值模擬過程中求解波動方程的重要問題。對于Lebedev網(wǎng)格有限差分技術(shù)，Lisitas等[16]給出了時間二階、空間2N階差分精度下的Lebedev網(wǎng)格穩(wěn)定性條件

(3)

式中： Δt為時間步長；vmax=max(vi)(i=1、2、3)表示三種波在TTI介質(zhì)中傳播的最大速度；Δx、Δy、Δz分別為對應(yīng)方向的空間步長。

由于Lebedev網(wǎng)格是在傳統(tǒng)交錯網(wǎng)格未定義的網(wǎng)格點(diǎn)上增加了相應(yīng)變量，因此這兩種網(wǎng)格的穩(wěn)定性條件基本相同，但是在理論上Lebedev網(wǎng)格穩(wěn)定性條件的適用性更廣泛[18-21]。

2.4　吸收邊界條件

在數(shù)值模擬過程中，由于計(jì)算空間有限，所以使模擬波場在傳播到計(jì)算邊界時會發(fā)生反射現(xiàn)象，增加了波場干擾，從而影響模擬效果。因此需要引入人工邊界，使波場在邊界附近的能量逐漸衰減。由于在特定情況下傳統(tǒng)的PML吸收邊界條件會因?yàn)榫哂兄笖?shù)增長解而存在不穩(wěn)定現(xiàn)象，因此本文采用多軸完全匹配層(M-PML)進(jìn)行邊界處理。M-PML吸收邊界條件由傳統(tǒng)的PML吸收邊界條件發(fā)展而來，傳統(tǒng)的PML吸收邊界條件方法通過引入衰減因子d、分裂變量而實(shí)現(xiàn)。M-PML吸收邊界條件方法是對傳統(tǒng)PML吸收邊界條件進(jìn)行一般化擴(kuò)展，并在多個方向上同時引入衰減因子。以垂直x軸的匹配層中的vx為例，有

(4)

式中：vx1、vx2分別為垂直于邊界和平行于邊界的vx；dx、dz為衰減因子，dx的取法與傳統(tǒng)PML中的衰減因子相同，dz=P1*dx(P1為穩(wěn)定性因子)。

同理，在垂直y軸的匹配層中存在穩(wěn)定性因子P2，穩(wěn)定性因子的取值通常根據(jù)經(jīng)驗(yàn)獲得，不同介質(zhì)模型對應(yīng)的穩(wěn)定性因子不同，本文的穩(wěn)定性因子取為P1=P2=0.32。特別地，當(dāng)穩(wěn)定性因子都為零時，M-PML吸收邊界退化為常規(guī)PML吸收邊界[17]。

3　正演模擬

3.1　單層介質(zhì)模型

為了驗(yàn)證Lebedev網(wǎng)格算法的有效性及穩(wěn)定性，通過模擬TTI介質(zhì)分析Lebedev網(wǎng)格有限差分技術(shù)的模擬效果與波場特征。隨著主頻升高會使圖像分辨率增強(qiáng)，但是過高的主頻則會產(chǎn)生嚴(yán)重的數(shù)值頻散現(xiàn)象。同時子波主頻的選取受空間和時間采樣間隔的影響，即：當(dāng)采樣間隔較小時需要選取較高的子波主頻以提高波形分辨率；當(dāng)采樣間隔較大時可以采用較低的子波主頻以壓制頻散效果。因此在保證模擬質(zhì)量的前提下，為了盡可能地提高圖像分辨率，采用主頻為100Hz的Ricker子波，震源位于模型中心，其他物性參數(shù)如表1所示，其中Cij表示彈性系數(shù)， φ為方位角， θ為極化角。

表1　TTI介質(zhì)物性參數(shù)

注：Cij的單位為109kg·m-1·s-2, 密度ρ的單位為g·cm-3

圖2為TTI介質(zhì)二維三分量波場快照和炮記錄。在圖中可以清晰地分辨TTI介質(zhì)中的qP波、qSV波、qSH波，其中x分量(圖2a)和z分量(圖2c)的縱波能量較強(qiáng)，y分量(圖2b)的縱波能量較弱，在TTI介質(zhì)中存在橫波分裂現(xiàn)象，因此三個分量的橫波在兩個對角線方向出現(xiàn)分叉；從炮記錄中可以看出，縱波速度最大，qSH波速度最小，這些特征皆符合波場傳播規(guī)律，因此Lebedev網(wǎng)格很好地模擬了地震波在TTI介質(zhì)中的傳播特征(圖2d～圖2f)；地震波沿對稱軸方向的傳播速度明顯小于垂直于對稱軸方向，這也從側(cè)面驗(yàn)證了Lebedev網(wǎng)格算法的有效性。

3.2　吸收邊界效果分析

為測試M-PML吸收邊界處理效果，采用單層TTI介質(zhì)進(jìn)行波場模擬。采用爆炸震源激發(fā)，震源位于模型中心，主頻為80Hz，其他物性參數(shù)如表1所示。圖3為未加M-PML吸收邊界條件的TTI介質(zhì)二維三分量波場快照和炮記錄。由圖可見，qP波在邊界處產(chǎn)生了強(qiáng)反射現(xiàn)象，由于震源的埋深大于四分之一地震波長，因此在邊界位置還產(chǎn)生了qP波虛反射(圖3b、圖3c)，嚴(yán)重干擾了數(shù)值模擬結(jié)果。圖4為加M-PML吸收邊界條件的TTI介質(zhì)二維三分量波場快照和炮記錄，可見邊界處的反射波和虛反射全部消失。圖5為PML與M-PML吸收邊界的x、z分量波場快照。 由圖可見，橫波在PML吸收邊界的左上方出現(xiàn)不穩(wěn)定波場反射(圖5a)，在M-PML吸收邊界條件的波場傳播較為穩(wěn)定，未出現(xiàn)雜亂反射(圖5b)。因此M-PML吸收邊界條件的處理效果較理想。

圖2　TTI介質(zhì)二維三分量波場快照(t=40ms)和炮記錄

圖3　未加M-PML吸收邊界條件的TTI介質(zhì)二維三分量波場快照(t=100ms)和炮記錄

圖5　PML(a)與M-PML吸收邊界(b)的x(左)、z分量(右)波場快照(t=80ms)

圖6　透鏡體介質(zhì)模型示意圖

模型尺寸為400m×400m，網(wǎng)格分為兩部分，白色區(qū)域?yàn)門TI介質(zhì)，黑色部分為透鏡體(各向同性介質(zhì))，x、z方向的空間采樣間隔均為1m，時間采樣間隔為0.1ms，采用爆炸震源激發(fā)，主頻為80Hz，震源位于中心

3.3　透鏡體模型

在沉積巖發(fā)育地區(qū)由于斷層、風(fēng)化和侵蝕等作用的影響，在地下介質(zhì)中可形成一種中間厚、邊緣薄的區(qū)域性地帶，因形似一枚透鏡，故稱作透鏡體。透鏡體一般是由礦物變質(zhì)形成的，因此研究透鏡體模型(圖6)可為礦產(chǎn)資源開發(fā)提供一定的模型基礎(chǔ)。采用非均勻TTI介質(zhì)進(jìn)行波場模擬，具體物性參數(shù)見表2。圖7為非均勻TTI介質(zhì)在40、50、60ms的x、y、z分量波場快照。由圖可見：在40ms時，縱波在透鏡體處發(fā)生了明顯的反射，可以觀察到透鏡體頂面的反射波(圖7a)；在50ms和60ms時，橫波在透鏡體處也發(fā)生反射，同時出現(xiàn)了縱波在透鏡體底面的反射波，在透鏡體的兩端出現(xiàn)繞射波(圖7b、圖7c)，因此基于Lebedev網(wǎng)格得到的波場快照符合地震波傳播規(guī)律。

表2　透鏡體模型物性參數(shù)

注：Cij單位為109kg·m-1·s-2, 密度ρ的單位為g·cm-3

圖7　非均勻TTI介質(zhì)在40ms(a)、50ms(b)、60ms(c)的x(左)、y(中)、z分量(右)波場快照

3.4　Lebedev網(wǎng)格與傳統(tǒng)交錯網(wǎng)格對比

為驗(yàn)證Lebedev網(wǎng)格的運(yùn)行優(yōu)勢，應(yīng)用雙層介質(zhì)模型(圖8)從運(yùn)行效率和差分精度兩方面對Lebedev網(wǎng)格與傳統(tǒng)交錯網(wǎng)格進(jìn)行對比、分析。圖9、圖10分別為由傳統(tǒng)交錯網(wǎng)格有限差分與Lebedev網(wǎng)格有限差分得到的雙層各向同性介質(zhì)在0.4s的x、z分量波場快照，運(yùn)行時間分別為497.107s和362.068s，故后者的運(yùn)行效率提高了27%。由圖可見：由于數(shù)值插值的影響，前者存在明顯的數(shù)值頻散，模擬精度較低，地震波在下層介質(zhì)中的波形顯示也較模糊(圖9)；后者基本不存在數(shù)值頻散，波形顯示也較清晰，模擬精度明顯提高(圖10)。

圖8　雙層介質(zhì)模型示意圖

設(shè)定雙層的各向同性介質(zhì)模型，模型尺寸為400m×400m，模型分為上、下兩層，下層介質(zhì)的速度和密度較大。x、z方向的空間采樣間隔均為10m，時間采樣間隔為1ms，采用爆炸震源激發(fā)，震源位于中心，主頻為30Hz

圖9　雙層各向同性介質(zhì)傳統(tǒng)交錯網(wǎng)格有限差分在0.4s的x(a)、z(b)分量波場快照

圖10　雙層各向同性介質(zhì)Lebedev網(wǎng)格有限差分在0.4s的x(a)、z(b)分量波場快照

4　結(jié)束語

文中基于Lebedev網(wǎng)格對TTI介質(zhì)二維三分量的應(yīng)力—速度方程進(jìn)行了高精度差分離散處理，并分別對單層和復(fù)雜TTI介質(zhì)模型進(jìn)行了正演模擬，得到以下認(rèn)識：①由于Lebedev網(wǎng)格對各向異性介質(zhì)的彈性波方程做離散時無需進(jìn)行波場插值，與傳統(tǒng)交錯網(wǎng)格有限差分方法相比，模擬精度更高；對TTI介質(zhì)進(jìn)行模擬時可以清晰地觀察到縱波、快橫波和慢橫波，并且快、慢橫波的偏振方向相反，在單炮記錄中觀測到的三種波的速度特征也符合波場傳播規(guī)律。②引入多軸完全匹配層(M-PML)吸收邊界條件后，在不影響模擬效果的情況下邊界反射現(xiàn)象被有效地壓制，因此M-PML吸收邊界是一種有效的邊界條件。③在復(fù)雜介質(zhì)條件下，使用Lebedev網(wǎng)格的模擬結(jié)果依然符合實(shí)際情況，因此Lebedev網(wǎng)格適用于各向異性介質(zhì)模型的波場數(shù)值模擬。

[1]Tsvankin I,Thomsen L.Inversion of reflection traveltimes for transverse isotropy.Geophysics,1995,60(4):1095-1107.

[2]Alterman Z,Karal F C.Propagation of elastic waves in layered media by finite difference methods.Bulletin of the Seismological Society of America,1968,58(1):367-398.

[3]Virieux J.SH wave propagation in heterogeneous media: velocity-stress finite-difference method.Geophysics,1984,49(11):1933-1942.

[4]Virieux J.P-SV wave propagation in heterogeneous media: velocity-stress finite-difference method.Geophysics,1986,51(4):889-901.

[5]Dablain M A.The application of high-order differencing to scalar wave equation.Geophysics,1986,51(1):54-66.

[6]Graves R W.Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences.Bulletin of the Seismological Society of America,1996,86(4):1091-1106.

[7]董良國,馬在田,曹景忠等.一階彈性波方程交錯網(wǎng)格高階差分解法.地球物理學(xué)報,2000,43(3):411-419.

Dong Liangguo,Ma Zaitian,Cao Jingzhong et al.Study on stability of staggered grid high order difference method for first order elastic wave equation.Chinese Journal of Geophysics,2000,43(3):411-419.

[8]裴正林.三維各向異性介質(zhì)中彈性波方程交錯網(wǎng)格高階有限差分法數(shù)值模擬.中國石油大學(xué)學(xué)報(自然科學(xué)版),2004,28(5):23-29.

Pei Zhenglin.Numerical simulation of elastic wave equation with staggered grid high-order finite difference method in three dimensional anisotropic medium.Journal of China University of Petroleum(Edition of Natural Science),2004,28(5):23-29.

[9]裴正林,董玉珊,彭蘇萍.裂隙煤層彈性波場方位各向異性特征數(shù)值模擬研究.石油地球物理勘探,2007,42(6):665-672.

Pei Zhenglin,Dong Yushan,Peng Suping.Numerical simulation of azimuthal anisotropy of elastic wave field in fractured seams.OGP,2007,42(6):665-672.

[10]裴正林.層狀各向異性介質(zhì)中橫波分裂和再分裂數(shù)值模擬.石油地球物理勘探,2006,41(1):17-25.

Pei Zhenglin.Numerical simulation of splitting and splitting of transverse wave in layered anisotropic media.OGP,2006,41(1):17-25.

[11]李東,董守華,趙小翠等.煤儲層雙相EDA介質(zhì)的地震波場模擬.地球物理學(xué)進(jìn)展,2011,26(2):654-663.

Li Dong,Dong Shouhua,Zhao Xiaocui et al.Simulation of seismic wave field in coal reservoir with dual phase EDA medium.Progress in Geophysics,2011,26(2):654-663.

[12]王宏偉,彭蘇萍,杜文鳳等.HTI構(gòu)造煤方位AVO正演.石油地球物理勘探,2014,49(6):1122-1130.

Wang Hongwei,Peng Suping,Du Wenfeng et al.HTI structural coal orientation AVO forward.OGP,2014,49(6):1122-1130.

[13]鎮(zhèn)晶晶,劉洋,李敏.幾種裂縫模型的波場傳播特征比較.石油地球物理勘探,2012,47(6):908-917.

Zhen Jingjing,Liu Yang,Li Min.Comparison of wave field characteristics among the fracture rock models.OGP,2012,47(6):908-917.

[14]Lebedev V I.Difference analogues of orthogonal decompositions,basic differential operators and some boundary problems of mathematical physics I.USSR Computational Mathematics and Mathematical Phy-sics,1964,4(3):69-92.

[15]Bernth H,Chapman C.A comparison of finite-difference grids for anisotropic elastic modeling.SEG Technical Program Expanded Abstracts,2010,29:2916-2920.

[16]Lisitsa V,Vishnevskiy D.Lebedev scheme for the numerical simulation of wave propagation in 3D anisotropic elasticity.Geophysical Prospecting,2010,58(4):619-635.

[17]林朋,彭蘇萍,盧勇旭等.基于旋轉(zhuǎn)交錯網(wǎng)格的雙相各向異性介質(zhì)二維三分量波場模擬.煤炭學(xué)報,2016,41(5):1203-1211.

Lin Peng,Peng Suping,Lu Yongxu et al.Study on 2D/3C wave propagation in two-phase anisotropic media using the rotated staggered-grid method.Journal of China Coal Society,2010,58(4):619-635.

[18]李娜,黃建平,李振春等.Lebedev網(wǎng)格改進(jìn)差分系數(shù)TTI介質(zhì)正演模擬方法研究.地球物理學(xué)報,2014,57(1):261-269.

Li Na,Huang Jianping,Li Zhenchun et al.Research on Lebedev mesh modified differential coefficient TTI forward modeling method.Chinese Journal of Geophysics,2014,57(1):261-269.

[19]李娜,李振春,黃建平等.Lebedev網(wǎng)格與標(biāo)準(zhǔn)交錯網(wǎng)格耦合機(jī)制下的復(fù)雜各向異性正演模擬.石油地球物理勘探,2014,49(1):121-131.

Li Na,Li Zhenchun,Huang Jianping et al.Complex anisotropic forward modeling based on Lebedev mesh and standard staggered mesh coupling mechanism.OGP,2014,49(1):121-131.

[20]黃建平,楊宇,李振春等.基于M-PML邊界的Lebedev網(wǎng)格起伏地表正演模擬方法及穩(wěn)定性分析.中國

石油大學(xué)學(xué)報(自然科學(xué)版),2016,40(4):47-56.

Huang Jianping,Yang Yu,Li Zhenchun et al.Forward modeling method and stability analysis of Lebedev grid based on M-PML boundary.Journal of China University of Petroleum(Edition of Natural Science),2016,40(4):47-56.

[21]楊宇,黃建平,雷建設(shè)等.Lebedev網(wǎng)格黏彈性介質(zhì)起伏地表正演模擬.石油地球物理勘探,2016,51(4):698-706.

Yang Yu,Huang Jianping,Lei Jianshe et al.Forward modeling of viscoelastic medium in Lebedev lattice.OGP,2016,51(4):698-706.

[22]楊慶節(jié),劉財,耿美霞等.交錯網(wǎng)格任意階導(dǎo)數(shù)有限差分格式及差分系數(shù)推導(dǎo).吉林大學(xué)學(xué)報,2014,44(1):375-385.

Yang Qingjie,Liu Cai,Geng Meixia et al.Staggered grid finite difference scheme and coefficients deduction of any number of derivatives.Journal of Jilin University(Earth Science Edition),2014,44(1):375-385.

[23]王洋,劉洪,張衡等.一種全局優(yōu)化的隱式交錯網(wǎng)格有限差分算法及其在彈性波數(shù)值模擬中的應(yīng)用.地球物理學(xué)報.2015,58(7):2508-2524.

Wang Yang,Liu Hong,Zhang Heng et al.A global optimized implicit staggered grid finite difference algorithm and its application in numerical simulation of elastic wave.Chinese Journal of Geophysics,2015,58(7):2508-2524.

[24]尹志恒,狄?guī)妥?魏建新等.裂隙參數(shù)對縱波能量衰減影響的物理模型研究.石油地球物理勘探,2012,47(5):728-734.

Yin Zhiheng,Di Bangrang,Wei Jianxin et al.P-wave attenuation by fracture parameter on physical models.OGP,2012,47(5):728-734.

[25]沈金松,詹林森,馬超.裂縫等效介質(zhì)模型對裂縫結(jié)構(gòu)和充填介質(zhì)參數(shù)的適應(yīng)性.吉林大學(xué)學(xué)報:地球科學(xué)版,2013,43(3):993-1003.

Shen Jinsong,Zhan Linsen,Ma Chao.The applicability of fracture effective model to diffeerent parameters of fracture geometric strutures and media filled in fracture pores.Journal of Jinlin University(Earth Science Edition),2013,43(3):993-1003.

[26]姚振岸,孫成禹,謝俊法等.黏彈TTI介質(zhì)旋轉(zhuǎn)交錯網(wǎng)格微地震波場模擬.石油地球物理勘探,2017,52(2):253-263.

Yao Zhen’an,Sun Chengyu,Xie Junfa et al.Simulation of micro-seismic wave field in viscous TTI medium rotating staggered.OGP,2017,52(2):253-263.