0 引言

叠前反演可揭示地下储层空间展布与含油气性等特征，在油气勘探、开发中具有重要作用^[1-3]。根据正演问题所用解析表达式的不同，叠前反演分为基于波动方程的反演和基于各种反射系数公式的反演^[4]。对于叠前三参数反演来说，受限于各种假设和条件^[5]，导致物理先验认识(AVA理论公式，如Aki公式^[6]、Gray公式^[7]、Fatti公式^[8]等)中的角度不易准确给定，加之大型稀疏矩阵的病态性，导致反演过程不稳定。为此，本文提出物理、数据先验认识融合(Harmonicon of Physical and Data Prior Knowledge, HPDPK)的叠前解耦分步反演方法。首先，对基于物理先验认识构建的正演框架进行非稀疏化，以保持反演的稳定性，同时也为后续解耦分步反演奠定基础；然后，以井资料作为数据先验认识，对正演框架进行修正，并利用修正后的正演框架对地震数据进行解耦，以获得稳定且准确的叠前地震属性数据；最后，利用井数据的低频信息作为先验约束，对叠前地震属性数据进行反演，获得地层弹性参数。

本文的反演方法如何获取高质量的叠前地震属性数据是关键。因常规方法利用物理先验认识直接提取的叠前地震属性数据不准确，不能与井合成记录匹配。作为改进，Levenberg-Marguarot方法在矩阵对角线引入单位矩阵和拉格朗日系数^[9-10]。后来，有人在贝叶斯概率化反演理论框架^[11]下，引入协方差^[12]以提高解耦精度。但在实际应用中，此方法对解耦性能的改善有限，根本原因在于无法获得符合工区实际情况的角度系数项，导致物理先验认识不适用。

为获取准确的角度项，本文将井控思想^[13-15]引入叠前三参数反演过程，提出物理、数据先验认识融合的地震数据解耦方法。通过地震数据与井数据的匹配，求取物理先验认识中的角度项，解决仅利用物理先验认识带来的解耦精度不足的问题。

1 方法原理 1.1 基于物理先验认识的非稀疏正演框架

以物理先验认识中的Gray公式为例介绍方法原理

$ R\left(\theta \right)={R}_{\lambda }\left(\frac{1}{2}-\frac{{\beta }^{2}}{{\alpha }^{2}}\right)\mathrm{s}\mathrm{e}{\mathrm{c}}^{2}\theta +\frac{{\beta }^{2}}{{\alpha }^{2}}\left(\mathrm{s}\mathrm{e}{\mathrm{c}}^{2}\theta -\\ \;\;\;\;\;\;\;\;4\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta \right){R}_{\mu }+\frac{1}{2}\left(1-\mathrm{t}\mathrm{a}{\mathrm{n}}^{2}\theta \right){R}_{\rho } $

(1)

式中：$ R\left(\theta \right) $为纵波反射系数，$ \theta $为界面处入射角与透射角的平均值；$ {R}_{\lambda }\mathrm{、}{R}_{\mu }\mathrm{、}{R}_{\rho } $分别为拉梅系数$ \lambda $、剪切模量$ \mu $和密度$ \rho $的反射系数；$ \beta /\alpha $为横波速度与纵波速度之比。对m个不同角度的反射系数，考虑子波$ w $后，可得到t时刻的叠前地震记录

$ d({\theta }_{i},t)=\left(\frac{1}{2}-\frac{{\beta }^{2}}{{\alpha }^{2}}\right)\mathrm{s}\mathrm{e}{\mathrm{c}}^{2}{\theta }_{\mathrm{i}}w\left(t\right){R}_{\mathrm{\lambda }}\left(t\right)+\frac{{\beta }^{2}}{{\alpha }^{2}}\left(\mathrm{s}\mathrm{e}{\mathrm{c}}^{2}{\theta }_{\mathrm{i}}-4\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}{\theta }_{\mathrm{i}}\right)w\left(t\right){R}_{\mathrm{\mu }}\left(t\right)+\frac{1}{2}\left(1-\mathrm{t}\mathrm{a}{\mathrm{n}}^{2}\theta \right)w\left(t\right){R}_{\mathrm{\rho }}\left(t\right) $

(2)

式中$ d({\theta }_{i},t) $为叠前地震记录，i=1, 2, …, m。将式(2)扩展至n个采样点，并写成矩阵形式

$\underbrace{\left[\begin{array}{c} d\left(\theta_1, t_1\right) \\ \vdots \\ d\left(\theta_1, t_n\right) \\ d\left(\theta_2, t_1\right) \\ \vdots \\ d\left(\theta_2, t_n\right) \\ \vdots \\ d\left(\theta_m, t_1\right) \\ \vdots \\ d\left(\theta_m, t_n\right) \end{array}\right]}_ \boldsymbol{d}= \underbrace{\left[ {\begin{array}{*{20}{l}} {{w_1}a\left( {{\theta _1}} \right)\;\;\;\;{w_1}b\left( {{\theta _1}} \right)\;\;\;\;{w_1}c\left( {{\theta _1}} \right)}\\ {\;\;\;\;\;\;\;\; \ddots \;\;\;\;\;\;\;\;\;\;\;\;\;\; \ddots \;\;\;\;\;\;\;\;\;\;\;\;\;\; \ddots }\\ {\;\;\;\;\;\;\;\;\;\;\;{w_n}a\left( {{\theta _1}} \right)\;\;\;\;{w_n}b\left( {{\theta _1}} \right)\;\;\;\;{w_n}c\left( {{\theta _1}} \right)}\\ {{w_1}a\left( {{\theta _2}} \right)\;\;\;\;{w_1}b\left( {{\theta _2}} \right)\;\;\;\;{w_1}c\left( {{\theta _2}} \right)}\\ {\;\;\;\;\;\;\;\; \ddots \;\;\;\;\;\;\;\;\;\;\;\;\;\; \ddots \;\;\;\;\;\;\;\;\;\;\;\;\;\; \ddots }\\ {\;\;\;\;\;\;\;\;\;\;\;{w_n}a\left( {{\theta _2}} \right)\;\;\;\;{w_n}b\left( {{\theta _2}} \right)\;\;\;\;{w_n}c\left( {{\theta _2}} \right)}\\ \ldots \\ {{w_1}a\left( {{\theta _m}} \right)\;\;\;\;{w_1}b\left( {{\theta _m}} \right)\;\;\;\;{w_1}c\left( {{\theta _m}} \right)}\\ {\;\;\;\;\;\;\;\; \ddots \;\;\;\;\;\;\;\;\;\;\;\;\;\; \ddots \;\;\;\;\;\;\;\;\;\;\;\;\;\; \ddots }\\ {\;\;\;\;\;\;\;\;\;\;\;{w_n}a\left( {{\theta _m}} \right)\;\;\;\;{w_n}b\left( {{\theta _m}} \right)\;\;\;\;{w_n}c\left( {{\theta _m}} \right)} \end{array}} \right]}_\boldsymbol{G} \underbrace{\left[\begin{array}{c}R_\lambda\left(t_1\right) \\ \vdots \\ R_\lambda\left(t_n\right) \\ R_\mu\left(t_1\right) \\ \vdots \\ R_\mu\left(t_n\right) \\ R_\rho\left(t_1\right) \\ \vdots \\ R_\rho\left(t_n\right)\end{array}\right]}_\boldsymbol{R}$

(3)

式中：$ {t}_{j}(j=\mathrm{1,2},\cdots ,n) $为时间序列；$ a\left(\theta \right)=\left(1/2-{\beta }^{2}/{a}^{2}\right)\mathrm{s}\mathrm{e}{\mathrm{c}}^{2}\theta $、$ b\left(\theta \right)=\left(\mathrm{s}\mathrm{e}{\mathrm{c}}^{2}\theta -4\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta \right){\beta }^{2}/{\alpha }^{2} $、$ c\left(\theta \right)=\left(1-\mathrm{t}\mathrm{a}{\mathrm{n}}^{2}\theta \right)/2 $为角度项；$ \boldsymbol{d} $为叠前地震数据；$ \boldsymbol{G} $为系数矩阵；$ \boldsymbol{R} $为反射系数。则式(3)表示为

$ {\boldsymbol{d}}_{mn\times 1}={\boldsymbol{G}}_{mn\times 3n}{\boldsymbol{R}}_{3n\times 1} $

(4)

从式(3)、式(4)不难看出，常规叠前反演希望从地震道集直接反演地层弹性参数。因不同参数相互耦合，且需求解大型稀疏矩阵$ \boldsymbol{G} $的逆矩阵，导致反演过程不稳定。同时，以下两个因素造成了矩阵$ \boldsymbol{G} $中由物理先验认识直接计算的$ a\left(\theta \right)\mathrm{、}b\left(\theta \right)\mathrm{、}c\left(\theta \right) $不准确性，从而加剧了反演的不稳定性：①式(3)、式(4)是以角度表示的，但地震数据是炮检距记录的，从炮检距域转化成角度域需要已知地下速度场。由于速度场难以准确获取，使从炮检距换算的角度存在较大偏差；②纵横波速度比不易确定。③叠前反演一般使用部分角度叠加数据，在反演中通常取角度范围的中间值。这些因素若用Zoeppritz方程进行反演，则存在求解困难且无法直接反演目标参数等问题。

在大数据领域，数据解耦是指将含有多个变量的数学方程变成用单个变量表示的方程组，从而简化计算。将此思路用于叠前三参数反演问题，如果消除参数之间的耦合作用，实现各个参数的独立存在，就可实现各弹性参数的稳定反演。对此，提出物理、数据先验认识融合的叠前解耦分步反演方法。首先根据褶积运算的交换律，将式(3)进行非稀疏化分解，得到如下正演框架

$ \underbrace{\left[\begin{array}{ccc} d\left(\theta_1, t_1\right) & \cdots & d\left(\theta_1, t_n\right) \\ d\left(\theta_2, t_1\right) & \cdots & d\left(\theta_2, t_n\right) \\ & \vdots & \\ d\left(\theta_m, t_1\right) & \cdots & d\left(\theta_m, t_n\right) \end{array}\right]}_\boldsymbol{d}=\underbrace{\left[\begin{array}{ccc} a\left(\theta_1\right) & b\left(\theta_1\right) & c\left(\theta_1\right) \\ a\left(\theta_2\right) & b\left(\theta_2\right) & c\left(\theta_2\right) \\ &\vdots & \\ a\left(\theta_m\right) & b\left(\theta_m\right) & c\left(\theta_m\right) \end{array}\right]}_\boldsymbol{A} \underbrace{\left[\begin{array}{ccc} R_\lambda\left(t_1\right) & \cdots & R_\lambda\left(t_n\right) \\ R_\mu\left(t_1\right) & \cdots & R_\mu\left(t_n\right) \\ R_\rho\left(t_1\right) & \cdots & R_\rho\left(t_n\right) \end{array}\right] \boldsymbol{w}}_{\boldsymbol{m} = \boldsymbol{R} \times \boldsymbol{w}} $

(5)

式(5)的矩阵形式如下

$ {\boldsymbol{d}}_{m\times n}={\boldsymbol{A}}_{m\times 3}{\boldsymbol{R}}_{3\times n}{\boldsymbol{w}}_{n\times n}={\boldsymbol{A}}_{m\times 3}{\boldsymbol{m}}_{3\times n} $

(6)

式中：A为与角度有关的系数矩阵，称为叠前角度项；$ \boldsymbol{w} $为子波矩阵；m为叠前地震属性。根据AVA近似式，纵波反射系数可表示为纵横波速度、密度、拉梅参数等反射系数的组合。通过不同的变换，可以将反射系数表示成不同参数的组合，这是在反射系数层面上对“解耦”的理解。同理，由式(6)可知，叠前地震记录可表示为不同参数的叠前地震属性的组合，这是从叠前地震数据的层面上对“解耦”的理解。通过变换，消除式(6)中A的作用，把叠前地震数据分解成各弹性参数对应的叠前地震属性数据，此过程称为叠前地震数据的解耦。最后，对各个叠前地震属性数据进行反演，获得对应的地层弹性参数。

对比式(3)和式(5)、式(4)和式(6)中的G和A可知，因G为稀疏矩阵，且矩阵维数较大，这给矩阵求逆带来了挑战。非稀疏化后，矩阵A的维数及病态性大大降低，为解耦分步反演奠定了基础。如何准确获取叠前地震属性数据是解耦分步反演的关键，具体解耦过程如下。

1.1.1 常规叠前数据解耦

对于式(6)，期望得到一组与叠前地震数据之间误差$ \boldsymbol{e} $平方和最小的预测数据对应的模型参数，其目标函数$ E $为

$ E={\boldsymbol{e}}^{\mathrm{T}}\boldsymbol{e}={(\boldsymbol{d}-\boldsymbol{A}\boldsymbol{m})}^{\mathrm{T}}(\boldsymbol{d}-\boldsymbol{A}\boldsymbol{m})\to \mathrm{m}\mathrm{i}\mathrm{n} $

(7)

对上式求导，可得

$ \boldsymbol{m}={\left({\boldsymbol{A}}^{T}\boldsymbol{A}\right)}^{-1}{\boldsymbol{A}}^{\mathrm{T}}\boldsymbol{d} $

(8)

从式(8)看出，利用叠前角度项$ \boldsymbol{A} $可将叠前地震数据分解为$ \lambda $、$ \mu $和$ \rho $表示的叠前地震属性数据。实际应用时，前面所提的种种原因使叠前角度项不准，导致式(8)提取的叠前地震属性精度较低，影响后续反演结果。

作为改进，引入协方差提高解耦精度^[16]。由贝叶斯理论^[17-18]，当$ \boldsymbol{m} $满足相关多元Gauss分布，噪声满足均值为零、方差为$ {\sigma }^{2} $的Gauss分布时，有

$ P\left(\boldsymbol{d}\left|\boldsymbol{m}\right.\right)=\frac{1}{\sigma \sqrt{2\mathrm{\pi }}}\mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{1}{2}\frac{{(\boldsymbol{A}\boldsymbol{m}-\boldsymbol{d})}^{\mathbf{T}}(\boldsymbol{A}\boldsymbol{m}-\boldsymbol{d})}{{\sigma }^{2}}\right] $

(9)

$ P\left(\boldsymbol{m}\right)=\frac{\mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{1}{2}{(\boldsymbol{m}-\overline{\boldsymbol{m}})}^{\mathrm{T}}{\boldsymbol{C}}_{\boldsymbol{m}}^{-1}(\boldsymbol{m}-\overline{\boldsymbol{m}})\right]}{{\left(2\mathrm{\pi }\right)}^{\frac{3}{2}}\sqrt{|\mathrm{d}\mathrm{e}\mathrm{t}{\boldsymbol{C}}_{\boldsymbol{m}}{|}^{3}}} $

(10)

式中：“det”表示行列式；$ {\boldsymbol{C}}_{\boldsymbol{m}} $为$ \boldsymbol{m} $的协方差，可由井数据统计得来^[19]；$ \overline{\boldsymbol{m}} $为$ \boldsymbol{m} $的均值。当$ \overline{\boldsymbol{m}}=0 $时，联合式(9)、式(10)可得$ \boldsymbol{m} $的后验概率为

$ \begin{array}{l}P\left(\boldsymbol{m}\right|\boldsymbol{d})\propto \mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{1}{2}\frac{{(\boldsymbol{A}\boldsymbol{m}-\boldsymbol{d})}^{\mathbf{T}}(\boldsymbol{A}\boldsymbol{m}-\boldsymbol{d})}{{\sigma }^{2}}\right]\times \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{1}{2}{\boldsymbol{m}}^{\mathrm{T}}{\boldsymbol{C}}_{\boldsymbol{m}}^{-1}\boldsymbol{m}\right)\end{array} $

(11)

取式(11)的对数，可得解耦目标函数

$ J\left(\boldsymbol{m}\right|\boldsymbol{d})=\frac{1}{2}{\boldsymbol{m}}^{\mathrm{T}}{\boldsymbol{C}}_{m}^{-1}\boldsymbol{m}+\frac{{(\boldsymbol{A}\boldsymbol{m}-\boldsymbol{d})}^{\mathbf{T}}(\boldsymbol{A}\boldsymbol{m}-\boldsymbol{d})}{2{\sigma }^{2}} $

(12)

令其梯度为零可得

$ \boldsymbol{m}={\left[{\boldsymbol{A}}^{\mathrm{T}}\boldsymbol{A}+{\sigma }^{2}{\boldsymbol{C}}_{m}^{-1}\right]}^{-1}{\boldsymbol{A}}^{\mathrm{T}}\boldsymbol{d} $

(13)

对比式(8)与式(13)，虽然协方差约束方法引入$ {\sigma }^{2}{\boldsymbol{C}}_{\boldsymbol{m}}^{-1} $改善解耦结果，但易受$ {\sigma }^{2} $的制约，且没从本质上考虑叠前角度项的影响。

1.1.2 物理、数据先验认识融合的叠前数据解耦

为获得精确的解耦结果，以测井数据为先验认识，利用HPDPK方法对叠前数据进行解耦。即利用井控计算的叠前角度项$ {\boldsymbol{A}}_{wc} $代替理论叠前角度项$ \boldsymbol{A} $进行数据解耦，过程为

$ \boldsymbol{d}=\boldsymbol{A}\boldsymbol{R}\boldsymbol{w}=\boldsymbol{A}\boldsymbol{m}\to \boldsymbol{d}={\boldsymbol{A}}_{wc}\boldsymbol{R}\boldsymbol{w}={\boldsymbol{A}}_{wc}\boldsymbol{m} $

(14)

具体步骤如下：

(1) 通过井震标定，将地震数据与井数据进行匹配。

(2) 利用井数据计算$ {R}_{\lambda }\mathrm{、}{R}_{\mu }\mathrm{、}{R}_{\rho } $，并与子波褶积得到合成叠前地震属性数据$ {\boldsymbol{m}}_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}} $。

(3) 利用d与$ {\boldsymbol{m}}_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}} $求解$ {\boldsymbol{A}}_{wc} $。求解过程为

$ \boldsymbol{e}=\boldsymbol{d}-\boldsymbol{A}{\boldsymbol{m}}_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}=0 $

(15)

$ \boldsymbol{L}={\boldsymbol{A}}^{\mathrm{T}}\boldsymbol{A}=\sum {\boldsymbol{A}}_{kl}^{2} $

(16)

即在$ \boldsymbol{e} $最小的条件下，求取使L最小的解^[20]。式中：k=1, 2, $ \cdots $, M; l=1, 2, $ \cdots $, N; M和N均为向量维度。

不失一般性，设地震数据$ \boldsymbol{d} $为M维向量，属性参数$ {\boldsymbol{m}}_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}} $为N维向量。当利用3个部分角度叠加地震数据进行三参数反演时，M=3、N=3。由式(15)、式(16)得目标函数

$ \varphi \left(\boldsymbol{A}\right)=\boldsymbol{L}+\sum\limits_{k=1}^{M}{\beta }_{k}{e}_{l}=\\\;\;\;\;\;\;\sum\limits_{k=1}^{M}\sum\limits_{l=1}^{N}{\boldsymbol{A}}_{kl}^{2}+\sum\limits_{k=1}^{M}{\beta }_{k}\left({\boldsymbol{d}}_{k}-\sum\limits_{l=1}^{N}{\boldsymbol{A}}_{kl}{{\boldsymbol{m}}_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}}_{l}\right) $

(17)

求式(17)的梯度并令其为

$ \begin{array}{l}\frac{\partial \varphi \left(\boldsymbol{A}\right)}{\partial \boldsymbol{A}}=\sum\limits_{k=1}^{M}\sum\limits_{l=1}^{N}2\frac{\partial {\boldsymbol{A}}_{kl}^{}}{\partial \boldsymbol{A}}{\boldsymbol{A}}_{kl}^{}-\sum\limits_{k=1}^{M}{\beta }_{k}\sum\limits_{l=1}^{N}\frac{\partial {\boldsymbol{A}}_{kl}^{}}{\partial \boldsymbol{A}}{\boldsymbol{m}}_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}l}\\ \;\;\;\;\;\;\;\;\;\; =2\sum\limits_{k=1}^{M}\sum\limits_{l=1}^{N}{\boldsymbol{A}}_{kl}^{}-\sum\limits_{k=1}^{M}{\beta }_{k}\sum\limits_{l=1}^{N}{\boldsymbol{m}}_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}l}=0\end{array} $

(18)

由式(18)可得

$ 2\boldsymbol{A}=\boldsymbol{\beta }{\boldsymbol{m}}_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}^{\mathrm{T}} $

(19)

式中$ \boldsymbol{\beta }=[{\beta }_{1}{\beta }_{2}\cdots {\beta }_{M}] $为加权因子。由式(14)和式(19)得

$ \boldsymbol{\beta }=2\boldsymbol{d}({\boldsymbol{m}}_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}^{\mathrm{T}}{\boldsymbol{m}}_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}{)}^{-1} $

(20)

将式(20)代入式(19)，消元后得到$ \boldsymbol{A} $的最小能量解$ {\boldsymbol{A}}_{wc} $

$ {\boldsymbol{A}}_{wc}=\boldsymbol{d}({\boldsymbol{m}}_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}^{\mathrm{T}}{\boldsymbol{m}}_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}{)}^{-1}{\boldsymbol{m}}_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}^{\mathrm{T}} $

(21)

根据式(21)可以得到符合实际工区的叠前角度项$ {\boldsymbol{A}}_{wc} $。

(4) 将$ {\boldsymbol{A}}_{wc} $代入式(8)，对叠前地震数据解耦可得

$ \boldsymbol{m}={\left[{\boldsymbol{A}}_{wc}^{\mathrm{T}}{\boldsymbol{A}}_{wc}\right]}^{-1}{\boldsymbol{A}}_{wc}^{\mathrm{T}}\boldsymbol{d} $

(22)

通过式(22)可得到精确的叠前地震属性数据。

1.2 叠前地震属性的井约束反演

鉴于地震反演的多解性，引入测井低频信息约束和稀疏约束^[21-22]，可得反演目标函数

$ F\left(\boldsymbol{r}\right)={‖\boldsymbol{m}-\boldsymbol{w}\boldsymbol{r}‖}_{2}^{2}+\gamma {‖\boldsymbol{r}-{\boldsymbol{r}}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}}‖}_{2}^{2}+\eta {‖\boldsymbol{r}‖}_{1} $

(23)

式中：r为$ {R}_{\lambda } $、$ {R}_{\mu } $、$ {R}_{\rho } $中的某个参数，视输入的叠前地震属性数据m而定；$ \gamma $和$ \eta $为正则化因子。式(23)第一项为合成地震记录与叠前地震属性数据间的残差；第二项为先验信息约束，其中$ {\boldsymbol{r}}_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}} $为测井低频信息；第三项是稀疏约束。

输入不同的叠前地震属性数据，利用梯度投影稀疏重构法^[23]对式(23)求解，得到$ {R}_{\lambda } $、$ {R}_{\mu } $、$ {R}_{\rho } $，再对其道积分即可得到各自对应的地层弹性参数。

2 数值处理与结果分析 2.1 模型数据验证

利用某工区W48井的测井曲线，根据Zoeppritz方程得到不同角度的纵波反射系数，与主频30 Hz的Ricker子波褶积得到1°~30°的叠前地震数据(图 1)。分别将1°~10°、11°~20°、21°~30°的数据部分叠加，得到小角度、中角度、大角度地震数据。

利用角度部分叠加数据作为d与W48井中$ \lambda \mathrm{、}\mu \mathrm{、}\rho $的合成属性数据$ {\boldsymbol{m}}_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}} $，通过式(21)得到$ {\boldsymbol{A}}_{wc} $，由Gray公式得到$ \boldsymbol{A} $。根据式(14)，由$ {\boldsymbol{A}}_{wc} $得到基于HPDPK方法的合成地震记录，由$ \boldsymbol{A} $得到基于Gray公式的合成地震记录。图 2a~图 2c分别为小角度、中角度和大角度地震数据，对比可知，Gray公式与Zoeppritz方程的正演结果在中角度和大角度数据中存在很大差异；而HPDPK方法的正演结果A_wc与Zoeppritz方程正演结果吻合度高，且很好地保留了AVA的特征信息，说明了HPDPK方法的可靠性。

2.2 实际数据应用与结果分析

根据式(21)和式(22)，将HPDPK方法应用于某工区实际地震资料1°~10°、11°~20°、21°~30°三个部分角度叠加数据和W48井数据，解耦获得$ \lambda \mathrm{、}\mu \mathrm{、}\rho $的叠前地震属性数据，并与利用物理先验认识直接解耦(式(8))和协方差约束方法(式(13))的结果进行对比。图 3展示了各方法解耦的λ叠前地震属性剖面，其中图 3a为利用物理先验认识的直接解耦结果，可见同相轴连续性差，属性结果与井合成属性数据$ {\boldsymbol{m}}_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}} $不匹配；图 3b为协方差约束方法的解耦结果，尽管同相轴连续性及$ {\boldsymbol{m}}_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}} $与解耦结果的对应关系有了一定程度的改善，但部分位置的匹配性仍然较差(如图中红色圆圈所示)；图 3c为HPDPK方法的解耦结果，同相轴的连续性最好，且与$ {\boldsymbol{m}}_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}} $良好匹配。

为更细致地观察各方法解耦结果与$ {\boldsymbol{m}}_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}} $的吻合情况，提取W48井的井旁属性数据进行对比(图 4)。图 4a~图 4c分别为$ \lambda \mathrm{、}\mu \mathrm{、}\rho $的叠前地震属性。由图 4可见：直接解耦结果与$ {\boldsymbol{m}}_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}} $匹配最差；协方差约束对解耦结果的改善有限，在图 4a的1.35、1.44 s处，图 4b的1.24、1.49 s处和图c的1.25、1.49 s处波谷不匹配；相比之下，HPDPK方法的解耦结果与$ {\boldsymbol{m}}_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}} $匹配最好。

为进一步验证HPDPK方法的稳定性，图 5给出了W344盲井的验证结果。由图可见，直接解耦的结果、协方差约束的解耦结果与W344盲井的$ {\boldsymbol{m}}_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}} $数据吻合度较差，如图 5a的1.21、1.41 s处，图 5b的1.31、1.37 s处和图 5c的1.39 s处箭头所示；而利用HPDPK方法的解耦结果仍与W344盲井的$ {\boldsymbol{m}}_{\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}} $保持良好的对应关系。

根据式(23)，利用HPDPK解耦的叠前地震属性数据进行反演，得到相应的地层弹性参数，如图 6~图 8所示。图 6是剪切模量的反演结果，可见反演结果的低值区（蓝色区域）对应测井曲线低值，高值区（红色区域）对应测井曲线高值，反演结果与井曲线吻合度较高。图 7、图 8分别是拉梅参数和密度反演结果。与图 6类似，反演结果均与井曲线吻合。

为更直观地观察HPDPK方法的反演精度，提取协方差约束方法和HPDPK方法在W48井处的反演曲线，并与原始井曲线对比，如图 9所示。由图可见，HPDPK方法的反演结果与测井曲线最相似，其中$ \lambda \mathrm{、}\mu \mathrm{、}\rho $与对应测井曲线的相似系数分别为89.1%、90.5%和87.4%；而协方差约束反演结果对应的相似系数分别只有78.6%、74.3%和76.5%。进一步将W48井、W344盲井和W47井处反演结果与测井曲线的相似系数绘制成表 1，可见较协方差约束方法，本方法反演的$ \lambda \mathrm{、}\mu \mathrm{、}\rho $精度分别提高了14.1%、13.6%和11.9%。

3 结论

针对叠前三参数反演存在的不稳定和精度不高问题，提出了物理、数据先验认识融合(HPDPK)的叠前解耦分步反演方法。先对AVA正演框架进行非稀疏化分解，增加了参数反演的稳定性；然后利用HPDPK方法从叠前部分叠加数据中解耦获得叠前地震属性数据，避免了角度不准而带来的巨大偏差，有利于提高地层弹性参数的反演精度；最后通过叠前地震属性数据反演得到对应的地层弹性参数。模型和实际数据试验结果表明，相比于业界常用的叠前反演方法，本文方法获得的叠前地震属性及反演结果更稳定，而且具有更高的精度。

本文方法将物理、数据先验认识融合，可以利用叠前数据进行任意敏感参数的反演，不限于模量、纵横波速度、泊松比等弹性参数，也可以自然扩展到AVAZ五维地震数据反演。