吳建魯　吳國忱　王　偉　何　京

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

1　引言

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

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

2　流固耦合介質正演模擬方法

2.1　流固耦合介質波動方程

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

(1)

其二維形式為

(2)

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

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

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

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

(5)

(6)

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

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

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

(7)

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

(8)

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

3　數(shù)值模擬結果

3.1　流固雙層模型

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

表1　兩種方法內存及運行時間對比

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

圖3　流固雙層模型

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

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

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

3.2　Marmousi 2模型

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

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

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

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

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

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

[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)格的變坐標系正演模擬方法在聲—彈介質中的應用.石油物探，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]宣領寬，張文平，明平劍等.基于不同時間步長時域非結構有限體積法模擬聲—彈性耦合問題.固體力學學報，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]黃金強，李振春和黃建平等.TTI介質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]高正輝，孫建國，孫章慶等.基于完全匹配層構建新方法的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]趙慧，劉穎，李予國.自適應有限元海洋大地電磁場二維正演模擬.石油地球物理勘探，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]黃一凡，胡祥云，韓波.自適應有限元法的大地電磁二維各向異性正反演.石油地球物理勘探，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.