陳可洋

(中國(guó)石油大慶油田有限責(zé)任公司勘探開發(fā)研究院)

三維快速高精度地震波正演數(shù)值模擬方法及其應(yīng)用

陳可洋

(中國(guó)石油大慶油田有限責(zé)任公司勘探開發(fā)研究院)

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

三維地震波動(dòng)方程 高階交錯(cuò)網(wǎng)格有限差分法 正演數(shù)值模擬 鑲邊法吸收邊界

0 引言

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

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

1 基本原理

1.1 三維高階交錯(cuò)網(wǎng)格有限差分算子的構(gòu)造

一般情況下，三維地震波動(dòng)方程的形式如下:

式中:

u—質(zhì)點(diǎn)的振動(dòng)位移;

v—介質(zhì)速度;

t、x、y、z—分代表時(shí)間和三個(gè)空間方向。

以x方向?yàn)槔?，定義三維二階導(dǎo)數(shù)的任意高階交錯(cuò)網(wǎng)格有限差分格式如下:

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

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

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

1.2 穩(wěn)定性條件和吸收邊界條件

經(jīng)推導(dǎo)，式(3)的穩(wěn)定性條件［3］與相應(yīng)的一階雙曲型情況相一致，其表達(dá)式如下:

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

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

2 應(yīng)用實(shí)例

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

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

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

圖1 三維似French速度型

圖2 0.175s時(shí)刻，三維任意方向波場(chǎng)快照切片

圖3 三維地震波數(shù)值模擬記錄

3 結(jié)論

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

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

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

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

3 陳可洋.標(biāo)量聲波波動(dòng)方程高階交錯(cuò)網(wǎng)格有限差分法［J］. 中國(guó)海上油氣，2009，21(4):232 －236.

4 楊微，陳可洋.加權(quán)吸收邊界條件的優(yōu)化設(shè)計(jì)［J］.石油物探，2009，48(3):244 －246，251.

5 陳可洋.邊界吸收中鑲邊法的評(píng)價(jià)［J］.中國(guó)科學(xué)院研究生院學(xué)報(bào)，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年獲大慶石油學(xué)院地球探測(cè)與信息技術(shù)專業(yè)碩士學(xué)位，現(xiàn)主要從事高精度彈性波正演數(shù)值模擬及逆時(shí)偏移成像方法研究。地址:(163712)黑龍江省大慶市讓胡路區(qū)大慶油田勘探開發(fā)研究院地震處理二室。電話:(0459)5508524，13504595794。E － mail:keyangchen@163.com

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