CN110579795A - 基于被动源地震波形及其逆时成像的联合速度反演方法 - Google Patents

Abstract

本发明涉及一种基于被动源地震波形及其逆时成像的联合速度反演方法，属于被动源地震探测的速度建模领域。主要包括如下步骤：1)对地震数据预处理并提取震源子波；2)计算逆时延拓波场，进行逆时成像和震源定位；3)计算理论地震记录与观测地震记录之间的差；4)建立反演速度的目标函数；5)基于伴随状态法计算目标函数对模型参数的梯度以及更新方向和步长；6)更新纵、横波速度模型；7)判断是否满足迭代终止条件，如果满足则输出定位结果及速度模型，否则返回步骤2。本发明的联合速度建模方法在震源位置未知情况下，可以获得高分辨率速度模型。

Description

基于被动源地震波形及其逆时成像的联合速度反演方法

技术领域

本发明涉及一种基于被动源地震波形及其逆时成像的联合速度反演方法，属于被动源地震探测的速度建模领域。

背景技术

在传统地震勘探中，由天然地震或地下微动诱发的地震波常被视为噪声干扰，在资料处理过程中会通过各种手段对其进行压制。但是由被动源诱发的地震波中同样包含大量的有用信息，可用于探测地下地质结构。被动源地震勘探相比主动源地震勘探具有诸如无需人工震源、施工方便、能连续观测等优势，在油气储层、大坝、地下油库的地质监测中具有较强的应用价值。

速度模型是诸多地震成像的基础，精确的速度模型可直接用于地下地质解释、储层预测等。被动源地震的能量往往很弱，地面监测方式易受噪声干扰，导致被动源地震资料的信噪比较低，地震旅行时往往难以拾取，这时，基于旅行时的定位和反演方法将不能应用。目前被广为研究的基于波形的反演方法，利用观测数据和模拟数据间的误差建立目标函数，其反演的速度模型具有高分辨率的优势。但在被动源地震勘探的反演问题中，由于速度模型、震源位置和发震时刻均为未知量，导致这种基于波形匹配的反演问题具有很强的非线性和不适定性，极容易发生周期跳跃，陷入局部极小解。

在被动源地震速度建模中，常规反演方法均需要已知震源位置。对于低信噪比资料，逆时成像是确定震源位置的有效方法。当速度模型准确时，纵波自相关、横波自相关和纵横波互相关三种成像条件的逆时成像结果将同时在震源位置处聚焦；而当速度模型不准确时，三者的逆时成像结果的聚焦程度和聚焦位置存在差异。因此利用这三种逆时成像对速度误差响应的不同，便可建立目标函数，迭代反演速度模型。但这种基于三种成像结果间差异的成像域速度反演方法获得的速度模型分辨率不高。

因此，本发明将成像域速度反演方法与分辨率较高的波形反演结合起来，形成一种联合的震源定位和速度建模的迭代反演方法，为被动源勘探速度建模提供一种高分辨率、高适用性的方法。

发明内容

本发明提出了一种基于被动源地震波形及其逆时成像的联合速度反演方法，先用被动地震逆时成像进行震源定位，再将不同成像条件逆时成像之间差异和地震记录波形匹配差异联合起来建立目标函数，利用梯度类最优化方法反演速度模型，从而研发了一种被动源地震的高分辨率速度建模方法，并同时给出准确震源位置。

本发明是采用以下的技术方案实现的：基于被动源地震波形及其逆时成像的联合速度反演方法，包括如下步骤：

1)获取观测数据。在地面布设三分量地震仪监测被动源地震数据；

2)对观测数据进行预处理。包括去噪、有效震源事件识别和提取，纵波与横波分离、取包络等处理；

3)提取地震子波。采用高阶统计方法，从观测资料中提取地震子波s；

4)计算目标函数中的权重系数。在反演迭代的初期，主要使用基于不同成像条件的逆时成像间差异的成像域反演方法，建立模型的背景速度场，并获得震源位置。当目标函数下降到一定程度时，逐渐将波形反演加入到反演过程中。在反演迭代后期，主要使用基于波形的反演方法提高反演速度模型的分辨率。因此，本发明采用下式计算权重系数：

其中，k是当前的迭代次数，N是总的迭代次数，J表示当前迭代目标函数值，J_max表示目标函数最大值，即初始目标函数值；

5)计算纵波和横波的逆时延拓波场。将分离的包络纵波和横波地震记录通过声波方程有限差分算法分别进行逆时延拓，获得逆时延拓的纵、横波场。声波方程可表示为：

其中，X表示模型空间位置(x,y,z)，T表示地震记录时间长度，t表示波传播时间，X_r表示模型空间中地震仪位置，D表示位于X_r的地震仪的地震记录，同时作为波场逆时延拓时的震源，v(X)是模型空间X点处的波传播速度，R(X,t)表示在模型空间X和时间t的逆时延拓波场；

6)震源逆时成像和定位。在逆时延拓过程中，对纵、横波逆时延拓波场分别应用纵波自相关成像条件、横波自相关成像条件、纵横波互相关成像条件进行震源成像(后文中分别使用 PP、SS和PS表示)。三种成像函数如下所示：

其中，R_p(X,t)和R_s(X,t)分别为地下空间位置X和传播时间t对应的纵波和横波逆时延拓波场，T_max为记录地震波场的最大时间，I_pp、I_ss和I_ps分别表示三种逆时成像条件的成像结果；7)构建惩罚函数。由于本发明中的被动地震成像域反演方法是基于不同成像条件的成像差异构建成像误差，因此需引入惩罚函数来构建各成像结果之间的残差，便于建立目标函数。惩罚函数定义为：

其中，I_i和F_i分别表示零延迟成像结果及其对应的惩罚函数(i＝PP，PS，SS)，α_i是无量纲参数，决定函数的惩罚范围，该惩罚函数对成像聚焦越好的地方施以越大的惩罚程度。I_i和 F_j之积表示两成像结果i与j之间的成像残差(i和j不同，可取PP，PS，SS)；

8)计算当前速度模型下的正传波场及其地震记录。利用给定的速度模型、步骤三提取的震源子波和步骤五的震源位置，计算正传波场。正传波场的计算采用下列方程：

其中，X表示模型空间位置，t表示时间，X_s表示模型空间中炮点位置，v(X)是模型空间X处的速度，p(X,t)是计算的正演模拟波场，s是从预处理后观测地震记录中提取的子波。地震记录的计算公式如下：

d_cal(r，t|s)＝R(X|X_r，X_s)p(X，t|X_s) (8)

其中，r表示地震仪号，t表示时间，X_s表示模型空间中炮点位置，X表示模型空间位置，X_r是地震仪在模型空间的位置，p(X,t)表示正演模拟波场，R表示将检波点处的波场值提取出来作为该时刻检波点的地震记录，d_cal(r,t)表示计算的地震记录；

9)计算地震记录残差。将计算的地震记录和观测的地震记录作差，得到地震记录的残差，计算公式如下：

其中，s是被动震源事件号，r是地震仪号，t是时间，d_cal(r,t)是计算的地震记录，d_obs(r,t)是观测的地震记录，d_err(r,t)是地震记录的残差；

10)构建目标函数。构建基于被动源地震勘探的三种成像条件逆时成像差异和波形匹配的联合反演目标函数如下：

其中，I_i和F_i分别表示逆时成像结果及其对应的惩罚函数(i＝PP，PS，SS)，s表示被动震源事件，ε_i表示各目标函数项的权重系数(i＝1～6)。d_cal(r,t)是计算的地震记录，d_obs(r,t)是观测的地震记录，α为步骤四中计算的权重系数。本发明通过最小化该目标函数来得到最优的速度模型；

11)判断是否满足迭代终止条件。如果满足终止条件，则结束反演流程，输出速度模型和震源定位结果。否则，继续进行以下步骤；

12)目标函数梯度计算。本发明利用梯度类最优化方法将模型扰动量δv和目标函数扰动量δJ近似表示为线性问题。对于波动方程反演问题，直接求解梯度计算量相当大，在这里使用伴随状态法近似地求解目标函数对模型参数的梯度，可表示为：

其中，表示逆时延拓波场对时间的二阶导数，a_i表示伴随波场(i＝p,s)，由伴随震源正演模拟获得；

13)伴随波场计算。伴随波场可由伴随震源正演模拟获得，表示为：

L_ia_i＝w_i (12)

L_i表示延拓算子，a_i和w_i分别表示伴随波场和伴随震源(i＝p,s)；

14)伴随震源计算。根据对步骤十中式10的推导，可得出如下伴随震源表达式：

其中，I_i和F_i分别表示逆时成像结果及其对应的惩罚函数，表示推导的中间变量(i＝PP， PS，SS)，ε_i表示目标函数各项的权重系数(i＝1～6)，d_err(r,t)是观测与计算地震记录之间的残差，α为步骤四中计算的权重系数；

15)计算目标函数梯度向量的共轭向量，其计算公式如下：

其中，d_k是第k次迭代目标函数梯度的共轭向量，g_k是第k次迭代目标函数梯度向量，β_k是使相邻两次共轭梯度正交的系数，其计算公式如下：

其中，g_k是梯度向量，d_k-1是共轭向量，符号<,>表示内积运算；

16)迭代更新速度模型。速度模型更新的计算公式如下：

v_k+1(X)＝v_k(X)+λ_kd_k(X) (17)

其中，v_k+1(X)是第k次迭代更新后在模型空间X位置的速度，λ是迭代更新的步长，可以通过一维线性搜索方法得到；

17)将更新后的速度模型作为新的初始速度模型，返回步骤四。

进一步地，对方程2和方程7采用时间2阶、空间12阶交错网格有限差分方法求解，同时采用完全匹配层吸收边界条件来压制模型边界处的反射波。

本发明的有益效果是：提高了被动源地震速度建模方法的适用性，将基于不同成像条件逆时成像结果的差异和基于波形的反演方法进行联合，可以在未知震源位置的情况下，同时获得准确的震源定位结果及速度模型。该发明在反演初期，主要利用基于三种逆时成像之间差异的反演方法，构建背景速度模型，得到较为精确的震源位置。而在反演后期，逐渐加入基于波形的反演方法，获取速度模型的高频成分，从而得到高分辨率速度模型。

附图说明

图1为本发明的流程图。

图2为本发明的理论纵波速度模型。

图3为本发明的理论横波速度模型。

图4为本发明的初始纵波速度模型。

图5为本发明的初始横波速度模型。

图6为基于不同逆时成像结果差异方法反演的纵波速度模型。

图7为基于不同逆时成像结果差异方法反演的横波速度模型。

图8为本发明的联合方法反演的纵波速度模型。

图9为本发明的联合方法反演的横波速度模型。

具体实施方式

为了使本发明的目的、技术方案更加清楚明白，下面结合附图和实例，对本发明作进一步详细说明。

本发明的流程图，如图1所示，包括如下步骤：

1)获取观测数据。在地面布设三分量地震仪监测被动源地震数据；

2)对观测数据进行预处理。包括去噪、有效震源事件识别和提取，纵波与横波分离、取包络等处理；

3)提取地震子波。采用高阶统计方法，从观测资料中提取地震子波s；

4)计算目标函数中的权重系数。在反演迭代的初期，主要使用基于不同成像条件逆时成像的差异的成像域反演方法，建立模型的背景速度场，并获得震源位置。当目标函数下降到一定程度时，逐渐将波形反演加入到反演过程中。在反演迭代后期，主要使用基于波形的反演方法提高反演速度模型的分辨率。因此，本发明采用下式计算权重系数：

其中，k是当前的迭代次数，N是总的迭代次数，J表示当前迭代目标函数值，J_max表示目标函数最大值，即初始目标函数值；

5)计算纵波和横波的逆时延拓波场。将分离的包络纵波和横波地震记录通过声波方程有限差分算法分别进行逆时延拓，获得逆时延拓的纵、横波场。声波方程可表示为：

其中，X表示模型空间位置(x,y,z)，T表示地震记录时间长度，t表示波传播时间，X_r表示模型空间中地震仪位置，D表示位于X_r的地震仪的地震记录，同时作为波场逆时延拓时的震源，v(X)是模型空间X点处的波传播速度，R(X,t)表示在模型空间X和时间t的逆时延拓波场；

6)震源逆时成像和定位。在逆时延拓过程中，对纵、横波逆时延拓波场分别应用纵波自相关成像条件、横波自相关成像条件、纵横波互相关成像条件进行震源成像(后文中分别使用 PP、SS和PS表示)。三种成像函数如下所示：

其中，R_p(X,t)和R_s(X,t)分别为地下空间位置X和传播时间t对应的纵波和横波逆时延拓波场，T_max为记录地震波场的最大时间，I_pp、I_ss和I_ps分别表示三种逆时成像条件的成像结果； 7)构建惩罚函数。由于本发明中的被动地震成像域反演方法是基于不同成像条件的成像差异构建成像误差，因此需引入惩罚函数来构建各成像结果之间的残差，便于建立目标函数。惩罚函数定义为：

其中，I_i和F_i分别表示零延迟成像结果及其对应的惩罚函数(i＝PP，PS，SS)，α_i是无量纲参数，决定函数的惩罚范围，该惩罚函数对成像聚焦越好的地方施以越大的惩罚程度。I_i和 F_j之积表示两成像结果i与j之间的成像残差(i和j不同，可取PP，PS，SS)；

8)计算当前速度模型下的正传波场及其地震记录。利用给定的速度模型、步骤三提取的震源子波和步骤五的震源位置，计算正传波场。正传波场的计算采用下列方程：

其中，X表示模型空间位置，t表示时间，X_s表示模型空间中炮点位置，v(X)是模型空间X处的速度，p(X,t)是计算的正演模拟波场，s是从预处理后观测地震记录中提取的子波。地震记录的计算公式如下：

d_cal(r，t|s)＝R(X|X_r，X_s)p(X，t|X_s) (8)

其中，r表示地震仪号，t表示时间，X_s表示模型空间中炮点位置，X表示模型空间位置，X_r是地震仪在模型空间的位置，p表示正演模拟波场，R表示将检波点处的波场值提取出来作为该时刻检波点的地震记录，d_cal(r,t)表示计算的地震记录；

9)计算地震记录残差。将计算的地震记录和观测的地震记录作差，得到地震记录的残差，计算公式如下：

其中，s是被动震源事件号，r是地震仪号，t是时间，d_cal(r,t)是计算的地震记录，d_obs(r,t)是观测的地震记录，d_err(r,t)是地震记录的残差；

10)构建目标函数。构建基于被动源地震勘探的三种成像条件逆时成像差异和波形匹配差值的联合目标函数如下：

其中，I_i和F_i分别表示逆时成像结果及其对应的惩罚函数(i＝PP，PS，SS)，s表示被动震源事件，ε_i表示各目标函数项的权重系数(i＝1～6)。d_cal(r,t)是计算的地震记录，d_obs(r,t)是观测的地震记录，α为步骤四中计算的权重系数。本发明通过最小化该目标函数来得到最优的速度模型；

11)判断是否满足迭代终止条件。如果满足终止条件，则结束反演流程，输出速度模型和震源定位结果。否则，继续进行以下步骤；

12)伴随震源计算。根据对步骤十中式10的推导，可得出如下伴随震源表达式：

其中，I_i和F_i分别表示逆时成像结果及其对应的惩罚函数，表示推导的中间变量(i＝PP， PS，SS)，ε_i表示各目标函数项的权重系数(i＝1～6)，d_err(r,t)是观测与计算地震记录之间的残差，α为步骤四中计算的权重系数；

13)伴随波场计算。伴随波场可由伴随震源正演模拟获得，表示为：

L_ia_i＝w_i (13)

L_i表示延拓算子，a_i和w_i分别表示伴随波场和伴随震源(i＝p,s)；

14)目标函数梯度计算。本发明利用梯度类最优化方法将模型扰动量δv和目标函数扰动量δJ近似表示为线性问题。对于波动方程反演问题，直接求解梯度计算量相当大，在这里使用伴随状态法近似地求解目标函数对模型参数的梯度，可表示为：

其中，表示逆时延拓波场对时间的二阶导数，a_i表示伴随波场(i＝p,s)，由伴随震源正演模拟获得；

15)计算目标函数梯度向量的共轭向量，其计算公式如下：

其中，d_k是第k次迭代目标函数梯度的共轭，g_k是第k次迭代目标函数梯度向量，β_k是使相邻两次共轭梯度正交的系数，其计算公式如下：

其中，g_k是梯度向量，d_k-1是共轭向量，符号<,>表示内积运算；

16)计算更新步长。本发明采用三点抛物线拟合法求取合适的更新步长，其表达式如下：

其中，λ₁和λ₂分别为两次不同的试探步长，J₁和J₂表示通过试探模型计算得到的目标函数，λ_k是计算得到的本次迭代最优化步长；

17)迭代更新速度模型。速度模型更新的计算公式如下：

v_k+1(X)＝v_k(X)+λ_kd_k(X) (18)

其中，v_k+1(X)是第k次迭代更新后在模型空间X位置处的速度，λ是迭代更新的步长，可以通过线性搜索的方法得到；

18)将更新后的速度模型作为新的初始速度模型，返回步骤四。

其中，公式2和7中的波场通过时间2阶、空间12阶精度的交错网格有限差分方法求解，采用完全匹配层边界条件来压制模型边界处的反射波。

其中，所述步骤五和八中，常使用由测井资料插值并平滑后的速度模型作为反演所用的初始速度模型。

下面结合具体实施方式，对于本发明的理论模型数据测试进行解释和说明。

为了进一步说明本方法的实现思路和实现过程并说明方法的有效性，用层状速度模型进行测试，并与单独使用成像域反演方法的反演结果进行比较。

S1：建立的真实纵波和横波速度模型，如图2和图3所示。该模型宽度为3000m，深度为3000m。由两层介质组成，上层纵波速度为2500m/s，横波速度为1250m/s，纵横波速度比为2.0；下层纵波速度为3000m/s，横波速度为1666m/s，纵横波速度比为1.8。在纵波模型的上层介质中含有一低速异常体，速度为2000m/s。采用正方形网格离散，网格大小为 10×10m²。

S2：观测系统，模拟被动源地震的地面监测方式。11个震源点均匀设置于深度为2500m处，震源间距为100m，如图2和图3中的星型所示。300个检波器均匀布设于地表，检波距为10m。

S3：利用方程7，通过时间2阶、空间12阶精度的弹性波交错网有限差分方法,，采用完全匹配层边界条件，正演数值模拟得到图2和图3所示速度模型及其震源的观测地震记录。

S4：利用式1，计算权重系数。

S5：建立初始速度模型。初始速度模型往往偏离真实模型。将图2所示的真实纵波速度模型上层速度增加5％，下层速度增加10％，把上层的异速体去掉，经平滑后形成初始纵波速度模型，如图4所示；将图3所示的真实横波速度理论模型的上层和下层速度均增高5％，经平滑后作为初始横波速度模型，如图5所示。

S6：用式2计算初始速度模型的逆时延拓波场。并通过式3、4和5计算逆时成像结果，并对震源进行定位，利用式6计算惩罚函数。

S7：由步骤六中得到的震源位置和当前速度模型，利用公式7计算正演波场。

S8：根据式9计算观测地震记录(S3中的模拟记录)和计算地震记录(S7中的模拟记录)的残差。

S9：根据式10和权重系数计算目标函数值。

S10：根据式11和12计算伴随震源函数，利用式13计算伴随波场。

S11：利用式14计算目标函数对速度的梯度向量。

S12：利用式15和16计算梯度向量的共轭向量。

S13：利用式17计算更新步长。

S13：利用式18更新速度模型。判断更新后的速度模型是否满足迭代终止条件，如果满足，则输出结果；如果不满足，则将更新后的速度模型返回S6进行迭代反演。

图6和图7分别为用基于不同逆时成像差异的方法反演的纵、横波速度模型，该结果良好的反演了速度模型的低频信息，但界面的分辨率不高，同时没有反演出纵波模型上层中的异速体。图8和图9分别为本发明的反演结果。与图6和图7对比，可以发现在该发明的反演结果中，界面的分辨率明显提高，纵波异速体被较好地反映出来。需要说明的是，在反演的纵波和横波模型的下层，由于震源到检波点的有效照明范围有限，照明不足的震源底部和两侧速度未被较好的反演出来，但该速度模型误差对震源定位不会产生影响。

以上所述仅为本发明的较佳实施例而己，并不以本发明为限制，凡在本发明的精神和原则之内所作的均等修改、等同替换和改进等，均应包含在本发明的专利涵盖范围内。

Claims (2)

1.一种基于被动源地震波形及其逆时成像的联合速度反演方法，其特征在于，包括如下步骤：

步骤一：在地面观测被动源地震三分量数据；

步骤二：对观测的被动源地震数据进行去噪、纵横波分离等预处理工作；

步骤三：提取震源子波；

步骤四：计算目标函数中的权重系数；

步骤五：计算纵波和横波的逆时延拓波场；

步骤六：通过逆时成像条件进行震源成像和定位；

步骤七：构建各震源成像结果对应的惩罚函数；

步骤八：由逆时成像获得的震源位置，计算当前模型下的正传波场及其对应的地震记录；

步骤九：计算观测地震记录和计算地震记录之间的差值；

步骤十：利用纵波自相关、横波自相关、纵波和横波互相关成像条件的成像结果差异和观测与计算的地震记录差值计算联合的目标函数；

步骤十一：判断是否满足迭代终止条件，如果满足终止条件，则结束反演流程，输出计算的速度模型和震源定位结果，否则，继续进行如下步骤；

步骤十二：根据目标函数计算伴随震源函数；

步骤十三：由伴随震源函数计算伴随波场；

步骤十四：基于伴随状态法计算目标函数对模型速度参数的梯度；

步骤十五：计算目标函数对模型速度参数梯度的共轭方向；

步骤十六：使用线性搜索法计算更新步长；

步骤十七：计算模型修正量并更新速度模型；

步骤十八：将更新后的速度模型作为新的初始速度模型，返回步骤四。

2.根据权利要求1所述的基于被动源地震波形及其逆时成像的联合速度反演方法，其特征在于，所述步骤十中，构建的基于不同成像条件的成像结果间的差异和观测地震记录与计算地震记录的差值建立的联合目标函数如下所示：

其中，I_i和F_i分别表示定位成像和其对应的惩罚函数(i＝PP，PS，SS)，s表示被动地震事件，ε_i表示各目标函数项的权重系数(i＝1～6)。d_cal(r,t)是由当前速度模型计算的地震记录，d_obs(r,t)是观测的地震记录，α为步骤四中计算的权重系数，可表示为：

其中，k是当前的迭代次数，N是总的迭代次数，J表示当前迭代目标函数值，J_max表示目标函数最大值，即初始目标函数值。本发明利用梯度类最优化算法，通过求解该目标函数的极小值来获取速度模型。