A rockfall disaster chain three-dimensional deformation monitoring method, system and medium

Abstract

This invention discloses a method, system, and medium for three-dimensional deformation monitoring of ice and rock avalanche disaster chains, belonging to the field of geological disaster remote sensing monitoring technology. The method includes the following steps: acquiring a multi-phase time-series digital elevation model (DEM) covering the target area; for each pair of adjacent time-phase DEMs: constructing a hierarchy through multi-scale decomposition, performing optical flow calculation under terrain gradient adaptive smoothing constraints to obtain an initial horizontal displacement field; using displacement amplitude to generate attention weights for cross-scale iterative fusion to obtain an optimized horizontal displacement field; registering the subsequent time-phase DEM based on this horizontal displacement field, and calculating the vertical displacement field through elevation differences; fusing horizontal and vertical displacements, and combining time intervals to generate a three-dimensional deformation field and a deformation rate field; and finally obtaining time-series three-dimensional deformation information reflecting the dynamic evolution of the disaster chain by sequentially processing all adjacent time-phases. This invention achieves high-precision, integrated time-series three-dimensional deformation monitoring of ice and rock avalanche disaster chains.

Description

Technical Field

[0001] This invention relates to the field of remote sensing monitoring technology for geological disasters, and in particular to a method, system and medium for monitoring three-dimensional deformation of ice and rock collapse disaster chains. Background Technology

[0002] High mountain and canyon areas, especially high-altitude slopes with well-developed glaciers, are prone to high-altitude geological disaster chains such as ice and rockfalls. High-precision, comprehensive three-dimensional deformation monitoring of these areas is crucial for early disaster identification and risk warning. However, the terrain in these regions is extremely complex, characterized by dramatic topographic relief, highly heterogeneous surface cover (such as a mixture of bedrock, snow, and glacial till), and frequent interference from clouds, fog, and shadows, posing significant challenges to existing monitoring technologies.

[0003] Currently, the main limitations of surface deformation monitoring methods applied to such areas are as follows: Synthetic Aperture Radar Interferometry (InSAR) technology is prone to monitoring failure or a significant decrease in accuracy in areas with extreme topographic gradients due to geometric distortion, overlay, and incoherence; the traditional Differential Digital Elevation (DOD) method can only acquire vertical deformation, completely ignoring the influence of horizontal displacement. However, in disasters such as landslides and collapses, horizontal displacement is often a more significant precursor signal, and ignoring it will lead to a serious underestimation or even misjudgment of the total deformation; the matching accuracy of pixel offset tracking technology based on optical images is highly dependent on image quality. Under the influence of variable lighting conditions, large-scale shadows, and frequent cloud and fog coverage in high mountain and canyon areas, the matching success rate and reliability are difficult to guarantee.

[0004] Therefore, existing technologies lack a robust method capable of simultaneously and accurately acquiring three-dimensional (east-west, north-south, and vertical) deformation fields of the Earth's surface under complex terrain. In particular, how to robustly calculate sub-pixel-level horizontal displacement from time-series digital elevation model (DEM) data and accurately separate the actual vertical deformation based on this has become a key technical bottleneck restricting the accurate monitoring and mechanism analysis of geological hazards in high mountain and canyon areas. The method described in claim 1 is precisely proposed to overcome this bottleneck. Summary of the Invention

[0005] The purpose of this invention is to provide a three-dimensional deformation monitoring method for ice and rock avalanche disaster chains, which realizes high-precision, integrated time-series three-dimensional deformation monitoring of ice and rock avalanche disaster chains.

[0006] To address the aforementioned technical problems, embodiments of the present invention provide a method for monitoring three-dimensional deformation of ice and rockfall disaster chains, comprising the following steps: Obtain multi-period time-series digital elevation models (DEMs) covering the target ice and rockfall disaster chain area; Multi-scale decomposition is performed on two adjacent DEMs to construct a multi-scale hierarchy from low resolution at the top to high resolution at the bottom. Starting from the top, the optical flow between corresponding levels of the two DEMs is calculated sequentially. Based on the gradient information of the DEM at that level, the smoothing constraint is weakened in the region with high terrain gradient and strengthened in the region with low terrain gradient to obtain the initial horizontal displacement field corresponding to each level. Based on the displacement amplitude of the pixels in each initial horizontal displacement field, attention weights for cross-scale fusion are generated; based on the attention weights, the initial horizontal displacement fields of adjacent layers are iteratively weighted and fused, and the optimized horizontal displacement field within the interval between the two adjacent DEM periods is output at the bottom layer. The horizontal displacement field includes east-west displacement components and north-south displacement components. Based on the optimized horizontal displacement field, the later DEM in two adjacent DEMs is spatially resampled to generate a corrected elevation model that is spatially aligned with the previous DEM; the elevation difference between the corrected elevation model and the previous DEM is determined to obtain the vertical displacement field reflecting the vertical elevation change during this period. Based on the horizontal displacement field, the vertical displacement field, and the time interval between the two DEM acquisitions, the three-dimensional deformation field and the corresponding deformation rate field within that time period are determined. Based on the three-dimensional deformation fields and corresponding deformation rate fields of all two adjacent DEM pairs, we obtain the temporal three-dimensional deformation information and deformation rate sequence that reflect the dynamic evolution of the ice and rock collapse disaster chain.

[0007] In some optional embodiments, the step of sequentially calculating the optical flow between corresponding layers of the two DEMs, and based on the gradient information of the DEM layer, reducing the smoothing constraint in regions with high terrain gradients and strengthening the smoothing constraint in regions with low terrain gradients to obtain the initial horizontal displacement field corresponding to each layer includes the following steps: The optical flow calculation is based on the Farneback dense optical flow algorithm and satisfies the assumptions of constant brightness and smooth motion; the assumption of constant brightness is approximated by the following first-order Taylor expansion: ; In the formula, The gradient of the DEM; This represents the east-west displacement component; dy This represents the displacement component in the north-south direction; The motion smoothing assumption is achieved by minimizing the following energy function: ; In the formula, E It is an energy function; For horizontal displacement field, , ; is the smoothing coefficient, where the smoothing constraint is weakened in areas with high terrain gradients and strengthened in areas with low terrain gradients. For the image domain.

[0008] In some optional embodiments, the step of generating attention weights for cross-scale fusion based on the displacement magnitudes of pixels in each initial horizontal displacement field, and iteratively weighting and fusing the initial horizontal displacement fields of adjacent layers based on the attention weights to output the optimized horizontal displacement field within the interval between two adjacent DEM periods at the bottom layer, includes the following steps: Calculate the displacement amplitude M of the initial horizontal displacement field of the current layer, where the displacement amplitude of each pixel is calculated using the following formula: ; In the formula, dy represents the east-west displacement component; dy represents the north-south displacement component. The displacement amplitude map M is Gaussian smoothed to obtain smoothed amplitude. The smoothed amplitude Normalization to The normalized weight formula is obtained from the interval as follows: ; In the formula, For normalized weights; To smooth the amplitude The maximum value in; For the normalized weights By applying a minimum weight constraint, we obtain the attention weights. W The formula is as follows: ; In the formula, This is the preset minimum weight value; Using the attention weight W, the optical flow field of the current level is... fused optical flow field with the previous adjacent level Weighted fusion is performed to obtain the fused optical flow field at the current level. The formula is as follows: ; In the formula, represents the fused optical flow field at the current level; W represents the attention weight. The optical flow field at the current level; This is the fused optical flow field of the adjacent layer above.

[0009] In some optional embodiments, the step of spatially resampling the later DEM between two adjacent DEMs based on the optimized horizontal displacement field, and calculating the elevation at the sub-pixel position after resampling using bilinear interpolation, includes the following steps: For mapped coordinates of sub-pixel positions The point, where x and y are the integer parts. Subpixel offset, interpolated elevation Calculate the elevation values ​​based on the elevation values ​​of the surrounding four integer grid points using the following formula: ; In the formula, For the second phase DEM in integer coordinates The elevation value at that location.

[0010] In some optional embodiments, determining the three-dimensional deformation field and the corresponding deformation rate field within the time period specifically includes: The horizontal resultant displacement field is calculated based on the horizontal displacement field, and the value of each pixel in the horizontal resultant displacement field is: ; In the formula, The horizontal displacement for each pixel; dy represents the east-west displacement component; dy represents the north-south displacement component. Based on the horizontal displacement field, the vertical displacement field, and the acquisition time interval Calculate the three-dimensional deformation rate field, which includes the east-west deformation rate. North-south deformation rate Vertical deformation rate ; The horizontal resultant velocity field is calculated based on the three-dimensional deformation rate field, and the value of each pixel in the horizontal resultant velocity field is: ; In the formula, The horizontal velocity of each pixel; The east-west deformation rate; The north-south deformation rate is represented by the value 't'.

[0011] In some optional embodiments, the method further includes preprocessing and registration of two adjacent DEMs, including the following steps: For invalid value regions in the digital elevation model, inverse distance weighted interpolation is used for filling; Gaussian filtering is applied to the filled digital elevation model for denoising; based on the denoised two-phase digital elevation models, the horizontal offset vector and overall elevation deviation between the two phases of digital elevation models are solved using a deviation correction model under terrain constraints; one phase of the digital elevation model is translated and its elevation is corrected according to the horizontal offset vector and overall elevation deviation to achieve registration of the two phases of digital elevation models.

[0012] In some optional embodiments, the method further includes optimizing the vertical displacement field, comprising the following steps: based on In principle, outliers in the vertical displacement field that exceed the range of "mean ± 3 times standard deviation" are removed; morphological closing operation filtering is performed on the vertical displacement field after removing outliers, and the closing operation is first dilation and then erosion.

[0013] Embodiments of the present invention also provide a three-dimensional deformation monitoring system for ice and rockfall disaster chains, comprising: Data acquisition module: used to acquire multi-period time-series digital elevation models (DEMs) covering the target ice and rockfall disaster chain area; The horizontal displacement field calculation module is used to perform multi-scale decomposition on two adjacent DEMs to construct a multi-scale hierarchy from a low-resolution top layer to a high-resolution bottom layer. Starting from the top layer, it sequentially calculates the optical flow between corresponding layers of the two DEMs. Based on the gradient information of the DEM layer, it weakens the smoothing constraint in areas with high terrain gradients and strengthens the smoothing constraint in areas with low terrain gradients to obtain the initial horizontal displacement field corresponding to each layer. Based on the displacement amplitude of the pixels in each initial horizontal displacement field, it generates attention weights for cross-scale fusion. Based on the attention weights, iteratively weighted and fused the initial horizontal displacement fields of adjacent layers, and outputs the optimized horizontal displacement field within the interval between the two adjacent DEMs at the bottom layer. The horizontal displacement field includes east-west displacement components and north-south displacement components. Vertical displacement field calculation module: used to spatially resample the later DEM in two adjacent DEMs based on the optimized horizontal displacement field, generate a corrected elevation model that is spatially aligned with the previous DEM; determine the elevation difference between the corrected elevation model and the previous DEM, and obtain the vertical displacement field reflecting the vertical elevation change during this period. Three-dimensional deformation synthesis module: used to determine the three-dimensional deformation field and the corresponding deformation rate field within the time period based on the horizontal displacement field, vertical displacement field and the acquisition time interval of the two DEMs; The time series analysis module is used to obtain time series three-dimensional deformation information and deformation rate sequence reflecting the dynamic evolution of the ice and rock collapse disaster chain based on the three-dimensional deformation field and the corresponding deformation rate field of all two adjacent DEM pairs.

[0014] Embodiments of the present invention also provide a computer device, including: at least one processor; and a memory communicatively connected to the at least one processor; wherein the memory stores instructions executable by the at least one processor, the instructions being executed by the at least one processor to enable the at least one processor to perform the above-described three-dimensional deformation monitoring method for ice and rockfall disaster chains.

[0015] Embodiments of the present invention also provide a computer-readable storage medium storing a computer program, which, when run by a processor, is capable of executing the above-described method for monitoring the three-dimensional deformation of ice and rock avalanche disaster chains.

[0016] The three-dimensional deformation monitoring method for ice and rockfall disaster chains provided by this invention has at least the following beneficial effects: This invention effectively resolves the contradiction between "missed large displacements" and "loss of local details" in single-scale methods for complex terrain areas by constructing a multi-scale Gaussian pyramid optical flow calculation framework, achieving high-precision synchronous monitoring of horizontal and vertical three-dimensional deformation of the land surface. It innovatively introduces terrain gradient as a priori constraint, suppressing noise-induced spurious deformation in flat areas and enhancing motion sensitivity in steep areas, significantly improving the robustness and anti-interference ability of displacement calculation. Furthermore, a cross-scale attention fusion mechanism based on optical flow amplitude is designed to intelligently weight and fuse multi-scale results, ensuring the overall consistency and detail fidelity of the displacement field. On this basis, geometric correction eliminates the interference of horizontal displacement on vertical deformation calculation, overcoming the error problem caused by neglecting horizontal displacement in traditional difference methods. Finally, a complete quantitative product containing three-dimensional displacement components, deformation rate field, and resultant indices is output, providing direct and reliable decision support for mechanism analysis, risk zoning, and early warning of high-altitude geological hazards, forming a complete high-precision three-dimensional deformation tracking solution for complex terrain. Attached Figure Description

[0017] The accompanying drawings, which are included to provide a further understanding of the invention and form part of this invention, illustrate exemplary embodiments of the invention and are used to explain the invention, but do not constitute an undue limitation of the invention. In the drawings:

[0018] Figure 1 This is a flowchart of a three-dimensional deformation monitoring method for ice and rockfall disaster chains according to an embodiment of the present invention; Figure 2 This is a schematic diagram of the research area provided according to an embodiment of the present invention; Figure 3 This is a comparison diagram of the results before and after fusing multi-scale attention according to an embodiment of the present invention. Figure 3 (a) and Figure 3 (b) Histogram of the distribution of horizontal displacement dimension Figure 3 (c) Histogram of vertical displacement dimension; Figure 4 This is a schematic diagram of three-dimensional deformation and deformation rate provided according to an embodiment of the present invention. Detailed Implementation

[0019] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention will be clearly and completely described below in conjunction with specific embodiments and corresponding drawings. Obviously, the described embodiments are only a part of the embodiments of this invention, and not all of them. Based on the embodiments of this invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this invention.

[0020] One embodiment of the present invention relates to a three-dimensional deformation monitoring method for ice and rock avalanche disaster chains. The implementation details of the three-dimensional deformation monitoring method for ice and rock avalanche disaster chains in this embodiment are described in detail below. The following content is only for the convenience of understanding and is not necessary for implementing this solution.

[0021] The specific process of the three-dimensional deformation monitoring method for ice avalanche disaster chains in this embodiment can be described as follows: Figure 1 As shown, it includes: Step 101: Obtain a multi-period time-series digital elevation model (DEM) covering the target ice and rockfall disaster chain area; DEM data can be acquired through satellite remote sensing, aerial photogrammetry, or ground surveying.

[0022] For invalid value regions in the digital elevation model, inverse distance weighted interpolation is used for filling; Gaussian filtering is applied to the filled digital elevation model for denoising; based on the denoised two-phase digital elevation models, the horizontal offset vector and overall elevation deviation between the two phases of digital elevation models are solved using a deviation correction model under terrain constraints; one phase of the digital elevation model is translated and its elevation is corrected according to the horizontal offset vector and overall elevation deviation to achieve registration of the two phases of digital elevation models.

[0023] For invalid values ​​in the DEM (such as shaded areas and cloud cover), inverse distance weighted interpolation (IDW) is used—the closer to the valid pixel, the greater the weight, ensuring that the filled values ​​conform to terrain continuity. For invalid pixels... Its fill value The calculation formula is:

[0024] ; In the formula, For the first [unit] within the search radius Elevation values ​​of one effective pixel; for With the Euclidean distance of effective pixels ; p This is the distance attenuation coefficient; n This represents the number of valid pixels within the search radius.

[0025] The DEM is smoothed by Gaussian kernel convolution, which suppresses high-frequency noise (such as random sensor errors) while preserving key terrain features such as glaciers.

[0026] Two-dimensional Gaussian kernel function The calculation formula is: ; DEM smoothed values Convolution of the Gaussian kernel with the original DEM: ; In the formula, The standard deviation is Gaussian (0.5~1.5 in the code). A value of 1.5 provides a stronger smoothing effect, making it suitable for DEMs with high noise levels. K For the size of the nucleus (usually taken as...) (This ensures that 99.7% of the weights are contained within the kernel).

[0027] When two DEMs are horizontally misaligned, the elevation difference between them ( ) and the slope of the terrain ( ), slope aspect ( The characteristics (magnitude, direction) of the offset vector and the offset vector have a clear mathematical relationship: By combining slope, aspect, and offset parameters, the calculation formula for the multi-scale deviation correction model under stable terrain constraints is obtained: In the formula, This represents the overall elevation deviation between the two DEMs. It is the magnitude of the horizontal offset vector. Related to the offset direction, The slope of the terrain (with due north as the azimuth).

[0028] Based on the solved offset parameters (horizontal offset vector, overall elevation deviation), the DEM is translated and its elevation is corrected to achieve registration.

[0029] Step 102: Perform multi-scale decomposition on two adjacent DEMs to construct a multi-scale hierarchy from low resolution at the top to high resolution at the bottom. Starting from the top, calculate the optical flow between the corresponding layers of the two DEMs in sequence. Based on the gradient information of the DEM at that layer, reduce the smoothing constraint in the region with high terrain gradient and strengthen the smoothing constraint in the region with low terrain gradient to obtain the initial horizontal displacement field corresponding to each layer. The optical flow calculation is based on the Farneback dense optical flow algorithm and satisfies the assumptions of constant brightness and smooth motion; the assumption of constant brightness is approximated by the following first-order Taylor expansion: ; In the formula, The gradient of the DEM; This represents the east-west displacement component; dy This represents the displacement component in the north-south direction; The motion smoothing assumption is achieved by minimizing the following energy function: ; In the formula, E It is an energy function; For horizontal displacement field, , ; is the smoothing coefficient, where the smoothing constraint is weakened in areas with high terrain gradients and strengthened in areas with low terrain gradients. For the image domain.

[0030] Step 103: Based on the displacement amplitude of the pixels in each initial horizontal displacement field, generate attention weights for cross-scale fusion; based on the attention weights, iteratively weighted fuse the initial horizontal displacement fields of adjacent layers, and output the optimized horizontal displacement field within the interval between the two adjacent DEMs at the bottom layer. The horizontal displacement field includes east-west displacement components and north-south displacement components. Calculate the displacement amplitude M of the initial horizontal displacement field of the current layer, where the displacement amplitude of each pixel is calculated using the following formula: ; In the formula, dy represents the east-west displacement component; dy represents the north-south displacement component. The displacement amplitude map M is Gaussian smoothed to obtain smoothed amplitude. The smoothed amplitude Normalization to The normalized weight formula is obtained from the interval as follows: ; In the formula, For normalized weights; To smooth the amplitude The maximum value in; For the normalized weights By applying a minimum weight constraint, we obtain the attention weights. W The formula is as follows: ; In the formula, This is the preset minimum weight value; Using the attention weight W, the optical flow field of the current level is... fused optical flow field with the previous adjacent level Weighted fusion is performed to obtain the fused optical flow field at the current level. The formula is as follows: ; In the formula, represents the fused optical flow field at the current level; W represents the attention weight. The optical flow field at the current level; This is the fused optical flow field of the adjacent layer above.

[0031] Step 104: Based on the optimized horizontal displacement field, spatially resample the later DEM in two adjacent DEM periods to generate a corrected elevation model spatially aligned with the previous DEM; determine the elevation difference between the corrected elevation model and the previous DEM to obtain the vertical displacement field reflecting the vertical elevation change during this period. Based on the optimized horizontal displacement field, the second DEM in the two-phase DEM is spatially resampled, and the elevation at the sub-pixel position after resampling is calculated using bilinear interpolation, including the following steps: For mapped coordinates of sub-pixel positions The point, where x and y are the integer parts. Subpixel offset, interpolated elevation Calculate the elevation values ​​based on the elevation values ​​of the surrounding four integer grid points using the following formula: ; In the formula, For the second phase DEM in integer coordinates The elevation value at that location.

[0032] For the invalid value regions in the digital elevation model, the inverse distance weighted interpolation method is used to fill them; Gaussian filtering is applied to the filled digital elevation model for noise reduction; based on the denoised two-phase digital elevation models, the deviation correction model under terrain constraints is used to solve the horizontal offset vector and overall elevation deviation between the two-phase digital elevation models. Based on the horizontal offset vector and the overall elevation deviation, the digital elevation model of one phase is translated and its elevation is corrected to achieve registration of the two phases of digital elevation models.

[0033] Step 105: Based on the horizontal displacement field, vertical displacement field and the acquisition time interval of the two DEMs, determine the three-dimensional deformation field and the corresponding deformation rate field within the time period. The horizontal resultant displacement field is calculated based on the horizontal displacement field, and the value of each pixel in the horizontal resultant displacement field is: ; In the formula, The horizontal displacement for each pixel; dy represents the east-west displacement component; dy represents the north-south displacement component. Based on the horizontal displacement field, the vertical displacement field, and the acquisition time interval Calculate the three-dimensional deformation rate field, which includes the east-west deformation rate. North-south deformation rate Vertical deformation rate ; The horizontal resultant velocity field is calculated based on the three-dimensional deformation rate field, and the value of each pixel in the horizontal resultant velocity field is: ; In the formula, The horizontal velocity of each pixel; The east-west deformation rate; The north-south deformation rate is represented by the value 't'.

[0034] Step 106: Based on the three-dimensional deformation fields and corresponding deformation rate fields of all two adjacent DEM pairs, obtain the time-series three-dimensional deformation information and deformation rate sequence that reflect the dynamic evolution of the ice and rock collapse disaster chain.

[0035] Embodiments of the present invention also provide a three-dimensional deformation monitoring system for ice and rock avalanche disaster chains. The implementation details of the three-dimensional deformation monitoring system for ice and rock avalanche disaster chains in this embodiment are described in detail below. The following content is only for the convenience of understanding the implementation details and is not necessary for implementing this solution.

[0036] Specifically, the data acquisition module is used to acquire multi-period time-series digital elevation models (DEMs) covering the target ice and rock avalanche disaster chain area; and to process each pair of adjacent DEMs in the multi-period time-series DEMs sequentially. The horizontal displacement field calculation module is used to perform multi-scale decomposition on two adjacent DEMs to construct a multi-scale hierarchy from a low-resolution top layer to a high-resolution bottom layer. Starting from the top layer, it sequentially calculates the optical flow between corresponding layers of the two DEMs. Based on the gradient information of the DEM layer, it weakens the smoothing constraint in areas with high terrain gradients and strengthens the smoothing constraint in areas with low terrain gradients to obtain the initial horizontal displacement field corresponding to each layer. Based on the displacement amplitude of the pixels in each initial horizontal displacement field, it generates attention weights for cross-scale fusion. Based on the attention weights, iteratively weighted and fused the initial horizontal displacement fields of adjacent layers, and outputs the optimized horizontal displacement field within the interval between the two adjacent DEMs at the bottom layer. The horizontal displacement field includes east-west displacement components and north-south displacement components. Vertical displacement field calculation module: used to spatially resample the later DEM in two adjacent DEMs based on the optimized horizontal displacement field, generate a corrected elevation model that is spatially aligned with the previous DEM; determine the elevation difference between the corrected elevation model and the previous DEM, and obtain the vertical displacement field reflecting the vertical elevation change during this period. Three-dimensional deformation synthesis module: used to determine the three-dimensional deformation field and the corresponding deformation rate field within the time period based on the horizontal displacement field, vertical displacement field and the acquisition time interval of the two DEMs; The time series analysis module is used to obtain time series three-dimensional deformation information and deformation rate sequence that reflect the dynamic evolution of the ice and rock collapse disaster chain by sequentially processing all adjacent two-period DEM pairs.

[0037] Example 1: This embodiment selects the Birch Glacier in the Lötschental Valley of Valais, Switzerland, where an icefall occurred on May 28, 2025, as the study area. Temporal DEM data was acquired using stereo photogrammetry, including two periods: May 10, 2025 (before the disaster) and May 29, 2025 (after the disaster), with a ground resolution of 10 meters. In the preprocessing stage, using the May 10 DEM as a baseline, a multi-scale deviation correction model under stable terrain constraints was used to register the May 29 DEM. Subsequently, Gaussian filtering (standard deviation 1.5) was applied to remove noise, and invalid value areas were filled using inverse distance weighted interpolation (search radius 5 pixels).

[0038] The study area of ​​this invention embodiment is as follows Figure 2As shown, against a backdrop of mountainous terrain encompassing both vegetated and snow-covered areas, the distribution of ice and rockfalls in the region is clearly displayed: the red areas correspond to ice avalanche occurrence points, the pink areas are rockfall areas, and the large gray area outlined in blue represents the disaster impact range of this ice and rockfall. The compass in the upper right corner also clearly indicates the region's location. Overall, the data intuitively presents the location of ice and rockfalls within the study area and the actual impact range of the disaster's spread, making it a practical application of the aforementioned monitoring methods.

[0039] In the 3D displacement field calculation, a 5-layer Gaussian pyramid was first constructed from two DEMs, with a top layer resolution of 64 meters and a bottom layer resolution of 10 meters. Starting from the top layer, the Lucas-Kanade optical flow algorithm was used to calculate the initial displacement layer by layer, and the result of the previous layer was used as the initial value for the next layer for iterative optimization. Based on the slope information calculated from the DEM, gradient constraint optimization was implemented: in the gentle region with a slope of less than 15°, the smoothness constraint was enhanced to suppress noise, and in the steep region with a slope of greater than 45°, the constraint was weakened to enhance motion sensitivity. Subsequently, attention weights were generated based on the optical flow amplitude (annualized rate) of each layer for cross-scale fusion: the high confidence region with an amplitude greater than 2 m / year had a weight of 0.9, the medium confidence region with an amplitude of 0.5-2 m / year had a weight of 0.6, and the low confidence region with an amplitude less than 0.5 m / year had a weight of 0.2. The results of the five layers were fused by weighted averaging, and the optimized horizontal displacement field (dx, dy) was finally output.

[0040] The comparison diagram of the results before and after fusing multi-scale attention in the embodiments of the present invention is shown below. Figure 3 As shown, the histograms include distribution histograms for the two horizontal displacement dimensions (a) and (b), and distribution histogram for the vertical displacement dimension (c). The horizontal axis corresponds to the displacement amount, and the vertical axis represents the displacement frequency. The figure uses different colored curves to distinguish between the results before fusion (red) and after fusion (blue). The mean, standard deviation, and normalized mean absolute error (NMAE) of both are also labeled. This intuitively demonstrates that after fusing multi-scale attention, the distribution of displacement results is more concentrated and the statistical errors (standard deviation, NMAE) are lower, clearly showing the improvement effect of this method on the accuracy and stability of displacement monitoring.

[0041] To obtain the vertical displacement, the second-phase DEM was first mapped and geometrically corrected using the horizontal displacement field (dx, dy). A corrected DEM, spatially aligned with the first-phase DEM, was then generated using bilinear interpolation. After calculating the elevation difference between the two DEMs to obtain the initial vertical displacement field, outliers were removed using the 3σ principle, and morphological closing operations (disk structural element radius 3) were applied for filtering. The final optimized vertical displacement field (dz) was then output. Experimental results show that, after processing using this method, the mean deviation and standard deviation of vertical deformation in the study area significantly decreased from 1.67±15.11 meters to 1.62±14.74 meters. The standard deviation of deformation estimation in all study areas improved by 2% to 53%, effectively verifying the reliability and superiority of this method for high-precision three-dimensional deformation monitoring in complex terrain areas.

[0042] The three-dimensional deformation and deformation rate of the embodiments of the present invention are as follows: Figure 4 As shown in the figure, the topographic displacement and vertical displacement distribution in the X direction (east-west) and Y direction (north-south) of the study area are displayed, and the annual average rate distribution in these three directions is also presented. Each sub-figure is presented as a topographic schematic diagram with color scales, and the color change corresponds to the magnitude of the displacement or rate. Combined with the ice avalanche areas marked in red, it intuitively reflects the degree and rate of deformation of the ice avalanche-related areas in different directions during the period, and clearly demonstrates the quantitative monitoring effect of the aforementioned monitoring method on three-dimensional deformation and deformation rate.

[0043] The steps of the various methods described above are only for clarity. In practice, they can be combined into one step or some steps can be split into multiple steps. As long as they include the same logical relationship, they are all within the protection scope of this invention. Adding insignificant modifications or introducing insignificant designs to the algorithm or process, without changing the core design of the algorithm and process, are also within the protection scope of this invention.

[0044] Another embodiment of the present invention relates to a computer-readable storage medium storing a computer program. When executed by a processor, the computer program implements the method embodiments described above.

[0045] That is, those skilled in the art will understand that all or part of the steps in the methods of the above embodiments can be implemented by a program instructing related hardware. This program is stored in a storage medium and includes several instructions to cause a device (which may be a microcontroller, chip, etc.) or processor to execute all or part of the steps of the methods described in the various embodiments of the present invention. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0046] Those skilled in the art will understand that the above embodiments are specific embodiments for implementing the present invention, and in practical applications, various changes can be made to them in form and detail without departing from the spirit and scope of the present invention.

Claims

1. A method for monitoring a rock avalanche disaster chain three-dimensional deformation, characterized in that, The method comprises: acquiring multi-period time-series digital elevation models DEMs covering a target ice-rock avalanche disaster chain area; performing multi-scale decomposition on two adjacent DEMs to construct a multi-scale hierarchy from a top layer with low resolution to a bottom layer with high resolution; starting from the top layer, calculating optical flow between corresponding levels of the two DEMs in sequence, obtaining initial horizontal displacement fields corresponding to each level based on gradient information of the DEM at the level, weakening smoothing constraints in areas with high terrain gradient and enhancing smoothing constraints in areas with low terrain gradient, and obtaining the initial horizontal displacement fields corresponding to each level; generating attention weights for cross-scale fusion based on displacement amplitudes of pixel points in the initial horizontal displacement fields; and iteratively weighting and fusing the initial horizontal displacement fields of adjacent levels based on the attention weights to output an optimized horizontal displacement field within an interval of the two adjacent DEMs at the bottom layer, wherein the horizontal displacement field comprises an east-west displacement component and a south-north displacement component; spatially resampling a later DEM in the two adjacent DEMs based on the optimized horizontal displacement field to generate a corrected elevation model spatially aligned with an earlier DEM; and determining an elevation difference between the corrected elevation model and the earlier DEM to obtain a vertical displacement field reflecting vertical elevation changes within the interval; determining a three-dimensional deformation field and a corresponding deformation rate field within the interval based on the horizontal displacement field, the vertical displacement field, and a time interval between the two DEMs; obtaining time-series three-dimensional deformation information and a deformation rate sequence reflecting a dynamic evolution process of the ice-rock avalanche disaster chain according to the three-dimensional deformation field and the corresponding deformation rate field of all pairs of adjacent two-period DEMs.

2. The ice-rock avalanche disaster chain three-dimensional deformation monitoring method according to claim 1, characterized in that, The step of sequentially calculating optical flow between corresponding levels of the two DEMs, obtaining initial horizontal displacement fields corresponding to each level based on gradient information of the DEM at the level, weakening smoothing constraints in areas with high terrain gradient, and enhancing smoothing constraints in areas with low terrain gradient, comprises the following steps: The optical flow calculation is based on a Farneback dense optical flow algorithm and meets a brightness constancy assumption and a motion smoothness assumption; the brightness constancy assumption is approximately expressed by the following first-order Taylor expansion: ； wherein is the gradient of the DEM; is the east-west displacement component; dy is the north-south displacement component; The motion smoothness assumption is realized by minimizing the following energy function: ； wherein, E is an energy function; is a horizontal displacement field, , ; is a smoothing coefficient, where the smoothing constraint is weakened in areas of high terrain gradient and strengthened in areas of low terrain gradient; is an image domain.

3. The ice-rock avalanche disaster chain three-dimensional deformation monitoring method according to claim 1, characterized in that, The step of generating attention weights for cross-scale fusion based on displacement amplitudes of pixel points in the initial horizontal displacement fields, and iteratively weighting and fusing the initial horizontal displacement fields of adjacent levels based on the attention weights to output an optimized horizontal displacement field within an interval of the two adjacent DEMs at the bottom layer, comprises the following steps: calculating displacement amplitudes M of the initial horizontal displacement field of the current level, wherein the displacement amplitude formula of each pixel is as follows: ； In the formula, dx is the east-west displacement component; and dy is the north-south displacement component. Gaussian smoothing is performed on the displacement amplitude map M to obtain a smoothed amplitude ; the smoothed amplitude is normalized to an interval, and a normalized weight formula is as follows: ； wherein is a normalized weight; is a smoothing amplitude is a maximum value in to the normalized weights applying a minimum weight constraint to obtain attention weights W , as follows: ； In the formula, is a preset minimum weight value; Using the attention weights W, the optical flow field of the current level is fused with the optical flow field of the previous adjacent level to obtain the fused optical flow field of the current level by weighted fusion, as follows: ​ ； In the formula, is the fused optical flow field of the current level; W is the attention weight; is the optical flow field of the current level; is the fused optical flow field of the previous adjacent level.

4. The ice-rock avalanche disaster chain three-dimensional deformation monitoring method according to claim 1, characterized in that, The step of spatially resampling a later DEM in the two adjacent DEMs based on the optimized horizontal displacement field, and calculating the elevation of a sub-pixel position after resampling by bilinear interpolation, comprises the following steps: For a point whose mapped coordinates are sub-pixel locations where x, y are integer parts, is the sub-pixel offset, the interpolated elevation is calculated according to the elevation values of the surrounding four integer grid points by the following formula: ； wherein is the elevation value of the second iteration DEM at integer coordinates x, y.

5. The ice-rock avalanche disaster chain three-dimensional deformation monitoring method according to claim 1, characterized in that, The step of determining a three-dimensional deformation field and a corresponding deformation rate field within the interval specifically comprises: calculating a horizontal combined displacement field according to the horizontal displacement field, wherein the value of each pixel in the horizontal combined displacement field is: ； wherein is the horizontal displacement of each pixel; is the east-west displacement component; dy is the north-south displacement component; from the horizontal displacement field, the vertical displacement field and the acquisition time interval a three-dimensional deformation rate field comprising an east-west deformation rate a north-south deformation rate a vertical deformation rate ; calculating a horizontal combined velocity field according to the three-dimensional deformation rate field, wherein the value of each pixel in the horizontal combined velocity field is: ； wherein is the horizontal resultant velocity for each pixel; is the east-west deformation rate; is the north-south deformation rate.

6. The ice-rock avalanche disaster chain three-dimensional deformation monitoring method according to claim 1, characterized in that, It also includes preprocessing and registration of adjacent DEMs, including the following steps: For invalid value regions in the digital elevation model, inverse distance weighted interpolation is used for filling; Gaussian filtering is applied to the filled digital elevation model for denoising; based on the denoised two-phase digital elevation models, the horizontal offset vector and overall elevation deviation between the two phases of digital elevation models are solved using a deviation correction model under terrain constraints; one phase of the digital elevation model is translated and its elevation is corrected according to the horizontal offset vector and overall elevation deviation to achieve registration of the two phases of digital elevation models.

7. The ice-rock avalanche disaster chain three-dimensional deformation monitoring method according to claim 1, characterized in that, It also includes optimizing the vertical displacement field, including the following steps: based on In principle, outliers in the vertical displacement field that exceed the range of "mean ± 3 times standard deviation" are removed; morphological closing operation filtering is performed on the vertical displacement field after removing outliers, and the closing operation is first dilation and then erosion.

8. A rockfall hazard chain three-dimensional deformation monitoring system, characterized in that, The system includes: Data acquisition module: used to acquire multi-period time-series digital elevation models (DEMs) covering the target ice and rockfall disaster chain area; The horizontal displacement field calculation module is used to perform multi-scale decomposition on two adjacent DEMs to construct a multi-scale hierarchy from a low-resolution top layer to a high-resolution bottom layer. Starting from the top layer, it sequentially calculates the optical flow between corresponding layers of the two DEMs. Based on the gradient information of the DEM layer, it weakens the smoothing constraint in areas with high terrain gradients and strengthens the smoothing constraint in areas with low terrain gradients to obtain the initial horizontal displacement field corresponding to each layer. Based on the displacement amplitude of the pixels in each initial horizontal displacement field, it generates attention weights for cross-scale fusion. Based on the attention weights, iteratively weighted and fused the initial horizontal displacement fields of adjacent layers, and outputs the optimized horizontal displacement field within the interval between the two adjacent DEMs at the bottom layer. The horizontal displacement field includes east-west displacement components and north-south displacement components. Vertical displacement field calculation module: used to spatially resample the later DEM in two adjacent DEMs based on the optimized horizontal displacement field, generate a corrected elevation model that is spatially aligned with the previous DEM; determine the elevation difference between the corrected elevation model and the previous DEM, and obtain the vertical displacement field reflecting the vertical elevation change during this period. Three-dimensional deformation synthesis module: used to determine the three-dimensional deformation field and the corresponding deformation rate field within the time period based on the horizontal displacement field, vertical displacement field and the acquisition time interval of the two DEMs; The time series analysis module is used to obtain time series three-dimensional deformation information and deformation rate sequence reflecting the dynamic evolution of the ice and rock collapse disaster chain based on the three-dimensional deformation field and the corresponding deformation rate field of all two adjacent DEM pairs.

9. A computer system, characterized by include: At least one processor; And a memory communicatively connected to the at least one processor; wherein the memory stores instructions executable by the at least one processor, the instructions being executed by the at least one processor to enable the at least one processor to perform the three-dimensional deformation monitoring method for ice and rockfall disaster chains as described in any one of claims 1 to 7.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a computer program, which, when executed by a processor, is capable of performing the three-dimensional deformation monitoring method for ice and rock avalanche disaster chains as defined in any one of claims 1 to 7.

Images