An InSAR image atmospheric phase compensation method based on self-supervised learning

Abstract

This invention discloses a self-supervised learning-based atmospheric phase compensation method for InSAR images, belonging to the field of differential interferometry (DI) technology for synthetic aperture radar. This method does not rely on external ground truth data. By establishing constraints for solving the atmospheric phase and integrating them into the network training process, it achieves self-supervised learning for atmospheric phase compensation. The method includes: constructing an atmospheric phase estimation network, which takes a time-series unwrapped phase image and related auxiliary data as input, and outputs the atmospheric phase estimation result or intermediate variables used to calculate the atmospheric phase; constructing differentiated constraint loss terms for reference points with different phase change modes based on the phase component composition of the interferometric phase; embedding physical prior constraints in both network structure and loss terms; combining various loss terms into a total loss function, and completing self-supervised training through iterative optimization to achieve atmospheric phase estimation and compensation.

Description

Technical Field

[0001] This invention belongs to the field of differential interferometry technology of synthetic aperture radar, specifically relating to an atmospheric phase compensation method for InSAR images based on self-supervised learning. Background Technology

[0002] Interferometric Synthetic Aperture Radar (InSAR) utilizes phase information from multiple SAR images to achieve high-precision measurement of surface deformation, offering advantages such as all-weather, all-day coverage, wide coverage, and non-contact measurement. However, spatiotemporal variations in atmospheric conditions cause nonlinear fluctuations in atmospheric refractive index, altering the propagation speed of radar electromagnetic waves and introducing additional atmospheric phase errors into the differential interferometric phase. This error is one of the main sources of error affecting the accuracy of deformation measurement and must be adequately compensated for.

[0003] Existing atmospheric phase compensation methods mainly fall into two categories: traditional parametric model methods and deep learning methods. Traditional parametric model methods typically establish a functional relationship between atmospheric phase and parameters such as slant range and elevation. These methods perform well in simple scenarios, but when there are drastic spatial variations in the atmosphere, the pre-defined functional form is difficult to accurately describe the complex spatial variation characteristics of atmospheric phase, resulting in limited compensation capabilities.

[0004] In recent years, deep learning methods have been introduced into atmospheric phase compensation, with the main idea being to achieve the mapping from interferometric phase images to compensated phase images through end-to-end training. However, existing deep learning methods typically employ fully supervised training, and their performance heavily relies on the completeness and accuracy of the training set. In practical applications, obtaining perfectly labeled data containing real atmospheric phase information is extremely difficult. The resulting incomplete training set will severely limit the model's generalization ability and compensation accuracy, leading to performance degradation in new scenarios. Summary of the Invention

[0005] This invention provides a self-supervised learning-based method for atmospheric phase compensation in InSAR images. This method does not rely on external ground truth data; it achieves self-supervised learning for atmospheric phase compensation by establishing constraints for solving the atmospheric phase and incorporating them into the network training process.

[0006] The technical solution of the present invention includes the following steps: S1. Construct an atmospheric phase estimation network. The network takes time-series unwrapped phase images and related auxiliary data as input and atmospheric phase estimation results or intermediate variables used to calculate atmospheric phase as output. The network adopts a deep network structure with feature extraction capabilities. S2. Based on the phase component composition of the interference phase, the differential constraint loss terms are constructed for reference points with different phase change modes by utilizing the differences in phase characteristics of each reference point in the phase image sequence. S3. Embed prior constraints into the network in the form of network structure or loss terms to limit the effectiveness of the network output; S4. Combine the prior constraint loss term in S3 with the loss terms corresponding to various reference points in S2 to form the total loss function. Perform self-supervised training on the atmospheric phase estimation network through iterative optimization to complete the estimation and compensation of atmospheric phase.

[0007] Furthermore, S1 includes the following steps: (1) Select appropriate phase image data as the input of the network. The input includes time-series unwrapped phase images. Select appropriate auxiliary data according to the task and perform preprocessing such as normalization and size adjustment on the input data. (2) Select an appropriate network structure to construct an atmospheric phase estimation network. The network adopts a deep network structure with feature extraction capability. The structure adopts a feature encoder-decoder framework and a spatiotemporal joint feature extraction module. The network output is the atmospheric phase estimation result or an intermediate variable used to calculate the atmospheric phase.

[0008] Furthermore, S2 includes the following steps: (1) Based on the spatiotemporal characteristics of the reference point phase, distinguish the reference points with different phase change modes; on this basis, for the reference point with atmospheric phase as the dominant phase component, construct the corresponding loss term by taking the difference between its unwrapped phase and the atmospheric phase estimate as the measure. (2) For a reference point whose phase is dominated by non-atmospheric change components, based on the spatiotemporal physical characteristics of its dominant phase components, a priori constraints corresponding to the characteristics of the phase change mode are applied to its phase after atmospheric compensation, and a corresponding loss term is constructed.

[0009] Furthermore, S3 includes the following steps: (1) Determine the prior constraints embedded in the form of network structure: Based on the physical properties of the network output, embed at least one physical constraint module at an appropriate position in the network to directly limit the validity of the output, including a numerical range limiting module that limits the output value to a preset range, or an invalid region zeroing module that sets the output to zero in the absence of an effective observation region. (2) Determine the prior constraints embedded in the form of loss terms: Based on the physical laws that the network output should follow, construct at least one physical prior constraint loss term, including a spatial continuity constraint loss term that encourages the network output to transition smoothly in the spatial dimension, or a temporal continuity constraint loss term that encourages the network output to change continuously in the temporal dimension.

[0010] Furthermore, S4 includes the following steps: (1) The loss terms constructed in S2 for different phase change mode reference points are weighted and summed with the physical prior constraint loss terms constructed in S3 to form the total loss function for network self-supervised training; (2) Set up an iterative optimization strategy, perform self-supervised training on the atmospheric phase estimation network based on the total loss function, update the network parameters through gradient descent or its improved method, and the training process may include a warm-up phase to determine the weighting coefficients of each loss term until the total loss function converges or meets the preset stopping condition. (3) Select the network model with the best performance during the training process as the final network and use the network to complete the atmospheric phase estimation. If the network in S1 is set to directly output the atmospheric phase estimation result, the result can be used directly. If the network in S1 is set to output intermediate variables, the intermediate variables need to be converted into atmospheric phase estimates through a preset physical relationship model, and the atmospheric phase estimates are removed from the original unwrapped phase to complete the atmospheric phase compensation.

[0011] Beneficial effects: 1. The method of this invention achieves atmospheric phase compensation based on a self-supervised learning framework, without the need for external ground truth labeled data. By embedding the physical constraints of atmospheric phase into the network training process, it effectively overcomes the problems of dependence on the completeness of the training set and insufficient generalization ability of existing fully supervised deep methods. 2. The method of this invention divides the reference points according to the differences in phase characteristics and constructs differentiated constraint loss terms for reference point sets with different phase change modes. At the same time, it incorporates physical priors through both structural embedding and loss embedding to enhance the physical rationality and numerical stability of the solution results and improve the applicability and compensation accuracy of the method in complex scenarios. Attached Figure Description

[0012] Figure 1 This is a schematic diagram of the atmospheric phase estimation network. Figure 2 This is a structural diagram of the numerical range limiting module; Figure 3 A schematic diagram illustrating the division of invalid regions; Figure 4 This is a schematic diagram of the self-supervised training iterative process; Figure 5 This is the overall flowchart of the present invention. Detailed Implementation

[0013] The present invention will now be described in detail with reference to the accompanying drawings and embodiments.

[0014] This invention provides an atmospheric phase compensation method based on self-supervised learning, the process of which is as follows: Figure 5As shown, the specific steps include the following: S1. Construct an atmospheric phase estimation network. The network takes time-series unwrapped phase images and related auxiliary data as input, and outputs atmospheric phase estimation results or intermediate variables used to calculate atmospheric phase. The network adopts a deep network structure with feature extraction capabilities. The specific process is as follows: S101. Select appropriate input data and perform preprocessing.

[0015] like Figure 1 As shown, in one embodiment of the present invention, a normalized slant range image, an effective point mask, and an unwrapped phase image sequence of length D are selected as network inputs. Each image is divided into H×W grids to estimate the refractive index change (i.e., using the refractive index change as an intermediate variable for atmospheric estimation). To enhance the network's ability to extract global spatial information, the spatial dimensions of the input images are adjusted to a preset multiple of the grid size (e.g., 8H×8W) through resampling.

[0016] S102. Construct the main network structure.

[0017] like Figure 1 As shown, in this embodiment of the invention, the main network adopts a deep structure with feature extraction capabilities, including a downsampling module, a residual module, and a convolutional module. The downsampling module consists of a 3D convolutional layer (3×3×3 kernel size, 1×2×2 stride, 1×1×1 padding) followed by a ReLU activation function, compressing the spatial size of the feature map layer by layer (height and width halved, depth unchanged). The residual module uses identity mapping residual connections. Each residual module contains a 3D convolutional layer with a 3×3×3 kernel, 1×1×1 stride, and 1×1×1 padding, followed by a ReLU activation function. The input is added to the convolutional output via an adder, thereby enhancing nonlinear expressive power and mitigating gradient vanishing. The convolutional module consists of a 3D convolutional layer with a 3×1×1 kernel, 1×1×1 stride, and no padding, used to adjust the number of feature map channels to match the target output. As described in S101, the network output is a gridded refractive index change, with output dimensions of (D, H, W), corresponding to the depth, length, and width dimensions of the output result, respectively.

[0018] S2. Phase component construction based on interferometric phase: Utilizing the differences in phase characteristics of each reference point in the phase image sequence, differentiated constraint loss terms are constructed for reference points with different phase change modes. The specific process is as follows: S201. Based on the spatiotemporal characteristics of the reference point phase, distinguish the reference points for different phase change modes.

[0019] In one embodiment of this invention for ground-based InSAR applications, the classification of reference points can be achieved by analyzing the coherence or stability characteristics of their phase over time series. For example, indicators such as amplitude deviation and phase stability can be used, or residuals can be calculated after preliminary fitting of the atmospheric phase using a parametric model, and reference points can be divided into deformed and non-deformed categories based on a preset threshold. Reference points of different categories will be used to subsequently construct differentiated constraint loss terms.

[0020] S202. Determine the conversion relationship from network output to atmospheric phase estimate.

[0021] If the network in S1 is set to directly output the atmospheric phase estimation result, the result can be used directly without additional conversion. If the network in S1 is set to output an intermediate variable, the intermediate variable needs to be converted into an atmospheric phase estimate value through a preset physical relationship model. In one embodiment of the present invention, the network output is the intermediate variable of the gridded refractive index change. According to the physical law of atmospheric delay, the propagation path of the radar signal from the radar to the reference point is divided into several segments by each grid. The atmospheric phase of any reference point can be expressed as the sum of the products of the path length of each segment and the corresponding grid refractive index change:

[0022] in The total number of grid cells. To reach the reference point The intercept of the propagation path in the k-th grid. Let be the refractive index change of the k-th grid (assuming the refractive index change within each grid is a constant). Using this physical relationship model, the atmospheric phase estimate for each reference point can be calculated from the gridded refractive index change output by the network.

[0023] S203. For the set of reference points with atmospheric phase as the dominant phase component, construct the first type of loss term.

[0024] In one embodiment of the present invention for a ground-based InSAR application, as described in S201, for non-deformation reference points, the dominant phase component is the atmospheric phase. For such reference points, the difference between the unwrapped phase and the atmospheric phase estimate should be as small as possible. Accordingly, the difference between the two is used as a metric to construct the loss term corresponding to this set. In one embodiment of the present invention, this loss term may be in the form of mean squared error:

[0025] in For this type of reference point set, For reference point The untangling phase, This is the atmospheric phase estimate for the reference point obtained in S202. This loss term drives the network to learn the spatial distribution characteristics of the atmospheric phase by minimizing the atmospheric phase estimation bias.

[0026] S204. For the set of reference points whose phase components are dominated by non-atmospheric variation components, construct a second type of loss term.

[0027] For such reference points, the unwrapped phase includes not only the atmospheric phase but also non-atmospheric phase components such as deformation. Given that these non-atmospheric phase components typically exhibit specific spatiotemporal variation patterns, it is not appropriate to directly construct constraints using the residuals between the unwrapped phase and the atmospheric phase estimates. Therefore, the atmospherically compensated phase for such reference points is calculated first:

[0028] in and These are the reference points in the d-th image. The unwrapped phase and atmospheric phase estimates are obtained. Then, based on the spatiotemporal physical characteristics of the dominant phase component, corresponding prior constraints are applied to the compensated phase to construct the loss term.

[0029] In one embodiment of the present invention for a ground-based InSAR application, as described in S201, for the deformation reference points, the dominant phase component is the deformation phase, and deformation typically exhibits smooth changes over a short period. Based on this physical characteristic, the second-order difference along the time dimension is calculated for the atmospherically compensated phase sequence of each reference point in the set, and its squared mean is used as the loss term:

[0030] in This is the set of reference points for this type. The loss term encourages the minimization of the second-order temporal difference of the compensated phase sequence, prompting the network to compensate for a reasonable, time-smoothly varying deformation phase in the deformation region, thereby achieving effective separation of the atmosphere and deformation components.

[0031] S205. If there is a set of reference points for other phase change modes, then construct corresponding constraint loss terms based on the physical characteristics of their dominant phase components. The specific form of each loss term can be flexibly set according to the spatiotemporal prior characteristics of the phase component, and is not limited to the two forms mentioned above.

[0032] S3. Embed prior constraints into the network in the form of network structure or loss terms to limit the effectiveness of the network output. The specific process is as follows: S301. Determine the prior constraints embedded in the form of a network structure.

[0033] Based on the physical properties of the network output, physical constraint modules are embedded at appropriate locations within the network to directly limit the effectiveness of the network output. In one embodiment of the invention, the network output is a meshed refractive index change, to which two types of structural embedding constraints are applied: One is the numerical range constraint. For example... Figure 2 As shown, a numerical range limiting module is connected in series at the end of the network. This module first uses an activation function (such as Sigmoid) to normalize the values ​​of the intermediate outputs of the network to the [0,1] interval, and then uses the upper and lower bounds of the refractive index change estimated from the data or set according to prior knowledge to linearly map the values ​​of each depth dimension (corresponding to each image) of the network output:

[0034] in This is the original output for the d-th depth dimension. and These represent the upper and lower bounds of the refractive index change in the corresponding depth dimension. After processing by this module, the refractive index change of each grid is strictly limited to a physically reasonable range.

[0035] The second is invalid region constraints. For example... Figure 3 As shown, based on the geometric coverage relationship of the radar signal propagation path, grids in the meshed region that are not traversed by any propagation path from the reference point to the radar are determined as invalid regions, and a binary mask is generated accordingly. (Dimensions are H×W). Masking operation at the network end. This forces the output value corresponding to the invalid region to be set to zero, ensuring that the network output only takes effect within the valid region that is physically supported by observations.

[0036] S302. Determine the prior constraints embedded in the form of loss terms.

[0037] Based on the physical laws that the network output should follow, a physical prior constraint loss term is constructed to encourage the network output to satisfy the corresponding physical characteristics. In one embodiment of the invention, a spatial continuity constraint loss term is constructed, using a total variational form:

[0038] in Let be the total variational function. For the first The network output corresponds to each image. The total variation measures the spatial smoothness of the output by calculating the sum of the absolute differences between adjacent pixels. This loss term encourages a smooth transition in the spatial dimension of the network output, suppresses local abrupt changes caused by observation noise, and enhances the stability of the solution. Furthermore, other forms of prior constraint loss terms, such as temporal continuity constraints, can be introduced as needed.

[0039] S4. Combine the prior constraint loss term from S3 with the loss terms corresponding to various reference points in S2 to form the total loss function. Iteratively optimize this function to perform self-supervised training on the atmospheric phase estimation network, thus completing the estimation and compensation of the atmospheric phase. The specific process is as follows: S401. Construct the total loss function.

[0040] The loss terms (main loss terms) constructed in S2 for different phase change mode reference points are weighted and summed with the physical prior constraint loss terms (auxiliary loss terms) constructed in S3 in the form of loss terms, thus forming the total loss function for network self-supervised training. In one embodiment of the present invention, the total loss term is:

[0041] in The loss term corresponds to the set of non-deformation reference points where atmospheric phase is the dominant component. This is the loss term corresponding to the set of deformation-type reference points where deformation phase is the dominant component. For spatial continuity constraint loss term, , , These are the weighting coefficients for each loss term. The weighting coefficients can be set according to the magnitude relationship of each loss term.

[0042] S402. Set the iterative optimization strategy and execute self-supervised training.

[0043] The atmospheric phase estimation network is self-supervised trained based on the total loss function, and the network parameters are updated using gradient descent or its improved methods (such as the Adam optimizer). In one embodiment of the present invention, as follows... Figure 4 As shown, the training process can include a warm-up phase, a formal training phase, and a closing phase.

[0044] The purpose of the warm-up phase is to confirm the weighting coefficients of the loss term. This phase can begin by setting auxiliary loss terms (such as...). The weight of ) is zero. After several iterations, the main loss term () is recorded. The magnitude relationship between the main loss term and the auxiliary loss term is used to determine reasonable weighting coefficients. For example, the spatial smoothing loss can be set about an order of magnitude lower than the main loss term to ensure that it plays a constraining role without dominating the gradient direction.

[0045] The purpose of the formal training phase is to obtain the optimal parameters of the network. This phase iteratively optimizes the network based on the complete total loss function. During training, the learning rate can be automatically adjusted according to the loss decrease: the learning rate is halved when the loss does not decrease for several consecutive rounds; training is terminated when the preset maximum number of iterations is reached or the loss does not decrease for several consecutive rounds, meeting the early stopping condition. The network parameters that achieve the optimal loss are automatically saved during training.

[0046] At the end of the training phase, the optimal parameters of the network obtained during the formal training are loaded and the network output is obtained as the final solution of the network for the atmospheric parameters.

[0047] S403, complete atmospheric phase compensation.

[0048] The trained network is used to estimate the atmospheric phase. If the network in S1 is set to directly output the atmospheric phase estimate, it is directly used as the basis for atmospheric phase compensation. If the network in S1 is set to output intermediate variables, the intermediate variables need to be converted into atmospheric phase estimates using the preset physical relationship model described in S202. In one embodiment of the invention, for the d-th image, the atmospheric phase estimate for each reference point is calculated using the gridded refractive index change output by the network through an atmospheric delay integral model. Finally, the atmospheric phase estimate is subtracted from the original unwrapped phase to obtain the atmospherically compensated phase, thus completing the atmospheric phase compensation.

[0049] In summary, the above are merely preferred embodiments of the present invention and are not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A self-supervised learning-based InSAR image atmospheric phase compensation method, characterized in that, Includes the following steps: S1. Construct an atmospheric phase estimation network. The network takes time-series unwrapped phase images and related auxiliary data as input and atmospheric phase estimation results or intermediate variables used to calculate atmospheric phase as output. The network adopts a deep network structure that includes feature extraction capabilities. S2. Based on the phase component composition of the interference phase, the differential constraint loss terms are constructed for reference points with different phase change modes by utilizing the differences in phase characteristics of each reference point in the phase image sequence. S3. Embed prior constraints into the network in the form of network structure or loss terms to limit the effectiveness of the network output; S4. Combine the prior constraint loss term in S3 with the loss terms corresponding to various reference points in S2 to form the total loss function. Perform self-supervised training on the atmospheric phase estimation network through iterative optimization to complete the estimation and compensation of atmospheric phase.

2. The method according to claim 1, characterized in that, S2 specifically includes: (1) Based on the spatiotemporal characteristics of the reference point phase, the reference points are divided into different phase change modes; for the set of reference points with atmospheric phase as the dominant phase component, the difference between the unwrapped phase and the atmospheric phase estimate is used as a metric to construct the corresponding loss term; (2) For a set of reference points whose dominant phase components are non-atmospheric change components, based on the spatiotemporal physical characteristics of their dominant phase components, a priori constraints corresponding to the characteristics of the phase change mode are applied to their phase after atmospheric compensation, and a corresponding loss term is constructed.

3. The method according to claim 2, characterized in that, For a set of reference points where atmospheric phase is the dominant phase component, the first type of loss term adopts the mean square error form: ； in For this type of reference point set, For reference point The untangling phase, For the network at reference point x i The atmospheric phase estimate is obtained; the loss term drives the network to learn the spatial distribution characteristics of the atmospheric phase by minimizing the atmospheric phase estimation bias.

4. The method according to claim 2, characterized in that, For a set of reference points whose phase components are dominated by non-atmospheric variation components, the atmospherically compensated phase of each reference point is first calculated: ； in and These are the reference points in the d-th image. The unwrapped phase and atmospheric phase estimates are obtained; then, the second-order difference is calculated along the time dimension of the atmospherically compensated phase sequence, and its squared mean is used as the second type of loss term: ； in For this type of reference point set, D The total number of image sequences; the loss term encourages the minimization of the second-order temporal difference of the compensated phase sequence, prompting the network to compensate for the temporally smooth atmospheric phase in the deformation region, thereby achieving effective separation of the atmosphere and deformation components.

5. The method according to claim 1, characterized in that, S3 specifically includes: (1) Determine the prior constraints embedded in the form of network structure: Based on the physical properties of the network output, embed at least one physical constraint module at an appropriate position in the network. The physical constraint module includes a numerical range limiting module that limits the output value to a preset range, or an invalid region zeroing module that sets the output to zero when there is no effective observation region. (2) Determine the prior constraints embedded in the form of loss terms: Based on the physical laws that the network output should follow, construct at least one physical prior constraint loss term, which includes a spatial continuity constraint loss term that encourages the network output to transition smoothly in the spatial dimension, or a temporal continuity constraint loss term that encourages the network output to change continuously in the temporal dimension.

6. The method according to claim 5, characterized in that, The calculation method of the numerical range limitation module is as follows: ； in This is the original output for the d-th depth dimension. and These represent the upper and lower bounds of the refractive index change for the corresponding depth dimension, respectively, with Sigmoid as the activation function. After processing by the module, the values ​​output by the network in each depth dimension are strictly limited to a physically reasonable range.

7. The method according to claim 5, characterized in that, The invalid region zeroing module processes the data as follows: Based on the model used and the observation data, invalid regions in the network output are identified, and the output values ​​corresponding to the invalid regions are set to zero, so that the network output only takes effect in valid regions with observation support.

8. The method according to claim 5, characterized in that, The spatial continuity constraint loss term is expressed in total variational form: ； in Let be the total variational function. For the first Network output corresponding to each image ，D The total number of image sequences is represented by the total variational function, which measures the spatial smoothness of the output by calculating the sum of the absolute values ​​of the differences between adjacent pixels. The loss term suppresses local abrupt changes caused by observation noise and enhances the spatial stability of the solution.

9. The method according to claim 1, characterized in that, S4 specifically includes: (1) Constructing the total loss function: The total loss function is obtained by weighted summing of the main loss terms constructed in S2 for different phase change mode reference points and the physical prior constraint auxiliary loss terms constructed in S3 in the form of loss terms. ； in The loss term corresponds to the set of non-deformation reference points where atmospheric phase is the dominant component. This is the loss term corresponding to the set of deformation-type reference points where deformation phase is the dominant component. For spatial continuity constraint loss term, , , Equations are weighting coefficients for each loss term, and these weighting coefficients are set according to the magnitude relationship of each loss term. (2) Set up an iterative optimization strategy and perform self-supervised training on the atmospheric phase estimation network based on the total loss function. The training process includes a warm-up phase, which is used to determine the weighting coefficients by recording the magnitude relationship of each loss term and to update the network parameters by the gradient descent method until the total loss function converges or meets the preset stopping condition. (3) Select the network model with the best performance during the training process as the final network to complete the estimation and compensation of atmospheric phase.

10. The method according to claim 1, characterized in that, The deep network structure adopts a feature encoder-decoder framework and a spatiotemporal joint feature extraction module. The network includes a three-dimensional convolutional downsampling module, a residual module, and a convolution module. It takes temporally unwrapped phase image sequences, normalized slant range images, and effective point masks as inputs and gridded refractive index changes as outputs.

Images