本发明属于微地震监测技术领域，特别是涉及一种微地震定位及层析成像方法。
背景技术：
水力压裂技术是非常规油气资源，如致密砂岩气、页岩气与煤层气等开发过程中重要的增产措施之一，该技术利用流体传压特性，通过液压致裂等开采手段将流体高压注入地层，使岩石破裂以提高地层渗透率，从而达到增产的目的。微地震监测技术通过观测、定位及分析微地震事件来监测生产活动及地下状态；该技术可用于水力压裂分析，在储层开发监测中发挥重要作用。然而微地震监测中，因为缺乏准确的地下速度模型描述，微地震事件的定位常存在误差。常规的微地震反演中采用的速度模型一般根据地震、声波测井、vsp等资料构建，并在压裂前采用射孔资料来校正速度模型。但是，射孔反演的速度对射孔附近的事件较为准确，而对距离射孔较远的事件来说，速度模型是不准确的。实际微地震监测中发现，水力压裂通常不发生射孔附近，而是距离射孔有几百米的距离。再者，射孔数据在压裂前采集，因此基于射孔资料的速度模型只能够表征压裂前的地层速度。实际上，在水力压裂过程中，高压流体会产生裂缝，同时通过渗流作用改变压裂缝附近地层的孔隙压力，新增裂缝和孔隙压力的改变都会引起地层速度的变化。因此，如果采用射孔资料矫正的速度模型进行微地震定位存在较大误差。
中国专利201410564587.9公开了“一种微地震震源位置和速度模型同时反演方法”，该方法为提高定位精度，考虑到速度模型横向变化的影响，采用三维速度模型描述地下介质，导致待反演的参数急剧增加，且速度模型的反演结果强烈依赖于地震射线对模型的覆盖程度。正因为如此，即使联合了普通层析方程、震源双差层析方程和检波器双差层析方程，未反演微地震发震时刻，反演的结果相对于传统双差层析方程有所提高，但是与真实结果仍然存在较大的误差，特别是速度模型和震源深度定位结果。再者，该方法采用两点间的射线追踪方法，所存在的明显不足是，当为适应增多的速度反演参数而采用更多的走时信息时，都需重新计算两点间的理论走时，计算太耗时。
技术实现要素：
本发明的目的是为克服现有技术所存在的不足而提供一种微地震定位及层析成像方法，本发明针对水平层状介质，利用地面和井中观测走时资料，能够同时反演震源参数和速度结构模型，这里震源参数是指震源位置和发震时刻，速度结构模型是指层速度和层分界面位置。
根据本发明提出的一种微地震定位及层析成像方法，其特征在于，包括如下具体步骤：
步骤1，读入初始模型及反演参数：所述初始模型及反演参数包括检波器的坐标值、拾取的微地震事件到各检波器的观测到时、初始震源参数及初始速度结构模型、待反演参数约束条件、循环迭代终止条件；所述初始震源参数是指震源位置和发震时刻，初始的速度结构模型是指层速度值和层分界面位置；
步骤2，采用基于界面元的最短路径算法计算微地震事件的理论走时和射线路径：具体是在层状介质中，射线传播的节点为界面离散点—界面元，在界面包围的层内，射线沿直线传播，因此网格剖分只在界面上进行，界面元就是射线传播过程中波前的面元控制点；
步骤3，计算当地下介质结构和微地震事件参数同时反演时，到时关于待反演参数的一阶偏导数；
步骤4，构建同时反演地下介质结构和微地震事件参数的层析方程组；
步骤5，共轭梯度法求解层析方程组，得到微地震震源参数和速度结构的更新量：具体是对步骤4中所构建的同时反演地下介质结构和微地震事件参数的层析方程组采用循环迭代的共轭梯度法进行求解，得到层速度、层界面位置、震源位置和发震时刻的更新量；
步骤6，更新并约束速度结构模型和震源参数：具体是将当前的反演参数值加上步骤5中与其对应的更新量，得到更新后的速度结构模型和震源参数，并用步骤1中输入的约束条件对更新后的参数进行约束；
步骤7，迭代终止判断：如果满足循环迭代终止条件，则输出当前的速度结构模型和震源参数；如果不满足，则返回步骤2继续；其中所述迭代终止条件为到时的均方根残差小于某一常数，或者参数更新量小于某一常数，或者循环迭代次数。
本发明提出的一种微地震定位及层析成像方法进一步的优选方案是：
本发明步骤3所述到时关于反演参数的一阶偏导数包含三部分：第一项为到时关于速度变化的一阶偏导数，第二项为到时关于层深度变化的一阶偏导数，第三项为到时关于微地震事件参数变化的的一阶偏导数：
(1)式中tij为第i个微地震震源到第j个检波器的到时，δtij是相应的走时扰动量；vk和dk分别是第k层内的速度分布和第k层界面的深度，δvk和δdk则分别是相应速度和层界面深度的扰动量；xik，k＝1、2、3、4分别代表第i个微地震事件在x、y、z方向的坐标值及发震时刻，δxik，k＝1、2、3、4是相应的微地震参数扰动量；
(1)式中第一项所述的到时关于速度变化的一阶偏导数采用公式(2)计算：
(2)式中l为射线在第k层内的射线长度；
(1)式中第二项所述的到时关于层界面深度变化的一阶偏导数为：
(3)式中向量wj,wj+1，和wz分别为入射单位向量、出射单位向量、及z轴的单位向量，vj和vj+1则分别是层界面两侧的速度值；
(1)式中第三项所述的到时关于微地震事件参数变化的一阶偏导数，由以下步骤计算，首先到时(t＝ta-to)关于发震时刻的一阶偏导数为：
到时关于震源位置三方向(x,y,z)变化的一阶偏导数由公式(5)计算：
(5)式中θh为射线与z轴夹角，
本发明步骤4所述构建同时反演地下介质结构和微地震事件参数的层析方程组是指：
将基于到时数据的非线性同时反演问题归结为带约束的阻尼最小二乘最优化问题，其反演问题的二范数目标函数为：
(6)式中：δt＝δt1,δt2......δtm是微地震事件的实际观测到时和理论数据的到时残差，m则为拾取的初至数量；δm＝[λ1δv；λ2δd；λ3δx；]，其中δv和δd分别是地下介质各层的层速度值和层分界面的相对改变量，δx为微地震事件的坐标和发震时刻的相对改变量，(λ1,λ2,λ3)则为各类参数的更新权重系数；
(6)式解的一般形式为：
本发明与现有技术相比其显著优点在于：
第一，本发明为解决微地震震源参数的速度模型估计不准确所带来的定位误差问题，而利用基于界面元的最短路径算法计算射线路径和走时，推导了到时关于层速度、层界面深度及震源参数的一阶偏导数，结合共轭梯度法求解带约束的阻尼最小二乘问题，发展了一种适用于水平层状介质，利用微地震初至到时数据同时进行速度结构反演和微地震定位的方法，计算速度快，速度结构反演和定位精度高，反演稳定性强。
第二，本发明采用水平层状介质描述地下介质，符合微地震监测的需求，有效减少反演速度结构参数的数目，降低速度结构反演结果对射线覆盖的依赖程度，而射线的覆盖程度主要由微地震事件的数量和分布位置来决定；同时进行微地震定位(震源位置和发震时刻)和层析成像(层速度和层界面深度)，提高微地震定位的精度。
第三，本发明提出的利用地面和井中微地震事件走时信息同时反演震源参数(震源位置和发震时刻)和速度结构模型(层速度和层分界面位置)的方法适用于微地震监测工程领域。
附图说明
图1为本发明提出的一种微地震定位及层析成像方法的流程示意图。
图2为本发明提出的界面元最短路径算法的示意图，其中圆点和三角形分别代表震源和检波器，v1、v2、v3和v4分别是第一、二、三及四层的层速度值。
图3为本发明提出的到时关于界面深度变化的偏导数计算的示意图，其中δd为界面深度变化量，a和b分别为界面深度扰动前和扰动后相应的射线路径，向量wj,wj+1，和wz分别为射线入射单位向量、射线出射单位向量、及z轴的单位向量，vj和vj+1则分别是层界面两侧的速度值，a和b分别代表射线路径a和b在入射段和出射段的长度差。
图4为本发明提出的地质模型与观测系统的示意图，其中三角形和圆点分别代表检波器和震源，黑色平面为层界面。
图5为本发明提出的层析成像结果的示意图。
图6为本发明提出的微地震定位结果中震源位置在xy面投影的示意图。
图7为本发明提出的微地震定位结果中震源位置在xz面投影的示意图。
图8为本发明提出的微地震定位结果中震源位置在yz面投影的示意图。
图9为本发明提出的微地震定位结果中发震时刻反演结果的示意图。
具体实施方式
下面结合附图和实施例对本发明的具体实施方式作进一步的详细说明。
本发明属于微地震监测技术领域，当无法准确获知地下介质的速度结构模型时，且水力压裂过程也会引起地下介质速度的变化，为避免速度模型估计误差的影响，提高微地震定位精度，选择同时进行微地震定位及层析成像的方法十分必要。微地震的特点是：微地震事件与接收器之间的距离通常在几十至几百米的范围，因此考虑地下介质为水平沉积岩层，利用观测到时与理论到时差同时反演地下介质的层速度、层界面位置、震源位置及发震时刻为本发明的基本条件。
如图1所示，本发明提出的一种微地震定位及层析成像方法，包括如下具体步骤：
步骤1，读入初始模型及反演参数：所述初始模型及反演参数包括检波器的坐标值、拾取的微地震事件到各检波器的观测到时、初始震源参数及初始速度结构模型、待反演参数约束条件、循环迭代终止条件；所述初始震源参数是指震源位置和发震时刻，初始的速度结构模型是指层速度值和层分界面位置；
步骤2，采用基于界面元的最短路径算法计算微地震事件的理论走时和射线路径：具体是在层状介质中，射线传播的节点为界面离散点—界面元，在界面包围的层内，射线沿直线传播，因此网格剖分只在界面上进行，界面元就是射线传播过程中波前的面元控制点；
如图2所示，步骤2的具体步骤进一步公开如下：
步骤2-1，模型参数化：
首先将模型参数化，依据层界面的形态位置划分为相对独立的层状区域，相邻区域由层界面相连接。图2中将模型划分为四层，从上而下依次为第一层、第二层、第三层及第四层，速度值分别为v1、v2、v3和v4。然后将模型边界及层界面离散化，图2中黑色圆点即为离散节点，震源和检波器也作为节点；
步骤2-2，最小走时的计算：
然后从震源位置开始进行波前扫描计算。震源位于第四层，即从第四层开始，进行当前层的波前追踪计算。当该层所有节点都计算完毕，波前停留在该层的界面离散节点上(即界面元)，程序同时保存界面元上的走时和相应的射线路径。然后在该界面元上走时最小的点继续进行新层(与当前计算层相邻的层，图中为第三层)的波前扫描计算。此时第三层作为当前层，当前层计算完毕之后，继续新区的计算，如此重复，直至模型中所有节点都计算完毕；
步骤2-3，输出检波器上的走时和射线路径：
最后依据检波器的位置挑选出对应的节点上保存的走时及相应的射线路径。当有多个检波器，不必重新计算，计算速度快；
步骤3，计算当地下介质结构和微地震事件参数同时反演时，到时关于待反演参数的一阶偏导数，具体计算公式及说明如下：
到时关于反演参数的一阶偏导数包含三部分：第一项为到时关于速度变化的一阶偏导数，第二项为到时关于层深度变化的一阶偏导数，第三项为到时关于微地震事件参数变化的的一阶偏导数：
(1)式中tij为第i个微地震震源到第j个检波器的到时，δtij是相应的走时扰动量；vk和dk分别是第k层内的速度分布和第k层界面的深度，δvk和δdk则分别是相应速度和层界面深度的扰动量；xik，k＝1、2、3、4分别代表第i个微地震事件在x、y、z方向的坐标值及发震时刻，δxik，k＝1、2、3、4是相应的微地震参数扰动量；
(1)式中第一项所述的到时关于速度变化的一阶偏导数采用公式(2)计算：
(2)式中l为射线在第k层内的射线长度；
(1)式中第二项所述的到时关于层界面深度变化的一阶偏导数为，计算示意图如图3：
(3)式中向量wj,wj+1，和wz分别为入射单位向量、出射单位向量、及z轴的单位向量，vj和vj+1则分别是层界面两侧的速度值；
(1)式中第三项所述的到时关于微地震事件参数变化的一阶偏导数，由以下步骤计算，首先到时(t＝ta-to)关于发震时刻的一阶偏导数为：
到时关于震源位置三方向(x,y,z)变化的一阶偏导数由公式(5)计算：
(5)式中θh为射线与z轴夹角，
步骤4，构建同时反演地下介质结构和微地震事件参数的层析方程组，具体构建方法如下：
将基于到时数据的非线性同时反演问题归结为带约束的阻尼最小二乘最优化问题，其反演问题的二范数目标函数为：
(6)式中：δt＝δt1,δt2......δtm是微地震事件的实际观测到时和理论数据的到时残差，m则为拾取的初至数量；δm＝[λ1δv；λ2δd；λ3δx；]，其中δv和δd分别是地下介质各层的层速度值和层分界面的相对改变量，δx为微地震事件的坐标和发震时刻的相对改变量，(λ1,λ2,λ3)则为各类参数的更新权重系数；
(6)式解的一般形式为：
步骤5，共轭梯度法求解层析方程组，得到微地震震源参数和速度结构的更新量：具体是对步骤4中所构建的同时反演地下介质结构和微地震事件参数的层析方程组采用循环迭代的共轭梯度法进行求解，得到层速度、层界面位置、震源位置和发震时刻的更新量；
步骤6，更新并约束速度结构模型和震源参数：具体是将当前的反演参数值加上步骤5中与其对应的更新量，得到更新后的速度结构模型和震源参数，并用步骤1中输入的约束条件对更新后的参数进行约束；
步骤7，迭代终止判断：如果满足循环迭代终止条件，则输出当前的速度结构模型和震源参数；如果不满足，则返回步骤2继续；其中所述迭代终止条件为到时的均方根残差小于某一常数，或者参数更新量小于某一常数，或者循环迭代次数。
实施例1为本发明提出的一种微地震定位及层析成像方法的算例，其具体计算结果进一步公开如下：
如图4所示，为本发明实施例1建立的模型，为四层的水平层状介质，速度从地表向下依次为2.0km/s、2.4km/s、2.8km/s和3.2km/s，层界面的深度分别为16m、30m和40m。采用微地震中常用的地表“米”子型和井中联合观测方式(检波器排列如图4中三角形所示)。图4中圆点代表设定的四个微地震事件，其空间坐标分别为s1(20,40,42)、s2(100,100,45)、s3(150,180,48)、s4(170,30,26)，激发时间分别为10ms、15ms、20ms和5ms。
图5是根据实施例1建立的模型同时反演时而得到的速度结构结果。从图5中可以看出，无论是层速度值还是层界面位置，均能很好的与真实模型吻合，除了第二、三层的层速度值存在较小的误差(0.06km/s)。
图6、图7、图8和图9分别是根据实施例1建立的模型同时反演时而得到的震源定位的结果，其中图6、图7、图8分别是反演的震源位置在xy、xz和yz面的投影示意图，图9为发震时刻反演结果示意图。对比发现，即使初始的震源位置和发震时刻偏离真实值较远，反演的震源位置参见图6、图7、图8和发震时刻参见图9都几乎与理论参数完全一致。
以上具体实施方式及实施例是对本发明提出的一种微地震定位及层析成像方法技术思想的具体支持，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在本技术方案基础上所做的任何等同变化或等效的改动，均仍属于本发明技术方案保护的范围。