Freeze-thaw disaster identification method and system in frozen soil region, and storage medium

Abstract

The invention relates to the technical field of remote sensing science and technology, in particular to a frozen soil region freeze-thaw disaster identification method and system and a storage medium, and the method comprises the steps: generating a target differential interference image set; carrying out homogeneous point filtering, phase unwrapping optimization and unwrapping error correction on the basis of the target differential interference image set to obtain an unwrapped phase image set; for each unwrapped phase diagram, performing uniform sub-block division, and based on an elevation-phase empirical relationship, estimating and removing an atmospheric delay phase existing in each sub-block through a linear fitting method to obtain a phase optimization diagram set after atmospheric influence is removed; performing sequential superposition on unwrapping phases corresponding to the interferograms in the phase optimization diagram set, and performing analysis by using the relationship between the average phase change rate after superposition and the radar wavelength to obtain earth surface deformation data; and unsupervised classification and screening are carried out based on the earth surface deformation data, and freeze-thaw change areas are optimized, so that a high-precision freeze-thaw recognition result is obtained.

Description

Technical Field

[0001] The present invention relates to the technical field of remote sensing science and technology, and particularly relates to a method and system for identifying freeze-thaw disasters in frozen soil areas. Background Art

[0002] In recent years, affected by global climate change, the degradation of frozen soil has become increasingly serious, mainly manifested as the deepening of the active layer, the melting of underground ice, etc. This change not only alters the regional ecological environment but also triggers a series of freeze-thaw disasters, having a profound impact on human society and the natural environment. Among them, freeze-thaw disasters refer to the abnormal changes of frozen soil caused by natural environmental changes and human activities, and then the phenomena or processes that damage human life and property, production and life, as well as resource environment. Such disasters mainly include the periodic frost heaving and thaw settlement of frozen soil, and various geological disasters such as thermokarst slumps, ice cones, debris flows, landslides, collapses, and soil erosion caused thereby.

[0003] In terms of deformation monitoring in frozen soil areas, traditional methods mainly include leveling, GPS observation, and burying frost heave gauges, etc. Although these methods have high measurement accuracy and reliable results, they mainly focus on single-point and small-area monitoring and are difficult to meet the current research needs of large-scale and high spatio-temporal resolution. With the development of remote sensing technology, Interferometric Synthetic Aperture Radar (InSAR) technology has become an important means for monitoring frozen ground surface deformation due to its advantages of high precision, all-weather, all-day, and wide coverage. In particular, multi-temporal InSAR (MT-InSAR) methods, such as Persistent Scatterer InSAR (PSI), Small Baseline Subset InSAR (SBAS-InSAR), and SqueeSAR, etc., by jointly screening high-coherence pixels (such as permanent scatterers and distributed scatterers), significantly improve the spatial sampling density while ensuring accuracy and effectively suppress the accumulation of phase unwrapping errors in complex surface environments. However, the application of InSAR technology in frozen soil areas still faces many challenges:

[0004] First of all, the influence of atmospheric delay effects on monitoring accuracy is particularly significant. Especially in frozen soil areas with complex climate environments and variable meteorological conditions, factors such as high-frequency temperature changes, precipitation, and snow cover will exacerbate the instability of atmospheric delay, making it difficult for traditional atmospheric correction methods to effectively eliminate errors, thus affecting the accuracy of deformation measurement. Although there are already atmospheric correction methods based on meteorological models, GNSS data, and time series analysis, there are still problems such as limited applicability and difficulty in obtaining data in frozen soil areas.

[0005] Secondly, the deformation process in permafrost areas is relatively complex, not only affected by temperature changes, but also closely related to factors such as water migration, underground ice content, and surface vegetation cover. Permafrost deformation is usually slow and shows significant seasonal characteristics, such as frost heave in winter and thaw settlement in summer, which makes it difficult for traditional InSAR methods to extract signals. Especially in time series analysis, how to effectively distinguish deformation signals caused by freeze-thaw cycles from other interference signals is still a challenge. In addition, uneven settlement and local sudden deformation of the surface in permafrost areas are often related to small-scale environmental factors, and traditional InSAR technology is difficult to monitor over a large area while taking into account high-resolution deformation extraction.

[0006] Finally, most current research focuses on deformation monitoring, while specialized identification methods for freeze-thaw disasters are still relatively lacking. Surface collapse, thermokarst landform evolution, and thaw-sinking disasters caused by permafrost degradation are often difficult to identify through deformation monitoring alone in the early stages. Existing InSAR analysis methods mainly focus on deformation trends, and the automatic identification, classification, and disaster warning capabilities of abnormal deformations still need to be improved.

[0007] Therefore, how to achieve high-precision deformation monitoring and accurate identification of freeze-thaw disasters is of great significance to improving the deformation monitoring capabilities and disaster warning levels in permafrost areas. Summary of the invention

[0008] In order to achieve high-precision deformation monitoring and accurate identification of freeze-thaw disasters, the purpose of the present invention is to provide a method and system for identifying freeze-thaw disasters in frozen soil areas. The technical solutions adopted are as follows:

[0009] In a first aspect, the present application discloses a method for identifying freeze-thaw disasters in a frozen soil area, the method comprising:

[0010] S1. Based on the acquired target SAR image set and DSM data, a target differential interferogram set is generated through main image selection and registration, differential interferometry processing based on time-space baseline threshold, terrain phase correction, and complex conjugate multiplication;

[0011] S2, performing homogeneous point filtering, phase unwrapping optimization, and unwrapping error correction based on the target differential interference atlas to obtain an unwrapped phase atlas;

[0012] S3. For each unwrapped phase map, the uniform sub-blocks are divided respectively, and based on the empirical relationship between elevation and phase, the atmospheric delay phase in each sub-block is estimated and removed by a linear fitting method to obtain a phase optimization atlas after removing the atmospheric influence;

[0013] S4. Superpose the unwrapped phases corresponding to the interference maps in the phase optimization map set in time series, and analyze and obtain the surface deformation data by using the relationship between the average phase change rate after superposition and the radar wavelength;

[0014] S5. Based on the surface deformation data, perform unsupervised classification, screening, and optimization of the freeze-thaw change area to obtain a high-precision freeze-thaw identification result.

[0015] Further, in step S1, based on the acquired target SAR image set and DSM data, through master image selection and registration, differential interferometry processing based on the spatio-temporal baseline threshold, terrain phase correction, and complex conjugate multiplication, a target differential interferometry map set is generated, including:

[0016] S11. Acquire the target SAR image set covering the study area and DSM data;

[0017] S12. Select the master image from the target SAR image set, and register the other images in the set with the master image to obtain a registered SAR image set;

[0018] S13. Perform differential interferometry processing on the registered SAR image set based on a preset spatio-temporal baseline threshold to obtain an initial differential interferometry map set;

[0019] S14. According to the time baseline and spatial baseline between the images, select the initial differential interferometry pairs from the initial differential interferometry map set, and simulate the terrain phase by using the DSM data to subtract the terrain phase from the initial differential interferometry pairs to obtain the target differential interferometry pairs;

[0020] S15. For each target differential interferometry pair, perform complex conjugate multiplication respectively to obtain the target differential interferometry map set.

[0021] Further, in step S2, based on the target differential interferometry map set, perform homogeneous point filtering, phase unwrapping optimization, and unwrapping error correction to obtain an unwrapped phase map set, including:

[0022] S21. Select homogeneous points with similar phase values based on the target differential interferometry map set, and perform filtering processing based on the selected homogeneous points to obtain a coherence coefficient map set reflecting the coherence between different pixel points;

[0023] S22. Based on the coherence coefficient map set, screen high-quality pixel points through a preset amplitude deviation threshold and time coherence coefficient threshold to obtain a set of high-coherence target points;

[0024] S23. Construct a phase unwrapping network based on the high-coherence target point set, and exclude the low-coherence edges in the phase unwrapping network by using the shortest path algorithm to achieve phase unwrapping optimization;

[0025] S24. Based on the minimum cost flow method, calculate the minimum cost path of each node in the phase unwrapping network, and accumulate the phase gradient along the path to restore the true phase value of each node;

[0026] S25. For the phase unwrapping error existing in the time dimension, correct the unwrapping error based on the closed phase information of redundant observations to obtain an unwrapped phase atlas.

[0027] Further, in step S25, the correcting the unwrapping error based on the closed phase information of redundant observations for the phase unwrapping error existing in the time dimension to obtain an unwrapped phase atlas includes:

[0028] S251. In the time dimension, define the unwrapped phase triangle closed loop of a given pixel in the following form:

[0029]

[0030] where C represents the matrix containing all closed loops, represents n unwrapped phases, U represents the integer cycle correction vector to be solved, and its size is n×1;

[0031] S252. When it is determined that the solved integer cycle correction vector U is not equal to 0, perform integer linear programming based on the branch cut method to estimate the parameter U and correct the result of phase unwrapping;

[0032] The unwrapped phase of each pixel after correction includes orbital error, terrain error, surface deformation, and atmospheric delay and noise phase components. Among them, the orbital error phase component is removed based on a bilinear polynomial, the terrain error phase component is removed based on the least squares method and the linear relationship between the vertical baseline and the terrain error, and the atmospheric effect related to elevation is removed based on the linear fitting method of elevation and phase.

[0033] Further, in step S3, the estimating and removing the atmospheric delay phase existing in each sub-block by using a linear fitting method based on the elevation-phase empirical relationship includes:

[0034] S31. Based on the elevation-phase empirical relationship, use a linear regression model to fit the phase data and elevation data in each sub-block to establish a linear relationship between the observed phase and elevation;

[0035] S32. Estimate the parameters to be solved in the linear regression model based on the least squares algorithm, and calculate and remove the residual phase based on the difference between the predicted phase and the actual observed phase of the fitted model.

[0036] Further, in the stage of sub-block splicing, the method further includes:

[0037] S33. Use a cosine weighting function to fuse the phases in the overlapping area to suppress the block-to-block mutation and maintain the spatial continuity of the phase field. The formula is as follows:

[0038]

[0039] where (x, y) represents the two-dimensional coordinates of the pixel in the image, and w k (x, y) represents the fusion weight value of the k-th sub-block at the pixel coordinates (x, y), and d x , d y respectively represent the horizontal distance and the vertical distance from the pixel (x, y) to the boundary of the k-th sub-block, and d represents the preset overlapping area width.

[0040] Further, in step S4, the average phase change rate after superposition is shown by the following formula:

[0041]

[0042] where Δt j represents the time interval of the j-th interferogram, represents the unwrapped phase of the j-th interferogram, and N represents the total number of interferograms superimposed;

[0043] In step S4, the relationship between the average phase change rate after superposition and the radar wavelength λ is shown by the following formula:

[0044]

[0045] where V disp represents the surface displacement deformation rate obtained by conversion.

[0046] Further, in step S5, the method of performing unsupervised classification, screening, and optimizing the freeze-thaw change area based on the surface deformation data to obtain a high-precision freeze-thaw identification result includes:

[0047] S51. Input the data to be classified and set the clustering parameters;

[0048] S52. Perform unsupervised clustering analysis based on the ISODATA clustering analysis algorithm to obtain the corresponding classification results;

[0049] S53, filtering out relevant data corresponding to the target category based on the classification result, and performing counting statistics based on the relevant data to determine the number of pixels included in each target category;

[0050] S54. Based on the counting statistics in step S53, the non-critical areas belonging to the designated first category and having a number of pixels less than a preset threshold are eliminated, and the critical areas belonging to the designated second category and having a number of pixels greater than the preset threshold are extracted.

[0051] In the second aspect, the present application discloses a freeze-thaw disaster identification system in a frozen soil area, the system comprising a data preprocessing module, a phase unwrapping module, an atmospheric correction module, a deformation analysis module, and a freeze-thaw identification optimization module, wherein:

[0052] The data preprocessing module is used to generate a target differential interferogram set based on the acquired target SAR image set and DSM data through main image selection and registration, differential interferometry processing based on time-space baseline threshold, terrain phase correction, and complex conjugate multiplication;

[0053] The phase unwrapping module is used to perform homogeneous point filtering, phase unwrapping optimization, and unwrapping error correction based on the target differential interference atlas to obtain an unwrapped phase atlas;

[0054] The atmospheric correction module is used to divide each unwrapped phase map into uniform sub-blocks, and based on the elevation-phase empirical relationship, estimate and remove the atmospheric delay phase in each sub-block through a linear fitting method to obtain a phase optimization atlas after removing the atmospheric influence;

[0055] The deformation analysis module is used to perform temporal superposition of the unwrapped phases corresponding to the interference patterns in the phase optimization atlas, and to obtain surface deformation data by analyzing the relationship between the average phase change rate after superposition and the radar wavelength;

[0056] The freeze-thaw identification optimization module is used to perform unsupervised classification, screening, and optimization of freeze-thaw change areas based on the surface deformation data to obtain high-precision freeze-thaw identification results.

[0057] In a third aspect, the present application discloses a computer storage medium, which is used to store computer execution instructions, and the computer execution instructions are used to execute the method for identifying freeze-thaw disasters in frozen soil areas as described in any of the preceding claims.

[0058] The present invention has the following beneficial effects:

[0059] 1) Generate high-quality differential interferograms through spatiotemporal baseline threshold control, which can reduce spatiotemporal decorrelation and noise interference and improve deformation monitoring accuracy;

[0060] 2) Effectively remove phase noise based on the same - point filtering and phase - unwrapping error correction, which can ensure the continuity and accuracy of phase information;

[0061] 3) Divide the study area by constructing uniform sub - blocks based on the regular grid - based block method, and correct the atmospheric delay phase of each sub - block through the linear fitting method, thereby improving the accuracy and calculation efficiency of atmospheric correction;

[0062] 4) Based on the ISO clustering method, preliminarily classify the deformation rate data through the ISODATA clustering algorithm, and optimize the clustering results through zonal statistics and threshold screening, so as to accurately identify the high - risk areas of freeze - thaw disasters and ensure spatial continuity. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] In order to more clearly illustrate the technical solutions and advantages in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required for the description of the embodiments or the prior art. Obviously, the following - described drawings are only some embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other drawings can also be obtained based on these drawings.

[0064] Figure 1 It is a flowchart of a method for identifying freeze - thaw disasters in frozen soil areas provided by an embodiment of the present invention;

[0065] Figure 2 It is a system structure diagram of a system for identifying freeze - thaw disasters in frozen soil areas provided by an embodiment of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0066] In order to further elaborate on the technical means and effects adopted by the present invention to achieve the intended invention purpose, the following, in combination with the drawings and preferred embodiments, details the specific implementation manners, structures, features, and effects of a method and system for identifying freeze - thaw disasters in frozen soil areas proposed according to the present invention. In the following description, different "one embodiment" or "another embodiment" do not necessarily refer to the same embodiment. In addition, the specific features, structures, or characteristics in one or more embodiments can be combined in any suitable form.

[0067] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by those of ordinary skill in the technical field to which the present invention belongs.

[0068] The following specifically describes the specific solutions of a method and system for identifying freeze - thaw disasters in frozen soil areas provided by the present invention with reference to the drawings.

[0069] Please refer to Figure 1, which shows a method flow chart of a method for identifying freeze-thaw disasters in frozen soil areas provided by an embodiment of the present invention. The method includes:

[0070] Step S1, based on the acquired target SAR image set and DSM data, generate a target differential interferogram set through main image selection and registration, differential interferometry processing based on spatio-temporal baseline threshold, terrain phase correction, and complex conjugate multiplication.

[0071] Step S2, perform homogeneous point filtering, phase unwrapping optimization, and unwrapping error correction based on the target differential interferogram set to obtain an unwrapped phase map set.

[0072] Step S3, for each unwrapped phase map, divide it into uniform sub-blocks respectively, and based on the elevation-phase empirical relationship, estimate and remove the atmospheric delay phase existing in each sub-block through a linear fitting method to obtain a phase optimization map set after removing the atmospheric influence.

[0073] Step S4, perform temporal superposition on the unwrapped phases corresponding to each interferogram in the phase optimization map set, and analyze to obtain surface deformation data by using the relationship between the average phase change rate after superposition and the radar wavelength.

[0074] Step S5, perform unsupervised classification, screening, and optimization of freeze-thaw change regions based on the surface deformation data to obtain a high-precision freeze-thaw identification result.

[0075] As can be seen from the above, a method for identifying freeze-thaw disasters in frozen soil areas disclosed in the present application generates a high-quality differential interferogram set through spatio-temporal baseline threshold control, which can reduce spatio-temporal decorrelation and noise interference and improve the accuracy of deformation monitoring; effectively remove phase noise based on homogeneous point filtering and unwrapping error correction, which can ensure the continuity and accuracy of phase information; divide the study area by constructing uniform sub-blocks based on the regular grid partitioning method, and correct the atmospheric delay phase of each sub-block through a linear fitting method, thereby improving the accuracy and calculation efficiency of atmospheric correction; based on the ISO clustering method, initially classify the deformation rate data through the ISODATA clustering algorithm, and optimize the clustering results through partition statistics and threshold screening, thereby accurately identifying high-risk regions of freeze-thaw disasters and ensuring spatial continuity.

[0076] In one embodiment, in step S1, the generating a target differential interferogram set through main image selection and registration, differential interferometry processing based on spatio-temporal baseline threshold, terrain phase correction, and complex conjugate multiplication based on the acquired target SAR image set and DSM data includes:

[0077] Step S11, acquire the target SAR image set covering the study area and DSM data.

[0078] Specifically, the present application uses the Sentinel-1 image set that covers the research area, has a time period close to one year (or a specified period according to research requirements, such as half a year, two years, etc.), and has a waveband type of C-band as the target SAR image set.

[0079] Specifically, the present application uses the AW3D DSM (ALOS World 3D Digital Surface Model) data with a resolution of 30m as the external reference DEM to simulate the topographic phase of the interferogram.

[0080] In one embodiment, the present application also performs multi-looking processing in the azimuth and range directions on the target SAR image set to reduce the noise levels in the azimuth and range directions of the images, balance the spatial resolution and signal-to-noise ratio, and provide more stable input data for subsequent interferometric processing. In a specific implementation case, the present application sets the multi-looking coefficient in the azimuth direction to 4 and the multi-looking coefficient in the range direction to 1, thereby significantly improving the signal-to-noise ratio in the azimuth direction while retaining the high resolution in the range direction, meeting the requirements for data quality in subsequent interferometric processing.

[0081] Step S12: Select a master image from the target SAR image set, and register the other images in the set with the master image to obtain a registered SAR image set.

[0082] Specifically, the present application selects a master image from the target SAR image set based on benchmarks such as the time baseline, spatial baseline, and imaging quality (such as coherence, cloud coverage ratio) of the images, and uses a corresponding configuration method (the present application does not limit the specific registration method, which can be selected according to actual application requirements) to register the other images in the set with the master image to obtain a registered SAR image set.

[0083] Step S13: Perform differential interferometric processing on the registered SAR image set based on a preset spatio-temporal baseline threshold to obtain an initial differential interferogram set.

[0084] Specifically, the selection of the spatio-temporal baseline threshold is crucial. If the spatio-temporal baseline threshold is set too small, the number of generated interferometric pairs will be small, resulting in too few SAR images participating in deformation monitoring and unreliable deformation information obtained. If the spatio-temporal baseline threshold is set too large, interferometric pairs with poor spatio-temporal coherence will participate in the calculation, seriously affecting the accuracy of deformation monitoring. Therefore, in order to obtain reliable deformation information while ensuring the accuracy of deformation monitoring, the spatio-temporal baseline threshold can be set based on the deformation rate of the research area, terrain complexity, radar system parameters (such as wavelength, incident angle), and historical interferometric data quality. The present application does not limit this.

[0085] Step S14: Select initial differential interferogram pairs from the initial differential interferogram set according to the temporal baseline and spatial baseline between the images, and simulate the topographic phase using the DSM data to subtract the topographic phase from the initial differential interferogram pairs to obtain target differential interferogram pairs.

[0086] Specifically, when simulating the topographic phase using DSM data, the DSM needs to be rasterized through a high-precision interpolation algorithm (such as inverse distance weighted interpolation or Kriging interpolation) to match its spatial resolution with that of the SAR image. Then, through geocoding, the DSM is projected into the SAR image coordinate system, and the topographic phase is calculated based on the geometric relationship between the line of sight (LOS) of the radar and the terrain slope to ensure that the topographic phase is accurately removed from the initial differential interferogram pairs, thereby obtaining target differential interferogram pairs that reflect real deformation information.

[0087] It should be noted that in the process of selecting initial differential interferogram pairs from the initial differential interferogram set: for N + 1 SAR images of the same study area, this application will generate M initial differential interferogram pairs based on the temporal baseline and spatial baseline between the images. Among them, assuming N is odd, the number of generated differential interferogram pairs M should satisfy: This ensures the effectiveness and quantity of the differential interferogram pairs to maximize the utilization of SAR image data within a reasonable range.

[0088] Step S15: For each target differential interferogram pair, perform complex conjugate multiplication respectively to obtain a target differential interferogram set.

[0089] Specifically, for each target differential interferogram pair, it is necessary to ensure that the two SAR images have been accurately registered, and then perform pixel-by-pixel complex conjugate multiplication based on the complex data of the two SAR images to finally generate the corresponding target differential interferogram set.

[0090] In the above embodiment, first, an initial interferogram set is generated through SAR image registration and differential interferometry processing, then high-quality initial interferogram pairs are selected by combining spatio-temporal baseline thresholds, and the influence of the topographic phase is removed by DSM data. Finally, a target differential interferogram set that only reflects surface deformation is generated through complex conjugate multiplication, thereby effectively improving the accuracy and reliability of deformation monitoring.

[0091] In one of the embodiments, in step S2, the process of performing homogeneous point filtering, phase unwrapping optimization, and unwrapping error correction on the target differential interferogram set to obtain an unwrapped phase map set includes:

[0092] Step S21: Select homogeneous points with similar phase values based on the target differential interferogram set, and perform filtering processing based on the selected homogeneous points to obtain a coherence coefficient map set that reflects the coherence between different pixel points.

[0093] Specifically, for each determined target differential interferogram, the present application will select co - homogenous points with similar phase values therefrom based on a fast co - homogenous point selection method, and perform co - homogenous point filtering analysis on the distributed pixels therein to improve the stability of the distributed scatterer phase, and generate a corresponding coherence coefficient map.

[0094] Step S22: Based on the coherence coefficient map set, high - quality pixel points are screened through a preset amplitude deviation threshold and a temporal coherence coefficient threshold to obtain a corresponding high - coherence target point set.

[0095] Specifically, the amplitude deviation threshold and the temporal coherence coefficient threshold can be set according to the specific surface characteristics of the research area, the quality of image data, and subsequent application requirements. The present application does not limit this. In one embodiment, based on on - site investigation of the research area, historical data statistics, and the specific accuracy requirements of the current monitoring task, the amplitude deviation threshold can be set to 0.6, and the temporal coherence coefficient threshold can be set to 0.3.

[0096] Step S23: Based on the high - coherence target point set, a phase unwrapping network is constructed, and the low - coherence edges in the phase unwrapping network are excluded by using the shortest - path algorithm to achieve phase unwrapping optimization.

[0097] Specifically, for the selected pixel points, the present application will construct a phase unwrapping network based on a triangulated irregular network. By taking each point in the high - coherence target point set as a node and the phase gradient between adjacent nodes as the edge weight, a weighted network is constructed. Then, the shortest - path algorithm is used to exclude the low - coherence edges in the phase unwrapping network to optimize the network structure.

[0098] Step S24: Based on the minimum - cost flow method, the minimum - cost path of each node in the phase unwrapping network is calculated, and the phase gradient is accumulated along the path to restore the true phase value of each node.

[0099] Step S25: For the phase unwrapping error existing in the time dimension, based on the closed - phase information of redundant observations, the unwrapping error is corrected to obtain an unwrapped phase map set.

[0100] In the above - mentioned embodiments, by selecting co - homogenous points and performing filtering processing to obtain a coherence coefficient map set reflecting the coherence of pixel points, then using preset thresholds to screen out a high - coherence target point set, constructing and optimizing a phase unwrapping network, and finally restoring the true phase value based on the minimum - cost flow method and correcting the phase unwrapping error in the time dimension, the accuracy and reliability of phase unwrapping are effectively improved, a high - quality unwrapped phase map set is obtained, and more accurate data support is provided for subsequent surface deformation monitoring.

[0101] In one embodiment, in step S25, for the phase unwrapping error existing in the time dimension, the unwrapping error is corrected based on the closed phase information of redundant observations to obtain an unwrapped phase atlas, including:

[0102] Step S251, in the time dimension, define the unwrapped phase triangle closed loop of a given pixel in the following form:

[0103]

[0104] where C represents the matrix containing all closed loops, represents n unwrapped phases, U represents the integer cycle correction vector to be solved, and its size is n×1.

[0105] Step S252, when it is determined that the solved integer cycle correction vector U is not equal to 0, perform integer linear programming based on the branch cut method to estimate the parameter U to correct the result of phase unwrapping.

[0106] Specifically, the branch cut method divides the phase discontinuous region by constructing branch cut lines in the phase unwrapping network, and uses the integer linear programming optimization algorithm to solve the integer cycle correction parameter under the condition of satisfying the phase continuity constraint, thereby correcting the 2π jump error generated in the phase unwrapping process.

[0107] After correction, the unwrapped phase of each pixel includes orbit error, terrain error, surface deformation, and atmospheric delay and noise phase components. Among them, the orbit error phase component is removed based on a bilinear polynomial, the terrain error phase component is removed based on the least squares method and the linear relationship between the vertical baseline and the terrain error, and the atmospheric effect related to elevation is removed based on the linear fitting method of elevation and phase.

[0108] In one embodiment, in step S3, based on the elevation-phase empirical relationship, estimate and remove the atmospheric delay phase existing in each sub-block by a linear fitting method, including:

[0109] Step S31, based on the elevation-phase empirical relationship, use a linear regression model to fit the phase data and elevation data in each sub-block, and establish a linear relationship between the observed phase and elevation.

[0110] It should be noted that this application will generate a grid space index matrix based on memory constraints, and evenly divide the entire image into M×N sub-blocks based on this index matrix. Among them, the size L x ×L y of each sub-block can be dynamically adjusted according to memory. In one embodiment, the size of each sub-block satisfies the following formula:

[0111] L x×L y ≤Available memory / Memory occupancy per single pixel · η;

[0112] where η is a safety factor with a value ranging from 0.6 to 0.8.

[0113] In addition, to avoid discontinuous phase unwrapping at the sub - block edges or the propagation of atmospheric delay correction errors, an overlapping area d is set between adjacent sub - blocks, and its width is jointly determined by the radar wavelength λ, the slant range R, the vertical baseline B ⊥ , and the digital elevation error Δh.

[0114] In one embodiment, the value range of the overlapping area d satisfies the following inequality:

[0115]

[0116] where α is the elevation sensitivity coefficient.

[0117] Specifically, the linear relationship between the observed phase and the elevation h can be referred to the following formula:

[0118]

[0119] where α is the elevation sensitivity coefficient, β is the constant term, and ε is the fitting residual.

[0120] Step S32: Estimate the parameters to be solved in the linear regression model based on the least - squares algorithm, and calculate and remove the residual phase based on the difference between the predicted phase of the fitted model and the actual observed phase.

[0121] Specifically, in this application, α and β are solved by minimizing the sum of squared errors between the actual observed phase and the model - predicted phase:

[0122]

[0123] where N represents the total number of sub - blocks, represents the observed phase of the i - th sub - block, and h i represents the elevation of the i - th sub - block.

[0124] After solving the parameters to be solved α and β, this application calculates the residual phase of each sub - block based on the following formula

[0125] In one embodiment, during the sub - block splicing stage, the method further includes:

[0126] Step S33: Use a cosine - weighted function to fuse the phase in the overlapping area to suppress the block - to - block mutation and maintain the spatial continuity of the phase field. The formula is as follows:

[0127]

[0128] Among them, (x, y) represents the two-dimensional coordinates of the pixel in the image, w k (x, y) represents the fusion weight value of the k-th sub-block at the pixel coordinates (x, y), d x , d y respectively represent the horizontal distance and the vertical distance from the pixel (x, y) to the boundary of the k-th sub-block, and d represents the preset width of the overlapping area.

[0129] Specifically, in the sub-block splicing stage, the present application uses a cosine weighting function to fuse the phases in the overlapping area. The above function specifically calculates the fusion weight value according to the distance (horizontal and vertical distances) from the pixel to the sub-block boundary to ensure the spatial continuity of the phase field and suppress the inter-block mutation.

[0130] In one embodiment, in step S4, the average phase change rate after superposition is shown by the following formula:

[0131]

[0132] Among them, Δt j represents the time interval of the j-th interferogram, represents the unwrapped phase of the j-th interferogram, and N represents the total number of the superimposed interferograms.

[0133] In step S4, the relationship between the average phase change rate after superposition and the radar wavelength λ is shown by the following formula:

[0134]

[0135] Among them, V disp represents the obtained ground displacement deformation rate.

[0136] It should be noted that the unwrapped phases corresponding to the interferograms in the phase optimization atlas are superimposed in time series. This process accumulates the phase information to form a time series phase superposition diagram, thereby improving the signal-to-noise ratio between the deformation information and the atmospheric error term. The atmospheric error phase after superposition is not the result of a simple multiple growth of the atmospheric phase error in each interferogram, but the result of a square root multiple growth of the number of interferograms, which improves the signal-to-noise ratio between the deformation information and the atmospheric error term.

[0137] In one embodiment, in step S5, the unsupervised classification, screening, and optimization of the freeze-thaw change area based on the surface deformation data to obtain a high-precision freeze-thaw identification result include:

[0138] Step S51: Input the data to be classified and set the clustering parameters.

[0139] Specifically, in this application, the maximum number of categories is set to 12, the minimum number of category samples is 100, and the sampling interval is 20, so as to control the complexity of the clustering process and the rationality of the category distribution, and avoid overfitting or underfitting phenomena.

[0140] Step S52: Perform unsupervised clustering analysis based on the ISODATA clustering analysis algorithm to obtain the corresponding classification results.

[0141] Specifically, in each iteration, the statistical features (such as mean, variance, etc.) of each category will be calculated according to the current clustering results, and the number of categories will be dynamically adjusted according to the preset merging and splitting rules until the convergence condition is met or the maximum number of iterations is reached.

[0142] Step S53: Screen out the relevant data corresponding to the target categories based on the classification results, and perform a count statistics based on the relevant data to determine the number of pixels included in each target category.

[0143] Specifically, assuming that the category numbers of the target categories are 0, 1, 2, 10, and 11, this application uses the following SQL query statement to remove the irrelevant data with category numbers from 3 to 9, and in the case of only retaining the target categories, the irrelevant areas can be removed: SELECT * FROM classification_results WHERE category NOT IN (3, 4, 5, 6, 7, 8, 9).

[0144] Specifically, to accurately quantify the pixel scale covered by each target category, this application uses the following SQL query statement to group the data of different target categories and count the number of pixels included in each target category: SELECT category, COUNT(*) AS pixel_count FROM classification_results GROUP BY category.

[0145] In one embodiment, to remove small regions that are scattered and have weak statistical significance, the present application also excludes categories with a pixel count less than 9000. The following SQL statement is specifically used to ensure the spatial continuity and rationality of the classification results while only retaining regions with sufficient representativeness: DELETE FROM classification_results WHERE category IN (SELECT category FROM classification_results GROUP BY category HAVING COUNT(*) < 9000).

[0146] Step S54, based on the count statistics in step S53, exclude non-critical regions belonging to the specified first category and having a pixel count less than the preset threshold, and extract critical regions belonging to the specified second category and having a pixel count greater than the preset threshold.

[0147] Specifically, to further optimize the screening criteria and ensure that no critical deformation information is missed, the present application further sets screening conditions to exclude non-critical regions with a category number between 3 and 10 (i.e., the category number of the first category) and a pixel count less than 5000, so as to ensure the integrity of valid information. The following SQL statement is specifically used to refine the classification results and remove small regions that may cause misjudgment: DELETE FROM classification_results WHERE category BETWEEN 3 AND 10 AND pixel_count <= 5000.

[0148] Furthermore, to ensure the complete retention of critical categories and ensure that important deformation regions are not accidentally deleted, the present application uses the following SQL statement to extract critical regions with a category number of 1, 2, 11, 12 (i.e., the category number of the second category) and all pixel counts greater than 5000: SELECT * FROM classification_results WHERE category IN (1, 2, 11, 12) OR pixel_count > 5000.

[0149] Please refer to Figure 2 , a freeze-thaw disaster identification system disclosed in the present application, the system includes a data preprocessing module, a phase unwrapping module, an atmospheric correction module, a deformation analysis module, and a freeze-thaw identification optimization module, wherein:

[0150] The data preprocessing module is used to generate a target differential interferogram set based on the acquired target SAR image set and DSM data through main image selection and registration, differential interferometry processing based on spatio-temporal baseline threshold, terrain phase correction, and complex conjugate multiplication.

[0151] The phase unwrapping module is used to perform homogeneous point filtering, phase unwrapping optimization, and unwrapping error correction based on the target differential interferogram set to obtain an unwrapped phase map set.

[0152] The atmospheric correction module is used to divide each unwrapped phase map into uniform sub-blocks, and estimate and remove the atmospheric delay phase existing in each sub-block through a linear fitting method based on the elevation-phase empirical relationship to obtain a phase optimization map set after removing the influence of the atmosphere.

[0153] The deformation analysis module is used to perform temporal superposition on the unwrapped phases corresponding to each interferogram in the phase optimization map set, and analyze and obtain surface deformation data by using the relationship between the average phase change rate after superposition and the radar wavelength.

[0154] The freeze-thaw identification optimization module is used to perform unsupervised classification, screening, and optimization of the freeze-thaw change area based on the surface deformation data to obtain a high-precision freeze-thaw identification result.

[0155] In one embodiment, the above modules are also used to implement a method for identifying freeze-thaw disasters in frozen soil areas as described in any one of the foregoing method embodiments, and this application does not make any limitations in this regard.

[0156] As can be seen from the above, a system for identifying freeze-thaw disasters in frozen soil areas disclosed in this application generates a high-quality differential interferogram set through spatio-temporal baseline threshold control, which can reduce spatio-temporal decorrelation and noise interference and improve the accuracy of deformation monitoring; effectively removes phase noise based on homogeneous point filtering and unwrapping error correction, which can ensure the continuity and accuracy of phase information; divides the study area by constructing uniform sub-blocks based on the regular grid partitioning method, and corrects the atmospheric delay phase of each sub-block through a linear fitting method, thereby improving the accuracy and calculation efficiency of atmospheric correction; based on the ISO clustering method, preliminarily classifies the deformation rate data through the ISODATA clustering algorithm, and optimizes the clustering results through zonal statistics and threshold screening, thereby accurately identifying high-risk areas of freeze-thaw disasters and ensuring spatial continuity.

[0157] Finally, this application also discloses a computer storage medium, which is used to store computer execution instructions, and the computer execution instructions are used to execute the method for identifying freeze-thaw disasters in frozen soil areas described in any one of the foregoing claims.

[0158] As can be seen from the above, a computer storage medium disclosed in this application can generate a high-quality differential interference atlas through spatio-temporal baseline threshold control, which can reduce spatio-temporal decorrelation and noise interference and improve the accuracy of deformation monitoring. Based on the same-point filtering and phase unwrapping error correction, phase noise can be effectively removed, ensuring the continuity and accuracy of phase information. Based on the regular grid partitioning method, the study area is divided by constructing uniform sub-blocks, and the atmospheric delay phase of each sub-block is corrected by the linear fitting method, thereby improving the accuracy and calculation efficiency of atmospheric correction. Based on the ISO clustering method, the deformation rate data is initially classified by the ISODATA clustering algorithm, and the clustering results are optimized through zonal statistics and threshold screening, so as to accurately identify the high-risk areas of freeze-thaw disasters and ensure spatial continuity.

[0159] It should be noted that the above sequence of the embodiments of the present invention is only for description and does not represent the superiority or inferiority of the embodiments. The processes depicted in the accompanying drawings do not necessarily require the specific order or continuous order shown to achieve the desired results. In some embodiments, multitasking and parallel processing are also possible or may be advantageous.

[0160] Each embodiment in this specification is described in a progressive manner, and the same or similar parts between the embodiments can be referred to each other. Each embodiment focuses on the differences from other embodiments.

Claims

1. A method for identifying freeze-thaw disasters in frozen soil areas, characterized in that The method includes: S1. Based on the acquired target SAR image set and DSM data, through master image selection and registration, differential interferometric processing based on spatio-temporal baseline threshold, topographic phase correction, and complex conjugate multiplication, generate a target differential interferogram set; S2. Based on the target differential interferogram set, perform homogeneous point filtering, phase unwrapping optimization, and unwrapping error correction to obtain an unwrapped phase map set; S3. For each unwrapped phase map, divide it into uniform sub-blocks respectively, and based on the elevation-phase empirical relationship, estimate and remove the atmospheric delay phase existing in each sub-block through linear fitting method to obtain a phase optimization map set after removing the atmospheric influence; S4. Temporally superimpose the unwrapped phases corresponding to each interferogram in the phase optimization map set, and analyze to obtain surface deformation data by using the relationship between the average phase change rate after superposition and the radar wavelength; S5. Based on the surface deformation data, perform unsupervised classification, screening, and optimization of freeze-thaw change areas to obtain a high-precision freeze-thaw identification result.

2. The method according to claim 1, characterized in that, In step S1, the generating of the target differential interferogram set based on the acquired target SAR image set and DSM data, through master image selection and registration, differential interferometric processing based on spatio-temporal baseline threshold, topographic phase correction, and complex conjugate multiplication, includes: S11. Acquire the target SAR image set covering the study area and DSM data; S12. Select a master image from the target SAR image set, and register other images in the set with the master image to obtain a registered SAR image set; S13. Perform differential interferometric processing on the registered SAR image set based on a preset spatio-temporal baseline threshold to obtain an initial differential interferogram set; S14. According to the temporal baseline and spatial baseline between images, select initial differential interferogram pairs from the initial differential interferogram set, and simulate the topographic phase by using the DSM data to subtract the topographic phase from the initial differential interferogram pairs to obtain target differential interferogram pairs; S15. For each target differential interferogram pair, perform complex conjugate multiplication respectively to obtain a target differential interferogram set.

3. The method according to claim 1, characterized in that, In step S2, the obtaining of the unwrapped phase map set by performing homogeneous point filtering, phase unwrapping optimization, and unwrapping error correction based on the target differential interferogram set includes: S21. Based on the target differential interferogram set, select homogeneous points with similar phase values, and perform filtering processing based on the selected homogeneous points to obtain a coherence coefficient map set reflecting the coherence between different pixel points; S22. Based on the coherence coefficient map set, screen high-quality pixel points through a preset amplitude deviation threshold and temporal coherence coefficient threshold to obtain a set of high-coherence target points; S23. Based on the set of high-coherence target points, construct a phase unwrapping network, and exclude low-coherence edges in the phase unwrapping network by using the shortest path algorithm to achieve phase unwrapping optimization; S24. Based on the minimum cost flow method, calculate the minimum cost path of each node in the phase unwrapping network, and accumulate the phase gradient along the path to restore the true phase value of each node; S25. For the phase unwrapping error existing in the time dimension, the unwrapping error is corrected based on the closed phase information of redundant observations to obtain an unwrapping phase atlas.

4. The method according to claim 3, wherein In step S25, the phase unwrapping error existing in the time dimension is corrected based on the closed phase information of the redundant observation to obtain an unwrapped phase atlas, including: S251. In the time dimension, define the unwrapped phase triangle closed loop of a given pixel in the following form: where C represents the matrix containing all closed loops, represents n unwrapped phases, U represents the integer ambiguity correction vector to be solved, and its size is n×1; S252, when it is determined that the integer correction vector U obtained by solving is not equal to 0, integer linear programming is performed based on the branch-cut method to estimate the parameter U to correct the result of the phase unwrapping; S253. The unwrapped phase of each pixel after correction includes orbit error, terrain error, surface deformation, atmospheric delay and noise phase components, wherein the orbit error phase component is removed based on a bilinear polynomial, the terrain error phase component is removed based on the least squares method and the linear relationship between the vertical baseline and the terrain error, and the atmospheric effect related to the elevation is removed based on a linear fitting method of elevation and phase.

5. The method according to claim 1, characterized in that, In step S3, the method of estimating and removing the atmospheric delay phase in each sub-block by a linear fitting method based on the elevation-phase empirical relationship includes: S31, based on the elevation-phase empirical relationship, a linear regression model is used to fit the phase data and elevation data in each sub-block to establish a linear relationship between the observed phase and elevation; S32. Estimate the parameters to be solved in the linear regression model based on the least squares algorithm, and calculate and remove the residual phase based on the difference between the predicted phase of the fitted model and the actual observed phase.

6. The method according to claim 5, wherein At the stage of sub-block splicing, the method further comprises: S33. The cosine weighting function is used to fuse the phases in the overlapping area to suppress the mutation between blocks and maintain the spatial continuity of the phase field. The formula is as follows: Among them, (x, y) represents the two-dimensional coordinates of the pixel in the image, w k (x, y) represents the fusion weight value of the k-th sub-block at the pixel coordinates (x, y), d x ,d y respectively represent the horizontal distance and the vertical distance from the pixel (x, y) to the boundary of the k-th sub-block, and d represents the preset width of the overlapping area.

7. The method according to claim 1, wherein In step S4, the average phase change rate after superposition is shown by the following formula: where, Δt j represents the time interval of the j-th interferogram, represents the unwrapped phase of the j-th interferogram, and N represents the total number of the superimposed interferograms; In step S4, the relationship between the average phase change rate after superposition and the radar wavelength λ is shown by the following formula: Among them, V disp represents the deformation rate of the ground displacement obtained by transformation.

8. The method according to claim 1, characterized in that, In step S5, unsupervised classification, screening, and optimization of freeze-thaw change areas are performed based on the surface deformation data to obtain high-precision freeze-thaw identification results, including: S51, input the data to be classified and set the clustering parameters; S52, performing unsupervised cluster analysis based on ISODATA cluster analysis algorithm to obtain corresponding classification results; S53, filtering out relevant data corresponding to the target category based on the classification result, and performing counting statistics based on the relevant data to determine the number of pixels included in each target category; S54. Based on the counting statistics in step S53, the non-critical areas belonging to the designated first category and having a number of pixels less than a preset threshold are eliminated, and the critical areas belonging to the designated second category and having a number of pixels greater than the preset threshold are extracted.

9. A freeze-thaw disaster identification system in a permafrost area, characterized in that, The system includes a data preprocessing module, a phase unwrapping module, an atmospheric correction module, a deformation analysis module, and a freeze-thaw identification optimization module, wherein: The data preprocessing module is used to generate a target differential interferogram set by performing main image selection and registration, differential interferometry processing based on a spatio-temporal baseline threshold, topographic phase correction, and complex conjugate multiplication on the acquired target SAR image set and DSM data; The phase unwrapping module is used to perform homogeneous point filtering, phase unwrapping optimization, and unwrapping error correction on the target differential interferogram set to obtain an unwrapped phase map set; The atmospheric correction module is used to divide each unwrapped phase map into uniform sub-blocks, and estimate and remove the atmospheric delay phase existing in each sub-block through a linear fitting method based on the elevation-phase empirical relationship to obtain an optimized phase map set after removing the atmospheric influence; The deformation analysis module is used to temporally stack the unwrapped phases corresponding to each interferogram in the optimized phase map set, and analyze and obtain surface deformation data by using the relationship between the average phase change rate after stacking and the radar wavelength; The freeze-thaw identification optimization module is used to perform unsupervised classification, screening, and optimization of the freeze-thaw change area based on the surface deformation data to obtain a high-precision freeze-thaw identification result.

10. A computer storage medium, characterized in that, The computer storage medium is used to store computer execution instructions, and the computer execution instructions are used to execute the freeze-thaw disaster identification method for frozen soil areas described in any one of claims 1 to 8.

Images