A Unified Algorithm for Channel Imbalance and Antenna Phase Center Position Calibration of a Single-Pass Multi-Baseline TomoSAR System

Abstract

The multi-baseline synthetic aperture radar (SAR) tomography (TomoSAR) system is employed in such applications as disaster remote sensing, urban 3-D reconstruction, and forest carbon storage estimation. This is because of its 3-D imaging capability in a single-pass platform. However, a high 3-D resolution of TomoSAR is based on the premise that the channel imbalance and antenna phase center (APC) position are precisely known. If this is not the case, the 3-D resolution performance will be seriously degraded. In this paper, a unified algorithm for channel imbalance and APC position calibration of a single-pass multi-baseline TomoSAR system is proposed. Based on the maximum likelihood method, as well as the least squares and the damped Newton method, we can calibrate the channel imbalance and APC position. The algorithm is suitable for near-field conditions, and no phase unwrapping operation is required. The effectiveness of the proposed algorithm has been verified by simulation and experimental results.

1. Introduction

In recent years, synthetic aperture radar (SAR) tomography (TomoSAR) has become a popular research topic due to its 3-D imaging capability [1,2,3,4,5,6]. TomoSAR has been successfully applied in many application contexts, such as forestry [7,8], 3D urban reconstruction [9,10], and glaciers [11]. The single-pass multi-baseline TomoSAR system has 3D resolution including the height resolving ability in a single-pass platform because there are multiple channels in the cross-track direction. However, the Rayleigh resolution in the height direction is very limited due to the limitation of the baseline length in a single-pass platform. As shown in Equation (25), the Rayleigh resolution of a TomoSAR system is about 35 m. If super-resolution performance is desired in the height direction, at a resolution of, say, 5 m, then a super-resolution algorithm must be introduced in that direction so as to distinguish multiple targets in small intervals within the Rayleigh resolution. Tebaldini [11] concludes that even a subwavelength accuracy of the antenna phase center (APC) position will hinder the focusing result in the height direction. The authors of [12] conclude that, when one wants to distinguish multiple point-like targets with different heights within a slant range-azimuth resolution cell, the requirements for phase stability or phase calibration accuracy are higher than those for traditional InSAR. Therefore, if we hope to obtain a super-resolution performance in the height direction, high requirements for the channel imbalance (also known as amplitude and phase inconsistency) and APC position calibration are required in a single-pass multi-baseline TomoSAR system.

Several calibration algorithms have been proposed in recent years. Pardini [13] calibrated the phase error due to APC position by an algorithm based on a minimum entropy criterion. In [14], Gocho compensates for the phase screen caused by APC position errors using eigenvalue decomposition and phase interpolation/extrapolation. However, these two methods cannot obtain the exact APC position, and they fail to take channel imbalance into account. A phase center double localization algorithm is proposed by Tebaldini [11,15] to obtain the APC position. However, this method loses its effectiveness when there is channel imbalance. A different method has been proposed in [16], which estimates the APC position based on raw data. Such processing seems infeasible for users with single-look complex (SLC) images only. Regarding the correction of amplitude and phase errors in multi-channel array systems, Kuoye Han [17] proposes a calibrator-based approach. On the basis of [17], Xiaolin Yang takes into account both channel imbalance and APC position in [18], but they failed to decouple the phase errors caused by the channel itself and the APC position. In general, these calibration algorithms above either do not consider both the channel imbalance and the APC position or are unable to decouple channel imbalance and the APC position error. In addition, a plane wave model is used in most of these algorithms. However, even in the far-field condition, the plane wave model will bring some non-negligible phase error when the ratio of the baseline length to line-of-sight distance is not sufficiently close to zero [19]. As an example, Figure 6 will show that the plane wave model brings a non-negligible phase error even in the far-field condition, where b = 0.6 m and r = 1625 m. Therefore, we use the Fresnel approximation, which assumes that spherical waves can be approximated by quadratic waves in our signal model.

In this paper, we propose a unified algorithm for the channel imbalance and APC position calibration of a single-pass multi-baseline TomoSAR system, which can not only calibrate channel imbalance but can also calibrate the APC position. Features of this algorithm are as follows: (1) the channel imbalance and APC position can be calibrated individually rather than confusing the phase error caused by the APC position error with the phase error of the channel itself; (2) the Fresnel approximation is used in the calibration signal model, which heightens the accuracy of the calibration signal model; and (3) there is no need for phase unwrapping.

This paper is structured as follows. Section 2 is devoted to establishing a signal model. The proposed calibration algorithm is described in Section 3. The effectiveness of the calibration algorithm is validated in Section 4 with a real data set acquired by the array InSAR system [20,21] developed by the Institute of Electronics, Chinese Academy of Sciences (IECAS). The discussion and conclusions follow in Section 5 and Section 6, respectively.

2. A Signal Model

Supposing that the TomoSAR system is working in side-looking mode, the TomoSAR acquisition geometry is depicted in Figure 1. Assuming that the number of the APC is N, axes x, y, and z are the cross-track, azimuth, and height directions, respectively. APC is supposed in the zero Doppler plane (see Figure 2), where s is the cross-range direction. Let $\left(x_n,z_n\right)$ denote the position of the nth APC; without loss of generality, we assume that APC1 is at the origin of the coordinate system, that is, APC1 is the reference APC and $\left(x_1,z_1\right)=\left(0,0\right)$. Applying a classical imaging algorithm to the raw SAR data collected in each channel, we obtain N 2-D SAR images, usually referred to as SLC images. After some sub-pixel accuracy coregistration to the reference channel, and under the Born weak-scattering approximation, the focused complex value of an azimuth-range pixel $\left(y_0,r_0\right)$ of the nth channel is [1]:

$${\stackrel{⌢}{\gamma }}_n\left(y_0,r_0\right)={\displaystyle \iint dydrf\left(y_0-y,r_0-r\right)}{\displaystyle \int ds\gamma \left(y,r,s\right)\exp \left[-j\frac{4\pi }{\lambda }{R}_n\left(r,s\right)\right]}$$

(1)

where $\lambda$ is the wavelength, $f\left(y_0-y,r_0-r\right)$ is the 2-D point spread function (PSF), $\gamma \left(y,r,s\right)$ is the function that models the 3-D scene scattering properties, and $R_n\left(r,s\right)$ represents the slant range between point targets located at $r$ and $s$ coordinates, and the nth APC. Under the Fresnel approximation, $R_n\left(r,s\right)$ can be written as

$$R_n\left(r,s\right)=\sqrt{{\left(r-b_{//n}\right)}^2+{\left(b_{\perp n}-s\right)}^2}\approx \left|r-b_{//n}\right|+\frac{{\left(b_{\perp n}-s\right)}^2}{2\left|r-b_{//n}\right|}$$

(2)

where $b_{//n}$ is the nth horizontal baseline, $b_{\perp n}$ is the nth orthogonal baseline, the relation between the baseline and APC position is shown as follows:

$$\{\begin{array}{l}{b}_{\perp n}=x_n\cos \theta +z_n\sin \theta \\ b_{//n}=x_n\sin \theta -z_n\cos \theta \end{array}$$

(3)

where $\theta$ is the off-nadir angle.

Based on the assumption of the point target, $f\left(y_0-y,r_0-r\right)$ can be regarded as the 2-D Dirac function. Following this, Equation (1) can be written as

$${\stackrel{⌢}{\gamma }}_n\left(y_0,r_0\right)={\displaystyle \int ds\gamma \left(y_0,r_0,s\right)\exp \left[-j\frac{4\pi }{\lambda }{R}_n\left(r_0,s\right)\right]}.$$

(4)

In the calibration processing, the ground control point (GCP) is usually placed in the non-layover area (e.g., the bare ground). That is, there is only one target in the s direction within a slant range-azimuth cell corresponding to the GCP. Therefore, Equation (4) can be rewritten as

$${\stackrel{⌢}{\gamma }}_n\left(y_0,r_0\right)=\exp \left[-j\frac{4\pi }{\lambda }{R}_n\left(r_0\right)\right]\gamma \left(y_0,r_0\right)$$

(5)

where

$$R_n\left(r_0\right)=\left|r_0-b_{//n}\right|+\frac{b_{\perp n}{}^2}{2\left|r_0-b_{//n}\right|}.$$

(6)

For notational convenience, the target position coordinates $r_0$ and $y_0$ in each symbol are omitted, and the focused complex value acquired by the nth channel is denoted as $g_n$. Considering that $r_0>>b_{//n}$, we have $R_n=r_0-b_{//n}+\frac{b_{\perp n}{}^2}{2r_0}$; therefore, Equation (5) can be rewritten as

$$g_n=\exp \left(-j\frac{4\pi }{\lambda }{R}_n\right)\gamma .$$

(7)

By combining the exponential term of the reference channel into the backscatter coefficient, and by taking into account channel imbalance, the calibration signal model of the multi-baseline TomoSAR system would be of the form:

$$\mathbf{g}=\mathbf{C}\mathsf{\alpha }{\gamma }^{\prime }+\mathbf{E}$$

(8)

where $\mathbf{g}={\left[g_1,g_2,...,g_N\right]}^T$ is the $N\times 1$ observation vector, the calibration matrix $\mathbf{C}=diag\left\{{\rho }_1{e}^{j{\phi }_1},...,{\rho }_n{e}^{j{\phi }_n},...,{\rho }_N{e}^{j{\phi }_N}\right\}$, where ${\rho }_n$ and ${\phi }_n$ are the amplitude and phase of the nth channel (note that, if all channels are exactly the same, then $\mathbf{C}=\mathbf{I}$), represents the channel imbalance, $\mathsf{\alpha }={\left[1,e^{-j\frac{4\pi }{\lambda }\left(R_2-R_1\right)},...,e^{-j\frac{4\pi }{\lambda }\left(R_N-R_1\right)}\right]}^T$, and ${\gamma }^{\prime }=\gamma \exp \left(-j\frac{4\pi }{\lambda }{R}_1\right)$, which is related to the backscatter coefficient of GCP.

3. The Calibration Algorithm

Let ${\theta }_m$ ($m=1,2,...,M$) be the off-nadir angle of the mth GCP, and $r_m$ be the slant range. Then,

$$\mathbf{g}=\mathbf{C}\mathsf{\alpha }\left(\mathsf{\psi },{\theta }_m\right){\gamma }_m+\mathbf{n}$$

(9)

where $\mathbf{g}={\left[g_1,g_2,...,g_N\right]}^T$is an observation vector of the mth GCP, ${\gamma }_m$ is the backscatter coefficient, $\mathsf{\alpha }\left(\mathsf{\psi },{\theta }_m\right)$ is an array manifold, and $\mathsf{\psi }={\left[x_1,x_2,...,x_N,z_1,z_2,...,z_N\right]}^T$ is a vector of the unknown APC position. By executing an eigenvalue decomposition of the covariance matrix $E\left[\mathbf{g}{\mathbf{g}}^T\right]$, and normalizing the first element of the eigenvector corresponding to the largest eigenvalue to 1, we obtain the array manifold estimation ${\mathsf{\alpha }}_{mea}\left({\theta }_m\right)$, and

$${\mathsf{\alpha }}_{mea}\left({\theta }_m\right)=\mathbf{C}\mathsf{\alpha }\left(\mathsf{\psi },{\theta }_m\right)+{\mathbf{n}}_m$$

(10)

where ${\mathbf{n}}_m\in {ℂ}^{N\times 1}$ is a random vector of additive noise. We shall assume that ${\mathbf{n}}_m$ has a zero- mean Gaussian distribution with a covariance ${\zeta }^2$.

Considering all M GCPs, the probability density function of the set of estimated array manifold is

$$\begin{array}{l}p\left({\mathsf{\alpha }}_{mea}\left({\theta }_1\right),{\mathsf{\alpha }}_{mea}\left({\theta }_2\right),...,{\mathsf{\alpha }}_{mea}\left({\theta }_M\right)|\mathbf{C},\mathsf{\psi },{\zeta }^2\right)\\ ={\left(\pi {\zeta }^2\right)}^{-MN}\exp \left\{-\frac{1}{{\zeta }^2}{\displaystyle \sum _{m=1}^M{\Vert {\mathsf{\alpha }}_{mea}\left({\theta }_m\right)-\mathbf{C}\mathsf{\alpha }\left(\mathsf{\psi },{\theta }_m\right)\Vert }^2}\right\}\end{array}.$$

(11)

It can be easily shown that the maximum likelihood estimates of $\mathbf{C}$ and $\mathsf{\psi }$ are the corresponding values that minimize the following cost function:

$$\underset{\mathbf{C},\mathsf{\psi }}{\min }{\Vert \mathbf{C}\mathbf{A}\left(\mathsf{\psi }\right)-{\mathbf{A}}_m\Vert }_F^2$$

(12)

where ${\Vert \cdot \Vert }_F$ denotes the Frobenius norm, ${\mathbf{A}}_m\triangleq \left[{\mathsf{\alpha }}_{mea}\left({\theta }_1\right),{\mathsf{\alpha }}_{mea}\left({\theta }_2\right),...,{\mathsf{\alpha }}_{mea}\left({\theta }_M\right)\right]$, and $\mathbf{A}\left(\mathsf{\psi }\right)\triangleq \left[\mathsf{\alpha }\left(\mathsf{\psi },{\theta }_1\right),\mathsf{\alpha }\left(\mathsf{\psi },{\theta }_2\right),...,\mathsf{\alpha }\left(\mathsf{\psi },{\theta }_M\right)\right]$. Since the calibration matrix $\mathbf{C}$ and the APC position contain $2N^2$ and $2N$ real unknowns, respectively, and since the M sources can provide $2MN$ independent measurements, $M\ge N+1$ is a necessary condition. This condition can also be found in [22]. According to this condition, at least 9 GCPs are needed to meet the requirement for our 8-channel array InSAR system.

Assuming that $\mathbf{A}\left(\mathsf{\psi }\right)$ is of full rank N, the cost equation, Equation (12), is a separable nonlinear least-squares optimization problem. When we keep $\mathsf{\psi }$ fixed, the least-square estimation of the calibration matrix $\mathbf{C}$ is

$$\mathbf{C}={\mathbf{A}}_m{\mathbf{A}}^H\left(\mathsf{\psi }\right){\left(\mathbf{A}\left(\mathsf{\psi }\right){\mathbf{A}}^H\left(\mathsf{\psi }\right)\right)}^{-1}.$$

(13)

Substituting Equation (13) into Equation (12), we obtain the maximum likelihood estimation for the APC position:

$${\widehat{\mathsf{\psi }}}_{ML}=\arg \min Tr\left(P_{{\mathbf{A}}^H\left(\mathsf{\psi }\right)}^{\perp }{\mathbf{A}}_m^H{\mathbf{A}}_m\right)=\arg \min f\left(\mathsf{\psi }\right)$$

(14)

where $P_{{\mathbf{A}}^H\left(\mathsf{\psi }\right)}^{\perp }=I-{\mathbf{A}}^H\left(\mathsf{\psi }\right){\left(\mathbf{A}\left(\mathsf{\psi }\right){\mathbf{A}}^H\left(\mathsf{\psi }\right)\right)}^{-1}\mathbf{A}\left(\mathsf{\psi }\right)$, and the minimization objective function $f\left(\mathsf{\psi }\right)\triangleq Tr\left(P_{{\mathbf{A}}^H\left(\mathsf{\psi }\right)}^{\perp }{\mathbf{A}}_m^H{\mathbf{A}}_m\right)$.

The estimation problem is now decoupled into two steps. The APC position estimation ${\widehat{\mathsf{\psi }}}_{ML}$ is obtained at first by solving the optimization problem Equation (14) before the estimation of the calibration matrix $C$ is derived by substituting ${\widehat{\mathsf{\psi }}}_{ML}$ estimated in the first step of Equation (13). We can then obtain the channel imbalance estimation by extracting the diagonal elements of the calibration matrix $\mathbf{C}$.

3.1. Step 1: APC Position Calibration

In this step, the APC position estimation ${\widehat{\mathsf{\psi }}}_{ML}$ is obtained by solving the optimization problem (Equation (14)). The minimization objective function $f\left(\mathsf{\psi }\right)$ may have many local minimums. However, when the nominal values of the APC position ${\mathsf{\psi }}_0$ are close enough to the true values, we can rewrite the minimization objective function as

$$f\left({\mathsf{\psi }}_0+\mathbf{p}\right)\approx f\left({\mathsf{\psi }}_0\right)+{\mathbf{g}}^t\left({\mathsf{\psi }}_0\right)\mathbf{p}+\frac{1}{2}{\mathbf{p}}^t\mathbf{H}\left({\mathsf{\psi }}_0\right)\mathbf{p}$$

(15)

where $\mathbf{g}$ and $\mathbf{H}$ are the gradient and Hessian of $f\left(\mathsf{\psi }\right)$, respectively. $\mathbf{p}$ is the search direction that minimizes the right hand of Equation (15).

In order to solve this optimization problem, the damped Newton method can be applied. The nominal values were set as the initial APC position ${\mathsf{\psi }}_0$, and the search direction ${\mathbf{p}}_{k+1}={\widehat{\mathsf{\psi }}}_{k+1}-{\widehat{\mathsf{\psi }}}_k=-{\mu }_k{\mathbf{H}}_k^{-1}{\mathbf{g}}_k$, where ${\mathbf{H}}_k\triangleq {\mathbf{H}|}_{\mathsf{\psi }={\widehat{\mathsf{\psi }}}_k}$ and ${\mathbf{g}}_k\triangleq {\mathbf{g}|}_{\mathsf{\psi }={\widehat{\mathsf{\psi }}}_k}$. The step length ${\mu }_k={\left(0.5\right)}^l$, where $l$ is the smallest nonnegative integer that satisfies $f\left({\widehat{\mathsf{\psi }}}_{k+1}\right)<f\left({\widehat{\mathsf{\psi }}}_k\right)$. Then the estimation of APC position is obtained when the damped Newton method converges or reaches the maximum number of iterations.

The solution is re-derived in the Appendix, and the required $\mathbf{g}$ and $\mathbf{H}$ are shown directly to be

$$\mathbf{g}=-2\mathrm{Re}\left\{\begin{array}{l}vecd\left({\mathbf{A}}_{\mathbf{x}}{\mathbf{P}}_{{\mathbf{A}}^H\left(\mathsf{\psi }\right)}^{\perp }{\mathbf{A}}_m^H{\mathbf{A}}_m{\mathbf{A}}^H\left(\mathsf{\psi }\right){\left(\mathbf{A}\left(\mathsf{\psi }\right)A^H\left(\mathsf{\psi }\right)\right)}^{-1}\right)\\ vecd\left({\mathbf{A}}_{\mathbf{z}}{\mathbf{P}}_{{\mathbf{A}}^H\left(\mathsf{\psi }\right)}^{\perp }{\mathbf{A}}_m^H{\mathbf{A}}_m{\mathbf{A}}^H\left(\mathsf{\psi }\right){\left(\mathbf{A}\left(\mathsf{\psi }\right){\mathbf{A}}^H\left(\mathsf{\psi }\right)\right)}^{-1}\right)\end{array}\right\}$$

(16)

$$\mathbf{H}=2\mathrm{Re}\left\{\mathbf{D}{\mathbf{P}}_{{\mathbf{A}}^H\left(\mathsf{\psi }\right)}^{\perp }{\mathbf{D}}^H\odot {\left({\mathbf{1}}_{2\times 2}\otimes \mathbf{E}\right)}^T\right\}$$

(17)

where

$${\mathbf{A}}_{\mathbf{x}}=j\frac{4\pi }{\lambda }\mathbf{A}\left(\mathsf{\psi }\right){\Lambda }_{\sin }-j\frac{4\pi }{\lambda }\mathbf{A}\left(\mathsf{\psi }\right)\odot \left[\begin{array}{l}{x}_1,z_1\\ x_2,z_2\\ \begin{array}{cc}⋮& ⋮\end{array}\\ x_N,z_N\end{array}\right]\left[\begin{array}{l}\frac{{\cos }^2{\theta }_1}{r_1},\frac{{\cos }^2{\theta }_2}{r_2},...,\frac{{\cos }^2{\theta }_M}{r_M}\\ \frac{\sin {\theta }_1\cos {\theta }_1}{r_1},\frac{\sin {\theta }_2\cos {\theta }_2}{r_2},...,\frac{\sin {\theta }_M\cos {\theta }_M}{r_M}\end{array}\right]$$

(18)

$${\mathbf{A}}_{\mathbf{z}}=-j\frac{4\pi }{\lambda }\mathbf{A}\left(\mathsf{\psi }\right){\Lambda }_{\cos }-j\frac{4\pi }{\lambda }\mathbf{A}\left(\mathsf{\psi }\right)\odot \left[\begin{array}{l}{x}_1,z_1\\ x_2,z_2\\ \begin{array}{cc}⋮& ⋮\end{array}\\ x_N,z_N\end{array}\right]\left[\begin{array}{l}\frac{\sin {\theta }_1\cos {\theta }_1}{r_1},\frac{\sin {\theta }_2\cos {\theta }_2}{r_2},...,\frac{\sin {\theta }_M\cos {\theta }_M}{r_M}\\ \frac{{\sin }^2{\theta }_1}{r_1},\frac{{\sin }^2{\theta }_2}{r_2},...,\frac{{\sin }^2{\theta }_M}{r_M}\end{array}\right]$$

(19)

$${\Lambda }_{\sin }=diag\left\{\left[\sin {\theta }_1,\sin {\theta }_2,...,\sin {\theta }_M\right]\right\}$$

(20)

$${\Lambda }_{\cos }=diag\left\{\left[\cos {\theta }_1,\cos {\theta }_2,...,\cos {\theta }_M\right]\right\}$$

(21)

$$D={\left[{\mathbf{A}}_{\mathbf{x}}^T,{\mathbf{A}}_{\mathbf{z}}^T\right]}^T$$

(22)

$$\mathbf{E}={\left(\mathbf{A}\left(\mathsf{\psi }\right){\mathbf{A}}^H\left(\mathsf{\psi }\right)\right)}^{-1}{\mathbf{A}}^H\left(\mathsf{\psi }\right){\mathbf{A}}_m^H{\mathbf{A}}_m{\mathbf{A}}^H\left(\mathsf{\psi }\right){\left(\mathbf{A}\left(\mathsf{\psi }\right){\mathbf{A}}^H\left(\mathsf{\psi }\right)\right)}^{-1}$$

(23)

where vecd(V) represents a vector formed from the diagonal elements of the matrix V, the symbol $\odot$ represents the matrix multiplication of the elements, $\otimes$ represents the Kronecker product, and $1_{p\times q}$ represents the $p\times q$ matrix with all entries equal to one.

3.2. Step 2: Channel Imbalance Calibration

After solving the APC position calibration problem, Equation (13) is solved in order to retrieve the channel imbalance. The diagonal element of the calibration matrix is the estimation of the channel imbalance. Non-diagonal elements are actually mutual coupling factors. In general, the mutual coupling factor is relatively small if the channels are spaced out far enough. In particular, the correction matrix will have some characteristics when the multi-channels array has some regular geometry. For example, the calibration matrix of the equidistant line array is the Toeplitz matrix.

3.3. Validation with Simulation Data

The simulation data is necessary since it makes it possible to directly compare the true channel imbalance and the APC position with those yielded by the calibration algorithm. This test would be quite hard to implement using real data, for which the true channel imbalance and APC position are, in general, not known with sufficient precision. The simulation data set consists of 8-channel 2-D focused SAR SLC images. Eight APCs corresponding to the eight-channel were distributed in the cross-track direction, and the longest baseline is 0.6 m. In this section, we will begin by showing a special case of the simulation data set in order to illustrate the effectiveness of the proposed algorithm. A Monte Carlo simulation is then carried out to assess the performance of the proposed algorithm in a statistical framework.

In the special case, eight APC trajectories—equally spaced out in the cross-track direction—are shown in Figure 3. Amplitude and phase inconsistency are set as AmpErr = [1, 1, 1, 1, 1, 1, 1, 1] and PhaseErr = [0, 0.3, 0.1, −0.2, 0.3, 0.1, 1, 0.4] rad. Other relevant system parameters are summarized in

Table 1. System parameters used in the simulation data set.

Simulation Data Parameters
Frequency	15 GHz (Ku-Band)
Bandwidth	500 MHz
PRF	2000 Hz
Airplane altitude	1000 m
Airplane Velocity	60 m/s
Azimuth beam width	2°
Depression angle	25°–41°

.

The observation scene includes $3\times 11$ non-layover GCPs, and four pairs of layover GCPs spaced out with super-resolution factors of 0.5, 1, 2, and 8, respectively. The super-resolution factor [4] is defined as the ratio of the Rayleigh resolution and layover scatters’ interval. The observation scene’s 2-D SAR image of the simulation data set is depicted in Figure 4. When taking a look at the 2-D SAR image in Figure 4, one can see that all non-layover GCPs are included in the imaging area and are well focused, and that four pairs of layover GCPs in red circles are appropriate superpositions, as expected.

After being coregistered to Channel 1 (the reference channel) with sub-pixel accuracy, the calibration algorithm proposed in this paper was carried out by utilizing all 33 non-layover GCPs. After calibration, we obtained the calibration values of the channel imbalance and APC position (see Figure 5). The true values of the channel imbalance and APC position are compared with the calibration values in order to verify the validity of the calibration algorithm. In Figure 5, the true values and calibration values are drawn in red and blue lines, respectively. Graphs in the left panel show the true and calibration values in contrast, while the differences (estimation errors) between them are plotted in the right panel. One can immediately note that the APC position estimation results are very accurate. The maximum and standard deviation are 0.16 mm and 0.105 mm, respectively. The maximum error of amplitude is only −30 dB. The maximum and standard deviation of the phase error are 0.12 rad and 0.06 rad, respectively. This seems a little disappointing because the maximum phase error is not as highly accurate as expected. An additional step, which can be helpful for improving the phase calibration accuracy, involves checking the same GCP in all channels or averaging the results of multiple GCPs, because we have obtained the exact APC position after the previous calibration.

In order to examine the effectiveness of this calibration values, a comparison experiment was performed on four pairs of layover GCPs. The true values and the estimation results of the target number and height are listed in

Table 2. Height resolution experiment on four pairs of layover GCPs.

Layover GCPs	Number of Targets			Height of Targets (m)
	True Value	Before Calibration	After Calibration	True Value	Before Calibration	After Calibration
1st pair	2	3	2	0, 56.9	9.7, 138.2, 166.9	0, 56.9
2nd pair	2	3	2	0, 32.7	0, 34, 143	0, 32.7
3rd pair	2	3	2	0, 12.4	8.4, 34.9, 78.1	0, 12.4
4th pair	2	4	2	0, 4.5	36, 126, 173.5, 176	0, 5.0

.

As shown in Table 2, each pair of layover GCPs is correctly identified as being two targets within the slant range-azimuth resolution cell after calibration. After calibration, the height of each target is also accurately estimated. Prior to calibration, however, the situation is disappointing, since not only the height estimation but even the estimate of the target number is wrong.

A precise performance assessment of the proposed algorithm has been carried out by means of Monte Carlo simulations. The systems parameters are the same as those in Table 1. The true APC positions, amplitude inconsistency, and phase inconsistency are modeled as follows:

$$\begin{array}{l}\mathrm{Amplitude}\mathrm{inconsistency}\left(\mathrm{dB}\right)~N\left(0,1\right)\\ \mathrm{Phase}\mathrm{inconsistency}\left(\mathrm{rad}\right)~U\left(-0.5,0.5\right)\\ \mathrm{APC}\mathrm{position}\left(\mathrm{m}\right):\mathbf{x}=[0,1,2,3,4,5,6,7]\times 0.6/7+\Delta \mathbf{x},\\ \mathrm{where}\Delta \mathbf{x}=\left[\Delta x_1,\Delta x_2,...,\Delta x_n,...,\Delta x_8\right]\mathrm{and}\Delta x_n~N\left(0,5\times {10}^{-3}\right)\\ \mathbf{z}=\left[z_1,z_2,...,z_n,...,z_8\right],\mathrm{where}{z}_n~N\left(0,10\times {10}^{-3}\right)\end{array}.$$

(24)

One hundred trials were carried out. For each trial, the mean and standard deviation of the calibration errors (amplitude and phase) are calculated. Then, results of the different trials are averaged. The mean and standard deviation of calibration errors are shown in

Table 3. Calibration errors.

	Mean $\mathit{\mu }$	Standard Deviation $\mathit{\sigma }$
Amplitude (dB)	−35.10	4.15
Phase (rad)	−0.0054	0.0577

. The calibration error of amplitude is below −30 dB in nearly all trials. As for the calibration error of the phase, the standard deviation is less than 0.06 rad. Regarding the calibration error of APC position, it is computed via averaging the RMSE of all trials, where $RMS{E}_i=\sqrt{\left[{\displaystyle \sum _{n=1}^8{\left(\Delta x_{n,i}\right)}^2+{\left(\Delta z_n{,}_i\right)}^2}\right]/8}$, and $\Delta x_{n,i}$ and $\Delta z_{n,i}$ represent the calibration errors of the APC position in the cross-track and height directions, respectively, in the ith trial. The simulation result shows that the calibration error of the APC position is less than 0.127 mm. Based on these simulation results, it can be concluded that the proposed algorithm performs well.

4. Experimental Results

4.1. Data Acquisition

In this section, we will validate the effectiveness of the calibration algorithm by using a real data set acquired by the array InSAR system. The array InSAR system was developed by IECAS in 2014, usually installed onboard a Y-12 aircraft. The radar system operates at 15 GHz (Ku-band), and has eight channels in the cross-track direction. Other system parameters are the same as those in Table 1. As a multi-baseline system, the array InSAR system can obtain 3-D images of the observed scene in a single-pass platform. The Rayleigh resolution [5] in the height direction of the array InSAR system can be derived approximately as

$${\rho }_h=\frac{\lambda R\sin \left(\theta \right)}{B_n}\simeq 35\mathrm{m}$$

(25)

where $\lambda$ represents the wavelength, the slant range $R$ is 1625 m, and the off-nadir angle $\theta$ is 38°. Super-resolution techniques must be used since many targets are under 35 m in height; if not, TomoSAR would be meaningless, especially in urban monitoring and mountain mapping tasks. The utilization of a super-resolution technique can achieve high resolution in the height direction, but it is based on the premise that the channel imbalance and APC position are precisely known. Therefore, calibration is particularly important.

It is clear in Figure 6 that, under the system parameters described above, the quadratic wave model is more accurate than the plane wave model, and that is why we use the Fresnel approximation in our signal model.

The real SAR data set used in this section was acquired by the array InSAR system in April 2015. The calibration site is located in Yuncheng county in the province of Shanxi; its SLC image of Channel 1 is shown in Figure 7. In this single-pass airborne TomoSAR campaign, images of nine GCPs and layover scenes, such as urban buildings, have been acquired. This provides valuable data for our subsequent calibration and validation experiments.

4.2. Validation Data Description

Three special scenes, including a flat ground, non-layover GCP, and layover building, were selected to verify the effectiveness of the unified calibration algorithm. The chosen flat ground is located at the calibration site with an altitude of about 550 m. The 3-D image of this flat ground should be like a horizontal plane. The selected GCPs are trihedral corner reflectors with a leg length of 20 cm, which are placed at a non-layover area; there should be only one target in the height slice for each non-layover GCP. As for the chosen building, serious layover phenomena appear due to the SAR side view principle. Note that there are a maximum of three layover points in a slant range-azimuth resolution cell of the chosen building. Figure 8 shows SLC images of the chosen scenes and corresponding optical images.

4.3. Calibration Algorithm Validation

The calibration processing of the channel imbalance and APC position was carried out as discussed in Section 3. Nine GCPs were used in the calibration site. The results shown in Figure 9 are the channel imbalance and APC position calibration results (blue line) of the array InSAR system campaign from April 2015. The nominal values of the APC position are depicted by a red line. Figure 9a,b show the nominal values and the calibration values of the APC position. The difference between the nominal values and calibration values may be due to the fact that the nominal values are obtained by rough measurement, which is not guaranteed to perfectly correspond to the electromagnetic phase center of the antenna. Figure 9c shows the amplitude difference between different channels. Channel 2 appears slightly smaller than the other channels. We found that the estimate amplitude consistency of each channel is similar to the amplitude calibration results in Figure 9c when the power of the same GCP is extracted in different channels. As the true phase inconsistency and APC position cannot be exactly known, special scenes were chosen to verify the effectiveness of the calibration results in Figure 9. As we shall see later, the 3D imaging results in Figure 10, Figure 11 and Figure 12 show the correctness of the calibration results.

In order to verify whether the calibration values are correct, TomoSAR focusing experiments were carried out for two non-layover GCPs, for flat ground and for a layover building; they are described above. In these experiments, a compressive sensing ()-based TomoSAR spectral estimation algorithm—SL1MMER [4]—was used in the imaging of the height direction. Note that SL1MMER can only achieve the best super-resolution power if the parameters are precisely known. Therefore, the more accurate the parameter calibration is, the better the imaging and resolving power in the height direction will be.

● GCP Results

Figure 10 shows the TomoSAR focusing results in the height slice of GCPs. The APC calibration only means that we only compute the estimation values of the APC position by means of a phase center double localization (PCDL) algorithm [11]. The channel imbalance calibration only means that we only compute the estimation results of the channel imbalance by means of the algorithm proposed in [18]. It is apparent that the TomoSAR focusing results of the unified calibration are more accurate in the target number and position estimation and have less disturbance points than the results of the APC calibration only or the results of the channel imbalance calibration only. The heights of the two chosen GCPs are 550.93 m and 552.79 m, respectively. The reconstructed heights obtained by utilizing the unified algorithm proposed in this paper are 551.06 m and 552.75 m, which is very close to the true heights. However, the results of the other methods are not satisfactory. This is not difficult to understand, as the phase error caused by the channel imbalance and APC position error would result in disturbance points and the target position offset.

● Flat Ground Validation

Figure 11 shows the 3D imaging results of the selected flat ground. The 3D imaging results of the unified calibration algorithm (see Figure 11c,f) are almost clustered on a plane, but the results of the APC calibration only and channel imbalance only are scattered, which is inconsistent with the actual flat terrain. It can therefore be said that the 3D imaging quality is improved by utilizing the system parameters calibrated by the proposed unified algorithm.

● Layover Building Validation

In this partition, the 3D imaging results of the layover building are shown in Figure 12. Figure 12a shows the building observation geometry. In this situation, Points A, B, and C with the slant range r1 will be superimposed within the same slant range-azimuth resolution cell in the 2D SAR image. The same is true for Points D and E. If the channel imbalance and APC position are calibrated precisely, these layover points will be perfectly reconstructed in the right height position utilizing a super-resolution algorithm. On the contrary, these layover points will not be properly reconstructed, such as the occurrence of disturbance points, location shifts, or even defocusing.

When the channel imbalance and APC position have been calibrated by the proposed algorithm, the 3D imaging result is shown in Figure 12b. In this figure, we see that Points A, B, and C are resolved and reconstructed in the corresponding position. This is similar to the other points (i.e., Points D and E) of this layover building. For comparison, the 3D imaging results of the APC calibration only and the channel imbalance calibration only are shown in Figure 12c,d, respectively. Compared with Figure 12c,d, there are fewer disturbance points in Figure 12b, especially in the area near the bottom of the building. This is because the imaging results of these areas are more sensitive to system parameters due to the very small difference in height between the building and the ground. It is therefore safe to say that the imaging results of Figure 12b is better and shows the validity of the unified calibration algorithm proposed in this paper.

To qualitatively exhibit the validity of the calibration algorithm we propose in this paper, the 3D focusing result of a residential area after channel imbalance and APC position calibration is shown below. The SAR 3D focusing result can be approached as follows. Applying a classical imaging algorithm to the raw SAR data collected by each channel, we obtain N 2D SLC images. After some sub-pixel accuracy coregistration to the reference channel and a unified calibration for the channel imbalance and APC position, a super-resolution algorithm, such as compressed sensing [3] or the direction of arrivals [12], is applied to each slant range-azimuth resolution cell. Following this, the 3D imaging result is obtained. The Google Earth image, SLC image, and 3D image of the residential area are shown in Figure 13. As shown in the SLC image, in the 2D SAR image, we cannot even correctly estimate the number of buildings owing to the serious layover phenomenon. When 3D imaging was carried out after precise calibration of the array parameters, the building footprints and the texture of the roof and façade were all clearly visible.

5. Discussion

Both the simulation and array InSAR system results demonstrated the effectiveness of the proposed unified calibration algorithm. The simulation results showed that the calibration accuracy of the phase inconsistency and APC position went up to 0.06 rad and 0.105 mm, respectively. The application to the array InSAR system, a single-pass multi-baseline TomoSAR system, demonstrated the effectiveness of the unified calibration algorithm. Both the point-like targets (e.g., GCPs) and distributed target (e.g., flat ground) presented more reliable 3D imaging results after calibration. In the simulation, the calibration accuracy mentioned above provided the super-resolution ability to distinguish two targets with a 4.5 m height interval (e.g., the 4th pair layover GCPs in Table 2) despite the system’s Rayleigh resolution of only 35 m. The resolution results in the height direction of the array InSAR system showed a similar super-resolution performance with a better separation of the building façade and ground (e.g., Figure 12b).

As discussed before, our proposed algorithm calibrated not only the channel imbalance but also the APC position. This is of great importance since both the channel imbalance (especially the phase inconsistency) and the APC position are the key parameters of the super-resolution algorithm. This proposed algorithm makes it possible to obtain a super-resolution performance in the height direction where the Rayleigh resolution is usually about 10–50 times below that in range or azimuth. In addition, the signal model uses the Fresnel approximation, which increases the accuracy of the calibration signal model. Although the proposed algorithm is of great accuracy, other tasks lie ahead: First, the question of how to choose GCPs for calibration, including the choice of size and position, needs to be studied, since the quality of GCPs will directly affect calibration results. In addition, the derivation of Cramer-Rao lower bound (CRLB) [23,24] and the quantitative assessment method require further study.

6. Conclusions

The multi-baseline TomoSAR system-like array InSAR system, whose antennas are placed in the cross-track direction, can obtain 3-D reconstruction images of the observation area in a single-pass flight campaign. In particular, when using the super-resolution algorithm in height focusing, we can obtain a high height resolution close to the azimuth and ground range resolutions. However, the high resolution in the height direction is based on the condition that the channel imbalance and APC position are precisely known. This is almost impossible, even though high-precision laboratory measurements have been carried out, since there are many factors that may lead to parameters changing in a flight campaign, such as installation error and the difference between the mechanical position and electromagnetic phase center. When a system error occurs, the 3-D resolution performance will be seriously degraded. In order to ensure the 3-D resolution performance of the TomoSAR system, especially the super-resolution performance in the height direction, the calibration processing of the channel imbalance and APC position is therefore necessary.

In our signal model, the Fresnel approximation is used. It is of great significance since the quadratic wave model is more accurate than the plane wave model, and most small aircraft platforms and UAVs work in relatively low flight height conditions. The importance of the proposed algorithm also lies in the fact that both the channel imbalance and APC position are accurately calibrated, which provides the basis for super-resolution in the height direction. This conclusion is supported by experiments based on a simulation data set and real data set acquired by the array InSAR system.