2018, Vol. 61
2. 吉林大学地球探测科学与技术学院, 长春 130026
2. College of Geo-exploration Science and Technology, Jilin University, Changchun 130026, China
人工地震是地球物理勘探中重要的方法之一，在保证精度和分辨率的情况下，探测深度甚至能够达到莫霍面.目前的石油、天然气等化石能源的确定以地震勘探为主要手段，随着国家重点研发计划“深地资源勘查开采”重点专项的开展，人工地震方法在地球深部矿产资源探查中也发挥重要的作用.地震资料的处理和成像主要通过地表获取的地震波反射数据来反映地下介质的结构，因此采集地震数据的数量和质量成为获取地下介质准确图像的关键环节.理想情况下地震波场采样应该是规则且全空间覆盖的，但由于一些实际环境和经济因素的考虑，如地表障碍物的存在、地形条件、仪器硬件问题引起的采集坏道，以及海上数据采集时电缆的羽漂等情况影响，实际数据在空间方向上是非规则且稀疏的，对最终准确成像产生严重的影响.随着勘探程度的逐渐加深，高精度的三维地震勘探以及宽方位角地震勘探方法快速发展，大大增加了采集数据的信息量，但是获取完整地震数据体仍然是不可实现的任务.通过对地震数据进行插值重建，可以提供更多的地球物理信息，减小地球物理反问题的不确定性，满足对复杂地质构造进行精细刻画的要求.因此，地震数据插值问题的求解是从资料处理初期的去噪，如：三维表层相关多次波消除技术(SRME)(Verschuur et al., 1992)到后期的偏移处理，如：克希霍夫叠前深度偏移(Gardner and Canning, 1994)、共方位角偏移(Mosher et al., 1997)、三维真振幅波动方程偏移(WEM)(Zhang et al., 2003)，再到油田开发阶段的储层监控，如：四维时移地震勘探(Morice et al., 2000)等重要理论成功应用的前提.因此，开发快速、有效的地震数据插值方法具有重要的实际意义.现存的地震数据插值方法主要分为两类：基于图像处理的方法和基于波动方程的方法.
(1) 基于图像处理最简单的插值方法是面元重置，即改变数据的真正位置.但这种方法会产生假象，并且最终导致成像失败.带有部分正常时差校正(NMO)的面元重置方法结合了简单的运动学特征，可以对炮点和检波点位置的微分时差进行校正，改善一定的处理效果，但是，这种方法仍然受过于简单的NMO假设所制约.另一种方法是利用褶积算子求解迭代最优化问题(Claerbout，1992)，这种方法在处理大的数据间隔时具有较高的效率，但是，时空方向同相轴的变化使该方法面临很大的挑战，尽管时空域变系数预测滤波器可以带来比较精确的处理结果，但是在处理三维数据时需要存储大量的滤波器系数，使得该方法的应用受到较大限制(Liu and Fomel, 2011).分形插值数据重建方法(李信富和李小凡，2008)能够很好地恢复地震数据的局部信息，但是有时存在计算效率的问题.基于非一致离散傅里叶变换(刘喜武等，2004；Xu et al., 2005；孟小红等，2008)以及Radon变换(王维红等，2007)的插值方法仅在对应变换方法对特定数据提供稀疏表征的前提下能够获得较好的效果，同时还依赖于迭代算法的有效性.
(2) 基于波动方程动力学关系的插值方法包括基于不同类型积分连续算子的方法，例如：炮检距连续(Fomel，2003)和炮连续(Mazzucchelli et al., 1999)等.积分连续算子可以直接应用于缺失数据的插值问题(辛可锋等，2002；管路平等，2009)，但在计算小距离连续时，由于积分孔径有限，处理效果有时不是很理想.基于干涉测量的地震数据插值方法(Wang et al., 2009)主要利用地震数据动力学的可预测特性，但方法往往受限于全波场数据运算时所产生的庞大计算量.
压缩感知(Compressive sensing，CS)是由Donoho(2006)等在信号处理领域正式提出的概念，根本条件是信号本身具有可压缩性，该理论将数据空间向模型空间进行投影，在模型空间的系数只有少数元素为非零(即具有稀疏的特性)的前提下，通过求取数学反问题解决数据缺失问题.基于压缩感知理论的地震数据插值方法主要分两个部分：稀疏变换和迭代插值方法.很多学者对利用不同稀疏变换的插值方法进行了研究，如：结合傅里叶变换(Xu et al., 2005；Zwartjes and Sacchi, 2007；高建军等，2009；唐刚和杨慧珠，2010)、小波变换(Sweldens, 1995; Sweldens and Schröder, 1996; Wang et al., 2011)、curvelet变换(刘国昌等，2011)以及Seislet变换(刘洋等，2009；Fomel and Liu, 2010; Liu and Fomel, 2010；刘财等，2013)的插值方法在处理地震数据重构问题时取得了较好的效果.另一方面，有许多解决地震数据插值相关反问题的迭代算法，例如：凸集投影法(POCS)(Abma and Kabir, 2006；Gao et al., 2010)，零范数稀疏优化方法(曹静杰等，2012)，迭代阈值法(IST)(Yang and Fomel, 2015)，非线性整形迭代法(Chen et al., 2015)，Bregman迭代法(Gou et al., 2014)，不准确Uzawa方法(Liang et al., 2014)以及两步法(Yu et al., 2015)等.但是，已有的迭代算法往往难以同时保证插值精度和计算效率.本文仅针对迭代插值方法进行研究，在稀疏变换的选取上并未有过多的讨论.
本文提出一种新的Bregman整形迭代算法，并与传统的POCS以及IST迭代方法进行对比，分析不同迭代算法的插值特性，同时针对本文迭代算法的特点，提出两种匹配的迭代控制准则来进一步优化计算效率.以傅里叶变换为例，在理论模型和实际数据处理中与POCS方法进行比较，为有效解决地震缺失数据插值问题提供一种可供选择的迭代策略.
1 理论基础 1.1 地震数据插值的数学反问题地震数据重建可以归纳为线性估计问题，并且通过迭代最佳算法进行求解(Fomel，2001).假设d为观测数据，m为拟重建的模型数据，利用正向模型算子L可以获得模型数据重建方程：
|
(1) |
其中，
缺失数据插值是一类特殊的数据重建问题，已知的数据已经在规则的观测系统网格里，只需要恢复在空面元内的缺失数据即可.通常情况，将选择算子K作为正向模型算子，K为对角矩阵，满足如下形式：
|
(2) |
其中，公式(2)中的I代表已知的地震数据位置，为N×N的单位矩阵，而O代表缺失的地震数据位置，为N×N的零矩阵，K为(M×N)×(M×N)的矩阵.此时公式(1)转化为数学欠定方程：
|
(3) |
为了得到合理的插值结果m，需要增加额外的约束条件.Claerbout(1992)提出了缺失数据插值问题的基本准则，即缺失数据的插值方法能够使恢复的数据在指定滤波之后有最小的能量.这一准则在数学上可以表示为如下方程：
|
(4) |
其中，D为特定的滤波器.结合公式(3)和(4)可以建立正则化条件下的插值问题：
|
(5) |
式中，ε为标量，表示惩罚因子.
如果选取D为梯度算子，则相当于利用Tikhonov正则化条件(Tikhonov，1963)进行问题求解.Fomel(2002)提出利用平面波分解滤波器(PWD)作为正则化条件的技术思路，由于使用了地震同相轴的倾角信息，能够获得较好的插值效果，但是该方法往往受数据中噪声的影响.利用正则化条件对数学欠定问题进行约束求解后，反演问题的解为局部平滑，但是局部平滑的插值结果往往并非所需要的解，因此需要提出新的约束条件.
压缩感知是由Candès(2006)和Donoho(2006)正式提出的新的采样理论，其能够突破传统香农采样定理的局限，可以在远小于奈奎斯特采样率的条件下获取信号的离散样本，保证信号的无失真重建，压缩感知理论的核心思想主要包括：信号的稀疏结构和不相关特性.
对于反演问题来说，如果模型空间是可压缩的(即具有极少的非零元素)，那么压缩感知理论可以将求解公式(1)转化为如下问题：
|
(6) |
其中，‖ · ‖ 0表示L0范数，即数据的非零个数.公式(6)为数学非凸优化问题，并且是NP-hard问题(Natarajan，1995)，在解该问题时局部最优解往往不是全局最优解，因此在数学上往往把公式(6)转化为凸优化问题：
|
(7) |
进行求解，此时局部最优解同时也是全局最优解.其中，‖ · ‖ 1表示L1范数，公式(6)到公式(7)的转化为凸松弛法中的基追踪算法(Basis Pursuit，BP)，Elad和Bruckstein(2002)证明了公式(6)和(7)的等效条件.
尽管公式(7)为凸优化问题，但是其仍然是约束问题，在数学上难以求解，通常将其进一步转化为非约束问题：
|
(8) |
其中，μ为惩罚参数.数学上有很多方法可以求解公式(8)，例如迭代重加权最小二乘算法(IRLS)(Chartrand and Yin, 2008).然而，地震数据m很少存在可直接压缩(只有少数元素为非零)的情况.为了将反问题(1)在压缩感知理论框架下进行求解，需要对m进行变量代换.如果m在某种变换域内具有可压缩的特性，即
|
(9) |
其中，T-1为稀疏反变换，p为变换域系数，同时p的非零数据个数远远小于m的非零数据个数.将公式(9)代入公式(7)，可得
|
(10) |
根据公式(8)，压缩感知框架下反问题的近似解为
|
(11) |
最终，获得的变换域系数p通过公式(9)即可计算所需的模型m.
应用压缩感知理论求解地震数据插值问题时，正向模型算子L为选择算子K，则对应如下约束最优化问题
|
(12) |
和非约束问题
|
(13) |
在求解公式(12)和(13)时，关键在于稀疏变换T的设计和迭代算法的选取.常用的稀疏变换有：傅里叶变换，小波变换，Seislet变换，曲波变换(Curvelet)变换，Radon变换等.各种变换的适用条件不同，在进行问题求解时需要首先分析数据的模式，然后选取对应的数学变换以获得最佳的稀疏表征.例如：如果地震数据中仅包含平面波，则傅里叶变换具有更好的波场稀疏表征能力，而小波变换主要用于分析分片平滑数据，地震数据往往不具有这个特征(时间和空间采样率不同)，因此传统的小波变换并不适用于地震波场的稀疏表征.本文并不针对稀疏变换进行压缩能力方面的研究，因此选取工业界常用的傅里叶变换作为测试迭代算法的稀疏变换.
1.2 Bregman整形迭代方法许多学者提出了不同的算法用于解决非约束最小二乘问题，如凸集投影(POCS)和迭代收缩阈值方法(IST).然而，传统方法很难在重建精度和计算速度之间达到较好的平衡.为了求解非约束问题(公式11)，许多学者提出了不同的方法.从公式
|
(14) |
出发，稀疏正变换算子T可以看做KT-1的近似逆算子，即级联算子TKT-1约等于单位算子I，那么存在简单的迭代方法
|
(15) |
公式(15)收敛于如下方程的解：
|
(16) |
公式(15)被Goldin(1986)称为“R方法”，其收敛的充分条件是公式(16)右侧是稀疏的.如果T被KT-1的转置算子所代替，那么公式(15)为“Landweber迭代”(Landweber，1951).由于公式(15)往往不满足所期待的模型解，Fomel(2008)提出了一种非线性整形正则化方法来约束模型空间的属性，通过定义一个整形算子S，公式(15)被修改为
|
(17) |
其中，整形算子S可以选取为非线性阈值算子(Donoho and Johnstone, 1994)，软阈值算子的表达式为：
|
(18) |
其中，p表示变换域系数，λ为选取的阈值参数，sgn为符号函数，(|p|-λ)+表示当(|p|-λ)>0时等于|p|-λ，当(|p|-λ) < 0时等于0.公式(18)的作用是去除小于阈值的变换域系数.当T被替换为KT-1的转置算子时，该迭代方程(17)为Daubechies等(2004)提出的收缩迭代方法，即迭代收缩阈值方法(IST).然而，迭代阈值方法中线性算子(KT-1)的转置算子很难求取.因此，广义条件下的整形迭代公式(17)更加适合解决公式(13)对应的插值目标函数(Chen et al，2015).
尽管整形迭代方法能够求解公式(13)，但是插值的数据过于平滑，容易引起插值数据与已知数据的明显差异，导致插值后地震同相轴不连续.主要原因在于公式(12)与(13)并不等价，除非存在零解的情况.Osher等(2005)提出解决公式(12)对应约束最优化问题的Bregman迭代方法，其通过计算一系列凸优化问题(Yin et al., 2008)进行约束问题求解，Bregman的迭代框架为
|
该方法已经用于求解地震数据插值问题(刘财等，2013；Gou et al., 2014).Bregman迭代方法的优势在于，公式(20)中的参数μ为常数，所以该方法确保了迭代方法的快速迭代.但是以往方法利用简单的迭代阈值方法求解公式(20)，由于未考虑选择算子K的影响，因此插值效果一般.本次研究利用整形迭代公式(17)作为计算公式(20)的迭代算子，并在Bregman迭代框架下求解插值目标函数(公式12)，可得如下方程：
|
对公式(22)进一步整理，则Bregman整形迭代方程为
|
方程(23)和(24)类似于全变分(TV)去噪方法(Osher et al., 2005)中Rudin-Osher-Fatemi(ROF)模型的“加回残差”方法，最终的插值结果为
|
(25) |
其中，N为迭代次数.需要注意的是，在Bregman整形迭代公式(23)和(24)中，数据空间dk和模型空间pk是两套并行的迭代体系，并且dk≠T-1pk，而POCS和IST方法只有一套独立的迭代方程，因此，Bregman整形迭代与POCS以及IST存在本质的不同.
1.3 Bregman整形迭代与POCS以及IST迭代的比较Abma和Kabir(2006)将凸集投影(POCS)方法引入到地震数据插值问题中，并且对传统的POCS方法进行了改进，其利用非线性阈值算子代替传统POCS方法中的线性投影算子，该方法具有POCS方法的特性，但是从严格的数学意义上，该方法不是标准POCS方法.尽管如此，为了与工业界习惯统一，本文仍将该方法称为POCS迭代，其在数据空间进行迭代，计算公式如下：
|
(26) |
迭代收缩阈值(IST)方法(Blumensath and Davies, 2009)可以有效地求解反问题公式(13).利用稀疏约束方法就可以求取非平稳函数的最优解.该方法是在模型空间进行迭代，计算公式如下：
|
(27) |
其中，(·)H为转置运算.如果T为正交矩阵，例如傅里叶变换，则T=(T-1)-1=(T-1)H，并且pk=Tdk，则公式(27)可以转化为数据空间迭代公式：
|
(28) |
对比公式(25)，(26)和(28)，可以看到，POCS迭代过程中选取固定的阈值函数S只能满足公式(12)中数据空间重建或者模型空间稀疏的两个目标函数之一，因此迭代速度较慢，Gao等(2010)使用变阈值函数一定程度提高了POCS算法的计算效率.对比公式(26)和(28)不难发现，IST迭代与Abma和Kabir(2006)提出的POCS迭代方法是等价的，区别仅在于是否在迭代结束后将已知位置数据恢复，IST迭代由于改变了已知数据，容易造成有效信息的损失.Bregman整形迭代通过选取大的固定阈值函数，能够保证模型空间迭代(公式24)的快速收敛，尽管Bregman迭代也会改变已知位置的数据，但是公式(23)能够保证已知位置的计算数据与原始数据的差异最小，使已知位置与缺失位置的同相轴连续，保证数据重建精度.
1.4 迭代控制准则迭代插值方法常常通过迭代次数来控制计算过程，但是过小的迭代次数可能导致迭代未达到收敛，而过大的迭代次数可能会引起计算资源的浪费.Gao等(2013)提出了一种可应用于实际数据处理的迭代控制目标函数：
|
(29) |
方程(29)表示第k次迭代解与第k-1次迭代解之间的差异，但是该方法有可能收敛到局部最小值.Bregman整形迭代方法也可以使用公式(29)作为控制准则，但是由于本文方法的特殊性，可以构建更加合理的迭代控制准则，避免局部最小值的发生.Bregman整形迭代方法在计算过程中不断逼近已知位置的数据值，当迭代次数k趋近无穷大时，理论上能够完整重建已知观测数据，因此可以根据这一特性提出新的迭代控制准则：
|
(30) |
方程(30)是通过数据之间能量的差异来定义目标函数，除此以外，还可以使用数据的相似性进一步修改迭代控制准则.可以通过全局相关系数来定义：
|
(31) |
式中，〈, 〉表示两个数据之间的点积.当两个信号完全相同时，相关系数γ=1，如果完全相反时，γ=-1.根据Cauchy-Schwartz不等式，在其他情况下，相关系数的绝对值会小于1.方程(31)所示的全局相关能够衡量两个数据的整体相似性，为了分析数据之间局部的差异，还可以应用局部相似性(Fomel，2007)进一步扩展迭代控制准则，本文选取公式(31)定义的控制准则进行分析.
2 理论模型分析为了验证Bregman整形迭代方法的插值效果，首先选取一个简单的平面波理论模型(图 1a)测试本文算法的有效性.该模型包含4条振幅不同的线性同相轴，随机去除50%地震道构建缺失数据模型，如图 1b所示.选取傅里叶变换对波场进行稀疏表征，通过Bregman整形迭代方法对缺失数据进行重建，将最大迭代次数设为50，通过方程(29)、(30)和(31)分别求取收敛曲线.可以看出方程(29)所示的曲线(图 2a)在初始的几次迭代中收敛比较迅速，但是该方法容易使得迭代收敛到局部最小值.方程(30)所示的曲线在开始的几次迭代过程中误差也不是单调递减，尽管随着迭代的进行，该能量误差最终趋于稳定，但是图 2b表征方程(30)所示算法也容易收敛到局部最小值.方程(31)所示的相似性曲线(图 2c)显示随着迭代次数的增加，插值结果与原始不规则数据的相似性单调增加，并收敛到全局最优解，从图中可以直接判断迭代近似等于33次时，收敛到最优解(图 2d).对比三种迭代控制准则，相似性准则对插值误差要求最为严格.
|
图 1 理论模型 (a)原始理论模型；(b)随机去除50%地震道的模型数据. Fig. 1 Synthetic model (a) Ideal synthetic data; (b) Data removed 50% of randomly selected traces. |
|
图 2 收敛曲线和本文方法数据重建结果 (a) Jk1收敛曲线；(b) Jk2收敛曲线；(c)相似性收敛曲线；(d)第33次迭代数据重建结果. Fig. 2 Convergence curves and the interpolated result by proposed method (a) Convergence curve of Jk1; (b) Convergence curve of Jk2; (c) Convergence curve of similarity norm; (d) The interpolated result at the 33th iteration. |
选取带随机噪声数据对插值方法进行测试.在图 1a所示的理论模型中加入随机噪声(图 3a)，随机去除50%地震道模拟带有随机噪声的缺失数据模型(图 3b).依然选取最大迭代次数为50，分析方程(29)、(30)和(31)计算的收敛曲线所表现出的不同特征关系.从收敛曲线(图 4a)看出方程(29)定义的收敛判断准则依然容易收敛到局部解，而方程(30)和(31)所示的收敛曲线都表现出较好的单调性(图 4b和图 4c)，相似性准则的稳定性更好.当迭代次数达到23左右时，本文算法可以收敛到全局最优解，最终的重建结果(图 4d)表明，插值数据后的数据具有合理的信噪比分布.
|
图 3 带有随机噪声的理论模型 (a)带有随机噪声的理论模型；(b)随机去除50%地震道的理论模型. Fig. 3 Synthetic model with random noise (a) Synthetic data with random noise; (b) Synthetic data with 50% traces removed. |
|
图 4 收敛曲线和本文方法数据重建结果 (a) Jk1收敛曲线；(b) Jk2收敛曲线；(c)相似性收敛曲线；(d)第23次迭代数据重建结果. Fig. 4 Convergence curves and the interpolated result by proposed method (a) Convergence curve of Jk1; (b) Convergence curve of Jk2; (c) Convergence curve of similarity norm; (d) The interpolated result at the 23th iteration. |
为了验证Bregman整形迭代方法的优势，本文选取凸集投影(POCS)方法进一步对比.分别使用两种数据重建方法对图 1b所示的缺失数据进行数据重建，利用公式：
|
(32) |
分别计算POCS插值方法和Bregman整形迭代插值方法数据插值结果与图 1a所示理论完整数据之间的能量差异，式中dc表示理论完整数据，得到图 5所示的误差曲线.从图中可以看出本文算法不论是在收敛速度还是重建数据的精度上都要略优于POCS方法，而且当处理带有噪声的数据插值问题(图 3b)时，本文算法的优势更加明显.仅在初始几次迭代时POCS迭代的误差要小于Bregman整形迭代方法，原因在于POCS插值方法每次迭代之后都会将已知位置数据恢复，而Bregman整形迭代算法并没有该步操作，但是初始迭代时的误差较大，POCS迭代并未达到收敛.选取两种方法都能够达到收敛的第33次迭代结果进行更直观地说明(图 6a和图 6b)，从插值结果比较，两种迭代方法均具有较好的重建效果，但是对比插值误差时则能够明显地看到图 5所示的两种插值结果的差异，为了更好地突出对比，插值误差剖面(图 6c和6d)选取了与插值结果(图 6a和6b)不同的显示参数.POCS方法获得的插值误差剖面(图 6c)表明，POCS重建的数据与理论数据仍然存在较大的误差，而本文算法计算的误差(图 6d)远小于POCS方法，证明Bregman整形迭代算法保证更好的插值精度.
|
图 5 Bregman整形迭代与POCS迭代误差对比曲线 虚线表示POCS插值方法，实线表示Bregman整形迭代插值方法. Fig. 5 The comparison of interpolation error between Bregman shaping iteration and POCS Dash line is POCS and solid line is Bregman shaping iteration. |
|
图 6 不同插值结果比较 (a) POCS迭代插值结果；(b) Bregman整形迭代插值结果(同图 2d)；(c) POCS迭代插值误差(显示参数不同于图 6a和图 6b)；(d) Bregman整形迭代插值误差(显示参数同图 6c). Fig. 6 The comparison of different interpolated results (a) POCS method; (b) Bregman shaping method (same as Fig. 2d); (c) Interpolation error of POCS method (using different plotting parameter from Fig. 6a and 6b); (d) Interpolation error of Bregman shaping iteration (using same plotting parameter as Fig. 6c). |
在实际资料处理中，选取墨西哥湾某深水地震数据(Crawley et al., 1999; Fomel, 2002)测试本文算法的适用性.图 7a为三维炮集数据，图中正上部为2.9 s位置的时间切片，左下图为第20炮的单炮记录，右下图为50道位置所示的共炮检距剖面.该数据最大道数为180道.从图中可以看出该数据体存在长周期多次波和由盐丘体导致的复杂散射波，并且该数据同相轴的振幅不均匀分布以及存在一定量的随机噪声干扰.为了测试和比较插值方法，随机去除40%地震道，如图 7b所示.为了对比本文算法的有效性，首先选取基于傅里叶变换的POCS算法对缺失数据(图 7b)进行插值重建，为了获得更好的插值效果，POCS迭代次数选取为100，处理结果如图 7c所示，从图中可以看出POCS插值方法可以重建大部分反射波信息，但是在空间方向50~150道，时间位置3.2 s附近的散射波恢复的较差，另外可以很明显地看到插值后的地震道信噪比更高，与原始数据差异较大.在应用Bregman整形迭代算法处理时，可以利用公式(30)所定义的迭代准则对本文算法的迭代过程进行更加精确地控制，迭代收敛曲线如图 8a所示，可以判断，基于傅里叶变换的Bregman整形迭代方法在第25次附近收敛，对应的插值重建结果如图 7d所示，可以看到本文方法对振幅存在非平稳变化的同相轴也能实现合理的插值.对比POCS方法，Bregman整形迭代算法能够更好地恢复散射波，同时插值后的地震道具有更加真实的信噪比，插值后的数据同相轴连续性更好.进一步根据公式(32)计算两种方法的插值误差曲线(图 8b)，可以看出，尽管POCS方法在前几次迭代的误差要小于本文方法，但是其并未达到收敛，另外即使利用更多的迭代也无法获得理想的插值误差，而本文算法在很少的迭代条件下就能够获得较小的插值误差.为了进一步在细节处对比两种方法，求取了数据重建结果与理想数据(图 7a)的差剖面，如图 9所示.对比差剖面结果也可以看出POCS方法比本文算法有更大的插值误差，尤其表现在散射波位置，进一步说明本文方法较传统方法的优越性.在实际数据处理中，两种方法均是对整个数据体进行三维傅里叶变换，如果分时窗处理，能够进一步提高插值的精度.
|
图 7 实际数据测试结果 (a)三维实际数据；(b)随机去除40%地震道的数据；(c) POCS迭代插值结果；(d) Bregman整形迭代插值结果. Fig. 7 Field data test (a) A 3D field dataset; (b) Data with 40% randomly missing traces; (c) Interpolated result using POCS; (d) Interpolated result using Bregman shaping iteration. |
|
图 8 迭代收敛曲线及插值误差曲线 (a) Bregman整形迭代收敛曲线；(b)不同方法的插值误差曲线(虚线为POCS方法，实线为Bregman整形迭代方法) Fig. 8 Iteration convergence curve and interpolation error curves (a) Convergence curve of Bregman shaping iteration; (b) Iteration error curves by using different methods (Dash line is POCS method and solid line is Bregman shaping iteration method) |
|
图 9 不同插值结果与原始数据的差剖面 (a) POCS迭代；(b) Bregman整形迭代. Fig. 9 Difference between reconstructed data and original data by using different methods (a) POCS method; (b) Bregman shaping interpolation. |
本文提出了一种基于稀疏变换的Bregman整形迭代插值算法用于解决地震数据缺失问题.该方法利用有效的稀疏变换表征地震波场，在压缩感知理论下，利用Bregman整形迭代法求解约束最优化问题，对缺失数据进行重建.同时，在迭代过程中提出两种与插值算法相匹配的迭代控制准则，通过求取能量曲线或相似性曲线获得最优化迭代次数，在保证插值结果准确性的前提下节约计算资源.对理论模型和实际地震数据进行处理和分析，并且对比工业POCS方法，结果表明本文迭代方法能够提供更合理的插值结果，不论误差精度还是收敛速度均优于工业POCS插值方法.另外，本文迭代方法可以直接与其他稀疏变换进行匹配，如果选取能够表征非平面波数据的稀疏变换，将会进一步提高本文迭代算法的缺失数据重建效果.
致谢第一作者感谢美国德州大学奥斯汀分校经济地质局Sergey Fomel教授的讨论和建议.感谢两位匿名审稿专家提出的宝贵意见.
Abma R, Kabir N.
2006. 3D interpolation of irregular data with a POCS algorithm. Geophysics, 71(6): E91-E97.
DOI:10.1190/1.2356088 |
|
Blumensath T, Davies M E.
2009. Iterative hard thresholding for compressed sensing. Applied and Computational Harmonic Analysis, 27(3): 265-274.
DOI:10.1016/j.acha.2009.04.002 |
|
Candès E J. 2006. Compressive sampling. //Proceedings of the International Congress of Mathematicians. Madrid, Spain: European Mathematical Society, 1433-1452.
|
|
Cao J J, Wang Y F, Yang C C.
2012. Seismic data restoration based on compressive sensing using the regularization and zero-norm sparse optimization. Chinese Journal of Geophysics (in Chinese), 55(2): 596-607.
DOI:10.6038/j.issn.0001-5733.2012.02.022 |
|
Chartrand R, Yin W T. 2008. Iteratively reweighted algorithms for compressive sensing. //IEEE International Conference on Acoustics, Speech and Signal Processing. Las Vegas, NV, USA: IEEE, 3869-3872.
http://ieeexplore.ieee.org/document/4518498/ |
|
Chen Y K, Zhang L L, Mo L W.
2015. Seismic data interpolation using nonlinear shaping regularization. Journal of Seismic Exploration, 24(5): 327-342.
|
|
Claerbout J F. 1992. Earth Soundings Analysis: Processing Versus Inversion. Boston: Blackwell Scientific Publications.
|
|
Crawley S, Clapp R, Claerbout J. 1999. Interpolation with smoothly nonstationary prediction-error filters. //69th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 1154-1157.
http://www.researchgate.net/publication/2349182_Interpolation_With_Smoothly_Nonstationary_Prediction-Error_Filters |
|
Daubechies I, Defrise M, De Mol C.
2004. An iterative thresholding algorithm for linear inverse problems with a sparsity constraint. Communications on Pure and Applied Mathematics, 57(11): 1413-1457.
DOI:10.1002/(ISSN)1097-0312 |
|
Donoho D L, Johnstone I M.
1994. Ideal spatial adaption by wavelet shrinkage. Biometrika, 81(3): 425-455.
DOI:10.1093/biomet/81.3.425 |
|
Donoho D L.
2006. Compressed sensing. IEEE Transactions on Information Theory, 52(4): 1289-1306.
DOI:10.1109/TIT.2006.871582 |
|
Elad M, Bruckstein A M.
2002. A generalized uncertainty principle and sparse representation in pairs of bases. IEEE Transactions on Information Theory, 48(9): 2558-2567.
DOI:10.1109/TIT.2002.801410 |
|
Fomel S.
2002. Applications of plane-wave destruction filters. Geophysics, 67(6): 1946-1960.
DOI:10.1190/1.1527095 |
|
Fomel S.
2003. Seismic reflection data interpolation with differential offset and shot continuation. Geophysics, 68(2): 733-744.
DOI:10.1190/1.1567243 |
|
Fomel S.
2007. Local seismic attributes. Geophysics, 72(3): A29-A33.
DOI:10.1190/1.2437573 |
|
Fomel S. 2008. Nonlinear shaping regularization in geophysical inverse problems. //78th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts2046-2051.
http://scitation.aip.org/getabs/servlet/GetabsServlet?prog=normal&id=SEGEAB000027000001002046000001&idtype=cvips&gifs=Yes |
|
Fomel S, Liu Y.
2010. Seislet transform and seislet frame. Geophysics, 75(3): V25-V38.
DOI:10.1190/1.3380591 |
|
Fomel S B. 2001. Three-dimensional seismic data regularization. Stanford: Stanford University.
|
|
Gao J J, Chen X H, Li J Y, et al.
2009. Study on reconstruction of seismic data based on nonuniform Fourier transform. Progress in Geophysics (in Chinese), 24(5): 1741-1747.
|
|
Gao J J, Chen X H, Li J Y, et al.
2010. Irregular seismic data reconstruction based on exponential threshold model of POCS method. Applied Geophysics, 7(3): 229-238.
DOI:10.1007/s11770-010-0246-5 |
|
Gao J J, Stanton A, Naghizadeh M, et al.
2013. Convergence improvement and noise attenuation considerations for beyond alias projection onto convex sets reconstruction. Geophysical Prospecting, 61(S1): 138-151.
|
|
Gardner G H F, Canning A. 1994. Effects of irregular sampling on 3-D prestack migration. //64th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 1553-1556.
http://www.researchgate.net/publication/249854823_effects_of_irregular_sampling_on_3-d_prestack_migration |
|
Goldin S V. 1986. Seismic Traveltime Inversion. Tulsa, Okla: SEG.
|
|
Gou F Y, Liu C, Liu Y, et al.
2014. Complex seismic wavefield interpolation based on the Bregman iteration method in the sparse transform domain. Applied Geophysics, 11(3): 277-288.
DOI:10.1007/s11770-014-0443-3 |
|
Guan L P, Tang Y X, Wang H Z.
2009. Common-offset plane-wave prestack time migration and demigration. Chinese Journal of Geophysics (in Chinese), 52(5): 1301-1309.
|
|
Landweber L.
1951. An iteration formula for Fredholm integral equations of the first kind. American Journal of Mathematics, 73(3): 615-624.
DOI:10.2307/2372313 |
|
Li X F, Li X F.
2008. Seismic data reconstruction with fractal interpolation. Chinese Journal of Geophysics (in Chinese), 51(4): 1196-1201.
|
|
Liang J W, Ma J W, Zhang X Q.
2014. Seismic data restoration via data-driven tight frame. Geophysics, 79(3): V65-V74.
DOI:10.1190/geo2013-0252.1 |
|
Liu C, Li P, Liu Y, et al.
2013. Iterative data interpolation beyond aliasing using seislet transform. Chinese Journal of Geophysics (in Chinese), 56(5): 1619-1627.
DOI:10.6038/cjg20130519 |
|
Liu G C, Chen X H, Guo Z F, et al.
2011. Missing seismic data rebuilding by interpolation based on Curvelet transform. Oil Geophysical Prospecting (in Chinese), 46(2): 237-246.
|
|
Liu X W, Liu H, Liu B.
2004. A study on algorithm for reconstruction of de-alias uneven seismic data. Chinese Journal of Geophysics (in Chinese), 47(2): 299-305.
|
|
Liu Y, Fomel S, Liu C, et al.
2009. High-order seislet transform and its application of random noise attenuation. Chinese Journal of Geophysics (in Chinese), 52(8): 2142-2151.
|
|
Liu Y, Fomel S.
2010. OC-seislet:seislet transform construction with differential offset continuation. Geophysics, 75(6): WB235-WB245.
DOI:10.1190/1.3479554 |
|
Liu Y, Fomel S.
2011. Seismic data interpolation beyond aliasing using regularized nonstationary autoregression. Geophysics, 76(5): V69-V77.
DOI:10.1190/geo2010-0231.1 |
|
Mazzucchelli P, Rocca F, Spagnolini U. 1999. Regularizing land acquisition using shot continuation operators: effects on amplitudes. //69th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 1995-1998.
http://www.researchgate.net/publication/240735663_Regularizing_land_acquisition_using_shot_continuation_operators_Effects_on_amplitudes |
|
Meng X H, Guo L H, Zhang Z F, et al.
2008. Reconstruction of seismic data with least squares inversion based on nonuniform fast Fourier transform. Chinese Journal of Geophysics (in Chinese), 51(1): 235-241.
|
|
Morice S, Ronen S, Canter P, et al. 2000. The impact of positioning differences on 4D repeatability. //70th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 1611-1614.
http://scitation.aip.org/getabs/servlet/GetabsServlet?prog=normal&id=SEGEAB000019000001001611000001&idtype=cvips&gifs=Yes |
|
Mosher C C, Foster D J, Hassanzadeh S. 1997. Common angle imaging with offset plane waves. //67th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 1379-1382.
http://scitation.aip.org/getabs/servlet/GetabsServlet?prog=normal&id=SEGEAB000016000001001379000001&idtype=cvips&prog=normal |
|
Natarajan B K.
1995. Sparse approximate solutions to linear systems. SIAM Journal on Computing, 24(2): 227-234.
DOI:10.1137/S0097539792240406 |
|
Osher S, Burger M, Goldfarb D, et al.
2005. An iterative regularization method for total variation-based image restoration. Multiscale Modeling & Simulation, 4(2): 460-489.
|
|
Sweldens W. 1995. Lifting scheme: A new philosophy in biorthogonal wavelet constructions. //Proceedings Volume 2569, Wavelet Applications in Signal and Image Processing Ⅲ. San Diego, CA, United States: SPIE, 68-79.
http://spie.org/Publications/Proceedings/Paper/10.1117/12.217619 |
|
Sweldens W, Schröder P. 1996. Building your own wavelets at home. //Klees R, Haagmans R. Wavelets in Computer Graphics. New York: ACM, 15-87.
http://www.springerlink.com/content/5672g2k5m21g7044 |
|
Tang G, Yang H Z.
2010. Seismic data compression and reconstruction based on Poisson Disk sampling. Chinese Journal of Geophysics (in Chinese), 53(9): 2181-2188.
DOI:10.3969/j.issn.0001-5733.2010.09.018 |
|
Tikhonov A N.
1963. Solution of incorrectly formulated problems and the regularization method. Soviet Mathematics Doklady, 4: 1035-1038.
|
|
Verschuur D J, Berkhout A J, Wapenaar C P A.
1992. Adaptive surface-related multiple elimination. Geophysics, 57(1): 1166-1177.
|
|
Wang W H, Pei J Y, Zhang J F.
2007. Prestack seismic data reconstruction using weighted parabolic Radon transform. Chinese Journal of Geophysics (in Chinese), 50(3): 851-859.
|
|
Wang Y B, Luo Y, Schuster G T.
2009. Interferometric interpolation of missing seismic data. Geophysics, 74(3): S137-S145.
|
|
Wang Y F, Cao J J, Yang C C.
2011. Recovery of seismic wavefields based on compressive sensing by an l1-norm constrained trust region method and the piecewise random subsampling. Geophysical Journal International, 187(1): 199-213.
DOI:10.1111/gji.2011.187.issue-1 |
|
Xin K F, Wang H Z, Wang C L, et al.
2002. Regularization of pre-stack seismic data. Oil Geophysical Prospecting (in Chinese), 37(4): 311-317.
|
|
Xu S, Zhang Y, Pham D, et al.
2005. Antileakage Fourier transform for seismic data regularization. Geophysics, 70(4): V87-V95.
DOI:10.1190/1.1993713 |
|
Yang P L, Fomel S.
2015. Seislet-based morphological component analysis using scale-dependent exponential shrinkage. Journal of Applied Geophysics, 118: 66-74.
DOI:10.1016/j.jappgeo.2015.04.003 |
|
Yin W T, Osher S, Goldfarb D, et al.
2008. Bregman iterative algorithms for l1 minimization with applications to compressed sensing. SIAM Journal on Imaging Sciences, 1(1): 143-168.
DOI:10.1137/070703983 |
|
Yu S W, Ma J W, Zhang X Q, et al.
2015. Interpolation and denoising of high-dimensional seismic data by learning a tight frame. Geophysics, 80(5): V119-V132.
DOI:10.1190/geo2014-0396.1 |
|
Zhang Y, Zhang G Q, Bleistein N.
2003. True amplitude wave equation migration arising from true amplitude one-way wave equations. Inverse Problems, 19(5): 1113-1138.
DOI:10.1088/0266-5611/19/5/307 |
|
Zwartjes P M, Sacchi M D.
2007. Fourier reconstruction of nonuniformly sampled, aliased seismic data. Geophysics, 72(1): V21-V32.
DOI:10.1190/1.2399442 |
|
曹静杰, 王彦飞, 杨长春.
2012. 地震数据压缩重构的正则化与零范数稀疏最优化方法. 地球物理学报, 55(2): 596–607.
DOI:10.6038/j.issn.0001-5733.2012.02.022 |
|
高建军, 陈小宏, 李景叶, 等.
2009. 基于非均匀Fourier变换的地震数据重建方法研究. 地球物理学进展, 24(5): 1741–1747.
|
|
管路平, 唐亚勋, 王华忠.
2009. 共偏移距道集平面波叠前时间偏移与反偏移. 地球物理学报, 52(5): 1301–1309.
|
|
李信富, 李小凡.
2008. 分形插值地震数据重建方法研究. 地球物理学报, 51(4): 1196–1201.
|
|
刘财, 李鹏, 刘洋, 等.
2013. 基于seislet变换的反假频迭代数据插值方法. 地球物理学报, 56(5): 1619–1627.
DOI:10.6038/cjg20130519 |
|
刘国昌, 陈小宏, 郭志峰, 等.
2011. 基于Curvelet变换的缺失地震数据插值方法. 石油地球物理勘探, 46(2): 237–246.
|
|
刘喜武, 刘洪, 刘彬.
2004. 反假频非均匀地震数据重建方法研究. 地球物理学报, 47(2): 299–305.
|
|
刘洋, FomelS, 刘财, 等.
2009. 高阶seislet变换及其在随机噪声消除中的应用. 地球物理学报, 52(8): 2142–2151.
|
|
孟小红, 郭良辉, 张致付, 等.
2008. 基于非均匀快速傅里叶变换的最小二乘反演地震数据重建. 地球物理学报, 51(1): 235–241.
|
|
唐刚, 杨慧珠.
2010. 基于泊松碟采样的地震数据压缩重建. 地球物理学报, 53(9): 2181–2188.
DOI:10.3969/j.issn.0001-5733.2010.09.018 |
|
王维红, 裴江云, 张剑锋.
2007. 加权抛物Radon变换叠前地震数据重建. 地球物理学报, 50(3): 851–859.
|
|
辛可锋, 王华忠, 王成礼, 等.
2002. 叠前地震数据的规则化. 石油地球物理勘探, 37(4): 311–317.
|
|