Quantitative signal subspace imaging

Multiple signal classification (MUSIC) was first introduced in signal processing [1] and has found multiple applications in imaging [2–5] in the last twenty years. MUSIC is one example of a sampling method. It forms an image through evaluation of a functional over an imaging region. MUSIC is often used for imaging point-like targets, but it has also been shown to be effective for imaging extended targets [6, 7]. Other examples of sampling methods are the linear sampling method [8, 9] and the factorization method [10]. Cheney [11] established an insightful equivalence between MUSIC and the factorization method for point targets. For a review on sampling and probe methods, we refer readers to [12] and references contained therein.

Sampling methods provide information about the support of the reflectivity, i.e. the location and shape of the scatterer. They typically do not provide quantitative information about the reflectivity of targets. However, the effective and efficient recovery of quantitative information about the reflectivity of targets is important, especially in remote sensing applications where on-the-fly reconstructions of the reflectivity are needed. Nonetheless, sampling methods are quite useful in quantitative imaging because they can be used as a priori information for an ℓ₂-minimization problem to determine the unknown quantities of interest [13–16]. Crocco et al [17] have used the linear sampling method as a preprocessing tool for quantitative imaging. Giorgi and Haddar [18] have also used a sampling method to recover the support of an extended target and then calculate the first eigenvalue of the far-field problem to obtain an approximation of the index of refraction.

In this paper we extend MUSIC for quantitative imaging using single frequency data collected on an array. We apply this quantitative imaging method to a single point target, multiple point targets, and an extended target. MUSIC uses only the noise subspace of the data to determine the support of small point scatterers. Our elementary extension of MUSIC makes use of the signal and noise subspaces together to construct an explicit imaging functional that recovers the support and magnitude of the reflectivity of point targets. For extended targets, this imaging functional recovers a scalar quantity that is related to the product of the volume and relative dielectric permittivity. By considering an alternate, yet still explicit, imaging functional, we are able to recover the full complex reflectivity of point targets after using the first imaging functional to locate them.

Moscoso et al [5] have shown that MUSIC exactly recovers the location of point targets when the data admit a special factorization. This exact reconstruction is true when targets lie on grid points in the mesh used to discretize the imaging region. When targets are off-grid, it is well known that the single frequency MUSIC image loses resolution dramatically, especially in the range direction (the coordinate orthogonal to the linear array corresponding to the main direction of propagation) [5, see figures 3 and 4]. Because our extension of MUSIC yields an explicit imaging functional, we are able to carry out a relatively simple resolution analysis. This resolution analysis provides valuable insight into this previously observed behavior. Specifically, we show that our imaging functional has cross-range (coordinate parallel to the linear array) resolution that is $O(\sqrt{{\epsilon}}\lambda L/a)$ and a range resolution that is $O(\sqrt{{\epsilon}}\lambda {L}^{2}/{a}^{2})$ with λ denoting wavelength, L denoting range, a denoting array aperture and related to the noise level. Note that these resolution estimates are simply the classical resolution estimates for array imaging multiplied by the factor $\sqrt{{\epsilon}}$. These resolution analysis results explain the subwavelength resolution obtained by MUSIC for large signal-to-noise ratios (SNRs) corresponding to 0 < ≪ 1. Moreover, they explain why the range resolution deteriorates faster since it is typical that L²/a² ≫ L/a. Finally, we determine that can be interpreted as a regularization parameter for the imaging functional which can be used to help address the off-grid target issue. It is in this way that we find that measurement noise surprisingly regularizes this imaging functional making the detection of targets more robust and effective. We show that this regularization-by-noise also plays a useful role in recovering the complex reflectivity of point targets.

The remainder of this paper is as follows. In section 2 we describe the single-frequency array imaging system to establish notation that we use throughout this paper. In section 3 we discuss imaging a single point target. Here we introduce our extension to MUSIC, give its resolution analysis, discuss its regularization-by-noise, and introduce an alternate imaging functional for recovering the complex reflectivity. In section 4 we apply the results from the previous section to the problem of multiple point targets. By artificially adjusting the noise level through a judicious choice of the user-defined parameter , we show that the quantitative imaging method is robust and effective for multiple point targets. We consider the quantitative information contained in measurements of scattering by an extended target in section 5. Although this quantitative information is not as directly useful as that for point targets, we show that it is robust even when the image of the support is degraded by noise. Section 6 contains our conclusions.

We consider a specific imaging system in which a multi-element, linear array with M transducers, each of which acts as both a source and a detector, operating at a single system frequency. The complete set of measurements yields an M × M matrix. The objective is to determine the support and reflectivity of targets in the imaging region from this measurement matrix. Moscoso et al [5] have shown that MUSIC for this measurement matrix can be cast as a linear system that admits a special factorization allowing for exact reconstructions in the noiseless case when point targets are collocated with grid points used to discretize the imaging region.

The measurements come from a suite of experiments, each of which consists of an individual array element emitting a signal and all array elements measuring the backscattered response. Figure 1 gives a schematic for a single experiment. Here we assume that the system bandwidth is narrow and therefore the signals may be approximated by single-frequency, time-harmonic waves. The complete suite of experiments includes all results from each array element acting as the emitter. The result is the measurement matrix, $B\in {\mathbb{C}}^{M\times M}$.

Zoom In Zoom Out Reset image size

Let x _i for i = 1, ..., M denote the positions of the elements of the array. We assume that the signal emitted by each element of the array is a spherical wave centered at that element's position,

$$\begin{equation}{G}_{0}(\boldsymbol{x}-{\boldsymbol{x}}_{i})=\frac{{\text{e}}^{\mathrm{i}k\vert \boldsymbol{x}-{\boldsymbol{x}}_{i}\vert }}{4\pi \vert \boldsymbol{x}-{\boldsymbol{x}}_{i}\vert },\quad i=1,\dots ,M.\end{equation} \tag{ 1 }$$

The wavenumber is k = 2π/λ with λ denoting the wavelength. For convenience, we introduce the vector of signals emitted by the elements of the array,

$$\begin{equation}\mathbf{g}(\boldsymbol{x})=\left[\begin{matrix}\hfill {G}_{0}(\boldsymbol{x}-{\boldsymbol{x}}_{1})\hfill \\ \hfill {G}_{0}(\boldsymbol{x}-{\boldsymbol{x}}_{2})\hfill \\ \hfill {\vdots}\hfill \\ \hfill {G}_{0}(\boldsymbol{x}-{\boldsymbol{x}}_{M})\hfill \end{matrix}\right].\end{equation} \tag{ 2 }$$

Suppose that a single point target with reflectivity ρ₀ is located at position y ₀ in the medium. The field scattered by this point target is a spherical wave centered at y ₀ whose amplitude is the product of the reflectivity and the field incident on the target. It follows that the field measured on element j of the array due to the signal emitted by element i of the array is given by [5, 6]

$$\begin{equation*}{b}_{ij}={\rho }_{0}{G}_{0}({\boldsymbol{x}}_{j}-{\boldsymbol{y}}_{0}){G}_{0}({\boldsymbol{y}}_{0}-{\boldsymbol{x}}_{i}).\end{equation*}$$

In terms of $\mathbf{g}({\boldsymbol{y}}_{0})={\Vert}\mathbf{g}({\boldsymbol{y}}_{0}){\Vert}\hat{\mathbf{g}}({\boldsymbol{y}}_{0})$ given in (2) where $\hat{\mathbf{g}}({\boldsymbol{y}}_{0})$ is a unit vector, the measurement matrix B is given by

$$\begin{equation}B={\rho }_{0}{\Vert}\mathbf{g}({\boldsymbol{y}}_{0}){{\Vert}}^{2}\hat{\mathbf{g}}({\boldsymbol{y}}_{0})\hat{\mathbf{g}}{({\boldsymbol{y}}_{0})}^{\text{T}},\end{equation} \tag{ 3 }$$

where the superscript T denotes the transpose (without complex conjugation).

From (3) we find that B is a rank-one matrix. The one-dimensional column space of B, given by $\text{span}\left\{\hat{\mathbf{g}}({\boldsymbol{y}}_{0})\right\}$, is called the signal subspace. The M − 1 dimensional subspace orthogonal to the signal subspace is called the noise subspace. We introduce the projections: ${P}_{\;\text{signal}}=\hat{\mathbf{g}}({\boldsymbol{y}}_{0})\hat{\mathbf{g}}{({\boldsymbol{y}}_{0})}^{\text{H}}$ and ${P}_{\;\text{noise}}=I-\hat{\mathbf{g}}({\boldsymbol{y}}_{0})\hat{\mathbf{g}}{({\boldsymbol{y}}_{0})}^{\text{H}}$ onto the signal and noise subspaces, respectively. Here, the superscript H denotes the Hermitian (complex conjugate transpose), and I denotes the M × M identity matrix.

In MUSIC one considers the vector corresponding to illuminating the medium at search point  x by all of the elements of the array. This vector of illuminations is g( x ) given in (2). In terms of g( x ), the MUSIC imaging functional is

$$\begin{equation}{F}_{\;\text{MUSIC}}(\boldsymbol{x})={\Vert}{P}_{\;\text{noise}}\mathbf{g}(\boldsymbol{x}){\Vert}.\end{equation} \tag{ 4 }$$

Here, ||⋅|| denotes the Euclidean norm of the vector. This non-negative functional is zero identically when x = y ₀ because g( y ₀) lies entirely in the signal subspace, and hence is orthogonal to the noise subspace. Otherwise, the evaluation of this imaging functional is non-zero because g( x ) has a component in the noise subspace. Therefore, by plotting 1/F_MUSIC( x ), we form an image that has a peak at y ₀ thereby identifying the location of the point target.

MUSIC uses the noise subspace to form an image of the target from which only the support of the target can be recovered. By also including the signal subspace, we can recover the reflectivity of the target in the following way. We first compute the singular value decomposition (SVD) of the measurement matrix, B = UΣV^H. According to (3) the SVD of the measurement matrix without noise due to a point target will yield only one non-zero singular value, σ₁ = |ρ₀|||g( y ₀)||². The corresponding singular vectors are ${\mathbf{u}}_{1}=\hat{\mathbf{g}}({\boldsymbol{y}}_{0}){\text{e}}^{\mathrm{i}({\theta }_{0}/2+\theta )}$ and ${\mathbf{v}}_{1}={\hat{\mathbf{g}}}^{\ast }({\boldsymbol{y}}_{0}){\text{e}}^{-\mathrm{i}({\theta }_{0}/2-\theta )}$ with ${\rho }_{0}=\vert {\rho }_{0}\vert {\text{e}}^{\mathrm{i}{\theta }_{0}}$ and θ denoting an arbitrary angle of rotation on the complex plane. To consider imaging with measurement noise, we introduce a dimensionless parameter in the M × M matrix,

$$\begin{equation}{{\Sigma}}_{{\epsilon}}^{+}={\sigma }_{1}^{-1}\mathrm{diag}(1,{{\epsilon}}^{-1},\dots ,{{\epsilon}}^{-1}),\end{equation} \tag{ 5 }$$

which is a pseudo-inverse for Σ. Instead of (4), we consider

$$\begin{equation}{F}_{{\epsilon}}(\boldsymbol{x})=\mathbf{g}{(\boldsymbol{x})}^{\text{H}}U{{\Sigma}}_{{\epsilon}}^{+}{U}^{\text{H}}\mathbf{g}(\boldsymbol{x}).\end{equation} \tag{ 6 }$$

We propose this functional for quantitative imaging of a point target to determine (i) the location y ₀ and (ii) the magnitude of reflectivity |ρ₀| since

$$\begin{align}{F}_{{\epsilon}}(\boldsymbol{x})& =\frac{1}{{\sigma }_{1}{\epsilon}}{\Vert}{P}_{\;\text{noise}}\mathbf{g}(\boldsymbol{x}){{\Vert}}^{2}+\frac{1}{{\sigma }_{1}}{\Vert}{P}_{\;\text{signal}}\mathbf{g}(\boldsymbol{x}){{\Vert}}^{2}\\ & =\frac{1}{\vert {\rho }_{0}\vert {\Vert}\mathbf{g}({\boldsymbol{y}}_{0}){{\Vert}}^{2}{\epsilon}}{\Vert}{P}_{\;\text{noise}}\mathbf{g}(\boldsymbol{x}){{\Vert}}^{2}+\frac{1}{\vert {\rho }_{0}\vert {\Vert}\mathbf{g}({\boldsymbol{y}}_{0}){{\Vert}}^{2}}{\Vert}{P}_{\;\text{signal}}\mathbf{g}(\boldsymbol{x}){{\Vert}}^{2}.\end{align} \tag{ 7 }$$

When x ≠ y ₀, the illumination vector g( x ) lies entirely in the noise subspace and so the second term in (7) evaluates to zero identically. Thus, F ( x ) behaves similarly to (4) when x ≠ y ₀. However, when x = y ₀ the first term of (7) evaluates to zero identically and the second term evaluates to 1/|ρ₀|. Therefore, when plotting 1/F ( x ), we form an image with peak at the location y ₀ whose value is the absolute value of the reflectivity |ρ₀| of the target.

The value of is chosen to be between the smallest singular value corresponding to the signal subspace and the largest singular value corresponding to the noise subspace. When the SNR is sufficiently large, choosing such a value of is relatively straight forward because the singular values corresponding to the signal subspace tend to be much larger than those corresponding to the noise subspace. Consequently, is set by the noise level in the measurements.

In figure 2 we show an image formed by plotting the 1/F ( x ) where F is given in (6) with = 10⁻⁶ on the xz-plane for a point target located at y ₀ = (kx₀, kL) = (0, 100) with magnitude of reflectivity |ρ₀| = 0.8. The measurement array with aperture ka = 80 has 41 elements equi-spaced on the interval −40 ⩽ kx ⩽ 40 on z = 0. We have added noise to measurements so that SNR = 4.5334  dB. This result shows that the image is peaked at the location of the target at which the magnitude is equal to |ρ₀|. Additionally, note that the region plotted about the point target, −0.02 ⩽ k(x − x₀) ⩽ 0.02 and −0.2 ⩽ k(z − L) ⩽ 0.2, is rather small compared to the wavelength indicating a high resolution achieved with this imaging method.

Zoom In Zoom Out Reset image size

3.1. Resolution analysis

We now determine the spatial resolution of the reconstructed image formed by evaluating the 1/F ( x ) where F ( x ) is given by (6). To keep the notation in this analysis simple, we consider only the xz-plane in what follows. We use a coordinate system in which the origin is situated at the center of the array that is aligned with the x-axis, the cross-range coordinate. The z-axis is the range coordinate. Extending this result to three dimensions is straightforward.

Theorem 3.1 (Resolution analysis). Let (x₀, L) with kL ≫ 1 denote the location of a point target with reflectivity ρ₀. Assuming that 0 < ≪ 1, the image formed by evaluating 1/F ( x ) where F ( x ) is given by (6) using an array with aperture size a (i) recovers |ρ₀|, (ii) has a cross-range resolution that is $O(\sqrt{{\epsilon}}\lambda L/a)$, and (iii) has a range resolution that is $O(\sqrt{{\epsilon}}\lambda {L}^{2}/{a}^{2})$.

Proof. Let x_i = −a/2 + (i − 1)a/(M − 1) for i = 1, ..., M denote the locations of the equi-spaced array elements over the aperture a. We write x _i = (x_i , 0), x _j = (x_j , 0), and y ₀ = (x₀, L). Since kL ≫ 1, we use the Fraunhofer approximation, $k\vert {\boldsymbol{x}}_{i}-{\boldsymbol{y}}_{0}\vert \sim kL+k{x}_{i}^{2}/(2L)-k{x}_{i}{x}_{0}/L$ as kL → ∞, and find

$$\begin{equation*}G({\boldsymbol{x}}_{i}-{\boldsymbol{y}}_{0})\sim \frac{{\text{e}}^{\mathrm{i}kL}}{4\pi L}{\text{e}}^{\mathrm{i}k{x}_{i}^{2}/(2L)-\mathrm{i}k{x}_{i}{x}_{0}/L}.\end{equation*}$$

Let ${\rho }_{0}=\vert {\rho }_{0}\vert {\text{e}}^{\mathrm{i}{\theta }_{0}}$ and

$$\begin{equation}{Q}_{0}={\text{e}}^{\mathrm{i}kL}\mathrm{diag}({\text{e}}^{\mathrm{i}k{x}_{1}^{2}/(2L)-\mathrm{i}k{x}_{0}{x}_{1}/L},\dots ,{\text{e}}^{\mathrm{i}k{x}_{M}^{2}/(2L)-\mathrm{i}k{x}_{0}{x}_{M}/L}).\end{equation} \tag{ 8 }$$

In terms of these quantities, the measurement matrix given in (3) is

$$\begin{equation}B=\left({\text{e}}^{\mathrm{i}{\theta }_{0}/2}{Q}_{0}{\mathbf{w}}_{0}\right)\frac{M\vert {\rho }_{0}\vert }{{(4\pi L)}^{2}}{\left({\text{e}}^{-\mathrm{i}{\theta }_{0}/2}{Q}_{0}^{\ast }{\mathbf{w}}_{0}\right)}^{\text{H}},\end{equation} \tag{ 9 }$$

with w₀ = M^−1/2 1_M where 1_M denotes the M-vector of ones.

Using (5) and ${\mathbf{u}}_{1}={\text{e}}^{\mathrm{i}{\theta }_{0}/2}{Q}_{0}{\mathbf{w}}_{0}$, we rewrite (6) as

$$\begin{align*}{F}_{{\epsilon}}(\boldsymbol{x})& =\frac{1}{{\epsilon}{\sigma }_{1}}\mathbf{g}{(\boldsymbol{x})}^{\text{H}}(I-{\mathbf{u}}_{1}{\mathbf{u}}_{1}^{\text{H}})\mathbf{g}(\boldsymbol{x})+\frac{1}{{\sigma }_{1}}\mathbf{g}{(\boldsymbol{x})}^{\text{H}}{\mathbf{u}}_{1}{\mathbf{u}}_{1}^{\text{H}}\mathbf{g}(\boldsymbol{x})\\ & =\frac{1}{{\sigma }_{1}}\left[\frac{1}{{\epsilon}}{\Vert}\mathbf{g}(\boldsymbol{x}){{\Vert}}^{2}-\left(\frac{1}{{\epsilon}}-1\right){\left\vert {\mathbf{w}}_{0}^{\text{T}}{Q}_{0}^{\text{H}}\mathbf{g}(\boldsymbol{x})\right\vert }^{2}\right].\end{align*}$$

According to the Fraunhofer approximation, the entries of g( x ) are

$$\begin{equation*}{g}_{j}(x,z)=\frac{{\text{e}}^{\mathrm{i}kz}}{4\pi z}\enspace {\text{e}}^{\mathrm{i}k{x}_{j}^{2}/(2z)}\enspace {\text{e}}^{-\mathrm{i}k{x}_{j}x/z},\quad j=1,\dots ,M,\end{equation*}$$

so ||g( x )||² = M/(4πz)², and

$$\begin{equation*}{\mathbf{w}}_{0}^{\text{T}}{Q}_{0}^{\text{H}}\mathbf{g}(x,z)=\frac{{\text{e}}^{\mathrm{i}k(z-L)}}{\sqrt{M}(4\pi z)}\sum\limits _{j=1}^{M}{\text{e}}^{\mathrm{i}k{x}_{j}^{2}(1/(2z)-1/(2L))}{\text{e}}^{-\mathrm{i}k{x}_{j}(x/z-{x}_{0}/L)}.\end{equation*}$$

Let

$$\begin{equation}A(x,z)=\sum\limits _{j=1}^{M}{\text{e}}^{\mathrm{i}k{x}_{j}^{2}(1/(2z)-1/(2L))}\enspace {\text{e}}^{-\mathrm{i}k{x}_{j}(x/z-{x}_{0}/L)}.\end{equation} \tag{ 10 }$$

It follows that

$$\begin{equation}{F}_{{\epsilon}}(x,z)=\frac{M}{{\epsilon}{\sigma }_{1}{(4\pi z)}^{2}}\left[1-\left(1-{\epsilon}\right)\frac{\vert A(x,z){\vert }^{2}}{{M}^{2}}\right].\end{equation} \tag{ 11 }$$

The second-degree Taylor polynomial approximation for A(x, z) is given by

$$\begin{equation}A(x,z)\approx M-\frac{\mathrm{i}k}{2{L}^{2}}{\mu }_{2}(z-L)-\frac{{k}^{2}}{2{L}^{2}}{\mu }_{2}{(x-{x}_{0})}^{2}-\left(\frac{{k}^{2}}{8{L}^{4}}{\mu }_{4}-\frac{\mathrm{i}k}{2{L}^{3}}{\mu }_{2}\right){(z-L)}^{2},\end{equation} \tag{ 12 }$$

with

$$\begin{equation}{\mu }_{2}=\sum\limits _{i=1}^{M}{x}_{i}^{2}={a}^{2}\frac{M(M+1)}{12(M-1)},\end{equation} \tag{ 13 }$$

and

$$\begin{equation}{\mu }_{4}=\sum\limits _{i=1}^{M}{x}_{i}^{4}={a}^{4}\frac{M(M+1)(3{M}^{2}-7)}{240(M-1)}.\end{equation} \tag{ 14 }$$

Using (12), we find that

$$\begin{align*}\frac{\vert A(x,z){\vert }^{2}}{{M}^{2}}& =1-\frac{{k}^{2}{\mu }_{2}}{M}\frac{{k}^{2}{(x-{x}_{0})}^{2}}{{(kL)}^{2}}\left(1-\frac{{k}^{2}{\mu }_{2}}{4M}\frac{{k}^{2}{(x-{x}_{0})}^{2}}{{(kL)}^{2}}\right)\\ & \quad +\left(\frac{{k}^{4}{\mu }_{2}^{2}}{4{M}^{2}}-\frac{{k}^{4}{\mu }_{4}}{4M}\right)\frac{{k}^{2}{(z-L)}^{2}}{{(kL)}^{4}}+O({(kL)}^{-5}).\end{align*}$$

Substituting this approximation into (11), we obtain

$$\begin{align}{F}_{{\epsilon}}(x,z)& =\frac{M}{{\epsilon}{\sigma }_{1}{(4\pi z)}^{2}}\left[{\epsilon}-(1-{\epsilon})\frac{{k}^{2}{\mu }_{2}}{M}\frac{{k}^{2}{(x-{x}_{0})}^{2}}{{(kL)}^{2}}\left(1-\frac{{k}^{2}{\mu }_{2}}{4M}\frac{{k}^{2}{(x-{x}_{0})}^{2}}{{(kL)}^{2}}\right)\right.\\ & \quad \left.+(1-{\epsilon})\left(\frac{{k}^{4}{\mu }_{2}^{2}}{4{M}^{2}}-\frac{{k}^{4}{\mu }_{4}}{4M}\right)\frac{{k}^{2}{(z-L)}^{2}}{{(kL)}^{4}}+O({(kL)}^{-5})\right].\end{align} \tag{ 15 }$$

Evaluating the reciprocal of approximation (15) on (x₀, L) and substituting σ₁ = M|ρ₀|/(4πL)², we find

$$\begin{equation*}\frac{1}{{F}_{{\epsilon}}({x}_{0},L)}=\vert {\rho }_{0}\vert .\end{equation*}$$

Thus, we recover the magnitude of the reflectivity at the target location.

To determine cross-range resolution, we evaluate approximation (15) on (x₀ + Δx, L) and find

$$\begin{align*}{F}_{{\epsilon}}({x}_{0}+{\Delta}x,L)& =\frac{1}{\vert {\rho }_{0}\vert }\left[1-\frac{(1-{\epsilon})}{{\epsilon}}\frac{{k}^{2}{\mu }_{2}}{M}\frac{{k}^{2}{\Delta}{x}^{2}}{{(kL)}^{2}}\left(1-\frac{{k}^{2}{\mu }_{2}}{4M}\frac{{k}^{2}{\Delta}{x}^{2}}{{(kL)}^{2}}\right)+O({(kL)}^{-5})\right].\end{align*}$$

We seek the value Δx* such that 1/F (x₀ + Δx*, L) = 1/(2F (x₀, L)) corresponding to the full-width/half-max (FWHM) and find that

$$\begin{equation*}1+\frac{(1-{\epsilon})}{{\epsilon}}\frac{{k}^{2}{\mu }_{2}}{M}\frac{{k}^{2}{\Delta}{x}^{\ast 2}}{{(kL)}^{2}}\left(1-\frac{{k}^{2}{\mu }_{2}}{4M}\frac{{k}^{2}{\Delta}{x}^{\ast 2}}{{(kL)}^{2}}\right)=0.\end{equation*}$$

In the limit as → 0⁺, we find upon substituting (13) that

$$\begin{equation}{\Delta}{x}^{\ast }\sim \pm 2\sqrt{3}\sqrt{\frac{M-1}{M+1}}\frac{\sqrt{{\epsilon}}L}{ka}=O\left(\frac{\sqrt{{\epsilon}}\lambda L}{a}\right).\end{equation} \tag{ 16 }$$

To determine range resolution, we evaluate approximation (15) on (x₀, L + Δz) and obtain

$$\begin{align*}{F}_{{\epsilon}}({x}_{0},L+{\Delta}z)& =\frac{{L}^{2}}{\vert {\rho }_{0}\vert {(L+{\Delta}z)}^{2}}\left[1+\frac{(1-{\epsilon})}{{\epsilon}}\left(\frac{{k}^{4}{\mu }_{2}^{2}}{4{M}^{2}}-\frac{{k}^{4}{\mu }_{4}}{4M}\right)\frac{{k}^{2}{\Delta}{z}^{2}}{{(kL)}^{4}}+O({(kL)}^{-5})\right].\end{align*}$$

We seek the FWHM value Δz* such that 1/F (x₀, L + Δz*) = 1/(2F (x₀, L)) and find that

$$\begin{equation*}\left[\frac{(1-{\epsilon})}{{\epsilon}}\left(\frac{{k}^{4}{\mu }_{2}^{2}}{4{M}^{2}}-\frac{{k}^{4}{\mu }_{4}}{4M}\right)\frac{1}{{(kL)}^{2}}-2\right]\frac{{\Delta}{z}^{\ast 2}}{{L}^{2}}-4\frac{{\Delta}{z}^{\ast }}{L}-1=0.\end{equation*}$$

In the limit as → 0⁺, we find upon substituting (13) and (14) that

$$\begin{equation}{\Delta}{z}^{\ast }\sim \pm \frac{24\sqrt{5}(M-1)}{\sqrt{5-74M+5{M}^{2}+120{M}^{3}-36{M}^{5}}}\frac{\sqrt{{\epsilon}}{L}^{2}}{k{a}^{2}}=O\left(\frac{\sqrt{{\epsilon}}\lambda {L}^{2}}{{a}^{2}}\right).\end{equation} \tag{ 17 }$$

This completes the proof.

3.2. Verifying resolution estimates

The resolution analysis of theorem 3.1 includes the familiar factors appearing in array imaging: λL/a for cross-range and λL²/a² for range. However, both factors are scaled by $\sqrt{{\epsilon}}$. To verify the resolution results in theorem 3.1, we numerically evaluate (6) where U and Σ are computed from the SVD of B given in (3), and determine the FWHM in range and cross-range. Note that B is computed using the exact Green's function (1). These numerical results are shown in figure 3.

Zoom In Zoom Out Reset image size

The left plot of figure 3 shows resolution results as a function of noise level for a point target with magnitude of reflectivity |ρ₀| = 0.8 located at kx₀ = 0 and kL = 100. The measurement array with aperture ka = 80 includes 41 elements sampling −40 ⩽ kx ⩽ 40 on kz = 0. The circle symbols correspond to the numerically computed values for the FWHM in range and the square symbols correspond to the numerically computed values for the FWHM in cross-range. The corresponding curves are least-squares linear fits to these data. For the range, the linear fit is log(Δz*) ≈ 0.5000 log() + 3.0910, and for the cross-range, the linear fit is log(Δx*) ≈ 0.5000 log() + 0.7885. These results are consistent with the $O(\sqrt{{\epsilon}})$ behavior of the range and cross-range resolutions.

The center plot of figure 3 shows resolution results as a function of the target range, kL for a point target with magnitude of reflectivity |ρ₀| = 0.8 located at kx₀ = 0. The measurement array with aperture ka = 80 includes 41 elements sampling −40 ⩽ kx ⩽ 40 on kz = 0. For the range, the linear fit is log(Δz*) ≈ 1.8878 log(kL) − 13.5770, and for the cross-range, the linear fit is log(Δx*) ≈ 0.9130 log(kL) − 11.4265. These results are consistent with the O((kL)²) behavior of the range and the O(kL) behavior of the cross-range resolutions.

The right plot of figure 3 shows resolution results as a function of the array aperture, ka for a point target with magnitude of reflectivity |ρ₀| = 0.8 located at kx₀ = 0 and kL = 100. The measurement array includes 41 elements sampling −ka/2 ⩽ kx ⩽ ka/2 on kz = 0. For the range, the linear fit is log(Δz*) ≈ −1.8412 log(ka) + 3.1718, and for the cross-range, the linear fit is log(Δx*) ≈ −0.8879 log(ka) − 3.3585. These results are consistent with the O((ka)⁻²) behavior of the range and the O((ka)⁻¹) behavior of the cross-range resolutions.

3.3. Off-grid target

A consequence of theorem 3.1 is given by the following corollary.

Corollary 3.3. In the limit as → 0⁺ corresponding to no noise and infinite precision arithmetic, 1/F ( x ) with F ( x ) given in (6) is zero everywhere except at the target location, (x₀, L) where it takes on the value |ρ₀|.

From the resolution point of view, corollary 3.3 indicates that when is small, this imaging method can yield a very high resolution quantitative image. However, the practical implication of corollary 3.3 is that recovering the location and magnitude of the reflectivity can be problematic. Suppose the target is 'off-grid' meaning that its location is not collocated with a grid point of the mesh used to reconstruct an image of the imaging region. Corollary 3.3 indicates that the target may not appear in the image. Moreover, if the location is not recovered exactly, the recovered magnitude of the reflectivity will have a large error that is an underestimate. This large quantitative error occurs because the functional producing the image falls off quickly from its peak value of |ρ₀| at the target location. Even if is set by finite precision arithmetic, the image may not reveal the target, or if it does, its quantitative recovery for |ρ₀| has a large error.

Measurement noise raises the value of . For that case, the $O(\sqrt{{\epsilon}})$ behavior of the resolution effectively regularizes (6). Thus, in the context of this imaging method, noise stabilizes this quantitative image reconstruction method. In fact, this regularization-by-noise opens the opportunity for a more robust method for recovering the location and magnitude of the reflectivity.

In figure 4 we show the evaluation of 1/F ( x ) on (x, L) to show its cross-range behavior (left) and on (x₀, z) to show its range behavior (right) for = 10⁻⁸ (blue curves) and = 10⁻⁵ (red curves). This figure shows how F ( x ) is regularized as increases thereby yielding a smoother result that peaks about the target location. In light of this regularization by , a method that seeks to determine the target location numerically will perform better with small-to-moderate noise than without noise since the corresponding imaging functional is smoother.

Zoom In Zoom Out Reset image size

To demonstrate how these ideas may work in practice, we consider the imaging region −0.5 ⩽ k(x − x₀) ⩽ 0.5 and −1 ⩽ k(z − L) ⩽ 1 centered about the target location (x₀, L). We plot 1/F ( x ) in this imaging region using a mesh of 100 × 100 equi-spaced grid points. Note that the exact target location is not a grid point on this mesh—it is displaced by half the grid size in both directions. Figure 5 shows results with = 10⁻¹⁶ (left), = 10⁻⁶ (center), and = 10⁻⁴ (right) using the same color scale for all three plots. One typically normalizes the images which does reveal a peak in each of these images. However, those peak values are quite different from one another. For = 10⁻¹⁶, 10⁻⁶ and 10⁻⁴ the peak values attained are 6.5381 × 10⁻¹¹, 0.3598, and 0.7903, respectively. By plotting the images without normalizing in figure 5, we show how affects the reconstructed image. Moreover, it gives a quantitative method for finding the exact location of the target since 1/F ( x ) attains its maximum there. This quantitative feature of this imaging method opens up the opportunity for using adaptive grid or numerical optimization methods for determining the target location, for example.

Zoom In Zoom Out Reset image size

3.4. Recovering the complex reflectivity

The imaging method described above only recovers the magnitude of the reflectivity, |ρ₀|, for the point target. However, the reflectivity for a point target is complex, in general. To see this, consider its scattering cross-section σ_s and total cross-section σ_t. These are related to one another through σ_t = σ_s + σ_a with σ_a denoting the absorption cross-section. Since scattering by a point target is isotropic, its scattering cross-section is σ_s = 4π|ρ₀|². By the optical theorem [19], its total cross-section is σ_t = 4π Im[ρ₀]/k. It follows that

$$\begin{equation*}{\rho }_{0}={\left[\frac{{\sigma }_{\text{s}}}{4\pi }-{\left(\frac{k{\sigma }_{\text{t}}}{4\pi }\right)}^{2}\right]}^{1/2}+\mathrm{i}\frac{k{\sigma }_{\text{t}}}{4\pi }.\end{equation*}$$

We now describe how to recover this complex reflectivity.

Using ${\rho }_{0}=\vert {\rho }_{0}\vert {\text{e}}^{\mathrm{i}{\theta }_{0}}$, we rewrite (3) as $B={\mathbf{u}}_{1}{\sigma }_{1}{\mathbf{v}}_{1}^{\text{H}}$ with σ₁ = |ρ₀|||g( y ₀)||²,

$$\begin{equation*}{\mathbf{u}}_{1}={\text{e}}^{\mathrm{i}({\theta }_{0}/2+\theta )}\enspace \hat{\mathbf{g}}({\boldsymbol{y}}_{0}),\end{equation*}$$

and

$$\begin{equation*}{\mathbf{v}}_{1}={\text{e}}^{-\mathrm{i}({\theta }_{0}/2-\theta )}\enspace {\hat{\mathbf{g}}}^{\ast }({\boldsymbol{y}}_{0}),\end{equation*}$$

with θ and the association of θ₀ to u₁ and v₁ arbitrary. Note that when we evaluate the functional

$$\begin{equation*}R(\boldsymbol{x})=\mathbf{g}{(\boldsymbol{x})}^{\text{T}}{\mathbf{v}}_{1}{\sigma }_{1}^{-1}{\mathbf{u}}_{1}^{\text{H}}\mathbf{g}(\boldsymbol{x}),\end{equation*}$$

on x = y ₀, we obtain

$$\begin{equation*}R({\boldsymbol{y}}_{0})=\frac{1}{\vert {\rho }_{0}\vert {\text{e}}^{\mathrm{i}{\theta }_{0}}}.\end{equation*}$$

Therefore, evaluating the reciprocal of R( y ₀) yields the complex reflectivity.

In light of this result, we consider the alternate imaging functional,

$$\begin{equation}{R}_{{\epsilon}}(\boldsymbol{x})=\mathbf{g}{(\boldsymbol{x})}^{\text{T}}V{{\Sigma}}_{{\epsilon}}^{+}{U}^{\text{H}}\mathbf{g}(\boldsymbol{x}),\end{equation} \tag{ 18 }$$

with ${{\Sigma}}_{{\epsilon}}^{+}$ given in (5). This imaging functional gives the complex reflectivity when evaluated at the position of the point target. However, this functional oscillates in a neighborhood of the target location. Therefore, it is not as useful as (6) for determining the location of a point target.

To recover the location and complex reflectivity of a point target, we propose the following two-stage procedure.

• (1)  
  Evaluate F given in (6) to determine the location and magnitude of the reflectivity of the point target.
• (2)  
  Evaluate R given in (18) on the determined location from (i) to determine the complex reflectivity.

We find that regularizes R just as it does for F ( x ) given in (6). Some example results are shown in figure 6. These results show that the 1/R becomes much smoother when increases. Again, the practical implications of this regularization is that a method to determine the complex reflectivity numerically will perform better with small-to-moderate noise than without noise since the corresponding imaging functional will be smoother.

Zoom In Zoom Out Reset image size

Consider N point targets located at y _n with complex reflectivities ρ_n for n = 1, ..., N. For this problem, the measurement matrix is

$$\begin{equation*}B=\sum\limits _{n=1}^{N}{\rho }_{n}\mathbf{g}({\boldsymbol{y}}_{n})\mathbf{g}{({\boldsymbol{y}}_{n})}^{\text{T}}.\end{equation*}$$

Let $\hat{G}$ denote the M × N matrix whose columns are given by $\hat{\mathbf{g}}({\boldsymbol{y}}_{n})$ for n = 1, ..., N. In terms of $\hat{G}$, the measurement matrix is

$$\begin{equation*}B=\hat{G}D{\hat{G}}^{\text{T}},\end{equation*}$$

with

$$\begin{equation}D=\mathrm{diag}\left({\rho }_{1}{\Vert}\mathbf{g}({\boldsymbol{y}}_{1}){{\Vert}}^{2},\dots ,{\rho }_{N}{\Vert}\mathbf{g}({\boldsymbol{y}}_{N}){{\Vert}}^{2}\right).\end{equation} \tag{ 19 }$$

The matrix $\hat{G}$ is not orthogonal in general. When we consider its QR-factorization, $\hat{G}=\hat{Q}\hat{R}$, we find that

$$\begin{equation}B=\hat{Q}\hat{R}D{\hat{R}}^{\text{T}}{\hat{Q}}^{\text{T}}.\end{equation} \tag{ 20 }$$

The singular values are related to the eigenvalues of $\hat{R}D{\hat{R}}^{\text{T}}$. In contrast to a single point target, we lose direct interpretation of the singular values for multiple point targets since $\hat{R}$ is unknown. Consequently, imaging through evaluation of F given in (6) allows for the recovery of point target locations, but does not provide direct quantitative information about the magnitude of their reflectivities.

We have found for imaging a single point target through evaluation of 1/F with F given in (6) that noise stabilizes the recovery of the location. This regularization-by-noise still holds for multiple point targets. However, this regularization yields an overall loss of resolution which may adversely affect resolving individual point targets. To see this we show in figure 7 the image produced through evaluation of 1/F for noisy data where SNR = 5.5146  dB. There are three point targets. The first point target is located at (kx₁, kL₁) = (−5.2, 96.0) with complex reflectivity ρ₁ = 2.9096i corresponding to σ_t1 = σ_s1 = 10. The second point target is located at (kx₂, kL₂) = (2.0, 102.0) with complex reflectivity ρ₂ = 3.5508i corresponding to σ_t2 = σ_s2 = 12. The third point target is located at (kx₃, kL₃) = (5.6, 98.0) with complex reflectivity ρ₃ = 2.2655i corresponding to σ_t3 = σ_s3 = 8. Although we find that the resulting image shown in figure 7 yields three peaks at the target locations, the regularization-by-noise leads to a blurred image of those targets, especially the second target whose range is largest among the three targets.

Zoom In Zoom Out Reset image size

Upon computing the SVD of the measurement matrix for the example shown in figure 7, we find that σ₁ = 5.7633 × 10⁻⁴, σ₂ = 1.7524 × 10⁻⁴, σ₃ = 7.9860 × 10⁻⁶, and σ_n < 2 × 10⁻⁹ for n ⩾ 4. As is typical, the singular values corresponding to the signal subspace are larger than those corresponding to the noise subspace. Therefore, we can often determine the dimension of the signal and noise subspaces by examining the singular values. Suppose the dimension of the signal subspace is r. When one can separate the singular values in this way, we propose imaging using F given in (6) with

$$\begin{equation}{{\Sigma}}_{{\epsilon}}^{+}={\sigma }_{1}^{-1}\mathrm{diag}(1,{\sigma }_{1}{\sigma }_{2}^{-1},\dots ,{\sigma }_{1}{\sigma }_{\text{r}}^{-1},{{\epsilon}}^{-1},\dots ,{{\epsilon}}^{-1}).\end{equation} \tag{ 21 }$$

With this imaging functional, one can adjust the value of that suitably resolves and regularizes the resulting image. By adjusting the value of , one can effectively change the resolution of the resulting image. But more importantly, adjusting the value of affects the smoothness of the peaks about the target locations allowing for numerically stable recovery of point target positions. Suppose one seeks to reconstruct an image over an imaging region using a fixed grid. One can adjust the value of so that the targets appear even if they are not collocated with a grid point (regularization) without too much blurring (resolution).

In figure 8, we plot 1/F with = 10⁻⁸ for the same noisy data used in figure 7. Next to each of the point targets, we have superimposed a plot of 1/F centered about the target location using a very fine mesh. Notice that in comparison to figure 7, the peaks in this image are more tightly focused about the target locations and blurring is significantly reduced. In the close-up images plotted on fine grids, the local maximum value is attained exactly at the target locations. These results indicate that even when targets are not collocated with grid points, adjusting the value of will lead to robust recovery of target locations. By using ≠ 0 we allow for robust recovery of targets regardless of the underlying grid used.

Zoom In Zoom Out Reset image size

Although we have found that imaging multiple point targets using F given in (6) leads to robust recovery of the target locations, we do not obtain the complex reflectivities. For that, we return to R given in (18). In (20) the matrix $\hat{R}D{\hat{R}}^{\text{T}}$ is symmetric, so its eigenvalues are real and its eigenvectors are orthogonal. However, the eigenvalues are not necessarily sign-definite nor sign-semidefinite. Hence, its diagonal factorization is

$$\begin{equation*}\hat{R}D{\hat{R}}^{\text{T}}=[\hat{W}\enspace \mathrm{exp}(\mathrm{i}\hat{{\Theta}}/2)]\hat{{\Sigma}}{[\hat{W}\enspace \mathrm{exp}(-\mathrm{i}\hat{{\Theta}}/2)]}^{\text{H}},\end{equation*}$$

where ${\hat{W}}^{\text{H}}\hat{W}=I$, $\hat{{\Sigma}}=\mathrm{diag}(\vert {\lambda }_{1}\vert ,\dots ,\vert {\lambda }_{N}\vert )$, and $\hat{{\Theta}}=\mathrm{diag}(\mathrm{arg}({\lambda }_{1}),\dots ,\mathrm{diag}({\lambda }_{N}))$. The measurement matrix B is then given by

$$\begin{equation*}B=[\hat{Q}\hat{W}\enspace \mathrm{exp}(\mathrm{i}\hat{{\Theta}}/2)]\enspace \hat{{\Sigma}}\enspace {[{\hat{Q}}^{\ast }\hat{W}\enspace \mathrm{exp}(-\mathrm{i}\hat{{\Theta}}/2)]}^{\text{H}}.\end{equation*}$$

In the SVD of the measurement matrix, B = UΣV^H, the first N columns of U are given by $\hat{Q}\hat{W}\enspace \mathrm{exp}(\mathrm{i}\hat{{\Theta}}/2)$, the first N columns of V are given by ${\hat{Q}}^{\ast }\hat{W}\enspace \mathrm{exp}(-\mathrm{i}\hat{{\Theta}}/2)$, and the first N diagonal entries of Σ correspond to the diagonal entries of $\hat{{\Sigma}}$. Suppose we evaluate R on the location of the nth target, y _n . Writing $\mathbf{g}({\boldsymbol{y}}_{n})={\Vert}\mathbf{g}({\boldsymbol{y}}_{n}){\Vert}\hat{Q}\hat{R}{\hat{\mathbf{e}}}_{n}$ with ${\hat{\mathbf{e}}}_{n}$ denoting the nth standard basis vector of ${\mathbb{R}}^{M}$, we find that

$$\begin{align*}{R}_{{\epsilon}}({\boldsymbol{y}}_{n})& ={\Vert}\mathbf{g}({\boldsymbol{y}}_{n}){{\Vert}}^{2}{\hat{\mathbf{e}}}_{n}^{\text{T}}{\hat{R}}^{\text{T}}{\hat{Q}}^{\text{T}}V{{\Sigma}}_{{\epsilon}}^{+}{U}^{\text{H}}\hat{Q}\hat{R}{\hat{\mathbf{e}}}_{n}\\ & ={\Vert}\mathbf{g}({\boldsymbol{y}}_{n}){{\Vert}}^{2}{\hat{\mathbf{e}}}_{n}^{\text{T}}{\hat{R}}^{\text{T}}\left[\hat{W}\enspace \mathrm{exp}(-\mathrm{i}\hat{{\Theta}}/2){{\Sigma}}^{-1}\enspace \mathrm{exp}(-\mathrm{i}\hat{{\Theta}}/2){\hat{W}}^{\text{H}}\right]\hat{R}{\hat{\mathbf{e}}}_{n}\\ & ={\Vert}\mathbf{g}({\boldsymbol{y}}_{n}){{\Vert}}^{2}{\hat{\mathbf{e}}}_{n}^{\text{T}}{\hat{R}}^{\text{T}}{\left[\hat{R}D{\hat{R}}^{\text{T}}\right]}^{-1}\hat{R}{\hat{\mathbf{e}}}_{n}\\ & ={\Vert}\mathbf{g}({\boldsymbol{y}}_{n}){{\Vert}}^{2}{\hat{\mathbf{e}}}_{n}^{\text{T}}{D}^{-1}{\hat{\mathbf{e}}}_{n}\\ & ={\rho }_{n}^{-1}.\end{align*}$$

We have used (19) to obtain the last step above. Therefore, evaluating 1/R with R ( x ) given in (18) at a point target location yields the complex reflectivity.

Using the results shown in figure 8, we recover the target locations. We evaluate 1/R at those locations and obtain approximations for the complex reflectivities. For target 1, the complex reflectivity is ρ₁ = 2.9096i and the estimated complex reflectivity is ρ₁ ≈ −0.000 506 76 + 2.9097i. For target 2, the complex reflectivity is ρ₂ = 3.5508i and the estimated complex reflectivity is ρ₂ ≈ 0.002 9675 + 3.5528i. For target 3, the complex reflectivity is ρ₃ = 2.2655i and the estimated complex reflectivity is ρ₃ ≈ −0.001 1034 + 2.2691i. Recall that the SNR = 5.5146  dB for this simulation.

To study the subwavelength resolution of this imaging method, we study imaging of two point targets in close proximity to one another. The measurement array of aperture ka = 80 included 41 elements sampling −40 ⩽ kx ⩽ 40 on kz = 0. The imaging region −1.2 ⩽ k(x − x₀) ⩽ 1.2 and −2.4 ⩽ k(z − L) ⩽ 2.4 is centered about (kx₀, kL) = (0, 100). We discretize this imaging region using mesh widths Δx = 0.1λL/a and Δz = 0.1λ(L/a)². Note that these mesh widths are related to the resolution estimates for active array imaging. We place two point targets at (k(x₁ − x₀), k(z₁ − L)) = (−0.45λL/a, −0.45λ(L/a)²) and at (k(x₂ − x₀), k(z₂ − L)) = (0.505λL/a, 0.505λ(L/a)²). Both of these locations are not aligned with grid points on the mesh.

In figure 9 we plot images of these two point targets formed through evaluation of 1/F ( x ) with F given in (6). Measurement noise has been added so that SNR = 4.5387  dB. The left plot in figure 9 is for = 10⁻⁶ and the right plot is for = 10⁻⁸. These results show that an appropriate choice of can separate two nearby scatterers. Note that the resulting quantities change with different choices of . Nonetheless, one can change the value of to accurately recover the location of targets. Then, one can evaluate 1/R ( x ) with R ( x ) given in (18) to recover the complex reflectivities.

Zoom In Zoom Out Reset image size

We now consider direct and inverse scattering by an extended target in 2D. To compute the measurement matrix, we use the method of fundamental solutions (MFS) introduced by Mathon and Johnston [20]. Then we use the Born approximation to study the quantitative subspace imaging method for an extended target. These results easily extend to 3D by using the appropriate Green's function for that problem.

5.1. MFS in 2D

Let Ω denote the support of the extended target, ∂Ω denote its boundary, and $\bar{{\Omega}}={\Omega}\cup \partial {\Omega}$ denote its closure. Additionally, let $E={\mathbb{R}}^{2} {\backslash}\bar{{\Omega}}$ denote the exterior. The wavenumber in E is k₀ and the wavenumber in Ω is k₁.

Let ${\boldsymbol{y}}_{n}^{\text{bdy}}$ for n = 1, ..., N denote points sampling the boundary of the extended target and let ν _n for n = 1, ..., N denote the corresponding unit outward normals. We introduce a small parameter ℓ and define ${\boldsymbol{y}}_{n}^{\text{int}}={\boldsymbol{y}}_{n}^{\text{bdy}}+\ell {\boldsymbol{\nu }}_{n}$ and ${\boldsymbol{y}}_{n}^{\text{ext}}={\boldsymbol{y}}_{n}^{\text{bdy}}-\ell {\boldsymbol{\nu }}_{n}$. Using these points we define the MFS approximation for the interior field according to

$$\begin{equation*}{u}^{\text{int}}(\boldsymbol{x})\approx \sum\limits _{n=1}^{N}{G}_{1}(\boldsymbol{x}-{\boldsymbol{y}}_{n}^{\text{int}}){c}_{n}^{\text{int}}\quad x\in {\Omega},\end{equation*}$$

with G₁ denoting Green's function for a medium with wavenumber k₁,

$$\begin{equation*}{G}_{1}(\boldsymbol{x}-{\boldsymbol{y}}_{n}^{\text{int}})=\frac{\mathrm{i}}{4}{H}_{0}^{(1)}({k}_{1}\vert \boldsymbol{x}-{\boldsymbol{y}}_{n}^{\text{int}}\vert ),\end{equation*}$$

and ${c}_{n}^{\text{int}}$ denoting coefficients to be determined. We define the MFS approximation for the scattered field according to

$$\begin{equation*}{u}^{\text{ext}}(\boldsymbol{x})\approx \sum\limits _{n=1}^{N}{G}_{0}(\boldsymbol{x}-{\boldsymbol{y}}_{n}^{\text{ext}}){c}_{n}^{\text{ext}}\quad x\in E,\end{equation*}$$

with G₀ denoting Green's function for a medium with wavenumber k₀,

$$\begin{equation*}{G}_{0}(\boldsymbol{x}-{\boldsymbol{y}}_{n}^{\text{int}})=\frac{\mathrm{i}}{4}{H}_{0}^{(1)}({k}_{0}\vert \boldsymbol{x}-{\boldsymbol{y}}_{n}^{\text{int}}\vert ),\end{equation*}$$

and ${c}_{n}^{\text{ext}}$ denoting coefficients to be determined. The field incident on the extended target is the signal emitted by array element i, which is G₀( x − x _i ). To determine coefficients ${c}_{n}^{\text{int}}$ and ${c}_{n}^{\text{ext}}$, we require that

$$\begin{equation}{G}_{0}({\boldsymbol{y}}_{n}^{\text{bdy}}-{\boldsymbol{x}}_{i})+{u}^{\text{ext}}({\boldsymbol{y}}_{n}^{\text{bdy}})={u}^{\text{int}}({\boldsymbol{y}}_{n}^{\text{bdy}})\quad n=1,\dots ,N,\end{equation} \tag{ 22 }$$

and

$$\begin{equation}{\partial }_{\nu }{G}_{0}({\boldsymbol{y}}_{n}^{\text{bdy}}-{\boldsymbol{x}}_{i})+{\partial }_{\nu }{u}^{\text{ext}}({\boldsymbol{y}}_{n}^{\text{bdy}})={\partial }_{\nu }{u}^{\text{int}}({\boldsymbol{y}}_{n}^{\text{bdy}})\quad n=1,\dots ,N,\end{equation} \tag{ 23 }$$

with ∂_ν denoting the normal derivative. Equations (22) and (23) give 2N linear equations for the 2N coefficients. Upon solution of this linear system, we compute entries of the measurement matrix through evaluation of

$$\begin{equation*}{b}_{ij}=\sum\limits _{n=1}^{N}{G}_{0}({\boldsymbol{x}}_{j}-{\boldsymbol{y}}_{n}^{\text{ext}}){c}_{n}^{\text{ext}},\quad j=1,\dots ,M.\end{equation*}$$

This computation is repeated for each i = 1, ..., M to compute the full measurement matrix.

5.2. Quantitative signal subspace imaging of extended targets

Suppose that the extended target is centered at position x ₀ = (x₀, L) and let Ω₀ denote its support when centered at the origin. Let ɛ_r denote the relative dielectric permittivity so that k₁ = k₀ ɛ_r. According to the Born approximation, the entries of the measurement matrix are given by

$$\begin{equation*}{b}_{ij}\approx {\int }_{{{\Omega}}_{0}}{G}_{0}({\boldsymbol{x}}_{i}-{\boldsymbol{x}}_{0}-\boldsymbol{y})\left[{k}^{2}\left({\varepsilon }_{\text{r}}(\enspace \boldsymbol{y})-1\right)\right]{G}_{0}(\enspace \boldsymbol{y}-{\boldsymbol{x}}_{j}+{\boldsymbol{x}}_{0})\mathrm{d}\boldsymbol{y}.\end{equation*}$$

Let y = (ξ, ζ). Assuming that kL ≫ 1, we use the leading behavior, ${H}_{0}^{(1)}(z)\sim \sqrt{2/(\pi z)}{\text{e}}^{-\mathrm{i}\pi /4}{\text{e}}^{\mathrm{i}z}$ for z → ∞ and

$$\begin{equation*}k\vert {\mathbf{x}}_{i}-{\mathbf{x}}_{0}-\mathbf{y}\vert \sim k\left(L-\frac{{x}_{0}{x}_{i}}{L}+\frac{{x}_{i}^{2}}{2L}\right)+k\zeta \left(1-\frac{{x}_{i}^{2}}{2L}\right)-k\xi \frac{{x}_{i}}{L},\quad kL\to \infty ,\end{equation*}$$

to obtain

$$\begin{equation}{G}_{0}({\mathbf{x}}_{i}-{\mathbf{x}}_{0}-\mathbf{y})\sim \frac{{\text{e}}^{-\mathrm{i}\pi /4}}{2\sqrt{2\pi kL}}\enspace {\text{e}}^{\mathrm{i}k(L+{x}_{i}^{2}/(2L)-{x}_{0}{x}_{i}/L)}\enspace {\text{e}}^{\mathrm{i}k\zeta (1-{x}_{i}^{2}/(2{L}^{2}))}\enspace {\text{e}}^{-\mathrm{i}k\xi {x}_{i}/L}.\end{equation} \tag{ 24 }$$

Using (24), we find that

$$\begin{align*}{b}_{ij}& \sim {\text{e}}^{\mathrm{i}k(L+{x}_{i}^{2}/(2L)-{x}_{0}{x}_{i}/L)}\left[\frac{-\mathrm{i}}{8\pi kL}v\left(\frac{k({x}_{i}+{x}_{j})}{L},2k-\frac{k({x}_{i}^{2}+{x}_{j}^{2})}{2{L}^{2}}\right)\right]{\text{e}}^{\mathrm{i}k(L+{x}_{j}^{2}/(2L)-{x}_{0}{x}_{j}/L)},\end{align*}$$

with

$$\begin{equation*}v({k}_{x},{k}_{z})={k}^{2}{\int }_{{{\Omega}}_{0}}\left[{\varepsilon }_{\text{r}}(\xi ,\zeta )-1\right]{\text{e}}^{-\mathrm{i}{k}_{x}\xi }\enspace {\text{e}}^{\mathrm{i}{k}_{z}\zeta }\enspace \mathrm{d}\xi \enspace \mathrm{d}\zeta .\end{equation*}$$

Let V denote the M × M matrix whose entries are given by

$$\begin{equation}{V}_{ij}=\frac{-\mathrm{i}}{8\pi kL}v\left(\frac{k({x}_{i}+{x}_{j})}{L},2k-\frac{k({x}_{i}^{2}+{x}_{j}^{2})}{2{L}^{2}}\right),\quad i,j=1,\dots ,M.\end{equation} \tag{ 25 }$$

It follows that the measurement matrix B is

$$\begin{equation*}B={Q}_{0}V{Q}_{0}^{\text{T}},\end{equation*}$$

with Q₀ given in (8). The matrix V contains information regarding the support of the extended target along with its material properties. It is symmetric since we can interchange x_i and x_j in (25) and the entries remain unchanged. Thus, its diagonal factorization is V = WΛW^H. However, it is not Hermitian, so the entries of Λ are not necessarily sign-definite nor sign-semidefinite. For this reason we write

$$\begin{equation*}{\Lambda}=\mathrm{exp}(\mathrm{i}{\Theta}/2)\hat{{\Sigma}}\mathrm{exp}(\mathrm{i}{\Theta}/2),\end{equation*}$$

with

$$\begin{equation*}{\Theta}=\mathrm{diag}(\mathrm{arg}({\lambda }_{1})-\pi /4,\dots ,\mathrm{arg}({\lambda }_{M})-\pi /4),\end{equation*}$$

and

$$\begin{equation*}\hat{{\Sigma}}=\frac{1}{8\pi kL}\mathrm{diag}(\vert {\lambda }_{1}\vert ,\dots ,\vert {\lambda }_{M}\vert ).\end{equation*}$$

Thus, we rewrite our model for the measurement matrix as

$$\begin{equation*}B={Q}_{0}W\enspace \mathrm{exp}(\mathrm{i}{\Theta}/2)\hat{{\Sigma}}\mathrm{exp}(\mathrm{i}{\Theta}/2){W}^{\text{H}}{Q}_{0}^{\text{T}}.\end{equation*}$$

In the SVD of the measurement matrix, B = UΣV^H, the columns of U are given by ${Q}_{0}W\enspace \mathrm{exp}(\mathrm{i}\hat{{\Theta}}/2)$, the columns of V are given by ${Q}_{0}^{\ast }W\enspace \mathrm{exp}(-\mathrm{i}\hat{{\Theta}}/2)$, and the diagonal entries of Σ correspond to the diagonal entries of $\hat{{\Sigma}}$. Just as we have done for point targets, we form an image of the extended target by plotting 1/F with F given in (6).

Note that in (25), the first argument of v is O(1/(kL)) and the second argument of v is 2k + O(1/(kL)²). Therefore, in the asymptotic limit kL → +∞, we find that

$$\begin{equation}V\sim \sum\limits _{n=0}^{\infty }{(kL)}^{-n}{V}_{n},\quad kL\to +\infty .\end{equation} \tag{ 26 }$$

Let

$$\begin{equation}{\nu }_{p,q}={\int }_{{{\Omega}}_{0}}{(-\mathrm{i}\xi )}^{p}{(\mathrm{i}\zeta )}^{q}V(\xi ,\zeta ){\text{e}}^{\mathrm{i}2k\zeta }\enspace \mathrm{d}\xi \enspace \mathrm{d}\zeta .\end{equation} \tag{ 27 }$$

It follows that ${V}_{0}={\nu }_{0,0}{\mathbf{1}}_{M}{\mathbf{1}}_{M}^{\text{T}}$. Using this leading behavior of asymptotic expansion (26), we find that the principal eigenvalue λ and corresponding eigenvector w satisfy,

$$\begin{equation*}\lambda \mathbf{w}=V\mathbf{w}\sim {V}_{0}\mathbf{w}={\nu }_{0,0}{\mathbf{1}}_{M}{\mathbf{1}}_{M}^{\text{T}}\mathbf{w},\end{equation*}$$

from which we find that λ ∼ Mν_0,0 and w ∼ M^−1/2 1_M = w₀. Because of asymptotic expansion (26), the other eigenvalues will be asymptotically smaller than this principal eigenvalue [6]. Let ${\nu }_{0,0}=\vert {\nu }_{0,0}\vert {\text{e}}^{\mathrm{i}{\theta }_{0}}$. By substituting this leading behavior for V into expression for B, we find

$$\begin{equation}B\sim ({\text{e}}^{\mathrm{i}{\theta }_{0}/2}\enspace {Q}_{0}{\mathbf{w}}_{0})\frac{M\vert {\nu }_{0,0}\vert }{8\pi kL}{({\text{e}}^{-\mathrm{i}{\theta }_{0}/2}\enspace {Q}_{0}^{\ast }{\mathbf{w}}_{0})}^{\text{H}},\quad kL\to \infty .\end{equation} \tag{ 28 }$$

Similar to (9) which is for a point target in 3D, we see that the leading behavior given in (28) corresponds to a point target in 2D whose reflectivity has magnitude |ν_0,0|. We find from (28) that the leading behavior of the first singular value of B is σ₁ ∼ M|ν_0,0|/(8πkL) as kL → +∞. After substituting (27) with p = q = 0, we find that

$$\begin{equation*}{\sigma }_{1}\sim \frac{M}{8\pi kL}\left\vert {\int }_{{{\Omega}}_{0}}{k}^{2}\left[{\varepsilon }_{\text{r}}(\xi ,\zeta )-1\right]{\text{e}}^{\mathrm{i}2k\zeta }\enspace \mathrm{d}\xi \enspace \mathrm{d}\zeta \right\vert ,\quad kL\to +\infty .\end{equation*}$$

Consider an extended target composed of a uniform medium so that ɛ_r is constant. For that case, we have

$$\begin{equation*}{k}^{2}\vert {\varepsilon }_{\text{r}}-1\vert \left\vert {\int }_{{{\Omega}}_{0}}{\text{e}}^{\mathrm{i}2k\zeta }\enspace \mathrm{d}\xi \enspace \mathrm{d}\zeta \right\vert \leqslant {k}^{2}\vert {\varepsilon }_{\text{r}}-1\vert \vert {{\Omega}}_{0}\vert ,\end{equation*}$$

where |Ω₀| denotes the area of the extended target. Therefore, we find

$$\begin{equation*}{\sigma }_{1}\leqslant \frac{M\vert {{\Omega}}_{0}\vert }{8\pi kL}{k}^{2}\vert {\varepsilon }_{\text{r}}-1\vert ,\quad kL\to +\infty .\end{equation*}$$

We can compute exact results for some examples. For the rectangular box extended target defined as Ω₀ = [−α_x /2, α_x /2] × [α_z /2, α_z /2], we have

$$\begin{equation*}{\sigma }_{1}^{\;\text{box}}\sim \frac{M}{8\pi kL}{k}^{2}{\alpha }_{x}{\alpha }_{z}\vert {\varepsilon }_{\text{r}}-1\vert \vert \text{sinc}\left(k{\alpha }_{z}\right)\vert ,\quad kL\to \infty .\end{equation*}$$

For the circular disk extended target, Ω₀ = {r < α₀}, we have

$$\begin{equation}{\sigma }_{1}^{\;\text{circle}}\sim \frac{M}{8\pi kL}\pi {(k{\alpha }_{0})}^{2}\vert {\varepsilon }_{\text{r}}-1\vert \frac{2{J}_{1}(2k{\alpha }_{0})}{2k{\alpha }_{0}},\quad kL\to \infty .\end{equation} \tag{ 29 }$$

The asymptotic behavior for the first singular value gives useful quantitative information. Provided that a signal is available in measurements, this singular value will be reliably accurate. On the other hand, the other singular values are more subject to corruption due to measurement noise which, in turn, will affect the reconstruction of the shape of the object.

5.3. Results for a circular disk target

We compute the measurement matrix for a circular disk extended target that is centered at (kx₀, kL) = (0, 100) with relative dielectric permittivity ɛ_r = 1.44 and radius kα₀ = 0.25 (left) and kα₀ = 1 (right). The measurement array of aperture ka = 80 included 41 elements sampling −40 ⩽ kx ⩽ 40 on kz = 0. For the MFS computation we used ${\boldsymbol{y}}_{n}^{\text{bdy}}=\left({x}_{0}+{\alpha }_{0}\enspace \mathrm{sin}(2\pi (n-1)/N),L+{\alpha }_{0}\enspace \mathrm{cos}(2\pi (n-1)/N)\right)$ for n = 1, ..., N with N = 64 and ℓ = 0.12α₀. The image formed by plotting 1/F (x, z) for = 10⁻¹⁵ on a log₁₀-scale is shown in figure 10. Noise has been added to the measurements so that SNR = 10.5172  dB. These imaging results show that 1/F (x, z) accurately recovers the location and support of the extended target.

Zoom In Zoom Out Reset image size

For the circular disk with kα₀ = 0.25, the computed first singular value is σ₁ = 0.001 3440 and the asymptotic approximation given in (29) is σ₁ = 0.001 3658. For the circular disk with kα₀ = 1, the computed first singular value is σ₁ = 0.012 737 and the asymptotic approximation is σ₁ = 0.013 005. These results show that the asymptotic approximation is accurate over a broad range of target sizes.

Figure 11 shows results for a circular disk target with kα₀ = 0.5 with different amounts of measurement noise. The upper left plot has SNR = ∞, the upper right plot has SNR = 10.5207  dB, the lower left plot has SNR = 5.5247  dB, and the lower right plot has SNR = 2.5211  dB. The value of is adapted to the measurement SNR so as to separate the signal from the noise subspace with ${{\Sigma}}_{{\epsilon}}^{+}$ defined as in (21). These results show that as SNR decreases, the image of the circular disk target degrades.

Zoom In Zoom Out Reset image size

For the SNR = 2.5211  dB result, we find that range resolution is nearly lost completely. Despite the fact that images degrade severely as the SNR decreases, we find that the quantitative estimate remains robust for different SNRs. For this circular disk with kα₀ = 0.5, the computed first singular value is σ₁ = 0.005 1974 and remains the same for the different SNRs. The asymptotic approximation is σ₁ = 0.004 9616. Using the computed σ₁ and the exact value for kα₀ in (29), we solve for ɛ_r and obtain ɛ_r = 1.4609 corresponding to a 1.5% error compared to the true value of ɛ_r = 1.44.

We have developed a quantitative imaging method for single-frequency array imaging. This method extends MUSIC which uses the noise subspace to reconstruct the support of targets by including also the signal subspace which provides quantitative information about the targets. This answers the open question raised in [11] of how to make use of the signal subspace in MUSIC.

We have applied this imaging method to a point target, multiple point targets, and an extended target. For a single point target, we have presented a resolution analysis. This resolution analysis provides valuable insight into the subwavelength resolution capabilities of this imaging method. It explicitly shows the difference between cross-range and range resolutions. Moreover, it shows the surprising result in which noise naturally regularizes the imaging method making it more robust and effective, especially when the target is not located on a grid point of the mesh discretizing the imaging region. For multiple point targets, we show that one can artificially change the noise level by changing the singular values associated with the singular vectors spanning the noise subspace to simultaneously change the resolution and regularization of the imaging method. For an extended target, we show that the first singular value is related to the product of the volume and relative dielectric permittivity of the target.

In our results we have found that the imaging method presented here is remarkably effective at determining the support and complex reflectivities of point targets. Moreover, the results verify that regularization-by-noise occurs for this imaging method. Although we have only discussed single-frequency array imaging, the method is general and can be applied to other imaging systems for which MUSIC has already proved to be effective. Because this imaging method is just an elementary extension of MUSIC, it should be useful for a variety of quantitative imaging problems.

A D Kim and C Tsogka are supported by the Air Force Office of Scientific Research (FA9550-21-1-0196). A D Kim is also supported by the National Science Foundation (DMS-1840265). Pedro González-Rodríguez is supported by the Spanish Ministerio de Ciencia e Innovación (PID2020-115088RB-100).

The data that support the findings of this study are available upon reasonable request from the authors.