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

(1. 長江大學 油氣資源與勘探技術教育部重點實驗室，湖北 武漢 430100； 2. 長江大學 地球物理與石油資源學院，湖北 武漢 430100； 3. 中石化西北油田分公司采油一廠，新疆 輪臺 841600)

0 前 言

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

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

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

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

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

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

(1)

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

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

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

(2)

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

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

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

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

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

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

速度水平分量：

速度垂直分量：

應力Pxx：

應力Qzz：

應力Sxz：

2 有限差分數(shù)值模擬中的關鍵技術及解決方案

2.1 模擬震源

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

雷克子波表達式(時間域)為：

(3)

通過傅立葉變換為頻率域：

(4)

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

2.2 邊界條件

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

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

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

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

1)透明邊界條件

(5)

2)Clayton-Engquist 邊界條件

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

(6)

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

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

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

由圖3可知，使用了邊界條件的波場快照在邊界處邊界反射波得到明顯的削弱，而沒有進行邊界處理的邊界反射波很明顯。通過引入恰當?shù)倪吔鐥l件，可以有效地削弱邊界反射。

2.3 穩(wěn)定性

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

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

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

(7)

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

(8)

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

(9)

因為vx

(10)

2.4 數(shù)值頻散問題

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

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

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

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

圖5 不同主頻波場快照

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

3 有限差分數(shù)值模擬實例

3.1 均勻介質(zhì)模型

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

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

模擬結果如圖7所示。

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

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

3.2 雙層介質(zhì)模型

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

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

模擬結果如圖9所示。

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

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

4 結 論

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

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

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

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

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

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

[3] 梁超, 呂力行, 童祥軒. 采空區(qū)頂板動力響應數(shù)值模擬研究[J]. 中國錳業(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ì)地震波場正演模擬方法研究[D]. 青島: 中國石油大學, 2008.

[9] 劉慶敏. 高階差分數(shù)值模擬方法研究[D]. 青島: 中國石油大學, 2007.

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