**Robust Interferometric Phase Estimation in InSAR via Joint Subspace Projection**

Hai Li and Renbiao Wu

*Tianjin Key Lab for Advanced Signal Processing, Civil Aviation University of China, Tianjin, P.R. China*

### **1. Introduction**

16 Will-be-set-by-IN-TECH

110 Recent Interferometry Applications in Topography and Astronomy

Lee, J. S.; Papathanassiou, K. P.; Ainsworth, T. L.; Grunes, M. R. & Reigber, A. (1998). "A new

Natsuaki, R. & Hirose A. (2011). "SPEC Method – A fine co-registration method for SAR

Pritt, M. & Shipman, j. (1994). "Least-squares two-dimensional phase unwrapping using

Reigber, A. & Moreia, J. (1997). "Phase unwrapping by fusion of local and global methods,"

Suksmono, A. B. & A. Hirose, A. (2002). "Adaptive noise reduction of InSAR images based on

Suksmono, A. B. and Hirose, A. (2006). "Progressive transform-based phase unwrapping

Tobita, M.; Fujiwara, S.; Murakami, M.; Nakagawa, H. & Rosen, P. A. (1999). "Accurate offset

Yamaki, R. & Hirose, A. (2007). "Singularity-Spreading Phase Unwrapping," *IEEE Transactions*

Yamaki, R. & Hirose, A. (2009). "Singular Unit Restoration in Interferograms Based on

*Geoscience and Remote Sensing Letters*, Vol.6, No.1, pp.18-22, January 2009.

*Geoscience and Remote Sensing*, vol. 36, pp. 1456-1464, Sept. 1998.

28-37, January 2011.

May 1994.

August 1997.

3, pp. 929-936, Mar. 2006.

*the Geodetic Society of Japan*, 45, 297-314, 1999.

*on Geoscience and Remote Sensing*, Vol.45, No.10, Oct. 2007.

2002.

technique for noise filtering of sar interferometric phase images,"*IEEE Transactions on*

interferometry," *IEEE Transactions Geoscience and Remote Sensing*, Vol.49, No.1, pp.

FFT's," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 32, no. 3, pp. 706-708,

in *International Geoscience and Remote Sensing Symposium (IGARSS) 1997*, pp. 869-871,

a complex-valued MRF model and its application to phase unwrapping problem," *IEEE Transactions on Geoscience and Remote Sensing*, Vol.40, No.3, pp.699-709, March

utilizing a recursive structure,"*IEICE Transactions on Communications*, vol. E89-B, no.

estimation between two SLC images for SAR interferometry" (in Japanese), *Journal of*

Complex-valued Markov Random Field Model for Phase Unwrapping," *IEEE*

Synthetic aperture radar interferometry (InSAR) is an important remote sensing technique to retrieve the terrain digital elevation model (DEM)[1][2]. Image coregistration, InSAR interferometric phase estimation (or noise filtering) and interferometric phase unwrapping[3][4][5][6] are three key processing procedures of InSAR. It is well known that the performance of interferometric phase estimation suffers seriously from poor image coregistration.

Image coregistration is an important preprocessing operation that aligns the pixels of one image to the corresponding pixels of another image. A review of recent as well as classic image registration methods can be found in Ref.[7]. Mutual information used for the registration of remote sensing imagery are presented in the literature[8]. In literature[9], the feature-based registration methods are presented. A new direct Fourier-transform-based algorithm for subpixel registration is proposed in Ref.[10]. In Ref.[11], a geometrical approach for image registration of SAR images is proposed, and the algorithm has been tested on several real data. A image registration method based on isolated point scatterers is proposed in Ref.[12].

Almost all the conventional InSAR interferometric phase estimation methods are based on interferogram filtering[13][14][15][16][17][18], such as pivoting mean filtering[13], pivoting median filtering[14], adaptive phase noise filtering[15], and adaptive contoured window filtering[18]. However, 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 are random in nature 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 this chapter, the interferometric phase estimation method based on subspace projection and its modified version were proposed. Theoretical analysis and computer simulation results show that the methods can provide accurate estimation of the terrain interferometric phase (interferogram) even if the coregistration error reaches one pixel. The remainder of this chapter is organized as follows. Section 2 presents the signal model of a single pixel pair

Robust Interferometric Phase Estimation in InSAR via Joint Subspace Projection 113

signal eigenvector) of ( )*i* **Rs** , respectively, and **β***n* is the noise eigenvector corresponding to

1 (( ) ) (( ) ) *H H i rs n n i rs J* **a**

The eigenvalues of *Cs*(*i*) result in a low dispersion due to the increase of the coregistration error. Actually, it is induced by the increase of noise eigenvalue, i.e., the signal component spreads into the noise space. Moreover, the 2-dimensional space will be fully occupied by signal component with a worse coregistration error of one pixel. At this instant, the noise eigenvalue is almost equal to the bigger one. The noise subspace dimension becomes zero. The degree of the signal component spreading to the noise subspace is smaller, the estimation of the InSAR interferometric phase is better(the subspace projection technique is used to estimate the InSAR interferometric phase ). On the contrary, the degree of the signal component spreading to the noise subspace is larger, the estimation of the InSAR interferometric phase is worse. The conclusions obtained from the preceding analysis are briefly summarized as follows: the degree of the signal component spreading to the noise subspace becomes larger and larger as the coregistration error increases when the example of formula (1) is used to build the data vector. In other words, the degree of dispersion is heavily impacted by the coregistration error. The simulation result shown in Fig.1

Fig. 1. Eigenspectra of the covariance matrix for accurate coregistration, coregistration errors

The cost function given by (5) can be used to estimate the InSAR interferometric phase when the SAR images are accurately coregistered. However, in the presence of coregistration

The minimization of 1*J* can provide the optimum estimate of the interferometric phase

the insignificant eigenvalue of ( ) *s* **C** *i* . From (4) we can note that (( ) )

*i* .

subspace, **β***n* is in the noise subspace, and (( ) )

to estimate the interferometric phase

demonstrates the above conclusions.

of 0.5 and 1 pixels, respectively.

The definition of cost function is,

*rs* and **β***rs* are the principal eigenvalue and the corresponding eigenvector (i.e., the

*i rs* **a β** are orthogonal to **β***n* , which is used

**β ββ a β** (5)

*i rs* **a β** is in the signal

*i* .

where

and the problem formulation. In Section 3, we discuss the interferometric phase estimation method based on subspace projection. In Section 4, the modified interferometric phase estimation method via subspace projection is presented in details. Finally, numerical and experemental results are presented in Section 5. Section 6 concludes the whole chapter.

### **2. Data model and problem formulation**

Assuming that the SAR images are accurately coregistered and the interferometric phases are flattened with a zero-height reference plane surface. The complex data vector, denoted as **s**( )*i* , of a pixel pair *i* (corresponding to the same ground area) of the coregistered SAR images can be formulated as follows[19],

$$\mathbf{s}(i) = \begin{bmatrix} \mathbf{s}\_1(i), \mathbf{s}\_2(i) \end{bmatrix}^\mathsf{T} = \mathbf{a}(\boldsymbol{\varphi}\_i) \odot \begin{bmatrix} \mathbf{x}\_1(i), \mathbf{x}\_2(i) \end{bmatrix}^\mathsf{T} + \mathbf{n}(i) = \mathbf{a}(\boldsymbol{\varphi}\_i) \odot \mathbf{x}(i) + \mathbf{n}(i) \tag{1}$$

where ( ) [1, ]*i j T i e* **a** is the spatial steering vector (i.e., the array steering vector) of the pixel *i* , 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 the 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 with zero-mean and the corresponding covariance matrix ( ) *s* **C** *i* is given by

$$\begin{aligned} \mathbf{C\_s}(i) &= E\{\mathbf{s}(i)\mathbf{s}^H(i)\} \\ &= \mathbf{a}(\boldsymbol{\varphi\_i})\mathbf{a}^H(\boldsymbol{\varphi\_i}) \odot E\{\mathbf{x}(i)\mathbf{x}^H(i)\} + \sigma\_n^2 \mathbf{I} \\ &= \sigma\_s^2(i)\mathbf{a}(\boldsymbol{\varphi\_i})\mathbf{a}^H(\boldsymbol{\varphi\_i}) \odot \mathbf{R\_s}(i) + \sigma\_n^2 \mathbf{I} \end{aligned} \tag{2}$$

$$\mathbf{R}\_s(i) = \begin{bmatrix} r\_{11}(i) \ \ r\_{12}(i) \\ r\_{21}(i) \ r\_{22}(i) \end{bmatrix} \tag{3}$$

where ( )*i* **Rs** is called the correlation coefficient matrix, **I** is a 22 identity matrix, *r i mn*( ) ( 0 () 1 *mn r i* , *n* 1,2 , and 1,2 *m* ) are the correlation coefficients between the satellites *m* and *n* , {} *E* denotes the statistical expectation, superscript *H* denotes vector conjugatetranspose, 2 ( ) *s i* is the echo power of the pixel *i* and 2 *n* is the noise power. In order to simplify the mathematical expressions, the denotation *i* (denoting the pixel) in the right side of the following expressions is omitted. In practice, the statistical covariance matrix of (2) can be adaptively estimated using the sample covariance matrix.

If the SAR images are accurately coregistered and the cross-correlation coefficients (i.e., the nondiagonal elements) of ( )*i* **Rs** are large enough, the number of the principal eigenvalues of the covariance matrix ( ) *s* **C** *i* is one; i.e., the dimensions of the signal subspace and the noise subspace are both one (in the absence of the layover). In this case, the eigendecomposition of the covariance matrix ( ) *s* **C** *i* is as follows:

$$\mathbf{C\_s}(i) = (\mathbb{A}\_{rs} + \sigma\_n^2)(\mathbf{a}(\rho\_i) \odot \mathbf{\mathcal{B}\_{rs}})(\mathbf{a}(\rho\_i) \odot \mathbf{\mathcal{B}\_{rs}})^H + \sigma\_n^2 \mathbf{\mathcal{B}\_n} \mathbf{\mathcal{B}\_n}^H \tag{4}$$

where *rs* and **β***rs* are the principal eigenvalue and the corresponding eigenvector (i.e., the signal eigenvector) of ( )*i* **Rs** , respectively, and **β***n* is the noise eigenvector corresponding to the insignificant eigenvalue of ( ) *s* **C** *i* . From (4) we can note that (( ) ) *i rs* **a β** is in the signal subspace, **β***n* is in the noise subspace, and (( ) ) *i rs* **a β** are orthogonal to **β***n* , which is used to estimate the interferometric phase *i* .

The definition of cost function is,

112 Recent Interferometry Applications in Topography and Astronomy

and the problem formulation. In Section 3, we discuss the interferometric phase estimation method based on subspace projection. In Section 4, the modified interferometric phase estimation method via subspace projection is presented in details. Finally, numerical and experemental results are presented in Section 5. Section 6 concludes the whole chapter.

Assuming that the SAR images are accurately coregistered and the interferometric phases are flattened with a zero-height reference plane surface. The complex data vector, denoted as **s**( )*i* , of a pixel pair *i* (corresponding to the same ground area) of the coregistered SAR

1 2 1 2 ( ) [ ( ), ( )] [ ( ), ( )] ( ) ( ) ( ) *i i i s i s i ( ) x i x i (i)*

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 the 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 with zero-mean and the

*i i* **s a na xn** (1)

is the spatial steering vector (i.e., the array steering vector) of the

2 2

**a a xx I aa R I**

*Ei i*

( ) ( ) { ( ) ( )} () ( ) ( ) ()

*i i*

11 12 21 22 ( ), ( ) ( ) ( ), ( ) *s riri*

*r ir i*

where ( )*i* **Rs** is called the correlation coefficient matrix, **I** is a 22 identity matrix, *r i mn*( ) ( 0 () 1 *mn r i* , *n* 1,2 , and 1,2 *m* ) are the correlation coefficients between the satellites *m* and *n* , {} *E* denotes the statistical expectation, superscript *H* denotes vector conjugate-

simplify the mathematical expressions, the denotation *i* (denoting the pixel) in the right side of the following expressions is omitted. In practice, the statistical covariance matrix of

If the SAR images are accurately coregistered and the cross-correlation coefficients (i.e., the nondiagonal elements) of ( )*i* **Rs** are large enough, the number of the principal eigenvalues of the covariance matrix ( ) *s* **C** *i* is one; i.e., the dimensions of the signal subspace and the noise subspace are both one (in the absence of the layover). In this case, the eigen-

> 2 2 ( ) ( )( ( ) )( ( ) )*H H rs n i rs i rs n n n* **C a s** *i*

*H s i is n*

*H H i i n*

2

**R** (3)

 **β a β ββ** (4)

*i* is the terrain

(2)

*n* is the noise power. In order to

pixel *i* , superscript *T* denotes the vector transpose operation,

( ) { ( ) ( )}

*i*

*i* is the echo power of the pixel *i* and 2

(2) can be adaptively estimated using the sample covariance matrix.

decomposition of the covariance matrix ( ) *s* **C** *i* is as follows:

*i Ei i*

**C ss s**

*H*

**2. Data model and problem formulation**

images can be formulated as follows[19],

corresponding covariance matrix ( ) *s* **C** *i* is given by

where ( ) [1, ]*i j T i e*

transpose, 2

( ) *s*

 **a**

$$J\_1 = \left(\mathbf{a}(\varphi\_i) \odot \mathfrak{B}\_{rs}\right)^H \mathfrak{B}\_n \mathfrak{B}\_n^H \left(\mathbf{a}(\varphi\_i) \odot \mathfrak{B}\_{rs}\right) \tag{5}$$

The minimization of 1*J* can provide the optimum estimate of the interferometric phase *i* .

The eigenvalues of *Cs*(*i*) result in a low dispersion due to the increase of the coregistration error. Actually, it is induced by the increase of noise eigenvalue, i.e., the signal component spreads into the noise space. Moreover, the 2-dimensional space will be fully occupied by signal component with a worse coregistration error of one pixel. At this instant, the noise eigenvalue is almost equal to the bigger one. The noise subspace dimension becomes zero. The degree of the signal component spreading to the noise subspace is smaller, the estimation of the InSAR interferometric phase is better(the subspace projection technique is used to estimate the InSAR interferometric phase ). On the contrary, the degree of the signal component spreading to the noise subspace is larger, the estimation of the InSAR interferometric phase is worse. The conclusions obtained from the preceding analysis are briefly summarized as follows: the degree of the signal component spreading to the noise subspace becomes larger and larger as the coregistration error increases when the example of formula (1) is used to build the data vector. In other words, the degree of dispersion is heavily impacted by the coregistration error. The simulation result shown in Fig.1 demonstrates the above conclusions.

Fig. 1. Eigenspectra of the covariance matrix for accurate coregistration, coregistration errors of 0.5 and 1 pixels, respectively.

The cost function given by (5) can be used to estimate the InSAR interferometric phase when the SAR images are accurately coregistered. However, in the presence of coregistration

Robust Interferometric Phase Estimation in InSAR via Joint Subspace Projection 115

in (7), we have,

2

1 1

*nsi* **β** ( 1,2, ,18 *l K* ) are the noise eigenvalue and the corresponding eigenvectors of

If the SAR images are accurately coregistered, the structure form of the joint correlation

*K K*

  ( ) 2 2 ( ) ( ) () ()

2

*i i*

*s m* ( *mi i i* 4, 3, , 4 ) in the diagonal of the ( )*i* **Rsi** are the

**R**

**s**

(9)

*s*

0 ( 4) ( 4)

coherence coefficient matrix (given by (3)) and the echo power of the pixel pair *m* , respectively. We can notice from (9) that when the SAR images are accurately coregistered, only the elements in **Rs** ( ) *m* are nonzero, while all the other elements of ( )*i* **Rsi** (or ( )*i* **Csi** ) are zero (assuming the complex reflectivity is independent from pixel to pixel and neglecting the noise). In other words, ( )*i* **Rsi** is a block diagonal matrix. However, if the SAR images are not accurately coregistered, the nonzero elements in the submatrices **Rs** ( ) *m* are diffused to other non-diagonal element positions of ( )*i* **Rsi** , as shown in Fig.3(b) and Fig.3(c). The dimensions of the joint correlation function matrix ( )*i* **Rsi** are 18×18. Fig.3(a) is the structure of ( )*i* **Rsi** for accuracy coregistration; Fig.3(b) is the structure of ( )*i* **Rsi** for the coregistration error of 0.5 pixel; Fig.3(c) is the structure of ( )*i* **Rsi** for the coregistration error of 1 pixel. From Fig.3 we can see that when the SAR images are not accurately coregistered, the same imaged ground area can be coregistered to different pixel positions; thus its correlation information appears in other non-diagonal element position of ( )*i* **Rsi** . The gray

level denotes the magnitude of each element with the white strongest and black zero.

( 4) ( 4) 0 0

**s**

*k k k H l lH rsi ni i rsi rsi n nsi nsi*

**ai β ai β ββ**

*csi*

*nsi* **β** ( 1,2, ,18 *l K* ) are in the noise subspace of ( )*i* **Csi** , and

18

of ( )*i* **Csi** , ( ) *k*

*rsi*

*i rsi* **ai β** ( 1,2, , *k K* ) are in the signal

*csi* **β** ( *k K* 1,2, , ) are the

of ( )*i* **Rsi** . 2

*rsi* **β** ( 1,2, , *k K* )

*n* and

(8)

*i ii i* (18×1).

steering vector of the pixel *i* is denoted by ( ) [ ( ), ( ), , ( )] *T T TT* **ai a a a**

18 () () () 2 () ()

**ββ ββ**

*k k kH l lH csi csi csi n nsi nsi*

( )( ( ) )( ( ) )

*k l*

Substituting ( )

( )*l*

subspace of ( )*i* **Csi** , ( )*l*

where **Rs** ( ) *m* and 2

*i rsi* **ai β** are orthogonal to ( )*l*

function matrix ( )*i* **Rsi** are given by

( ) ( ) *k*

*i* **ai** for 43 4 ( , ,, ) **α**

*EVD*

() ( ) ( ) ()

*i i*

*H*

*ii i*

**C ai ai R I si si**

1 1

*k l*

where *K* is the number of the principal eigenvalues of ( )*i* **Csi** , ( ) *k*

*nsi* **β** .

2

*s*

2

*i i*

**R**

*s*

0 ( 3) ( 3)

**s**

*i i*

**R**

are the eigenvectors corresponding to the principal eigenvalues ( ) *k*

eigenvectors corresponding to the principal eigenvalues ( ) *k*

( )*i* **Csi** ,respectively. From (8) we can note that ( ) ( ) *k*

( )

*i*

**si**

**R**

( )

*ii n K K*

error, the cross-correlation coefficients are smaller than 1 and the noise eigenvalue becomes large. In the completely misregistered case, the rank of ( )*i* **Rs** becomes 2, which means that the noise subspace vanishes in the eigenspace of ( ) *s* **C** *i* . So the cost function given by (5) can not be used to estimate the InSAR interferometric phase in the presence of coregistration error.

### **3. Interferometric phase estimation via subspace projection**

Considering the difficulties in accurate coregistration, we use not only the corresponding pixel pair *i* of the coarsely coregistered SAR image pair (as given in (1)) but also the neighboring pixel pairs centered on the pixel pair *i* to jointly construct the data vector. An example of the construction method for the data vector is shown in Fig.2, where a circle represents a SAR image pixel and *i* denotes the centric pixel pair (i.e., the desired pixel pair whose interferometric phase is to be estimated). We call this extended data vector **si**( )*i* the joint data vector. The number of the neighboring pixel pairs to construct the joint data vector is 8 as shown in Fig.2.

Fig.2. A construction method for the joint data vector.

The joint data vector **si**( )*i* as shown in Fig.2 can be written as:

$$\mathbf{s}(i) = \left[\mathbf{s}(i-4)^T, \mathbf{s}(i-3)^T, \dots, \mathbf{s}(i)^T, \dots, \mathbf{s}(i+4)^T\right]^T \tag{6}$$

The corresponding joint covariance matrix is given by

$$\begin{split} \mathbf{C}\_{\text{si}}(i) &= E\{\mathbf{si}(i)\mathbf{s}\mathbf{i}^H(i)\} \\ &= \mathbf{a}(\varphi\_{i-4}, \varphi\_{i-3}, \dots, \varphi\_{i+4}) \mathbf{a}^H(\varphi\_{i-4}, \varphi\_{i-3}, \dots, \varphi\_{i+4}) \odot \mathbf{R}\_{\text{si}}(i) + \sigma\_n^2 \mathbf{I} \end{split} \tag{7}$$

where 43 4 4 3 4 ( , , , ) [ ( ), ( ), , ( )] *T T TT* **α** *ii i i i i* **aa a** and ( )*i* **Rsi** are referred to as the joint steering vector and the joint correlation function matrix of the pixel pair *i* , respectively. In fact, most of the nature terrain can be approximated by a local plane. After we estimate the local slopes and correct the neighboring pixels, we can assume that the neighboring pixels have the almost identical terrain height[20]. That is, the spatial steering vectors of the pixel pairs in **si**( )*i* are assumed to be identical, i.e., 4 ( ) *i* **a** = 3 ( ) *i* **a** = = 4 ( ) [1, ]*i j T i e* **a** and 43 4 ( , , , ) [ ( ), ( ), , ( )] [1, ,1, ] *i i T T TT j j T ii i i i i e e* **α aa a** . The simplified joint

error, the cross-correlation coefficients are smaller than 1 and the noise eigenvalue becomes large. In the completely misregistered case, the rank of ( )*i* **Rs** becomes 2, which means that the noise subspace vanishes in the eigenspace of ( ) *s* **C** *i* . So the cost function given by (5) can not be used to estimate the InSAR interferometric phase in the presence of coregistration

Considering the difficulties in accurate coregistration, we use not only the corresponding pixel pair *i* of the coarsely coregistered SAR image pair (as given in (1)) but also the neighboring pixel pairs centered on the pixel pair *i* to jointly construct the data vector. An example of the construction method for the data vector is shown in Fig.2, where a circle represents a SAR image pixel and *i* denotes the centric pixel pair (i.e., the desired pixel pair whose interferometric phase is to be estimated). We call this extended data vector **si**( )*i* the joint data vector. The number of the neighboring pixel pairs to construct the joint data vector

( ) [ ( 4) , ( 3) , , ( ) , , ( 4) ] *T TT T T* **si s s s s** *ii i i i* (6)

>

43 4 43 4

joint steering vector and the joint correlation function matrix of the pixel pair *i* , respectively. In fact, most of the nature terrain can be approximated by a local plane. After we estimate the local slopes and correct the neighboring pixels, we can assume that the neighboring pixels have the almost identical terrain height[20]. That is, the spatial steering vectors of the pixel

*H*

( , , , ) ( , , , ) ()

**aa a** . The simplified joint

*ii i ii i n*

 *ii i i i i* **aa a** and ( )*i* **Rsi** are referred to as the

> *i* **a** = 3 ( )

*i*

**α α R I**

2

*i e*

 **a**

(7)

and

*i* **a** = = 4 ( ) [1, ]*i j T*

**si**

**3. Interferometric phase estimation via subspace projection**

Fig.2. A construction method for the joint data vector.

The corresponding joint covariance matrix is given by

( ) { ( ) ( )}

*iEi i*

pairs in **si**( )*i* are assumed to be identical, i.e., 4 ( )

**C si si**

**si**

The joint data vector **si**( )*i* as shown in Fig.2 can be written as:

*H*

43 4 ( , , , ) [ ( ), ( ), , ( )] [1, ,1, ] *i i T T TT j j T ii i i i i e e*

**α**

where 43 4 4 3 4 ( , , , ) [ ( ), ( ), , ( )] *T T TT* **α**

>

> >

error.

is 8 as shown in Fig.2.

steering vector of the pixel *i* is denoted by ( ) [ ( ), ( ), , ( )] *T T TT* **ai a a a** *i ii i* (18×1). Substituting ( ) *i* **ai** for 43 4 ( , ,, ) **α** *ii i* in (7), we have,

$$\begin{split} \mathbf{C}\_{\mathbf{s}\mathbf{i}}(i) &= \mathbf{a}\mathbf{i}(\boldsymbol{\varphi}\_{i})\mathbf{a}^{\mathbf{i}\cdot\mathbf{l}}(\boldsymbol{\varphi}\_{i}) \odot \mathbf{R}\_{\mathbf{s}\mathbf{i}}(i) + \sigma\_{n}^{2}\mathbf{I} \\ &\quad \underline{\mathbf{E}VD} \sum\_{k=1}^{K} \lambda\_{\text{csl}}^{(k)} \mathfrak{P}\_{\text{csl}}^{(k)} \mathfrak{P}\_{\text{csl}}^{(k)H} + \sum\_{l=1}^{18-K} \sigma\_{n}^{2} \mathfrak{P}\_{\text{nsl}}^{(l)} \mathfrak{P}\_{\text{nsl}}^{(l)H} \\ &= \sum\_{k=1}^{K} (\lambda\_{\text{rsl}}^{(k)} + \sigma\_{n}^{2}) (\mathbf{a}\mathbf{i}(\boldsymbol{\varphi}\_{i}) \odot \mathfrak{P}\_{\text{rsl}}^{(k)}) (\mathbf{a}\mathbf{i}(\boldsymbol{\varphi}\_{i}) \odot \mathfrak{P}\_{\text{rsl}}^{(k)})^{H} + \sum\_{l=1}^{18-K} \sigma\_{n}^{2} \mathfrak{P}\_{\text{nsl}}^{(l)} \mathfrak{P}\_{\text{nsl}}^{(l)H} \end{split} \tag{8}$$

where *K* is the number of the principal eigenvalues of ( )*i* **Csi** , ( ) *k csi* **β** ( *k K* 1,2, , ) are the eigenvectors corresponding to the principal eigenvalues ( ) *k csi* of ( )*i* **Csi** , ( ) *k rsi* **β** ( 1,2, , *k K* ) are the eigenvectors corresponding to the principal eigenvalues ( ) *k rsi* of ( )*i* **Rsi** . 2 *n* and ( )*l nsi* **β** ( 1,2, ,18 *l K* ) are the noise eigenvalue and the corresponding eigenvectors of ( )*i* **Csi** ,respectively. From (8) we can note that ( ) ( ) *k i rsi* **ai β** ( 1,2, , *k K* ) are in the signal subspace of ( )*i* **Csi** , ( )*l nsi* **β** ( 1,2, ,18 *l K* ) are in the noise subspace of ( )*i* **Csi** , and ( ) ( ) *k i rsi* **ai β** are orthogonal to ( )*l nsi* **β** .

If the SAR images are accurately coregistered, the structure form of the joint correlation function matrix ( )*i* **Rsi** are given by

$$\mathbf{R}\_{\text{si}}(i) = \begin{bmatrix} \sigma\_s^2(i-4)\mathbf{R}\_s(i-4) & 0 & \cdots & 0\\ 0 & \sigma\_s^2(i-3)\mathbf{R}\_s(i-3) \\\\ \vdots & \ddots & \vdots\\ 0 & & \sigma\_s^2(i+4)\mathbf{R}\_s(i+4) \end{bmatrix} \tag{9}$$

where **Rs** ( ) *m* and 2 ( ) *s m* ( *mi i i* 4, 3, , 4 ) in the diagonal of the ( )*i* **Rsi** are the coherence coefficient matrix (given by (3)) and the echo power of the pixel pair *m* , respectively. We can notice from (9) that when the SAR images are accurately coregistered, only the elements in **Rs** ( ) *m* are nonzero, while all the other elements of ( )*i* **Rsi** (or ( )*i* **Csi** ) are zero (assuming the complex reflectivity is independent from pixel to pixel and neglecting the noise). In other words, ( )*i* **Rsi** is a block diagonal matrix. However, if the SAR images are not accurately coregistered, the nonzero elements in the submatrices **Rs** ( ) *m* are diffused to other non-diagonal element positions of ( )*i* **Rsi** , as shown in Fig.3(b) and Fig.3(c). The dimensions of the joint correlation function matrix ( )*i* **Rsi** are 18×18. Fig.3(a) is the structure of ( )*i* **Rsi** for accuracy coregistration; Fig.3(b) is the structure of ( )*i* **Rsi** for the coregistration error of 0.5 pixel; Fig.3(c) is the structure of ( )*i* **Rsi** for the coregistration error of 1 pixel. From Fig.3 we can see that when the SAR images are not accurately coregistered, the same imaged ground area can be coregistered to different pixel positions; thus its correlation information appears in other non-diagonal element position of ( )*i* **Rsi** . The gray level denotes the magnitude of each element with the white strongest and black zero.

Robust Interferometric Phase Estimation in InSAR via Joint Subspace Projection 117

**Range(m)**

The joint subspace projection method mentioned above employs the projection of the joint signal subspace onto the corresponding joint noise subspace which is obtained from the eigendecomposition of the joint covariance matrix to estimate the terrain interferometric phase, and takes advantage of the coherence information of the neighboring pixel pairs to auto-coregister the SAR images, where the phase noise is reduced simultaneously. However, the noise subspace dimension of the covariance matrix changes with the coregistration error. For accurate estimating the InSAR interferometric phase, the noise subspace dimension of the covariance matrix must be known, and the performance of the method (i.e. subspace projection method) degrades when the noise subspace dimension is not estimated correctly. In this chapter, an modified joint subspace projection method for InSAR interferometric phase estimation is proposed. In this method, the benefit from the new formulation of joint data vector is that the noise subspace dimension of the covariance matrix is not affected by the coregistration error (i.e., the noise subspace dimension of the

as that of the covariance matrix with accurate coregistration). So the method does not need to calculate the noise subspace dimension before estimating the InSAR interferometric

For avoiding the trouble of calculating the noise subspace dimension before estimating the InSAR interferometric phase, the new formulation of joint data vector is used in the proposed method. An example to construct the new joint data vector **si**( )*i* is shown in Fig.5.

1222212222

*si si si si si si si si si si*

( 3), ( 2), ( 3), ( 6), ( 7), ( 4), ( 4), ( 5), ( 8), ( 9)]*T*

 **si** (11)

1 2 2 2 2 12 2 22

( ) [ ( 1), ( 6), ( 5), ( 2), ( 1), ( ), ( 4), ( 3), ( ), ( 1),

*i si si si si si sisi si sisi*

(a) (b) (c)

**4. Modified interferometric phase estimation via subspace projection**

Fig. 4. The interferograms obtained by the subspace projection method for the accurate coregistration (a), the coregistration error of 0.5 pixels (b) and the coregistration error of one

**Range(m)**

**0 100 200 300 400 500 600 700**

**0 100 200 300 400 500 600 700**

**Azimuth(m)**

 (0 1)

pixel is the same

**Range(m)**

**0 100 200 300 400 500 600 700**

**Azimuth(m)**

corresponding covariance matrix with the coregistration error

The joint data vector **si**( )*i* shown in Fig.5 can be written as

The corresponding covariance matrix ( )*i* **Csi** is given by

**4.1 Data modeling of the modified method**

**Azimuth(m)**

pixel (c).

phase.

Fig. 3. Structures of joint correlation function matrix for different coregistration errors: (a) accurate coregistration; (b) coregistration error of 0.5 pixels; (c) coregistration error of 1 pixel.

Comparing Fig.3 with (4), the conventional estimation approaches based on the single pixel pair will be not feasible if the coregistration error is very large, for example, larger than 0.5 pixel, while our approach can still achieve the optimum estimation due to the use of multiple neighboring pixel pairs.

From the literature[21], we can know that: for the accurate coregistration, the dimensions of the signal subspace and the noise subspace of the joint covariance matrix ( )*i* **Csi** ( 18 18 ) are both 9. For the coregistration error of 1 pixel, the dimensions of the signal subspace and the noise subspace of the joint covariance matrix ( )*i* **Csi** are changed to 12 and 6, respectively. The method can also provide accurate estimation of the terrain interferometric phase even if the coregistration error reaches one pixel only by changing the dimension of the noise subspace from 9 to 6.

As mentioned above, the signal subspace ( ) ( ) *k i rsi* **ai β** ( 1,2, , *k K* ) is orthogonal to the noise subspace (1) (2) ( ) {,, } *M K nsi nsi nsi span* **Nc ββ β** , which is used to estimate the interferometric phase *i* .

The definition of cost function is,

$$J\_2 = \sum\_{k=1}^{K} \sum\_{l=1}^{M-K} \left( \mathbf{a} \mathbf{i}(\varphi\_i) \odot \mathfrak{P}\_{rs}^{(k)} \right)^H \mathfrak{P}\_{ns}^{(l)} \mathfrak{P}\_{ns}^{(l)H} \left( \mathbf{a} \mathbf{i}(\varphi\_i) \odot \mathfrak{P}\_{rs}^{(k)} \right) \tag{10}$$

The minimization of 2*J* can provide the optimum estimate of the interferometric phase *i* .

Fig.s 4 shows the simulation results for various coregistration errors by the subspace projection method. Comparing these figures, we can observe that the large coregistration error has almost no effect on the interferograms obtained by the proposed method. We can see that the subspace projection method in this chapter is robust to large coregistration errors (up to one pixel).

(a) (b) (c)

pixel.

phase *i* .

multiple neighboring pixel pairs.

the noise subspace from 9 to 6.

The definition of cost function is,

errors (up to one pixel).

As mentioned above, the signal subspace ( ) ( ) *k*

noise subspace (1) (2) ( ) {,, } *M K*

2

*J*

1 1

The minimization of 2*J* can provide the optimum estimate of the interferometric phase

Fig.s 4 shows the simulation results for various coregistration errors by the subspace projection method. Comparing these figures, we can observe that the large coregistration error has almost no effect on the interferograms obtained by the proposed method. We can see that the subspace projection method in this chapter is robust to large coregistration

*k l*

*K MK*

Fig. 3. Structures of joint correlation function matrix for different coregistration errors: (a) accurate coregistration; (b) coregistration error of 0.5 pixels; (c) coregistration error of 1

Comparing Fig.3 with (4), the conventional estimation approaches based on the single pixel pair will be not feasible if the coregistration error is very large, for example, larger than 0.5 pixel, while our approach can still achieve the optimum estimation due to the use of

From the literature[21], we can know that: for the accurate coregistration, the dimensions of the signal subspace and the noise subspace of the joint covariance matrix ( )*i* **Csi** ( 18 18 ) are both 9. For the coregistration error of 1 pixel, the dimensions of the signal subspace and the noise subspace of the joint covariance matrix ( )*i* **Csi** are changed to 12 and 6, respectively. The method can also provide accurate estimation of the terrain interferometric phase even if the coregistration error reaches one pixel only by changing the dimension of

*nsi nsi nsi span* **Nc ββ β** , which is used to estimate the interferometric

( ) () () ( )

*k l lH H k i i rsi nsi nsi rsi*

**ai β ββ ai β** (10)

(() ) (() )

*i rsi* **ai β** ( 1,2, , *k K* ) is orthogonal to the

*i* .

Fig. 4. The interferograms obtained by the subspace projection method for the accurate coregistration (a), the coregistration error of 0.5 pixels (b) and the coregistration error of one pixel (c).

### **4. Modified interferometric phase estimation via subspace projection**

The joint subspace projection method mentioned above employs the projection of the joint signal subspace onto the corresponding joint noise subspace which is obtained from the eigendecomposition of the joint covariance matrix to estimate the terrain interferometric phase, and takes advantage of the coherence information of the neighboring pixel pairs to auto-coregister the SAR images, where the phase noise is reduced simultaneously. However, the noise subspace dimension of the covariance matrix changes with the coregistration error. For accurate estimating the InSAR interferometric phase, the noise subspace dimension of the covariance matrix must be known, and the performance of the method (i.e. subspace projection method) degrades when the noise subspace dimension is not estimated correctly. In this chapter, an modified joint subspace projection method for InSAR interferometric phase estimation is proposed. In this method, the benefit from the new formulation of joint data vector is that the noise subspace dimension of the 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). So the method does not need to calculate the noise subspace dimension before estimating the InSAR interferometric phase.

### **4.1 Data modeling of the modified method**

For avoiding the trouble of calculating the noise subspace dimension before estimating the InSAR interferometric phase, the new formulation of joint data vector is used in the proposed method. An example to construct the new joint data vector **si**( )*i* is shown in Fig.5.

The joint data vector **si**( )*i* shown in Fig.5 can be written as

$$\mathbf{s}(i) = \left[ s\_1(i-1), s\_2(i-6), s\_2(i-5), s\_2(i-2), s\_2(i-1), s\_1(i), s\_2(i-4), s\_2(i-3), s\_2(i), s\_2(i+1), \right] \tag{11}$$

$$s\_1(i+3), s\_2(i+2), s\_2(i+3), s\_2(i+6), s\_2(i+7), s\_1(i+4), s\_2(i+4), s\_2(i+5), s\_2(i+8), s\_2(i+9) \right]^T \tag{12}$$

The corresponding covariance matrix ( )*i* **Csi** is given by

Robust Interferometric Phase Estimation in InSAR via Joint Subspace Projection 119

*SAR* 1 *SAR* 2

*i* 1 *i* 6 *i* 5 *i* 2 *i* 1 *i i* 4 *i* 3 *i i* 1 *i* 3 *i*2 *i* 3 *i* 6 *i* 7 *i*4 *i* 4 *i* 5 *i* 8 *i* 9

( 1) ( 1) 0 0 0

**si**'( )*i i* 1 *i* 1 *i i i* 3 *i* 3 *i*4 *i* 4 *i* 6 *i* 5 *i* 2 *i* 4 *i* 3 *i* 1 *i*2 *i* 6 *i* 7 *i* 5 *i* 8 *i* 9

2

*s*

( 3) ( 3) ( 4) ( 4)

**R R**

**s s**

*i i i i*

( 4) ( ) ( 3)

2 2

*i*

*s s*

( 5) ( 2)

*i i*

2

*s*

0 ( 6) 00

2

*i*

*s*

0 0 ( 6) 0

00 0 ( 9)

From (14) we can see that the number of the principal eigenvalues of ( )*i* **Csi'** is 16 Ref.[21]. The noise subspace dimension of the covariance matrix ( )*i* **Csi'** ( 20 20 ) is 4, thus, the noise

2

*i*

*s*

2 2

*s s*

( 2)

*i i*

( 1)

*i*

2

*s*

2

*s*

( 7)

*i*

2 2

*s s*

( 8)

(15)

*i i*

( 5)

*i*

*i*1 *i i*3 *i*4

Fig. 6. Joint data vector for the accurate coregistration.

2 2 2

*s s s*

0 () ()

**s s**

*i i i i*

**R R**

2

*s*

*i*

**si'**

**R**

*i*2 *i* 1 *i i*1 *i* 6 *i*5 *i* 4 *i*3

*i* 2 *i*3 *i* 4 *i* 5 *i* 6 *i* 7 *i* 8 *i* 9

**si**( )*i*

$$\begin{aligned} \mathbf{C}\_{\text{si}}(i) &= E\{\mathbf{si}(i)\mathbf{si}^{H}(i)\} \\ &= \mathbf{a}\mathbf{i}(\varphi\_{i})\mathbf{a}\mathbf{i}^{H}(\varphi\_{i}) \odot \mathbf{R}\_{\text{si}}(i) + \sigma\_{n}^{2}\mathbf{I} \end{aligned} \tag{12}$$

where

( ) 1, , , , ,1, , , , ,1, , , , ,1, , , , *iiii iiii iiii iiii T jjjj jjjj jjjj jjjj i eeee eeee eeee eeee* **ai** and ( )*i* **Rsi** are referred to as the joint generalized steering vector and the joint correlation function matrix of the pixel pair *i*, respectively. The deduction of equation (12) is presented in appendix A.

Fig. 5. Formulation of the joint data vector.

In the following the noise subspace dimension of the joint covariance matrix ( )*i* **Csi** for different coregistration errors are discussed.

### **a. Accurate coregistration.**

For the accurate coregistration, the joint data vector, **si**( )*i* , that is shown in Fig.6, where circles represent SAR image pixels pair whose interferometric phase has to be estimated.

We can rearrange the elements (pixels) of **si**( )*i* to obtain ' **si** ( )*i* as shown in Fig.6, which does not change the eigenvalues of the corresponding covariance matrix[21]. The joint data vector ' **si** ( )*i* , shown in Fig.6, can be written as

$$\mathbf{s}\mathbf{i}'(i) = \left[s\_1(i-1), s\_2(i-1), s\_1(i), s\_2(i), s\_1(i+3), s\_2(i+3), s\_1(i+4), s\_2(i+4), s\_2(i-6), s\_2(i-5)\right] \tag{13}$$

$$s\_2(i-2), s\_2(i-4), s\_2(i-3), s\_2(i+1), s\_2(i+2), s\_2(i+6), s\_2(i+7), s\_2(i+5), s\_2(i+8), s\_2(i+9)\right]^T \tag{13}$$

The corresponding covariance matrix ( )*i* **Csi'** ( 20 20 ) is given by

$$\begin{aligned} \mathbf{C}\_{\text{si}^\circ}(i) &= E\{\mathbf{si}'(i)\mathbf{si}'(i)^H\} \\ &= \mathbf{a}\mathbf{i}(\varphi\_i)\mathbf{a}\mathbf{i}^H(\varphi\_i) \odot \mathbf{R}\_{\text{si}^\circ}(i) + \sigma\_n^2 \mathbf{I} \end{aligned} \tag{14}$$

*H*

*ii n*

*T jjjj jjjj jjjj jjjj*

**si**

( ) { ( ) ( )} ( ) ( ) ()

 **ai** and ( )*i* **Rsi** are referred to as the joint generalized steering vector and the joint correlation function matrix of the pixel pair *i*, respectively. The deduction of equation (12) is presented

*SAR* 1 *SAR* 2

*i* 1 *i* 6 *i* 5 *i* 2 *i* 1 *i i* 4 *i* 3 *i i* 1 *i* 3 *i*2 *i* 3 *i* 6 *i* 7 *i*4 *i* 4 *i* 5 *i* 8 *i* 9

In the following the noise subspace dimension of the joint covariance matrix ( )*i* **Csi** for

For the accurate coregistration, the joint data vector, **si**( )*i* , that is shown in Fig.6, where circles represent SAR image pixels pair whose interferometric phase has to be estimated.

We can rearrange the elements (pixels) of **si**( )*i* to obtain ' **si** ( )*i* as shown in Fig.6, which does not change the eigenvalues of the corresponding covariance matrix[21]. The joint data

2222222222

*si si si si si si si si si si*

( ) ( ) ()

*H*

*H*

*ii n*

**si'**

( 2), ( 4), ( 3), ( 1), ( 2), ( 6), ( 7), ( 5), ( 8), ( 9)]*T*

2

**ai ai R I** (14)

 *i*

 **si** (13)

1 2 121 2 1 2 2 2

( ) [ ( 1), ( 1), ( ), ( ), ( 3), ( 3), ( 4), ( 4), ( 6), ( 5),

*i si si sisisi si si si si si*

( ) { '( ) '( ) }

*iE i i*

The corresponding covariance matrix ( )*i* **Csi'** ( 20 20 ) is given by

**si'**

**C si si**

*iEi i*

*H*

( ) 1, , , , ,1, , , , ,1, , , , ,1, , , , *iiii iiii iiii iiii*

*i eeee eeee eeee eeee*

**si**

where

in appendix A.

*i*1 *i i*3 *i*4

Fig. 5. Formulation of the joint data vector.

different coregistration errors are discussed.

vector ' **si** ( )*i* , shown in Fig.6, can be written as

**a. Accurate coregistration.**

'

**C si si**

2

**ai ai R I** (12)

*i*2 *i* 1 *i i*1 *i* 6 *i*5 *i* 4 *i*3

*i*2 *i*3 *i* 4 *i* 5 *i* 6 *i* 7 *i* 8 *i* 9

 *i*

Fig. 6. Joint data vector for the accurate coregistration.

$$\mathbf{R}\_{\text{si}^{\prime}}(i) = \begin{bmatrix} \sigma\_{s}^{2}(i-1)\mathbf{R}\_{s}(i-1) & 0 & \cdots & 0 & \cdots & 0\\ 0 & \sigma\_{s}^{2}(i)\mathbf{R}\_{s}(i) & & & & \vdots\\ \vdots & \sigma\_{s}^{2}(i+3)\mathbf{R}\_{s}(i+3) & & & & &\\ & & \sigma\_{s}^{2}(i+4)\mathbf{R}\_{s}(i+3) & & &\\ & & & \sigma\_{s}^{2}(i-5) & & \cdots & 0 & \cdots & 0\\ & & & & & \sigma\_{s}^{2}(i-2)\\ \vdots & & & & & \sigma\_{s}^{2}(i-2)\\ \vdots & & & & & \sigma\_{s}^{2}(i-4) & &\\ & & & & & & \sigma\_{s}^{2}(i+1)\\ & & & & & & & \sigma\_{s}^{2}(i+2)\\ & & & & & & & & \sigma\_{s}^{2}(i+2)\\ & & & & & & & & \sigma\_{s}^{2}(i+2)\\ & & & & & & & & \sigma\_{s}^{2}(i+7)\\\\ 0 & \cdots & 0 & & & & & \sigma\_{s}^{2}(i+7)\\ & & & & & & & & \sigma\_{s}^{2}(i+5)\\\\ & & & & & & & & & \sigma\_{s}^{2}(i+5)\\\\ 0 & \cdots & 0 & \cdots & 0 & & & & \sigma\_{s}^{2}(i+9)\\\\ \end{bmatrix} \tag{15}$$

From (14) we can see that the number of the principal eigenvalues of ( )*i* **Csi'** is 16 Ref.[21]. The noise subspace dimension of the covariance matrix ( )*i* **Csi'** ( 20 20 ) is 4, thus, the noise

Robust Interferometric Phase Estimation in InSAR via Joint Subspace Projection 121

( 1) ( 1) 0 0 0

2

*s*

( 3) ( 3) ( 4) ( 4)

**R R**

**s s**

*i i i i*

( 5) ( ) ( 6)

2 2

*i*

*s s*

2

*s*

( 2) ( 1)

*i i*

0 ( 2) 0 0

2

*i*

*s*

0 0 ( 11) 0

0 0 0 ( 13)

 (0 1)

From the above discussion, we know the noise subspace dimension of the covariance matrix

 **pixel.**

pixel in the first satellite image), the joint data vector, **si**( )*i* ,is shown in Fig.8(a).

covariance matrix. The joint data vector **js**( )*i* , shown in Fig.8(b), can be written as

( ) { ( ) ( )}

*iEi i*

**C js js**

**js**

(i.e., the pixel of the image from the second satellite is shifted upwards compared to the

From the literature[21], we can know that the eigenspectrum (i.e., the distribution of eigenvalues) of a covariance matrix is invariant no matter how the elements of the corresponding data vector are permuted. So we can rearrange the elements of **si**( )*i* to obtain **js**( )*i* as shown in Fig.8(b) [21], which does not change the eigenvalues of the corresponding

1 2 121 2 1 2 2 2

*i si si sisisi si si si si si*

( ) [ ( 1), ( 1), ( ), ( ), ( 3), ( 3), ( 4), ( 4), ( 2), ( 2),

 **js** (19)

2 2 2 2 2 22 2 22

( ) ( ) ()

*ii n*

**js**

*H H*

*s i s i s i s i s As Bs i s i sCs D*

( 1), ( 5), ( 6), ( 7), ( ), ( ), ( 8), ( 9), ( ), ( )]*T*

2

**ai ai R I** (20)

 *i*

2

*i*

*s*

2 2

*i*

*s s*

( 7)

*i i*

> 2 2

*s s*

( 8) ( 9)

*i i*

( 10)

> 2 2

*i i*

*s s*

( 12)

(18)

pixel and its direction is upwards

2

*s*

*i*

**si'**

**R**

**c. Coregistration error of**

2

*s*

0 () ()

*i i i i*

**R R**

2 2

**s s**

*s s*

( )*i* **Csi** estimated from the joint data vector **si**( )*i* is 4.

When the azimuth coregistration error is

The corresponding covariance matrix ( )*i* **Cjs** is given by

 (0 1)

subspace dimension of the covariance matrix ( )*i* **Csi** estimated from the joint data vector **si**( )*i* is 4 Ref.[21].

### **b. Coregistration error of one pixel.**

When the azimuth coregistration error is one 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 joint data vector, **si**( )*i* , is shown in Fig.7.

We can rearrange the elements (pixels) of **si**( )*i* to obtain ' **si** ( )*i* as shown in Fig.7, which does not change the eigenvalues of the corresponding covariance matrix[21]. The joint data vector ' **si** ( )*i* , shown in Fig.7, can be written as

$$\mathbf{s}'(i) = \left[s\_1(i-1), s\_2(i-1), s\_1(i), s\_2(i), s\_1(i+3), s\_2(i+3), s\_1(i+4), s\_2(i+4), s\_2(i-2), s\_2(i+2)\right] \tag{16}$$

$$s\_2(i+1), s\_2(i+5), s\_2(i+6), s\_2(i+7), s\_2(i+10), s\_2(i+11), s\_2(i+8), s\_2(i+9), s\_2(i+12), s\_2(i+13)]^T \tag{17}$$

Fig. 7. Joint data vector for the coregistration error of one pixel. The corresponding covariance matrix ( )*i* **Csi'** is given by

$$\begin{aligned} \mathbf{C}\_{\text{si}^i}(i) &= E\{\mathbf{si}'(i)\mathbf{si}'(i)^H\} \\ &= \mathbf{ai}(\rho\_i)\mathbf{ai}^H(\rho\_i) \odot \mathbf{R}\_{\text{si}^i}(i) + \sigma\_n^2 \mathbf{I} \end{aligned} \tag{17}$$

subspace dimension of the covariance matrix ( )*i* **Csi** estimated from the joint data

When the azimuth coregistration error is one 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

We can rearrange the elements (pixels) of **si**( )*i* to obtain ' **si** ( )*i* as shown in Fig.7, which does not change the eigenvalues of the corresponding covariance matrix[21]. The joint data

22222 2 222 2

*si si si si si si si si si si*

( 1), ( 5), ( 6), ( 7), ( 10), ( 11), ( 8), ( 9), ( 12), ( 13)]*T*

 **si** (16)

*SAR* 1 *SAR* 2

*i* 1 *i* 2 *i* 1 *i* 2 *i* 3 *i i i* 1 *i* 4 *i*5 *i* 3 *i*6 *i*7 *i*10 *i*11 *i*4 *i*8 *i*9 *i*12 *i*13

**si**'( )*i i* 1 *i* 1 *i i i* 3 *i* 3 *i*4 *i* 4 *i*2 *i*2 *i* 1 *i*5 *i*6 *i*7 *i* 8 *i* 9

( ) { '( ) '( ) }

*iE i i*

**C si si**

( ) ( ) ()

*H*

*H*

*ii n*

**si'**

Fig. 7. Joint data vector for the coregistration error of one pixel.

**si'**

The corresponding covariance matrix ( )*i* **Csi'** is given by

*i*2 *i* 1 *i i*1

*i* 2 *i*3 *i* 4 *i* 5 *i* 6 *i* 7 *i* 8 *i* 9

*i*11 *i*12 *i*13

*i*10 *i*11 *i*12 *i*13

*i*10

2

**ai ai R I** (17)

 *i*

1 2 121 2 1 2 2 2

*i si si sisisi si si si si si*

( ) [ ( 1), ( 1), ( ), ( ), ( 3), ( 3), ( 4), ( 4), ( 2), ( 2),

first satellite image), the joint data vector, **si**( )*i* , is shown in Fig.7.

vector **si**( )*i* is 4 Ref.[21].

'

**si**( )*i*

**b. Coregistration error of one pixel.**

vector ' **si** ( )*i* , shown in Fig.7, can be written as

*i*1 *i i*3 *i*4

From the above discussion, we know the noise subspace dimension of the covariance matrix ( )*i* **Csi** estimated from the joint data vector **si**( )*i* is 4.

### **c. Coregistration error of**  (0 1)  **pixel.**

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 joint data vector, **si**( )*i* ,is shown in Fig.8(a).

From the literature[21], we can know that the eigenspectrum (i.e., the distribution of eigenvalues) of a covariance matrix is invariant no matter how the elements of the corresponding data vector are permuted. So we can rearrange the elements of **si**( )*i* to obtain **js**( )*i* as shown in Fig.8(b) [21], which does not change the eigenvalues of the corresponding covariance matrix. The joint data vector **js**( )*i* , shown in Fig.8(b), can be written as

$$\mathbf{j}(i) = \left[ s\_1(i-1), s\_2(i-1), s\_1(i), s\_2(i), s\_1(i+3), s\_2(i+3), s\_1(i+4), s\_2(i+4), s\_2(i-2), s\_2(i+2) \right] \tag{19}$$

$$s\_2(i+1), s\_2(i+5), s\_2(i+6), s\_2(i+7), s\_2(A), s\_2(B), s\_2(i+8), s\_2(i+9), s\_2(C), s\_2(D) \right]^T \tag{10}$$

The corresponding covariance matrix ( )*i* **Cjs** is given by

$$\begin{aligned} \mathbf{C}\_{\mathbf{js}}(i) &= E\{\mathbf{js}(i)\} \mathbf{s}^H(i) \| \\ &= \mathbf{a} \mathbf{i}(\varphi\_i) \mathbf{a} \mathbf{i}^H(\varphi\_i) \odot \mathbf{R}\_{\mathbf{js}}(i) + \sigma\_n^2 \mathbf{I} \end{aligned} \tag{20}$$

Robust Interferometric Phase Estimation in InSAR via Joint Subspace Projection 123

( 1) ( 1) 0 0 0

0 ( 2) 0 0

2

*s*

*i*

**js**

**R**

vector **si**( )*i* is also 4[21].

coregistration error

interferometric phase.

**4.2 Summary of the modified method**

method based on subspace projection.

 (0 1)

2

*s s s*

0 () ()

*i i i i*

**R R**

2 2

**s s**

2

*s*

( 3) ( 3) ( 4) ( 4)

**R R**

**s s**

*i i i i*

( 5) ( ) ( 6)

2 2

*i*

*s s*

2

*s*

( 2) ( 1)

*i i*

2

*i*

*s*

0 0 ( ) 0

0 0 0 ( )

We know the noise subspace dimension of the covariance matrix ( )*i* **Cjs** ( 20 20 ) is 4, thus, the noise subspace dimension of the covariance matrix ( )*i* **Csi** estimated from the joint data

From the results derived above, we can see that: the new formulation of joint data vector proposed in this chapter has the advantage that the noise subspace dimension of the corresponding covariance matrix is independent of the coregistration error. That is to say, the noise subspace dimension of the corresponding covariance matrix with the

Therefore, it is not required to calculate the noise subspace dimension, thus avoiding the trouble of calculating the noise subspace dimension before estimating the InSAR

In this section, we give the detailed steps for the modified interferometric phase estimation

**Step 1.** Coregister SAR images. The SAR images are coarsely coregistered using the

mitigate the complexity in image coregistration processing.

crosscorrelation information of the SAR image intensity or other strategies[1][2] after SAR imaging of the echoes acquired by each satellite. However, the allowable coregistration error of the proposed method can be very large (such as one pixel), which is useful in practice. The low coregistration accuracy requirement can greatly

2

*i*

*s*

2 2

*A B*

*s s*

( 7) ( )

*i*

2

*s*

pixel is the same as that of the accurate covariance matrix.

2

*s*

( 8)

*i*

2 2

*s s*

( )

*D*

(21)

*C*

( 9)

*i*

( ) *a*

Fig. 8. Joint data vector for the coregistration error of pixels.

*SAR* 1 *SAR* 2

'**si** ( )*i*

*i*1 *i i*3 *i*4

*SAR* 1 *SAR* 2

( 6)*b i* ( 1)*a i* ( 5)*b i*

( 10)*a i*

Fig. 8. Joint data vector for the coregistration error of

( 2)*b i*

( 2)*a i* ( 3)*a i* ( 1)*b* ( 2) *i b i*

*i*1 *i i*3 *i*4

*ai* ( 1)*a* ( 2) *i a i* ( 1)*a i*

( 3)*b* ( 6) *i b* ( 5) *i b* ( 4) *i b i*

( 2)*a i* ( 3)*a i* ( 4)*a i* ( 5)*a i* ( 6)*a i* ( 7)*a i* ( 8)*a i* ( 9)*a i*

( 6)*b i* ( 7)*b i* ( 8)*b i* ( 9) *b i*

( 1)*b i bi* ( 4)*b i* ( 5)*b* ( 2) *i b* ( 3) *i b i*

( 11)*a i* ( 12)*a i* ( 13)*a i*

( 1)*b i*

**si**( )*i i* 1 ( 2) *i i* 3 *i*4 *a i*

*ai* ( 1)*a i* ( 4)*a i* ( 5)*a i bi* ( 1)*b* ( 4) *i b* ( 3) *i b i*

*i* 1 *i* 2 *i* 1 *i* 2 *i* 3 *i i i* 1 *i* 4 *i*5 *i* 3 *i*6 *i*7 *A B i*4 *i*8 *i*9 *C D*

( ) *b*

pixels.

( 10)*a i*

( 2)*b i*

( ) *a*

*SAR* 2

( 10)*a i*

( 2)*b i*

*ai* ( 1)*a* ( 2) *i a i* ( 1)*a i*

( 1)*b i bi* ( 4)*b i* ( 5)*b* ( 2) *i b* ( 3) *i b i*

( 1)*b i*

( 10)*a* ( 6) *i a* ( 7) *i a i* ( 11)*a i* ( 3)*b* ( 2) *i b i* ( 6)*b i* ( 7)*b i*

*ai* ( 1)*a* ( 2) *i a i* ( 1)*a i*

( 3)*b* ( 6) *i b* ( 5) *i b* ( 4) *i b i*

( 2)*a i* ( 3)*a i* ( 4)*a i* ( 5)*a i* ( 6)*a i* ( 7)*a i* ( 8)*a i* ( 9)*a i*

( 6)*b i* ( 7)*b i* ( 8)*b i* ( 9) *b i*

( 1)*b i bi* ( 4)*b i* ( 5)*b* ( 2) *i b* ( 3) *i b i*

( 11)*a i* ( 12)*a i* ( 13)*a i*

( 1)*b i*

( 2)*a i* ( 3)*a i* ( 4)*a i* ( 5)*a i* ( 6)*a i* ( 7)*a i* ( 8)*a i* ( 9)*a i*

( 6)*b i* ( 7)*b i* ( 8)*b i* ( 9) *b i*

( 3)*b* ( 5) ( 4) *i b i b* ( 6) *i b i*

( 11)*a i* ( 12)*a i* ( 13)*a i*

*A B C D*

*A*

*SAR* 2

( 8)*a i* ( 9)*a i* ( 12)*a i* ( 13)*a i* ( 4)*b i* ( 5)*b i* ( 8)*b i* ( 9)*b i*

*i*2 *i* 1 *i i*1

*i*2 *i*3 *i*4 *i* 5 *i* 6 *i* 7 *i* 8 *i* 9

*B C D*

**js**( )*i i* 1 *i* 1 *i i i*3 *i* 3 *i*4 *i* 4 *i*2 *i*2 *i* 1 *i*5 *i*6 *i*7 *i* 8 *i* 9

We know the noise subspace dimension of the covariance matrix ( )*i* **Cjs** ( 20 20 ) is 4, thus, the noise subspace dimension of the covariance matrix ( )*i* **Csi** estimated from the joint data vector **si**( )*i* is also 4[21].

From the results derived above, we can see that: the new formulation of joint data vector proposed in this chapter has the advantage that the noise subspace dimension of the corresponding 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. Therefore, it is not required to calculate the noise subspace dimension, thus avoiding the trouble of calculating the noise subspace dimension before estimating the InSAR interferometric phase.

### **4.2 Summary of the modified method**

In this section, we give the detailed steps for the modified interferometric phase estimation method based on subspace projection.

**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. However, the allowable coregistration error of the proposed method can be very large (such as one pixel), which is useful in practice. The low coregistration accuracy requirement can greatly mitigate the complexity in image coregistration processing.

Robust Interferometric Phase Estimation in InSAR via Joint Subspace Projection 125

As shown by (8), the same signal subspace spanned by the principal eigenvectors

*csi* **β** ( *m K* 1,2, , )of ˆ ( )*i* **Csi** can be spanned by the Hadamard product vectors ˆ( ) ( ) *m*

*ii i rsi rsi rsi* **S ai c** *span*

**Step 4.** Projection of signal subspace onto noise subspace. The projection of the signal subspace onto the corresponding noise subspace is performed as follows:

( ) 1, , , , ,1, , , , ,1, , , , ,1, , , , *iiii iiii iiii iiii*

The cost function given by (29) is used to estimate the terrain interferometric phase

projection to coregistration errors by using simulated data and real data.

the minimization of 3*J* can provide the optimum estimate of the interferometric phase

Remark 2: The computational burden will be high if the minimization of 3*J* is obtained via

fast algorithm to compute the minimization of 3*J* is developed in Appendix B, where the

Using the above four steps, the terrain interferogram can be recovered after the pixel pairs of

In this section we demonstrate the robustness of the modified method via subspace

The simulated data are described as follows. We assume there are two formation-flying satellites in the cartwheel formation, and we select one orbit position for simulation, with an effective cross-track baseline of 281.46 m, an orbit height of 750 kilometers and an incidence angle of 45. We use a two-dimensional window to simulate the terrain and use the statistical model to generate the complex SAR image pairs[23]. The signal-to-noise ratio

Here, the number of the samples to estimate the covariance matrix is 7 (in range) 7 (in

Fig.s 9-12 compare the simulation results for various techniques and coregistration errors. Comparing Fig. 9,10 ,11and 12, we can observe that the large coregistration error heavily affects the interferograms obtained by pivoting median filtering, pivoting mean filtering and adaptive contoured window filtering. On the contrary, the large coregistration error has

*i eeee eeee eeee eeee*

ˆˆ ˆ (1) (2) ( ) () ,() , ,() *K*

 **β ai β ai β** (28)

>

. To reduce the computational burden, a

*i* is directly obtained by using the fast algorithm.

*i* . And

> *i* ,

**ai β ββ ai β** (29)

( ) () () ( )

*m l lH H m i i rsi nsi nsi rsi*

ˆ ˆˆ ˆ (() ) (() )

*T jjjj jjjj jjjj jjjj*

**ai** (30)

*rsi* **β** ( *m K* 1,2, , ).

By eigen-decomposing ˆ ( )*i* **Rsi** , we obtain *K* principal eigenvectors ˆ( ) *m*

20

*K K*

1 1

*i* in the principal phase interval [ , ]

*m l*

3

closed-form solution to the estimate of

(SNR) of the SAR images is 18 dB.

azimuth)=49.

the SAR images are processed separately.

**5. Numerical and experemental results**

*J*

ˆ( ) *m*

where

i.e., ˆ *i i* .

search of

*i rsi* **ai β** ( *m K* 1,2, , ), i.e.,

**Step 2.** Estimate the covariance matrix. The covariance matrix ( )*i* **Csi** can be estimated by using joint data vector **si**( )*i* shown in Fig.5. Under the assumption that the neighboring pixels have the identical terrain height and the complex reflectivity is independent from pixel to pixel[20][21], the covariance matrix ( )*i* **Csi** can be estimated by its sample covariance matrix ˆ ( )*i* **Csi** , i.e.,

$$\hat{\mathbf{C}}\_{\text{si}}(i) = \frac{1}{2K + 1} \sum\_{k=-L}^{L} \mathbf{s} \mathbf{i}(i+k) \mathbf{s} \mathbf{i}^H(i+k) \tag{22}$$

where 2 1 *L* is the number of i.i.d. samples from the neighboring pixel pairs.

Remark 1: It is easy to obtain enough i.i.d. samples for locally flat terrains. However, an imaging terrain in practice can not 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[20], the steering vector (i.e., the interferometric phase) variation due to the different terrain height from pixel to pixel can be compensated, which greatly enlarges the size of the sample window.

**Step 3.** Subspace estimation by Eigendecomposing. The estimated covariance matrix ˆ ( )*i* **Csi** of the dimensions 20 20 can be eigendecomposed into

$$
\hat{\mathbf{C}}\_{\text{si}}(i) = \sum\_{m=1}^{K} \hat{\lambda}\_{\text{csi}}^{(m)} \hat{\mathfrak{P}}\_{\text{csi}}^{(m)} \hat{\mathfrak{P}}\_{\text{csi}}^{(m)H} + \sum\_{l=1}^{20-K} \hat{\lambda}\_{\text{csi}}^{(l+K)} \hat{\mathfrak{P}}\_{\text{nsi}}^{(l)} \hat{\mathfrak{P}}\_{\text{nsi}}^{(l)H} \tag{23}
$$

where *K* is the number of the principal eigenvalues of ˆ ( )*i* **Csi** , ˆˆ ˆ ˆ ˆ (1) (2) ( ) ( 1) *K K* (20) *csi csi csi csi csi* , eigenvectors ˆ( )*l nsi* **β** ( 1,2, ,20 *l K* ) corresponding to the smaller eigenvalues ˆ( ) *l K csi* ( 1,2, ,20 *l K* ) span the noise subspace, i.e,

$$\mathbf{N}\_{\mathbf{c}} = \text{span}\left(\hat{\mathsf{B}}\_{n\bar{s}}^{(1)}, \hat{\mathsf{B}}\_{n\bar{s}}^{(2)}, \dots, \hat{\mathsf{B}}\_{n\bar{s}}^{(20-K)}\right) \tag{24}$$

whereas the larger eigenvectors ˆ( ) *m csi* **β** ( 1,2, , *m K* ) corresponding to the principal eigenvalues ˆ( ) *m csi* ( *m K* 1,2, , ) span the signal subspace, i.e.,

$$\mathbf{S}\_{\mathbf{c}} = \text{span}\left(\hat{\mathsf{B}}\_{\text{cs}}^{(1)}, \hat{\mathsf{B}}\_{\text{cs}}^{(2)}, \dots, \hat{\mathsf{B}}\_{\text{cs}}^{(K)}\right) \tag{25}$$

The noise power is often estimated by

$$
\hat{\sigma}\_n^2 = \frac{1}{20 - K} \sum\_{l=1}^{20 - K} \hat{\lambda}\_{\text{csi}}^{(l+K)} \tag{26}
$$

The joint correlation function matrix ˆ ( )*i* **Rsi** can be approximated as the amplitude (i.e., the absolute value) of the estimated covariance matrix ˆ ( )*i* **Csi** [20], i.e.,

$$
\hat{\mathbf{R}}\_{\rm si}(i) = \left| \hat{\mathbf{C}}\_{\rm si}(i) - \hat{\sigma}\_n^2 \mathbf{I} \right| \tag{27}
$$

By eigen-decomposing ˆ ( )*i* **Rsi** , we obtain *K* principal eigenvectors ˆ( ) *m rsi* **β** ( *m K* 1,2, , ). As shown by (8), the same signal subspace spanned by the principal eigenvectors ˆ( ) *m csi* **β** ( *m K* 1,2, , )of ˆ ( )*i* **Csi** can be spanned by the Hadamard product vectors ˆ( ) ( ) *m i rsi* **ai β** ( *m K* 1,2, , ), i.e.,

$$\mathbf{S}\_{\mathbf{c}} = \text{span}\left(\mathbf{a}\mathbf{i}(\rho\_i) \odot \hat{\mathfrak{P}}\_{\text{rs}}^{(1)}, \mathbf{a}\mathbf{i}(\rho\_i) \odot \hat{\mathfrak{P}}\_{\text{rs}}^{(2)}, \dots, \mathbf{a}\mathbf{i}(\rho\_i) \odot \hat{\mathfrak{P}}\_{\text{rs}}^{(K)}\right) \tag{28}$$

**Step 4.** Projection of signal subspace onto noise subspace. The projection of the signal subspace onto the corresponding noise subspace is performed as follows:

$$J\_3 = \sum\_{m=1}^{K} \sum\_{l=1}^{20-K} (\mathbf{a} \mathbf{i}(\phi\_i) \odot \hat{\mathfrak{P}}\_{rsi}^{(m)})^H \hat{\mathfrak{P}}\_{nsi}^{(l)} \hat{\mathfrak{P}}\_{nsi}^{(l)H} (\mathbf{a} \mathbf{i}(\phi\_i) \odot \hat{\mathfrak{P}}\_{rsi}^{(m)}) \tag{29}$$

where

124 Recent Interferometry Applications in Topography and Astronomy

**Step 2.** Estimate the covariance matrix. The covariance matrix ( )*i* **Csi** can be estimated by

1 ˆ ( ) ( )( ) 2 1 *L*

Remark 1: It is easy to obtain enough i.i.d. samples for locally flat terrains. However, an imaging terrain in practice can not 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[20], the steering vector (i.e., the interferometric phase) variation due to the different terrain height from pixel to pixel can be compensated, which greatly enlarges the

**Step 3.** Subspace estimation by Eigendecomposing. The estimated covariance matrix ˆ ( )*i* **Csi**

1 1 ˆ ˆ ˆ ˆ ˆ ˆ ˆ ( ) *K K*

where *K* is the number of the principal eigenvalues of

 , eigenvectors ˆ( )*l*

> *csi*

*m l*

20 ()()() ( ) () ()

ˆˆ ˆ (1) (2) (20 ) , ,, *K*

ˆˆ ˆ (1) (2) ( ) , ,, *K*

20 2 ( ) 1

*K l K*

1 ˆ ˆ 20

*n csi K l*

The joint correlation function matrix ˆ ( )*i* **Rsi** can be approximated as the amplitude (i.e., the

ˆ ˆ 2 () () ˆ*n* **RC I si si** *i i*

*m m mH l K l lH csi csi csi csi nsi nsi*

**Csi β β β β** (23)

( 1,2, ,20 *l K* ) span the noise subspace,

*nsi nsi nsi span* **Nc ββ β** (24)

*csi* **β** ( 1,2, , *m K* ) corresponding to the principal

*csi csi csi* **Sc** *span* **ββ β** (25)

(26)

(27)

*nsi* **β** ( 1,2, ,20 *l K* )

*k L i ik ik*

*K*

where 2 1 *L* is the number of i.i.d. samples from the neighboring pixel pairs.

of the dimensions 20 20 can be eigendecomposed into

( *m K* 1,2, , ) span the signal subspace, i.e.,

absolute value) of the estimated covariance matrix ˆ ( )*i* **Csi** [20], i.e.,

*i*

ˆ ( )*i* **Csi** , ˆˆ ˆ ˆ ˆ (1) (2) ( ) ( 1) *K K* (20) *csi csi csi csi csi*

corresponding to the smaller eigenvalues ˆ( ) *l K*

whereas the larger eigenvectors ˆ( ) *m*

The noise power is often estimated by

*csi*

estimated by its sample covariance matrix ˆ ( )*i* **Csi** , i.e.,

size of the sample window.

eigenvalues ˆ( ) *m*

i.e,

using joint data vector **si**( )*i* shown in Fig.5. Under the assumption that the neighboring pixels have the identical terrain height and the complex reflectivity is independent from pixel to pixel[20][21], the covariance matrix ( )*i* **Csi** can be

*H*

**Csi si si** (22)

$$\mathbf{a}\mathbf{i}(\boldsymbol{\phi}\_{i}) = \left[\mathbf{1}, \mathbf{e}^{j\phi\_{i}}, \mathbf{e}^{j\phi\_{i}}, \mathbf{e}^{j\phi\_{i}}, \mathbf{e}^{j\phi\_{i}}, \mathbf{1}, \mathbf{e}^{j\phi\_{i}}, \mathbf{e}^{j\phi\_{i}}, \mathbf{e}^{j\phi\_{i}}, \mathbf{e}^{j\phi\_{i}}, \mathbf{1}, \mathbf{e}^{j\phi\_{i}}, \mathbf{e}^{j\phi\_{i}}, \mathbf{e}^{j\phi\_{i}}, \mathbf{1}, \mathbf{e}^{j\phi\_{i}}, \mathbf{e}^{j\phi\_{i}}, \mathbf{e}^{j\phi\_{i}}, \mathbf{e}^{j\phi\_{i}}\right]^{T} \tag{30}$$

The cost function given by (29) is used to estimate the terrain interferometric phase *i* . And the minimization of 3*J* can provide the optimum estimate of the interferometric phase *i* , i.e., ˆ *i i* .

Remark 2: The computational burden will be high if the minimization of 3*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 3*J* is developed in Appendix B, where the closed-form solution to the estimate of *i* is directly obtained by using the fast algorithm.

Using the above four steps, the terrain interferogram can be recovered after the pixel pairs of the SAR images are processed separately.

### **5. Numerical and experemental results**

In this section we demonstrate the robustness of the modified method via subspace projection to coregistration errors by using simulated data and real data.

The simulated data are described as follows. We assume there are two formation-flying satellites in the cartwheel formation, and we select one orbit position for simulation, with an effective cross-track baseline of 281.46 m, an orbit height of 750 kilometers and an incidence angle of 45. We use a two-dimensional window to simulate the terrain and use the statistical model to generate the complex SAR image pairs[23]. The signal-to-noise ratio (SNR) of the SAR images is 18 dB.

Here, the number of the samples to estimate the covariance matrix is 7 (in range) 7 (in azimuth)=49.

Fig.s 9-12 compare the simulation results for various techniques and coregistration errors. Comparing Fig. 9,10 ,11and 12, we can observe that the large coregistration error heavily affects the interferograms obtained by pivoting median filtering, pivoting mean filtering and adaptive contoured window filtering. On the contrary, the large coregistration error has

Robust Interferometric Phase Estimation in InSAR via Joint Subspace Projection 127

**Range(m)**

Fig. 12. The interferograms obtained by the modified method for the accurate coregistration (a), the coregistration error of 0.5 pixels (b) and the coregistration error of one pixel (c).

In the following, we will verify the validity of the modified method with the ERS1/ERS2

Fig. 13 shows the interferograms generated from the ERS1/ERS2 real data. Fig.13(a) is the interferogram obtained by the conventional processing, and Fig.13(b) is that obtained by the

(a) (b) (c)

**Range(m)**

**0 100 200 300 400 500 600 700**

**0 100 200 300 400 500 600 700**

**Azimuth(m)**

**Range(m)**

modified method proposed in this chapter.

modified method for the ERS1/ ERS2 real data.

**6. Conclusions**

of the proposed new algorithms.

**0 100 200 300 400 500 600 700**

**Azimuth(m)**

(European Remote Sensing 1 and 2 tandem satellites) real data.

(a) (b)

Fig. 13. The interferograms obtained by (a) the conventional processing, and (b) the

In this chapter, the interferometric phase estimation method based on subspace projection and its modified method were presented. The interferometric phase estimation methods based on joint subspace projection can provide accurate estimate of the terrain interferometric phase (interferogram) even if the coregistration error reaches one pixel. Benefiting from the new formulation of joint data vector, the modified method does not need to calculate the noise subspace dimension, thus avoiding the trouble of calculating the noise subspace dimension before estimating the InSAR interferometric phase. A fast algorithm is developed to implement the modified method, which can significantly reduce the computational burden. Theoretical analysis and simulations demonstrate the efficiency

**Azimuth(m)**

almost no effect on the interferograms obtained by the modified method. We can see that the modified method in this chapter is robust to large coregistration errors (up to one pixel).

Fig. 9. The interferograms obtained by the pivoting median filtering for the accurate coregistration (a), the coregistration error of 0.5 pixels (b) and the coregistration error of one pixel (c).

Fig. 10. The interferograms obtained by the pivoting mean filtering for the accurate coregistration (a), the coregistration error of 0.5 pixels (b) and the coregistration error of one pixel (c).

Fig. 11. The interferograms obtained by the adaptive contoured window filtering for the accurate coregistration (a), the coregistration error of 0.5 pixels (b) and the coregistration error of one pixel (c).

Fig. 12. The interferograms obtained by the modified method for the accurate coregistration (a), the coregistration error of 0.5 pixels (b) and the coregistration error of one pixel (c).

In the following, we will verify the validity of the modified method with the ERS1/ERS2 (European Remote Sensing 1 and 2 tandem satellites) real data.

Fig. 13 shows the interferograms generated from the ERS1/ERS2 real data. Fig.13(a) is the interferogram obtained by the conventional processing, and Fig.13(b) is that obtained by the modified method proposed in this chapter.

Fig. 13. The interferograms obtained by (a) the conventional processing, and (b) the modified method for the ERS1/ ERS2 real data.

### **6. Conclusions**

126 Recent Interferometry Applications in Topography and Astronomy

almost no effect on the interferograms obtained by the modified method. We can see that the modified method in this chapter is robust to large coregistration errors (up to one pixel).

 (a) (b) (c) Fig. 9. The interferograms obtained by the pivoting median filtering for the accurate

 (a) (b) (c) Fig. 10. The interferograms obtained by the pivoting mean filtering for the accurate

(a) (b) (c)

Fig. 11. The interferograms obtained by the adaptive contoured window filtering for the accurate coregistration (a), the coregistration error of 0.5 pixels (b) and the coregistration

**Azimuth(m)**

coregistration (a), the coregistration error of 0.5 pixels (b) and the coregistration error of one

**Azimuth(m)**

pixel (c).

pixel (c).

**Azimuth(m)**

error of one pixel (c).

coregistration (a), the coregistration error of 0.5 pixels (b) and the coregistration error of one

In this chapter, the interferometric phase estimation method based on subspace projection and its modified method were presented. The interferometric phase estimation methods based on joint subspace projection can provide accurate estimate of the terrain interferometric phase (interferogram) even if the coregistration error reaches one pixel. Benefiting from the new formulation of joint data vector, the modified method does not need to calculate the noise subspace dimension, thus avoiding the trouble of calculating the noise subspace dimension before estimating the InSAR interferometric phase. A fast algorithm is developed to implement the modified method, which can significantly reduce the computational burden. Theoretical analysis and simulations demonstrate the efficiency of the proposed new algorithms.

Robust Interferometric Phase Estimation in InSAR via Joint Subspace Projection 129

2

**I**

*n*

1)

*n*

**I**

(A.4)

(A.2)

( ) ( ) , ( ) ( 4) , ( ) ( 3) , ( ) ( ) , ( ) ( 1)

( 3) ( ) , ( 3) ( 4) , ( 3) ( 3) , ( 3) ( ) , ( 3) ( 1)

21 22 22 22 22

( 1) ( ) , ( 1) ( 4) , ( 1) ( 3) , ( 1) ( ) , ( 1) ( 1)

*i ii i*

*j jj j*

21 22 22 22 22 21 22 2 2 22 22

*Es i s i Es i s i Es i s i Es i s i Es i s i r ii e r ii e r iie r ii e*

( ) ( ) , ( ) ( 4) , ( ) ( 3) , ( ) ( ) , ( ) (

12 12 12 12

1, ( , 4) , ( , 3) , ( , ) , ( , 1)

( 3, ) , ( 3, 4), 1, ( 3, ), ( 3, 1) ( , ) , ( , 4), ( , 3), 1, ( , 1) ( 1, ) , ( 1, 4), ( 1, 3), ( 1, ), 1

1, , , , 1, ( , 4), ( , 3), ( , ), ( , 1) , 1, 1, 1, 1 ( 4, ), 1, ( 4, 3), ( 4, ), ( 4, 1)

( ) 1, , , , *iiii*

*i eeee*

21 22 22 22

21 22 22 22 21 22 22 22

*r ii r ii r ii r ii*

( ) ( ) ( 3, ), ( 3, 4),1, ( 3, ), ( 3, 1)

So the covariance matrix ( )*i* **Csi** of the joint data vector **si**( )*i* , shown in Fig.5, can be given by

*H*

*ii n*

**si**

( ) { ( ) ( )} ( ) ( ) ()

*iEi i*

*H*

**si**

**9.1 Fast algorithm for optimal interferometric phase estimation**

If *U* , *V* and *W* are arbitrary complex column vectors, then[20]

**C si si**

*i i r i i r i i r i ir i i*

*T jjjj*

12 12 12 12

( , ), ( , 4), ( , 3), 1, ( , 1)

22 22 22 ( 1, ), ( 1, 4), ( 1, 3), ( 1, ),1 *i ir i i r i i r i i*

1, ( , 4), ( , 3), ( , ), ( , 1) ( 4, ), 1, ( 4, 3), ( 4, ), ( 4, 1)

*r ii r ii r ii r ii r i i r i i r i ir i i*

2

**ai ai R I** (A.5)

 *i*

**ρ** (A.3)

*eeee r ii r ii r ii r ii e r i i r i i r i ir i i*

*Es i s i Es i s i Es i s i Es i s i Es i s i Es is i Es is i Es is i Es is i Es is i*

22 22

3), ( 4, ), ( 4, 1)

*ri i ri i*

12 12 12 12 21 22 22 22

 2 1 22 22 22 21 22 22 22 21 22 22 22

*i ir i i r i ir i i r ii r ii r ii r ii r i ir i i r i i r i i*

( 3, ), ( 3, 4),1, ( 3, ), ( 3, 1) ( , ), ( , 4), ( , 3), 1, ( , 1) ( 1, ), ( 1, 4), ( 1, 3), ( 1, ),1

22 22 22

4) , ( 4) ( 3) , ( 4) ( ) , ( 4) ( 1)

*Es i s i Es i s i Es i s i*

1 1 1 2 1 2 12 12

*Es is i Es is i Es is i Es is i Es is i*

12 2 22 12 2 22

*Esis i s i s is i sis i s i s is i*

21 22 22 22 21 22 22 22

*r i ie r i i ri i ri i r iie r ii r ii r ii*

{[ ( ), ( 4), ( 3), ( ), ( 1)] [ ( ), ( 4), ( 3), ( ), ( 1)] }

*T*

( ) { ( ) ( )}

*iE i i*

**C ss ss ss**

2

*s*

( )

*i*

*i i i i*

*j j*

*j j*

*e e*

where

21 22

*Es i s i Es i s i*

*H*

( 4) ( ) , ( 4) (

21 22

*i*

*i i i*

*iiii*

*jjjj*

, 1, 1, 1, 1 ( )

*e i r*

*i*

2

*s*

21

*r*

  **I**

*j*

*j j j*

, 1, 1, 1, 1 , 1, 1, 1, 1

*ii n*

( ) ( ) ()

**ρ ρ R** 2

**ss**

**Rss**

**9. Appendix B**

*H*

( 4, ) , 1, ( 4,

21 22 22 22

*r i ie r i i r i i r i i*

2 2

*s*

*r i ie r i i*

### **7. Acknowledgement**

This work is supported in part by the National Natural Science Foundations of China under grant 60736009, 61071194 and 60979002, by the Fund of Civil Aviation University of China under grant ZXH2009D018 and 2011kyE06.

### **8. Appendix A**

### **8.1 Proof of equation (12)**

For easy discussion, we assume that the structure of joint data vector **ss**( )*i* is shown in Fig.A.1, where circles represent SAR image pixels and *i* denotes the desired pixel pair whose interferometric phase is to be estimated.

Fig. A.1. Formulation of the joint data vector **ss**( )*i* .

The joint data vector **ss**( )*i* , shown in Fig.A.1, can be written as

$$\mathbf{ss}(i) = \begin{bmatrix} s\_1(i), s\_2(i-4), s\_2(i-3), s\_2(i), s\_2(i+1) \end{bmatrix}^T \tag{A.1}$$

The corresponding covariance matrix ( )*i* **Css** is given by

 12 2 22 12 2 22 1 1 1 2 1 2 12 12 21 22 ( ) { ( ) ( )} {[ ( ), ( 4), ( 3), ( ), ( 1)] [ ( ), ( 4), ( 3), ( ), ( 1)] } ( ) ( ) , ( ) ( 4) , ( ) ( 3) , ( ) ( ) , ( ) ( 1) ( 4) ( ) , ( 4) ( *H T iE i i Esis i s i s is i sis i s i s is i Es is i Es is i Es is i Es is i Es is i Es i s i Es i s i* **C ss ss ss** 22 22 22 21 22 22 22 22 21 22 2 2 22 22 4) , ( 4) ( 3) , ( 4) ( ) , ( 4) ( 1) ( 3) ( ) , ( 3) ( 4) , ( 3) ( 3) , ( 3) ( ) , ( 3) ( 1) ( ) ( ) , ( ) ( 4) , ( ) ( 3) , ( ) ( ) , ( ) ( *Es i s i Es i s i Es i s i Es i s i Es i s i Es i s i Es i s i Es i s i Es is i Es is i Es is i Es is i Es is i* 21 22 22 22 22 12 12 12 12 21 22 2 1) ( 1) ( ) , ( 1) ( 4) , ( 1) ( 3) , ( 1) ( ) , ( 1) ( 1) 1, ( , 4) , ( , 3) , ( , ) , ( , 1) ( 4, ) , 1, ( 4, ( ) *i ii i i j jj j j s Es i s i Es i s i Es i s i Es i s i Es i s i r ii e r ii e r iie r ii e r i ie r i i i* 22 22 21 22 22 22 21 22 22 22 21 22 22 22 3), ( 4, ), ( 4, 1) ( 3, ) , ( 3, 4), 1, ( 3, ), ( 3, 1) ( , ) , ( , 4), ( , 3), 1, ( , 1) ( 1, ) , ( 1, 4), ( 1, 3), ( 1, ), 1 *i i i j j j ri i ri i r i ie r i i ri i ri i r iie r ii r ii r ii r i ie r i i r i i r i i* 2 12 12 12 12 21 22 22 22 2 2 1, , , , 1, ( , 4), ( , 3), ( , ), ( , 1) , 1, 1, 1, 1 ( 4, ), 1, ( 4, 3), ( 4, ), ( 4, 1) , 1, 1, 1, 1 ( ) , 1, 1, 1, 1 , 1, 1, 1, 1 *iiii i i i i n jjjj j j s j j eeee r ii r ii r ii r ii e r i i r i i r i ir i i e i r e e* **I** 2 1 22 22 22 21 22 22 22 21 22 22 22 ( 3, ), ( 3, 4),1, ( 3, ), ( 3, 1) ( , ), ( , 4), ( , 3), 1, ( , 1) ( 1, ), ( 1, 4), ( 1, 3), ( 1, ),1 ( ) ( ) () *n H ii n i ir i i r i ir i i r ii r ii r ii r ii r i ir i i r i i r i i i* **ss I ρ ρ R** 2 **I** (A.2)

where

128 Recent Interferometry Applications in Topography and Astronomy

This work is supported in part by the National Natural Science Foundations of China under grant 60736009, 61071194 and 60979002, by the Fund of Civil Aviation University of China

For easy discussion, we assume that the structure of joint data vector **ss**( )*i* is shown in Fig.A.1, where circles represent SAR image pixels and *i* denotes the desired pixel pair

*SAR* 1 *SAR* 2

**ss**( )*i i i* 4 *i* 3 *i i* 1

12 2 22 ( ) [ ( ), ( 4), ( 3), ( ), ( 1)]*T* **ss** *i sisi si sisi* (A.1)

*i i*1

*i* 4 *i*3

**7. Acknowledgement**

**8.1 Proof of equation (12)**

**8. Appendix A**

under grant ZXH2009D018 and 2011kyE06.

whose interferometric phase is to be estimated.

*i*

Fig. A.1. Formulation of the joint data vector **ss**( )*i* .

The joint data vector **ss**( )*i* , shown in Fig.A.1, can be written as

The corresponding covariance matrix ( )*i* **Css** is given by

$$\mathbf{p}(\wp\_{i}) = \left[\mathbf{1}, \mathbf{e}^{j\wp\_{i}}, \mathbf{e}^{j\wp\_{i}}, \mathbf{e}^{j\wp\_{i}}, \mathbf{e}^{j\wp\_{i}}\right]^{T} \tag{A.3}$$

$$\mathbf{R}\_{\mathsf{so}}(i) = \sigma\_s^2(i) \begin{bmatrix} 1, \ r\_{12}(i, i-4), r\_{12}(i, i-3), r\_{12}(i, i), r\_{12}(i, i+1) \\ r\_{21}(i-4, i), 1, r\_{22}(i-4, i-3), r\_{22}(i-4, i), r\_{22}(i-4, i+1) \\ r\_{21}(i-3, i), r\_{22}(i-3, i-4), 1, r\_{22}(i-3, i), r\_{22}(i-3, i+1) \\ r\_{21}(i, i), r\_{22}(i, i-4), r\_{22}(i, i-3), 1, \ r\_{22}(i, i+1) \\ r\_{21}(i+1, i), r\_{22}(i+1, i-4), r\_{22}(i+1, i-3), r\_{22}(i+1, i), 1 \end{bmatrix} \tag{A.4}$$

So the covariance matrix ( )*i* **Csi** of the joint data vector **si**( )*i* , shown in Fig.5, can be given by

$$\begin{aligned} \mathbf{C}\_{\text{si}}(i) &= E\{\mathbf{si}(i)\mathbf{si}^{H}(i)\} \\ &= \mathbf{a}\mathbf{i}(\boldsymbol{\wp}\_{i})\mathbf{a}\mathbf{i}^{H}(\boldsymbol{\wp}\_{i}) \odot \mathbf{R}\_{\text{si}}(i) + \boldsymbol{\sigma}\_{n}^{2}\mathbf{I} \end{aligned} \tag{A.5}$$

### **9. Appendix B**

### **9.1 Fast algorithm for optimal interferometric phase estimation**

If *U* , *V* and *W* are arbitrary complex column vectors, then[20]

Robust Interferometric Phase Estimation in InSAR via Joint Subspace Projection 131

5 5 1 1 2 2

*iiii i*

()()

(( ) ) ( )

( )

*b b be e be e*

*mn n n*

*mn n n*

2 cos( )

*mn n i*

( 0) *i*

[1] P. A. Rosen, S. Hensley, I. R. Joughin, F. K. Li, S. N. Madsen, E. Rodrìguez and R. M.

[2] R. Bamler and P. Hartl. Synthetic aperture radar interferometry. Inverse Problem,

[3] Wei Xu and Ian Cumming. A Region-Growing Algorithm for InSAR Phase Unwrapping.

[4] M. Costantini. A Novel Phase Unwrapping Method based on Network Programming.

[5] R. M. Goldstein, H. A. Zebker, C. L. Werner.Satellite radar interferometry : two-

dimensional phase unwrapping. Radio Sci ., 1988,23(4):713 -720.

*mn n n*

*jjjj j*

*e e e e bbbbb e*

55 5 5 11 1 1 22 2 2 55 5 5 11 1 1 22 2 2 55 5 5 11 1 1 22 2 2

*b b be be*

*mn n n*

*mn n n*

*b b be be*

*mn n n*

*n n*

*n n*

*b be*

*j*

where

*bbbbb e bbbbb e*

*i i*

*j j*

*i i*

*j j*

*i i*

*i* is directly obtained by using the fast

*k* ( *k* is an integer).

(B.8)

*j j j j*

( 0)

Goldstein. Synthetic Aperture Radar Interferometry. In: Proceeding of the IEEE,

*bbbbb e bbbbb*

5 1 2 ( ) *n n*

 and

(B.7)

*angle b*

1

*j*

*i*

*i i*

*j j*

so we can get

Since and *i* , thus

**10. References**

 .

5 5 1 1 2 2 ( ) *n n n n b b*

3

*J*

*H*

**β Bβ**

*H*

, Let

Using (B.5) and (B.6), the cost function of (B.3) can be rewritten as 4

1

1, , , ,

55 5 11 1 22 2

*mn n*

*b bb*

Obviously, the minimization of 3*J* can be obtained for 2 *i*

algorithm, which can significantly reduce the computational burden.

IEEE Trans. On GRS, 1999, 37(1): 124-134.

IEEE Trans. On GRS, 1998, 36(3): 813-831.

So the closed-form solution to the estimate of

2000,88(3):333-382.

1998,14: R1-R54.

() ()

*i i*

**β A β**

() ()

*i ni n*

$$\left(\left(\mathcal{U}\odot V\right)^{H}\mathcal{W}\cdot\mathcal{W}^{H}\left(\mathcal{U}\odot V\right)=\mathcal{U}^{H}\left[\left(\mathcal{W}\cdot\mathcal{W}^{H}\right)\odot\left(V^{\star}\cdot\left(V^{\star}\right)^{H}\right)\right]\mathcal{U}\tag{8.1}$$

Using the equation (B.1), we can rewrite the cost function of (29) as

$$\begin{split} J\_{3} &= \sum\_{m=1}^{K} \sum\_{l=1}^{20-K} \left( \mathbf{a} \mathbf{i} (\boldsymbol{\phi}\_{l}) \odot \hat{\boldsymbol{\mathfrak{P}}}\_{rsi}^{(m)} \right)^{H} \hat{\boldsymbol{\mathfrak{P}}}\_{rsi}^{(l)} \hat{\boldsymbol{\mathfrak{P}}}\_{nsi}^{(l)H} (\mathbf{a} (\boldsymbol{\phi}\_{l}) \odot \hat{\boldsymbol{\mathfrak{P}}}\_{rsi}^{(m)}) \\ &= \sum\_{m=1}^{K} \sum\_{l=1}^{20-K} \left\{ \mathbf{a} \mathbf{i}^{H} \left( \boldsymbol{\phi}\_{l} \right) \right\} \left( \hat{\boldsymbol{\mathfrak{P}}}\_{nsi}^{(l)} \hat{\boldsymbol{\mathfrak{P}}}\_{nsi}^{(l)} \right) \odot \left( \hat{\boldsymbol{\mathfrak{P}}}\_{rsi}^{(m)\*} \left( \hat{\boldsymbol{\mathfrak{P}}}\_{rsi}^{(m)\*} \right)^{H} \right) \right\} \mathbf{a} \mathbf{i} (\boldsymbol{\phi}\_{l}) \right\} \\ &= \mathbf{a} \mathbf{i}^{H} (\boldsymbol{\phi}\_{l}) \left\{ \sum\_{m=1}^{K} \sum\_{l=1}^{20-K} \left[ (\hat{\boldsymbol{\mathfrak{P}}}\_{nsi}^{(l)} \hat{\boldsymbol{\mathfrak{P}}}\_{nsi}^{(l)H}) \odot (\hat{\boldsymbol{\mathfrak{P}}}\_{rsi}^{(m)\*} (\hat{\boldsymbol{\mathfrak{P}}}\_{rsi}^{(m)\*})^{H}) \right] \right\} \mathbf{a} \mathbf{i} (\boldsymbol{\phi}\_{l}) \end{split} \tag{B.2}$$

Let 20 () () ( ) ( ) 1 1 ˆˆ ˆ ˆ [( ) ( ( ) )] *K K ll m m H H nsi nsi rsi rsi m l* **A ββ β β** . It can be easily proved that **A** ( 20 20 ) is a

Hermitian matrix. Then (B.2) can be rewritten as

$$\begin{split} \mathcal{I}\_{3} &= \mathbf{a}\mathbf{i}^{H}(\boldsymbol{\phi}) \left\{ \sum\_{m=1}^{K} \sum\_{l=1}^{20-K} \left[ (\hat{\mathbf{p}}\_{mi}^{(l)} \hat{\mathbf{p}}\_{mi}^{(l)})^{H} \right] \odot (\hat{\mathbf{p}}\_{mi}^{(m)\*} (\hat{\mathbf{p}}\_{mi}^{(m)\*})^{H}) \right] \mathbf{a}(\boldsymbol{\phi}) \\ &= \mathbf{a}\mathbf{i}^{H}(\boldsymbol{\phi}\_{1}) \mathbf{A} \mathbf{a}(\boldsymbol{\phi}\_{1}) \\ &= \left[ \mathbf{\hat{p}}^{H}(\boldsymbol{\phi}\_{1}), \mathbf{\hat{p}}^{H}(\boldsymbol{\phi}\_{1}), \mathbf{\hat{p}}^{H}(\boldsymbol{\phi}\_{1}), \mathbf{\hat{p}}^{H}(\boldsymbol{\phi}\_{1}) \right] \begin{bmatrix} \mathbf{A}\_{1} & \mathbf{A}\_{2} \\ \mathbf{A}\_{3} & \mathbf{A}\_{4} \end{bmatrix} \begin{bmatrix} \mathbf{\hat{p}}(\boldsymbol{\phi}\_{1}) \\ \mathbf{\hat{p}}(\boldsymbol{\phi}\_{1}) \\ \mathbf{\hat{p}}(\boldsymbol{\phi}\_{1}) \\ \mathbf{\hat{p}}(\boldsymbol{\phi}\_{1}) \end{bmatrix} \\ &= \sum\_{n=1}^{4} \mathbf{\hat{p}}^{H}(\boldsymbol{\phi}\_{i}) \mathbf{A}\_{n} \mathbf{\hat{p}}(\boldsymbol{\phi}\_{i}) \\ &= \mathbf{\hat{p}}^{H}(\boldsymbol{\phi}\_{i}) \left\{ \sum\_{n=1}^{4} \mathbf{A}\_{n} \right\} \mathbf{\hat{p}}(\boldsymbol{\phi}\_{i}) \end{split} \tag{6.3}$$

where

$$\mathbf{f}(\boldsymbol{\phi}\_{i}) = \begin{bmatrix} \mathbf{1}\_{\prime}e^{j\boldsymbol{\phi}\_{i}} \ \boldsymbol{e}^{j\boldsymbol{\phi}\_{i}} \ \boldsymbol{e}^{j\boldsymbol{\phi}\_{i}} \ \boldsymbol{e}^{j\boldsymbol{\phi}\_{i}} \end{bmatrix}^{T} \tag{\text{B.4}}$$

Let

$$\mathbf{B} = \sum\_{n=1}^{4} \mathbf{A}\_{n} = \begin{bmatrix} b\_{11} \ b\_{12} \ b\_{13} \ b\_{14} \ b\_{15} \\ b\_{21} \ b\_{22} \ b\_{23} \ b\_{24} \ b\_{25} \\ b\_{31} \ b\_{32} \ b\_{33} \ b\_{34} \ b\_{35} \\ b\_{41} \ b\_{42} \ b\_{43} \ b\_{44} \ b\_{45} \\ b\_{51} \ b\_{52} \ b\_{53} \ b\_{54} \ b\_{55} \end{bmatrix} \tag{B.5}$$

It can be easily proved that **B** is a Hermitian matrix, i.e.,

$$b\_{1n}^{\*} = b\_{n1} \left( n = 2, 3, 4, 5 \right) \tag{8.6}$$

Using (B.5) and (B.6), the cost function of (B.3) can be rewritten as

4 3 1 11 12 13 14 15 21 22 23 24 25 31 32 33 34 35 41 42 43 44 45 51 52 53 54 55 () () () () 1 1, , , , *i iiii i i i H i ni n H i i j jjjj j j j J bbbbb e bbbbb e e e e bbbbb e bbbbb e bbbbb e* **β A β β Bβ** 55 5 5 11 1 1 22 2 2 55 5 5 11 1 1 22 2 2 55 5 5 11 1 1 22 2 2 55 5 11 1 22 2 ()() (( ) ) ( ) ( ) 2 cos( ) *i i i i i i j j mn n n mn n n j j mn n n mn n n j j j j mn n n mn n n mn n i mn n b b be be b b be be b b be e be e b bb* (B.7)

Obviously, the minimization of 3*J* can be obtained for 2 *i k* ( *k* is an integer). Since and *i* , thus

$$\phi\_l = \begin{cases} -\pi - \mu & \text{ ( $\mu \le 0$ )}\\ \pi - \mu & \text{ ( $\mu > 0$ )} \end{cases} \tag{8.8}$$

So the closed-form solution to the estimate of *i* is directly obtained by using the fast algorithm, which can significantly reduce the computational burden.

### **10. References**

130 Recent Interferometry Applications in Topography and Astronomy

Using the equation (B.1), we can rewrite the cost function of (29) as

20

*K K*

ˆˆ ˆ ˆ ( ) [( ) ( (

**ai ββ β β**

*H H ll m m i nsi nsi rsi rsi*

1 1

20

*K K*

1 1

*iiii*

*m l*

*m l*

() () ( ) ( )

() ()

*i i*

() ()

*ini*

1

*n*

It can be easily proved that **B** is a Hermitian matrix, i.e.,

*n*

ˆˆ ˆ ˆ [( ) ( ( ) )]

*ll m m H H nsi nsi rsi rsi*

20

*K K*

1 1 20

*m l K K*

*m l*

Hermitian matrix. Then (B.2) can be rewritten as

*H*

4

1

*n*

**β**

*H*

**β A β**

4

1 () () *H i ni n*

 **A β**

1 1

**ai Aai**

3

3

*J*

20

*K K*

1 1

*m l*

Let

where

Let

*J*

( ) ( ) ( ) ( ( )) *HH H H H U V WW U V U WW V V U* (B.1)

() () ( ) (

() () ( ) ( )

ˆˆ ˆ ˆ ( ) [( ) ( ( ) )] ( )

*i i nsi nsi rsi rsi*

( ) ( ), ( ), ( ), ( ) ( )

*HHHH i*

( ) 1, , , , *iiii*

*i eeee*

*T jjjj*

11 12 13 14 15

*bbbbb bbbbb bbbbb bbbbb bbbbb*

31 32 33 34 35

41 42 43 44 45 51 52 53 54 55

21 22 23 24 25 4

**A A β ββββ A A β**

*H HH ll m m*

**ai ββ β β ai**

1 2 3 4

() () ( ) ( )

ˆ ˆˆ ˆ (() ) (() )

**ai β ββ ai β**

*HH H ll m m*

**A ββ β β** . It can be easily proved that **A** ( 20 20 ) is a

**ai ββ β β ai**

( ) () () ( )

*m l lH H m i i rsi nsi nsi rsi*

>

> > **ai**

*i*

(B.2)

(B.3)

 ) ) )] ( ) *H*

( )

**β**

*i*

( )

**β**

**β** (B.4)

**B A** (B.5)

1 1 *n n b b* ( 2,3,4,5 *n* ) (B.6)

*i i*

ˆˆ ˆ ˆ ( )[( ) ( ( ) )] ( )

*i i nsi nsi rsi rsi*

**7**

Tao Yu

*P. R. China*

**Airborne Passive Localization Method Based**

Airborne passive localization with single station is a key technique which is quickly developing lately. If the exact localization to target can be realized in single mobile platform, its significance is self-evident. Because the information content obtained by single station is

Phase interferometry is mainly in use for direction finding (DF). Moreover, Doppler shift may be used to determine the coordinate setting of target. As far as actual applications are concerned, these two methods are independently applied. Or, they firstly are respectively

Author researches the functional relation between Doppler shift and phase interferometry and presents some method which can syncretize these two methods directly from the equation of mathematical. Based on this fusion in physical relationship, we can obtain the Doppler shift by measuring phase. Or, we can determine the phase by measuring Doppler shift. This book describes these new methods as well as some applications in airborne

Section two describes the relationship between phase and Doppler shift. Firstly, the phase difference localization equation which is analogous to time difference localization equation is introduced by analyzing the phase in radial distance of target. And then, making use of the complementary relationship between the angle of advance at airborne platform and the angle of arrival for target, the approximate functional relation between Phase-difference and Doppler shift can be derived by combining DF formula based on phase interference and expression of Doppler shift. Where after, the calculating method of airborne Doppler rate is

Section three researches the problem of airborne Doppler direct range finding. Three solution methods are presented. First method is direct approximation for mathematical expression of Doppler rate. Second method is introducing differential transformation for Doppler shift equation based on angle rate. Third method is based on the specific value of Doppler rate between two adjacent detection nodes. On this basis, the airborne direct range finding formula based on phase measuring can be obtained by using the relationship

researched based on phase difference or Doppler frequency difference

between phase difference and Doppler shift.

lesser, the difficulty would be obviously bigger relative to multi-station observation.

detected. And then, the obtaining data are syncretized.

**1. Introduction**

passive localization system.

**on Doppler-Phase Interference Measuring**

*Shanghai Research Institute of Microwave Equipment, Shanghai*

Tao Yu

*Shanghai Research Institute of Microwave Equipment, Shanghai P. R. China*

### **1. Introduction**

132 Recent Interferometry Applications in Topography and Astronomy

[6] M. D. Pritt and J. S. Shipman.Least-Squares Two-Dimensional Phase Unwrapping Using

[7] B. Zitova and J. Flusser. Image registration methods: A survey. Image Vis. Comput.,

[8] A. Cole-Rhodes, K. L. Johnson, J. Le Moigne, and I. Zavorin. Multiresolution registration

[9] X. Dai and S. Khorram. A feature-based image registration algorithm using improved

[10] Harold S. Stone, Michael T. Orchard, Ee-Chien Chang, and Stephen A. Martucci. A

[11] Eugenio Sansosti, Paolo Berardino, Michele Manunta, Francesco Serafino, and

[12] Francesco Serafino. SAR Image Coregistration Based on Isolated Point Scatterers. IEEE

[13] P. H. Eichel, D. C. Ghiglia, et al. Spotlight SAR Interferometry for Terrain Elevation

[14] R. Lanari, G. Fornaro, et al. Generation of Digital Elevation Models by Using SIR-C/X-

[15] Jong-Sen Lee, Konstantinos P. Papathanassiou, et al. A New Technique for Noise

[16] Ireneusz Baran, Mike P. Stewart, Bert M. Kampes, Zbigniew Perski, and Peter Lilly. A

[17] Nan Wu, Da-Zheng Feng, and Junxia Li. A locally adaptive filter of interferometric

[18] Qifeng Yu, Xia Yang, Sihua Fu, Xiaolin Liu, and Xiangyi Sun. An Adaptive Contoured

[19] Hai Li, Guisheng Liao. An estimation method for InSAR interferometric phase based

[20] Hai Li, Zhenfang Li, Guisheng Liao, and Zheng Bao. An estimation method for InSAR

[21] Zhenfang Li, Zheng Bao, Hai Li, and Guisheng Liao. Image Auto-Coregistration and

[22] F. Lombardini*,* M. Montanari, and F. Gini. Reflectivity Estimation for Multibaseline

[23] F. Lombardini. Absolute Phase Retrieval in a Three-element Synthetic Aperture Radar

on MMSE criterion. IEEE Trans On GRS , 2010,48(3):1457-1469.

stochastic gradient. IEEE Trans. Image Process., 2003,12:1495-1511.

of remote sensing imagery by optimization of mutual information using a

chain-code representation combined with invariant moments. IEEE Trans. On GRS,

Fast Direct Fourier-Based Algorithm for Subpixel Registration of Images. IEEE

Gianfranco Fornaro. Geometrical SAR Image Registration. IEEE Trans. On GRS,

Mapping and Interferometric Change Detection. Sandia National Labs Tech.

SAR Multifrequency Two-Pass Interferometry: The Etna Case Study. IEEE Trans.

Filtering of SAR Interferometric Phase Images. IEEE Trans. On GRS, 1998,

Modification to the Goldstein Radar Interferogram Filter. IEEE Trans. On GRS,

Window Filter for Interferometric Synthetic Aperture Radar . IEEE GRS Letters,

interferometric phase combined with image auto-coregistration. Science in China,

InSAR Interferogram Estimation Using Joint Subspace Projection. IEEE Trans. On

Interferometric Radar Imaging of Layover Extended Sources. IEEE Trans On SP,

Interferometer. In: Proceeding of the 1996 CIE Int. Conf. of Radar, Beijing, China,

FFT's. IEEE Trans. On GRS ,1994, 32(3):706-708.

2003,21(11): 977-1000.

1999, 37(5):2351-2362.

2006, 44(10):2861-2870.

GRS Letters, 2006, 3(3):354-358.

On GRS, 1996, 34(5):1097-1114.

36(5):1456-1465.

2007, 4(1):23-26.

2003, 41(9):2114-2118.

Series F, 2006, 49(3): 386-396.

GRS ,2006,44(2):288-297.

2003, 51(6): 1508-1519.

1996, 309-312.

Trans. On GRS, 2001, 39(10):2235-2243.

Report, SAND93, December 1993, pp. 2539-2546.

phase images. IEEE GRS Letters, 2006, 3(1):73-77.

Airborne passive localization with single station is a key technique which is quickly developing lately. If the exact localization to target can be realized in single mobile platform, its significance is self-evident. Because the information content obtained by single station is lesser, the difficulty would be obviously bigger relative to multi-station observation.

Phase interferometry is mainly in use for direction finding (DF). Moreover, Doppler shift may be used to determine the coordinate setting of target. As far as actual applications are concerned, these two methods are independently applied. Or, they firstly are respectively detected. And then, the obtaining data are syncretized.

Author researches the functional relation between Doppler shift and phase interferometry and presents some method which can syncretize these two methods directly from the equation of mathematical. Based on this fusion in physical relationship, we can obtain the Doppler shift by measuring phase. Or, we can determine the phase by measuring Doppler shift. This book describes these new methods as well as some applications in airborne passive localization system.

Section two describes the relationship between phase and Doppler shift. Firstly, the phase difference localization equation which is analogous to time difference localization equation is introduced by analyzing the phase in radial distance of target. And then, making use of the complementary relationship between the angle of advance at airborne platform and the angle of arrival for target, the approximate functional relation between Phase-difference and Doppler shift can be derived by combining DF formula based on phase interference and expression of Doppler shift. Where after, the calculating method of airborne Doppler rate is researched based on phase difference or Doppler frequency difference

Section three researches the problem of airborne Doppler direct range finding. Three solution methods are presented. First method is direct approximation for mathematical expression of Doppler rate. Second method is introducing differential transformation for Doppler shift equation based on angle rate. Third method is based on the specific value of Doppler rate between two adjacent detection nodes. On this basis, the airborne direct range finding formula based on phase measuring can be obtained by using the relationship between phase difference and Doppler shift.

Thus, the path length difference *r* of radial distance between two antenna elements can be determined by measuring for phase difference. The localization equation based on phase difference can be obtained whose formed expression is completely analogous to the

1 2

(2-2)

(2-3)

1 2 ( ) 2

( ) 2

where: *NN N* 1 2 is the whole number of wavelength in path length difference;

Provided that incident wave from the some radiant is approximatively regarded as plane wave, the existing DF expression based on the principle of phase interferometry can be

( ) 2

*λ*

*d N* (2-4)

Or, by simply rearranging, the phase difference between two adjacent antennas can take the

2 2 sin

is arrival angle of target signal; *d* baseline length between two antennas.

Making use of the complementary relationship between the angle of advance and the angle of arrival, the approximate functional relation between phase-difference and Doppler shift can be obtained by combining the DF formula based on phase interference and the expression of Doppler shift. The analog calculation certificates the correctness of derived formula. The error analysis presents the measurement accuracy of Doppler shift obtained by phase-difference measuring. And calculation shows that the measuring error of Doppler shift obtained based on phase-difference measuring is lower if the detection error for integer

*N d*

*N*

*N N*

1 2

*rr r*

sin

 *r d*

localization equation based on time difference

1 2 phase difference.

approximatively obtained due to sine theorem

**2.2 Phase difference detecting of Doppler shift**

of wavelength in path difference is not considered.

form

where:

**2.2.1 Recapitulation**

In allusion to the problem to need solving integer ambiguity in available DF system based on phase interferometer, section four puts forward a DF method combining Doppler shift and phase interferometer. According to this method, the integer of wavelength in radial distant can be directly solved synthetically by making use of speed vector and Doppler shift as well as rate. Hence, the integer of wavelength in path length difference of radial distant between two adjacent antenna elements can be obtained. At the same time, the value less than one integer in path difference can be determined by using phase interferometer. As compared with available phase interferometer which firstly solves phase difference, this new method which firstly solves path difference has not phase ambiguity, nor it needs to limit baseline length.

Section fifthly presents a new airborne passive DF method using two orthogonal baselines based on the rate of direction cosine. This is a passive DF method which is only associated with Doppler frequency difference and is not associated with wavelength.

### **2. Relationship between phase and Doppler shift**

### **2.1 Phase-detecting for distance**

For single baseline interferometer as shown in fig.2-1, if the phase shift measured by descriminator, corresponding to the radial distance *ir* , is *i* , the expression of radial distance based on phase measuring is

$$r\_i = \mathcal{A}(N\_i + \frac{\phi\_i}{2\pi}) \tag{2-1}$$

where: is wavelength; *Ni* the whole number of wavelength.

Fig. 2-1. Single baseline array on airborne platform.

Thus, the path length difference *r* of radial distance between two antenna elements can be determined by measuring for phase difference. The localization equation based on phase difference can be obtained whose formed expression is completely analogous to the localization equation based on time difference

$$
\Delta r = r\_1 - r\_2
$$

$$
= \mathcal{A}(N\_1 - N\_2 + \frac{\phi\_1 - \phi\_2}{2\pi})
\tag{2-2}
$$

$$
= \mathcal{A}(\Delta N + \frac{\Delta \phi}{2\pi})
$$

where: *NN N* 1 2 is the whole number of wavelength in path length difference; 1 2 phase difference.

Provided that incident wave from the some radiant is approximatively regarded as plane wave, the existing DF expression based on the principle of phase interferometry can be approximatively obtained due to sine theorem

$$\begin{aligned} \sin \theta &= \frac{\Delta r}{d} \\\\ &= \frac{\mathcal{L}}{d} (\Delta N + \frac{\Delta \phi}{2\pi}) \end{aligned} \tag{2-3}$$

Or, by simply rearranging, the phase difference between two adjacent antennas can take the form

$$2\pi \Delta \text{N} + \Delta \phi = \frac{2\pi d}{\lambda} \sin \theta \tag{2-4}$$

where: is arrival angle of target signal; *d* baseline length between two antennas.

### **2.2 Phase difference detecting of Doppler shift**

### **2.2.1 Recapitulation**

134 Recent Interferometry Applications in Topography and Astronomy

In allusion to the problem to need solving integer ambiguity in available DF system based on phase interferometer, section four puts forward a DF method combining Doppler shift and phase interferometer. According to this method, the integer of wavelength in radial distant can be directly solved synthetically by making use of speed vector and Doppler shift as well as rate. Hence, the integer of wavelength in path length difference of radial distant between two adjacent antenna elements can be obtained. At the same time, the value less than one integer in path difference can be determined by using phase interferometer. As compared with available phase interferometer which firstly solves phase difference, this new method which firstly solves path difference has not phase ambiguity, nor it needs to

Section fifthly presents a new airborne passive DF method using two orthogonal baselines based on the rate of direction cosine. This is a passive DF method which is only associated

For single baseline interferometer as shown in fig.2-1, if the phase shift measured by

( ) 2

*y Txy* (, )

*r*

1 *d* 2 *x*

1 Flight direction of airborne platform

*i*

*i i r N* (2-1)

*i* , the expression of radial distance

with Doppler frequency difference and is not associated with wavelength.

is wavelength; *Ni* the whole number of wavelength.

2*r*

**2. Relationship between phase and Doppler shift**

descriminator, corresponding to the radial distance *ir* , is

1

Fig. 2-1. Single baseline array on airborne platform.

limit baseline length.

**2.1 Phase-detecting for distance**

based on phase measuring is

where:

*r*

Making use of the complementary relationship between the angle of advance and the angle of arrival, the approximate functional relation between phase-difference and Doppler shift can be obtained by combining the DF formula based on phase interference and the expression of Doppler shift. The analog calculation certificates the correctness of derived formula. The error analysis presents the measurement accuracy of Doppler shift obtained by phase-difference measuring. And calculation shows that the measuring error of Doppler shift obtained based on phase-difference measuring is lower if the detection error for integer of wavelength in path difference is not considered.

*T*

 1 *d* 2 *d* Flight direction of airborne platform 1 2 3 *x*

In order to verify the calculating formula of Doppler shift based on phase difference measuring, we make the analog calculation by replacing measured value with theoretical

in prescribed domain. Hence, the rest radial distance as well as path difference can be obtained according to circular function relationship. Simultaneity, the theoretical value of Doppler shift corresponding to every radial distance can be calculated by Doppler shift

On this basis, the Doppler shift based on phase measuring is computed by Eq.(2-7) and then

0

*f* (2-9)

0 100

, baseline

change

value. Firstly, we preset following parameter: radial distance 1*r* , wavelength

the relative calculation error is obtained by compared with the theoretical value

*di dai*

*di f f*

Without the notice, the adopted parameter is as follows: 1*r km* 100 , *v ms* 100 / ,

Fig.2-3 depicts the calculating error curve of Doppler shift *d*1 *f* when the baseline length is different. The mathematical simulation shows that calculating error is independent of flight

speed and directly related to wavelength and inversely related to radial distance.

where: subscript *a* expresses calculating value.

0.0375 *m* .

*f*

length *d* , the flight speed of flight device *v* . And then, we make the arrival angle

*y*

1*r*

Fig. 2-2. One-dimensional array with double baseline.

2*r* 3*r*

**2.2.3 Computational error**

formula.

*d* 5 ,  ### **2.2.2 Primitive formula**

As shown in fig.2-1, a single baseline array with two antenna elements is installed on airborne platform and the spacing of array is *d* . The direction of axis of baseline is parallel the axis of airborne platform. For the target *T* at stationary or low speed, the Doppler shift detected by airborne double channel measuring receiver in every element of single baseline interferometry is

$$
\lambda f\_{d1} = v \cos \beta\_1 \tag{2-5}
$$

where: *d*1 *f* is Doppler shift; *v* the flight speed of detection platform; the angle of advance between radial distance and travelling direction of flight device.

According to the geometric relationship as shown fig.2-1 and the analysis results in a previous section, results in

$$\begin{aligned} \cos \beta\_1 &= \sin(90^\circ - \beta\_1) \\\\ &= \sin \theta \approx \frac{\Delta r}{d} \\\\ &= \frac{\lambda}{d} (\Delta N + \frac{\Delta \phi}{2\pi}) \end{aligned} \tag{2-6}$$

Substituting (2-6) into the Doppler shift expression, we can obtain the calculating formula of Doppler shift based on phase difference measuring

$$\begin{split} f\_{d1} &= \frac{v}{d} (\Delta \mathbf{N} + \frac{\Delta \phi}{2\pi}) \\\\ &= \frac{v \Delta r}{d \lambda} \end{split} \tag{2-7}$$

According to above-mentioned derivation, we can directly written out the formula of Doppler frequency difference based on phase difference measuring by use of the array with double baseline in one-dimensional as shown in fig.2-2

$$\begin{aligned} \Delta f\_d &= \frac{\upsilon}{d} \Big[ (\Delta N\_1 + \frac{\Delta \phi\_1}{2\pi}) - (\Delta N\_2 + \frac{\Delta \phi\_2}{2\pi}) \Big] \\\\ &= \frac{\upsilon}{d} \Big[ (\Delta N\_1 - \Delta N\_2) + \frac{1}{2\pi} (\Delta \phi\_1 - \Delta \phi\_2) \Big] \end{aligned} \tag{2-8}$$

where: *NNN iii*1 ; *iii* 1 .

Fig. 2-2. One-dimensional array with double baseline.

### **2.2.3 Computational error**

136 Recent Interferometry Applications in Topography and Astronomy

As shown in fig.2-1, a single baseline array with two antenna elements is installed on airborne platform and the spacing of array is *d* . The direction of axis of baseline is parallel the axis of airborne platform. For the target *T* at stationary or low speed, the Doppler shift detected by airborne double channel measuring receiver in every element of single baseline

1 1

According to the geometric relationship as shown fig.2-1 and the analysis results in a

sin

Substituting (2-6) into the Doppler shift expression, we can obtain the calculating formula of

1 ( ) 2

*d v f N d*

According to above-mentioned derivation, we can directly written out the formula of Doppler frequency difference based on phase difference measuring by use of the array with

1 2

*v fN N d*

*v N N*

( )( ) 2 2

*v r d*

*N d*

0 1 1 cos sin(90 )

*r d*

1 2

>

> >

1 2 12

1 ( )( ) 2

( ) 2

(2-5)

the angle of

(2-6)

(2-7)

(2-8)

*f v d* cos

where: *d*1 *f* is Doppler shift; *v* the flight speed of detection platform;

advance between radial distance and travelling direction of flight device.

**2.2.2 Primitive formula**

previous section, results in

where: *NNN iii*1 ;

Doppler shift based on phase difference measuring

double baseline in one-dimensional as shown in fig.2-2

*iii* 1 .

*d*

*d*

interferometry is

In order to verify the calculating formula of Doppler shift based on phase difference measuring, we make the analog calculation by replacing measured value with theoretical value. Firstly, we preset following parameter: radial distance 1*r* , wavelength , baseline length *d* , the flight speed of flight device *v* . And then, we make the arrival angle change in prescribed domain. Hence, the rest radial distance as well as path difference can be obtained according to circular function relationship. Simultaneity, the theoretical value of Doppler shift corresponding to every radial distance can be calculated by Doppler shift formula.

On this basis, the Doppler shift based on phase measuring is computed by Eq.(2-7) and then the relative calculation error is obtained by compared with the theoretical value

$$
\omega\_f = \left| \frac{f\_{di} - f\_{dai}}{f\_{di}} \right| \times 100 \,\% \tag{2-9}
$$

where: subscript *a* expresses calculating value.

Without the notice, the adopted parameter is as follows: 1*r km* 100 , *v ms* 100 / , *d* 5 , 0.0375 *m* .

Fig.2-3 depicts the calculating error curve of Doppler shift *d*1 *f* when the baseline length is different. The mathematical simulation shows that calculating error is independent of flight speed and directly related to wavelength and inversely related to radial distance.

wavelength in path difference has error. According to the same reason, the larger measuring error of Doppler shift is also produced by the error of speed measuring in axle direction of

Fig.2-4 depicts the Doppler measuring error with different baseline length. The analysis shows that the measuring error in integer of wavelength has prodigious influence for measuring error of Doppler shift. Also, the change for baseline length is also considerable

Fig.2-5 illustrates the change curve versus wavelength for Doppler shift error when the integer of wavelength in path difference can be accurately detected. And the analyzer shows that the measuring error is not associated baseline and the influence produced by change of

0 10 20 30 40 50 60 70 80 90

d = 5 λ d = 10 λ d = 15 λ

Angle of arrival θ / (°)

baseline.

sensitive.

flight speed is also lesser.

150

Fig. 2-4. Measurement error of Doppler shift.

200

250

300

350

400

( )

σ / Hz

Error analysis

450

500

550

Fig. 2-3. Computational error of Doppler shift with different baseline lengths.

### **2.2.4 Precision analysis**

According to error estimation theory, the measuring error of Doppler shift produced by the measuring error of phase difference, flight speed and integer of wavelength in path difference is

$$\begin{split} \sigma &= \sqrt{\left(\frac{\partial f\_{d1}}{\partial \Delta \phi} \sigma\_{\phi}\right)^{2} + \left(\frac{\partial f\_{d1}}{\partial v} \sigma\_{v}\right)^{2} + \left(\frac{\partial f\_{d1}}{\partial \Delta N} \sigma\_{n}\right)^{2}} \\ &= \sqrt{\left(\frac{v}{2\pi d} \sigma\_{\phi}\right)^{2} + \left(\frac{\Delta r}{\lambda d} \sigma\_{v}\right)^{2} + \left(\frac{v}{d} \sigma\_{n}\right)^{2}} \end{split} \tag{2-10}$$

where: , *v* and *n* is respectively the mean-root-square error measuring phase difference, flight speed and integer of wavelength in path difference. In where, unit for is radian. Without the notice, the value of mean-root-square error is respectively as follows: 0 0 180 ( 1 ) , *v* 1 / *m s* , 1 *n* .

It can be seen from (2-10) that the Doppler measuring error produced only by measuring error of phase difference is just a constant term. If 0 180 (radian), value will be less than 2*Hz* making use of the parameter presented in previous section. The Doppler measuring error produced only by integer of wavelength in path difference is also a constant term. Because the speed is usually larger than baseline length for Airborne applications, the Doppler measuring error will be prodigious if the measuring for integer of

d = 5 λ d = 10 λ d = 15 λ

0 10 20 30 40 50 60 70 80 90

Angle of arrival θ / (°)

According to error estimation theory, the measuring error of Doppler shift produced by the measuring error of phase difference, flight speed and integer of wavelength in path

11 1

*ff f*

*v rv d dd*

difference, flight speed and integer of wavelength in path difference. In where, unit for

is radian. Without the notice, the value of mean-root-square error is respectively as follows:

It can be seen from (2-10) that the Doppler measuring error produced only by measuring

than 2*Hz* making use of the parameter presented in previous section. The Doppler measuring error produced only by integer of wavelength in path difference is also a constant term. Because the speed is usually larger than baseline length for Airborne applications, the Doppler measuring error will be prodigious if the measuring for integer of

*dd d v n*

2 2 2

*v N*

(radian),

(2-10)

value will be less

2 22

*v n*

*n* is respectively the mean-root-square error measuring phase

Fig. 2-3. Computational error of Doppler shift with different baseline lengths.

2

1 *n* .

error of phase difference is just a constant term. If 0 180

*v* 1 / *m s* ,

10-8

**2.2.4 Precision analysis**

difference is

where:

 , *v* and

0 0 180 ( 1 )

,

10-7

10-6

10-5

Computational error ε

f ( )/ %

10-4

10-3

10-2

10-1

wavelength in path difference has error. According to the same reason, the larger measuring error of Doppler shift is also produced by the error of speed measuring in axle direction of baseline.

Fig.2-4 depicts the Doppler measuring error with different baseline length. The analysis shows that the measuring error in integer of wavelength has prodigious influence for measuring error of Doppler shift. Also, the change for baseline length is also considerable sensitive.

Fig.2-5 illustrates the change curve versus wavelength for Doppler shift error when the integer of wavelength in path difference can be accurately detected. And the analyzer shows that the measuring error is not associated baseline and the influence produced by change of flight speed is also lesser.

Fig. 2-4. Measurement error of Doppler shift.

d = 5 λ d = 10 λ d = 15 λ

0 10 20 30 40 50 60 70 80 90

Angle of arrival θ / (°)

According to localization theory, the radial distance between measuring platform and measured target can be directly obtained based on Doppler changing rate equation. But, in fact, the localization method based on Doppler rate of change is not the classical method in current target localization for electronic warfare. A main reason is that measuring for change rate of Doppler shift is also relative difficult at present. At the same time, because Doppler rate of change is directly concerned with the tangential velocity, the problem related to direct range finding using Doppler rate equation is that the detecting system has to measure the angle of advance between the traveling direction of detection platform and the radial distance to target, along with measuring Doppler rate. Thereout, the localization equation can be solved. Hence, the localization method based on Doppler rate of change can be

In addition, for localization and tracker system, it is extraordinary valuable to obtain the Doppler changing rate of received signal in order to estimate the state and position of target. At present, Doppler rate of change is mainly obtained by estimating the frequency variation of received signal based on the principle that Doppler rate of change is the same in mathematical analysis as the rate of carrier frequency of received signal. These estimative algorithms are not only associated with modulation mode, but they are also complicated.

**2.3 Airborne Doppler changing rate obtaining by frequency difference or phase**

1

Fig. 2-6. Relative measurement error of Doppler shift.

completed with other localization method at present.

1.5

2

2.5

Relative measurement error

**difference**

**2.3.1 Recapitulation**

σr ( )/ %

3

3.5

4

4.5

Fig. 2-5. Measurement error of Doppler shift without integer error.

More, relative calculation error of Doppler shift is analyzed by perfect differential, that is

$$df\_{d1} = \frac{\partial f\_{d1}}{\partial \Delta \phi} d\Delta \phi\_i + \frac{\partial f\_{d1}}{\partial v} dv + \frac{\partial f\_{d1}}{\partial \Delta N} d\Delta N \tag{2-11}$$

When the error of observing variable is zero-mean and unaided reciprocally, the relative calculation error of Doppler shift is

$$\begin{aligned} \boldsymbol{\sigma}\_{r} &= \left| \frac{d f\_{di}}{f\_{di}} \right| \\\\ &= \frac{1}{f\_{di}} \left( \left| \frac{\partial f\_{di}}{\partial \Delta \phi\_{i}} \boldsymbol{\sigma}\_{\phi} \right| + \left| \frac{\partial f\_{d1}}{\partial \boldsymbol{v}} \boldsymbol{\sigma}\_{\boldsymbol{v}} \right| + \left| \frac{\partial f\_{d1}}{\partial \Delta N} \boldsymbol{\sigma}\_{\boldsymbol{n}} \right| \right) \\\\ &= \frac{\lambda d}{v \Delta r} \left( \left| \frac{\boldsymbol{v}}{2 \pi d} \boldsymbol{\sigma}\_{\phi} \right| + \left| \frac{\Delta r}{\lambda d} \boldsymbol{\sigma}\_{\boldsymbol{v}} \right| + \left| \frac{\boldsymbol{v}}{d} \boldsymbol{\sigma}\_{\boldsymbol{n}} \right| \right) \end{aligned} \tag{2-12}$$

The advantage making the expression which contains phase difference transform the expression which includes path difference *r* is that the analysis and calculation for integer and phase difference can be avoided. Fig.2-6 illustrates the change curve versus baseline length for relative calculation error of Doppler shift when the measuring for integer is without error. The calculation shows that the curvilinear change is basically not associated with flight speed and wavelength. The relative calculation error is inversely relative to baseline length. It can be seen that the relative calculation error of Doppler shift is less than 1.5% after the angle of arrival is bigger than 0 10

Fig. 2-6. Relative measurement error of Doppler shift.

### **2.3 Airborne Doppler changing rate obtaining by frequency difference or phase difference**

### **2.3.1 Recapitulation**

140 Recent Interferometry Applications in Topography and Astronomy

0 10 20 30 40 50 60 70 80 90

Angle of arrival θ / (°)

1 11

*f ff df d dv d N v N* (2-11)

  (2-12)

*d dd*

When the error of observing variable is zero-mean and unaided reciprocally, the relative

1 1 1

*dv r v vr d d d*

The advantage making the expression which contains phase difference transform the expression which includes path difference *r* is that the analysis and calculation for integer and phase difference can be avoided. Fig.2-6 illustrates the change curve versus baseline length for relative calculation error of Doppler shift when the measuring for integer is without error. The calculation shows that the curvilinear change is basically not associated with flight speed and wavelength. The relative calculation error is inversely relative to baseline length. It can be seen that the relative calculation error of Doppler shift is less than

*fff f vN*

*di d d v n*

>

*v n*

More, relative calculation error of Doppler shift is analyzed by perfect differential, that is

0

calculation error of Doppler shift is

5

10

15

20

( )

Error analysis σ / Hz

25

30

λ = 0.0375 m λ = 0.06 m λ = 0.3 m

Fig. 2-5. Measurement error of Doppler shift without integer error.

1

*di r di*

1.5% after the angle of arrival is bigger than 0 10

*df f*

*d i*

2

*di i*

According to localization theory, the radial distance between measuring platform and measured target can be directly obtained based on Doppler changing rate equation. But, in fact, the localization method based on Doppler rate of change is not the classical method in current target localization for electronic warfare. A main reason is that measuring for change rate of Doppler shift is also relative difficult at present. At the same time, because Doppler rate of change is directly concerned with the tangential velocity, the problem related to direct range finding using Doppler rate equation is that the detecting system has to measure the angle of advance between the traveling direction of detection platform and the radial distance to target, along with measuring Doppler rate. Thereout, the localization equation can be solved. Hence, the localization method based on Doppler rate of change can be completed with other localization method at present.

In addition, for localization and tracker system, it is extraordinary valuable to obtain the Doppler changing rate of received signal in order to estimate the state and position of target. At present, Doppler rate of change is mainly obtained by estimating the frequency variation of received signal based on the principle that Doppler rate of change is the same in mathematical analysis as the rate of carrier frequency of received signal. These estimative algorithms are not only associated with modulation mode, but they are also complicated.

1 1 sin

In fact, when the flight device is uniform motion, the left-hand component of equation is

2 2

 *d v f*

So we obtain the calculation formula of Doppler rate only based on measuring Doppler

*f f d d t*

1 *v d t* It can be seen that the time variation of detecting Doppler shift can be equivalently expressed by specific value between flight distance and flight speed of detection platform. Hence, we also prove that Doppler changing rate is only associated with flight speed and it

Again, the frequency difference of Doppler shift can is also obtained from one of radiation

1 2 1 2

 *dd d*

So Doppler changing rate can be determined by the actual measurement to frequency

As soon as we substitute the expression (2-8) of Doppler frequency difference based on phase difference measuring into (2-18), we can obtain the formula of Doppler changing rate

> 2 1 21 1 ( )( ) 2

*v f NN d* (2-19)

*tt t*

*ff f*

2

*d*

*ff f*

sin

1

*r*

*d*

*v v f r d* (2-17)

*d d v f f d* (2-18)

2 2

is tangential velocity.

In where: 1 *v v t* sin

frequency difference

There results:

frequency

difference of signal.

based on phase difference

namely basic expression of Doppler rate of change

According to the mathematical definition of Doppler rate

/

is wholly not associated with light speed.

where: *ti f* is measured value of signal.

On the basis of researching Doppler direct ranging, this section presents a method determining Doppler changing rate only by detecting Doppler shift value of received signal at some discrete node. Here, an analytic method solving Doppler changing rate only by detecting the frequency difference of Doppler shift or received signal is presented based on azimuth rate. This new method is not only straightforward in detection mode; it is also succinct on expression. More important characteristic is lying in:

### **2.3.2 Derivation**

According to the geometric relationship as shown in fig.2-1, we have as following the approximate representation of circular function

$$
\sin \beta\_1 \approx \frac{\sqrt{d^2 - \Delta r^2}}{d} \tag{2-13}
$$

$$