发布日期：2025-01-04 15:19 点击次数：127
0 引言
近几年，随着叠前地震反演技术和多波多分量技术的发展，储层流体识别方法的有效性和准确度不断提高。利用叠前振幅随偏移距或角度变化信息来进行流体识别是一种重要的方法[1-4]。印兴耀等[5-6]通过含有Gassmann流体项的弹性阻抗方程直接反演来提高流体识别的准确性；Zong等[7]建立了纵波模量和横波模量表示的孔隙弹性理论和反射系数近似方程，并发展了基于纵波模量和横波模量AVO反演的流体识别方法；张世鑫[8]发展了基于固液解耦的地震流体识别方法，克服了常规流体识别方法受孔隙度影响产生的流体识别假象，然而，这些方法仅使用了反射纵波地震数据，在岩性和流体的共同影响下，利用纵波得到的反演结果往往存在精度低、反演不稳定、边界刻画不清晰等问题。
与纵波速度相比，横波速度受流体影响较小，因此在反演中同时利用纵波和转换横波地震数据，可以降低反演的多解性，提高储层预测精度[9-11]。Larsen等[12]提出了纵波和转换波AVO联合加权叠加反演方法；陈天胜等[13]基于一种方向加速度优化算法，利用纵波和转换横波反演得到了稳定的速度比值；Du等[14]以弹性孔隙介质理论为基础，将流体因子引入Russell反射系数公式，实现了多波联合AVO反演，但常规流体因子由于受孔隙度和孔隙流体的耦合作用影响，流体识别存在的多解性问题仍未得到有效解决。
综合前人的研究结果，首先讨论Gassmann流体因子和流体体积模量随孔隙度的变化关系，然后针对复杂储层中Gassmann流体因子受孔隙流体和孔隙度耦合影响而存在的流体识别假象问题，推导基于流体体积模量、剪切模量和孔隙度的PP波和PS波Zoeppritz线性近似方程，建立一种能够稳定获取等效流体体积模量、剪切模量和孔隙度的纵横波同步联合反演方法，以期解决单纵波反演精度低、反演不稳定、边界刻画不清晰的问题。
1 理论方法
1.1 孔隙流体体积模量的敏感度分析
Domenico[15]假设流体的体积模量是气体和液体体积模量的体积加权函数
$
{K_{\rm{f}}} = {S_{\rm{w}}}{K_{\rm{w}}} + \left( {1 - {S_{\rm{w}}}} \right){K_{\rm{g}}}
$
(1)
式中：Kf为流体体积模量，GPa；Sw为含水饱和度，%；Kw为液体体积模量，GPa；Kg为气体体积模量，GPa。
Han等[16]提出了Gassmann流体项经验公式
$
f = G\left( \mathit{\Phi } \right){K_{\rm{f}}}
$
(2)
式中：f 为流体因子项，GPa; $ G(\mathit{\Phi }) = \frac{{{{\left( {1 - {K_{\rm{n}}}} \right)}^2}}}{\mathit{\Phi }}, {K_{\rm{n}}} = \frac{{{K_{{\rm{dry}}}}}}{{{K_{\rm{m}}}}}, G(\mathit{\Phi })$为岩石骨架矿物与孔隙度的综合增益函数；Φ 为孔隙度，%；Kdry为干岩石体积模量，GPa；Km为固体矿物基质体积模量，GPa；Kn为干岩石体积模量与基质体积模量的比值。
Nur等[17]给出了临界孔隙度模型
$
\left\{ {\begin{array}{*{20}{l}}
{{K_{{\rm{dry}}}} = {K_{\rm{m}}}\left( {1 - \frac{\mathit{\Phi }}{{{\mathit{\Phi }_{\rm{c}}}}}} \right)}\\
{{\mu _{{\rm{dry}}}} = {\mu _{\rm{m}}}\left( {1 - \frac{\mathit{\Phi }}{{{\mathit{\Phi }_{\rm{c}}}}}} \right)}
\end{array}} \right.
$
(3)
式中：μdry为干岩石剪切模量，GPa；μm为固体基质剪切模量（GPa）；Φc为临界孔隙度，%。对增益函数G (Φ)进行简化，得到
$
G\left( \mathit{\Phi } \right) = \frac{\mathit{\Phi }}{{\mathit{\Phi }_{\rm{c}}^2}}
$
(4)
假设固体骨架矿物相同，孔隙中包含水和气，则Gassmann流体因子和流体体积模量随含水饱和度和孔隙度的变化趋势（图 1）便可由式（1）与式（2）计算得到。由图 1可看出，Gassmann流体因子受孔隙流体和孔隙度的综合影响，即与孔隙度和含水饱和度呈非线性关系[图 1（a）]，而流体体积模量在一定孔隙度下随含水饱和度的变化趋势则为完全线性[图 1（b）]。因此，如果将流体体积模量与孔隙度分别单独作为变量，即可消除常规流体因子受孔隙度和孔隙流体的耦合影响，解决流体识别多解性的问题。
下载原图
图 1 流体因子随孔隙度与含水饱和度的变化趋势
Fig. 1 Trend of fluid factor with porosity and water saturation
1.2 Kf - μ - Φ 纵横波孔隙模量三参数AVO近似方程和反演
1.2.1 Kf - μ - Φ 纵横波孔隙模量三参数AVO近似方程推导
Russell等[18]建立了弹性孔隙介质理论和AVO技术之间的关系，推导了关于流体因子、剪切模量和密度的三参数AVO近似式
$
\begin{array}{*{20}{c}}
{{R_{{\rm{PP}}}} \approx \left[ {\left( {\frac{1}{4} - \frac{{\gamma _{{\rm{dry}}}^2}}{{4\gamma _{{\rm{sat}}}^2}}} \right){{\sec }^2}\theta } \right]\frac{{\Delta f}}{f} + \left[ {\frac{{\gamma _{{\rm{dry}}}^2}}{{4\gamma _{{\rm{sat}}}^2}}{{\sec }^2}\theta - } \right.}\\
{\left. {\frac{2}{{\gamma _{{\rm{sat}}}^2}}{{\sin }^2}\theta } \right]\frac{{\Delta \mu }}{\mu } + \left[ {\frac{1}{2} - \frac{{{{\sec }^2}\theta }}{4}} \right]\frac{{\Delta \rho }}{\rho }}
\end{array}
$
(5)
$
{R_{{\rm{ps}}}} \approx \frac{{{\gamma _{{\rm{san}}}}\tan \varphi }}{2}\left[ {\left( {\frac{2}{{\gamma _{{\rm{sat}}}^2}}{{\sin }^2}\theta - \frac{2}{{{\gamma _{{\rm{sat}}}}}}\cos \theta \cos \varphi } \right)\frac{{\Delta \mu }}{\mu } - \frac{{\Delta \rho }}{\rho }} \right]
$
(6)
式中：θ 为PP波的入射角与透射角的平均值，（°）；φ为PS波反射角与透射角的平均值，（°）；μ 为剪切模量，GPa；ρ 为密度，g/cm3；γsat为饱和岩石条件下纵横波速度比均值；γdry为干岩石情况下的纵横波速度比。
Gassmann流体项及剪切模量表达式为
$
f = {\left( {\rho v_{\rm{p}}^2} \right)_{{\rm{sat}}}} - \gamma _{{\rm{dry}}}^2{\left( {\rho v_{\rm{S}}^2} \right)_{{\rm{sat}}}}
$
(7)
$
\mu = \rho v_{\rm{S}}^2
$
(8)
式中：vP，vS分别为饱和岩石纵波、横波速度，m/s。由式（7）和式（8）的流体因子、剪切模量变化率、纵横波速度变化率与密度变化率的关系，可以得到
$
\frac{{\Delta f}}{f} = \frac{{{{\left( {v_{\rm{p}}^2} \right)}_{{\rm{sat}}}}\Delta \rho {{\left( {{v_{\rm{P}}}} \right)}_{{\rm{sat}}}}\Delta {v_{\rm{p}}} - \gamma _{{\rm{dry}}}^2\left[ {{{\left( {v_{\rm{S}}^2} \right)}_{{\rm{sat}}}}\Delta \rho + 2\rho {{\left( {{v_{\rm{S}}}} \right)}_{{\rm{sat}}}}\Delta {v_{\rm{S}}}} \right]}}{{{{\left( {\rho v_{\rm{P}}^2} \right)}_{{\rm{sat}}}} - \gamma _{{\rm{dry}}}^2{{\left( {\rho v_{\rm{S}}^2} \right)}_{{\rm{sat}}}}}} = \left( {\frac{{2\gamma _{{\rm{sat}}}^2}}{{\gamma _{{\rm{sat}}}^2 - \gamma _{{\rm{dry}}}^2}} + \frac{1}{4}} \right)\frac{{\Delta {v_{\rm{P}}}}}{{{v_{\rm{P}}}}} - \frac{{2\gamma _{{\rm{dry}}}^2}}{{\gamma _{{\rm{sat}}}^2 - \gamma _{{\rm{dry}}}^2}}\frac{{\Delta {v_{\rm{S}}}}}{{{v_{\rm{S}}}}}
$
(9)
$
\frac{{\Delta \mu }}{\mu } = 2\frac{{\Delta {v_{\rm{S}}}}}{{{v_{\rm{S}}}}} + \frac{{\Delta \rho }}{\rho } = 2\frac{{\Delta {v_{\rm{S}}}}}{{{v_{\rm{S}}}}} + \frac{1}{4}\frac{{\Delta {v_{\rm{P}}}}}{{{v_{\rm{P}}}}}
$
(10)
根据Gardner等[19]提出的经验关系式，得到纵波速度与密度变化率之间的关系
$
\frac{{\Delta \rho }}{\rho } = \frac{1}{4}\frac{{\Delta {v_{\rm{P}}}}}{{{v_{\rm{P}}}}}
$
(11)
结合式（9）—（11），可以推导得到
$
\frac{{\Delta \rho }}{\rho } = \frac{{\gamma _{{\rm{sat}}}^2 - \gamma _{{\rm{dry}}}^2}}{{9\gamma _{{\rm{sat}}}^2}}\frac{{\Delta f}}{f} + \frac{{\gamma _{{\rm{dry}}}^2}}{{9\gamma _{{\rm{sat}}}^2}}\frac{{\Delta \mu }}{\mu }
$
(12)
将式（12）分别代入式（5）和式（6），并可将其简化得到
$
{R_{{\rm{PP}}}} \approx \left[ {\frac{{2\left( {\gamma _{{\rm{sat}}}^2 - \gamma _{{\rm{dry}}}^2} \right)}}{{9\gamma _{{\rm{sat}}}^2}}{{\sec }^2}\theta + \frac{{\gamma _{{\rm{sat}}}^2 - \gamma _{{\rm{dry}}}^2}}{{18\gamma _{{\rm{sat}}}^2}}} \right]\frac{{\Delta f}}{f} + \left[ {\frac{{2\gamma _{{\rm{dry}}}^2}}{{9\gamma _{{\rm{sat}}}^2}}{{\sec }^2}\theta - \frac{2}{{\gamma _{{\rm{sat}}}^2}}{{\sin }^2}\theta + \frac{{\gamma _{{\rm{dry}}}^2}}{{18\gamma _{{\rm{sat}}}^2}}} \right]\frac{{\Delta \mu }}{\mu }
$
(13)
$
{R_{{\rm{PS}}}} \approx \frac{{{\gamma _{{\rm{sat}}}}\tan \varphi }}{2}\left[ {\left( {\frac{2}{{\gamma _{{\rm{sat}}}^2}}{{\sin }^2}\theta - \frac{2}{{{\gamma _{{\rm{sat}}}}}}\cos \theta \cos \varphi - \frac{{\gamma _{{\rm{dry}}}^2}}{{9\gamma _{{\rm{sat}}}^2}}} \right)\frac{{\Delta \mu }}{\mu } - \frac{{\gamma _{{\rm{sat}}}^2 - \gamma _{{\rm{dry}}}^2}}{{9\gamma _{{\rm{sat}}}^2}}\frac{{\Delta f}}{f}} \right]
$
(14)
根据Han等[16]总结得到的流体体积模量Kf与Gassmann流体项f 之间的关系，联立式（2）和式（4）可以得到
$
\frac{{\Delta f}}{f} = \frac{{\Delta \mathit{\Phi }}}{\mathit{\Phi }} + \frac{{\Delta {K_{\rm{f}}}}}{{{K_{\rm{f}}}}}
$
(15)
将式（15）代入式（13）和式（14），可以得到
$
\begin{array}{l}
{R_{{\rm{PP}}}} \approx \left[ {\frac{{2\left( {\gamma _{{\rm{sat}}}^2 - \gamma _{{\rm{dry}}}^2} \right)}}{{9\gamma _{{\rm{sat}}}^2}}{{\sec }^2}\theta + \frac{{\gamma _{{\rm{sat}}}^2 - \gamma _{{\rm{dry}}}^2}}{{18\gamma _{{\rm{sat}}}^2}}} \right]\frac{{\Delta {K_{\rm{f}}}}}{{{K_{\rm{f}}}}} + \\
\;\;\;\;\;\;\;\;\left[ {\frac{{2\gamma _{{\rm{dry}}}^2}}{{9\gamma _{{\rm{sat}}}^2}}{{\sec }^2}\theta - \frac{2}{{\gamma _{{\rm{sat}}}^2}}{{\sin }^2}\theta + \frac{{\gamma _{{\rm{dry}}}^2}}{{18\gamma _{{\rm{sat}}}^2}}} \right]\frac{{\Delta \mu }}{\mu } + \\
\;\;\;\;\;\;\;\;\left[ {\frac{{2\left( {\gamma _{{\rm{sat}}}^2 - \gamma _{{\rm{dry}}}^2} \right)}}{{9\gamma _{{\rm{sat}}}^2}}{{\sec }^2}\theta + \frac{{\gamma _{{\rm{sat}}}^2 - \gamma _{{\rm{dry}}}^2}}{{18\gamma _{{\rm{sat}}}^2}}} \right]\frac{{\Delta \mathit{\Phi }}}{\mathit{\Phi }}
\end{array}
$
(16)
$
\begin{array}{*{20}{c}}
{{R_{{\rm{PS}}}} \approx \frac{{{\gamma _{{\rm{sat}}}}\tan \varphi }}{2}\left[ {\left( {\frac{2}{{\gamma _{{\rm{sat}}}^2}}{{\sin }^2}\theta - \frac{2}{{{\gamma _{{\rm{sat}}}}}}\cos \theta \cos \varphi - \frac{{\gamma _{{\rm{dry}}}^2}}{{9\gamma _{{\rm{sat}}}^2}}} \right) \cdot } \right.}\\
{\left. {\frac{{\Delta \mu }}{\mu } - \frac{{\gamma _{{\rm{sat}}}^2 - \gamma _{{\rm{dry}}}^2}}{{9\gamma _{{\rm{sat}}}^2}}\frac{{\Delta {K_{\rm{f}}}}}{{{K_{\rm{f}}}}} - \frac{{\gamma _{{\rm{sat}}}^2 - \gamma _{{\rm{dry}}}^2}}{{9\gamma _{{\rm{bat}}}^2}}\frac{{\Delta \Phi }}{\Phi }} \right]}
\end{array}
$
(17)
式（16）和式（17）即为推导的PP波和PS波三参数Zoeppritz线性近似公式。
1.2.2 近似方程精度分析
为验证推导方程的精度[20]，分别设计了3种砂岩模型（表 1），分别用精确Zoeppritz方程[21]，Russell方程[18]与式（16），式（17）计算得到地层界面处PP波、PS波的反射系数，其中γdry2 = 2.333。当入射角不超过30°时，基于Kf -μ -Φ 的纵波近似式与Russell近似式精度相近，而转换波近似方程的精度则略高于Russell近似式精度[图 2（a）—（b），图 2（e）—（f）]。基于Kf - μ - Φ 的纵波近似式与Russell近似精度相近，而转换波近似方程精度在大于30°时变差[图 2（c）—（d）]。综合以上分析，新推导的孔隙模量三参数AVO近似方程在入射角0°~30°内具有较高的精度，因此在最佳入射角0°~30°时利用本文提出的Kf - μ - Φ 方程进行纵横波同步联合反演具有可行性。
下载原图
图 2 3种模型的PP波、PS波反射系数对比
Fig. 2 Comparison of PP wave and PS wave reflection coefficient of three kinds of models
下载CSV
表 1 砂岩模型参数
Table 1 Sandstone model parameters
1.3 基于贝叶斯理论的叠前同步联合反演
对PP波和PS波进行同步联合反演，考虑M层n 个入射角的情况，将式（16）和式（17）写成矩阵形式
$
{\left[ {\begin{array}{*{20}{c}}
{R_{{\rm{PP}}}^1}\\
\vdots \\
{R_{{\rm{PP}}}^n}\\
{R_{{\rm{PS}}}^1}\\
\vdots \\
{R_{{\rm{PS}}}^n}
\end{array}} \right]_{2n \times 1}} = {\left[ {\begin{array}{*{20}{c}}
{{A_1}}&{{B_1}}&{{C_1}}\\
\vdots & \vdots & \vdots \\
{{A_n}}&{{B_n}}&{{C_n}}\\
{{D_1}}&{{E_1}}&{{F_1}}\\
\vdots & \vdots & \vdots \\
{{D_n}}&{{E_n}}&{{F_n}}
\end{array}} \right]_{2n \times 3}}{\left[ {\begin{array}{*{20}{c}}
{\frac{{\Delta {K_{\rm{f}}}}}{{{K_{\rm{f}}}}}}\\
{\frac{{\Delta \mu }}{\mu }}\\
{\frac{{\Delta \Phi }}{\Phi }}
\end{array}} \right]_{3 \times 1}}
$
(18)
其中
$
A = \left[ {\frac{{2\left( {\gamma _{{\rm{sat}}}^2 - \gamma _{{\rm{dry}}}^2} \right)}}{{9\gamma _{{\rm{sat}}}^2}}{{\sec }^2}\theta + \frac{{\gamma _{{\rm{sat}}}^2 - \gamma _{{\rm{dry}}}^2}}{{18\gamma _{{\rm{sat}}}^2}}} \right],
$
$
B = \left[ {\frac{{2\gamma _{{\rm{dry}}}^2}}{{9\gamma _{{\rm{sat}}}^2}}{{\sec }^2}\theta - \frac{2}{{\gamma _{{\rm{sat}}}^2}}{{\sin }^2}\theta + \frac{{\gamma _{{\rm{dry}}}^2}}{{18\gamma _{{\rm{sat}}}^2}}} \right],
$
$
C = \left[ {\frac{{2\left( {\gamma _{{\rm{sat}}}^2 - \gamma _{{\rm{dry}}}^2} \right)}}{{9\gamma _{{\rm{sat}}}^2}}{{\sec }^2}\theta + \frac{{\gamma _{{\rm{sat}}}^2 - \gamma _{{\rm{dry}}}^2}}{{18\gamma _{{\rm{sat}}}^2}}} \right],
$
$
D = - \frac{{{\gamma _{{\rm{sat}}}}\tan \varphi }}{2}\frac{{\gamma _{{\rm{sat}}}^2 - \gamma _{{\rm{dry}}}^2}}{{9\gamma _{{\rm{sat}}}^2}},
$
$
E = \frac{{{\gamma _{{\rm{sat}}}}\tan \varphi }}{2}\left( {\frac{2}{{\gamma _{{\rm{sat}}}^2}}{{\sin }^2}\theta - \frac{2}{{{\gamma _{{\rm{sat}}}}}}\cos \theta \cos \varphi - \frac{{\gamma _{{\rm{dry}}}^2}}{{9\gamma _{{\rm{sat}}}^2}}} \right),
$
$
F = - \frac{{{\gamma _{{\rm{sat}}}}\tan \varphi }}{2}\frac{{\gamma _{{\rm{sat}}}^2 - \gamma _{{\rm{dry}}}^2}}{{9\gamma _{{\rm{sat}}}^2}}
$
根据地震褶积模型理论，增加子波矩阵至式（18），则叠前地震道集可以通过矩阵乘法表示为
$
\mathit{\boldsymbol{d}} = \left[ {\begin{array}{*{20}{c}}
{\mathit{\boldsymbol{d}}_{{\rm{PP}}}^1}\\
\vdots \\
{\mathit{\boldsymbol{d}}_{{\rm{PP}}}^n}\\
{\mathit{\boldsymbol{d}}_{{\rm{PS}}}^1}\\
\vdots \\
{\mathit{\boldsymbol{d}}_{{\rm{PS}}}^n}
\end{array}} \right] = {\left[ {\begin{array}{*{20}{c}}
{{\mathit{\boldsymbol{W}}_{\rm{P}}}{A_1}}&{{\mathit{\boldsymbol{W}}_{\rm{P}}}{B_1}}&{{\mathit{\boldsymbol{W}}_{\rm{P}}}{C_1}}\\
\vdots & \vdots & \vdots \\
{{\mathit{\boldsymbol{W}}_{\rm{P}}}{A_n}}&{{\mathit{\boldsymbol{W}}_{\rm{P}}}{B_n}}&{{\mathit{\boldsymbol{W}}_{\rm{P}}}{C_n}}\\
{{\mathit{\boldsymbol{W}}_{\rm{S}}}{D_1}}&{{\mathit{\boldsymbol{W}}_{\rm{S}}}{E_1}}&{{\mathit{\boldsymbol{W}}_{\rm{S}}}{F_1}}\\
\vdots & \vdots & \vdots \\
{{\mathit{\boldsymbol{W}}_{\rm{S}}}{D_n}}&{{\mathit{\boldsymbol{W}}_{\rm{S}}}{E_n}}&{{\mathit{\boldsymbol{W}}_{\rm{S}}}{F_n}}
\end{array}} \right]_{2n \times 3}}{\left[ {\begin{array}{*{20}{c}}
{\frac{{\Delta {K_{\rm{f}}}}}{{{K_{\rm{f}}}}}}\\
{\frac{{\Delta \mu }}{\mu }}\\
{\frac{{\Delta \mathit{\Phi }}}{\mathit{\Phi }}}
\end{array}} \right]_{3 \times 1}}
$
(19)
式中：WP和WS分别为PP波和PS波的子波矩阵。
式（19）可以写成
$
\mathit{\boldsymbol{d}} = \mathit{\boldsymbol{Gr}}
$
(20)
式中：d 为PP波和PS波地震记录矩阵；G 为含有子波的系数矩阵；r 为待反演的由等效体积模量、剪切模量和孔隙度组成的列向量系数矩阵，即
$
\mathit{\boldsymbol{G}} = \left[ {\begin{array}{*{20}{c}}
{{\mathit{\boldsymbol{W}}_{\rm{P}}}{A_1}}&{{\mathit{\boldsymbol{W}}_{\rm{P}}}{B_1}}&{{\mathit{\boldsymbol{W}}_{\rm{P}}}{C_1}}\\
\vdots & \vdots & \vdots \\
{{\mathit{\boldsymbol{W}}_{\rm{P}}}{A_n}}&{{\mathit{\boldsymbol{W}}_{\rm{P}}}{B_n}}&{{\mathit{\boldsymbol{W}}_{\rm{P}}}{C_n}}\\
{{\mathit{\boldsymbol{W}}_{\rm{S}}}{D_1}}&{{\mathit{\boldsymbol{W}}_{\rm{S}}}{E_1}}&{{\mathit{\boldsymbol{W}}_{\rm{S}}}{F_1}}\\
\vdots & \vdots & \vdots \\
{{\mathit{\boldsymbol{W}}_{\rm{S}}}{D_n}}&{{\mathit{\boldsymbol{W}}_{\rm{S}}}{E_n}}&{{\mathit{\boldsymbol{W}}_{\rm{S}}}{F_n}}
\end{array}} \right],
$
$
\mathit{\boldsymbol{r}} = {\left[ {\frac{{\Delta {K_{\rm{f}}}}}{{{K_{\rm{f}}}}}\frac{{\Delta \mu }}{\mu }\frac{{\Delta \mathit{\Phi }}}{\mathit{\Phi }}} \right]^T}
$
为提高反演的稳定性，本文采用基于贝叶斯理论[22]的PP波和PS波同步联合反演。贝叶斯反演的目的就是在给定测量数据d（带噪声n）的情况下估计模型参数r，把似然函数赋给噪音分布函数
$
p\left( {\mathit{\boldsymbol{d}},\mathit{\boldsymbol{r}}} \right) = p\left( n \right)
$
(21)
假设噪音服从高斯分布，噪音的平均值为n，方差为σ2，在这种情况下，概率函数为
$
p\left( {{n_i}} \right) = \frac{1}{{\sqrt {2{\rm{ \mathsf{ π} }}{\sigma ^2}} }}\exp \left[ {\frac{{ - 1}}{{2{\sigma ^2}}}{{\left( {{n_i} - \bar n} \right)}^2}} \right]
$
(22)
假定噪音的平均值n等于0，这一项就可以省略。如果噪音是不相关的，并且每种噪音样本的方差为常数，则似然函数可以写为
$
p\left( {\mathit{\boldsymbol{d}}\left| \mathit{\boldsymbol{r}} \right.} \right) = \frac{1}{{\sqrt {2{\rm{ \mathsf{ π} }}} {\sigma _n}}}\exp \left[ {\frac{{ - {{\left( {\mathit{\boldsymbol{d}} - \mathit{\boldsymbol{Gr}}} \right)}^T}\left( {\mathit{\boldsymbol{d}} - \mathit{\boldsymbol{Gr}}} \right)}}{{2\sigma _n^2}}} \right]
$
(23)
式中：σn2为噪音的方差。为了使反演结果稳定而进行的约束可以看作先验信息，通过讨论先验概率分布，可以更直观地理解对解施加的约束条件。用贝叶斯公式可将似然函数与柯西先验分布结合起来，则r
的后验概率密度变为
$
\begin{array}{*{20}{c}}
{p\left( {\mathit{\boldsymbol{r}},{\sigma _n}\left| {\mathit{\boldsymbol{d}},I} \right.} \right) \propto {K_1}{K_2}\prod\limits_{i = 1}^M {\left[ {\frac{1}{{1 + r_i^2/\sigma _{\rm{r}}^2}}} \right] \cdot } }\\
{\exp \left[ { - \frac{{{{\left( {\mathit{\boldsymbol{Gr}} - \mathit{\boldsymbol{d}}} \right)}^T}\left( {\mathit{\boldsymbol{Gr}} - \mathit{\boldsymbol{d}}} \right)}}{{2\sigma _n^2}}} \right]}
\end{array}
$
(24)
式中：σr2 表示模型的方差；${K_1}{K_2} = \frac{1}{{{{\left( {{\rm{ \mathsf{ π} }}{\sigma _{\rm{r}}}} \right)}^M}}}\frac{1}{{\sqrt {2\pi } {\sigma _n}}}$。
将（24）式代入到边缘化公式，取对数并进一步求导得到目标函数为
$
\nabla \mathit{\boldsymbol{L}} = \mathit{\boldsymbol{Qr}} - \frac{{\left( {{\mathit{\boldsymbol{G}}^T}\mathit{\boldsymbol{Gr}} - {\mathit{\boldsymbol{G}}^T}\mathit{\boldsymbol{d}}} \right)}}{{\tilde \sigma _n^2}}
$
(25)
其中
$
\tilde \sigma _n^2 = \frac{{{{\left( {\mathit{\boldsymbol{G\hat r}} - \mathit{\boldsymbol{d}}} \right)}^T}\left( {\mathit{\boldsymbol{G\hat r}} - \mathit{\boldsymbol{d}}} \right)}}{{\left( {N - 1} \right)}},
$
$
\mathit{\boldsymbol{Q}} = \left[ {\begin{array}{*{20}{c}}
{\frac{1}{{{{\left( {1 + \frac{{r_1^2}}{{\sigma _{\rm{r}}^2}}} \right)}^2}}}}&0&0&0\\
0&{\frac{1}{{{{\left( {1 + \frac{{r_2^2}}{{\sigma _{\rm{r}}^2}}} \right)}^2}}}}&0&0\\
0&0& \ddots &0\\
0&0&0&{\frac{1}{{{{\left( {1 + \frac{{r_M^2}}{{\sigma _{\rm{r}}^2}}} \right)}^2}}}}
\end{array}} \right]
$
令式（25）等于0，就得到问题的贝叶斯解
$
\left( {{\mathit{\boldsymbol{G}}^T}\mathit{\boldsymbol{G}} + \mu \mathit{\boldsymbol{Q}}} \right)\mathit{\boldsymbol{r}} = {\mathit{\boldsymbol{G}}^T}\mathit{\boldsymbol{d}}
$
(26)
式中：μ为数据保真项与模型约束项之间权重的调节因子。
反演式（26）存在一定的弱非线性，为得到稳定的求解结果，采用共轭梯度法进行求解。
2 应用实例
2.1 Marmousi2模型试算
为验证本方法的正确性，本文利用Marmousi2油藏模型[23]进行测试。图 3为利用纵、横波速度，密度以及孔隙度等参数，计算得到模型等效流体体积模量、剪切模量和孔隙度（Kf - μ - Φ）剖面。在此模型参数基础上，通过式（16）和式（17）正演得到反射系数，并与30 Hz雷克子波褶积生成纵波和转换波角度道集，同时加入信噪比为2的高斯噪音。图 4是CDP 300处PP波与PS波的角道集，其入射角为0°~30°。
下载原图
图 3 等效流体体积模量（a）、剪切模量（b）和孔隙度（c）模型
Fig. 3 Models of equivalent fluid bulk modulus(a), shear modulus(b)and porosity(c)
下载原图
图 4 CDP 300处PP波角道集（a）与PS波角道集（b）
Fig. 4 PP wave angle gathers(a)and PS wave angle gathers(b)at CDP 300
利用中心角度为5°，15°，25°部分的叠加数据，分别进行单纵波和纵横波同步联合反演。由图 5中CDP 300处单纵波、纵横波同步联合反演与实际模型单道对比结果可以得知，本文方法的反演结果误差相对较小，稳定性高于单纵波反演结果；图 6为单纵波和纵横波同步联合反演结果对比，其中，等效流体体积模量剖面能清楚地反映含油气砂岩特征，剪切模量剖面反映了砂岩与页岩的岩性差别，孔隙度剖面则准确反映了储层的孔隙特征。综合对比结果，采用本文纵横波同步联合反演方法得到的结果较单纵波反演结果具有更高的精确度，且能准确刻画地质体储层特征和含气特征。
下载原图
图 5 CDP 300处单纵波、纵横波同步联合反演结果与实际模型单道对比
Fig. 5 Single channel comparison of single P-wave and PP wave and SS wave with actual model at CDP 300
下载原图
图 6 单纵波和纵横波联合反演对比
Fig. 6 Comparison of single P-wave and PP-wave and SS wave joint inversion
2.2 实际数据应用
利用中国LJ地区A油田M测线的纵波和转换横波数据进行应用效果测试。该区位于凹陷南缓斜坡，受燕山运动和喜马拉雅运动的影响，形成了向北倾没的大型鼻状构造，受构造控制，该区油气藏类型主要为断块油气藏，兼有岩性、泥岩裂缝等特殊油气藏[24]。本区纵波地震资料主频为35 Hz，转换波地震资料主频为15 Hz。由于实际纵波和转换横波记录的反射振幅相差较大，为得到更可靠的联合反演结果，需要对纵波振幅和转换横波振幅做匹配。以纵波振幅为标准，每一角度计算一个校正系数[式（27）]，为了不改变整个角道集的AVO特征，将各个角度的校正系数求平均作为最终的校正系数，图 7为在时间域内匹配一致的纵波和转换波角道集。
下载原图
图 7 CDP 1000点PP波角道集（a）与PS波角道集（b）
Fig. 7 PP wave angle gathers(a)and PS wave angle gathers(b)at CDP 1000
$
{S'_{{\rm{dataPS}}}} = \frac{{\int_0^{{T_{{\rm{max }}}}} {\left| {{S_{{\rm{data PP }}}}\left( {t',\theta } \right)} \right|} {\rm{d}}t'}}{{\int_0^{{T_{{\rm{mat }}}}} {\left| {{S_{{\rm{dataps }}}}\left( {t',\theta } \right)} \right|} {\rm{d}}t'}}{S_{{\rm{dataPS}}}}\left( {t,\theta } \right)
$
(27)
利用反演得到的M测线目标层段局部放大的等效流体体积模量[图 8（a）]、剪切模量[图 8（b）]和孔隙度剖面[图 8（c）]。图中黑色线柱表示生产井段内的油层，紫色线柱表示生产井段内的水层，经过对比发现，油层对应等效流体体积模量剖面低值[图 8（a）中红色所示]，孔隙度剖面高值[图 8（c）中红色所示]，而水层等效流体体积模量值高于油层[图 8（a）中蓝色所示]，相对应的孔隙度剖面低值[图 8（c）中蓝色所示]。油水层所在砂岩与泥页岩区在剪切模量剖面上有较明显的区分[图 8（b）]。综上所述，反演结果中储层特征清晰，油、水层得到了较好的区分，且油层范围与已知生产井的信息较吻合。
下载原图
图 8 局部放大的纵横波联合反演的M测线目标层段反演结果
Fig. 8 Local amplified inversion results of M-line target interval in the joint inversion of PP wave and SS wave
3 结论
（1）新推导的孔隙模量三参数AVO近似方程在入射角0°~30°内具有较高的精度，因此在最佳入射角0°~30°时利用本文提出的方程进行纵横波同步联合反演具有可行性。
（2）综合理论分析与数据测试结果，本文提出的纵横波同步联合孔隙模量三参数AVO反演方法改善了单纵波反演结果不稳定和边界模糊性的问题，并在一定程度上压制了噪声，提高了反演精度和流体识别的可靠性。
致谢: 中国石化胜利油田分公司物探研究院孟宪军为本文提供了实际地震资料，谨此表示感谢。
