WO2024078134A1 - Excavation tunnel full-waveform inversion method based on multi-parameter constraint and structure correction - Google Patents

Excavation tunnel full-waveform inversion method based on multi-parameter constraint and structure correction Download PDF

Info

Publication number
WO2024078134A1
WO2024078134A1 PCT/CN2023/113748 CN2023113748W WO2024078134A1 WO 2024078134 A1 WO2024078134 A1 WO 2024078134A1 CN 2023113748 W CN2023113748 W CN 2023113748W WO 2024078134 A1 WO2024078134 A1 WO 2024078134A1
Authority
WO
WIPO (PCT)
Prior art keywords
correction
model
inversion
model update
waveform inversion
Prior art date
Application number
PCT/CN2023/113748
Other languages
French (fr)
Chinese (zh)
Inventor
李圣林
张平松
李洁
郭立全
邱实
Original Assignee
安徽理工大学
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 安徽理工大学 filed Critical 安徽理工大学
Publication of WO2024078134A1 publication Critical patent/WO2024078134A1/en

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/30—Analysis
    • 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
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/62—Physical property of subsurface

Definitions

  • the present invention relates to the field of mining exploration imaging technology, and in particular to a full waveform inversion method for tunnel excavation based on multi-parameter constraints and structural correction, and more specifically to a full waveform inversion method for elastic wave detection in tunnel excavation based on multi-parameter weighted constraints and spatial structural correction of wave velocity profiles.
  • Geological support technology is the basis for the safety of intelligent coal production. It is the basic data source for geological prediction, disturbance perception and risk assessment before, during and after tunneling construction. It is the prerequisite for the implementation of all key technologies of intelligent tunneling.
  • the full waveform inversion method can make full use of the kinematic and dynamic characteristics of seismic waves to obtain underground model parameter information. It has the advantages of high precision in imaging complex structures and good inversion effect of physical parameters. It has achieved good application results in surface seismic exploration and is the best choice for future mine seismic advance detection imaging. It can meet the geological support needs of intelligent tunneling.
  • the seismic advance detection mode is different from surface detection. Its observation system is highly restrictive. It only arranges a series of linearly arranged shot points and detection points on the center line behind the tunneling. For the abnormal body in front of the tunneling, it is similar to single offset detection. In addition, due to the limitation of detection space, the number of shot points and detection points is limited. Therefore, the detection mode has less data and smaller offset, which will lead to the enhancement of the multi-solution of full waveform inversion and the increase of physical parameter calculation error.
  • the present invention provides a full waveform inversion method for tunnel excavation based on multi-parameter constraints and structural correction, which combines multi-parameter weighted constraints and spatial structure correction of velocity profiles to effectively solve the problems of strong limitation of the observation system, small amount of detection data, small offset distance, and strong multi-solution of waveform inversion in tunnel advance detection, improves the inversion effect, and achieves high-precision imaging of abnormal geological structures ahead of excavation.
  • the present invention adopts the following technical solution:
  • a full waveform inversion method for tunneling based on multi-parameter constraints and structural correction includes the following steps:
  • Step 1 Input theoretical simulation data or measured conventional tunnel seismic advance data to construct a full waveform inversion initial model, and then use the multi-scale elastic wave full waveform inversion method to perform a single-scale inversion on the initial model;
  • Step 2 performing multi-parameter weighted constraint structure correction on the single-scale inversion result obtained in step 1 to obtain a preliminary correction result;
  • Step 3 Based on the preset restriction conditions, the one-dimensional velocity profile spatial structure correction and smoothing constraint are performed on the preliminary correction result to obtain the secondary correction result;
  • Step 4 Use the secondary correction result obtained in step 3 as the initial model of the next scale, and continue to perform full waveform inversion in the same manner as step 1;
  • Step 5 Repeat steps 2 to 4 until all scale inversions are completed to obtain the full waveform inversion result of the elastic wave of the tunnel advance detection.
  • step 2 multi-parameter weighted constraint structure correction is performed according to the following formula:
  • ⁇ m is the model update amount of a single iteration
  • ⁇ m′ is the model update amount after multi-parameter weighted structure correction
  • m is the initial model
  • i, j represent the positions of the grid nodes in the z and x directions, respectively
  • mi , jp are the parameter values of the i, j grid node positions of the initial model with a single parameter of longitudinal wave velocity, shear wave velocity or density
  • nz is the number of vertical grid points of the model, nx is the number of transverse grid points of the model
  • Vp is the longitudinal wave velocity;
  • Vs is the shear wave velocity;
  • Den is the density.
  • step 3 the preset restriction conditions are:
  • Constraint condition 1 is based on the value of the lower limit of spatial correction. Suppress and change to 0;
  • the second restriction condition is set according to actual needs.
  • the second restriction condition is based on the range of the structural correction area and the distance between areas.
  • the model update amount in the positive area on both sides of the negative area or in the negative area on both sides of the positive area will be suppressed to 1/5 of the original model update value.
  • step 3 includes:
  • Step 301 Select a vertical grid coordinate y i and extract a one-dimensional wave velocity profile along the lane axis.
  • x is the horizontal axis coordinate
  • Step 302 Calculate the one-dimensional wave velocity profile The corresponding model update amount With slope
  • Step 303 Update the model value below the spatial correction lower limit value Suppress it and turn it into 0 to get the corrected model update amount
  • Step 304 Update the model according to the corrected value According to constraint condition 2, each grid point on the one-dimensional velocity profile is judged and corrected one by one to obtain the corrected model update amount.
  • Step 305 Update the model according to the corrected value With slope According to the change of slope sign, each grid point on the one-dimensional wave velocity profile is judged one by one to obtain the corrected model update amount
  • Step 306 Update the model according to the corrected value Calculate the corrected one-dimensional velocity profile
  • Step 307 take the next vertical axis grid coordinate yi+1 , repeat steps 302 to 306, until the correction of the one-dimensional velocity profiles of all grids is completed, and obtain the secondary correction result of the one-dimensional velocity profile after spatial structure correction and smoothing constraint.
  • step 305 includes:
  • the model update sign of the grid points before and after the two points changes, and there is no point with a slope of 0 between the two points
  • the model update of the grid points in the area between the three adjacent points where the slope sign changes is corrected to the average of the model update of the first slope sign change point and the third slope sign change point, and the corrected model update is obtained.
  • the corrected model update amount Continue to calibrate.
  • the model update signs of all grid points in the range between two non-adjacent points where the slope sign changes do not change
  • the model update of the grid points in the area between the two slope sign change points is corrected to the average of the model update of all grid points in the area, and the corrected model update is obtained.
  • the inversion results in step 5 include: longitudinal wave inversion results, shear wave inversion results and density inversion results.
  • the present invention discloses a full waveform inversion method for tunneling based on multi-parameter constraints and structural correction, which has the following beneficial effects:
  • the present invention adopts an elastic wave full waveform inversion method, which utilizes the longitudinal wave velocity, shear wave velocity and density parameters in the observation record, and can obtain more accurate underground medium information than the acoustic wave full waveform inversion.
  • the present invention performs two corrections on the inversion results at a single scale by performing multi-parameter weighted constraint structure correction and one-dimensional velocity profile spatial structure correction and smoothing constraint, thereby achieving high-precision imaging of abnormal geological structures ahead of tunnel excavation.
  • the imaging results can accurately determine the location and occurrence of geological anomalies, and fill the gap in full waveform inversion technology in the field of mine advance detection.
  • the inversion strategy of the present invention not only provides a more accurate initial model for the next scale inversion, but also avoids the problem of excessive aggravation of false anomalies in multiple iterations. After the false anomalies are eliminated by structural correction, the direction of the full waveform inversion is constrained.
  • the present invention basically solves the problems of small amount of detection data, small offset distance and strong multi-solution of waveform inversion in tunnel advance detection, improves the inversion effect, and improves the accuracy of full waveform inversion by about 20%.
  • the present invention improves the inversion accuracy without increasing the amount of full waveform inversion calculations.
  • FIG1 is a conventional mine earthquake advance detection imaging result provided by the present invention.
  • FIG2 is a flow chart of a full waveform inversion method for detecting elastic waves in tunnels provided by the present invention
  • FIG. 3 is a schematic diagram of determining the update amount and slope change of a one-dimensional wave velocity profile model provided by the present invention
  • FIG4 is a schematic diagram of the structure of the initial model provided by the present invention.
  • FIG5 is a complex advanced geological theoretical model provided by the present invention.
  • FIG6 is an inversion result obtained by using the method of the present invention provided by the present invention.
  • FIG. 7 is a full waveform inversion result of multi-scale elastic waves in conventional time domain provided by the present invention.
  • an embodiment of the present invention discloses a full waveform inversion method for tunneling based on multi-parameter constraints and structural correction, comprising the following steps:
  • Step 1 Input theoretical simulation data or measured conventional tunnel seismic advance detection data to construct a full waveform inversion initial model (as shown in Figure 4), and then use the multi-scale elastic wave full waveform inversion method to perform a single-scale inversion on the initial model; the inversion includes three sets of parameters, namely, the longitudinal wave velocity model, the shear wave velocity model and the density model.
  • Step 2 performing multi-parameter weighted constraint structure correction on the single-scale inversion result obtained in step 1 to obtain a preliminary correction result;
  • Step 3 Based on the preset restriction conditions, the one-dimensional wave velocity profile spatial structure correction and smoothing constraint are performed on the preliminary correction result to obtain the secondary correction result;
  • Step 4 Use the secondary correction result obtained in step 3 as the initial model of the next scale, and continue to perform full waveform inversion in the same manner as step 1;
  • Step 5 Repeat steps 2 to 4 until all scale inversions are completed, and obtain the full waveform inversion results of the tunnel advance detection elastic wave.
  • the inversion results include: longitudinal wave inversion results, shear wave inversion results and density inversion results.
  • the seismic source is distributed in the middle of the coal seam, and multiple horizontal component and vertical component detectors are arranged linearly with the seismic source and also distributed in the middle of the coal seam to receive seismic signals.
  • the inversion initial model is set according to the uniform layered coal formation medium.
  • the conventional time domain multi-scale elastic wave full waveform inversion method Due to the special observation system conditions of tunnel advance detection, the conventional time domain multi-scale elastic wave full waveform inversion method has increased multi-solution problems, and it is difficult to further improve the inversion accuracy from the data perspective.
  • the present invention discovers the special structural features implied by the geological anomaly in the full waveform inversion results of tunnel advance detection, and constructs correction terms based on the structural features to constrain the direction of the full waveform inversion, and obtains the inversion results that meet the preset structural features, thereby improving the inversion effect.
  • the lateral velocity-depth curve of the vertical grid points inside the larger structure has a "wavy line” feature, and the model update sign is consistent within the range of the "wavy line”.
  • the present invention introduces multi-parameter weighted constraints and one-dimensional velocity profile spatial structure correction to perform multiple corrections on the inversion results, which can effectively solve the problems of small amount of detection data, small offset distance, and strong multi-solution of waveform inversion in tunnel advance detection, improve the inversion effect, and achieve high-precision imaging of abnormal geological structures ahead of excavation.
  • step 2 the correction method of step 2 and step 3 is described in detail.
  • step 2 multi-parameter weighted constraint structure correction is performed according to the following formula:
  • ⁇ m is the model update amount of a single iteration (i.e., the change between this iteration and the previous iteration model); ⁇ m′ is the model update amount after multi-parameter weighted structure correction;
  • m is the initial model;
  • i, j represent the positions of the grid nodes in the z and x directions, respectively;
  • mi , jp are the parameter values of the i, j grid node positions of the initial model with a single parameter of longitudinal wave velocity, shear wave velocity or density;
  • nz is the number of vertical grid points of the model, and nx is the number of horizontal grid points of the model;
  • Vp is the longitudinal wave velocity;
  • Vs is the shear wave velocity;
  • Den is the density
  • step 3 the preset restriction condition is:
  • Constraint 1 is based on the value of the lower limit of spatial correction.
  • the model update amount below the value of the lower limit of spatial correction is Suppress and become 0; in actual application, it is set independently according to the data processing effect, generally determined by experiments before formal inversion. For example, the value is set to be less than 5% of the initial model parameter.
  • Constraint 2 is based on the range of the structural correction area and the distance between areas. When there is a positive area on both sides of a negative area where the model update amount exists, or when there is a negative area on both sides of a positive area where the model update amount exists, and when the area range and the distance between areas meet Constraint 2, the model update amount in the positive area on both sides of the negative area or in the negative area on both sides of the positive area will be suppressed to 1/5 of the original model update value.
  • distance between areas refers to: when “there is a positive update area on both sides of the negative area”, the distance between the one-sided “positive area” and the middle “negative area”; or: when “there is a negative update area on both sides of the positive update area”, the distance between the one-sided "negative area” and the middle "positive area”.
  • Regular range refers to the range size of the "positive value region” on one side when “there is a region with positive update amount on both sides of the negative value region", or the range size of the "negative value region” on one side when “there is a region with negative update amount on both sides of the positive value region”.
  • the range of regions and the distance between regions are generally determined by experiments before formal inversion. For example, when the range of regions is: "There is a positive update region on both sides of the negative region", the range of the "positive region” on one side cannot be less than 1/2 of the range of the middle "positive region".
  • Distance between regions When there is a positive region on both sides of a negative region, the distance between the positive region on one side and the negative region in the middle cannot be greater than the size of the negative region in the middle.
  • step 3 includes:
  • Step 301 Select a vertical grid coordinate y i and extract a one-dimensional wave velocity profile along the lane axis.
  • x is the horizontal axis coordinate
  • Step 302 Calculate the one-dimensional wave velocity profile The corresponding model update amount With slope
  • the model update amount is: is the one-dimensional wave velocity profile at the e-th iteration, is the one-dimensional wave velocity profile at the e-1th iteration.
  • Slope is the one-dimensional wave velocity profile
  • the value of the j-th grid point in the horizontal direction is the one-dimensional wave velocity profile
  • the value of the j-1th grid point in the horizontal direction is the horizontal grid spacing of the model.
  • Step 303 Update the model value below the spatial correction lower limit value Suppress it and turn it into 0 to get the corrected model update amount
  • Step 304 Update the model according to the corrected value According to constraint condition 2, each grid point on the one-dimensional velocity profile is judged and corrected one by one to obtain the corrected model update amount.
  • Step 305 Update the model according to the corrected value With slope According to the change of slope sign, each grid point on the one-dimensional wave velocity profile is judged one by one to obtain the corrected model update amount
  • the corrected model update amount Continue to calibrate.
  • the model update signs of all grid points in the range between two non-adjacent points where the slope sign changes do not change
  • the model update of the grid points in the area between the two slope sign change points is corrected to the average of the model update of all grid points in the area, and the corrected model update is obtained.
  • search and determine according to the set area range limit means searching in order with 10 grid points as a round.
  • Step 306 Update the model according to the corrected value Calculate the corrected one-dimensional velocity profile
  • Step 307 take the next vertical axis grid coordinate yi+1 , repeat steps 302 to 306, until the correction of the one-dimensional velocity profiles of all grids is completed, and obtain the secondary correction result of the one-dimensional velocity profile after spatial structure correction and smoothing constraint.
  • the conventional time domain multi-scale elastic wave full waveform inversion results are used as comparative research objects.
  • FIG. 6 By comparing FIG. 6 with FIG. 7, it can be seen that the full waveform inversion results obtained by the present invention are closer to the real model.
  • the inversion results the boundaries on both sides of the collapse column and the fault fracture zone structure are clear and the internal parameters of the structure are well restored, and the small-scale fault lithology interface is also restored and revealed; the most obvious thing is that the false anomalies and disturbance interference in the results are basically completely suppressed, and the inversion accuracy has been greatly improved. Comparing the numerical difference between the inversion results and the real model, the full waveform inversion accuracy has been improved by about 20%.
  • each embodiment is described in a progressive manner, and each embodiment focuses on the differences from other embodiments.
  • the same or similar parts between the embodiments can be referred to each other.
  • the description is relatively simple, and the relevant parts can be referred to the method part.

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)
  • Radar Systems Or Details Thereof (AREA)

Abstract

Disclosed in the present invention is an excavation tunnel full-waveform inversion method based on multi-parameter constraint and structure correction. The method comprises the following steps: step 1, constructing a full-waveform inversion initial model, and then performing single-scale inversion on the initial model by using a multi-scale elastic-wave full-waveform inversion method; step 2, performing multi-parameter weighted constraint and structure correction on the single-scale inversion result, so as to obtain a preliminary correction result; step 3, on the basis of a preset limiting condition, performing one-dimensional wave-velocity profile space structure correction and smooth constraint on the preliminary correction result, so as to obtain a secondary correction result; step 4, using the secondary correction result as an initial model of the next scale, and continuing to perform full-waveform inversion; and repeating steps 2-4 until the inversion of all the scales is completed, so as to obtain a full-waveform inversion result. The present invention can effectively solve the problems of the observation system being highly restricted, the detection data volume being relatively small, the offset distance being relatively small, and the multiplicity of solution of waveform inversion being high in terms of advanced detection of an excavation tunnel, thus improving the inversion effect.

Description

基于多参数约束和结构校正的随掘巷道全波形反演方法Full waveform inversion method for tunneling based on multi-parameter constraints and structural correction 技术领域Technical Field
本发明涉及矿山勘探成像技术领域，尤其涉及一种基于多参数约束和结构校正的随掘巷道全波形反演方法，更具体的说是涉及一种基于多参数加权约束与波速剖面空间结构校正的随掘巷道探测弹性波全波形反演方法。The present invention relates to the field of mining exploration imaging technology, and in particular to a full waveform inversion method for tunnel excavation based on multi-parameter constraints and structural correction, and more specifically to a full waveform inversion method for elastic wave detection in tunnel excavation based on multi-parameter weighted constraints and spatial structural correction of wave velocity profiles.
背景技术Background technique
作为煤矿生产的两大核心环节之一，掘进智能化发展需求极为迫切，但巷道掘进过程中，煤与瓦斯突出、突水等灾害严重威胁着掘进生产安全和矿工的人身安全。地质保障技术是煤炭智能化生产安全保障的基础，是实现巷道掘进施工前、中、后地质预判、扰动感知与风险评估的基础数据来源，是一切智能掘进关键技术实施的前提保障。As one of the two core links of coal mine production, the demand for intelligent tunneling is extremely urgent. However, during tunneling, disasters such as coal and gas outbursts and water inrush seriously threaten tunneling production safety and the personal safety of miners. Geological support technology is the basis for the safety of intelligent coal production. It is the basic data source for geological prediction, disturbance perception and risk assessment before, during and after tunneling construction. It is the prerequisite for the implementation of all key technologies of intelligent tunneling.
但是，目前矿山地震超前探测成像解释结果多为“画弧”，假异常界面较多，成像精度较低(如图1所示)，难以满足巷道智能掘进地质保障需求。However, the current interpretation results of mine seismic advance detection imaging are mostly "arc drawing", with many false abnormal interfaces and low imaging accuracy (as shown in Figure 1), which makes it difficult to meet the geological support needs of intelligent tunneling.
全波形反演方法可以充分利用地震波的运动学和动力学特征来获取地下模型参数信息，具有复杂构造成像精度高、物性参数反演效果好等优点，在地表地震勘探中取得了很好的应用效果，是未来矿山地震超前探测成像的最优选择，可以满足巷道智能掘进地质保障需求。但地震超前探测模式与地表探测有所不同，其观测系统限制性较强，仅在掘进迎头后方中心线上布设一系列成线性排列的炮点与检波点，对于掘进前方异常体近似于单偏移距探测，并且，受探测空间限制，布设的炮点与检波点数量有限，因此，该探测模式数据量较少且偏移距偏小，会导致全波形反演的多解性增强，物性参数计算误差增大。The full waveform inversion method can make full use of the kinematic and dynamic characteristics of seismic waves to obtain underground model parameter information. It has the advantages of high precision in imaging complex structures and good inversion effect of physical parameters. It has achieved good application results in surface seismic exploration and is the best choice for future mine seismic advance detection imaging. It can meet the geological support needs of intelligent tunneling. However, the seismic advance detection mode is different from surface detection. Its observation system is highly restrictive. It only arranges a series of linearly arranged shot points and detection points on the center line behind the tunneling. For the abnormal body in front of the tunneling, it is similar to single offset detection. In addition, due to the limitation of detection space, the number of shot points and detection points is limited. Therefore, the detection mode has less data and smaller offset, which will lead to the enhancement of the multi-solution of full waveform inversion and the increase of physical parameter calculation error.
因此，如何提供一种适用于矿山随掘巷道探测模式下的高精度全波形反演方法是本领域技术人员亟需解决的问题。 Therefore, how to provide a high-precision full-waveform inversion method suitable for mine tunnel detection mode is an urgent problem to be solved by technical personnel in this field.
发明内容Summary of the invention
有鉴于此，本发明提供了一种基于多参数约束和结构校正的随掘巷道全波形反演方法，结合多参数加权约束和波速剖面空间结构校正，有效解决巷道超前探测中观测系统限制性较强、探测数据量较少且偏移距偏小、波形反演多解性强的问题，提升反演效果，实现掘进前方异常地质构造高精度成像。In view of this, the present invention provides a full waveform inversion method for tunnel excavation based on multi-parameter constraints and structural correction, which combines multi-parameter weighted constraints and spatial structure correction of velocity profiles to effectively solve the problems of strong limitation of the observation system, small amount of detection data, small offset distance, and strong multi-solution of waveform inversion in tunnel advance detection, improves the inversion effect, and achieves high-precision imaging of abnormal geological structures ahead of excavation.
为了实现上述目的，本发明采用如下技术方案：In order to achieve the above object, the present invention adopts the following technical solution:
一种基于多参数约束和结构校正的随掘巷道全波形反演方法，包括以下步骤：A full waveform inversion method for tunneling based on multi-parameter constraints and structural correction includes the following steps:
步骤1、输入理论模拟数据或实测常规巷道地震超前数据，构建全波形反演初始模型，再采用多尺度弹性波全波形反演方法对初始模型进行单一尺度反演；Step 1: Input theoretical simulation data or measured conventional tunnel seismic advance data to construct a full waveform inversion initial model, and then use the multi-scale elastic wave full waveform inversion method to perform a single-scale inversion on the initial model;
步骤2、对步骤1得到的单一尺度反演结果进行多参数加权约束结构校正，得到初步校正结果；Step 2, performing multi-parameter weighted constraint structure correction on the single-scale inversion result obtained in step 1 to obtain a preliminary correction result;
步骤3、基于预先设定的限制条件，对初步校正结果进行一维波速剖面空间结构校正和平滑约束，得到二次校正结果；Step 3: Based on the preset restriction conditions, the one-dimensional velocity profile spatial structure correction and smoothing constraint are performed on the preliminary correction result to obtain the secondary correction result;
步骤4、将步骤3得到的二次校正结果作为下一尺度的初始模型，继续按照步骤1的方式进行全波形反演；Step 4: Use the secondary correction result obtained in step 3 as the initial model of the next scale, and continue to perform full waveform inversion in the same manner as step 1;
步骤5、重复执行步骤2～步骤4，直至所有尺度反演完成，得到巷道超前探测弹性波全波形反演结果。Step 5: Repeat steps 2 to 4 until all scale inversions are completed to obtain the full waveform inversion result of the elastic wave of the tunnel advance detection.
进一步的，步骤2中，按照下式进行多参数加权约束结构校正：
Furthermore, in step 2, multi-parameter weighted constraint structure correction is performed according to the following formula:
其中，Δm为单次迭代模型更新量；Δm′为多参数加权结构校正后的模型更新量；m为初始模型；i,j分别表示z,x方向上网格节点的位置；mi,j p为纵波速度、横波速度或密度单一参数初始模型i,j网格节点位置的参数数值；nz为模型竖向网格点数，nx为模型横向网格点数；Vp为纵波速度；Vs为横波速度；Den为密度。 Among them, Δm is the model update amount of a single iteration; Δm′ is the model update amount after multi-parameter weighted structure correction; m is the initial model; i, j represent the positions of the grid nodes in the z and x directions, respectively; mi , jp are the parameter values of the i, j grid node positions of the initial model with a single parameter of longitudinal wave velocity, shear wave velocity or density; nz is the number of vertical grid points of the model, nx is the number of transverse grid points of the model; Vp is the longitudinal wave velocity; Vs is the shear wave velocity; Den is the density.
进一步的，步骤3中，预先设定的限制条件为：Furthermore, in step 3, the preset restriction conditions are:
根据实际需要设置限制条件一，限制条件一以空间校正下限value值为判定依据，将低于空间校正下限value值的模型更新量进行压制，变为0；Set constraint condition 1 according to actual needs. Constraint condition 1 is based on the value of the lower limit of spatial correction. Suppress and change to 0;
根据实际需要设置限制条件二，限制条件二以结构校正区域范围以及区域之间距离为判定依据，当模型更新量为负值区域两侧均存在一更新量为正值区域、或者当模型更新量为正值区域两侧均存在一更新量为负值区域，当区域范围以及区域之间距离满足限制条件二时，将负值区域两侧存在的正值区域内或者正值区域两侧存在的负值区域内的模型更新量进行压制，变为原来模型更新数值的1/5。The second restriction condition is set according to actual needs. The second restriction condition is based on the range of the structural correction area and the distance between areas. When there is a positive area on both sides of the negative area where the model update amount exists, or when there is a negative area on both sides of the positive area where the model update amount exists, when the area range and the distance between areas meet the second restriction condition, the model update amount in the positive area on both sides of the negative area or in the negative area on both sides of the positive area will be suppressed to 1/5 of the original model update value.
进一步的，步骤3包括：Furthermore, step 3 includes:
步骤301、选取一个纵轴网格坐标yi，沿巷道轴线提取一维波速剖面其中，x为横轴坐标；Step 301: Select a vertical grid coordinate y i and extract a one-dimensional wave velocity profile along the lane axis. Among them, x is the horizontal axis coordinate;
步骤302、计算一维波速剖面对应的模型更新量与斜率 Step 302: Calculate the one-dimensional wave velocity profile The corresponding model update amount With slope
步骤303、将低于空间校正下限value值的模型更新量进行压制，变为0，得到校正后的模型更新量 Step 303: Update the model value below the spatial correction lower limit value Suppress it and turn it into 0 to get the corrected model update amount
步骤304、根据校正后的模型更新量按照限制条件二对一维波速剖面上的每个网格点进行逐一判定和校正，得到校正后的模型更新量 Step 304: Update the model according to the corrected value According to constraint condition 2, each grid point on the one-dimensional velocity profile is judged and corrected one by one to obtain the corrected model update amount.
步骤305、根据校正后的模型更新量与斜率按照斜率符号变化情况对一维波速剖面上的每个网格点进行逐一判定，得到校正后的模型更新量 Step 305: Update the model according to the corrected value With slope According to the change of slope sign, each grid point on the one-dimensional wave velocity profile is judged one by one to obtain the corrected model update amount
步骤306、根据校正后的模型更新量计算校正后的一维波速剖面 Step 306: Update the model according to the corrected value Calculate the corrected one-dimensional velocity profile
步骤307、取下一个纵轴网格坐标yi+1，重复步骤302～步骤306，直至完成所有网格一维波速剖面的校正，得到一维波速剖面空间结构校正与平滑约束后的二次校正结果。Step 307, take the next vertical axis grid coordinate yi+1 , repeat steps 302 to 306, until the correction of the one-dimensional velocity profiles of all grids is completed, and obtain the secondary correction result of the one-dimensional velocity profile after spatial structure correction and smoothing constraint.
进一步的，步骤305包括： Further, step 305 includes:
当斜率符号发生变化的相邻三个点之间存在两个模型更新量为零的点，该两个点前后网格点模型更新量符号发生变化，且该两个点之间没有斜率为0的点，则对斜率符号发生变化的相邻三个点之间区域内的网格点的模型更新量进行校正，变为第一个斜率符号变化点与第三个斜率符号变化点两点模型更新量的平均值，得到校正后的模型更新量 When there are two points with zero model update between the three adjacent points where the slope sign changes, the model update sign of the grid points before and after the two points changes, and there is no point with a slope of 0 between the two points, the model update of the grid points in the area between the three adjacent points where the slope sign changes is corrected to the average of the model update of the first slope sign change point and the third slope sign change point, and the corrected model update is obtained.
再根据校正后的模型更新量继续校正，当斜率符号发生变化的不相邻的两个点之间范围内所有网格点模型更新量符号未发生改变，则对两个斜率符号变化点之间区域内的网格点的模型更新量进行校正，变为区域内所有网格点模型更新量的平均值，得到校正后的模型更新量 Then according to the corrected model update amount Continue to calibrate. When the model update signs of all grid points in the range between two non-adjacent points where the slope sign changes do not change, the model update of the grid points in the area between the two slope sign change points is corrected to the average of the model update of all grid points in the area, and the corrected model update is obtained.
进一步的，步骤5中的反演结果包括：纵波反演结果、横波反演结果和密度反演结果。Furthermore, the inversion results in step 5 include: longitudinal wave inversion results, shear wave inversion results and density inversion results.
经由上述的技术方案可知，与现有技术相比，本发明公开提供了一种基于多参数约束和结构校正的随掘巷道全波形反演方法，具有以下有益效果：It can be seen from the above technical solutions that, compared with the prior art, the present invention discloses a full waveform inversion method for tunneling based on multi-parameter constraints and structural correction, which has the following beneficial effects:
(1)本发明采用的是弹性波全波形反演方法，利用的是观测记录中的纵波速度、横波速度及密度参数，可得到比声波全波形反演更为准确的地下介质信息。(1) The present invention adopts an elastic wave full waveform inversion method, which utilizes the longitudinal wave velocity, shear wave velocity and density parameters in the observation record, and can obtain more accurate underground medium information than the acoustic wave full waveform inversion.
(2)本发明通过对单一尺度下的反演结果进行多参数加权约束结构校正和一维波速剖面空间结构校正和平滑约束，进行了两次校正，实现了巷道掘进前方异常地质构造的高精度成像，成像结果可以准确判断地质异常体的位置和产状，并填补了矿山超前探测领域全波形反演技术空白。(2) The present invention performs two corrections on the inversion results at a single scale by performing multi-parameter weighted constraint structure correction and one-dimensional velocity profile spatial structure correction and smoothing constraint, thereby achieving high-precision imaging of abnormal geological structures ahead of tunnel excavation. The imaging results can accurately determine the location and occurrence of geological anomalies, and fill the gap in full waveform inversion technology in the field of mine advance detection.
(3)本发明反演策略不仅为下一尺度反演提供了精度更高的初始模型，同时也避免了多次迭代假异常过度加重的问题，假异常被结构校正消除后，约束了全波形反演的方向。(3) The inversion strategy of the present invention not only provides a more accurate initial model for the next scale inversion, but also avoids the problem of excessive aggravation of false anomalies in multiple iterations. After the false anomalies are eliminated by structural correction, the direction of the full waveform inversion is constrained.
(4)本发明基本解决了巷道超前探测中探测数据量较少且偏移距偏小、波形反演多解性强的问题，提升了反演效果，全波形反演精度提升了20％左右。(4) The present invention basically solves the problems of small amount of detection data, small offset distance and strong multi-solution of waveform inversion in tunnel advance detection, improves the inversion effect, and improves the accuracy of full waveform inversion by about 20%.
(5)本发明在提高反演精度的同时，不会造成全波形反演计算量的增加。 (5) The present invention improves the inversion accuracy without increasing the amount of full waveform inversion calculations.
附图说明BRIEF DESCRIPTION OF THE DRAWINGS
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图获得其他的附图。In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on the provided drawings without paying creative work.
图1为本发明提供的常规矿山地震超前探测成像结果；FIG1 is a conventional mine earthquake advance detection imaging result provided by the present invention;
图2为本发明提供的随掘巷道探测弹性波全波形反演方法的流程图；FIG2 is a flow chart of a full waveform inversion method for detecting elastic waves in tunnels provided by the present invention;
图3为本发明提供的一维波速剖面模型更新量和斜率变化判定示意图；3 is a schematic diagram of determining the update amount and slope change of a one-dimensional wave velocity profile model provided by the present invention;
图4为本发明提供的初始模型的结构示意图；FIG4 is a schematic diagram of the structure of the initial model provided by the present invention;
图5为本发明提供的复杂超前地质理论模型；FIG5 is a complex advanced geological theoretical model provided by the present invention;
图6为本发明提供的采用本发明方法得到的反演结果；FIG6 is an inversion result obtained by using the method of the present invention provided by the present invention;
图7为本发明提供的采用常规时间域多尺度弹性波全波形反演结果。FIG. 7 is a full waveform inversion result of multi-scale elastic waves in conventional time domain provided by the present invention.
具体实施方式Detailed ways
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.
如图2所示，本发明实施例公开了一种基于多参数约束和结构校正的随掘巷道全波形反演方法，包括以下步骤：As shown in FIG2 , an embodiment of the present invention discloses a full waveform inversion method for tunneling based on multi-parameter constraints and structural correction, comprising the following steps:
步骤1、输入理论模拟数据或实测常规巷道地震超前探测数据，构建全波形反演初始模型(如图4所示)，再采用多尺度弹性波全波形反演方法对初始模型进行单一尺度反演；反演包括三组参数，即纵波速度模型、横波速度模型和密度模型。Step 1: Input theoretical simulation data or measured conventional tunnel seismic advance detection data to construct a full waveform inversion initial model (as shown in Figure 4), and then use the multi-scale elastic wave full waveform inversion method to perform a single-scale inversion on the initial model; the inversion includes three sets of parameters, namely, the longitudinal wave velocity model, the shear wave velocity model and the density model.
步骤2、对步骤1得到的单一尺度反演结果进行多参数加权约束结构校正，得到初步校正结果；Step 2, performing multi-parameter weighted constraint structure correction on the single-scale inversion result obtained in step 1 to obtain a preliminary correction result;
步骤3、基于预先设定的限制条件，对初步校正结果进行一维波速剖面空间结构校正和平滑约束，得到二次校正结果； Step 3: Based on the preset restriction conditions, the one-dimensional wave velocity profile spatial structure correction and smoothing constraint are performed on the preliminary correction result to obtain the secondary correction result;
步骤4、将步骤3得到的二次校正结果作为下一尺度的初始模型，继续按照步骤1的方式进行全波形反演；Step 4: Use the secondary correction result obtained in step 3 as the initial model of the next scale, and continue to perform full waveform inversion in the same manner as step 1;
步骤5、重复执行步骤2～步骤4，直至所有尺度反演完成，得到巷道超前探测弹性波全波形反演结果。反演结果包括：纵波反演结果、横波反演结果和密度反演结果。Step 5: Repeat steps 2 to 4 until all scale inversions are completed, and obtain the full waveform inversion results of the tunnel advance detection elastic wave. The inversion results include: longitudinal wave inversion results, shear wave inversion results and density inversion results.
本发明中，震源分布在煤层中部，多个水平分量和垂直分量检波器与震源成线性排列同样分布在煤层中部接收地震信号，按照均匀层状煤系地层介质情况设置反演初始模型。In the present invention, the seismic source is distributed in the middle of the coal seam, and multiple horizontal component and vertical component detectors are arranged linearly with the seismic source and also distributed in the middle of the coal seam to receive seismic signals. The inversion initial model is set according to the uniform layered coal formation medium.
由于巷道超前探测特殊的观测系统条件，常规时间域多尺度弹性波全波形反演方法多解性问题增强，很难从数据角度进一步改善反演精度。本发明从模型结构出发，发现了巷道超前探测全波形反演结果中地质异常体暗含的特殊结构特征，并根据结构特征构建校正项，约束全波形反演的方向，得到了满足预设结构特征的反演结果，提高了反演效果。Due to the special observation system conditions of tunnel advance detection, the conventional time domain multi-scale elastic wave full waveform inversion method has increased multi-solution problems, and it is difficult to further improve the inversion accuracy from the data perspective. Starting from the model structure, the present invention discovers the special structural features implied by the geological anomaly in the full waveform inversion results of tunnel advance detection, and constructs correction terms based on the structural features to constrain the direction of the full waveform inversion, and obtains the inversion results that meet the preset structural features, thereby improving the inversion effect.
发现的巷道超前探测弹性波全波形反演结果中地质异常体暗含的特殊结构特征为：The special structural features implied by the geological anomaly bodies in the full waveform inversion results of elastic waves found in the tunnel advance detection are:
(1)对于异常地质构造的恢复，密度参数的反演效果最好，而纵横波参数中假异常与扰动干扰压制效果更好；(1) For the restoration of abnormal geological structures, the inversion effect of density parameters is the best, while the suppression of false anomalies and disturbances in P- and S-wave parameters is better;
(2)纵、横波以及密度虽然在物性参数上具有较大的差异，但它们具有相同的地质结构特征；(2) Although P-wave, S-wave and density have great differences in physical parameters, they have the same geological structure characteristics;
(3)异常地质构造两侧边界附近均存在与自身属性特征相反的假异常区域；(3) There are false anomaly areas with opposite attribute characteristics to the abnormal geological structure near both sides of the boundary;
(4)垂直方向网格点横向速度-深度曲线中存在“W”或“M”型异常特征，且该“W”或“M”区域中间范围模型更新符号与两侧范围模型更新符号相反；(4) There are “W” or “M” type anomaly features in the lateral velocity-depth curve of the vertical grid points, and the model update sign in the middle range of the “W” or “M” area is opposite to the model update sign on both sides;
(5)体积较大构造内部垂直方向网格点横向速度-深度曲线存在“波浪线”特征，且该“波浪线”涉及范围内模型更新符号一致。(5) The lateral velocity-depth curve of the vertical grid points inside the larger structure has a "wavy line" feature, and the model update sign is consistent within the range of the "wavy line".
基于发现的上述特征，本发明引入了多参数加权约束和一维波速剖面空间结构校正对反演结果进行多次校正，能够有效解决巷道超前探测中探测数据量较少且偏移距偏小、波形反演多解性强的问题，提升反演效果，实现掘进前方异常地质构造高精度成像。Based on the above-mentioned characteristics discovered, the present invention introduces multi-parameter weighted constraints and one-dimensional velocity profile spatial structure correction to perform multiple corrections on the inversion results, which can effectively solve the problems of small amount of detection data, small offset distance, and strong multi-solution of waveform inversion in tunnel advance detection, improve the inversion effect, and achieve high-precision imaging of abnormal geological structures ahead of excavation.
下面，对步骤2和步骤3的校正方式进行详细说明。Next, the correction method of step 2 and step 3 is described in detail.
在一个具体实施例中，步骤2中，按照下式进行多参数加权约束结构校正：
In a specific embodiment, in step 2, multi-parameter weighted constraint structure correction is performed according to the following formula:
其中，Δm为单次迭代模型更新量(即本次迭代与上一次迭代模型之间的变化量)；Δm′为多参数加权结构校正后的模型更新量；m为初始模型；i,j分别表示z,x方向上网格节点的位置；mi,j p为纵波速度、横波速度或密度单一参数初始模型i,j网格节点位置的参数数值；nz为模型竖向网格点数，nx为模型横向网格点数；Vp为纵波速度；Vs为横波速度；Den为密度Wherein, Δm is the model update amount of a single iteration (i.e., the change between this iteration and the previous iteration model); Δm′ is the model update amount after multi-parameter weighted structure correction; m is the initial model; i, j represent the positions of the grid nodes in the z and x directions, respectively; mi , jp are the parameter values of the i, j grid node positions of the initial model with a single parameter of longitudinal wave velocity, shear wave velocity or density; nz is the number of vertical grid points of the model, and nx is the number of horizontal grid points of the model; Vp is the longitudinal wave velocity; Vs is the shear wave velocity; Den is the density
在一个具体实施例中，步骤3中，预先设定的限制条件为：In a specific embodiment, in step 3, the preset restriction condition is:
限制条件一，限制条件一以空间校正下限value值为判定依据，将低于空间校正下限value值的模型更新量进行压制，变为0；实际应用时，根据数据处理效果进行自主设定，一般通过正式反演前的实验确定。例如，将该值设定为小于初始模型参数的5％。Constraint 1: Constraint 1 is based on the value of the lower limit of spatial correction. The model update amount below the value of the lower limit of spatial correction is Suppress and become 0; in actual application, it is set independently according to the data processing effect, generally determined by experiments before formal inversion. For example, the value is set to be less than 5% of the initial model parameter.
限制条件二，限制条件二以结构校正区域范围以及区域之间距离为判定依据，当模型更新量为负值区域两侧均存在一更新量为正值区域、或者当模型更新量为正值区域两侧均存在一更新量为负值区域，当区域范围以及区域之间距离满足限制条件二时，将负值区域两侧存在的正值区域内或者正值区域两侧存在的负值区域内的模型更新量进行压制，变为原来模型更新数值的1/5。Constraint 2: Constraint 2 is based on the range of the structural correction area and the distance between areas. When there is a positive area on both sides of a negative area where the model update amount exists, or when there is a negative area on both sides of a positive area where the model update amount exists, and when the area range and the distance between areas meet Constraint 2, the model update amount in the positive area on both sides of the negative area or in the negative area on both sides of the positive area will be suppressed to 1/5 of the original model update value.
其中，“区域之间距离”是指：当“负值区域两侧均存在一更新量为正值区域”时，单侧“正值区域”距离中间“负值区域”的距离；或者是：“更新量为正值区域两侧均存在一更新量为负值区域”时，单侧“负值区域”距离中间“正值区域”的距离。 Among them, "distance between areas" refers to: when "there is a positive update area on both sides of the negative area", the distance between the one-sided "positive area" and the middle "negative area"; or: when "there is a negative update area on both sides of the positive update area", the distance between the one-sided "negative area" and the middle "positive area".
“区域范围”是指：“负值区域两侧均存在一更新量为正值区域”时，单侧“正值区域”的范围大小，或者是：“更新量为正值区域两侧均存在一更新量为负值区域”时，单侧“负值区域”的范围大小。"Region range" refers to the range size of the "positive value region" on one side when "there is a region with positive update amount on both sides of the negative value region", or the range size of the "negative value region" on one side when "there is a region with negative update amount on both sides of the positive value region".
区域范围以及区域之间距离一般通过正式反演前的实验确定。例如：“区域范围：“负值区域两侧均存在一更新量为正值区域”时，单侧“正值区域”的范围大小不能小于中间“正值区域”范围的1/2。The range of regions and the distance between regions are generally determined by experiments before formal inversion. For example, when the range of regions is: "There is a positive update region on both sides of the negative region", the range of the "positive region" on one side cannot be less than 1/2 of the range of the middle "positive region".
“区域之间距离”：当“负值区域两侧均存在一更新量为正值区域”时，单侧“正值区域”距离中间“负值区域”的距离不能大于中间“负值区域”的大小。"Distance between regions": When there is a positive region on both sides of a negative region, the distance between the positive region on one side and the negative region in the middle cannot be greater than the size of the negative region in the middle.
具体的，步骤3包括：Specifically, step 3 includes:
步骤301、选取一个纵轴网格坐标yi，沿巷道轴线提取一维波速剖面其中，x为横轴坐标；Step 301: Select a vertical grid coordinate y i and extract a one-dimensional wave velocity profile along the lane axis. Among them, x is the horizontal axis coordinate;
步骤302、计算一维波速剖面对应的模型更新量与斜率 Step 302: Calculate the one-dimensional wave velocity profile The corresponding model update amount With slope
其中，模型更新量：为第e次迭代时的一维波速剖面，为第e-1次迭代时的一维波速剖面。Among them, the model update amount is: is the one-dimensional wave velocity profile at the e-th iteration, is the one-dimensional wave velocity profile at the e-1th iteration.
斜率：为一维波速剖面横向第j个网格点的数值，为一维波速剖面横向第j-1个网格点的数值，dx为模型横向网格间距。Slope: is the one-dimensional wave velocity profile The value of the j-th grid point in the horizontal direction, is the one-dimensional wave velocity profile The value of the j-1th grid point in the horizontal direction, dx is the horizontal grid spacing of the model.
步骤303、将低于空间校正下限value值的模型更新量进行压制，变为0，得到校正后的模型更新量 Step 303: Update the model value below the spatial correction lower limit value Suppress it and turn it into 0 to get the corrected model update amount
步骤304、根据校正后的模型更新量按照限制条件二对一维波速剖面上的每个网格点进行逐一判定和校正，得到校正后的模型更新量 Step 304: Update the model according to the corrected value According to constraint condition 2, each grid point on the one-dimensional velocity profile is judged and corrected one by one to obtain the corrected model update amount.
步骤305、根据校正后的模型更新量与斜率按照斜率符号变化情况对一维波速剖面上的每个网格点进行逐一判定，得到校正后的模型更新量 Step 305: Update the model according to the corrected value With slope According to the change of slope sign, each grid point on the one-dimensional wave velocity profile is judged one by one to obtain the corrected model update amount
当斜率符号发生变化的相邻三个点之间存在两个模型更新量为零的点(如图3所示)，该两个点前后网格点模型更新量符号发生变化，且该两个点之间没有斜率为0的点，则对斜率符号发生变化的相邻三个点之间区域内的网格点的模型更新量进行校正，变为第一个斜率符号变化点与第三个斜率符号变化点两点模型更新量的平均值，得到校正后的模型更新量 When there are two points with zero model update amount between the three adjacent points where the slope sign changes (as shown in Figure 3), the model update amount of the grid points before and after the two points changes in sign, and there is no point with a slope of 0 between the two points, then the model update amount of the grid points in the area between the three adjacent points where the slope sign changes is corrected to the average of the model update amounts of the first slope sign change point and the third slope sign change point, and the corrected model update amount is obtained
再根据校正后的模型更新量继续校正，当斜率符号发生变化的不相邻的两个点之间范围内所有网格点模型更新量符号未发生改变，则对两个斜率符号变化点之间区域内的网格点的模型更新量进行校正，变为区域内所有网格点模型更新量的平均值，得到校正后的模型更新量 Then according to the corrected model update amount Continue to calibrate. When the model update signs of all grid points in the range between two non-adjacent points where the slope sign changes do not change, the model update of the grid points in the area between the two slope sign change points is corrected to the average of the model update of all grid points in the area, and the corrected model update is obtained.
为控制校正范围，需按照设定的区域范围限值进行搜索判定。也就是说控制最大校正范围，比如设置为10个网格点，“按照设定的区域范围限值进行搜索判定”就是按照10个网格点为一个轮次依次搜索。To control the correction range, it is necessary to search and determine according to the set area range limit. In other words, to control the maximum correction range, for example, if it is set to 10 grid points, "search and determine according to the set area range limit" means searching in order with 10 grid points as a round.
步骤306、根据校正后的模型更新量计算校正后的一维波速剖面 Step 306: Update the model according to the corrected value Calculate the corrected one-dimensional velocity profile
步骤307、取下一个纵轴网格坐标yi+1，重复步骤302～步骤306，直至完成所有网格一维波速剖面的校正，得到一维波速剖面空间结构校正与平滑约束后的二次校正结果。Step 307, take the next vertical axis grid coordinate yi+1 , repeat steps 302 to 306, until the correction of the one-dimensional velocity profiles of all grids is completed, and obtain the secondary correction result of the one-dimensional velocity profile after spatial structure correction and smoothing constraint.
为了进一步验证本发明基于多参数加权约束与波速剖面空间结构校正的随掘巷道探测弹性波全波形反演方法的有效性与高效性，将本发明提出的方案用于复杂超前地质理论模型中(图5)，得到反演结果如图6所示。In order to further verify the effectiveness and efficiency of the full waveform inversion method of elastic waves detected in tunnel excavation based on multi-parameter weighted constraints and wave velocity profile spatial structure correction, the scheme proposed in the present invention was applied to a complex advanced geological theoretical model (Figure 5), and the inversion results were obtained as shown in Figure 6.
将以常规时间域多尺度弹性波全波形反演结果(如图7所示)为对比研究对象，对比图6与图7可以看出，本发明得到的全波形反演结果更接近真实模型。反演结果中陷落柱与断层破碎带构造两侧边界清晰且构造内部参数恢复较好，而且小尺寸断层岩性分界面也有恢复和显现；最明显的是结果中的假异常与扰动干扰基本压制完全，反演精度得到了大幅度的提高。对比反演结果与真实模型的数值差异程度，全波形反演精度提升了20％左右。The conventional time domain multi-scale elastic wave full waveform inversion results (as shown in FIG. 7) are used as comparative research objects. By comparing FIG. 6 with FIG. 7, it can be seen that the full waveform inversion results obtained by the present invention are closer to the real model. In the inversion results, the boundaries on both sides of the collapse column and the fault fracture zone structure are clear and the internal parameters of the structure are well restored, and the small-scale fault lithology interface is also restored and revealed; the most obvious thing is that the false anomalies and disturbance interference in the results are basically completely suppressed, and the inversion accuracy has been greatly improved. Comparing the numerical difference between the inversion results and the real model, the full waveform inversion accuracy has been improved by about 20%.
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的装置而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。 In this specification, each embodiment is described in a progressive manner, and each embodiment focuses on the differences from other embodiments. The same or similar parts between the embodiments can be referred to each other. For the device disclosed in the embodiment, since it corresponds to the method disclosed in the embodiment, the description is relatively simple, and the relevant parts can be referred to the method part.
对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。 The above description of the disclosed embodiments enables one skilled in the art to implement or use the present invention. Various modifications to these embodiments will be apparent to one skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention will not be limited to the embodiments shown herein, but rather to the widest scope consistent with the principles and novel features disclosed herein.

Claims (6)

  1. 一种基于多参数约束和结构校正的随掘巷道全波形反演方法，其特征在于，包括以下步骤：A full waveform inversion method for tunneling based on multi-parameter constraints and structural correction, characterized by comprising the following steps:
    步骤1、输入理论模拟数据或实测常规巷道地震超前探测数据，构建全波形反演初始模型，再采用多尺度弹性波全波形反演方法对初始模型进行单一尺度反演；Step 1: Input theoretical simulation data or measured conventional tunnel seismic advance detection data to construct a full waveform inversion initial model, and then use the multi-scale elastic wave full waveform inversion method to perform a single-scale inversion on the initial model;
    步骤2、对步骤1得到的单一尺度反演结果进行多参数加权约束结构校正，得到初步校正结果；Step 2, performing multi-parameter weighted constraint structure correction on the single-scale inversion result obtained in step 1 to obtain a preliminary correction result;
    步骤3、基于预先设定的限制条件，对初步校正结果进行一维波速剖面空间结构校正和平滑约束，得到二次校正结果；Step 3: Based on the preset restriction conditions, the one-dimensional wave velocity profile spatial structure correction and smoothing constraint are performed on the preliminary correction result to obtain the secondary correction result;
    步骤4、将步骤3得到的二次校正结果作为下一尺度的初始模型，继续按照步骤1的方式进行全波形反演；Step 4: Use the secondary correction result obtained in step 3 as the initial model of the next scale, and continue to perform full waveform inversion in the same manner as step 1;
    步骤5、重复执行步骤2～步骤4，直至所有尺度反演完成，得到巷道超前探测弹性波全波形反演结果。Step 5: Repeat steps 2 to 4 until all scale inversions are completed to obtain the full waveform inversion result of the elastic wave of the tunnel advance detection.
  2. 根据权利要求1所述的一种基于多参数约束和结构校正的随掘巷道全波形反演方法，其特征在于，步骤2中，按照下式进行多参数加权约束结构校正：
    The method for full waveform inversion of tunneling based on multi-parameter constraints and structural correction according to claim 1 is characterized in that, in step 2, multi-parameter weighted constraint structural correction is performed according to the following formula:
    其中，Δm为单次迭代模型更新量；Δm′为多参数加权结构校正后的模型更新量；m为初始模型；i,j分别表示z,x方向上网格节点的位置；mi,j p为纵波速度、横波速度或密度单一参数初始模型i,j网格节点位置的参数数值；nz为模型竖向网格点数，nx为模型横向网格点数；Vp为纵波速度；Vs为横波速度；Den为密度。 Among them, Δm is the model update amount of a single iteration; Δm′ is the model update amount after multi-parameter weighted structure correction; m is the initial model; i, j represent the positions of the grid nodes in the z and x directions, respectively; mi , jp are the parameter values of the i, j grid node positions of the initial model with a single parameter of longitudinal wave velocity, shear wave velocity or density; nz is the number of vertical grid points of the model, nx is the number of transverse grid points of the model; Vp is the longitudinal wave velocity; Vs is the shear wave velocity; Den is the density.
  3. 根据权利要求1所述的一种基于多参数约束和结构校正的随掘巷道全波形反演方法，其特征在于，步骤3中，预先设定的限制条件为：The method for full waveform inversion of tunneling based on multi-parameter constraints and structural correction according to claim 1 is characterized in that in step 3, the preset constraint conditions are:
    根据实际需要设置限制条件一，限制条件一以空间校正下限value值为判定依据，将低于空间校正下限value值的模型更新量进行压制，变为0；Set constraint condition 1 according to actual needs. Constraint condition 1 is based on the value of the lower limit of spatial correction. Suppress and change to 0;
    根据实际需要设置限制条件二，限制条件二以结构校正区域范围以及区域之间距离为判定依据，当模型更新量为负值区域两侧均存在一更新量为正值区域、或者当模型更新量为正值区域两侧均存在一更新量为负值区域，当区域范围以及区域之间距离满足限制条件二时，将负值区域两侧存在的正值区域内或者正值区域两侧存在的负值区域内的模型更新量进行压制，变为原来模型更新数值的1/5。The second restriction condition is set according to actual needs. The second restriction condition is based on the range of the structural correction area and the distance between areas. When there is a positive area on both sides of the negative area where the model update amount exists, or when there is a negative area on both sides of the positive area where the model update amount exists, when the area range and the distance between areas meet the second restriction condition, the model update amount in the positive area on both sides of the negative area or in the negative area on both sides of the positive area will be suppressed to 1/5 of the original model update value.
  4. 根据权利要求2所述的一种基于多参数约束和结构校正的随掘巷道全波形反演方法，其特征在于，步骤3包括：The method for full waveform inversion of tunneling based on multi-parameter constraints and structural correction according to claim 2 is characterized in that step 3 comprises:
    步骤301、选取一个纵轴网格坐标yi，沿巷道轴线提取一维波速剖面其中，x为横轴坐标；Step 301: Select a vertical grid coordinate y i and extract a one-dimensional wave velocity profile along the lane axis. Among them, x is the horizontal axis coordinate;
    步骤302、计算一维波速剖面对应的模型更新量与斜率 Step 302: Calculate the one-dimensional wave velocity profile The corresponding model update amount With slope
    步骤303、将低于空间校正下限value值的模型更新量进行压制，变为0，得到校正后的模型更新量 Step 303: Update the model value below the spatial correction lower limit value Suppress it and turn it into 0 to get the corrected model update amount
    步骤304、根据校正后的模型更新量按照限制条件二对一维波速剖面上的每个网格点进行逐一判定和校正，得到校正后的模型更新量 Step 304: Update the model according to the corrected value According to constraint condition 2, each grid point on the one-dimensional velocity profile is judged and corrected one by one to obtain the corrected model update amount.
    步骤305、根据校正后的模型更新量与斜率按照斜率符号变化情况对一维波速剖面上的每个网格点进行逐一判定，得到校正后的模型更新量 Step 305: Update the model according to the corrected value With slope According to the change of slope sign, each grid point on the one-dimensional wave velocity profile is judged one by one to obtain the corrected model update amount
    步骤306、根据校正后的模型更新量计算校正后的一维波速剖面 Step 306: Update the model according to the corrected value Calculate the corrected one-dimensional velocity profile
    步骤307、取下一个纵轴网格坐标yi+1，重复步骤302～步骤306，直至完成所有网格一维波速剖面的校正，得到一维波速剖面空间结构校正与平滑约束后的二次校正结果。Step 307, take the next vertical axis grid coordinate yi+1 , repeat steps 302 to 306, until the correction of the one-dimensional velocity profiles of all grids is completed, and obtain the secondary correction result of the one-dimensional velocity profile after spatial structure correction and smoothing constraint.
  5. 根据权利要求4所述的一种基于多参数约束和结构校正的随掘巷道全波形反演方法，其特征在于，步骤305包括：The method for full waveform inversion of tunneling based on multi-parameter constraints and structural correction according to claim 4 is characterized in that step 305 comprises:
    当斜率符号发生变化的相邻三个点之间存在两个模型更新量为零的点，该两个点前后网格点模型更新量符号发生变化，且该两个点之间没有斜率为0的点，则对斜率符号发生变化的相邻三个点之间区域内的网格点的模型更新量进行校正，变为第一个斜率符号变化点与第三个斜率符号变化点两点模型更新量的平均值，得到校正后的模型更新量 When there are two points with zero model update between the three adjacent points where the slope sign changes, the model update sign of the grid points before and after the two points changes, and there is no point with a slope of 0 between the two points, the model update of the grid points in the area between the three adjacent points where the slope sign changes is corrected to the average of the model update of the first slope sign change point and the third slope sign change point, and the corrected model update is obtained.
    再根据校正后的模型更新量继续校正，当斜率符号发生变化的不相邻的两个点之间范围内所有网格点模型更新量符号未发生改变，则对两个斜率符号变化点之间区域内的网格点的模型更新量进行校正，变为区域内所有网格点模型更新量的平均值，得到校正后的模型更新量 Then according to the corrected model update amount Continue to calibrate. When the model update signs of all grid points in the range between two non-adjacent points where the slope sign changes do not change, the model update of the grid points in the area between the two slope sign change points is corrected to the average of the model update of all grid points in the area, and the corrected model update is obtained.
  6. 根据权利要求1所述的一种基于多参数约束和结构校正的随掘巷道全波形反演方法，其特征在于，步骤5中的反演结果包括：纵波反演结果、横波反演结果和密度反演结果。 According to the full waveform inversion method of tunnel excavation based on multi-parameter constraints and structural correction as described in claim 1, it is characterized in that the inversion results in step 5 include: longitudinal wave inversion results, shear wave inversion results and density inversion results.
PCT/CN2023/113748 2022-10-13 2023-08-18 Excavation tunnel full-waveform inversion method based on multi-parameter constraint and structure correction WO2024078134A1 (en)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN202211254264.0 2022-10-13
CN202211254264.0A CN115542393A (en) 2022-10-13 2022-10-13 Tunneling-following roadway full-waveform inversion method based on multi-parameter constraint and structural correction

Publications (1)

Publication Number Publication Date
WO2024078134A1 true WO2024078134A1 (en) 2024-04-18

Family

ID=84733508

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2023/113748 WO2024078134A1 (en) 2022-10-13 2023-08-18 Excavation tunnel full-waveform inversion method based on multi-parameter constraint and structure correction

Country Status (2)

Country Link
CN (1) CN115542393A (en)
WO (1) WO2024078134A1 (en)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115542393A (en) * 2022-10-13 2022-12-30 安徽理工大学 Tunneling-following roadway full-waveform inversion method based on multi-parameter constraint and structural correction

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20130097504A (en) * 2012-02-24 2013-09-03 서울대학교산학협력단 Method and apparatus of estimating underground structure using a plurality of weighting value
CN112698390A (en) * 2020-11-11 2021-04-23 中国石油天然气股份有限公司 Pre-stack seismic inversion method and device
CN113552625A (en) * 2021-06-21 2021-10-26 中国地质科学院地球物理地球化学勘查研究所 Multi-scale full waveform inversion method for conventional land-domain seismic data
CN114624779A (en) * 2020-12-14 2022-06-14 中国石油化工股份有限公司 Pre-stack multi-parameter inversion method for balanced model constraint
CN115453624A (en) * 2022-10-13 2022-12-09 安徽理工大学 Mine following excavation seismic signal imaging method
CN115524744A (en) * 2022-10-13 2022-12-27 安徽理工大学 Roadway full-waveform inversion method based on cut-off frequency optimization and regularization constraints
CN115542393A (en) * 2022-10-13 2022-12-30 安徽理工大学 Tunneling-following roadway full-waveform inversion method based on multi-parameter constraint and structural correction

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20130097504A (en) * 2012-02-24 2013-09-03 서울대학교산학협력단 Method and apparatus of estimating underground structure using a plurality of weighting value
CN112698390A (en) * 2020-11-11 2021-04-23 中国石油天然气股份有限公司 Pre-stack seismic inversion method and device
CN114624779A (en) * 2020-12-14 2022-06-14 中国石油化工股份有限公司 Pre-stack multi-parameter inversion method for balanced model constraint
CN113552625A (en) * 2021-06-21 2021-10-26 中国地质科学院地球物理地球化学勘查研究所 Multi-scale full waveform inversion method for conventional land-domain seismic data
CN115453624A (en) * 2022-10-13 2022-12-09 安徽理工大学 Mine following excavation seismic signal imaging method
CN115524744A (en) * 2022-10-13 2022-12-27 安徽理工大学 Roadway full-waveform inversion method based on cut-off frequency optimization and regularization constraints
CN115542393A (en) * 2022-10-13 2022-12-30 安徽理工大学 Tunneling-following roadway full-waveform inversion method based on multi-parameter constraint and structural correction

Also Published As

Publication number Publication date
CN115542393A (en) 2022-12-30

Similar Documents

Publication Publication Date Title
WO2024078134A1 (en) Excavation tunnel full-waveform inversion method based on multi-parameter constraint and structure correction
Cameron et al. Seismic velocity estimation from time migration
CN105301639A (en) Speed field inversion method and device based on VSP double-weight travel time tomography
CN115524744A (en) Roadway full-waveform inversion method based on cut-off frequency optimization and regularization constraints
CN111239819B (en) Direct envelope inversion method with polarity based on seismic channel attribute analysis
CN106526675B (en) Tomography spatial parameter extraction method
CN116774292B (en) Seismic wave travel time determining method, system, electronic equipment and storage medium
US20230384473A1 (en) Method and system for advanced detection and optimization of tunnel resistivity based on depth resolution
CN104570066A (en) Method for building seismic inversion low-frequency models
CN103105622B (en) Based on the homotype ripple time difference positioning method of database technology
CN105425286A (en) Earthquake time-travelling acquisition method and crosshole earthquake time-travelling tomography method based on the earthquake time-travelling acquisition method
CN102590857A (en) True surface relief prestack depth domain two-way wave imaging method
CN109655890A (en) A kind of shallow mid-deep strata joint chromatography inversion speed modeling method of Depth Domain and system
CN103217715B (en) Multiple dimensioned regular grid Static Correction of Tomographic Inversion method
CN107817524A (en) The method and apparatus of three-dimensional seismic tomography
US20230095632A1 (en) Interpretive-guided velocity modeling seismic imaging method and system, medium and device
CN115453624A (en) Mine following excavation seismic signal imaging method
CN111123404B (en) Data fusion method for roadway advanced detection by earthquake and direct current method
CN105319594A (en) Fourier domain seismic data reconstruction method on the basis of least-square parametric inversion
Zhang et al. Improvement of microseismic source location during cavern excavation in faulted rock mass using fast marching method
CN111208568B (en) Time domain multi-scale full waveform inversion method and system
CN110687591B (en) Method for determining physical property parameters of coal bed and surrounding rock based on waveform matching of prior data
CN113589369B (en) Large cave depot group seismic stability assessment method
CN112147700A (en) Low-frequency model construction method and system for speed abnormal area
CN104570100A (en) Multi-wavelet Kirchhoff seismic data migration method

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 23876339

Country of ref document: EP

Kind code of ref document: A1