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

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

Info

Publication number
CN110794460A
CN110794460A CN201911124128.8A CN201911124128A CN110794460A CN 110794460 A CN110794460 A CN 110794460A CN 201911124128 A CN201911124128 A CN 201911124128A CN 110794460 A CN110794460 A CN 110794460A
Authority
CN
China
Prior art keywords
stress
change
vsb
vsa
stress value
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN201911124128.8A
Other languages
English (en)
Inventor
陈卿
赵端
尹洪申
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China University of Mining and Technology CUMT
Original Assignee
China University of Mining and Technology CUMT
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China University of Mining and Technology CUMT filed Critical China University of Mining and Technology CUMT
Priority to CN201911124128.8A priority Critical patent/CN110794460A/zh
Publication of CN110794460A publication Critical patent/CN110794460A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/40—Seismology; Seismic or acoustic prospecting or detecting specially adapted for well-logging
    • G01V1/44—Seismology; Seismic or acoustic prospecting or detecting specially adapted for well-logging using generators and receivers in the same well
    • G01V1/48—Processing data
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/40—Seismology; Seismic or acoustic prospecting or detecting specially adapted for well-logging
    • G01V1/44—Seismology; Seismic or acoustic prospecting or detecting specially adapted for well-logging using generators and receivers in the same well
    • G01V1/48—Processing data
    • G01V1/50—Analysing data

Landscapes

  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • Remote Sensing (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明公开了一种应力值变化方向约束下的二维矿震全波形反演方法，属于煤矿井下二维矿震监测数据的全波形反演技术领域。本发明包括：步骤一确定目标区域网格剖分、步骤二将VSB进行二维FFT得m×m矩阵VFB、步骤三确定频域拟保留坐标长度D、将VFBC、VFAC中的自由变量用向量VFBCV、VFACV表示、步骤四确定应力值变化方向和变化位置矩阵SV、步骤五计算目标函数。本发明能提高非均匀岩体全波形反演的可靠性，为井下灾害预测预报提供基础数据支撑。

Description

应力值变化方向约束下的二维矿震全波形反演方法
技术领域
本发明涉及煤矿井下二维矿震监测数据的全波形反演技术领域，适用于提高非均匀岩体全波形反演的可靠性。
背景技术
矿震全波形反演可以获取目标区域的波速分布，进而获取井下地质构造、应力分布，目前的研究表明，通过波速反演可以实现井下灾害预测预报。但是，由于观测数据不完整，目前的井下波速反演的可靠性较低，且难以获得精确的反演结果。在弹性形变阶段，岩体所受应力的增大或减小将导致其中的弹性波波速相应增大或减小，因此可使用应力值变化方向约束波速反演过程，以获得更可靠的反演结果。
发明内容
为了克服上述现有技术的不足之处，本发明提供一种应力值变化方向约束下的二维矿震全波形反演方法，能提高非均匀岩体全波形反演的可靠性。
本发明是通过如下技术方案实现的：一种应力值变化方向约束下的二维矿震全波形反演方法，
步骤一：确定目标区域网格剖分，目标区域被剖分为m×m的矩形网格，m为奇数；应力变化前的波速分布矩阵为VSB，应力变化后的波速分布矩阵为VSA，波速分布矩阵中的元素与网格剖分中的每格波速一一对应，波速均为大于0的实数。
步骤二：将VSB进行二维FFT(Fast Fourier Transform)后所得m×m矩阵为VFB，记此过程为VFB＝FFT2(VSB)，逆过程为IFFT2，有VSA＝IFFT2(VFA)；根据二维FFT性质，VFB中元素均为复数，且幅度谱关于坐标轴对称，相位谱关于原点对称，其中心点坐标为
Figure BDA0002274929690000011
同样，得应力变化后波速分布矩阵的二维FFT矩阵VFA＝FFT2(VSA)。
步骤三：确定频域拟保留坐标长度D，即对于VFB、VFA，位于横坐标范围
Figure BDA0002274929690000012
或
Figure BDA0002274929690000013
纵坐标范围
Figure BDA0002274929690000014
或
Figure BDA0002274929690000015
的元素均为0，记频域截断后的VFB为VFBC，频域截断后的VFA为VFAC；同样的，VFBC、VFAC的幅度谱也关于坐标轴对称，相位谱也关于原点对称，因此VFBC、VFAC中的自由变量数目为4D2+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
Figure BDA0002274929690000021
其中(svxisv，svyisv)为VSB、VSA中发生应力变化的元素坐标，满足1≤svxisv≤m和1≤svyisv≤m且为整数，下标满足1≤nsv≤m2且为整数，svdisv为应力值变化方向，若应力值变大svdisv＝1，应力值变小svdisv＝-1，应力值不变svdisv＝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||2+||TruSeiA-TheoSeiA||2)
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中元素均为复数，且幅度谱关于坐标轴对称，相位谱关于原点对称，其中心点坐标为
Figure BDA0002274929690000041
同样，得应力变化后波速分布矩阵的二维FFT矩阵VFA＝FFT2(VSA)。
步骤三：确定频域拟保留坐标长度D，即对于VFB、VFA，位于横坐标范围
Figure BDA0002274929690000042
或纵坐标范围
Figure BDA0002274929690000044
或
Figure BDA0002274929690000045
的元素均为0，记频域截断后的VFB为VFBC，频域截断后的VFA为VFAC；同样的，VFBC、VFAC的幅度谱也关于坐标轴对称，相位谱也关于原点对称，因此VFBC、VFAC中的自由变量数目为4D2+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
Figure BDA0002274929690000051
其中(svxisv，svyisv)为VSB、VSA中发生应力变化的元素坐标，满足1≤svxisv≤m和1≤svyisv≤m且为整数，下标满足1≤nsv≤m2且为整数，svdisv为应力值变化方向，若应力值变大svdisv＝1，应力值变小svdisv＝-1，应力值不变svdisv＝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||2+||TruSeiA-TheoSeiA||2)
CHECK为应力分布检查函数，当VSB、VSA中应力值变化方向与SV中完全一致时，函数值为1，否则函数值为2。
本发明基于井下波速反演现状及应力值与波速间关系，采用的应力值变化方向约束下的二维矿震全波形反演方法，能提高非均匀岩体全波形反演的可靠性，获得更为精确的反演结果，从而能够实现更为准确的井下矿震灾害预测预报，提高矿井作业安全性。二维矿震全波形反演是典型的高维空间反问题，其模型空间维数等于描述目标区域的参数数量，对于本发明所述方法来说，即等于网格个数。因此，为使寻优算法在高维空间中能够可靠地收敛到真实解附近，本发明从两个方向上提出解决方案：一是通过二维傅里叶变换将模型由空域转换到频域，利用模型在频域的稀疏性降低模型维数从而降低模型空间大小；二是以测量到的应力值变化方向来约束反演过程以提高反演正确收敛的概率。
以上所述仅为本发明的示例性实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (2)

1.一种应力值变化方向约束下的二维矿震全波形反演方法，其特征在于：
步骤一：确定目标区域网格剖分，目标区域被剖分为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，位于横坐标范围
Figure FDA0002274929680000012
或
Figure FDA0002274929680000013
纵坐标范围或
Figure FDA0002274929680000015
的元素均为0，记频域截断后的VFB为VFBC，频域截断后的VFA为VFAC；同样的，VFBC、VFAC的幅度谱也关于坐标轴对称，相位谱也关于原点对称，因此VFBC、VFAC中的自由变量数目为4D2+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
Figure FDA0002274929680000016
其中(svxisv,svyisv)为VSB、VSA中发生应力变化的元素坐标，满足1≤svxisv≤m和1≤svyisv≤m且为整数，下标满足1≤nsv≤m2且为整数，svdisv为应力值变化方向，若应力值变大svdisv＝1，应力值变小svdisv＝-1，应力值不变svdisv＝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||2+||TruSeiA-TheoSeiA||2)
CHECK为应力分布检查函数，当VSB、VSA中应力值变化方向与SV中完全一致时，函数值为1，否则函数值为2。
2.根据权利要求1所述的应力值变化方向约束下的二维矿震全波形反演方法，其特征在于：所述步骤一中，若目标区域无法剖分为奇数行或奇数列，应补齐一行或一列以满足m为奇数。
CN201911124128.8A 2019-11-15 2019-11-15 应力值变化方向约束下的二维矿震全波形反演方法 Pending CN110794460A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911124128.8A CN110794460A (zh) 2019-11-15 2019-11-15 应力值变化方向约束下的二维矿震全波形反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911124128.8A CN110794460A (zh) 2019-11-15 2019-11-15 应力值变化方向约束下的二维矿震全波形反演方法

Publications (1)

Publication Number Publication Date
CN110794460A true CN110794460A (zh) 2020-02-14

Family

ID=69444890

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911124128.8A Pending CN110794460A (zh) 2019-11-15 2019-11-15 应力值变化方向约束下的二维矿震全波形反演方法

Country Status (1)

Country Link
CN (1) CN110794460A (zh)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105765408A (zh) * 2014-10-30 2016-07-13 伊迈格创新技术学院 用于分析位于地下矿场巷道上方的层的地质结构以及相对应力变化的方法及系统
CN106033127A (zh) * 2016-06-29 2016-10-19 中国石油化工股份有限公司 基于横波速度变化率的地应力方位地震预测方法
US20170335675A1 (en) * 2016-05-19 2017-11-23 Stan Lee Method To Predict Pore Pressure And Seal Integrity Using Full Wavefield Inversion
CN109557588A (zh) * 2018-11-16 2019-04-02 徐州工程学院 一种煤矿井下二维矿震波速反演降维方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105765408A (zh) * 2014-10-30 2016-07-13 伊迈格创新技术学院 用于分析位于地下矿场巷道上方的层的地质结构以及相对应力变化的方法及系统
US20170335675A1 (en) * 2016-05-19 2017-11-23 Stan Lee Method To Predict Pore Pressure And Seal Integrity Using Full Wavefield Inversion
CN106033127A (zh) * 2016-06-29 2016-10-19 中国石油化工股份有限公司 基于横波速度变化率的地应力方位地震预测方法
CN109557588A (zh) * 2018-11-16 2019-04-02 徐州工程学院 一种煤矿井下二维矿震波速反演降维方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
陈卿，等: "矿震波速反演模型空间降维方法", 《工矿自动化》 *
陈卿: "基于降维的矿震反演研究", 《中国博士学位论文全文数据库》 *

Similar Documents

Publication Publication Date Title
US20180292553A1 (en) Method and Apparatus for Separating Seismic Diffracted Wave
Tromp et al. Noise cross-correlation sensitivity kernels
CN114779330B (zh) 一种基于微震监测的采掘工作面主裂隙方位分析预测方法
US11119232B2 (en) System and method for real-time passive seismic event localization
Furumura et al. Large scale parallel simulation and visualization of 3D seismic wavefield using the Earth Simulator
CN109375252A (zh) 考虑不同发震构造最大可信地震的地震动参数评价方法
CN109375253A (zh) 基于全部发震构造最大可信地震的地震动参数评价方法
CN108051855B (zh) 一种基于拟空间域声波方程的有限差分计算方法
CN108108331A (zh) 一种基于拟空间域弹性波方程的有限差分计算方法
Homma et al. A physics-based Monte Carlo earthquake disaster simulation accounting for uncertainty in building structure parameters
CN115508908A (zh) 一种地震面波走时和重力异常联合反演方法与系统
Goebel et al. Detecting significant stress drop variations in large micro-earthquake datasets: A comparison between a convergent step-over in the San Andreas Fault and the Ventura thrust fault system, Southern California
Muratov et al. Grid-characteristic method on unstructured tetrahedral meshes
Voronina Recovering a tsunami source and designing an observational system based on an r-solution method
Kirezci et al. Probabilistic assessment of rogue wave occurrence in directional wave fields
Sun et al. Preliminary analysis of earthquake probability based on the synthetic seismic catalog
CN110794460A (zh) 应力值变化方向约束下的二维矿震全波形反演方法
He et al. Interseismic kinematics along the Tuolaishan-Lenglongling fault determined by Sentinel-1 InSAR observations
Lehmann et al. Synthetic ground motions in heterogeneous geologies: the HEMEW-3D dataset for scientific machine learning
Smaglichenko Modification of Gaussian elimination for the complex system of seismic observations
Magrin Multi-scale seismic hazard scenarios
Gallovic et al. High-frequency directivity in strong ground motion modeling methods
CN116243379B (zh) 一种基于震源机制与定位误差校准的强矿震预测方法
CN113221410A (zh) 一种逐元并行强地面运动快速数值模拟方法
Molodensky et al. Temporal variations in the tidal response of the medium in the vicinities of the sources of catastrophic earthquakes

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
RJ01 Rejection of invention patent application after publication
RJ01 Rejection of invention patent application after publication

Application publication date: 20200214