二阶同步提取变换的沉积旋回界面定位与追踪识别

李雪英, 王 鑫

(1.东北石油大学 地球科学学院, 黑龙江 大庆 163318; 2.黑龙江省油气藏形成机理与资源评价重点实验室, 黑龙江 大庆 163318)

0 引 言

沉积旋回研究是层序地层学中的一项重要工作,诸多学者已经对沉积旋回的规律做了大量研究[1-4]。地球物理观测数据反映了地层的岩性和物性,包含了大量与沉积旋回有关的信息[5]。通常情况下,由现场测井数据进行沉积旋回的划分,显然,这种方法对于井点附近的划分结果准确度高,但是对于无井点地区的分析,仅能通过已有地质规律推导,这样导致了该地段沉积旋回划分结果不尽人意。地震资料对于地下信息收集比较全面,也没有井点的制约,这种特性有利于更准确地研究地层变化。在地下各类地层介质中,泥岩与砂岩的占比比较大,也比较常见。而且通过砂岩泥岩重叠形式和厚度变化进行划分,可以作为重要指标指导判别沉积体系以及海侵-海退过程,也是使用时频分析方法判定沉积旋回模式的重要条件。时频分析方法判定沉积旋回模式的依据为:薄层时频响应机理中的升频降幅作用[6];旋回性薄互层中的小层厚度递增或递减变化,会相应地引起瞬时频谱峰值频率逐渐减小或增大。

通过时频分析方法识别沉积旋回的最早研究可以追溯到Бродов等提出的СФИ技术[7]。此后,加窗傅立叶变换[8]、小波变换[9-10]、广义S变换[11]和希尔伯特-黄变换[12]等已被开发用于确定沉积旋回。利用时频分析方法确定沉积旋回模式的过程中,越高的时频分辨率越有利于旋回模式的判别。为了实现理想的时频分析目标,一些先进的时频后处理方法成为近年来研究的重点。它主要包括RS[13]、SST[14]以及SET[15]等。经过推导与实验证明,二阶同步提取变换方法可以产生比RS、SST和SET方法能量聚焦更好的时频分析结果[16-17],这一特点非常符合沉积旋回判别的需求。

综上所述,笔者建立多套沉积旋回模型,以模拟地下各种情况下的组合形式,运用二阶同步提取变换对该模型的地震道数据进行时频分析处理,探究沉积旋回模型中,各单元顶界面附近的时频特征,给出沉积旋回顶界面的识别与追踪方法。

1 二阶同步提取变换原理

二阶同步提取变换(SET_2)是基于短时Fourier变换的时频分析后处理方法,其是对同步提取变换进行改良的算法;其主要思想是:首先对信号进行短时傅里叶变换分析,得到信号的时频谱,之后利用求取的二阶瞬时频率和delta函数,求得二阶同步提取算子(SEO_2),最后将算子与短时傅里叶变换结果相乘,求得二阶同步提取变换时频谱[20]。该思路的主要公式如下。

短时傅里叶变换公式为

式中:t——时间;

x——积分的自变量;

g(x-t)——高斯窗函数;

y(x)——输入信号;

Ge(t,ω)——短时傅里叶变换结果。

二阶估算频率求解公式[18]为

将delta函数与二阶估算频率结合,得到二阶同步提取算子,则二阶同步提取变换的表达式为

Te(t,ω)=Ge(t,ω)×δ(ω-ω1(t,ω)),

式中,δ——二阶同步提取算子。

二阶同步提取算子表达式为

2 时频特征与旋回界面定位方法

设计多套不同的沉积旋回单元模型,在泥岩背景下进行波动方程正演模拟,各模型的旋回单元厚度范围hr、厚度递变梯度hg、分割旋回单元的隔层厚度hc等参数,如表1所示。地震子波采用零相位的雷克子波,子波主频为39 Hz;通过正演模拟获取单炮纪录,并在单炮记录中提取零偏移距的地震道数据进行之后的研究。

表1 旋回单元组合模型参数

经过正演模拟之后各模型的时域波形及短时傅里叶变换(STFT)、二阶同步提取变换(SET2)结果,如图1～3所示。当模型中各旋回单元的厚度大,而且旋回单元中各小层厚度变化剧烈时,时域波形振幅在较厚小层处比较大,且波形较为稀疏,通过这个特征可以划分旋回模式及其厚度。当模型中各旋回单元厚度较小时(由模型c所示,其小层厚度均为入射波的1/8波长以下,均小于地震分辨率[19]),通过各模型的时域波形无法识别旋回特征。当模型中个旋回单元厚度较大,但是各小层厚度变化不大时(模型a),各时域波形在厚度较大的小层处振幅大,但是从整体来看,波形稀疏程度趋于一致,仅在旋回模型的顶界面附近出现强反射尖峰。

图1 模型a、b时频特征Fig.1 Time-frequency characterization of models a and b

图2 模型c、d时频特征Fig.2 Time-frequency characterization of models c and d

图3 模型e、f时频特征Fig.3 Time-frequency characterization of models e and f

综合分析短时傅里叶变换结果和二阶同步提取变换结果可知,二阶同步提取变换得到的时频谱能量更加聚焦,具有较高的时频分辨率。从图1～3可以看出,二阶同步提取变换频谱图中的每个频谱变化单元都对应着一个沉积旋回单元,而且对于不同旋回单元厚度和不同旋回单元厚度变化梯度,二阶同步提取变换都可以清晰反映出时频谱中高频到低频的变化特征,从这个方面可以对旋回单元的模式进行判别。短时傅里叶变换频谱图中在旋回单元较厚时频谱变化能够对应沉积旋回单元,然而当旋回单元较薄时,无法准确判断旋回单元。两种方法得到的时频谱都表现出了向模型中小层层厚增大的方向时移的特点。

通过时频分析结果分析总结沉积旋回单元顶界面的特征,当模型中沉积旋回单元的总厚度较大、各小层的厚度变化比较明显时,二阶同步提取变换的高频极值对应着旋回单元的顶界面(模型d、e、f中的SET),当厚度变化梯度比较小时,则需要结合二阶同步提取算子进行辅助判别,可以看到在顶界面附近,二阶同步提取算子呈现出高频突变现象,其拱形形态的尖峰对应着沉积旋回单元的顶界面。

薄层具有升频降幅的作用[6],在图1表示的沉积旋回模型中,由于薄层的升频作用,可以看出频率最高的位置总是指向薄层。薄层降幅作用导致薄层产生的高频信号能量弱,二阶同步提取变换对于弱信号具有较强的提取能力,可以很好地提取薄层产生的高频成分。由模型CSEO2可以看出,在旋回单元厚度阶跃的位置,厚层产生的强能量的低频信息及薄层产生的弱能量的高频信息均得到了清晰的刻画。在识别这类旋回单元时,需要结合二阶同步于是取变换的时频谱和二阶同步提取算子的结果进行判断,最终得出顶界面的位置。

综上分析可知,在二阶同步提取变换得到的时频谱中,沉积旋回单元的顶界面位置对应着时频谱的高频极值部分,可以根据这一性质对沉积旋回单元进行判别与划分,进而预测沉积旋回的厚度。

3 沉积旋回的噪声鲁棒性的判别

在实际的地震信号中,必定会有噪声的干扰,为了使文中提出的方法具有更广泛的适用性,探究二阶同步提取变换方法判断沉积旋回单元的抗噪声能力,选取一种前文建立的地质模型,加入一定能量的高斯白噪声使得该模型能够更好的模拟一般情况下底层情况,使得经验更具有一般性。

由前文可知,时域波形随着噪声能量的增大失去了本身的特征。对于二阶同步提取算子计算结果来说,该结果受噪声影响较大,但是通过综合比对二阶同步提取变换时频谱可知,算子计算结果仍然保留了顶界面处的高频突变现象,呈现为规则的拱形形态,并且该形态与噪声的频率分量有明显的差别,如图4所示,拱形形态在70 Hz左右。而从二阶同步提取变换的时频谱来看,加入各强度的噪声之后,噪声能量在时频谱中有一定的响应,但对于反映沉积旋回模型中各单元小层厚度变化的频谱能量并没有很大的干扰,如图5所示。

图4 二阶同步提取算子抗噪能力试验Fig.4 Noise immunity test of second-order synchronous extraction operator

图5 二阶同步提取变换抗噪能力试验Fig.5 Second-order simultaneous extraction transform noise immunity test

而实际地震探测中,噪声能量远不及实验中所加最大能量。综上所述,通过综合二阶同步提取变换以及二阶同步提取算子的计算结果,能够在噪声的干扰下较为准确的识别沉积旋回单元顶界面位置。

4 沉积旋回界面追踪方法

综合上述分析,沉积旋回单元顶界面在二阶同步提取变换时频谱中,必然表现出高频极值的特点,则可以通过这个高频极值来判别沉积旋回单元的顶界面位置。为了验证这个结论的合理性以及探究顶界面追踪方法,设计一个由三个正旋回单元组合的地质模型(图6),地质模型部分参数,如表2所示。

图6 三个正旋回单元组合厚度模型Fig.6 Combined thickness model of three positive cycle elements

表2 旋回单元组合地质模型部分参数

利用软件制作该模型的合成地震记录剖面(图7),然后使用二阶同步提取变换对正演结果的各地震道数据进行处理。

图7 由图6模型正演的地震剖面局部放大Fig.7 Local enlarged view of seismic profile modeled by model in figure 6

4.1 局限性分析

由图7可知,时域波形可以反映出沉积旋回单元的大致位置,当夹层厚度较大时,沉积旋回单元能够显示明显的分离界面,从时域波形的特征来看,能够明显看出沉积旋回单元的位置,但是单元界面与波形的峰谷对应程度不高;当夹隔层厚度较小时,沉积旋回单元界面的时域波形的特征不明显,无法准确判断各沉积旋回单元位置;另外时域波形易受噪声等因素的影响,很难指定精准的判别标准。由此,应通过时频分析结果结合时域波形综合讨论。

由图7地震道1～130可见,沉积旋回单元厚度不小于波长的一倍时:从时域波形来看,小层厚度较大的位置波形比较稀疏,反之小层厚度较小的位置波形比较密集。如图8中地震道波形98～131所示。时频谱中变化趋势清晰反映了沉积旋回单元的厚度变化趋势,并且各单元的顶界面均对应频率的极大值处,从时频谱中可准确判断沉积旋回单元顶界面位置及沉积旋回模式。

图8 图7剖面不同地震道的二阶同步提取变换结果及自动追踪位置与实际位置对比Fig.8 Second order synchronous extraction and transformation results of different seismic tracks in profile in fig.7 and comparison between automatic tracking position and actual position

如图7中地震道130～170可见,沉积旋回单元厚度在一半到一倍波长之间时:从时域波形来看,沉积旋回单元的波形特征仍比较清晰,但在图8中对应的时频谱中,顶界面的对应效果不佳,结合图6中的二阶同步提取算子计算结果则能够很好的矫正出顶界面的位置。如图7中地震道170～200所示,当沉积旋回单元厚度不大于一半波长时:时域波形趋于对称,顶底界面附近存在明显的强能量振幅;时频谱无法反应旋回单元频率变化趋势,但薄层位置依然对应高频极值,此时,应先根据时域波形限定旋回单元顶界面识别范围,进而追踪时频谱频率极大值位置定位旋回单元顶界面。

对于夹层来说,其在模型中起到了分割各沉积旋回单元的作用,由于夹层的存在且夹层厚度对于界面附近的频谱有一定的影响,需要注意夹层对于追踪结果的干扰。当沉积旋回单元的厚度大于波长的一半时,夹层对于判别旋回模式及识别顶界面没有影响。当沉积旋回单元厚度小于波长的一半并且夹层厚度大于波长的一半时,沉积旋回单元可以通过时域波形分辨;当夹层厚度小于波长的一半时,则需要直接从时频谱中进行判别。

4.2 沉积旋回单元顶界面追踪实验

根据上文总结的结论,制定出沉积旋回单元顶界面的追踪方法,如图9所示。分析时域波形的频率范围,给定阈值进行扫描;根据扫描范围内的时频谱频率极大值的位置,对沉积旋回顶界面进行初次识别追踪;将追踪结果结合地震剖面信息进行分析,对给定误差阈值的追踪结果结合二阶同步提取算计计算结果进行二次矫正,最终得出较为准确的沉积旋回单元顶界面位置。

图9 沉积旋回顶界面识别方法流程Fig.9 Sedimentation spinning back top interface identification method process

文中使用前文所建立模型进行实验验证。首先根据时域波形指定扫面范围为时频谱中大于最大振幅30%的值。首次识别结果在图8中以虚线划出,在图10b的地震剖面上以实现划出,可以看出识别结果能够较为准确的反映出沉积旋回单元的顶界面。但也可以在地震剖面图上看到由大幅度突变点,结合二阶同步提取算子计算结果(图10)对这些突变点附近的地震道顶界面进行二次矫正(图10中虚线所示),第二次矫正结果在图11b中以点划线表示,可以看出经过二次矫正之后,突变点附近的识别结果更贴近于真实情况。

图10 图7剖面二阶同步提取算子计算结果及位置对比Fig.10 Calculation results of second order synchronous extraction operators in profile in fig.7 and comparison of tracking positions

受噪声干扰下的合成地震剖面情况,如图11a所示。实验中加入的噪声强度为10 dB。从图中的识别追踪结果来看,受噪声干扰,追踪结果围绕真实界面有一定的波动情况,但仍然能较清晰的看出界面位置,经过二次矫正之后,识别结果更加贴近真实界面位置,相对误差基本能够保持在30%以内,追踪结果与实际位置基本吻合,如图11a中点划线所示。

综上,文中提出的方法具有较好的应用价值。

5 结 论

(1)当沉积旋回单元的厚度大于波长的一半时,夹层对于判别旋回模式及识别顶界面没有影响。当沉积旋回单元厚度小于波长的一半并且夹层厚度大于波长的一半时,沉积旋回单元可以通过时域波形分辨;当夹层厚度小于波长的一半时,则需要直接从时频谱中进行判别。

(2)二阶同步提取变换经由同步提取变换改进而得,该方法具有更高的时频分辨率,具有良好的抗干扰性能,鲁棒性更好。能够更好的提取弱能量信号,对于沉积旋回单元顶界面的薄层信号具有比一般时频分析方法更好的识别能力,提升了顶界面的定位精度。

(3)二阶同步提取变换计算的时频谱在沉积旋回单元顶界面位置出现稳定的频率极值,或出现一个具有稳定窄带拱形特征的频率极值;当二阶同步提取变换计算的时频谱在旋回单元顶界面频率极值能量弱无法识别时,结合二阶同步提取算子计算结果进行定位调整将提升定位精度。

(4)根据沉积旋回单元顶界面与同步提取变换结果频率极值之间的稳定关系,建立旋回单元顶界面追踪方法;通过模型结果分析,讨论方法的局限性,结合二阶同步提取算子计算结果建立人机交互的反馈调节机制,通过沉积旋回模型顶界面追踪实验证明了方法的可行性与实用性。