双曲Radon域综合预测反褶积和反馈循环方法压制多次波模型构建方法与流程

本发明涉及海洋地质工程领域，特别涉及双曲radon域综合预测反褶积和反馈循环方法压制浅地层多次波模型构建方法。

背景技术：

由于受到测量方式和海底浅地层沉积物结构的复杂性影响，浅地层剖面仪(浅剖)换能器接收阵不可避免地会接收到一些来自海面和浅地层层间的多次反射波，给浅地层有效一次波分离带来了困难，也严重地影响了浅地层成像的真实性和可靠性，进而误导了后期的浅地层底质解释。目前压制多次波的方法主要有两类：一类为基于有效波与多次波特征和性质之间差异的滤波方法，此类包括预测反褶积、fk变换、radon变换等(fosterandmosher,1992；yilmazandtaner,1994)；另一类为基于波动理论的预测相减法，此类包括波场外推法、反馈循环法、反散射级数法等(morleyandclaerbout,1983；wegleinetal.,1997)。预测反褶积原理是根据多次波具有周期性而一次波不具备周期性的差异分离出来的，适用于具有较大正速度梯度的介质，具有算法简单、计算成本低的特点。taner提出线性radon域预测反褶积方法，提高了压制中远偏移距多次波的能力，之后又将单道预测反褶积发展到二维，利用相邻地震道与不同滤波因子的褶积预测一次波，增加了对地层横向变化的适应能力(宋家文etal.,2014)。lokshtanov证实线性radon域预测反褶积也能够有效衰减多次波(lokshtanov,1995；1999)。赵昌垒等分析了线性radon预测反褶积的参数选取与适用性问题，并对南海深水地震资料多次波压制取得一定的效果(赵昌垒etal.,2013)。verschuur和berkhout提出压制多次波的反馈循环法，整个预测过程无需知道宏观速度场，从而增加了预测方法的适应性。随后两人又基于反馈模型理论提出2d级数展开压制表面多次波算法，后面发展提升了该理论，在无需对多个预测项求和的情况下，基于反馈模型理论提出迭代压制表面多次波算法，在某种程度上提高了计算精度(verschuur,1992；verschuuretal.,1995；berkhoutandverschuur,1997；verschuurandberkhout,1997；berkhout,1999)。基于该思想，kelamis、wang、vangroenestjin等先后对基于反馈循环法理论和应用进行了有益的完善和改进(kelamisandverschuur,2000；wang,2004；2007；vangroenestijnandverschuur,2009a；vangroenestijnandverschuur,2009b)。radon变换压制多次波方法来源于ryu提出的速度滤波概念，当多次波和一次波存在速度差时，多次波能量与一次波能量成像至不同的速度空间，在速度空间里，可将一次波和多次波分离。薛亚茹提出一种基于radon变换和正交多项式变换的多方向正交多项式变换压制多次波方法，该方法仅用一个曲率参数便可描述同相轴剩余时差参数，较大地提高了一次波和多次波的剩余时间分辨率(薛亚茹etal.,2012)。巩向博针对在速度各异及长偏距情况下，地震数据在radon域内能量仍不收敛的问题，提出了各向异性radon变换压制多次波方法，推导了由偏移距、慢度、非椭圆率三参数控制的积分曲线正反变换公式，避免了时间域稀疏脉冲迭代反演中的大矩阵运算，保持了较高的精度，并提高了计算效率(巩向博etal.,2014)。

然而，根据上面的描述，现有技术存在以下缺陷：

(1)当浅地层的声波速度逆转或者横向速度变换剧烈时，预测反褶积未能有效压制多次波，甚至会损伤有效波；

(2)反馈循环法在消除多次波时需要获知已知发射换能器的子波和地下介质结构，而地下介质结构往往较难准确获取，从而影响到是否精准预测并消除多次波；以及

(3)radon变换不适用于非水平均匀层状介质产生的多次波压制，以及对微屈多次波压制效果较差。

技术实现要素：

本发明要解决的技术问题是在全面分析预测反褶积法、反馈循环法和radon变化方法优缺点的基础上，结合获取的浅地层剖面数据特点，提供一种双曲radon域综合预测反褶积和反馈循环法压制多次波的模型构建方法，实现多次波与一次波的有效分离，为后期的海底地下介质解读提供更可靠的数据。

为实现上述发明目的，本发明要解决的技术问题包括如下几个方面：

(1)同时分别运用反馈循环法和预测反褶积法预测多次波，再将两者预测的多次波成分和预测误差进行双曲radon变换。

(2)针对预测反褶积预测误差中包含有效反射波和未彻底压制的长周期多次波，设计自适应滤波器，将有效反射波区域切除，保留多次波能量。

(3)当自适应滤波器输出值fm(τ,h)近似为0时，认为还存在有效反射波，则需返回再迭代处理，直至fm(τ,h)约等于1，结束迭代处理。之后，结合两种方法获得的多次波道数据模型，获得双曲radon域统一的多次波数据模型，再通过反双曲radon变换，最后从原始浅剖数据减去多次波数据即可完成多次波的综合压制过程。

为此，本发明提供的双曲radon域综合预测反褶积和反馈循环法压制多次波的模型构建方法，包括以下步骤：

(1)读入原始浅地层剖面数据，首先分别运用预测反褶积和反馈循环方法预测多次波成分；

(2)对以上两种方法预测的多次波进行双曲radon变换；

(3)针对预测反褶积预测误差中仍然存在长周期多次波，构建自适应滤波器，滤除一次有效波，进而获得运用两种方法之后的双曲radon域多次波能量模型；

(4)当自适应滤波器输出值fm(τ,h)近似为0时，认为还存在有效反射波，则需返回再迭代处理，直至fm(τ,h)约等于1，结束迭代处理。

(5)在双曲radon域中对多次波能量模型进行反radon变换；

(6)用原始浅地层剖面数据减去多次波模型数据，最终获得多次波压制后的高可靠性浅地层剖面数据。

在本发明的一个实施例中，预测反褶积的数学模型为：

在本发明的一个实施例中，反馈循环的数学模型为：

在本发明的一个实施例中，双曲radon变换的数学模型为：

在本发明的一个实施例中，自适应滤波器的数学模型为：

通过上述技术方案，本发明的有益效果是：

(1)针对浅海地区的浅地层剖面数据，由于受到海面自由波和海底浅地层层间多次波影响较严重，综合分析了预测反褶积法和反馈循环法的优缺点，利用两者的优势互补，提出一种双曲radon域综合预测反褶积和反馈循环法压制多次波的模型构建方法。

(2)设计出一种自适应滤波器，将预测反褶积预测误差中的有效波切除，进而得到成分更高的多次波能量模型，再利用原始数据减去多次波能量模型数据，最终达到有效压制多次波的目的。

附图说明

为了更清楚地说明本发明的技术特征，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的相关附图。

图1是预测反褶积压制多次波的流程图；

图2是预测反褶积和反馈循环综合压制多次波流程。

具体实施方式

为了使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，下面结合具体图示，进一步阐述本发明。

首先，本发明涉及如下技术术语：

预测反褶积法

预测反褶积在某种意义上可以说是一种广义的最小平方反褶积，它基于多次波具有周期性的特点，并通过设计预测误差滤波器消除原始数据中周期性出现的成分的过程(王兆旗,2007)。预测反褶积可应用于叠前记录数据，也可应用于叠后记录数据。

反馈循环法

反馈循环法是基于自由界面和层界面模型来预测多次波的，无需浅地层介质信息，即便多次波和有效反射波的时间差异不大，也能预测多次波(沈操,2003)。一般应用于叠前记录数据。

双曲radon域变换

双曲radon域变换应用地震共炮点或共中心点记录的双曲线特征，沿着双曲线路径首先对数据进行时间轴平方变换(在时间轴平方坐标下，时间轴平方与截距时间平方呈线性关系)，可使用频率域的算子进行radon正、反变换，变换后结果再进行时间轴开方变换即获得双曲radon变换结果(巩向博etal.,2016)。

本发明提供的双曲radon域综合预测反褶积和反馈循环方法压制多次波的模构建方法包括以下步骤：

(1)分别运用预测反褶积和反馈循环方法预测多次波成分；

(2)对以上两种方法预测的多次波进行双曲radon变换；

(3)针对预测反褶积预测误差中仍然存在长周期多次波，构建自适应滤波器，滤除一次有效波，进而获得运用两种方法之后的双曲radon域多次波能量模型；

(4)当自适应滤波器输出值fm(τ,h)近似为0时，认为还存在有效反射波，则需返回再迭代处理，直至fm(τ,h)约等于1，结束迭代处理。

(5)在双曲radon域中对多次波能量模型进行反radon变换；

(6)用原始浅地层剖面数据减去多次波模型数据，最终获得多次波压制后的高可靠性浅地层剖面数据。

参见图1至图2所示，本发明的一个具体实施方式如下：

(1)总体的技术路线

图2为双曲radon域综合预测反褶积和反馈循环方法压制多次波模型构建方法的原理图，首先将原始浅剖数据分别通过预测反褶积和反馈循环法进行预处理，利用反馈循环法在原始浅剖数据中粗略预测多次波的位置，作为多次波道数据模型，预测出来的多次波道数据模型无需与原始数据的多次波完全严格匹配。同时，结合预测反褶积的处理结果，由于预测误差中包含有效反射波和未彻底压制的长周期多次波，多次波能量主要分布在具有较大曲率参数的道上(曾忠玉,2013)。其次，将两者预测的多次波成分和预测误差进行双曲radon变换。然后通过本发明设计的自适应滤波器，将有效反射波区域切除，保留多次波能量。当fm(τ,h)近似为0时，如果还存在有效反射波，则需返回再迭代处理，直至fm(τ,h)约等于1，结束迭代处理。之后，结合两种方法获得的多次波道数据模型，获得双曲radon域统一的多次波数据模型，再通过反双曲radon变换，便可获得滤波后的多次波数据模型，即保留了多次波同相轴存在的浅剖数据。最后从原始浅剖数据减去多次波数据即可完成多次波的综合压制过程。

(2)预测反褶积理论

预测反褶积是基于多次波具有周期性的特点，并通过设计预测误差滤波器消除原始数据中周期性出现的成分的过程(王兆旗,2007)。该方法对在零偏移距处具有较好周期性的多次波压制效果显著，而对在中远偏移距处周期性的多次波压制效果不佳，所以此方法难以较好地压制整个偏移距范围内的多次波。

预测反褶积基本原理为，设预测步长为l的浅剖记录如式(9)所示，第二个

等号右边第一项表示预测误差(与一次反射波吻合)，第二项即表示预测的多次波(刘军,2008)。

又设d(t)为浅剖信号当前值，d(t-i)，i＝1,2,3,......为浅剖信号过去值，利用浅剖信号当前值和过去值综合预测其将来值即：

上式中t表示变换，针对线性时不变系统，预测值可用褶积模型来表示，即：

其中，c(t)＝{c(0),c(1),...}为线性预测因子。对于t+l时刻，预测误差为：

上式中d(t+l)是t+l时刻浅剖信号期望值，为实际输出值，即为多次波。选择合适的预测因子预测多次波，并且利用原始浅剖数据减去多次波成分即可得到有效反射一次波，如式(12)中的预测误差e(t+l)为有效一次波。预测反褶积方法压制多次波的流程如图1所示，首先输入浅地层剖面数据记录及参数，并取第k个时窗数据；其次求解出toeplitz矩阵；然后解线性方程组求滤波因子c(m)；最后求反褶积因子与时窗中每道记录的卷积，直到所有道数据处理完之后输出浅剖记录数据。

(3)反馈循环法理论

反馈循环法是基于自由界面和层界面模型来预测多次波的，无需浅地层介质信息，即便多次波和有效反射波的时间差异不大，也能预测多次波。在利用反馈循环法压制多次波时需满足的两个条件为：①需要修补缺失的近偏移距道数据；②需要估算声源子波。

反馈循环法的原理为，首先反演声波波场，将与自由界面有关的多次波预测出来，然后从原始波场中减去预测出的多次波(石颖etal.,2011)。反馈循环法是通过正演自由表面多次波波场将其削弱的，在此推导自由表面多次波压制的迭代算法。

由d(z0)＝d0(z0)+d0(z0)a(z0)d(z0)可得：

d0(z0)＝d(z0)-d0(z0)a(z0)d(z0)(13)

由式(13)结合d⁽⁰⁾＝d0得到迭代形式为：

式(14)为自由表面多次波滤波的迭代反演公式，n为迭代次数，d(z0)是原始浅剖数据(包含多次波)，d0⁽ⁿ⁾(z0)、d0⁽ⁿ⁺¹⁾(z0)分别为第n次和第n+1次迭代后除去多次波的浅剖数据估计，每次迭代过程，可取不同的自由界面算子a⁽ⁿ⁺¹⁾(z0)。

(4)综合压制方法

1)频域时不变双曲radon变换

由于浅剖的反射波为曲线且需要保证高分辨率高精度，因此可考虑运用双曲radon变换。再者，由于在时域进行双曲radon变换非常耗时，而且实用性不强，因此可考虑在频域上进行转换。时不变双曲radon变换可用下式表示(贾春梅etal.,2016)：

式(15)中，d(x,t)为时域浅剖数据，h为射线参数，τ为截距时间，zref为参考深度，m(τ,h)为双曲radon变换正剖面。

在式(15)两端分别进行傅里叶变换，可得：

其中θ可表示为：

针对浅剖单频分量，可用以下两个式子表示双曲radon正反变换：

m＝ld(18)

d＝l^hm(19)

上式中，l^h为l的共轭转置矩阵，l可表示为下式：

针对超定方程，式(18)可表示为以下最小平方形式：

m＝(l·l^h+λi)^-1ld(21)

当射线参数为等间隔采样时，l·l^h为toeplitz矩阵，可利用levinson递推法进行求解，具有较高的计算效率。

2)自适应滤波器

由于预测反褶积针对原始浅剖数据中长周期的多次波滤除效果不显著，经过预测反褶积处理后，有效反射波中仍然留下大部分的长周期多次波，因此需要设计一个自适应滤波器，将真正的有效一次反射波与长周期多次波实现分离，去除有效波能量，保留多次波能量。自适应滤波器可表示为：

式(22)中，d(τ,h)表示经过预测反褶积处理后包含多次波的浅剖数据，m(τ,h)表示多次波模型数据的双曲radon变换，η表示滤波器的剔除参数，通常取海底反射系数的绝对值，m表示滤波的平滑控制参数，通常取10以内的偶数。当fm(τ,h)近似为0时，认为不存在多次波，而当fm(τ,h)近似为1时，认为存在多次波。

3)综合压制方法流程

双曲radon域综合预测反褶积法与反馈循环法压制多次波的流程如图2所示。首先将原始浅剖数据分别通过预测反褶积和反馈循环法进行预处理，利用反馈循环法在原始浅剖数据中粗略预测多次波的位置，作为多次波道数据模型，预测出来的多次波道数据模型无需与原始数据的多次波完全严格匹配。同时，结合预测反褶积的处理结果，由于预测误差中包含有效反射波和未彻底压制的长周期多次波，多次波能量主要分布在具有较大曲率参数的道上(曾忠玉,2013)。其次，将两者预测的多次波成分和预测误差进行双曲radon变换。然后通过设计的自适应滤波器，将有效反射波区域切除，保留多次波能量。当fm(τ,h)近似为0时，认为还存在有效反射波，则需返回再迭代处理，直至fm(τ,h)约等于1，结束迭代处理。之后，结合两种方法获得的多次波道数据模型，获得双曲radon域统一的多次波数据模型，再通过反双曲radon变换，便可获得滤波后的多次波数据模型，即保留了多次波同相轴存在的浅剖数据。最后从原始浅剖数据减去多次波数据即可完成多次波的综合压制过程。

由此可见，本发明中解决了下列技术问题：

(1)针对浅海地区受海面自由波和浅地层层间多次波影响较严重的现象，综合分析了预测反褶积法和反馈循环法的优缺点，利用两者的优势互补，提出综合两者压制多次波的模型构建方案。

(2)针对预测反褶积预测误差中包含有效反射波和未彻底压制的长周期多次波，从设计自适应滤波器，将有效反射波区域切除，保留多次波能量，最后从原始数据中减去多次波能量数据获得高可靠性的浅地层剖面数据。

另外，本发明的技术特点如下：

(1)在充分分析浅海地区浅地层剖面数据特点，以及预测反褶积法和反馈循环法优缺点的基础上，构建综合预测反褶积法和反馈循环法压制多次波的模型。

(2)由于预测反褶积预测误差中包含有效反射波和未彻底压制的长周期多次波，从反面角度出发，经过多次迭代过程彻底消除有效反射波，进而获得成分纯度较高的多次波能量数据，最终从原始数据减去多次波能量数据得到高可靠性的浅地层剖面数据。

以上显示和描述了本发明的基本原理、主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等同物界定。

参考文献：

a.j.berkhout,1999.multipleremovalbasedonthefeedbackmodel.leadingedge.18,127-131.

a.j.berkhout,d.j.verschuur,1997.estimationofmultiplescatteringbyiterativeinversion,parti:theoreticalconsiderations.geophysics.62,1586-1595.

d.j.foster,c.c.mosher,1992.suppressionofmultiplereflectionsusingtheradontransform.geophysics.57,386-395.

p.g.kelamis,d.j.verschuur,2000.surface-relatedmultipleeliminationonlandseismicdata—strategiesviacasestudies.geophysics.65,719-734.

d.lokshtanov,1995.multiplesuppressionbysinglechannelandmultichanneldeconvolutioninthetau-pdomain.segtechnicalprogramexpandedabstracts.14,1482-1485.

d.lokshtanov,1999.multiplesuppressionbydata‐consistentdeconvolution.segtechnicalprogramexpandedabstracts.18,115-119.

l.morley,j.claerbout,1983.predictivedeconvolutioninshot‐receiverspace.geophysics.48,515-531.g.j.a.vangroenestijn,d.j.verschuur,2009a.estimatingprimariesbysparseinversionandapplicationtonear-offsetdatareconstruction.geophysics.74,a23-a28.

g.j.a.vangroenestijn,d.j.verschuur,2009b.estimationofprimariesandnear-offsetreconstructionbysparseinversion:marinedataapplications.geophysics.74,r119-r128.

d.j.verschuur,1992.adaptivesurface-relatedmultipleelimination.geophysics.57,1166-1177.

d.j.verschuur,a.j.berkhout,1997.estimationofmultiplescatteringbyiterativeinversion,partii:practicalaspectsandexamples.geophysics.62,1596-1611.

d.j.verschuur,a.j.berkhout,p.g.kelmis,1995.estimationofmultiplescatteringbyiterativeinversion,partii:examplesonmarineandlanddata.segtechnicalprogramexpandedabstracts.14,1470-1473.

y.h.wang,2004.multiplepredictionthroughinversion:afullydata-drivenconceptforsurface-relatedmultipleattenuation.geophysics.69,547-553.

y.h.wang,2007.multiplepredictionthroughinversion:theoreticaladvancementsandrealdataapplication.geophysics.72,v33-v39.

a.b.weglein,f.a.gasparotto,p.m.carvalho,1997.aninverse-scatteringseriesmethodforattenuatingmultiplesinseismicreflectiondata.geophysics.62,1975-1989.

o.yilmaz,m.t.taner,1994.discreteplane-wavedecompositionbyleast-mean-square-errormethod.geophysics.59,973-982.

曾忠玉,2013.复杂介质三维道集模型的纵波多分量nmo动校正速度研究,中国地震局地球物理研究所.巩向博,韩立国,李洪建,2014.各向异性radon变换及其在多次波压制中的应用.地球物理学报.57,2928-2936.

巩向博,韩立国,王升超,2016.混合域高分辨率双曲radon变换及其在多次波压制中的应用.石油物探.55,711-718.

贾春梅,姜国庆,刘志成,许璐,2016.频域稀疏双曲radon变换去噪方法.物探与化探.40,527-533.

刘军,2008.基于非高斯性最大化的预测反褶积技术,清华大学.

沈操,2003.基于波动方程的自由界面多次波压制,中国地质大学(北京).

石颖,刘洪,历玉英,2011.逆数据域表面多次波压制方法.吉林大学学报:地球科学版.41,271-276.

宋家文,d.j.verschuur,陈小宏,songjia-wen,d.j.verschuur,chenxiao-hong,2014.多次波压制的研究现状与进展.地球物理学进展.29,240-247.

王兆旗,2007.预测反褶积及其相关问题分析,长安大学.

薛亚茹,陈小宏,马继涛,2012.多方向正交多项式变换压制多次波.地球物理学报.55,3450-3458.

赵昌垒,叶月明,姚根顺,胡冰,庄锡进,2013.线性拉东域预测反褶积在海洋多次波去除中的应用.地球物理学进展.28,1026-1032.