Source estimation with incoherent waves in random waveguides

Source estimation with incoherent waves in random waveguides

Sebastian Acosta111Baylor College of Medicine, Houston, TX 77005., Ricardo Alonso222Departamento de Matemática, PUC–Rio, Brasil. and Liliana Borcea333Department of Mathematics, University of Michigan, Ann Arbor, MI 48109.

We study an inverse source problem for the acoustic wave equation in a random waveguide. The goal is to estimate the source of waves from measurements of the acoustic pressure at a remote array of sensors. The waveguide effect is due to boundaries that trap the waves and guide them in a preferred (range) direction, the waveguide axis, along which the medium is unbounded. The random waveguide is a model of perturbed ideal waveguides which have flat boundaries and are filled with known media that do not change with range. The perturbation consists of fluctuations of the boundary and of the wave speed due to numerous small inhomogeneities in the medium. The fluctuations are uncertain in applications, which is why we model them with random processes, and they cause significant cumulative scattering at long ranges from the source. The scattering effect manifests mathematically as an exponential decay of the expectation of the acoustic pressure, the coherent part of the wave. The incoherent wave is modeled by the random fluctuations of the acoustic pressure, which dominate the expectation at long ranges from the source. We use the existing theory of wave propagation in random waveguides to analyze the inverse problem of estimating the source from incoherent wave recordings at remote arrays. We show how to obtain from the incoherent measurements high fidelity estimates of the time resolved energy carried by the waveguide modes, and study the invertibility of the system of transport equations that model energy propagation in order to estimate the source.

Key words. Waveguides, random media, transport equations, Wigner transform.

AMS subject classifications. 35Q61, 35R60

1 Introduction

We study an inverse problem for the scalar (acoustic) wave equation, where we wish to estimate the source of waves from measurements of the acoustic pressure field at a remote array of receiver sensors. The waves propagate in a waveguide, meaning that they are trapped by boundaries and are guided in the range direction, the waveguide axis, along which the medium is unbounded. Ideally the boundaries are straight and the medium does not change with range. We consider perturbed waveguides filled with heterogeneous media, where the boundary and the wave speed have small fluctuations on scales similar to the wavelength. These fluctuations have little effect in the vicinity of the source, but they are important at long ranges because they cause significant cumulative wave scattering. We suppose that the array of receivers is far from the source, as is typical in applications in underwater acoustics, sound propagation in corrugated pipes, in tunnels, etc., and study how cumulative scattering impedes the inversion.

In most setups the fluctuations are uncertain, which is why we introduce a stochastic framework and model them with random processes. The inversion is carried in only one perturbed waveguide, meaning that the array measures one realization of the random pressure field, the solution of the wave equation in that waveguide. The stochastic framework allows us to study the chain of mappings from the uncertainty in the waveguide to the uncertainty of the array measurements and of the inversion results. The goal is to understand how to process the uncertain data and quantify what can be estimated about the source in a reliable (statistically stable) manner. Statistical stability means that the estimates do not change with the realization of the fluctuations of the waveguide, which are unknown.

The problem of imaging (localizing) sources in waveguides has been studied extensively in underwater acoustics [5, 21, 19, 1]. Typical imaging approaches are matched field and related coherent methods that match the measured with its mathematical model for search locations of the source. The model is based on wave propagation in ideal waveguides and the imaging is successful when is mostly coherent. The coherent part of is its statistical expectation with respect to realizations of the random waveguide, and the incoherent field is modeled by . As the waves propagate in the random waveguide they lose coherence due to scattering by the fluctuations of the boundary and the inhomogeneities in the medium. This manifests as an exponential decay in range of the expectation , and strengthening of the fluctuations .

Detailed studies of the loss of coherence of sound waves due to cumulative scattering are given in [20, 10, 13, 18, 14] for waveguides filled with randomly heterogeneous media and in [4, 17] for waveguides with random boundaries. These waveguides are two dimensional models of the ocean, and they may leak (radiate) in the ocean floor. The problem is similar in three dimensional acoustic waveguides with bounded cross-section. We refer to [6] for wave propagation in three dimensional waveguide models of the ocean which have unbounded cross-section and random pressure release top boundary, and to [3, 22] for three dimensional electromagnetic random waveguides. In all cases the analysis of loss of coherence is based on the decomposition of the wave field in an infinite set of monochromatic waves called waveguide modes, which are special solutions of the wave equation in the ideal waveguide. Finitely many modes are propagating waves, and we may associate them with plane waves that strike the boundary at different angles of incidence and are reflected repeatedly. The remaining infinitely many modes are evanescent and/or radiating waves. The cumulative scattering in the waveguide is modeled by fluctuations of the amplitudes of the modes. When scattering is weak, as is the case at moderate distances from the source, the amplitudes are approximately constant in range, and they are determined solely by the source excitation. Scattering builds up over long ranges and the mode amplitudes become random fields with exponentially decaying expectation on range scales called scattering mean free paths.

The mode dependence of the scattering mean free paths is analyzed in [4]. It turns out that the slow modes, which correspond to plane waves that strike the boundary at almost normal incidence, are most affected by scattering. These waves have long trajectories from the source to the array, and thus interact more with the boundary and medium fluctuations. We refer to [7] for an adaptive coherent imaging approach which detects which modes are incoherent and filters them out from the measurements in order to achieve statistically stable results. See also the results in [21, 19, 25]. However, when the array is farther from the source than the scattering mean free paths of all the modes, the data is incoherent and coherent imaging methods like matched field cannot work. In this paper we assume that this is the case and study an inversion approach based on a system of transport equations that models the propagation of energy carried by the modes. This system is derived in [20, 13, 4] and is used in [8] to estimate the location of a point source in random waveguides. Here we study the inverse problem in more detail and answer the following questions: (1) How can we obtain reliable estimates of the mode energies from the incoherent pressure field measured at the array? (2) What kind of information about the source can we recover from the transport equations? (3) Can we quantify the deterioration of the inversion results in terms of the range offset between the source and the array?

We begin in section LABEL:sect:setup with the mathematical formulation of the inverse problem, and recall in section LABEL:sect:stat the model of the random wave field derived in [20, 13, 4, 9]. The main results of the paper are in sections LABEL:sect:transpIm and LABEL:sect:inverse. We motivate there the inversion based on energy transport, and describe the forward mapping from the source to the expectation of the time resolved energy carried by the modes. We show how to calculate this energy from the incoherent array data, and describe how to invert approximately the transport equations. The results quantify the limited information that can be recovered about the source. We end with a summary in section LABEL:sect:summary.

2 Formulation of the problem

array at range
Figure 2.1: Schematic of the problem setup. The source emits a signal in a waveguide and the wave field is recorded at a remote array. The perturbed waveguide has fluctuating boundaries and is filled with a medium with fluctuating wave speed.

We limit our study to two dimensional waveguides with reflecting boundaries modeled by pressure release boundary conditions. This is for simplicity, but the results extend to other boundary conditions and to leaky and three dimensional waveguides, as discussed in section LABEL:sect:summary. We illustrate the setup in Figure LABEL:fig:Schem, and introduce the system of coordinates with range originating from the center of the source. The waveguide occupies the domain

where the cross-range takes values between the bottom and top boundaries modeled by and . The source has an unknown density which is compactly supported in , near , and emits a signal which is a pulse of support of order around , modulated by an oscillatory exponential


We introduce the bandwidth in the argument of the pulse to emphasize that the Fourier transform of the signal is supported in the interval around the central frequency ,


The array is a collection of receivers that are placed close together in the set

at range from the source, where is an interval called the array aperture. The receivers record the acoustic pressure field modeled by the solution of the acoustic wave equation


with pressure release boundary conditions


and initial condition for . Here is the sound speed.

The inverse problem is to determine the source density from the array data recordings . We model them by


using a recording time window centered at and of duration . We can take any continuous, compactly supported , but we assume henceforth that it equals one in the interval and tapers quickly to zero outside. We also approximate the array by a continuum aperture in the interval , and use the indicator function which equals one when and zero otherwise.

2.1 The random model of perturbed waveguides

In ideal waveguides the sound speed is modeled by a function that is independent of range and the boundaries are straight, meaning that and , a constant. The sound speed in the perturbed waveguide has fluctuations around and the boundaries and fluctuate around and . The fluctuations are small, with amplitude quantified by a positive dimensionless parameter . It is used in [20, 13, 4, 9] to analyze the pressure field at properly scaled long ranges where scattering is significant, in the asymptotic limit .

We take constant for simplicity, to write explicitly the mode decomposition, but the results extend easily to cross-range dependent . The perturbed sound speed is modeled by


where is a mean zero random process that is bounded almost surely, so that the right hand side in (LABEL:eq:fm5) remains positive. We assume that is stationary and mixing in range, meaning in particular that the auto-correlation


is absolutely integrable in the third argument over the real line. The process is normalized by and

where is the correlation length, the range offset over which the random fluctuations become statistically decorrelated. It compares to the central wavelength as . The scaling by the same of the cross-range in (LABEL:eq:fm5) means that the heterogeneous medium is isotropic, but we could have as well, without changing the conclusions. The amplitude of the fluctuations is scaled by which equals for a random medium and zero for a homogeneous medium.

We model similarly the boundary fluctuations


using two mean zero, stationary and mixing random processes and , that are bounded almost surely and have integrable autocorrelation and . We assume that , and are independent111If the random processes are not independent the moment formulae in this paper must be modified. Their derivation is a straightforward extension of the analysis in [20, 13, 4]. and use the same correlation length to simplify notation, but the results hold for any scales and of the order of . For technical reasons related to the method of analysis used in [4] we also assume that the processes and have bounded first and second derivatives, almost surely. Less smooth boundary fluctuations are considered in [17]. The boundary fluctuations are scaled by and which can be , or they may be set to zero to study separately the scattering effects of the medium and the boundary.

The theory of wave propagation in waveguides with long range correlations of the random fluctuations of is being developed [16], and our results are expected to extend (with modifications) to such settings. The case of turning waveguides with smooth and large variations of the boundaries, on scales that are comparable to , is much more difficult. The analysis of wave propagation in such waveguides is quite involved [24, 11, 2] and the mapping of random fluctuations of the sound speed to is not understood in detail, although it is considered formally in [22].

3 Cumulative scattering effects in the random waveguide

We write the solution of the wave equation (LABEL:eq:fm2)-(LABEL:eq:fm3) as


where is the wave field due to a point source at , emitting the signal defined in (LABEL:eq:fm1), and is the compact support of the source, which lies near . The points in (LABEL:eq:st1) are at range , for all .

It follows from [20, 13, 4, 9] that is a linear superposition of propagating and evanescent waves, called waveguide modes


The modes are special solutions of the wave equation in the ideal waveguide, and can be obtained with separation of variables. There are propagating modes


with index denoting forward going and backward going, and infinitely many evanescent modes


They are defined by the complete and orthonormal set of eigenfunctions of the symmetric linear operator with homogeneous Dirichlet boundary conditions at and , where . Because is constant we can write


and note explicitly how are associated with monochromatic plane waves that travel in the direction of the slowness vectors and strike the boundaries where they reflect according to Snell’s law. The mode wavenumbers are denoted by , and are determined by the square root of the eigenvalues of the operator


The of the propagating modes correspond to the first eigenvalues which are positive, where


and denotes the integer part.

We assume for simplicity that the bandwidth is not too large222In applications of imaging in open environments large bandwidths are desired for improved range resolution. In ideal waveguides good images can be formed with small bandwidths because the modes give different angle views of the support of the source. In random waveguides we may benefit from a large bandwidth, as explained in section LABEL:sect:I6. Such bandwidths may be divided in smaller sub-bands to which we can apply the analysis in this paper., so that there is the same number of propagating modes for all the frequencies of the pulse, and drop the dependence of on . We also suppose that there are no standing waves, meaning that are bounded below by a positive constant, for all .

The cumulative scattering effects in the random waveguide are modeled by the mode amplitudes and , which are random fields. In ideal waveguides the amplitudes are constant in range for


They depend on the cross-range in the support of the source, and the second equation in (LABEL:eq:st8) complies with the wave being outgoing. In random waveguides the mode amplitudes satisfy a coupled system of stochastic differential equations driven by the random fluctuations , and . They are analyzed in detail in [20, 13, 4, 9] and the result is that they are approximately the same as (LABEL:eq:st8)-(LABEL:eq:st10) for range offsets . This motivates the long range scaling


where cumulative scattering becomes significant. The evanescent modes may be neglected at such ranges333Note that although the evanescent modes do not appear explicitly in (LABEL:eq:st12), they affect the amplitudes of the propagating modes. This amplitude coupling is taken into account in the analysis in [20, 13, 4, 9] and thus in the results of this paper., and we use a further approximation that neglects the backward going waves to write


The forward scattering approximation holds for , and is justified by the fact that the backward mode amplitudes have very weak coupling with the forward ones, for autocorrelations of the fluctuations that are smooth enough in [20, 13, 4, 9]. We refer to [14] for the analysis of wave propagation that includes both the forward and backward going modes, but for the purpose of this paper it suffices to use (LABEL:eq:st12).

Let us write the amplitudes using the random propagator , which maps the amplitudes (LABEL:eq:st8) near the source at range , to those at the array


The propagator is analyzed in [20, 13, 4] in the asymptotic limit . It converges in distribution to a Markov diffusion with generator computed explicitly in terms of the autocorrelations of the random fluctuations. Thus, we can rewrite (LABEL:eq:st13) as


with symbol denoting approximate in distribution. It means that we can approximate the statistical moments of using the right hand side in (LABEL:eq:st14), with an error in the limit .

3.1 Data model

The data model follows from (LABEL:eq:fm4), (LABEL:eq:st1) and (LABEL:eq:st12)


where is the Fourier transform of the recording window and is given by (LABEL:eq:st13)-(LABEL:eq:st14). We take henceforth the bandwidth


which is small with respect to the center frequency. We ask that because the travel time of the modes is of order , and we need a pulse of much smaller temporal support in order to distinguish the arrival time of different modes. That is also needed for the statistical stability of the inversion, as we explain later. The choice is for convenience444For the results are similar, but higher powers of enter in the phase, and they change the shape of the pulse carried by the modes., because it allows us to linearize the phase in (LABEL:eq:model) as

with small error of order . When we use this approximation in (LABEL:eq:model) we see that in ideal waveguides where the modes propagate with range speed


In random waveguides only the expectation (coherent part) of propagates at speed (LABEL:eq:model4), but the energy of the mode is transported at different speed, as described in section LABEL:sect:I1. In any case, we note that the wavenumbers decrease monotonically with , so the first modes are faster as expected, because they take a more direct path from the source to the array. For example, in the case the slowness vectors of the plane waves associated with the first mode are almost parallel to the range direction, and the speed (LABEL:eq:model4) is approximately equal to . For the last modes the slowness vectors are almost orthogonal to the range direction and the speed is much smaller than .

It is natural to choose the duration of the recording window to be much longer than that of the pulse . We shall see in section LABEL:sect:transpIm that in fact we need to be at least of the order of the travel time of the waves in order for the incoherent imaging method to work. Thus, we let


with . We also assume that is a continuous function to simplify (LABEL:eq:model) slightly using the approximation


3.2 Loss of coherence

To compute the coherent part of the data model, we recall from [20, 13, 4, 9] the expectation of the limit propagator


where is the Kronecker delta symbol and the approximation is due to the fact that is much smaller than . Although the mean propagator is a diagonal matrix as in ideal waveguides, where it is the identity, its entries are exponentially damped in on scales , the scattering mean free path of the modes. There is also an anomalous phase accumulated on the mode dependent scales .

The scales and are defined in [7, equations (3.19),(3.28),(3.31)] and depend on the frequency and the autocorrelations , and of the fluctuations. Of particular interest in this paper are the scattering mean free paths because they give the range scale on which the modes randomize. The magnitude of the expectation (coherent part) of the mode amplitudes follows from (LABEL:eq:st14) and (LABEL:eq:st15)


where is the initial condition of at , equal to the amplitude (LABEL:eq:st8) in ideal waveguides. The exponential decay in (LABEL:eq:st16NN) is not caused by attenuation in the medium. The wave equation conserves energy, and we state in the next section that does not tend to zero. The meaning of the decay in (LABEL:eq:st16NN) is the randomization (loss of coherence) of the th mode due to scattering. It says that beyond scaled ranges the mode becomes incoherent i.e., the random fluctuations of its amplitude dominate its expectation.

The scattering mean free paths are given by


in terms of


where is the autocorrelation of the stationary process


This is for in (LABEL:eq:fm5) and (LABEL:eq:fm7), and for statistically independent random processes , and . The hat denotes the Fourier transform of the autocorrelations, which is non-negative by Bochner’s theorem.

To compare the scattering effects in the random medium with those at the boundary, we plot with solid line in Figure LABEL:fig:MatrixU for the case and Figure LABEL:fig:MatrixUBound for . In the first case we keep only the last term in (LABEL:eq:st22) and in the second case we keep the second term. The setup of the simulations is explained in the numerics section LABEL:sect:I3. Here we note two important facts displayed by the plots: The scales decrease monotonically with , and their mode dependence is much stronger in the random boundary case. This is intuitive once we recall that the first modes are waves that travel along a more direct path from the source to the array. These waves interact with the random boundary only once in a while and thus randomize on longer range scales than the slower modes. For example is more than a hundred times longer than in Figure LABEL:fig:MatrixUBound. The slow modes are waves that reflect repeatedly at the boundary and travel a long way in the waveguide as they progress slowly in range. They randomize on small range scales for both random boundary and medium scattering. However, the medium scatter leads to more dramatic loss of coherence as illustrated in Figure LABEL:fig:MatrixU, where all but the last modes have similar scattering mean free paths which are shorter than in Figure LABEL:fig:MatrixUBound.

The goal of this paper is to analyze what can be determined about the source of waves from measurements made at ranges , where


for all i.e., all the modes are incoherent. No coherent method can work in this regime, so we study an incoherent inversion approach based on the transport of energy theory summarized in the next two sections.

3.3 Statistical decorrelation

Since the wave equation is not dissipative, we have the conservation of energy relation [20, 13, 4, 9]


where the approximation is with an error as , due to the neglect of the backward going and evanescent waves. Thus, some second moments of the mode amplitudes remain finite, and can be used in inversion. To decide if we can estimate them reliably from the incoherent data, we need to know how the waves decorrelate. Statistical decorrelation means that the second moments of the amplitudes are equal approximately to the product of their expectations, which is negligible by (LABEL:eq:st16).

The two frequency analysis of the propagator is carried out in [13, 4], and the result is that the waves are decorrelated for frequency offsets Such small offsets are enough to cause the waves to interact differently with the random fluctuations over ranges , thus giving the statistical decorrelation. This result is important because it says that we can estimate those second moments of the amplitudes that do not decay in range by cross-correlating the Fourier transform of the data at nearby frequencies and and integrating over to obtain a statistically stable result. The bandwidth is much larger than by assumption (LABEL:eq:model2), and the statistical stability follows essentially from a law of large numbers, because we sum a large number of terms that are uncorrelated.

The second moments of the propagator at nearby frequencies are


where the bar denotes complex conjugate, is the Fourier transform of the Wigner distribution described below, and the scale is defined in terms of the autocorrelations , and (see [9, equation (6.26)]). These formulas follow from the calculations in [13, 4] which assume , and the law of iterated expectation with conditioning at , for . Denoting by the conditional expectation and using

because , we obtain

and (LABEL:eq:st18) follows from [13, 4] and the fact that in the support of the source .

The last term in (LABEL:eq:st18) corresponds to the coherent part of the mode amplitudes and it is negligible in our regime with . This is by (LABEL:eq:st15) and

Recalling the expression (LABEL:eq:st13) of the mode amplitudes in terms of the propagator, we see that (LABEL:eq:st18) states that the amplitudes of different modes are essentially uncorrelated. Therefore, the only second moments that remain large are the mean energies of the modes, which is why we use them in inversion.

3.4 The system of transport equations

The Wigner distribution defines the expectation of the energy of the th mode resolved over a time window of duration similar to the travel time, when the initial excitation is in the th mode. It satisfies the following system of transport equations derived in [20, 13, 4]


with initial condition


where is the Dirac delta distribution. The Fourier transform that appears in (LABEL:eq:st18) is defined by


where is the diagonal matrix


The matrix in (LABEL:eq:st19) models the transfer of energy between the modes, due to scattering. Its off-diagonal entries are defined in (LABEL:eq:st22)


and are non-negative, meaning that there is an outflow of energy from mode to the other modes. The energy lost by this mode is compensated by the gain of energy in the other modes, as stated by


4 Inversion based on energy transport equations

We now use the results summarized above to formulate our inversion approach. We give in section LABEL:sect:formod the forward model which maps the source density to the cross-correlations of the mode amplitudes. These are defined in section LABEL:sect:dataproc and are self-averaging with respect to different realizations of the random waveguide. Therefore, we can relate them to the Wigner distribution. The inversion method is studied in section LABEL:sect:inverse.

4.1 Data processing

The first question that arises is how to relate the incoherent array data to the moments (LABEL:eq:st18) of the propagator which are defined by the Wigner distribution. The answer lies in computing cross-correlations of the data projected on the eigenfunctions , as we now explain.

We denote by the Fourier transform of the measurements and by its projection on the eigenfunction


We are interested in its cross-correlation at lag and its inverse Fourier transform . The latter has the physical interpretation of energy carried by the th mode over the duration of a time window which we model with a bump function of dimensionless argument and order one support


Here has units of frequency, satisfying , so the integrand is compactly supported in the recording window . The scaling by of the argument of , the inverse Fourier transform of (LABEL:eq:inv1), is to be consistent with the travel time of the waves to the array, and the factors in front of the integral are chosen to get an order one


This expression is obtained by taking the inverse Fourier transform of (LABEL:eq:inv3), and the integral over is restricted by the support of to .

We relate below the expectation of to the Wigner distribution, and explain in Appendix LABEL:ap:A under which conditions is self-averaging, meaning that it is approximately equal to its expectation. The self-averaging is due to the rapid frequency decorrelation of over intervals of order , and the bandwidth assumption (LABEL:eq:model2). When we divide the frequency interval in smaller ones of order , we see that in (LABEL:eq:inv2) we are summing a large number of uncorrelated random variables. The self-averaging is basically by the law of large numbers, as long as is large. This happens for large enough arrays, for long recording times that scale as (LABEL:eq:durationT), and for times near the peak of .

The role of the projection (LABEL:eq:inv1) is to isolate in the data the effect of the th mode. We see from (LABEL:eq:model)-(LABEL:eq:model1) that


where we introduced the mode coupling matrix with entries


This coupling is an effect of the aperture of the array. The ideal setup is for an array with full aperture , because is the identity by the orthonormality of the eigenfunctions, and involves only the amplitude of the th mode. However, all the mode amplitudes enter the expression of when the array has partial aperture, and they are weighted by . The coupling matrix is diagonally dominant when the length of the aperture is not much smaller than the waveguide depth . This can be seen for example in the case of an array starting at the top boundary , where


We note in (LABEL:eq:inv4) that by choosing the support of the recording window as in (LABEL:eq:durationT), we can relate to the mode amplitudes in a frequency interval of order . This is important in the calculation of the cross-correlations , where the amplitudes must be evaluated at nearby frequencies. If we had a smaller , the cross-correlations would involve products of the amplitudes at frequency offsets that exceed . Such amplitudes are statistically uncorrelated and there is no benefit in calculating the cross-correlation.

4.2 The forward model

We show in Appendix LABEL:ap:A that


where the diagonal matrix defined in (LABEL:eq:cB) is evaluated at ,


are the Fourier coefficients of the unknown source density, and . Because the cross-correlations are self-averaging we can define the forward map from to the vector using equation (LABEL:eq:A4). We write it as


which is a simplification of (LABEL:eq:A4) based on the assumption that the recording window is well centered and sufficiently long to equal one at the times of interest. We also simplify the notation by dropping the argument of the wavenumbers , their derivatives and . The unknown source density appears in the model as the matrix of absolute values of its Fourier coefficients (LABEL:eq:A2). This is the most that we can expect to recover from the inversion.

5 Inversion

We have the following unknowns: the range , the matrix , and possibly the autocorrelations of the fluctuations. The question is what can be recovered from and how to carry the inversion. The range and some information about the autocorrelation of the fluctuations can be determined from the measurements of the travel times of . This is the easier part of the inversion and we discuss it first. The estimation of is more delicate and requires knowing and the autocorrelation of the fluctuations, so we can calculate the matrix . We discuss it in sections LABEL:sect:I2-LABEL:sect:I6. We illustrate the results with numerical simulations in section LABEL:sect:I3.

5.1 Arrival time analysis

If there were no random scattering effects i.e., no matrix , the