Methods and systems for non-invasively characterizing a heterogeneous medium using ultrasound

Abstract

In a first aspect, the present description relates to a system for non-invasively characterizing a heterogeneous medium using ultrasound, comprising at least one first array (10) of transducers configured to generate a series of incident ultrasound waves in a region of said heterogeneous medium and record the ultrasound waves which are backscattered by said region as a function of time, as well as a computing unit (42) which is coupled to said at least one first array of transducers and configured to: record an experimental reflection matrix defined between an emission basis at the input and the basis of the transducers at the output; determine, from said experimental reflection matrix, at least one first wideband reflection matrix defined in a focused base at the input and output; determine a first distortion matrix defined between said focused basis and an observation basis, said first distortion matrix resulting, directly or after a change of basis, from the term-by-term matrix product of said first wideband reflection matrix defined between said focused basis and an aberration correction basis, with the phase-conjugated matrix of a reference reflection matrix defined for a model medium, between the same bases; determine, from said first distortion matrix, at least one mapping of a physical parameter of said heterogeneous medium.

Claims

1. A method for the non-invasive ultrasound characterization of a heterogeneous medium, the method comprising: a step of generating, by means of at least a first array of transducers, a series of incident ultrasound waves in a region of said heterogeneous medium; a step of recording an experimental reflection matrix defined between an emission basis at input and the transducer basis at output; a step of determining, on the basis of said experimental reflection matrix, at least a first wideband reflection matrix, defined in a focused basis at input and at output, in a first spectral band; a step of determining a first distortion matrix defined in said first spectral band, between said focused basis and an observation basis, said first distortion matrix resulting, directly or after change of basis, from the term-by-term matrix product of said first wideband reflection matrix, defined between said focused basis and an aberration correction basis, with the phase conjugate matrix of a reference reflection matrix, defined for a model medium, between the same bases; a step of determining, on the basis of said first distortion matrix, at least one mapping of a physical parameter of said heterogeneous medium.

2. The method as claimed in claim 1, wherein the determination of said first wideband reflection matrix comprises: a change of basis at input and at output, from said experimental reflection matrix expressed in the frequency domain to form a reflection matrix, defined in the frequency domain, in a focused basis at input and at output; a coherent summation on said first given spectral band of said reflection matrix defined in a focused basis at input and at output to form said first wideband reflection matrix defined in a focused basis at input and at output.

3. The method as claimed in claim 1, wherein the determination of said first distortion matrix comprises: a change of basis, at input or at output, of said first wideband reflection matrix defined in a focused basis at input and at output, to obtain said first wideband reflection matrix determined between the focused basis and the aberration correction basis; the term-by-term matrix product of said first wideband reflection matrix defined between said focused basis and the aberration correction basis, with the phase conjugate matrix of the reference reflection matrix defined between the same bases.

4. The method as claimed in claim 1, wherein: said observation basis is the focused basis; said first distortion matrix is defined in a focused basis at input and at output, and is obtained either directly from said first wideband reflection matrix defined in a focused basis at input and at output, or by change of basis from the first distortion matrix defined between the focused basis and the aberration correction basis.

5. The method as claimed in claim 4, wherein the determination of at least one mapping of a physical parameter of said heterogeneous medium comprises determining, from said first distortion matrix defined in a focused basis at input and at output, and for at least one first isoplanatic domain identified, at least one first point spread function (PSF).

6. The method as claimed in claim 1, wherein the determination of at least one mapping of a physical parameter of said heterogeneous medium comprises: determining the invariants in said focused basis of said first distortion matrix, in order to identify at least one first isoplanatic domain in said focused basis; determining, for each first isoplanatic domain identified, a mapping of at least one first aberration law in said aberration correction basis.

7. The method as claimed in claim 6, wherein the determination of the invariants of said first distortion matrix comprises a singular value decomposition of at least one of the matrices of the group of matrices comprising: said first distortion matrix, said first distortion matrix normalized, a normalized correlation matrix of said first distortion matrix.

8. The method as claimed in claim 6, wherein the determination of at least one mapping of a physical parameter of said heterogeneous medium further comprises determining in said observation basis a first wideband reflection matrix corrected by said one or more first aberration laws.

9. The method as claimed in claim 8, further comprising: determining a corrected distortion matrix, defined between said focused basis and said aberration correction basis by the term-by-term matrix product of said corrected first wideband reflection matrix with the phase conjugate matrix of said reference reflection matrix.

10. The method as claimed in claim 9, further comprising: determining the invariants in said focused basis of said corrected distortion matrix, in order to identify at least one second isoplanatic domain in said focal plane; determining, for each second isoplanatic domain identified, a mapping of a second aberration law in said correction basis.

11. The method as claimed in claim 1, further comprising: a step of determining, on the basis of said experimental reflection matrix, at least one second wideband reflection matrix, defined in a focused basis at input and at output, in a second spectral band different from said first spectral band; a step of determining at least one second distortion matrix defined in said second spectral band between said focused basis and said observation basis, said second distortion matrix resulting, directly or after change of basis, from the term-by-term matrix product of said second wideband reflection matrix, defined between said focused basis and said aberration correction basis, with the phase conjugate matrix of a reference reflection matrix, defined for a model medium, between the same bases.

12. The method as claimed in claim 11, further comprising: determining the invariants in said focus basis of said at least one second distortion matrix, in order to identify at least one first isoplanatic domain in said focal plane; determining, for each isoplanatic domain identified, a mapping of at least one first aberration law in said correction basis; determining, in said observation basis, at least one second wideband reflection matrix corrected by said one or more first aberration laws; summing said corrected first wideband reflection matrix and said corrected at least one second wideband reflection matrix to form an ultra-wideband corrected reflection matrix.

13. The method as claimed in claim 1, wherein the determination of at least one mapping of a physical parameter of said heterogeneous medium comprises determining the speed of sound in the insonified region.

14. The method as claimed in claim 1, wherein the determination of at least one mapping of a physical parameter of said heterogeneous medium comprises determining a multiple scattering rate in the insonified region.

15. The method as claimed in claim 1, further comprising the identification and/or elimination of the specular reflection and/or multiple reflection components on the basis of said first distortion matrix, which is expressed at input and at output in the plane wave basis.

16. A system for the non-invasive ultrasound characterization of a heterogeneous medium comprising: at least one first array of transducers that is configured to generate a series of incident ultrasound waves in a region of said heterogeneous medium and to record the ultrasound waves backscattered by said region as a function of time; a computing unit, coupled to said at least one first array of transducers, and configured to: record an experimental reflection matrix defined between an emission basis at input and the transducer basis at output; determine, on the basis of said experimental reflection matrix, at least one first wideband reflection matrix, defined in a focused basis at input and at output, in a first given spectral band; determine a first distortion matrix defined in said first spectral band, between said focused basis and an observation basis, said first distortion matrix resulting, directly or after change of basis, from the term-by-term matrix product of said first wideband reflection matrix, defined between said focused basis and an aberration correction basis, with the phase conjugate matrix of a reference reflection matrix, defined for a model medium, between the same bases; determine, on the basis of said first distortion matrix, at least one mapping of a physical parameter of said heterogeneous medium.

Description

BRIEF DESCRIPTION OF THE FIGURES

(1) Other advantages and features of the technique presented above will become apparent from reading the detailed description below, provided with reference to the figures in which:

(2) FIGS. 1A-1C (already described) illustrate known emission/reception mechanisms for ultrasound imaging and quantification;

(3) FIG. 2 (already described) illustrates a technique for correcting aberrations in ultrasound imaging, according to the prior art, based on adaptive focusing;

(4) FIGS. 3A-3B (already described) schematically illustrate the concept of isoplanatic domains;

(5) FIG. 4 illustrates an exemplary ultrasound characterization system for implementing the ultrasound characterization methods according to the present description;

(6) FIG. 5 illustrates one step of the method according to the present description, aiming to determine, according to one exemplary embodiment, a first wideband focused reflection matrix on the basis of an experimental reflection matrix;

(7) FIG. 6 illustrates the basis change matrix operations for the projection of reflection matrices recorded in the focused base;

(8) FIGS. 7A, 7B illustrate one step of the method according to the present description, aiming to determine, according to one exemplary embodiment, a first distortion matrix on the basis of the wideband focused reflection matrix;

(9) FIG. 8 illustrates ultrasound images of a phantom without aberrations and with aberrations and the corresponding focused distortion matrices, determined at a given depth and corresponding to a line of the image;

(10) FIGS. 9A and 9B illustrate, respectively, the distribution of the singular values of the distortion matrix corresponding to the phantom with aberrations shown in FIG. 8 (FIG. 9A) and the phase of the associated first eigenvector in the Fourier plane (FIG. 9B);

(11) FIG. 10 illustrates an exemplary use of the distortion matrix (case of an aberrating layer of water on a phantom giving rise to a single isoplanatic domain) and more precisely illustrates ultrasound images without/with correction, the corresponding focused distortion matrices determined at a given depth and corresponding to a line of the image and a representation of the spread function with and without correction averaged over the insonified region and resulting from the focused distortion matrix;

(12) FIG. 11 illustrates an exemplary use of the distortion matrix (on a human breast phantom generating a plurality of isoplanatic domains) and more precisely illustrates a conventional ultrasound image and ultrasound images after correction of aberrations using the first vector of the distortion matrix and the second vector of the distortion matrix;

(13) FIGS. 12A-12C illustrate an exemplary use of the distortion matrix for the filtering of multiple reflections and more precisely, FIG. 12A, a diagram describing the experimental setup (Plexiglas plate placed between the ultrasound probe and the phantom), FIG. 12B, an ultrasound image for a homogeneous speed model (c=1800 m/s), before and after removal of specular components and multiple reflections from the image, examples of a focused distortion matrix for one line of the ultrasound image, before and after filtering, the same distortion matrices projected into the plane wave base, and FIG. 12C, a diagram explaining the signature of the specular and multiple reflection components in the plane wave base;

(14) FIG. 13 illustrates an exemplary use of the distortion matrix for a quantitative measurement of the speed of sound and illustrates more precisely an ultrasound image of the phantom studied, an example of a focused distortion matrix for a line of the image and a homogeneous speed model (c=1540 m/s), the corresponding distortion matrix between the focal plane at input and the Fourier plane, the width of the spread function normalized by the wavelength as a function of the speed of the waves c, the normalized first singular value of the distortion matrix as a function of this same speed;

(15) FIG. 14 illustrates an exemplary use of the distortion matrix for establishing a mapping of the speed of sound in the medium and more precisely illustrates a conventional ultrasound image of the phantom through the Plexiglas plate for a homogeneous speed model (c=1540 m/s), and the same image after filtering multiple reflections and optimizing the speed model for each line of the image, the profile of speed of sound integrated over depth;

(16) FIG. 15 illustrates an exemplary use of the distortion matrix for the local quantification of multiple scattering and more specifically shows a standard phantom ultrasound image, a diagram explaining the occurrence of multiple scattering in the distortion matrix, the focused distortion matrix for a line of the ultrasound image, an example of a spread function deduced from the preceding matrix showing an incoherent background linked to the multiple scattering component, the multiple scattering rate in each pixel of the ultrasound image;

(17) In the various embodiments which will be described with reference to the figures, similar or identical elements bear the same references.

DETAILED DESCRIPTION

(18) In the detailed description which follows, only some embodiments are described in detail in order to ensure clarity of the description, but these examples are not intended to limit the general scope of the principles that emerge from the present description.

(19) The various embodiments and aspects described in the present description may be combined or simplified in multiple ways. In particular, the steps of the various methods may be repeated, reversed, or performed in parallel, unless otherwise specified.

(20) FIG. 4 illustrates an exemplary ultrasound characterization system 40 for implementing the methods for the ultrasound characterization of a heterogeneous medium 20, according to the present description. The system 40 comprises at least one first array 10 of transducers 11, for example a linear or two-dimensional array; the transducers are, for example, ultrasound piezoelectric transducers which may take the conventional form of a rigid bar brought into contact with the medium 20. The array of transducers forms, for example, part of a probe device 41; the array of transducers is connected to a computing unit 42, which may itself be connected to a display device 43; the computing unit transmits and records electrical signals to/from each of the transducers 11. The ultrasound transducers then transform these electrical signals into ultrasound waves and vice versa. The computing unit 42 is configured for the implementation of calculating or processing steps, in particular for the implementation of steps of methods according to the present description.

(21) In FIG. 4, as in the remainder of the description, reference is made to an array of transducers for emission and reception, it being understood that, in a more general case, a plurality of arrays of transducers could be used simultaneously. They may both emit and receive, or else only emit for some and only receive for others.

(22) When, in the present description, reference is made to calculating or processing steps for the implementation in particular of method steps, it is understood that each calculating or processing step may be implemented by software, hardware, firmware, microcode or any appropriate combination of these technologies. When software is used, each calculating or processing step may be implemented by computer program instructions or software code. These instructions may be stored in or transmitted to a storage medium that is readable by a computer (or computing unit) and/or be executed by a computer (or computing unit) in order to implement these calculating or processing steps.

(23) The various experimental results shown in the present description (FIGS. 5 to 15) were taken from two ultrasound phantoms, i.e. synthetic samples reproducing the mechanical properties of biological tissues (c=1540 m/s) and generating an ultrasound speckle characteristic of same. The first phantom makes it possible to quantify the performance of the imaging system in terms of resolution and contrast. Specifically, it comprises near-pointlike targets located at various depths and lateral positions and a heterogeneous echogenic inclusion 5 mm in diameter (FIGS. 5 to 10 and 12 to 15). The second phantom reproduces the shape and properties of breast tissue and exhibits bright echoes mimicking the presence of microcalcifications within same (FIG. 11).

(24) The present description describes methods and systems for the non-invasive ultrasound characterization of a heterogeneous sample. These methods and systems are based on the determination of at least one first matrix called a “distortion matrix” in the remainder of the description. FIGS. 5 to 7 illustrate, according to some examples, steps of methods for determining the distortion matrix.

(25) As illustrated in FIG. 5, step 51, the non-invasive ultrasound characterization method comprises a step of generating, by means of an array 10 of transducers 11 (FIG. 4), a series of incident ultrasound waves in a region of said heterogeneous medium and a step of recording an experimental reflection matrix R.sub.ui(t) defined between an emission basis i at input and the transducer basis u at output; the emission basis i at input is for example the transducer basis u, the plane wave basis θ or the focused basis r, as described by means of FIGS. 1A-1C. In the example of FIG. 5, it is for example an experimental reflection matrix R.sub.uθ(t) recorded between the plane wave basis at input and the transducer basis at output.

(26) The method also comprises a step of determining, on the basis of said experimental reflection matrix R.sub.ui(t), at least one first wideband focused reflection matrix R.sub.rr(Δω.sub.1), defined in a focused basis r at input and at output, in a first spectral band Δω.sub.1. This step prior to the construction of the “distortion matrix” is illustrated according to one exemplary embodiment, by means of FIGS. 5 and 6.

(27) In the example of FIG. 5, a discrete Fourier transform operation 52 makes it possible to express the experimental reflection matrix in the frequency domain R.sub.uθ(ω). A step 53 of “wideband focusing” of the experimental reflection matrix comprises the projection of the reflection matrix R.sub.uθ(ω) into the focal basis at input and at output in order to obtain a focused reflection matrix R.sub.rr(ω), then a coherent summation on the spectral band Δω.sub.1. A wideband focused reflection matrix is then obtained, an example of which is illustrated (image 54) at a given depth z.

(28) FIG. 6 illustrates the projection of the reflection matrix R.sub.uθ(ω) into the focal basis r at input and at output. By way of example, image 63 represents the phase of the matrix R.sub.uθ(ω) at the frequency f=3.48 MHz, acquired in the first phantom, comprising near-pointlike targets and a heterogeneous echogenic inclusion of 5 mm in diameter, of which an ultrasound image is shown in FIG. 5 (image 56). This amounts to performing beamforming on emission (61) and on reception (62) in order to coherently sum the echoes coming from the same scatterer. Mathematically, these operations may be described using propagation matrices describing the outward and return paths of the waves in the medium. The transition from the incident plane wave basis to the focused basis is thus modeled by the matrix G′.sub.0(ω)=[G′.sub.0(r.sub.in, θ.sub.in, ω)], the elements of which are given by
G′.sub.0(r.sub.in,θ.sub.in,ω)=exp[jk(z.sub.in cos θ.sub.in+x.sub.in sin θ.sub.in)]  (1)
with r.sub.in=(x.sub.in, z.sub.in) the position vector of the focus points in the medium, θ.sub.in the angle of incidence of the emitted wave and k=ω/c the wavenumber. The path of the echoes reflected from the medium to the array of transducers may be described by the propagation matrix G.sub.0(ω)=[G.sub.0(r.sub.out, u.sub.out, ω)], the elements of which are given by

(29) $\begin{array}{cc}{G}_0\left(r_{\mathrm{out}},u_{\mathrm{out}},\omega \right)=\frac{\exp \left[\mathrm{jk}.Math.u_{\mathrm{out}}-r_{\mathrm{out}}.Math.\right]}{\sqrt{k.Math.u_{\mathrm{out}}-r_{\mathrm{out}}.Math.}}\frac{z_{\mathrm{out}}}{.Math.u_{\mathrm{out}}-r_{\mathrm{out}}.Math.}& \left(2\right)\end{array}$
with r.sub.out=z.sub.out) the position vector of the focus points in the medium and u.sub.out the position vector of the transducers. The last equation is a matrix translation of the |Rayleigh-Sommerfeld integral. Its first term is the 2D Green's function of the wave equation, or rather here its asymptotic form. The right-hand term is the direction cosine between the normal of the transducer array and the radial vector (u.sub.out−r.sub.out).

(30) The operations of beamforming on emission and on reception may be interpreted as backpropagation of the echoes to the scatterers associated therewith. In the frequency domain, this translates as a phase conjugation of the propagation matrices. The focused reflection matrix R.sub.rr may thus be deduced from the experimental reflection matrix via the following matrix product:
R.sub.rr(ω)=G*.sub.0(ω)×R.sub.uθ(ω)×G′.sub.0.sup.†(ω)  (3)
where the symbol × represents the matrix product and the exponent † denotes the conjugate transpose matrix. Of course, the projection of the reflection matrix into the focal bases at input and at output could also be carried out on the basis of an experimental reflection matrix expressed at input in another base, for example the transducer base. For this, it would be enough to use, in equation (3), the propagation matrix G.sub.0 instead of the propagation matrix G′.sub.0.

(31) The next step consists in summing the reflection matrices obtained on the spectral band Δω.sub.1 so as to obtain a wideband focused reflection matrix
R.sub.rr(τ)=Σ.sub.∫∈Δω.sub.1R.sub.rr(ω)e.sup.jωτ  (4)
τ is the echo time associated with the virtual array of transducers in the focused basis. Each coefficient R(r.sub.out, r.sub.in, τ) of the wideband focused reflection matrix corresponds to the analytical signal which would be recorded by the virtual transducer at r.sub.out immediately after the emission of a pulse from a virtual transducer at r.sub.in (diagram 55, FIG. 5). The diagonal of the matrix R.sub.rr(τ=0) in ballistic time τ=0 corresponds to the confocal signal with focusing on emission and on reception on the same point. It therefore directly gives a line of the ultrasound image (image 56) if the focused basis is limited to a focal plane located at a distance z from the real array of transducers. If, however, the focused basis contains all of the focus points of the volume to be imaged, the diagonal of R.sub.rr(τ=0) directly gives the entirety of the ultrasound image. In this latter case, the reflection matrix is a block matrix, only the diagonal blocks of which are non-zero. These correspond to the focused reflection matrices at each depth z of the medium. In the remainder of the description, in order to simplify the notation, the wideband reflection matrix will be called simply R.sub.rr, there will no longer be any mention of the time τ because it will always be set at ballistic time τ=0. However, for some applications, in particular for the correction of high-order aberrations, the study of multiple scattering or the measurement of the speed of sound, the reflection and distortion matrices may also be examined at times τ≠0.

(32) FIGS. 7A, 7B illustrate one step of the method according to the present description, aiming to determine, according to one exemplary embodiment, a first distortion matrix D on the basis of the wideband focused reflection matrix R.sub.rr.

(33) In order to estimate the aberrated component of the wavefronts, to start, an optimal basis for correcting the aberrations is chosen. For example, in the case of a thin aberrator located in proximity to the array of transducers, the transducer basis appears to be the most suitable. In the case chosen in this illustrative example of the characterization method, the heterogeneous medium is a translation invariant multilayer medium and the plane wave basis appears to be more suitable. This is for example the case with echography of the liver, where the sound waves have to pass through a succession of layers of muscle and adipose tissue. This example is considered hereinafter but the present description is in no way restricted to correction of aberrations from the Fourier base.

(34) Image 71 (FIG. 7A) shows an example of a wideband focused reflection matrix R.sub.rr recorded in the presence of an aberration generated by the insertion of a layer of water (c=1480 m/s) 10 mm thick between the array of transducers and the phantom (c=1540 m/s). First, a change of basis of the reflection matrix R.sub.rr at output (or at input) into the plane wave basis makes it possible to form the wideband reflection matrix R.sub.θr, such that
R.sub.θr=.sup.tG′.sub.0(ω.sub.c)×R.sub.rr  (5)
where the matrix G′.sub.0 at the central frequency ω.sub.c of the spectral band Δv.sub.1 is considered. Physically, each row of the matrix R.sub.θr thus obtained corresponds to the wavefront backscattered in the far field for each insonified point r.sub.in. Each wavefront has two components: a geometric component, linked to the diffraction of the wave between the transducer plane and the focal plane, and a distorted component, linked to the aberrations undergone by the wavefront going out and on return. In order to separate the two components, the matrix R.sub.θr is compared with a reference matrix that would be obtained in a model medium, for example a homogeneous medium without aberration (having only the geometric component). This reference matrix in the case of a homogeneous medium without aberrations is nothing other than the propagation matrix in a homogeneous medium G′.sub.0. A new matrix is then created (FIG. 7A, step 72), called the distortion matrix in the present description, from the Hadamard product or term-by-term matrix product between said first wideband reflection matrix R.sub.θr determined between the focused basis and the aberration correction base, and the phase conjugate matrix of the reference reflection matrix, G′*.sub.0(ω.sub.c), defined between the same bases at the central frequency ω.sub.c:
D.sub.θr=R.sub.θr∘G′.sub.0.sup.†(ω.sub.c)  (6)
where the symbol o represents the term-by-term matrix product (Hadamard matrix product) and the exponent † denotes the conjugate transpose matrix. To recall, the elements of the phase conjugate matrix G′.sub.0.sup.† of the matrix .sup.tG′.sub.0 are complex numbers having the same modulus as those of .sup.tG′.sub.0 but being of opposite phase. In terms of matrix coefficients, this equation is written as:
D(θ.sub.out,r.sub.in)=R(θ.sub.out,r.sub.in)×G′*.sub.0(r.sub.in,θ.sub.out,ω.sub.c)  (7)

(35) The operation performed above is illustrated by FIG. 7B. It consists in subtracting the deterministic geometric component from the measured wavefront and thus isolating its distorted component. Each column of the matrix D.sub.θr, contains this aberrated component of the wavefront seen from the far field for each focus point r.sub.in. Thus, in practice, to isolate the distorted component of the reflected wavefront, there is subtracted from the phase of the reflected field the phase of the field that would be obtained ideally in the absence of aberrations, for a model medium.

(36) In the case of an inhomogeneous model medium, the expression of the reference reflection matrix G′.sub.0 becomes more complicated than that given in equation (1). For example, in the case of a multilayer medium, the matrix G′.sub.0 may be calculated analytically by a matrix product between the different transmission matrices associated with the propagation of the wave in each of the layers. For a more complicated model medium, typically one that is not translation invariant (e.g. curved interface), the matrix G′.sub.0 may be determined by means of a numerical simulation or a semi-analytical calculation. Note that if the coefficients of the reference reflection matrix have an inhomogeneous modulus, the reference matrix may be normalized by imposing a constant modulus for each of these coefficients but retaining their phase.

(37) Image 73 in FIG. 7A shows the phase of the distortion matrix D.sub.θr deduced from the reflection matrix, the modulus of which is shown in image 71. The latter corresponds to the line of the ultrasound image 83 (FIG. 8) produced on the phantom in the presence of an aberrating layer of water. Diagram 74 illustrates the isolated aberrated wavefronts for each point of the focal plane and the way in which they are stored in the distortion matrix.

(38) The distortion matrix may also be studied in the focal plane. A “focal plane” distortion matrix D.sub.rr may be obtained on the basis of the distortion matrix expressed in the observation plane D.sub.θr via the following change of base:
D.sub.rr=Ĝ′*.sub.0(ω.sub.c)D.sub.θr  (8)
where Ĝ′.sub.0 is a normalized propagation matrix, the elements of which are given by

(39) $\begin{array}{cc}{\hat{G}}_0^{\prime }\left(r_{\mathrm{out}}^{\prime },{\theta }_{\mathrm{out}},\omega \right)=\frac{{\hat{G}}_0^{\prime }\left(r_{\mathrm{out}}^{\prime },{\theta }_{\mathrm{out}},\omega \right)}{{\hat{G}}_0^{\prime }\left(r_{\mathrm{out}}^{\prime }=0,{\theta }_{\mathrm{out}},\omega \right)}=\exp \left[{\mathrm{jkx}}_{\mathrm{out}}^{\prime }\sin {\theta }_{\mathrm{out}}\right]& \left(9\right)\end{array}$
Equation (8) is therefore rewritten in terms of matrix coefficients as follows:
D(r′.sub.out,r.sub.in)=Σ.sub.θ.sub.outD(θ.sub.out,r.sub.in)exp[−jk.sub.cx′.sub.out sin θ.sub.out]=R(r.sub.out−r.sub.in,r.sub.in)  (10)
where k.sub.c is the wavenumber at the central frequency ω.sub.c. In the present case where the aberration correction basis is the Fourier base, the focused distortion matrix D.sub.rr is simply deduced via spatial Fourier transform of D.sub.θr at output. However, it should be kept in mind that the relationship between the two matrices will be different [equation (8)] if the correction plane is no longer the Fourier plane.

(40) The columns of D.sub.rr, the modulus of which is shown in image 76, correspond to the field reflected in the image plane recentered on the focus point r.sub.in. The matrix D.sub.rr therefore gives the change in the reflection spread function at each point of the ultrasound image, the spread function being the spatial impulse response of the imaging system, that is to say the focal spot recentered on the focus point. As will be described below, the spread function makes it possible to quantify and characterize the aberrations and the multiple scattering caused by the sample upstream of the focal plane.

(41) By way of illustration, FIG. 8 illustrates an ultrasound image 81 of the phantom obtained without aberrations, by means of an experimental diagram recalled in diagram 80. Image 82 illustrates the focused distortion matrix determined at a given depth and corresponding to a line of image 81; image 83 is an ultrasound image of the phantom obtained in the presence of aberrations (aberrating layer of water 22, diagram 85) and image 84 illustrates the focused distortion matrix determined at a given depth and corresponding to a line of image 83. In the experiment without aberrations, all of the energy reflected by the medium appears on the central row of the distortion matrix (r.sub.in=r.sub.out): the width δ of the spread function is close to the limit imposed by diffraction (δ˜λz/D, with D the size of the array of transducers) which explains the excellent quality of the ultrasound image in terms of resolution and contrast. However, in the presence of the aberrating layer of water, the energy reflected spreads outside of the central row. The degradation of the spatial impulse response implies a significantly less well resolved ultrasound image. This example shows how the “focal plane” distortion matrix makes it possible to quantify the quality of focusing in the medium, independently of a qualitative analysis of the ultrasound image.

(42) To go further in the analysis of the distortion matrix D.sub.θr, it is advantageously possible to consider, on the one hand, the case of a reflection by the sample that is mainly specular and, on the other hand, the case of a reflection that is mainly random scattering. In practice, it is often an intermediate case which combines random scattering reflection and specular reflection. In both cases, the distortion matrix D.sub.θr exhibits correlations between its columns. These correlations correspond to the repetition of the same distortion pattern for the wavefronts coming from the same isoplanatic domain. For a sample causing a random scattering reflection, for example of speckle type (random distribution of unresolved scatterers), it will be possible to study the normalized correlation matrix of the distortion matrix D.sub.θr in order to extract these same correlations more efficiently. Whichever the case, as will be described in more detail below, searching for the invariants of the matrix D.sub.θr, for example by means of a singular value decomposition (SVD), makes it possible to extract the complex transmittance of the aberrator for each point of the focal plane and thus optimally correct the aberrations in each isoplanatic domain contained in the insonified region.

(43) Assuming a sample causing random scattering reflection, examples of the use of the distortion matrix to determine a mapping of one or more aberration laws in the observation base, in particular to establish a mapping of the local reflectivity of the insonified region, are described below.

(44) For this, the invariants of the distortion matrix are sought, in other words the aberration laws that are spatially invariant over the isoplanatic domains of the insonified region. Various methods are known to those skilled in the art for searching for the invariants of such a matrix, such as for example singular value decomposition (or “SVD”) or principal component analysis (“PCA”).

(45) Singular value decomposition is a powerful tool for extracting the correlations between the rows or columns of a matrix. Mathematically, the SVD of the matrix Dθ.sub.r, of size N×M, is written as follows:
D.sub.θr=U×Σ×V.sup.†  (11)
U and V are unit matrices of size N×M, the columns U.sub.i and V.sub.i of which correspond to the output and input eigenvectors. Each output eigenvector U.sub.i is defined in the Fourier plane identified by the angle θ.sub.out. Each input eigenvector V.sub.i is therefore defined in the focal plane identified by the vector r.sub.in. Σ is a matrix of size N×M, only the diagonal elements of which are non-zero:
Σ=diag(σ.sub.1,σ.sub.2, . . . ,σ.sub.N.sub.2  (12)
The diagonal elements of the matrix Σ are the singular values σ.sub.i of the matrix D.sub.θr which are real, positive and ranked in descending order:
σ.sub.1>σ.sub.2> . . . >σ.sub.N.sub.2  (13)
The coefficients D(θ.sub.out, r.sub.in) of the matrix D.sub.θr are therefore written as the following sum:
D(θ.sub.out,r.sub.in)=Σ.sub.i=1.sup.Lσ.sub.iU.sub.i(θ.sub.out)V*.sub.i(r.sub.in)  (14)
with L=min (M, N).

(46) SVD primarily decomposes a matrix into two subspaces: a signal subspace (a matrix characterized by substantial correlations between its rows and/or its columns) and a noise subspace (a random matrix without correlation between its rows and columns). The signal subspace is associated with the largest singular values while the noise subspace is associated with the smallest singular values. On the one hand, the SVD of D will therefore make it possible to filter the noise subspace which contains both the experimental noise and the incoherent contribution of the reflected field caused by multiple scattering events. On the other hand, each singular state of the signal subspace will make it possible to extract, according to the output eigenvector U.sub.i, the distortion undergone by the wave in the aberration correction basis for each isoplanatic domain of the image. In particular, in the case of a single isoplanatic domain, it is possible to show that the first output eigenvector U.sub.1 directly gives the complex transmittance A of the aberrator between the correction plane and the focal plane:
U.sub.1(θ.sub.out)=A(θ.sub.out)  (15)

(47) The following figures illustrate various applications of the distortion matrix, whether defined in the focal bases or between the focal basis and the aberration correction base.

(48) FIGS. 9A and 9B show the result of the SVD applied to the distortion matrix D.sub.θr shown in image 84 of FIG. 8. The spectrum of the singular values σ.sub.i of D.sub.θr (FIG. 9A) shows a dominance of the first singular value σ.sub.1. This indicates the presence of a single isoplanatic domain. The associated first eigenvector U.sub.1 (FIG. 9B) shows a parabolic phase, characteristic of a defocus defect explained by the presence of the layer of water in which the speed of sound (1490 m/s) is substantially different from that of the phantom (1540 m/s).

(49) To correct the wideband reflection matrix R.sub.rr (71, FIG. 7A), to start, in this example, the matrix R.sub.rr is projected into the aberration correction plane (here the Fourier plane):
R.sub.θθ=.sup.tG′.sub.0(ω.sub.c)×R.sub.rr×G′.sub.0  (16)
The aberrations are then corrected by applying the phase conjugate of U.sub.1 at input and output of the matrix R.sub.θθ:
R′.sub.θθ=exp(j×arg{U*.sub.1})∘R.sub.θθ∘exp(j×arg{U.sub.1.sup.†})  (17)
where R′.sub.θθ is the corrected reflection matrix and the function arg{ . . . } denotes the phase of the vector between curly brackets. Relationship (18) translates for the matrix coefficients as follows:
R′(θ.sub.out,θ.sub.in)=U*.sub.1(θ.sub.out)R(θ.sub.out,θ.sub.in)U*.sub.1(θ.sub.in)  (18)
Physically, the phase conjugation operation consists in re-emitting a wavefront modulated by a phase opposite that of the measured distortion and in applying a similar correction at output. This operation then makes it possible to compensate perfectly for the phase distortions accumulated by the wave on its outward and return paths. A corrected matrix R′.sub.rr may then be deduced from R′.sub.θθ via change of basis at input and output from the pupillary plane to the focal plane. A corrected confocal image I′ may be deduced from the diagonal (r.sub.in=r.sub.out) of the matrix R′.sub.rr:
I′(r.sub.in)=R′(r.sub.in,r.sub.in)  (19)
I′(r.sub.in) is then a reliable estimator of the reflectivity ρ(r.sub.in) of the sample.

(50) FIG. 10 illustrates the benefits of this correction by comparing the aberrated ultrasound image 101 and the ultrasound image 103 obtained according to the correction described by equation (17). A qualitative gain between these two images is observed in particular by virtue of the echogenic targets which appear much more contrasted and better resolved after correction of the aberrations.

(51) In echography, however, most of the time there are no echogenic targets for judging the resolution of the image. To quantify the latter, the focused distortion matrix may be used. That deduced from the starting reflection matrix, D.sub.rr, is presented in image 102 at a particular depth. That deduced from the corrected reflection matrix, D′.sub.rr, is presented in image 104 at the same depth. In the second case, the energy is much more concentrated around the central row of the focused distortion matrix than in the first case.

(52) On the basis of the focused distortion matrix, it is also possible to determine the spread function. This corresponds to the spatial Fourier transform of the aberration law measured in the aberration correction plane. It is spatially invariant over each isoplanatic domain. It may also be obtained from the singular value decomposition of the focused distortion matrix. If the matrix G′.sub.0 is unitary (which is the case if the Fourier and focused bases are conjugate with one another), then the matrix D.sub.rr has a singular value decomposition similar to D.sub.θr. The output eigenvectors, which will be called W.sub.i, are directly linked to the eigenvectors U.sub.i via a simple Fourier transform:
W.sub.i=Ĝ′*.sub.0(ω.sub.c)×U.sub.i  (20)

(53) Diagram 105 of FIG. 10 compares the modulus of the spread functions thus obtained before and after correction. The gain in terms of resolution is obvious and the diffraction limit is reached after correction of the aberrations.

(54) In the experiment performed, the measured aberrations are not dependent on the frequency since the speeds of sound in the phantom and in water vary little over the frequency band studied. In practice, this is not always the case and a frequency analysis of the distortion matrix may prove to be judicious in order to determine the space-frequency adaptive filter (or a space-time equivalent) that optimally corrects the aberrations. Specifically, the distortion matrix D may be studied on different spectral bands Δω.sub.i centered around distinct central frequencies ω.sub.c,i. An aberration law U.sub.1(ω.sub.c,i) may be extracted on each of these spectral bands and a corrected reflection matrix R′.sub.rr(ω.sub.c,i) may be deduced therefrom [equation (17)]. A coherent sum of the reflection matrices corrected independently on each of the spectral bands Δω.sub.i allows access to the ultra-wideband corrected reflection matrix R′.sub.rr:
R′.sub.rr=Σ.sub.ω.sub.c,iR′.sub.rr(ω.sub.c,i)  (21)
This frequency correction of aberrations is equivalent in beamforming to the convolution of the focused signals emitted and measured by the spatio-temporal filter U.sub.1(τ), where
U.sub.1(τ)=Σ.sub.ω.sub.c,iU.sub.1(ω.sub.c,i)e.sup.jω.sup.c,i.sup.τ  (22)

(55) Up to this point, the case of a translation invariant aberrating layer has been used, thus generating only one isoplanatic domain. A second exemplary use of the distortion matrix is now described in the case of a volumic aberrator that is not translation invariant, thus generating a plurality of isoplanatic domains.

(56) In this case, an iterative approach (post-processing) to the correction of aberrations is favored. Step #0 of the process consists in correcting the wideband reflection matrix R.sub.θr at output via the conjugate of the phase of the first eigenvector U.sub.1 of the distortion matrix D.sub.ur:
R.sub.θr.sup.(1)=exp(j×arg{U*.sub.1})∘R.sub.θr  (23)
where R.sub.θr.sup.(1) is the resulting corrected matrix. This initial step makes it possible to perform an overall correction for aberrations over the entirety of the image.

(57) Step #1 consists in recalculating a new distortion matrix D.sub.ru.sup.(1) deduced from R.sub.ru.sup.(1) defined this time between the Fourier plane at input and the focal plane at output. The random nature of the object means that it makes more sense to study the correlation matrix of D.sub.ru.sup.(1) in the Fourier plane, namely B.sup.(1)=.sup.tD.sub.rθ.sup.(1)D.sub.rθ.sup.(1)*. A singular value decomposition of the phase of this matrix B.sup.(1), exp[jarg{B.sup.(1)}], gives a new set of eigenvectors U.sub.i.sup.(1) that are associated with each isoplanatic domain of the image. The corresponding reflection matrix may therefore be corrected this time at input:
R.sub.rθ.sup.(2)=R.sub.rθ.sup.(1)∘exp(−jarg{U.sub.i.sup.(1)})  (24)

(58) The following steps consist in reproducing the same process by alternately correcting for residual aberrations at output (even iterations) and input (odd iterations). However, in each step, the correction is still performed with the first eigenvector U.sub.1.sup.(n) of the phase of the correlation matrix, exp[jarg{B.sup.(n)}], of the distortion matrix in the pupillary plane. Specifically, the choice of isoplanatic domain is made in step #1. Depending on whether the correction is at output or at input, the correlation matrix in the pupillary plane B.sup.(n) is given by: B.sup.(n)=D.sub.θr.sup.(n)D.sub.θr.sup.(n)† for even n (output) and B.sup.(n)=.sup.tD.sub.rθ.sup.(1)D.sub.rθ.sup.(1)* for odd n (input). In each step of the process, an image of the field of vision may be deduced from the diagonal of the reflection matrix R.sub.rr.sup.(n) expressed in the focal plane. In practice, a few iterations are sufficient to obtain an optimal correction for the selected isoplanatic domain. An image of the entirety of the field of vision may then be obtained by combining the corrections determined for each isoplanatic domain.

(59) FIG. 11 illustrates the benefits of the method according to the present description by means of an experiment on a breast phantom generating a plurality of isoplanatic domains and comprising microcalcifications at a depth of between 20 and 25 mm FIG. 11 shows a conventional ultrasound image 111, an ultrasound image 112 after correction of aberrations using the first eigenvector U.sub.1.sup.(1) of the normalized correlation matrix exp [jarg{B.sup.(1)}], and a closeup 113 on part of the preceding ultrasound image. Also shown is an ultrasound image 116 after correction of aberrations using the second vector U.sub.2.sup.(1) and, in image 117, a closeup on part of the preceding ultrasound image. It can be seen from these images that the first eigenvector U.sub.1.sup.(1) is associated with the isoplanatic domain containing the microcalcifications at a depth of between 20 and 25 mm Specifically, the ultrasound image (112) is more contrasted and better resolved at this depth. Conversely, the second eigenvector U.sub.2.sup.(1) optimally corrects image (117) at shallower depths (between 10 and 15 mm).

(60) Knowledge of the aberration laws at each point of the sample may be put to good use not only for imaging but also for physically focusing the ultrasound waves at any point of the sample. For example, on the experimental device of FIG. 2, it may be a matter of applying, to the array of transducers, a convolution of the initial focused signal via a spatio-temporal filter deduced from the eigenvectors of the distortion matrix. In the presence of a single isoplanatic area, it will be a matter of considering the Fourier transform of the first eigenvector U.sub.1 [equation (22)]. Note that if the aberration law is frequency stable, the spatio-temporal filter will simply consist in applying a delay law Δt that is directly proportional to the phase of the eigenvector U.sub.1 in the transducer base
Δt(u)=−arg{U.sub.1(u)}/ω.sub.c  (25)
In the presence of a plurality of isoplanatic areas, the aberration laws are determined by means of an iterative process (see previous paragraph). It will be a matter of bringing together, according to a vector W, all of the aberration corrections obtained at input of the sample:
W=exp(−j[arg{U.sub.i.sup.(1)}+Σ.sub.n=1 arg{U.sub.1.sup.(2n+1)}])  (26)

(61) In specular reflection mode and in the case of an ultrasound image contained in a single isoplanatic domain, it is possible to show that the first output eigenvector U.sub.1 of the distortion matrix gives the distortion caused by the aberrator cumulatively going out and on return:
U.sub.1(θ.sub.out)=A(θ.sub.out)A(θ.sub.in)δ(θ.sub.out+θ.sub.in)  (27)
with A(θ.sub.out) and A(θ.sub.in) the distortions undergone by the wavefront going out and on return and projected into the Fourier plane. The Dirac distribution δ in the preceding equation reflects the fact that, in specular reflection mode, a plane wave with an angle of incidence θ.sub.in is reflected as a plane wave at an angle θ.sub.out=−θ.sub.in at output. Additionally, it is possible to show that the input eigenvector V.sub.1 gives, for its part, direct access to the reflectivity ρ of the object:
V.sub.1(r.sub.in)=ρ(r.sub.in)  (28)
The reflection matrix may be corrected for aberrations by applying thereto, only once at input or output, the phase conjugate of the eigenvector U.sub.1:
R′.sub.ur=exp(−j×arg{U.sub.1})∘R.sub.ur  (29)
In the case of an ultrasound image containing a number P>1 of isoplanatic domains, the distortion matrix is of rank P. The input eigenvectors V.sub.i decompose the object in the focal plane over different isoplanatic domains and are associated with different aberration laws U.sub.i in the Fourier plane. The linear combination of the moduli of the eigenvectors V.sub.i weighted by the associated eigenvalues σ.sub.i.sup.2 finally gives access to an image of the field of vision corrected for the aberrations caused upstream
|ρ(r.sub.in)|.sup.2=Σ.sub.i=1.sup.Pσ.sub.i.sup.2|V.sub.i(r.sub.in)|.sup.2  (30)

(62) However, a reflector is never perfectly specular and its random scattering component will therefore not be corrected optimally by equation (29). Furthermore, a specular reflector may cause multiple reflections between its interfaces which may blur the ultrasound image. It is therefore necessary to be able to separate the specular and random scattering components of the field reflected by the medium. According to one or more exemplary embodiments, the characterization method therefore also comprises the identification and/or elimination of the specular component of the reflected field and/or of the multiple reflections arising between the array of transducers and the various interfaces or layers of the heterogeneous medium. To this end, the distortion matrix may be projected into the aberration correction basis both at input and at output. A “Fourier plane” distortion matrix D.sub.θθ may be obtained on the basis of the distortion matrix expressed in the observation plane D.sub.θr via the following change of base:
D.sub.θθ=D.sub.θr.sup.tG′.sub.0(ω.sub.c)  (31)

(63) The previous equation is therefore rewritten in terms of matrix coefficients as follows:
D(θ.sub.out,θ′.sub.in)=Σ.sub.r.sub.inD(θ.sub.out,r.sub.in)exp[jk(z.sub.in cos θ′.sub.in+x.sub.in sin θ′.sub.in)]=R(θ.sub.out,θ.sub.in+θ.sub.out)  (32)

(64) As illustrated in FIG. 13C, the specular and multiply reflected components appear for precise reflected and incident angle pairs, θ.sub.in and θ.sub.out, such that
θ.sub.in+θ.sub.out=Δθ  (33)
where Δθ corresponds to the angle between the specular reflector and the array of transducers. Δθ. It is for example zero if the specular reflector is parallel to the array of transducers. The specular and multiply reflected components may therefore be easily filtered and just the random scattering component (speckle) of the reflected field may be retained in order to precisely determine the aberrations caused by the Plexiglas plate and correct them. This discrimination of the random scattering component then allows direct access to the aberration laws to be applied at input and output in order to correct the reflection matrix and obtain an optimal image of the medium [equation (17)], which is not possible if the specular component predominates [equation (29)]. In addition, the filtering of the distortion matrix in the Fourier plane makes it possible to eliminate the multiple reflections which contaminate, for example, the ultrasound images at shallow depth (multiple echoes between the probe and the surface layers of the human body).

(65) FIGS. 12A-12C illustrate the use of the distortion matrix for the filtering of multiple reflections and its subsequent application to aberration correction and mapping the speed of sound. FIG. 12A shows the diagram of the experimental setup, namely a Plexiglas plate 121 placed between the ultrasound probe 41 and the phantom 20. The corresponding ultrasound image 122 is constructed using a homogeneous speed model (c=1800 m/s). This image is strongly degraded by the multiple reflections and the aberrations caused by the Plexiglas plate. FIG. 12B shows various forms of the distortion matrix for a depth of the ultrasound image (z=25 mm). Image 124 shows the modulus of the focused distortion matrix with a spread function that is particularly degraded by the presence of the Plexiglas plate. Image 125 shows the distortion matrix in the Fourier plane. A strong specular component for θ.sub.in+θ.sub.out=0 is clearly apparent. Image 126 shows the same distortion matrix after application of a specular filter which consists in eliminating all of the components of the field for which |θ.sub.in+θ.sub.out|<1/60 rad. Image 127 shows the focused distortion matrix after filtering of the specular component. Comparison with image 126 (before filtering) shows a distortion matrix, which while still strongly aberrated, has a random appearance typical of the random scattering mode (speckle). The central row of the distortion matrices thus filtered gives the ultrasound image 123 of FIG. 12A. Comparison with the raw ultrasound image 122 shows that the applied filter has succeeded in successfully eliminating the multiple reflections caused by the Plexiglas plate. Even though image 123 is still subject to aberrations, each of the targets of the phantom is now visible, which was not the case in the initial image 122.

(66) According to one or more exemplary embodiments, the determination of at least one mapping of a physical parameter of said heterogeneous medium comprises determining the speed of sound in the insonified region. This involves studying, for example, the change in the first singular value of the distortion matrix as a function of the speed model used to model the reference matrix. The optimal speed model is that which maximizes the first singular value of the distortion matrix. An alternative is to study the focused distortion matrix. Specifically, this makes it possible to study the change in the spread function as a function of the speed model used for the calculation of the reference matrix. The optimal speed model is that which minimizes the width of the wavelength-normalized spread function.

(67) Thus, FIG. 13 illustrates exemplary use of the distortion matrix for quantitative measurement of the speed of sound. FIG. 13 shows the ultrasound image 131 of the phantom already shown in image 81 (FIG. 8) in the absence of aberrations. Image 132 shows an example of a focused distortion matrix for a line of the image and a homogeneous speed model (c=1540 m/s) and image 133 shows the corresponding distortion matrix between the focal plane at input and the Fourier plane. Diagram 134 illustrates the width δ of the spread function normalized by the wavelength λ as a function of the speed of the waves c, and diagram 135 shows the normalized first singular value, {tilde over (σ)}.sub.1=σ.sub.1/Σ.sub.i=1.sup.Nσ.sub.i, of the distortion matrix as a function of this same speed.

(68) The ratio δ/λ in diagram 134 clearly shows at least around c=1539 m/s while the speed given by the phantom maker is 1540 m/s. An adequate speed model implies optimal focusing of the ultrasound waves on emission and on reception and a spread function which is correspondingly better resolved. The focused distortion matrix therefore provides a relevant criterion for quantitative measurement of the speed of sound in the phantom. Likewise, the first singular value {tilde over (σ)}.sub.1 of the distortion matrix in diagram 135 shows a maximum around a speed c=1542 m/s. Here again, the measured speed is close to that given by the maker. Specifically, the closer to reality the speed model used to calculate the distortion matrix is, the closer the distortion matrix is to a singular matrix of rank 1. This experiment shows how the distortion matrix offers various criteria for quantitative measurement of the speed of sound. Here, for proof of concept, it is limited to the case of a homogeneous medium but, of course, it is potentially possible to use the distortion matrix and its parameters δ/λ or {tilde over (σ)}.sub.1 as information feedback in order to construct a local mapping of the speed of sound in the medium. A first step for this will be presented in FIG. 14.

(69) Beyond the reflectivity of the medium or knowledge of the aberration laws at each point of a sample, the distortion matrix makes it possible to perform 3D tomography of the optical index of the sample studied. From the “focal plane” distortion matrix D.sub.rr, there is access to the spread function in the sample as well as to the isoplanatic domains. Instead of directly correcting the images, the idea is to change the reference medium on which is based the calculation of the distortion matrix from the reflection matrix. This leads to an iterative approach to the distortion matrix in which the reference medium is made to change in the direction of a more complex model (e.g. multilayer medium) so as to decrease the spatial extent of the spread function and increase the size of the isoplanatic domains. By going deeper and deeper into the sample, it is possible to gradually reconstruct a three-dimensional map of the speed of sound, while obtaining an image of the reflectivity of the medium that is all the more faithful the closer the reference medium is to reality.

(70) FIG. 14 illustrates an exemplary use of the distortion matrix for establishing a mapping of the speed of sound in the medium. FIG. 14A shows a conventional ultrasound image 141 of the phantom through the Plexiglas plate for a homogeneous speed model (c=1540 m/s). This ultrasound image is of poor quality, due to the aberrations and multiple reflections caused by the Plexiglas plate. By contrast, the reflectivity of the phantom is successfully revealed in ultrasound image 142. The latter was constructed first by filtering the multiple reflections in the image (FIG. 12) and then by optimizing the speed model for each line of the image (FIG. 13). These operations make it possible to obtain an optimum resolution only limited by diffraction. Only the low spatial frequencies of the targets are poorly recombined due to the filtering of the specular component. Diagram 143 shows the profile obtained for the speed of sound integrated over depth. This speed increases sharply on entering the Plexiglas layer (z=14 mm) before gradually decreasing as it leaves (z=17 mm). An algorithm then makes it possible to solve the inverse problem consisting in deducing from this integrated speed profile the exact profile of the speed of sound as a function of depth. The speed of sound measured in the Plexiglas is 2700 m/s and that in the phantom is, as expected, 1540 m/s.

(71) According to one or more exemplary embodiments, the determination of at least one mapping of a physical parameter of said heterogeneous medium comprises determining a multiple scattering rate in the insonified region. This involves studying, for example, the focused distortion matrix for which aberrations will have been corrected beforehand. The ratio of the average of the incoherent intensity (multiple scattering) obtained beyond the focal spot to the intensity at the focus point (central row of the distortion matrix) gives access to a local multiple scattering rate.

(72) FIG. 15 illustrates the advantage of the distortion matrix for the local quantification of multiple scattering. FIG. 15 shows a standard phantom ultrasound image 151. Diagram 152 explains the nature of the multiple scattering paths manifesting as an incoherent field in the distortion matrix. These occur upstream of the focal plane since they have the same times of flight as the waves that are singly scattered in the focal plane. A focused distortion matrix is shown in image 153 and corresponds to a line of the ultrasound image (z=75 mm). Its central row (r.sub.in=r.sub.out) clearly shows an overintensity characteristic of the single scattering contribution. However, it has an incoherent background of homogeneous intensity regardless of distance |r.sub.in−r.sub.out|. This incoherent background may be either caused by the noise in the ultrasound measurements, or linked to multiple scattering. To know its origin, the property of spatial reciprocity of waves may be examined. If R(r.sub.in, r.sub.out)=R(r.sub.out, r.sub.in), then it has to do with a signal of physical origin, hence multiple scattering. Otherwise, it can only be noise.

(73) The focused distortion matrix D.sub.rr provides access to the spread function 154 of the imaging system at any point of the field of vision. A multiple scattering rate γ(r.sub.in) may be measured locally on the basis of the ratio of the level of the incoherent background (|r.sub.out−r.sub.in|>>δ) to the signal level on the central row of D.sub.rr (r.sub.out=r.sub.in):

(74) $\begin{array}{cc}\gamma \left(r_{\mathrm{in}}\right)=\frac{I_M\left(r_{\mathrm{in}}\right)}{I_S\left(r_{\mathrm{in}}\right)+2I_M\left(r_{\mathrm{in}}\right)}=\frac{.Math.{.Math.D\left(r_{\mathrm{out}}\ne r_{\mathrm{in}},r_{\mathrm{in}}\right).Math.}^2.Math.}{{.Math.D\left(r_{\mathrm{out}}=r_{\mathrm{in}},r_{\mathrm{in}}\right).Math.}^2}& \left(34\right)\end{array}$
where the symbol . . . denotes an average along each row of D.sub.rr apart from the central elements (|r.sub.out−r.sub.in|>>δ). I.sub.M(r.sub.in)=|D(r.sub.out≠r.sub.in, r.sub.in)|.sup.2 is the multiple scattering intensity (incoherent background) at the point r.sub.in. I.sub.S (r.sub.in) is the single scattering intensity at the point r.sub.in. Equation (34) indicates that, on the central row (r.sub.out=r.sub.in), there is the following relationship: |D(r.sub.out=r.sub.in, r.sub.in)|.sup.2=I.sub.S(r.sub.in)+2I.sub.M(r.sub.in). The factor 2 comes from the phenomenon of coherent backscattering which consists of an overintensity by a factor of 2 of the backscattered wave on the virtual source in r.sub.in for the multiple scattering contribution. This phenomenon is caused by constructive interference between multiply scattered waves that have followed the same trajectory but in an opposite direction (reciprocal paths) [A. Aubry et al., Appl. Phys. Lett. 92, 124101, 2008].

(75) Image 155 of FIG. 15 shows the multiple scattering rate measured at each point of the ultrasound image 151. At shallow depths, this degree is non-negligible but linked to electronic noise (non-reciprocity of signals). Its increase at great depths is, on the other hand, clearly linked to the phenomenon of multiple scattering. In particular, the multiple scattering rate exceeds 50% beyond a depth of 60 mm.

(76) This measurement is confirmed by previous studies carried out in-vivo [A. Aubry et al., J. Acous. Soc. Am. 129, 225-233, 2011]. The contribution of the present invention resides in the fact that the focused distortion matrix allows local measurement of the multiple scattering rate whereas until now there was access only to an averaged value at each echo time or depth. [A. Aubry et al., J. Acous. Soc. Am. 129, 225-233, 2011].

(77) By taking times longer than the ballistic time [equation (4)], the spread function may also make it possible to follow the growth of the diffusive halo within the medium and derive therefrom a local measurement of the transport parameters of the multiply scattered wave (scattering coefficient or mean free path of transport). Studying the diffusive halo at short times gives access to a much finer spatial resolution than that obtained by the prior art in random scattering tomography [A. Aubry et al., Appl. Phys. Lett. 92, 124101, 2008].

(78) Although described through a number of detailed exemplary embodiments, the methods and systems for non-invasive ultrasound characterization comprise different variants, modifications and refinements which will be obviously apparent to those skilled in the art, it being understood that these different variants, modifications and refinements fall within the scope of the invention, as defined by the claims which follow.