Multiple scattering of ultrasound in weakly inhomogeneous media: Application to human soft tissues

Multiple scattering of ultrasound in weakly inhomogeneous media: Application to human soft tissues

Alexandre Aubry, Arnaud Derode Institut Langevin, ESPCI ParisTech
CNRS UMR 7587, Université Denis Diderot (Paris VII)
10 rue Vauquelin, 75005 Paris, France
July 13, 2019

Waves scattered by a weakly inhomogeneous random medium contain a predominant single scattering contribution as well as a multiple scattering contribution which is usually neglected, especially for imaging purposes. A method based on random matrix theory is proposed to separate the single and multiple scattering contributions. The experimental set up uses an array of sources/receivers placed in front of the medium. The impulse responses between every couple of transducers are measured and form a matrix. Single-scattering contributions are shown to exhibit a deterministic coherence along the antidiagonals of the array response matrix, whatever the distribution of inhomogeneities. This property is taken advantage of to discriminate single from multiple-scattered waves. This allows one to evaluate the absorption losses and the scattering losses separately, by comparing the multiple scattering intensity with a radiative transfer model. Moreover, the relative contribution of multiple scattering in the backscattered wave can be estimated, which serves as a validity test for the Born approximation. Experimental results are presented with ultrasonic waves in the MHz range, on a synthetic sample (agar-gelatine gel) as well as on breast tissues. Interestingly, the multiple scattering contribution is found to be far from negligible in the breast around 4.3 MHz.

©Copyright 2011 Acoustical Society of America. This article may be downloaded for personal use only. Any other use requires prior permission of the author and the Acoustical Society of America. The following article appeared in J. Acoust. Soc. Am. 129, 255-233 (2011), and may be found at

43.60.Gk, 43.20.Fn, 43.80.Ev, 43.35.Bf
preprint: Aubry et al.

I Introduction

Standard imaging techniques, such as ultrasonic echography angelsen (), radar stergiopoulos () or optical coherence tomography Schmitt (), are based on the same principle. One or several source(s) emit(s) a wave into the medium to be imaged. It is reflected by the inhomogeneities of the medium and the backscattered wave is measured by the same or other sensor(s). It contains two contributions:

A single scattering contribution (path in Fig.1): the incident wave undergoes only one scattering event before coming back to the sensor(s). This is the contribution which is used in imaging because there is a direct relation between the arrival time of the echo and the distance between the sensors and the scatterer, ( is the sound velocity). An image of the medium reflectivity can be built from measured signals.

A multiple scattering contribution (path in Fig.1): the wave undergoes several scattering events before reaching the sensor. Multiple scattering occurs when scatterers are strongly diffusive and/or highly concentrated. There is no correspondence between the arrival time and the position of a scatterer. Thus, classical imaging fails in multiple scattering media bly (); bordier (); karamata (); karamata2 ().

Standard imaging techniques rely on the single-scattering assumption (first Born approximation). However, there is no such thing as a purely single scattering medium. A multiple scattering contribution always exists, albeit negligible compared to single scattering. Naturally for imaging purposes, one tries to reduce the influence of multiple scattering, for instance by choosing an appropriate frequency domain where multiple scattering is not too strong aubry4 (). Focused beamforming with an array of transducers, or more generally synthetic aperture techniques stergiopoulos (), are also a way to enhance the single scattering contribution. It should be noted that even though multiple scattering is considered as the enemy of classical imaging techniques, studying it may bring additional information about the scattering structure. Indeed, a wave undergoing multiple scattering can be thought of as a random walker rossum (), with two essential parameters: the elastic mean-free path and the diffusion constant . Measuring these parameters is a way to characterize the microarchitecture of the scattering medium rytov (); rytov2 (). Yet in weakly inhomogeneous media where the first Born approximation is reasonably valid (especially human soft tissues probed by ultrasound in the MHz range), it is a challenge to study multiple scattering parameters because of the predominance of single scattering.

Recently, an original technique has been proposed to separate the single-scattered echo of a target drowned in a predominant multiple scattering background aubry5 (); aubry4 (). The method was based on a matrix approach. It has been successfully applied to target detection and imaging in highly scattering media. In this paper, we are also interested in discriminating single-scattering and multiple-scattering contributions from the total response of an unknown medium, based on matrix properties. However, the present approach is different from earlier works aubry5 (); aubry4 (), both in terms of method and applications. The situation we consider here is exactly the opposite: we want to extract the multiple scattering contribution from predominantly single-scattered waves in a weakly scattering media. Moreover, the method is based on a singular value decomposition applied to the antidiagonals of the array response matrix, and not to the array response matrix itself. The distinction between single and multiple scattering subspaces is then performed using random matrix theory tulino (); sengupta (), as it will be detailed in the next sections.

The interest of this work is twofold. First, once single and multiple scattering contributions are isolated, the proportion of multiple scattering within the wave response of the medium can be evaluated. This figure can be used as an indicator for the validity of the single-scattering (Born) approximation, which is the basis of classical imaging techniques. Note however that the purpose of this study is not to improve imaging of weakly scattering media. The second interest of this work is to provide a new tool for the characterization of weakly scattering media. More precisely, we will show how the multiple scattering contribution, once it is isolated, can be taken advantage of in order to estimate the scattering mean-free path independently from the absorption mean-free path , thus discriminating absorption and scattering losses.

The experimental results we present were obtained with pulsed ultrasonic waves firstly in a synthetic medium (agar-gelatine gel) around 3 MHz then in breast tissues around 4.3 MHz, but the principle of the technique can be applied to all fields of wave physics (e.g seismology, electromagnetism, acoustics etc.) for which the multi-element array technology is available and provides time-resolved measurements of the amplitude and the phase of the wave field.

Figure 1: Experimental setup: a 125-element linear array is placed in front of a random medium at a distance . The whole setup is immersed in a water tank.

Ii Transducers’ array configuration

We use a N-element ultrasonic array (here, ). The array is placed at a distance from the random scattering sample under investigation (see Fig. 1). The first step of the experiment consists in measuring the inter-element matrix. A sinusoidal burst of length , at the central frequency , is emitted from transducer into the scattering medium. Typical values here are s and MHz. The backscattered wave is recorded with the transducers of the same array, which yields a set of impulse responses ( denotes the receiver index). The operation is repeated for the emitting transducers. The responses form the impulse response matrix . Because of reciprocity, and is symmetric. A short-time Fourier analysis of the impulse response matrix is performed. The time signals are truncated into -long overlapping windows : with , anywhere else. The value of is chosen so that signals associated with the same scattering event(s) within the medium arrive in the same time window aubry3 (). Typical values here are s. A Fourier analysis of is achieved by means of a discrete Fourier transform. A response matrix is finally obtained at each time and frequency . The single and multiple scattering contributions can now be discriminated with the help of a matrix manipulation.

Iii The signature of single scattering

When studying the array response matrices , the predominance of single scattering manifests itself by the presence of a long-range deterministic coherence along the antidiagonals of the matrix, whatever the distribution of scatterers aubry4 (); aubry5 (); aubry3 (). As an example, Fig.2 shows the real part of one of the matrices , in the case of a synthetic medium (Agar-gelatine gel) which is enough weakly scattering for the Born approximation to be valid. Even if the inhomogeneities are randomly distributed, obviously exhibits some kind of coherence along its antidiagonals (i.e., for matrix elements such that ). This coherence is a typical signature of single scattering, and it vanishes when multiple scattering dominates. This has been thoroughly explained in aubry4 (); aubry5 (); aubry3 (), we briefly recall the main argument here.

Figure 2: Real part of matrix obtained in a gel (5% gelatine, 3% agar-agar) at time s and frequency 3.05 MHz. The source-sample distance was mm.

Generally, the can be written as the sum of single and multiple scattering contributions:


Under the paraxial approximation, the distance between the origin and an observer located slightly off-axis () is . As a result, the phase shift undergone by a wave travelling from a source with coordinates , scattered by a point and received in the plane of the source at reads :

with the wavenumber. The quadratic phase terms can be factorized since


Consider an ensemble of scatterers randomly distributed. As long as only the first and the last scattering of every scattering path are identical (which is naturally the case, if only single scattering takes place) the coefficients of the array response matrix at time and frequency will be proportional to:


with and the number of scatterers contained in the isochronous volume. depends on the reflectivity of the scatterer. and are random variables so is itself random. Interestingly, applying the factorization of Eq.2 to Eq.3, a deterministic relation arises along the antidiagonals of :


with the array pitch and denotes the distance between two array elements ( and ) on the same antidiagonal. Eq.4 implies that as long as there is only single scattering, there must be a form of coherence, a long-range deterministic relation, between the elements of the array response matrix, whatever the realisation of disorder. On the contrary, when multiple scattering occurs (except for recurrent scattering paths wiersma (), but this contribution is negligible in weakly scattering media), the elements cannot be factorized, and there is no such long-range deterministic coherence aubry4 (); aubry5 (); aubry3 ().

Iv Separation of single and multiple scattering

The key to separate single () and multiple () scattering contributions is the particular coherence of along its antidiagonals. In previous works aubry4 (); aubry5 (), was extracted from by projecting the antidiagonals of along the vector of Eq.4. But the simple form taken by Eq.4 results from a series of assumptions (paraxial approximation, pointlike array elements and scatterers) which do not all apply to our experimental configuration. In order to separate and , the method proposed in this paper is much less restrictive. We do not assume that Eq.4 exactly applies; we only assume that because of single scattering there must be a deterministic coherence between the antidiagonal elements of , but we do not suppose we know its exact form.

Under these conditions, the separation between and will essentially rely on a singular value decomposition (SVD) of the antidiagonals of . This separation is a three-step process:

Rotation of each matrix and construction of two sub-matrices and .

Filtering of matrices . is decomposed as the sum of two matrices: , where and contain respectively the single and multiple scattering signals.

Construction from and of the single and multiple scattering matrices and .

The first and third steps (rotation of data) have already been presented in previous works aubry4 (); aubry5 () and will be briefly recalled in Sec.IV.1 and Sec.IV.3. On the contrary, the second step (SVD of antidiagonals) constitutes the core of the method and differs completely from the previous approach aubry4 (); aubry5 (). The corresponding matrix operations are explained in details in Sec.IV.2.

iv.1 First step

A rotation of matrix data is achieved as depicted in Fig.3. It consists in building two matrices and from matrix :


where . Here and so is an even number.

Figure 3: Example of a matrix of dimension 17. The black points represent the elements of . The antidiagonals of are the columns of matrices and . Circles and squares represent the elements of and , respectively. Once single and multiple scattering contributions are separated, the final matrices and have elements (central square).

The columns of matrices and correspond to the antidiagonals of (see Fig.3). In the next subsection, we will no longer make the difference between matrices and because they are filtered in the same way. They will be called indifferently . is the dimension of . For matrix , ; for matrix , . Because of spatial reciprocity, is symmetric (). Thus, exhibits also a symmetry: each line of its upper part is identical to a line of its lower part. The symmetry axis is shown as a black line in Fig.3 and corresponds to the diagonal of the matrix . So, each column of the matrix contains only independent coefficients, even if its dimension is larger than .

iv.2 Second step

can be written as a sum of two matrices and , which correspond respectively to the single and multiple scattering contributions


Contrary to previous works aubry5 (); aubry4 (), the technique we propose in this paper consists in separating single and multiple scattering by achieving the singular value decomposition (SVD) of the matrix . The SVD decomposes a matrix into two subspaces: a signal subspace (a matrix characterized by an important correlation between its lines and/or columns) and a noise subspace (a random matrix without any correlations between its entries). When the SVD is applied to the matrix , the signal subspace (i.e., the largest singular values) corresponds to (the single scattering contribution characterized by a long-range correlation along its columns) and the noise subspace (i.e., the smallest singular values) corresponds to (the multiple scattering contribution).

The SVD of matrix is given by


and are square unitary matrices of dimension . Their respective columns and correspond to the singular vectors associated to the singular value . is a square diagonal matrix of dimension , containing the real positive singular values in a decreasing order (). Actually, has only non-zero singular values since it contains only independent lines, hence Eq.8 becomes:


The issue is to determine which rank of singular value separates the signal subspace (single scattering) from the noise subspace (multiple scattering). If Eq.4 were strictly true, the single scattering contribution would be of rank 1 and only the first singular space associated to the first singular value would correspond to the signal subspace. But when assumptions leading to Eq.4 do not strictly hold, is no longer of rank 1 and several singular spaces (associated to the largest singular values) are needed to fully describe the signal subspace. This happens for instance when scatterers are not pointlike or when the paraxial approximation does not hold. We have to define a threshold to discriminate the signal and noise subspaces, with the help of random matrix theory (RMT) tulino (); sengupta ().

By convention and for the sake of simplicity, the singular values are first normalized by their quadratic mean


For a random matrix of dimension (with ), whose entries are complex random variables, independently and identically distributed, the probability density function of the normalized singular values is given by sengupta ()


for and 0 otherwise, with


For random matrices of large dimensions, the singular value spectrum has a bounded support. In our case, is a square matrix of dimension . Yet, as it contains only independent lines, it is equivalent to a rectangular matrix. If it were truly random (which is expected to be the case of the multiple scattering contribution) its largest singular value should not exceed . This value is obtained from Eq.12, with and .

It should be noted that rigorously and are not large enough for the asymptotic law (Eq.11) to apply. Actually the first singular value obeys a complicated law, known as the Tracy-Widom distribution johnstone (), which is of unbounded support. The probability for to be larger than can be computed: it is found to be 0.08 for and . The presence of correlations between matrix entries also induce a deviation from Eq.11, as we will see later. For the sake of simplicity, we admit for now that within an acceptable probability of error, the singular values are upper bounded by .

According to Eq.7, is the sum of a matrix of rank (associated to single scattering) and a matrix of rank M (associated to multiple scattering). Sengupta and Mitra sengupta () have shown that the smallest singular values (linked to the noise subspace) exhibit the same distribution as singular values of a random matrix whose size is . Let denote the upper bound of the singular values distribution in the case of a random matrix of dimension ; we have:


From this property, one can propose a way to separate the signal and noise subspaces of . We first consider the first singular value upon normalization (Eq.10). If is larger than (Eq.13), it means that the first singular space is associated with the signal subspace. Then, we iterate the process and consider the second singular value. The are once again renormalized, considering only the singular values from :


The threshold value to consider this time is the one obtained for a random matrix of size , i.e (Eq.13). If , the second singular space is also linked to the single scattering contribution, and we iterate once more the process until rank for which . Finally, we obtain a threshold rank which allows to separate the signal () and noise () subspaces,


Ideally, should be devoid of multiple scattering. This is not strictly true, because multiple scattering signals are not strictly orthogonal to the single scattering subspace. Let and be the power of single and multiple scattering signals. In Appendix A, we show that the typical amplitude of the remaining multiple scattering contribution in is (). If we neglect this residual term, we have separated single and multiple scattering contributions: and . The whole separation process is summarized in Fig.4.

Figure 4: Principle of the separation between the single and multiple scattering contributions.

Note that a multiple scattering rate can be directly measured from the singular values of . The sum of the square of all the singular values corresponds to the total intensity backscattered by the medium towards the transducers’ array . Hence, a multiple scattering rate can be estimated from the singular values of A:


Until now, for simplicity we have implicitly assumed that along the antidiagonals of (or the columns of ) the matrix elements are completely decorrelated. However, experimentally short-range correlations may exist between elements, mostly because of mechanical coupling between neighboring transducers and of the coherence length of the diffuse wave-field aubry5 (); aubry3 (). Correlations between matrix elements can be taken into account theoretically sengupta (). Consequently, the actual probability density function is more complicated than the simple result of Eq.11, which modifies the upper bound aubry5 (); aubry3 (). In practice, has to be computed numerically, based on a acceptable probability of error (typically 1%). The details are given in Appendix B.

This technique of separation is based on the fact that the first singular value exceeds the value , otherwise there is no separation between single and multiple scattering and the whole signal is considered to be associated with multiple scattering. So, this approach is not well suited for strongly diffusive media, i.e random media for which the multiple scattering contribution is predominant aubry5 (); aubry4 ().

iv.3 Third step

The third step is the reverse of the first one. From and , two matrices and , of dimension , are built(see Fig.3) with a change of coordinates, back to the original system:

if is an integer,

if is not an integer,

contains the single scattering contribution (plus a residual multiple scattering contribution) and contains the multiple scattering contribution.

V Characterization of a weakly scattering medium

The experimental set up has already been described in Sec.II and is shown in Fig.1. The experiment takes place in a water tank. The ultrasonic array has elements. The emitted signal is a sinusoidal burst of length s at the central frequency (3 MHz). The sampling frequency is 20 MHz. Each array element is 0.39 mm in size and the array pitch is 0.417mm. The source-sample distance is mm. The first random medium of interest is a gel composed of 5% of gelatine and 3% of agar. In this kind of medium and frequency range, the single scattering contribution is by far predominant aubry (). The thickness of the scattering slab is 100 mm. Once the inter-element matrix is measured, the short-time Fourier analysis described in Sec.II yields a set of response matrices . Then, the separation of single and multiple scattering is achieved as described in Sec.IV.

Fig.5 shows a typical experimental result, taken at time s and frequency MHz. Note that the separation rank between the signal and noise subspaces is here , which confirms that Eq.4 doesnot strictly hold.. exhibits the deterministic coherence along the antidiagonals (Fig.5(a)) which is characteristic of single scattering. Obviously, is very close to the raw matrix (Fig.2), since single scattering is predominant. As to , it displays a random feature as expected for the multiple scattering contribution (Fig.5(b)). However, one cannot conclude that it originates in multiple scattered waves: it could also correspond to experimental noise.

Figure 5: Separation of single and multiple scattering contributions at time s and frequency MHz. (a) Real part of . (b) Real part of .

In order to establish the multiple-scattering origin of , we calculated the mean backscattered intensity as a function of the source-receiver distance and the arrival time


The symbol denotes an average over the quantities in the subscript, here frequency and all source/receiver couples separated by the same distance . In Fig.6, the spatial dependence of is compared with the total intensity, at a given time . Whereas the total intensity shows no preferred direction, exhibits a typical signature of multiple scattering: the coherent backscattering peak clearly arises around .

Figure 6: The multiple scattering intensity (black dots) and the total intensity (grey circles) are plotted versus , at time s. The intensity profiles have been renormalized with their maximum.

This phenomenon has been widely observed and studied in wave physics (optics Ishimaru1 (); wolf (); albada2 (); akkermans2 (); akkermans3 (); wiersma (); labeyrie1 (); labeyrie2 (), acoustics bayer (); tourin2 (); yamamoto (); rosny (); rosny3 (), seismology larose (); tiggelen (); margerin ()). It manifests itself as an enhancement, by a factor 2, in the backscattered intensity at the vicinity of the source (i.e ). Its physical origin lies in the constructive wave interference between reciprocal paths that have been scattered at least twice; it can only appear when multiple scattering occurs and the reciprocity symmetry is preserved. The intensity profile shown in Fig.6 is a spectacular evidence of multiple scattering and shows the efficiency of our technique for extracting the multiple-scattering waves among a predominant single scattering contribution.

Interestingly, though it is weak, the multiple scattering contribution can be taken advantage of in order to characterize the medium and determine separately the scattering losses and the absorption losses. When a wave propagates through a random medium, it loses progressively its coherence: after traveling over a distance , only a fraction of the initial energy still propagates in coherence with the initial wave. The parameter , called the extinction mean-free path, characterizes the extinction length of the coherent part of the wave. It comes from two distinct phenomena (scattering and intrinsic absorption of the medium) which are associated to two characteristic lengths: the elastic mean-free path and the absorption mean-free path , such that


Experimentally, can be determined by measurements of the ensemble-averaged field transmitted through a scattering layer page (); zhang (); scales1 (); scales2 (); derode4 (). However, this kind of experiments does not allow to distinguish from .

Figure 7: (a) Single-scattered intensity versus time. Experimental measurements (black circles) are fitted with the theoretical curve (continuous black line) considering an extinction length mm. (b) Multiply-scattered intensity versus time. Experimental measurements (white squares) are compared with thoretical results for (continuous black lines) and (dashed black lines) for different values of the mean free path, while keeping mm. All intensities have been normalized by the maximum of the single scattered intensity over time.

We focus on the single and multiple scattering intensities obtained at the source: and . They are plotted in Fig.7. Note that the intensity of the multiple scattering contribution is less than 1% of the single scattering contribution. Once and have been measured, we can fit both experimental curves with and as independent adjustable parameters. To that end, we need a theoretical model describing the spatial and temporal evolution of the mean intensity inside the random medium. In the literature, the mean intensity is often assumed to obey the diffusion equation akkermans (). The diffusion approximation is simple, but only valid in the long-time limit. Since we deal with a weakly scattering medium, the elastic mean free path is expected to be very large compared to the scattering path lengths (). Thus, the diffusion approximation does not apply to our problem. Instead, we used the radiative transfer equation (RTE) rossum (). Paasschens paasschens () proposed an exact solution of the RTE in time-domain and real space for an infinite 2D random medium (see Appendix C), Eq.C21). Based on this theoretical work, we have computed the exact expression of the single and double scattering intensities, and , as well as an approximate expression of the multiple scattering intensity , considering the medium as semi-infinite and two-dimensional.

The choice of a 2D model is justified as follows. Experimentally, the transducers are 10 mm in height, which is much larger than the average wavelength (0.5 mm). Moreover, a vertical cylindrical acoustic lens ensures that the emitted beam remains collimated in the plane. Similarly, in reception only waves propagating in the plane are recorded by the transducers. Thus the single scattering problem is clearly 2D. As to multiple scattering, for the same reason the only paths that can generate a signal on the receiving transducers are those whose first and last scatterer are in the plane. The gel sample being weakly scattering, is mostly dominated by double scattering (see Fig.7(b)). Thus even though the wave propagation in the gel sample is 3D, we have used the 2D solution for the RTE.

The detailed calculations of , and are shown in Appendix C. It appears that the single scattering intensity exhibits a temporal evolution which only depends on the extinction length . In the case of the gel studied here, the best fit of the experimental results yields mm (Fig.7(a)). Once is known, and can be determined by fitting with theory (Fig.7(b)) with only one adjustable parameter since . The scattering gel is found to be much more absorbing than scattering : mm, while mm. Fig.7(b) also displays the theoretical evolution of the double scattering contribution . For mm, and are nearly identical, which shows that the double scattering contribution clearly dominates the multiple scattering intensity in the gel sample. As the theoretical expression derived in Appendix C is exact for , the measured values of and are reliable.

In this example, the medium was a weakly scattering gel, with a ratio less than 1%. Yet the separation of single and multiple scattering contributions can also be achieved in real scattering media for which is closer to unity, as we present in the next section.

Vi Application to human soft tissues

Figure 8: (a) Experimental setup used for the investigation of multiple scattering in breast tissues. (b) Echographic image of the breast. The grey scale is in dB. (c) Multiple scattering rate as a function of time.

The same kind of experiment has been performed in a biological medium for which ultrasound is often used: the breast. The experimental set up is depicted in Fig.8(a). We use a -element ultrasonic array () with a 4.3 MHz central frequency and a 3.5-5MHz bandwidth; the array pitch is 0.33 mm. The emitted signal is a 0.7-s sinusoidal burst at MHz. The sampling frequency is 50 MHz. The experimental procedure is the same as in Sec.II. The separation of single and multiple scattering contributions can be performed as in Sec.IV, but an adjustment has to be made. Indeed, in our experimental configuration, the array of transducers is placed in the near-field of the scattering medium (). Consequently, the entries of matrix are not identically distributed: the variance (i.e the mean intensity ) of decreases significantly with the distance between the source and the receiver, as shown by Fig.9. This implies a different variance for each line of matrix , hence modifying its theoretical distribution of singular values. The upper bound can be computed numerically taking into account this non-uniform distribution of matrix elements (See Appendix D). However this is only possible at the first iteration (). Indeed, whereas the variance of can be estimated, the variance distribution of the subspaces of is unknown a priori. Unless we do the (strong) approximation that this variance is uniform, in which case Eq.13 could be applied, we cannot follow the procedure described in Sec.IV. Here, since clearly has a non-uniform variance, by precaution we choose a different strategy: the same upper bound is considered at each iteration of the single/multiple scattering separation process (see Sec.IV). Since , this precaution tends to overestimate the threshold and decrease the probability of error.

Fig.9 compares the multiple scattering, single scattering and total intensity profiles, at a given time . Contrary to the previous experiment in the gel sample (see Fig.6), the spatial intensity profiles, and , are not flat due to the near-field configuration of the experiment. Once again, exhibits a coherent backscattering peak on top of a flat incoherent intensity with an enhancement factor close to 2. Interestingly, this indicates that is not a noise contribution, but does originate from multiple wave scattering in the breast tissue, even though we operate at a frequency (4.3MHz) for which human soft tissues are usually treated as single-scattering.

Figure 9: The multiple scattering intensity (black dots), the single scattering intensity (grey circles) and the total intensity (black squares) are plotted versus , at time s. The intensity profiles have been renormalized with the maximum of the total intensity.

An ultrasound image of the breast has also been obtained with the same array, using 63-element subapertures (Fig.8(b)). As usual in ultrasound imaging, focusing in emission and reception is achieved by applying a set of 63 time delays to the signals transmitted/received by the array. The time delays are computed in order to focus at the desired region of interest, centred at coordinates and , assuming that the velocity of sound in soft tissues is known. In the case of breast, as for most soft biological tissues, it is close to that of water ( m/s). Here, 2666 focal planes ( 8 mm to 48 mm) and 63 values of ( to with the array pitch) have been used. At each depth , a new set of time-delays is calculated; this is more demanding that classical ultrasound imaging techniques, which generally use the same set of time-delays as long as is within the depth of field. A line of the resulting image represents, in gray level, the amplitude of the total echographic signal at the focal time , once focused beamforming has been applied to the 63 received signals. The resulting image displays the reflectivity of the medium under investigation. Ultrasound images of human tissues usually reveal the interfaces of inner organs, and often exhibit a speckled appearance due do random scattering by sub-wavelength inhomogeneities (cells, fibers, tissues etc.) Burckhardt (). Here the scanned area is particularly echogene between 30-40 s, corresponding to the depth range 22.5-30 mm. The typical ultrasound image of a human organ (here, the breast) is representative of the amplitude of backscatter at a given time which hopefully (under the single scattering assumption) corresponds to a given depth, but it does not allow us to distinguish single and multiple scattering contributions.

However, once the matrix has been recorded, not only can we build an echographic image but we can also isolate the multiple scattering contribution, and estimate the multiple scattering rate (Eq.16). has been averaged over the whole frequency spectrum and is displayed on Fig.8(c) as a function of time. Fig.8(c) complements the information brought by the echographic image. A relevant observation is that multiple scattering becomes predominant from s, that is to say beyond a depth of 34.5 mm. It means that the single scattering assumption, upon which the imaging process is based, is incorrect. It does not mean however that the image is totally wrong; a rate of means that half of the intensity received by one individual array element comes from multiple scattering. In classical array imaging, each line of the picture is constructed by focused beamforming in emission and reception. This procedure reduces the importance of multiple scattering in the final image because the single scattered contributions coming from a target in the focal zone add up coherently whereas the contributions from multiple scattering can be expected to be uncorrelated. With and assuming uncorrelated array elements, the multiple scattering rate becomes after beamforming. This is probably an underestimation since multiple scattering signals cannot be fully uncorrelated aubry5 (); aubry3 (). As a result, the proportion of multiple scattering in the final image around s is of the order of a few percents, which is still weak. Also note that at larger times, the technique we presented here would fail: after s the rate of multiple scattering becomes too large for the SVD to extract the single scattering contribution (), at least for some frequencies of the spectrum. As the results are averaged over the whole frequency spectrum (3.5-5 MHz), the multiple scattering rates that are presented here are still meaningful, until s. Beyond that time, the method we presented here would be inadequate to separate single and multiple scattering.

Vii Conclusion

The approach we developed here can separate single- and multiple-scattered waves in randomly heterogeneous media. It requires an array of transmitters/receivers and takes advantage of the persistence of a deterministic coherence of single scattering signals along the antidiagonals of the inter-element matrix. Once a singular value decomposition is applied, the single scattering contribution (signal subspace) is separated from the multiple scattering contribution (noise subspace) by using a criterion based on random matrix theory. Unlike previous works aubry5 (); aubry4 () this technique is particularly well-suited for weakly scattering media, for which single scattering dominates. In such media, the technique we presented here is not intended to enhance the quality of the ultrasound image, but rather to complement it, in two ways. Firstly, the experimental results indicate that this approach can be applied for characterization purposes: the separation of single and multiple scattering provides a way to measure the scattering and absorption mean-free paths independently. This idea was tested on a synthetic gel. Secondly, the technique was also applied in vivo to the case of breast imaging with ultrasonic waves around 4.3 MHz. The occurrence of multiple scattering has been established and its contribution to the backscattered wave-field is shown to be far from negligible. By measuring the relative amount of multiple scattering, the method serves as an experimental test for the first Born approximation (single scattering), which is usually made in such tissues.

The authors would like to thank Pr Mickael Tanter for fruitful discussions and Patricia Daenens for her technical help, as well as the groupe de recherches IMCODE of CNRS (GDR 2253).

Appendix A

Single and multiple scattering signals are not strictly orthogonal. Thus, can be expressed as the sum of the single scattering contribution plus a residual multiple scattering noise :


The aim of this appendix is to estimate the remaining part of multiple scattering which corrupts the signal subspace . In other words, we want to determine the mean intensity of coefficients ().

To that aim, we apply the first order perturbation theory cohen () to the autocorrelation matrix . Actually, can be decomposed as a sum of an unperturbed autocorrelation matrix (linked with single scattering) and a perturbation matrix due to multiple scattering :


As we deal with weakly scattering media, the single scattering contribution is predominant: . Thus, can be seen as a low perturbation of the autocorrelation matrix. Moreover, we can neglect the second order term in Eq.21. can be simplified into


Let us first focus on the reference state, i.e with no multiple scattering. In this case, we have


where denotes the singular value of . and are the singular vectors associated to .

The perturbation theory can now be applied to . At first order, this theory states that only the eigenvalues of are perturbed by multiple scattering. The singular vectors remain identical to those obtained in the unperturbed case (i.e. without multiple scattering), hence and . The first eigenvalues of can be written as cohen ():


is the unperturbed eigenvalue (that we would obtain without multiple scattering) and the term corresponds to the perturbation due to multiple scattering.

Using Eqs.24 & 22, one can try to express :


is a random complex variable with zero mean. Let us calculate its variance:


We assume that the matrix is random, hence the ensemble average of yields


where is the identity matrix and is the mean power of multiple scattering signals: . The combination of Eqs.26 & 27 yields:


Let us calculate the variance of (Eq.25):


By injecting Eq.28 in Eq.29, we obtain


Using the approximation , the last equation becomes:


where .

Using Eq.15 of the manuscript and Eq.23 and the fact that and , the matrix can be expressed as a function of :


Now, we can estimate the variance of the coefficients :


Because the singular vectors and are normalized, and . Using Eq.31, we finally obtain:


The mean power of the residual multiple scattering noise is . This residual noise is negligible compared to single scattering signals since .

Appendix B

This appendix describes the numerical process we use to determine the upper bound of the singular values distribution, when short-range correlations between matrix entries exist.

The first step is to estimate these correlations. This is done by measuring the correlation coefficients and between respectively the lines and columns of :


where the symbol denotes an average over the variables in the subscript, i.e the source/receiver pairs . The integer represents the distance between sources and receivers, in units of the array pitch .

The second step consists in building two correlation matrices (of size ) and (of size ) from the measured correlation coefficients and . The coefficients and of these correlation matrices are given by:


Once these correlation matrices are built, the next step consists in generating a random matrix which displays the same short-range correlations as matrix .

To that aim, a random matrix of dimension is first generated numerically. Its coefficients are totally independent and identically distributed. Then, a matrix is built from , such that


The matrix exhibits the same correlation properties as . The SVD of is then achieved and its first singular value is renormalized according to


By repeating this operation over 500 realizations, we can build a cumulative histogram of the first singular value . The distribution function of the first singular value, , can be estimated. The upper bound is then deduced from . An acceptable probability of error is first set (typically, ), hence the upper bound :


Appendix C

In this appendix, we derive an expression for the single and multiple-scattering contributions to the average backscattered intensity, based on radiative transfer theory. Radiative transfer theory describes the spatial and temporal dependence of the radiance (or specific intensity) in a random medium. Radiance is defined as energy flow, propagating in the direction , per unit normal area per unit solid angle per unit time. It follows a transport (Boltzmann) equation :


with the source term and the intensity, defined as the angular average of the radiance


Classically, the transport equation can be derived from the Bethe-Salpether equation frisch (); margerin2 (), neglecting all interference (coherent) effects. Here we adapt the theoretical developments of Paasschens paasschens (), which were derived for an infinite random medium, to our experimental configuration. The problem is solved in two dimensions.

The typical configuration is depicted in Fig.1 of the manuscript. The random medium is assumed to be semi-infinite, with a plane interface. It is characterized by an elastic mean-free path and an absorption length . The scattering events in the random medium are assumed to be isotropic, hence there is no difference between the elastic mean-free path and the transport mean-free path .

In our case, the medium is not statistically invariant under translation since it is semi-infinite: depends both on the observer () and the source (). The radiative transfer equation is modified into:


where an isotropic point source is now considered: . accounts for the semi-infinite nature size of the random medium:


This problem can be solved by considering separately the contributions to intensity from scattering events


The partial intensities satisfy


Let us first investigate the ballistic intensity . Following Paaschens paasschens (), the differential operators on the left hand side of Eq.47 can be integrated, to yield


The Dirac distribution describes the propagation of the energy pulse at the finite speed . accounts for the geometrical spreading of the ballistic radiance in two dimensions. is the ballistic path length within the scattering medium:

  • If , then .

  • If , then .

  • In other configurations (see Fig.10), the ballistic radiance propagates partly in free space, partly in the random medium: .

The angular average of Eq.49 yields the ballistic intensity :