地球物理学报  2014, Vol. 57 Issue (7): 2278-2290   PDF    
半径-斜率域加权反假频地震数据重建
黄小刚1,2, 王一博1, 刘伊克1, 常旭1    
1. 中国科学院地质与地球物理研究所工程地质力学重点实验室, 北京 100029;
2. 中国科学院大学, 北京 100049
摘要:为减小地震数据缺失给地震后续处理工作带来的影响,需要对地震数据进行插值重建.针对反假频插值重建这个难点问题,进行了相关研究,并由此提出了一种改进的R-P(半径-斜率)域加权反假频地震数据插值重建方法.该方法将F-K(频率-波数)谱变换到R-P域,在R-P域设计一个权函数并将其作用于每次的迭代插值过程.通过模型数据和实际数据的测试,证明了该方法具有较好的反假频插值重建能力.
关键词半径-斜率域     加权     反假频     重建    
Weighted anti-alias seismic data reconstruction in R-P domain
HUANG Xiao-Gang1,2, WANG Yi-Bo1, LIU Yi-Ke1, CHANG Xu1    
1. Key Laboratory of Engineering Geomechanics, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China;
2. University of Chinese Academy of Sciences, Beijing 100049, China
Abstract: In order to reduce the influence of the missing seismic traces on the follow-up seismic processes, seismic data reconstruction is required. To solve the problem of anti-alias seismic data reconstruction, some research about this was done and a new method called weighted anti-alias seismic data reconstruction in R-P (radius-slope) domain is proposed. The F-K (frequency-wave number) spectra are transformed to the R-P domain and a weighting function is designed. During every iteration process of the reconstruction, the R-P spectra are modified by the weighting function. By testing with synthetic and field data, the capacity of anti-alias reconstruction of this method is proved.
Key words: R-P     Weighted     Anti-alias     Reconstruction    
1 引言

在地震勘探野外工作中,由于各种原因,无法完 成对地震波场的充分采样(Abma and Kabir, 2006). 有时由于物理条件的限制,某些区域难以完成数据采集,例如高山及河流的存在区;另外,采集仪器的故障也可能导致同样的问题,例如一些检波器的故障会导致地震数据中存在若干坏道或空道.空间采样不足以及坏道的存在往往会影响到地震数据处理的效果,尤其是在波动方程偏移、表面相关多次波的消除、谱估计、时移地震数据处理等方面表现较为明显(Gao et al., 2010).利用有效数据恢复坏道或是加密地震数据显得尤为重要且必要.在这种背景下,地震道插值重建技术应运而生.

目前的地震数据重建方法有许多种,总体来说可分为两大类:基于波动方程的方法和基于信号处理的方法(曹静杰,2011).其中基于信号处理的方法又可以分为基于相干倾角、基于Fourier变换、基于线性预测、基于Radon变换、基于Curvelet变换和基于Seislet变换的方法等(曹静杰,2011唐刚,2010刘财等,2013).基于波动方程重建方法需要较大限度地利用地下信息,当地下信息难以获得或是不准确时,数据重建效果将不理想.另外,这种方法的计算量巨大从而难以投入实际生产(石颖,2010).Larner等(1981)首次提出了最大相干倾角插值的二维道内插方法.这种方法在时空域无法较好地处理多倾角的相交同相轴.Pieprzak等(1988)提出了一种处理多倾角同相轴的方法,不过这种方法应用于3D数据太过昂贵.基于Fourier变换重建方法的优点在于它是纯数据驱动的.基于线性预测的重建方法依据线性同相轴在f-x域的可预测性,将低频中抽取的自回归算子用于高频地震数据的重建.该类方法能够较好地处理假频问题,但仅限于规则采样的地震数据(Spitz,1991).为减小线性同相轴假 设的限制,发展出了基于Radon变换和基于Curvelet 变换的重建方法,在地震数据重建领域取得了较好 效果(Lu,1985Kabir, et al., 1995刘国昌等,2011).

地震数据重建中的难点之一在于反假频重建.对于一定的空间采样率,当地震信号的频率较低时,不容易产生假频;当频率较高时,假频就难以避免.不规则的地震道缺失会产生幅度较小的随机假频(Curry,2009),可视为频率-波数(F-K)域的随机噪声,通过一些常规的地震数据重建方法就能够成功地插值重建缺失道;当缺失道较为规则时,则可能形成幅度与有效频谱相当的规则假频,并且常与有效频谱混叠在一起.有时由于工作的需要,需给已有的数据隔道插入若干地震道,即内插.此时,要预先在待插道的位置设置零道.零道的插入会产生与有效频谱混叠在一起并且幅值与有效频谱相当的假频.我们可以将此内插问题等同于前述的含规则缺失道的地震数据插值重建问题.

为了进行反假频数据插值重建,一些学者进行了相关研究.目前看来,反假频重建方法大概有以下几类.第一类方法利用低频数据不易产生假频的性质,先设法重建低频数据,之后利用低频数据与高频数据的关系构建预测滤波算子,来重建高频数据,最终实现整个频段的数据重建.典型代表有Liu等(2004)刘喜武等(2004)提出的将最小加权范数内插(MWNI)与线性预测相结合的数据重建方法.Naghizadeh和Sacchi(2007)利用傅里叶方法重构无假频的低频谱,并由重构的低频分量提取预测滤波器预测高频分量.石颖等(20102013)给出了另一种由低频数据重建高频数据的方法,即最小加权范数插值与调制升频相结合的地震数据重建方法,利用最小加权范数插值方法重建无假频的低频数据,由重建的低频数据的AGC构建调制函数,并将调制函数作用于重建的低频数据,获得数据的高频成分.这一类方法能够一定程度地解决地震数据的反假频重建问题,但有时候由低频数据难以准确地构建高频数据.第二类方法就是稀疏促进的方法,将地震数据经过某种变换变换到稀疏域,在稀疏域施加一定的约束来压制假频,从而实现数据的反假频重建.Curry(2009)先将带假频的数据变换到频率-波数(F-K)域,再将F-K谱以极坐标的形式表示,形成半径-斜率(R-P)谱.在R-P谱中,他对谱施以半径相关的权重,获得了较好的反假频重建效果.Naghizadeh(2012)F-K谱中沿不同倾角进行能量扫描,找出若干个能量极值所对应的倾角,生成加权蒙板,作用于F-K谱,得到了较好的反假频重建效果.与第一类方法先重建低频成分,再由重建的低频成分重建高频成分不同的是,第二类方法插值重建始终是利用全频段的数据,这保证了高频成分的相对真实性.

凸集映射(POCS)方法是一种迅猛发展的图像 恢复方法(肖杰雄,2009Youla and Webb, 1982Oskoui and Stark, 1988杨李涛等,2010).早在1965年,Brègman在图像处理中首次提出了POCS的图像恢复方法(Brègman,1965).Abma和Kabir成功地将POCS方法应用于地震道插值重建(Abma and Kabir, 2006).该方法有几个优点:(1)可以恢复较大的数据缺失带,且能同时恢复地震信号和剖面中的随机噪声;(2)算法简单,易于实现;(3)可以利用FFT,计算效率较高(王淑琴,2010).目前,POCS地震道插值方法不仅能够对二维地震数据缺失道进行重建,也能用于三维地震数据,并且基于图形处理单 元(GPU),它的效率得到了显著提升(王淑琴,2010Wang et al., 2010).

Curry(2009)的方法利用了极坐标下有效频谱与假频容易区分的优势,取得了一定的反假频重建效果,但该方法需要通过共轭梯度法求解权函数,需要较大的计算量;Naghizadeh和Sacchi(2012)F-K谱沿倾角进行加权处理,能够较好地解决假频与有效频谱交叠的问题,但该方法求解反问题时没有利用POCS方法简洁和高效的优点.综合两者的优点,并考虑到POCS方法的优势,本文提出一种对极坐标表示下的谱亦即R-P(半径-斜率)谱沿斜率方向进行加权的POCS反假频地震数据重建方法,并通过合成数据和实际数据的测试,证明了方法的实用性. 2 方法原理

地震数据插值重建问题是一个欠定问题,理论上欠定问题有无穷多个解.对于欠定问题,可以根据先验信息如模型的稀疏性和光滑性等对问题施加一定的约束,从而获得满意的解.地震数据采集可以表示成:

其中,d 是完整的地震数据,亦即待重建得到的地震数据,S 是采样矩阵,d obs是得到的含有缺失道的地震数据.地震数据插值重建问题就是要由 d obs恢复出 d .Claerbout指出,缺失数据的重建就是使重建后的数据在某种滤波后具有最小的能量.由此,可以将重建问题表述为L-2范数下的问题:

其中,‖ · ‖2表示L-2范数.稀疏促进法对模型施以先验信息:模型在变换域是稀疏的,即模型在稀疏域的L-0范数最小,从而上述问题变为

其中,‖ · ‖0表示L-0范数,C 是某一稀疏变换算子,μ是一个调节前后两项权重的参数.(3)式是一个典型的NP-hard非凸优化问题,不易求解,所以通常用L-1范数来代替L-0范数.由此(3)式变成:

其中,‖ · ‖1表示L-1范数,μ依然是一个调节前后两项权重的参数.在解决此L-1范数约束下的稀疏反演问题上,POCS理论提供了较为方便的算法.

POCS算法是由Gerchberg和Saxton(1972)在图像处理中提出的阈值迭代算法,相应的凸约束集合的投影算子将解空间中的点投影到距离凸集表面最近的点上,经过有限次迭代,最终可以找到一个收敛于凸约束集合交集的解(李超峰等,2010).经过POCS方法重建的地震数据表达式为(Galloway and Sacchi, 2007):

式中,k是迭代次数,d k是第k次迭代后得到的插值重建结果,I 是单位阵,S 是采样矩阵,用于标记空道和有效道,T k是第k次迭代时使用的阈值算子,F 及 F -1分别是正反傅氏变换算子.

使用这种传统POCS方法能够较好地插值重建含有非规则缺失道的地震数据,但对于含有规则缺失道的地震数据,这种方法会失效.如图 1所示,(a、e)是含有非规则缺失道的地震数据及其F-K谱,(b、f)是采用传统POCS方法对A插值重建的结果及对应的F-K谱.传统POCS方法较好地重建了非规则的缺失道.(c、g)是含有规则缺失道的地震数据及其F-K谱,(d、h)是采用传统POCS方法对(c)插值重建的结果及对应的F-K谱.传统POCS方法对于含规则缺失道的数据基本无能为力.

图 1 传统POCS方法对含非规则缺失道及规则缺失道数据的插值重建效果
(a)含非规则缺失道的数据;(b)传统POCS对(a)的插值重建效果;(c)隔道缺失的地震数据; (d)传统POCS对(c)的插值重建效果;(e,f,g)和(h)分别为(a,b,c)和(d)对应的F-K谱.
Fig. 1 The reconstructed result of seismic data with r and omly missed traces and with evenly missed traces by traditional POCS method
(a)Data with r and omly missed traces;(b)The reconstructed result of(a)by traditional POCS method; (c)Data with evenly missed traces;(d)The reconstructed result of(c)by traditional POCS method. (e,f,g) and (h)are the F-K spectra of(a,b,c) and (d)respectively.

传统POCS方法对于含有规则缺失道的数据插值重建之所以失效是因为其含有幅值较强的假频,并且与有效频谱混叠在一起.为克服这一困难,恢复POCS方法的作用,需要对其进行改进.传统POCS方法能够有效插值重建含非规则缺失地震道的数据,利用了假频的幅值明显小于有效频谱幅值这一特点,但当缺失道规则时,这一特点不复存在.假频与有效频谱还有一重要区别在于有效频谱在F-K域都过原点或是延长线过原点.利用这个区别在极坐标表示的频率-波数谱,即R-P谱中可以非常方便地将假频和有效频谱分开.

对于含有缺失道的地震数据 d obs(x,t),通过二维FFT得到其F-K谱 D obs(f,k).对 D obs(f,k)施以直角坐标向极坐标的映射:

得到极坐标表示的频谱 D obs(r,p),即R-P(半径-斜率)谱.有效频谱在直角坐标下是过原点的直线,因此在极坐标下,它便变成了P为常数的直线;规则缺失道引起的假频在直角坐标下是不过原点的直线,在极坐标下,它便成了P不断变化的曲线.

在极坐标下,对不同的P进行能量扫描.能量扫描的具体表达式为

其中,abs表示复数取模.从扫描得到的所有 E(p)中选出数值最大的前m个值,记录其对应的P值.m的选择由数据窗内的同相轴倾角情况决定,对于一般的数据窗,m取值2~5较合适.利用最大的m个能量对应的P值进行权函数设计:

式中,b是沿着每条P线取窗的宽度,i是能量最大的m个P值对应编号,W t(r,p)是权重.当Pi和Pi+1相距较远时,权函数的两部分不重叠,两部分分别起作用.当Pi和Pi+1相距较小时,权函数的两部分可能发生重叠,将形成一个宽度大于b数值为1的重叠区,此时,权函数依然能够正常起作用.

在使用POCS方法进行插值重建的每次迭代过程中,将上述权函数作用于R-P谱,便能获得较好的反假频插值重建结果.

经过上述加权处理后,得到R-P域加权POCS反假频重建迭代公式为

式中,W t是权函数,R p及 R -1p是正反直角坐标-极坐标变换算子. 3 数据测试 3.1 线性同相轴的反假频插值重建

为去除其他因素的干扰(如噪声、振幅差异等),首先选取由褶积生成的含线性同相轴的模型数据来测试本方法的效果.如图 2a所示的数据是由主频40 Hz的雷克子波与反射系数序列褶积形成的.该记录有40道,每道样点数为200,采样时间间隔为4 ms,数据中包含有3条线性同相轴.图 2d是2a对应的频率-波数(F-K)谱.将2a隔道抽空,置为零值,得到2b所示的地震数据.由2b对应的F-K图 2e可见,在隔道抽空置零的过程中,带来了严重的假频,而且假频与有效频谱产生了交叠.若使用传统POCS方法,将无法插值重建出空道.图 2c为使用本文方法插值重建所得的结果,2f是其对应的F-K谱.图 3是插值重建的第7道(蓝线)与原始第7道(红线)的单道对比图,蓝线与红线基本重合.除了微弱的噪声外,插值重建的结果无论在振幅还是相位上都基本与原始数据吻合.为量化插值重建误差,引入能量误差参数Error.

图 2 R-P域加权反假频方法对线性同相轴的插值重建效果
(a)原始数据;(b)将奇数道置零后的数据;(c)本文方法插值重建所得的数据;(d,e,f)分别是(a,b)和(c)的F-K谱.
Fig. 2 The reconstructed result of the data with straight line events processed by the weighted anti-alias reconstruction in R-P domain method
(a)The original data;(b)data by zeroing the odd traces;(c)The reconstructed result by the proposed method; (d,e,f)are the F-K spectra of(a,b) and (c)respectively.

图 3 第7道插值前后单道波形对比(红线表示原始波形,蓝色表示插值所得波形) Fig. 3 The comparison of the original(red)7th trace and the 7th trace acquired by interpolation(blue)

式中 d interp表示重建所得的数据,d origin表示抽空前的原始数据,E(·)表示求能量.通 过计算得到图 2c图 2a的Error为0.4283%.对比F-K 谱可知,本方法有效地压制了假频.可见,本方法对于含有规则缺失道,带有严重假频的线性同相轴具有良好的插值重建效果.

图 4a所示为图 2b所对应的极坐标表示的频谱,即半径-斜率(R-P)谱.图 4a中有三条P值不变的垂线(P值分别在-0.2,0,0.2附近),分别对应图 2e中三条由坐标原点(f=0,k=0)向外辐射的射线,它们是有效频谱中的主倾角方向.图 4a中其余曲线的P值都在不断变化,基本为假频.由P值能量扫描得到了图 4b所示的结果,能量主峰分别在P=-0.2,0,0.2附近,进而生成了4c所示的权函数,黑线条处值为1,其余为0.在POCS每一次迭代过程中使用图 4c所示的权函数对R-P谱进行加权处理,便能够得到图 2c所示的插值重建结果.图 4d是使用权函数对R-P谱(图 4a)进行了一次处理后的R-P谱.

图 4(a)图 2b所对应的R-P谱;(b)由能量扫描所得的能量随斜率分布图;(c)由能量扫描生成的权函数图;(d)权函数c对a作用一次后的R-P Fig. 4(a)The R-P spectra of Fig. 1 b;(b)Energy-slope(E-P)plot produced by energy scanning; (c)The weighting function produced by(b);(d)The R-P spectra of a modified once by the weighting function

图 5,使用主频12 Hz和18 Hz的雷克子波褶积叠合形成的地震记录,采样间隔1 ms.对其每隔一道剔除三道后再随机剔除50%的地震道,所得结果如图 5a所示.使用Curry的方法对图 5a的数据进行插值重建,迭代10次后所得结果如图 5b所示.由图可见,虽然缺失道都得以较好恢复,但插值重建带了明显的虚假信息,表现为产生了大量平行于原有同相轴的微弱同相轴.而且,这种方法采用共轭梯度求取权函数,耗时较长.使用本文方法对5a进行插值重建,迭代10次,所得结果如图 5c所示.由图发现,除了缺失道得以较好地恢复外,图 5b中出现的那些平行于原有同相轴的微弱同相轴也基本被消除了.更重要的是,本方法避免了共轭梯度法求取权函数,节约了一些计算时间.

图 5(a)含有87.5%缺失道的数据(先每隔一道剔除三道,即先均匀地剔除75%,再随机地剔除50%地震道); (b)使用Curry的方法迭代10次所得结果;(c)使用本文方法迭代10次所得结果. Fig. 5(a)Seismic data with 87.5% missed traces(removing 75% traces in a regular fashion,then removing 50% traces r and omly);(b)Reconstruction result by Curry′s method;(c)Reconstruction result by the proposed method
3.2 双曲同相轴的反假频插值重建

双曲同相轴在地震道集中很常见,研究本方法对于双曲同相轴的效果很有必要.由主频15 Hz的雷克子波有限差分正演形成一个共炮点记录.该记录有151道,每道样点数为1200,采样时间间隔为1 ms,数据中包含有3条双曲同相轴.为研究本文方法对于随机缺失、规则缺失和随机缺失叠加规则缺失三种不同缺失情况的插值重建效果,进行了以下三种测试:

将上述合成数据随机抽空15%的地震道,置为零值,得到图 6a所示的地震数据.由图 6a对应的F-K图 7a可见,在随机抽空置零的过程中,带来了微弱的频谱噪声.使用传统POCS方法,就能插值重建出空道数据.现使用本文方法对图 6a进行分窗插值重建,窗的大小为:空间窗点数28,时间窗点数800,得到结果如图 6b所示.插值重建道与原始道保持了振幅和相位的较好一致性.插值误差见表 1.由单道波形对比图 8a可知,插值重建结果无论在振幅还是相位上都基本与原始数据吻合.可见,与传统POCS方法一样,本方法对含随机缺失地震道的地震数据依然有效.

图 6 R-P域加权反假频方法对含双曲同相轴的三种不同缺失道地震数据的插值重建效果
(a)随机缺失15%地震道的数据;(c)奇数道缺失的数据;(e)随机缺失15%的地震道后再将奇数道抽空的数据; (b,d)和(f)分别是使用R-P加权方法对(a,c)和(e)插值重建的结果.
Fig. 6 The reconstructed result of the data with hyperbolic events for three different kinds of missed traces processed
by the weighted anti-alias reconstruction in R-P domain method. (a)Data with 15% r and omly missed traces;(c)Data by zeroing the odd traces;(e)Data by zeroing the odd traces after r and omly zeroing 15% the traces;(b,d) and (f)are the reconstructed result by the proposed method for(a,c) and (e)respectively.

图 7(a、b、c、d、e)和(f)分别是图 6中(a、b、c、d、e)和(f)所对应的F-K Fig. 7(a,b,c,d,e) and (f)are the F-K spectra of(a,b,c,d,e) and (f)of Fig. 6 respectively

表 1 图 6中数据插值误差列表 Table 1 Error of the interpolation for the data in Fig. 6

将原始数据隔道抽空,置为零值,得到图 6c所示的地震数据.由图 6c对应的F-K图 7c可见,在隔道抽空置零的过程中,带来了严重的假频.若使用传统POCS方法,将无法插值重建出空道数据.使用本文方法对图 6c进行分窗插值重建,窗的大小设置同图 6a,得到结果如图 6d所示.插值重建道与原始道保持了振幅和相位的较好一致性.插值误差见表 1.由单道波形对比图 8b可知,插值重建结果在振幅和相位上都基本与原始数据吻合.对比图 6c图 6dF-K图 7c图 7d可知,本方法较有效地压制了假频.可见,本方法对由规则缺失道带来严重假频的双曲同相轴数据具有良好的反假频插值重建效果.

图 8 插值前后的单道波形对比
(a)图 6中(a)与(b)的第57道;(b)图 6中(c)与(d)的第77道;(c)图 6中(e)与(f)的 第61道,红线表示原始波形,蓝色表示插值所得波形.
Fig. 8 The comparison of the original(red)trace and the trace acquired by interpolation(blue)
(a)The 57th trace of(a) and (b)in Fig. 6;(b)The 77th trace of(c) and (d)in Fig. 6;(c)The 61th trace of(e) and (f)in Fig. 6.

将原始数据随机抽空15%的地震道,之后将数据隔道抽空,得到图 6e所示的地震数据.由图 6e对应的F-K图 7e可见,在随机抽空和隔道抽空置零的过程中,既带来了频谱噪声也带来了严重的规则假频.若使用传统POCS方法,只能重建出随机缺失的地震道,而无法重建出规则缺失的道.使用本文方法对图 6e进行分窗插值重建,窗的大小设置同图 6a,得到结果如图 6f所示.除了一些连续缺道较多的地方外,插值重建道与原始道都保持了振幅和相位的较好一致性.插值误差见表 1.由单道波形对比图 8c可知,插值重建结果无论在振幅还是相位上都基本与原始数据吻合.对比F-K图 7e图 7f可知,本方法既较好地压制了频谱噪声,又较有效地压制了假频.可见,本方法对同时含有规则缺失道和随机缺失道的双曲同相轴数据仍具有良好的反假频插值重建效果. 3.3 实际数据的反假频插值重建

图 9所示,某工区包含531道,每道时间样点数为1501,采样间隔为2 ms的叠后地震资料隔道抽空形成一个地震剖面.为使同相轴为线性或近线性的,对图 9的数据进行开窗处理,数据窗的大 小为:空间窗40个样点,时间窗1501个样点.对每 个数据窗使用本文方法进行反假频插值重建,得到的结果如图 10所示.由图 10看出:经过本文方法插值重建后,各缺失道均得到了较好的恢复,无论是同相轴的振幅还是相位都能够与附近道较好地匹配,保持了同相轴的较好连续性.利用公式(10)计算得插值重建的误差Error为13.9967%.图 11是插值重建前后的单道波形对比图.可见,经过插值重建后,除个别小振幅处的波形外,插值重建所得地震道波形基本与原始地震道吻合.为进一步对比插值重建效果,从图 9和10中同一位置分别取出了一个数据窗,如图 12a和12b所示.如图 12c的虚线右侧部分展示的是12b与抽空前的原始数据同一位置数据窗的差异剖面.为了在同一增益下对比误差道的幅度与原始道的幅度,将原始数据窗放在了图中,如图 12c中虚线左侧部分.可见,虽然存在一定的插值重建误差,但本文方法还是较好地压制了假频,成功地插值重建出了规则缺失的地震道,保持了同相轴的较好连续性.实际资料测试表明:对于带严重假频的实际资料,本方法依然具有良好的反假频插值重建效果.

图 9 奇数道置零后的叠加剖面 Fig. 9 The stacked data by zeroing the odd traces

图 10图 9使用R-P加权插值重建方法重建所得结果 Fig. 10 The reconstructed data by the weighted anti-alias reconstruction in R-P domain method for Fig. 9

图 11 第137道插值前后单道波形对比(红线表示原始波形,蓝色表示插值所得波形) Fig. 11 The comparison of the original(red)137th trace and the 137th trace acquired by interpolation(blue)

图 12(a)和(b)分别是图 9和10中同一位置的数据窗;(c)中虚线左侧是抽空前的原始数据窗,虚线右侧是插值的数据窗与原始数据窗做差的结果 Fig. 12(a) and (b)are windows from the same position of Fig. 8 and 10 respectively,the part on the left of the dotted line in(c)is the original windowed data and the right part is the difference between the reconstructed data and the original data
4 讨论

本文提出的R-P域加权反假频地震数据插值重建方法与之前的一些反假频方法相比,具有一些优点.第一,与先重建低频成分再由低频成分重建高频成分的反假频重建方法相比,本方法利用了全频段信息,保证了重建所得高频成分的相对真实性;第二,与Curry提出的在极坐标下沿径向加权的方法相比,本方法采用能量扫描来计算权函数,省去了通过共轭梯度法迭代求解权函数的过程,节省了计算量且实现更简便;第三,与Naghizadeh提出的对F-K谱沿倾角进行加权处理的方法相比,本方法利用了R-P域假频与有效频谱容易区分的优点,另外本方法基于POCS方法来求解反问题,具有相对的简单性和高效性.

本方法的一些参数选取对结果有一定影响.随着迭代次数的增加,处理效果会不断改善,但迭代次数的增加也会导致计算量的增大.基于傅氏变换的POCS方法要求同相轴是线性的,因而当待处理的数据同相轴为非线性时要做数据加窗处理.数据窗较大时,抗干扰能力较强,即受噪声的影响较小,但若窗太大,窗内同相轴难以保证是线性或近似线性的,起不到加窗的作用;窗较小,可以保证同相轴的良好线性性质,从而可以改善处理效果,但窗较小时抗干扰能力就较弱,噪声容易对结果造成影响.所以窗大小的选择要视情况而定.公式(8)中b值的选取 也对结果有一定影响.本文为了简单将介于-0.5b~0.5b 区间的权系数都置为1.为了减小噪声,改善结果,也可选用其他更优的权值设置方法.

缺失道的存在会引起频谱的变化,产生频谱噪声,因而可以将缺失道视为一种噪声,插值重建的过程也就成了去噪的过程.实际上,去噪与插值重建本质上是一样的,所以本文提出的反假频地震数据插值重建方法完全可以用于反假频去噪. 5 结论

通过研究,可以得出以下一些初步认识:

(1)极坐标表示的频率-波数谱即半径-斜率(R-P)谱具有方便区分有效频谱和假频的优点,有效频谱在R-P域中表现为P值不变的垂线,而假频则为P值不断改变的曲线.

(2)传统POCS方法能重建仅含非规则地震缺失道的数据,对带假频的数据无能为力;R-P域加权的反假频POCS地震数据插值重建方法不仅能够插 值重建仅含非规则缺失道的数据,还能够插值重建仅含规则缺失道带假频的数据,亦能够插值重建规则缺失叠加非规则缺失的地震数据,对于有效频谱与假频相交叠的情况,该方法也具有良好的插值重建效果.

(3)参数的选取对结果有一定影响,需要根据实际情况合理选择参数,改善插值重建效果.

致谢 本研究的顺利完成需要感谢国家自然科学基金(40930421,41074091)项目的支持.非常感谢审稿人的宝贵意见和建议.另外,在研究过程中,卢回忆、谢宋雷 等几位博士、硕士也给予了帮助和支持,在此一并感谢.
参考文献
[1] Abma R, Kabir N. 2006. 3D interpolation of irregular data with a POCS algorithm. Geophysics, 71(6): E91-E97, doi: 10.1190/1.2356088.
[2] Brègman L M. 1965. The method of successive projection for finding a common point of convex sets. Soviet Math, 6(3): 688-692.
[3] Cao J J. 2011. Seismic data compressed restoration and sparse inversion algorithm research (in Chinese). Beijing: Graduate University of Chinese Academy of Sciences.
[4] Curry W J. 2009. Interpolation with Fourier-Radial adaptive thresholding. SEG Expanded Abstracts, 3259-3263.
[5] Galloway E J, Sacchi M D. 2007. POCS method for seismic data reconstruction of irregularly sampled data.//Abstract of 2007-CSPG CSEG Conference. Calgary, Canada, 555.
[6] 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.
[7] Gerchberg R W, Saxton W O. 1972. A practical algorithm for the determination of phase from image and diffraction plane pictures. OPTIK, 35(2): 237-246.
[8] Kabir M M, Verschuur D J. 1995. Restoration of missing offsets by parabolic Radon transform. Geophysical Prospecting, 43(3): 347-368, doi: 10.1111/j.1365-2478.1995.tb00257.x.
[9] Larner K, Gibson B, Rothman D. 1981. Trace interpolation and the design of seismic surveys. Geophysics, 46(4): 407-415.
[10] Li C F, Liu H, Zhou F. 2010. Reconstruction of image super-resolution based on POCS algorithm. Microcomputer & Its Applications (in Chinese), (3): 20-22.
[11] Liu C, Li P, Liu Y, et al. 2013. Iterative data interpolation beyond aliasing using seislet transform. Chinese J. Geophys. (in Chinese), 56(5): 1619-1627, doi: 10.6038/cjg20130519.
[12] 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.
[13] Liu H, Liu B. 2004. Dealiasing band-limited uneven data using a linear prediction and adapted weighted least square. CPS/SEG Beijing 2004 international geophysical conference and exposition, 198-201.
[14] Liu X W, Liu H, Liu B. 2004. A study on algorithm for reconstruction of de-alias uneven seismic data. Chinese J. Geophys. (in Chinese), 47(2): 299-305.
[15] Lu L. 1985. Application of local slant-stack to trace interpolation. 57th Ann. Internat. Mtg., Soc. of Expl. Geophys., Expanded Abstracts, 560-562.
[16] Naghizadeh M, Sacchi M D. 2007. Multistep autoregressive reconstruction of seismic records. Geophysics, 72(6): v111-v118, doi: 10.1190/1.2771685.
[17] Naghizadeh M. 2012. Seismic data interpolation and denoising in the frequency-wavenumber domain. Geophysics, 77(2): v71-v80, doi: 10.1190/Geo2011-0172.1.
[18] Oskoui-Fard P, Stark H. 1988. Tomographic image reconstruction using the theory of convex projections. IEEE Transactions on Medical Imaging, 7(1): 45-58.
[19] Pieprzak A W, McClean J W. 1988. Trace interpolation of severely aliased events. 58th Ann. Internat. Mtg., Soc. of Expl. Geophys., Expanded Abstracts.
[20] Shi Y. 2010. Pragmatic multiple prediction and suppression investigation (in Chinese). Beijing: Graduate University of Chinese Academy of Sciences.
[21] Shi Y, Zhang Z, Wang J M, et al. 2013. Investigation on anti-aliasing regularization approach for seismic data. Progress in Geophys. (in Chinese), 28(1): 250-256, doi: 10.6038/pg20130126.
[22] Spitz S. 1991. Seismic trace interpolation in the F-X domain. Geophysics, 56(6): 785-794, doi: 10.1190/1.1443096.
[23] Tang G. 2010. Seismic data reconstruction and denoising based on compressive sensing and sparse representation (in Chinese). Beijing: Tsinghua University.
[24] Wang S Q, Gao X, Yao Z X. 2010. Accelerating POCS interpolation of 3D irregular seismic data with Graphics Processing Units. Computers & Geosciences, 36(10): 1292-1300, doi: 10.1016/j.cageo.2010.03.012.
[25] Wang S Q. 2010. Accelerating geophysical processing with Graphics Processing Units--application to 3D depth migration and POCS interpolation (in Chinese). Beijing: Graduate University of Chinese Academy of Sciences.
[26] Xiao J X. 2009. A POCS algorithm for super-resolution image reconstruction (in Chinese). Shanghai: Shanghai Jiaotong University.
[27] Yang L T, Lu L J, Fan Z Y. 2010, Video images super resolution reconstruction based on improved POCS algorithm. Microcomputer Applications (in Chinese), 26(7): 33-36.
[28] Youla D C, Webb H. 1982. Image restoration by the method of convex projections: part 1-theory. IEEE Transactions on Medical Imaging, 1(2): 81-94.
[29] 曹静杰. 2011. 地震数据压缩重建及稀疏反演算法研究[博士论文]. 北京: 中国科学院研究生院.
[30] 李超峰, 刘辉, 周峰. 2010. 基于POCS算法的图像超分辨率重建. 微型机与应用, (3): 20-22.
[31] 刘财, 李鹏, 刘洋等. 2013. 基于seislet变换的反假频迭代数据插值方法. 地球物理学报, 56(5): 1619-1627, doi: 10.6038/cjg20130519.
[32] 刘国昌, 陈小宏, 郭志峰等. 2011. 基于curvelet变换的缺失地震数据插值方法. 石油地球物理勘探, 46(2): 237-246.
[33] 刘喜武, 刘洪, 刘彬. 2004. 反假频非均匀地震数据重建方法研究. 地球物理学报, 47(2): 299-305.
[34] 石颖, 张振, 王建民等. 2013. 地震数据反假频规则化方法研究. 地球物理学进展, 28(1): 250-256, doi: 10.6038/pg20130126.
[35] 石颖. 2010. 多次波预测与压制的实用方法[博士论文]. 北京: 中国科学院研究生院.
[36] 唐刚. 2010. 基于压缩感知和稀疏表示的地震数据重建与去噪[博士论文]. 北京: 清华大学.
[37] 王淑琴. 2010. 利用图形处理单元加速地震数据处理研究——以三维深度偏移和POCS插值为例[博士论文]. 北京: 中国科学院研究生院.
[38] 肖杰雄. 2009. 基于POCS算法的超分辨率图像重建[硕士论文]. 上海: 上海交通大学.
[39] 杨李涛, 路林吉, 范征宇. 2010. 基于改进POCS算法的视频图像超分辨率重建. 微型电脑应用, 26(7): 33-36.