CN104076392A - 基于网格搜索和牛顿迭代的微震震源定位联合反演方法 - Google Patents

Abstract

本发明公开一种微地震反演算法。该算法是基于网格搜索和牛顿迭代的联合反演定位方法。本算法大体分为两个部分。第一部分是在微震事件所发生的区域内进行网格的划分，通过网格搜索法得到一个微地震事件发生的大概位置。第二部分是牛顿法。牛顿法反演时的初始迭代值正是第一部分中网格搜索法的结果。本方法解决了牛顿反演方法过分依赖初始迭代值的缺点。减少迭代次数。能够较好的使得目标函数最终收敛在理论真实值的附近。

Description

基于网格搜索和牛顿迭代的微震震源定位联合反演方法

技术领域

本发明是一种基于最优化理论的微地震定位反演方法。 

背景技术

微震是地震等级大于一级并且小于三级的地震，这样小的等级地震人们不易觉察到，只能用仪器监测。但是即使是很微小的地下震动，也会给地下的介质产生相应的激励，这样的激励可能改变地下介质的力学状态。比如，微震直接或者额间接的在很大程度上增加了矿井顶底板的破裂或者是瓦斯突出的概率。这样的微震给人们的生产生活带来极大的安全隐患。 

现如今，随着微地震监测技术的提高，已被广泛运用在工程领域中。在煤矿系统中，人们利用微震监测技术来判断地下介质的破裂情况，从而减少不必要的人员伤亡，和一系列不必要的财产损失，节约成本。在工程物探领域，人们可以利用微地震进行桩基检测，对建筑物的桩基进行质量检验，判断桩基内是否有裂缝。 

目前，反演定位的方法多种多样，在非线性反演中，有利用几何性质，通过导数求得结果的梯度法，牛顿法，共轭梯度法等。也有在全局随机搜索结果的蒙特卡罗法，还有通过自然规律进行的启发式的算法，如模拟退火法，遗传算法等等。 

这些算法都能在一定程度上有效的解决微地震反演的定位问题。但是，这些算法有存在着一些局限性，那就是对于初始值的给定有着极大的依赖，如果初始值的选取，在算法所能承受的范围之内，那么都过一系列的迭代，理论值会最终收敛在真实值的附近，但是如果初始值与真实值的偏差太大，这会给收敛的过程带来一定的影响，有时甚至决定最终的结果是否正确。 

所以，本算法的目的在于当对于牛顿法反演之前，通过网格搜索的方法，先给出地震源的大致位置，这样就会使得初始值与真实值的偏差不会太大，那么既能够得到真实的发震位置，又能进一步的提高微震反演的速度。 

发明内容

本发明是为了解决当前的反演算法对于初始值过度依赖而且加快反演速度的一种算法。 

基于网格搜索的牛顿反演定位方法的基本实现步骤如下： 

网格搜索部分 

步骤一：对于实行微地震观测的区域，在三维空间中，进行网格体的划分，确定每一 个维度的坐标值的值域，以及划分的格数。 

步骤二：建立代价函数t(x)。 

步骤三：将每一个网格上的点的坐标，代入到目标函数中去，得到相对应的目标函数值。 

步骤四：将每个点所对应目标函数的值代入到卡方分布的概率密度函数中，得到每个点对应的卡方分布密度函数值，并以此值作为权重。 

步骤五：得到每一个坐标点所对应的权重之后，将每一个维度的坐标进行加权求和。最终得到发震地点的大致位置x^(k)。 

牛顿法部分 

步骤六：令目标函数f(x)为在x^(k)位置上进行二阶Taylor级数展开。 

步骤七：对于目标函数进行了二级近似后，对新的目标函数求一次导数。 

步骤八：得到迭代公式，求得结果。 

进一步地，在所述步骤一中，对所述网格体的划分包括：先对划定的区域进行网格的划分，区域的大小，网格的大小，从而划分的点的个数，依据需要进行相应的设置。

进一步地，所述步骤二中，通过发震地点、储层速度，即地震波旅行时的函数关系来确定所述代价函数的建立，具体形式为: 

$t(x,y,z)=\frac{1}{{\mathrm{\ρ}}_T}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(t_i-\frac{1}{v}\×\sqrt{{(x_i-x)}^2+{(y_i-y)}^2+{(z_i-z)}^2})}^2$

其中t_i是地震波的旅行时，x_i,y_i,z_i,t_i,是每一个检波器所在的空间位置以及每一个检波器所得到的地震波的到时差，n是检波器的个数，ρ_T是一个常数。v是地下介质速度

所述步骤三中，每一个经过网格划分得到的点的位置都会对应一个(x,y,z)的空间坐标，将每一个点的(x,y,z)的空间坐标代入到所述步骤二中的目标函数，得到目标函数值。 

进一步地，所述步骤四中，将每个点所对应的所述目标函数的值代入到卡方分布的概率密度函数中，其中卡方分布的密度分布函数的具体形式： 

$p\left(t\right)=\frac{t^{\frac{n}{2}-\frac{3}{2}}{e}^{-\frac{3}{2}}}{\mathrm{\Γ}(\frac{n}{2}-\frac{1}{2})\×2^{\frac{n}{2}-\frac{1}{2}}}$

其中，n是检波器的个数，t是每一个点对应的代价函数值，得到每个点对应的卡方分布密度函数值，并以此值作为权重。Γ的具体形式为： 

$\mathrm{\Γ}\left(n\right)={\∫}_0^{+\∞}{e}^{-n}{x}^{n-1}\mathrm{dx},(n>0)$

公式中e为自然数。 

进一步地，所述步骤五中，得到每一个坐标点所对应的权重之后，将每一个维度的坐标进行加权求和，最终得到发震地点的大致位置x^(k)。 

进一步地，所述步骤六中，令所述目标函数f(x)在x^(k)位置上进行二阶Taylor级数展开，具体为： 

$\begin{array}{c}f\left(x\right)=f\left(x^{\left(k\right)}\right)+\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{\∂f\left(x^{\left(k\right)}\right)}{\∂x_i}\mathrm{\Δ}{x_i}^{\left(k\right)}+\frac{1}{2}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\frac{{\∂}^{2}f\left(x^{\left(k\right)}\right)}{\∂x_i\∂x_j}\mathrm{\Δ}{x_i}^{\left(k\right)}\mathrm{\Δ}{x_j}^{\left(k\right)}\\ =f\left(x^{\left(k\right)}\right)+\stackrel{\‾}{g}{\left(x^{\left(k\right)}\right)}^T\mathrm{\Δ}{x}^{\left(k\right)}+\frac{1}{2}\mathrm{\Δ}{x}^{\left(k\right)T}{H}^{\left(k\right)}\mathrm{\Δ}{x}^k\end{array}$

其中 $\stackrel{\‾}{g}{\left(x^{\left(k\right)}\right)}^T=\left[\begin{array}{cccc}\frac{\∂f\left(x^{\left(k\right)}\right)}{\∂x_1}& \frac{\∂f\left(x^{\left(k\right)}\right)}{\∂x_2}& \·\·\·& \frac{\∂f\left(x^{\left(k\right)}\right)}{\∂x_N}\end{array}\right]$ 为函数在x^(k)的梯度， 

Δx^(k)T＝[Δx₁ ^(k) Δx₂ ^(k) … Δx_N ^(k)]是模型参数的改正向量。其中，x是待求未知量组成的向量，N为未知数的个数。 

H^(k)为黑塞矩阵，它的具体数学表达式为： 

$H^{\left(k\right)}=\left(\begin{array}{cccc}\frac{{\∂}^{2}f\left(x^{\left(k\right)}\right)}{{\∂x}_1\∂x_1}& \frac{{\∂}^{2}f\left(x^{\left(k\right)}\right)}{\∂x_1\∂x_2}& \·\·\·& \frac{{\∂}^{2}f\left(x^{\left(k\right)}\right)}{\∂x_1\∂x_N}\\ \frac{{\∂}^{2}f\left(x^{\left(k\right)}\right)}{\∂x_2\∂x_1}& \frac{{\∂}^{2}f\left(x^{\left(k\right)}\right)}{\∂x_2\∂x_2}& \·\·\·& \frac{{\∂}^{2}f\left(x^{\left(k\right)}\right)}{\∂x_2\∂x_N}\\ \·& \·& & \·\\ \·& \·& & \·\\ \·& \·& & \·\\ \frac{{\∂}^{2}f\left(x^{\left(k\right)}\right)}{\∂x_N\∂x_1}& \frac{{\∂}^{2}f\left(x^{\left(k\right)}\right)}{\∂x_N\∂x_2}& \·\·\·& \frac{{\∂}^{2}f\left(x^{\left(k\right)}\right)}{\∂x_N\∂x_N}\end{array}\right)$

其中，N为未知数的个数。 

目标函数为： $f(x,y,z)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(t_i-\frac{1}{v}\×\sqrt{{(x_i-x)}^2+{(y_i-y)}^2+{(z_i-z)}^2})}^2.$

进一步地，所述步骤七中，对于目标函数进行了二级近似后，对新的目标函数求一次导数。由最终得 

$\mathrm{\Δ}{x}^{\left(k\right)}=-\frac{\stackrel{\‾}{g}\left(x^{\left(k\right)}\right)}{H^{\left(k\right)}}=-H^{\left(k\right)-1}\stackrel{\‾}{g}\left(x^{\left(k\right)}\right)$

步骤八：得到迭代公式 $\begin{array}{cc}{x}^{(k+1)}=x^{\left(k\right)}+\mathrm{\Δx}=x^{\left(k\right)}-H^{\left(k\right)-1}\stackrel{\‾}{g}\left(x^{\left(k\right)}\right)& (k=\mathrm{1,2},...),\end{array}$ 再令k＝k+1。 

其中k为迭代次数，H^(k)是目标函数在x^(k)处的黑塞矩阵，同理为目标函数在x^(k)处的梯度。 

设置迭代停止标准，当达到标准时，停止迭代，得到结果，如果没有带到标准返回到步骤六。 

本发明的技术方案，通过网格搜索的方法，先给出地震源的大致位置，这样就会使得初始值与真实值的偏差不会太大，继而通过牛顿法进行反演，既能够得到真实的发震位置，又能进一步的提高微震反演的速度。 

附图说明

图1为反演定位的流程框图；

具体实施方式

以下结合具体的方法实施过程对本发明的原理进行描述，所举实例只用于解释本发明，并非用于限定本发明的范围。 

一种基于网格搜索的牛顿反演定位方法，实施例可以为： 

步骤一：先对划定的区域进行网格的划分，区域的大小，网格的大小，从而划分的点的个数，依据需要进行相应的设置。 

步骤二：代价函数的建立通过发震地点，储层速度，即地震波旅行时的函数关系来确定。具体形式为: 

$t(x,y,z)=\frac{1}{{\mathrm{\ρ}}_T}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(t_i-\frac{1}{v}\×\sqrt{{(x_i-x)}^2+{(y_i-y)}^2+{(z_i-z)}^2})}^2$

其中t_i是地震波的旅行时，x_i,y_i,z_i,t_i是每一个检波器所在的空间位置以及每一个检波器所得到的地震波的到时差，n是检波器的个数，ρ_T是一个常数。v是地下介质速度。 

步骤三：每一个经过网格划分得到的点的位置都会对应一个(x,y,z)的空间坐标，将每一个点的(x,y,z)的空间坐标代入到步骤二中的目标函数，得到目标函数值。 

步骤四：将每个点所对应目标函数的值代入到卡方分布的概率密度函数中，其中卡方分布的密度分布函数的具体形式： 

$p\left(t\right)=\frac{t^{\frac{n}{2}-\frac{3}{2}}{e}^{-\frac{3}{2}}}{\mathrm{\Γ}(\frac{n}{2}-\frac{1}{2})\×2^{\frac{n}{2}-\frac{1}{2}}}$

其中，n是检波器的个数，t是每一个点对应的代价函数值，得到每个点对应的卡方分布密度函数值，并以此值作为权重。Γ的具体形式为： 

$\mathrm{\Γ}\left(n\right)={\∫}_0^{+\∞}{e}^{-n}{x}^{n-1}\mathrm{dx},(n>0)$

公式中e为自然数。 

步骤五：得到每一个坐标点所对应的权重之后，将每一个维度的坐标进行加权求和。最终得到发震地点的大致位置x^(k)。 

步骤六：令目标函数为f(x)在x^(k)位置上进行二阶Taylor级数展开。求得： 

$\begin{array}{c}f\left(x\right)=f\left(x^{\left(k\right)}\right)+\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{\∂f\left(x^{\left(k\right)}\right)}{\∂x_i}\mathrm{\Δ}{x_i}^{\left(k\right)}+\frac{1}{2}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\frac{{\∂}^{2}f\left(x^{\left(k\right)}\right)}{\∂x_i\∂x_j}\mathrm{\Δ}{x_i}^{\left(k\right)}\mathrm{\Δ}{x_j}^{\left(k\right)}\\ =f\left(x^{\left(k\right)}\right)+\stackrel{\‾}{g}{\left(x^{\left(k\right)}\right)}^T\mathrm{\Δ}{x}^{\left(k\right)}+\frac{1}{2}\mathrm{\Δ}{x}^{\left(k\right)T}{H}^{\left(k\right)}\mathrm{\Δ}{x}^k\end{array}$

其中 $\stackrel{\‾}{g}{\left(x^{\left(k\right)}\right)}^T=\left[\begin{array}{cccc}\frac{\∂f\left(x^{\left(k\right)}\right)}{\∂x_1}& \frac{\∂f\left(x^{\left(k\right)}\right)}{\∂x_2}& \·\·\·& \frac{\∂f\left(x^{\left(k\right)}\right)}{\∂x_N}\end{array}\right]$ 为函数在x^(k)的梯度，Δx^(k)T＝[Δx₁ ^(k) Δx₂ ^(k) … Δx_N ^(k)]是模型参数的改正向量。其中，N为未知数的个数。H^(k)为黑塞矩阵，它的具体数学表达式为： 

$H^{\left(k\right)}=\left(\begin{array}{cccc}\frac{{\∂}^{2}f\left(x^{\left(k\right)}\right)}{{\∂x}_1\∂x_1}& \frac{{\∂}^{2}f\left(x^{\left(k\right)}\right)}{\∂x_1\∂x_2}& \·\·\·& \frac{{\∂}^{2}f\left(x^{\left(k\right)}\right)}{\∂x_1\∂x_N}\\ \frac{{\∂}^{2}f\left(x^{\left(k\right)}\right)}{\∂x_2\∂x_1}& \frac{{\∂}^{2}f\left(x^{\left(k\right)}\right)}{\∂x_2\∂x_2}& \·\·\·& \frac{{\∂}^{2}f\left(x^{\left(k\right)}\right)}{\∂x_2\∂x_N}\\ \·& \·& & \·\\ \·& \·& & \·\\ \·& \·& & \·\\ \frac{{\∂}^{2}f\left(x^{\left(k\right)}\right)}{\∂x_N\∂x_1}& \frac{{\∂}^{2}f\left(x^{\left(k\right)}\right)}{\∂x_N\∂x_2}& \·\·\·& \frac{{\∂}^{2}f\left(x^{\left(k\right)}\right)}{\∂x_N\∂x_N}\end{array}\right)$

其中，N为未知数的个数。 

目标函数为： $f(x,y,z)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(t_i-\frac{1}{v}\×\sqrt{{(x_i-x)}^2+{(y_i-y)}^2+{(z_i-z)}^2})}^2.$

步骤七：对于目标函数进行了二级近似后，对新的目标函数求一次导数。 

由 最终得 

$\mathrm{\Δ}{x}^{\left(k\right)}=-\frac{\stackrel{\‾}{g}\left(x^{\left(k\right)}\right)}{H^{\left(k\right)}}=-H^{\left(k\right)-1}\stackrel{\‾}{g}\left(x^{\left(k\right)}\right)$

步骤八：得到迭代公式 

$\begin{array}{cc}{x}^{(k+1)}=x^{\left(k\right)}+\mathrm{\Δx}=x^{\left(k\right)}-H^{\left(k\right)-1}\stackrel{\‾}{g}\left(x^{\left(k\right)}\right)& (k=\mathrm{1,2},...),\end{array}$ 再令k＝k+1。 

其中k为迭代次数，H^(k)是目标函数在x^(k)处的黑塞矩阵，同理为目标函数在x^(k)处的梯度。 

设置迭代停止标准，当达到标准时，停止迭代，得到结果，如果没有带到标准返回到步骤六。 

下面分别采用正演模型、牛顿法、网格搜索法，以及采用本发明的基于网格搜索的牛顿法的反演结果进行反演运算，并对反演结果数据进行比较： 

 一、正演模型

正演模型建立了一个长1000米，宽500，高300米的区域，这是检波器所在的区域，微震源的发生空间位置没有限制。一共有12个检波器，它们的分布位置如表1：

表1 检波器位置

假设是个微震源的发生空间位置（如表2）为：

 二、牛顿法的反演结果

牛顿法反演结果如下（表3）：

 三、网格搜索法的反演结果

网格搜索法的反演结果如表4：

表4 网格搜索法结果

四、基于网格搜索法的牛顿法反演结果

基于网格搜索法的牛顿法反演结果如表5

表5 基于网格搜索法的牛顿法反演结果

下面对上述部分反演方法的结果比较及误差分析：

牛顿法的误差即相关计算如表6：

表6牛顿法的误差即相关计算

基于网格搜索的牛顿法误差分析及相关计算如表7：

表7基于网格搜索的牛顿法误差分析及相关计算

总结分析：

牛顿法的误差以及基于网格搜索的牛顿法的误差分析如表6、表7所示，我们可以看出，在某些位置，牛顿法的精度甚至要比改进的牛顿法更好，但是比较总体，改进的牛顿法的稳定度更好，因为牛顿法很大程度上精度取决于初始值，然而网格搜索法很好的解决了这个问题。而且，改进的牛顿法，最大的误差值只有不到牛顿法最大误差总体的三分之一，所以从这个方面分析，改进的牛顿法实质上可提高了整体计算的精确度的。

此外，比较两中牛顿法的迭代次数，改进的牛顿法很大程度上减少了迭代次数，这样也给反演的过程减少了不少的事件。 

综上，基于网格搜索的牛顿反演方法是一种有效改善牛顿反演方法的一种新型反演方法。 

 以上所述仅为实现本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。  

Claims (9)

1.一种基于网格搜索的牛顿反演定位方法，其特征在于，包括网格搜索部分和牛顿法部分，其中所述网格搜索部分包括： 

步骤一：对于实行微地震观测的区域，在三维空间中，进行网格体的划分，确定每一个维度的坐标值的值域，以及划分的格数； 

步骤二：建立代价函数t(x)； 

步骤三：将每一个网格上的点的坐标，代入到目标函数中去，得到相对应的目标函数值； 

步骤四：将每个点所对应目标函数的值代入到卡方分布的概率密度函数中，得到每个点对应的卡方分布密度函数值，并以此值作为权重； 

步骤五：得到每一个坐标点所对应的权重之后，将每一个维度的坐标进行加权求和，最终得到发震地点的位置x^(k)； 

所述牛顿法部分包括： 

步骤六：令目标函数为f(x)在x^(k)位置上进行二阶Taylor级数展开； 

步骤七：对于目标函数进行了二级近似后，对新的目标函数求一次导数； 

步骤八：得到迭代公式，求得结果。 

2.如权利要求1所述的方法，其特征在于，所述步骤一中，对所述网格体的划分包括：首先设定监测区域的大小，然后设定网格的大小，通过对划定的区域进行网格的划分，从而得到划分的点的个数。 

3.如权利要求1所述的方法，其特征在于，所述步骤二中，通过发震地点、储层速度，即地震波旅行时的函数关系来确定所述代价函数的建立，具体形式为: 

其中t_i是地震波的旅行时，x_i,y_i,z_i,t_i是每一个检波器所在的空间位置以及每一个检波器所得到的地震波的到时差，n是检波器的个数，ρ_T是一个常数。v是地下介质速度。 

4.如权利要求1所述的方法，其特征在于，所述步骤三中，每一个经过网格划分得到的点的位置都会对应一个(x,y,z)的空间坐标，将每一个点的(x,y,z)的空间坐标代入到所述步骤二中的目标函数，得到目标函数值。 

5.如权利要求1所述的方法，其特征在于，所述步骤四中，将每个点所对应的所述目标函数的值代入到卡方分布的概率密度函数中，其中卡方分布的密度分布函数的具体形式： 

其中，n是检波器的个数，t是每一个点对应的代价函数值，得到每个点对应的卡方分布密度函数值，并以此值作为权重。Γ的具体形式为： 

公式中e为自然数。 

6.如权利要求1所述的方法，其特征在于，所述步骤五中，得到每一个坐标点所对应的权重之后，将每一个维度的坐标进行加权求和，最终得到发震地点的大致位置x^(k)。 

7.如权利要求3所述的方法，其特征在于，所述步骤六中，令所述目标函数f(x)在x^(k)位置上进行二阶Taylor级数展开，具体为：， 

其中，x是待求未知量组成的向量， 

为函数在x^(k)的梯度， 

Δx^(k)T＝[Δx₁ ^(k) Δx₂ ^(k) … Δx_N ^(k)]是模型参数的改正向量， 

H^(k)为黑塞矩阵，它的具体数学表达式为： 

N为未知数的个数， 

所述目标函数为：

8.如权利要求7所述的方法，其特征在于，所述步骤七中，对于目标函数进行了二级近似后，对新的目标函数求一次导数，由

最终得

9.如权利要求8所述的方法，其特征在于，所述步骤八中，所述迭代公式为 再令k＝k+1。 

其中k为迭代次数，H^(k)是目标函数在x^(k)处的黑塞矩阵，同理为目标函数在x^(k)处的梯度。 

设置迭代停止标准，当达到标准时，停止迭代，得到结果，如果没有达到迭代标准返回到所述步骤六。