Multi-temporal InSAR remote sensing data ground surface monitoring and geological disaster early warning method

Abstract

The application discloses a kind of multi-period InSAR remote sensing data ground surface monitoring and geological disaster early warning method, it is related to ground surface monitoring technical field.The method includes simultaneously obtaining SAR satellite data, GNSS reference station coordinates, unmanned aerial vehicle LiDAR point cloud and weather data, constitutes multi-source data set;SAR data is interfered with processing, combined with DEM to eliminate terrain phase, annual deformation rate graph is obtained by phase unwrapping and time series analysis;GNSS, LiDAR and InSAR data are spatially registered and weighted fusion, and deformation rate field is generated;The field and slope, slope direction parameter are input into decision tree model to identify geomechanical deformation type;Future deformation trend is predicted using LSTM neural network, and the result is compared with early warning threshold to determine early warning level;Early warning level is matched to execute scheme, and early warning report containing visual map is generated.The application realizes closed-loop management from monitoring to early warning by multi-source data fusion and machine learning algorithm.

Description

Technical Field

[0001] This application relates to the field of surface monitoring technology, and more specifically, to a method for surface monitoring and geological disaster early warning using multi-phase InSAR remote sensing data. Background Technology

[0002] InSAR remote sensing data from multiple phases has been used for monitoring land surface changes and early warning of geological disasters, and its technological development has a rich history. Research began as early as the early 1980s, but due to the limited availability and low resolution of early radar satellite data, the application results were unsatisfactory. After 2007, with the maturation of the new generation of high-resolution SAR satellite technology internationally, domestic application exploration accelerated. Since 2019, integrated remote sensing methods centered on InSAR technology have been used for large-scale operational surveys in high-risk areas of geological disasters across the country, marking the transition of this technology from research to large-scale engineering application.

[0003] Existing technical solutions typically follow a relatively isolated processing paradigm. Their core process relies on a single type of SAR data (such as Sentinel-1 or LuTan-1) for routine registration, filtering, and computation to ultimately generate a deformation rate map. This approach has significant limitations: First, different technologies (such as InSAR and GNSS) often operate independently, lacking an effective deep fusion mechanism, resulting in a lack of complementary advantages and difficulty in eliminating errors. Second, deformation analysis is mostly based on static thresholds or simple models, lacking sufficient ability to predict deformation mechanisms, types, and future trends. Finally, the early warning output is relatively simplistic, making it difficult to achieve scientific risk classification and rapid, accurate response. These limitations result in existing technical solutions exhibiting characteristics of "data silos" and "shallow analysis."

[0004] In view of this, there is an urgent need for a method for surface monitoring and geological disaster early warning using multi-phase InSAR remote sensing data to solve the above-mentioned technical problems. Summary of the Invention

[0005] The main purpose of this application is to provide a method for surface monitoring and geological disaster early warning using multi-phase InSAR remote sensing data to solve the problem of insufficient prediction capability.

[0006] To achieve the above objectives, the first aspect of this application proposes a method for surface monitoring and geological disaster early warning using multi-phase InSAR remote sensing data, comprising: The target area is obtained by acquiring multiple SAR data from the same orbit collected by SAR satellites, three-dimensional coordinate point data of the reference station obtained by GNSS base station, digital surface model data generated by UAV scanning, and meteorological data for a preset time period. The multi-period SAR data from the multi-source data are combined to generate several interferometric pairs. The flat phase and terrain phase are eliminated sequentially using the DEM data in the digital elevation model to generate a differential interferogram. Then, the branch-cut phase unwrapping algorithm is used to solve the wrapped phase in the differential interferogram. Finally, the PS-InSAR and SBAS-InSAR dual algorithms are used to perform time series calculations to obtain the annual deformation rate map. Spatially align the three-dimensional coordinate point data, the digital surface model data, and the annual deformation rate map. Then, based on the weights assigned to the three-dimensional coordinate point data, the digital surface model data, and the annual deformation rate map respectively, calculate the fused deformation rate using a weighted formula to obtain the deformation rate field. The deformation rate field, slope and aspect parameters are input into a pre-trained decision tree classification model to obtain the geomechanical deformation type of the target area. Using a target LSTM neural network, the deformation rate of the target region within a preset future time period is predicted, thus obtaining the deformation rate prediction of the target region within the preset future time period. Based on the predicted deformation rate of the target area and the comparison between the geomechanical deformation type of the target area and the preset warning threshold, the warning level of the target area within a preset time period in the future is obtained; The warning level is matched with the preset warning level execution plan, and a report is generated to obtain the warning report of the execution plan.

[0007] Preferably, the step of obtaining multi-source data for the target area by acquiring multiple periods of SAR data from the same orbit collected by SAR satellites, three-dimensional coordinate point data of the reference station acquired by GNSS base stations, digital surface model data generated by UAV scanning, and meteorological data for a preset time period includes: The SAR satellite receiving module collects multiple SAR images of the same orbit at a preset revisit period to obtain the multi-phase SAR data. The three-dimensional coordinate point data of the base station are obtained through the GNSS base station module at a preset sampling rate; The target area is scanned along a preset flight path using the drone's LiDAR module, generating digital surface model data with a point cloud density greater than or equal to a preset density threshold. Meteorological data of the target area within a preset historical time period is collected through the meteorological data acquisition module; The multi-period SAR data, the three-dimensional coordinate point data, the digital surface model data, and the meteorological data are combined to obtain multi-source data for the target area.

[0008] Preferably, the steps of combining the multi-period SAR data from the multi-source data to generate several interferometric pairs, sequentially eliminating the flat land phase and terrain phase using DEM data from the digital elevation model to generate a differential interferogram, then using the branch-cut phase unwrapping algorithm to solve the wrapped phase in the differential interferogram, and using the PS-InSAR and SBAS-InSAR dual algorithms for time-series calculation to obtain the annual deformation rate map include: Using the InSAR time series solution algorithm, the multi-period SAR data are combined and paired according to a preset time baseline threshold to generate several interferometric pairs; Using DEM data from the digital elevation model, several interferometric pairs are sequentially processed to eliminate flat-ground phase effects and terrain phase effects, generating differential interferometric atlases. The phase unwrapping calculation of the differential interferogram is performed using the branch-cut phase unwrapping algorithm to obtain the unwrapped phase result; The PS-InSAR algorithm and the SBAS-InSAR algorithm are run respectively to perform time series analysis on the unwrapped phase results to obtain an annual deformation rate map. The PS-InSAR algorithm extracts the deformation of permanent scatterers in the unwrapped phase results, and the SBAS-InSAR algorithm performs deformation analysis on distributed scatterers in the unwrapped phase results.

[0009] Preferably, the step of spatially aligning the three-dimensional coordinate point data, the digital surface model data, and the annual deformation rate map, and then calculating the fused deformation rate using a weighted formula based on the weights assigned to the three-dimensional coordinate point data, the digital surface model data, and the annual deformation rate map to obtain the deformation rate field includes: Using transverse Mercator projection, the three-dimensional coordinate point data, the digital surface model data and the annual deformation rate map are spatially registered and aligned to obtain a data stack with unified spatial reference. Weights are assigned to the three-dimensional coordinate point data, the digital surface model data, and the annual deformation rate map in the unified data stack of the spatial reference, respectively. The fused deformation rate is then calculated using a weighted formula to obtain the deformation rate field. The formula for calculating the deformation rate field is as follows: ; Where V represents the deformation rate value of the deformation rate field. This represents the deformation rate value derived from InSAR technology. This represents the deformation value from GNSS measurements. This represents the terrain deformation value extracted from LiDAR data. , and These represent the weight coefficients of their respective data sources.

[0010] Preferably, the step of inputting the deformation rate field and the slope and aspect parameters into a pre-trained decision tree classification model to obtain the geomechanical deformation type of the target area includes: The slope parameters and aspect parameters are obtained using data from the digital elevation model. The deformation rate field, the slope parameter, and the aspect parameter are respectively input into a pre-trained decision tree classification model to obtain the geomechanical deformation type of the target area. The pre-trained decision tree classification model refers to a pre-trained decision tree classification model trained based on sample data of various historical geological disaster types.

[0011] Preferably, the step of using a target LSTM neural network to predict the deformation rate of the target region within a preset future time period, and obtaining the predicted deformation rate of the target region within the preset future time period, includes: Based on historical multi-source data within a preset time period of the target region, an LSTM neural network is trained to obtain the target LSTM neural network; Using the target LSTM neural network, the deformation rate of the target region within a preset time period is predicted to obtain the deformation rate prediction value.

[0012] Preferably, the step of comparing the predicted deformation rate of the target area and the geomechanical deformation type of the target area with a preset warning threshold to obtain the warning level of the target area within a preset future time period includes: Based on the geomechanical deformation type of the target area, the corresponding graded early warning threshold is called from the preset threshold rule library to obtain the graded early warning threshold range matching the geomechanical deformation type of the target area. Different geomechanical deformation types in the preset threshold rule library correspond to different threshold ranges. The deformation rate prediction of the target area is matched with the graded early warning threshold range to obtain the early warning level of the target area within a preset time period in the future.

[0013] Preferably, the step of matching the warning level with a preset warning level execution plan and generating a report to obtain the warning report of the execution plan includes: Construct a preset solution library, wherein the preset solution library includes an execution solution corresponding to each warning level; Based on the warning level, a matching scheme is performed in the preset scheme library to obtain the execution scheme for the target area; Based on the execution plan for the target area, a visual map product containing the spatial distribution of the warning area is generated, and a report is generated to obtain the warning report for the execution plan.

[0014] Secondly, this application provides a surface monitoring and geological disaster early warning system based on multi-phase InSAR remote sensing data, applied to the aforementioned method for surface monitoring and geological disaster early warning based on multi-phase InSAR remote sensing data, including: The acquisition unit is used to acquire SAR data from multiple periods of the same orbit collected by SAR satellites in the target area, three-dimensional coordinate point data of the reference station acquired by GNSS base station, digital surface model data generated by UAV scanning, and meteorological data for a preset time period, so as to obtain multi-source data of the target area. The deformation rate unit is used to combine the multi-period SAR data in the multi-source data to generate several interferometric pairs, use the DEM data in the digital elevation model to sequentially eliminate the flat phase and the terrain phase, generate a differential interferogram, and then use the branch-cut phase unwrapping algorithm to solve the wrapped phase in the differential interferogram. Finally, the PS-InSAR and SBAS-InSAR dual algorithms are used to perform time series calculations to obtain the annual deformation rate map. The rate field unit is used to spatially align the three-dimensional coordinate point data, the digital surface model data, and the annual deformation rate map. Then, based on the weights assigned to the three-dimensional coordinate point data, the digital surface model data, and the annual deformation rate map, a weighted formula is used to calculate the fused deformation rate to obtain the deformation rate field. The decision unit is used to input the deformation rate field and the slope and aspect parameters into a pre-trained decision tree classification model to obtain the geomechanical deformation type of the target area. The prediction unit is used to predict the deformation rate of the target region within a preset time period using a target LSTM neural network, thereby obtaining the deformation rate prediction of the target region within the preset time period. The threshold judgment unit is used to compare the predicted deformation rate of the target area and the geomechanical deformation type of the target area with a preset warning threshold to obtain the warning level of the target area within a preset time period in the future. The result unit is used to match the warning level with the preset warning level execution plan and generate a report to obtain the warning report of the execution plan.

[0015] Thirdly, this application provides a computer program that, when executed by a processor, implements the steps of the aforementioned method for surface monitoring and geological disaster early warning of multi-phase InSAR remote sensing data.

[0016] The technical solution provided in this application includes at least the following beneficial effects: In the surface monitoring and geological disaster early warning method using multi-phase InSAR remote sensing data described in this application, an integrated collaborative observation network of space (SAR satellites), air (UAVs), and ground (GNSS base stations) is constructed to achieve systematic acquisition and deep fusion of multi-source remote sensing data. A weighted fusion algorithm effectively suppresses errors from single data sources, significantly improving the accuracy and reliability of deformation monitoring. Furthermore, by introducing machine learning algorithms such as decision tree models and LSTM neural networks, intelligent identification of geological disaster types and accurate prediction of future deformation trends are achieved, enabling the analysis method to move from static description to dynamic forecasting and profoundly revealing the deformation mechanism. Finally, by establishing a dynamic early warning threshold library and a graded response mechanism associated with disaster types, closed-loop management from monitoring data to precise prevention and control actions is realized, generating a structured report integrating early warning levels, spatial distribution maps, and specific handling plans. This greatly improves the scientific nature, accuracy, and emergency response efficiency of geological disaster early warning, effectively overcoming the shortcomings of data silos and shallow analysis in existing technologies. Attached Figure Description

[0017] The accompanying drawings, which form part of this application, are used to provide a further understanding of the application and to make other features, objects, and advantages of the application more apparent. The illustrative embodiments and descriptions of this application are used to explain the application and do not constitute an undue limitation of the application. In the drawings: Figure 1 A flowchart of a method for surface monitoring and geological disaster early warning using multi-phase InSAR remote sensing data provided in this application. Detailed Implementation

[0018] To enable those skilled in the art to better understand the present application, the technical solutions in the embodiments of the present application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present application, and not all embodiments. Based on the embodiments in the present application, all other embodiments obtained by those of ordinary skill in the art without creative effort should fall within the scope of protection of the present application.

[0019] It should be noted that the terms "first," "second," etc., in the specification, claims, and accompanying drawings of this application are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate for the embodiments of this application described herein. Furthermore, the terms "comprising" and "having," and any variations thereof, are intended to cover non-exclusive inclusion; for example, a process, method, system, product, or apparatus that comprises a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0020] Furthermore, the terms "installation," "setup," "equipped with," "connection," "linked," and "socketing" should be interpreted broadly. For example, "connection" can be a fixed connection, a detachable connection, or an integral structure; it can be a mechanical connection or an electrical connection; it can be a direct connection or an indirect connection through an intermediate medium, or an internal connection between two devices, components, or parts. Those skilled in the art can understand the specific meaning of the above terms in this application according to the specific circumstances.

[0021] The following explanations of some terms used in this application are provided to aid in understanding the technical solution of this application: SAR: Synthetic Aperture Radar is an active microwave remote sensing technology that acquires surface phase and amplitude information by transmitting and receiving radar signals from satellites (such as Sentinel-1 and LuTan-1), and has all-weather imaging capabilities.

[0022] GNSS: Global Navigation Satellite Systems (such as GPS and BeiDou) provide three-dimensional coordinate positioning (sampling rate 1Hz) through base station modules, which is used for external compliance verification of InSAR deformation data to improve the accuracy and reliability of monitoring results.

[0023] LiDAR: LiDAR uses a module mounted on a drone to emit laser pulses to measure distance and generate high-precision point cloud data, which is used to obtain local terrain deformation information and assist in the verification of InSAR data.

[0024] DEM: Digital Elevation Model is raster data representing terrain elevation (ASTERGDEM, Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model). In InSAR processing, it is used to eliminate the flat-land phase effect caused by terrain and provide a basic terrain reference for deformation analysis.

[0025] Branch-cut phase unwrapping algorithm: The branch-cut method is a phase unwrapping algorithm that solves the phase ambiguity problem by identifying phase jump points and constructing branch-cut lines, avoiding the defects of the phase continuity assumption in traditional algorithms, thereby improving the unwrapping accuracy.

[0026] PS-InSAR: Permanent scatterer InSAR is used to monitor deformation at stable points such as buildings and bare rock, and extracts high-precision temporal deformation information. It is often used in urban areas or areas with exposed rock to supplement the shortcomings of SBAS-InSAR.

[0027] SBAS-InSAR: Small Baseline Set InSAR is a time-series InSAR technique suitable for processing non-permanent scatterers such as vegetation cover areas. It improves the stability of wide-area deformation monitoring by combining SAR data from multiple periods to monitor slow surface deformation.

[0028] Loop Closure: In InSAR data processing, loop closure verification is used to verify the accuracy of phase unwrapping. By checking the phase closure error in the interferogram, it ensures that the unwrapping error is controlled within the allowable range, thereby improving the reliability of deformation inversion.

[0029] UTM is a globally used map projection system, namely the transverse conformal cylindrical projection. It divides the Earth's surface (from 80°S to 84°N) into 60 longitude zones (6° each), establishing an independent Cartesian coordinate system for each zone, thereby minimizing terrain deformation within a limited area.

[0030] The Gamma correction model was originally a non-linear operation widely used in image processing, display technology, and computer graphics. It is mainly used to encode and decode light intensity or signal level to conform to the non-linear perception of light by the human visual system.

[0031] like Figure 1 As shown, in a first aspect, this application provides a method for surface monitoring and geological disaster early warning using multi-phase InSAR remote sensing data, including: S100: Acquire multi-source data of the target area by acquiring SAR data from multiple periods of SAR satellites in the same orbit, 3D coordinate point data of the reference station acquired by GNSS base station, digital surface model data generated by UAV scanning, and meteorological data for a preset time period.

[0032] Specifically, obtaining multi-source data for the target region may include the following steps: S101, through the SAR satellite receiving module, collects multiple SAR images of the same orbit according to a preset revisit period, and obtains multiple SAR data.

[0033] Specifically, the satellite orbit parameters are configured through the SAR satellite receiving module to ensure periodic observation of the same orbital area. The period can be set as needed, and multiple SAR images of the same orbit are collected according to a preset revisit cycle. In other words, the SAR satellite will repeatedly photograph the target area on Earth according to a fixed schedule, thereby forming multiple SAR data.

[0034] For example, multi-period SAR data can be Sentinel-1 data. The revisit period of Sentinel-1 satellite can be set to an interval of ≤12 days, the spatial resolution can be set to a resolution of ≥5m, and the commonly used polarization mode can be VV polarization.

[0035] S102 acquires the three-dimensional coordinate point data of the base station through the GNSS base station module at a preset sampling rate.

[0036] Specifically, multiple GNSS reference stations are deployed in the target monitoring area to form a deformation monitoring network. Next, the three-dimensional coordinate data (longitude, latitude, and elevation) of the reference stations are continuously recorded using a preset sampling frequency to obtain the three-dimensional coordinate point data of the GNSS reference stations, so that data correction can be performed using the three-dimensional coordinate points in subsequent steps.

[0037] For example, the sampling rate of the GNSS reference station is 1 Hz, and the correction uses "RTK technology correction".

[0038] S103, using the UAV LiDAR module, scans the target area along a preset flight path to generate digital surface model data with a point cloud density greater than or equal to a preset density threshold.

[0039] Specifically, the flight path of the UAV is planned for the target area, so that the UAV will fly in the target area along the preset flight path. Furthermore, the UAV LiDAR module on the UAV performs laser radar scanning, and collects the terrain of the target area according to the set point cloud collection purpose. Then, a digital surface model is generated based on the point cloud data. This digital surface model can reflect the terrain features of the target area and is used for terrain correction and local deformation verification.

[0040] For example, the point cloud density of the drone LiDAR is ≥50 points / m².

[0041] S104, The meteorological data of the target area within a preset historical time period is collected through the meteorological data acquisition module.

[0042] Specifically, through a meteorological data acquisition model, environmental parameters such as rainfall and temperature are continuously recorded. The acquisition frequency can be set as needed, for example, automatic acquisition at 1-hour intervals. In this way, meteorological data within a preset time period can be continuously recorded. This time period can be 3, 6, 9, or 12, with a quarter as a unit. Within a unit, meteorological environmental data is continuously collected from the target area. This meteorological data is used to analyze the impact of seasonal environmental factors on surface deformation and to provide environmental background parameters for the analysis of the causes of geological disasters.

[0043] S105, combine the multi-period SAR data, the three-dimensional coordinate point data, the digital surface model data, and the meteorological data to obtain multi-source data of the target area.

[0044] Specifically, in the aforementioned steps S101 to S104, the multi-period SAR data, the three-dimensional coordinate point data, the digital surface model data, and the meteorological data are obtained respectively. Next, after time alignment, they are combined to finally form a spatiotemporally correlated multi-source dataset, which is stored in the storage system of the preprocessing device to provide a data foundation for subsequent analysis.

[0045] S200: Combine the multi-period SAR data from the multi-source data to generate several interferometric pairs. Use the DEM data in the digital elevation model to sequentially eliminate the flat phase and terrain phase to generate a differential interferogram. Then, use the branch-cut phase unwrapping algorithm to solve the wrapped phase in the differential interferogram. Finally, use the PS-InSAR and SBAS-InSAR dual algorithms to perform time series calculations to obtain the annual deformation rate map.

[0046] Specifically, obtaining an annual deformation rate map may include the following steps: S201, using the InSAR time series calculation algorithm, the multi-period SAR data are combined and paired according to a preset time baseline threshold to generate several interferometric pairs.

[0047] Specifically, to eliminate potential random errors (such as atmospheric effects) in a single interferometric pair, an observation network containing a large number of interferometric pairs needs to be constructed using multiple images. The quality of this observation network is controlled by a "baseline threshold".

[0048] It should be noted that, using the InSAR time-series solution algorithm, based on a preset time baseline threshold (e.g., no more than 30 days), multiple SAR images are systematically combined and paired (i.e., combined according to preset time and spatial baseline thresholds). The first image is used as the geometric reference, and subsequent images are registered with it, for example, generating no fewer than 20 interferometric pairs. When combining interferometric pairs, the length of the spatial baseline must be controlled within a preset threshold to minimize the decorrelation effect caused by geometric differences, thereby ensuring the quality of the interferometric phase data and laying the foundation for high-precision deformation inversion.

[0049] S202 uses DEM data from the digital elevation model to sequentially perform flat-ground phase effect elimination and terrain phase effect elimination processing on several interferometric pairs to generate differential interferometric atlases.

[0050] Specifically, using known digital elevation models (DEMs), the phase contributions caused by surface undulations (topographic phase) and Earth curvature (flat phase) are simulated and subtracted.

[0051] For example, using external digital elevation model data (such as 30-meter resolution ASTER GDEM), non-deformation components in the interferometric phase are eliminated in two steps: first, the flat-area effect phase caused by the Earth's curvature and satellite orbital geometry is removed (subtracting large-scale systematic phase fringes caused by the reference ellipsoid); second, the terrain phase caused by topographic undulations is eliminated (subtracting the phase caused by actual topographic undulations, thus obtaining a differential interferogram mainly containing deformation, atmospheric, and noise components). The final generated clean differential interferogram set, in which each phase fringe represents a segment of deformation displacement in the radar line-of-sight direction.

[0052] S203, the differential interferogram is unwrapped using the branch-cut phase unwrapping algorithm to obtain the unwrapped phase result.

[0053] Specifically, because radar waves are periodic, with phase cycling in 2π (360 degrees) cycles, it can measure minute phase changes within a single cycle, but cannot directly distinguish whether this is a change within a single cycle or a cumulative change spanning multiple cycles. Therefore, phase entanglement occurs, and step S203 aims to untangle it and restore the complete number of cycles. For example, using the branching method, integration is performed along a continuous spatial path to restore the entangled relative phase value to an absolute, continuous deformation value, i.e., the untangled phase result. Furthermore, a "closed-loop" verification is used to ensure the correctness of the untangling process, i.e., integrating the solution along a closed path with an error controlled within 0.5 radians, ensuring the reliability of the result, thus obtaining the untangled phase result for each pixel. This untangled phase result is the direct input for calculating the deformation along the radar line-of-sight direction, and its physical meaning represents the precise distance change of the Earth's surface relative to the satellite.

[0054] S204. The PS-InSAR algorithm and SBAS-InSAR algorithm are run respectively to perform time series analysis on the unwrapped phase results to obtain the annual deformation rate map.

[0055] The PS-InSAR algorithm extracts the deformation of permanent scatterers in the unwrapped phase results, and the SBAS-InSAR algorithm performs deformation analysis on distributed scatterers in the unwrapped phase results.

[0056] An annual deformation rate map is a specialized map that can include different colors and uses a real geographic map (such as a satellite image or topographic map) as a base map. It can intuitively display the spatial distribution and intensity of slow surface deformation and is a routine and fundamental result map in the field of InSAR surface deformation monitoring.

[0057] Specifically, the deformation rate is extracted from the unwrapped phase sequence to obtain the annual deformation rate map.

[0058] For example, the PS-InSAR algorithm and the SBAS-InSAR algorithm are used to deal with different surface types: PS-InSAR (Permanent Scatterer Technology) identifies and tracks point targets, such as building corners and exposed rocks, that can stably reflect radar signals over long periods of time, and obtains their high-precision deformation time series.

[0059] SBAS-InSAR (Short Baseline Set Technology) processes distributed targets with unstable scattering characteristics, such as vegetated areas and farmland, by combining a large number of short spatiotemporal baseline interferometric pairs to extract their average deformation rate.

[0060] The results of the two algorithms are fused to generate a spatially continuous "annual deformation rate map". The value of each pixel in the map represents the deformation rate of the land surface at that point in one year, in millimeters per year.

[0061] It should be noted that PS-InSAR results provide high-precision point-like deformation information (such as building corners, exposed rocks, and other permanent scatterers). These points are extremely accurate, like "reference points" on the ground, but their spatial distribution is discontinuous. SBAS-InSAR results, on the other hand, provide spatially continuous planar deformation information, but this targets distributed scatterers such as vegetation and soil. The absolute accuracy of these points may be slightly lower than that of PS-InSAR points. Therefore, fusion is necessary. The fusion method treats PS-InSAR points as control points, while SBAS-InSAR results provide a continuous background field. The fusion algorithm uses PS-InSAR points as control points to constrain and calibrate the continuous deformation field generated by the SBAS-InSAR results, ultimately producing a deformation rate map that maintains spatial continuity and has high accuracy at key points.

[0062] S300, the three-dimensional coordinate point data, the digital surface model data, and the annual deformation rate map are spatially aligned, and then the fused deformation rate is calculated using a weighted formula based on the weights assigned to the three-dimensional coordinate point data, the digital surface model data, and the annual deformation rate map, to obtain the deformation rate field.

[0063] Specifically, obtaining the deformation rate field may include the following steps: S301, using transverse Mercator projection, spatially register and align the three-dimensional coordinate point data, the digital surface model data, and the annual deformation rate map to obtain a data stack with unified spatial reference.

[0064] Specifically, the aforementioned steps yielded data from satellites, ground base stations, and UAVs. These data sources have different initial coordinate systems, necessitating coordinate system unification and transformation to a single planar coordinate system before subsequent numerical calculations can be performed. The UTM (Unified Transverse Mercator) projection is a globally standardized zoned planar coordinate system that effectively controls projection distortion and ensures computational accuracy.

[0065] For example, the geodetic coordinates (latitude and longitude) obtained from GNSS measurements and the orbital coordinate system of InSAR data are converted into plane rectangular coordinates (X,Y) under a specified UTM projection zone through a mathematical model. Here, the mathematical model refers to the existing coordinate transformation formula based on ellipsoidal geodesy in map projection theory; this application has not adjusted or improved this formula, but essentially it systematically converts spherical coordinates (latitude and longitude) into plane coordinates. Next, the digital surface model data generated by LiDAR and the annual deformation rate map from InSAR are geometrically corrected (geometric correction refers to establishing a coordinate mapping relationship through corresponding feature points on the two images, thereby resampling one image to make it spatially perfectly aligned with the other), so that its grid pixels correspond to the geographical locations under the UTM coordinate system. Furthermore, through control point matching and interpolation algorithms, the registration error can be controlled within 1 meter; thus, a data stack with a completely unified spatial reference is obtained, where each pixel position corresponds to three types of deformation data from GNSS, LiDAR, and InSAR.

[0066] For example, the registration error can be controlled within 1 meter through control point matching and interpolation algorithms. This can be achieved by the following method: To achieve the accuracy requirement of a spatial registration error of no more than 1 meter for multi-source data, firstly, control points with clear positioning characteristics are extracted from the LiDAR digital surface model and the InSAR annual deformation rate map through a combination of manual identification and matching, and point correlation is performed using feature matching algorithms such as SIFT; then, a polynomial transformation model is established based on the successfully matched control point set, the model parameters are solved by the least squares method, and the residuals of each control point are calculated. After iteratively eliminating abnormal points in the residuals, the model is optimized; then, the coordinates of all pixels are resampled using a bilinear interpolation algorithm to transfer the high-precision positioning characteristics of the control points to the entire image area; finally, the registration results are verified using a checkpoint set independent of the modeling process to ensure that the root mean square error is strictly controlled within 1 meter.

[0067] S302, assign weights to the three-dimensional coordinate point data, the digital surface model data, and the annual deformation rate map in the unified data stack of the spatial reference, respectively, and calculate the fused deformation rate using a weighted formula to obtain the deformation rate field. The calculation formula for the deformation rate field is as follows: ; Where V represents the deformation rate value of the deformation rate field. This represents the deformation rate value derived from InSAR technology. This represents the deformation value from GNSS measurements. This represents the terrain deformation value extracted from LiDAR data. , and These represent the weight coefficients of their respective data sources.

[0068] Specifically, by assigning appropriate weights to different data sources, the optimal deformation estimate is calculated. For example, InSAR has a wide coverage but may be affected by the atmosphere, GNSS has high accuracy but only represents discrete points, and LiDAR has rich detail but only reflects local terrain changes. Weighted fusion aims to trust the more reliable data source and allow it to contribute more to the final result.

[0069] In the formula, a specific weight coefficient is assigned to each data point to reflect its relative reliability. For example... The weight is 0.6. The weight is 0.3. The weight is 0.1. The allocation of weights can be set based on experience, and this application does not impose any restrictions.

[0070] By performing the above calculations on all pixels in the entire image, a spatially continuous fused deformation rate field is finally generated after multi-source data correction and optimization.

[0071] S400, the deformation rate field and the slope and aspect parameters are input into the pre-trained decision tree classification model to obtain the geomechanical deformation type of the target area.

[0072] The S400 step is the core of intelligent diagnosis of geological disasters. It aims to use machine learning algorithms to comprehensively analyze multiple parameters and achieve intelligent identification from quantitative deformation data to qualitative disaster types. Its technical principle is to use a pre-trained decision tree classification model to simulate the comprehensive analysis logic of geological experts on parameters such as deformation rate and terrain features, so as to achieve automated classification of geological disaster types such as landslides.

[0073] Specifically, obtaining the geomechanical deformation type of the target area may include the following steps: S401, the slope parameters and aspect parameters are obtained through data from the digital elevation model.

[0074] Specifically, slope and aspect are key topographic factors controlling slope stability. Slope directly affects the stability of soil and rock under gravity, while aspect indirectly affects the mechanical properties of soil and rock by influencing sunlight, rainfall, and surface water runoff.

[0075] Furthermore, a Digital Elevation Model (DEM) is used as the base data source, typically obtained from satellite remote sensing or aerial photogrammetry, containing elevation information for each pixel. For example, based on the DEM data, slope values ​​(in degrees) and aspect values ​​(in degrees, 0° being true north) are calculated element-by-element using a spatial difference algorithm. Slopes greater than 25° are identified as high-risk steep slope areas, and north-facing slopes (e.g., 337.5°–22.5°) are identified as areas prone to water accumulation. Finally, slope and aspect raster maps are generated that perfectly match the spatial resolution and range of the deformation rate field, representing the slope and aspect parameters.

[0076] S402, the deformation rate field, the slope parameter and the aspect parameter are respectively input into the pre-trained decision tree classification model to obtain the geomechanical deformation type of the target area, wherein the pre-trained decision tree classification model represents a pre-trained decision tree classification model trained based on sample data of various historical geological disaster types.

[0077] Specifically, the decision tree classification model refers to classifying input data through hierarchical judgment rules. This decision tree classification model learns from a large number of historical samples of known disaster types to uncover the inherent correlation rules between different disaster types and parameters such as deformation rate, slope, and aspect.

[0078] For example, a decision tree classification model receives three features for each pixel: 1) Deformation rate value (V) in the deformation rate field; 2) Slope value (S) in the slope diagram; 3) The aspect value (A) in the aspect diagram.

[0079] Next, the decision tree classification model performs threshold judgments on the feature values ​​layer by layer starting from the root node. For example, it first judges whether the slope is greater than 30 degrees. If so, it proceeds to the next layer to judge whether the deformation rate is greater than 10 mm / year. Then it judges whether the slope direction is in the north direction. If all conditions are met, it is finally determined to be the disaster type corresponding to the end point of the path.

[0080] It should be noted that the pre-trained decision tree classification model was trained on a historical sample dataset containing 10 common landslide geomechanical deformation types. The sample set contains the true values ​​of features such as deformation rate, slope, and aspect for each type, along with their type labels. In the output, the decision tree classification model assigns a discrete classification label, i.e., the geomechanical deformation type, to each pixel. For example, when a pixel meets the following rules: "slope > 30°, deformation rate > 10 mm / year, and northward aspect," the pre-trained decision tree classification model will determine and output that pixel as a "soil-traction landslide."

[0081] S500: Using a target LSTM neural network, the deformation rate of the target region within a preset time period is predicted, thus obtaining the deformation rate prediction of the target region within the preset time period.

[0082] Specifically, obtaining the deformation rate prediction of the target region within a preset future time period may include the following steps: S501, Based on the historical multi-source data within a preset time period of the target region, train the LSTM neural network to obtain the target LSTM neural network.

[0083] S502, using the target LSTM neural network, the deformation rate of the target region within a preset time period is predicted to obtain the deformation rate prediction value.

[0084] Specifically, the historical multi-source data includes multiple periods of SAR data collected by SAR satellites in the same orbit in the target area, three-dimensional coordinate point data of the reference station obtained by GNSS base stations, digital surface model data generated by UAV scanning, and meteorological data for a preset time period. These historical multi-source raw deformation time-series data are used to train the initial LSTM neural network.

[0085] For example, for the initial training of the LSTM neural network: Training data includes time-aligned deformation rate time series acquired from InSAR, GNSS, LiDAR, and other technologies, as well as auxiliary features such as meteorological data. Each sample represents data from a time window (e.g., 30 time points).

[0086] Training objective: To learn the mapping relationship from historical sequences to deformation rates at one or more future time points.

[0087] Training process: includes data normalization, setting the number of iterations (500 times), learning rate (0.001), and using mean squared error (MSE) as the loss function for optimization.

[0088] Furthermore, the total number of training samples is no less than 1000 sets, and each set of samples contains deformation rate observations at 30 consecutive time points and their corresponding true deformation rate values ​​at future time points. The deformation rate data is normalized to convert all values ​​to the [0,1] interval to eliminate the influence of dimensions and accelerate model convergence. Next, the initial LSTM neural network is trained using historical multi-source raw deformation time series data. The training method is conventional and is not limited in this application. After training, the target LSTM neural network is obtained.

[0089] Subsequently, a target LSTM neural network is used to predict the deformation rate of the target region within a preset future time period, obtaining the predicted deformation rate value. For example, deformation rate data from the most recent 30 time points of the target region are extracted to form a complete time series window. The same standardized parameters as in the training phase are used to normalize the data to be predicted, ensuring that the input data has no missing or outlier values ​​and that the intervals between consecutive time points are uniform. The preprocessed time series data is input into the trained target LSTM neural network. Through linear transformations and nonlinear activation functions in each layer of the network, temporal features are extracted step by step. A rolling prediction method is used, first predicting the deformation rate at the first future time point, then adding this predicted value to the input sequence, and continuously predicting subsequent time points to generate the predicted deformation rate value for each time point within the next 6 months, forming a complete predicted time series curve.

[0090] S600, based on the predicted deformation rate of the target area and the comparison between the geomechanical deformation type of the target area and the preset warning threshold, the warning level of the target area within a preset time period is obtained.

[0091] Specifically, obtaining the warning level of the target area within a preset future time period may include the following steps: S601, based on the geomechanical deformation type of the target area, call the corresponding graded early warning threshold from the preset threshold rule library to obtain the graded early warning threshold range matching the geomechanical deformation type of the target area.

[0092] Among them, different geomechanical deformation types in the preset threshold rule base correspond to different threshold ranges.

[0093] Specifically, the preset threshold rule base generates specific early warning levels based on deformation prediction results and disaster type identification results through differentiated graded threshold judgments; the technology lies in setting exclusive graded early warning threshold ranges for each type of geomechanical deformation in the preset threshold rule base, thereby achieving precision and differentiation of early warning standards.

[0094] Furthermore, the construction principle of the differentiated threshold rule base: Type-specific threshold configuration: The preset threshold rule library establishes independent graded early warning threshold ranges for each type of geomechanical deformation.

[0095] For example, soil-induced landslides: yellow alert [5,10) mm / year, orange alert [10,15) mm / year, red alert [15,+∞) mm / year; Rockfall: Yellow alert [8,12) mm / year, Orange alert [12,18) mm / year, Red alert [18,+∞) mm / year; Push-type landslides: Yellow warning [6,9) mm / year, Orange warning [9,14) mm / year, Red warning [14,+∞) mm / year.

[0096] Threshold determination criteria: Threshold ranges for each type of disaster are set based on historical disaster-causing patterns, soil and rock mechanical properties, and movement characteristics, reflecting their unique instability critical conditions. These threshold ranges can be determined by comprehensively referencing industry standards, historical case statistics, geomechanical model inversion of the study area, and expert experience, thereby ensuring their scientific validity and regional applicability.

[0097] It should be noted that the threshold matching process includes receiving the geomechanical deformation type of the target area, and, based on the identified specific disaster type, calling the corresponding exclusive graded early warning threshold range from the rule base.

[0098] S602, the deformation rate prediction of the target area is matched and judged with the graded early warning threshold range to obtain the early warning level of the target area within a preset time period in the future.

[0099] Specifically, this step generates a targeted warning level by matching the predicted deformation rate of the target area with a specific graded warning threshold range.

[0100] In step S602, the categorized threshold comparison compares the predicted deformation rate with the threshold range specific to the disaster type to determine which threshold range the predicted deformation rate falls into, thereby determining the warning level; and finally generates the warning level for the target area within a preset time period in the future.

[0101] For example, the final output could be a warning level distribution map, which spatially indicates the warning level (e.g., no warning, yellow warning, orange warning, red warning) of different regions in a specific future period (e.g., the next 6 months), providing intuitive decision support for risk management.

[0102] S700, the warning level is matched with the preset warning level execution plan, and a report is generated to obtain the warning report of the execution plan.

[0103] Specifically, obtaining an early warning report for the execution plan may include the following steps: S701, build a preset scheme library.

[0104] The preset scheme library includes execution schemes corresponding to each warning level.

[0105] Specifically, the preset solution library can adopt a three-level classification structure, which corresponds completely to the warning levels: Yellow Alert Solution Library: For initial risks, it includes basic response solutions such as increasing monitoring frequency (e.g., once a week), inspecting key areas, and conducting preliminary risk assessments; Orange Alert Solution Library: For moderate risks, it includes enhanced response plans such as restricting activities in dangerous areas, deploying professional monitoring equipment, and activating emergency response plans. Red Alert Solution Library: Targeting high-risk situations, it includes emergency response plans such as emergency evacuation, engineering rescue, and 24 / 7 monitoring.

[0106] It should be noted that the standardization of the solution content includes: Spatial measures: Each plan includes specific spatial control requirements, such as the delineation of the warning area, the planning of evacuation routes, and the location of monitoring points; Timeframe: Clearly define the implementation timetable and response deadlines for each measure; Responsible entities: Designate the implementing departments and cooperating units for each measure; Resource allocation: Clarify the configuration plan for emergency supplies, professional personnel, and equipment.

[0107] S702, based on the warning level, a matching is performed in the preset scheme library to obtain the execution scheme for the target area.

[0108] Specifically, the execution plan for the target area includes: Level identification: Receives the warning level signal (yellow / orange / red) output from the threshold judgment unit; Solution retrieval: Automatically retrieves the corresponding standardized response plan template from the preset solution library based on the warning level; Parameter binding: Real-time parameters such as the coordinates of specific potential hazards, the scope of impact, and the risk level are injected into the solution template to generate a targeted execution plan.

[0109] The above method yields the execution plan for the target area.

[0110] S703, Based on the execution plan for the target area, generate a visualized map product containing the spatial distribution of the warning area and form a report to obtain the warning report of the execution plan.

[0111] Specifically, the generation of spatial visualization products includes: Heat map generation: Based on the GIS platform, the deformation rate field is rendered into a gradient color map to intuitively display the spatial distribution of risk; Warning area marking: Accurately mark the impact range of different warning levels on the base map, with boundary accuracy controlled within 50 meters; 3D situation display: Combined with digital elevation model, a 3D terrain deformation model is generated to show the development characteristics of geological disasters in a three-dimensional way.

[0112] The final early warning report for the implementation plan is obtained. This report includes integrated early warning levels, spatial distribution maps, and response plans, generated in a standardized PDF format. The report may include deformation data charts, on-site photos, historical comparison data, and other decision-making support information. Furthermore, the report is automatically pushed to relevant responsible personnel via a 4G module, simultaneously stored on a local server, and uploaded to the provincial monitoring platform.

[0113] Specifically, emergency notification triggers may include: Tiered notification mechanism: Different notification plans are activated according to the warning level, with a response mechanism such as within 10 minutes implemented for red warnings; Multiple channels for dissemination: Warning information is simultaneously released through multiple channels such as SMS, government platforms, and emergency broadcasts.

[0114] In one embodiment, after step S100, data preprocessing of multi-period SAR data is further included, specifically including: The format conversion unit of the data preprocessing device converts multi-period SAR data from CEOS format to ENVI format, preserving the original phase and amplitude information. The radiation correction unit inputs atmospheric profile data (obtained from the NASA website) and uses a Gamma correction model to eliminate the effects of atmospheric scattering. The geometric registration unit uses the first phase of SAR data as a reference and matches subsequent data based on the SIFT algorithm. Through iterative optimization, the registration error is made ≤1 / 8 pixel. The noise filtering unit uses a 5×5 window Goldstein (Goldstein refers to an adaptive filtering algorithm widely used in InSAR data processing, namely the Goldstein phase filtering algorithm) filter to smooth the noise in the interferogram.

[0115] In one embodiment, based on basic data of mining subsidence areas, collapse areas, and high-risk and key prevention areas for geological disasters such as debris flows and landslides in a certain city, surface deformation monitoring of 128 ground collapse points in the urban area was first carried out using InSAR technology based on Sentinel1-1 radar data. At the same time, internal consistency verification was carried out based on LuTan-1 radar data. Three deformation anomaly points were selected in each of the several counties under the jurisdiction of the city, and two deformation anomaly points were selected in the urban area of ​​the city for field investigation, and two external consistency verifications of surface deformation based on GNSS measurements were carried out.

[0116] All data were analyzed to establish a remote sensing monitoring database for geological disasters and a remote sensing monitoring database for the ecological environment of key mining areas; an integrated air-space-ground monitoring and early warning system platform for mine ecology and geological disasters in a certain city was built, providing multiple functions such as querying, managing and serving monitoring results.

[0117] In summary, this application has the following beneficial effects: Improving monitoring accuracy: By utilizing multi-source remote sensing, machine learning, and other technologies, and simultaneously inverting the overall surface deformation using both SBAS-InSAR and PS-InSAR time-series InSAR techniques, the accuracy of the results can be compared. This allows for more accurate acquisition of surface deformation information and improves the accuracy of automated disaster identification. For example, in the deformation monitoring of ground subsidence points in a certain city, this method was used to construct an accurate geological hazard catalog, laying the foundation for subsequent model training and testing. Compared with traditional methods, the monitoring accuracy is significantly improved.

[0118] Reduced monitoring costs: InSAR technology has the advantage of wide coverage. Compared with traditional geological disaster monitoring methods, such as manual inspections, it can monitor over a larger area, reducing the investment of manpower and resources, thereby reducing the cost of geological disaster monitoring and early warning.

[0119] Improving monitoring efficiency: InSAR-based methods can achieve parallel monitoring and analysis of multiple sub-regions through grid partitioning and automated data processing, greatly improving monitoring efficiency.

[0120] Enhancing early warning accuracy: By selecting appropriate susceptibility evaluation factors and constructing geological disaster early warning models using various machine learning algorithms, and evaluating the models using indicators such as accuracy, ROC curves, and AUC values ​​to determine the optimal model, the deviation between geological disaster early warning areas and early warning effects can be minimized. Furthermore, combining airborne LiDAR technology for microscopic verification further enhances the accuracy of early warnings.

[0121] A tiered early warning system is implemented: the displacement of dangerous rock masses is divided into different warning levels, achieving a clear classification of the risk of dangerous rock mass displacement. This classification method helps relevant departments and personnel quickly identify the current risk level, thereby taking corresponding countermeasures and improving the efficiency and effectiveness of emergency response.

[0122] Secondly, this application provides a surface monitoring and geological disaster early warning system based on multi-phase InSAR remote sensing data, applied to the aforementioned method for surface monitoring and geological disaster early warning based on multi-phase InSAR remote sensing data, including: The acquisition unit is used to acquire SAR data from multiple periods of the same orbit collected by SAR satellites in the target area, three-dimensional coordinate point data of the reference station acquired by GNSS base station, digital surface model data generated by UAV scanning, and meteorological data for a preset time period, so as to obtain multi-source data of the target area. The deformation rate unit is used to combine the multi-period SAR data in the multi-source data to generate several interferometric pairs, use the DEM data in the digital elevation model to sequentially eliminate the flat phase and the terrain phase, generate a differential interferogram, and then use the branch-cut phase unwrapping algorithm to solve the wrapped phase in the differential interferogram. Finally, the PS-InSAR and SBAS-InSAR dual algorithms are used to perform time series calculations to obtain the annual deformation rate map. The rate field unit is used to spatially align the three-dimensional coordinate point data, the digital surface model data, and the annual deformation rate map. Then, based on the weights assigned to the three-dimensional coordinate point data, the digital surface model data, and the annual deformation rate map, a weighted formula is used to calculate the fused deformation rate to obtain the deformation rate field. The decision unit is used to input the deformation rate field and the slope and aspect parameters into a pre-trained decision tree classification model to obtain the geomechanical deformation type of the target area. The prediction unit is used to predict the deformation rate of the target region within a preset time period using a target LSTM neural network, thereby obtaining the deformation rate prediction of the target region within the preset time period. The threshold judgment unit is used to compare the predicted deformation rate of the target area and the geomechanical deformation type of the target area with a preset warning threshold to obtain the warning level of the target area within a preset time period in the future. The result unit is used to match the warning level with the preset warning level execution plan and generate a report to obtain the warning report of the execution plan.

[0123] Thirdly, this application provides a computer program that, when executed by a processor, implements the steps of the aforementioned method for surface monitoring and geological disaster early warning of multi-phase InSAR remote sensing data.

[0124] It should be noted that the steps shown in the flowchart in the accompanying drawings can be executed in a computer system such as a set of computer-executable instructions, and although a logical order is shown in the flowchart, in some cases the steps shown or described may be executed in a different order than that shown here.

[0125] Obviously, those skilled in the art should understand that the various units or steps of this application described above can be implemented using general-purpose computing devices. They can be centralized on a single computing device or distributed across a network of multiple computing devices. Optionally, they can be implemented using computer-executable program code, thereby storing them in a storage device for execution by a computing device, or fabricating them separately as individual integrated circuit modules, or fabricating multiple modules or steps into a single integrated circuit module. Thus, this application is not limited to any particular combination of hardware and software.

[0126] The above description is merely a preferred embodiment of this application and is not intended to limit this application. Various modifications and variations can be made to this application by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this application should be included within the protection scope of this application.

Claims

1. A method for surface monitoring and geological disaster early warning using multi-phase InSAR remote sensing data, characterized in that, include: The target area is obtained by acquiring multiple SAR data from the same orbit collected by SAR satellites, three-dimensional coordinate point data of the reference station obtained by GNSS base station, digital surface model data generated by UAV scanning, and meteorological data for a preset time period. The multi-period SAR data from the multi-source data are combined to generate several interferometric pairs. The flat phase and terrain phase are eliminated sequentially using the DEM data in the digital elevation model to generate a differential interferogram. Then, the branch-cut phase unwrapping algorithm is used to solve the wrapped phase in the differential interferogram. Finally, the PS-InSAR and SBAS-InSAR dual algorithms are used to perform time series calculations to obtain the annual deformation rate map. Spatially align the three-dimensional coordinate point data, the digital surface model data, and the annual deformation rate map. Then, based on the weights assigned to the three-dimensional coordinate point data, the digital surface model data, and the annual deformation rate map respectively, calculate the fused deformation rate using a weighted formula to obtain the deformation rate field. The deformation rate field, slope and aspect parameters are input into a pre-trained decision tree classification model to obtain the geomechanical deformation type of the target area. Using a target LSTM neural network, the deformation rate of the target region within a preset future time period is predicted, thus obtaining the deformation rate prediction of the target region within the preset future time period. Based on the predicted deformation rate of the target area and the comparison between the geomechanical deformation type of the target area and the preset warning threshold, the warning level of the target area within a preset time period in the future is obtained; The warning level is matched with the preset warning level execution plan, and a report is generated to obtain the warning report of the execution plan.

2. The method for surface monitoring and geological disaster early warning using multi-phase InSAR remote sensing data as described in claim 1, characterized in that, The steps of obtaining multi-source data for the target area by acquiring multiple periods of SAR data from the same orbit collected by SAR satellites, three-dimensional coordinate point data of the reference station acquired by GNSS base stations, digital surface model data generated by UAV scanning, and meteorological data for a preset time period include: The SAR satellite receiving module collects multiple SAR images of the same orbit at a preset revisit period to obtain the multi-phase SAR data. The three-dimensional coordinate point data of the base station are obtained through the GNSS base station module at a preset sampling rate; The target area is scanned along a preset flight path using the drone's LiDAR module, generating digital surface model data with a point cloud density greater than or equal to a preset density threshold. Meteorological data of the target area within a preset historical time period is collected through the meteorological data acquisition module; The multi-period SAR data, the three-dimensional coordinate point data, the digital surface model data, and the meteorological data are combined to obtain multi-source data for the target area.

3. The method for surface monitoring and geological disaster early warning using multi-phase InSAR remote sensing data as described in claim 1, characterized in that, The steps of combining multi-period SAR data from the multi-source data to generate several interferometric pairs, sequentially eliminating flat and topographic phases using DEM data from the digital elevation model to generate differential interferograms, then using a branch-cut phase unwrapping algorithm to solve the wrapped phases in the differential interferograms, and finally using a dual algorithm of PS-InSAR and SBAS-InSAR to perform time-series calculations to obtain the annual deformation rate map include: Using the InSAR time series solution algorithm, the multi-period SAR data are combined and paired according to a preset time baseline threshold to generate several interferometric pairs; Using DEM data from the digital elevation model, several interferometric pairs are sequentially processed to eliminate flat-ground phase effects and terrain phase effects, generating differential interferometric atlases. The phase unwrapping calculation of the differential interferogram is performed using the branch-cut phase unwrapping algorithm to obtain the unwrapped phase result; The PS-InSAR algorithm and the SBAS-InSAR algorithm are run respectively to perform time series analysis on the unwrapped phase results to obtain an annual deformation rate map. The PS-InSAR algorithm extracts the deformation of permanent scatterers in the unwrapped phase results, and the SBAS-InSAR algorithm performs deformation analysis on distributed scatterers in the unwrapped phase results.

4. The method for surface monitoring and geological disaster early warning using multi-phase InSAR remote sensing data as described in claim 1, characterized in that, The steps of spatially aligning the three-dimensional coordinate point data, the digital surface model data, and the annual deformation rate map, and then calculating the fused deformation rate using a weighted formula based on the weights assigned to the three-dimensional coordinate point data, the digital surface model data, and the annual deformation rate map to obtain the deformation rate field include: Using transverse Mercator projection, the three-dimensional coordinate point data, the digital surface model data and the annual deformation rate map are spatially registered and aligned to obtain a data stack with unified spatial reference. Weights are assigned to the three-dimensional coordinate point data, the digital surface model data, and the annual deformation rate map in the unified data stack of the spatial reference, respectively. The fused deformation rate is then calculated using a weighted formula to obtain the deformation rate field. The formula for calculating the deformation rate field is as follows: ； Where V represents the deformation rate value of the deformation rate field. This represents the deformation rate value derived from InSAR technology. This represents the deformation value from GNSS measurements. This represents the terrain deformation value extracted from LiDAR data. , and These represent the weight coefficients of their respective data sources.

5. The method for surface monitoring and geological disaster early warning using multi-phase InSAR remote sensing data as described in claim 1, characterized in that, The step of inputting the deformation rate field, slope, and aspect parameters into a pre-trained decision tree classification model to obtain the geomechanical deformation type of the target area includes: The slope parameters and aspect parameters are obtained using data from the digital elevation model. The deformation rate field, the slope parameter, and the aspect parameter are respectively input into a pre-trained decision tree classification model to obtain the geomechanical deformation type of the target area. The pre-trained decision tree classification model refers to a pre-trained decision tree classification model trained based on sample data of various historical geological disaster types.

6. The method for surface monitoring and geological disaster early warning using multi-phase InSAR remote sensing data as described in claim 1, characterized in that, The step of using a target LSTM neural network to predict the deformation rate of the target region within a preset future time period, and obtaining the predicted deformation rate of the target region within the preset future time period, includes: Based on historical multi-source data within a preset time period of the target region, an LSTM neural network is trained to obtain the target LSTM neural network; Using the target LSTM neural network, the deformation rate of the target region within a preset time period is predicted to obtain the deformation rate prediction value.

7. The method for surface monitoring and geological disaster early warning using multi-phase InSAR remote sensing data as described in claim 1, characterized in that, The step of comparing the predicted deformation rate of the target area and the geomechanical deformation type of the target area with a preset warning threshold to obtain the warning level of the target area within a preset future time period includes: Based on the geomechanical deformation type of the target area, the corresponding graded early warning threshold is called from the preset threshold rule library to obtain the graded early warning threshold range matching the geomechanical deformation type of the target area. Different geomechanical deformation types in the preset threshold rule library correspond to different threshold ranges. The deformation rate prediction of the target area is matched with the graded early warning threshold range to obtain the early warning level of the target area within a preset time period in the future.

8. The method for surface monitoring and geological disaster early warning using multi-phase InSAR remote sensing data as described in claim 1, characterized in that, The step of matching the warning level with a preset warning level execution plan and generating a report to obtain the warning report of the execution plan includes: Construct a preset solution library, wherein the preset solution library includes an execution solution corresponding to each warning level; Based on the warning level, a matching scheme is performed in the preset scheme library to obtain the execution scheme for the target area; Based on the execution plan for the target area, a visual map product containing the spatial distribution of the warning area is generated, and a report is generated to obtain the warning report for the execution plan.

9. A surface monitoring and geological disaster early warning system based on multi-phase InSAR remote sensing data, characterized in that, A method for surface monitoring and geological disaster early warning using multi-phase InSAR remote sensing data as described in any one of claims 1 to 8, comprising: The acquisition unit is used to acquire SAR data from multiple periods of the same orbit collected by SAR satellites in the target area, three-dimensional coordinate point data of the reference station acquired by GNSS base station, digital surface model data generated by UAV scanning, and meteorological data for a preset time period, so as to obtain multi-source data of the target area. The deformation rate unit is used to combine the multi-period SAR data in the multi-source data to generate several interferometric pairs, use the DEM data in the digital elevation model to sequentially eliminate the flat phase and the terrain phase, generate a differential interferogram, and then use the branch-cut phase unwrapping algorithm to solve the wrapped phase in the differential interferogram. Finally, the PS-InSAR and SBAS-InSAR dual algorithms are used to perform time series calculations to obtain the annual deformation rate map. The rate field unit is used to spatially align the three-dimensional coordinate point data, the digital surface model data, and the annual deformation rate map. Then, based on the weights assigned to the three-dimensional coordinate point data, the digital surface model data, and the annual deformation rate map, a weighted formula is used to calculate the fused deformation rate to obtain the deformation rate field. The decision unit is used to input the deformation rate field and the slope and aspect parameters into a pre-trained decision tree classification model to obtain the geomechanical deformation type of the target area. The prediction unit is used to predict the deformation rate of the target region within a preset time period using a target LSTM neural network, thereby obtaining the deformation rate prediction of the target region within the preset time period. The threshold judgment unit is used to compare the predicted deformation rate of the target area and the geomechanical deformation type of the target area with a preset warning threshold to obtain the warning level of the target area within a preset time period in the future. The result unit is used to match the warning level with the preset warning level execution plan and generate a report to obtain the warning report of the execution plan.

10. A computer program, characterized in that, When the computer program is executed by the processor, it implements the steps of the method for surface monitoring and geological disaster early warning of multi-phase InSAR remote sensing data as described in any one of claims 1 to 8.

Images