Method and system for ultrasound characterization of media

Abstract

The invention relates to a method and a system for ultrasonically characterizing a medium. The method comprises the following steps: generating a series of incident ultrasonic waves; measuring a regular reflection matrix Rui (t) defined between the input transmit basis (i) and the output receive basis (u); determining a set of responses R of the medium obtained for a plurality of points in the region around the reference point at a plurality of frequencies f based on a regular reflection matrix Rui (t); determining a frequency correction rule, phi, suitable for the reference point and determined at the frequency f, based on the response of the medium at the plurality of points; for the plurality of frequencies f, a corrected response R'of the medium is determined by applying a frequency correction rule Phi to the response R of the medium around the reference point.

Description

Technical Field

[0001] This specification relates to methods and systems for ultrasonic characterization of heterogeneous media. These methods and systems use transducer arrays arranged to contact the medium to emit ultrasonic waves into the medium and measure the waves backscattered by the heterogeneous parts of the medium. Background Art

[0002] In the field of acoustic imaging, attempts are made to characterize a medium by actively probing the medium with ultrasonic waves. This is especially true for the principle of echo imaging in medical imaging.

[0003] Figure 10 A conventional focusing process for generating an ultrasonic image of a medium in the medium is shown. In this example, the medium includes an aberrant layer having a sound velocity different from that of other parts of the medium. This results in spatial distortion and temporal dispersion (reverberation) of the acoustic wavefront, leading to lateral and axial aberrations in the resulting ultrasonic image. These phenomena result in a decrease in resolution and contrast, as well as the appearance of reverberation artifacts, which are particularly inconvenient in, for example, medical examinations.

[0004] In the left figure, a transducer array arranged opposite the medium is used for acoustic transmission and imaging of the medium. The conventional method is to use a technique called beamforming, using focused transmission to perform acoustic transmission of the medium. A set of appropriate delays τ based on a uniform velocity model c 0 is applied to the signals emitted by each transducer so that the waves generated by each transducer undergo constructive interference at the target focal spatial position r in =(x in ,z). Due to the physical limitations of diffraction, the ultrasonic waves emitted through the aperture of the ultrasonic probe are concentrated in a region commonly referred to as the "focal spot". In addition, the waves passing through the aberrant layer are distorted and reflected several times, generating multiple echoes in the focal direction during transmission.

[0005] In the middle figure of this figure, the waves reflected at the focus return to the transducer array and then pass through the aberrant layer again, which further distorts the ultrasonic waves and results in an increase in echoes due to multiple reflections. Therefore, each scatterer in the medium generates multiple echoes, which generate multiple time pulses arriving at each transducer at different times, resulting in significant temporal dispersion of the ultrasonic signal. The beamforming process applied to such signals can be used to construct an ultrasonic image with significant axial distortion: as shown in the ultrasonic image in Figure 14 , the same scatterer appears at multiple depths. In addition, the non-uniform distribution of the sound velocity in the tissue being traversed affects the quality of the constructed image.

[0006] Figure 10The right figure in [Figure] shows how the signals received by each transducer for a single scatterer can be time-deconvolved and then time-reversed, so as to optimally focus the ultrasonic waves spatially and temporally on the scatterer involved. This time-reversal focusing technique still has limitations because it requires the medium to contain only a small number of scatterers. In ultrasonic echo imaging, the challenge is completely different because there is a medium with a large number of randomly distributed sub-resolved scatterers, and these large numbers of sub-resolved scatterers generate ultrasonic speckles.

[0007] Therefore, in conventional imaging, the common assumption of a homogeneous medium with a constant sound speed c is often not satisfied. 0 As a result, multiple internal reflections occur on the wave path towards the focus, generating reverberation. The result is a spatio-temporal distortion of the acoustic wavefront, which causes severe aberrations in the ultrasonic image and thus a decrease in its resolution and contrast. These aberrations can be harmful to ultrasonic characterization.

[0008] Document WO 2020 / 016250 proposes a technique for correcting aberrations in ultrasonic imaging based on post-processing operations on the reflection matrix of the medium. However, the method described in this document only considers IQ ultrasonic signals windowed in time around the expected ballistic time. Therefore, applying a phase shift to these signals to correct aberrations is equivalent to applying a time delay, and the magnitude of this delay must be lower than the time resolution of the ultrasonic signal. Therefore, the technique described in document WO2020 / 016250 is applicable to relatively low-order aberration corrections without time dispersion. Therefore, it cannot be used to correct reverberation or multiple scattering problems, which require determining different delay laws for each frequency component of the ultrasonic signal. Summary of the Invention

[0009] The object of the present disclosure is to improve known ultrasonic detection methods in order to correct aberrations in particular.

[0010] In a first aspect, the present disclosure relates to a method for ultrasonic characterization of a medium for medical analysis, the method comprising the following steps:

[0011] - Generating, by means of a transducer array, a series of incident ultrasonic waves in a region of the medium, the series of incident ultrasonic waves being the transmission basis (i); and

[0012] - Measuring a regular reflection matrix R(t) defined between the input transmission basis (i) and the output reception basis (u), the coefficients of the regular reflection matrix corresponding to the signals received by the transducer and caused by the ultrasonic waves reflected in the medium. ui The method further comprises the following steps:

[0013] The method further comprises the following steps:

[0014] - Determine (S130) a set of responses of the medium, these responses being for the sound speed model c 0 , based on the regular reflection matrix R ui (t), by, for a plurality of frequencies f of the signal received from the reflected ultrasonic wave, for a reference space position r p =(x p ,z p ), obtained by a focusing process for a plurality of points with a space position r=(x,z) in a region around a reference point;

[0015] - Determine (S140) a frequency correction law Φ by averaging or correlating the responses of the medium at different spatial position points (x, z) around the reference point, the frequency correction law being suitable for the reference point and determined at the frequency f;

[0016] - For a plurality of frequencies f, determine (S180) a corrected response R′ of the medium by applying the frequency correction law Φ to the response R of the medium around the reference point.

[0017] By means of these steps, the method advantageously enables local detection of the medium to obtain a local estimate of an appropriate frequency correction law for correcting the ultrasonic focusing process. This correction is used to reduce or eliminate, for example, the aberration caused by multiple reflections of waves generated by one or more aberration regions in the medium.

[0018] The estimate is calculated based on the measurements acquired and recorded in the regular reflection matrix. Therefore, the estimate can be calculated independently of the measurement acquisition phase, in particular by modifying various calculation parameters, so as to enable various ultrasonic characterization analyses either in real time or afterwards.

[0019] These correction calculations benefit from the local information extracted around the reference point and the local frequency information from the medium.

[0020] The method can be used in medical or veterinary imaging and all fields of ultrasonic imaging.

[0021] According to various embodiments of the method, any one of the following techniques can also be used.

[0022] According to a variant, the method further includes the following steps:

[0023] - Determine (S190) the intensity I of the echo imaging image point corresponding to the reference point by combining the corrected responses at a plurality of frequencies f of the reference point with a space position r p . c

[0024] According to a variant: ​

[0025] - Determine (S130) a set of responses R of the medium, by determining the response obtained by a focusing process between a first point at a spatial position r in =(x in , z) and a second point at a spatial position r out =(x out , z), where the first point corresponds to an input virtual transducer, the second point corresponds to an output virtual transducer, and the first and second points are identical (r in =r out ); and

[0026] Record all the responses R in a confocal reflection matrix R, the coefficients of which can be written as R = [R(x, z, f)];

[0027] - Perform the determination (S140) of a frequency correction law Φ, the coefficients of which are written as Φ = [φ(f, r p )], by correlating the responses at different spatial position points (x, z) around a reference point of the medium;

[0028] - Perform the determination (S180) of the corrected responses R' around the reference point, by applying the frequency correction law at each frequency f, by performing a term product between the confocal reflection matrix R and the phase conjugate of the frequency correction law Φ, i.e. by:

[0029]

[0030] where:

[0031] Record all the corrected responses R' in a corrected confocal reflection matrix R', the coefficients of which can be written as R' = [R'(x, z, f)];

[0032] . The symbol is the Hadamard product, such that:

[0033] R'(x, z, f) = R(x, z, f) φ * (f, r p ).

[0034] According to a variant, the method further comprises the following steps:

[0035] - Determine (S190) the intensity I of a point by combining the corrected responses R' of the echo imaging image point at a spatial position (x, z) at a plurality of frequencies f, i.e. by: c , i.e. by:

[0036] I c (x, z) = |∑ f R'(x, z, f)| 2 .

[0037] According to a variant:

[0038] Determining the frequency correction rule (S140) includes:

[0039] Based on the confocal reflection matrix R(z, f), constructing (S141) the correlation matrix C; and analyzing (S142) the correlation matrix C to determine the frequency correction rule Φ.

[0040] According to a variant:

[0041] The correlation matrix C is determined in the frequency domain by the following formula:

[0042] C(f, f′) = ∑ x，z R(x, z, f)R * (x, z, f′)

[0043] where: R is the confocal reflection matrix;

[0044] x, z are the coordinates of the points in the region around the reference point;

[0045] * is the conjugate operator.

[0046] According to a variant:

[0047] The correlation matrix C is determined in the basis of the image points at the spatial position (x, z) by the following formula:

[0048] C({x, z}, {x′, z′}) = ∑ f R(x, z, f)R * (x′, z′, f)

[0049] where: R is the confocal reflection matrix:

[0050] x, z are the coordinates of the image points in the region around the reference point;

[0051] * is the conjugate operator.

[0052] According to a variant:

[0053] Analyzing (S142) the correlation matrix C is the eigenvalue decomposition of the correlation matrix C, and the frequency correction rule Φ is the first eigenvector U of the correlation matrix C 1 .

[0054] According to a variant:

[0055] Analyzing (S142) the correlation matrix C is to solve the equation involving the correlation matrix C and the frequency correction rule Φ, and this equation solving corresponds to iterative time reversal or iterative phase reversal.

[0056] According to a variant:

[0057] An optimization algorithm that maximizes the confocal intensity in the region around the reference point of the ultrasound image is used to perform the determination (S140) of the frequency correction rule.

[0058] According to one variant:

[0059] The steps (S140, S180) of the correction process are iterated multiple times; and,

[0060] wherein, in each iteration, the corrected confocal reflection matrix R’(z, f) obtained in the previous iteration is used instead of the confocal reflection matrix R(z, f).

[0061] According to one variant:

[0062] In each iteration, the size of the region around the reference point with spatial position r p is reduced.

[0063] According to one variant:

[0064] Determining (S130) the confocal reflection matrix R(z, f) includes compensating for the time decay of the signal.

[0065] According to one variant:

[0066] - Determining (S130) a set of responses R of the medium includes determining the response obtained through the focusing process between a first point with spatial position ri n =(x in , z) and a second point with spatial position r out =(x out , z), where the first point corresponds to the input virtual transducer, the second point corresponds to the output virtual transducer, the first point and the second point are at the same depth z in the region, and the lateral positions x in and x out form the focusing basis (x) at each depth z; and

[0067] Recording all responses R in the focusing reflection matrix R xx (z, f), the coefficients of which can be written as R xx (z, f)=[R(x in , x out , z, f)];

[0068] - Determining (S140) the frequency correction rule Φ includes the following sub-steps implemented at each depth z and each frequency f:

[0069] - Determining (S150) the double reflection matrix R by the forward projection of the focusing reflection matrix R xx (z, f) onto the correction basis (c)c (z, f);

[0070] - Based on the double reflection matrix R c (z, f), calculate (S160) the frequency correction rule Φ, which is determined on the correction basis (c), Φ = [φ(c, f, r p )], so that the frequency correction rule Φ is a space-frequency correction rule;

[0071] - Determine (S170) the corrected double reflection matrix R' around the reference point c , and the coefficients of this matrix are written as R' c (z, f) = [R' c (x, c, f, z)], which is determined by performing the term product between the double reflection matrix R c (z, f) and the phase conjugate of the frequency correction rule Φ, that is, by:

[0072]

[0073] Where:

[0074] * The symbol refers to the phase conjugate operation;

[0075] The symbol is the Hadamard product, so that:

[0076] R' c (x, c, f, z) = R c (x, c, f, z) Φ * (x, c, f, z)

[0077] - Determine (S180) the corrected response R' of the medium around the reference point, including determining the corrected focused reflection matrix R' c (z, f) by the back-projection of the corrected double reflection matrix R' xx (z, f) on the focusing basis (x).

[0078] According to a variant, the method further includes the following steps:

[0079] - By combining the diagonal coefficients of the corrected focused reflection matrix R' xx (z, f) of the ultrasonic image points at spatial positions (x, z) at multiple frequencies f, determine (S190) the intensity I c of this point, that is:

[0080] I c (x, z) = |∑ f R' xx (x, x, z, f)|2 。

[0081] According to a variant:

[0082] The forward projection (150) is implemented by the matrix product between the transition matrix and the focusing reflection matrix R xx (z, f), i.e.:

[0083] R c (z, f) = P(z, f) × R xx (z, f)

[0084] where: P(z, f) = [P(c, x, z, f)] is the transition matrix between the focusing basis (x) and the correction basis (c) at depth z for each frequency f.

[0085] According to a variant:

[0086] The correction basis (c) is the input correction basis or the output correction basis.

[0087] According to a variant:

[0088] Calculating the frequency correction rule (S160) includes:

[0089] Based on the double reflection matrix R c (z, f), constructing (S161) the correlation matrix C; and

[0090] Analyzing (S162) the correlation matrix C to determine the frequency correction rule Φ.

[0091] According to a variant:

[0092] The correlation matrix C is determined in the correction basis (c) in the frequency domain by the following formula:

[0093]

[0094] where: R c is the double reflection matrix;

[0095] R ref is the model reflection matrix of the model medium, in which the sound speed is c 0 (the expected sound speed of the medium), and the plane reflector is located at depth z;

[0096] x, z are the coordinates of the points in the region around the point r p ;

[0097] * is the conjugate operator.

[0098] According to a variant:

[0099] In the basis of the image points at the spatial position (x, z), the correlation matrix C is determined by the following formula:

[0100]

[0101] where: R c is the double reflection matrix;

[0102] R ref is the model reflection matrix of the model medium, in which the sound speed is c 0 (the expected sound speed of the medium), and the plane reflector is located at the depth z;

[0103] x, z are the coordinates of the image points in the area around the point r p ;

[0104] * is the conjugate operator.

[0105] According to a variant:

[0106] The analysis (S162) that the correlation matrix C is the eigenvalue decomposition of the correlation matrix C, and the frequency correction rule Φ is the first eigenvector U of the correlation matrix C 1 .

[0107] According to a variant:

[0108] The analysis (S162) that the correlation matrix C is the solution of the equation involving the correlation matrix C and the frequency correction rule Φ, and this equation solution corresponds to iterative time reversal or iterative phase reversal.

[0109] According to a variant:

[0110] The frequency correction rule (S160) is calculated using an optimization algorithm that maximizes the confocal intensity in the area around the reference point of the ultrasonic image.

[0111] According to a variant:

[0112] The steps (S140, S180) of the correction process are iterated multiple times; and

[0113] where, in each iteration, the forward projection (S150) uses the corrected focusing reflection matrix R' xx (z, f) obtained during the back projection (S180) of the previous iteration, rather than the focusing reflection matrix R xx (z, f).

[0114] According to a variant:

[0115] In each iteration of the forward projection (S150), it alternates between the forward projection in the input correction basis and the forward projection in the output correction basis.

[0116] According to a variant:

[0117] In each iteration, the correction basis (c) of the forward projection (S150) is different.

[0118] According to a variant:

[0119] In each iteration, the size of the region around the reference point with a reduced spatial position r p is reduced.

[0120] According to a variant:

[0121] Determine (S130) the focusing reflection matrix R xx (z,f) including compensating for the temporal decay of the signal.

[0122] According to a second aspect, the present disclosure relates to an ultrasonic characterization system for analyzing a medium (for example, in the medical field), which is configured to implement the method as described above. The system according to this second aspect includes:

[0123] - A transducer array adapted to generate a series of incident ultrasonic waves in a region of the medium and, over time, measure the ultrasonic waves backscattered by the region; and

[0124] - A computing unit connected to the transducer array and adapted to implement the method according to the first aspect. BRIEF DESCRIPTION OF THE DRAWINGS

[0125] By the following detailed description, other advantages and features of the above technology will become apparent. This detailed description is presented in a non-limiting manner for illustrative purposes and with reference to the following drawings:

[0126] Figure 1 (a) to (h) thereof show transmission / reception sequences for ultrasonic imaging and ultrasonic characterization of a medium;

[0127] Figure 2 show conventional confocal imaging and matrix imaging, in which a focusing reflection matrix is synthesized, which contains the responses between individual virtual transducers synthesized at each depth by beamforming;

[0128] Figure 3 show an example of an ultrasonic characterization system for implementing the method according to the present disclosure;

[0129] Figure 4 show the definitions used in the method according to the present disclosure;

[0130] Figure 5is a flowchart of an ultrasonic characterization method according to the present disclosure, wherein a frequency correction rule is determined based on responses at multiple points in a medium at multiple frequencies, and the response of the medium is corrected by applying the frequency correction rule;

[0131] Figure 6 is a flowchart of an embodiment that performs the determination of the frequency correction rule of the method in Figure 5 by projecting a focused reflection matrix into a correction basis;

[0132] Figure 7 shows the specific calculation of the frequency correction rule by constructing and analyzing a correlation matrix; Figure 6 of

[0133] Figure 8 shows the specific calculation of the frequency correction rule by local ultrasonic image optimization; Figure 6 of

[0134] Figure 9 shows Figure 5 a flowchart of an embodiment of the steps for determining the frequency correction rule of the method in

[0135] Figure 10 which determines the frequency correction rule by constructing and analyzing a correlation matrix in the case of a confocal reflection matrix;

[0136] Figure 11 shows the aberration layer problem in ultrasonic imaging; Figure 10 shows a first embodiment of the method according to the present disclosure, which solves the reverberation problem in matrix form and generalizes the process described in

[0137] Figure 12 for the case of isolated scatterers to the "speckle" case;

[0138] Figure 13 shows an experimental medium of the ultrasonic phantom type on which the method according to the present disclosure is tested;

[0139] Figure 14 shows Figure 13 an ultrasonic image of the region of interest in the medium in

[0140] Figure 15 which shows reverberation artifacts that are particularly evident at the echo scatterers of the calibrated experimental medium; Figure 14 in region B1 of the image in

[0141] Figure 16An uncorrected ultrasound image is shown in Fig. (a), and an image obtained after correction using the method according to the present disclosure is shown in Fig. (b);

[0142] Figure 17 The same ultrasound image is shown, where three regions are selected to locally apply the method according to the present disclosure; Figure 14 The spatial frequency correction rules obtained in regions C1, C2, and C3 of the image in are shown;

[0143] Figure 18 The image obtained in is shown; Figure 17 In Figs. (a), an uncorrected ultrasound image is shown, in Fig. (b), an ultrasound image obtained after overall correction is shown, and in Fig. (c), ultrasound images obtained after correction in three regions C1, C2, and C3 are shown.

[0144] Figure 19 In each of the embodiments described with reference to the accompanying drawings, unless otherwise specified, similar or identical elements have the same reference numerals.

[0145] DETAILED DESCRIPTION DETAILED DESCRIPTION

[0146] In the following detailed description, for the sake of clarity, only certain embodiments are described in detail, but these examples are not intended to limit the overall scope of the principles that are obvious from the present disclosure.

[0147] In the various embodiments and aspects described in the present disclosure, they can be combined or simplified in various ways. In particular, unless otherwise specified, the steps in each process can be repeated, interchanged, and / or executed in parallel.

[0148] The present disclosure relates to methods and systems for ultrasonic characterization of a medium, particularly suitable for medical imaging of living or non-living tissues. The medium is, for example, a heterogeneous medium, and an attempt is made to characterize the heterogeneous medium to, for example, identify and / or characterize heterogeneities, provide accurate information about the medium under study, and detect unhealthy or damaged regions and / or tissues. For example, such data is very useful for medical applications, such as identifying lesions in the breast area, damaged muscle tissue, or deteriorating liver tissue. As is well known, these characterization techniques are non-invasive for the medium, and the medium is advantageously preserved, particularly in terms of its properties and integrity.

[0149] Common ultrasonic imaging modalities

[0150] In the field of ultrasonic imaging, it is often attempted to construct an image of the medium reflectivity based on the echoes backscattered by the inhomogeneities in the medium. This is the principle of ultrasonic scanners used in medical imaging, which are particularly used to view the internal anatomical structures of objects, humans or animals. For the purpose of simplification, in order to construct an ultrasonic image, the medium is considered to be homogeneous with a constant sound propagation speed c 0 .

[0151] Conventional ultrasonic methods typically use a piezoelectric transducer array, where the transducers are capable of emitting and / or receiving ultrasonic signals independently or almost independently, and each transducer is located at a position u in a strip supporting the array. The transducer array is placed opposite the medium to acoustically transmit the medium in different ways and construct a representative image thereof. The conventional method consists in using a technique called beamforming to acoustically transmit the medium with a focused emission. This method consists in applying a set of appropriate delays τ(u 0 , x in , z, c in ) to the signals emitted by each transducer, so that the wavelets generated by each transducer constructively interfere at the target focus at the spatial position (x 0 , z). Due to the physical limitations of diffraction, the ultrasonic waves are emitted through the apertures of the ultrasonic probe and are concentrated in a region with a lateral width of δx, usually called the "focus spot". in

[0152] During reception, a digital focusing step is also performed so that an ultrasonic image showing the characteristics of the medium under study can be constructed. By shifting the echoes collected by the transducer array in time, they are brought back to their original phase. The delay τ(u out , x out , z, c 0 ) is the same as the delay applied during transmission, and the variable u out indicates the position of each transducer. During the transmission phase, if the velocity model c 0 used corresponds to the actual situation of the medium under study, all the signals interfere at the position point (x in , z) at the ballistic time t = z / c 0 . During reception, by summing the signals at the echo time t = 2z / c 0 , the signals from the same point (x out = x in ​) interfere with the signal. This summation gives the final result of the received focusing. This confocal method with both transmitting and receiving focusing enables direct imaging of the reflectivity of the medium with a lateral resolution of δx and good contrast. However, this method is very time-consuming because it requires the transmitter to physically focus on each point of the medium, or at least at a given depth, on each line that will be constructed as an image representative of the medium.

[0153] Matrix method for ultrasonic imaging

[0154] In recent years, the present applicant has developed a matrix method for ultrasonic imaging. This method is based on constructing a canonical reflection matrix of the medium under study. This canonical reflection matrix is determined by experiments or calculations or digital simulations reproducing the experiments.

[0155] The first proposal for creating this canonical reflection matrix is to successively emit ultrasonic pulses from each transducer in the array, as Figure 1 shown in Figure (a), where the position of each transducer is marked with the coordinate u in . This generates a diverging cylindrical (or spherical) incident wave. This wave is reflected by the scatterers in the medium, and through each transducer, the backscattered field as a function of time is measured, as Figure 1 shown in Figure (b). By repeating this operation for each transducer successively used as a source, the canonical reflection matrix R uu (t) = [R(u out , u in , t)] is determined, which is composed of all the impulse responses R(u out , u in , t) between each transducer. This matrix thus contains a large amount of information describing the medium under study. However, this method assumes that the medium remains static during the measurement, which is particularly difficult in the case of a medical examination of a living patient. Moreover, since the medium is acoustically transmitted with a single transducer, the signal-to-noise ratio of the recorded signal is not good.

[0156] The second alternative for constructing the canonical reflection matrix is to acoustically transmit the medium using waves successively focused on multiple points in the medium, as is usually the case in a standard ultrasonic scanner. Figure 1 Figure (c) shows this transmit focusing illumination. The delay law is applied to each transmitted signal of this focusing. For each focused illumination, through all the position sensors u out , the reception shown in Figure 1 Figure (d) is measured, and each response forms the canonical reflection matrix R(u out , r in , t). However, this multi-point focusing is still too slow.

[0157] A third way to construct the regular reflection matrix consists in acoustically transmitting a series of plane waves through the medium. This method avoids the inherent problems in most of the above methods. Figure 1 Figure (e) shows the principle of this plane-wave illumination. The delay law τ' is applied to each emitted signal to form a wavefront that is inclined at an angle θ with respect to the transducer array. in At reception, as shown in Figure 1 Figure (f), for each incident plane wave θ in , through all u out position sensors, the domain R(u out , θ in , t) backscattered by the medium is measured. This set of responses forms the regular reflection matrix R uθ (t) = [R(u out , θ in , t)]. The above-mentioned dual-focusing method can be implemented digitally by time-shifting the measured signals before coherently summing them. This method enables ultrafast imaging and elastography, as described in the following references:

[0158] "Coherent plane-wave compounding for very high frame rateultrasonography and transient elastography", G. Montaldo et al. (IEEE Trans. Ultrasonics, Ferroelect., Freq. Control 56 489 - 506, 2009).

[0159] A third alternative for creating the regular reflection matrix consists in acoustically transmitting the medium with a diverging wave basis, as shown in Figure 1 Figure (g) and Figure 1 Figure (h), which allows a larger area of the acoustic domain to be illuminated than using plane waves. This technique is especially used for super-resolution imaging, as explained in the following references:

[0160] "Ultra fast imaging of the heart using Circular Wave SyntheticImaging with Phased Arrays", Couade et al, IEEE International UltrasonicsSymposium (2009).

[0161] Thus, conventional imaging consists in, as shown in Figure 2As shown in (a), for each pixel in the image, double focusing is performed at the same focal point (x in = x out ) during both transmission and reception.

[0162] As shown in figure (b) of Figure 2 , matrix imaging consists in decoupling the transmission and reception foci for the same echo time t. Thus, the focused reflection matrix R xx (τ) = [R(x in , x out , z, τ)] is determined, which is composed of the responses between an input virtual transducer and an output virtual transducer of the input and output synthetic focal spots corresponding to different spatial positions with spatial positions r in = (x in , z) and r out = (x out , z), and the desired depth z is determined by the echo time t: z = c 0 t / 2. Here, τ represents the time in the focusing basis such that τ = t - 2z / c 0 . The starting point of this time τ corresponds to the moment when the virtual source emits an ultrasonic pulse. Based on the regular reflection matrix R ui (t), the focused reflection matrix is obtained by focusing (e.g., using beamforming). This matrix provides more information characterizing the medium under study than the information obtained using the confocal image mode of conventional imaging. The focused reflection matrix can be synthesized in the time domain using the delay law (beamforming), or in the frequency domain by applying an appropriate phase shift.

[0163] In particular, in the following literature:

[0164] "Reflection Matrix Approach for Quantitative Imaging of ScatteringMedia", William Lambert, et al, Phys. Rev. X 10, 021048, (2020),

[0165] And then in the following literature:

[0166] "Ultrasound Matrix Imaging - Part I: The Focused Reflection Matrix, theF-Factor and the Role of Multiple Scattering", William Lambert, et al, IEEETrans. Med. Imag. 41, 3907 - 3920, (2022),

[0167] Describes a focused reflection matrix.

[0168] In these publications, it is considered at the ballistic time τ = 0 (t = 2z / c 0 ), a focused reflection array R between virtual transducers xx (z, τ). The diagonal coefficients of the matrix R xx (z, τ = 0) (x out = x in) enable the construction of a synthetic confocal image at depth z, while the off-diagonal coefficients provide information on potential aberrations and multiple scattering effects that may modify the quality of this image.

[0169] In the patent application WO2020 / 016250, a focused reflection matrix taking into account the ballistic time is considered and this matrix is studied to correct in particular the aberration problem. However, the corrections applied correspond to applying time delays, the amplitude of which must be less than the time resolution of the ultrasonic signal. Thus, the technique described in WO 2020 / 016250 is only applicable to relatively low-order aberration corrections without time dispersion. It cannot therefore be used to correct reverberation or multiple scattering problems, which require determining different delay laws for each frequency component of the ultrasonic signal.

[0170] Ultrasonic characterization system

[0171] Figure 3 Shows an example of an ultrasonic imaging system 1 for implementing an ultrasonic imaging method of a medium (e.g., a heterogeneous medium) according to the present disclosure. The system and method enable the formation of an ultrasonic image of at least a part (region of interest or field of view) of the medium.

[0172] System 1 includes:

[0173] - A detection system 20 or a probe 20;

[0174] - A computing unit 30 for computing an image based on signals from the probe 20;

[0175] - A control panel 40 connected to the computing unit 30, the control panel including for example keys 41 and a touchpad 42;

[0176] - A display device 50 for viewing the image and various elements or measurements.

[0177] The probe 20 is connected to the computing unit 30 by a cable 21 or by a wireless connection and is capable of emitting ultrasonic waves W into the medium M and receiving ultrasonic waves W from the medium M, the ultrasonic waves being originated from the reflection of the emitted ultrasonic waves on scattering particles (i.e., scatterers) within the medium.

[0178] The probe 20 may include an array 10, which includes a plurality of transducers 11. The array 10 is, for example, a linear or curved or two-dimensional or matrix array. The transducers 11 are capable of converting an electrical signal into vibrations and converting vibrations into an electrical signal. The transducers 11 are, for example, piezoelectric ultrasonic transducers, which may be in the form of a rigid strip that is in direct or indirect contact with the outer surface of the medium M to couple with the medium and vibrate and transmit and receive ultrasonic waves W. The array 10 of the transducers 11 of the probe 20 is then associated with the computing unit 30. The transducer array 10 may have one hundred or more transducers 11.

[0179] The computing unit 30 may have a housing 31, which includes receiving means for amplifying and / or filtering the signals received from the probe 20, and a converter (analog-to-digital converter and digital-to-analog converter) for converting the signals into data representative of the signals. The data may be stored in a memory in the computing unit 30 and / or be directly processed to calculate intermediate data (beamforming data or other data). The computing unit 30 may implement any process for constructing an image based on the signal data received from the probe 20 later, such as beamforming.

[0180] The calculated image may be:

[0181] - of the medium, usually in grayscale to view an image of an organ in the medium (B-mode image); and / or

[0182] - an image showing the velocity or flow in the medium (color image), which may be used, for example, to view blood vessels and associated flow in the medium; and / or

[0183] - an image showing the mechanical characteristics of the medium (such as elasticity according to the elastography mode, elasticity), which may be used, for example, to identify tumors within the medium.

[0184] The "connection" or "link" between the detection system 20, the computing unit 30, and the display device 50 refers to any type of wired connection of electrical or optical type, or any type of wireless connection using a protocol such as WIFi TM , or Bluetooth TM . These connections or links are either unidirectional or bidirectional. The associated display device 50 may be of any type, such as a touch or non-touch screen, connected or not.

[0185] The display device 50 is a screen for displaying the image calculated by the computing unit 30. The display device 50 may also display other information, such as an image scale, or configuration information for calculation or processing, or any measurement values or help information. The screen 50 may be hinged on a support arm 51 to improve its positioning for the user. The screen 50 is usually a large screen (at least 20 inches) so that it is easier for the user to view.

[0186] The control panel 40 is part of, for example, a system housing, and this part includes a panel housing having a substantially flat surface 40a that is inclined towards the user for single-handed operation. As Figure 3 shown, the control panel 40 may include a control screen 49 for displaying various configuration information.

[0187] The computing unit 30 is configured to perform computing and / or processing steps, in particular to perform the methods according to the present disclosure. Conventionally, as shown for the array 10 of transducers 11 located on the surface of the medium M, Figure 4 by adopting a first axis X and a second axis Z perpendicular thereto, a spatial reference system of the medium M is defined. For simplicity, the first axis X corresponds to the transverse direction along which the transducers 11 are aligned in the example of a linear array, and the second axis Z corresponds to the depth of the medium M relative to the array 10 of transducers 11. This definition can be adjusted according to the situation. Thus, for example, in the case of a two-dimensional array 10, this definition can be extended to a three-axis spatial reference system, or in the case of a curved array 10 to a polar coordinate reference system, or to any other reference system adapted to and / or depending on the structure and shape of the ultrasonic transducer array 10. Thus, in the following of the present disclosure, for simplicity of description, the Cartesian reference system XZ corresponding to the linear probe 20 will be used, but those skilled in the art can easily generalize and apply the results to any type of reference system.

[0188] In the following of the present disclosure, when referring to the array 10 of transducers 11 for transmitting and receiving, it should be understood that in a more general case, multiple transducer arrays can be used simultaneously. The transducers 11 can be both transmitters and receivers, or some are only transmitters and some are only receivers. Similarly, the array 10 can be composed of one (1) to N transducers 11 of the same or different types.

[0189] The array 10 of transducers 11 is, for example, used both as a transmitter and as a receiver, or is composed of multiple transducer sub-arrays, some dedicated to transmitting ultrasonic waves and some dedicated to receiving ultrasonic waves. A transducer array refers to at least one transducer, an aligned or non-aligned sequence of transducers, or a two-dimensional distribution of transducers (such as a transducer matrix), or any spatial distribution of transducers.

[0190] When specifically referring to computing or processing steps for implementing method steps in the present disclosure, it should be understood that each computing or processing step can be implemented by software, hardware, firmware, microcode, or any suitable combination of these or related technologies. When using software, each computing or processing step can be implemented by, for example, computer program instructions or code that can be interpreted or run. These instructions can be stored or sent to a computer (or computing unit) readable storage medium and / or run by a computer (or computing unit) to implement these computing or processing steps.

[0191] Analyze a point in a medium using a focusing reflection matrix

[0192] This disclosure describes methods and systems for ultrasonic characterization of a medium. In practice, the medium is assumed to be heterogeneous. These methods and systems are based on Figure 4 the definitions shown.

[0193] In the medium, define:

[0194] - A first point P1 with a spatial position r in in the spatial reference frame of the medium; and

[0195] - A first point P2 with a spatial position r out in the spatial reference frame of the medium.

[0196] Mark these spatial positions r in and r out in bold to indicate that these elements are position vectors, vectors taken in the spatial reference frame (X,Z) of the medium. Other representations and definitions of point positions are possible and available to any ultrasonic technician.

[0197] In this disclosure, the first point P1 has a lateral position x in . The second point P2 has a lateral position x out . The two points P1, P2 have the same desired depth z in = z out , also denoted as z, which is controlled by the echo time t under consideration. Thus, the spatial positions of points P1 and P2 are r in = (x in ,z) and r out = (x out ,z).

[0198] Select the two points P1 and P2 to be fairly close, i.e., separated by a few millimeters, e.g., twenty (20) millimeters or less.

[0199] The ultrasonic characterization method implemented by the computing unit 30 of the system 40 includes the following steps:

[0200] - Generate, via the array 10 of the transducer 11, a series of incident ultrasonic waves US in in a region of the medium, the series of incident ultrasonic waves being transmitted at base i; and

[0201] - Measure or construct a regular reflection matrix R ui (t) defined between the transmission base i at the input and the reception base u at the output, the coefficients of the regular reflection matrix corresponding to the signals received by the transducer caused by the ultrasonic waves reflected by the scatterers in the medium.

[0202] - Determine a focusing reflection matrix R that includes the response of the medium xx (z), where the response is based on the input virtual transducer TV at spatial position r in and the regular reflection matrix R in between the output virtual transducer TV at spatial position r out is calculated by focusing the regular reflection matrix R out (t). ui (t).

[0203] The obtained regular reflection matrix R ui (t) can be a "real" matrix, i.e., composed of real coefficients in the time domain, and the electrical signals recorded by each transducer are real. Alternatively, the matrix can be a "complex" matrix, i.e., composed of complex numerical values, for example, in the case of demodulation for phase beamforming and IQ beamforming.

[0204] The focusing reflection matrix R can be expressed in different ways xx . In the expressions of the present disclosure, the first point P1 at spatial position (x in , z) is used as a reference. These expressions can also be established with respect to the second point P2 at spatial position (x out , z), or with respect to the midpoint at spatial position ((x in + x out ) / 2, z) between P1 and P2, or with respect to any point designated as a reference point. Those skilled in the art will be able to make the necessary changes to the variables in the illustrated expressions.

[0205] Calculate the medium response by focusing based on the regular reflection matrix R ui (t).

[0206] The focusing reflection matrix R xx (z) response corresponds to the acoustic pressure domain calculated between all points at lateral positions x 0 and x in at the desired depth z in the medium for echo time t and assumed sound speed c out . In other words, the focusing reflection matrix R xx (z) is defined as: R xx (z) = [R(x in , x out , z)].

[0207] The parameters of the depth z and the sound speed c in the medium 0 affect the delay law used in the focusing process.

[0208] As in the foregoing of the present disclosure, in Figure 2 Figure (a) to Figure 2As described in FIG. (f), the input emission basis i is, for example, a wave basis in which each wave is generated by a single transducer in the transducer array 10, or a basis of plane waves having an angular inclination θ with respect to the X-axis, or a basis of virtual sources.

[0209] The receiving basis u is, for example, the basis of the transducers 11. Alternatively, other receiving bases can be used during reception, such as a plane wave basis (or spatial Fourier basis), or any other basis that allows, through a suitable beamforming process during reception, the measurement of the ultrasonic field in the intermediate plane between the ultrasonic probe and the considered focal plane.

[0210] Therefore, between the emission basis i and the receiving basis u, an ultrasonic generation step is performed. Thus, this ultrasonic generation step is defined for any type of focused or non-focused ultrasonic waves (such as plane waves).

[0211] In the measurement step, a regular reflection matrix R ui (t) is defined between the input reflection basis i and the output receiving basis u. This matrix contains, at time t, all the time responses of the medium measured by each transducer 11 with spatial coordinates u out for each emission i in . It should be understood that the elements marked with the subscript "in" refer to emission (i.e., input), while the elements marked with the subscript "out" refer to reception (i.e., output). The regular matrix can also be recorded and / or stored, for example, in the memory of a computing unit or on any other medium allowing persistent or temporary storage (whether removable or non-removable).

[0212] In the step of determining the focused reflection matrix R xx (z), this matrix can be obtained as follows:

[0213] - An input focusing process based on the regular reflection matrix R ui (t), which uses the forward flight time of the wave between the emission basis i and the input virtual transducer TV in and generates a so-called input focal spot around a first point P1 with spatial position r in , where the input focal spot corresponds to the input virtual transducer TV in ;

[0214] - An output focusing process based on the regular reflection matrix R ui (t), which uses the return flight time of the wave between the virtual output transducer TV out and the transducers of the receiving basis u and generates a so-called output focal spot around a second point P2 with spatial position r out , where the output focal spot corresponds to the virtual output transducer TV out .

[0215] These input and output focusing processes form an input-output focusing process, which will be referred to as the focusing process hereinafter, or more simply as focusing.

[0216] In other words, in this ultrasonic characterization method, the virtual input transducer TV in corresponds to an ultrasonic "virtual source" located at the spatial position r in in the medium, while the virtual output transducer TV out corresponds to an ultrasonic "virtual sensor" located at the spatial position r out in the medium. The virtual source and the virtual sensor are spatially separated by the difference Δr = r out – r in in their spatial positions. In this example, they are separated only along the horizontal axis by Δx = x out - x in . For the sound speed model c 0 , their desired depth is the parameter z used in the focusing law. Their actual depth is determined by the axial position (depth) of the isochronous volume, i.e., by the echo time t and the sound speed distribution c(r) in the medium. The lateral dimension of the virtual transducer is determined by the focal spot generated by focusing at this actual depth.

[0217] Moreover, it is possible to:

[0218] - either determine or calculate the focusing reflection matrix R xx (z) in the time domain, in which case the matrix can be explicitly represented using the time parameter t, i.e., as R xx (z,t), and

[0219] in practice, calculate the data of the focusing reflection matrix R xx (z,t) between two predetermined time points;

[0220] - either determine or calculate the focusing reflection matrix R xx (z) in the frequency domain, in which case the matrix can be explicitly represented using the pulse parameter ω corresponding to the frequency f, i.e., as R xx (z,ω), where ω = 2π / f, and

[0221] in practice, calculate the data of the focusing reflection matrix R c for the central pulse ω xx (z,ω) between two pulses ω with a frequency bandwidth (e.g., between the lower pulse ω- and the upper pulse ω+).

[0222] Thus, the following method calculations can be performed in either the time domain or the frequency domain.

[0223] In the case of calculation in the time domain, focusing is achieved by using input and output beamforming calculations to obtain the focusing reflection matrix R in between the input virtual transducer TV out and the output virtual transducer TV xx (z, t). The focusing reflection matrix R xx (z, t) can be determined by the following coefficients:

[0224]

[0225] where:

[0226] N in is the initial normalization coefficient;

[0227] N out is the second normalization coefficient;

[0228] R ui (t) is the confocal reflection matrix, where:

[0229] R(u out , i in , t) is the element of the confocal reflection matrix R in out ui (t) recorded by the transducer at spatial position u ui after transmission with index i in and r out applied to it;

[0230] A in (i in , x in z) and A out (u out , x out , z) are predetermined apodisation coefficients, for example, used to maintain a constant numerical aperture during transmission and reception.

[0231] For example, the first normalization coefficient N in can be defined by the following: Similarly, the second normalization coefficient N out can be defined by the following:

[0232] τ in (i in , x in , z) is the spatial position where each incident wave i in reaches the model medium with a sound speed c 0 at (xin The expected flight time of the first focus of (x, y, z).

[0233] τ out (u out , x out , z) is the expected flight time of the wave reflected from the second spatial position focus (x out , z) to the position transducer u out .

[0234] These delay times τ in and τ out are usually calculated by qualified personnel based on the established sound speed model. A relatively simplified assumption is to assume a homogeneous medium in which there is a constant sound speed c 0 . In this case, the flight time is directly obtained based on the distance between the probe transducer and the virtual transducer. Thus, the calculation of these delay times depends on the type of wave, the assumed sound speed, and the transducer array geometry.

[0235] N in The number of elements in the transmit base is, for example, greater than or equal to one (1), and advantageously greater than or equal to two (2). N out The number of elements in the receive base is, for example, greater than or equal to two (2).

[0236] The above beamforming formula is thus a double summation of the time responses recorded in the regular reflection matrix R ui : the first summation is according to the transmit base i reflecting the transmit focus, and the second summation is according to the receive base u regarding the receive focus. This calculation is performed for two points P1 and P2, that is, for the spatial coordinates of the points with expected spatial positions r in =(x in , z) and r out =(x out , z) respectively.

[0237] The result of this beamforming formula is thus the time signal for these two spatial coordinates (r in , r out ) or for these two lateral positions x in , x out .

[0238] Finally, the focused reflection matrix R xx (z, t) expressed in the time domain can be converted into the focused reflection matrix R xx (z, ω) in the frequency domain through Fourier transform, that is, through:

[0239]

[0240] This Fourier transform can be implemented by any type of discrete Fourier transform (whether normalized or non - normalized).

[0241] In a second case calculated in the frequency domain, the regular reflection matrix R ui (t), since it is constituted by the signals received by the transducer, can be converted, in the frequency domain, by Fourier transform, into the regular reflection matrix R ui (ω), that is, by:

[0242]

[0243] This Fourier transform can be implemented by any type of discrete Fourier transform (whether normalized or non - normalized).

[0244] Thus, the focused reflection matrix R xx (z) or R xx (z, ω) of the medium can be obtained by focusing through the following matrix calculation, which is roughly equivalent to focusing using the above - mentioned time beamforming, that is, through the following matrix product:

[0245]

[0246] Where:

[0247] The matrix R ui (ω) is the Fourier transform of the regular reflection matrix R ui (t);

[0248] The matrix G ux (z, ω) is the reception transition matrix suitable for transitioning from the reception basis (u) to the focused basis (x) at depth z and pulse ω;

[0249] The matrix P ix (z, ω) is the transmission transition matrix suitable for transitioning from the transmission basis (i)″i″ to the focused basis (x) at depth z and pulse ω;

[0250] The symbols * and respectively refer to the matrix conjugate and transpose conjugate operations.

[0251] Aberration correction process

[0252] The aim of method S100 according to the present disclosure is to correct the aberration in the ultrasonic characterization of the medium M, the aberration being, for example, due to structural changes in the medium that cause changes in the speed of sound and reflectivity.

[0253] Figure 5 Shows method S100 according to the present disclosure, which is implemented by the calculation unit 30 of system 40 and includes the following steps:

[0254] - Step S110, in which a series of incident ultrasonic waves US are generated in the region of the medium by the transducer array 10 of transducer 11 in , the series of incident ultrasonic waves being the emission basis i; and

[0255] - Step S120, in which the regular reflection matrix R ui (t) defined between the emission basis i at the input and the reception basis u at the output is measured or constructed, the coefficients of this regular reflection matrix corresponding to the signals received by the transducers and caused by the ultrasonic waves reflected by the scatterers in the medium.

[0256] These steps S110 and S120 correspond to the above - mentioned generation and measurement steps.

[0257] The method S100 according to the present disclosure further includes the following steps:

[0258] - Step S130, in which a set of responses R from the medium is determined.

[0259] This step S130 is similar to the step of determining the focused reflection matrix R xx (z, f), but this step is different in that the responses are not necessarily organized into a matrix. They may only correspond to confocal signals, i.e., only to the diagonal of the focused reflection matrix. Moreover, these responses are at multiple frequencies f, at spatial positions r p =(x p , z p ) in a region (analysis region) around a reference point at spatial position r=(x, z). The responses at multiple frequencies around the reference point enable the analysis of this region around the reference point, which allows local correction of the aberrations and / or chromatic dispersion and / or reverberation caused by the medium upstream of the focal plane.

[0260] In particular, the method includes the following steps:

[0261] - Step S130, in which a set of responses R of the medium is determined , the responses being obtained by a focusing process for multiple frequencies f of the signals received from the reflected ultrasonic waves, for multiple points at spatial position r 0 =(x ui , z p ) in a region around a reference point at spatial position r=(x, z), based on the regular reflection matrix R p , z p ) and for a sound - speed model c.

[0262] The method S100 according to the present disclosure then includes Calibration process , which includes the following steps:

[0263] - Step S140, in which the response at different spatial position points (x,z) around the reference point is averaged or correlated through a medium to Determine the frequency calibration rule Φ, which is a frequency correction rule suitable for the reference point and is determined at frequency f;

[0264] - Step S180, in which for a plurality of frequencies f, by applying the frequency correction rule Φ to the response R of the medium around the reference point, to Determine the calibrated response of the medium R'.

[0265] Thus, the method advantageously enables local detection of the medium and focusing the reflection matrix for aberration correction, especially by determining the correction rule for each point of the spatial position r p of the medium M, for each frequency f of the ultrasonic waves.

[0266] Moreover, the method may include the following steps:

[0267] - Step S190, by combining the corrected responses at points with spatial position r = (x,z) at a plurality of frequencies f, to determine The intensity I of the point in the ultrasound image corresponding to the reference point c .

[0268] In particular, the intensity is calculated by the quadratic sum of the corrected responses R' at points with spatial position r = (x,z) for at least a part or all of the previously calculated frequencies f. For example, the calculation may be limited to a predetermined frequency band for ultrasonic characterization of the medium.

[0269] 1 - First embodiment

[0270] The first embodiment of the method S100 of the present disclosure performs matrix calculation by implementing a focusing process between an input virtual transducer point and an output virtual transducer point, and the input virtual transducer and the output transducer point may be different from or the same as each other. Thus, a large number of responses from the medium can be calculated, enabling very detailed analysis and specific correction processing to characterize the medium M.

[0271] Figure 11 The figure shows the method according to the first embodiment of the present disclosure. A transducer array arranged to face the medium is used for acoustic transmission and imaging of a region of the medium with a random "speckle" reflectivity.

[0272] In the first figure (A) of this Figure 11 , using a technique called beamforming or focusing, by directing towards a spatial position of and The focused emission in the directions of multiple foci successively emits multiple wavefronts into the medium. The wavefronts pass through the aberration layer and create focal spots around the foci that are distorted in space and time. Here, the temporal dispersion of the focused wavefronts is due to the frequency dependence of the sound speed in the aberration layer and / or multiple reflected echoes caused by the aberration layer.

[0273] In the next three diagrams (B) of this figure, during the return journey, the waves reflected at each focus pass through the aberration layer again, which causes the ultrasonic waves to be distorted in space and time again and further increases the echoes caused by multiple reflections in the aberration layer. The time signals received by the transducer are thus very complex and all include numerous echoes related to multiple reflections.

[0274] As shown in the fifth diagram (C), by averaging or correlating the echoes caused by multiple foci in the region around a reference point with a spatial position of r p a time response that would be generated, for example, by a virtual coherent reflector is obtained. This calculation is used to determine the spatio - frequency correction rule (shown here with the transducer basis u).

[0275] The sixth diagram (D) in this figure explains how the de - convolved virtual response can be used as the optimized delay rule to be applied to correctly focus on each focus, which compensates for the aberration and / or dispersion and / or multiple - reflection problems.

[0276] The method of this embodiment is described in more detail below:

[0277] Determine a set of responses R of the medium Step S130 of... includes determining the response obtained through the focusing process between a first point with a spatial position of r in =(x in , z) and a second point with a spatial position of r out =(x out , z), where the first point corresponds to the input virtual transducer, the second point corresponds to the output virtual transducer, the first and second points are in the region around the reference point and at the same depth z, and the lateral positions x in and x out form the focusing basis (x) at each depth z.

[0278] Then the response R is recorded in the focusing reflection matrix R xx (z, f), and the coefficients of this matrix can be written as R xx (z, f)=[R(x in , x out , z, f)].

[0279] Determine the frequency correction rule Φ (S140)

[0280] As Figure 6As shown Determine the frequency calibration rule Φ Step S140 then includes the following sub-steps performed at each depth z and each frequency f:

[0281] - Step S150, in which the forward projection of the focusing reflection matrix R xx (z, f) on the correction basis (c) is performed, Determine Determine the double reflection matrix R c (z, f);

[0282] - Step S160, in which, based on the double reflection matrix R c (z, f), Calculate the frequency calibration rule Φ, the frequency correction law is determined on the correction basis (c), Φ = [φ(c, f, r p )], such that the frequency correction law Φ is a space-frequency correction law;

[0283] - Step S170, in which the calibrated double reflection matrix around the reference point is determined R′ c , the coefficients of this matrix are written as R′ c (z, f) = [R′ c (x, c, z, f)], which are determined by performing the element-wise product between the double reflection matrix R c (z, f) and the phase conjugate of the frequency correction law Φ, i.e., by:

[0284]

[0285] where:

[0286] The symbol * refers to the phase conjugate operation;

[0287] The symbol is the Hadamard product, which gives:

[0288] R′ c (x, c, z, f) = R c (x, c, z, f)Φ * (x, c, z, f).

[0289] Step S180 of determining the calibrated response R' of the medium around the reference point Then it includes the back-projection of the corrected double reflection matrix R′ c (z, f) on the focusing basis (x) to determine the corrected focusing reflection matrix R′ xx (z, f).

[0290] Thus, the method advantageously enables local probing of the medium and the correction of the focusing reflection matrix for aberrations, in particular by determining the correction law for each point of the spatial position r p of the medium M and for each frequency f of the ultrasonic waves.

[0291] The correction is performed in a correction basis c suitable for the aberration to be corrected. The correction basis c is an input correction basis or an output correction basis.

[0292] Examples of correction bases include:

[0293] - Plane wave basis or spatial Fourier basis;

[0294] - Transducer basis u;

[0295] - A basis corresponding to the assumed position of the aberration body in the medium, i.e., a plane between, for example, the transducer plane (transducer basis u) and the focal plane (focal basis x);

[0296] - A basis corresponding to a plane determined by optimization, for example, a plane determined by a correlation matrix whose first eigenvalue is a maximum.

[0297] Double reflection matrix R c (z, f)(S150)

[0298] According to an embodiment of the method of the present disclosure, the forward projection S150 is used to determine the double reflection matrix R c (z, f). The forward projection of this step S150 can be performed by:

[0299] The matrix product between the transition matrix P and the focal reflection matrix R xx (z, f), i.e.:

[0300] R c (z, f) = P(z, f) × R xx (z, f)

[0301] where: P(z, f) = [P(c, x, z, f)] is the transition matrix between the focal basis (x) and the correction basis (c) at depth z for each frequency f.

[0302] The transition matrix P depends on the correction basis c used.

[0303] For a correction basis corresponding to the plane wave basis (c = k), the transition matrix P is a Fourier transform operator.

[0304] For the array 10 of the linear transducers 11 for generating a two-dimensional image, the coefficients of the transition matrix P can be written as:

[0305] P(k x , x, z, f) = P(k x , x) = exp(-ik x x)

[0306] where: k xis the transverse component of the wave vector k associated with each plane wave.

[0307] For the array 10 of matrix type transducers 11 for generating three-dimensional images, the coefficients of the transition matrix P can be written as:

[0308] P(k || , ρ, z, f) = P(k || , ρ, z, f) = exp(-ik || ·ρ)

[0309] where:

[0310] k || is the transverse component of the wave vector k associated with each plane wave; and

[0311] ρ = (x, y) is the transverse position vector.

[0312] For the correction basis corresponding to the transducer basis (c = u), the coefficients of the transition matrix P correspond to the normal derivative of the Green's function connecting each focus at spatial position (x, z) to each transducer 11 at spatial position (u, 0).

[0313] For the array 10 of linear transducers 11 for generating two-dimensional images, the coefficients of the transition matrix P can be written as:

[0314]

[0315] where: is the gradient projected along the depth direction z; and

[0316] G 2D (u, r) is the two-dimensional Green's function that links each transducer u = (u, 0) to each point r of the medium M:

[0317]

[0318] where: k 0 = 2πf / c 0 is the wave number;

[0319]

[0320] For the array 10 of matrix type transducers 11 for generating three-dimensional images, the coefficients of the transition matrix P can be written as:

[0321]

[0322] where: G 3D (u, r) is the Green's function that links each transducer u = (u x , u y, 0) The two-dimensional Green's function linked to each point r = (ρ, z) of the medium M:

[0323]

[0324] Therefore, the coefficients of the transition matrix P can be written as:

[0325]

[0326] Spatial-frequency correction rule (S160)

[0327] According to Figure 7 The first embodiment of step S160 for calculating the spatial-frequency correction rule, this calculation step S160 includes:

[0328] Step S161, in which based on the double reflection matrix R c (z, f) determined in the forward projection step S140, construct the correlation matrix C; and

[0329] Step S162, in which the correlation matrix C is analyzed , to determine the spatial-frequency correction rule Φ.

[0330] Spatial - frequency calibration rule For correcting the double reflection matrix R c (z, f) in the correction basis c to obtain the corrected double reflection matrix R' c (z, f). The correction performed in the correction basis c is then applied to the focusing basis x through the backprojection step S180 of the corrected double reflection matrix R' c (z, f) to obtain the corrected focusing reflection matrix R' xx (z, f).

[0331] Correlation matrix (S161)

[0332] According to the first variant of step S161, by the following calculation of the elements of the correlation matrix C = C cc , the correlation matrix C is determined in the frequency domain in the correction basis c:

[0333]

[0334] Where: R c is the double reflection matrix;

[0335] R ref is the model reflection matrix of the model medium, in which the sound speed is c 0 (the expected sound speed of the medium), and the plane reflector is located at the depth z;

[0336] x, z are the points r pThe coordinates of points in the surrounding area;

[0337] * is the conjugate operator.

[0338] The reflection matrix R in the previous equation c Add the reference matrix R ref Can compensate for the geometric components of the reflection matrix predicted by the c0 sound speed model and isolate the wavefront components of distortion and reverberation. When viewed from the focal plane, this operation is equivalent to virtually moving each focus (x, z) to the origin of the reference frame. For a medium with a random reflectivity (ultrasonic speckle) and a sufficient number of foci (x, z) in the area around the reference point at spatial position r p The obtained correlation matrix is equivalent to the correlation matrix measured for a coherent virtual reflector in the correction basis, where the range of the coherent virtual reflector is defined by the focal spot created by the focusing process in the area. This virtual reflector forms a guide star for determining the optimized spatio-temporal focusing rule in the area around the reference point r p surrounding.

[0339] According to the second embodiment of step S161, the correlation matrix C is determined in the frequency domain in the correction basis c by the following calculation of the elements of cc correlation matrix C:

[0340]

[0341] where: D c is the double distortion matrix obtained by:

[0342]

[0343] It can be expressed by the following formula in calculating this matrix:

[0344]

[0345] where: R c is the double reflection matrix;

[0346] R ref is the model reflection matrix of the model medium, in which the sound speed is c 0 (the expected sound speed of the medium), and the plane reflector is located at a depth z;

[0347] x, z are the coordinates of points in the area around the point r p surrounding;

[0348] * is the conjugate operator.

[0349] This second variant is equivalent to the first variant. It reveals the double distortion matrix D cIntermediate calculations. Thus it is not as straightforward as the first variant, but allows for the direct extraction of the focusing law through the simple singular value decomposition of matrix D c to directly extract the focusing law.

[0350] Compared with prior work that only considers the distortion matrix windowed in time at a single center frequency, the originality of the distortion matrix considered here lies in its frequency dependence, which enables the obtaining of complex spatio-temporal focusing laws that go beyond simple time-delay laws. In addition to conventional aberrations, this also helps to compensate for reverberation and frequency dispersion problems.

[0351] According to the third variant of step S161, the correlation matrix C is determined in the basis of the image points at the spatial position (x, z) by the following calculation of the elements of the correlation matrix C = C rr :

[0352]

[0353] where: R c is the double reflection matrix:

[0354] R ref is the model reflection matrix of the model medium, in which the sound speed is c 0 (the expected sound speed of the medium), and the planar reflector is located at a depth z;

[0355] x, z are the coordinates of the image points in the region around the point r p ;

[0356] * is the conjugate operator.

[0357] Adding the reference matrix R c to the double reflection matrix R ref in the previous equation can compensate for the geometric components of the reflection matrix predicted by the c 0 sound speed model and isolate the distortion and reverberation components of the wavefront. When observed from the focusing plane, this operation is equivalent to virtually moving each focus (x, z) to the origin of the reference frame, thereby forming a set of coherent guide stars, the amplitude distribution of which depends on the focal spot support and is also modulated due to the random reflectivity of the medium. Calculating the correlation matrix C in the focusing basis enables the determination of the phase shift between each pair of guide stars in the medium and will be used later to determine how the phase can be readjusted and different configurations combined to generate coherent guide stars.

[0358] According to the fourth variant of step S161, the correlation matrix C is determined in the basis of the image points at the spatial position (x, z) by the following calculation of the elements of the correlation matrix C = C rr :

[0359]

[0360] where: D c is a double distortion matrix obtained by the following formula:

[0361]

[0362] It can be expressed in terms of calculating this matrix as:

[0363]

[0364] where: R c is a double reflection matrix;

[0365] R ref is the model reflection matrix of the model medium, in which the speed of sound is c 0 (the expected speed of sound of the medium), and the plane reflector is located at a depth z;

[0366] x, z are the coordinates of the points in the region around the point r p ;

[0367] * is the conjugate operator.

[0368] Analysis (S162)

[0369] According to the first variant of step S162, by eigenvalue decomposition of the correlation matrix C, the analysis of the correlation matrix C is performed, and the space-frequency correction rule Φ is the first eigenvector U of the correlation matrix C in the correction basis (c) 1 , that is, C = C cc .

[0370] Since the correlation matrix is a Hermitian matrix its eigenvalues are positive real numbers.

[0371] Thus, the correlation matrix C cc can be written as:

[0372]

[0373] Or in terms of matrix coefficients, written as:

[0374]

[0375] where: U p corresponds to the eigenvector of the correlation matrix C;

[0376] σ p corresponds to the positive real eigenvalues of the correlation matrix C cc arranged in descending order: σ 1 > σ 2 > … > σ N .

[0377] Then the spatial-frequency correction law Φ is obtained, which is equal to the first eigenvector, i.e., Φ(r p ) = U 1 , or its normalized version: Φ(r p ) = exp(j arg{U 1})), i.e., the coefficients of the spatial-frequency correction law have the same unit amplitude, but its phase is equal to the phase of U 1 (the symbol arg{X} refers to the phase of the vector X); or inverse filter type correction, Φ(r p ) = exp(j arg{U 1}) / |U 1 |. If the signal-to-noise ratio is poor (appropriate filter), the first option is preferred. But generally, the second option will be preferred so that the correction does not act as an amplitude filter, but only corrects the phase distortion. Finally, when the aberration medium unevenly attenuates certain components and / or frequencies in the domain and one wants to enhance these components and / or frequencies to finally obtain a more accurate estimate of the reflectivity, the third option is relevant.

[0378] According to the second variant of step S162, the correlation matrix C is analyzed by the singular value decomposition of the double distortion matrix D c defined above in the second variant of step S161. The eigenvalue decomposition of the correlation matrix C cc in the first variant of step S162 is actually equivalent to the singular value decomposition (SVD) of the double distortion matrix D c when its coefficients are organized according to the following definition:

[0379] D c = [D c ({c, f}, {x, z})]

[0380] The singular value decomposition is applied to a rectangular matrix and is applied to the double distortion matrix D c , which is written as:

[0381]

[0382] Or in terms of matrix coefficients, written as:

[0383]

[0384] where: U p = [U p (c, f)] corresponds to the singular vectors of the double distortion matrix D c in the correction basis, or equivalently, corresponds to the eigenvectors of the matrix C cc defined in the first variant of step S162;

[0385] V p = [V p (x, z)] corresponds to the double distortion matrix D c of the singular vectors;

[0386] λ p corresponds to the double distortion matrix D c of the singular values, which by definition is equal to the eigenvalue σ of the correlation matrix C as defined in the first variant of step S162 p of the square root:

[0387] Thus, the space - frequency correction rule Φ is equal to the first singular vector of the double distortion matrix D c i.e., Φ(r p ) = U 1 ; or its normalized version: Φ(r p ) = exp(jarg{U 1}), i.e., the coefficients of this space - frequency correction rule have unit magnitude, but its phase is equal to the phase of U 1 (the symbol arg{X} denotes the phase of vector X); or inverse filter type correction, Φ(r p ) = exp(jarg{U 1}) / |U 1 |.

[0388] Compared with the eigenvalue decomposition of the correlation matrix C cc , the advantage of the singular value decomposition of the double distortion matrix D c lies in the computational speed of the numerical algorithm for singular value decomposition.

[0389] This search for the space - frequency correction rule Φ is also equivalent to, by using the following expression corresponding to iterative time - reversal calculation:

[0390] Φ n+1 (r p ) = C cc × Φ n (r p ),

[0391] where: Φ 0 is a random wavefront, e.g., Φ 0 = [1…1] T ,

[0392] to iteratively solve the following equation:

[0393] aΦ(r p ) = C cc × Φ(r p )

[0394] where: × is the matrix product, and a is a constant.

[0395] Then, through the following formula:

[0396]

[0397] or its normalized version:

[0398]

[0399] or its inverse filter version:

[0400]

[0401] the spatio - frequency correction rule Φ is obtained.

[0402] For n → ∞, the iterative time - reversal algorithm converges to the same first singular vector U c of the matrix D 1 . In practice, using the iterative time - reversal algorithm instead of SVD may have advantages because it converges after several iterations, resulting in faster calculations.

[0403] According to the third variant of step S162, by using the following expression corresponding to the iterative phase - reversal calculation:

[0404] Φ n+1 (r p ) = exp(j arg{C cc × Φ n (r p )})

[0405] where: × is the matrix product;

[0406] and, Φ 0 is a random wavefront, for example, Φ 0 = [1…1] T ,

[0407] the following equation is solved iteratively:

[0408] Φ(r p ) = exp(j arg{C cc × Φ(r p )})

[0409] to analyze the correlation matrix C cc ,

[0410] then through:

[0411]

[0412] or its inverse filter version:

[0413]

[0414] to obtain the spatio - frequency correction law Φ.

[0415] Compared with the previous scheme, the advantage of the iterative phase - reversal algorithm is that it estimates the phase of the correction law Φ(r p ) more reliably, and thus finally compensates better for the phase distortion caused by the aberration body.

[0416] According to the fourth variant of step S162, by using the following expression:

[0417] W n+1 (r p ) = exp(j arg{C rr ×W n (r p )})

[0418] where: × is the matrix product,

[0419] and, W 0 is a random wavefront, for example W 0 = [1…1] T ,

[0420] iteratively solve the following equation:

[0421] W = exp(j arg{C rr ×W})

[0422] to analyze the correlation matrix C rr ,

[0423] where, × is the matrix product,

[0424] this gives the following vector W:

[0425]

[0426] This vector W = [W(x, z)] defined in the focusing basis contains the phases of each incoherent guide star synthesized by focusing in the region around the reference point r p .

[0427] The phase conjugate of this vector W can then be used to readjust the phases of each incoherent virtual star so that they can be coherently recombined, thereby obtaining an estimate of the spatio - frequency correction law Φ that is not affected by the random reflectivity of the medium. Mathematically, this operation is written as:

[0428] φ(c, f, r p ) = exp{j×arg{∑ x，z D c (x, c, z, f)W* (x, c, z, f)}}

[0429] with respect to the distortion matrix D c the SVD (second variant of step S162) or the iterative phase inversion algorithm (third variant of step S162), the advantage of this approach being that it converges towards a calibration rule that is as equiplanar as possible, i.e. efficient for every point in the region around the reference point r p in the region around the reference point r

[0430] Spatial-frequency calibration rule (S160)

[0431] According to the second embodiment shown in step S160 of calculating the spatial-frequency calibration rule Figure 8 by maximizing the intensity or the confocal intensity of the ultrasound image in the region Ω p around the spatial position reference point r p in the optimization algorithm S163, this calculation step S160 is carried out

[0432] In other words, by the following maximization, the spatial-frequency calibration rule Φ is determined:

[0433]

[0434] where: is the intensity of the ultrasound image in the region Ω p for the spatial-frequency calibration rule φ(c, f, r p ) applied to the input and output

[0435] For example, by the triple sum of the frequency value f, the input calibration basis c in and the output calibration basis c out this ultrasound image is determined

[0436] Using the definitions provided above, the following calculations can be obtained, for example:

[0437]

[0438] The above calculations include the double reflection matrix R cc (c in , c out , f), which matrix is obtained by forward projection (step S150) into the input calibration basis c in and the output calibration basis c out and the coefficients thereof can be expressed as:

[0439]

[0440] A method for iteratively determining the spatio - frequency correction law Φ using ultrasonic image calculations. Even if the ultrasonic image is limited to the area around the spatial position reference point r p The optimization algorithm iteration may still be computationally time - consuming.

[0441] However, the advantage of this method is that it more accurately determines the spatio - frequency correction law Φ because it takes into account the reciprocity of the aberration correction to be applied to the outward and return paths of the ultrasonic wave.

[0442] The corrected double - reflection matrix (S170)

[0443] In the embodiment shown in Figure 6 of the method according to the present disclosure, by performing the term - by - term product between the double - reflection matrix R c (z, f) and the phase conjugate of the spatio - frequency correction law Φ, in step Determine the calibrated double reflection matrix in S170 R′ c (z, f)=[R′ c (x, c, f, z)], that is, by:

[0444]

[0445] where: the exponent * refers to the phase - conjugate operation;

[0446] The symbol is the Hadamard product, that is, the term - by - term matrix product of the matrix coefficients, which gives:

[0447] R′ c (x, c, f, z)=R c (x, c, f, z)Φ * (x, c, f, z).

[0448] The corrected reflection matrix (S180)

[0449] In an embodiment of the method according to the present disclosure, in step S180, by back - projecting the corrected double - reflection matrix R′ c (z, f) onto the focusing basis (x), R′ Calibrated focused reflection matrix R′ xx (z, f) is determined.

[0450] This back - projection is performed by the matrix product between the transition matrix P defined above and the focusing reflection matrix R′ c (z, f), that is:

[0451]

[0452] where: the exponent refers to the cross - conjugate matrix operation.

[0453] Correction process iteration (L1)

[0454] According to the method of the present disclosure Figure 5 In the illustrated embodiment, as Figure 5 shown in the L1 loop, iterate the correction process steps multiple times (two or more than two times), namely the following steps:

[0455] - Step S140 of determining the frequency correction rule Φ, which may include determining the double reflection matrix R c (z, f) in step S150, step S160 of calculating the frequency correction rule Φ, and step S170 of determining the corrected double reflection matrix R′ c (z, f); and

[0456] - Step S180 of determining the corrected focusing reflection matrix R′ xx (z, f).

[0457] In each iteration, the forward projection in step S150 uses the corrected focusing reflection matrix R′ xx (z, f) obtained during the back projection in step S180 of the previous iteration, rather than the focusing reflection matrix R xx (z, f).

[0458] Thus, in each iteration, the spatial-frequency correction rule is improved to better take into account one or more aberrations in the medium M.

[0459] According to a first variant of this iterative process, in each iteration of step S150 of determining the double reflection matrix R c (z, f), different correction bases c are used to correct, for example, different aberrations located at different positions in the medium.

[0460] For example, the medium M can be discretized or modeled by successive layers along the depth direction z, and the correction bases c for each iteration correspond to the planes of these successive layers. In other words, corrections corresponding to multiple aberrations in the medium M are applied during the iteration.

[0461] For example, spatially, the medium M can be segmented into regions either in a predetermined manner or automatically based on a first ultrasound image of the medium M. Each iteration of the iterative process will perform corrections on the correction bases c corresponding to each region of the medium M.

[0462] According to a second variant of this iterative process, in each iteration of step S150 of determining the double reflection matrix R c (z, f), either the input correction base projected onto the reflection matrix or the forward projection of the output correction base projected onto the reflection matrix is used. In the latter case, the matrix Rxx (z, f) is projected onto the corrected basis as follows:

[0463]

[0464] where: The symbol denotes the transpose matrix operation.

[0465] In the iterative sequence, it is possible to alternate between using the input corrected basis and the output corrected basis. Thus, in each iteration, the spatio-temporal correction rule Φ is improved, and the correction of aberrations is improved.

[0466] According to the third variant of this iterative process, in each iteration of step S150 for determining the double reflection matrix R c (z, f), the region around the spatial position point r p is used, and its size becomes smaller and smaller during the iteration. In other words, in each iteration, the size of the calculation region around the reference point r p is reduced. The size of this region can be used to refer to the width in the x direction, or the depth in the z direction, or both, or any other convention for the size of this region, such as being suitable for the scanning mode of the medium M.

[0467] Thus, the correction rule Φ becomes increasingly suitable for the aberrations near the reference point with spatial position r p .

[0468] Confocal image (S190)

[0469] According to an embodiment of the ultrasonic characterization method of the present disclosure, the method further includes the following steps:

[0470] - Step S190, in which, based on the integral of the diagonal coefficients of the corrected focusing reflection matrix R′ xx (z, f) over the bandwidth of the ultrasonic signal, that is, by combining the corrected responses R′ of the ultrasonic image points at spatial positions (x, z) at multiple frequencies f, Determine that of this point intensity Degree I c (x, z) is calculated. For example, the calculation is as follows:

[0471]

[0472] The intensities determined previously at multiple points are used to construct the corrected confocal image of the medium M, which corresponds to a conventional ultrasonic image without aberration, reverberation, and sound speed frequency dispersion problems in the medium under study.

[0473] 2 - Second embodiment

[0474] The second embodiment of the method S100 of the present disclosure performs a more direct calculation, for example in order to determine the characterization of the medium M mainly confocal. This simplified embodiment can be used to simply determine the intensity I of an ultrasonic image point c , in order to more quickly determine the ultrasonic image of the region of interest in the medium M.

[0475] In this second embodiment, only the focusing process between the same input virtual transducer points and output virtual transducer points is determined or calculated. This is equivalent to only determining the diagonal components of the focusing reflection matrix R xx (z,f) and recording them in the confocal reflection matrix R(z,f), which greatly reduces the number of medium responses calculated.

[0476] Figure 12 This method according to the second embodiment of the present disclosure is shown. A transducer array arranged to face the medium is used for acoustic transmission and imaging of a region of the medium having a random "speckle" reflectivity.

[0477] In this Figure 11 first figure (A), using a technique called beamforming or focusing, multiple waves are successively transmitted into the medium by focusing emissions in the directions of multiple foci having spatial positions of and . The waves pass through the reverberation layer and reach the foci, and the reverberation layer causes multiple reflected echoes.

[0478] In the next three figures (B) of this figure, on the return journey, the waves reflected at each focus pass through the reverberation layer again, resulting in the echoes being increased again due to multiple reflections in the reverberation layer. The time signals received by the transducer are very complex and all include numerous echoes related to multiple reflections.

[0479] As shown in the fifth figure (C), by averaging or correlating the echoes caused by multiple foci in the region around a reference point having a spatial position of r p , a time response that would be generated by a virtual coherent reflector, for example, is obtained. This calculation is used to determine the frequency correction rule to be applied to the signal to compensate for reverberation in the ultrasonic image.

[0480] The sixth figure (D) in this figure explains how the deconvolved time return virtual response can be used as the optimization delay rule to be applied in order to correctly focus on each focus, which compensates for the time dispersion and / or multiple reflection problems.

[0481] The method of this second embodiment is described in more detail below.

[0482] This second embodiment of the method S100 has the following characteristics.

[0483] Step S130 of determining a set of responses R of the medium including determining a response obtained through a focusing process between a first point at a spatial position r in =(x in , z) and a second point at a spatial position r out =(x out , z), where the first point corresponds to an input virtual transducer, the second point corresponds to an output virtual transducer, and the first and second points are identical (r in =r out ); and

[0484] recording all responses R in a confocal reflection matrix R, the coefficients of which can be written as: R = [R(x, z, f)].

[0485] Thus, compared to the first embodiment, the confocal reflection matrix R now has only a single lateral position parameter x, rather than two independent lateral position parameters x in and x out .

[0486] Then, the step of determining the frequency correction law Φ (S140) is directly performed by correlating the responses at respective spatial position points (x, z) around a reference point in the medium, and the coefficients of the frequency correction law are written as: Φ = [φ(f, r p )].

[0487] Subsequently By applying the frequency correction law to each frequency f, Determine the calibrated response around the reference point R′ Step (S180) is directly performed by taking the term product between the confocal reflection matrix R and the phase conjugate of the frequency correction law Φ, i.e.:

[0488]

[0489] where:

[0490] All corrected responses R′ are recorded in a corrected confocal reflection matrix R′, the coefficients of which can be written as R′ = [R′(x, z, f)].

[0491] The symbol is the Hadamard product, which gives:

[0492] R′(x, z, f) = R(x, z, f)φ * (f, r p ).

[0493] With the aid of this system, the method advantageously enables local probing of the medium and aberration correction of the confocal reflection matrix, in particular for a spatial position r of the medium M pFor a reference point, for multiple frequencies f of ultrasonic waves, a calibration rule is directly determined.

[0494] Compared with the first embodiment, by correlating the received responses, that is, without using the calibration basis and the double reflection matrix, the frequency calibration rule Φ is calculated more directly and simply.

[0495] Frequency calibration rule (S140)

[0496] According to step S140 of determining the frequency calibration rule in Figure 9 the first embodiment shown, this step S140 includes:

[0497] Step S141, in which based on the confocal reflection matrix R(z, f), Construct the correlation matrix C ; and

[0498] Step S142, in which it is analyzed the correlation matrix C to determine the frequency calibration rule Φ.

[0499] This step S140 of determining the frequency calibration rule is thus similar to step S140 of calculating the frequency calibration rule in the first embodiment. It does not require a double reflection matrix and is directly implemented on the confocal reflection matrix. Thereby, the calculation is simplified, but as described below, the possible variations of these method steps are also similar to those in the first embodiment.

[0500] According to the first variant of step S141 , the correlation matrix C is determined in the frequency domain by the following calculation:

[0501] C(f, f′) = ∑ x，z R(x, z, f)R * (x, z, f′)

[0502] where: R is the confocal reflection matrix;

[0503] x, z are the coordinates of the points in the region around the reference point;

[0504] * is the conjugate operator.

[0505] According to the Second variant, the correlation matrix C is determined in the basis of the image points at the spatial position (x, z) by the following formula:

[0506] C({x, z}, {x′, z′}) = ∑ f R(x, z, f)R * (x′, z′, f)

[0507] where: R is the confocal reflection matrix;

[0508] x, z are the coordinates of the image points in the region around the reference point;

[0509] * is a conjugate operator.

[0510] of step S141 , the analysis of the correlation matrix C is the According to the first variant of step S142 , and the frequency correction rule Φ is the first eigenvector U of the correlation matrix C 1 .

[0511] Eigenvalue decomposition , the analysis of the correlation matrix C is the solution of an equation involving the correlation matrix C and the frequency correction rule Φ, and similar to that described above in the first embodiment, this equation solution corresponds to iterative time reversal or iterative phase reversal.

[0512] According to the second embodiment of step S140 for determining the frequency correction rule, similar to that described above in the first embodiment, the frequency correction rule is determined by an optimization algorithm that maximizes the confocal intensity of the ultrasonic image in the region around the reference point.

[0513] Correction process iteration (L1)

[0514] As According to the second variant of step S142 shown, as Figure 5 shown in the L1 loop in

[0515] - Step S140 for determining the frequency correction rule Φ; and

[0516] - Step S180 for determining the corrected confocal reflection matrix R′(z, f).

[0517] In each iteration, step S140 for determining the frequency correction rule uses the corrected confocal reflection matrix R′(z, f) obtained during step S180 for determining the corrected response in the previous iteration, rather than the confocal reflection matrix R(z, f).

[0518] Thus, in each iteration, the frequency correction rule is improved to better take into account one or more aberrations in the medium M.

[0519] In each iteration of step S140, a spatial position reference point r with an increasingly smaller size can be used during these iterations p around the region. In other words, in each iteration, the size of the calculation region around the reference point rp is reduced. The size of this region can be used to refer to the width in the x direction, or the depth in the z direction, or both, or any other convention for the size of this region, such as suitable for the scanning mode of the medium M.

[0520] The correction rule Φ becomes increasingly suitable for the aberrations near the reference point with a spatial position of r p .

[0521] Confocal image (S190)

[0522] According to an embodiment of the ultrasonic characterization method of the present disclosure, the method further includes the following steps:

[0523] - Step S190 , in which the corrected responses R' of the ultrasonic image points at spatial positions (x, z) are combined at multiple frequencies f, Figure 5 I c (x, z). For example, the calculation is as follows:

[0524]

[0525] The intensities determined previously at multiple points are used to construct a corrected confocal image of the medium M, which corresponds to a conventional ultrasonic image without aberration, reverberation, and sound speed frequency dispersion problems in the medium under study.

[0526] Results on a calibrated experimental medium (referred to as a "phantom")

[0527] Determine the intensity of this point Shows a calibrated experimental medium, the so-called "phantom", which generates ultrasonic speckles of two cylindrical inclusions with a reflectivity higher than that of the surrounding medium, and has a plurality of reflectors (nylon wires) arranged along two lines (one horizontal and the other vertical). In the medium, the sound speed far from the point reflector is approximately 1,540 m / s -1 . An acrylic layer is arranged on the phantom, and this acrylic layer corresponds to an aberration layer with a sound speed of approximately 2,750 m / s -1 . Then, an array 10 of a transducer 11 of an ultrasonic probe is arranged on this layer. The medium is acoustically transmitted with multi-channel plane waves at multiple angles from -40 to +40 degrees relative to the transducer plane. The ultrasonic waves in this medium pass through the medium at different sound speeds and experience multiple reflections (reverberation) between the interfaces of these media, and then reach the phantom medium scatterers.

[0528] Figure 13 Shows an ultrasonic image obtained in the region of interest in the experimental medium shown in Figure 14 . This ultrasonic image is obtained assuming that the sound speed in the volume of the medium being acoustically transmitted is 1,540 m / s -1 . This ultrasonic image shows the aberration and reverberation problems caused by the acrylic layer on the image of the point reflector. First, the image of each point reflector is stretched laterally, showing that the lateral resolution of the image is due to the sound speed model c 0is reduced due to the difference from the actual sound speed distribution in the medium. Secondly, the image of each point reflector is behind the ballistic image of the reflector and repeats in the depth z direction, showing reverberation caused by multiple reflection echoes in the plexiglass layer. Therefore, the ultrasonic image is greatly disturbed and has lateral and axial distortions. For medical imaging, these problems result in low contrast, low resolution, and artifacts in the ultrasonic image that seriously hinder the practitioner's diagnosis.

[0529] Figure 13 shows the frequency correction law Φ obtained in the medium region B1 according to the method of the present disclosure. In the second embodiment, the frequency correction law is obtained by iterative phase inversion of the frequency correlation matrix C = [[C(f,f′)]] of the confocal signals from each point in the medium. The figure first shows the modulus and phase of the frequency correction law as a function of frequency, and then shows the time representation of the frequency correction law obtained by inverse Fourier transform. The frequency correction law phase corresponds to the phase shift to be applied to each frequency component of the ultrasonic signal. Equivalently, the time correction law corresponds to the time delay law that must be convolved with the received ultrasonic signal to (partially) compensate for the reverberation problem and obtain a more satisfactory ultrasonic image. Figure 15 The time representation of the frequency correction law includes a first large-amplitude echo corresponding to the error in the sound speed aspect of the propagation model in the medium, and multiple subsequent echoes corresponding to multiple reflections in the plexiglass layer. The echo time of the first large-amplitude echo is shorter than the ballistic time. The frequency correction law allows the aberration to be corrected on average in the B1 region of the experimental medium, which corresponds to the "speckle" region in the medium. The law is then used for the region of interest to characterize the medium and especially to improve the corrected ultrasonic image of the region of interest, as visible in Figure (b) below

[0530] The left figure (a) in the figure is an ultrasonic image of an experimental medium obtained using the focusing method commonly used in the prior art, without aberration correction, assuming the sound speed c Figure 14 is 1,540 m.s

[0531] Figure 16 shows the improvement obtained by the method according to the present disclosure.

[0532] in the figure. The left figure (a) in the figure is an ultrasonic image of an experimental medium obtained using the focusing method commonly used in the prior art, without aberration correction, assuming the sound speed c 0 is 1,540 m.s -1 The image has low resolution and poor quality.

[0533] The right figure (b) in this figure is an image obtained using the method for correcting the sound velocity and multiple reflection aberrations. It can be seen that the spatial resolution of the image of the point reflector in the medium is significantly improved, and the multiple reflection echoes are (partially) eliminated. The background image also shows a greater amplitude difference between the virtual reflector and the surrounding speckles. The reverberation compensation process thus significantly improves the quality of the ultrasonic image, especially the contrast.

[0534] As shown in Figure 16 figure (b) of Figure 16 shows the region of interest of the experimental medium of Figure 17 , after applying the aberration and reverberation compensation according to the first embodiment, that is, by applying the spatio-frequency correction rule Φ = [Φ(k x , f)], the obtained ultrasonic image. This rule is obtained by the iterative phase inversion of the spatio-frequency correction matrix C kk = C({k x , f}, {k′ x , f′}), which is obtained by the correlation at all points in the average field of view. As shown in Figure 13 figure (b) of Figure 16 , this global correction only partially compensates for the aberration and reverberation. Different from Figure 16 figure (b) of

[0535] Figure 16 , three regions C1, C2, and C3 are identified in this image to locally and iteratively apply the method of the present disclosure.

[0535] Figure 16 shows the results of determining the spatio-frequency correction rule Φ in the first region C1 of Figure 18 , in the region C2 of Figure 17 , and then in the region C3 of Figure 17 . Thus, in the first row (a) of Figure 17 , the spatio-frequency correction rule Φ in the first region C1 is shown by its spectrum (|C x ×Φ|) as a function of the transverse component k kk of the wave vector and the frequency f (left figure) and its phase (arg[Φ) as a function of k x and the frequency f (right figure). Figure 18 The second row (b) of Figure 18 shows the spatio-frequency correction rule Φ in the second region C2 in the same format. Another spatio-frequency correction rule Φ is also determined for the third region C3 in the ultrasonic speckles.

[0536] Figure 18 then shows the improvement obtained by using the local spatio-temporal correction rule including the corrections C1, C2, and C3 according to the method of the present disclosure.

[0537] The left figure (a) in the figure is obtained using the focusing method commonly used in the prior art, without aberration correction, assuming a sound speed c 0 of 1,540 m / s -1 for the experimental medium. This image has low resolution and is severely damaged due to the reverberation phenomenon.

[0538] The middle figure (b) in the figure is the image obtained by the method that enables global correction of the sound speed and multiple reflection aberrations in the entire region of interest of the image. This correction corrects the aberration evenly across the entire region of interest, which is already a step forward. This image is similar to Figure 19 that shown.

[0539] The right figure (c) in the figure is the image obtained by the method, through local iterative calculation of the spatial-frequency correction rule for the regions C1, C2, and C3 shown as Figure 17 shown. Figure 18 The correction of the spatial-frequency correction rule significantly improves the spatial resolution of the image and largely eliminates the echoes caused by multiple reflections in the plexiglass layer. Finally, the relative amplitude between the point reflector and the surrounding speckles is much better, which shows the contrast improvement provided by local compensation of aberrations and reverberation.

[0540] Ultrasonic characterization system

[0541] In Figure 17 Figure 3 is shown an ultrasonic characterization system 1 of a medium M according to the present disclosure. It includes:

[0542] - An array 10 of transducers 11, which is adapted to generate a series of incident ultrasonic waves in a region of the medium and measure, over time, the ultrasonic waves backscattered by the region; and

[0543] - A calculation unit 30, which is connected to the transducer array and is adapted to implement a method including the following steps:

[0544] - S110, generating, through the array 10 of transducers 11, a series of incident ultrasonic waves US in in a region of the medium, the series of incident ultrasonic waves being the transmission basis i; and

[0545] - S120, measuring a regular reflection matrix R ui (t) defined between the input transmission basis i and the output reception basis u, the coefficients of the regular reflection matrix corresponding to the signals received by the transducers and caused by the ultrasonic waves reflected in the medium;

[0546] - S130, determining a set of responses R of the medium, the responses being for a sound speed model c 0 based on the regular reflection matrix Rui (t), by performing, for multiple frequencies f, a focusing process for multiple points at spatial positions r p =(x p , z p ) around a reference point at spatial position r = (z, z);

[0547] - S140, determining a frequency correction rule Φ based on the response of the medium at each spatial position point (x, z), the frequency correction rule being suitable for the reference point and determined at frequency f;

[0548] - S180, determining a corrected response R′ of the medium by applying the frequency correction rule Φ to the response R of the medium around the reference point for multiple frequencies f.

Claims

1. A method for ultrasonic characterization of a medium (S100), the method comprising the following steps: - generating (S110) a series of incident ultrasonic waves (US) in a region of interest of the medium through an array (10) of transducers (11) in ), the series of incident ultrasonic waves are the transmitting basis (i); and - measuring (S120) the regular reflection matrix R defined between the input transmitting basis (i) and the output receiving basis (u) ui (t), the coefficients of which correspond to the signals received by the transducer and caused by the ultrasonic waves reflected in the medium, The method is characterized in that the method further comprises a correction process, which comprises the following steps: - Determine (S130) a set of responses R of the medium, which are for the sound velocity model c0, based on the regular reflection matrix R ui (t), by taking multiple frequencies f of the signal received from the reflected ultrasonic wave, for the p =(x p ,z p ) is obtained by focusing a plurality of points whose spatial positions are r=(x, z) in the region around the reference point; - determining (S140) a frequency correction law Φ, adapted to said reference point and determined at a frequency f, by averaging or correlating said responses of said medium at different spatial positions (x, z) around said reference point; - for said plurality of frequencies f, determining (S180) a corrected response R' of said medium by applying said frequency correction law Φ to said response R of said medium around said reference point.

2. The method according to claim 1, further comprising the steps of: - By setting the spatial position to r p The corrected responses of the reference point at multiple frequencies f are combined to determine (S190) the intensity I of the echo imaging image point corresponding to the reference point c .

3. The method of claim 1 or claim 2, wherein: - The determining (S130) of a set of responses R of the medium includes determining a set of responses R of the medium through a spatial position r in =(x in ,z) and the first point with spatial position r out =(x out , z), wherein the first point corresponds to the input virtual transducer, the second point corresponds to the output virtual transducer, and the first point and the second point are the same (r in =r out );as well as All responses R are recorded in the confocal reflectance matrix R, the coefficients of which can be written as R = [R(x,z,f)]; - performing said determination (S140) of a frequency correction law Φ by relating the responses of said medium at different spatial positions (x, z) around said reference point, the coefficients of which are written as Φ=[φ(f,r p )]; - said determining (S180) of the corrected response R′ around said reference point is performed by applying said frequency correction law at each frequency f, by performing a term product between said confocal reflection matrix R and the phase conjugate of said frequency correction law Φ, i.e. by: in: All said corrected responses R′ are recorded in a corrected confocal reflection matrix R′, the coefficients of which can be written as R′=[R′(x,z,f)]; The symbol is the Hadamard product such that: R′(x,z,f)=R(x,z,f)φ * (f,r p )。 4. The method according to claim 3, further comprising the steps of: - Determine (S190) the intensity I of the echo imaging image point at the spatial position (x, z) by combining the corrected responses R′ at multiple frequencies f c , that is, by: I c (x,z)=|∑ f R′(x,z,f)| 2 。 5. The method of claim 3 or claim 4, wherein: The determining of the frequency correction rule (S140) includes: Based on the confocal reflection matrix R(z,f), construct (S141) a correlation matrix C; and The correlation matrix C is analyzed (S142) to determine the frequency correction law v.

6. The method of claim 5, wherein: The correlation matrix C is determined in the frequency domain by the following formula: C(f,f′)=∑ x,z R(x,z,f)R * (x,z,f′) Wherein: R is the confocal reflection matrix; x,z are the coordinates of points in the area around the reference point; * is the conjugate operator.

7. The method of claim 5, wherein: The correlation matrix C is determined in a basis of image points at spatial locations (x, z) by the following formula: C({x,z},{x′,z′})=∑ f R(x,z,f)R * (x′,z′,f) Where: R is the confocal reflection matrix: x,z are the coordinates of the image point in the area around the reference point; * is the conjugate operator.

8. The method according to any one of claims 5 to 7, wherein: The analyzing (S142) of the correlation matrix C is the eigenvalue decomposition of the correlation matrix C, and the frequency correction rule Φ is the first eigenvector U1 of the correlation matrix C.

9. The method according to any one of claims 5 to 7, wherein: The analyzing (S142) of the correlation matrix C is a solution to an equation involving the correlation matrix C and the frequency correction law Φ, and the solution to the equation corresponds to iterative time inversion or iterative phase inversion.

10. The method according to any one of claims 3 to 9, wherein: The determining (S140) of the frequency correction law is performed by an optimization algorithm that maximizes the confocal intensity of the ultrasound image in the region around the reference point.

11. The method according to any one of claims 3 to 10, wherein: Iterate the correction process steps (S140, S180) multiple times; and In each iteration, the corrected confocal reflection matrix R'(z, f) obtained in the previous iteration is used instead of the confocal reflection matrix R(z, f).

12. The method of claim 11, wherein: In each iteration, reduce the spatial position to r p The size of the area around the reference point.

13. The method of any one of claims 3 to 12, wherein: The determining (S130) of the confocal reflection matrix R(z,f) comprises compensating for the temporal decay of the signal.

14. The method of claim 1 or claim 2, wherein: - The determining (S130) of a set of responses R of the medium includes determining a set of responses R of the medium through a spatial position r in =(x in ,z) and the first point with spatial position r out =(x out , z), wherein the first point corresponds to the input virtual transducer, the second point corresponds to the output virtual transducer, the first point and the second point are located at the same depth z in the region, and the lateral positions x of the first point and the second point are in and x out forming a focusing basis (x) at each depth z; and All responses R are recorded in the focus reflection matrix R xx (z,f), the coefficients of this matrix can be written as R xx (z,f)=[R(x in ,x out ,z,f)]; - The determining (S140) of the frequency correction law Φ comprises the following sub-steps implemented at each depth z and each frequency f: - through the focusing reflection matrix R xx Forward projection of (z,f) on the correction basis (c) determines (S150) the double reflection matrix R c (z,f); -Based on the dual reflection matrix R c (z, f), calculate (S160) the frequency correction law Φ, the frequency correction law is determined on the correction basis (c), Φ = [φ (c, f, r p )], so that the frequency correction law Φ is a space-frequency correction law; - Determine (S170) the corrected double reflection matrix R′ around the reference point c , the coefficients of this matrix are written as R′ c (z,f)=[R′ c (x,c,f,z)], which is obtained by performing the double reflection matrix R c The product of (z, f) and the phase conjugate of the frequency correction law Φ is determined by: in: The * symbol refers to the phase conjugation operation; The symbol is the Hadamard product such that: R′ c (x,c,f,z)=R c (x,c,f,z)Φ * (x,c,f,z) - said determining (S180) said medium's corrected response R' around a reference point comprises, by means of said corrected double reflection matrix R' c The back projection of (z,f) on the focusing basis (x) determines the corrected focusing reflection matrix R′ xx (z,f).

15. The method of claim 14, further comprising the steps of: - By focusing the ultrasound image point at the spatial position (x, z) at multiple frequencies f, the corrected focus reflection matrix R' xx The diagonal coefficients of (z, f) are combined to determine (S190) the intensity I of the point c ,Right now: I c (x,z)=|∑ f R′ xx (x,x,z,f)| 2 。 16. A method as claimed in claim 14 or claim 15, wherein: By combining the transition matrix and the focusing reflection matrix R xx The matrix product between (z, f) implements the forward projection (150), that is: R c (z,f)=P(z,f)×R xx (z,f) Where: P(z,f)=[P(c,x,z,f)] is the transition matrix between the focusing basis (x) and the correction basis (c) at depth z at each frequency f.

17. The method of any one of claims 14 to 16, wherein: The correction basis (c) is an input correction basis or an output correction basis.

18. The method of any one of claims 14 to 17, wherein: The calculation frequency correction rule (S160) includes: Based on the dual reflection matrix R c (z,f), construct (S161) a correlation matrix C; and The correlation matrix C is analyzed (S162) to determine the frequency correction law Φ.

19. The method of claim 18, wherein: The correlation matrix C is determined in the frequency domain in the correction basis (c) by the following formula: Where: R c is the dual reflection matrix; R ref is the model reflection matrix for a model medium in which the speed of sound is c0, the expected speed of sound for the medium, and the planar reflector is located at depth z; x,z is point r p The coordinates of points in the surrounding area; * is the conjugate operator.

20. The method of claim 18, wherein: The correlation matrix C is determined in the basis of image points at spatial positions (x, z) by the following formula: Where: R c is the dual reflection matrix; R ref is the model reflection matrix for a model medium in which the speed of sound is c0, the expected speed of sound for the medium, and the planar reflector is located at depth z; x,z is point r p The coordinates of the image points in the surrounding area; * is the conjugate operator.

21. The method of any one of claims 18 to 20, wherein: The analyzing (S162) of the correlation matrix C is the eigenvalue decomposition of the correlation matrix C, and the frequency correction rule Φ is the first eigenvector U1 of the correlation matrix C.

22. The method of any one of claims 18 to 20, wherein: The analyzing (S162) of the correlation matrix C is a solution to an equation involving the correlation matrix C and the frequency correction law Φ, and the solution to the equation corresponds to iterative time inversion or iterative phase inversion.

23. The method of any one of claims 14 to 22, wherein: The frequency correction law is calculated using an optimization algorithm that maximizes the confocal intensity of the ultrasound image in the region around the reference point (S160).

24. The method of any one of claims 14 to 23, wherein: Iterate the steps of the correction process (S140, S180) multiple times; and Wherein, in each iteration, the forward projection (S150) uses the corrected focus reflection matrix R' obtained during the backprojection (S180) of the previous iteration xx (z,f), instead of the focus reflection matrix R xx (z,f).

25. The method of claim 24, wherein: In each iteration of the forward projection (S150), one alternates between forward projection in an input correction basis and forward projection in an output correction basis.

26. The method of claim 24, wherein: In each iteration, the correction basis (c) of the forward projection (S150) is different.

27. The method of claim 24, wherein: In each iteration, reduce the spatial position to r p The size of the area around the reference point.

28. The method of any one of claims 14 to 27, wherein: The determining (S130) of the focusing reflection matrix R xx (z,f) includes compensating for the time decay of the signal.

29. A system (1) for ultrasonic characterization of a medium (M), the system comprising: - a transducer array (10) adapted to generate a series of incident ultrasonic waves in a region of interest of the medium and to measure over time the ultrasonic waves backscattered by the region of interest; as well as - a computing unit (30) connected to the transducer array and adapted to implement a method according to any one of claims 1 to 28.

Images