阅读说明：本技术 应力值变化方向约束下的二维矿震全波形反演方法 (Two-dimensional mine earthquake full waveform inversion method under stress value change direction constraint ) 是由 陈卿 赵端 尹洪申 于 2019-11-15 设计创作，主要内容包括：本发明公开了一种应力值变化方向约束下的二维矿震全波形反演方法,属于煤矿井下二维矿震监测数据的全波形反演技术领域。本发明包括：步骤一确定目标区域网格剖分、步骤二将VSB进行二维FFT得m×m矩阵VFB、步骤三确定频域拟保留坐标长度D、将VFBC、VFAC中的自由变量用向量VFBCV、VFACV表示、步骤四确定应力值变化方向和变化位置矩阵SV、步骤五计算目标函数。本发明能提高非均匀岩体全波形反演的可靠性,为井下灾害预测预报提供基础数据支撑。(The invention discloses a two-dimensional mine earthquake full waveform inversion method under stress value change direction constraint, and belongs to the technical field of full waveform inversion of coal mine underground two-dimensional mine earthquake monitoring data. The invention comprises the following steps: determining target area grid subdivision, performing two-dimensional FFT on VSB to obtain an m multiplied by m matrix VFB, determining a frequency domain quasi-reserved coordinate length D, representing free variables in VFBC and VFAC by vectors VFBCV and VFACV, determining a stress value change direction and a change position matrix SV in the fourth step, and calculating a target function in the fifth step. The method can improve the reliability of the full waveform inversion of the inhomogeneous rock mass and provide basic data support for the prediction and forecast of the underground disasters.)

应力值变化方向约束下的二维矿震全波形反演方法

技术领域

本发明涉及煤矿井下二维矿震监测数据的全波形反演技术领域，适用于提高非均匀岩体全波形反演的可靠性。

背景技术

矿震全波形反演可以获取目标区域的波速分布，进而获取井下地质构造、应力分布，目前的研究表明，通过波速反演可以实现井下灾害预测预报。但是，由于观测数据不完整，目前的井下波速反演的可靠性较低，且难以获得精确的反演结果。在弹性形变阶段，岩体所受应力的增大或减小将导致其中的弹性波波速相应增大或减小，因此可使用应力值变化方向约束波速反演过程，以获得更可靠的反演结果。

发明内容

为了克服上述现有技术的不足之处，本发明提供一种应力值变化方向约束下的二维矿震全波形反演方法，能提高非均匀岩体全波形反演的可靠性。

本发明是通过如下技术方案实现的：一种应力值变化方向约束下的二维矿震全波形反演方法，

步骤一：确定目标区域网格剖分，目标区域被剖分为m×m的矩形网格，m为奇数；应力变化前的波速分布矩阵为VSB，应力变化后的波速分布矩阵为VSA，波速分布矩阵中的元素与网格剖分中的每格波速一一对应，波速均为大于0的实数。

步骤二：将VSB进行二维FFT(Fast Fourier Transform)后所得m×m矩阵为VFB，记此过程为VFB＝FFT2(VSB)，逆过程为IFFT2，有VSA＝IFFT2(VFA)；根据二维FFT性质，VFB中元素均为复数，且幅度谱关于坐标轴对称，相位谱关于原点对称，其中心点坐标为

同样，得应力变化后波速分布矩阵的二维FFT矩阵VFA＝FFT2(VSA)。

步骤三：确定频域拟保留坐标长度D，即对于VFB、VFA，位于横坐标范围

或

纵坐标范围

或

的元素均为0，记频域截断后的VFB为VFBC，频域截断后的VFA为VFAC；同样的，VFBC、VFAC的幅度谱也关于坐标轴对称，相位谱也关于原点对称，因此VFBC、VFAC中的自由变量数目为4D²+4D+1；

将VFBC、VFAC中的自由变量用向量VFBCV、VFACV表示，并记此过程为M_TO_V，则有VFBCV＝M_TO_V(VFBC)和VFACV＝M_TO_V(VFAC)，记逆过程为V_TO_M，则有VFBC＝V_TO_M(VFBCV)和VFAC＝V_TO_M(VFACV)。

步骤四：确定应力值变化方向和变化位置矩阵SV

其中(svx_isv，svy_isv)为VSB、VSA中发生应力变化的元素坐标，满足1≤svx_isv≤m和1≤svy_isv≤m且为整数，下标满足1≤nsv≤m²且为整数，svd_isv为应力值变化方向，若应力值变大svd_isv＝1，应力值变小svd_isv＝-1，应力值不变svd_isv＝0。

步骤五、计算目标函数

入参：TruSeiB(应力变化前的真实地震记录)、TruSeiA(应力变化后的真实地震记录)、VFBCV(应力变化前的波速分布向量)、VFACV(应力变化前的波速分布向量)、SV(应力值变化方向和变化位置矩阵)、SSP(震源、检波器位置，地震子波等参数)

返回值：TFV(目标函数值)

计算过程：

获得应力变化前、后的理论波速分布矩阵VSB、VSA

VSB＝IFFT2(V_TO_M(VFBCV))

VSA＝IFFT2(V_TO_M(VFACV))

获得应力变化前、后的理论地震记录：

记常密度声波方程正演为FW，应力变化前的理论地震记录为TheoSeiB，应力变化后的理论地震记录为TheoSeiA，有：

TheoSeiB＝FW(VSB，SSP)

TheoSeiA＝FW(VSA，SSP)

计算目标函数值：

TFV＝CHECK(VSB，VSA，SV)×(||TruSeiB-TheoSeiB||₂+||TruSeiA-TheoSeiA||₂)

CHECK为应力分布检查函数，当VSB、VSA中应力值变化方向与SV中完全一致时，函数值为1，否则函数值为2。

所述步骤一中，若目标区域无法剖分为奇数行或奇数列，应补齐一行或一列以满足m为奇数。

本发明的有益效果是：本发明基于井下波速反演现状及应力值与波速间关系，采用的应力值变化方向约束下的二维矿震全波形反演方法，能提高非均匀岩体全波形反演的可靠性，本发明能提高非均匀岩体全波形反演的可靠性，为井下灾害预测预报提供基础数据支撑，获得更为精确的反演结果。从而能够实现更为准确的井下矿震灾害预测预报，提高矿井作业安全性。二维矿震全波形反演是典型的高维空间反问题，其模型空间维数等于描述目标区域的参数数量，对于本发明所述方法来说，即等于网格个数。因此，为使寻优算法在高维空间中能够可靠地收敛到真实解附近，本发明从两个方向上提出解决方案：一是通过二维傅里叶变换将模型由空域转换到频域，利用模型在频域的稀疏性降低模型维数从而降低模型空间大小；二是以测量到的应力值变化方向来约束反演过程以提高反演正确收敛的概率。从数值实验结果上看，与不含应力值变化方向约束的方法相比，本发明所述方法所得结果与真实值间距离平均减少超过10％，方差降低超过80％，说明本方法切实可行。

附图说明

下面根据附图和实施例对本发明进一步说明。

图1是本发明震源与检波器布置俯视图；

图2是典型地震记录；

图3是本发明网格剖分示意图。

具体实施方式

下面将结合说明书附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅是本发明一部分实施例，而不是全部的实施例。以下对至少一个示例性实施例的描述实际上仅仅是说明性的，决不作为对本发明及其应用或使用的任何限制。基于本发明中的实施例，本领域普通技术人员在没有开展创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

对于本领域技术人员已知的技术、方法和设备可能不作详细讨论，但在适当情况下，所述技术、方法和设备应当被视为授权说明书的一部分。

如图1-图3所示的实施例：

鉴于井下波速反演现状及应力值与波速间关系，本发明提出一种使用应力值变化方向约束二维矿震全波形反演的方法。井下工作面附近震源、检波器及应力计布置俯视图如图1所示，震源、检波器和应力计布设于人员可达区域且位置已知，震源将震波射入目标岩体，震波在目标岩体中传播后被布设于巷道中的检波器接收而形成地震记录，典型地震记录如图2所示，图中横轴为时间纵轴为振幅。

若在某一时刻，图1中所示部分应力计检测到应力值变化，则这些应力值变化方向可作为同时反演应力变化前波速分布图像和应力变化后波速分布图像的约束，反演方法详述如下。

一种应力值变化方向约束下的二维矿震全波形反演方法，

步骤一：确定目标区域网格剖分，目标区域被剖分为m×m的矩形网格，m为奇数，若目标区域无法剖分为奇数行或奇数列，应补齐一行或一列以满足m为奇数。应力变化前的波速分布矩阵为VSB，应力变化后的波速分布矩阵为VSA，波速分布矩阵中的元素与网格剖分中的每格波速一一对应，波速均为大于0的实数。

步骤二：将VSB进行二维FFT(Fast Fourier Transform)后所得m×m矩阵为VFB，记此过程为VFB＝FFT2(VSB)，逆过程为IFFT2，有VSA＝IFFT2(VFA)；根据二维FFT性质，VFB中元素均为复数，且幅度谱关于坐标轴对称，相位谱关于原点对称，其中心点坐标为

同样，得应力变化后波速分布矩阵的二维FFT矩阵VFA＝FFT2(VSA)。

步骤三：确定频域拟保留坐标长度D，即对于VFB、VFA，位于横坐标范围

或纵坐标范围

或

的元素均为0，记频域截断后的VFB为VFBC，频域截断后的VFA为VFAC；同样的，VFBC、VFAC的幅度谱也关于坐标轴对称，相位谱也关于原点对称，因此VFBC、VFAC中的自由变量数目为4D²+4D+1；

将VFBC、VFAC中的自由变量用向量VFBCV、VFACV表示，并记此过程为M_TO_V，则有VFBCV＝M_TO_V(VFBC)和VFACV＝M_TO_V(VFAC)，记逆过程为V_TO_M，则有VFBC＝V_TO_M(VFBCV)和VFAC＝V_TO_M(VFACV)。

步骤四：确定应力值变化方向和变化位置矩阵SV

其中(svx_isv，svy_isv)为VSB、VSA中发生应力变化的元素坐标，满足1≤svx_isv≤m和1≤svy_isv≤m且为整数，下标满足1≤nsv≤m²且为整数，svd_isv为应力值变化方向，若应力值变大svd_isv＝1，应力值变小svd_isv＝-1，应力值不变svd_isv＝0。

步骤五、计算目标函数

入参：TruSeiB(应力变化前的真实地震记录)、TruSeiA(应力变化后的真实地震记录)、VFBCV(应力变化前的波速分布向量)、VFACV(应力变化前的波速分布向量)、SV(应力值变化方向和变化位置矩阵)、SSP(震源、检波器位置，地震子波等参数)

返回值：TFV(目标函数值)

计算过程：

获得应力变化前、后的理论波速分布矩阵VSB、VSA

VSB＝IFFT2(V_TO_M(VFBCV))

VSA＝IFFT2(V_TO_M(VFACV))

获得应力变化前、后的理论地震记录：

记常密度声波方程正演为FW，应力变化前的理论地震记录为TheoSeiB，应力变化后的理论地震记录为TheoSeiA，有：

TheoSeiB＝FW(VSB，SSP)

TheoSeiA＝FW(VSA，SSP)

计算目标函数值：

TFV＝CHECK(VSB，VSA，SV)×(||TruSeiB-TheoSeiB||₂+||TruSeiA-TheoSeiA||₂)

CHECK为应力分布检查函数，当VSB、VSA中应力值变化方向与SV中完全一致时，函数值为1，否则函数值为2。

本发明基于井下波速反演现状及应力值与波速间关系，采用的应力值变化方向约束下的二维矿震全波形反演方法，能提高非均匀岩体全波形反演的可靠性，获得更为精确的反演结果，从而能够实现更为准确的井下矿震灾害预测预报，提高矿井作业安全性。二维矿震全波形反演是典型的高维空间反问题，其模型空间维数等于描述目标区域的参数数量，对于本发明所述方法来说，即等于网格个数。因此，为使寻优算法在高维空间中能够可靠地收敛到真实解附近，本发明从两个方向上提出解决方案：一是通过二维傅里叶变换将模型由空域转换到频域，利用模型在频域的稀疏性降低模型维数从而降低模型空间大小；二是以测量到的应力值变化方向来约束反演过程以提高反演正确收敛的概率。

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