基于旋轉(zhuǎn)交錯(cuò)網(wǎng)格和FCT技術(shù)的高精度數(shù)值模擬

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

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

0 前 言

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

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

(2)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

其具體步驟如下：

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

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

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

4)反彌散通量的計(jì)算

a 計(jì)算k+1時(shí)刻的彌散通量P和Q；

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

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

2 模型試算

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

設(shè)計(jì)一個(gè)均勻各向同性介質(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，時(shí)間步長為0.5 ms，在模型(1 000 m，1 000 m)處使用震源信號(hào)主頻為30 Hz的Ricker子波激發(fā)。采用精度為o(△t2+△x4)的旋轉(zhuǎn)交錯(cuò)網(wǎng)格進(jìn)行數(shù)值模擬，得到X分量在t=0.2 s時(shí)的波場快照，如圖3所示，得到相應(yīng)時(shí)間切片相同位置x分量記錄振幅如圖4所示。

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

比較圖4(a)、(c)、(d)可得，利用通量傳輸校正技術(shù)后，數(shù)值頻散(尤其是橫波的數(shù)值頻散)均得到了有效地壓制，分辨率提高。能夠清晰地反映通量校正傳輸技術(shù)的校正結(jié)果與漫射因子和反射因子有關(guān)，比較圖4(c)、(d)、(f)，可得出壓制效果(d)>(f)>(c)，分析可得漫射因子的取值沒有特定標(biāo)準(zhǔn)，和數(shù)值頻散的大小有關(guān)，當(dāng)數(shù)值頻散較為嚴(yán)重時(shí)，漫射因子取較大值時(shí)，頻散壓制效果較為出色。當(dāng)漫射因子取值達(dá)到合適值時(shí)，頻散壓制效果不再隨著漫射因子的變化而發(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)交錯(cuò)網(wǎng)格不同參數(shù)波場快照對(duì)比

(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)交錯(cuò)網(wǎng)格不同參數(shù)振幅記錄對(duì)比

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

圖5 3種方法振幅記錄對(duì)比

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

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

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

表2 運(yùn)行效率和資源占用對(duì)比

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

設(shè)計(jì)了水平層狀介質(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，時(shí)間步長為1 ms，在模型(1 000 m，1 000 m)處使用震源信號(hào)主頻為30 Hz的Ricker子波激發(fā)，得到速度X分量在t=0.5 s時(shí)的波場快照，如圖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道記錄振幅進(jìn)行對(duì)比，如圖7所示。

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

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

3 結(jié) 論

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