Abstract
In this article, we propose a method to estimate the synthetic aperture radar interferometry (InSAR) interferometric phase based on the model of correlation weight joint pixel by using the joint subspace projection technique. In the method, the correlation weight joint data vector is constructed and the data vector can make the noise subspace dimension of the corresponding weight covariance matrix which is not affected by the coregistration error, thus avoiding the trouble of calculating the noise subspace dimension before estimating the InSAR interferometric phase. The method takes advantage of the coherence information of neighboring pixel pairs to autocoregister the SAR images and employs the projection of the joint signal subspace onto the corresponding joint noise subspace to estimate the terrain interferometric phase. The method can autocoregister the SAR images and reduce the interferometric phase noise simultaneously. Theoretical analysis and computer simulation results show that the method can provide accurate estimate of the interferometric phase (interferogram) even when the coregistration error reaches one pixel. The effectiveness of the method is verified via simulated data and real data.
Keywords:
SAR interferometry (InSAR); Joint subspace projection; Interferometric phase; Noise subspace; Correlation weight joint data vector1. Introduction
Synthetic aperture radar interferometry (InSAR) is an important remote sensing technique to retrieve the terrain digital elevation model [1,2]. Image coregistration [15], InSAR interferometric phase estimation (or noise filtering), and interferometric phase unwrapping [69] are key processing procedures of InSAR. It is well known that the performance of interferometric phase estimation suffers seriously from the inaccuracy of the image coregistration. Almost all the conventional InSAR interferometric phase estimation methods are based on interferogram filtering, such as pivoting mean filtering [10], pivoting median filtering [11], adaptive phase noise filtering [12]. and adaptive contoured window filtering [13]. The problem here is that when the quality of an interferogram is very poor due to a large coregistration error, it is very difficult for these methods to retrieve the true terrain interferometric phases. In fact, the interferometric phases we obtained are random quantities with their variances being inversely proportional to the correlation coefficients between the corresponding pixel pairs of the two coregistered SAR images [2]. Therefore, the terrain interferometric phases should be estimated statistically.
In the previous studies [14,15], we proposed a joint subspace projection method to estimate InSAR interferometric phase in the presence of large coregistration errors. However, the noise subspace dimension of the covariance matrix changes with the coregistration error [14]. For accurately estimating the InSAR interferometric phase, the noise subspace dimension of the covariance matrix must be known, and the performance of the method degrades when the noise subspace dimension is not estimated correctly [16]. In this article, we propose a new method based on a correlation weight joint subspace projection to estimate the terrain interferometric phase accurately in the presence of large coregistration errors. In this method, the benefit from the correlation weight joint data vector is that the noise subspace dimension of the weight covariance matrix is not affected by the coregistration error [i.e., the noise subspace dimension of the corresponding covariance matrix with the coregistration error μ(0 < μ ≤ 1) pixel is the same as that of the covariance matrix with accurate coregistration]. The key processing procedures of the approach are summarized as follows: after the coarse coregistration, the correlation weight joint data vector is constructed which can be used to estimate the corresponding weight covariance matrix. The noise subspace is obtained from the eigendecomposition of the estimated covariance matrix and the signal subspace is spanned by the vectors that are obtained by the Hadamard product of the principal eigenvectors (i.e., the signal eigenvectors) of the weight correlation function matrix (approximated by the magnitude of the weight covariance matrix) and the steering vector. The terrain interferometric phase estimation is then performed by the projection of the signal subspace onto the corresponding noise subspace, where the optimum estimation corresponds to minimizing the projection. For a pair of SAR images that are not coregistered accurately, the method can autocoregister them and accurately estimate the corresponding terrain interferometric phase.
The article is arranged as follows: the statistical model of the correlation weight joint data vector is given in Section 2; in Section 3, the processing steps of the proposed method are presented in details; the performance of the method is investigated with the simulated and the real data in Section 4; the conclusions are given in Section 5.
2. Statistical model of the correlation weight joint data vector
When the SAR images are accurately coregistered and the interferometric phases are flattened with a zeroheight reference plane surface, the complex data vector s(i) of a pixel pair i (corresponding to the same ground area) of the coregistered SAR images can be formulated as [16,17]:
where
denotes the spatial steering vector (i.e., the array steering vector) of the pixel i, where s_{1} and s_{2} are the complex SAR images, superscript T denotes the vector transpose operation, φ_{i} is the terrain interferometric phase to be estimated, ⊙ denotes the Hadamard product, x(i) is the complex magnitude vector (i.e., complex reflectivity vector of scene received by the satellites) of pixel i, and n(i) is the additive noise term. The complex data vector s (i) can be modeled as a joint complex circular Gaussian random vector [1,2] with zeromean and the corresponding covariance matrix C_{s}(i) is given by
where
is called the correlation coefficient matrix, I is a 2 × 2 identity matrix, r_{mn}(ii) (0 ≤ r_{mn}(ii) ≤ 1, n = 1,2, m = 1,2) are the correlation coefficients between the satellites m and n, E{} denotes the statistical expectation, superscript H denotes vector conjugatetranspose, σ_{s}^{2}(i) is the backscatter power of the pixel i and σ_{n}^{2} is the noise power.
If the SAR images are accurately coregistered and the crosscorrelation coefficients (i.e., the nondiagonal elements) of R_{s}(i) are high enough, the rank of the correlation coefficient matrix R_{s}(i) becomes 1, and the number of the principal eigenvalue of R_{s}(i) is one [14].
In this article, an estimation method for InSAR interferometric phase based on correlation weight joint subspace projection is proposed, which can estimate the terrain interferometric phase accurately in the presence of large coregistration errors. Benefitting from the correlation weight joint data vector, the method does not need to calculate the noise subspace dimension before estimating the InSAR interferometric phase.
When the azimuth coregistration error is ρ(0 < ρ < 1) pixel and its direction is upwards (i.e., the pixel of the image from the second satellite is shifted upwards compared to the pixel in the first satellite image), the formulation of the correlation weight joint data vector [18]si(iw) is shown in Figure 1, where circles represents SAR image pixels and i denotes the desired pixel pair whose interferometric phase is to be estimated,^{a} the approach to determine the data vector is presented in literature [18].
where
Figure 1. Formulation of the correlation weight joint data vector.
The corresponding weight covariance matrix can be given by
where R_{si}(i,w) is called the weight correlation function matrix of the pixel pair i.
In the following, the characteristics of the signal and the noise subspaces of the weight covariance matrix C_{si}(i,w) for different coregistration errors are discussed.
When the example of formula (1) is used to build the data vector, the correlation coefficient of the corresponding item of the data vector decreases with the increasing coregistration error. On the contrary, the correlation coefficient is not impacted by the coregistration error when the correlation weight data vector is used. That means the correlation coefficient [i.e., the crosscorrelation coefficients of R_{s}(a,w_{a}) (a = i – 1, i, i + 3, i + 4)] of the corresponding item of the correlation weight data vector with the coregistration error μ(0 < μ ≤ 1) pixel is as high as that of the corresponding item of the data vector with accurate coregistration. Figure 2 shows the correlation coefficient for various data vector versus the coregistration error to further demonstrate the above conclusions.
Figure 2. Correlation coefficient versus the coregistration error.
From the literature [14,19], we know that the rank of R_{s}(aw_{a})(a = i – 1, i, i + 3, i + 4) becomes 1 when the crosscorrelation coefficients (i.e., the nondiagonal elements) of R_{s}(aw_{a})(a = i – 1, i, i + 3, i + 4) are high enough. Therefore, the number of the principal eigenvalues of R_{si}(iw) is 4 according to (9), and the number of the principal eigenvalues of C_{si}(iw) is also equal to 4 [14]. That is the dimensions of the signal subspace and the noise subspace are both 4 in the presence of large coregistration errors. Figure 3 shows the eigenspectra of the weight covariance matrix for different coregistration errors to further demonstrate the above conclusions.
Figure 3. Eigenspectra of the weight covariance matrix for accurate coregistration, coregistration errors of 0.5 and 1 pixels, respectively.
In this case, the eigendecomposition of the weight covariance matrix C_{si}(iw) is as follows:
where β_{csi}^{(k)}(k = 1, 2, 3, 4) is the eigenvector corresponding to the principal eigenvalue λ_{csi}^{(k)}(k = 1, 2, 3, 4) of C_{si}(iw), β_{rsi}^{(k)}(k = 1, 2, 3, 4) is the eigenvector corresponding to the principal eigenvalue λ_{rsi}^{(k)}(k = 1, 2, 3, 4) of R_{si}(iw), σ_{n}^{2} and β_{nsi}^{(l)}(k = 1, 2, 3, 4) are the noise eigenvalue and the corresponding eigenvector of C_{si}(iw). From (12) we can note that ai(φ_{i}) ⊙ β_{rsi}^{(k)}(k = 1, 2, 3, 4) is in the signal subspace of C_{si}(iw), β_{nsi}^{(l)}(k = 1, 2, 3, 4) is in the noise subspace of C_{si}(iw), and the signal subspace are orthogonal to the noise subspace [14].
From the results derived above, we can see that the new formulation of the correlation weight joint data vector proposed in this article has the advantage that the noise subspace dimension of the corresponding weight covariance matrix is independent of the coregistration error. That is to say, the noise subspace dimension of the corresponding covariance matrix with the coregistration error μ(0 < μ ≤ 1) pixel is the same as that of the accurate covariance matrix (in this article, the noise subspace dimension of the corresponding covariance matrix with the accurate coregistration is 4). Therefore, it is not required to calculate the noise subspace dimension before estimating the InSAR interferometric phase, thus the introduced trouble will be avoided.
3. Processing procedures
In this section, we give the detailed steps for this proposed method.
Step 1. Coregister SAR images.
The SAR images are coarsely coregistered using the crosscorrelation information of the SAR image intensity or other strategies [1,2] after SAR imaging of the echoes acquired by each satellite.
Remark 1: The required image coregistration accuracy for the proposed method can be 1 pixel, which is much lower than the required accuracy (from 1/10 to 1/100 pixel) for conventional methods. The low coregistration accuracy requirement can greatly facilitate coregistration processing.
Step 2. Estimate the covariance matrix.
An example to construct the formulation of the correlation weight joint data vector js(i,w) is shown in Figure 4 when the coregistration error is μ(0 ≤ μ ≤ 1) pixels and its direction is unknown.
where
Figure 4. Formulation of the correlation weight data vector when the coregistration error is unknown.
The corresponding weight covariance matrix can be given by
where R_{js}(iw) is called the weight correlation function matrix of the pixel pair i. Under the assumption that the neighboring pixels have the identical terrain height and the complex reflectivity is independent from pixeltopixel [16], the weight covariance matrix C_{js}(iw) given by (16) can be estimated by its sample covariance matrix , i.e.,
where 2 K + 1 is the number of i.i.d. samples from the neighboring pixel pairs.
Remark 2: According to the Reed–Mallett–Brennan rule [20], the effective number of looks (i.e., the number of i.i.d. samples) that 2 K + 1 ≥ 2 M – 1 would make the estimation loss within 3 dB if the dimensions of the covariance matrix are M × M. The detailed analysis of the effective number of looks is presented in [21].
It is easy to obtain enough i.i.d. samples for locally flat terrains. However, an imaging terrain in practice cannot be relied upon to be so flat that the adjacent pixels have the identical terrain height. If the local terrain slope is available in advance or can be estimated [15], the steering vector (i.e., the interferometric phase) variation due to the different terrain height from pixeltopixel can be compensated, which greatly enlarges the size of the sample window.
Step 3. Subspace estimation by eigendecomposing.
The estimated covariance matrix of the dimensions 8 × 8 can be eigendecomposed into
where K (in this article, K is equal to 4) is the number of the principal eigenvalues of, , eigenvectors corresponding to the smaller eigenvalues span the noise subspace, i.e.,
whereas the larger eigenvectors corresponding to the principal eigenvalues span the signal subspace, i.e.,
The noise power is often estimated by [14]
The joint correlation function matrix can be approximated as the amplitude (i.e., the absolute value) of the estimated covariance matrix [14], i.e.,
By eigendecomposing , we obtain K principal eigenvectors . As shown by (12), the same signal subspace spanned by the principal eigenvectors of can be spanned by the Hadamard product vectors , i.e.,
Step 4. Projection of signal subspace onto noise subspace.
As mentioned above, the signal subspace is orthogonal to the noise subspace, which is used to estimate the interferometric phase φ_{i}.
Definition of cost function:
where
The minimization of J can provide the optimum estimate of the interferometric phase φ_{i}, i.e., .
Remark 3: The computational burden will be high if the minimization of J is obtained via search of φ_{i} in the principal phase interval [−π, +π]. To reduce the computational burden, a fast algorithm to compute the minimization of J is developed in Appendix, where the closedform solution to the estimate of φ_{i} is directly obtained by using the fast algorithm.
With the use of the above four steps, the terrain interferogram can be recovered after the pixel pairs of the SAR images are processed separately.
4. Performance investigation
In this section, we demonstrate the robustness of the method to coregistration errors by using two sets of simulated data and a real dataset.
We assume that there are two formationflying satellites in the cartwheel formation, and we select one orbit position for simulation, with an effective crosstrack baseline of 562.93 m, an orbit height of 750 km, and an incidence angle of 50°. We use a twodimensional Hann window to simulate the terrain and use the statistical model to generate the complex SAR image pairs [22]. The signaltonoise ratio of the SAR images is 16 dB.
Here, the number of the samples to estimate the covariance matrix is 7 (in range) × 7 (in azimuth) = 49. The variation of the standard deviation of the interferometric phase with the increasing coregistration error was computed by means of 1,000 Monte Carlo simulations.
Figures 5, 6, and 7 compare the simulation results for various techniques and coregistration errors. Comparing Figures 5, 6, and 7, we can observe that the large coregistration error heavily affects the interferograms obtained by pivoting mean filtering, pivoting median filtering, adaptive phase noise filtering, and adaptive contoured window filtering. On the contrary, the large coregistration error has almost no effect on the interferograms obtained by the proposed method. We can see that the proposed method in this article is robust to large coregistration errors (up to 1 pixel).
Figure 5. Accurate coregistration. Interferograms obtained by (a) pivoting mean filtering, (b) pivoting median filtering, (c) adaptive phase noise filtering, (d) adaptive contoured window filtering, and (e) the proposed method.
Figure 6. Image coregistration error of 0.5 pixels. Interferograms obtained by (a) pivoting mean filtering, (b) pivoting median filtering, (c) adaptive phase noise filtering, (d) adaptive contoured window filtering, and (e) the proposed method.
Figure 7. Image coregistration error of 1 pixel. Interferograms obtained by (a) pivoting mean filtering, (b) pivoting median filtering, (c) adaptive phase noise filtering, (d) adaptive contoured window filtering, and (e) the proposed method.
Now we demonstrate the robustness of the method to local misregistration [16] (note that by “local misregistration” we mean that the coregistration error is not same for every image pixel) using the simulated data above.
Figure 8 shows the interferograms obtained by various techniques for the local misregistration. From the simulation results, we can observe that the local misregistration heavily affects the interferograms obtained by pivoting mean filtering, pivoting median filtering, adaptive phase noise filtering, and adaptive contoured window filtering. However, the proposed method is robust to the increasing coregistration error. In other words, our method can accurately estimate the corresponding terrain interferometric phase in the presence of the local misregistration.
Figure 8. Local misregistration. Interferograms obtained by (a) pivoting mean filtering, (b) pivoting median filtering, (c) adaptive phase noise filtering, (d) adaptive contoured window filtering, and (e) the proposed method.
Figure 9 shows the variation of the standard deviation of the interferometric phase with the increasing coregistration error. We can see that the performance of our method is as good as that of the joint subspace projection method [14]. However, our method can avoid the trouble of calculating the noise subspace dimension before estimating the InSAR interferometric phase.
Figure 9. Standard deviation of the interferometric phase versus the coregistration error.
The second simulated data are the interferometric data pair of Mount Etna (the data are produced based on the SIRC/XSAR data acquired at Xband).^{b}
Figure 10 shows the interferograms generated from the Mount Etna data. Figure 10a is the interferogram obtained by conventional processing [16] (i.e., directly computing the interferometric phase pixelbypixel), and Figure 10b is the interferogram obtained by the approach proposed in this article.
Figure 10. Interferograms obtained by (a) conventional processing, and (b) proposed method for Mount Etna data.
Figure 11 confirms the effectiveness of the proposed method in processing of the ERS1/ERS2 (European Remote Sensing 1 and 2 tandem satellites, ERS1 orbit = 32585, ERS2 orbit = 12912, frame = 2781, 19971008/09, Zhangbei, China) real data.
Figure 11. The interferograms obtained by (a) the conventional processing, and (b) the proposed method for the ERS1/ERS2 real data.
Figure 11 shows the interferograms generated from the ERS1/ERS2 real data. Figure 11a is the interferogram obtained by conventional processing, and Figure 11b is the interferogram obtained by the approach proposed in this article.
5. Conclusions
We have proposed a new method to estimate the terrain interferometric phases from the InSAR image pair. Benefiting from the new formulation of correlation weight joint data vector, the method does not need to calculate the noise subspace dimension before estimating the InSAR interferometric phase, thus the introduced trouble will be avoided. The method is based on the projection of the joint signal subspace onto the corresponding joint noise subspace, and takes advantage of the coherence information of the neighboring pixel pairs to autocoregister the SAR images, where the phase noise is reduced simultaneously. A fast algorithm is developed to implement the method, which can significantly reduce the computational burden. The effectiveness of the method is verified via simulated data and real data.
Appendix
Fast algorithm to compute the optimum interferometric phase estimate
If U, V, and W are arbitrary complex column vectors, then [14]
Using Equation (26), we can rewrite the cost function of (24) as
Let . It can easily be proved that A(8 × 8) is a Hermitian matrix. Then (27) can be rewritten as
where
Let . It can easily be proved that B is a Hermitian matrix, i.e., b_{12}^{∗} = b_{21}, Let b_{12} = b_{12}e^{jμ}. Then,
Obviously, the minimization of J can be obtained for μ + ϕ_{i} = −π + 2kπ (k is an integer and μ = angle(b_{12})). Since –π ≤ μ < π and –π < ϕ_{i} < π, thus
Endnotes
^{a}As shown in Figure 1, the horizon is along range and the verticality is along azimuth.
^{b}Epsilon Nought, Radar Remote Sensing: http://epsilon.nought.de/ webcite.
Competing interests
The authors declare that they have no competing interests.
Acknowledgements
This paper was supported by the National Natural Science Foundations of China under grant 61071194, 60979002 and 61231017, by the Fund of Civil Aviation University of China under grant 2011kyE06, and the National Key Technology Research and Development Program of China (2011BAH24B12).
References

PA Rosen, S Hensley, IR Joughin, FK Li, SN Madsen, E Rodrìguez, RM Goldstein, Synthetic aperture radar interferometry. Proc. IEEE 88(3), 333–382 (2000)

R Bamler, P Hartl, Synthetic aperture radar interferometry. Inverse Probl. 14, R1–R54 (1998). Publisher Full Text

Q Lin, JF Vesecky, HA Zebker, Registration of Interferometric SAR Images. Proceedings of the International Geoscience and Remote Sensing Symposium’1992, Houston, Texas, USA, 1579–1581 (1992)

R Scheiber, A Moreira, Coregistration of interferometric SAR images using spectral diversity. IEEE Trans. Geosci. Remote Sens. 38(5), 2179–2191 (2000). Publisher Full Text

D Fernandes, G Waller, JR Moreira, Registration of SAR images using the chirp scaling algorithm. Proceedings of th International Geoscience and Remote Sensing Symposium’1996, Lincoln, NE, USA, 799–801 (1996)

X Wei, C Ian, A regiongrowing algorithm for InSAR phase unwrapping. IEEE Trans. Geosci. Remote Sens. 37(1), 124–134 (1999). Publisher Full Text

M Costantini, A novel phase unwrapping method based on network programming. IEEE Trans. Geosci. Remote Sens. 36(3), 813–831 (1998). Publisher Full Text

RM Goldstein, HA Zebker, CL Werner, Satellite radar interferometry: twodimensional phase unwrapping. Radio Sci. 23(4), 713–720 (1988). Publisher Full Text

MD Pritt, JS Shipman, Leastsquares twodimensional phase unwrapping using FFT’s. IEEE Trans. Geosci. Remote Sens. 32(3), 706–708 (1994). Publisher Full Text

PH Eichel, DC Ghiglia, CV Jakowatz Jr., PA Thompson, DE Wahl, Spotlight SAR interferometry for terrain elevation mapping and interferometric change detection. Sandia National Labs Technical Report, SAND93, 2539–2546 (1993)

R Lanari, G Fornaro, D Riccio, M Migliaccio, KP Papathanassiou, JR Moreira, M Schwabisch, L Dutra, G Puglisi, G Franceschetti, M Coltelli, Generation of digital elevation models by using SIRC/XSAR multifrequency twopass interferometry: the Etna case study. IEEE Trans. Geosci. Remote Sens. 34(5), 1097–1114 (1996). Publisher Full Text

JS Lee, KP Papathanassiou, TL Ainsworth, MR Grunes, A Reigber, A new technique for noise filtering of SAR interferometric phase images. IEEE Trans. Geosci. Remote Sens. 36(5), 1456–1465 (1998). Publisher Full Text

Y Qifeng, X Yang, F Sihua, X Liu, X Sun, An adaptive contoured window filter for interferometric synthetic aperture radar. IEEE Geosci. Remote Sens. Lett. 4(1), 23–26 (2007)

L Zhenfang, B Zheng, L Hai, L Guisheng, Image autocoregistration and InSAR interferogram estimation using joint subspace projection. IEEE Trans. Geosci. Remote Sens. 44(2), 288–297 (2006)

H Li, Z Li, G Liao, Z Bao, An estimation method for InSAR interferometric phase combined with image autocoregistration. Sci. China, Ser. F. 49(3), 386–396 (2006). Publisher Full Text

H Li, G Liao, An estimation method for InSAR interferometric phase based on MMSE criterion. IEEE Trans. Geosci. Remote Sens. 48(3), 1457–1469 (2010)

F Lombardini, M Montanari, F Gini, Reflectivity estimation for multibaseline interferometric radar imaging of layover extended sources. IEEE Trans. Signal Process. 51(6), 1508–1519 (2003). Publisher Full Text

L Hai, L Guisheng, A phase unwrapping method for multibaseline InSAR systems based on correlation coefficient weight data vector. Prog. Nat. Sci. (in Chinese) 18(3), 313–322 (2008)

L Hai, W Renbiao, InSAR interferometric phase estimation based on correlation weight subspace projection. 2011 IEEE CIE International Conference on Radar, Chengdu, China, 398–401 (2011)

IS Reed, JD Mallett, LE Brennan, Rapid convergence rate in adaptive arrays. IEEE Trans. Aerosp. Electron. Syst. 10, 853–863 (1974)

CH Gierull, IC Sikaneta, Estimating the effective number of looks in interferometric SAR data. IEEE Trans. Geosci. Remote Sens. 40(8), 1733–1742 (2002). Publisher Full Text

F Lombardini, Absolute phase retrieval in a threeelement synthetic aperture radar interferometer. Proceedings of the 1996 CIE International Conference on Radar, Beijing, China, 309–312 (1996)