A dynamic early warning method and system for loess liquefaction landslide

CN122531166APending Publication Date: 2026-08-07中国地质环境监测院(自然资源部地质灾害技术指导中心) +1
View PDF 0 Cites 0 Cited by

Patent Information

Authority / Receiving Office
CN · China
Patent Type
Applications(China)
Current Assignee / Owner
中国地质环境监测院(自然资源部地质灾害技术指导中心)
Filing Date
2026-04-29
Publication Date
2026-08-07

AI Technical Summary

Technical Problem

本发明的目的在于提供一种黄土液化滑坡的动态预警方法及系统,用于解决黄土边坡液化滑坡过程中难以及时准确区分土体结构破坏阶段与孔压-结构耦合液化阶段的动态预警问题

✦ Generated by Eureka AI based on patent content.

Smart Images

  • Figure CN122531166A_ABST
    Figure CN122531166A_ABST
Patent Text Reader

Abstract

The present application relates to the technical field of landslide analysis and early warning, in particular to a dynamic early warning method and system for loess liquefaction landslide, which obtains the resistivity time series and pore water pressure time series at different depths of the loess slope; separates the mutation component from the resistivity time series to calculate the resistivity drop rate at the current time; calculates the excess pore water pressure ratio at the current time according to the pore water pressure time series; forms a two-dimensional observation vector and a two-dimensional observation sequence by combining the resistivity drop rate and the excess pore water pressure ratio; calculates the identification state at the current time by using the hidden Markov model on the two-dimensional observation sequence; performs multi-level early warning according to the type of the identification state at the current time; and realizes accurate grading early warning for loess liquefaction landslide.
Need to check novelty before this filing date? Find Prior Art

Description

Technical Field

[0001] This invention relates to the field of landslide analysis and early warning technology, specifically a dynamic early warning method and system for loess liquefaction landslides. Background Technology

[0002] Loess is a unique type of soil widely distributed in Northwest China, characterized by its large porosity, low compaction, and sensitivity to water. Under external loads such as heavy rainfall or earthquakes, loess slopes are highly susceptible to landslides, primarily caused by liquefaction. Loess liquefaction landslides are characterized by rapid initiation, long distance traveled, and high destructive energy.

[0003] However, the instability process of loess liquefaction landslides is not dominated by a single mechanism, but rather undergoes a multi-stage evolution, from the gradual destruction of the soil skeleton structure to the rapid accumulation of pore water pressure, and then the coupling of these two factors driving the overall instability of the slope. The physical and mechanical mechanisms at different stages are fundamentally different: in the stage dominated by structural failure, the cementation between soil particles breaks down, and the bearing capacity of the soil skeleton decreases, but the pore pressure response is not yet significant; while in the pore pressure-structure coupled liquefaction stage, the rapid accumulation of excess pore water pressure and the continuous deterioration of the soil structure mutually reinforce each other, and the slope is already in a critical instability state, with a much higher degree of danger than in the former stage. Therefore, establishing an effective dynamic early warning system for loess liquefaction landslides is of significant practical importance for disaster prevention and mitigation. Summary of the Invention

[0004] (1) Technical problems to be solved The purpose of this invention is to provide a dynamic early warning method and system for loess liquefaction landslides, which can solve the problem of dynamic early warning in which it is difficult to distinguish between the soil structure failure stage and the pore pressure-structure coupled liquefaction stage in a timely and accurate manner during the process of loess slope liquefaction landslides.

[0005] (2) Technical solution To achieve the above objectives, on the one hand, the present invention provides a dynamic early warning method for loess liquefaction landslides, the method comprising: S1. Obtain the resistivity time series and pore water pressure time series at different depths of the loess slope.

[0006] S2. Perform signal decomposition on the resistivity time series to separate the baseline drift component, which characterizes the slow electrical changes caused by loess collapse, and the abrupt change component, which characterizes the sudden failure of the soil structure, and remove the baseline drift component; calculate the resistivity drop rate at the current moment based on the abrupt change component, where the resistivity drop rate is the relative magnitude of the resistivity decrease per unit time.

[0007] S3. Calculate the excess pore water pressure ratio at the current moment based on the pore water pressure time series. The excess pore water pressure ratio is the ratio of excess pore water pressure to the initial effective stress.

[0008] S4. The resistivity drop ratio and the ratio of ultrastatic pore water pressure are used to form a two-dimensional observation vector, and a two-dimensional observation sequence is formed by sliding sampling.

[0009] S5. Calculate the filtering probability of each hidden state at the current time using a Hidden Markov Model for the two-dimensional observation sequence, and use the hidden state corresponding to the maximum filtering probability as the identification state at the current time.

[0010] S6. If the current identification state is a structural damage-dominated state, then issue a first-level warning signal; if the current identification state is a pore pressure-structure coupling liquefaction state, then issue a second-level warning signal; otherwise, do not issue a warning signal.

[0011] Furthermore, the method for performing signal decomposition on the resistivity time series to separate the baseline drift component characterizing the slow electrical changes caused by loess subsidence and the abrupt change component characterizing the sudden failure of soil structure includes: The resistivity time series is decomposed to obtain a series of eigenmode function components arranged from high frequency to low frequency and a residual term.

[0012] Calculate the sample entropy or zero-crossing rate of each intrinsic mode function component. Identify intrinsic mode function components whose sample entropy is lower than the first entropy threshold or whose zero-crossing rate is lower than the first zero-crossing rate threshold as low-frequency components. Then, superimpose all low-frequency components with the residual term to obtain the baseline drift component that characterizes the slow electrical changes caused by loess subsidence.

[0013] The eigenmode function components whose remaining sample entropy is higher than or equal to the first entropy threshold and whose zero-crossing rate is higher than or equal to the first zero-crossing rate threshold are reconstructed to obtain the mutation component characterizing the sudden failure of the soil structure.

[0014] Furthermore, the method for decomposing the resistivity time series to obtain a series of eigenmode function components arranged from high frequency to low frequency and a residual term includes: Multiple sets of paired positive and negative adaptive white noise are added to the resistivity time series. Each set of paired positive and negative adaptive white noise constitutes a noise copy. Empirical mode decomposition is performed on each noise copy to obtain the first-order intrinsic mode function components corresponding to each noise copy.

[0015] The first-order intrinsic mode function components of all noise replicas are ensembled and averaged to obtain the first-order intrinsic mode function components of the resistivity time series. The first-order intrinsic mode function components are then subtracted from the resistivity time series to obtain the first-order residual components.

[0016] Add positive and negative paired adaptive white noise to the first-order residual component, and repeat the empirical mode decomposition and ensemble averaging operations to extract the second-order intrinsic mode function component, the third-order intrinsic mode function component, and so on up to the Nth-order intrinsic mode function component.

[0017] The residual components that cannot be further decomposed are taken as residual terms, and the extracted first-order intrinsic mode function components to the Nth-order intrinsic mode function components are arranged in descending order of frequency to obtain a series of intrinsic mode function components arranged from high frequency to low frequency.

[0018] Furthermore, the method for determining N includes: Perform a monotonicity test on the residual component of the current order obtained after each extraction. If the residual component of the current order is a monotonic function in its domain, stop decomposition and subtract one from the current order as the value of N.

[0019] If the current-order residual component is not a monotonic function, then calculate the permutation entropy of the current-order residual component, whereby the permutation entropy characterizes the complexity of the residual component.

[0020] Based on the expected energy distribution of resistivity mutation components in loess slopes during earthquakes or rainfall events that induce liquefaction landslides, a signal entropy threshold is pre-set. If the permutation entropy of the current-order residual component is lower than the signal entropy threshold, it is determined that the current-order residual component no longer contains identifiable structural damage mutation information, and the decomposition is stopped. The current order is reduced by one to obtain the value of N.

[0021] If the permutation entropy of the current order residual component is higher than or equal to the signal entropy threshold, then continue to execute the next order decomposition until the monotonicity or low entropy stopping condition is met.

[0022] Furthermore, the method of constructing a two-dimensional observation vector by combining the resistivity drop ratio with the excess pore water pressure ratio, and forming a two-dimensional observation sequence through sliding sampling, includes: Set a time window length and a sliding step size; at each sampling time, combine the calculated resistivity drop rate at the current sampling time with the excess pore water pressure ratio at the current sampling time into a two-dimensional observation vector and store it.

[0023] Each sampling time determined by the sliding step size is taken as the end time of the current window. Starting from the end time of the current window, the number of sampling points included in the time window length is traced back, and a two-dimensional observation vector sequence consisting of multiple consecutive two-dimensional observation vectors is extracted from the stored two-dimensional observation vectors.

[0024] Furthermore, the method for setting a time window length and a sliding step size includes: Based on the typical evolution time of loess slopes from initial soil structure failure to liquefaction instability under rainfall or seismic loads, an initial window length covering the typical evolution time is determined as the time window length.

[0025] Based on the data sampling period of resistivity drop rate and excess pore water pressure ratio, as well as the response speed requirements of the landslide early warning system, a sliding step size smaller than the time window length is set so that there is partial overlap between two adjacent sliding intercepted two-dimensional observation sequences.

[0026] The time window length and sliding step size are converted into window length points and sliding step size points, respectively, in terms of the number of sampling points. The window length points are the integer obtained by dividing the time window length by the sampling period, and the sliding step size points are the integer obtained by dividing the sliding step size by the sampling period.

[0027] Furthermore, the method of calculating the filtering probability of each hidden state at the current time using a Hidden Markov Model on the two-dimensional observation sequence, and using the hidden state corresponding to the maximum filtering probability as the identification state at the current time, includes: The three hidden states of the Hidden Markov Model are denoted as the steady state, the structural failure-dominated state, and the pore pressure-structure coupling liquefaction state, respectively. The filtering probabilities of each hidden state are initialized when the current time is zero. The filtering probability of the steady state is set to 1, and the filtering probabilities of the structural failure-dominated state and the pore pressure-structure coupling liquefaction state are both set to 0.

[0028] For the current time t, obtain the filtered probability of each hidden state at the previous time t-1, multiply the filtered probability of each hidden state at the previous time t-1 by the state transition probability of each hidden state at the current time t, and sum the products to obtain the predicted probability of each hidden state at the current time t.

[0029] Substitute the two-dimensional observation vector at the current time t into the pre-trained Gaussian mixture model to calculate the observation likelihood probability of observing the two-dimensional observation vector in each hidden state at the current time t.

[0030] Multiply the predicted probability of each hidden state at the current time t by the corresponding observation likelihood probability to obtain the unnormalized filtered probability of each hidden state at the current time t.

[0031] Calculate the sum of the unnormalized filtering probabilities of all hidden states at the current time t, and normalize the unnormalized filtering probability of each hidden state using the sum of the unnormalized filtering probabilities to obtain the filtering probability of each hidden state at the current time t.

[0032] Select the filter probability with the largest value from the filter probabilities of the three hidden states at the current time t, and take the hidden state corresponding to the maximum filter probability as the recognition state at the current time t.

[0033] Furthermore, the method of substituting the two-dimensional observation vector at the current time t into a pre-trained Gaussian mixture model to calculate the observation likelihood probability of the observed two-dimensional observation vector in each hidden state at the current time t includes: For each of the three hidden states, an independent Gaussian mixture model is configured. The Gaussian mixture model contains multiple Gaussian components, and each Gaussian component has a preset mean vector, covariance matrix, and mixing coefficients.

[0034] Substitute the two-dimensional observation vector at the current time t into each Gaussian component of the Gaussian mixture model of the current hidden state, and calculate the single Gaussian probability density value of the observed two-dimensional observation vector under each Gaussian component.

[0035] The weighted probability density value of each Gaussian component is obtained by multiplying the single Gaussian probability density value of each Gaussian component by the corresponding mixing coefficient.

[0036] The observation likelihood probability of the two-dimensional observation vector at time t is obtained by summing the weighted probability density values ​​of all Gaussian components in the current hidden state.

[0037] On the other hand, based on the same inventive concept, this invention also provides a dynamic early warning system for loess liquefaction landslides, the system comprising: The data sequence acquisition module is used to acquire resistivity time series and pore water pressure time series at different depths of loess slopes.

[0038] The resistivity signal processing module is used to decompose the resistivity time series, separate the baseline drift component, which characterizes the slow electrical changes caused by loess collapse, and the abrupt change component, which characterizes the sudden failure of the soil structure, and remove the baseline drift component; calculate the resistivity descent rate at the current moment based on the abrupt change component, where the resistivity descent rate is the relative magnitude of the resistivity decrease per unit time.

[0039] The module for calculating the excess pore water pressure ratio is used to calculate the excess pore water pressure ratio at the current moment based on the pore water pressure time series. The excess pore water pressure ratio is the ratio of the excess pore water pressure to the initial effective stress. A two-dimensional observation sequence construction module is used to construct a two-dimensional observation vector by combining the resistivity drop rate with the ratio of excess pore water pressure, and to form a two-dimensional observation sequence through sliding sampling.

[0040] The hidden state identification module is used to calculate the filtering probability of each hidden state at the current time using a hidden Markov model on the two-dimensional observation sequence, and to use the hidden state corresponding to the maximum filtering probability as the identification state at the current time.

[0041] The warning signal issuing module is used to issue a first-level warning signal if the current identification state is a structural damage-dominated state; to issue a second-level warning signal if the current identification state is a pore pressure-structure coupling liquefaction state; otherwise, no warning signal is issued.

[0042] (3) Beneficial effects Compared with the prior art, the beneficial effects of the present invention are: 1. By effectively separating the background baseline drift caused by loess subsidence and the resistivity abrupt change component caused by sudden damage to the soil structure, the resistivity drop rate feature can be accurately extracted. At the same time, by combining the excess pore water pressure ratio to construct a two-dimensional observation vector integrating multiple physics fields, the extraction accuracy and reliability of key physical features in the evolution of liquefaction landslides are significantly improved.

[0043] 2. A hidden Markov model is used to identify the dynamic state of a two-dimensional observation sequence. By utilizing its probabilistic modeling capability for non-stationary time series signals, the model can distinguish three evolution stages in real time under noise interference and signal uncertainty: stable state, structural failure-dominated state, and pore pressure-structure coupled liquefaction state. Based on this, corresponding graded early warning signals are output, providing a more reliable technical means for disaster prevention and mitigation of loess liquefaction landslides. Attached Figure Description

[0044] Figure 1 This is a flowchart of a dynamic early warning method for loess liquefaction landslides according to the present invention.

[0045] Figure 2 This is a schematic diagram of the module composition of a dynamic early warning system for loess liquefaction landslides according to the present invention. Detailed Implementation

[0046] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0047] Example 1: As Figure 1 As shown in the figure, this embodiment provides a dynamic early warning method for loess liquefaction landslides, the method comprising: S1. Obtain the resistivity time series and pore water pressure time series at different depths of the loess slope.

[0048] For example, in actual engineering deployments, resistivity sensor arrays and pore water pressure sensors are buried at different depths along the potential sliding surface of the loess slope. The resistivity and pore water pressure readings of each sensor are continuously recorded at fixed sampling intervals, thereby obtaining continuous resistivity and pore water pressure time series over time. The burial depth of the resistivity sensors should cover the loess layer above the potential sliding surface of the slope, typically including three depth ranges: shallow (1 to 2 meters below the surface), middle (3 to 5 meters), and deep (7 to 10 meters), to comprehensively capture early warning signals of soil structural failure at different depths. The pore water pressure sensors should be deployed correspondingly to the resistivity sensors on the same depth profile to ensure spatial comparability of the two physical quantities. The sampling interval should be reasonably set according to the response speed requirements of the early warning system. In heavy rainfall or earthquake events, the loess structural failure process may change significantly within seconds to tens of seconds; therefore, it is recommended that the sampling interval not exceed 10 seconds to ensure effective capture of sudden signals. Multi-point deployment at different depths is necessary because pore pressure accumulation and structural damage during loess liquefaction often occur preferentially in weak layers at specific depths. Multi-depth monitoring helps to capture the earliest signs of instability and provides a more complete spatial information basis for subsequent signal analysis.

[0049] S2. The method for performing signal decomposition on the resistivity time series to separate the baseline drift component characterizing the slow electrical changes caused by loess collapse and the abrupt change component characterizing the sudden failure of soil structure includes: Methods for decomposing the resistivity time series to obtain a series of eigenmode function components arranged from high frequency to low frequency and a residual term include: Multiple sets of paired positive and negative adaptive white noise are added to the resistivity time series. Each set of paired positive and negative adaptive white noise constitutes a noise copy. Empirical mode decomposition is performed on each noise copy to obtain the first-order intrinsic mode function components corresponding to each noise copy.

[0050] The first-order intrinsic mode function components of all noise replicas are ensembled and averaged to obtain the first-order intrinsic mode function components of the resistivity time series. The first-order intrinsic mode function components are then subtracted from the resistivity time series to obtain the first-order residual components.

[0051] Add positive and negative paired adaptive white noise to the first-order residual component, and repeat the empirical mode decomposition and ensemble averaging operations to extract the second-order intrinsic mode function component, the third-order intrinsic mode function component, and so on up to the Nth-order intrinsic mode function component.

[0052] The residual components that cannot be further decomposed are taken as residual terms, and the extracted first-order intrinsic mode function components to the Nth-order intrinsic mode function components are arranged in descending order of frequency to obtain a series of intrinsic mode function components arranged from high frequency to low frequency.

[0053] For example, during rainfall or earthquake-induced processes, the resistivity changes in loess simultaneously contain two distinct components: one is a long-period, slow-drifting background change caused by loess subsidence and a gradual increase in water content, known as the baseline drift component; the other is a rapid decrease in resistivity within a short period, caused by the sudden rupture of cementation between soil particles and the rapid collapse of the skeletal structure, known as the abrupt change component. If these two components are not separated, the baseline drift will mask structural damage characteristics, severely reducing the accuracy of early warning systems in identifying early structural damage signals.

[0054] Adaptive white noise with paired positive and negative noise is added multiple times to the resistivity time series. Each time a pair of positive and negative white noise is added, it forms a noise replica. The positive and negative white noise have equal amplitudes but opposite signs. Their amplitudes are adaptively determined according to a fixed proportion to the standard deviation of the current signal to be decomposed. This proportion is usually set to 0.2, meaning the standard deviation of the white noise is 0.2 times the standard deviation of the current signal to be decomposed, ensuring that the noise is sufficient to assist in decomposition but does not excessively interfere with the signal itself. It is recommended that there be no fewer than 100 noise replicas to ensure that the statistical residual of the noise after ensemble averaging is sufficiently small. Standard empirical mode decomposition is performed independently on each noise replica to extract the first-order intrinsic mode function (IMF) components of each replica. Then, the IMF components of all replicas are ensemble-averaged to eliminate the interference of artificially added noise, obtaining the first-order IMF components of the resistivity time series. The first-order IMF components are then subtracted from the resistivity time series to obtain the first-order residual components. Then, repeat the above operation on the first-order residual component, and extract the second-order, third-order and up to the Nth-order intrinsic mode function components in sequence, arrange them in descending order of frequency, and finally retain the residual components that cannot be further decomposed as residual terms.

[0055] The method for determining N includes: Perform a monotonicity test on the residual component of the current order obtained after each extraction. If the residual component of the current order is a monotonic function in its domain, stop decomposition and subtract one from the current order as the value of N.

[0056] If the current-order residual component is not a monotonic function, then calculate the permutation entropy of the current-order residual component, whereby the permutation entropy characterizes the complexity of the residual component. Based on the expected energy distribution of resistivity mutation components in loess slopes during earthquakes or rainfall events that induce liquefaction landslides, a signal entropy threshold is pre-set. If the permutation entropy of the current-order residual component is lower than the signal entropy threshold, it is determined that the current-order residual component no longer contains identifiable structural damage mutation information, and the decomposition is stopped. The current order is reduced by one to obtain the value of N.

[0057] If the permutation entropy of the current order residual component is higher than or equal to the signal entropy threshold, then continue to execute the next order decomposition until the monotonicity or low entropy stopping condition is met.

[0058] For example, after each extraction, the residual component of the current order is judged to constitute a monotonic function within its complete time domain. If all first-order differences of the residual component are non-negative or all are non-positive, that is, the residual component is monotonically non-decreasing or monotonically non-increasing throughout the time domain, then the residual component is considered to have degenerated into a trend term and no longer contains any oscillatory components. At this time, the decomposition is stopped, and N is taken as the current order minus one. The second stopping criterion is the permutation entropy test: if the residual component does not satisfy monotonicity, the permutation entropy of the residual component is calculated. The permutation entropy is calculated as follows: set the embedding dimension m (usually an integer between 3 and 5) and the time delay τ (usually 1), divide the residual component into continuous subsequences of length m according to the time delay τ, sort the elements in each subsequence in ascending order of value, record the symbol pattern corresponding to the permutation, count the relative frequency of each symbol pattern in all subsequences, denoted as the probability distribution, and calculate the permutation entropy according to the Shannon entropy formula: ;in, For the first The probability of occurrence of a certain symbol arrangement pattern For the embedding dimension is The permutation entropy is the total number of all possible permutation patterns. A higher permutation entropy value indicates that the residual components are more complex and contain more effective information; a lower permutation entropy value indicates that the residual components are closer to monotonic or periodic changes, and the identifiable structural damage information has been basically exhausted.

[0059] The signal entropy threshold is determined as follows: Based on existing monitoring data or indoor simulation test data of loess slopes in typical earthquake or heavy rainfall-induced liquefaction landslide events, the resistivity time series known to be in stable and structural failure states are decomposed respectively. The permutation entropy range corresponding to when it is confirmed that there is no longer any structural failure mutation information in each residual component is statistically analyzed. The upper limit of this range is taken as a reference for the signal entropy threshold, and then corrected by combining data from multiple historical events to finally determine the signal entropy threshold.

[0060] Calculate the sample entropy or zero-crossing rate of each intrinsic mode function component. Identify intrinsic mode function components whose sample entropy is lower than the first entropy threshold or whose zero-crossing rate is lower than the first zero-crossing rate threshold as low-frequency components. Then, superimpose all low-frequency components with the residual term to obtain the baseline drift component that characterizes the slow electrical changes caused by loess collapsing. The eigenmode function components whose remaining sample entropy is higher than or equal to the first entropy threshold and whose zero-crossing rate is higher than or equal to the first zero-crossing rate threshold are reconstructed to obtain the mutation component characterizing the sudden failure of the soil structure.

[0061] For example, after decomposition, it is necessary to classify the intrinsic mode function components of each order to determine which belong to baseline drift and which belong to abrupt changes. For each intrinsic mode function component, two statistics are calculated: sample entropy and zero-crossing rate. The sample entropy is calculated as follows: Given a template length *m* and a similarity tolerance *r* (usually *r* is 0.1 to 0.2 times the component standard deviation), the sample entropy is calculated by taking the negative logarithm of the ratio between the number of similarity logarithms between a template vector of length *m* and other position vectors within the tolerance *r*, and the number of similarity logarithms between template vectors of length *m+1*. Lower sample entropy indicates higher self-similarity and lower complexity of the component, corresponding to low-frequency, slowly varying baseline drift characteristics. The zero-crossing rate is calculated by counting the number of times the component sequence crosses from a positive value to a negative value or from a negative value to a positive value per unit time (or unit sampling points). A lower zero-crossing rate indicates slower component change, corresponding to low-frequency components.

[0062] When the slope is in a stable state (before being disturbed by rainfall or earthquakes), a representative reference resistivity time series is collected and fully adaptive noise ensemble empirical mode decomposition is performed. The sample entropy and zero-crossing rate of each intrinsic mode function component obtained by decomposition are calculated one by one. The maximum sample entropy and the maximum zero-crossing rate of each component that can be clearly attributed to the slow collapsibility background change (i.e., the frequency is significantly lower than the frequency of the soil structure failure signal) in this stage are used as the initial reference values ​​of the first entropy threshold and the first zero-crossing rate threshold, respectively. The initial reference values ​​are corrected by combining the calibration results of multiple historical monitoring events to determine the values ​​of the first entropy threshold and the first zero-crossing rate threshold.

[0063] The baseline drift component is removed; the resistivity descent rate at the current moment is calculated based on the abrupt change component, where the resistivity descent rate is the relative magnitude of the resistivity decrease per unit time. For example, the resistivity descent rate is defined as the relative magnitude of the resistivity decrease of the abrupt component per unit time, reflecting the rate characteristic of the rapid decrease in resistivity when the soil skeleton structure undergoes sudden failure. The calculation formula is as follows: ;in, The resistivity value of the abrupt change component at the current moment (unit: ), The resistivity value of the abrupt change component at the previous sampling time (unit: ), Sampling time interval (unit: seconds); The dimensions are When the soil skeleton structure experiences sudden failure, the resistivity of the abrupt component decreases rapidly within a short period, and the resistivity descent rate will show a significant positive value, which can be used as a characteristic quantity of the structural failure stage. If the resistivity of the abrupt component does not decrease, the resistivity descent rate is close to zero or negative, indicating that the structure has not yet experienced significant sudden failure at this moment. It should be noted that the calculation of the resistivity descent rate is based on the abrupt component after removing the baseline drift component.

[0064] S3. Calculate the excess pore water pressure ratio at the current moment based on the pore water pressure time series. The excess pore water pressure ratio is the ratio of excess pore water pressure to the initial effective stress. For example, the excess pore water pressure ratio is defined as the ratio of the current excess pore water pressure to the initial effective stress. It is used to quantify the proportion of the accumulation of pore water pressure in the soil relative to the initial effective stress of the soil. The calculation formula is as follows: ;in, The current static pore water pressure (unit: kPa) is equal to the pore water pressure value measured by the pore water pressure sensor at the current moment minus the hydrostatic pressure at the sensor's burial depth. The hydrostatic pressure is determined by the product of the groundwater level at the sensor's burial depth and the specific weight of water. The initial effective stress (unit: kPa) at the sensor burial depth is calculated by subtracting the hydrostatic pressure from the self-weight stress of the overlying soil. It is a known quantity determined by survey data before the sensor installation is completed and monitoring officially begins, and is used as a fixed parameter during the monitoring process. It is a dimensionless ratio.

[0065] The physical significance of the excess pore water pressure ratio lies in the following: when the excess pore water pressure ratio approaches 0, the pore pressure state of the soil is close to the hydrostatic pressure level, the effective stress is basically undisturbed, and the soil is in a stable state; when the excess pore water pressure ratio continues to increase and approaches 1, the pore water bears almost all the overburden pressure, the effective stress approaches zero, and the soil is about to completely lose its shear strength, which is a classic quantitative indicator of the critical liquefaction state. The excess pore water pressure ratio and the resistivity descent rate together constitute characteristic quantities that characterize two essentially different physical mechanisms in the evolution of loess liquefaction landslides. The joint observation of the two can significantly improve the accuracy of state identification.

[0066] S4. The method of constructing a two-dimensional observation vector by combining the resistivity drop ratio with the excess pore water pressure ratio, and forming a two-dimensional observation sequence through sliding sampling, includes: Methods for setting a time window length and a sliding step include: Based on the typical evolution time of loess slopes from initial soil structure failure to liquefaction instability under rainfall or seismic loads, an initial window length covering the typical evolution time is determined as the time window length.

[0067] Based on the data sampling period of resistivity drop rate and excess pore water pressure ratio, as well as the response speed requirements of the landslide early warning system, a sliding step size smaller than the time window length is set so that there is partial overlap between two adjacent sliding intercepted two-dimensional observation sequences.

[0068] The time window length and sliding step size are converted into window length points and sliding step size points, respectively, in terms of the number of sampling points. The window length points are the integer obtained by dividing the time window length by the sampling period, and the sliding step size points are the integer obtained by dividing the sliding step size by the sampling period.

[0069] For example, the time window length is determined based on the typical evolution time of a loess slope from initial soil structural failure to liquefaction instability under rainfall or seismic loads. Based on existing indoor tests and field case analyses of loess landslides, this evolution process typically lasts from several minutes to tens of minutes in heavy rainfall-induced scenarios and from tens of seconds to several minutes in earthquake-induced scenarios. The initial window length covering the typical evolution time can be determined by combining historical test data or numerical simulation results from specific monitoring projects. Taking a heavy rainfall-induced scenario as an example, if the typical evolution time is 20 minutes, the time window length can be set to 20 minutes. The sliding step size must be smaller than the time window length to ensure partial overlap between adjacent observation sequences, enabling the early warning system to continuously track state changes with high temporal resolution. An excessively large sliding step size will result in excessively long time intervals between adjacent state identifications, potentially missing critical nodes in rapid evolution; an excessively small sliding step size will increase the computational burden. Considering both response speed and computational efficiency, the sliding step size is usually set to an integer multiple of the sampling period, with its corresponding time length not exceeding one-tenth of the time window length. Taking a sampling period of 10 seconds and a time window length of 20 minutes as an example, the sliding step can be set to 10 seconds, meaning that a state recognition is triggered every time a new sampling moment is reached. Divide the time window length and sliding step by the sampling period to convert them into the number of sampling points, which are the window length points and the sliding step points, respectively. For example, when the time window length is 20 minutes and the sampling period is 10 seconds, the window length points are 120 sampling points.

[0070] At each sampling time, the calculated resistivity drop rate at the current sampling time is combined with the excess pore water pressure ratio at the current sampling time into a two-dimensional observation vector and stored. Each sampling time determined by the sliding step size is taken as the end time of the current window. Starting from the end time of the current window, the number of sampling points included in the time window length is traced back, and a two-dimensional observation vector sequence consisting of multiple consecutive two-dimensional observation vectors is extracted from the stored two-dimensional observation vectors.

[0071] For example, at each sampling time, the resistivity descent rate calculated at the current time is used. Compared with the pressure of superstatic pore water Before combining them in a fixed order, the two features must be standardized separately. The standardization method uses Z-score standardization, calculated using the offline training dataset as a benchmark. and Mean and standard deviation over all training samples: ;in, , On the training set respectively The mean and standard deviation; , On the training set respectively The mean and standard deviation of the features. After standardization, both features are dimensionless, with a mean of 0 and a standard deviation of 1 on the training set, making them comparable in magnitude. Standardization parameters Once determined during the offline training phase, the data is stored permanently. During the real-time monitoring phase, the real-time observations are standardized using the same parameters and are not updated further. After standardization, the two standardized features are combined in a fixed order to form a two-dimensional observation vector. The resistivity drop rate component reflects the degree of sudden damage to the soil skeleton structure, while the excess pore water pressure ratio component reflects the degree of pore water pressure accumulation. The combination of the two allows the observation vector to simultaneously capture two fundamentally different physical mechanisms in the evolution of loess liquefaction landslides: structural damage and pore pressure accumulation.

[0072] In actual operation, whenever new sampling data arrives and the cumulative number of arrived points meets the number of points corresponding to the sliding step size, the current sampling time is taken as the end time of the current window. The window length is then traced back by points, and a sequence of two-dimensional observation vectors consisting of multiple consecutive two-dimensional observation vectors is extracted from the stored historical records of two-dimensional observation vectors. This sequence is then used as the complete observation sequence input into the Hidden Markov Model at the current time. In the initial stage of monitoring, if the number of stored historical two-dimensional observation vectors is insufficient to meet the window length, the current observation sequence is constructed using all stored historical two-dimensional observation vectors.

[0073] S5. The method of calculating the filtering probability of each hidden state at the current time using a Hidden Markov Model for the two-dimensional observation sequence, and using the hidden state corresponding to the maximum filtering probability as the identification state at the current time includes: The three hidden states of the Hidden Markov Model are denoted as the steady state, the structural failure-dominated state, and the pore pressure-structure coupling liquefaction state, respectively. The filtering probabilities of each hidden state are initialized when the current time is zero. The filtering probability of the steady state is set to 1, and the filtering probabilities of the structural failure-dominated state and the pore pressure-structure coupling liquefaction state are both set to 0.

[0074] For the current time t, obtain the filtered probability of each hidden state at the previous time t-1, multiply the filtered probability of each hidden state at the previous time t-1 by the state transition probability of each hidden state at the current time t, and sum the products to obtain the predicted probability of each hidden state at the current time t.

[0075] For example, the stable state corresponds to the stage where the slope is generally safe, the soil skeleton is intact, and no abnormal accumulation of pore pressure is observed; the structural failure-dominated state corresponds to the stage where the cementation between soil particles begins to break, the skeleton bearing capacity decreases, the resistivity drop rate increases significantly, but the excess pore water pressure ratio has not yet accumulated in large quantities; the pore pressure-structure coupled liquefaction state corresponds to the stage where the excess pore water pressure ratio rises rapidly, mutually promotes the continuous deterioration of the skeleton structure, and the slope is already at the critical point of instability.

[0076] The parameters of the Hidden Markov Model (HMM) consist of three parts: the initial state probability vector, the state transition probability matrix, and the observation emission probability model. At time zero, the initial state probability vector sets the filtered probability of the steady state to 1, and the filtered probabilities of the structural failure-dominated state and the pore pressure-structure coupled liquefaction state to 0. This reflects the prior assumption that the slope is in a stable state when monitoring begins, which aligns with the practical situation in engineering where monitoring systems are typically deployed before any obvious anomalies appear on the slope.

[0077] The state transition probability matrix describes the probability of transitions between three hidden states at adjacent time points. The element in the i-th row and j-th column of the matrix represents the probability of transitioning from the i-th hidden state to the j-th hidden state at the next time point, and the sum of the elements in each row equals 1. The state transition probability matrix is ​​obtained by training a hidden Markov model on multiple sets of indoor triaxial liquefaction test data of loess and field landslide monitoring records, and iteratively optimizing it using the Baum-Welch expectation-maximization algorithm, reflecting the statistical transition laws between the three stages in the loess liquefaction evolution process.

[0078] Methods for calculating the observation likelihood probability of the observed two-dimensional observation vector in each hidden state at time t by substituting the two-dimensional observation vector at the current time t into a pre-trained Gaussian mixture model include: For each of the three hidden states, an independent Gaussian mixture model is configured. The Gaussian mixture model contains multiple Gaussian components, and each Gaussian component has a preset mean vector, covariance matrix, and mixing coefficients.

[0079] Substitute the two-dimensional observation vector at the current time t into each Gaussian component of the Gaussian mixture model of the current hidden state, and calculate the single Gaussian probability density value of the observed two-dimensional observation vector under each Gaussian component.

[0080] The weighted probability density value of each Gaussian component is obtained by multiplying the single Gaussian probability density value of each Gaussian component by the corresponding mixing coefficient.

[0081] The observation likelihood probability of the two-dimensional observation vector at time t is obtained by summing the weighted probability density values ​​of all Gaussian components in the current hidden state.

[0082] Multiply the predicted probability of each hidden state at the current time t by the corresponding observation likelihood probability to obtain the unnormalized filtered probability of each hidden state at the current time t.

[0083] Calculate the sum of the unnormalized filtering probabilities of all hidden states at the current time t, and normalize the unnormalized filtering probability of each hidden state using the sum of the unnormalized filtering probabilities to obtain the filtering probability of each hidden state at the current time t.

[0084] Select the filter probability with the largest value from the filter probabilities of the three hidden states at the current time t, and take the hidden state corresponding to the maximum filter probability as the recognition state at the current time t.

[0085] For example, the observation emission probability model uses a Gaussian mixture model, with an independent Gaussian mixture model configured for each hidden state. Taking the Gaussian mixture model in the steady state as an example, in the steady state, both the resistivity descent rate and the excess pore water pressure ratio (after standardization) are close to zero, and the distribution of the two-dimensional observation vector is relatively concentrated, so a Gaussian mixture model with fewer Gaussian components can be used for modeling. In the structural failure-dominated state, the standardized resistivity descent rate shows a significant positive value while the standardized excess pore water pressure ratio remains low, and the distribution shape of the two-dimensional observation vector reflects this characteristic. In the pore pressure-structure coupled liquefaction state, both features increase significantly and fluctuate greatly, so the number of components in the Gaussian mixture model can be increased accordingly to improve the fitting accuracy. The mean vector of each Gaussian component is a two-dimensional vector, the covariance matrix is ​​a two-dimensional positive definite matrix, and the mixing coefficient is a scalar. The sum of the mixing coefficients of all Gaussian components equals 1. The mean vector, covariance matrix, and mixing coefficients are all obtained by performing an expectation-maximization algorithm on the historical standardized two-dimensional observation vector sequence belonging to the corresponding state in the training samples. The number of Gaussian components in each hidden state Gaussian mixture model can be automatically determined during the training phase based on the Bayesian information criterion. That is, the Bayesian information criterion value is calculated for Gaussian mixture models with different numbers of components on the training data, and the number of components with the smallest Bayesian information criterion value is selected as the optimal configuration.

[0086] During the real-time operation phase, for the current time t, the standardized two-dimensional observation vector at the current time is... Substitute the Gaussian mixture model corresponding to each hidden state and calculate the observation likelihood probability as follows. For a given hidden state... Its Gaussian mixture model includes There are k Gaussian components, and the mean vector of the k-th Gaussian component is... The covariance matrix is The mixing coefficient is Then the observation vector at the current time In hidden state The observed likelihood probability is: ;in, This is a two-dimensional Gaussian probability density function. Due to the dimension of the observation vector... Its expression is: Substitute back, The formula simplifies to: ;in, Let be the determinant of the covariance matrix. This is the inverse of the covariance matrix. The quadratic form in the exponential term. The dimensionless squared Mahalanobis distance is used, and the dimension of the entire probability density function is the reciprocal of the observation space volume. Integrating this value yields the dimensionless probability. Since the input observation vector has been standardized, all components are dimensionless, and the covariance matrix has good numerical conditions, ensuring the stability of the numerical calculation.

[0087] For the current time t, the predicted probability of each hidden state is calculated by combining the filtered probability of the previous time t-1 and the state transition probability matrix. Let the filtered probabilities of the three hidden states at the previous time t-1 be... , , The state transition probability matrix from state Transition to state The probability is Then the hidden state at the current time t The predicted probability is: Multiply the predicted probability of each hidden state at time t by the corresponding observed likelihood probability to obtain the unnormalized filtered probability: ; Calculate the sum of the unnormalized filtered probabilities of all hidden states. ,by Normalize the unnormalized filter probability for each hidden state to obtain the filter probability at the current time t: After normalization, it satisfies , The maximum value among the filtering probabilities of the three hidden states is selected, and its corresponding hidden state is taken as the recognition state at the current time t.

[0088] It should be noted that all parameters of the above Hidden Markov Model, including the initial state probability vector, the state transition probability matrix, and the standardized parameters (…), are considered to be… The mean vector, covariance matrix, and mixing coefficients of each hidden-state Gaussian mixture model were determined during the offline training phase before the system was officially put into use. Training data sources included resistivity and pore water pressure monitoring data collected synchronously during indoor triaxial dynamic tests of loess liquefaction, as well as historical field monitoring data of loess landslides with known evolution stages. During training, the Baum-Welch expectation-maximization algorithm was used to iteratively optimize the hidden Markov model parameters until the log-likelihood function converged, thus ensuring that the model parameters fully reflected the statistical characteristics of the three evolution stages of loess liquefaction landslides. After training, the model parameters were permanently stored in the early warning system and were no longer updated during the real-time monitoring phase; only forward filtering inference operations were performed.

[0089] S6. If the current identification state is a structural damage-dominated state, then issue a first-level warning signal; if the current identification state is a pore pressure-structure coupling liquefaction state, then issue a second-level warning signal; otherwise, do not issue a warning signal.

[0090] For example, if the current identification state is dominated by structural failure, a first-level warning signal is output, indicating that the loess slope's framework structure has shown signs of localized damage, and the cementation between soil particles has begun to break but has not yet developed into full liquefaction. Immediate monitoring of slope deformation and seepage should be strengthened, and preparations should be made to activate the personnel evacuation plan. If the current identification state is a pore pressure-structure coupled liquefaction state, a second-level warning signal is output, indicating that the slope is in a critical liquefaction instability stage where a large accumulation of excess pore water pressure and continuous deterioration of the framework structure mutually reinforce each other. The risk of a large-scale liquefaction landslide is extremely high, and personnel evacuation and emergency response measures should be initiated immediately. If the current identification state is stable, no warning signal is output, and the system remains in normal monitoring mode.

[0091] Example 2: Based on the same inventive concept, such as Figure 2 As shown in the figure, this embodiment also provides a dynamic early warning system for loess liquefaction landslides, the system comprising: The data sequence acquisition module is used to acquire resistivity time series and pore water pressure time series at different depths of loess slopes.

[0092] The resistivity signal processing module is used to decompose the resistivity time series, separate the baseline drift component, which characterizes the slow electrical changes caused by loess collapse, and the abrupt change component, which characterizes the sudden failure of the soil structure, and remove the baseline drift component; calculate the resistivity descent rate at the current moment based on the abrupt change component, where the resistivity descent rate is the relative magnitude of the resistivity decrease per unit time.

[0093] The module for calculating the excess pore water pressure ratio is used to calculate the excess pore water pressure ratio at the current moment based on the pore water pressure time series. The excess pore water pressure ratio is the ratio of the excess pore water pressure to the initial effective stress.

[0094] A two-dimensional observation sequence construction module is used to construct a two-dimensional observation vector by combining the resistivity drop rate with the ratio of excess pore water pressure, and to form a two-dimensional observation sequence through sliding sampling.

[0095] The hidden state identification module is used to calculate the filtering probability of each hidden state at the current time using a hidden Markov model on the two-dimensional observation sequence, and to use the hidden state corresponding to the maximum filtering probability as the identification state at the current time.

[0096] The warning signal issuing module is used to issue a first-level warning signal if the current identification state is a structural damage-dominated state; to issue a second-level warning signal if the current identification state is a pore pressure-structure coupling liquefaction state; otherwise, no warning signal is issued.

[0097] It should be noted that the specific methods by which each module performs operations in the system described in the above embodiments have been described in detail in the embodiments related to the method, and will not be elaborated here.

[0098] Finally, it should be noted that although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A dynamic early warning method for loess liquefaction landslides, characterized in that, The method includes: Obtain time series of resistivity and pore water pressure at different depths of loess slope; The resistivity time series is decomposed to separate the baseline drift component, which characterizes the slow electrical changes caused by loess collapse, and the abrupt change component, which characterizes the sudden failure of soil structure. The baseline drift component is then removed. The resistivity drop rate at the current moment is calculated based on the abrupt change component. The resistivity drop rate is the relative magnitude of the resistivity decrease per unit time. The excess pore water pressure ratio at the current moment is calculated based on the pore water pressure time series. The excess pore water pressure ratio is the ratio of excess pore water pressure to the initial effective stress. The resistivity drop rate and the ratio of excess pore water pressure are used to form a two-dimensional observation vector, and a two-dimensional observation sequence is formed by sliding sampling. The filtering probability of each hidden state at the current time is calculated using a hidden Markov model for the two-dimensional observation sequence, and the hidden state corresponding to the maximum filtering probability is taken as the identification state at the current time. If the current identification state is a structural failure-dominated state, a first-level warning signal is issued; if the current identification state is a pore pressure-structure coupling liquefaction state, a second-level warning signal is issued; otherwise, no warning signal is issued.

2. The dynamic early warning method for loess liquefaction landslides according to claim 1, characterized in that, The method for performing signal decomposition on the resistivity time series to separate the baseline drift component characterizing the slow electrical changes caused by loess collapse and the abrupt change component characterizing the sudden failure of soil structure includes: The resistivity time series is decomposed to obtain a series of intrinsic mode function components arranged from high frequency to low frequency and a residual term; Calculate the sample entropy or zero-crossing rate of each intrinsic mode function component. Identify intrinsic mode function components whose sample entropy is lower than the first entropy threshold or whose zero-crossing rate is lower than the first zero-crossing rate threshold as low-frequency components. Then, superimpose all low-frequency components with the residual term to obtain the baseline drift component that characterizes the slow electrical changes caused by loess collapsing. The eigenmode function components whose remaining sample entropy is higher than or equal to the first entropy threshold and whose zero-crossing rate is higher than or equal to the first zero-crossing rate threshold are reconstructed to obtain the mutation component characterizing the sudden failure of the soil structure.

3. The dynamic early warning method for loess liquefaction landslides according to claim 2, characterized in that, The method for decomposing the resistivity time series to obtain a series of eigenmode function components arranged from high frequency to low frequency and a residual term includes: Multiple sets of positive and negative paired adaptive white noise are added to the resistivity time series. Each set of positive and negative paired adaptive white noise constitutes a noise copy. Empirical mode decomposition is performed on each noise copy to obtain the first-order intrinsic mode function components corresponding to each noise copy. The first-order intrinsic mode function components of all noise replicas are ensembled and averaged to obtain the first-order intrinsic mode function components of the resistivity time series. The first-order intrinsic mode function components are then subtracted from the resistivity time series to obtain the first-order residual components. Add positive and negative paired adaptive white noise to the first-order residual component, repeat the empirical mode decomposition and ensemble averaging operations, and extract the second-order intrinsic mode function component, the third-order intrinsic mode function component, and so on up to the Nth-order intrinsic mode function component in sequence. The residual components that cannot be further decomposed are taken as residual terms, and the extracted first-order intrinsic mode function components to the Nth-order intrinsic mode function components are arranged in descending order of frequency to obtain a series of intrinsic mode function components arranged from high frequency to low frequency.

4. The dynamic early warning method for loess liquefaction landslides according to claim 3, characterized in that, The method for determining N includes: Perform a monotonicity test on the residual component of the current order obtained after each extraction. If the residual component of the current order is a monotonic function in its domain, stop decomposing and subtract one from the current order to obtain the value of N. If the current-order residual component is not a monotonic function, then calculate the permutation entropy of the current-order residual component, whereby the permutation entropy characterizes the complexity of the residual component. Based on the expected energy distribution of resistivity mutation components in loess slopes during earthquakes or rainfall events that induce liquefaction landslides, a signal entropy threshold is pre-set. If the permutation entropy of the current residual component is lower than the signal entropy threshold, it is determined that the current residual component no longer contains identifiable structural damage mutation information, and the decomposition is stopped. The current order is reduced by one to obtain the value of N. If the permutation entropy of the current order residual component is higher than or equal to the signal entropy threshold, then continue to execute the next order decomposition until the monotonicity or low entropy stopping condition is met.

5. The dynamic early warning method for loess liquefaction landslides according to claim 1, characterized in that, The method of constructing a two-dimensional observation vector by the resistivity drop ratio and the excess pore water pressure ratio, and forming a two-dimensional observation sequence by sliding sampling, includes: Set a time window length and a sliding step size; at each sampling time, combine the calculated resistivity drop rate at the current sampling time with the excess pore water pressure ratio at the current sampling time into a two-dimensional observation vector and store it; Each sampling time determined by the sliding step size is taken as the end time of the current window. Starting from the end time of the current window, the number of sampling points included in the time window length is traced back, and a two-dimensional observation vector sequence consisting of multiple consecutive two-dimensional observation vectors is extracted from the stored two-dimensional observation vectors.

6. A dynamic early warning method for loess liquefaction landslides according to claim 5, characterized in that, The method for setting a time window length and a sliding step size includes: Based on the typical evolution time of loess slopes from initial soil structure failure to liquefaction instability under rainfall or seismic loads, an initial window length covering the typical evolution time is determined as the time window length. Based on the data sampling period of resistivity drop rate and excess pore water pressure ratio, as well as the response speed requirements of the landslide early warning system, a sliding step size smaller than the time window length is set so that there is partial overlap between two adjacent sliding intercepted two-dimensional observation sequences. The time window length and sliding step size are converted into window length points and sliding step size points, respectively, in terms of the number of sampling points. The window length points are the integer obtained by dividing the time window length by the sampling period, and the sliding step size points are the integer obtained by dividing the sliding step size by the sampling period.

7. A dynamic early warning method for loess liquefaction landslides according to claim 1, characterized in that, The method of calculating the filtering probability of each hidden state at the current time using a Hidden Markov Model for the two-dimensional observation sequence, and using the hidden state corresponding to the maximum filtering probability as the identification state at the current time, includes: The three hidden states of the Hidden Markov Model are denoted as the steady state, the structural failure-dominated state, and the pore pressure-structure coupling liquefaction state, respectively. The filtering probabilities of each hidden state are initialized when the current time is zero. The filtering probability of the steady state is set to 1, and the filtering probabilities of the structural failure-dominated state and the pore pressure-structure coupling liquefaction state are both set to 0. For the current time t, obtain the filtered probability of each hidden state at the previous time t-1, multiply the filtered probability of each hidden state at the previous time t-1 by the state transition probability from it to each hidden state at the current time t, and sum the product results to obtain the predicted probability of each hidden state at the current time t. Substitute the two-dimensional observation vector at the current time t into the pre-trained Gaussian mixture model to calculate the observation likelihood probability of observing the two-dimensional observation vector in each hidden state at the current time t. Multiply the predicted probability of each hidden state at the current time t by the corresponding observation likelihood probability to obtain the unnormalized filtered probability of each hidden state at the current time t. Calculate the sum of the unnormalized filtering probabilities of all hidden states at the current time t, and normalize the unnormalized filtering probability of each hidden state using the sum of the unnormalized filtering probabilities to obtain the filtering probability of each hidden state at the current time t. Select the filter probability with the largest value from the filter probabilities of the three hidden states at the current time t, and take the hidden state corresponding to the maximum filter probability as the recognition state at the current time t.

8. A dynamic early warning method for loess liquefaction landslides according to claim 1, characterized in that, The method of substituting the two-dimensional observation vector at the current time t into a pre-trained Gaussian mixture model to calculate the observation likelihood probability of the observed two-dimensional observation vector in each hidden state at the current time t includes: For each of the three hidden states, an independent Gaussian mixture model is configured. The Gaussian mixture model contains multiple Gaussian components, and each Gaussian component has a preset mean vector, covariance matrix and mixing coefficients. Substitute the two-dimensional observation vector at the current time t into each Gaussian component in the Gaussian mixture model of the current hidden state, and calculate the single Gaussian probability density value of the observed two-dimensional observation vector under each Gaussian component. Multiply the single Gaussian probability density value of each Gaussian component by the corresponding mixing coefficient to obtain the weighted probability density value of each Gaussian component. The observation likelihood probability of the two-dimensional observation vector at time t is obtained by summing the weighted probability density values ​​of all Gaussian components in the current hidden state.

9. A dynamic early warning system for loess liquefaction landslides, characterized in that, The system includes: The data sequence acquisition module is used to acquire resistivity time series and pore water pressure time series at different depths of loess slopes. The resistivity signal processing module is used to decompose the resistivity time series, separate the baseline drift component, which characterizes the slow electrical changes caused by loess collapse, and the abrupt change component, which characterizes the sudden failure of the soil structure, and remove the baseline drift component; and calculate the resistivity descent rate at the current moment based on the abrupt change component, wherein the resistivity descent rate is the relative magnitude of the resistivity decrease per unit time. The module for calculating the excess pore water pressure ratio is used to calculate the excess pore water pressure ratio at the current moment based on the pore water pressure time series. The excess pore water pressure ratio is the ratio of the excess pore water pressure to the initial effective stress. A two-dimensional observation sequence construction module is used to construct a two-dimensional observation vector by the resistivity drop rate and the ratio of excess pore water pressure, and to form a two-dimensional observation sequence by sliding sampling. The hidden state identification module is used to calculate the filtering probability of each hidden state at the current time through the hidden Markov model of the two-dimensional observation sequence, and use the hidden state corresponding to the maximum filtering probability as the identification state at the current time. The warning signal issuing module is used to issue a first-level warning signal if the current identification state is a structural damage-dominated state; to issue a second-level warning signal if the current identification state is a pore pressure-structure coupling liquefaction state; otherwise, no warning signal is issued.