王 婷，段焱文，蔡偉祥，易院平，汪 勇

(1. 長江大學(xué) 油氣資源與勘探技術(shù)教育部重點實驗室，湖北 武漢 430100； 2. 長江大學(xué) 地球物理與石油資源學(xué)院，湖北 武漢 430100； 3. 中石化重慶涪陵頁巖氣勘探開發(fā)有限公司，重慶 涪陵 408014; 4. 中冶集團武漢勘察研究院有限公司，湖北 武漢 430000)

0 前 言

波動方程數(shù)值求解法是一種首先離散網(wǎng)格化計算區(qū)域，然后通過數(shù)值求解描述地震波傳播的波動方程，最終得到波動方程較高精度的近似解析解，從而模擬地震波傳播的方法。有限差分法是一種快速有效求解波動方程的數(shù)值模擬方法，可以方便靈活地應(yīng)用于復(fù)雜的地質(zhì)模型，并且運動速度不受影響，邊界處理簡單。但是，當(dāng)在一個波長內(nèi)空間采樣點較少時，就會產(chǎn)生嚴重的網(wǎng)格頻散現(xiàn)象。譚未一等[1]采用規(guī)則網(wǎng)格有限差分方法模擬了粘彈VTI介質(zhì)中的波場，董良國等[2-4]將交錯網(wǎng)格與高階差分法有效結(jié)合，求解了橫向同性介質(zhì)中一階速度—應(yīng)力彈性波方程，并對其穩(wěn)定性進行了探究?；豇P斌等[5]對各向異性介質(zhì)彈性波進行了高階交錯網(wǎng)格有限差分數(shù)值模擬， Saenger等[6-7]首先提出了旋轉(zhuǎn)交錯網(wǎng)格有限差分算法，旋轉(zhuǎn)交錯網(wǎng)格與常規(guī)交錯網(wǎng)格有限差分不同，旋轉(zhuǎn)交錯網(wǎng)格有限差分是將所有物理量僅定義在兩個不同的位置。若是將應(yīng)力(或者質(zhì)點振動速度)定義于離散網(wǎng)格單元的中心，那么就將質(zhì)點振動速度(或者應(yīng)力)定義于離散網(wǎng)格單元的頂點，并且可將所有彈性模量都定義于相同位置(與應(yīng)力所在位置相同)。因此，在模擬非均勻介質(zhì)時，彈性常數(shù)就不需要進行平均處理，在模擬各向異性介質(zhì)時，彈性常數(shù)就不需要進行內(nèi)插處理[8]，Saenger等對此已經(jīng)進行了有力證明。嚴紅勇等[9]隨后在粘彈TTI介質(zhì)中進行了旋轉(zhuǎn)交錯網(wǎng)格高階有限差分數(shù)值模擬，進一步探究了旋轉(zhuǎn)交錯網(wǎng)格的特性。通量傳輸校正技術(shù)(Flux-corrected transport method，F(xiàn)CT)是一種在增加少量計算量的前提下，有效壓制數(shù)值模擬中造成數(shù)值頻散的技術(shù)。Book等[10]首次將流體動力學(xué)中的通量校正傳輸方法應(yīng)用于聲波方程求解，楊頂輝等[11]將通量校正傳輸方法與正演模擬相結(jié)合，應(yīng)用到彈性波方程正演模擬。均取得了較好的效果。張省等[12]在VTI介質(zhì)中將常規(guī)交錯網(wǎng)格與FCT技術(shù)相結(jié)合進行了數(shù)值模擬，有效壓制了數(shù)值頻散。本文在前人研究的基礎(chǔ)上，首次將通量傳輸校正技術(shù)與旋轉(zhuǎn)交錯網(wǎng)格相結(jié)合，實現(xiàn)了彈性波方程旋轉(zhuǎn)交錯網(wǎng)格高階有限差分數(shù)值模擬，提高了數(shù)值模擬穩(wěn)定性，并利用FCT技術(shù)有效壓制了粗網(wǎng)格情況下造成的數(shù)值頻散，提高了數(shù)值模擬精度。

1 方法原理

1.1 一階速度—應(yīng)力彈性波方程

在二維TI介質(zhì)中，用速度—應(yīng)力表示的彈性波方程(假設(shè)體力為零)[5]為：

(1)

其中，Vi為i(i=x或z)方向速度分量，τij為ij(ij=xx,zz或xz)方向應(yīng)力，ρ為密度，Cij為介質(zhì)的彈性常數(shù)。在各向同性情況下，C11=C33=λ+2μ，C13=λ，C44=μ。

1.2 時間上2M階差分近似

在本文中，使用旋轉(zhuǎn)交錯網(wǎng)格法來數(shù)值求解一階彈性波方程時，速度是在t+△t/2時刻，應(yīng)力是在 時刻進行計算的，我們可以利用泰勒公式展開式來提高時間差分精度，得到2M階精度的時間差分近似[4]：

(2)

其中，?2m-1Vx/?t2m-1為x方向速度對時間的2m-1階偏導(dǎo)數(shù)。

將公式(2)中速度對時間的奇數(shù)階高階導(dǎo)數(shù)用公式(1)中應(yīng)力對空間的導(dǎo)數(shù)表示，同理也可以將應(yīng)力對時間的奇數(shù)階高階導(dǎo)數(shù)用速度對空間的導(dǎo)數(shù)表示。

1.3 旋轉(zhuǎn)交錯網(wǎng)格有限差分格式

在常規(guī)交錯網(wǎng)格中，若將速度分量(或應(yīng)力分量)定義在不同整網(wǎng)格點，應(yīng)力分量(或速度分量)定義在不同半網(wǎng)格點，則此時彈性模量應(yīng)該定義在所有的半網(wǎng)格點和整網(wǎng)格點上，而在實際情況中，通常只將彈性模量定義在半網(wǎng)格點上，故需要對彈性模量進行插值。當(dāng)彈性模量變化很大時，普通交錯網(wǎng)格就容易出現(xiàn)計算不穩(wěn)定的問題。

在旋轉(zhuǎn)交錯網(wǎng)格中，若將應(yīng)力(或者質(zhì)點振動速度)位于離散網(wǎng)格單元的中心，質(zhì)點振動速度(或者應(yīng)力)位于離散網(wǎng)格單元的頂點，則此時所有彈性模量都可定義于相同網(wǎng)格上(與應(yīng)力所在位置相同)。因此，在對非均勻介質(zhì)進行數(shù)值模擬時，就不需要對彈性常數(shù)進行平均處理，在對各向異性介質(zhì)進行數(shù)值模擬時，就不需要對彈性常數(shù)進行內(nèi)插處理。圖1對二維情況下的常規(guī)交錯網(wǎng)格有限差分算法和旋轉(zhuǎn)交錯網(wǎng)格有限差分算法的空間位置進行了解釋。

(a)常規(guī)交錯網(wǎng)格 (b)旋轉(zhuǎn)交錯網(wǎng)格圖1 二維網(wǎng)格單元

旋轉(zhuǎn)交錯網(wǎng)格有限差分算法分兩步完成，首先是計算沿著網(wǎng)格對角線方向相關(guān)物理量的空間導(dǎo)數(shù)，然后線性疊加這些沿對角線方向計算得到的相關(guān)物理量結(jié)果，從而得到沿坐標軸方向的空間導(dǎo)數(shù)[7]。

將旋轉(zhuǎn)交錯網(wǎng)格沿坐標軸方向的空間導(dǎo)數(shù)與彈性波波動方程相結(jié)合，可以得到一階彈性波旋轉(zhuǎn)交錯網(wǎng)格差分格式，對此，Saenger進行了公式推導(dǎo)[6-7]。

1.4 旋轉(zhuǎn)交錯網(wǎng)格與常規(guī)交錯網(wǎng)格的穩(wěn)定性條件

分析一階速度—應(yīng)力彈性波方程組的時間二階，空間2N階精度旋轉(zhuǎn)交錯網(wǎng)格[10]和常規(guī)交錯網(wǎng)格[3]有限差分方程，可得到兩者數(shù)值解的穩(wěn)定性條件(見表1)。

表1穩(wěn)定性條件

1.5 通量傳輸校正技術(shù)

在有限差分法中，數(shù)值頻散包括空間頻散和時間頻散，但空間頻散占其主要方面。旋轉(zhuǎn)交錯網(wǎng)格屬于有限差分法，雖然與常規(guī)網(wǎng)格相比，數(shù)值頻散得到了改善，精度得到了提高，但仍存在有限差分法固有的數(shù)值頻散的問題，網(wǎng)格越大，數(shù)值頻散越嚴重。所以，提高差分方程階數(shù)來提高波場模擬精度和減小空間網(wǎng)格間距來壓制數(shù)值頻散是一種常規(guī)方法。但是，隨著網(wǎng)格間距的縮小，運算效率會降低。通量傳輸校正技術(shù)(FCT)適用于高階差分并有較高的計算精度,在采樣間隔較大時，可以在增加少量計算量的前提下有效壓制數(shù)值模擬中造成的數(shù)值頻散。通量傳輸校正技術(shù)雖然能有效的壓制數(shù)值頻散，但存在較難選取合適參數(shù)的問題，從而影響校正效果。作者將采用改進后的FCT與旋轉(zhuǎn)交錯網(wǎng)格相結(jié)合，借助旋轉(zhuǎn)交錯網(wǎng)格來實現(xiàn)對彈性波方程的有限差分數(shù)值模擬，在增加少量計算機內(nèi)存的情況下有效壓制粗網(wǎng)格造成的數(shù)值頻散。

FCT校正[12]是在假設(shè)所有的極值點都是由數(shù)值頻散引起的前提下，對參與計算的所有網(wǎng)格點進行擴散通量校正處理，再對非局部極值點進行逆擴散通量校正來補償損耗。FCT旋轉(zhuǎn)交錯網(wǎng)格有限差分方法的實現(xiàn)分兩個階段:旋轉(zhuǎn)交錯網(wǎng)格有限差分階段和通量校正階段。通量校正階段又分為兩步：平滑方程的解和反漫射矯正。

圖2 旋轉(zhuǎn)交錯網(wǎng)格FCT

其具體步驟如下：

1)用旋轉(zhuǎn)交錯網(wǎng)格有限差分格式表示彈性波波動方程。求解k+1時刻和k時刻的速度分量。

2)以速度vx為例，計算k時刻的彌散通量P和Q。

3)運用k+1時刻計算得到的彌散通量對所有網(wǎng)格點(包括非局部極值點)上的速度vx進行修正。

4)反彌散通量的計算

a 計算k+1時刻的彌散通量P和Q；

c 進行反彌散通量處理來得到消除數(shù)值頻散后的正確波場值。

反彌散通量的計算不是對所有的網(wǎng)格點進行處理，僅針對進行漫射通量校正前未產(chǎn)生頻散的網(wǎng)格點(即非局部極值點)，通過修正被平滑過的非局部極值點方程的解，恢復(fù)不需要進行漫反射位置(未產(chǎn)生數(shù)值頻散位置)的真實波場。具體計算公式可參考文獻[12]。

2 模型試算

2.1 均勻介質(zhì)VTI介質(zhì)

設(shè)計一個均勻各向同性介質(zhì)模型，其大小為2 000 m×2 000 m，縱波速度vp=3 000 m/s，橫波速度vs=1 732 m/s，密度ρ=1 500 kg/m3，空間網(wǎng)格間距8 m和空間網(wǎng)格間距10 m，時間步長為0.5 ms，在模型(1 000 m，1 000 m)處使用震源信號主頻為30 Hz的Ricker子波激發(fā)。采用精度為o(△t2+△x4)的旋轉(zhuǎn)交錯網(wǎng)格進行數(shù)值模擬，得到X分量在t=0.2 s時的波場快照，如圖3所示，得到相應(yīng)時間切片相同位置x分量記錄振幅如圖4所示。

由圖4(a)和(b)可以看出，在網(wǎng)格大小8 m的情況下，P波波場無明顯的數(shù)值頻散，而S波由于速度較低，其波場模擬結(jié)果數(shù)值頻散較為嚴重，存在明顯的"拖尾"現(xiàn)象。在網(wǎng)格大小10 m的情況下，P波波場有輕微的數(shù)值頻散，而S波波場數(shù)值頻散比小網(wǎng)格嚴重?？梢?，空間網(wǎng)格較小的情況下，頻散得到一定的改善。

比較圖4(a)、(c)、(d)可得，利用通量傳輸校正技術(shù)后，數(shù)值頻散(尤其是橫波的數(shù)值頻散)均得到了有效地壓制，分辨率提高。能夠清晰地反映通量校正傳輸技術(shù)的校正結(jié)果與漫射因子和反射因子有關(guān)，比較圖4(c)、(d)、(f)，可得出壓制效果(d)>(f)>(c)，分析可得漫射因子的取值沒有特定標準，和數(shù)值頻散的大小有關(guān)，當(dāng)數(shù)值頻散較為嚴重時，漫射因子取較大值時，頻散壓制效果較為出色。當(dāng)漫射因子取值達到合適值時，頻散壓制效果不再隨著漫射因子的變化而發(fā)生顯著變化。

(a)△x=△z=8，η1=0，η2=0；(b)△x=△z=10，η1=0，η2=0；(c)△x=△z=8，η1=0.01，η2=0.02；(d)△x=△z=8，η1=0.02，η2=0.03；(e)△x=△z=10，η1=0.01，η2=0.02；(f)△x=△z=10，η1=0.02，η2=0.03

圖3旋轉(zhuǎn)交錯網(wǎng)格不同參數(shù)波場快照對比

(a)△x=△z=8，η1=0，η2=0；(b)△x=△z=10，η1=0，η2=0；(c)△x=△z=8，η1=0.01，η2=0.02；(d)△x=△z=8，η1=0.02，η2=0.03；(e)△x=△z=10，η1=0.01，η2=0.02；(f)△x=△z=10，η1=0.02，η2=0.03

圖4旋轉(zhuǎn)交錯網(wǎng)格不同參數(shù)振幅記錄對比

為了進一步探究在得到相同計算精度的情況下使用FCT和不使用FCT的運行效率和資源占用情況，本文在未使用FCT情況下縮小空間網(wǎng)格間距，若要達到使用FCT(△x2=△z2=10，η1=0.02,η2=0.03)的精度(模型相同)，需要將網(wǎng)格空間大小取到4 m，得到波形對比圖如圖5所示，a實線為旋轉(zhuǎn)交錯網(wǎng)格應(yīng)用FCT后的數(shù)值模擬結(jié)果，b實線為旋轉(zhuǎn)交錯網(wǎng)格數(shù)值模擬結(jié)果，c虛線為常規(guī)交錯網(wǎng)格的精確解析解。

圖5 3種方法振幅記錄對比

從圖5可以看出，三者基本一致，旋轉(zhuǎn)交錯網(wǎng)格

應(yīng)用FCT技術(shù)的模擬結(jié)果略優(yōu)于旋轉(zhuǎn)交錯網(wǎng)格的數(shù)值模擬，旋轉(zhuǎn)交錯網(wǎng)格略優(yōu)于常規(guī)交錯網(wǎng)格。

資源占用量和效率如表2所示。從表2中數(shù)組個數(shù)和數(shù)組大小可以得到，使用FCT所用的存儲空間(數(shù)組個數(shù)×數(shù)組大小)比未使用FCT的存儲空間小，達到相同效果所需計算機運行時間也略小。因此，我們可以看到使用FCT技術(shù)確實提高了效率和節(jié)省了存儲空間。

表2 運行效率和資源占用對比

2.2 水平層狀介質(zhì)模型

設(shè)計了水平層狀介質(zhì)模型(兩種介質(zhì)分界面位于深度1 000 m處)，模型總大小為2 000 m×2 000 m(vp1=3 000 m/s，vs1=1 732 m/s，ρ1=2 000 kg/m3，vp2=4 200 m/s，vs2=2 425 m/s，ρ2=2 500 kg/m3)，空間網(wǎng)格間距△x=△z=8 m，時間步長為1 ms，在模型(1 000 m，1 000 m)處使用震源信號主頻為30 Hz的Ricker子波激發(fā)，得到速度X分量在t=0.5 s時的波場快照，如圖6(a)所示，取漫射因子η1為0.02，反漫射因子η2為0.03，波場快照如圖6(b)所示。

①為第一層介質(zhì)入射S波；②為第一層介質(zhì)入射P波；③為第一層介質(zhì)PP透射波；④為第一層介質(zhì)PP反射波；⑤為第一層介質(zhì)PS透射波；⑥為第一層介質(zhì)PS轉(zhuǎn)換波

圖6水平層狀介質(zhì)波場快照

抽取第125道記錄振幅進行對比，如圖7所示。

圖7 第125道記錄振幅對比(△x1=△z1=8)

觀察圖5和圖6可得，在層狀介質(zhì)模型中進行正演模擬，通量傳輸校正技術(shù)也可以有效壓制數(shù)值頻散，但由于通量校正傳輸方法進行擴散通量校正處理是對所有網(wǎng)格點進行處理，造成了非局部極值點的損失，隨后再對有損失的非局部極值點進行逆散通量校正處理來進行補償損失，所以也造成了一定的振幅偏差。圖7可以十分清晰的看出使用FCT后的效果得到顯著提高。

3 結(jié) 論

旋轉(zhuǎn)交錯網(wǎng)格有限差分數(shù)值模擬方法是一種有效的地震波場數(shù)值模擬方法，可以將地震波在地下介質(zhì)中的傳播過程及規(guī)律清晰明了地給呈現(xiàn)出來，模擬精度較高，穩(wěn)定性相較于常規(guī)網(wǎng)格得到了加強，實現(xiàn)過程較簡單容易。通量傳輸校正技術(shù)可以有效壓制數(shù)值模擬過程中由于網(wǎng)格差分產(chǎn)生的數(shù)值頻散，特別是在大網(wǎng)格情況下很大地提高了計算精度，清晰地呈現(xiàn)出波在地下介質(zhì)中傳播過程，反射因子和漫反射因子的取值沒有唯一標準，取值大小直接影響通量傳輸校正技術(shù)的效果，當(dāng)頻散較為嚴重(網(wǎng)格取值較大)時，漫射因子可以取較大值，當(dāng)產(chǎn)生輕微頻散時，可以取較小值。當(dāng)取值的變化不再影響頻散的壓制效果時，說明取值合適。本文首次將通量傳輸校正技術(shù)與旋轉(zhuǎn)交錯網(wǎng)格相結(jié)合，對均勻介質(zhì)和層狀介質(zhì)模型進行正演模擬，不僅有效壓制數(shù)值頻散，達到較高模擬精度，還減少了增加網(wǎng)格點來提高計算精度所產(chǎn)生的計算量，提高了運算速度。