频变子波压缩处理方法、装置、计算机设备和存储介质与流程

1.本发明涉及油气探勘中地震资料处理技术领域，特别涉及一种频变子波压缩 处理方法、装置、计算机设备和存储介质。


背景技术：

2.提高地震资料的分辨率是地球物理勘探的难点，也是地球物理一直探讨的话 题。高分辨率地震资料有助于分辨更薄的地层、更精细地描述小断块、小断层以 及油气藏边缘储层的特点。目前存在的垂向分辨率提高技术包括反褶积、反q滤 波、谱白化以及谱蓝化等。
3.麻省理工地球物理分析组在60年代提出了反褶积，基于robinson褶积模型， 首先求反算子，利用反算子求取反射系数序列。peacork和treitel两人于1969年， 开创了预测反褶积；1975年，基于最大熵分析方法的伯格反褶积被提出；徐伯勋 等人于1992年提出了自适应预测反褶积算法；1994年，王承曙提出了多道预测 反褶积。脉冲反褶积作为另一种反褶积技术，是估计反射系数、提高地震资料分 辨率最常规也是最成熟的方法。但是脉冲反褶积主要有两个限制条件：需要假设 子波是最小相位的；需要假设反射系数是白噪；针对这两个问题，后人对脉冲反 褶积做了相应的改进：基于混合相位未知脉冲的最小平方反褶积技术、基于扫描 子波相位的混合相位反褶积技术、分形反褶积方法以及最小熵反褶积方法等等。 velis提出了随机稀疏脉冲反褶积技术，该方法需要同时估计反射系数的位置和反 射系数的幅值。为了规避地震子波最小相位及反射系数白噪的假设，ulrych提出 了同态反褶积；2003年，高少武等采用l1模、parsimony准则和最小熵准则选择 最优子波对同态反褶积进行改进。凌云等在同态理论的基础上开发了零相位同态 反褶积。wiggins于1985年提出了最小熵反褶积，这种方法突出了厚层反射，可 以从地震振幅异常得到尽可能多的有用信息。最小熵反褶积的算法核心是用一个 线性算子来最大程度的增强地震数据的尖脉冲特征，使反射信号的熵最小。由于 子波随着地层深度的增加而发生振幅和相位的变化，非稳态子波反褶积应运而生。margrave于1998年最开始提出了时变反褶积，随后在2011年其通过非稳态反褶 积利用gabor变换估计反射系数，突破了传统的时间域或频率域反褶积模式。以 吸收衰减模型为基础的吸收衰减补偿技术，对大地滤波引起的振幅衰减和相位畸 变进行补偿和校正，补偿效果较依赖于q值精度和资料与模型的匹配度。吸收补 偿技术常用的方法是反q滤波，即在合理估算地下q场分布的前提下，按照理论 衰减模型(如kjartansson模型、futterman模型)，对各频率成分进行振幅补偿 和相位校正。但由于其在振幅补偿方面本质上不稳定，现阶段依然有许多科研人 员在解决该问题上攻坚克难。地震波在传播过程中其振幅呈指数衰减，在传播到 深部位置时，有效信号的能量便淹没于背景噪声之中，而一般的反q滤波方法是 对整体能量的指数抬升，因此深部噪声的能量得到急剧放大，从而导致不稳定现 象。即便在无噪情况下，由于计算精度问题，数值误差的产生及指数型放大也会 带入不稳定现象。
4.然而以上提高地震资料分辨率方法，在补偿地震资料高频成分的同时，也会 对地
震资料的低频成分造成一定破坏，从而不能有效拓宽地震频带。而地震资料 的垂向分辨率与地震频带息息相关，如何更加有效地保护地震数据的低频成分， 成为提高地震资料分辨率方面一个至关重要的环节。


技术实现要素：

5.基于此，有必要针对上述技术问题，提供一种频变子波压缩处理方法、装置、 计算机设备和存储介质。
6.一种频变子波压缩处理方法，包括：
7.从三维地震工区中提取零相位统计子波；
8.基于所述的零相位统计子波的频谱特征，确定保低频频变尺度因子的关键参 数；
9.基于所述的保低频频变尺度因子，结合三次样条插值技术，构建宽频子波；
10.基于所述的零相位统计子波和宽频子波，计算保低频拓频滤波器；
11.基于所述的保低频拓频滤波器对所述三维地震工区中的每一道地震数据进 行频率域拓频处理，得到拓频处理后的地震数据。
12.在一个实施例中，所述从三维地震工区中提取零相位统计子波的步骤包括：
13.从三维地震工区中，每隔k道提取一道地震数据，并计算该道地震数据的自 相关；
14.基于地震数据自相关的振幅谱与子波的能量谱的关系，对所述提取地震数据 自相关的振幅谱的开方进行平均，进而得到统计子波的振幅谱；
15.基于统计子波为零相位的假设，根据所述统计子波的振幅谱计算得到统计子 波的频谱。
16.在一个实施例中，所述保低频频变尺度因子的计算公式为：
[0017][0018]
其中，a(f)表示保低频频变尺度因子；f
l
和fr分别表示左参考频率和右参 考频率，且f
l
＜fr，一般的，f
l
可选择零相位统计子波的峰值频率，fr可选择f
l
2； a
max
是一个大于1的数，表示最大尺度因子；f
nyq
表示奈奎斯特频率。
[0019]
在一个实施例中，所述基于所述的保低频频变尺度因子，结合三次样条插值 技术，构建宽频子波的步骤包括：
[0020]
基于所述的保低频频变尺度因子，与原频率坐标轴进行点乘操作，得到频变 拉伸处理后的新频率坐标轴；
[0021]
利用三次样条插值技术，将所述的零相位统计子波的频谱，从新频率坐标轴 重采样到原频率坐标轴，即可得到宽频子波的频谱。
[0022]
在一个实施例中，所述基于所述宽频子波频谱和所述统计子波的频谱，计算 得到保低频拓频滤波器的步骤中，所述保低频拓频滤波器的计算式为：
[0023]
[0024]
其中，表示统计子波的频谱的共轭；表示如果向量表示如果向量中的元素比λ大，则不变，反之，则取λ；λ表示正则化系数，λ＝μ
·
χ，μ表 示白噪系数，χ表示向量中元素最大的值。
[0025]
一种频变子波压缩处理装置，包括：
[0026]
零相位统计子波提取模块，用于从三维地震工区中提取零相位统计子波；
[0027]
保低频频变尺度因子确定模块，用于基于所述的零相位统计子波的频谱特征， 确定保低频频变尺度因子的关键参数；
[0028]
宽频子波构建模块，用于基于所述的保低频频变尺度因子，结合三次样条插 值技术，构建宽频子波；
[0029]
保低频拓频滤波器计算模块，用于基于所述的零相位统计子波和宽频子波， 计算保低频拓频滤波器；
[0030]
地震数据拓频处理模块，用于基于所述的保低频拓频滤波器对所述三维地震 工区中的每一道地震数据进行频率域拓频处理，得到拓频处理后的地震数据。
[0031]
在一个实施例中，所述零相位统计子波提取模块包括：
[0032]
地震数据提取单元，用于从三维地震工区中，每隔k道提取一道地震数据， 并计算该道地震数据的自相关；
[0033]
子波振幅谱获取单元，用于基于地震数据自相关的振幅谱与子波的能量谱的 关系，对所述提取地震数据自相关的振幅谱的开方进行平均，进而得到统计子波 的振幅谱；
[0034]
统计子波的频谱获取单元，用于基于统计子波为零相位的假设，根据所述统 计子波的振幅谱计算得到统计子波的频谱。
[0035]
在一个实施例中，所述宽频子波构建模块包括：
[0036]
新频率坐标轴获取单元，包括基于所述的保低频频变尺度因子，与原频率坐 标轴进行点乘操作，得到频变拉伸处理后的新频率坐标轴；
[0037]
宽频子波的频谱获取单元，包括利用三次样条插值技术，将所述的零相位统 计子波的频谱，从新频率坐标轴重采样到原频率坐标轴，即可得到宽频子波的频 谱。
[0038]
一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，其 特征在于，所述处理器执行所述计算机程序时实现以下步骤：
[0039]
从三维地震工区中提取零相位统计子波；
[0040]
基于所述的零相位统计子波的频谱特征，确定保低频频变尺度因子的关键参 数；
[0041]
基于所述的保低频频变尺度因子，结合三次样条插值技术，构建宽频子波；
[0042]
基于所述的零相位统计子波和宽频子波，计算保低频拓频滤波器；
[0043]
基于所述的保低频拓频滤波器对所述三维地震工区中的每一道地震数据进 行频率域拓频处理，得到拓频处理后的地震数据。
[0044]
一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理 器执行时实现以下步骤：
[0045]
从三维地震工区中提取零相位统计子波；
[0046]
基于所述的零相位统计子波的频谱特征，确定保低频频变尺度因子的关键参 数；
[0047]
基于所述的保低频频变尺度因子，结合三次样条插值技术，构建宽频子波；
[0048]
基于所述的零相位统计子波和宽频子波，计算保低频拓频滤波器；
[0049]
基于所述的保低频拓频滤波器对所述三维地震工区中的每一道地震数据进 行频率域拓频处理，得到拓频处理后的地震数据。
[0050]
上述频变子波压缩处理方法、装置、计算机设备和存储介质，通过构建保低 频频变尺度因子，使得保低频频变尺度因子在左右参考频率区间内为一个接近正 弦且随频率变化的函数，进而使得频变子波压缩的过程中不仅可以补偿子波的高 频成分，还可以有效保护子波的低频成分，最终通过结合统计子波与频变压缩后 的宽频子波构建保低频拓频滤波器，对地震数据进行逐道滤波，进而可以获取高 分辨率的地震资料。该发明在叠后地震资料提高分辨率方面具有较高的实用价值。
附图说明
[0051]
图1a为一个实施例中频变子波压缩处理方法的流程示意图；
[0052]
图1b为一个实施例中频变子波压缩处理方法的流程示意图；
[0053]
图2为一个实施例中频变子波压缩处理装置的结构框图；
[0054]
图3为一个实施例中计算机设备的内部结构图；
[0055]
图4为另一个实施例中频变子波压缩处理方法的流程示意图；
[0056]
图5为一个实施例中的验证本技术的理论数据示意图；
[0057]
图6为一个实施例中的理论数据统计子波与宽频子波对比示意图；
[0058]
图7为一个实施例中的理论数据的统计子波与宽频子波的振幅谱对比示意图；
[0059]
图8为一个实施例中的宽频子波合成地震记录与原合成地震记录提高分辨率 后的对比示意图；
[0060]
图9为一个实施例中的实际地震数据统计子波与宽频子波对比图；
[0061]
图10为一个实施例中的实际数据统计子波振幅谱与宽频子波振幅谱对比示 意图；
[0062]
图11为一个实施例中的实际地震剖面与经过保低频频变子波压缩提高分辨 率后的地震剖面对比示意图；
[0063]
图12为一个实施例中的际地震数据的振幅谱与经过保低频频变子波压缩提 高分辨率后的地震数据的振幅谱对比示意图。
具体实施方式
[0064]
为了使本技术的目的、技术方案及优点更加清楚明白，以下结合附图及实施 例，对本技术进行进一步详细说明。应当理解，此处描述的具体实施例仅仅用以 解释本技术，并不用于限定本技术。
[0065]
实施例一
[0066]
一种频变子波压缩处理方法，如图1a所示，其包括：
[0067]
步骤110，从三维地震工区中提取零相位统计子波。
[0068]
本步骤中，从三维地震工区中提取零相位统计子波的频谱。
[0069]
具体地，根据三维叠后地震工区的大小，每隔k道抽取一道叠后地震数据dj， 其中j＝1,2,
…
m，m表示从三维叠后地震工区中抽取的总道数；
[0070]
计算抽取的每一道地震数据的自相关aj，其中j＝1,2,
…
m；
[0071]
对抽取的每一道地震数据的自相关做傅里叶变换fft(aj)，并根据地震数据 自相关的振幅谱等于子波的能量谱这一关系，通过如下公式计算统计子波的振幅 谱：
[0072][0073]
表示统计子波的振幅谱，real(*)表示对复数*取实部，imag(*)表示对 复数*取虚部。
[0074]
假设统计子波为零相位，利用如下公式得到统计子波的频谱：
[0075]
步骤120，基于所述的零相位统计子波的频谱特征，确定保低频频变尺度因 子的关键参数。
[0076]
本步骤中，基于所述的零相位统计子波的频谱特征，确定保低频频变尺度因 子的关键参数，进而根据保低频频变尺度因子的关键参数，构建保低频频变尺度 因子。
[0077]
本实施例中，保低频频变尺度因子的计算式为：
[0078][0079]
其中，其中，a(f)表示保低频频变尺度因子；f
l
和fr分别表示左参考频率 和右参考频率，且f
l
＜fr，一般的，f
l
可选择零相位统计子波的峰值频率，fr可 选择f
l
2；a
max
是一个大于1的数，表示最大尺度因子；f
nyq
表示奈奎斯特频率。
[0080]
保低频频变尺度因子相对常规尺度因子的优势在于：保低频频变尺度因子是 随频率进行变化的；在低频段(0～f
l
)，尺度因子恒为1，即不对频率轴进行拉伸， 确保宽频子波的低频成分与提取的统计子波相一致，从而保护拓频数据的低频成 分；在高频段(fr～f
nyq
)，尺度因子恒为a
max
，即对频率轴进行a
max
倍的拉伸，从 而可以拓宽宽频子波的高频成分；在中频段(f
l
～fr)，尺度因子用一个类正弦函 数表示，确保从1到a
max
是平缓变化的。
[0081]
该公式的构建思路如下：假设时间域子波和频率域子波是对应的傅里叶变换 对则根据傅里叶尺度变换性质得到：
[0082][0083]
其中，φ表示相位谱，可以看出当尺度因子a》1时，地震子波被压缩，其频 率域的振幅谱拉伸，高频成分增加，分辨能力增加；当0＜a《1时，地震子波被 拉伸，子波频谱压缩，高频成分降低，分辨能力降低。通过利用傅里叶尺度变换 性质，可以有效压缩子波，进而提高地震数据的分辨率。然而常规子波压缩变换 中，尺度因子a是一个常数，不随频率f发生变化，因此造成拓频后的地震数据 的低频成分受到压制。
[0084]
针对此问题，通过构建保低频频变尺度因子公式(3)，可以有效保护地震 数据的低频成分。当0≤f＜f
l
时，频变尺度因子a＝1，即不对地震数据的低频成 分进行尺度校正，
从而有效保护低频成分；当f
l
＜f≤f
nyq
，频变尺度因子a＝a
max
， 即对地震数据的高频成分进行尺度校正，从而在高频段进行有效拓频；当f
l
≤f≤fr时，尺度校正因子的值在(1,a
max
)之间，且随频率变化，如公式(3)所 示是一个类似正弦的函数，该函数保证了频变尺度因子a(f)是个光滑的函数，防 止在频变子波压缩时产生异常的高频噪声。
[0085]
步骤130，基于所述的保低频频变尺度因子，结合三次样条插值技术，构建 宽频子波。
[0086]
在一个实施例中，所述基于所述的保低频频变尺度因子，结合三次样条插值 技术，构建宽频子波的步骤包括：基于所述的保低频频变尺度因子，与原频率坐 标轴进行点乘操作，得到频变拉伸处理后的新频率坐标轴；利用三次样条插值技 术，将所述的零相位统计子波的频谱，从新频率坐标轴重采样到原频率坐标轴， 即可得到宽频子波的频谱。
[0087]
本实施例中，首先通过下述式子计算子波压缩变换后的频率向量：
[0088]fnew
＝f
·aꢀꢀ
(4)
[0089]
其中，f
new
表示子波压缩变换后的频率向量；f表示原始频率向量，
[0090]
f＝[0 δf 2δf
ꢀ…ꢀfnyq
]，δf表示频率间隔；
[0091]
a＝[a(0) a(δf) a(2δf)
ꢀ…ꢀ
a(f
nyq
)]。
[0092]
采用三次样条插值技术，将所述子波压缩变换后的频率向量重采样到所述原 始频率向量下，得到压缩变换后的宽频子波频谱。
[0093]
具体地，子波压缩变换后的频率向量为f
new
，其对应的宽频子波的频谱是为了得到原始频率向量下的宽频子波频谱，需要利用三次样条插值技术将子波压 缩变换后的频率向量f
new
对应的频谱p0重采样到原始频率向量f下，得到作为压缩变换后的宽频子波频谱。
[0094]
步骤140，基于所述的零相位统计子波和宽频子波，计算保低频拓频滤波器。
[0095]
本步骤中，基于所述的零相位统计子波的频谱和宽频子波，计算保低频拓频 滤波器。
[0096]
本实施例中，所述保低频拓频滤波器的计算式为：
[0097][0098]
其中，表示统计子波的频谱的共轭；表示如果向量表示如果向量中的元素比λ大，则不变，反之，则取λ；λ表示正则化系数，防止分母过小， 出现计算异常，λ＝μ
·
χ，μ表示白噪系数，通常取0.01，χ表示向量中 元素最大的值。
[0099]
该保低频拓频滤波器相对传统拓频滤波器的优势在于：在相同白噪系数下， 中高频段的拓频能力更强，并且拓频地震数据的振幅谱与期望的宽频子波振幅谱 更接近。
[0100]
公式(5)所示的拓频滤波器的推导过程如下：
[0101]
根据褶积理论，地震记录可以看成是地震子波与地层反射系数褶积而成：
[0102]
d＝w0*r
ꢀꢀ
(6)
[0103]
其中，d表示地震数据，w0表示地震子波，r表示反射系数。考虑到时间域 褶积等价于频率域相乘，因此：
[0104]
[0105]
其中，表示地震数据频谱，表示地震子波频谱，表示反射系数频谱。
[0106]
应该理解的是，反射系数序列是一系列的尖脉冲，地震记录的带限特性主要 来源于地震子波，地震子波的主频高低影响着地震记录的分辨率。假设要得到高 分辨率的地震资料，其对应的宽频地震子波的频谱为因此：
[0107][0108]
其中，表示高分辨率地震资料的频谱。结合(7)式和(8)式得到：
[0109][0110]
上式表明，对地震子波进行处理就是对原始地震数据进行相应的滤波处理。 因此，如果对地震子波进行处理可以改变地震记录的分辨率，关键是求得滤波器 h：
[0111][0112]
针对公式(10)的求解，可以在方程左右两边同时乘以子波频谱的共轭，即针对公式(10)的求解，可以在方程左右两边同时乘以子波频谱的共轭，即因此得到
[0113]
步骤150，基于所述保低频拓频滤波器对所述三维地震工区中的每一道地震 数据进行拓频滤波，得到拓频滤波后的地震数据。
[0114]
在一个实施例中，所述基于所述保低频拓频滤波器对三维地震工区中的每一 道地震数据进行拓频滤波，得到拓频滤波后的地震数据的步骤包括：对所述三维 地震工区中的每一道所述地震数据进行傅里叶变换，得到傅里叶变换后的地震数 据；将傅里叶变换后的地震数据与所述保低频拓频滤波器相乘，得到拓频滤波后 的地震数据；对拓频滤波后的地震数据进行反傅里叶变换，得到提高分辨率后的 地震数据。
[0115]
本步骤中，针对三维地震工区中的每一道地震数据进行傅里叶变换，并与滤 波器h相乘，再进行反傅里叶变换，即可获得提高分辨率后的地震数据。
[0116]
本技术本扩充了常规频变子波压缩提高分辨率方法，通过构建保低频频变尺 度校正因子，并结合傅里叶尺度变换性质与三次样条插值技术，获取宽频子波， 进而构建保低频拓频滤波器，对地震数据进行拓频滤波，最终可以获得高分辨率 的地震资料。
[0117]
在一个实施例中，如图1b所示，步骤110包括：
[0118]
步骤111，从三维地震工区中，每隔k道提取一道地震数据，并计算该道地 震数据的自相关。
[0119]
从三维地震工区中，每隔k道提取一道地震数据，计算提取出来的每一道所 述地震数据的自相关。
[0120]
具体地，三维叠后地震工区中包括多道的地震数据，本步骤中，从三维叠后 地震工区中，每隔k道提取一道地震数据，并对提取出来的地震数据进行自相关 计算。
[0121]
本实施例中，根据三维叠后地震工区的大小，每隔k道抽取一道叠后地震数 据dj，其中j＝1,2,
…
m，m表示从三维叠后地震工区中抽取的总道数。
[0122]
计算抽取的每一道地震数据的自相关aj，其中j＝1,2,
…
m。
[0123]
随后，对每一道所述地震数据的自相关做傅里叶变换，得到傅里叶变换后的 地震数据。
[0124]
本步骤中，对抽取的每一道地震数据的自相关做傅里叶变换即为fft(aj)。
[0125]
步骤112，基于地震数据自相关的振幅谱与子波的能量谱的关系，对所述提 取地震数据自相关的振幅谱的开方进行平均，进而得到统计子波的振幅谱。
[0126]
本步骤中，基于地震数据自相关的振幅谱与子波的能量谱的关系，对所述傅 里叶变换后的地震数据自相关的振幅谱的开方进行求平均，进而得到统计子波的 振幅谱。
[0127]
具体地，本实施例中，参考地球物理中基本假设，假设当反射系数为白噪声 时，子波的能量谱可以近似用地震数据的能量谱表示。因此，通过震数据自相关 的振幅谱与子波的能量谱的关系，能够对对所述傅里叶变换后的地震数据进行振 幅谱转换处理，计算得到统计子波的振幅谱。
[0128]
本实施例中，根据地震数据自相关的振幅谱等于子波的能量谱这一关系，通 过如下公式计算统计子波的振幅谱：
[0129][0130]
其中，表示统计子波的振幅谱，real(*)表示对复数*取实部，imag(*) 表示对复数*取虚部。
[0131]
步骤113，基于统计子波为零相位的假设，根据所述统计子波的振幅谱计算 得到统计子波的频谱。
[0132]
本实施例中，假设统计子波为零相位，利用如下公式得到统计子波的频谱：
[0133][0134]
实施例二
[0135]
请结合图4，本一种频变子波压缩处理方法包括：步骤一：从三维地震工区 中每隔k道抽取一道地震数据；步骤二：计算抽取的每道地震数据的自相关；步 骤三：根据能量谱与振幅谱的关系，计算统计子波的振幅谱；步骤四：基于零相 位假设，获取统计子波的频谱；步骤五：根据傅里叶尺度变换性质，构建保低频 频变尺度因子；步骤六：计算子波压缩变换后的频率向量；步骤七：利用三次样 条插值技术，重采样获取宽频子波的频谱；步骤八：计算保低频拓频滤波器；步 骤九：针对三维地震数据中每一道进行拓频滤波。
[0136]
结合图4，本技术的具体实施过程为：
[0137]
图4为本发明的流程图，通过给定地震波速度、地震波主频确定薄层最大的 时间厚度，结合反射系数奇偶分解理论以及褶积模型理论，构建冗余的楔形字典， 在模型趋势约束和l1范数最小化约束下，利用基追踪分解可以获取一组最稀疏 的解，进而转化为反射系数后，再通过道积分便可以获取高分辨率的波阻抗反演 结果。其具体实施过程如下：
[0138]
步骤一：根据三维叠后地震工区的大小，每隔k道抽取一道叠后地震数据dj， 其中j＝1,2,
…
m，m表示从三维叠后地震工区中抽取的总道数；
[0139]
步骤二：计算抽取的每一道地震数据的自相关aj，其中j＝1,2,
…
m；
[0140]
步骤三：对抽取的每一道地震数据的自相关做傅里叶变换fft(aj)，并根据 地震数据自相关的振幅谱等于子波的能量谱这一关系，通过如下公式计算统计子 波的振幅谱：
[0141]
[0142]
表示统计子波的振幅谱，real(*)表示对复数*取实部，imag(*)表示对复数 *取虚部。
[0143]
步骤四：假设统计子波为零相位，利用如下公式得到统计子波的频谱：
[0144][0145]
步骤五：根据傅里叶尺度变换性质，构建保低频频变尺度因子：
[0146][0147]
其中，a(f)表示保低频频变尺度因子；f
l
和fr分别表示左参考频率和右参考 频率，且f
l
＜fr，一般的，f
l
可选择零相位统计子波的峰值频率，fr可选择f
l
2； a
max
是一个大于1的数，表示最大尺度因子；f
nyq
表示奈奎斯特频率。
[0148]
保低频频变尺度因子相对常规尺度因子的优势在于：保低频频变尺度因子是 随频率进行变化的；在低频段(0～f
l
)，尺度因子恒为1，即不对频率轴进行拉伸， 确保宽频子波的低频成分与提取的统计子波相一致，从而保护拓频数据的低频成 分；在高频段(fr～f
nyq
)，尺度因子恒为a
max
，即对频率轴进行a
max
倍的拉伸，从 而可以拓宽宽频子波的高频成分；在中频段(f
l
～fr)，尺度因子用一个类正弦函 数表示，确保从1到a
max
是平缓变化的。
[0149]
步骤六：计算子波压缩变换后的频率向量：
[0150]fnew
＝f
·aꢀꢀ
(4)
[0151]
其中f
new
表示子波压缩变换后的频率向量；f表示原始频率向量，
[0152]
f＝[0 δf 2δf
ꢀ…ꢀfnyq
]，δf表示频率间隔；
[0153]
a＝[a(0) a(δf) a(2δf)
ꢀ…ꢀ
a(f
nyq
)]。
[0154]
步骤七：子波压缩变换后的频率向量为f
new
，其对应的宽频子波的频谱是为了得到原始频率向量下的宽频子波频谱，需要利用三次样条插值技术将子波压 缩变换后的频率向量f
new
对应的频谱p0重采样到原始频率向量f下，得到作为压缩变换后的宽频子波频谱。
[0155]
步骤八：计算保低频拓频滤波器h：
[0156][0157]
表示统计子波的频谱的共轭；表示如果向量中 的元素比λ大，则不变，反之，则取λ；λ表示正则化系数，防止分母过小，出 现计算异常，λ＝μ
·
χ，μ表示白噪系数，通常取0.01，χ表示向量中元 素最大的值。
[0158]
该保低频拓频滤波器相对传统拓频滤波器的优势在于：在相同白噪系数下， 中高频段的拓频能力更强，并且拓频地震数据的振幅谱与期望的宽频子波振幅谱 更接近。
[0159]
步骤九：针对三维地震工区中的每一道地震数据进行傅里叶变换，并与滤波 器h相乘，再进行反傅里叶变换，即可获得提高分辨率后的地震数据。
[0160]
为了验证本发明的效果，利用图5中的上图给出的反射系数曲线和主频为 30hz的
雷克子波制作合成地震记录，如图5的下图所示。其中上图为反射系数， 下图为合成的理论数据，通过将上图的反射系数与主频为30hz的雷克子波褶积 获取。
[0161]
图6为基于图5中的合成地震数据，通过自相关获取的统计子波(黑线)与 通过保低频频变子波压缩获取的宽频子波(灰线)的对比图，可以看出，相对统 计子波，宽频子波更尖锐，旁瓣更小。
[0162]
图7为图6中合成地震数据的统计子波与宽频子波的振幅谱对比，可以看出， 相对统计子波(黑线)，宽频子波(灰线)的主频更高，频带更宽，并且在低频 段与统计子波保持一致，因此验证了保低频频变尺度因子的有效性。
[0163]
图8给出了宽频子波合成地震记录(黑线)与原合成地震记录经保低频频变 子波压缩提高分辨率后的结果(灰虚线)的对比，可以看出原合成地震记录经过 保低频频变子波压缩提高分辨率后，与宽频子波合成地震记录吻合度十分高，因 此理论数据验证了该方法的可靠性。
[0164]
图9给出了验证本发明的实际地震数据利用自相关获取的统计子波(黑线) 与利用保低频频变子波压缩获取的宽频子波(灰线)对比图。可以看出，相对统 计子波，宽频子波更尖锐，旁瓣更小。
[0165]
图10为实际地震数据的统计子波(黑线)与宽频子波(灰线)的振幅谱对 比，可以看出，相对统计子波，宽频子波的主频更高，频带更宽，并且在低频段 与统计子波保持一致，因此验证了保低频频变尺度因子的有效性。
[0166]
图11为实际地震剖面(上图)与经过保低频频变子波压缩提高分辨率后的 地震剖面(下图)对比图，可以看出经过保低频频变子波压缩提高分辨率后，地 震剖面的垂向分辨率得到了一定的提升，相对原始地震剖面，展示了更加丰富的 细节，并且也能有效凸出弱反射。
[0167]
图12为实际地震数据的振幅谱(黑线)与经过保低频频变子波压缩提高分 辨率后的地震数据的振幅谱(灰线)对比。可以看出，经过保低频频变子波压缩 提高分辨率后，振幅谱的高频成分得到一定的补偿，同时在低频段，与原始地震 数据的振幅谱保持一致，因此实际数据验证了保低频频变尺度因子的有效性。
[0168]
实施例三
[0169]
一种频变子波压缩处理装置，如图2所示，包括：
[0170]
零相位统计子波提取模块210，用于从三维地震工区中提取零相位统计子波；
[0171]
保低频频变尺度因子确定模块220，用于基于所述的零相位统计子波的频谱 特征，确定保低频频变尺度因子的关键参数；
[0172]
宽频子波构建模块230，用于基于所述的保低频频变尺度因子，结合三次样 条插值技术，构建宽频子波；
[0173]
保低频拓频滤波器计算模块240，用于基于所述的零相位统计子波和宽频子 波，计算保低频拓频滤波器；
[0174]
地震数据拓频处理模块250，用于基于所述的保低频拓频滤波器对所述三维 地震工区中的每一道地震数据进行频率域拓频处理，得到拓频处理后的地震数据。
[0175]
在一个实施例中，所述零相位统计子波提取模块包括：
[0176]
地震数据提取单元，用于从三维地震工区中，每隔k道提取一道地震数据， 并计算
该道地震数据的自相关；
[0177]
子波振幅谱获取单元，用于基于地震数据自相关的振幅谱与子波的能量谱的 关系，对所述提取地震数据自相关的振幅谱的开方进行平均，进而得到统计子波 的振幅谱；
[0178]
统计子波的频谱获取单元，用于基于统计子波为零相位的假设，根据所述统 计子波的振幅谱计算得到统计子波的频谱。
[0179]
在一个实施例中，所述保低频频变尺度因子的计算公式为：
[0180][0181]
其中，a(f)表示保低频频变尺度因子；f
l
和fr分别表示左参考频率和右参 考频率，且f
l
＜fr，一般的，f
l
可选择零相位统计子波的峰值频率，fr可选择f
l
2； a
max
是一个大于1的数，表示最大尺度因子；f
nyq
表示奈奎斯特频率。
[0182]
在一个实施例中，所述宽频子波构建模块包括：
[0183]
新频率坐标轴获取单元，包括基于所述的保低频频变尺度因子，与原频率坐 标轴进行点乘操作，得到频变拉伸处理后的新频率坐标轴；
[0184]
宽频子波的频谱获取单元，包括利用三次样条插值技术，将所述的零相位统 计子波的频谱，从新频率坐标轴重采样到原频率坐标轴，即可得到宽频子波的频 谱。
[0185]
在一个实施例中，所述保低频拓频滤波器的计算式为：
[0186][0187]
其中，表示统计子波的频谱的共轭；表示如果向量表示如果向量中的元素比λ大，则不变，反之，则取λ；λ表示正则化系数，λ＝μ
·
χ，μ表 示白噪系数，χ表示向量中元素最大的值。
[0188]
关于频变子波压缩处理装置的具体限定可以参见上文中对于频变子波压缩 处理方法的限定，在此不再赘述。上述频变子波压缩处理装置中的各个模块可全 部或部分通过软件、硬件及其组合来实现。上述各模块可以硬件形式内嵌于或独 立于计算机设备中的处理器中，也可以以软件形式存储于计算机设备中的存储器 中，以便于处理器调用执行以上各个模块对应的操作。
[0189]
在一个实施例中，提供了计算机设备，其内部结构图可以如图3所示。该计 算机设备包括通过系统总线连接的处理器、存储器、网络接口、显示屏和输入装 置。其中，该计算机设备的处理器用于提供计算和控制能力。该计算机设备的存 储器包括非易失性存储介质、内存储器。该非易失性存储介质存储有操作系统和 计算机程序。该内存储器为非易失性存储介质中的操作系统和计算机程序的运行 提供环境。该计算机设备的网络接口用于与其他计算机设备连接通信。该计算机 程序被处理器执行时以实现一种频变子波压缩处理方法。该计算机设备的显示屏 可以是液晶显示屏或者电子墨水显示屏，该计算机设备的输入装置可以是显示屏 上覆盖的触摸层，也可以是计算机设备外壳上设置的按键、轨迹球或触控板，还 可以是外接的键盘、触控板或鼠标等。
[0190]
本领域技术人员可以理解，图3中示出的结构，仅仅是与本技术方案相关的 部分结构的框图，并不构成对本技术方案所应用于其上的计算机设备的限定，具 体的计算机设备可以包括比图中所示更多或更少的部件，或者组合某些部件，或 者具有不同的部件布置。
[0191]
实施例四
[0192]
一种计算机设备，包括存储器和处理器，该存储器存储有计算机程序，该处 理器执行计算机程序时实现以下步骤：
[0193]
从三维地震工区中提取零相位统计子波；
[0194]
基于所述的零相位统计子波的频谱特征，确定保低频频变尺度因子的关键参 数；
[0195]
基于所述的保低频频变尺度因子，结合三次样条插值技术，构建宽频子波；
[0196]
基于所述的零相位统计子波和宽频子波，计算保低频拓频滤波器；
[0197]
基于所述的保低频拓频滤波器对所述三维地震工区中的每一道地震数据进 行频率域拓频处理，得到拓频处理后的地震数据。
[0198]
在一个实施例中，处理器执行计算机程序时还实现以下步骤：
[0199]
从三维地震工区中，每隔k道提取一道地震数据，并计算该道地震数据的自 相关；
[0200]
基于地震数据自相关的振幅谱与子波的能量谱的关系，对所述提取地震数据 自相关的振幅谱的开方进行平均，进而得到统计子波的振幅谱；
[0201]
基于统计子波为零相位的假设，根据所述统计子波的振幅谱计算得到统计子 波的频谱。
[0202]
在一个实施例中，处理器执行计算机程序时还实现以下步骤：
[0203]
基于所述的零相位统计子波的频谱特征，确定保低频频变尺度因子的关键参 数，其中，所述保低频频变尺度因子的计算公式为：
[0204][0205]
其中，a(f)表示保低频频变尺度因子；f
l
和fr分别表示左参考频率和右参 考频率，且f
l
＜fr，一般的，f
l
可选择零相位统计子波的峰值频率，fr可选择f
l
2； a
max
是一个大于1的数，表示最大尺度因子；f
nyq
表示奈奎斯特频率。
[0206]
保低频频变尺度因子相对常规尺度因子的优势在于：保低频频变尺度因子是 随频率进行变化的；在低频段(0～f
l
)，尺度因子恒为1，即不对频率轴进行拉伸， 确保宽频子波的低频成分与提取的统计子波相一致，从而保护拓频数据的低频成 分；在高频段(fr～f
nyq
)，尺度因子恒为a
max
，即对频率轴进行a
max
倍的拉伸，从 而可以拓宽宽频子波的高频成分；在中频段(f
l
～fr)，尺度因子用一个类正弦函 数表示，确保从1到a
max
是平缓变化的。
[0207]
在一个实施例中，处理器执行计算机程序时还实现以下步骤：
[0208]
基于所述的保低频频变尺度因子，与原频率坐标轴进行点乘操作，得到频变 拉伸处理后的新频率坐标轴；
[0209]
利用三次样条插值技术，将所述的零相位统计子波的频谱，从新频率坐标轴 重采
样到原频率坐标轴，即可得到宽频子波的频谱。
[0210]
在一个实施例中，处理器执行计算机程序时还实现以下步骤：所述基于所述 宽频子波频谱和所述统计子波的频谱，计算得到保低频拓频滤波器的步骤中，所 述保低频拓频滤波器的计算式为：
[0211][0212]
其中，表示统计子波的频谱的共轭；表示如果向量表示如果向量中的元素比λ大，则不变，反之，则取λ；λ表示正则化系数，λ＝μ
·
χ，μ表 示白噪系数，χ表示向量中元素最大的值。
[0213]
该保低频拓频滤波器相对传统拓频滤波器的优势在于：在相同白噪系数下， 中高频段的拓频能力更强，并且拓频地震数据的振幅谱与期望的宽频子波振幅谱 更接近。
[0214]
实施例五
[0215]
一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执 行时实现以下步骤：
[0216]
从三维地震工区中提取零相位统计子波；
[0217]
基于所述的零相位统计子波的频谱特征，确定保低频频变尺度因子的关键参 数；
[0218]
基于所述的保低频频变尺度因子，结合三次样条插值技术，构建宽频子波；
[0219]
基于所述的零相位统计子波和宽频子波，计算保低频拓频滤波器；
[0220]
基于所述的保低频拓频滤波器对所述三维地震工区中的每一道地震数据进 行频率域拓频处理，得到拓频处理后的地震数据。
[0221]
在一个实施例中，计算机程序被处理器执行时还实现以下步骤：
[0222]
从三维地震工区中，每隔k道提取一道地震数据，并计算该道地震数据的自 相关；
[0223]
基于地震数据自相关的振幅谱与子波的能量谱的关系，对所述提取地震数据 自相关的振幅谱的开方进行平均，进而得到统计子波的振幅谱；
[0224]
基于统计子波为零相位的假设，根据所述统计子波的振幅谱计算得到统计子 波的频谱。
[0225]
在一个实施例中，计算机程序被处理器执行时还实现以下步骤：
[0226]
基于所述的零相位统计子波的频谱特征，确定保低频频变尺度因子的关键参 数，其中，所述保低频频变尺度因子的计算公式为：
[0227][0228]
其中，a(f)表示保低频频变尺度因子；f
l
和fr分别表示左参考频率和右参 考频率，且f
l
＜fr，一般的，f
l
可选择零相位统计子波的峰值频率，fr可选择f
l
2； a
max
是一个大于1的数，表示最大尺度因子；f
nyq
表示奈奎斯特频率。
[0229]
保低频频变尺度因子相对常规尺度因子的优势在于：保低频频变尺度因子是 随
频率进行变化的；在低频段(0～f
l
)，尺度因子恒为1，即不对频率轴进行拉伸， 确保宽频子波的低频成分与提取的统计子波相一致，从而保护拓频数据的低频成 分；在高频段(fr～f
nyq
)，尺度因子恒为a
max
，即对频率轴进行a
max
倍的拉伸，从 而可以拓宽宽频子波的高频成分；在中频段(f
l
～fr)，尺度因子用一个类正弦函 数表示，确保从1到a
max
是平缓变化的。
[0230]
在一个实施例中，计算机程序被处理器执行时还实现以下步骤：
[0231]
基于所述的保低频频变尺度因子，与原频率坐标轴进行点乘操作，得到频变 拉伸处理后的新频率坐标轴；
[0232]
利用三次样条插值技术，将所述的零相位统计子波的频谱，从新频率坐标轴 重采样到原频率坐标轴，即可得到宽频子波的频谱。
[0233]
在一个实施例中，计算机程序被处理器执行时还实现以下步骤：
[0234]
基于所述宽频子波频谱和所述统计子波的频谱，计算得到保低频拓频滤波器 的步骤中，所述保低频拓频滤波器的计算式为：
[0235][0236]
其中，表示统计子波的频谱的共轭；表示如果向量表示如果向量中的元素比λ大，则不变，反之，则取λ；λ表示正则化系数，λ＝μ
·
χ，μ表 示白噪系数，χ表示向量中元素最大的值。
[0237]
该保低频拓频滤波器相对传统拓频滤波器的优势在于：在相同白噪系数下， 中高频段的拓频能力更强，并且拓频地震数据的振幅谱与期望的宽频子波振幅谱 更接近。
[0238]
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是 可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一非 易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法 的实施例的流程。其中，本技术所提供的各实施例中所使用的对存储器、存储、 数据库或其它介质的任何引用，均可包括非易失性和/或易失性存储器。非易失性 存储器可包括只读存储器(rom)、可编程rom(prom)、电可编程rom (eprom)、电可擦除可编程rom(eeprom)或闪存。易失性存储器可包括 随机存取存储器(ram)或者外部高速缓冲存储器。作为说明而非局限，ram 以多种形式可得，诸如静态ram(sram)、动态ram(dram)、同步dram (sdram)、双数据率sdram(ddrsdram)、增强型sdram(esdram)、 同步链路(synchlink)dram(sldram)、存储器总线(rambus)直接ram (rdram)、直接存储器总线动态ram(drdram)、以及存储器总线动态 ram(rdram)等。
[0239]
以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实 施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的 组合不存在矛盾，都应当认为是本说明书记载的范围。
[0240]
以上所述实施例仅表达了本技术的几种实施方式，其描述较为具体和详细， 但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普 通技术人员来说，在不脱离本技术构思的前提下，还可以做出若干变形和改进， 这些都属于本技术的保护范围。因此，本技术专利的保护范围应以所附权利要求 为准。