3D Point Cloud Reconstruction Using Inversely Mapping and Voting from Single Pass CSAR Images

Abstract

3D reconstruction has raised much interest in the field of CSAR. However, three dimensional imaging results with single pass CSAR data reveals that the 3D resolution of the system is poor for anisotropic scatterers. According to the imaging mechanism of CSAR, different targets located on the same iso-range line in the zero doppler plane fall into the same cell while for the same target point, imaging point will fall into the different positions at different aspect angles. In this paper, we proposed a method for 3D point cloud reconstruction using projections on 2D sub-aperture images. The target and background in the sub-aperture images are separated and binarized. For a projection point of target, given a series of offsets, the projection point will be mapped inversely to the 3D mesh along the iso-range line. We can obtain candidate points of the target. The intersection of iso-range lines can be regarded as voting process. For a candidate, the more times of intersection, the higher the number of votes, and the candidate point will be reserved. This fully excavates the information contained in the angle dimension of CSAR. The proposed approach is verified by the Gotcha Volumetric SAR Data Set.

1. Introduction

Traditional synthetic aperture radar (SAR) can obtain projection of the target on the two-dimensional plane, which loses the elevation information and limits the wide application of SAR imagery. 3D reconstruction of the observed object plays an important role in the SAR image applications [1,2]. In the field of computer vision, it is a fundamental problem to reconstruct 3D structures from single image [3,4] or image sequences taken at different viewpoints [5,6,7]. It is known that under different observation angles, different features of the target can be obtained. When the number of observation angles increases, more information of the target can be obtained, which is more conducive to 3D reconstruction [8].

Circular SAR (CSAR) moving along a circular trajectory [9], makes it possible to provide features of the target from different aspect angles. For isotropic scatterers, a single circular flight can achieve very high resolution and 3D imaging [10]. However, the 3D resolution of the anisotropic scatterers is very poor. The impulse response function (IRF) and consequently the retrieved information is distorted by undesired sidelobes [11,12]. It is almost impossible to achieve 3D imaging for anisotropic scatterers.

It has been deeply investigated for digital elevation map (DEM) generation using CSAR images by radargrammetry [13,14]. The key step is to measure the parallax between the two corresponding image points. However, due to the anisotropic nature of the object [15], it is difficult to find the corresponding points in the sub-aperture images under different aspect angles. Meanwhile the elevation information provided by DEM is insufficient for a detailed description of the target. Then from the work in [16,17], after polarimetric processing, the result provides rough shape of the vehicle outline from multi-view layover. Interferometric SAR (IFSAR) is a method of sparse 3D image reconstruction [18]. For two SAR images obtained at two elevation angles with little difference, it extracts height of the target by phase difference of the two images. These ways can give a better description of the target. Above methods all need to use the polarization information or multiple trajectory data. They have certain requirements for the type of data acquisition, and increase the difficulty of data processing.

In fact, we believe that the dimension of aspect angle can provide a rich amount of information. The research idea of this paper is to make full use of the image features under different aspect angles using single pass of CSAR images and try to avoid image registration step.

For an elevated target at the scene, the shape and location of projection of the target are related to the transient velocity vector of the SAR platform and the target itself. To be specific, if the target locates at the scene center, the projection of the target will be in different extent due to the height of the target. If the target is away from the scene center and is near to the radar platform, the projection becomes little. If the 3D coordinates of the target in the actual scene can be found by projection in each sub-aperture image, it is feasible to retrieve 3D shape of the target from CSAR images.

We proposed a 3D reconstruction algorithm to build the 3D point cloud of the object based on the projections in sub-aperture images provided by CSAR in this paper. Firstly, the target and background are separated by image segmentation and binarized. There is a certain offset between the projection point and the target in 3D space. Then based on different offsets, a projection point in the sub-aperture image is inversely mapped into the 3D mesh along the iso-range line. The above operations are performed on all sub-aperture images and the intersection times are accumulated in the 3D mesh. This process is defined as “voting” for each grid point, which is consistent with the idea in [19]. If the number of votes is greater than a certain threshold, the point is considered to be existing in 3D space. Otherwise it does not exist. Finally we can obtain the 3D point cloud from CSAR images.

Here the separation of objects and scene will affect the result of 3D reconstruction [20]. Due to the inherent characteristics of SAR image, such as speckle noise, target segmentation processing technology of SAR image is quite different from that of optical image [21]. In this paper, based on the characteristics of CSAR image, the target and background are segmented. The projection direction of the target in each sub-aperture image is perpendicular to the irradiation direction of SAR. Thus, each sub-aperture image can be viewed as the superposition of the inherent background and the aspect angle-dependent projection [22]. Sub-aperture images are reshaped to a vector and form a composite matrix. Sparse and low rank matrix decomposition techniques are used to separate the similar and unique parts of each sub-aperture image, that is, the background and foreground are separated. We keep the sparse matrix and separate the foreground target. Robust principal component analysis (RPCA) which is a robust and efficient method of solving the separation problem [23], is adopted in this paper.

The proposed method realizes the 3D reconstruction of the target without additional information, which can describe the characteristics of the target more precisely and make full use of the information contained in the angle dimension of CSAR. It transforms the inevitable image matching step in radargrammetry into the voting process, which is more efficiency.

The remainder of this letter is organized as follows. In Section 2.1, RPCA is explained. In Section 2.2, we describe the basic principle of the proposed method. In Section 3, the proposed method is applied on real Gotcha mission data. Furthermore, in Section 4, we give some discussions about the result.

2. Methodology

2.1. Robust Principal Component Analysis

The problem solved by low rank decomposition model is to decompose real data A and interference data E from observation data D. In practical application, the observation matrix in the model can be understood as the data we have got in reality. Through the decomposition algorithm, the corresponding real data and interference data can be obtained, i.e., $D=A+E$. The real data A need to satisfy the low rank constraint, while the interference data E needs to satisfy the sparse constraint. Using the similarity between sub-aperture images, each image is reshaped into a vector to form a composite matrix, which has a relatively low rank. In this way, the low rank part can be extracted from the constructed observation matrix. The low-rank matrix contains the targets with the same height, which equals to the reference imaging height, such as the inherent scene and speckle noise. The sparse matrix contains the targets with other heights, which is related to aspect angle of SAR. Therefore, the objects and changes in the image appear in the sparse part. Speckle noise is an inherent feature of SAR image and is evenly distributed. Thus, it has the same sparsity as the scene. The SAR image can be denoised at the same time. RPCA [24] is adopted in this paper because of its stability and effectiveness.

The recovery of matrix can be described by the following optimization problem:

$$\underset{A,E}{min}rank\left(A\right)+\lambda {\parallel E\parallel }_{0}subjectA+E=D$$

(1)

where $\lambda$ is the weight of the sparse part of the objective function. ${\parallel •\parallel }_0$ Represents a $l_0$-norm operation, which solves the number of non-zero elements in a matrix.

In the above formula, the rank function and $l_0$-norm are not convex, which is NP hard. The kernel norm of matrix is a reasonable relaxation of the solution of matrix rank function and $l_0$-norm is relaxation of $l_1$-norm.

For these relaxations, the final form of RPCA is derived as follows:

$$\underset{A,E}{min}{\parallel A\parallel }_{\ast }+\lambda {\parallel E\parallel }_{1}subjectA+E=D$$

(2)

Equation (2) can be solved by an augmented Lagrangian method [25]. Constructing Lagrange function as follows:

$$L(A,E,Y)={\parallel A\parallel }_{\ast }+\lambda {\parallel E\parallel }_1+<Y,D-A-E>$$

(3)

where Y is the Lagrange multiplier.

The constrained problem is transformed into an unconstrained problem, and a penalty term is added.

$$L(A,E,Y,\mu )={\parallel A\parallel }_{\ast }+{\lambda \parallel E\parallel }_1+<Y,D-A-E>+\frac{\mu }{2}·{\parallel D-A-E\parallel }_F^2$$

(4)

where $\mu$ is positive scalar quantity. Then coordinate descent method is applied to solve the Y, A, E and $\mu$ which satisfy the conditions. In each iteration, the extreme value is calculated along one coordinate axis direction (e.g., Y), while all other coordinate axes (e.g., A, E and $\mu$) are fixed. $\rho$ is the updated factor of $\mu$ in each iteration.

2.2. Inversely Mapping and Voting

In this paper, sub-aperture images of CSAR are formed on a reference ground plane by back-projection algorithm (BPA) [26]. In this way, all the sub-aperture images are formed in the same ground coordinate system.

Consider that a point $T(x_T,y_T,z_T)$ is above the reference ground plane where the height of the imaging plane is $z=z_p$. As is shown in Figure 1, if the elevation of the imaging plane is not consistent with the actual elevation of the target, the two dimensional coordinate of the target on the sub-aperture image will deviate from the actual position, resulting in geometric deformation of the full-aperture image [27].

Here we take viewing angle B for example to explain the imaging geometry of the target point, shown in Figure 2.

According to the imaging mechanism of SAR [28], which is described as follows:

$$\left\{\begin{array}{c}\hfill \overrightarrow{V_s}·\overrightarrow{ST}=0\\ \hfill \left|ST\right|=\left|SP\right|\end{array}\right.$$

(5)

where S denotes the zero doppler point corresponding to the viewing angle B. $\overrightarrow{V_S}$ is velocity of SAR platform. $\theta$ denotes the incidence angle. The projection point P is located at the intersection of zero doppler plane, the iso-range line L and imaging plane. Therefore when target is located on the plane which is perpendicular to the velocity direction of SAR platform, its doppler frequency is zero.

With Equation (5), we can project any point T from the 3D object space to the imaging plane along iso-range line L. In Figure 2, Line $LT$ is the intersection line between the imaging plane and the zero doppler plane. It indicates that the intersection line $LT$ is perpendicular to the instantaneous velocity of the center of synthetic aperture. The aspect angle $\phi$ refers to the angle between intersection line $LT$ and the positive direction of X-axis. Unless explicitly noted during the rest of this article, aspect angle is used to refer to different SAR imaging positions. The slope of intersection line $LT$ of a sub aperture image is expressed as follows:

$${\overrightarrow{k}}_{LT}\ast \overrightarrow{V_S}=0$$

(6)

Compared with the distance from the radar to the target, iso-range line L is approximately considered as a straight line $L^{\prime }$. Then the height of target $\Delta h$ can be obtained by Equation (7).

$$\Delta h=\Delta r\ast tan\theta$$

(7)

where $\Delta h$ is the height difference between the target and the reference ground plane. $\Delta r$ is the offset between the target and the projection point.

Above is the process of generating projection point from target. Then we will map the projection points into 3D space against the projection process step by step.

The calibration rules of candidate point: Firstly, along the intersection line $LT$, for any offset $\Delta x$ in x-direction, we can get the only $\Delta y$ in y-direction. The offset between imaging point and candidate point can be solved by Equation (8). Furthermore, the height value can be obtained by Equation (7). Finally corresponding to different offset, the projection point can be inversely mapped to a set of candidate points $(P_1,P_2,P_3,\cdots P_n)$ in 3D object space. The inversely mapping process of projection point P is shown in Figure 3. The projection point Q is inversely mapped along iso-range line to a set of candidate points $(Q_1,Q_2,Q_3,\cdots Q_n)$ under another aspect angle. Two iso-range lines intersect at $P_3$ and $Q_3$. The times of intersection of iso-range lines can be regarded as the process of voting for candidate points. Each candidate point in the 3D object space obtained from the inversely mapping process of the two-dimensional image is marked as “1” in the 3D object space. The inversely mapping process and voting process are carried out under multiple aspect angles. In Figure 3, for a candidate point, the more iso-range lines intersect, the higher the number of votes. Finally, the candidate point whose number of votes is greater than a certain threshold will be reserved, which is considered as a point of the target in 3D space.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \Delta r& =\sqrt{{(\Delta x)}^2+{(\Delta y)}^2}\hfill \end{array}\end{array}$$

(8)

where $\Delta x$ and $\Delta y$ is the component of the offset in the X-axis and Y-axis, respectively.

Then the procedure of 3D point cloud reconstruction method is described in detail. The algorithm flow chart is shown in Figure 4.

Step 1: The CSAR data is divided into several sub apertures. Back projection girds are utilized to form SAR images.

Step 2: The target and background of each sub-aperture image is separated by RPCA. Then the image is binarized. The projection part is set to 1, and the background part is set to 0.

Step 3: Under a certain aspect angle, the projection direction of the sub-aperture image is consistent. Under the constraint of the line of iso-range line, any point in the binarized “1” region is reversely positioned after given offset and projected into the 3D mesh. A series of candidate points are generated.

Step 4: Under a certain aspect angle, traverse all the points in the projection part.

Step 5: The sub-aperture images of all aspect angles are processed by Step 3 to Step 4.

Step 6: Counting the voting results of each point in the 3D grid. If the voting results are greater than a certain threshold F, the candidates are preserved.

3. Experiment and Results

In this section, to test the reconstruction performance of the proposed algorithm, the Gotcha volumetric SAR data is used. It is collected at the X-band with a 640-MHz bandwidth [29]. The spotlighted scene is located at a parking lot. In this experiment, we only use the eighth-pass data.

BPA is used to process the phase history data and form the raw sub-images sequence. Each sub-aperture covers 1${}^{\circ }$. The pixel spacing of the sub-aperture images is $0.2\times 0.2$ m. In this paper, the height of the imaging plane is 0 m.

Quantitative analyses of the reconstruction result and computational burden with comparative results of the conventional CSAR approach [16] using single pass CSAR data are demonstrated here. These experiments were conducted with the MATLAB R2018b software on a computer with an Intel Core 2.2-GHz processor and 64.0 GB of physical memory.

3.1. Results Obtained by the Proposed Method

Based on the aforementioned RPCA model, the targets and the background are separated from the raw sub-images sequence. The low-rank matrix contains the targets with the same height, which equals to the reference imaging height. The sparse matrix contains the targets with other heights. It is supposed that the sparse portion involved targets, whereas the low-rank portion indicated the scene and inherent noise. The full scene of entire aperture after RPCA is shown in Figure 5.

Then the entire aperture is divided into 36 sub-apertures and every 10${}^{\circ }$ sub-aperture images after RPCA are non-coherently accumulated. In total, 36 group images can be obtained through the way of dividing sub aperture in this experiment. Two parameters $\mu$ and $\rho$ are set to perform RPCA processing. In Equation (4), $\mu$ is the weight of penalty item and $\rho$ is the updated factor of $\mu$ in each iteration. $\mu =1/{\parallel D\parallel }_2$, $\rho =1.5$. The weight on the sparse term is set as $\lambda =1/max(m,n)$ for $m\times n$ matrix D.

Under a certain aspect angle, some strong scattering points of the target form projection points in the 2D plane. As is shown in Figure 6, it is difficult to separate the target from the scene because the image is affected by the speckle noise. The projection of the extracted target on the plane is flaky, so the iso-range lines generated by the projection points in the reverse mapping will be very dense. The structure of the final target is blocky and hard to be recognised. However, in Figure 7, the aspect-independent object can be easily obtained from the image.

Then the cars marked by yellow box, red box and blue box in Figure 7 are selected to restore 3D point cloud from the SAR image with the proposed method.

Some targets are anisotropic which means that they can only be illuminated under some aspect angles. So the threshold of voting number is set at 20 in this paper. Meanwhile the grid spacing in elevation direction is 0.2 m. The models of the three selected vehicles are different, which are shown in Figure 8a, Figure 9a and Figure 10a.

The 3D reconstruction results of the selected targets are shown in Figure 8b, Figure 9b and Figure 10b. Compared with the optical image of the car, 3D point cloud and the heading are restored successfully. The structures of the three kinds of vehicles can be distinguished and are consistent with the models shown in the optical image. The shape of the vehicles can also be seen in the result of 3D reconstruction. Furthermore, we can distinguish different types of vehicles according to the reconstructed image using only single pass CSAR data.

Table 1. Quantitative evaluation of vehicle C.

Dimensions	Length (m)	Width (m)	Height (m)
actual value	4.98	1.86	1.42
restored value	4.88	1.74	1.44
error	0.10	0.12	−0.02

,

Table 2. Quantitative evaluation of vehicle B.

Dimensions	Length (m)	Width (m)	Height (m)
actual value	4.75	1.74	1.41
restored value	4.60	1.60	1.40
error	0.15	0.14	0.01

and

Table 3. Quantitative evaluation of vehicle F.

Dimensions	Length (m)	Width (m)	Height (m)
actual value	4.45	1.77	1.44
restored value	4.40	1.80	1.60
error	0.05	−0.03	−0.16

show the actual vehicles’ size parameters and the estimated results obtained by using this method. It can be seen from the tables that compared with the actual parameters, the estimated geometric parameter error of the vehicle is less than 0.2 m, using only single pass CSAR data. On this basis, without additional information, a more detailed description of the target is achieved, where the shape and contour information of the target are provided.

3.2. Comparison with 3D Imaging Method Using Single Pass of CSAR Data

Traditional 3D Back projection(3D BP) imaging [10] performs interpolation operation and phase compensation operation for each azimuth pulse and grid point by grid point. The number of full aperture pulses is specified as $N_a$ and the size of 3D imaging grid of CSAR is specified as $N_x\times N_y\times N_z$. 3D reconstruction image is obtained by coherently accumulating the 3D imaging results of each azimuth pulse. The final result can be obtained by testing of generalized likelihood ratio (GLRT) [30] after 3D imaging of full aperture. The maximum value of the corresponding grid position in the three-dimensional imaging results of each sub aperture can be selected as the final imaging result at the grid.

3D imaging results of the three vehicles are shown in Figure 11. As mentioned before, the retrieved information is distorted by undesired sidelobes. It is hardly impossible to achieve 3D CSAR imaging using only one single pass CSAR data for anisotropic scatterers. Therefore the geometric size of the target can not be obtained from the 3D imaging results using one single pass CSAR data.

The following is the comparison of computational cost of the two methods:

For an $N_x\times N_y$ CSAR image, the size of 3D imaging grid of CSAR is specified as $N_x\times N_y\times N_z$. The operation time of one interpolation and phase compensation of BP algorithm is specified as $t_1$.

Traditional 3D BP algorithm: the traditional 3D imaging BP algorithm needs to do interpolation and phase compensation operation $N_x\times N_y\times N_z\times N_a$ times, and the time consumed by the 3D BP algorithm is:

$$T_1=N_x\times N_y\times N_z\times N_a\times t_1$$

(9)

The proposed method is to decompose 3D BP imaging algorithm into two parts: 2D BPA and restoring height information from 2D sub-aperture images. Therefore the time consumed by our algorithm needs to consider two parts: one is the 2D BP algorithm, and the other is the process of inversely mapping and voting process for the 3D mesh $N_x\times N_y\times N_z$ on the sub-aperture images. The full aperture is divided into $N_{sub}$ sub-apertures. The times of execution of the second part is $N_x\times N_y\times N_z\times N_{sub}$. The operation time of inversely mapping and voting process is specified as $t_2$. The time consumed by the proposed algorithm is:

$$T_2=N_x\times N_y\times N_a\times t_1+N_x\times N_y\times N_z\times N_{sub}\times t_2$$

(10)

Therefore, the time-consuming ratio of the two methods is

$$\begin{array}{cc}\hfill \frac{T_2}{T_1}& =\frac{N_x\times N_y\times N_a\times t_1+N_x\times N_y\times N_z\times N_{sub}\times t_2}{N_x\times N_y\times N_z\times N_a\times t_1}\hfill \\ \hfill & =\frac{N_x\times N_y\times N_a\times t_1+N_x\times N_y\times N_z\times N_{sub}\times t_2}{N_x\times N_y\times N_a\times t_1+N_x\times N_y\times (N_z-1)\times N_a\times t_1}\hfill \end{array}$$

(11)

where $t_2<t_1$, $N_{sub}\ll N_a$.

According to the above analysis, we transform the ($N_x\times N_y\times N_z\times N_a$) times interpolation and phase compensation operation along elevation direction into ($N_x\times N_y\times N_z\times N_{sub}$) times inversely mapping and voting process with two multiplication operations and one addition operation, which greatly reduces the amount of calculation. Hence, compared to 3D BP imaging method, the proposed method can realize the 3D point cloud reconstruction using single pass CSAR image with much higher efficiency.

The time-consuming comparison between the traditional 3D BP imaging of single pass CSAR image and the proposed method is shown in

Table 4. Time-consuming comparison between 3D BP imaging and the proposed method.

Time-Consuming	3D Imaging (s)	The Proposed Method (s)
vehicle C	8111	1023
vehicle B	7608	970
vehicle F	8135	1037

. It can be seen from the results that the time consumption of this method is much lower than that of 3D BP imaging method.

4. Discussion

The proposed method utilized the relationship between the projection point of circular SAR image sequence and the position of target in 3D space to reconstruct 3D point cloud of the target. The candidate points of target are located on the iso-range lines of each projection point. The intersection of iso-range lines can be regarded as a voting process. For a candidate point, the more times of intersection, the higher the number of votes, and the candidate point will be reserved.

According to the results of 3D reconstruction, not only the geometric parameters of the target can be obtained, such as length, width and height, but also different types of vehicles can be distinguished according to the reconstructed 3D point cloud. Under the constraint of projection direction, the inverse mapping of projection points is carried out in the projection area, which achieves more accurate and fast reconstruction. In the case of no additional information, such as polarimetric scattering information, single circular SAR data is used to realize 3D reconstruction. In this way, the information contained in the aspect angle dimension is mined to the greatest extent.

It can be seen that there is noise interference in SAR measured data in Figure 6, which makes it difficult to segment the target region. In this paper, RPCA is used to segment the target and background preliminarily which achieves the effect of denoising at the same time, shown in Figure 7. Then it is conducive to preserve the details of the subsequent three-dimensional reconstruction. Compared with the optical image and the actual geometric parameters of the target, the proposed method is effective.

However, we also noticed that the contour part is not complete when only single pass CSAR data is used for processing, shown in Figure 10b. Since the vehicle is close to the lawn, the projection information of the vehicle itself is affected, the wrong points appear in the 3D reconstruction point cloud of the vehicle.

In this paper, we only select independent vehicle targets with a gap for the experiment, so as to ensure that they are not affected by other vehicles. In the future research, we will consider the situation of high target density. In this case, the projection information between targets will overlap. Therefore, it may be necessary to add other types of auxiliary information to help expand the application scope of the proposed method.