Quantitative subsurface imaging in strongly scattering media

Pedro González-Rodríguez; Arnold D. Kim; Miguel Moscoso; Chrysoula Tsogka

doi:10.1364/OE.26.027346

1. Introduction

Subsurface imaging in strongly scattering media has several applications in biomedical optics and geophysical remote sensing. This imaging problem is severely ill-posed due to the inherent lack of measurements available and inherent mathematical instabilities associated with the inverse problem. The main challenge with this problem is extracting useful information about targets in the medium from measurements of backscattered light. Because light is strongly scattered by the medium, effectively none of the backscattered light corresponds to direct information about targets. Consequently, there is an overall loss of data fidelity which, in turn, affects image resolution and stability to noise.

Much of the work done on subsurface imaging has to do with processing measurements and filtering out unwanted parts to increase data fidelity and overall signal-to-noise ratio. For example, polarization gating has been shown to distinguish between less scattered light and diffusely scattered light [1–3]. Other examples include using the time [4–6], frequency [7], and angle characteristics of backscattered light [8]. All these techniques discriminate between weakly and strongly scattered light giving better image resolution. However, there is less known about how well a particular imaging algorithm will perform in reconstructing subsurface images from backscattered light measurements. For example, it may be the case that an appropriate imaging method may alleviate the need for sophisticated technology requirements or may better identify the issues surrounding limitations inherent in an imaging system.

In this work, we introduce a method for subsurface imaging of targets at shallow depths. This imaging method uses time-dependent backscattered light measurements to reconstruct point-like targets in the medium. Rather than reconstructing images using multiple source-detector pairs, we show that collecting a series of one-dimensional reconstructions at null source-detector distances works better. The key point here is that this method is able to reconstruct quantitatively accurate, high depth resolution images using only diffusive light at moderate sampling rates. This imaging method has an inherent depth limitation due to the exponential decay of the signal due to scattering and absorption by the surrounding medium. The method proposed here reconstructs targets at depths of the order of a scattering mean-free path as OCT does [9], but uses only diffusive light. Consequently, this imaging method opens opportunities for developing novel imaging modalities. Imaging deeper into the tissue, beyond the ballistic limit, can be enhanced by including a few angle-resolved measurements as it was shown in [8]. To keep the measurement process as simple as possible, we use angle-averaged measurements here to form the images.

The method has two stages. In the first stage, the location of the targets is determined using the MUSIC algorithm. MUSIC is an abbreviation of MUltiple SIgnal Classification and the algorithm was first introduced by Schmidt in [10]. It has been used successfully in signal processing [11] and imaging [12, 13]. We also refer to [14] for an analysis of MUSIC for single-snapshot spectral estimation. More recently, a modified MUSIC algorithm was proposed for diffuse optical tomography by Dileep et al. [15]. The algorithm proposed here is different from that work in that it is adapted to null source-detector distances. In this case, only a single data vector is available for imaging at each measurement location. This is a challenging setup for MUSIC that relies on multiple measurement vectors collected at several detectors. We show, however, that using a Prony-type idea [22] we can form a suitable data matrix from which the locations of the targets can be determined with MUSIC. In the second stage, those target locations are used to simplify the determination of the optical quantities associated to them. This second stage is done through the solution of a relatively simple, low dimensional least-squares problem.

One key to the effectiveness of this imaging algorithm is due to the use of an asymptotic model for backscattered light measurements developed by Rohde and Kim [16]. This asymptotic theory introduces a slightly modified diffusion approximation that is valid at null source-detector distances and at shallow penetration depths. It is this diffusion approximation that is used in the two-stage imaging method described above.

The remainder of this paper is as follows. In Section 2, we give a description of the asymptotic model used in the reconstruction method. Included in this section are some initial results that use a standard ℓ₂ method for image reconstruction. Those poor results motivate the use of a more sophisticated imaging algorithm that is described in Section 3. In Section 4, we state the limitations of this imaging method and the sources for error. We show numerical results of the new imaging method in Section 5. Section 6 contains our conclusions.

2. Modeling measurements

We use radiative transfer theory to describe light backscattered by a half space containing small, point-like inhomogeneities located superficially below the surface of a strongly scattering medium. Using the asymptotic solution of the radiative transfer equation in strongly scattering media, we derive a model for measurements. This model is further simplified by linearizing the effect of those inhomogeneities on the data using the Born approximation. We show that the minimum ℓ₂-norm solution for image reconstruction is ineffective for this problem thereby motivating the need for a more sophisticated imaging method.

2.1. Radiative transfer

The specific intensity I(ŝ, r, t) gives the power flowing in direction ŝ at position r, and time t. The radiative transfer equation

\frac{1}{c} \frac{\partial I}{\partial t} + \hat{s} \cdot \nabla I + μ_{a} I + μ_{s} I = μ_{s} \int_{4 π} f (\hat{s} \cdot {\hat{s}}^{'}) I ({\hat{s}}^{'}, r, t) d {\hat{s}}^{'}

governs I in a medium that absorbs and scatters light. Here, c is the speed of light in the medium, μ_a is the absorption coefficient, and μ_s is the scattering coefficient. The scattering phase function is denoted by f. It gives the fraction of light incident in direction ŝ′ that is scattered in direction ŝ. We have assumed that f is spherically symmetric so that it depends only on ŝ · ŝ′.

Consider the half space z > 0 with boundary z = 0. The refractive index inside the half space is different from that outside of it. We seek to solve (1) subject to the initial condition

{I |}_{t = 0} = 0,

and boundary condition

{I |}_{z = 0} = T_{0} δ (\hat{s} - \hat{z}) h (x, y, t) + R [I]

prescribed over the hemisphere of directions that point into the medium corresponding to ŝ · ẑ > 0. Here, T₀ is the transmission coefficient for the pulsed beam, denoted by h (x, y, t), incident normally on the boundary. In addition, R[I] is the internal reflection of light incident on the boundary from inside the half space due to the refractive index mismatch. The initial-boundary value problem comprised of (1) subject to initial condition (2), boundary condition (3), and requiring that I → 0 as z → ∞ gives a complete mathematical description of the direct problem.

For measurements, we consider the time-resolved diffuse reflectance

r (x, y, t) = - \int_{NA} T^{out} [I] (\hat{s}, x, y, 0, t) \hat{s} \cdot \hat{z} d \hat{s} .

Here, NA denotes the numerical aperture of the detector. It is a subset of directions satisfying ŝ · ẑ < 0, corresponding to light exiting the half space. The operation T^out[I] denotes the light incident on z = 0 from within the half space that is transmitted out of the half space across the index mismatched boundary. These measurements give the spatial-temporal distribution of light backscattered by the medium. Because we only consider backscattered light, the inverse problem that forms the images is severely ill-posed.

2.2. Asymptotic solution in strongly scattering media

The strong scattering limit corresponds to when multiple scattering in the medium is dominant and absorption is weak so that μ_s ≫ μ_a. It is well understood that the diffusion approximation accurately describes light that has propagated deep into a strongly scattering medium [17,18]. This approximation states that $u (r, t) ~ \int_{4 π} I (\hat{s}, r, t) d \hat{s}$ satisfies the diffusion equation

\frac{1}{c} \frac{\partial u}{\partial t} + μ_{a} u - \nabla \cdot (D \nabla u) = 0,

where D denotes the diffusion coefficient.

The challenge in using the diffusion approximation is that approximate boundary conditions and sources must be used. For this reason, the conventional wisdom is that the diffusion approximation is not valid near boundaries or sources. Rohde and Kim [16] have recently developed an asymptotic theory for the radiative transfer equation in the strong scattering limit. This theory includes a so-called boundary layer solution that corrects the error made by the diffusion approximation near boundaries and sources. Using this asymptotic theory, we determine that the diffusion coefficient is D = 1/[3μ_s(1 − g)] with g denoting the anisotropy factor or mean cosine of the scattering angle. Very often, μ_a is included in the diffusion coefficient so that D = 1/[3μ_s(1 − g) + μ_a], but because μ_s ≫ μ_a is assumed here, its role in D is negligible in this asymptotic theory. The asymptotic approximation for the time-resolved diffuse reflectance given by (4) is

r ~ r_{0} h + r_{1} \partial_{n} {u |}_{z = 0},

with h denoting the same pulsed beam introduced in (3), and ∂_nu|_z=0 denoting the normal derivative of the solution of (5) (evaluated on z = 0) subject to the initial condition

{u |}_{t = 0} = 0,

and boundary condition

{u |}_{z = 0} = c_{0} h (x, y, t) .

The coefficients, r₀, r₁, and c₀ are constants that depend on fundamental quantities in the problem.

The main difference in this diffusion approximation, from what is typically used, is the Dirichlet boundary condition (8) instead of a Robin-type boundary condition, which is derived through a boundary layer that fulfills boundary conditions and sources in the radiative transfer equation. Effectively, the boundary layer solution relieves the diffusion approximation from having to completely resolve boundary conditions and sources in the problem on its own, thereby allowing for this simpler boundary condition (see [16] for details).

According to (6), each measurement is a linear combination of the known pulsed beam, and the normal derivative of the solution of the diffusion equation evaluated at the boundary. Consequently, the term proportional to ∂_nu|_z=0 is the only one that carries information about the interior of the domain. Suppose that the background optical properties of the medium are known. For that case, we can determine the coefficients r₀ and r₁, and since h is known we can remove its contribution from the measurements. Even if the coefficients are not known exactly, one can effectively remove the contribution due to h from the measurements using r−(〈h, r〉/〈h, h〉) h, where 〈·, ·〉 denotes the inner product. This calculation removes the projection of the measurements onto h, which effectively leaves only the portion of the data due to ∂_nu|_z=0. Thus, in what follows, we take as measurements

\tilde{r} = \partial_{n} {u |}_{z = 0},

with u satisfying (5) subject to initial condition (7) and boundary condition (8).

In what follows, we study perturbations of D and develop a method for recovering those perturbations from measurements. One can consider perturbations of μ_a instead. However, we discuss below in Section 4 that this method is not able to recover both D and μ_a simulatenously because it is assumed throughout that μ_s ≫ μ_a. It follows that any perturbations of μ_a are negligible compared to those of D, especially in the presence of noise.

2.3. Born approximation

We will consider measurements in the frequency domain corresponding to the Fourier transform

U (r, ω) = \frac{1}{2 π} \int u (r, t) e^{- i ω t} d t,

where U satisfies the equation

\nabla \cdot (D \nabla U) - (μ_{a} + i ω / c) U = 0

in the half space z > 0, subject to

{U |}_{z = 0} = c_{0} H (x, y, ω),

with H(x, y, ω) denoting the Fourier transform of the pulsed beam, h(x, y, t). Let U₀ denote the solution to the homogeneous problem, with constant background diffusion coefficient D₀,

\nabla^{2} U_{0} - k^{2} U_{0} = 0

in z > 0 and subject to boundary condition (12), where k² = (μ_a + iω/c)/D₀.

We model inhomogeneities in the medium, denoted by δD, as a perturbation to the constant background diffusion coefficient, so D = D₀ + δD. Substituting this into (11) yields, after rearranging terms,

\nabla^{2} (U - U_{0}) - k^{2} (U - U_{0}) = - D_{0}^{- 1} \nabla \cdot (δ D \nabla U),

subject to boundary condition

{(U - U_{0}) |}_{z = 0} = 0 .

The half-space Green’s function G for this problem satisfies

\nabla^{2} G - k^{2} G = - D_{0}^{- 1} δ (r - r^{'})

in z > 0, subject to

{G |}_{z = 0} = 0 .

It is given by

G (r, r^{'}) = \frac{e^{- k r}}{4 π D_{0} r} - \frac{e^{- k r^{*}}}{4 π D_{0} r^{*}},

with

r = \sqrt{{(x - x^{'})}^{2} + {(y - y^{'})}^{2} + {(z - z^{'})}^{2}}

and

r^{*} = \sqrt{{(x - x^{'})}^{2} + {(y - y^{'})}^{2} + {(z + z^{'})}^{2}}

. In terms of G, the solution of (14) is given by

U - U_{0} = \int_{z^{'} > 0} G (r, r^{'}) [\nabla \cdot (δ D \nabla U)] d r^{'} .

The operations and quantities within the brackets in (19) are to be evaluated with respect to r′.

We now apply the Born approximation, which gives

U \approx U_{0} + \int_{z^{'} > 0} G (r, r^{'}) [\nabla \cdot (δ D \nabla U_{0})] d r^{'} .

Eq. (20) gives an explicit representation for the perturbation made by δD on the solution U in terms of the half-space Green’s function G(r, r′). Since G is the only function of r in (20), we find that the perturbed measurements in the frequency domain, R̃ = ∂_zU|_z=0, are given by

\tilde{R} \approx \partial_{z} {U_{0} |}_{z = 0} - \int_{z^{'} > 0} [\nabla G_{z} (r, r^{'}) \cdot \nabla U_{0}] δ D d r^{'},

where we have used partial integration to obtain this result. In (21), G_z denotes the partial derivative of G with respect to z.

We now write U₀ in terms of the surface Green’s function G̃ satisfying

\nabla^{2} \tilde{G} - k^{2} \tilde{G} = 0

in z > 0, subject to boundary condition

{\tilde{G} |}_{z = 0} = δ (x - x^{'}) δ (y - y^{'}) .

It follows that U₀ is given by the convolution U₀ = G̃ * c₀H, which is an integral over the z = 0 plane. Using Green’s identities, we find that G̃ is related to the half-space Green’s function G given in (18) according to G̃ = −D₀G_z|_z′=0⁺.

Finally, we assume an infinitely narrow beam incident on the boundary at point (x_s, y_s, 0). Thus, U₀ is the solution to the homogeneous problem (13) subject to boundary condition U₀|_z=0 = δ(x − x_s)δ(y − y_s). In other words, U₀ = G̃ in (21) and the measurements become

\tilde{R} \approx \partial_{z} {U_{0} |}_{z = 0} + D_{0} \int_{z^{'} > 0} [\nabla G_{z} (r_{d}, r^{'}) \cdot \nabla G_{z} (r^{'}, r_{s})] δ D d r^{'} .

Note that R̃ = R̃(r_d, r_s) depends on the source and detector positions r_s = (x_s, y_s, 0) and r_d = (x_d, y_d, 0), respectively.

2.4. Linear system

Our aim is to reconstruct a small point-like diffusion perturbation δD assuming that we can measure the difference ∂_z(U − U₀)|_z=0 due to the perturbation for multiple sources and detectors. Under the Born approximation the measurements given by (24) are linear in the unknown perturbation δD. Thus, the inverse problem reduces to solving a system of linear equations. We assume that the images are sparse, meaning that only a small number M of point-like diffusion perturbations are present.

Specifically, if we discretize the region of interest using N grid points r′_n we can write (24) as the linear system

𝒜_{s} x = b_{s}

for a fixed source r_s. In (25), 𝒜_s is the model matrix that maps the unknown vector x = δD ∈ ℝ^N to the data vector b_s ∈ ℂ^D. If we use N_d detectors and n_f frequencies to form the images, then D = N_d · n_f. The components (b_s)_j of b_s correspond to the difference ∂_z(U − U₀)|_z=0/D₀ measured at detector i at frequency ω_l [j ≡ j(i, l) = i + (l − 1)N_d, with i = 1, . . ., N_d, and l = 1, . . ., n_f, is double indexed]. Finally, the elements

𝒜_{s} (i + (l - 1) N_{d}, n) = \nabla G_{z} (r_{i}, r_{n}; ω_{l}) \cdot \nabla G_{z} (r_{n}, r_{s}; ω_{l})

of the model matrix come from the discretization of the integral in (24). In (26), i = 1, . . ., N_d, l = 1, . . ., n_f, and n = 1, . . ., N.

2.5. Initial imaging results

Typically, the linear system (25) is underdetermined and, hence, there are infinitely many δD distributions that match the data vector b_s. If there is only one source of illumination, the usual method to form an image is to pick the minimum ℓ₂-norm solution. If data from multiple sources are available, one can average over the different solutions obtained for each source. In the top row of Fig. 1 we show the distribution of targets δD that we seek to reconstruct. They take values between 0.1cm and 0.8cm. Note that in the right image there are two nearby targets located at the same cross-range. The bottom row of Fig. 1 shows the reconstructions for these three cases choosing the minimum ℓ₂-norm solution. The background medium characteristics are μ_a0 = 0.05 cm −1, μ_s0 = 7.5 cm −1, and g = 0.8. In the reconstructions we have used N_s = 21 sources and N_d = 21 detectors located on the surface. Both are 0.3 cm apart. We have also used information from n_f = 16 frequencies, with central frequency 10 GHz and sampling rate 0.13 GHz. This means a 2 GHz bandwidth. In all the experiments shown in the paper the data is contaminated by additive noise so that the signal-to-noise ratio (SNR) ranges between 30dB and 35dB. These results show that the usual minimum ℓ₂-norm approach fails to provide good images.

Fig. 1 Top row: Exact distributions of perturbations δD to be reconstructed. Bottom row: minimum ℓ₂-norm reconstructions using all source-detector pairs with a 2 GHz system bandwidth.

Download Full Size | PDF

Given the large distance between detectors, we regularize the problem by considering only measurements corresponding to null source-detector distances. Thus, we use only data received at the detector that is collocated with the source. Short source-detector distances in diffuse optical imaging have been proposed in [6,19,20] to improve the images. The reconstructions in Fig. 2 are obtained by solving N_d one-dimensional problems where, at each measurement location, we seek to recover the perturbation δD as a function of depth only. This strategy needs less computational resources (and less measurements) but does not improve the quality of the images. Indeed by comparing the images in the bottom row of Fig. 1 with those in Fig. 2 we observe that the minimum ℓ₂-norm reconstructions are similar whether we use data from multiple sources and multiple detectors or only null source-detector separation data.

Fig. 2 Minimum ℓ₂-norm reconstructions using collocated sources and detectors for the same cases as shown in Fig. 1 with a 2 GHz system bandwidth.

Download Full Size | PDF

3. Image reconstruction method

The starting point of the imaging method proposed here is the remark pointed out above; null source-detector separation data are sufficient to image the obstacles or inhomogeneities when cross-range resolution of the order of the scattering mean free path or larger is sufficient. In this case, given the relatively large distance between detectors, these data are optimal in the sense that they carry most of the information needed to find out how deep the inhomogeneities are. The proposed imaging method consists of two steps. First, we recover the depth or range of the obstacles below each measurement location using the MUSIC algorithm with null source-detector separation data. By scanning the sample along the cross-range direction we determine the support of all the obstacles. Once the support is found, we estimate in the second step the values of the inhomogeneities using least-squares restricted to the recovered support. This two step strategy performs remarkably well as it is shown in Section 5 through numerical experiments.

Next, we describe in more detail the reduced, one dimensional inverse problem for a source-detector pair at a fixed location. We explain how to apply MUSIC to find the depth or range of the unknown inhomogeneities. The method is similar to the one given in [21] for the one-dimensional reflectivity imaging problem for the wave equation. The cross-range dependence of the inhomogeneities is recovered by moving the source-detector pair across the sample.

For a fixed source-detector location, the multifrequency data vector b = [b₁, b₂, . . ., b_{n_f}] has components

b_{l} = \int_{z > 0} \nabla G_{z} (r_{d}, r; ω_{l}) \cdot \nabla G_{z} (r, r_{s}; ω_{l}) δ D (r) d r,

with r_d = r_s since we use null source-detector distances. We consider measurements at n_f = 2S equidistant frequencies ω_l = ω₀ ± lΔω, l = 1, . . ., S, within a small bandwidth B = S Δω ≪ ω₀.

Given a one dimensional grid with N points in depth z_n = z₀ + (n − 1)Δz, n = 1, . . ., N, we denote by δD_n the unknown perturbation in the diffusion coefficient at position z_n. Thus, discretizing (27), our one-dimensional discrete model is given by

b_{l} \approx \sum_{n = 1}^{N} G_{z z} (z_{0}; z_{n}; ω_{l}) G_{z z} (z_{n}, z_{0}; ω_{l}) δ D_{n} Δ z,

for l = 1, . . ., n_f. In (28), z₀ = 0 and G_zz(z₀, z_n, ω_l) is the second-order partial derivative of the one dimensional Green’s function

G (z, z^{'}; ω_{l}) = \frac{e^{- k_{l} | z - z^{'} |}}{4 π D_{0} | z - z^{'} |} - \frac{e^{- k_{l} | z + z^{'} |}}{4 π D_{0} | z + z^{'} |}

with respect to z evaluated at z = z_n. In (29),

k_{l} = \sqrt{(μ_{a 0} + i ω_{l} / c) / D_{0}}

.

We seek to recover the vector δD = [δD₁, δD₂, . . ., δD_N] from the multifrequency data vector b = [b₁, b₂, . . ., b_{n_f}], with n_f = 2S ≪ N. We assume that the vector δD is sparse, meaning that only a few components are different from zero. Eq. (28) can be re-written as

𝒜 x = b,

with x = δD, and

𝒜 (l, n) = G_{z z}^{2} (z_{0}, z_{n}; ω_{l})

for n = 1, . . ., N, and l = 1, . . ., n_f. Since 𝒜 ∈ ℂ^{n_f ×N} with n_f ≪ N, the linear system (30) is underdetermined meaning that there infinitely many obstacle configurations that match the data. If we pick the one whose ℓ₂-norm is minimal the results are bad as it was is illustrated in Fig. 2.

The key idea so as to form better images is to determine first the support of the unknown vector x by using MUSIC. The MUSIC algorithm, however, needs multiple measurement vectors to compute the left nullspace space of a suitable data matrix, usually referred to as the noise subspace. Since we only have one data vector b ∈ ℂ^2S, the use of MUSIC is not straightforward. Still, following the Prony-type idea [22] we form the S × S data matrix

ℬ = (\begin{matrix} b_{1} & b_{2} & \dots & b_{S} \\ b_{2} & b_{3} & \dots & b_{S + 1} \\ \dots & \dots & \dots & \dots \\ b_{S} & b_{S + 1} & \dots & b_{2 S} \end{matrix})

to extract valuable information from the single vector b. The main observation is that if the illuminations are sufficiently diverse we can find the support of the inhomogeneities with high precision using the MUSIC algorithm; the support is just the zero set of the orthogonal projections of the column vectors of the matrix 𝒜 onto the noise subspace.

Accordingly, we compute the SVD of the data matrix

ℬ = U Σ V^{*} = \sum_{j = 1}^{S} σ_{j} u_{j} v_{j}^{*} .

If the unknown vector x = [δD₁, δD₂, . . ., δD_N] has M nonzero components, with M < S, and the data is noiseless, then there are M nonzero singular values σ₁ > σ₂ > · · · > σ_M > 0. The remaining singular values σ_j, j = M + 1, . . ., S, are essentially zero. The first M singular vectors u_j span the the range or column subspace of ℬ, also referred to as the signal subspace, while the remaining N − M singular vectors u_j span the noise subspace.

Then, the support of δD can be found from the peaks of the imaging functional

ℐ_{n}^{MUSIC} = \frac{{‖ {\tilde{a}}_{n} ‖}_{ℓ_{2}}}{\sum_{j = M + 1}^{N} {| {\tilde{a}}_{n}^{*} u_{j} |}^{2}}, n = 1, \dots, N,

with imaging vectors ã_n ∈ ℂ^S defined as

{\tilde{a}}_{n} = {[G_{z z}^{2} (z_{0}, z_{n}, ω_{1}), G_{z z}^{2} (z_{0}, z_{n}, ω_{2}), \dots, G_{z z}^{2} (z_{0}, z_{n}, ω_{S})]}^{t} .

Once the support is recovered, the problem becomes overdetermined, and the nonzero values of δD can be easily found by solving the corresponding linear system (30) restricted to the given support with a least squares method.

4. Limitations and sources of error

We make here a few remarks concerning the limitations of the proposed method.

The first limitation concerns the range and is due to the exponential decay of the Green’s function with depth. For an infinite SNR, this exponential decay is not an issue, but for a finite noise level we can only image at depths for which the reflectance measurements remain above the noise. This typically restricts the depth to be of the order of the scattering mean free path.

The second remark is about the resolution in cross-range, or in other words, the minimal spacing between consecutive measurements in cross-range. The proposed method uses null source-detector separation measurements and a one-dimensional approximation, which mainly assumes that what is recorded at a particular source-detector cross-range is only due to perturbations at that cross-range. The one-dimensional approximation is no longer valid if we consider small measurement separations in cross-range. The typical cross-range distance for which the modeling error using this assumption can be neglected is of the order of the scattering mean free path or larger.

A finite SNR also places additional limitations on the value of the perturbation that can be detected. In particular, the perturbation must be large enough to be detected over the noise in measurements. Otherwise, the MUSIC algorithm will not be able to distinguish it from noise. In the examples used here, we have considered perturbations between δD = 0.1 cm and 0.8 cm to ensure that they are detectable. Otherwise, the imaging method would not be able to recover these targets at all.

Finally, using the proposed method we can reconstruct either scattering or absorption perturbations, but not both simultaneously. This is so, because scattering and absorption coefficients typically differ by two orders of magnitude. Thus, absorption perturbations can be successfully reconstructed only when the scattering perturbations are negligible and the SNR is large enough.

5. Results

We consider here the same distribution of targets as shown in the top row of Fig. 1. We use N_s = 21 sources and detectors that are 0.3 cm apart, and n_f = 16 frequencies with central frequency 10 GHz and sampling rate 0.13 GHz corresponding to a 2 GHz bandwidth. The SNR is between 30dB and 35dB.

To recover in a first step the supports of the targets, we use the MUSIC algorithm described in Section 3. The results are shown in the top row of Fig. 3. The supports are not recovered exactly. In fact, we do not recover the two nearby targets in the right image of this figure. The reconstruction of the actual values of δD are obtained by solving the linear system (30) restricted to the support determined by MUSIC. The results of this second step are shown in the bottom row of Fig. 3. Both the locations and the quantitative values of δD are quite accurate, and much better than the ones shown in Fig. 1 and Fig. 2 obtained with a minimum ℓ₂-norm approach.

Fig. 3 Top row: Results from the first stage of the imaging algorithm showing MUSIC reconstructions of the support of δD for the same cases shown in Fig. 1. Bottom row: Results from the second stage of the imaging algorithm in which δD quantities are recovored over the support results from the first stage by a least-squares approximation. Here, we have used a 2 GHz system bandwidth.

Download Full Size | PDF

The resolution in depth is determined by the bandwidth B used to form the images, and it is of the order of c/B (in [21], numerical results suggest a 0.1 c/B range resolution for MUSIC). For a 2 GHz bandwidth, this means a poor c/B = 15 cm depth resolution. Hence, in the previous numerical experiments, we could not resolve two obstacles at the same cross-range. Certainly, this limitation can be overcome with a larger bandwidth as it is illustrated in Fig. 4, where we use a 2 THz bandwidth with sampling rate 0.13 THz. More precisely, in Fig. 4 we present the reconstructions corresponding to the target distribution shown in the right image of Fig. 1 that has two nearby targets located at the same cross-range. In the left plot of Fig. 4, we show the minimum ℓ₂-norm solution using all pairs of sources and detectors. In the middle plot, we show the minimum ℓ₂-norm solution using null source-detector separation measurements. Although there is a small improvement with respect to Figs. 1 and 2, where a smaller 2 GHz bandwidth was used, these reconstructions are very similar to those ones. Finally, the right plot of Fig. 4 shows the image obtained with the proposed, two-stage method using a 2 THz bandwidth. As expected, we observe that the image improves by increasing the bandwidth and the two nearby perturbations are now well resolved by the proposed method. Moreover, the reconstructed δD is almost exact, with an ℓ_∞-norm error equal to 2.710⁻¹² and in ℓ₂-norm error equal to 3.310⁻¹².

Fig. 4 Reconstructions of δD using a broader bandwidth of 2 THz with sampling rate 0.13 THz. The left plot shows the minimum ℓ₂-norm solution using all source-detector pairs. The middle plot shows the minimum ℓ₂-norm solution using null source-detector separation data. The right plot shows the MUSIC + least-squares reconstruction.

Download Full Size | PDF

6. Conclusions

We have presented a method for imaging small obstacles near the surface of a diffusive medium. The proposed method is effective and simple as it only uses backscattered light gathered at a few points on the surface. It relies on null source-detector separation measurements, and assumes that the obstacles are well separated in cross-range. This implies that one-dimensional reconstructions can be carried out to accurately resolve the inhomogeneities in depth if large enough bandwidths are available. Under the same conditions we have shown that minimum ℓ₂-norm reconstructions completely fail. The proposed method is limited to penetration depths of the order of one scattering mean-free path.

Funding

Air Force Office of Scientific Research (Grant: FA9550-17-1-0238), MINECO grant FIS2016-77892-R.

Acknowledgments

C. Tsogka and A. D. Kim are supported by the Air Force Office of Scientific Research (Grant: FA9550-17-1-0238). P. Gonzalez-Rodriguez and M. Moscoso are supported by the MINECO grant FIS2016-77892-R.

References

1. J. M. Schmitt, A. H. Gandjbakche, and R. F. Bonner, “Use of polarized light to discriminate short-path photons in a multiply scattering medium,” Appl. Opt. 31, 6535–6546 (1992). [CrossRef] [PubMed]

2. S. G. Demos and R. R. Alfano, “Optical polarization imaging,” Appl. Opt. 36, 150–155 (1997). [CrossRef] [PubMed]

3. M. Moscoso, J. B. Keller, and G. Papanicolaou, “Depolarization and blurring of optical images by biological tissue,” J. Opt. Soc. America A 18, 948–960 (2001). [CrossRef]

4. K. M. Yoo, B. B. Das, and R. R. Alfano, “Imaging of translucent object hidden in highly scattering medium from the early portion of the diffuse component of a transmitted ultrafast laser pulse,” Opt. Lett. 17, 958–960 (1992). [CrossRef] [PubMed]

5. B. B. Das, K. M. Yoo, and R. R. Alfano, “Ultrafast time-gated imaging in thick tissues: a step toward optical mammography,” Opt. Lett. 18, 1092–1094 (1993). [CrossRef] [PubMed]

6. D. Contini, A. D. Mora, L. Spinelli, A. Farina, A. Torricelli, R. Cubeddu, F. Martelli, G. Zaccanti, A. Tosi, G. Boso, F. Zappa, and A. Pifferi, “Effect of time-gated detection in diffuse optical imaging at short source-detector separation,” J. Phys. D: Appl. Phys. 48, 045401 (2015). [CrossRef]

7. S. G. Demos, H. B. Radousky, and R. R. Alfano, “Deep subsurface imaging in tissues using spectral and polarization filtering,” Opt. Express 7, 23–28 (2000). [CrossRef] [PubMed]

8. P. González-Rodríguez, A. D. Kim, and M. Moscoso, “Robust depth selectivity in mesoscopic scattering regimes using angle-resolved measurements,” Opt. Lett. 38, 787–789 (2013). [CrossRef] [PubMed]

9. D. Huang, E. A. Swanson, C. P. Lin, J. S. Schuman, W. G. Stinson, W. Chang, and J. G. Fujimoto, “Optical coherence tomography,” Science 254, 1178–1181 (1991). [CrossRef] [PubMed]

10. R. O. Schmidt, “Multiple emitter location and signal parameter estimation,” IEEE Trans. Antennas Propag. 34, 276–280 (1986). [CrossRef]

11. M. H. Hayes, Statistical Digital Signal Processing and Modeling (John Wiley & Sons, 1996).

12. A. Devaney, E. Marengo, and F. Gruber, “Time-reversal-based imaging and inverse scattering of multiply scattering point targets,” J. Acoust. Soc. Am. 118, 3129–3138 (2005). [CrossRef]

13. R. Griesmaier and C. Schmiedecke, “A multifrequency MUSIC algorithm for locating small inhomogeneities in inverse scattering,” Inverse Probl. 33, 035015 (2017). [CrossRef]

14. W. Liao and A. Fannjiang, “A MUSIC for single-snapshot spectral estimation: stability and super-resolution,” Appl. Comput. Harmon. Anal. 40, 33–67 (2016). [CrossRef]

15. B. V. P. Dileep, T. Das, and P. K. Dutta, “Modified CS-MUSIC for diffuse optical tomography using joint sparsity,” Optik 158, 1478–1490 (2018). [CrossRef]

16. S. B. Rohde and A. D. Kim, “Backscattering of continuous and pulsed beams,” Multiscale Model. Simul. 15, 1356–1375 (2017). [CrossRef]

17. L. V. Wang and H.-I. Wu, Biomedical Optics: Principles and Imaging (John Wiley & Sons, 2012).

18. A. Ishimaru, Wave Propagation and Scattering in Random Media (Wiley-IEEE, 1997).

19. L. Spinelli, F. Martelli, S. Del Bianco, A. Pifferi, A. Torricelli, R. Cubeddu, and G. Zaccanti, “Absorption and scattering perturbations in homogeneous and layered diffusive media probed by time-resolved reflectance at null source-detector separation,” Phys. Rev. E 74, 021919 (2006). [CrossRef]

20. A. Pifferi, A. Torricelli, L. Spinelli, D. Contini, R. Cubeddu, F. Martelli, G. Zaccanti, A. Tosi, A. D. Mora, F. Zappa, and S. Cova, “Time-resolved diffuse reflectance using small source-detector separation and fast single-photon gating,” Phys. Rev. Lett. 100, 138101 (2008). [CrossRef] [PubMed]

21. M. Moscoso, A. Novikov, G. Papanicolaou, and C. Tsogka, “Robust multifrequency imaging with MUSIC,” submitted, 2018.

22. R. Prony, “Essai experimental et analytique sur les lois de la dilatabilite de fluides elastiques et sur celles del la force expansive de la vapeur de l’alkool a differentes temperatures,” J. L’Ecole Polytech.1, 24–76 (1795).

Quantitative subsurface imaging in strongly scattering media

Abstract

1. Introduction

2. Modeling measurements

2.1. Radiative transfer

2.2. Asymptotic solution in strongly scattering media

2.3. Born approximation

2.4. Linear system

2.5. Initial imaging results

3. Image reconstruction method

4. Limitations and sources of error

5. Results

6. Conclusions

Funding

Acknowledgments

References

Cited By

Figures (4)

Equations (35)

Optics Express