**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

<sup>1</sup> (( ) ) (( ) ) *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 ( ) *<sup>s</sup>* **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 **β***<sup>n</sup>* , 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, ]*<sup>i</sup> <sup>j</sup> <sup>T</sup> <sup>i</sup> e* **a** is the spatial steering vector (i.e., the array steering vector) of the pixel *i* , superscript *T* denotes the vector transpose operation, *<sup>i</sup>* 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 ( ) *<sup>s</sup>* **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, <sup>2</sup> ( ) *<sup>s</sup> i* is the echo power of the pixel *i* and <sup>2</sup> *<sup>n</sup>* 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 ( ) *<sup>s</sup>* **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 ( ) *<sup>s</sup>* **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 ( ) *<sup>s</sup>* **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 **β***<sup>n</sup>* , 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 ( ), ( ) ( ) ( ), ( ) *<sup>s</sup> 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 ( ) *<sup>s</sup>* **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 <sup>s</sup>** *i*

 

*H s i is n*

*H H i i n*

 

2

**R** (3)

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

 

 

*<sup>i</sup>* is the terrain

(2)

*<sup>n</sup>* 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 <sup>2</sup>

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

decomposition of the covariance matrix ( ) *<sup>s</sup>* **C** *i* is as follows:

 

*i Ei i*

**C ss <sup>s</sup>**

*H*

 

**2. Data model and problem formulation** 

images can be formulated as follows[19],

corresponding covariance matrix ( ) *<sup>s</sup>* **C** *i* is given by

where ( ) [1, ]*<sup>i</sup> <sup>j</sup> <sup>T</sup> <sup>i</sup> e*

transpose, <sup>2</sup>

( ) *<sup>s</sup>* 

 **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*

*<sup>s</sup> 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 ( )*<sup>i</sup>* **Csi** , ( ) *<sup>k</sup>*

*rsi* 

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

 

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

of ( )*<sup>i</sup>* **Rsi** . <sup>2</sup>

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

*<sup>n</sup>* 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 ( )*<sup>i</sup>* **Csi** , ( )*<sup>l</sup>*

where **Rs** ( ) *<sup>m</sup>* and <sup>2</sup>

*<sup>i</sup> rsi* **ai β** are orthogonal to ( )*<sup>l</sup>*

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

( ) ( ) *<sup>k</sup>* 

*<sup>i</sup>* **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 ( )*<sup>i</sup>* **Csi** , ( ) *<sup>k</sup>*

*nsi* **β** .

2

*s*

2

*i i*

**R**

*s*

0 ( 3) ( 3)

**s**

*i i*

**R**

are the eigenvectors corresponding to the principal eigenvalues ( ) *<sup>k</sup>*

eigenvectors corresponding to the principal eigenvalues ( ) *<sup>k</sup>*

( )*<sup>i</sup>* **Csi** ,respectively. From (8) we can note that ( ) ( ) *<sup>k</sup>*

( )

*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 ( ) *<sup>s</sup>* **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* **<sup>α</sup>** *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 ( ) *<sup>i</sup>* **a** = <sup>3</sup> ( ) *<sup>i</sup>* **<sup>a</sup>** <sup>=</sup> <sup>=</sup> <sup>4</sup> ( ) [1, ]*<sup>i</sup> <sup>j</sup> <sup>T</sup> <sup>i</sup> e* **a** and 43 4 ( , , , ) [ ( ), ( ), , ( )] [1, ,1, ] *i i T T TT j j <sup>T</sup> 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 ( ) *<sup>s</sup>* **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

> *<sup>i</sup>* **a** = <sup>3</sup> ( )

*i*

**α α R I**

2

*<sup>i</sup> e*

 **a** 

(7)

and

 

*<sup>i</sup>* **<sup>a</sup>** <sup>=</sup> <sup>=</sup> <sup>4</sup> ( ) [1, ]*<sup>i</sup> <sup>j</sup> <sup>T</sup>*

**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 <sup>T</sup> ii i i i i e e*

 

**α**

where 43 4 4 3 4 ( , , , ) [ ( ), ( ), , ( )] *T T TT* **<sup>α</sup>**

 

> 

> >

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 ( ) *<sup>i</sup>* **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 ( )*<sup>i</sup>* **Csi** , ( ) *<sup>k</sup> csi* **β** ( *k K* 1,2, , ) are the eigenvectors corresponding to the principal eigenvalues ( ) *<sup>k</sup> csi* of ( )*<sup>i</sup>* **Csi** , ( ) *<sup>k</sup> rsi* **β** ( 1,2, , *k K* ) are the eigenvectors corresponding to the principal eigenvalues ( ) *<sup>k</sup> rsi* of ( )*<sup>i</sup>* **Rsi** . <sup>2</sup> *<sup>n</sup>* and ( )*l nsi* **β** ( 1,2, ,18 *l K* ) are the noise eigenvalue and the corresponding eigenvectors of ( )*<sup>i</sup>* **Csi** ,respectively. From (8) we can note that ( ) ( ) *<sup>k</sup> <sup>i</sup> rsi* **ai β** ( 1,2, , *k K* ) are in the signal subspace of ( )*<sup>i</sup>* **Csi** , ( )*<sup>l</sup> nsi* **β** ( 1,2, ,18 *l K* ) are in the noise subspace of ( )*i* **Csi** , and ( ) ( ) *<sup>k</sup> <sup>i</sup> rsi* **ai β** are orthogonal to ( )*<sup>l</sup> 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** ( ) *<sup>m</sup>* and <sup>2</sup> ( ) *<sup>s</sup> 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)]*<sup>T</sup>*

 **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 ( ) ( ) *<sup>k</sup> <sup>i</sup> 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 ( ) ( ) *<sup>k</sup>*

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 <sup>β</sup>** (10)

(() ) (() )

*<sup>i</sup> 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 <sup>i</sup>* <sup>1</sup> *<sup>i</sup>* <sup>1</sup> *i i <sup>i</sup>* <sup>3</sup> *<sup>i</sup>* <sup>3</sup> *<sup>i</sup>*<sup>4</sup> *<sup>i</sup>* <sup>4</sup> *<sup>i</sup>* <sup>6</sup> *<sup>i</sup>* <sup>5</sup> *<sup>i</sup>* <sup>2</sup> *<sup>i</sup>* <sup>4</sup> *<sup>i</sup>* <sup>3</sup> *<sup>i</sup>* <sup>1</sup> *<sup>i</sup>*<sup>2</sup> *<sup>i</sup>* <sup>6</sup> *<sup>i</sup>* <sup>7</sup> *<sup>i</sup>* <sup>5</sup> *<sup>i</sup>* <sup>8</sup> *<sup>i</sup>* <sup>9</sup>

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 <sup>T</sup> jjjj jjjj jjjj jjjj <sup>i</sup> 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*

*<sup>T</sup> 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)]*<sup>T</sup>*

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*

*<sup>i</sup> 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), ( ), ( )]*<sup>T</sup>*

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)]*<sup>T</sup>*

 **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 <sup>i</sup>* <sup>1</sup> *<sup>i</sup>* <sup>1</sup> *i i <sup>i</sup>* <sup>3</sup> *<sup>i</sup>* <sup>3</sup> *<sup>i</sup>*<sup>4</sup> *<sup>i</sup>* <sup>4</sup> *<sup>i</sup>*<sup>2</sup> *<sup>i</sup>*<sup>2</sup> *<sup>i</sup>* <sup>1</sup> *<sup>i</sup>*<sup>5</sup> *<sup>i</sup>*<sup>6</sup> *<sup>i</sup>*<sup>7</sup> *<sup>i</sup>* <sup>8</sup> *<sup>i</sup>* <sup>9</sup>

( ) { '( ) '( ) }

*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)*<sup>b</sup> <sup>i</sup>* ( 1)*<sup>a</sup> <sup>i</sup>* ( 5)*<sup>b</sup> i*

( 10)*<sup>a</sup> i*

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

( 2)*<sup>b</sup> i*

( 2)*<sup>a</sup> i* ( 3)*<sup>a</sup> i* ( 1)*<sup>b</sup>* ( 2) *i <sup>b</sup> i*

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

*ai* ( 1)*<sup>a</sup>* ( 2) *<sup>i</sup> <sup>a</sup> <sup>i</sup>* ( 1)*<sup>a</sup> <sup>i</sup>*

( 3)*<sup>b</sup>* ( 6) *i <sup>b</sup>* ( 5) *i <sup>b</sup>* ( 4) *i <sup>b</sup> <sup>i</sup>*

( 2)*<sup>a</sup> i* ( 3)*<sup>a</sup> i* ( 4)*<sup>a</sup> i* ( 5)*<sup>a</sup> i* ( 6)*<sup>a</sup> <sup>i</sup>* ( 7)*<sup>a</sup> i* ( 8)*<sup>a</sup> i* ( 9)*<sup>a</sup> i*

( 6)*<sup>b</sup> i* ( 7)*<sup>b</sup> i* ( 8)*<sup>b</sup> i* ( 9) *<sup>b</sup> <sup>i</sup>*

( 1)*<sup>b</sup> <sup>i</sup> bi* ( 4)*<sup>b</sup> i* ( 5)*<sup>b</sup>* ( 2) *i <sup>b</sup>* ( 3) *i <sup>b</sup> i*

( 11)*<sup>a</sup> i* ( 12)*<sup>a</sup> i* ( 13)*<sup>a</sup> i*

( 1)*<sup>b</sup> i*

**si**( )*i <sup>i</sup>* <sup>1</sup> ( 2) *i <sup>i</sup>* <sup>3</sup> *<sup>i</sup>*<sup>4</sup> *<sup>a</sup> <sup>i</sup>*

*ai* ( 1)*<sup>a</sup> <sup>i</sup>* ( 4)*<sup>a</sup> <sup>i</sup>* ( 5)*<sup>a</sup> <sup>i</sup> bi* ( 1)*<sup>b</sup>* ( 4) *<sup>i</sup> <sup>b</sup>* ( 3) *<sup>i</sup> <sup>b</sup> <sup>i</sup>*

*i* 1 *i* 2 *i* 1 *i* 2 *<sup>i</sup>* <sup>3</sup> *i i i* 1 *i* 4 *i*5 *i* 3 *i*6 *i*7 *A <sup>B</sup> i*4 *i*8 *i*9 *<sup>C</sup> <sup>D</sup>*

( ) *b*

pixels.

( 10)*<sup>a</sup> i*

( 2)*<sup>b</sup> i*

( ) *a*

*SAR* 2

( 10)*<sup>a</sup> i*

( 2)*<sup>b</sup> i*

*ai* ( 1)*<sup>a</sup>* ( 2) *<sup>i</sup> <sup>a</sup> <sup>i</sup>* ( 1)*<sup>a</sup> <sup>i</sup>*

( 1)*<sup>b</sup> <sup>i</sup> bi* ( 4)*<sup>b</sup> i* ( 5)*<sup>b</sup>* ( 2) *i <sup>b</sup>* ( 3) *i <sup>b</sup> i*

( 1)*<sup>b</sup> i*

( 10)*<sup>a</sup>* ( 6) *i <sup>a</sup>* ( 7) *i <sup>a</sup> i* ( 11)*<sup>a</sup> <sup>i</sup>* ( 3)*<sup>b</sup>* ( 2) *i <sup>b</sup> i* ( 6)*<sup>b</sup> i* ( 7)*<sup>b</sup> i*

*ai* ( 1)*<sup>a</sup>* ( 2) *<sup>i</sup> <sup>a</sup> <sup>i</sup>* ( 1)*<sup>a</sup> <sup>i</sup>*

( 3)*<sup>b</sup>* ( 6) *i <sup>b</sup>* ( 5) *i <sup>b</sup>* ( 4) *i <sup>b</sup> <sup>i</sup>*

( 2)*<sup>a</sup> i* ( 3)*<sup>a</sup> i* ( 4)*<sup>a</sup> i* ( 5)*<sup>a</sup> i* ( 6)*<sup>a</sup> <sup>i</sup>* ( 7)*<sup>a</sup> i* ( 8)*<sup>a</sup> i* ( 9)*<sup>a</sup> i*

( 6)*<sup>b</sup> i* ( 7)*<sup>b</sup> i* ( 8)*<sup>b</sup> i* ( 9) *<sup>b</sup> <sup>i</sup>*

( 1)*<sup>b</sup> <sup>i</sup> bi* ( 4)*<sup>b</sup> i* ( 5)*<sup>b</sup>* ( 2) *i <sup>b</sup>* ( 3) *i <sup>b</sup> i*

( 11)*<sup>a</sup> i* ( 12)*<sup>a</sup> i* ( 13)*<sup>a</sup> i*

( 1)*<sup>b</sup> i*

( 2)*<sup>a</sup> i* ( 3)*<sup>a</sup> i* ( 4)*<sup>a</sup> i* ( 5)*<sup>a</sup> i* ( 6)*<sup>a</sup> <sup>i</sup>* ( 7)*<sup>a</sup> i* ( 8)*<sup>a</sup> i* ( 9)*<sup>a</sup> i*

( 6)*<sup>b</sup> i* ( 7)*<sup>b</sup> i* ( 8)*<sup>b</sup> i* ( 9) *<sup>b</sup> <sup>i</sup>*

( 3)*<sup>b</sup>* ( 5) ( 4) *i <sup>b</sup> <sup>i</sup> <sup>b</sup>* ( 6) *i <sup>b</sup> i*

( 11)*<sup>a</sup> i* ( 12)*<sup>a</sup> i* ( 13)*<sup>a</sup> i*

*A B C D*

*A*

*SAR* 2

( 8)*<sup>a</sup> i* ( 9)*<sup>a</sup> i* ( 12)*<sup>a</sup> <sup>i</sup>* ( 13)*<sup>a</sup> i* ( 4)*<sup>b</sup> i* ( 5)*<sup>b</sup> i* ( 8)*<sup>b</sup> i* ( 9)*<sup>b</sup> <sup>i</sup>*

*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 <sup>i</sup>* <sup>1</sup> *<sup>i</sup>* <sup>1</sup> *i i <sup>i</sup>*<sup>3</sup> *<sup>i</sup>* <sup>3</sup> *<sup>i</sup>*<sup>4</sup> *<sup>i</sup>* <sup>4</sup> *<sup>i</sup>*<sup>2</sup> *<sup>i</sup>*<sup>2</sup> *<sup>i</sup>* <sup>1</sup> *<sup>i</sup>*<sup>5</sup> *<sup>i</sup>*<sup>6</sup> *<sup>i</sup>*<sup>7</sup> *<sup>i</sup>* <sup>8</sup> *<sup>i</sup>* <sup>9</sup>


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* **<sup>β</sup>** ( *m K* 1,2, , )of <sup>ˆ</sup> ( )*<sup>i</sup>* **Csi** can be spanned by the Hadamard product vectors ˆ( ) ( ) *<sup>m</sup>*

*ii i rsi rsi rsi* **S ai <sup>c</sup>** *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

*<sup>i</sup> eeee eeee eeee eeee*

ˆˆ ˆ (1) (2) ( ) () ,() , ,() *<sup>K</sup>*

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

> 

. To reduce the computational burden, a

*<sup>i</sup>* is directly obtained by using the fast algorithm.

*<sup>i</sup>* . And

> *i* ,

**ai β ββ ai <sup>β</sup>** (29)

( ) () () ( )

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

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

*<sup>T</sup> jjjj jjjj jjjj jjjj*

**ai** (30)

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

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

20

*K K*

1 1

*<sup>i</sup>* 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*

ˆ( ) *<sup>m</sup>*

where

i.e., ˆ *i i* .

search of

*<sup>i</sup> 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 <sup>ˆ</sup> ( )*<sup>i</sup>* **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 <sup>ˆ</sup> ( )*<sup>i</sup>* **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 <sup>ˆ</sup> ( )*<sup>i</sup>* **Csi** , ˆˆ ˆ ˆ ˆ (1) (2) ( ) ( 1) *K K* (20) *csi csi csi csi csi* , eigenvectors ˆ( )*<sup>l</sup> 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 ˆ( ) *<sup>m</sup> csi* **β** ( 1,2, , *m K* ) corresponding to the principal eigenvalues ˆ( ) *<sup>m</sup> 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 <sup>ˆ</sup> ( )*<sup>i</sup>* **Rsi** can be approximated as the amplitude (i.e., the absolute value) of the estimated covariance matrix <sup>ˆ</sup> ( )*<sup>i</sup>* **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 <sup>ˆ</sup> ( )*<sup>i</sup>* **Rsi** , we obtain *K* principal eigenvectors ˆ( ) *<sup>m</sup> rsi* **β** ( *m K* 1,2, , ). As shown by (8), the same signal subspace spanned by the principal eigenvectors ˆ( ) *<sup>m</sup> csi* **<sup>β</sup>** ( *m K* 1,2, , )of <sup>ˆ</sup> ( )*<sup>i</sup>* **Csi** can be spanned by the Hadamard product vectors ˆ( ) ( ) *<sup>m</sup> <sup>i</sup> 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

<sup>1</sup> <sup>ˆ</sup> ( ) ( )( ) 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 <sup>ˆ</sup> ( )*<sup>i</sup>* **Csi**

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

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

 , eigenvectors ˆ( )*<sup>l</sup>*

> *csi*

*m l*

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

ˆˆ ˆ (1) (2) (20 ) , ,, *<sup>K</sup>*

ˆˆ ˆ (1) (2) ( ) , ,, *<sup>K</sup>*

20 2 ( ) 1

*K l K*

 

<sup>1</sup> <sup>ˆ</sup> <sup>ˆ</sup> <sup>20</sup>

*<sup>n</sup> csi <sup>K</sup> <sup>l</sup>*

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

<sup>ˆ</sup> <sup>ˆ</sup> <sup>2</sup> () () <sup>ˆ</sup>*<sup>n</sup>* **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 <sup>ˆ</sup> ( )*<sup>i</sup>* **Csi** [20], i.e.,

*i*

 

<sup>ˆ</sup> ( )*<sup>i</sup>* **Csi** , ˆˆ ˆ ˆ ˆ (1) (2) ( ) ( 1) *K K* (20) *csi csi csi csi csi*

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

 

whereas the larger eigenvectors ˆ( ) *<sup>m</sup>*

The noise power is often estimated by

*csi* 

estimated by its sample covariance matrix <sup>ˆ</sup> ( )*<sup>i</sup>* **Csi** , i.e.,

size of the sample window.

eigenvalues ˆ( ) *<sup>m</sup>*

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 *<sup>i</sup>* . 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 *<sup>i</sup>* 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 *<sup>i</sup>* 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*

*<sup>i</sup> 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*

*<sup>T</sup> jjjj*

12 12 12 12

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

<sup>22</sup> <sup>22</sup> <sup>22</sup> ( 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*

**<sup>ρ</sup>** (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

 <sup>2</sup> 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** <sup>2</sup>

**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** <sup>2</sup> 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** <sup>2</sup> **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 <sup>i</sup>* <sup>4</sup> *<sup>i</sup>* <sup>3</sup> *i <sup>i</sup>* <sup>1</sup>

12 2 22 ( ) [ ( ), ( 4), ( 3), ( ), ( 1)]*<sup>T</sup>* **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) *<sup>i</sup>*

[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 <sup>e</sup>*

*i i*

*j j*

*i i*

*j j*

*i i*

 

*<sup>i</sup>* 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 ( ) *<sup>n</sup> n*

 and

(B.7)

*angle b*

1

*j*

*i*

*i i*

 

*j j*

so we can get

Since and *<sup>i</sup>* , 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 *<sup>i</sup>*

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*

**<sup>β</sup> <sup>A</sup> <sup>β</sup>**

() ()

*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* **<sup>A</sup> ββ β β** . It can be easily proved that **<sup>A</sup>** ( 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 <sup>e</sup>* **<sup>β</sup> <sup>A</sup> <sup>β</sup> β 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 *<sup>i</sup> k* ( *k* is an integer). Since and *<sup>i</sup>* , 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 *<sup>i</sup>* 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 () () *<sup>H</sup> i ni n* 

 **<sup>A</sup> <sup>β</sup>**

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*

*<sup>i</sup> eeee*

*<sup>T</sup> 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 <sup>β</sup> ββββ A A <sup>β</sup>**

*H HH ll m m*

**ai ββ β β ai**

1 2 3 4

() () ( ) ( )

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

**ai β ββ ai β**

*HH H ll m m*

 

**<sup>A</sup> ββ β β** . It can be easily proved that **<sup>A</sup>** ( 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*

( )

**β**

**<sup>β</sup>** (B.4)

**B A** (B.5)

1 1 *n n b b* ( 2,3,4,5 *<sup>n</sup>* ) (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 ( ) <sup>2</sup>

( ) <sup>2</sup>

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

( ) <sup>2</sup>

*λ*   

*<sup>d</sup> <sup>N</sup>* (2-4)

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

<sup>2</sup> <sup>2</sup> 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 *<sup>i</sup>* , 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

( ) <sup>2</sup> 

*y Txy* (, )

*r*

1 *d* 2 *x*

<sup>1</sup> Flight direction of airborne platform

*i*

*i i r N* (2-1)

*<sup>i</sup>* , 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

*<sup>f</sup>* (2-9)

<sup>0</sup> <sup>100</sup>

, baseline

change

value. Firstly, we preset following parameter: radial distance <sup>1</sup>*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 *<sup>d</sup>*<sup>1</sup> *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*<sup>1</sup> *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*<sup>1</sup> ; *iii* <sup>1</sup> .

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

<sup>1</sup> ( ) <sup>2</sup> 

*d <sup>v</sup> f N <sup>d</sup>*

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

*<sup>v</sup> fN N <sup>d</sup>*

*<sup>v</sup> N N*

( )( ) 2 2

*v r d*

*N d*

0 1 1 cos sin(90 )

*r d*

1 2

 

> 

> >

1 2 12

<sup>1</sup> ( )( ) <sup>2</sup>

( ) <sup>2</sup>

(2-5)

the angle of

(2-6)

(2-7)

(2-8)

*f v <sup>d</sup>* cos

where: *d*<sup>1</sup> *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*<sup>1</sup> ;

Doppler shift based on phase difference measuring

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

*iii* <sup>1</sup> .

*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 <sup>1</sup>*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 *<sup>d</sup>*<sup>1</sup> *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: , *<sup>v</sup>* and *<sup>n</sup>* 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 ) , *<sup>v</sup>* 1 / *m s* , 1 *<sup>n</sup>* .

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 <sup>0</sup> 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*

*<sup>n</sup>* is respectively the mean-root-square error measuring phase

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

2

1 *<sup>n</sup>* .

error of phase difference is just a constant term. If <sup>0</sup> 180

*<sup>v</sup>* 1 / *m s* ,

10-8

**2.2.4 Precision analysis** 

difference is

where:

 

 , *<sup>v</sup>* 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** 

σ<sup>r</sup> ( )/ %

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 <sup>0</sup> 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 <sup>d</sup> <sup>d</sup> 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 <sup>r</sup> di*

1.5% after the angle of arrival is bigger than <sup>0</sup> 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

 *<sup>d</sup> <sup>v</sup> <sup>f</sup>*

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 <sup>1</sup> ( )( ) <sup>2</sup>

*<sup>v</sup> f NN <sup>d</sup>* (2-19)

*tt t*

*ff f*

2

*d*

*ff f*

sin 

1

*r*

*<sup>d</sup>*

*v v <sup>f</sup> r d* (2-17)

*d d <sup>v</sup> <sup>f</sup> <sup>f</sup> <sup>d</sup>* (2-18)

2 2

is tangential velocity.

In where: 1 *v v <sup>t</sup>* 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}
$$

$$
\cos \beta\_1 \approx \frac{\Delta r}{d} \tag{2-14}
$$

By applying differential transformation, we can make a deformation to the Doppler shift equation

$$\begin{aligned} \lambda f\_d &= v \cos \beta\_1 = \frac{v}{\alpha} \frac{\partial \sin \beta\_1}{\partial t} \\\\ &= -\frac{v \Delta r}{\alpha d \sqrt{d^2 - \Delta r^2}} \frac{\partial \Delta r}{\partial t} \\\\ &= \frac{v \Delta r \lambda \Delta f\_d}{\alpha d \sqrt{d^2 - \Delta r^2}} \end{aligned} \tag{2-15}$$

where: *dd d* 2 1 *ff f* is the frequency difference of Doppler shift; / *t* angular velocity.

Making use of circular function and by rearranging, we get

$$
\Delta \mathfrak{o} f\_d = \frac{v}{d} \Delta f\_d \mathfrak{c} \mathfrak{t} \mathfrak{g} \mathcal{B}\_1 \tag{2-16}
$$

Further, the Eq.(2-16) can be deformed by use of angular velocity *<sup>t</sup> <sup>v</sup> r* and expression *f v <sup>d</sup>* cos

$$\frac{v^2 \sin^2 \beta\_1}{\lambda r\_1} = \frac{v}{d} \Delta f\_d \tag{2-17}$$

In where: 1 *v v <sup>t</sup>* sin is tangential velocity.

In fact, when the flight device is uniform motion, the left-hand component of equation is namely basic expression of Doppler rate of change

$$\stackrel{\bullet}{f}\_d = \frac{v^2 \sin^2 \beta}{\lambda r\_1}$$

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

$$\stackrel{\bullet}{f}\_d = \frac{\upsilon}{d} \Delta f\_d \tag{2-18}$$

According to the mathematical definition of Doppler rate

$$\stackrel{\bullet}{f}\_d = \Delta f\_d \ne \Delta t$$

There results:

142 Recent Interferometry Applications in Topography and Astronomy

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

1. The Doppler rate of change can be directly obtained when the wavelength of measured

2. It is wholly not associated with light speed. And this is completely helpful for realizing high-accuracy measuring. The existing analyzing indicates that the result of error

According to the geometric relationship as shown in fig.2-1, we have as following the

<sup>1</sup> sin 

<sup>1</sup> cos

*<sup>r</sup>*

By applying differential transformation, we can make a deformation to the Doppler shift

1

*<sup>v</sup> f v <sup>t</sup>*

 

*vr f dd r*

cos 

*d*

*d d*

where: *dd d* 2 1 *ff f* is the frequency difference of Doppler shift;

Further, the Eq.(2-16) can be deformed by use of angular velocity

Making use of circular function and by rearranging, we get

2 2

*vr r t dd r*

> 1

2 2

*d*

*d r*

2 2

1

sin

*<sup>d</sup>* (2-13)

*<sup>d</sup>* (2-14)

 

*<sup>v</sup> <sup>f</sup> f ctg <sup>d</sup>* (2-16)

 *<sup>t</sup> <sup>v</sup> r*

(2-15)

/ *t* angular

and expression

analysis will become very bad if the resolving result is relative to light speed.

succinct on expression. More important characteristic is lying in:

signal is yet unknown.

approximate representation of circular function

**2.3.2 Derivation** 

equation

velocity.

*f v <sup>d</sup>* cos 
$$\frac{v}{d} = \frac{1}{\Delta 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 is wholly not associated with light speed.

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

$$\begin{aligned} \Delta f\_d &= f\_{d1} - f\_{d2} \\ &= f\_{t1} - f\_{t2} = \Delta f\_t \end{aligned}$$

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

So Doppler changing rate can be determined by the actual measurement to frequency difference of signal.

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 based on phase difference

$$\stackrel{\bullet}{f}\_d = \left(\frac{v}{d}\right)^2 \left[ (\Delta \text{N}\_2 - \Delta \text{N}\_1) + \frac{1}{2\pi} (\Delta \phi\_2 - \Delta \phi\_1) \right] \tag{2-19}$$

The mathematical simulation shown yet that relative calculation error are independent of flight speed and wavelength. The relative calculation error curve is smoother as compared

> r = 50 km r = 100 km r = 300 km

∆t = 5 s

0 10 20 30 40 50 60 70 80 90

Angle of advance β / (°)

Fig. 2-8. Relative calculation error curve of Doppler changing rate against advancing angle

According to mathematical definition, we can realize airborne single-station ranging only based on Doppler shift measurement by approximatively dealing with Doppler changing rate equation. On this basis, the expressing problem of average for Doppler changing rate in a time interval is researched. On condition that airborne measuring station is uniform flight along linear motion, the analysis shows that the average value of Doppler changing rate is directly related to the product of tangential velocity at two terminals and is inversely relative to radial distance at end position in a time interval. And the analog calculation verifies that the better ranging result can be obtained only from this average expression. As contrasted to existing method, the ranging method derived in this text requires neither to

detect Doppler changing rate directly nor to use other localization methods.

0

in different radial distance

**3.1.1 Recapitulation** 

**3. Airborne range finding method** 

**3.1 Principle of Doppler direct range finding** 

0.2

0.4

0.6

0.8

1

1.2

( )

Computational error ε / %

1.4

1.6

to the curve given by existing literature. And the behavior is basically identical.

### **2.3.3 Analog calculation**

In order to verify the accuracy of analysis formula of Doppler rate based on measuring frequency difference, the mathematical simulation is applied with measured value replaced by theoretical value. The wavelength , radial distance <sup>1</sup>*r* , the flight speed of flight device *v* and time interval *t* are set beforehand and making the angle of advance 1 change in prescribed domain. Hence, the rest radial distance and angle of advance can be obtained according to geometric relationship as shown in fig.2-1. Thus, the theoretical value of Doppler shift and Doppler rate of change can be calculated corresponding to every radial distance.

On this basis, Doppler changing rate is computed by Eq.(2-18) and then the relative calculation error is obtained by compared with the theoretical value.

Without the notice, the adopted parameter is as follows: 0.25 *m* , *r km* 100 , *v ms* 100 / , *t s* 5 .

Fig. 2-7. Error curve of Doppler rate against advancing angle in different detecting times.

Fig. 2-7 depicts the relative calculation error curve against the angle of advance for Doppler changing rate when the detecting time interval is different. It can be seen that the error will be conspicuously augmentation if the time interval is too long. Fig.2-8 depicts the relative calculation error is inversely relative to radial distance.

The mathematical simulation shown yet that relative calculation error are independent of flight speed and wavelength. The relative calculation error curve is smoother as compared to the curve given by existing literature. And the behavior is basically identical.

Fig. 2-8. Relative calculation error curve of Doppler changing rate against advancing angle in different radial distance

### **3. Airborne range finding method**

### **3.1 Principle of Doppler direct range finding**

### **3.1.1 Recapitulation**

144 Recent Interferometry Applications in Topography and Astronomy

In order to verify the accuracy of analysis formula of Doppler rate based on measuring frequency difference, the mathematical simulation is applied with measured value replaced

prescribed domain. Hence, the rest radial distance and angle of advance can be obtained according to geometric relationship as shown in fig.2-1. Thus, the theoretical value of Doppler shift and Doppler rate of change can be calculated corresponding to every radial

On this basis, Doppler changing rate is computed by Eq.(2-18) and then the relative

0 10 20 30 40 50 60 70 80 90

Angle of advance β / (°)

Fig. 2-7. Error curve of Doppler rate against advancing angle in different detecting times.

Fig. 2-7 depicts the relative calculation error curve against the angle of advance for Doppler changing rate when the detecting time interval is different. It can be seen that the error will be conspicuously augmentation if the time interval is too long. Fig.2-8 depicts the relative

, radial distance <sup>1</sup>*r* , the flight speed of flight device *v*

∆t = 5 s ∆t = 3 s ∆t = 0.5 s 0.25 *m* , *r km* 100 ,

1 change in

calculation error is obtained by compared with the theoretical value.

Without the notice, the adopted parameter is as follows:

and time interval *t* are set beforehand and making the angle of advance

**2.3.3 Analog calculation** 

*v ms* 100 / , *t s* 5 .

0

calculation error is inversely relative to radial distance.

0.1

0.2

0.3

0.4

0.5

( )

Computational error ε / %

0.6

0.7

0.8

distance.

by theoretical value. The wavelength

According to mathematical definition, we can realize airborne single-station ranging only based on Doppler shift measurement by approximatively dealing with Doppler changing rate equation. On this basis, the expressing problem of average for Doppler changing rate in a time interval is researched. On condition that airborne measuring station is uniform flight along linear motion, the analysis shows that the average value of Doppler changing rate is directly related to the product of tangential velocity at two terminals and is inversely relative to radial distance at end position in a time interval. And the analog calculation verifies that the better ranging result can be obtained only from this average expression. As contrasted to existing method, the ranging method derived in this text requires neither to detect Doppler changing rate directly nor to use other localization methods.

2 22 <sup>2</sup> ( ) 

 *t d d d vt v f t <sup>r</sup> f f*

In fact, Doppler changing rate equation only describes the change of Doppler shift in a certain moment or at a certain point. It can not shows the average change of Doppler shift within a time interval. Later-day research has proved that the speed value derived according to basic physical definition is adverse to one from plane geometry. For speed change within a time interval, the speed value obtaining by plane geometry relationship is

> *f v <sup>d</sup>* (cos cos )

1 1

1

*t*

*v*

 1 2 

The result shows that the average of Doppler changing rate within a certain time is directly related to the product of tangential velocity at two-terminal of flight distance,that is

*<sup>d</sup> t t f v v*

*v*

*d*

cos ( ) cos

1

*<sup>f</sup> <sup>f</sup> <sup>t</sup>*

*d*

sin

   

2 1

111

1

sin

*t*

is associated with the rotary radius vector *r* , that is

(cos cos sin sin cos )

   

and sin

(3-4)

 

*<sup>t</sup> r v* (3-7)

*t r* (3-8)

 2 1 

 as

(3-5)

(3-6)

, we have after

(3-3)

**3.1.3 Improvement of measuring accuracy** 

The Doppler frequency difference can be formulated as

shown in fig.3-1, using approximate expression cos 1

*df v*

Substituting (3-5) into (3-2) it follows that

*v*

*v*

According to the relationship between interior angle and exterior angle

 

more exact.

rearranging for (3-4)

The angular velocity

As a result

### **3.1.2 Basic range finding formula**

Provided that airborne platform is uniform motion along straight line as shown in fig.3-1, the Doppler changing rate equation between detection platform and measured target is

$$\stackrel{\bullet}{f}\_d = \frac{\upsilon\_t^2}{\varkappa \cdot r} \tag{3-1}$$

From the viewpoint of basic mathematical definition, during a time interval *t* , the Doppler changing rate can be approximatively expressed by measured value of Doppler frequency difference between two detection nodes

$$\begin{aligned} \stackrel{\bullet}{f}\_d &= \frac{\Delta f\_d}{\Delta t} \\\\ &= \frac{f\_{d2} - f\_{d1}}{\Delta t} \end{aligned} \tag{3-2}$$

where: *t* is a time interval.

Fig. 3-1. Geometric relationship used for mobile ranging system.

Synthetically applying above-named two equations and speed vector equation <sup>222</sup> *r t vvv* as well as the relationship between radial speed and Doppler shift *v f r d* , the basic range finding formula can be obtained

$$r = \frac{\nu\_t^2 \Delta t}{\lambda \Delta f\_d} = \frac{(\nu^2 - \lambda^2 f\_d^2) \Delta t}{\lambda \left| \Delta f\_d \right|} \tag{3-3}$$

### **3.1.3 Improvement of measuring accuracy**

146 Recent Interferometry Applications in Topography and Astronomy

Provided that airborne platform is uniform motion along straight line as shown in fig.3-1, the Doppler changing rate equation between detection platform and measured target is

> 2 *<sup>t</sup> <sup>d</sup> <sup>v</sup> <sup>f</sup>*

From the viewpoint of basic mathematical definition, during a time interval *t* , the Doppler changing rate can be approximatively expressed by measured value of Doppler frequency

*d*

*<sup>f</sup> <sup>f</sup> <sup>t</sup>*

*d*

*r*

2 1

*d d*

*f f t*

Flight direction of airborne platform

<sup>2</sup>

*S*

as well as the relationship between radial speed and Doppler shift *v f r d*

Synthetically applying above-named two equations and speed vector equation <sup>222</sup> *r t vvv*

, the basic range

Fig. 3-1. Geometric relationship used for mobile ranging system.

(3-1)

(3-2)

**3.1.2 Basic range finding formula** 

difference between two detection nodes

where: *t* is a time interval.

1

finding formula can be obtained

2

1

*r*

In fact, Doppler changing rate equation only describes the change of Doppler shift in a certain moment or at a certain point. It can not shows the average change of Doppler shift within a time interval. Later-day research has proved that the speed value derived according to basic physical definition is adverse to one from plane geometry. For speed change within a time interval, the speed value obtaining by plane geometry relationship is more exact.

The Doppler frequency difference can be formulated as

$$
\lambda \Delta f\_d = \upsilon (\cos \beta\_2 - \cos \beta\_1) \tag{3-4}
$$

According to the relationship between interior angle and exterior angle 2 1 as shown in fig.3-1, using approximate expression cos 1 and sin , we have after rearranging for (3-4)

$$\begin{aligned} \Delta \Delta f\_d &= \upsilon \left[ \cos \left( \beta\_1 + \Delta \beta \right) - \cos \beta\_1 \right] \\\\ &= \upsilon \left( \cos \beta\_1 \cos \Delta \beta - \sin \beta\_1 \sin \Delta \beta - \cos \beta\_1 \right) \\\\ &\approx -\upsilon \Delta \beta \sin \beta\_1 \end{aligned} \tag{3-5}$$

Substituting (3-5) into (3-2) it follows that

$$
\hat{f}\_d = \frac{\Delta f\_d}{\Delta t}
$$

$$
\approx -\frac{v\Delta\beta\sin\beta\_1}{\lambda\Delta t}\tag{3-6}
$$

$$
\approx -\frac{v\_{f1}\alpha}{\lambda}
$$

The angular velocity is associated with the rotary radius vector *r* , that is

$$
tau = \upsilon\_t \tag{3.7}$$

As a result

$$\frac{\Delta f\_d}{\Delta t} \approx -\frac{v\_{t1}v\_{t2}}{\Delta r} \tag{3-8}$$

The result shows that the average of Doppler changing rate within a certain time is directly related to the product of tangential velocity at two-terminal of flight distance,that is

10 20 30 40 50 60 70 80 90

10 20 30 40 50 60 70 80 90

/ ° ( )

Angle of advance β<sup>1</sup>

Fig. 3-3. Improvement of calculation accuracy after mean processing.

/ ° ( )

 Mean value improvement Primitive formula

d = 100 m d = 1000 m d = 3000 m

Angle of advance β<sup>1</sup>

Fig. 3-2 Relative calculation error curve against advancing angle for radial distance with

10-3

0

0.5

Computational error ε / (%)

1

1.5

different displacement.

10-2

10-1

Computational error ε / (%)

100

101

directly related to geometric mean value at two-terminal of flight distance. And it is inversely related to radial distance at terminal position. Right now, the ranging finding equation whose accuracy can be improved can be obtained after expressing tangential velocity in term of speed vector equation and Doppler speed equation

$$\begin{split} r &= \frac{v\_{t1} \cdot v\_{t2} \cdot \Delta t}{\lambda (f\_{d1} - f\_{d2})} \\\\ &= \frac{\sqrt{v^2 - \lambda^2 f\_{d1}^2} \sqrt{v^2 - \lambda^2 f\_{d2}^2} \Delta t}{\lambda (f\_{d1} - f\_{d2})} \end{split} \tag{3-9}$$

### **3.1.4 Analog verification**

In order to validate the accuracy of range finding, the analog calculation is done by replacing measured value with theoretical value. Preassigning following parameter: wavelength , radial distance *r* , flight speed *v* , flight time *t* or spacing *d* , and making 1 continuously change within preassigned interval, the radial distance and angle of advance in other nodes can be calculated in turn according to the geometric relationship as shown in fig. 3-1. Hence, the theoretical value of Doppler shift *di f* corresponding to every radial distance can be computed according to Doppler shift equation.

On this basis, the range can be calculated according to (3-9). The relative calculation error can be obtained by compared with the theoretical value.

Without the notice, the adopted parameter is as follows: *r km* 100 , *v ms* 100 / , *t s* 5 , *t dv* (in where: 1000 *d m* ).

Fig.3-2 is the relative calculation error curve again the angle of advance for ranging calculation with different flight time, time interval in graphs has been transformed into mobile distance *d* . It can be seen that the calculation value obtained from (3-9) has more accuracy if the time interval is less.

Analog calculation proves that the relative calculation error is not associated with wavelength and movement speed of airborne platform. Fig.3-3 presents the comparison of error result between (3-3) and (3-9). Apparently, the calculation accuracy after mean processing has biggish improvement.

In processing analog calculation, we must advert that Doppler frequency difference in denominator of ranging formula must take absolute value. The correct result can be also obtained without absolute value if the operation order of Doppler shift in two detection nodes can be determined according to derivation process presented in 3.1.2 section.

### **3.2 Direct ranging method based on angle rate**

### **3.2.1 Recapitulation**

An airborne Doppler direct ranging method with single baseline is presented. Based on the principle that the angle rate can be determined by Doppler frequency difference between

directly related to geometric mean value at two-terminal of flight distance. And it is inversely related to radial distance at terminal position. Right now, the ranging finding equation whose accuracy can be improved can be obtained after expressing tangential

> 2 22 2 22 1 2 1 2

In order to validate the accuracy of range finding, the analog calculation is done by replacing measured value with theoretical value. Preassigning following parameter:

1 continuously change within preassigned interval, the radial distance and angle of advance in other nodes can be calculated in turn according to the geometric relationship as shown in fig. 3-1. Hence, the theoretical value of Doppler shift *di f* corresponding to every

On this basis, the range can be calculated according to (3-9). The relative calculation error

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

Fig.3-2 is the relative calculation error curve again the angle of advance for ranging calculation with different flight time, time interval in graphs has been transformed into mobile distance *d* . It can be seen that the calculation value obtained from (3-9) has more

Analog calculation proves that the relative calculation error is not associated with wavelength and movement speed of airborne platform. Fig.3-3 presents the comparison of error result between (3-3) and (3-9). Apparently, the calculation accuracy after mean

In processing analog calculation, we must advert that Doppler frequency difference in denominator of ranging formula must take absolute value. The correct result can be also obtained without absolute value if the operation order of Doppler shift in two detection

An airborne Doppler direct ranging method with single baseline is presented. Based on the principle that the angle rate can be determined by Doppler frequency difference between

nodes can be determined according to derivation process presented in 3.1.2 section.

( )

*v f v f t f f*

 

, radial distance *r* , flight speed *v* , flight time *t* or spacing *d* , and making

(3-9)

*d d d d*

velocity in term of speed vector equation and Doppler speed equation

radial distance can be computed according to Doppler shift equation.

can be obtained by compared with the theoretical value.

*t s* 5 , *t dv* (in where: 1000 *d m* ).

accuracy if the time interval is less.

processing has biggish improvement.

**3.2.1 Recapitulation** 

**3.2 Direct ranging method based on angle rate** 

*r*

**3.1.4 Analog verification** 

wavelength

1 2 1 2

*vv t*

*f f*

*t t d d*

( )

Fig. 3-2 Relative calculation error curve against advancing angle for radial distance with different displacement.

Fig. 3-3. Improvement of calculation accuracy after mean processing.

1 *d* 2 Flight direction of airborne platform

2 2 *ctg*

1

sin

*d*

*f*

*dv*

*t*

1

2 2 1

The simulation analysis shows that the frequency difference in denominator must introduce

As soon as the Doppler shift expression based on phase difference measurement is substituted into the Doppler ranging equation derived based on angle change rate, according to the geometric relationship as shown fig.2-2, we can obtain the airborne ranging

2

1 2

*N N*

1 1 ( )( ) 2 2

1( ) <sup>2</sup>

1 1 2 1

*d N d*

1 2

*d d*

*dv f v f*

( )

*d d*

1 2

1 2

 

*rd r* (3-15)

*<sup>d</sup> d f* (3-16)

(3-17)

(3-18)

*T*

Fig. 3-4. Geometric relationship used for airborne Doppler direct ranging.

 sin 

1

*r*

*<sup>t</sup> v r* into (3-16), we obtain

2*r*

1*r*

<sup>1</sup>

By rearrangement, the equation can be expressed

By substituting 1 

absolute value sign.

formula only using phase shift measuring

1

*r*

Because of

*r*

two antennas, the radial distance from measured target to detection platform can be derived making use of angular velocity that is obtained by differential deformation after introducing sine angle rate of change into Doppler shift equation. The analog calculation verifies that derived formula is correct. Moreover, it is proofed that the direct ranging method based on angle rate of change is equivalent to the one based on Doppler rate of change. As a comparison, the direct ranging method based on angle rate of change has following characteristics:


### **3.2.2 Primitive formula**

As shown in fig.3-4, 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 fixed or low speed target *T* , the Doppler shift detected by airborne double channel measuring receiver in no.1 antenna element is

$$
\mathcal{A}f\_{d1} = \upsilon \cos \beta\_1 \tag{3-10}
$$

By differential transmutation, the Doppler shift is expressed in terms of the sine rate of angle of advance

$$
\lambda f\_{d1} = \frac{v}{o\nu} \frac{\mathbf{d} \sin \beta\_1}{\mathbf{d}t} \tag{3-11}
$$

where: 1 *<sup>t</sup> v r* .

Regarding the incident signal from target as parallel wave, according to the geometric relationship in diagram, we approximatively have

$$
\sin \beta\_1 \approx \sqrt{d^2 - \Delta r^2} \Big/ d\tag{3-12}
$$

Again

$$\stackrel{\bullet}{r}\_i = \upsilon\_{ri} = \mathcal{A} f\_{di} \tag{3-13}$$

Resulting in

$$\begin{aligned} \lambda f\_{d1} &= \frac{v}{o} \frac{d \sin \beta\_1}{dt} \\\\ &= \frac{v}{o} \frac{\Delta r}{d} \frac{\lambda \Delta f\_d}{\sqrt{d^2 - \Delta r^2}} \end{aligned} \tag{3-14}$$

where: *dd d* 1 2 *ff f* .

Fig. 3-4. Geometric relationship used for airborne Doppler direct ranging.

Because of

150 Recent Interferometry Applications in Topography and Astronomy

two antennas, the radial distance from measured target to detection platform can be derived making use of angular velocity that is obtained by differential deformation after introducing sine angle rate of change into Doppler shift equation. The analog calculation verifies that derived formula is correct. Moreover, it is proofed that the direct ranging method based on angle rate of change is equivalent to the one based on Doppler rate of change. As a comparison, the direct ranging method based on angle rate of change has following

As shown in fig.3-4, 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 fixed or low speed target *T* , the Doppler shift detected by

1 1

dsin <sup>d</sup>

By differential transmutation, the Doppler shift is expressed in terms of the sine rate of angle

Regarding the incident signal from target as parallel wave, according to the geometric

*r v i ri di* 

dsin <sup>d</sup>

*v r f d d r*

2 2

1

2 2

*d*

1

(3-10)

*<sup>v</sup> <sup>f</sup> <sup>t</sup>* (3-11)

*d rd* (3-12)

*f* (3-13)

(3-14)

*f v <sup>d</sup>* cos

3. It needs not assume that the detecting platform must be uniform motion.

airborne double channel measuring receiver in no.1 antenna element is

1

*d*

<sup>1</sup> sin 

1

*d*

*<sup>v</sup> <sup>f</sup> <sup>t</sup>*

relationship in diagram, we approximatively have

characteristics:

of advance

where: 1 *<sup>t</sup> v r* .

Again

Resulting in

where: *dd d* 1 2 *ff f* .

1. Only requiring a measurement 2. Independent of time measurement

**3.2.2 Primitive formula** 

$$\text{ctg}\,\beta = \Delta r \Big/ \sqrt{d^2 - \Delta r^2} \tag{3-15}$$

By rearrangement, the equation can be expressed

$$
\hbar \, \mathrm{od} \, \sin \beta\_1 = \lambda \Delta \xi\_d \,\tag{3-16}
$$

By substituting 1 *<sup>t</sup> v r* into (3-16), we obtain

$$\begin{aligned} r\_1 &= \frac{d\upsilon\_t \sin \beta\_1}{\lambda \Delta f\_d} \\\\ &= \frac{d\left[\upsilon^2 - \left(\lambda f\_{d1}\right)^2\right]}{\lambda \upsilon \left|\Delta f\_d\right|} \end{aligned} \tag{3-17}$$

The simulation analysis shows that the frequency difference in denominator must introduce absolute value sign.

As soon as the Doppler shift expression based on phase difference measurement is substituted into the Doppler ranging equation derived based on angle change rate, according to the geometric relationship as shown fig.2-2, we can obtain the airborne ranging formula only using phase shift measuring

$$r\_1 = \frac{d\_1 \left[ 1 - \frac{\lambda^2}{d\_1^2} (\Delta N\_1 + \frac{\Delta \phi\_1}{2\pi})^2 \right]}{\lambda \left| \frac{1}{d\_1} (\Delta N\_1 + \frac{\Delta \phi\_1}{2\pi}) - \frac{1}{d\_2} (\Delta N\_2 + \frac{\Delta \phi\_2}{2\pi}) \right|} \tag{3-18}$$

0 10 20 30 40 50 60 70 80 90

Angle of arrival θ / (°)

As shown in fig. 3-6, provided that the detection platform which is uniform motion along straight line detects the signal of target by fixed cycle and carries out at least three continuing measurement respectively to earth-fixed target. The expression of Doppler

*r* ( 1,2,3) *<sup>i</sup>* (3-21)

(3-22)

(3-23)

Fig. 3-5. Computational error of ranging formula with different baseline length.

**3.3.2 Ratio of Doppler changing rate between two adjacent detecting nodes** 

 2 

*<sup>v</sup> <sup>f</sup>*

*di*

The ratio of Doppler changing rate between two adjacent detecting nodes is

*q*

1

*ti*

*i*

 *<sup>d</sup> <sup>t</sup> <sup>t</sup> <sup>d</sup>*

sin sin sin sin

 

*i i i ti i i i ti r vv r vv*

Namely, in the case of uniform motion, the ratio of radial distance between two adjacent detecting nodes equals the ratio of tangential velocity. Substituting Eq.(3-23) into (3-22)

*f r v*

*<sup>r</sup> <sup>v</sup> <sup>f</sup>*

The ratio of radial distance between two adjacent detecting nodes according to sine theorem

2 2 1 2 2 <sup>2</sup> <sup>1</sup> <sup>1</sup>

1 1 ( 1)

 

 

0

0.01

changing rate at per node is

is

gives

0.02

0.03

0.04

Computational error ε / (%)

0.05

0.06

0.07

0.08

### **3.2.3 Error between computational and theoretical value**

Based on the equivalence in computation, the ranging formula (3-18) can be transformed into as follows form

$$r = \frac{d\_1 \left[ 1 - \frac{\Delta r\_1^2}{d\_1^2} \right]}{\left| \frac{\Delta r\_1}{d\_1} - \frac{\Delta r\_2}{d\_2} \right|}\tag{3-19}$$

Preassigning following parameter: radial distance *r* and baseline length *d* , and making the angle of arrival <sup>0</sup> 1 90 continuously change within preassigned interval. The other radial distance can be in turn solved making use of cosine law. Then, the path difference is obtained. With that, the radial distance of target can be calculated by (3-19) and the relative calculation error between computational and preassigning value can be obtained by comparison

$$
\omega = \frac{\left| X - X\_a \right|}{X} \times 100\,\%\tag{3-20}
$$

where: *X* and *Xa* respectively express computational and preassigning value.

The benefit using (3-19) in analog calculation is that the analysis for integer of wavelength and phase difference value can be avoided.

Fig.3-5 shows the error curve of radial distance between computational and preassigning value in which the total length of baseline is different when the angle of arrival is linear variation from <sup>0</sup> 0 to <sup>0</sup> 90 . Its basic feature is that the longer is the baseline length, the greater is the computational error, namely, the computational error is inversely proportional to radial distance.

The basic parameters adopting in simulation analysis is: 1*r km* 50 and *dd d* 1 2 .

### **3.3 Direct ranging method based on Doppler rate of change**

### **3.3.1 Recapitulation**

Under the condition that detection platform is uniform motion along straight line and continuing multipoint detection, on the one hand, the ratio of Doppler changing rate between two adjacent detecting nodes can be expressed as the cubic of specific value of tangential velocities making use of circular function relationship. On the other hand, making use of the formula of Doppler rate of change based on measuring Doppler frequency difference, the ratio of Doppler changing rate between two adjacent detecting nodes can be also expressed as the specific value of Doppler frequency difference. On this basis, a direct ranging formula based on Doppler shift as well as frequency difference can be obtained from the identical relation of speed vector.

Based on the equivalence in computation, the ranging formula (3-18) can be transformed

1 *<sup>r</sup> <sup>d</sup> d*

Preassigning following parameter: radial distance *r* and baseline length *d* , and making the

radial distance can be in turn solved making use of cosine law. Then, the path difference is obtained. With that, the radial distance of target can be calculated by (3-19) and the relative

100%

The benefit using (3-19) in analog calculation is that the analysis for integer of wavelength

Fig.3-5 shows the error curve of radial distance between computational and preassigning value in which the total length of baseline is different when the angle of arrival is linear variation from <sup>0</sup> 0 to <sup>0</sup> 90 . Its basic feature is that the longer is the baseline length, the greater is the computational error, namely, the computational error is inversely proportional

Under the condition that detection platform is uniform motion along straight line and continuing multipoint detection, on the one hand, the ratio of Doppler changing rate between two adjacent detecting nodes can be expressed as the cubic of specific value of tangential velocities making use of circular function relationship. On the other hand, making use of the formula of Doppler rate of change based on measuring Doppler frequency difference, the ratio of Doppler changing rate between two adjacent detecting nodes can be also expressed as the specific value of Doppler frequency difference. On this basis, a direct ranging formula based on Doppler shift as well as frequency difference can be obtained

The basic parameters adopting in simulation analysis is: 1*r km* 50 and *dd d* 1 2 .

**3.3 Direct ranging method based on Doppler rate of change** 

*X Xa*

*r r d d*

*r*

where: *X* and *Xa* respectively express computational and preassigning value.

continuously change within preassigned interval. The other

*<sup>X</sup>* (3-20)

between computational and preassigning value can be obtained by

(3-19)

**3.2.3 Error between computational and theoretical value** 

1

into as follows form

angle of arrival <sup>0</sup>

calculation error

to radial distance.

**3.3.1 Recapitulation** 

comparison

90

and phase difference value can be avoided.

from the identical relation of speed vector.

Fig. 3-5. Computational error of ranging formula with different baseline length.

### **3.3.2 Ratio of Doppler changing rate between two adjacent detecting nodes**

As shown in fig. 3-6, provided that the detection platform which is uniform motion along straight line detects the signal of target by fixed cycle and carries out at least three continuing measurement respectively to earth-fixed target. The expression of Doppler changing rate at per node is

$$\stackrel{\bullet}{f}\_{di} = \frac{v\_{ti}^2}{\lambda r\_i} \quad \text{( $i = 1, 2, 3$ )}\tag{3-21}$$

The ratio of Doppler changing rate between two adjacent detecting nodes is

$$q = \frac{\stackrel{\bullet}{f}\_{d2}}{\stackrel{\bullet}{f}\_{d1}} = \frac{r\_1}{r\_2} \frac{v\_{t2}^2}{v\_{t1}^2} \tag{3-22}$$

The ratio of radial distance between two adjacent detecting nodes according to sine theorem is

$$\frac{r\_{i+1}}{r\_i} = \frac{\sin \beta\_i}{\sin \beta\_{i+1}} = \frac{v \sin \beta\_i}{v \sin \beta\_{i+1}} = \frac{v\_{\text{fi}}}{v\_{t(i+1)}} \tag{3-23}$$

Namely, in the case of uniform motion, the ratio of radial distance between two adjacent detecting nodes equals the ratio of tangential velocity. Substituting Eq.(3-23) into (3-22) gives

$$q = \frac{\stackrel{\bullet}{f}\_{d2}}{\stackrel{\bullet}{f}\_{d1}} = \frac{\upsilon\_{t2}^3}{\upsilon\_{t1}^3} \tag{3-24}$$

Then, separately substituting the Doppler shift which is associated with radial velocity and the Doppler changing rate which is associated with tangential velocity as well as their ratio

> 22 1 2 1 22 ( ) ( 1)

> > 2 2 2 1

( )

*d d*

( 1)

0 10 20 30 40 50 60 70 80

/ (°)

Angle of advance β<sup>1</sup>

Fig. 3-7. Computational error of ranging formula in equally spaced detection.

2

*d*

211 2 1

*d ff f u vf*

*dd d d*

2

( )

1

 

*f f <sup>r</sup>*

( 1)

*u f*

*f f rf u dd d* (3-28)

(3-29)

2

Hence, we can obtain the calculating formula of radial distance

2

= 500 m

= 1000 m

= 1500 m

3 2 1 2 3

*<sup>d</sup> <sup>f</sup> u q d f* .

0

0.5

1

1.5

Computational error ε

r ( )/ %

2

2.5

3

2 1

*d d*

d1 = d2

d1 = d2

d1 = d2

into (3-27) gives

where:

Hence, when detection platform is uniform motion, the ratio of Doppler change rate between two adjacent detecting nodes will equal the cubic of specific value of tangential velocities on two adjacent nodes. Introducing the Doppler change rate expression (2-18) based on Doppler frequency difference, the specific value of Doppler change rate between two adjacent detecting nodes can be also written as

$$q = \frac{d\_1}{d\_2} \frac{\Delta f\_{d2}}{\Delta f\_{d1}}\tag{3-25}$$

Fig. 3-6. Geometric model used for analyzing Doppler passive localization of moving single station.

### **3.3.3 Solution of radial distance**

According to the velocity component at every node in mobile path of the platform, we have identical relation about speed

$$\upsilon^2 = \upsilon\_{r2}^2 + \upsilon\_{t2}^2 = \upsilon\_{r1}^2 + \upsilon\_{t1}^2 \tag{3-26}$$

Rearranging yield

$$
v\_{r2}^2 - v\_{r1}^2 = v\_{t1}^2 - v\_{t2}^2\tag{3-27}$$

Then, separately substituting the Doppler shift which is associated with radial velocity and the Doppler changing rate which is associated with tangential velocity as well as their ratio into (3-27) gives

$$
\lambda (f\_{d2}^2 - f\_{d1}^2) = r\_2 \stackrel{\bullet}{f\_{d2}} (\mu^{-1} - 1)\tag{3-28}
$$

 $\mu = \sqrt[3]{q^2} = \sqrt[3]{\left(\frac{d\_1}{d\_2} \frac{\Delta f\_{d2}}{\Delta f\_{d1}}\right)^2}$ .

Hence, we can obtain the calculating formula of radial distance

where:

154 Recent Interferometry Applications in Topography and Astronomy

*f v*

*<sup>v</sup> <sup>f</sup>*

Hence, when detection platform is uniform motion, the ratio of Doppler change rate between two adjacent detecting nodes will equal the cubic of specific value of tangential velocities on two adjacent nodes. Introducing the Doppler change rate expression (2-18) based on Doppler frequency difference, the specific value of Doppler change rate between

> 1 2 2 1 *d d*

> > 3

*r*

*r*

22 2 2 2 *rt rt* 22 11 *vv v v v* (3-26)

2 2 22 *vv vv rr tt* 2 112 (3-27)

Flight direction of airborne platform

2

<sup>1</sup> 3

*r*

*T*

Fig. 3-6. Geometric model used for analyzing Doppler passive localization of moving single

According to the velocity component at every node in mobile path of the platform, we have

*<sup>d</sup> <sup>f</sup> <sup>q</sup> <sup>d</sup> <sup>f</sup>* (3-25)

 *<sup>d</sup> <sup>t</sup> <sup>t</sup> <sup>d</sup>*

*q*

two adjacent detecting nodes can be also written as

3

*r* 2

1

2 *d*

2

1 *d*

1

**3.3.3 Solution of radial distance** 

identical relation about speed

Rearranging yield

station.

(3-24)

$$\begin{aligned} r\_2 &= \frac{\lambda (f\_{d2}^2 - f\_{d1}^2)}{\left(\mu^{-1} - 1\right) f\_{d2}} \\\\ &= \frac{\lambda d\_2 \Delta f\_{d1} (f\_{d1} + f\_{d2})}{\left(\mu^{-1} - 1\right) v \Delta f\_{d2}} \end{aligned} \tag{3-29}$$

Fig. 3-7. Computational error of ranging formula in equally spaced detection.

wavelength is very short for high-frequency signal. This moment, not only the antenna element must be made to do very small, but also very high demand is put forward for antenna arrangement. It will bring about coupled between antennas and bring down antenna gain. At the same time, higher demand will be required for measurement accuracy of interferometer. For algorithm resolving phase ambiguity with multi-baseline, the

The study presented in this section shows that applying a single baseline array in airborne platform can realize the high-accuracy DF without the ambiguity of viewing angle after combining with Doppler information. Only from the viewpoint of principle of measurement, selection for baseline length is arbitrary. So it is more suitable for carrying out

The analysis discovers that the integer of wavelength in radial distance can be directly obtained compositely making use of the velocity vector equation and Doppler shift as well as Doppler changing rate equation. From this, the integer of wavelength in path length difference of radial distance between two adjacent antenna elements can be determined. Further, the value less than a wavelength in path length difference can be obtained by phase difference measurement. As compared with now existing interferometry firstly determining phase difference, this sort of direction finding method associating with Doppler and phase difference and firstly determining path length difference does not exist phase ambiguity nor

According to derivation in section three, in fact, we can obtain two relations from the

1 2 11 ( ) ( 1)

1 2 22 <sup>1</sup> ( ) (1 )

int int

*<sup>r</sup> f f <sup>N</sup>*

*f f rf dd d <sup>u</sup>*

2 2

*d d*

2 2

*d d*

( 1)

( 1)

1 1 2

2 1 2

( ) int int

*<sup>r</sup> <sup>f</sup> f u <sup>N</sup>*

Substituting these expressions into the existing DF expression based on the principle of phase interferometry, we can obtain the airborne DF formula based on phase shift and

1

2

*f u*

 

*d*

*f u*

 

*d*

*f f rf u dd d* (4-1)

(4-2)

(4-3)

(4-4)

2 2

2 2

Thus, we can obtain integer values of wavelength about two radial distances

1

2

computing amount is heavy due to demanding multidimensional integer search.

detection operation in a broadband.

require restricting baseline length.

**4.2 Analytical derivation** 

identical relation about speed

Doppler shift as well as its rate

Fig. 3-8. Computational error of ranging formula in unequal spaced detection.

According to calculation error formula between calculated value and theoretical value, fig.3- 7 shows the computational error curve between computational and theoretical value when the platform is movement in equal interval. Fig.3-8 is not equal interval. It can be seen that the formula has best accuracy when two adjacent flight distances is equal. If the distance is not equal, the error is larger. There is divergence phenomenon when the angle of advance goes to <sup>0</sup> 90 . The simulation calculation shows that the error is not connected with the flight speed of mobile platform and wavelength of measured signal.

The basic parameters adopting in simulation analysis is: *v ms* 100 / , <sup>1</sup>*r km* 100 and 0.25 *m* .

### **4. Airborne DF method without ambiguity based on Doppler as well as rate**

### **4.1 Recapitulation**

The phase interferometry is a direction finding method with better measurement accuracy. It is widely used for active and passive detection system. But for single baseline phase interferometry, there is the contradiction between accuracy of direction finding and maximum angle without ambiguity. To solve this problem, existing method is to utilize multi-baseline system including the method combining long baselines with short ones and algorithm resolving phase ambiguity with multi-baseline.

In actual application, the method combining long baselines with short ones have two limitations. In fact, corresponding baselines will also become extremely small since wavelength is very short for high-frequency signal. This moment, not only the antenna element must be made to do very small, but also very high demand is put forward for antenna arrangement. It will bring about coupled between antennas and bring down antenna gain. At the same time, higher demand will be required for measurement accuracy of interferometer. For algorithm resolving phase ambiguity with multi-baseline, the computing amount is heavy due to demanding multidimensional integer search.

The study presented in this section shows that applying a single baseline array in airborne platform can realize the high-accuracy DF without the ambiguity of viewing angle after combining with Doppler information. Only from the viewpoint of principle of measurement, selection for baseline length is arbitrary. So it is more suitable for carrying out detection operation in a broadband.

The analysis discovers that the integer of wavelength in radial distance can be directly obtained compositely making use of the velocity vector equation and Doppler shift as well as Doppler changing rate equation. From this, the integer of wavelength in path length difference of radial distance between two adjacent antenna elements can be determined. Further, the value less than a wavelength in path length difference can be obtained by phase difference measurement. As compared with now existing interferometry firstly determining phase difference, this sort of direction finding method associating with Doppler and phase difference and firstly determining path length difference does not exist phase ambiguity nor require restricting baseline length.

### **4.2 Analytical derivation**

156 Recent Interferometry Applications in Topography and Astronomy

= 1100 m

= 1000 m

0 10 20 30 40 50 60 70 80 90

/ (°)

Angle of advance β<sup>1</sup>

According to calculation error formula between calculated value and theoretical value, fig.3- 7 shows the computational error curve between computational and theoretical value when the platform is movement in equal interval. Fig.3-8 is not equal interval. It can be seen that the formula has best accuracy when two adjacent flight distances is equal. If the distance is not equal, the error is larger. There is divergence phenomenon when the angle of advance goes to <sup>0</sup> 90 . The simulation calculation shows that the error is not connected with the flight

The basic parameters adopting in simulation analysis is: *v ms* 100 / , <sup>1</sup>*r km* 100 and

**4. Airborne DF method without ambiguity based on Doppler as well as rate** 

The phase interferometry is a direction finding method with better measurement accuracy. It is widely used for active and passive detection system. But for single baseline phase interferometry, there is the contradiction between accuracy of direction finding and maximum angle without ambiguity. To solve this problem, existing method is to utilize multi-baseline system including the method combining long baselines with short ones and

In actual application, the method combining long baselines with short ones have two limitations. In fact, corresponding baselines will also become extremely small since

Fig. 3-8. Computational error of ranging formula in unequal spaced detection.

speed of mobile platform and wavelength of measured signal.

algorithm resolving phase ambiguity with multi-baseline.

2.5

3

3.5

4

4.5

Computational error ε

0.25 *m* .

**4.1 Recapitulation** 

r ( )/ %

5

5.5

6

6.5

7

d1

d1

= 1000 m, d2

= 1100 m, d2

According to derivation in section three, in fact, we can obtain two relations from the identical relation about speed

$$
\mathcal{A}(f\_{d1}^2 - f\_{d2}^2) = r\_1 \stackrel{\bullet}{f\_{d1}} (\mu - 1) \tag{4-1}
$$

$$
\mathcal{A}(f\_{d1}^2 - f\_{d2}^2) = r\_2 \stackrel{\bullet}{f}\_{d2} (1 - \frac{1}{u}) \tag{4-2}
$$

Thus, we can obtain integer values of wavelength about two radial distances

$$N\_1 = \text{int}\left[\frac{r\_1}{\mathcal{A}}\right] = \text{int}\left[\frac{f\_{d1}^2 - f\_{d2}^2}{\stackrel{\bullet}{f\_{d1}}(\mu - 1)}\right] \tag{4-3}$$

$$N\_2 = \text{int}\left[\frac{r\_2}{\mathcal{A}}\right] = \text{int}\left[\frac{(f\_{d1}^2 - f\_{d2}^2)\mu}{\stackrel{\bullet}{f}\_{d2}(\mu - 1)}\right] \tag{4-4}$$

Substituting these expressions into the existing DF expression based on the principle of phase interferometry, we can obtain the airborne DF formula based on phase shift and Doppler shift as well as its rate

As shown in fig.5-1, three antenna units arrange in L form in horizontal plane. Two baselines are placed at right angles to each other. In which, the direction of one baseline is parallel to axes of aerocraft. On condition that measured target is motionless or low speed

Approximately upon condition that the incident wave is parallel, according to the direction

1 1

*t td*

1 2 1

*r r d*

*f f <sup>d</sup>*

1

Hence, making use of two antenna units whose baseline is parallel to the axis of aerial vehicle can obtain the sinusoidal triangle function with regard to relative azimuth depending on the

1 cos sin

1

*T*

*r*

2

<sup>1</sup> *S* 1 *d* 2 *S* Flight direction of airborne platform

*r*

<sup>3</sup> 1

2

1

*r*

1 2

1

1 1

*t r t d*

1 1

*d*

*f d*

( )

*d d*

*<sup>i</sup> f v* (5-1)

1 sin

*v r* *r*

(5-2)

(5-3)

angular velocity; <sup>1</sup> *d*

motion, the Doppler shift received by airborne receiver in every antenna unit is

*di* cos

cos

Doppler frequency difference and on the angular velocity and on the base length

1

where: *d dd* <sup>112</sup> *fff* is Doppler frequency difference; <sup>1</sup>

baseline length parallel to the axis of aerial vehicle..

3

Fig. 5-1. Airborne passive direction finding array with L form.

3 *S*

2 *d*

1 *r*

1 cosine change rate, we have

$$\begin{split} \sin \theta &= \frac{\Delta r}{L} \\ &= \frac{\lambda}{L} \Bigg( \text{int} \Bigg[ \frac{f\_{d1}^2 - f\_{d2}^2}{\bullet} \Bigg[ -\text{int} \Bigg[ \frac{(f\_{d1}^2 - f\_{d2}^2) \mu}{\bullet} \Bigg] - \text{int} \Bigg[ \frac{(f\_{d1}^2 - f\_{d2}^2) \mu}{f\_{d2} (\mu - 1)} \Bigg] + \frac{\Delta \phi}{2\pi} \Bigg) \end{split} \tag{4-5}$$

### **5. Airborne passive DF with orthogonal baseline**

### **5.1 Recapitulation**

At present, main methods applicable to airborne DF have the amplitude comparison and phase interferometry, etc. The measuring precision of direction-finding system based on amplitude comparison have always had to suffer the biggish influence conduced by incompatible from antenna and reception channel of measuring receiver. Moreover, phase interferometry needs to solve phase ambiguity.

This section presents a Doppler DF method applicable to airborne based on the direction cosine change rate.

If three antenna units are divided into two set and two baselines are placed at right angles to each other, in which the direction of one baseline is parallel to the actual flight direction of air vehicle, the sine and cosine function of target bearing respectively in two baseline directions can be simultaneously obtained according to the analysis principle of the direction cosine change rate. The angle measurement formula only based on Doppler frequency difference can be derived after eliminating the unknown parameters including angular velocity and wavelength by the specific value of two circular functions. The analog calculation shows that the relative calculation error is in direct proportion to the baseline length provided that the incident wave is parallel in derivation. Furthermore, the derived formula has irregularity in airborne axis direction. But the error analysis depicts that the measurement accuracy is in direct proportion to the baseline length. Moreover, the measurement accuracy when the azimuth angle is minor can be usefully enhanced by changing the specific value between two baselines. Since the new method is not associated with wavelength, the direction finding only based on Doppler frequency difference will be more adapted to passive sounding as compared with phase interference method.

### **5.2 Derivation**

If we can use the direction cosine change rate for single baseline array, the incidence angle of measured signal can be expressed as the function depending on the Doppler frequency difference and on the angular velocity and on the baseline length and on the wavelength. According to this result, the study discovery that the airborne direction finding only based on Doppler frequency difference can be realized by making use of the direction cosine change rate for orthogonal double-baseline. Firstly, a planar array with L-shape is structured by use of three antennas. Then, the sine and cosine circular function with regard to the incidence angle of measured signal can be simultaneously obtained due to the direction cosine change rate. Moreover, the unknown angular velocity and the wavelength can be eliminated by the specific value of the sine and cosine circular function. Thus, the derived tangent angle is only associated with known Doppler frequency difference and the baseline length.

sin

interferometry needs to solve phase ambiguity.

**5.1 Recapitulation** 

cosine change rate.

**5.2 Derivation** 

*r L*

**5. Airborne passive DF with orthogonal baseline** 

22 22 12 12

*dd dd*

 

*f f f fu*

<sup>2</sup> ( 1) ( 1)

(4-5)

1 2

At present, main methods applicable to airborne DF have the amplitude comparison and phase interferometry, etc. The measuring precision of direction-finding system based on amplitude comparison have always had to suffer the biggish influence conduced by incompatible from antenna and reception channel of measuring receiver. Moreover, phase

This section presents a Doppler DF method applicable to airborne based on the direction

If three antenna units are divided into two set and two baselines are placed at right angles to each other, in which the direction of one baseline is parallel to the actual flight direction of air vehicle, the sine and cosine function of target bearing respectively in two baseline directions can be simultaneously obtained according to the analysis principle of the direction cosine change rate. The angle measurement formula only based on Doppler frequency difference can be derived after eliminating the unknown parameters including angular velocity and wavelength by the specific value of two circular functions. The analog calculation shows that the relative calculation error is in direct proportion to the baseline length provided that the incident wave is parallel in derivation. Furthermore, the derived formula has irregularity in airborne axis direction. But the error analysis depicts that the measurement accuracy is in direct proportion to the baseline length. Moreover, the measurement accuracy when the azimuth angle is minor can be usefully enhanced by changing the specific value between two baselines. Since the new method is not associated with wavelength, the direction finding only based on Doppler frequency difference will be

more adapted to passive sounding as compared with phase interference method.

associated with known Doppler frequency difference and the baseline length.

If we can use the direction cosine change rate for single baseline array, the incidence angle of measured signal can be expressed as the function depending on the Doppler frequency difference and on the angular velocity and on the baseline length and on the wavelength. According to this result, the study discovery that the airborne direction finding only based on Doppler frequency difference can be realized by making use of the direction cosine change rate for orthogonal double-baseline. Firstly, a planar array with L-shape is structured by use of three antennas. Then, the sine and cosine circular function with regard to the incidence angle of measured signal can be simultaneously obtained due to the direction cosine change rate. Moreover, the unknown angular velocity and the wavelength can be eliminated by the specific value of the sine and cosine circular function. Thus, the derived tangent angle is only

*d d*

*<sup>L</sup> fu fu*

( ) int int

As shown in fig.5-1, three antenna units arrange in L form in horizontal plane. Two baselines are placed at right angles to each other. In which, the direction of one baseline is parallel to axes of aerocraft. On condition that measured target is motionless or low speed motion, the Doppler shift received by airborne receiver in every antenna unit is

$$\mathcal{A}f\_{di} = \upsilon \cos \theta\_i \tag{5-1}$$

Approximately upon condition that the incident wave is parallel, according to the direction cosine change rate, we have

$$\begin{aligned} \frac{\partial \cos \theta\_1}{\partial t} &= \frac{\partial}{\partial t} \left( \frac{\Delta r\_1}{d\_1} \right) \\ &\qquad \bullet \quad \text{\textquotedblleft} \\ &= \frac{r\_1 - r\_2}{d\_1} \\ &= \frac{\lambda}{d\_1} (f\_{d1} - f\_{d2}) \end{aligned} \tag{5-2}$$

Hence, making use of two antenna units whose baseline is parallel to the axis of aerial vehicle can obtain the sinusoidal triangle function with regard to relative azimuth depending on the Doppler frequency difference and on the angular velocity and on the base length

$$\begin{split} \sin \theta\_1 &= -\frac{1}{o\rho\_\theta} \frac{\partial \cos \theta\_1}{\partial t} \\ &= -\frac{1}{o\rho\_\theta} \frac{\partial}{\partial t} \left( \frac{\Delta r\_1}{d\_1} \right) \\ &= -\frac{\lambda \Delta f\_{d1}}{o\rho\_\theta d\_1} \end{split} \tag{5-3}$$

where: *d dd* <sup>112</sup> *fff* is Doppler frequency difference; <sup>1</sup> 1 sin *v r* angular velocity; <sup>1</sup> *d* baseline length parallel to the axis of aerial vehicle..

Fig. 5-1. Airborne passive direction finding array with L form.

Fig.5-2 and fig.5-3 depict the relative calculation error curve with different baseline and radial distance. When azimuth angle tends to zero, the derived formula has irregularity. Obviously, the farther is the object distance, or the shorter is the baseline length, the smaller is the relative calculation error of formula. There occurs this occurrence provided that the

d1 = d 2

d1 = d 2

d1 = d 2

r 1 = 10 km

r 1 = 50 km

r 1 = 100 km

= 5 m

= 10 m

= 20 m

0 10 20 30 40 50 60 70 80 90

0 10 20 30 40 50 60 70 80 90

/ (°)

Azimuth angle θ<sup>1</sup>

Azimuth angle θ1 / (°)

incident wave is parallel in derivation.

10-5

10-5

10-4

10-3

10-2

Relative error ε / (%)

10-1

100

101

10-4

10-3

Basic parameters:

Fig. 5-2. Relative calculation error: different baseline lengths.

Basic parameters:

d1 = d2 = 5 m

Fig. 5-3. Relative calculation error: different radial distances.

r 1 = 100 km

10-2

Relative error ε / (%)

10-1

100

101

Further, making use of two antenna units whose baseline is at right angles to the axis of aerial vehicle can obtain the cosine triangle function with regard to relative azimuth depending on the Doppler frequency difference and on the angular velocity and on the base length

$$\begin{aligned} \cos \theta\_1 &= \sin(90 - \theta\_1) \\\\ &= -\frac{1}{a\rho\_\theta} \frac{\partial \cos(90 - \theta\_1)}{\partial t} \\\\ &= -\frac{1}{a\rho\_\theta} \frac{\partial}{\partial t} \left(\frac{\Delta r\_2}{d\_2}\right) = -\frac{\lambda \Delta f\_{d2}}{a\rho\_\theta d\_2} \end{aligned} \tag{5-4}$$

where: *d dd* 2 13 *f f f* is Doppler frequency difference; 2 *d* baseline length at right angles to the axis of aerial vehicle.

By way of the specific value, the tangential triangle function only depending on the Doppler frequency difference and on the baseline length can be obtained

$$\begin{aligned} \, &tg \theta\_1 = \frac{\sin \theta\_1}{\cos \theta\_1} \\\\ &= \frac{d\_2}{d\_1} \frac{\Delta f\_{d1}}{\Delta f\_{d2}} \end{aligned} \tag{5-5}$$

The relative azimuth between aerial vehicle and measured target is

$$\theta\_1 = \text{tg}^{-1}\left[\frac{d\_2}{d\_1}\frac{\Delta f\_{d1}}{\Delta f\_{d2}}\right] \tag{5-6}$$

This analytic function in form is analogous to the formula of the amplitude comparison and phase interference method. Moreover, it is also not connected with wavelength.

### **5.3 Simulation analysis**

We make the analog verification by replacing measured value with theoretical value. Firstly, presetting following parameter: radial distance 1*r* and base length *<sup>i</sup> d* as well as wavelength and speed, the theoretical value of the rest radial distance and azimuth angle can be computed by making azimuth angle continuous change in specified interval. Hence, we can obtain the theoretical value of the Doppler shift corresponding to each radial distance. After that, the value of azimuth angle can be calculated by Eq.(5-6) and the relative calculation error can be obtained by comparison with the theoretical value.

Because the simulation analysis of the relative calculation error is not associated with both wavelength and speed, there is not the specified value of wavelength and speed.

Further, making use of two antenna units whose baseline is at right angles to the axis of aerial vehicle can obtain the cosine triangle function with regard to relative azimuth depending on the Doppler frequency difference and on the angular velocity and on the base

1 cos(90 )

*t*

 

where: *d dd* 2 13 *f f f* is Doppler frequency difference; 2 *d* baseline length at right angles to

By way of the specific value, the tangential triangle function only depending on the Doppler

sin cos

1

1

2 1 1 2

1 2 1

*d*

This analytic function in form is analogous to the formula of the amplitude comparison and

We make the analog verification by replacing measured value with theoretical value. Firstly, presetting following parameter: radial distance 1*r* and base length *<sup>i</sup> d* as well as wavelength and speed, the theoretical value of the rest radial distance and azimuth angle can be computed by making azimuth angle continuous change in specified interval. Hence, we can obtain the theoretical value of the Doppler shift corresponding to each radial distance. After that, the value of azimuth angle can be calculated by Eq.(5-6) and the relative calculation

Because the simulation analysis of the relative calculation error is not associated with both

wavelength and speed, there is not the specified value of wavelength and speed.

1 2 

*d*

*<sup>d</sup> <sup>f</sup> tg d f* (5-6)

*d d*

*d f d f* 1

2 2 2 2

*r f td d*

 

  *d*

(5-4)

(5-5)

1 1

cos sin(90 )

frequency difference and on the baseline length can be obtained

The relative azimuth between aerial vehicle and measured target is

error can be obtained by comparison with the theoretical value.

1

phase interference method. Moreover, it is also not connected with wavelength.

1

1

*tg*

length

the axis of aerial vehicle.

**5.3 Simulation analysis** 

Fig.5-2 and fig.5-3 depict the relative calculation error curve with different baseline and radial distance. When azimuth angle tends to zero, the derived formula has irregularity. Obviously, the farther is the object distance, or the shorter is the baseline length, the smaller is the relative calculation error of formula. There occurs this occurrence provided that the incident wave is parallel in derivation.

Fig. 5-2. Relative calculation error: different baseline lengths.

Fig. 5-3. Relative calculation error: different radial distances.

Basic parameters:

= 0.1 Hz

 = 100 km v = 300 m/s λ = 0.15 m

σf

r 1

Fig. 5-4. Measurement errors: different baseline lengths.

σf

r 1

d2 = 5 m v = 300 m/s λ = 0.15 m

Basic parameters:

= 0.1 Hz

= 100 km

0 10 20 30 40 50 60 70 80 90

0 10 20 30 40 50 60 70 80 90

/ (°)

Azimuth angle θ<sup>1</sup>

Fig. 5-5. Measurement errors: different ratio between two baseline lengths.

/ (°)

d1 = d2

d1 = d2

d1 = d2

> d1 : d2 = 10

d1 : d2 = 5

d1 : d2 = 1

= 5 m

= 15 m

= 30 m

Azimuth angle θ<sup>1</sup>

10-1

100

101

Measurement error

σ / (°)

102

103

100

101

Measurement error

σ / (°)

102

103

104

105

### **5.4 Error analysis**

On condition that the locating error of baseline is neglected, making differential for Doppler frequency difference gives

$$\frac{\partial \theta\_1}{\partial \Delta f\_{d1}} = \frac{1}{1 + A^2} \frac{d\_2}{d\_1 \Delta f\_{d2}} \tag{5-7}$$

$$\frac{\partial \mathcal{O}\_1}{\partial \Delta f\_{d2}} = \left| \frac{1}{1 + A^2} \frac{d\_2 \Delta f\_{d1}}{d\_1 \Delta f\_{d2}^2} \right| \tag{5-8}$$

where: <sup>2</sup> <sup>1</sup> 1 2 *d d <sup>d</sup> <sup>f</sup> <sup>A</sup> d f*

According to error estimation theory, the overall error produced by the measurement of Doppler frequency difference is

$$
\sigma = \sigma\_f \sqrt{\sum\_{i=1}^{n=2} \frac{\widehat{\mathcal{O}\theta\_i}}{\widehat{\mathcal{O}\Delta} \xi\_{di}}} \tag{5.9}
$$

where: *<sup>f</sup>* is the mean-root-square error measuring Doppler frequency difference.

Fig.5-4 shows that the error curve falls exponentially. When azimuth angle tends to zero, there is obviously the blind zone for direction finding. And the accuracy is best when azimuth angle is close to right angle. The analysis shows that the measuring error can usefully decrease if the length of two baselines is simultaneously increased. For example, as soon as the length of two baselines can be increased to 15 meter, the maximum angle measuring error is less than <sup>0</sup> 2 after azimuth angle is bigger than <sup>0</sup> 10 .

Fig.5-5 depicts that the measurement error when the azimuth angle is less than <sup>0</sup> 30 can be quickly reduced by changing the specific value between two baselines.

Although the relative calculation error is not connected with wavelength, the calculation shows that the measurement accuracy is connected with wavelength. Fig.5-6 shows that the measurement accuracy is inversely proportional to the wavelength.

Moreover, though the relative calculation error is not connected with the flight speed, the calculation shows that the measurement accuracy is connected with the flight speed. Fig.5-7 shows that the measurement accuracy is in direct proportional to the flight speed.

### **5.5 Summarize briefly**

Since the new method presented in this text is not associated with wavelength, the direction finding only based on Doppler frequency difference will be more adapted to passive sounding as compared with phase interference method.

The multi-objective passive localization and wideband operation from single airborne observer is a key and difficulty subject for modern electron reconnaissance. Because the

On condition that the locating error of baseline is neglected, making differential for Doppler

1 2 2 1 12 1 1 

*d*

2 2 2 1 2 1 1 

*d f*

*<sup>f</sup> <sup>A</sup> d f* (5-7)

*<sup>f</sup> A d <sup>f</sup>* (5-8)

*<sup>i</sup> di <sup>f</sup>* (5-9)

*d d*

1 2 1

 *d d d*

According to error estimation theory, the overall error produced by the measurement of

2

*<sup>n</sup> <sup>i</sup>*

1

 

*f*

*<sup>f</sup>* is the mean-root-square error measuring Doppler frequency difference.

Fig.5-4 shows that the error curve falls exponentially. When azimuth angle tends to zero, there is obviously the blind zone for direction finding. And the accuracy is best when azimuth angle is close to right angle. The analysis shows that the measuring error can usefully decrease if the length of two baselines is simultaneously increased. For example, as soon as the length of two baselines can be increased to 15 meter, the maximum angle

Fig.5-5 depicts that the measurement error when the azimuth angle is less than <sup>0</sup> 30 can be

Although the relative calculation error is not connected with wavelength, the calculation shows that the measurement accuracy is connected with wavelength. Fig.5-6 shows that the

Moreover, though the relative calculation error is not connected with the flight speed, the calculation shows that the measurement accuracy is connected with the flight speed. Fig.5-7

Since the new method presented in this text is not associated with wavelength, the direction finding only based on Doppler frequency difference will be more adapted to passive

The multi-objective passive localization and wideband operation from single airborne observer is a key and difficulty subject for modern electron reconnaissance. Because the

shows that the measurement accuracy is in direct proportional to the flight speed.

 

measuring error is less than <sup>0</sup> 2 after azimuth angle is bigger than <sup>0</sup> 10 .

quickly reduced by changing the specific value between two baselines.

measurement accuracy is inversely proportional to the wavelength.

sounding as compared with phase interference method.

**5.4 Error analysis** 

where: <sup>2</sup> <sup>1</sup>

where:

**5.5 Summarize briefly** 

*<sup>d</sup> <sup>f</sup> <sup>A</sup> d f*

1 2 *d d*

Doppler frequency difference is

frequency difference gives

Fig. 5-4. Measurement errors: different baseline lengths.

Fig. 5-5. Measurement errors: different ratio between two baseline lengths.

direction finding method presented in this text is only associated with Doppler frequency difference, that is, which is not associated with the setting up the wavelength and baseline

At the same time, this method is also very applicable to the electron reconnaissance for multi-objective from single airborne observer because the information of frequency

Available localization method based on Doppler as well as its rate of change uses the coordinate variation in rectangular coordinate system for analyzing the angle of advance in polar coordinate system and uses geometric projection for analyzing the speed vector. Because the Doppler shift is the function of position and motion state of target, there are 2n unknown number in n-dimensional plane. For the localization of single station, in order to resolve the position and motion state of target, 2n nonlinear equations must be set up according to the measured value obtained during n measuring cycle. The solution process is

The method of Doppler direct ranging changes such status. Recently, research done by author shows that the direct range finding for target with analytic form can be realized making use of Doppler frequency difference during a time interval from the viewpoint of the basic definition of Doppler rate. More, the direct ranging formula whose characteristic is

Though the auctorial research result has proves that the single station which is moving can realize fast ranging for fixed or slow speed target making use of Doppler shift as well as its rate of change. However, the measurement accuracy for carrier frequency in passive

restricts the application of Doppler passive localization technique in engineering. As a rule, the value of Doppler frequency difference may be less than Hz since the baseline length of antenna array is shorter in single moving platform. This conduces that the demand for

In allusion to such difficulty, author researches the phase detection method for Doppler shift and presents the airborne passive phase interference ranging method based on the principle of Doppler direct ranging. The other one which is being explored by author is the DF method combining Doppler with phase interference in which the path length difference is firstly determined by directly solving the integer of wavelength in path length difference based on measuring Doppler as well as its rate. Though these methods are also imperfect, elementary result has shown that the measuring method of Doppler-phase interference will help to development of airborne passive localization technique as well as correlation

As soon as the phase shift can directly integrate with the Doppler frequency shift in physical relationship, more novel method of localization can be appeared. The mutual equivalent conversion between phase shift and Doppler shift means that the method obtained by a kind of measurement technology can be quickly expanded to another measurement technology.

*<sup>f</sup> f* based on actual technique level. This measuring ability

length, this method is very applicable to wideband operation.

measurement is advantageous for signal sorting and recognition.

comparative complicated and the analytic result can not be obtained.

better can be obtained by applying angle rate.

measurement accuracy of frequency difference is very rigor.

localization is just <sup>4</sup> 10

technique.

**6. Conclusions** 

Fig. 5-6. Measurement errors: different wavelengths.

Fig. 5-7. Measurement errors: different flight speeds.

direction finding method presented in this text is only associated with Doppler frequency difference, that is, which is not associated with the setting up the wavelength and baseline length, this method is very applicable to wideband operation.

At the same time, this method is also very applicable to the electron reconnaissance for multi-objective from single airborne observer because the information of frequency measurement is advantageous for signal sorting and recognition.

### **6. Conclusions**

164 Recent Interferometry Applications in Topography and Astronomy

Basic parameters:

= 5 m

= 0.1 Hz

= 100 km

v = 300 m/s

σf

r 1

d1 = d2

Fig. 5-6. Measurement errors: different wavelengths.

σf

r 1

Fig. 5-7. Measurement errors: different flight speeds.

d = d1 = d2

λ = 0.15 m

Basic parameters:

= 5 m

= 0.1 Hz

= 100 km

0 10 20 30 40 50 60 70 80 90

0 10 20 30 40 50 60 70 80 90

/ (°)

Azimuth angle θ<sup>1</sup>

/ (°)

λ = 0.15 m λ = 0.5 m λ = 1 m

v = 50 m/s v = 150 m/s v = 300 m/s

Azimuth angle θ<sup>1</sup>

100

100

101

Measurement error

σ / (°)

102

103

101

Measurement error

σ / (°)

102

103

Available localization method based on Doppler as well as its rate of change uses the coordinate variation in rectangular coordinate system for analyzing the angle of advance in polar coordinate system and uses geometric projection for analyzing the speed vector. Because the Doppler shift is the function of position and motion state of target, there are 2n unknown number in n-dimensional plane. For the localization of single station, in order to resolve the position and motion state of target, 2n nonlinear equations must be set up according to the measured value obtained during n measuring cycle. The solution process is comparative complicated and the analytic result can not be obtained.

The method of Doppler direct ranging changes such status. Recently, research done by author shows that the direct range finding for target with analytic form can be realized making use of Doppler frequency difference during a time interval from the viewpoint of the basic definition of Doppler rate. More, the direct ranging formula whose characteristic is better can be obtained by applying angle rate.

Though the auctorial research result has proves that the single station which is moving can realize fast ranging for fixed or slow speed target making use of Doppler shift as well as its rate of change. However, the measurement accuracy for carrier frequency in passive localization is just <sup>4</sup> 10 *<sup>f</sup> f* based on actual technique level. This measuring ability restricts the application of Doppler passive localization technique in engineering. As a rule, the value of Doppler frequency difference may be less than Hz since the baseline length of antenna array is shorter in single moving platform. This conduces that the demand for measurement accuracy of frequency difference is very rigor.

In allusion to such difficulty, author researches the phase detection method for Doppler shift and presents the airborne passive phase interference ranging method based on the principle of Doppler direct ranging. The other one which is being explored by author is the DF method combining Doppler with phase interference in which the path length difference is firstly determined by directly solving the integer of wavelength in path length difference based on measuring Doppler as well as its rate. Though these methods are also imperfect, elementary result has shown that the measuring method of Doppler-phase interference will help to development of airborne passive localization technique as well as correlation technique.

As soon as the phase shift can directly integrate with the Doppler frequency shift in physical relationship, more novel method of localization can be appeared. The mutual equivalent conversion between phase shift and Doppler shift means that the method obtained by a kind of measurement technology can be quickly expanded to another measurement technology.

McCormick W S, Tsui J B Y, BakkieV L. A Noise Insensitive Solution to an Ambiguity

Messer H,Singal G. On the achievable DF accuracy of two kinds of active interferometers.

Poisel,Richard A. (2008). Electronic warfare target location methods. Translated by Qu Xiaoxu. *Electronic Industry Press,* ISBN 978-7-121-06326-8, Beijing, China

Sun Zhongkang, Guo Fucheng, Feng Daowang. (2008). Passive location and tracking

Sundaram K R, Ranjan K M. Modulo conversion method for estimation the direction of

Wang Qiang. A Research on Passive Location and Tracking Technology of a Single Airborne

Wang Zhi-rong. On High Precision DOA Estimation. *Radio Engineering of China*,.Vol.37,

Wei Xing & Wan Jian-wei & Huang Fu-kan. Study of Passive Location System Based on

Xu Yao-wei & Sun Zhong-kang. Passive Location of Fixed Emitter Using Phase Rate of

Yang Yue-lun. Summary of New Airborne Passive Detection and Location Technologies,

Yu Chun1ai, Wan Jianwei, Zhan Ronghui. An Estimation Algorithm for Doppler Frequency

Yu Tao, An Analytic Method for Doppler Changing Rate Based on Frequency

Yu Tao. Airborne Direction Finding Method Based on Doppler-Phase Measuremen. *Frontiers* 

Yu Yang. The Improvement of Airborne DF Adjust Method Based on Four-channel

Zhang Gangbing & Liu Yu & Liu Zongmin. Unwrapping Phase Ambiguity Algorithm Based

Intelligence Information Security. Beijing, China, December 17-19, 2010 Yu Tao. A new airborne passive DF method only based on frequency difference. *Advanced* 

*materials research*, Vols.219-220, (2011), pp.846-850. ISSN 1022-6680

*Technology*, Vol.30, No.10, (2008), pp. 2303-2306, ISSN 1009-5896

Vol.25, No.5, (1989), pp.729-732, ISSN 0018-9251

ISSN 0018-9251

05817-8, Beijing, China.

PR.China. (Nov. 2004).

ISSN 1004-7859

1001-506X

9167

3584

pp.1391-1396, ISSN 0018-9251

No.11, (2007), pp.24-25, ISSN 1003-3106

(2007),pp.24-26,39, ISSN 1674-2230

No.5, (2008), pp.665-669, ISSN 1005-2615

Problem in Spectral Estimation. *IEEE Trans on Aerospace and Electronic Systems*,

*IEEE Trans. on Aerospace and Electronic Systems*, Vol.32, No.3, (1996), pp.1158-1164,

technology by single observer. *National Defense Industry Press*, ISBN 978-7-118-

arrival. *IEEE Trans. on Aerospace and Electronic Systems*, Vol.36, No.4, (2000),

Observer. National University of Deefnse Technology, Changsha, Hunan.

Multi Base-line Interferometers. *Modern Radar*,.Vol.29, No.5, (2007), pp.22-25,35,

Change. *Systems Engineering and Electronics*, Vol.21, No.3, (1999), pp. 34-37, ISSN

*Shipboard Electronic Countermeasure*, Vol.34, No.3, (2011), pp.14-17+56, ISSN 1673-

Rate-of-Change with PCM Coherent Pulse Train. *Journal of Electronics & Information* 

Measurement[C]. 2010 International Conference on Communications and

*of Electrical and Electronic Engineering in China*, Vol.5, No.4, pp.493-495, ISSN 1673-

Amplitude Comparison. *Electronic Information Warfare Technology*, Vol.22, No.3,

on Baseline Ratio. *Journal of Nanjing University of Aeronautics & Astronautics*, Vol.40,

Hence, the research for passive localization by applying the functional relation between phase shift and Doppler shift may be a more meaningful work.

### **7. References**


Hence, the research for passive localization by applying the functional relation between

Abatzoglou T. Fast maximum likelihood joint estimation of frequency and frequency rate.

An Xiao-jun. Research on DF Based on Improved Phase Interferometer. *Radio Engineering of* 

Diao ming, Wang Yue. Research of passive location based on the Doppler changing rate.

Djuric Petar M and Kay Steven M. Parameter estimation of chirp signals. *IEEE Trans. on* 

Ernest J. Ambiguity resolution in interferometer. *IEEE Trans. Aerospace and Electronic* 

Feng Daowang & Zhou Yiyu & Li Zonghua. A Fast and Accurate Estimation for Doppler

Gong Xiang-yi & Yuan Jun-quan & Su Ling-hua. A Multi-pare Unwrap Ambiguity of

Hua Yang. Research into the Single-station Passive Location for Electronic Reconnaissance

Li Yong & Zhao Guo-wei & Li Tao. Algorithm of Solving Interferometer Phase Difference

Li Zonghua, Xiao Yuqin, Zhou Yiyu, Sun Zhongkang. Single-observer passive location and

Li Zonghua,Xiao Yuqin,Zhou Yiyu,et al. Single-observer passive location and tracking

Liao Ping & Yang Zhong-hai & Jiang Dao-an. A Single Observer Passive Location Algorithm

Lin Yi-meng & Liu Yu & Zhang Ying-nan. Algorithm of Direction Finding for Broadband

Lu Xiaomei. Overview on the Technology of the Single Airborne Station Passive Location, *Shipboard Electronic Countermeasure*, Vol.26, No.3, (2003), pp:20-23, ISSN 1673-9167

*and Electronics*, Vol.26, No.5, (2004), pp. 613-616, ISSN 1001-506X

*Electronics*, Vol.26, No.5, (2004), pp. 613-616, ISSN 1001-506X

*Information Technology*, Vol.28, No.1, (2006), pp.55-59, ISSN 1009-5896 Gong Xiang-yi & Huang Fu-kan & Yuan Jun-quan. A New Algorithm for Estimation of

*IEEE Trans. On Aerospace and Electronic Systems*, AES-Vol.22, No.6, (1986), pp.708-

*Systems Engineering and Electronics*, Vol28, No.5, (2006), pp. 696-698, ISSN 1001-506X

*Acoustics*, *Speech and Signal Processing*, Vol.38, No.12, (1990), pp. 2118-2126, ISSN

Rate-of-Change with the Coherent Pulse Train. *Signal Processing,* Vol.20, No.1,

Interferometer Array for Estimation of Direction of Arrival. *Journal of Electronics &* 

Direction of Arrival Based on the Second-Order Difference of Phase of Interferometer Array. *Acta Electronica Sinica*,.Vol.33, No.3, (2005), pp.444-446, ISSN

UAV. *Shipboard Electronic Countermeasure*, Vol.32, No.2, (2008), pp.14-17, ISSN 1673-

Ambiguity by Airborne Single Position. *Chinese Journal of Sensors and Actuators*,

tracking algorithms using frequency and spatial measurements. *Systems Engineering* 

algorithms using frequency and spatial measurements. *Systems Engineering and* 

Based on Probability in Multi-target Environment, *Telecommunication Engineering*.

Digital Signal. *Journal of Nanjing University of Aeronautics & Astronautics*, Vol.37,

phase shift and Doppler shift may be a more meaningful work.

*China*, Vol.39, No.03, (2009), pp.59-61, ISSN 1003-3106

*Systems*, Vol.117, No.6, (1981), pp.766-780, ISSN 0018-9251

Nol.19, No.6, (2006), pp.2600-2602, ISSN 1004-1699

Vol.46, No.1, (2006), pp.45-49, ISSN 1001-896X

No.3, (2005), pp.335-340, ISSN 1005-2615

**7. References** 

715, ISSN 0018-9251

(2004), pp.40-43, ISSN 1003-0530

1053-587X

0372-2112

9167


**8** 

**Experiences in Boreal Forest Stem Volume** 

Jan Askne1 and Maurizio Santoro2 *1Chalmers University of Technology,* 

*2Gamma Remote Sensing,* 

*1Sweden 2Switzerland* 

**Estimation from Multitemporal C-Band InSAR** 

During the last two decades synthetic aperture radar (SAR) and interferometric synthetic aperture radar (InSAR) have become important tools for airborne and satellite remote sensing. SAR delivers a high resolution image of the radar backscatter and by means of InSAR a second image from almost the same orbit is combined with the first one to provide the relative phase difference and stability/coherence of the backscattered signal. Principles of InSAR can be found in various review articles (Bamler and Hartl, 1998; Massonnet and Feigl, 1998; Rott, 2009). This presentation will focus on the use of InSAR observations for

Satellite methods are important to determine forest above-ground biomass on a global scale in order to follow changes in a consistent manner over long periods of time as part of climate modeling. Methods to determine forest stem volume are also important for the economical management of forest areas. For these applications many remote sensing techniques have been investigated. Our interest will focus on a radar method. In contrast to methods based on optical images, a semi-empirical, physically-based relation between measurements and stem volume/biomass can be established because microwaves penetrate the forest canopy to a certain extent and the backscattered signal is therefore modulated by

With the early C-band SAR systems, the European Remote Sensing Satellites ERS-1 and ERS-2, and in particular the observations during the tandem period with one day interval between the observations, large amounts of InSAR data became available, which allowed the possibility to retrieve forest parameters. One argument against C-band for forest applications has been that the backscatter is originating from the top layer of the forest with little relation to the major components of the biomass, such as the stem and then saturate for low biomass values (Imhoff, 1995). However, this changed with repeat pass InSAR data, and coherence can under certain conditions provide high accuracy stem volume estimates up to the stem volumes available in the test areas studied below, i.e. up to 539 m3/ha or approximately 265 tons/ha. Microwave penetration into forests is related to the vegetation properties including the gaps in the vegetation. Since the density of boreal forests is often low this is a fundamental property for the application of C-band and higher frequencies,

investigation of forest stem volume or above ground biomass.

**1. Introduction** 

the forest structural properties.

