A Novel Phase Compensation Method for Urban 3D Reconstruction Using SAR Tomography

Abstract

Synthetic aperture radar (SAR) tomography (TomoSAR) has been widely used in the three-dimensional (3D) reconstruction of urban areas using the multi-baseline (MB) SAR data. For urban scenarios, the MB SAR data are often acquired by repeat-pass using the spaceborne SAR system. Such a data stack generally has long time baselines, which result in different atmospheric disturbances of the data acquired by different tracks. These factors can lead to the presence of phase errors (PEs). PEs are multiplicative noise for observation data, which can cause diffusion and defocus in TomoSAR imaging and seriously affect the extraction of target 3D information. In this paper, we combine the methods of the block-building network (BBN) and phase gradient autofocus (PGA) to propose a novel phase compensation method called BBN-PGA. The BBN-PGA method can effectively and efficiently compensate for PEs of the MB SAR data over a wide area and improve the accuracy of 3D reconstruction of urban areas. The applicability of this proposed BBN-PGA method is proved by using simulated data and the spaceborne MB SAR data acquired by the TerraSAR-X satellite over an area in Barcelona, Spain.

1. Introduction

With the development of the economy and society, the construction and management of the urban environment are particularly important. The acquisition of spatial information for urban environments, such as the 3D structure of buildings, plays an important role in urban planning, infrastructure construction, disaster monitoring, etc. [1]. Synthetic aperture radar (SAR) tomography (TomoSAR) technology is a three-dimensional imaging technology developed from SAR technology. After years of development, TomoSAR technology has been widely used in 3D imaging in urban and forest areas [2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19].

TomoSAR imaging depends on the phase of multi-baseline (MB) SAR data, which is related to the range between scatterers and sensor [20]. In the actual situation, the MB data can be obtained by airborne or spaceborne SAR systems. The phase errors (PEs) are unavoidable due to the baseline measurement errors or the atmospheric delays [21,22,23]. For airborne MB data, time decorrelation is weak and the effect of atmospheric delay is small, but due to the instability of aircraft movement, the baseline measurement error is relatively large [21,22]. For the spaceborne MB SAR data, the satellite motion is relatively stable and the baseline measurement error is small, but the atmospheric delay effect is relatively large due to the long time baseline [20]. In the case of urban areas, the MB SAR data are often acquired by spaceborne SAR systems. Therefore, PEs are mainly caused by atmospheric delays. PEs can cause high sidelobes or defocus for TomoSAR imaging, which will damage the extraction of target information [22]. Therefore, the PE compensation operation must be performed on the MB SAR data before TomoSAR imaging.

Many methods have been proposed to achieve PE compensation of MB SAR data. These methods can be classified according to the data acquisition platform and can be roughly divided into two categories: the PE compensation methods of spaceborne and airborne MB SAR data.

For spaceborne MB SAR data, PEs are mainly caused by atmospheric delays, and the main application scenarios of the data are urban and mountainous areas. Many methods have been studied to compensate for the atmospheric phase. Classical methods proposed in the earlier literature include persistent scatterer interferometric SAR (PS-InSAR) [24], Small Baseline Subset (SBAS) [25] and other techniques, which are developed on the differential interferometric SAR technique. These methods can obtain estimates of the atmospheric phase. Therefore, these methods can be used for PE compensation of TomoSAR imaging. However, the processing process of these methods is cumbersome, and the atmospheric phase estimation is only obtained by simple space-time filtering, which is less accurate [26]. In [26], the author proposed to use the phase gradient autofocus (PGA) method to estimate and compensate for atmospheric phases. However, the PGA method considers the atmospheric phase to be spatially consistent and can only be applied to small areas. In [11], a two-tier network construction (TTNC) method was proposed to compensate the atmospheric phase to achieve 3D reconstruction of urban buildings. However, this method involves phase unwrapping and large matrix operations, and the computational complexity is very high. Aiming at the drastically changing atmosphere in mountainous areas, the literature [27] proposed a data-driven method to estimate the change in atmospheric delay phase using a regression-kriging framework and used 32-strip SAR images to prove the availability of this method. However, this method requires the use of the PS-InSAR method for initial atmospheric phase estimation. Therefore, the processing process of this method is heavy and complicated, and the performance of this method depends on the accuracy of the initial atmospheric phase estimation.

For airborne MB SAR data, PEs are mainly caused by baseline errors, and the main application scenarios of the data is forest areas. The literature [28] introduced a method to minimize the entropy of the tomographic profile to compensate for PEs. However, this method will introduce unwanted shifts in the elevation direction. To solve the unwanted shifts, the literature [20] proposed a modified method based on the constrained optimization of phase derivatives. In [23], the authors proposed a two-step calibration method to compensate for space-varying PEs in range and azimuth directions. In [29], the authors proposed a phase center double localization method to compensate for PEs without any prior information. However, this method is not suitable for processing the MB SAR data with low coherence images. In [30,31], the PGA method was used to compensate for PEs caused by baseline errors. In [32], the authors combined the network construction (NC) method with the PGA method to propose the NC-PGA method, which can realize the phase error compensation in a large area.

The PGA method was originally used in SAR imaging and has strong robustness in phase error compensation [33]. Moreover, it has been used with TomoSAR technology to compensate for PEs of the spaceborne and airborne MB SAR data [26,30,31]. However, since the PEs caused by atmospheric delays or baseline errors are space-variant, the PGA is not suitable for handling large areas. In [32], the NC-PGA method was proposed to overcome the shortcoming of the PGA method, which can realize the estimation and compensation of spatially variable PEs in a large area. The literature [32] only verified the performance of NC-PGA for PE compensation of the airborne MB SAR data in the forest area application. Since both the network construction (NC) method and the PGA method can be used for the spaceborne MB SAR data, it is not difficult to conclude that the NC-PGA method can also be used for PE compensation of the spaceborne MB SAR data. However, the density of permanent scatterers (PSs) in urban areas is large, and the matrix dimension required for processing by the NC method is large, so the computational complexity is very high. Therefore, for the 3D reconstruction of urban areas, the NC-PGA method is not suitable for large-scale area processing. To address these problems, this paper proposes a BBN-PGA method to compensate for PEs in urban area tomographic imaging. The BBN-PGA method adopted the block-building network (BBN) method to replace the NC method to estimate the elevations of PSs in urban areas and combine the PGA method to compensate for PEs. The main contributions are as follows.

• The BBN-PGA method can obtain estimated elevations of PSs in the whole test area and bring these elevations into the subsequent PGA processing to compensate for PEs of the whole area;
• The BBN method divides the entire test area into several sub-blocks for individual network construction to form multiple permuted sub-networks, and each adjacent sub-network overlaps. Compared with the NC method, this method can reduce the matrix dimension and computational complexity of the algorithm while ensuring the estimation accuracy, greatly improving the processing speed of the algorithm;
• In this paper, the applicability and effectiveness of the NC-PGA method are verified for the first time by using the spaceborne MB SAR data in urban areas. It is also proved that the BBN-PGA method proposed in this paper, as an improved method of NC-PGA, can be used for the PE compensation of airborne and spaceborne MB SAR data and can be applied to the 3D imaging of forests and urban areas.

The rest of this paper is formed as follows. Section 2 gives the tomographic SAR model with PEs and introduces the workflow of the BBN-PGA method. Section 3 gives the experimental results of the proposed BBN-PGA method using simulated data and the TerraSAR-X data stack. A further discussion of the BBN-PGA method is given in Section 4. Finally, conclusions are drawn in Section 5.

2. Methods

2.1. TomoSAR Signal Model with Phase Errors

Figure 1 shows the geometry of the spaceborne TomoSAR imaging system. Let us assume that a SAR satellite acquired N SAR images in the same test area after N repeat flights. One SAR image can be selected as the master image, and the other SAR images are slave images. We can calculate the orthogonal baseline information of N images by using the master image as a reference. As shown in Figure 1, B is the orthogonal baseline span (i.e., tomographic aperture), and $\Delta B$ is the baseline interval for adjacent images. The processing steps of TomoSAR technology mainly include co-registration of images, deramping, PE compensation and tomographic imaging [2]. After the co-registration of images and deramping steps, the TomoSAR signal model can be expressed as [7]:

$$g_n=\underset{\Delta S}{\int }\gamma \left(s\right)exp\left(j2\pi {\xi }_{n}s\right)ds$$

(1)

where $g_{n=1,2\dots ,N}$ is the complex value at a pixel of the image, $\gamma \left(s\right)$ is the reflectivity function along elevation s and $\Delta S$ is the sampling range of s. ${\xi }_n=2b_n/\lambda R$ is the spatial frequency, where $\lambda$ is the wavelength, $b_n$ is the orthogonal baseline of the nth image and R is the slant range between the sensor and target.

For spaceborne MB SAR data, PEs are mainly caused by atmospheric delays and have spatially correlated and temporally uncorrelated distribution characteristics [20,26]. Moreover, the PE is commonly a low-frequency phase distortion, mainly due to spatial inhomogeneity of atmospheric water vapor concentrations [24]. Assuming that ${\phi }_n$ is the PE of the nth SAR image, and considering the presence of noise, Equation (1) can be rewritten as:

$$\mathbf{g}=\Phi \odot \mathbf{A}\Gamma +\mathbf{e}$$

(2)

where $\mathbf{g}={[g_1,\cdots ,g_n]}^{\mathrm{T}}$ is the complex vector, $\Phi ={[exp\left(j{\phi }_1\right),\cdots ,exp\left(j{\phi }_N\right)]}^{\mathrm{T}}$ is the vector of PEs, $\mathbf{A}$ is the $N\times M$ steering matrix with $A_{nm}=exp\left(j2\pi {\xi }_n{s}_m\right)$ (M is the number of discrete samplings of $\Delta S$), $\Gamma ={[\gamma \left(s_1\right),\cdots ,\gamma \left(s_M\right)]}^{\mathrm{T}}$ is the reflectivity vector, $\mathbf{e}$ is the noise vector, ⊙ represents the Hadamard product and ${(·)}^{\mathrm{T}}$ represents transpose operation.

2.2. Phase Error Compensation Using the BBN-PGA Method

This subsection introduces the workflow of the BBN-PGA method. The workflow of the BBN-PGA method is shown in Figure 2. BBN is an improved method of the NC method to estimate elevations of PSs. Then, these elevations can be used in PGA to achieve PE compensation for each subarea.

2.2.1. Applying BBN Method in the Test Area

As an improved method of the NC method, the BBN method reduces the complexity and improves the processing efficiency by dividing the test area into sub-blocks and constructing the network separately in each sub-block. As shown in Figure 2, the BBN method mainly include four steps, as follows:

Firstly, PS selection. The BBN method requires the PSs in the whole area to construct a network. The PSs can be selected based on the amplitude dispersion of the pixel [24], which is defined by:

$$\mu =\sqrt{\frac{1}{N}\sum _{n=1}^N{\left|g_n\right|}^2-{\left(\frac{1}{N}\sum _{n=1}^N\left|g_n\right|\right)}^2}/\phantom{\sqrt{\frac{1}{N}\sum _{n=1}^N{\left|g_n\right|}^2-{\left(\frac{1}{N}\sum _{n=1}^N\left|g_n\right|\right)}^2}\frac{1}{N}\sum _{n=1}^N\left|g_n\right|}\frac{1}{N}\sum _{n=1}^N\left|g_n\right|$$

(3)

A threshold $T_{\mu }$ should be set to judge PSs. The pixel with $\mu <T_{\mu }$ will be selected as a PS.

In the urban environment, the density of PSs is relatively large. Therefore, the number of PSs that can be selected for a large test area is very large, which will make the calculation complexity of the network construction method very large and the calculation time long. To solve this problem, the BBN method proposes the idea of limiting the number of PSs, that is, limiting the maximum number of PSs to p in an area of a certain size. Assuming that there are a total of P PSs in a certain area, the amplitude dispersion value of the p PSs selected in the BBN method must meet the following conditions:

$${\mu }_1\le {\mu }_2\le \cdots \le {\mu }_p\le {\mu }_{p+1}\le \cdots \le {\mu }_{\mathrm{P}}\le T_{\mu }$$

(4)

Secondly, sub-block division. To reduce the computational complexity, the BBN method not only limits the number of PSs but also proposes the idea of sub-block segmentation and overlapping network construction. We can evenly divide the test area into multiple sub-blocks with a certain size $w_a\times w_r$ (for example, 1 km × 1 km), and considering that the elevations of PSs need to be jointly calibrated between each sub-block, adjacent sub-blocks have overlapping areas (the width is $w_o$). As shown in sub-plot b. in Figure 2, the test area is divided into four sub-blocks by red lines, and the gray areas are the overlapping areas of adjacent sub-blocks. Therefore, the size of these sub-blocks is $(w_a+w_o)\times (w_r+w_o)$.

Thirdly, network construction. The Delaunay triangulation network (DTN) is a robust network to estimate the information of PSs [11,24]. Therefore, we can use the PSs in a sub-block area to construct the DTN. The sub-plot c. in Figure 2 shows the result of network construction. The networks of each sub-block are constructed separately and are not connected, but the networks of adjacent sub-blocks have overlapping parts.

Finally, estimate the elevations of PSs. We can use the NC method introduced in [11] and [32] to the DTN of each sub-block to estimate the elevations of PSs in each sub-block. The NC method mainly includes the estimation of relative elevations at arcs of the DTN, the screening of arcs and the estimation of absolute elevations of PSs. For a detailed introduction of these steps, please refer to the literature [11,32], which is not be repeated in this paper. After obtaining the elevations of PSs in each sub-block area, it is necessary to use the elevations of the PSs in the overlapping area of two adjacent sub-blocks to perform elevation calibration, so that the reference elevation of PSs in the entire test area is consistent. We assume that there are two adjacent sub-blocks denoted as $B_i$ and $B_j$. There are I and J PSs in $B_1$ and $B_2$, respectively, and there are p PSs in the overlapping area of $B_1$ and $B_2$. Then, the elevations of PSs in sub-block $B_2$ can be calibrated according to the elevations of PSs in $B_1$ as:

$$s_j=s_j+\frac{1}{p}\sum _{k=1}^p\left(s_k^{B_1}-s_k^{B_2}\right),j=1,\dots ,J$$

(5)

where $s_j$ is the elevation of jth PS in the sub-block $B_2$, and $s_k^{B_1}$ and $s_k^{B_2}$ denote the elevations of the kth PS in the overlapping area in sub-block $B_1$ and sub-block $B_2$, respectively.

2.2.2. Applying PGA for Each Subarea

After obtaining the elevations of PSs using the BBN method, these elevations can be used in the PGA method to remove the linear phase term (caused by the elevation) of the PS and better complete the estimation and compensation of PEs. As shown in Figure 2, the test area should be first divided into multiple subareas. After that, the PE of each subarea can be compensated by the PGA method. The size of subareas can be set according to the actual situation, generally $50\times 50$ m to $100\times 100$ m can be selected.

Assuming that the subarea contains K PSs, the signal of the kth PS in the nth SAR image is:

$$g_{k,n}={\gamma }_{k}exp\left(j{\phi }_n+j2\pi {\xi }_n{s}_k\right)+e_{k,n}$$

(6)

where ${\gamma }_k$ is the complex reflectivity, $s_k$ represents the elevation of the kth PS and $e_{k,n}$ denotes the noise component [30,31]. The BBN-PGA method considers the phase error to be spatially consistent within the subarea. Therefore, the phase error ${\phi }_n$ in Equation (6) does not have changes along azimuth and range directions and is a fixed value within the subarea.

Using the PGA method to estimate PEs can be roughly divided into following steps.

• Remove the linear phase term. The linear phase $2\pi {\xi }_n{s}_k$ (shown in Equation (6)) must be removed, then the gradient of ${\phi }_n$ can be estimated. After this operation, the $g_{k,n}$ can be expressed as:

  $$g_{k,n}^{\prime }={\gamma }_{k}exp\left(j{\phi }_n\right)+e_{k,n}^{\prime }$$

  (7)

  where $e_{k,n}^{\prime }=e_{k,n}·exp(-j2\pi {\xi }_n{s}_k)$ can be processed as white noise.
• Estimate the phase gradient. The phase gradient ${\widehat{\dot{\phi }}}_n$ of ${\phi }_n$ can be estimated using the maximum-likelihood estimation (MLE) kernel, the expression of which is:

  $${\widehat{\dot{\phi }}}_n=arg\left\{\sum _{k=1}^K\left[{\left({g^{\prime }}_{k,n-1}\right)}^{*}·{g^{\prime }}_{k,n}\right]\right\},2\le n\le N$$

  (8)

  where ${(·)}^{*}$ represents the conjugate operation, and $arg(·)$ denotes the operation of phase extraction.
• Estimate phase errors. The PEs ${\phi }_n$ can be acquired by integrating the phase gradient ${\widehat{\dot{\phi }}}_n$, the expression of which is:

  $${\widehat{\phi }}_n=\left\{\begin{array}{cc}0,\hfill & n=1\hfill \\ {\displaystyle \sum _{i=2}^n}{\widehat{\dot{\phi }}}_i,\hfill & n=2,3,\dots ,N\hfill \end{array}\right.$$

  (9)
  Then, the estimated PEs can be removed from the data of the subarea.

To improve the accuracy and robustness of the phase error compensation, the PGA method adopts an iterative method for processing. A threshold $T_{pga}$ can be set to determine the convergence of the iteration, which can be expressed as [33]:

$$\sum _{n=1}^N{\left|{\widehat{\phi }}_n^a-{\widehat{\phi }}_n^{a-1}\right|}^2<T_{pga}$$

(10)

where a is the number of iterations, ${\widehat{\phi }}_n^a$ is the estimated PEs at the ath iteration and the $T_{pga}$ is commonly set to ${10}^{-3}$.

2.3. TomoSAR Imaging

After the PE compensation, we can use the TomoSAR imaging method to extract the height of scatterers by using the calibrated MB SAR data. In this paper, beamforming (BF) is used to realize tomographic imaging. The basic idea of the BF algorithm is to design a spatial filter with the finite impulse response, so that the spatial frequency corresponding to the strong scatterer with a height of s can pass through without distortion, and the response amplitude of other spatial frequencies can be weakened as much as possible [34,35]. The power of spatial spectrum estimated by using BF method as [2,36]:

$${\widehat{P}}_{BF}\left(s\right)=\frac{{\mathbf{a}}^H\left(s\right)\mathbf{Ra}\left(s\right)}{N^2}$$

(11)

where $\mathbf{a}\left(s\right)={\left[exp\left(j2\pi {\xi }_{1}s\right),exp\left(j2\pi {\xi }_{2}s\right),\cdots ,exp\left(j2\pi {\xi }_{N}s\right)\right]}^{\mathrm{T}}$ is the steering vector, ${(·)}^H$ denotes the Hermitian operator and $\widehat{\mathbf{R}}$ is the covariance matrix and can be estimated over L looks as:

$$\widehat{\mathbf{R}}=\frac{1}{L}\sum _{l=1}^L{\mathbf{g}}_l{{\mathbf{g}}_l}^H$$

(12)

The ${\mathbf{g}}_{l=1,\cdots ,L}$ can be selected from the pixels around a fixed pixel. [32,36]. Then, the heights of scatterers can be estimated by detecting the peaks of ${\widehat{P}}_{BF}$.

3. Experiments and Results

In this section, we conduct experiments with simulated and real data to evaluate the performance of the BBN-PGA method and demonstrate its effectiveness and practicality. Furthermore, we conducted a comparative analysis using PGA [30,31], TTNC [11], and NC-PGA [32] methods to prove the advantages of BBN-PGA methods.

3.1. Simulated Experiments and Results

In this experiment, the parameters and baseline configurations of the TerraSAR-X data stack (shown in Section 3.2.1) can be used to construct the simulated data. In the simulated experiment, we simulate a $500\times 500$ pixels area. The characteristics of the simulated scene are shown in

Table 1. The characteristics of simulated scene.

Location	Elevation [m]	Size [pixels]
Top left	80	$150\times 150$
Top right	25	$150\times 150$
Center	5–128	$150\times 100$
Bottom left	100	$150\times 150$
Bottom right	45	$150\times 150$

. The simulated complex data uses the signal model of Equation (6), and the signal-to-noise ratio (SNR) is set to 5 dB. Figure 3 shows the elevations and mean amplitude of the simulated data.

As discussed in Section 2, the PEs caused by atmospheric delay exhibit low-frequency distortion in space and are incoherent in the temporal dimension [24]. Therefore, in this simulation experiment, we consider that PEs exhibit a linear change in the azimuth and range directions. The PE ${\phi }_n$ can be expressed as:

$${\phi }_n(x,r)=c_1·{\alpha }_{1,n}+\frac{c_2·{\alpha }_{2,n}}{N_x}·x+\frac{c_3·{\alpha }_{3,n}}{N_r}·r$$

(13)

where ${\alpha }_{1,n}$, ${\alpha }_{2,n}$ and ${\alpha }_{3,n}$ are uniformly distributed random numbers in $[-0.5,0.5]$; $c_1$, $c_2$ and $c_3$ are constants that control the range of PE fluctuations; x and r denote the azimuth and range coordinates of the pixel, respectively; and $N_x$ and $N_r$ denote the number of pixels in the azimuth and range directions, respectively. In this experiment, the PEs are simulated using Equation (13) by setting $c_1=\pi$ radian and $c_2=c_3=2\pi$ radian.

To compensate for PEs, we first use the BBN method to obtain the elevations of PSs. In the experiment, wet set $T_u=0.23$ to select PSs and restricted it to keep 20 PSs in an area of size $50\times 50$ pixels. Figure 4 shows the simulated result of the BBN method, including selected PSs, the networks of sub-blocks and estimated elevations of PSs. As shown in Figure 4b, the simulated test area is divided into 4 sub-blocks to build 4 independent networks, and adjacent networks have overlapping parts. Comparing Figure 3a and Figure 4c, it can be seen that the estimated elevations are basically the same as true elevations.

After obtaining the elevations of PSs, we can take these estimated elevations into the PGA method and set the size of the subarea to $100\times 100$ pixels for PE estimation and compensation. Then, we use the BF method for TomoSAR imaging and extract the elevation of each pixel. Figure 5 shows the comparison of the elevation estimated results of the simulated test area. Figure 6 gives the joint distribution of the true elevation and the estimated elevation. We also calculate the bias, root mean square error (RMSE), and the coefficient of determination ($R^2$) of the estimated elevations from the true elevations to evaluate the performance of the proposed BBN-PGA method. $R^2$ is a number that can be used to represent the goodness of fit of the estimated value to the reference value [37]. The evaluation results are given in Table 2. From the simulated results, some conclusions can be summarized as follows.

• The existence of PEs will seriously damage the extraction of target information, resulting in a large deviation of the estimated elevation (as shown in Figure 5b and Figure 6b);
• The simulated results of the PGA method are significantly improved. The elevations of most pixels are estimated relatively correctly. However, there are many outliers (with larger errors of the estimated elevation) at the edges of the test area (as shown in Figure 5c and Figure 6c). This is because the PGA method assumes that the PEs have spatial consistency and cannot achieve high-precision compensation of the space-variant PEs;
• As can be seen from Figure 5d–f and Figure 6d–f, the estimated elevations of the three methods TTNC, NC-PGA and BBN-PGA are in good agreement with the true elevations, which proves that the three methods have good PE compensation capability. Combined with the evaluation results given in Table 2, the PE compensation capabilities of the NC-PGA and BBN-PGA methods are better than the TTNC method. This is because the NC-PGA and BBN-PGA methods are better than the TTNC method in terms of bias, RMSE and $R^2$.
• The PE compensation capability of methods NC-PGA and BBN-PGA is almost the same. The elevation estimation accuracy after compensating the PEs by the two methods is almost the same as that without PEs. This also proves that the two methods, NC-PGA and BBN-PGA, can well compensate the space-variant PEs.

Figure 5. Comparison of elevation estimated results of the simulated test area: (a) PEs free data; (b) the data with PEs; the data after PE compensation by using the (c) PGA, (d) TTNC, (e) NC-PGA and (f) BBN-PGA methods.

Figure 5. Comparison of elevation estimated results of the simulated test area: (a) PEs free data; (b) the data with PEs; the data after PE compensation by using the (c) PGA, (d) TTNC, (e) NC-PGA and (f) BBN-PGA methods.

Figure 6. Joint distribution of the true and estimated elevations: (a) PEs free data; (b) the data with PEs; the data after PE compensation by using the (c) PGA, (d) TTNC, (e) NC-PGA and (f) BBN-PGA methods.

Figure 6. Joint distribution of the true and estimated elevations: (a) PEs free data; (b) the data with PEs; the data after PE compensation by using the (c) PGA, (d) TTNC, (e) NC-PGA and (f) BBN-PGA methods.

Table 2. Evaluation results of simulated experiments.

Method	Bias (m)	RMSE (m)	${\mathbf{R}}^2$
PEs free	−0.002	2.099	0.9961
With PEs	−44.93	65.47	–
PGA	−2.497	20.83	0.6147
TTNC	−2.597	5.889	0.9689
NC-PGA	−0.017	2.141	0.9959
BBN-PGA	−0.056	2.161	0.9959

Table 2. Evaluation results of simulated experiments.

Method	Bias (m)	RMSE (m)	${\mathbf{R}}^2$
PEs free	−0.002	2.099	0.9961
With PEs	−44.93	65.47	–
PGA	−2.497	20.83	0.6147
TTNC	−2.597	5.889	0.9689
NC-PGA	−0.017	2.141	0.9959
BBN-PGA	−0.056	2.161	0.9959

Furthermore, we compare the computational complexity of different methods by analyzing their running time.

Table 3. Running time of different PE compensation methods.

Method	Quantity (s)
PGA	9.1
TTNC	120.08
NC-PGA	87.61
BBN-PGA	5.92

gives the running time of different PE compensation methods for processing the entire test area (a desktop with an Intel Xeon Gold 5218 16-core processor). The BBN-PGA method takes the shortest time, and the PE compensation process is completed in only 5.92 s, which is about 14 times faster than the NC-PGA method and nearly 20 times faster than the TTNC method.

The above simulated experimental results demonstrate that the BBN-PGA method can realize the compensation of space-variant PEs, and on the basis of ensuring that the performance is similar to the NC-PGA method, the computational complexity is reduced and the computational efficiency is greatly improved.

3.2. Real Data Experiments and Results

3.2.1. Study Area and Data Stack

The MB SAR data is acquired by the TerraSAR-X satellite. It has 24 stripmap SAR images acquired from July 2008 to November 2009 over Barcelona, Spain. Figure 7 shows the study area in this experiments. The size of the study area is $1000\times 2000$ pixels and its spatial dimensions are approximately $1.8\times 2.4$ km.

Table 4. The parameters of the TerraSAR-X data.

Parameters	Quantity
Orbit height	513 km
Wavelength	3.11 cm
Center incidence angle	$35.{32}^{\circ }$
Center slant range	618 km

gives the basic parameters of the TerraSAR-X data. Figure 8 shows the spatial and temporal baselines of the MB SAR data. The span of the spatial baselines is about 246.7 m, and the span of the temporal baseline is 374 days.

3.2.2. Inversion of Urban Height

In this experiment, we set the sub-block size to $500\times 500$ pixels and the overlapping area width to 50 pixels for the BBN-PGA method to construct the sub-block network to compensate for PEs of the real data stack. After compensating for PEs, we use the BF method given by Equation (11) for tomographic imaging and extract the target height by detecting the peak. Figure 9 shows the estimated heights of the study area. The results can be summarized as follows.

• Before PE compensation, the height estimation results of urban buildings are very poor, and there are many height outliers.
• After using PGA to compensate for PEs, it can be seen that the height estimation results in the central area of the study area have been improved, and the heights of buildings are reasonable. However, there are still many outliers (as shown in Figure 9b). Additionally, the estimated heights of the buildings near the edge area are poor (marked by the red ellipse in Figure 9b). This also proves that the PGA method is not suitable for large area processing.
• After using BBN-PGA to compensate for PEs, the distribution of the estimated building heights in the study area is reasonable, and there are very few outliers. This proves that the proposed BBN-PGA method compensates well for PEs of the real data over a large area.
  Due to the high computational complexity of the TTNC and NC-PGA methods, the computation for a larger area takes a long time and occupies a large amount of memory. Therefore, in this experiment, a test area with a size of $1000\times 1000$ pixels (as shown in the red box in Figure 7b) is selected to compare the performance of the TTNC, NC-PGA and BBN-PGA methods. Figure 10 shows the estimated heights of the test area after PE compensation. The estimated heights of the three methods are similar, which can reflect that the PE compensation capabilities of the three methods are similar. To evaluate the performance of the three methods more quantitatively, we chose the high-rise building 1 (as shown in Figure 10a) for further comparative analysis.

Figure 9. The estimated heights of the study area: (a) without PE compensation; with PE compensation using the (b) PGA and (c) BBN-PGA methods.

Figure 9. The estimated heights of the study area: (a) without PE compensation; with PE compensation using the (b) PGA and (c) BBN-PGA methods.

A height model can be simulated to evaluate the accuracy of height estimation. The height model contains the true height of scatterers covering building 1. The model height increases linearly from the bottom to the top of building 1. Figure 11 shows the joint distribution of the model height and the estimated height. It can be seen that the estimated heights of most scatterers agree well with the model height. However, there are some scatteters in the red elliptical area in Figure 11 that do not match the linear model height, which is caused by the overlapping of the roof and the wall of building 1. We calculate the bias, RMSE, and $R^2$ of the estimated heights from the model heights. It should be noted that the overlapped region no longer satisfies the linear height model (which matches the actual situation). Therefore, we remove the overlapped part shown in Figure 11 when calculating the evaluation indicators to ensure the accuracy of the evaluation results. We also use the same computer (a desktop with an Intel Xeon Gold 5218 16-core processor) to count the running time of the three methods, TTNC, NC-PGA and BBN-PGA, to complete the PE compensation in the test area.

Table 5. Evaluation results of real data experiments.

Method	Bias (m)	RMSE (m)	${\mathbf{R}}^2$	Running Time (s)
TTNC	−1.34	7.84	0.9111	3191.3
NC-PGA	0.26	4.35	0.9416	3081.8
BBN-PGA	0.31	3.64	0.9510	27.6

shows the elevation results of this experiments. It can be seen that in the term of PE compensation, the methods NC-PGA and BBN-PGA have similar performance and are superior to the TTNC method. This is because the evaluation indicators of NC-PGA and BBN-PGA are less different and better than the TTNC method. In terms of processing speed, the BBN-PGA method takes the shortest running time, and its processing speed is more than 110 times faster than the NC-PGA and TTNC methods.

The above simulated and real data experimental results demonstrate that the BBN-PGA method can realize the compensation of space-variant PEs over urban areas, and its performance is better than that of the PGA and TTNC methods. Moreover, the proposed BBN-PGA method can reduce computational complexity and greatly improve computational efficiency while keeping the performance similar to the NC-PGA method.

4. Discussion

The computational complexity and generality of the BBN-PGA and NC-PGA methods are discussed in this section.

4.1. Computational Complexity

We contrast the NC-PGA and BBN-PGA methods through an analysis of the most computationally complex step. The NC-PGA method uses the weighted least square (WLS) estimators to estimate elevations of selected PSs, which can be expressed as [38]:

$$\mathbf{X}={\left({\mathbf{D}}^T\mathbf{WD}\right)}^{-1}{\mathbf{D}}^T\mathbf{WS}$$

(14)

where $\mathbf{X}={[x_1,\cdots ,x_P]}^{\mathrm{T}}$ is a $P\times 1$ vector that contains the absolute elevations of selected PSs (P is the number of selected PSs), $\mathbf{S}=[s_1,\cdots ,s_Q]$ is a $Q\times 1$ vector that contains the relative elevations at arcs (Q is the number of the retained arcs in DTN), $\mathbf{W}$ is a $Q\times Q$ is a diagonal matrix with weighted values diagonally (the weighted values can be determined according to the literature [11]) and $\mathbf{D}$ is a $Q\times P$ matrix, that is [32]:

$$\mathbf{D}=\left[\begin{array}{cccccc}1& \cdots & -1& \cdots & 0& \cdots \\ ⋮\\ \cdots & 1& 0& \cdots & -1& \cdots \\ ⋮\end{array}\right]$$

(15)

where $-1$ represents the start PS, and 1 represents the end PS of one arc in the DTN, respectively. Equation (14) involves the matrix multiplication and matrix inversion operations of the large matrices D and W, and its computational complexity is very high. Moreover, the dimensions of matrices D and W are related to the number of PSs and arcs. Therefore, the computational complexity depends on the number of PSs (P) and the number of arcs (Q) in the DTN.

We assume that the NC-PGA method selects $P_1$ PSs in the entire test area, and the DTN has $Q_1$ arcs, then, the computational complexity of the WLS estimator in the NC-PGA method is about:

$$O(P_1{Q}_1+2P_1{Q_1}^2+2Q_1{P_1}^2+{P_1}^3)$$

(16)

For the same test area, it is assumed that the BBN-PGA method selects $P_2$ PSs, and the corresponding number of DTN arcs is $Q_2$. We assume that the BBN method divides the test area into L sub-blocks evenly, and each sub-block contains $P_2/L$ PSs and about $Q_2/L$ arcs, then, the computational complexity of the WLS estimator in the BBN-PGA method is about:

$$O(\frac{P_2{Q}_2}{L}+\frac{2P_2{Q_2}^2+2Q_2{P_2}^2+{P_2}^3}{L^2})$$

(17)

Therefore, the computational complexity of the BBN method is not only related to the number of PSs and arcs but also inversely proportional to the number of sub-blocks. In addition, considering the overlapping area between sub-blocks, the computational complexity of the WLS estimator in the BBN-PGA method is slightly larger than Equation (17). Comparing Equations (16) and (17), when $P_1=P_2$, the computational complexity of BBN-PGA is about L to $L^2$ times lower than that of NC-PGA. Moreover, BBN-PGA also limits the number of PSs, that is, $P_2<P_1$, so the computational complexity of BBN-PGA is much lower than that of NC-PGA.

The analysis of the results of the experiments with simulated and real data is given below to verify the above discussion.

• In the simulated experiments, we set $T_{\mu }=0.23$ to select PSs. The NC-PGA method preserves 9793 PSs and $29,297$ arcs in DTN. While the BBN-PGA method retains $[720,600,600,500]$ PSs and $[2135,1778,1780,1478]$ arcs in the four sub-block DTNs, respectively. Figure 12 show the DTN of NC-PGA and BBN-PGA methods. The running time of the WLS step for NC-PGA and BBN-PGA is 68.69 s and 0.17 s, respectively. In the processing speed of the WLS step, the BBN-PGA method is about 405 times faster than the NC-PGA method, which is consistent with the above analysis. Table 3 shows the overall running time comparison of NC-PGA and BBN-PGA. The results show that BBN-PGA takes 5.92 s and NC-PGA takes 87.61 s. The BBN-PGA method is about 14 times faster than the NC-PGA method.
• In the real data experiments, we set $T_{\mu }=0.23$ to select PSs. The NC-PGA method preserves 32,937 PSs and 98,663 arcs in DTN. While the BBN-PGA method retains $[2327,2115,1882,1913]$ PSs and $[6949,6314,5598,5713]$ arcs in the four sub-block DTNs, respectively. Figure 13 show the DTN of NC-PGA and BBN-PGA methods. The running time of the WLS step for NC-PGA and BBN-PGA is 2971.83 s and 3.99 s, respectively. In the processing speed of the WLS step, the BBN-PGA method is about 744 times faster than the NC-PGA method, which is consistent with the above analysis. Table 5 shows the overall running time comparison of NC-PGA and BBN-PGA. The results show that BBN-PGA takes 27.6 s and NC-PGA takes 3081.8 s. The BBN-PGA method is about 110 times faster than the NC-PGA method.

To sum up, the BBN method reduces the number of PSs and arcs involved in the operation of Equation (14) by limiting the number of PSs and dividing the sub-blocks, which can significantly reduce the computational complexity and improve the processing efficiency.

4.2. Generality

In this paper, we demonstrate that the NC-PGA method can achieve the estimation and compensation of PEs caused by atmospheric delays. Furthermore, in [32], the authors have proved the good ability of the NC-PGA method in compensating for PEs caused by baseline errors. Therefore, the NC-PGA method has the capability of estimating and compensating the PEs of the spaceborne and airborne MB SAR data stacks and is suitable for TomoSAR processing in urban and forest areas. As an improved method of NC-PGA, in this paper, the proposed BBN-PGA method was proven to have a performance close to that of NC-PGA in the PE compensation of spaceborne MB SAR data acquired over urban areas. We have reason to believe that the BBN-PGA method also has PE estimation and compensation ability for airborne MB SAR data and can be applied to forest areas.

Here, we conduct experiments with airborne MB SAR data to verify the generality of the BBN-PGA method. The MB SAR data was acquired by the ONERA SETHI airborne system in the frame of the European Space Agency’s campaign TropiSAR [39]. Figure 14 shows the SAR image of the test forest area. The test area is located in French Guiana, and its introduction can be found in [20,32]. Figure 15 shows tomographic imaging results along a line (marked by the red line in Figure 14). Obviously, the TomoSAR profiles of the data after PE compensation using the BBN-PGA method are well focused, and the sidelobes of profiles are well suppressed, which are highly similar to those of NC-PGA.

In summary, the BBN-PGA method and the NC-PGA method have similar PE compensation performance. Both methods can achieve good compensation for PEs of spaceborne and airborne MB SAR data stacks and are suitable for tomographic processing in urban and forest areas. Therefore, both BBN-PGA and NC-PGA methods have good practicability and generality.

5. Conclusions

To realize the 3D reconstruction of urban areas, a novel BBN-PGA method is proposed in this paper. The BBN-PGA method combines the BBN method and the PGA method to compensate for PEs in a large urban area using spaceborne MB SAR data. The BBN method can estimate the accurate elevation of PSs in the whole area, which can be used in PGA to complete PE compensation. Furthermore, the BBN method reduces the computational complexity by limiting the number of PSs and dividing the sub-blocks, which greatly improves the processing efficiency of the method.

This paper has demonstrated the effectiveness and applicability of the proposed BBN-PGA method by using simulated data and the TerraSAR-X data stack. The experimental results show that the BBN-PGA method can well realize PE compensation of the spaceborne MB SAR data, so that the TomoSAR technology can better and more accurately extract the height of urban buildings. Compared with the PGA and TTNC methods, the PE compensation performance of the BBN-PGA method is better, and it can achieve more accurate estimated heights of scatterers. Compared with the NC-PGA method, the performance of the BBN-PGA method is similar to that of the NC-PGA method, but the computational complexity of the BBN-PGA method is lower, and its processing efficiency is greatly improved (as shown in Table 3 and Table 5, the processing speed of BBN-PGA is about 14 times and 110 times faster than that of NC-PGA). Furthermore, this paper is the first time the capability of the NC-PGA method in PE compensation of spaceborne MB SAR data has been verified. We also demonstrate that the proposed BBN-PGA method can achieve PE (caused by atmospheric delays or baseline errors) compensation for both spaceborne and airborne MB SAR data for tomographic processing in urban and forested areas. These experiments are sufficient to demonstrate the effectiveness, practicability and generality of the proposed BBN-PGA method in PE compensation.

Of course, the analysis of the BBN-PGA method in this paper is not comprehensive enough. In the future, we will focus on the impact of the size of sub-blocks and the number of PSs on the performance of the BBN-PGA method to improve the robustness of this method, so that the BBN-PGA method can be applied to more scenes such as mountains and glaciers.