CN113534243B - 一种被动源Marchenko成像方法及系统 - Google Patents

Abstract

本发明涉及一种被动源Marchenko成像方法及系统，属于地震勘探成像领域，首先采用互相关法对采集到的被动源地震数据进行地震干涉处理，得到各个检波器位置的虚拟源地震数据；然后进行去噪处理以及直达波切除和子波反褶积的预处理，得到反射格林函数响应；采用Marchenko方法对反射格林函数响应进行消除多次波处理，得到消波后的反射格林函数响应；在非精确的速度模型上进行正演，得到直达波响应；将消波后的反射格林函数响应和直达波响应作为输入，重构出上行格林函数和下行格林函数，计算得到反射系数；排列后得到成像结果。无需输入精确速度模型，不会因多次波造成虚假构造，有效提升了成像效果。

Description

一种被动源Marchenko成像方法及系统

技术领域

本发明涉及地震勘探成像领域，特别是涉及一种被动源Marchenko成像方法及系统。

背景技术

地震勘探是地质资源勘查过程中常用的高精度地球物理勘查技术，其利用地震波的传播特征来获取地下构造的空间分布。在常规的地震勘探数据采集过程中，检波器接收到的地震波主要由重锤、炸药等震源激发产生。近年来，随着地震干涉理论的发展，从杂乱无章的地下背景噪声中提取有效的反射地震信号变成了可能。传统的被动源地震勘探成像方法不需要主动激发地震波，只需要将检波器布设于地表并进行长时间采集，就能获得所需要的地震背景噪声记录。背景噪声记录经过地震干涉处理后可以得到虚拟源反射地震记录，进而可以通过常规的地震勘探数据处理方法对虚拟源数据进行处理和成像，这种方法被称为被动源地震勘探。

传统的被动源地震勘探成像方法首先基于互相关法进行反射波场的重构，然后基于传统的主动源地震数据处理方法对地震数据进行处理得到叠加剖面及偏移成像剖面，进而将数据处理的结果用于地下构造的解释。但是，这种方法对速度模型要求很高，由于速度模型不精确导致在最终的成像结果中含有由于多次波造成的虚假构造，成像精度低、效果差，不利于对地下构造的地质解释。

因此，目前亟需一种全新的被动源地震勘探成像方法，以解决传统成像方法需要精确的速度模型作为输入、成像结果因多次波产生虚假构造而导致的成像精度低、成像效果差的问题。

发明内容

本发明的目的是提供一种被动源Marchenko成像方法及系统，结合被动源地震勘探、Marchenko多次波消除和Marchenko方法成像的优点，无需精确的速度模型作为输入，且最大程度保留被动源信号中的有用信息，成像结果也不受多次波的影响造成虚假构造，有效提升了被动源的成像精度和成像效果。

为实现上述目的，本发明提供了如下方案：

一种被动源Marchenko成像方法，包括：

采用互相关法对采集到的被动源地震数据进行地震干涉处理，得到各个检波器位置的虚拟源地震数据；

对所述虚拟源地震数据进行去噪处理，得到去噪后的虚拟源地震数据；

对所述去噪后的虚拟源地震数据分别进行直达波切除和子波反褶积的预处理，得到反射格林函数响应；所述反射格林函数响应包含一次反射波和多次波；

采用Marchenko方法对所述反射格林函数响应进行消除自由表面多次波和层间多次波处理，得到消波后的反射格林函数响应；所述消波后的反射格林函数响应只包含一次反射波；

在非精确的速度模型上进行正演，得到地下虚拟震源位置到地表检波器之间的直达波响应；

将所述消波后的反射格林函数响应和所述直达波响应作为所述Marchenko方法的输入，重构出虚拟检波器位置到地表检波器间的上行格林函数和下行格林函数，并计算得到所有网格成像点的反射系数；

对所述所有网格成像点的反射系数进行排列，得到被动源地震成像结果。

可选的，所述采用互相关法对采集到的被动源地震数据进行地震干涉处理，得到各个检波器位置的虚拟源地震数据，具体包括：

利用公式(1)计算所述虚拟源地震数据：

其中，

分别表示A、B位置处的观测波场，

表示x_B为震源、x_A为检波器的格林函数，

表示功率谱，

表示实部，c和ρ分别代表传播速度和密度，*代表共轭和褶积。

可选的，所述对所述虚拟源地震数据进行去噪处理，得到去噪后的虚拟源地震数据，具体包括：

采用小波变换及F-K滤波方法压制所述虚拟源地震数据中的低频段面波；

利用多域迭代去噪思想，在所述虚拟源地震数据中的共炮点域与共检波点域进行曲波阈值迭代去噪和二维多级中值滤波去噪，并在去噪后的差剖面中提取有效反射波信息，得到所述去噪后的虚拟源地震数据。

可选的，所述对所述去噪后的虚拟源地震数据分别进行直达波切除和子波反褶积的预处理，得到反射格林函数响应；所述反射格林函数响应包含一次反射波和多次波，具体包括：

采用浅地表调查的方式计算所述去噪后的虚拟源地震数据中震源点到各个检波器的直达波到达时间，然后将所述直达波及之前的数据置0；

对于数值模拟的情形，在被动源的噪声源模拟阶段，将Sinc波作为地震子波进行输入，将所述Sinc波作为震源模拟出来的随机位置的多个噪声源依次进行激发，满足互相关法重构反射波格林函数场所需要的条件，将此时得到的虚拟源反射波场作为所述反射格林函数响应；

对于野外实际数据的情形，通过多维反褶积来消除子波的影响，得到包含一次反射波和多次波的所述反射格林函数响应。

可选的，所述采用Marchenko方法对所述反射格林函数响应进行消除自由表面多次波和层间多次波处理，得到消波后的反射格林函数响应；所述消波后的反射格林函数响应只包含一次反射波，具体包括：

利用Marchenko方程对所述反射格林函数响应进行消除自由表面多次波和层间多次波，利用公式(2)计算所述消波后的反射格林函数响应：

其中，U^-表示消波后的反射格林函数响应，x″₀为检波器位置，x'₀和x₀为两个不同的震源位置，t表示时间，t₂表示时间下限，

表示采集边界，d为微分符号，dt'是t'的微分，v^-、

分别为聚焦函数，r为反射系数，R表示地表接收的反射响应。

可选的，所述在非精确的速度模型上进行正演，得到地下虚拟震源位置到地表检波器之间的直达波响应，具体包括：

在非精确的速度模型上进行正演，通过构建Torque作业管理系统将所有直达波正演计算任务平均分配到所述Torque作业管理系统的各个计算核心上，待所有所述直达波正演计算任务计算完毕后，得到地下虚拟震源位置到地表检波器之间的所述直达波响应。

可选的，所述通过构建Torque作业管理系统将所有直达波正演计算任务平均分配到所述Torque作业管理系统的各个计算核心上，具体包括：

假设所述Torque作业管理系统具有m个计算核心，成像区域中n个地下虚拟震源位置对应着n个所述直达波正演计算任务，且n>m；

所述Torque作业管理系统首先将m个所述直达波正演计算任务对应分配到m个所述计算核心中，当m个所述直达波正演计算任务计算结束后，所述Torque作业管理系统继续将剩余的n-m个所述直达波正演计算任务依次分配到m个所述计算核心中，直到所述直达波正演计算任务分配完毕为止；每一个所述计算核心每次仅处理一个所述直达波正演计算任务。

可选的，所述将所述消波后的反射格林函数响应和所述直达波响应作为所述Marchenko方法的输入，重构出虚拟检波器位置到地表检波器间的上行格林函数和下行格林函数，并计算得到所有网格成像点的反射系数，具体包括：

将所述消波后的反射格林函数响应和所述直达波响应作为所述Marchenko方法的输入，重构出虚拟检波器位置到地表检波器间的上行格林函数和下行格林函数；

通过计算所述上行格林函数和所述下行格林函数的互相关函数和多维反褶积，得到地下介质面向目标成像区域所有网格成像点的反射系数：

R₀＝C(Γ+εI)^-1 (4)

其中，x_i、x'_i和x″₀分别表示三个不同的检波器位置，ω表示角频率，

表示采集边界，G^-(x_i,x″₀,ω)表示上行格林函数，G⁺表示下行格林函数，R₀为反射系数，d为微分符号，C表示互相关函数，Γ表示点扩散函数，I表示单位矩阵，ε为常量。

可选的，所述对所述所有网格成像点的反射系数进行排列，得到被动源地震成像结果，具体包括：

对所述所有网格成像点的反射系数按照从左到右和从上到下的顺序进行排列，得到所述被动源地震成像结果。

一种被动源Marchenko成像系统，包括：

虚拟源地震数据获取模块，用于采用互相关法对采集到的被动源地震数据进行地震干涉处理，得到各个检波器位置的虚拟源地震数据；

虚拟源地震数据去噪模块，用于对所述虚拟源地震数据进行去噪处理，得到去噪后的虚拟源地震数据；

反射格林函数响应获取模块，用于对所述去噪后的虚拟源地震数据分别进行直达波切除和子波反褶积的预处理，得到反射格林函数响应；所述反射格林函数响应包含一次反射波和多次波；

多次波消除模块，用于采用Marchenko方法对所述反射格林函数响应进行消除自由表面多次波和层间多次波处理，得到消波后的反射格林函数响应；所述消波后的反射格林函数响应只包含一次反射波；

直达波响应获取模块，用于在非精确的速度模型上进行正演，得到地下虚拟震源位置到地表检波器之间的直达波响应；

上下行格林函数构建模块，用于将所述消波后的反射格林函数响应和所述直达波响应作为所述Marchenko方法的输入，重构出虚拟检波器位置到地表检波器间的上行格林函数和下行格林函数，并计算得到所有网格成像点的反射系数；

被动源地震成像结果获取模块，用于对所述所有网格成像点的反射系数进行排列，得到被动源地震成像结果。

根据本发明提供的具体实施例，本发明公开了以下技术效果：

本发明提出了一种被动源Marchenko成像方法，采用互相关法对采集到的被动源地震数据进行地震干涉处理，得到各个检波器位置的虚拟源地震数据；对虚拟源地震数据进行去噪处理，得到去噪后的虚拟源地震数据，去除无用的噪声，得到信噪比较高的虚拟源地震数据；对去噪后的虚拟源地震数据分别进行直达波切除和子波反褶积的预处理，得到反射格林函数响应，反射格林函数响应包含一次反射波和多次波，预处理过程有效解决了被动源地震数据基于Marchenko方程的高精度成像问题，实现了基于纯被动源数据对地下构造的精确成像；采用Marchenko方法对反射格林函数响应进行消除自由表面多次波和层间多次波处理，得到消波后的反射格林函数响应，消波后的反射格林函数响应只包含一次反射波，从而有效去除了多次波，仅保留一次反射波，解决了多次波会造成虚假构造的问题；在非精确的速度模型上进行正演，得到地下虚拟震源位置到地表检波器之间的直达波响应；将消波后的反射格林函数响应和直达波响应作为Marchenko方法的输入，重构出虚拟检波器位置到地表检波器间的上行格林函数和下行格林函数，并计算得到所有网格成像点的反射系数；对所有网格成像点的反射系数进行排列，得到被动源地震成像结果。

结合被动源地震勘探、Marchenko多次波消除和Marchenko方法成像的优点，无需精确的速度模型作为输入，成像结果也不会受多次波的影响造成虚假构造，且最大程度保留被动源信号中的有用信息，有效提升了被动源成像精度和成像效果，为地震勘探提供了新的被动源地震数据反演成像思路，可以广泛应用于自然保护区、城市及其周边等不适合主动源放炮的地区，具备生态保护意义和广阔的应用前景。

附图说明

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

图1为本发明实施例1提供的被动源Marchenko成像方法的流程图；

图2为本发明实施例1提供的被动源Marchenko成像方法的原理图；

图3为本发明实施例1提供的速度模型的示意图；

图4为本发明实施例1提供的只包含一次反射波的消波后的反射格林函数响应的示意图；

图5为本发明实施例1提供的基于Marchenko方法重构出的(1000，500)位置处的格林函数的示意图；

图6为本发明实施例1提供的互相关函数法得到的被动源Marchenko成像结果示意图；

图7为本发明实施例1提供的多维反褶积法得到的被动源Marchenko成像结果示意图；

图8为本发明实施例2提供的被动源Marchenko成像系统的结构框图。

具体实施方式

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

本发明的目的是提供一种被动源Marchenko成像方法，在纯被动源数据情形下，将被动源背景噪声数据应用在Marchenko成像中，基于Marchenko方程，对重构得到的反射地震响应进行层间多次波和自由表面多次波的消除，进而得到符合Marchenko成像所需要条件的、仅保留一次反射波的反射格林函数响应即格林函数，然后将得到的反射格林函数响应和直达波响应以及一个非精确的速度模型作为输入，实现纯被动源数据情形下，无需精确速度模型的Marchenko成像，消除了传统成像方法中的虚假构造，提升了成像精度和成像效果。

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

实施例1

首先，对本发明涉及的概念进行解释说明：

地震勘探：基于地震波动理论，对检波器接收到的波动进行处理，可以对地下介质构造进行成像。

被动源：相对于主动源的地震波主要由重锤、炸药等震源激发产生，被动源无需主动激发地震波，故此时检波器接收到的地震波是被动产生的。

Marchenko成像：基于Marchenko方法推导出的、对地下介质构造进行成像的地震勘探成像方法。

多次波：由于传统的地震成像方法都是基于波的单次散射假设，所以当地震波在地层间进行反射的时候，容易产生多次波，这类波在传统的地震成像中容易产生构造假象。多次波包括层间多次波和自由表面多次波。

格林函数：点源的脉冲响应。

地震勘探是地质资源勘查过程中常用的高精度地球物理勘查技术，其利用地震波的传播特征来获取地下构造的空间分布。在常规的地震勘探数据采集过程中，检波器接收到的地震波主要由重锤、炸药等震源激发产生。近年来，随着地震干涉理论的发展，从杂乱无章的地下背景噪声中提取有效的反射地震信号变成了可能。地震干涉理论由Claerbout于1968年提出，Calvert于2004基于地震干涉理论提出“虚震源法”，可以将震源重置到检波点位置，该方法在勘探地震数据上得到了广泛应用。这种全新的地震信号采集方式不需要主动激发地震波，只需要将检波器布设于地表并进行长时间采集，就能获得所需要的地震背景噪声记录。背景噪声记录经过地震干涉处理后可以得到虚拟源反射地震记录，进而可以通过常规的地震勘探数据处理方法对虚拟源数据进行处理和成像，这种方法被称为被动源地震勘探。被动源地震勘探作为一种新的勘探方式，由于不需要人工激发地震波，具有绿色、安全、经济的优点。

传统的被动源地震勘探成像方法首先基于互相关法进行反射波场的重构，然后基于传统的主动源地震数据处理方法对地震数据进行处理得到叠加剖面及偏移成像剖面，进而将数据处理的结果用于地下构造的解释。现有的方法技术中，只能用传统的地震成像方法对重构后的数据进行成像，但是，传统的地震成像方法对速度模型要求很高，由于速度模型不精确导致在最终的成像结果中含有由于多次波造成的虚假构造，不利于对地下构造的地质解释。

基于此背景，本发明提出了一种地震勘探的被动源成像方法，尤其是基于纯被动源地震数据的互相关地震干涉法被动源反射波场重构、反射波场多次波消除和数据驱动的面向复杂地下构造目标的高精度Marchenko成像。结合被动源地震勘探、Marchenko多次波消除和Marchenko方法成像的优点，基于纯被动源数据的Marchenko成像，即被动源Marchenko成像方法，其不需要精确的速度模型作为输入，且成像结果不受多次波的影响且最大程度保留了被动源信号中的有用信息，为地震勘探提供新的被动源地震数据反演成像思路。

现有基于Marchenko方程的成像方法中，对于格林函数的重构需要输入一个地表主动源放炮得到的地震数据。结合本发明提出的对重构出的反射格林函数响应数据的子波反褶积、直达波切除等预处理方法，解决了被动源地震数据基于Marchenko方程的高精度成像问题，达到了无需精确的速度模型，也可以基于纯被动源数据对地下构造进行精确成像。

如图1所示，本实施例示出了一种被动源Marchenko成像方法，包括：

S1、采用互相关法对采集到的被动源地震数据进行地震干涉处理，得到各个检波器位置的虚拟源地震数据；具体包括：

利用公式(1)计算所述虚拟源地震数据：

其中，

分别表示A、B位置处的观测波场，

表示x_B为震源、x_A为检波器的格林函数，

表示功率谱，

表示实部，c和ρ分别代表传播速度和密度，*代表共轭和褶积。

S2、对所述虚拟源地震数据进行去噪处理，得到去噪后的虚拟源地震数据；所述去噪处理的实质是提取有用反射波并压制无用信号。所述虚拟源地震数据实际上是由原始被动源地震数据经过互相关法处理后重构得到的被动源地震数据。重构是指将地下的杂乱无章的背景噪声信号变成有用的虚拟源地震记录，只有重构得到虚拟源地震记录后，才可以进行Marchenko成像。

步骤S2具体包括以下步骤：

S21、采用小波变换及F-K滤波方法压制所述虚拟源地震数据中的低频段面波。

被动源地震数据属于地下噪声源传播到地表检波器的透射波，所以被动源中的反射波信号属于低频信号，频段范围与面波相似。为此，本实施例采用小波变换及F-K滤波方法压制虚拟源地震数据中的低频段面波。作为一种时频分解方法，小波变换是一种有效的时变滤波方法。本实施例借助于F-K滤波方法，将小波变换成功地应用于地震信号的滤波之中。具体方法是：将一维信号变换到二维频率-时间平面上，通过改变不同的尺度因子，将信号的不同频率分量刻画出来，再通过改变位移参数来描述信号的时间特性。本发明将F-K滤波方法与小波变换结合应用，是一种十分有效的滤波去噪方法。

S22、对于虚拟源地震记录中包含的较强观测噪声，利用多域迭代去噪思想，在所述虚拟源地震数据中的共炮点域与共检波点域进行曲波阈值迭代去噪和二维多级中值滤波去噪，并在去噪后的差剖面中提取有效反射波信息，得到所述去噪后的虚拟源地震数据。

S3、对所述去噪后的虚拟源地震数据分别进行直达波切除和子波反褶积的预处理，得到反射格林函数响应；所述反射格林函数响应包含一次反射波和多次波。本实施例中，预处理过程有效解决了被动源地震数据基于Marchenko方程的高精度成像问题，实现了基于纯被动源数据对地下构造的精确成像。具体包括以下步骤：

本实施例中，直达波切除预处理就是选取直达波到到达时间的数据，消除直达波的振幅和相位的过程，直达波切除预处理具体包括：

S31、采用浅地表调查的方式计算所述去噪后的虚拟源地震数据中震源点到各个检波器的直达波到达时间，然后将所述直达波及之前的数据置0；

本实施例中，子波反褶积预处理具体包括：

S32、对于数值模拟的情形，在被动源的噪声源模拟阶段，将Sinc波作为地震子波进行输入，将所述Sinc波作为震源模拟出来的随机位置的多个噪声源依次进行激发，满足互相关法重构反射波格林函数场所需要的条件，将此时得到的虚拟源反射波场作为所述反射格林函数响应；

S33、对于野外实际数据的情形，通过多维反褶积来消除子波的影响，得到包含一次反射波和多次波的所述反射格林函数响应。

需要说明的是，本发明中，被动源Marchenko成像方法的核心是基于互相关地震干涉法对具有大量背景噪声的被动源地震数据进行反射波场重构，对重构出的反射波场进行多次波消除以及面向地下构造目标的Marchenko成像。由于需要将经过子波反褶积之后的反射波场数据作为Marchenko方法多次波消除和成像的输入，因此，本实施例提出了基于数值模拟和实际应用两种实现途径的子波反褶积处理方法：

(1)如步骤S32所示，对于数值模拟的情形，在被动源的噪声源模拟阶段，将Sinc波作为地震子波进行输入，此时的Sinc波在频谱上的形态是平的，所以可以认为将Sinc波作为震源模拟出来的随机位置的多个噪声源按照一定的次序激发，满足互相关法重构反射波格林函数场所需要的条件，此时得到的虚拟源反射波场可以认为是反射格林函数响应。

(2)如步骤S33所示，对于野外实际数据的情形，由于基于互相关法重构出的格林函数受到噪声源震源子波的影响，所以可以通过多维反褶积方法来消除子波的影响，得到包含一次波、层间多次波的反射格林函数响应。

然后，将此时得到的反射格林函数响应作为Marchenko方法多次波消去和成像的输入，来分别实现迭代求取只包含一次反射波的消波后的反射格林函数响应，以及地下目标区域的上行格林函数和下行格林函数，最终基于互相关函数法和多维反褶积法来实现面向目标地下构造的被动源Marchenko成像，如图2所示。

S4、采用Marchenko方法对所述反射格林函数响应进行消除自由表面多次波和层间多次波处理，得到消波后的反射格林函数响应；所述消波后的反射格林函数响应只包含一次反射波，从而有效去除了自由表面多次波和层间多次波等多次波，仅保留了一次反射波，进而解决了多次波会造成虚假构造的问题。

步骤S4具体包括以下步骤：

利用Marchenko方程对所述反射格林函数响应进行消除自由表面多次波和层间多次波，利用公式(2)计算所述消波后的反射格林函数响应：

其中，U^-表示消波后的反射格林函数响应，x″₀为检波器位置，x'₀和x₀为两个不同的震源位置，t表示时间，t₂表示时间下限，

表示采集边界，d为微分符号，dt'是t'的微分，v^-、

分别为聚焦函数，r为反射系数，R表示地表接收的反射响应。

S5、在非精确的速度模型上进行正演，得到地下虚拟震源位置到地表检波器之间的直达波响应。由于这个过程由于需要进行大量的正演计算，所以本实施例通过构建Torque作业管理系统来将正演计算任务平均分配到集群或者并行机器的物理核心上，以大幅度提高运算效率，从而提高成像效率。具体包括：

在非精确的速度模型上进行正演，通过构建Torque作业管理系统将所有直达波正演计算任务平均分配到所述Torque作业管理系统的各个计算核心上，待所有所述直达波正演计算任务计算完毕后，得到地下虚拟震源位置到地表检波器之间的所述直达波响应。

本实施例中，通过构建Torque作业管理系统将所有直达波正演计算任务平均分配到所述Torque作业管理系统的各个计算核心上，具体包括：

假设所述Torque作业管理系统具有m个计算核心，成像区域中n个地下虚拟震源位置对应着n个所述直达波正演计算任务，且n>m；

所述Torque作业管理系统首先将m个所述直达波正演计算任务对应分配到m个所述计算核心中，当m个所述直达波正演计算任务计算结束后，所述Torque作业管理系统继续将剩余的n-m个所述直达波正演计算任务依次分配到m个所述计算核心中，直到所述直达波正演计算任务分配完毕为止；每一个所述计算核心每次仅处理一个所述直达波正演计算任务。

所述Torque作业管理系统将n次所述直达波正演计算任务依次、平均地分配到所有计算核心中进行波动方程正演计算。

S6、将所述消波后的反射格林函数响应和所述直达波响应作为所述Marchenko方法的输入，重构出虚拟检波器位置到地表检波器间的上行格林函数和下行格林函数，并计算得到所有网格成像点的反射系数；具体包括：

S61、将所述消波后的反射格林函数响应和所述直达波响应作为所述Marchenko方法的输入，重构出虚拟检波器位置到地表检波器间的上行格林函数和下行格林函数；

S62、通过计算所述上行格林函数和所述下行格林函数的互相关函数和多维反褶积，得到地下介质面向目标成像区域所有网格成像点的反射系数：

R₀＝C(Γ+εI)^-1 (4)

其中，x_i、x'_i和x″₀分别表示三个不同的检波器位置，ω表示角频率，

表示采集边界，G^-(x_i,x″₀,ω)表示上行格林函数，G⁺表示下行格林函数，R₀为反射系数，d为微分符号，C表示互相关函数，Γ表示点扩散函数，I表示单位矩阵，ε为常量。

本实施例中，在得到上行格林函数和下行格林函数后，可以通过计算上、下行格林函数的互相关函数或者多维反褶积得到成像区域中虚拟检波器位置的反射系数的值。其中，地下介质面向目标成像区域所有网格成像点表征的就是虚拟检波器的位置。如步骤S7所示，然后再将无数个反射系数值，按一定顺序和网格间距排列进行成图，得到的图像就是最终的面向目标的、不包含多次波造成虚假构造的被动源地震成像结果。

S7、对所述所有网格成像点的反射系数进行排列，得到被动源地震成像结果。具体包括：

对所述所有网格成像点的反射系数按照从左到右和从上到下的顺序进行排列，得到所述被动源地震成像结果。

应说明，对反射系数进行排序时，与图像的像素点相同。本实施例优选从左到右和从上到下的顺序进行排序，但排序的方式不是固定和唯一的，可根据实际情况而定，会得到不同的成像结果。

为了更加详细地、形象地说明本发明的具体过程，下面采用具体实例进行举例说明：

采用如图3所示的速度模型进行测试，速度模型的水平方向上的长度范围是-6000m～6000m，垂直方向上的深度范围是0～2000m。目标成像区域为水平方向的-1000m～1000m以及垂直方向上400m～1400m。在速度模型上进行噪声源数值模拟，预先在地表布设了601个压力速度检波器，以进行被动源地震数据的采集，接收时长为50min，相邻的检波器的间距为10m，检波器的位置布设在离地表10m的深度。地下的噪声震源总共1500个，这些震源包围着检波器排列。按照本实施例1中方法各个步骤，首先利用互相关法重构出检波器位置的虚拟源地震数据，然后依次对该虚拟源地震数据进行去噪处理以及直达波切除、子波反褶积的预处理，然后再基于Marchenko方法进行多次波消除处理，得到如图4所示的只包含一次反射波的消波后的反射格林函数响应。然后，通过基于输入的非精确速度模型，利用Marchenko方法重构出虚拟检波器位置到地表检波器间的上行格林函数和下行格林函数，(1000m，500m)位置处的格林函数如图5所示。最后，将得到的上行格林函数和下行格林函数基于互相关函数法和多维反褶积法进行Marchenko成像，得到图6所示的互相关函数法的被动源Marchenko成像结果示意图，以及图7所示的多维反褶积法的被动源Marchenko成像结果示意图。

需要理解的是，上述的速度模型涉及的水平方向和垂直方向的模型范围，目标成像区域范围，以及检波器的数量、接收时长、间距、布设深度等，还有震源数量和格林函数的具体位置，这些具体数值仅仅是举例说明，这些值并不是固定的、唯一的，不应作为对本发明的保护范围的限定，还可以设置为其他值，可根据实际情况自行确定。但按照上述方法采用任何数值都应该在本发明的保护范围之内。

还需要说明的是，基于本发明提出的被动源Marchenko成像方法，由于对重构之后的反射格林函数响应数据不包含自由表面多次波和层间多次波，因此，还可以通过传统的地震成像方法进行成像，例如逆时偏移成像法等，逆时偏移成像法为现有技术，此处不再赘述，但无论最后采用何种成像方法实现成像，都应在本发明保护范围之内。

实施例2

如图8所示，本实施例示出了一种被动源Marchenko成像系统，该系统具体包括：

虚拟源地震数据获取模块M1，用于采用互相关法对采集到的被动源地震数据进行地震干涉处理，得到各个检波器位置的虚拟源地震数据；

虚拟源地震数据去噪模块M2，用于对所述虚拟源地震数据进行去噪处理，得到去噪后的虚拟源地震数据；

反射格林函数响应获取模块M3，用于对所述去噪后的虚拟源地震数据分别进行直达波切除和子波反褶积的预处理，得到反射格林函数响应；所述反射格林函数响应包含一次反射波和多次波；

多次波消除模块M4，用于采用Marchenko方法对所述反射格林函数响应进行消除自由表面多次波和层间多次波处理，得到消波后的反射格林函数响应；所述消波后的反射格林函数响应只包含一次反射波；

直达波响应获取模块M5，用于在非精确的速度模型上进行正演，得到地下虚拟震源位置到地表检波器之间的直达波响应；

上下行格林函数构建模块M6，用于将所述消波后的反射格林函数响应和所述直达波响应作为所述Marchenko方法的输入，重构出虚拟检波器位置到地表检波器间的上行格林函数和下行格林函数，并计算得到所有网格成像点的反射系数；

被动源地震成像结果获取模块M7，用于对所述所有网格成像点的反射系数进行排列，得到被动源地震成像结果。

本说明书中各个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。本说明书中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

Claims (8)

1.一种被动源Marchenko成像方法，其特征在于，包括：

采用互相关法对采集到的被动源地震数据进行地震干涉处理，得到各个检波器位置的虚拟源地震数据；

对所述虚拟源地震数据进行去噪处理，得到去噪后的虚拟源地震数据；

对所述去噪后的虚拟源地震数据分别进行直达波切除和子波反褶积的预处理，得到反射格林函数响应；所述反射格林函数响应包含一次反射波和多次波；具体包括：

采用浅地表调查的方式计算所述去噪后的虚拟源地震数据中震源点到各个检波器的直达波到达时间，然后将所述直达波及之前的数据置0；

对于数值模拟的情形，在被动源的噪声源模拟阶段，将Sinc波作为地震子波进行输入，将所述Sinc波作为震源模拟出来的随机位置的多个噪声源依次进行激发，满足互相关法重构反射波格林函数场所需要的条件，将此时得到的虚拟源反射波场作为所述反射格林函数响应；

对于野外实际数据的情形，通过多维反褶积来消除子波的影响，得到包含一次反射波和多次波的所述反射格林函数响应；

采用Marchenko方法对所述反射格林函数响应进行消除自由表面多次波和层间多次波处理，得到消波后的反射格林函数响应；所述消波后的反射格林函数响应只包含一次反射波；

在非精确的速度模型上进行正演，得到地下虚拟震源位置到地表检波器之间的直达波响应；

将所述消波后的反射格林函数响应和所述直达波响应作为所述Marchenko方法的输入，重构出虚拟检波器位置到地表检波器间的上行格林函数和下行格林函数，并计算得到所有网格成像点的反射系数；具体包括：

将所述消波后的反射格林函数响应和所述直达波响应作为所述Marchenko方法的输入，重构出虚拟检波器位置到地表检波器间的上行格林函数和下行格林函数；

通过计算所述上行格林函数和所述下行格林函数的互相关函数和多维反褶积，得到地下介质面向目标成像区域所有网格成像点的反射系数：

R₀＝C(Γ+εI)^-1

其中，x_i、x′_i和x″₀分别表示不同的检波器位置，ω表示角频率，

表示采集边界，G^-(x_i，x″₀，ω)表示上行格林函数，G⁺表示下行格林函数，R₀为反射系数，d为微分符号，C表示互相关函数，Γ表示点扩散函数，I表示单位矩阵，ε为常量；

对所述所有网格成像点的反射系数进行排列，得到被动源地震成像结果。

2.根据权利要求1所述的被动源Marchenko成像方法，其特征在于，所述采用互相关法对采集到的被动源地震数据进行地震干涉处理，得到各个检波器位置的虚拟源地震数据，具体包括：

利用公式(1)计算所述虚拟源地震数据：

其中，

分别表示A、B位置处的观测波场，

表示x_B为震源、x_A为检波器的格林函数，

表示功率谱，

表示实部，c和ρ分别代表传播速度和密度，*代表共轭和褶积。

3.根据权利要求1所述的被动源Marchenko成像方法，其特征在于，所述对所述虚拟源地震数据进行去噪处理，得到去噪后的虚拟源地震数据，具体包括：

采用小波变换及F-K滤波方法压制所述虚拟源地震数据中的低频段面波；

利用多域迭代去噪思想，在所述虚拟源地震数据中的共炮点域与共检波点域进行曲波阈值迭代去噪和二维多级中值滤波去噪，并在去噪后的差剖面中提取有效反射波信息，得到所述去噪后的虚拟源地震数据。

4.根据权利要求1所述的被动源Marchenko成像方法，其特征在于，所述采用Marchenko方法对所述反射格林函数响应进行消除自由表面多次波和层间多次波处理，得到消波后的反射格林函数响应；所述消波后的反射格林函数响应只包含一次反射波，具体包括：

利用Marchenko方程对所述反射格林函数响应进行消除自由表面多次波和层间多次波，利用公式(2)计算所述消波后的反射格林函数响应：

其中，U^-表示消波后的反射格林函数响应，x″₀为检波器位置，x′₀和x₀为两个不同的震源位置，t表示时间，t₂表示时间下限，

表示采集边界，d为微分符号，dt′是t′的微分，v^-、

分别为聚焦函数，r为反射系数，R表示地表接收的反射响应。

5.根据权利要求1所述的被动源Marchenko成像方法，其特征在于，所述在非精确的速度模型上进行正演，得到地下虚拟震源位置到地表检波器之间的直达波响应，具体包括：

在非精确的速度模型上进行正演，通过构建Torque作业管理系统将所有直达波正演计算任务平均分配到所述Torque作业管理系统的各个计算核心上，待所有所述直达波正演计算任务计算完毕后，得到地下虚拟震源位置到地表检波器之间的所述直达波响应。

6.根据权利要求5所述的被动源Marchenko成像方法，其特征在于，所述通过构建Torque作业管理系统将所有直达波正演计算任务平均分配到所述Torque作业管理系统的各个计算核心上，具体包括：

假设所述Torque作业管理系统具有m个计算核心，成像区域中n个地下虚拟震源位置对应着n个所述直达波正演计算任务，且n>m；

所述Torque作业管理系统首先将m个所述直达波正演计算任务对应分配到m个所述计算核心中，当m个所述直达波正演计算任务计算结束后，所述Torque作业管理系统继续将剩余的n-m个所述直达波正演计算任务依次分配到m个所述计算核心中，直到所述直达波正演计算任务分配完毕为止；每一个所述计算核心每次仅处理一个所述直达波正演计算任务。

7.根据权利要求1所述的被动源Marchenko成像方法，其特征在于，所述对所述所有网格成像点的反射系数进行排列，得到被动源地震成像结果，具体包括：

对所述所有网格成像点的反射系数按照从左到右和从上到下的顺序进行排列，得到所述被动源地震成像结果。

8.一种被动源Marchenko成像系统，其特征在于，包括：

虚拟源地震数据获取模块，用于采用互相关法对采集到的被动源地震数据进行地震干涉处理，得到各个检波器位置的虚拟源地震数据；

虚拟源地震数据去噪模块，用于对所述虚拟源地震数据进行去噪处理，得到去噪后的虚拟源地震数据；

反射格林函数响应获取模块，用于对所述去噪后的虚拟源地震数据分别进行直达波切除和子波反褶积的预处理，得到反射格林函数响应；所述反射格林函数响应包含一次反射波和多次波；具体包括：

采用浅地表调查的方式计算所述去噪后的虚拟源地震数据中震源点到各个检波器的直达波到达时间，然后将所述直达波及之前的数据置0；

对于数值模拟的情形，在被动源的噪声源模拟阶段，将Sinc波作为地震子波进行输入，将所述Sinc波作为震源模拟出来的随机位置的多个噪声源依次进行激发，满足互相关法重构反射波格林函数场所需要的条件，将此时得到的虚拟源反射波场作为所述反射格林函数响应；

对于野外实际数据的情形，通过多维反褶积来消除子波的影响，得到包含一次反射波和多次波的所述反射格林函数响应；

多次波消除模块，用于采用Marchenko方法对所述反射格林函数响应进行消除自由表面多次波和层间多次波处理，得到消波后的反射格林函数响应；所述消波后的反射格林函数响应只包含一次反射波；

直达波响应获取模块，用于在非精确的速度模型上进行正演，得到地下虚拟震源位置到地表检波器之间的直达波响应；

上下行格林函数构建模块，用于将所述消波后的反射格林函数响应和所述直达波响应作为所述Marchenko方法的输入，重构出虚拟检波器位置到地表检波器间的上行格林函数和下行格林函数，并计算得到所有网格成像点的反射系数；具体包括：

将所述消波后的反射格林函数响应和所述直达波响应作为所述Marchenko方法的输入，重构出虚拟检波器位置到地表检波器间的上行格林函数和下行格林函数；

通过计算所述上行格林函数和所述下行格林函数的互相关函数和多维反褶积，得到地下介质面向目标成像区域所有网格成像点的反射系数：

R₀＝C(Γ+εI)^-1

其中，x_i、x′_i和x″₀分别表示不同的检波器位置，ω表示角频率，

表示采集边界，G^-(x_i，x″₀，ω)表示上行格林函数，G⁺表示下行格林函数，R₀为反射系数，d为微分符号，C表示互相关函数，Γ表示点扩散函数，I表示单位矩阵，ε为常量；

被动源地震成像结果获取模块，用于对所述所有网格成像点的反射系数进行排列，得到被动源地震成像结果。

Images