一种频率域三维不规则地震数据重建方法
【技术领域】
[0001] 本发明涉及的是一种缺失道地震数据重建的方法，具体是一种频率域三维不规则 地震数据重建方法。 技术背景
[0002] 在野外数据采集过程中，由于采集设备、项目成本以及地形条件的限制，地震数据 沿空间方向上通常低于奈奎斯特采样定理要求进行不规则采样，从而导致采集到的地震数 据不规则、不完整，出现空间假频，难以满足后续处理的要求，例如速度分析、自由表面多次 波消除、偏移归位等，在这种情况下，需要发展较好的地震数据重建方法，重建出缺失的地 震道，得到满足处理要求的地震数据，并且该方法也可以指导野外数据采集，压缩地震数据 采集量。
[0003]基于曲波变换的数据重建方法不需要地下结构的先验信息，能够重建出规则缺失 和不规则缺失的地震道，且计算速度快，精度高，得到了广泛的应用。薛念等采用阈值迭代 法，实现了基于曲波变换的二维地震数据重建(基于Curvelet变换的地震数据插值和去噪 [硕士论文]，西南交通大学，2010)，刘国昌等采用凸集投影算法，在曲波域实现了基于P0CS 算法的二维地震数据重建(基于curvelet变换的缺失地震数据插值方法[J].石油地球物理 勘探，2011，46(2)，237-245。），该方法更能够重建非线性同相轴，更加适合具有各项异性特 征的地震波场，但没有深入研究基于曲波变换的三维数据重建方法，而且对阈值参数选择 也没有进一步研究，仅仅采用指数阈值参数公式，其重建后的效果还是有限。而后张华等在 此基础上实现了基于曲波变换的三维地震数据重建(基于jitter采样和曲波变换的三维地 震数据重建[J].地球物理学报，2013,56(5): 1637-1649。），通过对时间切片进行逐次重建， 取得较好的重建效果，并且探讨了阈值参数选择策略。王本锋等提出了不同阈值下的曲波 重建方法(P0CS联合改进的Jitter采样理论曲波域地震数据重建方法[J].地球物理学报， 2015,50(1):20-27。），然而该方法在时间域进行重建，且只采用硬阈值参数进行处理，当重 建数据同时包含高和低能量振幅同相轴时，由于P0CS算法在前几次迭代过程中，低于阈值 参数以下的同相轴系数被滤除，因此这些方法在重建低振幅有效波同相轴的效果有限，不 能完全恢复出地震信号的低能量同相轴及提高局部异常区的信噪比，并且迭代次数多，重 建效率低，制约了该方法的进一步应用。
【发明内容】
[0004]本发明的目的是为了能够快速高精度的重建出缺失道三维地震数据，而提供了一 种频率域三维不规则地震数据重建方法。
[0005]本发明技术方案：
[0006]本发明提出频率域三维不规则地震数据重建方法，首先利用傅立叶变换将时间域 三维地震数据变换到频率域，然后采用凸集投影算法思想和引入能够刻画地震数据局部化 特征的曲波变换。在迭代过程中提出新的按指数规律衰减的阈值参数，采用软阈值算子，对 每个频率切片单独进行重建，降低了迭代次数和提高了重建精度，从而达到重建三维不规 则地震数据的目的。
[0007] 一种频率域三维不规则地震数据重建方法，地震数据的重建问题描述为由一组不 完整数据通过线性算子的作用恢复出完整的数据，假设如下线性正演模型：
[0008] y〇bs=Mf
[0009] 这里yobseR"代表采集的地震数据;feRN，且N2n，表示待重建的无假频数据;Me RnXN表示对角矩阵，其元素1和0分别表示已知地震道和未知地震道，假设系数X是f在曲波域 C中的稀疏表示，则上述方程为：
[0010]
[0011] 这里上标Η代表共辄转置矩阵，当从采集的数据yQbs中重建无假频数据f时，由于X 是稀疏的，从而可以采用稀疏算法求解该欠定方程。
[0012] 在稀疏促进反演后，重建信号由；7 = 决定，同时
[0013]
[0014] 在这个表达式中，i代表估计值，Li范数定义为K ，x[i]是向量X中第i个 元素，通过求解上述方程，原始无假频的数据就重建出来；
[0015]在地震数据重建过程中，采用凸集投影算法，具体重建步骤如下：
[0016]步骤1:在时间域中输入具有缺失道的三维地震数据，然后采用傅立叶变换，将重 建数据从时间域变换到频率域；
[0017]步骤2:采用曲波变换对频率域三维地震数据进行稀疏变换，得到曲波域中曲波系 数，根据曲波系数的大小选择合适的阈值参数公式M(i=l，2,3，· · ·，Ν，其中N为迭代次 数)及其公式；
[0018] 步骤3:在曲波域中，采用软阈值算子进行处理，也即大于阈值人,的曲波系数减去 一个阈值大小的值，小于阈值的曲波系数加上一个阈值大小的值，而其它的曲波系数置 零；
[0019] 步骤4:将阈值化后的曲波系数做反曲波变换得到时间域地震数据，然后再用缺失 地震数据的未缺失地震道填充到反变换后的地震数据中去；
[0020] 步骤5:最后将得到的地震数据代入步骤(2)，重新进行迭代，直到运行Ν次结束，再 对迭代Ν次后的地震数据做反傅立叶变换即得到最终的重建结果。
[0021] 进一步，所述曲波变换的定义为：
[0022]
[0023] 式中：表示曲波函数，j，l，k分别表示尺度，方向和位置参数，f(x)为地震数 据，其频率域定义式为：
[0024]
[0025] 经过变换后得到的曲波系数，可用(：{」}{1}仏1氺2)表示其结构，其中」表示尺度，1 表示方向，（lu，k2)表示j尺度1方向上的矩阵系数。
[0026] 进一步，所述软阈值算子，其表达式如下：
[0027]
[0028] Tk表示软阈值算子，Sk-i表示第k-Ι次迭代得到的重建数据的曲波系数，满足 L=G7(^Ft和F-\表示关于时间变量t的正反傅立叶变换，74._,表示第k_l次迭代得到的 时间域重建数据，其中表示N维阈值集合，· · ·，λΝ}，且满足λρλ；^ > · · · >λΝ，Ν表示最大迭代次数。
[0029] 进一步，阈值参数公式，其表达式如下：
[0030]
[0031]其中Max为|CFy|的最大值，即曲波变换系数绝对值的最大值。Ν为总的迭代次数，ε为接近零的小值，与数据中噪声的能量有关。
[0032]本发明的优点：本发明通过利用傅立叶变换将时间域三维地震数据变换到频率 域，并采用凸集投影算法和能够刻画地震数据局部化特征的曲波变换进行数据重建，从而 避免了低能量有效波同相轴系数被滤除，提高了缺失地震数据重建精度，同时在迭代过程 中提出新的按指数规律衰减的阈值参数，采用软阈值算子，对每个频率切片单独进行重建， 可以减少迭代次数，降低算法的计算复杂度，大幅度地提高了计算效率，节省了运算时间。
【附图说明】
[0033]图1表示本发明实施例中频率域三维不规则地震数据重建方法流程图。
[0034] 图2不同阈值参数重建信噪比曲线对比图。
[0035] 图3为二维随机欠采样及其重建结果对比图。
[0036]图4为重建前后二维频谱对比图。
[0037]图5为信噪比与采样率关系曲线对比图。
[0038]图6为含噪数据重建过程图。
【具体实施方式】
[0039]以下实施案例用于说明本发明，但不用来限制本发明的范围。
[0040] 实施例1
[0041] 实现该方法的步骤主要包括，重建方程的构建，时间域变换到频率域，凸集投影重 建算法，阈值参数处理等。具体步骤如下：
[0042] 步骤1:重建方程的构建。地震数据的重建问题可以描述为由一组不完整数据通过 线性算子的作用恢复出完整的数据，假设如下线性正演模型
[0043] y〇bs=Mf
[0044] 这里yobseR"代表采集的地震数据;feRN，且N2n，表示待重建的无假频数据;Me RnXN表示对角矩阵，其元素1和0分别表示已知地震道和未知地震道，假设系数X是f在曲波域 C中的稀疏表示，则上述方程可以写成
[0045] y〇bs=Axiiα=μ€η
[0046]这里上标H代表共辄转置矩阵，当从采集到的数据yQbs中重建无假频数据f时，由于 X是稀疏的，从而使得该欠定方程有解。
[0047] 在稀疏促进反演后，重建信号由7 = (:〃 ^决定，同时
[0048]
r
[0049] 在这个表达式中，^代表估计值，Li范数定义为|χ| ，x[i]是向量x中第i 个元素，通过求解上述欠定方程，原始完整而规则的地震数据就可以重建出来。
[0050] 步骤2:时间域变换到频率域。本发明采用傅立叶变换对时间域地震数据进行计 算，其离散傅立叶变换公式如下：
[0051]
[0052] 式中：m、n均为整数。m= 0，1，"·Μ-1;n= 0，1，"·Ν-1 为频谱分析的时窗长 度。
[0053]步骤3:凸集投影重建算法。针对传统重建算法重建精度低和计算速度慢的缺点， 提出在频率域中进行重建的算法，并采用软阈值算子和新的按指数规律衰减的阈值参数公 式，其算法迭代表达式为：
[0054]
[0055]其中，fk表示第k次迭代得到的时间域重建数据，表示原始采集到数据y〇bs(t，x， y)，满足/n，Ft和F、表示关于时间变量t的正反傅立叶变换，c和(Γ1表示正反曲波变 换，其曲波变换的定义为：
[0056]
[0057] 式中：表示曲波函数，j，l，k分别表示尺度，方向和位置参数，f(x)为地震数 据，其频率域定义式为：
[0058]
[0059]经过变换后得到的曲波系数，可用(：{」}{1}仏1氺2)表示其结构，其中」表示尺度，1 表示方向，（lu，k2)表示j尺度1方向上的矩阵系数。
[0060] Tk表示软阈值算子，其元素满足：
[0061]
[0062]Sk-i表示第k-Ι次迭代得到的重建数据