基于三维稀疏聚焦域反演的噪声型数据一次波反演方法与流程

本发明属于地球物理勘探技术领域，具体涉及一种噪声型被动源数据的一次波反演方法，特别涉及一种基于三维稀疏聚焦域反演的噪声型数据一次波反演方法。

背景技术：

在地球物理勘探中，被动源地震勘探是一个新兴领域。它的提出与应用极大地的丰富了人们对地震波传播规律及性质的认识和研究，被认为是地球物理勘探的一大进步。不同于人工激发震源的主动源勘探，被动源地震勘探不需要人工震源，而是将地下的微振动或背景噪声作为震源。因此，被动源地震勘探能够极大的节约生产成本。

噪声型被动源地震探勘是被动源地震勘探中极为重要的一种勘探方式，通常背景噪声被视为噪声型被动源数据的震源。通过将检波器放置在地表，接收来自地下的波场从而形成被动源数据地震记录。

针对被动源数据地震记录，互相关算法是最常用的一种处理被动源数据的算法。该算法能够将被动源数据进行转化，形成地表激发，地表接收的虚拟炮集记录，从而获得地下的构造信息。同时，该虚拟炮集同传统主动源地震记录一样，不仅包含一次波信息也包含表层相关多次波信息。但是，这种方法依赖于地下均匀照明的假设，在激发角度和强度方面需要具有一致性。在实际采集中，这种假设很难实现，因此互相关算法很难提供具有准确振幅的近地表响应，所以获得结果也无法为后续的地震数据处理解释以及地下成像提供准确的地震资料。

为了克服上述问题，vangroenestijn等人(2010.estimationofprimariesbysparseinversionfrompassiveseismicdata.geophysics75(4).sa61–sa69)改进了针对传统主动源地震记录的一次波估计方法epsi(estimationofprimariesbysparseinversion)，使其能够使用于被动源地震数据多次波的压制。改进后的epsi不需要满足地下均匀照明的假设就能获得具有准确振幅的一次波响应。然而，该方法需要人为设置时窗，事先设置准确的一次波响应初始值。程浩等人(2015.estimatingprimariesfrompassiveseismicdata.explorationgeophysics46(2).184-191)提出将l1范数稀疏反演与epsi相结合，避免了传统epsi需要人为设置时窗的缺点，同时也提高了一次波响应估计结果的精度。然而，现有的这些方法基本只适用于脉冲型被动源数据，同时也无法解决被动源数据远偏移距振幅不准的问题。gabriel等人(2015.closed-loopsurface-relatedmultipleeliminationanditsapplicationtosimultaneousdatareconstruction.geophysics80(6).v189-v199)提出closed-loopsrme，该方法用一种更为稳定的目标函数进行主动源地震记录一次波的估计，但是该方法并不能直接应用于噪声型被动源地震数据的处理。

噪声型被动源数据中存在大量的相干噪声，会对一次波反演产生严重的干扰，因而这些方法如果直接应用在噪声型被动源数据上难以获得理想的效果。因此，针对噪声型被动源数据一次波估计的研究也鲜有报道。

技术实现要素：

本发明的目的就在于针对上述现有技术的不足，对closed-loopsrme及其相关方法进行大量针对性的改进，提出一种基于三维稀疏聚焦域反演的噪声型数据一次波反演方法。

本发明的目的是通过以下技术方案实现的：

首先，对噪声型被动源数据的特征进行分析，通过与主动源数据进行对比，总结相应的特征，针对噪声型被动源数据的具体特征提出相应的解决方案；然后，噪声型被动源数据的相干噪声和常规噪声有所不同，其形态和地震记录中的同相轴颇为类似，因此很难使用常规的去噪方法进行去噪，为此我们考虑将聚焦变换应用到噪声型被动源数据的处理中；最后对closed-loopsrme方法的目标函数进行改进并引入三维稀疏变换作为约束条件，同时将聚焦变换作为算子的形式引入到噪声型被动源数据多次波压制的过程中去，从而构建噪声型被动源数据一次波估计以及噪声压制的一体化进程。

一种基于三维稀疏聚焦域反演的噪声型数据一次波反演方法，包括以下步骤：

a、输入噪声型被动源数据；

b、根据地下全波场、一次波和多次波的物理关系，用线性算子相乘的形式表示地下波场传播的物理过程

式中，l为线性算子的集合，pnoise和pnoise，0分别代表输入波场列向量和噪声型被动源数据的一次波波场列向量；代表矩阵和向量的kronecker乘积；blockdiag代表在频率域生成对角线的块状矩阵；ft和分别代表能够实现时频域互相转换的傅里叶正反变换算子；i代表单位算子；a代表表层相关算子；

c、随后，引入l1范数构建噪声型被动源数据一次波估计反演体系，避免被动源数据一次波估计过程陷入局部最小值；

pnoise，0←argmin||pnoise，0||1s.t.||pnoise-lpnoise，0||2≤σ

通过谱梯度投影法(spgl1)求解上述最优化方程，得到一个初始的噪声型被动源数据一次波估计值，σ为输入数据的l2范数；

d、将一次波估计值pnoise，0带入下式

a＝argmin||p-pnoiss，0(i+ap)||2

求取第一次迭代过程中表层相关算子的估计值；

e、构建三维稀疏变换算子s以及聚焦去噪算子d,然后将算子结合为

其中，s是三维稀疏变换算子(2维curvelet变换以及1维wavelet变换的组合)，d是聚焦去噪算子，通过该算子能够实现噪声型被动源数据在聚焦域的去噪和远偏移距处的波形重构；

f、将构建好的三维稀疏聚焦算子t引入到反演过程中。将初始的一次波估计值pnoise，0置零，同时将求取的表层相关算子a的估计值带入到下式中去使用spgl1算法进行求解，得到噪声型被动源数据一次波估计值

g、重复步骤d，得到下一次迭代的表层相关算子的估计值；

h、将新的表层算子的估计值带入到步骤f中，直到得到的残差规定值，所得的一次波估计结果即为噪声型被动源数据的一次波估计值。

与现有技术相比，本发明的有益效果在于：

本发明将基于三维l1范数稀疏约束的closed-loopsrme方法引入到噪声型被动源数据一次波估计过程中去。同时，针对噪声型被动源数据的特征进行了大量针对性的改进，在反演过程中引入了聚焦变换进行噪声型被动源数据进行噪声压制和同相轴恢复，显著的提高了一次波估计结果的信噪比和精度。

本发明有以下特点：

1、基于三维稀疏聚焦域反演的噪声型数据一次波反演方法是一种数据驱动的一次波估计方法，无需知道地下模型先验信息，仅通过输入数据就能够实现一次波估计的目的。

2、在一次波估计过程引入了聚焦变换，在压制噪声方面有极好的效果，信噪比有了很大的提高。同时，借助于聚焦变换的作用，远偏移距处的同相轴信息得到正确的恢复，有利于后续的地震数据处理与解释。

3、本方法中的三维稀疏变换和聚焦变换的组合，能够结合稀疏性和聚焦性的特点，通过将输入数据转换至稀疏聚焦域，使信息得到更好的表示，进一步提高了一次波估计结果的精度。

4、基于三维稀疏聚焦域反演的噪声型数据一次波反演方法与常规的基于三维l1范数稀疏变换的closed-loopsrme方法不同，能够更好地适应输入数据中含有噪声的情况，尤其是含有这种和同相轴类似的相干噪声的情况。因而在噪声型被动源数据中有良好的表现，能够获得信噪比高的、更加精确的噪声型被动源数据一次波估计结果。

附图说明

图1基于三维稀疏聚焦域反演的噪声型数据一次波反演方法流程图；

图2地下全波场、一次波波场和多次波波场的地下路径物理传播关系图；

图3l1范数求解图；

图4对表层相关算子进行时间域稀疏约束；

图5a普通含噪声地震记录；

图5b对图5a做聚焦变换；

图5c在聚焦域进行噪声切除；

图5d将切除后的结果从聚焦域变换时空域；

图6a简单模型数值算例中的模型；

图6b噪声震源分布范围；

图6c噪声震源形态；

图6d原始噪声型被动源输入数据；

图6e主动源地震数据原始记录；

图6f第一次迭代获得的噪声型被动源数据一次波估计结果；

图6g第一次迭代获得的表层算子估计值；

图6h、图6i、图6j和图6k每次迭代获得的噪声型被动源数据一次波估计值；图6l、图6m、图6n和图6o每次迭代获得的表层算子估计值；

图6p最终获得的噪声型被动源数据一次波估计结果；

图6q对主动源地震记录估计得到主动源一次波结果；

图7a盐丘模型噪声型被动源数据数值算例中的模型；

图7b震源分布范围；

图7c震源形态；

图7d噪声型被动源数据输入记录；

图7e主动源地震数据原始记录；

图7f第一次迭代获得的噪声型被动源数据一次波估计结果；

图7g第一次迭代获得的表层相关算子估计值；

图7h最终的噪声型被动源数据一次波估计结果；

图7i真实的一次波；

图7j输入噪声型被动源数据零偏移距剖面；

图7k噪声型被动源数据一次波估计结果零偏移距剖面；

图7l主动源数据零偏移距剖面；

图7m主动源数据真实一次波结果零偏移距剖面。

具体实施方式

下面结合附图和实例对本发明进一步的详细说明。

本发明利用在大规模尺度反演的过程中引入三维稀疏变换来估计一次波，同时，在反演过程中利用聚焦变换压制噪声型被动源数据的噪声并重构其远偏移距的振幅。

基于三维稀疏聚焦域反演的噪声型数据一次波反演方法，包括以下步骤：

a、输入噪声型被动源数据原始记录。

b、根据地下全波场、一次波和多次波的物理关系，用线性算子相乘的形式

表示地下波场传播的物理过程

式中，pnoise和pnoise，0分别代表输入波场列向量和噪声型被动源数的一次波波场列向量；代表矩阵和向量的kronecker乘积；blockdiag代表在频率域生成对角线的块状矩阵；ft和分别代表能够实现时频互相转换的傅里叶正反变换算子；i代表单位矩阵算子；a代表表层相关算子。

对于噪声型被动源数据地下波场的传播关系，我们采用线性算子相乘的形式来表示地下全波场、一次波波场和多次波波场的物理传播过程，如图2所示，二阶多次波可以由一次波和一阶多次波褶积生成，也就是说全体多次波波场可以表示为一次波波场和多次波波场相乘的形式，基于此我们构建出ax＝b的标准正演形式。为了便于计算，我们用线性算子实现在频率域波场作用的过程，通过ft实现时间域波场到频率域波场的转化，随后再将频率域获得的结果通过转换回时间域，极大的提高了反演过程的稳定性和简便性。

c、引入l1范数至目标函数反演过程，避免被动源数据一次波估计过程陷入局部最小值。引入l1范数来该反演过程可以理解为：在反演过程中引入一次波估计结果的l1范数作为约束条件，可以稳定反演过程，控制求解方向。同时，由于l1范数集合和目标函数构成的集合都是凸集合。因此，对求解过程只会有全局最值，如图3所示，这样可以避免在一次波估计过程中陷入局部最小值

pnoise，0←argmin||pnoise，0||1s.t.||pnoise-lpnoise，0||2≤σ

通过基于l1范数的谱梯度投影法求解上述最优化方程，得到一个初始的噪声型被动源数据一次波估计值pnoise，0。

d、将第一次迭代获得的一次波估计值pnoise，0带入下式

a＝argmin||p-pnoise，0(i+ap)||2

求取第一次迭代过程中表层相关算子的估计值；只有获得一个正确的表层相关算子才能正确的表示地下传播过程中全波场，一次波波场，多次波波场之间的关系。因为如果只是一次波波场和全波场进行褶积，得到的结果在振幅和相位方面都是不准的，它能够纠正由于褶积产生的振幅和相位的差异。因此，通过在时间域对表层算子进行稀疏约束，如图4所示，选取与一次波对应的同相轴是至关重要。

e、构建三维稀疏变换算子s以及聚焦去噪算子d,然后将算子结合为

这里的三维稀疏变换(s)采用二维curvelet变换和一维小波变换结合的形式，将二维curvelet变换应用到炮-检波器域，同时在时间域上应用一维小波变换，这样避免了三维curvelet变换耗时过长的缺点。聚焦变换d的引入，可以压制数据中的噪声。通过图5a-图5d可以清楚的看到聚焦变换对噪声的压制作用，覆盖的数据(图5a)变换至聚焦域(图5b)然后切除(图5c)最终变换回时间域(图5d)，证明噪声型被动源数据的噪声完全可以通过聚焦变换的引入来实现。

f、将构建好的三维稀疏聚焦算子t引入到反演过程中。将初始的一次波估计值pnoise，0置零,同时将求取的表层相关算子a的估计值带入到下式中去

使用spgl1算法进行求解，得到噪声型被动源数据一次波估计值。σ值的选择是至关重要的，一般选择为输入数据l2范数的1％～10％左右，如果σ过大，就会产生导致反演过程不足，就无法得到完整的一次波估计结果；如果σ值给的过小，就会在反演过程中引入假象，这会极大地降低反演过程的准确度。

g、重复步骤d。

h、重复步骤f，直到满足要求的σ值，最终获得的一次波估计值即为最终的噪声型被动源数据一次波估计结果。

本发明基于三维稀疏聚焦域反演的噪声型数据一次波反演方法是通过matlab的平台实现的。

实施例1

基于三维稀疏聚焦域反演的噪声型数据一次波反演方法，包括以下步骤：

a、输入噪声型被动源数据记录(图6d)；该数据是由5000*1800的二维模型(图6a)正演产生，网格间距5m，地表均匀分布250个检波器，检波器间距为20m；99个噪声型震源(图6b)在地下1200m-1500m深度范围随机位置激发(图6c)，记录总时间800s。同时利用相同的模型正演主动源地震记录，本算例的目的是为了检测聚焦稀疏域能否在恢复远偏移距振幅的同时压制相干噪声；

b、根据地下全波场、一次波和多次波的物理关系，用线性算子相乘的形式

表示地下波场传播的物理过程

式中，pnoise和pnoiss，0分别代表输入波场数据列向量和噪声型被动源数据的一次波波场数据列向量；代表矩阵和向量的kronecker乘积；blockdiag代表在频率域生成对角线的块状矩阵；ft和分别代表能够实现傅里叶正反变换的傅里叶正反变换算子；i代表单位矩阵；a代表表层算子矩阵；

c、引入l1范数构建一次波反演过程，避免被动源数据一次波估计过程陷入局部最小值，同时能够稳定一次波估计过程。

pnoiss，0←argmin||pnoiss，0||1s.t.||pnoiss-lpnoiss，0||2≤σ

先将pnoiss，0和表层算子a初始值设为0，然后带入上式求解第一次获得pnoiss，0的初始迭代值(图6f)，σ取原数据l2范数的10％；

d、将第一次迭代获得的一次波估计值pnoiss，0带入下式

a＝argmin||p-p0(i+ap)||2

对得到的结果进行时间域的时窗稀疏获得准确的表层相关算子a(图6g)；

e、构建三维稀疏变换算子s以及聚焦去噪算子d,然后将算子结合为

这里的三维稀疏变换(s)采用二维curvelet变换和一维小波变换结合的形式，将二维curvelet变换应用到炮-检波器域，同时在时间域切片上应用一维小波变换，这样避免了三维curvelet变换耗时过长的缺点。聚焦变换d的引入，实现了在反演过程中去噪的目的；

f、将构建好的三维稀疏聚焦算子t引入到反演过程中，同时将初始的一次波估计值pnoise，0置零,同时将求取的表层相关算子a的估计值带入到下式中去

主循环中的σ选择为全数据l2范数的7％，然后使用spgl1算法进行求解，得到噪声型被动源数据一次波估计值；

g、将获得的被动源数据一次波估计值带入步骤d并重复步骤d。算法结后束会获得一系列的表层相关算子的估计值(图6l-图6o)，可以看到随着迭代次数的增加，表层相关算子的分辨率越来越高；

h、重复步骤f，获得一系列一次波估计值(图6h-图6k)，直到满足要求的σ值，即可得到最终的噪声型被动源数据一次波估计结果(图6p)。通过与主动源原始地震记录(图6e)以及相应的一次波记录(图6q)的波形进行对比，发现获得噪声型被动源数据的一次波结果成功的压制了噪声，恢复了远偏移距的振幅；同时由黑色箭头所示，多次波获得了极大的压制，一次波(白色箭头所示)获得了完整的保留，这对后续的地震成像和解释具有重大的意义。

实施例2

基于三维稀疏聚焦域反演的噪声型数据一次波反演方法，包括以下步骤：

a、输入噪声型被动源数据记录(图7d)；该数据是由5400m*1200大小的盐丘模型正演产生，网格间距5m，地表均匀分布250个检波器，检波器间距为20m；99个噪声型震源(图7b)在地下900m-100m深度范围随机位置激发(图7c)，记录总时间800s。该模型是多次波压制领域著名的模型(图7a)，多用于检测多次波压制的效果，所以使用该模型可以有效的检测本方法的合理性和准确性，主要用于验证该方法对一次波估计的精确度。因此，用同样的模型(图7a)来正演主动源地震记录(图7e)及真实一次波记(图7i)，用于评估噪声型被动源一次波估计效果。在(图7d)中，白色箭头所示为一次波同相轴，但应该注意到在远偏移距处一次波的振幅出现误差，这是由于非均匀照明导致的。图7d中黑色箭头所示为多次波；

b、根据地下全波场、一次波和多次波的物理关系，用线性算子相乘的形式表示地下波场传播的物理过程

式中，pnoiss和pnoiss，0分别代表输入波场数据列向量和噪声型被动源数据的一次波波场数据列向量；代表矩阵和向量的kronecker乘积；blockdiag代表在频率域生成对角线的块状矩阵；ft和分别代表能够实现傅里叶正反变换的傅里叶正反变换算子；i代表单位矩阵；a代表表层算子矩阵；

c、引入l1范数构建一波次反演过程，避免被动源数据一次波估计过程

陷入局部最小值的情况，同时能够稳定一次波估计过程。

pnoise，0←argmin||pnoise，0||1s.t.||pnoise-lpnoise，0||2≤σ

先将pnoiss，0和表层算子a初始值设为0，然后带入上式求解第一次获得pnoiss，0的初始迭代值，图7f所示。

d、将第一次迭代获得的一次波估计值pnoiss，0带入下式

a＝argmin||p-p0(i+ap)||2

获得表层相关算子估计值之后，再通过时间域的稀疏获取准确的表层相关算子(图7g)；

e、构建三维稀疏变换算子s以及聚焦去噪算子d,然后将算子结合为

这里的三维稀疏变换(s)采用二维curvelet变换和一维小波变换结合的形式，将二维curvelet变换应用到炮-检波器域，同时在时间域切片上应用一维小波变换；

f、将构建好的三维稀疏聚焦算子t引入到反演过程中，将初始的一次波估计

值pnoise，0置零,同时将求取的表层相关算子a的估计值带入到下式中去。

使用spgl1算法进行求解，得到噪声型被动源数据一次波估计值。σ值选择为输入数据l2范数的1％～10％左右；

g、重复步骤d；

h、重复步骤f，直到满足要求的σ值，即可得到最终的噪声型被动源数据一次波估计结果(图7h)。跟真实一次波(图7i)相比，尽管输入数据中有极大的噪声干扰，但是在0.5s-1.0s的一次波同相轴得到了正确的估计，同时位于远偏移距处的同相轴得到了正确的纠正。通过与主动源记录及其估计的一次波记录，发现在远偏移距的振幅以及传播时间获得了准确的重构；同时噪声型被动源记录中的相干噪声得到了极大的压制。在一次波估计方面，通过将输入数据零偏剖面(图7j)及其一次波零偏剖面(图7k)、主动源数据(图7l)及其一次波零偏剖面(图7m)对比，验证本方法能够实现噪声型被动源数据一次波估计、噪声压制和振幅恢复的一体化进程。