CN110687593A - 二维小波域矿震监测数据反演方法 - Google Patents

二维小波域矿震监测数据反演方法 Download PDF

Info

Publication number
CN110687593A
CN110687593A CN201910971036.7A CN201910971036A CN110687593A CN 110687593 A CN110687593 A CN 110687593A CN 201910971036 A CN201910971036 A CN 201910971036A CN 110687593 A CN110687593 A CN 110687593A
Authority
CN
China
Prior art keywords
vector
target area
twc
coefficients
wavelet
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.)
Granted
Application number
CN201910971036.7A
Other languages
English (en)
Other versions
CN110687593B (zh
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 CN201910971036.7A priority Critical patent/CN110687593B/zh
Publication of CN110687593A publication Critical patent/CN110687593A/zh
Application granted granted Critical
Publication of CN110687593B publication Critical patent/CN110687593B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

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/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/282—Application of seismic models, synthetic seismograms
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis

Landscapes

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

Abstract

本发明公布一种二维小波域矿震监测数据反演方法，属于煤矿矿震监测技术领域。步骤一，目标区域网格剖分；步骤二，建立目标区域波速分布的1行m2列的Haar小波域描述向量TW，Haar小波与消失矩为1的Daubechies小波等价；步骤三，确定目标区域分解尺度DeSc；步骤四，进一步降低TW的维数，记降维后的TW为TWC；步骤五，目标函数；①，由反演向量TWC得波速分布TS；②，计算理论地震记录；③，返回目标函数值。通过本发明的反演方法，可以改善由于观测数据不完整、正演模型与实际情况不符等原因造成的目标函数多峰值形态问题，便于寻优算法正确收敛到全局最优解附近，能够提高非均匀岩体波速分布全波形反演的可靠度和速度。

Description

二维小波域矿震监测数据反演方法
技术领域
本发明涉及煤矿矿震监测技术领域，具体是一种二维小波域矿震监测数据反演方法，用于提高非均匀岩体波速分布全波形反演的可靠度和速度。
背景技术
波速分布与井下应力分布及强震分布强相关，因此，通过实时监测目标区域的波速分布可以实现矿井灾害的预测预报。井下工作面附近震源与检波器布置俯视图如图1所示，震源、检波器布设位置已知，震源将震波射入目标岩体，震波在目标岩体中传播后被布设于巷道中的检波器接收而形成地震记录，典型地震记录如图2所示，图中横轴为时间，纵轴为振幅。
获取目标岩体波速分布在数学上等价于一个最优化问题：
min||dTrue-dTheo||2，寻优对象为目标岩体波速分布mTheo，dTrue为实测地震记录矩阵，dTheo为理论地震记录矩阵，波速分布与地震记录之间通过映射d＝F(m)建立关系，F代表常密度声波方程正演、变密度声波方程正演等正演方法。
但是，由于观测数据不完整、正演模型与实际情况不符等原因，目标函数 ||dTrue-dTheo||2通常呈多峰值形态，导致寻优算法无法正确收敛到全局最优解附近。
发明内容
为解决目标函数||dTrue-dTheo||2不易全局收敛的问题，考虑到矿震监测的连续性，本发明提供一种二维小波域矿震监测数据反演方法。
本发明通过以下技术方案实现：一种二维小波域矿震监测数据反演方法，
步骤一，目标区域网格剖分；
将目标区域剖分为m×m的网格，每网格中的波速为定值，变量m满足： m＝2n，n为大于0的整数，记目标区域波速分布为矩阵TS；
步骤二，建立目标区域波速分布的1行m2列的Haar小波域描述向量TW， Haar小波与消失矩为1的Daubechies小波等价；
步骤三，确定目标区域分解尺度DeSc，以获得向量TW的格式，其中 1<DeSc<n且为整数；
尺度DeSc包含尺度系数、水平、垂直、斜向小波系数，其余尺度仅包含小波系数，每组系数长度标于TW格式下方，TW总长度为22n＝m2；
步骤四，定义尺度系数拟保留变量数
Figure BDA0002231262910000024
确定尺度 1拟保留变量数DePaToKpMinS，要求2≤DePaToKpMinS<DePaToKpMaxS且为整数；
任一尺度TS保留的变量数PaNumToKepTS由公式：
获得，其中，1<TS<DeSc；
公式中min表示取两者中的较小值，round表示按照四舍五入规律取整数，由此向量TW的维数得以进一步降低，记降维后的TW为TWC；
步骤五，目标函数；
入参：TruSei1，TruSei1是真实地震记录1；TruSei2，TruSei2是真实地震记录2；TWC1，TWC1是待反演向量1；TWC2，TWC2是待反演向量2； SenSouPar，SenSouPar是检波器与震源参数；
返回值：TFV；
①，由反演向量TWC得波速分布TS；
TWC1、TWC2是TWC的一个实例，用于生成波速分布TS的实例TS1、TS1，
记由TWC获得TW的过程为：TW＝WAVC_TO_WAV(TWC)，
记由TW获得TS的过程为：TS＝WAV_TO_V(TW)；
WAV_TO_V为以Haar小波为基的小波逆变换；
WAVC_TO_WAV等价于将TWC中的每组系数转换为TW中对应的系数；
记TWC中的某一组系数为ZC＝[ZC1,ZC2,L,ZCic,L,ZCnc]，
记TW中某一组系数为Z＝[Z1,Z2,L,Ziw,L,Znw]。
记ZC的横坐标向量为
Figure BDA0002231262910000022
斜率向量为KC＝[KC1,KC1,L,KCic,L,KCnc]，其中KC1＝1，
Figure BDA0002231262910000023
则Z中元素与ZC中元素关系为：Z1＝ZC1，
当
Figure BDA0002231262910000031
时，
Figure BDA0002231262910000032
②，计算理论地震记录；
记理论地震记录1为TheoSei1，理论地震记录2为TheoSei2，常密度声波方程正演为FW，则TheoSei1＝FW(TS1,SenSouPar)， TheoSei2＝FW(TS2,SenSouPar)；
③，返回目标函数值；
目标函数值计算公式为：
TFV＝||TheoSei1-TruSei1||2+||TheoSei2-TruSei2||2+λ||TS1-TS2||2
寻优目标是找到最优的TWC1和TWC2使得TFV最小。
优选的：步骤一中，若目标区域无法被剖分为m×m的网格，则先对目标区域进行填充，使目标区域满足剖分条件之后，再将目标区域分为m×m的网格；填充区域的波速根据实际情况确定。
优选的：步骤五中，λ为可变参数，计算方法如下：
记L1＝floor(log10(||TheoSei1-TruSei1||2+||TheoSei2-TruSei2||2))
记L2＝floor(log10(||TS1-TS2||2))，其中，floor表示向下取整；
若L1＝L2，则λ＝1，若L1>L2，则λ＝10L1-L2，若L1<L2，则λ＝10L2-L1。
与现有技术相比，本发明的有益效果是：改善由于观测数据不完整、正演模型与实际情况不符等原因造成的目标函数多峰值形态问题，便于寻优算法正确收敛到全局最优解附近，能够提高非均匀岩体波速分布全波形反演的可靠度和速度。
附图说明
图1是震源与检波器布置俯视图；
图2是典型地震记录；
图3是目标区域网格剖分示意图；
图4是实施例中TW格式；
图5是实施例中TWC格式；
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。
一种二维小波域矿震监测数据反演方法，
步骤一，目标区域网格剖分，结合图1所示；
将目标区域剖分为m×m的网格，每网格中的波速为定值，变量m满足：
m＝2n，n为大于0的整数，记目标区域波速分布为矩阵TS；
若目标区域无法被剖分为m×m的网格，则先对目标区域进行填充，使目标区域满足剖分条件之后，再将目标区域分为m×m的网格；填充区域的波速根据实际情况确定。比如在图3中，若需要在左侧填充一列波速，则填充的这一列波速可设置成与已有波速最左侧一列波速一致。若需要在下方填充一行波速，则填充的这一行波速可设置成与已有波速最下一行波速一致。
步骤二，建立目标区域波速分布的1行m2列的Haar小波域描述向量TW， Haar小波与消失矩为1的Daubechies小波等价；
步骤三，确定目标区域分解尺度DeSc，以获得向量TW的格式，其中 1<DeSc<n且为整数；
尺度DeSc包含尺度系数、水平、垂直、斜向小波系数，其余尺度仅包含小波系数，每组系数长度标于TW格式下方，TW总长度为22n＝m2；
步骤四，定义尺度系数拟保留变量数
Figure BDA0002231262910000041
确定尺度 1拟保留变量数DePaToKpMinS，要求2≤DePaToKpMinS<DePaToKpMaxS且为整数；
任一尺度TS保留的变量数PaNumToKepTS由公式：
Figure BDA0002231262910000042
获得，其中，1<TS<DeSc；
公式中min表示取两者中的较小值，round表示按照四舍五入规律取整数，由此向量TW的维数得以进一步降低，记降维后的TW为TWC；
步骤五，目标函数；
入参：TruSei1，TruSei1是真实地震记录1；TruSei2，TruSei2是真实地震记录2；TWC1，TWC1是待反演向量1；TWC2，TWC2是待反演向量2；
SenSouPar，SenSouPar是检波器与震源参数；
返回值：TFV；
①，由反演向量TWC得波速分布TS；
TWC1、TWC2是TWC的一个实例，用于生成波速分布TS的实例TS1、TS1，
记由TWC获得TW的过程为：TW＝WAVC_TO_WAV(TWC)，
记由TW获得TS的过程为：TS＝WAV_TO_V(TW)；
WAV_TO_V为以Haar小波为基的小波逆变换；
WAVC_TO_WAV等价于将图5所示TWC中的每组系数(即D1′、V1′、 H1′、……、A′DeSc)转换为图4所示TW中对应的系数(即D1、V1、H1、……、 ADeSc)；
记TWC中的某一组系数为ZC＝[ZC1,ZC2,L,ZCic,L,ZCnc]，
记TW中某一组系数为Z＝[Z1,Z2,L,Ziw,L,Znw]。
记ZC的横坐标向量为
Figure BDA0002231262910000051
斜率向量为KC＝[KC1,KC1,L,KCic,L,KCnc]，其中KC1＝1，
Figure BDA0002231262910000052
则Z中元素与ZC中元素关系为：Z1＝ZC1，
当时，
Figure BDA0002231262910000054
②，计算理论地震记录；
记理论地震记录1为TheoSei1，理论地震记录2为TheoSei2，常密度声波方程正演为FW，则TheoSei1＝FW(TS1,SenSouPar)， TheoSei2＝FW(TS2,SenSouPar)；
③，返回目标函数值；
目标函数值计算公式为：
TFV＝||TheoSei1-TruSei1||2+||TheoSei2-TruSei2||2+λ||TS1-TS2||2
寻优目标是找到最优的TWC1和TWC2使得TFV最小；
λ为可变参数，计算方法如下：
记L1＝floor(log10(||TheoSei1-TruSei1||2+||TheoSei2-TruSei2||2))
记L2＝floor(log10(||TS1-TS2||2))，其中，floor表示向下取整；
若L1＝L2，则λ＝1，若L1>L2，则λ＝10L1-L2，若L1<L2，则λ＝10L2-L1。
本实施例改善由于观测数据不完整、正演模型与实际情况不符等原因造成的目标函数多峰值形态问题，便于寻优算法正确收敛到全局最优解附近，能够提高非均匀岩体波速分布全波形反演的可靠度和速度。本发明利用模型小波域表示的稀疏性，将反演模型空间(即反演对象的所有可能值构成的空间)由空域(即图 3中每个网格中的所有可能的波速构成的空间)转换到小波域，并通过折线拟合进一步降低了模型小波域系数的数量，极大降低了待寻优模型空间的维数和大小，使反演计算更快收敛到全局最优解附近以获得正确反演结果。
尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由所附权利要求及其等同物限定。

Claims (3)

1.一种二维小波域矿震监测数据反演方法，
步骤一，目标区域网格剖分；
将目标区域剖分为m×m的网格，每网格中的波速为定值，变量m满足：m＝2n，n为大于0的整数，记目标区域波速分布为矩阵TS；
步骤二，建立目标区域波速分布的1行m2列的Haar小波域描述向量TW，Haar小波与消失矩为1的Daubechies小波等价；
步骤三，确定目标区域分解尺度DeSc，以获得向量TW的格式，其中1<DeSc<n且为整数；
尺度DeSc包含尺度系数、水平、垂直、斜向小波系数，其余尺度仅包含小波系数，每组系数长度标于TW格式下方，TW总长度为22n＝m2；
步骤四，定义尺度系数拟保留变量数DePaToKpMaxS＝22(n-DeSc)，确定尺度1拟保留变量数DePaToKpMinS，要求2≤DePaToKpMinS<DePaToKpMaxS且为整数；
任一尺度TS保留的变量数PaNumToKepTS由公式：
Figure FDA0002231262900000011
获得，其中，1<TS<DeSc；
公式中min表示取两者中的较小值，round表示按照四舍五入规律取整数，由此向量TW的维数得以进一步降低，记降维后的TW为TWC；
步骤五，目标函数；
入参：TruSei1，TruSei1是真实地震记录1；TruSei2，TruSei2是真实地震记录2；TWC1，TWC1是待反演向量1；TWC2，TWC2是待反演向量2；
SenSouPar，SenSouPar是检波器与震源参数；
返回值：TFV；
①，由反演向量TWC得波速分布TS；
TWC1、TWC2是TWC的一个实例，用于生成波速分布TS的实例TS1、TS1，
记由TWC获得TW的过程为：TW＝WAVC_TO_WAV(TWC)，
记由TW获得TS的过程为：TS＝WAV_TO_V(TW)；
WAV_TO_V为以Haar小波为基的小波逆变换；
WAVC_TO_WAV等价于将TWC中的每组系数转换为TW中对应的系数；
记TWC中的某一组系数为ZC＝[ZC1,ZC2,L,ZCic,L,ZCnc]，
记TW中某一组系数为Z＝[Z1,Z2,L,Ziw,L,Znw]。
记ZC的横坐标向量为
Figure FDA0002231262900000021
斜率向量为KC＝[KC1,KC1,L,KCic,L,KCnc]，其中KC1＝1，则Z中元素与ZC中元素关系为：Z1＝ZC1，
当
Figure FDA0002231262900000023
时，
Figure FDA0002231262900000024
②，计算理论地震记录；
记理论地震记录1为TheoSei1，理论地震记录2为TheoSei2，记常密度声波方程正演为FW，则TheoSei1＝FW(TS1,SenSouPar)，TheoSei2＝FW(TS2,SenSouPar)；
③，返回目标函数值；
目标函数值计算公式为：
TFV＝||TheoSei1-TruSei1||2+||TheoSei2-TruSei2||2+λ||TS1-TS2||2
寻优目标是找到最优的TWC1和TWC2使得TFV最小。
2.根据权利要求1所述的二维小波域矿震监测数据反演方法，其特征在于：
步骤一中，若目标区域无法被剖分为m×m的网格，则先对目标区域进行填充，使目标区域满足剖分条件之后，再将目标区域分为m×m的网格；填充区域的波速根据实际情况确定。
3.根据权利要求1所述的二维小波域矿震监测数据反演方法，其特征在于：
步骤五中，λ为可变参数，计算方法如下：
记L1＝floor(log10(||TheoSei1-TruSei1||2+||TheoSei2-TruSei2||2))
记L2＝floor(log10(||TS1-TS2||2))，其中，floor表示向下取整；
若L1＝L2，则λ＝1，若L1>L2，则λ＝10L1-L2，若L1<L2，则λ＝10L2-L1。
CN201910971036.7A 2019-10-12 2019-10-12 二维小波域矿震监测数据反演方法 Active CN110687593B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910971036.7A CN110687593B (zh) 2019-10-12 2019-10-12 二维小波域矿震监测数据反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910971036.7A CN110687593B (zh) 2019-10-12 2019-10-12 二维小波域矿震监测数据反演方法

Publications (2)

Publication Number Publication Date
CN110687593A true CN110687593A (zh) 2020-01-14
CN110687593B CN110687593B (zh) 2021-08-17

Family

ID=69112486

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910971036.7A Active CN110687593B (zh) 2019-10-12 2019-10-12 二维小波域矿震监测数据反演方法

Country Status (1)

Country Link
CN (1) CN110687593B (zh)

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20140039799A1 (en) * 2012-08-06 2014-02-06 Christine E. Krohn Seismic Inversion for Formation Properties and Attentuation Effects
CN105022091A (zh) * 2015-08-07 2015-11-04 中国矿业大学 一种无预测速的远场震源快速定位方法
CN105785436A (zh) * 2016-03-17 2016-07-20 北京矿冶研究总院 一种矿用微震监测方法
CN108549100A (zh) * 2018-01-11 2018-09-18 吉林大学 基于非线性高次拓频的时间域多尺度全波形反演方法
CN109557588A (zh) * 2018-11-16 2019-04-02 徐州工程学院 一种煤矿井下二维矿震波速反演降维方法
CN110018517A (zh) * 2019-05-07 2019-07-16 西安石油大学 一种多尺度地面微地震逆时干涉定位方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20140039799A1 (en) * 2012-08-06 2014-02-06 Christine E. Krohn Seismic Inversion for Formation Properties and Attentuation Effects
CN105022091A (zh) * 2015-08-07 2015-11-04 中国矿业大学 一种无预测速的远场震源快速定位方法
CN105785436A (zh) * 2016-03-17 2016-07-20 北京矿冶研究总院 一种矿用微震监测方法
CN108549100A (zh) * 2018-01-11 2018-09-18 吉林大学 基于非线性高次拓频的时间域多尺度全波形反演方法
CN109557588A (zh) * 2018-11-16 2019-04-02 徐州工程学院 一种煤矿井下二维矿震波速反演降维方法
CN110018517A (zh) * 2019-05-07 2019-07-16 西安石油大学 一种多尺度地面微地震逆时干涉定位方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
ODILE GAUTHIER 等: "Two-dimensional nonlinear inversion of seismic waveforms:Numerical results", 《GEOPHYSICS》 *
孙兴林 等: "基于Matlab的矿震信号小波分析", 《煤矿安全》 *
李铁等: "基于震源机制解的矿井采动应力场反演与应用 ", 《岩石力学与工程学报》 *
陈卿 等: "矿震波速反演模型空间降维方法", 《工矿自动化》 *

Also Published As

Publication number Publication date
CN110687593B (zh) 2021-08-17

Similar Documents

Publication Publication Date Title
CN111239802B (zh) 基于地震反射波形和速度谱的深度学习速度建模方法
EP2869096B1 (en) Systems and methods of multi-scale meshing for geologic time modeling
CN104122585B (zh) 基于弹性波场矢量分解与低秩分解的地震正演模拟方法
CN110031895B (zh) 一种基于图像缝合的多点地质统计学随机反演方法及装置
CN108508482A (zh) 一种地下裂缝地震散射响应特征模拟方法
CN104977607A (zh) 利用变步长网格声波波场模拟的时间域全波形反演方法
CN109239773A (zh) 一种高阶模式瑞雷波的重建方法
CN106249297A (zh) 基于信号估计的水力压裂微地震震源定位方法及系统
CN115508908A (zh) 一种地震面波走时和重力异常联合反演方法与系统
CN104570095A (zh) 一种基于Radon变换消除斜缆虚反射的方法
CN110687593B (zh) 二维小波域矿震监测数据反演方法
CN109145467B (zh) 一种适用于台风区域的沙波运移预测方法
Hu et al. Numerical simulation of barotropic tides around Taiwan
CN104123449B (zh) 复杂山地区域的分区局部变加密不等距双重网格剖分方法
Petrova et al. Space-time evolution of random wave groups with high waves based on the quasi-determinism theory
CN103886129B (zh) 将测井数据离散到储层网格模型的方法和装置
Ikutama et al. Source modeling for predicting ground motions and permanent displacements very close to the fault trace
CN109164488B (zh) 一种梯形网格有限差分地震波场模拟方法
CN109100791B (zh) 基于纵横向空间约束的速度反演方法
CN110187386B (zh) 一种自动快速识别地质构造的dtw地震体属性分析方法
CN106353799A (zh) 一种纵横波联合层析速度反演方法
CN109033588A (zh) 一种基于空间传播的不确定性量化方法
CN110609328B (zh) 一种基于井震结合的含油饱和度预测方法
CN109752758B (zh) 一种地震数据分解方法、系统及存储介质与终端
CN111126439A (zh) 一种煤矿井下二维矿震成像训练数据集生成方法

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
CB03 Change of inventor or designer information

Inventor after: Ding Enjie

Inventor after: Liu Zhongyu

Inventor after: Chen Qing

Inventor after: Liu Yafeng

Inventor before: Chen Qing

Inventor before: Liu Zhongyu

Inventor before: Ding Enjie

Inventor before: Liu Yafeng

CB03 Change of inventor or designer information
GR01 Patent grant
GR01 Patent grant