具体实施方式

下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本申请保护的范围。

为使本申请的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本申请作进一步详细的说明。

实施例一

如图1所示，为本申请实施例一中的基于照明方向约束的被动源干涉偏移成像方法流程示意图，主要由被动源干涉、照明信息提取、照明方向约束偏移成像三步骤构成。下面详细介绍本申请实施例一的三个步骤。

第一步，对单点脉冲型被动源的地震记录进行干涉处理，得到被动源的反褶积函数，通过反褶积函数，可以获得不准确的干涉虚炮集记录。

对于脉冲型被动源，Wapenaar&Fokkema(2006)介绍了在声波介质中，在被动源连续均匀分布(各自位于x_s处，在本实施例一中标记x_s表征被动源)的情况下，重构两个检波器x_A、x_B(其中，x_A表征虚震源)之间格林函数的关系式。本申请实施例一参考该方法，当两个检波器位于自由表面记录被动源地震时，重构出的两个检波器间的格林函数即为地下界面的反射响应：

其中R代表实部，{}^*代表复共轭，ω是角频率，^代表波场在频率域中，ρ和c_p表示介质中被动源位置处的密度和声波速度。

公式1中，观测波场由在x_s处因垂直点应力产生的垂直质点速度格林函数与震源函数的傅里叶变换的乘积表示。公式1右边对应于在时间域中的互相关。通过对在x_s处的被动源进行积分表示对在地下介质中分布的被动源的叠加。积分结果由期望得到的反射响应(表示在虚震源x_B处激发，在x_A处接收的垂直质点速度波场)与震源函数的能量谱的乘积表示。

若想通过公式1得到正确的反射响应估计，则要求具有相同谱的被动源均匀分布，从所有可能的角度对检波器进行照明。但在实际情况中，被动源的分布相对比较稀疏或成簇分布。在这种情况下，本申请实施例一将公式1离散后近似：

其中，代表在x_s处的单个被动源的互相关函数：

由公式2和3可进一步得知，对于在x_s处的单个被动源，两个检波器x_A和x_B之间的格林函数为：

但是，公式4中包含了震源函数的能量谱该项随震源函数的变化而变换，是一个复杂的频率函数。在实际数据处理中，很难对震源函数能量谱进行估计，因此会在一定程度上影响干涉成像的效果。然而检波器x_A和x_B处接收波场和的反褶积为：

其中为的自相关。由公式5可知，被抵消掉了，因此在进行反褶积干涉的时候，就不需要再对震源函数进行估计，从而提高干涉成像的分辨率和信噪比。

第二步，对反褶积函数进行倾斜叠加处理，并提取零时刻处的峰值对应的射线参数，得到有效照明信息方向。

在被动源分布不均匀的情况下，被动源干涉成像的困难之处在于缺失了用来重构正确反射响应的部分能量。由于不能对干涉假象进行有效压制，导致了对反褶积函数的偏移成像不理想。从被动源辐射出，经自由表面反射后形成的散射波场携带有大量有关地下界面的信息。对此，本申请实施例一将已知虚震源位置的反褶积函数的两处直达波射线路径均提取出来，以对后续的偏移处理进行约束。

图2展示了如何重构两个检波器之间的反射响应。其中，三角形表示检波器，五角星表示被动源。可以看出在稳相点区域的震源辐射出的直达波到达一个检波器后，经过镜像反射到达另一个检波器，而射线路径中均包含虚震源所在的检波器x_s。在横向均匀介质中，对于每一个被动源-虚震源对，均存在一个特定的射线参数来确定该镜像射线。这个射线参数定义了来自虚震源处镜像反射射线的方向。为了找到稳相点区域的射线路径，只需要提取出这个射线参数，而不需要知道被动源x_s的位置。

在常规的线性倾斜叠加中，地震数据可以由算子形式方程d＝Lm表示，其中d是原始时间域地震数据，m是radon模型，L是由观测系统和模型参数定义的radon变换算子。传统最小二乘倾斜叠加方法通过对代价函数求最小极值得到：

其中λ是控制方程解粗糙程度的阻尼因子。radon模型的最小二乘解为m＝(L^TL+λI)^-1L^Td，其中，I为单位矩阵。

由于地震数据是在有限的空间内进行采集，只能在有限的时空域内进行求和，且对所有频率都采用相同的阻尼因子，导致radon模型的能量团无法较好的聚焦，因此该方法求解得到的结果分辨率较低。为了提高radon模型的分辨率，在本申请实施例一中，采用在时间域进行迭代的迭代阈值收缩算法：

m_k＝T_α{m_k-1+2tF^-1[(L^TL)^-1L^T(F[d]-LF[m_k-1])]} (7)

其中k为迭代次数，m_k为迭代k次后的时间域radon模型；F和F^-1分别代表正反傅里叶变换；t为控制收敛速度的步长，通常在0～1之间选择；T_a是2D模型收敛算子，a为控制每次迭代中阈值收缩的算子，通常0<α<1。

对反褶积函数的射线参数分布的分析通过对在时间域中的反褶积函数进行高分辨率线性倾斜叠加来实现。在t＝0时：

其中p是射线参数向量，x_H,A代表x_A的水平坐标，是虚震源x_B的反褶积函数的射线参数分布。其中优势射线参数界定了虚震源波场的有效照明信息方向：

第三步，基于反褶积函数和有效照明信息方向，对虚炮集记录进行偏移成像处理，得到偏移成像结果，完成被动源干涉偏移成像。

从反褶积函数中利用特定的射线参数获得正确的成像结果，需要进行基于方向约束的偏移成像处理。在本实施例一中，成像条件由正传波场和反传波场的互相关来表示，分别从虚震源和检波器位置激发。本实施例一使用高斯波束的高频渐进来重构介质中的格林函数，通过将各个方向的高斯波束叠加来近似表示介质中的波场。每一个高斯波束由介质中紧挨波束的任意位置x的射线中心坐标s(x)和n(x)来表示，现有理论中，在3D介质中，点x_B和任意点x之间的格林函数由每一个高斯波束沿不同方向(由方位角θ和极角ф进行描述)的积分表示。表达式为：

其中射线中心坐标s和n定义了波束通过其附近的观测位置x。表示高斯波束的初始振幅。高斯波束的行为可以通过其宽度和曲率来控制。这些参数在检波器位置处被确定。

在被动源对检波器进行均匀照明的情况下，在偏移处理中，正传波场应该向所有角度辐射。但是，对于偏移单个被动源的反褶积函数来说，正传波场应该被约束于有效照明信息方向。通过利用之前描述的照明分析结果，对反褶积函数的正传波场进行照明方向约束。利用虚震源位置处的介质速度c_p(x_B)，将反褶积函数的坐标从射线参数转化为角度坐标：其中射线参数的水平慢度坐标对虚震源在x_B处，在x处估计的格林函数进行方向约束近似，是通过对反褶积函数的归一化射线参数分布进行加权实现的：

该公式可以通过将照明方向约束于射线参数分布达到的最大值处被简化为：

因此，由被动源辐射出的直达波方向上的高斯波束所构建。

在t时刻的正传波场(或称为下行波场)在虚震源位置x_B处产生，由公式8的格林函数近似所表示：

其中代表被动源的震源函数。震源函数可以根据被动源的瞬态特征从波场的直达波中被估计出，或可以将地震记录的自相关近似作为震源函数。

对于检波器波场的反传，我们使用公式6中的高斯波束近似构建了反褶积函数的渐进形式，并将克希霍夫积分调整至x_A处的检波器边界：

因此，在t时刻的反传波场(或称为上行波场)通过在x处叠加所有方向上的反褶积函数的高斯波束函数来计算：

上行波场包括反褶积函数提供的震源信号的自相关。

反褶积函数的反传估计与正传波场有所不同，反传波场的构建不需要进行特殊的方向约束。

两个波场和的零时延互相关即为偏移的成像条件

其中是由被动源x_s激发，在x_B处的虚震源照明得到的局部成像结果。将每个位置的虚震源得到的成像结果进行叠加，即可得到最终的成像结果：

得到的结果说明地下局部介质可以通过利用单个被动源x_s中包含的约束射线参数信息进行准确地成像。

实施例二

如图3所示，为本申请实施例二中的基于照明方向约束的被动源干涉偏移成像系统结构示意图，包括反褶积模块、照明分析模块和成像模块。

其中，反褶积模块用于对单点脉冲型被动源的地震记录进行干涉处理，得到反褶积函数，反褶积函数用于得到虚炮集记录；照明分析模块用于对反褶积函数进行倾斜叠加处理，并提取零时刻处的峰值对应的射线参数，得到有效照明信息方向；成像模块用于根据根据反褶积函数和有效照明信息方向，对虚炮集记录进行偏移成像处理，得到偏移成像结果。

进一步的，成像模块由正传波场单元、反传波场单元和偏移成像单元组成，具体的，正传波场单元用于根据反褶积函数和有效照明信息方向，得到正传波场；反传波场单元用于基于反褶积函数构建高斯波束函数，并叠加高斯波束函数，得到反传波场；偏移成像单元用于根据正传波场和反传波场，得到偏移成像结果。

以上所述的实施例仅是对本申请优选方式进行的描述，并非对本申请的范围进行限定，在不脱离本申请设计精神的前提下，本领域普通技术人员对本申请的技术方案做出的各种变形和改进，均应落入本申请权利要求书确定的保护范围内。