Method for monitoring high-position loose body activity based on multi-source SAR data

Abstract

The present application relates to synthetic aperture radar technical field, disclose a kind of based on multi-source SAR data monitoring high loose body activity method, to solve the problem of poor accuracy of existing high loose body activity identification mode, scheme mainly includes: constructing spatial analysis unit;Obtain multi-source SAR data and adaptive preprocessing;Maximum backscattering feature extraction and activity quantification;High loose body activity identification and multi-scale feature analysis;Activity high loose body multidimensional risk comprehensive evaluation.The present application is combined by multi-orbit cooperation, fine feature extraction and multidimensional model evaluation, which provides an effective solution to the problem of early identification of high loose body in high mountain and valley area, risk assessment inaccuracy, significantly improves the accuracy and timeliness of monitoring and early warning.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of synthetic aperture radar, and particularly relates to a method for monitoring activity of high-position loose body based on multi-source SAR data. BACKGROUND

[0002] High-position loose body is a kind of special geological disaster body distributed in high-altitude and high-steep slope areas, such as moraine, collapse accumulation and landslide accumulation. They are usually in a critical stable state, and are prone to instability under strong earthquakes, extreme rainfall or freeze-thaw action, forming chain disasters such as mudslides and landslides, with the characteristics of large disaster scale, long movement distance, strong suddenness and huge destructive power, which poses a serious threat to downstream major projects, infrastructure and people's life and property safety.

[0003] Early identification and risk assessment of the activity of high-position loose body is a key prerequisite for realizing advanced warning and active prevention and control of geological disasters. However, the existing technical means have significant limitations in dealing with such disasters:

[0004] Firstly, the traditional geological investigation means are difficult to achieve detailed investigation in a large range and high frequency due to the extremely harsh natural environment and accessibility, not only with high risk and huge cost, but also with strong subjectivity, which makes it difficult to find micro-deformation and early hidden dangers below the ground surface.

[0005] Secondly, the monitoring method based on optical remote sensing is easily affected by frequent clouds, rain, snow and night conditions in high mountain areas, with serious data missing and discontinuous observation period, which cannot realize effective observation in all-weather and all-season, and easily misses the key surface change information.

[0006] Moreover, although the Synthetic Aperture Radar (SAR) technology has the ability of all-weather millimeter-level deformation monitoring, the existing SAR analysis method has exposed fundamental defects in practical application, and the most prominent problem is the geometric distortion effect of single-track SAR data under complex terrain conditions. The radar side-looking imaging characteristics cause systematic shadow and overlap phenomena in the loose body distribution area of a specific slope direction, forming a large area of monitoring blind area. This geometric limitation makes many potential dangerous areas cannot be effectively covered, which seriously affects the integrity of the monitoring. At the same time, the existing method generally uses rough statistical indicators such as regional average to extract the backscattering characteristics. Although this processing method can suppress noise, it also flattens the spatial difference within the loose body. In fact, the instability of the loose body often begins with subtle changes in the local area, and these early signs are systematically ignored under the traditional analysis framework. More importantly, the existing method relies on simple change detection of empirical threshold, lacks a quantitative model that links SAR observations to the physical state of loose bodies, and fails to understand the physical process reflected by the change of backscattering coefficient. During the activity of the loose body, the evolution of its surface roughness, water content and internal structure will change the scattering mechanism of radar waves. This complex physical response relationship needs to be accurately described through a quantitative model. Without a quantitative model that links SAR observations to the physical state of loose bodies, it is impossible to scientifically and accurately assess the risk level of loose bodies and identify the critical state of loose body instability. SUMMARY

[0007] The present application aims to solve the problem of poor accuracy in existing high loose body activity recognition methods, and proposes a method for monitoring high loose body activity based on multi-source SAR data.

[0008] The technical solution adopted by the present application to solve the above technical problems is:

[0009] The method for monitoring high loose body activity based on multi-source SAR data comprises:

[0010] Performing hydrological analysis based on the digital elevation model of the target area to divide a plurality of watershed units, and taking each watershed unit as a basic spatial analysis framework;

[0011] Obtaining ascending SAR data and descending SAR data of the target area within a preset time window and covering the same time period, and sequentially performing thermal noise removal, radiation scaling, terrain illumination correction, speckle filtering and adaptive preprocessing of multi-temporal image registration on the ascending SAR data and descending SAR data to generate backscattering coefficient time series data;

[0012] For each pixel point in each watershed unit, in the pre-processed backscattering coefficient time series data, a weighted maximum backscattering coefficient is calculated from all valid observations in the preset time window based on the data quality factor and the incidence angle weight of the corresponding time phase; and a normalized difference index for quantifying the activity change degree of loose bodies of each pixel point is constructed based on the weighted maximum backscattering coefficients of the current monitoring period and the historical reference period.

[0013] At the pixel scale, the activity abnormal pixel points are identified from all pixel points according to the normalized difference index and based on an adaptive threshold method; at the object scale, spatially adjacent activity abnormal pixel points are aggregated into independent activity patches through a region growing algorithm, and the morphological features and average activity intensity of each activity patch are calculated; at the watershed scale, the activity intensity comprehensive index of the watershed unit is calculated according to the spatial distribution, area and intensity characteristics of all activity patches in the watershed unit.

[0014] For the watershed unit with the activity intensity comprehensive index greater than the preset comprehensive index threshold, the SAR change intensity index, the source richness index and the terrain susceptibility index are calculated respectively, the comprehensive risk index is calculated according to the SAR change intensity index, the source richness index and the terrain susceptibility index, and the risk level of high-position loose bodies is divided and warned based on the value range of the comprehensive risk index.

[0015] Further, the adaptive preprocessing process specifically includes:

[0016] Thermal noise removal: remove the system thermal noise according to the following formula:

[0017] ;

[0018] Wherein, represents the system thermal noise power, represents the Boltzmann constant, represents the system noise temperature, represents the system bandwidth, represents the noise coefficient;

[0019] Radiometric calibration: perform radiometric calibration according to the following formula:

[0020] ;

[0021] Wherein, represents the backscattering coefficient, represents the original value of the SAR data, represents the estimated noise amplitude, represents the amplitude correction factor, represents the calibration offset constant;

[0022] Terrain illumination correction: Terrain illumination correction is performed according to the following terrain illumination correction model:

[0023] ;

[0024] wherein, denotes the terrain illumination corrected backscattering coefficient, denotes the backscattering coefficient on the reference ellipsoid, denotes the reference incidence angle, denotes the local incidence angle, denotes the terrain factor, , denotes an empirical coefficient, denotes the surface roughness, denotes an exponential function;

[0025] Speckle filtering: The filter window size is dynamically adjusted according to the local coefficient of variation:

[0026] ;

[0027] wherein, denotes the filter window size, denotes the minimum size of the filter window, denotes the maximum size of the filter window, denotes the local coefficient of variation, denotes an adjustment parameter;

[0028] Multi-temporal image registration: A three-stage registration strategy based on orbital parameters, stable ground feature control point matching and mutual information optimization is adopted to perform multi-temporal image registration.

[0029] Further, the calculation formula of the weighted maximum backscattering coefficient is as follows:

[0030] ;

[0031] ;

[0032] wherein, denotes the weighted maximum backscattering coefficient of the pixel point in the time window , denotes the maximum backscattering coefficient of the pixel point in the time window , denotes an enhancement coefficient, denotes a spatial weight, denotes the neighborhood centered on the pixel point , denotes the pixel points in the neighborhood maximum backscatter coefficient within a time window, denotes a neighborhood mean, denotes a positive operator, denotes an effective number of observations within a time window, denotes a backscatter coefficient of a pixel at the th observation, denotes a data quality factor of a pixel at the th observation, denotes a local incidence angle, i.e. the angle between the radar beam and the surface normal at the pixel at the th observation, denotes an optimal incidence angle, denotes an incidence angle weight width parameter, denotes an exponential function.

[0033] Further, the constructing a normalized difference index for quantifying the degree of change in mobility of each pixel of the loose material comprises:

[0034] constructing the normalized difference index:

[0035] ;

[0036] wherein, denotes the normalized difference index of a pixel, denotes a weighted maximum backscatter coefficient of a pixel within an early time phase, denotes a weighted maximum backscatter coefficient of a pixel within a late time phase, denotes a trend enhancement factor, denotes a reference backscatter coefficient; calculating the trend of the backscatter coefficient by linear regression analysis on a plurality of historical period data: ;

[0037] wherein, denotes the annual trend slope of the backscatter coefficient of a pixel,

[0038] denotes the number of historical periods, denotes the time of the th time phase,

[0039] th observation, ​​​​​​​​​Indicates the time average. This represents the mean backscattering coefficient. Indicates the first The weighted maximum backscattering coefficient for each phase;

[0040] Applying medium-range filtering to remove isolated noise points:

[0041] ;

[0042] in, Represents the pixels after spatial filtering Normalized difference index This represents the mean value operator. Represents the pixels in the neighborhood Normalized difference index Represented by pixels A rectangular neighborhood of 3 rows by 3 columns centered on the center.

[0043] Further, based on the normalized difference index and an adaptive thresholding method, abnormal activity pixels are identified from all pixels, including:

[0044] For each pixel to be detected An adaptive thresholding method is used to determine the detection threshold:

[0045] ;

[0046] in, Represents pixels Adaptive detection threshold, Represents pixels The average of the normalized difference index of all pixels within the local neighborhood window. The significance level coefficient is represented by the coefficient. Represents pixels The standard deviation of the normalized difference index of all pixels within a local neighborhood window. This indicates the number of valid pixels within the locally moving window. Represents pixels The normalized difference index.

[0047] Calculate each pixel to be detected In the early phase and later phases The absolute change in the backscattering coefficient between Simultaneously, each pixel to be detected is calculated separately. In the early phase within the corresponding local neighborhood window and later phases coefficient of variation of backscattering coefficient and ;

[0048] When pixel The normalized difference index is greater than its corresponding adaptive detection threshold, and simultaneously satisfies as well as At that time, the pixel The pixel was identified as having abnormal activity.

[0049] Furthermore, spatially adjacent active anomalous pixels are aggregated into independent active patches using a region growing algorithm, including:

[0050] Calculate the fitness score of each pixel as a seed point for region growing, and select the initial seed point for region growing based on the fitness score. The formula for calculating the fitness score is as follows:

[0051] ;

[0052] in, Represents pixels As a fitness score for regional growth seed points Represents pixels Confidence level of pixels identified as having abnormal activity. Represents pixels Normalized difference index Represents an exponential function. Represents pixels Euclidean distance to the nearest inactive area edge, Represents the distance scale parameter. Indicates terrain gradient, Indicates the maximum terrain gradient;

[0053] For the pixels to be judged that are adjacent to the currently active patch Calculate its neighboring pixels within the current active patch. The comprehensive similarity metric value is used to determine which pixel is merged when the comprehensive similarity metric value is greater than a preset merging threshold. The formula for calculating the comprehensive similarity metric, which is then incorporated into the currently active patch, is as follows:

[0054] ;

[0055] in, Represents pixels With pixels The comprehensive similarity measure, Represents pixels Normalized difference index Represents pixels Normalized difference index Represents the spectral scale parameter, Represents pixels With pixels Spatial distance between them Represents spatial scale parameters. Indicators representing terrain continuity , Represents pixels The slope value, Represents pixels The slope value, Indicates the threshold for slope difference. These represent the corresponding weight coefficients;

[0056] Calculate the morphological compactness of the current active patch, and retain active patches whose morphological compactness is greater than a preset compactness threshold as valid active patches. The formula for calculating the morphological compactness is as follows:

[0057] ;

[0058] in, Indicates the compactness of the shape. Indicates the area of ​​the active patch. Indicates the perimeter of the active plaque. It represents pi (π).

[0059] Furthermore, the comprehensive index of activity intensity for this watershed unit is calculated using the following formula:

[0060] ;

[0061] in, This represents a comprehensive index of activity intensity for a watershed unit. Represents the first unit within the watershed. The area of ​​each active patch, Indicates the number of active patches within a watershed unit. Indicates the first The average of the normalized difference index of all pixels within an active patch. No. Topographic weighting factors for each active patch Indicators of clustering degree. Indicates the aggregation enhancement coefficient. Represents the total area of ​​a watershed unit. The topographic complexity of a watershed unit is represented by its representation. Indicates the terrain complexity suppression coefficient;

[0062] The formula for calculating the terrain weight factor is as follows:

[0063] ;

[0064] wherein, represents the elevation of the current calculation point, represents the most suitable elevation for the development of high loose bodies, represents the control coefficient of the distribution width of the elevation term, represents the slope of the current calculation point, represents the reference slope, represents the exponential adjustment coefficient of the slope term, represents the local relief of the current calculation point, represents the adjustment coefficient of the relief term, represents the exponential function.

[0065] Further, the calculation formula of the SAR variation intensity index is as follows:

[0066] ;

[0067] wherein, represents the SAR variation intensity index, represents the area of the th active patch in the basin unit, represents the number of active patches in the basin unit, represents the average value of the normalized difference index of all pixel points in the th active patch, represents the total area of the basin unit, represents the distance weight coefficient, represents the shortest distance from the th active patch to the potential motion channel, represents the characteristic distance constant, represents the spatial distribution mode factor, represents the natural constant;

[0068] The calculation formula of the material source richness index is as follows:

[0069] ;

[0070] wherein, represents the material source richness index, represents the volume of loose bodies, represents the material density correction coefficient, represents the vegetation coverage, represents the time decay coefficient, represents the time since the last activity, represents the exponential function;

[0071] The calculation formula of the terrain susceptibility index is as follows:

[0072] ;

[0073] ;

[0074] wherein, denotes the terrain susceptibility index, denotes the source area susceptibility score, denotes the motion channel susceptibility score, denotes the accumulation area susceptibility score, denotes the corresponding weight coefficient respectively, denotes the elevation of the current calculation point, denotes the lower limit elevation, denotes the most suitable elevation, denotes the growth rate control parameter, denotes the slope of the current calculation point, denotes the control coefficient of the elevation term distribution width, denotes the lithology factor.

[0075] Further, the calculation formula of the comprehensive risk index is as follows:

[0076] ;

[0077] wherein, denotes the comprehensive risk index, denotes a nonlinear mapping function for mapping the linear weighted sum to the [0, 1] interval, denotes the SAR variation intensity index, denotes the source richness index, denotes the terrain susceptibility index, denotes the corresponding weight coefficient respectively.

[0078] Further, the method further comprises:

[0079] For the high loose body with the comprehensive risk index greater than the preset risk index, the instability probability is calculated according to the following formula:

[0080] ;

[0081] wherein, denotes the instability probability of the high loose body in the evaluation period denotes the probability of occurrence of a disaster triggering event in the evaluation period denotes the probability of occurrence of a disaster triggering event in the evaluation period denotes the probability of occurrence of a disaster triggering event in the evaluation period represents the conditional probability of high loose body movement and transformation into a disaster under the condition of a triggering event, represents the historical average triggering rate of the evaluation period, represents the evaluation period, represents the comprehensive risk index corresponding to the high loose body, represents the half-saturation constant, represents the natural constant.

[0082] The beneficial effects of the present application are: the method for monitoring the activity of high loose bodies based on multi-source SAR data provided by the present application effectively solves the problem of terrain shielding existing in single orbit data by establishing an adaptive fusion mechanism for ascending and descending orbit SAR data and dynamically allocating weights based on terrain geometric relationships, achieving full-coverage, non-missing coverage of different slope high loose body distribution areas in complex high mountain terrain, and significantly improving the integrity of the monitoring range; using the maximum backscattering coefficient as the core feature and constructing a normalized difference index for quantification, the method can effectively capture the scattering characteristic changes caused by early weak activities such as crack development and local slip within the loose body, while suppressing the systematic errors caused by imaging condition differences, greatly enhancing the discovery ability and identification accuracy of early hidden dangers and weak activity signals; the constructed multi-scale analysis framework realizes progressive analysis from micro abnormal pixel identification, to meso activity patch aggregation, to macro basin comprehensive evaluation; by fusing three dimensions of current activity intensity, source potential energy and terrain susceptibility, a quantitative evaluation model of comprehensive risk index is established, realizing a fundamental change from traditional empirical and qualitative judgment to modeling and quantitative risk assessment based on physical mechanism, providing accurate scientific basis for risk classification and accurate early warning. Through the organic combination of multi-orbit cooperation, fine feature extraction and multi-dimensional model evaluation, the present application provides an effective solution to the problems of early identification of high loose bodies in high mountain and canyon areas and inaccurate risk assessment, significantly improving the accuracy and timeliness of monitoring and early warning. BRIEF DESCRIPTION OF DRAWINGS

[0083] Figure 1 The flowchart of the method for monitoring the activity of high loose bodies based on multi-source SAR data provided by the present application is shown. DETAILED DESCRIPTION

[0084] Due to the limitations and inefficiency of the current high loose body activity identification method, the optical remote sensing technology is severely restricted by weather conditions, and the single SAR analysis method has systematic defects, so the existing technology cannot systematically and all-weather obtain complete data, even if the data is obtained, the analysis method is too rough and empirical, and cannot interpret and quantify the early activity of loose bodies from the physical mechanism level, resulting in incomplete activity monitoring, inaccurate identification, and timely early warning.

[0085] Based on this, the technical scheme of the present application is proposed. In the present application, first, the ascending track and descending track SAR data are used simultaneously to make up for the geometric limitations of single track observation, such as shadow and overlap, reduce the observation blind area, increase the effective observation number, reduce the risk of data loss caused by geometric distortion, and ensure the continuity of data in time and space. Then, through steps such as thermal noise removal, radiation calibration, and terrain correction, system errors and terrain interference in SAR data are eliminated; data quality factors and incidence angle weights are introduced to calculate the weighted maximum backscattering coefficient, highlighting reliable observation values and suppressing noise. By comparing the weighted maximum backscattering coefficient of the current monitoring period with that of the historical reference period, a normalized difference index is constructed to sensitively capture changes in surface scattering characteristics, such as changes in scattering intensity caused by loose body displacement, and quantify activity. Then, through multi-scale anomaly detection and aggregation, progressive analysis from pixels to objects to basins is performed, which not only retains detailed information but also suppresses false anomalies through spatial aggregation. Finally, the integrity of geological processes is reflected in the basin unit. At the pixel scale, adaptive threshold-based anomaly pixels are identified to avoid false judgments of complex terrain caused by fixed thresholds; at the object scale, region growing algorithm is used to aggregate adjacent anomaly pixels into active patches, extract shape and intensity features, and suppress isolated noise; at the basin scale, combined with basin unit analysis, the activity intensity comprehensive index is calculated to reflect the spatial aggregation of loose body activity. Finally, the comprehensive risk index is calculated by combining SAR change intensity, source richness, and terrain susceptibility to realize risk level division and early warning and improve the reliability of risk assessment.

[0086] The technical solutions in the embodiments will be clearly and completely described below with reference to the drawings in the embodiments. Obviously, the described embodiments are only a part of the embodiments of the present application, not all embodiments.

[0087] Figure 1 A flowchart of a method for monitoring high loose body activity based on multi-source SAR data is shown. Please refer to Figure 1 , which includes the following steps:

[0088] Step 1, construct a spatial analysis unit:

[0089] Based on the digital elevation model of the target area, hydrological analysis is performed to divide multiple basin units, and each basin unit is taken as a basic spatial analysis framework.

[0090] It can be understood that the spatial distribution and evolution of high loose body are strictly controlled by topographic conditions, therefore, the present embodiment adopts a watershed unit as the basic spatial analysis framework. This selection is based on the basic understanding that the migration of loose body material follows gravity-driven hydrological processes. The watershed boundary is extracted through hydrological analysis of the digital elevation model, and the boundary is refined in combination with high-resolution remote sensing images to ensure that the analysis unit is consistent with the natural geomorphic unit. This spatial division method enables subsequent activity analysis to be carried out within the framework of geological disaster dynamics. The specific slope unit division can be achieved by means of existing flow concentration segmentation algorithms, such as the existing patent CN113850822A, which is relatively mature and will not be described in detail here.

[0091] Step 2, obtaining multi-source SAR data and adaptive preprocessing:

[0092] Obtain ascending SAR data and descending SAR data of the target area within a preset time window and covering the same time period, and sequentially perform thermal noise removal, radiation scaling, terrain illumination correction, speckle filtering, and adaptive preprocessing of multi-temporal image registration on the ascending SAR data and descending SAR data to generate backscattering coefficient time series data.

[0093] In the extremely complex high mountain terrain environment, obtaining high-quality SAR images that truly reflect the scattering characteristics of the ground surface is the premise of quantitative monitoring. In actual application, the present embodiment can guarantee data quality through the following systematic process:

[0094] (1) The present embodiment selects C-band data of Sentinel-1 satellite as the main information source, and its 5.6 cm wavelength achieves an ideal balance between penetration ability and surface sensitivity, which can obtain internal information of loose body shallow layer and maintain sufficient sensitivity to surface changes. The 6-day revisit period of the satellite system can effectively capture the progressive change process of loose body, and the 10-meter resolution ground range multi-view product significantly reduces the coherent speckle noise through multi-view processing, providing a high-quality data basis for fine analysis.

[0095] (2) The present embodiment focuses on the observation time in the winter window period, which has important physical significance. During December to February of the next year, the environmental conditions in the high mountain area present unique stability: the continuous low temperature makes the water in the loose body remain in solid state, eliminating the disturbance of liquid water content change on the dielectric constant; the stable snow cover forms a uniform scattering background; the vegetation enters the dormant period, greatly reducing the complexity of volume scattering. The combined effect of these factors makes the change of winter SAR backscattering coefficient more truly reflect the evolution of loose body internal structure rather than environmental noise interference.

[0096] (3) To overcome the limitation of single observation geometry, the embodiment synchronously utilizes ascending and descending SAR data. The ascending SAR data is imaged from the east side during the south-to-north flight, forming a favorable incidence angle of about 39 degrees to the east-facing slope; the descending SAR data is observed from the west side during the north-to-south flight, optimizing the monitoring capability to the west-facing slope. The geometric complementary relationship of the two kinds of orbit data ensures the complete coverage of various slope loose bodies.

[0097] (4) In the data preprocessing link, the embodiment constructs a systematic adaptive preprocessing procedure, aiming to extract the true ground scattering information from the original SAR data. Specifically, it includes:

[0098] ① Firstly, thermal noise removal is performed to restore the weak signal in high-altitude areas. The embodiment removes the system thermal noise according to the following formula:

[0099] ;

[0100] wherein, represents the system thermal noise power, represents the Boltzmann constant, represents the system noise temperature, represents the system bandwidth, represents the noise coefficient.

[0101] By accurately removing the thermal noise component, the signal quality of the loose body weak scattering area is significantly improved.

[0102] ② Subsequently, radiometric calibration is performed to convert the dimensionless numerical value into the backscattering coefficient with physical meaning, according to the following formula:

[0103] ;

[0104] wherein, represents the backscattering coefficient, represents the original value of the SAR data, represents the estimated noise amplitude, represents the amplitude correction factor, represents the calibration offset constant.

[0105] The backscattering coefficient is a dimensionless physical quantity, representing the ratio of the electromagnetic wave energy reflected by the ground target to the incident electromagnetic wave energy per unit area, and its value is usually expressed in decibels (dB). This coefficient is the core physical property of each pixel point in the SAR image, directly reflecting the dielectric and geometric properties of the ground target. Through radiometric calibration, SAR data of different time phases and different orbits are ensured to be comparable.

[0106] The terrain illumination correction is the most challenging step in SAR processing in mountainous areas. The radiation distortion caused by complex terrain seriously affects the authenticity of the scattering coefficient. In this embodiment, an improved terrain illumination correction model is used for terrain illumination correction:

[0107] ;

[0108] wherein, represents the backscattering coefficient after terrain illumination correction, represents the backscattering coefficient on the reference ellipsoid, represents the reference incidence angle, represents the local incidence angle, which is the angle between the radar beam and the local ground normal direction, and changes with the terrain, represents the terrain factor, which reflects the surface scattering characteristics. For rough surfaces such as high loose bodies, its value is determined by an empirical relationship: , represents an empirical coefficient, represents the surface roughness, represents an exponential function.

[0109] The above adaptive terrain illumination correction method effectively compensates for the radiation non-uniformity caused by terrain.

[0110] The speckle noise suppression adopts a terrain adaptive filtering strategy, which dynamically adjusts the filter window size according to the local coefficient of variation:

[0111] ;

[0112] wherein, represents the filter window size, represents the minimum size of the filter window, represents the maximum size of the filter window, represents the local coefficient of variation, which is the ratio of the standard deviation to the mean of the pixel values in the local window. A low local coefficient of variation indicates that the backscattering in the region is relatively uniform, and a high local coefficient of variation indicates that there are significant scattering differences in the region, represents an adjustment parameter.

[0113] A larger filter window is used in flat terrain areas to fully suppress noise, and a smaller filter window is used in broken terrain areas to preserve spatial details, achieving an optimal balance between noise suppression and detail preservation.

[0114] ⑤ Multi-temporal image registration employs a three-level registration strategy based on orbit parameters, stable ground feature control point matching, and mutual information optimization. Coarse registration based on orbit parameters provides pixel-level alignment, control point matching of stable ground features achieves sub-pixel accuracy, and mutual information optimization ultimately achieves a registration accuracy of 0.1 pixels. This step-by-step refinement method is particularly suitable for complex mountainous environments with local deformations.

[0115] Step 3: Maximum backscattering feature extraction and activity quantification:

[0116] For each pixel within each watershed unit, in the preprocessed backscattering coefficient time series data, the weighted maximum backscattering coefficient is calculated from all valid observations within the preset time window based on the data quality factor and incident angle weight of its corresponding time phase; based on the weighted maximum backscattering coefficients of the current monitoring period and the historical reference period, a normalized difference index is constructed to quantify the degree of activity change of loose bodies at each pixel.

[0117] It is understandable that early activity in high-altitude loose bodies often manifests as subtle changes in local structures, such as fracture development or local slippage, which can alter radar scattering mechanisms. While traditional regional statistical methods can suppress noise through averaging, they inevitably mask these crucial local signals.

[0118] Based on this, in this embodiment, the formula for calculating the weighted maximum backscattering coefficient is as follows:

[0119] ;

[0120] ;

[0121] in, Represents pixels In the time window The weighted maximum backscattering coefficient within the range, Represents pixels In the time window The maximum backscattering coefficient within, Indicates the enhancement coefficient. Indicates spatial weights, Represented by pixels The neighborhood centered on, Represents the pixels in the neighborhood In the time window The maximum backscattering coefficient within, Represents the neighborhood mean. This represents the positive operator. Indicates time window The number of valid observations within the period Indicates the first Pixels at the next observation The backscattering coefficient, Indicates the first Pixels at the next observation Data quality factor Represents the local incident angle, i.e., the first... radar beam and pixels during the second observation The angle between the direction of the surface normal and the surface. Indicates the optimal angle of incidence. This represents the incident angle weighted width parameter. This represents an exponential function.

[0122] Considering the potential for short-term disturbances and data quality variations within the winter observation window, this embodiment employs a weighted maximum value extraction algorithm. This algorithm introduces a data quality factor and an incident angle weighting function to ensure that the extracted maximum value reflects both the true scattering peak and the influence of observation geometry. Furthermore, recognizing that the activity of high-altitude loose bodies often exhibits spatial clustering characteristics, this embodiment introduces a spatial neighborhood enhancement technique. By analyzing the distribution of scattering values ​​above the mean within the neighborhood surrounding the target pixel, the ability to identify local anomalies is enhanced.

[0123] In this embodiment, the construction of the normalized difference index is crucial for quantifying changes in the activity of loose materials. Specifically, it includes the following steps:

[0124] Step 31: Construct the normalized difference index:

[0125] ;

[0126] in, Represents pixels Normalized difference index Represents pixels In the early phase The weighted maximum backscattering coefficient within the range, Represents pixels In the later phase The weighted maximum backscattering coefficient within the range, Indicates the trend enhancement factor. This represents the reference scattering coefficient.

[0127] In this embodiment, the normalized difference index is constructed by considering not only relative changes but also the influence of absolute changes, thereby improving the ability to distinguish between activities of different intensities.

[0128] Step 32: Calculate the trend of scattering coefficient variation by performing linear regression analysis on data from multiple historical periods.

[0129] ;

[0130] wherein, denotes the slope of the annual variation trend of the backscattering coefficient of the pixel point , denotes the number of historical periods, denotes the time of the th time phase, denotes the time average, denotes the backscattering coefficient average, denotes the weighted maximum backscattering coefficient of the th time phase.

[0131] The variation trend of the backscattering coefficient is calculated by linear regression analysis on multiple historical period data, thereby improving the reliability of the identification.

[0132] Step 33, applying median filtering to remove isolated noise points:

[0133] ;

[0134] wherein, denotes the normalized difference index of the pixel point after spatial filtering, denotes the median operator, denotes the normalized difference index of the pixel point in the neighborhood, denotes the 3-row-by-3-column rectangular neighborhood centered on the pixel point .

[0135] In actual applications, median filtering can be applied first to remove isolated noise points, and then morphological opening operation is used to remove small noise patches and retain active regions with a certain spatial scale. This processing effectively suppresses the influence of random noise while maintaining the integrity of the real active region.

[0136] Step 4, high loose body activity identification and multi-scale feature analysis:

[0137] At the pixel scale, active abnormal pixel points are identified from all pixel points according to the normalized difference index and based on an adaptive threshold method; at the object scale, spatially adjacent active abnormal pixel points are aggregated into independent active patches through a region growing algorithm, and the morphological features and average activity intensity of each active patch are calculated; at the watershed scale, the activity intensity comprehensive index of the watershed unit is calculated according to the spatial distribution, area, and intensity characteristics of all active patches in the watershed unit.

[0138] Specifically, the activity characteristics of high-altitude loose bodies exhibit significant spatial aggregation and scale dependence, necessitating the establishment of a multi-scale analysis system from micro to macro levels. This embodiment constructs a progressive analysis framework at three levels: pixel, object, and watershed, enabling accurate identification and quantitative assessment of activity. The specific steps include:

[0139] Step 41: For the pixel scale, based on the normalized difference index and an adaptive thresholding method, identify active abnormal pixels from all pixels, specifically including:

[0140] Step 411: For each pixel to be detected An adaptive thresholding method is used to determine the detection threshold:

[0141] ;

[0142] in, Represents pixels Adaptive detection threshold, Represents pixels The average of the normalized difference index of all pixels within the local neighborhood window. The significance level coefficient is represented by the coefficient. Represents pixels The standard deviation of the normalized difference index of all pixels within a local neighborhood window. This indicates the number of valid pixels within the locally moving window. Represents pixels The normalized difference index.

[0143] Step 412: Calculate each pixel to be detected. In the early phase and later phases The absolute change in the backscattering coefficient between Simultaneously, each pixel to be detected is calculated separately. In the early phase within the corresponding local neighborhood window and later phases coefficient of variation of backscattering coefficient and .

[0144] Step 413, when pixel The normalized difference index is greater than its corresponding adaptive detection threshold, and simultaneously satisfies as well as At that time, the pixel The pixel was identified as having abnormal activity.

[0145] The above scheme dynamically determines the detection threshold by analyzing the statistical distribution of the normalized difference index in the local window, however, relying only on the statistical threshold may misjudge the changes caused by non-geological factors as loose body activities. Therefore, the embodiment introduces a physical constraint condition based on the scattering mechanism. Real loose body activities must be accompanied by an increase in backscattering coefficient and spatial heterogeneity, and only changes that meet both statistical significance and physical reasonableness are identified as potential activity regions.

[0146] Step 42, for the object scale, the identified activity abnormal pixel points need to be aggregated into activity patches with geological significance, and this process is realized by an improved region growing algorithm, specifically including:

[0147] Step 421, calculate the fitness score of each pixel point as a region growing seed point, and select the initial seed point of region growing according to the fitness score, and the calculation formula of the fitness score is as follows:

[0148] ;

[0149] Wherein, represents the pixel point as a region growing seed point, represents the pixel point identified as an activity abnormal pixel point, represents the pixel point the normalized difference index, represents the exponential function, represents the pixel point the Euclidean distance to the nearest non-active region edge, represents the distance scale parameter, represents the terrain gradient, represents the maximum terrain gradient.

[0150] The above scoring system can ensure that the seed point selection is in a position with high activity intensity, far from the edge and relatively flat terrain, avoiding misjudgment of terrain mutations as activity centers.

[0151] Step 422, for the pixel point adjacent to the current activity patch, calculate the comprehensive similarity measure value of the pixel point and the adjacent pixel points in the current activity patch, when the comprehensive similarity measure value is greater than a preset merging threshold, the pixel point is merged into the current activity patch, and the calculation formula of the comprehensive similarity measure value is as follows:

[0152] ;

[0153] Wherein, pixel point pixel point a comprehensive similarity measure value of pixel points, pixel point a normalized difference index of pixel points, pixel point a normalized difference index of pixel points, a spectral scale parameter, pixel point a spatial distance between pixel points, pixel point a spatial scale parameter, a terrain continuity index, , pixel point a slope value of pixel points, pixel point a slope value of pixel points, a slope difference threshold value, respectively represent corresponding weight coefficients;

[0154] The above region generation merging criterion considers not only the spectral similarity but also the terrain continuity constraint in the pixel homogenization during the region growing process, thereby improving the merging accuracy.

[0155] Step 423, calculate the morphological compactness of the current active patch, and retain the active patch with a morphological compactness greater than a preset compactness threshold value as an effective active patch, and the calculation formula of the morphological compactness is as follows:

[0156] ;

[0157] wherein, the morphological compactness, the area of the active patch, the perimeter of the active patch, the circumference ratio.

[0158] Specifically, the formed active patch is verified for its reliability through morphological analysis. The compactness index can effectively distinguish the true active region and noise interference. The true loose body activity usually forms a relatively regular active patch, while the change caused by noise presents an elongated or scattered morphology.

[0159] Step 43, for the watershed scale, the embodiment constructs a high-position loose body activity intensity comprehensive index, which comprehensively considers the number, area, intensity and spatial distribution characteristics of the active patch, and the calculation formula is as follows:

[0160] ;

[0161] wherein, represents the activity intensity comprehensive index of the watershed unit, represents the area of the i-th activity patch in the watershed unit, represents the number of activity patches in the watershed unit, represents the average value of the normalized difference index of all pixel points in the i-th activity patch, represents the terrain weight factor of the i-th activity patch, reflecting the prior probability of loose body activity at the location of the patch, represents the inter-aggregation degree index, which takes a higher value when the activity patches present an aggregated distribution, represents the aggregation enhancement coefficient, represents the total area of the watershed unit, represents the terrain complexity of the watershed unit, represents the terrain complexity inhibition coefficient. wherein, the calculation of the terrain weight factor is based on the terrain control law of high loose body distribution, to ensure that the activity signal in the typical high loose body development area is appropriately enhanced, and the calculation formula is as follows:

[0162] ;

[0163] ;

[0164] wherein, represents the elevation of the current calculation point, represents the most suitable elevation for the development of high loose bodies, represents the control coefficient of the distribution width of the elevation term, represents the slope of the current calculation point, represents the reference slope, represents the exponential adjustment coefficient of the slope term, represents the local terrain relief of the current calculation point, represents the adjustment coefficient of the terrain relief term, represents the exponential function.

[0165] The multi-scale analysis framework in this embodiment ensures the sensitivity to local subtle changes and realizes accurate grasp of the overall activity trend from fine recognition at the pixel level to comprehensive evaluation at the watershed level, thereby providing reliable technical support for risk assessment of high loose bodies.

[0166] Step 5, multi-dimensional risk comprehensive evaluation of active high loose bodies:

[0167] ​​For the basin unit whose activity intensity comprehensive index is greater than the preset comprehensive index threshold, the SAR change intensity index, the material source richness index and the terrain susceptibility index are calculated respectively, the comprehensive risk index is calculated according to the SAR change intensity index, the material source richness index and the terrain susceptibility index, and the risk level of the high-level loose body is divided and warned based on the value range of the comprehensive risk index.

[0168] The activity basin unit whose activity intensity comprehensive index identified by the multi-scale feature analysis is greater than the preset comprehensive index threshold only reflects the current state of the high-level loose body. To realize reliable advanced identification and risk assessment, a more comprehensive comprehensive evaluation system needs to be established. The high-level loose body from the initial activity to the final instability is a complex geomorphological dynamic process involving multiple coupling actions of material conditions, terrain environment and triggering factors. A single SAR change index is difficult to fully describe this complexity. Based on this, the embodiment constructs a three-dimensional coupled risk evaluation framework to comprehensively evaluate the risk degree of the high-level loose body from the SAR change intensity, the material source richness and the terrain susceptibility.

[0169] Among them, the SAR change intensity index reflects the current activity state of the loose body, and is the core index of risk evaluation. When calculating the index, the embodiment not only considers the characteristics of the identified activity patch itself, but also pays special attention to the influence of its spatial distribution pattern on the risk, and the specific calculation formula is as follows:

[0170] ;

[0171] Among them, represents the SAR change intensity index, represents the area of the i th activity patch in the basin unit, represents the number of activity patches in the basin unit, represents the average value of the normalized difference index of all pixel points in the i th activity patch, represents the total area of the basin unit, represents the distance weight coefficient, represents the shortest distance from the i th activity patch to the potential motion channel, represents the characteristic distance constant, represents the spatial distribution pattern factor, represents the natural constant.

[0172] ​​​The material source richness index is used to evaluate the potential disaster scale, which is the material basis of the hazard assessment. The total amount of loose accumulation directly determines the possible disaster scale, while the mobility of the material affects the instability level. This embodiment combines the static material conditions with the dynamic evolution process, which more truly reflects the actual threat level of the material source. The calculation formula is as follows:

[0173] ;

[0174]

[0175] The terrain susceptibility index measures the natural conditions of disaster occurrence, which is the environmental constraint of the hazard assessment. The formation of high loose body disaster is a complete process from the material source area to the movement area and then to the accumulation area, and the terrain conditions of each stage jointly determine the occurrence probability and influence range of the disaster. This embodiment fully considers the characteristics of high mountain landform process, and the calculation formula is as follows:

[0176] ;

[0177] ;

[0178]

[0179] Based on the above three-dimensional indexes, this embodiment constructs a comprehensive hazard index of high loose body:

[0180] ; ​​​​​​​​​​​​​​​​​​​​​

[0181] wherein, denotes the comprehensive risk index, denotes a non-linear mapping function for mapping the linear weighted sum to the interval [0, 1], denotes the SAR variation intensity index, denotes the source richness index, denotes the terrain susceptibility index, denotes the corresponding weight coefficient, respectively.

[0182] In this embodiment, the natural breakpoint method combined with expert experience is used to divide the risk level, and the comprehensive risk index is divided into four levels: when , it is divided into extremely high risk; when , it is divided into high risk; when , it is divided into medium risk; when , it is divided into low risk.

[0183] In this embodiment, for the high-risk high-level loose body with a comprehensive risk index greater than the preset risk index, a failure probability estimation model is further established:

[0184]

[0185] wherein, denotes the failure probability of the high-level loose body within the evaluation period , denotes the probability of occurrence of a disaster triggering event within the evaluation period , denotes the conditional probability of the high-level loose body moving and transforming into a disaster under the triggering event occurrence condition, denotes the historical average triggering rate of the evaluation period, denotes the evaluation period, denotes the comprehensive risk index corresponding to the high-level loose body, denotes the half-saturation constant, denotes the natural constant. In summary, the method for monitoring the activity of high-level loose bodies based on multi-source SAR data provided in this embodiment has the following advantages:

[0186] (1) Breaks through the full-coverage monitoring technical bottleneck under complex terrain conditions

[0187] (2) The monitoring accuracy is high

[0188] ​​The prior art mostly uses single-track SAR data or DInSAR technology, and there are a large number of monitoring blind areas in the complex terrain of high mountain and valley area, especially the high loose body activity information in the radar shadow and overlap area is completely lost. The embodiment realizes all-around monitoring of various slope high loose bodies by terrain self-adaptive weighted fusion of ascending and descending track data, and establishes an intelligent fusion model based on observation geometry and terrain complexity.

[0189] (2) Realize fine identification from overall statistics to local anomaly

[0190] The traditional method uses overall statistical indicators such as regional average value or median value to represent the backscattering characteristics. Although this processing method can suppress noise, it also masks the local signal of early activity of high loose bodies, which leads to missing the best early warning opportunity. The embodiment innovatively proposes a maximum backscattering coefficient aggregation technology, which can sensitively capture the scattering anomalies caused by local structural changes such as cracks and cavities in loose bodies by extracting the peak response in the time window and combining spatial neighborhood enhancement. This method improves the detection sensitivity of early activity signs by an order of magnitude, and can identify potential unstable areas 2-3 months before macroscopic deformation, truly realizing the technical goal of "early identification".

[0191] (3) Establish a quantitative evaluation technology system of high loose body activity

[0192] Existing researches stay at the level of qualitative judgment or simple threshold classification, and cannot accurately evaluate the relative risk and evolution trend of different loose bodies, which seriously restricts the scientificity of disaster prevention decision. The embodiment constructs a three-dimensional quantitative evaluation model including change intensity, change frequency and change duration, and for the first time establishes a mathematical mapping relationship between the time series change of backscattering coefficient and the stability state of loose body by introducing historical trend factor and spatiotemporal consistency constraint. This breakthrough makes the risk of high loose body can be compared and sorted by a unified quantitative index, providing a scientific basis for the optimal allocation of limited disaster prevention resources, and also laying a technical foundation for establishing a regional dynamic early warning system.

[0193] The application will be further described in detail through specific examples.

[0194] I. Selection of research area and data preparation

[0195] The validation study of the method in this embodiment selects the hinterland of the Hengduan Mountains on the southeast edge of the Qinghai-Tibet Plateau as the target region. This region represents a typical environment with the most developed high-level geological disasters, with an altitude rising sharply from 3000 meters to 5500 meters, and extremely complex terrain. Influenced by the deep Indian monsoon, the annual precipitation varies between 600 and 1000 millimeters, and is highly concentrated in the summer, providing sufficient hydrodynamic conditions for the seasonal activity of loose bodies. More importantly, strong geological action since the Quaternary has formed a large amount of loose materials such as moraine and collapse deposits in this area, providing an ideal natural experiment field for verifying the method in this embodiment.

[0196] The data preparation strictly follows the technical specifications proposed in this embodiment. Based on the 30-meter resolution digital elevation model, the study area is finely divided into 912 watershed units through hydrological analysis, with an average unit area of 2.30 square kilometers. This division ensures the internal consistency of the analysis framework and the natural geomorphic process. SAR data collection focuses on the five consecutive winter window periods from 2018 to 2023, i.e. the key period from early December to the end of February of the next year. This period is chosen based on its unique environmental advantages - the continuous freezing of the ground eliminates the interference of water changes, the dormancy of vegetation reduces the complexity of volume scattering, and the stable snow cover provides a homogeneous scattering background. A total of 150 high-quality images of the Sentinel-1 satellite were obtained over the five years, and the combination of ascending and descending orbit observation modes ensured full coverage of the complex terrain. All raw data have undergone systematic preprocessing procedures, from thermal noise removal to radiation calibration, from terrain illumination correction to speckle suppression. Each link strictly controls the quality, and finally obtains time series data of backscattering coefficient with clear physical meaning and good spatial continuity.

[0197] II. Maximum backscattering coefficient extraction and time series construction

[0198] The weighted maximum value extraction strategy, the core innovation of this embodiment, has shown excellent performance in practical application. Taking the winter data processing from 2022 to 2023 as an example, although the 12 effective images obtained during this period have temporal differences, the stable and reliable maximum backscattering intensity is successfully extracted through careful design of the weighted processing. The quality factor is dynamically adjusted according to the local signal-to-noise ratio, with a variation range of 0.75 to 0.98 in the study area, and an average value of 0.89, indicating that the overall quality of the data is excellent. The incidence angle weight function is optimally designed with 35 degrees as the center, effectively compensating for the systematic bias caused by different observation geometries.

[0199] The superiority of this strategy is fully demonstrated in specific cases. The monitoring results of a large-scale collapse accumulation body in the northern study area show that the backscattering coefficient obtained by the traditional average value method is stable at around -12.5 decibels, with weak changes and lack of spatial details. In contrast, the maximum value extraction method successfully captures the local strong scattering signal of -9.41 decibels, clearly reflecting the abnormal changes in the internal structure. Through systematic processing of five winter data, the time series constructed reveals rich information about the evolution of loose bodies. The 23 active areas in the study area show a consistent trend of scattering enhancement, with an average annual enhancement rate of more than 0.5 decibels. Especially noteworthy is the collapse accumulation body numbered BD-01, whose backscattering coefficient has been continuously enhanced from -14.2 decibels in 2018 to -10.1 decibels in 2023. This significant long-term trend indicates that the internal stability is undergoing fundamental changes, providing a clear signal for early warning.

[0200] III. Intelligent fusion of ascending and descending track data and full coverage monitoring

[0201] The complex terrain conditions in the study area fully demonstrate the necessity and effectiveness of the fusion technology of ascending and descending track data. Terrain analysis reveals a key issue: the east-facing slope accounts for 41.3% of the total area, and the west-facing slope accounts for 38.7%. This nearly balanced slope distribution means that any single track observation will miss nearly half of the potential dangerous areas. The adaptive weight fusion algorithm in this embodiment achieves intelligent data integration by accurately analyzing the geometric relationship between terrain slope and radar line of sight. In the east-facing slope, the ascending track data is given a high weight of 0.87 due to its favorable observation angle; while in the west-facing slope, the weight of the descending track data is increased to 0.85; the transition area achieves smooth fusion through gradual weight transition. Among them, the monitoring case of the watershed unit numbered FU-156 vividly demonstrates the key role of the fusion technology. The unit has complex terrain, spanning the east and west two main slopes, and is distributed with multiple loose accumulation bodies of different sizes. When using only ascending track data, the west slope area is almost completely in the monitoring blind area, and only two obvious activities on the east slope can be identified. After intelligent fusion processing, the monitoring capability is qualitatively improved - not only the monitoring accuracy of the east slope is significantly improved, but more importantly, the hidden activity area of the watershed unit numbered WS-23 on the west slope is successfully discovered, with a backscattering anomaly value of -8.7 decibels. Subsequent field verification confirms that there is indeed a landslide body developing slowly at this location, fully proving the irreplaceability of multi-track fusion for achieving full coverage monitoring.

[0202] IV. Calculation of normalized difference index and activity identification

[0203] On the basis of obtaining high-quality timing data, the embodiment uses an enhanced normalized difference index to quantitatively depict the activity change of loose bodies. Compared with traditional indexes, the enhanced design not only retains the sensitivity of relative change by introducing a trend enhancement factor, but also incorporates the contribution of absolute change amount, significantly improving the ability to distinguish different intensity activities. In practical application, the trend enhancement factor is optimized to 0.4, and the reference scattering coefficient is selected as the average value of the stable region-15 decibels. These parameters have been verified through multiple tests to ensure the stability and reliability of the index calculation. The core challenge of activity recognition lies in the scientific determination of the threshold. The adaptive threshold strategy of the embodiment fully considers the control effect of the terrain, realizing differentiated and accurate recognition. In the alpine region above an altitude of 4500 meters, the background variability increases due to strong freeze-thaw action, so the local threshold is increased to 0.28 to avoid misjudgment; in the middle-altitude zone of 3500 to 4500 meters, a moderate threshold of 0.19 is used to balance sensitivity and accuracy; and in the relatively stable low-altitude accumulation area, a lower threshold of 0.12 is used to improve the detection ability of weak signals. This refined threshold setting initially identifies 8734 abnormal pixels in the study area. Through the double constraint verification based on the scattering physical mechanism, which requires an increase in backscattering enhancement and spatial heterogeneity, 5421 pixels are finally confirmed as real activity signals, effectively excluding 3313 false abnormalities caused by environmental interference. Spatial aggregation analysis organizes these active pixels into 67 active patches with clear geological significance, providing a reliable foundation for comprehensive evaluation.

[0204] V. Multi-scale feature analysis and comprehensive evaluation

[0205] After activity recognition, the embodiment realizes system evaluation from micro to macro through the constructed multi-scale analysis framework. The pixel-scale analysis takes seed point optimization as the starting point, and comprehensively considers activity intensity, spatial position and terrain characteristics to ensure that regional growth starts from the most representative position. The seed point SP-07 with the highest score exhibits typical activity center characteristics: the normalized difference index is as high as 0.43, the position 850 meters away from the edge ensures its representativeness, and the terrain gradient of only 0.15 indicates that it is in a relatively stable local platform. Starting from these preferred seed points, the regional growth algorithm aggregates discrete active pixels into active patches with reasonable morphology by balancing spectral similarity, spatial proximity and terrain continuity. Among the 67 patches formed, the compactness of 43 patches exceeds 0.3, showing a regular spatial morphology, which is consistent with the physical characteristics of real loose body activity.

[0206] The comprehensive evaluation at the basin scale fully demonstrates the systematic advantages of the present embodiment. Taking the basin unit numbered FU-045 as an example, the total area of the 5 identified active patches in this basin unit reaches 45.2 hectares. By integrating the area, intensity, terrain weight and spatial distribution pattern of each patch, the comprehensive index of activity intensity is calculated to be 0.71, clearly indicating that the basin is in a high activity state. Further three-dimensional risk assessment deepens the understanding of risk: the current SAR change intensity index is 0.68, reflecting the active current state; the source abundance index is 0.54, indicating that there is considerable potential disaster material; the terrain susceptibility index is as high as 0.72, revealing the terrain conditions conducive to disaster occurrence. Through the weight integration and nonlinear mapping processing determined by the analytic hierarchy process, the final comprehensive risk index is 0.79, clearly belonging to the extremely high risk level. The instability probability model based on the Poisson process further quantifies the risk level, and the estimated 10.3% annual instability probability provides a clear quantitative basis for disaster prevention decision-making.

[0207] Six, verification results and application effects

[0208] The effectiveness of the method of the present embodiment was convincingly confirmed in the system verification in the 2023 flood season. Among the 27 high-risk and extremely high-risk loose bodies identified in the early stage, 24 showed clear signs of activity during the monitoring period - including new cracks, local collapse, debris flow initiation and other phenomena, with a verification accuracy rate of 88.9%, fully demonstrating the reliability of the method. Case No. HL-2023-007 is particularly representative. In the monitoring in December 2022, the normalized difference index of this loose body suddenly mutated to 0.45, exceeding the local threshold by 137%, causing the system to pay close attention. Subsequent monitoring showed that its activity intensity index quickly rose from 0.42 in 2021 to 0.76 in early 2023, and the comprehensive risk index reached an extremely high level of 0.82. On July 15, 2023, the loose body as expected occurred local instability, forming a debris flow of 78,000 cubic meters, which was highly consistent with the model's estimated scale of 82,000 cubic meters, again verifying the accuracy of the method.

[0209] Through the establishment of a business operation system in the study area, the present embodiment has realized the transition from experimental verification to normalized application. Compared with the traditional manual survey method, the present embodiment not only improves the identification accuracy from about 60% to nearly 90%, but also realizes an average early warning period of 8.5 months, which saves valuable decision-making and action time for geological disaster prevention and control. These practical achievements fully demonstrate the important application value and broad popularization prospects of the present embodiment in early identification and risk prevention of high-level geological disasters.

Claims

1. A method for monitoring high-level loose material activity based on multi-source SAR data, characterized in that, The method comprises: performing hydrological analysis based on a digital elevation model of the target region to divide a plurality of watershed units, and taking each watershed unit as a basic spatial analysis framework; acquiring ascending SAR data and descending SAR data of the target region within a preset time window and covering the same time period, and sequentially performing thermal noise removal, radiation scaling, terrain illumination correction, speckle filtering, and adaptive preprocessing of multi-temporal image registration on the ascending SAR data and the descending SAR data to generate backscattering coefficient time series data; for each pixel point in each watershed unit, based on a data quality factor and an incident angle weight of the corresponding phase in the preprocessed backscattering coefficient time series data, a weighted maximum backscattering coefficient is calculated from all effective observation values within the preset time window; and based on the weighted maximum backscattering coefficients of the current monitoring period and the historical reference period, a normalized difference index for quantifying the activity change degree of each pixel point of loose body is constructed; at a pixel scale, activity abnormal pixel points are identified from all pixel points according to the normalized difference index and based on an adaptive threshold method; at an object scale, spatially adjacent activity abnormal pixel points are aggregated into independent activity patches through a region growing algorithm, and the morphological features and average activity intensity of each activity patch are calculated; and at a watershed scale, an activity intensity comprehensive index of the watershed unit is calculated according to the spatial distribution, area, and intensity features of all activity patches in the watershed unit; for the watershed unit with the activity intensity comprehensive index greater than a preset comprehensive index threshold, a SAR change intensity index, a material source richness index, and a terrain susceptibility index are calculated respectively, a comprehensive risk index is calculated according to the SAR change intensity index, the material source richness index, and the terrain susceptibility index, and a high-level loose body risk level is divided and warned based on the value range of the comprehensive risk index; the construction of the normalized difference index for quantifying the activity change degree of each pixel point of loose body comprises: constructing the normalized difference index: ； wherein, represents a normalized difference index of pixel points represents a weighted maximum backscatter coefficient of pixel points in an early phase represents a weighted maximum backscatter coefficient of pixel points in a late phase represents a trend enhancement factor, represents a reference scatter coefficient;​​​ calculating the change trend of the scattering coefficient through linear regression analysis on a plurality of historical period data: ； wherein, represents the slope of the annual trend of the backscatter coefficient of the pixel point , represents the number of historical periods, represents the time of the th phase, represents the time average, represents the backscatter coefficient average, represents the weighted maximum backscatter coefficient of the th phase; applying median filtering to remove isolated noise points: ； wherein, denotes the normalized difference index of the pixel point after spatial filtering, denotes the median operator, denotes the normalized difference index of the pixel points in the neighborhood, denotes the 3x3 rectangular neighborhood centered on the pixel point .

2. The method for monitoring high seated loose material activity based on multi-source SAR data according to claim 1, characterized in that, the adaptive preprocessing process specifically comprises: thermal noise removal: removing system thermal noise according to the following formula: ； wherein, represents the system thermal noise power, represents the Boltzmann constant, represents the system noise temperature, represents the system bandwidth, represents the noise figure; radiation scaling: performing radiation scaling according to the following formula: ； wherein denotes the backscatter coefficient, denotes the original value of the SAR data, denotes the estimated noise amplitude, denotes the amplitude correction factor, denotes the scaling offset constant; terrain illumination correction: performing terrain illumination correction according to the following terrain illumination correction model: ； wherein denotes the backscatter coefficient after topographic illumination correction, denotes the backscatter coefficient on a reference ellipsoid, denotes the reference incidence angle, denotes the local incidence angle, denotes the topographic factor, , denotes an empirical coefficient, denotes the surface roughness, denotes an exponential function; speckle filtering: dynamically adjusting the filter window size according to the local coefficient of variation: ； wherein, denotes the filter window size, denotes the minimum size of the filter window, denotes the maximum size of the filter window, denotes the local coefficient of variation, denotes the adjustment parameter; multi-temporal image registration: adopting a three-level registration strategy based on orbit parameters, stable feature control point matching, and mutual information optimization to perform multi-temporal image registration.

3. The method for monitoring high seated loose material activity based on multi-source SAR data according to claim 1, characterized in that, The calculation formula of the weighted maximum backscattering coefficient is as follows: ； ； wherein denotes the pixel point the weighted maximum backscatter coefficient within the time window , denotes the pixel point the maximum backscatter coefficient within the time window , denotes the enhancement coefficient denotes the spatial weight denotes the neighborhood centered at the pixel point , denotes the maximum backscatter coefficient within the time window for the pixel points in the neighborhood, denotes the neighborhood mean denotes the take positive operator denotes the effective number of observations within the time window , denotes the backscatter coefficient of the pixel point at the th observation, denotes the data quality factor of the pixel point at the th observation, denotes the local incidence angle, i.e. the angle between the radar beam and the surface normal at the pixel point at the th observation, denotes the optimal incidence angle denotes the incidence angle weight width parameter denotes the exponential function.

4. The method for monitoring high seated loose material activity based on multi-source SAR data according to claim 1, characterized in that, identifying activity abnormal pixel points from all pixel points according to the normalized difference index and based on an adaptive threshold method comprises: For each pixel point to be detected An adaptive threshold method is used to determine the detection threshold: ； wherein, denotes the adaptive detection threshold value of the pixel point , denotes the average value of the normalized difference index of all pixel points within the local neighborhood window of the pixel point , denotes the saliency level coefficient, denotes the standard deviation of the normalized difference index of all pixel points within the local neighborhood window of the pixel point , denotes the number of valid pixels within the local moving window, denotes the normalized difference index of the pixel point . calculating the backscattering coefficient of each pixel to be detected at the early phase and the late phase the absolute variation of the backscattering coefficient between the early phase ; and calculating the backscattering coefficient of each pixel to be detected at the early phase and the late phase the coefficient of variation of the backscattering coefficient of each pixel to be detected at the early phase and the late phase When the normalized difference index of the pixel point is greater than its corresponding adaptive detection threshold, and simultaneously satisfies and , the pixel point is identified as an active abnormal pixel point.

5. The method for monitoring high seated loose material activity based on multi-source SAR data according to claim 1, characterized in that, aggregating spatially adjacent activity abnormal pixel points into independent activity patches through a region growing algorithm comprises: The fitness score of each pixel point as a region growing seed point is calculated, and an initial seed point of region growing is selected according to the fitness score, and the calculation formula of the fitness score is as follows: ； wherein, represents a pixel point a fitness score of a region growing seed point, represents a pixel point a confidence of being identified as an active abnormal pixel point, represents a pixel point a normalized difference index, represents an exponential function, represents a pixel point a Euclidean distance to the nearest non-active region edge, represents a distance scale parameter, represents a terrain gradient, represents a maximum terrain gradient; For the pixel point to be judged adjacent to the current active patch , a comprehensive similarity measure value with the adjacent pixel points in the current active patch is calculated , and when the comprehensive similarity measure value is greater than a preset merging threshold, the pixel point is merged into the current active patch, and the calculation formula of the comprehensive similarity measure value is as follows: ； wherein, denotes a pixel point and a pixel point a combined similarity measure value, denotes a normalized difference index of a pixel point a normalized difference index of a pixel point a normalized difference index of a pixel point a normalized difference index of a pixel point denotes a spectral scale parameter, denotes a spatial distance between a pixel point and a pixel point a spatial distance between a pixel point denotes a spatial scale parameter, denotes a terrain continuity index, , denotes a slope value of a pixel point a slope value of a pixel point a slope value of a pixel point a slope value of a pixel point denotes a slope difference threshold, denotes a corresponding weight coefficient, respectively; The morphological compactness of the current active patch is calculated, and the active patch with the morphological compactness greater than a preset compactness threshold is reserved as an effective active patch, and the calculation formula of the morphological compactness is as follows: ； wherein, represents the compactness of the shape, represents the area of the active plaque, represents the perimeter of the active plaque, represents the circumference.

6. The method for monitoring high seated loose material activity based on multi-source SAR data according to claim 1, characterized in that, The calculation formula of the active intensity comprehensive index of the flow region unit is as follows: ； in, This represents a comprehensive index of activity intensity for a watershed unit. Represents the first unit within the watershed. The area of ​​each active patch, Indicates the number of active patches within a watershed unit. Indicates the first The average of the normalized difference index of all pixels within an active patch. No. Topographic weighting factors for each active patch Indicators of clustering degree. Indicates the aggregation enhancement coefficient. Represents the total area of ​​a watershed unit. The topographic complexity of a watershed unit is represented by its representation. Indicates the terrain complexity suppression coefficient; The calculation formula of the terrain weight factor is as follows: ； wherein, represents the elevation of the current calculation point, represents the most suitable elevation for the development of high loose bodies, represents the control coefficient of the distribution width of the elevation term, represents the slope of the current calculation point, represents the reference slope, represents the exponential adjustment coefficient of the slope term, represents the local relief of the current calculation point, represents the adjustment coefficient of the relief term, represents the exponential function.

7. The method for monitoring high seated loose material activity based on multi-source SAR data according to claim 1, characterized in that, The calculation formula of the SAR change intensity index is as follows: ； in, Indicators representing the intensity of SAR change. Represents the first unit within the watershed. The area of ​​each active patch, Indicates the number of active patches within a watershed unit. Indicates the first The average of the normalized difference index of all pixels within an active patch. Represents the total area of ​​a watershed unit. This represents the distance weighting coefficient. Indicates the first The shortest distance from an active plaque to a potential movement channel Represents the characteristic distance constant. Indicates spatial distribution pattern factor, Represents the natural constant; The calculation formula of the material source richness index is as follows: ； wherein, denotes the material source richness indicator, denotes the volume of the loose body, denotes the material density correction factor, denotes the vegetation coverage, denotes the time decay coefficient, denotes the time since the last activity, denotes the exponential function; The calculation formula of the terrain susceptibility index is as follows: ； ； wherein, denotes a terrain susceptibility index, denotes a source area susceptibility score, denotes a transport pathway susceptibility score, denotes a accumulation area susceptibility score, denotes a corresponding weight factor, respectively, denotes an elevation of the current calculation point, denotes a lower limit elevation, denotes an optimum elevation, denotes a growth rate control parameter, denotes a slope of the current calculation point, denotes a control factor for the elevation term distribution width, denotes a lithology factor.

8. The method for monitoring high seated loose material activity based on multi-source SAR data according to claim 1, characterized in that, The calculation formula of the comprehensive risk index is as follows: ； wherein, denotes the comprehensive risk index, denotes a non-linear mapping function for mapping the linear weighted sum to the interval [0, 1], denotes the SAR change intensity indicator, denotes the source richness indicator, denotes the terrain susceptibility indicator, denote the corresponding weight coefficients, respectively.

9. The method for monitoring high seated loose material activity based on multi-source SAR data according to claim 1, characterized in that, The method further comprises: For the high loose body with the comprehensive risk index greater than a preset risk index, the instability probability is calculated according to the following formula: ； wherein, represents the probability of instability of the high-level loose material within the assessment period, represents the probability of a disaster-triggering event occurring within the assessment period, represents the probability of a disaster-triggering event occurring within the assessment period, represents the probability of a disaster-triggering event occurring within the assessment period, represents the conditional probability of the high-level loose material moving and translating into a disaster given the occurrence of a triggering event, represents the historical average triggering rate for the assessment period, represents the historical average triggering rate for the assessment period, represents the comprehensive risk index corresponding to the high-level loose material, represents the half-saturation constant, represents the natural constant.

Images