2. 中国石油大学(北京)油气资源与探测国家重点实验室, 北京 102249
2. State Key Laboratory of Petroleum Resources and Prospecting, China University of Petroleum, Beijing 102249, China
因障碍物、禁采区及海洋拖缆羽状漂移等各种因素的影响, 采集的地震数据沿空间方向不规则, 该不规则性会影响自由表面多次波去除、AVO分析、反演及偏移等处理结果的精度, 因此有必要对采集到的地震数据进行插值重建[1]。插值重建方法可以大致归纳为4类[2-4]:基于稀疏变换的插值重建方法[5]; 基于预测滤波的插值重建方法[6]; 基于波动方程的插值方法[7]和基于降秩的插值方法[8-10]。
随着压缩感知理论的问世[11], 基于稀疏变换的插值重建方法, 特别是基于稀疏Curvelet变换的插值方法得到广泛关注[12-13]。高精度的凸集投影(POCS)方法[2-3, 14-15]是基于稀疏变换的重要插值方法之一, 该方法常用于图像重构[16], 1991年, MENKE[17]首次将其引入地球物理领域的插值重建方面。目前, 国内专家学者研究了基于压缩感知策略的地震数据插值重建问题[18-21]。大多数插值重建方法属于迭代方法的范畴[12-13, 22-23], 因此建立一个合理的质量控制准则(QCC)成为高效地震数据插值重建的关键, 质量控制准则可在保证插值精度的同时, 适时地结束迭代, 达到提高插值重建计算效率的目的。GAO等[24]给出了两种不同的质量控制准则来监测插值迭代的收敛过程:第一种方法利用了完整的数据体, 仅适用于插值重建理论方法研究; 第二种方法利用了相邻的迭代更新解, 但容易陷入局部极值, 影响插值重建的计算精度[25]。
本文基于迭代硬阈值方法及投影算子推导了精确高效POCS方法, 在此基础上, 提出了一种新的质量控制准则, 该质量控制准则仅利用了观测数据的信息, 具有更好的适应性。基于稀疏Curvelet变换及POCS插值重建方法, 利用两套模拟数据和实际数据验证了该质量控制准则的有效性及计算的高效性。
1 方法原理 1.1 POCS插值方法空间方向不规则的观测地震数据dobs可以近似表征为:
| $ {\mathit{\boldsymbol{d}}_{\rm{obs}}} = {\rm{ }}\mathit{\boldsymbol{Rd}}{_0} $ | (1) |
式中:R为对角型采样算子; d0为需重建的完整数据体。由于受到观测地震数据频带有限及数据缺失等因素的影响, 方程(1)的求解是不适定的, 因此基于压缩感知理论, 利用稀疏变换及稀疏约束构建目标泛函为:
| $ \mathit{\Phi} \left( {{\rm{ }}\mathit{\boldsymbol{x}}{\rm{ }}} \right) = \left\| {\mathit{\boldsymbol{d}}{_{\rm{obs}}}{\rm{-}}R{C^{\rm{T}}}x} \right\|_2^2 + \lambda {\left\| \mathit{\boldsymbol{x}} \right\|_0} $ | (2) |
式中:x为Curvelet系数向量; CT为Curvelet逆变换(C为Curvelet变换); λ为正则化因子, 用于权衡Curvelet系数稀疏性与拟合残差。方程(2)可以由迭代硬阈值方法进行迭代求解[26-27]:
| $ \mathit{\boldsymbol{x}}{_{k + 1}} = {T_{{\lambda _k}}}[{\rm{ }}\mathit{\boldsymbol{x}}{_k} + {({\rm{ }}\mathit{\boldsymbol{RC}}{^{\rm{T}}})^{\rm{T}}}({\rm{ }}\mathit{\boldsymbol{d}}{_{\rm{obs}}}-{\rm{ }}\mathit{\boldsymbol{RC}}{^{\rm{T}}}){\rm{ }}\mathit{\boldsymbol{x}}{_k}] $ | (3) |
式中:xk为Curvelet域内第k次的更新解; Tλk为硬阈值函数; λk为阈值[3]。地震数据插值重建的目的是得到数据域重建结果, 因此将变换域的更新结果转换到数据域得:
| $ \mathit{\boldsymbol{d}}{_{k + 1}} = {\rm{ }}\mathit{\boldsymbol{C}}{^{\rm{T}}}\mathit{\boldsymbol{x}}{_{k + 1}} $ | (4) |
插值迭代收敛之前, 在观测位置处dk+1的值偏离观测数据dobs, 为了提高插值重建的效率, 将dk+1投影到观测平面{d|dobs=Rd}上, 得:
| $ \begin{array}{l} {{\mathit{\boldsymbol{\bar d}}}_{k + 1}} = {\rm{ }}\mathit{\boldsymbol{Rd}}{_{{\rm{obs}}}} + \left( {{\rm{ }}\mathit{\boldsymbol{I}}{\rm{- }}\mathit{\boldsymbol{R}}{\rm{ }}} \right){\rm{ }}\mathit{\boldsymbol{d}}{_{k + 1}} = {\rm{ }}\mathit{\boldsymbol{d}}{_{{\rm{obs}}}} + \left( {{\rm{ }}\mathit{\boldsymbol{I}}{\rm{- }}\mathit{\boldsymbol{R}}{\rm{ }}} \right)\cdot\\ \;\;\;\;\;\;\;\;{\mathit{\boldsymbol{C}}^{\rm{T}}}{T_{{\lambda _k}}}\{ {\rm{ }}\mathit{\boldsymbol{C}}[{\rm{ }}\mathit{\boldsymbol{d}}{_{\rm{obs}}} + \left( {{\rm{ }}\mathit{\boldsymbol{I}}{\rm{-}}\mathit{\boldsymbol{R}}{\rm{ }}} \right){\rm{ }}\mathit{\boldsymbol{d}}{_k})]\} \end{array} $ | (5) |
其中dk+1为观测数据置入之后的重建结果。记dk=dobs+(I-R)dk, 则POCS公式可以归纳总结为:
| $ \begin{array}{l} {\mathit{\boldsymbol{d}}_{k + 1}} = {\rm{ }}{\mathit{\boldsymbol{C}}^{\rm{T}}}{T_{{\lambda _k}}}({\rm{ }}\mathit{\boldsymbol{C\bar d}}{_k})\\ {{\mathit{\boldsymbol{\bar d}}}_{k + 1}} = {\rm{ }}\mathit{\boldsymbol{d}}{_{\rm{obs}}} + \left( {{\rm{ }}\mathit{\boldsymbol{I}}{\rm{-}}\mathit{\boldsymbol{R}}{\rm{ }}} \right){\rm{ }}\mathit{\boldsymbol{d}}{_{k + 1}} \end{array} $ | (6) |
实际上POCS插值方法是一个两步法:稀疏变换域内更新解以及数据域内置入观测数据。该方法具有较高的插值重建精度及计算效率。
1.2 新型质量控制准则由于POCS插值方法是一种迭代方法, 有必要建立一个合理的质量控制准则来适时地终止迭代, 进而提高插值重建的计算效率。GAO等[24]给出了两种质量控制准则:
| $ J_1^k\mathit{\boldsymbol{ = }}{\left\| {\mathit{\boldsymbol{\bar d}}{_k}-{\rm{ }}\mathit{\boldsymbol{d}}{_0}} \right\|^2}/{\left\| {{\mathit{\boldsymbol{d}}_0}} \right\|^2} $ | (7) |
| $ J_2^k\mathit{\boldsymbol{ = }}{\left\| {\mathit{\boldsymbol{\bar d}}{_k}-{\rm{ }}{{\mathit{\boldsymbol{\bar d}}}_{k{\rm{-}}1}}} \right\|^2}/{\left\| {{{\mathit{\boldsymbol{\bar d}}}_k}} \right\|^2} $ | (8) |
为了更好地阐释每一种质量控制准则的物理意义, 作如图 1所示的示意图。完整数据体的位置集合Ω可以分解为两部分:观测位置Ω1及缺失位置Ω2, 满足:Ω1+Ω2=Ω。插值重建迭代收敛之前, dk在观测位置Ω1和缺失位置Ω2处, 均偏离完整数据体d0。然而, 由于观测数据置入的影响, 在观测位置Ω1处dk与完整数据体d0一致, 差异仅在于缺失位置Ω2处。方程(7)中, J1表征利用d0进行归一化的缺失位置Ω2处dk和d0之间的差异, 由于利用了完整数据体d0的信息, J1仅适用于理论插值重建方法研究。另外, J1可视为参考解来评价其它的质量控制准则。公式(8)中, J2表征利用dk进行归一化的缺失位置Ω2处dk和dk-1之间的差异, J2可能会导致插值重建陷入局部极值。上述两种质量控制准则均由缺失位置Ω2处的数据残差来定义, 然而利用观测位置Ω1处的数据信息可能具有更好的适应性。基于此, 本文研究了基于POCS方法的新型质量控制准则, 仅利用观测位置Ω1处的信息, 具有更好的适应性, 且与参考值J1有相似的效果:
| $ J_3^k\mathit{\boldsymbol{ = }}{\left\| {{\mathit{\boldsymbol{d}}_k}-{\rm{ }}{{\mathit{\boldsymbol{\bar d}}}_k}} \right\|^2}/{\left\| {{\mathit{\boldsymbol{d}}_{{\rm{obs}}}}} \right\|^2} $ | (9) |
|
图 1 地震数据位置集合分解示意 |
式中:J3表征利用观测数据dobs进行归一化的观测位置Ω1处dk和dk之间的差异; dk为观测数据置入之前数据域第k次更新解; dk为观测数据置入之后的数据域第k次更新解。dk和dk在缺失位置Ω2处一致, 差异在于观测位置Ω1处(观测位置Ω1处, dk为观测数据真实值, dk为重构值)。J3与J1相似, 不同点在于前者仅利用了观测数据及观测数据置入前后的更新解, 在实际应用中具有更好的灵活性。利用提出的新型质量控制准则, 可以适时地结束插值迭代过程, 在保证插值精度的同时, 进一步提高插值重建的计算效率。
2 数值实验首先, 利用两套模拟数据验证本文提出的质量控制准则的可行性及有效性; 其次, 利用一套实际资料进一步验证本文提出的质量控制准则能够在保持插值精度的同时, 提高插值计算效率的有效性。
2.1 模拟算例图 2a为层状介质模型对应的模拟数据体(模拟数据1), 包含201道, 每道1 001个采样点, 空间采样率与时间采样率分别为12.5 m和2 ms。图 2b为利用改进的Jittered欠采样方法进行50%地震道缺失得到的不规则采样数据。采用基于稀疏Curvelet变换的POCS插值方法进行插值重建, 根据经验, 最大迭代次数设置为50次, 利用不同的质量控制准则J1, J2, J3, 并进行归一化, 得到的收敛曲线如图 3所示。J1为理论的质量控制准则, 可以作为参考解, 合理且逐渐下降至收敛, 其收敛性可以保证插值方法具有较高的精度。J2为利用相邻迭代更新解定义的质量控制准则, 较大地偏离参考解, 可能使插值迭代陷入局部极小值, 进而影响插值重建结果的精度。然而, J3利用观测位置Ω1处的信息定义质量控制准则, 其与参考解J1具有较好的一致性, 具有更广的实际应用前景。图 3验证了质量控制准则J3的可行性与有效性, 即当J3满足预先设定的允许误差时, 插值迭代可以适时地终止, 在保证插值精度的同时提高插值计算效率。
|
图 2 完整数据体(a)与不规则采样数据体(b) |
|
图 3 不同质量控制准则对应的归一化收敛曲线(模拟数据1) |
图 4a为利用基于稀疏Curvelet变换的POCS插值重建方法得到的插值结果, 可以看出, 其与原始数据具有较高的一致性。图 4b为插值结果与原始数据之间的重构误差, 可见, 图上除了有微弱的边界效应外, 残差较小。该边界效应是由边界处数据的不规则性产生, 可以通过扩边的方式降低这种边界效应, 但会增加计算量, 具体量化分析留待后续研究。图 3和图 4验证了基于稀疏变换的POCS插值方法的有效性, 且提出的质量控制准则与参考值的一致性较好。
|
图 4 插值重建结果(a)与重建误差(b) |
为了进一步验证本文提出的质量控制准则的有效性, 将其应用于含有干涉及弱反射的模拟数据(模拟数据2)中, 图 5a和图 5b分别为完整数据体及不规则采样数据体, 其中图 5b缺失40%数据。利用不同的质量控制准则得到的归一化收敛曲线如图 6所示, 其中最大迭代次数为50次。图 6与图 3相似, 均显示出本文提出的质量控制准则J3与理论解J1具有较好的一致性, 但J2与理论解J1差别较大。图 6进一步验证了本文提出的质量控制准则J3的可行性与有效性, 其与理论值J1具有较好的一致性, 同时也指出了J2的缺点。插值重建结果如图 7a所示, 其与原始完整数据体具有较好的一致性, 且残差较小, 可以忽略不计(图 7b), 插值重建结果与插值残差验证了基于稀疏变换的POCS插值方法的有效性。
|
图 5 完整数据体(a)与40%缺失不规则采样数据体(b) |
|
图 6 不同质量控制准则对应的归一化收敛曲线(模拟数据2) |
|
图 7 插值结果(a)与插值误差(b) |
上述两个模拟算例验证了本文提出的质量控制准则J3的有效性与可行性, 其与理论值J1具有较好的一致性, 但又与J1不同, J3仅利用了观测数据的信息, 在实际数据插值重建中具有更好的灵活性, 应用范围更广泛, 在保证插值重建精度的同时, 提高了插值重建计算效率。
2.2 实际数据算例对不规则实际数据进行插值重建, 理论质量控制准则J1不适用, 因为J1利用了完整数据体的信息, 但是本文提出的质量控制准则J3是一个较好的选择, 而且与J1具有相似的性态。对于不规则缺失实际数据插值重建, 下面仅比较质量控制准则J2和J3的应用效果, 进一步验证质量控制准则J3的优势。
图 8a为利用改进的Jittered欠采样方法得到的不规则实际数据, 数据的不规则性会影响自由表面多次波去除、偏移与反演的精度。为了插值重建完整的数据体, 利用基于稀疏Curvelet变换的POCS插值方法得到的重构数据体如图 8b所示, 数据插值中最大迭代次数为50。由图 8b可看出, 插值重建后地震数据的横向连续性得到提高, 该重建的地震数据可以作为后续地震数据处理的完整数据体。质量控制准则J2和J3对应的归一化收敛曲线如图 8c所示, 可以看出, J3对应的收敛曲线缓慢下降, 具有一定的合理性, 然而J2可能会使插值迭代陷入局部极值。为了利用J3进一步提高插值重建的计算效率, 当满足预先设定的容许残差(J3＜0.005)时, 插值迭代被适时终止, 插值结果如图 8d所示。此时, 迭代次数为30次(如图 8c中纵向虚线所示), 可以看出, 相对50次迭代而言, 计算效率提高了40%, 但没有明显的降低插值重建精度。实际上, 图 8b与图 8d具有较好的一致性, 从而验证了本文提出的新型质量控制准则J3的有效性, 在保证插值重建精度的同时, 可以有效地提高插值重建的计算效率。从图 8b和8d中还可以看出, 插值结果中有方向各异的微弱假象, 该假象由稀疏变换域系数截断引起, 因此需研究更多的插值重建策略, 进一步衰减假象、提高插值重建的精度。
|
图 8 实际数据插值重建结果 a 不规则实际数据; b 迭代50次后插值重建结果; c 质量控制准则J2和J3对应的归一化收敛曲线; d 当J3＜0.005时的插值重建结果 |
本文基于POCS插值重建方法推导过程, 提出了一种新型质量控制准则J3来监测插值重建过程, 可以适时地终止插值迭代过程, 在保证插值重建精度的同时, 进一步提高插值重建的计算效率。阐述并对比分析了不同类型的质量控制准则, 提出的质量控制准则J3与理论值J1具有较好的一致性, 而J3仅利用了观测数据的信息, 在实际地震数据插值重建监测过程中具有更好的灵活性与适用性, 应用范围更广。模拟与实际数据处理验证了新型质量控制准则J3的有效性。
| [1] |
GAN S, WANG S, CHEN Y, et al. Dealiased seismic data interpolation using seislet transform with low-frequency constraint[J]. IEEE Geoscience and Remote Sensing Letters, 2015, 12(10): 2150-2154. DOI:10.1109/LGRS.2015.2453119 |
| [2] |
WANG B F, WU R S, GENG Y, et al. Dreamlet-based interpolation using POCS method[J]. Journal of Applied Geophysics, 2014, 109(10): 256-265. |
| [3] |
WANG B F, WU R S, CHEN X, et al. Simultaneous seismic data interpolation and denoising with a new adaptive method based on dreamlet transform[J]. Geophysical Journal International, 2015, 201(2): 1180-1192. |
| [4] |
WANG B F, LI J, CHEN X. A novel method for simultaneous seismic data interpolation and noise removal based on the L0 norm constraint[J]. Journal of Seismic Exploration, 2015, 24(2): 187-204. |
| [5] |
NAGHIZADEH M, SACCHI M D. Beyond alias hierarchical scale curvelet interpolation of regularly and irregularly sampled seismic data[J]. Geophysics, 2010, 75(6): WB189-WB202. DOI:10.1190/1.3509468 |
| [6] |
SPITZ S. Seismic trace interpolation in the F-X domain[J]. Geophysics, 1991, 56(6): 785-794. DOI:10.1190/1.1443096 |
| [7] |
RONEN J. Wave-equation trace interpolation[J]. Geophysics, 1987, 52(7): 973-984. DOI:10.1190/1.1442366 |
| [8] |
MA J. Three-dimensional irregular seismic data reconstruction via low-rank matrix completion[J]. Geophysics, 2013, 78(5): V181-V192. DOI:10.1190/geo2012-0465.1 |
| [9] |
GAO J, SACCHI M D, CHEN X. A fast reduced-rank interpolation method for prestack seismic volumes that depend on four spatial dimensions[J]. Geophysics, 2013, 78(1): V21-V30. |
| [10] |
GAO J, STANTON A, SACCHI M D. Parallel matrix factorization algorithm and its application to 5D seismic reconstruction and denoising[J]. Geophysics, 2015, 80(6): V173-V187. DOI:10.1190/geo2014-0594.1 |
| [11] |
DONOHO D L. Compressed sensing[J]. IEEE Transactions on Information Theory, 2006, 52(4): 1289-1306. |
| [12] |
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 |
| [13] |
HERRMANN F J, WANG D, HENNENFENT G, et al. Curvelet-based seismic data processing:a multiscale and nonlinear approach[J]. Geophysics, 2008, 73(1): A1-A5. |
| [14] |
ABMA R, KABIR N. 3D interpolation of irregular data with a POCS algorithm[J]. Geophysics, 2006, 71(6): E91-E97. |
| [15] |
WANG B F, CHEN X H, LI J Y, 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 Observations and Remote Sensing, 2016, 9(1): 228-235. DOI:10.1109/JSTARS.2015.2496374 |
| [16] |
YOULA D C, WEBB H. Image restoration by the method of convex projections:part 1-theory[J]. IEEE Transactions on Medical Imaging, 1982, 1(2): 81-94. |
| [17] |
MENKE W. Applications of the POCS inversion method to interpolating topography and other geophysical fields[J]. Geophysical Research Letters, 1991, 18(3): 435-438. DOI:10.1029/90GL00343 |
| [18] |
冯飞, 王征, 刘成明, 等. 基于Shearlet变换稀疏约束地震数据重建[J]. 石油物探, 2016, 55(5): 682-691. FENG F, WANG Z, LIU C M, et al. Seismic data reconstruction based on sparse constraint in the Shearlet domain[J]. Geophysical Prospecting for Petroleum, 2016, 55(5): 682-691. DOI:10.3969/j.issn.1000-1441.2016.05.007 |
| [19] |
马继涛, 王建花, 刘国昌. 基于频率域奇异值分解的地震数据插值去噪方法研究[J]. 石油物探, 2016, 55(2): 205-213. MA J T, WANG J H, LIU G C. Seismic data noise attenuation and interpolation using singular value decomposition in frequency domain[J]. Geophysical Prospecting for Petroleum, 2016, 55(2): 205-213. DOI:10.3969/j.issn.1000-1441.2016.02.006 |
| [20] |
陈小春, 陈辉, 喻勤, 等. 反假频POCS数据规则化及其在偏移成像中的应用[J]. 石油地球物理勘探, 2017, 52(1): 13-19. CHEN X C, CHEN H, YU Q, et al. Seismic data interpolation with anti-aliasing POCS method and its application in seismic migration imaging[J]. Oil Geophysical Prospecting, 2017, 52(1): 13-19. |
| [21] |
王本锋, 陈小宏, 李景叶, 等. POCS联合改进的Jitter采样理论曲波域地震数据重建[J]. 石油地球物理勘探, 2015, 50(1): 20-28. WANG B F, CHEN X H, LI J Y, et al. Seismic data reconstruction based on POCS and improved Jittered sampling in the curvelet domain[J]. Oil Geophysical Prospecting, 2015, 50(1): 20-28. |
| [22] |
TROPP J A, GILBERT A C. Signal recovery from random measurements via orthogonal matching pursuit[J]. IEEE Transactions on Information Theory, 2007, 53(12): 4655-4666. DOI:10.1109/TIT.2007.909108 |
| [23] |
YI L, TAKAHASHI K, SATO M. A fast iterative interpolation method in f-k domain for 3-D irregularly sampled GPR data[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remonte Sensing, 2016, 9(1): 9-17. DOI:10.1109/JSTARS.2015.2489260 |
| [24] |
GAO J, STANTON A, NAGHIZADEH M, et al. Convergence improvement and noise attenuation considerations for beyond alias projection onto convex sets reconstruction[J]. Geophysical Prospecting, 2012, 61(S1): 138-151. |
| [25] |
CHEN P, HUANG J, ZHANG X. A primal-dual fixed point algorithm for convex separable minimization with applications to image restoration[J]. Inverse Problems, 2013, 29(2): 025011. |
| [26] |
BLUMENSATH T, DAVIES M E. Iterative thresholding for sparse approximations[J]. Journal of Fourier Analysis and Applications, 2008, 14(5): 629-654. |
| [27] |
BLUMENSATH T, DAVIES M E. Iterative hard thresholding for compressed sensing[J]. Applied and Computational Harmonic Analysis, 2009, 27(3): 265-274. DOI:10.1016/j.acha.2009.04.002 |