基于Lebedev網格的TTI介質二維三分量正演模擬

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

(中國礦業大學(北京)煤炭資源與安全開采國家重點實驗室，北京 100083)

1　引言

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

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

2　方法原理

2.1　TTI介質二維三分量應力—速度方程

(1)

式中：v、σ分別代表速度和應力；Cij為TTI介質彈性張量系數矩陣的元素。

2.2　TTI介質Lebedev網格有限差分格式

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

(2a)

(2b)

圖1　Lebedev網格變量分布示意圖

2.3　穩定條件

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

(3)

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

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

2.4　吸收邊界條件

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

(4)

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

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

3　正演模擬

3.1　單層介質模型

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

表1　TTI介質物性參數

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

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

3.2　吸收邊界效果分析

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

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

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

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

圖6　透鏡體介質模型示意圖

模型尺寸為400m×400m，網格分為兩部分，白色區域為TTI介質，黑色部分為透鏡體(各向同性介質)，x、z方向的空間采樣間隔均為1m，時間采樣間隔為0.1ms，采用爆炸震源激發，主頻為80Hz，震源位于中心

3.3　透鏡體模型

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

表2　透鏡體模型物性參數

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

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

3.4　Lebedev網格與傳統交錯網格對比

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

圖8　雙層介質模型示意圖

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

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

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

4　結束語

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

[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]董良國,馬在田,曹景忠等.一階彈性波方程交錯網格高階差分解法.地球物理學報,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]裴正林.三維各向異性介質中彈性波方程交錯網格高階有限差分法數值模擬.中國石油大學學報(自然科學版),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]裴正林,董玉珊,彭蘇萍.裂隙煤層彈性波場方位各向異性特征數值模擬研究.石油地球物理勘探,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]裴正林.層狀各向異性介質中橫波分裂和再分裂數值模擬.石油地球物理勘探,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介質的地震波場模擬.地球物理學進展,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構造煤方位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]鎮晶晶,劉洋,李敏.幾種裂縫模型的波場傳播特征比較.石油地球物理勘探,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]林朋,彭蘇萍,盧勇旭等.基于旋轉交錯網格的雙相各向異性介質二維三分量波場模擬.煤炭學報,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網格改進差分系數TTI介質正演模擬方法研究.地球物理學報,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網格與標準交錯網格耦合機制下的復雜各向異性正演模擬.石油地球物理勘探,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網格起伏地表正演模擬方法及穩定性分析.中國

石油大學學報(自然科學版),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]楊宇,黃建平,雷建設等.Lebedev網格黏彈性介質起伏地表正演模擬.石油地球物理勘探,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]楊慶節,劉財,耿美霞等.交錯網格任意階導數有限差分格式及差分系數推導.吉林大學學報,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]王洋,劉洪,張衡等.一種全局優化的隱式交錯網格有限差分算法及其在彈性波數值模擬中的應用.地球物理學報.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]尹志恒,狄幫讓,魏建新等.裂隙參數對縱波能量衰減影響的物理模型研究.石油地球物理勘探,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]沈金松,詹林森,馬超.裂縫等效介質模型對裂縫結構和充填介質參數的適應性.吉林大學學報:地球科學版,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介質旋轉交錯網格微地震波場模擬.石油地球物理勘探,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.