High-resolution High-dimensional Imaging of Urban Building Based on GaoFen-3 SAR Data（in English）

1. Introduction

Synthetic Aperture Radar (SAR), unlike traditional optical observation methods, has all-day and all-weather surveillance capability, making it widely used in fields such as land and resource surveying and natural disaster monitoring^[1]. Traditional SAR imaging, on the other hand, can only obtain an azimuth-range two-Dimensional (2-D) image and cannot accurately reflect the three-Dimensional (3-D) scattering characteristics of the target, which limits the SAR image’s further application. TomoSAR (Synthetic Aperture Radar Tomography) is an extension of SAR imaging technology. It extends the synthetic aperture principle in the elevation direction, allowing it to obtain azimuth-range-elevation information from multiple 2-D complex-valued images and thus achieve 3-D imaging^[2]. Reigber et al.^[3,4] presented the first airborne TomoSAR imaging results in 1999, achieving 3-D imaging of forest regions using spectral estimation methods. Fornaro and Serafino^[5] carried out the spaceborne experiment of TomoSAR imaging in 2006, using long-term baseline data, and confirmed the feasibility and effectiveness of applying spectral estimation technology to elevation reconstruction. Reale et al.^[6] demonstrated in 2011 that high-resolution data combined with advanced interference-processing techniques could effectively reconstruct the 3-D structure of buildings. Shahzad and Zhu^[7] proposed a new building facade reconstruction method in 2013. It obtains the 3-D point cloud of the Bellagio Hotel in Las Vegas, confirming the TomoSAR point cloud’s superiority in building a dynamic city model. Wang et al.^[8] proposed an iterative reweighted alternating direction multiplier method for fast TomoSAR imaging in 2017. Wang and Zhu^[9] proposed a TomoSAR imaging method based on kernel principal component analysis in 2018, which uses minimal cost to separate multiple scatterers distributed along the elevation direction in the same azimuth-range resolution cell. Qin et al.^[10] introduced a machine learning-based building target recognition and extraction algorithm in 2019 to address the problems of low elevation resolution, layover, and low efficiency of building target feature extraction in TomoSAR. It increased the efficiency of feature extraction of observation targets and validated the method’s effectiveness using airborne array 3-D SAR data. D-TomoSAR (Differential Synthetic Aperture Radar Tomography) is a TomoSAR extension. It is possible to achieve azimuth-range-elevation-time four-Dimensional (4-D) imaging by using multi-baseline observation data. It not only solved the problem of high dislocation and ambiguity in TomoSAR, but it also obtained highly accurate target deformation information. Lombardini first proposed the D-TomoSAR concept in 2003^[11]. Fornaro et al.^[12] demonstrated in 2007 that D-TomoSAR could be regarded as an effective alternative to traditional permanent scatterer monitoring techniques, allowing effective monitoring of large-scale scene deformation. Fornaro et al.^[13] proposed and applied a technique for separable relevant time series of interfering scatterers to D-TomoSAR in Rome in 2008. It demonstrated that this technique could be used to monitor the deformation of urban complex scenes. Zhu et al.^[14] achieved D-TomoSAR imaging of Las Vegas using TerraSAR-X data in 2009, resulting in an urban 4-D map. Fornaro et al.^[15] used ERS data for D-TomoSAR imaging of Rome’s Grotta Perfetta area in 2010, demonstrating that the D-TomoSAR method can effectively distinguish different scatterers in the same resolution cell, overcome the limitations of differential interference, and thus improve the monitoring capacity of urban infrastructure deformation. Zhu and Bamler^[16] proposed a “time-distorted” method and implemented D-TomoSAR imaging of urban areas based on TerraSAR-X data in 2011, effectively obtaining the linear and seasonal deformation velocity of urban buildings. Siddique et al.^[17] used 50 TerraSAR-X images in 2015 to combine the D-TomoSAR method and permanent scatterer interference technology to obtain spatial-temporal inversion results of high-rise buildings. Wang and Liu^[18] proposed a generalized D-TomoSAR imaging system model and a quasi-maximum likelihood-based algorithm that inversed multiple deformation velocities, including linear and seasonal deformation velocities, in 2020. At the moment, TomoSAR and D-TomoSAR imaging technologies have shown great application potential in 3-D reconstruction and long-term deformation monitoring of urban buildings and infrastructures.

Compressive Sensing (CS) is an important sparse signal processing technology that can achieve high-quality recovery for sparse signals from fewer samples^[19-21]. The observed scene’s elevation distribution must be sparse for CS-Tomo-SAR to work. Because urban areas are mostly made up of man-made structures, their elevation distributions satisfy the sparsity condition. As a result, CS has a promising application in 3-D and 4-D urban imaging. Initially, the elevation resolution of TomoSAR is far lower than the azimuth or range direction due to the orbit limitation of modern meter-resolution spaceborne SAR systems. As a result, there is an urgent need for a super-resolution algorithm to solve this problem. Zhu and Bamler^[22] introduced a TomoSAR imaging method based on CS theory in 2010. It is discovered that, when compared to traditional spectral estimation methods, the proposed algorithm can achieve super-resolving elevation distribution reconstruction. In 2012, Zhu and Bamler^[23] applied the proposed SL1MMER algorithm to Terra-SAR-X spaceborne data processing and obtained a high-resolution 3-D image of the Bellagio Hotel in Las Vegas, validating the SL1MMER algorithm’s super-resolution capability. Weiss et al.^[24] proposed an adaptive CS algorithm for Tomo-SAR in 2015, which accurately identifies two scatterer positions in the same resolution cell. In 2017, Li et al.^[25] investigated the SPICE-based TomoSAR imaging method and used 8 Terra-SAR-X strip map images to create a high-precision 3-D reconstruction of a building in Genhe, Inner Mongolia. Zhu and Bamler^[26] applied the CS technique to D-TomoSAR imaging and demonstrated its superiority over conventional spectral estimation methods on multi-scatterer separation in elevation direction in 2010. It is also mentioned that CS can automatically determine the number of scatterers, making it ideal for 3-D and 4-D imaging of spaceborne SAR systems. In 2010, Zhu and Bamler^[27] used the SL1MMER algorithm to obtain D-TomoSAR imaging results of the Las Vegas Convention Center, validating the advantages and capabilities of CS in 4-D imaging. Leng et al.^[28] (2014) demonstrated the elevation reconstruction and deformation monitoring results of Barcelona by applying the least absolute shrinkage and selection operator CS algorithm to building regions.

The Gaofen-3 (GF-3) satellite is China’s first C-band multi-polarization 1-meter resolution SAR satellite. On August 10, 2016, it was launched by a CZ-4C carrier rocket^[29]. The “China High-resolution Earth Observation Project” includes only one civil microwave remote sensing imaging satellite, the GF-3. It has high-resolution, wider swath, multi-imaging mode, and long-period operation capabilities and can perform all-day and all-weather global ocean and land monitoring^[30]. The GF-3 is now successfully used in high-precision surveying and mapping, natural disaster monitoring, and other fields^[31,32]. However, because high-dimensional imaging capability was not considered during system design, the collected GF-3 SAR images have Spatio-Temporal de-coherence, which poses some challenges for their future applications to InSAR, D-InSAR, TomoSAR, and D-TomoSAR. Yu et al.^[33] used the GF-3 data data to conduct interferometric experiments in the areas surrounding Dengfeng city in Henan Province, China, in 2019. They verified the interference capability of the GF-3 and demonstrated that it could be used to extract surface deformation information by comparing experimental results with the Sentinel-1 satellite. Huang et al.^[34] successfully extracted and obtained the digital elevation model of the surveillance region using interferometric experiments based on the GF-3 data in July 2021.

In this paper, we used seven complex-valued SAR images of the GF-3 to perform TomoSAR and D-TomoSAR imaging experiments with CS technology, resulting in high-resolution 3-D and 4-D images of the buildings surrounding Yanqi Lake in Beijing, China. Based on a Chinese satellite, we achieved high-quality 3-D reconstruction and high-precision deformation monitoring, which provides technical support for subsequent interference series applications and multi-dimensional high-resolution imaging of the GF-3.

The remainder of this paper is structured as follows. Section 2 provides a brief overview of the TomoSAR and D-TomoSAR imaging models, as well as the CS solution for the first two models. The GF-3 dataset, which is used in this paper, is introduced in Section 3. Section 4 performs TomoSAR and D-TomoSAR imaging experiments on simulated data to demonstrate the effectiveness of CS in high-resolution 3-D imaging and high-precision deformation monitoring. In Section 5, we performed TomoSAR and D-TomoSAR imaging studies of the buildings around Yanqi Lake in Beijing, using seven complex-valued SAR images of the GF-3, and obtained 3-D and 4-D images of representative buildings and large-scale areas. Section 6 contains the conclusions.

2. Imaging Model

2.1 TomoSAR imaging model

TomoSAR performs aperture synthesis along the elevation direction using multiple registered 2-D complex-valued SAR images of the same area, obtaining the target’s 3-D scattering information^[3,4,35,36]. Fig. 1 depicts the imaging geometry of TomoSAR. Let N denote the number of data acquisition baselines and ${b_n}\left( {n = 1,2,\cdots,N} \right)$ denote the elevation aperture. The focused complex-valued measurement at the ${n_{{\text{th}}}}$ SAR image can be expressed for a given azimuth-range pixel as

Figure  1.  TomoSAR imaging geometry

DownLoad: Full-Size Img PowerPoint

where ${\xi _n} = - 2{b_n}/\left( {\lambda r} \right)$ denotes the elevation frequency, $\lambda$ denotes the wavelength, r denotes the slant range, $\gamma \left( s \right)$ denotes the complex reflection function along the elevation direction s, and $\Delta s$ denotes the elevation span. After discretizing the complex reflectivity function $\gamma \left( s \right)$ along the elevation direction s, the imaging model in Eq. (1) can be approximated as

where L denotes the number of points along the elevation direction, ${\boldsymbol{g}} = {\left[ {{g_1},{g_2},\cdots,{g_N}} \right]^{\rm{T}}}$ denotes the measurement vector, ${\boldsymbol{R}} = \exp \left( { - {\rm{j}}2\pi {\xi _n}{s_l}} \right)$ denotes the observation matrix constructed based on the TomoSAR imaging geometry, and ${\boldsymbol{\gamma}} = \left[ \gamma \left( {{s_1}} \right), \gamma \left( {{s_2}} \right),\cdots,\gamma \left( {{s_L}} \right) \right]^{\rm{T}}$ denotes the discrete complex reflection function along ${s_l}\left( {l = 1,2,\cdots,L} \right)$ . From Eq. (2), the TomoSAR imaging model can be considered an irregular sampling-based discrete Fourier transform of $\gamma \left( s \right)$ . A SAR measurement can thus be thought of as a spectral parameter of the target complex reflection function along the elevation direction. For the non-parametric spectral analysis problem, the theoretical elevation resolution ${\rho _s}$ is determined by the elevation aperture $\Delta b$ . In the case of dense sampling in the elevation, ${\rho _s}$ can be calculated by

2.2 D-TomoSAR imaging model

When compared to TomoSAR, D-TomoSAR adds another time dimension. It combines two apertures into elevation and deformation velocity directions to achieve joint resolution, yielding a 4-D image of the observed target^[11]. For N complex-valued images, when the elevation aperture position is ${b_n}$ and time baseline is ${t_n}$ , the focused measurement of the ${n_{{\text{th}}}}$ acquisition can be expressed as

where ${\eta _n} = - 2\pi {t_n}/\lambda$ denotes the temporal frequency and $V\left( s \right)$ denotes the deformation velocity. The model in Eq. (4) also can be written as

where $\Delta v$ denotes the deformation velocity span of the observed target and $\delta \left( \cdot \right)$ denotes the spectral distribution induced by the deformation. Let ${a_\gamma }\left( {s,v} \right) = \gamma \left( s \right)\delta \left( {v - V\left( s \right)} \right)$ , then we can rewrite Eq. (5) as:

The above model can be considered a 2-D Fourier transform of $a\left( {s,v} \right)$ on the elevation-deformation plane. Thus, its projection in the vertical direction is $\gamma \left( s \right)$ ^[37]. After discretizing s and v in Eq. (6), the D-TomoSAR imaging model can be expressed as

where ${\boldsymbol{g}} = {\left[ {{g_1},{g_2},\cdots,{g_N}} \right]^{\rm{T}}}$ denotes the measurement vector, ${\boldsymbol{R}} = \exp \left( { - {\rm{j}}2\pi \left( {{\xi _n}{s_l} + {\eta _n}{v_q}} \right)} \right)$ is the D-TomoSAR observation matrix, ${s_l}\left( {n = 1,2,\cdots,L} \right)$ denotes the discrete elevation distribution, ${v_q}\left( {q = 1,2,\cdots,Q} \right)$ denotes the discrete distribution in the deformation direction, and $\gamma$ consists of a discretized $a\left( {s,v} \right)$ . Then the deformation resolution ${\rho _v}$ can be calculated with the time aperture size of $\Delta {\text{t}}$ by

2.3 CS reconstruction

The observed target in the city region is primarily man-made architecture, with a sparse elevation distribution, i.e., finite scatterers in each azimuth-range resolution cell. Thus, for the models in Eqs. (2) and (7), when the measurement matrix R satisfies the Restricted Isometry Property (RIP) condition, we can achieve CS-based TomoSAR and D-TomoSAR imaging by solving the following equation:

where $\beta$ denotes the regularization parameters, determined by the noise level and the number of samples. The CS algorithm quickly recovers high-quality signals from acquired samples^[38,39]. This paper employs the CS algorithm for TomoSAR and D-TomoSAR imaging due to its benefits.

3. Introduction to the GF-3 Dataset

The GF-3 dataset used in this paper is made up of seven complex-valued images, the parameters of which are listed in Tab. 1. The elevation aperture size of the seven images in this dataset is 1417 m. The time baseline spans 464 d, from June 2018 to September 2019. Fig. 2 depicts the spatial-temporal baseline distribution and the specific parameters are listed in Tab. 2. The SAR image acquired on March 1, 2019, serves as the master image in this paper, with the remaining six serving as slave images. By the way, the slave images’ spatial-temporal baseline position is calculated relative to the master image. The simulations in Section 4 of this paper will be run using the parameters listed in Tab. 1. Section 5 will carry out experiments using the GF-3 datasets.

Table  1.  Parameters of the GF-3 dataset

Parameter	Value		Parameter	Value
Spatial baseline span	1417.4 m		Number of scenes	7
Temporal baseline span	464 d		Azimuth resolution	0.3626 m
Slant range	1052747 m		Range resolution	0.765692 m
Wavelength	0.056 m		Elevation resolution	20.6174 m
Incident angle	47.2330015°		Information resolution	21.8 mm/year

 | Show Table

DownLoad: CSV

Figure  2.  Spatio-temporal baseline distribution of the GF-3 dataset

DownLoad: Full-Size Img PowerPoint

Table  2.  Spatio-temporal baseline parameters of the GF-3dataset

Number	Time of acquisition	Spatial baseline (m)	Temporal baseline (d)
1	2018.06.13	–459.108	–261
2	2019.01.31	–628.551	–29
3	2019.03.01	0	0
4	2019.03.30	–724.517	29
5	2019.07.24	692.863	145
6	2019.08.22	–38.211	174
7	2019.09.20	–510.491	203

 | Show Table

DownLoad: CSV

4. Simulations

We set two scatterers in the elevation direction using the experimental parameters listed in Tab. 1, generated seven scenes of simulated data, and added noise with SNR = 20 dB. TomoSAR and D-TomoSAR imaging use conventional spectral estimation methods and the CS algorithm, respectively^[40]. This section will present the simulation results of three spectral estimation algorithms, namely BeamForming (BF)^[41], Adaptive beamforming (Capon)^[42], and MUltiple SIgnal Classification (MUSIC)^[43,44], which will be compared to the results of CS to demonstrate the advantages of the CS algorithm in TomoSAR and D-TomoSAR imaging. The TomoSAR imaging results of two scattering points in the elevation directions are shown in Fig. 3. The elevation position distribution is represented by the x-coordinate, and the scattering point’s amplitude value is represented by the y-coordinate. The distances between the two scatterers were set to 11 m and 50 m, respectively. When the distance is less than 11 m, the three spectral estimation methods are unable to accurately separate the two scatterers, resulting in a failed reconstruction. At the same time, the CS algorithm is still capable of effectively distinguishing two scatterers to achieve elevation super-resolution imaging. When the distance is 50 m, the reconstruction result of the BF algorithm has significant and irregular sidelobes. Based on the results, the reconstruction images of the Capon and MUSIC algorithms have lower sidelobes than BF, which improves the quality of the elevation reconstruction. When compared to the other three spectral estimation algorithms, CS can effectively suppress sidelobes and noise while also improving elevation resolution. The TomoSAR imaging results for a pixel with three scatterers in the elevation direction are shown in Fig. 4. The distance between three scatterers is 20 m, which is comparable to the imaging results obtained with two scatterers. As can be seen, the CS algorithm works well for multiple scatterer separation.

Figure  3.  TomoSAR reconstructed reflectivity profiles of two scattering points along the elevation direction (left image: the distance between two scatterers is 11 m; right image: the distance between two scatters is 50 m)

DownLoad: Full-Size Img PowerPoint

Figure  4.  TomoSAR reconstructed reflectivity profiles of three scattering points along the elevation direction (the distance between three scatterers is 20 m)

DownLoad: Full-Size Img PowerPoint

Following that, several D-TomoSAR simulations were run. It is also assumed that the azimuth-range resolution cell contains two scatterers, located at –10 m and 10 m, respectively, and that the distance between them is less than the elevation resolution. The deformation velocities are set to 4 mm/year and –7 mm/year. We obtained the corresponding simulated data and added an additive noise with SNR = 20 dB based on the parameters in Tab. 1. Fig. 5 depicts the D-TomoSAR simulation results of various algorithms. The x-coordinate represents the target elevation position, and the y-coordinate represents the scatterer’s linear deformation variable. According to Fig. 5, three spectral estimation methods could not distinguish the two scatterers in the elevation direction due to algorithm limitations, particularly the BF algorithm, which is difficult to break through the theoretical resolution. The CS algorithm can accurately distinguish two scatterers separated by 20 m, which is consistent with the above results and further validates CS’s super-resolution ability in 3-D imaging. According to Tab. 1, the theoretical resolution of deformation is about 21.8 mm/year. While only the CS algorithm is capable of accurately estimating deformation rates of 4 mm/year and –7 mm/year, it confirms its superiority over spectral estimation techniques in deformation monitoring.

Figure  5.  D-TomoSAR simulation results (elevation position of two scatterers are –10 m and 10 m; deformation velocity of two scatterers are 4 mm/year and –7 mm/year, respectively)

DownLoad: Full-Size Img PowerPoint

5. Results of the GF-3 Dataset

5.1 Representative building inversion

We present the TomoSAR and D-TomoSAR reconstruction results based on the GF-3 dataset in this section. The Haihuayundu Eco-agriculture Ltd and the Yanqi Lake International Convention and Exhibition Center in Beijing were chosen as two representative buildings in the observation scene for TomoSAR and D-TomoSAR imaging. The optical image of the Haihuayundu Eco-agriculture Ltd is shown in Fig. 6(a), and the area framed by the red dashed line is the target of interest. The 2-D SAR image of the Haihuayundu Eco-agriculture Ltd is shown in Fig. 6(b). The height information of five buildings in the area is obtained by TomoSAR reconstruction, as shown in Fig. 7(a). The heights of the five individual buildings are all 35 m, which is consistent with the actual height, as shown in Fig. 7(a). It demonstrated the efficacy of the CS-TomoSAR imaging technology and demonstrated that it could be used for high-precision 3-D imaging of the GF-3 dataset. The reconstruction result of D-TomoSAR based on the CS algorithm is shown in Fig.7(b). The five buildings have different deformations, which may be related to the seasonal thermal expansion and contraction of building materials. It should be monitored daily to avoid the risks associated with it. A 3-D point cloud of the Haihuayundu Eco-agriculture Ltd, which more accurately reflects the 3-D scattering structure of the buildings, is shown in Fig. 8 to provide a visual representation of the interested area.

Figure  6.  The Haihuayundu Eco-agriculture Ltd

DownLoad: Full-Size Img PowerPoint

Figure  7.  Elevation and deformation velocity maps of the Haihuayundu Eco-agriculture Ltd

DownLoad: Full-Size Img PowerPoint

Figure  8.  3-D point cloud of the Haihuayundu Eco-agriculture Ltd

DownLoad: Full-Size Img PowerPoint

The optical and SAR images of the Yanqi Lake International Convention and Exhibition Center, with an actual height of about 30 m, are shown in Fig. 9. The CS-TomoSAR imaging result is shown in Fig. 10(a), and its 3-D point cloud image is shown in Fig. 11. The reconstructed image clearly reflects the 3-D structure of the building, and the height of the reconstructed building corresponds to the actual information. It also demonstrates that we can perform high-precision 3-D reconstruction of complex buildings using the GF-3 SAR data. The D-TomoSAR result of the Yanqi Lake International Convention and Exhibition Center based on CS technology is shown in Fig. 10(b). According to Fig. 10(b), the building’s deformation velocity ranges between 10 mm/year and –10 mm/year. The linear deformation rate of the left and right halves of the building is opposite, indicating that the ground beneath the building may be changing between uplift and collapse, which should be monitored more closely.

Figure  9.  Beijing Yanqi Lake International Convention and Exhibition Center

DownLoad: Full-Size Img PowerPoint

Figure  10.  Elevation and deformation velocity maps of Beijing Yanqi Lake International Convention and Exhibition Center

DownLoad: Full-Size Img PowerPoint

Figure  11.  3-D point cloud of Beijing Yanqi Lake International Convention and Exhibition Center

DownLoad: Full-Size Img PowerPoint

5.2 3-D and 4-D imaging of large-scale scene

The TomoSAR and D-TomoSAR imaging results of a large-scale scene from the GF-3 dataset are shown in this section. The optical and SAR images of Dingxiumeiquan Town are shown in Fig. 12. The TomoSAR imaging results based on CS are shown in Fig. 13(a). The heights of the buildings in the area are all between 15 and 20 m. Because the actual buildings have 4 to 6 floors, the results are consistent. The D-TomoSAR imaging results are shown in Fig. 13(b). It demonstrates that building deformation in the lower left and upper right areas of Dingxiumeiquan Town is approximately –10 mm/year, most likely due to construction in these areas causing ground subsidence. The deformation rate in the center of Dingxiumeiquan Town is approximately 5 mm/year and is relatively stable. These experiments demonstrate that by utilizing the GF-3 dataset, we were able to achieve the high-quality 3-D reconstruction of large-scale scenes as well as high-precision deformation monitoring, thereby confirming the GF-3 satellite’s application potential in urban sensing and monitoring.

Figure  12.  The Dingxiumeiquan Town

DownLoad: Full-Size Img PowerPoint

Figure  13.  Elevation and deformation velocity maps of the Dingxiumeiquan Town

DownLoad: Full-Size Img PowerPoint

6. Conclusions

In this paper, we conduct TomoSAR and D-TomoSAR imaging research using a GF-3 SAR complex-valued image dataset, obtain 3-D and 4-D images of two representative buildings in Beijing’s Yanqi Lake area, and show 3-D and 4-D SAR images of a large-scale area. These studies validate the potential of the Chinese GF-3 satellite for TomoSAR and D-TomoSAR imaging applications, as well as provide technical support for the GF-3’s subsequent expansion into interference series applications.

In the future, our team will continue relevant research and collect the GF-3 data from multiple areas and scenes. We will also investigate new high-resolution TomoSAR and D-TomoSAR imaging algorithms to realize large-scale 3-D and 4-D reconstruction of urban complex scenes, as well as further investigate GF-3’s enormous potential in interferometric series applications.