A high-cold canyon GB-SAR high-precision deformation information extraction method and system

Abstract

The application provides a high-cold canyon GB-SAR high-precision deformation information extraction method and system, and belongs to the technical field of remote sensing detection. A GB-SAR time series interference data set of a research area is acquired, and the wrapped phase is corrected; the corrected phase is subjected to partition compensation in combination with a DEM elevation layering strategy; a high-density effective feature point set is constructed by screening PS candidate points and SBAS candidate points by using a polarized scattering power enhanced coherence coefficient field; a triangular grid is constructed based on the feature points, time series unwrapped phase is obtained by developing triangular grid phase unwrapping, and finally, a linear deformation model is constructed to complete inversion, deformation rate and cumulative deformation information of each feature point are output, and surface deformation extraction is realized. The application can solve the problem of poor applicability of the traditional uniform atmospheric assumption model, improve the scattering point recognition ability in the vegetation area, and improve the monitoring precision and reliability of the high-cold canyon deformation monitoring.

Description

Technical Field

[0001] This invention belongs to the field of remote sensing technology, and in particular relates to a method and system for high-precision deformation information extraction from high-altitude cold canyons using GB-SAR. Background Technology

[0002] The high-altitude canyon region is characterized by dramatic topographic relief, large elevation differences, and complex micro-topography. Local atmospheric temperature, humidity, and pressure parameters exhibit significant stratified variations with elevation, and atmospheric turbulence effects are prominent. Traditional deformation monitoring techniques based on GB-SAR (Ground-Based Synthetic Aperture Radar) mostly employ a single uniform atmospheric model for phase correction, failing to consider the unique elevation-stratified atmospheric distortion characteristics of high-altitude canyons. This makes it difficult to accurately characterize elevation-related nonlinear atmospheric phase errors, resulting in high residual noise and limited phase correction accuracy, which severely affects the accuracy of subsequent deformation inversion. Meanwhile, in the high-altitude canyon region, vegetation and bare rock areas are interspersed, resulting in significant differences in the temporal stability of surface scattering. Vegetated areas exhibit rapid coherence decay and weak effective scattering signals. Traditional deformation monitoring methods, which rely on a single coherence coefficient threshold to select feature points, lack sufficient ability to identify effective scattering targets in weakly coherent vegetation areas. This leads to problems such as sparse effective feature points, uneven spatial distribution, and low regional sampling coverage, making it impossible to achieve refined deformation sampling across the entire high-altitude canyon and insufficient to support high-precision deformation monitoring of complex canyon slopes and other terrain features. Furthermore, the complex monitoring environment of high-altitude canyons results in prominent interferometric phase noise and sparse, discrete distribution of effective feature points. Conventional triangular mesh phase unwrapping algorithms fail to distinguish the phase quality differences among feature points, easily introducing low-quality, high-noise invalid points into the triangular mesh. This causes phase unwrapping errors to continuously spread and accumulate within the mesh, leading to problems such as phase unwrapping ambiguity and jumps, reducing the temporal phase unwrapping accuracy of high-altitude canyons. Ultimately, this results in poor stability and low reliability of deformation monitoring results for canyon slopes and other terrain features, failing to meet the high-precision deformation monitoring requirements of high-altitude canyon terrain. Summary of the Invention

[0003] To address the aforementioned problems in existing technologies, this invention proposes a high-precision deformation information extraction method and system for GB-SAR in high-altitude canyons. First, a GB-SAR time-series interferometric dataset of the study area is acquired. A three-dimensional atmospheric phase field is constructed based on an atmospheric tropospheric refractive index model to correct atmospheric distortion of the entangled phase. Then, a DEM (Digital Elevation Model) elevation stratification strategy is used to perform zonal compensation on the corrected phase. Finally, the coherence coefficient field is enhanced using polarization scattering power, fusing candidate PS (Persistent Scatterer, i.e., coherent point target) points selected by amplitude dispersion index with SBAS (Small Baseline) points selected by enhanced coherence coefficient. A high-density effective feature point set is constructed by selecting candidate points (i.e., a small baseline set). Based on the feature points, a triangular mesh is constructed and phase unwrapping is carried out using phase quality factor weighted constraints to obtain high-quality temporally unwrapped phase. Finally, a linear deformation model is constructed to complete the temporal deformation inversion, outputting the deformation rate and cumulative deformation information of each feature point. This enables high-precision surface deformation extraction under complex terrain of high-altitude canyons, thereby solving the problem of poor applicability of traditional uniform atmospheric assumption models. At the same time, it improves the ability to identify scattering points in vegetation areas and enhances the monitoring accuracy and reliability of deformation monitoring in high-altitude canyons.

[0004] To achieve the above objectives, the present invention specifically adopts the following technical solution: a method for high-precision deformation information extraction from high-altitude cold canyons using GB-SAR, comprising the following steps: S1. Obtain the GB-SAR time-series interferometric dataset for the high-altitude canyon study area. This dataset includes the time-series SAR (Synthetic Aperture Radar) amplitude field, entangled phase field, coherence coefficient field, polarization scattering power field, elevation matrix (DEM), and phase quality field. S2. Calculate the three-dimensional atmospheric phase field based on the atmospheric tropospheric refractive index model. Correct the entangled phase field based on this three-dimensional atmospheric phase field to obtain the corrected entangled phase field. S3. Divide the study area into multiple elevation-layered regions. Perform phase compensation on the corrected entangled phase field in each elevation-layered region to obtain an optimized entangled phase field. S4. Enhance the coherence coefficient field based on the polarization scattering power and the amplitude-dispersion... S5. Based on the dispersion index, PS candidate points are selected, and SBAS candidate points are selected based on the enhanced coherence coefficient threshold. The two types of candidate points are fused to obtain the effective feature point set and the number of effective feature points. S6. Based on the effective feature point set and the number of effective feature points, a triangular unit index matrix is ​​constructed. A weight vector is constructed using the phase mass field. The triangular mesh phase unwrapping function is called to obtain the temporal unwrapped phase of the effective feature points. S7. A linear deformation coefficient matrix is ​​constructed. The temporal unwrapped phase of the effective feature points is inverted to obtain the deformation rate sequence and cumulative deformation matrix of each effective feature point, thus completing the high-precision surface deformation information extraction of the high-altitude canyon.

[0005] In step S2, an elevation-related atmospheric phase is constructed based on an atmospheric tropospheric refractive index model. The wound phase field is then corrected in the wound phase domain to obtain the corrected wound phase field, specifically: S21: Using the reference thermodynamic temperature, subtract the product of the preset standard atmospheric thermodynamic temperature lapse rate and the elevation matrix DEM to obtain the thermodynamic temperature T at the elevation location; based on the atmospheric tropospheric refractive index model N=k1×P / T+k2×e×RH / T 2 Calculate the tropospheric refractive index modulus N, where k1 is the empirical coefficient of the first model, k2 is the empirical coefficient of the second model, P is the atmospheric static pressure, RH is the relative humidity, and e is the saturated water vapor pressure. S22: Multiply the atmospheric tropospheric refractive index modulus N by the atmospheric refractive index modulus conversion factor, and then multiply it by the phase distance conversion factor to obtain the absolute atmospheric phase; superimpose the phase offset on the absolute atmospheric phase, perform a modulo operation with 2π, and then subtract π from the result of the modulo operation to obtain the entangled domain atmospheric phase; S23: The atmospheric phase of the entangled domain is copied and extended along the time dimension to form a three-dimensional atmospheric phase field; the difference between the entangled phase field and the three-dimensional atmospheric phase field is obtained to obtain a phase difference field; the phase difference field is superimposed with a phase offset and modulo operation is performed with 2π, and then the result of the modulo operation is subtracted from π to obtain the corrected entangled phase field.

[0006] Specifically, S3 is: S31: After expanding the elevation matrix DEM into a one-dimensional array, calculate the maximum and minimum values ​​of the one-dimensional array, and divide the difference between the maximum and minimum values ​​by the preset number of elevation layers to obtain the equal elevation layer interval. S32: Construct an outer loop with sequence number I, and values ​​ranging from 1 to the preset elevation layer number. For each I, calculate the product of the difference between I and 1 and the isoelevation layer interval, and then add it to the minimum value in S31 to obtain the lower limit elevation; calculate the product of I and the isoelevation layer interval, and then add it to the minimum value in S31 to obtain the upper limit elevation; filter the pixels in the elevation matrix DEM that are greater than the lower limit elevation and less than or equal to the upper limit elevation to form a binary mask; remove the intervals in the binary mask where the number of effective pixels is lower than a preset threshold. S33: Construct an inner loop with the sequence number t, and the value is from 1 to the number of preset time-series interferograms. For each t, ​​read the two-dimensional phase frame data corresponding to index t in the corrected wound phase field, use the binary mask as an index to extract the corresponding data in the two-dimensional phase frame data as a set of phase sampling points, and then calculate the mean of the set of phase sampling points. S34: Calculate the difference between the set of phase sampling points and the mean of the set of phase sampling points, add the phase offset, and perform a modulo operation with 2π. Then subtract π from the result of the modulo operation to obtain the corrected set of phase sampling points. Update the two-dimensional phase frame data using the corrected set of phase sampling points and assign it to the corrected winding phase field. After the inner and outer loops are completed, use the corrected winding phase field as the optimized winding phase field.

[0007] Specifically, S4 is: S41: Solve for the minimum and maximum power values ​​of the polarization scattering power field along the two-dimensional spatial dimension. Subtract the minimum power value from the polarization scattering power field and divide it element-wise by the normalized denominator to obtain the normalized power. The normalized denominator is the difference between the maximum and minimum power values. Set a polarization weight coefficient, multiply the polarization weight coefficient by the normalized power, and then add a unit constant to construct a polarization gain factor. Perform a pixel-wise multiplication operation between the coherence coefficient field and the polarization gain factor to obtain the enhanced coherence coefficient field. S42: Obtain a three-dimensional temporal amplitude field by taking the modulus of the temporal SAR amplitude field, calculate the mean and standard deviation of the three-dimensional temporal amplitude field along the third dimension, divide the standard deviation and the mean element by element to obtain the amplitude dispersion index, mark the pixels with the amplitude dispersion index less than the preset amplitude dispersion index threshold as PS candidate points and generate a binary mask of PS candidate points. S43: The average coherence coefficient matrix is ​​obtained by averaging the enhanced coherence coefficient field along the third dimension. The average coherence coefficient matrix is ​​compared point by point with a preset coherence threshold. When the value of the average coherence coefficient matrix is ​​greater than the preset coherence threshold, the corresponding position mask is marked as a valid identifier, thereby obtaining the coherence coefficient mask. S44: Perform a logical OR fusion on the binary mask of the PS candidate points and the coherence coefficient mask to obtain the effective region mask. Traverse the effective region mask, extract the planar coordinates corresponding to all effective points to obtain the effective feature point set, and calculate the number of effective feature points in the effective feature point set.

[0008] Specifically, S5 is: S51: Extract the horizontal coordinate sequence and vertical coordinate sequence from the set of effective feature points and concatenate them horizontally to generate an effective point coordinate matrix. Based on the effective point coordinate matrix, perform Delaunay triangulation to obtain a triangular unit index matrix. S52: Initialize the time-series unwrapped phase matrix such that the first dimension of the time-series unwrapped phase matrix is ​​the number of effective feature points and the second dimension is the preset number of time-series interferograms; construct an outer loop with the sequence number t1 and a value from 1 to the preset number of time-series interferograms; construct a wrapped phase sequence and a phase quality factor sequence with a length equal to the number of effective feature points, and initialize them to zero. S53: Construct an inner loop with index p, where the value ranges from 1 to the number of effective feature points. For each p, read the horizontal coordinate sequence index and vertical coordinate sequence index of the p-th effective point. From the optimized winding phase field, store the horizontal coordinate sequence index, vertical coordinate sequence index, and t1 sampled phase at the p-th position of the winding phase sequence. From the phase quality field, store the horizontal coordinate sequence index, vertical coordinate sequence index, and t1 sampled quality factor at the p-th position of the phase quality factor sequence. S54: After the inner loop with sequence number p ends, the phase quality factor of the phase quality factor sequence is compared with the preset phase quality threshold point by point. Points with phase quality factors higher than the preset phase quality threshold are assigned 1, and the remaining points are assigned 0, to obtain a numerical point weight vector. S55: Based on the effective point coordinate matrix, triangular unit index matrix, entanglement phase sequence and numerical point weight vector, call the triangular mesh phase unwrapping function, take the output of the triangular mesh phase unwrapping function as the temporal unwrapping phase corresponding to the outer loop t1, and when the outer loop with the sequence number t1 ends, obtain the temporal unwrapping phase of the effective feature point.

[0009] The triangular mesh phase unwrapping function takes the effective point coordinate matrix, triangular element index matrix, wrapped phase sequence, and numerical point weight vector as input, specifically: S551: Obtain the number of rows in the effective point coordinate matrix and the number of rows in the triangular unit index matrix, construct a phase integer compensation sequence with a length equal to the number of rows in the effective point coordinate matrix and initialize it to zero; S552: Construct a loop with index i1 and value from 1 to the row number of the triangular unit index matrix. For each i1, extract the entire row of data of the i1th triangular unit index matrix to obtain a one-dimensional vertex index vector. Based on the one-dimensional vertex index vector, obtain the first vertex index, the second vertex index, and the third vertex index respectively. S553: ​​Obtain the numerical point weight vector values ​​corresponding to the first vertex index, the second vertex index, and the third vertex index. If any of the numerical point weight vector values ​​corresponding to the first vertex index, the second vertex index, and the third vertex index is less than the preset weight threshold, return to S552 and execute the next loop; if the numerical point weight vector values ​​corresponding to the first vertex index, the second vertex index, and the third vertex index are all greater than or equal to the preset weight threshold, execute S554. S554: Calculate the phase difference between the winding phase sequence corresponding to the second vertex index and the winding phase sequence corresponding to the first vertex index and use it as the first phase difference. Divide the first phase difference by 2π, round it down, and then multiply it by 2π to use it as the first phase compensation value. Add the first phase compensation value and the first phase difference to obtain the first normalized phase gradient. Calculate the phase difference between the winding phase sequence corresponding to the third vertex index and the winding phase sequence corresponding to the first vertex index and use it as the second phase difference. Divide the second phase difference by 2π, round it down, and then multiply it by 2π to obtain the second phase compensation value. Superimpose the second phase compensation value and the second phase difference to obtain the second normalized phase gradient. S555: Calculate the winding phase sequence corresponding to the first vertex index, subtract the winding phase sequence corresponding to the second vertex index, add the first phase difference, divide by 2π and round to obtain the phase integer compensation amount corresponding to the second vertex; calculate the winding phase sequence corresponding to the first vertex index, subtract the winding phase sequence corresponding to the third vertex index, add the second phase difference, divide by 2π and round to obtain the phase integer compensation amount corresponding to the third vertex. S556: After the loop with sequence number i1 is completed, the phase integer compensation sequence is updated. The updated phase integer compensation sequence is multiplied by 2π and then superimposed on the wrapped phase sequence as the output of the triangular mesh phase unwrapping function.

[0010] Specifically, S6 is: S61: Construct an equally spaced time reference sequence starting from 1 and ending at a preset number of time-series interferograms; construct a linear deformation coefficient matrix, where the first column of the linear deformation coefficient matrix is ​​a constant term basis vector and the second column is the equally spaced time reference sequence; construct a deformation rate sequence with a length equal to the number of effective feature points and initialize it to zero; construct a cumulative deformation matrix with a dimension equal to the number of effective feature points multiplied by the preset number of time-series interferograms and initialize it to zero; S62: Establish an outer loop with sequence number p1 and values ​​from 1 to the number of effective feature points. For each p1, extract the p1th time-series unwrapped phase and transpose it to obtain the time-series phase observation vector. Divide the time-series phase observation vector by the linear deformation coefficient matrix to obtain the point deformation parameter vector. Extract the second component of the point deformation parameter vector to obtain the phase linear change rate. Multiply the phase linear change rate by the radar wavelength and divide it by the first deformation coefficient to obtain the p1th deformation rate sequence. S63: Establish an inner loop with sequence number t2 and values ​​from 1 to the preset number of time-series interferograms. For each t2, calculate the product of the second component of the point deformation parameter vector and the t2th equally spaced time reference sequence, and then superimpose the first component of the point deformation parameter vector as the fitting phase. Multiply the fitting phase by the radar wavelength and divide it by the second deformation coefficient to obtain the cumulative deformation matrix elements under the p1th outer loop and the t2th inner loop. S64: After the inner loop with sequence number t2 and the outer loop with sequence number p1 are completed, the deformation rate sequence and the cumulative deformation matrix are obtained.

[0011] This invention also proposes a high-precision deformation information extraction system for high-altitude cold canyons using GB-SAR, which is used to implement the aforementioned method for high-precision deformation information extraction from high-altitude cold canyons using GB-SAR. The high-precision deformation information extraction system for high-altitude cold canyons using GB-SAR includes a data acquisition module, a data analysis module, and a result output module. The data acquisition module is connected to the data analysis module, and the data analysis module is connected to the result output module. The data acquisition module is used to acquire GB-SAR time-series interferometric datasets of the high-altitude cold canyon study area. The data analysis module is used to obtain high-precision surface deformation information of the high-altitude cold canyon. The result output module is used to output the high-precision surface deformation information of the high-altitude cold canyon.

[0012] The beneficial technical effects of this invention compared with the prior art are as follows: (1) Based on the physical model of atmospheric tropospheric refractive index, the entanglement phase correction is realized, and then the phase compensation is carried out by the elevation layering and partitioning of DEM and other elevations to gradually weaken the elevation-related atmospheric phase distortion caused by the dramatic elevation fluctuations in the high-altitude valleys, thus solving the problem of poor applicability of the traditional uniform atmospheric assumption model. (2) The coherence coefficient field is enhanced by polarization scattering power, and PS candidate points are screened by combining amplitude dispersion index and SBAS candidate points after enhancement and fusion, thereby increasing the effective feature point density of vegetation-covered areas in alpine canyons, enhancing the ability to identify scattering points in vegetation areas, increasing the spatial sampling density of alpine canyon deformation monitoring, and improving the reliability of alpine canyon deformation monitoring. (3) Based on the phase mass field, the point weight vector is generated, and the triangular unit containing low quality points is automatically removed in the triangular mesh. The phase integer compensation is calculated by the phase gradient of the triangular unit vertex, which suppresses the spread of unwrapping error caused by noise and improves the monitoring accuracy of deformation monitoring in high-altitude canyons. Attached Figure Description

[0013] To more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are merely exemplary, and those skilled in the art can derive other embodiments based on the provided drawings without creative effort.

[0014] Figure 1 This is a simplified flowchart of a high-precision deformation information extraction method for high-altitude cold canyons using GB-SAR, according to the present invention. Figure 2 This is a simplified flowchart of the method for obtaining the effective feature point set and the number of effective feature points according to the present invention. Figure 3 This is a flowchart illustrating the workflow of the triangular mesh phase unwrapping function of the present invention. Figure 4 This is a simplified structural diagram of a GB-SAR high-precision deformation information extraction system for high-altitude cold canyons according to the present invention. Detailed Implementation

[0015] The technical solutions of the embodiments of the present invention 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. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0016] The concepts involved in this application will first be described with reference to the accompanying drawings. It should be noted that the following descriptions of various concepts are only for the purpose of making the content of this application easier to understand and do not constitute a limitation on the scope of protection of this application; furthermore, the embodiments and features in the embodiments of this application can be combined with each other unless otherwise specified. This application will now be described in detail with reference to the accompanying drawings and embodiments.

[0017] Refer to the instruction manual. Figure 1 The present invention proposes a method for high-precision deformation information extraction from high-altitude cold canyons using GB-SAR, comprising the following steps: S1. Obtain the GB-SAR time-series interferometric dataset for the high-altitude canyon study area. The dataset includes the time-series SAR amplitude field, entangled phase field, coherence coefficient field, polarization scattering power field, elevation matrix DEM, and phase quality field. Except for the elevation matrix DEM, which is two-dimensional, all other datasets are three-dimensional. S2. Calculate the three-dimensional atmospheric phase field based on the atmospheric tropospheric refractive index model. Correct the entangled phase field based on the three-dimensional atmospheric phase field to obtain the corrected entangled phase field. S3. Divide the study area into multiple elevation-layered regions. Perform iterative phase compensation on the corrected entangled phase field in each elevation-layered region. After the maximum number of iterations, obtain the optimized entangled phase field that can suppress atmospheric disturbances in the vertical stratification of the canyon. S4. Fully utilize polarization information... To improve the low coherence coefficient and poor signal-to-noise ratio in high-undulation canyon areas, the coherence coefficient field is enhanced by polarization scattering power. PS candidate points are selected based on the amplitude dispersion index, and SBAS candidate points are selected based on the enhanced coherence coefficient threshold. The two types of candidate points are fused to obtain the effective feature point set and the number of effective feature points. S5: A triangular element index matrix is ​​constructed based on the effective feature point set and the number of effective feature points. A weight vector is constructed using the phase mass field, and the triangular mesh phase unwrapping function is called to obtain the temporal unwrapped phase of the effective feature points. S6: A linear deformation coefficient matrix is ​​constructed, and the temporal unwrapped phase of the effective feature points is inverted. The residual topographic phase and deformation phase are separated to obtain the deformation rate sequence and cumulative deformation matrix of each effective feature point, completing the high-precision extraction of surface deformation information in high-altitude canyons.

[0018] In S2, an elevation-related atmospheric phase is constructed based on the atmospheric tropospheric refractive index model. The wound phase field is then corrected in the wound phase domain to obtain the corrected wound phase field, specifically: S21: The terrain of the high-altitude canyon is dramatically undulating, and the air temperature varies significantly with elevation. Based on a baseline thermodynamic temperature of 276K, an international standard atmospheric temperature lapse model is used. The two-dimensional thermodynamic temperature field T corresponding to each elevation pixel in the study area is obtained by subtracting the standard atmospheric thermodynamic temperature lapse rate of 0.0065K / m from the baseline thermodynamic temperature and multiplying it by the pixel-by-pixel product of the elevation matrix DEM. Based on the atmospheric tropospheric refractive index model N=k1×P / T+k2×e×RH / T 2 Calculate the tropospheric refractive index modulus N, where k1 is the empirical coefficient of the first model (77.7 K / hPa) and k2 is the empirical coefficient of the second model (373454 K). 2 P is atmospheric static pressure (580 hPa), RH is relative humidity (0.6, dimensionless), and e is saturated water vapor pressure (approximately 7.5 hPa). S22: Electromagnetic waves experience a propagation delay when passing through the troposphere, which is determined by the atmospheric refractive index modulus. Multiplying the two-dimensional tropospheric refractive index modulus field N by the atmospheric refractive index modulus conversion factor (10...) -6 The atmospheric propagation delay field is obtained; then, the delay field is multiplied by the phase-distance conversion coefficient (4π / 0.056) to calculate the two-dimensional absolute atmospheric phase field. To map the continuous absolute atmospheric phase to the entangled phase interval, the phase offset (π) is superimposed on the absolute atmospheric phase, a modulo 2π operation is performed, and then π is subtracted from the result to constrain the phase to the interval [-π, π], thus obtaining the two-dimensional entangled domain atmospheric phase. S23: Since the GB-SAR temporal interferometric dataset is a three-dimensional structure, the atmospheric phase of the two-dimensional winding domain cannot be directly calculated with the original three-dimensional winding phase field. Therefore, the atmospheric phase of the two-dimensional winding domain is copied and extended along the temporal dimension to construct a three-dimensional atmospheric phase field that is completely matched with the dimension of the original dataset. The original winding phase field is subtracted from the three-dimensional atmospheric phase field pixel by pixel to obtain the phase difference field. Similarly, the phase offset (π) is superimposed on the phase difference field, and after performing a modulo 2π operation, π is subtracted to finally obtain the corrected winding phase field after atmospheric phase compensation.

[0019] S3 specifically refers to: S31: After expanding the elevation matrix DEM into a one-dimensional array, calculate the maximum and minimum values ​​of the one-dimensional array. Divide the difference between the maximum and minimum values ​​by the preset number of elevation layers (e.g., 4 layers) to obtain the equal elevation layer interval step size, and complete the calculation of the basic parameters of the elevation layer. S32: Construct an outer loop with index I and a traversal range from 1 to a preset elevation layer number. For each I, calculate the product of the difference between I and 1 and the equal elevation layer interval, and then add it to the minimum value in S31 to obtain the lower limit elevation. Calculate the product of I and the equal elevation layer interval, and then add it to the minimum value in S31 to obtain the upper limit elevation. Filter pixels in the elevation matrix DEM that are greater than the lower limit elevation and less than or equal to the upper limit elevation to form a binary mask. To avoid unreliable statistical values ​​due to insufficient pixels, perform effectiveness screening on each elevation layer mask, removing elevation intervals where the number of effective pixels in the mask is lower than a preset threshold (preferably 10), and discarding layers with insufficient samples. S33: Construct an inner loop with sequence number t and a traversal range from 1 to the preset number of time-series interferograms (e.g., 22). For each t, ​​read the two-dimensional phase frame data corresponding to index t in the corrected wound phase field, use a binary mask as an index to extract the corresponding data in the two-dimensional phase frame data as a set of phase sampling points, and then calculate the arithmetic mean of the set of phase sampling points to characterize the statistical level of the phase under the elevation layer. S34: Calculate the difference between the set of phase sampling points and the mean of the set of phase sampling points, then add the phase offset, and perform a modulo operation with 2π. Subtract π from the result of the modulo operation to obtain the corrected set of phase sampling points. Use the corrected set of phase sampling points to write back and update the pixel values ​​of the corresponding two-dimensional phase frame intra-mask region, and assign the updated phase frame to the corrected wound phase field. Complete all elevation layer outer loops and all temporal frame intra-loops in sequence, and finally output the optimized wound phase field.

[0020] Instruction manual attached Figure 2 The basic steps of S4 are summarized below. S4 is as follows: S41: To fully utilize polarization scattering information, the polarization scattering power field is normalized. The minimum and maximum power values ​​of the polarization scattering power field are calculated along the two-dimensional spatial dimension. The minimum power value is subtracted from the polarization scattering power field, and then the normalized power is divided element-wise by the normalized denominator to obtain the normalized power. The normalized denominator is the difference between the maximum and minimum power values. A polarization weighting coefficient (0.55) is set, and the polarization weighting coefficient is multiplied by the normalized power and then superimposed with a unit constant to construct the polarization gain factor. The coherence coefficient field is multiplied pixel-wise by the polarization gain factor to obtain the enhanced coherence coefficient field, thereby achieving weighted enhancement of the coherence of high scattering intensity pixels. S42: Obtain the three-dimensional temporal amplitude field by taking the modulus of the temporal SAR amplitude field. Calculate the mean and standard deviation of the three-dimensional temporal amplitude field along the third dimension, i.e., the time dimension. Divide the standard deviation from the mean element-wise to obtain the amplitude dispersion index. This index characterizes the degree of fluctuation of the pixel backscattering intensity over time; the smaller the index, the more stable the scattering. Mark pixels with an amplitude dispersion index less than the preset amplitude dispersion index threshold (0.05) as PS candidate points and generate a binary mask for PS candidate points. S43: The mean coherence coefficient matrix is ​​obtained by averaging the enhanced coherence coefficient field along the third dimension. The mean coherence coefficient matrix is ​​compared point by point with the preset coherence threshold (0.22). When the value of the mean coherence coefficient matrix is ​​greater than the preset coherence threshold, the corresponding position mask is marked as valid, and the other positions are marked as invalid, thus obtaining the coherence coefficient mask. S44: The binary mask of PS candidate points and the coherence coefficient mask are fused by logical OR operation. As long as the PS candidate condition or coherence condition is met, it is determined to be a valid region, and a comprehensive valid region mask is obtained. The valid region mask is traversed, and the planar spatial coordinates corresponding to all marked valid pixels are extracted to form a set of valid feature points. The total number of feature points in the set of valid feature points is counted to provide a reliable point target basis for subsequent phase unwrapping and deformation time series inversion.

[0021] S5 specifically refers to: S51: Extract the horizontal and vertical coordinate sequences from the effective feature point set and perform horizontal splicing to generate an effective point coordinate matrix. Based on the effective point coordinate matrix, perform two-dimensional Delaunay triangulation to obtain a triangular element index matrix. Remove invalid triangular elements that exceed the research boundary or have malformed narrow triangles in the triangulation results to complete the triangular mesh topology construction. S52: Initialize the time-series unwrapped phase matrix, such that the first dimension of the time-series unwrapped phase matrix is ​​the number of effective feature points, and the second dimension is the preset number of time-series interferograms; construct an outer loop with the index t1 and a value from 1 to the preset number of time-series interferograms, construct a wrapped phase sequence and a phase quality factor sequence with a length equal to the number of effective feature points, and initialize them to zero; the wrapped phase sequence is used to store the wrapped phase observation values ​​of all feature points of the current interferogram, and the phase quality factor sequence characterizes the reliability of the phase observation at each point.

[0022] S53: Construct a nested inner loop with index p, traversing from 1 to the number of valid feature points. For each p, read the x-coordinate sequence index and y-coordinate sequence index of the p-th valid point. From the optimized winding phase field, store the x-coordinate sequence index, y-coordinate sequence index, and t1 sampled phase at the p-th position of the winding phase sequence. From the phase quality field, store the x-coordinate sequence index, y-coordinate sequence index, and t1 sampled quality factor at the p-th position of the phase quality factor sequence. After traversing all feature points, obtain the complete point phase and quality information of the current interferogram. S54: After the inner loop with sequence number p ends, the phase quality factor of the phase quality factor sequence is compared with the preset phase quality threshold (0.25) point by point. Points with phase quality factors higher than the preset phase quality threshold are assigned 1, and the rest are assigned 0, to obtain a numerical point weight vector used to suppress low quality point interference during phase unwrapping. S55: Based on the effective point coordinate matrix, triangular unit index matrix, entangled phase sequence and numerical point weight vector, the triangular mesh phase unwrapping function is called. The output of the triangular mesh phase unwrapping function is used as the temporal unwrapping phase corresponding to the outer loop t1. When the outer loop with the sequence number t1 ends, the temporal unwrapping phase of the effective feature point is obtained.

[0023] In step S55 above, the temporal phase unwrapping is achieved using a triangular mesh phase unwrapping function, as detailed in the attached manual. Figure 3 This demonstrates the working steps of the triangular mesh phase unwrapping function. The function takes the effective point coordinate matrix, the triangular element index matrix, the wrapped phase sequence, and the numerical point weight vector as inputs. Specifically: S551: Obtain the number of rows in the valid point coordinate matrix and the number of rows in the triangular unit index matrix, construct a phase integer compensation sequence with a length equal to the number of rows in the valid point coordinate matrix and initialize it to zero; S552: Construct a loop with index i1, where the value ranges from 1 to the row number of the triangular unit index matrix. For each i1, extract the entire row of data from the i1th triangular unit index matrix to obtain a one-dimensional vertex index vector containing the three vertex numbers. Based on the one-dimensional vertex index vector, obtain the first vertex index, the second vertex index, and the third vertex index respectively. By accurately locating the three feature points of the triangular unit, the positional basis can be provided for the subsequent phase operation and quality judgment of the single triangular unit. S553: ​​Based on the extracted three vertex indices, match and read the corresponding weight values ​​in the numerical point weight vector to determine the validity of the current triangular unit. Set the preset weight threshold to 0.1 and execute the unit filtering mechanism: if the weight value of any of the three vertices is less than 0.1, the triangular unit is determined to contain low-coherence, high-noise invalid points, which are prone to causing phase unwrapping error propagation. The current unit is skipped directly, and the process returns to S552 to execute the next round of the loop; if the weight values ​​of all three vertices are greater than or equal to 0.1, the point data of the triangular unit is determined to be reliable and the phase quality meets the standard. It can participate in phase gradient calculation and integer compensation operation, and the process proceeds to S554. S554: Based on the phase data of the three vertices of a qualified triangular element, a bidirectional phase gradient normalization calculation is performed. Taking the first vertex as the reference vertex, the phase difference between the second and first vertices is first calculated to obtain the first phase difference. This phase difference is then scaled by 2π, and by dividing by 2π, rounding, and then multiplying by 2π, a first phase compensation value is obtained to compensate for decimal errors. This compensation value is then superimposed on the original phase difference to eliminate the jump deviation caused by phase winding, resulting in the first normalized phase gradient. Similarly, the phase difference between the third and first vertices is calculated to obtain the second phase difference. Through the same compensation operation logic, the second normalized phase gradient is obtained, achieving the standardization correction of the phase gradient within the triangular element. S555: Calculate the winding phase sequence corresponding to the first vertex index, subtract the winding phase sequence corresponding to the second vertex index, add the first phase difference, divide by 2π and round to obtain the phase integer compensation amount corresponding to the second vertex; calculate the winding phase sequence corresponding to the first vertex index, subtract the winding phase sequence corresponding to the third vertex index, add the second phase difference, divide by 2π and round to obtain the phase integer compensation amount corresponding to the third vertex. S556: After the loop with index i1 is completed, the phase integer compensation sequence is updated. The updated phase integer compensation sequence is multiplied by 2π and then superimposed on the wrapped phase sequence as the output of the triangular mesh phase unwrapping function.

[0024] S6 specifically refers to: S61: Construct an equally spaced time reference sequence starting from 1 and ending at a preset number of time-series interferograms. This sequence is used to characterize the relative observation time corresponding to each interferogram, providing a time dimension independent variable for linear deformation time-series fitting. Construct a linear deformation coefficient matrix, assigning all values ​​of 1 to the first column of the matrix as constants to form the model constant term basis vector, used to fit the initial fixed phase shift of the points. Substitute the equally spaced time reference sequence into the second column to form a time-varying basis vector, used to characterize the linear change trend of the phase over time. Construct a deformation rate sequence with a length equal to the number of effective feature points and initialize it to zero, used to store the final linear annual deformation rate of each point. Construct a cumulative deformation matrix with a dimension equal to the number of effective feature points multiplied by the preset number of time-series interferograms and initialize it to zero, used to store the cumulative surface deformation point by point and time by time, ensuring that the data storage space for subsequent point-by-point traversal operations is independent and the initial state is unbiased. S62: Establish an outer loop with index p1, where the value ranges from 1 to the number of effective feature points. For each p1, extract the p1th time-series unwrapped phase and transpose it to generate a one-dimensional time-series phase observation vector, which serves as the dependent variable for linear fitting. Perform a left division operation on the time-series phase observation vector using the preset linear deformation coefficient matrix, and obtain the point deformation parameter vector through least-squares fitting. This vector contains two core parameters: a fixed phase offset and a temporal phase change rate. Extract the second component of the deformation parameter vector, which represents the linear phase change rate of that feature point, characterizing the continuous phase change trend over time. By combining the physical parameters observed by radar, the physical quantity conversion is completed. The linear rate of phase change is multiplied by the radar wavelength and then divided by the preset first deformation coefficient (365×1000 / 4π) to complete the conversion of phase change to surface deformation. At the same time, the time scale and length scale are normalized to obtain the linear deformation rate of the point line of sight in millimeters per year. The value is then synchronously updated to the corresponding position in the deformation rate sequence. The first deformation coefficient is introduced with a 365-day time factor to realize the scale conversion of the phase change rate to the annual average deformation rate.

[0025] S63: Establish an inner loop with sequence number t2 and values ​​ranging from 1 to the preset number of time-series interferograms. For each t2, calculate the product of the second component of the point deformation parameter vector and the t2th equally spaced time reference sequence, and then superimpose the first component of the point deformation parameter vector as the fitted phase. Multiply the fitted phase by the radar wavelength and divide by the second deformation coefficient (1000 / 4π) to accurately convert the phase fitting value into line-of-sight cumulative deformation in millimeters, obtaining the cumulative deformation matrix elements under the p1th outer loop and the t2th inner loop. Here, the second deformation coefficient is the basic conversion coefficient from phase to millimeter cumulative deformation, and does not need to include time scale conversion.

[0026] S64: After completing the entire traversal of the inner loop (t2) and the outer loop (p1), the linear deformation parameter fitting, physical quantity conversion, and temporal cumulative deformation assignment for all valid feature points in the entire domain are completed. The final output is a complete deformation rate sequence and a cumulative deformation matrix. The deformation rate sequence represents the long-term stable annual average linear deformation characteristics of each feature point, and the cumulative deformation matrix represents the surface cumulative deformation state at each observation time relative to the initial time, providing a complete temporal deformation dataset for subsequent spatiotemporal deformation analysis, anomalous subsidence identification, and deformation residual calculation.

[0027] This invention also proposes a high-precision deformation information extraction system for high-altitude cold canyons using GB-SAR, for performing the aforementioned high-precision deformation information extraction method for high-altitude cold canyons using GB-SAR. (See attached specification.) Figure 4 The GB-SAR high-precision deformation information extraction system for high-altitude cold canyons includes a data acquisition module, a data analysis module, and a result output module; the output end of the data acquisition module is connected to the input end of the data analysis module, and the output end of the data analysis module is connected to the input end of the result output module.

[0028] The data acquisition module is used to collect GB-SAR time-series interferometric datasets for the high-altitude canyon study area. The time-series interferometric datasets contain multiple time-series interferograms, phase quality information, radar system parameters, observation geometric parameters, and basic geographic information of the study area. It completes the import of raw observation data, format parsing, and data preprocessing, providing a basic data source for subsequent phase analysis. The high-altitude canyon area has dramatic topographic relief and significant spatial and temporal differences in the atmosphere. The module supports batch reading of long-term continuous observation time-series data and removal of abnormal data.

[0029] The data analysis module is used to perform phase error suppression, feature point optimization, triangular mesh phase unwrapping, and temporal deformation parameter inversion on the received GB-SAR time-series interferometric dataset step by step, and calculate high-precision surface deformation information of the high-altitude canyon. The data analysis module includes an atmospheric phase correction unit, an elevation layer compensation unit, a feature point screening unit, a phase unwrapping unit, and a temporal deformation inversion unit, and the above units are coupled in sequence.

[0030] The atmospheric phase correction unit is used to separate and eliminate large-scale atmospheric phase errors in time-series interferograms, targeting the characteristics of vertically stratified atmospheric disturbances in high-altitude canyons, thereby reducing phase noise introduced by uneven vertical distribution of temperature and water vapor. The elevation layered compensation unit is used to combine canyon topographic elevation information to compensate for phase deviation caused by topographic undulation and elevation residual, and further optimize the interferometric phase field. The feature point screening unit is used to screen stable and effective feature points based on indicators such as coherence coefficient and phase signal-to-noise ratio, construct an effective feature point set, and remove noise points with low coherence and easy loss of coherence. The phase unwrapping unit is used to construct a Delaunay triangular mesh for the optimized wrapped phase field. Combined with the point phase quality weight, it realizes phase integer compensation and phase unwrapping of the triangular mesh interferogram per cycle and outputs the effective feature point time-series unwrapped phase. The temporal deformation inversion unit is used to construct a linear deformation least squares model based on the temporal unwrapped phase, invert the linear deformation rate of each effective feature point and the cumulative deformation at each time step, and complete the conversion of phase observation values ​​to physical deformation values.

[0031] The results output module is used to receive high-precision surface deformation information of high-altitude canyons from the data analysis module, realize the storage, visualization rendering, and data format export of linear deformation rate sequences and time-series cumulative deformation matrices; it supports raster mapping of deformation results and output of point time-series curves, which facilitates the identification and spatiotemporal evolution analysis of deformation hazards on high-altitude canyon slopes and rock masses.

[0032] The embodiments and / or implementation methods described above are merely preferred embodiments and / or implementation methods for implementing the technology of the present invention, and are not intended to limit the implementation methods of the technology of the present invention in any way. Any person skilled in the art can make some modifications or alterations to other equivalent embodiments without departing from the scope of the technical means disclosed in the content of the present invention, but they should still be regarded as the technology or embodiments that are substantially the same as the present invention.

[0033] This document uses specific examples to illustrate the principles and implementation methods of this application. The descriptions of the above embodiments are only for the purpose of helping to understand the methods and core ideas of this application. The above descriptions are only preferred embodiments of this application. It should be noted that due to the limitations of written expression, while there are objectively infinite specific structures, those skilled in the art can make several improvements, modifications, or changes without departing from the principles of this application, and can also combine the above technical features in an appropriate manner. These improvements, modifications, changes, or combinations, or the direct application of the inventive concept and technical solution to other situations without modification, should all be considered within the scope of protection of this application.

Claims

1. A method for high-precision deformation information extraction from high-altitude cold canyons using GB-SAR, characterized in that, Includes the following steps: S1. Obtain the GB-SAR time-series interferometric dataset of the high-altitude canyon study area. The dataset includes the time-series SAR amplitude field, the entangled phase field, the coherence coefficient field, the polarization scattering power field, the elevation matrix DEM, and the phase quality field. S2. Calculate the three-dimensional atmospheric phase field based on the atmospheric tropospheric refractive index model, and correct the winding phase field based on the three-dimensional atmospheric phase field to obtain the corrected winding phase field. S3. Divide the study area into multiple elevation-layered regions, and perform phase compensation on the corrected winding phase field in each elevation-layered region to obtain an optimized winding phase field. S4. The coherence coefficient field is enhanced based on the polarization scattering power. PS candidate points are selected based on the amplitude dispersion index and SBAS candidate points are selected based on the enhanced coherence coefficient threshold. The two types of candidate points are fused to obtain the effective feature point set and the number of effective feature points. S5. Construct a triangular unit index matrix based on the set of effective feature points and the number of effective feature points, construct a weight vector using the phase mass field, and call the triangular mesh phase unwrapping function to obtain the temporal unwrapped phase of the effective feature points; S6. Construct a linear deformation coefficient matrix, and invert the temporal unwrapped phase of the effective feature points to obtain the deformation rate sequence and cumulative deformation matrix of each effective feature point, thus completing the high-precision extraction of surface deformation information in the high-altitude canyon.

2. The method for high-precision deformation information extraction from high-altitude cold canyons using GB-SAR according to claim 1, characterized in that, In step S2, an elevation-related atmospheric phase is constructed based on an atmospheric tropospheric refractive index model. The wound phase field is then corrected in the wound phase domain to obtain the corrected wound phase field, specifically: S21: Using the reference thermodynamic temperature, subtract the product of the preset standard atmospheric thermodynamic temperature lapse rate and the elevation matrix DEM to obtain the thermodynamic temperature T at the elevation location; based on the atmospheric tropospheric refractive index model N=k1×P / T+k2×e×RH / T 2 Calculate the tropospheric refractive index modulus N, where k1 is the empirical coefficient of the first model, k2 is the empirical coefficient of the second model, P is the atmospheric static pressure, RH is the relative humidity, and e is the saturated water vapor pressure. S22: Multiply the atmospheric tropospheric refractive index modulus N by the atmospheric refractive index modulus conversion factor, and then multiply it by the phase distance conversion factor to obtain the absolute atmospheric phase; superimpose the phase offset on the absolute atmospheric phase, perform a modulo operation with 2π, and then subtract π from the result of the modulo operation to obtain the entangled domain atmospheric phase; S23: The atmospheric phase of the entangled domain is copied and extended along the time dimension to form a three-dimensional atmospheric phase field; the difference between the entangled phase field and the three-dimensional atmospheric phase field is obtained to obtain a phase difference field; the phase difference field is superimposed with a phase offset and modulo operation is performed with 2π, and then the result of the modulo operation is subtracted from π to obtain the corrected entangled phase field.

3. The method for high-precision deformation information extraction from high-altitude cold canyons using GB-SAR according to claim 1, characterized in that, Specifically, S3 is... S31: After expanding the elevation matrix DEM into a one-dimensional array, calculate the maximum and minimum values ​​of the one-dimensional array, and divide the difference between the maximum and minimum values ​​by the preset number of elevation layers to obtain the equal elevation layer interval. S32: Construct an outer loop with serial number I and a value from 1 to the preset elevation layer number. For each I, calculate the product of the difference between I and 1 and the equal elevation layer interval, and then add it to the minimum value in S31 to obtain the lower limit elevation. Calculate the product of I and the equal elevation layer interval, and then add it to the minimum value in S31 to obtain the upper limit elevation; Pixels in the elevation matrix DEM that are greater than the lower limit elevation and less than or equal to the upper limit elevation are selected to form a binary mask; intervals in the binary mask where the number of effective pixels is less than a preset threshold are removed. S33: Construct an inner loop with the sequence number t, which takes values ​​from 1 to the preset number of time-series interferograms. For each t, ​​read the two-dimensional phase frame data corresponding to index t in the corrected wound phase field, use the binary mask as an index to extract the corresponding data in the two-dimensional phase frame data as a set of phase sampling points, and then calculate the mean of the set of phase sampling points. S34: Calculate the difference between the set of phase sampling points and the mean of the set of phase sampling points, add the phase offset, and perform a modulo operation with 2π. Then subtract π from the result of the modulo operation to obtain the corrected set of phase sampling points. Update the two-dimensional phase frame data using the corrected set of phase sampling points and assign it to the corrected winding phase field. After the inner and outer loops are completed, use the corrected winding phase field as the optimized winding phase field.

4. The method for high-precision deformation information extraction from high-altitude cold canyons using GB-SAR according to claim 1, characterized in that, Specifically, S4 is: S41: Solve for the minimum and maximum power values ​​of the polarization scattering power field along the two-dimensional spatial dimension. Subtract the minimum power value from the polarization scattering power field and divide it element-wise by the normalized denominator to obtain the normalized power. The normalized denominator is the difference between the maximum power value and the minimum power value. A polarization weighting coefficient is set, and the polarization weighting coefficient is multiplied by the normalized power and then superimposed with a unit constant to construct a polarization gain factor; the coherence coefficient field is multiplied pixel by pixel with the polarization gain factor to obtain the enhanced coherence coefficient field; S42: Obtain a three-dimensional temporal amplitude field by taking the modulus of the temporal SAR amplitude field, calculate the mean and standard deviation of the three-dimensional temporal amplitude field along the third dimension, divide the standard deviation and the mean element by element to obtain the amplitude dispersion index, mark the pixels with the amplitude dispersion index less than the preset amplitude dispersion index threshold as PS candidate points and generate a PS candidate point binary mask. S43: The average coherence coefficient matrix is ​​obtained by averaging the enhanced coherence coefficient field along the third dimension. The average coherence coefficient matrix is ​​compared point by point with the preset coherence threshold. When the value of the average coherence coefficient matrix is ​​greater than the preset coherence threshold, the corresponding position mask is marked as a valid identifier, thereby obtaining the coherence coefficient mask. S44: Perform a logical OR fusion on the binary mask of the PS candidate points and the coherence coefficient mask to obtain the effective region mask. Traverse the effective region mask, extract the planar coordinates corresponding to all effective points to obtain the effective feature point set, and calculate the number of effective feature points in the effective feature point set.

5. The method for high-precision deformation information extraction from high-altitude cold canyons using GB-SAR according to claim 1, characterized in that, Specifically, S5 is: S51: Extract the horizontal coordinate sequence and vertical coordinate sequence from the set of effective feature points and concatenate them horizontally to generate an effective point coordinate matrix. Based on the effective point coordinate matrix, perform Delaunay triangulation to obtain a triangular unit index matrix. S52: Initialize the time-series unwrapped phase matrix such that the first dimension of the time-series unwrapped phase matrix is ​​the number of effective feature points and the second dimension is the number of preset time-series interferograms; construct an outer loop with the sequence number t1 and a value from 1 to the number of preset time-series interferograms; construct a wrapped phase sequence and a phase quality factor sequence with a length equal to the number of effective feature points, and initialize them to zero. S53: Construct an inner loop with sequence number p, and values ​​from 1 to the number of effective feature points. For each p, read the horizontal coordinate sequence index and vertical coordinate sequence index of the p-th effective point. From the optimized winding phase field, the horizontal coordinate sequence index, the vertical coordinate sequence index, and the t1 sampled phase are stored together at the p-th position of the winding phase sequence. From the phase quality field, the horizontal coordinate sequence index, the vertical coordinate sequence index, and the t1 sampled quality factor are stored together at the p-th position of the phase quality factor sequence. S54: After the inner loop with sequence number p ends, the phase quality factor of the phase quality factor sequence is compared with the preset phase quality threshold point by point. Points with phase quality factors higher than the preset phase quality threshold are assigned 1, and the remaining points are assigned 0, to obtain a numerical point weight vector. S55: Based on the effective point coordinate matrix, triangular unit index matrix, entanglement phase sequence and numerical point weight vector, call the triangular mesh phase unwrapping function, take the output of the triangular mesh phase unwrapping function as the temporal unwrapping phase corresponding to the outer loop t1, and when the outer loop with the sequence number t1 ends, obtain the temporal unwrapping phase of the effective feature point.

6. The method for high-precision deformation information extraction from high-altitude cold canyons using GB-SAR according to claim 5, characterized in that, The triangular mesh phase unwrapping function takes the effective point coordinate matrix, triangular element index matrix, wrapped phase sequence, and numerical point weight vector as input, specifically: S551: Obtain the number of rows in the effective point coordinate matrix and the number of rows in the triangular unit index matrix, construct a phase integer compensation sequence with a length equal to the number of rows in the effective point coordinate matrix and initialize it to zero; S552: Construct a loop with index i1 and value from 1 to the row number of the triangular unit index matrix. For each i1, extract the entire row of data of the i1th triangular unit index matrix to obtain a one-dimensional vertex index vector. Based on the one-dimensional vertex index vector, obtain the first vertex index, the second vertex index, and the third vertex index respectively. S553: ​​Obtain the numerical point weight vector values ​​corresponding to the first vertex index, the second vertex index, and the third vertex index. If any of the numerical point weight vector values ​​corresponding to the first vertex index, the second vertex index, and the third vertex index is less than the preset weight threshold, return to S552 and execute the next loop; if the numerical point weight vector values ​​corresponding to the first vertex index, the second vertex index, and the third vertex index are all greater than or equal to the preset weight threshold, execute S554. S554: Calculate the phase difference between the winding phase sequence corresponding to the second vertex index and the winding phase sequence corresponding to the first vertex index and use it as the first phase difference. Divide the first phase difference by 2π, round it down, and then multiply it by 2π to get the first phase compensation value. Add the first phase compensation value and the first phase difference to obtain the first normalized phase gradient. Calculate the phase difference between the winding phase sequence corresponding to the third vertex index and the winding phase sequence corresponding to the first vertex index and use it as the second phase difference. Divide the second phase difference by 2π, round it down, and then multiply it by 2π to obtain the second phase compensation value. Add the second phase compensation value and the second phase difference to obtain the second normalized phase gradient. S555: Calculate the winding phase sequence corresponding to the first vertex index, subtract the winding phase sequence corresponding to the second vertex index, add the first phase difference, divide by 2π and round to obtain the phase integer compensation amount corresponding to the second vertex. The winding phase sequence corresponding to the first vertex index is calculated, the winding phase sequence corresponding to the third vertex index is subtracted, the second phase difference is added, and then the result is divided by 2π and rounded to obtain the phase integer compensation amount corresponding to the third vertex. S556: After the loop with sequence number i1 is completed, the phase integer compensation sequence is updated. The updated phase integer compensation sequence is multiplied by 2π and then superimposed on the winding phase sequence as the output of the triangular mesh phase unwrapping function.

7. The method for high-precision deformation information extraction from high-altitude cold canyons using GB-SAR according to claim 1, characterized in that, Specifically, S6 is: S61: Construct an equally spaced time reference sequence starting from 1 and ending at a preset number of time-series interferograms; construct a linear deformation coefficient matrix, wherein the first column of the linear deformation coefficient matrix is ​​a constant term basis vector and the second column is the equally spaced time reference sequence; Construct a deformation rate sequence with a length equal to the number of effective feature points and initialize it to zero; construct a cumulative deformation matrix with a dimension equal to the number of effective feature points multiplied by a preset number of temporal interferograms and initialize it to zero; S62: Establish an outer loop with sequence number p1 and values ​​from 1 to the number of effective feature points. For each p1, extract the p1th time-series unwrapped phase and transpose it to obtain the time-series phase observation vector. Divide the time-series phase observation vector by the linear deformation coefficient matrix to obtain the point deformation parameter vector. Extract the second component of the point deformation parameter vector to obtain the phase linear change rate. Multiply the phase linear change rate by the radar wavelength and divide it by the first deformation coefficient to obtain the p1th deformation rate sequence. S63: Establish an inner loop with sequence number t2 and values ​​from 1 to the preset number of time-series interferograms. For each t2, calculate the product of the second component of the point deformation parameter vector and the t2th equally spaced time reference sequence, and then superimpose the first component of the point deformation parameter vector as the fitting phase. Multiply the fitting phase by the radar wavelength and divide it by the second deformation coefficient to obtain the cumulative deformation matrix elements under the p1th outer loop and the t2th inner loop. S64: After the inner loop with sequence number t2 and the outer loop with sequence number p1 are completed, the deformation rate sequence and the cumulative deformation matrix are obtained.

8. A high-precision deformation information extraction system for high-altitude cold canyons using GB-SAR, used to implement the high-precision deformation information extraction method for high-altitude cold canyons using GB-SAR as described in any one of claims 1-7, characterized in that, The high-precision GB-SAR deformation information extraction system for high-altitude cold canyons includes a data acquisition module, a data analysis module, and a result output module. The data acquisition module is connected to the data analysis module, and the data analysis module is connected to the result output module. The data acquisition module is used to acquire GB-SAR time-series interferometric datasets of the high-altitude cold canyon study area. The data analysis module is used to obtain high-precision surface deformation information of the high-altitude cold canyon. The result output module is used to output the high-precision surface deformation information of the high-altitude cold canyon.

Images