Method for detecting historical activity of high-position loose body based on multi-source DEM data

Abstract

The present application relates to the technical field of radar remote sensing, and discloses a high-position loose body historical activity detection method based on multi-source DEM data, aiming to solve the problem that the prior art cannot accurately and quantitatively evaluate the historical activity of high-position loose bodies in a glacier-covered area, and the scheme mainly comprises the following steps: automatically dividing a watershed unit and screening a glacier-affected area based on DEM data; eliminating systematic errors by fusing multi-source DEM data, correcting geographical positioning deviation, correcting height distortion deviation, and correcting orbit mode deviation; establishing a penetration depth physical model, using C / X band height difference optimization parameters to correct the radar penetration depth in the glacier area; calculating the height change rate by using robust regression; calculating the activity comprehensive index by comprehensively considering multiple factors, and identifying hot regions and time periods by probability estimation and spatiotemporal clustering. The present application realizes automatic and high-precision quantitative evaluation of the historical activity of high-position loose bodies, and is particularly suitable for geological disaster warning and risk assessment 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 DEM data obtained by different sensors at different times have differences in coordinate system, spatial resolution, height datum, etc., and the systematic deviation seriously affects 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. 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 mainly 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 realize high-precision surface change detection. SUMMARY

[0007] The present application aims to solve the technical problem that the existing technology cannot realize 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 scheme adopted by the present application to solve the above technical problems is:

[0009] The method for detecting the historical activity of high loose bodies based on multi-source DEM data comprises:

[0010] The method for detecting the historical activity of high loose bodies based on multi-source DEM data comprises:

[0011] 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;

[0012] 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 subjected to system bias correction;

[0013] 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 and snow covered area in the candidate basin unit, and the C-band DEM data is spatially corrected to obtain the final DEM data subjected to penetration depth correction;

[0014] The final DEM data at different time periods is subjected to difference calculation, and the height change trend of each basin unit is fitted by a formula based on robust regression to obtain height change rate spatio-temporal distribution data reflecting the surface dynamics;

[0015] Based on the height change rate spatio-temporal distribution data, the historical activity comprehensive index 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.

[0016] Further, the calculation formula of the ice coverage rate is as follows:

[0017] ;

[0018] 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 and snow area, then , otherwise 0, represents the basin unit, represents the total number of grid units in the basin unit.

[0019] Further, the spatial resolution unification adopts a cubic convolution interpolation method, and the formula is as follows:

[0020] ;

[0021] wherein, represents an output pixel elevation value to be solved, and represent continuous coordinate values in a target coordinate system, represents an elevation value of a known input pixel, and represent discrete matrix indexes of an input grid, represents a weight kernel function of the cubic convolution interpolation, represents a normalized distance in a horizontal direction, represents a normalized distance in a vertical direction, and respectively represent weight values of the normalized distances in the horizontal direction and the vertical direction;

[0022] The weight kernel function is defined as follows:

[0023] When , ;

[0024] When , ;

[0025] When , ;

[0026] wherein, represents a normalized distance, i.e., a distance from a center point of a target pixel to a center point of a neighboring original pixel, represents an interpolation coefficient.

[0027] Further, the geographic positioning bias correction corresponds to the formula as follows:

[0028] ;

[0029] wherein, represents an elevation difference value of two DEMs, represents a terrain slope, represents a magnitude of a horizontal offset vector, represents a direction of an offset vector, represents a terrain slope direction, i.e., an orientation angle of a slope surface, represents a ratio of a vertical bias to a tangent of an average slope;

[0030] The , the formula is as follows:

[0031] ;

[0032] 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 aspect of the th sample point.

[0033] Further, the corresponding formula of the elevation distortion bias correction is as follows:

[0034] ;

[0035] wherein, denotes the pixel elevation value, denotes the elevation distortion bias estimation value, i.e. the systematic elevation error at the elevation , and denote the minimum elevation and the maximum elevation of the study area respectively, denotes the polynomial coefficient, denotes the polynomial order, denotes the polynomial order index;

[0036] The Akaike information criterion is used to determine the polynomial order, and the formula is as follows:

[0037] ;

[0038] 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 represents the error sum of squares of the model fitting.

[0039] Further, the corresponding formula of the satellite orbit model bias correction is as follows:

[0040] ;

[0041] wherein, denotes the orbit-related system bias, i.e. the elevation error estimation value at the orbit coordinate , denotes the pixel coordinate along the orbit direction, represents the coordinate along the cross-track direction, represents the trend surface coefficient, and represents the trend surface order, and represents the polynomial order index;

[0042] Orbit coordinate system coordinates obtained by a rotation transformation, the formula is as follows:

[0043] ;

[0044] wherein, and represents the coordinates of the pixel under the original map projection coordinate system, represents the orbit direction angle.

[0045] Further, the penetration depth physical model is as follows:

[0046] ;

[0047] wherein, represents the radar penetration depth, that is, the penetration depth under the condition of the elevation , snow density , and represents the snow surface elevation, represents the snow density, represents the sea level reference penetration depth, represents the radar frequency, represents the speed of light, represents the loss angle tangent, represents the snow thickness, represents the elevation gradient coefficient, represents the reference elevation, represents the constant term, represents a natural exponential function;

[0048] The calculation formula of the penetration depth difference is as follows:

[0049] ;

[0050] wherein, represents the penetration depth difference 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;

[0051] The parameter optimization formula of the penetration depth physical model is as follows:

[0052] ;

[0053] 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.

[0054] Furthermore, the formula for the robust regression is as follows:

[0055] ;

[0056] 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.

[0057] Furthermore, the formula for calculating the comprehensive index of historical activity of the high-level loose body is as follows:

[0058] ;

[0059] in, This represents a comprehensive index indicating the historical activity of high-level loose aggregates. Indicates the evaluation factor number. This represents the total number of evaluation factors. Indicates the first The weights of each evaluation factor Indicates the first Factor values ​​of each evaluation factor Indicates the first Confidence level of each assessment factor;

[0060] The formula for estimating the activity probability is as follows:

[0061] ;

[0062] in, denotes an activity probability, denotes an intercept term, denotes a regression coefficient of the evaluation factor, denotes a standardized factor value of the evaluation factor, denotes a natural exponential function.

[0063] Further, the identification of the active hot spot area and the time period adopts a spatiotemporal clustering analysis model, and the spatiotemporal clustering analysis model is as follows:

[0064] ;

[0065] wherein, denotes a spatiotemporal clustering statistic, denotes a spatial window, denotes a time window, denotes a number of high-activity points in the window, denotes an expected value based on a Poisson distribution, denotes a total number of points, denotes a window likelihood function value, denotes a total likelihood function value.

[0066] The present application has the beneficial effects that: the high loose body historical activity detection method based on multi-source DEM data provided by the present application realizes high-precision, quantitative and systematic detection of the historical activity of high loose bodies in high mountain glacier areas by constructing a complete multi-source DEM data processing flow. Specifically, through the collaborative use and systematic correction of multi-source DEM data, the observation blind area of a single data source in a harsh environment is effectively overcome, forming a continuous and reliable time series elevation data set; by establishing a progressive systematic 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

[0067] Figure 1 is a structural schematic diagram provided for an embodiment;

[0068] Figure 2 is a geometric relationship schematic diagram between DEM and offset DEM provided for an embodiment;

[0069] Figure 3 is a system bias correction effect comparison schematic diagram provided for an embodiment;

[0070] Figure 4 The typical flow basin unit elevation time sequence change curve schematic diagram provided for the embodiment is shown in the following figure:

[0071] Figure 5 The high loose body activity multi-factor comprehensive evaluation thunderstorm schematic diagram provided for the embodiment is shown in the following figure:

[0072] Figure 6 The ROC curve schematic diagram of the activity detection model provided for the embodiment is shown in the following figure:

[0073] Figure 7 The detection result confusion matrix of the activity detection model provided for the embodiment is shown in the following table:

[0074] Figure 8 The space-time clustering analysis result schematic diagram provided for the embodiment is shown in the following figure. DETAILED DESCRIPTION

[0075] In order to overcome the problems of the existing monitoring technology, such as environmental limitation, poor observation continuity, and serious signal incoherence 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 present application is proposed.

[0076] In the present application, first, the flow basin unit is divided. The movement of high loose body is controlled by terrain, such as collapse and debris flow, and the material moves along the slope and finally converges in the channel. Therefore, the flow basin 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 for 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.

[0077] 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.

[0078] Then, system bias correction is performed. There are systematic errors between different DEM data that are not caused by real terrain changes. The present application adopts a three-step progressive correction strategy: geographic positioning bias correction, height distortion bias correction, and orbit pattern bias. Specifically, at the terrain undulations, a small shift in horizontal position will appear as a huge height difference. The geographic positioning bias correction establishes a mathematical relationship between the height difference, slope, and aspect, and inversely calculates the horizontal shift amount and direction and corrects it. In high-altitude areas, ground control points are scarce, and DEM errors are often systematically related to altitude. The height distortion bias correction fits a polynomial relationship between height and error, and subtracts it as a function of altitude. The orbit and attitude of the satellite can cause a wave-like systematic error along the flight direction. The orbit pattern bias eliminates the sensor's own distortion by converting the coordinates to the orbit coordinate system and fitting a two-dimensional trend surface.

[0079] Then, radar penetration depth correction is performed. Radar signals can penetrate snow and ice layers, causing the measured height to be 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 penetration depth and snow density, water content (loss tangent), frequency, and other parameters based on the attenuation law of electromagnetic waves in snow and ice media, and uses the C-band and X-band DEM data obtained simultaneously by SRTM. Since the X-band penetration ability is much weaker than the C-band, the height difference between the two at the same location mainly reflects the difference in penetration depth. By this observable difference value, the key parameters in the physical model are inversely optimized and solved, thereby realizing spatialized penetration depth correction and modifying the radar DEM to the true snow and ice surface.

[0080] Then, multi-temporal DEM difference analysis is performed. The activity of high loose bodies must be manifested as changes in surface elevation, such as elevation decrease caused by glacier retreat and elevation increase caused by collapse accumulation. The present application performs pixel-level subtraction calculation, i.e. difference, on DEMs of different periods that have been strictly corrected, directly obtains the elevation change, and uses a robust regression method to fit the time series to obtain a reliable elevation annual change rate.

[0081] Finally, high loose body activity assessment is performed. Single elevation change is not enough to comprehensively assess activity, and it needs to be combined with multiple geomorphic activity indicators 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 to convert the comprehensive index into an occurrence probability, making the evaluation result more statistically meaningful, and identifies hotspots and time periods that are significantly aggregated in the spatial and temporal dimensions through scan statistics, revealing the group occurrence rules and key risk periods of disasters.

[0082] The technical solutions in this embodiment will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments.

[0083] Figure 1 A flowchart illustrating a method for detecting the historical activity of high-altitude loose bodies based on multi-source DEM data is shown. Please refer to [link / reference]. Figure 1 The method includes the following steps:

[0084] Step 1: Watershed Unit Division:

[0085] Based on the digital elevation model (DEM) data, watershed units naturally defined by ridgelines and valleys are automatically delineated; based on glacier boundary data, the glacier coverage of each watershed unit is calculated, and candidate watershed units significantly affected by glaciers are selected according to the set glacier coverage threshold.

[0086] In high-altitude glacial regions, the formation, transport, and deposition of high-level loose material are strictly controlled by topography. Source areas are typically located on steep periglacial slopes, transport paths are constrained by valley systems, and deposition areas are distributed at valley mouths or in areas with gentler terrain. Therefore, a watershed unit naturally defined by ridgelines (watersheds) and valley lines (catchment lines) is a closed and independent natural system that can fully reflect the complete disaster process of source-transport-deposition. It can accurately reflect the transport paths and accumulation patterns of surface materials, ensuring that the analysis results conform to the laws of geomorphological evolution and the principles of disaster dynamics. Automated watershed unit division can be achieved using mature algorithms based on confluence segmentation, such as Chinese Patent Application Publication No. CN113850822A, which will not be elaborated upon in this embodiment.

[0087] After obtaining the watershed units, it is necessary to identify candidate watershed units significantly affected by glaciers. This embodiment uses glacier coverage as a screening criterion. In terms of technical implementation, the vector format glacier boundary data (such as RGI data) is first converted into raster data with the same resolution as the DEM. For raster units that cross watershed boundaries, the area weighting method is used to accurately calculate their contribution: a weight value between 0 and 1 is assigned based on the actual coverage area of ​​the glacier within that raster, and then the weights are accumulated to obtain the glacier coverage of the entire watershed unit.

[0088] The formula for calculating glacier coverage is as follows:

[0089] ;

[0090] in, Indicates glacier coverage. and This indicates the row and column number of a cell within a watershed unit in the raster data. Represents a pixel The binary function of glacier cover, if like the element Located within the glacier area, ,otherwise 0, Represents a watershed unit. This indicates the total number of grid cells within a watershed unit.

[0091] In practical applications, a glacier coverage rate greater than 10% can be set as a screening criterion. This ensures that the selected watersheds are significantly affected by glacial activity, such as moraine supply and meltwater triggering, while also preventing the omission of areas with relatively small glacier areas that still pose potential disaster risks. For example, although some watersheds may not have a high overall glacier coverage rate, the retreat of local glaciers in their upstream or periphery may still expose large amounts of loose material, forming potential sources of disaster.

[0092] Step 2, Multi-source DEM data processing:

[0093] Acquire multi-source DEM data for the study area, including C-band and X-band DEM data of SRTM DEM, ALOSAW3D30 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 dataset.

[0094] Since a single remote sensing data source cannot simultaneously meet all requirements for temporal coverage, spatial resolution, measurement accuracy, and anti-interference capability, this embodiment employs multi-platform, multi-sensor, and multi-temporal DEM data for complementary advantages. However, these raw data differ fundamentally in imaging mechanisms, spatial references, acquisition time, and accuracy characteristics, making them unsuitable for direct high-precision differential comparison. Therefore, a rigorous preprocessing procedure is necessary to unify them to the same benchmark, forming a spatiotemporally comparable, format-consistent standard dataset to lay the foundation for subsequent precise calibration and analysis. This embodiment eliminates systematic errors caused by inconsistent data standards through a unified coordinate system, data format, and spatial resolution, ensuring direct comparability of data from different periods and sources.

[0095] Three types of DEM data are used in this implementation: First, SRTM DEM data (acquired in February 2000), which is global elevation data acquired by a radar system carried by the space shuttle. This data contains both C-band (5.7 GHz) and X-band (9.7 GHz) radar measurements with a spatial resolution of 30 meters. SRTM data is complete in coverage and stable in quality, and can be used as a reference for ground elevation in 2000. Second, ALOS AW3D30 DEM data (acquired from 2006 to 2011), generated from optical stereo pairs of the Japanese ALOS satellite, also with a spatial resolution of 30 meters. Unlike radar data, optical measurements directly obtain the surface elevation of ground objects, and there is no signal penetration problem in snow-covered areas, making 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 precision is extremely high, reaching centimeters, providing a reliable reference for the accuracy verification of other DEM data.

[0096] After obtaining multi-source DEM data, the first step is to unify the data format and coordinate system. All original data is first converted to GeoTIFF format, using 32-bit floating-point type storage to ensure the accuracy of 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 checks are performed. Verify whether the data is complete and covers the study area, whether the elevation values are within a reasonable range (such as high mountain areas are usually between 2000-8000 meters), and remove obvious outliers. Then, all data is 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, use UTM 47N projection zone, and the ellipsoid reference is unified to WGS84. This step ensures the strict consistency of different data sources in spatial position.

[0097] The second step is the unification of spatial resolution. In order to perform accurate differential analysis, all data need to be resampled to the same spatial resolution. The resampling strategy depends on the relationship between the original resolution and the target resolution: for fine data with resolution higher than 30 meters, bilinear interpolation is used for down-sampling. This method determines the new elevation value by calculating the weighted average of the four original pixels around the target pixel, which can both smooth local details and maintain the overall characteristics of the terrain. For coarse data with resolution lower than 30 meters (such as MODIS snow products with 500-meter resolution), up-sampling is needed. In this implementation, cubic convolution interpolation is used:

[0098] ;

[0099] where, is the output pixel elevation value to be solved, and are the continuous coordinate values in the target coordinate system, is the known input pixel elevation value, and are the discrete matrix indices of the input grid, is the weight kernel function of cubic convolution interpolation, is the normalized distance in the horizontal direction, is the normalized distance in the vertical direction, and are the weight values of the normalized distance in the horizontal and vertical directions, respectively;

[0100] The weight kernel function is defined as follows:

[0101] When , ;

[0102] When , ;

[0103] When , ;

[0104] where, is the normalized distance, i.e., the distance from the center of the target pixel to the center of the adjacent original pixel, is the interpolation coefficient.

[0105] The above steps effectively overcome the inherent defects of single data source in high mountain area, 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 available elevation information can be obtained in complex environment. Through strict format, coordinate and resolution triple unification, the data from different sources are integrated into a set of highly standardized multi-temporal DEM data. This eliminates the most basic systematic differences between data, so that DEM data from SRTM, ALOS and ICESat can be directly and accurately compared at the pixel level, which makes it possible to detect sub-millimeter level surface changes in the subsequent.

[0106] Step 3, system bias correction:

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

[0108] 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 systematic bias. If these biases are not corrected, they will produce false signals far exceeding the true changes in the differential analysis. The present application innovatively decomposes the complex total system error 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.

[0109] In this embodiment, the total system error can be decomposed into geographic positioning bias, elevation distortion bias related to elevation, and orbit pattern bias related to sensor orbit geometry. The three biases are different in magnitude and influence, and the geographic positioning 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 orbit-related bias is eliminated. Improper correction sequence will lead to error transmission and amplification.

[0110] Among them, the geographic positioning 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 , The magnitude of the horizontal offset vector is represented by DEM, which represents the surface line of the digital elevation model at the reference period. Indicates the slope of the terrain. This indicates the false elevation gain generated in the ridgeline region. This indicates the false elevation loss generated in the valley line area. This represents the elevation difference between the two DEM phases. This embodiment identifies and corrects this deviation by analyzing the relationship between the elevation difference and the terrain slope; the corresponding formula is as follows:

[0111] ;

[0112] in, This represents the elevation difference between the two DEM phases. Indicates the slope of the terrain. This represents the magnitude of the horizontal offset vector. Indicates the direction of the offset vector. Indicates the slope aspect, that is, the angle at which the slope faces. It represents the ratio of vertical deviation to the average slope tangent.

[0113] in, The solution is obtained using the least squares method, as shown in the following formula:

[0114] ;

[0115] 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 The 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.

[0116] Elevation distortion is the variation of DEM error with altitude. This phenomenon is particularly pronounced in mountainous regions, mainly due to the sparse distribution of ground control points at high altitudes. For example, SRTM data may exhibit a systematic bias of 5-10 meters above 4000 meters. This invention uses a polynomial model to describe this elevation dependency:

[0117] ;

[0118] in, Indicates the pixel elevation value. This represents the estimated elevation distortion deviation, i.e., at elevation... Systematic elevation error at the location, and respectively represent the minimum and maximum elevation of the study area, represent the polynomial coefficients, represent the polynomial order, represent the polynomial order index.

[0119] where the polynomial order is determined using the Akaike information criterion (AIC) as follows:

[0120] ;

[0121] where, 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, the total number of sample points used to calculate the average elevation difference in each layer of the stable area (such as bedrock) after stratification by elevation (such as 100-meter intervals), is the residual variance, representing the sum of squared errors of the model fitting.

[0122] 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 within each layer. ② Use the weighted least squares method to fit the polynomial coefficients. The weight is proportional to the number of effective sample pixels included in each layer, ensuring that the data-rich layer contributes more. ③ Use the Akaike information criterion (AIC) to adaptively determine the optimal polynomial order and select the polynomial order with the smallest AIC. ④ Use cross-validation (such as 70% training set and 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.

[0123] Orbit pattern bias is related to satellite orbit characteristics and is caused by specific geometric characteristics of the sensor (such as periodic changes in attitude, geometric distortion, and directional differences in atmospheric refraction). It manifests itself as different systematic error patterns along the satellite flight orbit direction and perpendicular to the orbit (cross-orbit) direction. For example, the SRTM DEM often exhibits long-wave characteristics with a wavelength of about 50-100 kilometers and an amplitude of 2-3 meters along the orbit direction. To correct such biases, this embodiment establishes a two-dimensional polynomial trend surface model in the satellite orbit coordinate system:

[0124] ;

[0125] where, represents the orbit-related systematic bias, i.e., the elevation error estimate at the coordinates in the orbit coordinate system, represents the coordinate of the pixel along the orbit direction, represents the coordinate along the cross-orbit direction, represents the trend surface coefficients, and denotes the order of the trend surface, and denotes the polynomial order index.

[0126] wherein the orbital coordinate system coordinates is obtained by a rotation transformation, the formula is as follows:

[0127] ;

[0128] wherein, and denotes the coordinates of the pixel under the original map projection coordinate system, denotes the orbital direction angle.

[0129] In practical applications, the orbital mode deviation correction implementation steps are as follows: ①convert the geographic coordinates of all pixels in the study area into orbital coordinates; ②fit the trend surface coefficients 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 subtract the deviation, and use the iterative method for correction. Calculate the standard deviation of the residual error after each iteration; ④stop iteration when the residual error standard deviation is less than the preset threshold (usually 0.5 meters) or reaches the maximum iteration number (usually 5 times).

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

[0131] Step 4, radar penetration depth correction:

[0132] 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.

[0133] The core principle of this step is the physical interaction of radar wave and snow-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.

[0134] This embodiment is based on the propagation characteristics of electromagnetic waves in snow and ice medium, and establishes a physical model of penetration depth:

[0135] ;

[0136] wherein, represents the radar penetration depth, i.e. the height of the snow surface measured by radar, , snow density , and the penetration depth under the condition of represents the snow surface elevation measured by radar, represents the snow density, represents the reference penetration depth of sea level, 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 a 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 a stable area, represents a natural exponential function.

[0137] The accurate calibration of the parameters of the penetration depth physical model relies on the establishment of 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.

[0138] 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 lacking field observations.

[0139] Specifically, first, the penetration depth difference is calculated:

[0140] ;

[0141] 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.

[0142] 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:

[0143] ;

[0144] 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, denotes a parameter set denotes a predicted penetration depth difference, denotes a regularization coefficient, denotes a regularization term.

[0145] 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).

[0146] 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, and perform pixel-by-pixel correction on the C-band DEM to obtain the final DEM data after penetration depth correction, whose elevation value represents the true snow and ice surface elevation.

[0147] Through the above correction process, the systematic deviation caused by the radar penetration effect can be effectively eliminated, and the elevation value truly reflects the position of the ground surface, avoiding the misjudgment of the penetration effect as ground surface activity, and fundamentally ensuring the accuracy and reliability of the final activity evaluation results. 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 foundation for accurate evaluation of glacier change and loose body activity.

[0148] Step 5, multi-temporal DEM difference analysis:

[0149] Difference calculation is performed on the final DEM data of different time periods, and the formula based on robust regression is used 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 ground surface dynamics.

[0150] After the system preprocessing, bias correction and radar penetration depth correction of multi-source DEM data, the embodiment reveals the dynamic changes of the 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 body. 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 activities such as debris flow. Therefore, accurate quantification of the spatio-temporal change pattern of ground elevation provides the most direct quantitative evidence for evaluating the historical activity of high loose body.

[0151] The traditional two-period DEM simple difference method is easily disturbed by accidental errors and outliers, affecting the reliability of 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 more robust estimates of annual average change rate by fitting the overall trend of elevation observations over time, but also effectively identify nonlinear dynamic processes, such as accelerated ablation of glaciers 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 and short-term cloud shadows), so as to obtain more reliable and more consistent with the actual physical process of elevation change rate.

[0152] In the embodiment, the formula of the robust regression is as follows:

[0153] ;

[0154] Wherein, and represent the annual change rate of elevation, represent the number of time series observations, represent the elevation observation value of the th period, represent the observation time of the th period, represent the reference elevation, represent the reference time point, represent the standard deviation of elevation observation, represent the Huber loss function.

[0155] Compared with 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, instantaneous snow, vegetation, etc., and avoid its excessive influence on trend estimation. It can provide a change trend that is more in line with the actual physical process, for example, it can identify the accelerating ablation trend of glaciers or the mutation signal of collapse events, without being misled by individual abnormal observations. This can accurately distinguish different activity intensity areas of loose bodies in the future, and provide a scientific basis for risk classification and precise prevention and control.

[0156] Step 6, high loose body activity assessment:

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

[0158] The ultimate goal of this embodiment is to quantitatively evaluate the historical activity degree of the high loose body in the high mountain glacier area by comprehensively analyzing the multi-temporal DEM difference results and auxiliary information. This evaluation not only needs to identify the area where the ground surface has a significant elevation change, but also needs to deeply understand the geomorphological meaning of these changes to judge whether they indicate potential disaster risks, thereby providing a scientific basis for risk management. Therefore, this embodiment constructs a multi-level evaluation system, from single-factor change analysis, multi-factor comprehensive evaluation to spatio-temporal pattern recognition, gradually improving the reliability and practicality of the evaluation results.

[0159] First, in order to fully capture the active characteristics of the high loose body, this embodiment establishes a comprehensive evaluation model based on the evidence weight method, which is used to calculate the historical activity comprehensive index of the high loose body, as follows:

[0160] ;

[0161] Among them, represents the historical activity comprehensive index of the high loose body, 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.

[0162] In the above formula, different evaluation factors have different contributions to activity, and their overall impact 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 characteristics of the region. Each evaluation factor is standardized to an interval to ensure that factors of different dimensions can be compared.

[0163] The core evaluation factors recommended by this embodiment and their initial weights and feature threshold suggestions are as follows: ① Elevation change rate factor , reflecting the intensity of loose body erosion / deposition, weight 0.4, threshold > 2 m / year for high activity; ② Ice lake expansion rate factor , indicating the increase of upstream glacial meltwater and the risk of breaching, weight 0.2, threshold > 20% / 10 years for high risk; ③ Glacier retreat factor , reflecting the speed of ice loss and exposure of loose material source, weight 0.2, end retreat > 200 m / 10 years for significant; ④ Channel erosion factor , indicating the intensity of channel bed incision or side erosion of activities such as debris flow, weight 0.2, width increase > 50% for active. Confidence is evaluated according to data coverage, temporal consistency and processing accuracy (e.g. high quality data = 1.0, medium = 0.7, low quality = 0.4).

[0164] Secondly, to further quantify the possibility of disaster occurrence, this embodiment constructs a logistic regression model based on historical disaster samples to predict the activity probability, i.e. considering activity as a probability event, establishing a statistical model based on historical disaster samples, converting the comprehensive index into a more intuitive occurrence probability, the formula is as follows:

[0165] ;

[0166] wherein, represents the activity probability, represents the intercept term, which is obtained by training the historical disaster samples, represents the regression coefficient of the th evaluation factor, reflecting the importance of each factor, represents the standardized factor value of the th evaluation factor, using z-score standardization, represents the natural exponential function.

[0167] Based on the activity probability, the level can be divided, for example: activity probability ≥ 0.7, determined as high activity, 0.3 ≤ activity probability < 0.7, determined as medium activity, and activity probability < 0.3, determined as low activity.

[0168] Finally, the occurrence of disasters is not random in space and time, but shows clustering. By identifying these hotspots, the group of disasters and the key risk period can be revealed. This embodiment identifies the clustering pattern (hot area and time period) of disaster activities in space and time by using the space-time scanning statistic, identifies the most likely cluster and secondary cluster with statistical significance, and the space-time cluster analysis model is as follows:

[0169] ;

[0170] 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 the Poisson distribution, represents the total number of points, represents the window likelihood function value, represents the total likelihood function value.

[0171] The technical implementation of space-time cluster 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 space-time combinations. ② Calculate the likelihood ratio of each scanning window . ③ Through Monte Carlo simulation (such as 999 random rearrangements of space-time positions), construct the empirical distribution of the likelihood ratio statistic, and calculate the empirical value E corresponding to the observed maximum likelihood ratio. ④ Identify space-time clusters with statistical significance (E < 0.05), distinguish between “most likely clusters” (with the largest likelihood ratio value) and “secondary clusters”, and generate a space-time risk hotspot distribution map.

[0172] This step constructs a comprehensive and objective quantitative evaluation system by integrating a variety of key geomorphic process indicators, significantly improving the scientificity and reliability of the evaluation results. Not only includes quantitative comprehensive index, but also easy-to-understand risk probability and clear space-time risk hotspot distribution map, this multi-dimensional output greatly facilitates the interpretation and application of the results, and provides intuitive decision-making basis for managers and decision-makers with different knowledge backgrounds.

[0173] To verify the effectiveness of the method of the present application, the southeast region of Tibet, China (94°-98°E, 27°-30°N) is selected as the test area to verify the method provided in this embodiment. The region is located at the intersection of the eastern section of the Himalayas and the Hengduan Mountains, with an altitude span of 1500-7500 meters, and a large number of oceanic glaciers are distributed. The annual precipitation in the 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:

[0174] (1) Basin unit division

[0175] The 30-meter SRTM DEM is processed by filling the concave points, a total of 3847 points, with an average depth of 0.3 meters. The D8 single flow direction algorithm is used to calculate the cumulative flow, which takes about 35 minutes. Set the cumulative flow threshold to 1000 (catchment area 0.9 km²) to extract the river network, and obtain a total length of 12,450 km of river network system. Combined with 5-meter resolution Google Earth images, manual correction is performed to correct 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%).

[0176] (2) Multi-source DEM data processing

[0177] The following core data is obtained and preprocessed: ① SRTM DEM (February 2000): C-band and X-band, resolution 30 meters, vertical accuracy ±16 meters. ② ALOS AW3D30 DEM (2006-2011): resolution 30 meters, vertical accuracy ±5 meters. ③ 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 product. The preprocessing includes: using GDAL tools to unify to GeoTIFF format (32-bit floating point); projection conversion to UTM 47N (WGS84 ellipsoid); using cubic convolution interpolation method for spatial resolution unification.

[0178] (3) System bias correction

[0179] ① Geolocation bias correction: 23,456 pixels in stable bedrock area with slope of 15°-45° were selected as control points. Iterative solution obtained SRTM horizontal offset of 12.3 meters / 47° and ALOS offset of 8.7 meters / 135°; SRTM vertical bias coefficient of -1.2 and ALOS of 0.8. After correction, the elevation difference standard deviation of stable area was reduced from 18.6 meters to 6.2 meters, with an improvement of 66.7%.

[0180] ② Elevation distortion bias correction: 60 layers were divided according to 100-meter elevation. AIC criterion determined that 3-order polynomial was adopted, and the fitting polynomial coefficients were: =-2.3, =8.7, =-12.4, =5.1. After correction, the bias mean of each elevation band was controlled within ±0.5 meters. In southeast Tibet, the elevation correlation bias between SRTM DEM and ALOS AW3D30 DEM was significant, and each 1000-meter change in altitude corresponded to about 40 meters of systematic bias.

[0181] ③ Orbit pattern bias correction: SRTM orbit angle was 10.2°, and ALOS was -8.4°. The fitting trend surface coefficients were =0.8, =-0.003, =0.001, =0.00002. After correction, the residual error standard deviation was reduced to 0.48 meters.

[0182] The system bias correction results of the embodiment are shown in Figure 3 . Figure 3 The left side is a schematic diagram of elevation bias standard deviation and its systematic bias trend. Among them, the columnar filling represents the change of elevation bias standard deviation at different processing stages, and the dashed curve represents the systematic bias trend of each 1000-meter altitude. The values (including 18.6, 6.2, 2.4, 0.48) shown in the figure are the elevation bias standard deviation values corresponding to each stage, which shows that the elevation accuracy is gradually improved with the advancement of correction steps.

[0183] Figure 3 The right side is a bias comparison column chart at different correction stages. From left to right, they are the bias of the original, after geolocation correction, after elevation distortion correction, and 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 elevation bias and gradually improve the accuracy of the measurement results.

[0184] (4) Radar penetration depth correction

[0185] Combined with field snow pit measurements (snow density 300-450 kg / m³, loss tangent: dry snow 0.0008-0.0012, wet snow 0.008-0.015, snow depth 0.5-8.0 meters) and application of the penetration depth physical model. The calculations show that at an altitude of 4000 meters, the C-band penetration is 4.3±0.7 meters, and the X-band 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 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 is 0.7±0.3 meters. The C-X band elevation optimization model parameters are used for 28450 glacier pixels. After correction, the glacier area elevation 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.

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

[0187]

[0188] (5) Multi-temporal DEM difference analysis

[0189] Select 29 ICESat data sufficient basins: as shown in Figure 4 The elevation change curves of A, B, C, D four basin units from 2000 to 2023 based on SRTM, ICESAT-1 and ICESAT-2 data. Figure 4 The specific change rate of a number of typical basin units is 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 4 The Huber regression converges after 15 iterations. Set the significance threshold (SRTM-ALOS difference standard deviation 11.4 meters x 2 = 22.8 meters; SRTM-ICESat 10.0 meters x 2 = 20.0 meters). Among the 661 basins, 412 (62.3%) were detected to have significant elevation changes: the average decrease of the glacier area was 18.7±8.2 meters (2000-2010), and the average decrease of the non-glacier area was only -0.3±3.1 meters, which verified the effectiveness of the correction. The results show that the average elevation of these areas gradually decreases over time, indicating that the glacier-covered area is continuously melting, increasing the risk of high loose body disasters.

[0190] (6) High loose body activity assessment

[0191] = 0.85 (change rate - 25.3 m / 10 years), = 0.92 (glacial lake 0.05→0.15 km2), = 0.73 (terminal retreat 385 m), = 0.67 (gully width 15→25 m), and the comprehensive index is 0.81 (high activity).

[0192] = -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%.

[0193]

[0194] The multi-factor comprehensive evaluation results of the embodiment are shown in Figure 5 . Figure 5 The main factors affecting the geological activity and their weights in the evaluation system are listed, including: topographic gradient (as a supplementary index), gully erosion (weight: 0.11), and 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 gradual change strip, which is used to correspond to different levels of geological activity. Figure 5 The effectiveness and application results of the technical scheme of the present application in comprehensively integrating multi-source indicators, quantifying weights, and generating a geological activity index are illustrated, which embodies the systematization and discrimination of the evaluation method.

[0195] ​​​(7) Result verification and technical effect

[0196] Comprehensive evaluation finally identified 78 historical active points: high activity 25 (30.3%, concentrated in the glacier terminus / ice lake periphery), medium activity 35 (46.1%, glacier retreat area), low activity 18 (23.6%, sporadic distribution). Compared with the actual disaster record from 2010 to 2020: 49 out of 60 disasters were successfully identified, with an accuracy rate 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%), three new ice lakes were formed in the main ditch (total area 0.28 km²), and the elevation showed an accelerating downward trend, indicating that the future disaster risk will continue to rise. Performance evaluation and verification effect of active detection model Figures 6 to 8

[0197] Figure 6 Fig. 7 shows the ROC curve of the active detection model, that 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, reflecting the good classification accuracy and reliability of the model in identifying geological activity.

[0198] Figure 7 Fig. 8 shows the confusion matrix of the active detection model. The confusion matrix compares the model's predicted results with the actual observed results in tabular form, 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; 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. This directly verifies that the active detection model constructed by the present application has high accuracy and high reliability.

[0199] Figure 8 Fig. 9 shows a spatiotemporal clustering analysis result diagram, which shows the clustering structure obtained by analyzing the spatiotemporal characteristics of the target region using the method of the present application, including a main cluster and multiple sub-clusters (including sub-cluster 1, sub-cluster 2, and sub-cluster 3) to represent the hierarchical features of the data in spatial distribution and temporal evolution.

[0200] ​​In summary, the high loose body historical activity detection method based on multi-source DEM data provided by the embodiment realizes accurate and quantitative detection of the historical activity of high 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 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 reduction 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 the above 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 watershed unit to correct the spatialization of the C-band DEM data, and the C-band DEM data after penetration depth correction is obtained; The processed DEM data at different time periods are subjected to difference calculation, and the height change trend of each watershed unit is fitted by a formula based on robust regression to obtain height change rate spatiotemporal distribution data reflecting the surface dynamics; Based on the height change rate spatiotemporal distribution data, the comprehensive index of historical activity of high-level 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 evaluation result of the historical activity of high-level loose body; 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 parameter optimization formula 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.

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 of the glacier coverage binary function, 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 and vertical directions. 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 body 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 periods of 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 deviation to the average slope tangent; Solving by least squares 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 aspect of the th 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 the pixel elevation value, represents the elevation twist bias estimate, i.e. the systematic elevation error at the elevation represents the elevation twist bias estimate, i.e. the systematic elevation error at the elevation and respectively represent the minimum and maximum elevation of the study area, represents the polynomial coefficient, represents the polynomial order, represents the 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 the 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 estimate of the elevation error at the orbit coordinate system coordinate denotes the trend surface coefficients, denotes the coordinates of the pixel along the orbit direction, denotes the coordinates 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 image element 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 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.

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 for calculating the comprehensive index of historical activity of high-level loose body is as follows: ； in, This represents a comprehensive index indicating the historical activity of high-level loose aggregates. Indicates the evaluation factor number. This represents the total number of evaluation factors. Indicates the first The weights of each evaluation factor Indicates the first Factor values ​​of each evaluation factor Indicates the first The confidence levels of the assessment factors, including elevation change rate factor, glacial lake expansion rate factor, glacial retreat factor, and gully erosion factor; The formula for estimating the activity probability 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.

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 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, 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