多孔岩石中通常存在一种或多种流体，当地震波穿过含流体的多孔岩石时会引起孔隙中的流体流动，进而产生能量的损耗，引起速度频散和能量衰减。了解地震波在孔隙介质中传播的频散和衰减特征，对于烃类勘探、二氧化碳封存监测、水资源管理、地热能开发等至关重要^[1]。

流体流动可以分为宏观、介观和微观三种尺度。Biot^[2]建立了宏观尺度流体流动的地震波方程，为双相介质理论研究奠定了基础，但只考虑宏观尺度流体流动的模型^[3-4]无法解释地震勘探中实际观测到的波频散和衰减现象^[5]。Mavko和Nur^[6]提出的微观“喷射流”机制也不适用于地震频带。Dvorkin和Nur^[7]将Biot理论与喷射流机制相结合，提出了BISQ理论，在高频范围内能取得较好的效果，但仍无法解释地震频带内地震波的频散和衰减。近年来，研究表明介观尺度下的地震波诱导孔隙内流体流动是地震频带内波速度频散和衰减的主要原因^{[5, 8]}。White^[9]首次引入介观尺度的耗散机制并提出气液呈同心球体分布的球状斑块模型。Johnson^[10]对White模型进行了修正，推导出了含任意两种不同流体的多孔介质的动态体积模量表达式。此后研究人员提出了各种模型，以分析由孔隙结构^[11-12]和流体分布^[13-14]引起的波频散和衰减特性。

介观尺度的孔隙介质模型可以基本解释在地震频带内波的衰减和频散效应，但只考虑流体介观流动的单尺度模型具有一定的局限性。为解决单尺度模型的不足，许多研究者建立了双尺度模型和三尺度模型。双尺度模型具体可分为宏观与微观双尺度模型^{[7, 15]}、宏观与介观双尺度模型^{[8, 16-19]}、介观与微观双尺度模型^[20]三类。宏观与微观双尺度模型以BISQ模型^[7]为代表，此后有研究者研究了不同压力和裂缝条件下波在孔隙介质中的频散和衰减情况^[15]，结果显示宏观与微观双尺度模型描述波在高频时的频散和衰减较为准确。宏观与介观双尺度模型以双重孔隙介质模型^{[8, 11]}为代表，此后许多研究者进行了宏观与介观双尺度模型的研究^[16-19]，该种模型可以较准确地预测波在地震频带内的频散和衰减。介观与微观双尺度模型则是通过考虑孔隙介质中介观与微观两种尺度的裂缝，来描述波在孔隙介质内的频散和衰减情况^[20]。三尺度模型^[21-25]涵盖了三种尺度的流体流动，在全频带范围内都能取得较好的效果。但三尺度模型参数较为复杂，应用较为困难，并且目前的三尺度模型均未考虑束缚水的影响。

在油气藏形成过程中，由于油、水、气对岩石的润湿性差异和毛细管力的作用，运移的油气不可能把岩石孔隙中的水完全驱替出去，会有一定量的水残存在岩石孔隙中。由于特殊的分布和存在状态，这一部分水相对于岩石骨架几乎是不流动的，因而被称为束缚水。杨佳明等^[26]在Biot方程基础上采用黏弹性理论描述当地震波传播时束缚水与骨架间的耗散机制，但该模型只考虑了宏观尺度的流体流动，且模型内只含有一种流体。Zhang和He^[16]考虑到束缚水的影响，同时也考虑了宏观和介观尺度的流体流动，建立了部分饱和孔隙介质中的纵波方程。但该方程未考虑毛细管压力对介观尺度流体流动的影响，实际上毛细管压力会影响介观尺度流体斑块特性进而影响地震波的衰减^[27]。

为考虑毛细管压力对介观尺度流体流动的影响，本文将毛细管参数引入Johnson体积模量计算公式，并将改进后的体积模量计算公式代入部分饱和孔隙介质中的纵波方程，得到了部分饱和孔隙介质中的改进纵波方程。部分饱和孔隙介质中的改进纵波方程既考虑宏观和介观尺度的流体流动，也兼顾束缚水与毛细管压力的影响。为更加直观地表明地震波的传播情况，本文对部分饱和孔隙介质中的改进纵波方程进行有限差分求解，推导出该方程的高阶有限差分格式；采用Xie等^[28]提出的余弦型混合吸收边界，给出了边界区域的有限差分格式；最终得到了高质量的波场快照和叠后地震记录。通过对不同毛细管压力下的波场快照进行分析，可以得知毛细管压力对波速度的影响。对薄层状模型的叠后地震记录进行抽道分析，发现穿过部分饱和孔隙介质的地震波主频明显降低，说明本文正演算法正确，同时也说明部分饱和孔隙介质中的改进纵波方程能很好地反映流体流动导致的地震波衰减。对厚层状模型的叠后地震记录进行抽道分析，发现毛细管压力较大时，振幅衰减更大。

1 部分饱和孔隙介质中改进纵波方程的建立

Zhang和He^[16]基于White模型，在Biot方程的基础上建立了部分饱和孔隙介质中的纵波方程：

$ \rho \frac{{{\partial ^2}\theta }}{{\partial {t^2}}} + {\rho _{{f_2}}}\frac{{{\partial ^2}\varepsilon }}{{\partial {t^2}}} = H{\mathbf{∇} ^2}\theta + 2\gamma D{\mathbf{∇} ^2}\varepsilon 。$

(1)

$ {\rho _{{f_2}}}\frac{{{\partial ^2}\theta }}{{\partial {t^2}}} + m\frac{{{\partial ^2}\varepsilon }}{{\partial {t^2}}} = 2\gamma D{\mathbf{∇} ^2}\theta + 2D{\mathbf{∇} ^2}\varepsilon - \frac{{{\eta ^2}}}{\kappa }\frac{{\partial \varepsilon }}{{\partial t}}。$

(2)

其中:

$ \begin{aligned} & \rho=(1-\phi) \rho_s+\phi \mathrm{S}_1 \rho_{f_1}+\phi \mathrm{S}_2 \rho_{f_2}, \theta=\mathbf{∇} \cdot {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over u} }, \\ & \varepsilon=\mathbf{∇} \cdot {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over w} }, {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over w} }=\phi({\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over u} }-{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over U} }, H=\lambda+2 \mu=K+\frac{4}{3} \mu, \\ & \gamma=1-\frac{\bar{K}_m}{\bar{K}_s}, D=\frac{\bar{K}_s}{2}\left[\gamma+\frac{S_2 \phi}{K_{f_2}}\left(\bar{K}_s-K_{f_2}\right)\right]^{-1}, \\ & \bar{K}_s=K_b+ \\ & \frac{\left(1-K_b / K_s\right)^2}{\left(\phi S_1 /\left(1-\phi S_2\right)\right)\left(1 / K_{f_1}-1 / K_s\right)+\left(1 / K_s\right)\left(1-K_b / K_s\right)}, \\ & \bar{K}_m=\bar{K}_s\left(1-\left(\frac{\phi S_1}{\phi_c\left(1-\phi S_2\right)}\right)^n\right), m=s \frac{\rho_{f_2}}{\phi S_2} . \end{aligned} $

式中：ρ是孔隙岩石等效质量密度；ρ_s是岩石基质颗粒的质量密度；ρ_f1和ρ_f2分别是流体1和流体2的质量密度；ϕ是孔隙介质的孔隙度；S₁和S₂分别是流体1和流体2的饱和度；θ和ε分别是固相和流相的体应变； ${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over w} } $为流体2相对于固体的位移；${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over u} } $为含流体1的等效固相位移；${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over U} } $为流体2在绝对坐标系下的位移；n是计算固体骨架体积模量时的经验系数；m是动态耦合系数；s是与孔隙结构有关的常数；η₂是流体2的黏度；κ是岩石渗透率；t是时间；K_b和K_s分别是干燥岩石骨架和致密岩石基质的体积模量；K_f1和K_f2分别是流体1和流体2的体积模量；K_s是致密等效固体的体积模量；K_m是等效固体骨架的体积模量；H为部分孔隙介质的平面波模量；λ和μ为拉梅常数；K为部分饱和孔隙介质的体积模量；ϕ_C是临界孔隙度，对砂岩取0.4。

在上述方程中，部分饱和孔隙介质的体积模量K由下式给出：

$ K{\rm{ }} = {K_{{\rm{BGH}}}} - {\rm{ }}\frac{{{K_{{\rm{BGH}}}} - {K_{{\rm{BGW}}}}}}{{1 - \xi + \xi \sqrt {1 - i\omega \tau /{\xi ^2}} }}。$

(3)

式中：ω为地震波角频率；i为虚数单位；参数τ和ξ的计算方法详见文献[10]；K_BGH为部分饱和孔隙介质体积模量的高频极限，满足式(4)：

$ \frac{1}{{{K_{{\rm{BGH}}}} + \frac{4}{3}\mu }} = \frac{{{S_1}}}{{{K_{{\rm{BG}}}}\left( {{K_{{f_1}}}} \right) + \frac{4}{3}\mu }} + {\rm{ }}\frac{{{S_2}}}{{{K_{{\rm{BG}}}}\left( {{K_{{f_2}}}} \right) + \frac{4}{3}\mu }}。$

(4)

其中

$ {K_{{\rm{BG}}}}\left( {{K_i}} \right) = \frac{{{K_s} + \left[ {\phi {\rm{ }}\left( {{K_s}/{K_i}} \right) - \phi - 1} \right]{K_b}}}{{1{\rm{ }} - \phi - \left( {{K_b}/{K_s}} \right) + \phi {\rm{ }}\left( {{K_s}/{K_i}} \right)}},i = {f_1},{f_2}。$

(5)

式中：K_BGW为部分饱和孔隙介质体积模量的低频极限，${K_{{\rm{BGW}}}} = {K_{{\rm{BG}}}}\left( {{K_{\rm{W}}}} \right) $，K_w为两种混合流体的体积模量，将K_W代入式(5)即可得到K_BGW。在式(3)中，两种混合流体的体积模量K_W的计算采用的是Wood公式，即：

$ {K_{\rm{W}}} = \frac{{{K_{{f_1}}}\cdot{K_{{f_2}}}}}{{{S_1}{K_{{f_2}}} + {S_2}{K_{{f_1}}}}}。$

(6)

然而，Pagageorgiou等^[29]指出，当不相混溶的两种流体混合时，两种混合流体的等效体积模量会受毛细管压力的影响。假设流体1为润湿相，流体2为非润湿相，则混合流体的体积模量计算公式为：

$ K{'_{\rm{W}}} = {\rm{ }}\frac{{{S_1}{K_{{f_1}}} + \alpha {S_2}{K_{{f_2}}}}}{{{S_1} + \alpha {S_2}}},{\rm{ }}1 \le \alpha \le \frac{{{K_{{f_1}}}}}{{{K_{{f_2}}}}}。$

(7)

式(7)即为改进后的混合流体体积模量。其中α为毛细管压力参数，α越小毛细管压力越大，α=1时毛细管压力最大；$\alpha = \frac{{{K_{{f_1}}}}}{{{K_{{f_2}}}}} $时，可视为无毛细管压力影响，式(7)与式(6)相同。

将改进后的混合流体体积模量计算公式代入式(5)可得到考虑毛细管压力影响的部分饱和孔隙介质的体积模量低频极限K′_BGW，其表达式为：

$ \begin{array}{l} K{'_{{\rm{BGW}}}} = {K_{{\rm{BG}}}}K{'_W} = {K_{{\rm{BG}}}}\frac{{{S_1}{K_{{f_1}}} + \alpha {S_2}{K_{{f_2}}}}}{{{S_1} + \alpha {S_2}}} = \\ \frac{{{K_s} + \left[ {\phi \left( {{K_s}\frac{{\left( {{S_1} + \alpha {S_2}} \right)}}{{{S_1}{K_{{f_1}}} + \alpha {S_2}{K_{{f_2}}}}}} \right) - \phi - 1} \right]{K_b}}}{{1 - \phi - \left( {{K_b}/{K_s}} \right) + \phi \left( {{K_s}\frac{{\left( {{S_1} + \alpha {S_2}} \right)}}{{{S_1}{K_{{f_1}}} + \alpha {S_2}{K_{{f_2}}}}}} \right)}} \end{array}。$

(8)

将K′_BGW和K′_W代入式(3)(参数τ和ξ中也涉及K_BGW和K_W，详见附录A)，得到考虑毛细管压力影响的部分饱和孔隙介质的体积模量K′：

$ K{\rm{ }}' = {K_{{\rm{BGH}}}} - {\rm{ }}\frac{{{K_{{\rm{BGH}}}} - K{'_{{\rm{BGW}}}}}}{{1 - \xi ' + \xi {\rm{ }}'\sqrt {1 - i\omega \tau '/\xi {'^2}} }}。$

(9)

将式(9)代入方程(1)，即可得到考虑毛细管压力影响的部分饱和孔隙介质改进纵波方程：

$ \rho \frac{{{\partial ^2}\theta }}{{\partial {t^2}}} + \rho {f_2}\frac{{{\partial ^2}\varepsilon }}{{\partial {t^2}}} = H{\mathbf{∇} ^2}\theta + 2\gamma D{\mathbf{∇} ^2}\varepsilon 。$

(10)

$ {\rho _{{f_2}}}\frac{{{\partial ^2}\theta }}{{\partial {t^2}}} + m\frac{{{\partial ^2}\varepsilon }}{{\partial {t^2}}} = 2\gamma D{\mathbf{∇} ^2}\theta + 2D{\mathbf{∇} ^2}\varepsilon - \frac{{{\eta ^2}}}{\kappa }\frac{{\partial \varepsilon }}{{\partial t}}。$

(11)

其中$H' = K' + \frac{4}{3}\mu $。

Tisato and Quintal^[30]通过实验测得了1~100 Hz内含两种流体的孔隙介质中快纵波的衰减量，孔隙介质物性参数如表 1所示。为探究毛细管压力对部分饱和孔隙介质体积模量的影响，本文采用表 1参数绘制了不同毛细管压力作用下K′随频率变化的曲线。图 1中深红色、绿色、黄色和蓝色线分别代表$\alpha = 1{\rm{ }}.{\rm{ }}0 \times 1{\rm{ }}{0^4} $、$\alpha = 1{\rm{ }}.{\rm{ }}5 \times 1{\rm{ }}{0^4} $、$\alpha = 1{\rm{ }}.{\rm{ }}8 \times 1{\rm{ }}{0^4} $和$\alpha = 2{\rm{ }}.{\rm{ }}2 \times 1{\rm{ }}{0^4} $时部分饱和孔隙介质体积模量随频率变化的曲线。观察图 1(a)可知部分饱和孔隙介质体积模量实部随频率的增大而增大，且毛细管压力越大(α越小)，部分饱和孔隙介质体积模量实部越大。观察图 1(b)可以发现部分饱和孔隙介质体积模量虚部随频率的增大而减小，且毛细管压力越大，部分饱和孔隙介质体积模量虚部越小。

为验证部分饱和孔隙介质中改进纵波方程的正确性，本文采用平面波分析法求解该方程得到了快纵波逆品质因子的表达式。根据表 1参数，分别计算了改进后方程(α=1.8×10⁴)与White模型、Johnson模型、改进前方程预测的快纵波衰减，并将其与实验测得的衰减量进行对比(见图 2)。从图 2可以看出，改进后方程预测的衰减与实验室测量数据最为吻合，验证了部分饱和孔隙介质改进纵波方程的正确性。

2 部分饱和孔隙介质中改进纵波方程的有限差分解法

在验证部分饱和孔隙介质改进纵波方程正确性后，对其进行正演模拟，推导其时间二阶、空间高阶的差分格式。

首先对空间模型网格化，将网格进行等步长划分，x方向和z方向的空间步长分别为Δx和Δz，时间采样步长为Δt。$\varepsilon^n(i, j) $ 和$\theta^n(i, j) $ 分别表示(i, j)点在n时刻的流相体应变和固相体应变的波场值。考虑二维情况，在复数域求解，推导时间二阶、空间2N阶的高阶有限差分解。

令$\theta = {\rm{ }}{\theta _{re}} + i{\theta _{im}},\varepsilon = {\rm{ }}{\varepsilon _{re}} + i{\varepsilon _{im}},H' = H{'_{re}} + iH{'_{im}} $，i为虚数单位，则式(10)和式(11)可写为：

$ \begin{array}{l} \rho \frac{{{\partial ^2}{\theta _{{\rm{re}}}}}}{{\partial {t^2}}} + {\rm{i}}\rho \frac{{{\partial ^2}{\theta _{{\rm{im}}}}}}{{\partial {t^2}}}{\rm{ }} + {\rho _{{f_2}}}\frac{{{\partial ^2}{\varepsilon _{{\rm{re}}}}}}{{\partial {t^2}{\rm{ }}}}{\rm{ }} + {\rm{i}}{\rho _{{f_2}}}\frac{{{\partial ^2}{\varepsilon _{{\rm{im}}}}}}{{\partial {t^2}}}{\rm{ }} = {\rm{ }}\left( {H{'_{{\rm{re}}}} + {\rm{i}}H{'_{{\rm{im}}}}} \right)\cdot\\ \left( {\frac{{{\partial ^2}{\theta _{{\rm{re}}}}}}{{\partial {x^2}}} + \frac{{{\partial ^2}{\theta _{{\rm{re}}}}}}{{\partial {z^2}}}{\rm{ }} + {\rm{i}}\frac{{{\partial ^2}{\theta _{{\rm{im}}}}}}{{\partial {x^2}}}{\rm{ }} + {\rm{i}}\frac{{{\partial ^2}{\theta _{{\rm{im}}}}}}{{\partial {z^2}}}{\rm{ }}} \right) + 2\gamma D\left( {\frac{{{\partial ^2}{\varepsilon _{{\rm{re}}}}}}{{\partial {x^2}}}{\rm{ }} + \frac{{{\partial ^2}{\varepsilon _{{\rm{ re}}}}}}{{\partial {z^2}}}{\rm{ }} + {\rm{i}}\frac{{{\partial ^2}{\varepsilon _{{\rm{im}}}}}}{{\partial {x^2}}}{\rm{ }} + {\rm{i}}\frac{{{\partial ^2}{\varepsilon _{{\rm{im}}}}}}{{\partial {z^2}}}{\rm{ }}} \right) \end{array}。$

(12)

$ \begin{array}{l} {\rho _{{f_2}}}\frac{{{\partial ^2}{\theta _{{\rm{re}}}}}}{{\partial {t^2}}} + {\rm{i}}{\rho _{{f_2}}}\frac{{{\partial ^2}{\theta _{{\rm{im}}}}}}{{\partial {t^2}}} + m\frac{{{\partial ^2}{\varepsilon _{{\rm{re}}}}}}{{\partial {t^2}}} + {\rm{i}}m\frac{{{\partial ^2}{\varepsilon _{{\rm{im}}}}}}{{\partial {t^2}}} = 2\gamma D\left( {\frac{{{\partial ^2}{\theta _{{\rm{re}}}}}}{{\partial {x^2}}} + \frac{{{\partial ^2}{\theta _{{\rm{re}}}}}}{{\partial {z^2}}} + {\rm{i}}\frac{{{\partial ^2}{\theta _{{\rm{im}}}}}}{{\partial {x^2}}} + {\rm{i}}\frac{{{\partial ^2}{\theta _{{\rm{im}}}}}}{{\partial {z^2}}}} \right) + \\ 2D\left( {\frac{{{\partial ^2}{\varepsilon _{{\rm{re}}}}}}{{\partial {x^2}}} + \frac{{{\partial ^2}{\varepsilon _{{\rm{re}}}}}}{{\partial {z^2}}} + {\rm{i}}\frac{{{\partial ^2}{\varepsilon _{{\rm{im}}}}}}{{\partial {x^2}}} + {\rm{i}}\frac{{{\partial ^2}{\varepsilon _{{\rm{im}}}}}}{{\partial {z^2}}}} \right) \end{array}。$

(13)

将式(12)和式(13)进行整理，并将虚部和实部进行分离可得式(14)—(17)：

$ \rho \frac{{{\partial ^2}{\theta _{{\rm{re}}}}}}{{\partial {t^2}}} + {\rho _{{f_2}}}\frac{{{\partial ^2}{\varepsilon _{{\rm{re}}}}}}{{\partial {t^2}{\rm{ }}}} = H{'_{{\rm{re}}}}\left( {\frac{{{\partial ^2}{\theta _{{\rm{re}}}}}}{{\partial {x^2}}} + \frac{{{\partial ^2}{\theta _{{\rm{re}}}}}}{{\partial {z^2}}}} \right) - {\rm{i}}H{'_{{\rm{im}}}}\left( {\frac{{{\partial ^2}{\theta _{{\rm{im}}}}}}{{\partial {x^2}}} + \frac{{{\partial ^2}{\theta _{{\rm{im}}}}}}{{\partial {z^2}}}} \right) + 2\gamma D\left( {\frac{{{\partial ^2}{\varepsilon _{{\rm{re}}}}}}{{\partial {x^2}}}{\rm{ }} + \frac{{{\partial ^2}{\varepsilon _{{\rm{ re}}}}}}{{\partial {z^2}}}} \right) $

(14)

$ \rho \frac{{{\partial ^2}{\theta _{{\rm{im}}}}}}{{\partial {t^2}}} + {\rho _{{f_2}}}\frac{{{\partial ^2}{\varepsilon _{{\rm{im}}}}}}{{\partial {t^2}}} = H{'_{{\rm{im}}}}\left( {\frac{{{\partial ^2}{\theta _{{\rm{re}}}}}}{{\partial {x^2}}} + \frac{{{\partial ^2}{\theta _{{\rm{re}}}}}}{{\partial {z^2}}}} \right) + H{'_{{\rm{re}}}}\left( {\frac{{{\partial ^2}{\theta _{{\rm{im}}}}}}{{\partial {x^2}}} + \frac{{{\partial ^2}{\theta _{{\rm{im}}}}}}{{\partial {z^2}}}} \right) + 2\gamma D\left( {\frac{{{\partial ^2}{\varepsilon _{{\rm{im}}}}}}{{\partial {x^2}}} + \frac{{{\partial ^2}{\varepsilon _{{\rm{im}}}}}}{{\partial {z^2}}}} \right) $

(15)

$ {\rho _{{f_2}}}\frac{{{\partial ^2}{\theta _{{\rm{re}}}}}}{{\partial {t^2}}} + m\frac{{{\partial ^2}{\varepsilon _{{\rm{re}}}}}}{{\partial {t^2}{\rm{ }}}} = 2\gamma D\left( {\frac{{{\partial ^2}{\theta _{{\rm{re}}}}}}{{\partial {x^2}}} + \frac{{{\partial ^2}{\theta _{{\rm{re}}}}}}{{\partial {z^2}}}} \right) + 2D\left( {\frac{{{\partial ^2}{\varepsilon _{{\rm{re}}}}}}{{\partial {x^2}}}{\rm{ }} + \frac{{{\partial ^2}{\varepsilon _{{\rm{ re}}}}}}{{\partial {z^2}}}} \right) - \frac{{{\eta _2}}}{\kappa }\frac{{\partial {\varepsilon _{{\rm{re}}}}}}{{\partial t}} $

(16)

$ {\rho _{{f_2}}}\frac{{{\partial ^2}{\theta _{{\rm{im}}}}}}{{\partial {t^2}}} + m\frac{{{\partial ^2}{\varepsilon _{{\rm{im}}}}}}{{\partial {t^2}}} = 2\gamma D\left( {\frac{{{\partial ^2}{\theta _{{\rm{im}}}}}}{{\partial {x^2}}} + \frac{{{\partial ^2}{\theta _{{\rm{im}}}}}}{{\partial {z^2}}}} \right) + 2D\left( {\frac{{{\partial ^2}{\varepsilon _{{\rm{im}}}}}}{{\partial {x^2}}} + \frac{{{\partial ^2}{\varepsilon _{{\rm{im}}}}}}{{\partial {z^2}}}} \right) - \frac{{{\eta _2}}}{\kappa }\frac{{\partial {\varepsilon _{{\rm{im}}}}}}{{\partial t}} $

(17)

由于：

$ \begin{array}{*{20}{l}} {\frac{{{\partial ^2}\theta }}{{\partial {t^2}}} = \frac{{{\theta ^{n + 1}}(i,j) + {\theta ^{n - 1}}(i,j) - 2{\theta ^n}(i,j)}}{{\Delta {t^2}}},}\\ {\frac{{\partial \varepsilon }}{{\partial t}} = \frac{{{\varepsilon ^{n + 1}}(i,j) - {\varepsilon ^{n - 1}}(i,j)}}{{2\Delta t}},\frac{{{\partial ^2}\varepsilon }}{{\partial {t^2}}} = \frac{{{\varepsilon ^{n + 1}}(i,j) + {\varepsilon ^{n - 1}}(i,j) - 2{\varepsilon ^n}(i,j)}}{{\Delta {t^2}}},}\\ {{\mathbf{∇} ^2}\theta = \frac{{{\partial ^2}\theta }}{{\partial {x^2}}} + \frac{{{\partial ^2}\theta }}{{\partial {z^2}}} = \sum\limits_{m = - N}^N {\frac{{{c_m}^{(N)}{\theta ^n}(i + m,j)}}{{\Delta {x^2}}}} {\rm{ + }}\sum\limits_{m = - N}^N {\frac{{{c_m}^{(N)}{\theta ^n}(i,j + m)}}{{\Delta {z^2}}},} }\\ {{\mathbf{∇} ^2}\varepsilon = \frac{{{\partial ^2}\varepsilon }}{{\partial {x^2}}} + \frac{{{\partial ^2}\varepsilon }}{{\partial {z^2}}} = \sum\limits_{m = - N}^N {\frac{{{c_m}^{(N)}{\varepsilon ^n}(i + m,j)}}{{\Delta {x^2}}}} {\rm{ + }}\sum\limits_{m = - N}^N {\frac{{{c_m}^{(N)}{\varepsilon ^n}(i,j + m)}}{{\Delta {z^2}}}。} } \end{array} $

故式(14)—(17)可变为：

$ \begin{array}{l} \rho \frac{{{\theta _{{\rm{re}}}}^{n + 1}(i,j) + {\theta _{{\rm{re}}}}^{n - 1}(i,j) - 2{\theta _{{\rm{re}}}}^n(i,j)}}{{\Delta {t^2}}} + {\rho _{{f_2}}}\frac{{{\varepsilon _{{\rm{re}}}}^{n + 1}(i,j) + {\varepsilon _{{\rm{re}}}}^{n - 1}(i,j) - 2{\varepsilon _{{\rm{re}}}}^n(i,j)}}{{\Delta {t^2}}} = \\ H{'_{{\rm{re}}}}\left( {\sum\limits_{m = - N}^N {\frac{{{\rm{ }}{c_m}^{(N)}{\theta _{{\rm{re}}}}^n(i{\rm{ }} + m,j)}}{{\Delta {x^2}}}} {\rm{ + }}\sum\limits_{m = - N}^N {\frac{{{c_m}^{(N)}{\theta _{{\rm{re}}}}^n(i{\rm{ }},j + m)}}{{\Delta {z^2}}}} } \right)\\ H{'_{{\rm{im}}}}\left( {\sum\limits_{m = - N}^N {\frac{{{\rm{ }}{c_m}^{(N)}{\theta _{{\rm{im}}}}^n(i{\rm{ }} + m,j)}}{{\Delta {x^2}}}} {\rm{ + }}\sum\limits_{m = - N}^N {\frac{{{c_m}^{(N)}{\theta _{{\rm{im}}}}^n(i{\rm{ }},j + m)}}{{\Delta {z^2}}}} } \right)\\ 2\gamma D\left( {\sum\limits_{m = - N}^N {\frac{{{\rm{ }}{c_m}^{(N)}{\varepsilon _{{\rm{re}}}}^n(i{\rm{ }} + m,j)}}{{\Delta {x^2}}}} {\rm{ + }}\sum\limits_{m = - N}^N {\frac{{{c_m}^{(N)}{\varepsilon _{{\rm{re}}}}^n(i{\rm{ }},j + m)}}{{\Delta {z^2}}}} } \right) \end{array}。$

(18)

$ \begin{array}{l} \rho \frac{{{\theta _{{\rm{im}}}}^{n + 1}(i,j) + {\theta _{{\rm{im}}}}^{n - 1}(i,j) - 2{\theta _{{\rm{im}}}}^n(i,j)}}{{\Delta {t^2}}} + {\rho _{{f_2}}}\frac{{{\varepsilon _{{\rm{im}}}}^{n + 1}(i,j) + {\varepsilon _{{\rm{im}}}}^{n - 1}(i,j) - 2{\varepsilon _{{\rm{im}}}}^n(i,j)}}{{\Delta {t^2}}} = \\ {\rm{ }}H{'_{{\rm{im}}}}\left( {\sum\limits_{m = - N}^N {\frac{{{\rm{ }}{c_m}^{(N)}{\theta _{{\rm{re}}}}^n(i{\rm{ }} + m,j)}}{{\Delta {x^2}}}} {\rm{ + }}\sum\limits_{m = - N}^N {\frac{{{c_m}^{(N)}{\theta _{{\rm{re}}}}^n(i{\rm{ }},j + m)}}{{\Delta {z^2}}}} } \right)\\ H{'_{{\rm{re}}}}\left( {\sum\limits_{m = - N}^N {\frac{{{\rm{ }}{c_m}^{(N)}{\theta _{{\rm{im}}}}^n(i{\rm{ }} + m,j)}}{{\Delta {x^2}}}} {\rm{ + }}\sum\limits_{m = - N}^N {\frac{{{c_m}^{(N)}{\theta _{{\rm{im}}}}^n(i{\rm{ }},j + m)}}{{\Delta {z^2}}}} } \right)\\ 2\gamma D\left( {\sum\limits_{m = - N}^N {\frac{{{\rm{ }}{c_m}^{(N)}{\varepsilon _{{\rm{im}}}}^n(i{\rm{ }} + m,j)}}{{\Delta {x^2}}}} {\rm{ + }}\sum\limits_{m = - N}^N {\frac{{{c_m}^{(N)}{\varepsilon _{{\rm{im}}}}^n(i{\rm{ }},j + m)}}{{\Delta {z^2}}}} } \right) \end{array}。$

(19)

$ \begin{array}{l} {\rho _{{f_2}}}\frac{{{\theta _{{\rm{re}}}}^{n + 1}(i,j) + {\theta _{{\rm{re}}}}^{n - 1}(i,j) - 2{\theta _{{\rm{re}}}}^n(i,j)}}{{\Delta {t^2}}} + m\frac{{{\varepsilon _{{\rm{re}}}}^{n + 1}(i,j) + {\varepsilon _{{\rm{re}}}}^{n - 1}(i,j) - 2{\varepsilon _{{\rm{re}}}}^n(i,j)}}{{\Delta {t^2}}} = \\ 2\gamma D\left( {\sum\limits_{m = - N}^N {\frac{{{\rm{ }}{c_m}^{(N)}{\theta _{{\rm{re}}}}^n(i{\rm{ }} + m,j)}}{{\Delta {x^2}}}} {\rm{ + }}\sum\limits_{m = - N}^N {\frac{{{c_m}^{(N)}{\theta _{{\rm{re}}}}^n(i{\rm{ }},j + m)}}{{\Delta {z^2}}}} } \right) + \\ 2D\left( {\sum\limits_{m = - N}^N {\frac{{{\rm{ }}{c_m}^{(N)}{\varepsilon _{{\rm{re}}}}^n(i{\rm{ }} + m,j)}}{{\Delta {x^2}}}} {\rm{ + }}\sum\limits_{m = - N}^N {\frac{{{c_m}^{(N)}{\varepsilon _{{\rm{re}}}}^n(i{\rm{ }},j + m)}}{{\Delta {z^2}}}} } \right) - \frac{{{\eta _2}}}{\kappa }\frac{{{\varepsilon _{{\rm{re}}}}^{n + 1}(i,j) - {\varepsilon _{re}}^{n - 1}(i,j)}}{{2\Delta t}} \end{array}。$

(20)

$ \begin{array}{l} {\rho _{{f_2}}}\frac{{{\theta _{{\rm{im}}}}^{n + 1}(i,j) + {\theta _{{\rm{im}}}}^{n - 1}(i,j) - 2{\theta _{{\rm{im}}}}^n(i,j)}}{{\Delta {t^2}}} + m\frac{{{\varepsilon _{{\rm{im}}}}^{n + 1}(i,j) + {\varepsilon _{{\rm{im}}}}^{n - 1}(i,j) - 2{\varepsilon _{{\rm{im}}}}^n(i,j)}}{{\Delta {t^2}}} = \\ 2\gamma D\left( {\sum\limits_{m = - N}^N {\frac{{{\rm{ }}{c_m}^{(N)}{\theta _{{\rm{im}}}}^n(i{\rm{ }} + m,j)}}{{\Delta {x^2}}}} {\rm{ + }}\sum\limits_{m = - N}^N {\frac{{{c_m}^{(N)}{\theta _{{\rm{im}}}}^n(i{\rm{ }},j + m)}}{{\Delta {z^2}}}} } \right) + \\ 2D\left( {\sum\limits_{m = - N}^N {\frac{{{\rm{ }}{c_m}^{(N)}{\varepsilon _{{\rm{im}}}}^n(i{\rm{ }} + m,j)}}{{\Delta {x^2}}}} {\rm{ + }}\sum\limits_{m = - N}^N {\frac{{{c_m}^{(N)}{\varepsilon _{{\rm{im}}}}^n(i{\rm{ }},j + m)}}{{\Delta {z^2}}}} } \right) - \frac{{{\eta _2}}}{\kappa }\frac{{{\varepsilon _{{\rm{im}}}}^{n + 1}(i,j) - {\varepsilon _{{\rm{im}}}}^{n - 1}(i,j)}}{{2\Delta t}} \end{array}。$

(21)

求解式(18)—(21)可得：

$ \left\{ {\begin{array}{*{20}{c}} {{\theta _{{\rm{re}}}}{{^n}^{ + 1}}(i,j) = \frac{{{\rho _{{f_2}}}B - (m + \frac{{{\eta _2}}}{{2\kappa }}\Delta t)A}}{{{\rho ^2}_{{f_2}} - m\rho - \frac{{{\eta _2}}}{{2\kappa }}\rho \Delta t}}}\\ {{\theta _{{\rm{im}}}}{{^n}^{ + 1}}(i,j) = \frac{{{\rho _{{f_2}}}{B^*} - (m + \frac{{{\eta _2}}}{{2\kappa }}\Delta t){A^*}}}{{{\rho ^2}_{{f_2}} - m\rho - \frac{{{\eta _2}}}{{2\kappa }}\rho \Delta t}}}\\ {{\varepsilon _{{\rm{re}}}}{{^n}^{ + 1}}(i,j) = \frac{{\rho B - {\rho _{{f_2}}}A}}{{m\rho + \frac{{{\eta _2}}}{{2\kappa }}\rho \Delta t - {\rho ^2}_{{f_2}}}}}\\ {{\varepsilon _{{\rm{im}}}}{{^n}^{ + 1}}(i,j) = \frac{{\rho {B^*} - {\rho _{{f_2}}}{A^*}}}{{m\rho + \frac{{{\eta _2}}}{{2\kappa }}\rho \Delta t - {\rho ^2}_{{f_2}}}}} \end{array}} \right. $

(22)

其中：

$ \begin{aligned} & A=\left[H_{\mathrm{re}}^{\prime}\left(\sum\limits_{m=-\mathrm{N}}^N \frac{c_m^{(\mathrm{N})} \theta_{\mathrm{re}}^n(i+m, j)}{\Delta x^2}+\sum\limits_{m=-\mathrm{N}}^N \frac{c_m^{(\mathrm{N})} \theta_{\mathrm{re}}^n(i, j+m)}{\Delta z^2}\right)-H_{\mathrm{im}}^{\prime}\left(\sum\limits_{m=-\mathrm{N}}^N \frac{c_m^{(\mathrm{N})} \theta_{\mathrm{im}}^n(i+m, j)}{\Delta x^2}+\sum\limits_{m=-\mathrm{N}}^N \frac{c_m^{(\mathrm{N})} \theta_{\mathrm{im}}^n(i, j+m)}{\Delta z^2}\right)+\right. \\ & \left.2 \gamma D\left(\sum\limits_{m=-\mathrm{N}}^N \frac{c_m^{(\mathrm{N})} \varepsilon_{\mathrm{re}}^n(i+m, j)}{\Delta x^2}+\sum\limits_{m=-\mathrm{N}}^N \frac{c_m^{(\mathrm{N})} \varepsilon_{\mathrm{re}}^n(i, j+m)}{\Delta z^2}\right)\right] \Delta t^2-\rho\left(\theta_{\mathrm{re}}^{n-1}(i, j)-2 \theta_{\mathrm{re}}^n(i, j)\right)-\rho_{f_2}\left(\varepsilon_{\mathrm{re}}^{n-1}(i, j)-2 \varepsilon_{\mathrm{re}}^n(i, j)\right), \\ & A^*=\left[H_{\mathrm{im}}^{\prime}\left(\sum\limits_{m=-\mathrm{N}}^N \frac{c_m^{(N)} \theta_{\mathrm{re}}^n(i+m, j)}{\Delta x^2}+\sum\limits_{m=-\mathrm{N}}^N \frac{c_m^{(\mathrm{N})} \theta_{\mathrm{re}}^n(i, j+m)}{\Delta z^2}\right)+H_{\mathrm{re}}^{\prime}\left(\sum\limits_{m=-\mathrm{N}}^N \frac{c_m^{(\mathrm{N})} \theta_{\mathrm{im}}^n(i+m, j)}{\Delta x^2}+\sum\limits_{m=-\mathrm{N}}^N \frac{c_m^{(N)} \theta_{\mathrm{im}}^n(i, j+m)}{\Delta z^2}\right)+\right. \\ & \left.2 \gamma D\left(\sum\limits_{m=-\mathrm{N}}^N \frac{c_m^{(N)} \varepsilon_{\mathrm{im}}^n(i+m, j)}{\Delta x^2}+\sum\limits_{m=-\mathrm{N}}^N \frac{c_m^{(N)} \varepsilon_{\mathrm{im}}^n(i, j+m)}{\Delta z^2}\right)\right] \Delta t^2-\rho\left(\theta_{\mathrm{im}}^{n-1}(i, j)-2 \theta_{\mathrm{im}}^n(i, j)\right)-\rho_{f_2}\left(\varepsilon_{\mathrm{im}}^{n-1}(i, j)-2 \varepsilon_{\mathrm{im}}^n(i, j)\right), \\ & B=\left[2 \gamma D\left(\sum\limits_{m=-\mathrm{N}}^N \frac{c_m^{(N)} \theta_{\mathrm{re}}^n(i+m, j)}{\Delta x^2}+\sum\limits_{m=-\mathrm{N}}^N \frac{c_m^{(\mathrm{N})} \theta_{\mathrm{re}}^n(i, j+m)}{\Delta z^2}\right)+2 D\left(\sum\limits_{m=-\mathrm{N}}^N \frac{c_m^{(\mathrm{N})} \varepsilon_{\mathrm{re}}^n(i+m, j)}{\Delta x^2}+\sum\limits_{m=-\mathrm{N}}^N \frac{c_m^{(N)} \varepsilon_{\mathrm{re}}^n(i, j+m)}{\Delta z^2}\right)\right] \Delta t^2+ \\ & \frac{\eta_2}{2 \kappa_{\mathrm{re}}^{n-1}}(i, j) \Delta t-\rho_{f_2}\left(\theta_{\mathrm{re}}^{n-1}(i, j)-2 \theta_{\mathrm{re}}^n(i, j)\right)-m\left(\varepsilon_{\mathrm{re}}^{n-1}(i, j)-2 \varepsilon_{\mathrm{re}}^n(i, j)\right), \\ & B^*=\left[2 \gamma D\left(\sum\limits_{m=-\mathrm{N}}^N \frac{c_m^{(N)} \theta_{\mathrm{im}}^n(i+m, j)}{\Delta x^2}+\sum\limits_{m=-\mathrm{N}}^N \frac{c_m^{(N)} \theta_{\mathrm{im}}^n(i, j+m)}{\Delta z^2}\right)+2 D\left(\sum\limits_{m=-\mathrm{N}}^N \frac{c_m^{(N)} \varepsilon_{\mathrm{im}}^n(i+m, j)}{\Delta x^2}+\sum\limits_{m=-\mathrm{N}}^N \frac{c_m^{(N)} \varepsilon_{\mathrm{im}}^n(i, j+m)}{\Delta z^2}\right)\right] \Delta t^2+ \\ & \frac{\eta_2}{2 \kappa_{\mathrm{im}}^{n-1}}(i, j) \Delta t-\rho_{f_2}\left(\theta_{\mathrm{im}}^{n-1}(i, j)-2 \theta_{\mathrm{im}}^n(i, j)\right)-m\left(\varepsilon_{\mathrm{im}}^{n-1}(i, j)-2 \varepsilon_{\mathrm{im}}^n(i, j)\right) 。\end{aligned} $

式(22)即为部分饱和孔隙介质时间二阶、空间高阶的差分格式，其中$\theta_{\mathrm{re}}^{n+1}(i, j) $ 和$\theta_{\mathrm{im}}^{n+1}(i, j) $ 分别为固相在(i, j)点n+1时刻体应变的实部和虚部；$\varepsilon _{{\rm{re}}}^{n + 1}(i,j) $ 和$\varepsilon_{\mathrm{im}}^{n+1}(i, j) $ 分别为流相在(i, j)点n+1时刻体应变的实部和虚部。$c_m^{(N)} $ 为差分系数，$c_m^{(N)}=c_{-m}^{(N)} $ ，差分系数具体值见文献[31]。

在使用有限差分技术进行数值模拟时，需要满足数值模拟的稳定性要求。根据规则网格的稳定性条件计算公式^[30]，满足本文正演模拟需求的时间2阶、空间2N阶的稳定性条件为：

$ {v_{{\rm{max}}}}\Delta t\sqrt {\frac{1}{{{{\left( {\Delta x} \right)}^2}}} + \frac{1}{{{{\left( {\Delta z} \right)}^2}}}} \le {\left( {\frac{1}{{\sum\limits_{m = 1}^{{N_1}} {c_{_{2m - 1}}^{^{(N)}}} }}} \right)^{1/2}}。$

(23)

式中：N₁为不超过N的最大奇数；v_max代表模型的速度最大值，当v_max为含流体储层中第一纵波速度时，其表达式见文献[16]。

3 吸收边界条件

目前的人工边界处理方法主要包括完美匹配层(Perfectly matched layer, PML)方法和混合型吸收边界条件(Absorbing boundary conditions, ABCs)方法。PML边界通常需要在边界上使用几十甚至上百个匹配层才能实现更好的吸收效果，这会产生巨大的计算和内存负担。2017年，Liu等^[32]将混合型ABCs边界与传统分裂PML进行比较，发现混合型ABCs边界可以在更短的运算时间内实现更好的吸收效果。2020年，Xie等^[28]将传统混合型ABCs边界的过渡区域的线性加权模式改为余弦型，使得中心波场方程和边界方程更好地耦合。本文采用Xie等提出的余弦型混合吸收边界，吸收部分饱和孔隙介质中边界处的伪反射。

混合ABCs边界如图 3所示，白色区域为中心波场区域，其差分格式如式(22)所示。灰色区域为混合吸收边界区域，共M层。混合吸收边界可被图 3中蓝线划分为8个边界区域，本文直接给出左边界和左上边界的吸收边界条件及其差分格式，其余边界条件可以类推，具体推导过程可见文献[28]。

混合边界左边界的吸收边界条件为

$ \frac{{{\partial ^2}\varphi }}{{\partial {x^2}}} + \frac{{cos{\alpha _1}cos{\alpha _2}}}{{{v^2}}}\frac{{{\partial ^2}\varphi }}{{\partial {t^2}}} - \frac{{\left( {cos{\alpha _1} + cos{\alpha _2}} \right)}}{v}{\rm{ }}\frac{{{\partial ^2}\varphi {\rm{ }}}}{{\partial t\partial x}}{\rm{ }} = 0 $

(24)

混合边界左边界差分格式为

$ \begin{array}{l} {\varphi _l} = {\varphi ^{n + 1}}(i,j) = (1 - {\sigma _l})[2{\varphi ^n}(i,j) - {\varphi ^{n - 1}}(i,j) + \\ {v^2}\Delta {t^2}\sum\limits_{n = 1}^N {c_m^{(N)}} {\rm{ }}(\frac{{{\varphi ^n}(i - m,j) + {\varphi ^n}(i + m,j)}}{{\Delta {x^2}}} + \frac{{{\varphi ^n}(i,j - m){\rm{ }} + {\rm{ }}{\varphi ^n}(i,j + m)}}{{\Delta {z^2}}}) + {v^2}\Delta {t^2}c_0^{(N)}{\rm{ }}(\frac{{{\varphi ^n}(i{\rm{ }},j)}}{{\Delta {x^2}}}{\rm{ }} + \frac{{{\varphi ^n}(i{\rm{ }},j)}}{{\Delta {z^2}}})] + \\ {\sigma _l}[({\rm{ }}2 - {\lambda _x} - {\mu _x}){\varphi ^n}(i,j) + ({\lambda _x} + 2{\mu _x}){\varphi ^n}(i + 1{\rm{ }},j) - {\mu _x}{\varphi ^n}(i + 2{\rm{ }},j) + ({\lambda _x} - 1{\rm{ }}){\varphi ^{n - 1}}(i,j) - {\lambda _x}{\varphi ^{n - 1}}(i + 1{\rm{ }},j)] \end{array}。$

(25)

混合边界左上边界的吸收边界条件为：

$ (\frac{{{\partial ^2}\varphi }}{{\partial {x^2}}} + \frac{{{\partial ^2}\varphi }}{{\partial {z_2}}}{\rm{ }}) + \frac{{2cos{\alpha _1}cos{\alpha _2}}}{{{v^2}}}\frac{{{\partial ^2}\varphi }}{{\partial {t^2}}} - \frac{{\left( {cos{\alpha _1} + cos{\alpha _2}} \right)}}{v}(\frac{{{\partial ^2}\varphi }}{{\partial t\partial x}}{\rm{ }} + \frac{{{\partial ^2}\varphi }}{{\partial t\partial z}}{\rm{ }}) = 0。$

(26)

混合边界左上边界差分格式为:

$ \begin{array}{l} {\varphi _l} = {\varphi ^{n + 1}}(i,j) = (1 - {\sigma _l})[2{\varphi ^n}(i,j) - {\varphi ^{n - 1}}(i,j) + \\ {v^2}\Delta {t^2}\sum\limits_{n = 1}^N {c_m^{(N)}} {\rm{ }}(\frac{{{\varphi ^n}(i - m,j) + {\varphi ^n}(i + m,j)}}{{\Delta {x^2}}} + \frac{{{\varphi ^n}(i,j - m){\rm{ }} + {\rm{ }}{\varphi ^n}(i,j + m)}}{{\Delta {z^2}}}) + {v^2}\Delta {t^2}c_0^{(N)}{\rm{ }}(\frac{{{\varphi ^n}(i{\rm{ }},j)}}{{\Delta {x^2}}}{\rm{ }} + \frac{{{\varphi ^n}(i{\rm{ }},j)}}{{\Delta {z^2}}})] + \\ {\sigma _l}\{ 0.5[({\rm{ }}2 - {\lambda _x} - {\mu _x}){\varphi ^n}(i,j) + ({\lambda _x} + 2{\mu _x}){\varphi ^n}(i + 1{\rm{ }},j) - {\mu _x}{\varphi ^n}(i + 2{\rm{ }},j) + ({\lambda _x} - 1{\rm{ }}){\varphi ^{n - 1}}(i,j) - {\lambda _x}{\varphi ^{n - 1}}(i + 1{\rm{ }},j)] + \\ 0.5[({\rm{ }}2 - {\lambda _z} - {\mu _z}){\varphi ^n}(i,j) + ({\lambda _z} + 2{\mu _z}){\varphi ^n}(i{\rm{ }},j + 1) - {\mu _z}{\varphi ^n}(i{\rm{ }},j + 2) + ({\lambda _z} - 1{\rm{ }}){\varphi ^{n - 1}}(i,j) - {\lambda _z}{\varphi ^{n - 1}}(i{\rm{ }},j + 1)]\} \end{array}。$

(27)

其中：

$ \left\{ {\begin{array}{*{20}{l}} {{\lambda _x} = \frac{{cos{\alpha _1} + cos{\alpha _2}}}{{cos{\alpha _1}cos{\alpha _2}}} \cdot v \cdot \frac{{\Delta t}}{{\Delta x}}}\\ {{\mu _x} = \frac{1}{{cos{\alpha _1}cos{\alpha _2}}} \cdot {v^2} \cdot \frac{{\Delta {t^2}}}{{\Delta {x^2}}}}\\ {{\lambda _z} = \frac{{cos{\alpha _1} + cos{\alpha _2}}}{{cos{\alpha _1}cos{\alpha _2}}} \cdot v \cdot \frac{{\Delta t}}{{\Delta z}}}\\ {{\mu _z} = \frac{1}{{cos{\alpha _1}cos{\alpha _2}}} \cdot {v^2} \cdot \frac{{\Delta {t^2}}}{{\Delta {z^2}}}} \end{array}} \right. $

α₁和α₂是两个可选的角度参数(通常为零)。φ_l为第l层任意场分量的值，l = 1, 2, 3 …M, M为吸收边界总层数。σ_l是余弦型加权算子，${\sigma _l} = 1 - cos{\rm{ }}(\frac{\pi }{2}\cdot({\rm{ }}1 + \frac{{1 - l}}{M})) $。v是快纵波速度，求解方式见文献[16]。

分别将$\theta _{{\rm{re}}}^{n + 1}{\rm{ }}(i,j) $、$\theta _{{\rm{im}}}^{n + 1}{\rm{ }}(i,j) $、$\varepsilon _{{\rm{re}}}^{n + 1}{\rm{ }}(i,j) $、$\varepsilon _{{\rm{im}}}^{n + 1}{\rm{ }}(i,j) $当作${\varphi ^{n + 1}}(i,j) $，代入式(25)即可得到部分饱和孔隙介质的左吸收边界差分格式，其余边界可以类推。

为验证吸收边界效果，对均匀含油部分饱和孔隙介质进行正演模拟，获取波场快照。图 4为采用以上边界条件计算的均匀含油部分饱和孔隙介质中t=140 ms时的纵波波场快照，激发源设在网格中心，计算空间大小为1 000 m×1 000 m，空间步长为5 m×5 m，时间步长为0.2 ms，主频25 Hz，差分精度为16阶，吸收边界厚度为20层，介质参数见表 2，在该介质中α取值范围为1~2.25。在大多数情况下毛细管压力参数α取值一般接近$\frac{{{K_{{f_1}}}}}{{{K_{{f_2}}}}} $，即介质中存在毛细管压力影响，但毛细管压力并不占据主导地位，因此本文选择毛细管压力参数α为2.1。图 5为没有采用任何边界条件处理计算得到的均匀部分饱和孔隙介质中的纵波波场快照。比较图 4和图 5可以看出，图 4中快纵波在人工边界的伪反射几乎完全被吸收，而图 5中快纵波边界处伪反射强烈，说明上述边界条件能较好地压制边界的伪反射。

4 数值模拟

为验证部分饱和孔隙介质中改进纵波方程及其正演算法的正确性，本文设计三个模型进行正演模拟。模型1为均匀部分饱和孔隙介质，其物性参数如表 2所示。图 6为采用16阶空间精度差分的正演过程中不同时刻(70、90、110 ms)的固相和流相波场快照。计算空间大小为1 000 m×1 000 m，空间步长为5 m×5 m，时间步长为0.2 ms，激发源设在网格中心，主频25 Hz。

从图 6(a)、图 6(b)、图 6(c)中可以清楚地观察到，部分饱和孔隙介质中(孔隙中两种流体分别是油和水)，在固相波场同时存在相位相反的快纵波和慢纵波，且其速度与理论计算结果一致。在流相波场中(见图 6(d)、图 6(e)、图 6(f))同样可以观察到两类纵波，只是流相波场中快纵波的振幅比慢纵波的振幅要小很多，这与介质特性有关。综合对比图 6(a)—(f)可以发现快纵波的相位在固相部分和流相部分是相同的，慢纵波的相位在固相部分和流相部分是相反的，这与Biot理论一致，进一步说明了本文正演结果的正确性。观察图 6(a)—(f)还可发现慢纵波的频散现象比快纵波严重，这是由慢纵波自身的易频散特性决定的，这也与理论结果一致。

为考虑毛细管压力参数对波场的影响，在保持表 2中其他参数不变的前提下，计算了不同毛细管压力参数下的波场快照。图 7为t=100 ms时的波场快照，其中图 7(a)—(c)分别为α=1、α=1.6和α=2.25($\alpha = \frac{{{K_{{f_1}}}}}{{{K_{{f_2}}}}} $)时的固相波场快照，图 7(d)—(f)分别为α=1、α=1.6、α=2.25的流相波场快照。图 7中红色虚线和蓝色虚线位置固定，观察图 7(a)—(c)中快纵波与红色虚线关系，可以发现α=1时快纵波最外围已明显穿过红色虚线，α=2.25时快纵波最外围刚刚到达红色虚线，α=1.6时快纵波位置位于α=1和α=2.25中间，这说明毛细管压力越大(α越小)，快纵波速度越大。观察图 7(a)—(f)中快纵波与蓝色虚线关系，发现毛细管压力对慢纵波无明显影响。

为探究地震波在部分饱和孔隙介质中的衰减情况，设置了模型2，即如图 8所示的层状模型，其中介质1和介质4为页岩，厚度为500 m，介质2为致密砂岩，介质3为含油砂岩，介质2和介质3厚度为20 m。模型2物性参数如表 3所示，在介质3中α取值范围为1~3.214，本文选择介质3中毛细管压力参数α为3。

对模型2进行正演模拟，在地表设置激发点和接收点，以纵波源作为激发震源，采用自激自收方式接收反射波。激发主频为30 Hz，共有601道接收，道间距为10 m，采样间隔为0.2 ms，记录长度为450 ms，得到的叠后地震记录如图 9所示。模型2的目的层反射波即200~300 ms之间的地震波同相轴。图 10为从图 9中抽取的第50道和第300道地震记录，其中蓝色线为第50道地震记录，即未穿过含油储层的单道地震记录，红色线为第300道地震记录，即穿过含油储层的地震记录。观察图 10可发现含油地震道振幅大于未含油地震道，这是由于上下介质波阻抗差异导致的，含油地震道上下介质的波阻抗差大于未含油地震道。波阻抗的差异会影响反射系数，进而导致含油地震道的振幅大于未含油地震道的振幅。对第50道和第300道地震记录进行傅里叶变换得到图 11，观察图 11可知红色曲线高于蓝色曲线，这也是因为波阻抗的差异导致的。对图 11中振幅进行归一化可得到图 12。从图 12中可明显看出地震波在含油储层中传播时其主频小于未含油时的情况。为突出频谱变化，将两者相减(含油地震道振幅谱减未含油地震道振幅谱)得到图 13。

由图 13可见，在小于地震波主频的频率范围内含油地震道的振幅大于未含油地震道的振幅，而在大于地震波主频的频率范围内含油地震道的振幅小于未含油地震道的振幅。此现象表明地震波在穿过油气区后其主频发生变化，主频向低频方向移动，地震波低频能量相对增强，高频能量相对减弱。这是因为不同频率的地震波使得流体流动的剧烈程度不同，从而导致不同频率的波耗散也不同，也即油气(流体)对不同频率的波吸收程度不同，对频率高的波吸收程度大，对频率低的波吸收程度小，从而导致地震波的高频成分衰减大，低频成分衰减小，主频降低。这说明本文部分饱和孔隙介质中的改进纵波方程能够反映地震频带内流体对波的吸收衰减的影响，同时也说明本文正演算法正确。

为进一步分析毛细管压力对地震波衰减的影响，设置了模型3。模型3示意图可参考模型2，不同的是将介质2和介质3的厚度扩展至200 m，其目的是将顶底界面的反射波分开，避免顶底反射波相互干涉对振幅造成的影响，其余参数与模型2相同。对模型3进行正演模拟，其观测系统与模型2正演模拟时相同，得到的叠后地震记录如图 14所示。图 14(a)为储层毛细管压力参数α=1时的叠后地震记录，即毛细管压力最大时的叠后地震记录；图 14(b)为储层毛细管压力参数$\alpha = 3{\rm{ }}.{\rm{ }}2{\rm{ }}1{\rm{ }}4{\rm{ }}(\alpha = \frac{{{K_{{f_1}}}}}{{{K_{{f_2}}}}}) $时的叠后地震记录，即不存在毛细管压力影响时的叠后地震记录。分别抽取图 14(a)和图 14(b)的第300道地震记录得到图 15。

图 15中深绿色曲线为α=1时的地震记录，深红色曲线为α=3.214时的地震记录。在图 15中200~300 ms为顶界面的反射波，300~450 ms为底界面反射波。综合对比顶、底界面反射波，能够明显看出，在穿过油气区域之后(底界面反射波)，毛细管压力越大(α=1时)，地震波的振幅衰减越显著。

5 结语

本文将毛细管压力参数融入Johnson体积模量计算公式，并将改进后的体积模量计算公式代入部分饱和孔隙介质的纵波方程，从而构建出考虑毛细管压力的部分饱和孔隙介质中的改进纵波方程。为探究该方程对地震波传播特征的描述，本文推导了方程的高阶有限差分格式，利用有限差分法对方程进行了正演模拟，得到了高精度的波场快照，通过对同一时刻不同毛细管压力作用下的波场快照进行分析可知毛细管压力越大，快纵波速度越大。通过对薄层状模型正演所得叠后地震记录的抽道分析，发现穿过含油部分饱和孔隙介质的地震波主频显著降低，这与实际含油气储层的主频降低现象一致。通过对厚层状模型的叠后地震记录分析可知毛细管压力越大，地震波振幅衰减越明显。因此，改进后的部分饱和孔隙介质中的纵波方程能够反映流体流动引发的地震波衰减，同时也表明本文正演算法的正确性，为后续利用地震波的衰减预测储层流体提供了依据。

附录A

考虑毛细管压力的部分饱和孔隙介质体积模量K′的表达式为

$ K' = {K_{{\rm{BGH}}}} - {\rm{ }}\frac{{{K_{{\rm{BGH}}}} - K{'_{{\rm{BGW}}}}}}{{1 - \xi ' + \xi {\rm{ }}'\sqrt {1 - i\omega \tau '/\xi {'^2}} }}。$

(A-1)

其中:

$ \tau ' = {\left[ {\frac{{{K_{{\rm{BGH}}}} - K{'_{{\rm{BGW}}}}}}{{{K_{{\rm{BGH}}}}G}}} \right]^2}，$

(A-2)

$ \xi {\rm{ }}' = {\rm{ }}\frac{{\left( {{K_{{\rm{BGH}}}} - K{'_{{\rm{BGW}}}}} \right)}}{{2K{'_{{\rm{BGW}}}}}}\cdot\frac{{\tau '}}{{T'}}。$

(A-3)

上式中

$ G = {\rm{ }}\frac{{\kappa \cdot{K_{{\rm{BGH}}}}}}{{{\eta _1}\sqrt {{D_1}} + {\eta _2}\sqrt {{D_2}} }} \times {\rm{ }}{\left\{ {\frac{{{R_2} + {Q_2}\left[ {{K_1} + \left( {{\rm{ }}4{\rm{ }}/{\rm{ }}3} \right)N} \right]{\rm{ }} - {\rm{ }}\left( {{R_1} + {Q_1}} \right)\left[ {{K_2} + \left( {{\rm{ }}4{\rm{ }}/{\rm{ }}3} \right)N} \right]{\rm{ }}}}{{\phi {\rm{ }}{S_1}{K_1}\left[ {{K_2} + \left( {4{\rm{ }}/{\rm{ }}3{\rm{ }}} \right)N} \right] + \phi {\rm{ }}{S_2}{K_2}\left[ {{K_1} + \left( {4{\rm{ }}/{\rm{ }}3{\rm{ }}} \right)N} \right]}}} \right\}^2} \times \frac{S}{V}，$

(A-4)

$ N = \mu ，$

(A-5)

$ {K_1} = {K_{{\rm{BG}}}}\left( {{K_{{f_1}}}} \right)，$

(A-6)

$ {K_2} = {K_{{\rm{BG}}}}\left( {{K_{{f_2}}}} \right)，$

(A-7)

$ {R_1} = {\rm{ }}\frac{{{\phi ^2}{K_s}}}{{1 - \phi - {K_b}/{K_s} + \phi {K_s}/{K_{{f_1}}}}}，$

(A-8)

$ {R_2} = {\rm{ }}\frac{{{\phi ^2}{K_s}}}{{1 - \phi - {K_b}/{K_s} + \phi {K_s}/{K_{{f_2}}}}}，$

(A-9)

$ {Q_1} = \left[ {\left( {1 - \phi - {K_b}/{K_s}} \right)\phi {\rm{ }}{K_s}} \right]/{\rm{ }}\left[ {1 - \phi - {K_b}/{K_s} + \phi {\rm{ }}{K_s}/{K_{{f_1}}}} \right]，$

(A-10)

$ {Q_2} = \left[ {\left( {1 - \phi - {K_b}/{K_s}} \right)\phi {\rm{ }}{K_s}} \right]/{\rm{ }}\left[ {1 - \phi - {K_b}/{K_s} + \phi {\rm{ }}{K_s}/{K_{{f_2}}}} \right]，$

(A-11)

$ {D_1} = \frac{{\kappa {K_{{f_1}}}}}{{{\eta _1}\phi }}{\left( {1 + \frac{{{K_{{f_1}}}}}{{\phi {\rm{ }}\left( {{K_b} + \frac{4}{3}N} \right)}} \times \left\{ {1 + \frac{1}{{{K_s}}}\left[ {\frac{4}{3}N\left( {1 - \frac{{{K_b}}}{{{K_s}}}} \right) - {K_b} - \phi \left( {{K_b} + \frac{4}{3}N} \right)} \right]} \right\}} \right)^{ - 1}}，$

(A-12)

$ {D_2} = \frac{{\kappa {K_{{f_2}}}}}{{{\eta _2}\phi }}\left( {1 + \frac{{{K_{{f_2}}}}}{{\phi {\rm{ }}\left( {{K_b} + \frac{4}{3}N} \right)}} \times \left\{ {1 + \frac{1}{{{K_s}}}\left[ {\frac{4}{3}N\left( {1 - \frac{{{K_b}}}{{{K_s}}}} \right) - {K_b} - \phi \left( {{K_b} + \frac{4}{3}N} \right)} \right]} \right\}} \right)，$

(A-13)

$ T^{\prime}=\frac{K_{\mathrm{BGW}}^{\prime} \phi^2}{30 \kappa R_b^3}\left\{\begin{array}{l} {\left[3 \eta_2 g_2^{\prime 2}+5\left(\eta_1-\eta_2\right) \mathrm{g}_1^{\prime} g_2^{\prime}-3 \eta_1 g_1^{\prime 2}\right] R_a^5-15 \eta_2 g_2^{\prime}\left(g_2^{\prime}-g_1^{\prime}\right) R_a^3 R_b^2+} \\ 5 g_2^{\prime}\left[3 \eta_2 g_2^{\prime 2}-\left(2 \eta_2+\eta_1\right) \mathrm{g}_1^{\prime}\right] R_a^2 R_b^3-3 \eta_2 g_2^{\prime 2} R_b^5 \end{array}\right\}, $

(A-14)

$ g{'_1} = \frac{{\left( {1 - \frac{{{K_b}}}{{{K_s}}}} \right)\left( {\frac{1}{{K{'_{\rm{W}}}}} - \frac{1}{{{K_f}_{_1}}}} \right)}}{{1 - \frac{{{K_b}}}{{{K_s}}} - \phi \frac{{{K_b}}}{{{K_s}}} + \phi \frac{{{K_b}}}{{K{'_{\rm{W}}}}}}}，$

(A-15)

$ g{'_2} = \frac{{\left( {1 - \frac{{{K_b}}}{{{K_s}}}} \right)\left( {\frac{1}{{K{'_{\rm{W}}}}} - \frac{1}{{{K_f}_{_2}}}} \right)}}{{1 - \frac{{{K_b}}}{{{K_s}}} - \phi \frac{{{K_b}}}{{{K_s}}} + \phi \frac{{{K_b}}}{{K{'_{\rm{W}}}}}}}。$

(A-16)

在式(A - 4)中S是流体1和流体2之间的表面积, V是总体积, 因此有$S = 4{\rm{ }}\pi R_a^2 $, $V = \frac{4}{3}\pi R_b^3 $, $\frac{S}{V} = \frac{{3R_a^2}}{{R_b^3}} $。R_a和R_b分别为同心球体的内径和外径。