High-Resolution Interferometric Synthetic Aperture Imaging in scattering media

TL;DR

This paper presents a high-resolution interferometric synthetic aperture imaging method that effectively mitigates wave distortion in scattering media, enabling sharp imaging through a phase retrieval process.

Contribution

It introduces a novel high-resolution interferometric approach based on CINT that preserves image sharpness despite scattering distortions.

Findings

Mitigates wave distortion using empirical cross-correlations.

Achieves robust estimation of the Fourier modulus of reflectivity.

Enables high-resolution imaging via phase retrieval.

Abstract

The goal of synthetic aperture imaging is to estimate the reflectivity of a remote region of interest by processing data gathered with a moving sensor which emits periodically a signal and records the backscattered wave. We introduce and analyze a high-resolution interferometric method for synthetic aperture imaging through an unknown scattering medium which distorts the wave. The method builds on the coherent interferometric (CINT) approach which uses empirical cross-correlations of the measurements to mitigate the distortion, at the expense of a loss of resolution of the image. The new method shows that, while mitigating the wave distortion, it is possible to obtain a robust and sharp estimate of the modulus of the Fourier transform of the reflectivity function. A high-resolution image can then be obtained by a phase retrieval algorithm.

Equations

$$\text{mbox}data=\{R_n(t)\text{mbox}fort\in (0,T),n=0,\dots ,N\}.$$

$$\displaystyle\frac{1}{c^{2}}\Big{[}1+\sigma\mu\Big{(}\frac{{\itbf x}}{\ell_{c}}\Big{)}+\rho({\itbf x})\Big{]}\partial_{t}^{2}u_{n}(t,{\itbf x})-\Delta u_{n}(t,{\itbf x})$$

$$u_n(t,\text{itbf}x)$$

$${\cal R}({\itbf x}-{\itbf x}^{\prime})=\mathbb{E}\big{[}\mu({\itbf x})\mu({\itbf x}^{\prime})\big{]},$$

$$\mathcal{R}(0)=1\text{mbox}and{\int }_{{\mathbb{R}}^2}d\text{itbf}x\mathcal{R}(\text{itbf}x)=1.$$

$$\mathcal{R}(\text{itbf}x)=\exp (-\pi |\text{itbf}x{|}^2).$$

$$R_n(t)$$

$$=\frac{1}{2\pi }{\int }_{\mathbb{R}}d\omega e^{-i\omega (t-nT)}\hat{s}(\omega )k^2(\omega ){\int }_{{\mathbb{R}}^2}d\text{itbf}y\rho (\text{itbf}y){\hat{G}}_{\mu }^2(\omega ,\text{itbf}y,{\text{itbf}x}_n)+W_n(t),$$

$$\hat{s}(\omega )={\int }_{\mathbb{R}}dts(t)e^{i\omega t},s(t)=\frac{1}{2\pi }{\int }_{\mathbb{R}}d\omega e^{-i\omega t}\hat{s}(\omega ).$$

$$\mathbb{E}\big{[}\overline{\widehat{W}_{n}(\omega)}\widehat{W}_{n^{\prime}}(\omega^{\prime})\big{]}=\sigma_{\rm W}^{2}\delta(\omega-\omega^{\prime})\delta_{nn^{\prime}},$$

$$\displaystyle\widehat{G}(\omega,{\itbf x},{\itbf y})=\frac{i}{4}H_{0}^{(1)}\big{(}k(\omega)|{\itbf x}-{\itbf y}|\big{)}\approx\frac{\exp\big{(}ik(\omega)|{\itbf x}-{\itbf y}|+i\frac{\pi}{4}\big{)}}{2^{3/2}\sqrt{\pi k(\omega)|{\itbf x}-{\itbf y}|}},$$

$${\mathcal{I}}_{\mathrm{SAR}}({\text{itbf}y}^S)$$

$$\displaystyle=\Big{|}\frac{1}{2\pi}\sum_{n=0}^{N}\int_{\mathbb{R}}d\omega\,\widehat{R}_{n}(\omega)\overline{\widehat{F}_{n}(\omega,{\itbf y}^{S})}e^{-i\omega nT}\Big{|}^{2},$$

$$F_n(t,{\text{itbf}y}^S)=\frac{1}{2\pi }{\int }_{\mathbb{R}}d\omega e^{-i\omega t}{\hat{F}}_n(\omega ,{\text{itbf}y}^S),{\hat{F}}_n(\omega ,{\text{itbf}y}^S)=\hat{s}(\omega ){\hat{G}}^2(\omega ,{\text{itbf}y}^S,{\text{itbf}x}_n).$$

$${\mathcal{I}}_{\mathrm{CINT}}({\text{itbf}y}^S)=$$

$$\displaystyle\times e^{i\omega nT-i\omega^{\prime}n^{\prime}T}\exp\Big{(}-\frac{|{\itbf x}_{n}-{\itbf x}_{n^{\prime}}|^{2}}{2X^{2}}-\frac{(\omega-\omega^{\prime})^{2}}{2\Omega^{2}}\Big{)}.$$

$$\mathcal{I}({\text{itbf}y}^S,{{\text{itbf}y}^S}^{\prime })=$$

$$\displaystyle\times e^{i\omega nT-i\omega^{\prime}n^{\prime}T}\exp\Big{(}-\frac{|{\itbf x}_{n}-{\itbf x}_{n^{\prime}}|^{2}}{2X^{2}}-\frac{(\omega-\omega^{\prime})^{2}}{2\Omega^{2}}\Big{)}.$$

$${\mathcal{I}}_{\mathrm{CINT}}({\text{itbf}y}^S)=\mathcal{I}({\text{itbf}y}^S,{\text{itbf}y}^S).$$

$${\cal I}_{\rm HCINT}(\widetilde{\itbf y}^{S})=\int_{\mathbb{R}^{2}}d{{\itbf y}^{S}}\,{\cal I}\Big{(}{{\itbf y}^{S}}+\frac{\widetilde{\itbf y}^{S}}{2},{{\itbf y}^{S}}-\frac{\widetilde{\itbf y}^{S}}{2}\Big{)}.$$

$${\hat{\mathcal{I}}}_{\mathrm{HCINT}}(\mathbf{\kappa })$$

$$\displaystyle=\int_{\mathbb{R}^{2}}d{\itbf y}^{S}\int_{\mathbb{R}^{2}}d{\widetilde{\itbf y}^{S}}\,{\cal I}\Big{(}{{\itbf y}^{S}}+\frac{\widetilde{\itbf y}^{S}}{2},{{\itbf y}^{S}}-\frac{\widetilde{\itbf y}^{S}}{2}\Big{)}e^{-i\bm{\kappa}\cdot\widetilde{\itbf y}^{S}}.$$

$${\sigma }^2\frac{L^3}{{\ell }_c^3}\ll \frac{{\lambda }_o^2}{{\sigma }^2{\ell }_{c}L}\ll 1.$$

$${\hat{G}}_{\mu }(\omega ,\text{itbf}x,\text{itbf}y)$$

$${\cal T}_{\mu}({{\itbf x}},{{\itbf y}})=\frac{\sigma|{{\itbf x}}-{{\itbf y}}|}{2c}\int_{0}^{1}dh\,\mu\Big{(}\frac{{{\itbf y}}+h({{\itbf x}}-{{\itbf y}})}{\ell_{c}}\Big{)}$$

$$\mathbb{E}\big{[}{\cal T}_{\mu}({\itbf x}_{n},{\itbf y}){\cal T}_{\mu}({\itbf x}_{n^{\prime}},{\itbf y}^{\prime})\big{]}=\tau^{2}{\cal C}\Big{(}\frac{|{\itbf x}_{n}-{\itbf x}_{n^{\prime}}|}{\ell_{c}}\Big{)},\qquad{\cal C}(r)=\frac{1}{r}\int_{0}^{r}dh\,e^{-\pi h^{2}}.$$

$$\tau =\frac{\sigma \sqrt{{\ell }_{c}L}}{2c},$$

$$\omega \tau \sim {\omega }_o\tau \gg 1,$$

$$\displaystyle\mathbb{E}\big{[}\exp\big{(}2i\omega_{1}{\cal T}_{\mu}({\itbf x}_{n_{1}},{\itbf y}_{1})\big{)}\big{]}=$$

$$\displaystyle\mathbb{E}\big{[}\exp\big{(}2i\omega_{1}{\cal T}_{\mu}({\itbf x}_{n_{1}},{\itbf y}_{1})-2i\omega_{2}{\cal T}_{\mu}({\itbf x}_{n_{2}},{\itbf y}_{2})\big{)}\big{]}=$$

$$\displaystyle-4\omega_{1}\omega_{2}\tau^{2}\Big{[}1-{\cal C}\Big{(}\frac{|{\itbf x}_{n_{1}}-{\itbf x}_{n_{2}}|}{\ell_{c}}\Big{)}\Big{]}\Big{\}},$$

$$\displaystyle\mathbb{E}\big{[}\exp\big{(}2i\omega_{1}{\cal T}_{\mu}({\itbf x}_{n_{1}},{\itbf y}_{1})-2i\omega_{2}{\cal T}_{\mu}({\itbf x}_{n_{2}},{\itbf y}_{2})-2i\omega_{3}{\cal T}_{\mu}({\itbf x}_{n_{3}},{\itbf y}_{3})+2i\omega_{4}{\cal T}_{\mu}({\itbf x}_{n_{4}},{\itbf y}_{4})\big{)}\big{]}$$

$$\displaystyle=\exp\bigg{\{}-2\tau^{2}\Big{[}\sum_{j=1}^{4}\omega_{j}^{2}+2\omega_{1}\omega_{4}{\cal C}\Big{(}\frac{|{\itbf x}_{n_{1}}-{\itbf x}_{n_{4}}|}{\ell_{c}}\Big{)}+2\omega_{2}\omega_{3}{\cal C}\Big{(}\frac{|{\itbf x}_{n_{2}}-{\itbf x}_{n_{3}}|}{\ell_{c}}\Big{)}$$

$$\displaystyle\quad-2\omega_{1}\omega_{2}{\cal C}\Big{(}\frac{|{\itbf x}_{n_{1}}-{\itbf x}_{n_{2}}|}{\ell_{c}}\Big{)}-2\omega_{1}\omega_{3}{\cal C}\Big{(}\frac{|{\itbf x}_{n_{1}}-{\itbf x}_{n_{3}}|}{\ell_{c}}\Big{)}-2\omega_{2}\omega_{4}{\cal C}\Big{(}\frac{|{\itbf x}_{n_{2}}-{\itbf x}_{n_{4}}|}{\ell_{c}}\Big{)}$$

$$\displaystyle\quad-2\omega_{3}\omega_{4}{\cal C}\Big{(}\frac{|{\itbf x}_{n_{3}}-{\itbf x}_{n_{4}}|}{\ell_{c}}\Big{)}\Big{]}\bigg{\}}.$$

Full text

1 Introduction

In synthetic aperture imaging, a sensor mounted on a moving platform (e.g., plane, satellite) emits a signal $s(t-nT)$ at time instants $nT$, counted by $n=0,1,\ldots,N$, and records the backscattered wave, the “response” $R_{n}(t)$. The goal of imaging is to estimate the reflectivity $\rho({\itbf x})$ of a remote region of interest from

[TABLE]

The signal $s(t)$ is either a broad-band pulse defined by an envelope of small temporal support of order $1/B\ll T$, where $B$ is the bandwidth, modulated at carrier frequency $\omega_{o}\gg B$, or a chirp that can be compressed to a pulse with data processing [5]. We assume the former. The variables $t$ and $nT$ are referred to as the “fast time” and the “slow time”, respectively. The trajectory of the platform can be arbitrary, but for simplicity, and without loss of generality, we suppose that it is straight and the motion is uniform, so that the signal emission is from the regularly spaced positions ${\itbf x}_{n}$, for $n=0,\ldots,N$ (see Fig. 1). The line segment connecting ${\itbf x}_{0}$ to ${\itbf x}_{N}$ is called the synthetic aperture and its length $a=|{\itbf x}_{N}-{\itbf x}_{0}|$ is the aperture size.

The data (1) have only two degrees of freedom, so it is not possible to estimate a reflectivity function $\rho({\itbf x})$ in three dimensions. Here we consider the problem in two dimensions, but the results extend to imaging in three dimensions on a surface with known topography. Again, for simplicity, we center the aperture above the remote region ${\cal D}$ of interest, and we introduce the system of coordinates ${\itbf x}=(x_{\parallel},x_{\perp})$ with origin at the center of ${\cal D}$, the “range coordinate” $x_{\parallel}$ measured along the main direction of propagation, orthogonal to the aperture, and “cross-range coordinate” $x_{\perp}$ measured along the aperture.

The classic synthetic aperture imaging method assumes that the medium between the sensor and the imaging region is non-scattering, with smooth and known wave speed. Without going into technical details, the imaging function is given roughly by the sum over $n$ of the returns $R_{n}(t)$ evaluated (synchronized) at the roundtrip travel time between ${\itbf x}_{n}$ and the imaging point [4, 5]. When this point lies in the support of $\rho({\itbf x})$, denoted by $\mbox{supp}(\rho)$, the synchronized returns add constructively and the imaging function is large. Therefore, the set $\mbox{supp}(\rho)$ can be estimated from the imaging function displayed in ${\cal D}$ above some user defined threshold value. In the particular case of a homogeneous medium with constant wave speed $c$, and for an idealized reflectivity supported at two points, well known resolution formulas state that these points can be distinguished if they are separated by a distance of order $c/B$ in range and $\lambda_{o}L/a$ in cross-range, where $\lambda_{o}=2\pi c/\omega_{o}$ is the carrier wavelength and $L$ is the range offset from the aperture.

We are interested in imaging in heterogeneous media with microstructure, as sketched in Fig. 1, where the wave speed fluctuates about a known reference profile, which we take equal to the constant $c$ for simplicity. The fluctuations have small amplitude and occur on a length scale that is much smaller than $L$. They are unknown and cannot be estimated as part of imaging from the band limited data (1). Thus, there is uncertainty in the wave propagation, which motivates modeling the wave speed as a random perturbation of the reference $c$. This model introduces a stochastic framework where we can quantify the robustness of imaging methods with respect to the uncertainty of the microstructure. Robust images cannot be obtained by empirical averaging over many realizations of the random wave speed, because the imaging experiment occurs in a single medium. However, with careful data processing, it is possible to get images that are practically insensitive to the particular realization i.e., are statistically stable.

The cumulative scattering effect of the microstructure (the wave distortion) depends in a complicated way on the amplitude of the fluctuations of the wave speed as well as the relation between the length scale of the fluctuations, the carrier wavelength and the travel distance. When the distortion is strong, the classic synthetic aperture imaging method [4, 5] gives noisy images that are difficult to interpret and unreliable (statistically unstable). Different data processing is needed to mitigate the wave distortion, based on the empirical cross-correlations of the measurements. The coherent interferometric (CINT) method [1, 2] forms an image using such cross-correlations. It is known [13, 14, 16, 17] that scattering causes statistical decorrelation of the time-harmonic components of the wave field over a frequency offset $\Omega_{d}$ called “decoherence frequency” and a spatial offset $X_{d}$ called “decoherence length”. These scales depend on the statistics of the fluctuations of the wave speed, not the particular realization, and CINT takes them into account by calculating the empirical cross-correlations in a time window of duration $1/\Omega$ and for sensor locations that are within a distance $X$ of each other. The image is then formed by superposing the empirical cross-correlations synchronized relative to the imaging point with travel time delays calculated in the reference medium. There is a trade-off between the resolution of the image and its robustness to the uncertainty of the microstructure, which is quantified by the threshold parameters $X$ and $\Omega$ [1, 2]. The smaller these are, the less sensitive is the CINT image to the microstructure, as long as the aperture and the bandwidth are large enough. However, the range resolution is of the order $c/\Omega$ and the cross-range resolution is of the order of $\lambda_{o}L/X$, so robustness comes at the cost of loss of resolution. The optimal choice is $X\approx X_{d}<a$ and $\Omega\approx\Omega_{d}<B$, and in practice this can be determined by optimizing a measure of quality of the image [2].

CINT has been used for imaging with arrays of sensors [1, 2] (see also references therein) and with synthetic apertures [11]. A modification of CINT introduced recently in [3] in the context of imaging a constellation of point sources with a passive array of receivers shows that it is possible to localize the sources with resolution that is comparable to that in the homogeneous medium. More precisely, nearby point sources within a blurry peak of the CINT function can be resolved with resolution $c/B$ in range and $\lambda_{o}L/a$ in cross-range, up to an overall (rigid body) translation and rotation of the constellation. The algorithm in [3] involves a point search (it is targeted toward imaging a constellation of points) and has prohibitive computational cost for many sources. In this paper we extend the ideas in [3] to synthetic aperture imaging of a general reflectivity $\rho({\itbf x})$, based on a new HCINT imaging function, where the acronym stands for high-resolution CINT. We show that this function is on one hand robust to the uncertainty of the fluctuations of the wave speed and on the other hand it allows a precise estimate of the modulus of the Fourier transform of $\rho({\itbf x})$. A high-resolution image can then be obtained from this estimate using phase retrieval [6, 7, 8, 9, 15].

Our mathematical analysis of HCINT is based on a geometrical optics model of wave propagation through random media. This simple model accounts for wavefront distortion and allows an explicit quantification of robustness i.e., calculation of the variance of the imaging function. The HCINT method is not model specific, and it can also be analyzed with more complex wave propagation models like in [12].

The paper is organized as follows: We begin in section 2 with the mathematical formulation of the problem and the expression of the three imaging functions: classic synthetic aperture imaging, CINT and HCINT. The model of wave propagation in the random medium is described in section 3 and it is used to analyze the imaging functions in section 4. The estimation of the modulus of the Fourier coefficients of the reflectivity function and the subsequent imaging based on phase retrieval is in section 5. We present numerical results in section 6 and end with a summary in section 7.

2 Formulation of the problem and the imaging functions

The wave $u_{n}(t,{\itbf x})$ emitted from the location ${\itbf x}_{n}=(L,x_{n\perp})$ in the aperture satisfies

[TABLE]

for $n=0,\ldots,N$, time $t\in\mathbb{R}$ and position ${\itbf x}\in\mathbb{R}^{2}$, with the initial condition

[TABLE]

Here $c$ is the constant reference wave speed, $\mu$ models the random fluctuations and $\rho$ is the unknown reflectivity, assumed compactly supported away from the aperture. The random process $\mu$ is statistically homogeneous, with mean zero and integrable autocovariance

[TABLE]

normalized so that

[TABLE]

We assume henceforth, for convenience and without loss of generality, the Gaussian autocovariance

[TABLE]

The dimensionless parameter $\sigma$ in (2) is the standard deviation of the random fluctuations and the length scale $\ell_{c}$ is the correlation length.

The inverse problem is to estimate the reflectivity $\rho$ from the data (1). We study its solution using the three imaging functions given in sections 2.1–2.3. Their expression is based on two standard approximations: (1) the “start-stop” approximation [4, 5] which assumes that the sensor movement during the roundtrip travel time to the imaging region ${\cal D}$ is negligible; (2) the single scattering (Born) approximation which assumes that the reflectivity $\rho$ is not too strong. The data model is

[TABLE]

where $k(\omega)=\omega/c$ is the wavenumber and we denote with “hat” the Fourier transform with respect to time, defined with the convention

[TABLE]

The propagation through the random medium is modeled in (7) with the Green’s function $\widehat{G}_{\mu}$ of the Helmholtz equation with wave speed c\big{[}1+\sigma\mu({\itbf x}/\ell_{c})\big{]}^{-1/2} and $W_{n}$ denotes additive noise. For convenience, we let $W_{n}(t)$ be Gaussian, white in time $t$ and uncorrelated in $n$, with mean zero and covariance given in the frequency domain by

[TABLE]

where $\delta_{nn^{\prime}}$ stands for the Kronecker symbol. The bar is used throughout to denote complex conjugate.

2.1 Synthetic aperture radar (SAR) imaging

The fluctuations of the wave speed are neglected in the standard synthetic aperture imaging method, meaning that the wave propagation is modeled by the Green’s function in the reference medium

[TABLE]

where $H_{0}^{(1)}$ is the Hankel function of the first kind and of order [math] and the approximation is for a large distance $|{\itbf x}-{\itbf y}|$ with respect to the wavelength $\lambda=2\pi/k(\omega)$.

The imaging function

[TABLE]

is the superposition of the data (1) convolved (matched filtered) with

[TABLE]

Note that $F_{n}$ is the data model for a point reflector at the search (imaging) point ${\itbf y}^{S}$ in the reference medium. The matched filtering is called henceforth “backpropagation” to the imaging point ${\itbf y}^{S}$. At long range $|{\itbf x}_{n}-{\itbf y}^{S}|\sim L\gg a$, where “$\sim$” means of the order of, and neglecting constant amplitude factors, the backpropagation amounts to evaluating the response $R_{n}(t)$ at the roundtrip travel time $2|{\itbf x}_{n}-{\itbf y}^{S}|/c$. The index SAR in (10) is the acronym for synthetic aperture radar, the most common application of the synthetic aperture imaging modality.

2.2 The CINT imaging function

The CINT imaging function is obtained by backpropagating selected empirical cross-correlations of the measured responses, for nearby pairs of sensor locations ${\itbf x}_{n},{\itbf x}_{n^{\prime}}$ and frequencies $\omega,\omega^{\prime}$,

[TABLE]

We use Gaussian sensor offset and frequency windows for convenience in the calculations, with standard deviations $X$ and $\Omega$ accounting for the decorrelation of the wave components due to scattering, as explained in the analysis in section 4.

2.3 The HCINT imaging function

Instead of (12), consider the CINT-like imaging function defined for two nearby imaging points ${\itbf y}^{S}$ and ${{\itbf y}^{S}}^{\prime}$, as proposed in a passive array imaging context in [3],

[TABLE]

We call it the “two-point CINT” imaging function and note that it is a generalization of the CINT imaging function (12) since we have

[TABLE]

The HCINT function is defined by the integral of (13) over the center locations

[TABLE]

We will also use its Fourier transform

[TABLE]

We will show in sections 4–5 that ${\cal I}_{\rm CINT}$ gives a statistically stable but low resolution image of the reflectivity $\rho$, whereas the HCINT imaging function (16) gives an estimate of the modulus of the Fourier transform of $\rho$. This estimate can then be used in a phase retrieval algorithm [6, 7, 8, 9, 15] to get a high-resolution image of the reflectivity.

In practice, HCINT may be used to improve the resolution by zooming a region of interest in the support of the CINT image. The integrals over ${{\itbf y}^{S}}$ and $\widetilde{\itbf y}^{S}$ in (16) can be evaluated with numerical quadrature formulas in such a small zoom region.

3 Random travel time model

In this section we briefly review the geometrical optics model of wave propagation through the random medium with wave speed $c[1+\sigma\mu({\itbf x}/\ell_{c})]^{-1/2}$. Its derivation is given in [10, Section 12.1] under the high-frequency scaling assumption $\lambda_{o}\ll\ell_{c}<L$ and the weak fluctuations assumption $\sigma^{2}\ll(\ell_{c}/L)^{3}$, so that scattering does not effect the amplitude of the wave and the rays remain straight. It is only the travel time calculated along the straight ray that is randomized. We are interested in a long range $L\gg\ell_{c}$, where the random travel time fluctuations have Gaussian statistics (even if $\mu$ is not Gaussian), by the central limit theorem. To showcase the effect of the random medium, we assume that these fluctuations are large (wave front is strongly distorted), which amounts to having

[TABLE]

The Green’s function is

[TABLE]

where $\widehat{G}$ is given in (9) and

[TABLE]

models the random fluctuations of the travel time, given by the line integral of the random process $\mu$ along the straight ray connecting ${\itbf x}$ and ${\itbf y}$. For points ${\itbf y},{\itbf y}^{\prime}$ in the neighborhood of the origin, satisfying $|{\itbf y}-{\itbf y}^{\prime}|<\ell_{c}$, and for ${\itbf x}_{n},{\itbf x}_{n^{\prime}}$ in the aperture, the process ${\cal T}_{\mu}$ has Gaussian statistics with mean zero and covariance function

[TABLE]

Here we used the assumption (6) and introduced the time scale

[TABLE]

which quantifies the standard deviation of the random fluctuations of the travel time. Note that

[TABLE]

by the assumption (17), so the phase of the Green’s function (18) has very large fluctuations.

For arbitrary four points $({\itbf x}_{n_{j}})_{j=1,\ldots,4}$ in the aperture, indexed by $0\leq n_{j}\leq N$, four points $({\itbf y}_{j})_{j=1,\ldots,4}$ in the search (imaging) region ${\cal D}$ with diameter smaller than $\ell_{c}$, and for frequencies $(\omega_{j})_{j=1,\ldots,4}$, we have by the Gaussian property of ${\cal T}_{\mu}$ that

[TABLE]

and

[TABLE]

We conclude from equations (18), (22) and (23) that

[TABLE]

Physically, this means that the wave front is strongly distorted (randomized) due to scattering, so that averaging the Green’s function over realizations of the wave speed gives a negligible result. We also get from (24) and the assumption on the bandwidth

[TABLE]

that the second moments (25) are negligible unless ${\itbf x}_{n_{1}}$ and ${\itbf x}_{n_{2}}$ are nearby. For such points we can expand the covariance in the phase of (24) around the origin and obtain the simpler formula

[TABLE]

The scales of decay in the sensor and frequency offsets

[TABLE]

quantify the statistical decorrelation of the wave in the random medium and are called the decoherence length and decoherence frequency.

4 Analysis of the imaging functions

We now use the random travel time model of wave propagation to analyze the three imaging functions given in sections 2.1–2.3. We compare the results to those in the reference homogeneous medium and account for the effects of the additive noise, as well.

In the analysis we choose a probing pulse with Gaussian envelope

[TABLE]

and with Fourier transform

[TABLE]

The successive positions of the sensor are close to each other

[TABLE]

so we approximate the sums over the index $n$ of ${\itbf x}_{n}$ by integrals over the aperture and take the Gaussian appodization $\exp(-x_{\perp}^{2}/a^{2})$

[TABLE]

The aperture size $a$ is assumed smaller than the range $L$, so we can use the paraxial approximation

[TABLE]

for all points ${\itbf x}=(L,x_{\perp})$ in the aperture and points ${\itbf y}^{S}=(y_{\parallel}^{S},y_{\bot}^{S})$ in the imaging region ${\cal D}$. These choices lead to explicit expressions of the imaging functions, but do not play an essential role in the conclusions.

4.1 Analysis of the SAR imaging function

The index $n$ is no longer needed in the continuum aperture approximation (31), so we change slightly the notation

[TABLE]

From now on, $W(t,x_{\perp})$ is a Gaussian process, white in time $t$ and in space $x_{\perp}$, with mean zero and covariance given in the frequency domain by

[TABLE]

The expression (10) of the imaging function becomes

[TABLE]

with

[TABLE]

and

[TABLE]

4.1.1 Homogeneous medium

In the absence of the fluctuations of the wave speed and noise, the expression of the imaging function (34) would be

[TABLE]

where $k_{o}=k(\omega_{o})=2\pi/\lambda_{o}$ and

[TABLE]

This is obtained with straightforward calculations from equations (9), (34–36) and the paraxial approximation (32).

In the particular case of a point reflector at location ${\itbf y}_{\star}$, the imaging function is proportional to the square modulus of the “point-spread function” ${\cal K}_{a,B}$,

[TABLE]

This and the expression (38) show that the resolution in the cross-range direction is $\sim\lambda_{o}L/a$ and in the range direction is $\sim c/B$, as reported in the literature [4].

4.1.2 Random medium

To explain the behavior of the SAR imaging function in the random medium, we describe here its expectation and covariance in the noiseless case (i.e., when there is no additive noise). The effect of the additive noise is analyzed in the next section.

We obtain from equations (9), (32–36) and the moment formula (27) that

[TABLE]

with the same constant ${\rm C}$ as in (37) and with $\widetilde{a}$ and $\widetilde{B}$ defined by

[TABLE]

If the random medium is strong enough so that the decoherence parameters (28) satisfy $X_{d}<a$ and/or $\Omega_{d}<B$, we obtain from (40) that the amplitude of expectation of the imaging function is reduced and there is loss of resolution. This is evident in the case of a single point scatterer, where

[TABLE]

Comparing with (39) and using the definition (38), we note that the peak amplitude of (42) is smaller by the factor $(\widetilde{a}/a)^{2}(\widetilde{B}/B)^{2}$ and the resolution is reduced to ${\lambda_{o}L}/{\widetilde{a}}$ and $c/\widetilde{B}$ in the cross-range and range direction. These resolution limits are the scales of decay of the kernel ${\cal K}_{\widetilde{a},\widetilde{B}}$ defined in (38), with $a$ and $B$ replaced by $\widetilde{a}$ and $\widetilde{B}$.

The calculation of the covariance of ${\cal I}_{\rm SAR}$ involves the fourth order moments (25), for distinct points and frequencies satisfying

[TABLE]

Here we introduced the center frequencies $\omega,\omega^{\prime}\sim\omega_{o}$ and the frequency offsets which satisfy $|\widetilde{\omega}|,|\widetilde{\omega}^{\prime}|\lesssim 2B$. Similarly, ${\itbf x},{\itbf x}^{\prime}$ are center points in the aperture i.e., with range coordinate $L$ and cross-range coordinates in the interval $(-a/2,a/2)$. The spatial offsets are $\widetilde{{\itbf x}}=(0,\widetilde{x}_{\perp})$ and $\widetilde{{\itbf x}}^{\prime}=(0,\widetilde{x}_{\perp}^{\prime})$, with $|\widetilde{x}_{\perp}|,|\widetilde{x}_{\perp}^{\prime}|\lesssim a$.

If the diameter of the search region ${\cal D}$ is smaller than $\ell_{c}$ and if $X_{d}<a$ and $\Omega_{d}<B$, then we find the covariance

[TABLE]

This indicates that the image displays bright and dark spots, so-called “speckle” of size of the order of the correlation radius $\lambda_{o}L/\widetilde{a}$ in the cross-range direction and $c/\widetilde{B}$ in the range direction, the scales of decay of the kernel ${\cal K}_{\widetilde{a},\widetilde{B}}$. Moreover, if we let ${\itbf y}^{S}={{\itbf y}^{S}}^{\prime}$ in this expression, we obtain that the variance is equal to the square mean:

[TABLE]

Therefore, the coefficient of variation (i.e., the ratio of the standard deviation over the mean) of the SAR image near its peak values is large,

[TABLE]

i.e., the SAR image has strong random fluctuations, it is not statistically stable.

4.1.3 Additive noise

The presence of the Gaussian additive noise introduces an additional speckle pattern in the image, modeled by

[TABLE]

This is independent of the fluctuations of the wave speed in the random medium. We now describe the mean and correlation radius of (45), where the latter gives the typical noise induced speckle size.

We obtain after straightforward calculations, using definitions (33), (9), (11) and the paraxial approximation (32) that the mean of the speckle pattern is uniform

[TABLE]

and the covariance is

[TABLE]

where we have used the fourth moment property satisfied by the Gaussian process

[TABLE]

By letting ${\itbf y}^{S}={{\itbf y}^{S}}^{\prime}$ in the last equation we obtain the variance

[TABLE]

The decay of the covariance (47) shows that the noise induced speckle pattern has correlation radius $\lambda_{o}L/a$ in the cross range direction and $c/B$ in the range direction, which means that the image displays bright and dark spots with these typical sizes. This is in addition to the random fluctuations induced by scattering in the random medium.

4.2 Analysis of the two-point CINT imaging function

We now describe the mean and variance of the two-point CINT imaging function

[TABLE]

which is related to CINT by (14). Its use in the high-resolution CINT imaging method is discussed in section 5.

4.2.1 The mean

The expression of the mean is obtained from definitions (13), (9), (11) and equations (32–36). The calculation is the same in the homogeneous and the random medium, except that in the latter case we use the moment formula (27) and in the former case there is no need for the expectation if there is no noise.

We obtain that in the noiseless case

[TABLE]

with the same constant ${\rm C}$ as in (37) and with kernels

[TABLE]

The first kernel gives the resolution in the central spatial variable, modeled by the decay of (52) on the scale $\lambda_{o}L/\widetilde{X}$ in the cross-range direction and $c/\widetilde{\Omega}$ in the range direction. If the medium were homogeneous, $\widetilde{X}$ and $\widetilde{\Omega}$ would depend on the window parameters $X$ and $\Omega$ in the definition (13) of ${\cal I}({\itbf y}^{S},{{\itbf y}^{S}}^{\prime})$, the aperture $a$ and the bandwidth $B$ of the probing pulse as follows

[TABLE]

In the random medium, they also depend on the decoherence length $X_{d}$ and frequency $\Omega_{d}$ defined in (28),

[TABLE]

The second kernel (53) gives the resolution in the spatial offset, modeled by the decay on the scale $\lambda_{o}L/a$ in the cross-range direction and $c/B$ in the range direction. These scales coincide with the resolution limits of the SAR imaging function in the homogeneous medium.

We conclude that the mean of the two-point CINT image displays excellent resolution in the spatial offset locations ${\itbf y}^{S}-{{\itbf y}^{S}}^{\prime}$, and reduced resolution in the mid-point locations $({\itbf y}^{S}+{{\itbf y}^{S}}^{\prime})/2$. The latter is the same as the resolution of the CINT imaging function, obtained from (14) by setting ${\itbf y}^{S}={{\itbf y}^{S}}^{\prime}$ in equation (51),

[TABLE]

In the particular case of a single point scatterer, the CINT imaging function has the simple expression

[TABLE]

and its resolution can be compared easily to that of the mean SAR image in (42). We note in particular that the smaller the window parameter $\Omega$ and $X$ are, the worse the resolution of CINT is. We show in the next section that to ensure the statistical stability of the two-point CINT (and therefore of the CINT) image, the window parameters should satisfy

[TABLE]

If scattering in the random medium is so weak that the wave components remain correlated across the aperture and bandwidth, in the sense that $X_{d}>a$ and $\Omega_{d}>B$, then we can remove the windowing in (13) to obtain

[TABLE]

and \mathbb{E}\big{[}{\cal I}_{\rm SAR}({\itbf y}^{S})\big{]} is approximately given by its expression (37) in homogeneous medium.

We are interested in strong scattering in the random medium, where $X_{d}\ll a$ and $\Omega_{d}\ll B$, and where the windowing in (13) is needed. The CINT image is then much blurrier than what the SAR image gives in homogeneus medium. This is the cost of statistical stability, as explained next.

4.2.2 The variance

The calculation of the variance of the two-point CINT imaging function uses the fourth order moment formula (25) for points and frequencies satisfying

[TABLE]

The center frequencies are $\omega,\omega^{\prime}\sim\omega_{o}$ and the frequency offsets satisfy $|\widetilde{\omega}|,|\widetilde{\omega}^{\prime}|\lesssim\Omega$. Similarly, ${\itbf x},{\itbf x}^{\prime}$ are center points in the aperture i.e., with range coordinate $L$ and cross-range coordinates of the order of $a$ and the spatial offsets $\widetilde{{\itbf x}}=(0,\widetilde{x}_{\perp})$ and $\widetilde{{\itbf x}}^{\prime}=(0,\widetilde{x}_{\perp}^{\prime})$ satisfy $|\widetilde{x}_{\perp}|,|\widetilde{x}_{\perp}^{\prime}|\lesssim 2X$. We are interested in the choice (58) of the window parameters, where the two-point CINT imaging function is statistically stable, as shown below. Definitions (28), (21) and the assumption (22) give that $X<X_{d}\ll\ell_{c}$, and after long but straightforward calculations we obtain that in this regime the variance takes the simple form

[TABLE]

We conclude that the two-point CINT imaging function gives statistically stable results when the window parameters satisfy the relation (58) because then the coefficient of variation is smaller than one:

[TABLE]

The optimal choice of the window parameters reflects the trade-off between the stability and resolution and corresponds to $X\lesssim X_{d}$ and $\Omega\lesssim\Omega_{d}$, as stated in [1, 2].

Remark 4.1**.**

We assumed a search (imaging) region ${\cal D}$ of radius less than $\ell_{c}$ in order to simplify the expressions of the second- and fourth-order moments of the Green’s function, and therefore the mean and variance of the two-point CINT imaging function. We have seen that the reflectivity function can be localized and imaged by CINT at the scales $\lambda_{o}L/X$ and $c_{o}/\Omega$ in the cross-range and range directions. This gives a consistent and relevant result because the second condition in (17) and definition (28) ensure that with $X\sim X_{d}$ and $\Omega\sim\Omega_{d}$ we have $\lambda_{o}L/X\ll\ell_{c}$ and $c/\Omega\ll\ell_{c}$ .

4.2.3 Additive noise

The effect of the additive noise on the two-point CINT imaging function is modeled by the expression

[TABLE]

Using definitions (33), (9), (11) and the paraxial approximation (32), we obtain the mean

[TABLE]

with the same constant ${\rm C}_{\rm W}$ as in (46). The covariance is calculated using the Gaussian property (48) of the noise and the result is

[TABLE]

When we set ${\itbf y}^{S}={{\itbf y}^{S}}^{\prime}$ and $\itbf z^{S}={\itbf z^{S}}^{\prime}$ in equations (62–63) we obtain that the additive noise effect on the classic CINT image consists of speckle with uniform mean

[TABLE]

and with covariance

[TABLE]

Therefore, the noise induced speckle size in the CINT image is of the order of $\lambda_{o}L/a$ in the cross-range direction and $c/B$ in the range direction. These are smaller than the cross-range resolution $\lambda_{o}L/\widetilde{X}$ and range resolution $c/\widetilde{\Omega}$ obtained in the previous section. Thus, if the noise is weak (i.e., ${\rm C}_{\rm W}$ is small), then it does not affect the CINT image. If it is moderate, then it is possible to remove the induced speckle, up to the uniform mean, using a low-pass filter on the CINT image.

However, the noise induced speckle plays a role in the two-point CINT image, because its typical size is of the same order as the resolution in the offset spatial variables. We discuss this point further in the next section.

5 High-resolution CINT imaging

With the window parameters chosen optimally, as explained above, so that the two-point CINT image is approximated by its mean, we obtain from equation (51) and definitions (15–16) that the high-resolution CINT image in the noiseless case is

[TABLE]

where $\bm{\kappa}=(\kappa_{\parallel},\kappa_{\bot})$, $\bm{\kappa}_{o}=(k_{o},0)$ and $\widehat{\rho}$ is the Fourier transform of the unknown reflectivity function. Note that

[TABLE]

the Fourier transform of the reflectivity $\rho$ modulated in range at wavenumber $k_{o}$,

[TABLE]

Equation (66) shows that we can estimate $|\widehat{\rho}_{k_{o}}(\bm{\kappa})|$ at wave vectors $\bm{\kappa}=(\kappa_{\parallel},\kappa_{\bot})$ with $|\kappa_{\parallel}|\lesssim B/c$ and $|\kappa_{\bot}|\lesssim a/(\lambda_{o}L)$. In the spatial domain, this corresponds to sampling $\rho$ on a grid of size $\lambda_{o}L/a$ in the cross-range direction and $c/B$ in the range direction. An estimate of $\rho$ on such a grid can be obtained from (66) using phase retrieval, as explained in section 5.1

In principle, there may be another way of estimating $\rho$ from the two-point CINT image, without phase retrieval. In definition (15–16) of HCINT we integrate ${\cal I}({\itbf y}^{S},{{\itbf y}^{S}}^{\prime})$ over the center points $({\itbf y}^{S}+{{\itbf y}^{S}}^{\prime})/2$ and then take the Fourier transform with respect to the offset ${\itbf y}^{S}-{{\itbf y}^{S}}^{\prime}$. We could consider instead the function

[TABLE]

and obtain from (51) and the definitions (52–53) that

[TABLE]

This can be used to determine the phase of $\widehat{\rho}$. For instance, from (70) we get

[TABLE]

so we could estimate $\nabla_{\bm{\kappa}}{\rm arg}\widehat{\rho}(\bm{\kappa}-2\bm{\kappa}_{o})$ from \nabla_{\widetilde{\bm{\kappa}}}{\rm arg}\,\widehat{\mathfrak{I}}(\bm{\kappa}-2\bm{\kappa}_{o},\widetilde{\bm{\kappa}})\big{|}_{\widetilde{\bm{\kappa}}={\bf 0}} and then integrate to get the phase. Our experience is that this approach is not stable, because it requires an estimate of $\arg\widehat{\rho}$ and its gradient over the whole domain of $\bm{\kappa}$, while $\widehat{\rho}$ may be very small in some regions of this domain.

The numerical results in section 6 are based on the phase retrieval method described next.

5.1 Phase retrieval

The goal of phase retrieval [6, 7, 8, 9, 15] is to determine a function $\eta({\itbf x})$ with Fourier transform $\widehat{\eta}(\bm{\kappa})$, such that

[TABLE]

with $|\widehat{\rho}_{k_{o}}(\bm{\kappa})|$ obtained from (66–67). Note that the restriction of the wave vector $\bm{\kappa}$ to the domain in (71) ensures that the exponential in the right-hand side of (66) is of order one and thus can be safely inverted. If the reflectivity function is known to be nonnegative valued, then we can use a phase retrieval algorithm with positivity constraints

[TABLE]

which is known to give a good reconstruction with resolution $\lambda_{o}L/a$ in the cross-range direction and $c/B$ in the range direction, up to a global shift and a symmetry with respect to the origin [6, 7]. The global shift uncertainty can be compensated by the standard CINT image, with a precision given by the classic CINT resolution. That is to say, the estimated reflectivity

[TABLE]

can be centered in the zoom area within a peak of CINT, where we seek to improve the resolution of the image.

Remark 5.1**.**

One could also apply a phase retrieval algorithm with support constraint [8] determined from the standard CINT image. In fact, it should be possible to apply the two-step approach suggested in [9] to extract a complex valued $\rho$ from $|\widehat{\rho}_{k_{o}}|$ and the low resolution standard CINT image. Moreover, if the support of $\rho$ is small, i.e. a few points, then one can use a phase retrieval algorithm with sparsity contraints [15]. The latter is more sensitive to noise than the phase retrieval with positivity constraints used in our numerical simulations.

5.2 Additive noise

The additive noise contribution to the HCINT imaging function is, in the mean,

[TABLE]

where $|\cal D|$ is the area of the imaging region. Note that (74) peaks at $\bm{\kappa}=0$, just as (66), and has similar decay. This makes it difficult to filter out the noise effect, which will impede the high-resolution imaging via phase retrieval when the noise is strong enough.

6 Numerical results

In this section we use numerical simulations to illustrate the performance of the HCINT imaging method. We consider a reflectivity $\rho$ supported at four identical point-like scatterers, as in Fig. 2. The data (1) are generated with the random travel time model***In the analysis we assumed that $\ell_{c}\ll L$ in order to use the central limit theorem and obtain Gaussian random travel time fluctuations ${\cal T}_{\mu}$. In the numerical simulations we use a Gaussian $\mu$, which means that ${\cal T}_{\mu}$ is Gaussian, no matter the ratio $L/\ell_{c}$. described in section 3, using a Gaussian zero-mean random process $\mu$ for the fluctuations of the wave speed, with correlation length $\ell_{c}=L$ and $\sigma=0.06$. The additive noise has standard deviation $\sigma_{\rm W}$ equal to $0\%$, $20\%$ or $40\%$ of the maximal amplitude of the returns. The bandwidth is $B/\omega_{o}=1/5$ and $N=60$. The window parameters of the CINT function are $\Omega=B/5$ and $X=a/5$. The imaging functions are calculated as defined in sections 2.1–2.3 and the phase-retrieval is done with the simple (error-reduction) algorithm proposed in [8]. There are better phase retrieval algorithms but this simple one was sufficient for our purpose.

The SAR imaging function gives a good estimate of the support of $\rho$ in the absence of noise (Fig. 3(a)) but it performs poorly in the presence of medium perturbations (Fig. 4(a)–5(a)). The standard CINT function gives a robust but low-resolution image, in the absence or in the presence of medium perturbations and additive noise (Fig. 2(b)–5(b)). The HCINT imaging function is displayed in Fig. 2(c)–5(c) and gives a high-resolution image, in the absence or in the presence of medium perturbations (Fig. 2(d)–5(d)).

However, HCINT is sensitive to strong additive noise. The central peak of ${\cal I}_{\rm HCINT}({\itbf y}^{S})$ is enhanced as predicted in section 5.2 and as seen in Fig. 5(c), where the central peak dominates the others, compared to the other (c) labeled pictures. Consequently, the phase-retrieval algorithm cannot determine the correct amplitudes of the four peaks of the function $\rho$ (Fig. 5(d)). By reducing the amplitude of the main peak of ${\cal I}_{\rm HCINT}({\itbf y}^{S})$ by $20\%$ we get a better image Fig. 5(f). In practice, there may be a way to automate this procedure, based on an estimate of the noise level. In any case, the phase-retrieval algorithm is known to be sensitive to noise.

As we stated earlier, there is ambiguity in the estimated reflectivity in the sense that $\rho({\itbf x})$ and $\rho({\itbf x}-{\itbf x}_{\star})$ cannot be distinguished, for arbitrary ${\itbf x}_{\star}$. In the figures we display the results given by the phase retrieval algorithm, but the reconstructions could be shifted by hand to the center of the peak of the CINT image. Since $\rho({\itbf x})$ and $\rho(-{\itbf x})$ have the same modulus of the Fourier transform, it is impossible to distinguish them with phase retrieval. This is not an issue in Fig. 2(d)–5(d) because the true reflectivity is symmetric with respect to the origin.

7 Summary

We introduced a novel interferometric imaging method for high-resolution synthetic aperture imaging of the reflectivity of a remote region, when the waves propagate through scattering random media. The method builds on the coherent interferometric (CINT) approach which uses empirical cross-correlations calculated over carefully chosen data sets in order to mitigate the distortion of the wave caused by scattering. This mitigation comes at the expense of the resolution. The algorithm introduced in this paper is based on a modified version of the CINT method, where the imaging scene is sampled at pairs of points. It shows how to use this modified imaging function to estimate the modulus of the Fourier transform of the unknown reflectivity function. The image of the reflectivity is obtained from this estimate using a phase retrieval algorithm, and the resolution is comparable to that of imaging through known and non-scattering media. The imaging method inherits the robustness of CINT with respect to the uncertainty of the random medium. However, the phase retrieval part is sensitive to noise.

In its current form, the imaging method is computationally intensive. However, the principle should be of interest for synthetic aperture radar and other imaging modalities, and it should be possible to find more efficient implementations that limit the numerical overburden and can be used in three-dimensional setups.

Acknowledgements

This research is supported in part by the Air Force Office of Scientific Research under award number FA9550-18-1-0131.