Modeling measurements for quantitative imaging subsurface targets ^††thanks: This work was supported by the AFOSR Grant FA9550-24-1-0191. A. D. Kim also acknowledges support by NSF Grant DMS-1840265.

In two dimensions, the air-soil interface is $z=h(x)$. Let $x_{p}=-L/2+(p-1)L/P$ for $p=1,\dots,P$ denote a grid of $P$ points that uniformly samples the interval $-L/2\leq x<L/2$ with meshwidth $L/P$. We then use the points $\mathbf{r}^{\text{int}}_{p}=(x_{p},h(x_{p}))$ for $p=1,\dots,P$ to sample the air-soil interface. At $(x_{p},h(x_{p}))$ we use as the unit normal to this air-soil interface, $\hat{\mathbf{n}}_{p}=(-h^{\prime}(x_{p}),1)/\sqrt{1+(h^{\prime}(x_{p}))^{2}}$. Note that this choice always points into Region 0.

We assume the target boundary $\partial\Omega$ is a closed curve that can be parameterized according to $(\xi(t),\zeta(t))$ for $0\leq t\leq 2\pi$ with $(\xi(2\pi),\zeta(2\pi))=(\xi(0),\zeta(0))$. Let $t_{q}=2\pi(q-1)/Q$ for $q=1,\dots,Q$ denote a grid of $Q$ points that uniformly samples $0\leq t<2\pi$ with meshwidth $2\pi/Q$. We then use the points $\mathbf{r}_{q}^{\text{tar}}=(\xi(t_{q}),\zeta(t_{q}))$ for $q=1,\dots,Q$ to sample the target boundary. At $(\xi(t_{q}),\zeta(t_{q}))$ we use as the unit normal to this target boundary, $\hat{\boldsymbol{\nu}}_{q}=(\zeta^{\prime}(t_{q}),-\xi^{\prime}(t_{q}))/\sqrt{(\xi^{\prime}(t_{q}))^{2}+(\zeta(t_{q}))^{2}}$. Note that this choice always points into Region 1.

We introduce two, user-defined parameters: $\delta^{\text{int}}>0$ and $\delta^{\text{tar}}>0$. The MFS approximation for scattered field $u_{0}$ in Region 0 is then

with $G_{0}(\mathbf{r})=\mathrm{i}\frac{1}{4}H_{0}^{(1)}(k_{0}|\mathbf{r}|)$ denoting the free-space Green’s function with wavenumber $k_{0}$ and $\mathbf{r}$ denoting any point in Region 0. Note that the source positions: $\mathbf{r}^{\text{int}}_{p}-\delta^{\text{int}}\hat{\mathbf{z}}$ for $p=1,\dots,P$ lie outside of Region 0 (inside Region 1), so (4) exactly solves (1a) in Region 0. The expansion coefficients $a_{p}$ for $p=1,\dots,P$ are to be determined.

The MFS approximation for scattered field $u_{1}$ in Region 1 is given by

with $G_{1}(\mathbf{r})=\mathrm{i}\frac{1}{4}H_{0}^{(1)}(k_{1}|\mathbf{r}|)$ denoting the free-space Green’s function with wavenumber $k_{1}$ and $\mathbf{r}$ denoting any point in Region 1. Note that the source points $\mathbf{r}^{\text{int}}_{p}+\delta^{\text{int}}\hat{\mathbf{z}}$ for $p=1,\dots,P$ and $\mathbf{r}^{\text{tar}}_{q}-\delta^{\text{tar}}\hat{\boldsymbol{\nu}}_{q}$ for $q=1,\dots,Q$ lie outside of Region 1 (inside of Region 0 and Region 2, respectively), so (5) exactly solves (1b) in Region 1. The expansion coefficients $b_{p}$ for $p=1,\dots,P$ and $c^{\text{ext}}_{q}$ for $q=1,\dots,Q$ are to be determined.

Finally, the MFS approximation for the field interior to the target $u_{2}$ in Region 2 is

with $G_{2}(\mathbf{r})=\mathrm{i}\frac{1}{4}H_{0}^{(1)}(k_{2}|\mathbf{r}|)$ denoting the free-space Green’s function with wavenumber $k_{2}$ and $\mathbf{r}$ denoting any point in Region 2. Since the source points $\mathbf{r}^{\text{tar}}_{q}+\delta^{\text{tar}}\hat{\boldsymbol{\nu}}_{q}$ for $q=1,\dots,Q$ lie outside of Region 2 (inside Region 1), (6) exactly solves (1c) in Region 2. The expansion coefficients $c^{\text{int}}_{q}$ for $q=1,\dots,Q$ are to be determined.

MFS approximations (4), (5), and (6) are all just superpositions of free-space Green’s functions whose source positions lie outside of the regions where they are evaluated. Additionally, each of these free-space Green’s functions satisfy the correct outgoing conditions. Therefore, (4), (5), and (6) satisfy the appropriate outgoing conditions in their respective regions.

In (4), (5), and (6), there are $2P+2Q$ undetermined coefficients: $a_{p}$ and $b_{p}$ for $p=1,\dots,P$, and $c_{q}^{\text{ext}}$ and $c_{q}^{\text{int}}$ for $q=1,\dots,Q$. To determine these undetermined expansion coefficients, we apply a collocation method for boundary conditions (2) and (3).

For boundary conditions (2) we require that

and

for $p=1,\dots,P$. Note that $\partial_{n}=\hat{\mathbf{n}}_{p}\cdot\nabla$ in the notation above. Equations (7) and (8) give $2P$ conditions.

For boundary conditions (3), we require

and

for $q=1,\dots,Q$. Note that $\partial_{\nu}=\hat{\boldsymbol{\nu}}_{q}\cdot\nabla$ in the notation above. Equations (9) and (10) give $2Q$ equations.

Let $A_{1}$, $A_{2}$, $A_{3}$ and $A_{4}$ be $P\times P$ matrices whose entries are given by

Let $B_{1}$ and $B_{2}$ be the $P\times Q$ matrices whose entries are given by

Let $C_{1}$ and $C_{2}$ be the $Q\times P$ matrices whose entries are given by

Finally, let $S_{1}$, $S_{2}$, $S_{3}$, and $S_{4}$ denote the $Q\times Q$ matrices whose entries are given by

Because $\delta^{\text{int}}\neq 0$ and $\delta^{\text{tar}}\neq 0$, these matrix entries never correspond to evaluation of the free-space Green’s function where it is singular.

Using these matrices, we find that the $2P+2Q$ equations given by (7), (8), (9), and (10) as the following block linear system,

Here, $\mathbf{a}$ and $\mathbf{b}$ are the $P$-vectors whose components are $a_{p}$ and $b_{p}$ for $p=1,\dots,P$, respectively, and $\mathbf{c}^{\text{ext}}$ and $\mathbf{c}^{\text{int}}$ are the $Q$-vectors whose components are $c^{\text{ext}}_{q}$ and $c^{\text{int}}_{q}$ for $q=1,\dots,Q$, respectively. The right-hand side is made up of $P$-vectors $\mathbf{y}_{1}$ and $\mathbf{y}_{2}$ whose components are given by

for $p=1,\dots,P$.

Block system (15) can be solved numerically leading to the determination of the coefficients in (4), (5), and (6). With those expansion coefficients determined, those approximations can be evaluated anywhere within the regions to which they are defined. Here, we show an example computation.

In Fig. 2 we show MFS results for the fields scattered due to a point source located at $\mathbf{r}_{\text{src}}=(-25,75)$ cm, at frequency $f=5.1$ GHz and with $c=30$ cm GHz. For soil, we have used $\epsilon_{r1}=9$ and for the target we have used $\epsilon_{r2}=5$.

The air-soil interface is one realization of a Gaussian-correlated random rough surface with RMS height $h_{\text{RMS}}=0.4$ cm and correlation length $\ell=8$ cm generated using the method described by Tsang et al [11]. In Fig. 2, we have used $L=400$ cm and $P=512$ for the air-soil interface. For the target, we consider the “kite” whose boundary is given by $\mathbf{r}^{\text{tar}}(t)=(3,-14)+(3.0\cos t+1.8\cos 2t-0.65,3.4\sin t)$ and $Q=128$. We have set $\delta^{\text{int}}=\delta^{\text{tar}}=0.1$.

The example shown in Fig. 2 illustrates that the MFS approximation captures the complex interactions of the fields with the rough air-soil interface and the target. These complex interactions involve a variety of spatial scales and the computed solution smoothly transitions over these scales. Because the MFS is an easy to implement method and accurately accounts for the underlying physics of multiple scattering by the rough air-soil interface and the target, we use it for simulating measurements.

GP-SAR measurements are multi-frequency signals taken over several locations along the flight path. Let $\omega_{m}$ for $m=1,\dots,M$ denote the $M$ frequencies used to sample the system bandwidth and let $\mathbf{r}_{n}$ for $n=1,\dots,N$ denote the $N$ locations along the flight path where these multi-frequency signals are emitted and measured. For each location, we set $\mathbf{r}_{\text{src}}=\mathbf{r}_{n}$ and solve the direct scattering problem described above for frequency $\omega_{m}$. We take as measurements evaluations of $u_{0}(\mathbf{r}_{n};\omega_{m})$ corresponding to the scattered field $u_{0}$ evaluated at frequency $\omega_{m}$ at position $\mathbf{r}_{n}$. Through repeated solution of the direct scattering problem over all $M$ frequencies and all $N$ source/measurement positions, we obtain the measurement matrix $D\in\mathbb{C}^{M\times N}$. Measurements are typically corrupted by additive noise, so we model the entries of the data matrix as the sum,

with $\eta_{mn}$ denoting additive Gaussian white noise.

Note that we have considered only frequency-flat, isotropic point sources. Additionally, we use point-wise evaluations of the scattered field as measurements. Different sources and measurements can be used through an appropriate modification of the direct scattering problem. However, the measurement matrix we use here is fundamental in that other sources and measurements can be computed from it.

A key challenge in this inverse scattering problem is the ground bounce signals in the measurements. Ground bounce signals correspond to the portion of measurements that are reflections of the incident field by the air-soil interface. They contribute more to the overall power in measurements than signals scattered by the subsurface target. However, those signals carry little or no information about the target. Therefore, we must remove ground bounce signals from measurements to be able to image the subsurface target.

When the air-soil interface is known, one can model these reflections and use those to remove them from measurements. Otherwise, one needs to consider alternative methods. In what follows, we describe how to use principal component analysis (PCA) to approximately remove ground bounce signals from measurements.

For PCA, we first compute the singular value decomposition (SVD) of the data matrix so that $D=U\Sigma V^{H}$ where $[\cdot]^{H}$ denotes the Hermitian (conjugate) transpose. We assume the non-zero singular values $\sigma_{1},\sigma_{2},\dots,\sigma_{r}$ along the diagonal of $\Sigma$ are in descending order. The key idea with PCA is that ground bounce signals are stronger than the scattered signals, so they are associated with the largest singular values. Hence, we consider

instead of $D$. Here, $J$ is a user-defined truncation, $\mathbf{u}_{j}$ and $\mathbf{v}_{j}$ are the $j$th columns of $U$ and $V$, respectively. Provided that $J$ is chosen well, the resulting processed data matrix $\tilde{D}_{J}$ is amenable to imaging methods such as Kirchhoff migration which we explain below.

Using PCA to approximately remove ground bounce signals is completely data-driven. It does not require a priori knowledge of the electromagnetic properties of soil. Rather, it is based off of the assumption that the ground bounce signals dominate the measurements over the scattered signals by the target.

Kirchhoff migration (KM), also known as back-propagation or reverse time-migration, is a so-called sampling method [12]. Rather than having to solve an optimization problem to reconstruct an image of the target, KM requires only the evaluation of an imaging function over some region of interest, which we call the imaging region. Let $\mathbf{y}$ denote a point in this imaging region. The KM imaging function is given by

with $\tilde{d}_{mn}$ denoting the entry of $\tilde{D}_{J}$, $a_{mn}(\mathbf{y})$ denoting what we call the illuminations and $[\cdot]^{\ast}$ denoting the complex conjugate.

To define the illuminations, we assume we know the mean elevation $L$ of the flight path over the air-soil interface. Additionally, we assume we know the relative refractive index for soil denoted by $\epsilon_{r1}$. Let $\mathbf{r}_{n}=(x_{n},L)$, $\mathbf{y}=(\xi,\zeta)$, and $k_{m}=\omega_{m}/c$. The illuminations are given by

The illuminations given in (20) correspond to the phase accumulated from scattering by a point target at location $\mathbf{y}$ below a flat air-soil interface on $z=0$ using the Fresnel approximation. Additionally, illuminations (20) assume only one interaction between this point target and the flat air-soil interface, whereas the full boundary value problem includes infinitely-many of these interactions.

We form an image by plotting $|I^{\text{KM}}|$ over the imaging region. When KM is effective, this plot reveals spatial information about the target. Note that this implementation of KM starts with measurements defined by (17) which include additive measurement noise. We then compute (18), the processed data matrix $\tilde{D}$ resulting from PCA which approximately removes ground bounce signals. Finally, this KM implementation evaluates (19) with illuminations (20) defined for a flat air-soil interface.

In this formulation of KM, we assume a priori knowledge of a ground permittivity $\epsilon_{r1}$ and the mean height of the flight path over the air-soil interface $L$. Other than those two parameters, no other parameters are needed to evaluate (19).

We now show examples of images produced using the methods described above on data simulated using the MFS. For these examples, we have used $41$ frequencies uniformly sampling the bandwidth between $3.5$ GHz and $5.5$ GHz (at $50$ MHz steps) similarly to what was done by Garcia-Fernandez et al [4]. We have set the elevation above the mean of the air-soil interface to be $L=75$ cm and use $35$ spatial locations unformly sampling the synthetic aperture of size $a=102$ cm (at $3$ cm steps).

We consider the air-soil interface to be one realization of a Gaussian-correlated random rough surface with RMS height $0.4$ cm and correlation length $8$ cm. We showed recently that random rough surfaces with these parameters exhibit enhanced backscattering indicating significant multiple scattering caused by the surface roughness. As suggested by Daniels [13], we set $\epsilon_{r1}=9$ for soil and $\epsilon_{r2}=2.3$ for the target. For the MFS computations to generate simulated measurements, we used $P=512$ points to sample the surface over the interval of length $L=400$ cm.

In Fig. 3 we show results for the same “kite” target used in Fig. 2. We have added measurement noise to the simulated measurements so that SNR $\approx 25$ dB. The singular values of the data matrix $D$ are shown in Fig. 3(a). We observe from those results that the first two singular values are significantly larger than the remaining singular values with large gaps between them. In Figs. 3(b)–(d) we show KM images applied to the original data matrix $D$ (Fig. 3(b)), the PCA processed data matrix $\tilde{D}_{1}$ with one term removed (Fig. 3(c)), and the PCA processed data matrix $\tilde{D}_{2}$ with two terms removed (Fig. 3(d)). These KM images show plots of the absolute value of (19) normalized by its maximum value.

These results show how PCA effectively removes the ground bounce signals that interfere with the identification and location of the target. In Fig. 3(b) we find that the image is concentrated near $z=0$ corresponding to the air-soil interface. With the contribution by the first singular value removed, the image shown in Fig. 3(c) identifies the target but also has artifacts located near the air-soil interface. After removing the contributions by the first two singular values, we find in Fig. 3(d) that the resulting image identifies the target. Specifically, it appears to identify the single point on the boundary of the target that is closest to the synthetic aperture.

In Fig. 4 we consider the KM images for differently shaped targets: (a) a circle of radius $3.5$ cm centered at $(3,-14)$ cm, (b) an ellipse with major axis $3.5$ cm and minor axis $2.5$ cm centered at $(3,-14)$ cm, (c) a star whose boundary is given by $(3,-14)+(2.5+0.6\cos 5t,-3\sin 5t)$ cm for $0\leq t\leq 2\pi$, and (d) a rectangle with width $7$ cm and height $5$ cm centered at $(3,-15)$ cm. All of the targets have relative dielectric constant $\epsilon_{r2}=2.3$. All of these images have used PCA processed data $\tilde{D}_{2}$ and these KM images show plots of the absolute value of (19) normalized by its maximum value.

The maximum value of each of the KM images in Fig. 4(a)–(d) are plotted as a red “+” symbol. All of these results indicate that the KM image concentrates on a single point near the boundary of the target facing the synthetic aperture. There are some differences of the spread of the KM image about that point most notably for the rectangle. Nonetheless, these KM images do no provide too much spatial information other than this point. This behavior of KM images is due to the limitations in measurements. Specifically, SAR measurements (17) only subtend a relatively small portion of all scattering angles. Moreover, monostatic SAR measurements are restricted to just the backscattering or retroreflected direction. These space/angle limitations correspond to a severely limited aperture imaging problem. We cannot expect to recover very much spatial information about the target. The images shown in Fig. 4 suggest that one cannot obtain much more than a single point associated with a target. This point can be used to identify a target and locate it, but not much else.

The first simplifying assumption we make is to consider only first-order interactions between the target and the air-soil interface. Even when there is multiple scattering in the direct scattering problem, methods for the inverse scattering problem based only on first-order interactions are effective. It is understood that these first-order interaction approximations carry sufficient phase information regarding boundaries where discontinuities in wave fields occur.

Assuming only first-order interactions between the target and the air-soil interface, we can consider the following procedure to approximate the solution of the direct scattering problem.

1. 1.

   For a source located above the air-soil interface, compute reflection and transmission by the rough air-soil interface with no target below it.

2. 2.

   Use the field transmitted across the air-soil interface to determine the field incident on the target.

3. 3.

   Solve the scattering problem by the target using that incident field.

4. 4.

   Compute the field scattered by the target that is transmitted back across the air-soil interface from soil into air and evaluate that result at the receiver.

Using block system (15) for the MFS implementation we have used here, this procedure corresponds to first solving

then solving

and finally solving

With these results, we approximate measurements through evaluation of (4) using $a_{p}\approx a_{p}^{(0)}+a_{p}^{(1)}$.

In Fig. 5 we show the absolute difference ($\log_{10}$-scale) between the MFS solution of the full boundary value problem described in Section III and the first-order target-interface interaction approximation described here. For these results, we have used one realization of a Gaussian-correlated random rough surface with RMS height $0.4$ cm and correlation length $8$ cm. The relative dielectric constant for soil is $\epsilon_{r1}=9$ and for the target is $\epsilon_{r2}=2.3$. We used a point source at location $\mathbf{r}_{\text{src}}=(-25,75)$ cm at frequency $5.1$ GHz.

This result is characteristic of this first-order target-interface interaction approximation. The absolute difference is nearly three orders of magnitude smaller above the air-soil interface than below it. This non-uniformity is to be expected since the multiple interactions between the target and interface that are neglected in this approximation are most pronounced below the air-soil interface. However, from the point of view of modeling measurements taken above the air-soil interface, this result suggest that this approximation is accurate and useful simplifying the model of SAR measurements.

In Fig. 6 we show comparisons of the KM images produced using measurements computed (a) from the full direct scattering problem and (b) from the first-order target-interface interaction approximation. The air-soil interface is the same realization of a Gaussian-correlated random rough surface with RMS height $0.4$ cm and correlation length $8$ cm for both data sets. The relative dielectric constant in soil is $\epsilon_{r1}=9$ and in the the target is $\epsilon_{r2}=2.3$. Additive measurement noise has been added so that SNR $\approx 25$ dB. Both of these images are produced using $\tilde{D}_{2}$ from their respective data matrices. These results are nearly indistinguishable from one another indicating that the first-order approximation is valid for modeling measurements. The location where the maximum of the KM image is attained for the first-order approximation is exactly the same as for for the full direct scattering problem data.

The shower curtain effect is a phenomenon related to imaging through a scattering medium [14, 15]. To explain it consider a source and a detector separated by a fixed distance with some multiple scattering region of finite extent (i.e. the shower curtain) somewhere in between. The location of the shower curtain relative to the source affects the image quality. Blurring is stronger when the shower curtain is closer to the receiver than the source. When the shower curtain is closer to the source, blurring is reduced in comparison. Scattering by the shower curtain causes blurring through introducing angular diversity. When the receiver is close to the shower curtain, it measures a substantial fraction of large-angle scattered fields which, in turn, introduces a mixture of phases that cause blurring. In contrast, when the receiver is far from the shower curtain, it does not measure those large-angle scattered fields thereby reducing blur.

When we use the first-order target-interface approximation described above, scattering by the target is a secondary source below the rough air-soil interface which acts as the shower curtain. Since a buried landmine is typically situated at distances on the order of the wavelength below the air-soil interface, its distance to it is shorter than the distance from the UAV to it. Considering the shower curtain effect, we assume that the roughness of the air-soil interface can be neglected when computing scattering by the target. For that case, we solve (21) and (22) as before, but replace (23) with

where the entries of the blocks $A_{1}^{\text{flat}}$, $A_{2}^{\text{flat}}$, $A_{3}^{\text{flat}}$, and $A_{4}^{\text{flat}}$ are computed as (11) but for a flat air-soil interface in which $\mathbf{r}_{p}^{\text{int}}=(x_{p},0)$ for $p=1,\dots,P$. We call this approximation the first-order/flat interface approximation.

In Fig. 7 we show results using this first-order/flat interface approximation for the same problem as in Fig. 5. The absolute difference ($\log_{10}$-scale) between the MFS solution for the full boundary value problem and this first-order/flat interface approximation is shown in Fig. 7(a). In comparison with Fig. 5 we see that the error for this approximation in the region above the air-soil interface is not markedly different and suggests its validity for modeling measurements.

To test its validity for modeling measurements for imaging, we show in Fig. 7(b) the resulting KM image normalized by its maximum value using $\tilde{D}_{2}$ of the data simulated using this approximation. The rough air-soil interface is the same one used for the results shown in Fig. 6. Additionally, we have added measurement noise to the approximate measurements so that SNR $\approx 25$ dB just as with the results shown in Fig. 6. This image formed using the first-order/flat interface approximation to simulate measurements is nearly indistinguishable with those shown in Fig. 6. Moreover, the location where the KM image attains its maxiumum value is exactly the same as that in Fig. 6(a). In the context of faithfully reproducing KM imaging results, these results show that this first-order/flat interface approximation is sufficiently accurate.

All of the KM images we have shown exhibit a concentration about a single location where the KM image attains maximum value. The image for the rectangular target shown in Fig. 4(d) is the most different in that the concentration of the image about that point seems to follow the flat boundary. As we discussed in Section IV, this GP-SAR problem is a severely limited aperture imaging problem so we do not expect to recover much spatial information about the subsurface target. In light of these severe limitations, we consider simplifying a model of measurements by restricting our attention to a point-like target.

We know from elementary scattering theory that the field scattered by a particle in the far-field is proportional to the free-space Green’s function modified by a scattering amplitude [14]. It is the scattering amplitude that contains information about the particle such as its size, shape, and material properties. Limitations in GP-SAR measurements result in a severely limited sample of the scattering amplitude, but the KM images suggest the underyling Green’s function is robustly obtained.

Suppose we approximate the target below a rough air-soil interface by a point target located at $\mathbf{r}_{0}^{\text{tar}}$ below a flat air-soil interface. To consider a point target model within the MFS implementation we have used here, we solve (21) as before, but instead of solving (22) and (23), we instead solve

with

denoting the field exciting the point target. The components of the $P$-vectors $\mathbf{f}_{1}$ and $\mathbf{f}_{2}$ are given by

for $p=1,\dots,P$, respectively. We call the scalar parameter $\varrho$ the reflectivity which gives the scattering strength of the point target.

In Fig. 8 we show the KM image using the first-order/flat interface approximation with a point target located at $\mathbf{r}_{0}^{\text{tar}}=(0,-10.36)$ cm. This target location corresponds to the location where the KM image attains its maximum value using the full direct scattering problem measurements. The reflectivity has been arbitrarily set to $\varrho=8$. Additive noise has been added to these approximate measurements so that SNR $\approx 25$ dB.

The KM image shown in Fig. 8 is nearly indistinguishable with all other images shown here. This result strongly indicates that this first-order/flat interface approximation with a point target model captures the spatial information contained in measurements based off of the full direct scattering problem. Alternatively, these results suggest that any additional sophistication in modeling measurements is lost when considering KM images of targets.

Recall that entry $d_{mn}$ of the $M\times N$ data matrix $D$ gives the measurements at frequency $\omega_{m}$ at spatial location $\mathbf{r}_{n}$. We have discussed the following simplifying approximations above.

1. (i)

   First-order target-interface interaction

2. (ii)

   Flat interface for scattering by the target

3. (iii)

   Point target model

The results above show that within the context of KM imaging these approximations are valid in that they faithfully capture the qualitative features of KM images using data simulated from solutions of the full boundary value problem. These approximations lead to the following model for GP-SAR measurements of a subsurface target:

with the matrix $R^{\text{rough}}$ containing ground bounce signals by the rough air-soil interface, and the matrix $S^{\text{flat}}(\mathbf{r}_{0})$ containing signals scattered by a point target at location $\mathbf{r}_{0}$ with reflectivity $\varrho=1$ situated below a flat air-soil interface.

Within the context of this model, we can interpret the PCA and KM imaging method we have discussed above as follows. After computing the SVD of the measurements: $D=U\Sigma V^{H}$, we compute

with

denoting the error in removing the ground bounce signal using the contributions made by the first $J$ singular values. When $E$ is small, it acts effectively as additive noise to $\varrho S(\mathbf{r}_{0})$. The evaluation of $I^{\text{KM}}(\mathbf{y})$ given in (19) with illuminations (20) is designed to “unwind” the phase accumulated through propagation of the scattered field by a point-like target below a flat air-soil interface at the search location $\mathbf{y}$. Consequently, when the search point is in a neighborhood of the point target position, the phases of the data and illuminations match and lead to the absolute value $|I^{\text{KM}}|$ developing a peak.

By systematically testing and validating each of the simplifying assumptions above, we have established the validity of model (28). A key point to this model is the recognition of the limitations in GP-SAR measurements and understanding how those limitations affect the overall information content in those measurements. This model provides valuable insight into how the PCA and KM imaging method we have used here works.

Although our study uses only simulations of the direct scattering problem for scalar waves in two dimensions, the principles behind these basic assumptions extend naturally to more complicated three-dimensional scattering problems. For that case, the entries of matrices $R$ and $S(\mathbf{r}_{0})$ correspond to the respective field values for the three-dimensional problem. Moreover, measurement model (28) can be modified to incorporate more sophistication when that is necessary. For example, we can include additional $S$ matrices for multiple subsurface targets. If the point target model is not valid for a related, but different problem, one can replace the matrix $S$ with a more sophisticated scattering model.

We have recently proposed a dispersive point target model where the reflectivity $\varrho$ varies with frequency [16]. That target model is easily incorporated in (28) as

where $\varrho_{m}$ corresponds to the reflectivity at frequency $\omega_{m}$. From this dispersive point target model, we easily extend it further to include anisotropy in the reflectivity using

with $\odot$ denoting the Hadamard (element-wise) matrix product. The entries of the $M\times N$ matrix $F$ are the reflectivities $\varrho_{mn}$ which now vary in frequency and space. These two extensions open opportunities for quantitative imaging methods that seek to determine reflectivities rather than seeking to determine spatial information about targets. Given the inherent space/angle limitations in GP-SAR measurements, these quantitative imaging problems may be more practically useful than seeking to find more spatial information about targets.

To investigate the validity of (32), we consider $R^{\text{rough}}$ and $S^{\text{flat}}(\mathbf{r}_{0}^{\text{tar}})$ that we used to evaluate the point target assumption in Fig. 8. We then consider the measurements used for Fig. 7 that used the first-order target-interface approximation and scattering by the kite-shaped target below a flat air-soil interface. To determine the reflectivity matrix $F$, we subtract $R^{\text{rough}}$ from those measurements and then divide that result element-wise by the entries of $S^{\text{flat}}(\mathbf{r}_{0}^{\text{tar}})$. A plot of the absolute values of this result are shown in Fig. 9(a), which we have labelled as “model.” We do the same computation using the measurements from the full direct scattering problem and those results are shown in Fig. 9(b), which we have labelled as “data.” These results show that the measurements computed from full direct scattering problem can be modeled accurately using (32) which includes the reflectivity matrix $F$.

The results shown in Fig. 9 make explicit use of $R^{\text{rough}}$ which tacitly implies that the rough air-soil interface is known. We now assume it is not known and apply PCA as we have done above by truncating the contributions made by the first two singular values to obtain $\tilde{D}_{2}$. For that case, we estimate the reflectivity matrix $F$ by dividing $\tilde{D}_{2}$ element-wise by the entries of $S^{\text{flat}}(\mathbf{r}_{0}^{\text{tar}})$. A plot of the absolute values of this result are shown in the right plot of Fig. 10. This result shows that this estimate has large, spurious oscillations due to the error $E$ introduced in (30). Although $E$ only acts as additive measurement noise for imaging using KM, its effect on estimating $F$ here is much stronger. In light of these strong oscillatory errors, we consider the spatial average of absolute values:

Then we normalize this spatial average by its $2$-norm and denote it by $\widehat{|\varrho|}(f)$. A plot of $\widehat{|\varrho|}(f)$ for the results shown in Fig. 9 and this PCA result are shown in the right plot of Fig. 10. Even though these normalized frequency spectra do not quantitatively agree with one another, they share several characteristic features, especially within the sub-band between $4$ GHz and $5$ GHz.

In Fig. 4 we showed that the imaging method we use here does not produce images that show clear differences for markedly different target shapes. In Fig. 11 we show $\widehat{|\varrho|}(f)$ for each of these targets. These results show that these frequency spectra are markedly different for the different targets. We have recently studied how these normalized frequency spectra can be used to classify subsurface targets [17]. From these results we see that there is characterizing information contained in GP-SAR measurements. By considering (32), we are able to propose novel quantitative methods that enable target classification.