Terrain segmentation based regional adaptive multi-scale insar atmospheric delay correction method

Abstract

The application discloses a regional adaptive multi-scale InSAR atmospheric delay correction method based on terrain segmentation, first, the main watershed is extracted by using DEM data and an improved terrain segmentation algorithm, the research region is adaptively divided, and multiple sub-regions with a spatial scale more suitable for the atmospheric delay distribution characteristics are obtained, and on this basis, the atmospheric delay phase in each sub-region is estimated by using a linear model iteration, so that multi-scale atmospheric delay correction is realized, and accurate estimation of local atmospheric delay is helpful.

Description

Technical Field

[0001] This invention belongs to the field of atmospheric delay phase correction technology of differential interferograms in InSAR deformation measurement technology. It relates to a regional adaptive multi-scale InSAR atmospheric delay correction method based on terrain segmentation. Specifically, it is a method that uses an improved terrain segmentation algorithm to adaptively divide the study area based on DEM elevation data, taking into account the temporal and spatial characteristics of the correlation and differences between atmospheric delay distribution and regional terrain. On this basis, a linear model is used to estimate the atmospheric delay phase at different times and scales, thereby realizing InSAR atmospheric delay phase correction. Background Technology

[0002] Tropospheric atmospheric delay has always been a major source of error in InSAR deformation monitoring results, mainly related to the differences in atmospheric water vapor content and distribution. Zebker HA et al. (Zebker HA, Rosen PA, Hensley S. Atmospheric effects in interferometric synthetic aperture radar surface deformation and topographic maps. Journal of Geophysical Research, 1997, 102(B4)) showed that even a 20% change in atmospheric water vapor content between two SAR images in a single interferometric processing can lead to a deformation error of 10–14 cm, severely impacting the accuracy and reliability of InSAR deformation monitoring. In mountainous areas with undulating terrain, the variation of water vapor with surface elevation due to stratification effects is very significant, severely affecting the accuracy of InSAR deformation monitoring. Therefore, effective atmospheric delay correction is crucial for improving the measurement accuracy of InSAR technology. Empirical function models (such as linear models) are currently one of the widely used atmospheric correction methods. By simulating the relationship between terrain elevation and phase, they can effectively estimate the atmospheric phase delay caused by vertical stratification effects at small spatial scales. However, when differences in atmospheric pressure, temperature, relative humidity, and other factors lead to regional variations in atmospheric delay on a large spatial scale, a single empirical function relationship will not be applicable to the needs of global atmospheric phase delay estimation. To address this issue, Bekaert D et al. (Bekaert D, Hooper A, and Wright TA, Spatially Variable Power Law Tropospheric Correction Technique for InSAR Data. Journal of Geophysical Research: Solid Earth, 2015, 120(2): 1345-56) proposed an atmospheric delay estimation method based on local windows. This method divides the study area into multiple small windows to accommodate the spatial scale variations in atmospheric delay characteristics.

[0003] However, existing region segmentation methods often use regular grids of a specified size or divide the space at equal scales based on the geometric extent of the study area. While these methods can mitigate the spatial variation of atmospheric delay to some extent, they fail to consider the differences in the correlation between the spatial distribution characteristics of atmospheric delay and topography. In practical applications, climate conditions in complex mountainous environments often vary with topography and watershed characteristics, leading to spatiotemporal heterogeneity in atmospheric delay phase. Conventional region segmentation methods easily introduce multiple atmospheric delays with different scales within a single window, causing biases in model parameter estimation and introducing errors during the correction process. Therefore, how to consider the correlation and differences between atmospheric delay distribution characteristics and topography, adopt reasonable methods to divide the study area, and achieve accurate estimation of multi-scale atmospheric phase are pressing challenges to be solved when using functional model methods to estimate InSAR atmospheric delay phase. Summary of the Invention

[0004] Technical problem solved: To address the above-mentioned technical problems, this invention provides a regional adaptive multi-scale InSAR atmospheric delay correction method based on terrain segmentation, which can effectively solve the multi-scale problems caused by the lack of adaptability of the window to terrain features and the spatial heterogeneity of atmospheric delay phase in traditional atmospheric correction methods.

[0005] Technical solution: A regional adaptive multi-scale InSAR atmospheric delay correction method based on terrain segmentation, comprising the following steps:

[0006] S1. Data acquisition and preprocessing: Acquire SAR image data and DEM data of the study area and perform preprocessing to obtain unwrapped phase data and DEM data in radar coordinate system;

[0007] S2. Study Area Division: Using DEM data in radar coordinate system, the main watershed within the study area is extracted based on terrain segmentation algorithm to adaptively divide the study area.

[0008] S3. Iterative estimation of model parameters: Select observation points in the image according to the set step size, remove observation points in the deformed area using the deformation rate mask, and use stable observation points to participate in the model parameter estimation.

[0009] S4. Atmospheric Delay Phase Smoothing: The estimated model parameter values ​​corresponding to each sub-region are assigned to the centroid pixel of that region as known points for interpolation. The atmospheric model parameters are fitted and interpolated to all pixels in the entire image using the B-spline function to ensure a smooth transition of atmospheric delay phase between different windows.

[0010] S5. Atmospheric Delay Correction and Deformation Result Calculation: Subtract the atmospheric delay phase processed in step S4 from the initial interferometric data, and use the atmospherically corrected interferometric data to perform InSAR deformation calculation, finally obtaining the atmospherically corrected InSAR deformation result.

[0011] Preferably, the specific process of step S2 is as follows: First, the morphological gradient of the DEM data is calculated. Then, morphological opening and closing reconstruction operations are used to perform multi-level filtering on the gradient image at different structural scales. The morphological filtering results are then weighted and averaged to obtain the gradient result after multi-level filtering. The morphological gradient of the original DEM is filtered and reconstructed to optimize the gradient image. The expression is shown in the following formula:

[0012]

[0013] Where g(x,y)' is the final gradient result, and · are the opening and closing reconstruction operators, respectively, i is the order in which the opening and closing reconstruction operations are performed, and the structure element b i The radius r increases as i increases, and r i =10i, where n is the total number of reconstruction operations, w i Let w be the weight. i =exp(1-i), where g(x,y) is the morphological gradient image.

[0014] Furthermore, the expression for the morphological gradient image g(x,y) is as follows:

[0015]

[0016] Where f(x,y) is the original image, and b is the disk-shaped structuring element. and These represent the dilation and erosion operations in grayscale morphology, respectively.

[0017] Preferably, step S3 is as follows: using the unwrapped phase data as a reference, the atmospheric delay phase is fitted with a linear model in each of the divided sub-regions, and the model parameters of each sub-region are estimated using the least squares method, thereby obtaining the multi-scale atmospheric delay phase in different regions. M-1 interferometric pairs were generated using N SLC images. The topographic and orbital phase components in each differential interferogram have been removed. n observation points are sampled from a given sub-region. The relationship between the atmospheric delay phase and surface elevation in the linear model is shown in the following equation:

[0018]

[0019] in, Let k be the atmospheric delay phase at a pixel within the i-th sub-region of the N-th SLC image. i,n and c i,n These are the slope and intercept parameters in the linear model, h, respectively. i The surface elevation corresponding to this pixel; the interferometric phase after atmospheric phase correction can be expressed as:

[0020]

[0021] in, Let i be the residual phase matrix of each observation point in the i-th sub-region after atmospheric correction. E is the Kronecker product operator. i Let n×1 be an all-1 vector. For each interference phase matrix, S represents the observed interferometric phase value at each sampling point in each interferogram. i Let A be the observation matrix used for model parameter estimation within this sub-region, and let A be the M×(N-1) design matrix describing the interferometric pair mesh. Without loss of generality, the atmospheric delay of the main image is set to 0, thereby removing its corresponding column from the design matrix. i The model parameter matrix for each SLC scene is expressed as follows:

[0022]

[0023] The least squares method is used to solve for the model parameters in each divided region of each SLC scene, thereby calculating the atmospheric phase in each scene's differential interferogram. After obtaining the preliminary estimate of the atmospheric delay phase in the interferogram, it is subtracted from the original phase unwrapping data, and the residual phase is used as new input data to iteratively estimate the model parameters. This process of estimating model parameters is repeated until the iteration stops, and the final model parameters in that region are the sum of the model parameters estimated in all iterations. That is, it satisfies:

[0024]

[0025]

[0026] in, and This represents the final result of the iterative estimation of the linear model parameters in the i-th sub-region of the M-th interferogram, where j is the iteration number. and These are the linear model parameter values ​​obtained during the j-th iteration. and By substituting the relationship between atmospheric delay phase and surface elevation in the linear model, the atmospheric delay phase of the sub-region in each SLC scene can be obtained.

[0027] Furthermore, the condition for stopping the iteration is: the change in the root mean square value of the phase residual in the sub-region during the two iterations is less than 1 mm or the maximum number of iterations is reached.

[0028] Beneficial effects: The regional adaptive multi-scale InSAR atmospheric delay correction method based on terrain segmentation of this invention first uses DEM data and an improved terrain segmentation algorithm to extract the main watershed, and adaptively divides the study area to obtain multiple sub-regions with spatial scales that better match the atmospheric delay distribution characteristics. On this basis, a linear model is used to iteratively estimate the atmospheric delay phase in each sub-region, thereby realizing multi-scale atmospheric delay correction, which helps to accurately estimate local atmospheric delay. Attached Figure Description

[0029] Figure 1 This is a flowchart of the regional adaptive multi-scale InSAR atmospheric delay correction method based on terrain segmentation according to the present invention;

[0030] Figure 2 The image shows the annual average deformation rate of the study area obtained under the following conditions: without atmospheric correction, with atmospheric correction using a linear model method, and with the terrain segmentation-based regional adaptive multi-scale InSAR atmospheric delay correction method of this invention. Detailed Implementation

[0031] The present invention will now be described in detail with reference to the accompanying drawings and specific embodiments:

[0032] Example 1

[0033] like Figure 1 As shown, the regional adaptive multi-scale InSAR atmospheric delay correction method based on terrain segmentation includes the following steps:

[0034] S1. Data Acquisition and Preprocessing:

[0035] SAR image data and DEM data of the study area are acquired. Preprocessing steps such as fine orbit correction, image registration, differential interferometry, and phase unwrapping are performed using professional software such as GAMMA and ISCE2 to obtain unwrapped phase data and DEM data in radar coordinate system. Deformation calculation is performed using software such as StaMPS to obtain the uncorrected surface deformation rate of the study area, and a rate threshold is set to obtain the deformation zone mask file.

[0036] S2. Study Area Division:

[0037] Reasonable division of the study area is a prerequisite for accurate estimation of model parameters. Using DEM data in radar coordinate system, the study area is adaptively divided based on the terrain segmentation algorithm to extract the main watersheds within the study area. The terrain segmentation algorithm uses elevation maxima as peaks and minima as valleys, extracting the corresponding embankments at the intersection of two local minima to serve as the boundaries between the two areas. The set of these boundaries is the ridgeline of the area, also known as the watershed, which is one of the commonly used algorithms in image segmentation, also called the watershed algorithm. However, in practical applications, directly using the standard watershed transform often results in oversegmentation due to too many local extrema. Therefore, this invention improves the watershed algorithm by first calculating the morphological gradient of the DEM data, then using morphological opening and closing reconstruction operations at different structural scales to perform multi-level filtering on the gradient image, and then weighting the morphological filtering results to obtain the multi-level filtered gradient result. The morphological gradient of the original DEM is then filtered and reconstructed to optimize the gradient image. Morphological gradient is defined as the difference between the dilation and erosion results of the original image. It can describe the drastic degree of elevation change in the DEM and preserve the edge contour of the target. Its calculation method is as follows:

[0038]

[0039] Where g(x,y) is the morphological gradient image, f(x,y) is the original image, and b is the disk-shaped structuring element. and These represent the dilation and erosion operations in grayscale morphology, respectively. At this point, the morphological gradient still contains many local extrema details and noise, making it impossible to directly apply watershed transform. Therefore, morphological opening and closing reconstruction operations are needed to reconstruct the gradient image to eliminate local extrema caused by irregular elevations and preserve contour extrema. Morphological reconstruction is used to restore the boundaries preserved after the opening and closing operations. By combining opening and closing reconstruction operations, both maxima and minima within the region can be corrected simultaneously. However, due to differences in texture and noise scales, using only a single specified structuring element is insufficient to achieve similar filtering effects in different scenes. Therefore, this invention utilizes structuring elements of different radii to perform opening and closing reconstructions separately, and then performs a weighted average of the results to obtain better gradient filtering results. That is:

[0040]

[0041] Where g(x,y)' is the final gradient result, ο and · are the opening and closing reconstruction operators respectively, i is the order of the opening and closing reconstruction operations, and the structuring element b i The radius r increases as i increases, and r i =10i, where n is the total number of reconstruction operations, w i Let w be the weight.i =exp(1-i), where the weight decreases as the radius increases. This is to prevent excessively large structuring elements from causing terrain contour shifts. Subsequently, the corrected gradient image can be used for watershed transformation to identify the main boundaries in the target region for adaptive region division.

[0042] S3. Iterative estimation of model parameters:

[0043] After dividing the study area, a step size (10 or 20) is set to obtain evenly distributed observation points. Then, a deformation region mask file is used to cover deformation regions with an annual average deformation rate exceeding 30 mm / y, ensuring that atmospheric delay estimation uses only stable points and preventing the introduction of errors. Using unwrapped phase data as a reference, a linear model is applied to fit the atmospheric delay phase in each sub-region after division. The model parameters for each region are estimated using the least squares method, thus obtaining the multi-scale atmospheric delay phase for different regions. The relationship between atmospheric delay phase and surface elevation in the linear model can be expressed by the following formula:

[0044]

[0045] in Let k be the atmospheric delay phase at a pixel within the i-th sub-region of the N-th SLC image. i,n and c i,n These are the slope and intercept parameters in the linear model, h, respectively. i Let be the surface elevation corresponding to that pixel. The former characterizes the relationship between atmospheric delay and elevation, while the latter includes atmospheric phase components independent of elevation. Assuming that N SLC images are used to generate M-1 interferometric pairs, and the topographic and orbital phase components in each differential interferogram have been removed, and n observation points are sampled from a defined region, then the interferometric phase after atmospheric phase correction can be expressed as follows:

[0046]

[0047] in The residual phase matrix for each observation point after atmospheric correction for the i-th sub-region is given. E is the Kronecker product operator. i Let n×1 be an all-1 vector. For each interference phase matrix, This represents the observed interference phase value at each sampling point in each interferogram. i Let A be the observation matrix used for model parameter estimation within this sub-region, and let A be the M×(N-1) design matrix describing the interferometric pair mesh. Without loss of generality, the atmospheric delay of the main image is set to 0, thereby removing its corresponding column from the design matrix. iThe model parameter matrix for each SLC scene is expressed as follows:

[0048]

[0049] Solving equation (5) using the least squares method yields the model parameters for each segmented region in each SLC scene, thus calculating the atmospheric phase in each scene's differential interferogram. After obtaining the preliminary estimated atmospheric delay phase of the interferogram, it is subtracted from the original unwrapped data, and the residual phase is used as the original unwrapped phase in the second parameter iteration estimation. The linear model parameters are estimated and the phase residual is calculated again using the above method to obtain the second round of atmospheric residual phase estimation results. This process is iterated until the change in the root mean square of the phase residual in the sub-region during the two iterations is less than 1 mm (the phase needs to be converted to deformation) or the maximum number of iterations is reached (generally 3). That is, the following conditions are met:

[0050]

[0051]

[0052] in, and This represents the final result of the iterative estimation of the linear model parameters within the i-th partitioned region of the M-th interferogram, where j is the iteration number. and These are the linear model parameter values ​​obtained during the j-th iteration. By estimating the model parameters within each divided region, we can obtain the different model parameters for the atmospheric phase corresponding to each region, thereby achieving atmospheric delay phase correction at different scales.

[0053] S4. Atmospheric Delay Phase Smoothing: To ensure good spatial continuity of the estimated atmospheric phase between each divided region, the model parameters of each sub-region estimated in the previous step are assigned to the centroid pixels of these sub-regions. Then, these pixels are used as known interpolation points, and the model parameters are interpolated to each pixel in the image using the B-spline function, thereby making the overall estimated atmospheric phase smoother and preventing phase jumps between different sub-regions. After obtaining the model parameters corresponding to each pixel, the atmospheric delay phase at each pixel can be calculated using equation (3), thus obtaining the atmospheric phase of the entire image.

[0054] S5. Atmospheric Delay Correction and Deformation Result Calculation: The atmospheric delay phase, processed in the previous step, is subtracted from the original interferometric unwrapped data to obtain the atmospheric delay-corrected phase result. Then, specialized software is used to convert the phase information into deformation information, yielding an accurate atmospherically corrected deformation result.

[0055] Example 2

[0056] The study area comprises 94 Sentinel-1 images covering the China-Nepal border north of Kathmandu, Nepal, from September 2019 to September 2022. DS-InSAR technology was used for surface deformation monitoring. The terrain of the study area is mainly plateau and mountainous, with significant elevation variations, resulting in severe atmospheric delay interference in the interferograms. This embodiment uses both a traditional linear model atmospheric correction method and a terrain-segmented regional adaptive InSAR atmospheric delay correction method to perform atmospheric delay correction, and compares the correction effects of different methods on atmospheric delay.

[0057] Figure 2 The diagrams show the annual average deformation rate without atmospheric correction, and after atmospheric correction using both the linear model method and the method of this invention. It can be seen that compared to the result without atmospheric correction, the linear model method introduces a more significant terrain-related phase error in some areas. This is largely due to the failure to consider the local characteristics of atmospheric delay and the use of a global window for parameter estimation. In contrast, the atmospheric correction using the method of this invention effectively reduces the local atmospheric delay phase trend.

[0058] Theoretically, removing atmospheric delay noise can make the phase more continuous and smooth, and the smaller the residual between the estimated atmospheric phase and the original unwrapped phase, the more accurate the atmospheric phase estimation result. Therefore, in this embodiment, the mean standard deviation of all phase unwrapped data and the mean root mean square error of the estimated atmospheric delay phase are used to measure the effectiveness of different atmospheric correction methods. The mean standard deviation and mean root mean square error calculated by the two methods are shown in Table 1 below:

[0059] Table 1. Mean Standard Deviation of Atmospheric Corrected Phase

[0060]

[0061] As shown in the table, the mean standard deviation of the original phase unwrapped data is 0.9921. Atmospheric correction using the traditional linear model method actually increases the mean standard deviation of the phase to 1.0232, while the mean standard deviation of the phase obtained by the terrain-segmented regional adaptive InSAR atmospheric delay correction method is only 0.8351. Furthermore, the root mean square error also shows that the method of this invention is significantly better than the traditional linear model method. These statistical results indicate that the traditional linear model method not only fails to effectively correct atmospheric delay but also easily introduces errors, causing the mean standard deviation to increase rather than decrease. The terrain-segmented regional adaptive InSAR atmospheric delay correction method, on the other hand, can reasonably divide the study area through watersheds, better adapting to the changes in atmospheric delay at local spatial scales, thus obtaining a more accurate atmospheric delay correction effect.

[0062] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A regional adaptive multi-scale InSAR atmospheric delay correction method based on terrain segmentation, characterized in that, Includes the following steps: S1. Data acquisition and preprocessing: Acquire SAR image data and DEM data of the study area and perform preprocessing to obtain unwrapped phase data and DEM data in radar coordinate system; S2. Study Area Division: Using DEM data in radar coordinate system, the main watershed within the study area is extracted based on terrain segmentation algorithm to adaptively divide the study area. S3. Iterative estimation of model parameters: Observation points are selected in the image according to the set step size. Observation points in the deformed area are removed using a deformation rate mask. Stable observation points are used for model parameter estimation. The specific process is as follows: Using phase unwrapping data as a reference, the atmospheric delay phase is fitted with a linear model in each sub-region after division. The model parameters of each sub-region are estimated using the least squares method, thereby obtaining the multi-scale atmospheric delay phase in different regions. M-1 interferometric pairs are generated using N SLC images. The topographic phase and orbital phase components in each differential interferogram have been removed. n observation points are sampled from a certain sub-region. The relationship between atmospheric delay phase and surface elevation in the linear model is shown in the following formula: ， in, For the Nth SLC image, the first... Atmospheric delay phase at a certain pixel within a sub-region and These are the slope and intercept parameters in the linear model, respectively. The surface elevation corresponding to this pixel; the interferometric phase after atmospheric phase correction is expressed as: ， in, For the first The residual phase matrix of each observation point after atmospheric correction for each sub-region For the Kronecker product operator, Let n×1 be an all-1 vector. For each interference phase matrix, This represents the observed interference phase value at each sampling point in each interferogram. Let A be the observation matrix used for model parameter estimation within this sub-region, and let A be the M×(N-1) design matrix describing the interferometric pair mesh. Without loss of generality, the atmospheric delay of the main image is set to 0, thereby removing its corresponding column from the design matrix. The model parameter matrix for each SLC scene is expressed as follows: ， The least squares method is used to solve for the model parameters in each segmented region of each SLC scene, thereby calculating the atmospheric phase in the differential interferogram of each scene. After obtaining the preliminary estimate of the atmospheric delay phase in the interferogram, it is subtracted from the original phase unwrapping data, and the residual phase is used as new input data to iteratively estimate the model parameters. The above model parameter estimation process is repeated until the iteration stops. Finally, the model parameters in this region are the sum of the model parameters estimated in all iterations, that is, satisfying: ， in, and For the Mth interferogram The final result of iterative estimation of linear model parameters within each sub-region. For the number of iterations, and It is the first The linear model parameter values ​​obtained during the next iteration; and Substituting the relationship between atmospheric delay phase and surface elevation in the linear model, the atmospheric delay phase of the sub-region in each SLC scene can be obtained. S4. Atmospheric Delay Phase Smoothing: The estimated model parameter values ​​corresponding to each sub-region are assigned to the centroid pixel of that region as known points for interpolation. The atmospheric model parameters are fitted and interpolated to all pixels in the entire image using the B-spline function to ensure a smooth transition of atmospheric delay phase between different windows. S5. Atmospheric Delay Correction and Deformation Result Calculation: Subtract the atmospheric delay phase processed in step S4 from the initial interferometric data, and use the atmospherically corrected interferometric data to perform InSAR deformation calculation, finally obtaining the atmospherically corrected InSAR deformation result.

2. The regional adaptive multi-scale InSAR atmospheric delay correction method based on terrain segmentation according to claim 1, characterized in that, The specific process of step S2 is as follows: First, calculate the morphological gradient of the DEM data. Then, perform multi-level filtering on the gradient image using morphological opening and closing reconstruction operations at different structural scales. Next, perform a weighted average of the morphological filtering results to obtain the gradient result after multi-level filtering. Filter and reconstruct the morphological gradient of the original DEM to optimize the gradient image. The expression is shown in the following formula: ， in, For the final gradient result, and These are the opening and closing refactoring operators, respectively. The order in which opening and closing reconstruction operations are performed, structuring elements The radius r varies with Increase and grow, have n represents the total number of reconstruction operations. As the weight, satisfying , This is a morphological gradient image.

3. The regional adaptive multi-scale InSAR atmospheric delay correction method based on terrain segmentation according to claim 2, characterized in that, Morphological gradient image The expression is as follows: ， in, For the original image, It is a disk-shaped structural element. and These represent the dilation and erosion operations in grayscale morphology, respectively.

4. The regional adaptive multi-scale InSAR atmospheric delay correction method based on terrain segmentation according to claim 1, characterized in that, The condition for stopping the iteration is: the change in the root mean square value of the phase residual in the sub-region during the two iterations is less than 1 mm or the maximum number of iterations is reached.

Images