一种煤矿井下二维矿震波速反演降维方法与流程

本发明涉及矿井安全的技术领域，尤其涉及一种煤矿井下二维矿震波速反演降维方法。

背景技术：

波速分布与井下应力分布以及强震分布强相关，因此通过实时监测目标区域的波速分布可以实现矿井灾害的预测预报。现有的波速反演方法可以分为线性方法和非线性方法。线性方法收敛速度快而且计算量小，但是受限于初值选取，容易陷入局部最优而无法获得全局最优解。非线性方法与初值选取无关，适应性更好，但是收敛速度慢而且计算量大，实践中不便应用。

技术实现要素：

为解决现有技术存在的局限和缺陷，本发明提供一种煤矿井下二维矿震波速反演降维方法，包括：

使用直线将目标区域等间隔划分成m行n列的正方形网格，其中m和n为奇数，其中待反演区域为pt，模型空间mt的维数为mtnum＝m×n；

获得目标区域的待反演向量v＝(v11，v12，…，v1n，v21，v22，…，vij，…，vmn)，所述目标区域的每个正方形的介质均匀，所述目标区域的波速为vij，i∈[1，m]，j∈[1，n]；

根据所述待反演向量获得所述目标区域的每个正方形的波速上限向量和波速下限向量，所述波速上限向量和所述波速下限向量分别为：

获得降维之后的待反演变量mod_angvector＝(modvector，angvector)，其中

所述待反演区域pt的频域为pf，所述频域pf的维数为mfnum＝mtnum＝m×n，所述待反演变量mod_angvector的维数为mfcnum，根据modvector和angvector形成的频域为pfc，频域pf的中心点坐标为频域pf的最远点坐标为

获得所述待反演变量mod_angvector的上限向量为其中modvector和angvector的每个元素上限构成的向量分别为：

获得所述待反演变量mod_angvector的下限向量为其中modvector和angvector的每个元素下限构成的向量分别为：

当||tobs-ttheo||2最小时，根据所述上限向量和所述下限向量获得待反演变量mod_angvector，其中tobs为测得到时差向量，ttheo为理论到时差向量。

可选的，所述根据所述待反演向量获得所述目标区域的每个正方形的波速上限向量和波速下限向量的步骤之后，所述获得降维之后的待反演变量的步骤之前，包括：

对所述待反演向量、所述波速上限向量、所述波速下限向量进行归一化处理，获得归一化之后的待反演向量、波速下限向量、波速上限向量分别为：

v＝(v11，v12，…，v1n，v21，v22，…，vij，…，vmn)

根据归一化之后的待反演向量、波速上限向量、波速下限向量获得归一化系数为：

其中，min(v^min)为向量v^min中最小的元素，max(v^max)为向量v^max中最大的元素，min(v^min)为向量v^min中最小的元素，max(v^max)为向量v^max中最大的元素。

可选的，还包括：

根据所述归一化系数获得如下公式：

vij＝k^norm·vij(2)

根据上述公式(2)获得波速归一化公式为：

v＝k^norm·v(3)。

可选的，所述频域pf之中拟保留的变量数目为

mfcnum＝4d²+4d+1。

可选的，所述理论到时差向量为：

ttheo＝ray_trace(fft_2d^-1(map_to_vector^-1(mod_angvector)))(4)

其中，fft_2d^-1为根据所述待反演区域pt获得所述频域pf的逆过程，map_to_vector^-1为根据所述频域pfc获得待反演变量mod_angvector的逆过程，ray_trace为根据射线追踪获得理论到时差向量ttheo。

本发明具有下述有益效果：

本发明提供的煤矿井下二维矿震波速反演降维方法，包括：将目标区域等间隔划分成多个正方形网格，获得正方形网格的波速上限向量和波速下限向量，获得拟保留的待反演变量和变量范围，降低反演模型空间维数，在降维之后的模型空间内寻找最优解。本发明提供的技术方案可以解决反演过程中的初值依赖问题，从而降低了反演结果收敛至错误位置的可能性。本发明提供的技术方案在降维之后的模型空间内寻找最优解，从而降低了反演寻优的模型空间大小，提高了反演结果的可靠度以及反演计算的速度。

附图说明

图1为本发明实施例一提供的目标区域的网格划分示意图。

图2(a)为本发明实施例一提供的波速分布图像。

图2(b)为本发明实施例一提供的波速分布图像的第一幅度谱。

图2(c)为本发明实施例一提供的波速分布图像的第二幅度谱。

图3为本发明实施例一提供的频域pf中心点附近的坐标分布示意图。

具体实施方式

为使本领域的技术人员更好地理解本发明的技术方案，下面结合附图对本发明提供的煤矿井下二维矿震波速反演降维方法进行详细描述。

实施例一

本实施例提供一种煤矿井下二维矿震波速反演降维方法，从而降低反演寻优的模型空间大小，提高反演结果的可靠度以及反演计算速度。本实施例提供的降维方法如下：

本实施例首先获得目标区域的网格剖分。图1为本发明实施例一提供的目标区域的网格划分示意图。如图1所示，本实施例使用直线将目标区域等间隔剖分成m行n列的正方形网格，m和n均为奇数，若剖分完成后m或n为偶数，可将最后一行或最后一列删除以形成奇数行和奇数列，或者复制最后一行或最后一列将图像补齐为奇数行和奇数列。

本实施例将剖分后的待反演区域设定为pt，模型空间mt的维数mtnum＝m×n，模型空间是指每个网格中所有可能的波速组成的波速分布图像全体。假定每个正方形的介质均匀，设定其波速为vij，其中i∈[1,m]，j∈[1,n]，则待反演向量为v＝(v11,v12,l,v1n,v21,v22,l,vij,l,vmn)。

本实施例获得网格的波速上下限，并对其进行归一化，获得波速归一化公式。具体来说，本实施例根据实际情况，获得每个正方形的波速上限向量和波速下限向量，所述波速上限向量和所述波速下限向量分别为：

所述波速上限向量和所述波速下限向量限制了待反演向量v的范围，即为了方便反演程序的设计，本实施例需要将待反演变量v、波速下限向量v^min、波速上限向量v^max归一化到统一的上下限范围之内。本实施例提供的技术方案可以解决反演过程中的初值依赖问题，从而降低了反演结果收敛至错误位置的可能性。

本实施例获得归一化后的待反演变量、波速下限向量、波速上限向量分别为：

v＝(v11,v12,l,v1n,v21,v22,l,vij,l,vmn)

其中，

本实施例进行如下设定：向量v^min中最小的元素为min(v^min)，向量v^max中最大的元素为max(v^max)，向量v^min中最小的元素为min(v^min)，向量v^max中最大的元素为max(v^max)，获得归一化系数为：

因此，本实施例可以根据所述归一化系数获得

vij＝k^norm·vij(2)

也就是说，本实施例获得波速归一化公式为：

v＝k^norm·v(3)

本实施例获得拟保留的待反演变量和变量范围，从而降低反演模型空间维数。具体来说，若认为图1是一幅二维图像，则v的每一个元素都可以认为是图1的每一个像素的灰度，对图1进行离散傅氏变换的快速算法(fastfouriertransformation，fft)可以获得对应的频域图像。图2(a)为本发明实施例一提供的波速分布图像，图2(b)为本发明实施例一提供的波速分布图像的第一幅度谱，图2(c)为本发明实施例一提供的波速分布图像的第二幅度谱。可以看出，图2(a)为波速分布图像，图2(b)为45度视图的幅度谱，图2(c)为俯视图的幅度谱。

显然，从频域上描述一幅图像并不需要完全保留mtnum维的模型，只需对图像进行频域截断，仅保留相对低频成分就足以保留图像的绝大部分信息。因此，本实施例提供的模型空间降维的过程即是使用理想低通滤波器对待反演区域pt进行理想低通滤波的过程。本实施例提供的模型空间降维过程如下：

本实施例将待反演区域pt的频域表示为pf，则频域pf的维数同样为mfnum＝mtnum＝m×n。图3为本发明实施例一提供的频域pf中心点附近的坐标分布示意图。获得频域pf的中心点坐标为则频域pf中心点附近部分坐标如图3所示。

获得频域pf中拟保留的最远点坐标为由于幅度谱关于中心点偶对称，相位谱关于中心点奇对称，而且中心点的相位恒为0，因此频域pf中拟保留的变量数目mfcnum＝4d²+4d+1。

获得降维后的待反演变量为mod_angvector＝(modvector,angvector)，待反演变量为模向量和幅角向量的组合，其中

向量modvector之中元素的下标表示该模在图3中的坐标；

向量angvector之中元素的下标表示该幅角在图3中坐标。

由于恒为0，若不计则mod_angvector的维数等于mfcnum。

获得mod_angvector上限向量为其中：

为modvector中每个元素上限构成的向量；

为angvector中每个元素上限构成的向量，而且恒为0。

获得mod_angvector下限向量为其中：

为modvector中每个元素下限构成的向量；

为angvector中每元素下限构成的向量，而且恒为0。

模的上限向量元素中的最大值和下限向量元素中的最小值取决于v^max和v^min，或者根据波速分布情况进行设定。幅角的上限向量元素中的最大值取π或者根据实际情况进行设定，幅角的下限向量元素中的最小值取-π或者根据实际情况进行设定。

本实施例将仅由modvector和angvector生成的频域表示为pfc，其中位于拟保留坐标范围之外的模和幅角设置为0。本实施例中，由频域pfc获得mod_angvector的过程为map_to_vector，逆过程为map_to_vector^-1，则有mod_angvector＝map_to_vector(pfc)，pfc＝map_to_vector^-1(mod_angvector)。由频域pf获得频域pfc的过程为low_pass，则有pfc＝low_pass(pf)。由待反演区域pt获得频域pf的过程为fft_2d，逆过程为fft_2d^-1，则有pf＝fft_2d(pt)，pt＝fft_2d^-1(pf)。

本实施例在降维后的模型空间内寻优以求解问题。具体来说，本实施例获得测得到时差向量为tobs以及理论到时差向量为ttheo。对于一个已知的波速分布，理论到时差向量ttheo可以由射线追踪获得，此过程为ray_trace，则反演目标区域波速等价于一个有约束最优化问题。本实施例可以将上述问题转化为寻找一个最优的mod_angvector，使得||tobs-ttheo||2最小，其中理论到时差为：

ttheo＝ray_trace(fft_2d^-1(map_to_vector^-1(mod_angvector)))(4)

本实施例中，约束上限向量为约束下限向量为最优化方法不唯一，可以根据需要选择遗传算法(geneticalgorithm，ga)、粒子群算法(particleswarmoptimization，pso)等常用算法。

本实施例提供的煤矿井下二维矿震波速反演降维方法，包括：将目标区域等间隔划分成多个正方形网格，获得正方形网格的波速上限向量和波速下限向量，获得拟保留的待反演变量和变量范围，降低反演模型空间维数，在降维之后的模型空间内寻找最优解。本实施例提供的技术方案可以解决反演过程中的初值依赖问题，从而降低了反演结果收敛至错误位置的可能性。本实施例提供的技术方案在降维之后的模型空间内寻找最优解，从而降低了反演寻优的模型空间大小，提高了反演结果的可靠度以及反演计算的速度。

可以理解的是，以上实施方式仅仅是为了说明本发明的原理而采用的示例性实施方式，然而本发明并不局限于此。对于本领域内的普通技术人员而言，在不脱离本发明的精神和实质的情况下，可以做出各种变型和改进，这些变型和改进也视为本发明的保护范围。