A method and system for early dynamic identification and monitoring of large-scale regional landslides

Abstract

The application discloses a large-range regional landslide early dynamic identification monitoring method and system, relates to the landslide dynamic identification technical field, and comprises the following steps: acquiring topographic monitoring data and radar monitoring data of a target grid area in a current monitoring period; based on the topographic monitoring data, performing physical validity screening and vector restoration processing on the radar monitoring data to obtain gravity flow direction displacement rates of different grid units in the target grid area; based on the gravity flow direction displacement rates, material looseness degrees and stress concentration degrees of the grid units, calculating a path penetration resistance map and corresponding optimal resistance dissipation path sets; and determining that the target grid area has a landslide risk in the case that corresponding indexes of the path penetration resistance map and the optimal resistance dissipation path sets meet preset conditions. The application improves the accuracy of early identification of hidden and penetrating landslide disasters.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of landslide dynamic identification, and particularly relates to a large-range regional landslide early dynamic identification monitoring method and system. BACKGROUND

[0002] The southwest mountainous area of China has complex terrain and geology, and is a region where landslide disasters frequently occur. At present, the main method for checking large-range landslide hazards is to obtain surface deformation information by using synthetic aperture radar interferometry (InSAR) technology. However, the existing monitoring methods face the following technical bottlenecks: first, the geometric projection deviation and noise interference are serious, the radar line-of-sight (LOS) observation is only one-dimensional projection of three-dimensional surface movement, and the traditional interpolation does not consider the anisotropy of the gravity flow direction of landslide movement, which easily introduces non-structural noise, and thus cannot accurately restore the real movement vector along the slope surface.

[0003] Secondly, it is difficult to identify hidden deep stress concentration. For the “pushing back and resisting forward” type of rock landslide, the precursors before sliding often show stress concentration and sliding surface penetration in deep rock-soil mass, and the surface deformation may not be significant, or it is difficult to distinguish temporary shallow creep caused by rainfall from critical instability.

[0004] The existing method lacks the comprehensive analysis capability of combining surface deformation, hydrological excitation caused terrain material loose degree and terrain gravity constraint, and cannot quantify the full-path penetration trend of the sliding surface in the slope body, which easily causes missed reports or false reports of hidden landslides. SUMMARY

[0005] In order to solve the technical problems in the related art that the comprehensive analysis capability of combining surface deformation, hydrological excitation caused terrain material loose degree and terrain gravity constraint is lacking, the full-path penetration trend of the sliding surface in the slope body cannot be quantified, and missed reports or false reports of hidden landslides are easily caused, the present application provides a large-range regional landslide early dynamic identification monitoring method and system.

[0006] The technical solutions adopted are as follows:

[0007] Obtain terrain monitoring data and radar monitoring data of a target grid area in a current monitoring period;

[0008] Based on the terrain monitoring data, perform physical validity screening and vector restoration processing on the radar monitoring data to obtain gravity flow direction displacement rates of different grid units in the target grid area, and the gravity flow direction displacement rate is used to represent the real movement rate of the surface along the terrain gravity flow direction;

[0009] The path penetration resistance map and the corresponding optimal resistance dissipation path set are calculated based on the gravity flow direction displacement rate, the material loose degree of each grid unit, and the stress concentration degree.

[0010] In a case where the path penetration resistance map and the corresponding optimal resistance dissipation path set satisfy a preset condition, it is determined that the target grid region has a landslide risk.

[0011] In a possible implementation of the present application, the radar monitoring data includes a deformation rate value and a direction unit vector of each grid unit under the radar line of sight.

[0012] Based on the terrain monitoring data, the radar monitoring data is subjected to physical validity screening and vector restoration processing to obtain the gravity flow direction displacement rate of different grid units in the target grid region, including:

[0013] Based on the digital elevation model and the terrain monitoring data, the terrain gravity flow direction vector of the target grid region is constructed.

[0014] The projection coefficient of the terrain gravity flow direction vector in the direction unit vector is calculated, and the radar monitoring data is subjected to physical validity screening processing based on the projection coefficient to obtain the first monitoring data after removing abnormal values.

[0015] The first monitoring data of each grid unit is subjected to vector restoration processing based on the projection coefficient and the deformation rate value to obtain the gravity flow direction displacement rate of different grid units in the target grid region.

[0016] In a possible implementation of the present application, the radar monitoring data is subjected to physical validity screening processing based on the projection coefficient to obtain the first monitoring data after removing abnormal values, including:

[0017] It is determined whether the projection sign and the projection coefficient in the radar monitoring data of each grid unit are consistent.

[0018] If they are not consistent, and the absolute value of the radar monitoring data is greater than a preset measurement error threshold, it is determined that the radar monitoring data is non-structural noise, and the data of the current grid unit is marked as an invalid value, so as to perform physical validity screening processing on the radar monitoring data to obtain the first monitoring data after removing abnormal values.

[0019] In a possible implementation of the present application, the path penetration resistance map and the corresponding optimal resistance dissipation path set are calculated based on the gravity flow direction displacement rate, the material loose degree of each grid unit, and the stress concentration degree, including:

[0020] The water level load data of the target grid region is obtained. Based on the water level load data and the gravity flow displacement rate, the displacement response coefficient is calculated. The displacement response coefficient is used to characterize the material looseness of the grid element.

[0021] For any grid cell, based on the difference in displacement response coefficients between the current grid cell and downstream neighbor cells in the preset gravity downstream neighbor index table, as well as the gravity flow direction displacement rate, the flow direction differential blocking index is calculated. The flow direction differential blocking index is used to characterize the stress concentration of the locking effect in the subsurface region.

[0022] Each grid cell is used as a node. For any node in the directed geoimpedance graph, a directed edge is established for the adjacent nodes in the preset gravity downstream neighborhood index table to construct the directed geoimpedance graph.

[0023] Based on the flow direction differential blocking index, displacement response coefficient, and Euclidean distance between two adjacent nodes in the geological impedance directional graph, the single-segment slip resistance value along the gravity flow direction of each directional edge is calculated.

[0024] By using a preset path search algorithm and single-segment slip resistance value, the shortest path integral value from any grid cell to the preset target node is searched on the geological impedance directed graph, resulting in a path penetration resistance graph and its corresponding optimal resistance dissipation path set.

[0025] In one possible implementation of this application, a flow-direction differential blocking index is calculated based on the difference in displacement response coefficients between the current grid cell and downstream neighbor cells in a preset gravity downstream neighbor index table, and the gravity flow direction displacement rate, including:

[0026] Calculate the average displacement response coefficient of the downstream neighboring cells in the preset gravity downstream neighboring index table of the current grid cell to obtain the average displacement.

[0027] The flow direction differential blocking index is calculated based on the difference between the displacement response coefficient and the average displacement of the current grid cell, and the product of the gravity flow direction displacement rate.

[0028] In one possible implementation of this application, the single-segment slip resistance value along each directional edge in the direction of gravity flow is calculated based on the flow direction differential blocking index, displacement response coefficient, and Euclidean distance between two adjacent nodes in the directed geoimpedance graph, including:

[0029] Calculate the weighted sum of the flow direction differential blocking index and the displacement response coefficient;

[0030] Based on the ratio of the Euclidean distance between two adjacent nodes in the directed geoimpedance graph to the weighted sum, the single-segment slip resistance value of each directed edge along the gravity flow direction is calculated. The single-segment slip resistance value is inversely proportional to the flow direction differential blocking index and the displacement response coefficient.

[0031] In one possible implementation of this application, by using a preset path search algorithm and a single-segment slip resistance value, the shortest path integral value from any grid cell to a preset target node is searched on the geological impedance directed graph to obtain a path penetration resistance graph and its corresponding optimal resistance dissipation path set, including:

[0032] Using a preset path search algorithm, the optimal path from any grid cell to the preset target node is searched on the geological impedance directed graph.

[0033] Calculate the minimum single-segment slip resistance value of each grid cell on the optimal path to obtain the shortest path integral value, and use the shortest path integral value as the path penetration resistance value of each grid cell.

[0034] Based on the path penetration resistance value and the optimal path corresponding to each grid cell, a path penetration resistance map and its corresponding optimal resistance dissipation path set are generated.

[0035] In one possible implementation of this application, if the indicators corresponding to the path connectivity resistance map and the optimal resistance dissipation path set meet preset conditions, it is determined that there is a landslide risk in the target grid area, including:

[0036] The path penetration resistance diagrams from the previous and current time points are subjected to time-series difference calculations to obtain the penetration resistance attenuation rate.

[0037] The degree of overlap between the optimal resistance dissipation path set at the current moment and the previous moment in terms of spatial distribution is calculated to obtain the degree of overlap of the dominant slip zone.

[0038] If the penetration resistance attenuation rate and the overlap of the dominant slip zone meet the preset conditions, it is determined that there is a landslide risk in the target grid area.

[0039] In one possible implementation of this application, when the penetration resistance attenuation rate and the overlap of dominant slip zones meet preset conditions, it is determined that there is a landslide risk in the target grid area, including:

[0040] If the penetration resistance attenuation rate exceeds the preset attenuation threshold and the overlap of the dominant slip zone exceeds the preset locking threshold, then the penetration resistance attenuation rate and the overlap of the dominant slip zone meet the preset conditions, and the target grid area is determined to have a landslide risk.

[0041] To achieve the above objectives, a large-scale regional landslide early dynamic identification and monitoring system is also provided. The system includes a memory, a processor, and a computer program stored in the memory and running on the processor. When the processor executes the computer program, it implements the steps of any of the methods described above.

[0042] This application has, but is not limited to, the following technical effects:

[0043] By acquiring topographic and radar monitoring data of the target grid area within the current monitoring period, and performing physical validity screening and vector reconstruction processing on the radar monitoring data based on the topographic monitoring data, the gravity flow displacement rate of different grid units is obtained. This determines the true movement rate along the topographic gravity flow direction. Furthermore, based on the gravity flow displacement rate, as well as the material looseness and stress concentration of each grid unit, a path penetration resistance map and the corresponding optimal resistance dissipation path set are determined. If the indicators corresponding to the path penetration resistance map and the optimal resistance dissipation path set meet preset conditions, it is determined that there is a landslide risk in the target grid area. In this application, the topographic monitoring data is first quantified to determine the true movement rate along the topographic gravity flow direction. Then, based on the gravity flow displacement rate, as well as the material looseness and stress concentration of each grid unit, the path penetration resistance map and the optimal resistance dissipation path set are calculated. This allows the monitoring results to reflect the overall topological connectivity of the landslide and effectively distinguish the essential difference between local soil loosening and overall slip surface penetration. Furthermore, by combining surface deformation, material looseness, and topographic gravity constraints, the full-path continuity trend of the slip surface inside the slope is quantified, which improves the accuracy of early identification of hidden and continuous landslide disasters and reduces the underreporting or false reporting of hidden landslides. Attached Figure Description

[0044] Figure 1 This is a flowchart illustrating the first embodiment of the method for early dynamic identification and monitoring of landslides in a large area according to this application.

[0045] Figure 2 This is a schematic diagram of the overall implementation process of the large-scale regional landslide early dynamic identification and monitoring method in this application;

[0046] Figure 3 This is a schematic diagram of the device structure of the hardware operating environment involved in the embodiments of this application. Detailed Implementation

[0047] It should be understood that the specific embodiments described herein are merely illustrative of this application and are not intended to limit this application.

[0048] This application provides a method for early dynamic identification and monitoring of landslides in a large area. In the first embodiment of this method, refer to...Figure 1 The methods include:

[0049] Step S10: Obtain terrain monitoring data and radar monitoring data of the target grid area within the current monitoring period.

[0050] As an example, the current monitoring cycle can be one hour, two hours, etc., without specific limitations. First, determine the digital elevation model (DEM) of the monitoring area. The target grid area can be defined based on the spatial resolution of the DEM data, defining the elevation grid of the monitoring area. The grid Include OK The grid cells of a column, where the coordinate index of any grid cell is denoted as... ,in The grid structure remains unchanged throughout all subsequent monitoring cycles, serving as a container for data alignment.

[0051] As an example, terrain monitoring data can be terrain data of the monitored target grid area, used to calculate the terrain gravity flow vector. .

[0052] As an example, radar monitoring data could be a terrain displacement rate map / radar line-of-sight displacement rate map in the current period along the radar line of sight. And the unit vector of the line-of-sight direction during radar satellite imaging. The radar line-of-sight displacement rate map is generated by an external InSAR processing system and includes the line-of-sight deformation rate scalar value / deformation rate value for each point within the monitoring area. .

[0053] As an example, the method for early dynamic identification and monitoring of landslides in large areas can also be applied to a system for early dynamic identification and monitoring of landslides in large areas. A schematic diagram of the overall execution flow of this system is shown below. Figure 2 As shown, the specific modules and their functions are as follows:

[0054] 1. Static benchmark construction and initialization module;

[0055] Function definition: Responsible for environment initialization during system cold start and constructing geometric and physical constraint benchmarks that can be reused throughout the entire lifecycle.

[0056] Core processing: Calculating the terrain gravity flow vector based on the DEM (Diagram of Earth's surface). ), and establish a gravity downstream neighborhood index ( Extracting topographic drainage baselines ( ) and calculate historical hydrological mean ( ).

[0057] Output: Static geometric constraint field, environment regularization parameters.

[0058] 2. Physical constraint data standardization module;

[0059] Function definition: Responsible for the access of dynamic observation data, screening of physical validity and vector reconstruction, and elimination of non-structural noise.

[0060] Core processing: Utilizing gravity flow vectors to analyze radar line-of-sight data / radar monitoring data ( Projection consistency verification is performed, anti-gravity noise and geometric blind spots are removed, and interpolation is used to generate a displacement rate map along the gravity flow direction. ), which is a velocity diagram composed of the gravity flow displacement rates of multiple grid cells.

[0061] Output: A spatially continuous real displacement rate field.

[0062] 3. Multiphysics decoupling and feature extraction module;

[0063] Functional definition: Responsible for separating material properties and stress state from apparent deformation, and identifying local kinematic blockages.

[0064] Core processing: Introducing water level load data ( ) Calculate the normalized displacement response coefficient / displacement response coefficient ( The flow direction differential blocking index is calculated based on the flow direction gradient. ).

[0065] Output: Material looseness field, local stress concentration field.

[0066] 4. Topological potential energy inversion and path planning module;

[0067] Functional definition: Responsible for constructing a geological impedance network and quantifying the overall continuity trend of the slip surface in three-dimensional space.

[0068] Core processing: Constructing a directed geoimpedance graph ( ), calculate the single-segment slip resistance ( The shortest path algorithm is used to calculate the total cumulative breakthrough resistance / path breakthrough resistance diagram. ) and extracting the optimal resistance dissipation path set ( ).

[0069] Output: Global resistance distribution and potential slip channel trajectory.

[0070] 5. Critical dynamic evolution identification module;

[0071] Functional definition: Responsible for the joint determination of temporal evolution analysis and critical instability state.

[0072] Core processing: Calculating the penetration resistance attenuation rate ( ) and the degree of overlap with the dominant slip zone ( ), and perform a double threshold AND gate decision.

[0073] Output: Early warning signal, landslide boundary prediction.

[0074] As an example, the data flow logic of the entire system is as follows:

[0075] Startup Phase (Phase 0): The system powers on, runs the static baseline construction module, reads the DEM, and generates static data. It is also resident in memory.

[0076] Dynamic Input (Phase 1): Entering the monitoring cycle External radar data access With hydrological data .

[0077] Data Standardization (Phase 2): The data standardization module calls static... right Perform cleaning and restoration, output .

[0078] Feature decoupling (Phase 3): Feature extraction module receives... and Output intermediate physical quantities and .

[0079] Potential energy inversion (Phase 4): The path planning module combines static topological constraints and intermediate physical quantities to calculate the core state variables. and .

[0080] State identification (Phase 5): The identification module reads the state of the previous cycle from the historical database. If there is no historical data (cold start), only the current state is stored; if there is, the differential index is calculated and a warning is issued.

[0081] Closed-loop update (Phase 6): Change the current state Write to the database, overwrite or archive, wait. cycle.

[0082] Step S20: Based on terrain monitoring data, perform physical validity screening and vector restoration processing on radar monitoring data to obtain the gravity flow displacement rate of different grid cells in the target grid area. The gravity flow displacement rate is used to characterize the actual movement rate of the ground surface along the terrain gravity flow direction.

[0083] As an example, this step aims to address the issues of geometric projection distortion, blind spots, and non-structural noise interference in radar line-of-sight observation data. Since Synthetic Aperture Radar (InSAR) can only acquire a one-dimensional projection of the three-dimensional surface motion along the radar line-of-sight, and this projection value is significantly affected by terrain undulations, directly using the raw observation data cannot accurately reflect the gravity-driven slope sliding state. Therefore, in this embodiment, a static geometric constraint model based on the terrain is first established. Then, the physical validity of the dynamically input radar monitoring data is screened and vectorized. Finally, a spatially continuous and standardized displacement rate map along the gravity flow direction, composed of individual grid cells, is output, providing a basic input for subsequent dynamic analysis.

[0084] As an example, the gravity flow displacement rate is used to characterize the actual movement rate of the surface along the gravity flow direction of the terrain. The gravity flow displacement rate map is the actual movement velocity matrix of the slope after vector reconstruction. Each element value in this matrix corresponds to the gravity flow displacement rate value of each grid cell.

[0085] The radar monitoring data includes the deformation rate value and direction unit vector of each grid cell under the radar line of sight; step S20 of the early dynamic identification and monitoring of landslides in a large area also includes steps S21 to S23, including:

[0086] Step S21: Based on the digital elevation model and terrain monitoring data, construct the terrain gravity flow vector of the target grid area.

[0087] As an example, the terrain gravity flow vector can be a unit vector representing the theoretical direction of matter movement determined by the terrain. Specifically, for a grid... Each grid cell in The slope angle was calculated using DEM elevation data. and slope angle Based on slope and aspect, a terrain gravity flow direction vector is constructed in a three-dimensional Cartesian coordinate system. It is a unit vector whose direction indicates the theoretical direction of the movement of surface material along the maximum slope line under the action of gravity alone; after the vector field is calculated, it is stored in the system's static database as a geometric benchmark for subsequent steps to judge the physical rationality of the observation data.

[0088] Step S22: Calculate the projection coefficient of the terrain gravity flow direction vector in the direction unit vector. Based on the projection coefficient, perform physical validity screening on the radar monitoring data to obtain the first monitoring data after removing outliers.

[0089] As an example, the projection coefficient is the topographic gravity flow vector. In the radar line of sight The projection component on the surface is used to characterize the sensitivity of the radar observation direction to the actual slope motion, for each grid cell in the target grid region. Projection coefficient The calculation method can be:

[0090]

[0091] in, Indicates the direction vector of terrain gravity flow. Represents a unit vector of direction.

[0092] Step S22 includes:

[0093] Determine whether the projection symbol and projection coefficient in the radar monitoring data of each grid cell are consistent;

[0094] If there is a discrepancy, and the absolute value of the radar monitoring data is greater than the preset measurement error threshold, the radar monitoring data is determined to be unstructured noise, and the data of the current grid cell is marked as invalid. This is to perform physical validity screening on the radar monitoring data and obtain the first monitoring data after removing outliers.

[0095] As an example, in the physical validity screening process, geometric blind zones are first identified, and then noise in the radar monitoring data of each grid cell is removed. The geometric blind zone identification process is as follows:

[0096] Since the projection coefficient approaches zero when the radar line of sight is perpendicular to the direction of slope movement, subsequent division to restore the value would lead to numerical explosion. Therefore, a projection angle threshold is set. (As an example, the threshold value can be...) to (Dimensionless values ​​between). Where, if the condition is satisfied... If the grid cell is located in the radar geometric observation blind zone, it will be placed in... The corresponding numerical value is marked as invalid (NaN).

[0097] The process of anti-gravity noise removal is as follows:

[0098] In landslide monitoring scenarios, the actual sliding motion should be downward along the slope, meaning the projection symbol along the radar line of sight should be consistent with... The signs of the observed values ​​are consistent with those of the projection coefficients. If the sign of the observed value is opposite to that of the projection coefficient, and the absolute value of the observed value exceeds the system's preset measurement error threshold, then the error is considered to be due to a mismatch between the observed value and the projection coefficient. (As an example, the threshold can be set to) to If the data point is not found to be non-structural noise caused by atmospheric turbulence or unwrapping errors, then the data in that grid cell is marked as invalid (NaN). If the observed value is opposite to the direction of gravity projection and its absolute value is less than the measurement error threshold, the velocity of that point is forcibly corrected to 0 instead of being marked as invalid. After the above processing, a sparse observation matrix after removing outliers is obtained, which is the first monitoring data.

[0099] Step S23: Based on the projection coefficient and deformation rate value, perform vector restoration processing on the first monitoring data of each grid cell to obtain the gravity flow displacement rate of different grid cells in the target grid region.

[0100] As an example, the purpose of vector restoration processing is to restore the one-dimensional scalar in the radar line-of-sight direction to the true motion vector magnitude along the slope and fill the gaps caused by screening in order to construct a spatially continuous velocity field.

[0101] As an example, for the target grid area Each one in the array is marked as a valid grid cell. The gravitational flow displacement rate is calculated using projection relationships. The calculation formula is:

[0102]

[0103] This calculation uses projection coefficients to restore the observations to the true slope slip component.

[0104] For regions of invalid values ​​(NaN) in the matrix (caused by geometric blind spots or noise removal), a spatial interpolation algorithm is used for numerical filling. One implementation method is inverse distance weighted (IDW) interpolation: for each invalid raster, valid raster values ​​within its specified radius neighborhood are searched; specifically, the search neighborhood is the terrain gravity flow vector. Within the sector, a weighted average is calculated using the reciprocal of the distance as the weight, which is then used as the estimate for that point.

[0105] Finally, based on the gravitational flow displacement rate of different grid cells, the current monitoring period is output. Displacement rate diagram along the direction of gravity flow The matrix covers the entire target grid area. It has spatial continuity, and all values ​​represent the actual motion rate along the direction of gravity flow in the terrain.

[0106] Step S30: Based on the gravity flow displacement rate, the material looseness of each grid cell, and the stress concentration, the path penetration resistance map and its corresponding optimal resistance dissipation path set are calculated. The path penetration resistance map is used to characterize the overall stability of the potential slip surface below the corresponding position of the grid cell at the current moment.

[0107] As an example, the path penetration resistance diagram can be the total energy / resistance matrix required to penetrate from a certain point to the toe of the slope. The smaller the path penetration resistance diagram, the easier it is for penetration failure to occur, which indicates that the overall stability of the potential slip surface is worse at the current moment. The optimal resistance dissipation path set can be the set of optimal paths from all points in the entire domain to the endpoint at the current moment.

[0108] As an example, the purpose of this step is to:

[0109] By decoupling external environmental excitation from soil material properties, the physical response characteristics of the in-situ medium are obtained. Then, the gradient difference concept in fluid mechanics is used to quantify the kinematic blocking effect. Finally, the stress concentration state is transformed into the passage cost in the topological network, and the penetration behavior of potential slip surfaces along the path of least resistance under gravity is simulated by graph theory algorithm, thereby identifying hidden critical instability regions.

[0110] Step S30 includes steps S31 to S35:

[0111] Step S31: Obtain water level load data for the target grid region. Based on the water level load data and gravity flow displacement rate, calculate the displacement response coefficient. The displacement response coefficient is used to characterize the material looseness of the grid element.

[0112] As an example, water level load data includes the current cumulative hydrological load and the historical average hydrological load, both of which can be obtained from historical meteorological databases. For example, the historical average hydrological load could be the average annual rainfall or the average groundwater level over the past 10 years. This data serves as a static scalar constant and is used as a regularization factor in subsequent normalization calculations.

[0113] As an example, based on the current monitoring cycle Displacement rate diagram along the direction of gravity flow and the current cumulative hydrological load For the target grid region Each grid cell in Calculate the displacement response coefficient The calculation formula is as follows:

[0114]

[0115] In the formula, This is an environmental regularization coefficient used to maintain the numerical stability of the denominator when the current hydrological load approaches zero (e.g., during the dry season), preventing the calculation results from diverging. As one implementation method, The value range is set to to Between. Based on the displacement response coefficient of each grid cell, a displacement response coefficient map of the entire target grid region is generated. The calculated The matrix characterizes the surface displacement response efficiency under unit hydrological load excitation; the larger the value, the looser the rock and soil structure at that location or the more sensitive it is to water.

[0116] Step S32: For any grid cell, based on the difference in displacement response coefficients between the current grid cell and downstream neighbor cells in the preset gravity downstream neighbor index table, and the gravity flow direction displacement rate, the flow direction differential blocking index is calculated. The flow direction differential blocking index is used to characterize the stress concentration of the locking effect in the subsurface region.

[0117] As an example, the preset gravity downstream neighborhood index table is based on the terrain gravity flow direction vector. The constructed downstream gravity neighborhood index table is pre-built with a global grid. For the target grid area any unit in Index table Stores the current grid cell All in the 8-neighborhood of Points to the coordinates of adjacent cells within the sector.

[0118] As an example, the purpose of this step is to identify kinematically incompatible regions, i.e., potential stress concentration areas, by analyzing the spatial gradient of the motion state.

[0119] Step S32 includes:

[0120] Calculate the average displacement response coefficient of the downstream neighboring cells in the preset gravity downstream neighboring index table of the current grid cell to obtain the average displacement.

[0121] As an example, for each grid cell Call the static gravity downstream neighborhood index table First, retrieve all downstream adjacent units contained in the index table. Read the displacement response coefficient And calculate the arithmetic mean of the response coefficients of these downstream units, denoted as the displacement average. .

[0122] The flow direction differential blocking index is calculated based on the difference between the displacement response coefficient and the average displacement of the current grid cell, and the product of the gravity flow direction displacement rate.

[0123] As an example, the response difference between the current cell and its downstream neighbor is calculated, and combined with the current motion magnitude, the flow-direction differential congestion index is calculated. The calculation formula is:

[0124] In this operation, An operation (equivalent to the modified linear unit ReLU) is used to implement unidirectional filtering: the difference is positive only when the response intensity of the current unit is significantly greater than that of the downstream unit (i.e., there is a squeezing effect along the direction of gravity), indicating blockage; otherwise, if the downstream movement is faster (tension effect), the difference is set to zero, because the tension zone usually does not have the mechanical characteristics of lock-up collapse. Simultaneously, the natural logarithm function is introduced. Nonlinear compression of the displacement rate modulus is applied to balance the weights of regions with different deformation orders of magnitude, preventing the numerical values ​​in the high-rate region from masking structural features. The calculation results form a flow-direction differential blockage index diagram. The high-value areas in this figure indicate strong stress concentration zones with deep locking effects. Among them, if the downstream neighborhood index table... If it is empty (i.e., located at the boundary), then let Equal to the current grid (Assuming downstream is unobstructed and there is no obstruction gradient), or mark this point. =0.

[0125] Step S33: Using each grid cell as a node, for any node in the directed geoimpedance graph, establish directed edges for adjacent nodes in the preset gravity downstream neighborhood index table to construct the directed geoimpedance graph.

[0126] As an example, step S33 uses graph theory to transform discrete physical properties into a continuous path planning space, simulating the process of the slip surface extending along the weak zone or breaking through the locked segment. Based on this, a directed geoimpedance graph is constructed using each grid cell as a node. The node set in this graph directly corresponds to the target grid region. All mesh cells in the array. Edge construction follows gravity flow constraints: for any node... Only index the downstream neighborhood of its gravity. Directed edges are established between adjacent nodes recorded in the data.

[0127] Step S34: Based on the flow direction differential blocking index, displacement response coefficient, and Euclidean distance between two adjacent nodes in the geological impedance directional graph, calculate the single-segment slip resistance value of each directional edge along the gravity flow direction.

[0128] Step S34 includes:

[0129] Calculate the weighted sum of the flow direction differential blocking index and the displacement response coefficient;

[0130] Based on the ratio of the Euclidean distance between two adjacent nodes in the directed geoimpedance graph to the weighted sum, the single-segment slip resistance value of each directed edge along the gravity flow direction is calculated. The single-segment slip resistance value is inversely proportional to the flow direction differential blocking index and the displacement response coefficient.

[0131] As an example, for each directed edge in the diagram, calculate its single-segment slip resistance value. This resistance value represents the energy cost required for the slip surface to penetrate between two adjacent nodes, and its value is set to be inversely proportional to the material looseness of the target node (characterized by the normalized displacement response coefficient) and the stress concentration (characterized by the flow differential blockage index). and Before calculation, the values ​​are standardized to the [0,1] interval using Min-Max, and then a weighted sum is performed. The calculation formula is as follows:

[0132]

[0133] In the formula, For grid cells Its downstream adjacent unit The Euclidean distance between the center points. and This is a weighting balancing coefficient used to adjust the contribution ratio of material properties and stress state in drag calculation. As one implementation method, both are normalized to... A constant within the interval. For example, a very small positive number (e.g., (), used to prevent the denominator from being zero. The formula represents the weighted sum. Its physical meaning is that when the target area is loosely textured or the current area has extremely high stress, the resistance to the expansion of the slip surface decreases, and the cost of penetration is reduced.

[0134] Step S35: Using a preset path search algorithm and single-segment slip resistance value, search the shortest path integral value from any grid cell to the preset target node on the geological impedance directed graph to obtain the path penetration resistance graph and its corresponding optimal resistance dissipation path set.

[0135] As an example, a global optimization algorithm is used to quantify the overall difficulty of penetrating the slip surface when any location within the entire domain is taken as the trailing edge point of a landslide.

[0136] As an example, the preset target node can be any endpoint in the defined global endpoint set. Specifically, based on hydrological analysis of the digital elevation model, gully lines and river areas are extracted and defined as topographic drainage baselines, serving as the global endpoint set for all potential slip paths.

[0137] Step S35 includes:

[0138] Using a preset path search algorithm, the optimal path from any grid cell to the preset target node is searched on the geological impedance directed graph.

[0139] Calculate the minimum single-segment slip resistance value of each grid cell on the optimal path to obtain the shortest path integral value, and use the shortest path integral value as the path penetration resistance value of each grid cell.

[0140] Based on the path penetration resistance value and the optimal path corresponding to each grid cell, a path penetration resistance map and its corresponding optimal resistance dissipation path set are generated.

[0141] As an example, the terrain drainage baseline The covered node set is the global target node set, and the cumulative impedance of these target nodes is initialized to zero in the geological impedance directed graph. Above, a single-source shortest path search algorithm / preset path search algorithm (as a preferred implementation, a variant of Dijkstra's algorithm or A* algorithm) is used to calculate the shortest path from each grid cell. Starting from the beginning, find the integral value of the shortest path along the directed edges to any global target node. This integral value represents the cumulative resistance to penetration along the entire path. The calculation process can be represented as finding an optimal path. This minimizes the sum of the single-segment slip resistance values ​​along the path:

[0142]

[0143] in, This represents the single-segment slip resistance value of the k-th directed edge, and the final output is the current monitoring period. Path penetration resistance diagram The cumulative penetration resistance along the entire path of each grid cell constitutes a path penetration resistance map. Each value in this matrix represents the overall stability of the potential slip surface below the corresponding location at the current moment; the lower the value, the more likely penetration failure will occur. Simultaneously, the optimal path node sequence corresponding to each grid cell is recorded, forming the optimal resistance dissipation path set. .

[0144] Step S40: If the indicators corresponding to the path penetration resistance map and the optimal resistance dissipation path set meet the preset conditions, it is determined that there is a landslide risk in the target grid area.

[0145] As an example, this step aims to address the technical challenge of distinguishing between steady-state creep and critical instability in large-scale landslide monitoring. Critical instability of a landslide is dynamically manifested as the irreversible loss of shear strength of the potential slip surface (energy dissipation) and the spatial fixation of the dominant slip channel (geometric locking). Based on the cumulative through-path resistance map and the optimal resistance dissipation path set, the critical precursors of a landslide are identified by monitoring the synchronicity of resistance decay and channel locking.

[0146] Step S40 includes:

[0147] The path penetration resistance attenuation rate is obtained by performing time-series difference calculation on the path penetration resistance diagrams of the previous and current time points.

[0148] As an example, determining the current monitoring period Does it have historical reference data? If so... If so, the current path penetration resistance map and the optimal resistance dissipation path set are stored in the system storage unit as the initial reference, and the current cycle process ends. If Then retrieve the previous monitoring cycle from the storage unit. Path penetration resistance diagram and the optimal resistance dissipation path set .

[0149] For each grid cell Perform the following calculations:

[0150] Calculate the attenuation rate of the penetration resistance This index characterizes the relative rate of loss of potential slip surface penetration resistance, and is calculated using the following formula:

[0151]

[0152] in, This represents the path penetration resistance value at the corresponding time. For example, the smallest positive number (e.g.) (), used to prevent the denominator from being zero. When When the value is positive and large, it indicates that the slip surface beneath the region is experiencing rapid resistance decay. It is the baseline value of the path penetration resistance after the sliding window smoothing process at the previous moment.

[0153] The degree of overlap between the optimal resistance dissipation path set at the current moment and the previous moment in terms of spatial distribution is calculated to obtain the degree of overlap of the dominant slip zone.

[0154] As an example, for mesh cells Advantageous slip zone overlap The calculation method can be:

[0155]

[0156] in, This represents the set of optimal resistance dissipation paths at time t (the current time). This represents the set of optimal drag dissipation paths at time t-1 (the previous time step), and this index characterizes the spatial stability of potential slip channels. It extracts the set of grid coordinates traversed by the optimal path corresponding to the current grid cell. and the corresponding set from the previous moment. The overlap ratio of dominant slip zones was calculated using the Jaccard coefficient.

[0157] in, Indicates the number of elements in the set. When Approaching When the time is right, it indicates that the spatial location of the potential slip channel is highly fixed and no longer wanders randomly. This is a typical geometric feature of stress concentration leading to the penetration of the shear band.

[0158] If the penetration resistance attenuation rate and the overlap of the dominant slip zone meet the preset conditions, it is determined that there is a landslide risk in the target grid area.

[0159] As an example, the attenuation rate of the penetration resistance and the overlap of the dominant slip zone are compared with their corresponding thresholds. Based on the comparison results, it is determined whether the preset conditions are met, and thus, the possibility of a landslide risk in the target grid area is determined. Through joint judgment logic, temporary acceleration interference caused by environmental factors such as rainfall is eliminated, and irreversible instability trends are accurately identified.

[0160] The step of determining that there is a landslide risk in the target grid area when the penetration resistance attenuation rate and the overlap of the dominant slip zone meet preset conditions includes:

[0161] If the penetration resistance attenuation rate exceeds the preset attenuation threshold and the overlap of the dominant slip zone exceeds the preset locking threshold, then the penetration resistance attenuation rate and the overlap of the dominant slip zone meet the preset conditions, and the target grid area is determined to have a landslide risk.

[0162] As an example, set a preset attenuation threshold. (As an example, the value is...) to ) and preset locking threshold (As an example, the value is...) to A point-by-point assessment of the global grid cells is performed, and the region is determined to be in a critical unstable state if and only if both of the following conditions are met simultaneously:

[0163] This indicates that the drag on the slip surface is being rapidly lost.

[0164] This indicates that the position of the sliding channel has been spatially locked.

[0165] The system generates early warning signals for grid cells that meet the above conditions and outputs the corresponding optimal resistance dissipation path as the predicted landslide boundary and sliding direction for disaster prevention departments to carry out emergency response. Finally, the current moment's... and Updated to the storage unit for analysis in the next monitoring cycle.

[0166] This application provides a method for early dynamic identification and monitoring of landslides in large-scale areas. It acquires topographic monitoring data and radar monitoring data of the target grid area within the current monitoring period. Based on the topographic monitoring data, it performs physical validity screening and vector reconstruction processing on the radar monitoring data to obtain the gravity flow displacement rate of different grid units. This determines the true movement rate along the topographic gravity flow direction. Furthermore, based on the gravity flow displacement rate, the material looseness, and stress concentration of each grid unit, it determines the path penetration resistance map and the corresponding optimal resistance dissipation path set. If the indicators corresponding to the path penetration resistance map and the optimal resistance dissipation path set meet preset conditions, it is determined that the target grid area has a landslide risk. In this application, the topographic monitoring data is first quantified to determine the true movement rate along the topographic gravity flow direction. Then, based on the gravity flow displacement rate, the material looseness, and stress concentration of each grid unit, the path penetration resistance map and the optimal resistance dissipation path set are calculated. This ensures that the monitoring results reflect the overall topological connectivity of the landslide and effectively distinguish the essential difference between local soil loosening and overall slip surface penetration. Furthermore, by combining surface deformation, material looseness, and topographic gravity constraints, the full-path continuity trend of the slip surface inside the slope is quantified, which improves the accuracy of early identification of hidden and continuous landslide disasters and reduces the underreporting or false reporting of hidden landslides.

[0167] Reference Figure 3 , Figure 3 This is a schematic diagram of the device structure of the hardware operating environment involved in the embodiments of this application.

[0168] like Figure 3As shown, the large-area landslide early dynamic identification and monitoring device may include: a processor 1001, a memory 1003, and a communication bus 1002. The communication bus 1002 is used to realize the connection and communication between the processor 1001 and the memory 1003.

[0169] Optionally, the large-scale landslide early dynamic identification and monitoring equipment may also include a user interface, a network interface, a camera, RF (Radio Frequency) circuitry, sensors, a WiFi module, etc. The user interface may include a display screen and an input submodule such as a keyboard; optional user interfaces may also include standard wired or wireless interfaces. The network interface may include standard wired or wireless interfaces (such as a Wi-Fi interface).

[0170] Those skilled in the art will understand that Figure 3 The structure of the large-area landslide early dynamic identification and monitoring equipment shown in the figure does not constitute a limitation on the large-area landslide early dynamic identification and monitoring equipment. It may include more or fewer components than shown in the figure, or combine certain components, or have different component arrangements.

[0171] like Figure 3 As shown, the memory 1003, serving as a storage medium, may include an operating system, a network communication module, and a large-scale landslide early dynamic identification and monitoring program. The operating system is a program that manages and controls the hardware and software resources of the large-scale landslide early dynamic identification and monitoring equipment, supporting the operation of the large-scale landslide early dynamic identification and monitoring program and other software and / or programs. The network communication module is used to enable communication between the various components within the memory 1003, as well as communication with other hardware and software in the large-scale landslide early dynamic identification and monitoring system.

[0172] exist Figure 3 In the large-scale landslide early dynamic identification and monitoring device shown, the processor 1001 is used to execute the large-scale landslide early dynamic identification and monitoring program stored in the memory 1003 to implement the steps of the large-scale landslide early dynamic identification and monitoring method described above.

[0173] The specific implementation method of the large-area landslide early dynamic identification and monitoring equipment in this application is basically the same as the embodiments of the above-mentioned large-area landslide early dynamic identification and monitoring method, and will not be repeated here.

[0174] It should be noted that, in this document, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or system that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or system. Unless otherwise specified, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or system that includes that element.

[0175] The sequence numbers of the embodiments in this application are for descriptive purposes only and do not represent the superiority or inferiority of the embodiments.

[0176] Through the above description of the embodiments, those skilled in the art can clearly understand that the methods of the above embodiments can be implemented by means of software plus necessary general-purpose hardware platforms. Of course, they can also be implemented by hardware, but in many cases the former is a better implementation method. Based on this understanding, the technical solution of this application, in essence, or the part that contributes to the prior art, can be embodied in the form of a software product. This computer software product is stored in a storage medium (such as ROM / RAM, magnetic disk, optical disk) as described above, and includes several instructions to cause a terminal device (which may be a mobile phone, computer, server, air conditioner, or network device, etc.) to execute the methods described in the various embodiments of this application.

[0177] The above are merely preferred embodiments of this application and do not limit the scope of this application. Any equivalent structural or procedural transformations made based on the description and drawings of this application, or direct or indirect applications in other related technical fields, are similarly included within the scope of protection of this application.

[0178] It should be noted that the order of the embodiments described above is merely for descriptive purposes and does not represent the superiority or inferiority of the embodiments. The processes depicted in the accompanying drawings do not necessarily require a specific or sequential order to achieve the desired result. In some embodiments, multitasking and parallel processing are also possible or may be advantageous.

[0179] The various embodiments in this specification are described in a progressive manner. The same or similar parts between the various embodiments can be referred to each other. Each embodiment focuses on describing the differences from other embodiments.

Claims

1. A method for early dynamic identification and monitoring of landslides in a large area, characterized in that, The method includes: Acquire terrain monitoring data and radar monitoring data for the target grid area within the current monitoring period; Based on the terrain monitoring data, the radar monitoring data is subjected to physical validity screening and vector restoration processing to obtain the gravity flow displacement rate of different grid cells in the target grid area. The gravity flow displacement rate is used to characterize the actual movement rate of the surface along the terrain gravity flow direction. Based on the gravity flow displacement rate, the material looseness of each mesh element, and the stress concentration, a path penetration resistance map and its corresponding optimal resistance dissipation path set are calculated, specifically including: Obtain water level load data for the target grid region. Based on the water level load data and the gravity flow displacement rate, calculate the displacement response coefficient, which is used to characterize the material looseness of the grid cells. For any grid cell, based on the difference in displacement response coefficients between the current grid cell and downstream neighbor cells in the preset gravity downstream neighbor index table, and the gravity flow direction displacement rate, the flow direction differential blocking index is calculated. The flow direction differential blocking index is used to characterize the stress concentration of the locking effect in the subsurface region. Each of the aforementioned grid cells is used as a node. For any given node, a directed edge is established for the adjacent nodes in the preset gravity downstream neighborhood index table to construct a directed graph of geological impedance. Based on the flow direction differential blocking index, displacement response coefficient, and the Euclidean distance between two adjacent nodes in the geological impedance directional graph, the single-segment slip resistance value along the gravity flow direction of each directional edge is calculated. Using a preset path search algorithm and the single-segment slip resistance value, the shortest path integral value from any grid cell to the preset target node is searched on the geological impedance directed graph to obtain the path penetration resistance graph and its corresponding optimal resistance dissipation path set. The path penetration resistance graph is used to characterize the overall stability of the potential slip surface below the corresponding position of the grid cell at the current moment. If the indicators corresponding to the path penetration resistance map and the optimal resistance dissipation path set meet the preset conditions, it is determined that there is a landslide risk in the target grid area.

2. The method for early dynamic identification and monitoring of landslides in large-scale areas as described in claim 1, characterized in that, The radar monitoring data includes the deformation rate value and direction unit vector of each grid cell under the radar line of sight; The process of performing physical validity screening and vector reconstruction on the radar monitoring data based on the terrain monitoring data to obtain the gravity flow displacement rate of different grid cells in the target grid region includes: Based on the digital elevation model and the terrain monitoring data, the terrain gravity flow vector of the target grid area is constructed; Calculate the projection coefficient of the terrain gravity flow direction vector onto the unit vector in the direction, and based on the projection coefficient, perform physical validity screening on the radar monitoring data to obtain the first monitoring data after removing outliers; Based on the projection coefficient and deformation rate value, the first monitoring data of each grid cell is vector-reconstructed to obtain the gravity flow displacement rate of different grid cells in the target grid region.

3. The method for early dynamic identification and monitoring of landslides in large-scale areas as described in claim 2, characterized in that, The step of performing physical validity screening on the radar monitoring data based on the projection coefficient to obtain the first monitoring data after removing outliers includes: Determine whether the projection symbol in the radar monitoring data of each grid cell is consistent with the projection coefficient; If there is a discrepancy, and the absolute value of the radar monitoring data is greater than the preset measurement error threshold, the radar monitoring data is determined to be unstructured noise, and the data of the current grid cell is marked as invalid. This is to perform physical validity screening on the radar monitoring data and obtain the first monitoring data after removing outliers.

4. The method for early dynamic identification and monitoring of large-scale landslides as described in claim 1, characterized in that, The calculation of the flow direction differential blocking index based on the difference in displacement response coefficients between the current grid cell and downstream neighbor cells in the preset gravity downstream neighbor index table, and the gravity flow direction displacement rate, includes: Calculate the average displacement response coefficient of the downstream neighboring cells in the preset gravity downstream neighboring index table of the current grid cell to obtain the average displacement. The flow direction differential blocking index is calculated based on the difference between the displacement response coefficient of the current grid cell and the average displacement value, and the product of the gravity flow direction displacement rate.

5. The method for early dynamic identification and monitoring of landslides in large-scale areas as described in claim 1, characterized in that, The single-segment slip resistance value along each directional edge in the direction of gravity flow is calculated based on the flow direction differential blocking index, displacement response coefficient, and Euclidean distance between two adjacent nodes in the directed geoimpedance graph, including: Calculate the weighted sum of the flow direction differential blocking index and the displacement response coefficient; Based on the ratio of the Euclidean distance between two adjacent nodes in the geoimpedance directed graph to the weighted sum, the single-segment slip resistance value of each directed edge in the direction of gravity flow is calculated, wherein the single-segment slip resistance value is inversely proportional to the flow direction differential blocking index and the displacement response coefficient.

6. The method for early dynamic identification and monitoring of landslides in large-scale areas as described in claim 1, characterized in that, The process involves searching the shortest path integral value from any grid cell to a preset target node on the geological impedance directed graph using a preset path search algorithm and the single-segment slip resistance value, thereby obtaining a path penetration resistance graph and its corresponding optimal resistance dissipation path set, including: Using a preset path search algorithm, the optimal path from any grid cell to the preset target node is searched on the geological impedance directed graph. Calculate the minimum single-segment slip resistance value of each grid cell on the optimal path to obtain the shortest path integral value, and use the shortest path integral value as the path penetration resistance value of each grid cell. Based on the path penetration resistance value and optimal path corresponding to each grid cell, a path penetration resistance map and its corresponding optimal resistance dissipation path set are generated.

7. The method for early dynamic identification and monitoring of landslides in large-scale areas as described in claim 1, characterized in that, The determination that there is a landslide risk in the target grid area when the corresponding indicators of the path penetration resistance map and the optimal resistance dissipation path set meet preset conditions includes: The path penetration resistance diagrams from the previous and current time points are subjected to time-series difference calculations to obtain the penetration resistance attenuation rate. The degree of overlap between the optimal resistance dissipation path set at the current moment and the previous moment in terms of spatial distribution is calculated to obtain the degree of overlap of the dominant slip zone. If the penetration resistance attenuation rate and the overlap of the dominant slip zone meet the preset conditions, it is determined that there is a landslide risk in the target grid area.

8. The method for early dynamic identification and monitoring of landslides in large-scale areas as described in claim 7, characterized in that, The determination that the target grid area has a landslide risk when the penetration resistance attenuation rate and the overlap of the dominant slip zone meet preset conditions includes: If the penetration resistance attenuation rate exceeds a preset attenuation threshold and the overlap of the dominant slip zone exceeds a preset locking threshold, then it is determined that the penetration resistance attenuation rate and the overlap of the dominant slip zone meet preset conditions, and the target grid area is judged to have a landslide risk.

9. A large-scale regional landslide early dynamic identification and monitoring system, characterized in that, The system includes a memory, a processor, and a computer program stored in the memory and running on the processor, wherein the processor executes the computer program to implement the steps of the method as described in any one of claims 1 to 8.

Images