0 引言

火成岩是一种特殊的岩性体，具有高速度、非均质、原生孔隙度大等特点。火成岩与围岩有多种接触关系，火成岩反射波场和转换波场较复杂，致使火成岩发育区的地震波场混乱，对于下伏围岩的地震响应特征具有很强屏蔽作用，导致火成岩下伏储层反射信号能量较弱，严重制约了火成岩发育区的储层预测。

目前去除地震强屏蔽方法多集中在多子波分解与重构技术。对于信号分解，匹配追踪(Matching Pursuit, MP)算法在信号的稀疏表示方面效果很好，是多子波分解的核心算法。Mallat等^[1]利用基于Gabor原子形成的完备原子库，通过不断迭代，求取原信号与原子库子波的最佳稀疏表示——MP算法；Chakraborty等^[2]首次把MP算法应用于地震信号分析，其中MP原子库(即子波类型和子波分辨率)的选择十分重要。该项技术主要是利用不同频率、振幅、相位的各类原子与地震道匹配^[3-8]，分解成一个原子库，其中最大能量原子对应于强背景，在重构中舍去，从而达到去屏蔽的效果。

多子波分解与重构技术往往存在以下两个问题^[9]：首先，对于分解的子波实际意义不明确，地震剖面上最强子波不一定都代表火成岩，也可能是其他岩体；其次，地震道子波分解是对单道处理，无法考虑地层横向关系，对于复杂的火成岩发育区，火成岩存在形态和厚度变化，反射特征并不一致，容易造成匹配过度和重构剖面连续性差。

金成志等^[9]在2017年首次结合Wheeler变换和主成分分析(Principal Component Analysis, PCA)技术去除地震强屏蔽，可了解强屏蔽层之下的砂体展布特征，为去除强屏蔽提供了新思路。长、短旋回波形特征的主要区别在于相位和频率，因此，采用PCA提取反射波形的空间一致性相位信息分离长、短旋回，可去除地震强屏蔽的影响。由于“长旋回反射相位特征一致”的认识仅在相对地质年代域中成立，所以在利用PCA分解长、短旋回波形时，需要进行Wheeler变换。

De Groot等^[10]在2006年正式提出三维Wheeler域变换技术；De Bruin等^[11]于2007年提出Wheeler域与时间域地震数据交互解释的技术思路。至此，Wheeler域变换技术被正式推广至油气勘探领域。目前，应用地震数据构建Wheeler数据体的方法主要有三种^[12]：第一种是由De Groot等^[10]、Nordlund等^[13]提出的基于地震层位解释的Wheeler数据体生成方法，通过追踪尽可能多的层位，解释层序边界构建Wheeler数据体。但利用人工拾取同相轴构建Wheeler数据体效率低、工作量大。第二种方法是Lomask^[14]、Fomel^[15]提出的将数据体拉平进而生成Wheeler数据体的方法，其优势是不需要层位拾取，但是在存在不整合面和断层等的复杂地质情况下，无法正确拉平同相轴。第三种方法由Stark^[16-17]提出，首先通过地震瞬时相位解缠生成相对地质时间(RGT)体，然后在RGT体的基础上构建Wheeler数据体。这种方法的精度取决于瞬时相位解缠的精度，难于取得理想的应用效果。

由于火成岩相带类型多、地震波振幅横向变化快、反射特征复杂，因此在利用长、短旋回波形分析法去除火成岩强屏蔽时，往往存在以下问题^{[9, 12]}：①需要对火成岩发育区多个层序界面进行精细层位解释，以形成准确的地层切片体^[12]；②若Wheeler域中火成岩反射同相轴无法拉平，通过PCA提取的长旋回波形存在一定偏差，从而影响长、短旋回分离后的短旋回波形的有效精度；③局部区域发育的火成岩爆发相、喷溢相的界面往往与长旋回界面不平行，难以准确剥离和剔除火成岩强反射；④在地下断层发育、地层出现尖灭、同相轴连续性差等复杂情况下，经Wheeler正、反变换后，重构的地震信号可能出现局部失真，并且重构过程较复杂。

为了剥离火成岩强反射，可根据解释的火成岩层位信息设置相应的时窗宽度，通过搜索纵、横向地震波形特征准确获得火成岩层位，实现局部层拉平，再利用子体波形PCA剥离火成岩强反射。基于此，本文采用三维多项式曲面拟合代替Wheeler变换实现局部层拉平，形成了基于三维多项式曲面拟合的子体波形PCA火成岩强屏蔽剥离技术，避开了Wheeler变换带来的问题。同时，在剥离过程中无需引入振幅阈值控制，对于火成岩反射振幅横向变化快的问题，只要在设置的子体窗口范围内火成岩反射振幅变化不大，就可以提取该子体窗口的火成岩强反射，并通过滑动子体窗口，在横向剥离火成岩强反射，提高了方法的实用性。

1 三维多项式曲面拟合局部层拉平

首先对解释的火成岩层位进行三维多项式曲面拟合，通过空间搜索获取精确的层位，从而有效解决Wheeler变换中的问题①；针对问题②和③，如果存在局部发育的火成岩相带，无论其界面是否与本地区长旋回界面平行，根据解释的火成岩层位信息，设置相应的时窗宽度，通过搜索纵、横向波形特征拟合火成岩的准确层位，从而实现局部层拉平，最终基于子体波形PCA分解方法剥离火成岩强反射，从而有效避开了问题④。具体过程如下。

利用三维地震数据体拟合三维多项式曲线时，需要选取一个分析区域，首先划分一个初始窗口，然后假设某一时窗W的时间取值范围为t∈[-l, l]，横向窗口范围为x∈[-m, m]、y∈[-n, n]，对该时窗沿x、y方向建立离散坐标系

$ \begin{gathered} \boldsymbol{D}=\{(x, y) \mid x \in[-m, m], \\ y \in[-n, n] ; x, y \in \bf{Z}\} \end{gathered} $

(1)

式中：x、y分别为Crossline、Inline方向的坐标；Z为整数集合；D为研究区的二维空间。

时窗W的中心点时间拟合多项式为

$ T(x, y)=\sum\limits_{i=0}^{5} C_{i} G_{i}=C_{0}+\sum\limits_{i=1}^{5} C_{i} G_{i} $

(2)

式中：T(x，y)为由拟合得到的点(x，y)的强反射同相轴时间；C_i为第i个正交多项式的拟合系数；G_i为D的二元二次正交多项式，即

$ \left\{\begin{array}{l} G_{0}=1 \\ G_{1}=x \\ G_{2}=x y \\ G_{3}=x y \\ G_{4}=y^{2}-\frac{1}{3} m(m+1) \\ G_{5}=y^{2}-\frac{1}{3} n(n+1) \end{array}\right. $

(3)

根据多道互相关性最强的原则确定最佳C_i。对W内数据S(x，y，t)可以计算归一化多道互相关系数，即

$ R\left(C_{0}, C_{1}, \cdots, C_{5}\right)=\sum\limits_{t=-l}^{l} \frac{A-B}{C D} $

(4)

其中

$ \left\{\begin{array}{l} A=\left\{\sum\limits_{x=-m}^{m} \sum\limits_{y=-n}^{n} S\{x, y, [T(x, y)+t]\}\right\}^{2} \\ B=\sum\limits_{x=-m}^{m} \sum\limits_{y=-n}^{n} S^{2}\{x, y, [T(x, y)+t]\} \\ C=(2 m+1)(2 n+1)-1 \\ D=\sum\limits_{t=-l}^{l} \sum\limits_{x=-m}^{m} \sum\limits_{y=-n}^{n} S^{2}\{x, y, [T(x, y)+t]\} \end{array}\right. $

(5)

为了不影响拟合效果，求取其他系数时，需要假设C₀的值为0，可以确保拟合的信号时间固定在窗口中心。

以C₁的扫描为例，选取初始窗口，以初始窗口边界的样点时间作为扫描范围，通过不同的样点时间得到不同的C₁；由式(4)求出相应的相关系数，其中最大相关系数对应的C₁即为最终系数。由于G_i为二元二次正交多项式，故可以用独立扫描的方式求得C_i，从而简化计算。具体过程为：①固定C₀(可取为0)，先扫描C₁，此时C₂~C₅均为0，比较多道互相关值确定C₁；②保持已确定的多道互相关系数，再扫描C₂，依此类推。扫描过程中窗口的形状随多项式系数的变化而变化，求取系数的过程就是搜索信号的过程，由扫描确定的窗口就是期望信号的窗口。

2 子体波形PCA分解剥离火成岩强屏蔽

PCA的核心思想是将数据降维，它能够有效地提取数据中的“主要”信息，降低了复杂数据的维度，从而去除冗余成分，挖掘数据中的隐藏信息。PCA方法原理简单，已用于众多领域。因为PCA能将一组相关性较强的变量转化为一组相关性较弱的矩阵，最大特征值对应的即为主要成分，而小特征值对应的往往是相关性较差的噪声。因而PCA常用于去除随机噪声和信息聚类，本文采用基于三维多项式曲面拟合局部层拉平数据体，火成岩反射层位经过层拉平，在横向上波形特征基本一致，通过PCA将其作为背景信息提取。

具体实现过程为：首先以一个地震道为中心，取若干相邻道形成一个PCA单元，再通过层位时间信息时窗约束进行层拉平；其次通过PCA波形分析提取局部层拉平后的相似波形，将此波形从PCA单元中剔除，达到剥离火成岩强屏蔽的目的。上述基于子体波形的PCA方法称为子体波形PCA分解方法。具体原理阐述如下。

设X =(s₁，s₂，…，s_q)为一个p×q的子窗口内的地震数据矩阵，其中s_i为某道地震数据，子窗口内的道数为q，每道的时间样点数为p。经过线性变换将原始数据X转换到另一个变量空间Y，Y中的各变量彼此不相关，即

$ \boldsymbol{Y}_{r \times q}=\boldsymbol{K}_{r \times p} \boldsymbol{X}_{p \times q} $

(6)

其中r＜p为Y的维度。式(6)表示原始地震数据的“提炼”过程，通常利用奇异值分解(Singular Value Decomposition，SVD)精准地提炼有用信息，获得矩阵K _r×p，即

$ \boldsymbol{X}_{p \times q}=\boldsymbol{V}_{p \times p} \mathit{\pmb{\Sigma}}_{p \times q} \boldsymbol{U}_{q \times q}^{\mathrm{T}} $

(7)

式中：V和U分别为左奇异矩阵和右奇异矩阵，并且V^TV = U^TU = I；Σ_p×q为特征值λ_g(g=1，…，q)组成的对角矩阵，且λ₁>λ₂>…>λ_q，其中特征值越大，包含的数据信息特征越多。此时如果选取前r个特征值重构，便得到$\boldsymbol{K}_{r \times p}\left(\boldsymbol{K}_{r \times p}=\boldsymbol{V}_{p \times r}^{\mathrm{T}}\right)$

$ \hat{\boldsymbol{X}}_{p \times q}=\boldsymbol{V}_{p \times r} \boldsymbol{Y}_{r \times q} $

(8)

式中$\hat{\boldsymbol{X}}_{p \times q}$即为重构数据，保留了信号的主成分。

3 模型测试及实际资料应用

图 1为惠民凹陷YHM地区火成岩地质模型及其正演模拟记录。由图可见：火山岩溢流相厚度小，横向分布范围广，产状近似水平(图 1a)，顶、底形成强反射，波形、振幅基本保持不变(图 1b)；火山岩通道相表现为空白反射或杂乱反射(图 1b)；火山岩爆发相、喷溢相的横向厚度变化较大，中间厚，向两边逐渐变薄(图 1a)，多个界面的反射叠合在一起形成复波，同相轴呈凸起状，并且淹没了下伏地层反射(图 1b)。

图 2为图 1b的Wheeler正变换结果。由图可见，Wheeler正变换的时间段为100~600ms，主要包含5组同相轴。在Wheeler域通过基于PCA分解方法提取强振幅、相位特征相同的并与长旋回对应的子体波形，在此基础上再进行Wheeler反变换，得到火成岩强反射(图 3a)，最终得到剔除两套火成岩强反射的结果(图 3b)。可见：①在Wheeler域火山岩溢流相和火山岩沉积相界面下方同相轴被拉平，但火山岩爆发相、喷溢相反射同相轴并没有完全拉平，因此剥离的火山岩强反射存在一定偏差，影响剥离火成岩强屏蔽效果(图 3a)。②火山岩溢流相横向分布范围广，产状近似水平，因此被很好地剔除，只是在断层附近有一定残余；火成岩爆发相和喷溢相的位置仍存留一定剩余波形，同时在断层左侧火成岩不发育位置提取了假火成岩反射，说明仅通过振幅阈值约束往往难以判定火成岩反射的可靠性，需要解释人员结合研究区地质条件分析进行人工干预(图 3b)。

图 4为三维多项式曲面拟合得到的火成岩强反PCA子体的纵向时窗宽度为60ms，横向窗口以分析道为中心往左、右两边各取2道射层在时窗中心点处的时间剖面。由图可见，获得的火成岩反射层横向连续性强，并且在断层位置处较精准地提取了火成岩空间位置，在此基础上可对火成岩反射层进行局部层拉平。

通过分析子体波形PCA分解的各分量发现，在火成岩反射层局部拉平的子体窗口内，主要能量集中在最大特征值对应的分量上，且横向波形特征一致，该分量即为局部拉平后的火成岩强反射，而被火成岩淹没的其他反射层信号分布在其余特征值对应的分量中。因此，选取最大特征值对应的分量进行重构，可以剥离火成岩强反射(图 5a)，并得到剔除火成岩强反射的剖面(图 5b)。可见，本文方法剥离的火成岩强反射横向连续性强(图 5a)，能够较准确地剔除火成岩强反射界面的综合响应，保留火成岩下方的地层反射信息，从而更好地削弱火成岩对下伏地层反射波的屏蔽作用(图 5b)。

为了验证基于多项式曲面拟合子体波形PCA分解方法的实际应用效果，对惠民凹陷YHM地区的实际地震资料进行处理。该地区火成岩十分发育，火成岩顶、底界面较大的波阻抗差异导致出现强反射同相轴，强反射同相轴严重屏蔽了火成岩下方储层的弱信号，降低了储层预测精度。

W1井综合录井图(图 6a)显示：W1井T3界面上方火成岩为玄武岩，T3界面下方沙三段3砂组储层发育。W3井综合录井图(图 6b)显示：W3井T3界面上方火成岩为凝灰岩，T3界面下方沙三段3砂组储层不发育。由W1-W2-W3井连井地震剖面(图 7)可见，T3反射层(图 7绿色虚线)上方火成岩发育，为中强连续反射(图 7蓝框位置)，主要目的层沙三段3砂组(图 7黄色虚线)在T3反射层下方50ms处。

图 8为对图 7剥离火成岩强屏蔽的结果。由图可见：①多子波分解法剥离的火成岩强反射(图 8a)横向连续性较差，其使用的子波库为雷克子波库，剔除火成岩强反射剖面(图 8b)残留部分火成岩信息。②基于Wheeler正、反变换的长、短旋回PCA分解法在断层位置处剥离了火成岩强反射(图 8c)，并错误地识别了火成岩范围(蓝框区域)，导致剔除的火成岩强反射不准确(图 8d)；基于长、短旋回PCA分解法无法准确提取局部火成岩同相轴(红框区域)，如果扩大提取的时间范围，会损失有效信息。③三维多项式曲面拟合得到的火成岩强反射层在时窗中心点处的时间剖面(图 8e)在横向上光滑、连续，对应层位准确。④采用基于三维多项式曲面拟合的子体波形PCA分解法剥离的火成岩强反射考虑了火成岩发育区的地震反射特征，实现了局部层拉平，再进行子体波形PCA分解，准确地提取了火成岩强反射(图 8f)，剔除火成岩强反射剖面(图 8g)的横向连续性强，较准确地剔除了强反射界面的影响，从而在一定程度上增强了下伏地层反射。⑤利用基于时频分析的地震波能量衰减补偿方法^[18]对图 8g进行高频能量衰减补偿处理(图 8h)，可提高对地震弱反射信息的刻画能力。图 9为图 7蓝框区域及图 8h蓝框区域的局部放大，可见：W1井T3界面上方的火成岩强反射(图 9a红圈处)淹没了沙二段1砂组以及1砂组与火成岩之间的一套薄互层砂组(图 9a红色虚线框)，导致在原始地震剖面上无法识别这两套砂组；在剔除火成岩强反射及能量衰减补偿处理后，这两套砂组对应的同相轴均呈弱反射信号特征，能够进行横向追踪和识别(图 9b红色虚线框)，同时准确地刻画了3砂组，砂体边界向右延伸了175m(图 9中蓝色箭头处)。实际应用表明：本文方法能够削弱火成岩对下伏地层反射的屏蔽作用，提高对地震弱反射信息的刻画能力。

4 结论

本文采用三维多项式曲面拟合实现火成岩强反射局部层拉平，结合子波相位特性，利用子体波形PCA分解方法剥离火成岩强反射，避开了Wheeler变换带来的问题，形成了一种非常适用的火成岩强屏蔽剥离技术。

(1) 该项技术在剥离火成岩强反射过程中不需要引入振幅阈值的控制，通过分析研究区火成岩反射特征，合理设置子体窗口，有效解决了火成岩反射振幅横向变化快的问题；同时，通过子体窗口的连续滑动，解决了火成岩强反射横向剥离问题，克服了多子波分解方法和基于Wheeler变换的长、短旋回PCA分解方法横向连续性差的缺陷。

(2) 该项技术在实际应用过程中，通过对解释的火成岩层位进行三维多项式曲面拟合，利用空间搜索获取精确层位，自动生成局部层拉平数据体，为子体波形PCA分解提供了可靠的基础数据，简化了实现步骤，实用性大大增强。

(3) 模型测试和实际资料应用表明，该方法剥离的火成岩强屏蔽横向连续性强，能够较准确地剔除火成岩强反射界面的综合响应，在一定程度上削弱了火成岩对下伏地层反射波的屏蔽作用，提高了对地震弱反射的识别能力，为弱信号的能量补偿和后续处理奠定了基础。