基于一种改进凸集投影方法的地震数据同时插值和去噪

曹静杰， 王本锋

1 石家庄经济学院， 石家庄 050031 2 中国石油大学(北京)，油气资源与勘探测国家重点实验室， 北京 102249



基于一种改进凸集投影方法的地震数据同时插值和去噪

曹静杰1， 王本锋2

1 石家庄经济学院， 石家庄 050031 2 中国石油大学(北京)，油气资源与勘探测国家重点实验室， 北京 102249

基于稀疏反演的地震插值方法是一种重要的插值方法，然而大多数这类方法只针对无噪声数据或者高信噪比数据插值.实际上，地震数据含有各种噪声，使得插值问题变得更加困难.凸集投影方法是一种高效的插值算法，但是对于含噪声数据的插值效果不理想，针对含噪声数据提出的加权凸集投影方法能够实现同时插值和去噪，但是除了最小阈值需要认真选取外，增加一个权重因子来实现去噪功能.本文由迭代阈值算法推导出加权凸集投影方法，证明其是解无约束优化问题的一种方法，加权因子可以看作拟合误差项的系数.本文还提出了一种改进的凸集投影方法，与原始凸集投影方法相比该方法不需要增加任何计算量，只要通过阈值的选择来进行插值和去噪.数值模拟证明了该算法的计算效率，并且对含噪声数据能够实现较好的插值效果；先插值后去噪的结果证明了同时去噪和插值算法的可靠性和稳定性.

稀疏反演； 去噪； 插值； 反演

1 引言

由于采集成本、地形、坏道等因素的限制，野外采集的地震数据经常不满足香农采样定理(陆基孟， 1993). 不完整的数据会产生采集脚印，严重地影响叠前时间偏移(Symes，2007)、三维表面多次波消除(Berkhout and Verschuur，1997)和全波形反演(Bunks et al.，1995)的效果，因此需要采用插值方法来得到完备的波场信息. 地震数据插值不仅能够提供可靠的波场数据，还可以指导采集设计，节约采集成本(Yilmaz，1987).目前，地震插值方法主要包括预测滤波方法(Spitz， 1991)，基于稀疏变换的方法(Liu and Sacchi，2004；Zwartjes，2005；Herrmann and Hennenfent，2008；Xu et al.，2010)，基于波动方程的方法(Ronen，1987；Stolt，2002)和矩阵降秩方法(Oropeza and Sacchi，2011；Kreimer and Sacchi，2012；Ma，2013). 预测滤波方法专门针对规则不完备采样数据的插值(Spitz，1991)；基于波动方程的方法需要地下的参数信息并且计算量较大. 随着压缩感知理论引入到图像和三维数据的填充，矩阵/张量完备方法也被用于地震数据插值(Oropeza and Sacchi，2011；Kreimer and Sacchi，2012；Ma， 2013)，该方法将插值问题看作图像填充问题，基于完整的矩阵或张量的秩最小化来实现地震数据插值. 基于稀疏变换的方法利用各种变换，比如 Fourier 变换(Liu and Sacchi，2004；Zwartjes，2005；Xu et al.， 2010)，Radon 变换(Trad et al.，2002)，seislet 变换(Fomel and Liu，2010)，Dreamlet变换(耿瑜等，2012)，curvelet 变换(Cao et al.，2011；Wang et al., 2011；曹静杰等，2012；Mansour et al.，2013)等作为稀疏变换来压缩地震数据，由于其高效性和稳定性，该方法得到了很大的发展.

野外地震数据经常包含不同类型的噪声，这些噪声可以分为两类：随机噪声和相干噪声.目前的插值方法很少直接研究含噪声数据的插值.对于含随机噪声的采样不足数据，一般是先插值后去噪. Hennenfent 和 Herrmann (2006)提出了一种基于非规则采样curvelet 变换和阈值法的地震数据插值和去噪；Oropeza 和Sacchi (2011)提出了一种同时插值和去噪的重加权替换方法，通过加权来消除随机噪声，Gao 等(2013)提出一种加权凸集投影(POCS)方法来同时去噪和插值.

由于地震数据插值需要处理大量数据， 因此属于大规模计算问题， 需要快速稳定的插值算法. 目前凸集投影(POCS)算法是工业界一种重要的插值方法(Abma and Kabir，2006)，该方法计算效率很高，容易实现，并且具有很好的插值效果.但是对含噪声数据的插值，由于采样的噪声被带入到插值结果，因此张华和陈小宏(2013)证明该方法的结果不太理想. 加权POCS算法作为其改进能够实现去噪功能，但是除了需要调节最小的阈值外，增加的权重因子同样也需要认真选取.

本文针对含随机噪声数据的插值问题进行研究.我们从一般的0范数正则化和优化问题出发，由迭代阈值方法推导出了加权POCS方法，证明了加权POCS方法是解无约束优化问题的一种方法，权重因子可以看作是拟合误差项的系数(Oropeza and Sacchi，2011；Gao et al.，2013).同时，我们还提出了一种改进的POCS方法来实现同时插值和去噪，该方法不需要增加任何的计算量，能够适合含噪声数据的插值.数值模拟表明，改进的POCS方法在保持POCS方法计算效率高的同时，能够很好地重建无噪声的数据和含噪声的数据.通过和先插值后去噪的结果比较，证明了同时插值和去噪的稳定性.

2 地震数据的同时插值和去噪

假设地震数据的采集过程可以表示为以下的数学模型：

Φd+n=dobs，

(1)

其中Φ是采样算子矩阵，d是完整的波场数据，n为可加性噪声，此处假设为随机白噪声，dobs为采集的地震数据. 如果完整的波场数据d在某个变换域是稀疏的，即x=ΨTd是某个变换域的系数， 是稀疏的或可压缩的，其中ΨT是一个正交变换或者紧框架，则公式(1)可以表示为

ΦΨx+n=dobs.

(2)

定义A=ΦΨ，则公式 (2)的解可由以下优化问题得到

(3)

其中λ是用来平衡拟合误差项和x的稀疏性约束项的正则参数，R(x)是关于x的算子，是对解的稀疏性的描述. 地震数据插值属于大规模计算问题，因此需要快速稳定的算法. 在众多稀疏反演方法中，迭代阈值类算法是一种有效的解法. 如果R(x)=‖x‖1，则可以用迭代软阈值法解问题(3)(Daubechiesetal.，2004)； 如果R(x)=‖x‖0，则迭代硬阈值法是解问题(3)的一个很好的策略(BlumensathandDavies，2008).理论上，存在唯一的λ使得问题 (3) 能得到最优解(王彦飞，2007). 但是，迭代阈值类方法属于局部优化方法，为了防止得到局部解，在初始迭代时参数λ应当选择较大的值，然后逐渐减小到合适的值.

现有的基于问题(3)的大多数插值方法只是研究无噪声数据或高信噪比数据的插值，实际上问题(3) 同样适合于含有随机噪声数据的插值，如果参数λ选择合适就能够实现同时插值和去噪.参数λ和迭代阈值类算法中的阈值关系紧密，直接影响阈值算法的结果.迭代阈值类算法一般都是在变换域中进行阈值运算，本质上，变换域中的阈值处理是一种去噪过程，因此基于阈值类算法的插值可以看作变换域的去噪. 参数λ与噪声的能量有关，较大的噪声对应于较大的λ. 如果λ选择合适，则能够达到插值效果和去噪的最佳平衡. 然而不同数据的噪声能量是不同的，因此一般很难选择精确的λ.

Oropeza 和 Sacchi (2011)提出一种基于多道奇异谱分析的重加权替换方法，其中地震数据被重新排列为块Hankel/Toeplitz矩阵，通过降低Hankel/Toeplitz矩阵的秩来进行插值，秩的作用和λ的作用类似. 该算法同样通过逐渐降低秩来实现插值，因此多道奇异谱分析方法的思想和迭代阈值法的思想是一致的.

3 加权POCS方法和改进的POCS方法

凸集投影法是地震插值领域最常用的方法之一，该方法计算效率高，结果稳定，非常适合大规模计算，但该方法对含噪声数据的插值效果不理想，为此Gao等(2013) 提出一种加权POCS方法来插值含噪声数据.下面我们对加权POCS方法进行推导，证明其是解无约束优化问题的一种方法，权重因子可以看作是拟合误差项的系数.Yang等(2012)证明了对于插值问题，迭代硬阈值法可以推导出POCS方法.我们借鉴其中的思想推导出加权POCS方法.POCS方法可以用硬阈值迭代，也可以用软阈值迭代，取决于选用的正则算子是0范数还是1范数. 本文采用硬阈值进行推导，软阈值的推导过程和硬阈值的推导过程是一样的.迭代硬阈值法的迭代公式为

(4)其中N是迭代次数，xk为第k个迭代解，

是硬阈值运算，τ是阈值.τk随着迭代次数的增大而减小.τk的下降方式对插值结果有一定的影响，但不是关键的，关键的参数是最小阈值(正则参数)的大小，对于含噪声数据，应该仔细地选择最小阈值，对于无噪声数据，最小阈值可以很小.τk的下降方式可以采用线性下降或者指数下降，本文采用指数下降方式. 根据公式(4)，迭代硬阈值方法逐渐更新变换域的系数xk，当得到最终的迭代解xN时，只要通过反变换就能得到时间-空间域的插值结果.POCS方法的迭代公式为

(5)

其中I是元素全部为1的矩阵，该方法通过更新时间-空间域的解uk来得到最优解. 如果POCS法的初始解为0，则POCS方法第一次迭代的解为dobs，因此POCS方法的第一次迭代就得到一个好的初始解. 虽然POCS方法插值的信噪比在前几次迭代提高的很快，但是多次迭代后其最终解和迭代硬阈值法的解是一样的(Yangetal.，2012).

加权POCS方法利用加权的采样数据和加权的迭代解来构造新的迭代解(Gaoetal.，2013)，其中的权重因子是一个关键的参数. 下面给出其推导过程，定义

(6)

=Hτk(xk-αAT(Axk-dobs))，

(7)

其中τk为第k个阈值参数. 对公式(7)两边乘以Ψ，并且定义dk=Ψxk，可以得到

dk+1=Ψxk+1

=ΨHτk(xk+α(ΦΨ)T(dobs-Φdk))

=ΨHτk(ΨTdk+α(ΦΨ)T(dobs-Φdk))

=ΨHτk(ΨTdk+αΨT(ΦTdobs-ΦTΦdk))

=ΨHτk(ΨT(dk+αdobs-αΦdk))

=ΨHτk(ΨT(αdobs+(I-αΦ)dk)).

(8)

关于(8)的推导可以参考文献Yang等(2012). 定义uk=αdobs+(I-αΦ)dk，则公式 (8)可以变为

dk+1=ΨHτk(ΨT(uk)).

(9)

因此，

uk+1=αdobs+(I-αΦ)dk+1

=αdobs+(I-αΦ)dk+1

=αdobs+(I-αΦ)ΨHτk(ΨT(uk))

=αdobs+(I-αΦ)ΨHτk(ΨT(uk)).

(10)这即为Gao等(2013)文中的加权POCS的迭代公式.

根据上述推导过程，权重因子α实际上为问题(6)中拟合误差项的系数. 若问题(6)两边除以α，则问题 (6)可以变为

(11)

很明显，正则算子的系数λ变为λ/α.因此，加权POCS方法是解无约束优化问题的方法，权重因子可以看作为拟合误差项的系数. 与之不同的是，原始POCS方法实际上是解等式约束优化问题的方法， 这是两者最大的不同之处，关于原始POCS方法的理论推导将在其他文章中阐述.

由公式(5)可知，POCS方法得到的解uN是dobs+(I-Φ)dN，其中dN=Ψ(xN)是时间-空间域的迭代解. 这里的(I-Φ)dN的作用是将采样点位置的元素置为0，而未被采样的点的位置的元素置为dN，则uN的采样点位置的解为观测数据dobs，未被采样的位置的解为dN.如果观测数据dobs含有噪声， 这样就会将噪声带入插值结果，因此POCS方法对于含噪声数据是不合适的.加权POCS算法的解在采样位置为αdobs+(1-α)dN，在未被采样的位置为dN，采样位置的解仍然受到含噪声的观测数据dobs的影响，此外该方法中除了需要认真选择最小阈值外，还需要权重因子α来得到好的结果.由公式(6)可知，实际上权重因子可以消去，只要选择合适的λ或最小阈值就可以实现同时插值和去噪，如果最小阈值选择合适，并且返回dN为最终解，则可以实现同时插值和去噪，因此我们得到一种改进的POCS算法.该算法的迭代公式为

dk+1=ΨHτk(ΨT(dobs+(I-Φ)dk)),k=1,2,…,N，

(12)其中的迭代解为dk，属于时间-空间域的解，当多次迭代后，返回dN为最终解.下面我们称原始POCS为POCS(uN)，改进的POCS方法为POCS(dN)，对这两个算法进行比较，验证改进算法对含噪声数据和不含噪声数据插值的结果.

4 数值模拟

本节前两个例子用来证明改进POCS 方法对无噪声数据的插值效果；后两个例子用来证明该方法对含噪声数据插值的效果.我们用先插值后去噪的结果来验证同时插值和去噪的可靠性.本文采用 curvelet 变换为稀疏变换(Candès et al.，2006).

4.1 无噪声模拟数据的插值

图1 (a)原始模拟数据；(b)随机采样的模拟数据Fig.1 The original synthetic data (a) and randomly sampled synthetic data (b)

图2 (a) POCS (uN) 方法插值结果；(b) 图2a与原始数据的误差Fig.2 Interpolation of POCS (uN) method (a) and difference between Fig.2a and the original data (b)

图3 (a) POCS(dN) 方法的插值结果；(b) 图3a与原始数据的误差Fig.3 Interpolation of POCS(dN) method (a) and difference between Fig.3a and the original data (b)

图4 (a) 原始炮数据； (b) 随机采样的炮数据Fig.4 The original field data (a) and randomly sampled field data (b)

图5 (a) POCS(uN) 插值结果； (b) 图5a与原始数据的误差Fig.5 Interpolation of POCS (uN) method (a) and difference between Fig.5a and the original data (b)

图6 (a) POCS(dN) 方法的插值结果；(b) 图6a与原始数据的误差Fig.6 Interpolation of POCS (dN) method (a) and difference between Fig 6a and the original data (b)

图7 (a) 噪声模拟数据； (b) 含噪声采样数据Fig.7 The noisy synthetic data (a) and randomly sampled noisy synthetic data (b)

图8 (a) POCS(uN)方法插值结果； (b) 图8a与图1a的误差Fig.8 Interpolation of POCS(uN) method (a) and difference between Fig.8a and Fig.1a (b)

图9 (a) POCS(uN)插值结果的第3道； (b) POCS(uN)插值结果的第5道Fig.9 The 3rd trace of POCS(uN) interpolation (a) and the 5th trace of POCS (uN) interpolation (b)

图10 (a) POCS (dN) 方法的插值结果； (b) 图10a与原始数据的误差Fig.10 Interpolation of POCS (dN) method (a) and difference between Fig.10a and Fig.1a (b)

图11 (a) POCS (uN)插值结果，最小阈值为 0.01； (b) 图11a去噪结果Fig.11 Interpolation of POCS (uN) with least threshold 0.01 (a) and denoising result of Fig.11a (b)

图12 (a) 含噪声炮数据， 信噪比为5.1492 dB； (b) 采样的数据， 信噪比为2.34 dBFig.12 The noisy field data， SNR=5.1492 dB (a) and randomly sampled noisy field data， SNR=2.34 dB (b)

图13 (a) POCS(uN)方法插值结果； (b) 图13a的残差Fig.13 Interpolation of POCS (uN) method (a) and residual of Fig.13a (b)

图14 (a) POCS(dN)方法的插值结果； (b) 图14a的残差Fig.14 Interpolation of POCS(dN) method (a) and residual of Fig.14a (b)

图15 (a) POCS(uN)插值结果，最小阈值1； (b) 图15a的去噪结果Fig.15 Interpolation of POCS (uN) method with least threshold 1 (a) and denoising result of Fig.15a (b)

表1 无噪声模拟数据POCS(uN)方法和POCS(dN)方法的插值结果比较Table 1 Interpolation of POCS (uN) and POCS(dN) for noiseless synthetic data

这里POCS(uN)方法比POCS(dN)方法的结果稍微好一点，这是因为在采样位置，POCS(uN)的插值结果被采样dobs代替了，但是这仅仅对无噪声数据是可行的.对于含噪声数据这样得到的解是不可接受的，下面的含噪声数据插值的例子证明了这一点.4.3 含噪声模拟数据的插值图7a是含有零均值随机噪声的模拟数据，该数据的信噪比为4.1617 dB，图7b是图7a的采样结果，信噪比为1.6023 dB，采样道数为128.两种算法的最小阈值为0.17，迭代步数N为 60. 图8a是POCS (uN)的插值结果，由于仍然含有较严重噪声，因此图8a的结果不能令人满意，图8b是图8a与原始数据(图 1a)的误差.图9a给出了POCS(uN)插值结果的第3道，该道为采样道，图9b是POCS(uN)插值的第5道，这是一个未被采样的道，可以看到采样道的插值结果不是很好. 图10a是POCS(dN) 的插值结果，图10b是图10a与原始数据的误差. 两种算法的CPU时间，信噪比和相对误差见表3.

表2 无噪声真实数据POCS(uN)方法和POCS(dN)方法插值结果比较Table 2 Interpolation of POCS(uN) and POCS(dN) for noiseless field data

为了验证同时插值和去噪的可靠性，我们首先用POCS(uN)方法对图7b进行插值，然后对插值结果去噪. POCS(uN)方法的最小阈值为0.01，最大迭代步数N为60. 插值结果见图11a，信噪比为5.0561 dB.对图11a采用阈值去噪的结果见图11b，信噪比为12.6132 dB.很明显，图10a和图11b的信噪比相当，因此当参数选择合适时，同时插值和去噪是稳定的.

表3 含噪声合成数据POCS (uN)方法和POCS(dN)方法插值结果比较Table 3 Interpolation of POCS(uN) and POCS(dN) for noisy synthetic data

4.4 含噪声真实地震数据的插值

图12a是图4a添加随机白噪声的结果，信噪比为5.1492 dB，图12b是图12a随机采样的结果，信噪比为2.34 dB. 选择的最小阈值为16，最大迭代步数N为 60. 图13a是POCS(uN)的插值结果，图13b是图13a的残差.图14a是 POCS(dN)的插值结果，图14b是图14a的残差. 两种算法的CPU 时间，信噪比SNR和相对误差见表4. 由图13a可知，POCS(uN)方法的插值结果是不可接受的，但是POCS(dN)方法可以得到较满意的结果. 为了验证同时插值和去噪的可靠性，图15a是用POCS(uN)方法的插值结果，所选的最小阈值为1，信噪比为4.9458 dB，图15b是图15a阈值法的去噪结果，信噪比为8.9228 dB.上述结果存在一些噪声，这是因为有效信号的部分curvelet系数和噪声的curvelet系数耦合在一起了，仅仅通过阈值类去噪算法无法消除这种现象.

根据以上的数值模拟，POCS(dN)方法的计算效率和结果得到了验证，其不仅仅适合去无噪声数据的插值，而且适合含噪声数据的插值.

表4 含噪声真实数据POCS (uN)方法和POCS(dN) 方法插值结果比较Table 4 Interpolation of POCS (uN) and POCS(dN) for noisy field data

5 结论

本文由迭代阈值类方法推导出了加权POCS方法，证明其是解无约束优化问题的一种方法，权重因子可以看作拟合误差项的系数，同时也放大了最小阈值. 本文还提出了一种改进的 POCS方法来实现同时插值和去噪，该方法只含有一个超参数(阈值)，通过选择合适的阈值就能够实现插值和去噪. 数值模拟表明，新的方法比原始POCS方法不增加任何计算量，在保持高效计算的同时能够很好地插值无噪声数据和含噪声数据. 通过和先插值后去噪结果的比较，证明了同时插值和去噪的可靠性.

基于迭代阈值类算法进行同时插值和去噪的关键是选择合适的正则参数λ.λ和噪声的能量有直接关系，虽然可以针对每个数据的噪声能量来估计λ，但是这样做效率较低，另外本文采用最大迭代步数来控制收敛，这个迭代步数是根据我们的经验选择的，更实用的方法是选择一种自适应的停机准则来控制迭代停止.为了便于我们的推导，我们采用二维数据进行验证，本文的算法和结果对高维数据同样适用.本算法不仅仅适用于curvelet变换，其他的稀疏变换也适用于该算法框架.未来的研究方向之一是实现含噪声高维地震数据的插值和基于其他小波类变换的波场插值.关于POCS方法是解等式约束优化问题的方法的证明将在其他文章中阐述.

Abma R, Kabir N. 2006. 3D interpolation of irregular data with a POCS algorithm.Geophysics, 71(6): E91-E97.Berkhout A J, Verschuur D J. 1997. Estimation of multiple scattering by iterative inversion, Part I: Theoretical considerations.Geophysics, 62(5): 1586-1595.Blumensath T, Davies M E. 2008. Iterative thresholding for sparse approximations.JournalofFourierAnalysisandApplications, 14(5-6): 629-654. Bunks C, Saleck F M, Zaleski S, et al. 1995. Multiscale seismic waveform inversion.Geophysics, 60(5): 1457-1473.

Candès E, Demanet L, Donoho D, et al. 2006. Fast discrete curvelet transforms.SIAMJournalonMultiscaleModelingandSimulation, 5(3): 861-899.

Cao J J, Wang Y F, Zhao J T, et al. 2011. A review on restoration of seismic wavefields based on regularization and compressive sensing.InverseProblemsinScienceandEngineering, 19(5): 679-704.

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.ChineseJournalofGeophysics(in Chinese), 55(2): 596-607, doi: 10.6038/j.issn.0001 5733.2012.02.22.Daubechies I, Defrise M, de Mol C. 2004. An iterative thresholding algorithm for linear inverse problems with a sparsity constraint.CommunicationsonPureandAppliedMathematics, 57(11): 1413-1457.

Fomel S, Liu Y. 2010. Seislet transform and seislet frame.Geophysics, 75(3): V25-V38.

Gao J J, Stanton A, Naghizadeh M, et al. 2013. Convergence improvement and noise attenuation considerations for beyond alias projection onto convex sets reconstruction.GeophysicalProspecting, 61(S1): 138-151.

Geng Y, Wu R S, Gao J H. 2012. Dreamlet compression of seismic data.ChineseJournalofGeophysics(in Chinese), 55(8): 2705-2715, doi: 10.6038/issn.0001-5733.2012.08.023.

Hennenfent G, Herrmann F J. 2006. Seismic denoising with nonuniformly sampled curvelets. Computing in Science and Engineering, 8: 16-25. Herrmann F J, Hennenfent G. 2008. Non-parametric seismic data recovery with curvelet frames.GeophysicalJournalInternational, 173(1): 233-248.Kreimer N, Sacchi M D. 2012. A tensor higher-order singular value decomposition for prestack seismic data noise reduction and interpolation.Geophysics, 77(3): V113-V122.

Liu B, Sacchi M D. 2004. Minimum weighted norm interpolation of seismic records.Geophysics, 69(6): 1560-1568.

Lu J M. 1993. Principle of Seismic Exploration (in Chinese). Dongying: Petroleum University Press. Ma J W. 2013. Three-dimensional irregular seismic data reconstruction via low-rank matrix completion.Geophysics, 78(5): V181-V192.Mansour H, Herrmann F J, Yilmaz O. 2013. Improved wavefield reconstruction from randomized sampling via weighted one-norm minimization.Geophysics, 78(5): V193-V206.

Oropeza V, Sacchi M D. 2011. Simultaneous seismic data denoising and reconstruction via multichannel singular spectrum analysis.Geophysics, 76(3): V25-V32.

Ronen J. 1987. Wave-equation trace interpolation.Geophysics, 52(7): 973-984.

Spitz S. 1991. Seismic trace interpolation in the f-x domain.Geophysics, 56(6): 785-794. Stolt R H. 2002. Seismic data mapping and reconstruction.Geophysics, 67(3): 890-908. Symes W W. 2007. Reverse time migration with optimal checkpointing.Geophysics, 72(5): SM213-SM221. Trad D O, Ulrych T J, Sacchi M D. 2002. Accurate interpolation with high-resolution time-variant Radon transforms.Geophysics, 67(2): 644-656. Wang Y F. 2007. Computational Methods for Inverse Problems and Their Applications (in Chinese). Beijing: Higher Education Press. 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.GeophysicalJournalInternational, 187(1): 199-213.

Xu S, Zhang Y, Lambaré G. 2010. Antileakage Fourier transform for seismic data regularization in higher dimensions.Geophysics, 75(6): WB113-WB120.Yang P L, Gao J H, Chen W C. 2012. Curvelet-based POCS interpolation of nonuniformly sampled seismic records.JournalofAppliedGeophysics, 79: 90-99.

Yilmaz O. 1987. Seismic Data Processing, Investigations in Geophysics, Vol. 2. Tulsa, Okla: Society of Exploration Geophysicists. Zhang H, Chen X H. 2013. Seismic data reconstruction based on jittered sampling and curvelet transform.ChineseJournalofGeophysics(in Chinese), 56(5): 1637-1649, doi: 10.6038/cjg20130521.

Zwartjes P M. 2005. Fourier reconstruction with sparse inversion [Ph. D. Thesis]. Delft: Delft University of Technology.

附中文参考文献

曹静杰, 王彦飞, 杨长春. 2012. 地震数据压缩重构的正则化与零范数稀疏最优化方法. 地球物理学报, 55(2): 596-607, doi: 10.6038/j.issn.0001 5733.2012.02.22.

耿瑜, 吴如山, 高静怀. 2012. 基于Dreamlet变换的地震数据压缩理论与方法. 地球物理学报, 55(8): 2705-2715, doi: 10.6038/issn.0001-5733.2012.08.023.

陆基孟. 1993. 地震勘探原理. 东营: 石油大学出版社.

王彦飞. 2007. 反演问题的计算方法及其应用. 北京: 高等教育出版社.

张华, 陈小宏. 2013. 基于jitter采样和曲波变换的三维地震数据重建. 地球物理学报, 56(5): 1637-1649, doi: 10.6038/cjg20130521.

(本文编辑 胡素芳)

An improved projection onto convex sets method for simultaneous interpolation and denoising

CAO Jing-Jie1， WANG Ben-Feng2

1ShijiazhuangUniversityofEconomics，Shijiazhuang050031，China2StateKeyLaboratoryofPetroleumResourcesandProspecting，ChinaUniversityofPetroleum，Beijing102249，China

Seismic interpolation can provide a complete wavefield from incomplete wavefield information and is crucial for some processings such as multiple elimination, wave equation migration and denoising. Currently, most interpolation methods are suitable for high signal-to-noise ratio data. However, filed seismic data often contain different kinds of noise, especially the random noise. The existence of noise can deteriorate the result of interpolation. The weighted projection onto convex sets(POCS) method is a method for random noise contained data interpolation. However, there is no specific explanation of this method except the iteration formula. This paper analyzed the algorithm theory of the weighted POCS method, and proved that it is a method for unconstrained optimization problem, the weighting factor can be regarded as the coefficient of the fitting error. Meanwhile, we also proposed an improved POCS method for simultaneous interpolation and denoising. This algorithm requires no additional computational effort compared with the original POCS method.Seismic interpolation especially the high dimensional interpolation belongs to large-scale computational problems, thus fast algorithms should be chosen to improve the efficiency of seismic interpolation. Sparse transform methods change the interpolation problem into a sparse optimization based on the assumption that seismic data can be expressed sparsely in some transform domain. POCS method is such a method with high efficiency and robust results. However, this method can not get acceptable results for random noise contained data interpolation. The weighted POCS method is a method for random noise contained data interpolation on the basis of the POCS method. But there is no specific explanation of this method up to now. This paper deduced the weighted POCS method from the iterative thresholding method, we proved that the weighted POCS method is a method for unconstrained optimization, the weighted factor can be seen as the coefficient of the fitting error. However, the original POCS method is aiming at an equality constrained optimization. This is the essential difference between them. Meanwhile, we also proposed an improved POCS method with the same computational complexity as the original POCS method. The improved POCS method can realize interpolation and denoising simultaneously by adjustment of the least thresholds.In order to test the interpolation and denoising capability of the improved POCS method, four data sets were used to make numerical experiments. The first experiment is a synthetic noiseless data interpolation. The dimension of the synthetic data is 256×256 with trace distance 25 m and time sampling rate 4 ms. The original POCS method is chosen as the benchmark. Both the computational efficiency and numerical result are comparable for them. The second data is field data with 115 traces. The trace distance is 25 m and the time sampling rate is 4 ms. This experiment also proved that both of them can get comparable results for noiseless data interpolation. The third example is a noise contained data interpolation. The original POCS method can not get acceptable result no matter how to choose the least thresholds, while the improved POCS method can provide good result by choosing the least thresholds carefully. The fourth example also proved that the original POCS method failed in the noise contained field data interpolation, and the improved POCS method is suitable for noise contained data interpolation. Numerical experiment on first-interpolation-then-denoising proved the rationality and robustness of simultaneous interpolation and denoising.We deduced the weighted POCS method from the iterative thresholding method, and proved that it is a method for unconstrained optimization, the weighting factor is the coefficient of the fitting error. However, the original POCS method is a method for equality constrained optimization. This is the essential difference between them. We also proposed an improved POCS method to realize simultaneous interpolation and denoising. This method just need to adjust the least thresholds to realize noiseless data and random noise contained data interpolation. The computational effort of the improved POCS method is the same as the original POCS method but with more broad applications.

Sparse inversion； Denoise； Interpolation； Inversion

国家自然科学基金项目(41204075)，河北省自然科学基金(D2014403007)和河北省普通高等学校青年拔尖人才项目(BJ2014049)共同资助．

曹静杰，男，副教授，1982年生，2011年于中国科学院地质与地球物理研究所获博士学位，主要从事地震信号处理和地球物理反演算法研究. E-mail：cao18601861@163.com

10.6038/cjg20150826.

10.6038/cjg20150826

P631

2014-09-17，2015-05-26收修定稿

曹静杰，王本锋. 2015. 基于一种改进凸集投影方法的地震数据同时插值和去噪.地球物理学报，58(8):2935-2947,

Cao J J, Wang B F. 2015. An improved projection onto convex sets method for simultaneous interpolation and denoising.ChineseJ.Geophys. (in Chinese),58(8):2935-2947,doi:10.6038/cjg20150826.