2018,
Vol. 42
2. 中国石化石油工程地球物理有限公司胜利分公司, 山东东营 257088
2. Shengli Branch, Geophysical Corporation, SINOPEC, Dongying 257088, China
随着油气勘探开发的不断深入, 地震勘探目标逐渐转向复杂构造[1]、地层[2]和岩性圈闭油气藏[3], 其后期数据处理技术[4-5]对地震勘探数据的规则性及完整性提出了更高的要求。实际数据采集时, 高精度的仪器设备保障了时间方向的无假频高密集采样, 但空间上受野外环境(某些难以逾越的障碍物和禁采区)、复杂地质条件以及废道、废炮剔除等多种复杂因素影响, 导致地震数据缺失。对此, 众多学者提出并发展了基于波动方程、预测滤波理论、各种数学变换域、矩阵降秩理论以及近年来较为热门的压缩感知理论的重构方法等。霍志周等[6]仅对前3种方法的发展做了较详细总结, 且没有系统地给出这些方法使用条件。如何在众多的方法之中进行最优的选择是在实际进行数据重构时需要面临的关键问题之一。笔者对地震道缺失在频率波数域的表现特征进行深入分析的基础上, 对现有的重构方法进行归类和分析, 总结各种方法的适用条件及其优缺点, 给出在不同条件下重构方法的选择策略, 并就不同的重构方法进行模型实验和对比分析。
1 地震数据缺失在频率波数域的表现特征实质上, 缺失地震数据属于不规则地震数据。高建军等[7]根据不规则地震数据的特点, 将其分为4类:
(1) 稀疏地震道数据。空间规则采样, 道间距较大, 地震道数据完整。
(2) 非均匀道缺失地震道数据。空间规则采样, 地震道随机缺失。
(3) 不规则地震道数据。空间不规则采样, 道间距随机, 地震道数完整。
(4) 混合不规则地震道数据。空间不规则采样, 且地震道不足。
一般将第一类称之为规则采样不足, 第二类、第三类和第四类称之为非规则采样不足。地震数据的规则或非规则采样不足都会对地震数据的后续处理产生不利的影响。
图 1展示了规则采样不足地震道及其频谱特征。图 1(a)为一完整的空间规则采样的多层地震数据模型, 道间距为25 m, 采样频率为2 ms, 图 1(b)为空间稀疏规则采样的地震数据, 道间距为50 m, 采样频率也是2 ms, 图 1(c)和图 1(d)分别为其对应的归一化频谱(f-k)图。如图 1(d)所示, 地震数据经稀疏采样后, 频谱图中出现了部分回折的空间假频。实际上地震数据的稀疏采样也可以看成是数据的规则采样不足。
|
图 1 空间规则稀疏采样及f-k谱 Fig.1 Regular sparse sampling and corresponding f-k spectrums |
图 2展示了非规则采样不足情况下随机缺失不同程度的地震数据的频谱图对比。图 2(a)为完整的多层地震数据模型, 道间距为25 m, 图 2(b)、(c)分别为图 2(a)随机缺失20%和50%的地震数据模型, 图 2(d)~(f)分别对应前述3种模型的归一化频谱(f-k)图。如图 2(d)所示, 地震数据在波数域是稀疏的, 其能量主要集中在少量波数成分上。由图 2(e)可以看出, 采样随机缺失20%的地震数据时, 频谱的能量向各个波数进行分散, 出现空间假频, 随着缺失数据的增加到50%, 如图 2(f)所示, 频谱能量更加分散, 空间假频现象更加严重。
|
图 2 缺失地震数据及其f-k谱 Fig.2 Missing seismic data and corresponding f-k spectrums |
通过对比分析得出如下结论:地震数据的规则或非规则采样不足本质上均会造成其频谱能量在波数域分散, 影响后续的处理结果, 尤其是偏移成像时, 如果偏移剖面上存在空间假频及频散现象, 会使横向分辨率变差, 导致解释人员无法正确解释地下构造。
2 地震数据重构方法 2.1 规则采样不足地震道重构方法规则采样不足地震道重构主要是针对空间规则稀疏采样的情况, 研究的重点是地震道抗假频插值重构, 基本方法有基于倾角扫描的重构方法、基于预测滤波理论的重构方法以及基于波动方程的重构方法。
2.1.1 基于倾角扫描的重构方法倾角扫描是最早使用的地震插值技术, Pieprzak等[8]提出在地震数据上开时窗, 在每段时窗中沿某个倾角进行等角度扫描, 找出相似或一致的特征, 最后根据这些特征恢复缺失地震道。假设每段时空窗的地震记录是由有效波(定义在3~5道连续的地震道内存在线性同相轴)与不相关的干扰波组成, 第i个时窗的地震数据在x-ω域的表现形式为
| $ {\rm{[}}{\mathit{\boldsymbol{G}}_\mathit{i}}\left( \mathit{\omega } \right)\left] {\rm{ = }} \right[\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{H}}_{\mathit{i}{\rm{, }}\;\mathit{j}}}\left( \mathit{\omega } \right)\left] {\rm{\cdot}} \right[{\mathit{\boldsymbol{S}}_\mathit{j}}\left( \mathit{\omega } \right)\left] {\rm{ + }} \right[{\mathit{\boldsymbol{R}}_\mathit{i}}\left( \mathit{\omega } \right){\rm{]}}{\rm{.}} $ | (1) |
式中, Gi(w)为列向量, 地震记录在第i个时窗的频谱; SHi, j(w)为线性时移矩阵(由对角元素为ejiωdt的子矩阵组成, 其初始值为扫描的角度信息); Sj(w)为第j个倾角的有效波的频谱, 其维数为倾角j与频率数的乘积; Ri(w)为不满足线性同相轴假设的干扰波的频谱。求解由各时窗频谱组成的线性方程组就可以得到重构的地震数据Sj(w)。
局部倾角的自动拾取需要有一定的先验信息(最大、最小反射同相轴倾角等), 其微小误差都会对扫描倾角之外的反射同相轴造成不利的影响，因此该方法一般仅用于单道数据。
2.1.2 基于预测滤波理论的重构方法预测误差滤波方法是重构含有空间假频的规则采样数据的常用方法[9]。对于N道等间距地震记录, 含L组同相轴(每组同相轴的视速度不同), 第i组同相轴表示为Ai(t)(i=1, 2, …, L), pi为相邻两道间的时移, 则地震道gk在频率域的描述为
| $ {\mathit{\boldsymbol{G}}_\mathit{k}}\left( \mathit{f} \right){\rm{ = }}\sum\limits_{\mathit{i} = 1}^\mathit{L} {{\mathit{\boldsymbol{A}}_\mathit{i}}} \left( \mathit{f} \right)\mathit{\boldsymbol{Z}}_\mathit{i}^{\mathit{k}{\rm{ - 1}}}\left( \mathit{f} \right), \;\;\mathit{k}{\rm{ = 1, 2, }} \cdots {\rm{, }}\mathit{N}{\rm{.}} $ | (2) |
式中, Ai(f)是第i组同相轴Ai(t)的傅氏变换; Zi(f)=ei2πifp, Ai(f)在频率f处的相移。向前、向后一步预测算子Pi(f)可以应用最小二乘原理由下式求出：
| $ \left\{ \begin{array}{l} {\mathit{\boldsymbol{G}}_\mathit{k}}\left( \mathit{f} \right) \approx \sum\limits_{\mathit{i} = 1}^\mathit{L} {{\mathit{\boldsymbol{P}}_\mathit{i}}} \left( \mathit{f} \right){\mathit{\boldsymbol{G}}_{\mathit{k}{\rm{ - }}\mathit{i}}}\left( \mathit{f} \right){\rm{, }}\;\;\;\mathit{k}{\rm{ = }}\mathit{L}{\rm{ + 1, }} \cdots {\rm{, }}\mathit{N}{\rm{.}}\\ {\mathit{\boldsymbol{G}}_\mathit{k}}\left( \mathit{f} \right) \approx \sum\limits_{\mathit{i} = 1}^\mathit{L} {\mathit{\boldsymbol{P}}_\mathit{i}^*} \left( \mathit{f} \right){\mathit{\boldsymbol{G}}_{\mathit{k}{\rm{ + }}\mathit{i}}}\left( \mathit{f} \right){\rm{, }}\;\;\;\mathit{k}{\rm{ = 1}}, \cdots {\rm{, }}\mathit{N}{\rm{ - 1}}{\rm{.}} \end{array} \right. $ | (3) |
式中, *表示复共轭。
如果对Gk(f)做一阶内插, 得到由偶数道分量和奇数道分量组成的Gk(f)(k=1, 2, …, 2N-1), 以及内插后在频率f处一步预测算子P′i(f), 其表达式与式(2)和式(3)类似。偶数道分量则可以由奇数道分量和P′i(f)构成的矩阵预测得到。
为了提高计算效率, Sacchi[10]又发展了f-k域道内插值算法, 对应的f-k域非假频区能够得到高保真的内插结果, 但在假频区会产生错误的内插结果。李国发[11]结合f-k域与f-x域道内插的优缺点, 形成速度快、精度高、抗空间假频的内插技术。为适应不均匀采样数据的处理, Naghizadeh等[12]和Liu等[13]又结合傅里叶变换使其应用范围拓展到非规则地震数据重建。
基于预测滤波理论的重构方法能够较好地处理含有空间假频的规则采样数据, 在处理非规则采样数据的处理时, 需要结合稀疏反演理论等方法。
2.1.3 基于波动方程的重构方法基于波动方程的地震数据重构方法以积分算子为基础, 应用求和方式代替积分方程进行计算, Bolondi[14]指出此方法能够恢复并提高常规处理中损失的空间分辨率。该方法主要包括DMO(倾角时差校正)方法、偏移距域延拓法和炮域延拓法, 其中以DMO方法的研究较多, 以此为例简要介绍其重构原理。对于采样不足的N道地震数据记录dj, 在波数域由积分算子和正演模型表示为
| $ {\mathit{\boldsymbol{d}}_\mathit{j}}\left( \mathit{k} \right){\rm{ = }}\sum\limits_\mathit{n} {\mathit{\boldsymbol{D}}_\mathit{j}^ + } {\rm{(}}\mathit{k}{\rm{ - }}\mathit{n}{\mathit{k}_{\rm{s}}}{\rm{)}}\mathit{\boldsymbol{m}}{\rm{(}}\mathit{k}{\rm{ - }}\mathit{n}{\mathit{k}_{\rm{s}}}{\rm{)}}{\rm{.}} $ | (4) |
式中, Dj+为逆DMO(DMO-1)算子; m为正演模型; ks为采样波数(ks=2π/x, x为采样间隔)。由式(4)可知将DMO作为插值算子, 对dj做DMO处理, 得到正演模型, 再做正常时差校正(NMO)、DMO和叠加(stacking)处理, 通过局部的时间倾角和共轭梯度迭代等方法实现数据规则化重建, Ronen[15]应用该方法实现数据规则化重建。为了加快运算速度, Canning等[16]将地震数据的时间坐标轴变换到对数尺度, 在f-k域逐步进行DMO和DMO-1, 但对数拉伸后数据量增加, 受计算机内存限制计算效率低。管路平等[17]应用DMO积分算法分步校正方位角, 解决了含假频的叠前地震数据规则化问题, 但受积分孔径影响较大。
基于波动方程的地震数据重建方法可以充分利用地下介质的相关速度分布信息(叠加速度、偏移速度、均方根速度等), 且由于积分算子可以处理各种几何关系的采样, 在地震数据处理中可以灵活使用, 提高重构精度。但是, 当地下信息不准确或未知时, 重建结果会受到严重影响, 如果地震数据为稀疏采样, 重建后还会出现空间假频, 而且该方法计算量大、耗时长, 在实际生产中难以推广应用。
2.2 非规则采样不足地震道重构方法非规则采样不足地震数据重构的核心是抗假频数据规则化和插值重建。目前发展的方法主要有基于数据变换的重构方法、基于矩阵降秩理论的重构方法以及基于压缩感知理论的重构方法。
2.2.1 基于数学变换的重构方法基于数学变换的地震数据重构方法主要包括Fourier变换、Radon变换、Wavelet变换、curvelet变换等, 其中Fourier变换是其他数学变换方法的基础。
对于一个完整的地震数据g可分为存在和缺失两部分, g可以看做是模型空间, 将采样不足地震道的插值问题描述为最小平方问题, 即
| $ \begin{array}{l} {\rm{0}} \approx {\left\| {\mathit{\boldsymbol{d}}{\rm{ - }}\mathit{\boldsymbol{Mg}}} \right\|_{\rm{2}}}{\rm{, }}\\ {\rm{0}} \approx \mathit{\varepsilon }{\left\| {\mathit{\boldsymbol{Cg}}} \right\|_{\rm{2}}}{\rm{.}} \end{array} $ | (5) |
式中, d为插值前地震数据; M为mask矩阵, 即对应缺失数据的位置元素为0, 其他为1;C为稀疏变换算子, 约束插值后的地震数据满足一定的条件, ≈表示两端误差最小。式(5)的理论解为
| $ \left\langle {\mathit{\boldsymbol{g}}} \right\rangle {\rm{ = }}{\mathit{\boldsymbol{C}}^{{\rm{ - 1}}}}{{\rm{(}}{\mathit{\boldsymbol{C}}^{{\rm{ - }}\mathit{T}}}{\mathit{\boldsymbol{M}}^{\rm{T}}}\mathit{\boldsymbol{M}}{\mathit{\boldsymbol{C}}^{{\rm{ - 1}}}}{\rm{ + }}\mathit{\varepsilon }\mathit{\boldsymbol{I}}{\rm{)}}^{{\rm{ - 1}}}}{\mathit{\boldsymbol{C}}^{{\rm{ - }}\mathit{T}}}{\mathit{\boldsymbol{M}}^{\rm{T}}}\mathit{\boldsymbol{d}}{\rm{.}} $ | (6) |
Fourier变换对地球物理领域的频谱分析、f-k滤波去噪和数据重构发挥了重要作用。Duijindam[18]提出了直接应用傅里叶理论和参数反演重构任意不规则采样数据, 在地震数据有限带宽的条件下, 无需经过任何地球物理假设; 为了有效去除空间假频, 刘喜武等[19]又提出离散傅里叶变换(DFT)加权范数规则化算法, 将带限地震数据重建归结为最小二乘反演问题, 将功率谱支持和谱的形状结合到重建过程中; 霍志周等[20]则将Fourier变换与贝叶斯参数反演方法结合, 使用预条件共轭梯度法实现抗假频重建, 加快收敛速度, 提高计算效率。
Radon变换是一种积分变换, 根据积分路径分为线性Radon变换、双曲Radon变换、抛物Radon变换。线性Radon变换一般用于垂直地震剖面(VSP)、直达波等具有线性特征的地震数据同相轴的波场分离; 双曲Radon变换的积分路径与CMP道集内地震数据的同相轴的双曲特征更为相近[21], 能够提高重构的分辨率, 但其在时空域需要计算一个大矩阵的逆, 计算效率低; 抛物Radon变换算子具有时不变性, 可以进行傅里叶变换并对每个频率进行单独求解, 计算效率较线性和双曲Radon变换有很大的提升, 但受固有分辨率的限制, 也不能得到理想的处理效果。后来, 王维红等[22]提出加权抛物线Radon变换重构方法, 提高重构地震数据分辨率的同时也提高了计算效率和精度, 但权系数受地震剖面的品质影响较大; 薛亚茹等[23]将高阶3D Radon变换与正交多项式面拟合结合, 进一步完善3D Radon变换, 提高了分辨率。
Wavelet变换和curvelet变换都属于多尺度变换, 其目的是获得时间与频率的相互关系。Wavelet变换可以对非平稳信号进行描述和压缩, 郭念民等[24]实现了非规则采样不足地震数据的插值重构, 但对含有弯曲同相轴的地震数据需要做分频处理。curvelet变换是在脊波变换理论的基础上提出的, 具有良好的各向异性, 能够高效描述地震数据的边缘特征。Naghizedh等[25]利用其稀疏性, 实现地震数据重建和弯曲同相轴的重构和恢复, 为了加快迭代收敛速度, 提高对含噪声的地震数据插值效果, Wang等[26]采用改进的基于curvelet变换的加权凸集投影方法, 使其更适用于实际地震数据的处理。
基于数学变换的重构方法利用地震数据在各数字域中的稀疏性实现规则或非规则地震数据重构, 不需要各种地质、地球物理等先验信息, 但在处理非规则地震数据时抗假频性能较差。
2.2.2 基于矩阵降秩理论的重构方法基于降秩理论的地震数据重构方法主要有Cadzow滤波法、奇异谱分析法(SSA)和多道奇异谱分析法(MSSA)等。
Cadzow滤波法[27]最早由Cadzow在1988年提出, 基本原理与MSSA降秩法类似, 首先将地震数据变换到频率域, 对每一个频率切片的数据构建(块)Hankel矩阵, 应用截断SVD(奇异值分解)进行降秩运算, 将降秩后的矩阵沿反对角线取平均, 完成数据重构。传统的截断SVD算法计算量大, 难以实现降秩处理, Oropeza等[28]用随机SVD算法取代截断SVD, 对低维数据重构有一定效果; 魏小强等[29]结合POCS(凸集投影法)和随机SVD算法对MSSA重构方法进行改进, 实现非规则地震数据的插值。
降秩理论重构方法与地震数据重构物理意义贴合密切, 但计算量随数据矩阵阶数的增加而增大, 缺少抗假频机制, 且无法对规则采样的稀疏地震数据进行加密处理。因此, 该方法在含假频的非规则地震数据重构和规则地震数据加密中的应用还有待进一步的研究。
2.2.3 基于压缩感知理论的重构方法压缩感知(compressive sensing, CS)是由Donoho等[30]提出的寻找线性欠定系统稀疏解的技术, 打破了传统的Nyquist采样理念。该理论认为, 如果信号是稀疏的, 或者在某个变换域可以被稀疏表示, 可以通过设计一个与稀疏变换基不相关的测量矩阵来观测该信号, 从而将高维信号投影到低维空间上, 得到少量观测值, 利用稀疏促进算法求解最优化问题实现数据重构。Herrmann[31]等首先将该理论应用到地震数据的恢复中, 提出了CRSI(curvelet-based recovery by sparsity-promoting inversion)方法, 主要包括3部分:稀疏表示、压缩采样和稀疏促进算法设计。对于一个完整信号f∈RN, 设f在Ψ变换域的稀疏表示为
| $ \mathit{\boldsymbol{f}}{\rm{ = }}\mathit{\boldsymbol{ \boldsymbol{\varPsi} u}}{\rm{.}} $ | (7) |
式中, u为变换域系数向量。测量完整信号f到不完整信号y∈RM(M≪N)的映射, 矩阵形式表示为
| $ \mathit{\boldsymbol{y}}{\rm{ = }}\mathit{\boldsymbol{ \boldsymbol{\varPhi} f}}{\rm{ }}{\rm{.}} $ | (8) |
式中, Φ∈RM×N(M≪N)为采样矩阵。将式(7)带入式(8)得
| $ \mathit{\boldsymbol{y}}{\rm{ = }}\mathit{\boldsymbol{ \boldsymbol{\varPhi} \boldsymbol{\varPsi} u}}{\rm{ = }}\mathit{\boldsymbol{ \boldsymbol{\varTheta} u}}{\rm{.}} $ | (9) |
式中, Θ∈RM×N为测量矩阵; y为测量得到的高维原始信号在低维空间上的投影, 因此式(9)为欠定方程, 压缩感知理论通过设计合理的稀疏促进算法求解该方程的唯一解。Ma[32]结合Bregman和ATV(anisotropic total variation)迭代算法, 构造联合迭代算法框架, 加快了迭代初期的收敛速度, 同时避免了迭代后期噪声干扰的影响。
对于稀疏的地震信号, 压缩感知理论作为一种新的数据重构方式, 需要的先验信息最少, 也不需要像抛物线Radon变换、DMO等方法一样预先去除面波或者预知地下介质的相关信息(倾斜层或速度模型)等, 如果结合数学变换对地震信号进行稀疏表示, 则比传统的重构方法更加有效。
3 地震数据规则化重构方法选择策略应用地震数据规则化重构方法时需考虑具体的实现条件和重构需求, 如地震数据空间采样的规则性、稀疏程度以及先验信息的已知程度等, 具体的重构策略如图 3所示。
|
图 3 地震数据重构选择策略 Fig.3 Optimal strategy of reconstruction method |
对于不含空间假频的地震数据, 分为已知地层倾角信息、已知速度信息以及未知先验信息等3种情况进行讨论:
(1) 如果已知地层倾角信息, 可以选择基于倾角扫描的重构方法, 但使用的倾角信息必须要准确, 否则会影响扫描倾角之外的反射同相轴, 因此该方法一般只用于单道数据处理。
(2) 如果已知的是速度信息, 可以选择基于波动方程的重构方法, 在f-k域逐步进行DMO和DMO-1, 利用DMO积分算法分步校正方位角, 通过局部的时间倾角和共轭梯度迭代等方法实现数据规则化重建。该方法计算量大、耗时长, 如果地下介质的相关速度分布信息不准确时, 重建结果会受到严重影响, 因此在实际生产中推广较难。
(3) 无论地层倾角和速度信息已知与否, 均可以采用预测误差滤波理论、基于各种数学变换的重构方法。预测误差滤波理论有f-x域、τ-x域、f-k域道内插值; 基于各种数学变换的方法, 如变权系数的抛物线Radon变换、小波变换等。
对于含有空间假频的地震数据, 可以考虑选择预测误差滤波方法和结合稀疏反演的各种数学变换方法:
(1) 预测误差滤波算法。f-x域、τ-x域、f-k域等道内插值方法, 重构效率高, 适用于地震道加密处理。
(2) 结合稀疏反演的各种数学变换方法。数学变换方法可以对规则数地震数据进行快速高精度重构, 但是抗假频能量差, 因此需要结合稀疏反演方法进行抗假频处理。例如, 利用curvelet变换的稀疏性, 通过L2范数约束求解正则化反问题, 利用粗尺度无假频的曲波系数约束细尺度有假频的曲波系数, 实现地震数据重构; 应用快速广义Fourier变换(FGFT)识别非平稳信号在每段时间频率空间波数的变化, 结合最小二乘拟合法恢复最优FGFT系数, 实现地震数据重构。
3.2 非规则采样数据重构策略对于非规则采样地震数据, 如果不含空间假频, 可以选择基于各种数学变换的方法或基于矩阵降秩理论等即可得到准确的重构数据, 无需任何先验信息。
(1) 基于各种数学变换的方法。应用加权抛物线Radon变换、非离散小波变换及基于曲波框架和类傅里叶变换的鲁棒性迭代稀疏约束插值算法等。此类方法处理不含假频的地震数据, 计算速度快, 效率高, 且对输入数据要求较少, 被广泛应用于生产中。
(2) 基于矩阵降秩理论的方法。MSSA方法结合POCS(凸集投影法)和随机奇异值分解法等, 对不含噪声的地震数据处理效果较好。该方法提出时只能用于非规则地震数据处理, 且受噪声影响较大, 如果用于实际资料处理还需要进行深入研究。
对于含空间的假频的非规则采样数据, 处理难度大大增加, 也是实际生产中存在的不可避免的问题, 从先验信息已知和未知两方面进行讨论:
(1) 如果已知速度等先验信息, 可以选择基于波动方程的重构方法, 但需要结合数学变换方法, 例如, 将偏移距域延拓与seilet变换联合, 应用凸集映射或Bregman迭代算法精确描述复杂波场信息解, 解决地震数据插值中的无约束问题。该方法的缺点是对地下介质的相关速度分布信息较依赖。
(2) 如果先验信息未知, 则可以选择基于预测滤波理论结合稀疏反演理论的重构方法或数学变换结合压缩感知理论的重构方法。因为预测滤波理论方法适用于规则含假频的地震数据重构, 可以结合傅里叶变换将其扩展到非规则地震数据重构; 基于数学变换的重构方法适用于非规则无假频的地震数据重构, 扩展到含假频的情况需要结合抗假频的压缩感知理论。例如, 应用多步自回归预测滤波方法, 结合傅里叶变换构造和自适应f-x域插值算法, 通过求解一个全局最小二乘约束问题获得光滑非平稳预测滤波误差系数, 直观地选择规则化系数; 在曲波域可以联合Bregman算法与ATV最小化算法, 在保留数据特征和边缘信息的同时降低噪声和假频干扰, 还可以结合最速下降和Bregman两种迭代算法, 加快迭代初期的收敛速度, 实现地震数据的重构。
4 模型测试实际生产中含假频的地震数据对后期的数据处理影响比较大, 且速度等先验信息一般也较难获取, 针对先验信息未知的含假频的规则和非规则采样不足地震数据, 对图 1和图 2所示的缺失地震数据进行模型测试, 比较不同重构策略的处理效果和适用性。
4.1 规则采样不足对于含假频的规则采样不足的地震数据可以采用预测误差滤波方法和结合稀疏反演的各种数学变换方法。这里的预测滤波方法选择不受空间假频影响的f-x域道内插值方法, 数学变换选择目前研究应用较多且重构精度较高的曲波变换。
图 4为两种方法重构的结果对比, 图 4(a)、(b)分别为f-x域道内插值方法和曲波变换+稀疏反演插值方法的地震数据，图 4(c)、(d)分别为两种重构方法分别对应的f-k频谱。由图 4(a)和图 4(b)可以看出, 两种方法都能够对含假频的规则采样不足的地震数据进行较好的插值, f-x域道内插值方法重构出的地震数据产生的噪声较多, 实际应用时还需要进行相应的去噪处理; 对比图 4(c)和图 4(d)可以看出, f-x域道内插值方法可以较好地集中低波数的频率, 而对高波数的频率较分散, 曲波变换由于结合了抗假频的稀疏反演方法, 能够更好地克服假频现象, 将波数域的高频和低频能量分别集中在少数波数成份上, 更好地实现完整地震数据重构。
|
图 4 规则采样不足地震数据重构效果 Fig.4 Reconstructed image for regular model |
表 1对比了这两种方法进行重构时使用的时间和重构结果的峰值信噪比, 可以看出f-x域道内插值方法重构的速度更快一些, 但是curvelet变换+稀疏反演的重构方法的峰值信噪比要更高一些。
|
表 1 f-x域道内插值和曲波变换+稀疏反演插值重构结果对比 Table 1 Interpolation of f-x and curvelet+sparse inversion for regular model |
对比可知, 在对规则稀疏地震数据进行插值时, 需要根据实际需要选取不同的方法, 即如果数据量比较大, 对时间要求较紧时应尽量选择f-x域道内插值方法; 如果对精度要求较高且时间较宽裕, 可以选择curvelet变换+稀疏反演的重构方法。
4.2 非规则采样不足对于含假频的非规则采样不足的地震数据可以采用基于预测滤波理论结合稀疏反演理论和数学变换结合压缩感知理论的方法, 这里的预测滤波方法和数学变换同样选取f-x域道内插值和curvelet变换。稀疏反演理论和压缩感知理论都具有抗假频的性质, 因此还同时对比了数学变换结合稀疏反演理论的重构方法。
图 5为f-x域道内插值+稀疏反演、曲波变换插值+稀疏反演和曲波变换+压缩感知3种方法对非规则采样不足(20%)的采样地震数据进行重构的结果。对比图 5(a)~(c)可以看出, f-x域道内插值+稀疏反演的重构方法尽管在重构过程中产生的噪声较多, 但对缺失的地震道能够基本恢复; curvelet变换+稀疏反演的重构方法虽然在重构过程中产生的噪声较少, 但对缺失的地震道恢复的效果较差, 从图 5(b)也可以看出对缺失道中较弱的信号重构效果较差; curvelet变换+压缩感知的重构方法在重构过程中产生的噪声最少, 缺失地震道中较弱的信号也能够较好地恢复出来。对比频谱图 5(d)~(f)可以看出, f-x域道内插值重构方法结合稀疏反演可以基本恢复低波数的频率成分, 而对高波数的频率恢复的较差; curvelet变换+稀疏反演/压缩感知两种重构方法都能够使低波数和高波数的能量得到很好的集中, 由于压缩感知理论能够用较少的数据即可高精度重构恢复原始完整信号, 因此其抗假频性能更好, 对波数域的能量集中得更好。也可以看出, 压缩感知理论作为一种更加有效的重构方法, 是未来地震数据高精度重构的一个发展方向。
|
图 5 非规则采样不足(20%)地震数据重构效果 Fig.5 Reconstructed image for irregular model with 20% missing data |
表 2对比了3种方法进行重构时使用的时间和峰值信噪比。从表中可以看出, f-x域道内插值+稀疏反演的重构方法重构的速度更快一些, 但重构精度较低; curvelet变换+压缩感知的重构方法的峰值信噪比虽然较高, 但其重构所用时间较长; curvelet变换+稀疏反演的重构方法耗时也较长, 重构效果较差。基于curvelet变换的重构方法耗时长是由于需要将地震数据变换到curvelet域进行插值和规则化处理然后再进行反变换, 因此该方法的发展方向之一是减少迭代次数, 提高迭代效率。
|
|
表 2 f-x域道内插值+稀疏反演、曲波变换+稀疏反演和曲波变换+压缩感知插值重构结果对比 Table 2 Interpolation of f-x +sparse inversion, curvelet+sparse inversion and curvelet+CS for regular model |
通过对上述3种方法的分析可知, 非规则采样不足的地震数据重构可以选择预测误差方法(f-x域道内插值)+稀疏反演和数学变换(curvelet变换)+压缩感知, 如果对时间要求较紧且对结果的精度要求不高时可以选择预测误差方法(f-x域道内插值)+稀疏反演重构方法, 如果对精度要求较高且时间较宽裕, 则可以选择数学变换(curvelet变换)+压缩感知重构方法。
5 结论(1) 地震数据在波数域是稀疏的, 其能量主要集中于少数波数成分, 地震数据的规则或非规则采样不足均会造成频谱能量在波数域分散, 导致空间假频出现。对于非规则采样不足的地震数据随缺失地震数据的增加, 频谱能量逐渐向各个波数分散, 加重空间假频现象。因此, 地震数据重构方法的重点是抗假频重建。
(2) 地震数据规则化重构策略是根据重构条件和使用需求对不同的重构方法进行合理选择。对规则采样不足的地震数据进行重构, 如果不含空间假频可以采用预测误差滤波理论和基于各种数学变换的重构方法; 如果含空间假频, 且对计算效率要求比较高, 可以选择基于预测滤波理论的方法, 如果对精度要求比较高, 则可以选择基于数学变换与抗假频的稀疏反演组合的方法。对非规则采样不足的地震数据进行重构, 如果不含空间假频, 可以选择基于各种数学变换或基于矩阵降秩理论等方法即可得到准确的重构数据; 如果含空间假频, 且对计算效率要求比较高, 可以选择基于预测滤波理论与稀疏反演结合的方法; 如果对精度要求较高, 则可以选择数学变换与压缩感知理论进行结合。
| [1] |
王岩泉, 边伟华, 刘宝鸿, 等. 辽河盆地火成岩储层评价标准与有效储层物性下限[J]. 中国石油大学学报(自然科学版), 2016, 40(2): 13-22. WANG Yanquan, BIAN Weihua, LIU Baohong, et al. Evaluation criterion and cut-off value of igneous rock reservoirs in Liaohe Basin[J]. Journal of China University of Petroleum(Edition of Natural Science), 2016, 40(2): 13-22. |
| [2] |
马立民, 李志鹏, 林承焰, 等. 东营凹陷沙四下盐湖相沉积序列[J]. 中国石油大学学报(自然科学版), 2014, 38(6): 24-31. MA Limin, LI Zhipeng, LIN Chengyan, et al. Sedimentary sequences of salt-lake facies in Lower Es4 of Dongying Depression[J]. Journal of China University of Petroleum(Edition of Natural Science), 2014, 38(6): 24-31. |
| [3] |
杨田, 操应长, 王艳忠, 等. 东营凹陷沙三中亚段浊积岩低渗透储层有效性评价[J]. 中国石油大学学报(自然科学版), 2016, 40(4): 1-11. YANG Tian, CAO Yingchang, WANG Yanzhong. Effectiveness evaluation of low permeability turbidite reservoirs in the middle of the third member of Shahejie Formation in Dongying Sag[J]. Journal of China University of Petroleum(Edition of Natural Science), 2016, 40(4): 1-11. |
| [4] |
徐秀刚, 李振春, 叶月明, 等. 起伏地表条件下波动方程法共聚焦点成像技术[J]. 中国石油大学学报(自然科学版), 2010, 34(2): 38-42. XU Xiugang, LI Zhenchun, YE Yueming, et al. Wave equation common focus-point migration based on irregular surface[J]. Journal of China University of Petroleum(Edition of Natural Science), 2010, 34(2): 38-42. |
| [5] |
薛亚茹, 杨静, 钱步仁. 基于高阶稀疏Radon变换的预测多次波自适应相减方法[J]. 中国石油大学学报(自然科学版), 2015, 39(1): 43-49. XUE Yaru, YANG Jing, QIAN Buren. Multiples prediction and adaptive subtraction with high-order sparse Radon transform[J]. Journal of China University of Petroleum(Edition of Natural Science), 2015, 39(1): 43-49. |
| [6] |
霍志周, 熊登, 张剑锋. 地震数据重建方法综述[J]. 地球物理学进展, 2013, 28(4): 1749-1756. HUO Zhizhou, XIONG Deng, ZHANG Jianfeng. The overview of seismic data reconstruction methods[J]. Progress in Geophys, 2013, 28(4): 1749-1756. |
| [7] |
高建军, 陈小宏, 王芳芳, 等. 不规则地震道数据规则化重建方法研究[J]. 地球物理学进展, 2011, 26(3): 983-991. GAO Jianjun, CHEN Xiaohong, WANG Fangfang, et al. Study on regularized reconstruction of uneven seismic traces[J]. Progress in Geophys, 2011, 26(3): 983-991. |
| [8] |
PIEPRZAK A, MCCLEAN J. Trace interpolation of severely aliased events[C/OL]//1988 SEG Annual Meeting. Anaheim, California, October 30-November 3, 1988: SEG Technical Program Expanded Abstracts, c1988: 658-660. [2016-05-10]. http://library.seg.org/doi/pdfplus/10.1190/1.1892282.
|
| [9] |
SPITZ S. Seismic trace interpolation in the F-X domain[J]. Geophysics, 1991, 67: 794-890. |
| [10] |
SACCHI M, ULRYCH T J. High-resolution velocity gathers and offset space reconstruction[J]. Geophysics, 1995, 60(4): 1169-1177. DOI:10.1190/1.1443845 |
| [11] |
李国发. f-k域与f-x域联合实现道内插[J]. 石油地球物理勘探, 1995, 30(5): 693-701. LI Guofa. Joint trace interpolation in f-k and f-x domains[J]. Oil Geophysical Prospecting, 1995, 30(5): 693-701. |
| [12] |
NAGHIZADEH M, SACCHI M D. f-x adaptive seismic trace interpolation[J]. Geophysics, 2009, 74(1): V9-V16. DOI:10.1190/1.3008547 |
| [13] |
LIU Yang, SERGEY F. Seismic data interpolation beyond aliasing using regularized nonstationary auto regression[J]. Geophysics, 2011, 76(5): V69-V77. DOI:10.1190/geo2010-0231.1 |
| [14] |
BOLONDI G, LOINGER E, ROCCA F. Offset continuation of seismic sections[J]. Geophysical Prospection, 1982, 30: 813-828. DOI:10.1111/gpr.1982.30.issue-6 |
| [15] |
RONEN J. Wave-equation trace interpolation[J]. Geophysics, 1987, 52: 973-984. DOI:10.1190/1.1442366 |
| [16] |
CANNING A, GARDNER G. Regularizing 3-D data sets with DMO[J]. Geophysics, 1996, 61: 1108-1114. |
| [17] |
GUAN L P, TANG Y X, WANG H Z. Common-offset plane-wave prestack time migration and demigration[J]. Chinese Journal of Geophysics, 2009, 52(5): 1301-1309. |
| [18] |
DUIJNDAM AJW, SCHONEWILLE MA, HINDRIKS C O H. Reconstruction of band-limited signals, irregularly sampled along one spatial direction[J]. Geophysics, 1999, 64(2): 524-538. DOI:10.1190/1.1444559 |
| [19] |
刘喜武, 刘洪, 刘彬. 反假频非均匀地震数据重建方法[J]. 地球物理学报, 2004, 47(2): 299-305. LIU Xiwu, LIU Hong, LIU Bin. A study on algorithm for reconstruction of de-alias uneven seismic data[J]. Chinese Journal of Geophysics, 2004, 47(2): 299-305. |
| [20] |
霍志周, 熊登, 张剑锋. 预条件共轭梯度法在地震数据重建方法中的应用[J]. 地球物理学报, 2013, 56(4): 1321-1330. HUO Zhizhou, XIONG Deng, ZHANG Jianfeng. Application of the preconditioned conjugate gradient method to reconstruction of seismic data[J]. Chinese Journal of Geophysics, 2013, 56(4): 1321-1330. |
| [21] |
THORSON J R, CLAERBOUT J F. Velocity stack and slant stack stochastic inversion[J]. Geophysics, 1985, 50(12): 2727-2741. DOI:10.1190/1.1441893 |
| [22] |
王维红, 裴江云, 张剑锋. 加权抛物Radon变换叠前地震数据重建[J]. 地球物理学报, 2007, 50(3): 851-859. WANG Weihong, PEI Jiangyun, ZHANG Jianfeng. Prestack seismic data reconstruction using weighted parabolic transform[J]. Chinese Journal of Geophysics, 2007, 50(3): 851-859. |
| [23] |
薛亚茹, 陈健升, 钱步仁. 高阶3D Radon变换及其在数据重建应用[J]. 中国石油大学学报(自然科学版), 2016, 40(3): 69-76. XUE Yaru, CHEN Jiansheng, QIAN Buren. High order 3D Radon transform and its application in data reconstruction[J]. Journal of China University of Petroleum(Edition of Natural Science), 2016, 40(3): 69-76. |
| [24] |
郭念民, 李海山, 冯雪梅, 等. 非抽样离散小波变换叠前地震数据重建方法[J]. 石油地球物理勘探, 2014, 49(3): 508-516. GUO Nianmin, LI Haishan, FENG Xuemei, et al. Pre-stack seismic data reconstruction based on the undecimated wavelet transform[J]. Oil Geophysical Prospecting, 2014, 49(3): 508-516. |
| [25] |
NAGHIZADEH M, SACCHI M D. Beyond alias hierarchical scale curvelet interpolation of regularly and irregularly sampled seismic data[J]. Geophysics, 2010, 75: 189-202. |
| [26] |
WANG B, CHEN X, LI J, et al. An improved weighted projection onto convex sets method for seismic data interpolation and denoising[J]. IEEE Journal of Selected Topics in Applied Earth Observation and Remote Sensing, 2016, 9(1): 228-235. DOI:10.1109/JSTARS.2015.2496374 |
| [27] |
CADZOW J. Signal enhancement a composite property mapping algorithm[J]. IEEE Transactions on Acoustics, Speech and Signal Processing, 1988, 36(1): 49-62. DOI:10.1109/29.1488 |
| [28] |
OROPEZA V, SACCHI MD. Simultaneous seismic data denoising and reconstruction via multichannel singular spectrum analysis[J]. Geophysics, 2011, 76(3): V25-V32. DOI:10.1190/1.3552706 |
| [29] |
魏小强, 雷秀丽, 马庆珍. 基于多道奇异谱分析的三维地震数据规则化方法[J]. 石油地球物理勘探, 2014, 49(5): 846-850. WEI Xiaoqiang, LEI Xiuli, MA Qingzhen. 3D seismic data regularization based on multichannel singular spectrum analysis[J]. Oil Geophysical Prospecting, 2014, 49(5): 846-850. |
| [30] |
DONOHO DL, HUO X. Uncertainty principles and ideal atomic decomposition[J]. IEEE Transactions on Information Theory, 2001, 47(7): 2845-2862. DOI:10.1109/18.959265 |
| [31] |
HERRMANN F J, HENNENFENT G. Non-parametric seismic data recovery with curvelet frames[J]. Geophysical Journal International, 2008, 173(1): 233-248. DOI:10.1111/gji.2008.173.issue-1 |
| [32] |
MA J W. Improved iterative curvelet thresholding for compressed sensing[J]. IEEE Transactions on Instrumentation and measurement, 2011, 60(1): 126-136. DOI:10.1109/TIM.2010.2049221 |