Karst conduit three-dimensional trajectory tracking method and system based on array micro-motion signal interference imaging

Abstract

The application discloses a karst pipeline three-dimensional trajectory tracking method and system based on array microseismic signal interference imaging. The method collects environmental microseismic signals by arranging multiple microseismic detector sites on the ground; continuous microseismic waveform data are acquired and a cross-power spectral density matrix is constructed to extract characteristic microseismic dispersion information representing underground fluid pulsation; then, a three-dimensional interference focusing imaging algorithm is used to calculate the superposition strength of the seismic source signal and construct a three-dimensional energy distribution model; finally, a strip-shaped high-energy anomaly area is identified, and the continuous three-dimensional spatial trajectory and spatial positioning reliability distribution of the karst pipeline are reconstructed. Compared with the prior art, the application greatly improves the signal-to-noise ratio of target identification by using the unique fluid pulsation characteristics, and realizes high-precision and quantifiable evaluation of the true three-dimensional trajectory reconstruction of the strike, dip angle and burial depth of the medium-deep karst pipeline.

Description

Technical Field

[0001] This invention belongs to the field of engineering surveying technology, specifically relating to a method and system for three-dimensional trajectory tracking of karst pipelines based on array-type micro-motion signal interferometric imaging. Background Technology

[0002] Karst regions often feature underground systems of caves, fissures, and conduits of varying sizes. Among these, karst conduits, characterized by clear water flow channels, are the most important hydrogeological structural units in karst areas and a major potential source of geological hazards such as ground subsidence and uneven foundation settlement. Therefore, identifying the three-dimensional spatial distribution patterns of karst conduits is of paramount engineering value for project planning and site selection, groundwater resource assessment, and geological hazard prevention.

[0003] Currently, the main methods for detecting karst conduits include:

[0004] Currently used geophysical exploration methods include high-density direct current (DC) resistivity tomography (HDT), transient electromagnetic methods (TEM), and ground-penetrating radar (GPR). Although GPR has high resolution, the high-frequency electromagnetic waves experience severe polarization attenuation in water-bearing clay overburden (red clay commonly found in karst areas), and the penetration depth is usually no more than 20 meters, making it unable to cover underground rivers buried more than 100 meters deep. On the other hand, the spatial resolution of DC resistivity tomography and transient electromagnetic methods decreases exponentially with increasing detection depth, often only detecting a large, blurry area of ​​low resistivity anomaly, and completely failing to characterize the specific three-dimensional pipe morphology of underground rivers.

[0005] Other methods include micromotion exploration, which uses natural environmental vibrations as a passive seismic source and utilizes the surface wave component of environmental noise to invert underground velocity structures. Domestic patents, such as CN108646302A (A Method for Detecting Underground Targets Based on Micromotion Signals), describe a method for velocity tomography using micromotion signals. However, existing micromotion exploration methods generally suffer from the following drawbacks:

[0006] (1) Existing micro-motion exploration methods are mainly aimed at one-dimensional velocity structure inversion, which can only obtain the vertical velocity profile directly below the measuring point and lack the ability to locate and characterize the three-dimensional spatial anomalies in the underground.

[0007] (2) Existing methods are not good at identifying linear targets such as karst conduits and cannot distinguish the difference in the dispersion curve response between conduits and layered low-velocity strata;

[0008] (3) Existing methods lack a specific extraction mechanism for the characteristic micro-motion response excited by fluid pulsation in karst pipes, and fail to make full use of the active "sound source" characteristics of the pipe fluid;

[0009] (4) Existing methods only provide two-dimensional profiles or simple planar distribution maps, and cannot output continuous three-dimensional trajectory models with practical engineering value.

[0010] In summary, existing technologies for the detailed three-dimensional detection of concealed karst pipelines suffer from major drawbacks, including insufficient depth, weak three-dimensional positioning capabilities, poor adaptability to urban environments, and high costs. Therefore, a method is needed that can achieve high-precision quantitative tracking of the three-dimensional trajectory of karst pipelines under passive source conditions. Summary of the Invention

[0011] To address the shortcomings of existing technologies, such as insufficient depth, weak 3D positioning capability, poor adaptability to urban environments, and high cost, the present invention aims to provide a 3D trajectory tracking method for karst pipelines based on array-type micro-motion signal interferometry imaging, comprising the following steps:

[0012] S1. Deploy multiple three-component micro-motion detector sites on the surface of the target exploration area to synchronously and continuously collect environmental micro-motion signals.

[0013] S2. Preprocess the environmental micro-motion signal to obtain continuous micro-motion waveform data for each station;

[0014] S3. Based on the continuous micro-motion waveform data, construct the cross-power spectral density matrix between the station pairs, and extract the characteristic micro-motion dispersion information that characterizes the fluid pulsation in the underground karst pipeline through spatial coherence function analysis;

[0015] S4. Using the aforementioned characteristic micro-motion dispersion information as input, a three-dimensional interferometric focusing imaging algorithm is employed to calculate the superposition intensity of the seismic source signal for a preset underground three-dimensional spatial grid unit, thereby constructing a three-dimensional energy distribution model of the target detection area;

[0016] S5. Identify the banded high-energy anomaly zone characterizing the activity of the karst conduit in the three-dimensional energy distribution model, and reconstruct the continuous three-dimensional spatial trajectory of the karst conduit by combining the spatial distribution information and micro-motion waveform characteristics of each detector station, and calculate the corresponding spatial positioning confidence distribution.

[0017] Furthermore, the step of S5 reconstructing the continuous three-dimensional spatial trajectory of the karst conduit specifically includes: S2.1. Performing a depth-by-depth horizontal slice analysis on the three-dimensional energy distribution model and extracting the energy maxima points on each depth slice;

[0018] S2.2. Perform a spatial continuity test on the energy maximum point, fit the principal axis direction of the banded high-energy anomaly area on each depth slice, and generate the planar projection trajectory of the karst conduit;

[0019] S2.3. Divide the multiple three-component micro-motion geophone stations into multiple local geophone groups, and calculate the spatial tilt angle of the karst pipe based on the imaging travel time difference of the target position of the planar projection trajectory in different local geophone groups.

[0020] S2.4. Polarization analysis is performed using the three-component particle motion vector of the environmental micro-motion signal to extract the incoming wave direction angle, and then combined with the spatial tilt angle to determine the three-dimensional spatial orientation of the karst conduit.

[0021] Furthermore, in the step of deploying multiple three-component micro-motion detector sites on the surface of the target exploration area, the array's deployment geometry includes a multi-ring concentric nested array or a mutually orthogonal L-shaped nested array.

[0022] Furthermore, the preprocessing step includes: performing mean removal and instrument response correction processing on the original environmental micro-motion signal;

[0023] A time-frequency domain joint normalization method is used to suppress transient strong noise events. In this method, bit normalization is performed on the waveform in the time domain, and spectral whitening is performed in the frequency domain.

[0024] The processed data is divided into windows according to the preset time window length, and abnormal time windows that are strongly polluted are removed based on the root mean square threshold of the signal energy in each window.

[0025] Furthermore, the three-dimensional interferometric focusing imaging algorithm includes: establishing the underground three-dimensional spatial grid unit in the target detection area; obtaining the shear wave velocity structure model of the target area; for each three-dimensional spatial grid unit, calculating the theoretical travel time from each station to the unit based on the shear wave velocity structure model; aligning the cross-correlation functions of all station pairs according to the time shift of the theoretical travel time difference and then superimposing them; and using the superimposed energy amplitude as the initial source signal superposition intensity of the three-dimensional spatial grid unit.

[0026] Furthermore, the step of calculating the superposition intensity of the source signals also includes a phase-weighted superposition step, specifically including: calculating the phase consistency index of the cross-correlation function of all station pairs at the time-shift alignment time, wherein the phase consistency index is defined as: the ratio of the superposition amplitude after normalization of the cross-correlation function of all station pairs to the total number of station pairs participating in the superposition.

[0027] The final superposition intensity of the source signals of the three-dimensional spatial grid cell is obtained by multiplying the initial source signal superposition intensity by a preset power of the phase consistency index.

[0028] Furthermore, the steps of reconstructing the continuous three-dimensional spatial trajectory and spatial positioning confidence distribution of the karst pipeline specifically include: selecting multiple different characteristic frequency bands to independently execute the three-dimensional interferometric focusing imaging algorithm to obtain the three-dimensional energy distribution model of multiple independent frequency bands; calculating adaptive weighting coefficients based on the imaging resolution index and signal-to-noise ratio level within the detection depth range corresponding to each frequency band; using the adaptive weighting coefficients to perform weighted summation and fusion of the three-dimensional energy distribution models of the multiple independent frequency bands to output the final continuous three-dimensional spatial trajectory model of the karst pipeline; calculating the spatial positioning confidence distribution based on the beam directivity index of energy focusing within the three-dimensional spatial grid cell and the consistency deviation between the imaging results of the multiple independent frequency bands; wherein, a confidence ellipsoid is used to characterize the spatial positioning uncertainty of each trajectory point in the continuous three-dimensional spatial trajectory model.

[0029] It also includes a karst conduit three-dimensional trajectory tracking system based on array-type micro-motion signal interferometry imaging, comprising:

[0030] The signal acquisition module is used to acquire environmental micro-motion signals synchronously and continuously collected by multiple three-component micro-motion detector sites deployed on the surface of the target exploration area.

[0031] The data preprocessing module is used to preprocess the environmental micro-motion signals to obtain continuous micro-motion waveform data of each station;

[0032] The coherence analysis module is used to construct the cross power spectral density matrix between station pairs based on the continuous micro-motion waveform data, and to extract characteristic micro-motion dispersion information that characterizes fluid pulsation in underground karst pipes through spatial coherence function analysis.

[0033] The three-dimensional interferometric imaging module is used to calculate the superposition intensity of the source signal for a preset underground three-dimensional spatial grid unit by taking the characteristic micro-motion dispersion information as input and using the three-dimensional interferometric focusing imaging algorithm, and constructing a three-dimensional energy distribution model of the target detection area.

[0034] The trajectory reconstruction module is used to identify the banded high-energy anomaly zone characterizing the activity of karst conduits in the three-dimensional energy distribution model, and to reconstruct the continuous three-dimensional spatial trajectory of the karst conduits by combining the spatial distribution information and micro-waveform characteristics of each detector station.

[0035] Furthermore, it includes a processor and a memory, the memory storing a computer program, the processor executing the computer program to implement the steps of the method.

[0036] The beneficial effects of this invention are as follows: Compared with the prior art, this invention adopts an anti-interference non-destructive passive source detection method, and significantly improves the target recognition signal-to-noise ratio through unique fluid pulsation feature extraction and multi-band joint imaging technology, breaking through the limitations of conventional one-dimensional exploration, and realizing high-precision true three-dimensional trajectory reconstruction of the direction, dip angle and burial depth of medium and deep karst pipelines. Attached Figure Description

[0037] Figure 1 This is the overall flowchart of the present invention;

[0038] Figure 2 This is a schematic diagram of the geometric layout of the multi-ring nested three-component detector array of the present invention;

[0039] Figure 3 This is a schematic diagram of the array micro-motion signal interference focusing imaging principle of the present invention;

[0040] Figure 4 This is a schematic diagram of the three-dimensional energy distribution model and pipeline trajectory extraction of the present invention;

[0041] Figure 5 This is a structural block diagram of the karst pipeline three-dimensional trajectory tracking device based on array-type micro-motion signal interferometry imaging according to the present invention. Detailed Implementation

[0042] To make the implementation objectives, technical solutions, and features of this invention application clearer, the technical solutions of this invention application will be clearly and completely described below with reference to the accompanying drawings of the embodiments. Obviously, the described embodiments are only a part of the examples of this invention application, and not all of the implementation cases. The embodiments of this invention application described and shown in the accompanying drawings can generally be arranged and designed in different configurations. Example

[0043] This embodiment provides three-dimensional trajectory tracking for karst conduits in standard plain areas, such as... Figure 1 As shown, the method in this embodiment mainly includes the following steps:

[0044] Step S1:

[0045] S1. Three-component micro-motion geophone stations are geometrically deployed in a multi-ring nested circular array on the surface of the target exploration area to synchronously and continuously acquire environmental micro-motion signals. Specifically, two concentric circular arrays are deployed: the inner circle (radius R1=20m) has 6 three-component geophone stations evenly deployed, with a spacing of approximately 20.9m between each station; the outer circle (radius R2=60m) has 12 three-component geophone stations evenly deployed, with a spacing of approximately 30.9m between each station; a central reference station is deployed at the center of the inner circle. The overall array aperture D=120m, and the target detection depth H=80m, satisfying the design constraint 0.3H(24m)≤D(120m)≤1.5H(120m); each three-component geophone station is equipped with one broadband three-component velocity geophone. The vertical component Z of the geophone is orthogonally installed to the north and east of the horizontal component using a leveling bubble, and is firmly coupled to the bedrock or compacted soil layer with a special steel cone.

[0046] All data acquisition systems at all sites achieve strict clock synchronization via a GPS timing module, with a synchronization accuracy better than 0.1ms. Acquisition parameters are set as follows: sampling rate 1000Hz, continuous acquisition duration 12 hours, recording on-site environmental conditions such as temperature and wind speed during the acquisition period.

[0047] S2. Micro-motion signal preprocessing to obtain continuous micro-motion waveform data from each station; the preprocessing steps include: performing mean removal and instrument response correction processing on the original environmental micro-motion signals;

[0048] A time-frequency domain joint normalization method is used to suppress transient strong noise events. In this method, bit normalization is performed on the waveform in the time domain, and spectral whitening is performed in the frequency domain.

[0049] The processed data was divided into windows according to a preset time window length, with a time window length of 50 seconds (covering 10 to 250 cycles in the target frequency range of 0.5Hz to 5Hz). The normalized continuous micro-motion waveform data was windowed with a 75% overlap rate, resulting in approximately 1680 short-time records. Abnormal time windows contaminated by strong interference were removed based on the root mean square threshold of the signal energy within each window. The root mean square (RMS) value of the energy of each of the three components was calculated for each time window, and abnormal time windows where the RMS of any component exceeded three times the global median were removed. Finally, approximately 1400 valid time windows were retained for subsequent analysis.

[0050] S3. Spatial coherence function analysis and characteristic dispersion information extraction: The following analysis is performed on the retained effective time window:

[0051] 1. Construction of cross-power spectral density matrix: For the vertical component Z continuous micro-motion waveform data of all N=19 stations, calculate the N×N cross-power spectral density matrix C(f) at each frequency point f. The matrix element C_{ij}(f) represents the cross-power spectral density between the i-th station and the j-th station at frequency f. The cross-power spectral density matrices of all effective time windows are superimposed and averaged.

[0052] 2. Extended Spatial Autocorrelation Analysis: All stations are grouped according to the station spacing r. In this embodiment, they are divided into four spacing groups: r≈10m, r≈20m, r≈30m, and r≈60m. The average spatial coherence coefficient as a function of frequency ρ(f,r) is calculated for each group. Based on the Bessel function theory, the phase velocity c(f) in ρ(f,r)=J0[2πf·r / c(f)] is fitted to obtain the dispersion curve.

[0053] 3. Identification of fluid pulsation characteristics in karst pipes: The pipe humming effect excited by the fluid in the pipe wall in karst pipes manifests as an abnormal decrease in phase velocity in a specific frequency range (usually 1Hz to 3Hz) on the dispersion curve. Combined with the azimuth clustering characteristics of the dominant propagation direction in the corresponding frequency band in the FK spectrum analysis (frequency-wavenumber spectrum, with the strongest coherence in the direction perpendicular to the pipe direction), the characteristic frequency range f_c1 to f_c2 of the fluid pulsation excitation component in the pipe is located. In this embodiment, the characteristic frequency range is identified as 1.2Hz to 2.8Hz, corresponding to the response of karst pipes at depths of 40m to 90m underground.

[0054] S4. Three-dimensional interferometric focusing imaging

[0055] like Figure 3 As shown, based on the feature dispersion information extracted in step S3, three-dimensional interferometric focusing imaging is performed:

[0056] (1) Velocity model establishment: Using the dispersion curve obtained in step S103, a genetic algorithm is used to invert the Rayleigh wave dispersion curve to obtain a layered model of shear wave velocity (Vs) within the depth range of 0 to 100m in the target area. In this embodiment, the velocity structure obtained by inversion is as follows: 0 to 5m, Vs=180m / s; 5 to 20m, Vs=280m / s; 20 to 50m, Vs=450m / s; 50 to 100m, Vs=720m / s (intact limestone); the stratum where the karst conduit is located shows a relatively low velocity anomaly in the velocity structure;

[0057] (2) Establishment of three-dimensional spatial grid unit: Establish regular three-dimensional spatial grid unit in underground three-dimensional space, with a horizontal range covering the area inside the array and the area within the outer circle (140m×140m), and a depth range of 5m to 150m. The size of the three-dimensional spatial grid unit is uniformly set to 2m×2m×2m (approximately 1 / 2 to 2 / 3 of the target pipe diameter of 3 to 5m).

[0058] (3) Phase-weighted interferometric focusing superposition: For each three-dimensional spatial grid cell V(x,y,z) in the three-dimensional mesh, perform the following calculations:

[0059] a. Based on the velocity model, calculate the S-wave travel time τᵢ from the i-th station to the three-dimensional spatial grid cell V;

[0060] b. For all site pairs (i,j), align the preprocessed cross-correlation function by a time shift of (τᵢ - τⱼ);

[0061] c. Superimpose the cross-correlation functions of all stations after time-shift alignment, and calculate the amplitude E(V) of the superimposed waveform at zero time delay, which is used as the focused imaging value of the three-dimensional spatial grid cell;

[0062] d. Simultaneously calculate the phase consistency index P(V) of the cross-correlation function of each station at the alignment time (i.e., the ratio of the normalized superimposed amplitude of the cross-correlation function of each station to the total number of samples), and use it as the phase weighted superposition (PWS) weight;

[0063] e. The final focused imaging value is E_PWS(V) = E(V) × P(V)^γ, where γ is the phase weighting exponent, and in this embodiment, γ=2;

[0064] (4) Local normalization: Divide the three-dimensional energy distribution model E_PWS(V) by the local background energy on each depth slice (calculate the local mean using 5×5 three-dimensional spatial grid cells as the window) to obtain the normalized three-dimensional energy distribution model, such as Figure 4 As shown.

[0065] S5 Pipeline Planar Trajectory Recognition

[0066] like Figure 4 As shown, a depth-by-depth horizontal slice analysis is performed on the three-dimensional energy distribution model:

[0067] (1) Extract local energy maxima for each horizontal slice at depth z, set the energy threshold to 60% of the slice peak value, and identify suspected pipeline anomalies;

[0068] (2) Determine the spatial continuity of anomalous bodies in each depth slice: It is required that there are related anomalous bodies in three vertically adjacent slices (within a total depth range of 6m) to eliminate isolated noise interference;

[0069] (3) For the anomalies that pass the continuity test, fit the main axis direction of the high-energy anomaly zone in the anomaly band on each depth slice to obtain the estimated value of the direction angle of the pipeline plane projection at each depth.

[0070] (4) Project the center points of the anomalies at different depths onto the ground surface to form the planar projection trajectory of the pipeline. In this embodiment, a continuous karst pipeline axial high-energy anomaly zone with a direction of NW45° and a planar projection length of about 85m was identified in the range of depth from 45m to 90m.

[0071] Step S6: Joint constraint inversion of the pipeline's 3D orientation and inclination. Based on the planar trajectory in Step S5, perform joint constraint inversion:

[0072] (1) Local detector grouping: All 19 sites are divided into 4 local detector groups according to the four quadrants NE, NW, SW and SE. Each local detector group contains at least 6 sites.

[0073] (2) Independent focusing imaging of local detector groups: Three-dimensional interferometric focusing imaging of each local detector group in each azimuth is performed independently to obtain three-dimensional energy distribution models in four azimuths;

[0074] (3) Calculation of travel time difference: At the pipeline axis position determined in step S5, compare the peak energy depth of the pipeline target in the focused imaging of local detector groups in different azimuths, calculate the depth difference between local detector groups in different azimuths, and use the geometric relationship (depth difference divided by the horizontal spacing of local detector groups) to calculate the longitudinal inclination angle δ of the pipeline; the calculation result in this embodiment is: the pipeline is inclined downward in the northwest at an inclination angle of about 12° along the NW direction;

[0075] (4) Particle motion polarization analysis: Extract the particle motion ellipse parameters (principal axis direction, ellipticity) of each three-component detector station in the characteristic frequency f_c range along the pipeline direction determined in step S105. Based on the principle that the horizontal direction of the principal axis of the Rayleigh wave particle motion ellipse corresponds to the main surface wave source direction, accurately determine the azimuth angle of the surface wave from the pipeline direction.

[0076] (5) Least square joint fitting: The tilt angle and azimuth angle calculated by the local detector group travel time difference decomposition are used as initial constraints and combined with the polarization analysis results for joint least square fitting to finally determine the three-dimensional direction angle (NW305°) and tilt angle (12°, tilting downward in the NW direction) of the pipeline.

[0077] Step S7: Multi-band weighted fusion and 3D model output

[0078] (1) Multi-band imaging: In step S104, three-dimensional interferometric focusing imaging is independently completed for the three frequency bands: low frequency band (0.5Hz to 1.5Hz, corresponding to a detection depth of 80m to 150m), mid frequency band (1.2Hz to 3Hz, corresponding to a detection depth of 30m to 90m) and high frequency band (2.5Hz to 5Hz, corresponding to a detection depth of 10m to 40m);

[0079] (2) Adaptive weight calculation: For each frequency band, calculate the imaging signal-to-noise ratio (the ratio of abnormal peak energy of the pipe to the background noise energy) in the target pipe area, and the resolution index (the lateral resolution is about λ / 2, where λ is the Rayleigh wave wavelength of the corresponding frequency) in the target pipe area; use the product of the two as the fusion weight of each frequency band.

[0080] (3) Multi-band fusion: The normalized three-dimensional energy distribution models of the three frequency bands are weighted and summed according to the above weights to obtain the final fused three-dimensional imaging result;

[0081] (4) Construction of three-dimensional trajectory model: such as Figure 4 As shown, a continuous set of three-dimensional trajectory points is extracted along the center of the anomaly maximum value of the pipeline in the fused three-dimensional energy distribution model. The distance between adjacent trajectory points is set to 5m, and they are connected to form a three-dimensional trajectory polyline. For each trajectory point, a confidence ellipsoid 305 is constructed with the weighted average of the lateral resolution of each participating frequency band as the semi-axis length.

[0082] (5) Output results: The final output is a three-dimensional trajectory model (including the three-dimensional coordinates, azimuth, tilt and confidence ellipsoid parameters of each trajectory point), as well as the energy density map of each depth slice and the data volume file of the fused three-dimensional energy distribution model.

[0083] Example 2: A 3D Trajectory Tracking System for Karst Pipes Based on Array-based Micro-motion Signal Interferometric Imaging. Based on the same inventive concept as Example 1 above, this example provides a tracking system, the structure of which is shown in the attached figure. Figure 5 As shown. The system includes a signal acquisition module, used to acquire environmental micro-motion signals synchronously and continuously collected by multiple three-component micro-motion detector stations deployed on the surface of the target exploration area;

[0084] The data preprocessing module is used to preprocess the environmental micro-motion signals to obtain continuous micro-motion waveform data of each station;

[0085] The coherence analysis module is used to construct the cross power spectral density matrix between station pairs based on the continuous micro-motion waveform data, and to extract characteristic micro-motion dispersion information that characterizes fluid pulsation in underground karst pipes through spatial coherence function analysis.

[0086] The three-dimensional interferometric imaging module is used to calculate the superposition intensity of the source signal for a preset underground three-dimensional spatial grid unit by taking the characteristic micro-motion dispersion information as input and using the three-dimensional interferometric focusing imaging algorithm, and constructing a three-dimensional energy distribution model of the target detection area.

[0087] The trajectory reconstruction module is used to identify the banded high-energy anomaly zone characterizing the activity of karst conduits in the three-dimensional energy distribution model, and to reconstruct the continuous three-dimensional spatial trajectory of the karst conduits by combining the spatial distribution information and micro-waveform characteristics of each detector station.

[0088] It also includes a processor and a memory, the memory storing a computer program, the processor executing the computer program to implement the steps of the method.

[0089] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for three-dimensional trajectory tracking of karst pipes based on array-type micro-motion signal interferometric imaging, characterized in that, Includes the following steps: S1 deploys multiple three-component micro-motion detector sites on the surface of the target exploration area to synchronously and continuously collect environmental micro-motion signals. S2 preprocesses the environmental micro-motion signal to obtain continuous micro-motion waveform data for each station; S3 constructs the cross-power spectral density matrix between station pairs based on the continuous micro-motion waveform data, and extracts characteristic micro-motion dispersion information that characterizes fluid pulsation in underground karst pipes through spatial coherence function analysis; S4 takes the aforementioned characteristic micro-motion dispersion information as input, uses a three-dimensional interferometric focusing imaging algorithm, calculates the superposition intensity of the source signal for the preset underground three-dimensional spatial grid unit, and constructs a three-dimensional energy distribution model of the target detection area; S5 identifies banded high-energy anomaly zones characterizing the activity of karst conduits in the three-dimensional energy distribution model, and reconstructs the continuous three-dimensional spatial trajectory of the karst conduit by combining the spatial distribution information and micro-motion waveform characteristics of each detector station, and calculates the corresponding spatial positioning confidence distribution; the step of S5 reconstructing the continuous three-dimensional spatial trajectory of the karst conduit specifically includes: S2.1 Perform depth-by-depth horizontal slicing analysis on the three-dimensional energy distribution model to extract the energy maxima points on each depth slice; S2.2 Perform spatial continuity test on the energy maximum point, fit the principal axis direction of the banded high-energy anomaly area on each depth slice, and generate the planar projection trajectory of the karst conduit; S2.3 Divide the multiple three-component micro-motion geophone stations into multiple local geophone groups, and calculate the spatial tilt angle of the karst pipe based on the imaging travel time difference of the target position of the planar projection trajectory in different local geophone groups. S2.4 Polarization analysis is performed using the three-component particle motion vector of the environmental micro-motion signal to extract the incoming wave direction angle, and then combined with the spatial tilt angle to determine the three-dimensional spatial orientation of the karst conduit. The three-dimensional interferometric focusing imaging algorithm includes: establishing the underground three-dimensional spatial grid unit in the detection area; obtaining the shear wave velocity structure model of the target area; for each three-dimensional spatial grid unit, calculating the theoretical travel time from each station to the unit based on the shear wave velocity structure model; aligning the cross-correlation functions of all station pairs according to the time shift of the theoretical travel time difference and then superimposing them; and using the superimposed energy amplitude as the initial source signal superposition intensity of the three-dimensional spatial grid unit. The step of calculating the superposition intensity of the source signals also includes a phase-weighted superposition step, which specifically includes: calculating the phase consistency index of the cross-correlation function of all station pairs at the time-shift alignment time, wherein the phase consistency index is defined as the ratio of the superposition amplitude after normalization of the cross-correlation function of all station pairs to the total number of station pairs participating in the superposition; multiplying the initial superposition intensity of the source signals by a preset power of the phase consistency index to obtain the final superposition intensity of the source signals of the three-dimensional spatial grid cell; The steps for reconstructing the continuous three-dimensional spatial trajectory and spatial positioning confidence distribution of the karst pipeline specifically include: selecting multiple different characteristic frequency bands and independently executing the three-dimensional interferometric focusing imaging algorithm to obtain the three-dimensional energy distribution model of multiple independent frequency bands; calculating adaptive weighting coefficients based on the imaging resolution index and signal-to-noise ratio level within the detection depth range corresponding to each frequency band; using the adaptive weighting coefficients to perform weighted summation and fusion of the three-dimensional energy distribution models of the multiple independent frequency bands to output the final continuous three-dimensional spatial trajectory model of the karst pipeline; calculating the spatial positioning confidence distribution based on the beam directivity index of energy focusing within the three-dimensional spatial grid cell and the consistency deviation between the imaging results of the multiple independent frequency bands; wherein, a confidence ellipsoid is used to characterize the spatial positioning uncertainty of each trajectory point in the continuous three-dimensional spatial trajectory model.

2. The method for three-dimensional trajectory tracking of karst pipes based on array-type micro-motion signal interferometry imaging according to claim 1, characterized in that, In the step of deploying multiple three-component micro-motion detector sites on the surface of the target exploration area, the array layout geometry includes a multi-ring concentric nested array or a mutually orthogonal L-shaped nested array.

3. The method for three-dimensional trajectory tracking of karst pipes based on array-type micro-motion signal interferometric imaging according to claim 1, characterized in that, The preprocessing steps include: performing mean removal and instrument response correction processing on the original environmental micro-motion signal; A time-frequency domain joint normalization method is used to suppress transient strong noise events. In this method, bit normalization is performed on the waveform in the time domain, and spectral whitening is performed in the frequency domain. The processed data is divided into windows according to the preset time window length, and abnormal time windows that are strongly polluted are removed based on the root mean square threshold of the signal energy in each window.

4. A three-dimensional trajectory tracking system for karst pipes based on array-type micro-motion signal interferometric imaging, characterized in that, The system is used to perform the method according to any one of claims 1 to 3, the system comprising: The signal acquisition module is used to acquire environmental micro-motion signals synchronously and continuously collected by multiple three-component micro-motion detector sites deployed on the surface of the target exploration area. The data preprocessing module is used to preprocess the environmental micro-motion signals to obtain continuous micro-motion waveform data of each station; The coherence analysis module is used to construct the cross power spectral density matrix between station pairs based on the continuous micro-motion waveform data, and to extract characteristic micro-motion dispersion information that characterizes fluid pulsation in underground karst pipes through spatial coherence function analysis. The three-dimensional interferometric imaging module is used to calculate the superposition intensity of the source signal for a preset underground three-dimensional spatial grid unit by taking the characteristic micro-motion dispersion information as input and using the three-dimensional interferometric focusing imaging algorithm, and constructing a three-dimensional energy distribution model of the target detection area. The trajectory reconstruction module is used to identify the banded high-energy anomaly zone characterizing the activity of karst conduits in the three-dimensional energy distribution model, and to reconstruct the continuous three-dimensional spatial trajectory of the karst conduits by combining the spatial distribution information and micro-waveform characteristics of each detector station.

5. The karst pipe three-dimensional trajectory tracking system based on array-type micro-motion signal interferometric imaging according to claim 4, characterized in that, The system further includes a processor and a memory, the memory storing a computer program, the processor executing the computer program to implement the steps of the method according to any one of claims 1 to 3.

Images