Full waveform inversion method for tunneling based on multi-parameter constraints and structural correction

Abstract

The present invention discloses a full waveform inversion method for tunneling while excavating based on multi-parameter constraints and structural correction, comprising the following steps: step 1, constructing a full waveform inversion initial model, and then using a multi-scale elastic wave full waveform inversion method to perform a single-scale inversion on the initial model; step 2, performing a multi-parameter weighted constraint structural correction on the single-scale inversion result to obtain a preliminary correction result; step 3, based on a pre-set restriction condition, performing a one-dimensional wave velocity profile spatial structure correction and a smoothing constraint on the preliminary correction result to obtain a secondary correction result; step 4, using the secondary correction result as the initial model of the next scale, and continuing to perform full waveform inversion; repeating steps 2 to 4 until all scale inversions are completed to obtain a full waveform inversion result. The present invention can effectively solve the problems of strong observation system restriction, small amount of detection data, small offset distance, and strong waveform inversion multi-solution in advance detection of tunneling while excavating, and improve the inversion effect.

Description

Technical Field

[0001] 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 Art

[0002] 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.

[0003] However, the current interpretation results of mine seismic advance detection imaging are mostly "arc drawing", with many false abnormal interfaces and low imaging accuracy (such as Figure 1 As shown in the figure, it is difficult to meet the geological support requirements for intelligent tunneling.

[0004] 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 complex structure imaging 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 a small amount of data and a small offset, which will lead to the enhancement of the multi-solution of full waveform inversion and the increase of physical parameter calculation error.

[0005] Therefore, how to provide a high-precision full-waveform inversion method suitable for mine tunnel detection mode is an urgent problem that technicians in this field need to solve. Summary of the invention

[0006] 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.

[0007] In order to achieve the above object, the present invention adopts the following technical solution:

[0008] A full waveform inversion method for tunneling based on multi-parameter constraints and structural correction includes the following steps:

[0009] 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;

[0010] 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;

[0011] 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;

[0012] 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;

[0013] 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.

[0014] Furthermore, in step 2, multi-parameter weighted constraint structure correction is performed according to the following formula:

[0015]

[0016] Wherein, Δ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; m i,j p is the parameter value of the initial model i, j grid node position of the single parameter of P-wave velocity, S-wave velocity or density; nz is the number of vertical grid points of the model, nx is the number of horizontal grid points of the model; V p is the longitudinal wave velocity; V s is the shear wave velocity; Den is the density.

[0017] Furthermore, in step 3, the preset restriction conditions are:

[0018] 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;

[0019] 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.

[0020] Furthermore, step 3 includes:

[0021] Step 301, select a vertical axis grid coordinate y i , extract the one-dimensional wave velocity profile along the tunnel axis Among them, x is the horizontal axis coordinate;

[0022] Step 302: Calculate the one-dimensional wave velocity profile The corresponding model update amount With slope

[0023] 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

[0024] Step 304: Update the model according to the corrected value According to constraint condition 2, each grid point on the one-dimensional wave velocity profile is judged and corrected one by one to obtain the corrected model update amount.

[0025] 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 velocity profile is judged one by one to obtain the corrected model update amount

[0026] Step 306: Update the model according to the corrected value Calculate the corrected one-dimensional velocity profile

[0027] Step 307: Take the next vertical axis grid coordinate y i+1 , repeat steps 302 to 306 until the correction of the one-dimensional velocity profiles of all grids is completed, and the secondary correction result of the one-dimensional velocity profile after spatial structure correction and smoothing constraint is obtained.

[0028] Further, step 305 includes:

[0029] 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.

[0030] Then update the amount of the corrected model 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.

[0031] Furthermore, the inversion results in step 5 include: longitudinal wave inversion results, shear wave inversion results and density inversion results.

[0032] 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:

[0033] (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.

[0034] (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.

[0035] (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.

[0036] (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%.

[0037] (5) The present invention improves the inversion accuracy without increasing the amount of full waveform inversion calculations. BRIEF DESCRIPTION OF THE DRAWINGS

[0038] 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.

[0039] Figure 1 The conventional mine earthquake advance detection imaging results provided by the present invention;

[0040] Figure 2 A flow chart of the full waveform inversion method for elastic wave detection in tunneling provided by the present invention;

[0041] Figure 3 A schematic diagram of determining the update amount and slope change of the one-dimensional wave velocity profile model provided by the present invention;

[0042] Figure 4 A schematic diagram of the structure of the initial model provided by the present invention;

[0043] Figure 5 The complex advanced geological theoretical model provided by the present invention;

[0044] Figure 6 The inversion result obtained by the method of the present invention is provided by the present invention;

[0045] Figure 7 The present invention provides a conventional time domain multi-scale elastic wave full waveform inversion result. DETAILED DESCRIPTION

[0046] 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.

[0047] like Figure 2 As shown, the 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:

[0048] Step 1: Input theoretical simulation data or measured conventional tunnel seismic advance detection data to construct the full waveform inversion initial model (such as Figure 4 As shown in Figure 3, the initial model is inverted at a single scale using the multi-scale elastic wave full waveform inversion method. The inversion includes three sets of parameters, namely the P-wave velocity model, the S-wave velocity model and the density model.

[0049] 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;

[0050] 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;

[0051] 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;

[0052] 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.

[0053] In the present invention, the seismic source is distributed in the middle of the coal seam, and a plurality of 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 initial inversion model is set according to the medium conditions of the uniformly layered coal-bearing strata.

[0054] 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.

[0055] 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:

[0056] (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;

[0057] (2) Although P-wave, S-wave and density have great differences in physical parameters, they have the same geological structural characteristics;

[0058] (3) There are false anomaly areas with opposite attribute characteristics to the abnormal geological structure near both sides of the boundary;

[0059] (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;

[0060] (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".

[0061] 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.

[0062] Next, the correction method of step 2 and step 3 is described in detail.

[0063] In a specific embodiment, in step 2, multi-parameter weighted constraint structure correction is performed according to the following formula:

[0064]

[0065] Wherein, Δm is the single iteration model update amount (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; m i,j p is the parameter value of the initial model i, j grid node position of the single parameter of P-wave velocity, S-wave velocity or density; nz is the number of vertical grid points of the model, nx is the number of horizontal grid points of the model; V p is the longitudinal wave velocity; V s is the shear wave velocity; Den is the density

[0066] In a specific embodiment, in step 3, the preset restriction condition is:

[0067] 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.

[0068] Constraint 2: Constraint 2 is based on the range of the structural correction area and the distance between areas. When there is a positive update area on both sides of a negative model update area, or when there is a negative update area on both sides of a positive model update area, 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.

[0069] 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".

[0070] "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".

[0071] 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 "positive region" in the middle.

[0072] "Distance between regions": When there is a positive update 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.

[0073] Specifically, step 3 includes:

[0074] Step 301, select a vertical axis grid coordinate y i , extract the one-dimensional wave velocity profile along the tunnel axis Among them, x is the horizontal axis coordinate;

[0075] Step 302: Calculate the one-dimensional wave velocity profile The corresponding model update amount With slope

[0076] 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.

[0077] 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.

[0078] 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

[0079] Step 304: Update the model according to the corrected value According to constraint condition 2, each grid point on the one-dimensional wave velocity profile is judged and corrected one by one to obtain the corrected model update amount.

[0080] 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 velocity profile is judged one by one to obtain the corrected model update amount

[0081] When there are two points with zero model update between the three adjacent points where the slope sign changes (such as Figure 3 As shown in the figure, if the sign of the model update amount 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 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.

[0082] 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.

[0083] 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.

[0084] Step 306: Update the model according to the corrected value Calculate the corrected one-dimensional velocity profile

[0085] Step 307: Take the next vertical axis grid coordinate y i+1 , repeat steps 302 to 306 until the correction of the one-dimensional velocity profiles of all grids is completed, and the secondary correction result of the one-dimensional velocity profile after spatial structure correction and smoothing constraint is obtained.

[0086] In order to further verify the effectiveness and efficiency of the full waveform inversion method of elastic wave detection in tunneling based on multi-parameter weighted constraints and wave velocity profile spatial structure correction, the scheme proposed in the present invention is applied to a complex advanced geological theoretical model ( Figure 5 ), and the inversion result is as follows Figure 6 shown.

[0087] The conventional time domain multi-scale elastic wave full waveform inversion results (such as Figure 7 The comparison is shown in Figure 2. Figure 6 and Figure 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 of the collapse column and the fault fracture zone on both sides are clear and the internal parameters of the structure are well restored. In addition, the lithological interface of the small-scale fault is also restored and revealed. The most obvious thing is that the false anomalies and disturbances in the results are basically suppressed completely, and the inversion accuracy has been greatly improved. Comparing the numerical difference between the inversion results and the real model, the accuracy of the full waveform inversion has been improved by about 20%.

[0088] 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.

[0089] 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

1. A full waveform inversion method for tunneling based on multi-parameter constraints and structural correction, characterized in that: The following steps are involved: 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; 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 to obtain the full waveform inversion result of the elastic wave of the tunnel advance detection.

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 structure correction is performed according to the following formula: Wherein, Δ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; m i,j p is the parameter value of the initial model i, j grid node position of the single parameter of P-wave velocity, S-wave velocity or density; nz is the number of vertical grid points of the model, nx is the number of horizontal grid points of the model; V p is the longitudinal wave velocity; V s is the shear wave velocity; Den is the density.

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 pre-set constraints are: 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; 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. 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 includes: Step 301, select a vertical axis grid coordinate y i , extract the one-dimensional wave velocity profile along the tunnel axis Among them, 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 wave 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 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 y i+1 , repeat steps 302 to 306 until the correction of the one-dimensional velocity profiles of all grids is completed, and the secondary correction result of the one-dimensional velocity profile after spatial structure correction and smoothing constraint is obtained.

5. 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 includes: 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 update the amount of the corrected model 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. The method of full waveform inversion of tunneling based on multi-parameter constraints and structural correction according to claim 1 is characterized in that: The inversion results in step 5 include: longitudinal wave inversion results, shear wave inversion results and density inversion results.

Images