CN113885079A - 基于弹性波场解耦的高精度多方位逆时震源成像方法 - Google Patents

Abstract

本发明公开了一种基于弹性波场解耦的高精度多方位逆时震源成像方法,涉及地球物理震源成像技术领域技术领域，该方法以多个方位台网观测到的地震记录作为边界条件，通过求解弹性波方程，计算反传弹性波场，并利用矢量亥姆霍兹分解对反传波场进行纵横波场分离，再根据传统零延迟互相关成像条件，将分离的P波和S波波场的能量在局部时间窗内相乘求和，得到该时刻震源区域破裂能量的分布。本发明使用纵横波信息以及多方位成像条件，不仅能够压制偏移假象，而且能够更加准确的刻画震源区的破裂特征；生成了高精度四维破裂成像结果，能够精细刻画能量在不同深度处的传播特征，对于确定震源位置分布和研究地下隐伏断层至关重要。

Description

基于弹性波场解耦的高精度多方位逆时震源成像方法

技术领域

本发明涉及地球物理震源成像技术领域，尤其涉及一种基于弹性波场解耦的高精度多方位逆时震源成像方法。

背景技术

目前，震源成像是研究震源破裂机制及区域构造的关键，四维震源成像的结果直接反映了震源区域破裂过程的时空演化情况。常规的震源成像的方法包括震源时间函数方法和直接波形方法，前者取决于提取的震源时间函数的精度，后者则需要断层的先验信息来假设断层模型，使得震源成像结果的精度过于依赖初始条件。为了解决这一难题，采用勘探地震学中的逆时偏移成像的原理，震源逆时成像方法应运而生，其使用地震台站接收到的地震记录，对其进行反向延拓并采用互相关成像条件实现震源高精度成像。该方法不仅能够定位震源的位置，还可以可对震源区域的破裂过程进行精细刻画和描述。

震源逆时成像方法主要针对远震数据对震源破裂过程进行成像，成像结果通常仅能反映震源深度平面上的破裂特征，无法精细地刻画深度剖面上的破裂特征，这对于研究构板块间的地震存在非常大的局限性。在实际应用中，震源逆时成像方法需要对单方位角台阵进行成像，导致成像结果存在明显的向台网方向移动的假象，使得成像结果无法正确反映震源的破裂特征。频率域方法是压制假象的一个常用方法，即对来自不同震源和噪音的非相干信号在频域内搜索最相干的能量体。频率域方法能够改善成像质量，但其实现更加复杂，需要更高的计算成本。

因此，对于近震台站接收到的地震数据，急需研发一种适应于近震数据的高精度震源成像方法。

发明内容

为解决现有震源成像方法中存在的成像结果无法正确反映震源的破裂特征的技术问题，本发明公开了一种基于弹性波场解耦的高精度多方位逆时震源成像方法，该方法以多个方位台网观测到的地震记录作为边界条件，通过求解弹性波方程，计算反传弹性波场，并利用矢量亥姆霍兹分解对反传波场进行纵横波场分离，根据传统零延迟互相关成像条件，将分离的P波和S波波场的能量在局部时间窗内相乘求和，得到该时刻震源区域破裂能量的分布。由于近震台站方位角覆盖较大，借助多方位角成像结果相乘叠加的原理，可大幅度提高震源成像的精度。

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

一种基于弹性波场解耦的高精度多方位逆时震源成像方法，具体包括以下步骤：

步骤S1：获取输入数据；

步骤S2：对观测数据进行数据预处理，检查每个台站记录的波形，移除噪音记录、时间中断记录以及同一台站多通道的记录；

步骤S3：对处理后的数据进行频率域滤波，据台站记录的地震波形，选取合适的滤波频段，对所有台站记录波形进行滤波处理；

步骤S4：对台网内不同台站进行方位角划分，根据区域台站方位分布情况，遵循各方位内台站数据尽可能相近的原则，将近震台站划分为多个方位；

步骤S5：进行弹性波场反传，对于划分的每个方位，以台站记录的近震台网地震数据作为边界条件，通过求解弹性波波动方程计算地下的反传波场；

步骤S6：对延拓后的弹性波场进行纵横波分离，利用矢量亥姆霍兹分解，通过求解泊松方程对地下的反传波场进行纵横波分离分解为真振幅矢量P波波场和S波波场；

步骤S7：根据互相关成像条件，利用真振幅矢量P波波场和S波波场对震源区域进行破裂成像；

步骤S8：利用得到的四维震源区域能量分布图，研究该区域震源的破裂的特征，通过震源区域主要能量的分布情况，确定点源的位置分布并高精度的刻画和描述整个地震破裂的动态传播过程。

进一步地，步骤S1中，输入数据包括：研究区域纵波速度模型v_P(x)、横波速度模型v_S(x)、研究区域密度模型ρ(x)、近震台网地震记录d_obs(x_r,t)，x为地下模型空间的位置坐标，x_r代表台站位置。

进一步地，步骤S5中，利用研究区域纵波速度模型v_P(x)、横波速度模型v_S(x)、研究区域密度模型ρ(x)进行弹性波场反传；

弹性波波动方程如下：

其中，t_max为地震记录持续时间；σ(x,t)为应力场；v(x,t)为偏振速度波场；

为一阶时间导数，d_obs(x_r,t_max-t)为观测数据的逆时表达形式；δ(x)为狄克尔delta函数。

L为空间偏导数矩阵，表示为：

且

分别为x、y、z方向的一阶空间偏导；

C为模型参数矩阵，表示为：

进一步地，步骤S6中，矢量亥姆霍兹分解式表示为：

v＝v^P+v^S (4)

泊松方程表示为：

真振幅矢量P波波场v^P(x,t)和S波波场v^S(x,t)如下：

其中，w为满足泊松方程的辅助波场，

为空间梯度算子，

为散度算子，

为旋度算子。

进一步地，步骤S7中，能量成像条件如下：

其中，I(x,t)为震源区域能量分布，|v^P(x,t')|²为P波能量，|v^S(x,t')|²为S波能量，dt为能量叠加时间窗口。

本发明的有益效果是，使用多方位的近震地震数据，通过求解弹性波动方程，计算地下反传波场，使用矢量亥姆霍兹分解将反传波场分解为真振幅矢量P波波场和S波波场，然后采用能量互相关成像条件，对震源区域进行四维的破裂成像。

相比于传统声波介质逆时震源成像方法，本方法使用纵横波信息以及多方位成像条件，不仅能够压制偏移假象，而且能够更加准确的刻画震源区的破裂特征。相比于频率域的反投影震源成像方法，本方法不需要在频域内逐步搜索能量相干体，具有较高的计算效率。相较于时间域震源成像方法，本方法最大的优势在于生成高精度四维破裂成像结果，特别是在深度剖面上分辨率较高，能够精细刻画能量在不同深度处的传播特征，对于确定震源位置分布和研究地下隐伏断层至关重要。

附图说明

图1为Ridgecrest地震序列及近震台站分布图；

图2为弹性波逆时震源成像所需的P波速度场；

图3为Mw7.1主震三分量地震数据的处理结果；

图4为Mw7.1主震CI台网FDR台站记录波形预处理后的结果；

图5为Mw6.4前震(左)和Mw7.1主震(右)的三维震源成像结果；

图6为Mw6.4前震在8km和12km深度切片上的震源成像叠加结果；

图7为Mw6.4前震在深度剖面上的震源成像结果；

图8为Mw7.1主震在8km和12km深度切片上的震源成像叠加结果；

图9为Mw7.1主震在深度剖面上的震源成像结果。

具体实施方式

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

常规的反投影震源成像方法主要使用远震数据，而且存在明显随台网变化偏移的假象，使得成像结果无法正确反映震源区域破裂特征。为获得精确的近震震源破裂成像结果，并压制单方位台网成像中的偏移假象，本发明公开了一种基于弹性波场解耦的高精度多方位逆时震源成像方法，用于准确刻画震源破裂的四维时空演化过程，具体包括以下步骤：

一、获取输入数据。

特别的，输入数据包括：研究区域纵波速度模型v_P(x)、横波速度模型v_S(x)、研究区域密度模型ρ(x)、近震台网地震记录d_obs(x_r,t)，x为地下模型空间的位置坐标，x_r代表台站位置。

二、对观测数据d_obs(x_r,t)进行数据预处理，检查每个台站记录的波形，移除噪音记录、时间中断记录以及同一台站多通道的记录。

三、对处理后的数据进行频率域滤波，据台站记录的地震波形，选取合适的滤波频段，对所有台站记录波形进行滤波处理。

四、对台网内不同台站进行方位角划分，根据区域台站方位分布情况，遵循各方位内台站数据尽可能相近的原则，将近震台站划分为多个方位。

五、进行弹性波场反传，对于划分的每个方位，以台站记录的近震台网地震数据作为边界条件，通过求解弹性波波动方程计算地下的反传波场。

特别的，利用研究区域纵波速度模型v_P(x)、横波速度模型v_S(x)、研究区域密度模型ρ(x)进行弹性波场反传，对于划分的每个方位，以台站记录的近震台网地震记录d_obs(x_r,t_max-t)作为边界条件，通过求解弹性波波动方程。

弹性波波动方程如下：

其中，t_max为地震记录持续时间；σ(x,t)为应力场；v(x,t)为偏振速度波场；

为一阶时间导数，d_obs(x_r,t_max-t)为观测数据的逆时表达形式；δ(x)为狄克尔delta函数。

L为空间偏导数矩阵，表示为：

且

分别为x、y、z方向的一阶空间偏导；

C为模型参数矩阵，表示为：

六、对延拓后的弹性波场进行纵横波分离，利用矢量亥姆霍兹分解，通过求解泊松方程对地下的反传波场进行纵横波分离，分解为真振幅矢量P波波场v^P(x,t)和S波波场v^S(x,t)。

特别的，矢量亥姆霍兹分解式表示为：

v＝v^P+v^S (4)

泊松方程表示为：

真振幅矢量P波波场v^P(x,t)和S波波场v^S(x,t)如下：

其中，w为满足泊松方程的辅助波场，

为空间梯度算子，

为散度算子，

为旋度算子。

七、根据互相关成像条件，利用真振幅矢量P波波场v^P(x,t)和S波波场v^S(x,t)对震源区域进行破裂成像。

特别的，能量成像条件如下：

其中，I(x,t)为震源区域能量分布，|v^P(x,t')|²为P波能量，|v^S(x,t')|²为S波能量，dt为能量叠加时间窗口。

八、利用得到的四维震源区域能量分布图，研究该区域震源的破裂的特征，通过震源区域主要能量的分布情况，确定点源的位置分布并高精度的刻画和描述整个地震破裂的动态传播过程。

将本发明公开的方法应用于2019年Ridgecrest地震序列中，取得了较理想的成像结果。

具体地：

2019年，Ridgecrest地震序列及近震台站分布如图1所示，图中每个三角形代表一个台站，可以看出近震台站方位角覆盖广，支持多方位的震源成像。弹性波逆时震源成像所需的研究区域P波速度模型如图2所示，相应的S波速度和密度由经验公式计算获得。Mw7.1主震三分量地震数据的处理结果如图3所示。Mw7.1主震CI台网FDR台站记录波形预处理后的结果如图4所示，可以看出挑选的地震记录的质量都比较好，滤波频段的选取能够反映主震高频能量。Mw6.4前震(左)和Mw7.1主震(右)的三维震源成像结果(灰色等值面)如图5所示。Mw6.4前震及Mw7.1主震空间平面和深度剖面震源成像结果如图7-9所示，可以看出成像结果不仅可以在不提供源信息的基础上精确定位点源位置，而且高精度的刻画出震源区域在空间平面及深度剖面上的破裂特征，震源成像结果所描述的Mw6.4前震和Mw7.1主震在空间平面上的破裂过程与SCEDC余震目录确定的断层几何形状基本一致，即Mw6.4前震先沿西北方向断层破裂后转向西南方向传播，Mw7.1主震沿西北-东南方向主要断层传播。

本发明公开了一种基于弹性波场解耦的高精度多方位逆时震源成像方法，即利用推导的互相关成像条件对分离的矢量弹性波场，在近震台站的多个方位角同时进行震源成像，对多个方位角的成像结果相乘叠加，得到多方位的高精度的震源成像结果，该反演流程既能压制震源假象，又可以对震源区域特别是在深度剖面上的四维破裂过程进行高精度地刻画，该方法计算成本较低，具有较好的应用前景。

当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。

Claims (5)

1.一种基于弹性波场解耦的高精度多方位逆时震源成像方法，其特征在于，具体包括以下步骤：

步骤S1：获取输入数据；

步骤S2：对观测数据进行数据预处理，检查每个台站记录的波形，移除噪音记录、时间中断记录以及同一台站多通道的记录；

步骤S3：对处理后的数据进行频率域滤波，据台站记录的地震波形，选取合适的滤波频段，对所有台站记录波形进行滤波处理；

步骤S4：对台网内不同台站进行方位角划分，根据区域台站方位分布情况，遵循各方位内台站数据尽可能相近的原则，将近震台站划分为多个方位；

步骤S5：进行弹性波场反传，对于划分的每个方位，以台站记录的近震台网地震数据作为边界条件，通过求解弹性波波动方程计算地下的反传波场；

步骤S6：对延拓后的弹性波场进行纵横波分离，利用矢量亥姆霍兹分解，通过求解泊松方程对地下的反传波场进行纵横波分离分解为真振幅矢量P波波场和S波波场；

步骤S7：根据互相关成像条件，利用真振幅矢量P波波场和S波波场对震源区域进行破裂成像；

步骤S8：利用得到的四维震源区域能量分布图，研究该区域震源的破裂的特征，通过震源区域主要能量的分布情况，确定点源的位置分布并刻画和描述整个地震破裂的动态传播过程。

2.如权利要求1所述的一种基于弹性波场解耦的高精度多方位逆时震源成像方法，其特征在于，步骤S1中，输入数据包括：研究区域纵波速度模型v_P(x)、横波速度模型v_S(x)、研究区域密度模型ρ(x)、近震台网地震记录d_obs(x_r,t)，x为地下模型空间的位置坐标，x_r代表台站位置。

3.如权利要求1所述的一种基于弹性波场解耦的高精度多方位逆时震源成像方法，其特征在于，步骤S5中，利用研究区域纵波速度模型v_P(x)、横波速度模型v_S(x)、研究区域密度模型ρ(x)进行弹性波场反传；

弹性波波动方程如下：

其中，t_max为地震记录持续时间；σ(x,t)为应力场；v(x,t)为偏振速度波场；

为一阶时间导数；d_obs(x_r,t_max-t)为观测数据的逆时表达形式；δ(x)为狄克尔delta函数；

L为空间偏导数矩阵，表示为：

且

分别为x、y、z方向的一阶空间偏导；

C为模型参数矩阵，表示为：

4.如权利要求1所述的一种基于弹性波场解耦的高精度多方位逆时震源成像方法，其特征在于，步骤S6中，矢量亥姆霍兹分解式表示为：

v＝v^P+v^S；

泊松方程表示为：

真振幅矢量P波波场v^P(x,t)和S波波场v^S(x,t)如下：

其中，w为满足泊松方程的辅助波场，

为空间梯度算子，

为散度算子，

为旋度算子。

5.如权利要求1所述的一种基于弹性波场解耦的高精度多方位逆时震源成像方法，其特征在于，步骤S7中，能量成像条件如下：

其中，I(x,t)为震源区域能量分布，|v^P(x,t')|²为P波能量，|v^S(x,t')|²为S波能量，dt为能量叠加时间窗口。

Images