三維快速高精度地震波正演數值模擬方法及其應用

陳可洋

(中國石油大慶油田有限責任公司勘探開發研究院)

三維快速高精度地震波正演數值模擬方法及其應用

陳可洋

(中國石油大慶油田有限責任公司勘探開發研究院)

如何有效提高三維地震波正演數值模擬精度和計算效率一直是勘探地球物理學研究的重要問題。為了克服常規中心有限差分法較難快速提高差分精度的缺陷和一階雙曲型波動方程內存占用多、計算量大、引入變量較多的困難，采用高階交錯網格有限差分法直接求解三維地震波動方程，推導的高階差分格式計算形式簡單，可以推廣于求解任意偶數階時空導數，同時給出其穩定性條件。在人工邊界處，對比了鑲邊法和常規旁軸近似法兩種吸收邊界條件。從三維似French模型的正演結果看出，采用的高階交錯網格差分算法在快速有效地提高數值模擬精度的同時，大大提高了計算效率，同時結合鑲邊法吸收邊界條件還可有效壓制邊界反射，提高整個計算域內波場的信噪比。圖3參5

三維地震波動方程 高階交錯網格有限差分法 正演數值模擬 鑲邊法吸收邊界

0 引言

針對當前高精度地震勘探的要求，地震勘探方法必須考慮地下三維空間內非均勻介質對地震資料采集的影響。三維地震波正演數值模擬方法因此成為準確認識地震波場傳播規律(保留幾何學、運動學、動力學特征)、指導地震觀測系統的優化設計和檢驗地震資料處理與解釋方法準確性的一種重要手段。只有準確地研究復雜的地質構造和油氣儲集體所對應的地震波場特征，才能有效地進行構造和儲層的識別與劃分。傳統的基于褶積模型的正演方法僅考慮了縱向上介質的變化，無法完整地描述三維空間的局部構造或非均勻性介質變化產生的復雜波場響應［1，2］。

目前，國內外對三維地震波的數值模擬方法進行了大量研究，逐步將二維方法推廣應用于三維情況，主要包括單程波正演方法(如隱式有限差分法、Fourier法、傅里葉有限差分法、顯式短算子方法等)和雙程波正演方法(顯式有限差分法、隱式有限差分法、有限元法、精細積分法、偽譜法、Hartley變換法等)，其中單程波正演方法是在頻率－空間域進行交互處理，在每一步波場遞推過程中，均需引入正反傅立葉變換，因而計算量較大。而雙程波正演方法是在時空域進行計算的，因此其計算量較小，其中使用最多的是有限差分法，常規中心有限差分法較難快速提高有限差分精度，如果將標量地震波動方程轉化為一階雙曲型方程來計算，則需要引入幾個輔助變量，這將增加計算量和計算的復雜度。另外，對于大工區的三維地震波正演而言，其計算機內存的占用量是很龐大的，計算效率也必然較低。為了克服上述這些難題，本文在總結前人研究成果的基礎之上，提出了將高階交錯網格有限差分法直接引入到求解三維標量地震波動方程的新思路，并結合鑲邊法吸收邊界條件，以期快速提高局部地震波場的數值模擬精度和計算效率。

1 基本原理

1.1 三維高階交錯網格有限差分算子的構造

一般情況下，三維地震波動方程的形式如下:

式中:

u—質點的振動位移;

v—介質速度;

t、x、y、z—分代表時間和三個空間方向。

以x方向為例，定義三維二階導數的任意高階交錯網格有限差分格式如下:

式(2)中，Lx=?/?x為x方向的空間微分算子，代表后向差分算子代表前向差分算子，二階導數的差分算子就是將前向差分算子和后向差分算子組合得到，另外，(i，j，k，n)中的每一個量依次代表(x，y，z，t)的每一個方向的離散位置變量，am為高階交錯網格有限差分系數，N為差分階數(正整數)。

將式(1)按照式(2)進行差分離散，得到的三維地震波動方程的時間2階、空間4N－2階交錯網格有限差分法計算公式如下:

分析公式(3)可知，本文方法適用于非均勻網格的三維地震波正演數值模擬。此外可以看出，任意空間方向二階導數的高階交錯網格有限差分格式具有差分規律，在相同差分階數N情況下，常規中心網格有限差分法的差分精度為2N，而本文方法的差分精度為4N－2。由此可以看出，本文方法的差分精度與差分階數N是近似4倍的關系，且為常規方法差分精度的2倍。另外，求解三維地震波動方程只需要引入三個不同時間層、相同的遞推變量(即同一變量在三個不同時間層的不同表示)，而采用相同情況下的一階雙曲型方程則至少需要引入四個交錯時間層、不相同的遞推變量。因此，本文方法可以大大提高三維地震波正演數值模擬的計算效率。另外，文中的高階差分方法還可以推廣應用于求解任意偶數階導數，仍以x方向為例，其計算通式為:其中2m為x方向導數的階數。式(4)在高階時間差分近似情況下較為常用(此時通常是將時間的高階導數轉化為空間高階導數來實現)。

1.2 穩定性條件和吸收邊界條件

經推導，式(3)的穩定性條件［3］與相應的一階雙曲型情況相一致，其表達式如下:

其中，S=Max{v2/Δx2，v2/Δy2，v2/Δz2}，在高階時間差分近似條件下，Δt的上限值可以適當放寬。

為了能夠削弱或消除計算邊界處的反射波，同時保證邊界計算過程的穩定，在人工截斷邊界處，采用了常規旁軸近似吸收邊界條件［4］(不需要進行外側鑲邊，僅依賴于邊界附近節點處不同時間層的值，并采用二階近似的單程地震波動方程來進行邊界點值的預測，其計算效率最高，但受到邊界吸收角度的限制)和一定厚度的外側鑲邊法吸收邊界條件［5］(即在3D模型的邊界外側，鑲上20個網格節點的阻尼條帶，使得邊界反射波在該條帶內多次吸收衰減，其精度較高，但效率偏低)兩種方法，并對比這兩種吸收邊界條件對三維空間有效波場響應的影響。

2 應用實例

采用的速度模型類似于French模型，其中包含一水平層界面、一個傾斜斷層、一個隆起構造、一個凹陷構造(圖1，數值越大代表層界面離地表越深)，模型總大小為500m×500m×500m，三個方向空間網格大小均為5m，震源位于模型地表中央位置，其主頻為60Hz，模型中上層介質速度為2000m/s，下層介質速度為2500m/s，時間步長為0.5ms，滿足穩定性條件式(5)，時空差分精度為(Δt2+Δx10)(此時N=3)。檢波器布置于地表下方50m深度處，自激自收方式合成的水平疊加剖面(這里考慮了速度差異界面反射系數的大小)，在三維模型的邊界處分別采用常規旁軸近似吸收邊界條件和一定厚度的外側鑲邊法吸收邊界條件。

圖2分別為0.175s時刻，x、y、z三方向的中間位置且平行于yoz、xoz、xoy三個平面的波場快照切片，圖中的①和②與圖1中波場快照切片位置(虛線)相對應，邊界處采用了外側鑲邊法吸收邊界條件。對比速度模型和波場快照可以比較容易地識別出直達波、傾斜界面的反射波、水平界面的反射波以及隆起構造的反射波這四種波，可以看出，地震波是在三維空間中進行傳播的，僅考慮二維空間的傳播問題是不準確的。另外，在地震波到達模型邊界時無邊界反射波形成，在波場傳播過程中無任何頻散現象，這表明本文算法精度較高。

圖3(a)和圖3(b)分別為采用鑲邊法吸收邊界條件情況下的縱觀測系統和非縱觀測系統接收到的單炮模擬記錄。分析圖3(a)和圖3(b)可知，邊界反射波的能量較弱，而有效波能量強，同相軸清晰，結合速度模型可以較準確地識別出各種有效波的來源，從而可以進行三維空間地震波場的傳播規律研究。圖3(c)、圖3(d)、圖3(e)、圖3(f)分別為采用常規旁軸近似吸收邊界條件和鑲邊法吸收邊界條件情況下接收到的兩條正交的自激自收剖面(類似于水平疊加剖面，剖面位置如圖1中的虛線位置所示，①代表x方向剖面位置，②代表y方向剖面位置)，比較x方向剖面(圖3(c)和圖3(d))和y方向剖面(圖3(e)和圖3(f))可知:采用旁軸近似吸收邊界條件時，地震波場不清晰，比較雜亂，信噪比很低;而采用鑲邊法吸收邊界條件后，有效波特征明顯，同相軸易追蹤。再結合三維速度模型和三維地震波的波場快照就可以很容易地識別出各種復雜地震波場的成因問題。

圖1 三維似French速度型

圖2 0.175s時刻，三維任意方向波場快照切片

圖3 三維地震波數值模擬記錄

3 結論

本文提出了用高階交錯網格有限差分法直接求取三維地震波動方程的新思路，并詳細推導得到了三維地震波動方程的高精度離散方程，給出了計算所需的穩定性條件和吸收邊界條件。從計算量、計算效率、計算復雜度上對比分析了本文方法和一階雙曲型方法，得出本文方法在三維快速高精度實現正演數值模擬方面的優點為:①計算速度快;②占用內存小;③計算格式簡單有規律，且計算復雜度低;④計算量小。基于這些優點，本文方法可以用來快速模擬野外三維地震資料的采集過程以及觀測系統的優化設計，同時還可以推廣應用于三維二階各向異性介質彈性波的傳播數值模擬問題。

從分別采用外側鑲邊法吸收邊界條件和常規旁軸近似吸收邊界條件的計算結果可以看出，前者邊界吸收效果好，三維模擬記錄中無任何頻散現象，能夠較清晰地識別出各種構造所形成的復雜反射波，而后者邊界反射波與有效波相互疊加，造成了信噪比降低，有效波場模糊、難識別。因此，采用本文方法并結合鑲邊法吸收邊界條件就可以快速有效地模擬三維地震波的傳播規律。

1 陳可洋，劉洪林，楊微，等.隨機介質模型的改進方法及應用［J］.大慶石油地質與開發.2008，27(5):124－126，131.

2 陳可洋.三維隨機建模方法及其波場模擬分析［J］.勘探地球物理進展，2009，32(5):315－320.

3 陳可洋.標量聲波波動方程高階交錯網格有限差分法［J］. 中國海上油氣，2009，21(4):232 －236.

4 楊微，陳可洋.加權吸收邊界條件的優化設計［J］.石油物探，2009，48(3):244 －246，251.

5 陳可洋.邊界吸收中鑲邊法的評價［J］.中國科學院研究生院學報，2010，27(2):170 －175.

3D FAST AND HIGH－RESOLUTION SEISMIC－WAVE FORWARD NUMERICALSIMULATION AND ITS APPLICATION

CHEN Keyang(Research Institute of Exploration and Development，PetroChina Daqing Oilfield Company)．

How to effectively increase both accuracy and calculation efficiency of 3D seismic－wave forward numerical simulation is an important problem in geophysical prospect．But for conventional central finite－ difference method，it is difficult to fast improve difference accuracy;and for one－stage dual－curve wave equation，there are some defects of occupying much memory，large amount of calculation and introducing much variable．In this study，a method of higher－order staggered－grid finite difference is adopted to directly solve a 3D seismic wave equation and there are some advantages:(1)simple calculation form;(2)it may also be applied to solving a random even－order space－time derivative;(3)it can provide with some stable conditions．Moreover，edging method is correlated to conventional paraxial approximation to adsorb in boundary condition．It is shown from the forward result of 3D quasi－French model that:(1)the higher－order staggered－grid finite difference method can not only fast and effectively improve simulation accuracy but also increase calculation efficiency;and(2)combined with edging method，the higher－order staggered－grid finite difference method can effectively impose boundary reflection and improve signal－to－noise ratio of wave field within whole calculation domain．

3D seismic wave equation，higher－order staggered－grid finite difference method，forward numerical simulation，edging method

陳可洋，男，1983年出生，助理工程師;2009年獲大慶石油學院地球探測與信息技術專業碩士學位，現主要從事高精度彈性波正演數值模擬及逆時偏移成像方法研究。地址:(163712)黑龍江省大慶市讓胡路區大慶油田勘探開發研究院地震處理二室。電話:(0459)5508524，13504595794。E － mail:keyangchen@163.com

(修改回稿日期 2010－12－30 編輯 陳 玲)NATURALGAS EXPLORATION＆DEVELOPMENT．v．34，no．3 ，pp．12－15，7/25/2011