Elevation Extraction from Spaceborne SAR Tomography Using Multi-Baseline COSMO-SkyMed SAR Data

Abstract

SAR tomography (TomoSAR) extends SAR interferometry (InSAR) to image a complex 3D scene with multiple scatterers within the same SAR cell. The phase calibration method and the super-resolution reconstruction method play a crucial role in 3D TomoSAR imaging from multi-baseline SAR stacks, and they both influence the accuracy of the 3D SAR tomographic imaging results. This paper presents a systematic processing method for 3D SAR tomography imaging. Moreover, with the newly released TanDEM-X 12 m DEM, this study proposes a new phase calibration method based on SAR InSAR and DEM error estimation with the super-resolution reconstruction compressive sensing (CS) method for 3D TomoSAR imaging using COSMO-SkyMed Spaceborne SAR data. The test, fieldwork, and results validation were executed at Zipingpu Dam, Dujiangyan, Sichuan, China. After processing, the 1 m resolution TomoSAR elevation extraction results were obtained. Against the terrestrial Lidar ‘truth’ data, the elevation results were shown to have an accuracy of 0.25 ± 1.04 m and a RMSE of 1.07 m in the dam area. The results and their subsequent validation demonstrate that the X band data using the CS method are not suitable for forest structure reconstruction, but are fit for purpose for the elevation extraction of manufactured facilities including buildings in the urban area.

1. Introduction

3D SAR tomography (TomoSAR) [1,2,3,4] and its closely associated 4D SAR differential tomography (Diff-TomoSAR) [5,6,7,8,9] take advantage of multi-baseline SAR data stacks to create an important innovation in SAR interferometry, allowing for the sensing of complex scenes with multiple scatterers mapped into the same SAR cell. The idea of tomographic imaging was first introduced to the field of SAR research in the 1990s [10,11,12] in order to overcome the limitations (e.g., 3D information extraction) of 2D SAR imaging. The initial experiment was carried out in a laboratory under ideal experimental conditions [13] and subsequently by using airborne systems [1]. In 2010, researchers such as Zhu and Bamler of the DLR Laboratory in Germany used a tomographic SAR inversion method based on L1 norm compression sensing to separate scattering particles distributed along the vertical dimension within the same cell using spaceborne SAR data. Subsequently, this method has been applied to 3D ultra-high-resolution tomographic SAR imaging in complex urban environments, while the compressive sensing method based on the L1 norm has also been applied to differential tomography SAR [14,15,16]. In this way, classical 2D InSAR can be considered as a simple parametric case of 3D TomoSAR. D-TomoSAR (4D SAR imaging) [6,16,17] exploits the strengths of both TomoSAR and PSI, which inverts the motion of the scatterers with 3D TomoSAR reconstruction [16], which can retrieve the motion [18,19] and elevation information of multiple scatterers in a SAR pixel cell using a spectral analysis method [20]. In addition to 3D shape reconstruction, they permit a solution for monitoring deformation in complex urban/infrastructure areas [2,4] as well as recent cryospheric ice investigations [21], emerging tomographic remote sensing applications including forest scenarios [3,22,23] (e.g., tree height and biomass estimation), sub-canopy topographic mapping, and even search, rescue, and surveillance applications under forest.

However, these scenes are characterized and influenced by DEM uncertainty, the temporal decorrelation of scatterers, orbital, tropospheric and ionospheric phase distortion, and an open issue regarding possible height blurring and accuracy losses for TomoSAR applications, particularly in densely vegetated mountainous rural areas and polar icy regions of the Earth, alongside the polar region of many planets (the Moon, Mars, and Mercury) and their icy satellites (moons) such as Europa, Ganymede, Callisto, and so on. To address the problems inherent in these techniques, baseline estimation should be conducted in advance and phase calibration is critical and indispensable in 3D SAR tomographic imaging. Here, COSMO-SkyMed spaceborne SAR data were applied to various applications in the world. Currently, there have only been rare reports on TomoSAR applications using these particular SAR satellite data. This paper presents a systematic processing method for 3D SAR tomography imaging using COSMO-SkyMed spaceborne SAR. In addition, a new phase calibration method is presented based on SAR interferometry (InSAR) and DEM error estimation, and compensation with the newly released geographically extremely limited areas of TanDEM-X 12 m DEM for 3D TomoSAR imaging is also shown. The test, fieldwork, results validation, new method discussion, CS algorithm adaptability analysis, and results analysis were conducted at Zipingpu Dam, Dujiangyan, Sichuan, China. In Section 2.1, the mathematical principles are described, whilst in Section 2.2, the SAR interferometry phase (InSAR) calibration with the DEM error is described. In Section 3, the results are shown for the dam test site whereas in Section 4, these are discussed and our conclusions are presented in the final section.

2. Methodology

2.1. Principles

The two-dimensional imaging principle for SAR assumes that the target is a two-dimensional plane target, while the actual target exists in a three-dimensional space. Two-dimensional SAR imaging results from the projection of the three-dimensional structure and targets into a two-dimensional plane (azimuth-range plane). Three-dimensional space is divided into a series of equidistant cylindrical surfaces along the azimuth axis, but SAR fails to distinguish targets on the same cylindrical surface because these targets are compressed into the same pixel at the same distance. This is the cylindrical symmetric ambiguity problem in SAR two-dimensional imaging. Since high azimuth resolution imaging can be obtained by the synthesis aperture method using a small sized radar antenna, two-dimensional high-resolution imaging in the elevation plane can be realized if a two-dimensional synthetic aperture is formed in the elevation direction with the help of high-resolution range direction imaging to achieve a true three-dimensional radar imaging. This is the basic idea of three-dimensional SAR imaging, TomoSAR. A normal monostatic imaging SAR system consists of a side-looking transmitter and receiver mounted on a moving platform such as an airplane or satellite. TomoSAR builds up a synthetic elevation aperture from a stack of N complex SAR datasets of the same area taken at different times and slightly different orbit positions. The native 3D reference frame of a SAR sensor and the parameters of SAR maps in the three directions are defined as below (x azimuth, r range, s elevation), which is displayed in the TomoSAR imaging geometry in Figure 1. According to previous research [1,2,3,4], the workflow of SAR tomography and D-TomoSAR is shown in Figure 2. First, all complex images were co-registered into SAR stacks. Then, the atmospheric and ionospheric corrections were performed. Next, deramping and phase error compensation were executed. Finally, the results can be obtained after 3D TomoSAR imaging and post-processing.

Theoretically, the range-azimuth resolution cell in a sequence of images ($M$ images) corresponds to the same feature with pixels indexed by ($s$), and its complex image single-look-complex (SLC) value $Q\left(m\right)$ which can be written as below [2,24],

$$Q\left(m\right)={{\displaystyle \int }}_{s_{min}}^{s_{max}} \gamma \left(s\right)exp\left(-j\frac{4\pi f_0}{c}{R}_m\left(s\right)\right)ds, m=1,2,\dots ,M,$$

(1)

where $\gamma \left(s\right)$ refers to the complex scattering coefficient; s represents the elevation; $f_0$ is the frequency; $c$ is the speed of light; and $R_m\left(s\right)$ is the slant range. After removing the phase term caused by the reference slant range $R_m\left(0\right)$ at the phase center of each image, called deramping [2,24], the mth image can be written as $\mathrm{g}\left(m\right)$.

$$\begin{gathered} \mathrm{g}\left(m\right)=exp\left(j\frac{4\pi f_0}{c}{R}_m\left(0\right)\right)\cdot Q\left(m\right), \\ ={{\displaystyle \int }}_{s_{min}}^{s_{max}}\gamma \left(s\right)exp\left(-j\frac{4\pi f_0}{c}\left(R_m\left(s\right)-R_m\left(0\right)\right)\right)ds, \end{gathered}$$

(2)

After calculation and simplification [2,25], Equation (2) be expressed in the following form,

$$\mathrm{g}\left(m\right)={{\displaystyle \int }}_{s_{min}}^{s_{max}}\gamma \left(s\right)exp\left(j2\pi {\xi }_{m}s\right)ds,$$

(3)

$${\xi }_m=\frac{2b_{\perp m}}{\lambda r},$$

(4)

where ${\xi }_m$ is the spatial frequency corresponding to the height $s$ from the normal slant-range direction (NSR direction or here called $s,$ elevation direction), and $b_{\perp m}$ refers to the perpendicular baseline. As seen from Equation (3), the complex values $\mathrm{g}\left(m\right)$, $m=1,2,\cdot \cdot \cdot ,M$ are the target discrete cell samples of electromagnetic scattering characteristics $\gamma \left(s\right)$ along the NSR direction at ${\xi }_m$ after SLC image deramping. In other words, SAR tomographic imaging is essentially the use of spectral discrete sampling reconstruction of the original signal problem, and it is necessary to obtain the spatial frequency ${\xi }_m$ first.

Based on Equation (2), the deramping operation plays a critical role in SAR tomography. It is precisely because of the removal of the central slant phase that the intrinsic relationship between the observed data and the target NSR information is established.

The reference slant range used in the deramping step can be obtained in two ways: (1) calculating the reference slant range using the radar transmission center time delay and velocity recorded by the radar; and (2) calculating the reference slant range based on some reference topography and radar position [1,2,3,4]. The reference slant range used in the deramping can be the center distance of each image, or the distance between the radar antenna phase center of each image and a reference terrain (e.g., known coarse resolution DEM data). Although the most direct deramping method is to remove the slant range to the phase center via the radar-recorded electromagnetic propagation delay, the use of its deramping introduces additional atmospheric phase errors because the echo delay is affected by atmospheric interference, which also needs to use an offset co-registration interpolation method to calculate the reference slant range. A reference terrain is more commonly used because it only requires the calculation of the reference slant range based on the radar position and the terrain [24,26]. Therefore, a coarse DEM is generally applied via flattening of the simulated phase by a DEM. In addition, for scenes with lower relief, one can also use ellipse models that have been improved based on the average elevation of the scene [24,26]. According to previous research [27], the phase simulation by DEM is expressed as:

$${\phi }_{simu\_m}=-\frac{4\pi }{\lambda }{b}_{\parallel m},$$

(5)

where $b_{\parallel m}$ is the parallel baseline; therefore, according to Equation (5), it can be observed that parallel baseline estimation is extremely important for TomoSAR deramping. Considering that the parallel baseline estimation is common in InSAR [28,29,30], it is not described here because this is not the key point of this study.

As noted in InSAR [28,29,30], the InSAR calculation method (interferometry and then flattening) is shown in Equation (6). $I_{mast}$ is the master SLC data; $I_{slave}$ is the slave SLC data; $(){ }^{*}$ means the conjugated calculation. In this way, we define that SLC data ($I$) are flattened, followed by Equation (7), which is the basic processing operation in software codes in many InSAR software.

$${\left\{I_{mast}·{\left(I_{slave}\right)}^{*}\right\}}_{flattening}=I_{mast}·{\left(I_{slave}\right)}^{*}· exp\left(-j·{\phi }_{sim{u}_m} \right)=I_{mast}·{\left( I_{slave}·exp\left(j·{\phi }_{sim{u}_m} \right) \right)}^{*},$$

(6)

$${\left\{I_{slave}\right\}}_{flattening}=I_{slave}·exp\left(j·{\phi }_{simu\_m}\right) ,$$

(7)

If Equation (2) is multiplied by $exp\left(-j\frac{4\pi }{\lambda }R\right)$, $R$ is the central phase slant range of the master image, then Equation (2) becomes Equation (8).

$$\begin{gathered} p\left({\xi }_m\right)=exp\left(-j\frac{4\pi }{\lambda }R\right).\mathrm{g}\left(m\right)=exp\left(-j\frac{4\pi }{\lambda }R\right)exp\left(j\frac{4\pi f_0}{c}{R}_m\left(0\right)\right)\cdot Q\left(m\right)= \\ exp\left(j·-\frac{4\pi }{\lambda }\left[R-R_m\left(0\right)\right]\right)· Q\left(m\right), \end{gathered}$$

(8)

From Equation (2), R_m(0) is the central slant range of each image (slave image) and if ${\phi }_m$ is defined as

$${\phi }_m=-\frac{4\pi }{\lambda }\left[R-R_m\left(0\right)\right],$$

(9)

Equation (2) becomes Equation (10).

$$p\left({\xi }_m\right)=exp\left(j· {\phi }_m\right)· Q\left(m\right),$$

(10)

It can be seen that ${\phi }_m$ is the interferometric phase. If the DEM error is ignored and ${\phi }_m$ is replaced by the DEM simulation phase, Equation (10) is identical to Equation (7). This is because $Q\left(m\right)=I\left(m\right)$, and they are all the same (mth) SLC data. In this way, based on Equation (2), the flattening of the SLC data ($I_{flattening}$) is shown in Equation (11).

$$\begin{gathered} I_{flattening}\left(m\right)=I\left(m\right)·exp\left(j·{\phi }_{simu\_m}\right)=Q\left(m\right)·exp\left(j·{\phi }_{simu\_m}\right) \\ =exp\left(j·{\phi }_{simu\_m}\right)· Q\left(m\right)=p\left({\xi }_m\right)=exp\left(-j\frac{4\pi }{\lambda }R\right)\mathrm{g}\left(m\right) \\ ={{\displaystyle \int }}_{s_{min}}^{s_{max}}exp\left(-j\frac{4\pi }{\lambda }R\right)\gamma \left(s\right)exp\left(j2\pi {\xi }_{m}s\right)ds, \end{gathered}$$

(11)

If ${\gamma }^{\prime }\left(s\right)$ is defined in Equation (12), the flattening of the SLC data in Equation (11) becomes Equation (13). Therefore, the flattening of SLC data ($I_{flattening}$) can be directly used in TomoSAR processing, followed by the use of core Equations (8)–(13). ${\gamma }^{\prime }\left(s\right)$ is the inversion results of TomoSAR processing along the s direction (NSR direction).

$${\gamma }^{\prime }\left(s\right)=exp\left(-j \frac{4\pi }{\lambda }R\right) \gamma \left(s\right),$$

(12)

$$I_{flattening}\left(m\right)=p\left({\xi }_m\right)={{\displaystyle \int }}_{s_{min}}^{s_{max}}{\gamma }^{\prime }\left(s\right) exp\left(j2\pi {\xi }_{m}s\right)ds,$$

(13)

2.2. SAR Interferometry Phase (InSAR) Calibration with DEM Error

The phase calibration step is indispensable for TomoSAR imaging. In this paper, we employed the newly released TanDEM-X 12 m DEM available for very small geographical areas, alongside a new method to correct the phase via the SAR interferometric phase with the DEM error. The 12 m TanDEM-X DEM data were obtained through a data grant from DLR for 3D TomoSAR imaging over Zipingpu Dam, Dujiangyan, Sichuan Province, China, and its location is shown in Figure 3. A hill shaded map by GMT5 is shown in Figure 4.

DEM data, like all other spatial datasets, have errors (systematic and random errors as well as blunders). Therefore, it is critical to validate DEM products before using them for TomoSAR processing. The SRTM 30 m (absolute vertical accuracy ≤16 m) [31] and ICESAT GLA14 data (absolute vertical accuracy ≤1 m) [32] are used to assess the TanDEM-X 12 m DEM data. The vertical comparison accuracy (the height difference between the DEM and the ‘truth’ DEM) of the TanDEM-X 12 m DEM data in Dujiangyan can be summarized as follows: against the ICESat GLAS14 elevation data, the TanDEM-X 12 m DEM had an accuracy of 0.6 ± 16.1 m (mean ± standard deviation), and the RMSE was 16.1 m; against the SRTM 30 m data, the TanDEM-X 12 m DEM had an accuracy of −0.4 ± 15.9 m and the RMSE was 15.9 m. The standard deviation and RMSE of the difference between the 12 m TanDEM-X DEM data and SRTM data was a little bit larger; this might be because the 12 m TanDEM-X DEM data had errors, noise, and a lot of changes (e.g., new bridges, a new dam, earthquake deformation, landslide deformation, and so on in this area) after 2000 when the SRTM data were collected. The 12 m TanDEM-X DEM data had random errors (mainly from high frequency noise) that needed filtering before using them. As the TanDEM-X 12 m DEM included holes, the SRTM 30 m data were used to replace the holes within the Environment for Visualizing Images (ENVI). After that, the DEM data were edited using the noise removal and smoothing tool in ENVI. Then, the new improved TanDEM-X 12 m DEM was used for TomoSAR imaging.

In this study, a high-resolution DEM (e.g., the quality improved TanDEM-X 12 m DEM) was employed as the reference DEM in the co-registration step of TomoSAR processing. Then, SAR interferometry (InSAR) was applied. The interferograms between the slave image (10/08/2016) and the master image (25/07/2016) are shown in Figure 6 and Figure 7, which demonstrate the elimination of the DEM phase. The fringes could be clearly seen on the building, road, and dam areas in these figures, while the fringes were not clear in the mountain tree areas as the coherence was very low; this was caused by foreshortening, shadows, and the SAR layover in these complex mountain tree areas.

In addition, as the DEM simulation phase was used for deramping, the deramping phase in Equation (9) changed to that of Equation (14). As the simulation phase by DEM is based on an external DEM, and DEM data have errors, the DEM error phase ${\phi }_{DEMerror}$ must be taken into consideration in TomoSAR processing.

$${\phi }_m=-\frac{4\pi }{\lambda }\left[R-R_m\left(0\right)\right]={\phi }_{simu\_m}+{\phi }_{DEMerror},$$

(14)

As is well-known, SLC complex data have errors, which include an atmospheric error (water vapor and ionosphere), orbit error, deformation (linear and nonlinear), and noise. The real complex data can be written as an ideal complex data plus all of the error phase ${\phi }_e$. ${\phi }_{e\_master}$ is the error phase of the master data and ${\phi }_{e\_slave}$ is the error phase of the slave data. In this way, the interferometric phase in Equation (14) becomes Equation (16). Based on Equations (8), (9), and (13), Equation (10) becomes Equation (17). In addition, the real flattening complex data are shown in Equation (15), then Equation (17) becomes Equation (18). It is known that ${\phi }_{InSAR}={\phi }_{DEMerror}+{\phi }_{e\_master}-{\phi }_{e\_slave}$ is the interferometric phase after DEM flattening, as shown in Equation (19). This is because the InSAR phase is the DEM error phase, plus the master error phase, minus the slave error phase after flattening. Therefore, based on Equations (15) and (19), the TomoSAR processing Equation (18) becomes Equation (20). It can be seen from these equations that the InSAR phase can be used for phase calibration. With phase errors and the DEM error phase, Equation (20) should be used for the TomoSAR processing, as shown below. In Equation (20), $I_{flattening\_real}\left(m\right)$ is the flattening complex data using an external DEM, and ${\phi }_{InSAR}$ is the interferometric phase after flattening by this DEM. ${\gamma }^{\prime }\left(s\right)$ in Equation (12) is the inversion results of TomoSAR processing along the s direction.${\xi }_m$ in Equation (4) is the spatial frequency, corresponding to the height $s$ from the normal-slant-range direction ($s$ elevation direction).

$$I_{flattening\_real}\left(m\right)=exp\left(j·{\phi }_{simu\_m}\right)·Q\left(m\right),$$

(15)

$${\phi }_m=-\frac{4\pi }{\lambda }\left[R-R_m\left(0\right)\right]={\phi }_{simu\_m}+{\phi }_{DEMerror}+{\phi }_{e\_master}-{\phi }_{e\_slave},$$

(16)

$$\begin{gathered} p\left({\xi }_m\right)=exp\left(j· {\phi }_m\right)· Q\left(m\right)=exp\left(j·{\phi }_{simu\_m}+j·{\phi }_{DEMerror}+j·{\phi }_{e_{master}}-j· \\ {\phi }_{e\_slave}\right)· Q\left(m\right)={{\displaystyle \int }}_{s_{min}}^{s_{max}}\gamma \prime \left(s\right)exp\left(j2\pi {\xi }_{m}s\right)ds, \end{gathered}$$

(17)

$$\begin{gathered} I_{flattening\_real}\left(m\right) exp\left(j·{\phi }_{DEMerror}+j·{\phi }_{e_{master}}-j·{\phi }_{e\_slave} \right)=p\left({\xi }_m\right) \\ ={{\displaystyle \int }}_{s_{min}}^{s_{max}}\gamma \prime \left(s\right)exp\left(j2\pi {\xi }_{m}s\right)ds \end{gathered}$$

(18)

$${\phi }_{InSAR}={\phi }_{DEMerror}+{\phi }_{e\_master}-{\phi }_{e\_slave},$$

(19)

$$I_{flattening\_real}\left(m\right) exp\left(j·{\phi }_{InSAR}\right)=p\left({\xi }_m\right)={{\displaystyle \int }}_{s_{min}}^{s_{max}}\gamma \prime \left(s\right)exp\left(j2\pi {\xi }_{m}s\right)ds,$$

(20)

After this phase calibration, the data are ready for super-resolution reconstruction via Capon and compressive sensing. However, the height reference is based on the strong scattering center after INSAR calibration. Hence, it is difficult to obtain a real height reference after this method. If geocoding is needed, control points are needed, or the result can be referenced to the DEM data after DEM error estimation and DEM error compensation based on Equation (23).

$$\begin{gathered} I_{flattening\_real}\left(m\right) exp\left(j·{\phi }_{InSAR}\right)·exp\left(-j·{\phi }_{DEMerror}\right)={{\displaystyle \int }}_{s_{min}}^{s_{max}}{\gamma }^{\prime }\left(s\right)· \\ exp\left(-j·{\phi }_{DEMerror}\right)exp\left(j2\pi {\xi }_{m}s\right)ds, \end{gathered}$$

(21)

$${\gamma }^{″}\left(s\right)={\gamma }^{\prime }\left(s\right)·exp\left(-j·{\phi }_{DEMerror}\right),$$

(22)

$$\begin{gathered} I_{flattening\_real}\left(m\right) exp\left(j·({\phi }_{InSAR}-{\phi }_{DEMerror})\right)= \\ {{\displaystyle \int }}_{s_{min}}^{s_{max}}{\gamma }^{″}\left(s\right)exp\left(j2\pi {\xi }_{m}s\right)ds, \end{gathered}$$

(23)

As shown in the above equations, ${\phi }_{InSAR}-{\phi }_{DEMerror}={\phi }_{e\_master}-{\phi }_{e_{slave}}$, and ${\phi }_{e\_master}-{\phi }_{e_{slave}}$ represent the difference of phase errors, which includes an atmospheric error term (water vapor and ionosphere), orbit error, deformation (linear and nonlinear), and noise, hence why this step is called phase calibration.

Moreover, as shown in Equation (22) for all of the slave images, and in the right part of Equation (23), the DEM error phase $exp\left(-j·{\phi }_{DEMerror}\right)$ goes to $\gamma \prime \left(s\right)$ and thus forms ${\gamma }^{″}\left(s\right)$, as shown in Equation (22) for all of the slave images. The DEM error phase does not influence the magnitude of $\gamma \prime \left(s\right)$, which means that the magnitude of $\gamma \prime \left(s\right)$ is the same as the magnitude of ${\gamma }^{″}\left(s\right)$. The magnitude of ${\gamma }^{″}\left(s\right)$ is the inversion results of TomoSAR processing. As a result, the InSAR with the DEM error phase method can be used for phase calibration. According to Equation (23), the accuracy of the DEM error estimation will influence the inversion results. After PS-InSAR (or SBAS InSAR), the DEM error can be obtained. Then, the DEM error of the whole image can be obtained by interpolation. Therefore, the DEM error was finally used to convert the unknown strong scattering center to the DEM data as the reference.

2.3. The Compressive Sensing Method

As SAR tomography is a semi-discrete problem ($\gamma \left(s\right)$ is continuous), it is necessary to discretise $\gamma \left(s\right)$ for actual processing so that it can be converted into a discrete problem for solution. Let $\Delta s$ be the sampling interval of the NSR direction (elevation s) in the discretization process in the sampling interval $\left[S_{min}, S_{max}\right]$, then the total samples ($N$) are obtained. Next, the SAR tomography model can be written in the following form.

$$\mathrm{g}\cong \mathsf{\Phi }\gamma +e,$$

(24)

where $\mathrm{g}$ is an $M\times 1$ dimensional observation vector,

$$\mathrm{g}={\left[\mathrm{g}\left({\xi }_1\right),\mathrm{g}\left({\xi }_2\right),\cdot \cdot \cdot ,\mathrm{g}\left({\xi }_M\right)\right]}^T,$$

(25)

The ${\xi }_m$ is the spatial frequency corresponding to the height s in the NSR direction, $\mathsf{\Phi }$ is a $M\times N$ dimensional matrix,

$$\mathsf{\Phi }={\left[{\varphi }_1,{\varphi }_2,\cdot \cdot \cdot ,{\varphi }_M\right]}^T,$$

(26)

$${\varphi }_M={\left[exp\left(j2\pi {\xi }_m{s}_1\right),exp\left(j2\pi {\xi }_m{s}_2\right),\cdot \cdot \cdot ,exp\left(j2\pi {\xi }_m{s}_N\right)\right]}^T,$$

(27)

$\gamma$ is an $N\times 1$ dimensional unknown signal vector, which is what we need to achieve at the end.

$$\gamma ={\left[\gamma \left(s_1\right),\gamma \left(s_2\right),\cdot \cdot \cdot ,\gamma \left(s_N\right)\right]}^T,$$

(28)

$e$ is an $M\times 1$ dimensional noise vector.

$$e={\left[e\left({\xi }_1\right)e\left({\xi }_2\right)\cdot \cdot \cdot e\left({\xi }_m\right)\right]}^T,$$

(29)

According to the model in Equation (24), SAR tomography essentially reconstructs the unknown signal $\gamma$ from the measured data $\mathrm{g}$. Since the dimension of the measured data is much smaller than the dimension $M<<N$ of the unknown signal $\gamma$, according to the classical signal reconstruction theory, a direct solution of Equation (24) is a pathological problem (cannot be solved in the classical signal reconstruction theory as the number of unknown parameters are larger than the number of the measurements). There are usually only a few strong scatterers in the same azimuth-distance resolution unit (one pixel), that is, $\gamma$ is sparse in the height field. In practice, there is some unavoidable clutter in addition to the main strong scatterers, whose intensities are far less than those of strong scatterers [33,34,35]. Therefore, generally speaking, $\gamma$ is compressible in the elevation region. In this case, the sparse basis is the Dirac basis, the sparse basis matrix $\mathsf{\Psi } = \mathrm{I}$, $\mathrm{I}$ is the $N\times N$ dimension unit matrix, $\gamma = \mathsf{\Psi } \gamma = \gamma$. In view of this, we can solve the problem of SAR tomography in the theoretical framework of compressive sensing. The result is

$$\underset{\gamma }{min}{\gamma }_1 \mathrm{subject} \mathrm{to} \mathrm{g} =\mathsf{\Phi }\gamma ,$$

(30)

For a given $\mathrm{K}$, since $\gamma$ is compressible in the elevation region, that is, the sparse basis matrix $\mathsf{\Psi } = \mathrm{I}$, $\gamma = \mathsf{\Psi }, \gamma = \gamma ,$ quadratic relaxation can provide a good approximation of the sparsest solution. When noise is considered, Equation (30) becomes

$$M=O\left(K log N\right),$$

(31)

$$\underset{\gamma }{min}{\gamma }_1 \mathrm{subject} \mathrm{to} \Vert \mathrm{g}-\mathsf{\Phi }\gamma \Vert \le {\sigma }_{e },$$

(32)

where $M$ is the smallest number of SLC data for CS method; ${\sigma }_e$ is the standard deviation of the noise. If $K$ is not known and measurement noise exists, Equation (30) can be approximated by:

$$\widehat{\gamma }=\arg \underset{\gamma }{min}(\left|\left|\mathrm{g}-\mathsf{\Phi }\gamma \right|\right|\begin{array}{c}2\\ 2\end{array}+{\lambda }_k{\left|\left|\gamma \right|\right|}_1),$$

(33)

where $\widehat{\gamma }=\arg \underset{\gamma }{min}()$ means the optimum estimate of $\gamma$; $\widehat{\gamma }$ is the best optimum estimate of $\gamma$, when $\left|\left|\mathrm{g}-\mathsf{\Phi }\gamma \right|\right|\begin{array}{c}2\\ 2\end{array}+{\lambda }_k{\left|\left|\gamma \right|\right|}_1$ is the smallest and $\gamma \approx \widehat{\gamma }$; ${\lambda }_k$ is a Lagrange multiplier depending on the number of samples $N$ [36] and the noise level ${\sigma }_e$. Equation (33) consists of an L2 norm residual and an L1 norm regularizer and Equation (33) can be interpreted as a Bayesian estimate with an exponential prior favoring sparse solution. By scaling down via the L1 and L2 norm minimization method, model selection, parameter estimation, and least-squares estimation, all the results can be achieved.

3. Results

3.1. Test Sites and COSMO-SkyMed Spaceborne SAR

The densely vegetated mountainous rural areas at Zipingpu Dam, Dujiangyan, Sichuan, China were employed as a test area, as shown in Figure 8a,b. Figure 8c displays a small test subarea over the Zipingpu Dam for TomoSAR imaging. In 3D SAR imaging (SAR tomography), 14 COSMO-SkyMed Spotlight ascending data stacks are used, as shown in

Table 1. The ascending COSMO-SkyMed Spotlight data stacks, the master image is 25/07/2016.

ID	Incidence Angle	Time	Track Type	Baseline (m)
1	37.66	03/06/2016	Ascending	−373.44
2	37.66	11/06/2016	Ascending	−289.25
3	37.66	19/06/2016	Ascending	743.28
4	37.66	23/06/2016	Ascending	920.38
5	37.66	05/07/2016	Ascending	−431.28
6	37.66	09/07/2016	Ascending	−629.15
7	37.66	25/07/2016	Ascending	0
8	37.66	06/08/2016	Ascending	820.31
9	37.66	10/08/2016	Ascending	203.46
10	37.66	22/08/2016	Ascending	−435.71
11	37.66	26/08/2016	Ascending	−140.25
12	37.66	07/09/2016	Ascending	186.17
13	37.66	11/09/2016	Ascending	611.54
14	37.66	23/09/2016	Ascending	−330.33

. Figure 8 reveals the TDX-DEM 12 m and SRTM 30 m data as well as the orbit information. In this test, the reference terrain-TanDEM-X 12 m data (Figure 8d) were used for deramping with all of the TomoSAR height results being based on the TanDEM-X DEM12 m surface. Work on the dam was started on 29/03/2001 and completed on 30/09/2005. However, the SRTM data were acquired in February 2000. Therefore, Figure 8f demonstrates the height difference and shape of the dam. Moreover, after PS-InSAR, the DEM errors of each PS point are estimated. Then, the DEM error of the whole image can be obtained by interpolation. The DEM error map of 03/06/2016 via PS-InSAR is shown in Figure 9. After estimating the DEM errors, the phase of this DEM error will be compensated via simple phase addition and subtraction for TomoSAR processing based on Equation (23).

3.2. TomoSAR Results and Fieldwork for Validation

The results derived from the compressive sensing [33] method in the sub-test area after InSAR phase calibration and DEM referenced correction is shown in Figure 10. All results were referenced to the DEM elevation, which can be applied for geocoding using SAR observation geometry. After CS imaging, the errors were filtered, and targets selected by the improved constant false alarm rate (CFAR) method for our application. First, the mean value and standard deviation were calculated, and all peaks were found along the height direction; then, when the amplitude of the peaks was larger than three times the standard deviation plus the mean value, they were selected as targets. Moreover, the results illustrated a strong backscatter curve line with the few targets around it, confirming that the X band data could hardly penetrate trees. The simulation using the COSMO-SkyMed X-band satellite parameters with the baseline distribution of the dataset used demonstrated that the theoretical vertical precision of the extraction elevation was 0 m.

Finally, after geocoding (adding the TanDEM-X 12 m DEM height), the 1 m TomoSAR elevation extraction results of the X-band COSMO-SkyMed Spotlight data can be obtained, as shown in Figure 10. Compared to the TanDEM-X 12 m DEM data in Figure 10a, the results of the 1 m TomoSAR in Figure 10b exhibited better resolution. As shown in Figure 10b, there are two line steps on the dam in the 1 m TomoSAR imaging results, which can also be observed in the Lidar and photos in Figure 11a–c.

As the X band has little vegetation penetration capability and no penetration for the X band at all for the man-made dam, only one point of the TomoSAR imaging result for each SAR pixel along the height direction in the dam area could be validated using Lidar data. Furthermore, as the TomoSAR results were referenced to the TANDEM-X 12 m data, the results already have a set of absolute coordinates after geocoding (adding DEM height). If the high accuracy of the position is needed, the Lidar or other control data can be used as control points for geocoding.

Fieldwork to acquire terrestrial Lidar data and photos were collected at Zipingpu Dam, Dujiangyan, China between 24/09/2017 and 29/09/2017. The photos of the Zipingpu dam are shown in Figure 11a,b. The Lidar data collected by the ground-based V-Line 3D Terrestrial Laser Scanner (TLS) RIEGL VZ-1000 at Zipingpu Dam were used for validation, as shown in Figure 11c.

Finally, and Figure 12 and

Table 2. The results of the validation statistics.

Basic Stats	Area	Min (m)	Max (m)	Mean (m)	Stdev σ (m)	RMSE (m)
TomoSAR-Lidar	Dam and mountain trees	−8.5	8.6	0.11	2.81	2.82
TomoSAR-Lidar	Dam	−5.7	5.3	0.25	1.04	1.07

indicate the absolute difference map between the X-band TomoSAR imaging result and the Lidar data over the field site at Dujiangyan and the height difference statistical results. The map shows that the compressive sensing result had a good match with Lidar data after geocoding. As shown in Table 2, the maximum difference was 8.6 m, the minimum difference was −8.5 m, the mean difference was 0.11 m, the standard deviation was 2.81 m, and the RMSE was 2.82 m. After masking out the tree and mountain areas, the maximum difference of the dam area was 5.3 m, the minimum difference was −5.7 m, the mean difference was 0.25 m, the standard deviation was 1.04 m, and the RMSE was 1.07 m (over the dam area). These results demonstrate that the compressive sensing result showed a good match with the Lidar data after geocoding. Therefore, the TomoSAR imaging algorithm appears to work satisfactorily and result in a 1 m DEM. These high resolution TomoSAR results in the dam area are outstanding for application because they match very well with the Lidar data and have high accuracy.

Finally, we show several azimuth lines along the range direction (points of each line are stored in Google Earth KML format, the red colour points in the picture) of our TomoSAR results in Google Earth, shown in Figure 13 below.

4. Discussion

First, all of the SAR SLC data were co-registered in the radar coordinate system to sub-pixel accuracy for TomoSAR processing. This was based on the master image orbit information and external DEM. Thereafter, a precise orbit was used to estimate the orbit baseline. Moreover, because DEM deramping establishes the spectrum relationship (based on the perpendicular baseline) between the observed SAR images and 3D SAR tomographic reconstruction inversion parameters, 3D SAR inversion is possible. The reference slant range used in the deramping process can be the central range of each image using the radar transmission time delay and velocity, or the distance between the radar antenna phase center of each image and a reference terrain (e.g., external DEM data). It should be noted that the only differences between the two methods are the zero-point benchmark and phase errors. According to our experiment, the coarse DEM (TanDEM-X 12 m) was used for deramping without additional atmospheric phase errors in the accurate reference slant range calculation. However, even after pre-processing of the high resolution DEM, the DEM still had elevation errors, which may also have an impact on the TomoSAR results. Therefore, the phase calibration step is indispensable for TomoSAR imaging, eventually ensuring the accuracy of the inversion results. After PS-InSAR (or SBAS InSAR), the DEM error was estimated and obtained. Then, the phase calibration method based on SAR interferometry (InSAR) and DEM error estimation and compensation for 3D TomoSAR imaging was used. Based on this method, the accuracy of the DEM error estimation will affect the inversion results. However, the high accuracy of the DEM error estimation in a constructed facility area can guarantee good performance for the phase calibration.

In the Dujiangyan test area, it was difficult to obtain favorable compressive sensing results using the original baseline, while compressive sensing worked successfully after using the improved baseline with the new phase calibration. The compressive sensing results in the dam area had a good match with the Lidar data. Although there were some differences except for the 0 m difference, they were around 0.18 m difference on average, which might be caused by other errors such as DEM uncertainty, orbital, tropospheric and ionospheric phase distortion, and thermal noise.

The CS algorithm is generally executed using SLC data (without averaging for noise filtering), which is automatically performed in the stack dimension using the SAR model and sparsity driven estimation technique to characterize a minimal number (<4) of targets in the measured signal. Therefore, CS is well-adapted to discrete scatterers or targets as well as 3D and 4D high resolution tomographic reconstruction. The results demonstrate that the X band had little penetration capability, and thus cannot be used for forest structure reconstruction. However, it can be used to extract the canopy top, the shape of the constructed structures (dam, buildings, and manufactured facilities), and the top of the 3D terrain, which is optimal for high resolution DSM extraction and target detection.

5. Conclusions

In this paper, based on our theoretical study and mathematical derivation, a systematic TomoSAR algorithm and the methods were described, demonstrated, tested, and analyzed to achieve 3D tomographic SAR imaging results using COSMO-SkyMed X band data over Zipinpu Dam, Dujiangyan, Sichuan, China. We demonstrated that in 3D SAR tomography, phase calibration plays a crucial role in the data processing. With the newly released TanDEM-X 12 m DEM and suitable pre-processing, this study proposes a new phase calibration method based on SAR interferometry (InSAR) and DEM error estimation and compensation for 3D TomoSAR imaging. Then, the super-resolution reconstruction CS method was studied, which demonstrated that the X band data with the CS method are not suitable for forest structure reconstruction, but fit for the reconstruction of manufactured facilities (dam) and buildings in the urban area. Finally, the fieldwork and validation were described, and the results revealed that our methods worked properly using COSMO-SkyMed X band data at Zipinpu Dam, Dujiangyan, Sichuan, China.