任意偶數(shù)階精度交錯網(wǎng)格聲波方程數(shù)值模擬

吳從輝， 李紅星， 張 智

(1.東華理工大學(xué)核工程與地球物理學(xué)院，江西 南昌 330013;2.桂林理工大學(xué)地球科學(xué)學(xué)院，廣西 桂林 541004)

地震波場數(shù)值模擬是地震學(xué)和地震勘探的重要基礎(chǔ)，在地震學(xué)和地震勘探各個工作環(huán)節(jié)中都起到重要作用(裴正林等，2004)。通過正演模擬技術(shù)，可以模擬地震波在各種地質(zhì)結(jié)構(gòu)中的傳播規(guī)律，可為地震資料處理解釋中各種方法的應(yīng)用提供標(biāo)準(zhǔn)的先驗數(shù)據(jù)體，檢測方法的實用性(趙海波等，2011;張華等，2012)。

有限差分方法因其計算效率高、占用計算內(nèi)存小等特點，成為當(dāng)前最常用的數(shù)值模擬方法(張英杰等，2007;呂玉增等，2011)。交錯網(wǎng)格技術(shù)相對于規(guī)則網(wǎng)格在復(fù)雜模型模擬和模擬精度等方面都有較大的優(yōu)勢(李紅星等，2008;李國平等2010)，交錯網(wǎng)格引入的頻散明顯低于規(guī)則網(wǎng)格，因此地震波場數(shù)值模擬中交錯網(wǎng)格有限差分法得到廣泛應(yīng)用。

為得到高精度的正演模擬結(jié)果，本文基于一階應(yīng)力-速度聲波方程組，運用任意偶數(shù)階精度的交錯網(wǎng)格有限差分格式進(jìn)行了各向同性介質(zhì)模型中的正演模擬。

1 波動方程的有限差分格式建立

1.1 一階應(yīng)力-速度聲波方程組

根據(jù)彈性力學(xué)的分析，可得到動力學(xué)方程組和應(yīng)變應(yīng)力方程組。二維情況下一階應(yīng)力-速度聲波方程組(張伯軍，2010)如下:

式中K = ρv2，K 為地質(zhì)體壓強，ρ 表示地質(zhì)體密度，vx表示在x 方向上的速度分量，vz表示在z 方向上的速度分量。

1.2 一階應(yīng)力-速度聲波方程組交錯網(wǎng)格計算實現(xiàn)

一階應(yīng)力-速度聲波方程組各個物性參量在交錯網(wǎng)格中的位置設(shè)置如下:

空間網(wǎng)格中，參數(shù)P 放置在規(guī)則的i，j 節(jié)點上。而參數(shù)vx在x 方向上位于j +1/2 的半節(jié)點上，在z方向上位于規(guī)則的i 節(jié)點上;參量vz在x 方向上位于規(guī)則的j 節(jié)點上，在z 方向上位于i+1/2 半節(jié)點上(圖1)。

時間網(wǎng)格中，參數(shù)P 在規(guī)則的k 節(jié)點上，vx與vz均設(shè)置在k +1/2 的半節(jié)點上(圖2)。

圖1 空間維度上的網(wǎng)格點分布Fig.1 The grid point distribution in space domain

圖2 時間維度上的網(wǎng)格點分布Fig.2 The grid point distribution in time domain

1.3 任意偶數(shù)階交錯網(wǎng)格差分計算格式實現(xiàn)

本文從Taylor 級數(shù)展開式出發(fā)，實現(xiàn)了任意偶數(shù)階精度交錯網(wǎng)格下的差分格式算法(劉洋等，1998;陳可洋，2009)。設(shè)函數(shù)U(x)為存在2L +1 階導(dǎo)數(shù)的函數(shù)，其在空間i+1/2 半節(jié)點處，分別對i 與i +1 點進(jìn)行泰勒展開處理，得到Ui+1和Ui在i+1/2處的泰勒展開式，即可以得到式(2)和式(3)。

將式(2)與式(3)相減得:

類似上述推導(dǎo)過程，可得任意Ui+m- Ui-m+1(m= 1，2，3…L)的推導(dǎo)結(jié)果:

其中m = 1，2，3，…L。由上式可得式(6)，用于計算

解方程式得:

式中am便是任意偶數(shù)階精度交錯網(wǎng)格下的一階導(dǎo)數(shù)差分系數(shù)，于是:

1.4 一階應(yīng)力-速度聲波方程組任意偶數(shù)階精交錯網(wǎng)格差分算法實現(xiàn)

由式(1)給出的二維情況下一階應(yīng)力- 速度方程組，可知需對vx，vz，P 三個參量進(jìn)行計算。采用時間二階精度的差分方式對三個參量的時間導(dǎo)數(shù)進(jìn)行計算:

由式(1)結(jié)合式(9)及式(10)給出的一階導(dǎo)數(shù)差分計算公式，便可實現(xiàn)時間二階空間任意偶數(shù)階精度交錯網(wǎng)格下的聲波方程差分格式(11)。

2 數(shù)值模擬中的幾個重要問題

2.1 震源函數(shù)

震源函數(shù)有多種取法，如高斯子波、雷克子波等。在波場數(shù)值模擬中，震源函數(shù)的選取對模擬結(jié)果有重要的影響(張海燕等，2007;董良國等，2004)。本文采用雷克子波作為震源函數(shù)。

2.2 穩(wěn)定性條件

對文中一階應(yīng)力-速度聲波方程組任意偶數(shù)階精度交錯網(wǎng)格有限差分方程，通過平面諧潑的分析(裴正林等，2003)，得出其數(shù)值解的穩(wěn)定性條件:

其中，Δt 為時間域網(wǎng)格步長，Δx、Δz 表示空間域網(wǎng)格步長，vp表示介質(zhì)的縱波速度，am為相應(yīng)的交錯網(wǎng)格差分算子系數(shù)。由式(12)得出交錯網(wǎng)格高階有限差分的穩(wěn)定性隨著差分精度的提高而不斷增高(Cao et al.，1997)。

2.3 衰減邊界條件

基于在模擬區(qū)外圍附加吸收層的思想，Cerjan等(1985)提出了衰減邊界條件。在模擬區(qū)域內(nèi)，波按照正常的波動傳播;當(dāng)波到達(dá)外圍衰減區(qū)時，波場乘以一個衰減因子G，使其按指數(shù)級逐漸衰減。其中:式中N 為衰減區(qū)總網(wǎng)格數(shù)，i 為衰減區(qū)的網(wǎng)格編號(1 ≤i ≤N)，a 為衰減系數(shù)。圖3 為衰減邊界示意圖。

圖3 衰減邊界示意圖Fig.3 Pictorial view of the damping boundary condition

李信富等(2007)研究表明，衰減因子參數(shù)的選擇應(yīng)采取折中方法，即在取某一較小衰減系數(shù)的前提下，設(shè)置較少的衰減帶網(wǎng)格數(shù)使邊界的能量反射最小。經(jīng)對模擬速度和精度的綜合考慮，本文a 值選用0.15，N 值為50。那么在衰減區(qū)波場計算公式為:

3 精度選取測試及邊界效果對比

3.1 差分精度選取測試

理論上說，差分算子的階數(shù)越高，波場模擬結(jié)果也越精確。但階數(shù)過高，會造成一定的計算困難(左瑩，2007)。那么在效果與效率平衡的基礎(chǔ)上，選取合適的差分精度是很有必要的。

為確定合適的差分精度，本文進(jìn)行了差分精度選取測試。測試模型為均勻各向同性介質(zhì)模型，模型介質(zhì)參數(shù):v =1 000 m/s，ρ =1 g/m3。模型大小為1 000 ×1 000 m，x 和z 方向上均有201個網(wǎng)格點，空間網(wǎng)格步長dx，dz 為5 m，dt 為0.4 ms。震源設(shè)置在模型中心點處，震源函數(shù)頻率為25 Hz。采用時間二階精度，空間不同偶數(shù)階精度進(jìn)行模擬試算。得到空間不同偶數(shù)階精度下的波場模擬快照(圖4)。

通過對比分析，易得出隨著空間差分階數(shù)的提高，模擬精度也逐漸增高?？臻g二階精度時，存在嚴(yán)重的頻散現(xiàn)象;空間四階精度時，頻散現(xiàn)象明顯減弱?？臻g八階精度時，幾乎無頻散現(xiàn)象出現(xiàn);階數(shù)再增加，模擬精度變化不明顯。

因此，為得到高精確的正演模擬結(jié)果，空間差分精度不得低于八階。本文的波場模擬都是采用空間八階精度進(jìn)行的。

圖4 空間不同偶數(shù)階精度下的波場模擬快照Fig.4 Wave field snapshots at diffierent even-order accurate in space

3.2 衰減邊界效果對比

測試模型為均勻各向同性介質(zhì)模型。震源放置在模型中心點處，震源函數(shù)頻率為25 Hz。選取640 ms 時未使用衰減邊界條件和使用衰減邊界條件的波場快照進(jìn)行邊界效果對比(圖5)。

對比可知，未使用衰減邊界條件的波場快照(圖5a)，邊界反射明顯;使用衰減邊界條件的模型波場快照(圖5b)，未出現(xiàn)邊界反射現(xiàn)象。結(jié)果表明，衰減邊界條件可以有效消除邊界反射。

圖5 未使用衰減邊界和使用衰減邊界時波前快照對比圖Fig.5 Wave field snapshots without using the attenuation boundary condition and wave field snapshots with using it

4 模型正演模擬

分別對不同類型的模型進(jìn)行了正演模擬。模擬過程均采用時間二階空間八階差分精度。

4.1 水平層狀介質(zhì)模型正演模擬

如圖6a 所示水平層狀介質(zhì)模型，模型大小為1 000 ×1 000 m，水平分界面在埋深500 m 處。上層介質(zhì)參數(shù):v1=1 000 m/s，ρ1=1 g/m3;下層介質(zhì)參數(shù):v2=1 500 m/s，ρ2=1 g/m。模型剖分為201 ×201個離散節(jié)點，空間網(wǎng)格步長為5 m，時間網(wǎng)格步長為0.4 ms。

震源位于水平方向500 m，埋深10 m 處，震源函數(shù)頻率為25 Hz。對模型進(jìn)行正演模擬，得到不同時刻的波場快照。

圖6b 和圖6c 分別為水平層狀介質(zhì)模型560 ms 和600 ms 時刻的波場快照。波場快照中，界面反射波和透射波清晰可見，邊界反射被有效吸收，未出現(xiàn)頻散現(xiàn)象。圖6d 為水平層狀介質(zhì)模型的共炮點記錄，直達(dá)波及水平分界面z =500 m 處的反射波同相軸清晰。

4.2 傾斜介質(zhì)模型正演模擬

介質(zhì)模型依舊采用201 ×201 的網(wǎng)格剖分，空間網(wǎng)格步長為5 m，時間網(wǎng)格步長0.4 ms，第一層介質(zhì)模型參數(shù):v1=1 000 m/s，ρ1=1 g/m3;第二層介質(zhì)模型的參數(shù)為:v2=1 500 m/s，ρ2=1 g/m3;第三層的介質(zhì)模型參數(shù)為:v3=2 000 m/s，ρ3=1 g/m3。

圖6 水平層狀介質(zhì)模型、模型波場快照和共炮點記錄Fig.6 The model of horizontal stratified medium、Wave field snapshots and Synthetic seismogram of horizontal stratified medium

圖7 傾斜介質(zhì)模型、模型波場快照和共炮點記錄Fig.7 The model of layered medium with tilted interface，Wave field snapshots and Synthetic seismogram of layered medium with tilted interface

傾斜介質(zhì)模型中傾斜分界面從左邊界埋深300 m 處延伸到右邊界埋深550 m 處(圖7a);水平分界面在埋深z =550 m 處。震源位于水平方向600 m，埋深10 m 處，震源函數(shù)頻率為25 Hz。同樣采用空間八階時間二階差分精度進(jìn)行波場模擬。

圖7b 和圖7c 分別為傾斜介質(zhì)模型420 ms 和600 ms 時刻的波場快照，可清晰反映出波在模型中的傳播規(guī)律。邊界反射被有效吸收，但存在輕微頻散現(xiàn)象。圖7d 為傾斜介質(zhì)模型的共炮點記錄，圖中直達(dá)波同相軸和界面反射波同相軸清晰明顯。邊界反射被有效吸收，整體效果良好。

4.3 斷層介質(zhì)模型正演模擬

如圖8a 所示斷層介質(zhì)模型，斷層分界面中，垂直段在x=600 m 處，從埋深380 m 延伸到520 m。斷層水平段一部分在埋深380 m 處，水平坐標(biāo)x=0至x=600 之間;水平段另一部分在埋深520 m 處，水平坐標(biāo)x=600 至x =1 000 m 之間。第二、三層介質(zhì)的水平分界面在z =600 m 處。第一層模型的介質(zhì)參數(shù)為:v1=1 000 m/s，ρ1=1 g/m3;第二層的模型介質(zhì)參數(shù)為:v2=1 500 m/s，ρ2=1 g/m3;第三層的模型介質(zhì)參數(shù)為:v3=2 000 m/s，ρ3=1 g/m3。

震源位于水平方向500 m，埋深10 m 處，震源函數(shù)頻率為25 Hz。采用同樣的差分精度對斷層介質(zhì)模型進(jìn)行模擬試算。

圖8b 和圖8c 分別為420 ms 和600 ms 時刻的模型波場快照，充分反映出斷層介質(zhì)模型中波的傳播規(guī)律。圖中可清晰分辨出界面反射波和透射波，以及斷層的斷點上產(chǎn)生的繞射波;邊界反射吸收效果良好，未出現(xiàn)頻散現(xiàn)象。圖8d 為斷層介質(zhì)模型的共炮點記錄。圖中直達(dá)波、斷層界面的反射波及水平界面反射波同相軸清晰。斷層上下水平段界面對應(yīng)的反射波同相軸末端，都伴隨著明顯的繞射波同相軸，與實際情況相符合。

圖8 斷層介質(zhì)模型、模型波場快照和共炮點記錄Fig.8 The fault modeling，Wave field snapshots and Synthetic seismogram in the fault modeling

5 結(jié)論

高偶數(shù)階交錯網(wǎng)格有限差分?jǐn)?shù)值模擬方法可以清晰模擬出波在地下介質(zhì)中的傳播規(guī)律。本文方法模擬精度高，計算效率快，能有效壓制頻散噪音，邊界反射吸收效果良好，充分證實了該數(shù)值模擬方法的準(zhǔn)確性和可行性，具有較高的應(yīng)用價值?？梢酝茝V應(yīng)用于三維聲波高精度正演模擬中，幫助我們深入學(xué)習(xí)和研究地震波的傳播規(guī)律，也為地震勘探工作提供了實用的參考價值。

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

董良國，李培明.2004.地震波傳播數(shù)值模擬中的頻散問題［J］.天然氣工業(yè)，24(6):53-56.

李國平，程利敏.2010.非均勻介質(zhì)交錯網(wǎng)格高階差分地震波數(shù)值模擬［J］.石油天然氣學(xué)報，32(4):220-225.

李紅星，劉財，吳志成，等.2008. 基于BISQ 機制的雙相孔隙介質(zhì)三維數(shù)值模擬研究［J］.石油天然氣學(xué)報，30(2):77-80.

李信富，李小凡，張美根.2007. 偽譜法彈性波場數(shù)值模擬中的邊界條件［J］.地球物理學(xué)進(jìn)展，22(5):1375-1379.

劉洋，李承楚，牟永光.1998.任意偶數(shù)階精度數(shù)值模擬算法［J］. 石油地球物理勘探，33(1):1-10.

呂玉增，熊彬，薛霆虓.2011. 地球物理數(shù)據(jù)處理基礎(chǔ)［M］. 北京:地質(zhì)出版社:90-106.

裴正林，牟永光.2003.非均勻介質(zhì)地震波傳播交錯網(wǎng)格高階有限差分法模擬［J］.石油大學(xué)學(xué)報，27(6):17-21.

裴正林，牟永光.2004.地震波傳播數(shù)值模擬［J］. 地球物理學(xué)進(jìn)展，19(4):933-941.

張伯軍.2010.彈性動力學(xué)簡明教程［M］.北京:科學(xué)出版社:35-64.

張海燕，李慶忠.2007.幾種常用解析子波的特征分析［J］.石油地球物理勘探，42(6):651-657.

張華，龔育齡.2012.花崗巖裂隙儲層各向異性高階有限差分?jǐn)?shù)值模擬［J］.東華理工大學(xué)學(xué)報:自然科學(xué)版，35(4):416-421.

張英杰，樣長春，吳萍.2007.地震波有限差分模擬綜述［J］. 地球物理學(xué)進(jìn)展，22(2):487-491.

趙海波，陳百軍，李奎周，等.2011.黏彈性介質(zhì)VSP 記錄模擬及在估計Q 值研究中應(yīng)用［J］.地球物理學(xué)報，54(2):329-335.

左瑩.2007. 基于高階交錯網(wǎng)格的有限差分地震波場數(shù)值模擬［D］.西安:長安大學(xué).

Cao S H，Stewart G.1997.Attenuating boundary conditions for numerical modeling of acoustic wave propagation［J］. Geophysics，63(1):231-243.

Cerjan C，Kosloff R，Reshef M.1985. A nonreflecting boundary condition for discrete acoustic and elastic wave equations［J］. Geophysics，50(4):705-708.