Historical activity detection method for high-order loose body based on multi-source DEM (Digital Elevation Model) data

Abstract

The invention relates to the technical field of radar remote sensing, discloses a high-order loose body historical activity detection method based on multi-source DEM data, and aims to solve the problem that the historical activity of a high-order loose body in a glacier coverage area is difficult to accurately and quantitatively evaluate in the prior art. The scheme mainly comprises the steps that drainage basin units are automatically divided based on DEM data, and glacier influence areas are screened; multi-source DEM data are fused, and system errors are eliminated through geographic positioning deviation correction, elevation distortion deviation correction and track mode deviation correction; establishing a penetration depth physical model, and correcting the radar penetration depth in the glacier area by using the C / X wave band elevation difference optimization parameter; calculating an elevation change rate by adopting robust regression; and calculating an activity comprehensive index by integrating multiple factors, and identifying a hot spot region and a time period through probability estimation and spatio-temporal clustering. According to the method, automatic and high-precision quantitative evaluation on the historical activity of the high-position loose body is realized, and the method is particularly suitable for geological disaster early warning and risk evaluation in high mountain areas.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of radar remote sensing, and particularly relates to a high loose body historical activity detection method based on multi-source DEM data. BACKGROUND

[0002] With the continuous intensification of global warming, the geological environment of high mountain glacier regions is undergoing significant changes. The rapid retreat of glaciers exposes a large amount of loose material that was originally covered by ice bodies, forming widely distributed unconsolidated or weakly consolidated debris accumulations, i.e., high loose bodies, in high-altitude areas. These loose accumulations are extremely vulnerable to triggering by external forces such as rainfall, snowmelt, or earthquakes, and can easily transform into chain disasters such as outburst floods, debris flows, or landslides, posing a serious threat to downstream settlements, infrastructure, and ecological environment. Therefore, accurately identifying and assessing the historical activity characteristics of high loose bodies is of great scientific and practical significance for understanding the development law of disasters, conducting risk assessment, and implementing early warning.

[0003] Traditional field geological survey methods face severe challenges in high mountain glacier regions. High mountain glacier regions usually have harsh environments, rugged terrain, poor traffic accessibility, and other characteristics, making manual surveys not only costly and time-consuming, but also difficult to achieve systematic monitoring on a large scale, with obvious shortcomings in terms of spatial and temporal density and continuity of data acquisition.

[0004] The advent of remote sensing technology has provided a new technical approach for geological disaster identification in high mountain areas, but existing remote sensing monitoring methods still have obvious limitations. First, optical remote sensing image interpretation is easily disturbed by frequent cloud and fog weather in high mountain areas, resulting in a very limited effective observation time window and making it difficult to obtain continuous and usable time series image data. Second, although synthetic aperture radar interferometry technology has the ability to penetrate clouds and fog, in areas with steep terrain and dramatic surface changes where high loose bodies are distributed, the radar signal is severely out of phase, making it impossible to obtain reliable ground deformation information. In addition, analysis methods based on single time phase digital elevation models can only reflect static topographic features and cannot capture the dynamic evolution process of loose bodies under the action of glacier retreat, freeze-thaw cycles, and other actions, nor can they quantitatively assess their historical activity.

[0005] In recent years, multi-temporal DEM difference technique has made important progress in the field of glacier mass balance monitoring. By comparing and analyzing surface elevation data at different times, the spatial and temporal variation of glacier thickness and volume can be quantitatively revealed. However, these studies mainly focus on monitoring the changes of the glacier body, and have not been systematically applied to the field of high loose body activity detection. More importantly, the comprehensive utilization of multi-source DEM data faces many technical problems that have not been solved: the differences in coordinate system, spatial resolution, and height datum between DEM data obtained by different sensors at different times seriously affect the accuracy of difference analysis. In particular, the signal penetration effect of radar DEM in the glacier snow area is particularly prominent. For example, the penetration depth of the widely used SRTM C-band data in dry snow can reach 4-10 meters, and has significant spatial variability. However, the existing correction methods mostly use simple empirical constants, which do not fully consider this spatial variability feature, resulting in limited correction accuracy. In addition, the existing technology lacks systematic analysis of error propagation path and strict evaluation of result uncertainty, making the reliability of loose body activity detection based on DEM difference poor.

[0006] In summary, the existing technology has the following defects in detecting the historical activity of high loose bodies in high mountain glacier areas: first, the reliability and continuity of data acquisition are insufficient, which makes it difficult to meet the needs of time series analysis; second, the processing accuracy of multi-source DEM data is limited, especially the lack of effective correction scheme for the systematic deviation of multi-source DEM and radar penetration effect in complex high mountain terrain conditions, which makes it impossible to achieve high-precision surface change detection. SUMMARY

[0007] The present application aims to solve the technical problem that the existing technology cannot achieve high-precision continuous detection of the historical activity of high loose bodies in high mountain glacier areas, and proposes a method for detecting the historical activity of high loose bodies based on multi-source DEM data.

[0008] The technical solution adopted by the present application to solve the above technical problems is: The method for detecting the historical activity of high loose bodies based on multi-source DEM data comprises: The method for detecting the historical activity of high loose bodies based on multi-source DEM data comprises: The multi-source DEM data of the research area is acquired, including C-band and X-band DEM data of SRTM DEM, ALOS AW3D30 DEM data, and ICESat-1 and ICESat-2 laser altimetry data; the multi-source DEM data is preprocessed, including data format unification, coordinate system unification and spatial resolution unification, to obtain a standardized multi-temporal DEM data set; The standardized multi-temporal DEM data is sequentially subjected to geographic positioning bias correction, height distortion bias correction and satellite orbit pattern bias correction, to obtain high-precision multi-temporal DEM data after system bias correction; Based on the propagation characteristics of electromagnetic waves in ice and snow medium, a penetration depth physical model is established; the height difference of C-band and X-band DEM data in the high-precision multi-temporal DEM data is used to calculate the penetration depth difference, and the parameter optimization of the penetration depth physical model is performed; the optimized penetration depth physical model is applied to the ice-covered area in the candidate basin unit to perform spatialization correction on the C-band DEM data, to obtain the final DEM data after penetration depth correction; The final DEM data at different time periods is subjected to difference calculation, and a formula based on robust regression is used to fit the height change trend of each basin unit, to obtain height change rate spatio-temporal distribution data reflecting the surface dynamics; Based on the height change rate spatio-temporal distribution data, the comprehensive index of historical activity of high-level loose body is calculated in the candidate basin unit, the activity probability is estimated, and the activity hot area and time period are identified, to obtain the quantitative evaluation result of the historical activity of high-level loose body.

[0009] Further, the calculation formula of the ice coverage rate is as follows: ; Wherein, represents the ice coverage rate, and represents the row number and column number of the pixel in the grid data in the basin unit, represents the ice coverage binary function of the pixel , if the pixel is located in the ice area, then , otherwise 0, represents the basin unit, represents the total number of grid units in the basin unit.

[0010] Further, the spatial resolution unification adopts a cubic convolution interpolation method, and the formula is as follows: ; Wherein, denotes the elevation value of the output pixel to be solved, and denotes the continuous coordinate value in the target coordinate system, denotes the elevation value of the known input pixel, and denotes the discrete matrix index of the input grid, denotes the weight kernel function of cubic convolution interpolation, denotes the normalized distance in the horizontal direction, denotes the normalized distance in the vertical direction, and denote the weight values of the normalized distances in the horizontal and vertical directions, respectively; the weight kernel function is defined as follows: when , ; when , ; when , ; wherein, denotes the normalized distance, i.e. the distance from the center point of the target pixel to the center point of the adjacent original pixel, denotes the interpolation coefficient.

[0011] Further, the corresponding formula of the geographic positioning bias correction is as follows: ; wherein, denotes the elevation difference of the two DEMs, denotes the terrain slope, denotes the magnitude of the horizontal offset vector, denotes the direction of the offset vector, denotes the terrain slope direction, i.e. the orientation angle of the slope surface, denotes the ratio of the vertical bias to the average slope tangent; solved by the least squares method , the formula is as follows: ; wherein, denotes the sample point index, denotes the total number of sample points participating in the calculation, denotes the weight value of the th sample point, denotes the elevation difference of the th sample point, denotes the slope of the th sample point, denotes the slope direction of the i-th sample point.

[0012] Further, the formula of the elevation distortion bias correction is as follows: ; wherein, denotes the pixel elevation value, denotes the elevation distortion bias estimation value, i.e. the systematic elevation error at the elevation , and and denote the minimum and maximum elevation of the study area, respectively, denotes the polynomial coefficient, denotes the polynomial order, denotes the polynomial order index. The polynomial order is determined using the Akaike information criterion (AIC), which is given by: ; wherein, denotes the AIC index, denotes the total number of stable area elevation difference sample points used to fit the elevation distortion bias model, is the residual variance, which denotes the error sum of squares of the model fitting.

[0013] Further, the formula of the satellite orbit model bias correction is as follows: ; wherein, denotes the orbit-related systematic bias, i.e. the elevation error estimation value at the orbit coordinate , and denotes the pixel coordinate along the orbit direction, denotes the coordinate along the cross-orbit direction, denotes the trend surface coefficient, and denote the trend surface order, and denote the polynomial order index. The orbit coordinate is obtained by a rotation transformation, which is given by: ; wherein, and denote the pixel coordinate in the original map projection coordinate system, denotes the orbit direction angle.

[0014] Further, the formula of the penetration depth physical model is as follows: ​ ; in, Indicates the radar penetration depth, i.e., at elevation. Snow density Penetration depth under certain conditions Indicates the snow surface elevation. Indicates snow density, Indicates the reference penetration depth at sea level. Indicates the radar frequency. Represents the speed of light. Indicates the loss tangent. Indicates the thickness of the snow cover. Represents the elevation gradient coefficient. Indicates reference elevation. Represents a constant term. Represents the natural exponential function; The formula for calculating the penetration depth difference is as follows: ; in, This indicates the difference in penetration depth between C-band and X-band DEM data. and These represent the elevation values ​​of the C-band and X-band DEM data, respectively. and These represent the penetration depths of C-band and X-band DEM data, respectively. The parameter optimization formula for the physical model of penetration depth is as follows: ; in, Represents the set of pixels in a stable region. This represents the set of parameters to be optimized. This represents the elevation difference between the observed C-band and X-band DEM data. Indicates based on parameter set The difference in predicted penetration depth Represents the regularization coefficient. This represents the regularization term.

[0015] Furthermore, the formula for the robust regression is as follows: ; in, and This represents the annual rate of change in elevation. Represents the number of time series observations. Indicates the first Periodic elevation observations, Indicates the first Observation period Indicates reference elevation. represents a reference time point, represents a standard deviation of height observation, represents a Huber loss function.

[0016] Further, the calculation formula of the high loose body historical activity comprehensive index is as follows: ; wherein, represents a high loose body historical activity comprehensive index, represents an evaluation factor serial number, represents a total number of evaluation factors, represents a weight of the evaluation factor, represents a factor value of the evaluation factor, represents a confidence degree of the evaluation factor; The formula of the activity probability estimation is as follows: ; wherein, represents an activity probability, represents an intercept term, represents a regression coefficient of the evaluation factor, represents a standardized factor value of the evaluation factor, represents a natural exponential function.

[0017] Further, the identification of the activity hotspot area and period adopts a spatiotemporal clustering analysis model, and the formula of the spatiotemporal clustering analysis model is as follows: ; wherein, represents a spatiotemporal clustering statistic, represents a spatial window, represents a time window, represents a number of high activity points in the window, represents an expected value based on a Poisson distribution, represents a total number of points, represents a window likelihood function value, represents a total likelihood function value.

[0018] The application has the beneficial effects that: the high loose body historical activity detection method based on multi-source DEM data provided by the application realizes high-precision, quantification and systematization detection of the historical activity of high loose bodies in the high mountain glacier area by constructing a complete multi-source DEM data processing flow. Specifically, by using and systematizing the multi-source DEM data, the observation blind area of a single data source in a harsh environment is effectively overcome, and a continuous and reliable time series elevation data set is formed; by establishing a progressive system bias correction framework and a spatialized radar penetration depth model, the comprehensive accuracy of the DEM data is improved from the meter level to the sub-meter level, thereby fundamentally ensuring the reliability of the difference result; by constructing a multi-element fusion quantitative evaluation system, the traditional qualitative discrimination is converted into an objective quantitative activity index and probability, thereby providing a scientific basis directly supporting decision-making for precise prevention and control of geological disaster risks. BRIEF DESCRIPTION OF DRAWINGS

[0019] Figure 1 A structural schematic diagram provided for an embodiment; Figure 2 A geometric relationship schematic diagram between DEM and offset DEM provided for an embodiment; Figure 3 A system bias correction effect comparison schematic diagram provided for an embodiment; Figure 4 A typical watershed unit elevation time series change curve schematic diagram provided for an embodiment; Figure 5 A high loose body activity multi-factor comprehensive evaluation radar schematic diagram provided for an embodiment; Figure 6 An ROC curve schematic diagram of an activity detection model provided for an embodiment; Figure 7 A detection result confusion matrix of an activity detection model provided for an embodiment; Figure 8 A time and space clustering analysis result schematic diagram provided for an embodiment. DETAILED DESCRIPTION

[0020] In order to overcome the problems of the existing monitoring technology, such as large environmental limitation, poor observation continuity, serious signal incoherence and the like in the high mountain glacier area, and realize long-term, continuous and reliable surface change monitoring data acquisition in a large range of high loose body distribution area, the technical scheme of the application is proposed.

[0021] In the present application, first, the watershed unit division is carried out. The movement of high-position loose bodies is controlled by the terrain, such as collapse, debris flow, and the material moves along the slope and finally converges in the channel. Therefore, the watershed unit divided based on the ridge line and the valley line is the most ideal spatial unit for analyzing its activity. The present application uses DEM to carry out hydrological analysis, automatically extracts the watershed (ridge line) and the catchment network (valley line), forms an independent slope-channel system, ensures that the subsequent analysis is carried out in a unit with clear physical meaning, and the result conforms to the law of geomorphic evolution.

[0022] Then, the multi-source DEM data processing is carried out. The DEM data obtained by different sensors and at different times have differences in coordinates, formats and resolutions, and cannot be directly compared, and must be unified to the same reference framework. The present application converts all data to the same file format and map projection, ensures that each pixel is strictly aligned in space, and adjusts all data to the same resolution through resampling to eliminate the comparison error caused by different pixel sizes.

[0023] Then, the system deviation correction is carried out. There are systematic errors between different DEM data which are not caused by real terrain changes. The present application adopts a three-step progressive correction strategy: geographic positioning deviation correction, height distortion deviation correction and orbit mode deviation. Specifically, at the topographic relief, a small horizontal position shift will appear as a huge height difference, the geographic positioning deviation correction corrects the horizontal shift amount and direction by establishing the mathematical relationship between the height difference and the slope and slope direction; in high-altitude areas, ground control points are rare, and DEM errors are often systematically related to altitude, the height distortion deviation correction fits a polynomial relationship between height and error as a function of altitude to eliminate it; the orbit and attitude of the satellite will cause a wave-like systematic error along the flight direction, and the orbit mode deviation eliminates the sensor's own distortion by converting the coordinates to the orbit coordinate system and fitting a two-dimensional trend surface.

[0024] Then, the radar penetration depth correction is carried out. The radar signal can penetrate the snow and ice layer, resulting in the measured height being the reflection surface height at a certain depth inside the snow and ice, rather than the true snow and ice surface height. The present application establishes a physical model of the penetration depth and the snow density, water content (loss tangent), frequency and other parameters based on the attenuation law of electromagnetic waves in snow and ice medium, and uses the C-band and X-band DEM data obtained by SRTM at the same time. Since the X-band penetration ability is much weaker than the C-band, the height difference of the two at the same position mainly reflects the difference in penetration depth. Through this observable difference value, the key parameters in the physical model are solved in reverse, so as to realize the spatialized penetration depth correction and correct the radar DEM to the true snow and ice surface.

[0025] Then, multi-temporal DEM difference analysis is performed. The activity of high loose body must be manifested as the change of ground elevation, such as the elevation decrease caused by glacier retreat and the elevation increase caused by collapse accumulation. The present application performs pixel-level subtraction calculation, i.e. difference, on DEMs in different periods after strict correction, directly obtains the elevation change amount, adopts robust regression method to fit the time series, and obtains reliable elevation annual change rate.

[0026] Finally, the activity of high loose body is evaluated. Single elevation change is not enough to comprehensively evaluate the activity, and multiple geomorphic activity signs need to be combined for comprehensive judgment. The present application combines the elevation change rate with other key geomorphic process indicators, forms a comprehensive activity index through weighted fusion, trains a logistic regression model based on historical disaster samples, converts the comprehensive index into occurrence probability, makes the evaluation result more statistically significant, and identifies the hot spot areas and time periods that are significantly aggregated in the time and space dimensions through scan statistics, and reveals the group occurrence rule and key risk period of disasters.

[0027] 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, rather than all the embodiments.

[0028] Figure 1 A flowchart of a high loose body historical activity detection method based on multi-source DEM data is shown, please refer to Figure 1 The method comprises the following steps: Step 1, division of basin unit: The basin unit naturally defined by the ridge line and the valley line is automatically divided based on the digital elevation model DEM data; the ice coverage rate of each basin unit is calculated according to the ice boundary data, and the candidate basin unit significantly affected by the glacier is screened according to the set ice coverage rate threshold.

[0029] In the high mountain glacier area, the formation, migration and accumulation of high loose body are strictly controlled by the terrain. The source area is usually located on the steep ice edge slope, the migration path is constrained by the valley system, and the accumulation area is distributed at the valley mouth or gentle place. Therefore, the basin unit naturally defined by the ridge line (watershed) and the valley line (converging line) is a closed and independent natural system that can fully reflect the complete disaster process of source-migration-accumulation, can truly reflect the migration path and convergence mode of surface matter, and ensure that the analysis result conforms to the geomorphic evolution rule and disaster dynamics principle. The automatic division of the basin unit can be realized by using the mature algorithm based on flow segmentation, such as the Chinese patent with the application publication number CN113850822A, which will not be described herein.

[0030] After obtaining the basin units, it is necessary to identify the candidate basin units significantly affected by glaciers. In this embodiment, the glacier coverage is used as a screening index. In technical implementation, first, the vector format glacier boundary data (such as RGI data) is converted into raster data with the same resolution as the DEM. For the grid cells across the basin boundary, the area weight method is used to accurately calculate the contribution: according to the actual coverage area proportion of the glacier in the grid, a weight value between 0 and 1 is given, and then the glacier coverage of the entire basin unit is accumulated.

[0031] The calculation formula of the glacier coverage is as follows: ; Wherein, represents the glacier coverage, and represents the row number and column number of the pixel in the grid data, represents the glacier coverage binary function of the pixel , if the pixel is located in the glacier area, then , otherwise 0, represents the basin unit, represents the total number of grid cells in the basin unit.

[0032] In practical application, the glacier coverage greater than 10% can be set as the screening standard, which can not only ensure that the screened basins are significantly affected by glaciers, such as moraine supply and ice melt water triggering, but also will not miss the areas with relatively small glaciers but still have potential disaster risks. For example, although the overall glacier coverage of some basins is not high, the local glacier retreat in the upstream or edge of the basin may still expose a large amount of loose material, forming a potential disaster source.

[0033] Step 2, multi-source DEM data processing: Obtain multi-source DEM data of the study area, including C-band and X-band DEM data of SRTM DEM, ALOS AW3D30 DEM data, and ICESat-1 and ICESat-2 laser altimetry data; preprocess the multi-source DEM data, including data format unification, coordinate system unification and spatial resolution unification, to obtain a standardized multi-temporal DEM data set.

[0034] Since no single remote sensing data source can meet all requirements of temporal coverage, spatial resolution, measurement accuracy and anti-interference ability at the same time, the embodiment adopts DEM data of multi-platform, multi-sensor and multi-temporal to complement each other. However, these original data have essential differences in imaging mechanism, spatial reference, acquisition time and accuracy characteristics, and cannot be directly used for high-precision differential comparison. Therefore, it is necessary to unify them to the same benchmark through strict preprocessing procedures to form a set of standard data sets with spatial and temporal comparability and consistent format, laying a foundation for subsequent precise correction and analysis. Through unified coordinate system, data format and spatial resolution, the embodiment eliminates systematic errors caused by inconsistent data standards and ensures that data of different periods and different sources are directly comparable.

[0035] The embodiment adopts three main types of DEM data: first, SRTM DEM data (acquired in February 2000), which is global elevation data acquired by a radar system carried by a space shuttle. The data contains C-band (5.7 GHz) and X-band (9.7 GHz) radar measurement results, with a spatial resolution of 30 meters. SRTM data is complete in coverage and stable in quality, and can be used as the ground elevation benchmark in 2000. Second, ALOS AW3D30 DEM data (acquired in 2006-2011), generated from optical stereo pairs of Japanese ALOS satellite, also with a spatial resolution of 30 meters. Unlike radar data, optical measurement directly obtains the surface elevation of ground objects, and there is no signal penetration problem in snow-covered areas, which makes it an important reference for correcting radar penetration errors. The third type is ICESat laser altimetry data, including ICESat-1 (2003-2009) and ICESat-2 (2019-2023) data. Although the spatial sampling of laser altimetry is relatively sparse (the laser point spacing of ICESat-1 is about 172 meters), its vertical accuracy is extremely high, up to centimeters, providing a reliable benchmark for the accuracy verification of other DEM data.

[0036] After obtaining the multi-source DEM data, the first step is to unify the data format and coordinate system. All original data are first converted to GeoTIFF format, using 32-bit floating-point type storage to ensure the accuracy of the elevation values. During the conversion process, the meta-information of the original data needs to be completely preserved, including projection parameters, geographic range, and invalid value identification, etc. For data containing multiple bands (such as the C-band and X-band of SRTM), they need to be extracted and stored as independent files respectively. After completing the format conversion, data quality inspection is carried out. Verify whether the data is complete and covers the study area, whether the elevation values are within a reasonable range (such as the high mountain area is usually between 2000-8000 meters), and eliminate obvious outliers. Subsequently, all data are projected to the UTM coordinate system corresponding to the study area. For example, if the center of the study area is located at 95 degrees east longitude, the UTM 47N projection zone is used, and the ellipsoid reference is unified to WGS84. This step ensures the strict consistency of different data sources in spatial position.

[0037] The second step is to unify the spatial resolution. In order to carry out accurate differential analysis, all data need to be resampled to the same spatial resolution. The resampling strategy is determined according to the relationship between the original resolution and the target resolution: for fine data with a resolution higher than 30 meters, bilinear interpolation method is used for down-sampling. This method determines the new elevation value by calculating the weighted average value of the four original pixels around the target pixel, which can smooth local details and maintain the overall characteristics of the terrain. For coarse data with a resolution lower than 30 meters (such as 500-meter resolution MODIS snow product), up-sampling is needed. This embodiment uses cubic convolution interpolation method: ; wherein, represents the output pixel elevation value to be solved, and represent the continuous coordinate values in the target coordinate system, represents the elevation value of the known input pixel, and represent the discrete matrix index of the input grid, represents the weight kernel function of cubic convolution interpolation, represents the normalized distance in the horizontal direction, represents the normalized distance in the vertical direction, and represent the weight values of the normalized distance in the horizontal and vertical directions, respectively; the weight kernel function is defined as follows: when , ; when , ; When , ; Wherein, represents a normalized distance, i.e. represents an interpolation coefficient.

[0038] The above steps effectively overcome the inherent defects of single data source in high mountain areas, such as cloud and fog obstruction, signal incoherence or sparse sampling, by integrating DEM data of three different principles of radar, optical and laser, form a robust data foundation of complementary advantages and mutual redundancy, and ensure that continuous and available elevation information can still be obtained in complex environment. Through strict format, coordinate and resolution three unification, the data from different sources are integrated into a set of highly standardized multi-temporal DEM data set. This eliminates the most basic systematic differences between data, enabling DEM data from SRTM, ALOS and ICESat to be directly and accurately compared at the pixel level, making it possible to detect sub-millimeter changes in the surface. In the process of resolution unification, the interpolation algorithm is optimized for different situations, especially the cubic convolution interpolation method considering 16 neighborhood pixels, which can better maintain the continuity and sharpness of the terrain features compared with the nearest neighbor method, significantly reducing the smoothing error introduced by the resampling process, and preserving the key terrain information for subsequent high-precision differential analysis.

[0039] Step 3, system bias correction: The standardized multi-temporal DEM data are sequentially subjected to geographic positioning bias correction, height distortion bias correction and satellite orbit pattern bias correction to obtain high-precision multi-temporal DEM data after system bias correction.

[0040] Due to differences in sensor physical characteristics, satellite orbit geometry, processing algorithms, etc., the errors of DEM data from different sources are not randomly distributed, but have obvious regularity of systematic bias. If these biases are not corrected, they will produce false signals far exceeding the true changes in differential analysis. The present application innovatively decomposes the complex total system bias into three relatively independent components with clear physical origin, and establishes the corresponding mathematical model for three-step progressive correction, thereby stripping the system error layer by layer and restoring the true terrain signal.

[0041] In this embodiment, the total system error can be decomposed into three components: the geolocation bias caused by horizontal position offset, the elevation distortion bias related to elevation, and the track pattern bias related to sensor track geometry. The three biases are different in magnitude and impact. The geolocation bias, which is sensitive to terrain and has the greatest impact, must be corrected first, then the elevation-related bias is processed, and finally the large-scale track-related bias is eliminated. Improper correction sequence will lead to error propagation and amplification.

[0042] wherein the geolocation bias is the most common system error, which is manifested as the horizontal displacement of the same terrain feature in different DEMs. In the terrain undulating area, this horizontal displacement will be converted into false elevation difference. Please refer to Figure 2 In Figure 2 , represents the amplitude of the horizontal offset vector, DEM represents the ground surface line of the reference period digital elevation model, represents the terrain slope, represents the false elevation gain generated in the ridge line area, represents the false elevation loss generated in the valley line area, represents the elevation difference of two period DEMs. This embodiment identifies and corrects this bias by analyzing the relationship between the elevation difference and the terrain slope, and the corresponding formula is as follows: ; wherein, represents the elevation difference of two period DEMs, represents the terrain slope, represents the amplitude of the horizontal offset vector, represents the direction of the offset vector, represents the terrain aspect, i.e. the orientation angle of the slope surface, represents the ratio of the vertical bias to the average slope tangent.

[0043] wherein, solved by least squares method, the formula is as follows: ; wherein, represents the sample point index, represents the total number of sample points participating in the calculation, represents the weight value of the th sample point, represents the elevation difference of the th sample point, represents the slope of the th sample point, represents the aspect of the th sample point.

[0044] Elevation distortion bias is the variation of DEM error with elevation. This phenomenon is particularly evident in high mountainous areas, mainly due to the sparse distribution of ground control points in high altitude areas. For example, SRTM data may have a systematic bias of 5-10 meters above 4000 meters. The present application uses a polynomial model to describe this elevation-dependent relationship: ; wherein, represents the elevation value of a pixel, represents the estimated value of the elevation distortion bias, i.e. the systematic elevation error at the elevation , and represent the minimum elevation and maximum elevation of the study area, respectively, represents the polynomial coefficient, represents the polynomial order, represents the polynomial order index.

[0045] wherein the polynomial order is determined using the Akaike Information Criterion, as follows: ; wherein, represents the AIC index, represents the total number of stable area elevation difference sample points used to fit the elevation distortion bias model, for example, after stratification by elevation (e.g. 100 meter intervals), the total number of sample points used to calculate the average elevation difference in the stable area (e.g. bedrock) of each layer, is the residual variance, representing the sum of squared errors of the model fitting.

[0046] In practical applications, the steps for implementing the elevation distortion bias correction are as follows: ① stratify the study area by 100 meter elevation intervals and calculate the average elevation difference in the stable area of each layer. ② Fit the polynomial coefficients using weighted least squares. The weight is proportional to the number of effective sample pixels in each layer, ensuring that the data-rich layer contributes more. ③ Use the Akaike Information Criterion to adaptively determine the optimal polynomial order, selecting the polynomial order with the smallest AIC. ④ Use cross-validation (e.g. 70% training set, 30% validation set) to prevent model overfitting, and the number of stable area pixels participating in the statistics for each elevation layer should not be less than 1000.

[0047] Orbital pattern bias is related to satellite orbit characteristics, resulting from specific sensor geometric characteristics (such as periodic changes in attitude, geometric distortion, and directional differences in atmospheric refraction), and is manifested as different systematic error patterns along the satellite flight orbit direction and perpendicular to the orbit (cross-track) direction. For example, SRTM DEM often exhibits a long-wave feature with a wavelength of about 50-100 kilometers and an amplitude of 2-3 meters along the orbital direction. To correct such bias, the present embodiment establishes a two-dimensional polynomial trend surface model in the satellite orbit coordinate system: ; wherein, represents the orbit related system bias, i.e. the elevation error estimate at the orbit coordinate system coordinate , represents the coordinate of the pixel along the orbit direction, represents the coordinate along the cross-orbit direction, represents the trend surface coefficients, and represents the trend surface order, and represents the polynomial order index.

[0048] wherein, the orbit coordinate system coordinate is obtained by a rotation transformation, the formula is as follows: ; wherein, and represents the coordinate of the pixel under the original map projection coordinate system, represents the orbit direction angle.

[0049] In practical application, the orbit mode bias correction implementation steps are as follows: ①convert the geographic coordinates of all pixels in the study area into the orbit coordinate system; ②fit the trend surface coefficients by using the elevation difference data of the stable terrain area in the study area; ③apply the fitted trend surface model to the entire DEM data to deduct the bias, and use the iterative method for correction. Calculate the standard deviation of the residual error after each iteration; ④stop iteration when the standard deviation of the residual error is less than the preset threshold (usually 0.5 meters) or reaches the maximum iteration number (usually 5 times).

[0050] Through the above three steps of correction, the errors of different physical origins are modeled and eliminated respectively, the systematic bias between multi-source DEM data is separated and suppressed from the root, and the incompleteness of mixed correction is avoided. The systematic bias between different source DEMs can be reduced to sub-meter level, providing a reliable basis for subsequent high-precision differential analysis.

[0051] Step 4, radar penetration depth correction: Based on the propagation characteristics of electromagnetic waves in ice and snow medium, a penetration depth physical model is established; the elevation difference between C-band and X-band DEM data in the high-precision multi-temporal DEM data is used to calculate the penetration depth difference, and the parameter optimization of the penetration depth physical model is carried out; the optimized penetration depth physical model is applied to the ice-covered area in the candidate watershed unit, and the C-band DEM data is spatialized and corrected to obtain the final DEM data after penetration depth correction.

[0052] The core principle of this step is the physical interaction of radar wave with snow and ice medium and dual-frequency differential detection. Unlike optical sensors, radar signals have a certain penetration ability for snow and ice, and their return signals come from the inside of snow and ice rather than the real ground surface. This results in a systematic lower elevation measured by radar DEM than the real ground surface elevation, and the difference is the penetration depth. The traditional method considers the penetration depth as a constant, ignoring its spatial variability with snow and ice physical properties (density, water content) and terrain (elevation), resulting in limited correction accuracy.

[0053] This embodiment is based on the propagation characteristics of electromagnetic waves in snow and ice medium, and establishes a physical model of penetration depth: ; Wherein, represents the radar penetration depth, i.e. the difference between the radar DEM and the true ground surface elevation, , snow density , and terrain (elevation) under the condition of represents the snow surface elevation measured by radar altimetry, represents the snow density, represents the sea level reference penetration depth, C band takes 4.5 meters, X band takes 0.5 meters, represents the radar frequency, C band is 5.7 GHz, X band is 9.7 GHz, represents the speed of light, represents the loss tangent, dry snow takes 0.001, wet snow takes 0.01, represents the snow thickness, which can be estimated by meteorological model or measured, represents the elevation gradient coefficient, represents the reference elevation, usually the average elevation of the study area, represents the constant term, which is calibrated by stable area, represents the natural exponential function.

[0054] The accurate calibration of the parameters of the penetration depth physical model relies on establishing a reliable reference benchmark in stable regions (i.e., regions without surface elevation changes and without ice and snow coverage). The elevation differences in these regions should be mainly attributed to systematic errors of the radar penetration effect itself rather than real surface changes. In practical applications, a multi-criteria comprehensive discrimination method can be used to identify stable regions: ① terrain stability, the slope is less than 40°, avoiding DEM errors caused by terrain shadows and overlap effects; ② no glacier coverage, based on multi-temporal (at least 5 years of time span) high-resolution optical image analysis, confirming that there is no glacier in the region; ③ no seasonal snow, the normalized difference snow index (NDSI) is less than 0.4 (calculated based on clear sky optical images), excluding the influence of seasonal snow; ④ land cover type, the land cover type is bare rock or sparse vegetation (verified by GlobeLand30 and other data). The regions that meet all the above conditions at the same time are marked as stable regions, which are used for subsequent calibration and verification of the parameters of the penetration depth model.

[0055] In addition, this step also utilizes the C-band and X-band data simultaneously acquired by SRTM. Since the X-band (9.7 GHz) has a higher frequency, the penetration depth is much smaller than the C-band (5.7 GHz). Therefore, the elevation difference of the two at the same location mainly reflects the difference in penetration depth, which contains the information of snow physical parameters such as snow density and thickness. By establishing the functional relationship between the difference and the snow parameters, the penetration depth correction value can be estimated in the areas where there is a lack of in-situ observations.

[0056] Specifically, first, the penetration depth difference is calculated: ; wherein, represents the penetration depth difference of the C-band and X-band DEM data, and represent the elevation values of the C-band and X-band DEM data, respectively, and represent the penetration depths of the C-band and X-band DEM data, respectively.

[0057] Then, the parameter optimization of the penetration depth physical model is realized by minimizing the deviation of the predicted penetration depth difference from the actual observed elevation in the stable regions, and regularization is introduced to prevent overfitting. The parameter optimization formula is as follows: ; wherein, represents the set of stable region pixels, represents the set of parameters to be optimized, represents the observed elevation difference of the C-band and X-band DEM data, represents the predicted penetration depth difference based on the parameter set ​denotes a regularization coefficient, denotes a regularization term.

[0058] In practical applications, the optimization process adopts a three-stage strategy to ensure robustness and accuracy. The specific steps are as follows: First, perform a rough search, and perform a grid scan within the physically reasonable range of parameters to quickly locate the approximate area of the optimal solution. Then, perform fine optimization, and use the Levenberg-Marquardt algorithm to start from the coarse search result and iteratively solve until convergence. Finally, perform uncertainty evaluation, and evaluate the stability and confidence interval of parameter estimation through Bootstrap resampling (1000 times).

[0059] Substitute the optimal parameter set obtained by optimization into the penetration depth physical model to obtain the spatialized penetration depth model. Apply the spatialized penetration depth model to the ice-covered area in the candidate drainage unit to perform per-pixel correction of the C-band DEM, and obtain the final DEM data after penetration depth correction, whose elevation value represents the true snow and ice surface elevation.

[0060] Through the above correction process, the systematic deviation caused by radar penetration effect can be effectively eliminated, and the elevation value truly reflects the position of the ground surface, avoiding misjudgment of the penetration effect as ground surface activity, and fundamentally ensuring the accuracy and reliability of the final active evaluation result. Practice shows that the elevation accuracy of the corrected DEM in the glacier area can be improved from meter level to sub-meter level, providing a reliable data basis for accurately evaluating glacier change and loose body activity.

[0061] Step 5, multi-temporal DEM difference analysis: Perform difference calculation on the final DEM data of different time periods, and use a formula based on robust regression to fit the elevation change trend of each drainage unit to obtain the spatial and temporal distribution data of the elevation change rate reflecting the dynamic of the ground surface.

[0062] After completing the system preprocessing, bias correction and radar penetration depth correction of multi-source DEM data, the embodiment reveals the dynamic change of ground elevation through multi-temporal DEM difference analysis. This step is the core technical means for monitoring glacier ablation, collapse and identifying the historical activity of high loose bodies. The basic principle is to compare the ground elevation models obtained in different periods after strict correction, and calculate the elevation change at a specific location. In the high mountain glacier area, such elevation change is mainly caused by key geomorphic processes, such as ice loss caused by glacier retreat, ice-rock collapse events, and erosion and accumulation caused by debris flow and other activities. Therefore, accurately quantifying the spatio-temporal change pattern of ground elevation provides the most direct quantitative evidence for evaluating the historical activity of high loose bodies.

[0063] The traditional two-period DEM simple difference method is susceptible to accidental errors and outliers, affecting the reliability of the change trend estimation. To overcome this limitation, the embodiment uses a robust regression method based on multi-temporal (not less than 3 periods) DEM data to estimate the annual change rate of elevation. This method can not only obtain a more robust estimate of the annual average change rate by fitting the overall trend of elevation observations over time, but also effectively identify nonlinear dynamic processes, such as accelerated glacier ablation or sudden collapse events. The core is to introduce the robust estimation theory, so that the fitted change trend is not sensitive to outliers in the data (such as errors caused by temporary snow, short-term cloud shadows), so as to obtain more reliable and more realistic physical process of elevation change rate.

[0064] In this embodiment, the formula of the robust regression is as follows: ; Wherein, and denote the annual change rate of elevation, denote the number of time series observations, denote the elevation observation value of the period, denote the observation time of the period, denote the reference elevation, denote the reference time point, denote the standard deviation of elevation observation, denote the Huber loss function.

[0065] Compared with the simple two-period difference, the robust regression method based on multi-temporal can effectively suppress the interference of abnormal elevation values caused by residual noise, temporary snow, vegetation, etc., and avoid its excessive influence on trend estimation. It can provide a change trend that is more consistent with the actual physical process, for example, it can identify the accelerated ablation trend of glaciers or the mutation signal of collapse events, without being misled by individual abnormal observations. This provides a scientific basis for accurately distinguishing different activity intensity of loose body regions in the future, and for risk classification and precise prevention and control.

[0066] Step 6, high loose body activity assessment: Based on the spatial and temporal distribution data of the elevation change rate, the historical activity comprehensive index of high loose body is calculated in the candidate watershed unit, the activity probability is estimated, and the activity hotspot area and period are identified, to obtain the quantitative assessment result of the historical activity of high loose body.

[0067] The final goal of the embodiment is to quantitatively evaluate the historical activity of high loose bodies in the high mountain glacier area by comprehensively analyzing the multi-temporal DEM difference results and auxiliary information. The evaluation not only needs to identify the areas where the ground has significant elevation changes, but also needs to deeply understand the geomorphological meaning of these changes to judge whether they indicate potential disaster risks, so as to provide a scientific basis for risk management. Therefore, the embodiment constructs a multi-level evaluation system, from single factor change analysis, multi-factor comprehensive evaluation to spatiotemporal pattern recognition, to gradually improve the reliability and practicality of the evaluation results.

[0068] Firstly, in order to fully capture the active characteristics of high loose bodies, the embodiment establishes a comprehensive evaluation model based on the evidence weight method, which is used to calculate the historical activity comprehensive index of high loose bodies, as follows: ; Among them, represents the historical activity comprehensive index of high loose bodies, represents the evaluation factor serial number, represents the total number of evaluation factors, represents the weight of the th evaluation factor, represents the factor value of the th evaluation factor, represents the confidence of the th evaluation factor.

[0069] In the above formula, different evaluation factors have different contributions to the activity, and the overall influence is quantified by weighted synthesis. While ensuring the interpretability of the model, the weights of each evaluation factor can be flexibly adjusted according to the regional characteristics. Each evaluation factor is standardized to the interval to ensure that factors of different dimensions can be compared.

[0070] The recommended core evaluation factors and their initial weights and feature threshold values are as follows: ① elevation change rate factor , reflecting the erosion / deposition intensity of loose bodies, weight 0.4, threshold > 2m / year for high activity; ② ice lake expansion rate factor , indicating the increase of upstream glacier meltwater and the risk of outburst, weight 0.2, threshold > 20% / 10 years for high risk; ③ glacier retreat factor , reflecting the speed of ice loss and exposed loose material source, weight 0.2, end retreat > 200m / 10 years for significant; ④ gully erosion factor , indicating the gully bed incision or lateral erosion intensity of activities such as debris flow, weight 0.2, width increase > 50% for active. Confidence According to the data coverage, the phase consistency and the processing precision evaluation assignment (such as: high-quality data =1.0, medium =0.7, low-quality =0.4).

[0071] Secondly, in order to further quantify the possibility of disaster occurrence, the embodiment constructs a logistic regression model based on historical disaster samples to predict the activity probability, that is, the activity is regarded as a probability event, a statistical model is established based on historical disaster samples, and the comprehensive index is converted into a more intuitive occurrence probability, and the formula is as follows: ; wherein, represents the activity probability, represents the intercept term, which is obtained by training the historical disaster samples, represents the regression coefficient of the first evaluation factor, reflecting the importance of each factor, represents the standardized factor value of the first evaluation factor, using z-score standardization, represents the natural exponential function.

[0072] Based on the activity probability, the grade can be divided, for example: the activity probability is greater than or equal to 0.7, which is determined as high activity, 0.3 is less than or equal to the activity probability and is less than 0.7, which is determined as medium activity, and the activity probability is less than 0.3, which is determined as low activity.

[0073] Finally, the occurrence of disaster is not random in space and time, but shows aggregation. By identifying these hot spots, the group occurrence rule and key risk period of disaster can be revealed. In order to identify the aggregation mode (hot area and time period) of disaster activity in space and time, the embodiment adopts the space-time scanning statistic to identify the most possible cluster and secondary cluster with statistical significance, and the space-time cluster analysis model is as follows: ; wherein, represents the space-time cluster statistic, represents the spatial window, represents the time window, represents the number of high activity points in the window, represents the expected value based on Poisson distribution, represents the total number of points, represents the window likelihood function value, represents the total likelihood function value.

[0074] The technical implementation of the spatio-temporal clustering analysis includes the following key steps: ① Construct a cylinder with a spatial circle as the bottom surface and a time span as the height, and move the position and size of the cylinder to cover all possible spatio-temporal combinations. ② Calculate the likelihood ratio of each scanning window . ③ Construct the empirical distribution of the likelihood ratio statistic through Monte Carlo simulation (e.g., 999 random rearrangements of spatio-temporal positions), and calculate the empirical value E corresponding to the observed maximum likelihood ratio. ④ Identify spatio-temporal clusters with statistical significance (E < 0.05), distinguish between "most likely clusters" (with the largest likelihood ratio value) and "secondary clusters", and generate a spatio-temporal risk hotspot distribution map.

[0075] This step constructs a comprehensive and objective quantitative evaluation system by integrating multiple key geomorphic process indicators, significantly improving the scientificity and reliability of the evaluation results. Not only does it include quantitative comprehensive indexes, but also risk probabilities that are easy to understand and explicit spatio-temporal risk hotspot distribution maps. This multi-dimensional output greatly facilitates the interpretation and application of the results, providing intuitive decision-making basis for managers and decision-makers with different knowledge backgrounds.

[0076] To verify the effectiveness of the method, the southeastern Tibetan Plateau region (94°-98°E, 27°-30°N) was selected as the test area to verify the method provided in this embodiment. This region is located at the junction of the eastern section of the Himalayas and the Hengduan Mountains, with an elevation span of 1500-7500 meters and a large number of oceanic glaciers. The annual precipitation in this region is 800-2000 millimeters, and the glacier equilibrium line height is about 4800-5200 meters. Historical records show that there have been more than 200 mudslides in this region, 70% of which are related to glacier activity, providing an ideal natural experiment field for method verification. The specific process includes the following steps: (1) Basin unit division After filling the concave points of the 30-meter SRTM DEM, a total of 3847 points were processed, with an average depth of 0.3 meters. The D8 single flow direction algorithm was used to calculate the cumulative flow, which took about 35 minutes. The cumulative flow threshold was set to 1000 (catchment area 0.9 km²) to extract the river network, obtaining a total length of 12,450 km of river network system. Combined with 5-meter resolution Google Earth images, manual correction was performed on 1238 boundaries. Calculate the glacier coverage of each basin unit, and set the screening standard of glacier coverage >10%, finally obtain 661 candidate basin units (area range 20.6-651.8 km², average 87.3 km²; average slope 28.7°±8.4°, average elevation 4285±752 meters, average glacier coverage 31.2%).

[0077] (2) Multi-source DEM data processing The following core data were acquired and pre-processed: ① SRTM DEM (February 2000): C-band and X-band, resolution 30 m, vertical accuracy ± 16 m. ② ALOS AW3D30 DEM (2006-2011): resolution 30 m, vertical accuracy ± 5 m. ③ ICESat laser altimetry data: ICESat-1 (2003-2009) and ICESat-2 ATL06 (2019-2023). ④ Auxiliary data: RGI 6.0 glacier inventory, GlobeLand30 land cover data, MODIS snow cover product. Pre-processing includes: uniform conversion to GeoTIFF format (32-bit floating point) using GDAL tools; projection conversion to UTM 47N (WGS84 ellipsoid); spatial resolution unification using cubic convolution interpolation method.

[0078] (3) Systematic bias correction ① Geolocation bias correction: 23,456 pixels in stable bedrock area with slope 15°-45° were selected as control points. Iterative solution obtained SRTM horizontal offset 12.3 m / 47°, ALOS offset 8.7 m / 135°; SRTM vertical bias coefficient -1.2, ALOS 0.8. After correction, the standard deviation of elevation difference in stable area decreased from 18.6 m to 6.2 m, improved by 66.7%.

[0079] ② Elevation distortion bias correction: 60 layers according to 100 m elevation. AIC criterion determines to use 3-order polynomial, fitting polynomial coefficients: =-2.3, =8.7, =-12.4, =5.1. After correction, the mean bias of each elevation band is controlled within ± 0.5 m. In southeast Tibet, the elevation-related bias between SRTM DEM and ALOS AW3D30 DEM is significant, and each 1000 m of altitude change corresponds to about 40 m of systematic bias.

[0080] ③ Orbit pattern bias correction: SRTM orbit angle 10.2°, ALOS -8.4°. Fitting trend surface coefficients =0.8, =-0.003, =0.001, =0.00002. After correction, the standard deviation of residual error decreased to 0.48 m.

[0081] The results of the systematic bias correction of this embodiment are shown in Figure 3 . Figure 3The left side is a diagram of the standard deviation of height deviation and its systematic deviation trend. The columnar filling represents the change of the standard deviation of height deviation at different processing stages, and the dashed curve represents the systematic deviation trend of each 1000 meters of altitude. The values shown in the figure (including 18.6, 6.2, 2.4, and 0.48) are the corresponding height deviation standard deviation values at each stage, which shows that the height accuracy gradually improves as the correction steps proceed.

[0082] Figure 3 The right side is a deviation comparison column chart at different correction stages. From left to right, they are the original deviation, the deviation after geographic positioning correction, the deviation after height distortion correction, and the deviation after orbit pattern correction. The column height gradually decreases, which directly shows that the multi-stage correction method provided by the present application can effectively reduce the height deviation and gradually improve the accuracy of the measurement results.

[0083] (4) Radar penetration depth correction Combined with the field snow pit measurement (snow density 300-450 kg / m³, loss tangent: dry snow 0.0008-0.0012, wet snow 0.008-0.015, snow thickness 0.5-8.0 meters), the penetration depth physical model is applied. The calculation shows that at an altitude of 4000 meters, the C-band penetration is 4.3±0.7 meters, and the X-band penetration is 0.5±0.2 meters; at an altitude of 5000 meters, the C-band penetration is 4.8±0.8 meters, and the X-band penetration is 0.6±0.2 meters; at an altitude of 6000 meters, the C-band penetration is 5.2±0.9 meters, and the X-band penetration is 0.7±0.3 meters. The C-X band height difference optimization model parameters are calculated using 28450 glacier pixels. After correction, the glacier area height deviation is significantly improved from -4.7±3.2 meters to -0.2±1.1 meters. Table 1 shows the radar penetration depth correction parameters at different altitudes.

[0084] Table 1 Radar penetration depth correction parameters at different altitudes

[0085] (5) Multi-temporal DEM difference analysis Select 29 ICESat data sufficient basins: as shown in Figure 4 The height change curves of basins A, B, C, and D from 2000 to 2023 based on SRTM, ICESAT-1, and ICESAT-2 data. Figure 4 The specific change rates of a plurality of typical basin units are shown in the form of a list: the change rate of basin A is -2.3±0.4 meters / year, the change rate of basin B is -1.8±0.3 meters / year, the change rate of basin C is -3.1±0.5 meters / year, and the change rate of basin D is -1.2±0.2 meters / year. Figure 4The uncertainty ranges of the values are also marked in the figure. Huber regression converges after 15 iterations. The significance threshold is set as 2 times the standard deviation of the difference between SRTM and ALOS (11.4 m; 22.8 m) and SRTM and ICESat (10.0 m; 20.0 m). Among the 661 basins, 412 (62.3%) are detected with significant elevation changes: -18.7 ± 8.2 m (2000-2010) in the glacierized areas and -0.3 ± 3.1 m in the non-glacierized areas, which verifies the effectiveness of the correction. The results show that the average elevation of these areas gradually decreases over time, indicating that the glacier-covered areas are continuously melting, increasing the risk of high loose body disasters.

[0086] (6) High loose body activity assessment ① Multi-factor comprehensive assessment: In this embodiment, the AHP method is used to determine the weight: the elevation change rate is 0.42, the ice lake expansion is 0.28, the glacier retreat is 0.19, and the channel erosion is 0.11. Taking basin C as an example, the factor value of each evaluation factor is: = 0.85 (change rate -25.3 m / 10 years), = 0.92 (ice lake 0.05→0.15 km²), = 0.73 (end retreat 385 m), = 0.67 (gully width 15→25 m), and the comprehensive index is 0.81 (high activity).

[0087] ② Activity probability estimation: Based on 87 historical disaster points (100 m buffer) and 261 stable points, a logistic regression model is trained. The intercept term and regression coefficients are: = -1.73, = 2.14, = 1.86, = 1.52, = 0.98. The model AUC = 0.87, and the optimal threshold is 0.65. 10-fold cross-validation: average accuracy 84.3%, recall rate 78.5%.

[0088] ③ Spatio-temporal clustering analysis: Set the spatial radius to 1-10 km and the time window to 1-5 years. The main cluster (center: 96.5°E, 28.3°N; radius 6 km; period 2005-2010; relative risk RR = 3.7; Monte Carlo 999 times, E = 0.003) and 3 secondary significant clusters (E < 0.05) are identified.

[0089] The multi-factor comprehensive assessment results of this embodiment are shown in Figure 5 . Figure 5The main factors affecting geological activity and their weights in the evaluation system are listed, including: terrain gradient (as a supplementary indicator), gully erosion (weight: 0.11), glacier retreat (weight: 0.19). Figure 5 The comprehensive evaluation results of three typical basins are also shown: Basin C is evaluated as "high activity" with a comprehensive index of 0.81, Basin A is evaluated as "medium activity" with a comprehensive index of 0.58, and Basin D is evaluated as "low activity" with a comprehensive index of 0.35. Figure 5 The value range of the comprehensive index (1.0 to 0.2) is also marked in the form of a gradient strip, which corresponds to different levels of geological activity. Figure 5 The effectiveness and application results of the technical scheme of the present application in the aspects of comprehensive multi-source index, weight quantification and generation of geological activity index are described in general, which embodies the systematization and discrimination of the evaluation method.

[0090] (7) Result verification and technical effect The comprehensive evaluation finally identifies 78 historical active points: 25 high activity points (30.3%, concentrated in the glacier terminal / ice lake periphery), 35 medium activity points (46.1%, glacier retreat area), and 18 low activity points (23.6%, sporadic distribution). Compared with the actual disaster records from 2010 to 2020: 49 out of 60 disasters were successfully identified, with an identification accuracy of 81.6%, a missed detection rate of 18.4%, and a false alarm rate of 15.8%, indicating high reliability of the method. Typical case (Basin A): From 2000 to 2010, the glacier area decreased from 19.11 km² to 6.99 km² (-63.4%), 3 new ice lakes were formed in the main gully (total area 0.28 km²), and the elevation showed an accelerating downward trend, indicating that the disaster risk will continue to rise in the future. The performance evaluation and verification effect of the active detection model is shown in Figures 6 to 8 .

[0091] Among them, Figure 6 The ROC curve of the active detection model is shown, which is the receiver operating characteristic curve. The ROC curve represents the performance of the model at different decision thresholds, where the horizontal axis is the false positive rate and the vertical axis is the true positive rate. The marked point on the curve indicates the optimal decision threshold (0.65) determined after model optimization. This threshold point corresponds to a combination of a high true positive rate and a low false positive rate, which reflects the good classification accuracy and reliability of the model in identifying geological activity.

[0092] Figure 7A detection result confusion matrix of the activity detection model is shown, and the confusion matrix compares the model prediction results with the actual observation results in tabular form, specifically including: the number of instances that are actually active and correctly predicted by the model as active is 49, accounting for 81.7% of the total number of actual active instances; the number of instances that are actually active but incorrectly predicted by the model as stable is 11, accounting for 18.3% of the total number of actual active instances; the number of instances that are actually stable but incorrectly predicted by the model as active is 12, accounting for 2.0% of the total number of actual stable instances; and the number of instances that are actually stable and correctly predicted by the model as stable is 589, accounting for 98.0% of the total number of actual stable instances. The activity detection model constructed by the application is intuitively verified to have high accuracy and high reliability.

[0093] Figure 8 A spatio-temporal clustering analysis result schematic diagram is shown, which shows the clustering structure obtained by performing spatio-temporal feature analysis on the target region using the method of the application, including a main cluster and multiple secondary clusters (including secondary cluster 1, secondary cluster 2, and secondary cluster 3) to represent the hierarchical features of the data in spatial distribution and temporal evolution.

[0094] In summary, the high-position loose body historical activity detection method based on multi-source DEM data provided in the embodiment realizes accurate and quantitative detection of the historical activity of high-position loose bodies by constructing a systematic data processing chain and a multi-level evaluation model. Through multi-source data fusion, physical model correction, intelligent algorithm analysis, and spatio-temporal pattern recognition, a complete high-position loose body historical activity detection technology system is formed, which has important application value in the field of early identification and risk assessment of geological disasters in glacial landforms, and provides reliable technical support for engineering planning and disaster prevention and mitigation in high mountain and canyon areas.

Claims

1. A method for detecting historical activity of high loose body based on multi-source DEM data, characterized in that, The method comprises: The method comprises: Obtain multi-source DEM data of the study area, including C-band and X-band DEM data of SRTM DEM, ALOS AW3D30 DEM data, and ICESat-1 and ICESat-2 laser altimetry data; pre-process the multi-source DEM data, including data format unification, coordinate system unification, and spatial resolution unification, to obtain a standardized multi-temporal DEM data set; The standardized multi-temporal DEM data is sequentially subjected to geographic positioning bias correction, height distortion bias correction, and satellite orbit pattern bias correction to obtain high-precision multi-temporal DEM data after system bias correction; Based on the propagation characteristics of electromagnetic waves in ice and snow medium, a penetration depth physical model is established; the height difference between C-band and X-band DEM data in the high-precision multi-temporal DEM data is used to calculate the penetration depth difference, and the parameter optimization of the penetration depth physical model is performed; the optimized penetration depth physical model is applied to the ice-covered area in the candidate watershed unit to correct the spatialization of the C-band DEM data, and the final DEM data after penetration depth correction is obtained; Differential calculation is performed on the final DEM data at different time periods, and a formula based on robust regression is used to fit the height change trend of each watershed unit to obtain height change rate spatiotemporal distribution data reflecting the surface dynamics; Based on the height change rate spatiotemporal distribution data, the historical activity comprehensive index of high-level loose body is calculated in the candidate watershed unit, the activity probability is estimated, and the activity hotspot area and time period are identified to obtain the quantitative evaluation result of the historical activity of high-level loose body.

2. The method for detecting high-situated loose body historical activity based on multi-source DEM data according to claim 1, characterized in that, The formula for calculating the ice coverage rate is as follows: ； wherein, represents the glacier coverage ratio, and represents the row number and the column number of the pixel in the grid data within the basin unit, represents the pixel a glacier coverage binary function of the pixel, if the pixel is located within the glacier area, then , otherwise 0, represents the basin unit, represents the total number of grid units within the basin unit.

3. The method for detecting high-situated loose material historical activity based on multi-source DEM data according to claim 1, characterized in that, The spatial resolution unification adopts a cubic convolution interpolation method, and the formula is as follows: ； wherein, represents the output pixel elevation value to be solved, and represents the continuous coordinate value in the target coordinate system, represents the elevation value of the known input pixel, and represents the discrete matrix index of the input raster, represents the weight kernel function of the cubic convolution interpolation, represents the normalized distance in the horizontal direction, represents the normalized distance in the vertical direction, and respectively represent the weight values of the normalized distances in the horizontal direction and the vertical direction. The weight kernel function is defined as follows: When Time, ; When time, ; When time, ; wherein, represents the normalized distance, i.e. the distance of the center point of the target pixel to the center point of the adjacent original pixel, represents the interpolation coefficient.

4. The method for detecting high-situated loose material historical activity based on multi-source DEM data according to claim 1, characterized in that, The formula corresponding to the geographic positioning bias correction is as follows: ； wherein, represents the elevation difference between two DEMs, represents the terrain slope, represents the magnitude of the horizontal offset vector, represents the direction of the offset vector, represents the terrain aspect, i.e. the orientation angle of the slope surface, represents the ratio of the vertical offset to the average slope tangent; Solving by least squares The formula is as follows: ； in, Indicates the sample point index. This represents the total number of sample points involved in the calculation. Indicates the first The weight values ​​of each sample point Indicates the first Elevation difference of each sample point Indicates the first The slope of each sample point Indicates the first The slope aspect of each sample point.

5. The method for high-situated loose body historical activity detection based on multi-source DEM data according to claim 1, characterized in that, The formula corresponding to the height distortion bias correction is as follows: ； wherein, represents an image element elevation value, represents an elevation twist bias estimate, i.e. systematic elevation error at elevation represents an elevation twist bias estimate, i.e. systematic elevation error at elevation and respectively represent the minimum and maximum elevation of the study area, represents a polynomial coefficient, represents a polynomial order, represents a polynomial order index; The Akaike information criterion is used to determine the polynomial order, and the formula is as follows: ； wherein, represents the AIC index, represents the total number of stable area elevation difference sample points for fitting the elevation distortion bias model, is the residual variance, which represents the error sum of squares of model fitting.

6. The method for high-situated loose body historical activity detection based on multi-source DEM data according to claim 1, characterized in that, The formula corresponding to the satellite orbit pattern bias correction is as follows: ； wherein denotes the orbit related system bias, i.e. the elevation error estimate at the orbit coordinate system coordinate , denotes the pixel coordinate along the orbit direction, denotes the pixel coordinate along the cross-orbit direction, denotes the trend surface coefficients, and denotes the trend surface order, and denotes the polynomial order index; Orbital coordinate system coordinates Obtained by a rotation transformation, the formula being as follows: ； wherein, and denotes the coordinates of the pixel in the original map projection coordinate system, denotes the track direction angle.

7. The method for high-situated loose body historical activity detection based on multi-source DEM data according to claim 1, characterized in that, The penetration depth physical model is as follows: ； wherein represents the radar penetration depth, i.e. the penetration depth at an altitude , the snow density under the condition of the penetration depth, represents the snow surface altitude, represents the snow density, represents the sea level reference penetration depth, represents the radar frequency, represents the speed of light, represents the loss tangent, represents the snow depth, represents the altitude gradient coefficient, represents the reference altitude, represents the constant term, represents the natural exponential function; The formula for calculating the penetration depth difference is as follows: ； wherein, represents the difference in penetration depth of the C-band and X-band DEM data, and respectively represent the elevation values of the C-band and X-band DEM data, and respectively represent the penetration depths of the C-band and X-band DEM data; The formula for parameter optimization of the penetration depth physical model is as follows: ； wherein, denotes a stable region pixel set, denotes a set of parameters to be optimized, denotes an observed elevation difference of C-band and X-band DEM data, denotes a predicted penetration depth difference based on the set of parameters denotes a predicted penetration depth difference, denotes a regularization coefficient, denotes a regularization term.

8. The method for high-situated loose body historical activity detection based on multi-source DEM data according to claim 1, characterized in that, The formula of the robust regression is as follows: ； in, and This represents the annual rate of change in elevation. Represents the number of time series observations. Indicates the first Periodic elevation observations, Indicates the first Observation period Indicates reference elevation. Indicates a reference time point. Indicates the standard deviation of elevation observations. This represents the Huber loss function.

9. The method for high-situated loose body historical activity detection based on multi-source DEM data according to claim 1, characterized in that, The formula for calculating the historical activity comprehensive index of high-level loose body is as follows: ； wherein, denotes a high-rank loose body history activity comprehensive index, denotes an evaluation factor serial number, denotes a total number of evaluation factors, denotes a weight of the evaluation factor, denotes a factor value of the evaluation factor, denotes a confidence of the evaluation factor; The formula for the activity probability estimation is as follows: ； wherein, represents the activity probability, represents the intercept term, represents the regression coefficient of the th evaluation factor, represents the standardized factor value of the th evaluation factor, represents the natural exponential function.

10. The method for high-situated loose body historical activity detection based on multi-source DEM data according to claim 1, characterized in that, The identification of the activity hotspot area and time period adopts a spatiotemporal clustering analysis model, and the formula of the spatiotemporal clustering analysis model is as follows: ； wherein, denotes a spatio-temporal clustering statistic, denotes a spatial window, denotes a temporal window, denotes the number of highly active points within a window, denotes an expected value based on a Poisson distribution, is a total number of points, denotes a window likelihood function value, denotes a total likelihood function value.

Images