基于频率域多阶微分融合的地震频带拓宽方法研究

郭　欣,雍学善,高建虎,刘伟方,李胜军

(1.中国石油天然气股份有限公司勘探开发研究院,北京100083;2.中国石油天然气股份有限公司勘探开发研究院西北分院,甘肃兰州730020;3.中国石油天然气集团公司油藏描述重点实验室,甘肃兰州730020)



基于频率域多阶微分融合的地震频带拓宽方法研究

郭欣1,雍学善2,3,高建虎2,3,刘伟方2,3,李胜军2

(1.中国石油天然气股份有限公司勘探开发研究院,北京100083;2.中国石油天然气股份有限公司勘探开发研究院西北分院,甘肃兰州730020;3.中国石油天然气集团公司油藏描述重点实验室,甘肃兰州730020)

摘要:针对地震资料频带较窄的问题,提出了一种基于频率域多阶微分融合的地震频带拓宽方法,该方法利用微分运算可揭示信号微观特征的性质来拓展高频信息。微分算子在振幅谱上与频率呈平稳的递增关系,微分阶数大小决定了信号高频的提高程度。在测井反射系数振幅谱低频趋势的约束下,将不同阶微分提取的信息进行多次迭代融合,并将其变换到时间域,获得了宽频的地震记录,其高截频可以提升1倍左右,有效地降低了地震子波带限的影响。利用该方法对薄互层数学模型、尖灭体物理模型和实际地震数据进行了试处理,均取得了较好的效果,验证了本方法对于拓展地震资料频带的可行性及有效性。

关键词:多阶微分;频率域;迭代融合;拓宽频带;经验模态分解

宽频带地震资料是精细地震解释的基础,而实际地震资料频带较窄,分辨率较低,难以满足精细地震解释的需求。拓宽地震资料频带的方法有很多,包括反褶积[1-8]、谱白化[9-10]、反Q滤波[11-12]、高频补偿[13]及谱整形[14-20]等,不同方法具有不同的优势和特点,也解决了不同时期地震资料高分辨处理面临的一些难题。但是,随着岩性地层油气藏勘探开发工作的深入,现有的反褶积、谱白化等方法技术因其各自的局限,难以满足薄互层识别和分辨的需求。

微分可以揭示信号的细节或高频信息。赵圣亮等[21]提出了一种采用微分来补偿地震频率的方法,处理后地震剖面的各层能量基本上保持一致且提高了资料的分辨率。云美厚等[22]提出了多级微分频谱整形滤波器的设计思路,针对地震采集仪的模数转换前置频谱进行整形,论述了2,4,6阶微分整形滤波器的理论整形效果,并指出,当微分阶数n≥2时,地震记录频带可扩展到原来的两倍或更宽。然而,孤立地采用某一阶微分并不能真实提高地震资料的分辨率,只是提高了视分辨率。只有将多阶微分的结果进行合理地融合(加权相加),才可能真正地拓展频带宽度。SAJID等[23]对地震记录在时间域进行了2,4,6阶微分变换,并在相同相位条件下相加,得到了高分辨率的地震剖面,其结果优于谱白化算法以及Gabor反褶积的结果。但上述方法受到时间域运算的限制,为了防止高阶微分产生异常高频,只能进行少数几阶微分计算,而且不同地震道拓宽频带后频谱差异较大,可能会影响到地震剖面同相轴的连续性。

本文提出的基于频率域多阶微分融合的频带拓宽方法利用了微分算子在频率域平稳、直观、灵活的特性,将主要运算放在频率域,可将微分阶数增加到十几阶,并拓展到分数(小数)阶。在反射系数振幅谱低频趋势的约束下,对不同阶微分运算提取的频率成分进行迭代融合,获得了较宽频带的地震剖面,较好地压制了高阶微分产生的异常高频。

1方法原理

1.1频率域多阶微分的特点和优势

在地震勘探中,由于地表激发和接收因素的限制,使得地震子波延续时间较长、频带较窄,在剖面上对薄层或薄互层的分辨能力不足。在信号处理中,微分表达信号的变化率,可以将波形中隐含的高频弱信息放大,凸显其细节。如果将地震信息看作物理量中的振幅、其1阶微分表达速度,2阶微分表达加速度,多阶微分虽然还没有明确的物理名称,但却能追溯高频信息通过“同频共振”所表达的薄层或薄互层的地层信息。因此,在地震剖面上被子波干涉所掩盖的薄层或薄互层信息,就可以通过高阶微分将其挖掘和揭示出来。

为了形象阐述信号微分的特点,建立了薄互层模型并进行微分处理,结果见图1。图1a中,左图表示薄互层模型,其中白色为泥岩,黄色为砂岩,假设速度相同(4000m/s),砂泥岩密度比为1∶2,地层厚度由上下两端向中心逐渐变薄,变化间隔为2m,具体为36,34,…,2,4,…,36m,将其转化到时间域。提取模型的反射系数,并与主频为50Hz的雷克子波褶积得到合成地震记录(图1a中的右图)。

图1　薄互层模型及其微分地震记录a 薄互层模型及其合成地震记录; b 时间域实现不同阶微分记录; c 频率域实现不同阶微分记录(压制异常高频后)

对合成地震记录在时间域进行1～4和11阶微分处理(如图1b),微分阶数每增加1阶,相位旋转90°;地震记录的奇数阶微分反映阻抗(或层速度),而偶数阶微分反映地层界面(或反射系数);随着微分阶数的增加,调谐厚度减小,所反映的薄层信息被逐步凸显出来,但当微分阶数达到11阶时,高频异常造成的“毛刺”现象严重影响了波形。

地震记录s(t)的n阶微分振幅谱为:

(1)

式中:f为频率;S(f)为地震记录的振幅谱;S(n)(f)为地震记录的n阶微分振幅谱(简称微分谱)。离散微分算子在频率域并非严格线性变化(图2),其低频近似线性变化,而高频趋于平缓。离散微分算子的这种性质避免了过高地提升高频噪声,同时注意到:随着微分阶数的增加,低频逐渐被压制,而高频逐步得到抬升。

图2　脉冲信号的微分振幅谱

频率域微分运算的优势:①微分阶数在频率域以频率的指数项表达,可以是分数(小数),为多阶微分迭代融合提供了灵活性,而在时间域微分阶数只能为整数;②微分算子在振幅谱上为光滑曲线,改造子波运算更加直观简单;③高阶微分运算产生的高频异常造成的毛刺现象,在时间域难以消除,而在频率域可以在高截频后加衰减函数,以减小高频提升过高的影响(图1c中第11阶微分记录)。

1.2多阶微分融合

提高地震资料分辨率的关键在于保护低频并拓展高频,而单阶微分只能提高信号的主频。如果将不同阶微分的频谱进行合理地融合,可以有效拓展地震记录的频带宽度。

将多阶微分谱直接相加,影响最终频带范围的主要因素为微分阶数,但微分阶数是一个抽象概念,不能有效控制拓宽频带的范围,但微分阶数跟微分后信号的主频可以建立联系,如图3。

图3　微分阶数与微分信号主频的关系曲线

对于某一地震记录,微分阶数越大,微分后信号的主频越高,即微分信号的主频与微分阶数呈单调正相关关系,可以通过曲线拟合得函数φ,使得fm=φ(n),其中fm为微分后信号的主频,则

(2)

函数关系φ由地震子波的实际情况来确定,表1 为单道模型的微分信号主频与微分阶数的对应关系。

表1　微分信号主频与微分阶数的关系

由公式(1)和公式(2)可知:

(3)

若对其加权叠加,可得到叠合谱:

(4)

式中:ai为加权系数;fmin为最小微分信号主频;fmax为最高微分信号主频。为了获取加权系数值,需要建立微分谱融合的标准及融合过程所采用的算法。

首先介绍微分谱融合的标准。在理想情况下,根据褶积理论,地震记录消除子波的影响后为反射系数,可采用测井提取的反射系数振幅谱作为微分融合的迭代标准。地震反射系数并非白噪声序列,而具有偏蓝色特征。为了能从反射系数振幅谱中提取低频趋势项(即反射系数序列的非白噪成分谱),本文采用了WALDEN[24]提出的一阶自回归滑动平均(ARMA)模型方法来拟合,其中AR参数θ和MA参数φ均可以通过标准方法来求得,则其非白噪成分为:

(5)

将其转化到频率域为Rn-w,ARMA曲线如图4a红色虚线。

其次,在融合过程中建立线性方程。已知融合标准为Rn-w,则公式(4)改进为:

(6)

其中,B(·)为求包络函数,在本文中求包络算法采用了先提取微分谱局部极大值点,然后采用经验模态分解求余项[25]的方法。其步骤为:先提取微

分谱序列的局部极大值,并对其进行三次样条插值,将新的序列作为输入信号进行经验模态分解,求取残余分量为低频趋势项,即微分谱的包络。为了满足算法的使用条件,校正某些微分谱的不规则性(非中间高、两端低的形态),基于高斯函数的特点(曲线平滑且权重集中于位置参数附近)对微分谱进行高斯函数约束:

(7)

式中:fm为微分信号的主频。高斯函数的位置参数为微分信号的主频,尺度参数c经试验后确定为拓频后频宽的一半。微分谱的包络如图4b黑色曲线所示。

公式(7)结合公式(6)得到线性方程:

(8)

根据所求取的系数ai经加权叠加得到叠合谱为宽频的振幅谱,结合处理前地震记录的相位谱经反傅氏变换即可获得宽频的地震记录。在实际应用中,应根据处理的需求来调节频宽的大小,即调节参数fmin和fmax,其中fmin为原始地震记录的主频。

图4　多阶微分融合原理a 反射系数振幅谱及其拟合的ARMA曲线; b 微分谱及其包络

图5　拓宽频带后的振幅谱准确性验证过程a 原始合成地震记录振幅谱; b处理后的叠合振幅谱(叠合谱); c 反射系数振幅谱; d 叠合谱与反射系数振幅谱的交会曲线

图5显示了拓宽频带后的振幅谱准确性验证过程,对拓频后地震记录振幅谱及反射系数振幅谱在所圈区间[fmin,fmax]进行交会分析(图5d),两者呈良好的线性关系,拟合度达到了0.964,说明拓频后的振幅谱具有较高的可信度。将拓宽频带后的振幅谱结合地震记录相位谱,经反傅氏变换到时间域(图6)。图6a为合成地震

记录,可识别最小层厚为10m;图6b为本文方法处理过的地震记录,参数fmin=50Hz,对应微分阶数为0,fmax=150Hz,对应微分阶数为16。频谱拓宽后的记录在中心薄层处可分辨4m左右的薄层,分辨能力提高2.5倍,厚层反射的主瓣宽度变窄,分辨能力提高。厚层处相位的增加不可避免,这是由子波旁瓣引起的,如果地震信号频带足够宽,且无噪声存在,这种现象将弱化。本方法将微分融合的标准设定为反射系数谱的低频趋势,目标是识别薄层,相位虽然增加,但振幅值变小。为了减小子波旁瓣影响,我们将地震记录进行道积分处理,然后将其振幅谱经归一化后与叠合谱相加,从而增加低频成分(图6c)。在实际应用中,薄层区域的子波旁瓣基本上可以被强同相轴压制,而在厚层处可被测井合成地震记录识别出来。

图6　处理前后的合成地震记录a 原始合成地震记录; b fmin=50Hz,fmax=150Hz时,拓频后的地震记录; c 合成记录道积分谱归一化后与叠合谱相加经反傅氏变换后的地震记录

2应用效果分析

2.1物理模型试验

采用物理模型模拟碳酸盐岩地层,模型设置了

多个密度异常体,代表溶洞;同时存在一个岩性尖灭体(图7a黄色部分)。图7b为经过处理后的叠加剖面,剖面为从南到北显示,起始道对应CDP490点,终止道对应CDP990点。

模型数据的主频为10Hz,频带范围:5～15Hz;处理采用参数fmin=10Hz,fmax=25Hz,对应最高微分阶数约为6.3,处理后频带范围为5～30Hz。此模型关注的重点是岩性尖灭体的识别。图7b中,由于其频带的限制,子波延续时间较长,相互干涉,在地层尖灭处分辨能力不够;图7c是经过多阶微分融合方法处理后的剖面,同相轴的连续性变强,在岩性尖灭处的分辨率明显提高;为了与模型更好对比,对处理后的结果进行道积分(图7d),道积分剖面上的相对阻抗信息与模型具有较好的一致性。

2.2实际资料应用

实际资料采用柴达木盆地某三维工区抽取的一条二维剖面。应用本文方法前,对叠前时间偏移的叠加剖面做零相位化预处理。剖面的主频为23Hz,频带范围10～35Hz;经过试验选择参数fmin=10Hz,fmax=60Hz,对应最高微分阶数为7.3。本文方法不涉及相位的变化,属于纯振幅运算,为了防止高频过高带来的子波旁瓣的影响,将地震记录道积分的振幅谱归一化后与叠合谱相加,补偿低频成分,减小子波旁瓣影响。

图7　物理模型的处理效果a 物理模型; b 叠后剖面; c 多阶微分融合处理后的剖面; d 多阶微分融合处理后的道积分剖面

图8a是叠前时间偏移剖面,剖面中心由于地震子波频带较窄,薄层反射较弱;图8b是在图8a剖面基础上采用频率域多阶微分融合方法处理后的剖面,剖面中心位置的反射同相轴增强并保持了较好的连续性和信噪比,且与合成地震记录具有较高的一致性,处理前井震相关系数为0.90,处理后相关系变成0.86,虽然相关系数略有下降，但剖面分辨率提高明显，且井震标定较好。

图8　处理前后的剖面a 叠前时间偏移剖面; b 采用多阶微分处理后的剖面

图9a显示了井旁道地震记录拓频前后的振幅谱变化,拓频前的频带范围为10～35Hz,处理后频带得到拓宽,频带范围为10～65Hz。为了弱化子波旁瓣的影响,对地震记录作道积分,将其振幅谱归一化后加入叠合谱,使低频成分得到保留且向低频方向拓展。图9b为拓频后每隔10道相加的振幅谱,可以明显看到振幅谱微弱的偏蓝化现象,同时由于道积分谱的加入,低频成分得以补充。

图9　处理前后振幅谱变化a 井旁道地震记录拓频前后振幅谱对比; b 间隔10道采样多道地震记录拓频前后振幅谱对比

3结束语

研究发现,本文提出的基于频率域多阶微分融合的频带拓展方法具有以下优点:

1) 频率域多阶微分算子具有平稳、直观、灵活的特点,可以实现信号的分数阶微分,且消除异常高频干扰较为简单。

2) 通过信号不同阶微分谱的迭代融合,可以较好地恢复地震记录振幅谱,并且跟反射系数振幅谱具有较高的相关性。

3) 多阶微分结果在频率域经过系数校正融合后,可以将高截频拓展1倍左右,获得宽频带的地震信息,可有效增强识别和分辨薄互层的能力。

4) 尖灭体物理模型和实际资料试处理表明,频率域多阶微分融合方法可以有效拓宽地震资料频带,大幅提高地震资料的分辨率。

参考文献

[1]郭廷超,曹文俊,陶长江,等.时变谱模拟反褶积方法研究[J].石油物探,2015,54(1):36-42

GUO T C,CAO W J,TAO C J,et al.Research on time-varying spectral modeling deconvolution method[J].Geophysical Prospecting for Petroleum,2015,54(1):36-42

[2]王红丽,丁在宇,桂德军,等.震源子波反褶积在海洋地震资料处理中的应用[J].石油物探,2013,52(1):49-54

WANG H L,DING Z Y,GUI D J,et al.Application of seismic source wavelet deconvolution in marine seismic survey[J].Geophysical Prospecting for Petroleum,2013,52(1):49-54

[3]梁东辉,陈生昌.基于L0范数稀疏约束的地震数据反褶积[J].石油物探,2014,53(4):397-403

LIANG D H,CHEN S C.Deconvolution of seismic data based onL0 norm sparse constraint[J].Geophysical Prospecting for Petroleum,2014,53(4):397-403

[4]王嘉松,曹桂荣.分形脉冲反褶积方法[J].地球物理学报,1998,41(1):99-108

WANG J S,CAO G R.The method of fractal impulse deconvolution[J].Chinese Journal of Geophysics,1998,41(1):99-108

[5]蒋立,谭佳,张绪建,等.同态反褶积的子波稳定性方法研究[J].石油物探,2012,51(3):239-243

JIANG L,TAN J,ZHANG X J,et al.Research on the methods for wavelet stability of homomorphic deconvolution[J].Geophysical Prospecting for Petroleum,2012,51(3):239-243

[6]唐博文,赵波,吴艳辉.一种实现谱模拟反褶积的新途径[J].石油地球物理勘探,2010,45(增刊1):66-70

TANG B W,ZHAO B,WU Y H.A new way to realize spectral modeling deconvolution[J].Oil Geophysical Prospecting,2010,45(S1):66-70

[7]李金山,姜秀萍,王修田.主能量脉冲反褶积[J].中国海洋大学学报,2011,41(7/8):119-126

LI J S,JIANG X P,WANG X T.The main energy spiking deconvolution[J].Periodical of Ocean University of China,2011,41(7/8):119-126

[8]王秀槐.纯振幅反褶积[J].石油地球物理勘探,1986,21(1):47-53

WANG X H.Pure-amplitude deconvolution[J].Oil Geophysical Prospecting,1986,21(1):47-53

[9]范小东,曾灏,刘益强.地震资料时空变谱白化处理[J].石油地球物理勘探,1995,30(4):550-555

FAN X D,ZENG H,LIU Y Q.Time-space variant spectrum whitening of seismic data[J].Oil Geophysical Prospecting,1995,30(4):550-555

[10]陈传仁,周熙襄.小波谱白化方法提高地震资料的分辨率[J].石油地球物理勘探,2000,35(6):703-709

CHEN C R,ZHOU X X.Improving resolution of seismic data using wavelet spectrum whitening[J].Oil Geophysical Prospecting,2000,35(6):703-709

[11]赵岩,刘洋,胡光义,等.基于衰减补偿的地震资料高分辨率处理方法[J].石油物探,2014,53(1):38-45

ZHAO Y,LIU Y,HU G Y,et al.High resolution processing method of seismic data based on attenuation compensation[J].Geophysical Prospecting for Petroleum,2014,53(1):38-45

[12]郭建,王咸彬,胡中平,等.Q补偿技术在提高地震分辨率中的应用-以准噶尔盆地Y1井区为例[J].石油物探,2007,46(5):509-513

GUO J,WANG X B,HU Z P,et al.Application ofQcompensation technique in improving seismic resolution: case study of Y1 well area,Junggar basin[J].Geophysical Prospecting for Petroleum,2007,46(5):509-513

[13]丁进杰,戴永寿,张亚南.基于高频补偿方法提高地震资料分辨率的初步研究[J].地球物理学进展,2013,28(6):3214-3221

DING J J,DAI Y S,ZHANG Y N.Improving seismic data resolution based on high frequency compensation[J].Progress in Geophysics,2013,28(6):3214-3221

[14]邬达理.谱整形提高分辨率处理技术及应用效果分析[J].石油物探,2011,50(1):33-37

WU D L.Spectrum-shaping based resolution improvement processing technology and the analysis on its application[J].Geophysical Prospecting for Petroleum,2011,50(1):33-37

[15]孙成禹.谱模拟方法及其在提高地震资料分辨率中的应用[J].石油地球物理勘探,2000,35(1):27-35

SUN C Y.Spectrum modeling method and its application to seismic resolution improvement[J].Oil Geophysical Prospecting,2000,35(1):27-35

[16]孙成禹.地震资料的双约束频谱补偿方法[J].石油物探,2000,39(1):35-41

SUN C Y.Dual constraint frequency spectrum compensation of seismic data[J].Geophysical Prospecting for Petroleum,2000,39(1):35-41

[17]GREGORY B,JEREMY N.Increasing seismic resolution using spectral blueing and colored inversion:Cannonball field,Trinidad[J].Expanded Abstracts of 74thAnnual Internat SEG Mtg,2004:1794-1797

[18]ZHOU B Z,IAIN M,HATHERLY P.Tuning seismic resolution by frequency shift[J].Expanded Abstracts of 17thASEG Geophysical Conference and Exhibition,2004:1-4

[19]NEIL H,SVEN T.Frequency extension,resolution,and sparse inversion[J].Expanded Abstracts of 83thAnnual Internat SEG Mtg,2013:3345-3349

[20]FONS K,STEFFEN B.Broadband seismic data-the importance of low frequencies[J].Geophysics,2013,78(2):WA3-WA14

[21]赵圣亮,王晓明,王建华.一种简易的频率补偿方法[J].石油地球物理勘探,1994,29(2):231-235

ZHAO S L,WANG X M,WANG J H.A simple

method compensating frequency[J].Oil Geophysical Prospecting,1994,29(2):231-235

[22]云美厚,田景文,张国富.多级微分频谱整形滤波器的理论研究[J].石油仪器,1998,12(6):22-24

YUN M H,TIAN J W,ZHANG G F.Theoretical study on multiple-order differential spectrum shaping filter[J].Petroleum Instruments,1998,12(6):22-24

[23]SAJID M,GHOSH D.A fast and simple method of spectral enhancement[J].Geophysics,2014,79(3):V75-V80

[24]WALDEN A T,HOSKEN J W J.An investigation of the spectral properties of primary reflection coefficients[J].Geophysical Prospecting,1985,33(3):400-435

[25]HUANG N E,SHEN Z,LONG S R,et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J].Royal Society of London Proceedings,1998,454(1971):903-995

[26]张韵华,奚梅成,陈效群.数值计算方法与算法[M].北京:科学出版社,2009:130-140

ZHANG Y H,XI M C,CHEN X Q.Numerical calculation method and algorithm[M].Beijing:Science Publishing Company,2009:130-140

(编辑:朱文杰)

Seismic bandwidth broadening method based on multi-order differential fusion in frequency domain

GUO Xin1,YONG Xueshan2,3,GAO Jianhu2,3,LIU Weifang2,3,LI Shengjun2

(1.ResearchInstituteofPetroleumExploration&Development,Beijing100083,China;2.ResearchInstituteofPetroleumExploration&Development-Northwest,Petrochina,Lanzhou730020,China;3.CNPCReservoirDescriptionKeyLaboratory,Lanzhou730020,China)

Abstract:Due to the limited bandwidth of seismic wavelet,the information outside of the wavelet bandwidth is difficult to be distinguished and identified,which leads to the low resolution.To the issue,we present a multi-order differential fusion method in frequency domain to broaden the seismic bandwidth.Differential operators have linear features of increasing monotonically and can improve high frequency and suppress low frequency of signals,namely it has the attribute of frequency spectral decomposition and higher differential order represents higher frequency.We can extract the reflection coefficient from logging data and calculate the low frequency change trend of its amplitude spectrum,which are non-white noise components.Then we add up the differential signal through multiple iterations under the constraint of the non-white noise components.A new broadband seismic record is achieved after back to the time domain from the frequency domain,which has the same trend with the spectrum of reflection coefficient and its high cut-off frequency can be doubled,so it can effectively weaken the influence of the limited bandwidth of the seismic wavelet.It demonstrates that this method can broaden the bandwidth well and improve the resolution of seismic data through the test processing of thin interbedded mathematical model,physical model with pinchout bodies and real seismic data.

Keywords:multi-order differential,frequency domain,iterative fusion,bandwidth broadening,empirical mode decomposition

文章编号:1000-1441(2016)02-0271-09

DOI:10.3969/j.issn.1000-1441.2016.02.013

中图分类号:P631

文献标识码:A

基金项目:国家科技重大专项“大型油气田与煤层气开发”专题(2016ZX05007-006)资助。

作者简介:郭欣(1989—),男,硕士在读,主要从事地震资料处理及储层预测方面的研究。

收稿日期:2015-08-19;改回日期:2015-10-29。

This research is financially supported by the National Science and Technology Major Project of China (Grant No.2016ZX05007-006).