20世纪80年代, 勘探工作者发现地震数据的反射振幅随炮检距的变化可以用于油气检测, 从而发展了利用纵波AVO属性进行构造解释和岩性识别的AVO分析技术^[1]。在AVO分析技术的基础上发展起来的AVO反演方法, 可从地震叠前资料中获得反映地下储层信息的纵、横波速度和密度等信息, 为储层预测与烃类检测提供量化依据。随着AVO反演方法逐渐趋于成熟, 人们开始利用地震、测井和岩心等数据的互补特性, 从井间地震数据中反演地层岩石的弹性参数和物性参数^[2-3]。与单纯的AVO反演相比, 这种多源信息联合反演的方法不仅可以减少反演多解性的问题, 还能够同时获得多种储层参数，因此受到了学术界和工业界的广泛关注。

多源信息联合反演研究多于一个的目标函数在给定区域上的最优化, 求取目标函数最小或最大值所对应的一组模型参数, 并把这组模型参数看成是满足这些条件的“最佳解”。显然, 目标函数不同, 最佳解也就不一样。为了求解方便, 绝大多数方法都假设噪声信号服从高斯分布^[4-7], 利用统计性方法和确定性方法求解目标函数。前者将观测数据和模型参数都视为随机变量, 用统计的方法来确定解所服从的概率分布, 得到的是满足观测数据的解具有多大概率; 而后者认为观测数据和模型参数都是确定量, 因而用数理方程或代数方法求解反演问题, 根据观测数据求得模型参数。BOSCH等^[8]在贝叶斯框架下将地震似然函数、岩石物理似然函数和先验信息结合, 构建服从高斯分布的后验概率密度函数, 利用统计性方法(蒙特卡罗采样)获得储层的孔隙度和波阻抗。统计性方法存在不完全依赖于初始猜测、反演过程理论上不会陷入局部极值的优点。但其缺点也十分明显, 就是计算工作量大、效率低, 难以完全满足大规模地球物理反演的要求。BOSCH等^[9-10]在贝叶斯框架下利用确定性方法(最小二乘法)对信噪比相对较高的地震叠后资料进行了孔隙度和波阻抗同步反演。李志勇等^[11]在BOSCH等研究成果的基础上尝试利用AVO正演模拟和Gasmann方程建立储层物性参数与地震叠前观测数据之间的联系, 利用最小二乘法同步反演获得纵横波速度、密度、孔隙度、含水饱和度和泥质含量。研究发现, 利用最小二乘准则可以方便地计算最速下降方向和修正步长, 并且具有不错的收敛速度。但对于信噪比相对较低的叠前地震数据, 使用最小二乘准则会使误差数据得到相同的权重, 严重影响反演结果的稳定性。另外, 由于速度模型与埋藏深度通常存在正比关系, 而储层物性参数服从0~1的截断分布, 因此，建立储层物性参数与弹性参数之间的岩石物理模型时, 误差同样不完全服从高斯分布^[12]。LI等^[13]采用自适应加权系数来调整地震似然函数、岩石物理似然函数和先验信息之间的权重, 提高了多参数同步反演方法的稳定性, 但无法从根本上消除脉冲噪声信号的影响。

对于脉冲噪声信号, TAYLOR等^[14]使用L₁范数建立目标函数, 得到了比传统的L₂范数更鲁棒的解。然而, L₁范数在残差接近于0的时候存在奇异点, 容易导致反演算法的不稳定。在TAYLOR等研究成果的基础上, 很多学者研究了基于L₁范数的变形算法。如BUBE等^[15]采用迭代加权算法来最小化L₁/L₂混合范数。GUITTON等^[16]、吕晓春等^[17]采用Huber范数建立目标函数, 在残差较小时采用L₂范数, 在残差较大时采用L₁范数, 避免了L₁范数的奇异性带来的优化困难。Huber范数在一定程度上减小了脉冲噪声对反演算法的影响, 但由于阈值点为无穷大, 也没有从根本上消除脉冲噪声的影响。随着脉冲噪声比例的增加, Huber范数将会失效。考虑到Biweight范数能够对大于设定阈值的噪声信号予以忽略, JI^[18]将Biweight范数引入到地震反演中, 提高了算法对脉冲噪声的抑制能力。刘洋等^[19]提出了针对非高斯噪声的地震叠前非高斯反演概念和思想, 构造了能同时压制高斯噪声和脉冲噪声的混合范数并将其作为目标函数。岳碧波等^[20]假设地震噪声信号服从非高斯α稳定分布, 提出了基于非高斯α稳定分布的最小L_p范数地震反演方法。另外, 考虑到广义极值概率密度分布函数在一定情况下具有逼近任意概率密度分布函数的能力, ZHANG等^[21]提出了基于广义极值分布的AVO反演方法。

误差协方差矩阵是L₂范数的一个重要统计性质, 其正确估计是进行数据完美匹配的前提。由于协方差矩阵对脉冲噪声的存在非常敏感, 因此脉冲噪声可能导致对协方差矩阵估计的错误, 解决这一问题最直接的办法就是采用鲁棒性估计方法修正协方差矩阵。本文假设地震噪声信号和岩石物理模型误差服从高斯分布且带有加性脉冲噪声, 引入Hampel三截尾函数构建重加权算子来对协方差矩阵进行修正。Hampel三截尾函数结合了Huber范数和Biweight范数的优点, 不仅能够减小脉冲噪声对协方差矩阵估计的影响, 而且还能将超过阈值的噪声信号完全排除, 提高算法对高斯噪声和脉冲噪声的抑制能力。

地震资料能够提供纵、横波速度和密度等信息, 利用这些信息可以进一步估算孔隙度、饱和度和泥质含量等储层物性参数。地震勘探中通常使用Gassmann方程来定量描述低频条件下储层物性参数发生变化时介质有效弹性模量的相应变化。采用基于Gassmann方程的流体替换技术可以得到不同孔隙流体条件下砂岩储层的体积模量和剪切模量^[22]:

(1)

(2)

式中:K_sat和μ_sat分别是流体饱和岩石的有效体积模量和有效剪切模量; K_dry和μ_dry分别是干岩石的有效体积模量和有效剪切模量, 可以采用改进的Hashin-Strikman弹性模量上边界计算得到^[23]; K₀是岩石基质的有效体积模量; K_fl是孔隙流体的有效体积模量; φ是孔隙度。

对于任意岩石, 给定各成分体积含量之后, 其等效模量将处于上、下界限之间, 但其精确的数值将依赖于几何细节。由于Gassmann-Biot理论假设岩石具有一个均匀的矿物模量, 因此混合矿物常常用一种“平均矿物”来描述。本文采用最简单、但不一定是最好的界限。以砂、泥岩为例, 如果岩石的矿物成分已知, 那么利用Voigt-Reuss-Hill平均可以计算岩石基质的等效体积模量K₀^[24]:

(3)

式中:C是泥质含量, K_s和K_c分别是砂岩和泥岩的体积模量。

同理, 当岩石孔隙所含流体为多相流体时, 如果流体的成分已知, 那么利用Reuss模型^[25]可以计算孔隙流体的等效体积模量K_fl:

(4)

式中:S_w是含水饱和度; K_hc是油或气的体积模量; K_brine是地层水的体积模量。

建立岩石物理模型的另一个参数是饱和岩石的密度ρ_sat, 该参数由岩石骨架密度ρ₀和充填流体的密度ρ_fl按体积比例加权平均得到:

(5)

式中:ρ_sat, ρ₀, ρ_fl分别为流体饱和岩石、岩石骨架和充填流体的密度。

在计算得到流体饱和岩石的体积模量K_sat, 剪切模量μ_sat和密度ρ_sat之后, 流体饱和岩石的纵波、横波速度公式可以表示为:

(6)

(7)

文献[11]通过假设地震噪声信号和岩石物理模型误差服从高斯分布, 给出了基于最小二乘法的储层弹性和物性参数同步反演目标函数:

(8)

式中:d是地震角道集向量; x是储层弹性参数向量; G是地震叠前正演模型算子; C_d是地震观测数据的协方差矩阵; y是储层物性参数向量; f(y)是建立储层物性参数y与弹性参数x联系的岩石物理方程; C_x|y是以储层物性参数y为条件、关于弹性参数x的协方差矩阵; C_y是储层物性参数y的协方差矩阵; y是储层物性参数y的先验信息向量。

通过研究(x, y)为何值时使S(x, y)达到最小可以得到最佳解。然而, 公式(8) 中的协方差矩阵C_d和C_x|y对脉冲噪声十分敏感, 容易导致反演算法的不稳定。本文假设地震噪声信号和岩石物理模型误差服从高斯分布且带有加性脉冲噪声, 利用Hampel三截尾函数构建重加权算子对协方差矩阵进行修正。新的反演目标函数表示为:

(9)

式中:ω_d和ω_x|y是基于Hampel三截尾函数的重加权算子^[26], 可表示为:

(10)

(11)

(12)

式中:ω(e)=φ(e)/e, φ(e)=∂ρ(e)/∂e, ρ(e)是Hampel三截尾函数; ε, Δ₁和Δ₂是Hampel三截尾函数的阈值, 用来控制函数对脉冲噪声的抑制能力, 阈值越小, 其抗脉冲噪声的能力越强。图 1为Hampel函数和加权算子示意图。可以看出:当|e_i|＜ε时, Hampel函数为二次函数; 当ε≤|e_i|＜Δ₂时, Hampel函数为一次函数; 当误差|e_i|＞Δ₂时, Hampel函数为常数。由此可见, Hampel函数结合了Huber函数和Biweight函数的优点:当误差较小时, 使用L₂范数求解; 当存在少数脉冲噪声时, 使用L₁范数求解; 当存在大量脉冲噪声时, 直接对噪声信号进行截断, 从根本上排除了脉冲噪声对反演的影响。

图 1(Fig. 1)

图 1 Hampel函数(a)和由Hampel函数构建的加权函数(b)

显然, 阈值的合理设置是压制高斯和脉冲噪声信号的关键。在噪声信号服从高斯分布的假设下, 如果噪声信号大于某一给定阈值T, 则其概率可表示为:

(13)

式中, σ表示高斯噪声信号的标准差。相应地, p{|e|>ε}, p{|e|>Δ₁}和p{|e|>Δ₂}分别表示|e|大于ε, Δ₁和Δ₂的概率。ZOU等^[27]提出根据数据的95.0%, 97.5%和99.0%置信区间来设置ε, Δ₁和Δ₂, 计算获得相应的阈值:

(14)

(15)

(16)

式中:σ(k)是迭代残差的瞬时标准差。一种常用的方法是利用递归公式对瞬时标准差估计进行滤波:

(17)

由于σ(k)容易受到脉冲噪声的影响, 这里采用一种鲁棒性的估计方法:

(18)

式中:mad{·}表示平均绝对中值偏差, 用来抑制脉冲噪声对方差估计的影响。

将公式(9) 的梯度用泰勒公式展开, 可以近似地得到:

(19)

式中:∇²S(x_n, y_n)是S(x, y)的Hessen矩阵, ∇S(x_n+1, y_n+1)和∇S(x_n, y_n)是S(x, y)的梯度向量。令∇S(x_n+1, y_n+1)=0, 可以得到:

(20)

式中:

(21)

(22)

式中:F^T是储层弹性参数x对物性参数y的梯度。为了方便求解, 本文将公式(20) 中的∇²S(x_n, y_n)化简为上三角矩阵。在公式(20) 的左右两边分别左乘储层物性参数和弹性参数联合的协方差矩阵C_{x, y}= 的逆矩阵C_{x, y}^-1= , 可以得到:

(23)

式中:D=F(y-Y)+f(y)-x。通过迭代求解方程组(23), 可以同步反演得到储层弹性参数x和物性参数y。

首先由测井曲线解释得到地层岩石的孔隙度、含水饱和度和泥质含量等储层物性参数(图 2a至图 2c), 然后根据岩石物理模型计算纵、横波速度和密度等弹性参数(图 2d至图 2f), 最后利用叠前正演模型合成地震角道集(图 3), 测试其含脉冲噪声时本文方法的反演效果。地震子波采用30Hz的零相位雷克子波, 道集的最小炮检距为150m, 最大炮检距为4000m, 炮检距间隔为350m。在反演过程中, 初始模型和先验模型可以通过测井曲线滤波得到, 物性参数的初始模型和先验模型相同。

图 2(Fig. 2)

图 2 利用测井曲线解释得到的储层参数(一) a孔隙度; b含水饱和度; c泥质含量; d纵波速度; e横波速度; f密度

图 3(Fig. 3)

图 3 利用测井数据合成的地震角道集(一)

当地震角道集含有30%的随机脉冲噪声时(图 4), 分别采用最小二乘法和本文方法进行反演, 结果如图 5和图 6所示。由图 5和图 6可以看出, 两者均能识别反射层位, 但最小二乘法的反演结果不稳定, 锯齿现象明显。这是因为使用最小二乘法会使残差得到相同的权重, 脉冲噪声的存在使地震数据得到偏大的权重(C_d偏小), 岩石物理约束和先验信息的权重相对减小, 弹性参数和物性参数的反演结果主要来自于地震数据。本文利用Hampel三截尾函数设置加权矩阵, 自适应调整C_d的大小, 减小了地震脉冲噪声在目标函数中的权重, 因而提高了算法对脉冲噪声的抑制能力。

图 4(Fig. 4)

图 4 加入30%脉冲噪声信号的地震角道集(一)

图 5(Fig. 5)

图 5 利用最小二乘法反演得到的储层参数(一) a孔隙度; b含水饱和度; c泥质含量; d纵波速度; e横波速度; f密度

图 6(Fig. 6)

图 6 利用本文方法反演得到的储层参数(一) a孔隙度; b含水饱和度; c泥质含量; d纵波速度; e横波速度; f密度

首先由测井曲线解释得到岩石的孔隙度、含水饱和度和泥质含量等储层物性参数(图 7a至图 7c); 然后根据岩石物理模型计算纵、横波速度和密度等弹性参数(图 7d至图 7f, 红色曲线), 并在弹性参数中加入脉冲噪声(图 7d至图 7f, 黑色曲线); 最后利用含脉冲噪声的弹性参数合成地震角道集(图 8), 并在地震角道集中加入30%的随机脉冲噪声(图 9), 分别测试岩石物理模型含脉冲噪声与岩石物理模型和地震角道集都含脉冲噪声两种情况下本文方法的反演效果。地震子波采用30Hz的零相位雷克子波, 道集的最小炮检距为150m, 最大炮检距为4000m, 炮检距间隔为350m。在反演过程中, 初始模型和先验模型可以通过测井曲线滤波得到, 物性参数的初始模型和先验模型相同。

图 7(Fig. 7)

图 7 利用测井曲线解释得到的储层参数(二) a孔隙度; b含水饱和度; c泥质含量; d纵波速度; e横波速度; f密度

图 8(Fig. 8)

图 8 利用测井数据合成的地震角道集(二)

图 9(Fig. 9)

图 9 加入30%随机脉冲噪声后的地震角道集(二)

当岩石物理模型含脉冲噪声、地震角道集不含脉冲噪声时(图 8), 分别采用最小二乘法和本文方法进行反演, 结果如图 10和图 11所示。由图 10和图 11可以看出, 两者均能识别反射层位, 但最小二乘法的弹性参数反演结果与真实值相差较大。这是因为脉冲噪声的存在使岩石物理模型得到偏大的权重(C_x|y偏小), 地震角道集和先验信息的权重相对减小, 弹性参数的反演结果主要来自于岩石物理约束。本文利用Hampel三截尾函数设置加权矩阵, 自适应调整C_x|y的大小, 减小了岩石物理脉冲噪声在目标函数中的权重, 提高了算法对脉冲噪声的抑制能力。

图 10(Fig. 10)

图 10 利用最小二乘法反演得到的储层参数(二) a孔隙度; b含水饱和度; c泥质含量; d纵波速度; e横波速度; f密度

图 11(Fig. 11)

图 11 利用本文方法反演得到的储层参数(二) a孔隙度; b含水饱和度; c泥质含量; d纵波速度; e横波速度; f密度

当岩石物理模型和地震角道集都含脉冲噪声时(图 9), 分别采用最小二乘法和本文方法进行反演, 结果如图 12和图 13所示。由图 12和图 13可以看出, 两者均能识别反射层位, 但使用最小二乘法得到的弹性参数反演结果不稳定且偏离真实值。这是因为最小二乘法会使残差数据得到相同的权重, 增加了脉冲噪声在目标函数中的权重。同时, 地震数据和岩石物理约束条件以及先验信息之间的相对关系变得不稳定, 反演结果既有可能来自于地震数据, 也有可能来自于岩石物理约束。本文方法能够自适应调整地震数据、岩石物理约束和先验信息之间的相对关系, 提高算法的稳定性, 因此弹性参数和物性参数的反演结果都与真实值吻合较好。

图 12(Fig. 12)

图 12 利用最小二乘法反演得到的储层参数(二) a孔隙度; b含水饱和度; c泥质含量; d纵波速度; e横波速度; f密度

图 13(Fig. 13)

图 13 利用本文方法反演得到的储层参数(二) a孔隙度; b含水饱和度; c泥质含量; d纵波速度; e横波速度; f密度

利用苏里格气田实际数据对本文方法的应用效果进行了研究。研究区目的层是上古生界二叠系底部H₈—S₁砂层组, 属沼泽背景下的辫状河沉积, 天然气成藏条件复杂。同时, 构造对天然气聚集不起主要控制作用, 储层经过强烈的成岩作用改造, 有效孔隙以次生孔隙为主, 具有低孔、低渗和非均质性严重的特点。图 14为苏里格气田A测线地震剖面, 该测线共有601道, 道间距为25m, 时间采样率为2ms。该剖面经过6口井, 各井所在CDP号分别为1665, 1691, 1792, 1856, 2007和2044, 图中黑色曲线为H₈层位。

图 14(Fig. 14)

图 14 苏里格气田叠后地震剖面

图 15是A测线地震角道集的方差分布图, 可以看出, 沿着该测线方向的方差变化较大。在这种情况下, 地震数据、岩石物理约束和先验信息之间的相对权重关系十分复杂, 基于噪声信号高斯分布的假设往往会导致反演结果不稳定。图 16a是位于CDP1856处的实际地震角道集, 可以看出, 1970ms处角道集的振幅远远大于上覆地层的振幅。图 16b是利用最小二乘法反演结果生成的地震角道集, 图 16c是实际角道集与合成角道集的残差分布图, 可以看出, 当上覆地层的地震残差很小时, 1970ms处的地震残差仍然很大。显然, 图 16c中的角道集残差已经不再满足噪声信号服从高斯分布的假设, 假设地震残差服从高斯分布且带有加性脉冲噪声更为合理。

图 15(Fig. 15)

图 15 地震角道集的方差分布

图 16(Fig. 16)

图 16 CDP1856处地震角道集 a实际角道集; b合成角道集; c实际角道集与合成角道集的残差

图 17和图 18分别是采用最小二乘法和本文方法进行反演得到的结果。可以看出, 当岩石物理模型和地震角道集都含有脉冲噪声时, 最小二乘法的反演结果并不理想。一方面, 岩石物理模型的脉冲噪声使反演得到的弹性参数偏离初始模型的趋势; 另一方面, 地震角道集的脉冲噪声使反演得到的物性参数变得不可靠。本文方法通过设置加权矩阵自适应调整地震角道集、岩石物理约束和先验信息之间的相对权重, 在一定程度上能够压制脉冲噪声的影响, 提高算法的稳定性。

图 17(Fig. 17)

图 17 CDP1856处最小二乘法反演结果 a纵波速度; b横波速度; c密度; d孔隙度; e含水饱和度; f泥质含量

图 18(Fig. 18)

图 18 CDP1856处本文方法反演结果 a纵波速度; b横波速度; c密度; d孔隙度; e含水饱和度; f泥质含量

沿着H₈层位曲线, 弹性参数与物性参数的初始模型和先验模型可以通过6口井的内插、外推模拟得到。图 19是利用本文方法反演得到的纵波速度、横波速度剖面, 图 20是利用本文方法反演得到的密度、孔隙度、含水饱和度和泥质含量剖面。可以看出, 弹性参数和物性参数的反演结果稳定性较好, 且反演结果与测井曲线较为吻合。

图 19(Fig. 19)

图 19 纵波速度(a)和横波速度(b)反演结果

图 20(Fig. 20)

图 20 密度(a)、孔隙度(b)、含水饱和度(c)和泥质含量(d)反演结果

储层弹性和物性参数同步反演是提取岩性参数的一条重要途径。随着勘探目标由构造油气藏转向岩性油气藏, 同步反演的应用越来越广泛。本文假设地震噪声信号和岩石物理模型误差服从高斯分布且带有加性脉冲噪声, 利用Hampel三截尾函数构建加权矩阵对反演目标函数进行约束, 通过反演残差自适应调整目标函数中的加权矩阵, 提高了同步反演算法对高斯噪声和脉冲噪声的适应能力。

当观测数据不完全服从高斯分布时, 多源信息约束反演本身是严重不稳定的。利用Hampel三截尾函数构建加权矩阵能在一定程度上自适应调节各个目标函数之间的权重, 消除地震和岩石物理噪声的统计性误差造成的反演不稳定性, 最终得到稳定的储层弹性和物性参数反演结果。