基于等效交錯網(wǎng)格的流固耦合介質(zhì)地震波模擬

吳建魯　吳國忱　王　偉　何　京

(①中國石油大學(xué)(華東)，山東青島 266580； ②海洋國家實(shí)驗(yàn)室海洋礦產(chǎn)資源評價與探測技術(shù)功能實(shí)驗(yàn)室，山東青島 266071； ③中海石油(中國)有限公司天津分公司渤海石油研究院,天津塘沽 300459)

1　引言

流固耦合介質(zhì)在現(xiàn)實(shí)世界中普遍存在，包括海洋環(huán)境、河流、壩體和船舶等。地震波在流固耦合介質(zhì)中的數(shù)值模擬，一直是研究熱點(diǎn)和難點(diǎn)。Hou等[1]將流固耦合介質(zhì)地震波數(shù)值模擬方法分為兩種： ①單一方程法，即在流相和固相中采用相同的地震波波動方程，利用一階速度—應(yīng)力方程針對流固耦合介質(zhì)進(jìn)行地震波數(shù)值模擬[2-7]，不必考慮界面的幾何特征，實(shí)現(xiàn)方法簡單，但是由于方程中的變量較多(各向同性條件下二維情況為5個，三維情況為9個)，計算效率相對較低； ②雙方程耦合方法，即分別用聲波方程和彈性波方程描述地震波在液相和固相介質(zhì)中的傳播，并利用連續(xù)性條件將兩種介質(zhì)聯(lián)系起來，但需要精確考慮流固耦合介質(zhì)的界面信息。Stephen[5，8]和Sochackid等[9]利用流固界面的應(yīng)力、應(yīng)變連續(xù)性條件采用有限差分方法對聲學(xué)和彈性介質(zhì)的純位移方程耦合進(jìn)行了正演模擬，雖然短時間內(nèi)波場傳播十分穩(wěn)定，但隨著傳播時間的增加，界面耦合近似誤差積累逐漸導(dǎo)致地震波傳播不穩(wěn)定； 另外，界面處耦合差分方法很難適用于高階格式，極大地限制了空間差分步長的選擇。Lee等[10]結(jié)合流固界面連續(xù)性條件利用單元網(wǎng)格有限差分?jǐn)?shù)值法模擬了水平界面情況下的地震波傳播，同樣該方法很難推廣到高階差分格式，限制了實(shí)際應(yīng)用。曲英明等[11]和Qu等[12]利用一階速度—應(yīng)力方程將該方法推廣到海底起伏模型情況，克服了起伏界面模擬時的階梯散射波場干擾。Komatitsch等[13]利用譜元法數(shù)值模擬了水平及起伏海底條件的地震波傳播，模擬結(jié)果相對穩(wěn)定； Carcione等[14]利用偽譜法分析了海底不同形式波的傳播特征； Zhang[15]結(jié)合有限差分和有限元方法進(jìn)行了單元格子法流固界面地震波數(shù)值模擬方法研究，該方法對起伏不規(guī)則海底界面適用性相對較強(qiáng)； 宣領(lǐng)寬等[16]考慮了四邊形單元的雙線性特性并結(jié)合變步長思想進(jìn)一步提高了數(shù)值模擬精度。Choi等[17]在二階聲壓波動理論和彈性位移波動理論框架下采用頻率域有限元方法模擬了流固耦合介質(zhì)地震波傳播，Bae等[18，19]將該方法應(yīng)用于不同形式下的全波形反演，進(jìn)一步驗(yàn)證了方法的穩(wěn)定性和有效性。Yu等[20，21]通過聲壓與應(yīng)力參數(shù)的定量關(guān)系將傳統(tǒng)的一階速度—應(yīng)力方程變換為一階速度—應(yīng)力—壓力方程，并用于海上各向同性和VTI介質(zhì)的逆時偏移。

上述方法中采用一階速度—應(yīng)力波動方程模擬精度相對較高且流固界面模擬穩(wěn)定，但是在現(xiàn)有計算環(huán)境下進(jìn)行實(shí)際應(yīng)用時，勢必會面臨計算量大、計算效率相對較低的問題[17，22-24]。基于有限元、譜元和單元格子法采用雙方程耦合方法可以有效地模擬流固邊界處的地震波傳播，但是該方法需要進(jìn)行精確的網(wǎng)格剖分，尤其針對復(fù)雜三維情況下的剖分需要特別留意，以免引入不必要的誤差[25-27]。Bartolo等[28，29]借助交錯網(wǎng)格差分思想提出了針對二階聲壓標(biāo)量波方程和二階純位移波動方程的時間域等效交錯差分正演模擬，該方法不但和一階交錯差分具有同樣的差分精度，而且含有的參數(shù)變量較少，節(jié)省了大量內(nèi)存，提高了正演模擬效率。本文借助等效交錯網(wǎng)格思想并充分考慮密度參數(shù)空間變化對地震波傳播的影響，在流相介質(zhì)和固相介質(zhì)中分別采用非均質(zhì)情況二階聲壓標(biāo)量波、二階純位移控制方程，為提高地震波在流固相介質(zhì)中傳播的精度和穩(wěn)定性，在流固耦合界面處采用一階位移—應(yīng)力彈性波動方程作為流相和固相之間的轉(zhuǎn)換過渡層，并詳細(xì)說明了過渡層與上、下介質(zhì)進(jìn)行時間和空間差分的耦合方法。利用該方法對不同模型進(jìn)行測試并與一階交錯網(wǎng)格速度—應(yīng)力方程的結(jié)果進(jìn)行對比，證明了方法的準(zhǔn)確性和穩(wěn)定性。

2　流固耦合介質(zhì)正演模擬方法

2.1　流固耦合介質(zhì)波動方程

在各向同性非均質(zhì)情況下，流體中剪切應(yīng)力為零，一般用聲壓參數(shù)和位移矢量進(jìn)行刻畫，流相中的物理參數(shù)一般包括兩個：模量(K)和密度(ρ)，通過一階速度—壓力方程可得非均質(zhì)情況下聲壓標(biāo)量波的波動方程

(1)

其二維形式為

(2)

式中：P為壓力；K為液相的體積模量；S(x，t)為震源函數(shù)；ρ為密度。

在非均質(zhì)條件下各向同性二階純位移波動方程為

(3)

其二維形式為

(4)

式中：u(x，t)=[U(x，z，t)，W(x，z，t)]為位移向量；λ和μ為拉梅系數(shù)；F=(fx，fz)為震源函數(shù)。

2.2　等效交錯網(wǎng)格思想

以二維為例，如圖1所示，傳統(tǒng)的交錯網(wǎng)格差分格式是在空間和時間上應(yīng)力(τxx，τzz，τxz)和速度(vx，vz)交錯更新，每次都需要保存至少兩個時刻的所有參數(shù)值，而等效交錯網(wǎng)格在利用其周圍應(yīng)力或應(yīng)變信息時，又可以采用其本身的差分格式來代替，從而實(shí)現(xiàn)僅含有純壓力(應(yīng)力)或純位移交錯差分思想。式(4)中x方向空間四階、時間二階的差分格式為

圖1　傳統(tǒng)交錯和等效交錯網(wǎng)格差分格式示意圖

(5)

(6)

式中：Δx和Δz分別為x和z方向的網(wǎng)格間隔；p為介質(zhì)的彈性參數(shù)或彈性參數(shù)的組合。式(2)和式(4)中的其他偏微分項(xiàng)的差分格式與式(5)和式(6)類似，不再贅述。

2.3　基于交錯網(wǎng)格流固界面過渡層

借助交錯網(wǎng)格可以穩(wěn)定數(shù)值模擬地震波在流固界面?zhèn)鞑サ奶攸c(diǎn)，在流固界面處采用一階位移—應(yīng)力波動方程差分模擬，與上覆流相中的聲壓標(biāo)量波方程和下伏固相中的純位移彈性波方程的等效交錯網(wǎng)格格式耦合差分，如圖2所示。假設(shè)流固界面為水平界面，Δt為時間間隔，NfΔz為流固界面的深度，從深度(Nf-4)Δz到(Nf+5)Δz為常規(guī)交錯網(wǎng)格的區(qū)域。假設(shè)已知(n-1)Δt時刻從0到NfΔz之間的壓力分量、從(Nf-10)Δz到(Nf+8)Δz之間的應(yīng)力分量和從(Nf-6)Δz到NzΔz處位移分量(Nz為模型z方向網(wǎng)格數(shù))，首先應(yīng)用應(yīng)力—位移關(guān)系(式(8))和彈性波位移方程(式(4))分別更新從(Nf-10)Δz到(Nf+4)Δz之間和從(Nf+5)Δz到NzΔz之間nΔt時刻的位移分量，再利用應(yīng)力—位移關(guān)系式(式(8))更新從(Nf-3)Δz到(Nf+8)Δz之間nΔt時刻的應(yīng)力分量，利用式(2)和壓力—應(yīng)力連續(xù)方程(式(7))分別更新從0到(Nf-4)Δz之間和從(Nf-3)Δz到NfΔz之間nΔt時刻的壓力分量，最后利用式(7)更新(Nf-10)Δz到(Nf-4)Δz之間nΔt時刻的應(yīng)力分量，從而完成nΔt時刻所有變量的更新工作。利用上述流程可完成所有時刻的變量計算工作，進(jìn)而完成正演模擬工作。

(7)

圖2　基于交錯網(wǎng)格過渡層與上、下等效交錯網(wǎng)格耦合過程示意圖

(8)

流固界面處的常規(guī)交錯網(wǎng)格和上、下等效交錯網(wǎng)格之間的銜接部分同樣采用交錯網(wǎng)格差分格式，實(shí)現(xiàn)了模型整個區(qū)域都采用交錯網(wǎng)格差分思想，一方面保證了地震波在流、固相之間傳播的穩(wěn)定性，同時也保證了數(shù)值模擬的精度。本文在流相和固相中采用二階常規(guī)PML吸收邊界條件，在過渡層中采用一階常規(guī)PML吸收邊界條件。本文方法可以推廣到更高階的差分格式。

3　數(shù)值模擬結(jié)果

3.1　流固雙層模型

首先，通過流固雙層模型數(shù)值模擬驗(yàn)證本文正演方法的優(yōu)勢。圖3為流固雙層速度模型示意圖，五角星代表炮點(diǎn)位置，距離水面60m，三角形代表聲壓分量接收點(diǎn)位置，距離水面120m。另外在海底接收水平和垂直位移分量，x與z方向的網(wǎng)格間隔均為6m，選用主頻為20Hz的雷克子波激發(fā)，時間采樣間隔為0.2ms，記錄時間為1.2s。表1為本文等效交錯網(wǎng)格(ESG)方法與傳統(tǒng)交錯網(wǎng)格差分(SSG)方法運(yùn)行內(nèi)存和計算時間的對比，可見ESG方法極大地降低了正演模擬的內(nèi)存需求，約為SSG方法的三分之一。另一方面，由于耦合情況下的ESG方法只需在液相和固相中分別使用純聲壓和位移彈性波動方程，提高了計算效率，ESG方法相比于SSG方法運(yùn)行時間縮短了四分之一。

表1　兩種方法內(nèi)存及運(yùn)行時間對比

同時為驗(yàn)證本文等效交錯方法(ESG)數(shù)值模擬的準(zhǔn)確性，對圖3中的模型進(jìn)行細(xì)網(wǎng)格剖分以便求取更為精確的數(shù)值模擬結(jié)果(可近似為解析解)，空間網(wǎng)格間距為1m， 時間間隔為0.05ms， 模型大小和炮點(diǎn)信息不變。圖4為粗、細(xì)網(wǎng)格模擬的聲壓分量、水平位移及垂直位移記錄對比，兩種網(wǎng)格獲得的模擬結(jié)果基本相同。為進(jìn)一步研究兩種方法模擬結(jié)果的波形一致性，抽取x方向900m處的單道進(jìn)行比較(圖5)。由于正演模擬中不可避免的空間頻散和時間頻散造成了二者振幅和相位的細(xì)微差異，但在誤差允許的范圍之內(nèi)，說明了等效交錯網(wǎng)格方法的準(zhǔn)確性。圖6為等效交錯方法模擬的0.5s時刻流固雙層模型的波場快照，液相為壓力波場，固相為位移波場。流相和固相中的反射、透射波以及首波的模擬結(jié)果十分清晰，從另一個方面也說明了本文針對流固耦合介質(zhì)地震波數(shù)值模擬方法的準(zhǔn)確性。

圖3　流固雙層模型

圖4　　　　雙層模型6m(左)和1m(右)網(wǎng)格間距ESG方法模擬結(jié)果對比

圖5　　　　雙層模型6m(黑色)和1m(紅色)網(wǎng)格間距ESG方法模擬的x方向900m處地震記錄

圖6　雙層模型ESG方法模擬的0.5s波場快照

3.2　Marmousi 2模型

為進(jìn)一步驗(yàn)證和說明本文方法對復(fù)雜模型的適用性和穩(wěn)定性，采用復(fù)雜的Marmousi 2模型(圖7)進(jìn)行測試。需要特別指出的是，為保證模型的模擬精度，本文未考慮原始模型中海底松軟低速層的影響。模型網(wǎng)格數(shù)為2401×1201，兩個方向的網(wǎng)格間距均為5m，水層的厚度為900m，在水面下50m激發(fā)主頻為20Hz的雷克子波，檢波器的設(shè)置在水下60m。圖8為炮點(diǎn)位于x=6000m處利用等效交錯網(wǎng)格(ESG)和傳統(tǒng)交錯網(wǎng)格(SSG)方法模擬的地震記錄，兩種方法都采用時間二階和空間四階差分格式。圖9為檢波器位于x=5000m處單道記錄的對比，可以發(fā)現(xiàn)，兩種方法模擬的地震信號基本一致。圖10為Marmousi 2模型2.0s時刻的波場快照，可見本文提出的流固介質(zhì)等效交錯網(wǎng)格方法在復(fù)雜模型中對變速構(gòu)造的模擬結(jié)果穩(wěn)定，且與傳統(tǒng)交錯網(wǎng)格方法的空間精度相同。

圖7　Marmousi 2模型

圖8　ESG方法(a)和SSG方法(b)模擬的Marmousi 2模型單炮記錄

圖9　　　　ESG方法和SSG方法模擬的Marmousi 2模型單道記錄(a)及其局部放大(b)

圖10　Marmousi 2模型ESG方法模擬的2.0s波場快照

4　分析與結(jié)論

本文借助等效交錯網(wǎng)格的思想，給出了流固耦合介質(zhì)正演模擬方法，流相介質(zhì)由聲壓標(biāo)量波方程控制，固相介質(zhì)由二階純位移彈性波方程控制，在流固界面附近采用傳統(tǒng)的一階位移—應(yīng)力交錯網(wǎng)格作為過渡層，從而實(shí)現(xiàn)流固介質(zhì)地震波正演模擬。

(1)與傳統(tǒng)的交錯網(wǎng)格正演模擬相比，本文提出的流固耦合介質(zhì)正演模擬在保持與傳統(tǒng)交錯網(wǎng)格方法相同模擬精度的前提下，盡可能地減少了正演過程中的變量，降低了內(nèi)存需求，提高了計算效率。

(2)通過簡單流固雙層模型和復(fù)雜Marmousi 2模型試算，不僅驗(yàn)證了等效交錯網(wǎng)格和傳統(tǒng)交錯網(wǎng)格具有相同的模擬精度，而且說明了本文所述方法的準(zhǔn)確性和穩(wěn)定性。

(3)該方法同樣可以推廣到三維情況下，并且相對于傳統(tǒng)的交錯網(wǎng)格有限差分方法，本文的流固耦合地震波數(shù)值模擬方法對于海上的地震偏移和全波形反演等的研究和應(yīng)用有較好的適用性。

(4)當(dāng)海底界面為非水平時，則必須增加過渡層的層數(shù)，使得海底界面在過渡層之內(nèi)即可。

[1]Hou G，Wang J and Layton A.Numerical methods for f-structure interaction — a review.Computer Physics Communications，2012，12(2):337-377.

[2]Virieux J.P-SV wave propagation in heterogeneous media：Velocity-stress finite-difference method.Geophysics，1986，51(4):889-901.

[3]Burns D R and Stephen R A.Three-dimensional numerical modeling of geoacoustic scattering from seafloor topography.The Journal of the Acoustical Society of America，1990，88(5)：2338-2345.

[4]Dougherty M E and Stephen R A.Seismo/acoustic propagation through rough seafloors.The Journal of the Acoustical Society of America，1991，90(5):2637-2651.

[5]Stephen R A.A review of finite difference methods for seismo-acoustics problems at the seafloor.Reviews of Geophysics，1988，26(3)：445-458.

[6]Van Vossen R，Robertsson J and Chapman C.Finite-difference modeling of wave propagation in a f-solid configuration.Geophysics，2002，67(2)：618-624.

[7]De Basabe J D and Sen M K.A comparison of finite-difference and spectral-element methods for elastic wave propagation in media with a f-solid interface.Geophysical Journal International，2015，200(1)：278-298.

[8]Stephen R A.A comparison of finite difference and reflectivity seismograms for marine models.Geophy-sical Journal International，1983，72(1)：39-57.

[9]Sochacki J S.Interface conditions for acoustic and elastic wave propagation.Geophysics，1991，56(2)：168-181.

[10]Lee H，Lim S and Min D.2D time-domain acoustic-elastic coupled modeling：a cell-based finite-difference method.Geosciences Journal，2010，13(4):407-414.

[11]曲英銘，黃建平，李振春等.基于單元交錯網(wǎng)格的變坐標(biāo)系正演模擬方法在聲—彈介質(zhì)中的應(yīng)用.石油物探，2015，54(5)：582-591.

Qu Yingming，Huang Jianping and Li Zhenchun et al.The application of variable coordinate system forward modeling using cell-based staggered grids in acoustic-elastic coupled medium.GPP，2015，54(5)：582-591.

[12]Qu Y，Huang J and Li Z et al.A hybrid grid method in an auxiliary coordinate system for irregular f-solid interface modeling.Geophysical Journal International，2017，208(3)：1540-1556.

[13]Komatitsch D，Barnes C and Tromp J.Wave propagation near a fluid-solid interface：A spectral-element approach.Geophysics，2000，65(2):623-631.

[14]Carcione J M and Helle H B.The physics and simulation of wave propagation at the ocean bottom.Geophysics，2004，69(3)：825-839.

[15]Zhang J.Wave propagation across f-solid interfaces：a grid method approach.Geophysical Journal International，2004，159(1)：240-252.

[16]宣領(lǐng)寬，張文平，明平劍等.基于不同時間步長時域非結(jié)構(gòu)有限體積法模擬聲—彈性耦合問題.固體力學(xué)學(xué)報，2013，34(2)：158-168.

Xuan Lingkuan，Zhang Wenping and Ming Pingjian et al.Numerical simulation of acoustic-elastic interaction based on muti-time-step unstructured finite volume time domain method.Chinese Journal of Solid Mechanics.2013，34(2)：158-168.

[17]Choi Y，Min D J and Shin C.Two-dimensional waveform inversion of multi-component data in acoustic-elastic coupled media.Geophysical Prospecting，2008，56(6)：863-881.

[18]Bae H S，Shin C and Cha Y et al.2D acoustic-elastic coupled waveform inversion in the Laplace domain.Geophysical Prospecting，2010，58(6)：997-1010.

[19]Bae H S，Pyun S and Chuang W et al.Frequency-domain acoustic-elastic coupled waveform inversion using the Gauss-Newton conjugate gradient method.Geophysical Prospecting，2011，60(3)：413-432.

[20]Yu P，Geng J and Li X et al.Acoustic-elastic coupled equation for ocean bottom seismic data elastic reverse time migration.Geophysics，2016，81(5)：S333-S345.

[21]Yu P，Geng J and Li X et al.Separating quasi-P-wave in transversely isotropic media with a vertical symmetry axis by synthesized pressure applied to ocean-bottom seismic data elastic reverse time migration.Geophysics，2016，81(6)：C295-C307.

[22]黃金強(qiáng)，李振春和黃建平等.TTI介質(zhì)Lebedev網(wǎng)格高階有限差分正演模擬及波型分離.石油地球物理勘探，2017，52(5)：915-927.

Huang Jinqiang，Li Zhenchun and Huang Jianping et al.Lebedev grid high-order finite-difference modeling and elastic wave-mode separation for TTI media.OGP，2017，52(5)：915-927.

[23]高正輝，孫建國，孫章慶等.基于完全匹配層構(gòu)建新方法的2.5維聲波近似方程數(shù)值模擬.石油地球物理勘探，2016，51(6)：1128-1133.

Gao Zhenghui，Sun Jianguo and Sun Zhangqing et al.2.5D acoustic approximate equation numerical simulation with a new construction method of the perfectly matched layer.OGP，2016，51(6)：1128-1133.

[24]梁文全，王彥飛，楊長春.聲波方程數(shù)值模擬矩形網(wǎng)格有限差分系數(shù)確定法.石油地球物理勘探，2017，52(1)：56-62.

Liang Wenquan，Wang Yanfei and Yang Changchun.Acoustic wave equation modeling with rectangle grid finite difference operator and its linear time space domain solution.OGP，2017，52(1)：56-62.

[25]Xu K and McMechan G.2D frequency-domain elastic full-waveform inversion using time-domain modeling and a multi step-length gradient approach.Geophy-sics，2014，79(2)：R41-R53.

[26]趙慧，劉穎，李予國.自適應(yīng)有限元海洋大地電磁場二維正演模擬.石油地球物理勘探，2014，49(3)：578-585.

Zhao Hui，Liu Ying and Li Yuguo.Adaptive finite element forward modeling for two-dimensional marine magnetotelluric fields.OGP，2014，49(3)：578-585.

[27]黃一凡，胡祥云，韓波.自適應(yīng)有限元法的大地電磁二維各向異性正反演.石油地球物理勘探，2016，51(4)：809-820.

Huang Yifan，Hu Xiangyun and Han Bo.Forward and inverse modeling of the magnetotelluric field in 2D anisotropic media with an adaptive finite element method.OGP，2016，51(4)：809-820.

[28]Bartolo L D，Dors C and Mansur W.A new family of finite-difference schemes to solve the heterogeneous acoustic wave equation.Geophysics,2012,77(5)：T187-T199.

[29]Bartolo L D，Dors C and Mansur W.Theory of equi-valent staggered-grid schemes：application to rotated and standard grids in anisotropic media.Geophysical Prospecting，2015，63(5)：1097-1125.