Earthquake continuous monitoring and predicting system
By building a continuous earthquake monitoring and prediction system, using multi-dimensional signal processing and geological structure inversion technology, combined with artificial intelligence analysis, the earthquake precursor problem of difficult to extract deep seismic wave signals is solved, and accurate prediction and early warning of earthquakes is achieved.
Patent Information
- Application Number
- CN202311768055.2
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2023-12-20
- Publication Date
- 2025-07-08
AI Technical Summary
Existing seismic monitoring methods cannot effectively detect deep underground geological activities, resulting in inaccurate earthquake prediction. Especially when deep underground seismic wave signals are small and susceptible to interference, it is difficult to extract precursor information of earthquakes.
By building a continuous earthquake monitoring and prediction system, seismic monitoring beacons are used to collect multi-dimensional physical information, perform noise reduction preprocessing and grid signal calculation, combine geological structure inversion and artificial intelligence analysis, extract the characteristics of earthquake precursor signals, and use the seismic prediction model for visual interpretation.
It realizes fine monitoring and prediction of deep seismic wave signals, improves the accuracy of earthquake prediction, provides more accurate early warning information for earthquake occurrence, and supports rescue work.
Smart Images

Figure CN120276019A_ABST
Abstract
Description
Technical Field
[0001] The present invention relates to the technical field of earthquake prediction, and more specifically, to an earthquake continuous monitoring and prediction system. Background Art
[0002] Most of the earthquakes harmful to humans occur in the deeper underground at a depth of 10 - 30 km, and currently, humans have insufficient ability to detect the deeper underground greater than 10 km. In recent years, with the progress of ground wave detection technology, it is possible to judge the trend of plate movement by capturing tiny vibration signals in the deeper underground.
[0003] The occurrence of a destructive earthquake is often the sudden outbreak of stress accumulation with a long time span. The plate movement situations around the world are different, yet interrelated and interacting. For each earthquake zone, there are its unique causal relationships. In the past, due to the inability to effectively sense the detailed geological activities in the deeper underground, only macro - analysis, the time difference between the propagation of P - waves and S - waves can be used to give a short - term warning of an earthquake. Or try to analyze the plate energy accumulation cycle by combining historical earthquake data that has occurred, or estimate the possibility of subsequent earthquakes through seismic wave signals in the same region with higher similarity. The main reason why humans cannot accurately predict the occurrence of earthquakes is the lack of understanding of the root cause mechanism of earthquake formation and the deep geological movement.
[0004] The earth's core thermal energy release from the earth's core is conducted through the mantle to the crust for heat dissipation. Through the flow of mantle materials, the crustal plates are continuously pushed to move without stopping. At the macroscopic level, geological movements are occurring all the time, and earthquakes are also occurring all the time. The seismic wave records monitored on the earth's surface are noisy and diverse. These repetitive, phase - changing, and decaying waveforms are likely to cover the tiny but important fracture information in the deep stress area, making it difficult to analyze and interpret in depth. If only relying on historical records, surface measurements, and comparison of periodic similar seismic waveforms to estimate the occurrence of earthquakes, it is very difficult and has a certain degree of uncertainty. Some research institutions have demonstrated the feasibility of an artificial intelligence system for monitoring and warning of earthquake occurrence by constructing a deep - learning earthquake model. The main large - scale earthquake records in China over thousands of years and the full - process records of earthquake information in the past century, as the basic data for artificial intelligence training, will play a positive role, which also creates the possibility of predicting earthquakes through a complex artificial intelligence system in the future.
[0005] If we observe geological movements over a time span of thousands of years, the subduction and compression uplift of continental plates are obvious and intense. If we consider the mantle as a fluid, with the hot and cold mixing at the Mohorovicic discontinuity, the Earth's crust is like an eggshell that is constantly crushed, broken, and then fused. The crust is cold and brittle, and at the boundary between the crust and the mantle, it is the place where hot and cold meet, and the strata mix and rub violently. The deformation of large areas of rock layers caused by plate movements is the source of stress for the most destructive shallow earthquakes. Similarly, if we observe over a time span of days, the plate movements in the world's major earthquake-prone areas can be considered gentle. Earthquakes occur based on basic plate movements. Regional earthquake prediction should be based on the geological characteristics and macroscopic movement trends of the region to demarcate stress blocks, combined with the release degree of the stress area, and conduct verification monitoring for the stress accumulation and release stages of the stress blocks. Especially in areas with a large population in high earthquake-risk zones, more proactive detection means should be adopted to gain a detailed understanding of potential earthquake blocks.
[0006] Existing earthquake monitoring methods often only include incomplete information such as P-waves and S-waves with relatively shallow depths, and there are few other information dimensions available for earthquake monitoring. At the same time, the ability to resist interference signals for weak signals is insufficient, making it difficult to match or effectively extract the information signals of earthquake precursors.
[0007] Based on the known information about the underground stress structure points, the macroscopic layout of geological faults, and past earthquake data in a certain place, through continuous monitoring of ground waves and other physical conditions, by constructing a geological structure fluid movement model and a stress block evolution stage model, conducting artificial intelligence data analysis, determining the stress release stage, extracting a series of information change characteristics before an earthquake occurs, and monitoring these characteristics, it is possible to achieve relatively accurate earthquake prediction. In view of this, we propose a continuous earthquake monitoring and prediction system. Summary of the Invention
[0008] The purpose of the present invention is to provide a continuous earthquake monitoring and prediction system to solve the problems raised in the above background technology.
[0009] To achieve the above purpose, the present invention provides a continuous earthquake monitoring and prediction system, including the following steps:
[0010] S1: The data center receives real-time signals transmitted by earthquake monitoring beacons and obtains multi-dimensional physical information on stratum changes;
[0011] S2: The data processing center performs noise reduction preprocessing on the information;
[0012] S3: Grid the multi-beacon signals and cross-calculate to determine the source points of underground vibration signals and the corresponding fluid patterns of the geology;
[0013] S4: The data processing center performs continuous geological inversion of the geological structure by obtaining ground wave and other multi-dimensional physical change information;
[0014] S5: Match the series of earthquake precursor signals or the initial signal of the shock wave with a threshold with the earthquake omen signal;
[0015] S6: Mark and give an alarm to the signals with a higher degree of matching in the geological inversion system, and then interpret them using the visualization of the earthquake prediction model.
[0016] Preferably, in S1, the earthquake monitoring beacon collects multi-dimensional physical transformation information through a signal processing algorithm.
[0017] Preferably, the signal processing algorithm is specifically as follows:
[0018] Add two groups of white noise n a (t) and -n a (t) with equal numerical magnitudes and opposite signs to the original signal x(t), and obtain:
[0019]
[0020] where n i (t) represents the a-th additive Gaussian white noise sequence, and x1(t) and x2(t) represent the noisy signals after adding positive and negative white noise for the a-th time;
[0021] Perform EMD decomposition on x1(t) and x2(t) to obtain IMF components x1(t) and x2(t):
[0022]
[0023] where M a,b1 and M a,b2 are respectively the b-th IMF component decomposed after adding positive and negative Gaussian white noise for the a-th time, c a,b1 (t) and c a,b2 (t) are the residual components of the EMD decomposition, and B is the number of IMFs;
[0024] Repeat the above steps A times, perform ensemble average operation on the corresponding IMFs above, and the CEEMD decomposition obtains the b-th IMF component c b as:
[0025]
[0026] where M ab is the b-th IMF component decomposed after adding Gaussian white noise for the a-th time;
[0027] The signal s(t) output after CEEMD decomposition is also processed by a non - linear bistable system to output an enhanced signal:
[0028]
[0029]
[0030]
[0031]
[0032]
[0033] r4 = 2l(u(y n + r3)- v(y n + r3) 3 + s n+1 )
[0034] where y(t) is the output enhanced signal, s(t) is the signal output after CEEMD decomposition, u and v are system parameters, is additive white Gaussian noise with a mean of 0 and a variance of σ 2 , y n (t) and y n+1 (t) are the n - th sampling value and the (n + 1)-th sampling value of the output enhanced signal respectively, s n and s n+1 are the n - th sampling value and the (n + 1)-th sampling value of the signal output after CEEMD decomposition respectively, r1, r2, r3 and r4 are all substitution values, and l is the integration step size.
[0035] Preferably, stress plates and stress regions are constructed based on the results of geological inversion in S4 to build a geological model. The main process of geological inversion is based on minimizing the difference between seismic actual data and the geological model, that is, reducing the residual between observed data and model prediction by iteratively updating the geological model.
[0036] Preferably, the process of judging the stress accumulation in the stress region is as follows:
[0037] 1), Data pre - processing and feature extraction:
[0038] The original stress data σ is standardized:
[0039] where μσ and σσ are the mean and standard deviation of the stress data respectively. Features such as the peak value, mean value, variance of the stress are extracted and used as the input of the model;
[0040] 2), LSTM model construction:
[0041] Input layer: Let \(I(t)\) be the stress characteristic data at time \(t\);
[0042] LSTM layer: The following iterative formula is adopted:
[0043] \(f_t=\sigma(W_f\cdot[H_{t - 1},I(t)]+b_f)\)
[0044] \(i_t=\sigma(W_i\cdot[H_{t - 1},I(t)]+b_i)\)
[0045]
[0046]
[0047] \(o_t=\sigma(W_o\cdot[H_{t - 1},I(t)]+b_o)\)
[0048] \(H_t = o_t\times\tanh(C_t)\)
[0049] Output layer: Let \(O(t)\) be the output at time \(t\), representing the predicted stress accumulation state;
[0050] 3), Loss function and optimization:
[0051] Use the training data to optimize the weight and bias parameters of the model, use the mean square error as the loss function, and adopt an optimization algorithm for parameter update:
[0052]
[0053] where \(Y(i)\) is the true stress accumulation state at time \(i\);
[0054] Use the optimizer: Initial learning rate: \(\alpha\), estimated first moment: \(m_t\), estimated second moment: \(v_t\), decay rates: \(\beta_1,\beta_2\). Update rule:
[0055] \(m_t=\beta_1m_{t - 1}+(1 - \beta_1)g_t\)
[0056] \(v_t=\beta_2v_{t - 1}+(1 - \beta_2)g_t^2\)
[0057]
[0058]
[0059]
[0060] where \(\theta\) is the parameter to be updated and \(g_t\) is the gradient;
[0061] 4), Model validation and evaluation:
[0062] Check the performance of the model through a separate validation dataset to prevent overfitting and find the best model and hyperparameters;
[0063] Calculate the evaluation metrics of the model to quantify the performance of the model on the validation set;
[0064]
[0065] If the performance of the model on the validation dataset is poor, it is necessary to return to the model construction stage for adjustment in combination with the actual situation;
[0066] 5), Result application and geological inversion integration:
[0067] Use the trained LSTM model to analyze stress data and obtain relevant features of the predicted stress accumulation value;
[0068] Analyze the predicted stress accumulation state, stress peak, and stress mean value, and analyze them in time and space to better understand their distribution and evolution trend. Use the stress-related features obtained from the prediction and analysis to strengthen the input feature layer of the geological inversion algorithm, use the predicted stress value to correct or constrain the update rule of the geological parameter model, or use it as additional information to guide the inversion in the subsequent inversion process.
[0069] Preferably, the system can also estimate and judge the plate convergence area, movement trend, and stress accumulation stage of the formation based on macroscopic geological conditions, low-frequency microseismic data, surface displacement, and plate movement history data; among them, the composition of the stress accumulation stage label algorithm is as follows:
[0070] 1), Input data:
[0071] Let Xt represent the input data set at time t, where: σt is the stress data, Et and νt are the elastic modulus and Poisson's ratio of the rock, and dt is the plate displacement estimated based on surface displacement measurement and plate movement history data;
[0072] Xt = [σt, Et, νt, dt]
[0073] 2), GRU model:
[0074] Update gate and reset gate:
[0075] zt = σ(Wz × Xt + Uz × ht-1 + bz)
[0076] rt = σ(Wr × Xt + Ur × ht-1 + br)
[0077] Among them:
[0078] Wz and Wr are the weights of the input data to the update and reset gates, Uz and Ur are the weights of the previous hidden state to the update and reset gates, bz and br are the biases, and σ is the activation function;
[0079] 3), Candidate hidden layer:
[0080]
[0081] Where: W and U are weights, b is the bias, and ⊙ represents element-wise multiplication;
[0082] Updated hidden state:
[0083] ht = (1 - zt) ⊙ ht-1 + zt ⊙ ht
[0084] Output layer:
[0085] 1). Degree of stress accumulation: y stress = σ(Ws × ht + bs), where: Ws is the weight and bs is the bias.
[0086] 2). Accumulation stage of stress: y phase = softmax(Wp × ht + bp), where: Wp is the weight and bp is the bias;
[0087] 4), Loss function:
[0088] Let be the actual label of the degree of stress accumulation, be the actual label of the stress accumulation stage, and the loss function L is:
[0089]
[0090] Where: MSE is the mean squared error and CrossEntropy is the cross-entropy loss.
[0091] Preferably, the geological inversion process in S4 belongs to a preferred method after the stress block is delimited, and its specific process is as follows:
[0092] 1) Initial model construction:
[0093] According to the existing geological information and vibration data, initialize the geological velocity model;
[0094] Use the ray tracing algorithm to simulate the propagation of seismic waves in the initial model and determine its main stress area or boundary;
[0095] 2) Inversion iteration:
[0096] Set the discrete grid of space and time:
[0097] xi = i·Δx, yj = j·Δy, zk = k·Δz, tn = n·Δt
[0098] The derivatives of time and space are approximated using finite difference operators;
[0099]
[0100] where Γ is the absorption coefficient, representing the energy loss in wave propagation, and s is the source term;
[0101] In each time step, the above discrete equation is used to update the pressure wave field;
[0102] 3). Residual calculation:
[0103] Calculate the residual between the observed data and the forward model: R = Dobs - Dcalc, where Dobs is the observed data and Dcalc is the data calculated from the current model;
[0104] 4). Model update:
[0105] Using inverse problem techniques, update the model based on the residual, and the update formula is:
[0106]
[0107] where m new and m old are the model parameters after and before the update respectively, α is the step size, is the residual gradient with respect to the model parameters;
[0108] 5). Convergence check and analysis of geological inversion results:
[0109] Check whether the magnitude of the model update or the reduction of the residual meets the predefined stopping criterion;
[0110] Analyze the accuracy and reliability of the geological inversion results. Uncertainty quantification and resolution analysis are required. Quantifying uncertainty and resolving resolution are key components of the inverse problem solution;
[0111] Select a prior distribution P(m) for the model parameter m, which contains the determined formation information and the range of ground wave spectra. Combining the likelihood function P(Dobs∣m) that describes the probability of the observed data Dobs given the model parameter m, use Bayes' theorem to calculate the posterior distribution:
[0112]
[0113] The posterior distribution P(m∣Dobs) provides the uncertainty information of the updated parameters;
[0114] To ensure the convergence and stability of the algorithm, an appropriate optimization method can be selected, and the step size and convergence criterion can be carefully chosen. The step size α can be dynamically selected through line search or trust region methods;
[0115] The convergence criterion can be based on the reduction of the residual or the magnitude of parameter updates. When implementing, first select an appropriate Bayesian framework and optimization algorithm according to prior information and data errors, and then perform the inversion process by iteratively updating parameters, evaluating the posterior distribution, and calculating the resolution matrix. In each iteration, it is necessary to check the convergence and stability of the algorithm and may need to adjust the optimization parameters and convergence criterion;
[0116] The data processing center in S4 performs inversion of geological structures based on multi-dimensional physical change information obtained from the processing results of signal processing algorithms. In some areas that require detailed investigation, mobile artificial vibration signal sources can be set to detect the formation, and reflection seismic tomography can be used. That is, a ground vibration signal with a determined energy is generated on the ground. When these ground waves encounter faults or structural blocks, they will produce reflections and diffractions, and the distributed beacon grid will receive a series of signals from underground. Under the condition of quantitative movement of the ground signal and multiple signal verifications, the underground formation structure can be detected more accurately.
[0117] Preferably, after the basic formation and stress blocks are determined, the continuous low-frequency vibration signals in S5 are manifested as extrusion friction, rock layer fracture, and intermittent rock deformation vibration signals.
[0118] The series of earthquake precursor signals in S5 are the extraction markers of the regular characteristics of the vibration signals related to the formation movement of the plot before a major earthquake. In a complete earthquake signal record, except for the energy release with strong displacement and vibration such as the main shock, the continuous tiny vibration signals record the development changes and movement trends of its stress area. This time span is relatively long. Using the long attention algorithm of artificial intelligence, the movement of the formation plot is calculated by association, and the characteristics of the series of earthquake signals before the earthquake are extracted under similar geological types, providing support for subsequent earthquake prediction;
[0119] The formation displacement algorithm based on the ground displacement measurement of seismic monitoring beacons and underground inertial measurement data is composed as follows:
[0120] 1) Calculate the velocity v_imu(t) and displacement s_imu(t) from the acceleration a_imu_filtered(t):
[0121] V_imu(t) = ∫₀ᵗ a_imu_filtered(τ)dτ + C1
[0122] S_imu(t) = ∫₀ᵗ v_imu(τ)dτ + C2
[0123] Among them, C1 and C2 are integration constants, which can be determined according to the initial values. For multiple IMUs, the average value of their displacement data is: s_imuavg(t) = (1 / n)∑i = 1n s_imu_i(t), where n is the number of IMUs;
[0124] 2) Displacement ratio comparison and analysis
[0125] Let the surface GPS (or Beidou) displacement be s_gps(t), and compare s_imuavg(t) and s_gps(t) to obtain the displacement difference δ(t): δ(t) = s_imuavg(t) - s_gps(t);
[0126] 3) Subsurface stress block displacement judgment
[0127] For the IMU data s_imu_i(t) at different depths, where i represents different depths, calculate the difference Δs_i(t) between them: Δs_i(t) = s_imu_i(t) - s_imuavg(t)
[0128] If all Δs_i(t) are less than a certain threshold, it is determined that the displacement of the subsurface measurement point to the surface area formation block is uniform displacement; otherwise, there is a non-uniform moving layer from the subsurface formation to the surface.
[0129] Preferably, during the continuous monitoring of the earthquake by the seismic beacon, small earthquakes or extremely small subsurface friction displacement signals will continuously occur in the monitoring area. Extract features from the vibration signals before these earthquakes, and at the same time extract some pre-earthquake features from similar terrains or recorded large earthquakes. After extracting the relevant features, monitor these signal features. When the degree of feature matching detected is relatively high, it is predicted that an earthquake will occur. Use algorithms such as deep learning to verify with subsequent facts, including the prediction and verification of aftershocks after the main earthquake, and continuously improve the prediction accuracy;
[0130] Among them, the feature extraction algorithm process is as follows:
[0131] 1) Input the historical s_historical(t) and the seismic signal time series s_current(t) of the current monitoring area, and output: the feature vector F.
[0132] 2) Extraction steps:
[0133] Use the Short-Time Fourier Transform (STFT) or Wavelet Transform on s_historical(t) and s_current(t) to obtain the time-frequency representation of the signal;
[0134] S_current(f,t) = STFT(s_current(t))
[0135] Shistorical(f,t) = STFT(shistorical(t))
[0136] Considering the importance of signal correlation over a long time span, especially in the context of earthquake prediction, the autocorrelation and cross-correlation characteristics of signals can help identify regular patterns and potential cyclic behaviors over a long time span:
[0137] (1) Extract the following features from the time-frequency representation: center of energy, bandwidth, instantaneous frequency, duration of the vibration signal, skewness and kurtosis of the frequency distribution, and total energy;
[0138] (2) Autocorrelation: R(τ) = ∫s(t)s(t + τ)dt, where τ is the delay.
[0139] The autocorrelation function measures the similarity between a signal and its time-delayed version. For seismic signals, it can help identify periodic or repetitive patterns;
[0140] (3) Cross-correlation: Rxy(τ) = ∫x(t)y(t + τ)dt
[0141] From historical records and multiple recent seismic signals, cross-correlation can help measure the similarity between these signals;
[0142] Energy change over a long time span: E(t, Tw) = ∫t - Twts2(u)du, the energy of the signal s(t) within the time window Tw. Evaluating the energy change of the signal within a longer time window can provide information about signal persistence and change trends, and this trend information provides data support for modeling the macroscopic movement of the strata;
[0143] (5) Low-frequency energy ratio: The low-frequency component contains potential information caused by crustal deformation and plate movement. The ratio of low-frequency energy to total energy over the entire frequency range can be calculated to obtain this information;
[0144] (6) Statistical characteristics of historical signals: Statistics such as the mean, standard deviation, maximum, and minimum extracted from historical seismic signals. These statistics can help understand the similarities and differences between the current signal and historical patterns;
[0145] 3). For historical earthquake data, label it according to the intensity of the earthquake and the correlation of the foreshock signals, and use a supervised learning method to train a classifier based on the features extracted from shistorical(t) and the generated labels;
[0146] 4). Use the features extracted from scurrent(t) to predict its classification. If the prediction is accurate, keep the features; otherwise, adjust the feature extraction method and repeat steps 1)-3). The process is continuously optimized to improve the accuracy of the prediction.
[0147] Preferably, the visual interpretation of the earthquake prediction model in S6 is achieved by converting known data into a visual gridded geological 3D model, providing more intuitive decision support for earthquake prediction.
[0148] Compared with the prior art, the present invention has the following beneficial effects:
[0149] In this continuous earthquake monitoring and prediction system, by building an inversion model and performing continuous inversion verification with the collected seismic wave signals (the reference physical information is not limited to ground waves, but also includes electromagnetic signals), more detailed geological structure information can be obtained. The earthquake prediction system conducts earthquake rehearsal based on the earthquake inversion model, and compares the correlation between the characteristic signals extracted from the monitoring rehearsal earthquake and the monitored signals, thereby warning of the occurrence of earthquakes, so as to provide people with more accurate warning information and provide support for rescue work after the earthquake. BRIEF DESCRIPTION OF THE DRAWINGS
[0150] Figure 1 It is a flowchart of the overall process of the present invention. DETAILED DESCRIPTION
[0151] In order to make the purpose, technical solution and advantages of the present invention more clear, the technical solution in the embodiment of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiment of the present invention. Obviously, the described embodiment is only a part of the embodiment of the present invention, not all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.
[0152] A continuous earthquake monitoring and prediction system comprises the following steps:
[0153] S1: The data center receives real-time signals transmitted by seismic monitoring beacons and obtains multi-dimensional physical information of stratum changes;
[0154] S2: The data processing center performs noise reduction preprocessing on the information;
[0155] S3: Gridding multiple beacon signals, cross calculation to determine the underground vibration signal source point and the corresponding geological fluid mode;
[0156] S4: The data processing center performs continuous geological inversion of geological structures by acquiring ground waves and other multi-dimensional physical change information;
[0157] S5: Match the series of earthquake precursor signals or the initial seismic wave signals with thresholds to the earthquake omen signals;
[0158] S6: Mark and give an alarm to the signals with a higher degree of matching in the geological inversion system, and then interpret them using the visualization of the earthquake prediction model.
[0159] Preferably, in the S1, the earthquake monitoring beacon collects multi-dimensional physical transformation information through a signal processing algorithm.
[0160] It should be noted that stress release is often the fracture or frictional movement of the rock layer, and obvious changes in electromagnetic signals may occur. The beacon buried deeper can not only monitor the changes in displacement and electromagnetic signals, but also be accompanied by certain air pressure changes (such as the flowing gas in the rock pores). These changes in physical information in other dimensions except the ground wave can be used as evidence for earthquake precursors and occurrences, which can enhance the earthquake recognition ability.
[0161] Preferably, the signal processing algorithm is specifically as follows:
[0162] Add two groups of white noises n a (t) and -n a (t) with equal numerical magnitudes and opposite signs to the original signal x(t), and obtain:
[0163]
[0164] Among them, n i (t) represents the a-th additive Gaussian white noise sequence, and x1(t) and x2(t) represent the noisy signals after adding positive and negative white noises for the a-th time;
[0165] Perform EMD decomposition on x1(t) and x2(t) to obtain the IMF components x1(t) and x2(t):
[0166]
[0167] Among them, M a,b1 and M a,b2 are respectively the b-th IMF component decomposed after adding positive and negative Gaussian white noises for the a-th time, c a,b1 (t) and c a,b2 (t) are the residual components of the EMD decomposition, and B is the number of IMFs;
[0168] Repeat the above steps A times, perform ensemble average operation on the corresponding IMFs above, and the CEEMD decomposition obtains the b-th IMF component c b as:
[0169]
[0170] Among them, M ab is the b-th IMF component decomposed after adding Gaussian white noise for the a-th time;
[0171] The signal s(t) output after CEEMD decomposition is also processed through a non-linear bistable system to output an enhanced signal:
[0172]
[0173]
[0174]
[0175]
[0176]
[0177] r4 = 2l(u(y n + r3) - v(y n + r3) 3 + s n+1 )
[0178] Among them, y(t) is the output enhanced signal, s(t) is the signal output after CEEMD decomposition, u and v are system parameters, is additive Gaussian white noise with a mean of 0 and a variance of σ 2 , y n (t) and y n+1 (t) are respectively the n-th sampling value and the (n + 1)-th sampling value of the output enhanced signal, s n , s n+1 are respectively the n-th sampling value and the (n + 1)-th sampling value of the signal output after CEEMD decomposition, r1, r2, r i and r4 are all substitution values, and l is the integration step size.
[0179] Preferably, stress plates and stress regions are constructed through the results of geological inversion in S4, and a geological model is built. The main process of geological inversion is based on minimizing the difference between the actual seismic data and the geological model, that is, by iteratively updating the geological model to reduce the residual between the observed data and the model prediction.
[0180] Preferably, the process of judging the stress accumulation in the stress region is as follows:
[0181] 1), Data preprocessing and feature extraction:
[0182] The original stress data σ is standardized:
[0183] Among them, μσ and σσ are the mean and standard deviation of the stress data respectively. Extract features such as the peak value, mean value, variance of the stress, etc., and use them as the input of the model;
[0184] 2) Construction of the LSTM model:
[0185] Input layer: Let I(t) be the stress feature data at time t;
[0186] LSTM layer: The following iterative formula is adopted:
[0187] ft = σ(Wf · [Ht-1, I(t)] + bf)
[0188] it = σ(Wi · [Ht-1, I(t)] + bi)
[0189]
[0190]
[0191] ot = σ(Wo · [Ht-1, I(t)] + bo)
[0192] Ht = ot × tanh(Ct)
[0193] Output layer: Let O(t) be the output at time t, representing the predicted stress accumulation state;
[0194] 3) Loss function and optimization:
[0195] Use the training data to optimize the weight and bias parameters of the model. Use the mean squared error as the loss function and adopt an optimization algorithm (such as Adam) to update the parameters:
[0196]
[0197] Among them, Y(i) is the true stress accumulation state at time i;
[0198] Use an optimizer (such as Adam): Initial learning rate: α, estimated first moment: mt, estimated second moment: vt, decay rates: β1, β2.
[0199] Update rule:
[0200] mt = β1mt-1 + (1 - β1)gt
[0201] vt = β2vt-1 + (1 - β2)gt2
[0202]
[0203]
[0204]
[0205] Among them, θ is the parameter to be updated, and gt is the gradient;
[0206] 4), Model verification and evaluation:
[0207] Check the performance of the model through a separate validation dataset (data not involved in training) to prevent overfitting and find the best model and hyperparameters;
[0208] Calculate the evaluation metrics of the model, such as MAE (Mean Absolute Error) or RMSE (Root Mean Squared Error), to quantify the performance of the model on the validation set;
[0209]
[0210] If the performance of the model on the validation dataset is poor, it is necessary to return to the model construction stage for adjustment in combination with the actual situation;
[0211] 5), Result application and integration with geological inversion:
[0212] Use the trained LSTM model to analyze stress data and obtain features related to the predicted stress accumulation value;
[0213] Analyze the predicted stress accumulation state, stress peak, and stress mean value, and conduct analysis in terms of time and space to better understand its distribution and evolution trend. Use the stress-related features obtained from the prediction and analysis to strengthen the input feature layer of the geological inversion algorithm, use the predicted stress value to correct or constrain the update rule of the geological parameter model, or use it as additional information to guide the inversion in the subsequent inversion process;
[0214] The stress stage is a key description of the earthquake process. The occurrence of all types of earthquakes is the result of stress accumulation reaching a certain stage and then being released. Judging which stage the stress area is in, namely stress accumulation, stress criticality, and stress release, helps to judge the time of subsequent earthquake occurrence.
[0215] This LSTM-based method can better capture the time dependence of stress and its relationship with other geological parameters, provide richer information for geological inversion, and thus enhance the inversion effect.
[0216] Based on the above algorithm, iterative optimization is carried out. When the maximum number of iterations is reached, the iteration stops and the inversion result of the geological structure is output. It should be noted that the above is only a basic framework algorithm. In practical applications, each step may involve decomposition into more complex steps and implementation decisions. The inversion result is continuously corrected and needs to be combined with geological information, seismograph data and other information for comprehensive interpretation.
[0217] Preferably, the system can also estimate and judge the plate convergence area, movement trend and stress accumulation stage of the strata based on macroscopic geological conditions, low-frequency micro-seismic data, surface displacement and plate movement history data; among them, the composition of the stress accumulation stage label algorithm is as follows:
[0218] 1), Input data:
[0219] Let Xt represent the input data set at time t, where: σt is the stress data, Et and νt are the elastic modulus and Poisson's ratio of the rock, and dt is the plate displacement estimated based on surface displacement measurement and plate movement history data;
[0220] Xt = [σt, Et, νt, dt]
[0221] 2), GRU model:
[0222] Update gate and reset gate:
[0223] zt = σ(Wz × Xt + Uz × ht-1 + bz)
[0224] rt = σ(Wr × Xt + Ur × ht-1 + br)
[0225] Among them:
[0226] Wz, Wr are the weights from the input data to the update and reset gates, Uz, Ur are the weights from the previous hidden state to the update and reset gates, bz, br are the biases, and σ is the activation function;
[0227] 3), Candidate hidden layer:
[0228]
[0229] Among them: W and U are the weights, b is the bias, and ⊙ represents element-wise multiplication;
[0230] Updated hidden state:
[0231] ht = (1 - zt) ⊙ ht-1 + zt ⊙ ht
[0232] Output layer:
[0233] 1). Degree of stress accumulation: y_stress = σ(Ws×ht + bs), where: Ws is the weight and bs is the bias.
[0234] 2). Stress accumulation stage: y_phase = softmax(Wp×ht + bp), where: Wp is the weight and bp is the bias;
[0235] 4). Loss function:
[0236] Let be the actual label of the stress accumulation degree, be the actual label of the stress accumulation stage, and the loss function L is:
[0237]
[0238] where: MSE is the mean squared error and CrossEntropy is the cross-entropy loss.
[0239] Preferably, the geological inversion process in S4 belongs to a preferred method after stress block delineation, and its specific process is as follows:
[0240] 1) Initial model construction:
[0241] Initialize the geological velocity model according to the existing geological information and vibration data;
[0242] Use the ray tracing algorithm to simulate the propagation of seismic waves in the initial model and determine its main stress area or boundary;
[0243] 2) Inversion iteration:
[0244] Set the discrete grids in space and time:
[0245] xi = i·Δx, yj = j·Δy, zk = k·Δz, tn = n·Δt
[0246] Use the finite difference operator to approximate the derivatives in time and space;
[0247]
[0248] where Γ is the absorption coefficient, representing the energy loss in wave propagation, and s is the source term;
[0249] In each time step, use the above discrete equation to update the pressure wave field;
[0250] 3). Residual calculation:
[0251] Calculate the residual between the observed data and the forward model: R = Dobs - Dcalc, where Dobs is the observed data and Dcalc is the data calculated by the current model;
[0252] 4). Model update:
[0253] Use inverse problem techniques, such as conjugate gradient method or other optimization algorithms, to update the model according to the residuals. The update formula is:
[0254]
[0255] where \(m_{new}\) and \(m_{old}\) are the model parameters after and before the update respectively, \(\alpha\) is the step size, is the residual gradient with respect to the model parameters;
[0256] 5). Convergence check and analysis of geological inversion results:
[0257] Check whether the size of the model update or the reduction of the residuals meets the predefined stopping criteria;
[0258] Analyze the accuracy and reliability of the geological inversion results. Uncertainty quantification and resolution analysis are required. Quantifying uncertainty and resolving resolution are key components of the inverse problem solution;
[0259] Select a prior distribution \(P(m)\) for the model parameter \(m\). This distribution contains the determined formation information and the range of ground wave spectra. Combining the likelihood function \(P(D_{obs}\mid m)\) which describes the probability of the observed data \(D_{obs}\) given the model parameter \(m\), use Bayes' theorem to calculate the posterior distribution:
[0260]
[0261] The posterior distribution \(P(m\mid D_{obs})\) provides the uncertainty information of the updated parameters;
[0262] To ensure the convergence and stability of the algorithm, an appropriate optimization method can be selected, and the step size and convergence criteria can be carefully chosen. The step size \(\alpha\) can be dynamically selected by line search or trust region methods;
[0263] The convergence criteria can be based on the reduction of the residuals or the size of the parameter update. When implementing, first select an appropriate Bayesian framework and optimization algorithm according to the prior information and data errors, and then perform the inversion process by iteratively updating the parameters, evaluating the posterior distribution, and calculating the resolution matrix. In each iteration, it is necessary to check the convergence and stability of the algorithm, and it may be necessary to adjust the optimization parameters and convergence criteria;
[0264] The data processing center in S4 performs inversion of the geological structure based on the multi-dimensional physical change information obtained from the processing results of the signal processing algorithm. In some areas that require detailed investigation, a mobile artificial vibration signal source can be set to detect the formation, and reflection seismic tomography can be carried out, that is, a ground vibration signal with a determined energy is created on the ground. When these ground waves encounter faults or structural blocks, they will produce reflections and diffractions. The distributed beacon grid will receive a series of signals from underground. Under the condition of quantitative movement of the ground signal and multiple signal verifications, the underground formation structure can be detected more accurately.
[0265] Preferably, after the basic formation and stress blocks are determined in S5, the continuous low-frequency vibration signals are manifested as extrusion friction, rock layer fracture, and intermittent rock deformation vibration signals. These are the pre-earthquake omen signals that this invention needs to focus on capturing. For example, obvious formation displacement can often be detected before a strike-slip earthquake occurs.
[0266] The series of earthquake precursor signals in S5 are the extraction markers of the regular characteristics of the vibration signals related to the formation movement of the plot before a major earthquake. In a complete earthquake signal record, except for the displacement, vibration, and strong energy release such as the main shock, the continuous small vibration signals record the development changes and movement trends of its stress area. This time span is relatively long. Using the long attention algorithm of artificial intelligence, the movement of the formation plot is calculated by correlation, and the characteristics of the series of earthquake signals before the earthquake are extracted under similar geological types, and support is provided for subsequent earthquake prediction;
[0267] The formation displacement algorithm based on the ground displacement measurement of the seismic monitoring beacon and the underground inertial measurement data is composed as follows:
[0268] 1) Calculate the velocity v_imu(t) and displacement s_imu(t) from the acceleration a_imu_filtered(t):
[0269] V_imu(t) = ∫₀ᵗ a_imu_filtered(τ)dτ + C1
[0270] S_imu(t) = ∫₀ᵗ v_imu(τ)dτ + C2
[0271] Among them, C1 and C2 are integration constants and can be determined according to the initial values. For multiple IMUs (underground inertial measurement devices), the average value of their displacement data is: s_imu_avg(t) = (1 / n) ∑ᵢ₌₁ⁿ s_imu_i(t), where n is the number of IMUs;
[0272] 2) Displacement comparison and analysis
[0273] Let the surface GPS (or Beidou) displacement be s_gps(t). Compare s_imuavg(t) and s_gps(t) to obtain the displacement difference δ(t): δ(t) = s_imuavg(t) - s_gps(t);
[0274] 3) Displacement judgment of formation stress blocks
[0275] For IMU data s_imui(t) at different depths, where i represents different depths, calculate the difference Δsi(t) between them: Δsi(t) = s_imui(t) - s_imuavg(t)
[0276] If all Δsi(t) are less than a certain threshold, it is determined that the displacement of the formation block from the underground measurement point to the surface area is a uniform displacement; otherwise, there is a non-uniform moving layer from the underground formation to the surface.
[0277] Preferably, during the continuous monitoring of the earthquake by the seismic beacon, pre-seismic events, small earthquakes, or extremely small formation friction displacement signals will continuously occur in the monitoring area. Extract features from these vibration signals before the earthquakes, and at the same time extract some pre-seismic features from similar terrains or records of large earthquakes that have occurred. After extracting the relevant features, monitor these signal features. When the degree of feature matching detected is relatively high, it is predicted that an earthquake will occur. Use algorithms such as deep learning to verify through subsequent facts, including the prediction and verification of aftershocks after the main earthquake, and continuously improve the prediction accuracy;
[0278] Among them, the feature extraction algorithm process is as follows:
[0279] 1) Input the historical shistorical(t) and the seismic signal time series scurrent(t) of the current monitoring area, and output: the feature vector F.
[0280] 2) Extraction steps:
[0281] Use the Short-Time Fourier Transform (STFT) or Wavelet Transform on shistorical(t) and scurrent(t) to obtain the time-frequency representation of the signals;
[0282] Scurrent(f,t) = STFT(scurrent(t))
[0283] Shistorical(f,t) = STFT(shistorical(t))
[0284] The importance of considering the signal correlation over a long time span, especially in the context of earthquake prediction. The autocorrelation and cross-correlation features of signals can help identify regular patterns and potential cyclic behaviors over a long time span:
[0285] (1) Extract the following features from the time-frequency representation: the energy center, bandwidth, instantaneous frequency, duration of the vibration signal, skewness and kurtosis of the frequency distribution, and the total energy;
[0286] (2) Autocorrelation: R(τ) = ∫s(t)s(t + τ)dt, where τ is the delay.
[0287] The autocorrelation function measures the similarity between a signal and its time-delayed version. For earthquake signals, it can help identify periodic or repetitive patterns;
[0288] (3) Cross-correlation: Rxy(τ) = ∫x(t)y(t + τ)dt
[0289] Collect multiple earthquake signals from historical records and recent records. Cross-correlation can help measure the similarity between these signals;
[0290] Energy change over a long time span: E(t, Tw) = ∫t - Twts2(u)du, which is the energy of the signal s(t) within the time window Tw. Evaluating the energy change of the signal within a longer time window can provide information about the signal persistence and change trend. This trend information provides data support for modeling the macroscopic movement of the strata;
[0291] (5) Low-frequency energy ratio: The low-frequency component contains potential information caused by crustal deformation and plate movement. The ratio of low-frequency energy to total energy within the entire frequency range can be calculated to obtain this information;
[0292] (6) Statistical properties of historical signals: Statistics such as the mean, standard deviation, maximum value, and minimum value extracted from historical earthquake signals. These statistics can help understand the similarities and differences between the current signal and historical patterns;
[0293] 3). For historical earthquake data, label it according to the intensity of the earthquake and the correlation of the foreshock signals. Use a supervised learning method to train a classifier based on the features extracted from shistorical(t) and the generated labels;
[0294] 4). Use the features extracted from scurrent(t) to predict its classification. If the prediction is accurate (e.g., successfully predicting small earthquakes or similar historical events), keep the features; otherwise, adjust the feature extraction method and repeat steps 1)-3). The process is continuously optimized to improve the prediction accuracy;
[0295] This earthquake continuous monitoring and prediction system can also monitor long-distance signals. Long-distance signals refer to waveform signals generated by geological movement behaviors that have a long time span, that is, the correlation characteristics of extremely low-frequency signals, which can often reflect the macroscopic movement of the land block. Although the time span is long, it has strong correlation and logic, and can corroborate and explain the evolution of stress blocks.
[0296] The long-distance signal monitoring algorithm is as follows:
[0297] 1. Data preparation:
[0298] 1) Feature calculation: For each time point t, calculate the feature, Ft = FeatureExtractor(s(t))
[0299] where: s(t) is the earthquake signal at time t, and Ft is the feature vector at time t.
[0300] 2) Data division: TrainSet, TestSet = DataSplit(F, α);
[0301] where: F is the entire feature data set, and α is the proportion of the training set.
[0302] 2. Consideration of feature correlation:
[0303] Use feature selection or feature engineering methods to consider the correlation between different features. For example, use the Pearson correlation coefficient or other correlation measurement methods for analysis: Correlation = Pearson(Ft1, Ft2);
[0304] where: Ft1 and Ft2 are feature vectors at different times t1 and t2.
[0305] Based on this correlation analysis, features that have the greatest impact on the model can be selected, and redundant or features that have little relationship with the target can be discarded.
[0306] 3. Model establishment:
[0307] 1) Input layer: I = InputLayer(Ft);
[0308] 2) Long attention layer:
[0309] For a given time point t, the long attention mechanism will consider all previous time points and assign weights to them.
[0310] At = LongAttention(I, Ft-1, Ft-2,..., F0)
[0311] where: At is the attention output at time t after considering all historical data.
[0312] 3) Fully Connected Layer: Ot = FullyConnected(At)
[0313] Where: Ot is the output of the fully connected layer at time t.
[0314] 4) Output Layer: Pt = σ(OutputLayer(Ot))
[0315] Where: σ is the sigmoid function, and Pt is the seismic probability at time t.
[0316] 3. Training the Model:
[0317] Use a loss function L (such as cross - entropy) and an optimizer Optimizer to update the weights of the model.
[0318] L = Loss(P, Y)
[0319] Model = Optimizer(L)
[0320] Where: P is the prediction of the model, and Y is the true label.
[0321] 4. Validation and Testing:
[0322] After the model training is completed, it is necessary to use the test set separated from the entire dataset before to evaluate the performance of the model.
[0323] 1) Model Evaluation:
[0324] Performance = Evaluate(Model, TestSet)
[0325] Where: Model is the trained earthquake prediction model, TestSet is the independent dataset used to evaluate the model, and Performance is the performance of the model on the test set.
[0326] 2) Performance Metrics: The evaluation metrics include:
[0327] The proportion of samples correctly predicted by the model:
[0328] Among the samples correctly predicted as positive, the proportion of those actually being positive:
[0329] Among all samples actually being positive, the proportion correctly predicted as positive
[0330] The harmonic mean of Precision and Recall etc.;
[0331] 3) Result analysis and model tuning:
[0332] For the monitoring and warning of impending earthquakes, manual intervention is required for tuning. While not missing any real earthquake events, the accuracy of the results is ensured. Manual tuning helps with the subsequent evolution and upgrade of the model, and based on the deviation between the forecast results and the actual situation, timely reshaping and optimization are carried out in aspects such as feature engineering, model structure, and hyperparameters.
[0333] 5. Real-time data update:
[0334] Among them, when new earthquake data is received, the model needs to be updated.
[0335] (a) Calculation of new features: Fnew = FeatureExtractor(snew(t))
[0336] (b) Model update: Model = Update(Model, Fnew)
[0337] Where: Fnew are the features of the newly collected data. Within a certain time interval (e.g., daily or weekly), the model is retrained or fine-tuned using the new data.
[0338] 6. Real-time monitoring and warning:
[0339] For the real-time earthquake signal sreal-time(t): Preal-time = Model(sreal-time(t))
[0340] If: Preal-time > threshold, then a warning is issued.
[0341] Preferably, the visual interpretation of the earthquake prediction model in S6 is achieved by converting the known data into a visual grid-based geological 3D model, providing more intuitive decision-making support for earthquake prediction. Visualizing the model is an important way to understand the operation of the geological model system and know the latest progress of model simulation. Using 3D software to construct the geological model can achieve an intuitive visual effect and enhance the interaction between the model system and people. The earthquake prediction model uses intelligent learning algorithms, inputs all the known seismic wave data of the region, conducts manual calibration training on a large amount of earthquake data, updates the data in a timely manner, and promotes the model to tend towards the verification of earthquake results or artificial earthquake signals. Based on the dynamic changes of the geological model obtained by the inversion mechanism, the occurrence of earthquakes is pre-enacted in advance, which helps with the prediction of earthquake occurrence when earthquake precursor information appears;
[0342] Provides visualization of real-time and historical earthquake data, including digital earth models, macroscopic geological data, earthquake source distribution, earthquake magnitude, predicted stress accumulation block diagram, etc. Allows expert users to fine-tune earthquake prediction models, such as modifying parameters, selecting different algorithms, or combining multiple algorithms.
[0343] Copy the latest copy of the system. Without affecting the earthquake monitoring function, users can import historical earthquake data and perform necessary preprocessing through the module, such as model selection, parameter setting and simulation experiments. You can also try to operate multiple prediction models, such as the LSTM model and GRU classifier described above, and allow users to set parameters as needed. The prediction system performs earthquake deduction prediction with accelerated time based on the imported data and the selected model, and presents the prediction results through dynamic charts. Users can select historical time periods, compare actual earthquake events with prediction results, and evaluate the accuracy and robustness of the model. Provide a special toolbox that allows expert users to make in-depth adjustments and optimizations to the model to cope with complex earthquake environments. When the module predicts a possible large-scale earthquake event in a certain area, it prompts personnel to pay attention.
[0344] The above shows and describes the basic principles, main features and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited by the above embodiments. The above embodiments and descriptions are only preferred examples of the present invention and are not intended to limit the present invention. Without departing from the spirit and scope of the present invention, the present invention may have various changes and improvements, which fall within the scope of the present invention. The scope of protection of the present invention is defined by the attached claims and their equivalents.
Claims
1. An earthquake continuous monitoring and prediction system, characterized in that, It includes the following steps: S1: The data center receives the real-time signals transmitted by the seismic monitoring beacons and obtains the multi-dimensional physical information of the formation changes; S2: The data processing center performs noise reduction preprocessing on the information; S3: Grid the multi-beacon signals, and cross-calculate to determine the underground vibration signal source points and the corresponding fluid patterns of the geology; S4: The data processing center performs continuous geological inversion of the geological structure by obtaining ground wave and other multi-dimensional physical change information; S5: Match the series of seismic precursor signals or the initial shock wave signals with thresholds with the earthquake omen signals; S6: Mark and give an alarm to the signals with a higher degree of matching in the geological inversion system, and then interpret them using the visualization of the earthquake prediction model.
2. The earthquake continuous monitoring and prediction system according to claim 1, wherein: In S1, the seismic monitoring beacon collects the multi-dimensional physical transformation information through a signal processing algorithm.
3. The earthquake continuous monitoring and prediction system according to claim 1, characterized in that: The signal processing algorithm is specifically as follows: Add two sets of white noise \(n\) with the same magnitude but opposite signs, following the standard normal distribution, to the original signal \(x(t)\). a (t) and \(-n\) a (t), and obtain: where n i (t) represents the a-th additive white Gaussian noise sequence, and x1(t) and x2(t) represent the noisy signals after adding positive and negative white noises for the a-th time; Perform EMD decomposition on x1(t) and x2(t) to obtain the IMF components x1(t) and x2(t): Among them, M a,b1 and M a,b2 are respectively the b-th IMF component decomposed after adding positive and negative Gaussian white noises for the a-th time, c a,b1 (t) and c a,b2 (t) are the residual components of EMD decomposition, and B is the number of IMFs; Repeat the above steps A times, perform ensemble averaging on the corresponding IMFs above, and the b-th IMF component c is obtained by CEEMD decomposition b which is Among them, M ab is the b-th IMF component obtained by decomposition after adding Gaussian white noise for the a-th time; The signal s(t) output after CEEMD decomposition is also processed by a non-linear bistable system to output an enhanced signal: r4 = 2l(u(y n + r3) - v(y n + r3) 3 + s n+1 ) Among them, y(t) is the output enhanced signal, s(t) is the signal output after CEEMD decomposition, u and v are system parameters, and is additive white Gaussian noise with a mean of 0 and a variance of σ 2 , and y n (t) and y n+1 (t) are the n-th sampling value and the (n + 1)-th sampling value of the output enhanced signal respectively, and s n and s n+1 are the n-th sampling value and the (n + 1)-th sampling value of the signal output after CEEMD decomposition respectively, and r1, r2, r i and r4 are all substitution values, and l is the integration step size.
4. The earthquake continuous monitoring and prediction system according to claim 1, wherein: Construct stress plates and stress regions through the results of geological inversion in S4, build a geological model, and the main process of geological inversion is based on minimizing the difference between the actual seismic data and the geological model, that is, reducing the residuals between the observed data and the model prediction by iteratively updating the geological model.
5. The earthquake continuous monitoring and prediction system according to claim 1, characterized in that: The process of judging the stress accumulation in the stress region is as follows: 1). Data preprocessing and feature extraction: Normalize the original stress data σ: Among them, μσ and σσ are the mean and standard deviation of the stress data respectively, extract features such as the peak value, mean value, variance of the stress, etc., and use them as the input of the model; 2). LSTM model construction: Input layer: Let I(t) be the stress feature data at time t; LSTM layer: Use the following iterative formula: ft=σ(Wf·[Ht-1,I(t)]+bf) it=σ(Wi·[Ht-1,I(t)]+bi) ot=σ(Wo·[Ht-1,I(t)]+bo) Ht=ot×tanh(Ct) Output layer: Let O(t) be the output at time t, representing the predicted stress accumulation state; 3). Loss function and optimization: Use the training data to optimize the weight and bias parameters of the model, use the mean square error as the loss function, and use an optimization algorithm to update the parameters: Among them, Y(i) is the true stress accumulation state at time i; Use the optimizer: Initial learning rate: α, estimated first moment: mt, estimated second moment: vt, decay rates: β1, β2. Update rule: mt=β1mt-1+(1-β1)gt vt=β2vt-1+(1-β2)gt2 Among them, θ is the parameter to be updated, and gt is the gradient; 4). Model verification and evaluation: Check the performance of the model through a separate validation data set, prevent overfitting, and find the best model and hyperparameters; Calculate the evaluation metrics of the model to quantify the performance of the model on the validation set; If the performance of the model on the validation data set is not good, it is necessary to combine the actual situation and go back to the model construction stage for adjustment; 5) Result application and geological inversion integration: Using the trained LSTM model to analyze stress data and obtain the features related to the predicted stress accumulation value; Analyze the predicted stress accumulation state, stress peak value, and stress mean value, conduct analysis in terms of time and space to better understand their distribution and evolution trend. Utilize the stress-related features obtained from prediction and analysis to strengthen the input feature layer of the geological inversion algorithm, use the predicted stress value to correct or constrain the update rule of the geological parameter model, or serve as additional information to guide the inversion in the subsequent inversion process.
6. The earthquake continuous monitoring and prediction system according to claim 5, characterized in that: This system can also estimate and judge the plate convergence area, movement trend, and stress accumulation stage of the formation based on macroscopic geological conditions, low-frequency micro-seismic data, surface displacement, and plate movement history data. Among them, the composition of the stress accumulation stage labeling algorithm is as follows: 1) Input data: Let Xt represent the input data set at time t, where: σt is the stress data, Et and νt are the elastic modulus and Poisson ratio of the rock, and dt is the plate displacement estimated based on surface displacement measurement and plate movement history data; Xt = [σt, Et, νt, dt] 2) GRU model: Update gate and reset gate: zt = σ(Wz × Xt + Uz × ht-1 + bz) rt = σ(Wr × Xt + Ur × ht-1 + br) Where: Wz, Wr are the weights from input data to the update and reset gates, Uz, Ur are the weights from the previous hidden state to the update and reset gates, bz, br are the biases, and σ is the activation function; 3) Candidate hidden layer: Where: W and U are the weights, b is the bias, and ⊙ represents element-wise multiplication; Updated hidden state: ht = (1 - zt) ⊙ ht-1 + zt ⊙ ht Output layer: 1). Degree of stress accumulation: y stress = σ(Ws × ht + bs), where: Ws is the weight, bs is the bias. 2). Stress accumulation stage: y phase = softmax(Wp × ht + bp), where: Wp is the weight, bp is the bias; 4) Loss function: Let be the label of the actual stress accumulation degree, be the actual label in the stress accumulation stage, and the loss function L is as follows: Where: MSE is the mean square error, and CrossEntropy is the cross-entropy loss.
7. The earthquake continuous monitoring and prediction system according to claim 6, characterized in that: The geological inversion process in S4 belongs to an optimized method after stress block delineation, and its specific process is as follows: 1) Initial model construction: According to the existing geological information and vibration data, initialize the geological velocity model; Use the ray tracing algorithm to simulate the propagation of seismic waves in the initial model and determine its main stress area or boundary; 2) Inversion iteration: Set the discrete grids in space and time: xi = i·Δx, yj = j·Δy, zk = k·Δz, tn = n·Δt Use the finite difference operator to approximate the derivatives in time and space; Where, Γ is the absorption coefficient, representing the energy loss in wave propagation, and s is the source term; In each time step, use the above discrete equation to update the pressure wave field; 3). Residual calculation: Calculate the residual between the observed data and the forward model: R = Dobs - Dcalc, where Dobs is the observed data and Dcalc is the data calculated by the current model; 4). Model update: Using inverse problem techniques, the model is updated based on the residuals, and the update formula is: where \(m_{new}\) and \(m_{old}\) are the model parameters after and before the update respectively, and \(\alpha\) is the step size, is the residual gradient with respect to the model parameters; 5). Convergence check and analysis of geological inversion results: Check whether the reduction in the size of the model update or the reduction in residuals meets the predefined stopping criteria; Analyze the accuracy and reliability of the geological inversion results. Uncertainty quantification and resolution analysis are required. Quantifying uncertainty and resolving resolution are key components of the inverse problem solution; Select a prior distribution P(m) for the model parameters m. This distribution contains the determined formation information and the ground wave frequency spectrum range. Combining the likelihood function P(Dobs∣m) describes the probability of the observed data Dobs given the model parameters m. Use Bayes' theorem to calculate the posterior distribution: The posterior distribution P(m∣Dobs) provides information on the uncertainty of the updated parameters; To ensure the convergence and stability of the algorithm, an appropriate optimization method can be selected, and the step size and convergence criteria can be carefully chosen. The step size α can be dynamically selected by line search or trust region methods; The convergence criterion can be based on the reduction in residuals or the size of parameter updates. When implementing, first select an appropriate Bayesian framework and optimization algorithm based on prior information and data errors, and then perform the inversion process by iteratively updating parameters, evaluating the posterior distribution, and calculating the resolution matrix. In each iteration, it is necessary to check the convergence and stability of the algorithm and may require adjusting the optimization parameters and convergence criteria; The data processing center in S4 performs geological structure inversion based on the multi-dimensional physical change information obtained from the processing results of signal processing algorithms. In some areas that require detailed investigation, mobile artificial vibration signal sources can be set up to detect the formation. Reflection seismic tomography, that is, generating a ground vibration signal with a determined energy on the ground. When these ground waves encounter faults or structural blocks, they will produce reflections and diffractions. The distributed beacon grid will receive a series of signals from underground. The underground formation structure can be detected more accurately under the condition of quantitative movement of ground signals and multiple signal verifications.
8. The earthquake continuous monitoring and prediction system according to claim 7, characterized in that: After the basic formation and stress blocks are determined in S5, its continuous low-frequency vibration signals are manifested as extrusion friction, rock layer fracture, and intermittent rock deformation vibration signals. The series of earthquake precursor signals in S5 are the extraction and marking of the regular characteristics of the vibration signals related to the formation movement of the plot before a major earthquake. In a complete earthquake signal record, except for the displacement and strong vibration energy release such as the main shock, its continuous small vibration signals record the development changes and movement trends of its stress area. This time span is relatively long. Applying the long attention algorithm of artificial intelligence, calculate the movement of the formation plot through correlation, and extract the characteristics of the series of earthquake signals before the earthquake under similar geological types, and provide support for subsequent earthquake prediction; The formation displacement algorithm based on the ground displacement measurement of seismic monitoring beacons and underground inertial measurement data is composed as follows: 1) Calculate the velocity v imu(t) and displacement s imu(t) from the acceleration a imufiltered(t): V imu(t)=∫0t a imufiltered(τ)dτ+C1 S imu(t) = ∫₀ᵗ v imu(τ)dτ + C2 Among them, C1 and C2 are integration constants and can be determined according to the initial values. For multiple IMUs, the average value of their displacement data is: s imuavg(t) = (1 / n)∑ᵢ₌₁ⁿ s imui(t), where n is the number of IMUs; 2) Displacement ratio analysis Let the surface GPS (or Beidou) displacement be s gps(t), and the displacement difference δ(t) is obtained by comparing s imuavg(t) and s gps(t): δ(t) = s imuavg(t) - s gps(t); 3) Formation stress block displacement judgment For the IMU data s imui(t) at different depths, where i represents different depths, calculate the difference Δsi(t) between them: Δsi(t) = s imui(t) - s imuavg(t) If all Δsi(t) are less than a certain threshold, it is determined that the displacement of the formation block from the underground measurement point to the surface area is uniform displacement; otherwise, there is a non-uniform moving layer from the underground formation to the surface.
9. The earthquake continuous monitoring and prediction system according to claim 8, wherein: During the continuous monitoring of the earthquake by the said seismic beacon, small earthquakes or extremely small formation friction displacement signals will continuously occur in the monitoring area. Extract features from the vibration signals before these earthquakes, and at the same time extract some pre-earthquake features from similar terrains or recorded large earthquakes that have occurred. After extracting the relevant features, monitor these signal features. When the matching degree of the monitored features is relatively high, it is predicted that an earthquake will occur. Use algorithms such as deep learning and verify with subsequent facts, including the prediction and verification of aftershocks after the main earthquake, and continuously improve the prediction accuracy; Among them, the feature extraction algorithm process is as follows: 1) Input the historical shistorical(t) and the time series of seismic signals scurrent(t) in the current monitoring area, and output: the feature vector F. 2) Extraction steps: Use the Short-Time Fourier Transform (STFT) or Wavelet Transform to obtain the time-frequency representation of the signals for shistorical(t) and scurrent(t); Scurrent(f,t) = STFT(scurrent(t)) Shistorical(f,t) = STFT(shistorical(t)) Considering the importance of the signal correlation over a long time span, especially in the context of earthquake prediction, the autocorrelation and cross-correlation features of the signals can help identify regular patterns and potential cyclic behaviors over a long time span: (1) Extract the following features from the time-frequency representation: energy center, bandwidth, instantaneous frequency, duration of the vibration signal, skewness and kurtosis of the frequency distribution, and total energy; (2) Autocorrelation: R(τ) = ∫s(t)s(t + τ)dt, where τ is the delay. The autocorrelation function measures the similarity between a signal and its time-delayed version. For seismic signals, it can help identify periodic or repetitive patterns; (3) Cross - correlation: $R_{xy}(\tau)=\int x(t)y(t + \tau)dt$ From multiple seismic signals in historical and recent records, cross - correlation can help measure the similarity between these signals; Energy change over a long time span: $E(t,T_w)=\int_{t - T_w}^t s^2(u)du$, the energy of the signal $s(t)$ within the time window $T_w$. Evaluating the energy change of the signal within a longer time window can provide information about the signal's persistence and change trend, and this trend information provides data support for modeling the macroscopic movement of the strata; (5) Low - frequency energy ratio: The low - frequency component contains potential information caused by crustal deformation and plate movement. The ratio of low - frequency energy to total energy over the entire frequency range can be calculated to obtain this information; (6) Statistical characteristics of historical signals: Statistics such as the mean, standard deviation, maximum, and minimum extracted from historical seismic signals. These statistics can help understand the similarities and differences between the current signal and historical patterns; 3). For historical earthquake data, label it according to the intensity of the earthquake and the correlation of the foreshock signals, and use a supervised learning method to train a classifier based on the features extracted from $s_{historical}(t)$ and the generated labels; 4). Use the features extracted from $s_{current}(t)$ to predict its classification. If the prediction is accurate, keep the features; otherwise, adjust the feature extraction method and repeat steps 1) - 3). The process is continuously optimized to improve the accuracy of the prediction.
10. The earthquake continuous monitoring and prediction system according to claim 9, characterized in that: The visual interpretation of the earthquake prediction model in S6 is achieved by converting the known data into a visual grid - like geological 3D model, providing more intuitive decision - making support for earthquake prediction.
Citation Information
Cited By
High-resolution earthquake risk information analysis method and system
CN121634268A