A method for intelligently monitoring geological disasters based on satellite remote sensing

Abstract

The present application relates to the technical field of geological monitoring, and particularly relates to a method for intelligently monitoring geological disasters based on satellite remote sensing, which comprises: obtaining corresponding ground cover information through image registration and orthographic correction processing; identifying terrain boundaries, geological structure features, and geological attribute information through the ground cover information; identifying surrounding hidden trouble points through a terrain feature extraction method; extracting stable scatterers and constructing a deformation monitoring network through an image fusion enhancement method in cooperation with a permanent scatterer interferometry method, and calculating phase difference values in a high-risk area; obtaining ground three-dimensional deformation field information in a monitoring area through phase difference value unwrapping processing; establishing a geological disaster risk assessment model through spatiotemporal clustering analysis; and quantitatively assessing geological disaster risks in the monitoring area through the established risk assessment model. The present application has the beneficial effect of realizing comprehensive monitoring and accurate early warning of geological disaster hidden dangers around a substation.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of geological monitoring, and particularly relates to a method for intelligently monitoring geological disasters based on satellite remote sensing. BACKGROUND

[0002] With the acceleration of urbanization, the environment around important infrastructures such as substations is increasingly complex, and the risk of geological disasters is increasingly prominent. Geological disasters such as landslides, collapses and ground subsidence not only threaten personal safety and property, but also may cause power supply interruption, seriously affecting social and economic stability and development. In recent years, the rapid development of remote sensing technology has provided a new means for geological disaster monitoring. The combination of high-resolution multispectral satellite images and synthetic aperture radar (SAR) images makes the extraction of surface cover change and geological features more accurate. Through the analysis of multi-temporal images, potential geological disaster hazard areas can be effectively identified, so as to conduct risk assessment and early warning. In this process, the fusion of visible light image intelligent analysis and geological disaster hazard risk assessment faces multiple technical challenges. After obtaining high-resolution multispectral satellite images, how to realize accurate surface cover information extraction through image registration and orthorectification, and accurately classify ground objects. When extracting surface texture information and identifying terrain boundaries, how to effectively combine edge detection and segmentation technology to improve the identification accuracy of potential geological disaster hazard points. Finally, when fusing optical images and synthetic aperture radar images for surface deformation monitoring, the interference processing accuracy of permanent scatterer interferometric radar needs to be improved to obtain more accurate three-dimensional deformation field information of the surface. The solution of these interrelated technical problems is crucial for accurate geological disaster risk assessment. SUMMARY

[0003] In view of the above or existing problems in the prior art, the present application is proposed.

[0004] Therefore, the purpose of the present application is to provide a method for intelligently monitoring geological disasters based on satellite remote sensing, which can accurately identify potential geological disaster hazard areas, conduct risk assessment and early warning.

[0005] To solve the above technical problems, the present application provides the following technical scheme: a method for intelligently monitoring geological disasters based on satellite remote sensing, comprising obtaining corresponding surface cover information through image registration and orthorectification processing;

[0006] Through the surface cover information, the terrain boundaries and geological structure features are identified;

[0007] Through the terrain boundaries and geological structure features, the geological attribute information of lithology and strata occurrence is extracted;

[0008] Through the terrain feature extraction method, the surrounding hazard points are identified;

[0009] Through the image fusion enhancement method, the permanent scatterer interference method is matched, the stable scatterer is extracted, and a deformation monitoring network is constructed, and the phase difference value in the high-risk area is calculated;

[0010] Through phase difference value unwrapping processing, combined with track parameters and atmospheric phase, the surface three-dimensional deformation field information in the monitoring area is inverted, and the time series analysis method is used to calculate the cumulative deformation amount;

[0011] Through the space-time clustering analysis, the geological disaster risk assessment model is established;

[0012] Through the establishment of the risk assessment model, the geological disaster risk in the monitoring area is quantitatively evaluated, and the risk grade is divided to generate early warning information.

[0013] As a preferred scheme of the method for intelligently monitoring geological disasters based on satellite remote sensing, the image registration and orthographic correction processing comprises: acquiring multi-time period multi-spectral image data and synthetic aperture radar image of the area around the substation by using a high-resolution multi-spectral satellite.

[0014] The multi-spectral image data comprises visible light, near-infrared and short-wave infrared bands.

[0015] The synthetic aperture radar image comprises radar echo signals of the area around the substation and corresponding radar image data.

[0016] As a preferred scheme of the method for intelligently monitoring geological disasters based on satellite remote sensing, the surface cover information comprises: extracting surface texture information by using a texture feature extraction method, and identifying terrain boundaries and geological structure features by using an edge detection and segmentation method.

[0017] The texture feature extraction method comprises calculating texture features by using a gray level co-occurrence matrix to obtain surface texture information, and the surface texture information is a surface texture feature map.

[0018] The edge detection and segmentation method comprises detecting an enhanced texture feature map to obtain terrain boundary line segments, and performing image segmentation on the edge detection result by using a watershed algorithm according to the terrain boundary line segments, setting a minimum region area threshold, merging regions smaller than the threshold, and obtaining a geological structure feature region.

[0019] As a preferred scheme of the method for intelligently monitoring geological disasters based on satellite remote sensing, the watershed algorithm comprises gradient calculation, minimum value marking, flood filling, catchment point meeting, post-processing, minimum region area threshold and morphological operation.

[0020] As a preferred embodiment of the method for intelligent monitoring of geological hazards based on satellite remote sensing of the present invention, the geological attribute information of lithology and stratigraphic attitude includes extracting shape index and directional features from the geological structural feature area, calculating the circularity and main direction angle of each segmented area, combining the orientation information of the topographic boundary line segment, identifying the lithology type and stratigraphic dip angle, determining the stratigraphic attitude, and obtaining geological attribute information.

[0021] As a preferred embodiment of the method for intelligent monitoring of geological disasters based on satellite remote sensing of the present invention, the terrain feature extraction method includes extracting terrain slope and aspect information from geological attribute information;

[0022] The surrounding potential hazards include high-risk areas prone to geological disasters such as landslides, collapses, and ground subsidence, identified through topographic slope and aspect information, forming a distribution map of potential hazards around the substation.

[0023] As a preferred embodiment of the method for intelligent monitoring of geological disasters based on satellite remote sensing of the present invention, the image fusion enhancement method includes: using discrete wavelet transform to fuse optical images and synthetic aperture radar images, taking the average value of low-frequency coefficients, selecting high-frequency coefficients using the absolute value maximum method, extracting high-frequency and low-frequency information through multi-scale decomposition, using different weight coefficients for fusion reconstruction to obtain a fused high-resolution image, and georegistering the fused high-resolution image with the hazard point distribution map;

[0024] The permanent scatterer interferometry method includes selecting points in high-risk areas with amplitude deviation index less than a preset threshold and coherence greater than a preset threshold as candidate permanent scatterers based on the registered fused image and the distribution map of potential hazard points. The candidate permanent scatterers are then screened, and the screened stable scatterers are subjected to time-series interferometry processing.

[0025] As a preferred embodiment of the satellite remote sensing intelligent monitoring method for geological disasters of the present invention, the phase difference unwrapping process includes: using a least squares phase expansion algorithm to unwrap the phase difference to obtain a continuous phase field; based on the continuous phase field, atmospheric phase correction is performed using a combination of tropospheric delay model and ionospheric delay model to obtain corrected phase information; the transformation matrix between line-of-sight deformation and three-dimensional deformation components is solved using the least squares method, converting the corrected phase information into a three-dimensional deformation field of the Earth's surface; the three-dimensional deformation field of the Earth's surface is then decomposed at multiple scales, the long-term trend term is extracted, and the cumulative deformation is calculated.

[0026] As a preferred embodiment of the satellite remote sensing intelligent monitoring method for geological disasters of the present invention, the spatiotemporal clustering analysis includes extracting the spatial distribution characteristics and temporal evolution patterns of deformation anomaly areas, combining historical disaster data statistics, meteorological factors, and geological attribute information.

[0027] As a preferred embodiment of the satellite remote sensing intelligent monitoring method for geological disasters of the present invention, the risk assessment model includes: quantitatively assessing the risks of landslides, collapses, and ground subsidence geological disasters within the monitoring area, classifying risk levels, and generating multiple types of hazard warning information, which include risk type, level, scope of impact, and recommended measures.

[0028] The beneficial effects of this invention are as follows: By combining high-resolution multispectral satellite imagery and SAR imagery, this invention can more accurately extract land cover changes and geological features, effectively identify potential geological hazard areas, and also help to conduct risk assessment and early warning of geological hazards, improve the monitoring capabilities of disasters such as landslides, collapses and ground subsidence, thereby protecting personal safety, property and power supply, and promoting the stability and development of the social economy. Attached Figure Description

[0029] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort. Wherein:

[0030] Figure 1 This is a flowchart illustrating a method for intelligent monitoring of geological hazards based on satellite remote sensing. Detailed Implementation

[0031] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings.

[0032] Many specific details are set forth in the following description in order to provide a full understanding of the invention. However, the invention may also be practiced in other ways different from those described herein, and those skilled in the art can make similar extensions without departing from the spirit of the invention. Therefore, the invention is not limited to the specific embodiments disclosed below.

[0033] Secondly, the term "one embodiment" or "embodiment" as used herein refers to a specific feature, structure, or characteristic that may be included in at least one implementation of the present invention. The phrase "in one embodiment" appearing in different places in this specification does not necessarily refer to the same embodiment, nor is it a single or selective embodiment that is mutually exclusive with other embodiments.

[0034] Example

[0035] Reference Figure 1This is one embodiment of the present invention, which provides a method for intelligent monitoring of geological disasters based on satellite remote sensing, which can accurately identify, assess risks, and provide early warnings for areas with potential geological disaster hazards.

[0036] Specifically, the corresponding land cover information is obtained through image registration and orthorectification.

[0037] By using land cover information, we can identify topographic boundaries and geological structural features;

[0038] Geological attribute information on lithology and stratigraphic occurrence is extracted by analyzing topographic boundaries and geological structural features.

[0039] Identify potential hazards in the surrounding area using terrain feature extraction methods;

[0040] By combining image fusion enhancement with permanent scatterer interferometry, stable scatterers are extracted and a deformation monitoring network is constructed to calculate the phase difference value in high-risk areas.

[0041] By unwinding the phase difference and combining it with orbital parameters and atmospheric phase, the three-dimensional deformation field information of the surface in the monitoring area is obtained by inversion, and the cumulative deformation is calculated by time series analysis.

[0042] A geological hazard risk assessment model was established through spatiotemporal clustering analysis.

[0043] By establishing a risk assessment model, the geological disaster risks in the monitoring area are quantitatively assessed, risk levels are classified, and early warning information is generated.

[0044] Furthermore, image registration and orthorectification processing includes acquiring multi-time-segment multispectral image data and synthetic aperture radar imagery of the area surrounding the substation using high-resolution multispectral satellites;

[0045] Multispectral image data includes visible light, near-infrared, and short-wave infrared bands;

[0046] Synthetic aperture radar imagery includes radar echo signals from the area surrounding the substation, as well as corresponding radar image data.

[0047] Specifically, high-resolution multispectral satellite sensors are used to image the area surrounding the substation over multiple time periods. Multispectral imagery data encompassing visible, near-infrared, and short-wave infrared bands is acquired quarterly. Simultaneously, a synthetic aperture radar system is used to collect radar echo signals from the area, generating corresponding radar images. For the acquired multispectral and radar images, the SIFT (Scale Invariant Eigentransform) algorithm is used to extract image feature points. This includes: performing multi-scale analysis on the input image to create a series of spatial images at different scales; comparing each pixel between adjacent scales to find locations with local maxima in scale space as keypoints; fitting candidate keypoints to remove edge-response and low-contrast keypoints; calculating the gradient histogram of the neighborhood around each keypoint and assigning one or more principal directions to each keypoint; calculating the gradient direction histogram within the neighborhood of each keypoint and quantizing it into a fixed-length vector to form a descriptor; and then using R... The ANSAC (Random Sample Consensus) algorithm removes mismatched points, establishes a geometric transformation model between images, achieves accurate registration of multi-source heterogeneous images, and eliminates geometric deviations between images. This includes: randomly selecting a minimum number of matching point pairs (at least 3 pairs are needed for affine transformations) to estimate transformation parameters; using the estimated transformation parameters to check whether the remaining matching point pairs conform to the transformation, with matching point pairs being called inliers; iterative optimization repeating the above steps multiple times, selecting different matching point pairs each time, to find the transformation parameter set with the most inliers; and finally, model construction using the most inlier set to re-estimate the transformation parameters, resulting in the final geometric transformation model. Based on SRTM (shuttle radar topography mission) digital elevation model data and RPCs (rational polynomial coefficients) models, orthorectification is performed on the registered multispectral and radar images to eliminate geometric distortions caused by terrain undulations and imaging perspective, resulting in orthorectified images with a unified geographic coordinate system. Orthorectification includes: loading SRTM / DEM data to ensure its coverage matches the multispectral and radar images to be rectified; extracting RPC coefficients from image metadata, which describe the sensor's imaging geometry; calculating the terrain height corresponding to each pixel using DEM data—this step is necessary because terrain height affects image geometry; inversely projecting each pixel into the ground coordinate system using the RPC model, taking into account the terrain effect as well as the sensor's attitude and position; and resampling the inverse projection results to ensure they fall on a standard grid to generate the orthorectified image. Resampling methods include nearest neighbor, bilinear interpolation, and cubic convolution interpolation. The final output orthophoto is topographically and visually corrected, eliminating geometric distortions caused by topographic undulations and imaging perspective. All pixels are located on the same horizontal plane and have a unified geographic coordinate system.Atmospheric correction and radiometric calibration were performed on the orthorectified multispectral images to eliminate atmospheric effects and sensor errors, improving the spectral information quality of the images. Principal component analysis was then used to reduce the dimensionality of the multi-band images and extract key feature information. Land cover classification features were constructed using multi-temporal orthorectified image data. A support vector machine classification algorithm based on radial basis function kernels was employed, with grid search determining the optimal penalty and kernel parameters. The reduced image features were then used for training and classification to identify the land cover types in the area surrounding the substation, generating a high-precision land cover classification map. High-resolution multispectral satellite sensors acquired images within a 10-kilometer radius of the substation, with a spatial resolution of 0.5 meters, including blue, green, red, and near-infrared bands, acquired quarterly. Simultaneously, an X-band synthetic aperture radar system operating at 9.6 GHz was used to acquire radar images. The SIFT algorithm was applied to the acquired multi-source images, with a feature point threshold of 1000. After feature point extraction, the RANSAC algorithm was used to remove mismatched points, with an inlier ratio threshold of 0.8. Using the SRTM digital elevation model with a vertical accuracy of 16 meters and a horizontal resolution of 30 meters, orthorectification was performed using the RPC model to obtain orthorectified imagery in the WGS84 geographic coordinate system. The FLAASH atmospheric correction model was applied to the orthorectified multispectral imagery, selecting a mid-latitude summer atmospheric model with a visibility setting of 40 kilometers for atmospheric correction and radiometric calibration. Principal component analysis was used to reduce the dimensionality of the corrected multiband imagery, retaining principal components with a cumulative contribution rate of 95%. Land cover classification features including the Normalized Difference Vegetation Index (NDVI), Normalized Difference Water Index (NDWI), and Normalized Difference Building Index (NDBI) were constructed. Support vector machines with radial basis function kernels were used for classification. The optimal penalty parameter C was determined to be 100, and the kernel parameter γ to be 0.01, determined through five-fold cross-validation and grid search. Finally, a land cover classification map with a spatial resolution of 0.5 meters was generated, including four land cover categories: water bodies, vegetation, buildings, and bare land.

[0048] SIFT is a feature extraction method used to detect and describe key points in images at different scales and orientations. The main parameters of SIFT include scale threshold, edge threshold, and key point threshold, which are used to determine the stability and uniqueness of key points.

[0049] RANSAC is an iterative method used to estimate the parameters of a mathematical model from a set of data containing outliers. The main parameters of RANSAC include the number of iterations, the interior point threshold, and the model fit threshold. These parameters determine the robustness and accuracy of the algorithm.

[0050] The SRTM digital elevation model provides high-resolution terrain data for the entire globe. SRTM data is typically provided in the form of a digital elevation model (DEM), with parameters including resolution and accuracy.

[0051] RPCs are a set of polynomial coefficients used for geometric correction of satellite images. RPCs parameters include coefficients of linear, quadratic, and cubic terms used to convert image pixel coordinates into ground coordinates.

[0052] The NDVI formula is as follows:

[0053]

[0054] Wherein, NIR is the reflectance in the near-infrared band, and R is the reflectance in the red band;

[0055] The NDWI formula is as follows:

[0056]

[0057] Where G is the reflectivity of the green light band;

[0058] The NDBI formula is as follows:

[0059]

[0060] Wherein, MIR is the reflectance of the mid-wave infrared band, and SWIR1 is the reflectance of the short-wave infrared band.

[0061] The FLAASH atmospheric correction model can remove atmospheric effects from remote sensing imagery to obtain more accurate surface reflectance. Based on a physical atmospheric radiative transfer model, the FLAASH model considers the effects of gases, aerosols, and clouds in the atmosphere. Its formula is as follows:

[0062] I(λ)=I0(λ)·τ(λ)+(1-τ(λ))·L(λ)

[0063] Where I(λ) is the observed radiance (including atmospheric path radiation), I0(λ) is the radiance at the top of the atmosphere (i.e., surface reflected radiation), τ(λ) is atmospheric transparency, representing the attenuation of radiation by the atmosphere, and L(λ) is atmospheric path radiation.

[0064] Furthermore, the land cover information includes extracting land texture information through texture feature extraction methods, and identifying terrain boundaries and geological structure features by combining edge detection and segmentation methods;

[0065] Texture feature extraction methods include calculating texture features using a gray-level co-occurrence matrix to obtain surface texture information, which is a surface texture feature map;

[0066] The edge detection and segmentation method involves detecting the enhanced texture feature map to obtain terrain boundary segments; based on the terrain boundary segments, the watershed algorithm is applied to segment the image of the edge detection results, a minimum region area threshold is set, and regions smaller than the threshold are merged to obtain the geological structure feature regions.

[0067] Furthermore, the watershed algorithm includes gradient calculation, minimum value marking, flood filling, watershed meeting, post-processing, minimum region area thresholding, and morphological operations.

[0068] Furthermore, the geological attribute information of lithology and stratigraphic attitude includes extracting shape index and directional features from the geological structural feature regions, calculating the circularity and main direction angle of each segmented region, combining the orientation information of the topographic boundary line segments, identifying lithology categories and stratigraphic dip angles, determining stratigraphic attitude, and obtaining geological attribute information.

[0069] Specifically, based on the acquired land cover information, texture features are calculated using a gray-level co-occurrence matrix. By setting the sliding window size to 9×9 pixels, statistics such as the second moment of angles, contrast, entropy, and correlation in four directions (0°, 45°, 90°, and 135°) are calculated to obtain a land texture feature map. Adaptive histogram equalization is then applied to the texture feature map to improve image contrast and enhance texture details. For the enhanced texture feature map, the Canny edge detection operator is used for edge detection. The Canny edge detection operator processes the texture feature map, including Gaussian filtering for noise reduction, calculating gradient strength and direction, non-maximum suppression, and dual-threshold detection. The OTSU algorithm is used to adaptively calculate high and low thresholds to detect terrain boundary segments. The OTSU algorithm adaptively determines the high and low thresholds required by the Canny edge detection operator by maximizing the inter-class variance between the foreground and background. A distance transform is used to calculate a marker map, which indicates seed points or initial segmentation regions in the image. This map measures the distance from each pixel to the nearest edge point. Then, a watershed algorithm is applied to segment the image based on the edge detection results. The image is treated as terrain, where pixel intensity represents altitude. The segmentation goal is to find "valleys" and "peaks," i.e., region boundaries. A minimum region area threshold of 100 square pixels is set to filter out small, unimportant regions, merging those smaller than the threshold to obtain geological structural feature regions. Shape indices and directional features are extracted from the segmented geological structural feature regions. The circularity and principal orientation angle of each segmented region are calculated. Combined with the orientation information of the terrain boundary lines, a support vector machine classifier is used to identify lithology and stratigraphic dip angles, determining the stratigraphic attitude and obtaining geological attribute information. A Support Vector Machine (SVM) classifier is trained using known geological sample data. A Radial Basis Function (RBF) kernel is employed, and optimal parameters are determined through grid search. This process includes ensuring the geological sample data has been preprocessed (e.g., cleaning, standardization / normalization) to suit its training purpose. The data should be divided into features (X) and labels (y), where features are the training dataset and labels are the corresponding categories or target variables. An SVM classifier is selected, specifying the use of the RBF kernel. The RBF kernel is well-suited for handling non-linearly separable datasets because it maps the original feature space to a higher-dimensional space, where a suitable hyperplane can be more easily found to partition the data. Two key hyperparameters of the SVM require tuning: the regularization parameter C and the width parameter gamma of the RBF kernel. The choice of these parameters affects the model's generalization ability. To find the optimal parameter combination, a grid search is used to traverse all possible combinations, and cross-validation is employed to evaluate the effectiveness of each combination. Based on geological principles and known regional geological background information, the classification results are post-processed to correct abnormal classifications that do not conform to geological patterns, thereby improving the accuracy of geological attribute information.Texture features of land cover information are calculated using the gray-level co-occurrence matrix (GLCM). A 9×9 pixel sliding window is used to calculate the second moment, contrast, entropy, and correlation of the angles at 0°, 45°, 90°, and 135°, resulting in a 16-dimensional texture feature vector. Adaptive histogram equalization is applied to the texture feature map, with an image patch size of 8×8 pixels and a contrast limit threshold of 0.02 to enhance texture details. Edge detection is performed using the Canny edge detection operator, and the OTSU algorithm is used to adaptively calculate high and low thresholds. Typically, the high threshold is approximately 120 and the low threshold is approximately 60. A marker map is calculated using distance transform, with the Euclidean distance from the pixel to the boundary used as the transform value. Then, the watershed algorithm is applied for image segmentation, with a minimum region area threshold of 100 square pixels. Shape index and directional features are extracted from the segmented regions, and roundness (4π×area / perimeter²) and principal direction angle (ellipse fitting based on the region's second moment) are calculated. Support vector machine (SVM) classifiers were used to identify lithology and stratigraphic attitude. The training samples included 5 common lithologies and 3 stratigraphic attitudes. A radial basis function (RBF) kernel was used, and a penalty factor C of 10 and a kernel parameter γ of 0.01 were determined through grid search. Based on regional geological background information, such as the known stratigraphic strike of 45°E, the classification results were post-processed to correct outlier classifications that deviated from the main strike by more than 30°, improving the accuracy of geological attribute information to over 90%.

[0070] The main parameters of the Canny edge detection operator include the standard deviation of the Gaussian filter, the gradient threshold, and the suppression threshold.

[0071] The OTSU algorithm is an automatic image thresholding method, and its formula is as follows:

[0072]

[0073] The main parameter of the OTSU algorithm is the threshold T, which divides the image histogram into foreground and background, minimizing the intra-class variance between them.

[0074] Furthermore, terrain feature extraction methods include extracting terrain slope and aspect information from geological attribute information;

[0075] The surrounding potential hazards include high-risk areas prone to geological disasters such as landslides, collapses, and ground subsidence, identified through topographic slope and aspect information, forming a distribution map of potential hazards around the substation.

[0076] Specifically, based on the elevation data in the geological attribute information, the terrain slope is calculated using the finite difference method. The slope value is obtained by the ratio of the elevation difference between adjacent grid cells to the horizontal distance. The slope is divided into five levels: 0-5°, 5-15°, 15-25°, 25-35°, and >35°, for comprehensive assessment of potential risk areas. The terrain aspect is calculated using the elevation data. The aspect angle is determined by the elevation gradient of the grid cells in the east-west and north-south directions, obtaining aspect values ​​from 0 to 360 degrees. These values ​​are then divided into eight 45° sector areas, corresponding to the eight aspect categories: north, northeast, east, southeast, south, southwest, west, and northwest. Combining slope and aspect information with geological attributes, a random forest algorithm is used to establish a geological hazard risk assessment model. Data cleaning and feature engineering are performed on historical hazard data obtained from geological departments to train model parameters and calculate the probability of landslides, collapses, and ground subsidence. This includes: obtaining location information of historical landslides, collapses, and ground subsidence events from geological departments; acquiring DEM data of the study area and calculating slope and aspect information from it; collecting geological attribute data such as soil type, rock strata properties, and groundwater level; checking for missing values ​​in the dataset and deciding how to handle them (deletion, imputation, etc.); identifying and handling outliers to avoid their negative impact on model training and ensuring the consistency and reliability of all data sources. Reliability; Extract slope and aspect information from DEM data, as these two features are crucial for geological hazard risk assessment; Extract useful features from geological attribute data, such as soil permeability and rock compressive strength; Determine which features are most important for predicting geological hazards through correlation analysis or other feature selection methods; Divide the dataset into training and test sets, and select the random forest algorithm as the geological hazard risk assessment model; Adjust the hyperparameters of the random forest, such as the number of trees and maximum depth, using methods such as cross-validation, and train the random forest model using the training set data; Evaluate the model's performance on the test set, such as accuracy, recall, and F1 score; For each test sample, the model will output the probability of landslides, collapses, and ground subsidence. Based on the calculation results of the risk assessment model, a risk threshold is set, and the 90th percentile is selected as the high-risk threshold to identify high-risk areas. A kernel density estimation method is used to perform spatial cluster analysis on high-risk points. Using GIS software or other visualization tools, a heat map is drawn based on the density estimation results to show the distribution of potential hazards around the substation. This includes: using a selected bandwidth to estimate the density of high-risk points and obtaining the density value at each point; using the density values ​​obtained from KDE, high-density areas can be identified, and these areas are potential hazard clusters. The density threshold can be used to determine which areas belong to the same cluster, thereby identifying spatial hazard clusters.Combining the substation's location information and the distribution of surrounding infrastructure, a secondary assessment of high-risk areas is conducted, considering the impact of human activities on the geological environment, such as excavation and filling projects. Weighted overlay analysis is used to adjust the risk assessment results, resulting in a final map of potential hazard points around the substation. Based on a 10-meter resolution digital elevation model, the finite difference method is used to calculate the terrain slope. A 3×3 pixel window is used to calculate the maximum slope value of the center pixel, classifying the slope into five levels: 0-5°, 5-15°, 15-25°, 25-35°, and >35°, assigning risk values ​​from 1 to 5. The same elevation data is used to calculate the terrain aspect, and the Zevenbergen & Thorne algorithm is used to calculate the aspect angle, obtaining aspect values ​​from 0 to 360 degrees. These are then divided into eight 45° sector regions, corresponding to the eight aspect categories: North, Northeast, East, Southeast, South, Southwest, West, and Northwest, each assigned a coding value from 1 to 8. By combining slope and aspect information with geological attributes such as lithology and dip angle, a random forest model was constructed using 500 decision trees, with a maximum tree depth of 10 and a minimum leaf node sample size of 5. Data from 100 historical disaster sites over the past 20 years, provided by the geological department, was used to divide the data into an 80% training set and a 20% test set. After cross-validation, the probabilities of landslides, collapses, and ground subsidence were obtained. A 90th percentile was set as the high-risk threshold. Spatial clustering analysis was performed on high-risk points using kernel density estimation (with a bandwidth of 1000 meters), generating a 100×100 meter grid heatmap of hazard point distribution. Incorporating engineering activity information within the 500-meter buffer zone of the substation, such as excavation depth and fill height, adjustment coefficients of 0.8-1.2 were set. Weighted overlay analysis was used to fine-tune the risk assessment results, ultimately forming a 1:10000 scale map of hazard point distribution around the substation.

[0077] The Zevenbergen & Thorne algorithm includes slope calculation, flow direction determination, water flow accumulation, threshold setting, river network extraction, and river network refinement.

[0078] The slope calculation formula is as follows:

[0079]

[0080] Among them, z i,j Δx is the elevation of cell (i, j), and Δy is the size of the cell.

[0081] Furthermore, the image fusion enhancement method includes: using discrete wavelet transform to fuse optical images and synthetic aperture radar images; taking the average value of low-frequency coefficients; using the absolute value maximum method to select high-frequency coefficients; extracting high-frequency and low-frequency information through multi-scale decomposition; using different weight coefficients for fusion reconstruction to obtain the fused high-resolution image; and georegistering the fused high-resolution image with the hazard point distribution map.

[0082] The permanent scatterer interferometry method includes selecting points in the high-risk area with amplitude deviation index less than a preset threshold and coherence greater than a preset threshold as candidate permanent scatterers based on the registered fused image and the distribution map of potential hazard points. The candidate permanent scatterers are then screened, and time-series interferometry is performed on the screened stable scatterers. Specifically, a deformation monitoring network is constructed using a spatial Delaunay triangulation network, and the network adjustment is performed using the least squares adjustment method to obtain the relative phase difference between adjacent monitoring points in the high-risk area.

[0083] Specifically, discrete wavelet transform is used to fuse optical and synthetic aperture radar (SAR) images. A three-level decomposition is performed using the db4 wavelet basis function. Low-frequency coefficients are averaged, while high-frequency coefficients are selected using the absolute value maximization method. High-frequency and low-frequency information is extracted through multi-scale decomposition, and fusion reconstruction is performed using different weighted coefficients to obtain the fused high-resolution image. This includes: performing three-level wavelet decomposition on the two original images (optical and SAR images); for the low-frequency coefficients of each level, the average of the corresponding low-frequency coefficients from the two images is taken; for the high-frequency coefficients of each level, the absolute value maximization method is used (selecting the coefficient with the larger absolute value as the final fusion coefficient); after fusing the low-frequency and high-frequency coefficients, the information extracted from multi-scale decomposition is used to reconstruct the fusion coefficients through inverse discrete wavelet transform (IDWT) to obtain the fused high-resolution image. The fused high-resolution image is then geo-registered with the hazard point distribution map using an affine transformation model and nearest neighbor resampling method to ensure accurate spatial correspondence between the two. Based on the registered fused image and the distribution map of potential hazard points, points within high-risk areas with an amplitude deviation index less than 0.3 and a coherence greater than 0.7 were selected as candidate permanent scatterers. The phase-to-baseline relationship of the scatterers was fitted using the least squares method, and spurious scattering points with residuals greater than a preset threshold were eliminated. Time-series interferometry was performed on the selected stable scatterers. By constructing interferogram pairs and using the minimum cost flow algorithm for phase unwrapping, the cumulative deformation and deformation rate of each monitoring point were obtained.A deformation monitoring network is constructed using a spatial Delaunay triangulation network. Network adjustment is performed using the least squares adjustment method to obtain the relative phase difference between adjacent monitoring points in high-risk areas, thereby determining the deformation trend and potential risks in the region. This process includes: selecting a certain number of monitoring points, which should be evenly distributed across the monitoring area, with increased density in high-risk areas to ensure sufficient monitoring accuracy; installing appropriate monitoring equipment, such as GNSS receivers or other deformation monitoring instruments, at the selected monitoring points; connecting these monitoring points using the Delaunay triangulation algorithm to form a triangular network without intersecting boundaries, ensuring that the maximum gap circle of each triangle does not contain any other monitoring points; after the triangulation network is constructed, the three vertices within each triangle constitute an interconnected monitoring unit, facilitating subsequent deformation analysis; and periodically or continuously collecting phase measurement data from each monitoring point. These data can then be used to... The process involves: analyzing carrier phase observations from GNSS receivers; calculating the baseline vector for each triangle side (i.e., the line connecting two monitoring points); calculating the relative phase difference for each baseline vector using the least squares method, which reflects the relative motion between the two points; performing adjustment calculations on the entire monitoring network using the least squares method to reduce measurement errors and improve the accuracy of monitoring results; analyzing the relative phase difference within each triangle based on the adjustment results to determine the deformation trend within the region; identifying potential deformation trends by observing changes in the relative phase difference between monitoring points over different time periods; assessing the potential risk level based on the deformation trend; for example, if monitoring points in a certain area show a consistent direction of movement, it may indicate a deformation risk in that area; summarizing the analysis results and compiling a detailed deformation monitoring report to provide decision support for relevant departments and to take necessary preventive or emergency response measures. Discrete wavelet transform is used to fuse 10-meter resolution optical images and 30-meter resolution synthetic aperture radar images. A three-level decomposition is performed using the db4 wavelet basis function, with low-frequency coefficients averaged and high-frequency coefficients selected using the absolute value maximization method. The fused reconstruction yields a high-resolution 10-meter resolution image. The fused imagery was georegistered with a 1:10000 scale map of potential hazard points using a 6-parameter affine transformation model and nearest-neighbor resampling, achieving a registration accuracy better than 0.5 pixels. Points with an amplitude deviation index less than 0.3 and a coherence greater than 0.7 within high-risk areas were selected as candidate permanent scatterers, with an average of 50-100 points per square kilometer. The phase-to-baseline relationship of the scatterers was fitted using the least squares method, with a residual threshold of 0.5 radians set to eliminate spurious scatterers exceeding the threshold. Time-series interferometry was performed on the selected stable scatterers, constructing 190 interferogram pairs using 20 SAR images. A minimum-cost flow algorithm was used for phase unwrapping, and the cumulative deformation and annual average deformation rate for each monitoring point were calculated.A deformation monitoring network was constructed using a spatial Delaunay triangulation network with a maximum side length of 500 meters. The network adjustment was performed using the least squares adjustment method to obtain the relative phase difference between adjacent monitoring points in high-risk areas, with an accuracy better than 0.1 radians. This was used to determine the deformation trend and potential risks in the region.

[0084] Furthermore, the phase difference unwrapping process includes using a least-squares phase expansion algorithm to unwrap the phase difference, obtaining a continuous phase field. Based on the continuous phase field, atmospheric phase correction is performed using a combination of tropospheric delay models and ionospheric delay models to obtain corrected phase information. The transformation matrix between line-of-sight deformation and three-dimensional deformation components is solved using the least-squares method, converting the corrected phase information into a three-dimensional surface deformation field. The three-dimensional surface deformation field is then decomposed into multiple scales, extracting long-term trend terms and calculating cumulative deformation.

[0085] Specifically, a least-squares phase expansion algorithm is used to unwrap the phase difference, with an upper limit of 50 iterations and a convergence threshold of 0.1 radians. Phase jumps are eliminated through iterative calculations to obtain a continuous phase field. Based on satellite orbit parameters and an atmospheric phase model, atmospheric phase correction is performed using a combination of tropospheric delay and ionospheric delay models. Baseline refinement is then performed using ground GPS control point data to improve deformation measurement accuracy. This includes: acquiring accurate satellite orbit files, such as precise ephemeris SP3 files, which contain satellite position information at different times; correcting the satellite's orbital position based on the information in the orbit files to reduce the impact of orbital errors on deformation measurements; collecting data from ground GPS control points, including receiver positions, reception times, and received satellite signals; preprocessing the collected data, such as removing gross errors and filtering; using a tropospheric delay model to estimate the impact of tropospheric delay based on meteorological parameters such as atmospheric pressure, temperature, and humidity along the satellite's line of sight; correcting carrier phase observations by subtracting the tropospheric delay value; and using an ionospheric delay model to estimate the impact of ionospheric delay based on electron content along the satellite's line of sight. Most ionospheric effects are eliminated by combining dual-frequency observations, or by using single-frequency observations combined with an external ionospheric model for correction. The baseline vector, i.e., the relative position between two monitoring stations, is calculated using the corrected carrier phase observations. The baseline vectors in the entire monitoring network are adjusted using the least squares method or other optimization algorithms to reduce measurement errors and improve accuracy. Ground GPS control points with known locations are used as constraints to further improve the accuracy of the baseline vectors. The actual coordinates of the control points are incorporated into the adjustment calculation to correct the baseline vectors, thereby obtaining more accurate monitoring point coordinates. Through these steps, the accuracy of deformation measurement can be significantly improved, specifically in the following aspects: correction of tropospheric and ionospheric delays reduces the impact of atmospheric effects on the received signal; baseline refinement and ground control point constraints improve the accuracy of the baseline vectors; and network adjustment calculations enhance the overall stability of the monitoring network, making the monitoring results more reliable. The transformation matrix between line-of-sight deformation and three-dimensional deformation components was solved using the least squares method. The corrected phase information was converted into a three-dimensional deformation field on the surface, and the east-west, north-south, and vertical displacement components of each scattering point within the monitoring area were calculated. Multi-scale decomposition was performed using db4 wavelets to extract long-term trend terms. Time series processing was then applied to the three-dimensional deformation field to calculate the cumulative deformation, with a threshold of 30 mm / year. Based on the cumulative deformation and geological background information, a fuzzy comprehensive evaluation method was used to classify potential hazard points into high, medium, and low risk levels. A distribution map of hazard points was generated, and high-risk areas were identified as potential geological hazard points.Phase difference unwrapping was performed on 300,000 scattering points within the monitoring area using a least-squares phase expansion algorithm. An upper limit of 50 iterations and a convergence threshold of 0.1 radians were set, resulting in a continuous phase field after an average of 15 iterations. A tropospheric delay model was constructed using ECMWF meteorological data, and ionospheric delay correction was performed using a global ionospheric TEC map. Baseline refinement was performed using data from 10 ground GPS control points, improving deformation measurement accuracy to the millimeter level. A 100×100 km monitoring grid was constructed, and the transformation matrix between line-of-sight deformation and three-dimensional deformation components was solved using the least-squares method. The east-west, north-south, and vertical displacement components of each scattering point were calculated, with average accuracies of 3 mm, 4 mm, and 2 mm, respectively. A 4-level db4 wavelet decomposition was performed on 20 SAR images to extract long-term trend terms and calculate cumulative deformation, setting a threshold of 30 mm / year. Based on the cumulative deformation and geological background information, the fuzzy comprehensive evaluation method is used to classify potential hazard points into levels. A membership function is set with a weight vector of [0.4, 0.3, 0.3], corresponding to the deformation, geological conditions and topographic factors, to generate a 1:10000 scale distribution map of hazard points. High-risk areas account for 5% of the total area and are identified as potential geological hazard points.

[0086] Furthermore, spatiotemporal clustering analysis includes extracting the spatial distribution characteristics and temporal evolution patterns of deformation anomaly areas, combined with historical disaster data statistics, meteorological factors, and geological attribute information.

[0087] Specifically, the DBSCAN algorithm is used to spatially cluster identified geological hazard points, with a cluster radius of 500 meters and a minimum point threshold of 10. This yields the spatial distribution characteristics of the deformation anomaly areas, determining their geometric shape and size. The process includes: preparing a dataset containing the latitude and longitude coordinates of hazard points (which can be obtained from high-risk points in previous steps) and converting these coordinates into points in two-dimensional space; using the DBSCAN algorithm from the Scikit-learn library to cluster the converted coordinates, with a cluster radius (eps) of 500 meters and a minimum point threshold (min_samples) of 10; and analyzing the cluster labels output by the DBSCAN algorithm to obtain the spatial distribution characteristics of the deformation anomaly areas, including calculating the bounding box and area information of each cluster. The STL decomposition method is then used to perform temporal evolution analysis on the deformation anomaly areas, setting a seasonal period of 12 months. By extracting trend, periodic, and random terms, the long-term trend and seasonal fluctuation patterns of the deformation anomaly areas are obtained. Based on historical disaster data, a Gaussian kernel function was used for kernel density estimation, and the Silverman criterion was adopted for bandwidth selection. The frequency distribution of disasters within the region was calculated. Combined with meteorological data such as rainfall and temperature, a disaster sensitivity index was constructed. Integrating geological attribute information such as lithology, stratigraphic dip, and fault distribution, the weights of each factor were determined using the analytic hierarchy process (AHP) to construct a geological vulnerability index. Through multiple regression analysis, the spatial distribution characteristics, temporal evolution patterns, historical disaster frequencies, meteorological sensitivity index, and geological vulnerability index of the deformation anomaly area were used as input variables to establish a geological disaster risk assessment model. The model calculated the disaster risk level within the study area, classifying it into low, medium, and high levels with thresholds of 0.3 and 0.7, respectively. This included obtaining data on deformation anomaly areas within the study area, extracting their spatial distribution characteristics and temporal evolution patterns as one input variable; obtaining historical geological disaster data within the study area and statistically analyzing the frequency of disasters as another input variable; and obtaining meteorological data within the study area and calculating the sensitivity index of meteorological factors to geological disasters as a third input variable. Geological environmental data of the study area was acquired, and the vulnerability index of the geological environment was calculated as the fourth input variable of the model. These four types of variables were then input into a pre-established multiple regression analysis model to calculate the geological hazard risk value for each location within the study area. Based on preset risk level thresholds, the calculated risk values ​​were divided into three levels: low, medium, and high, resulting in a geological hazard risk level distribution map of the study area. If the risk level of a location is high, focused monitoring and early warning will be implemented at that location; if the risk level is medium, the patrol frequency at that location will be increased; and if the risk level is low, routine monitoring will be maintained.Within a 10km × 10km study area, the DBSCAN algorithm was used to spatially cluster 1000 identified geological hazard hazard points. A cluster radius of 500 meters and a minimum point threshold of 10 were set, resulting in 15 deformation anomaly areas, with the largest area being 2.5 square kilometers and the smallest 0.3 square kilometers. Temporal evolution analysis was performed on these anomaly areas using the STL decomposition method. Using 36 months of monitoring data and a 12-month seasonal cycle, the long-term trend of an average annual deformation rate of -15 mm / year and seasonal fluctuations of ±5 mm were extracted. Using data from 200 historical hazard points over the past 50 years, a Gaussian kernel function was used for kernel density estimation, with a bandwidth calculated to be 1.2 km using the Silverman criterion, generating a 100m × 100m grid of hazard occurrence frequency distribution maps. Combined with meteorological data such as an average annual rainfall of 1200 mm and an average annual temperature of 15℃, hazard sensitivity indicators were constructed. By integrating information on five lithological types, stratigraphic dip angles ranging from 15° to 45°, and the distribution of three major fault zones within the region, the analytic hierarchy process (AHP) was used to determine the weights of 0.4, 0.3, and 0.3, respectively, to calculate a geological vulnerability index within the range of 0-1. After standardizing all the above indicators, they were input into a multiple regression model, yielding regression coefficients of 0.3, 0.25, 0.2, 0.15, and 0.1, respectively, to establish a geological hazard risk assessment model. Based on the model calculation results, the study area was divided into low, medium, and high risk zones, with thresholds of 0.3 and 0.7, respectively, with the high-risk zone accounting for 15% of the total area.

[0088] Furthermore, the risk assessment model includes a quantitative assessment of the risks of geological hazards such as landslides, collapses, and ground subsidence within the monitoring area, classifying risk levels and generating multiple types of hazard warning information, including risk type, level, scope of impact, and recommended measures.

[0089] Specifically, based on the established risk assessment model, a 30m×30m sliding window with a step size of 15m is used to analyze the topographic data within the monitoring area, calculating parameters such as slope, aspect, and topographic relief. Combined with geological attributes such as lithology and dip angle, sensitivity indices for landslides, collapses, and ground subsidence are obtained. Using multi-source remote sensing data and ground monitoring data, time series decomposition is performed using wavelet transform to extract long-term trends and periodic changes. Combined with triggering factors such as rainfall and groundwater level, hazard indices for landslides, collapses, and ground subsidence are calculated. The weights of the sensitivity and hazard indices are determined using the analytic hierarchy process (AHP), a judgment matrix is ​​constructed, eigenvectors are calculated, and consistency checks are performed. A comprehensive risk index for landslides, collapses, and ground subsidence is obtained through weighted superposition, and risk levels are classified according to preset thresholds. Based on the distribution and importance of facilities surrounding the substation, risk levels are locally adjusted, and differentiated early warning thresholds and response strategies are formulated. Based on the comprehensive risk index and risk level, combined with the geographic information around the substation, multiple types of hazard warning information are generated, including risk type, level, impact range, and recommended measures. This information is updated daily and in real-time in case of emergencies. The risk warning information is visualized and dynamically updated through a geographic information system platform. Within a 10km×10km monitoring area, a 30m×30m sliding window with a step size of 15m is used to analyze DEM data, calculating parameters such as slope range of 0-60°, eight azimuth categories of slope aspect, and topographic relief of 0-200m. Combined with five lithological types and stratum dip angle information of 15°-45°, a fuzzy comprehensive evaluation method is used to obtain the landslide, collapse, and ground subsidence sensitivity index within the range of 0-1. Using data from 12 Sentinel-1 SAR scenes and 50 ground GPS monitoring points, a 4-level decomposition was performed using db4 wavelet analysis to extract the long-term trend of the annual average deformation rate from -20 mm / year to +5 mm / year and its seasonal variation of ±10 mm. Combined with triggering factors such as annual average rainfall of 1500 mm and groundwater level changes of ±2 m, a hazard index within the 0-1 range was calculated. A 3×3 judgment matrix was constructed using the analytic hierarchy process (AHP), with a sensitivity index weight of 0.6 and a hazard index weight of 0.4, and a consistency ratio (CR) of 0.05. Weighted summation yielded a comprehensive risk index within the 0-1 range, and risk levels were classified into low, medium, and high levels based on thresholds of 0.3 and 0.7. The risk level within the 500 m buffer zone of the substation was increased by one level, and the warning threshold was lowered by 20%. The system ultimately generates early warning information that includes risk type, level, scope of impact, and recommended measures. It visualizes the risk warning map at a scale of 1:10,000 through a WebGIS platform, updates it every 24 hours, and triggers a real-time update mechanism when the daily rainfall exceeds 50 mm.

[0090] In summary, this invention, by combining high-resolution multispectral satellite imagery and SAR imagery, more accurately extracts land cover changes and geological features, effectively identifies potential geological hazard areas, and also helps in risk assessment and early warning of geological hazards, improving the monitoring capabilities for disasters such as landslides, collapses, and ground subsidence, thereby protecting personal safety, property, and power supply, and promoting social and economic stability and development.

[0091] It is important to note that the constructions and arrangements of this application shown in several different exemplary embodiments are merely illustrative. Although only a few embodiments are described in detail in this disclosure, those who consult this disclosure will readily understand that many modifications are possible (such as variations in installation arrangement, use of materials, color, orientation, etc.) without substantially departing from the novel teachings and advantages of the subject matter described in this application. For example, an element shown as integrally formed may be composed of multiple parts or elements, the position of the element may be inverted or otherwise changed, and the nature or number or position of the discrete elements may be altered or changed. Therefore, all such modifications are intended to be included within the scope of the invention. The order or sequence of any process or method steps may be changed or rearranged according to alternative embodiments. In the claims, any "support plus function" clause is intended to cover the structure performing the function described herein, and not only structurally equivalent but also equivalent in structure. Other substitutions, modifications, alterations, and omissions may be made in the design, operation, and arrangement of the exemplary embodiments without departing from the scope of the invention. Therefore, the invention is not limited to the particular embodiments but extends to a variety of modifications that still fall within the scope of the appended claims.

[0092] Furthermore, in order to provide a concise description of exemplary embodiments, not all features of actual embodiments may be omitted.

[0093] It should be understood that numerous specific implementation decisions can be made during the development of any practical implementation, such as in any engineering or design project. Such development efforts may be complex and time-consuming, but for those of ordinary skill in the art who benefit from this disclosure, the development effort will be a routine task in design, manufacturing, and production without requiring extensive experimentation.

[0094] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A method for intelligent monitoring of geological hazards based on satellite remote sensing, characterized in that: include, The corresponding land cover information is obtained through image registration and orthorectification. By using land cover information, we can identify topographic boundaries and geological structural features; Geological attribute information on lithology and stratigraphic occurrence is extracted by analyzing topographic boundaries and geological structural features. Identify potential hazards in the surrounding area using terrain feature extraction methods; By combining image fusion enhancement with permanent scatterer interferometry, stable scatterers are extracted and a deformation monitoring network is constructed to calculate the phase difference value in high-risk areas. By unwinding the phase difference and combining it with orbital parameters and atmospheric phase, the three-dimensional deformation field information of the surface in the monitoring area is obtained by inversion, and the cumulative deformation is calculated by time series analysis. A geological hazard risk assessment model was established through spatiotemporal clustering analysis. By establishing a risk assessment model, the geological disaster risk in the monitoring area is quantitatively assessed, risk levels are classified, and early warning information is generated. The image fusion enhancement method includes: using discrete wavelet transform to fuse optical images and synthetic aperture radar images; taking the average value of low-frequency coefficients; using the absolute value maximum method to select high-frequency coefficients; extracting high-frequency and low-frequency information through multi-scale decomposition; using different weight coefficients for fusion reconstruction to obtain a fused high-resolution image; and georegistering the fused high-resolution image with the hazard point distribution map. The permanent scatterer interferometry method includes selecting points in high-risk areas with amplitude deviation index less than a preset threshold and coherence greater than a preset threshold as candidate permanent scatterers based on the registered fused image and the distribution map of potential hazard points, screening the candidate permanent scatterers, and performing time-series interferometry processing on the screened stable scatterers. The phase difference unwrapping process includes: unwrapping the phase difference using a least squares phase expansion algorithm to obtain a continuous phase field; performing atmospheric phase correction based on the continuous phase field using a combination of tropospheric delay model and ionospheric delay model to obtain corrected phase information; solving the transformation matrix between line-of-sight deformation and three-dimensional deformation components using the least squares method to convert the corrected phase information into a three-dimensional deformation field of the Earth's surface; performing multi-scale decomposition of the three-dimensional deformation field of the Earth's surface to extract long-term trend terms and calculate cumulative deformation. The spatiotemporal clustering analysis includes extracting the spatial distribution characteristics and temporal evolution patterns of deformation anomaly areas, combined with historical disaster data statistics, meteorological factors, and geological attribute information; The risk assessment model includes a quantitative assessment of the risks of landslides, collapses, and ground subsidence geological hazards within the monitoring area, classifying risk levels and generating multiple types of hazard warning information, including risk type, level, scope of impact, and recommended measures.

2. The method for intelligent monitoring of geological hazards based on satellite remote sensing as described in claim 1, characterized in that: The image registration and orthorectification processing includes acquiring multi-time period multispectral image data and synthetic aperture radar images of the area surrounding the substation using high-resolution multispectral satellites. The multispectral image data includes visible light, near-infrared, and short-wave infrared bands; The synthetic aperture radar imagery includes radar echo signals from the area surrounding the substation, as well as corresponding radar image data.

3. The method for intelligent monitoring of geological hazards based on satellite remote sensing as described in claim 2, characterized in that: The land cover information includes extracting land texture information through texture feature extraction methods, and identifying terrain boundaries and geological structure features by combining edge detection and segmentation methods; The texture feature extraction method includes calculating texture features using a gray-level co-occurrence matrix to obtain surface texture information, wherein the surface texture information is a surface texture feature map. The edge detection and segmentation method includes detecting the enhanced texture feature map to obtain terrain boundary line segments; based on the terrain boundary line segments, applying the watershed algorithm to segment the image of the edge detection results, setting a minimum region area threshold, merging regions smaller than the threshold, and obtaining geological structure feature regions.

4. The method for intelligent monitoring of geological hazards based on satellite remote sensing as described in claim 3, characterized in that: The watershed algorithm includes gradient calculation, minimum value marking, flood filling, watershed meeting, post-processing, minimum region area thresholding, and morphological operations.

5. The method for intelligent monitoring of geological hazards based on satellite remote sensing as described in claim 4, characterized in that: The geological attribute information of lithology and stratigraphic attitude includes extracting shape index and directional features from geological structural feature areas, calculating the roundness and main direction angle of each segmented area, combining the orientation information of topographic boundary line segments, identifying lithology type and stratigraphic dip angle, determining stratigraphic attitude, and obtaining geological attribute information.

6. The method for intelligent monitoring of geological hazards based on satellite remote sensing as described in claim 5, characterized in that: The terrain feature extraction method includes extracting terrain slope and aspect information from geological attribute information; The surrounding potential hazards include high-risk areas prone to geological disasters such as landslides, collapses, and ground subsidence, identified by the terrain slope and slope aspect information, forming a distribution map of potential hazards around the substation.

Images