二維聲波方程交錯(cuò)網(wǎng)格有限差分?jǐn)?shù)值模擬研究

楊曉柳，竺 俊，姚 銘，桂志先

(1. 長(zhǎng)江大學(xué) 油氣資源與勘探技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430100； 2. 長(zhǎng)江大學(xué) 地球物理與石油資源學(xué)院，湖北 武漢 430100； 3. 中石化西北油田分公司采油一廠，新疆 輪臺(tái) 841600)

0 前 言

在地震勘探領(lǐng)域中，人們常常利用地震波場(chǎng)數(shù)值模擬技術(shù)研究地震波的傳播機(jī)理。地震數(shù)值模擬就是提前給定地下介質(zhì)模型和相應(yīng)物理參數(shù)，模擬地震波在地下介質(zhì)中的傳播過(guò)程，并對(duì)地面或地下各個(gè)觀測(cè)點(diǎn)觀測(cè)到的數(shù)值地震記錄進(jìn)行計(jì)算分析的一種地震模擬方法。地震數(shù)值模擬方法在地震勘探和天然地震領(lǐng)域得到了普遍地推廣使用，在地震數(shù)據(jù)采集、處理和解釋的各個(gè)環(huán)節(jié)都發(fā)揮著重要作用。

地震波場(chǎng)數(shù)值模擬方法主要包括有限差分法、有限元法、偽譜法等[1]，其中有限差分法是利用 Taylor 展開式來(lái)近似代替波動(dòng)方程中的微分，使微分形式的波動(dòng)方程變成差分形式，來(lái)用于求解波動(dòng)方程。他具有占有內(nèi)存小、計(jì)算速度快，編程相對(duì)簡(jiǎn)單等優(yōu)點(diǎn)[1]，是最常用的地震波場(chǎng)數(shù)值模擬方法之一。本文通過(guò)建立地質(zhì)模型，利用Matlab語(yǔ)言進(jìn)行編程，完成了聲波方程交錯(cuò)網(wǎng)格有限差分算法的數(shù)值模擬。通過(guò)模擬聲波在不同介質(zhì)中的傳播過(guò)程，對(duì)比分析不同時(shí)刻、不同介質(zhì)的波場(chǎng)快照，深入理解地震波在不同介質(zhì)中的傳播規(guī)律。有限差分算法在模擬精度和模擬計(jì)算速率方面都展現(xiàn)出一定優(yōu)勢(shì)。綜上所述，地震波場(chǎng)數(shù)值模擬對(duì)于人們深入研究地震波的傳播規(guī)律，準(zhǔn)確解釋實(shí)際勘探中的地震資料等均具有重要的理論和實(shí)際意義。

1 聲波方程及其交錯(cuò)網(wǎng)格有限差分格式

在實(shí)際的地震勘探過(guò)程中，最常用到的數(shù)值資料是縱波資料。因此討論利用聲波方程來(lái)進(jìn)行地震波的模擬。有限差分法是基于差分原理的一種數(shù)值計(jì)算方法。其基本思想是：將求解區(qū)域網(wǎng)格化，應(yīng)用差分原理，用差商來(lái)近似的代替微商，將求解連續(xù)函數(shù)的問(wèn)題轉(zhuǎn)化為求解網(wǎng)格節(jié)點(diǎn)上的差分方程組的問(wèn)題，得到數(shù)值解。在波動(dòng)方程網(wǎng)格離散化的過(guò)程中，可以利用交錯(cuò)網(wǎng)格的差分形式。

1.1 均勻各向同性介質(zhì)二維聲波方程

均勻各向同性介質(zhì)二維聲波方程可表示為：

(1)

其中，u(x,y,z,t)為地表記錄的壓力波場(chǎng)或位移波場(chǎng)，v(x,z)為聲波速度場(chǎng)。

1.2 非均勻介質(zhì)中二維聲波方程

非均勻各向同性介質(zhì)中二維聲波方程的一階應(yīng)力—速度方程形式如下：

(2)

其中，vx，vz是質(zhì)點(diǎn)速度，u是法線應(yīng)力，ρ是密度，vp是縱波速度。

1.3 交錯(cuò)網(wǎng)格有限差分

在對(duì)介質(zhì)模型進(jìn)行離散化處理的過(guò)程中，網(wǎng)格是一種常用手段。對(duì)波動(dòng)方程進(jìn)行網(wǎng)格離散，可以利用交錯(cuò)網(wǎng)格的差分形式。交錯(cuò)網(wǎng)格就是把速度和應(yīng)力分配到兩套不同的網(wǎng)格中，這樣可以使速度、應(yīng)力得到很好的耦合[2]。利用交錯(cuò)網(wǎng)格有限差分法對(duì)一階速度—應(yīng)力波動(dòng)方程進(jìn)行求解時(shí)，應(yīng)力、速度等分量在模型交錯(cuò)網(wǎng)格節(jié)點(diǎn)中的位置分布如圖1所示。

圖1 交錯(cuò)網(wǎng)格剖分

其中三角形為τxx，τzz，四邊形為τxz，圓形位于個(gè)立方體中心為vx，六邊形為τz。

利用Taylor級(jí)數(shù)展開法，把波動(dòng)方程中的導(dǎo)數(shù)用網(wǎng)格節(jié)點(diǎn)上的函數(shù)值的差商近似代替，進(jìn)行離散；推導(dǎo)出了時(shí)間域2階、空間域2N階精度的交錯(cuò)網(wǎng)格離散差分格式。

速度水平分量：

速度垂直分量：

應(yīng)力Pxx：

應(yīng)力Qzz：

應(yīng)力Sxz：

2 有限差分?jǐn)?shù)值模擬中的關(guān)鍵技術(shù)及解決方案

2.1 模擬震源

正演模擬中有一個(gè)重要的部分就是加載震源。需要加載震源，震源的加載方式以及類型等因素，都能影響正演模擬的結(jié)果。本文選擇的是雷克子波作為加載的震源，來(lái)模擬激發(fā)地震波。雷克子波是零相位的，零相位子波可以達(dá)到分辨率的極限[1]。

雷克子波表達(dá)式(時(shí)間域)為：

(3)

通過(guò)傅立葉變換為頻率域：

(4)

其中，初始時(shí)刻是t0，主頻是fM。我選用頻率為40 Hz，然后采用時(shí)間域的雷克子波的公式做成圖，如圖2所示，該圖為本文在進(jìn)行地震波場(chǎng)數(shù)值模擬時(shí)才采用的雷克子波波形。

2.2 邊界條件

在進(jìn)行數(shù)值模擬的過(guò)程中，計(jì)算機(jī)的內(nèi)存是有限的，無(wú)法做到模擬一個(gè)無(wú)限大的介質(zhì)區(qū)域，都會(huì)存在一定的范圍限制，即存在一定的邊界，這就是我們所說(shuō)的人為邊界。而在實(shí)際情況下，真正的地下介質(zhì)類似于半無(wú)限大的空間。因此，在有限的計(jì)算區(qū)域內(nèi)，地震波在傳播過(guò)程中若遇到人工邊界就會(huì)產(chǎn)生邊界反射，這些來(lái)自邊界上的邊界反射波會(huì)對(duì)計(jì)算區(qū)域的波場(chǎng)值造成干擾，產(chǎn)生邊界效應(yīng)，嚴(yán)重降低地震波數(shù)值模擬結(jié)果的準(zhǔn)確性。因此在進(jìn)行波場(chǎng)模擬的時(shí)候。引入恰當(dāng)?shù)倪吔鐥l件是十分重要的[3]。

圖2 主頻為40 Hz的雷克子波的特征

為了盡量降低人為邊界所產(chǎn)生的影響，我們可以引入恰當(dāng)?shù)倪吔鐥l件來(lái)更加真實(shí)的模擬無(wú)限的實(shí)際空間，盡最大可能去削弱邊界反射，改善正演模擬的效果。

許多學(xué)者進(jìn)行了相關(guān)方面的研究，提出多種構(gòu)造邊界的方法。Clayton和Engquist[4]在1977年提出了一組包含不同階數(shù)的吸收邊界條件(簡(jiǎn)稱 CE)。Reynolds[5]在1978年通過(guò)分解波動(dòng)方程提出了透明邊界條件(簡(jiǎn)稱TBC)。Berenger[6]在1994年針對(duì)電磁波的傳播，提出了PML(完全匹配層)吸收邊界條件，并從理論上證明了利用PML邊界條件可以吸收任一頻率、任一方向的電磁波。Frank，Hastings 等[7]在1995年將PML的算法思想成功地應(yīng)用于對(duì)彈性波的模擬中。

1)透明邊界條件

(5)

2)Clayton-Engquist 邊界條件

Clayton和Engquist 提出的一組包含不同精度的吸收邊界條件，其中如下的二階近似邊界條件應(yīng)用最為普遍[4]。

(6)

為了測(cè)試邊界條件的作用，本文在此利用一個(gè)均勻介質(zhì)模型，模型的區(qū)域大小為1 000 m×1 000 m，震源位于(500 m，500 m)處，計(jì)算網(wǎng)格大小為5 m×5 m，時(shí)間采樣率為0.5 ms，總采樣時(shí)間為0.2 s，速度采用3 500 m/s，子波頻率為40 Hz。

采用交錯(cuò)網(wǎng)格有限差分算子對(duì)擁有邊界條件跟無(wú)邊界條件的情況分別進(jìn)行模擬，得到t=200 ms時(shí)的波場(chǎng)快照，如圖3所示。

圖3 均勻介質(zhì)模型下有無(wú)邊界條件在t=200 ms的波場(chǎng)快照

由圖3可知，使用了邊界條件的波場(chǎng)快照在邊界處邊界反射波得到明顯的削弱，而沒(méi)有進(jìn)行邊界處理的邊界反射波很明顯。通過(guò)引入恰當(dāng)?shù)倪吔鐥l件，可以有效地削弱邊界反射。

2.3 穩(wěn)定性

由于波動(dòng)方程的有限差分格式一般都是按時(shí)間遞推的，因此前一時(shí)間波函數(shù)值的舍入誤差必然會(huì)對(duì)后一時(shí)間的波場(chǎng)產(chǎn)生影響。這就需要對(duì)誤差傳播和積累的情況進(jìn)行了解，使誤差不至于隨時(shí)間的推進(jìn)而快速累積，破壞整個(gè)數(shù)值解，影響數(shù)值模擬的效果，嚴(yán)重時(shí)還可能導(dǎo)致計(jì)算溢出，使得數(shù)值模擬無(wú)法繼續(xù)進(jìn)行。因此穩(wěn)定性問(wèn)題成為有限差分格式無(wú)法忽略的一個(gè)基本問(wèn)題。根據(jù)Lax等價(jià)定理[8]，在保證穩(wěn)定性的同時(shí)，也就保證了差分格式的收斂性。

穩(wěn)定性問(wèn)題是有限差分法正演模擬中不可忽略的問(wèn)題之一。每進(jìn)行一次波場(chǎng)模擬，都要分析穩(wěn)定性的問(wèn)題。一般來(lái)說(shuō)，隨著時(shí)間的推移，如果誤差不增加，那么這個(gè)差分格式可以認(rèn)為是符合穩(wěn)定性條件的。

針對(duì)均勻介質(zhì)，顯式的穩(wěn)定性條件可以由譜分析得到：

(7)

其中，vz是縱波速度，可見(jiàn)上式與橫波速度vx無(wú)關(guān)。令Δx=Δz，則穩(wěn)定性條件為：

(8)

常規(guī)差分方程其穩(wěn)定性條件為：

(9)

因?yàn)関x

(10)

2.4 數(shù)值頻散問(wèn)題

數(shù)值頻散是有限差分法最大的缺點(diǎn)。在用有限差分進(jìn)行數(shù)值模擬的過(guò)程中，常常會(huì)觀察到，一些同相軸的形狀隨著時(shí)間的推移逐漸變化，有的還逐漸分散，由一根同相軸分裂為許多根同相軸，這種現(xiàn)象是由于波傳播相速度和群速度不一致所導(dǎo)致的，我們稱該現(xiàn)象為數(shù)值頻散現(xiàn)象[9]。數(shù)值頻散在實(shí)質(zhì)上是由于在數(shù)值計(jì)算過(guò)程中離散網(wǎng)格而產(chǎn)生的一種偽波動(dòng)，這種頻散在差分方程求解波動(dòng)方程的過(guò)程中是無(wú)法避免的。

在有限差分正演模擬過(guò)程中，當(dāng)一個(gè)波場(chǎng)內(nèi)采樣點(diǎn)太少時(shí)，用差分逼近微分的過(guò)程中就會(huì)產(chǎn)生較大的誤差，就會(huì)產(chǎn)生數(shù)值頻散[10]。在給定地震子波頻率的情況下，若速度越低，一個(gè)波場(chǎng)內(nèi)的網(wǎng)格點(diǎn)數(shù)就越少，頻散就越嚴(yán)重；在給定空間網(wǎng)格間距的情況下，子波的頻率越高，一個(gè)波場(chǎng)內(nèi)的網(wǎng)格點(diǎn)數(shù)就會(huì)越少，數(shù)值頻散也越嚴(yán)重。因此，在進(jìn)行數(shù)值模擬時(shí)，提高差分階數(shù)，減少網(wǎng)格步長(zhǎng)和降低地震子波頻率是減弱頻散現(xiàn)象最有效的方法。

本文通過(guò)建立一個(gè)均勻介質(zhì)模型，通過(guò)改變網(wǎng)格步長(zhǎng)和子波頻率，利用有限差分算法進(jìn)行數(shù)值模擬，得到波場(chǎng)快照如圖4～5所示。

圖4 網(wǎng)格步長(zhǎng)為5 m(左)與網(wǎng)格步長(zhǎng)為6 m(右)波場(chǎng)快照

圖5 不同主頻波場(chǎng)快照

圖4與圖5表明隨著網(wǎng)格步長(zhǎng)和子波主頻的增加，頻散現(xiàn)象越來(lái)越嚴(yán)重。即可以通過(guò)適當(dāng)減小網(wǎng)格步長(zhǎng)和子波主頻來(lái)壓制頻散現(xiàn)象。

3 有限差分?jǐn)?shù)值模擬實(shí)例

3.1 均勻介質(zhì)模型

模型大小為1 000 m×1 000 m,空間步長(zhǎng)為5 m，速度采用3 500 m/s(圖6)，震源位于模型(500 m，500 m)處，采用雷克子波作為震源，其主頻為40 Hz，滿足彈性波數(shù)值模擬的穩(wěn)定性條件，對(duì)聲波方程利用交錯(cuò)網(wǎng)格有限差分進(jìn)行數(shù)值模擬，總采樣時(shí)間為0.2 s。均勻介質(zhì)模型如圖6所示。

圖6 均勻介質(zhì)模型

模擬結(jié)果如圖7所示。

圖7 均勻介質(zhì)模型在不同時(shí)刻得到的波場(chǎng)快照

圖7可以看出，波場(chǎng)快照中的同相軸是圓形的，說(shuō)明在均勻各向同性介質(zhì)中，點(diǎn)源激發(fā)的波前面是一個(gè)圓，這與理論也是吻合的。并且隨著時(shí)間的增大，波前面的面積逐漸增大，說(shuō)明地震波從震源中心向外傳播。

3.2 雙層介質(zhì)模型

模型大小為2 000 m×2 000 m，空間步長(zhǎng)為5 m，上層介質(zhì)速度為35 000 m/s，下層介質(zhì)速度為4 500 m/s(圖8)，震源位于模型(1 000 m，500 m)處，采用雷克子波作為震源，其主頻為40 Hz。對(duì)聲波方程利用交錯(cuò)網(wǎng)格有限差分進(jìn)行數(shù)值模擬，總采樣時(shí)間為0.2 s。雙層介質(zhì)模型如圖8所示。

圖8 雙層介質(zhì)模型

模擬結(jié)果如圖9所示。

圖9 雙層介質(zhì)模型不同時(shí)刻的波場(chǎng)快照

圖9中可以看出，在未遇到界面前，地震波在均勻介質(zhì)中的波前面一個(gè)圓。當(dāng)?shù)卣鸩ㄈ肷涞絻煞N不同介質(zhì)的分界面時(shí)，會(huì)產(chǎn)生透射波和反射波，符合波動(dòng)理論。

4 結(jié) 論

1)針對(duì)數(shù)值模擬過(guò)程中存在的邊界反射問(wèn)題，可以通過(guò)引入恰當(dāng)?shù)倪吔鐥l件來(lái)解決，盡最大可能去削弱邊界反射，改善正演模擬的效果。

2)空間采樣間隔越大，計(jì)算速度越快，相應(yīng)的數(shù)值頻散現(xiàn)象越嚴(yán)重，因此，要合理選擇空間步長(zhǎng)，以便同時(shí)兼顧抑制頻散和提高計(jì)算效率。

3)震源子波頻率對(duì)模擬效果也有一定的影響。子波頻率越高，分辨率也會(huì)越高，但數(shù)值頻散現(xiàn)象也越嚴(yán)重，因此需要合理選擇震源子波的主頻，在保證具有較高分辨率的同時(shí)可以充分壓制數(shù)值頻散。

4)對(duì)不同時(shí)刻波場(chǎng)快照分析說(shuō)明，利用聲波波動(dòng)方程計(jì)算地震波場(chǎng)具有很好的運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)特征。為深入了解地震波在地下介質(zhì)中地震波傳播規(guī)律的研究提供理論依據(jù)。

[1] 馮英杰, 楊長(zhǎng)春, 吳萍. 地震波有限差分綜述[J]. 地球物理學(xué)進(jìn)展, 2007, 22(2): 487-491.

[2] 王毓軍. 二維地震波場(chǎng)有限差分正演模擬及合成地震記錄的實(shí)現(xiàn)[D]. 桂林: 桂林理工大學(xué), 2011.

[3] 梁超, 呂力行, 童祥軒. 采空區(qū)頂板動(dòng)力響應(yīng)數(shù)值模擬研究[J]. 中國(guó)錳業(yè), 2016, 34(4): 61-67.

[4] Robert C, Born E. Absorbing boundary conditions for acoustic and clastic wave equation[J]. Bulletin of the Scismalagical Society of America, 1977, 66(11): 1529-1540.

[5] Reynolds A C. Boundary conditions for the numerical solution of wave propagation problems[J]. Geophysics, 1978, 43(6): 1099-1110.

[6] Jean-pierre B. A perfectly matched layer for the absorption of electromagnetic waves[J]. Jouenal of Computational physics, 1994(114): 185-200.

[7] Hastings F D, Schneider J B. Application of the perfectly matched layer[J]. Journal AcousricalSociety of America, 1996, 100(5): 3051-3068.

[8] 倪長(zhǎng)寬. 復(fù)雜介質(zhì)地震波場(chǎng)正演模擬方法研究[D]. 青島: 中國(guó)石油大學(xué), 2008.

[9] 劉慶敏. 高階差分?jǐn)?shù)值模擬方法研究[D]. 青島: 中國(guó)石油大學(xué), 2007.

[10] 李勝軍. 地震波數(shù)值模擬中的頻散壓制方法分析[J]. 石油物探, 2008, 47(5): 444-449.