Kirchhoff migration without phases

0.54 MB·

by  Patrick Bardsley

KIRCHHOFF MIGRATION WITHOUT PHASES
PATRICKBARDSLEYANDFERNANDOGUEVARAVASQUEZ
6
1
0
2 Abstract. We present a simple, frequency domain, preprocessing step to
Kirchhoffmigrationthatallowsthemethodtoimagescattererswhenthewave
n
field phase information is lost at the receivers, and only intensities are mea-
a
sured. Theresultingimagingmethoddoesnotrequireknowingthephasesof
J
theprobingfieldormanipulatingthephaseofthewavefieldatthereceivers.

8 Inaregimewherethescatteredfieldissmallcomparedtotheprobingfield,the
problemofrecoveringthefull-waveformscatteredfieldfromintensitydatacan

] be formulated as an embarrassingly simple least-squares problem. Although
A
this only recovers the projection (on a known subspace) of the full-waveform
N scatteredfield,weshowthat,forhighfrequencies,thisprojectiongivesKirch-
hoffimagesasymptoticallyidenticaltotheimagesobtainedwithfullwaveform
.
h data. OurmethodcanalsobeusedwhenthesourceismodulatedbyaGauss-
t ianprocessandautocorrelationsaremeasuredatanarrayofreceivers.
a
m AMS classification numbers: 35R30, 78A46
Keywords: intensity-only imaging, correlation-based imaging, migration
[

1 1. Introduction
v
7 Imaging scatterers in a homogeneous medium from full waveform data is well
6
understood. Themediumisprobedwithwavesemanatingfromoneormoresources
6
and the reflections from scatterers in the medium are recorded at one or more
2
0 receivers. An image of the scatterers can be formed from these recordings by using
. classicimagingmethodssuchas Kirchhoffortraveltime migration(seee.g.[3]); or
1
MUSIC (see e.g. [6]). Both Kirchhoff migration and MUSIC rely on full-waveform
0
6 measurements at the receivers to form an image. Here we work in the frequency
1 domain and we assume that phase information is lost and only intensities can
: be measured at the receivers. To be more precise, if u (ω) is the wave field at
v (cid:98)r
i receiver r and angular frequency ω, we can only measure the intensity u(cid:98)r(ω)2.
X | |
Intensity measurements arise in a variety of physical problems, e.g. when the
r response time of a receiver is much larger than a typical wave period. Such is the
a
case in optical coherence tomography [23, 22] and diffraction tomography [14, 15].
In these situations, intensities are much easier and more cost-effective to measure
than full waveform data.
The setup we analyze consists of one source and N receivers, all of known lo-
cation. The receivers can only record intensities and only the source intensity is
known. Ifthescatteredfieldissmallcomparedtotheprobingfield(atthereceivers)
then the scattered field projected onto a known subspace can be found from the
intensity data by solving an underdetermined real least squares problem of size
Mathematics Department, University of Utah, 155 S 1400 E RM 233, Salt Lake City
UT 84112-0090.
E-mail addresses: [email protected], [email protected].
1
2 KIRCHHOFF MIGRATION WITHOUT PHASES
N 2N (per frequency sample). This system is underdetermined because we are
×
trying to use the intensity data, i.e. N real measurements, to recover the scattered
field, i.e. N complex or 2N real numbers. Fortunately, a stationary phase argu-
ment shows that the error made by projecting the scattered field does not affect
Kirchhoff imaging (for high frequencies). Moreover the least-squares problem is
typicallywell-conditionedanditssolutionisembarrassinglysimple: itmerelycosts
about 2N complex operations (additions or multiplications). Hence our method is
comparable in computational cost to Kirchhoff migration. Well-known resolution
studies for Kirchhoff migration can also be used for our method.
1.1. Related work for imaging with intensities. One way of imaging with
intensitiesiscalled“phaseretrieval”andconsistsoffirstrecoveringthephasesfrom
theintensitydataandthenusingthereconstructedfieldtoimage. Examplesofthis
approach include using intensity measurements at two different planes to recover
full waveform measurements at a single plane [14], iterative approaches e.g. [16,
10,19,8]andusingdifferentialidentitiestorelateintensitymeasurementswithfull
waveformdata(e.g. [29,17]). Othermethodstreatintensity-onlymeasurementsas
noisy measurements of full waveform data (e.g. [9]) or use optimization techniques
to fit assumed models of scatterers to measured intensity data (e.g. [27]).
The problem of imaging a few point scatterers can be reformulated as a convex
optimization problem involving low rank matrices [5, 4, 31]. Alternatively the
polarization identity 4Re(u∗v) = u + v 2 u v 2,u,v CN, and linear
(cid:107) (cid:107) − (cid:107) − (cid:107) ∈
combinations of single source experiments can be used to recover dot products of
two single source experiments [20]. Notice that recovering the dot product u∗v
from the polarization identity requires multiplying u and v by 1 or ı. This
± ±
requiresmanipulationofthesourcephases,bye.g. introducingdelays. MUSICcan
thenbeusedtoimagewiththisquadraticfunctionalofthefullwaveformdata[20].
Insteadofdirectlycontrollingthesourcephases,in[1]weusetwosourcesthatsend
exactly the same signal from two different locations. The problem of recovering
the full-waveform scattered field as measured at one receiver location is formulated
as a least-squares problem, which is analyzed in [1]. The least-squares systems are
typically2N 2N andthescatteredfieldcanberecovereduptoaonedimensional
×
nullspacethatdoesnotaffecttheKirchhoffimages. Incontrast,thepresentmethod
requires less measurements (only N), no pairwise illuminations (remark 1) and the
least-squares systems are trivial to solve and (usually) well-conditioned.
1.2. Imaging with correlations. Correlations are used for imaging when the
sources are not well known. For example in seismic imaging the sources may have
unknown locations [24, 26, 25] and may even consist of ambient noise [11, 12, 13].
Fortunately, the correlation of recordings at two points contains information about
the Green function of the medium between the two points [11], which can be used
to form an image of the medium. Correlations are also used in radar imaging [7],
since the operating frequencies make it impractical to measure the phases at the
receivers. Infactevenstochasticsignalscanbeusedinplaceofdeterministicsignals
for the probing fields [28, 30].
As in [1] we observe that autocorrelations (i.e. correlating the recorded signal
withitself)areequivalenttointensitymeasurements(bytheWiener-Khinchinthe-
orem, see e.g. [18]). Therefore the method we present here can be applied to the
casewhereonlytheautocorrelation(orthepowerspectrum)ofthesourceisknown
KIRCHHOFF MIGRATION WITHOUT PHASES 3
x
1
A
x~
r
ρ(y~)
x
3
x
2
x~
s
Figure 1. Physical setup for an array of receivers with a wave
A
source located at x(cid:126) . The scatterer is characterized by the com-
s
pactly supported function ρ(y(cid:126)) (in yellow).
and autocorrelations (or power spectra) are measured at the receivers. Since cor-
relations are robust to additive noise, we expect our method to work in low signal
to noise ratio situations, as we illustrate with numerical experiments.
1.3. Contents. Thephysicalsetupandnotationsweusearedescribedin 2.1and
§
we briefly review Kirchhoff migration in 2.2. Using the Born approximation, we
§
formulate the problem of recovering the full wave scattered field at the array from
intensity-onlymeasurementsasalinearleastsquaresproblem 2.3. In 3weanalyze
§ §
and solve the least squares problem and show that its solution can be used with
Kirchhoff migration. We extend this imaging method to stochastic illuminations
andautocorrelationmeasurementsin 4. Numericalexperimentsforanopticregime
§
are provided in 5 and we conclude with a discussion in 6.
§ §
2. Wave propagation and intensity-only measurements
HereweintroducethesetupweworkwithandbrieflyrecallKirchhoffmigration.
Hereinafter we use the Fourier transform convention for functions of time t:
(cid:90) ∞ 1 (cid:90) ∞
(1) f(cid:98)(ω)= dtf(t)eıωt, f(t)= dωf(cid:98)(ω)e−ıωt, for f(t), f(cid:98)(ω) L2(R).
2π ∈
−∞ −∞
2.1. Wave propagation in a homogeneous medium. The physical setup we
considerisdepictedinfigure1. Weprobethemediumwithapointsourcelocatedat
x(cid:126) . Waves are recorded on an array of receivers x(cid:126) =(x ,0) for r =1,...,N,
s r r
where Rd 0 and d = 2,3 is the dimension. We us∈eAthe notation x for
the firsAt d⊂ 1 c×om{po}nents of a vector x(cid:126) Rd. For simplicity we consider a linear
− ∈
array in 2D or a square array in 3D, i.e. =[ a/2,a/2]d−1 0 , however other
A − ×{ }
shapes may be considered. We impose only mild conditions on the positions of
the source and receivers, in particular that the source is not in the array. We
assume the medium contains scatterers with reflectivity ρ(y(cid:126)) with supp(ρ) = ,
R
and background wave velocity c .
0
The total field arriving at the receiver location x(cid:126) from frequency modulation
r
f(cid:98)(ω) at the source location x(cid:126)s is
(2) u(cid:98)(x(cid:126)r,x(cid:126)s,ω)=G(cid:98)(x(cid:126)r,x(cid:126)s,ω)f(cid:98)(ω),
where G(cid:98) is the Green’s function for the (inhomogeneous) medium in the frequency
domain. We assume the scatterers are weak so that multiple scatterings may be
4 KIRCHHOFF MIGRATION WITHOUT PHASES
neglected and by the Born approximation
(cid:90)
(3) G(cid:98)(x(cid:126)r,x(cid:126)s,ω) G(cid:98)0(x(cid:126)r,x(cid:126)s,ω)+k2 dy(cid:126)ρ(y(cid:126))G(cid:98)0(x(cid:126)r,y(cid:126),ω)G(cid:98)0(y(cid:126),x(cid:126)s,ω),
≈
R
where G(cid:98)0 is the Green’s function for the Helmholtz equation:

 4ıH0(1)(k|x(cid:126) −y(cid:126)|), d=2,
(4) G(cid:98)0(x(cid:126),y(cid:126),ω)=
 exp4(πıkx(cid:126)|x(cid:126) −y(cid:126)y(cid:126)|), d=3.
| − |
Here k =ω/c is the wave number.
0
We express the total fields received on the array with linear algebra notation as
u(x(cid:126)r,x(cid:126)s,ω)=eTr(cid:0)g0(x(cid:126)s,ω)+p(x(cid:126)s,ω)(cid:1)f(cid:98)(ω), for r =1,...,N,
where the vector g is the vector of direct arrivals (or incident field) at the array
0
(cid:104) (cid:105)T
(5) g0(x(cid:126)s,ω)= G(cid:98)0(x(cid:126)1,x(cid:126)s,ω),G(cid:98)0(x(cid:126)2,x(cid:126)s,ω), ,G(cid:98)0(x(cid:126)N,x(cid:126)s,ω) ,
···
and the array response vector (or scattered field at the array) is
(cid:90)
(6) p(x(cid:126)s,ω)=k2 dz(cid:126)g0(z(cid:126),ω)G(cid:98)0(x(cid:126)s,z(cid:126),ω)ρ(z(cid:126)).
R
2.2. Kirchhoff migration. When full waveform measurements are available, i.e.
u(x(cid:126) ,x(cid:126) ,ω) is known for r = 1,...,N, the scattered field p can be obtained from
r s
the total field at the array, g0 and f(cid:98)(ω). The scatterers in the medium can be
imaged with Kirchhoff migration applied to p, which for a single frequency ω has
the form:
(7) ΓKM(cid:2)p,ω(cid:3)(y(cid:126))=G(cid:98)0(x(cid:126)s,y(cid:126),ω)g0(y(cid:126),ω)∗p(x(cid:126)s,ω).
Herey(cid:126) isapointintheimage. TheKirchhoffmigrationfunctionalhasbeenstudied
extensively (see e.g. Bleistein et al. [3] for a review). In the cross-range (the
directionparalleltothearray),wecanexpectaresolutionofλL/awhereλ=2π/k
is the wavelength, L is the array to scatterer distance and a is array aperture.
This is the Rayleigh resolution limit. To obtain resolution in the range direction
(the direction perpendicular to the array), Γ [p,ω] needs to be integrated over
KM
a frequency band , e.g. = [ ω , ω ] [ω ,ω ]. In this case we can
max min min max
B B − − ∪
expect the resolution to be c /(ω ω ).
0 max min
−
2.3. Intensity-only measurements. Using the illumination f(cid:98)(ω) at the source
location x(cid:126) , the intensity-only measurement of the wave field at x(cid:126) is:
s r
∈A
(cid:104) (cid:105)
(8) |u(x(cid:126)r,x(cid:126)s,ω)|2 =|f(cid:98)(ω)|2eTr (g0+p)(cid:12)(g0+p) ,
wheretheoperator denotesthecomponentwiseorHadamardproductoftwovec-
tors and e N is(cid:12)the standard orthonormal basis of RN. Our objective is to find
{ r}r=1
asmuchaswecanaboutpfromthevectorofmeasurements[u(x(cid:126) ,x(cid:126) ,ω)2] .
r s r=1,...,N
| |
This is done by linearization, so we need to assume that the scattered field p is
small compared to the direct arrival g at the array.
0
Assumption 1. The position of the receivers, the source and the reflectivity are
such that eTp eTg , r =1,...,N.
| r |(cid:28)| r 0|
KIRCHHOFF MIGRATION WITHOUT PHASES 5
This assumption is satisfied e.g. if the reflectivity is sufficiently small and the
source x(cid:126) is near the receiver array (as is shown in figure 1). With assumption 1
s
wecanneglectquadratictermsinptoapproximatetheintensitymeasurements(8)
by a vector d(x(cid:126) ,ω) defined by
s
|u(x(cid:126)r,x(cid:126)s,ω)|2 ≈eTrd(x(cid:126)s,ω)≡|f(cid:98)(ω)|2eTrRe[g0(cid:12)(g0+2p)].
This is not, strictly speaking, a linear system for p CN since u Re(u) is not a
linear mapping from CN to CN. However we can w∈rite an underd→etermined linear
system for the real and imaginary parts of g +p as follows
0
(cid:20) (cid:21)
Re(g +2p)
(9) |f(cid:98)(ω)|2M(x(cid:126)s,ω) Im(g00+2p) =d(x(cid:126)s,ω),
where the matrix M(x(cid:126) ,ω) RN×2N is given by
s
∈
(cid:2) (cid:0) (cid:1) (cid:0) (cid:1)(cid:3)
(10) M(x(cid:126) ,ω)= diag Re(g ) diag Im(g ) .
s 0 0
We give in the next section an explicit solution to the least squares problem (9).
3. Migrating a least-squares estimate of the scattered field
The first step in our imaging method consists of a cheap least-squares prepro-
cessing step that gives an approximation to the array response vector ( 3.1). The
§
second step is to migrate this approximation with standard Kirchhoff migration
3.2. Crucially we show in theorem 1 that the mistake we make by using this
§
approximation of the array response vector does not affect the Kirchhoff images.
3.1. Recovering a projection of the array response vector. We start by
finding a simple and explicit expression to the pseudoinverse of the matrix M that
we obtained from linearizing the problem of finding the real and imaginary parts
of the array response vector p. This can be used to recover from the data d the
orthogonal projection of [Re(p)T,Im(p)T]T onto a known N dimensional subspace
(that depends only on g ). Moreover the process is well conditioned.
0
First notice that the matrix M is full-rank. Indeed a simple calculation gives
thatMMT =diag(g g ). Thismatrixisclearlyinvertiblebecauseitisadiagonal
0(cid:12) 0
matrix with the moduli of 2D or 3D Green functions on the diagonal. Hence the
Moore-Penrose pseudoinverse M† can be written explicitly
(cid:20) (cid:0) (cid:1)(cid:21)
diag Re(g )
(11) M† =MT(MMT)−1 = diag(cid:0)Im(g00)(cid:1) diag(g0(cid:12)g0)−1.
WecanuseM† toseewhatinformationaboutpwecanrecoverfromtherighthand
side d in the least-squares problem (9). Since M has an N dimensional nullspace,
we can only expect to recover the orthogonal projection of [Re(p)T,Im(p)T]T onto
range(MT) = (null(M))⊥. This projection has a simple form when we write it in
CN, as can be seen in the next proposition.
Proposition 1. Provided f(cid:98)(ω)2 =0, the intensity measurements d determine
| | (cid:54)
(12) p p+(g )−1 g p,
(cid:101)≡ 0 (cid:12) 0(cid:12)
where the inverse of a vector is understood componentwise. Moreover p can be
(cid:101)
obtained in about 2N complex operations from d with
(13) p(cid:101)=|f(cid:98)(ω)|−2(g0)−1(cid:12)d−g0.
6 KIRCHHOFF MIGRATION WITHOUT PHASES
Proof. Since we use the first (resp. last) N rows of M† to recover the real (resp.
imaginary) part of a vector in CN, it is convenient to consider the matrix
(cid:2)I iI(cid:3)M† =diag(g )diag(g g )−1 =diag(g )−1,
0 0(cid:12) 0 0
where I is the N N identity matrix. To see what information about p we can
×
recover from the right hand side d in the least-squares system (9) we can evaluate:
(cid:20) (cid:21)
f(cid:98)(ω)−2(cid:2)I iI(cid:3)M†d=diag(g )−1M Re(g0+2p)
| | 0 Im(g0+2p)
=diag(g )−1[Re(g )2+Im(g )2+2Re(g )Re(p)+2Im(g )Im(p)]
0 0 0 0 0
=g +diag(g )−1(g p+g p)
0 0 0(cid:12) 0(cid:12)
=g +p+g (g )−1 p=g +p.
0 0(cid:12) 0 (cid:12) 0 (cid:101)
Hence we can get p from the intensity data d with essentially N complex multipli-
(cid:101)
cations and N complex additions.
(cid:3)
Anaturalquestiontoaskiswhetherwecanobtainpinastablemannerfromd.
(cid:101)
ThiscanbeansweredbylookingattheconditioningofM,i.e. theratioofthelargest
singular value σ of M to σ , the smallest one. These are easily obtained from the
1 N
squarerootsoftheeigenvaluesofthediagonalmatrixMMT =diag(g g ). Hence
0(cid:12) 0
theconditioningofMistheratioofthelargesttothesmallestmodulioftheentries
of g :
0
(14) cond(M)=mmainxrr(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)HH0(0(11))((kk||x(cid:126)x(cid:126)(cid:12)rr−−x(cid:126)x(cid:126)ss||))(cid:12)(cid:12)(cid:12)(cid:12), for d=2,
maxr(cid:12)(cid:12)x(cid:126)r−x(cid:126)s(cid:12)(cid:12), for d=3.
minr(cid:12)x(cid:126)r x(cid:126)s(cid:12)
−
In figure 2 we show the condition number of M(x(cid:126) ,ω) plotted over an optical
s
frequencyband. Theexperimentalsetupisthatgivenin 5. Theconditionnumber
§
(14) is clearly independent of frequency for d = 3 and for d = 2 we have the
approximation for high frequencies:
max x(cid:126) x(cid:126) 1/2
condM(x(cid:126) ,ω)= r| r− s| (1+ (1/ω)), as ω .
s min x(cid:126) x(cid:126) 1/2 O →∞
r r s
| − |
This approximation follows from the Hankel function asymptotic (see e.g. [21])
(cid:114)
2
H(1)(t)= exp[ı(t π/4)](1+ (1/t)), as t .
0 πt − O →∞
ThustheconditioningofMisdeterminedbytheratiooflargesttosmallestsource-
to-receiver distances.
3.2. Kirchhoff migration. We now show that migrating the recovered data p
(cid:101)
(12) using Γ gives essentially the same image as migrating the true data p. We
KM
establish this result by means of a stationary phase argument but in order to do
this, we need the following assumption on the location of the source x(cid:126) .
s
KIRCHHOFF MIGRATION WITHOUT PHASES 7
3
2
N
σ
/
σ1
1
500 600 700
FrequencyinTHz
Figure 2. Condition number of M(x(cid:126) ,ω) for d = 2 (red) and
s
d=3 (blue) for the setup given in 5
§
Assumption 2 (Geometric imaging conditions). For a scattering potential with
support contained inside an image window , we assume x(cid:126) satisfies
s
W
x(cid:126) x(cid:126) x(cid:126) y(cid:126)
r− s = r− ,
x(cid:126) x(cid:126) (cid:54) x(cid:126) y(cid:126)
r s r
| − | | − |
for r =1,...,N, and y(cid:126) .
∈W
We interpret this assumption as a restriction on the placement of our source
location x(cid:126) as follows. Fix a receiver position x(cid:126) and consider the cone
s r
(cid:26) (cid:27)
y(cid:126) x(cid:126)
(x(cid:126) )= α − r :α>0,y(cid:126) .
K r y(cid:126) x(cid:126) ∈W
r
| − |
As long as x(cid:126) / (x(cid:126) ), then we have that (x(cid:126) x(cid:126) )/x(cid:126) x(cid:126) =(y(cid:126) x(cid:126) )/y(cid:126) x(cid:126)
s r s r s r r r
∈K − | − |(cid:54) − | − |
for any y(cid:126) , i.e. assumption 2 holds for x(cid:126) . Ensuring this is satisfied for all
r
∈ W
receiver locations x(cid:126) for r = 1,...,N, we require x(cid:126) / N (x(cid:126) ). In figure 3 we
r s ∈ ∪r=1K r
illustratethisassumption. Here, thedarkblueregiondepictsthecone (x(cid:126) )while
r
K
the union of cones N (x(cid:126) ) is depicted by the light blue region. Assumption 2
∪r=1K r
simply requires x(cid:126) to be outside the light blue region.
s
A
x~
r
(x~ )
r
K
N (x~ )
∪r=1K r
W
Figure 3. Illustration of assumption 2. If x(cid:126) is outside of the
s
light blue region then (x(cid:126) x(cid:126) )/x(cid:126) x(cid:126) =(x(cid:126) y(cid:126))/x(cid:126) y(cid:126) for
r s r s r r
− | − |(cid:54) − | − |
all x(cid:126) and y(cid:126) .
r
∈A ∈W
Theorem 1. Provided assumption 2 holds, the image of the reconstructed array
response vector is
Γ (cid:2)p+(g )−1 g p,ω(cid:3)(y(cid:126)) Γ (cid:2)p,ω(cid:3)(y(cid:126)).
KM 0 (cid:12) 0(cid:12) ≈ KM
8 KIRCHHOFF MIGRATION WITHOUT PHASES
Proof. We begin by approximating the Kirchhoff migration functional (7) by an
integral over the array :
A
Γ (cid:2)(g )−1 g p(cid:3)(y(cid:126))
KM 0 (cid:12) 0(cid:12)
=G(cid:98)0(x(cid:126)s,y(cid:126),ω)g0(y(cid:126),ω)∗(cid:2)(g0(x(cid:126)s,ω))−1(cid:12)g0(x(cid:126)s,ω)(cid:12)p(x(cid:126)s,ω)(cid:3)
(15) (cid:90) (cid:90)
k2 dz(cid:126) dx C(x(cid:126) ,x(cid:126) ,z(cid:126),y(cid:126))
r s r
∼ ×
R A
exp(ık(2x(cid:126) x(cid:126) x(cid:126) z(cid:126) z(cid:126) x(cid:126) x(cid:126) y(cid:126) y(cid:126) x(cid:126) )).
r s r s r s
| − |−| − |−| − |−| − |−| − |
Here denotes equal up to a constant and C(x(cid:126) ,x(cid:126) ,z(cid:126),y(cid:126)) is a smooth real valued
s r
∼
function that collects the various G(cid:98)0 geometric spreading terms.
Now we apply the method of stationary phase (see e.g. [3]) to the integral over
. In the large wavenumber limit k , dominant contributions to the integral
A → ∞
come from stationary points of the phase, i.e. points x(cid:126) that satisfy
r
(2x(cid:126) x(cid:126) x(cid:126) z(cid:126) z(cid:126) x(cid:126) x(cid:126) y(cid:126) y(cid:126) x(cid:126) )=0.
∇x(cid:126)r | r− s|−| r− |−| − s|−| r− |−| − s|
This expression is equivalent to
(cid:18) (cid:19)
x(cid:126) x(cid:126) 1 x(cid:126) y(cid:126) x(cid:126) z(cid:126)
(16) r− s = r− + r− .
x(cid:126) x(cid:126) 2 x(cid:126) y(cid:126) x(cid:126) z(cid:126)
r s r r
| − | | − | | − |
If (16) holds then we must have
(cid:12)(cid:12) x(cid:126)r x(cid:126)s (cid:12)(cid:12)2 1(cid:12)(cid:12) x(cid:126)r y(cid:126) x(cid:126)r z(cid:126) (cid:12)(cid:12)2
(cid:12) − (cid:12) = (cid:12) − + − (cid:12)
(cid:12) x(cid:126) x(cid:126) (cid:12) 4(cid:12) x(cid:126) y(cid:126) x(cid:126) z(cid:126) (cid:12)
(17) | r− s| | r− | | r− |
x(cid:126) y(cid:126) x(cid:126) z(cid:126)
1= r− r− .
⇐⇒ x(cid:126) y(cid:126) · x(cid:126) z(cid:126)
r r
| − | | − |
Since (x(cid:126) y(cid:126))/x(cid:126) y(cid:126) and (x(cid:126) z(cid:126))/x(cid:126) z(cid:126) are both unit vectors, it follows from
r r r r
− | − | − | − |
the Cauchy-Schwarz equality that (17) holds only if z(cid:126) =y(cid:126). Thus stationary points
must satisfy
x(cid:126) x(cid:126) x(cid:126) y(cid:126)
r− s = r− ,
x(cid:126) x(cid:126) x(cid:126) y(cid:126)
r s r
| − | | − |
where y(cid:126) . By assumption 2 there are no such stationary points and therefore,
∈W
neglectingboundaryeffects,thisintegralvanishesfasterthananypolynomialpower
of ω (see e.g. [2]). (cid:3)
Remark 1. We used a similar idea in [1] to show that with multiple sources, a
single receiver, and a specific pairwise illumination scheme it is possible to image
with sole knowledge of the intensities of the wave fields at the receiver and of the
probing fields. We approached the problem by estimating the array response vec-
tor (a vector in CN) with 2N (or more) real measurements, which are essentially
the measured intensities for 2N or more different pairs of sources. The results of
3.1 and 3.2 can be modified by reciprocity to apply to the setup we considered
§ §
in [1]. Hence images similar to those in [1] can be obtained without the pairwise
illumination scheme and the number of required illuminations is reduced from 3N
to N.
KIRCHHOFF MIGRATION WITHOUT PHASES 9
4. Stochastic illuminations and autocorrelations
Our imaging method can also be used when the source is driven by a stationary
stochastic process (for which we only assume knowledge of the autocorrelation or
power spectra) and only empirical autocorrelations are measured at the receiver
locations.
To be more precise, the source at x(cid:126) is driven by f(t), a stationary mean zero
s
Gaussian process with autocorrelation function
(18) f(t)f(t+τ) =F(τ).
(cid:104) (cid:105)
Here denotes expectation with respect to realizations of f and we recall that
(cid:104)·(cid:105)
F(τ)=F( τ). In the time domain, the field recorded at x(cid:126) is
r
−
(cid:90)
1
u(x(cid:126)r,x(cid:126)s,t)= 2π dωe−ıωtG(cid:98)(x(cid:126)r,x(cid:126)s,ω)f(cid:98)(ω), for r =1,...,N,
where we assume G(cid:98) is given by the Born approximation (3).
Themeasurementsatthereceiverlocationsx(cid:126) aretheempiricalautocorrelations:
r
1 (cid:90) T
(19) ψ(x(cid:126) ,x(cid:126) ,τ)= dtu(x(cid:126) ,x(cid:126) ,t)u(x(cid:126) ,x(cid:126) ,t+τ) for r =1,...,N,
r s 2T r s r s
−T
where T is a fixed acquisition time. As shown by Garnier and Papanicolaou [11],
these measurements are independent of the acquisition time T and ergodic as we
summarize in the following proposition.
Proposition2. Assumef(t)isastationarymeanzeroGaussianprocesssatisfying
(18). The expectation of the empirical autocorrelations (19) is independent of the
acquisition time T:
ψ(x(cid:126) ,x(cid:126) ,τ) =Ψ(x(cid:126) ,x(cid:126) ,τ),
r s r s
(cid:104) (cid:105)
where
(cid:90)
1
(20) Ψ(x(cid:126)r,x(cid:126)s,τ)= 2π dωe−ıωτF(cid:98)(ω)eTr [g(x(cid:126)s,ω)(cid:12)g(x(cid:126)s,ω)],
with g g +p. Furthermore, (19) is ergodic, i.e.
0
≡
T→∞
(21) ψ(x(cid:126) ,x(cid:126) ,τ) Ψ(x(cid:126) ,x(cid:126) ,τ).
r s r s
−−−−→
Proof. The proof is a straight-forward application of [11, Proposition 4.1]. (cid:3)
Theergodicity(21)ofthispropositionguaranteesthatforsufficientlylargeacqui-
sitiontimeT ,theautocorrelationψ(x(cid:126) ,x(cid:126) ,τ)isclosetoanintensitymeasurement,
r s
i.e.
(cid:104) (cid:105)
ψ(cid:98)(x(cid:126)r,x(cid:126)s,ω)−T−→−−∞→Ψ(cid:98)(x(cid:126)r,x(cid:126)s,ω)=F(cid:98)(ω)eTr (g0+p)(cid:12)(g0+p) .
Proceeding analogously as in 2.3, we neglect the quadratic term in p:
§
Ψ(cid:98)(x(cid:126)r,x(cid:126)s,ω)≈F(cid:98)(ω)eTrRe(g0(cid:12)(g0+2p)).
Thecollectionofautocorrelationsforr =1,...,N canbeexpressed,approximately,
as
(cid:104) (cid:105) (cid:20)Re(g +2p)(cid:21)
Ψ(cid:98)(x(cid:126)r,x(cid:126)s,ω) r=1,...,N ≈d(x(cid:126)s,ω)≡F(cid:98)(ω)M(x(cid:126)s,ω) Im(g00+2p) ,
where M(x(cid:126) ,ω) RN×2N is given by (10). Therefore, the techniques developed in
s
∈
3 can be applied to image from the autocorrelation measurements (19).
§
10 KIRCHHOFF MIGRATION WITHOUT PHASES
5. Numerical experiments
We now provide 2D numerical experiments of our proposed imaging method.
The physical scalings we use correspond to an optic regime. We use the back-
ground wave velocity of c =3 108 m/s and central frequency of about 590 THz
0
×
which gives a central wavelength λ of about 509 nm. Our receiver array is a
0
A
linear array centered at the origin and consists of 501 receivers located at coor-
dinates x(cid:126) = (0, 5+(r 1)(10/500))mm for r = 1,...,501. This corresponds
r
− −
to using a 1 cm linear array of receivers spaced approximately 20µm apart. We
place the wave source at coordinate x(cid:126) = (5, 7.5)mm to guarantee assumption 2
s
−
is satisfied. We begin with experiments in the deterministic setting ( 5.1) followed
§
by an experiment in the stochastic setting ( 5.2). Lastly, we investigate situations
§
whereassumptions1and/or2areviolatedandourmethodisnotexpectedtowork
( 5.3). For all experiments, we assume 3D wave propagation for simplicity so that
§
G(cid:98)0 is given by (4) for d=3.
5.1. Deterministic illuminations. Usingtheilluminationf(cid:98)(ω) 1,wegenerate
≡
the intensity data d(x(cid:126) ,ω) using the Born approximation:
s
d(x(cid:126) ,ω)=(g +p) (g +p),
s 0 0
(cid:12)
withg andpdefinedby(5)and(6)respectively,for100uniformlyspacedfrequen-
0
cies in the frequency band [430,750] THz. This corresponds to obtaining intensity
datafor100differentmonochromaticilluminationswithwavelengthsλ [400,700]
∈
nm, equally spaced in the frequency band. Note that the quadratic term p p is
(cid:12)
present in our data, but using assumption 1 we proceed assuming d is well approx-
imated by the linear system (9).
We recover the approximate array response vector p=(g )−1 d g for each
(cid:101) 0 (cid:12) − 0
frequencyω intheangularfrequencyband ,where(2π)−1 =[430,750]THz. An
B B
image is then formed using the Kirchhoff migration functional integrated over :
B
(cid:90)
Γ [p](y(cid:126))= dωΓ [p,ω](y(cid:126)),
KM (cid:101) KM (cid:101)
B
where Γ is defined in (7). Here we consider image points y(cid:126) = (50mm+
KM
∈ W {
iλ /2.5,jλ /2.5), for i,j = 25,...,25 .
0 0
− }
For our first experiment, we consider a point reflector located at coordinate y(cid:126) =
(50,0)mm with refractive index perturbation ρ(y(cid:126))=1 10−15 (roughly equivalent
×
to a reflector of area (λ )2 and reflectivity 5978). The migrated images of the
0
true array response vector p and the recovered array response vector p are shown
(cid:101)
in figure 4a. Although we are significantly undersampling both in frequency and
on the array (recall the spacing between receivers is approx. 20µm λ /2), the
0
(cid:29)
images still exhibit the cross-range (Rayleigh) resolution estimate λ L/a 5λ
0 0
≈
and range resolution estimate c / 1λ . Our second experiment (figure 4b)
0 0
|B| ≈
uses two point reflectors located at coordinates y(cid:126) =(50mm 3λ , λ ) and y(cid:126) =
1 0 0 2
− −
(50mm+6λ ,5λ ) each with ρ(y(cid:126) ) = 1 10−15. We show an extended scatterer
0 0 i
×
(a disk) in figure 5. The disk is generated as a set of point reflectors, each with
ρ(y(cid:126) )=1 10−15 separated by λ /4.
i 0
×
5.2. Stochastic illumination. Here we image with power spectrum data d gen-
erated from a stochastic illumination as in 4. Since we work in an optic regime,
§