CN105259575A - 快速3d自由表面多次波预测方法 - Google Patents

Abstract

本发明属于地震勘探技术中地震信号处理领域，具体公开了一种快速3D自由表面多次波预测方法。该方法采用混合频率-时间域抛物线稀疏反演方法来估计横测线MCG的同相轴顶点位置，实现3D自由表面多次波预测。该方法对于每一个曲率参数，能沿着同相轴顶点的时间和空间方向施加稀疏约束。在已知道位置相同的横测线MCG中，基于快速迭代收缩阈值算法的混合频率-时间域抛物线稀疏反演方法对每一个频率值，只需计算一次矩阵求逆。相比于传统的基于频率域抛物线稀疏反演的3D自由表面多次波预测方法，本发明方法在降低计算量的同时，能进一步提高多次波预测的精度。

Description

快速3D自由表面多次波预测方法

技术领域

本发明属于地震勘探技术中地震信号处理领域，具体涉及一种基于混合频率-时间域抛物线稀疏反演的快速3D自由表面多次波预测方法。

背景技术

SRME(SurfaceRelatedMultipleElimination)广泛用来压制海洋地震勘探中的自由表面多次波。2DSRME并不能有效地适应地下介质在横测线方向存在倾角或拖缆采集存在羽角漂移的情况(Dragoset,B.,E.Verschuur,I.MooreandR.Bisley,2010,Aperspectiveon3Dsurface-relatedmultipleelimination:Geophysics,75,75A245-75A261.)。考虑到地震波场的3D传播特性，3DSRME比2DSRME能更好地预测自由表面多次波的旅行时。然而，直接应用3DSRME需要地震波场在横测线和纵测线方向有足够大的采样密度。实际3D地震采集系统在炮线之间或拖缆线之间有比较大的间距，这会限制3DSRME方法的直接和广泛应用。

目前，已经提出多个方法来解决3DSRME方法面临的数据稀疏采样问题。这些方法包括在应用3DSRME之前进行数据重建的方法、即时内插方法和横测线多次波贡献道集(MultipleContributionGather/MCG)重建方法等。基于双曲线稀疏反演的横测线MCG重建方法(vanDedem,E.J.,2002,3Dsurface-relatedmultipleprediction:Ph.D.thesis,DelftUniversityofTechnology.)假设横测线MCG中的同相轴满足双曲线时距曲线规律，每一个同相轴的顶点位置表征多次波的真实旅行时。双曲线稀疏反演方法通过对每一个曲率参数，沿着同相轴顶点的时间和空间方向施加稀疏约束，来估计同相轴顶点位置。后来，抛物线稀疏反演方法(Hokstad,K,.andR.Sollie,2006,3Dsurface-relatedmultipleeliminationusingparabolicsparseinversion:Geophysics,71,V145-V152.)被提出来，并假设横测线MCG中的同相轴满足抛物线时距曲线规律。抛物线稀疏反演方法在频率域进行Radon变换，具有比双曲线稀疏反演方法更高的计算效率。频率域抛物线稀疏反演方法将数据空间和模型空间的关系表达为：

$d(y_k,\mathrm{\ω})=\underset{u=1}{\overset{N_q}{\Σ}}\underset{v=1}{\overset{N_y}{\Σ}}m(q_u,y_v^0,\mathrm{\ω})e^{{\mathrm{iwq}}_u{\left(y_k-y_v^0\right)}^2},=1,2,...,N_d,---\left(1\right)$

或

d＝Lm，(2)

$d={\left[\begin{array}{cccc}d(y_1,\mathrm{\ω})& d(y_2,\mathrm{\ω})& ...& d(y_{N_d},\mathrm{\ω})\end{array}\right]}^T,---\left(3\right)$

$m={\left[\begin{array}{ccccccc}m(q_1,y_1^0,\mathrm{\ω})& ...& m(q_1,y_{Ny}^0,\mathrm{\ω})& ...& m(q_{Nq},y_1^0,\mathrm{\ω})& ...& m(q_{Nq},y_{Ny}^0,\mathrm{\ω})\end{array}\right]}^T,---\left(4\right)$

$L=\left[\begin{array}{ccccccc}{e}^{{\mathrm{iwq}}_1{\left(y_1-y_1^0\right)}^2}& \mathrm{...}& e^{{\mathrm{iwq}}_1{\left(y_1-y_{N_y}^0\right)}^2}& \mathrm{...}& e^{{\mathrm{iwq}}_{N_q}{\left(y_1-y_1^0\right)}^2}& \mathrm{...}& e^{{\mathrm{iwq}}_{N_q}{\left(y_1-y_{N_y}^0\right)}^2}\\ e^{{\mathrm{iwq}}_1{\left(y_2-y_1^0\right)}^2}& \mathrm{...}& e^{{\mathrm{iwq}}_1{\left(y_2-y_{N_y}^0\right)}^2}& \mathrm{...}& e^{{\mathrm{iwq}}_{N_q}{\left(y_2-y_1^0\right)}^2}& \mathrm{...}& e^{{\mathrm{iwq}}_{N_q}{\left(y_2-y_{N_y}^0\right)}^2}\\ e^{{\mathrm{iwq}}_1{\left(y_{N_d}-y_1^0\right)}^2}& \mathrm{...}& e^{{\mathrm{iwq}}_1{\left(y_{N_d}-y_{N_y}^0\right)}^2}& \mathrm{...}& e^{{\mathrm{iwq}}_{N_q}{\left(y_{N_d}-y_1^0\right)}^2}& \mathrm{...}& e^{{\mathrm{iwq}}_{N_q}{\left(y_{N_d}-y_{N_y}^0\right)}^2}\end{array}\right],---\left(5\right)$

其中，d表示数据空间(横测线MCG)，m表示模型空间，L表示Radon变换算子，N_q表示总的曲率参数个数，N_y表示总的顶点个数，N_d表示横测线MCG中总的地震道数。q_u表示曲率参数，表示顶点位置，y_k表示第k道的位置到相应同相轴顶点位置的偏移距。

频率域抛物线稀疏反演方法的优化问题为：

$\arg \underset{m}{min}\left|\right|d-Lm|{|}_2^2+\mathrm{\λ}|\left|m\right|{|}_1,---\left(6\right)$

其中，λ表示正则化因子。频率域抛物线稀疏反演方法并不能在时间方向施加稀疏约束，相比于时间域双曲线稀疏反演方法，会降低多次波预测的精度。

另外，频率域抛物线稀疏反演方法采用迭代重加权最小二乘算法来估计横测线MCG的同相轴顶点位置。迭代重加权最小二乘算法在预测每一道中的多次波时，需要在每一次迭代对于每一个频率值，进行一次矩阵-矩阵相乘和矩阵求逆运算，计算量较大。

发明内容

本发明的目的在于提出一种快速3D自由表面多次波预测方法，该方法利用基于快速迭代收缩阈值算法的混合频率-时间域抛物线稀疏反演方法来估计横测线MCG的同相轴顶点位置，并对同相轴顶点位置进行振幅和相位校正后，沿横向相加得到3D自由表面多次波的预测道，在所有横测线MCG处理完毕后得到所有道的3D自由表面多次波预测结果。

为了实现上述目的，本发明采用如下技术方案：

快速3D自由表面多次波预测方法，包括如下步骤：

a对3D地震数据缺失的近偏移距进行外推，并对缺失的地震炮集进行内插，然后构建相应的横测线MCG；对横测线MCG进行同相轴顶点估计的数学模型表示为：

$\stackrel{\‾}{d}=F^{-1}\[LF\[\stackrel{\‾}{m}\]\],$

其中，F表示沿时间方向的傅立叶正变换，F^-1表示相应的逆傅立叶变换，表示时间域的同相轴顶点位置道集，表示横测线MCG，L表示Radon变换算子；

b设置变量初始值；需要设置初始值的变量包括收缩阈值s_α、阻尼因子β、最大迭代次数N、迭代步长l_t、曲率参数个数N_q和大小、顶点参数个数N_y和大小；

c根据曲率参数个数N_q和大小、顶点参数个数N_y和大小构建Radon变换算子L，并计算矩阵 $\stackrel{\‾}{L}=L^H{({\mathrm{LL}}^H+\mathrm{\β}I)}^{-1};$

其中，L^H表示Radon变换算子L的共轭转置矩阵，I表示单位矩阵；

d利用步骤c得到的矩阵对横测线MCG进行同相轴顶点估计，并得到3D自由表面多次波的预测道M_p；

e判断横测线MCG是否全部处理完毕；如果否，返回步骤d；如果全部处理完毕，输出3D自由表面多次波的预测结果。

进一步，步骤d中3D自由表面多次波的预测道M_p的具体计算过程如下：

d1设置迭代数n＝0，利用步骤c得到的矩阵计算同相轴顶点位置的初始值 ${\stackrel{\‾}{m}}^{\left(0\right)}=F^{-1}\[\stackrel{\‾}{L}F\[\stackrel{\‾}{d}\]\];$

d2令n＝n+1,对上一步估计的顶点位置进行收缩阈值操作：

${\tilde{m}}^{\left(n\right)}=T_{\mathrm{\α}}\left({\stackrel{\‾}{m}}^{\left(n\right)}\right),$

其中，T_α为收缩阈值算子，定义为：

$T_{\mathrm{\α}}{\left(\stackrel{\‾}{m}\right)}_{i,j,k}={\left(\right|{\stackrel{\‾}{m}}_{i,j,k}|-s_{\mathrm{\α}}C)}_{+}\mathrm{sgn}\left({\stackrel{\‾}{m}}_{i,j,k}\right),$

i＝1,2,…,N_q,j＝1,2,…,N_y,k＝1,2,…,N_t

其中，表示第n步估计的顶点位置的收缩阈值结果，表示第n步估计的顶点位置，N_t表示时间采样点个数，表示向量的下标为(i,j,k)的元素，

${\left(a\right)}_{+}=\left\{\begin{array}{cc}a,& a\&GreaterEqual;0\\ 0,& a=0\end{array}\right.,\mathrm{sgn}\left(a\right)=\left\{\begin{array}{cc}1,& a>0\\ 0,& a=0\\ -1,& a<0\end{array}\right.;$

d3对收缩阈值结果进行更新：

${\tilde{y}}^{\left(n\right)}=\left\{\begin{array}{cc}{\tilde{m}}^{\left(1\right)},& n=1\\ {\tilde{m}}^{\left(n\right)}+\left(\right(t_s^{\left(n\right)}-1)/t_s^{(n+1)})({\tilde{m}}^{\left(n\right)}-{\tilde{m}}^{(n-1)}),& n>1\end{array}\right.,$

其中，表示更新后的收缩阈值结果，序列表达为

d4求取更新后的同相轴顶点位置：

${\stackrel{\&OverBar;}{m}}^{(n+1)}={\tilde{y}}^{\left(n\right)}+l_t{F}^{-1}\&lsqb;L^H{({\mathrm{LL}}^H+\mathrm{\&beta;}I)}^{-1}(F\&lsqb;\stackrel{\&OverBar;}{d}\&rsqb;-LF\&lsqb;{\tilde{y}}^{\left(n\right)}\&rsqb;)\&rsqb;;$

d5判断迭代次数n是否达到最大迭代次数N；如没达到，返回步骤d2；如果达到，利用同相轴顶点位置的估计结果，按下面公式得到3D自由表面多次波的预测道M_p：

$M_p=\underset{v=1}{\overset{N_y}{\&Sigma;}}\underset{u=1}{\overset{N_q}{\&Sigma;}}F(q_u,\mathrm{\&tau;})\stackrel{\&OverBar;}{m}(q_u,y_v),F(q_u,\mathrm{\&tau;})=\sqrt{\frac{2\mathrm{\&pi;}\mathrm{\&tau;}}{{\mathrm{\&omega;q}}_u}}{e}^{-i\mathrm{\&pi;}/4},$

其中，F(q_u,τ)表示振幅和相位校正因子，表示同相轴顶点位置的估计结果，q_u表示曲率参数，y_v表示顶点位置，τ表示旅行时，ω表示频率，i表示虚数单位。

本发明具有如下优点：

本发明方法采用混合频率-时间域抛物线稀疏反演方法来估计横测线MCG的同相轴顶点位置，实现3D自由表面多次波预测。该方法对于每一个曲率参数，能沿着同相轴顶点的时间和空间方向施加稀疏约束。在已知道位置相同的横测线MCG中，基于快速迭代收缩阈值算法的混合频率-时间域抛物线稀疏反演方法对每一个频率值，只需计算一次矩阵求逆。相比于传统的基于频率域抛物线稀疏反演的3D自由表面多次波预测方法，本发明方法在降低计算量的同时，能进一步提高多次波预测的精度。

附图说明

图1为本发明中快速3D自由表面多次波预测方法的流程示意图；

图2a为原始数据的共炮点道集图；

图2b为真实一次波的共炮点道集图；

图3a为横测线MCG示意图；

图3b为本发明中快速3D自由表面多次波预测方法估计的同相轴顶点示意图，其中，由图3a顶部的黑色方框标明的地震道作为输入数据；

图3c为传统频率域抛物线稀疏反演方法估计的同相轴顶点示意图，其中，由图3a顶部的黑色方框标明的地震道作为输入数据；

图4a为直接利用3DSRME方法得到的预测多次波示意图；

图4b为利用本发明快速3D自由表面多次波预测方法得到的预测多次波示意图；

图4c为传统的频率域抛物线稀疏反演方法得到的预测多次波示意图；

图4d为2DSRME方法得到的预测多次波示意图；

图5a为图2a中2050毫秒到2396毫秒的放大显示结果示意图；

图5b为图4a中2050毫秒到2396毫秒的放大显示结果示意图；

图5c为图4b中2050毫秒到2396毫秒的放大显示结果示意图；

图5d为图4c中2050毫秒到2396毫秒的放大显示结果示意图；

图5e为图4d中2050毫秒到2396毫秒的放大显示结果示意图；

图6a为利用图4a中的预测多次波进行自适应相减的结果示意图；

图6b为利用图4b中的预测多次波进行自适应相减的结果示意图；

图6c为利用图4c中的预测多次波进行自适应相减的结果示意图；

图6d为利用图4d中的预测多次波进行自适应相减的结果示意图。

具体实施方式

本发明的基本思想为：将横测线MCG的同相轴顶点估计问题表征为混合频率-时间域抛物线稀疏反演问题，构建相应的数学模型，在频率域进行Radon正反变换，并构建对每一个曲率参数，沿同相轴顶点的时间和空间方向施加稀疏约束的优化问题，同时引入快速迭代收缩阈值算法降低优化问题求解的计算复杂度。

具体地，下面结合附图以及具体实施方式对本发明作进一步详细说明：

结合图1所示，快速3D自由表面多次波预测方法，包括如下步骤：

a对3D地震数据缺失的近偏移距进行外推，并对缺失的地震炮集进行内插，然后构建相应的横测线MCG；本发明中方法的数学模型为：

$\stackrel{\&OverBar;}{d}=F^{-1}\&lsqb;LF\&lsqb;\stackrel{\&OverBar;}{m}\&rsqb;\&rsqb;,$

其中，F表示沿时间方向的傅立叶正变换，F^-1表示相应的逆傅立叶变换，表示时间域的同相轴顶点位置道集，表示横测线MCG，L表示Radon变换算子；

估计同相轴顶点的优化问题可以表示为：

$\arg \underset{\stackrel{\&OverBar;}{m}}{min}\left|\right|\stackrel{\&OverBar;}{d}-F^{-1}\&lsqb;LF\&lsqb;\stackrel{\&OverBar;}{m}\&rsqb;\&rsqb;|{|}_2^2+\mathrm{\&lambda;}|\left|\stackrel{\&OverBar;}{m}\right|{|}_1,$

其中，λ表示正则化因子。

为了求解上式中给定的优化问题，本发明方法采用快速迭代收缩阈值算法估计同相轴顶点位置实现3D自由表面多次波预测，即：

在每一步迭代，快速迭代收缩阈值算法采用收缩阈值算子对每一个曲率值，沿同相轴顶点的时间和空间方向促进所估计同相轴顶点位置的稀疏性，得到同相轴顶点位置的收缩阈值结果，并将当前步和前一步的收缩阈值结果进行线性组合来对当前步的收缩阈值结果进行更新，然后采用更新后的收缩阈值结果来对同相轴顶点位置进行更新。

b设置变量初始值；需要设置初始值的变量包括收缩阈值s_α、阻尼因子β、最大迭代次数N、迭代步长l_t、曲率参数个数N_q和大小、顶点参数个数N_y和大小；

c根据曲率参数个数N_q和大小、顶点参数个数N_y和大小构建Radon变换算子L，并计算矩阵 $\stackrel{\&OverBar;}{L}=L^H{({\mathrm{LL}}^H+\mathrm{\&beta;}I)}^{-1};$

其中，L^H表示Radon变换算子L的共轭转置矩阵，I表示单位矩阵；

d利用步骤c得到的矩阵对横测线MCG进行同相轴顶点估计，并得到3D自由表面多次波的预测道M_p；具体计算过程如下：

d1设置迭代数n＝0，利用步骤c得到的矩阵计算同相轴顶点位置的初始值 ${\stackrel{\&OverBar;}{m}}^{\left(0\right)}=F^{-1}\&lsqb;\stackrel{\&OverBar;}{L}F\&lsqb;\stackrel{\&OverBar;}{d}\&rsqb;\&rsqb;;$

d2令n＝n+1,对上一步估计的顶点位置进行收缩阈值操作：

${\tilde{m}}^{\left(n\right)}=T_{\mathrm{\&alpha;}}\left({\stackrel{\&OverBar;}{m}}^{\left(n\right)}\right),$

其中，T_α为收缩阈值算子，定义为：

$T_{\mathrm{\&alpha;}}{\left(\stackrel{\&OverBar;}{m}\right)}_{i,j,k}={\left(\right|{\stackrel{\&OverBar;}{m}}_{i,j,k}|-s_{\mathrm{\&alpha;}}C)}_{+}\mathrm{sgn}\left({\stackrel{\&OverBar;}{m}}_{i,j,k}\right),$

i＝1,2,…,N_q,j＝1,2,…,N_y,k＝1,2,…,N_t

其中，表示第n步估计的顶点位置的收缩阈值结果，表示第n步估计的顶点位置，N_t表示时间采样点个数，表示向量的下标为(i,j,k)的元素，

${\left(a\right)}_{+}=\left\{\begin{array}{cc}a,& a\&GreaterEqual;0\\ 0,& a=0\end{array}\right.,\mathrm{sgn}\left(a\right)=\left\{\begin{array}{cc}1,& a>0\\ 0,& a=0\\ -1,& a<0\end{array}\right.;$

d3对收缩阈值结果进行更新：

${\tilde{y}}^{\left(n\right)}=\left\{\begin{array}{cc}{\tilde{m}}^{\left(1\right)},& n=1\\ {\tilde{m}}^{\left(n\right)}+\left(\right(t_s^{\left(n\right)}-1)/t_s^{(n+1)})({\tilde{m}}^{\left(n\right)}-{\tilde{m}}^{(n-1)}),& n>1\end{array}\right.,$

其中，表示更新后的收缩阈值结果，序列表达为

d4求取更新后的同相轴顶点位置：

${\stackrel{\&OverBar;}{m}}^{(n+1)}={\tilde{y}}^{\left(n\right)}+l_t{F}^{-1}\&lsqb;L^H{({\mathrm{LL}}^H+\mathrm{\&beta;}I)}^{-1}(F\&lsqb;\stackrel{\&OverBar;}{d}\&rsqb;-LF\&lsqb;{\tilde{y}}^{\left(n\right)}\&rsqb;)\&rsqb;;$

d5判断迭代次数n是否达到最大迭代次数N；如没达到，返回步骤d2；如果达到，利用同相轴顶点位置的估计结果，按下面公式得到3D自由表面多次波的预测道M_p：

$M_p=\underset{v=1}{\overset{N_y}{\&Sigma;}}\underset{u=1}{\overset{N_q}{\&Sigma;}}F(q_u,\mathrm{\&tau;})\stackrel{\&OverBar;}{m}(q_u,y_v),F(q_u,\mathrm{\&tau;})=\sqrt{\frac{2\mathrm{\&pi;}\mathrm{\&tau;}}{{\mathrm{\&omega;q}}_u}}{e}^{-i\mathrm{\&pi;}/4},$

其中，F(q_u,τ)表示振幅和相位校正因子，表示同相轴顶点位置的估计结果，q_u表示曲率参数，y_v表示顶点位置，τ表示旅行时，ω表示频率，i表示虚数单位；

e判断横测线MCG是否全部处理完毕；如果否，返回步骤d；如果全部处理完毕，输出3D自由表面多次波的预测结果。

在仿真实验中，利用3D模型数据验证本发明中方法的有效性：

图2a为原始数据的共炮点道集图，图2b为真实一次波的共炮点道集图。图2a和图2b中横坐标InlineCoordinate表示纵测线坐标，单位为米(m)，纵坐标Time表示时间，单位为毫秒(ms)。为预测图2a中第6道的3D自由表面多次波，图3a给出了相应的横测线MCG。图3a中，横坐标CrosslineCoordinate表示横测线坐标，单位为米(m)，纵坐标Time表示时间，单位为毫秒(ms)。在图3a顶部以黑色方框标识的地震道作为稀疏反演方法的输入数据。本发明方法和频率域抛物线稀疏反演方法得到的同相轴顶点估计结果显示在图3b和图3c中。图3b和图3c中横坐标ApexPosition表示顶点位置，单位为米(m)，纵坐标Time表示时间，单位为毫秒(ms)。从图3b和图3c可以看出，本发明方法比传统的频率域抛物线稀疏反演方法得到时间分辨率更高的同相轴顶点估计结果。

直接应用3DSRME方法处理不存在稀疏采样问题的3D模型数据，得到图4a中的3D预测多次波。横测线间距为100米的3D数据作为本发明方法和传统频率域抛物线稀疏反演方法的输入数据，分别得到图4b和图4c中的3D预测多次波。2DSRME方法得到的2D预测多次波显示在图4d中。图4a至图4d中横坐标InlineCoordinate表示纵测线坐标，单位为米(m)，纵坐标Time表示时间，单位为毫秒(ms)。将图2a以及图4a至图4d中从2050毫秒到2396毫秒的部分显示在图5a至图5e中。图5a至图5e中横坐标InlineCoordinate表示纵测线坐标，单位为米(m)，纵坐标Time表示时间，单位为毫秒(ms)。从图5a至图5e可以看出，3D预测多次波的旅行时比2D预测多次波的旅行时更加准确。将图4a至图4d中的预测多次波从原始数据中自适应减去，得到的一次波估计结果分别显示在图6a至图6d中。图6a至图6d中横坐标InlineCoordinate表示纵测线坐标，单位为米(m)，纵坐标Time表示时间，单位为毫秒(ms)。由于2DSRME方法并不能有效地预测多次波的旅行时，因此图6d中存在明显的残余多次波(由白色箭头标明)。在图4c中，频率域抛物线稀疏反演方法并不能有效地预测黑色箭头标明的多次波的振幅。因此，图6c中存在明显的残余多次波(由黑色箭头标明)。对于图6a至图6d中的一次波，相应的信噪比为31.68、31.35、30.27和29.51。利用本发明方法得到的3D自由表面多次波进行自适应相减，比利用传统的频率域抛物线稀疏反演方法得到的3D自由表面多次波进行自适应相减，得到更高的信噪比。对于图4b和图4c中预测的3D自由表面多次波，相应的计算时间为82.3秒和185.2秒。本发明方法比传统的频率域抛物线稀疏反演方法具有更高的计算效率。

当然，以上说明仅仅为本发明的较佳实施例，本发明并不限于列举上述实施例，应当说明的是，任何熟悉本领域的技术人员在本说明书的教导下，所做出的所有等同替代、明显变形形式，均落在本说明书的实质范围之内，理应受到本发明的保护。

Claims (2)

1.快速3D自由表面多次波预测方法，其特征在于，包括如下步骤：

a对3D地震数据缺失的近偏移距进行外推，并对缺失的地震炮集进行内插，然后构建相应的横测线MCG；对横测线MCG进行同相轴顶点估计的数学模型表示为：

$\stackrel{\&OverBar;}{d}=F^{-1}\&lsqb;LF\&lsqb;\stackrel{\&OverBar;}{m}\&rsqb;\&rsqb;,$

其中，F表示沿时间方向的傅立叶正变换，F^-1表示相应的逆傅立叶变换，表示时间域的同相轴顶点位置道集，表示横测线MCG，L表示Radon变换算子；

b设置变量初始值；需要设置初始值的变量包括收缩阈值s_α、阻尼因子β、最大迭代次数N、迭代步长l_t、曲率参数个数N_q和大小、顶点参数个数N_y和大小；

c根据曲率参数个数N_q和大小、顶点参数个数N_y和大小构建Radon变换算子L，并计算矩阵 $\stackrel{\&OverBar;}{L}=L^H{({\mathrm{LL}}^H+\mathrm{\&beta;}I)}^{-1};$

其中，L^H表示Radon变换算子L的共轭转置矩阵，I表示单位矩阵；

d利用步骤c得到的矩阵对横测线MCG进行同相轴顶点估计，并得到3D自由表面多次波的预测道M_p；

e判断横测线MCG是否全部处理完毕；如果否，返回步骤d；如果全部处理完毕，输出3D自由表面多次波的预测结果。

2.根据权利要求1所述的快速3D自由表面多次波预测方法，其特征在于，步骤d中3D自由表面多次波的预测道M_p的具体计算过程如下：

d1设置迭代数n＝0，利用步骤c得到的矩阵计算同相轴顶点位置的初始值 ${\stackrel{\&OverBar;}{m}}^{\left(0\right)}=F^{-1}\&lsqb;\stackrel{\&OverBar;}{L}F\&lsqb;\stackrel{\&OverBar;}{d}\&rsqb;\&rsqb;;$

d2令n＝n+1,对上一步估计的顶点位置进行收缩阈值操作：

${\tilde{m}}^{\left(n\right)}=T_{\mathrm{\&alpha;}}\left({\stackrel{\&OverBar;}{m}}^{\left(n\right)}\right),$

其中，T_α为收缩阈值算子，定义为：

$T_{\mathrm{\&alpha;}}{\left(\stackrel{\&OverBar;}{m}\right)}_{i,j,k}={\left(\right|{\stackrel{\&OverBar;}{m}}_{i,j,k}|-s_{\mathrm{\&alpha;}}C)}_{+}\mathrm{sgn}\left({\stackrel{\&OverBar;}{m}}_{i,j,k}\right),$

i=1，2，…，N_q，j=1，2，…，N_y，k=1，2，…，N_t

其中，表示第n步估计的顶点位置的收缩阈值结果，表示第n步估计的顶点位置，N_t表示时间采样点个数，表示向量的下标为(i,j,k)的元素， $\begin{array}{cc}{\left(a\right)}_{+}=\{\begin{array}{cc}a,& a\&GreaterEqual;0\\ 0,& a=0\end{array},& \mathrm{sgn}\left(a\right)=\left\{\begin{array}{cc}1,& a>0\\ 0,& a=0\\ -1,& a<0\end{array}\right.\end{array};$

d3对收缩阈值结果进行更新：

${\tilde{y}}^{\left(n\right)}=\{\begin{array}{cc}{\tilde{m}}^{\left(1\right)},& n=1\\ {\tilde{m}}^{\left(n\right)}+\left(\right(t_s^{\left(n\right)}-1)/t_s^{(n+1)})({\tilde{m}}^{\left(n\right)}-{\tilde{m}}^{(n-1)}),& n>1\end{array},$

其中，表示更新后的收缩阈值结果，序列表达为

d4求取更新后的同相轴顶点位置：

${\stackrel{\&OverBar;}{m}}^{(n+1)}={\tilde{y}}^{\left(n\right)}+l_t{F}^{-1}\&lsqb;L^H{({\mathrm{LL}}^H+\mathrm{\&beta;}I)}^{-1}(F\&lsqb;\stackrel{\&OverBar;}{d}\&rsqb;-LF\&lsqb;{\tilde{y}}^{\left(n\right)}\&rsqb;)\&rsqb;;$

d5判断迭代次数n是否达到最大迭代次数N；如没达到，返回步骤d2；如果达到，利用同相轴顶点位置的估计结果，按下面公式得到3D自由表面多次波的预测道M_p：

$\begin{array}{cc}{M}_p=\underset{v=1}{\overset{N_y}{\&Sigma;}}\underset{u=1}{\overset{N_q}{\&Sigma;}}F(q_u,\mathrm{\&tau;})\stackrel{\&OverBar;}{m}(q_u,y_v),& F(q_u,\mathrm{\&tau;})=\sqrt{\frac{2\mathrm{\&pi;}\mathrm{\&tau;}}{{\mathrm{\&omega;q}}_u}}{e}^{-i\mathrm{\&pi;}/4}\end{array},$

其中，F(q_u,τ)表示振幅和相位校正因子，表示同相轴顶点位置的估计结果，q_u表示曲率参数，y_v表示顶点位置，τ表示旅行时，ω表示频率，i表示虚数单位。