Earthquake comprehensive prediction method

Abstract

The invention relates to the technical field of earthquake monitoring, and discloses a comprehensive earthquake prediction method, which comprises the following steps of: acquiring a continuous micro-motion data sequence to reconstruct an empirical green function sequence, and dividing the empirical green function sequence into a compression phase subset and a stretching phase subset by utilizing a synchronously acquired theoretical earth tide stress phase; respectively extracting a causal branch waveform corresponding to the positive time axis and a non-causal branch waveform corresponding to the negative time axis, and determining a weight index for inhibiting the non-structural interference by calculating a reciprocity residual error between a causal wave velocity offset and a non-causal wave velocity offset; correcting a wave velocity offset mean value by using the weight index to obtain a stress sensitivity coefficient reflecting the stress state of the earth crust medium; according to the method, a physical criterion is established by using an elastic wave propagation reciprocity principle, decoupling of medium attribute change and noise source spatial drift is realized, and visual wave velocity offset induced by environmental interference is eliminated.

Description

Technical Field

[0001] This invention relates to a comprehensive earthquake prediction method, belonging to the field of earthquake monitoring technology. Background Technology

[0002] Currently, the mainstream approach in earthquake monitoring is to reconstruct empirical Green's functions using background noise cross-correlation techniques and then monitor crustal wave velocity perturbations based on Green's function wake waves.

[0003] Shallow crustal media are sensitive to surface environmental loads. Fluctuations in surface air pressure, changes in water level, and non-stationary spatial migration of noise sources generate apparent velocity shifts in monitoring sequences. These environmental loads induce interference signals that cover weak tectonic stress anomalies, making it difficult for observation systems to distinguish between real changes in medium properties and random disturbances in the background field. For example, the method for synthesizing destructive earthquake acceleration waveforms from fixed stations, authorized by CN119165530B, relies on empirical Green's function superposition to fit source parameters. While this scheme compensates for the lack of strong earthquake records, it is limited by existing source model experience and does not address the real-time decoupling of environmental load response and tectonic stress response. Due to the lack of symmetric constraints on wavefield propagation reciprocity, spatial drift of noise sources introduces apparent velocity deviations in the simulation sequence. Linear superposition methods cannot identify the nonlinear instability characteristics of the critical rupture stage of the crustal medium, making it difficult to achieve physical calibration of stress saturation state. Extending the superposition time of the reconstructed waveform smooths random errors, resulting in smoothing loss on the time axis of the warning signal, making it impossible to capture transient abrupt changes before fault rupture. Increasing the density of monitoring sites cannot achieve decoupling of environmental load response and tectonic stress response at the physical mechanism level.

[0004] Therefore, how to construct a demodulation mechanism based on solid tidal load as the physical reference, and use wave field topological constraints to achieve physical isolation between environmental interference and tectonic response, thereby improving the reliability of identifying the risk of crustal instability, has become the technical problem to be solved by this invention. Summary of the Invention

[0005] To address the problems mentioned in the background art, the technical solution of the present invention is as follows: A comprehensive earthquake prediction method, comprising the following steps: Step S101: Obtain continuous micro-motion data sequences from multiple monitoring sites in the seismic monitoring network; Step S102: Perform cross-correlation operation on the continuous micro-motion data sequence to reconstruct the empirical Green's function sequence characterizing the medium features between monitoring sites; Step S103: Simultaneously acquire the theoretical solid tidal stress phase of the monitoring area within the time window corresponding to the continuous micro-motion data sequence; Step S104: Based on the theoretical solid tidal stress phase, the time window is divided into a compression phase interval and a tension phase interval, and the empirical Green's function sequence is divided into a compression phase subset and a tension phase subset. Step S105: Extract the causal branch waveforms corresponding to the positive time axis and the non-causal branch waveforms corresponding to the negative time axis from the empirical Green's functions of the compressed phase subset and the stretched phase subset, respectively. Step S106: Using moving window cross-spectral analysis, calculate the causal wave velocity offset of the causal branch between the compressed phase subset and the stretched phase subset, and calculate the non-causal wave velocity offset of the non-causal branch between the compressed phase subset and the stretched phase subset. Step S107: Calculate the absolute difference between the causal wave velocity offset and the non-causal wave velocity offset, which serves as the reciprocity residual characterizing the degree of noise source drift interference. Step S108: Substitute the reciprocity residual into the preset Gaussian kernel function to generate an effectiveness weight index for suppressing unconstructed wave velocity interference. Step S109: Use the effectiveness weight index to perform weighted correction on the mean values ​​of causal wave velocity offset and non-causal wave velocity offset to obtain the stress sensitivity coefficient that reflects the stress state of the crustal medium. Step S110: Monitor the slope characteristics of the stress sensitivity coefficient evolution over time, and when the slope characteristics continuously exceed the preset linear threshold, output earthquake early warning results in combination with the spatial distribution connectivity characteristics between monitoring sites.

[0006] Preferably, step S109 specifically includes: step S201, synchronously acquiring time-series data of surface atmospheric pressure in the monitoring area; step S202, extracting the environmental load component orthogonal to the theoretical solid tidal stress phase from the surface atmospheric pressure time-series data; step S203, establishing a response operator for the environmental load component to the change in medium wave velocity, and calculating the background contribution value of the environmental load component to the medium wave velocity; step S204, using the background contribution value to perform phase compensation on causal wave velocity offset and non-causal wave velocity offset, so as to eliminate non-structural wave velocity offset caused by atmospheric pressure fluctuations.

[0007] Preferably, the method includes the following steps: Step S301, continuously acquiring causal wave velocity offset sequences over multiple periods to construct a response curve with the theoretical solid tidal stress phase as the independent variable; Step S302, performing a discrete Fourier transform on the response curve to extract higher-order harmonic features relative to the principal frequency of the theoretical solid tidal stress; Step S303, determining the nonlinear distortion index of the medium based on the energy ratio of the higher-order harmonic features relative to the fundamental frequency signal. The calculation formula is as follows: ,in, The nonlinear distortion index of the medium. The amplitude of the fundamental frequency signal. For the first The amplitude of the first harmonic component, The preset highest harmonic order; Step S304, when the nonlinear distortion index of the medium... When the numerical values ​​show increased asymmetry and the energy distribution shifts to higher frequencies, it is determined that the crustal medium has entered a sub-instability stage.

[0008] Preferably, step S102 of reconstructing the empirical Green's function sequence characterizing the medium characteristics between monitoring sites further includes: step S401, calculating the cross-spectral coherence coefficient of the reconstructed Green's function relative to the preset reference waveform for each time period; step S402, establishing a nonlinear weight mapping relationship based on the cross-spectral coherence coefficient, and assigning quality weight indicators to the Green's function components in different time windows; step S403, using the quality weight indicators to perform weighted superposition on the empirical Green's function sequence to suppress transient interference introduced by the non-diffused wave field.

[0009] Preferably, step S110, which monitors the slope characteristics of the stress sensitivity coefficient evolving over time and outputs earthquake early warning results when the slope characteristics continuously exceed a preset linear threshold, further includes: step S501, clustering the paths of dual stations at different azimuth angles in the earthquake monitoring network; step S502, calculating the stress sensitivity coefficients corresponding to the paths in different azimuth angle groups respectively, and constructing a stress sensitivity tensor model; step S503, extracting the polarization direction features in the stress sensitivity tensor model, and identifying the correlation between the polarization direction features and the known fault strikes within the monitoring area.

[0010] Preferably, in step S105, both the causal branch waveform and the non-causal branch waveform are extracted from the tail wave segment of the empirical Green's function, and the start time of the tail wave segment is set to be within the range of 2 to 4 times the arrival time of the direct wave.

[0011] Preferably, after obtaining the continuous micro-motion data sequence in step S101, the method further includes performing filtering processing on the continuous micro-motion data sequence using a bandpass filter of 0.1Hz to 2.0Hz to extract the dominant frequency component modulated by solid moisture stress.

[0012] Preferably, the earthquake early warning result is output by combining the spatial distribution connectivity characteristics between monitoring sites, specifically including: step S801, extracting the stress sensitivity coefficients corresponding to multiple adjacent monitoring site pairs; step S802, verifying the connectivity characteristics of the stress sensitivity coefficients in spatial distribution, and outputting the earthquake early warning result when the stress sensitivity coefficients of multiple adjacent monitoring site pairs simultaneously show a sudden increase and exhibit spatial geometric connectivity.

[0013] Preferably, in step S108, the effectiveness weight index is generated using the following Gaussian kernel formula. : ,in, As an effectiveness weighting indicator, For reciprocal residuals, The preset bandwidth parameter is based on the background wavefield coherence between monitoring sites.

[0014] Preferably, in step S110, when monitoring the slope characteristics, the stress saturation state of the crustal medium is identified by calculating the second derivative of the stress sensitivity coefficient array of the most recent 120 tidal cycles.

[0015] Compared with the prior art, the beneficial effects of the present invention are: 1. In comprehensive earthquake prediction, a response demodulation mechanism based on solid tidal load is established. By utilizing the theoretical tidal stress periodicity and the interaction law of the crustal medium, the empirical Green's function sequence is reconstructed and processed according to the phase of the compression phase and the tension phase. This enables the extracted wave velocity offset characteristic quantity to form a physical correlation with the deep crustal stress state, automatically suppresses non-stationary environmental noise unrelated to the tidal period, and reduces the risk of false alarms caused by random fluctuations in the background field.

[0016] 2. By introducing the reciprocity residual constraint of the causal and attribution branches of the Green's function, and using the physical symmetry of wave field propagation as the signal criterion, the consistency of positive and negative time axis offsets is compared to achieve physical decoupling between changes in medium properties and spatial drift of noise sources. This avoids the apparent wave velocity shift induced by the migration of surface noise source locations in the cross-correlation function, ensuring the physical authenticity of stress sensitivity monitoring indicators in complex source environments. An instability identification mechanism based on the harmonic energy distribution of the response curve is constructed. The high-order harmonic characteristics of the wave velocity response curve relative to the tidal dominant frequency are extracted, the nonlinear distortion index of the medium is calculated, and the physical information of the crustal medium transitioning from the linear elastic stage to the nonlinear interaction stage of microfractures is captured. This enables the calibration of the critical state of the fault sub-instability stage and extends the early warning period for catastrophic rupture events.

[0017] 3. Implement anisotropic scanning of stress sensitivity based on azimuth clustering, fit the differences in stress response along different monitoring paths, reconstruct the stress sensitivity tensor model, extract polarization direction features, establish a correlation criterion between stress accumulation direction and known fault strike, transform risk assessment from single stress amplitude monitoring into a spatially directional fault stability evaluation, and improve the efficiency of identifying rupture risks in specific directions; adopt a dynamic weighted integration mechanism based on cross-spectral coherence to calculate the coherence coefficient between the reconstructed waveform and the reference waveform within a short time window, assign nonlinear weight indices to the Green's function components at different time periods, and use the weight function to suppress non-diffusion wave field interference introduced by human activities or extreme weather, maintain the real-time nature of early warning, and ensure the resilience of long-term monitoring sequences. Attached Figure Description

[0018] Figure 1 This is a schematic diagram of the earthquake comprehensive prediction method based on reciprocity residual correction according to the present invention. Figure 2 This is a schematic diagram illustrating the statistical trend of the calculation error of this invention as a function of the tidal period window length; Figure 3 This is a diagram of the cloud-edge-device layered collaborative earthquake prediction system architecture of the present invention. Detailed Implementation

[0019] The following disclosure is intended to explain and illustrate the present invention, and not to limit the scope of protection of the present invention.

[0020] This invention provides a comprehensive earthquake prediction method. It acquires continuous micro-motion data sequences from a distributed seismic monitoring network and performs cross-correlation calculations to reconstruct an empirical Green's function sequence characterizing the medium characteristics between monitoring sites. Using synchronously acquired theoretical solid tidal stress phases, the empirical Green's function sequence is divided into compressional and tensile subsets. By extracting causal and non-causal waveforms and calculating reciprocity residuals, an effectiveness weight index is generated to correct the mean wave velocity offset, thereby obtaining a stress sensitivity coefficient reflecting the stress state of the crustal medium. Finally, by monitoring the temporal evolution slope characteristics of the stress sensitivity coefficient, the stress saturation state of the crustal medium is identified, and earthquake early warning results are output. Since there is a contradiction between background noise interference and weak precursor signals in earthquake prediction, and the migration of noise source locations can cause a large-scale apparent wave velocity offset, leading to a risk of false alarms in the early warning system, this method acquires continuous micro-motion data sequences from multiple monitoring sites in the seismic monitoring network and utilizes... Hz to A bandpass filter of Hz is used to filter the continuous micro-motion data sequence to extract the dominant frequency component modulated by solid tidal stress. Cross-correlation operation is performed on the processed sequence to reconstruct the empirical Green's function sequence characterizing the medium characteristics between monitoring sites. This procedure uses natural micro-motion signals as seismic sources to obtain Green's function waveforms that reflect the medium properties.

[0021] Based on the fact that spatial migration of non-stationary noise sources can induce spurious phase fluctuations in the Green's function, the system introduces a decoupling mechanism based on wavefield reciprocity symmetry constraints. This mechanism synchronously acquires the theoretical solid tidal stress phase within the time window corresponding to the continuous micro-motion data sequence of the monitoring area. The time window is divided into a compression phase interval and a tension phase interval based on the theoretical solid tidal stress phase. The empirical Green's function sequence is further divided into compression and tension subsets. Causal branch waveforms corresponding to the positive time axis and non-causal branch waveforms corresponding to the negative time axis are extracted from the empirical Green's functions of the compression and tension subsets. Both causal and non-causal branch waveforms are extracted from the tailwave segment of the empirical Green's function, with the start time of the tailwave segment set to the arrival time of the direct wave. Doubled Within the interval, the system uses a moving window cross-spectral analysis method to calculate the causal wave velocity offset of the causal branch between the compressed and stretched phase subsets, and calculates the non-causal wave velocity offset of the non-causal branch between the compressed and stretched phase subsets. The system calculates the absolute difference between the causal and non-causal wave velocity offsets as the reciprocity residual characterizing the degree of noise source drift interference. To suppress unstructured disturbances, the system utilizes the formula Generate effectiveness weight index ,in, As an effectiveness weighting indicator, For reciprocal residuals, This is a preset bandwidth parameter based on the background wavefield coherence between monitoring sites. This weighting index reflects the physical symmetry of the signal. When it increases This rapidly reduces the contribution of noise pollution during periods of noise source migration.

[0022] Obtaining the theoretical solid tidal stress phase At that time, the latitude and longitude coordinates and elevation data of the monitoring site are input into the preliminary reference earth model, and the tidal stress tensor components and bandwidth parameters corresponding to the sampling time are calculated using the physical tidal prediction program. The reciprocity residuals of the statistical monitoring array over a period of no less than 28 calendar days during the background field stationary period are obtained. The sequence distribution characteristics are determined, specifically by calculating the sequence standard deviation and... Set the weight index to 1.5 times the standard deviation. Under stable spatial distribution of noise sources, the system maintains a value above 0.95, corrects the mean wave velocity offset, and reduces the contribution of non-structural background fluctuations. To eliminate non-structural wave velocity offsets caused by atmospheric pressure fluctuations, the system simultaneously acquires time-series data of surface atmospheric pressure in the monitoring area. Extracting time-series data of surface atmospheric pressure A linear response operator for the change in medium wave velocity to the environmental load component is established by considering the environmental load component that is orthogonal to the phase of the theoretical solid tidal stress. The system executes a linear response operator based on multiple linear regression. The online calibration procedure utilizes the latest Time series data of surface atmospheric pressure over a natural day Using the corresponding wave velocity offset sequence as the input variable and the Wiener filtering process as the response variable, the linear response operator is determined by solving the derivative response matrix in the frequency domain using the least squares algorithm. The parameter distribution is such that the system triggers the operator at midnight every day. The sliding update is used to correct the transfer function deviation induced by periodic fluctuations in shallow geological water content, ensuring that the compensation accuracy for the background contribution of atmospheric pressure load is maintained within a certain range. The above calculates the background contribution of environmental load components to the medium wave velocity, and uses this background contribution to perform phase compensation on the causal and non-causal wave velocity offsets. This is time-series data of surface atmospheric pressure. It is a linear response operator.

[0023] When deploying a seismic monitoring network, the physical distance between adjacent monitoring sites Based on the reference wave velocity of the crustal medium With dominant frequency center value Settings, Maintenance No less than three times the dominant wavelength Satisfying the background noise cross-correlation spread wave field assumption, Depend on and The ratio is determined to suppress the effect of near-field non-scattered wave components on the causal wave velocity offset. Non-causal wave velocity offset Computational interference, time series data of surface atmospheric pressure The contribution of the medium wave velocity background is determined by the linear response operator. Compensation, Operator A 14-day sliding window multiple linear regression method is employed, updating the data daily at midnight. The derivative response matrix is ​​solved in the frequency domain and processed using Wiener filtering to eliminate non-tectonic wave velocity shifts caused by periodic fluctuations in shallow geological water content and atmospheric pressure loads. This yields a stress sensitivity coefficient reflecting the stress state of the crustal medium. The system utilizes validity weighting indicators. The stress sensitivity coefficient is obtained by performing a weighted correction on the mean of the causal and non-causal wave velocity offsets. Monitoring stress sensitivity coefficient The slope characteristics that evolve over time are calculated by... The second derivative of the stress sensitivity coefficient array for each tidal cycle identifies the stress saturation state of the crustal medium. When the second derivative continuously exceeds a preset linear threshold and exhibits spatial geometric connectivity among multiple adjacent monitoring point pairs, an earthquake early warning result is output. In the sub-instability stage where the medium is close to rupture, the system executes an instability identification procedure based on response waveform distortion, continuously acquiring wave velocity offset sequences over multiple cycles to construct a response curve with the theoretical solid tidal stress phase as the independent variable. The system then performs a discrete Fourier transform on the response curve to extract higher-order harmonic features relative to the dominant frequency of the theoretical solid tidal stress. Based on the energy ratio of the higher-order harmonic features relative to the fundamental frequency signal, the nonlinear distortion index of the medium is determined. Using the formula Calculate the distortion index, where, The nonlinear distortion index of the medium. The amplitude of the fundamental frequency signal. For the first The amplitude of the first harmonic component, The highest harmonic order is preset, and the nonlinear distortion index of the medium is... When the crustal medium exhibits enhanced asymmetry and a shift in energy distribution towards higher frequencies, it is determined that the medium has entered a sub-instability stage. This procedure captures the critical physical information of the medium's transition from a linear stage to a nonlinear interaction stage.

[0024] Monitoring stress sensitivity coefficient The slope characteristics evolved over time, and the stress sensitivity coefficient for the most recent 120 tidal cycles was calculated. Second derivative of an array Identify the stress saturation state of the crustal medium and determine the linear mechanical threshold of seismic risk. Based on the historical static period Given a sample array, calculate the sample mean. with standard deviation set up This makes the judgment benchmark correspond to the inherent background fluctuation level of the tectonic zone, and the nonlinear distortion index of the medium. According to the formula Calculation, where The amplitude of the fundamental frequency signal. For the first First harmonic component amplitude, To preset the highest order of harmonics, when When the asymmetry of the energy distribution increases and shifts to higher frequencies above 2.0 Hz, the crustal medium is determined to have entered a sub-instability stage, and an earthquake early warning result is output. To improve monitoring stability, the system calculates the cross-spectral coherence coefficient of the reconstructed Green's function relative to the preset reference waveform for each time period. Establish a cross-spectral coherence coefficient The nonlinear weight mapping relationship assigns quality weight indices to the Green's function components within different time windows, and performs weighted superposition to suppress transient interference introduced by the non-diffused wavefield. The coherence coefficient is used to identify the principal stress accumulation direction of the fault rupture by executing a stress sensitivity tensor model construction procedure based on azimuth angle clustering, targeting the directional characteristics of fault rupture. This process groups the paths of dual stations at different azimuth angles within the monitoring array according to... Spatial clustering is performed using step size, and the stress sensitivity coefficient corresponding to each group path is calculated separately. The least squares method is used to fit the observed values ​​at different azimuth angles to a surface in order to solve for the second-order tensor describing the stress state of the medium. The procedure determines the principal stress polarization direction based on the spatial projection of tensor eigenvalues, and identifies the spatial orientation of stress concentration by monitoring the evolution trend of the angle between the principal stress polarization direction and the known fault strike within the monitoring area. This procedure ensures that the early warning results have a definite physical and geometric basis when the dominant tectonic stress direction deviates. The stress sensitivity coefficient, For the stress sensitivity tensor components.

[0025] Example 1: In the scenario of deploying a distributed seismic monitoring network at the edge of a highly dynamic sea area, continuous micro-motion data sequences are controlled by background noise from ocean waves with spatial migration characteristics. The fluctuation of the tides causes changes in the geometric centroid of the noise source, resulting in a magnitude of [insert magnitude here] in the reconstructed empirical Green's function. Apparent wave velocity migration, the values ​​of which cover stress perturbation induced in the deep crust. Variations in wave velocity cause the monitoring system to generate false earthquake risk reports due to non-stationary fluctuations from external signal sources; to eliminate systematic biases introduced by non-stationary noise sources, the system utilizes... Hz to A Hz bandpass filter extracts the dominant frequency components from a continuous micro-motion data sequence. Based on the synchronous calculation theory of solid tidal stress phase, the empirical Green's function sequence is divided into a compressional phase subset and a tensile phase subset. The causal branch wave velocity offset is obtained using the moving window cross-spectral analysis method. Non-causal branch wave velocity offset Calculate the reciprocity residuals This characterizes the degree of wavefield symmetry breaking and applies the reciprocity residual. Substitution based on parameters Calibrate the Gaussian kernel formula to generate an effectiveness weight index ,in, For causal branch wave velocity offset, This is a non-causal branch wave velocity offset. For reciprocal residuals, For preset bandwidth parameters, As an effectiveness weighting index, by reducing the weight of periods polluted by noise source drift, the system reduces non-structural background wave velocity fluctuations to a lower level. the following.

[0026] The system executes a multi-source load phase orthogonal decoupling procedure, extracts the phase orthogonal components of the theoretical solid tidal stress from the surface atmospheric pressure data, establishes a linear response operator, calculates the contribution of environmental loads to the background wave velocity of the medium, and performs phase compensation to obtain the stress sensitivity coefficient reflecting the stress state of the crustal medium. The system monitors the recent Tidal cycle stress sensitivity coefficient The second derivative characteristics of the array were used to observe the nonlinear uplift trend of the crustal medium evolution curve and confirm the nonlinear distortion index of the medium. If the energy distribution exceeds a preset threshold and shifts to higher frequencies, the polarization direction determined by the stress sensitivity tensor model is used to pinpoint the risk area, and a high signal-to-noise ratio earthquake early warning result is output before the mainshock occurs. The stress sensitivity coefficient, The nonlinear distortion index of the medium.

[0027] Example 2: In a controlled experiment used to verify the resolution of tectonic stress identification, data was collected from observation areas deployed in the fault sub-instability zone. A broadband seismograph acquired raw, continuous micro-motion data sequences. The data used in the experiment came from measured waveform data acquired by this array. The frequency response range of the monitoring equipment was [missing information]. Hz to Hz, sampling rate set to Hz, which selects a balance between the analytical accuracy of the waveform time axis and the computational storage load in order to extract Hz to When the dominant frequency component of Hz satisfies the Nyquist sampling theorem and avoids spectral aliasing, the system actively superimposes energy ratio in the input sequence as follows: The system simulates random Gaussian noise in dB and a set of spatially non-stationary wave micro-motion sources to simulate environmental disturbances; it performs cross-correlation calculations to reconstruct the empirical Green's function sequence, completes the division of the compression phase subset and the tension phase subset based on the theoretical solid tidal stress phase, and extracts the causal branch wave velocity offset. Non-causal branch wave velocity offset By calculating reciprocity residuals Determine the effectiveness weight index Preset bandwidth parameters Set as This value is calibrated by normal fitting of the reciprocity residual distribution of the observation array during the background field stationary period, aiming to identify unconstructed apparent velocity shifts. Table 1 shows the comparison results of the monitoring indicators of the sample group of this invention and the control sample group under different noise source shift intensities.

[0028] Table 1: Comparison of monitoring indicators between the present invention sample group and the control sample group under different noise source offset intensities. in For causal branch wave velocity offset, This is a non-causal branch wave velocity offset. For reciprocal residuals, As an effectiveness weighting indicator, The stress sensitivity coefficient is shown in Table 1. When the spatial offset of the noise source is... Increase to At that time, the mean original wave velocity offset corresponding to the comparison sample group showed a monotonically large increase due to the apparent wave velocity effect, causing the tectonic stress signal to be covered by environmental fluctuations. However, in the sample group of this invention, the reciprocity residual increased with the change in wave velocity. cross Boundary, effectiveness weight index The nonlinear decay and tendency to zero of the data from the contaminated period suppresses the contribution of the data to the final mean, thus reducing the stress sensitivity coefficient after correction. It remains constant under different interference intensities. to The experiment achieves physical decoupling between changes in crustal medium properties and spatial drift of noise sources within a stable range; and to address the nonlinear response characteristics of the medium entering the fracturing stage, the experiment incorporates tectonic load intensity to simulate stress saturation processes, and the system calculates the most recent... Stress sensitivity coefficient per tidal cycle Extracting stress sensitivity coefficients from arrays The second derivative characteristics over time, when the tectonic stress level is at the elastic limit of the crustal medium. The following is the stress sensitivity coefficient The evolution slope remains linear and stable when the stress level exceeds... Stress sensitivity coefficient observed after performance inflection point The second derivative is from Increase to Meanwhile, the nonlinear distortion index of the medium The gradient effect is enhanced; experimental data show that as the stress level evolves from the linear region to the saturation region, the nonlinear distortion index of the medium increases. From the benchmark value Increase to And its energy distribution is oriented in the spectrum towards Migration to higher frequency bands above Hz.

[0029] Example 3: This example combines Figures 1 to 3 This section explains a comprehensive earthquake prediction method, such as... Figure 1 As shown, the external entity—the seismic monitoring network—inputs the continuously collected micromotion data sequence into the data preprocessing and filtering module. After bandpass filtering in the 0.1Hz-2.0Hz range, the dominant frequency components are extracted and transmitted to the cross-correlation and phase diversity module. Simultaneously, the external entity—the environment and theoretical model—provides this module with the theoretical solid tidal stress phase to assist in the division of compressive and tensile phases. The generated compressive / tensile phase subsets are then fed into the wave velocity migration and reciprocity analysis module for moving window cross-spectral analysis. This module outputs causal and non-causal wave velocity migrations and reciprocity residuals. The weight calculation and coefficient correction module, combined with bandwidth parameters included in the preset parameter library, is then processed. Parameters provided by linear threshold In addition to the externally input ground atmospheric pressure time series data, the stress sensitivity coefficient S is generated by weighting with a Gaussian kernel function. Finally, the evolution monitoring and early warning judgment module analyzes the coefficient based on slope characteristics and spatial connectivity, generates earthquake early warning results, and sends them to the external entity: the early warning release terminal.

[0030] like Figure 2As shown in the chart, the relationship between the calculation error percentage and the window length and number of tidal cycles is illustrated. The horizontal axis represents the window length range, covering 24 to 240 cycles, and the vertical axis represents the percentage of calculation error. The bar chart with diagonal lines shows that as the window length increases from 24 cycles to 120 cycles, the calculation error decreases and reaches its lowest value at 120 cycles. However, when the window length continues to increase to 240 cycles, the calculation error gradually increases. Figure 3 As shown, the system is divided into a perception edge layer, an intelligent cloud computing layer, and an application decision layer from bottom to top. The perception edge layer includes a wideband seismometer configured for the earthquake monitoring array, an acquisition and transmission unit, and a barometric pressure sensor and a temperature sensor configured for the environmental load sensing array. It is responsible for reporting the raw data stream to the intelligent cloud computing layer. The intelligent cloud computing layer deploys a multi-source data lake to store waveform data, barometric pressure data, solid tide models, and intelligent operation and maintenance modules to perform status monitoring, automatic parameter calibration, log auditing, and core business computing engines. It integrates Green's function reconstruction services, reciprocity residual verification services, and stress sensitivity coefficient calculation services. The early warning signals and analysis results generated by the calculations are transmitted to the application decision layer, and it includes an earthquake early warning release terminal, a fault stress visualization screen, and an emergency response linkage interface.

[0031] Example 4: When a distributed seismic monitoring network is deployed in a compressional zone at the edge of a tectonic plate, the stress sensitivity coefficient of the crustal medium... Influenced by the coupling effect of seasonal meteorological loads and tectonic stress accumulation, the system obtains the evolution characteristics of the medium through a dynamic calibration procedure of observation window length and stress saturation threshold, and calculates the stress sensitivity coefficient. Sequence stationarity gain at different time scales The gain Defined as the ratio of stress response signal energy to background fluctuation variance, in areas covering at least Systematic traversal under the physical context of a complete tidal cycle One tidal cycle to The value range of each tidal cycle is identified when the length is set to... One tidal cycle Tidal Difference and Tidal energy produces a gain superposition in cross-spectral coherence analysis, affecting the stress sensitivity coefficient. The calculation error is minimized, and the system determines accordingly. Each tidal period serves as the reference window length for reconstructing the second derivative features.

[0032] To determine the threshold for identifying the stress saturation state of the crustal medium, the system utilizes an adaptive calibration procedure based on residual probability distribution to retrieve the second derivative of the stress sensitivity coefficient of the monitored area during historical tectonic static stability periods from memory. Sample array, calculate the average of the sample array. with standard deviation According to the formula Determine the mechanical linear threshold The system will use a mechanical linear threshold. Anchored to the inherent background fluctuation level of this structural region, when the second derivative of the stress sensitivity coefficient is measured in real time... continuous All sampling points exceeded the mechanical linear threshold. Furthermore, the nonlinear distortion index of the medium When the crustal medium exhibits a synchronous trend of migrating to higher frequencies, it is determined that the crustal medium has entered a sub-instability stage of stress saturation from the linear elastic stage. The stress sensitivity coefficient, It is the second derivative. For mechanical linear threshold, The mean, Standard deviation The nonlinear distortion index is used for the medium; when performing phase orthogonal decoupling of multi-source loads, the system acquires time-series data of surface atmospheric pressure in the monitoring area. By using the Hilbert transform to extract the orthogonal components of the theoretical solid tidal stress phase, a linear response operator for the medium wave velocity offset to the orthogonal components is established. By calculating the background contribution value and removing atmospheric pressure phase interference, a stress sensitivity coefficient with a high signal-to-noise ratio can be obtained. sequence, This is time-series data of surface atmospheric pressure. For linear response operators, This is the stress sensitivity coefficient.

[0033] Example 5: When a new monitoring array is connected to the regional network and initialization is performed, the system establishes a monitoring benchmark through a site feature calibration procedure based on historical wavefield statistical distribution, and obtains the monitoring site's position during the background field's static and stable period. Hourly continuous micro-motion data sequences were used to determine the benchmark power spectral density reflecting the background noise level of the local geological environment using a probability distribution fitting method. The system will compare the real-time measured noise power with the reference power spectral density. The comparison confirms that the deviation is within the preset range. Within the tolerance range, time-series data of medium wave velocity offset versus ground atmospheric pressure are established accordingly. The initial linear response matrix is ​​obtained, thus providing a spatially targeted response operator for phase orthogonal decoupling of multi-source loads. As the reference power spectral density, This is time-series data of surface atmospheric pressure.

[0034] When a ground-to-atmospheric pressure sensor link fails or a signal is lost at a monitoring site, the system switches to autonomous constrained monitoring mode using a pre-stored medium elasticity model in its memory, and extracts the most recent data from that monitoring site. A monthly average response coefficient array for each natural day is used to generate a virtual load sequence to replace the real-time measured values ​​using a linear interpolation algorithm. The system utilizes virtual payload sequences A correlation scan was performed with the theoretical solid tidal stress phase, and the stress sensitivity coefficient was maintained after removing the background contribution of environmental interference. To ensure computational continuity, this procedure guarantees that earthquake early warning results will be produced at a frequency no less than once per tidal cycle even when sensor data is interrupted. Secondly, the identification error of the second derivative of the crustal medium under stress saturation remains at... within, For virtual payload sequence, This is the stress sensitivity coefficient.

[0035] Example 6: In a scenario where a seismic monitoring array is deployed in a geological structural zone with scattering characteristics, the system executes a station spacing calibration procedure based on wave theory to obtain the reference wave velocity of the medium within the monitoring area. Combined with the center value of the dominant frequency after filtering Through formula Calculate the dominant wavelength of the micro-motion signal The physical distance between the two monitoring sites Set at no less than Within the range, this geometric constraint procedure is used to satisfy the diffused wavefield assumption in the empirical Green's function reconstruction process, aiming to suppress the influence of near-field non-scattered wave components on causal wave velocity offsets. Non-causal wave velocity offset Interference in the calculation process, For reference wave speed, The dominant frequency center value, The main wavelength, For physical spacing, For causal branch wave velocity offset, This is the non-causal branch wave velocity offset.

[0036] When the system performs a continuous monitoring task based on empirical Green's function tail waves in the active region, it executes a time resolution optimization procedure for the moving window, increasing the width of the moving window. Set as the theoretical solid tidal period to The ratio range is set by... Window overlap rate and stress sensitivity coefficient The sliding update monitors the coherence coefficient generated during the calculation process. State, when the coherence coefficient is observed Below the threshold When adjusting the reconstruction of the empirical Green's function sequence under the operating conditions, the duration of the micro-motion data superposition is increased until the coherence coefficient is reached. Rebound to the preset level.

[0037] It will be apparent to those skilled in the art that the present invention is not limited to the details of the exemplary embodiments described above, and that the present invention can be implemented in other specific forms without departing from the spirit or essential characteristics of the present invention.

[0038] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention.

Claims

1. A method of earthquake comprehensive prediction, characterized in that, The method comprises the following steps: Step S101, acquiring continuous microseismic data sequences of multiple monitoring sites in a seismic monitoring network; Step S102, performing cross-correlation operation on the continuous microseismic data sequences to reconstruct an empirical Green function sequence representing medium characteristics between the monitoring sites; Step S103, synchronously acquiring theoretical earth tide stress phases of the monitoring region within a time window corresponding to the continuous microseismic data sequences; Step S104, dividing the time window into compression phase intervals and stretching phase intervals according to the theoretical earth tide stress phases, and dividing the empirical Green function sequence into compression phase subsets and stretching phase subsets; Step S105, extracting causal branch waveforms corresponding to the positive time axis and acausal branch waveforms corresponding to the negative time axis from the empirical Green functions in the compression phase subsets and the stretching phase subsets, respectively; Step S106, calculating causal wave velocity offset of the causal branch between the compression phase subsets and the stretching phase subsets by using moving window cross-spectrum analysis, and calculating acausal wave velocity offset of the acausal branch between the compression phase subsets and the stretching phase subsets; Step S107, calculating an absolute difference value between the causal wave velocity offset and the acausal wave velocity offset as a reciprocity residual representing the degree of noise source drift interference; Step S108, substituting the reciprocity residual into a preset Gaussian kernel function to generate an effectiveness weight index for suppressing non-tectonic wave velocity interference; Step S109, performing weighted correction on the mean values of the causal wave velocity offset and the acausal wave velocity offset by using the effectiveness weight index to obtain a stress-sensitive coefficient reflecting the stress state of the crust medium; Step S110, monitoring the slope characteristics of the stress-sensitive coefficient over time, and outputting a seismic warning result when the slope characteristics continuously exceed a preset linear threshold combined with the spatial distribution connectivity characteristics between the monitoring sites.

2. The method according to claim 1, wherein, Step S109 specifically comprises: Step S201, synchronously acquiring ground atmospheric pressure time series data of the monitoring region; Step S202, extracting an environmental load component in the ground atmospheric pressure time series data orthogonal to the theoretical earth tide stress phase; Step S203, establishing a response operator of medium wave velocity change to the environmental load component to calculate a background contribution value of the environmental load component to the medium wave velocity; and Step S204, performing phase compensation on the causal wave velocity offset and the acausal wave velocity offset by using the background contribution value to eliminate non-tectonic wave velocity offset caused by atmospheric pressure fluctuation.

3. The method according to claim 1, characterized in that, The step S102 of reconstructing the empirical Green function sequence representing medium characteristics between the monitoring sites further comprises: Step S401, calculating a cross-spectrum coherence coefficient of the Green function reconstructed in each time period relative to a preset reference waveform; Step S402, establishing a nonlinear weight mapping relationship based on the cross-spectrum coherence coefficient to assign a quality weight index to the Green function components in different time windows; and Step S403, performing weighted superposition on the empirical Green function sequence by using the quality weight index. Step S301, continuously acquire the causal wave velocity offset sequence in multiple periods, and construct a response curve with the theoretical solid tide stress phase as the independent variable; step S302, perform discrete Fourier transform on the response curve, and extract high-order harmonic characteristics relative to the theoretical solid tide stress main frequency; step S303, determine the medium nonlinear distortion index according to the energy proportion of the high-order harmonic characteristics relative to the base frequency signal , and the calculation formula is: , wherein, is the medium nonlinear distortion index, is the base frequency signal amplitude, is the amplitude of the mth harmonic component, is the amplitude of the mth harmonic component, is the preset highest order of the harmonic; and step S304, when the numerical value of the medium nonlinear distortion index presents asymmetric enhancement and energy distribution migration to the high frequency band, it is determined that the crust medium enters the sub-instability stage.

4. The method according to claim 1, wherein, ​ 5. The method according to claim 1, wherein, The step S110 of monitoring the slope feature of the stress sensitivity coefficient evolution over time and outputting the earthquake warning result when the slope feature continuously exceeds the preset linear threshold further comprises: a step S501 of clustering and grouping the double-station paths of different azimuth angles in the earthquake monitoring network; a step S502 of calculating the corresponding stress sensitivity coefficients under different azimuth angle grouping paths respectively, and constructing a stress sensitivity tensor model; and a step S503 of extracting the polarization direction feature in the stress sensitivity tensor model, and identifying the correlation degree of the polarization direction feature relative to the strike of the known fault in the monitoring area.

6. The method according to claim 1, wherein, In the step S105, the causal branch wave and the acausal branch wave are both intercepted from the coda paragraph of the empirical Green function, and the starting time of the coda paragraph is set to be within the interval of 2 to 4 times the direct wave arrival time.

7. The method according to claim 1, wherein, After obtaining the continuous microseismic data sequence in the step S101, further comprising performing filtering processing on the continuous microseismic data sequence by using a band-pass filter of 0.1 Hz to 2.0 Hz, so as to extract the dominant frequency component modulated by the solid tide stress.

8. The method according to claim 1, wherein, The earthquake warning result is output in combination with the spatial distribution connectivity feature of the monitoring sites, specifically comprising: a step S801 of extracting the stress sensitivity coefficients corresponding to a plurality of adjacent monitoring site pairs; a step S802 of verifying the spatial distribution connectivity feature of the stress sensitivity coefficients, and outputting the earthquake warning result when the stress sensitivity coefficients of the plurality of adjacent monitoring site pairs synchronously appear sudden growth and present spatial geometric connectivity.

9. The method according to claim 1, wherein, In step S108, the validity weight indicator is generated by using the following Gaussian kernel formula : wherein, is the validity weight indicator, is the reciprocity residual, is a preset bandwidth parameter based on the background wave field coherence between monitoring sites.

10. The method of claim 1, wherein, In the step S110 of monitoring the slope feature, the stress saturation state of the crust medium is identified by calculating the second derivative of the stress sensitivity coefficient array of the last 120 tidal periods.

Images