P-Band UAV-SAR 4D Imaging: A Multi-Master Differential SAR Tomography Approach

Abstract

Due to its rapid deployment, high-flexibility, and high-accuracy advantages, the unmanned-aerial-vehicle (UAV)-based differential synthetic aperture radar (SAR) tomography (D-TomoSAR) technique presents an attractive approach for urban risk monitoring. With its sufficiently long spatial and temporal baselines, it offers elevation and velocity resolution beyond the dimensions of range and azimuth, enabling four-dimensional (4D) SAR imaging. In the case of P-band UAV-SAR, a long spatial-temporal baseline is necessary to achieve high enough elevation-velocity dimensional resolution. Although P-band UAV-SAR maintains temporal coherence, it still faces two issues due to the extended spatial baseline, i.e., low spatial coherence and high sidelobes. To tackle these problems, we introduce a multi-master (MM) D-TomoSAR approach, contributing three main points. Firstly, the traditional D-TomoSAR signal model is extended to a MM one, which improves the average coherence coefficient and the number of baselines (NOB) as well as suppresses sidelobes. Secondly, a baseline distribution optimization processing is proposed to equalize the spatial–temporal baseline distribution, achieve more uniform spectrum samplings, and reduce sidelobes. Thirdly, a clustering-based outlier elimination method is employed to ensure 4D imaging quality. The proposed method is effectively validated through computer simulation and P-band UAV-SAR experiment.

1. Introduction

Differential synthetic aperture radar (SAR) tomography (D-TomoSAR) is a rapidly evolving technique for four-dimensional (4D) SAR imaging [1]. It surpasses the range and azimuth dimensions by providing deformation velocity measurements while achieving elevation-dimensional imaging. In fact, it is a three-dimensional (3D) + one-dimensional (1D) imaging technology, where 1D represents the deformation information. With its high-resolution and high-accuracy advantages, it has vast application potential in urban management and risk monitoring [2].

D-TomoSAR is a progression of TomoSAR, which has established the foundation of multi-dimensional SAR imaging [3]. In prior research, airborne platforms greatly contributed to the feasibility verification of three-dimensional (3D) imaging [4] through SAR tomography (TomoSAR). The combination of differential interferometric SAR (D-InSAR) and TomoSAR, D-TomoSAR was introduced in [5], extending SAR imaging to 4D. Since the core of (D-)TomoSAR is spectrum estimation [6,7], spaceborne platforms have provided a wealth of data and propelled the study of sparse spectrum estimation algorithms [8,9,10,11,12]. Particularly, the introduction of compressed sensing (CS) theory has significantly enhanced imaging quality [13,14,15]. Moreover, imaging dimensions have even been extended to 5 or more [11,16,17]. Utilizing the advantages of large mapping width, spaceborne D-TomoSAR is a common solution for urban monitoring [18,19].

Due to the rapid deployment and high flexibility [20,21], the unmanned-aerial-vehicle (UAV)-based D-TomoSAR technique offers a new way for 4D urban mapping. In addition, benefiting from the advantage in penetration ability and temporal coherence maintaining [22,23], the P-band SAR is used in the vertical structure extraction of the forests [24,25], biomass mapping [26,27], moisture estimation of the soils [28], ground deformation surveying [29], etc. Therefore, it is worth exploring the application of UAV-based P-band D-TomoSAR. The benefits of such a system can be summarized as follows.

1.

Rapid emergency response. Compared to spaceborne SAR systems, which often have several-day or longer revisit periods, UAV-SAR systems have the advantage of being rapidly deployable due to their flexible trajectory planning capabilities [30].

2.

Baseline control. Compared with the uncontrolled baseline of the spaceborne SAR formed by orbit perturbation [31], the UAV much more easily controls their trajectories. The spatial baselines can be achieved from the utilization of the integrated navigation system (INS) and real-time kinematic (RTK) technique [32].

3.

High-quality SAR imaging. Benefiting from the long wavelength, the UAV-based P-band SAR is not sensitive to motion errors [33]. With the help of the INS and RTK techniques, the centimeter-level UAV trajectory measurement accuracy can easily meet the requirements of high-quality SAR imaging [34].

4.

Temporal coherence maintaining. Benefiting from the low working frequency, the UAV-based P-band SAR has obvious advantages in temporal coherence maintaining [22]. It creates favorable conditions for long-period urban monitoring.

5.

No atmospheric phase error. Because of its low flight altitude [35], it is not necessary to consider the atmospheric phase error problem in UAV-SAR, which is an unavoidable problem in the spaceborne case [36].

However, despite the promising capabilities of the UAV-based P-band D-TomoSAR, its application remains infrequent. The resolutions of elevation and deformation velocity are proportional to the wavelength and inverse to the lengths of spatial and temporal baselines. Due to the long wavelength of P-band SAR, long baselines are necessary to achieve better resolutions for elevation and deformation velocity measurements. This demand for a long baseline can introduce several problems, including:

1.

Low spatial coherence. The multi-pass working mode is usually adopted in D-TomoSAR to acquire the multi-baseline (MB) interferograms [2]. There will be both spatial and temporal baselines between the SAR images. Image coherence is mainly affected by two aspects: temporal and spatial coherence [22,37,38]. Although the P-band SAR is easier to maintain temporal coherence [22], it still faces spatial decoherence under long baselines (quantitatively analyzed in Section 2.2). In the traditional single-master (SM) D-TomoSAR model, it is difficult for the fixed master image to balance the spatial coherence with all other images. Even if the master trajectory is optimally located at the center of all trajectories, there will still be significant spatial decoherence of the far-away trajectories’ images. The average coherence coefficient (ACC)’s reduction of the interferograms in the D-TomoSAR data stack (DDS) will introduce serious phase noise and affect the D-TomoSAR processing.

2.

High sidelobes. The essence of D-TomoSAR is two-dimensional (2D) (spatial and temporal) sparse spectrum estimation [6,7]. Specifically, the distribution of the spatial-temporal baselines determines the sampling form of the 2D spectrum [39,40]. The number of baselines (NOB) is equivalent to the sample number, which will affect the spectrum estimation effect [7,41]. However, the limited UAV flights make the NOB insufficient. The sparse problem is inevitable, and strong sidelobes will occur. What is more, the long baseline makes the sparsity problem more prominent. With the influence of phase noise, the problem of high sidelobes will be further aggravated.

To address the above problems, a multi-master (MM) D-TomoSAR approach is proposed in this paper. The three main contributions of this approach are summarized as follows.

1.

Signal model extension. A fixed master image cannot balance the spatial coherences with all other images. It will reduce the ACC and NOB and cause high sidelobes. To address this issue, we propose a solution through the idea of MM D-TomoSAR using a non-fixed master image. Although the bistatic-like MM model has been introduced previously in [39,40,42,43], it is a limited one, and only bistatic-like SAR images can be paired to generate interferograms to avoid atmospheric phase errors. In contrast, the proposed model effectively utilizes the advantage that UAV-SAR is free from atmospheric phase errors. Any two SAR images can generate interferograms to construct the DDS. In this way, all images can be paired with their adjacent ones to obtain higher spatial coherence. The improved NOB and ACC (quantitatively analyzed in Section 3.2) will contribute to sparse spectrum estimation and suppress the sidelobes.

2.

Baseline distribution optimization. The essence of D-TomoSAR lies in 2D sparse spectrum estimation, and a more even sampling can be beneficial [7]. We extend the baseline sign reassignment process introduced in [39] to the 2D case. It is achieved through a step-by-step vector-accumulation-based baseline sign reassignment. The sign reassignment in each step aims to minimize the magnitude of the accumulated vectors. It can achieve a more uniform 2D spectrum sampling, facilitating the subsequent spectrum estimation.

3.

Outlier elimination. To further mitigate the impact of false estimations caused by sidelobes, an outlier elimination method is proposed in this paper. The rationale behind this method is the clustering behavior of the scatterers in a local region. Any estimated scatterers that do not belong to these clusters can be identified as outliers and removed. The clustering is accomplished through a 2D threshold-based discrimination process, which ensures the accuracy of elevation-velocity estimation.

Based on the above improvements, the proposed MM D-TomoSAR methodology can be implemented through three stages. In the first stage, the DDS can be constructed based on the MB interferograms utilizing the MM approach. The second stage entails optimizing the DDS through a step-by-step baseline sign reassignment. The third stage involves performing the elevation-velocity joint estimation. The outliers are eliminated by means of adaptive-threshold-based clustering.

The rest of this paper is organized as follows. In Section 2, the traditional SM D-TomoSAR approach with its problems is introduced. In Section 3, the proposed MM D-TomoSAR approach with its signal model, advantages, and processing flow is explained in detail. In Section 4, the simulation and P-band UAV-SAR experiments validate the proposed approach. The discussion of the proposed method is in Section 5. The conclusion and future work are described in Section 6.

2. SM D-TomoSAR

2.1. SM Signal Model

The implementation of D-TomoSAR is similar to that of MB-InSAR [44]. It attains the imaging in elevation-velocity dimensions by synthesizing spatial and temporal apertures. The temporal aperture is generated by the multi-pass observation mode. As illustrated in Figure 1a, multiple UAV flights at varying altitudes and times provide 2D apertures. Leveraging the capability of elevation resolution, layover points such as A and $A^{\prime }$ can be discerned, and their deformation velocities can be quantified.

Assuming N images are obtained through N UAV flights, the deformation of a certain target at an elevation s is commonly assumed to be a linear motion with a velocity of v. The value $g_n$ of an azimuth-range pixel for the nth SAR image is the integral of the reflected signal along the elevation and velocity dimensions, as described in [3]:

$$g_n=\underset{\Delta v}{\int }\underset{\Delta s}{\int }\gamma (s,v)exp(-j2\pi ({\xi }_{n}s+{\eta }_{n}v))dsdv$$

(1)

where $\Delta s$ and $\Delta v$ are the range of possible elevations and velocities. $\gamma (s,v)$ represents the reflectivity function along elevation and velocity. ${\xi }_n$ and ${\eta }_n$ are the elevation and velocity frequencies, respectively; ${\xi }_n=-2b_n/\left(\lambda r\right)$ and ${\eta }_n=2t_n/r$, where $\lambda$ and r are the wavelength and slant range, respectively; $b_n$ and $t_n$ are the spatial and temporal baselines between the nth and the master images.

Based on Equation (1), the measurement $g_n$ can be considered to be a two-dimensional Fourier transform (FT) of $\gamma (s,v)$. Thus, the estimation of elevation s and velocity v can be performed through spectrum analysis. It should be noted that the estimated elevations in this context refer to the normal slant-range (NSR) direction that is determined by the view angle. To obtain the target’s actual height (perpendicular to the ground), a conversion through the view angle is required. Similarly, the measured deformation velocity is in the line-of-sight (LOS) direction. The measurements corresponding to different spatial–temporal baselines can be seen as discrete samples of the spectrum (as depicted in Figure 1b). In the presence of noise $\mathbf{\epsilon }$, Equation (1) after discretization can be written as [3]:

$$\mathit{g}=\mathit{G}·\mathbf{\gamma }+\mathbf{\epsilon }$$

(2)

where $\mathit{g}={[g_1,g_2,\cdots ,g_N]}^T$ is the measurements vector with N elements. $\mathit{G}$ is an $N\times K$ mapping matrix with $G_{n,k}=exp[-j2\pi ({\xi }_n{s}_k+{\eta }_n{v}_k)]$. $\mathbf{\gamma }$ is the discrete elevation-velocity reflectivity matrix after vectorization with $K=P\times Q$ elements, and $\mathbf{\gamma }={[\gamma (s_1,v_1),\cdots ,\gamma (s_1,v_Q),\cdots ,\gamma (s_P,v_1),\cdots ,\gamma (s_P,v_Q)]}^T$.

In D-TomoSAR, the estimation of $\mathbf{\gamma }$ is crucial. As depicted in Figure 1b, the 2D spectrum has limited samples due to insufficient N, resulting in sparse problems during spectrum analysis. Various algorithms can be adopted to address sparse spectrum estimation, such as the nonparametric spectrum analysis class represented by the translated singular value decomposition (TSVD) algorithm [45,46,47,48] and the CS class represented by the iterative shrinkage thresholding algorithm (ISTA) [48,49,50,51].

2.2. The Problems of the SM Signal Model

As for its performance of D-TomoSAR, the elevation resolution ${\rho }_s$ and velocity resolution ${\rho }_v$ are both determined by the wavelength $\lambda$ and baselines, and ${\rho }_s=\lambda r/(2\Delta b)$ and ${\rho }_v=\lambda /(2\Delta t)$. $\Delta b$ and $\Delta t$ are the maximum spatial (as shown in Figure 1a) and temporal baselines, respectively. As for the distribution of the 2D baselines, the selection of the master image is crucial. All the SAR images should be registered to the master one for interference (as the orange arrows in Figure 1a) and generate the interferograms to construct the DDS.

For P-band SAR, due to the long wavelength, the need for a long baseline is prominent. For clarity, Figure 2a shows the required baselines of the P-band SAR under different elevation resolution demands. The UAV is set to 150 m altitude with 65${}^{\circ }$ view angle. The working frequency is 400 MHz. It can be seen that it requires ∼100 m baseline to obtain the elevation resolution within 1.5 m. At this time, the baseline length is even comparable to the flight altitude. Under such a long baseline, there will be two problems that occur in P-band UAV-SAR for differential tomography processing.

One is the low spatial coherence introduced by a long baseline. Although the P-band SAR has the advantage of maintaining temporal coherence [22], it still faces spatial decoherence in long-baseline cases. Figure 2b shows the spatial coherence coefficients curve under different baselines. The range resolution and terrain slope are set to 0.5m and 15°, respectively. The other parameters are the same as in Figure 2a. It can be seen that the spatial coherence coefficient descends significantly with the growth of baseline length. The coherence coefficient of 100 m baseline can drop to ∼0.7. Considering the impact of noise, the image coherence will be much lower. It can impact the quality of the interferograms which are used to construct the DDS. Additionally, the fixed master image in the SM model (described in Section 3.2) can not balance the spatial coherences with all the other images. So, the ACC will decline so as to bring in phase noise to the DDS. It will also reduce the estimation accuracy of the elevation and deformation velocity.

Another is the high sidelobe introduced by the limited NOB. Since the essence of D-TomoSAR is 2D spectrum estimation [6,7], the spectrum sampling form will affect the estimation effect. With the limited NOB due to insufficient UAV flights, the samples in the spatial–temporal spectrum are also sparse. So, the problem of high sidelobes will occur unavoidably. What is more, this problem will be more prominent under the long baseline. The above problems will all further impact the 4D imaging effect.

3. MM D-TomoSAR

3.1. MM Signal Model

As discussed earlier, the SM signal model with a fixed master image is limited in long-baseline cases. Previous research [39,40,42,43] has proposed a bistatic-like MM signal model for atmospheric phase error elimination, which can be adopted for image pairing in D-TomoSAR. This model allows for an unfixed master image. Building upon this idea and leveraging the fact that UAV-SAR is free from the atmospheric phase, we propose an extended MM signal model (abbreviated as MM model henceforth).

The schematic of the proposed MM model is shown in Figure 3a. Herein, any SAR image can be used as the master one and paired with others. The image of any spatial-location trajectory can generate interferograms with high spatial coherence by pairing them with adjacent ones. Through this processing, among the N images, there are up to $C_N^2$ interferograms can be extracted. The NOB can also be improved. This phenomenon is equivalent to the amplification of the samples in the two-dimensional spectrum (as exemplified in Figure 3b), thereby boosting the processing of spectrum estimation. Based on the SM model in Section 2.1, the proposed model is given as:

$$\mathbf{\psi }=\mathsf{\Psi }·\mathbf{\gamma }+\mathbf{\epsilon }$$

(3)

where $\mathbf{\psi }$ is the measurement vector with $C_N^2$ elements, which corresponds to the $C_N^2$ interferograms. $\mathsf{\Psi }$ is a $C_N^2\times K$ mapping matrix, and $\mathbf{\gamma }$ is the discrete reflectivity vector with K elements that need to be estimated. $\mathbf{\epsilon }$ still represents the noise.

Among the above matrices (vectors), $\mathbf{\psi }$ and $\mathsf{\Psi }$ both consist of multiple submatrices, and each submatrix corresponds to each master image, i.e,

$$\left\{\begin{array}{c}\mathbf{\psi }={[{\mathbf{\psi }}_1,\cdots ,{\mathbf{\psi }}_i,\cdots ,{\mathbf{\psi }}_{N-1}]}^T\hfill \\ {\mathbf{\psi }}_i=[{\psi }_{i_{i+1}},\cdots ,{\psi }_{i_j},\cdots ,{\psi }_{i_N}]\hfill \end{array}\right.$$

(4)

$$\left\{\begin{array}{c}\mathsf{\Psi }={[{\mathsf{\Psi }}_1,\cdots ,{\mathsf{\Psi }}_i,\cdots ,{\mathsf{\Psi }}_{N-1}]}^T\hfill \\ {\mathsf{\Psi }}_i=\left[\begin{array}{ccccc}{\mathsf{\Psi }}_{i_{i+1},1}& \cdots & {\mathsf{\Psi }}_{i_{i+1},k}& \cdots & {\mathsf{\Psi }}_{i_{i+1},K}\\ ⋮& \ddots & ⋮& \ddots & ⋮\\ {\mathsf{\Psi }}_{i_j,1}& \cdots & {\mathsf{\Psi }}_{i_j,k}& \cdots & {\mathsf{\Psi }}_{i_j,K}\\ ⋮& \ddots & ⋮& \ddots & ⋮\\ {\mathsf{\Psi }}_{i_N,1}& \cdots & {\mathsf{\Psi }}_{i_N,k}& \cdots & {\mathsf{\Psi }}_{i_N,K}\end{array}\right]\hfill \end{array}\right.$$

(5)

where ${\psi }_{i_j}=S_i·{S_j}^{\ast }·{\varphi }_{flat,i_j}^{\ast }$. $S_i$ is an azimuth-range pixel in the ith SAR image. $S_j$ represents the jth one. ${[·]}^{\ast }$ stands for the conjugate operation, and ${\varphi }_{flat,i_j}$ is the flat phase.

In the mapping matrix $\mathsf{\Psi }$, its elements ${\mathsf{\Psi }}_{i_j,k}$ are mainly determined by elevation and velocity frequencies. Different from the fixed view angle of the master image in the SM model, the view angles and the NSR directions vary with the master images here. So, a certain target will have different elevations after being projected to the respective NSR direction of different master images. Similarly, the LOS-direction deformation velocities of different master images are also different. Consequently, there will be a problem in the combination of the submatrices (${\mathsf{\Psi }}_i$) and the subvectors (${\mathbf{\psi }}_i$).

To tackle this problem, a solution based on scale zooming is proposed. The elevations of the targets that are perpendicular to the ground (abbreviated as ‘sky-direction elevations’) are fixed with respect to a certain reference ground. The elevation in the NSR direction can be converted from the sky-direction elevation through the view angles. Therefore, the different estimated elevations in respective NSR directions can be unified through the conversion. It is essentially a view-angle-based scale zooming process. Similarly, the different LOS-direction deformation velocities can also be unified through scale zooming. The deformation velocity in the sky-direction is converted from the LOS direction. For brevity, the sky-direction elevation and velocity are both abbreviated as “elevation” and “velocity” in later, respectively, and there is

$$\left\{\begin{array}{cc}\hfill {\displaystyle {\mathsf{\Psi }}_{i_j,k}}& {\displaystyle =exp(-j2\pi ({\xi }_{i_j}{h}_k+{\eta }_{i_j}{z}_k))}\hfill \\ \hfill {\displaystyle {\xi }_{i_j}}& {\displaystyle =\frac{-2b_{i_j}}{\lambda r_{i}sin{\theta }_i}}\hfill \\ \hfill {\displaystyle {\eta }_{i_j}}& {\displaystyle =\frac{2t_{i_j}cos{\theta }_i}{\lambda }}\hfill \end{array}\right.$$

(6)

where $h_k$ and $z_k$ are the kth elevation and velocity in the discrete reflectivity vector. $b_{i_j}$ and $t_{i_j}$ represent the spatial and temporal baseline between the ith and jth images. For the ith master image, its slant range $r_i$ and view angle ${\theta }_i$ are adopted in the calculation of the baseline parameters. and $1/sin{\theta }_i$ and $cos\theta$ are the scale zooming factors of the elevation and velocity, respectively. Corresponding to the above scale zooming processing, the reflectivity vector $\mathbf{\gamma }$ in Equation (3) is adjusted to $\mathbf{\gamma }=[\gamma (h_1,z_1),\cdots ,\gamma (h_1,z_Q),\cdots ,\gamma (h_P,z_1),\cdots ,$$\gamma (h_P,z_Q){]}^T$; and $K=P\times Q$.

In addition, the proposed MM model does not change the essence of D-TomoSAR. It is still a sparse spectrum estimation problem. Only the spectrum sampling form has been changed. The commonly used algorithms can all be adopted to solve this model.

3.2. The Advantages of the MM Signal Model

Although the MM model does not change the essence of D-TomoSAR as a sparse spectrum estimation, it still brings some benefits to processing. The two principal benefits of the proposed model are summarized as follows.

One is the increase in the NOB. In the SM model, $N-1$ interferograms can be extracted from N SAR images, but the MM model can generate $C_N^2$ interferograms. Figure 4a illustrates the comparison of the NOB between the two models for different numbers of images. As the number of images increases, the MM model’s advantage in the NOB becomes more significant. A larger NOB leads to better elevation estimation accuracy [7,41]. Furthermore, as D-TomoSAR involves sparse spectrum estimation, the increase in the NOB is equivalent to increasing the spectrum samples (as shown in Figure 1b and Figure 3b). It is beneficial to spectrum estimation by alleviating the sparsity problem.

Another is the improvement of the ACC. The ACC can indicate the degree of phase noise and impact the effectiveness of D-TomoSAR processing. Although the proposed MM model cannot enhance the spatial coherence of the long-baseline interferograms, it can increase the number of short baselines with higher coherence. This augmentation of the shorter baselines can improve the ACC. To clarify, consider the example of short baselines defined as the baselines between adjacent images. In the SM model, the reference master image is limited to only two short baselines. However, in the proposed MM model, $N-1$ short baselines can be generated from the N images. Assuming that the short baselines between images are evenly distributed within a 100 m range, using the coherence coefficient analysis results (Figure 2b), the ACC of all the baselines based on different models can be obtained in Figure 4b. It is apparent that the MM model’s ACC (the red curve) is superior to that of the SM model (the blue curve). This advantage becomes more prominent as N increases. It will contribute to the elevation-velocity joint estimation.

3.3. MM D-TomoSAR Processing Flow

3.3.1. Overview

Based on the proposed MM model, an MM D-TomoSAR processing approach is introduced. Aiming at obtaining high-accuracy elevations and velocities, the proposed approach adds the processes of baseline distribution optimization and outliers elimination. The processing flow of the proposed method is summarized as follows (shown in Figure 5).

Stage-1 involves DDS construction. Image registration is performed to extract the MB interferograms through traversal pairing. The DDS is composed of these interferograms. Stage-2 involves the DDS optimization. Baseline distribution optimization is performed through sign reassignment. It can achieve a more balanced distribution of the baselines and help to uniformly sample the spectrum. According to the optimization results, the spatial and temporal baseline parameters in the DDS are adjusted. Stage-3 involves imaging effect optimization. Elevation-velocity joint estimation is performed through spectrum analysis, and a clustering-based outliers elimination is adopted to guarantee the precision of the estimated elevations and velocities. This process is accomplished using a straightforward 2D threshold-based discrimination.

3.3.2. DDS Construction

The DDS is composed of several interferograms of different spatial and temporal baselines. In the MM model, the interferograms are generated from the registered images by traversal pairing. So, for interferogram extraction, image registration is the premise. Because the MM model increases the NOB and the amount of the interferograms significantly, pair-by-pair image registration will obviously reduce the processing efficiency. For example, among the N, there are $C_N^2$ interference pairs that can be obtained. According to the commonly used image registration algorithms based on similarity measurement [52,53,54,55], the largely added operations of sample window sliding will inevitably increase its time cost. With the increase in N, the cost will grow further.

Aiming at this problem, a simple solution is introduced here. This is, the registration parameters are extracted through only one fixed reference image, and registration parameters of other image pairs are obtained through recursion. In this way, among N images, only $N-1$ sets of registration parameters (i.e., pixel offsets) need to be calculated based on similarity measurement. The rest $C_N^2-(N-1)$ sets of pixel offsets (POs) can be obtained through recursion. Because the processing of recursion is more efficient than sample window sliding, this solution can maximize the registration efficiency under the MM model. In this way, the paired images in the MM model can all be registered with each other. The MB interferograms can also be extracted based on the registered images. Then, the DDS can be constructed based on these interferograms after flat phase removal.

In addition, the specific processing flow of the MM image registration is described in Appendix A.1 and will not be further detailed here.

3.3.3. DDS Optimization

On the basis of the constructed DDS, spectrum estimation can be performed for elevation and velocity estimation, but the baseline sign reassignment is added here before the spectrum estimation for DDS optimization. This processing aims at balancing the distribution of the spatial and temporal baselines. It will equalize the 2D spectrum sampling and contribute to the subsequent elevation-velocity joint estimation.

Based on the sign reassignment processing introduced in [39], the baseline sign reassignment is extended to the 2D case in this paper. The proposed method converts the sign reassignment into a vector accumulation processing under the constraint of modulus minimization.

In this method, the 2D baselines obtained in the MM model are distinguished by the dimension. Then, two sets of baseline sequences are normalized by the respective maximum value. The elements in the normalized sequences are paired to compose the 2D baseline vectors. These vectors can be regarded as the plane vector in the interval of $[-1,1]$. When the modulus of the accumulated vectors reaches the minimum, the distribution of the baseline reaches the most uniform. Assuming that $F_i$ is the assigned sign of the baseline vector ${\mathit{B}}_{i_j}$, the objective function for baseline sign adjustment is

$$\underset{F_{i_j}\in \{-1,+1\}}{\mathrm{minimize}}{∥\sum _{n=1}^{C_N^2}{F}_{i_j}·{\mathit{B}}_{i_j}∥}_2$$

(7)

Since $F_i$ has only two values of 1 and −1, the baseline sign reassignment can be achieved by exhaustive search in the case of the small NOB, but under the MM model, the sharply improved NOB increases the computation amount, and exhaustive search is no longer applicable.

To solve the third problem, a concise reassignment method based on step-by-step adjusting is introduced here. As shown in Figure 6, it can be performed as follows. First, among the sorted baseline vectors, a sign (1 or −1) for the 1st baseline is set as the starting. Second, the sign of the second baseline is adjusted to minimize the modulus of the accumulated vector of first and second. Then, continue the previous processing to the rest baselines step by step, and each step should minimize the current accumulated vector’s modulus (as Equation (7)). When the last step adjustment is reached, the baseline sign reassignment is also completed. In this way, a more uniform 2D spectrum sampling can be achieved. It will contribute to the subsequent elevation-velocity joint estimation.

Additionally, for the normalization and construction process of the 2D baseline vectors, refer to Appendix A.2, and there will be no further detailed explanation here.

3.3.4. Imaging Effect Optimization

Based on the optimized DDS, the signs of the obtained InSAR phases are adjusted according to the baseline signs. Correspondingly, the elevation and velocity frequency parameters are also adjusted as

$${\psi }_{i_j}=exp[j{F}_{i_j}·\angle (S_i·{S_j}^{\ast }·{\varphi }_{flat,i_j}^{\ast })]$$

(8)

$${\mathrm{\Psi }}_{i_j,k}=exp[-j{F}_{i_j}·2\pi (\frac{-2b_{i_j}}{\lambda r_{i}sin{\theta }_i}{h}_k+\frac{2t_{i_j}cos{\theta }_i}{\lambda }{z}_k)]$$

(9)

where ${\psi }_{i_j}$ and ${\mathrm{\Psi }}_{i_j,k}$ are the elements of $\mathbf{\psi }$ and $\mathsf{\Psi }$ in Equation (3). $\angle (·)$ represents the angle-taking operation. $F_{i_j}$ is the sign flag. When the assigned baseline sign is consistent with the original one, $F_{i_j}$ is 1. Otherwise, it is −1.

As described in Section 3.1, the essence of the proposed MM model (as Equation (3)) is still spectrum estimation. Based on the estimated reflectivity vector $\mathbf{\gamma }$, the layover points with their elevations and deformation velocities can be obtained. The commonly used algorithms can all be directly adopted for elevation-velocity joint estimation, such as the TSVD [45,46,47,48], ISTA [48,49,50,51], etc. The MM model does not increase the difficulty of model solving. Instead, benefiting from the increased number of samples and the more uniform sampling, the sidelobes will be effectively suppressed. The specific estimation processing will not be explained here in detail.

Although the MM model has the advantage of sidelobe suppression, there will still inevitably be sidelobes in the estimated results due to the influence of phase noise and insufficient sampling. To optimizate the imaging effect and obtain high-accuracy elevation and deformation velocity, an outliers elimination method is proposed here.

This method is based on the similarity of the local regions. The observation geometries and the building structures in a small local region are similar. So, for a small pixel window, the elevations of the scattering points in all pixels (whether there is a layover or not) shall have the clustering characteristics. Similarly, the deformation velocities also have similarities in local regions, and there is rarely mutation of a single pixel. For example, if there is only the overlap between the building facade and the ground in a local pixel window, the estimated points can almost be classified into two clusters. One represents the ground, which has lower elevations and a type of deformation (such as land subsidence). The other represents the building facade, which has higher elevations, and another type of deformation (such as thermal expansion). In this case, the estimated points that do not belong to these two clusters can be determined as outliers.

Herein, the clustering is realized through a simple way, i.e., the 2D threshold-based discrimination. Specifically, all estimated points in a local window can make up a pixel set. For the central pixel in this window, the estimated points with certain elevations and velocities can be identified one by one. These elevations and velocities of the non-outlier points should be close to those of other points in this set. If not, this point can be identified as an outlier and be eliminated.

As shown in Figure 7, the outlier elimination can be performed as follows. Firstly, a pixel window is set in the image. Except for the center pixel, all the estimated elevations and deformation velocities of the other pixels can make up the elevation and velocity point sets. Similarly, the estimated points in the central pixel can also compose two point sets. Secondly, an index as the number of intra-cluster points (NICP) is set here. It is calculated based on these point sets through the 2D threshold-based discrimination, and the threshold can be determined according to the observation geometries and the noise level. Third, the discrimination of outliers can be performed based on the obtained NICPs. The scattering points with small NICPs will be set as outliers and eliminated from the central pixel.

Through the above way, the outliers in the D-TomoSAR imaging results can be eliminated, avoiding the false points introduced by sidelobes and ensuring the accuracy of the estimated elevations and deformation velocities.

In addition, the specific processing flow with the formulas of the outlier elimination is described in Appendix A.3. We will not go into more detail here.

4. Experiment

To verify the proposed MM D-TomoSAR approach, the computer simulation and UAV-SAR experiment are both adopted. Among them, the simulation is based on point targets to verify the sidelobe suppression effect of the MM model. The P-band UAV-SAR experiment lasts for 5 days, which is used to verify the actual 4D imaging effect of the proposed method.

4.1. Computer Simulation

The simulation parameters are shown in

Table 1. The system parameters of the simulation.

Parameter	Value	Unit
Working Frequency	400	MHz
View Angle	65	°
Flight Number	26	-
Flight Altitude	100–200	m
Altitude Interval	4	m
SNR	5	dB

. There are 26 UAV flights evenly distributed within the altitudes of 100–200 m to form the 100 m spatial baseline, but the distribution of the observation time is uneven and within the time span of 15 days (360 h). Under such long enough baselines, a high enough elevation (within 1.5 m) and velocity (about 1 mm/h) resolution can be achieved.

As for the simulated targets, there are three sets of points has been adopted. Each set contains two layover points with the same scattering amplitudes but different elevations and deformation velocities. The elevation and velocity differences between the two points change with sets to provide a full comparison. In the simulation, the two points with the same elevations but different velocities, different elevations but the same velocities, and different elevations and velocities have all been adopted. The specific parameters of the two points are shown in

Table 2. The parameters of the simulated targets.

	Set-1		Set-2		Set-3
	Point-1	Point-2	Point-1	Point-2	Point-1	Point-2
Elevation	0 m	5 m	0 m	0 m	0 m	5 m
Velocity	0 mm/h	0 mm/h	0 mm/h	10 mm/h	0 mm/h	10 mm/h

. In it, the elevation and velocity of one point (denoted by ‘Point-1’) are fixed to 0, but the elevation and velocity of another one (denoted by ‘Point-2’) changes. The elevation difference is set to 5 m, and the velocity difference is set to 10 mm/h. All the simulations are conducted under the SNR of 5 dB.

For comparison, there are four methods that have been adopted, including two types of signal models (the SM and MM model) with two elevation-velocity joint estimation algorithms (TSVD and ISTA). The adopted signal models and algorithms of the four methods are shown in

Table 3. The comparison methods used in the computer simulation (and the UAV-SAR experiment).

	Method-1	Method-2	Method-3	Method-4
Signal Model	SM	SM	MM	MM
Algorithm	TSVD	ISTA	TSVD	ISTA

.

Figure 8 shows the comparison of the spatial and temporal baselines of the two types of models, whereas Figure 8a shows the distribution of the simulated UAV flights in the time-altitude coordinates. Each circle represents a flight. Figure 8b,c are the obtained spatial and temporal baselines of the SM and MM models, respectively. It can be seen that the NOB of the MM model has an obvious advantage (raised from 25 to 325). Additionally, the red circles in Figure 8c exactly represent the baselines in Figure 8b. It makes the contrast clearer.

Figure 9 shows the comparison of the estimated elevation-velocity planes. These planes have all been normalized based on the respective maximum peak values. These peaks represent the possible targets. The greater the peak value, the higher the possibility that there is a target at this elevation-velocity coordinate, where the red and white boxes represent the real positions of Point-1 and -2, respectively. The centers of these boxes are exactly the actual coordinates. Based on these estimated planes, the comparisons of the two aspects are as follows.

As for the comparison in the estimation algorithms, Figure 9a (Method-1) and Figure 9b (Method-2) both adopted the SM model, but the ISTA-based Method-2 has obviously lower sidelobes. The comparison of Figure 9c (Method-3) and Figure 9d (Method-4) shows a similar effect which proves the sidelobe suppression advantage of the ISTA-based methods. In addition to the above comparison of Set-1 simulation, the results of Set-2 (Figure 9e–h) and -3 (Figure 9i–l) provide a consistent conclusion.

As for the comparison of the signal models, Figure 9a (Method-1) and Figure 9c (Method-3) both adopted the TSVD algorithm, but the MM model-based Method-3 has lower sidelobes. Most of the sidelobes in Figure 9a are effectively suppressed in Figure 9c. The comparison of Figure 9b (Method-2) and Figure 9d (Method-4) shows a similar effect. Benefiting from the ISTA algorithm, the sidelobes in Figure 9b are much lower than Figure 9a, but the MM model-based Method-4 further suppresses the remaining sidelobes in Figure 9d. Similarly, the comparison of Set-2 (Figure 9e–h) and -3 (Figure 9i–l) provide a consistent conclusion. It can be seen that the sidelobe suppression advantage of the MM model is obvious whether the CS-class (ISTA) algorithm is adopted or not. Even with the traditional TSVD algorithm, the MM model can still obtain the estimated planes with lower sidelobes.

To quantitatively evaluate the sidelobe suppression effect, an index of mainlobe energy percentage (MEP) is adopted here. it represents the ratio of the mainlobe’s energy to the sum of all the lobes’ energies (whether they are main- or the sidelobes), and this index is provided as a percentage. The energy of the mainlobe is calculated by the integral of all the strengths from the mainlobe’s peak to the first trough. The higher the strength of the sidelobes, the lower the MEP, and otherwise, the opposite. So, it can be used as an effective index for the sidelobe suppression effect evaluation.

The evaluated MEPs of the 3 simulation sets with different methods are shown in

Table 4. The mainlobe energy percentage (MEP) of different methods in the simulation.

Methods	Method-1	Method-2	Method-3	Method-4
Simulation Set-1	5.25%	71.77%	7.96%	100%
Simulation Set-2	5.63%	75.08%	9.29%	97.23%
Simulation Set-3	5.87%	45.56%	9.32%	90.83%

. For the comparison of the estimation algorithms, the advantage of the ISTA-based methods (Method-2 and -4) in sidelobe suppression are prominent. Compared to the TSVD-based methods, the MEP has an improvement of approximately 10 times. For the comparison of the signal models, the proposed MM-model-based methods (Method-3 and -4) have significant sidelobe suppression ability. The MM model can improve the MEP by 1.5∼2 times in the comparison of the same estimation algorithms. So, the proposed method has the advantage of sidelobe suppression. It is consistent with the previous analysis.

Additionally, the elevation and velocity estimation accuracy have been evaluated (as shown in

Table 5. The evaluated elevation and deformation velocity accuracies in the simulation.

	Method-1	Method-2	Method-3	Method-4
Elevation Accuracy	0.35 m	0.17 m	0.28 m	0.17 m
Velocity Accuracy	0.47 mm/h	0.33 mm/h	0.35 mm/h	0 mm/h

) to compare the elevation-velocity joint estimation effect. It should be noted that the sidelobe suppression advantage of the MM model has been proved. The false points caused by the sidelobes are not involved in this evaluation. The elevation and velocities of the mainlobes (in the red and white boxes) close to the actual coordinates are selected as the estimated results. Taking the actual elevation and velocities as the references, the root square mean errors (RSMEs) of the estimated results can been calculated. The average RSMEs of the three sets of simulations are taken as the evaluated accuracies. From the comparison, it can be seen that under the same estimation algorithms, the MM-model-based methods have better accuracy than those of the SM model. Additionally, because the baseline is long and the Rayleigh resolution is high enough, the super-resolution advantage of the ISTA algorithm is not obvious. The layover targets can all be distinguished, but the ISTA-based methods also have greater advantages in the estimation accuracy than the others. So, Method-4 has the best accuracy among these methods.

4.2. P-Band UAV-SAR Experiment

4.2.1. Experimental Condition

In addition to the computer simulation, a P-band UAV-SAR experiment has also been conducted for verification. The main parameters of the SAR system are shown in

Table 6. The main parameters of the P-band UAV-SAR system in this experment.

Parameter	Value	Unit
Working Frequency	400	MHz
Resoultion	$0.5\left(R\right)\times 0.5\left(A\right)$	m
Polarization	HH	–
View Angle	65	°
Flight Number	19	–
Flight Altitude	90∼180	m
Altitude Interval	5	m

. This experiment is conducted in Huang Songyu township, Pinggu district, Beijing. The terrain of the observation area gradually rises from southwest to northeast. The local digital elevation model (DEM) of this region is shown in Figure 10a, which is acquired by TanDEM-X. It is under the 90 m grid with ∼2 m elevation accuracy [56]. Although the sparse DEM grid loses the terrain details, it can still be used as a reference for the local terrain trends.

An eight-rotor UAV (as Figure 10b) was taken as the platform. Based on the high-precision INS and RTK, it can accurately cruise according to the preset flight routes with trajectory control accuracy within 0.5 m. The measurement accuracy of the UAV trajectories is better than 3 cm, which can meet the requirements of P-band SAR imaging [21]. As for the designed UAV routes, there are 19 UAV flights evenly distributed within the altitudes of 90–180 m with the direction from west to east. The real UAV tracks are shown in Figure 10c. The whole experiment lasted for 5 days, and 4 flights were conducted every day except the last one. The specific observation time was unevenly distributed within the 5 days (as shown in Figure 10d). Additionally, it should be noted that in order to facilitate the 2D spectrum estimation processing, the correlation of the spatial–temporal baseline is minimized as much as possible. So, the flight altitude is not positively related to the observation time, and the time–altitude coordinates of each flight are shown in Figure 11a.

The center of the observation area is a small rotorcraft airport (as shown in Figure 10e). Two reflectors are arranged in the open area of the airport to provide the evaluation reference. One reflector (denoted by Reflector-1, as shown in Figure 10f and the red circle in Figure 10e) has the ability of deformation. It can move along the triaxial sliding tracks. The displacement control accuracy can reach 0.2 mm. During the experiment, according to the observation time (Figure 10d) of each flight, Reflector-1 adjusts its elevation-direction position. Assuming $D_i$ is the elevation-direction displacement of the ith flight, there is $D_i=V_d\times T_i$, where $T_i$ is the observation time of the ith flight, and $V_d$ is the deformation velocity. Here, it is set to 5 mm/h. Considering the long wavelength of the P-band SAR, the size (edge length) of the reflector is large and up to 1.3 m. In order to support such a heavy reflector, its base is also large and made of metal, which will also reflect the radar waves. The elevation difference between the reflector and the base is 1.5 m. Another reflector (denoted by Reflector-2, as shown in Figure 10g and the blue circle in Figure 10e) is fixed on the ground. It does not move at all during the 5 days., and the size of it is consistent with Reflector-1 to provide a comparison.

There are 19 SAR images that have been acquired through the back projection algorithm (BPA) [33]. The 10th image is shown in Figure 10h. Corresponding to the red and blue circles in Figure 10h, the local amplifications of Reflector-1 and -2 are shown in Figure 10i and Figure 10j, respectively. It can be seen that these two reflectors are both well focused. It is worth noting that the imaging result in Figure 10i is an accumulation result of two targets: Reflector-1 and its base. Due to the large size and metal material of this base, it will also reflect the radar waves. The slant range difference between Reflector-1 and the base relative to the UAV is within the range resolution (0.5 m). They will overlap in the SAR image. So, they can be regarded as two layover targets: one is fixed (the base), and the other one has deformation (Reflector-1). Although it is difficult to distinguish the two targets in 2D SAR images, this problem will be effectively solved through the subsequent D-TomoSAR processing (explained later).

4.2.2. Experimental Results

Similar to the previous computer simulation, four methods have been adopted for comparison. They include two types of signal models (the SM and MM model) with two elevation-velocity joint estimation algorithms (TSVD and ISTA). The four methods are consistent with Table 3.

Figure 11a shows the time–altitude information of the 19 UAV flights. It can be seen that there is no obvious correlation between the altitudes and observation times, which can be equivalent to the random sampling of the 2D spectrum. Further, the temporal–spatial baselines under SM and MM models are shown in Figure 11b and Figure 11c, respectively, where the red circles in Figure 11c exactly represent the 2D baselines in Figure 11b. It is clear that the NOB of the MM model rapidly increased (from 18 to 171). This is equivalent to a 9.5-time increase in the number of spectrum samples.

Based on the different (SM and MM) models, two sets of DDSs have been constructed. In the SM-based DDS, the 10th image is selected as the master one, and all the other images are registered to it. So, 18 interferograms with their coherence coefficient maps (CCMs) can be obtained. For the MM-based DDS, the master image is unfixed. There are 171 interferograms with their CCMs. To compare the quality of the DDSs of different models, their average CCMs have been obtained as Figure 12. They are calculated through the average of all (18 or 171) CCMs under respective models. From the comparison between Figure 12a,b, it can be seen that the coherence coefficient in the MM-model-based CCM is higher the that of the SM model. Especially for regions (such as the black boxes) with flat terrains, this advantage is more obvious. The increase in the average coherence coefficient reflects the improvement of the DDS quality, which will contribute to the subsequent D-TomoSAR processing.

After the DDS construction and the 2D baseline distribution optimization, the elevation-velocity joint estimation is performed. It should be explained that some low-quality pixels have been deleted here. The determination of these pixels is based on the image amplitude and the CCMs, which aims at avoiding their impact on the overall imaging effect.

Figure 13 shows the estimated elevations and deformation velocities based on the four methods, where Figure 13a,e,i,m are the estimated elevations of Method-1, -2, -3 and -4. Figure 13c,g,k,o are the corresponding deformation velocities. It should be noted that these 2D figures only show the elevations and deformation velocities of the strongest scattering points estimated in each pixel unit. From the comparison, it can be seen that there is no obvious difference in the overall terrain trends among the four methods. The estimated elevations rise slightly from southwest to northeast, which is consistent with the reference DEM (Figure 10a). As for the estimated deformation velocities, they are all almost close to 0 (except for Reflector-1). It is because the geologic structure of the observation area where the airport is located is relatively stable, and no obvious deformation occurs in the limited experiment duration, but Reflector-1 can also provide a reliable reference for the verification of deformation velocity estimation.

For a detailed comparison of the local regions, Figure 13b,f,j,n show the local magnifications of the black boxes. This region is farmland with flat terrain. It can be seen that there is obvious burr noise in the estimated elevations in Figure 13b (Method-1) and Figure 13f (Method-2), which are both based on the SM model. Similar to elevation estimation results, the estimated deformation velocities in Figure 13d (Method-1) and Figure 13h (Method-2) both have burr noise. Although the noise level of the ISTA-based Method-2 is slightly lower than that of the TSVD-based Method-1, it is still difficult to avoid the impact on the estimation accuracy, but the MM model-based methods show less noise in Figure 13j,l,n,p. Compared with the SM model, it has the obvious advantage in noise suppression. Moreover, the MM-model-based Method-4 has slightly less noise, benefiting from the sidelobe suppression ability of the ISTA algorithm.

To further compare the four methods, the estimated elevation-velocity planes of Reflector-1 are shown in Figure 14a–d. It can be seen that two peaks appear here. Considering that the Reflector-1 and its base are overlapping here, the two peaks represent the two targets, respectively. Based on the known elevation-velocity information, two reference boxes are set here. Where the red box represents the base with the center coordinate of (0 m, 0 mm/h). The white box represents Reflector-1 with the center coordinate of (1.5 m, 5 mm/h). On the whole, the estimated peaks are almost close to the actual target positions, but there are differences in the strengths of the sidelobes. Through the comparison between Figure 14a (Method-1) and Figure 14c (Method-3), most of the sidelobes in Figure 14a have been suppressed in Figure 14c, which is due to the use of the MM model. Benefiting from the ISTA algorithm, the amount of sidelobes has decreased in Figure 14b (Method-2) and Figure 14d (Method-4), and the strengths of these sidelobes have also declined. Through the comparison between Figure 14c (Method-3) and Figure 14d (Method-4), Method-4 shows a slight advantage due to the ISTA algorithm. Additionally, it should be explained that the contrast effect here is not as prominent as that in the previous computer simulation. The reason is that the large-size reflector has a high SNR, which ensures the estimation effect of all methods.

Figure 14e–h show the estimated elevation-velocity planes of Reflector-2. For this reflector, it is fixed on the ground. So, the center coordinate of the reference (red) boxes is fixed to (0 m, 0 mm/h). For such a stable reference target, all the four methods have shown better estimation performance than those of Reflector-1. On the whole, the estimated peaks are almost close to the actual positions. Except for Figure 14e (Method-1), most of the sidelobes of other methods have been effectively suppressed. Combined with the previous analysis of other regions, the MM model-based methods are more reliable. Among them, the ISTA-based Method-4 has the most robust estimation performance.

To quantitatively evaluate the sidelobe suppression effect, MEP is also adopted here. The evaluated MEPs of the two reflectors with different methods are shown in

Table 7. The evaluated MEP of different methods based on the reflectors in the UAV-SAR experiment.

Methods	Method-1	Method-2	Method-3	Method-4
Reflector-1	5.81%	81.50%	9.32%	100%
Reflector-2	5.42%	100%	9.19%	100%

. For the comparison of Reflector-1, the sidelobes are obvious due to the layover problem. The proposed MM model improves the MEP by 1.60 and 1.23 times of the TSVD- and ISTA-based methods, respectively. For the comparison of Reflector-2, all the four methods have shown better estimation performance than those of Reflector-1. By comparing Method-1 and -3, the proposed MM model has improved the MEP by 1.69 times. The above comparison shows the advantages of the MM-model-based methods in sidelobe suppression. This conclusion is consistent with the previous analysis.

Additionally, the elevation and velocity estimation accuracy of the four methods has been evaluated. The reference data for the elevation accuracy evaluation is the TerraSAR-DEM, as shown in Figure 10a. Although the grid of the reference DEM is sparse, it can still be used as a reference for local terrain trends. In order to avoid the impact of insufficient terrain details, only the flat regions in the black boxes in Figure 12) are selected for evaluation. Moreover, the relative elevation accuracy is selected as the evaluation result. It can avoid the impact of the absolute elevation error of the reference DEM. As for the evaluation of the deformation velocity, the average value of each set of estimated deformation velocities has been counted. They are all almost close to 0. Therefore, taking 0 as the reference, the RSMEs of the estimated deformation velocities have been calculated. It can reflect the dispersion degree of the velocities, which is equivalent to the noise level. The evaluation results are shown in

Table 8. The evaluated elevation and deformation velocity accuracies in the UAV-SAR experiment.

	Method-1	Method-2	Method-3	Method-4
Elevation Accuracy	3.25 m	2.97 m	1.52 m	1.45 m
Velocity Accuracy	0.81 mm/h	0.69 mm/h	0.19 mm/h	0.17 mm/h

. It can be seen that the MM-model-based methods have an obvious advantage in elevation-velocity joint estimation. Benefiting from the ISTA algorithm, Method-4 has the highest estimation accuracy.

Based on the estimated elevations and velocities of Method-4, the D-TomoSAR point cloud has been obtained. The 3D point cloud is shown in Figure 15a, where the color represents the elevation. The 4D point cloud colored by the deformation velocities is shown in Figure 15b. The color of the whole image is relatively consistent (except for Reflector-1), which represents that this region is relatively stable during the experiment. Figure 15c,d are the fusion images of the point cloud and the optical image in 2D view. It can be seen that the main features of the two types of images match well. It also proves the accuracy of estimated elevations.

5. Discussion

Based on the above experimental results and analysis, we summarize the advantages of the proposed MM-model-based D-TomoSAR approach. Its potential shortcomings and solutions are also given as follows.

5.1. The Advantages

The advantages of the proposed method can be summarized as two aspects: sidelobe suppression and noise suppression.

Sidelobe suppression is in terms of the elevation-velocity joint estimation of a single pixel. The proposed method achieve sidelobe suppression via two main points. On the one hand, the proposed MM model can increase the NOB through the traverse combination of the SAR images. Because that D-TomoSAR processing is essentially a spectrum estimation problem in the elevation-velocity dimensions, the NOB can be equivalent to the number of spectrum samples. The increased NOB (and spectrum samples) can contribute to spectrum estimation. On the other hand, the proposed method also adopted the DDS optimization through baseline sign reassignment. It can make the sampling more balanced and also benefit sidelobes suppression. So, for the estimated elevation-velocity plans, the proposed method usually has lower sidelobes (as the estimated elevation-velocity planes in Figure 9 and Figure 14).

Noise suppression is in terms of the imaging effect of the whole image. The proposed method reduces the noise via three main points. Firstly, the proposed method has the advantage of sidelobe suppression (as mentioned above). So, it can reduce the false points in the elevation-velocity joint estimation, which often manifest as noise. Secondly, the proposed method constructs the DDS based on the MB interferograms through a traverse combination of the SAR images. In this way, the number of short baselines will significantly improve. Short-baseline image pairs often have higher spatial coherence and lower phase noise. So, the noise in the DDS and imaging result can also been reduced. Thirdly, the proposed method adopted a clustering-based outlier elimination. It can eliminate the outliers in the imaging results and also reduce the noise. So, the estimated elevations and deformation velocities usually have less noise (as shown in Figure 13).

5.2. The Shortcomings

The processing efficiency of the proposed method may be lower than that of the traditional SM-model-based methods. The limited processing efficiency is mainly shown in two aspects: image registration and elevation-velocity joint estimation.

For image registration, the proposed method adopts a traversal combination way for MB interferograms extraction. This approach can increase the NOB, alleviate the sparsity issues, and contribute to sidelobe suppression. However, improving the NOB requires more sets of image pairs to be registered. It will increase the computation amount of image registration processing.

For the elevation-velocity joint estimation, the proposed MM model required a larger size of matrices in the model solving. The size of the matrices is influenced by the NOB. A larger NOB corresponds to a larger matrix size. So, the MM model involves large-scale matrix operations, which increases the computational amount and reduces the processing efficiency.

5.3. The Solutions for the Shortcomings

The solutions for the above shortcoming introduced by the image registration and elevation-velocity joint estimation have also been provided here.

For image registration, we have provided a solution via the recursion-based image registration strategy (as Section 3.3.2). That is, the registration parameters are extracted through only one fixed reference image, and registration parameters of other image pairs are obtained through recursion. This strategy significantly reduces the window sliding processing in image registration and replaces them with recursive processing. Because the recursion processing is more efficient than window sliding, this solution can maximize the registration efficiency under the MM model.

For the elevation-velocity joint estimation, hardware acceleration based on the graphics processing unit (GPU) can be adopted. The spectrum estimation in D-TomoSAR mainly involves matrix operations. The MM model only changes the size and elements of the matrices, without changing the processing steps of the estimation algorithms. So, accelerating the efficiency of matrix solving can reduce time costs. GPU can efficiently handle matrix processing [57]. The time cost introduced by the large-size matrices of the MM model can be reduced in this way.

6. Conclusions

This paper proposes an MM D-TomoSAR approach which is suitable for P-band UAV-SAR. It solves the problems of low spatial coherence and high sidelobes introduced by long baseline via three main improvements. First, the traditional D-TomoSAR signal model is extended to an MM one. The rapidly increased NOB relieves the sparse problem in spectrum estimation. The increased number of short baselines improves the ACC and the quality of the DDS. Second, a baseline distribution optimization processing is proposed, which is achieved by a step-by-step sign reassignment processing of the baseline vectors. It can equalize the distribution of the spatial–temporal baselines to achieve more uniform spectrum samplings. Thirdly, to ensure the accuracy of the elevation-velocity joint estimation, a clustering-based outlier elimination method is introduced. It is realized through adaptive-threshold-based discrimination, which can effectively detect and eliminate outliers. The proposed approach is validated through computer simulation and P-band UAV-SAR experiment, which shows its advantages in sidelobe suppression and high-accuracy 4D imaging.

Additionally, it should be explained that no obvious deformation has been found in the UAV-SAR experiment (except for the reflector). It is because of the limited experimental duration and the stable observation scene. Therefore, future work can be divided into three aspects. First, the UAV-SAR observation shall be carried out in areas prone to deformation (such as mountains, mines, and bridges) for a longer time to verify its actual application effect. Second, combining the advantages of persistent scatterer (PS)-InSAR [58] that can eliminate the atmospheric phase errors, this method can be extended to the spaceborne mode to achieve large-width 4D imaging. Thirdly, due to high sidelobes being the common issue in D-TomoSAR, the proposed method also has application potential in other working bands. Future work will also extend this method in other bands and conduct validations.