趙 威，李偉波，桂志先，周 游，于曉東

(長江大學 油氣資源與勘探技術教育部重點實驗室，湖北 武漢 430100)

0 前 言

地震波正演模擬技術一直被人們視為一種比較直觀、形象地了解和認識波在地下介質中傳播特性的重要手段。所謂正演模擬，就是在首先假定地下介質模型和相應的物理參數已知的情況下，模擬地震波在地下各種介質中的傳播規(guī)律，從而根據其計算出的地面或地下各觀測點所應觀測到的數值地震記錄來進行分析地下模型的一種地震模擬方法。地震數值模擬由來已久、發(fā)展迅速，目前已經有多種地震數值模擬方法在地震勘探和地震學中得到廣泛而有效的應用。這些地震波場數值模擬方法主要可以歸納為3大類，即幾何射線法、積分方程法和波動方程法[1]。

波動方程法數值模擬本質上是求解地震波動方程，模擬的地震記錄包含了地震波傳播的所有信息，能夠很好的了解地下情況，從而被廣泛應用。波動方程的求解方法有解析法和數值法兩種，而波動方程的解很難以用解析形式表達出來，因此，數值法成為求解波動方程的主要方法。波動方程數值模擬的方法主要有有限差分法、有限元法和偽譜法等。相比其他兩種方法，有限差分作為一種最常用的數值模擬方法，能夠較準確的描述地下介質分布狀況及地震波在各種介質中的傳播規(guī)律。具有計算效率高，占用內存小等特點，但由于其本身時間和空間都采用差分離散形式，往往會出現(xiàn)數值頻散問題[2]。怎樣減少頻散，提高模擬精度始終是有限差分數值模擬的關鍵問題。為此，學者們做了大量工作，當信號在一個波長范圍內，具有多個離散節(jié)點是，頻散現(xiàn)象基本消失。一般隨著網格間距增大，頻散問題越嚴重；反之若減小網格間距相應會增加網格點數，將會大大增加計算量，造成計算時間過長。既能控制模擬精度，又能節(jié)省計算時間的方法顯得尤為重要。Boris和Book等人首先提出了通量校正傳輸的方法，它最早起源于求解流體力學連續(xù)方程中，隨后才開始被運用于聲波方程數值模擬的求解；國內楊頂輝等[3]將其應用于彈性波動方程正演數值模擬頻散處理，取得一定效果。為了更好的對比分析未加FCT數值模擬和加了FCT 方法的區(qū)別，本文基于聲波和彈性波兩方面，借助于交錯網格的模擬優(yōu)勢，對各向同性介質波場進行了有限差分數值模擬，并對其相應的波場快照和計算效率進行對比分析。

1 波動方程有限差分算法原理

1.1 聲波方程形式

二維二階標量聲波方程的一般形式為：

(1)

經過化簡變形可得到一階壓力—速度方程組形式：

(2)

針對(1)式左右兩邊使用二階中心差商公式可得到相應的差分形式：

(3)

其相應的一階壓力—速度方程差分形式為：

Δx=Δz=h

(4)

1.2 一階速度—應力彈性波方程

在二維各向同性介質中，不考慮體力的影響下，彈性波方程可以化為一階速度—應力的形式：

(5)

方程組(5)中，Vx、Vz表示x、z方向速度分量，ρ表示介質的密度，σxx、σxz、σzz分別表示應力分量，C11、C13、C14為介質的彈性參數。

為了研究方便，我們采用時間二階、空間四階交錯網格有限差分形式[4]，交錯網格的網格化和常規(guī)有所區(qū)別，交錯網格引進了半網格，速度和應力分別在整網格和半網格上計算，

時間的2M階差分近似：

(6)

上面兩式相減，得2M階精度的差分近似式：

(7)

當M取1時，上式即為時間二階精度差分格式，帶入示例并整理得(8)：

(8)

空間差分格式近似：

把時間二階差分精度公式和空間四階公式代入方程組(5)，可得一階速度—應力彈性波方程

離散格式：

(10)

2 通量校正FCT方法

通量校正傳輸法(FCT)最初被廣泛應用于流體力學中沖擊波前陣面極大濃度梯度過程中，滿足連續(xù)守恒方程式，守恒定律保證了其方程式和物理定律的一致性。此后，隨著應用的推廣被引用于聲波方程中，在地震波在模型中的傳播方程中有效的壓制了在構造交錯網格差分格式時計算產生的數值頻散，即網格交界面邊界在同一時刻有流入和流出等量的物理通量，用這種方法可以有效的避免了截斷誤差的累積。FCT方法主要包含兩個步驟即：交錯網格有限差分計算和通量的校正[5]。在其通量校正的階段，用改變波場的方法來壓制偽波動。但由于在實際操作過程中并沒有事先證明散射存在的位置，所以FCT首先假設所有極值均由數值頻散引起，從而對全部的網點進行通量校正，在校正的過程中同時對反漫射進行修正，進而使其結果更加準確。以速度分量vx為例，其具體步驟如下。

1)建立各項同性介質條件下各項同性非均勻介質彈性波動方程，進而得到一階速度—應力交錯網格有限差分方程。由FCT算法的特性可知，在對速度分量vx、vz進行計算后，只要速度分量被校正，其應力分量也會被隨之校正。

2)k-1時刻漫射通量的計算：

(11)

其中，η1取值范圍為0～1，為漫射因子，其值可以是一個常數，根據實際情況可以由數值實驗得到。

3)k+1時刻漫射通量的計算。

其中，η2取值范圍與η1一致，為反漫射因子，通常由于反漫射運算需要補償人為加入漫射引起的誤差和有限差分運算帶來的振幅損失，所以一般取值略大于漫射因子η1。

4)對k-1時刻漫射用通量P對速度分量進行平滑修正。

(13)

5)求解k+1時刻的反漫射通量。

(14)

6)利用反漫射通量進行修正，得到校正之后的速度分量：

(15)

其中：

(16)

3 數值模擬頻散分析

有限差分數值模擬是將微分方程離散為差分方程進行求解的過程。當空間網格間距或者時間網格大小選取不合適就會引起波形畸變，地震波在傳播過程中，波前形狀就會發(fā)生變化，并且逐漸散開，這就是所謂的數值頻散現(xiàn)象。其本質是離散化波動方程產生的偽波動，這是有限差分方法所固有的本質特征[6]。目前壓制有限差分數值頻散的思路主要三個方面：一是通過控制網格步長，減少網格步長可以明顯提高數值頻散，但網格點數會增加，計算量會顯著增加；二是提高差分階數，但差分階數達到一定時，再提高效果不明顯且計算量很大；最后就是利用通量校正傳輸技術等。

3.1 均勻介質模型

圖1為均勻介質聲波方程有限差分數值模擬的波場快照。所選網格步長如圖1a，從波場快照來看，很明顯出現(xiàn)了數值頻散，會影響模擬精度；圖1b～d是通過減少、縮減網格步長后的波場快照，b圖可以看出消頻效果明顯，但隨著不斷縮減后期效果并不顯著。圖1e為FCT消頻散后的波場模擬快照，可以看出FCT消頻散雖然所用時間較長，但模擬效果較好；與此同時，比較圖1d、e達到相同的精度時，通量校正傳輸的方法所需時間還是相較前者要少，波場清晰度更高，能夠起到一定的效果。

圖1 標量聲波方程有限差分模擬波場快照

3.2 層狀介質模型

圖2、3分別為均勻介質模型和層狀介質模型的時間二階、空間四階交錯網格有限差分數值模擬波場快照，通過對比通量校正前后的波場快照，可以明顯得到：傳統(tǒng)交錯網格有限差分在空間步長較大的情況下，波的模擬過程中出現(xiàn)了出現(xiàn)了偽波動，產生了數值頻散現(xiàn)象，如圖2a、3a；而結合FCT通量校正后的波場快照可以看出(如圖2b、3b)，波的擾動明顯減少，波場更加穩(wěn)定平滑，很好的消除頻散；再者，傳統(tǒng)方法要達到相同精度的波場快照所需時間遠遠超過FCT。總之，F(xiàn)CT消頻散能夠明顯改善數值模擬結果，能夠有很好的壓制頻散的效果。

(a)未消頻散X、Z分量波場快照；(b)消頻散后X、Z分量波場快照(均勻介質模型下)

圖2時間二階、空間四階交錯網格有限差分數值模擬波場快照

(a)未消頻散X、Z分量波場快照；(b)消頻散后X、Z分量波場快照(層狀介質模型下)

圖3時間二階、空間四階交錯網格有限差分數值模擬波場快照

4 結 論

本文在前人已有的成果基礎上，推導了二維聲波方程、時間二階空間四階彈性波有限差分計算公式，并結合通量校正輸出技術(FCT)，給出了具體的實施方法步驟；將FCT運用到基于各向同性介質、層狀介質的波動方程數值模擬中。通過對比數值模擬波場快照，可以看出該方法可以有效壓制各向同性介質有限差分正演模擬過程中因網格剖分而引起的數值頻散。雖然加入通量校正傳輸技術會需要增加計算量，但相比減小步長，增加網格點數來消除數值頻散，對于相同模擬精度下，該方法能夠大大減少計算量；與此同時，提高精度的同時，對于較大采樣間隔的網格劃分能節(jié)省計算量，提高相應計算速度。因此可以總結為，F(xiàn)CT是在同性介質下，地震波數值模擬時壓制頻散問題的一種有效可能的方法。

參考文獻：

[1] 李信富，李小凡，張美根.地震波數值模擬方法研究綜述[J].防災減災工程學報，2007，27(2)：241-245.

[2] 李賓.橫向各向同性介質有限差分法波場模擬方法研究[D].青島：中國石油大學，2009.

[3] 楊頂輝，滕吉文.各向異性介質中三分量地震記錄的FCT有限差分模擬[J].石油地球物理勘探，1997,32(2):181.

[4] 岳曉鵬，白超英，岳崇旺.高階交錯網格有限差分彈性波場模擬的精度分析[J].煤田地質與勘探，2017,45(1):125-130.

[5] 李立平，何兵壽.TTI介質彈性波FCT有限差分數值模擬[J].地球物理學進展，2017,3(4):1584-1590.

[6] 吳國忱，王華忠.波場模擬中的數值頻散分析與校正策略[J].地球物理學進展，2005，20(1)：58-65.