0 引言

石油和天然气为重要的战略资源，随着国民经济的发展，油气需求量逐渐增大，油气勘探目标逐渐从构造油气藏转向岩性油气藏^[1]，薄层及薄互层逐渐成为勘探工作的重点。砂体厚度小、横向连通性差、难以确定沉积尖灭位置等因素增大了薄层和薄互层的识别难度^[2]，特别是当目标储层上、下存在高阻层时，较强的波阻抗差在地震剖面上呈强反射^[3]，屏蔽了目标砂层的弱反射信息，难以确定薄储层的尖灭线位置。地层尖灭线是控制油气成藏的重要因素，对于油气勘探、开发评价具有重要意义^[4]。在一定的地域范围内，通常地层的结构特征变化不大，在地震剖面上同相轴的横向连续性较好^[5]，其反射子波的振幅、频率、相位等属性保持相对一致。因此根据这种特征利用匹配追踪稀疏理论拾取和剥离强反射干扰，进而削弱强反射对薄层砂体有效反射信息的屏蔽作用。

Mallat等^[6]提出的匹配追踪算法(Matching Pursuit)根据信号局部结构特征构建超完备的字典原子库，在一组能够适应各种时频特性的原子上将信号展开，最终利用有限个原子稀疏表示原始信号，实现信号的自适应分解。Liu等^[7-8]提出了动态匹配追踪算法的概念，在每次迭代过程中参与內积运算的时频原子随信号的局部瞬时属性发生变化，提高了信号的分解效率。匹配追踪算法在各领域得到广泛应用。宋维琪等^[9]将复数子波匹配追踪算法应用于储层预测领域，可较好地识别薄层砂体信息；杨昊等^[10]利用匹配追踪算法自动提取薄层顶、底反射系数，定量预测薄层厚度；张繁昌等^[11]研究了基于匹配追踪瞬时谱计算的砂体尖灭线识别方法，并用于检测三角洲沉积砂体的尖灭线；李海山等^[12]采用贪婪匹配追踪算法分离煤层强反射；李坤等^[13]利用匹配追踪算法识别流体、刻画储层边界，减少了反演流程，且反演结果与测井资料保持一致。

基于以上研究背景，本文采用动态匹配追踪算法识别强反射，通过得到的最佳匹配强反射的子波剥离强反射。考虑到传统的瞬时频率运算结果存在负频率现象且易受噪声影响，即当计算结果出现负频率时，只能进行时频原子频率属性的全局搜索，并不是真正意义上的动态搜索^[14]。为此，本文引进局部频率计算模式用于约束动态匹配追踪算法最优原子的搜索范围，其频率计算结果合理，可以确定最佳匹配强反射的地震子波。通过构建含高阻层的地质模型和实际资料验证了方法的可行性。

1 匹配追踪强反射识别与分离方法 1.1 基于局部频率约束的动态匹配追踪算法

传统的匹配追踪算法是一种贪婪算法^[15]，通过创建超完备的时频原子库，对时移、频率和相位参数进行全局扫描，以寻找最佳相关的匹配原子，但该算法每一次迭代的內积运算量大，对信号的分解效率极低。动态匹配追踪算法将信号的瞬时特征作为先验信息约束匹配原子的搜索范围^[16]，以先验信息创建动态子波库，使每一次迭代的內积运算量大为减小，极大地提高了分解效率。传统方法以瞬时频率确定最优原子搜索范围，而瞬时频率是瞬时相位对时间的导数，计算结果存在“负频率”现象，且受噪声干扰大，计算结果极不稳定。当频率计算结果不合理时，只能对时频原子的频率属性进行全局搜索，并不是真正意义上的动态搜索。为此，本文将局部频率算法用于动态匹配追踪算法，根据整形光滑理论，利用差分算子综合某数据点邻域内的频率值，计算得到该数据点的频率属性信息，局部频率计算结果合理、稳定，抗噪性较强, 因此可以直接用于确定最优时频原子的搜索邻域。

连续信号x(t)的复地震道为^[17]

$ c(t) = x(t) + {\rm{i}}h(t) = A(t){{\rm{e}}^{{\rm{i}}\varphi (t)}} $

(1)

式中：h(t)为x(t)的Hilbert变换；A(t)为地震道包络；φ(t)为信号的瞬时相位。

传统的瞬时角频率ω(t)是瞬时相位φ(t)的时间变化率^[18]，即

$ \begin{array}{l} \omega (t) = {\varphi ^\prime }(t) = {\mathop{\rm Im}\nolimits} \left[ {\frac{{{c^\prime }(t)}}{{c(t)}}} \right] = \frac{{x(t){h^\prime }(t) - {x^\prime }(t)h(t)}}{{{x^2}(t) + {h^2}(t)}}\\ \;\;\;\;\;\;\;\; = \frac{{u(t)}}{{v(t)}} \end{array} $

(2)

引入Fomel^[19]的局部角频率ω_loc的计算策略

$ {\omega _{{\rm{loc}}}} = {\left[ {{\lambda ^2}\mathit{\boldsymbol{I}} + \mathit{\boldsymbol{S}}\left( {\mathit{\boldsymbol{V}} - {\lambda ^2}\mathit{\boldsymbol{I}}} \right)} \right]^{ - 1}}\mathit{\boldsymbol{SU}} $

(3)

式中：λ为平滑尺度因子；I为单位向量；S为整形正则化算子^[20]；U和V分别为u(t)和v(t)组成的向量矩阵。

为了验证局部频率计算结果的合理性，设计两个用于对比的瞬时频率与局部频率信号(图 1)。图 2为线性调谐信号与实际地震信号的瞬时频率。由图可见，由于瞬时频率是相位对时间的导数，计算结果出现很多“负频率”现象，利用负频率奇异点不能直接确定动态匹配追踪频率参数搜索范围，需要进行全局扫描，因此不是真正意义上的动态搜索方式。图 3为线性调谐信号与实际地震信号的局部频率。由图可见，线性调谐信号的局部频率计算结果与理论值一致(图 3a)，实际地震信号的局部频率计算结果较为合理(图 3b)，在理论范围内不存在“负频率”现象，可以用于确定最佳匹配原子的频率参数。

首先确定最大包络振幅的时间位置u₀，即为此次迭代搜索的中心，进而得到该时间位置u₀处的局部角频率ω₀和相位φ₀，以此创建匹配搜索的先验信息点α₀(u₀, ω₀, φ₀)，并设定相邻控制参数的间隔Δα，得到最佳匹配原子的搜索范围γ′=[α₀-nΔα, α₀+nΔα](n为正整数, α为时频原子的三个控制参数u、ω和φ, Δα为相邻控制参数α的间隔)，即得到迭代搜索的动态字典原子库D′={g_γ′(t)}(D′∈D, g_γ′(t)为时频原子), D为过完备字典原子库。

假定第1次迭代需要扫描动态字典原子库中的m个时频原子

$ \begin{array}{l} \mathop g\limits_{i = 1, \cdots , {m^{\gamma '}}} = \left\{ {{u_0}, {\omega _i} \in \mathit{\boldsymbol{U}}\left[ {\omega \left( {{u_0}} \right), {\delta _\omega }} \right]} \right., \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;{\varphi _i} \in \mathit{\boldsymbol{U}}\left[ {\varphi \left( {{u_0}} \right), {\delta _\varphi }} \right]\} \end{array} $

(4)

式中：u₀为每个最佳匹配子波的搜索中心位置处的时间；ω(u₀)和φ(u₀)分别为搜索中心位置处信号的瞬时角频率和相位信息；U[ω(u₀), δ_ω]、U[φ(u₀), δ_φ]分别为角频率和相位的搜索邻域, δ_ω、δ_φ为相应参数搜索半径。

第1次迭代与信号最佳相关的时频原子g_γ1(t)满足

$ \left| {\left\langle {f, {g_{{\gamma _1}}}} \right\rangle } \right| = \mathop {\sup }\limits_{i = 1, \cdots , m} \left| {\left\langle {f, {g_{{\gamma _i}}}} \right\rangle } \right| $

(5)

则信号f被表示为沿g_γ1(t)方向的分量和残差分量，即

$ f = \left\langle {f, {g_{{\gamma _1}}}} \right\rangle {g_{{\gamma _1}}} + {{\rm{R}}^1}f $

(6)

式中R¹f为第一次迭代后的残差信号。按上述方式依次迭代下去，直至达到预设的迭代次数和残差阈值，则迭代停止。若经过n次迭代后，完成对信号的分解过程，最终将信号表示为n个最佳匹配原子的线性组合

$ f = \sum\limits_{i = 1}^n {\left\langle {{{\rm{R}}^{i - 1}}f, {g_{{\gamma _i}}}} \right\rangle } {g_{{\gamma _i}}} + {{\rm{R}}^n}f $

(7)

式中：R^i-1f为i-1次迭代后的残差信号；Rⁿf为n次迭代后的残差信号，当i=1时, R⁰f为原始信号。

利用9个Morlet小波合成理论信号，基于局部频率约束的匹配追踪算法迭代分解理论信号，图 4为无噪声理论信号的迭代分解过程。由图可见，每一次迭代优选出一个最佳匹配的原子，经过9次迭代后，重构信号与理论信号高度吻合，残差曲线为零值。为验证本文提出的匹配追踪算法的抗噪性，对理论合成信号进行加噪，在信噪比为5:1的情况下，对加噪信号迭代分解，图 5为加噪理论信号的迭代分解过程。由图可见，基于局部频率约束的匹配追踪算法对信号的迭代分解过程稳定，经过9次迭代后，得到9个最佳匹配原子重构理论信号，残差曲线基本为噪声信号，证明该算法具有一定抗噪性。

1.2 基于局部频率约束的动态匹配追踪算法

在对原始数据进行强反射识别与分离时，可以利用层位约束或一定的时窗范围内信号包络振幅最大位置拾取强反射具体位置。通过求取强反射时间位置的复地震道的局部频率、瞬时相位信息，构建局部反射特征波形库，利用式(5)搜索最佳相关的子波作为强反射子波。在对实际数据进行强反射剥离前，需要先确定强反射剥离的强度^[21]，即

$ \frac{1}{M}\left\| {{d_N} - \xi \times {d_{{\rm{strong}}}}} \right\|_2^2 = \frac{1}{N}\left\| {{d_M}} \right\|_2^2 $

(8)

式中：d_N为强反射时窗内地震数据能量, d_M为背景时窗内地震数据能量, N和M分别代表其时窗内地震数据的采样点数；ξ为强反射压制参数；d_strong为拾取的强反射地震数据。

图 6为去除强反射前、后记录。根据压制强反射的能量均衡原则(式8)，认为使去除强反射后的平均能量(红色虚线框)与背景能量(蓝色虚线框)基本保持一致的ξ为最佳强反射压制系数，能够达到较为理想的强反射压制效果(图 6b)。图 7为M区不同ξ的单道处理结果。为了满足能量均衡原则(式8)，即使去除强反射后的平均能量(红色虚线框内)与地震数据背景能量(蓝色虚线框内)保持一致，通过对比选取ξ=0.6为最佳强反射压制系数，以此系数对实际地震数据进行强反射剥离，能够较为理想地压制强反射并突显目标储层弱反射信息(图 7d)。

2 模型试算

将本文提出的基于优化局部频率约束的动态匹配追踪强反射识别与分离方法用于模型数据的处理，以验证方法可靠性。构建含高阻层的理论地质模型(图 8)，利用28Hz的Morlet小波合成理论模型的正演地震记录(图 9)，可见由于受高速层的影响，顶部薄砂体的有效反射信息被强反射淹没，难以识别上覆砂体真正的形态和位置。图 10为由本文方法识别和剥离强反射后的地震记录。由图可见，明显地削弱了强干扰的影响，能够恢复上覆薄砂岩的地震响应。为了更清晰地看出强反射对砂体特征识别的影响，图 11和图 12分别为去除强反射前、后的地震道集。由图可见，剥离强反射后突显了砂体形态特征，有助于对储层进行精细描述，为后续地震数据处理提供了有效的数据基础(图 12)。

3 实际资料处理

A区位于济阳坳陷，该区存在多级不整合，利于发育与不整合有关的油藏。实际地震资料存在典型的不整合反射特征，为准确描述岩性圈闭提供了基础数据。由于该区储层倾角小、储层厚度小且目的层上、下存在高阻抗体，地震反射特征不能真实反映不整合圈闭的地质特征。为了准确地拾取砂体反射特征信息，有必要识别与分离强反射。

图 13为原始地震记录。由图可见，砂岩尖灭反射信息淹没于上覆高速层强反射中，难以识别真正的砂体尖灭点。考虑到强反射记录具有很好的横向连续性，利用本文提出方法将T2和Tr强反射从工区实际数据中剥离，从剥离后的地震记录(图 14)上可见，砂体的反射信息得到突显。从去除强反射前、后的地震道集上也能明显地看到：由于Tr强反射与砂体反射干涉、叠加，无法准确识别砂体尖灭位置(图 15)；剥离Tr和T2强反射后砂体尖灭反射信息更清楚(图 16)。因此，高速层强反射识别与剥离技术为砂体尖灭识别提供了可靠的数据。

4 结束语

为削弱目标储层上、下高阻层形成的强反射干扰影响，本文提出了基于局部频率约束的动态匹配追踪强反射识别与分离方法。考虑到瞬时频率计算结果存在“负频率”现象，认为其不适用于作为时频原子的频率搜索参数，本文提出局部频率计算模式，计算结果较为合理，将其作为构建动态子波库的控制参数。利用强反射位置处的局部频率和瞬时相位信息作为先验信息，构建动态反射特征波形库，以此搜索最佳匹配强反射子波，通过单道测试确定强反射压制系数，以保证剥离强反射后的均衡能量与背景能量一致，有助于拾取目标储层的反射信息，最终实现储层的精细描述。通过构建含高阻层的模型数据验证了基于局部频率约束的动态匹配追踪强反射识别和分离方法的可行性，并将该方法应用于实际地震数据，削弱了强反射干扰影响，突显了弱反射信息。