A mine slope absolute displacement field monitoring and error self-correction method and system

Abstract

This invention relates to the field of mine slope deformation monitoring technology, specifically to a method and system for monitoring and self-correcting the absolute displacement field of mine slopes. The method includes the following steps: processing the atmospheric-corrected GB-SAR relative displacement, projection coefficient matrix, and motion pattern partitioning based on Kriging synergetics interpolation and adaptive weighted fusion to obtain the GNSS-calibrated absolute displacement field; processing the displacement period characteristic parameters and the surface displacement gradient field using multi-feature fuzzy comprehensive evaluation to obtain the comprehensive index field of local slope stability; and combining the GNSS warning level, the GNSS-calibrated absolute displacement field, and the comprehensive local stability index, using graded warning and displacement field source analysis to complete high-precision absolute displacement field monitoring and error self-correction. This invention enables high-precision monitoring of surface absolute displacement fields, accurate identification of the causes of displacement anomalies, and error self-correction.

Description

Technical Field

[0001] This invention relates to the field of mine slope deformation monitoring technology, specifically to a method and system for monitoring the absolute displacement field and self-correcting errors of mine slopes. Background Technology

[0002] Mine slope deformation monitoring is a core component of ensuring safe mine production. Its core requirements are high-precision, full-area coverage of displacement field acquisition and reliable landslide hazard early warning. Among existing monitoring technologies, Global Navigation Satellite Systems (GNSS) offer high absolute displacement measurement accuracy, but can only achieve point-based monitoring, resulting in incomplete point coverage and an inability to reflect the planar deformation characteristics of slopes. While ground-based synthetic aperture radar (GB-SAR) can achieve planar coverage monitoring and capture the overall deformation trend of slopes, it is susceptible to atmospheric refraction and signal interference, its measurement accuracy relies on external correction, and it can only acquire line-of-sight (LOS) displacement, failing to directly provide three-dimensional displacement information, leading to inaccurate planar displacement.

[0003] To address the limitations of single sensors, existing technologies attempt to fuse GNSS and GB-SAR data, but several key technical shortcomings remain: First, atmospheric correction often employs fixed atmospheric models or single-moment correction methods, failing to consider the spatiotemporal dynamic heterogeneity of the atmosphere caused by dust and large diurnal temperature variations in mining environments, easily resulting in residual atmospheric errors. Second, when converting line-of-sight displacement to three-dimensional displacement using GB-SAR, a general imaging geometric model is often used, ignoring the uneven terrain features and structural surfaces of mine slopes, leading to low conversion accuracy. Third, the data fusion process uses a globally unified weighting strategy, failing to consider the differences in motion patterns across different areas of the slope, thus failing to achieve precise complementarity of data advantages. Fourth, early warning methods often use single displacement or rate thresholds, failing to deduct periodic fluctuations such as precipitation and freeze-thaw cycles, and failing to set dynamic thresholds based on geological conditions, easily leading to false alarms and missed alarms. Fifth, the lack of a complete error tracing and self-correction closed-loop mechanism makes it impossible to accurately identify the causes of displacement anomalies (sensor failure, mining interference, environmental interference, or actual deformation), hindering dynamic optimization of monitoring accuracy. Therefore, given the complex geological and production scenarios of mine slopes, there is an urgent need for a fusion monitoring method that can achieve high-precision surface absolute displacement field monitoring, accurately identify the causes of displacement anomalies, and has error self-correction capabilities, in order to solve the core problems of incomplete coverage, insufficient accuracy, unreliable early warning, and lack of error correction mechanisms in existing technologies. Summary of the Invention

[0004] To address the shortcomings of existing methods and the needs of practical applications, and to solve the aforementioned problems, this invention provides a method for monitoring the absolute displacement field and self-correcting errors of mine slopes, comprising the following steps:

[0005] Based on GNSS anchor point-based atmospheric phase screen iterative modeling and correction, GB-SAR time-series differential phase map data is processed to obtain atmospherically corrected GB-SAR relative displacement sequences. Through geometric constraint calibration of line-of-sight to three-dimensional displacement vector transformation, GNSS three-dimensional coordinate time-series data and GB-SAR pixel geometric information are analyzed to obtain the projection coefficient matrix from GB-SAR LOS displacement to true three-dimensional displacement. Based on Kriging co-interpolation and adaptive weighted fusion, the atmospherically corrected GB-SAR relative displacement, projection coefficient matrix, and motion mode partitioning are processed to obtain the GNSS-calibrated absolute displacement field. Using multi-feature fuzzy comprehensive evaluation, displacement period characteristic parameters and surface displacement gradient fields are processed to obtain a comprehensive index field of slope local stability. Combining GNSS warning levels, the GNSS-calibrated absolute displacement field, and the comprehensive local stability index, high-precision absolute displacement field monitoring and error self-correction are achieved using graded warning and displacement field source analysis.

[0006] Optionally, the atmospheric phase screen iterative modeling and correction based on GNSS anchor points, which processes GB-SAR time-series differential phase map data to obtain atmospherically corrected GB-SAR relative displacement sequences, includes the following steps:

[0007] GNSS anchor points are selected, and atmospheric delay at the anchor points is estimated. An initial atmospheric phase screen is constructed, and the atmospheric delay and the initial atmospheric phase screen are combined to obtain the atmospherically corrected GB-SAR relative displacement sequence.

[0008] Optionally, the geometric constraint calibration through line-of-sight to three-dimensional displacement vector transformation, analyzing GNSS three-dimensional coordinate time series data and GB-SAR pixel geometric information to obtain the projection coefficient matrix from GB-SAR LOS displacement to true three-dimensional displacement, includes the following steps:

[0009] Extract the three-dimensional displacement of GNSS monitoring points and the local geometric features of the slope; construct an initial projection coefficient model, and use the three-dimensional displacement and the local geometric features of the slope to perform geometric constraint calibration and coefficient optimization on the initial projection coefficient model to obtain the projection coefficient matrix from GB-SAR LOS displacement to the true three-dimensional displacement.

[0010] Optionally, the step of processing the atmospheric-corrected GB-SAR relative displacement, projection coefficient matrix, and motion mode partitions based on Kriging co-interpolation and adaptive weighted fusion to obtain the GNSS-calibrated absolute displacement field includes the following steps:

[0011] The relative displacement of GB-SAR after atmospheric correction is transformed into three dimensions using the projection coefficient matrix; based on the transformation results and motion mode partitioning, Kriging co-interpolation fusion is performed according to the partition variogram function and adaptive weights to obtain the absolute displacement field of GNSS calibration.

[0012] Optionally, extracting the motion pattern partition includes the following steps:

[0013] Phase statistical features are extracted and the areas affected by mining are marked. Motion pattern clustering analysis is performed by combining the phase statistical features and the marking results to generate a motion pattern partition map.

[0014] Optionally, the step of using multi-feature fuzzy comprehensive evaluation to process the displacement period feature parameters and the surface displacement gradient field to obtain the comprehensive index field of local slope stability includes the following steps:

[0015] An evaluation index system is constructed using displacement period characteristic parameters and surface displacement gradient field; index weights are determined by combining the analytic hierarchy process (AHP) and mine geological constraints; a fuzzy evaluation matrix is ​​constructed based on the evaluation index system and index weights, and a comprehensive index field for local slope stability is generated through fuzzy comprehensive evaluation.

[0016] Optionally, extracting the displacement period characteristic parameters includes the following steps:

[0017] Gross error identification and multi-criteria robust filtering are used to obtain denoised GNSS coordinate time series data; displacement trend term and period term are separated from the GNSS coordinate time series data, and then displacement period feature parameters are extracted.

[0018] Optionally, analyzing the surface displacement gradient field includes the following steps:

[0019] Coherent point targets are screened by region to obtain a preliminary coherent point target set; a spatial adjacency network of coherent point targets is constructed based on the preliminary coherent point target set, and then the relative displacement difference between adjacent coherent point targets is quantified; the surface deformation gradient field is analyzed through the relative displacement of the neighborhood.

[0020] Optionally, the process of combining GNSS early warning levels, GNSS-calibrated absolute displacement fields, and local stability comprehensive indicators, and utilizing graded early warning and displacement field source analysis to achieve high-precision absolute displacement field monitoring and error self-correction, includes the following steps:

[0021] By combining GNSS early warning levels, absolute displacement field deformation, and local stability comprehensive indicators, graded early warning of mine slopes can be achieved; through displacement field source analysis, displacement anomalies and error sources can be identified, and error self-correction can be achieved.

[0022] This invention lays the foundation for high-precision relative displacement by iteratively correcting atmospheric phase screens to remove complex atmospheric interference; it achieves a breakthrough in GB-SAR dimensions through geometric constraint calibration, providing the core basis for three-dimensional displacement conversion; it leverages the complementary advantages of point and surface data through Krifin integration to obtain a high-precision absolute displacement field across the entire region; it quantifies the stability surface domain based on multi-feature fuzzy evaluation to accurately locate potential hazard areas; and it combines hierarchical early warning and source tracing analysis to complete monitoring, early warning, and error self-correction, ensuring data reliability and landslide early warning accuracy, and providing scientific support for mine slope safety management.

[0023] Secondly, to efficiently execute the mine slope absolute displacement field monitoring and error self-correction method provided by this invention, this invention also provides a mine slope absolute displacement field monitoring and error self-correction system, including a processor, an input device, an output device, and a memory. The processor, input device, output device, and memory are interconnected. The memory stores a computer program containing program instructions. The processor is configured to call the program instructions to execute the mine slope absolute displacement field monitoring and error self-correction method as described in the first aspect of this invention. The mine slope absolute displacement field monitoring and error self-correction system of this invention has a compact structure and stable performance, and can stably execute the mine slope absolute displacement field monitoring and error self-correction method provided by this invention, further improving the overall applicability and practical application capability of this invention. Attached Figure Description

[0024] Figure 1 This is a flowchart of a method for monitoring the absolute displacement field and self-correcting errors of a mine slope, provided in an embodiment of the present invention.

[0025] Figure 2 This is a framework diagram of a mine slope absolute displacement field monitoring and error self-correction system provided in an embodiment of the present invention. Detailed Implementation

[0026] Specific embodiments of the present invention will now be described in detail. It should be noted that the embodiments described herein are for illustrative purposes only and are not intended to limit the invention. In the following description, numerous specific details are set forth in order to provide a thorough understanding of the invention. However, it will be apparent to those skilled in the art that these specific details are not necessary to practice the invention. In other instances, well-known circuits, software, or methods have not been specifically described to avoid obscuring the invention.

[0027] Throughout this specification, references to an embodiment, example, or illustration mean that a particular feature, structure, or characteristic described in connection with that embodiment or example is included in at least one embodiment of the invention. Therefore, phrases appearing in various places throughout the specification, such as "in one embodiment," "in an embodiment," "an example," or "an illustration," do not necessarily refer to the same embodiment or example. Furthermore, specific features, structures, or characteristics can be combined in any suitable combination and / or sub-combination in one or more embodiments or examples. Moreover, those skilled in the art will understand that the illustrations provided herein are for illustrative purposes and are not necessarily drawn to scale.

[0028] Please see Figure 1 To address the aforementioned problems, this invention provides a method for monitoring the absolute displacement field and self-correcting errors of mine slopes, such as... Figure 1 As shown, in one embodiment, the method includes the following steps:

[0029] S1. Iterative modeling and correction of atmospheric phase screen based on GNSS anchor points: Process GB-SAR time-series differential phase map data to obtain atmospherically corrected GB-SAR relative displacement sequence.

[0030] The atmospheric phase screen iterative modeling and correction based on GNSS anchor points processes GB-SAR time-series differential phase map data to obtain atmospherically corrected GB-SAR relative displacement sequences, including the following steps:

[0031] S11. Select GNSS anchor points and estimate the atmospheric delay at the anchor points.

[0032] The stability of the anchor points directly determines the reliability of the atmospheric correction benchmark. Therefore, it is necessary to strictly follow the three principles of being far away from mining-related influences, having uniform geological conditions, and ensuring unobstructed signals. 3-5 candidate anchor points should be selected in the non-mining-related influence zone more than 500m outside the mine slope. The candidate points should avoid areas where large equipment is operating, under power transmission lines, and areas of loose deposits.

[0033] To cover the atmospheric fluctuation characteristics of the entire mining production cycle, 30 consecutive days of raw GNSS data were selected for stability analysis. First, the 3σ criterion was used to eliminate outliers in the raw data. Then, the three-dimensional coordinate time series of each candidate point was calculated, and the displacement standard deviation of each coordinate component (x, y, z directions) was calculated. Furthermore, a multi-dimensional threshold method was used for stability determination: the displacement standard deviation of a single coordinate component < 2 mm, and the three-dimensional composite displacement standard deviation < 3 mm, while the daily average displacement change between two adjacent days < 0.5 mm. Only candidate points meeting all these conditions were included in the stable anchor point set.

[0034] Then, the carrier phase data of the GNSS anchor points are preprocessed, including satellite orbit error correction, clock error correction, and multipath effect suppression (using the Melbourne-Wübbena combination method to weaken the multipath effect). Next, atmospheric delay components are separated: the ionospheric delay is calculated using a dual-frequency combination method (L1 / L2 carrier phase without ionospheric combination); the tropospheric delay uses the Saastamoinen model as the base model, combined with localized corrections based on time-series meteorological data (temperature, pressure, and humidity observations) of the mining area. Considering the high dust levels in the mines, a visibility parameter is introduced to adjust the wet delay coefficient in the model, reducing the error in delay estimation caused by refractive index changes due to dust. Finally, the separated ionospheric delay is superimposed with the corrected tropospheric delay to obtain the total atmospheric delay value at each anchor point, forming an atmospheric delay time series.

[0035] S12. Construct an initial atmospheric phase screen, and combine the atmospheric delay and the initial atmospheric phase screen to obtain the atmospherically corrected GB-SAR relative displacement sequence.

[0036] First, each GB-SAR time-series differential phase map is preprocessed to remove pixels that fail to unwrap phase (pixels with phase values ​​exceeding the reasonable range of [-π, π]), and median filtering is used to smooth phase noise. Then, the pixel phase values ​​corresponding to the GNSS anchor point positions in each GB-SAR differential phase map are extracted. Considering the spatial matching error between the anchor point and GB-SAR pixels in the mining scene, the average phase value of the 3×3 pixels surrounding the anchor point projection position is selected as the GB-SAR phase value corresponding to the anchor point to reduce the impact of single-point pixel noise.

[0037] Next, based on the total atmospheric delay value at the anchor point and the corresponding GB-SAR phase value, a linear mapping relationship between the two is established. The mapping coefficient is obtained by fitting using the least squares method to quantify the degree of atmospheric delay interference with the phase. Considering the significant spatial heterogeneity of the mine atmosphere (e.g., the temperature difference between the top and bottom of the slope can reach 5-8℃, and the dust distribution is uneven), an inverse distance weighted interpolation method is used to generate the initial atmospheric phase screen. During interpolation, the weight index is adjusted according to the complexity of the mine terrain (the weight index is 2.5 for areas with large terrain undulations and 2.0 for flat areas) to ensure that the interpolation results can accurately reflect the spatial gradation characteristics of the atmospheric phase. The final generated initial atmospheric phase screen must be completely matched with the GB-SAR monitoring area, and the pixel resolution must be consistent with the original GB-SAR data.

[0038] Furthermore, the initial atmospheric phase screen is stripped frame by frame from the corresponding GB-SAR time-series differential phase map to obtain preliminary atmospheric phase data. Then, the relative displacement sequence at each anchor point is inverted based on the GB-SAR imaging geometric parameters (azimuth, elevation, and range resolution).

[0039] Simultaneously, atmospheric influence was subtracted from the measured displacement data of the GNSS anchor points (the total atmospheric delay value was converted into the displacement influence), yielding the true relative displacement sequence at the anchor points. Residual analysis was performed on the relative displacement sequences obtained by the two methods, calculating the mean and standard deviation of the residuals at each time step. If the mean residual > 0.5 mm or the standard deviation > 0.8 mm, it indicates that the initial atmospheric phase screen did not fully capture the spatiotemporal characteristics of atmospheric changes. The model parameters need to be corrected to reconstruct the phase screen, adjust the weighting index of the inverse distance weighted interpolation, or supplement atmospheric delay data from adjacent time steps to optimize the mapping relationship. The initial atmospheric phase screen construction and phase stripping process should then be repeated.

[0040] The residual index is recalculated after each iteration until the iteration converges. The atmospheric phase data obtained after convergence is then converted into atmospherically corrected GB-SAR relative displacement sequence using the phase-displacement conversion formula (combined with GB-SAR operating frequency and observation distance to calculate the conversion coefficient).

[0041] This invention uses GNSS anchor points in stable mining areas as dynamic atmospheric benchmarks, achieves adaptive correction of atmospheric phase through iterative modeling, and incorporates optimized model parameters based on mining meteorological data. This solves the problem of GB-SAR phase interference caused by atmospheric refraction in complex mining atmospheric environments, providing high-precision relative displacement baseline data for subsequent absolute displacement field inversion and avoiding deformation misjudgment caused by atmospheric errors.

[0042] S2. By calibrating the geometric constraints of the line-of-sight to three-dimensional displacement vector transformation, analyze the GNSS three-dimensional coordinate time series data and GB-SAR pixel geometric information to obtain the projection coefficient matrix from GB-SAR LOS displacement to the real three-dimensional displacement.

[0043] GB-SAR can only acquire line-of-sight (LOS) displacement and cannot reflect the three-dimensional deformation characteristics of the slope, while GNSS can directly acquire three-dimensional absolute displacement. This step utilizes the digital elevation model (DEM) of the mine slope and the geometric parameters of GB-SAR imaging to construct the conversion relationship between LOS displacement and three-dimensional displacement. The projection coefficient matrix of the conversion model is calibrated using the three-dimensional displacement data of GNSS monitoring points. At the same time, the constraints are optimized by combining the actual geometric characteristics of the mine slope, such as slope and aspect, to improve the conversion accuracy.

[0044] Specifically, the geometric constraint calibration through line-of-sight to three-dimensional displacement vector transformation, analyzing GNSS three-dimensional coordinate time series data and GB-SAR pixel geometric information to obtain the projection coefficient matrix from GB-SAR LOS displacement to true three-dimensional displacement includes the following steps:

[0045] S21. Extract the three-dimensional displacement of GNSS monitoring points and the local geometric features of the slope.

[0046] Precise calculations were performed using professional GNSS data processing software. A local coordinate system was established based on a stable anchor point set. A relative positioning mode was used for the slope monitoring points (the baseline calculation adopted the L1 / L2 dual-frequency carrier phase combination method to reduce the influence of the ionosphere and multipath). The calculation yielded the three-dimensional coordinate time series of each monitoring point once per hour. Then, the three-dimensional displacement time series (Δx, Δy, Δz) of each point was obtained by calculating the coordinate difference between adjacent time points, where Δx is the east-west displacement, Δy is the north-south displacement, and Δz is the vertical displacement.

[0047] Furthermore, a 3×3 window analysis method was used to extract the local geometric parameters corresponding to each GNSS monitoring point and GB-SAR pixel: the slope α was calculated using the maximum slope drop method, which is the angle corresponding to the maximum elevation difference between the central pixel and the eight adjacent pixels within the window, reflecting the local steepness of the slope; the slope aspect β was calculated using the azimuth method, which is the angle between the direction of the maximum slope drop and the due north direction (range 0°-360°), reflecting the orientation of the slope.

[0048] To adapt to the characteristics of the development of the structural surface of the mine slope, the topographic curvature (reflecting the degree of convexity and concavity of the slope) at each location was also extracted. Combined with the structural surface orientation data in the geological survey report, the dominant direction of slope deformation was identified. For example, when the slope aspect is consistent with the orientation of the structural surface, shear deformation along the slope aspect is likely to occur; when the slope has obvious concave curvature, the risk of tensile deformation perpendicular to the slope aspect is higher, providing key geometric basis for the subsequent constraint optimization of the projection coefficient matrix.

[0049] S22. Construct an initial projection coefficient model, and use the three-dimensional displacement and the local geometric features of the slope to perform geometric constraint calibration and coefficient optimization on the initial projection coefficient model to obtain the projection coefficient matrix from GB-SAR LOS displacement to the real three-dimensional displacement.

[0050] Based on the imaging geometry principle of GB-SAR, the spatial relationship between the LOS direction and the three-dimensional coordinate system (East-North-Sky coordinate system) is clarified: the LOS direction is the straight line direction from the GB-SAR antenna to the monitored pixel, and its spatial attitude is determined by the GB-SAR system parameters (azimuth angle α_sar, elevation angle β_sar).

[0051] Establish a linear transformation equation between the LOS displacement ΔLOS and the three-dimensional displacement (Δx, Δy, Δz): ΔLOS = A·[Δx, Δy, Δz] TA is a 3×1 projection coefficient matrix used to quantify the projection weights of each component of the three-dimensional displacement in the LOS direction. The calculation of the initial projection coefficient matrix A requires combining GB-SAR imaging parameters and the local slope aspect β: First, calculate the angle θ between the LOS direction and the vertical direction (skyward direction) (calculated from the GB-SAR elevation angle β_sar, observation distance D, and antenna height H, θ=arctan(D / H)-β_sar); then correct the horizontal projection weights according to the slope aspect β, and finally obtain the initial projection coefficient matrix A=[cosθ·sinβ,cosθ·cosβ,sinθ].

[0052] To address the issue of localized occlusion on mine slopes, occlusion detection (based on line-of-sight analysis of DEM) is performed on GB-SAR pixels. Occluded pixels are marked and their initial projection coefficients are temporarily set to 0. These pixels are then corrected during subsequent optimization to avoid conversion errors caused by occlusion.

[0053] Furthermore, the three-dimensional displacement time series (Δx, Δy, Δz) of each GNSS monitoring point is substituted into the initial transformation equation time by time to calculate the theoretical LOS displacement ΔLOS_sim at each time. At the same time, the measured LOS displacement ΔLOS_obs of the corresponding pixels of each GNSS monitoring point is extracted from the atmospherically corrected GB-SAR relative displacement sequence. Then, a geometrically constrained objective function is constructed: with minimizing the sum of squared residuals between theoretical and measured LOS displacements as the core objective, combined with the local slope constraint of the slope: for areas with slope α > 30° (steep slope areas in mines, where the vertical slope displacement accounts for a higher proportion), the weight coefficient of the vertical slope displacement component (Δx·sinβ+Δy·cosβ) is increased (weight coefficient = 1+0.02α, where α is in degrees), to enhance the deformation adaptability of steep slope areas.

[0054] The optimized projection coefficient matrix A is solved by iterative weighted least squares method: In each iteration, observations with residual absolute value > 1 mm are given a lower weight (weight coefficient = 1 / residual absolute value) to weaken the influence of outliers; the number of iterations is set to 10, and the iteration can be terminated in advance if the residual sum of squares converges during the iteration.

[0055] Furthermore, the accuracy of the projection coefficient matrix is ​​verified using root mean square error (RMSE) and residual distribution checks: if RMSE > 0.3 mm or the residuals are non-normally distributed, structural surface constraints are added. Combining the structural surface orientation γ from the geological survey, a penalty term for the structural surface directional displacement component (Δx·sinγ+Δy·cosγ) is introduced to correct the horizontal weight distribution in the projection coefficient matrix. The above optimization process is repeated until RMSE ≤ 0.3 mm and the residuals are normally distributed, and the final projection coefficient matrix is ​​output.

[0056] This invention uses the local slope and aspect extracted from the DEM of the mine slope as geometric constraints in the coefficient calibration process, and dynamically optimizes the coefficient matrix using GNSS measured data to achieve accurate conversion of GB-SAR LOS displacement to three-dimensional displacement.

[0057] S3. Based on Kriging co-interpolation and adaptive weighted fusion, the atmospheric-corrected GB-SAR relative displacement, projection coefficient matrix, and motion mode partition are processed to obtain the absolute displacement field of GNSS calibration.

[0058] The deformation and movement patterns of different areas of the mine slope vary significantly (e.g., stable zone, uniform deformation zone, accelerated deformation zone, and abrupt deformation zone). By statistically analyzing the atmospheric-corrected GB-SAR time series phase, multi-dimensional feature parameters are extracted. Combined with the mine's mining plan and the distribution of goaf areas, a clustering algorithm is used to divide the slope movement pattern into zones, providing a basis for zone adaptation for subsequent data fusion.

[0059] Specifically, extracting the motion pattern partition includes the following steps:

[0060] S311. Extract phase statistical features and mark the area affected by mining.

[0061] First, the temporal phase data of each coherent point target is extracted from the atmospherically corrected GB-SAR relative displacement sequence, and then the phase statistical characteristic parameters are analyzed, including:

[0062] Displacement rate: Calculated using a weighted linear fitting method, with the coherence coefficient of the phase data as the weight. Higher coherence coefficients result in higher weights, preventing low-coherence data from lowering the fitting accuracy. The fitting formula is as follows: ,in Let be the phase value at time t. To fit the displacement rate, This is the initial phase;

[0063] Rate stability: The standard deviation σ_v of the fitted rate sequence is calculated. σ_v < 0.2 mm / d is considered to be stable rate, and σ_v > 0.5 mm / d is considered to be drastic rate fluctuation. This parameter can reflect the stability of deformation in the area where the coherent point target is located and is adapted to the rate fluctuation characteristics under mining disturbance.

[0064] Phase residuals: Calculate the root mean square (RMSE) of the linear fitting residuals. Before calculating the residuals, the residual sequence needs to be smoothed using the 5-point moving average method to eliminate random noise. RMSE < 0.1 rad indicates that the phase fitting effect is good, while RMSE > 0.3 rad requires a re-examination of the phase data quality.

[0065] Number of mutations: A dynamic threshold is set based on the displacement standard deviation of coherent point targets in the stable region. The calibrated coherent point targets in the stable region are selected, and the 3σ value of their displacement sequence is calculated as the mutation threshold. σ is the displacement standard deviation of coherent point targets in the stable region. The number of times the displacement exceeds the threshold in the time series is counted. Each mutation must meet the confirmation condition of exceeding the threshold for two consecutive monitoring times to avoid misstatistics caused by single-point data anomalies.

[0066] After all feature parameters are extracted, a 4-dimensional feature vector [displacement rate, rate stability, phase residual, number of mutations] is finally formed for each coherent point target.

[0067] Furthermore, an in-depth analysis of the mining plan was conducted, reviewing the mining plan documents for the next 1-3 years and extracting core parameters, including the advancing direction, rate, mining strata, and mining intensity of the mining face. Based on these parameters, empirical formulas were used to calculate the theoretical boundary of the mining impact area: the horizontal boundary of the mining impact area is 50-100m ahead of the mining face in the advancing direction (the greater the mining intensity, the farther the boundary); the lateral boundaries are 30-50m beyond the width of the working face; the vertical boundary is from the surface above the mining strata to 20m below the mining strata, forming the preliminary theoretical range of the mining impact area.

[0068] Then, the labeling is optimized by combining the distribution data of the goaf: based on the goaf detection report, a three-dimensional model of the goaf is constructed to clarify the spatial coordinates, volume and burial depth of the goaf; for the surface area above the goaf, the boundary of the mining impact range is expanded because the goaf is prone to uneven surface subsidence and has higher deformation sensitivity; for the edge area of ​​the goaf, it is marked as the area with strong mining impact and is given the highest weight in subsequent clustering.

[0069] The third step is to verify historical deformation data: extract historical GB-SAR relative displacement data for the past year, analyze the deformation characteristics within the theoretical mining-affected area, and if the average displacement rate of coherent point targets in the area is >1 mm / d and the rate stability σ_v is >0.3 mm / d, which matches the characteristics of uniform or accelerated deformation that easily occurs in the mining-affected area, then the area is confirmed as a mining-affected area; if the displacement rate of coherent point targets in the area is <0.2 mm / d and there is no obvious rate fluctuation, then the theoretical boundary is adjusted and it is classified as a non-mining-affected area.

[0070] The final output is a mining-affected area marking map, which clearly divides the area into three categories: strong mining-affected area, weak mining-affected area, and non-mining-affected area. Among them, the non-mining-affected area is mostly a stable area, and its deformation characteristics (displacement rate < 0.2 mm / d, number of mutations = 0) can be used as the stable benchmark prior information for cluster analysis.

[0071] S312. Combine the phase statistical features and labeling results to perform motion pattern clustering analysis and generate a motion pattern partition map.

[0072] Due to the significant differences in the dimensions of the phase statistical feature parameters, the min-max normalization method is used to map each feature value to the [0,1] interval; then, the core parameters of weighted K-means clustering are set: using the normalized 4-dimensional feature vector as input, combined with the common deformation types of mine slopes, the number of clusters is set to 4, and the expected core features of each category are clarified:

[0073] ①Stable mode: Displacement rate < 0.2 mm / d, velocity stability σ_v < 0.2 mm / d, phase residual < 0.1 rad, number of abrupt changes = 0;

[0074] ② Uniform deformation mode: displacement rate 0.2-2 mm / d, velocity stability σ_v < 0.3 mm / d, phase residual < 0.2 rad, number of abrupt changes ≤ 1;

[0075] ③Accelerated deformation mode: displacement rate > 2 mm / d, velocity stability σ_v > 0.3 mm / d, phase residual 0.1-0.3 rad, number of abrupt changes 1-3;

[0076] ④ Abrupt deformation mode: displacement rate > 5 mm / d, velocity stability σ_v > 0.5 mm / d, phase residual > 0.3 rad, number of abrupt changes ≥ 3;

[0077] Furthermore, based on the marking results of mining-affected areas, sample weights were set: sample weight for areas with strong mining impact = 1.2, sample weight for areas with weak mining impact = 1.0, and sample weight for areas without mining impact = 0.8, thus strengthening the clustering priority of mining-dominated deformation areas.

[0078] The third step involves iterative cluster optimization: When initializing cluster centers, typical samples from each region are selected first to avoid initial center deviation leading to cluster failure; during the iteration process, Euclidean distance is used to calculate the similarity between samples and cluster centers, and the formula for updating cluster centers is as follows: , where C_k is the cluster center of the kth class, w_i is the weight of the i-th sample, and x_i is the feature vector of the i-th sample;

[0079] Furthermore, the silhouette coefficient is used to evaluate the clustering effect. The silhouette coefficient ranges from [-1, 1]. A value > 0.5 indicates that the clustering effect is good and the samples of each category are significantly distinguishable. If the silhouette coefficient < 0.3, the number of clusters or the sample weights are adjusted, and the clustering process is re-executed. Finally, the motion pattern category label of each coherent point target is output, forming a coherent point target-motion pattern association table.

[0080] Furthermore, a high-precision zoning map of mine slope motion patterns is generated based on seed point selection, constrained region growth, and post-zoning processing. First, seed points are refined: seed points are selected from the coherent point target-motion pattern association table. The selection criteria are: ① core samples belonging to the motion pattern (Euclidean distance between the sample and the corresponding cluster center < 0.2); ② coherent point target coherence coefficient ≥ 0.7, ensuring reliable seed point deformation characteristics; ③ uniform spatial distribution of seed points, avoiding insufficient growth due to sparse seed points. Then, region growth with terrain constraints is performed: using the selected coherent point target motion pattern category as the seed point, growth constraints are constructed based on the high-precision DEM data of the mine slope. In this embodiment, the slope difference threshold between adjacent pixels is set to 10°. If the slope difference between adjacent pixels is >10°, it is considered a terrain boundary and growth is stopped. During the growth process, the feature similarity between the pixel to be grown and the seed point is calculated. If the similarity is >0.8, the pixel is classified into the same motion pattern as the seed point. For blank areas without coherent point target coverage, the kriging interpolation method is used to obtain its feature parameters first, and then growth classification is performed. During interpolation, the variogram parameters are adjusted in combination with the mining-affected area marking results. The third step involves post-processing optimization of the partitions: ① Small area merging: Isolated partitions are removed and merged into adjacent main partitions; ② Boundary smoothing: Morphological closing operations are used to smooth partition boundaries, eliminating jagged edges and ensuring that partition boundaries match topographic features (such as slope shoulders and structural surfaces); ③ Cross-regional verification: Combining mining activity plans and goaf distribution data, the rationality of the partitions is verified. For example, areas heavily affected by mining should be mainly classified as uniform or accelerated deformation modes, while non-mining-affected areas should be mainly classified as stable modes. If contradictions exist, the growth constraints are locally adjusted (such as reducing the slope difference threshold to 8°) and the area is regrown. Finally, a motion mode partition map is generated. The map must include the partition range, boundary lines, typical coherence point target sample point locations, and motion mode labels for the four major motion modes, along with partition statistics (area, proportion, and average displacement rate of each mode area), providing a clear partition adaptation basis for subsequent data fusion (different fusion weights are used for different motion mode areas) and stability evaluation.

[0081] GB-SAR has the advantage of area coverage, but its accuracy is affected by motion patterns. GNSS has the advantage of absolute displacement accuracy, but its data is distributed in a point-like manner. Fusing the data from both can achieve a high-precision absolute displacement field combining point and area data. This step is based on motion pattern partitioning, using Kriging synergetics interpolation (considering the spatial correlation of the data), and combining it with synergetics theory (considering the synergistic effect of GB-SAR and GNSS data) to assign adaptive weights to regions with different motion patterns (e.g., higher weight for GNSS in accelerated deformation regions, and higher weight for GB-SAR in stable regions), thus achieving data fusion.

[0082] Specifically, the process of processing the atmospheric-corrected GB-SAR relative displacement, projection coefficient matrix, and motion mode partitions based on Kriging co-interpolation and adaptive weighted fusion to obtain the GNSS-calibrated absolute displacement field includes the following steps:

[0083] S321. Use the projection coefficient matrix to perform three-dimensional displacement transformation on the atmospherically corrected GB-SAR relative displacement.

[0084] Based on the geometric parameters of GB-SAR imaging and the projection coefficient matrix A, the inverse transformation equation from LOS displacement to three-dimensional displacement Δs_gb3d=A is constructed. -1 · ,in Δs_gb3d represents the atmospheric-corrected GB-SAR relative displacement, while Δs_gb3d represents the transformed three-dimensional displacement, which includes three components: east-west Δx, north-south Δy, and vertical Δz. Considering the dramatic topographic relief of mine slopes, a local slope α correction term is introduced during the transformation process. For steep slope areas with a slope α > 30°, a slope influence coefficient (coefficient = 1 + 0.01α) is added to the calculation of the vertical displacement component to compensate for the transformation deviation caused by the steep terrain.

[0085] S322. Based on the transformation results and motion mode partitioning, Kriging co-interpolation fusion is performed according to the partitioning mutation function and adaptive weights to obtain the absolute displacement field calibrated by GNSS.

[0086] As the core of Kriging interpolation, the variogram function's parameter settings need to be deeply adapted to the deformation characteristics of different motion modes of the mine slope. First, based on the motion mode zoning map, the variogram function type is accurately matched. In the stable zone (motion mode is stable), the deformation is gentle and the spatial correlation is strong, so a spherical variogram function is selected; in the accelerated deformation zone (motion mode is accelerated / abrupt deformation), the deformation difference is significant and the spatial correlation is weak, so an exponential variogram function is selected; in the uniform deformation zone, a spherical-exponential hybrid variogram function is used to take into account both its stable deformation and local small fluctuations. Second, the initial parameters are set: the initial parameters for the spherical variogram function in the stable zone are a range a = 20m (based on the spatial distribution density of coherent point targets in the stable zone, ensuring coverage of 3-5 adjacent coherent point targets), and sill value. Gold value The initial parameters for the exponential variogram of the accelerated deformation zone are: range a = 12m (the spatial correlation range of the deformation concentration zone is smaller) and sill value. Gold value 2Then, a surface displacement gradient field is introduced for parameter optimization: for regions with high gradient values ​​(areas with concentrated deformation), the range of the variogram in the corresponding region is reduced by 20%-30% to enhance the interpolation sensitivity of local deformation details; for regions with low gradient values ​​(areas with uniform deformation), the range is appropriately increased by 10%-15% to improve the spatial continuity of interpolation. Finally, the weighted least squares method is used to fit the variogram of the target displacement data of coherent points in each region, with the objective function being to minimize the sum of squared residuals, ultimately forming a dedicated variogram parameter set for each motion mode region.

[0087] Furthermore, based on synergetics theory, a three-level weight allocation system is constructed, consisting of data credibility calculation, motion mode weight adjustment, and dynamic correction, to fully explore the synergistic advantages of GB-SAR and GNSS data.

[0088] First, the reliability of the two types of data is accurately calculated, with the core being error quantization and normalization: ① Calculation of GNSS data reliability ω_gnss: First, quantify the GNSS displacement error σ_gnss, which consists of three parts: observation error, solution error, and environmental interference error, synthesized using the root mean square method; the total error σ_total is the sum of the maximum and minimum values ​​of the errors of the two types of data; the final reliability ω_gnss = 1 - σ_gnss / σ_total, with a value range of [0,1]. The larger the value, the more reliable the GNSS data. ② Calculation of GB-SAR three-dimensional displacement conversion error σ_gb: This includes projection coefficient matrix error (σ_A), phase unwrapping error (σ_unw), and conversion model error (σ_mod), satisfying: Similarly, the GB-SAR data reliability ω_gb=1-σ_gb / σ_total is calculated, and ω_gnss+ω_gb=1 is satisfied.

[0089] Then, the weights are adjusted based on the motion pattern partitioning to adapt to the sensor data characteristics of different deformation areas in the mine:

[0090] ① Stable region: Deformation is weak and spatial continuity is strong, GB-SAR is more advantageous in area coverage. The weights are adjusted to ω_gb=0.7 and ω_gnss=0.3 to enhance the ability of area data to characterize global trends; ② Accelerated deformation region: Deformation rate is fast and local abrupt changes are easy to occur. GNSS absolute displacement accuracy is higher (able to accurately capture single-point abrupt changes). The weights are adjusted to ω_gb=0.3 and ω_gnss=0.7 to highlight the constraining effect of high-precision point data; ③ Uniform deformation region: Deformation is stable and has both spatial continuity and single-point reliability. The balanced weights ω_gb=0.5 and ω_gnss=0.5 are adopted to achieve complementary advantages of the two types of data.

[0091] Finally, a dynamic correction mechanism is introduced. Based on the surface displacement gradient field, an additional GNSS data weight of 0.1 is added to the gradient high value region (deformation concentration region) (e.g., ω_gnss is adjusted to 0.8 in the gradient high value region of accelerated deformation region), and a GB-SAR data weight of 0.1 is added to the gradient low value region (deformation uniform region) (e.g., ω_gb is adjusted to 0.8 in the gradient low value region of stable region). At the same time, the consistency of the two types of data is monitored in real time, and the 3D displacement difference between GNSS and GB-SAR at the same monitoring point is calculated. If the difference is >0.5mm, the weight of the data with lower confidence is temporarily reduced by 0.1 to ensure that the weight allocation can dynamically adapt to changes in data quality.

[0092] Based on this, the absolute displacement field calibrated by GNSS was obtained through Kriging synergetics fusion. Based on the zonal variogram, with GNSS control points as the core, the monitoring area was divided into uniform grids matching the resolution of the mine slope DEM. Interpolation calculations were performed separately for regions with different motion modes: ordinary Kriging interpolation was used in the stable region, universal Kriging interpolation was used in the accelerated deformation region, and co-Kriging interpolation was used in the uniform deformation region. During the interpolation process, the interpolation error of each grid was recorded in real time. For grid points with large errors, the variogram parameters were readjusted and interpolation was performed a second time until the error met the requirements, and finally, the preliminary GNSS interpolated absolute displacement field was obtained. Furthermore, for each grid point, the preliminary GNSS interpolated displacement value S_gnss and GB-SAR three-dimensional displacement value S_gb are extracted and fused according to adaptive weights. The fusion formula is S_fuse=ω_gnss·S_gnss+ω_gb·S_gb. The final output GNSS calibration absolute displacement field contains time-by-time three-dimensional absolute displacement data of the entire monitoring area, combining the high precision of GNSS and the area coverage advantage of GB-SAR.

[0093] S4. Using multi-feature fuzzy comprehensive evaluation, the displacement period characteristic parameters and surface displacement gradient field are processed to obtain the comprehensive index field of local slope stability.

[0094] Extracting the displacement period characteristic parameters includes the following steps:

[0095] S411. Obtain the denoised GNSS coordinate time series data by using gross error identification and multi-criteria robust filtering.

[0096] In this embodiment, the interference period is first accurately marked. Based on the mine production log (which must include the detonation time, explosive quantity, and impact range of the blasting operation, the route, time period, and load of large transport vehicles, as well as the specific time window for equipment maintenance), a time + space dual-dimensional marking rule is adopted: in terms of time, the 30 minutes before and after the blasting operation is marked as the interference period, and the large transport is marked as the fixed daily travel time and GNSS monitoring points within 50m along the transport route; in terms of space, combined with the distance between the monitoring point coordinates and the interference source (blasting point, transport route), the monitoring point data with a distance of <100m is additionally marked as highly suspicious interference data.

[0097] Then, a multi-criteria gross error identification system was implemented to form a comprehensive gross error screening system: Criterion 1 is the displacement mutation threshold criterion, which combines the statistical characteristics of historical displacement data of the mine in recent years (calculating the mean and standard deviation of displacement change rate in different regions and seasons), and adopts regional dynamic threshold setting. The threshold for dangerous areas such as slope toe and crack development is set at 2.5 mm / hour, and the threshold for areas such as slope top and stable platform is set at 3 mm / hour. When the displacement change rate exceeds the threshold of the corresponding area at a certain moment, it is judged as a suspicious gross error; Criterion 2 is the residual threshold criterion, which uses the stable anchor point set as a benchmark, first performs Kalman filtering preprocessing on the anchor point GNSS data, extracts the filter residual sequence, and calculates the 3σ value of the residual as the global residual threshold (σ is the standard deviation of the anchor point residual sequence). When the filter residual of the monitoring point exceeds this threshold, it is marked as a suspicious gross error; Criterion 3 is the interference period association criterion, which temporarily judges the GNSS data in the marked interference period as a suspicious gross error, regardless of whether its displacement or residual exceeds the standard, and then weakens its influence through subsequent weight adjustment.

[0098] Next, a robust Kalman filter algorithm adapted to the mining scenario is used for gross error removal and denoising. Weight levels are assigned based on the intensity of interference type: data during blasting periods are weighted at 0.2, data during large-scale transportation periods at 0.5, data during power transmission line electromagnetic interference periods at 0.6, and data during normal, interference-free periods at 0.8. The filtering state equation and observation equation are adjusted for mining adaptation. The state equation incorporates geological stiffness parameters of the mine slope to constrain the reasonable range of displacement change rate, avoiding deformation information loss due to excessive smoothing. The observation equation incorporates a multipath effect correction term based on the surrounding terrain of the monitoring point, such as the distribution of deposits and vegetation. The filtering process is iterative. After each iteration, the weight iteration coefficient of suspicious gross errors is calculated. When the weight coefficient is <0.1, the data point is directly removed until the filtering residual converges, finally outputting the denoised GNSS coordinate time-series data.

[0099] S412. Separate the displacement trend term and period term from the GNSS coordinate time series data, and then extract the displacement period feature parameters.

[0100] First, a precise fitting of the trend term is performed. Considering that the displacement of the mine slope may have linear and nonlinear trends (such as linear deformation in the early stage of mining and accelerated nonlinear deformation in the later stage), the piecewise least squares method is used to fit the denoised GNSS displacement sequence: taking the mining face advancement node and major precipitation event as time boundary points, the entire displacement time series is divided into multiple continuous analysis periods. Linear or quadratic polynomial fitting is performed for each period. The piecewise curves obtained by fitting are superimposed to form a complete displacement trend term. This trend term intuitively reflects the deformation trend of the slope under the long-term effects of mining, geological stress and other factors.

[0101] The fitted residual sequence was then preprocessed and its stationarity was tested: First, the residual sequence was smoothed using a 5-point moving average method to remove random noise interference, resulting in a smoothed residual sequence; The stationarity of the residual sequence was determined using the ADF test. During the test, the optimal lag order was determined based on the time scale of the mining data using the AIC criterion. If the P-value of the ADF test was <0.05, the sequence was considered stationary and periodic terms could be extracted directly; if the P-value was ≥0.05, it indicated that the sequence had a unit root, and the residual sequence needed to be differencing by first or second order until the sequence was stationary (the ADF test was performed again after differencing to ensure that the P-value was <0.05). Considering that the displacement of mine slopes is affected by environmental factors such as precipitation and freeze-thaw cycles, and that the periodic characteristics are non-stationary (e.g., the amplitude of precipitation cycle deformation varies from year to year), Empirical Mode Decomposition (EMD) combined with mine environmental data is used to extract periodic terms. First, EMD decomposition is performed on the stationary residual sequence to obtain multiple intrinsic mode functions (IMFs) and one residual term, with each IMF representing a fluctuation component of a different frequency. Combined with the temporal environmental data of the mining area (daily precipitation, weekly average temperature, and start and end times of freeze-thaw periods), the correlation coefficient between each IMF and the environmental data is calculated (using Pearson correlation analysis), and IMFs with an absolute correlation coefficient ≥ 0.6 are selected as candidate periodic terms. To avoid mode aliasing affecting the accuracy of period extraction, the selected candidate periodic terms are further decomposed and optimized using the EnsembleEMD (EEMD) method, adding Gaussian white noise with an amplitude of 0.2 times the standard deviation of the original sequence, and repeating the decomposition 200 times to eliminate mode aliasing. Finally, the main periodic components were determined by power spectrum analysis (using the Welch method, with a frequency resolution of 0.001 Hz and the Hanning window as the window function): power spectrum calculation was performed on the optimized candidate periodic terms, the frequency corresponding to the peak value of the power spectrum was identified, the frequency was converted into the period length (T=1 / frequency), the rationality of the period was verified in combination with the characteristics of the mining environment, and the main periodic components that need to be retained were finally determined.

[0102] Furthermore, individual feature fitting is performed on the main periodic terms (such as annual and quarterly periodic terms). Considering that the periodic deformation caused by the mining environment may have a decreasing or increasing trend in amplitude, a variable parameter sinusoidal model is used for fitting: For each main periodic term, the model Δs_period,i=A_i(t)·sin(2πt / T_i+φ_i) is constructed, where A_i(t) is the amplitude that changes with time, and is fitted with a linear function A_i(t)=a_i·t+b_i, where a_i is the rate of change of amplitude, reflecting the strengthening or weakening trend of periodic deformation, T_i is the period length, and φ_i is the phase.

[0103] In this embodiment, the Levenberg-Marquardt algorithm is used for parameter optimization, with the objective function being the minimization of the sum of squares of the fitted residuals. During the iteration process, parameter constraints are set as follows: the absolute value of the amplitude change rate a_i ≤ 0.02 mm / day (in combination with historical mine data to avoid unreasonable large amplitude changes), and the phase φ_i ∈ [0, 2π) to ensure that the physical meaning of the parameters is reasonable. The iteration convergence condition is set as the parameter change between two adjacent iterations < 0.001, or the change in the sum of squares of the residuals < 0.01. The upper limit of the number of iterations is set to 100. If convergence is not achieved, the initial parameters are adjusted and refitted.

[0104] Then, the fitting parameters of all major periodic terms are integrated to form a complete set of periodic feature parameters: for each periodic term, the period length T_i, initial amplitude b_i, amplitude change rate a_i, phase φ_i, and goodness of fit R are recorded. 2 Simultaneously, the contribution weight of each periodic term is calculated to clarify the degree of influence of different periodic components on slope displacement.

[0105] Finally, a complete displacement component decomposition model is constructed: the fitted trend term Δs_trend, the fitted sequence of each major periodic term Δs_period (all Δs_period,i superimposed), and the residual random term Δs_random after EMD decomposition are integrated to obtain Δs_total=Δs_trend+Δs_period+Δs_random, thereby obtaining the displacement periodic characteristic parameter set of trend term parameters (fitting equation coefficients, goodness of fit), periodic characteristic parameters (T_i, A_i(t), φ_i, contribution weight of each period), and random term statistical parameters (mean, standard deviation).

[0106] The analysis of the surface region displacement gradient field includes the following steps:

[0107] S421. Select relevant target points by region to obtain a preliminary set of relevant target points.

[0108] First, the original GB-SAR coherence coefficient map sequence is preprocessed to eliminate mine scene-specific interference: median filtering (3×3 window) is used to smooth the random noise in the coherence coefficient map (mine dust and equipment electromagnetic interference can easily cause coherence coefficient jumps), and then histogram equalization is used to enhance the dynamic range of the coherence coefficient and improve the distinction between low coherence and high coherence regions.

[0109] Subsequently, a regional dynamic threshold screening strategy was implemented. The threshold setting needs to be precisely adapted to the geological and deformation characteristics of the mine slope. In the example, the dangerous area (slope toe, crack development area, and above the goaf) is the core monitoring area for landslide hazards. It is necessary to prioritize ensuring the distribution density of coherent point targets to capture subtle deformations. Therefore, the coherence coefficient threshold is lowered to 0.6. Even if some low coherence points are included, invalid points can be removed through subsequent network optimization. The stable area (slope top platform, non-mining affected area) has weak deformation. It is necessary to focus on the quality of coherent point targets to avoid error accumulation. Therefore, the threshold is set to 0.7, and only reliable points with high coherence are retained.

[0110] To further improve the purity of the coherent point target set, it is necessary to combine multi-source data to remove low-quality points: based on slope DEM data, identify densely vegetated areas and loose deposit areas, mark coherent points in these areas and remove them directly; combine GB-SAR time-series phase data to calculate the phase stability of each candidate coherent point target (phase standard deviation < 0.2 rad is considered stable), and remove points with drastic phase fluctuations.

[0111] The final preliminary set of coherent point targets needs to be correlated and verified with atmospherically corrected GB-SAR relative displacement data: extract the displacement sequence of each coherent point target, and use the 3σ criterion to remove coherent point targets with a displacement jump rate greater than 5%, ensuring that the displacement data of the included coherent point targets has been freed from atmospheric interference, so as to provide a high-quality sample basis for subsequent network construction and gradient extraction.

[0112] S422. Construct a spatial adjacency network of coherent point targets based on the preliminary coherent point target set, and then quantify the relative displacement differences of adjacent coherent point targets.

[0113] First, the accurate calculation of the three-dimensional spatial coordinates of the coherent point targets is completed. Based on the high-precision DEM of the mine slope, the planar coordinates (x, y) of each coherent point target on the DEM are extracted. Combined with the elevation information (z) of the DEM, the three-dimensional coordinates (x, y, z) of the coherent point targets are obtained. Considering the characteristics of the mine terrain with severe undulations, the quadratic polynomial interpolation method is used to optimize the DEM elevation extraction accuracy and avoid elevation errors caused by the DEM sampling interval.

[0114] Secondly, the Delaunay triangulation method, adapted to the irregular terrain of the mine, is used to construct the adjacency network. By automatically generating the optimal triangular mesh, it is ensured that each coherent point target forms a stable connection with 3-6 adjacent coherent point targets, and the interior angles of the triangles are all greater than 30°, avoiding the distortion of neighborhood relationships caused by narrow triangles. During the construction process, terrain constraints are also introduced, that is, the elevation difference between adjacent coherent point targets must not exceed 5m. If the elevation difference exceeds the standard, the connection is forcibly disconnected and a suitable adjacent coherent point target is found again.

[0115] Finally, network structure optimization and isolated point handling are implemented: First, the number of connections for each coherent point target is counted. Coherent point targets with less than 2 connections are identified as isolated points, as these points cannot reflect the spatial transmission characteristics of displacement and need to be removed. For coherent point targets with 2 connections (edge ​​points), the partitioning window is expanded. If the number of connections still cannot be increased, they are marked as valid edge points, and their weights are reduced in subsequent gradient calculations. At the same time, abnormal connections in the network are checked (such as spatial distance between adjacent coherent point targets > 20m, exceeding the spatial correlation range of local deformation of the mine slope), and abnormal connections are forcibly disconnected and the local network is reconstructed.

[0116] Furthermore, from the atmospherically corrected GB-SAR relative displacement sequence, the LOS direction displacement value of each coherent point target at each monitoring time is extracted, and linear interpolation is used to fill in the data for a small number of missing times. For displacement data during special periods such as mine blasting and rainstorms, the 3σ criterion is first used to remove outliers to avoid outliers affecting the calculation of neighborhood difference. Next, a set of neighborhood coherent point target pairs is constructed: based on the optimized spatial adjacency network of coherent point targets, all effective adjacent coherent point targets of each coherent point target are traversed to form neighborhood coherent point target pairs (i,j), and the spatial distance d_ij of each pair of coherent point targets is recorded. To avoid duplicate calculations, an i < j numbering rule is set to ensure that each neighborhood relationship is calculated only once. The third step is to calculate the relative displacement difference in the neighborhood: For each pair of coherent point targets (i,j) in the neighborhood, calculate the displacement difference Δs_ij(t)=s_i(t)-s_j(t) at each time step, where s_i(t) and s_j(t) are the LOS displacement values ​​of coherent point targets i and j at time t, respectively; during the calculation, the sign of the difference needs to be recorded synchronously. A positive difference indicates that the displacement of coherent point target i is greater than that of coherent point target j, and a negative difference indicates the opposite. The sign of the difference will be used to determine the deformation gradient direction in the subsequent calculation. The fourth step is to construct the neighborhood relative displacement matrix: using the coherent point target numbers as rows and columns, construct an N×N symmetric matrix (N is the total number of coherent point targets). The matrix elements are the time-by-time relative displacement sequences of the corresponding coherent point target pairs. For non-adjacent coherent point target pairs, the matrix elements are set to NaN (invalid values). At the same time, a separate spatial distance matrix is ​​constructed to record the spatial distance d_ij of each neighborhood coherent point target pair, providing basic data for subsequent deformation gradient calculation. The final output neighborhood relative displacement matrix needs to pass a consistency check: for different adjacent coherent point targets of the same coherent point target, calculate the correlation of their relative displacements. If the correlation coefficient is ≥0.7, it is considered consistent, ensuring the rationality of neighborhood displacement transmission.

[0117] S423. Analyze the surface deformation gradient field through the relative displacement of the neighborhood.

[0118] First, the deformation gradient of a single coherent point target is calculated: For each coherent point target, the relative displacement difference Δs_ij(t) and spatial distance d_ij of all its effective neighboring coherent point target pairs (i,j) are extracted. At each time step, the gradient component G_ij(t) = |Δs_ij(t)| / d_ij of the coherent point target in each neighborhood direction is calculated. The direction of the gradient component is the direction of the displacement difference (i.e., from coherent point target -j to coherent point target -i; if Δs_ij(t) is positive, then it is negative). Since a single coherent point target... There may be multiple neighborhood gradient components. The weighted average method is used to calculate the comprehensive deformation gradient of the coherent point target: the product of the coherence coefficients of each neighborhood coherent point target pair is used as the weight, and the weighted average is used to obtain the comprehensive gradient magnitude G_i(t) and gradient principal direction (the vector synthesis direction of all gradient components) of the coherent point target -i at time t. For coherent point targets whose gradient magnitude exceeds the reasonable range (G_i(t) > 5mm / m, set in combination with historical deformation data of the mine), they are identified as anomalous points, and their gradient values ​​are replaced by the average gradient of the three adjacent coherent point targets.

[0119] Then, surface deformation gradient field interpolation and post-processing are performed: using the comprehensive deformation gradient magnitude and main direction of each coherent point target as sample points, the Kriging interpolation method is used to extend the discrete coherent point target gradient data to a uniform grid of the entire monitoring surface area, so as to obtain the surface deformation gradient magnitude field and direction field at each time step; the interpolated gradient field is smoothed by Gaussian filtering to eliminate interpolation noise, while retaining the detailed features of the high gradient value area.

[0120] The final output of the surface displacement gradient field should include a gradient magnitude distribution map and a direction distribution map. The high gradient magnitude area can be directly marked as the deformation concentration area, and the convergence or divergence characteristics of the gradient direction can reflect the deformation transmission trend, providing core surface deformation basic data for subsequent local stability evaluation.

[0121] The local stability of mine slopes is influenced by multiple factors (displacement period characteristics, deformation gradient, motion pattern, etc.), and the degree of influence of each factor is ambiguous (e.g., high gradient, large period, and large amplitude are not clearly defined). This step adopts a multi-feature fuzzy comprehensive evaluation method, combined with mine geological weights, to perform fuzzy processing and comprehensive scoring on each evaluation index, thereby obtaining a comprehensive index field for local slope stability and achieving a surface-domain quantitative evaluation of stability.

[0122] Specifically, the method of using multi-feature fuzzy comprehensive evaluation to process displacement period feature parameters and surface displacement gradient field to obtain a comprehensive index field of local slope stability includes the following steps:

[0123] S431. An evaluation index system is constructed using displacement period characteristic parameters and surface displacement gradient field.

[0124] In this embodiment, the evaluation metrics include:

[0125] Periodic amplitude A (reflecting the intensity of periodic deformation): The annual average periodic amplitude A_year and the seasonal average periodic amplitude A_season are selected as core parameters. Considering the differences in periodic deformation sensitivity in different areas of the mine, an additional monthly periodic amplitude component A_month is added to the loose accumulation area at the toe of the slope. The calculation method is A=(A_year+A_season) / 2+k·A_month, where k is the regional correction coefficient. k=0.3 in the loose area and k=0.1 in the hard rock area. This index can quantify the cumulative effect of repeated deformation caused by environmental factors. The larger the amplitude, the weaker the slope's resistance to periodic deformation.

[0126] Deformation gradient G (reflecting the degree of deformation concentration): Based on the surface displacement gradient field, the gradient magnitude of the principal direction is selected for each evaluation unit. Considering the characteristics of mine structural surface development, if the evaluation unit is located within ±15° of the structural surface strike, the gradient value needs to be multiplied by an amplification factor of 1.2 (structural surfaces are prone to becoming deformation concentration channels). This index directly points to the potential sliding zone of the landslide, and areas with high gradient values ​​are often areas with weak stability.

[0127] Displacement amplitude S (reflecting the cumulative deformation scale): Based on the absolute displacement field calibrated by GNSS, a combined calculation method is adopted, which combines the cumulative displacement over the past 3 months with the predicted displacement rate trend, i.e., S = S_3month + v·30 (v is the average displacement rate over the past month, in mm / d). This method takes into account both the current cumulative deformation and the short-term development trend, avoiding the lag caused by judging solely based on historical data, and adapting to the dynamic deformation characteristics under the continuous impact of mining activities.

[0128] Movement pattern category M (reflecting the deformation development trend): Based on the movement pattern zoning map, a quantitative value is assigned. The basic assignment rule is: stable mode = 1, uniform deformation = 2, accelerated deformation = 3, and abrupt deformation = 4. To strengthen the weight of mining impact, the uniform / accelerated deformation modes in the mining-affected area are assigned values ​​of 2.5 / 3.5 respectively, while the abrupt deformation mode remains unchanged at 4. This indicator can intuitively distinguish the danger level of deformation. The larger the value, the higher the landslide risk.

[0129] After all indicators are extracted, they need to be normalized (mapped to the [0,1] interval) to eliminate the difference in units and lay the foundation for subsequent fuzzy evaluation.

[0130] S432. Determine the index weights by combining the analytic hierarchy process (AHP) and mine geological constraints.

[0131] A three-dimensional weight determination system combining the Analytic Hierarchy Process (AHP), mine geological constraints, and dynamic adjustment is adopted to ensure that the weight allocation aligns with the dominant stability factors in different regions.

[0132] First, a complete hierarchical structure is constructed: the target layer consists of the evaluation weights for local slope stability, the criterion layer consists of the three core influencing factors of mines: lithology type, degree of mining impact, and structural surface development level, and the indicator layer consists of the above four evaluation indicators.

[0133] Then, the expert consultation and judgment matrix was constructed: a review group was formed by inviting experts in mine geological exploration, monitoring technology, and safety management. Based on the specific geological report of the mine (including lithological distribution map, structural plane occurrence, and goaf location), the importance of each factor in the criterion layer and each indicator in the index layer was compared pairwise. For example, in areas with weak lithology, the degree of development of structural planes has a significant impact on stability. Therefore, the weight of structural planes (0.4) in the criterion layer is higher than that of mining influence (0.3) and lithology itself (0.3). Thus, it was deduced that the deformation gradient (structural plane dominance) has the highest weight in the index layer.

[0134] Third, the consistency test and weight calculation of the judgment matrix: The weight vector of the judgment matrix is ​​solved by the sum-product method, and the consistency ratio (CR) test is passed. If the consistency ratio is not met, experts are organized to readjust the comparison results. In the example, the final output of the regional weight scheme is as follows: ① Weak lithology + developed structural planes + strong mining influence area (such as loose roof and side slopes of coal mines): deformation gradient weight 0.35, displacement amplitude weight 0.3, period amplitude weight 0.2, movement mode weight 0.15 (deformation concentration and cumulative effect are the core risk sources); ② Hard lithology + undeveloped structural planes + non-mining influence area (such as integral rock slopes of metal mines): deformation gradient weight 0.25, displacement amplitude weight 0.25, period amplitude weight 0.25, movement mode weight 0.25 (balanced influence of various factors); ③ Intermediate transition area (medium lithology + relatively developed structural planes + weak mining influence area): deformation gradient weight 0.3, displacement amplitude weight 0.28, period amplitude weight 0.22, movement mode weight 0.2 (taking into account deformation concentration and development trend).

[0135] S433. Based on the evaluation index system and the index weights, construct a fuzzy evaluation matrix, and generate a comprehensive index field for local slope stability through fuzzy comprehensive evaluation.

[0136] Based on the fuzzy characteristics of mine slope stability evaluation, a targeted membership function is constructed for each evaluation index through piecewise linear membership functions and regional adaptation correction. This clarifies the degree of membership of each index value to the four evaluation levels of stable, relatively stable, unstable, and extremely unstable, ultimately forming an m×n fuzzy evaluation matrix (where m is the number of evaluation indicators and n is the number of evaluation levels).

[0137] First, deformation gradient G (unit: mm / m): Based on the statistics of historical landslide cases in mines, a regional membership function is constructed: For weak lithology zones: when G < 0.5, the membership is stable 0.9, relatively stable 0.1, unstable 0, and extremely unstable 0; when 0.5 ≤ G < 1.0, the membership is stable 0.2, relatively stable 0.6, unstable 0.2, and extremely unstable 0; when 1.0 ≤ G < 2.0, the membership is stable 0, relatively stable 0.2, unstable 0.6, and extremely unstable 0.2; when G ≥ 2.0, the membership is stable 0, relatively stable 0, unstable 0.1, and extremely unstable 0.9; the threshold for hard lithology zones is generally increased by 50% (e.g., G < 0.75 is stable), to suit the strong deformation resistance of hard rock.

[0138] Second, periodic amplitude A (unit: mm): Soft rock area: A < 2 is stable 0.9; 2 ≤ A < 5 is relatively stable 0.7; 5 ≤ A < 8 is unstable 0.6; A ≥ 8 is extremely unstable 0.9; Hard rock area: A < 5 is stable 0.9; 5 ≤ A < 10 is relatively stable 0.7; 10 ≤ A < 15 is unstable 0.6; A ≥ 15 is extremely unstable 0.9.

[0139] Third, displacement amplitude S (unit: mm): Mining-affected area: S < 10 is considered stable 0.9; 10 ≤ S < 20 is considered relatively stable 0.6; 20 ≤ S < 30 is considered unstable 0.7; S ≥ 30 is considered extremely unstable 0.9; Non-mining-affected area threshold floating: S < 20 is considered stable 0.9, and so on.

[0140] Fourth, the motion mode category M: When M=1 (stable), the membership is 0.9 for stable and 0.1 for relatively stable; when M=2 (uniform speed), the membership is 0.7 for relatively stable and 0.3 for unstable; when M=2.5 (uniform speed), the membership is 0.4 for relatively stable and 0.6 for unstable; when M=3 (acceleration), the membership is 0.7 for unstable and 0.3 for extremely unstable; when M=3.5 (accelerated speed), the membership is 0.3 for unstable and 0.7 for extremely unstable; when M=4 (abrupt change), the membership is 0.9 for extremely unstable. After constructing the membership function, a typical region with known stability states needs to be selected for verification to ensure the accuracy of the function in distinguishing stability levels. If the standard is not met, the threshold parameters should be adjusted and the function reconstructed.

[0141] Furthermore, a weighted average fuzzy synthesis + partitioned interpolation optimization is used to achieve accurate conversion from evaluation unit scores to surface domain index fields.

[0142] The first step is evaluation unit division and data matching: Based on the high-precision DEM of the mine slope, the monitoring area is divided into uniform evaluation units (matching the absolute displacement field resolution). Each unit is associated with the membership vector of the four indicators at the corresponding location and the indicator weight vector of the region, ensuring that the weights are accurately matched with the geological characteristics of the unit.

[0143] The second step is to calculate the fuzzy comprehensive score: a weighted average fuzzy synthesis operator is used, which satisfies: B=W·R, where W is a 1×4 index weight vector, R is a 4×4 fuzzy evaluation matrix, and B is a 1×4 comprehensive membership vector; the stability comprehensive score is calculated through the B vector: S_score=0.9·B1+0.65·B2+0.35·B3+0.1·B4 (B1-B4 are the membership degrees of stable, relatively stable, unstable, and extremely unstable, respectively, and the coefficients are the benchmark scores of the corresponding levels), and the score range is mapped to [0,1], with 1 being the most stable and 0 being the most unstable.

[0144] The third step is surface interpolation and post-processing: using the comprehensive score of each evaluation unit as sample points, a preliminary comprehensive index field of surface stability is generated using the Kriging interpolation method; during interpolation, the variogram parameters are adjusted in conjunction with the motion mode zoning to improve the interpolation accuracy of deformation-sensitive areas; the interpolation results are post-processed and optimized: ① Small area smoothing: median filtering with a 3×3 window is used to eliminate isolated high / low score points and avoid noise interference; ② Boundary correction: the boundary of the index field is made to match the boundary of the geological zoning (lithology, structural plane) and the mining-affected area to eliminate cross-regional score abrupt changes; ③ Mining verification: for the area surrounding the mining working face, the score is lowered by 0.1-0.2 (mining exacerbates instability) to ensure that the score matches the actual impact.

[0145] The fourth step is stability level classification and verification: Based on the comprehensive score, four stability levels are defined: 0-0.3: Extremely unstable (emergency remediation required immediately); 0.3-0.6: Unstable (requires intensified monitoring and specialized investigation); 0.6-0.8: Relatively stable (routine monitoring sufficient); 0.8-1.0: Stable (no additional intervention required). Historical verification points are selected for accuracy verification to achieve the required level determination accuracy. The final output is a comprehensive area stability index field containing stability level zones, comprehensive score contour lines, and hazard area markers, providing precise spatial positioning data for landslide hazard investigation.

[0146] S5. Combining GNSS early warning levels, GNSS-calibrated absolute displacement field, and local stability comprehensive indicators, high-precision absolute displacement field monitoring and error self-correction are achieved by utilizing graded early warning and displacement field source analysis.

[0147] Landslides on mine slopes typically exhibit abrupt changes in both displacement and displacement rate. Early warning systems based on a single parameter (displacement or velocity) are prone to false alarms (e.g., large short-term displacement but stable velocity, which may indicate periodic fluctuations) or missed alarms (e.g., a sudden change in velocity but displacement not reaching the threshold). This step, based on the geological conditions of the mine slope (lithology, structural planes) and historical landslide cases, sets a dual-parameter early warning threshold for displacement and velocity, extracts the early warning level sequence from GNSS monitoring points, and implements tiered early warning systems.

[0148] First, a geological weighting system was constructed, and the weights of core influencing factors were determined using the analytic hierarchy process (AHP). A three-level hierarchical structure was established: target layer (early warning threshold), criterion layer (lithology, degree of development of structural planes, slope height), and scheme layer (different geological grades). Experts in mining geology and monitoring were invited to conduct pairwise comparisons of the importance of each factor in the criterion layer to construct a judgment matrix. The rationality of the judgment matrix was verified through consistency checks, and finally, the weights of each geological factor were calculated (typical weight allocation: lithology 0.4, structural planes 0.35, slope height 0.25).

[0149] Secondly, based on historical cases, a basic threshold is fitted: hierarchical clustering is performed according to the weight of geological factors, and the cases are divided into three typical scenarios: weak lithology - developed structural planes - medium lithology of high slopes - relatively developed structural planes - hard lithology of medium slopes - undeveloped structural planes - low slopes; for each scenario, a logarithmic fitting method is used to construct a displacement-velocity threshold curve, and the inflection point of the curve is extracted as the basic threshold.

[0150] Next, a dynamic correction mechanism is introduced: the threshold is lowered by 20% when the rainfall in the past month is >100mm, and by 15% during the freeze-thaw period, based on recent mine meteorological data and the impact range of mining activities (the threshold is lowered by 25% within 50m of the mining face). Simultaneously, a strict periodic fluctuation deduction process is implemented: from the set of periodic characteristic parameters, the periodic term fitting sequence Δs_period(t) and the periodic rate sequence v_period(t) of each GNSS monitoring point are matched to ensure that the timestamps of Δs_obs (measured displacement) are completely synchronized with Δs_period(t), and v_obs (measured rate) are completely synchronized with v_period(t). Abnormal displacement and abnormal rate are calculated using Δs_abn = Δs_obs - Δs_period(t) and v_abn = v_obs - v_period(t), completely eliminating false anomalies caused by normal periodic fluctuations such as rainfall and freeze-thaw cycles, and avoiding misjudging seasonal deformation as landslide warning signals.

[0151] Furthermore, based on the actual needs of mine safety management and with reference to relevant standards in the "Safety Regulations for Metal and Non-metal Mines," a four-level early warning system is established, clearly defining the judgment criteria and response requirements for each level:

[0152] Blue Alert (Attention Level): The judgment condition is that v_abn>v1 / 2 and Δs_abn≤Δs1 / 2 (e.g., in the case of weak lithology, v_abn>1mm / d and Δs_abn≤5mm). The response requirement is to increase the GNSS sampling frequency in the area and strengthen the frequency of manual inspections.

[0153] Yellow alert (warning level): The judgment condition is that v_abn>v1 and Δs_abn≤Δs1 (in weak scenarios, v_abn>2mm / d and Δs_abn≤10mm). The response requirement is to activate the regional monitoring emergency plan, set up warning signs, and prohibit non-monitoring personnel from entering.

[0154] Orange alert (warning level): The judgment conditions are v_abn>v1×1.5 and Δs_abn>Δs1 (v_abn>3mm / d and Δs_abn>10mm in weak scenarios). The response requirements are to organize the evacuation of surrounding workers, activate the slope treatment special team, and formulate an emergency response plan.

[0155] Red Alert (Emergency Level): The judgment conditions are v_abn>v1×2 and Δs_abn>Δs1×1.5 (v_abn>4mm / d and Δs_abn>15mm in weak scenarios). The response requirements are to activate the mine's Level I emergency response, seal off the entire slope monitoring area, and carry out 24-hour real-time monitoring and landslide emergency response.

[0156] During the early warning level matching process, a continuous time confirmation mechanism is introduced: when a single moment meets a certain level condition, it is temporarily marked as a suspected early warning; only when three consecutive monitoring moments (e.g., three consecutive hours under regular monitoring, or 30 consecutive minutes under high-frequency monitoring) meet the level condition can it be officially determined as the corresponding early warning level, avoiding false alarms caused by single-point data anomalies. Subsequently, the early warning levels of each monitoring point are extracted in chronological order to generate a GNSS early warning level sequence. The sequence content includes core information such as monitoring point number, monitoring time, abnormal displacement, abnormal velocity, early warning level, and response suggestions. Finally, the early warning level sequence is post-processed: a consistency check between adjacent time levels is adopted. If the level of a certain moment differs from the levels before and after by more than two levels (e.g., jumping directly from blue to red), anomaly investigation is carried out in combination with the sensor operating status at that moment (whether there is signal obstruction or equipment failure) and mining production activities (whether there is blasting or large equipment operation). The sequence is updated only after confirmation, ensuring the reliability of the early warning sequence and providing high-quality point-like early warning data for subsequent global hierarchical early warning fusion.

[0157] Furthermore, the method of combining GNSS early warning levels, GNSS-calibrated absolute displacement fields, and local stability comprehensive indicators, and utilizing graded early warning and displacement field source analysis to achieve high-precision absolute displacement field monitoring and error self-correction, includes the following steps:

[0158] S51. By combining GNSS early warning levels, absolute displacement field deformation, and local stability comprehensive indicators, graded early warning of mine slopes can be achieved.

[0159] In this embodiment, the global early warning rules are refined and dynamically adapted: ① Red Alert (Emergency Level): The judgment conditions must meet three constraints: there are GNSS points with red alerts at three consecutive times (to avoid false alarms at a single time), the stability index of the corresponding grid is <0.3, and the displacement rate is >5mm / d. If the area is a strong mining-affected zone, the displacement rate threshold is lowered by 20%, and the response requirements are to activate the slope emergency monitoring system and update the displacement data every 5 minutes; ② Orange Alert (Early Warning Level): two consecutive orange alert GNSS points + stability index <0.6 + displacement rate >3mm / d, the threshold for the rate in the mining-affected area is lowered by 15%, and the response requires supplementary manual monitoring of the crack width on the structural surface, twice a day; ③ Yellow alert (warning level): two consecutive yellow alert GNSS points + stability index <0.8 + displacement rate >1mm / d, the threshold for abnormal weather periods such as rainstorms and freeze-thaw cycles is lowered by 10%, and the response requires increasing the frequency of encrypted meteorological data collection, once per hour; ④ Blue alert (attention level): only blue alert GNSS points exist or the stability index is 0.6-0.8, and the response requires generating a regional deformation trend report weekly.

[0160] Furthermore, consistency verification of early warnings is required: select typical areas with known early warning status, verify that the accuracy of early warning level determination meets the standard, and if it does not meet the standard, adjust the threshold parameters and rematch to finally obtain a global hierarchical early warning zoning map.

[0161] S52. By analyzing the displacement field, identify displacement anomalies and error sources to achieve error self-correction.

[0162] Specifically, through a three-level tracing process of anomaly area definition, multi-dimensional source investigation, and cross-validation, the core causes of displacement anomalies are accurately identified, with the tracing priority being sensor failure > environmental interference > mining interference > actual deformation.

[0163] Precise definition of anomalous regions: Based on the absolute displacement field calibrated by GNSS, a threshold and spatial clustering methods are used to define anomalous regions. First, a displacement rate anomalous threshold is set, and grid points with rates exceeding the threshold are marked. Then, the DBSCAN clustering algorithm is used to aggregate discrete anomalous points into continuous anomalous regions and remove isolated anomalous regions.

[0164] Sensor fault tracing and troubleshooting: Based on real-time records of sensor operating status, establish fault judgment criteria. For GNSS monitoring points: if the PDOP value is >6, the number of satellites is <4, or the signal obstruction duration is >10 minutes at a certain moment, and the displacement of the corresponding grid point is abnormal, it is judged as a GNSS signal obstruction fault. For GB-SAR data: if the GB-SAR coherence coefficient in the abnormal area is <0.6 and the phase unwrapping residual is >0.3rad, it is judged as a GB-SAR signal interference fault. Fault judgment needs to be cross-verified through data from adjacent sensors.

[0165] Source tracing and investigation of mining-induced disturbances: Based on the mining plan and the distribution of goaf areas, calculate the horizontal distance d between the abnormal area and the mining face. If d < 50m (strong mining-induced disturbance area), and the deformation rate of the abnormal area is positively correlated with the working face advance rate (v_mining) (correlation coefficient ≥ 0.7), then it is determined to be a mining-induced disturbance.

[0166] Using empirical formulas for deformation caused by mining (k is the lithology correction coefficient, k=0.8 for weak rock and k=0.3 for hard rock) Quantify the impact of mining. If the relative error between the calculated impact of mining and the abnormal displacement is ≤15%, the cause of mining interference can be further confirmed.

[0167] Environmental disturbance source tracing and investigation: Based on the time series meteorological data of the mine, establish the correlation rules between environmental disturbance and deformation. Rainstorm disturbance: Abnormal deformation occurs within 24-72 hours after daily precipitation > 50 mm, and the deformation trend coincides with the soil saturation time after precipitation. Freeze-thaw disturbance: Abnormal deformation occurs during the freeze-thaw period (when the temperature fluctuates around 0℃), and the deformation cycle is consistent with the freeze-thaw cycle. Environmental disturbance needs to be verified by periodic characteristic matching.

[0168] Confirmation of actual deformation anomalies: After excluding all the above-mentioned interfering factors, and combining the comprehensive local stability index field of the slope, if the stability index of the abnormal area is <0.6 (unstable or below level), and the deformation trend shows an accelerating characteristic (displacement rate increases with time, rate of change >0.2mm / d), then... 2 The slope was determined to be abnormally deformed (potential landslide hazard); additional geological survey data is needed for cross-verification (such as the presence of structural surfaces or geologically weak points above mined-out areas in the abnormal area) to ensure the accuracy of the determination.

[0169] Based on the source analysis results, differentiated correction strategies are implemented for different causes, and a dynamic update mechanism for the model parameters of the preceding steps is established, forming a closed loop of source tracing, correction, and optimization.

[0170] Sensor fault area correction: Neighborhood interpolation plus trend constraint method is used for correction. For GNSS fault points: 3-5 effective GNSS monitoring points are selected in the surrounding area. Kriging interpolation is performed based on the partition variogram (spherical function for stable areas and exponential function for accelerated deformation areas). The interpolation error needs to be controlled within ≤0.3mm. For GB-SAR interference areas: GB-SAR three-dimensional displacement sequences of adjacent stable areas of the abnormal area are extracted. Linear trend extrapolation method is used to supplement the data of the abnormal period. The extrapolation residual needs to be verified by GNSS data (RMSE≤0.5mm). After correction, the spatial continuity of the abnormal area needs to be checked. When the displacement difference between adjacent grid points is >1mm, the 5-point moving average method is used for smoothing.

[0171] Mining-induced disturbance correction: Based on the mining-induced impact Δs_mining quantified during the source tracing stage, a mining-induced disturbance deduction model is constructed: Δs_corrected = Δs_obs - Δs_mining × w_mining (w_mining is the mining-induced impact weight; w_mining = 1.0 when d < 30m, and w_mining = 0.6 when 30m ≤ d < 50m). After deduction, the correction effect needs to be verified: the difference between the corrected displacement rate and the rate of the same lithology area in the non-mining zone is ≤ 0.5mm / d, ensuring effective stripping of mining-induced disturbance.

[0172] Environmental interference area correction: Combine the periodic feature parameter set to fit the environmental deformation component Δs_env (the superposition of periodic terms such as rainstorm, freeze-thaw, etc.) of the abnormal area, and achieve interference subtraction by Δs_corrected=Δs_obs-Δs_env; after subtraction, calculate the residual between the correction sequence and the stable trend term. If RMSE≤0.4mm, the correction is deemed qualified. If it does not meet the standard, the periodic feature parameters are re-optimized (such as adjusting the noise amplitude of EMD decomposition).

[0173] Statistical analysis was performed on the intermediate error data of each step. Atmospheric correction residuals: if the mean residual value is >0.5mm, the number of iterations of the atmospheric phase screen was adjusted (increased from 3 to 5) and the inverse distance weighted interpolation weight index was increased (±0.3). Projection transformation error: if RMSE is >0.3mm, the projection coefficient matrix of GNSS monitoring points was recalibrated and structural surface constraints were added. Fusion error: if the RMSE after fusion is >0.3mm, the adaptive weights of different motion mode regions were adjusted (±0.1). After parameter optimization, iterative verification was required until the intermediate errors of each step met the accuracy requirements (mean atmospheric correction residual value ≤0.5mm, projection transformation RMSE ≤0.3mm, fusion RMSE ≤0.3mm), thus achieving self-adaptive optimization of the model.

[0174] Please see Figure 2 In this embodiment, to efficiently execute the mine slope absolute displacement field monitoring and error self-correction method provided by the present invention, the present invention also provides a mine slope absolute displacement field monitoring and error self-correction system, including: an input device, an output device, a processor, and a memory. The input device, output device, processor, and memory are interconnected. The memory contains program instructions used for the steps of the mine slope absolute displacement field monitoring and error self-correction method. The mine slope absolute displacement field monitoring and error self-correction system of the present invention has a compact structure and stable performance, and can stably execute the mine slope absolute displacement field monitoring and error self-correction method of the present invention, further improving the overall applicability and practical application capability of the present invention.

[0175] In summary, this invention lays the foundation for high-precision relative displacement by iteratively correcting atmospheric phase screens to remove complex atmospheric interference; it achieves a breakthrough in GB-SAR dimensions through geometric constraint calibration, providing the core basis for three-dimensional displacement conversion; it leverages the complementary advantages of point and surface data through Krifin integration to obtain a high-precision absolute displacement field across the entire region; it quantifies the stability surface domain based on multi-feature fuzzy evaluation, accurately locating potential hazard areas; and it combines hierarchical early warning and source tracing analysis to complete monitoring, early warning, and error self-correction, ensuring data reliability and landslide early warning accuracy, thus providing scientific support for mine slope safety management.

[0176] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention, and they should all be covered within the scope of the present invention.

Claims

1. A method for monitoring absolute displacement field of mine slope and self-correcting errors, characterized in that, Includes the following steps: Based on the atmospheric phase screen iterative modeling and correction of GNSS anchor points, GB-SAR time-series differential phase map data is processed to obtain atmospheric-corrected GB-SAR relative displacement sequence. By using geometric constraint calibration to transform line-of-sight to three-dimensional displacement vector, we analyze GNSS three-dimensional coordinate time series data and GB-SAR pixel geometric information to obtain the projection coefficient matrix from GB-SAR LOS displacement to the real three-dimensional displacement. Based on Kriging co-interpolation and adaptive weighted fusion, the atmospheric-corrected GB-SAR relative displacement, projection coefficient matrix, and motion mode partitioning are processed to obtain the absolute displacement field of GNSS calibration. By using multi-feature fuzzy comprehensive evaluation, the displacement period characteristic parameters and surface displacement gradient field are processed to obtain the comprehensive index field of local slope stability. By combining GNSS early warning levels, GNSS-calibrated absolute displacement field, and local stability comprehensive indicators, and utilizing graded early warning and displacement field source analysis, high-precision absolute displacement field monitoring and error self-correction can be achieved. The atmospheric phase screen iterative modeling and correction based on GNSS anchor points processes GB-SAR time-series differential phase map data to obtain atmospherically corrected GB-SAR relative displacement sequences, including the following steps: Select GNSS anchor points and estimate the atmospheric delay at the anchor points; An initial atmospheric phase screen is constructed, and the atmospheric delay and the initial atmospheric phase screen are combined to obtain the atmospherically corrected GB-SAR relative displacement sequence; The geometric constraint calibration, which involves transforming the line-of-sight to a three-dimensional displacement vector, analyzes GNSS three-dimensional coordinate time-series data and GB-SAR pixel geometric information to obtain the projection coefficient matrix from GB-SAR LOS displacement to the true three-dimensional displacement. This includes the following steps: Extract the three-dimensional displacement of GNSS monitoring points and the local geometric features of the slope; An initial projection coefficient model is constructed, and the geometric constraints of the initial projection coefficient model are calibrated and the coefficients are optimized using the three-dimensional displacement and the local geometric features of the slope to obtain the projection coefficient matrix from the GB-SAR LOS displacement to the real three-dimensional displacement. The process of processing the atmospheric-corrected GB-SAR relative displacement, projection coefficient matrix, and motion mode partitions based on Kriging co-interpolation and adaptive weighted fusion to obtain the GNSS-calibrated absolute displacement field includes the following steps: The relative displacement of atmospherically corrected GB-SAR is transformed into three-dimensional displacement using the projection coefficient matrix; Based on the transformation results and motion mode partitioning, Kriging co-interpolation fusion is performed according to the partitioning variogram function and adaptive weights to obtain the absolute displacement field calibrated by GNSS. Extracting the motion pattern partition includes the following steps: Extract phase statistical features and mark the area affected by mining. Motion pattern clustering analysis is performed by combining the phase statistical features and labeling results to generate a motion pattern partitioning map; The method combines GNSS early warning levels, GNSS-calibrated absolute displacement field, and local stability comprehensive indicators, and utilizes graded early warning and displacement field source analysis to achieve high-precision absolute displacement field monitoring and error self-correction, including the following steps: By combining GNSS early warning levels, absolute displacement field deformation, and local stability comprehensive indicators, graded early warning of mine slopes can be achieved. By analyzing the displacement field, displacement anomalies and error sources can be identified, enabling error self-correction.

2. The method according to claim 1, wherein, The method of using multi-feature fuzzy comprehensive evaluation to process displacement period feature parameters and surface displacement gradient fields to obtain a comprehensive index field of local slope stability includes the following steps: An evaluation index system is constructed using displacement period characteristic parameters and surface displacement gradient field. The weights of the indicators were determined by combining the analytic hierarchy process (AHP) with the geological constraints of the mine. Based on the evaluation index system and the index weights, a fuzzy evaluation matrix is ​​constructed, and a comprehensive index field for local slope stability is generated through fuzzy comprehensive evaluation.

3. The method according to claim 2, wherein, Extracting the displacement period characteristic parameters includes the following steps: Gross error identification and multi-criteria robust filtering are used to obtain denoised GNSS coordinate time series data; The displacement trend term and period term are separated from the GNSS coordinate time series data, and then the displacement period feature parameters are extracted.

4. The method according to claim 2, wherein, The analysis of the surface region displacement gradient field includes the following steps: By filtering coherent point targets by region, a preliminary set of coherent point targets is obtained; Based on the preliminary coherent point target set, a coherent point target spatial adjacency network is constructed, thereby quantifying the relative displacement differences between adjacent coherent point targets; The deformation gradient field of the surface region is analyzed by relative displacement of the neighborhood.

5. A mine slope absolute displacement field monitoring and error self-correction system, characterized in that, The mine slope absolute displacement field monitoring and error self-correction system includes: an input device, an output device, a processor, and a memory. The input device, output device, processor, and memory are interconnected. The memory includes program instructions, which are used to execute the mine slope absolute displacement field monitoring and error self-correction method according to any one of claims 1-4.

Images