# Self-Portrait of the Focusing Process in Speckle: III. Tailoring Complex Spatio-Temporal Focusing Laws To Overcome Reverberations in Reflection Imaging

Elsa Giraudat, Flavien Bureau, William Lambert, Mathias Fink, Alexandre Aubry

<figure><img src="image_1.png" /></figure>

To cite this version:

Elsa Giraudat, Flavien Bureau, William Lambert, Mathias Fink, Alexandre Aubry. Self-Portrait of the Focusing Process in Speckle: III. Tailoring Complex Spatio-Temporal Focusing Laws To Overcome Reverberations in Reflection Imaging. 2026. hal-05497442

HAL Id: hal-05497442

https://hal.science/hal-05497442v1

Preprint submitted on 6 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Self-Portrait of the Focusing Process in Speckle:

III. Tailoring Complex Spatio-Temporal Focusing Laws To
Overcome Reverberations in Reflection Imaging

Elsa Giraudat,¹ Flavien Bureau,¹ William Lambert,² Mathias Fink,¹ and Alexandre Aubry*¹

¹Institut Langevin, ESPCI Paris, PSL University, CNRS, 75005 Paris, France

²SuperSonic Imagine, Aix-en-Provence, France

(Dated: February 6, 2026)
Abstract

This is the third article in a series of three dealing with the exploitation of speckle for imaging purposes. In complex media, a fundamental limit is the multiple scattering phenomenon that completely blurs the imaging process in depth. Matrix imaging can provide a relevant framework for solving this problem. As it proved to be an adequate tool for probing reverberations in speckle [E. Giraudat et al., Part I], we will show how it can be used to tailor complex spatio-temporal focusing laws to monitor the interference between the multiply-reflected paths and the ballistic component of the wave-field. To do so, we extend the distortion matrix concept to the frequency domain. An iterative phase reversal process operated from the space-time Fourier space is then used to compensate for reverberations and optimize both the axial and transverse resolution of the confocal image. Here, we first present an experimental proof-of-concept consisting in imaging a tissue-mimicking phantom through a reverberating plate before outlining the potential and the limits of this strategy for transcranial ultrasound and beyond.

Multiple scattering of waves in heterogeneous media is a fundamental problem with important applications, ranging from seismic wave propagation through the Earth's crust to deep tissue imaging in microscopy or trans-cranial ultrasound for brain imaging. Conventional focusing and imaging techniques based on Born approximation generally fail in strongly scattering media because of the multiple scattering events undergone by the incident wavefront. To compensate for those wave-front distortions induced by the medium heterogeneities, numerous adaptive focusing methods have been developed [1, 2]. These approaches generally consist in an application of time delays or phase shifts to each transmission and reception channel of the imaging device. While they are efficient, to a certain extent, to compensate for aberrations induced by forward multiple scattering trajectories [3, 4], they cannot overcome the problem of reverberations that require the tailoring of more complex spatio-temporal focusing laws. Each frequency component of the back-scattered echoes shall be actually compensated independently in order to make multiply-scattered waves associated with a broad time-of-flight distribution interfere constructively at any point inside the medium [5].

In this paper, we focus on the problem of reverberations in ultrasound which is particu-
larly prominent in transcranial imaging. Standard ultrasound actually relies on a focused
beamforming process in which only the ballistic component of the wave-field is used to
construct an image of the medium reflectivity (Fig. 1). This shortcoming can be the source
of multiple artifacts in complex media. Figure 1d illustrates this last statement by displaying
the ultrasound image obtained on a head phantom reproducing the acoustic characteristics of
an adult skull and containing a spherical target for evaluating the aberrations induced by the
skull layer. The reverberations generated at the skull-tissue interfaces induce multiple ghost
images of the bright scatterer. Furthermore, although they are not central to our discussion
below, it should also be noted that the multiple specular reflections at the skull interfaces
also seriously degrade image quality at shallow depths.

In the past, the problem of ultrasound focusing through the skull has been addressed in different ways. On the one hand, a first strategy consists in modeling the skull to estimate the time delays that need to be applied in order to correctly focus ultrasound. However, this requires prior knowledge of the geometry and acoustic properties of the skull, as well as the ability to predict wave propagation in this complex medium. The characteristics of the skull can be obtained in several ways, whether it be with X-rays [6–10], MRI [10–14] or using ultrasound data itself [15–19]. Coupling different imaging modalities can be quite tedious, as the success of the correction depends on the correct positioning of the probe. Moreover, these methods have been gradually refined to take into account refraction [20–23] and shear waves [24]. The more the approach is sophisticated, the more it is sensitive: a small error in estimating the thickness or mechanical properties of the skull can lead to a significant error in estimating the delay laws. Last but not least, a significant computation time is required for methods involving simulations [25, 26], which constitutes a barrier to the practical use of these methods [27, 28].

An alternative strategy has been to employ adaptive focusing techniques inspired from
astronomy. The broadband nature of ultrasound has led to the emergence of time reversal
processing to properly harness the multiply-scattered waves through the skull. However, such
approaches require either an active source inside the brain [12, 29] such as a miniaturized
transducer or the use of artificially induced ultrasonic stars such as a bubble that can
be generated by cavitation [30] or vaporization [31]. Such a process can then be coupled
with beam steering techniques to scan the focal spot around the initial source location [32].
However, it relies on an angular memory effect that turns to be very limited through the
<figure><img src="image_3.png" /><figcaption>FIG. 1. **Artifacts caused by reverberations in ultrasound.** a A medium is insonified by a plane wave. After passing through a reverberating layer, the plane wave is split into a tail of multiply-reflected waves. b The wave reflected by a strong diffuser in the medium produces a set of echoes, that are multiplied again when they pass back through the reverberating layer. c During the image reconstruction stage using a homogeneous velocity model, the multiple reflected echoes, associated with longer times-of-flight, are back-propagated to overestimated depths, giving rise to ghost images, or replicas, of the strong scatterer. Only a dispersive focusing law could correct for these artifacts. d, Example of a confocal image obtained on an ultrasound head phantom containing a bright target, appearing tripled due to reverberations induced by the skull.</figcaption></figure>

skull, especially for multiply-reflected waves. It also means that an active or passive source shall be used for each isoplanatic volume inside the region of interest. The constraints of in vivo imaging are such that a non-invasive reflection correction method that does not require the presence of a guide star is required.

To meet these criteria, it has been shown that spatial correlations of back-scattered echoes induced by a focused incident wave-front can be used to extract aberration laws [33]. To do so, an average over different incident focal spots is performed to synthesize a virtual guide star from different realizations of disorder. However, this experimental demonstration was restricted to the compensation of wave distortions and not to reverberations. Moreover, it
rely on an hypothesis of local isoplanicity that implies a repetition of the whole process for
each point in the field of view. It therefore requires a particularly time-consuming iterative
focusing process that makes this approach impractical for trans-cranial imaging.

To avoid this tedious focusing process, the reflection matrix approach is particularly
relevant [34, 35]. Once recorded, this matrix contains all the information available on the
medium and any iterative focusing process can now be performed in post-processing in
real time. In this paper, we will show how the time-focused reflection matrix and the
frequency-dependent distortion matrix introduced in the first paper of the series [36] can be
used to address each frequency component of the reflection matrix independently and thus go
beyond aberration correction involving a simple application of time delays to recorded signals.
In the first part, we propose to focus on the frequency degrees of freedom. A method for
compensating temporal dispersion phenomena is presented and illustrated using the academic
experiment involving an ultrasound phantom placed behind a reverberating plate [36]. In
the second part, the method is refined to address each angular component of the ultrasonic
field independently. Finally, in the last part, this method is applied to the more ambitious
case of a head phantom with a particularly heterogeneous skull bone. The potential and
limits of matrix imaging for trans-cranial applications are established before showing the
interest of the proposed method in other fields of wave physics such as optical microscopy
and seismology.

# I. DISPERSIVE CORRECTION OF REVERBERATIONS

In this work, we propose an approach to compensate for reverberations, particularly in
the difficult case of random scattering media encountered in medical ultrasound in absence
of any guide star. To begin with, we propose to focus on the frequency degrees of freedom of
the wave-field. In the first paper of the series [36], we saw that the temporal reflection matrix
can be used to synthesize a coherent virtual reflector in speckle. Here, we will show that this
virtual guide star can be used to tailor a dispersive focusing law in order to compensate for
reverberations.
<figure><img src="image_4.png" /><figcaption>FIG. 2. Extraction of a frequency dispersion law in speckle. **a** Acquisition of the reflection matrix **R**<sub>uθ</sub>(t): (i) each incident plane wave is multiply-reflected by the reverberating layer before being scattered by sub-resolved heterogeneities; (ii) the reflected waves also undergo reverberations before being recorded by the ultrasonic probe. **b** Confocal delay-and-sum beamforming process at each point **r** in the medium. The temporal dispersion of the echoes dispersed is captured by investigating the confocal signal I<sub>0</sub>(**r**, τ) as a function of time lapse τ with respect to the expected ballistic time. **c** An IPR analysis conducted on the correlation matrix between the backscattered echoes at each frequency leads to the synthesis of a coherent virtual source and the access to the associated reverberated signal H(τ) that accumulates the axial aberrations of the incident and reflected paths. **d** Time reversal of this signal yields a wave-front that can compensate for the axial aberrations of the reflection matrix, leading to a temporal recompression of the echoes associated with each point **r** of the medium.</figcaption></figure>

**A. Proof-of-concept experiment**

We are primarily interested in the academic experiment described in Ref. 36, in which a
plexiglas plate is inserted between the probe and the ultrasound phantom. The canonical
reflection matrix **R**<sub>uθ</sub>(t) of the medium is acquired using the sequence described in Table 1
of the same paper (Fig. 2a). The reverberating Plexiglas layer displays a much higher sound
velocity (c<sub>r</sub> ~ 2750 m.s<sup>-1</sup>) than the phantom (c = 1540 m.s<sup>-1</sup>). This academic experiment
mimics a transcranial imaging configuration. The Plexiglass layer induces strong defocusing and reverberation artifacts (Fig. 3a) similar to the ones produced by a skull (Fig. 1d).

As described in the first paper of the series [36], the time-focused de-scan reflection matrix $R_{\Delta r}(\tau)$ is constructed by a delay-and-sum beamforming process applied to $R_{u\theta}(t)$, with an homogeneous velocity model corresponding to the speed-of-sound $c$ in the phantom (Fig. 2b). Since we are first restricting the reverberation compensation to a time correction of temporal dispersion of echoes, the study of the focused reflection matrix is restricted to the confocal signal:

$$ \mathcal{I}_0(\mathbf{r}_{\text{in}}, \tau) = R(\Delta \mathbf{r} = 0, \mathbf{r}_{\text{in}}, \tau) \qquad (1) $$

This confocal signal is the result of a standard beamforming process with identical focusing points at input and output ($\mathbf{r}_{\text{in}} = \mathbf{r}_{\text{out}}$, Fig. 2b) but letting the echo time evolve beyond the expected ballistic time. The confocal signal taken at time lapse origin $\tau = 0$ corresponds to the raw confocal image displayed in Fig. 3b. It exhibits an amplitude drop with depth. This

<figure><img src="image_5.png" /><figcaption>FIG. 3. **Dispersive compensation of reverberations in speckle.** a Sketch of an incident focused wave-field in presence of a reverberating layer. The ballistic and multiply-reflected components focus at different depths shifted with respect to their respective isochronous volume. b Confocal image obtained for the homogeneous model $c_0 = 1540$ ms in the phantom - Plexiglas experiment. c Modulus and phase of the frequency dispersion law $\bar{H}(f)$ (Eq. 5) extracted by IPR in the speckle on the blue area in panel b. d Normalized real part of the temporal dispersion law $H(\tau)$ (Eq. 6), highlighting a first peak at $\tau_{0,0}^* \sim -12.4$ µs associated with the ballistic echo and secondary peaks associated with multiple reflections. e Corrected confocal image $\mathcal{I}_1(\mathbf{r}, \tau = 0)$ obtained by applying a dispersive focusing law (Eqs. 8 and 9). f Sketch of the dispersive focusing law tailored and applied in post-processing to make interfere constructively each multiply-reflected path with the ballistic component of the wave-field.</figcaption></figure>
attenuation is due to the decay of energy experienced by ultrasonic waves due to scattering and absorption. As shown in Appendix A, this phenomenon can be detrimental to the assessment of reverberations. A time gain compensation is thus applied to circumvent wave attenuation.

## B. Dispersive Phase Law

A reverberating layer implies that each frequency component of the wave-field undergoes a different phase aberration. A temporal Fourier transform of the confocal signal is thus performed in order to proceed with a frequency analysis of the wave-field:

$$ \bar{\mathcal{I}}_0(\mathbf{r}, f) = \int d\tau \mathcal{I}_0(\mathbf{r}, \tau) e^{-i2\pi f \tau}. \quad (2) $$

For lateral aberration correction [34], a broadband aberration phase law was extracted by considering the correlation between each spatial frequency component of the wave-fronts induced by virtual sources belonging to the same isoplanatic area. For axial aberration correction, a correlation between each temporal frequency component of the wave-fronts shall be considered instead. The difficulty here consists in finding an homogeneous speckle area that can provide the self-portrait of the wave-field. This area should be large enough in order to encompass a sufficiently high number of resolution cells to guarantee the convergence of the correlation matrix towards its covariance. However, it should avoid the presence of bright scatterers that would bias the estimation of a dispersive phase law. As a result, we here selected the area $A$ of virtual sources depicted by the blue rectangle in Fig. 3b whose lateral and axial dimensions are $w_x = 22$ mm by $w_z = 10$ mm, respectively.

The correlation between each frequency component of the wave-field in the selected area is stored in a spectral correlation matrix $C = [C(f, f')]$ whose coefficients are computed as follows:

$$ C(f, f') = \sum_{\mathbf{r} \in A} \bar{\mathcal{I}}_0(\mathbf{r}, f) \bar{\mathcal{I}}_0^*(\mathbf{r}, f'). \quad (3) $$

As described in the first paper of the series, an iterative phase reversal process can be applied to the correlation matrix to synthesize a coherent guide star from each virtual source in the patch $A$. The associated dispersive phase law is obtained using the following iterative process,

$$ \mathbf{W}^{(n+1)} = \exp \left( i \arg \left\{ \mathbf{C} \times \mathbf{W}^{(n)} \right\} \right), \quad (4) $$
with $\mathbf{W}^{(0)}$, an arbitrary wave-front that we will take arbitrarily as constant: $\mathbf{W}^{(0)} = (1 \cdots 1)^{\top}$. The phase $\phi_H$ of the invariant obtained at the end of this process, $\exp(i\phi_H) = \lim_{n \to \infty} \mathbf{W}^{(n)}$, is an estimator of the frequency dispersion law induced by the reverberating layer (Fig. 3c). As we will see further, this frequency dispersion law corresponds to the set of phase shifts to be applied to each of the frequency components of the echoes in order to make each multiple reflection path interfere constructively behind the reverberating plates. To check the consistency of this law over the frequency spectrum, a simple matrix product between $\mathbf{C}$ and $\mathbf{W}$ provides an estimator of the transfer function of the reverberating layer:

$$ \overline{\mathbf{H}} = \mathbf{C} \times \exp(i\phi_H). \qquad (5) $$

The amplitude of the resulting vector $\overline{\mathbf{H}}$ provides the frequency spectrum on which the IPR invariant is effective (Fig. 3c). Contrary to a SVD process that would result into a quasi-monochromatic phase law [36], the IPR process leads to a correction effective over the whole frequency bandwidth, here 3.5-4.5 MHz. The transfer function $\overline{\mathbf{H}}(f)$ also displays an amplitude modulation that can be related to the Fabry-Pérot resonances supported by the reverberating layer (Fig. 3c).

The inverse Fourier transform of $\overline{\mathbf{H}}$ gives access to the temporal impulse response associated with the reverberating layers (Fig. 2c):

$$ H(\tau) = \int df \overline{\mathbf{H}}(f) e^{i2\pi f\tau}. \qquad (6) $$

This temporal dependence of this impulse response $H(\tau)$ is displayed in Fig. 3d. It shows a main peak associated with the ballistic wave located at $\tau_{0,0}^* \sim -12.4 \ \mu s$. This time shift of the ballistic wave is the signature of the mismatch between the wave velocity distribution and the velocity model $c$ [36]. This impulse response also consists of secondary echoes of lesser intensity that can be attributed to multiple reflections induced by reverberating plate but also by the acoustic lens that lies at the surface of the probe [36]. They are the temporal counterparts of the Fabry-Pérot resonances highlighted by the transfer function $\overline{\mathbf{H}}(f)$ displayed in Fig. 3c.

The exact arrival times of these echoes observed in Fig. 3d can be determined analytically [36]. By considering the case of a double-layered medium, the iterative phase reversal process is shown to provide the reflected wave-field that would be obtained for a virtual reflector at the focusing point $\mathbf{r}^* = (x, z^{(0)})$ of the ballistic echo, with $z^{(0)} = z + L(1 - c_r/c)L$.
The arrival time of each multiply-reflected reflected echo depends on the number of reflections undergone by the wave on the way out (l) and on the way back (m):

$$ \tau_{l,m}^* = 2\frac{L}{c}\left(\frac{c_r}{c} - \frac{c}{c_r}\right) - \frac{2(l+m)L}{c_r} \quad (7) $$

The temporal dispersion law here corresponds to the accumulated reverberations on the forward and return paths. It is not possible here to separate the temporal dispersion experienced by the incident beam from that experienced by the reflected wave. The expected reverberation time, $\tau_r = 2L/c_r$ of the plexiglass plate ($c_r \sim 1000$ m.s$^{-1}$, L = 6 mm) and of the acoustic lens is of 4.3 µs and 2µs. Hence, we probably observe a superimposition of multiply-reflected echoes induced by each layer in Fig. 3.

### C. Time Compression

Because of the geometry of the reverberating plates and a wave velocity model $c_0$ matching with the phantom speed-of-sound, one can show that the aberration induced by the reverberation plate is spatially-invariant by translation. Thus, as a first approximation, the frequency dispersion law extracted in the speckle area depicted in Fig. 3a can be considered as isoplanatic. This is why, a global image correction can be performed using this law. To do so, the phase conjugate of this law is applied to the confocal signal $\bar{I}_0(\mathbf{r}, f)$ over the whole field-of-view (Fig. 2d), according to:

$$ \bar{I}_1(\mathbf{r}, f) = \bar{I}_0(\mathbf{r}, f) \exp[-i\phi_H(f)]. \quad (8) $$

A corrected confocal signal, $I_1(\mathbf{r}, \tau)$, is then recovered by an inverse Fourier transform:

$$ I_1(\mathbf{r}, \tau) = \int df \bar{I}_1(\mathbf{r}, f) e^{i2\pi fd\tau}. \quad (9) $$

Alternatively, the correction could have been performed in the time domain by convolving the reflection matrix with the time-reversed impulse response $H(\tau)$ (Fig. 2e).

The corrected confocal image, $I'(\mathbf{r}, \tau = 0)$, is shown in Fig. 3e. Its comparison with the initial image (Fig. 3b) shows the benefit of a dispersive phase correction. The targets now appear with much better contrast, comparable to that observed in a case without a reverberating layer [36]. The axial and transverse resolutions of their main lobe are also close to those obtained in an ideal case. This improvement in image quality can be assessed
quantitatively by plotting the transverse section of the images along lines of scatterers (Fig. 4d). It is made possible primarily by the time shift applied to the ballistic component, which adjusts the associated isochronous volume to its focusing depth $z^{*(0)}$. It confirms that even when considering a purely temporal correction, we are able to act on the transverse spread of the focal spot when this one is mainly a defocus. The image improvement is also induced by a time shift of multiple reflections to make them interfere constructively with the ballistic echo at the lapse time $\tau = 0$ (Fig. 3f).

However, we can see that this recombination is only partial and that axial secondary lobes associated with reverberations subsist and that their level is actually increased along the vertical line of scatterers (Fig. 4e). The reverberation plate not only induces an isoplanatic axial distortion of the wave-field. It also generates higher-order aberrations that requires a sharper tailoring of local spatio-temporal focusing laws. This is the aim of the next part.

## II. SPATIO-TEMPORAL CORRECTION OF REVERBERATIONS

Temporal and spatial aberration correction methods offer complementary approaches. In general, it may be necessary to combine them in order to shape a spatio-temporal focusing law adapted to the medium and correct the transverse and axial distortions of the image. With this in mind, we now propose to combine the frequency correction method presented in Sec. I with a transverse correction by generalizing the distortion matrix concept introduced in previous works [34, 35, 37].

### A. Frequency Analysis of Wave Distortions

To implement this, the time focused reflection matrix [36] is now considered at each depth $z$ of the medium: $\mathbf{R}_{xx}(z, \tau) = [R(x_{\text{out}}, x_{\text{in}}, z, \tau)]$. It contains the cross-talk between virtual sources $\mathbf{r}_{\text{in}} = (x_{\text{in}}, z)$ and the virtual detectors $\mathbf{r}_{\text{out}} = (x_{\text{out}}, z)$ located at the same depth. In the first paper of the series [36], the associated de-scan matrix, $\mathbf{R}_{\Delta x}(z, \tau) = [R(\Delta x, x_{\text{in}}, z, \tau)$ (with $\Delta x = x_{\text{out}} - x_{\text{in}}$), allowed us to simultaneously probe the temporal dispersion and spatial spread of the focal spots. Here, we show how to harness its spatio-temporal degrees of freedom in order to compensate for the residual distortions exhibited by the confocal image in Fig. 3e.
<figure><img src="image_6.png" /><figcaption>FIG. 4. **Local spatio-temporal correction of aberrations in the phantom-plexiglass experiment.** a Confocal image $I_0$ obtained for the homogeneous model $c_0 = 1540 \text{ m.s}^{-1}$ in the phantom - Plexiglas experiment. b Image $I_1$ (Eqs. 8-9) obtained after global frequency compensation of reverberations. c Image $I_2$ (Eqs. 14-15) after spatio-temporal refocusing of ultrasound data. d Transverse cross-section of each image at the depth corresponding to the horizontal line of scatterers in the phantom. e Longitudinal cross-section of each image at the lateral position corresponding to the vertical line of scatterers in the phantom. In panels d and e, the original image $I_0$ is displayed as a red continuous line while the corrected images $I_1$ and $I_2$ correspond to the yellow and green continuous lines, respectively.</figcaption></figure>

To do so, a correction basis shall be first determined. Ideally, this basis should maximize the spatial correlations between the distorted wave-fields. In the present case, the lateral invariance of the reverberating layer implies a maximal isoplanicity of wave distortions in
<figure><img src="image_7.png" /><figcaption>FIG. 5. Extraction of a spatio-temporal focusing law in speckle. a Acquisition of the reflection matrix $\mathbf{R}_{u\theta}(t)$: (i) each incident plane wave is multiply-reflected by the reverberating layer before being scattered by sub-resolved heterogeneities; (ii) the reflected waves also undergo reverberations before being recorded by the ultrasonic probe. b A delay-and-sum beamforming process decoupling the input and output focal spots is applied to $\mathbf{R}_{u\theta}(t)$ to compute a time focused reflection matrix $\mathbf{R}_{xx}(z, \tau)$ at each depth $z$. A de-scan operation followed by a 2D Fourier transform yields the frequency-dependent distortion matrix $\mathbf{D}_{kx}(f)$ in the plane wave basis. c A correlation matrix is computed between the distorted wave-fields associated with each virtual source $\mathbf{r}_{in}$ belonging to the selected area $\mathcal{A}$ centered around the point $\mathbf{r}_{\mathcal{A}}$. An IPR analysis conducted at each spatio-temporal frequency leads to the synthesis of a coherent virtual source and the estimation of the reverberated wave-field $\bar{T}(k_x, f)$ in the Fourier space, or equivalently, to $T(\Delta x, \tau)$ in real space. d Time reversal of this spatio-temporal wave-field yields a focusing law that can compensate for reverberations, leading to a spatio-temporal recompression of the echoes associated with each point $\mathbf{r}$ of the medium.</figcaption></figure>

the plane wave basis. A spatial Fourier transform is therefore applied to the de-scan matrix $\mathbf{R}_{\Delta x}(z, f)$ (Fig. 5b):

$$ \mathbf{D}_{kx}(z, f) = \mathbf{F}^{\top} \times \bar{\mathbf{R}}_{\Delta x}(z, f), \quad (10) $$
using the spatial Fourier transform operator F:

$$F(x, k_x) = \exp (ik_x x), \qquad (11)$$

where $k_x$ represents the transverse spatial frequency of the associated plane wave. The result of Eq. 10 is a distortion matrix whose columns contain the angular decomposition of the wave-front induced by each virtual source $\mathbf{r}_{\text{in}} = (x_{\text{in}}, z)$ (Fig. 5b). It contains the signature of the wave-front distortions undergone by each reflected wave-field when going through the reverberating layer.

In previous works [34, 38], this distortion matrix was considered at the single ballistic time and its singular value decomposition was used to extract an aberration phase law at the central frequency, or equivalently, a set of time delays in the temporal domain. Here, the distortion matrix is frequency-dependent and we will now show how this polychromatic feature can be used to extract a spatio-temporal focusing law. To do so, we extend to the spatial frequency domain [35] the iterative phase reversal process only operated in the temporal frequency domain in the first part of the manuscript. To do so, a correlation matrix, $\mathbf{C}_A = [C_A(\{k_x, f\}, \{k'_x, f'\})]$, shall be built between each spatio-temporal frequency component $(k_x, f)$ of the distorted wave-fields over a given area $A$ in the field-of-view:

$$C_A(\{k_x, f\}, \{k'_x, f'\}) = N_A^{-1} \sum_{\mathbf{r}_{\text{in}} \in A} D(k_x, f, \mathbf{r}_{\text{in}}) D^*(k'_x, f', \mathbf{r}_{\text{in}}) \exp \{-i[\phi_H(f) - \phi_H(f')]\}, \qquad (12)$$

where $N_A$ is the number of resolution cells contained in $A$. The choice of $A$ is subject to the following dilemma. On the one hand, it should encompass a sufficiently large number of resolution cells in speckle to guarantee the convergence of $\mathbf{C}_A$ towards its ensemble average, the covariance matrix $\langle \mathbf{C}_A \rangle$. On the other hand, it should be as small as possible to grasp high-order aberrations associated with a small isoplanatic length. To this end, the field of view has been divided into a set of spatial windows $A$ shown in Fig. 6a. Two types of areas can be distinguished: extended averaging areas in the speckle and more restricted areas on strong reflectors. Note that the phase term $\exp \{-i[\phi_H(f) - \phi_H(f')]\}$ in Eq. 12 enables a prior compensation of the frequency dispersion law determined at the previous step. This operation is made to facilitate the convergence of Eq. 12 towards $\langle \mathbf{C}_A \rangle$ by maximizing the correlation between each frequency component of the wave-field.
<figure><img src="image_8.png" /><figcaption>FIG. 6. Extraction of spatio-temporal focusing laws in the phantom through a reverberating layer. a The field of view is divided into a set of areas for estimating the aberration law, as shown in white. b-d The laws obtained for the three colored areas are shown as examples. For each, from left to right: amplitude and phase of the law $\overline{T}(k_x, f)$, and amplitude of its spatio-temporal Fourier transform.</figcaption></figure>

## B. Space-Time Transmittance

The transmittance $\overline{T}_A = [T_A(k_x, f)]$ of the reverberating layer seen from each point $\mathbf{r}$ in the field-of-view is extracted using the iterative phase reversal process (Eq. 4) applied to the spatio-temporal frequency correlation matrix $C_A$ computed over the corresponding patch $A$ (Eq. 12):

$$ \overline{\overline{T}}(\mathbf{r} \in A) = C_A \times \lim_{n \to \infty} \mathbf{W}^{(n)} \quad (13) $$

Figures 6b₁-d₁ and 6b₂-d₂ present the amplitude $|\overline{T}_A(k_x, f)|$ and phase $\phi_A(k_x, f)$, respectively, of the transmittance $\overline{T}_A(k_x, f)$ obtained over the different colored areas depicted in Fig. 6a. The associated focal spots, $T_A(\Delta x, t) = \int df \int dk_x \overline{T}_A(k_x, f)e^{ik_x x}e^{-i2\pi ft}$, are displayed in Figs. 6b₃-d₃. Figs. 6b₃ and d₃ display multiply-reflected echoes that were not perfectly
corrected by the global dispersive law $\phi_H(f)$ at the first step of the process. The focal spot displayed in Fig. 6c_3 shows a time shift compared to the ballistic time and can be explained by the subsistence of a defocus after the global dispersion correction. The variations observed between these focal spots (Figs 6b_3-d_3) and the corresponding phase laws (Figs 6b_2-d_2) highlight the spatial fluctuations of the wave-distortions across the field-of-view. This justifies a posteriori the need for a local estimation of the space-time transmittance. The reverberation phase laws also show transverse phase fluctuations that demonstrate the interest of going beyond a dispersive focusing law and the need for a spatio-temporal correction of aberrations. As previously revealed by the global frequency response of the reverberating layer (Fig. 3c), the frequency spectra $|\bar{T}(k_x, f)|$ of each aberration phase law (Figs. 6b_1-d_1) exhibit amplitude fluctuations that are reminiscent of the Fabry-Pérot resonances induced by the reverberating layer.

### C. Spatio-Temporal Focusing

A phase conjugation of the transmittance matrix estimator at each frequency can be used to tailor spatio-temporal focusing laws. The compensation of reverberations can then be performed as follows:

$$R'(\Delta x, r_{\text{in}} \in \mathcal{A}, \tau) = \int dk_x \int df \exp(-i[\phi_H(f) + \phi_{\mathcal{A}}(k_x, f)]) R(k_x, r_{\text{in}}, f) e^{ik_x \Delta x} e^{-i2\pi f \tau} \quad (14)$$

This operation is therefore equivalent to a phase reversal experiment with both a recompression of the confocal signal both in time and space (Fig. 5d). This is confirmed by the novel confocal image, $\mathcal{I}_2(\mathbf{r}_{\text{in}})$, that can be deduced from the confocal component ($\Delta x = 0$) of the corrected de-scan matrix at the ballistic time ($\tau = 0$):

$$\mathcal{I}_2(\mathbf{r}_{\text{in}}) = R'(\Delta x = 0, \mathbf{r}_{\text{in}}, \tau = 0) \quad (15)$$

This result shown in Fig. 4c is compared with the initial image (Fig. 4a) and the intermediate one (Fig. 4b) obtained after application of a global dispersive focusing law at the first step of correction. The additional local space-time correction provides a better contrast and a better axial resolution through a sharper recombination of multiple reflections, as highlighted by the arrows pointing on bright scatterers in Fig. 4. To be more quantitative, Fig. 4d displays a transverse cross-section of the confocal image along the horizontal line of scatterers in the
phantom. Compared to the initial image shown in red, the compensation of reverberations
leads to a remarkable contrast enhancement of 10 dB for the point-like scatterers throughout
the field-of-view. Fig. 4e displays the longitudinal cross-section of the confocal image along
the vertical line of scatterers in the phantom. Compared to the dispersive correction that
previously increased the level of secondary lobes along this line of scatterers (orange line in
Fig. 4e), the spatio-temporal correction restores and even slightly improves the initial contrast
while drastically improving the transverse resolution of the image.

Nevertheless, we can see that the recombination of multiple reflections is only partial and
that side lobe effects associated with reverberations remain in Fig. 4c. This is partly due to
the phase-matched nature of the filter chosen for the correction [39]. If we wanted to obtain
a perfect image, we would need to be able to whiten the frequency spectrum of the image
which is altered by the Fabry Pérot resonances of the reverberating plate (Fig. 6). This is
the principle of an inverse filter but this operation would be extremely sensitive to noise and
here hampered by out-of-focus echoes.

### III. APPLICATION TO TRANSCRANIAL IMAGING

Now that the distortion matrix concept has been extended to the frequency domain, the
method is applied to the more realistic case of a head phantom. In a previous study [35], matrix
imaging has been used to determine phase shift laws that compensate for the aberrations
suffered by direct echoes through the skull of the head phantom. Here we show how the
approach developed above can compensate for the temporal dispersion induced by multiple
reflections at the skull interfaces.

#### A. Experiment on an Adult Head Phantom

Spatio-temporal matrix imaging is now applied to an adult head ultrasound phantom. This
phantom (True Phantom Solutions) allows us to reproduce *in vitro* the conditions of brain
imaging in adults (Fig. 7a). It consists of a material that mimics the acoustic characteristics
of the human skull and an internal tissue that mimics those of biological soft tissues (Tab. I).
Highly reflective spherical inclusions made of bone-mimicking material are placed inside the
phantom. They are evenly spaced at 10 mm and arranged in a three-dimensional cross
pattern. Their diameter increases with depth: 0.2, 0.5, 1, 2, and 3 mm (Fig. 7a).

<table><caption>TABLE I. Acoustic properties of the head phantom.</caption><thead><tr><th>phantom area</th><th>speed-of-sound [m.s<sup>-1</sup>]</th><th>density [g cm<sup>-3</sup>]</th><th>attenuation at 2.25 MHz [dB cm<sup>-1</sup>]</th></tr></thead><tbody><tr><td>cortical bone</td><td>3000 ± 30</td><td>2.31</td><td>6.4 ± 0.3</td></tr><tr><td>trabecular bone</td><td>2800 ± 50</td><td>2.03</td><td>21 ± 2</td></tr><tr><td>cerebral tissues</td><td>1400 ± 10</td><td>0.99</td><td>1.0 ± 0.2</td></tr><tr><td>skin tissues</td><td>1400 ± 10</td><td>1.01</td><td>1.7 ± 0.2</td></tr></tbody></table>

The ultrasonic acquisition of the reflection matrix is performed as follows (see also Tab. II). A linear phased array probe is placed on the temporal window of the head phantom where the skull thickness is approximately 6 mm, as shown in Fig. 7a. Impedance matching was performed using ultrasonic contact gel. The reflection matrix was recorded in the plane wave basis to allow for a better signal-to-noise ratio [40]. A sequence of 81 plane waves is emitted, assuming a sound velocity of $c_0 = 1400 \text{ m.s}^{-1}$, in accordance with the expected velocity in the soft tissues of the phantom. The resulting confocal image is shown in Fig. 7b. It suffers from significant aberrations, as highlighted by the three targets in the field of view, which appear with significant reverberation echo tails.

## B. Local Space-Time Correction of Skull Reverberations

Following the approach described above, an iterative phase reversal process (Eq. 4) is applied to the correlation matrix of wave distortions expressed in the $(k, f)$-space (Eq. 12). Transmittance phase distribution $\overline{T}_A(k_x, f)$ are extracted in the colored areas shown in Fig. 7b and surrounding the bright scatterers. Even though the frequency bandwidth is limited to 1-3 MHz due to skull attenuation, the phase distribution $\phi_A(k_x, f)$ of the estimated transmittance is particularly complex whether it be along the temporal or spatial frequency dimension. This complexity can be assessed by computing the Strehl ratio $S$ that provides the ratio between the energy at focus with and without aberrations and reverberations:

$$S = \frac{\left| \langle \overline{T}(k_x, f) \rangle_{k_x, f} \right|^2}{\left| \langle \overline{T}(k_x, f) \rangle \right|^2_{k_x, f}} \qquad (16)$$
<figure><img src="image_9.png" /><figcaption>FIG. 7. **Experiment on head phantom.** a A linear phased array probe is placed at the temporal window of the head phantom. The blue beam indicates the field of view captured by the probe. b Original confocal image, under a velocity model $c_0 = 1400 \text{ m.s}^{-1}$. c Image obtained after local space-time correction using the signal from strong scatterers. d, e Phase distribution of the transmittance $\overline{T}(k_x, f)$ estimated over the blue and yellow areas displayed in panel b. The amplitude of $\overline{T}(k_x, f)$ is encoded with a transparency mask. f, g Associated spatio-temporal focal spots $T(\Delta x, \tau)$. h, i Depth dependence of the confocal signal $I(\mathbf{r}, \tau = 0)$ at the target location before (red line) and after reverberation compensation (green line).</figcaption></figure>

where the symbol $\langle \cdots \rangle$ stands for an average over the subscripts, here the spatial and temporal frequencies ($k_x, f$). The transmittances displayed in Figs. 7d and e display Strehl ratios $S$ of $7.6 \times 10^{-3}$ and $5.7 \times 10^{-3}$, respectively. These extremely low values illustrate the detrimental impact of skull heterogeneities on the focusing process. The corresponding focal spots $T_A(\Delta x, \tau)$ are also displayed in Figs. 7f and g. In each case, the focal spot displays a tail of multiply-reflected echoes lasting for more than $10 \ \mu\text{s}$ (Figs. 7f and g). This time dispersion is a manifestation of the strong reverberation phenomena taking place between the two skull interfaces.
<table><caption>TABLE II. Probe and acquisition sequence for recording the canonical reflection matrix $R_{u\theta}(t)$ in the head phantom experiment.</caption><thead><tr><th></th><th>parameter</th><th>value</th></tr></thead><tbody><tr><td rowspan="4">probe</td><td>type</td><td>linear, phased array</td></tr><tr><td>number of transducers</td><td>N<sub>u</sub> = 96</td></tr><tr><td>inter-element distance</td><td>δu = 0.2 mm ~ λ/2</td></tr><tr><td>aperture</td><td>Δu = 19.2 mm</td></tr><tr><td></td><td>central frequency</td><td>f<sub>c</sub> = 2.8 MHz</td></tr><tr><td rowspan="10">Acquisition</td><td>imaging platform</td><td>Aixplorer®, Supersonic Imagine</td></tr><tr><td>speed-of-sound model</td><td>c<sub>0</sub> = 1400 m s<sup>-1</sup></td></tr><tr><td rowspan="3">plane wave sequence</td><td>maximal angle</td><td>θ<sub>max</sub> = 40°</td></tr><tr><td>angular pitch</td><td>δθ = 1°</td></tr><tr><td>number</td><td>N<sub>θ</sub> = 81</td></tr><tr><td>emitted signal</td><td>two half-period pulse at f<sub>c</sub></td></tr><tr><td>sampling frequency</td><td>f<sub>s</sub> = 5 MHz</td></tr><tr><td>recording time</td><td>Δt = 217 µs</td></tr></tbody></table>

As before, the estimated reverberation phase laws (Figs. 7d and e) are used to correct the confocal image (Eqs. 14-15). The resulting image is displayed in Fig. 7c. These complex focusing laws compensate for most of the reverberations that initially polluted the confocal image (Fig. 7b). The new confocal image thus constitutes a more accurate estimator of the brain phantom reflectivity. As shown by the confocal signal at the target location displayed in Fig. 7h and i, the image of the bright scatterers shows a drastic gain of in terms of contrast (~10 dB at $z = 56$ mm, see Fig. 7i).

### C. Limit in the Speckle Regime

Even though the difference between Figs. 7b and 7g is spectacular, we have not been able to extract a spatio-temporal focusing law in any speckle area of the brain for several reasons:

* the extremely weak signal-to-noise ratio due to skull attenuation

* the limited size of the purely speckle areas not affected by bright scatterers in the
space-time ranges $(\Delta x, \tau)$

* the strong anisoplanetism of transcranial aberrations due to the curved shape of the skull layer, its heterogeneous composition, and its irregular thickness.

The latter aspect can be illustrated by comparing the phase distributions of the skull transmittance seen from two different locations inside the brain phantom (Figs. 7d and e). The two transmittances actually exhibit a weak similarity index of 0.13. This limited isoplanicity is critical in speckle because it imposes a limit on the size of the correlation area $A$. This is a fundamental problem because our transmittance estimator suffers from a bias $|\delta\bar{T}(k_x, f)|^2$ scaling as:

$$|\delta\bar{T}(k_x, f)|^2 \propto \frac{N_T}{N_{in}}. \qquad (17)$$

where $N_T \sim S^{-1}$ is the number of coherence grains exhibited by the transmittance $\bar{T}$. As the reverberation laws are quite complex in the present case ($N_T \sim 200$), the required number $N_{in}$ of focusing points is way too large to meet an even partial condition of isoplanicity. This is why we were unable to observe convergence of our method towards a satisfying focusing law in speckle in the present case.

## IV. DISCUSSION

Even though this last experiment has outlined some of its limits for trans-cranial imaging, the proposed method constitutes a breakthrough in the field of ultrasound imaging. Admittedly, it can be seen as the digital counterpart of the various adaptive focusing methods already developed in the past [33, 41]. However, to our knowledge, these approaches had never been applied to the compensation of reverberations. Furthermore, they were particularly difficult to implement in practice because they required the iterative time reversal process to be repeated many times on all the speckle grains in the area under consideration. On the contrary, matrix formalism allows for compensation of aberrations and reverberations in post-processing. This feature offers much greater flexibility in terms of correction basis. It also allows an iteration of the adaptive focusing process, which is particularly critical when trying to correct for such high-order aberrations as reverberations. In fact, the number of iterations required to converge on the focusing laws is no longer limited by their physical
execution, and the adaptive focusing process can be repeated for any correction plane and
focal point of interest.

With regards to the isoplanatic limit, a transition to 3D imaging could help us to address
a larger number of speckle grains in the same isoplanatic area. A previous study on this
head phantom has actually shown that isoplanicity is greater in the (x, y) plane than in the
(x, z) plane [35]. Moreover, brain imaging is widely performed nowadays with contrast agents
(microbubbles) that are exploited to beat diffraction in ultrasound localization microscopy [42].
Interestingly, microbubbles can also be used as ultrasound guide stars for adaptive focusing.
While current aberration correction methods have been restricted so far to the compensation
of wave distortions undergone by the ballistic wave-front [43], the frequency-dependent
distortion matrix approach can be exploited for tailoring spatio-temporal focusing laws
that would enable the constructive interference of multiply-reflected paths with the ballistic
component of the wave-field. At last, speckle dynamics provided by blood flow can also be of
interest to multiply the realization of disorder for speckle [44] and ensure the convergence of
the iterative phase reversal process despite a limited isoplanatic area.

Beyond trans-cranial imaging, the extension of matrix imaging to the temporal degrees of
freedom can be fruitful for the inspection of other organs in ultrasound in which reverberation
is a long-standing problem. For instance, the abdominal layer induce strong multiple
reflections that often prevent imaging the liver at shallow depth [45]. Beyond ultrasound, the
proposed approach can be applied to other field of wave physics in which the multi-element
technology allows a reflection matrix to be measured. In seismology, multiples are known to
often pollute the migration process and have been tackled by Marchenko redatuming [46] in
multi-layered configurations. Matrix imaging can, nevertheless be a more robust and adaptive
tool than the existing approaches that require a precise knowledge of the seismic wave velocity
inside the medium [47, 48]. Complex areas such as volcanoes or fault zones are typical areas
in which conventional migration or more sophisticated tools such as full waveform inversion
have failed so far. Extending the distortion matrix concept to the frequency domain would
lead to the design of spatio-temporal focusing laws that would allow a better monitoring of
magma storage areas or fluid injection in depth.

At an opposite scale, adaptive optics has been limited so far to a spatial correction of
aberrations in microscopy [1]. Yet, light propagation in tissues is distorted by forward
multiple scattering events [49]. A frequency-dependent analysis would be thus extremely
rewarding for a proper realignment of these complex trajectories. While such a process is
illusory experimentally with conventional adaptive optics, matrix imaging offers an adequate
framework to perform this task in post-processing [3, 4]. Several studies have actually shown
how a multi-wavelength and multi-angle illumination scheme coupled to an interferometric
measurement of the reflected wave-field can lead to a measurement of a multi-spectral
distortion matrix [50–52]. As demonstrated in this study, this quantity is actually the ideal
starting point for the design of spatio-temporal focusing laws able to control the interference
of multiple scattering paths in depth.

To conclude, in this work, we have tackled the issue of reverberations in ultrasound
imaging. The aim of this work was to propose a correction method that broadens the range
of aberrations that can be compensated for. The spatial and temporal degrees of freedom
offered by broadband transducer arrays have been exploited using a matrix formalism to
probe the spatial and temporal dispersion of ultrasonic echoes and then to find the focusing
law adapted to the medium heterogeneities. Extending the distortion matrix concept to the
frequency domain and taking full advantage of the iterative phase reversal process in the
(f, k)-space, the method was successfully applied to compensate for defocus and multiple
reflections induced by a reverberating layer on a tissue-mimicking phantom in an academic
experiment. The interest of this method has then been explored for trans-cranial imaging
through an experiment on a head phantom. While the approach has shown its ability for
designing adaptive spatio-temporal focusing laws in presence of bright scatterers, it has also
shown some limits in speckle areas due to the strong heterogeneities of the skull that limit
the range of the memory effect. Nevertheless, the proposed method paves the way towards
numerous applications ranging from ultrasound to optical microscopy or seismology in which
deep imaging is long-standing goal and reverberations a fundamental problem.
**Acknowledgments.** The authors are grateful for the funding provided by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement 819261, REMINISCENCE project, AA). W.L. acknowledges financial support from the SuperSonic Imagine company.

**Author contributions.** A.A. and M.F. initiated the project. A.A. supervised the project. F.B. and W.L. performed the phantom experiment. E.G. developed the post-processing tools. E.G. and A.A. performed the theoretical analysis. E.G. and A.A. prepared the manuscript. E.G., F.B., W.L., M.F. and A.A. discussed the results and contributed to finalizing the manuscript.

## Appendix A: Time gain compensation

<figure><img src="image_10.png" /><figcaption>FIG. 8. **Wave attenuation compensation.** a Raw confocal image. b Linear fit of the natural logarithm of the average confocal intensity (besides bright scatterers) as a function of depth leading to an estimation of the attenuation coefficient $\beta$. c,c Estimated impulse response $H(\Delta x, \tau)$ obtained in the speckle area before (c) and after (d) attenuation compensation, respectively.</figcaption></figure>

In a statistically homogeneous disordered medium, the mean confocal intensity, $\langle |\mathcal{I}_0(x, z)|^2 \rangle$, shall scale as $\exp(-2z/\ell_{\text{ext}})$, with $\ell_{\text{ext}}$, the extinction length. To retrieve such an exponential
decay, the speckle-like fluctuations of the confocal image (Fig. 8a) should be priorly smoothed out by lateral averaging. The resulting mean confocal intensity is displayed in log-scale as a function of effective depth z (Fig. 8b). Its evolution in speckle (i.e outside of bright scatterers) can be fitted by an exponential curve (red dashed line in Fig. 8b) whose decay provides an estimation of $\ell_{ext} \simeq 117$ mm. The overall fitting curve can be used to normalize at each time the raw data such that:

$$ \mathbf{R}'_{u\theta}(t) = \mathbf{R}_{u\theta}(t) \exp [c_0 t/(2\ell_{ext})] \quad (A1) $$

Figure 8c and d compares the impulse response $H(\Delta x, \tau)$ associated with the reverberating layer without and with the prior time gain compensation of Eq. A1. In absence of time gain compensation, the amplitude of echoes decreases with echo time. It gives rise to a strong bias on our estimation of the impulse response (Fig. 8c) with a shift towards negative time lapse of the ballistic echo time $\tau_{0,0}^*$ and the underestimation of multiply-reflected echoes that arrive at larger times-of-flight. After time gain compensation (Fig. 8d), the direct echo is retrieved at the expected ballistic time and multiple reflections are satisfactorily grasped by the estimated impulse response.

References.

[1] M. J. Booth, Adaptive optical microscopy: the ongoing quest for a perfect image, Light Sci. Appl. 3, e165-e165 (2014).

[2] R. Ali, T. Brevett, L. Zhuang, H. Bendjador, A. S. Podkowa, S. S. Hsieh, W. Simson, S. J. Sanabria, C. D. Herickhoff, and J. J. Dahl, Aberration correction in diagnostic ultrasound: A review of the prior field and current directions, Z. Med. Phys., 267 (2023).

[3] S. Kang, Y. Kwon, H. Lee, S. Kim, J. H. Hong, S. Yoon, and W. Choi, Tracing multiple scattering trajectories for deep optical imaging in scattering media, Nat. Commun. 14, 6871 (2023).

[4] U. Najar, V. Barolle, P. Balondrade, M. Fink, A. C. Boccara, M. Fink, and A. Aubry, Harnessing forward multiple scattering for optical imaging deep inside an opaque medium, Nat. Commun. 15, 7349 (2024).

[5] A. P. Mosk, A. Lagendijk, G. Lerosey, and M. Fink, Controlling waves in space and time for imaging and focusing in complex media, Nat. Photonics 6, 283 (2012).

[6] J.-F. Aubry, M. Tanter, M. Pernot, J.-L. Thomas, and M. Fink, Experimental demonstration
of noninvasive transskull adaptive focusing based on prior computed tomography scans, J. Acoust. Soc. Am. **113**, 84 (2003).

[7] S. Pichardo, V. W. Sin, and K. Hynynen, Multi-frequency characterization of the speed of sound and attenuation coefficient for longitudinal transmission of freshly excised human skulls, Phys. Med. Biol. **56**, 219 (2011).

[8] U. Vyas, P. Ghanouni, C. H. Halpern, J. Elias, and K. B. Pauly, Predicting variation in subject thermal response during transcranial magnetic resonance guided focused ultrasound surgery: Comparison in seventeen subject datasets, Med. Phys. **43**, 5170 (2016).

[9] S. A. Leung, T. D. Webb, R. R. Bitton, P. Ghanouni, and K. Butts Pauly, A rapid beam simulation framework for transcranial focused ultrasound, Sci. Rep. **9**, 7965 (2019).

[10] T. D. Webb, S. A. Leung, P. Ghanouni, J. J. Dahl, N. J. Pelc, and K. B. Pauly, Acoustic Attenuation: Multifrequency Measurement and Relationship to CT and MR Imaging, IEEE Trans. Ultras. Ferroel. Freq. Cont. **68**, 1532 (2021).

[11] J. Sun and K. Hynynen, Focusing of therapeutic ultrasound through a human skull: A numerical study, J. Acoust. Soc. Am. **104**, 1705 (1998).

[12] K. Hynynen and F. A. Jolesz, Demonstration of Potential Noninvasive Ultrasound Brain Therapy Through an Intact Skull, Ultras. Med. Biol. **24**, 275 (1998).

[13] P. Su, S. Guo, S. Roys, F. Maier, H. Bhat, E. Melhem, D. Gandhi, R. Gullapalli, and J. Zhuo, Transcranial MR Imaging-Guided Focused Ultrasound Interventions Using Deep Learning Synthesized CT, Am. J. Neuroradiol. **41**, 1841 (2020).

[14] S. A. Leung, D. Moore, Y. Gilbo, J. Snell, T. D. Webb, C. H. Meyer, G. W. Miller, P. Ghanouni, and K. Butts Pauly, Comparison between MR and CT imaging used to correct for skull-induced phase aberrations during transcranial focused ultrasound, Sci. Rep. **12**, 13407 (2022).

[15] A. Wydra, E. Malyarenko, K. Shapoori, and R. G. Maev, Development of a practical ultrasonic approach for simultaneous measurement of the thickness and the sound speed in human skull bones: A laboratory phantom study, Phys. Med. Biol. **58**, 1083 (2013).

[16] M. Hajian, R. Gaspar, and R. G. Maev, Accurate 3-D Profile Extraction of Skull Bone Using an Ultrasound Matrix Array, IEEE Trans. Biomed. Eng. **64**, 2858 (2017).

[17] H. Nguyen Minh, J. Du, and K. Raum, Estimation of Thickness and Speed of Sound in Cortical Bone Using Multifocus Pulse-Echo Ultrasound, IEEE Trans. Ultrason. Ferroel. Freq. Cont. **67**, 568 (2020).
[18] G. Renaud, P. Clouzet, D. Cassereau, and M. Talmant, Measuring anisotropy of elastic wave velocity with ultrasound imaging and an autofocus method: Application to cortical bone, Phys. Med. Biol. **65**, 235016 (2020).

[19] B. Du, J. Wang, H. Zheng, C. Xiao, S. Fang, M. Lu, and R. Mao, A novel transcranial ultrasound imaging method with diverging wave transmission and deep learning approach, Comput. Methods Programs Biomed. **186**, 105308 (2020).

[20] B. D. Lindsey and S. W. Smith, Refraction Correction in 3D Transcranial Ultrasound Imaging, Ultrason. Imag. **36**, 35 (2014).

[21] K. Shapoori, J. Sadler, A. Wydra, E. V. Malyarenko, A. N. Sinclair, and R. G. Maev, An Ultrasonic-Adaptive Beamforming Method and Its Application for Trans-skull Imaging of Certain Types of Head Injuries; Part I: Transmission Mode, IEEE Trans. Biomed. Eng. **62**, 1253 (2015).

[22] C. Jiang, Y. Li, B. Li, C. Liu, F. Xu, K. Xu, and D. Ta, Ray Theory-Based Transcranial Phase Correction for Intracranial Imaging: A Phantom Study, IEEE Access **7**, 163013 (2019).

[23] M. Mozaffarzadeh, E. Verschuur, M. D. Verweij, V. Daeichin, N. De Jong, and G. Renaud, Refraction-Corrected Transcranial Ultrasound Imaging Through the Human Temporal Window Using a Single Probe, IEEE Trans. Ultrason. Ferroelec. Freq. Cont. **69**, 1191 (2022).

[24] G. T. Clement, P. J. White, and K. Hynynen, Enhanced ultrasound transmission through the human skull using shear mode conversion, J. Acoust. Soc. Am. **115**, 1356 (2004).

[25] T. Wang and Y. Jing, Transcranial ultrasound imaging with speed of sound-based phase correction: a numerical study, Phys. Med. Biol. **58**, 6663 (2013).

[26] M. Mozaffarzadeh, C. Minonzio, N. de Jong, M. D. Verweij, S. Hemm, and V. Daeichin, Lamb Waves and Adaptive Beamforming for Aberration Correction in Medical Ultrasound Imaging, IEEE Trans. Ultrason. Ferroelec. Freq. Cont. **68**, 84 (2021).

[27] J.-F. Aubry, O. Bates, C. Boehm, K. Butts Pauly, D. Christensen, C. Cueto, P. Gélat, L. Guasch, J. Jaros, Y. Jing, R. Jones, N. Li, P. Marty, H. Montanaro, E. Neufeld, S. Pichardo, G. Pinton, A. Pulkkinen, A. Stanziola, A. Thielscher, B. Treeby, and E. Van 'T Wout, Benchmark problems for transcranial ultrasound simulation: Intercomparison of compressional wave models, J. Acoust. Soc. Am. **152**, 1003 (2022).

[28] C. Angla, B. Larrat, J.-L. Gennisson, and S. Chatillon, Transcranial ultrasound simulations: A review, Med. Phys. **50**, 1051 (2023).
[29] J.-L. Thomas and M. Fink, Ultrasonic beam focusing through tissue inhomogeneities with a time reversal mirror: Application to transskull therapy, IEEE Trans. Ultrason. Ferroel. Freq. Cont. 43, 1122 (1996).

[30] M. Pernot, G. Montaldo, M. Tanter, and M. Fink, "Ultrasonic stars" for time-reversal focusing using induced cavitation bubbles, Appl. Phys. Lett., 4 (2014).

[31] K. J. Haworth, J. B. Fowlkes, P. L. Carson, and O. D. Kripfgans, Towards aberration correction of transcranial ultrasound using acoustic droplet vaporization, Ultrasound Med. Biol. 34, 435 (2008).

[32] G. Clement and K. Hynynen, Micro-receiver guided transcranial beam steering, IEEE Trans. Ultrason. Ferroelec. Freq. Cont. 49, 447 (2002).

[33] G. Montaldo, M. Tanter, and M. Fink, Time Reversal of Speckle Noise, Phys. Rev. Lett. 106, 054301 (2011).

[34] W. Lambert, L. A. Cobus, T. Frappart, M. Fink, and A. Aubry, Distortion matrix approach for ultrasound imaging of random scattering media, Proc. Nat. Acad. Sci. USA 117, 14645 (2020).

[35] F. Bureau, J. Robin, A. Le Ber, W. Lambert, M. Fink, and A. Aubry, Three-dimensional ultrasound matrix imaging, Nat. Commun. 14, 6793 (2023).

[36] E. Giraudat, F. Bureau, W. Lambert, M. Fink, and A. Aubry, Self-portrait of the focusing process in speckle: I. spatio-temporal imaging of wave packets in complex media, arXiv:2602.05907 (2026).

[37] A. Badon, V. Barolle, K. Irsch, A. C. Boccara, M. Fink, and A. Aubry, Distortion matrix concept for deep optical imaging in scattering media, Sci. Adv. 6, aay7170 (2020).

[38] W. Lambert, L. A. Cobus, J. Robin, M. Fink, and A. Aubry, Ultrasound matrix imaging—Part II: The distortion matrix for aberration correction over multiple isoplanatic patches, IEEE Trans. Med. Imag. 41, 3921 (2022).

[39] E. Herrin and T. Goforth, Phase-matched filters: Application to the study of rayleigh waves, Bull. Seismol. Soc. Am. 67, 1259 (1977).

[40] G. Montaldo, M. Tanter, J. Bercoff, N. Benech, and M. Fink, Coherent plane-wave compounding for very high frame rate ultrasonography and transient elastography, IEEE Trans. Ultras. Ferroel. Freq. Cont. 56, 489 (2009).

[41] S. Flax and M. O'Donnell, Phase-aberration correction using signals from point reflectors and
diffuse scatterers: basic principles, IEEE Trans. Ultras. Ferroel. Freq. Cont. **35**, 758 (1988).

[42] C. Errico, J. Pierre, S. Pezet, Y. Desailly, Z. Lenkei, O. Couture, and M. Tanter, Ultrafast ultrasound localization microscopy for deep super-resolution vascular imaging, *Nature* **527**, 499 (2015).

[43] J. Robin, C. Demené, B. Heiles, V. Blanvillain, L. Puke, F. Perren, and M. Tanter, In vivo adaptive focusing for clinical contrast-enhanced transcranial ultrasound imaging in human, *Phys. Med. Biol.* **68**, 025019 (2023).

[44] B.-F. Osmanski, G. Montaldo, M. Tanter, and M. Fink, Aberration correction by time reversal of moving speckle noise, *IEEE Tran. Ultrason. Ferroel. Freq. Cont.* **59**, 1575 (2012).

[45] G. Ferraioli, R. G. Barr, A. Berzigotti, I. Sporea, V. W. Wong, T. Reiberger, T. Karlas, M. Thiele, A. C. Cardoso, O. T. Ayonrinde, L. Castera, C. F. Dietrich, H. Iijima, D. H. Lee, W. Kemp, C. P. Oliveira, and S. K. Sarin, Wfumb guidelines/guidance on liver multiparametric ultrasound. part 2: Guidance on liver fat quantification, *Ultrasound Med. Biol.* **50**, 1088 (2024).

[46] K. Wapenaar, J. Thorbecke, J. van der Neut, F. Broggini, E. Slob, and R. Snieder, Marchenko imaging, *GEOPHYSICS* **79**, WA39–WA57 (2014).

[47] R. Touma, T. Blondel, A. Derode, M. Campillo, and A. Aubry, A distortion matrix framework for high-resolution passive seismic 3-D imaging: application to the San Jacinto fault zone, *California, Geophys. J. Int.* **226**, 780 (2021).

[48] E. Giraudat, A. Burtin, A. L. Ber, M. Fink, J.-C. Komorowski, and A. Aubry, Matrix imaging as a tool for high-resolution monitoring of deep volcanic plumbing systems with seismic noise, *Commun. Earth Environ.* **5**, 509 (2024).

[49] C. Dunsby and P. M. W. French, Techniques for depth-resolved imaging through turbid media including coherence-gated imaging, *J. Phys. D: Appl. Phys.* **36**, R207–R227 (2003).

[50] P. Balondrade, V. Barolle, N. Guigui, E. Auriant, N. Rougier, C. Boccara, M. Fink, and A. Aubry, Multi-Spectral Reflection Matrix for Ultra-Fast 3D Label-Free Microscopy, *Nat. Photon.* **18**, 1097 (2024).

[51] Y. Zhang, M. Dinh, Z. Wang, T. Zhang, T. Chen, and C. W. Hsu, Deep imaging inside scattering media through virtual spatiotemporal wavefront shaping (2023), arXiv:2306.08793.

[52] Y.-R. Lee, D.-Y. Kim, Y. Jo, M. Kim, and W. Choi, Exploiting volumetric wave correlation for enhanced depth imaging in scattering medium, *Nat. Commun.* **14**, 1878 (2023).