Methods, devices, equipment, media, and products for ground subsidence monitoring based on passive source seismic data.

Abstract

This invention discloses a ground subsidence monitoring method, device, equipment, medium, and product based on passive source seismic data, relating to the field of geological disaster monitoring technology. The method first performs cross-correlation analysis on corresponding passive source seismic data collected in the most recent period and multiple consecutive periods for each pair of nodes. Based on the analysis results, multiple delay time estimates and correlation coefficients are obtained through time-domain segmentation and correlation processing, and frequency-domain processing and correlation calculations. Then, based on the calculation results, the corresponding valid data are selected. Finally, based on the valid data and geographical location of each pair of nodes, the distribution of seismic wave velocity changes in the subsurface medium of the target area in the most recent period is analyzed, and ground subsidence early warning actions are triggered accordingly. This effectively improves the accuracy and timeliness of ground subsidence early warning and has advantages such as continuous monitoring, simplicity, efficiency, and controllable cost, making it easy to apply and promote in practice.

Description

Technical Field

[0001] This invention belongs to the field of geological disaster monitoring technology, specifically relating to a ground subsidence monitoring method, device, equipment, medium, and product based on passive source seismic data. Background Technology

[0002] Against the backdrop of rapid urban expansion, ground collapses, such as those on urban roads, are becoming increasingly frequent. Road collapses severely disrupt transportation on major thoroughfares, leading to traffic congestion and accidents. This not only affects the smooth flow of traffic but also causes vehicle damage and pedestrian injuries or even deaths, posing a serious threat to urban development and the safety of people's lives and property, and impacting the sustainable development of the social economy.

[0003] To address the hazards caused by ground / road collapses, various research and engineering institutions have tried numerous methods and achieved certain results. Among them, commonly used monitoring methods for road collapse problems include total station measurement, GPS (Global Positioning System) monitoring, satellite remote sensing, laser scanning, and video monitoring.

[0004] Total station surveying is a traditional and classic method for monitoring the earth's surface. Utilizing its high-precision measurement capabilities, a total station can measure the horizontal and vertical angles and distances at different points on the earth's surface, thereby determining the shape and elevation of the surface. This technique is commonly used in engineering surveying, but it can also be applied to monitoring road subsidence. By regularly conducting total station measurements on the surface around roads, minute changes in the surface, including subsidence and uplift, can be monitored, providing reliable data support for timely intervention. However, it is labor-intensive and time-consuming.

[0005] GPS monitoring technology can observe the movement and deformation of the earth's surface in real time by measuring its position and elevation through GPS receivers installed on the surface. By continuously monitoring changes in the position of the GPS receivers, surface deformation, including road subsidence or movement, can be detected. This technology has advantages such as high real-time performance and wide coverage, enabling effective monitoring of surface changes over a large area. However, its data accuracy may be limited by the influence of satellite signals, especially in densely built-up or high-rise urban areas. Furthermore, weather conditions may affect GPS signal reception, impacting monitoring accuracy.

[0006] Satellite remote sensing technology acquires images and data of the Earth's surface using sensors on satellites or aerial platforms to monitor changes in the land surface. This method can be used to detect signs of surface subsidence, cracks, and landslides around roads, enabling timely analysis and early warning. Remote sensing technology offers advantages such as rapid data acquisition and wide coverage, providing crucial support for urban infrastructure management. However, due to the periodic nature of remote sensing data acquisition, real-time monitoring may not be possible. Secondly, it requires specialized image interpretation techniques and geoscientific knowledge, and the interpretation process can be complex. Furthermore, cloud cover can affect the quality of remote sensing images and the accuracy of interpretation.

[0007] Laser scanning technology is also widely used in the field of land surface monitoring. It uses laser beams to measure the shape and elevation of the land surface, enabling the acquisition of a high-precision three-dimensional model. Regular laser scanning of the surface around roads allows for the detection of minute surface changes, enabling timely repair and reinforcement. Laser scanning technology offers advantages such as high precision and efficiency, making it a crucial technique in land surface monitoring. However, the purchase and maintenance costs of laser scanning equipment are relatively high, the requirements for identifying and interpreting ground features are demanding, it may be affected by surface cover, and the processing and analysis of laser scanning data requires specialized software and technical support, making the operation relatively complex.

[0008] Video surveillance is a simple yet effective method for monitoring land surfaces. Installing cameras around roads and periodically taking photos or videos of the surface allows for real-time monitoring of changes. By comparing photos or videos from different points in time, issues such as subsidence and cracks can be detected and addressed promptly. Video surveillance offers advantages such as low cost and ease of installation, making it widely applicable in urban infrastructure management. However, it also has drawbacks, including susceptibility to weather conditions, the need for a continuous energy supply, and complex data processing.

[0009] The aforementioned monitoring methods all fall under the category of surface monitoring, primarily targeting surface morphology and displacement. However, since ground / road subsidence often begins with changes in the underground structure and develops gradually, surface monitoring methods suffer from a "post-event" drawback—they only show results when the subsidence is imminent or has already occurred. Therefore, exploring a new ground / road subsidence monitoring method with features such as early diagnosis, continuous observation, simplicity, efficiency, and cost-effectiveness is crucial for overcoming the shortcomings of existing methods and enriching ground / road subsidence monitoring technology solutions. This is of great significance for improving the forecasting and early warning capabilities for urban ground / road subsidence disasters.

[0010] Ground / road subsidence is essentially a chain reaction resulting from the gradual accumulation of deterioration effects on the ground / road foundation. The most common factor is the action of water (e.g., leaking pipes or groundwater erosion). Under the influence of water, soil and rock undergo a series of significant physical property changes. For example, fluid diffusion in soil pores significantly increases local pore water pressure. Water flow erodes the cementing substances between soil particles (such as clay minerals), weakening the friction and interlocking between particles. The direct effect is a sharp decrease in the effective stress borne by the soil skeleton due to the increased pore water pressure, directly affecting the soil's shear modulus. Simultaneously, water flow dissolves carbonate / ferrous cementing materials, causing the gravel skeleton to degenerate from an "interlocked state" to a "suspended state," resulting in a simultaneous decrease in both bulk modulus and shear modulus. Furthermore, some fine particles are carried away by the water flow (i.e., undercutting), leading to a decrease in the average density of the medium.

[0011] The propagation speed of seismic waves in underground media is directly physically related to the elastic modulus and density of the media. This relationship can be described as follows: changes in wave velocity are essentially an external manifestation of the dynamic equilibrium between soil stiffness and mass. Taking shear waves (S-waves) as an example, their propagation speed v... s It is determined by both the shear modulus μ (which reflects the soil's ability to resist shear deformation) and the density ρ, i.e. The increase or decrease in shear modulus plays a dominant role in wave velocity. When the shear modulus μ of the underground medium increases significantly due to compaction, increased stress, or cementation (e.g., more than 30% after roadbed compaction), even if the density ρ increases simultaneously (usually by less than 5%), the wave velocity will still rise significantly due to the dominant contribution of the modulus. Conversely, in catastrophic scenarios (such as water-soil interaction caused by underground pipeline leaks), the infiltration of high-pressure fluids into the soil will increase pore water pressure, weakening the effective stress of the soil-rock skeleton, causing the shear modulus μ to decrease exponentially; at the same time, the scouring effect of water flow carrying away fine particles will lead to a slight decrease in density ρ. Under the combined effect of the above two factors, the wave velocity exhibits a systematic decrease, which can reach 1% to 2%.

[0012] Therefore, the sensitive attenuation of wave velocity in the ground / road subgrade is a physical precursor to the deterioration of soil mechanical properties, providing a crucial basis for early warning of ground / road subsidence risks. However, obtaining information about these velocity changes is a challenging task.

[0013] Currently, traditional methods for obtaining ground seismic wave velocities from time-series data mainly include the refraction wave method, reflection wave method, active source surface wave method, borehole wave velocity testing, cross-bore wave velocity testing, and cross-bore tomography. Specifically, the refraction wave method utilizes the refraction phenomenon of seismic waves at the interfaces of strata with different velocities to invert stratum velocities through the first arrival travel time. However, this method cannot detect low-velocity layers (such as loose zones) because low-velocity layers create a "blind zone," resulting in low resolution (typically >10m). The reflection wave method involves exciting seismic waves and receiving reflected signals from strata interfaces (interfaces with different wave impedances), constructing a subsurface profile using travel time and amplitude. However, this method has a significant shallow blind zone (lack of near-surface reflected waves), and the absorption by loose surface layers leads to high-frequency signal attenuation and a decrease in the dominant frequency. Furthermore, data processing requires complex static correction and noise suppression. The active source surface wave method excites Rayleigh waves by hammering or vibrating, extracts dispersion curves to invert shear wave velocity profiles, and directly reflects soil shear stiffness. However, the lateral resolution of this method is limited by wavelength (i.e., better at shallow depths and worse at deep depths), and it is easily affected by traffic vibrations in urban environments. Furthermore, the inversion results depend on the selection of the initial model. Borehole wave velocity testing involves placing detectors at different depths within the borehole and directly measuring the travel time of P-waves / S-waves by excitation at the ground surface to calculate in-situ soil dynamic parameters (such as shear modulus). However, this method has limited vertical resolution (usually >2m) and only reflects local information near the borehole, making it difficult to represent the condition of a large-scale subgrade. Cross-hole wave velocity testing involves transmitting and receiving seismic waves between adjacent boreholes to accurately calculate inter-layer wave velocities and achieve high-resolution (0.5–1 m) layering. However, this method requires the cooperation of multiple boreholes, is costly, and has low operational efficiency (advancing only tens of meters per hour), making it only suitable for detailed investigation of small, key areas. Cross-hole tomography is based on multiple sets of cross-hole travel time data, generating two-dimensional / three-dimensional velocity structures through ray tracing or waveform inversion. It can finely identify cavities and fractures (resolution up to 0.3 m). However, this method requires complex equipment and a huge amount of computation, data processing requires specialized software support, and the detection range is strictly limited by the borehole layout.

[0014] Passive-source seismic data has gradually become a promising research area due to its ability to overcome the limitations of dependence on seismic events. This is thanks to the possibility of retrieving seismic Green's functions from the cross-correlation functions of passive-source seismic wavefield records acquired from different locations within the study area. In fact, even when conditions do not allow for complete reconstruction, constructing seismic Green's functions using cross-correlation of passive-source seismic data has proven to be robust. Furthermore, since passive-source seismic data are continuously recorded and independent of seismic sources, these cross-correlation functions can be analogous to records from continuously repeating dual-source seismic systems placed at each station, and can be used to extract observations of seismic velocity variations.

[0015] Therefore, how to use continuous passive source seismic data to monitor the changes in the seismic wave velocity field of the underground medium in the target area, and to carry out ground subsidence monitoring and early warning based on these changes, so as to make up for the deficiencies of existing ground subsidence monitoring schemes, enrich ground subsidence monitoring technology schemes, and improve the ability to predict and warn of urban ground subsidence disasters, is a topic that urgently needs to be studied by those skilled in the art. Summary of the Invention

[0016] The purpose of this invention is to provide a ground subsidence monitoring method, device, computer equipment, computer-readable storage medium, and computer program product based on passive source seismic data, in order to solve the problem that existing ground subsidence monitoring schemes only issue warnings after ground subsidence has occurred due to the lack of "early diagnosis" features.

[0017] To achieve the above objectives, the present invention adopts the following technical solution:

[0018] Firstly, a ground subsidence monitoring method based on passive source seismic data is provided, including:

[0019] The system receives passive source seismic data collected in real time by multiple seismic monitoring nodes, wherein the multiple seismic monitoring nodes are discretely arranged in the target area.

[0020] For each pair of nodes among the multiple seismic monitoring nodes, based on a first sliding window with an overlap rate of a first percentage, cross-correlation analysis is periodically performed on the corresponding two passive source seismic data acquired in the most recent period to obtain N in the most recent period. cur A cross-correlation function, and then the N cur The cross-correlation functions are superimposed and normalized to obtain the corresponding current cross-correlation function d1(t), and the cross-correlation function d1(t) is also obtained by superimposing and normalizing the N most recent consecutive cross-correlation functions. ref The corresponding cross-correlation functions for each period are superimposed and normalized to obtain the corresponding reference cross-correlation function d0(t), where N cur Represents a positive integer greater than or equal to 10, and the duration of the current most recent period is N times the duration of the first sliding window. cur Times, N ref t represents a positive integer greater than or equal to 10, and t represents time.

[0021] For each pair of nodes, based on a second sliding window with an overlap rate of a second percentage, the corresponding reference cross-correlation function d0(t) is divided into N... w The first piece of data, and the corresponding current cross-correlation function d1(t) divided into N w The second piece of data, and on the N w First data and the N wEach data point in the second data set undergoes root mean square energy normalization, windowing, and cosine shrinking to obtain the corresponding N. w The first processed data and N w The second processed data, wherein the duration of the second sliding window is shorter than the duration of the first sliding window, and the second percentage is higher than the first percentage, N w Represents positive integers;

[0022] For each pair of nodes, apply the Fourier transform to obtain the corresponding N... w The first processed data and N w The frequency domain function corresponding to each piece of processed data in the second processed data is obtained, and the corresponding N is obtained based on the transformation result. w Mutual spectrum;

[0023] For each pair of nodes, based on all corresponding frequency domain functions and N w Based on the cross spectrum and a preset effective frequency range, the energy density is calculated and the corresponding N is obtained. w ×J mutually coherent function values ​​and a weighting coefficient corresponding to each mutually coherent function value, where J represents the total number of sampling frequencies in the effective frequency range;

[0024] For each pair of nodes, based on the corresponding N w Given J weighted coefficients, the corresponding N is obtained using the weighted least squares inversion method. w A slope, and based on that N w The slope is calculated to obtain the corresponding N. w Each delay time estimate, and based on the corresponding N w The corresponding N is obtained by calculating ×J mutually coherent function values. w One correlation coefficient;

[0025] For each pair of nodes, based on the preset delay time estimation threshold and correlation coefficient threshold, the corresponding N will be... w All second data in the second data set that satisfy the corresponding delay time estimate value being less than or equal to the delay time estimate threshold and the corresponding correlation coefficient being greater than or equal to the correlation coefficient threshold are considered as the corresponding valid data.

[0026] Based on the valid data and geographical location of each pair of nodes, the distribution of changes in seismic wave velocity in the subsurface medium of the target area during the most recent period is analyzed.

[0027] The ground collapse early warning action for the target area is triggered based on the distribution of changes in seismic wave velocity.

[0028] Based on the above-mentioned invention, a novel scheme for monitoring ground subsidence in a target area using continuous passive source seismic data is provided. This scheme involves receiving passive source seismic data collected in real-time by multiple seismic monitoring nodes. First, for each pair of nodes, cross-correlation analysis is periodically performed on corresponding passive source seismic data sets collected in the most recent period and the most recent consecutive periods. Based on the analysis results, multiple delay time estimates and correlation coefficients are obtained through time-domain segmentation and correlation processing, as well as frequency-domain processing and correlation calculations. Then, the corresponding valid data are selected based on the calculation results. Finally, based on the valid data of each pair of nodes and their geographical location, the distribution of seismic wave velocity changes in the subsurface medium of the target area in the most recent period is analyzed, triggering a ground subsidence early warning action. This overcomes the shortcomings of existing ground subsidence monitoring schemes, which lack "early diagnosis" features and only issue warnings after ground subsidence has occurred. This effectively improves the accuracy and timeliness of ground subsidence early warnings and offers advantages such as continuous monitoring, simplicity, efficiency, and controllable cost, making it suitable for practical application and widespread adoption.

[0029] In one possible design, for each pair of nodes among the plurality of seismic monitoring nodes, based on a first sliding window with an overlap rate of a first percentage, cross-correlation analysis is periodically performed on the corresponding two passive source seismic data acquired in the current most recent period to obtain N in the current most recent period and corresponding cur The cross-correlation functions include:

[0030] For a pair of nodes among the multiple earthquake monitoring nodes, periodically acquire the corresponding two passive source earthquake data collected in the most recent period;

[0031] The two sets of passive source seismic data of a certain pair of nodes are preprocessed to obtain two sets of preprocessed data. The preprocessing includes filtering, anomaly suppression, energy normalization and / or spectral whitening.

[0032] Based on a first sliding window with an overlap rate of a first percentage, cross-correlation analysis is performed on the two preprocessed data sets to obtain N of the pair of nodes in the current most recent period. cur There are N cross-correlation functions, where N cur Represents a positive integer greater than or equal to 10, and the duration of the current most recent period is N times the duration of the first sliding window. cur times.

[0033] In one possible design, when the preprocessing includes filtering, the filtering is used to extract a preferred frequency band signal with a frequency band range of 2.0 to 40 Hz.

[0034] In one possible design, for each pair of nodes, the corresponding N is obtained based on the transformation result. wEach cross spectrum includes:

[0035] For a given pair of nodes among the multiple earthquake monitoring nodes, the corresponding nth cross spectrum is obtained based on the transformation result. Where n represents less than or equal to N w positive integers, F 0,n (f) represents the frequency domain function corresponding to the nth slice of first processed data in time sequence. F represents the frequency domain function corresponding to the nth slice of second-processed data in time sequence. 1,n (f) is the conjugate form of f, where f represents the frequency.

[0036] In one possible design, for each pair of nodes, based on all corresponding frequency domain functions and N w Based on the cross spectrum and a preset effective frequency range, the energy density is calculated and the corresponding N is obtained. w ×J mutually coherent function values ​​and the corresponding weight coefficients for each mutually coherent function value, including:

[0037] For a given pair of nodes, based on all corresponding frequency domain functions and N w Based on the cross spectrum and the preset effective frequency range, the energy density is calculated according to the following formula, and the corresponding N is obtained. w ×J mutually coherent function values ​​and the corresponding weight coefficient for each mutually coherent function value:

[0038]

[0039] In the formula, j represents the index of the sampling frequency in the effective frequency range and is a positive integer less than or equal to J, J represents the total number of sampling frequencies in the effective frequency range, and C n,j Indicates the cross spectrum with the nth cross spectrum X n (f) and the cross-coherence function value corresponding to the j-th sampling frequency in the effective frequency range, X n,j This means substituting the j-th sampling frequency into the n-th cross spectrum X. n The function value F obtained after (f) 0,n,j This means substituting the j-th sampling frequency into the frequency domain function F. 0,n The function value F obtained after (f) 1,n,j This means substituting the j-th sampling frequency into the frequency domain function F. 1,n The function value obtained after (f), w n,j Represents the value of the mutually coherent function C. n,j The corresponding weighting coefficients, the effective frequency range is 5 to 20 Hz.

[0040] In one possible design, for each pair of nodes, based on the corresponding N wGiven J weighted coefficients, the corresponding N is obtained using the weighted least squares inversion method. w A slope, and based on that N w The slope is calculated to obtain the corresponding N. w Each delay time estimate, and based on the corresponding N w The corresponding N is obtained by calculating ×J mutually coherent function values. w The correlation coefficients include:

[0041] For a given pair of nodes, based on the corresponding N w Given ×J weighted coefficients, the corresponding N is obtained using the weighted least squares inversion method according to the following formula. w One slope:

[0042]

[0043] In the formula, k n Indicates the cross spectrum with the nth cross spectrum X n The slope corresponding to (f), f j φ represents the frequency value of the j-th sampling frequency. j This represents the phase information corresponding to the j-th sampling frequency;

[0044] Based on the N of the aforementioned pair of nodes w The slope is used to calculate N for a given pair of nodes according to the following formula. w One estimated delay time:

[0045] δt n =k n ÷(2×π)

[0046] In the formula, δt n Indicates the cross spectrum with the nth cross spectrum X n (f) The estimated delay time corresponding to;

[0047] Based on the N of the aforementioned pair of nodes w Given ×J mutually coherent function values, N for a given pair of nodes is calculated using the following formula. w One correlation coefficient:

[0048]

[0049] In the formula, C n Indicates the cross spectrum with the nth cross spectrum X n (f) Correlation coefficient.

[0050] Secondly, a ground collapse monitoring device based on passive source seismic data is provided, comprising a seismic data receiving unit, a cross-correlation analysis unit, a data fragmentation processing unit, a data frequency domain processing unit, a first calculation processing unit, a second calculation processing unit, an effective data filtering unit, a wave velocity change analysis unit, and a monitoring and early warning triggering unit that are connected in sequence.

[0051] The earthquake data receiving unit is used to receive passive source earthquake data collected in real time by multiple earthquake monitoring nodes, wherein the multiple earthquake monitoring nodes are discretely arranged in the target area.

[0052] The cross-correlation analysis unit is used to periodically perform cross-correlation analysis on corresponding passive source seismic data acquired in the most recent period for each pair of nodes among the plurality of seismic monitoring nodes, based on a first sliding window with an overlap rate of a first percentage, to obtain N corresponding to the current most recent period. cur A cross-correlation function, and then the N cur The cross-correlation functions are superimposed and normalized to obtain the corresponding current cross-correlation function d1(t), and the cross-correlation function d1(t) is also obtained by superimposing and normalizing the N most recent consecutive cross-correlation functions. ref The corresponding cross-correlation functions for each period are superimposed and normalized to obtain the corresponding reference cross-correlation function d0(t), where N cur Represents a positive integer greater than or equal to 10, and the duration of the current most recent period is N times the duration of the first sliding window. cur Times, N ref t represents a positive integer greater than or equal to 10, and t represents time.

[0053] The data sharding processing unit is used to divide the corresponding reference cross-correlation function d0(t) into N parts for each pair of nodes, based on a second sliding window with an overlap rate of a second percentage. w The first piece of data, and the corresponding current cross-correlation function d1(t) divided into N w The second piece of data, and on the N w First data and the N w Each data point in the second data set undergoes root mean square energy normalization, windowing, and cosine shrinking to obtain the corresponding N. w The first processed data and N w The second processed data, wherein the duration of the second sliding window is shorter than the duration of the first sliding window, and the second percentage is higher than the first percentage, N w Represents positive integers;

[0054] The data frequency domain processing unit is used to apply Fourier transform to each pair of nodes to obtain the corresponding N... w The first processed data and N wThe frequency domain function corresponding to each piece of processed data in the second processed data is obtained, and the corresponding N is obtained based on the transformation result. w Mutual spectrum;

[0055] The first computational processing unit is configured to, for each pair of nodes, perform computation based on all corresponding frequency domain functions and N. w Based on the cross spectrum and a preset effective frequency range, the energy density is calculated and the corresponding N is obtained. w ×J mutually coherent function values ​​and a weighting coefficient corresponding to each mutually coherent function value, where J represents the total number of sampling frequencies in the effective frequency range;

[0056] The second calculation processing unit is used to, for each pair of nodes, based on the corresponding N w Given J weighted coefficients, the corresponding N is obtained using the weighted least squares inversion method. w A slope, and based on that N w The slope is calculated to obtain the corresponding N. w Each delay time estimate, and based on the corresponding N w The corresponding N is obtained by calculating ×J mutually coherent function values. w One correlation coefficient;

[0057] The effective data filtering unit is used to, for each pair of nodes, filter the corresponding N nodes based on a preset delay time estimation threshold and a correlation coefficient threshold. w All second data in the second data set that satisfy the corresponding delay time estimate value being less than or equal to the delay time estimate threshold and the corresponding correlation coefficient being greater than or equal to the correlation coefficient threshold are considered as the corresponding valid data.

[0058] The wave velocity change analysis unit is used to analyze the distribution of seismic wave velocity changes in the subsurface medium of the target area in the current most recent period based on the effective data and geographical location of each pair of nodes.

[0059] The monitoring and early warning triggering unit is used to trigger a ground collapse early warning action for the target area based on the distribution of the seismic wave velocity change.

[0060] Thirdly, the present invention provides a computer device comprising a memory, a processor, and a transceiver connected in sequence for communication, wherein the memory is used to store a computer program, the transceiver is used to send and receive messages, and the processor is used to read the computer program and execute the ground subsidence monitoring method as described in the first aspect or any possible design in the first aspect.

[0061] Fourthly, the present invention provides a computer-readable storage medium storing instructions that, when executed on a computer, perform the ground subsidence monitoring method as described in the first aspect or any possible design within the first aspect.

[0062] Fifthly, the present invention provides a computer program product, including a computer program or instructions, which, when executed by a computer, implement the ground subsidence monitoring method as described in the first aspect or any possible design in the first aspect.

[0063] The beneficial effects of the above scheme are:

[0064] (1) This invention creatively provides a new scheme for monitoring ground subsidence in a target area based on continuous passive source seismic data. After receiving passive source seismic data collected in real time by multiple seismic monitoring nodes, cross-correlation analysis is performed periodically on the two corresponding passive source seismic data collected in the current most recent period and the current consecutive periods for each pair of nodes. Based on the analysis results, multiple delay time estimates and correlation coefficients are obtained through time domain segmentation and correlation processing and frequency domain processing and correlation calculation. Then, the corresponding effective data are obtained by filtering based on the calculation results. Finally, based on the effective data and geographical location of each pair of nodes, the distribution of seismic wave velocity changes in the underground medium of the target area in the current most recent period is analyzed, and ground subsidence early warning action is triggered accordingly. This can make up for the defects of existing ground subsidence monitoring schemes, which are not characterized by "early diagnosis" and thus issue early warnings only after ground subsidence. It effectively improves the accuracy and timeliness of ground subsidence early warning and has the advantages of continuous monitoring, simplicity, efficiency and controllable cost.

[0065] (2) It can also enrich the technical solutions for monitoring ground subsidence, improve the forecasting and early warning capabilities of urban ground subsidence disasters, and facilitate practical application and promotion. Attached Figure Description

[0066] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0067] Figure 1 This is a flowchart illustrating the ground subsidence monitoring method based on passive source seismic data provided in an embodiment of this application.

[0068] Figure 2 A schematic diagram of the structure of a ground subsidence monitoring device based on passive source seismic data provided in this application embodiment.

[0069] Figure 3 A schematic diagram of the structure of a computer device provided in an embodiment of this application. Detailed Implementation

[0070] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the present invention will be briefly introduced below in conjunction with the accompanying drawings and descriptions of the embodiments or the prior art. Obviously, the following description of the structure of the accompanying drawings is only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort. It should be noted that the description of these embodiments is for the purpose of helping to understand the present invention, but does not constitute a limitation of the present invention.

[0071] It should be understood that although the terms "first" and "second", etc., may be used herein to describe various objects, these objects should not be limited by these terms. These terms are only used to distinguish one object from another. For example, the first object may be referred to as the second object, and similarly, the second object may be referred to as the first object, without departing from the scope of the exemplary embodiments of the invention.

[0072] It should be understood that the term "and / or" that may appear in this document is merely a description of the relationship between related objects, indicating that three relationships can exist. For example, A and / or B can mean: A exists alone, B exists alone, or A and B exist simultaneously. Another example is A, B and / or C, which can mean that any one of A, B, and C or any combination thereof exists. The term " / and" that may appear in this document describes another relationship between related objects, indicating that two relationships can exist. For example, A / and B can mean: A exists alone or A and B exist simultaneously. In addition, the character " / " that may appear in this document generally indicates that the related objects before and after it are in an "or" relationship.

[0073] Example:

[0074] like Figure 1 As shown, the ground subsidence monitoring method based on passive source seismic data provided in the first aspect of this embodiment can be executed, but is not limited to, by a computer device with certain computing resources and a communication connection to the seismic monitoring node. For example, it can be executed by an electronic device such as a platform server, a personal computer (PC, referring to a multi-purpose computer of a size, price, and performance suitable for personal use; desktop computers, laptops, mini-laptops, tablets, and ultrabooks are all considered personal computers), a smartphone, a personal digital assistant (PDA), or a wearable device. Figure 1 As shown, the ground subsidence monitoring method may include, but is not limited to, the following steps S1 to S9.

[0075] S1. Receive passive source seismic data collected in real time by multiple seismic monitoring nodes, wherein the multiple seismic monitoring nodes are discretely arranged in the target area.

[0076] In step S1, the target area may be, but is not limited to, urban roads or public squares. The earthquake monitoring nodes are existing station equipment used to collect earthquake data, capable of routinely acquiring passive source earthquake data generated by the presence of passive seismic sources in real time. Specifically, the earthquake monitoring nodes can be connected to local equipment via existing wired or wireless IoT technologies. Multiple earthquake monitoring nodes may be, but are not limited to, discretely arranged in a matrix within the target area. The sampling rate of the earthquake monitoring nodes can be designed based on a comprehensive analysis of the monitoring target requirements and IoT data pressure, for example, designed to sample once every 4 milliseconds. The local equipment also needs to record the unique identifier of each earthquake monitoring node and its pre-measured geographical location. Furthermore, the local equipment also needs to obtain the reference seismic wave velocity (e.g., 200 m / s) and the approximate range of the target area (e.g., 100 m) in the target area.

[0077] S2. For each pair of nodes among the multiple seismic monitoring nodes, based on a first sliding window with an overlap rate of a first percentage, periodically perform cross-correlation analysis on the corresponding two passive source seismic data acquired in the current most recent period to obtain N in the current most recent period and corresponding... cur A cross-correlation function, and then the N cur The cross-correlation functions are superimposed and normalized to obtain the corresponding current cross-correlation function d1(t), and the cross-correlation function d1(t) is also obtained by superimposing and normalizing the N most recent consecutive cross-correlation functions. ref The corresponding cross-correlation functions for each period are superimposed and normalized to obtain the corresponding reference cross-correlation function d0(t), where N cur Represents a positive integer greater than or equal to 10, and the duration of the current most recent period is N times the duration of the first sliding window. cur Times, N ref t represents a positive integer greater than or equal to 10, and t represents time.

[0078] In step S2, the aforementioned period reflects the frequency of ground subsidence monitoring and early warning. It can be, but is not limited to, 1 hour. That is, for a pair of nodes among the multiple seismic monitoring nodes, the corresponding two passive source seismic data collected in the most recent hour can be acquired every hour, and cross-correlation analysis can be performed on them. To improve the quality of the analyzed data and results, preferably, for each pair of nodes among the multiple seismic monitoring nodes, based on a first sliding window with an overlap rate of a first percentage, cross-correlation analysis is periodically performed on the corresponding two passive source seismic data collected in the most recent period to obtain the N corresponding to the most recent period. cur Each cross-correlation function includes, but is not limited to, the following steps S21 to S23.

[0079] S21. For a pair of nodes among the plurality of earthquake monitoring nodes, periodically acquire the corresponding two passive source earthquake data collected in the most recent period.

[0080] In step S21, the pair of nodes corresponds one-to-one with the two sets of passive source seismic data.

[0081] S22. Preprocess the two sets of passive source seismic data of the pair of nodes respectively to obtain two sets of preprocessed data, wherein the preprocessing includes, but is not limited to, filtering, anomaly suppression, energy normalization and / or spectral whitening.

[0082] In step S22, the two sets of preprocessed data also correspond one-to-one with the pair of nodes. The specific processes of the filtering, abnormal interference suppression, energy normalization, and spectral whitening are existing conventional techniques and will not be described in detail here. Specifically, when the preprocessing includes filtering, the filtering is used to extract a preferred frequency band signal with a frequency range of 2.0–40 Hz; the aforementioned frequency band range can be predetermined by conducting multiple tests for urban ground / road subsidence monitoring.

[0083] S23. Based on a first sliding window with an overlap rate of a first percentage, perform cross-correlation analysis on the two preprocessed data sets to obtain the N of the pair of nodes in the current most recent period. cur There are N cross-correlation functions, where N cur Represents a positive integer greater than or equal to 10, and the duration of the current most recent period is N times the duration of the first sliding window. cur times.

[0084] In step S23, the first sliding window is used to segment the two preprocessed data respectively, so as to obtain N based on the segmentation. curPerform the cross-correlation analysis on the segmented data respectively to obtain N of the pair of nodes in the current most recent period. cur A cross-correlation function. The specific process of the cross-correlation analysis is a prior art and will not be described in detail here. For example, the duration of the first sliding window is 60 seconds, the first percentage is 0%, and N is... cur The value is 60, meaning that for a given pair of nodes, 60 cross-correlation functions corresponding to the most recent period can be obtained.

[0085] In step S2, the superposition and normalization process includes two actions: superposition and normalization. The specific processes of the aforementioned superposition and normalization are existing technologies and will not be described in detail here. The current cross-correlation function d1(t) reflects the current underground situation of the target area, and the reference cross-correlation function d0(t) reflects the underground historical background of the target area. Specifically, N... ref For example, if the value is 240 (i.e., the most recent ten days), then the number of all cross-correlation functions is 240 × 60 = 14400, which is far greater than N. cur .

[0086] S3. For each pair of nodes, based on a second sliding window with an overlap rate of the second percentage, the corresponding reference cross-correlation function d0(t) is divided into N... w The first piece of data, and the corresponding current cross-correlation function d1(t) divided into N w The second piece of data, and on the N w First data and the N w Each data point in the second data set undergoes root mean square energy normalization, windowing, and cosine shrinking to obtain the corresponding N. w The first processed data and N w The second processed data, wherein the duration of the second sliding window is shorter than the duration of the first sliding window, and the second percentage is higher than the first percentage, N w Represents a positive integer.

[0087] In step S3, the second sliding window is used to segment the reference cross-correlation function d0(t) and the current cross-correlation function d1(t) respectively, so as to obtain the estimated time delay between the two cross-correlation functions through subsequent frequency domain analysis. The duration of the second sliding window, the overlap rate, and the positive integer N are also considered. wThe choice of N typically depends on the frequency content of the cross-correlation function being considered; in this embodiment, the duration of the second sliding window can be calculated based on the reference seismic wave velocity in the target area (e.g., 200 m / s) and the approximate size of the target area (e.g., 100 m), for example, the duration of the second sliding window is 100 m ÷ 200 m / s ÷ 10 = 50 ms, and the second percentage is, for example, 10%. w The first piece of data and the N w The second data segment corresponds one-to-one in time sequence, the N w The first processed data of the chip and the N w The first piece of data corresponds one-to-one, the N w The second processed data and the N w The second piece of data corresponds one-to-one, therefore the N w The first processed data of the chip and the N w The data processed in the second step will also correspond one-to-one in time sequence. The aforementioned root mean square energy normalization, windowing, and cosine tapering processing includes the following three actions: root mean square energy normalization, windowing (a key technique in signal processing to weight finite-length signals using window functions to reduce spectral leakage and improve the accuracy of frequency domain analysis), and cosine tapering (i.e., Cosine taper, a signal processing technique mainly used to apply cosine window functions at both ends of a signal to reduce edge effects and spectral leakage; this technique is usually used before Fast Fourier Transform, especially when processing finite-length data). The specific processes of these processing methods are existing technologies and will not be elaborated here. In addition, in the aforementioned cosine tapering processing, the cosine tapering parameter can be set to 80% for example.

[0088] S4. For each pair of nodes, apply the Fourier transform to obtain the corresponding N... w The first processed data and N w The frequency domain function corresponding to each piece of processed data in the second processed data is obtained, and the corresponding N is obtained based on the transformation result. w Mutual spectrum.

[0089] In step S4, the cross-power spectrum, short for cross-power density spectrum, is used to describe the statistical correlation between two different signals in the frequency domain. Specifically, for each pair of nodes, the corresponding N is obtained based on the transformation result. w The cross spectrum includes, but is not limited to, the following steps: for a pair of nodes among the plurality of earthquake monitoring nodes, obtain the corresponding nth cross spectrum based on the transformation result. Where n represents less than or equal to N w positive integers, F 0,n(f) represents the frequency domain function corresponding to the nth slice of first processed data in time sequence. F represents the frequency domain function corresponding to the nth slice of second-processed data in time sequence. 1,n The conjugate form of (f), where f represents the frequency. Therefore, N... w The cross spectrum and the N w The second set of data will correspond one-to-one in time sequence. Furthermore, the specific process of the Fourier transform is a current technical method and will not be elaborated upon here.

[0090] S5. For each pair of nodes, based on all corresponding frequency domain functions and N w Based on the cross spectrum and a preset effective frequency range, the energy density is calculated and the corresponding N is obtained. w ×J mutually coherent function values ​​and a weighting coefficient corresponding to each mutually coherent function value, where J represents the total number of sampling frequencies in the effective frequency range.

[0091] In step S5, cross-spectral X is considered. n (f) can be written in the following form:

[0092]

[0093] In the formula, e represents the base of the natural logarithm function, and i represents the imaginary unit. This represents the phase information with frequency f as the independent variable. It also considers that the time delay between two signals can be obtained from the phase of their corresponding cross-spectrums, which has the following linear relationship with frequency:

[0094]

[0095] In the formula, j represents the index of the sampling frequency in the effective frequency range and is a positive integer less than or equal to J, J represents the total number of sampling frequencies in the effective frequency range, and f j k represents the frequency value of the j-th sampling frequency. n Represents cross-spectral X n (f) corresponds to the slope, δt n Represents cross-spectral X n (f) The corresponding estimated delay time. Furthermore, when estimating the cross-spectral time delay, it is necessary to analyze the similarity between the two signals. This similarity can be quantitatively assessed through the cross-coherence function values ​​between their energy densities, and the time delay between the two signals is estimated through the slope of a linear regression. In the linear regression calculation, weighted coefficient alignment can be introduced for weighting. Therefore, this embodiment requires calculating the corresponding N values ​​between the energy densities for each pair of nodes. w×J mutually coherent function values ​​and the corresponding weight coefficients for each mutually coherent function value. Specifically, for each pair of nodes, based on all corresponding frequency domain functions and N... w Based on the cross spectrum and a preset effective frequency range, the energy density is calculated and the corresponding N is obtained. w ×J mutually coherent function values ​​and the corresponding weight coefficients for each mutually coherent function value, including but not limited to: for a certain pair of nodes, based on all corresponding frequency domain functions and N w Based on the cross spectrum and the preset effective frequency range, the energy density is calculated according to the following formula, and the corresponding N is obtained. w ×J mutually coherent function values ​​and the corresponding weight coefficient for each mutually coherent function value:

[0096]

[0097] In the formula, j represents the index of the sampling frequency in the effective frequency range and is a positive integer less than or equal to J, J represents the total number of sampling frequencies in the effective frequency range, and C n,j Indicates the cross spectrum with the nth cross spectrum X n (f) and the cross-coherence function value corresponding to the j-th sampling frequency in the effective frequency range, X n,j This means substituting the j-th sampling frequency into the n-th cross spectrum X. n The function value F obtained after (f) 0,n,j This means substituting the j-th sampling frequency into the frequency domain function F. 0,n The function value F obtained after (f) 1,n,j This means substituting the j-th sampling frequency into the frequency domain function F. 1,n The function value obtained after (f), w n,j Represents the value of the mutually coherent function C. n,j The corresponding weighting coefficients, the effective frequency range is 5–20 Hz. Based on the aforementioned weighting coefficients w n,j The calculation formula shows that the weighting coefficient w n,j By taking into account both the amplitude of the cross spectrum and the strength of their mutual interference, the accuracy of the weight allocation can be ensured.

[0098] S6. For each pair of nodes, based on the corresponding N w Given J weighted coefficients, the corresponding N is obtained using the weighted least squares inversion method. w A slope, and based on that N w The slope is calculated to obtain the corresponding N. w Each delay time estimate, and based on the corresponding N w The corresponding N is obtained by calculating ×J mutually coherent function values. w One correlation coefficient.

[0099] In step S6, the N w The estimated delay time value and the N w Each slope corresponds one-to-one with the N slope. w The slope and the N w Each cross spectrum corresponds one-to-one, therefore the N w The estimated delay time value and the N w The second piece of data will also correspond one-to-one in terms of timing; and the N mentioned above. w The correlation coefficients and the N w Each cross spectrum corresponds one-to-one, therefore the N w The correlation coefficients and the N w The second set of data will also correspond one-to-one in time sequence. Specifically, for each pair of nodes, based on the corresponding N... w Given J weighted coefficients, the corresponding N is obtained using the weighted least squares inversion method. w A slope, and based on that N w The slope is calculated to obtain the corresponding N. w Each delay time estimate, and based on the corresponding N w The corresponding N is obtained by calculating ×J mutually coherent function values. w The correlation coefficient includes, but is not limited to, the following steps S61 to S63.

[0100] S61. For a given pair of nodes, based on the corresponding N w Given ×J weighted coefficients, the corresponding N is obtained using the weighted least squares inversion method according to the following formula. w One slope:

[0101]

[0102] In the formula, k n Indicates the cross spectrum with the nth cross spectrum X n The slope corresponding to (f), f j φ represents the frequency value of the j-th sampling frequency. j This represents the phase information corresponding to the j-th sampling frequency.

[0103] In step S61, φ j Based on Calculated.

[0104] S62. Based on the N of the aforementioned pair of nodes w The slope is used to calculate N for a given pair of nodes according to the following formula. w One estimated delay time:

[0105] δt n =k n ÷(2×π)

[0106] In the formula, δt n Indicates the cross spectrum with the nth cross spectrum X n (f) corresponds to the estimated delay time.

[0107] S63. Based on the N of the aforementioned pair of nodes w Given ×J mutually coherent function values, N for a given pair of nodes is calculated using the following formula. w One correlation coefficient:

[0108]

[0109] In the formula, C n Indicates the cross spectrum with the nth cross spectrum X n (f) Correlation coefficient.

[0110] S7. For each pair of nodes, based on the preset delay time estimation threshold and correlation coefficient threshold, in the corresponding N... w All second data in the second data set that satisfy the condition that the corresponding delay time estimate is less than or equal to the delay time estimate threshold and the corresponding correlation coefficient is greater than or equal to the correlation coefficient threshold are considered as the corresponding valid data.

[0111] In step S7, the valid data is used to ensure accurate acquisition of observed seismic wave velocity changes for the corresponding node pairs. For example, the delay time estimation threshold is 5 milliseconds, and the correlation coefficient threshold is 70%.

[0112] S8. Based on the valid data and geographical location of each pair of nodes, analyze the distribution of seismic wave velocity changes in the underground medium of the target area during the most recent period.

[0113] In step S8, specifically, for each pair of nodes, dt / t (dt represents the differential form with respect to time t) can be routinely analyzed based on the corresponding effective data to obtain the corresponding dv / v, which is the observed value of the corresponding seismic wave velocity change. Then, based on the observed value of the seismic wave velocity change of each pair of nodes and their geographical location, the distribution of the seismic wave velocity change of the underground medium in the target area during the current most recent period can be routinely obtained.

[0114] S9. Trigger a ground collapse early warning action for the target area based on the distribution of the seismic wave velocity changes.

[0115] In step S9, for example, based on the distribution of seismic wave velocity changes, if the observed value of the seismic wave velocity change corresponding to any pair of nodes is found to meet the preset early warning conditions, then a ground collapse early warning action is executed for the target area (e.g., issuing an alarm message instructing the closure of the target area, etc.).

[0116] Therefore, based on the ground subsidence monitoring method described in steps S1 to S9 above, a new scheme for monitoring ground subsidence in a target area based on continuous passive source seismic data is provided. Specifically, after receiving passive source seismic data collected in real time by multiple seismic monitoring nodes, cross-correlation analysis is performed periodically on two corresponding passive source seismic data sets collected in the most recent period and several consecutive periods for each pair of nodes. Based on the analysis results, multiple delay time estimates and correlation coefficients are obtained through time-domain segmentation and correlation processing, and frequency-domain processing and correlation calculations. Then, the corresponding valid data are selected based on the calculation results. Finally, based on the valid data and geographical location of each pair of nodes, the distribution of seismic wave velocity changes in the subsurface medium of the target area in the most recent period is analyzed, and a ground subsidence early warning action is triggered accordingly. This can overcome the shortcomings of existing ground subsidence monitoring schemes, which lack "early diagnosis" features and only issue warnings after ground subsidence has occurred. It effectively improves the accuracy and timeliness of ground subsidence early warnings and has advantages such as continuous monitoring, simplicity, efficiency, and controllable cost, making it easy to apply and promote in practice.

[0117] like Figure 2 As shown, the second aspect of this embodiment provides a virtual device for implementing the ground subsidence monitoring method described in the first aspect, including a seismic data receiving unit, a cross-correlation analysis unit, a data fragmentation processing unit, a data frequency domain processing unit, a first calculation processing unit, a second calculation processing unit, an effective data filtering unit, a wave velocity change analysis unit, and a monitoring and early warning triggering unit that are connected in sequence.

[0118] The earthquake data receiving unit is used to receive passive source earthquake data collected in real time by multiple earthquake monitoring nodes, wherein the multiple earthquake monitoring nodes are discretely arranged in the target area.

[0119] The cross-correlation analysis unit is used to periodically perform cross-correlation analysis on corresponding passive source seismic data acquired in the most recent period for each pair of nodes among the plurality of seismic monitoring nodes, based on a first sliding window with an overlap rate of a first percentage, to obtain N corresponding to the current most recent period. cur A cross-correlation function, and then the N cur The cross-correlation functions are superimposed and normalized to obtain the corresponding current cross-correlation function d1(t), and the cross-correlation function d1(t) is also obtained by superimposing and normalizing the N most recent consecutive cross-correlation functions. ref The corresponding cross-correlation functions for each period are superimposed and normalized to obtain the corresponding reference cross-correlation function d0(t), where N cur Represents a positive integer greater than or equal to 10, and the duration of the current most recent period is N times the duration of the first sliding window. cur Times, Nref t represents a positive integer greater than or equal to 10, and t represents time.

[0120] The data sharding processing unit is used to divide the corresponding reference cross-correlation function d0(t) into N segments for each pair of nodes, based on a second sliding window with an overlap rate of a second percentage. w The first piece of data, and the corresponding current cross-correlation function d1(t) divided into N w The second piece of data, and on the N w First data and the N w Each data point in the second data set undergoes root mean square energy normalization, windowing, and cosine shrinking to obtain the corresponding N. w The first processed data and N w The second processed data, wherein the duration of the second sliding window is shorter than the duration of the first sliding window, and the second percentage is higher than the first percentage, N w Represents positive integers;

[0121] The data frequency domain processing unit is used to apply Fourier transform to each pair of nodes to obtain the corresponding N... w The first processed data and N w The frequency domain function corresponding to each piece of processed data in the second processed data is obtained, and the corresponding N is obtained based on the transformation result. w Mutual spectrum;

[0122] The first computational processing unit is configured to, for each pair of nodes, perform computation based on all corresponding frequency domain functions and N. w Based on the cross spectrum and a preset effective frequency range, the energy density is calculated and the corresponding N is obtained. w ×J mutually coherent function values ​​and a weighting coefficient corresponding to each mutually coherent function value, where J represents the total number of sampling frequencies in the effective frequency range;

[0123] The second calculation processing unit is used to, for each pair of nodes, based on the corresponding N w Given J weighted coefficients, the corresponding N is obtained using the weighted least squares inversion method. w A slope, and based on that N w The slope is calculated to obtain the corresponding N. w Each delay time estimate, and based on the corresponding N w The corresponding N is obtained by calculating ×J mutually coherent function values. w One correlation coefficient;

[0124] The effective data filtering unit is used to, for each pair of nodes, filter the corresponding N nodes based on a preset delay time estimation threshold and a correlation coefficient threshold. wAll second data in the second data set that satisfy the corresponding delay time estimate value being less than or equal to the delay time estimate threshold and the corresponding correlation coefficient being greater than or equal to the correlation coefficient threshold are considered as the corresponding valid data.

[0125] The wave velocity change analysis unit is used to analyze the distribution of seismic wave velocity changes in the subsurface medium of the target area in the current most recent period based on the effective data and geographical location of each pair of nodes.

[0126] The monitoring and early warning triggering unit is used to trigger a ground collapse early warning action for the target area based on the distribution of the seismic wave velocity change.

[0127] The working process, working details and technical effects of the aforementioned device provided in the second aspect of this embodiment can be found in the ground subsidence monitoring method described in the first aspect, and will not be repeated here.

[0128] like Figure 3 As shown, the third aspect of this embodiment provides a computer device for executing the ground subsidence monitoring method as described in the first aspect, including a memory, a processor, and a transceiver connected in sequence. The memory stores a computer program, the transceiver sends and receives messages, and the processor reads the computer program to execute the ground subsidence monitoring method as described in the first aspect. Specifically, the memory may include, but is not limited to, random-access memory (RAM), read-only memory (ROM), flash memory, first-in-first-out (FIFO) memory, and / or first-in-last-out (FILO) memory, etc.; the processor may include, but is not limited to, a microprocessor of the STM32F105 series. Furthermore, the computer device may also include, but is not limited to, a power module, a display screen, and other necessary components.

[0129] The working process, working details and technical effects of the aforementioned computer equipment provided in the third aspect of this embodiment can be found in the ground subsidence monitoring method described in the first aspect, and will not be repeated here.

[0130] This fourth aspect of the embodiment provides a computer-readable storage medium storing instructions comprising the ground subsidence monitoring method as described in the first aspect. Specifically, the computer-readable storage medium stores instructions that, when executed on a computer, perform the ground subsidence monitoring method as described in the first aspect. The computer-readable storage medium refers to a data storage medium, which may include, but is not limited to, floppy disks, optical disks, hard disks, flash memory, USB flash drives, and / or Memory Sticks. The computer may be a general-purpose computer, a special-purpose computer, a computer network, or other programmable devices.

[0131] The working process, working details and technical effects of the aforementioned computer-readable storage medium provided in the fourth aspect of this embodiment can be found in the ground subsidence monitoring method described in the first aspect, and will not be repeated here.

[0132] This fifth aspect of the embodiment provides a computer program product, including a computer program or instructions, which, when executed by a computer, implement the ground subsidence monitoring method as described in the first aspect. The computer may be a general-purpose computer, a special-purpose computer, a computer network, or other programmable device.

[0133] Finally, it should be noted that the above description is merely a preferred embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A ground subsidence monitoring method based on passive source seismic data, characterized in that, include: The system receives passive source seismic data collected in real time by multiple seismic monitoring nodes, wherein the multiple seismic monitoring nodes are discretely arranged in the target area. For each pair of nodes among the multiple seismic monitoring nodes, based on a first sliding window with an overlap rate of a first percentage, cross-correlation analysis is periodically performed on the corresponding two passive source seismic data acquired in the most recent period to obtain N in the most recent period. cur A cross-correlation function, and then the N cur The cross-correlation functions are superimposed and normalized to obtain the corresponding current cross-correlation function d1(t), and the cross-correlation function d1(t) is also obtained by superimposing and normalizing the N most recent consecutive cross-correlation functions. ref The corresponding cross-correlation functions for each period are superimposed and normalized to obtain the corresponding reference cross-correlation function d0(t), where N cur Represents a positive integer greater than or equal to 10, and the duration of the current most recent period is N times the duration of the first sliding window. cur Times, N ref t represents a positive integer greater than or equal to 10, and t represents time. For each pair of nodes, based on a second sliding window with an overlap rate of a second percentage, the corresponding reference cross-correlation function d0(t) is divided into N... w The first piece of data, and the corresponding current cross-correlation function d1(t) divided into N w The second piece of data, and on the N w First data and the N w Each data point in the second data set undergoes root mean square energy normalization, windowing, and cosine shrinking to obtain the corresponding N. w The first processed data and N w The second processed data, wherein the duration of the second sliding window is shorter than the duration of the first sliding window, and the second percentage is higher than the first percentage, N w Represents positive integers; For each pair of nodes, apply the Fourier transform to obtain the corresponding N... w The first processed data and N w The frequency domain function corresponding to each piece of processed data in the second processed data is obtained, and the corresponding N is obtained based on the transformation result. w Mutual spectrum; For each pair of nodes, based on all corresponding frequency domain functions and N w Based on the cross spectrum and a preset effective frequency range, the energy density is calculated and the corresponding N is obtained. w ×J mutually coherent function values ​​and a weighting coefficient corresponding to each mutually coherent function value, where J represents the total number of sampling frequencies in the effective frequency range; For each pair of nodes, based on the corresponding N w Given J weighted coefficients, the corresponding N is obtained using the weighted least squares inversion method. w A slope, and based on that N w The slope is calculated to obtain the corresponding N. w Each delay time estimate, and based on the corresponding N w The corresponding N is obtained by calculating ×J mutually coherent function values. w One correlation coefficient; For each pair of nodes, based on the preset delay time estimation threshold and correlation coefficient threshold, the corresponding N will be... w All second data in the second data set that satisfy the corresponding delay time estimate value being less than or equal to the delay time estimate threshold and the corresponding correlation coefficient being greater than or equal to the correlation coefficient threshold are considered as the corresponding valid data. Based on the valid data and geographical location of each pair of nodes, the distribution of seismic wave velocity changes in the subsurface medium of the target area during the most recent period is analyzed. The ground collapse early warning action for the target area is triggered based on the distribution of changes in seismic wave velocity.

2. The ground subsidence monitoring method according to claim 1, characterized in that, For each pair of nodes among the multiple seismic monitoring nodes, based on a first sliding window with an overlap rate of a first percentage, cross-correlation analysis is periodically performed on the corresponding two passive source seismic data acquired in the most recent period to obtain N in the most recent period. cur The cross-correlation functions include: For a pair of nodes among the multiple earthquake monitoring nodes, periodically acquire the corresponding two passive source earthquake data collected in the most recent period; The two sets of passive source seismic data of a certain pair of nodes are preprocessed to obtain two sets of preprocessed data. The preprocessing includes filtering, anomaly suppression, energy normalization and / or spectral whitening. Based on a first sliding window with an overlap rate of a first percentage, cross-correlation analysis is performed on the two preprocessed data sets to obtain N of the pair of nodes in the current most recent period. cur There are N cross-correlation functions, where N cur Represents a positive integer greater than or equal to 10, and the duration of the current most recent period is N times the duration of the first sliding window. cur times.

3. The ground subsidence monitoring method according to claim 2, characterized in that, When the preprocessing includes filtering, the filtering is used to extract a preferred frequency band signal with a frequency band range of 2.0 to 40 Hz.

4. The ground subsidence monitoring method according to claim 1, characterized in that, For each pair of nodes, the corresponding N is obtained based on the transformation result. w Each cross spectrum includes: For a given pair of nodes among the multiple earthquake monitoring nodes, the corresponding nth cross spectrum is obtained based on the transformation result. Where n represents less than or equal to N w positive integers, F 0,n (f) represents the frequency domain function corresponding to the nth slice of first processed data in time sequence. F represents the frequency domain function corresponding to the nth slice of second-processed data in time sequence. 1,n (f) is the conjugate form of f, where f represents the frequency.

5. The ground subsidence monitoring method according to claim 4, characterized in that, For each pair of nodes, based on all corresponding frequency domain functions and N w Based on the cross spectrum and a preset effective frequency range, the energy density is calculated and the corresponding N is obtained. w ×J mutually coherent function values ​​and the corresponding weight coefficients for each mutually coherent function value, including: For a given pair of nodes, based on all corresponding frequency domain functions and N w Based on the cross spectrum and the preset effective frequency range, the energy density is calculated according to the following formula, and the corresponding N is obtained. w ×J mutually coherent function values ​​and the corresponding weight coefficient for each mutually coherent function value: In the formula, j represents the index of the sampling frequency in the effective frequency range and is a positive integer less than or equal to J, J represents the total number of sampling frequencies in the effective frequency range, and C n,j Indicates the cross spectrum with the nth cross spectrum X n (f) and the cross-coherence function value corresponding to the j-th sampling frequency in the effective frequency range, X n,j This means substituting the j-th sampling frequency into the n-th cross spectrum X. n The function value F obtained after (f) 0,n,j This means substituting the j-th sampling frequency into the frequency domain function F. 0,n The function value F obtained after (f) 1,n,j This means substituting the j-th sampling frequency into the frequency domain function F. 1,n The function value obtained after (f), w n,j Represents the value of the mutually coherent function C. n,j The corresponding weighting coefficients, the effective frequency range is 5 to 20 Hz.

6. The ground subsidence monitoring method according to claim 5, characterized in that, For each pair of nodes, based on the corresponding N w Given J weighted coefficients, the corresponding N is obtained using the weighted least squares inversion method. w A slope, and based on that N w The slope is calculated to obtain the corresponding N. w Each delay time estimate, and based on the corresponding N w The corresponding N is obtained by calculating ×J mutually coherent function values. w The correlation coefficients include: For a given pair of nodes, based on the corresponding N w Given ×J weighted coefficients, the corresponding N is obtained using the weighted least squares inversion method according to the following formula. w One slope: In the formula, k n Indicates the cross spectrum with the nth cross spectrum X n The slope corresponding to (f), f j φ represents the frequency value of the j-th sampling frequency. j This represents the phase information corresponding to the j-th sampling frequency; Based on the N of the aforementioned pair of nodes w The slope is used to calculate N for a given pair of nodes according to the following formula. w One estimated delay time: δt n =k n ÷(2×π) In the formula, δt n Indicates the cross spectrum with the nth cross spectrum X n (f) The estimated delay time corresponding to; Based on the N of the aforementioned pair of nodes w Given ×J mutually coherent function values, N for a given pair of nodes is calculated using the following formula. w One correlation coefficient: In the formula, C n Indicates the cross spectrum with the nth cross spectrum X n (f) Correlation coefficient.

7. A ground subsidence monitoring device based on passive source seismic data, characterized in that, It includes a seismic data receiving unit, a cross-correlation analysis unit, a data fragmentation processing unit, a data frequency domain processing unit, a first calculation processing unit, a second calculation processing unit, an effective data filtering unit, a wave velocity change analysis unit, and a monitoring and early warning triggering unit, which are connected in sequence. The earthquake data receiving unit is used to receive passive source earthquake data collected in real time by multiple earthquake monitoring nodes, wherein the multiple earthquake monitoring nodes are discretely arranged in the target area. The cross-correlation analysis unit is used to periodically perform cross-correlation analysis on two corresponding passive source seismic data acquired in the most recent period for each pair of nodes among the plurality of seismic monitoring nodes, based on a first sliding window with an overlap rate of a first percentage, to obtain N corresponding to the current most recent period. cur A cross-correlation function, and then the N cur The cross-correlation functions are superimposed and normalized to obtain the corresponding current cross-correlation function d1(t), and the cross-correlation function d1(t) is also obtained by superimposing and normalizing the N most recent consecutive cross-correlation functions. ref The corresponding cross-correlation functions for each period are superimposed and normalized to obtain the corresponding reference cross-correlation function d0(t), where N cur Represents a positive integer greater than or equal to 10, and the duration of the current most recent period is N times the duration of the first sliding window. cur Times, N ref t represents a positive integer greater than or equal to 10, and t represents time. The data sharding processing unit is used to divide the corresponding reference cross-correlation function d0(t) into N segments for each pair of nodes, based on a second sliding window with an overlap rate of a second percentage. w The first piece of data, and the corresponding current cross-correlation function d1(t) divided into N w The second piece of data, and on the N w First data and the N w Each data point in the second data set undergoes root mean square energy normalization, windowing, and cosine shrinking to obtain the corresponding N. w The first processed data and N w The second processed data, wherein the duration of the second sliding window is shorter than the duration of the first sliding window, and the second percentage is higher than the first percentage, N w Represents positive integers; The data frequency domain processing unit is used to apply Fourier transform to each pair of nodes to obtain the corresponding N... w The first processed data and N w The frequency domain function corresponding to each piece of processed data in the second processed data is obtained, and the corresponding N is obtained based on the transformation result. w Mutual spectrum; The first computational processing unit is configured to, for each pair of nodes, perform computation based on all corresponding frequency domain functions and N. w Based on the cross spectrum and a preset effective frequency range, the energy density is calculated and the corresponding N is obtained. w ×J mutually coherent function values ​​and a weighting coefficient corresponding to each mutually coherent function value, where J represents the total number of sampling frequencies in the effective frequency range; The second calculation processing unit is used to, for each pair of nodes, based on the corresponding N w Given J weighted coefficients, the corresponding N is obtained using the weighted least squares inversion method. w A slope, and based on that N w The slope is calculated to obtain the corresponding N. w Each delay time estimate, and based on the corresponding N w The corresponding N is obtained by calculating ×J mutually coherent function values. w One correlation coefficient; The effective data filtering unit is used to, for each pair of nodes, filter the corresponding N nodes based on a preset delay time estimation threshold and a correlation coefficient threshold. w All second data in the second data set that satisfy the corresponding delay time estimate value being less than or equal to the delay time estimate threshold and the corresponding correlation coefficient being greater than or equal to the correlation coefficient threshold are considered as the corresponding valid data. The wave velocity change analysis unit is used to analyze the distribution of seismic wave velocity changes in the subsurface medium of the target area in the current most recent period based on the effective data and geographical location of each pair of nodes. The monitoring and early warning triggering unit is used to trigger a ground collapse early warning action for the target area based on the distribution of the seismic wave velocity change.

8. A computer device, characterized in that, The device includes a memory, a processor, and a transceiver connected in sequence, wherein the memory is used to store a computer program, the transceiver is used to send and receive messages, and the processor is used to read the computer program and execute the ground subsidence monitoring method as described in any one of claims 1 to 6.

9. A computer-readable storage medium, characterized in that... The computer-readable storage medium stores instructions that, when executed on a computer, perform the ground subsidence monitoring method as described in any one of claims 1 to 6.

10. A computer program product, comprising a computer program or instructions, characterized in that, When the computer program or the instructions are executed by the computer, they implement the ground subsidence monitoring method as described in any one of claims 1 to 6.

Images