CN105589100A - 一种微地震震源位置和速度模型同时反演方法 - Google Patents

Abstract

本发明提供了一种微地震震源位置和速度模型同时反演方法，属于微震监测中震源定位领域。本方法包括：(1)拾取微地震事件的初至时间；(2)利用声波测井资料，构建初始速度模型，然后获取微地震事件的初始震源位置和发震时刻；(3)基于两点射线追踪计算微地震事件的理论走时，构建联合层析方程组；(4)求解所述联合层析方程组，计算出参数更新量δu、Δx和Δτ；(5)判断是否满足终止条件，如果是，则转入步骤(6)，如果否，则利用所述参数更新量对速度模型、震源位置和发震时刻进行更新，然后返回步骤(3)；(6)输出最终的速度模型、震源位置和发震时刻。

Description

一种微地震震源位置和速度模型同时反演方法

技术领域

本发明属于微震监测中震源定位领域，具体涉及一种微地震震源位置和速度模型同时反演方法，利用微地震事件走时信息同时反演震源位置和速度模型。

背景技术

微地震监测技术是20世纪90年代发展和活跃起来的一项新的非常规技术，它可以应用于油气藏开发、矿山监测、地质灾害监测等诸多领域。微地震监测中震源的位置、发震时刻、震源强度和速度模型都是未知的，确定这些参数是微地震监测的首要任务。

微地震定位时，由于震源位置和速度模型都是未知量，两者耦合在一起。微地震检测时事件的定位精度很大程度取决于速度模型的精度。常规的定位方法都是利用基于测井资料且经过射孔校正后的一维速度模型来进行定位。当地下速度模型存在横向变速时基于一维速度模型的定位精度就会降低，因此就需要对震源位置和速度模型同时进行反演。

Thuber(1983，1992)提出局部地震层析方法，该方法利用局部子空间反演技术将震源参数和速度参数分离，从而实现同时反演三维速度模型和震源位置。国内外地震学家将其广泛应用到地震定位和速度反演中。

Zhang(2003)在双差定位法(Waldhauser和Fllsworth，2000)基础上发展出了双差层析方法，该有效地吸取了双差定位的优势且兼顾反演速度。zhang方法的反演效果要优于局部地震层析方法。

但是，Zhang(2003)的方法在构造层析方程时射线路径在浅层基本重合，不利于浅层速度反演。

发明内容

本发明的目的在于解决上述现有技术中存在的难题，提供一种微地震震源位置和速度模型同时反演方法，引入检波器双差层析方程，联合普通层析方程、震源双差层析方程，三者构建联合层析方程。与zhang的双差层析算法相比，将有利于改善浅层速度反演精度，可以进一步提高微地震反演速度模型的精度以及定位精度。

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

一种微地震震源位置和速度模型同时反演方法，包括：

(1)，拾取微地震事件的初至时间；

(2)，利用声波测井资料，构建初始速度模型，然后获取微地震事件的初始震源位置和发震时刻；

(3)，基于两点射线追踪计算微地震事件的理论走时，构建联合层析方程组；

(4)，求解所述联合层析方程组，计算出参数更新量δu、Δx和Δτ；

(5)，判断是否满足终止条件，如果是，则转入步骤(6)，如果否，则利用所述参数更新量对速度模型、震源位置和发震时刻进行更新，然后返回步骤(3)；

(6)，输出最终的速度模型、震源位置和发震时刻。

所述步骤(3)是这样实现的：

建立联合层析方程组，包括普通层析方程、震源双差层析方程和检波器双差层析方程，其中普通层析方程如公式(2)所示：

$r_k^i=\mathrm{\Δ}{\mathrm{\τ}}^i+\underset{l=1}{\overset{3}{\mathrm{\Σ}}}\frac{\∂T_k^i}{\∂x_l^i}\mathrm{\Δ}{x}_l^i+{\∫}_i^k\mathrm{\δuds}---\left(2\right)$

其中，为微地震事件i到接收点k的计算旅行时和拾取旅行时之间的残差；δu为慢度模型更新量；分别为微地震事件i震源位置x、y、z方向的更新量；Δτⁱ为微地震事件i的发震时刻更新量；震源偏导数为震源处的慢度向量，为：

$\frac{\∂T_k^i}{\∂x_l^i}=-\frac{1}{V}\left(\frac{{\mathrm{dx}}_l^i}{\mathrm{ds}}\right)---\left(3\right)$

所述震源双差层析方程如公式(4)所示：

${\mathrm{dr}}_k^{\mathrm{ij}}=r_k^i-r_k^j=\mathrm{\Δ}{\mathrm{\τ}}^i+\underset{l=1}{\overset{3}{\mathrm{\Σ}}}\frac{\∂T_k^i}{\∂x_l^i}\mathrm{\Δ}{x}_l^i+{\∫}_i^k\mathrm{\δuds}-\mathrm{\Δ}{\mathrm{\τ}}^j-\underset{l=1}{\overset{3}{\mathrm{\Σ}}}\frac{\∂T_k^j}{\∂x_l^j}\mathrm{\Δ}{x}_l^i-{\∫}_j^k\mathrm{\δuds}---\left(4\right)$

其中，分别为微地震事件i、j到接收点k的计算旅行时和拾取旅行时之间的残差；Δτ^j为微地震事件j的发震时刻更新量；为震源双差时间，是记录两个事件的同一检波的旅行时残差之差：

${\mathrm{dr}}_k^{\mathrm{ij}}=r_k^i-r_k^j={(t_k^i-t_k^j)}^{\mathrm{obs}}-{(t_k^i-t_k^j)}^{\mathrm{cal}}---\left(5\right)$

所述检波器双差层析方程如公式(6)所示：

${\mathrm{dt}}_{\mathrm{mn}}^i=r_m^i-r_n^i={\∫}_i^m\mathrm{\δuds}-{\∫}_i^n\mathrm{\δuds}---\left(6\right)$

其中，为检波器双差时间，记录同一事件的两个检波器的旅行时残差之差：

${\mathrm{dr}}_{\mathrm{mn}}^i=r_m^i-r_n^i={(t_m^i-t_n^i)}^{\mathrm{obs}}-{(t_m^i-t_n^i)}^{\mathrm{cal}}---\left(7\right)$

所述步骤(4)中是利用阻尼LSQR方法求解所述联合层析方程组。

所述步骤(5)中利用所述参数更新量对速度模型、震源位置和发震时刻进行更新是这样实现的：

将步骤(2)中的参数分别加上步骤(4)得到的与其对应的更新量δu、Δx和Δτ。

所述步骤(5)中的终止条件如下：

平均旅行时残差，即绝对值算术平均值小于某一常数或参数更新量小于某一常数。即分别给这四个参数设置阀值，这四个常数为人为给定的很小的量，如可以用0.00001。

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

1)，反演稳定性更强。

2)，与传统双差定位方法相比，微地震震源定位精度高；

3)，与双差层析方法相比，反演结果有利于改善浅层速度反演精度。

附图说明

图1a双差层析定位结果中的XY视图

图1b双差层析定位结果中的XZ视图

图1c双差层析定位结果中的YZ视图

图2a本方法定位结果中的XY视图

图2b本方法定位结果中的XZ视图

图2c本方法定位结果中的YZ视图

图3准确速度模型(深度400m)

图4双差层析反演速度模型(深度400m)

图5本方法反演速度模型(深度400m)

图6本方法的步骤框图。

具体实施方式

下面结合附图对本发明作进一步详细描述：

联合层析方程方程组

1)普通层析方程：微地震事件i到接收点k的实际观测旅行时

$t_k^i={\mathrm{\τ}}^i+T_k^i={\mathrm{\τ}}^i+{\∫}_i^k\mathrm{uds}---\left(1\right)$

其中，u为慢度场(速度的倒数)，τⁱ为发震时刻。微地震定位时，震源坐标(x₁，x₂，x₃)、发震时刻和慢度场都是未知的参数。

检波器接收到的到达时(P波或者S波)，利用试验的震源位置，发震时刻和初始速度模型(先验信息)来计算理论到达时两者旅行时的残差为

计算旅行时和观测旅行时残差可以用震源和速度模型的扰动线性近似为

$r_k^i=\mathrm{\Δ}{\mathrm{\τ}}^i+\underset{l=1}{\overset{3}{\mathrm{\Σ}}}\frac{\∂T_k^i}{\∂x_l^i}\mathrm{\Δ}{x}_l^i+{\∫}_i^k\mathrm{\δuds}---\left(2\right)$

其中，震源偏导数为震源处的慢度向量，为

$\frac{\∂T_k^i}{\∂x_l^i}=-\frac{1}{V}\left(\frac{{\mathrm{dx}}_l^i}{\mathrm{ds}}\right)---\left(3\right)$

2)震源双差层析方程：对于相同检波器k，接收到不同的微地震事件i，j，旅行时残差的差，可以得到如下方程(Zhang，2003)

${\mathrm{dr}}_k^{\mathrm{ij}}=r_k^i-r_k^j=\mathrm{\Δ}{\mathrm{\τ}}^i+\underset{l=1}{\overset{3}{\mathrm{\Σ}}}\frac{\∂T_k^i}{\∂x_l^i}\mathrm{\Δ}{x}_l^i+{\∫}_i^k\mathrm{\δuds}-\mathrm{\Δ}{\mathrm{\τ}}^j-\underset{l=1}{\overset{3}{\mathrm{\Σ}}}\frac{\∂T_k^j}{\∂x_l^j}\mathrm{\Δ}{x}_l^i-{\∫}_j^k\mathrm{\δuds}---\left(4\right)$

其中，为震源双差时间(Waldhauser和Fllsworth，2000)，为记录两个事件的同一检波的旅行时残差之差。

${\mathrm{dr}}_k^{\mathrm{ij}}=r_k^i-r_k^j={(t_k^i-t_k^j)}^{\mathrm{obs}}-{(t_k^i-t_k^j)}^{\mathrm{cal}}---\left(5\right)$

3)检波器双差层析方程：对于同一微地震事件i，两个不同检波器，接收点k的另一微地震事件j，也有相似的表达式；两者相减可以得到如下方程

${\mathrm{dt}}_{\mathrm{mn}}^i=r_m^i-r_n^i={\∫}_i^m\mathrm{\δuds}-{\∫}_i^n\mathrm{\δuds}---\left(6\right)$

其中，为检波器双差时间，记录同一事件的两个检波器的旅行时残差之差。

${\mathrm{dr}}_{\mathrm{mn}}^i=r_m^i-r_n^i={(t_m^i-t_n^i)}^{\mathrm{obs}}-{(t_m^i-t_n^i)}^{\mathrm{cal}}---\left(7\right)$

普通层析方程，震源双差层析方程和检波器双差层析方程，即方程(2)、(4)和(6)，联立构成最终的联合层析方程。另外，可以为三种层析方程给定不同的加权系数w₁，w₂，w₃，来控制三者的权重比例。

如图6所示，本发明方法的具体实施步骤如下：

(1)，拾取微地震事件的初至时间；

(2)，利用声波测井资料，构建初始速度模型(水平层状)；通过人为给定或者其它的微地震定位方法(例如Gaiger法、网格搜索法等)，获取微地震事件的初始震源位置和发震时刻；

(3)，基于两点射线追踪计算微地震事件的理论走时；按照公式(2)、(4)和(6)构造联合层析方程(包括普通层析方程、震源双差层析方程和检波器双差层析方程)；

(4)，利用阻尼LSQR方法求解联合层析方程，计算出参数更新量(δu、Δx和Δτ)，进而更新速度模型、震源位置和发震时刻；

(5)，重复步骤(3)(4)，直到满足终止条件，即旅行时残差小于某一常数或参数更新量小于某一常数；

(6)，将最终的速度模型、震源位置和发震时刻等参数作为输出结果输出。

下面利用为一个实施例，用来检验本发明的有效性，具体如下：

图1a-图1c为双差层析定位结果，分别为XY视图、XZ视图和YZ视图，图2a-图2c为本发明的定位结果，分别为XY视图、XZ视图和YZ视图；对比图1a-图1c和图2a-图2c可以看出，本发明反演得到的震源位置更精确。图4为双差层析反演得到的速度模型(深度400m)，图5为本发明反演得到的速度模型(深度400m)，对比真实速度模型(如图3所示)可以看出，本发明在反演速度模型的精细程度上也要优于双差层析方法。

本发明涉及微震监测中震源定位方法，是在初始速度模型不准确时，提高定位精度同时反演速度模型的方法。本发明利用联合层析方程来同时反演出微地震震源位置、发震时刻和速度模型；联合层析方程包括了普通层析方程，震源双差层析方程和检波器双差层析方程。由于同时采用了走时残差和走时残差的差作为优化变量，提高了震源定位精度和速度模型的反演精度，尤其是在速度模型不准确的情况下，也提高了反演结果的稳定性。

上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。

Claims (5)

1.一种微地震震源位置和速度模型同时反演方法，其特征在于：所述方法包括：

(1)，拾取微地震事件的初至时间；

(2)，利用声波测井资料，构建初始速度模型，然后获取微地震事件的初始震源位置和发震时刻；

(3)，基于两点射线追踪计算微地震事件的理论走时，构建联合层析方程组；

(4)，求解所述联合层析方程组，计算出参数更新量δu、Δx和Δτ；

(5)判断是否满足终止条件，如果是，则转入步骤(6)，如果否，则利用所述参数更新量对速度模型、震源位置和发震时刻进行更新，然后返回步骤(3)；

(6)，输出最终的速度模型、震源位置和发震时刻。

2.根据权利要求1所述的微地震震源位置和速度模型同时反演方法，其特征在于：所述步骤(3)是这样实现的：

建立联合层析方程组，包括普通层析方程、震源双差层析方程和检波器双差层析方程，其中普通层析方程如公式(2)所示：

$r_k^i={\mathrm{\Δ\τ}}^i+\underset{l=1}{\overset{3}{\mathrm{\Σ}}}\frac{{\∂}_k^i}{{\∂x}_l^i}{\mathrm{\Δx}}_l^i+{\∫}_i^k\mathrm{\δuds}---\left(2\right)$

其中，为微地震事件i到接收点k的计算旅行时和拾取旅行时之间的残差；δu为慢度模型更新量；分别为微地震事件i震源位置x、y、z方向的更新量；Δτⁱ为微地震事件i的发震时刻更新量；震源偏导数为震源处的慢度向量，为：

$\frac{{\∂T}_k^i}{{\∂x}_l^i}=-\frac{1}{V}\left(\frac{{\mathrm{dx}}_l^i}{\mathrm{ds}}\right)---\left(3\right)$

所述震源双差层析方程如公式(4)所示：

${\mathrm{dr}}_k^{\mathrm{ij}}=r_k^i-r_k^j={\mathrm{\Δ\τ}}^i+\underset{l=1}{\overset{3}{\mathrm{\Σ}}}\frac{{\∂T}_k^i}{{\∂x}_l^i}{\mathrm{\Δx}}_l^i+{\∫}_i^k\mathrm{\δuds}-{\mathrm{\Δ\τ}}^j-\underset{l=1}{\overset{3}{\mathrm{\Σ}}}\frac{{\∂T}_k^j}{{\∂x}_l^j}{\mathrm{\Δx}}_l^j-{\∫}_j^k\mathrm{\δuds}---\left(4\right)$

其中，分别为微地震事件i、j到接收点k的计算旅行时和拾取旅行时之间的残差；Δτ^j为微地震事件j的发震时刻更新量；为震源双差时间，是记录两个事件的同一检波的旅行时残差之差：

${\mathrm{dr}}_k^{\mathrm{ij}}=r_k^i-r_k^j={(t_k^i-t_k^j)}^{\mathrm{obs}}-{(t_k^i-t_k^j)}^{\mathrm{cal}}---\left(5\right)$

所述检波器双差层析方程如公式(6)所示：

${\mathrm{dr}}_{\mathrm{mn}}^i=r_m^i-r_n^i={\∫}_i^m\mathrm{\δuds}-{\∫}_i^n\mathrm{\δuds}---\left(6\right)$

其中，为检波器双差时间，记录同一事件的两个检波器的旅行时残差之差：

${\mathrm{dr}}_{\mathrm{mn}}^i=r_m^i-r_n^i={(t_m^i-t_n^i)}^{\mathrm{obs}}-{(t_m^i-t_n^i)}^{\mathrm{cal}}---\left(7\right).$

3.根据权利要求2所述的微地震震源位置和速度模型同时反演方法，其特征在于：所述步骤(4)中是利用阻尼LSQR方法求解所述联合层析方程组。

4.根据权利要求3所述的微地震震源位置和速度模型同时反演方法，其特征在于：所述步骤(5)中利用所述参数更新量对速度模型、震源位置和发震时刻进行更新是这样实现的：

将步骤(2)中的参数分别加上步骤(4)得到的与其对应的更新量δu、Δx和Δτ。

5.根据权利要求4所述的微地震震源位置和速度模型同时反演方法，其特征在于：所述步骤(5)中的终止条件如下：

平均旅行时残差，即绝对值算术平均值小于某一常数或参数更新量小于某一常数。