CN102841374B - 基于扫描面正演的伪三维快速微地震正演方法 - Google Patents

Abstract

一种基于扫描面正演的伪三维快速微地震正演方法，包括：建立三维横向各向同性的各向异性模型，坐标系为(X,Y,Z)，在各向异性模型的内部设置有微地震震源，在各向异性模型的上表面上设置有多个地震检波器；对微地震震源和每个地震检波器建立二维垂直面，使得微地震震源和所述每个地震检波器位于二维垂直面上；将二维垂直面垂直投影到各向异性模型的水平坐标系(X,Y)上，根据微地震震源和地震检波器在水平坐标系(X,Y)下的水平坐标，得到方位角θ；根据方位角θ和坐标变换公式对水平坐标系(X,Y)进行旋转，得到第一新坐标系(X¹,Y¹)；利用坐标系(X,Y,Z)的Z轴和第一新坐标系(X¹,Y¹)的X¹轴形成第二新坐标系(X¹,Z)；在第二新坐标系(X¹,Z)下，计算从微地震震源发出的微地震波到每个地震检波器点的初至旅行时。

Description

基于扫描面正演的伪三维快速微地震正演方法

技术领域

本发明涉及石油天然气地震勘探领域，具体地讲，涉及一种基于扫描面正演的伪三维快速微地震正演方法，主要应用于石油地震勘探的微地震监测。

背景技术

微地震压裂监测技术是近年来快速发展起来的解决低渗透油气藏开发的重要技术。尽管20世纪70年代已经进行过微地震压裂的监测，但那时只是作为实验性质。随着国家对能源开发的进一步重视及需求，微地震监测与研究对致密砂岩、页岩等非常规油气藏的开发起到至关重要的作用。我国每年至少有上万口井需要做微地震压裂监测。微地震技术是通过对相邻井中的地震检波器接收到的来自压裂井在压裂过程中的微地震信号（目前主要指微地震初至旅行时）进行分析，来描述压裂过程中裂缝生长的几何分布以及流体运移特征（例如裂缝的高度、长度以及方位）。这些信息可以优化压裂设计以及提高油气藏管理，从而提高油气田的产能。

在微地震正演研究中存在两个关键的技术部分：微地震模型建立以及微地震初至旅行时正演计算。由于真实地下油气藏区域（甚至整个地壳内部）的介质属性表现为复杂的各向异性特征，因此原有的微地震各向同性模型不能精确地描述油气藏区域的物理属性。如果不能开发出微地震各向异性模型而采用原有的各向同性模型，则很难有效地利用实际资料中的初至旅行时资料（这些资料记录来自具有各向异性的特征的油气藏）。基于各向同性模型的计算结果可能会导致错误的微地震震源位置以及裂缝空间几何分布。

基于三维各向异性介质模型的微地震初至旅行时正演计算决定了微地震震源定位以及裂缝分布计算的效率。通常，三维微地震波初至旅行时正演计算存在计算效率慢的问题。如何提高微地震波初至旅行时的运算速度且保持相应的精度是微地震研究中的一个重要研究课题。

因此，本发明提出一种基于扫描面正演的各向异性伪三维快速微地震正演的方法。

发明内容

微地震波初至旅行时正演计算是当前进行微地震研究的基础和前提。有效的微地震初至旅行时计算能够保证微地震后期研究（例如微地震的震源定位）的顺利进行。另外，基于三维各向异性介质的微地震波旅行时正演计算存在计算效率慢的问题。因此，提出一种基于扫描面正演的伪三维快速微地震正演的技术，地下地质体被假设为各向异性介质。这项技术具有计算微地震初至旅行时精确、运算速度快以及计算稳定等优点，能够提高对油气藏的开发与利用。本发明涉及微地震模型的建立、坐标系变换以及数值模拟等方法。

根据本发明的一方面，提供一种基于扫描面正演的伪三维快速微地震正演方法，所述方法包括以下步骤：建立三维横向各向同性的各向异性模型，各向异性模型的坐标系为(X,Y,Z)，各向异性模型被划分为多个正方体网格，在各向异性模型的内部设置有微地震震源，在各向异性模型的上表面上布设有多个地震检波器；对微地震震源和每个地震检波器建立二维垂直面，二维垂直面与各向异性模型的上表面垂直，使得微地震震源和所述每个地震检波器位于所述二维垂直面上；将二维垂直面垂直投影到各向异性模型的水平坐标系(X,Y)上，根据微地震震源和地震检波器在水平坐标系(X,Y)下的水平坐标，得到方位角θ；根据方位角θ和坐标变换公式对水平坐标系(X,Y)进行旋转，得到第一新坐标系(X¹,Y¹)，使得微地震震源和地震检波器位于第一新坐标系(X¹,Y¹)的X¹坐标轴上；利用坐标系(X,Y,Z)的Z轴和第一新坐标系(X¹,Y¹)的X¹轴形成第二新坐标系(X¹,Z)；在第二新坐标系(X¹,Z)下，计算从微地震震源发出的微地震波到每个地震检波器点的初至旅行时。

所述各向异性模型可包括关于模型尺寸、网格尺寸、每个网格点的各向异性参数、微地震震源位置以及地震检波器位置的信息。

所述各向异性参数可包括微地震波的垂直相速度以及Thomsen参数ε和δ。

坐标变换公式可以是： $\left(\begin{array}{c}{X}^1\\ Y^1\end{array}\right)=\left(\begin{array}{cc}\cos \left(\mathrm{\θ}\right)& \sin \left(\mathrm{\θ}\right)\\ -\sin \left(\mathrm{\θ}\right)& \cos \left(\mathrm{\θ}\right)\end{array}\right)\·\left(\begin{array}{c}X\\ Y\end{array}\right).$

可采用程函方程有限差分方法来计算从微地震震源发出的微地震波到每个地震检波器点的初至旅行时。

采用程函方程有限差分方法来计算微地震波到每个地震检波器点的初至旅行时的步骤可包括：(1)将二维垂直面划分为多个正方形网格，微地震震源和地震检波器位于二维垂直面的网格点上；(2)将二维垂直面上的微地震震源所在的垂直线作为垂直边界线，将二维垂直面上的微地震震源所在的水平线作为水平边界线；(3)根据微地震震源的发生时间、网格宽度以及微地震波的速度计算微地震波到达水平边界线以及垂直边界线上的各个网格点的旅行时；(4)根据微地震震源的发生时间以及所述多个正方形网格中分别位于垂直边界线以及水平边界线上的与微地震震源相邻的两个网格点的旅行时计算微地震波传播的斜率；(5)根据微地震波传播的斜率，通过反正切函数得到微地震波传播的相角；(6)根据微地震波传播的相角、微地震波的垂直相速度以及Thomsen参数ε和δ，计算微地震波的相速度；(7)根据微地震波的相速度、正方形网格的宽度、微地震震源的发生时间以及正方形网格中与微地震震源相邻的所述两个网格点的旅行时，计算微地震波到二维垂直面上的所述多个正方形网格中的另一网格点的旅行时，所述另一网格点与所述两个网格点相邻并且不在水平边界线和垂直边界线上；(8)按照与步骤(7)相同的方式，在二维垂直面上从微地震震源出发从下到上且从左到右进行迭代计算，对于二维垂直面上的不在水平边界线和垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型顶部的网格点的初至旅行时为止；(9)按照与步骤(7)相同的方式，在二维垂直面上从微地震震源出发从下到上且从右到左进行迭代计算，对于二维垂直面上的不在水平边界线和垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型顶部的网格点的初至旅行时为止；(10)按照与步骤(7)相同的方式，在二维垂直面上从微地震震源出发从上到下且从左到右进行迭代计算，对于二维垂直面上的不在水平边界线和垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型底部的网格点的初至旅行时为止；(11)按照与步骤(7)相同的方式，在二维垂直面上从微地震震源出发从上到下且从右到左进行迭代计算，对于二维垂直面上的不在水平边界线和垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型底部的网格点的初至旅行时为止，由此完成二维垂直面的第一次扫描；(12)按照与步骤(7)相同的方式，从二维垂直面上的垂直边界线与各向异性模型顶部相交的网格点出发从上到下且从左到右进行迭代计算，以先前计算的各向异性模型顶部网格点的初至旅行时作为初始值，对于二维垂直面上的不在垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型底部的网格点的初至旅行时为止；(13)按照与步骤(7)相同的方式，从二维垂直面上的垂直边界线与各向异性模型顶部相交的网格点出发从上到下且从右到左进行迭代计算，以先前计算的各向异性模型顶部网格点的初至旅行时作为初始值，对于二维垂直面上的不在垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型底部的网格点的初至旅行时为止，由此完成二维垂直面的第二次扫描；(14)按照与步骤(7)相同的方式，从二维垂直面上的垂直边界线与各向异性模型底部相交的网格点出发从下到上且从左到右进行迭代计算，以先前计算的各向异性模型底部网格点的初至旅行时作为初始值，对于二维垂直面上的不在垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型顶部的网格点的初至旅行时为止；(15)按照与步骤(7)相同的方式，从二维垂直面上的垂直边界线与各向异性模型底部相交的网格点出发从下到上且从右到左进行迭代计算，以先前计算的各向异性模型底部网格点的初至旅行时作为初始值，对于二维垂直面上的不在垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型顶部的网格点的初至旅行时为止，由此完成二维垂直面的第三次扫描，其中，在二维垂直面的每一次扫描过程中，计算得到的每一个网格点上的旅行时与先前扫描计算所得到的所述网格点上的旅行时进行比较，保留最小旅行时，由此得到微地震波到达地震检波器的初至旅行时。

在上述步骤中，通过二维垂直面的三次扫描，最后取迭代计算的最小旅行时，可进一步提高微地震正演精度。

可将微地震震源的发生时间作为水平边界线上的网格点的初始旅行时，迭代计算微地震波到达水平边界线上的网格点的旅行时，其中，对于位于水平边界线上的特定网格点，根据已经计算出的水平边界线上的与所述特定网格点相邻的网格点的旅行时、网格宽度以及微地震波的速度计算微地震波到达水平边界线上的所述特定网格点的旅行时。

可将微地震震源的发生时间作为垂直边界线上的网格点的初始旅行时，以迭代计算微地震波到达垂直边界线上的网格点的旅行时，其中，对于位于垂直边界线上的特定网格点，根据已经计算出的垂直边界线上的与所述特定网格点相邻的网格点的旅行时、网格宽度以及微地震波的速度计算微地震波到达垂直边界线上的所述特定网格点的旅行时。

附图说明

通过结合附图，从下面的实施例的描述中，本发明这些和/或其它方面及优点将会变得清楚，并且更易于理解，其中：

图1是根据本发明的基于扫描面正演的各向异性伪三维快速微地震正演方法的流程图；

图2示出了根据本发明的三维横向各向同性的各向异性模型的微地震震源与地震检波器关系的示意图；

图3示出了通过微地震震源S(x₀,y₀,z₀)和地震检波器R(x₁,y₁,z₁)的二维垂直面；

图4是连接微地震震源与地震检波器的投影点S₁和R₁的方位角的示意图；

图5是微地震震源和地震检波器在第一新坐标系(X¹，Y¹)下的位置的示意W图；

图6是微地震震源和地震检波器在第二新坐标系(X¹,Z)下的位置的示意图；

图7是旅行时初始化的示意图；

图8示出微地震波初至旅行时的几何示意图；

图9示出三维五层的各向异性模型的示意图；

图10是各向异性模型内的初至旅行时的等值线图；

图11是地表测线得到的初至旅行时的曲线。

具体实施方式

在根据本发明的基于扫描面正演的各向异性伪三维快速微地震正演方法，通过对微地震震源以及每一个地震检波器接收几何排列（位于地表或井中）进行分组，得到多条二维垂直面。对每个二维垂直面，通过垂直投影将垂直面投影到水平面，这时在水平面上为一条直线。进一步进行旋转坐标系变换，然后利用各向异性程函方程有限差分法计算得到微地震初至旅行时。在计算过程中，在模型顶部与模型底部设计了三次扫描计算策略，确保得到的旅行时为最小。

下面参照图1详细描述根据本发明的基于扫描面正演的各向异性伪三维快速微地震正演方法。图1是根据本发明的基于扫描面正演的各向异性伪三维快速微地震正演方法的流程图。

参照图1，在步骤101，建立三维横向各向同性的各向异性模型，坐标系为(X,Y,Z)。各向异性模型被划分为多个正方体网格，在各向异性模型的内部设置有微地震震源，在各向异性模型的上表面(表示地表)上布设有多个地震检波器(每个地震检波器位于各向异性模型的上表面上的相应网格点上)。所述各向异性模型包括模型尺寸、网格尺寸、每个网格点的各向异性参数、微地震震源位置以及地震检波器位置等信息。各向异性参数包括微地震波的垂直相速度v₀以及Thomsen参数（ε和δ）。

图2示出了根据本发明的三维横向各向同性的各向异性模型的微地震震源(微震点)与地震检波器关系的示意图。在图2中，S(x₀,y₀,z₀)表示微地震震源，R(x₁,y₁,z₁)表示地震检波器。微地震震源发出的微地震波被地震检波器接收和检测。

在步骤102，对微地震震源以及每一个地震检波器（位于地表或井中）建立二维垂直面，二维垂直面与各向异性模型的上表面垂直，使得微地震震源以及地震检波器位于该二维垂直面上。图3示出了一个通过微地震震源S(x₀,y₀,z₀)和地震检波器R(x₁,y₁,z₁)的二维垂直面，即，微地震震源S(x₀,y₀,z₀)和地震检波器R(x₁,y₁,z₁)位于该二维垂直面上。

在步骤103，将二维垂直面垂直投影到各向异性模型的水平坐标系(X,Y)上，微地震震源和地震检波器在水平坐标系(X,Y)上的投影点坐标分别为S₁(x₀,y₀)和R₁(x₁,y₁)，然后根据微地震震源和地震检波器的水平坐标，得到方位角θ。

方位角 $\mathrm{\θ}=\arctan \left(\frac{y_1-y_0}{x_1-x_0}\right)$

图4是连接微地震震源与地震检波器的投影点S₁和R₁的方位角的示意图。

在步骤104，根据方位角θ和坐标变换公式对各向异性模型的水平坐标系(X,Y)进行旋转，得到第一新坐标系(X¹,Y¹)，使得微地震震源和地震检波器位于第一新坐标系(X¹,Y¹)的X¹坐标轴上。坐标变换公式如下：

$\left(\begin{array}{c}{X}^1\\ Y^1\end{array}\right)=\left(\begin{array}{cc}\cos \left(\mathrm{\θ}\right)& \sin \left(\mathrm{\θ}\right)\\ -\sin \left(\mathrm{\θ}\right)& \cos \left(\mathrm{\θ}\right)\end{array}\right)\·\left(\begin{array}{c}X\\ Y\end{array}\right)---\left(1\right)$

经过坐标变换，微地震震源和地震检波器在第一新坐标系(X¹,Y¹)下的坐标分为和图5是微地震震源和地震检波器在第一新坐标系(X¹,Y¹)下的位置的示意图。

在步骤105，由于微地震震源和地震检波器位于第一新坐标系(X¹,Y¹)的X¹轴上，所以利用坐标系(X,Y,Z)的Z轴以及第一新坐标系(X¹,Y¹)的X¹轴进一步形成新的垂直坐标系(X¹,Z)，即第二新坐标系(X¹,Z)。在第二新坐标系(X¹,Z)下，微地震震源和地震检波器的坐标变为和图6是为微地震震源和地震检波器在第二新坐标系(X¹,Z)下的位置的示意图。

在步骤106，在第二新坐标系(X¹,Z)下，计算从微地震震源发出的微地震波到每个地震检波器点的初至旅行时。可采用程函方程有限差分方法来计算从微地震震源发出的微地震波到每个地震检波器点的初至旅行时。程函方程有限差分方法相比其它方法计算速度快而且稳定。

下面详细描述采用程函方程有限差分方法来计算从微地震震源发出的微地震波到每个地震检波器点的初至旅行时的步骤。

具体地，采用程函方程有限差分方法来计算初至旅行时的步骤可包括：(1)将二维垂直面划分为多个正方形网格，微地震震源和地震检波器位于二维垂直面的网格点上；(2)将二维垂直面上的微地震震源所在的垂直线作为垂直边界线，将二维垂直面上的微地震震源所在的水平线作为水平边界线；(3)根据微地震震源的发生时间、网格宽度以及微地震波的速度计算微地震波到达水平边界线以及垂直边界线上的各个网格点的旅行时；(4)根据微地震震源的发生时间以及所述多个正方形网格中分别位于垂直边界线以及水平边界线上的与微地震震源相邻的两个网格点的旅行时计算微地震波传播的斜率；(5)根据微地震波传播的斜率，通过反正切函数得到微地震波传播的相角；(6)根据微地震波传播的相角、微地震波的垂直相速度以及Thomsen参数ε和δ，计算微地震波的相速度；(7)根据微地震波的相速度、正方形网格的宽度、微地震震源的发生时间以及正方形网格中与微地震震源相邻的所述两个网格点的旅行时，计算微地震波到二维垂直面上的所述多个正方形网格中的另一网格点的旅行时，所述另一网格点与所述两个网格点相邻并且不在水平边界线和垂直边界线上；(8)按照与步骤(7)相同的方式，在二维垂直面上从微地震震源出发从下到上且从左到右进行迭代计算，对于二维垂直面上的不在水平边界线和垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型顶部的网格点的初至旅行时为止；(9)按照与步骤(7)相同的方式，在二维垂直面上从微地震震源出发从下到上且从右到左进行迭代计算，对于二维垂直面上的不在水平边界线和垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型顶部的网格点的初至旅行时为止；(10)按照与步骤(7)相同的方式，在二维垂直面上从微地震震源出发从上到下且从左到右进行迭代计算，对于二维垂直面上的不在水平边界线和垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型底部的网格点的初至旅行时为止；(11)按照与步骤(7)相同的方式，在二维垂直面上从微地震震源出发从上到下且从右到左进行迭代计算，对于二维垂直面上的不在水平边界线和垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型底部的网格点的初至旅行时为止，由此完成二维垂直面的第一次扫描；(12)按照与步骤(7)相同的方式，从二维垂直面上的垂直边界线与各向异性模型顶部相交的网格点出发从上到下且从左到右进行迭代计算，以先前计算的各向异性模型顶部网格点的初至旅行时作为初始值，对于二维垂直面上的不在垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型底部的网格点的初至旅行时为止；(13)按照与步骤(7)相同的方式，从二维垂直面上的垂直边界线与各向异性模型顶部相交的网格点出发从上到下且从右到左进行迭代计算，以先前计算的各向异性模型顶部网格点的初至旅行时作为初始值，对于二维垂直面上的不在垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型底部的网格点的初至旅行时为止，由此完成二维垂直面的第二次扫描；(14)按照与步骤(7)相同的方式，从二维垂直面上的垂直边界线与各向异性模型底部相交的网格点出发从下到上且从左到右进行迭代计算，以先前计算的各向异性模型底部网格点的初至旅行时作为初始值，对于二维垂直面上的不在垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型顶部的网格点的初至旅行时为止；(15)按照与步骤(7)相同的方式，从二维垂直面上的垂直边界线与各向异性模型底部相交的网格点出发从下到上且从右到左进行迭代计算，以先前计算的各向异性模型底部网格点的初至旅行时作为初始值，对于二维垂直面上的不在垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型顶部的网格点的初至旅行时为止，由此完成二维垂直面的第三次扫描，其中，在二维垂直面的每一次扫描过程中，计算得到的每一个网格点上的旅行时与先前扫描计算所得到的所述网格点上的旅行时进行比较，保留最小旅行时，由此得到微地震波到达地震检波器的初至旅行时。

在上述步骤中，通过二维垂直面的三次扫描，最后取迭代计算的最小旅行时，可进一步提高微地震正演精度。

可将微地震震源的发生时间作为水平边界线上的网格点的初始旅行时，以迭代计算微地震波到达水平边界线上的网格点的旅行时；对于位于水平边界线上的特定网格点(相应地，微地震波相角为0°)，根据已经计算出的水平边界线上的与所述特定网格点相邻的网格点的旅行时、网格宽度以及微地震波的速度计算微地震波到达水平边界线上的所述特定网格点的旅行时。

同样地，可将微地震震源的发生时间作为垂直边界线上的网格点的初始旅行时，以迭代计算微地震波到达垂直边界线上的网格点的旅行时；对于位于垂直边界线上的特定网格点(相应地，微地震波相角为90°)，根据已经计算出的垂直边界线上的与所述特定网格点相邻的网格点的旅行时、网格宽度以及微地震波的速度计算微地震波到达垂直边界线上的所述特定网格点的旅行时。

图7是旅行时初始化的示意图。如图7所示，假设微地震震源位于T₀网格点，T₀为微地震震源的发生时间。二维垂直面的网格宽度h以及微地震波的速度已知，则可以利用等式T₁＝T₀+hs(0°)计算第一初始旅行时T₁，依次类推可以计算得到所有L1线(水平边界线)上网格点的旅行时（实心三角形表示）。s(0°)指水平方向上的慢度。可以利用T₂＝T₀+hs(90°)计算第二初始旅行时T₂，依次类推可以计算得到所有C0线(垂直边界线)上网格点的旅行时（实心圆点表示）。s(90°)指垂直方向上的慢度。以上各个网格点的旅行时只是作为初始值，可在随后的网格点扫描过程中被迭代更新。

下面结合附图和公式具体描述各个网格点的旅行时的迭代计算。前面的T₀以及计算的T₁和T₂可作为迭代计算的旅行时初始值。

假定微地震波的各向异性参数（即，微地震波的垂直相速度v₀以及Thomsen参数ε和δ）在每一个正方形网格单元内为恒定，并且假设微地震波按照相速度（V(α)）以平面波形式向外传播，其中，α为微地震波相角。一旦相角α被确定，则可以计算相速度V(α)，然后进一步计算微地震波到达每个地震检波器的初至旅行时。

在横向各向同性介质的弱-中等各向异性情况下，微地震波的相速度V(α)的计算表达式为：

V(α)=v₀(1+δsin²α+(ε-δ)sin⁴α)    (2)

可以看出，微地震波沿每个传播方向的相速度是不同的。

二维各向同性介质的程函方程如下：

$\frac{1}{V^2}={\left(\frac{\∂\mathrm{\τ}}{\∂x}\right)}^2+{\left(\frac{\∂\mathrm{\τ}}{\∂z}\right)}^2---\left(3\right)$

其中，τ表示时间。

从二维各向同性介质的程函方程可以得到二维各向异性介质的程函方程的近似形式如下：

$\frac{1}{V{\left(\mathrm{\α}\right)}^2}={\left(\frac{\∂\mathrm{\τ}}{\∂x}\right)}^2+{\left(\frac{\∂\mathrm{\τ}}{\∂z}\right)}^2---\left(4\right)$

在程函方程有限差分方法中，根据微地震波到位于二维垂直面的三个网格点的旅行时计算微地震波到二维垂直面的与所述三个网格点相邻的另一网格点的未知旅行时(所述另一网格点不在水平边界线和竖直边界线上)。

如图7所示，t₀、t₁以及t₂为已知的微地震波到二维垂直面的三个网格点的已知旅行时，微地震波到另外一个网格点旅行时t为待求，二维垂直面的正方形网格宽度为h。

通过下式计算微地震波传播的斜率：

然后通过反正切函数得到微地震波传播的相角α。将相角α代入到相速度计算公式（2），从而得到相速度V(α)。进一步将相速度V(α)代入下式：

$t=t_0+\sqrt{2{(h/V\left(\mathrm{\α}\right))}^2-{(t_2-t_1)}^2}---\left(6\right)$

由此，可以得出微地震波到另外一个网格点的旅行时t。

然后按照相同的方式，在二维垂直面上从微地震震源出发从下到上且从左到右进行迭代计算，对于二维垂直面上的不在水平边界线和垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型顶部的网格点的初至旅行时为止；按照相同的方式，在二维垂直面上从微地震震源出发从下到上且从右到左进行迭代计算，对于二维垂直面上的不在水平边界线和垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型顶部的网格点的初至旅行时为止；按照相同的方式，在二维垂直面上从微地震震源出发从上到下且从左到右进行迭代计算，对于二维垂直面上的不在水平边界线和垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型底部的网格点的初至旅行时为止；按照相同的方式，在二维垂直面上从微地震震源出发从上到下且从右到左进行迭代计算，对于二维垂直面上的不在水平边界线和垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型底部的网格点的初至旅行时为止，由此完成二维垂直面的第一次扫描。

然后按照相同的方式，从二维垂直面上的垂直边界线与各向异性模型顶部相交的网格点出发从上到下且从左到右进行迭代计算，以先前计算的各向异性模型顶部网格点的初至旅行时作为初始值，对于二维垂直面上的不在垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型底部的网格点的初至旅行时为止；按照相同的方式，从二维垂直面上的垂直边界线与各向异性模型顶部相交的网格点出发从上到下且从右到左进行迭代计算，以先前计算的各向异性模型顶部网格点的初至旅行时作为初始值，对于二维垂直面上的不在垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型底部的网格点的初至旅行时为止，由此完成二维垂直面的第二次扫描。

然后按照相同的方式，从二维垂直面上的垂直边界线与各向异性模型底部相交的网格点出发从下到上且从左到右进行迭代计算，以先前计算的各向异性模型底部网格点的初至旅行时作为初始值，对于二维垂直面上的不在垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型顶部的网格点的初至旅行时为止；按照相同的方式，从二维垂直面上的垂直边界线与各向异性模型底部相交的网格点出发从下到上且从右到左进行迭代计算，以先前计算的各向异性模型底部网格点的初至旅行时作为初始值，对于二维垂直面上的不在垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型顶部的网格点的初至旅行时为止，由此完成二维垂直面的第三次扫描。

在二维垂直面的每一次扫描过程中，计算得到的每一个网格点上的旅行时与先前扫描计算所得到的所述网格点上的旅行时进行比较，保留最小旅行时，由此得到微地震波到达地震检波器的初至旅行时。

如上所述，在第一次扫描(迭代计算)中，可利用T₀、T₁和T₂分别作为t₀、t₁和t₂的初始值进行计算。另外，对于位于水平边界线或垂直边界线上的特定网格点(相应地，微地震波相角为0°或90°)，可根据已经计算出的水平边界线或垂直边界线上的与所述特定网格点相邻的网格点的旅行时、网格宽度以及微地震波的速度计算微地震波到达水平边界线或垂直边界线上的所述特定网格点的旅行时；对于不在水平边界线和竖直边界线上的网格点，可按照迭代公式(6)计算计算微地震波到达所述网格点的旅行时，直到完成整个二维垂直面的扫描为止。

将本发明应用到三维五层各向异性介质模型中，建立三维横向各向同性的各向异性模型，各层模型各向异性参数如下表所示：

表1模型各向异性参数

层数	v0	ε	δ
				1	2000	0	0
2	2400	0.05	0.05
				3	2800	0.1	0.1
4	3000	0.1	0.1
				5	3500	0.15	0.15

经过坐标变换得到第二新坐标系(X¹,Z)下的模型尺寸为3000×2000，即模型长3000米、深2000米（参见图9所示的三维五层各向异性模型）。网格尺寸为2米，微地震震源位置为网格点（1000，980），如图9中的星号所示。地表每个网格点设置有一个地震检波器，虚线表示测线。经过本发明提出的微地震正演方法，能够快速、准确地计算得到初至旅行时。

图10是各向异性模型内的初至旅行时的等值线图，图11是地表测线得到的初至旅行时的曲线。

模型实验表明，根据本发明的基于扫描面正演的伪三维快速微地震正演技术不仅保持了微地震初至旅行时计算的精确，而且具有运算速度快以及计算稳定等优点，极大地提高了微地震监测的定位精度。

本发明适用于微地震监测软件系统的自主开发。页岩气藏具有较强的各向异性特征，并对定位精度有着重要的影响，本发明能够通过基于各向异性的正演方法，对微地震以后的微地震实时监测及事后精细处理、观测系统设计等具有较强的实际指导作用。

虽然本发明是参照其示例性的实施例被具体描述和显示的，但是本领域的普通技术人员应该理解，在不脱离由权利要求限定的本发明的精神和范围的情况下，可以对其进行形式和细节的各种改变。

Claims (5)

1.一种基于扫描面正演的伪三维快速微地震正演方法，包括以下步骤：

建立三维横向各向同性的各向异性模型，各向异性模型的坐标系为(X,Y,Z)，各向异性模型被划分为多个正方体网格，在各向异性模型的内部设置有微地震震源，在各向异性模型的上表面上设置有多个地震检波器；

对微地震震源和每个地震检波器建立二维垂直面，二维垂直面与各向异性模型的上表面垂直，使得微地震震源和所述每个地震检波器位于所述二维垂直面上；

将二维垂直面垂直投影到各向异性模型的水平坐标系(X,Y)上，根据微地震震源和地震检波器在水平坐标系(X,Y)下的水平坐标，得到方位角θ；

根据方位角θ和坐标变换公式对水平坐标系(X,Y)进行旋转，得到第一新坐标系(X¹,Y¹)，使得微地震震源和地震检波器位于第一新坐标系(X¹,Y¹)的X¹坐标轴上；

利用坐标系(X,Y,Z)的Z轴和第一新坐标系(X¹,Y¹)的X¹轴形成第二新坐标系(X¹,Z)；

在第二新坐标系(X¹,Z)下，计算从微地震震源发出的微地震波到每个地震检波器点的初至旅行时，

其中，坐标变换公式是： $\left(\begin{array}{c}{X}^1\\ Y^1\end{array}\right)=\left(\begin{array}{cc}\cos \left(\mathrm{\θ}\right)& \sin \left(\mathrm{\θ}\right)\\ -\sin \left(\mathrm{\θ}\right)& \cos \left(\mathrm{\θ}\right)\end{array}\right)\·\left(\begin{array}{c}X\\ Y\end{array}\right),$

其中，采用程函方程有限差分方法来计算从微地震震源发出的微地震波到每个地震检波器点的初至旅行时，

其中，采用程函方程有限差分方法来计算微地震波到每个地震检波器点的初至旅行时的步骤包括：

(1)将二维垂直面划分为多个正方形网格，微地震震源和地震检波器位于二维垂直面的网格点上；

(2)将二维垂直面上的微地震震源所在的垂直线作为垂直边界线，将二维垂直面上的微地震震源所在的水平线作为水平边界线；

(3)根据微地震震源的发生时间、网格宽度以及微地震波的速度计算微地震波到达水平边界线以及垂直边界线上的各个网格点的旅行时；

(4)根据微地震震源的发生时间以及所述多个正方形网格中分别位于垂直边界线以及水平边界线上的与微地震震源相邻的两个网格点的旅行时计算微地震波传播的斜率；

(5)根据微地震波传播的斜率，通过反正切函数得到微地震波传播的相角；

(6)根据微地震波传播的相角、微地震波的垂直相速度以及Thomsen参数ε和δ，计算微地震波的相速度；

(7)根据微地震波的相速度、正方形网格的宽度、微地震震源的发生时间以及正方形网格中与微地震震源相邻的所述两个网格点的旅行时，计算微地震波到二维垂直面上的所述多个正方形网格中的另一网格点的旅行时，所述另一网格点与所述两个网格点相邻并且不在水平边界线和垂直边界线上；

(8)按照与步骤(7)相同的方式，在二维垂直面上从微地震震源出发从下到上且从左到右进行迭代计算，对于二维垂直面上的不在水平边界线和垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型顶部的网格点的初至旅行时为止；

(9)按照与步骤(7)相同的方式，在二维垂直面上从微地震震源出发从下到上且从右到左进行迭代计算，对于二维垂直面上的不在水平边界线和垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型顶部的网格点的初至旅行时为止；

(10)按照与步骤(7)相同的方式，在二维垂直面上从微地震震源出发从上到下且从左到右进行迭代计算，对于二维垂直面上的不在水平边界线和垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型底部的网格点的初至旅行时为止；

(11)按照与步骤(7)相同的方式，在二维垂直面上从微地震震源出发从上到下且从右到左进行迭代计算，对于二维垂直面上的不在水平边界线和垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型底部的网格点的初至旅行时为止，由此完成二维垂直面的第一次扫描；

(12)按照与步骤(7)相同的方式，从二维垂直面上的垂直边界线与各向异性模型顶部相交的网格点出发从上到下且从左到右进行迭代计算，以先前计算的各向异性模型顶部网格点的初至旅行时作为初始值，对于二维垂直面上的不在垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型底部的网格点的初至旅行时为止；

(13)按照与步骤(7)相同的方式，从二维垂直面上的垂直边界线与各向异性模型顶部相交的网格点出发从上到下且从右到左进行迭代计算，以先前计算的各向异性模型顶部网格点的初至旅行时作为初始值，对于二维垂直面上的不在垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型底部的网格点的初至旅行时为止，由此完成二维垂直面的第二次扫描；

(14)按照与步骤(7)相同的方式，从二维垂直面上的垂直边界线与各向异性模型底部相交的网格点出发从下到上且从左到右进行迭代计算，以先前计算的各向异性模型底部网格点的初至旅行时作为初始值，对于二维垂直面上的不在垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型顶部的网格点的初至旅行时为止；

(15)按照与步骤(7)相同的方式，从二维垂直面上的垂直边界线与各向异性模型底部相交的网格点出发从下到上且从右到左进行迭代计算，以先前计算的各向异性模型底部网格点的初至旅行时作为初始值，对于二维垂直面上的不在垂直边界线上的其它网格点，根据已经计算出的与所述其它网格点相邻的三个网格点的旅行时计算所述其它网格点的旅行时，直到计算出二维垂直面上的位于各向异性模型顶部的网格点的初至旅行时为止，由此完成二维垂直面的第三次扫描，

其中，在二维垂直面的每一次扫描过程中，计算得到的每一个网格点上的旅行时与先前扫描计算所得到的所述网格点上的旅行时进行比较，保留最小旅行时，由此得到微地震波到达地震检波器的初至旅行时。

2.根据权利要求1所述的伪三维快速微地震正演方法，其中，所述各向异性模型包括关于模型尺寸、网格尺寸、每个网格点的各向异性参数、微地震震源位置以及地震检波器位置的信息。

3.根据权利要求2所述的伪三维快速微地震正演方法，其中，所述各向异性参数包括微地震波的垂直相速度以及Thomsen参数ε和δ。

4.根据权利要求1所述的伪三维快速微地震正演方法，其中，将微地震震源的发生时间作为水平边界线上的网格点的初始旅行时，迭代计算微地震波到达水平边界线上的网格点的旅行时，

其中，对于位于水平边界线上的特定网格点，根据已经计算出的水平边界线上的与所述特定网格点相邻的网格点的旅行时、网格宽度以及微地震波的速度计算微地震波到达水平边界线上的所述特定网格点的旅行时。

5.根据权利要求1所述的伪三维快速微地震正演方法，其中，将微地震震源的发生时间作为垂直边界线上的网格点的初始旅行时，以迭代计算微地震波到达垂直边界线上的网格点的旅行时，

其中，对于位于垂直边界线上的特定网格点，根据已经计算出的垂直边界线上的与所述特定网格点相邻的网格点的旅行时、网格宽度以及微地震波的速度计算微地震波到达垂直边界线上的所述特定网格点的旅行时。