Down-Looking Airborne Radar Imaging Performance: The Multi-Line and Multi-Frequency

: The paper proposes an analytical study regarding airborne radar imaging performances and accounts for a down-looking radar system moving along parallel lines far, in terms of probing wavelength, from the investigated domain and collecting multi-frequency and multi-monostatic data. The imaging problem is formulated in a constant depth plane by exploiting the Born approximation. Hence, a linear inverse scattering problem is faced by considering both the Adjoint and the Truncated Singular Value Decomposition reconstruction schemes. Analytical and simulated results are provided to state how the achievable performances depend on the measurement conﬁguration. These results are of practical usefulness because, in operative conditions, it is unfeasible to plan a ﬂight grid made up by a high number of closely (in terms of probing wavelength) spaced lines. Hence, the understanding of how the availability of under-sampled data affects the radar imaging allows for a trade-off between operative data collection constrains and reliable reconstructions of the scenario under test.


Introduction
Radar systems mounted on airborne platforms like small unmanned aerial vehicles (UAVs) are worth being considered because UAV allows for a not trivial simplification of the measurement logistic phase. UAV makes, indeed, possible an easy access to areas, which are hard to be reached by human operators or dangerous for them, and it permits the survey of wide areas without efforts. Accordingly, different kind of UAV-based downlooking radar systems have been recently proposed [1][2][3][4][5][6][7][8], someone of them working at high frequency and applied for landmine detection [2][3][4][5]7] and cultural heritage [8].
Whatever is the considered radar system and the application of interest, a common issue is the design of the measurement configuration as a function of the flight parameters, i.e., flight altitude, number of measurement lines, and number of measurement points along each line. The radar imaging result depends, indeed, on the amount and quality of data as well as on the adopted data processing strategy [9,10]. Moreover, when a UAV radar system is employed, one must account for that, in practice, it is unfeasible to plan a flight grid with a high number of lines spaced by a small, in terms of probing wavelength, spatial offset. Therefore, it is clear that an open practical issue is to understand how the imaging capabilities depend on the measurement setup.
A contribution to this topic has been given in [11], where analytical and numerical results have been presented for the case of single frequency, multi-lines, multi-monostatic data. Specifically, in [11] a radar moving along parallel lines at a certain altitude and collecting single frequency data at nadir has been considered. In addition, the imaging has been formulated in the 2D (x-y) plane under the Born approximation [12], and the Truncated Singular Values Decomposition (TSVD) approach [13] has been adopted as inversion procedure.
In this paper, the same measurement configuration and assumptions of [11] are made but, differently from [11], multi-frequency data are exploited. Moreover, the imaging capabilities of the Adjoint and the TSVD reconstruction approaches [13] are investigated. The analysis is performed by accounting for the Point Spread Function (PSF) as a tool to foresee the effect of the measurement configuration on the achievable spatial resolution along x and y directions. Analytical results are firstly derived. After, numerical results, referred to a wide-band high frequency radar system, are presented to show the benefits introduced by the availability of multi-frequency data and to compare the performance of the two considered inversion approaches.
The paper is organized as follows. Section 2 deals with the formulation of the imaging problem and recalls the inversion procedures. Section 3 provides analytical results stating how multi-frequency multi-monostatic well and under-sampled data affect the imaging capabilities in terms of spatial resolution and occurrence of grating lobes. Numerical results concerned with a point-like target are presented in Section 4, while the case of an extended target is tackled in Section 5. Conclusions are given in Section 6.

Imaging Problem
Let us consider the imaging scenario depicted in Figure 1. A radar system, mounted on a UAV platform, collects data over the planar surface M = (−a, a) × (−b, b) by traveling along n parallel measurement lines, with n = 1, . . . , N ny . The measurement lines are parallel to the x-axis, evenly spaced along the y-axis and at a constant altitude h above the imaging surface. This latter is the spatial domain D = (−a', a') × (−b', b') located at z = 0. The radar operates in monostatic mode and illuminates the scene at nadir (i.e., along z-axis) at each measurement point r = xx + yŷ + hẑ. The transmitting and receiving antennas are modeled as wide-band linear dipoles, oriented along the y-axis and operating at a central frequency belonging to the C-band. The targets are located into the imaging domain D. Under the above assumptions and considering the Born Approximation, the scattering phenomenon is described in the frequency domain by the following linear integral equation: In Equation (1), E S (r, ω) denotes the collected component of the backscattered field at the measurement point r and at the angular frequency ω; k 0 = ω √ µ 0 ε 0 is the free-space propagation constant (µ 0 and ε 0 being the free-space magnetic permeability and dielectric permittivity, respectively). Moreover, χ(r ) = ε t (r ) ε 0 − 1 is the unknown contrast function, ε t (r ) being the dielectric permittivity function at the generic point r = x x + y ŷ in D; G = f s (r, r , ω) is the dyadic Green function in free space; while E f s i (r, r , ω) is the incident filed, whose expression is: where I 0 l is the transmitting dipole current moment and Z 0 is the free-space impedance. The start and stop approximation is assumed, i.e., the radar is considered motionless (at a fixed measurement point r) in the time interval between the transmission of the probing signal, i.e., the incident field, and the reception of the corresponding backscattered field.
Let the radar operate in far-field [12], by substituting Equation (2) in Equation (1) and exploiting the far field expression of the Greens function, the scattering equation is rewritten as: where the dependency onŷ has been omitted.
is the linear scattering operator mapping the unknown space into data space, where L 2 denotes the space of the square integrable functions and Ω = [ω min , ω max ] is the radar frequency band, ω min and ω max being the minimum and maximum angular frequencies, respectively.
To face the imaging problem as defined by Equation (3), two reconstruction strategies are considered: the Adjoint procedure [13] and the Truncated Singular Value Decomposition (TSVD) inversion scheme [13]. These approaches allow for a closed form expression of the contrast function given in terms of the SVD of the discretized version of the operator L. Specifically, after discretization, the imaging is faced as the solution of the matrix inverse problem: where E S is the N-dimensional data vector, χ is the Q-dimensional unknown vector, Q being the number of points in D, and L is the N × Q-dimensional matrix obtained by discretizing the integral operator L. Let {u n , σ n , v n } be the singular system of the matrix L, we have [13]: Lv n = σ n u n and L + u n = σ n v n (5) where L + denotes the adjoint matrix, {σ n } are the singular values in decreasing order, while {v n } and {u n } define the orthonormal basis of the object space and the data space, respectively. By exploiting the SVD representation of L, the estimated contrast function is expressed as:χ according to the Adjoint procedure, and aŝ according to the TSVD approach. In Equations (6) and (7), ·, · is the scalar product; the parameter N in Equation (6) is the amount of collected data, while the parameter T in Equation (7) is the regularization threshold [13]. The modulus of the contrast function as defined by Equation (6) or Equation (7) is the tomographic image of the scenario under test.

Imaging Performance Analysis
This Section aims at providing analytical results relating the imaging capabilities to the measurement configuration. These results are propaedeutic to the outcomes of the numerical analysis presented in the following sections. In detail, with respect to the imaging scenario sketched in Figure 1, this Section discusses about the amount of independent data, i.e., the minimum number of samples required to represent the backscattered field properly, as well as the achievable spatial resolution limits and the occurrence of grating lobes.
The analysis deals with the far field condition and, as previously done in [11] for the single frequency case, it starts with the assumption that both the frequency range Ω and the measurement domain M are sampled continuously. Thereafter, it takes into account that a dense sampling is reasonable for Ω and the along-track direction, i.e., the x-axis for the scenario at hand (see Figure 1), whereas it is practically unfeasible for the across direction, i.e., the y-axis. The sampling measurement step along y-axis, indeed, must be compliant with the positioning accuracy reachable by the on-board Global Navigation Units (GNUs), which control the UAV flight trajectory. Consequently, a coarse measurement step is expected along the y-axis. The discrete sampling case is distinguished by considering: (i) data sampled in agreement with the Nyquist criterion (i.e., well sampled data) and (ii) under sampled data.

Continuous Sampling
Based on the far field condition, we assume k 0 h 1, h |x − x | and h |y − y |. In this case, by applying the paraxial approximation [14] the distance term R appearing in the exponential term of Equation (3) is expressed as: and Equation (3) is rewritten as: The kernel of the integral in Equation (9) is separable and we can tackle the dependence on the spatial variables independently. Let us start by considering the variable x and the equation: , Equation (10) is rewritten as: where k x = 2k 0 x h . Equation (11) is the Fourier Transform (FT) of the contrast function χ, which has the spatial limited support [−a , a ].
According to the Nyquist criterion, the sampling step of the FT of a spatial limited function in [−a , a ] is ∆ = π a [13]. On the other hand, it is easy to verify that in the multimonostatic and multi-frequency case herein at hand, the variable k x takes values into the range − 4πa λ min h , 4πa λ min h , where λ min is the minimum wavelength. Hence, it follows that the number of samples required to represent E S (x, ω) is given by [15,16]: By repeating the same reasoning for the spatial variable y', the minimum number of samples to represent E S (y, ω) is: Therefore, the total number of multi-monostatic and multi-frequency data required to represent the scattered field, while assuring no loss of information, is: It is worth pointing out that the Nyquist sample number, as given by Equation (14), defines the amount of not redundant samples necessary to represent the backscattered field. However, it does not provide a practical rule to properly set the number of measurement points on M, i.e., along x and y axes, as well as the number of frequencies in Ω. Indeed, although x, y and ω are independent variables, there is an infinite number of their combinations which satisfies the Nyquist criterion for sampling k x and k y . On the other hand, an evenly sampling of the measurement domain M, along x and y, as well as of the frequency range Ω, can provide an insufficient amount of independent samples of the backscattered field even if the gathered number of data is equal or higher than the Nyquist number.

Discrete Sampling: Resolution Limits and Grating Lobes
Let us assume that the backscattered field is measured at a finite number of angular frequencies ω i with i = 1, . . . , N f in the range Ω and at a finite number of evenly spaced points (x m , y n ) on M with m = 1, . . . , N mx and n = 1, . . . , N ny . Moreover, M x and M y denote the amount of multi-monostatic and multi-frequency data collected along the x-axis and the y-axis, respectively. As in the previous Section, we consider the dependence on the two spatial variables independently and start by considering the x variable. Hence, let

Nyquist Sampling Configuration
Let us consider M x = NDF 1D x defined in Equation (12). In this case, the samples E s,p are the Fourier harmonics of χ(x ) and the solution of the inverse problem in Equation (15) has the following expression: being In order to appreciate the spatial resolution, we use Equation (16) to calculate the Point Spread Function (PSF) referred to a point like target located a x = 0. By means of some matematical passages, which are not reported for sake of brevity, it is obtained the following expression [11]: The PSF in Equation (17) depends on x by means of the quantity B x , it is maximum at x = 0 and it has its first null is at: This means that the spatial resolution limit along the x-axis is dictated by the maximum frequency of the considered range and depends on the flight altitude h.
Note that the modulus of the PSF in Equation (17) has period π with respect to the variable B x 2 . Hence, grating lobes occurr at B x 2 = ±qπ with q ∈ N − {0}.
All the above considerations are valid also for y , provided that M x is replaced by M y , ∆k x by ∆k y = 8πb M y λ min h and B x by B y = − 8πb M y λ min h y . Hence, with respect to y , the PSF is given by while the 2D PSF is expressed as: It is worth pointing out that, while the spatial resolution limits along x and y do not depend on M x and M y , these latter affect the imaging capabilities due to the possible occurrence of grating lobes. In the case at hand, M x and M y are given by Equations (12) and (13), respectively, hence ∆k x,opt = π a and ∆k y,opt = π b while the parameters B x and B y get the values: Consequently, since grating lobes occurr when B x 2 = ±qπ and B y 2 = ±qπ with q ∈ N − {0}, they are at x = ±2qa and y = ±2qb , which are points outside the investigated domain. In other words, as it is obvious, the Nyquist sample number assures that no aliasing effects occurr.

Undersampling Configuration
Let us consider that the data are collected according to the Nyquist criterion along the UAV travelling direction, i.e., the x-axis, while under-sampled data are collected along the across direction, i.e., the y-axis. As said at the beginning of this Section, this is compliant with the fact that a large amount of measurement points is feasible along the x-axis, while practical flight constrains make unfeasible to consider a large number of dansly spaced measurement lines. Hence, in practice, an often encountered condition is the availabilty of an amount of data such that M y < NDF 1D y . Let us express the sampling step as ∆k y = π b , b being an auxiliary spatial parameter such that ∆k y = π b > ∆k y,opt = π b The amount of backscattered field samples is given by M y = 8bb λ min h . In this case, the resolution limits do not change because they do not depend on M y ; conversely, grating lobes occur at: and alsiaing issues arise.

Performance Analysis
This Section investigates how the imaging capabilities of the Adjoint procedure and the TSVD regularization scheme (see Section II) depend on the adopted measurement configuration. The analysis deals with spatial resolution limits and occurrence of grating lobes. In particular, a point-like target centered in (x t = 0, y t = 0), i.e., at the origin of the reference system given in Figure 1, is considered. The imaging domain, D, is a square located at z = 0 m whose side is 3 m long, i.e., a = b = 1.5 m. The data are gathered in the frequency range Ω = (3.5, 4.5) GHz by a UAV mounted radar. The latter covers a measurement domain M, which is a square with sides 3 m long, and it travels at height h along multiple measurement lines N ny parallel to the x-axis. Three examples are considered and are referred to as Example A, B and C.

Example A (Height of the Platform Equal to 15 m)
Let us assume that the flight altitude is h = 15 m. According to the above definition of D, M and Ω and the analytical results given in Section 3, we have NDF 1D x = NDF 1D y =18, the spatial resolution is expected to be 0.17 m along x and y, while grating lobes occur when an amount of not redundant data lower than NDF 2D = 18 × 18 = 324 is gathered.
However, the problem arises about the effective location of the measurement points on M and the number of frequencies in Ω. Under the assumption that the measurement domain M and the frequency range Ω are evenly sampled, we consider the five cases summarized in Table 1, which also provides the total amount of collected data, N, and the value of TSVD threshold, T. This latter is set in correspondence of the beginning of the fast decay of the singular values. For the considered cases, this means to filter out all the singular values whose normalized amplitude is lower than −10 dB, see Figure 2.  Figures 3a-c and 4a-c show the normalized modulus of the PSF obtained by using Equation (6) (Adjoint PSF) and Equation (7) (TSVD PSF), respectively. In addition, Figure 5a,b compares the PSF cuts along y directions. These figures corroborate that both the reconstruction approaches allow for comparable performance in Case 1 and Case 3, while grating lobes occur at the sides of the investigated domain in Case 2, even if they have negligible amplitude for the TSVD approach. Moreover, both the approaches reach a spatial resolution well approximating the expected one. In both the spatial directions, indeed, the spatial resolution achievable by Adjoint and TSVD are 0.18 m and 0.16 m, respectively.
Based on the provided results, we consider Case 3 as the best trade-off in terms of imaging performance and amount of data to be processed. In this case, indeed, both the approaches provide images not affected by grating lobes. Accordingly, it is reasonable of assuming that the measurement parameters of Case 3 allows us to collect the required amount of not redundant data. The other outcome of the present analysis is that, for the example at hand, the number of spatial measurement has a larger effect on the grating lobes issue.
In Case 4 and Case 5, the number of measurement points along x and of frequencies are fixed as in Case 3, while the number of measurement lines along y is progressively decreased in order to investigate how the undersampling expected in this direction affects the imaging results. Figure 6a-c is referred to Case 4 while Figure 7a-c accounts for Case 5. As expected, the spatial resolution does not change in both the spatial directions and grating lobes appearing along y. Moreover, whatever the inversion procedure, their number increases as the measurement lines decrease, and their location is well approximated by Equation (22) with M y = N vs NDF 1D x , N vs being the number of about constant singular values, i.e., the singular values whose normalized amplitude is not less than −6 dB.       Figures 6 and 7 show also that the TSVD inversion allows for improved reconstruction capabilities since the grating lobes have lower amplitude than in the Adjoint reconstruction. This result is compliant with those given in [16,17] and it is due to the not steplike decreasing behaviour of the singular values. It is also worth pointing out that the loss of data, due to reduced number of measurement lines, cannot be compensated for by increasing the number of frequencies. In particular, by setting = 5, the imaging results  Figures 6 and 7 show also that the TSVD inversion allows for improved reconstruction capabilities since the grating lobes have lower amplitude than in the Adjoint reconstruction. This result is compliant with those given in [16,17] and it is due to the not step-like decreasing behaviour of the singular values. It is also worth pointing out that the loss of data, due to reduced number of measurement lines, cannot be compensated for by increasing the number of frequencies. In particular, by setting N f = 5, the imaging results are similar to those given in Figures 6 and 7 and are not shown for the sake of brevity.

Example B (Height of the Platform Equal to 2.5 m)
Let us investigate what happens if the assumptions h |x − x | and h |y − y | are removed and assume the flight altitude h = 2.5 m. In this case, the PSF analytical expression given in Equation (20) is not valid in principle but, in practice, can be still considered to have indication about the expected performance and to have hints about the design of the measurement configuration. Specifically, in the case at hand, the spatial resolution expected on the basis of the analysis reported in Section 3 is 0.03 m along x and y, while grating lobes are foreseable when an amount of independent data lower than NDF 2D = 108 × 108 = 11664 are gathered.
Based on the previous results, we have considered the cases summarized in Table 2.  Figures 8-10 show the 2D PSF reconstructed by means the Adjoint and the TSVD procedures and their cuts along y for N mx = 63, N f = 3 and N ny = {63, 5, 3}. As for the Example A and the following ones, the TVSD threshold has been set in such a way to consider all the singular values before the fast decay. This implies that the threshold T is equal to the amount of available data N, when few measurement lines are accounted for. Figures 8-10 corroborate that, as discussed in Section 3, the spatial resolution improves by decreasing the flight altitude; it does not change when the amount of data decreases, and its value is well estimated by Equation (18). In fact, the spatial resolution is about 0.03 m in both the directions and for both the considered reconstruction approaches. Moreover, as for the flight altitude h = 15 m, when few measurement lines are considered the TSVD allows for a better reconstruction, in terms of amplitude of grating lobes, compared to the Adjoint procedure, even if the reconstruction provided by both the approaches are affected by artefacts occurring along y.   It is worth noting that now the artefacts can be reduced by setting N f = 5, as one can observe by Figure 11a,b, which show the Adjoint PSF and TSVD PSF for the Case 4 in Table 2, and comparing them with Figures 8c and 9c, respectively.

Example C
Let us consider the reconstruction of a point-like target when the spacing among the measurement lines y is fixed and equal to 0.6 m. The target is still located at the center of the investigated domain, which is the same of the previous examples, while the data are collected at the flight altitude h = 2.5 m. The measurement points along the x-axis are N mx = 63; the measurement lines are N ny = 5 and N ny = 3, the frequencies in Ω are N f = 5 (see Table 3).  Figure 12a,b show the Adjoint PSF and the TSVD PSF referred to Case 1, while those referred to Case 2 are given in Figure 13a,b.
It is worth noting that the resolutions along y are different w.r.t. the previous cases, because the size of the measurement domain along y changes. Indeed, for Case 1 the parameter b is 1.2 m, while for Case 2 it is b = 0.6 m. The resolution values along y axis are 0.35 m and 0.7 m for Case 1 and Case 2, respectively; they are compliant with the analytical values given by Equation (18). Moreover, concerning grating lobes, the TSVD approach achieves again better performance w.r.t. the Adjoint inversion procedure.
Indeed, Figures 12 and 13 show that the Adjoint PSFs have grating lobes and artefacts whose amplitude is larger compared to the TSVD PSFs.

Numerical Example
This section aims at assessing the performances of the two reconstruction approaches in the case of extended objects. In particular, we consider an extended object located on a ground having relative permittivity equal to 4. The object is a parallelepiped, having size 0.50 × 0.50 × 0.05 m 3 and relative permittivity equal to 6, and it is located about at the center of the investigated domain (the simulated scenario is shown in Figure 14). Specifically, the center of the object is placed at x = −0.10 m, y = 0 m and z = 0.025 m.
Time domain data have been simulated by means of gprMax 3D, which is a FDTD simulation software commonly used to simulate time-domain GPR data [18]. Five straight measurement lines are considered along the y-axis. They are 3 m long along the x-axis, 0.6 m spaced from y = −1.2 m to y = 1.2 m, and 2.5 m far from the air-soil interface. A dipole antenna oriented along the y-axis, whose central frequency is 4 GHz, is taken as primary source and only the y-component of the backscattered field is collected at N mx = 75 evenly spaced measurement points along each line. The simulated data have been transformed into the frequency domain by setting the useful frequency range equal to (3.5 − 4.5) GHz and considering N f = 5 frequencies evenly spaced of 0.5 GHz.  Figure 14 shows a sketch of the reference scenario, while Figure 15a,b show Adjoint and TSVD reconstructions in the plane y = 0 m, where the black dashed contour represents the true target. As for the examples reported in Section 4, the TSVD threshold is such to consider the singular values before the fast decay, i.e., T = 1795 (see Figure 15d). According to the analytical results, the TSVD approach allows for a better spatial resolution and reconstructs accurately shape and size of the object into the (x − y) plane. On the other hand, few ghosts appear and are due to the grating lobes occurring for the low number of measurement lines. In the TSVD approach, the amplitude of the grating lobes is lower than the amplitude of the object with the exception of those arising on the side of the investigated domain. Conversely, in the Adjoint-based reconstruction, the grating lobes have amplitude comparable to the object and this implies that the y-size of the object is significantly overestimated.

Discussion and Conclusions
The performance of UAV radar imaging systems depends on multiple system factors, i.e., radar operative parameters, measurement setup and data processing strategy. The paper has presented an analytical and numerical study aimed at assessing the reconstruction capabilities of down looking airborne radar systems able to collect multi-frequency data. The study regards the performance achievable when the imaging is faced under the Born Approximation, which is the usual assumption for facing a radar imaging problem, f.i. in applicative framework of Ground Penetrating Radar [9].
The first part of the presented analysis deals with data collected in far-field condition and satisfying the Nyquist criterion. Under these assumptions, a closed expression of the Point Spread Function (PSF) has been derived and considered as analytical tool to estimate the spatial resolution limits and investigate how the imaging capabilities depend on the amount of available data. It is worth pointing out that the analytical results presented in Section 3 do not depend on the adopted inversion approach and provide a general indication about the achievable imaging capabilities.
In addition, a numerical analysis has been provided to account for the reduced amount of data commonly available, and the far-field condition may not satisfied. The numerical analysis has been carried out by using two widespread inversion strategies, i.e., the Adjoint and the TSVD approaches. Moreover, first, it accounts for a single point-like target and, then, an extended object, whose scattered field has been simulated by using a full-wave electromagnetic software widely exploited by the GPR community. When the measurement configuration is in agreement with the theoretical hypotheses, the numerical results corroborate that the imaging results achieved by means of Adjoint and TSVD inversion approaches exactly match with the resolutions and the grating lobes appearance of the analytical PSF. Conversely, whereas the resolution limits do not change, the occurrence of grating lobes depends on the adopted inversion scheme and the presented results state that the TSVD approach allows for improved results compared to the Adjoint, especially when the flight altitude decreases. Herein, the Adjoint inversion has been considered because it is at the basis of the classical imaging methods (see [10]), while the TSVD inversion was been taken into account due to the fact that the singular values are characterized by a slow dynamic before the fast exponential decay [13]. Other regularization inverse schemes might be used especially when the singular values have a smooth decay [13]; this could be an interesting topic of future work.
As a final remark, we underline that the presented analysis accounts for data collected in ideal conditions, i.e., does not account for the uncertainties about the actual location of the measurement points, which is one of the main problems arising when consider experimental data. Accordingly, the manuscript provides indications about the imaging capabilities achievable by means of a UAV mounted radar system in the best case. As future research activity, we are planning to assemble a radar system on-board a UAV platform, design a measurement campaign in a laboratory-like controlled scenario and exploit the Carrier-phase differential global positioning system (CDGPS) technique [19] to accurately estimate the UAV positions during the data acquisition step. The final goal of this future activity will be the evaluation of the discrepancies between the results achieved by the theoretical analysis and those achieved by processing real datasets.

Conflicts of Interest:
The authors declare no conflict of interest.