在地震勘探领域, 提高地震资料分辨率一直是国内外学者重点关注的研究内容之一, 但是, 通常主要关注地震资料中的高频信息。随着勘探技术的不断发展, 地球物理学家越来越认识到低频信息对于地震勘探的重要性。地震资料中低频信息的传播距离比高频信息的传播距离远得多, 其携带的有效信息也更加丰富。因此, 可以利用低频信息传播距离远、穿透能力强的性质, 有效提高复杂地质构造的成像质量^[1]。对油气资源地震勘探开发而言, 低频信息更为重要, 这是因为在全波形反演中, 低于10Hz的频率成分发挥了更为重要的作用。如果勘探地震资料中含有丰富的低频信息, 将有效提升速度建模和成像的精度。

反褶积作为提高地震资料分辨率最常使用的处理方法, 既可以压缩地震子波、衰减多次波, 还可以用于估计反射系数。最先提出的是时间序列褶积模型, 该模型将地震道看作许多地震子波叠加而成的复合波, 后者由反射系数序列与地震子波褶积而成, 之后多种预测反褶积方法相继被提出^[2]。WIGGINS^[3]提出了对子波的相位无任何要求的最小熵反褶积算法, 通过调整反子波使输出时序方差的模最大, 最终得到反射系数。GRAY^[4]、GODFREY^[5]以及WALDEN^[6]相继对该方法进行了深入研究, 优化了最小熵反褶积算法。LARUE等^[7]和VAN等^[8]提出了两个随机变量相关性最低化条件下的反褶积算法, 该方法通过时频转换实现反褶积处理。ROSA等^[9]将谱模拟算法引入反褶积处理中, 进一步提高了地震资料的分辨率。MASOOMZADEH等^[10]和郭树祥等^[11]分别提出了空间谱白化和分频去噪的处理方法, 有效保护了地震资料中的低频能量。WOODBURN等^[12]利用零相位匹配算子对零相位子波的低频信号进行频率补偿, 提高了低频信号的振幅谱, 实现了对低频成分的补偿。此外, 管路平等^[13]采用预测滤波方法在频率域对地震资料的频谱进行外推, 将丢失的低频信息恢复出来。PEDERSEN等^[14]提出基于测井数据的多变量插值方法, 将测井数据、地震参数和速度数据有效结合, 通过估计地震资料的层速度、层深、地层厚度等, 提高低频信息预测结果的准确性, 该方法明显改善了2~8Hz频率范围的波阻抗信息。韩立国等^[15]、ZHANG等^[16]和张彬彬等^[17]分别提出了基于压缩感知补偿地震低频信息的不同方法。陈祖庆等^[18]和夏红敏等^[19]分别在压缩感知框架下实现了反射系数谱反演和地震数据的插值重构。本文基于前人在低频补偿方面的研究成果, 在压缩感知框架下, 利用反射系数反演, 研究了对地震资料进行低频端和高频端双向拓频以提高地震资料分辨率的方法。

1 反射系数反演基本理论

压缩感知(compressed sensing, CS)是一种利用信号的稀疏性, 寻求欠定方程组稀疏解的方法, 该方法理论由DONOHO^[20]和BARANIUK^[21]等提出。一般情况下, 压缩感知理论由信号的稀疏表示、测量矩阵及稀疏促进算法3部分构成。本文进行了基于反射系数反演的频带双向拓展的提高地震资料分辨率的方法研究。

1.1 压缩感知框架下的反射系数反演方法

褶积模型的基本思想是地震子波经过地下各反射界面时都会产生反射波, 将各个反射波叠加起来就可得到地震记录。含噪地震记录如下:

$ g(t) = w(t) * r(t) + n(t) $

(1)

式中: g(t)为含噪地震记录; w(t)为地震子波; r(t)为地下反射系数; n(t)为随机噪声。通过Fourier变换将(1)式变换到频率域:

$ \mathit{\boldsymbol{G}} = \mathit{\boldsymbol{WR}} + \mathit{\boldsymbol{N}} = \mathit{\boldsymbol{WF}}r + \mathit{\boldsymbol{N}} $

(2)

式中: F是傅里叶变换算子, 为稀疏表示矩阵; G, W, R, N分别是g(t), w(t), r(t), n(t)在频率域的稀疏简化表示矩阵。

在压缩感知理论框架下, 先构造一个合适的测量矩阵, 然后利用信号的稀疏特性, 通过求解稀疏正则化的凸优化问题将完整的地震数据较为准确地重构出来。褶积模型中的反射系数在时域具有稀疏性, 符合压缩感知理论应用的前提条件^[17]。

在压缩感知理论框架指导下, 将(2)式的欠定方程问题转化为求解稀疏规则算子约束下的最小误差能量(J)问题^[19]:

$ \begin{array}{c} \mathit{\boldsymbol{J}} = \frac{1}{2}\left\| {\mathit{\boldsymbol{WR}} - \mathit{\boldsymbol{G}}} \right\|_2^2 + \lambda \left\| r \right\|_1^1 = \frac{1}{2}\left\| {\mathit{\boldsymbol{WF}}r - } \right.\\ \left. \mathit{\boldsymbol{G}} \right\|_2^2 + \lambda \left\| r \right\|_1^1 \end{array} $

(3)

式中: ‖·‖₂²为L₂范数约束的误差函数; W是根据地震子波构造的测量矩阵; ‖·‖₁¹为L₁范数约束的准确解; λ是正则化参数, 用于调节两个范数的权重以确保稀疏性与数据保真之间的平衡, λ越大, 稀疏规则算子占的权重越大, 求得的反射系数稀疏性越好。(3)式第1项采用最小误差能量约束, 在求解过程中不断向完整的地震资料收敛, 第2项采用稀疏规则算子约束, 以获取稀疏性更好的反射系数r。

为了稳定求解, 将(3)式变换如下:

$ \begin{array}{c} \mathit{\boldsymbol{J}} = {(\mathit{\boldsymbol{WF}})^{\rm{T}}}\mathit{\boldsymbol{WF}}r - r{(\mathit{\boldsymbol{WF}})^{\rm{T}}}\mathit{\boldsymbol{G}} - {\mathit{\boldsymbol{G}}^{\rm{T}}}(\mathit{\boldsymbol{WF}}r) + \\ {\mathit{\boldsymbol{G}}^{\rm{T}}}\mathit{\boldsymbol{G}} + \lambda \left\| r \right\|_1^1 \end{array} $

(4)

式中: T表示转置(实数)或共轭转置(复数)。为了求取最优解, (4)式中的J对r的求导如下:

$ \nabla \mathit{\boldsymbol{J}}(r) = {(\mathit{\boldsymbol{WF}})^{\rm{T}}}\mathit{\boldsymbol{WF}}r - {(\mathit{\boldsymbol{WF}})^{\rm{T}}}\mathit{\boldsymbol{G}} + \lambda {\rm{sign}}(r) $

(5)

式中: sign(r)为符号函数, 当r＞0时, sign(r)＞0;当r=0时, sign(r)=0;当r＜0时, sign(r)＜0。

可以采用贪婪算法^[18]或者迭代阈值算法等求解(5)式, 下面介绍迭代阈值法求解的具体算法流程。

初始化: 赋初始值r₀, 外部循环最大迭代次数L_out, 内部循环最大迭代次数L_in, 残差为ε

迭代开始

外循环: 判断‖G－WFr^k‖₂>ε∩k≤L_out

内循环: j=0 to L_in－1, r^j+1=T_λk[r^j+1/α(WF)^T(G－WFr^j), λ/α]

执行: k=k+1;r^k=r^L

输出: $\tilde r$

上述算法流程中, T_λ为阈值函数; α, λ为参数; α必须大于WF的最大奇异值的平方。

1.2 基于反射系数反演的频带双向展宽

反演得到反射系数$\tilde r$后, 采用傅里叶变换算子F将$\tilde r$变为$\mathit{\boldsymbol{\tilde R}}$, 为下一步做准备:

$ \mathit{\boldsymbol{\tilde R}} = \mathit{\boldsymbol{F}}\mathit{\tilde r} $

(6)

$ \mathit{\boldsymbol{\tilde R}}$是采用有限带宽地震反演得到的全带宽反射系数频谱, 用于重构地震资料、实现频带双向展宽。频带双向展宽的过程如下。

1) 选取基于反射系数的频谱$\mathit{\boldsymbol{\tilde R}}$双向重构原数据频谱G的两个起始点N₁, N₂, 其中N₁ < N₂, N₁对应于向低频拓展的起始频率f₁, N₂对应于向高频拓展的起始频率f₂。

2) 提取N₁点之前$\mathit{\boldsymbol{\tilde R}}$的频谱值、N₁, N₂之间G的频谱值以及N₂之后$\mathit{\boldsymbol{\tilde R}}$的频谱值, 将上述频谱值组合为频带双向展宽补偿后的频谱$\mathit{\boldsymbol{\tilde G}}$, 即:

$ \mathit{\boldsymbol{\tilde G}} = {\mathit{\boldsymbol{F}}^{ - 1}}(\mathit{\boldsymbol{L}}(\mathit{\boldsymbol{\tilde R}}) + \mathit{\boldsymbol{M}}(\mathit{\boldsymbol{G}}) + \mathit{\boldsymbol{H}}(\mathit{\boldsymbol{\tilde R}})) $

(7)

式中; L, M和H分别为提取低频、中频和高频频谱值的算子; F^-1为反傅里叶变换算子。

3) 对$\mathit{\boldsymbol{\tilde G}}$进行反傅里叶变换, 即可得到频带双向展宽补偿后的地震资料。

2 数值试算 2.1 模型试算

为了验证利用反射系数反演对地震资料双向拓频方法的正确性、稳定性和抗噪性, 分别构建了无噪和含噪的单道反射系数模型, 在模型中直接设定了相应的反射系数值, 子波采用主频为30Hz的雷克子波, 子波与反射系数模型褶积得到合成地震记录, 该合成地震记录采样点数为501, 采样间隔为2ms, 对于该地震合成记录采用本文方法反演反射系数, 将反射系数与子波褶积得到反演地震记录, 上述模型及反演结果分别如图 1和图 2所示。对比可知, 反演反射系数与原始反射系数在时间和极性上完全对应, 幅值方面略有差异, 通过能量匹配可以消除幅值差异的影响, 从而得到更好的反演结果。含噪模型测试结果验证了本文方法的抗噪性。基于反射系数反演对原始地震资料进行双向拓频频谱分析, 得到了无噪情况下一系列频谱对比结果, 如图 3所示。图 3a是反演得到的反射系数频谱, 即通过有限带宽地震资料反演得到的全带宽反射系数频谱; 图 3b是原始地震记录频谱; 图 3c是原始地震记录频谱与反演反射系数频谱的对比结果; 图 3d是对原始地震记录频谱进行低频补偿后得到的频谱, 图 3e是对原始地震记录频谱进行高频补偿后得到的频谱, 图 3f是对原始地震记录频谱进行双向拓频补偿后得到的频谱, 对比拓频前、后的地震记录和对应的频谱可知, 本文方法能很好地对地震资料进行双向频带拓展, 实现地震分辨率的提升, 模型试验结果验证了本文方法的正确性、稳定性和抗噪性。

2.2 实际资料试算

在压缩感知框架下基于反射系数反演对地震资料进行双向拓频, 从本质上看是纯数据驱动的, 即双向拓频的补偿能量出自数据本身, 其中的稀疏反演步骤旨在得到更为准确的反射系数, 实现地震资料的重构, 从而获得具有更高分辨率特性的地震资料。

利用江苏油田某地区的实际地震资料进行测试, 主要包括2个步骤: ①提取地震子波; ②利用压缩感知理论反演得到全带宽的反射系数, 并通过全带宽的反射系数对地震记录进行双向拓频。测试参数如下: λ=0.9;α值通过(WF)^TWF矩阵对应的最大特征值来确定, 不同数据对应的观测矩阵不同; 迭代次数为100次, 试验证明当迭代次数足够大时, 其迭代出的反射系数趋于稳定, 不存在异常变化; λ用于调节反演得到的反射系数稀疏度。

图 4a和图 4b分别为江苏油田某线双向拓频前的原始地震剖面和采用本文方法双向拓频后的地震剖面。整体而言, 图 4b的地震剖面质量得到提高, 同相轴明显变窄, 层位关系更加清晰, 原本不清晰的同相轴能量也得到了增强, 分辨率明显提高, 反射波组关系有所改善, 能够更好地反映地下层组间的接触关系。

图 5为双向拓频前、后的频谱对比结果, 其中图 5a和图 5b分别为原始地震资料第116道的频谱和拓频后的频谱, 图 5c和图 5d分别为原始地震资料第363道的频谱和拓频后的频谱。以最高能量的1/8界定有效频带(图 5中横线所对应的频带范围), 可以看出, 低频拓展超过6Hz, 高频拓展约20Hz, 有效频带展宽超过26Hz, 在低频端和高频端两个方向上拓宽了地震记录优势频率的带宽。拓展补偿后的地震资料的能量数量级与原始地震资料一致, 本文方法可以有效减小地震记录同相轴的旁瓣幅值, 取得了良好的补偿效果, 验证了本文方法的有效性。

图 6为江苏油田某二维线原始地震剖面与双向拓频数据分频(1~7Hz频率段)地震剖面, 利用本文方法进行双向拓频处理后, 图 6b所示的地震剖面中低频信息丰富有效。图 7为江苏油田某二维线原始地震剖面与双向拓频数据分频(35~60Hz频率段)地震剖面, 对比分析可以看出, 利用本文拓频方法处理后, 图 7b所示的地震剖面中高频信息能量得到增强, 特征清晰, 进一步验证了本文方法的有效性。

3 结束语

本文研究了在压缩感知框架下基于反射系数反演对地震资料进行双向拓频的方法。模型和实际资料测试结果验证了方法的有效性, 低频和高频拓宽结果均真实可靠, 双向拓频不仅增强了地震弱信号的能量强度, 而且在没有改变原始地震剖面特征的同时, 使地震资料有效频带拓宽了20Hz以上, 提高了地震资料识别薄储层的能力, 改善了地震资料品质, 为后续地震解释提供了高分辨率的地震资料, 有利于识别地下地质构造和寻找油气。该方法适合在相关探区推广应用。在现有双向拓频方法的基础上, 后续将在时频域进行有针对性的均衡补偿, 以期取得更好的应用效果。