三維粘彈性介質地震波場數值模擬技術研究

陳玉璽，王 璐，章 陽

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

三維粘彈性介質地震波場數值模擬技術研究

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

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

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

非完全彈性；交錯網格高階有限差分；差分算子；數值模擬

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

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

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

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

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

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

2 交錯網格高階有限差分法

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

圖1 三維粘彈性介質交錯網格示意圖

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

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

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

Δt——時間網格間距；

vmax——模型最大速度；

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

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

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

其中，R為理論反射系數，δ為PML厚度[2]。在8個角點區域，d(x)≠0，d(y)≠0，d(z)≠0；在12條棱區域，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個面區域，垂直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。這樣整個模型區域都引入PML吸收邊界，地震波傳播經過模型邊界到達PML吸收層時，能量逐漸減弱，因此不會產生任何邊界反射。

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

其中：

其中：

3 三維復雜粘彈性介質模型

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

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

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

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

4 結論

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

圖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]牟永光，裴正林.三維復雜介質地震數值模擬[M].北京：石油工業出版社，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]董良國，馬在田，曹景忠.一階彈性波方程交錯網格高階差分解法穩定性研究[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-），男（漢族），山東臨沂人，中國海洋大學在讀碩士研究生，研究方向：油氣田與煤田地球物理勘探方法與技術。