基于匹配追踪时频谱计算的砂体尖灭线检测方法

汪瑞良， 张文珠， 刘徐敏， 董国辉， 李志晔

(中海石油(中国)有限公司 深圳分公司，深圳 518000)

基于匹配追踪时频谱计算的砂体尖灭线检测方法

汪瑞良， 张文珠， 刘徐敏， 董国辉， 李志晔

(中海石油(中国)有限公司 深圳分公司，深圳 518000)

匹配追踪算法能够实现信号的自适应分解。首先研究了快速动态匹配追踪算法，实现了匹配追踪高分辨率时频谱计算，然后通过构建楔形砂体模型深入剖析了薄砂储层的时域和频域反射特征，发展了利用时频谱分量来指示砂体尖灭点位置的方法。该方法综合了薄层调谐能量和薄层反射向高频移动的特征，采用相对高频的瞬时谱分量的异常高值识别薄砂尖灭点，相比单独使用振幅属性的尖灭点检测方法精度有所提高。模型试算验证了利用时频谱分量来指示砂体尖灭点位置的方法识别薄层砂体尖灭点的有效性,在实际资料的薄砂尖灭识别应用中，该方法取得了较好的效果，有效地证实了其实用性和可靠性。

匹配追踪； 尖灭线检测； 薄层砂体； 时频分析

0 引言

随着勘探形势的日益复杂，复杂岩性圈闭边界的精确落实以及砂泥岩尖灭线检测已经成为国内、外研究难点。由于地震信号在地下复杂介质传播过程中受到大地滤波作用，地震信号是一个带限信号，地震资料的垂直分辨率较低。当储层厚度小于调谐厚度(约1/4波长)时，干涉作用导致相邻地层界面的有效反射形成单个复合波，因此，时间域振幅信息难以准确有效地识别薄砂储层的空间尖灭位置。

为了克服地震资料分辨率的局限性，高分辨率地震采集和处理方法被不断提出并得到了广泛应用(宽频地震采集、时频谱白化、吸收衰减补偿(Gabor反褶积)以及稀疏反褶积等)，这为利用地震资料识别薄层提供了较好的数据基础[1-6]。基于地震属性的薄层分析方法在地震解释领域发展迅速，其主要包括振幅类、频率类、相位类、相干类以及几何类等属性。作为一种数据驱动的解释手段，地震属性薄层分析需明确属性与薄层结构之间的映射关系[7]。李国发等[8]研究了基于模型的薄互层振幅属性分析及应用，总结出了薄互层砂体累计厚度和反射振幅之间的定量关系。薄层调谐能量对应的位置十分接近真正的尖灭点，且此位置在瞬时谱剖面会形成亮点，容易被识别和追踪，因此可利用瞬时谱特征指示尖灭位置。

地震信号本质上是一种非平稳信号，即信号所含的频率分量随着时间的变化而变化，地层厚度变化和吸收衰减效应等都会引起频率异常现象。针对非平稳地震信号的时频分析方法主要包括短时傅立叶变换(STFT)、连续小波变换(CWT)、S变换(ST)以及Wigner-Ville分布等，这些常规的时频分析方法虽然计算效率高，但不能满足识别砂体尖灭线的精度要求。相比传统的时频分析方法，基于匹配追踪的时频分析能够较好地适应信号本质结构特征，对地层的识别能力更强，可用于对薄砂岩储层尖灭的识别[9-12]。何胡军等[13]通过基于匹配追踪算法的子波分解技术识别薄互层储层；相比张繁昌等[14]研究的基于匹配追踪瞬时谱计算的于三角洲砂体尖灭线识别，笔者结合实际地质沉积模式构建地层薄砂互层模型，对复杂沉积模式下的尖灭线识别进行了更全面地研究分析。

笔者在前人研究的基础上，首先研究了基于Morlet小波的快速匹配追踪算法，并通过理论信号和地震信号，验证了几种常规时频分析方法和匹配追踪时频谱计算结果的时频分辨率特征差异。其次，通过构建楔形模型和薄砂互层模型研究了薄层砂体的振幅和频率响应特征，分析发现薄层调谐振幅位置已十分接近真正砂体尖灭点；在此基础上，综合薄层反射位置的频率域响应信息，利用时频谱分量来指示砂体尖灭点位置，相比单独使用振幅属性的尖灭点检测精度有所提高。最后，通过对某海域实际资料砂体沉积尖灭线的检测，验证利用时频谱分量来指示砂体尖灭点位置的方法的有效性。

1 砂体尖灭线识别原理

1.1 匹配追踪时频分析

时频分析方法主要可分为三大类：①窗口限制的线性时频表征方法；②非线性时频分布方法；③反演谱分解。第一类线性时频表征的分辨率主要取决于窗口长度的选择， STFT仅具有单一的时频分辨特征，而CWT、ST等在实现过程中具有自适应的分辨能力；第二类非线性时频分布由于不涉及时间窗口的选择，因此不存在分辨率的问题，如Wigner-Ville分布(WVD)具有较高的时频分辨率，但其存在交叉项干扰的问题。匹配追踪算法(MP)可实现信号的自适应分解，通过将信号在超完备时频原子库中进行投影，将地震信号表示为匹配原子的线性组合形式。

时频原子库是一系列时频原子的集合，也称为字典，常用的时频原子是对窗函数g(t)进行拉伸、平移与调制得到。由高斯窗函数产生的时频原子具有最高的联合时频分辨率，Morlet小波作为一种高斯窗函数调制的时频原子被广泛应用于地震勘探中，Morlet小波字典具体表达式如式(1)所示。

D={mγ(t)}γ∈Γ={mγ=(u,ω,φ)(t)}γ∈Γ

(1)

式中：γ=(u,ω,φ)为时频原子的调制参数；u为中心时间延迟；ω为调制频率；φ为调制相位。字典D中的原子是非正交的且是过完备的，信号f在D中的分解是稀疏的，其分解可以等价为式(2)所示的优化问题，通过控制迭代次数δ和迭代阈值ε实现信号的自适应投影分解。

(2)

式中：a为稀疏表示系数向量;f为原始地震信号向量。

匹配追踪算法是一种贪婪迭代算法，能够很好地对式(2)进行求解，即估计信号f在字典D中原子上的投影。首次迭代后信号f可以被分解为沿mγ1方向及垂直方向分量之和：

f=[f,mγ1]mγ1+R1f

(3)

式中：[f,mγ1]为信号f与原子mγ1的内积，即信号f沿mγ1方向的分量；R1f为投影近似后的残差信号，即信号f沿与mγ1垂直方向的分量。为使残差信号尽可能小，选择最佳原子应满足式(4)。

|[f,mγ1]|=max{|[f,mγn|,n=1，2，L}

(4)

利用同样的方法继续对残差进行分解，假设迭代算法已经进行了N次，得到残差信号为RNf，信号f可以分解为：

(5)

式中：mγn(u,ω,φ)为时频匹配原子；an为时频原子的振幅系数。

针对匹配追踪算法巨大计算量的问题，考虑原始信号的瞬时振幅、瞬时频率和瞬时相位等先验信息[14]，进一步约束频率和相位的搜索半径，从而提高匹配追踪的计算效率。因此每一次搜索方式可以简单表述为式(6)。

γn={u0,ω∈U[ω(u0),δω],φ∈U[φ(u0),δφ]}

(6)

式中：u0={t0:A(t0)=max[A(t)]}为最大振幅包络对应的时间；A(t)、φ(t)和ω(t)分别表示地震信号的瞬时振幅、瞬时相位及瞬时频率信息；U[ω(u0),δω]、U[φ(u0),δφ]为频率和相位搜索邻域；δω、δφ为相应参数搜索半径。按照上述快速匹配追踪原理即可实现地震信号的自适应分解，借助Wigner-Ville分布的高时频分辨率特征，则其时频能量谱可以表示为匹配原子WVD分布的叠加形式

(7)

时频谱Wf(t,ω)不但继承了WVD分布的高分辨率特征，而且消除了WVD分布存在的交叉项干扰的问题。

我们设计了由10个Morlet小波组成的理论信号，针对匹配追踪算法的稀疏重构能力进行了测试，图1中黑色曲线Syn为理论信号，红线Rec为匹配重构信号，Res为匹配残差信号，1st-10th分别表示逐次迭代原子的结果。由此可见，匹配重构信号与理论信号基本保持一致，可较好地识别出合成信号中不同频率、不同相位的Morlet小波原子。

图1 理论信号匹配追踪信号重构过程Fig.1 Reconstruction process of theoretic signal matching pursuit

图2 地震信号时频分析分辨率特征对比Fig.2 Characteristics comparison of different spectrum analysis methods(a)信号；(b)MP；(c)CWT；(d)STFT

为了验证基于Morlet小波的动态匹配追踪WVD时频表征方法的应用效果，通过对1D地震信号MP时频谱与连续小波变换(CWT)及短时傅立叶变换(STFT)的时频谱对比分析，由图2可以看出，MP时频谱的时间和频率分辨率均要高于常规时频分析方法。MP高分辨率时频谱分解方法为薄层砂岩尖灭点识别奠定了数据基础，有助于提高砂泥岩尖灭线的识别精度。

1.2 砂体尖灭线识别原理

为了研究薄层反射信号的振幅和频率响应特征，我们采用25 Hz的雷克子波对楔形地层模型进行了正演(图3、图4)。

由图4(a)和图4(b)可见，单楔体情况下，反射最大振幅(均方根振幅)，随砂体厚度的增大呈现先增大后减小最后达到稳定的变化特征，瞬时频率则呈现先减小后增大最后趋于稳定的变化特征。此外，当薄层厚度小于调谐厚度时(图4中黑线和红线指示位置，砂体厚度为子波波长的1/4)，最大振幅与薄层厚度呈正相关关系，瞬时频率与薄层厚度呈负相关关系。

由图4(c)和图4(d)可见，薄互层的振幅和频率随砂体总厚度的变化趋势与单楔形模型基本一致。由于地震资料受到“子波叠印”的影响[15]，仅仅利用地震振幅信息难以识别出可靠的砂岩尖灭位置；瞬时频率和瞬时相位属性在实际地震资料中不稳定存在[16]，因此均难以提供较为准确和可靠的砂岩尖灭点识别效果。结合上述，要综合薄层位置的反射振幅和频率信息发展时频域砂体尖灭线识别方法，利用时频谱分量来指示砂体尖灭点位置。时频域地震信号比振幅属性的尖灭点检测精度要高，它克服了单独使用瞬时频率信息进行尖灭点识别的不稳定问题。

图3 楔形体模型及合成地震记录Fig.3 Wedge models and corresponding synthetics(a)单楔形理论模型；(b)单楔形合成地震；(c)双楔形理论模型；(d)双楔形合成记录

图4 楔形体最大振幅和瞬时频率响应分析Fig.4 Maximum amplitude and instantaneous frequency analysis of wedge models(a)单楔形最大振幅；(b)单楔形瞬时频率；(c)双楔形最大振幅；(d)双楔形瞬时频率

基于地震振幅的薄层尖灭点识别方法，仅能够识别至调谐能量位置处，而基于时频谱的尖灭线识别方法，综合利用了调谐尺度内频率随地层厚度减小而增大的信息，时频谱中薄层反射信号的能量团向高频端移动(调谐尺度内)，通过选取高频瞬时谱分量即可提高尖灭点的识别精度。由于地震频带宽度是局限的，不存在无限增大的现象，因此地震频带宽度越宽，薄层尖灭点的识别精确越高。同时，借助高分辨率匹配追踪时频分析方法计算地震信号时频谱，为薄层尖灭点的精确识别提供了可靠的数据基础。

2 模型试算

为了验证基于时频谱分析的砂体尖灭线识别方法的可行性与稳定性，根据我国某南部探区中的地质沉积模式构建了2套纵波阻抗地层尖灭模型。模型顶部为泥岩地层，泥岩地层下方包含两套薄砂互层，泥岩和薄砂互层之间发育一套高阻地层(图5和图6)。模型大小为横向上CDP道数为450，纵向上采样时长为300 ms，采样间隔为1 ms，上覆砂体尖灭位置在CDP100位置。

图6 不同岩性充填纵波阻抗模型以及砂体尖灭点检测结果Fig.6 Impedance models filled with different lithology and corresponding pinch-out point detection results(a)理论地质模型IV；(b)模型IV合成地震记录；(c)模型IV去强反射道集；(d)理论地质模型V ；(e)模型V合成地震记录；(f)模型V去强反射道集；(g)理论地质模型VI；(h)模型VI合成地震记录；(i)模型VI去强反射道集

图5(a)、图5(d)、图5(g)为不同厚度泥岩隔层纵波阻抗模型I、II、III，泥岩隔层厚度依次为3 ms～13 ms、2 ms～8 ms及1 ms～3 ms，砂体厚度变化范围为1 ms～12 ms，高速层厚度为10 ms～13 ms。利用主频为25 Hz的零相位Ricker子波基于褶积模型合成地震记录如图5(b)、图5(e)、图5(h)所示，由图5(b)、图5(e)、图5(h)中可见，砂岩尖灭点位置受到地震波不同程度的模糊作用，仅依靠地震剖面中的相位突变点难以精确识别砂岩尖灭位置。此外，随着泥岩隔层厚度的不断减小，尖灭点的识别精度不断下降(相位突变点：CDP145- CDP150- CDP165)。图5(c)、图5(f)、图5(i)为基于时频谱分析的尖灭点识别结果，分析图5(c)、图5(f)、图5(i)可得，由于时频谱分析综合了薄层的振幅和频率响应信息，尖灭点识别位置相比相位突变点更加可靠，并且精度更高(能量异常点：CDP130- CDP143- CDP155)，相比振幅相位突变点检测到的尖灭点位置向前逼近了约10-15个采用点。

图6(a)、图6(d)、图6(g)为不同岩性充填纵波阻抗模型IV、V、VI，岩性的充填差异主要体现在砂岩和泥岩的纵波阻抗差异上(砂岩为低阻地层)，模型VI中较差砂质代表砂岩和泥岩的纵波阻抗差较小。图6(c)、图6(e)、图6(h)为对应模型IV、V、VI的合成地震记录。分析图6(c)、图6(e)、图6(h)发现，地震剖面上“砂体尖灭点”位置依次在CDP140、CDP150、CDP150，与真实尖灭点位置CDP100仍然存在较大的差距。图6(c)、图6(f)、图6(i)为砂体尖灭点的时频谱检测结果，由图6(c)、图6(f)、图6(i)可见，匹配追踪算法可以有效检测出调谐尺度内比较可靠的尖灭点位置(CDP135、CDP135、CDP135)，并且相比地震剖面中的相位突变点识别精度有所提高。

3 实际资料处理

图7(a)为原始地震剖面I，其中绿色箭头指示位置为地震剖面中相位突变点；图7(b)为相位突变处地震道集的时频谱分析结果，由图7可见MP-WVD时频谱的时间和频率分辨率明显要高于CWT时频谱，因此MP-WVD可以更好地为薄层砂体尖灭点识别提供有效的数据基础。

对比图8(a)和图8(b)可得，两种方法均能有效的识别出尖灭点大致位置；但由于CWT谱时频分辨能力的局限性，尖灭点在地震剖面中的时空展布仍然存在较大的误差，而MP-WVD方法识别结果的分辨能力却有较大程度地提高。

图9(a)为工区内的另一地震剖面，尖灭点无明显相位突变现象，图9(b)为其尖灭点检测结果，从图9(b)中可看出，MP时频谱落实了沉积砂体的尖灭位置(图9(a)中绿色箭头所指处)出现了明显的能量异常现象，该识别效果与地质沉积背景相一致。

针对3D地震资料，开展了各期沉积体边界平面展布预测。图10为各期沉积体边界平面展布范围预测结果。由于该实际工区多期前积砂体受到砂岩厚度、泥岩隔层厚度、砂岩品质以及砂体连通性等因素的影响，存在多种尖灭点类型，不同尖灭点类型需要优选不同频率成分进行分析。因此在落实砂体边界平面展布时，综合了多个优势频率瞬时谱的能量。图10中所示的3期沉积体边界预测结果与地质沉积环境相吻合，从而验证了本文方法的实用性和可靠性。

图7 原始地震剖面I以及尖灭点位置时频谱分析Fig.7 Orignal seismic section I and spectrum analysis corresponding to pinch-out trace(a)原始地震数据I；(b)尖灭点位置时频谱

图8 地震剖面I的CWT谱和MP-WVD谱的尖灭点识别结果Fig.8 Pinch-out detection results of CWT and MP-WVD spectrums corresponding to orignal seismic section I(a)CWT检测结果；(b)MP检测结果

图9 原始地震数据II及时频谱检测结果Fig.9 Orignal seismic section II and spectrum analysis corresponding to pinch-out trace(a)原始地震数据II；(b)MP时频谱检测结果

图10 各期沉积体边界平面展布范围预测结果Fig.10 Predicted sediments distribution map in different stages(a)砂体尖灭线I；(b)砂体尖灭线II；(c)砂体尖灭线III

4 结论

通过对比分析了常规时频分析方法与匹配追踪时频谱计算方法的时频分辨率特征，发现匹配追踪时频谱具有最高的时频分辨率。针对薄层砂体的尖灭线检测，利用地震资料振幅本身均难以达到理想的识别效果。笔者发展了基于瞬时谱的薄层砂岩尖灭线检测方法，该方法综合利用了薄层反射的时间域和频率域响应信息，利用时频谱分量来指示砂体尖灭点位置，比振幅属性的尖灭点检测方法精度更高；同时，该方法建立在高分辨率时频谱计算的基础上，薄层尖灭点的识别精确又有较大的提高。此外，值得关注的是，由于薄层砂体沉积环境的复杂性，薄层砂体边界的平面预测往往需要综合优势频带内多频率瞬时谱分量，才能更好地落实薄层砂体边界的沉积全貌。

[1] CHAI XINTAO, WANG SHANGXU, YUAN SANYI, et al. Sparse reflectivity inversion for nonstationary seismic data [J]. Geophysics, 2014, 79(3): V93-V105.

[2] NGUYEN T. High resolution seismic reflectivity inversion[D]. Houston: University of Houston, 2008.

[3] ZHANG RUI, CASTAGNA JOHN. Seismic sparse-layer reflectivity inversion using basis pursuit decomposition[J]. Geophysics, 2011, 76(6):147-158.

[4] CHOPRA S., CASTAGNA J. P. , PORTNIAGUINE O. Thin-bed reflectivity inversion[C]. 76th SEG Annual International Meeting, 2006：2057-2061.

[5] 贺锡雷，黄德济，贺振华. 薄互层反射系数序列时～频特征研究[J]. 物探化探计算技术，2009，31(3):227-232

HE X L, HUANG D J, HE Z H. Research on time-frequency characteristics of reflection coefficients of thin interbeds[J]. Computing Techniques for Geophysical and Geochemical Exploration, 2009，31(3):227-232.(In Chinese)

[6] 孙雷鸣，曾维辉，方中于. 地震薄层反射系数谱反演算法研究及应用[J]. 物探化探计算技术，2014，36(4):462-470.

SUN L M，ZENG W H,FANG Z Y.Thin-bed reflectivity inversion and seismic application[J]. Computing Techniques for Geophysical and Geochemical Exploration, 2014，36(4):462-470. (In Chinese)

[7] 李国发，岳英，国春香,等. 基于模型的薄互层地震属性分析及其应用[J]. 石油物探，2011，50(2): 144-149.

LI G F,YUE Y,GUO C X, et al. Seismic attributes analysis based on model in thin interbedded layers and its application[J].Geophysical Prospecting for Petroleum，2011，50(2): 144-149. (In Chinese)

[8] 李国发，岳英，熊金良,等. 基于三维模型的薄互层振幅属性实验研究[J]. 石油地球物理勘探，2011，46(1): 115-120.

LI G F,YUE Y,XIONG J L, et al. Experimental study on seismic amplitude attribute of thin interbed based on 3D model[J]. Oil Geophysical Prospecting，2011，46(1): 115-120. (In Chinese)

[9] MARFURT S, ZHANG Z. Matchging-pursuit with time frequency dictionaries[J]. IEEE Transactions on Signal Processing, 1993, 41(12)：3397-3415.

[10] 张繁昌，李传辉，印兴耀. 基于动态匹配子波库的地震数据快速匹配追踪[J]. 石油地球物理勘探，2010，45(5):667-673.

ZHANG F C,LI C H , YIN X Y. Seismic data fast matching pursuit based on dynamic matching wavelet library[J].Oil Geophysical Prospecting，2010，45(5):667-673. (In Chinese)

[11] 张繁昌，李传辉. 地震信号复数域高效匹配追踪分解[J]. 石油地球物理勘探，2013，48(02):171-175.

ZHANG F C,LI C H.Complex domain efficient matching pursuit decomposition of seismic signals[J]. Geophysical Prospecting for Petroleum, 2013，48(02):171-175. (In Chinese)

[12] 刘杰，张忠涛，刘道理,等. 强反射背景下沉积体边界检测及流体识别方法[J]. 石油物探，2016，55(1):142-149.

LIU J Z,ZHANG Z T,LIU D L，et al.Sediment boundary identification and fluid detection for the seismic data with strong background reflections[J]. Geophysical Prospecting for Petroleum,2016，55(1):142-149. (In Chinese)

[13] 何胡军，王秋语，程会明. 基于匹配追踪算法子波分解技术在薄互层储层预测中的应用[J]. 物探化探计算技术，2010，32(6):641-644.

HE H J, WANG Q Y, CHENG M H.The application of wavelet decomposition technique based on matching pursuit algorithm in thin interbedded reservoir prediction[J]. Computing Techniques for Geophysical and Geochemical Exploration, 2010，32(6):641-644. (In Chinese)

[14] 张繁昌，李传辉，印兴耀. 三角洲砂岩尖灭线的地震匹配追踪瞬时谱识别方法[J]. 石油地球物理勘探，2012，47(1): 82-88.

ZHANG F C,LI C H,YIN X Y.Delta fringe line recognition based on seismic matching pursuit instantaneous spectral characteristics[J].Oil Geophysical Prospecting，2012，47(1): 82-88. (In Chinese)

[15] 刘豪, 辛仁臣.油气储层地震综合预测技术与应用[M].北京：中国地质大学出版社，2013.

LIU H,XIN R C.Seismic integrative prediction of reservoir and hydrocarbon[M].Binjing:China University of Geoscience Press,2013.(In Chinese)

[16] 高静怀，陈凤，陈树民.利用地震瞬时谱属性进行薄互层分析[J].煤田地质与勘探，2005,33(5)：67-70.

GAO J H,CHEN F,CHEN S M.Using seismic instantaneous attributes to analyze thin interbeds[J]. Coal Geology & Exploration , 2005,33(5)：67-70. (In Chinese)

Themethodofthinsandpinch-outboundarydetectionviaT-Fspectrumbasedonmatchingpursuit

WANG Ruiliang, ZHANG Wenzhu, LIU Xumin, DONG Guohui, LI Zhiye

(CNOOC China LIMITED-Shenzhen, Shenzhen 518000, China)

The adaptive decomposition of seismic signal can be performed by Matching pursuit algorithm. The fast dynamic matching pursuit algorithm is studied, and the computation of high resolution spectrum of seismic data is accomplished in this study. We analyzes the reflection characteristics in time domain and frequency domain through the wedge models-based forwarding, and then develops the method of detection of pinch-out boundary of thin sand reservoir using the time-frequency spectrum component. The energy of thin layers' reflection moves to higher frequency, which is taken into consideration to detect the thin sand pinch-out point. Compared with the conventional method that only use the amplitude of seismic, the proposed method improves the accuracy of pinch-out point detection. Tests on models proves the validity of the proposed method for identifying the sharp vanishing point of thin sand layers. Finally, good results of pinch-out boundary detection in real case application strengthen the practicability and reliability of our method.

matching pursuit; pinch-out line detection; thin sand bed; time-frequency analysis

2017-08-24 改回日期： 2017-09-05

“十二五”国家科技重大专项(2011ZX05023-002-007)

汪瑞良(1960-)，男，教授级高级工程师，主要从事地球物理方法研究及其在油气勘探开发中的应用，E-mail：wangrl@cnooc.com.cn。

1001-1749(2017)06-0799-09

P 631.4

A

10.3969/j.issn.1001-1749.2017.06.13