A Method for Optimizing the Phase of Distributed Targets in Temporal InSAR

Abstract

The present invention discloses a method for optimizing the phase of a distributed target in temporal InSAR, which relates to the technical field of synthetic aperture radar interferometry data processing, and includes the following steps: obtaining N scenes of InSAR images; applying the CMP algorithm to select the homogeneous pixel set Ω of each pixel p, using the selected homogeneous pixel set Ω to estimate the complex covariance matrix, pre-estimating the empirical coherence coefficient based on the InSAR dataset information, jointly determining the parameter s by the empirical coherence coefficient and the homogeneous pixel set Ω, performing iterative correction on the coherence amplitude matrix through the coherence deviation corrector and the homogeneous pixel set Ω, and combining the coherence amplitude matrix after deviation correction with the EMI phase evaluator to obtain the finally optimized phase sequence. The present invention significantly reduces the deviation of the coherence matrix and obtains a high-precision reconstructed phase sequence, especially in the case of a small number of homogeneous pixel sets.

Description

Technical Field

[0001] The present invention relates to the technical field of synthetic aperture radar interferometry data processing, and particularly to a method for optimizing the phase of distributed targets in time-series InSAR. Background Art

[0002] Among various multi-temporal interferometric synthetic aperture radar (MT-InSAR) techniques, persistent scatterer InSAR (PSInSAR) is a high-precision geodetic method for monitoring the displacement of the Earth's surface. The basic idea of this technique is to identify and analyze persistent scatterers (PSs) with stable phases during the entire observation period. These scatterers usually correspond to artificial targets. Therefore, a major limitation of PSInSAR is the low spatial density of measurement points, especially in non-urban environments.

[0003] Distributed scatterers (DSs) are used to increase the density of measurement points. Different from PSs, DSs usually cover multiple small targets with the same scattering behavior in one resolution cell. These targets have geometric and temporal decorrelations and are usually noisier than PSs. SqueeSAR performs maximum likelihood estimation to reduce the decorrelation of DSs, which is called phase estimation. In addition, several phase estimation methods have been proposed, among which the maximum likelihood estimation of interferometric phase based on eigen-decomposition (EMI) is an effective method that takes into account both estimation accuracy and computational efficiency. By estimating the optimized single-master image phase sequence from all available interferograms, the decorrelation in this process is significantly reduced. Then, using traditional PSInSAR processing, the estimated phase sequence can be used to invert geophysical signals. The EMI phase estimation algorithm is similar to other phase estimation algorithms and uses a variant of the coherence matrix as the weight. The influence of the coherence matrix error on its accuracy is significant.

[0004] Due to the difficulty in accurately estimating the coherence matrix, the assumption that the sample coherence matrix is close to the true one is broken, and the performance of EMI drops severely. Previous studies have improved the accuracy of DS phase estimation by reducing the coherence error. However, due to the limited correction from simple bias correction estimators, the improvement is insufficient. In order to reduce the loss of deformation monitoring accuracy, there is still a lack of high-precision phase estimation. Summary of the Invention

[0005] The present invention provides a method for optimizing the phase of distributed targets in time-series InSAR, aiming to solve the drawbacks existing in the prior art.

[0006] To achieve the above object, the present invention provides the following technical solution: A method for optimizing the phase of distributed targets in time-series InSAR, comprising the following steps:

[0007] Obtain N InSAR images and construct an InSAR dataset;

[0008] Select a homogeneous pixel set Ω for each pixel p of the image through the CMP algorithm;

[0009] Use the complex covariance matrix estimator and the selected homogeneous pixel set Ω to estimate the complex covariance matrix And for the complex covariance matrix Take the modulus to obtain the coherence amplitude matrix

[0010] Pre-estimate the empirical coherence coefficient based on the information of the InSAR dataset

[0011] Through the empirical coherence coefficient And the homogeneous pixel set Ω to determine the empirical constant s;

[0012] Obtain the coherence coefficient bias corrector by combining the empirical constant s with the natural logarithm moment of the pixel

[0013] Through the coherence coefficient bias corrector And the homogeneous pixel set Ω perform iterative correction on the coherence amplitude matrix ;

[0014] Input the coherence amplitude matrix after bias correction Into the EMI phase evaluator to obtain the final optimized phase sequence.

[0015] Preferably, the formula of the complex covariance matrix estimator is:

[0016]

[0017] Among them, L is the number of homogeneous pixels, s(k) represents the normalized complex vector, and H represents the conjugate transpose.

[0018] Preferably, the empirical coherence coefficient Is pre-estimated as:

[0019]

[0020] Among them,

[0021]

[0022] Is the thermal decoherence, where the signal-to-noise ratio SNR depends on the system parameters and is set to a constant value of 12 dB;

[0023] Is the geometric decoherence, where B ⊥ And B ⊥max Are the vertical baseline and the critical baseline respectively;

[0024] is the time decoherence, where B T and B Tmax are the baseline and the decoherence rate in time, respectively.

[0025] Preferably, the natural logarithm moment of the pixel is approximated as:

[0026]

[0027] where s is an empirical constant and L is the number of homogeneous pixels.

[0028] Preferably, the coherence coefficient corrector is:

[0029]

[0030] where

[0031]

[0032] Preferably, the empirical constant s is estimated by the following formula:

[0033]

[0034] where round the result of it to the largest integer not greater than it.

[0035] Preferably, the specific formula for inputting the coherence amplitude matrix after bias correction into the EMI phase evaluator is as follows:

[0036]

[0037] where, replace the coherence amplitude matrix after correction is the estimated consistent time series phase, and the symbol ° represents the Hadamard product and use the Lagrange multiplier method to solve the

[0038] Preferably, EMI performs matrix inversion operation. For the case where the coherence amplitude matrix is not positive definite, pseudo-inverse processing is adopted.

[0039] The present invention has the following beneficial effects compared with the prior art:

[0040] The new estimator can simultaneously incorporate the information of the interference coherence coefficient and the multi-look number, and achieve significant bias correction for each element of the coherence matrix. Then, the coherence matrix after bias correction is combined with EMI to obtain an optimized phase sequence. Due to the low bias of the coherence matrix, this method obtains a phase sequence with high accuracy, especially in the case of a small number of SHP sets. Description of the Drawings

[0041] Figure 1 is the expectation of the sampling coherence coefficient under different sample numbers L of the present invention The difference graph from the true coherence coefficient γ

[0042] Figure 2 is the corrected coherence coefficient under different sample numbers L of the present invention The difference graph from the true coherence coefficient γ

[0043] Figure 3 is the algorithm flow chart of the present invention

[0044] Figure 4 is the coherence matrix estimated based on a 5×5 pixel window of the present invention: (a) true data; (b) sampling coherence matrix; (c) bias-corrected coherence matrix

[0045] Figure 5 is for the estimated coherence matrix under different multi-look numbers of the present invention: (a) mean value; (b) standard deviation

[0046] Figure 6 is the standard deviation of the phase sequence residuals estimated based on different coherence matrices of the present invention

[0047] Figure 7 is for the volcanic research area of the present invention: (a) optical image; (b) average intensity map from 40 Sentinel-1 images

[0048] Figure 8 is for the present invention: (a) original single-look interferogram, with the reconstructed phase using: (b) sampling coherence coefficient, (c) general corrected coherence coefficient, (d) corrected coherence coefficient proposed by the present invention

[0049] Figure 9 is for the characteristics of the original and reconstructed interferometric phases in the present invention: (a) horizontal profile; (b) vertical profile

[0050] Figure 10 is the average value of the temporal coherence coefficient under different SHP numbers of the present invention Detailed Embodiments

[0051] The following further describes the detailed embodiments of the present invention in conjunction with the drawings. The following embodiments are only used to more clearly illustrate the technical solutions of the present invention and cannot be used to limit the protection scope of the present invention

[0052] The present invention provides a method for optimizing the phase of distributed targets in temporal InSAR. First, a theoretical analysis of the coherence matrix error is carried out. On this basis, by considering the factors of the interference coherence coefficient and the number of looks, a method for systematically reducing the coherence coefficient error is proposed. The coherence matrix after bias correction is combined with EMI to obtain the optimized phase sequence. Tests are carried out on simulation and real Sentinel-1 data to verify the effectiveness of the algorithm.

[0053] Among them, for any pixel 2 in the search window and the reference pixel 1 in the CMP algorithm, the similarity is measured by the following formula:

[0054]

[0055] Where

[0056]

[0057] s j is the complex signal of the j-th SAR image. Based on the similarity obtained above, according to the prior threshold, the homogeneous pixel set of pixel 1 is determined.

[0058] In the traditional method, for the single-look complex images registered for N scenes, the normalized complex vector s = [s1, s2,..., sN] T is usually assumed to follow a zero-mean N-variable complex circular normal distribution. The complex covariance matrix CCM can be estimated as:

[0059]

[0060] Ω is the set of statistically homogeneous pixels SHP with L adjacent pixels, L is the number of homogeneous pixels, s(k) represents the normalized complex vector, and H represents the conjugate transpose.

[0061] Although can be used for phase unwrapping and deformation monitoring, it is affected by spatio-temporal decorrelation. Phase estimation is to minimize the influence of decorrelation by using the statistical characteristics of the data. As a recently proposed estimator, EMI has been proven to be the current technique with the best accuracy. Replace with the corrected coherence amplitude matrix θ = [0, θ2,..., θ N T is the estimated consistent time series phase. The specific formula for combining the bias-corrected coherence amplitude matrix with the EMI phase estimator is as follows:

[0062]

[0063] Among them, the symbol represents the Hadamard product, and the Lagrange multiplier method with the minimum eigenvalue is used to solve​

[0064] As in previous studies, a similar mathematical model is also applicable to Equation (2). The most fundamental difference is the weight matrix which assigns different weights to different interferograms. However, Equation (1) is a suboptimal method, and the estimated coherence magnitude is biased. In addition, the matrix inversion operation amplifies the bias. These factors seriously impair the accuracy of EMI.

[0065] Amplitude deviation statistic:

[0066] The empirical coherence coefficient is an element of the coherence magnitude matrix whose probability density function pdf is:

[0067]

[0068] where γ is the true coherence coefficient magnitude, L is the number of samples, and F is the generalized hypergeometric function.

[0069] Expected value The analytical expression of

[0070]

[0071] is derived as:

[0072] Figure 1 where Γ is the gamma function. represents the relationship between the coherence coefficient deviation and the true coherence coefficient and the number of samples L. It can be seen that

[0073] is biased towards higher values, and the deviation increases with the decrease of γ and L. Therefore, the correction should consider the influence of different γ and L.

[0074] Specifically, the deviation correction of the consistency amplitude:

[0075]

[0076] The above equation can also be called the logarithmic moment E((ln(γ 2 ))) s .

[0077] The improved phase estimation algorithm in a temporal InSAR distributed target phase optimization method is described in detail below, which includes the following steps:

[0078] Step 1: Obtain N scenes of InSAR images and construct an InSAR dataset.

[0079] Step 2: Select a set of homogeneous pixels Ω for each pixel p of the image through the CMP algorithm.

[0080] Step 3: Use the complex covariance matrix estimator and the selected set of homogeneous pixels Ω to estimate the complex covariance matrix and take the modulus of the complex covariance matrix to obtain the coherence amplitude matrix

[0081] Step 4: Pre-estimate the empirical coherence coefficient based on the information of the InSAR dataset

[0082] Step 5: Determine the empirical constant s through the empirical coherence coefficient and the set of homogeneous pixels Ω.

[0083] Step 6: Obtain the coherence coefficient bias corrector by combining the empirical constant s with the natural logarithm moment of the pixel

[0084] Step 7: Perform iterative correction on the coherence amplitude matrix using the coherence coefficient bias corrector and the set of homogeneous pixels Ω.

[0085] Step 8: Input the coherence amplitude matrix after bias correction into the EMI phase evaluator to obtain the final optimized phase sequence.

[0086] Specifically, use the complex covariance matrix estimator, that is, formula (1), and the selected set of homogeneous pixels Ω and substitute them into formula (1) to estimate the complex covariance matrix

[0087] Specifically, by taking the inverse of formula (5), a coherence coefficient bias corrector is obtained:

[0088] The coherence coefficient bias corrector is:

[0089]

[0090] where

[0091]

[0092] As we know it is affected by both γ and L. Here, these two factors are also considered. According to numerical calculations, the estimated bias is as Figure 2As shown. It can be seen that the value of s has an important influence on the estimation deviation. Generally speaking, in order to obtain a more accurate result, when the true coherence coefficient is low and the number of samples is small, a larger value should be assigned to s. When s = 1, this estimator has a simple form and has been tested in previous studies.

[0093] Because the integral of formula (5) is difficult to implement in practice, the logarithmic moment can be replaced by the statistical average of a large number of independent samples. Using the SHP set here, the logarithmic moment of the pixel can be approximated as:

[0094]

[0095] where s is an empirical constant and L is the number of homogeneous pixels.

[0096] Based on the theoretical relationship between formulas (5) and (6), the empirical constant s can be estimated by the following formula:

[0097]

[0098] where [] represents the floor operation. Since the true coherence is unknown, is its substitute. Considering that has a relatively large variance, taking it as will lead to an inaccurate result. Here, an empirical coherence coefficient is pre-estimated as:

[0099]

[0100] where,

[0101]

[0102] is the thermal decoherence where the signal-to-noise ratio SNR depends on the system parameters and is set to a constant 12 dB; is the geometric decoherence, where B ⊥ and B ⊥max are the vertical baseline and the critical baseline respectively; is the temporal decoherence, where B T and B Tmax are the temporal baseline and the decoherence rate respectively. B ⊥max and B Tmax are set to 1100 meters and 200 days respectively. The empirical coherence coefficient is fixed throughout the interferogram. Therefore, the computational burden is not too heavy. In addition, it has different values in different interferogram pairs, and all elements in the matrix can be made consistent by using their corresponding average values for correction.

[0103] The proposed coherent error correction method can be called a general strategy and can be applied to other DS phase estimators. In practice, Equation (8) is estimated using a homogeneous pixel set, and it is difficult to meet the assumption of a large number of samples due to the limited number of samples. Here, this problem can be alleviated by iteratively performing Step 7. Considering the balance between accuracy improvement and increased computation time, the iteration is usually performed twice. In addition, EMI requires matrix inversion operations. For the case where the coherent amplitude matrix is not positive definite, a pseudo-inverse is used to solve the occasional problems.

[0104] Example:

[0105] 1. Test on the synthetic InSAR dataset:

[0106] Simulated data was synthesized according to the coherence matrix. Consider an exponential decorrelation model:

[0107] γ = (γ0 - γ ∝ )·γ thermal ·γ geom ·γ temp +γ ∝ (12)

[0108] where γ0 and γ ∝ are the short-term and long-term coherences, and are set to 0.7 and 0.03. γ thermal , γ geom and γ temp are thermal decorrelation, geometric decorrelation, and temporal decorrelation respectively. The signal-to-noise ratio was set to 12 dB, resulting in γ thermal = 0.92. The vertical baseline was normally distributed with a mean of zero and a standard deviation of 50 m, and the critical baseline was 1100 m. Similar to Sentinel 1, the revisit period was set to 12 days, and the decorrelation rate was 200 days. The interferometric phase was only related to the velocity and was 1 mm / year. Based on the coherence matrix, 40 images were simulated using a similar method. The size of the simulated image was 100×100 pixels. Since is an average coherence coefficient in Equation (11), SNR, B ⊥max and B Tmax were set to 12 dB, 1100 m, and 200 days respectively.

[0109] Figure 4 is the result of the coherence matrix estimation. Comparing the true coherence matrix 4(a) and the sampled coherence matrix 4(b), it can be clearly seen that there is an overestimation. By performing the bias correction method proposed in the present invention, the coherence matrix has been significantly improved.

[0110] For further quantitative evaluation, the mean and standard deviation of all pixels of N(N - 1) / 2 interferograms were calculated and shown in Figure 5Among them. Considering that the number of homogeneous pixels of DS usually has a threshold, for example, L>20, two search windows are set here, which are 5×5 and 7×7 respectively. Observe Figure 5 (a), and it can also be seen that the sampling coherence coefficient is overestimated. When the number of looks becomes smaller, the overestimation is more serious. When using the traditional simple bias correction method, it can be seen that the bias is reduced. Finally, by using the method proposed in the present invention, a coherence coefficient close to the true value is obtained, realizing significant coherence coefficient bias correction. This is mainly because the adopted s takes into account the true coherence coefficient and the magnitude of the multi-look number of looks. In addition, by comparing the standard deviation image 5(b), it can be seen that the proposed method obtains a lower value. The more accurate mean and the lower standard deviation together indicate the effectiveness of the proposed method in correcting the coherence matrix bias.

[0111] Based on different coherence amplitude matrices, the EMI algorithm is used for phase estimation. We statistically calculated the standard deviation of the residuals, that is, the difference between the estimated phase sequence and the true value, as Figure 6 shown. The Cramér-Rao lower bound CRLB provides the theoretically achievable highest precision. Generally, for different methods, when the multi-look number of looks is larger, the accuracy of phase estimation is higher. This is because a more accurate complex covariance matrix is obtained, including more accurate amplitude and more accurate phase. By comparing the estimation results of different coherence coefficients, it can be found that the coherence coefficient corrected in the present invention has better accuracy than the coherence coefficient corrected by the existing method and is close to the true coherence coefficient, and the average residual is reduced by 0.11 rad. This shows that the coherence amplitude matrix plays an important role in the estimated phase sequence. In addition, compared with the traditional sampling coherence coefficient , the phase estimated by the method of the present invention has significantly higher accuracy, and the average residual is reduced by 0.47 rad, which is equivalent to the deformation deviation of 2.4 mm in the C-band of Sentinel-1, which may lead to obvious deformation differences.

[0112] 2. Test on real InSAR datasets

[0113] The 40 scenes of Sentinel-1 images of the Alcedo volcano are used to test the performance of the algorithm. Figure 7 The optical and average SAR intensity data are shown, revealing the diversity of the covered ground objects in the scene. The Sentinel-1 images are obtained in wideband mode and VV polarization mode, and are obtained in descending orbit mode from January 2021 to May 2022. The test area is limited to 650×2400 pixels. The CMP algorithm is used to select homogeneous pixels in a 7×13 window, and pixels exceeding 20 are selected for DS phase estimation.

[0114] The performance of the method was visually evaluated by examining the original and reconstructed interferograms with the longest time baseline of 288 days. The results are as Figure 8 shown.

[0115] To better compare the details, a crater area delimited by a rectangular box was enlarged for display. Comparing Figure 8 (a) and Figure 8 (b), it can be seen that the traditional EMI phase estimation method has a significant noise suppression effect. When using the existing deviation-corrected coherence coefficient a smoother interferometric phase map was obtained, but there were still some noise points. Using the corrected coherence coefficient proposed in the present invention it can be seen that the noise points were further reduced.

[0116] To further visually evaluate the reconstructed interferogram, Figure 8 the two profiles in Figure 9 (a) were analyzed. shows the profiles of the original and reconstructed interferograms of the horizontal and vertical black solid lines. The profile line of the existing deviation-corrected coherence coefficient contains more outliers. Obviously, the coherence coefficient corrected by the proposed method

[0117] produces a smoother profile with fewer outliers. The better performance demonstrates the effectiveness of the proposed method, which is mainly attributed to the effectiveness of the coherence matrix deviation correction.

[0118]

[0119] where, represents the interferometric phase of the coherence matrix . θ i and θ j are the phases reconstructed by the phase estimation algorithm. Figure 10 represents the average value of the temporal coherence coefficients for different numbers of SHPs. It can be found that the temporal coherence coefficient of the existing deviation-corrected coherence coefficient is higher than that of the sampling coherence coefficient and the improvement is most obvious when the number of samples is between 40 and 55. When the number of samples is low, due to the obvious coherence coefficient deviation, the improvement of the temporal coherence coefficient is insufficient. When using the corrected coherence coefficient proposed in the present invention the highest temporal coherence coefficient was obtained, showing the best performance in the reconstructed phase.

[0120] The present invention proposes an improved temporal InSAR distributed target phase estimation method. In particular, the new estimator can simultaneously incorporate the information of the interferometric coherence coefficient and the multi-look number, and achieve significant bias correction for each pixel of the image and each element of the coherence matrix. Then, the bias-corrected coherence matrix is combined with EMI to obtain an optimized phase sequence. Since the bias of the coherence matrix is reduced, this method obtains a phase sequence with high-precision reconstruction, especially in the case where the number of homogeneous pixels is small. Experiments on simulated data and Sentinel-1 real data verify the effectiveness of this method.

[0121] The above-described embodiments are only preferred specific embodiments of the present invention, and the protection scope of the present invention is not limited thereto. Any simple changes or equivalent replacements of the technical solutions that can be obviously obtained by those skilled in the art within the technical scope disclosed by the present invention shall fall within the protection scope of the present invention.

Claims

1. A method for optimizing the phase of distributed targets in temporal InSAR, characterized in that, It includes the following steps: Obtain N scenes of InSAR images and construct an InSAR dataset; Select a homogeneous pixel set Ω for each pixel p of the image through the CMP algorithm; Estimate the complex covariance matrix using a complex covariance matrix estimator and a selected set of homogeneous pixels Ω and the complex covariance matrix Take the modulus to obtain the coherence amplitude matrix Pre-estimate the empirical coherence coefficient in advance based on the information of the InSAR dataset The empirical coherence coefficient Is pre-estimated as: Among them, is thermal decoherence, where the signal-to-noise ratio SNR depends on system parameters and is set to a constant 12 dB; is geometric decorrelation, where B ⊥ and B ⊥max are the vertical baseline and the critical baseline, respectively; is temporal decoherence, where B T and B Tmax are the time baseline and the decoherence rate, respectively; Through the empirical coherence coefficient and the homogeneous pixel set Ω to determine the empirical constant s; The coherence coefficient bias corrector is obtained by combining the empirical constant s with the natural logarithm moments of the pixels The natural logarithm moment of the pixel is approximated as: where s is an empirical constant and L is the number of homogeneous pixels; The coherence coefficient corrector is: Among them, Through a coherence coefficient deviator Perform iterative correction on the coherence amplitude matrix with the homogeneous pixel set Ω ; The coherent amplitude matrix after deviation correction is input into the EMI phase evaluator to obtain the finally optimized phase sequence.

2. The method for optimizing the phase of distributed targets in temporal InSAR according to claim 1, characterized in that, The formula for the complex covariance matrix estimator is: where L is the number of homogeneous pixels, s(k) represents a normalized complex vector, and H represents conjugate transpose.

3. The method for optimizing the phase of distributed targets in temporal InSAR according to claim 1, characterized in that, The empirical constant s is estimated by the following formula: Among them, rounds the result to the largest integer not greater than it.

4. The method for optimizing the phase of distributed targets in temporal InSAR according to claim 1, characterized in that, The coherent amplitude matrix after deviation correction The specific formula for inputting into the EMI phase evaluator is as follows: Among them, replace with the corrected coherent amplitude matrix θ = [0, θ2,..., θ N T is the estimated consistent time series phase, the symbol represents the Hadamard product and the Lagrange multiplier method is used to solve the one with the minimum eigenvalue ​ 5. The method for optimizing the phase of distributed targets in temporal InSAR according to claim 1, characterized in that, EMI performs matrix inversion operations. For the case where the coherence amplitude matrix is not positive definite, the pseudo-inverse is used for processing.

Images