An Inversion Method for S-Wave Velocity Structure Based on Dense Micromotion Observations

Abstract

This invention discloses an inversion method for S-wave velocity structure based on dense micromotion observations, belonging to the field of geophysical exploration and engineering geological exploration technology. The method includes: acquiring micromotion data; performing segmented preprocessing and cross-correlation calculations; constructing a phase consistency metric function based on the instantaneous phase information of the cross-correlation function for each time period; generating adaptive weighting coefficients; and then constructing an empirical Green's function through weighted superposition. Subsequently, time-frequency analysis is performed to obtain a time-frequency energy map; preliminary extraction of candidate dispersion velocities for each frequency is performed, and their comprehensive confidence level is calculated; after verification, a reliable dispersion curve is obtained. Based on the dispersion curve, a model parameter vector is established, and a joint objective function is constructed. A differential evolution algorithm is used for parallel global optimization inversion to obtain a one-dimensional S-wave velocity structure model below each station pair. The one-dimensional models are interpolated and fused to generate a two-dimensional S-wave velocity structure profile. This invention has the advantages of high inversion stability, high computational efficiency, and a high degree of automation.

Description

Technical Field

[0001] This invention belongs to the field of geophysical exploration and engineering geological exploration technology, and more specifically, relates to an inversion method for S-wave velocity structure based on dense micromotion observation. Background Technology

[0002] Micromotion detection is a passive source geophysical method that utilizes continuous vibration signals (environmental noise) generated by natural and human activities to detect underground structures. This method requires no artificial seismic source and has advantages such as low cost, convenient construction, and minimal environmental interference, making it particularly suitable for areas where active seismic sources are limited, such as cities, mining areas, and cultural heritage protection zones. Its core principle is to extract the dispersion characteristics (i.e., the relationship between phase velocity or group velocity and frequency) of the dominant surface wave (such as Rayleigh wave) from the micromotion signal, and then invert the shear wave velocity structure of the underground medium, thereby revealing geological targets such as stratigraphic interfaces, lithological changes, cavities, and mining subsidence areas.

[0003] With the development of sensor technology and data acquisition capabilities, long-term micro-motion observation using dense arrays has become an important means of obtaining high-resolution information on underground structures. However, the influx of massive amounts of data poses a severe challenge to traditional micro-motion data processing workflows.

[0004] A search revealed patent publication number CN120122178A, published on June 10, 2025, which discloses a data processing method for micro-motion detection array deployment. This method generates field stations in batches using RTK measured coordinates and uses the DBSCAN clustering algorithm to spatially average the spatial autocorrelation coefficients of station pairs, thereby picking up dispersion data and performing shear wave velocity inversion. While this type of method improves the efficiency of field station deployment and optimizes the calculation of autocorrelation coefficients through cluster analysis, it still mainly relies on the central station model at the algorithm level. Failure of the central equipment may affect the overall acquisition. Its inversion stage often employs an adaptive gradient strategy, which is highly dependent on the initial model. Furthermore, when processing large-scale dense array data, it is prone to getting trapped in local minima, affecting the stability of the inversion results.

[0005] Patent publication number CN111164462A, published on May 15, 2020, discloses an artificial source surface wave exploration method. This method obtains underground structural information by artificially generating surface waves and performing multi-channel acquisition and dispersion analysis. This method is a typical active source surface wave detection technology. While it has some application value in engineering surveys, its data acquisition and processing relies on artificial seismic sources and fixed survey lines, making it unsuitable for urban environments or large-scale, long-term continuous monitoring scenarios. Furthermore, it lacks optimization design for the automatic extraction and inversion of surface wave dispersion under conditions of subtle noise.

[0006] Meanwhile, existing methods for automatically extracting background noise dispersion curves, such as patent publication number CN112861721A, published on May 28, 2021, disclose a method for automatically extracting background noise dispersion curves, including: acquiring background noise data recorded by seismic stations, preprocessing the background noise data, cross-correlating and superimposing the background noise data to obtain an empirical Green's function, and then performing a vector wavenumber domain transformation on the empirical Green's function to obtain a dispersion spectrum; obtaining candidate dispersion regions based on the dispersion spectrum and the Res-Unet++ network model; and extracting and classifying the dispersion points of the candidate dispersion regions using the gradient method and the chasing method to obtain the target dispersion curve. This method focuses primarily on segmenting and extracting the dispersion spectrum through neural networks, while traditional dispersion extraction criteria often rely solely on single indicators such as the maximum energy point or peak strength of the time-frequency graph. Traditional single-indicator discrimination methods are prone to mispicking in low signal-to-noise ratio regions, leading to non-physical jumps in the dispersion curve. In addition, traditional inversion objective functions typically fit only a single phase velocity and use equal-weighted errors, resulting in a high degree of ambiguity in the inversion results.

[0007] With the rapid development of dense array deployment technology and long-term continuous recording capabilities, the acquired micro-motion data has seen significant improvements in both spatial coverage and temporal length, providing a data foundation for high-resolution, large-scale detailed imaging of underground structures. However, this also brings new technical challenges: 1) The data volume is huge, and traditional manual or semi-automatic dispersion extraction methods are inefficient and cannot meet the needs of engineering practice for rapid processing; 2) Significant differences in signal quality exist between different stations and at different times, necessitating an effective automated quality control mechanism to filter reliable data; 3) Faced with thousands of station pairs, traditional inversion methods are computationally time-consuming, have strong dependence on the initial model, and lack stability and automation of the inversion results.

[0008] Therefore, there is an urgent need for an inversion method for S-wave velocity structure based on dense micro-motion observation, which can automatically and efficiently process micro-motion data and stably and reliably invert the underground S-wave velocity structure. Summary of the Invention

[0009] 1. The problem to be solved The purpose of this invention is to provide an inversion method for S-wave velocity structure based on dense micro-motion observations, aiming to address the shortcomings of existing technologies in terms of automation, processing efficiency, result stability, and end-to-end integration. Specifically, it utilizes massive, multi-source, low signal-to-noise ratio raw data, and constructs a high-quality empirical Green's function (EGF) through an innovative and adaptive signal quality assessment and fusion method. Then, it extracts key features with high confidence (dispersion curves), and finally obtains the highly reliable target structure (S-wave velocity) through multi-constraint stable inversion.

[0010] 2. Technical Solution To solve the above problems, the technical solution adopted by the present invention is as follows: This invention provides an inversion method for S-wave velocity structure based on dense micromotion observations, comprising the following steps: S1. Constructing the empirical Green's function: S11. Acquire micro-motion data collected by a dense array of stations deployed in the detection area; the micro-motion data is a long time series. S12. Perform segmented preprocessing on the micro-motion data, and calculate the cross-correlation function for the signals of each station pair in each time period; S13. Based on the instantaneous phase information of the cross-correlation function of each time period, construct a phase consistency measurement function and generate adaptive weighting coefficients accordingly. S14. Using the weighting coefficients, the cross-correlation functions of each time period are weighted and superimposed to construct an empirical Green's function; S2. Extract the dispersion curve: S21. Perform time-frequency analysis on the empirical Green's function obtained in step S14 to obtain the time-frequency energy map; S22. Based on the time-frequency energy map, preliminary extraction of candidate values ​​for dispersion velocities corresponding to each frequency is performed, and the comprehensive confidence level of each dispersion point is calculated; the candidate values ​​for dispersion velocities include at least candidate values ​​for phase velocity and candidate values ​​for group velocity. S23. Perform continuity verification on the initially extracted dispersion curve, remove outliers that do not meet the continuity requirements, and obtain a reliable dispersion curve. S3 and S-wave velocity structure batch inversion: S31. Based on the dispersion curves of multiple station pairs extracted in step S23, establish the layered medium model parameter vector and construct a joint inversion objective function that integrates the relative error between observed and theoretical dispersion curves, the comprehensive confidence weight, and the model smoothing constraint. S32. The objective function is optimized nonlinearly and globally using the differential evolution algorithm, and a one-dimensional S-wave velocity structure model below each station is obtained by parallel inversion. S4. Velocity profile construction: The one-dimensional S-wave velocity structure models of all station pairs obtained in step S32 are interpolated and fused according to their spatial positions to generate a two-dimensional S-wave velocity structure profile of the detection area.

[0011] As one possible implementation, in step S12, the segmented preprocessing includes: performing detrending, mean removal, and bandpass filtering on the original micro-motion data of each station, and performing RMS normalization segment by segment.

[0012] As one possible implementation, in step S13, the method for constructing the phase consistency metric function is as follows: The continuously acquired micro-motion signals are divided into N time periods. Cross-correlation is performed on the micro-motion signals from station i and station j within each time period to obtain the corresponding cross-correlation function. Where τ represents the time delay, and n = 1, 2, ..., N; Cross-correlation function for each time period Perform analytical signal transformation to obtain its corresponding instantaneous phase information. ; At the same time delay τ, the instantaneous phase of all time periods is statistically analyzed to construct a phase consistency metric function. :

[0013] Where k is the sensitivity adjustment parameter, with a value range of [1,2]; reference phase It is the vector average of the instantaneous phases over all time periods; In the above formula: Reference phase , where i is the imaginary unit, calculates the statistically significant average phase by vector averaging the instantaneous phases of all N time periods; it provides an objective dynamic benchmark for measuring the dispersion of the phase in each time period, and can effectively reflect the overall phase characteristics of the station at a specific delay.

[0014] Quantification of phase deviation To measure the difference between the phase of the nth time period and the average phase, the periodicity of the cosine function is used to handle the circular statistical problem of the phase, accurately capturing the consistency of signal fluctuations.

[0015] Exponential mapping and sensitivity adjustment: A negative exponential function maps the mean deviation to the [0,1] interval, and a parameter k is introduced to adjust the sensitivity, ranging from [1,2]. A larger value indicates greater sensitivity to the penalty for phase inconsistency. This transforms the abstract statistical characteristics of phase into an intuitive reliability score; higher consistency results in a higher score. The closer a score is to 1, the worse the consistency, and the score rapidly decays towards 0.

[0016] As one possible implementation, in step S13, the method for generating the adaptive weighting coefficients is as follows:

[0017] in: Phase consistency measurement function The mean of the neighborhood around the delay τ.

[0018] Let be the phase standard deviation at delay τ for all time periods.

[0019] The instantaneous phase at delay τ in the nth time interval is compared with the reference phase. The absolute deviation.

[0020] α is a nonlinear enhancement parameter, ranging from [1,4]; it acts like a magnifying glass, amplifying the weight difference between high-quality signals and low-quality noise. If the overall coherence is poor at this time delay point, the weight will decrease rapidly.

[0021] β is the phase dispersion adjustment parameter, with a value range of [0.5, 2.0]; it controls the degree to which dispersion reduces the weight. The larger the dispersion (the more unstable the signal), the smaller the exponential term, thus effectively suppressing signal segments with severe transient interference.

[0022] γ is a phase deviation adjustment parameter for a single time period, with a value range of [0.1, 1.0]; it allows the formula to treat specific time periods individually. Even if the overall quality of the delay point is good, if the phase of a specific time period n deviates from the reference phase, this term will increase the denominator, thereby fine-tuning and reducing the weight contribution of that specific time period.

[0023] ε is the stability control parameter, taking values ​​

[10] . -6 10 -3 This prevents numerical calculation anomalies in cases of extremely high coherence (denominator tending to 0), ensuring the mathematical robustness of large-scale, long-term data during automated processing.

[0024] As one possible implementation, in step S14, the weighting coefficients are used to perform a weighted superposition of the cross-correlation functions of each time period to obtain the empirical Green's function between station i and station j. :

[0025] Where ψ(·) is a nonlinear magnitude mapping function used to adjust the magnitude of the cross-correlation function; satisfying: ;in Represents the cross-correlation function The mean.

[0026] The phase coupling term, used to characterize the phase consistency across time intervals, is defined as follows: .

[0027] Traditional phase-weighted superposition (PWS) typically uses only one phase coherence factor:

[0028] Empirical Green's Function of this Invention The method differs significantly from traditional formulas: first, the instantaneous phase of each cross-correlation function is extracted; then, the phase consistency across all time periods at the same delay is statistically analyzed; further, a phase dispersion adjustment parameter β, a single-time-period phase deviation adjustment parameter γ, and a stability control parameter ε are added to the weights; finally, a nonlinear amplitude mapping function ψ( ) and phase coupling terms.

[0029] As one possible implementation, in step S22, the comprehensive confidence level... q Calculated using the following formula:

[0030] in, Indicates frequency; This represents the maximum peak amplitude of the cross-correlation function at this frequency.

[0031] All frequency points The statistical average is used for amplitude normalization.

[0032] This represents the time-frequency energy value corresponding to this frequency; the time-frequency energy at this frequency point obtained through FTAN time-frequency analysis.

[0033] This represents the maximum energy value in the time-frequency energy graph.

[0034] The standard deviation of the cross-correlation peak values ​​at different time periods at this frequency reflects the stability of the signal.

[0035] For all frequency points The statistical mean and the average standard deviation are used for the normalization of stability indicators.

[0036] , The feature weight coefficients satisfy the following conditions: .

[0037] η is a nonlinear mapping adjustment parameter that controls the sensitivity of confidence evaluation. A higher value will make the confidence score more polarized, that is, make high quality points closer to 1 and low quality points more quickly approach 0.

[0038] The overall confidence level As a frequency-related weight, it directly participates in the construction of the joint inversion objective function, which is used to enhance the constraint effect of high-reliability frequency scatter points on the inversion results and suppress the influence of low-reliability frequency scatter points.

[0039] Traditional dispersion extraction often relies on two simple criteria: selecting only the maximum energy point on the time-frequency graph, or using only peak strength as the screening standard. These methods are essentially single-index discrimination. However, the comprehensive confidence q of this invention integrates three types of information: the maximum peak amplitude of cross-correlation; the standard deviation of peak values ​​in different time periods; and the time-frequency energy of FTAN.

[0040] As one possible implementation, in step S23, the rule for continuity verification is: a frequency point is considered a reliable point and is retained only if the dispersion curve is effectively extracted on at least M consecutive frequency points, where M is a preset continuity threshold.

[0041] As one possible implementation, in step S23, the dispersion curve is subjected to continuity constraint screening. The dispersion curve is only allowed to appear continuously at at least A adjacent frequency points, where A is a preset continuity threshold, usually 3-7. The larger the value, the stronger the constraint.

[0042] As one possible implementation, in step S31, the layered medium model parameter vector m is:

[0043] in, This represents the S-wave velocity of the d-th layer. This represents the thickness of the d-th layer, where d is the total number of layers in the model.

[0044] As one possible implementation, in step S31, the joint inversion objective function for:

[0045] in, This represents the k-th discrete frequency point.

[0046] q denoted as the overall confidence level of the k-th discrete frequency point.

[0047] and These are the observed phase velocity and group velocity, respectively.

[0048] and These are the theoretical phase velocity and group velocity calculated from model m.

[0049] L is the model smoothing matrix, which is a second-order difference matrix used to quantify the changes in velocity and thickness between adjacent strata.

[0050] λ is a regularization parameter used to balance data fitting error and model smoothness. Adjusting this value can prevent overfitting or non-physical solutions from being generated during inversion. When λ is large, the inversion results tend to be smoother, which helps to identify large-scale layered structures. The λ inversion results will fit the observation data better, but may introduce more unstable fluctuations.

[0051] The penalty model changes drastically to prevent overfitting, noise, and non-physical solutions.

[0052] Common inversion objective functions are often written as a single error term, fitting only the phase velocity, for example:

[0053] In the formula v obs For the observed phase velocity, v cal This is the theoretical phase velocity obtained through calculation.

[0054] However, the joint inversion objective function of this invention In the middle: the fitting of single phase velocity was changed to joint fitting of phase velocity and group velocity; the ordinary error was changed to relative error; the equal weighted error was changed to the introduction of comprehensive confidence as frequency-related weight; and a model smoothing regularization term was added.

[0055] As one possible implementation, in step S32, the dispersion curves of multiple station pairs are inverted in parallel using the differential evolution algorithm. The objective function of each station pair is solved independently, and finally a one-dimensional S-wave velocity structure model below the path of each station pair is obtained.

[0056] 3. Beneficial effects Compared with the prior art, the beneficial effects of the present invention are as follows: (1) This invention proposes a comprehensive confidence index that integrates cross-correlation peak amplitude, peak stability, and time-frequency energy, and constructs a multi-objective objective function with joint constraints of phase velocity and group velocity, confidence weighting, and smoothing regularization. This multi-dimensional quality control and joint constraint mechanism effectively suppresses overfitting noise and significantly improves the reliability of inversion results in the detection of complex underground structures.

[0057] (2) The present invention uses the differential evolution algorithm for global optimization, which significantly reduces the dependence on the initial model and realizes batch parallel inversion of multiple stations.

[0058] (3) This invention achieves batch extraction and inversion of dispersion curves of large-scale stations by using parallel computing strategy, fast dispersion extraction algorithm (such as short-time Fourier transform) and differential evolution algorithm suitable for global search, which greatly improves data processing efficiency.

[0059] (4) In the construction of the empirical Green's function, this invention introduces a multi-parameter adaptive weighting mechanism based on instantaneous phase consistency, phase dispersion, and single-time phase deviation, which significantly improves the signal-to-noise ratio. In the dispersion extraction, a comprehensive confidence index integrating cross-correlation peak amplitude, peak stability, and time-frequency energy is proposed, and combined with continuity verification, automatic quality control is achieved. In the inversion, a multi-objective function with joint constraints of phase velocity and group velocity, confidence weighting, and smoothing regularization is constructed, which effectively suppresses noise interference, avoids the inversion from falling into local extrema, and improves the stability of the results.

[0060] (5) The entire process of this invention is designed as an automated production line, minimizing human intervention. It is particularly suitable for processing massive micro-motion data collected by dense arrays, and can quickly generate intuitive two-dimensional velocity profiles. It has a clear ability to identify low-velocity anomalies such as goaf and karst caves, and has strong engineering interpretability.

[0061] (6) The method framework of the present invention is clear, and each module is relatively independent. The parameters (such as filter frequency band, number of inversion layers, regularization intensity, etc.) can be adjusted according to the actual data quality and detection requirements. It is easy to expand and integrate into different micro-motion data processing systems. Attached Figure Description

[0062] Figure 1 This is a flowchart of the overall process for the inversion method of S-wave velocity structure based on dense micro-motion observations according to the present invention.

[0063] Figure 2 This is a schematic diagram illustrating the process of constructing the empirical Green's function and extracting the dispersion curve in this invention.

[0064] Figure 3 This is a schematic diagram of the joint constraint inversion model.

[0065] Figure 4 This is a schematic diagram illustrating the enhancement of the dispersion curve in the embodiment.

[0066] Figure 5 The diagram shows the dispersion inversion structure of a certain station in the embodiment.

[0067] Figure 6 This is a two-dimensional S-wave velocity profile of measuring line a in the embodiment.

[0068] Figure 7 This is a two-dimensional S-wave velocity profile of measuring line b in the embodiment.

[0069] Figure 8 This is a schematic diagram of the borehole depth in a goaf area, as shown in the example. Detailed Implementation

[0070] The present invention will be further described below with reference to specific embodiments.

[0071] Example 1: Application of goaf detection in a mining area Using the detection of a goaf in a certain XX as an application scenario, this embodiment demonstrates the specific implementation process of the present invention.

[0072] The overall processing flow is as follows Figure 1 As shown, the method comprises four main steps: constructing an empirical Green's function, extracting dispersion curves, batch inversion of S-wave velocity structure, and constructing velocity profiles. These steps are interconnected via data streams to form a complete automated processing flow.

[0073] Specifically: S1. Constructing the empirical Green's function: S11. Acquire micro-motion data collected by a dense array of stations deployed in the detection area: 247 densely arranged micro-motion geophones were deployed along the pre-defined survey line a, and 258 stations were deployed along survey line b. Continuous observations were conducted for 3 hours at each station, with a sampling interval of 4 ms. The raw data recorded at each station were detrended and mean-removed to eliminate baseline drift. Bandpass filtering of 0.1–10 Hz was applied to preserve the effective frequency band of the surface waves.

[0074] S12. Perform segmented preprocessing on the micro-motion data, and calculate the cross-correlation function for the signals of each station pair within each time period: The long record was divided into N=36 5-minute intervals, and RMS amplitude normalization was performed on each interval to suppress strong transient interferences such as earthquakes and vehicle traffic. For all possible station pairs (approximately 30,000 pairs), the cross-correlation function within each time interval was calculated. .

[0075] S13. Based on the instantaneous phase information of the cross-correlation function of each time period, construct a phase consistency measurement function and generate adaptive weighting coefficients accordingly: Instantaneous phase extracted by analytical signal transformation And calculate the baseline average phase for the entire time period. Construct a phase consistency metric function The sensitivity adjustment parameter k is set to 1.5, such as... Figure 2The diagram illustrates the process from segmenting the original micro-motion data, calculating cross-correlation, extracting the instantaneous phase, to weighted superposition to generate the Empirical Green's Function (EGF). It also demonstrates time-frequency analysis based on EGF, calculation of comprehensive confidence levels, and extraction of reliable dispersion curves through continuity verification.

[0076] Where k is the sensitivity adjustment parameter, with a value range of [1,2]; reference phase It is the vector average of the instantaneous phases over all time periods; In the above formula: Reference phase , where i is the imaginary unit, calculates the statistically significant average phase by vector averaging the instantaneous phases of all N time periods; it provides an objective dynamic benchmark for measuring the dispersion of the phase in each time period, and can effectively reflect the overall phase characteristics of the station at a specific delay.

[0077] Quantification of phase deviation To measure the difference between the phase of the nth time period and the average phase, the periodicity of the cosine function is used to handle the circular statistical problem of the phase, accurately capturing the consistency of signal fluctuations.

[0078] Exponential mapping and sensitivity adjustment: A negative exponential function maps the mean deviation to the [0,1] interval, and a parameter k is introduced to adjust the sensitivity, ranging from [1,2]. A larger value indicates greater sensitivity to the penalty for phase inconsistency. This transforms the abstract statistical characteristics of phase into an intuitive reliability score; higher consistency results in a higher score. The closer a score is to 1, the worse the consistency, and the score rapidly decays towards 0.

[0079] based on Construct the weighting coefficients for each time period:

[0080] in, Phase consistency measurement function The mean of the neighborhood around the delay τ.

[0081] Let be the phase standard deviation at delay τ for all time periods.

[0082] The instantaneous phase at delay τ in the nth time interval is compared with the reference phase. The absolute deviation.

[0083] Nonlinear enhancement parameters Phase dispersion adjustment parameter Single-time phase deviation adjustment parameters Stable control parameters .

[0084] S14. Using the weighting coefficients, the cross-correlation functions of each time period are weighted and superimposed to construct an empirical Green's function: finally, an empirical Green's function with a high signal-to-noise ratio is generated by superposition. :

[0085] Where ψ(·) is a nonlinear magnitude mapping function used to adjust the magnitude of the cross-correlation function; satisfying: ;in Represents the cross-correlation function The mean.

[0086] The phase coupling term, used to characterize the phase consistency across time intervals, is defined as follows: .

[0087] S2. Extract the dispersion curve: S21. Perform time-frequency analysis on the empirical Green's function obtained in step S14, that is, perform short-time Fourier transform on the EGF of each station pair to generate a time-frequency energy map.

[0088] S22. Based on the time-frequency energy map, preliminary candidate values ​​for dispersion velocities corresponding to each frequency are extracted, and the overall confidence level of each dispersion point is calculated: by tracing the path of maximum energy, the phase velocity dispersion curve is preliminarily calculated. The overall confidence level of each station for the EGF is calculated.

[0089] Indicates frequency; This represents the maximum peak amplitude of the cross-correlation function at this frequency.

[0090] All frequency points The statistical average is used for amplitude normalization.

[0091] This represents the time-frequency energy value corresponding to this frequency; the time-frequency energy at this frequency point obtained through FTAN time-frequency analysis.

[0092] This represents the maximum energy value in the time-frequency energy graph.

[0093] In this embodiment, the weight is set as follows: =0.4, =0.4, =0.2, adjust parameter The value is 8.

[0094] reserve q A value of >=0.80 is used for subsequent inversion weighting.

[0095] S23. Perform continuity verification on the initially extracted dispersion curve (in this embodiment, the effective points are continuous across 5 adjacent frequencies), eliminate abnormal jump points, and finally obtain a smooth and reliable dispersion curve, such as... Figure 4 As shown.

[0096] S3 and S-wave velocity structure batch inversion: S31. Based on the dispersion curves of multiple station pairs extracted in step S23, establish the layered medium model parameter vector m: Set up a 5-layer horizontal layered initial model space (the speed and thickness range are set according to the geological data of the work area).

[0097]

[0098] Table 1 is a table of recommended initial model space settings for 5 layers for a specific mining area detection scenario (such as a goaf) in the example: Subsequently, a joint inversion objective function was constructed that integrates the relative error between observed and theoretical dispersion curves, the comprehensive confidence weight, and the model smoothing constraint:

[0099] in, This represents the k-th discrete frequency point.

[0100] q denoted as the overall confidence level of the k-th discrete frequency point.

[0101] and These are the observed phase velocity and group velocity, respectively.

[0102] and These are the theoretical phase velocity and group velocity calculated from model m.

[0103] L is the model smoothing matrix, which is a second-order difference matrix used to quantify the changes in velocity and thickness between adjacent strata.

[0104] λ is a regularization parameter, set to 0.5, used to balance data fitting error and model smoothness. Adjusting this value prevents overfitting or non-physical solutions from being generated during inversion. When λ is larger, the inversion results tend to be smoother, which helps to identify large-scale layered structures. The λ inversion results will fit the observation data better, but may introduce more unstable fluctuations.

[0105] The penalty model changes drastically to prevent overfitting, noise, and non-physical solutions.

[0106] Table 1. Suggested initial model space settings for 5 layers in a mining area detection scenario (such as a goaf) in this embodiment.

[0107] S32. The differential evolution algorithm is used to perform nonlinear global optimization on the above objective function, and the one-dimensional S-wave velocity structure model below each station is obtained by parallel inversion: Leveraging the parallel capabilities of the computing cluster, the dispersion curves of hundreds of station pairs are simultaneously inverted, automatically outputting a one-dimensional Vs model for each path. For example... Figure 3 As shown, the model focuses on the relative error between observed dispersion data (phase velocity and group velocity) and theoretical dispersion data. A comprehensive confidence level is introduced as a frequency-related weight, and a model smoothing regularization term is added. Global optimization is performed using a differential evolution algorithm to output a one-dimensional S-wave velocity structure model. An example using a set of stations is shown below. Figure 5 As shown.

[0108] S4. Velocity Profile Construction: The one-dimensional model obtained from the inversion is fused according to spatial location using Kriging interpolation to generate two-dimensional S-wave velocity profiles along survey lines a and b. (See details below.) Figure 6 , Figure 7 As shown.

[0109] Drilling data is direct, real-time geological information obtained from field boreholes. It has extremely high vertical resolution and stratigraphic boundary accuracy, mainly including station coordinates, elevation coordinates, and physical conditions such as lithology, hardness, groundwater distribution, and the presence of mined-out areas at different depths. In the exploration process, it is regarded as the surface truth and is used to calibrate inversion models and verify the accuracy of geophysical inferences. A schematic diagram of some borehole depths in a mined-out area of ​​a certain XX is shown in this example. Figure 8 As shown.

[0110] The two-dimensional S-wave velocity profile of survey line a clearly shows multiple transversely continuous low-velocity anomaly zones, with wave velocities mainly concentrated between 600-1100 m / s. The spatial distribution of the low-velocity anomaly zones obtained by inversion is consistent with... Figure 8 The location of the goaf in the known two-dimensional boreholes is highly consistent with that of the goaf, proving that the method has extremely high identification accuracy for underground anomalies such as goafs.

[0111] The above description is merely an embodiment of this specification and does not limit the patent scope of this invention. Various modifications and variations can be made to this specification by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this specification should be included within the scope of the claims of this specification.

Claims

1. An inversion method for S-wave velocity structure based on dense micromotion observations, characterized in that: Includes the following steps: S1. Constructing the empirical Green's function: Acquire micro-motion data collected by a dense array of stations deployed in the detection area; perform segmented preprocessing on the micro-motion data, and calculate the cross-correlation function for the signals of each station pair in each time period; construct a phase consistency measurement function based on the instantaneous phase information of the cross-correlation function in each time period, and generate adaptive weighting coefficients accordingly; use the weighting coefficients to perform weighted superposition of the cross-correlation functions in each time period to construct the empirical Green's function. S2. Extracting dispersion curves: Perform time-frequency analysis on the empirical Green's function obtained in step S1 to obtain a time-frequency energy map; based on the time-frequency energy map, initially extract candidate values ​​of dispersion velocity corresponding to each frequency, and calculate the comprehensive confidence level of each dispersion point; perform continuity verification on the initially extracted dispersion curves, remove outliers that do not meet the continuity requirements, and obtain reliable dispersion curves. S3, Batch Inversion of S-wave Velocity Structure: Based on the dispersion curves of multiple station pairs extracted in step S2, a layered medium model parameter vector is established, and a joint inversion objective function is constructed that integrates the relative error between observed and theoretical dispersion curves, the comprehensive confidence weight, and the model smoothing constraint. The objective function is nonlinearly optimized globally using the differential evolution algorithm, and a one-dimensional S-wave velocity structure model below each station pair is obtained through parallel inversion. S4. Velocity profile construction: The one-dimensional S-wave velocity structure models of all station pairs obtained from step S3 are interpolated and fused according to their spatial positions to generate a two-dimensional S-wave velocity structure profile of the detection area.

2. The inversion method for S-wave velocity structure based on dense micro-motion observations according to claim 1, characterized in that: In step S1, the segmented preprocessing includes: performing detrending, mean removal, and bandpass filtering on the raw micro-motion data of each station, and performing RMS normalization segment by segment.

3. The inversion method for S-wave velocity structure based on dense micro-motion observations according to claim 1, characterized in that: In step S1, the method for constructing the phase consistency metric function is as follows: The continuously acquired micro-motion signals are divided into N time periods. Cross-correlation is performed on the micro-motion signals from station i and station j within each time period to obtain the corresponding cross-correlation function. Where τ represents the time delay, and n = 1, 2, ..., N; Cross-correlation function for each time period Perform analytical signal transformation to obtain its corresponding instantaneous phase information. ; At the same time delay τ, the instantaneous phase of all time periods is statistically analyzed to construct a phase consistency metric function. : Where k is the sensitivity adjustment parameter, with a value range of [1,2]; reference phase It is the vector average of the instantaneous phase over all time periods.

4. The inversion method for S-wave velocity structure based on dense micro-motion observations according to claim 3, characterized in that: In step S1, the method for generating the adaptive weighting coefficients is as follows: in, Phase consistency measurement function The mean of the neighborhood around the delay τ; The phase standard deviation at delay τ for all time periods; The instantaneous phase at delay τ in the nth time interval is compared with the reference phase. The absolute deviation; α is the nonlinear enhancement parameter, and its value ranges from [1,4]. β is the phase dispersion adjustment parameter, with a value range of [0.5, 2.0]. γ is a single-time phase deviation adjustment parameter, with a value range of [0.1, 1.0]. ε is the stability control parameter, taking values ​​[10]. -6 10 -3 ].

5. The inversion method for S-wave velocity structure based on dense micro-motion observations according to claim 4, characterized in that: In step S1, the empirical Green's function between station i and station j for: Wherein, ψ(·) is a nonlinear amplitude mapping function used to adjust the amplitude of the cross-correlation function.

6. The inversion method for S-wave velocity structure based on dense micro-motion observations according to claim 1, characterized in that: In step S2, the comprehensive confidence level q Calculated using the following formula: in, Indicates frequency; This represents the maximum peak amplitude of the cross-correlation function at this frequency; All frequency points The statistical average; This is the time-frequency energy value corresponding to that frequency; This represents the maximum energy value in the time-frequency energy graph. This represents the standard deviation of the peak cross-correlation values ​​at different time intervals at this frequency. For all frequency points The statistical average; , The feature weight coefficients satisfy the following conditions: ; η is a nonlinear mapping adjustment parameter that controls the sensitivity of confidence evaluation.

7. The inversion method for S-wave velocity structure based on dense micro-motion observations according to claim 6, characterized in that: In step S2, the rule for continuity verification is: a frequency point is considered a reliable point and is retained only if the dispersion curve is effectively extracted on at least M consecutive frequency points, where M is a preset continuity threshold.

8. The inversion method for S-wave velocity structure based on dense micro-motion observations according to claim 1, characterized in that: In step S3, the layered medium model parameter vector m is in, This represents the S-wave velocity of the d-th layer. This represents the thickness of the d-th layer, where d is the total number of layers in the model.

9. The inversion method for S-wave velocity structure based on dense micro-motion observations according to claim 8, characterized in that: In step S3, the joint inversion objective function for: in, This represents the k-th discrete frequency point; q The overall confidence level of the k-th discrete frequency point; and These are the observed phase velocity and group velocity, respectively. and These are the theoretical phase velocity and group velocity calculated from model m; L is the model smoothing matrix; λ is a regularization parameter used to balance data fitting error and model smoothness.

10. The inversion method for S-wave velocity structure based on dense micro-motion observations according to claim 9, characterized in that: In step S3, the differential evolution algorithm is used to invert the dispersion curves of multiple station pairs in parallel. The objective function of each station pair is solved independently, and finally a one-dimensional S-wave velocity structure model below the path of each station pair is obtained.

Images