0 引言

现今多分量地震记录采集与处理技术日趋成熟，多分量地震数据解释已不再仅仅追求精细描述构造。在多波反演技术的帮助下，岩性反演、流体描述和断裂描述以及各向异性解释等已取得较大进展，然而，相比于多波地震数据采集和处理方法，多波地震数据解释方法进展较慢，难以体现多波地震技术的优越性^[1-6]。主要原因为：现有转换波解释技术尚不完善；没有很好地应用多波记录蕴含的地震信息；难以构建精确的多波地震模型；欠缺分析转换波数据的岩石物理依据；多分量地震资料与测井资料联合反演、AVO分析等方面存在诸多问题^[7-10]，需要深入探讨理论、方法。

地震属性是指分析或运算地震数据而获得的突显地层地震响应特征的数据。地震属性能揭示蕴含在常规地震数据中的隐晦信息，如，地层含流体时地震振幅变强、频率变低，地层存在裂缝时地震波相干性变弱等^[11-12]。基于地质因素差异，如果通过研究独立属性分析地质问题，往往因地震数据分辨率、信噪比等的限制而会出现偏差。振幅、频率和相位是最重要的地震属性，运用这些属性能够重建地震记录。然而，对依据纵波地震属性基于多波地震资料建立有效的多波属性解释方法，为高质量储层评价与流体识别提供有力保障，仍缺乏足够重视^[13-14]。

随着计算机网络与地震探测方法的迅速发展，可充分提取和利用地震记录中潜藏的信息，日益增多的地震属性提取方法也为地球物理解释提供了新工具或手段。如，通过图像分割对多类属性开展聚类分析从而优选最佳属性、以多谱成像原理为基础融合多属性，这些算法的共性在于基于非线性统计方法以及综合多属性信息提取最优属性^[15-17]。实际应用证明，非线性方法比其他方法的鲁棒性更强，多属性融合方法考虑了各种属性的异常信息，可充分体现地下地质体特征。当前，地震属性提取手段日渐丰富，属性提取算法正逐步转向非线性领域，纵波多属性分析（模糊聚类、深度学习）已广泛用于地震流体识别，将来地震属性分析落脚点将逐渐向地质方向靠拢^[18]。常规地震属性提取与分析多基于纵波地震资料，并未综合利用储层纵、横波信息。此外，常规属性求解方法只截出纵波波形的一种或几种特性，由于地球物理的多解性问题，仅依赖一种地震属性难以准确预测储层^[19-20]。基于多波地震数据有效识别地层岩性和含油气性，需要从地震信号中求取敏感地震属性，由于不同地层的敏感地震属性不同，因此在实际应用中需融合多种属性才能突出储层特性。

基于上述研究背景，本文提出了一种多尺度属性融合技术，以增强多波复合属性的油气敏感性。对于实际二维地震剖面，首先，构造一个高斯金字塔，以生成不同分辨率的各种复合属性。其次，融合不同尺度的所有属性形成有效的多尺度增强属性。再次，借鉴图像叠加理论，将不同的多尺度增强属性混合叠加，可有效突出储层发育区域，且保留了各种属性间的差异，能更好地描述砂体的含油气性。

1 多波地震复合属性构建原理

振幅、频率和相位是最基本的地震属性，利用这些属性能够重建地震记录、估算地层厚度等。本文利用希尔伯特变换获得多波瞬时振幅、瞬时相位和瞬时频率等地震特征属性^[21]，通过分析多波模型与实际数据的作用机理，研究油、水层的差异。具体表现为：当地震波经过油气储层时，纵波振幅增强、频率降低；转换波振幅基本不变、频率降低；瞬时相位不受地震信号强度的影响，能够清楚地反映弱反射信号，可获得更多的深层结构信息。此外，利用瞬时相位对时间的导数识别地震反射界面的效果较好^[14]。为提高多波属性的应用效果，本文提出了一种基于多波地震资料的新复合属性，该属性与振幅呈正相关、与频率呈负相关，不但与地层厚度存在较好的线性映射关系，满足厚度预测的一般要求，且通过指数运算，可放大多波属性间的差异，增强了多波复合属性识别油气的能力。值得注意的是，在求解多波复合属性时，由于传统的瞬时频率多呈无地质意义的负值，且多波的瞬时频率差异并不明显，故采用最小二乘法求解多波瞬时频率，放大了纵波与转换波的频率差异。首先依据实际井旁道数据，基于希尔伯特变换求解多波三瞬属性，分析储层（1.09~1.11 s）和非储层的差异。

实际地震记录$ \boldsymbol{X}\left(\tau \right) $的希尔伯特变换为

$ \boldsymbol Y\left(t\right)=\frac{1}{{{\rm{ \mathsf{ π} }} }}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\frac{\boldsymbol{X}\left(\tau \right)}{t-\tau }\mathrm{d}\tau $

(1)

式中：$ t $为时间；$ \tau $为时间延迟。$ \boldsymbol Y $$ \left(t\right) $的解析函数$ \boldsymbol Z\left(t\right) $、瞬时振幅$ \boldsymbol A\left(t\right) $、瞬时频率$ f\left(t\right) $和瞬时相位$ \boldsymbol \varphi \left(t\right) $分别为

$ \left\{\begin{array}{l}\boldsymbol Z\left(t\right)=\boldsymbol A\left(t\right){\mathrm{e}}^{\mathrm{i}t}\\ \boldsymbol A\left(t\right)=\sqrt{{\boldsymbol X}^{2}\left(t\right)+{\boldsymbol Y}^{2}\left(t\right)}\\ \boldsymbol \varphi \left(t\right)=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\left[\frac{\boldsymbol Y\left(t\right)}{\boldsymbol X\left(t\right)}\right]\\ \boldsymbol f\left(t\right)=\frac{\mathrm{d}\boldsymbol \varphi \left(t\right)}{\mathrm{d}t}\end{array}\right. $

(2)

图 1为实际井旁道的三瞬属性。由图可见，在储层位置（红色椭圆），PP波和PS波的瞬时振幅存在明显差异（图 1a），瞬时频率的差异不明显（图 1b），瞬时相位很稳定，不受油气影响（图 1c）。为突出多波三瞬属性的差异，引入基于整形正则化算子的瞬时频率求解算法^[22]，假定连续信号$ \boldsymbol X\left(t\right) $的复地震道为

$ \boldsymbol{S}\left(t\right)=\boldsymbol{X}\left(t\right)+\mathrm{i}\boldsymbol{h}\left(t\right)=\boldsymbol{A}\left(t\right){\mathrm{e}}^{\mathrm{i}\boldsymbol \varphi \left(t\right)} $

(3)

式中$ \boldsymbol{h}\left(t\right) $是$ \boldsymbol{X}\left(t\right) $的希尔伯特变化。

$ \boldsymbol{f}\left(t\right) $是$ \boldsymbol{\varphi }\left(t\right) $的时间变化率，即

$ \begin{array}{l}\boldsymbol{f}\left(t\right)={\boldsymbol{\varphi }}^{\mathrm{\text{'}}}\left(t\right)=\mathrm{I}\mathrm{m}\left[\frac{{\boldsymbol{S}}^{\mathrm{\text{'}}}\left(t\right)}{\boldsymbol{S}\left(t\right)}\right]\\\ \ \ \ \ \ \ =\frac{\boldsymbol{X}\left(t\right){\boldsymbol{h}}^{\mathrm{\text{'}}}\left(t\right)-{\boldsymbol{X}}^{\mathrm{\text{'}}}\left(t\right)\boldsymbol{h}\left(t\right)}{{\boldsymbol{X}}^{2}\left(t\right)+{\boldsymbol{h}}^{2}\left(t\right)}=\frac{\boldsymbol{l}\left(t\right)}{\boldsymbol{o}\left(t\right)}\end{array} $

(4)

考虑地震信号的连续相位，使用最小二乘法进一步求解$ \boldsymbol{f}\left(t\right) $

$ \boldsymbol{f}\left(t\right)={\left[{p}^{2}\boldsymbol{I}+{\left(\boldsymbol {DO}\right)}^{\mathrm{T}}\left(\boldsymbol {DO}\right)\right]}^{-1}{\left(\boldsymbol {DO}\right)}^{\mathrm{T}}\boldsymbol{I} $

(5)

式中：I和O分别为l(t)和o(t)组成的矩阵，I为单位矩阵；$ D $为整形算子；p为权系数，通常取$ \boldsymbol O $中数据最大值的1%~5%。

图 2为基于整形正则化算子求解的瞬时频率。由图可见，在储层位置（红色椭圆）PP波和PS波的瞬时频率差异较非储层位置明显，表明多波属性能够有效突出储层特征。

基于此，假定复合属性γ₁为

$ {\gamma }_{1}={A}_{\mathrm{h}}{{\mathrm{e}}^{-}}^{\alpha {f}_{\mathrm{h}}} $

(6)

式中：$ {A}_{\mathrm{h}} $为多波资料的瞬时振幅；$ {f}_{\mathrm{h}} $为多波资料的瞬时频率；α为频率控制系数。由式（6）可知，$ {\gamma }_{1} $与多波振幅呈正相关、与频率呈负相关，且与地层厚度存在较强的线性映射关系，满足厚度预测的一般要求，通过指数运算，放大了多波属性间的差异，增强了多波复合属性识别油气的能力。

图 3为A井、B井井旁道的γ₁。由图可见，γ₁对储层较敏感，在储层位置（红色椭圆）PP波和PS波属性差异远高于非储层位置。

实际多波资料的作用机制表明：当地震资料经过油气储层时，PP波振幅增强、频率降低；PS波振幅基本不变、频率降低；瞬时相位不受地震信号强度影响，能够清楚地反映弱反射信号，可获得更多的深层结构信息。为此，将γ₁乘以瞬时频率f(t)($ \boldsymbol{f}\left(t\right) $的标量形式)得到复合属性$ {\gamma }_{2} $

$ {\gamma }_{2}={\gamma }_{1}f\left(t\right) $

(7)

γ₁通过求逆的方式求解$ {f}_{\mathrm{h}} $，并引入正则化算子平滑，相当于削弱了底层的影响（瞬时相位的响应），并通过指数形式表征整体地震响应。基于瞬时相位求导获得的瞬时频率f(t)突出了结构信息，其提供的较高分辨率放大了底层的弱反射信息，能够得到不同尺度的多波属性，可更好地对比多波属性间的差异。

图 4为A井、B井井旁道的γ₂。由图可见，γ₂的分辨率较高，清晰地指示了储层位置（1.11 s）。且储层位置PP波的γ₂远高于非储层位置。

2 基于金字塔的多尺度增强技术

在图像处理中，多尺度分析方法是指以不同分辨率级别表示图像，每个分辨率级别都包含想要缩放的属性特征。该方法已广泛用于图像处理，如图像编码、压缩、分析和检测。金字塔表示是多尺度分析的经典方法之一，它产生了尺度空间理论和其他多尺度信号表示。高斯金字塔与拉普拉斯金字塔是典型的多尺度分解算法^[23-24]，本文利用高斯金字塔挖掘多波复合属性的不同尺度特征。针对二维地震剖面$ {\boldsymbol{S}}_{}\left(M,N\right)\in {\boldsymbol{R}}^{M\times N} $，S是尺寸为$ M\times N $的矩形支撑，$ {\boldsymbol{R}}^{M\times N} $是尺寸为$ M\times N $的实数空间集。首先设置金字塔的尺度为0，表示全分辨率尺度；然后，利用高斯核模糊$ {\boldsymbol{S}}_{0}\left(M,N\right) $得到尺度1——$ {\boldsymbol{S}}_{1}\left(M,N\right) $，通过下采样2倍计算$ {\boldsymbol{S}}_{1}\left(M,N\right) $，其余尺度以类似方式产生。高斯金字塔构造公式为

$ {\boldsymbol{S}}_{i}(M,N)=\sum\limits_{\forall k,l\in \boldsymbol{G}}g(k,l){f}_{i-1}(2M+k,2N+l) $

(8)

式中：$ i=\mathrm{1,2},\cdots n,n $是高斯金字塔分解的尺度；$ g(k,l) $为定义在区域G的高斯核，即

$ g(k,l)=\frac{1}{2\mathrm{{\rm{ \mathsf{ π} }} }{\sigma }^{2}}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{{k}^{2}+{l}^{2}}{2{\sigma }^{2}}\right)\ \ \forall (k,l)\in \boldsymbol{G} $

(9)

式中：$ k,l $是矩形尺寸；$ {\sigma }^{2} $是平滑因子，一般取固定值。

对于高斯金字塔，输入地震剖面的每个样本将有K个相应的属性样本，即每个尺度对应一个样本。基于算数平均或加权平均方法融合K个样本获得一个代表性的样本。值得注意的是，尺度为0的属性图表征由传统方法计算的属性，其中金字塔中较高尺度的属性将具有较低的分辨率，同时保留了大规模特征。另一方面，金字塔中的较低尺度产生的属性细节更丰富但存在噪声。为了生成单一的增强属性，文中设计了一种融合机制，该机制要求属性映射大小相同。因此，在融合之前，通过插值将所有映射大小调整为尺度为0。与传统的空间滤波方法相比，多尺度融合步骤可以解释为在尺度空间进行滤波。如，想得到高分辨率尺度特征，可以使用较大权重的加权平均融合方案降低金字塔尺度，还可以使用非线性算子进行多尺度融合，如秩滤波器。图 6为多尺度属性增强流程。

3 多波属性联合渲染技术

最初人们用黑色和白色描述地震记录，因此难以看清很多地震曲线。为此，Balch^[25]在1971年首次使用彩色工具描绘地震资料，可使人们看到更多的地震细节信息。从此，地震属性彩色展示成为地球物理解释的有力工具，但当时的彩色展示只能逐一描述多个地震属性，无法清晰地体现总体趋势，且难以有效突出地下特殊异常体的存在^[26]。为更好地利用地震资料中的隐晦信息，人们引入RGB理论对地震资料颜色渲染，可有效地突出储层地质异常^[27]。本文提出的复合属性主要依赖于地震数据的振幅、频率，数据的量纲保持大体一致，在此基础上，将基于多波的多尺度属性映射至RGB颜色空间，并借鉴图像分析思路叠加显示，可有效地突出储层发育区域。

基于RGB理论的颜色渲染技术是将输入的地震属性通过映射变换到RGB空间，其中的每个值对应颜色空间中的数值

$ {I}_{\mathrm{R}\mathrm{G}\mathrm{B}}\left(\boldsymbol{x}\right)=R\left[{I}_{\mathrm{R}}\left(\boldsymbol{x}\right),{I}_{\mathrm{G}}\left(\boldsymbol{x}\right),{I}_{\mathrm{B}}\left(\boldsymbol{x}\right)\right] $

(10)

式中：$ {I}_{\mathrm{R}\mathrm{G}\mathrm{B}}\left(\boldsymbol{x}\right) $为变换后的颜色值，x=(x, y, z)为三维地震记录的某一点；$ R $[•]为RGB变换算法，$ {I}_{\mathrm{R}} $、$ {I}_{\mathrm{G}} $、$ {I}_{\mathrm{B}} $为红、绿、蓝三基色，通过红、绿、蓝的不同分量表征与坐标轴距离，能展示地震资料的具体颜色特征，即

$ \begin{array}{l}{\boldsymbol{I}}_{\mathrm{R}\mathrm{G}\mathrm{B}}\left(\boldsymbol{x}\right)=\left\{R\left({I}_{\mathrm{R}},{I}_{\mathrm{G}},{I}_{\mathrm{B}}\right)\left|0\le {I}_{\mathrm{R}}\le 255,\right.\right.\\ \ \ \ \ \ \ \ \ \ \ \left.0\le {I}_{\mathrm{G}}\le \mathrm{255,0}\le {I}_{\mathrm{B}}\le 255\right\}\end{array} $

(11)

$ {I}_{\mathrm{R}\mathrm{G}\mathrm{B}}\left(\boldsymbol{x}\right) $综合利用了多波复合属性的细节，与常规单个地震属性相比，地质特征更鲜明、细节更丰富。

对于色彩饱和度，本文基于最大和最小值法，将提取的多波属性归一化

$ {C}_{q}\left(\boldsymbol{x}\right)=255\times \frac{{R}_{q}\left(\boldsymbol{x}\right)-{{{R}_{q,\min}}}}{{{R}_{q,{\max}}}-{{R}_{q,}}_{\min}} $

(12)

避免了多波属性赋值差异造成的一些有效响应信息被隐藏或丢失。式中：$ {C}_{q}\left(\boldsymbol{x}\right) $为归一化多波属性，q =R，G，B表示红、绿、蓝；$ {R}_{q}\left(\boldsymbol{x}\right) $为提取的多波复合属性，$ {{{R}_{q,\min}}} $与$ {{R}_{q,{\max}}} $分别为$ {R}_{q}\left(\boldsymbol{x}\right) $的最小值与最大值。

将$ {C}_{\mathrm{R}} $、$ {C}_{\mathrm{G}} $、$ {C}_{\mathrm{B}} $映射至RGB空间，可分别反映同一层面的检测信息。基于对应像素坐标不发生偏移的原理，本文直接进行叠加融合（图 7），即

$ I=0.229{C}_{\mathrm{R}}+0.587{C}_{\mathrm{G}}+0.114{C}_{\mathrm{B}} $

(13)

4 模型试算

为了验证基于多波地震的多尺度属性融合方法的效果，依据A区实际地质数据构建了复杂地质模型（图 8），可见模型内围岩速度与最深处含油气砂岩速度相差不大。图 9为基于数学理论构建的多种复合属性，其中多波融合地震记录是通过RGB理论融合所得，具体方法参见文献[26]。图 10为基于高斯金字塔的图 9f、图 9i增强，图 11为由RGB得到的多波融合属性与叠加属性。由图可见：基于振幅、频率差异构建的多种复合属性从不同角度描述了真实储层的发育位置，为描述深层含油气砂岩提供了有力依据，为多波联合反演奠定了基础（图 9）；高斯金字塔技术进一步增强了多波复合属性（图 10）；基于RGB联合渲染技术得到的多波联合属性识别的储层特征与图 8基本一致，但是由于未考虑实际噪声影响及多波匹配差异等，因此最终融合效果不好（图 11）。

为对比、分析多波属性融合方法的效果，抽取模型的单道属性结果进行对比（图 12）。可见，在储层发育处（红色方框处）多波融合属性（图 12e）显著突出了储层差异，PP波属性和PS波属性较一致，图 12e与PP波、PS波属性整体差异不大，具体表现为：①由γ₁（图 12a）及γ₁增强（图 12b）清晰地展示了PP波和PS波响应的较大差异，但未能明显突出储层发育区域，多波复合属性融合了PP波和PS波的差异，突出了储层响应（图 12b）。②γ₂（图 12c）及γ₂增强（图 12d）有效突出了地层细节信息（蓝色线框处形成了完整的地震界面）。③图 12e充分考虑了PP波和PS波的差异，放大了储层响应特征，且具有较高的分辨率。图 9~图 12表明，对于含油气砂岩与围岩差异较小的复杂二维模型，依据RGB的多波地震属性融合能够检测储层，并可识别反射界面位置。

5 实际数据试算

四川盆地致密气藏储量大，勘探潜力大，沙溪庙组致密气藏整体勘探、开发效果好。通过分析可知，A区以中—细砂岩为主，重点发育河流—湖泊沉积体系，平面上砂体纵向交错，横向连续分布，为生成油气储层夯实了牢靠的物质基础。A区主要目的层（1.11 s附近）由河道砂体叠置而成，具有低PP波阻抗、高PS波阻抗的特征。

图 13~图 16分别为多波地震记录、基于多波地震记录得到的多种复合属性、基于高斯金字塔的多波γ₁、多波γ₂增强以及基于RGB的多波融合属性与叠加属性，其中图 16是由RGB理论融合所得^[26]。由图可见：①仅仅依据多波振幅差异难以识别砂体（图 13c），且图 13c的砂体振幅弱于PP波记录（图 13a）。②多波γ₁（图 14c）、γ₂（图 14f）有效识别了储层发育位置，且清晰地刻画了河道边界，此外γ₁、γ₂的横向连续性很好，能够较好地描述河道展布。③对比多波γ₁ (图 15a)、多波γ₂ (图 15b)增强与图 14c、图 14f可见，图 15a、图 15b有效增强了多波复合属性，突显了储层差异。④RGB融合技术整体考虑了多波复合属性性质，可突出河道砂体信息，RGB叠加技术综合利用多波复合属性差异（图 16b），可更好地预测储层，并为确定水平井位置提供了参考依据。

图 17为实际单道多波属性对比。由图可见，在储层发育位置（红色方框处）多波融合属性（图 17e）显著突出了储层差异，具体表现为：①PP波γ₁和PS波γ₁差异较大（图 17a），而多波复合属性能够融合纵波和转换波的差异，可有效突出储层差异；②多波复合属性在储层发育位置响应明显，并通过增强压制了非储层影响，突出了储层响应（图 17b）；③γ₂（图 17c）及γ₂增强（图 17d）有效突出了地层细节信息（蓝色线框处形成了完整的地震界面）；④图 17e充分考虑了PP波和PS波的差异，放大了储层响应特征，且具有较高的分辨率，并极大地压制了非储层影响。

6 结论

针对传统地震属性提取方法基本依托于纵波资料，缺乏多波资料属性提取技术的理论依据，为此，本文基于多波记录的波形差异分析和岩石物理研究，联合信号处理手段，基于数据驱动的多波融合属性检测油气。首先依据实际井旁道数据，基于希尔伯特变换求解多波三瞬属性，分析储层和非储层的差异。综合图像检测手段，开发了基于多波地震的多种复合属性计算方法，提出了一种基于多波地震的多尺度属性融合技术，在拓展多波分析手段的同时，准确提取了多波反射信息，降低了纵波油气预测中的多解性。实际数据试算表明，所提方法能分辨靶区河道砂体的储层分布，获得了较好效果。

尚需指出，文中方法存在以下局限性：①针对多波资料，需进行高精度多波资料匹配，以保证多波复合属性计算结果的精度；②研究区的多波复合属性与多波振幅呈正相关、与频率呈负相关，当复合属性与振幅或频率呈负相关时，应适当调整复合属性公式。