Landslide sliding zone identification method and identification system based on three-dimensional coordinate solution

Abstract

The present application provides a landslide sliding zone identification method and identification system based on three-dimensional coordinate solution, relates to the surface deformation measurement technical field, and is characterized in that: a GNSS monitoring point is arranged in a landslide area, a reference point is arranged in an external area of the landslide, carrier phase and pseudo-range observation data of at least two sets of satellite navigation systems are collected, single-difference and double-difference preprocessing of the observation data is carried out, the three-dimensional coordinates of the monitoring point are solved, the three-dimensional coordinate time sequence of the monitoring point is generated in combination with the fixed whole-week ambiguity, the three-dimensional coordinate time sequence is converted to a local coordinate system of the slope body, the three-dimensional displacement vector, displacement rate and displacement gradient are solved, the triangular mesh unit is established by using a triangular subdivision algorithm, the deformation difference, rate deviation and gradient change amplitude of adjacent mesh are calculated, the isolated abnormal unit is removed, and the landslide sliding zone boundary is identified based on the continuous abnormal mesh profile, so that the defects that the traditional monitoring can only obtain the deformation of discrete points and cannot accurately define the sliding zone range are solved.

Description

Technical Field

[0001] This invention relates to the field of surface deformation measurement technology, specifically a method for identifying landslide slip zones based on three-dimensional coordinate calculation. Background Technology

[0002] Landslide disaster monitoring is an important part of geological disaster prevention and control. Global Navigation Satellite System (GNSS) positioning technology has been widely used in the field of landslide surface deformation monitoring due to its advantages such as all-weather operation, automation, and continuous observation.

[0003] In the prior art, CN109443188A discloses a two-layer multidimensional landslide monitoring method. This method constructs a first-layer displacement monitoring network between a base station and monitoring stations, and a second-layer attitude monitoring network between the monitoring stations. It uses a multi-mode receiver to collect raw observation data, combines the carrier double-difference observation equation and the extended Kalman filter algorithm to solve the coordinates and deformation of the monitoring stations, and solves the attitude angle of the landslide body based on the baseline vector between the monitoring stations. This method can simultaneously obtain the displacement information and attitude changes of the landslide body, reflecting the overall deformation trend of the landslide in multiple dimensions, and has good application value in engineering practice. However, its deformation is limited to two scales: discrete points and the overall macro scale. It can only qualitatively judge the overall deformation rate and overall deflection state of the landslide. It does not establish an automatic identification mechanism for concentrated sliding deformation areas, making it difficult to automatically delineate the precise spatial distribution boundary of the sliding zone and its evolution over time from the deformation data. This makes it difficult to support the spatial zoning assessment of landslide risk and targeted engineering treatment.

[0004] The information disclosed in the background section is only intended to enhance the understanding of the background of this disclosure, and therefore may include information that does not constitute prior art known to those skilled in the art. Summary of the Invention

[0005] The purpose of this invention is to provide a method for identifying landslide slip zones based on three-dimensional coordinate calculation, so as to solve the problems mentioned in the background art.

[0006] To achieve the above objectives, the present invention provides the following technical solution: The landslide slip zone identification method based on three-dimensional coordinate calculation includes the following steps: Step 1: Deploy multiple GNSS monitoring points in the landslide area and at least one reference point outside the landslide area. Collect multi-mode GNSS raw observation data from each GNSS monitoring point. The multi-mode GNSS raw observation data includes carrier phase observations and pseudorange observations from no less than two satellite navigation systems. Step 2: Combining the known fixed reference coordinates of the benchmark point with all the original observation data of multi-mode GNSS, perform joint differential calculation on all GNSS monitoring points and complete the carrier phase ambiguity fixation to obtain the time-series three-dimensional coordinates of each GNSS monitoring point. Step 3: Transform the three-dimensional coordinate time series of each GNSS monitoring point to the local coordinate system of the slope. Based on the coordinate difference between adjacent monitoring times, solve the three basic deformation characteristics of each GNSS monitoring point: three-dimensional displacement vector, displacement rate, and displacement gradient. Step 4: Divide the landslide area into grid units according to the spatial distribution of each GNSS monitoring point, determine the main sliding direction of the landslide by statistically analyzing the average direction of the three-dimensional displacement vector of each GNSS monitoring point, obtain the energy index to characterize the rock and soil stiffness of the corresponding grid unit through dynamic penetration test at the center point of each grid unit, and screen the unstable grid units by combining the surface crack development parameters of each grid unit. Step 5: Along the main sliding direction of the landslide, calculate the three types of difference features between the unstable grid cell and its adjacent downstream grid cell based on the three types of basic deformation characteristics, namely, the displacement difference amplitude, deformation rate deviation, and deformation gradient change amplitude, so as to construct a discrimination rule for identifying the sliding zone. Based on the discrimination rule, analyze the three types of difference features of the neighboring grid cells, identify the boundary of the sliding zone, and mark the corresponding grid cells.

[0007] Furthermore, the specific method for performing joint differential calculations on all GNSS monitoring points and fixing the carrier phase ambiguity to obtain the time-series three-dimensional coordinates of each GNSS monitoring point is as follows: Simultaneously extract carrier phase observations and pseudorange observations from each GNSS monitoring point and the reference point. Select one satellite as the reference satellite from all satellites observed by both, and use the remaining satellites as non-reference satellites. The difference between the carrier phase observations of the GNSS monitoring point and the reference point to the same satellite is calculated to obtain the inter-station single difference value of the carrier phase of the corresponding satellite; the difference between the pseudorange observations of the GNSS monitoring point and the reference point to the same satellite is calculated to obtain the inter-station single difference value of the pseudorange of the corresponding satellite. The difference between the carrier phase inter-station single difference value of each non-reference satellite and the carrier phase inter-station single difference value of the reference satellite is calculated to obtain the carrier phase inter-station-inter-satellite double difference value. The pseudo-range station-to-station single difference value of each non-reference satellite is calculated by comparing it with the pseudo-range station-to-station single difference value of the reference satellite to obtain the pseudo-range station-to-satellite double difference value. By using the pseudo-range inter-station-inter-satellite double difference and the reference coordinates of the benchmark point, the pseudo-range inter-station-inter-satellite double difference observation equation is constructed, and the approximate three-dimensional coordinates of the GNSS monitoring point are obtained by solving it. Using the approximate three-dimensional coordinates as constraints, and combining the carrier phase inter-station-inter-satellite double difference, a carrier phase inter-station-inter-satellite double difference observation equation is constructed. The integer ambiguity in the equation is searched, matched and fixed by the LAMBDA algorithm. Substitute the fixed integer ambiguity back into the carrier phase inter-station-inter-satellite double-difference observation equation to calculate the three-dimensional coordinate increment of the GNSS monitoring point relative to the reference point. By combining the known reference coordinates of the benchmark point with the three-dimensional coordinate increments, the three-dimensional coordinates of the GNSS monitoring points at each monitoring time are solved, and the coordinates are arranged in the order of monitoring time to generate a time sequence of the three-dimensional coordinates of each GNSS monitoring point.

[0008] Furthermore, based on the coordinate differences between adjacent monitoring times, the method for solving the three-dimensional displacement vector, displacement rate, and displacement gradient of each GNSS monitoring point is as follows: Three-dimensional coordinate time series data of two adjacent monitoring times of the same GNSS monitoring point are selected. The coordinates of the previous monitoring time are used as the initial reference and the coordinates of the next monitoring time are used as the deformed coordinates. The displacement change components of the GNSS monitoring point in three-dimensional space are obtained by coordinate difference calculation. The displacement components of each dimension are integrated to construct a three-dimensional displacement vector. Based on the time interval between adjacent monitoring times and the magnitude of the corresponding three-dimensional displacement vector, the displacement change per unit time is calculated to obtain the displacement rate of the GNSS monitoring point. Taking a single GNSS monitoring point as the center, and combining the three-dimensional displacement data of the surrounding adjacent GNSS monitoring points, the displacement spatial change rate at the location of the GNSS monitoring point is solved to obtain the corresponding displacement gradient; Traverse all GNSS monitoring points within the region, and solve for the three-dimensional displacement vector, displacement rate, and displacement gradient of each GNSS monitoring point according to the above calculation logic.

[0009] Furthermore, the method for dividing the landslide area into grid cells based on the spatial distribution of each GNSS monitoring point is as follows: Based on the planar coordinates of each GNSS monitoring point in the local coordinate system of the slope, the spatial distribution range of each GNSS monitoring point is determined. Using the planar position of each GNSS monitoring point as a node, the Delaunay triangulation algorithm is used to construct multiple irregular triangular grid cells covering the landslide area.

[0010] Furthermore, the method for obtaining the energy index of soil and rock stiffness for the corresponding grid unit through dynamic penetration testing is as follows: The center point of each grid cell was selected as the test site, and in-situ dynamic penetration tests were carried out. The drop weight, standard drop distance, probe cross-sectional area, and penetration depth were recorded during the test. Based on the mechanical mechanism of dynamic penetration tests, the single impact input energy was calculated according to the principle of work done by the drop weight. The single impact input energy is the product of the drop weight and the standard drop distance. The drop weight is obtained by multiplying the drop weight by the gravitational acceleration. The total impact input energy was obtained by accumulating the impact input energy of multiple impacts. The total input energy was divided by the product of the probe cross-sectional area and the actual penetration depth to obtain the impact input energy per unit penetration depth and per unit probe cross-sectional area, which is used as the energy index characterizing the soil stiffness of the corresponding grid cell.

[0011] Furthermore, the energy index characterizing the stiffness of the soil and rock is obtained for each grid cell. For grid cells with energy indices higher than a preset stiffness threshold, a screening rule is constructed based on the surface crack development parameters of each grid cell to screen unstable grid cells, as follows: Preset threshold values ​​for surface fissure development parameters, including threshold values ​​for average fissure width, fissure extension length, and fissure development density. If all surface fissure development parameters of a grid cell exceed the corresponding surface fissure development parameter threshold, it is determined to be an unstable grid cell. Traverse all grid cells in the landslide area and perform cell-by-cell discrimination according to the above screening rules to identify unstable grid cells.

[0012] Furthermore, the three types of difference characteristics among neighboring grid cells are calculated: displacement difference amplitude, deformation rate deviation, and deformation gradient change amplitude. The specific steps are as follows: Traverse all neighboring mesh elements and calculate the mean of the three-dimensional displacement vector, displacement rate, and displacement gradient for each mesh element as the corresponding three-dimensional displacement vector, displacement rate, and displacement gradient. Extract the three-dimensional displacement vector, displacement rate, and displacement gradient for each pair of adjacent mesh elements. Calculate the absolute difference between the magnitudes of the two sets of three-dimensional displacement vectors as the displacement difference amplitude between the two sets of mesh elements. Calculate the vector deviation between the two sets of displacement rates as the deformation rate deviation between the two sets of mesh elements. Calculate the difference between the magnitudes of the two sets of displacement gradients as the deformation gradient change amplitude between the two sets of mesh elements.

[0013] Furthermore, the method for constructing the discrimination rules for identifying the sliding band is as follows: For each unstable grid cell, a profile line parallel to the main slip direction is drawn through its center point. All grid cells along the profile line are extracted and arranged in order of elevation from high to low to form a process analysis grid sequence. Using the unstable grid cells within the analysis grid sequence as boundary nodes, a unified elevation datum is selected, and the elevation difference of each boundary node relative to the unified elevation datum is calculated. Simultaneously, three types of difference features are calculated between each boundary node and its adjacent downstream grid cells: displacement difference amplitude, deformation rate deviation, and deformation gradient change amplitude. The elevation differences, impact input energy, and corresponding three types of features obtained from all boundary nodes based on historical data are used to construct a training dataset. Using elevation difference and impact input energy in the training dataset as input features, and three types of difference features—displacement difference amplitude, deformation rate deviation, and deformation gradient change amplitude—as output labels, a regression algorithm is used to train three corresponding feature threshold prediction models. A preset threshold for the soil and rock stiffness energy index is set. Grid cells with an energy index higher than the threshold are designated as high-energy grid cells. Each high-energy grid cell is considered a boundary node. A series of analysis grids along the main sliding direction of the landslide are constructed. In the analysis grid sequence constructed at any boundary node, the elevation difference and impact input energy of the boundary node are used as inputs to the trained feature threshold prediction model to output three feature thresholds: displacement difference threshold, velocity deviation threshold, and gradient abrupt change threshold corresponding to the boundary node and its adjacent downstream grid cells. By monitoring the displacement difference amplitude, deformation rate deviation, and deformation gradient change amplitude of the grid cell corresponding to the boundary node and its downstream grid cell in real time, if the grid cell corresponding to the boundary node simultaneously meets the condition that all three deformation parameters exceed the threshold, the grid cell is determined to be a deformation anomalous grid. In the analysis grid sequence along the path of all deformation anomalous grids, the grid cells in the opposite direction to the main sliding direction are all determined to be the area where the sliding zone is located. Based on the outer contour of the grid cells in the area where the sliding zone is located, a closed spatial boundary is generated to delineate the distribution range of at least one landslide sliding zone, thus completing the delineation of the spatial boundary of the sliding zone. All grid cells belonging to the landslide sliding zone are uniformly marked.

[0014] Additionally, a landslide slip zone identification system based on three-dimensional coordinate calculation is provided. This system is used to execute any of the aforementioned landslide slip zone identification methods based on three-dimensional coordinate calculation, including: The reference module is used to set up multiple GNSS monitoring points in the landslide area and at least one reference point outside the landslide area to collect multi-mode GNSS raw observation data from each GNSS monitoring point. The multi-mode GNSS raw observation data includes carrier phase observations and pseudorange observations from no less than two satellite navigation systems. The data processing module is used to combine the known fixed reference coordinates of the benchmark point with all multi-mode GNSS raw observation data to perform joint differential calculation on all GNSS monitoring points, and to fix the carrier phase ambiguity, so as to obtain the time-series three-dimensional coordinates of each GNSS monitoring point. The feature calculation module is used to uniformly transform the three-dimensional coordinate time series of each GNSS monitoring point to the local coordinate system of the slope. Based on the coordinate difference between adjacent monitoring times, it solves the three-dimensional displacement vector, displacement rate and displacement gradient of each GNSS monitoring point respectively. The grid selection module is used to divide the landslide area into grid units according to the spatial distribution of each GNSS monitoring point, determine the main sliding direction of the landslide by statistically analyzing the average direction of the three-dimensional displacement vector of each GNSS monitoring point, obtain the energy index to characterize the rock and soil stiffness of the corresponding grid unit through dynamic penetration test at the center point of each grid unit, and select unstable grid units by combining the surface crack development parameters of each grid unit. The difference identification module is used to calculate three types of difference features between the unstable grid cell and its adjacent downstream grid cell based on three basic deformation features along the main sliding direction of the landslide: displacement difference amplitude, deformation rate deviation, and deformation gradient change amplitude. This is used to construct a discrimination rule for identifying the sliding zone. Based on this discrimination rule, the module analyzes the three types of difference features of the neighboring grid cells, identifies the boundary of the sliding zone, and marks the corresponding grid cells.

[0015] Compared with the prior art, the beneficial effects of the present invention are: This invention collects raw observation data from multi-mode GNSS systems containing at least two satellite navigation systems, performs joint differential calculations and ambiguity fixation, and can still stably output three-dimensional coordinate time-series data of each GNSS monitoring point even in scenarios where the landslide area is obscured by terrain and the number of visible satellites is insufficient. This allows for the calculation of three-dimensional displacement vectors, displacement rates, and displacement gradients. Subsequently, the landslide area is divided into grid cells based on the spatial distribution of each GNSS monitoring point. Three types of difference characteristics are obtained between adjacent grid cells: displacement difference amplitude, deformation rate deviation, and deformation gradient change amplitude. The main sliding direction of the landslide is determined by the average direction of the three-dimensional displacement vector. By analyzing the relationship between the energy index and the three types of difference characteristics, a discrimination rule is constructed along the main sliding direction of the landslide. Furthermore, by monitoring the three types of difference characteristics of adjacent grid cells in real time, the boundary of the sliding zone can be identified and its location marked, achieving automatic identification and spatial positioning capabilities of the sliding zone. Attached Figure Description

[0016] Figure 1 This is a schematic diagram of the overall method flow of the present invention; Figure 2 This is a schematic diagram of the overall system structure of the present invention. Detailed Implementation

[0017] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to specific embodiments.

[0018] It should be noted that, unless otherwise defined, the technical or scientific terms used in this invention should have the ordinary meaning understood by one of ordinary skill in the art to which this invention pertains. The terms "first," "second," and similar terms used in this invention do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Terms such as "comprising" or "including" mean that the element or object preceding the word encompasses the elements or objects listed following the word and their equivalents, without excluding other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. Terms such as "upper," "lower," "left," and "right" are used only to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0019] Example: Please see Figure 1 The present invention provides a technical solution: The landslide slip zone identification method based on three-dimensional coordinate calculation includes the following steps: Step 1: Deploy multiple GNSS monitoring points in the landslide area and at least one reference point outside the landslide area. Collect multi-mode GNSS raw observation data from each GNSS monitoring point. The multi-mode GNSS raw observation data includes carrier phase observations and pseudorange observations from no less than two satellite navigation systems. Significant spatial differences exist in the displacement amplitude, displacement direction, and deformation evolution characteristics of soil and rock masses at different locations on the landslide slope. Single-point monitoring data cannot fully characterize the deformation state of the entire landslide area. Therefore, multiple sets of multi-mode GNSS monitoring points are reasonably deployed throughout the entire landslide area to achieve dynamic monitoring of the entire landslide slope with full coverage and no blind spots, and to obtain deformation data of discrete monitoring points throughout the entire area.

[0020] For example, transverse monitoring lines are set up along the rear edge (tension zone), middle (main sliding zone), and front edge (shear / compression zone) of the main sliding direction of the landslide. The plane spacing of each GNSS monitoring point is controlled between 10 and 30m to ensure that the measuring points completely cover the entire slope of the landslide without large monitoring blank areas. In order to achieve high-precision differential positioning, at least one benchmark point is set up in a stable area outside the landslide body to ensure that its own coordinates do not change over time.

[0021] Traditional single-mode satellite navigation monitoring methods are susceptible to interference from multiple system and propagation errors, such as satellite orbit deviation, satellite clock error, upper-level ionospheric signal refraction delay, and lower-level tropospheric water vapor delay. In particular, landslide monitoring areas often have mountain obstructions, vegetation cover, and canyon terrain, making single-system observation data prone to discontinuities, jumps, and offsets. Therefore, this method adopts a multi-mode GNSS fusion observation approach. By accessing carrier phase and pseudorange observations from multiple satellite navigation systems, complementary calibration of multi-source data is achieved. Differential calculations weaken these major error sources, improving positioning accuracy from meter-level to millimeter-level. Landslide monitoring areas are often accompanied by complex terrain features such as mountain obstruction, tall vegetation cover, steep cliffs, and deep canyons. This leads to problems such as signal obstruction, observation link interruption, and data loss for a single satellite navigation system, seriously affecting the continuity of monitoring. At the same time, to further improve the stability of differential calculations and reduce the cumulative effect of observation errors, it is necessary to collect multi-mode GNSS raw observation data from each GNSS monitoring point. Specifically, the multi-mode GNSS raw observation data includes carrier phase and pseudorange observations from no less than two satellite navigation systems. Through complementary fusion of multi-system data, the integrity of monitoring data and the reliability of calculation results are ensured.

[0022] Step 2: Combining the known fixed reference coordinates of the benchmark point with all the original multi-mode GNSS observation data, perform joint differential calculation on all GNSS monitoring points and complete the carrier phase ambiguity fixation to obtain the time-series three-dimensional coordinates of each GNSS monitoring point. The specific method for performing joint differential calculations on all GNSS monitoring points and fixing the carrier phase ambiguity to obtain the time-series three-dimensional coordinates of each GNSS monitoring point is as follows: Simultaneously extract carrier phase observations and pseudorange observations from each GNSS monitoring point and the reference point. Select one satellite as the reference satellite from all satellites observed by both, and use the remaining satellites as non-reference satellites. The difference between the carrier phase observations of the GNSS monitoring point and the reference point to the same satellite is calculated to obtain the inter-station single difference value of the carrier phase of the corresponding satellite; the difference between the pseudorange observations of the GNSS monitoring point and the reference point to the same satellite is calculated to obtain the inter-station single difference value of the pseudorange of the corresponding satellite. The difference between the carrier phase inter-station single difference value of each non-reference satellite and the carrier phase inter-station single difference value of the reference satellite is calculated to obtain the carrier phase inter-station-inter-satellite double difference value. The pseudo-range station-to-station single difference value of each non-reference satellite is calculated by comparing it with the pseudo-range station-to-station single difference value of the reference satellite to obtain the pseudo-range station-to-satellite double difference value. The pseudorange single-point positioning observation equations of multiple satellites are solved simultaneously using the least squares method to obtain the initial approximate coordinates; The specific process of constructing and solving the pseudo-range inter-station-inter-satellite double-difference observation equation using the pseudo-range inter-station-inter-satellite double-difference values ​​and the reference coordinates of the benchmark point is as follows: Calculate pseudorange inter-station to inter-satellite geometric distances: in, Represents the three-dimensional spatial coordinates of the satellite. This represents the three-dimensional coordinates of the GNSS monitoring point to be determined. Indicates the reference coordinates of the benchmark point. This represents the geometric straight-line distance from the reference point to the satellite. This represents the geometric straight-line distance from the GNSS monitoring point to the satellite; Calculate the inter-station geometric single difference of pseudorange: In the formula, Indicates the geometric single difference between stations; For non-reference satellites, the inter-station geometric single difference is calculated based on the formula for calculating pseudorange: In the formula, This represents the geometric single difference between stations of a non-reference star. This represents the geometric straight-line distance from the GNSS monitoring point to the non-reference satellite. This represents the geometric straight-line distance from the reference point to a non-reference star; For the reference satellite, the inter-station geometric single difference is calculated based on the formula for calculating pseudorange: In the formula, This represents the geometric single difference between stations of the reference star. This represents the geometric straight-line distance from the GNSS monitoring point to the reference satellite. This represents the geometric straight-line distance from the reference point to the reference star; Calculate the inter-station-inter-satellite geometric distance double difference based on the non-reference satellite's inter-station geometric single difference and the reference satellite's inter-station geometric single difference: In the formula, This represents the double difference in the inter-station-inter-satellite geometric distance between the non-reference satellite and the reference satellite; Construct the linear equation for the measured pseudo-range station-to-satellite double-difference observations: In the formula, This represents the measured pseudorange inter-station / inter-satellite double difference. This represents the sum of pseudo-distance inter-station and inter-satellite double-difference observation noise and residuals; Establish a system of equations for all non-reference stars, and rearrange it into matrix-vector form: In the formula, This represents the coordinate correction value of the monitoring station, that is, the difference between the actual coordinates to be determined and the initial approximate coordinates of the monitoring station. The coefficient matrix represents the partial derivatives, which are formed by combining the partial derivatives of the inter-station and inter-satellite geometric distance differences at the initial approximate coordinates after performing a first-order Taylor expansion. This represents the pseudorange double-difference observation constant term, which is the measured pseudorange inter-station-inter-satellite double-difference value minus the theoretical inter-station-inter-satellite geometric distance double-difference value calculated from the initial approximate coordinates; Finally, the least squares method is used to solve the problem. The obtained coordinate correction increments, superimposed with the initial approximate coordinates, yield approximate three-dimensional coordinates of the GNSS monitoring point, where, This indicates the transpose of the coefficient matrix; Using this approximate three-dimensional coordinate system as a constraint, and combining the carrier phase inter-station-inter-satellite double difference values, the carrier phase inter-station-inter-satellite double difference observation equation is constructed as follows: In the formula, This represents the inter-station and inter-satellite double difference in carrier phase. This indicates the carrier wavelength corresponding to the satellite navigation signal. This represents the carrier phase double-difference integer ambiguity, i.e., the integer parameter to be solved. The sum of carrier phase observation noise and residuals is represented by the LAMBDA algorithm. The integer ambiguities in the equation are searched, matched and fixed. The fixed integer ambiguities are then substituted back into the carrier phase inter-station-inter-satellite double-difference observation equation. The three-dimensional coordinate increments of the GNSS monitoring points relative to the reference point are obtained by least squares solution. Then, by superimposing the known reference coordinates of the reference point, the three-dimensional coordinates of the GNSS monitoring points at each monitoring time are solved. The coordinates are then arranged in the order of monitoring time to generate the three-dimensional coordinate time sequence of each GNSS monitoring point.

[0023] Step 3: Transform the three-dimensional coordinate time series of each GNSS monitoring point to the local coordinate system of the slope. Based on the coordinate difference between adjacent monitoring times, solve the three basic deformation characteristics of each GNSS monitoring point: three-dimensional displacement vector, displacement rate, and displacement gradient. Since the coordinate axes of the global geodetic coordinate system do not conform to the slope sliding direction, the displacement values ​​cannot directly distinguish whether it is sliding along the slope, lateral displacement, or slope heave and settlement. Therefore, the three-dimensional coordinate time series of each GNSS monitoring point is uniformly transformed to the local coordinate system of the slope, so that the coordinate changes can characterize the slope bulging, collapse, tension crack deformation. The three-dimensional coordinate time series data of two adjacent monitoring times of the same GNSS monitoring point are selected. The coordinates of the previous monitoring time are used as the initial reference and the coordinates of the next monitoring time are used as the coordinates after deformation. The displacement change components of the GNSS monitoring point in three-dimensional space are obtained by calculating the coordinate difference. The displacement components of each dimension are integrated to construct a three-dimensional displacement vector. The three-dimensional components correspond to the main sliding displacement, lateral shear displacement, and vertical slope expansion and contraction displacement, respectively, which characterize the three-dimensional spatial deformation characteristics of the landslide. Based on the time interval between adjacent monitoring times and the magnitude of the corresponding three-dimensional displacement vector, the displacement change per unit time is calculated to obtain the displacement rate of the GNSS monitoring point. The displacement rate characterizes the deformation per unit time and represents different dangerous stages such as landslide stability creep, accelerated sliding, and sudden sliding. Centered on a single GNSS monitoring point, and combined with the three-dimensional displacement data of neighboring GNSS monitoring points, the displacement spatial change rate at the location of the GNSS monitoring point is calculated to obtain the corresponding displacement gradient. Since the displacement inside the slip zone is large and the gradient is gentle, while the displacement at the edge of the slip zone drops sharply and the gradient amplitude increases sharply under the critical state of the landslide, this displacement gradient is calculated to distinguish the boundary of the slip zone from the stable slope. Traverse all GNSS monitoring points within the region, and solve for the three basic deformation characteristics of each GNSS monitoring point—three-dimensional displacement vector, displacement rate, and displacement gradient—according to the above calculation logic.

[0024] Step 4: Divide the landslide area into grid units according to the spatial distribution of each GNSS monitoring point, determine the main sliding direction of the landslide by statistically analyzing the average direction of the three-dimensional displacement vector of each GNSS monitoring point, obtain the energy index to characterize the rock and soil stiffness of the corresponding grid unit through dynamic penetration test at the center point of each grid unit, and screen the unstable grid units by combining the surface crack development parameters of each grid unit. Extract the two-dimensional plane coordinates of all GNSS monitoring points in the local coordinate system of the slope, traverse all the plane coordinates of the measuring points to obtain the maximum and minimum values ​​of the abscissa and the maximum and minimum values ​​of the ordinate, and use these extreme values ​​to delineate the rectangular outer boundary, and use the rectangular outer boundary as the spatial distribution range of each GNSS monitoring point, so that the blank slope surface without measuring points is also completely covered by the grid cells. Using the planar positions of each GNSS monitoring point as nodes, the Delaunay triangulation algorithm is used to construct multiple irregular triangular mesh units covering the landslide area. The reason for using the Delaunay triangulation algorithm is that a single GNSS measuring point can only reflect the deformation of a single point, which is discrete data. The triangular mesh unit divides the entire landslide slope into continuous units. Each triangular mesh unit is constrained by three measured measuring points, which transforms discrete point data into spatially continuous data, realizing the visualization and quantitative analysis of landslide deformation across the entire area. In addition, the three vertices of a single triangle have their own three-dimensional displacement vector, displacement rate, and displacement gradient. The deformation difference amplitude and gradient abrupt change comparison between adjacent mesh units can be directly calculated using the triangular unit as the smallest calculation unit. The two-dimensional plane coordinates of all GNSS monitoring points in the local coordinate system of the slope, extracted at the latest monitoring time, are all constructed as two-dimensional displacement vectors. The average displacement vector is obtained by vector synthesis and averaging of all two-dimensional displacement vectors. The azimuth of the average displacement vector relative to true north is taken as the main sliding direction of the landslide. This direction represents the overall displacement trend of the landslide and will serve as the spatial reference benchmark for constructing the analysis profile along the path in subsequent steps. The logic of determining the main sliding direction of the landslide by the average direction of the three-dimensional displacement vectors of each GNSS monitoring point is as follows: the essence of a landslide is the shear sliding deformation of the slope along the weak structural surface. The overall sliding direction is consistent with the direction of the maximum sliding force of the slope and is the core factor that determines the distribution of slope deformation. By statistically averaging the displacement vectors of all monitoring points, the directional offset caused by local random deformation and boundary effects is offset, so that the main sliding direction has the representativeness of the whole area. In this scheme, the energy index, as a quasi-static background field characterizing the inherent mechanical properties of the soil and rock mass in each grid unit, needs to be acquired independently and periodically. The specific acquisition method is as follows: select the center point of each grid unit as the test point, conduct in-situ dynamic penetration tests, and record the hammer mass, standard drop distance, probe cross-sectional area, and penetration depth during the test; combine the mechanical mechanism of dynamic penetration test, calculate the single impact input energy according to the principle of work done by the hammer's own weight falling, the single impact input energy is the product of the hammer's weight and the standard drop distance, and the hammer's weight is obtained by multiplying the hammer's mass and gravitational acceleration; accumulate the impact input energy of multiple impacts to obtain the total impact input energy, divide the total input energy by the product of the probe cross-sectional area and the actual penetration depth to obtain the impact input energy per unit penetration depth and per unit probe cross-sectional area, which serves as the energy index characterizing the stiffness of the soil and rock. Among these, the dynamic penetration test is a test of the perturbation of the soil and rock mass, and the data is acquired periodically every 3 to 6 months. When the slip zone is located deep within the slope, surface cracks often reflect the deformation response of the shallow soil and rock mass. The distribution of cracks alone cannot directly and completely characterize the entire boundary of the deep slip zone. When surface cracks develop to a large scale, they manifest as tension crack zones, shear crack zones, or feather-like crack clusters. These cracks are the projected boundaries of the deep slip zone on the surface. Therefore, the significance of cracks is used as the basis for judging the unstable grid unit. At the same time, in the landslide mechanics mechanism, the slip zone is not a uniform and continuous weak surface, but is locked by a series of concave and convex bodies with different strengths. Among them, the concave and convex bodies that bear the main resistance and play a controlling role in slope stability are called critical concave and convex bodies. When the surface crack development parameters (average crack width, extension length, development density) exceed a certain value, it indicates that the high-stiffness critical concave and convex bodies have broken and cracked. The slip zone will be locally connected at this location. At this time, the deformation difference characteristics will show a numerical abrupt change, which is the physical law of the rapid decrease in bearing capacity when the high-stiffness medium breaks. Based on the above understanding of the mechanism, it can be known that grid cells with higher energy indices correspond to key concave-convex features on the slip surface, while exceeding the limit of crack development parameters indicates that the key concave-convex features have been destroyed and the slip surface has been partially connected. Therefore, it is necessary to continuously monitor grid cells with higher energy indices, and when the surface cracks develop to a large scale, they should be marked as unstable grid cells. The specific screening method is as follows: The energy index representing the stiffness of soil and rock is extracted from each grid cell. A threshold for the energy index of soil and rock stiffness is preset. Grid cells with an energy index higher than the threshold of the energy index of soil and rock stiffness are regarded as high-energy grid cells. For grid cells with an energy index higher than the preset stiffness threshold, each high-energy grid cell is regarded as a boundary node. Combined with the surface crack development parameters at each monitoring time, a surface damage feature screening rule is constructed. A preset threshold for surface crack development parameters is established, which includes the average crack width, crack extension length, and crack development density. The crack development density is the total length of cracks within a unit grid area. The specific method for obtaining and calibrating these threshold parameters is as follows: historical displacement rate data from three monitoring points within each high-energy grid cell are extracted. The statistical mean and standard deviation of the displacement rate of that grid cell are calculated. The mean is multiplied by three times the standard deviation as the rate mutation threshold. When the displacement rate of a high-energy grid cell exceeds this threshold, a damage verification mechanism is triggered for that grid cell. This involves conducting an in-situ dynamic penetration test on that grid cell to obtain the dynamic penetration energy index. If this dynamic penetration energy index is lower than the soil stiffness energy... A threshold is set to determine whether a grid cell has changed from a high-energy grid to a non-high-energy grid, indicating that the key type of concavity / convexity at the grid cell has undergone macroscopic fracturing. The surface crack development parameters corresponding to the grid cell are recorded, including the average crack width, crack extension length, and crack development density. A fracture crack feature sample set is constructed from the surface crack development parameters corresponding to at least five grid cells. The sample mean and standard deviation of the surface crack development parameters are calculated based on the fracture crack feature sample set. The result of subtracting twice the sample standard deviation from the sample mean is used as the corresponding surface crack development parameter threshold. For example, in this embodiment, the average crack width threshold is set to 15 mm, the crack extension length threshold is set to 3 m, and the crack development density threshold is set to 0.3 m / m². At any given monitoring time, if the surface crack development parameters of a grid cell all exceed the corresponding surface crack development parameter threshold, the grid cell is determined to be an unstable grid cell. Among them, in the method of setting the threshold of the rock and soil stiffness energy index, the rock and soil types, densities, weathering degrees and particle compositions of different landslides vary greatly, and the corresponding energy index magnitudes can range from several times to even an order of magnitude. The location where the surface crack development exceeds the limit is essentially the direct response of the key concave-convex body to the failure of the locking function on the surface. The energy index of such locations directly quantifies the inherent stiffness level of the key concave-convex body that undertakes the core anti-sliding function. At the same time, the core characteristic of the key concave-convex body is that its stiffness is significantly higher than that of the surrounding ordinary rock and soil. Therefore, for grid cells where the average crack width, extension length and development density all significantly exceed the standard, the previously obtained energy index is extracted, and its minimum value is used as the threshold of the rock and soil stiffness energy index.

[0025] Step 5: Along the main sliding direction of the landslide, calculate the three types of difference features between the unstable grid cell and its adjacent downstream grid cell based on the three types of basic deformation characteristics, namely, the displacement difference amplitude, deformation rate deviation, and deformation gradient change amplitude, so as to construct a discrimination rule for identifying the sliding zone. Based on the discrimination rule, analyze the three types of difference features of the neighboring grid cells, identify the boundary of the sliding zone, and mark the corresponding grid cells. By traversing all pairwise adjacent triangular mesh cells, the three-dimensional displacement vector, displacement rate, and displacement gradient—three types of deformation parameters—are extracted from each pair of adjacent mesh cells. Multi-dimensional deformation difference quantification is then performed. The specific principles and calculation logic are as follows: The three-dimensional displacement vector represents the cumulative total deformation of the mesh cell from the start of monitoring to the current time. Significant differences exist in the deformation state of the soil and rock on both sides of the landslide sliding zone boundary, and the cumulative deformation is prone to abrupt changes. Therefore, by calculating the absolute difference in the magnitude of the three-dimensional displacement vectors of adjacent mesh cells, the cumulative total deformation difference of the soil and rock in adjacent areas is quantified, and this absolute difference is defined as the displacement difference amplitude of adjacent mesh cells. The displacement rate represents the deformation amplitude of the soil and rock per unit time, which can intuitively reflect the landslide... The real-time sliding activity of the slope surface is assessed. The landslide slip zone often exhibits continuous creep and accelerated sliding characteristics, with a displacement rate significantly higher than the surrounding stable area. By calculating the vector deviation of the displacement rates of adjacent grid cells, the degree of differentiation in the real-time sliding state of the soil in adjacent grid areas can be accurately determined, yielding the corresponding deformation rate deviation. The displacement gradient, representing the spatial rate of change of slope displacement, characterizes the intensity of the spatial transition in slope deformation. Deformation within the slip zone is continuous and gradual, resulting in a smaller displacement gradient value. However, the edge of the slip zone represents abrupt deformation interfaces, where the displacement gradient exhibits a sharp increase. By calculating the difference in the magnitude of the displacement gradient between adjacent grid cells, the local abrupt deformation transition area of ​​the slope can be effectively identified. This difference is defined as the magnitude of the deformation gradient change between adjacent grid cells.

[0026] For each unstable grid cell, a profile line parallel to the main slip direction is drawn through its grid center point. All grid cells along the profile line are extracted and arranged in order of elevation from high to low to form a process analysis grid sequence. Using the unstable grid cells within the analysis grid sequence along the landslide as boundary nodes, a unified elevation benchmark is selected, such as the average elevation of all grid center points on the entire landslide slope or the highest point of the slope top as a unified benchmark. The elevation difference of each boundary node relative to this unified elevation benchmark is calculated. Simultaneously, three types of difference features are calculated between each boundary node and its adjacent downstream grid cells: displacement difference amplitude, deformation rate deviation, and deformation gradient change amplitude. The elevation differences, impact input energy, and corresponding three types of features obtained from all boundary nodes based on historical data are used to construct a training dataset. The unstable grid corresponds to the location on the sliding surface where the key concave-convex body is in the failure stage and the locking effect fails. Before the locking effect fails, the key concave-convex body bears the main anti-sliding force, and the soil on the downstream side is blocked and constrained, with small displacement, low deformation rate, and gentle displacement gradient. When the key concave-convex body fails, the locking force is released instantly, the sliding constraint of the soil on the downstream side is released, and the displacement and velocity will increase stepwise, and the displacement gradient will also change abruptly. Therefore, using the elevation difference and impact input energy in the dataset as input features, and the displacement difference amplitude, deformation rate deviation, and deformation gradient change amplitude as output labels, respectively, a regression algorithm is used to train three corresponding feature threshold prediction models. This enables the prediction of the critical amplitude of three types of deformation change that should occur on the downstream side when the key concave-convex body fails at any location of the landslide grid cell, which is the adaptive judgment threshold for the corresponding location. By analyzing the relationship between the energy index and three types of differential characteristics, a discrimination rule is constructed along the main sliding direction of the landslide. The logic of generating the feature threshold prediction model is as follows: the sliding deformation and stress transmission of the landslide develop gradually along the main sliding direction, and the penetration of the sliding zone and the damage of the locking concave and convex bodies also occur sequentially along the sliding path. Constructing a profile along the main sliding direction can capture the evolution law of deformation amount, deformation rate and deformation gradient along the process, and locate the boundary position of deformation abrupt change. The adjacent downstream grid cell is defined with the main sliding direction as a reference, the direction of sliding forward is downstream, and the lower elevation side is the downstream adjacent grid.

[0027] At the monitoring time, the three types of differences between each pair of adjacent grid cells are compared one by one. All neighboring grid cells are traversed. The three-dimensional displacement vector, displacement rate and displacement gradient mean of each grid cell are calculated as the three-dimensional displacement vector, displacement rate and displacement gradient of the corresponding grid cell. The displacement difference amplitude, deformation rate deviation and deformation gradient change amplitude between grid cells are extracted respectively. Each high-energy grid cell is used to construct a process analysis grid sequence along the main sliding direction of the landslide. In the analysis grid sequence constructed by any high-energy grid cell, the elevation difference and impact input energy of the high-energy grid cell are used as inputs to the trained feature threshold prediction model to output the three feature thresholds: displacement difference threshold, rate deviation threshold and gradient change threshold corresponding to the high-energy grid cell and its adjacent downstream grid cells. By real-time monitoring of the displacement difference amplitude, deformation rate deviation, and deformation gradient change amplitude of the high-energy grid cell and its downstream grid cell, if the high-energy grid cell simultaneously meets the threshold condition of all three types of difference characteristics, the grid cell is determined to be a deformation anomalous grid. In the analysis grid sequence along the slip direction of all deformation anomalous grid cells, the grid cells in the opposite direction to the main slip direction are all determined to be the slip zone region. The reason for identifying all grid cells in the opposite direction of the main sliding direction as the location of the sliding zone is that, in engineering, the sliding zone refers to a region of rock and soil that is continuously distributed along the sliding surface and undergoes overall sliding, rather than a discrete single-point failure location. The boundary of the landslide sliding zone is controlled by a series of key high-stiffness concave-convex bodies on the sliding zone. Only these bodies, whose own strength is much higher than that of the surrounding rock and soil, will form an observable deformation boundary before and after failure. That is, the displacement rate of the sliding zone will increase stepwise with the downstream grid cells. Grid cells of ordinary stiffness correspond to homogeneous rock and soil, and the deformation along the main sliding direction is a continuous and gradual transition. Observing the features extracted based on the three basic deformation features of three-dimensional displacement vector, displacement rate and displacement gradient will include a large number of small rate fluctuations in normal gradual transition areas, resulting in a large error in the sliding zone identification result. Tracing back to the upstream area covering the entire profile along the opposite direction of the main sliding direction as the location of the sliding zone is consistent with its continuous distribution. Connectivity analysis is performed on the global deformation anomaly grid cells to divide spatially adjacent and contiguous anomaly grid cells into the same anomaly connected domain, and isolated single-point anomaly grid cells without adjacent relationships are removed. For each contiguous anomaly grid connected domain, the outer edge segments of the outermost triangular grid cells of the connected domain are extracted, and all outer edge segments are spliced ​​and fitted to generate a closed polygon contour. This closed contour is the spatial distribution boundary of the landslide sliding zone, and finally the distribution range of at least one landslide sliding zone is delineated, completing the spatial boundary delineation of the sliding zone. All triangular grid cells falling inside the closed contour are uniformly marked.

[0028] Please see Figure 2 The present invention further provides a landslide slip zone identification system based on three-dimensional coordinate calculation, used to execute the landslide slip zone identification method based on three-dimensional coordinate calculation described in any of the above claims, comprising: The reference module is used to set up multiple GNSS monitoring points in the landslide area and at least one reference point outside the landslide area to collect multi-mode GNSS raw observation data from each GNSS monitoring point. The multi-mode GNSS raw observation data includes carrier phase observations and pseudorange observations from no less than two satellite navigation systems. The data processing module is used to combine the known fixed reference coordinates of the benchmark point with all multi-mode GNSS raw observation data to perform joint differential calculation on all GNSS monitoring points, and to fix the carrier phase ambiguity, so as to obtain the time-series three-dimensional coordinates of each GNSS monitoring point. The feature calculation module is used to uniformly transform the three-dimensional coordinate time series of each GNSS monitoring point to the local coordinate system of the slope. Based on the coordinate difference between adjacent monitoring times, it solves the three-dimensional displacement vector, displacement rate and displacement gradient of each GNSS monitoring point respectively. The grid selection module is used to divide the landslide area into grid units according to the spatial distribution of each GNSS monitoring point, determine the main sliding direction of the landslide by statistically analyzing the average direction of the three-dimensional displacement vector of each GNSS monitoring point, obtain the energy index to characterize the rock and soil stiffness of the corresponding grid unit through dynamic penetration test at the center point of each grid unit, and select unstable grid units by combining the surface crack development parameters of each grid unit. The difference identification module is used to calculate three types of difference features between the unstable grid cell and its adjacent downstream grid cell based on three basic deformation features along the main sliding direction of the landslide: displacement difference amplitude, deformation rate deviation, and deformation gradient change amplitude. This is used to construct a discrimination rule for identifying the sliding zone. Based on this discrimination rule, the module analyzes the three types of difference features of the neighboring grid cells, identifies the boundary of the sliding zone, and marks the corresponding grid cells.

[0029] The above formulas are all dimensionless calculations. The formulas are derived from software simulations based on a large amount of collected data to obtain the most recent real-world results. The preset parameters in the formulas are set by those skilled in the art according to the actual situation.

[0030] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented in software, the above embodiments can be implemented, in whole or in part, as a computer program product. Those skilled in the art will recognize that the units and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented by electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution.

[0031] The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; they may be located in one place or distributed across multiple network units. Some or all of the units can be selected to achieve the purpose of this embodiment, depending on actual needs.

[0032] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.

Claims

1. A method for identifying landslide slip zones based on three-dimensional coordinate calculation, characterized in that, The specific steps include: Step 1: Deploy multiple GNSS monitoring points in the landslide area and at least one reference point outside the landslide area. Collect multi-mode GNSS raw observation data from each GNSS monitoring point. The multi-mode GNSS raw observation data includes carrier phase observations and pseudorange observations from no less than two satellite navigation systems. Step 2: Combining the known fixed reference coordinates of the benchmark point with all the original observation data of multi-mode GNSS, perform joint differential calculation on all GNSS monitoring points and complete the carrier phase ambiguity fixation to obtain the time-series three-dimensional coordinates of each GNSS monitoring point. Step 3: Transform the three-dimensional coordinate time series of each GNSS monitoring point to the local coordinate system of the slope. Based on the coordinate difference between adjacent monitoring times, solve the three basic deformation characteristics of each GNSS monitoring point: three-dimensional displacement vector, displacement rate, and displacement gradient. Step 4: Divide the landslide area into grid units according to the spatial distribution of each GNSS monitoring point, determine the main sliding direction of the landslide by statistically analyzing the average direction of the three-dimensional displacement vector of each GNSS monitoring point, obtain the energy index to characterize the rock and soil stiffness of the corresponding grid unit through dynamic penetration test at the center point of each grid unit, and screen the unstable grid units by combining the surface crack development parameters of each grid unit. Step 5: Along the main sliding direction of the landslide, calculate the three types of difference features between the unstable grid cell and its adjacent downstream grid cell based on the three types of basic deformation characteristics, namely, the displacement difference amplitude, deformation rate deviation, and deformation gradient change amplitude, so as to construct a discrimination rule for identifying the sliding zone. Based on the discrimination rule, analyze the three types of difference features of the neighboring grid cells, identify the boundary of the sliding zone, and mark the corresponding grid cells. The center point of each grid cell was selected as the test site, and in-situ dynamic penetration tests were carried out. The drop weight, standard drop distance, probe cross-sectional area, and penetration depth were recorded during the test. Based on the mechanical mechanism of dynamic penetration tests, the single impact input energy was calculated according to the principle of work done by the drop weight. The single impact input energy is the product of the drop weight and the standard drop distance. The drop weight is obtained by multiplying the drop weight by the gravitational acceleration. The total impact input energy was obtained by accumulating the impact input energy of multiple impacts. The total input energy was divided by the product of the probe cross-sectional area and the actual penetration depth to obtain the impact input energy per unit penetration depth and per unit probe cross-sectional area, which is used as the energy index characterizing the soil stiffness of the corresponding grid cell.

2. The landslide slip zone identification method based on three-dimensional coordinate calculation according to claim 1, characterized in that: The specific method for performing joint differential calculations on all GNSS monitoring points and fixing the carrier phase ambiguity to obtain the time-series three-dimensional coordinates of each GNSS monitoring point is as follows: Simultaneously extract carrier phase observations and pseudorange observations from each GNSS monitoring point and the reference point. Select one satellite as the reference satellite from all satellites observed by both, and use the remaining satellites as non-reference satellites. The difference between the carrier phase observations of the same satellite by the GNSS monitoring point and the reference point is calculated to obtain the inter-station single difference value of the carrier phase of the corresponding satellite. The pseudorange observations of the same satellite by the GNSS monitoring point and the reference point are calculated to obtain the single difference value between the pseudoranges of the corresponding satellite between stations. The difference between the carrier phase inter-station single difference value of each non-reference satellite and the carrier phase inter-station single difference value of the reference satellite is calculated to obtain the carrier phase inter-station-inter-satellite double difference value. The pseudo-range station-to-station single difference value of each non-reference satellite is calculated by comparing it with the pseudo-range station-to-station single difference value of the reference satellite to obtain the pseudo-range station-to-satellite double difference value. By using the pseudo-range inter-station-inter-satellite double difference and the reference coordinates of the benchmark point, the pseudo-range inter-station-inter-satellite double difference observation equation is constructed, and the approximate three-dimensional coordinates of the GNSS monitoring point are obtained by solving it. Using the approximate three-dimensional coordinates as constraints, and combining the carrier phase inter-station-inter-satellite double difference, a carrier phase inter-station-inter-satellite double difference observation equation is constructed. The integer ambiguity in the equation is searched, matched and fixed by the LAMBDA algorithm. Substitute the fixed integer ambiguity back into the carrier phase inter-station-inter-satellite double-difference observation equation to calculate the three-dimensional coordinate increment of the GNSS monitoring point relative to the reference point. By combining the known reference coordinates of the benchmark point with the three-dimensional coordinate increments, the three-dimensional coordinates of the GNSS monitoring points at each monitoring time are solved, and the coordinates are arranged in the order of monitoring time to generate a time sequence of the three-dimensional coordinates of each GNSS monitoring point.

3. The landslide slip zone identification method based on three-dimensional coordinate calculation according to claim 2, characterized in that: The method for calculating the three-dimensional displacement vector, displacement rate, and displacement gradient of each GNSS monitoring point based on the coordinate difference between adjacent monitoring times is as follows: Three-dimensional coordinate time series data of two adjacent monitoring times of the same GNSS monitoring point are selected. The coordinates of the previous monitoring time are used as the initial reference and the coordinates of the next monitoring time are used as the deformed coordinates. The displacement change components of the GNSS monitoring point in three-dimensional space are obtained by coordinate difference calculation. The displacement components of each dimension are integrated to construct a three-dimensional displacement vector. Based on the time interval between adjacent monitoring times and the magnitude of the corresponding three-dimensional displacement vector, the displacement change per unit time is calculated to obtain the displacement rate of the GNSS monitoring point. Taking a single GNSS monitoring point as the center, and combining the three-dimensional displacement data of the surrounding adjacent GNSS monitoring points, the displacement spatial change rate at the location of the GNSS monitoring point is solved to obtain the corresponding displacement gradient; Traverse all GNSS monitoring points within the region, and solve for the three-dimensional displacement vector, displacement rate, and displacement gradient of each GNSS monitoring point according to the above calculation logic.

4. The landslide slip zone identification method based on three-dimensional coordinate calculation according to claim 1, characterized in that: The method for dividing the landslide area into grid cells based on the spatial distribution of each GNSS monitoring point is as follows: Based on the planar coordinates of each GNSS monitoring point in the local coordinate system of the slope, the spatial distribution range of each GNSS monitoring point is determined. Using the planar position of each GNSS monitoring point as a node, the Delaunay triangulation algorithm is used to construct multiple irregular triangular grid cells covering the landslide area.

5. The landslide slip zone identification method based on three-dimensional coordinate calculation according to claim 1, characterized in that: The energy index characterizing the stiffness of soil and rock is obtained for each grid cell. For grid cells with energy indices higher than a preset stiffness threshold, a screening rule is constructed based on the surface crack development parameters of each grid cell to screen unstable grid cells, as follows: Preset threshold values ​​for surface fissure development parameters, including threshold values ​​for average fissure width, fissure extension length, and fissure development density. If all surface fissure development parameters of a grid cell exceed the corresponding surface fissure development parameter threshold, it is determined to be an unstable grid cell. Traverse all grid cells in the landslide area and perform cell-by-cell discrimination according to the above screening rules to identify unstable grid cells.

6. The landslide slip zone identification method based on three-dimensional coordinate calculation according to claim 1, characterized in that: The specific steps for calculating three types of difference characteristics among neighboring grid cells are as follows: Displacement difference amplitude, deformation rate deviation, and deformation gradient change amplitude. Traverse all neighboring mesh elements and calculate the mean of the three-dimensional displacement vector, displacement rate, and displacement gradient for each mesh element as the corresponding three-dimensional displacement vector, displacement rate, and displacement gradient. Extract the three-dimensional displacement vector, displacement rate, and displacement gradient for each pair of adjacent mesh elements. Calculate the absolute difference between the magnitudes of the two sets of three-dimensional displacement vectors as the displacement difference amplitude between the two sets of mesh elements. Calculate the vector deviation between the two sets of displacement rates as the deformation rate deviation between the two sets of mesh elements. Calculate the difference between the magnitudes of the two sets of displacement gradients as the deformation gradient change amplitude between the two sets of mesh elements.

7. The landslide slip zone identification method based on three-dimensional coordinate calculation according to claim 6, characterized in that: The method for constructing the discrimination rules for identifying sliding bands is as follows: For each unstable grid cell, a profile line parallel to the main slip direction is drawn through its center point. All grid cells along the profile line are extracted and arranged in order of elevation from high to low to form a process analysis grid sequence. Using the unstable grid cells within the analysis grid sequence along the path as boundary nodes, a unified elevation datum is selected, and the elevation difference of each boundary node relative to the unified elevation datum is calculated. Simultaneously calculate three types of difference features between each boundary node and its adjacent downstream grid cells: displacement difference amplitude, deformation rate deviation, and deformation gradient change amplitude. Construct a training dataset using the elevation difference, impact input energy, and corresponding three types of features obtained from all boundary nodes based on historical data. Using elevation difference and impact input energy in the training dataset as input features, and three types of difference features—displacement difference amplitude, deformation rate deviation, and deformation gradient change amplitude—as output labels, a regression algorithm is used to train the corresponding feature threshold prediction models for the three types of features. The elevation difference and impact input energy of the boundary node are input into the trained feature threshold prediction model to output the displacement difference threshold, velocity deviation threshold and gradient mutation threshold corresponding to the boundary node and its adjacent downstream grid cells.

8. The landslide slip zone identification method based on three-dimensional coordinate calculation according to claim 7, characterized in that: The method for identifying the sliding zone boundary and marking the corresponding mesh element is as follows: A threshold for the rock and soil stiffness energy index is preset. Grid cells with an energy index higher than the threshold are designated as high-energy grid cells. Each high-energy grid cell is used to construct a path analysis grid sequence along the main sliding direction of the landslide. In the analysis grid sequence constructed by any high-energy grid cell, the elevation difference and impact input energy of the high-energy grid cell are used as inputs to the trained feature threshold prediction model to output three feature thresholds: displacement difference threshold, velocity deviation threshold, and gradient abrupt change threshold corresponding to the high-energy grid cell and its adjacent downstream grid cells. By monitoring the displacement difference amplitude, deformation rate deviation, and deformation gradient change amplitude of the high-energy grid cell and its downstream grid cell in real time, if the high-energy grid cell simultaneously meets the three deformation parameter exceeding threshold conditions, the grid cell is determined to be a deformation anomalous grid. In the analysis grid sequence along the path of all deformation anomalous grid cells, the grid cells in the opposite direction to the main sliding direction are all determined to be the area where the sliding zone is located. Based on the outer contour of the grid cells in the area where the sliding zone is located, a closed spatial boundary is generated, and the distribution range of at least one landslide sliding zone is delineated, thus completing the delineation of the spatial boundary of the sliding zone. All grid cells belonging to the landslide sliding zone are uniformly marked.

9. A landslide slip zone identification system based on three-dimensional coordinate calculation, characterized in that: The system is used to execute the landslide slip zone identification method based on three-dimensional coordinate calculation as described in any one of claims 1-8: The reference module is used to set up multiple GNSS monitoring points in the landslide area and at least one reference point outside the landslide area to collect multi-mode GNSS raw observation data from each GNSS monitoring point. The multi-mode GNSS raw observation data includes carrier phase observations and pseudorange observations from no less than two satellite navigation systems. The data processing module is used to combine the known fixed reference coordinates of the benchmark point with all multi-mode GNSS raw observation data to perform joint differential calculation on all GNSS monitoring points, and to fix the carrier phase ambiguity, so as to obtain the time-series three-dimensional coordinates of each GNSS monitoring point. The feature calculation module is used to uniformly transform the three-dimensional coordinate time series of each GNSS monitoring point to the local coordinate system of the slope. Based on the coordinate difference between adjacent monitoring times, it solves the three-dimensional displacement vector, displacement rate and displacement gradient of each GNSS monitoring point respectively. The grid selection module is used to divide the landslide area into grid units according to the spatial distribution of each GNSS monitoring point, determine the main sliding direction of the landslide by statistically analyzing the average direction of the three-dimensional displacement vector of each GNSS monitoring point, obtain the energy index to characterize the rock and soil stiffness of the corresponding grid unit through dynamic penetration test at the center point of each grid unit, and select unstable grid units by combining the surface crack development parameters of each grid unit. The difference identification module is used to calculate three types of difference features between the unstable grid cell and its adjacent downstream grid cell based on three basic deformation features along the main sliding direction of the landslide: displacement difference amplitude, deformation rate deviation, and deformation gradient change amplitude. This is used to construct a discrimination rule for identifying the sliding zone. Based on this discrimination rule, the module analyzes the three types of difference features of the neighboring grid cells, identifies the boundary of the sliding zone, and marks the corresponding grid cells.

Images