地球物理学进展  2014, Vol. 29 Issue (4): 1752-1757   PDF    
引用本文
石颖, 张振, 李婷婷, 刘淑芬, 曾科. λ-f域加权抛物Radon变换地震数据重建方法研究[J]. 地球物理学进展, 2014, 29(4): 1752-1757, doi: 10.6038/pg20140435  
SHI Ying, ZHANG Zhen, LI Ting-ting, LIU Shu-fen, ZENG Ke. Investigation on seismic data reconstruction by λ-f domain weighted parabolic Radon transform method. Progress in Geophysics, 2014, 29(4): 1752-1757, doi: 10.6038/pg20140435   
λ-f域加权抛物Radon变换地震数据重建方法研究
石颖1,2, 张振1, 李婷婷1, 刘淑芬1, 曾科1     
1. 东北石油大学 地球科学学院, 大庆 163318;
2. 黑龙江省普通高校科技创新团队“断层变形、封闭性及与流体运移”, 大庆 163318
收稿日期: 2014-3-15 修回日期: 2014-6-11 .
基金项目：黑龙江省教育厅科学技术研究项目（12511025）资助.
作者简介：石颖，女，1976年生，博士，教授，博士生导师.主要从事地震资料处理方面的研究.（E-mail：shiyingdqpi@163.com）
摘要：大多数地震处理技术都要求地震数据具备完整性和规则性，然而，由于诸多因素的影响，地震勘探所采集到的资料普遍存在数据缺失，需要对地震数据进行重建.本文在传统Radon变换重建的基础上提出一种λ-f域加权抛物Radon变换的地震数据重建方法.通过引入新变量λ，消除了Radon变换算子对频率的依赖，使得Radon变换算子及算子的逆只需计算一次，显著提高了计算效率.同时，在λ-f域抛物Radon变换迭代计算过程中引入变化的权系数，更好地实现了λ-f域的能量聚焦.理论模型及实际数据试算表明，文中方法对地震数据重建精度较高，单道对比吻合较好.
关键词： λ-f域     抛物Radon变换     地震数据重建     权系数    
Investigation on seismic data reconstruction by λ-f domain weighted parabolic Radon transform method
SHI Ying1,2, ZHANG Zhen1, LI Ting-ting1, LIU Shu-fen1, ZENG Ke1    
1. Geoscience Institute of Northeast Petroleum University, Daqing 163318, China;
2. Science and Technology Innovation Team on Fault Deformation, Sealing and Fluid Migration, Daqing 163318, China
Abstract: Most seismic processing techniques require integrality and regularity of seismic data. However, as the seismic data often suffers from data missing during seismic exploration for the influence of many factors, we have to reconstruct the seismic data. In this paper, we propose to do seismic data reconstruction by λ-f domain weighted parabolic Radon transform method which is on basis of conventional Radon transform. The new variable λ is introduced to make the Radon transform operator frequency-independent. The transform operator and inverse operator could be computed only once, which can significantly improve the computational efficiency. During the iteration process of λ-f domain parabolic Radon transform, changing weighted coefficients are introduced to improve energy focusing problem in λ-f domain. Theoretical models and real data computation show that the accuracy of the seismic reconstruction could be improved with a better single-channel comparison by using the method in this paper.
Key words: λ-f domain     parabolic Radon transform     seismic data reconstruction     weighted coefficients    

0 引 言

由于采集环境、勘探成本等诸多因素的影响，地震数据在空间方向常常出现缺失、稀疏或不规则等情况，而稀疏或不规则采样的地震数据对多次波压制(刘伊克等，2004；李翔和胡天跃，2009；马继涛等，2009)及偏移成像(陈必远和马在田，1994;张文生等，2001；金胜汶等，2002；刘定进和印兴耀，2007；王西文等，2010)影响较大，也由此影响后续的地震解释精度，为此，高精度的地震数据重建(刘喜武等，2004；李信富和李小凡，2008；孟小红等，2008；高建军等，2011；陆艳洪等，2012；刘红艳等，2014)在数据处理中显得尤为重要.目前，基于变换域的数据重建方法主要包括Radon变换、Fourier变换和Curvelet变换等.Duijndam和Schonewille(1999)提出了利用非均匀快速傅里叶变换与最小范数准则相结合的方法进行数据规则化，随后，Hindriks和Duijndam(2000)将该方法拓展应用于3D地震数据处理.基于Fourier变换，Liu和Sacchi(2004)提出了一种的最小加权范数插值方法，将数据重建描述为最小二乘问题，但并没有给出较好的反假频方法.离散抛物Radon变换由Hampson(1986)提出，并在地震资料多次波压制、波场分离和数据重建领域等得到广泛应用，也有一些学者从分辨率，计算效率或反假频等多种角度出发，提出了时域、频域以及混合域的Radon变换计算方法，其中，Trad等(2002)提出了一种时变的高分辨率Radon域重建方法，但由于时变反演无法应用Levinson递推快速求解方法，极大地降低了计算效率.与时域方法需对大矩阵求逆的问题比较，频域方法可在地震信号有限带宽内对所有频率成分解耦地处理小矩阵的复杂计算问题.为此，Kabir和Verschuur(1995)利用抛物Radon变换实现了近偏移距缺失数据插值重建，但该方法迭代次数较多，计算效率低；Sacchi和Ulrych(1995)提出了频率域的高分辨率抛物Radon变换，但重建效果不理想；Schonewille和Duijndam(2001)使用最小二乘法在频域实现了地震数据的重建和规则化，并详细讨论了Radon域的采样准则及其计算效率问题；熊登等(2009)提出了混合域高分辨率抛物Radon变换的方法，对频率域抛物Radon变换引入时变的稀疏权，进而改善了高分辨Radon变换的计算效率.Abbad(2011)提出了一种改进的高效λ-f域抛物Radon变换，该方法显著提高了计算效率，但是并未在地震数据重建中取得应用.

本文提出一种基于λ-f域加权抛物Radon变换的地震数据规则化方法，该方法将加权思想(王维红等，2007)引入到λ-f域抛物Radon变换中.采用λ-f域抛物Radon变换消除Radon变换算子对频率的依赖，提高计算效率.通过引入变化的权系数对缺失道在数据域进行加权，增强变换域的聚焦性和成像精度.对理论模型及实际地震资料的CMP道集进行试算，重建数据效果较为理想.抽取单道数据对比分析表明重建数据的相位、振幅保持良好.

1 抛物Radon变换理论

抛物Radon正变换可表述为

式中，d(xn，t)为时空域数据，N为道数，m(q，τ)为Radon域数据，q为曲率，τ为Radon域零偏移距截距时间.

将(1)式进行傅里叶变换，其频域形式为

式中，M(q，f)和D(xn，f)分别为m(q，τ)和d(xn，t)的傅里叶变换，f为频率.抛物Radon正变换可在每个频率分量f上解耦完成.

由方程(2)，抛物Radon反变换的近似形式可表述为

式中，Nq为q值的个数.

方程(3)可写为矩阵形式

其中，

频率域抛物Radon变换的优点是频率域解耦，计算效率很高，但是在方程求解过程中对于信号频带内的每一个频率分量，都要计算一次Radon算子 L ni(f)及 L ni(f)的逆，增加了较大的计算量.

2 λ-f域抛物Radon变换理论

为了消除抛物Radon算子 L ni(f)对频率f的依赖性，引入新的变量λ=qf，则方程(3)式可写为

方程(8)可写成矩阵形式为

其中，Radon算子 L(λ)表达式如下：

采用最小二乘方法定义(9)式对应的λ-f域抛物Radon正变换为

式中，(λ)为Radon域的数据，I 为单位矩阵，α通常取主对角线上数值的1%，其作用是在矩阵求逆过程中作为稳定因子，增强求解的稳定性.对于正演的理论模型可以不加稳定因子直接计算，而对随机噪音较大的实际数据则需加入稳定因子.由于 L H(λ)L(λ)+αI是一个Toeplitz-Hermit型矩阵，可采用Levinson递推法计算，以节省存储空间，提高计算效率.

方程(9)和(11)即为λ-f域抛物Radon正反变换的公式，相比于传统抛物Radon变换，该方法的Radon算子 L ni(λ)是Nx×Nλ维的复数矩阵，与频率无关，所以只需计算一次，计算效率有很大程度的提高.同时，对于时空域具有一定抛物时差(q值为常数)的同相轴，该方法与传统Radon变换方法映射出的效果不同.传统抛物Radon变换由于曲率q值固定，截距时间τ由时空域抛物线顶点的时间位置确定，因此抛物同相轴在传统抛物Radon域中会映射为一个点.而由于q=λ/f，在λ-f域抛物Radon变换中，λ和f表现为正比例线性关系，因此，抛物同相轴在λ-f域表现为过原点以q为斜率的直线.

由定义λ=qf定义可知，当λ取值一定时，频率f的取值越小，曲率q的值越大，因此，相对于高频同相轴，λ-f域抛物Radon变换对低频同相轴有更大的曲率扫描范围，该方法可以实现低频数据能量在λ-f域的有效映射.而低频数据对改善地震数据的信噪比有重要作用，传统抛物Radon变换可能因低频同相轴无法满足采样定理，使得低频数据能量在Radon域不能聚焦.

3 λ-f域抛物Radon变换的离散采样和权系数

λ-f域抛物Radon变换的采样间隔Δλ和λ参数的范围是控制假频、稳定性、分辨率以及计算效率方面的重要参数.本文在计算时，采用

其中，xmax和xmin分别为数据的最大和最小偏移距.

当λ>0时，在方程(8)中Radon算子相位误差必须满足

因此，

同理，当λ<0时

所以，λ参数的范围为

另外，为改善Radon域重建剖面能量聚焦的收敛性，本文对待重建的地震道引入了变化的权系数，以改变缺失道在正变换过程中的参与程度.若用 W 表示权系数矩阵，则式(11)可写为

式(17)和(9)就构成了λ-f域加权抛物Radon变换的正反变换对.式(17)中引入权系数矩阵 W 后，并未破坏矩阵的Toeplitz结构，仍可使用原来的Levinson递推法求解，其计算量相对于λ-f域抛物Radon变换也未发生改变.权系数的选取准则为对于原始数据的规则道，其权系数的值始终为1.0.对于待重建的空道，其初始的权系数取为0，而后随着迭代次数的累积逐渐变化，依次取为0.4，0.7，1.0，自第4次迭代开始均取为1.0，一般迭代4次即可获得满意的效果.

图 1a为合成的含有三个同相轴的单炮地震记录，共100道，道间距20 m，时间采样间隔4 ms，图 1b为图 1a近偏移距缺失10道数据的剖面.对图 1a和图 1b所示的地震数据采用λ-f域抛物Radon正变换，得到图 2a和2b所示的λ-f域剖面.对图 1b所示的缺失道地震数据采用文中提出的λ-f域加权抛物Radon变换方法计算，得到图 2c所示的λ-f域剖面.对比图 2b和图 2c可知，不加权的方法计算缺失道地震数据会在λ-f域产生假象.而对缺失地震道数据应用加权方法后，其在λ-f 域显示了很好的聚焦性，与完整地震数据对应的λ-f 域变换的振幅基本一致.分析可知，本文方法对缺失多道的近偏移距数据的重建效果仍然较好.

图 1 合成的单炮记录 (a)合成地震剖面；(b)近偏移距缺失10道的剖面. Fig. 1 Synthetic single shot record (a)Synthetic seismic section;(b)Seismic section with 10 traces missing at the near offset.

图 2 λ-f域抛物Radon变换与λ-f域加权抛物Radon变换效果对比 (a)图 1a所示数据的λ-f域抛物Radon正变换剖面；(b)图 1b所示数据的λ-f域抛物Radon正变换剖面； (c)对图 1b所示数据应用文中方法重建的变换域剖面(迭代一次). Fig. 2 Result comparison between λ-f l domain parabolic Radon transform and λ-f l domain weighted parabolic Radon transform (a)λ-f l domain parabolic Radon forward transform of data showed in Fig. 1a;(b)λ-f l domain parabolic Radon forward transform of data showed in Fig. 1b;(c)Transformed domain section of data reconstructed by the method proposed in this paper(one iteration).
4 算例分析

4.1 理论模型1

基于波动方程的表面多次波预测方法SRME(Surface-related multiple elimination)要求全波场数据，为此，数据重建的效果直接影响着地震资料多次波预测和压制的效果.为测试文中提出的λ-f域加权抛物Radon变换重建地震数据的效果，合成了含表面多次波的模型数据，如图 3a所示.完整的理论模型数据共100道，其道间距为20 m，采样间隔为4 ms.人为的从模型数据第1道、第20道、第40道、第60道、第80道开始各剔除5道数据，如图 3b所示.对其采用文中所述的λ-f域加权抛物Radon变换方法进行重建，共迭代4次，重建后的炮记录如图 3c所示，重建效果较为理想，同相轴自然连续.由此可见，该方法对含多次波的近、中、远偏移距数据重建均有较好的适应性.

图 3 含多次波的单炮地震数据重建 (a)含有多次波的原始单炮记录；(b)待重建的单炮记录；(c)数据重建结果. Fig. 3 Seismic data reconstruction of shot record with multiples (a)Original shot record with multiples;(b)Shot record to be reconstructed;(c)Reconstruction result.

为了精确地对比数据重建的精度，抽取单道数据进行分析，图 3a中原始完整炮集记录的第3道数据和图 3c中应用文中方法重建后炮记录的第3道数据分别如图 4a和4b所示.可以观察到重建单道数据与原始单道数据在旅行时、振幅和相位方面均保持了较好的一致性.

图 4 单道重建效果对比 (a)图 3a所示原始炮记录的第3道数据; (b)图 3c所示重建后炮记录的第3道数据. Fig. 4 Result comparison of single trace reconstruction (a)The third trace of original shot record showed in Fig. 3a; (b)The third trace of shot record after reconstruction showed in fig3c.
4.2 理论模型2

本文抽取了Marmousi模型的CMP道集进行数据重建测试，图 5a为原始第224个CMP道集所示的地震数据，图 5b是该CMP道集前五道置零后的地震剖面.利用本文方法计算可得到图 5c所示地震数据剖面，可见重建的地震数据和原始地震数据相似度较高，重建效果较为理想.测试也表明文中方法不仅对理论模型一中的水平层状介质地震反射数据有效，对复杂介质的地震反射数据重建也较为有效.

图 5 Marmousi模型CMP224道集 (a)原始CMP道集；(b)近偏移距缺失5道的地震剖面；(c)数据重建的结果. Fig. 5 224th CMP gather of Marmousi model (a)Original CMP gather;(b)Seismic section with 5 traces missing at the near offset;(c)The reconstruction result.
4.3 实际地震数据

为进一步验证文中提出方法的有效性，抽取某海洋实际地震数据的第280个CMP道集，该CMP道集共65道，道间距为50 m，采样间隔为2 ms，近偏移距缺失5道，如图 6a所示，对其使用λ-f域加权抛物Radon变换的方法进行重建，结果如图 6b所示，重建效果较好，同相轴自然连续.

图 6 实际地震数据CMP道集重建 (a)近偏移距缺失的实际CMP道集；(b)数据重建结果. Fig. 6 CMP gather reconstruction of real seismic data (a)The real CMP gather with missing near offset;(b)Reconstruction result.
5 结 论

本文提出一种基于λ-f l域的加权抛物Radon变换地震数据重建方法，通过引入变量l λ，消除了常规Radon变换算子对频率的依赖关系.在迭代过程中，对待重建的地震道引入变化的权系数，在很大程度上降低了假频的影响，提高了变换域数据的重建精度.文中在λ-f l域抛物Radon正变换的求解过程中，利用了与地震数据无关的Toeplitz结构算子改善计算效率.理论模型及实际地震数据的试算表明，基于λ-f l域加权抛物Radon变换方法重建精度较高，重建结果与原始数据的振幅、相位保持一致，重建效果较为理想.

致 谢 感谢审稿专家和编辑部老师的指导和帮助.

参考文献
[1] Abbad B, Ursin B, Porsani M J. 2011. A fast, modified parabolic Radon transform[J].   Geophysics, 76(1): V11-V24.
[2] Chen B Y, Ma Z T. 1994. New technic for 3D prestack migration[J]. Chinese J. Geophys. (in Chinese), 37(3): 400-407.
[3] Duijndam A, Schonewille M A. 1999. Nonuniform fast Fourier transform[J].   Geophysics, 64(2): 539-551.
[4] Gao J J, Chen X H, Wang F F, et al. 2011. Study on regularized reconstruction of uneven seismic traces[J]. Progress in Geophys. (in Chinese), 26(3): 983-991.
[5] Hampson D. 1986. Inverse velocity stacking for multiple elimination[J]. Journal of the Canadian Society of Exploration Geophysicists, 22(1): 44-55.
[6] Hindriks K, Duijndam A. 2000. Reconstruction of 3-D seismic signals irregularly sampled along two spatial coordinates[J].   Geophysics, 65(1): 253-263.
[7] Jin S W, Xu S Y, Wu R S. 2002. Wave equation based prestack depth migration using generalized screen propagator[J]. Chinese J. Geophys. (in Chinese), 45(5): 684-690.
[8] Kabir M M N, Verschuur D J. 1995. Restoration of missing offsets by parabolic Radon transform[J]. Geophys. Prosp.  , 43(3): 347-368.
[9] Li X, Hu T Y. 2009. Surface-related multiple removal with inverse scattering series method[J]. Chinese J. Geophys. (in Chinese), 52(6): 1633-1640.
[10] Li X F, Li X F. 2008. Seismic data reconstruction with fractal interpolation[J]. Chinese J. Geophys. (in Chinese), 51(4): 1196-1201.
[11] Liu B, Sacchi M D. 2004. Minimum weighted norm interpolation of seismic records[J].   Geophysics, 69(6): 1560-1568.
[12] Liu D J, Yin X Y. 2007. A method of Fourier finite-difference preserved-amplitude prestack depth migration[J]. Chinese J. Geophys. (in Chinese), 50(1): 268-276.
[13] Liu H Y, Chen Y K, Li X F. 2014. Fractal interpolation method for reconstruction and resampling of seismic data[J]. Progress in Geophys. (in Chinese), 29(2): 518-522.
[14] Liu X W, Liu H, Liu B. 2004. A study on algorithm for reconstruction of de-alias uneven seismic[J]. Chinese J. Geophys. (in Chinese), 47(2): 299-305.
[15] Liu Y K, Sun H C, Chang X. 2004. Multiple removal by wavepath migration[J]. Chinese J. Geophys.   (in Chinese), 47(4): 697-701.
[16] Lu Y H, Lu W K, Zhai Z J. 2012. An edge-preserving seismic data interpolation method[J]. Chinese J. Geophys. (in Chinese), 55(3): 991-997.
[17] Ma J T, Sen M K., Chen X H. 2009. Multiple attenuation using inverse data processing in the plane-wave domain[J]. Chinese J. Geophys. (in Chinese), 52(3): 808-816.
[18] 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[J]. Chinese J. Geophys. (in Chinese), 51(1): 235-241.
[19] Sacchi M D, Ulrych T J. 1995. High-resolution velocity gathers and offset space reconstruction[J].   Geophysics, 60(4): 1169-1177.
[20] Schonewille M A, Duijndam A J W. 2001. Parabolic Radon transform, sampling and efficiency[J].   Geophysics, 66(2): 667-678.
[21] Trad D O, Ulrych T J, Sacchi M D. 2002. Accurate interpolation with high-resolution time-variant Radon transform[J].   Geophysics, 67(2): 644-656.
[22] Wang W H, Pei J Y, Zhang J F. 2007. Prestack seismic data reconstruction using weighted parabolic Radon transform[J]. Chinese J. Geophys. (in Chinese), 50(3): 851-859.
[23] Wang X W, Liu W Q, Wang Y C, et al. 2010. Research and application of common reflection angle domain pre-stack migration[J]. Chinese J. Geophys. (in Chinese), 53(11): 2732-2738.
[24] Xiong D, Zhao W, Zhang J F. 2009. Hybrid-domain high-resolution parabolic Radon transform and its application to demultiple[J]. Chinese J. Geophys. (in Chinese), 52(4): 1068-1077.
[25] Zhang W S, Guan Q, Song H B. 2001. Prestack depth migration by hybrid method with high precision and its parallel implementation[J]. Chinese J. Geophys. (in Chinese), 44(4): 542-551.
[26] 陈必远, 马在田. 1994. 三维叠前偏移新技术[J].   地球物理学报, 37(3): 400-407.
[27] 高建军, 陈小宏, 王芳芳,等. 2011. 不规则地震道数据规则化重建方法研究[J].   地球物理学进展, 26(3): 983-991.
[28] 金胜汶, 许士勇, 吴如山. 2002. 基于波动方程的广义屏叠前深度偏移[J].   地球物理学报, 45(5): 684-690.
[29] 李翔, 胡天跃. 2009. 逆散射级数法去除自由表面多次波[J].   地球物理学报, 52(6): 1633-1640.
[30] 李信富, 李小凡. 2008. 分形插值地震数据重建方法研究[J].   地球物理学报, 51(4): 1196-1201.
[31] 刘定进, 印兴耀. 2007. 傅里叶有限差分法保幅叠前深度偏移方法[J].   地球物理学报, 50(1): 268-276.
[32] 刘红艳, 陈宇坤, 李信富. 2014. 地震道重建和重采样的分形插值方法研究[J].   地球物理学进展, 29(2): 518-522.
[33] 刘喜武, 刘洪, 刘彬. 2004. 反假频非均匀地震数据重建方法研究[J].   地球物理学报, 47(2): 299-305.
[34] 刘伊克, Sun H C, 常旭. 2004. 基于波射线路径偏移压制多次波[J].   地球物理学报, 47(4): 697-701.
[35] 陆艳洪, 陆文凯, 翟正军. 2012. 一种边缘保持的地震数据插值方法[J].   地球物理学报, 55(3): 991-997.
[36] 马继涛, Mrinal K. Sen, 陈小宏. 2009. 平面波域反数据处理压制多次波方法研究[J].   地球物理学报, 52(3): 808-816.
[37] 孟小红, 郭良辉, 张致付,等. 2008. 基于非均匀快速傅里叶变换的最小二乘反演地震数据重建[J].   地球物理学报, 51(1): 235-241.
[38] 王维红, 裴江云, 张剑锋. 2007. 加权抛物Radon变换叠前地震数据重建[J].   地球物理学报, 50(3): 851-859.
[39] 王西文, 刘文卿, 王宇超,等. 2010. 共反射角叠前偏移成像研究及应用[J].   地球物理学报, 53(11): 2732-2738.
[40] 熊登, 赵伟, 张剑锋. 2009. 混合域高分辨率抛物Radon变换及在衰减多次波中的应用[J].   地球物理学报, 52(4): 1068-1077.
[41] 张文生, 关泉, 宋海斌. 2001. 高精度混合法叠前深度偏移及其并行实现[J].   地球物理学报, 44(4): 542-551.