Full waveform inversion method for tunnel based on cutoff frequency optimization and regularization constraints

Abstract

The present invention discloses a tunnel full waveform inversion method based on cutoff frequency optimization and regularization constraint, including: constructing a multi-scale cutoff frequency optimization scheme for elastic waves in a tunnel advance detection mode; determining an optimal regularization constraint weight factor; using input seismic data and the optimal regularization constraint weight factor, as well as a multi-scale cutoff frequency optimization scheme, performing a single-scale inversion of a tunnel elastic wave full waveform inversion method based on total variation regularization constraint; performing a wave velocity structure correction on a single-scale inversion result; using the wave velocity structure correction result as an initial model for the next scale to continue full waveform inversion until all scale inversions in the multi-scale cutoff frequency optimization scheme are completed, and obtaining a tunnel advance detection elastic wave full waveform inversion result. The present invention can better alleviate the nonlinearity and ill-posedness problems of full waveform inversion, improve the inversion effect, and achieve high-precision imaging of abnormal geological structures ahead of tunneling.

Description

Technical Field

[0001] The present invention relates to the field of mining geophysical exploration imaging technology, and in particular to a tunnel full waveform inversion method based on cutoff frequency optimization and regularization constraints, and more specifically to a tunnel advance detection elastic wave full waveform inversion method based on multi-scale cutoff frequency optimization and regularization constraints under wave velocity structure correction. Background Art

[0002] Intelligent coal mining is the core technical support for adapting to the development trend of the modern industrial technological revolution, achieving the "dual carbon" goals, ensuring national energy security, and realizing high-quality development of the coal industry. As one of the two core links in coal mine production, the demand for the development of intelligent tunneling is extremely urgent. However, during tunnel excavation, 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 ensuring the safety of intelligent coal production. It is the basic data source for realizing geological prediction, disturbance perception and risk assessment before, during and after tunnel excavation construction. It is the prerequisite and guarantee for the implementation of all key technologies for intelligent tunneling. However, the current interpretation results of mine seismic advance detection imaging are mostly "arcing", 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.

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

[0004] 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 conventional time domain multi-scale elastic wave full waveform inversion results cannot effectively guide the safe production of coal mine excavation. False anomalies and disturbance noise are relatively serious, and there is still a problem of inversion discontinuity for abnormal geological structures, which urgently needs optimization and improvement.

[0005] Therefore, how to provide a tunnel full waveform inversion method based on cutoff frequency optimization and regularization constraints to solve the problem of discontinuous inversion of abnormal geological structures ahead of tunnel excavation and achieve high-precision imaging 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 elastic waves in advance detection of tunnels based on multi-scale cutoff frequency optimization and regularization constraints under wave velocity structure correction, which solves the problem of discontinuous inversion of geological structures in the inversion results of the full-waveform inversion method with structural correction provided by the inventor, and better alleviates the nonlinearity and instability of the full-waveform inversion, effectively solves the problems of small amount of detection data, small offset distance, and strong multi-solution of waveform inversion in advance detection of tunnels, improves the inversion effect, and realizes 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 tunnel full waveform inversion method based on cutoff frequency optimization and regularization constraints includes the following steps:

[0009] S1: Combined with the initial scale cutoff frequency, the longitudinal wave velocity, the shear wave velocity and the density parameters are comprehensively considered to construct the optimization scheme of the elastic wave multi-scale cutoff frequency in the tunnel advance detection mode;

[0010] S2: Determine the optimal regularization constraint weight factor based on the full waveform inversion initial model and the input seismic data;

[0011] S3: using the input seismic data and the optimal regularization constraint weight factor described in S2, and the multi-scale cutoff frequency optimization scheme described in S1, a single-scale inversion of the tunnel elastic wave full waveform inversion method based on total variation regularization constraint is performed;

[0012] S4: Perform velocity structure correction on the single-scale inversion result obtained in S3;

[0013] S5: Use the velocity structure correction result obtained in S4 as the initial model for the next scale and continue to perform full waveform inversion in the same manner as S3;

[0014] S6: Repeat S4 to S5 until all scale inversions in the multi-scale cutoff frequency optimization scheme described in S1 are completed, and the full waveform inversion result of the tunnel advance detection elastic wave is obtained.

[0015] Preferably, in S1, the preferred scheme for elastic wave multi-scale cutoff frequency in the lane advance detection mode is constructed as follows:

[0016]

[0017] Where: f n With f n+1 is the cutoff frequency of two adjacent scales; Δf is the cutoff frequency interval of adjacent scales; is the minimum cosine value of the reflection angle generated by the incident wave in the target layer, and the maximum target layer depth of the inversion is z max , the maximum offset distance is h; is the comprehensive background velocity of coal-bearing strata; f is the main frequency of the full waveform inversion source; f1 is the initial scale cutoff frequency.

[0018] Preferably, S2 specifically includes:

[0019] According to the known geological data of the tunnel, the initial model of full waveform inversion is constructed, and the theoretical simulation data or the measured conventional tunnel seismic advance detection data are input;

[0020] The optimal regularization constraint weight factor is obtained experimentally based on the input seismic data.

[0021] Preferably, S3 specifically includes:

[0022] The single-scale inversion of the tunnel elastic wave full waveform inversion method based on total variation regularization constraints is performed according to the following equation:

[0023]

[0024] Where: E(m) is the objective function; λ is the weight coefficient of the regularization term; m represents the P-wave and S-wave velocity models and the density model; u cal and u obs Respectively represent the simulated and observed particle displacement data; m k represents the kth iteration model parameter, α k-1 represents the step size, δm k-1 represents the search direction, i.e. the model update amount; L(p N ,q N ) is the total variation regularization constraint update term, p N ,q N is the regularization constraint operator.

[0025] Preferably, the single-scale inversion of the tunnel elastic wave full waveform inversion method based on total variation regularization constraint in S3 specifically includes:

[0026] S301: Calculate the objective function value according to E(m);

[0027] S302: Using the objective function value obtained in S301, use the parabola fitting method and the L-BFGS method to calculate the step length α k-1 and search direction δm k-1 ;

[0028] S303: Based on step size α k-1 and search direction δm k-1 Calculate the total variation regularization constraint update term L(p N ,q N );

[0029] S304: According to formula m k =(m k-1 +α k-1 δm k-1 )-λ·L(p N ,q N ) Iteratively update the model, λ is the weight coefficient of the regularization term, m k represents the kth iteration model parameter, m k-1 represents the k-1th iteration model parameter, and the single-scale inversion result of the tunnel elastic wave full waveform inversion method based on total variation regularization constraint is obtained.

[0030] Preferably, the S4 specifically includes:

[0031] S401: The single-scale inversion result obtained in S3 is subjected to multi-parameter weighted constraint structure correction processing according to the following equation, and the result after multi-parameter weighted constraint structure correction is obtained:

[0032]

[0033] Where: Δm i,j p is the single iteration model update amount; Δm′ i,j p is the model update amount after multi-parameter weighted structure correction; m i,j p is the multi-parameter initial model; p = V p ,V s , Den represents the medium longitudinal wave velocity, shear wave velocity, and density physical parameters, nx is the number of horizontal grids in the model, nz is the number of vertical grids in the model, and i, j represent the positions of grid nodes in the z and x directions respectively;

[0034] S402: Perform one-dimensional velocity profile spatial structure correction and smoothing constraint on the result of multi-parameter weighted constraint structure correction obtained in S401.

[0035] Preferably, the S402 specifically includes:

[0036] S4021: Set the restriction conditions of the spatial structure correction of the one-dimensional wave velocity profile, including the lower limit value of the spatial correction, the range of the structural correction area, and the distance between the areas;

[0037] S4022: Select a vertical grid coordinate y i, extract the one-dimensional wave velocity profile along the tunnel axis Among them, x is the horizontal axis coordinate;

[0038] S4023: Calculate one-dimensional wave velocity profile The corresponding model update amount With slope

[0039] S4024: 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

[0040] S4025: Update the amount based on the corrected model 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.

[0041] S4026: Update the amount based on the corrected model 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

[0042] S4027: Update the amount based on the corrected model Calculate the corrected one-dimensional velocity profile

[0043] S4028: Take the next vertical grid coordinate y i+1 , repeat S4023 to S4027 until the correction of the one-dimensional velocity profiles of all grids is completed, and the result after the velocity structure correction is obtained.

[0044] Preferably, the S4026 specifically includes:

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

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

[0047] The present invention also provides a computer-readable medium, which stores instructions. When the instructions are executed on the readable medium, the readable medium executes the tunnel full waveform inversion method based on cutoff frequency optimization and regularization constraints.

[0048] The present invention also provides an electronic device, including a memory and a processor, wherein the memory stores a computer program that can be run on the processor, and when the processor executes the computer program, the tunnel full waveform inversion method based on cutoff frequency optimization and regularization constraints is implemented.

[0049] It can be seen from the above technical solutions that, compared with the prior art, the present invention discloses a tunnel full waveform inversion method based on cutoff frequency optimization and regularization constraints, which has the following beneficial effects:

[0050] (1) The present invention provides 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;

[0051] (2) The present invention realizes 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;

[0052] (3) The present invention provides an optimal multi-scale cutoff frequency selection strategy, which saves memory space and computing time while ensuring inversion accuracy;

[0053] (4) The present invention adopts total variation regularization constraints to solve the problem of discontinuous inversion of geological structures in the inversion results of the structure-corrected full waveform inversion method, and better alleviates the nonlinearity and ill-posedness problems of the full waveform inversion;

[0054] (5) 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 25%. BRIEF DESCRIPTION OF THE DRAWINGS

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

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

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

[0058] Figure 3 The complex advanced geological theoretical model provided by the present invention;

[0059] Figure 4 The conventional time domain multi-scale elastic wave full waveform inversion result provided by the present invention;

[0060] Figure 5 The full waveform inversion result of elastic waves in the tunnel based on wave velocity structure correction provided by the present invention;

[0061] Figure 6 A flow chart of a tunnel full waveform inversion method based on cutoff frequency optimization and regularization constraints provided by the present invention;

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

[0063] Figure 8 The invention provides an inversion result obtained by using the method of the invention. DETAILED DESCRIPTION

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

[0065] In a first aspect, an embodiment of the present invention discloses a tunnel full waveform inversion method based on cutoff frequency optimization and regularization constraints to construct a complex advanced geological theoretical model (such as Figure 2-Figure 3The invention is described by taking the example of FIG. 1 as an example. In addition to the embodiments, the invention is also applicable to other models and actual data. The invention is based on the conventional time domain multi-scale elastic wave full waveform inversion method. The inversion includes three sets of parameters, namely the longitudinal wave velocity model, the shear wave velocity model and the density model.

[0066] Figure 2 The initial model shown is a real model of longitudinal waves. The earthquake source is distributed in the middle of the coal seam. Multiple horizontal component and vertical component geophones are arranged linearly with the earthquake source and are also distributed in the middle of the coal seam to receive seismic signals. The layout of the observation system is as follows: Figure 3 The geophone and source are shown in the figure. The initial inversion model is set according to the uniform layered coal-bearing stratum medium.

[0067] The conventional time domain multi-scale elastic wave full waveform inversion results (such as Figure 4 The inversion technology used is the most practical and advanced technology in surface seismic exploration. However, its results cannot effectively guide the safe production of coal mine excavation. False anomalies and disturbance noise are relatively serious, and optimization and improvement are urgently needed. Figure 4 (a) is the result of longitudinal wave inversion, Figure 4 (b) is the shear wave inversion result, Figure 4 (c) is the density inversion result.

[0068] 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. A tunnel elastic wave full waveform inversion method based on wave velocity structure correction is proposed, which constrains the direction of the full waveform inversion and obtains the inversion results that meet the preset structural features.

[0069] The tunnel elastic wave full waveform inversion method based on wave velocity structure correction provided includes two technical methods, namely:

[0070] (1) Structural correction method based on multi-parameter weighted constraints;

[0071] (2) One-dimensional velocity profile spatial structure correction and smoothing constraint method.

[0072] The proposed tunnel elastic wave full waveform inversion method based on wave velocity structure correction is applied to complex advanced geological theoretical models ( Figure 3 ), and the inversion result is as follows Figure 5 shown. Figure 5 (a) is the result of longitudinal wave inversion, Figure 5 (b) is the shear wave inversion result, Figure 5 (c) is the density inversion result.

[0073] contrast Figure 4 and Figure 5 It can be seen that the full waveform inversion results of tunnel elastic waves based on wave velocity structure correction are closer to the real model. In the inversion results, the boundaries of the collapse column and the fault fracture zone structure are clear, the internal parameters of the structure are well restored, and the small-scale fault lithology interface is also restored and revealed; the most obvious is that the false anomalies and disturbances in the results are basically suppressed completely, and the inversion accuracy has been greatly improved. Then, the numerical difference between the inversion results and the real model is compared, and the full waveform inversion accuracy is improved by about 20%.

[0074] However, the only drawback is that there is an inversion discontinuity in the abnormal geological structure in the full waveform inversion results of tunnel elastic waves based on wave velocity structure correction, which is caused by errors in the structure correction process.

[0075] Therefore, the present invention provides a full waveform inversion method for elastic waves of advance detection in tunnels based on multi-scale cutoff frequency optimization and regularization constraints under wave velocity structure correction, comprising the following steps:

[0076] S1: Combined with the initial scale cutoff frequency, the multi-scale cutoff frequency optimization scheme of elastic waves in the tunnel advance detection mode is constructed by comprehensively considering the longitudinal wave velocity, shear wave velocity and density parameters;

[0077] S2: Determine the optimal regularization constraint weight factor based on the full waveform inversion initial model and the input seismic data;

[0078] S3: Using the seismic data input in S2 and the optimal regularization constraint weight factor, as well as the multi-scale cutoff frequency optimization scheme in S1, a single-scale inversion of the tunnel elastic wave full waveform inversion method based on total variation regularization constraint is performed;

[0079] S4: Perform velocity structure correction on the single-scale inversion result obtained in S3;

[0080] S5: Use the velocity structure correction result obtained in S4 as the initial model for the next scale and continue to perform full waveform inversion in the same manner as S3;

[0081] S6: Repeat S4 to S5 until all scale inversions in the multi-scale cutoff frequency optimization scheme in S1 are completed, and the full waveform inversion result of the elastic wave of the tunnel advance detection is obtained.

[0082] In the present invention, the so-called multi-scale inversion is to divide the seismic data into different scales by low-pass filter, so as to realize the successive inversion of multiple scales. In the preferred scheme of multi-scale cutoff frequency, the cutoff frequency is the frequency limit of low-pass filter filtering, and each scale corresponds to a different cutoff frequency.

[0083] For multi-scale full waveform inversion, it is very important to choose an appropriate frequency scale. Too few frequency scales will lead to unstable inversion, incomplete inversion information, and easy to fall into local minima; too many frequency scales will not increase the accuracy of inversion, but will only increase the amount of redundant calculations, waste memory space and calculation time. Therefore, how to choose an optimal set of multi-scale frequency band ranges while taking into account the stability and efficiency of calculation is a key issue in time domain multi-scale inversion.

[0084] In one embodiment, the preferred scheme for the elastic wave multi-scale cutoff frequency in the lane advance detection mode in S1 is constructed as follows:

[0085]

[0086] Where: f n With f n+1 is the cutoff frequency of two adjacent scales; Δf is the cutoff frequency interval of adjacent scales; is the minimum cosine value of the reflection angle generated by the incident wave in the target layer, and the maximum target layer depth of the inversion is z max , unlike surface seismic exploration, tunnel advance detection is similar to single offset detection for the front of excavation, with a maximum offset of h; is the comprehensive background velocity of coal-bearing strata; f is the main frequency of the full waveform inversion source; f1 is the initial scale cutoff frequency, that is, the frequency value of the cutoff frequency of the first scale in the multi-scale full waveform inversion process; N represents the number of layers of the geological model.

[0087] Different from the selection strategy of full waveform inversion cutoff frequency in conventional acoustic wave media, the selection of elastic wave multi-scale cutoff frequency in tunnel advance detection mode needs to be adjusted according to the following situations:

[0088] (1) The physical parameters of P-wave velocity, S-wave velocity and density are different, that is, the wave number corresponding to the same frequency is different. In other words, in the elastic wave multi-scale full waveform inversion, the cutoff frequencies selected for P-wave velocity, S-wave velocity and density should also be different. However, in the elastic wave full waveform inversion, it is impossible to perform scale filtering on the P-wave velocity and S-wave density separately. Therefore, it is necessary to comprehensively consider the P-wave, S-wave and density to obtain a comprehensive Δf for calculation.

[0089] (2) Since coal-bearing strata are typical layered strata, the background velocity co used should be a comprehensive value.

[0090] (3) Within a certain frequency range, when the initial scale cutoff frequency is higher, the structural features of the structure are better portrayed and the boundaries are clearer; when the initial scale cutoff frequency is lower, the internal parameter details of the structure are better restored. It is necessary to comprehensively consider the structural characteristics and parameter details to determine the initial cutoff frequency.

[0091] In one embodiment, S2 specifically includes:

[0092] According to the known geological data of the tunnel, the initial model of full waveform inversion is constructed, and the theoretical simulation data or the measured conventional tunnel seismic advance detection data are input;

[0093] According to the input seismic data, the optimal regularization constraint weight factor is obtained experimentally.

[0094] In one embodiment, adding a regularization term to the objective function to constrain the inversion process is a way to alleviate the nonlinearity and ill-posedness of full waveform inversion, and is an important way to constrain the full waveform inversion results. The total variation regularization method can achieve noise suppression by solving the nonlinear minimization problem of the total variation constraint while retaining the original structural characteristics of the structure; and the total variation regularization is a sparse constrained full waveform inversion method, which can highlight and retain the boundary information of the inversion result.

[0095] The tunnel elastic wave full waveform inversion method based on total variation regularization constraint has limitations in use, which are:

[0096] (1) The velocity structure correction method cannot directly participate in every iteration of the full waveform inversion. It needs to be combined with the multi-scale inversion to perform structural correction on the results of the single-scale inversion, and use the structurally corrected model as the initial model of the next scale to continue the full waveform inversion. The tunnel elastic wave full waveform inversion method based on total variation regularization constraints needs to participate in every iteration of the full waveform inversion.

[0097] (2) The weight factor of the regularization constraint has a great influence on the constraint effect of the regularization term, but this factor does not have a specific optimal value. Therefore, it is necessary to obtain the optimal weight factor through experiments before the formal inversion.

[0098] In S3, single-scale inversion of tunnel elastic wave full waveform inversion method based on total variation regularization constraint is performed according to the following equation:

[0099]

[0100] Where: E(m) is the objective function; λ is the weight coefficient of the regularization term; m represents the P-wave and S-wave velocity models and the density model; u cal and u obs Respectively represent the simulated and observed particle displacement data; m k represents the kth iteration model parameter, α k-1 represents the step size, δm k-1 represents the search direction, i.e. the model update amount; L(p N ,q N ) is the total variation regularization constraint update term, pN ,q N is the regularization constraint operator.

[0101] In this embodiment, the specific execution steps of the single-scale inversion of the tunnel elastic wave full waveform inversion method based on total variation regularization constraint include:

[0102] S301: Calculate the objective function value according to E(m);

[0103] S302: Using the objective function value obtained in S301, use the parabola fitting method and the L-BFGS method to calculate the step length α k-1 and search direction δm k-1 ;

[0104] S303: Based on step size α k-1 and search direction δm k-1 Calculate the total variation regularization constraint update term L(p N ,q N );

[0105] S304: According to formula m k =(m k-1 +α k-1 δm k-1 )-λ·L(p N ,q N ) Iteratively update the model, λ is the weight coefficient of the regularization term, m k represents the kth iteration model parameter, m k-1 represents the k-1th iteration model parameter, and the single-scale inversion result of the tunnel elastic wave full waveform inversion method based on total variation regularization constraint is obtained.

[0106] Among them, the total variation regularization constraint update term L(p N ,q N ):

[0107] ① Set the maximum number of computational iterations N for the total variation regularization constraint update term during a single iteration of full waveform inversion.

[0108] ②

[0109] Where p and q are regularization constraint operators, N is the maximum number of calculation stacking times of the total variation regularization constraint update term in a single iteration of full waveform inversion, i = 1…nz, j = 1…nx, nx is the number of horizontal grids in the model, and nz is the number of vertical grids in the model.

[0110] and pass and Obtained from multiple superposition calculations of e=1…N.

[0111] Where:

[0112] Where:

[0113] Where:

[0114] Where:

[0115] Where:

[0116] In the formula, let h 1 =1, is a zero matrix, h is a regularized update operator, and C and D are regularized iteration operators.

[0117] In one embodiment, S4 specifically includes:

[0118] S401: The single-scale inversion result obtained in S3 is subjected to multi-parameter weighted constraint structure correction processing according to the following equation, and the result after multi-parameter weighted constraint structure correction is obtained:

[0119]

[0120] Where: Δm i,j p is the single iteration model update amount; Δm′ i,j p is the model update amount after multi-parameter weighted structure correction; m i,j p is the multi-parameter initial model; p = V p ,V s , Den represents the medium longitudinal wave velocity, shear wave velocity, and density physical parameters, nx is the number of horizontal grids in the model, nz is the number of vertical grids in the model, and i, j represent the positions of grid nodes in the z and x directions respectively;

[0121] S402: Perform one-dimensional velocity profile spatial structure correction and smoothing constraint on the result of multi-parameter weighted constraint structure correction obtained in S401.

[0122] In this embodiment, S402 specifically includes:

[0123] S4021: Set the constraint conditions for the spatial structure correction of the one-dimensional velocity profile, including the lower limit value of the spatial correction, the range of the structural correction area, and the distance between the areas; Setting the constraint conditions required for the spatial structure correction of the one-dimensional velocity profile is related to the quality of the correction effect during the application of the present invention, and is also related to the degree of human intervention in the inversion results. There is no specific setting value, and it needs to be set accordingly according to the actual situation. The pre-set constraint conditions are specifically:

[0124] Constraint 1: Setting the lower limit value of spatial correction;

[0125] The spatial correction lower limit value is used as the basis for judgment, and the model update amount below the preset lower limit value is Suppress and directly change to 0. The lower limit value of spatial correction needs to be set according to actual needs. For example, the value is set to be less than 5% of the initial model parameter.

[0126] Constraint 2: Setting of structural correction area range and distance between areas;

[0127] Based on the structural correction area range and the distance between areas, when there is a positive area on both sides of the negative area, or when there is a negative area on both sides of the positive area, when the area range and the distance between areas meet the setting 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 is suppressed to 1 / 5 of the original model update value. Before applying this embodiment, the area range and the distance between areas need to be set according to actual needs.

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

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

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

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

[0132] S4022: Select a vertical grid coordinate y i , extract the one-dimensional wave velocity profile along the tunnel axis Among them, x is the horizontal axis coordinate;

[0133] S4023: Calculate one-dimensional wave velocity profile The corresponding model update amount With slope

[0134] Among them, the model update amount is: is the one-dimensional wave velocity profile at the kth iteration, is the one-dimensional wave velocity profile at the k-1th iteration.

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

[0136] S4024: 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

[0137] S4025: Update the amount based on the corrected model 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.

[0138] S4026: Update the amount based on the corrected model 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

[0139] When there are two points with zero model update between the three adjacent points where the slope sign changes (such as Figure 3As 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.

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

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

[0142] S4027: Update the amount based on the corrected model Calculate the corrected one-dimensional velocity profile

[0143] S4028: Take the next vertical grid coordinate y i+1 , repeat S4023 to S4027 until the correction of the one-dimensional velocity profiles of all grids is completed, and the result after the velocity structure correction is obtained.

[0144] The embodiments of the present invention take the constructed complex advanced geological theoretical model as an example to illustrate the invention content. In addition to the embodiments, the invention is also applicable to other models and actual data. The present invention is based on the conventional time domain multi-scale elastic wave full waveform inversion method to carry out innovation. The inversion includes three sets of parameters, namely, the P-wave velocity model, the S-wave velocity model and the density model.

[0145] The second aspect of the embodiment of the present invention further discloses a computer-readable medium, which stores instructions. When the instructions are executed on the readable medium, the readable medium executes the tunnel full waveform inversion method based on cutoff frequency optimization and regularization constraints disclosed in the first aspect.

[0146] The third aspect of an embodiment of the present invention further discloses an electronic device, including a memory and a processor, wherein the memory stores a computer program that can be run on the processor, and when the processor executes the computer program, the tunnel full waveform inversion method based on cutoff frequency optimization and regularization constraints disclosed in the first aspect is implemented.

[0147] 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 3 ), and the inversion result is as follows Figure 8 shown. Figure 8 (a) is the result of longitudinal wave inversion, Figure 8 (b) is the shear wave inversion result, Figure 8 (c) is the density inversion result.

[0148] contrast Figure 4 , Figure 5 and Figure 8 It can be seen that the full waveform inversion result obtained by the embodiment of the present invention is closer to the real model. In the inversion result, the boundaries on both sides of the structure are clear and the internal parameters of the structure are well restored. The false anomalies and disturbance interference in the inversion result are basically completely suppressed, and the inversion accuracy has been greatly improved. Moreover, the existing discontinuity phenomenon of geological structure inversion is also eliminated. Then, by comparing the numerical difference between the inversion result and the real model, the accuracy of the full waveform inversion is improved by about 25% compared with the conventional method.

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

[0150] 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 tunnel full waveform inversion method based on cutoff frequency optimization and regularization constraints, characterized in that: The following steps are involved: S1: Combined with the initial scale cutoff frequency, the longitudinal wave velocity, the shear wave velocity and the density parameters are comprehensively considered to construct the optimization scheme of the elastic wave multi-scale cutoff frequency in the tunnel advance detection mode; S2: Determine the optimal regularization constraint weight factor based on the full waveform inversion initial model and the input seismic data; S3: using the input seismic data and the optimal regularization constraint weight factor described in S2, and the multi-scale cutoff frequency optimization scheme described in S1, a single-scale inversion of the tunnel elastic wave full waveform inversion method based on total variation regularization constraint is performed; S4: Perform velocity structure correction on the single-scale inversion result obtained in S3; S5: Use the velocity structure correction result obtained in S4 as the initial model for the next scale and continue to perform full waveform inversion in the same manner as S3; S6: Repeat S4 to S5 until all scale inversions in the multi-scale cutoff frequency optimization scheme described in S1 are completed, and the full waveform inversion result of the tunnel advance detection elastic wave is obtained.

2. A tunnel full waveform inversion method based on cutoff frequency optimization and regularization constraint according to claim 1, characterized in that: In S1, the optimization scheme of elastic wave multi-scale cutoff frequency in lane advance detection mode is constructed as follows: Where: f n With f n+1 is the cutoff frequency of two adjacent scales; Δf is the cutoff frequency interval of adjacent scales; is the minimum cosine value of the reflection angle generated by the incident wave in the target layer, and the maximum target layer depth of the inversion is z max , the maximum offset distance is h; is the comprehensive background velocity of coal-bearing strata; f is the main frequency of the full waveform inversion source; f1 is the initial scale cutoff frequency.

3. The tunnel full waveform inversion method based on cutoff frequency optimization and regularization constraint according to claim 1 is characterized in that: The S2 specifically includes: According to the known geological data of the tunnel, the initial model of full waveform inversion is constructed, and the theoretical simulation data or the measured conventional tunnel seismic advance detection data are input; The optimal regularization constraint weight factor is obtained experimentally based on the input seismic data.

4. The tunnel full waveform inversion method based on cutoff frequency optimization and regularization constraint according to claim 1 is characterized in that: The S3 specifically includes: The single-scale inversion of the tunnel elastic wave full waveform inversion method based on total variation regularization constraints is performed according to the following equation: Where: E(m) is the objective function; λ is the weight coefficient of the regularization term; m represents the P-wave and S-wave velocity models and the density model; u cal and u obs Respectively represent the simulated and observed particle displacement data; m k represents the kth iteration model parameter, α k-1 represents the step size, δm k-1 represents the search direction, i.e. the model update amount; L(p N ,q N ) is the total variation regularization constraint update term, p N ,q N is the regularization constraint operator.

5. A tunnel full waveform inversion method based on cutoff frequency optimization and regularization constraint according to claim 4, characterized in that: The single-scale inversion of the tunnel elastic wave full waveform inversion method based on total variation regularization constraint in S3 specifically includes: S301: Calculate the objective function value according to E(m); S302: Using the objective function value obtained in S301, use the parabola fitting method and the L-BFGS method to calculate the step length α k-1 and search direction δm k-1 ; S303: Based on step size α k-1 and search direction δm k-1 Calculate the total variation regularization constraint update term L(p N ,q N ); S304: According to formula m k =(m k-1 +α k-1 δm k-1 )-λ·L(p N ,q N ) Iteratively update the model, λ is the weight coefficient of the regularization term, m k represents the kth iteration model parameter, m k-1 represents the k-1th iteration model parameter, and the single-scale inversion result of the tunnel elastic wave full waveform inversion method based on total variation regularization constraint is obtained.

6. The tunnel full waveform inversion method based on cutoff frequency optimization and regularization constraint according to claim 1 is characterized in that: The S4 specifically includes: S401: The single-scale inversion result obtained in S3 is subjected to multi-parameter weighted constraint structure correction processing according to the following equation, and the result after multi-parameter weighted constraint structure correction is obtained: Where: Δm i,j p is the single iteration model update amount; Δm′ i,j p is the model update amount after multi-parameter weighted structure correction; m i,j p is the multi-parameter initial model; p = V p ,V s , Den represents the medium longitudinal wave velocity, shear wave velocity, and density physical parameters, nx is the number of horizontal grids in the model, nz is the number of vertical grids in the model, and i, j represent the positions of grid nodes in the z and x directions respectively; S402: Perform one-dimensional velocity profile spatial structure correction and smoothing constraint on the result of multi-parameter weighted constraint structure correction obtained in S401.

7. A tunnel full waveform inversion method based on cutoff frequency optimization and regularization constraint according to claim 6, characterized in that: The S402 specifically includes: S4021: Set the restriction conditions of the spatial structure correction of the one-dimensional wave velocity profile, including the lower limit value of the spatial correction, the range of the structural correction area, and the distance between the areas; S4022: Select a vertical grid coordinate y i , extract the one-dimensional wave velocity profile along the tunnel axis Among them, x is the horizontal axis coordinate; S4023: Calculate one-dimensional wave velocity profile The corresponding model update amount With slope S4024: 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 S4025: Update the amount based on the corrected model 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. S4026: Update the amount based on the corrected model 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 S4027: Update the amount based on the corrected model Calculate the corrected one-dimensional velocity profile S4028: Take the next vertical grid coordinate y i+1 , repeat S4023 to S4027 until the correction of the one-dimensional velocity profiles of all grids is completed, and the result after the velocity structure correction is obtained.

8. A tunnel full waveform inversion method based on cutoff frequency optimization and regularization constraint according to claim 7, characterized in that: The S4026 specifically 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.

9. A computer readable medium storing instructions, characterized in that: When the instructions are executed on the readable medium, the readable medium is caused to execute the method according to any one of claims 1 to 8.

10. An electronic device comprising a memory and a processor, wherein the memory stores a computer program that can be run on the processor, wherein: When the processor executes the computer program, the method according to any one of claims 1 to 8 is implemented.

Images