陳玉璽，王 璐，章 陽

（1.中國海洋大學，山東青島266100；2.江西地震局，江西南昌330039）

三維粘彈性介質(zhì)地震波場數(shù)值模擬技術(shù)研究

陳玉璽*1，王 璐1，章 陽2

（1.中國海洋大學，山東青島266100；2.江西地震局，江西南昌330039）

為了適應油氣田勘探的發(fā)展需求，提高勘探精度，地震勘探由二維向三維發(fā)展。傳統(tǒng)的地震波數(shù)值模擬一般假設地層介質(zhì)是均勻、完全彈性的，但在實際情況中，由于介質(zhì)是非完全彈性、各向異性、物性參數(shù)可連續(xù)變化的，因此地震波的能量在傳播過程中是逐漸衰減的。在非均勻介質(zhì)情況下進行地震波的數(shù)值模擬，可以更真實地反映地震波的傳播和能量變化規(guī)律。推導出三維交錯網(wǎng)格高階有限差分算法，求取優(yōu)化差分算子，實現(xiàn)三維粘彈性波數(shù)值模擬，獲得地震記錄，并分析其波場特征和能量變化情況。

非完全彈性；交錯網(wǎng)格高階有限差分；差分算子；數(shù)值模擬

地震波數(shù)值模擬技術(shù)是研究復雜油氣田地球物理勘探最基本也是最重要的一種方法，是反演的基礎(chǔ)。數(shù)值模擬的精度也決定了反演的精度，對于我們認識地質(zhì)構(gòu)造，油氣田勘探，識別裂隙型油藏具有重要的指導作用。三維復雜粘彈性介質(zhì)的數(shù)值模型是地震波數(shù)值模擬的主要工作，也是當今地震勘探的重點。

1 三維粘彈性波方程的建立

目前使用最廣泛的粘彈性介質(zhì)模型為佛各特（Viogt）模型[1]，它將粘彈性介質(zhì)的應力分為2部分：完全彈性應力和粘性應力。

假設λ和μ為完全彈性拉梅常數(shù)，λ′和μ′為粘性拉梅常數(shù)，則粘彈性波動方程中的拉梅常數(shù)等效為將完全彈性波方程中的拉梅常數(shù)等效替換成粘彈性波方程中的拉梅常數(shù)，就得到了粘彈性波方程。本文主要研究一階速度-應力形式的三維粘彈性波方程，在均勻各向同性介質(zhì)中（假設體力為0），其表達式為式（1）。

其中，λ和μ是彈性拉梅常數(shù)，λ′和μ′為是粘性拉梅常數(shù)，我們知道根據(jù)前人對粘彈性介質(zhì)理論的研究有：由這幾個式子可以得出式（1）中的拉梅常數(shù)：

其中，ρ是介質(zhì)的密度，vp和vs是介質(zhì)的縱波和橫波速度，Qp和Qs是縱、橫波品質(zhì)因子，w是圓頻率。

2 交錯網(wǎng)格高階有限差分法

交錯網(wǎng)格法在常規(guī)網(wǎng)格的基礎(chǔ)上引入半網(wǎng)格點，可以有效地處理速度分量和應力分量之間的耦合關(guān)系，提高數(shù)值模擬的精度和穩(wěn)定性。有限差分法中用離散化的差商近似代替連續(xù)的微商，這樣就會產(chǎn)生無法避免的數(shù)值頻散，而高階有限差分法可以有效地降低由于網(wǎng)格離散化造成的數(shù)值頻散，提高數(shù)值模擬的精度。本文用交錯網(wǎng)格高階有限差分法對三維粘彈性波方程進行數(shù)值模擬，其網(wǎng)格剖分如圖1所示。

圖1 三維粘彈性介質(zhì)交錯網(wǎng)格示意圖

其中將正應力分量σxx、σyy、σzz定義在點( )i,j,k上；剪應力分量σxy定義在點上，σxz定義在點上，σyz定義在點上；速度分量vx和速度變化率分量定義在點上，vy和定義在點上，vz和定義在點上。

在地震波場數(shù)值模擬中，有幾個問題是需要著重考慮的，首先是穩(wěn)定性條件，本文采用裴正林等提出的穩(wěn)定性條件，其表達式為[3]：

式中：Δx、Δy、Δz——空間x、y、z三個方向的網(wǎng)格間距；

Δt——時間網(wǎng)格間距；

vmax——模型最大速度；

am——高階差分系數(shù)，本文選取優(yōu)化差分算子（以時間2階，空間10階為例）：a1=1.250，a2=-0.1200，a3=0.0314，a4=-0.0092，a5=0.0018。

其次是吸收邊界條件，實際地震勘探中，地震波往往是在無限空間的地下地層中傳播。而在用計算機數(shù)值模擬時，我們所選擇的理論物理模型通常是在有限區(qū)域范圍內(nèi)，這樣便產(chǎn)生了人工邊界反射問題[4]。本文采用PML吸收邊界條件實現(xiàn)對三維粘彈性波方程的邊界處理[5]。PML的核心思想是在計算區(qū)域外面構(gòu)造有限厚度的吸收層，吸收或衰減向外傳播的波，同時在計算模型區(qū)域與吸收層的鏈接處是透明的，使得產(chǎn)生最小可能的虛假反射來解決模擬時的人工邊界問題。

三維模型需要在6個面、12條棱和8個角點處分別添加PML吸收邊界。在x、y、z三個方向均加載衰減項d(x)、d(y)、d(z)，其計算公式為：

其中，R為理論反射系數(shù)，δ為PML厚度[2]。在8個角點區(qū)域，d(x)≠0，d(y)≠0，d(z)≠0；在12條棱區(qū)域，xoy平面和 yoz平面相交的棱域，d(x)≠0，d(y)=0，d(z)≠0，xoy平面與xoz平面相交的棱域，d(x)=0，d(y)≠0，d(z)≠0，xoz平面與yoz平面相交的棱域，d(x)≠0，d(y)≠0，d(z)=0；在6個面區(qū)域，垂直x軸的2個平面，d(x)≠0，d(y)=0，d(z)=0，垂直y軸的2個平面，d(x)=0，d(y)≠0,d(z)=0,垂直z軸的2個平面，d(x)=0,d(y)=0,d(z)≠0。這樣整個模型區(qū)域都引入PML吸收邊界，地震波傳播經(jīng)過模型邊界到達PML吸收層時，能量逐漸減弱，因此不會產(chǎn)生任何邊界反射。

以vx和σxx為例，根據(jù)PML的分裂思想，則三維粘彈性波方程的精度為Ο(Δ t2+Δx2N)的引入PML吸收邊界條件的交錯網(wǎng)格高階差分格式為[6]：

其中：

其中：

3 三維復雜粘彈性介質(zhì)模型

在三維French速度模型的基礎(chǔ)上[7]，加入縱橫波品質(zhì)因子Qp、Qs從而構(gòu)建三維French粘彈性介質(zhì)速度模型。模型包含1個水平層面、1個傾斜界面和2個隆起構(gòu)造，整個模型σ=0.25，ρ=2000kg/m3。模型大小為1000m × 1000m × 1000m，空 間 網(wǎng) 格 步 長Δx=Δy=Δz=5m，時間步長Δt=0.5ms，震源子波選擇主頻為30Hz的雷克子波，計算精度為O(Δt2+Δx10)，圖2為相應各個平面的二維垂直切片的縱波速度模型。

圖2 各平面上的二維垂直剖面的縱波速度

圖3、圖4、圖5分別為t=500ms時，不同平面上的波場快照，分析各個平面上的粘彈性波波場快照切片，由于三維波場之間的相互干涉和相互影響，形成的波場十分復雜。xoz平面上，隆起構(gòu)造的側(cè)面反射波和透射波明顯，能量強，yoz平面上同樣能觀察到能量很強的側(cè)面反射波和透射波，xoy平面上能觀察到2個隆起構(gòu)造產(chǎn)生的垂直反射波信息和多次側(cè)面反射波信息，還產(chǎn)生能量較強的散射波。隨著傳播深度的不斷增大，由于存在粘滯吸收衰減作用，能量產(chǎn)生一定的衰減，而且記錄到的散射波的能量明顯比反射波弱。

圖3 t=500ms時刻y=500m處xoz平面上的粘彈性波波場快照

4 結(jié)論

以French速度模型為基礎(chǔ)，加入縱橫波品質(zhì)因子Qp、Qs構(gòu)建三維粘彈性介質(zhì)速度模型，以粘彈性介質(zhì)為例，根據(jù)其應力-應變關(guān)系，結(jié)合完全彈性波方程和表征介質(zhì)吸收衰減特性的品質(zhì)因子Q，推導出三維粘彈性波方程，并進一步給出其一階速度-應力關(guān)系式。實驗發(fā)現(xiàn)，地震波在粘彈性介質(zhì)中傳播時存在明顯的能量衰減現(xiàn)象，地震波的振幅明顯衰減，反射波同相軸能量明顯減弱，同時隨著地震波傳播深度的不斷增加，能量衰減越發(fā)嚴重。反射波波形發(fā)生畸變，而且吸收衰減作用越明顯，地震波波形畸變越嚴重。并且這種粘滯吸收衰減作用對橫波影響比對縱波強。同時本文研究的三維粘彈性交錯網(wǎng)格高階有限差分法模擬技術(shù)不僅能夠更詳細全面地描述地震波的傳播及其全波場特征，而且該算法穩(wěn)定，無明顯的頻散現(xiàn)象，計算精度高，對以后實際的三維地震勘探技術(shù)會有一定的幫助。

圖4 t=500ms時刻x=500m處yoz平面上的粘彈性波波場快照

圖5 t=500ms時刻z=500m處xoy平面上的粘彈性波波場快照

[1]Stanke F E,Bumidge R.Comparison of Spatial and Ensemble Averaging Methods Applied to Wave Propagation in Finely Lay?ered Media[C]//60th Ann.Internat.Mtg.，Soc.expl.Geophys.Ex? panded Abstracts，1990:1062-1065.

[2]董敏煜.地震勘探[M].東營：中國石油大學出版社，2006.

[3]牟永光，裴正林.三維復雜介質(zhì)地震數(shù)值模擬[M].北京：石油工業(yè)出版社，2005.

[4]陸基孟.地震勘探原理[M].東營：中國石油大學出版社，2009.

[5]Collino F,Tsogka C.Application of the Perfectly Matched Ab?sorbing layer Model to the Linear Elastodynamic Problem in Anisotropic Heterogeneous Media[J].Geophysics，2001，66(1): 294-307.

[6]董良國，馬在田，曹景忠.一階彈性波方程交錯網(wǎng)格高階差分解法穩(wěn)定性研究[J].地球物理學報，2000，43(6):856-864.

[7]Saenger E H，Gold N，Shapiro S A.Modeling the Propagation of Elastic Waves Using a Modified Finite-difference Grid[J]. Wave Motion，2000(31):77-92.

P631

A

1004-5716(2016)12-0026-04

2016-03-04

2016-03-07

陳玉璽（1990-），男（漢族），山東臨沂人，中國海洋大學在讀碩士研究生，研究方向：油氣田與煤田地球物理勘探方法與技術(shù)。