地震反演驱动的改进匹配追踪煤层识别方法

潘辉， 印兴耀，2*， 李坤，2， 裴松

1 中国石油大学(华东)地球科学与技术学院， 青岛 266580 2 海洋国家实验室海洋矿产资源评价与探测技术功能实验室， 青岛 266580

0 引言

我国含油气盆地的油气地震勘探中，煤系地层的存在给勘探开发带来了极大的阻力，煤系地层中煤层纵向上多层叠置，多套煤层与砂岩直接接触形成砂煤薄互层，含煤储层地震响应信息夹杂煤层反射信息，煤层反射干扰了储层的振幅和频率信息，增加了地震解释的不确定性(李海山等, 2014; Liu and Ghosh, 2017).一方面，通常在地震剖面上看到的亮点被错误地解释为油气存在的指示；另一方面，极低纵波阻抗引起的强烈振幅可能会掩盖煤层附近储层的反射(Guo and Wang, 2004; Liu et al., 2014).如何在煤层地震响应微弱的基础上建立有效的地震数据目标处理方法，为优质储层评价和流体识别奠定数据基础是一个亟待解决的问题.

目前，传统的煤层识别方法可分为两类：(1)基于测井资料的煤层识别方法(沈怀磊等,2010)；(2)基于地震资料的煤层识别方法(张玉忠和杨永波,2009;甘嫦华和喻岳钰,2010; 张珊珊等,2015).针对以上基于地震资料的煤层识别方法，都是通过地震资料选择对煤层识别最有效的地震属性或其组合，并借助相关处理和解释方法提高煤层解释的精度.笔者认为煤层反射识别要得到理想的结果，关键在于首先要针对工区特点选取对煤层敏感的地球物理参数.其次，由于地球物理反演的多解性问题，仅靠一个地球物理参数无法准确预测煤层.本文提出的煤层指示因子优点在于综合了地震、岩石物理与信号领域，达到直观、精确识别煤层的目的.

近年来，稀疏表示理论的研究日趋深入，在地震随机噪声压制、反褶积、地震信号时频分析、地震数据重建等地震信号处理领域得到了广泛的应用(田琳和胡津健，2021).针对地震信号的稀疏表示方法，国内外学者对匹配追踪-MP(L0范数)算法作了大量应用与研究(Mallat and Zhang, 1993; Liu et al., 2004; Liu and Marfurt, 2005)，其主要体现在非平稳地震信号的时频分析(Castagna et al., 2003; Partyka et al., 1999)、稀疏地震反演(Nguyen and Castagna, 2010; 印兴耀等, 2020; 李坤等, 2018)、地震稀疏分离去噪处理(Hou et al., 2018; Li and Sacchi, 2021)、瞬时谱提取(Wang, 2007, 2010; 张繁昌等, 2012; 李坤等, 2016)、高分辨率反射率反演(Rodriguez et al., 2010)、尖灭点识别及滤波(赵天姿等, 2008; 张繁昌等, 2012)、煤层强反射识别与分离(李海山等, 2014; 许璐等, 2019).匹配追踪(Matching Pursuit, MP)算法能在一定程度上提高地震分辨率，但面临的主要问题是算法最优时效果可能不是最优，因此有必要明确地质目标，以便用匹配追踪算法重建的地震数据能够对预期的地质体是可分辨的，故本文提出地震反演驱动的改进匹配追踪煤层识别方法.

基于以上研究背景，针对地震资料中储层地震响应信息夹杂煤层反射信息，且煤层地震响应特征不够明显的情况，常规的匹配追踪煤层反射分离在煤层单一、厚度变化较小的情况下是比较精确的；然而当地层广泛发育多套薄煤层，波形不规则，且横向连续性较差时，煤层反射特征不够明显，很难准确识别煤层波形.故本文采用基于煤层指示因子信息约束的改进匹配追踪方法对煤层反射进行识别与分离，首先基于岩石物理分析构建煤层敏感指示因子，进而利用匹配追踪算法反演煤层敏感因子，精确拾取主要煤层所在位置为后续匹配追踪算法提供先验信息.其次考虑到复数道动态匹配追踪算法计算效率偏低、所求取的瞬时频率结果存在无物理意义的“负频率”且易受噪声影响(Boashash,1992)，即当计算结果出现异常值时，只能全局搜索时频原子频率属性，这并不是真正意义上的动态搜索.为此，本文引入基于EMD字典的改进匹配追踪算法，以缓解现有匹配追踪算法计算效率低的问题(潘辉等，2021)，利用煤层反射先验信息约束改进匹配追踪搜索范围，进而将煤层反射记录与原始地震记录分离，减少振幅“假亮点”及其对下伏地层造成的砂体连续分布的假象，通过含煤储层地质模型和实际数据处理，验证了该方法的可行性.

1 煤层指示因子地震匹配追踪反演方法

1.1 基于岩石物理分析的煤层指示因子构建方法

岩石物理分析提供了煤层性质和地震属性之间的联系，在进行大规模煤层数据识别之前，研究敏感的煤层岩石物理参数，对于模拟每一个薄煤层的地震响应是至关重要的，使用我国某油田的数据来对煤层敏感岩石物理参数进行识别和排序，为后续煤层匹配追踪和识别奠定了坚实的基础.

首先，基于理论地震岩石物理模型估算岩石弹性参数；然后，分析岩石弹性参数与地层煤质含量间的量化关系；其次，利用测井数据进行煤层指示因子的交会图分析，计算不同弹性参数对煤层的指示系数，构建出最敏感的煤层指示因子为：

(1)

式中，γ表示基于统计学理论计算得到的弹性参数的煤层指示系数，L表示煤层发育段的实际测井数据，M表示非煤层发育位置的实际测井数据，mean(·)表示均值运算符,只对深度方向上进行平均，std(·)表示标准差运算符.

如图1纵波阻抗和纵波速度的交会图所示，其中煤层为黑色对应于低纵波阻抗、砂岩为黄色对应于高纵波阻抗、红色泥岩的纵波阻抗相对高于煤层且低于砂岩.根据计算的敏感指示系数，结合岩石物理交会，观察到纵波阻抗对煤层较为敏感,在纵波阻抗低于一定阈值后可以清楚地区分煤层和砂泥岩.

在岩石物理交会图分析的基础上，根据公式(1)计算不同弹性参数对煤层的指示系数，寻找出最敏感的煤层指示因子.

如图2所示，纵坐标为根据公式(1)算出的γ，表示基于统计学理论计算得到的弹性参数的煤层指示系数，横坐标为各种煤层地震弹性参数，由图1a、b清晰看出相较于其他地震弹性物理参数，纵波阻抗的煤层指示系数远远高于其他值，这是结合实际工区计算结果，不同工区具有不同敏感参数.

图2 单井煤层敏感性分析示意图(a) A井煤层敏感性分析; (b) B井煤层敏感性分析; (c) 横坐标示意图.Fig.2 Coal seam sensitivity analysis diagram of single well(a) Sensitivity analysis ofcoal seam in well A； (b) Sensitivity analysis of coal seam in well B； (c) Horizontal coordinate diagram.

为进一步验证纵波阻抗在本工区的适用性，对实际工区多井进行煤层敏感性分析见图3，其中纵坐标为煤层指示系数γ，横坐标为各种煤层地震弹性参数，不同颜色表示实际不同测井数据.

从图3可以看出在实际工区多井数据分析中对比不同岩石弹性参数纵波阻抗煤层敏感指示系数γ最高.在地震学上，煤层是低速、低密的夹层，基于纵波阻抗包含速度、密度等岩石物性参数的信息，故结合岩石物理分析，提出将某一阈值下纵波阻抗作为煤层指示因子用于煤层预测.

图3 多井煤层敏感性分析示意图Fig.3 Coal seam sensitivity analysis diagram of multi-well

1.2 基于匹配追踪的煤层指示因子反演方法

在实际地震资料处理中测井数据较少，对于大工区的地震数据，基于测井数据获得的纵波阻抗远远不够，往往需要借助地震反演获得纵波阻抗.由于地震数据是带限的且受到随机噪声的干扰，导致地震反问题存在较强的不适定性(杨培杰和印兴耀,2008；李坤等,2016；肖爽等，2020).稀疏正则化方法能够较好地缓解地震反问题的病态性，以获取模型参数的近似解(Tarantola,2005；Tikhonov,1963).本文利用L0范数正则化方法，假设地层反射系数具有稀疏性，构建了基于匹配追踪的反演目标泛函，且加入了模型参数的低频先验信息约束项，以增强反演稳定性(Alemie and Sacchi,2011；Zong et al.,2015; Yin and Zhang,2014; Yin et al.,2016；杨俊等，2020).匹配追踪稀疏反演是在传统反演方法的基础上，对模型参数进行稀疏约束度量(Chopra et al., 2006;Wen et al.,2015;刘晓晶等,2015;李坤等,2018).通过控制反射系数个数、迭代次数和迭代阈值，经过有限步迭代搜索获得目标最优解.将模型约束的地震反演目标函数与匹配追踪优化算法相结合，得到的稀疏目标泛函为：

(2)

通过求解方程(2)获得地层相对反射系数r，进而地层的绝对波阻抗IP(t)可由下面积分方程得到：

(3)

1.3 基于煤层指示因子预测煤层空间展布

基于上文岩石物理的研究成果，提出以煤层硬阈值约束下的地层波阻抗作为煤层指示因子Ccoal，预测煤层空间展布，即：

Ccoal=IP(t),{IP(t)

(4)

其中e为针对研究区煤层岩性特点，基于岩石物理分析得出的煤层阈值范围.

此时煤层指示因子Ccoal只有两种情况非零值与零值，为了便于计算，故令：

(5)

进而确定煤层指示因子等于1的时间位置为t(Clocation)i(i=1,2,3,…,n)，n表示每道地震记录确定的位置个数.为通过煤层指示因子时间位置对煤层顶底板进行区分，对煤层指示因子时间位置求导得t_ci=t′(Clocation)i(i=1,2,3,…,n)，则可以对煤层顶底板进行区分，假定煤层顶板为t(Clocation)roof，煤层底板为t(Clocation)floor，则本文将煤层位置分布范围确立为δt=[t(Clocation)-nΔt,t(Clocation)+nΔt]，其中，煤层地震响应长度为n=t(Clocation)roof-t(Clocation)floor，Δt为地震记录的采样间隔.因而，在对煤层信号进行匹配追踪识别时，利用煤层时间位置t(Clocation)i(i=1,2,3,…,n)作为动态最优搜索的先验信息，实现煤层信号的快速高效识别.

为更清晰地阐述煤层先验位置计算过程，具体流程如下：

(1)利用快速匹配追踪反演纵波阻抗IP，令煤层指示因子Ccoal=inverse_IP;

(2)基于岩石物理分析，设置煤层硬阈值e，寻找符合条件的煤层指示因子位置; Location1 =find(Ccoal

(3)计算煤层指示因子位置导数，便于分辨煤层顶底板，t_c=diff(Location1)；

(4)寻找煤层底部边界位置，Location2= find(tc～=1);

(5)计算煤层数量响应，Coal_number=length(Location2);

(6)寻找煤层顶部边界位置，Location3=location1(Location2(coal_number-1)+1);

(7)计算煤层响应时间长度，n= Location2-Location3+1;

(8)故最终约束煤层匹配追踪搜索范围为sigamt=[Location3-n×delat Location2+n×delat]、delat=0.002 s.

图4是单道煤层指示因子示意图，黑色为煤层指示先验信息，红色为所设计的与之对应高斯窗处理后的先验信息，在实际资料处理中，为保证实际资料的纵向连续性，对提取出的煤层指示先验信息进行了加窗处理.从图中能清晰的分辨出每个煤层响应的顶底板时间位置t_ci(i=1,2,3,…,7)，从而获得每个煤层响应的地震时间厚度，进而约束匹配追踪搜索邻域.

图4 单道地震信号煤层指示因子示意图Fig.4 Coal seam indicator diagram of single seismic signal

图5是叠后地震数据，图6是与之相对应的煤层指示因子剖面，由图可见，通过与实际测井数据进行对比，煤层指示因子可以明显区分煤层所在位置，且指示的煤层横向连续性较高、边界刻画清晰，无井区域也具有一定的煤层指示效果.故本文提出将煤层硬阈值约束下纵波阻抗作为煤层指示因子具有一定可行性.

图5 叠后地震剖面Fig.5 Post-stack seismic section

图6 煤层指示因子剖面Fig.6 Indication factor section of coal seam

图7是根据煤层指示因子刻画的煤层三维空间展布示意图，由图可见，通过与实际测井数据进行对比，煤层指示因子可以明显区分煤层所在位置，且指示的煤层横向连续性较高、边界刻画清晰，无井区域也具有一定的煤层指示效果，表明该参数能够较好的刻画三维空间展布.

图7 煤层指示因子三维空间展布示意图Fig.7 Three-dimensional spatial distribution of coal seam indicator factors

2 基于EMD字典的改进匹配追踪算法

2.1 基于经验模态分解构建EMD字典

本文通过地震反演获得煤层位置的先验信息，将其引入匹配追踪算法中约束煤层原子搜索范围，可缓解现有匹配追踪算法运行效率低的问题；同时，本文将经验模态分解(EMD)引入匹配追踪算法，给出了EMD字典构建方法(潘辉等，2021).具体地，将EMD视为一种基于过完备时频字典的稀疏分解方法(Huang et al., 1998; Liu and Chen, 2019; Huang et al., 2020；邬蒙蒙等，2020；贺月等，2021)，在计算匹配煤层原子主频的过程中，引进基于连续相位的阻尼最小二乘反演求解瞬时频率(Taner et al., 1979; Shatilo, 2010; McGowan and Kuc, 1982)，并加入整形正则化算子(Fomel, 2007)对数据进行平滑处理，有效避免了连续相位计算过程中奇异值的出现，改善了匹配煤层原子瞬时频率的求取方法，得到了稳定的瞬时频率(Vincent and Bengio, 2002; Durka et al., 2005; Gribonval, 2011)，同时显著提升了计算效率.

本文给出的EMD字典是根据EMD算法分解获得的所有可能的固有模式函数(Instinsic Mode Function, IMF)的集合来定义(Hou and Shi, 2013)：

Dimf(t)=a(t)cosφ(t):a(t),φ′(t)∈V(φ,λ),φ′(t)

≥0，

(6)

其中Dimf(t)为EMD分解后的一个内蕴模式函数，a(t)为瞬时振幅函数，瞬时相位函数φ(t)为递增函数，其导数φ′(t)可以表示瞬时频率.由于IMF函数的前提假设是平稳的，故信号的能量和瞬时频率不随时间发生剧烈变化，因此a(t)和φ′(t)需要满足属于由瞬时相位φ(t)定义的谐波列向量所组成的子空间V(φ,λ)，一般来说，将V(φ,λ)构造为过完备的傅里叶基是最有效的.V(φ,λ)为：

(7)

D={ωγ(t)}γ ∈Γ={ωγ=(tc,fc,φc)(t)}γ ∈Γ,ωγ(t)

=a(t)cosφ(t)，

(8)

其中，D表示过完备字典原子库，ωγ(t)表示过完备字典原子库的母小波，ωγ=(tc,fc,φc)(t)为调制后的原子，Γ表示上述时间范围及其对应的频率和相位的矩阵，a(t)表示振幅，φ(t)表示瞬时相位，γ=(tc,fc,φc)为控制煤层信号的参数集合，tc表示为煤层指示因子约束下最优煤层时频原子的中心时间，fc表示煤层指示因子约束下最优煤层时频原子的瞬时频率，φc表示煤层指示因子约束下最优煤层时频原子的瞬时相位.

本文搜索方法地每次迭代可以简单表示为：

γi={tc∈U[δt],fc∈U[f(tc),δf],φc∈U[φ(tc),

δφ]} (i=1,2,3,…,L)，

(9)

式中，L表示每一次迭代搜寻所确定的最优匹配子波的个数；δt为上文约束的煤层时频原子匹配范围；其中，U[f(tc),δf]和U[φ(tc),δφ]分别在煤层指示因子约束下频率与相位的搜索邻域集合，f(tc)和φ(tc)为煤层时频原子中心时间的瞬时相位与瞬时频率，δf和δφ为煤层指示因子约束下相应参数搜索半径.与全局搜索相比，这种动态搜索极大提高了计算效率；且本文每次迭代确定的时间位置是在煤层指因子约束下将信号能量所有极大值处作为每次迭代的初始时间位置，这样既保证了匹配追踪的精度，同时也提高了计算效率；为本文匹配追踪算法应用于三维煤层数据处理提供了可能.

2.2 基于连续相位的改进瞬时频率求解方法

针对传统瞬时频率的负值现象及计算过程鲁棒性差的问题，本文引入了连续相位，并借助阻尼最小二乘法，加入整形正则化算子(Fomel, 2007)对数据进行平滑处理，改善了匹配原子瞬时频率的求取方法，得到了稳定的瞬时频率，并从中获得匹配原子主频.

连续信号x(t)的复地震道为：

S(t)=x(t)+jh(t)=A(t)ejφ(t)，

(10)

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

传统的瞬时频率f(t)是瞬时相位φ(t)的变化率，即：

(11)

将连续相位替代瞬时相位φ(t)，引入阻尼最小二乘法，并加入整形正则化算子求取瞬时频率：

f(t)=[p2I+(ZO)T(ZO)]-1(ZO)Tl，

(12)

式中：p为权系数，一般取O中元素最大值的1%～5%即可；I为单位向量；Z为整形正则化算子；l和O分别为l(t)和o(t)组成的向量矩阵.

为了验证该瞬时频率计算结果的合理性，设计了1个理论地震信号用于对比瞬时频率的求解差异(图8).图8a—f分别为合成地震信号、传统瞬时频率、基于连续相位求导的瞬时频率、连续相位常规阻尼最小二乘法所求瞬时频率、局部频率以及基于连续相位改进阻尼最小二乘法所求瞬时频率，由图8b—c对比可看出，由于传统瞬时频率是相位对时间的导数，计算结果存在较多的“负频率”.而基于连续相位求导的瞬时频率，虽然避免了出现负值的现象，但在0.3～0.45 s等处瞬时频率出现异常.图8d、f为基于连续相位直接求逆与改进阻尼所求瞬时频率，可见阻尼最小二乘法解决了瞬时频率异常现象，并且改进阻尼所求瞬时频率比常规阻尼更加平滑，更加突显频率异常点；图8e为局部频率，由图8e—f对比可看出，基于连续相位求解的瞬时频率能更加符合真实频率值，且更加突显高频成分.改进阻尼最小二乘法所求瞬时频率，相比局部频率能更加符合真实频率值，且更加突显高频成分.

图8 地震信号频率求解策略(a) 合成地震信号； (b) 瞬时频率； (c) 连续相位求导的瞬时频率； (d) 连续相位直接求逆所求瞬时频率； (e) 局部频率； (f) 连续相位改进阻尼所求瞬时频率.Fig.8 Frequency solving strategy of seismic signal(a) Synthesis seismic signal； (b) Instantaneous frequency； (c) Instantaneous frequency obtained by derivative of continuous phase； (d) Instantaneous frequency obtained by inversion of continuous phase； (e) Local frequency； (f) Instantaneous frequency of continuous phase modified damping.

图9为理论信号含信噪比(1∶1)噪声的情况下进行抗噪性测试，从图9的结果，可以看出，基于本文的瞬时频率求法能够很好地处理加噪地震数据，并且相对于局部频率曲线更加光滑，更加突显频率异常点，证明本文方法有很好的抗噪性.

图9 加噪地震信号(信噪比1∶1)频率求解策略(a) 加噪地震信号； (b) 瞬时频率； (c) 连续相位求导的瞬时频率； (d) 连续相位直接求逆所求瞬时频率； (e) 局部频率； (f) 连续相位改进阻尼所求瞬时频率.Fig.9 Frequency solving strategy of noisy seismic signal(a) Noisy seismic signal； (b) Instantaneous frequency； (c) Instantaneous frequency obtained by derivative of continuous phase； (d) Instantaneous frequency obtained by inversion of continuous phase； (e) Local frequency； (f) Instantaneous frequency of continuous phase modified damping.

本文提出的地震反演驱动的改进匹配追踪数学原理为：

(13)

其中，S为含煤地震信号，ωi为匹配分解的EMD原子，n为实际地震资料中所包含的随机噪声，ai为煤层匹配原子的相应振幅值，D为地震反演信息约束下构建的超完备EMD字典.本文匹配算法采用了批量筛选煤层原子的策略，首先需要为每个匹配过程设置硬迭代阈值δ，一般0≤δ≤1，同时满足煤层指示因子约束与迭代阈值的原子将被列为候选匹配集合，具体表述为：

|Si|≥δ|kmax|,Si=〈ωi,Rs〉,

(14)

式中，Si为前一次迭代的残差与第i个EMD原子的内积值；Rs为迭代残差向量，kmax为每次迭代内积向量中的最大绝对值，假设地震信号经过N次迭代后，基于匹配的M(M>N)个原子构成原子矩阵JN=[ω1,ω2,…ωM]，然后利用阻尼最小二乘算法对匹配的煤层原子振幅进行逐次修正，即：

aN=[(JN)T(JN)+σ2I]-1(JN)TS，

(15)

式中，aN为第N次迭代后的修正振幅，IM×M是单位矩阵，σ2为阻尼因子，S为残差信号.假设SN是第N次迭代的残存信号，搜索得到的最优匹配原子是ωn+1，则SN可表示为：

SN=S-JNaN，

(16)

其中ωn+1满足：

|〈SN,ωn+1〉|=supi∈(1,2,…N)|〈SN,ωi〉|，

(17)

其中〈〉表示内积，supi∈(1,2,…N)表示上确界，即选择ωn+1∈D，使得ωn+1与SN最为匹配.

(18)

3 地震反演驱动的改进匹配追踪煤层识别

本文通过地震反演的煤层指示因子精确拾取主要煤层所在位置为改进匹配追踪算法提供先验信息，约束基于EMD字典匹配追踪算法的搜索范围，进而将煤层反射信息在EMD字典中匹配出来，对煤层反射进行定量表征与识别，可有效消煤层反射对目的层有效反射消息的屏蔽作用,并大大提高了计算效率.煤层指示参数约束下改进的匹配追踪煤层反射识别示意图(图10)所示.图10a中黑色为煤层指数先验信息，红色为所设计的与之对应高斯窗；图10b中黑色为煤层指示因子指示的煤层地震数据，绿色为与之对应的高斯窗处理后数据，红色为利用改进匹配追踪算法匹配出的煤层地震数据；图10c中蓝色为实际地震记录，红色为与之对应的分离煤层后地震记录.基于煤层指数先验信息提供的数据连续性较差，故设计了与之对应的高斯窗对其进行平滑处理，从本文方法提取的煤层地震记录，以及原始煤层记录与分离煤层后地震记录对比可以发现，在煤层先验位置处对原始地震数据有较好的分离效果，煤层反射信息得到了很好的压制，从图10c中可以看出去除煤层处理后地震记录仍与周边记录保持较好的横向连续性，且去除后地震记录(红线)与去除前地震记录(蓝线)相位正负未发生变化，波形相似，并未发生波形极性反转现象.

图10 基于煤层指示约束的匹配追踪反射识别示意图Fig.10 Schematic of matching pursuit reflection recognition based on coal seam indication constraint

详细基于煤层指示因子约束的匹配追踪反射识别流程(图11)如下：

(1)分析煤系储层中煤层敏感岩石物理参数, 构建煤层敏感指示因子；

(2)提出基于快速匹配追踪的煤层指示因子地震反演算法；

(3)引入经验模态分解思想，构建基于EMD字典的改进匹配追踪；

(4)研发地震反演驱动的改进匹配追踪煤层识别方法.

图11 基于煤层指示因子约束的匹配追踪反射识别流程Fig.11 Matching pursuit reflection recognition process based on coal seam indicator factor constraint

为了验证上述匹配追踪识别煤层方法的可行性，我们首先建立不同砂泥煤耦合方式进行了模型试算，设计了如图12、图13、图14不同地质背景下的含煤层理论模型，其中砂岩速度为4160 m·s-1、煤层速度为2850 m·s-1、泥岩速度为3885 m·s-1.图(a)为理论纵波阻抗模型，图(b)为理论反射系数序列，其中因煤层产生的反射系数出现在160ms附近；图(c)为正演的25 Hz零相位雷克子波，图(d)正演合成记录；由图(e)为利用本文算法匹配出的煤层反射记录；图(f)为煤层反射分离后地震道集；图(g)为处理前后单道记录对比，值得注意的是，图(g)中红线为煤层反射分离后的地震记录，黑线为对应含煤地震记录.图(d)中的合成记录与图(e)中匹配追踪记录可见两种地震记录波形基本吻合，图(f)、图(g)中处理道集前后对比可知，在不同地质背景下本文基于EMD字典的匹配追踪技术对煤层反射具有很好的分离效果，较好地恢复了弱小反射层的反射信息.

图12 单套煤层反射模型及匹配追踪识别结果(上下砂岩且速度一致)(a) 含煤层波阻抗模型； (b) 反射系数序列； (c) 雷克子波； (d) 合成地震记录；(e) 匹配出的煤层反射； (f) 分离煤层反射后记录； (g) 分离煤层反射前后单道记录对比.Fig.12 Single set of coal seam reflection model and matching pursuit recognition results (upper and lower sandstone and consistent velocity)(a) Impedance model containing coal seam； (b) Reflection coefficient sequence； (c) Rayke wavelet； (d) Synthetic seismogram； (e) Matching coal seam reflection； (f) Records after separation of coal seam reflection； (g) Comparison of single record before and after reflection in separated coal seam.

图13 单套煤层反射模型及匹配追踪识别结果(上覆砂岩，下覆泥岩)(a) 含煤层波阻抗模型； (b) 反射系数序列； (c) 雷克子波； (d) 合成地震记录； (e) 匹配出的煤层反射； (f) 分离煤层反射后记录； (g) 分离煤层反射前后单道记录对比.Fig.13 Single set of coal seam reflection model and matching tracking recognition results (overlying sandstone, underlying mudstone)(a) Impedance model containing coal seam； (b) Reflection coefficient sequence； (c) Rayke wavelet； (d) Synthetic seismogram； (e) Matching coal seam reflection； (f) Records after separation of coal seam reflection； (g) Comparison of single record before and after reflection in separated coal seam.

图12为上下砂岩一致背景下含煤层反射界面理论模型，可看出匹配追踪技术对煤层反射有很好的分离效果，匹配追踪去煤层后的地震记录与不含煤层模型得到的地震记录基本吻合.

图13为上覆砂岩、下覆泥岩背景下含煤层反射界面理论模型，可看出在砂煤耦合时，匹配追踪技术对煤层反射有很好的分离效果，去煤层后反射地震记录与无煤层背景下地震记录基本一致.

图14为4套煤层复合波模型及匹配追踪识别结果，从图14e、f可见，煤层复合波引起的反射匹配识别效果依然保持较好，煤层反射分离后，弱小反射层的反射信息依旧得到了较好的恢复.

4 模型试算

为了进一步验证煤层指示因子约束下基于EMD字典的匹配追踪与识别方法的可行性和实用性.依据工区实际地质情况构建煤层反射背景下砂泥煤岩薄互层理论地质模型(图15)，具体参数如表1所示，模型大小为横向上1150CDP，纵向上采样时长1118 ms，图16为通过快速匹配追踪算法反演的纵波阻抗结果，图17为理论模型的低频模型，图18为理论模型的煤层指示因子剖面图，图19是通过25 Hz的 Ricker子波合成理论模型的正演模拟记录，从图16、图17、图18与图19可看出由于受到广泛发育的薄煤层影响，薄砂体的反射信息被煤层反射所吞没，难以分辨出砂体的真实形态和位置，而煤层指示因子准确的刻画了煤层位置分布.图20为单道地震记录基于煤层指示因子约束的匹配追踪反射识别示意图，展示了基于煤层指示因子约束下识别的煤层地震记录与去除后的煤层地震记录.图21为由本文方法剥离煤层反射后的地震记录，图22为由本文方法识别的煤层反射地震记录，图23为理论模型正演的煤层地震记录，图21和图22的残差信号如图24所示，对比图19、图21、图22、图23和图24，清晰看出本文方法识别的煤层反射记录与理论模型正演记录基本一致，且明显地减弱了煤层干扰的影响，较好的恢复薄砂岩的地震响应，且两者残差为零.图15—图24验证了对于复杂砂泥煤岩薄互层纵波阻抗二维模型，基于EMD字典的匹配追踪具有一定的煤层识别效果，对比分离煤层记录与理论煤层合成记录，表明对于广泛分布的复杂砂泥煤模型，本文方法可以达到基本识别每一个煤层，并且具有较好的分离效果,表明对于复杂砂泥煤耦合模型匹配追踪具有较好实用性.对比分离煤层前后记录可以发现，匹配追踪反射分离在不改变围岩反射信息的前提下有效恢复了弱小反射信息，并且匹配识别的煤层记录与理论煤层合成记录一致，分离后煤层地震记录与无煤层地震记录基本吻合.

表1 理论地质模型参数Table 1 Theoretical geological model parameters

图15 理论地质模型模型背景(绿色)为大套泥岩(平均速度为3524 m·s-1，密度为2.643 kg·m-3)，中部(黑色)为广泛发育的薄煤层(平均速度为3195 m·s-1，密度为2.274 kg·m-3)，薄层砂岩(红色)(平均速度为4162 m·s-1，密度为2.399 kg·m-3).Fig.15 Theoretical geological modelThe background of the model (green) is a large set of mudstone (average velocity 3524 m·s-1, density 2.643 kg·m-3), the middle (black) is a widely developed thin coal seam (average velocity 3195 m·s-1, density 2.274 kg·m-3), and the thin sandstone (red) (average velocity 4162 m·s-1, density 2.399 kg·m-3).

图16 理论模型的纵波阻抗反演结果Fig.16 P-wave impedance inversion results of the theoretical model

图17 理论模型的低频模型Fig.17 Low frequency model of the theoretical model

图18 理论模型的煤层指示因子剖面Fig.18 Coal seam indicator factor profile of theoretical model

图19 理论模型的正演地震记录Fig.19 Forward seismogram of the theoretical model

图20 单道地震记录煤层信号匹配追踪反射识别示意图Fig.20 Single-channel seismic record coal seam signal matching tracking reflection identification diagram

图21 由本文方法剥离煤层反射后的地震记录Fig.21 Seismic records after coal seam stripping by the proposed method

图22 由本文方法识别的煤层反射地震记录Fig.22 Coal seam seismic records identified by this method

图23 理论模型正演的煤层地震记录Fig.23 Coal seam seismic records forwarded by theoretical model

图24 识别的煤层记录与理论模型合成的煤层记录残差Fig.24 Coal seam records identified and coal seam records residual synthesized by theoretical model

为验证该算法的抗噪能力，对加噪地震记录(信噪比1∶1)进行测试，图25为理论加噪模型(信噪比1∶1)的纵波阻抗反演结果，图26为理论加噪模型(信噪比1∶1)的低频模型，图27为理论加噪模型(信噪比1∶1)的煤层指示因子剖面，图28是通过25 Hz的Ricker子波合成理论加噪模型(信噪比1∶1)的正演地震记录；从图25、图26、图27与图28对比表明，在含噪声的情况下砂体依然受到广泛发育的薄煤层影响，煤层反射屏蔽了薄砂层的反射信息，难以分辨出砂体的真实形态和位置，而煤层指示因子依然较为准确的刻画了煤层位置分布.图29是由本文方法剥离煤层反射后的地震记录，图30是由本文方法识别的煤层反射地震记录，图31理论加噪模型(信噪比1∶1)正演的煤层地震记录，图32 是识别的煤层记录与理论加噪模型(信噪比1∶1)合成的煤层记录残差，对比图28、图29、图30、图31和图32，表明本文方法在加噪情况下，仍有不错的适用性，识别的煤层反射记录与理论加噪模型(信噪比1∶1)正演记录基本一致，同时本文方法识别的煤层记录随机噪声有所降低，且依然较为明显地减弱了煤层干扰的影响，较好的恢复薄砂岩的地震响应，两者残差基本为随机噪声.

图25 理论加噪模型(信噪比1∶1)的纵波阻抗反演结果Fig.25 P-wave impedance inversion results of the theoretical noise model (SNR 1∶1)

图26 理论加噪模型(信噪比1∶1)的低频模型Fig.26 Low frequency model of theoretical noise model (SNR 1∶1)

图27 理论加噪模型(信噪比1∶1)的煤层指示因子剖面Fig.27 Coal seam indicator factor profile of theoretical noise model (SNR 1∶1)

图28 理论加噪模型(信噪比1∶1)的正演地震记录Fig.28 Forward seismic records of the theoretical noise model (SNR 1∶1)

图29 由本文方法剥离煤层反射后的地震记录Fig.29 Seismic records after coal seam stripping by the proposed method

图30 由本文方法识别的煤层反射地震记录Fig.30 Coal seam seismic records identified by this method

图31 理论加噪模型(信噪比1∶1)正演的煤层地震记录Fig.31 Coal seam seismic records of the theoretical noise model (SNR 1∶1)

图32 识别的煤层记录与理论加噪模型(信噪比1∶1)合成的煤层记录残差Fig.32 Identification of coal seam records and residuals of coal seam records synthesized by theoretical noise model (SNR 1∶1)

5 实际资料处理

研究区存在广泛发育的薄煤层，且砂煤泥耦合方式复杂，地震数据受薄煤层反射的干扰严重，容易引起岩性识别的虚假亮点，尤其在泥煤密集反射段；基于常规地震反演的储层敏感参数预测误差较大，且易受煤层反射的影响，难以落实储层的空间展布，为了精确提取砂体的反射信息，有必要对薄煤层的反射进行识别与分离.

研究区煤层套数约18～25套，单套煤层厚度范围为0.5～2 m，砂岩约16～17套，单套砂岩厚度范围为1～30 m，多套薄煤层与薄砂岩耦合，由于煤层与砂体反射信息相互干扰、叠加，无法精确分辨砂体位置.图33、图34、图35分别显示了过井分离煤层反射前后的地震记录以及识别的煤层地震记录，红色线框处为4个主要煤层所在位置.由图33可见，该工区存在广泛发育的薄煤层且煤层横向连续性较差，砂岩反射信息淹没于煤层反射中，难以识别真正的砂体，利用本文方法将复杂薄煤层反射从实际地震数据中分离，从分离后的地震记录(图34)可看到红色线框处的主要煤层得到了有效分离，移除后的煤层振幅相对较低，分离煤层反射的地震记录具有很好的横向连续性，地震能量得到了平衡，并且原先淹没在煤层反射中的有效砂体反射能够到突显.识别得到的煤层反射记录如图35所示，通过与实际测井数据进行对比可见，煤层反射信息得到了有效的匹配识别.分离煤层后的地震能量得到了平衡，被淹没的砂体反射信息得到了突出.明对发育多套薄煤层的地震记录，本文方法剥离煤层反射具有较好的效果，这样可以利用无煤层反射的地震剖面预测储层和识别流体，为获得目的层信息提供了有力的保证.

图33 原始地震记录Fig.33 Original seismic record

图34 本文方法煤层反射分离后地震记录Fig.34 Seismic records after coal seam reflection separation by this method

图35 本文方法识别的煤层地震记录Fig.35 Coal seam seismic records identified by this method

由图34可见，煤层反射分离后，微弱反射信息(绿色线框处)相对明显，但仅从分离前后地震剖面对比缺乏说服力，为了定量描述煤层反射分离后砂体在属性上的表征特征，本文对去除煤层前后单道地震记录匹配追踪时频谱分解，来进一步分析煤层反射层去除前、后的效果，如图36所示，煤层反射分离后，煤层的反射信息得到压制，砂体的弱小信息得到凸显(红色线框处)，进一步验证了本文方法的可行性.

图36 煤层去除前后匹配追踪时频分析Fig.36 Time-frequency analysis of matching pursuit before and after coal seam removal

为测试计算效率，对本文方法处理剖面时间进行了计算(该剖面344道，每道1001个采样点，采样间隔2 ms)，同时，在相同的阈值约束下本文利用阻抗提供先验位置信息约束常规复数道动态匹配追踪算法对煤层反射进行识别与分离，历时513.64 s，而本文算法处理时间为13.43 s，表明本文方法运行时间相较于实际大幅度减小，适用于三维地震数据处理，提高了算法的实用性.

为进一步说明本文方法优势所在，在此本文利用煤层指示因子作为先验信息对常规复数道动态匹配追踪识别煤层算法进行了改进，并对图33所示的实际资料进行了处理.图37和图38分别显示了过井分离煤层反射后的地震记录以及识别的煤层地震记录，对比分离前后的地震记录，常规方法取得了一定的效果，但在识别的煤层波形上可以清晰看到本文方法的效果更优，匹配出煤层的波形更为完整.

图37 常规方法煤层反射分离后地震记录Fig.37 Seismic records after coal seam reflection separation by conventional methods

图38 常规方法识别的煤层地震记录Fig.38 Coal seam seismic records identified by conventional methods

6 结论

本文针对常规匹配追踪贪婪算法很难有效识别复杂薄煤层反射信息，且效率较低，无法应用于大规模野外数据处理，应用前景受到限制.研究煤层相关敏感参数的基础上，综合测井与地震资料，发展了基于匹配追踪算法的煤层指示因子反演方法，提出了一种地震反演驱动的改进匹配追踪煤层反射识别分离的方法，准确提取煤层反射信息的同时有效提高匹配追踪计算效率.将该方法应用于实际资料处理中，能分辨出由煤层反射引起的靶区储层砂体连续分布的假象，有效突出煤层反射面附近的储层信息，这证明了本文方法的可行性和抗噪性，且证明了其高效性.这使得将MP应用于3D地震数据成为可能，煤层分离后的数据使得解释被煤层振幅覆盖的储层成为可能，为后续的工作提供了强有力的保证.总结该方法的适用性和应用局限性，主要存在以下两点需要注意：

(1)本文研究方法仍有较大的发展前景，适用于煤系地层以外，也适用于地层不整合强屏蔽、砂体顶部因泥岩脱水形成的高速层干扰、油页岩覆盖强屏蔽、灰岩覆盖强屏蔽等情形，须提出合适的指示因子作为先验信息约束.

(2)对于复杂地下沉积模式，分离参数设置是有所不同的，当非有效反射信息较强时，应将分离参数增强同时减小匹配次数，当非有效反射信息较弱时，应减小分离参数并同时增强匹配次数来提高分离精度.