Pulse-echo ultrasound attenuation tomography

Objective. We present the first fully two-dimensional attenuation imaging technique developed for pulse-echo ultrasound systems. Unlike state-of-the-art techniques, which use line-by-line acquisitions, our method uses steered emissions to constrain attenuation values at each location with multiple crossing wave paths, essential to resolve the spatial variations of this tissue property. Approach. At every location, we compute normalized cross-correlations between the beamformed images that are obtained from emissions at different steering angles. We demonstrate that their log-amplitudes provide the changes between attenuation-induced amplitude losses undergone by the different incident waves. This allows us to formulate a linear tomographic problem, which we efficiently solve via a Tikhonov-regularized least-squares approach. Main results. The performance of our tomography technique is first validated in numerical examples and then experimentally demonstrated in custom-made tissue-mimicking phantoms with inclusions of varying size, echogenicity, and attenuation. We show that this technique is particularly good at resolving lateral variations in tissue attenuation and remains accurate in media with varying echogenicity. Significance. Based on a similar principle, this method can be easily combined with computed ultrasound tomography in echo mode for speed-of-sound imaging, paving the way towards a multi-modal ultrasound tomography framework characterizing multiple acoustic tissue properties simultaneously.


Introduction
Ultrasound waves propagating through tissues undergo energy losses mainly due to absorption and scattering (Szabo 2014).This attenuation significantly differs between tissue types and can reveal disease-related changes in tissue composition and structure (Duck 1990).For example, the controlled attenuation parameter (CAP), which estimates the average ultrasound attenuation in the liver, is a widely used noninvasive tool to assess diffuse hepatic steatosis, i.e., fat accumulation in the liver (Eddowes et al. 2019, Cao et al. 2022, An et al. 2022).Despite being a powerful diagnostic technique, CAP is designed to provide a single value corresponding to a small volume within the liver (Sasso et al. 2010) and cannot thus characterize the inherent tissue heterogeneity.This hinders its extension to clinical applications involving focal lesions or other organs, such as the breast, where attenuation has proven relevant for differentiating malignant and benign lesions (Calderon et al. 1976, D'Astous & Foster 1986, Nam et al. 2013).
Imaging techniques that quantify the spatial distribution of ultrasound attenuation in tissue can offer a versatile tool to overcome the limitations of CAP.Several approaches have been suggested using either time-domain (Ghoshal & Oelze 2012) or frequencydomain formulations (Kanayama et al. 2013, Pawlicki & O'Brien 2013, Coila & Lavarello 2018, Brandner et al. 2021, Kim & Varghese 2008).Although computationally more demanding, the latter are more common as they provide efficient ways to compensate for system-and transducer-dependent effects.Frequency-domain techniques retrieve the attenuation in tissue by analyzing the spectral information of backscattered radiofrequency signals.They are classified into three types depending on the information they exploit (Labyed & Bigelow 2011), which can be either frequency-dependent amplitude variations (spectral difference and spectral-log difference methods), the downshift of the center frequency (spectral shift methods), or a combination of the two (hybrid methods).Independent of the type or domain, nearly all these techniques consider conventional line-by-line scanning to relate data variations with depth to attenuation.That is, they use simplified one-dimensional formulations where ultrasound waves propagate only in the axial direction, limiting the amount of data used to constrain the attenuation at each location and, thus, limiting the spatial resolution and accuracy of reconstructed images (Labyed & Bigelow 2011, Samimi & Varghese 2017).
To account for the two-dimensional (2-D) nature of attenuation in tissue, Coila & Lavarello (2018) reformulated the spectral-log difference technique as a regularized inverse problem.Instead of estimating the attenuation at each depth independently, their method includes 2-D spatial constraints via regularization.They show that this can substantially improve the trade-off between the spatial resolution and accuracy of state-of-the-art techniques; however, regularization strategies rely on subjective choices rather than actual 2-D data constraints, which would require interrogating each tissue location with ultrasound waves propagating along multiple directions.An attempt in this direction was suggested by Treece et al. (2007) to correct shadowing and enhancement artifacts in B-mode images.They developed a method that estimates the lateral variations in the cumulative attenuation from pairs of beams steered at opposite angles.Although successful at restoring B-mode images, their technique cannot characterize the local attenuation in tissue since this information is not actually required for artifact correction.
In this work, we present a fundamentally new 2-D approach to quantify the spatial distribution of tissue attenuation using pulse-echo ultrasound systems.We suggest using a set of steered emissions in order to measure the changes between the amplitudes of echoes that are detected using different steering angles.We can relate these to the attenuation-induced amplitude losses undergone by the incident waves, allowing us to formulate a tomographic technique to reconstruct tissue attenuation.Computed ultrasound tomography in echo mode (CUTE) uses a similar approach to relate tissue speed of sound to the echo phase shifts observed when probing tissue at different angles (Jaeger et al. 2015, Stähli, Kuriakose, Frenz & Jaeger 2020, Jaeger et al. 2022) and has demonstrated unprecedented spatial and contrast resolution in both tissuemimicking phantoms (Stähli, Frenz & Jaeger 2020) and in-vivo (Stähli et al. 2023).In the following, we first introduce the theoretical foundations of the proposed attenuation imaging technique in section 2 and discuss its practical implementation in section 3. We then validate this approach using numerical wave propagation simulations in section 4 and finally demonstrate its experimental feasibility in section 5 using calibrated tissuemimicking phantoms with varying acoustic properties.

Beamformed ultrasound images
In pulse-echo ultrasound, we use a transducer array to emit pressure waves through the tissue and acquire the backscattered wavefield at each transducer element.If we neglect multiple scattering and denote χ(x) the scattering function of the medium defined in Ω ⊂ R 2 , we can generally express the backscattered wavefield sensed at x r in the temporal-frequency domain as where ω = 2πf is the angular frequency of the waves, p e (x, ω) refers to the emitted pressure field, and G(x r , x, ω) is the Green's function, i.e., the impulse response of the medium between locations x and x r (Jensen 1991).In this study, we consider insonifying the tissue with steered plane waves, i.e., where A denotes the amplitude of the plane wave, and the wave number vector k defines its propagation direction.Since tissue is attenuating, this vector is considered complex, i.e., k = k r + jk i .The real part k r = ω/c(ω) encodes the wave phase velocity c, whereas the imaginary part k i = α(ω) describes the ultrasound attenuation (in units of Np/m) exhibiting the frequency power law relation α(ω) = α 0 ω y in tissues (Szabo 2014).Here α 0 is the attenuation coefficient, and y is the power law exponent typically ranging from 1 to 2. In the following, we omit the dependency on ω for a condensed notation.
We can generate an ultrasound image I(r) per plane-wave emission by applying delay-and-sum beamforming on recorded signals p(x r ).This process typically assumes a lossless homogeneous medium (indicated with the subscript 0) to focus the received energy into an arbitrary focal point r via the operation where is the receive point-spread function (Lambert et al. 2020).Here * indicates the complex conjugate operation, and the summation is performed over all receiving elements.Note that k 0 is related to the beamforming process and is thus a real vector.Although we typically integrate the equation (3) over the frequency range of ultrasound pulses to obtain high-resolution images, the next subsection considers the monochromatic formulation for simplicity.

Relationship between attenuation and cross-correlations of beamformed images
Assume we emit two plane waves steered at directions k 1 and k 2 sequentially and reconstruct the corresponding images I 1 (r) and I 2 (r), respectively.As a first approximation, we consider tissue as a diffuse scattering medium, meaning that scatterers are spatially uncorrelated and satisfy ⟨χ * (x 2 )χ(x 1 )⟩ = ⟨|χ| 2 ⟩δ(x 1 − x 2 ), where ⟨•⟩ denotes an ensemble average and δ is the Dirac distribution (Mallart & Fink 1991, Lambert et al. 2020).Then, the ensemble-averaged cross-correlation between the two images C 12 (r) := ⟨I * 2 (r)I 1 (r)⟩ is given by A i is the amplitude of the emitted plane wave along k i direction.
The point-spread functions in equation (3) restrict the range of locations significantly contributing to any arbitrary image position r 0 .Thus, without loss of generality, we can focus our analysis of C 12 to a subvolume Ω 0 ⊂ Ω centered around r 0 , where x = r 0 + ξ ∈ Ω 0 .In this case, C 12 at r 0 reduces to When the difference between the plane-wave steering angles ∆φ is sufficiently small to satisfy where the loss per wavelength is assumed small, i.e., k i ≪ k r (Szabo 1994), the exponential inside the integral in equation ( 6) approximates to 1, and the crosscorrelation becomes Thus, in a first approximation, when plane-wave emissions satisfy equation ( 7), only the propagation time differences related to the incident waves contribute to the phase of the cross-correlation.This approximation led to the first development of CUTE for estimating the speed-of-sound distribution in tissue (Jaeger et al. 2015), and we propose a similar approach for attenuation imaging in this work.In order to isolate the attenuation-induced amplitude terms from equation ( 8), we suggest dividing C 12 with the auto-correlation Then, the log-amplitude of the normalized cross-correlation is where we assumed A 1 ≈ A 2 following the requirement of small steering angle differences previously discussed.Thus, the normalized cross-correlation operation is equivalent to placing a virtual receiver at r 0 that measures the changes between attenuation-induced amplitude losses undergone by the incident waves.Equation (11) establishes the forward problem of the proposed attenuation imaging technique in this work.Note that if we use C 22 for the normalization instead, only a sign change occurs in equation ( 11); so, we can alternatively express the forward problem as which may be useful to reduce the observational noise via the averaging of both logamplitude terms.The validity of this forward problem is subject to condition (7), which depends on the size of point-spread functions in equation (3) and, thus, on the focusing quality of beamformed images.For instance, for steering angle differences of ∆φ = 2.5 • (see next section), we would require a ||ξ|| approximately smaller than 4λ, which is considerably larger than the size of diffraction-limited point-spread functions arising when aberrations are not present (Chen et al. 2020).This condition may not hold when the assumed and true impulse responses in equation (3) differ substantially.In this case, the normalization operation cannot accurately isolate attenuation-induced amplitude losses of the incident waves, and we may expect an increased noise level in the log-amplitude measurements.We will discuss this point further in section 5.1.

Inversion
In attenuation imaging, it is common to consider ultrasound waves propagating as straight rays.Under this simplification, the term on the right-hand side of equation ( 11) becomes where L i (r 0 ) refers to the ray path connecting the transducer array and position r 0 along direction k i , and dl is the differential arc length along this ray.Upon inserting equation ( 13) into equation ( 11) or ( 12), we can express the forward problem in matrix notation as where d ∈ R N is a vector containing normalized cross-correlation log-amplitude measurements for each tissue location and pair of incident waves, m ∈ R M contains the attenuation values α at each location of the tissue (discretized here using a rectilinear grid), and F ∈ R N ×M is the forward operator whose rows contain differences between pairs of ray paths for each measurement location.Due to the limited data coverage, the forward operator F is ill-conditioned, and the inverse problem requires regularization.To ensure close-to-real-time solutions critical to medical ultrasound, we use first-order Tikohonov regularization, for which the inversion has the closed-form solution Here, the superscript T stands for the matrix transpose operation, D = [D x D z ] T contains the first-order finite-difference operators along the lateral (x) and axial (z) directions penalizing non-smooth solutions, and λ is the regularization parameter, which we optimize using the L-curve method (Hansen & O'Leary 1993).Note that equation (15) assumes normally distributed and uncorrelated noise, equal for all observations.As observed in pulse-echo speed-of-sound tomography, reconstructions can benefit from anisotropic regularization, particularly in layered media (Sanabria et al. 2018, Stähli, Kuriakose, Frenz & Jaeger 2020).This can be implemented by reformulating λD T D as The L-curve approach cannot optimize more than one regularization parameter without additional constraints.Thus, we keep a fixed λ x /λ z ratio to apply this approach for anisotropic regularization.

Uncertainty assessment
Uncertainties of attenuation estimates m in equation ( 15) are described by means of the posterior covariance matrix Γ post := F T F + λD T D −1 .For example, its diagonal elements provide the variances of m and can be used to visualize how reliable our attenuation estimates are in each location.We also consider here the resolution matrix R, which relates the true medium with the estimated one.Assuming that our forward model is perfectly accurate, we can relate the observations to a true tissue model satisfying d = Fm true and replace this in equation ( 15) to obtain The columns of R provide the relationship between reconstructed attenuation values at different locations, i.e., the tomographic point-spread functions (Korta Martiartu et al. 2020).Note that since our problem is linear, both Γ post and R are independent of m and can be precomputed to provide fast images and uncertainty estimates in practice.

Practical implementation
The proposed attenuation tomography technique is conceptually similar to the speedof-sound imaging technique CUTE.Its practical implementation thus closely follows the one typically used in CUTE.We describe its most important aspects here and refer the reader to Stähli, Kuriakose, Frenz & Jaeger (2020) for a more detailed justification of the choice of acquisition and processing parameters.

Ultrasound data acquisition and beamforming
While the theoretical derivations presented so far are valid regardless of the ultrasound probe type, this work focuses on linear transducer arrays for simplicity.We consider an acquisition sequence consisting of steered plane-wave emissions at angles φ ranging from -27.5 • to 27.5 • , with an angular step of 0.5 • .Typical linear probes have an interelement pitch close to the wavelength; thus, this angle range is limited by grating lobes becoming dominant at larger angles.The small angular step is required for the coherent compounding described in the next subsection.For each emission, we collect the backscattered signals p(x r , t) arriving at each transducer element and transform them into the analytic signals p(x r , t) using the Hilbert transform.Section 2.1 described the beamforming process in the temporal-frequency domain to simplify the theoretical derivations.In practice, however, we implement the delay-and-sum algorithm in the time domain and reconstruct the image I φ (r) corresponding to the emission along φ as using the focusing law where r = (x, z), c 0 refers to the assumed tissue speed of sound, and a is the apodization factor limiting the receiving angular aperture.Note that this factor was omitted in section 2.1 for simplicity, although it could have been included in equation ( 4).The apodization factor does not explicitly affect the derivation of the forward problem (11).However, it shapes the receive point-spread function and may therefore increase the size of Ω 0 , affecting the fulfillment of condition ( 7).Here, we use a factor a that limits the receiving angular aperture to ±30 • , which was empirically found to provide a good trade-off between the focusing quality and observational noise.Furthermore, we use a rectilinear grid for reconstructing the images I φ (r), with the lateral and axial grid spacing equal to the element pitch and c 0 ∆t/2, respectively, where ∆t refers to the sampling period of p(x r , t).

Normalized cross-correlation log-amplitude measurements
The quality of the images I φ (r) is degraded due to the clutter arising from, e.g., electronic cross-talk, grating lobes, or multiple scattering that is unaccounted for in the beamforming.To improve this, we apply coherent compounding and collimate the images synthetically on emission along the predefined angles ranging from −25 • to 25 • , with an angular step of 2.5 • (figure 1(a)).We use a Gaussian weighting function with the mean equal to the angle of interest and a standard deviation of 3 • / √ 2. This value is a compromise between reducing clutter and keeping incident beams close to plane waves, as assumed in the forward problem (11).
After coherent compounding, we compute zero-lag cross-correlations between pairs of images of successive angles (Loupas et al. 1995), and we normalize them with the corresponding auto-correlations following the equation ( 12) (figure 1(b)).We use a crosscorrelation kernel of 1 mm by 1 mm size around each grid point.While increasing its size reduces noise in log-amplitude measurements, it also reduces the spatial resolution of reconstructed attenuation images.Thus, it must be carefully chosen to optimize this trade-off.
Due to the limited aperture of the transducer array, the emitted waves do not illuminate the entire image area.We identify the nonilluminated areas by applying geometrical considerations and exclude measurements at these locations from our observations (see white areas in figure 1(b)).Finally, we interpolate the measurements into a coarser rectilinear grid with a fixed lateral and axial grid spacing of 0.5 mm to reduce redundancies, and we use the same grid for reconstructing tissue attenuation (figure 1(c)).

Calibration
The forward problem in equation ( 11) assumes negligible differences between the amplitudes of the cross-correlated plane waves; however, this does not hold in practice due to, e.g., the directivity pattern of transducer elements.To effectively remove transducer-dependent effects, we calibrate our observations with data collected in a reference phantom with known acoustic properties, as is common in state-of-the-art techniques (Yao et al. 1990, Labyed & Bigelow 2010, Pawlicki & O'Brien 2013, Coila & Lavarello 2018).With this calibration, the inverse problem in equation ( 15) becomes where d ref contains the normalized cross-correlation log-amplitudes measured in the reference phantom with attenuation m ref .

Simulation studies
This section uses numerically computed data to validate the normalized cross-correlation log-amplitude measurements and discuss the performance of the proposed technique in phantoms with varying contrast in attenuation and echogenicity.We use the k-Wave open-source toolbox (Treeby & Cox 2010) to simulate 2-D ultrasound wave propagation in lossy, nondispersive media.We consider a 256-element transducer array with an interelement pitch of 0.2 mm, emitting a 7-cycle Gaussian-modulated tone burst with a center frequency of 5 MHz.The media in our examples have a constant density of 1000 kg/m 2 , a constant speed of sound of 1540 m/s, a homogeneous background with an attenuation coefficient of α 0 = 0.5 dB/cm/MHz, and a constant power-law exponent of y = 1.For calibration data, we use a homogeneous medium with the same speed of sound and α 0 = 0.2 dB/cm/MHz.Diffuse scattering is mimicked by introducing random density perturbations with zero-mean Gaussian distribution at each point of the simulation grid, which has a grid spacing of 25 µm in both lateral and axial directions.All studies in this section use five realizations of the random scattering media and an assumed velocity of 1540 m/s for beamforming.Attenuation reconstructions are obtained with isotropic regularization unless otherwise indicated.

Validating normalized cross-correlation log-amplitude measurements
In section 2.2, we showed that the normalized cross-correlation operation is similar to placing virtual receivers inside the medium measuring attenuation-induced differences in the amplitudes of the incident waves.To demonstrate the validity of this conclusion, we consider a medium containing a circular inclusion of α 0 = 1.0 dB/cm/MHz and perform two numerical experiments.In the first experiment, we place the receivers inside the medium at every 12th simulation-grid point (0.3 mm spacing) along both directions.Moreover, we consider plane-wave emissions only at the steering angles of the coherent compounding described in section 3.2.For each recorded signal, we select the first arrivals and isolate attenuation-related amplitude information using calibration signals.Finally, we compute log-amplitude differences between pairs of emissions of successive angles, shown in figure 2(a).We refer to this data as the ground truth.In the second experiment, we consider a pulse-echo ultrasound scenario and follow the approach described in section 3 to measure the log-amplitudes of normalized cross-correlations, shown in figure 2(b).These measurements capture the ground-truth log-amplitude information remarkably well (mean absolute error (MAE): 3 × 10 −3 Np), confirming the accuracy of our theoretical derivations.Measured values at each location reflect the differences between the cumulative losses undergone by the incident waves during their propagation.Thus, they are zero at the transducer location and tend to increase with the propagation distance.Below the inclusion, we observe tails of positive and negative values.At these locations, only the propagation path of one of the incident waves crosses the inclusion, leading to large amplitude differences.Interestingly, the ground-truth data shows narrow areas of positive and negative values surrounding each tail, probably caused by the contributions of higher-order Fresnel zones occurring for waves with finite frequencies (Korta Martiartu et al. 2020).Similar areas are observed for edge waves at the margins of nonilluminated areas.These effects become less noticeable in crosscorrelation log-amplitude measurements due to coherent compounding, which acts as a collimator.This will benefit the reconstructions with our proposed technique since the forward modeling in equation ( 13) does not account for such wave phenomena.Cross-correlation log-amplitudes also show granular noise that resembles the small-scale fluctuations of tissue echogenicity.This noise is caused by the relatively low number of realizations used for ensemble averaging, which cannot completely remove the imprint of the scatterer distribution.We will discuss the role of this averaging in more detail in section 5.2 and proceed to analyze the performance of our attenuation imaging technique in the following.

Attenuation reconstruction: media with a circular inclusion
The inverse problem formulated in section 2.3 reconstructs the frequency-dependent quantity α(r) = α 0 (r)ω y(r) .Yet, throughout this work, we show reconstructed images in terms of the attenuation coefficient α 0 (r) to allow direct comparisons with the groundtruth medium properties.We remove the frequency-dependent term from α(r) by using the center frequency of the signals and the known power exponent y = 1.First, we want to understand the quality of our reconstructions relative to the highest attainable quality from this type of data.To this end, we compare images retrieved from both the ground-truth data (figure 3(a)) and normalized cross-correlation log-amplitude measurements (figure 3(b)).The corresponding lateral and axial profiles are shown in figures 3(c) and 3(d), respectively.Overall, we observe an excellent agreement between both results, with a mean absolute percentage error (MAPE) between images of 5 % and the same full width at half maximum (FWHM) in the lateral (0.90 cm) and axial (1.13 cm) profiles.However, the granular noise in crosscorrelation data (see figure 2) leads to larger errors at the bottom part of the image, slightly diminishing the contrast-to-noise (CNR) ratio.Here, the data coverage is relatively poor, and reconstructed values become highly uncertain (figure 3(f)).In both cases, the reconstructed inclusion appears elongated in the axial direction, showing that our technique provides a lower axial than lateral resolution.This is confirmed by the tomographic point-spread function shown in figure 3(e) and is caused by the limited angular aperture of the steered emissions.Such an aperture moreover leads to X-shaped artifacts crossing the inclusion, also seen in the point-spread function.Since our data cannot constrain the axial resolution well, the regularization enforces smooth features along this direction, spreading the inclusion according to the partial volume effect (Jaeger et al. 2022).At the lateral margins of reconstructed inclusions, we observe lower-value areas, which are particularly noticeable for the ground-truth data.There are two factors contributing to this: (i) data from higher-order Fresnel zones, which are neglected by our ray-based forward modeling, and (ii) the relationship between attenuation estimates at different locations shown by the point-spread function in figure 3(e).The latter is intrinsic to our tomographic technique, whereas the first factor is expected to become less relevant at higher frequencies (Korta Martiartu et al. 2020).This first factor, which is more significant for the ground-truth data (see figure 2), is moreover responsible for the line artifacts appearing in figure 3(a).
Reconstructed values within the inclusion are sensitive to the true inclusion size and the regularization.This is shown in figure 4, where we change either the diameter of the circular inclusion or the regularization strength.Compared to the result in ).This is caused by the smoothing effect of the regularization, which reduces the values of structures that are smaller than the resolution limit of our technique.The CNR, however, is the same in figures 4(a) and 4(b) since the former contains more artifacts, especially below the inclusion.Here, the large, highly attenuating area significantly decreases the signal-to-noise ratio of log-amplitude measurements.When we use a weaker regularization than in figure 3, the attenuation of the inclusion becomes more accurate (figure 4(c)), suggesting that the spatial resolution has improved.Yet, reconstruction artifacts contaminate the image, and the CNR deteriorates.A stronger regularization smooths out the noise at the bottom part of the image but underestimates the attenuation values within the inclusion, having no benefits for the CNR (figure 4(d)).The next subsection discusses different regularization strategies particularly suited for layered media.

Attenuation reconstruction: layered media
Previous results demonstrate that our technique provides images with lower axial than lateral resolution.To understand how this affects the retrieval of layered structures, we consider a homogeneous medium containing a 1-cm-thick horizontal layer with α 0 = 1.0 dB/cm/MHz.Figure 5(a) shows the reconstructed attenuation image using the same regularization as in figure 3.Although we recover the layered structure at the correct location, the regularization substantially affects its thickness and attenuation, resulting in a relatively poor CNR.Here, the axial resolution is lower than in figure 3, where it has probably been enhanced by the lateral variations in the attenuation.As previously discussed, we can improve this with a weaker regularization, shown in figure 5(b).Attenuation estimates within the layer become more accurate, but there is almost no CNR gain due to the strong artifacts appearing in areas with poor data coverage.These examples illustrate the limits of applying the same regularization strength in the axial and lateral directions.
In layered media, we expect smoother features along the lateral than axial direction, and this knowledge can be incorporated into our problem via anisotropic regularization.Figure 5(c) shows the reconstructed attenuation image with the same regularization strength as in figure 5(b) along the axial direction and 50 times stronger along the lateral direction.This approach effectively reduces most artifacts while keeping the axial resolution intact, allowing us to recover the layered structure more accurately and improve the CNR.

Attenuation reconstruction: media with echogenicity contrast
In section 2.2, we showed that the normalization step is essential to remove the contribution of tissue echogenicity from cross-correlation log-amplitude measurements.To illustrate this better, we consider the two media shown in figure 6(a).Both have the same constant attenuation, but medium II contains a circular inclusion with +6 dB echogenicity contrast.We simulate this by increasing the standard deviation of the Gaussian distribution used to introduce random density perturbations in the medium.Despite the inclusion, normalized cross-correlation log-amplitude measurements are practically identical in both cases (figure 6(b), MAE: 10 −4 Np); thus, the echogenicity contrast has no influence on the images reconstructed with our technique (figure 6(c)).While the amplitudes of cross-correlations do contain variations due to the inclusion (figure 6(d)), the normalization with the auto-correlations effectively removes this information from log-amplitude measurements, as theoretically expected from equation ( 11).These results are excellent in a context where state-of-the-art spectral-difference methods fail to provide accurate attenuation estimates in tissue with varying echogenicity (Kim & Varghese 2008, Pawlicki & O'Brien 2013).

Experimental results
This section aims to demonstrate the experimental feasibility of our attenuation imaging technique using tissue-mimicking phantoms with similar complexity as in previous examples.We also analyze here the influence of incorrect speed-of-sound assumptions on the performance of the technique and discuss the role of ensemble averaging.We acquire ultrasound data using the SuperSonic ® MACH ® 30 ultrasound system (Hologic ® -Supersonic Imagine ® , Aix en Provence, France) with the L18-5 linear probe.As in our numerical simulations, the probe contains 256 transducer elements with an interelement pitch of 0.2 mm.Acquired signals have a center frequency of 6.5 MHz with a 0.6 fractional bandwidth and a sampling frequency of 40 MHz.Data is gathered in custom-made calibrated CIRS phantoms (Computerized Imaging Reference Systems, Inc. Norfolk, VA, USA) with inclusions of varying size, echogenicity, and attenuation.Their speed-of-sound values vary less than 5 m/s with respect to the background medium, which has a speed of sound of 1518 m/s.We use this value for delay-and-sum beamforming unless otherwise indicated.We take five measurements per position for the ensemble averaging and ten measurements in a homogeneous area of the phantoms for calibration.
Figure 7 shows the results for phantoms containing cylindrical inclusions.Unlike in our numerical simulations, plane-wave emissions lead to electronic cross-talk noise in the shallow areas, which we exclude to avoid meaningless observations (see the second row in figure 7).Since our calibration data is acquired in the same phantom, observed experimental data do not contain information about amplitude losses caused by the background attenuating media.We mainly observe inclusion-related tails whose values increase in magnitude with increasing contrast in attenuation.Data noise is larger than in our numerical examples, especially at the margins of the image area, where the focusing quality of beamformed images is lower.Consequently, optimal regularization parameter values are at least one order of magnitude larger than in previous examples.As shown in figure 4(d), this leads to images with a lower axial resolution where inclusions are more axially elongated, reducing their attenuation estimates accordingly.A stronger regularization forces stronger relations between reconstructed values at different locations.As a result, we observe more pronounced X-shaped artifacts centered at the reconstructed inclusions, increasing attenuation estimates of the background medium above and below the inclusions.Similarly, the lower-value areas at the lateral edges of inclusions also become more noticeable, as shown by the reconstruction of, e.g., the largest attenuating inclusion in phantom III.Despite these regularization artifacts, which are also observed in figure 4(d), our results show reasonable average CNR values, consistent with the ground-truth inclusion properties.As in our numerical examples, the relatively good lateral resolution allows us to reconstruct structures that are 0.5 cm small (see phantom I).Moreover, our attenuation imaging technique accurately recovers the relative variations in the attenuation values of the inclusions (e.g., see phantom II and III) and remains insensitive to variations in echogenicity (see phantom IV).
Similarly, figure 8 shows the experimental results in phantoms containing a horizontal-layer inclusion with contrast in either attenuation or echogenicity.Following our observations in section 4.3, we use anisotropic regularization satisfying λ x = 50λ z .We keep this λ x /λ z ratio constant in order to optimize the regularization parameter values using the L-curve method (see supplementary figure 2).
Although the experimental noise leads to more smoothing than in our numerical examples, and this particularly affects the already limited axial resolution of our technique, we still recover the attenuating inclusion with relatively good CNR and in the correct location (see phantom V).However, the inclusion becomes more diffused in the axial direction, increasing the values of the reconstructed background medium.For phantom VI, our attenuation reconstruction is approximately homogeneous with a variability of 0.05 dB/cm/MHz.There are a few artifacts in the top right part of the image (low and high attenuating areas), but the overall influence of the echogenicity contrast remains relatively low.This example thus confirms the robustness of our technique against variations in echogenicity, as already observed in the previous section.

Sensitivity to beamforming speed of sound
Previous results used the true speed-of-sound value for delay-and-sum beamforming.This ensures optimally focused point-spread functions, allowing us to fulfill the condition (7) required for the accuracy of our forward problem.In practice, however, the true tissue speed of sound usually deviates from the assumed speed of sound for beamforming, and aberrations will degrade the focusing quality of beamformed images.To understand how this affects our results, figure 9 compares attenuation images reconstructed with the beamforming speed of sound c 0 ranging from 1480 m/s to 1560 m/s.As an example, we consider the phantom III in figure 7. The results confirm that cross-correlation log-amplitudes measurements become noisier the larger the mismatch between the true and assumed speed of sound is (see the top row in figure 9).This noise increases with depth as aberrations become stronger and deteriorates the CNR of reconstructed attenuation images.Interestingly, we observe an increase in attenuation estimates with increasing c 0 .This could indicate that the mislocation of echoes (i.e., our virtual receivers) magnifies the existing gradient with depth in log-amplitude measurements, thereby introducing an apparent shift in attenuation estimates.Importantly, the reconstructed contrast of inclusions appears insensitive to the choice of c 0 , meaning that our technique reliably captures the relative attenuation variations of the medium.

The role of ensemble averaging
Our forward problem relates the attenuation in tissue to ensemble-averaged crosscorrelations.We perform this averaging by taking five repeated measurements at the same location, slightly moving the ultrasound probe.Since this approach can be prone to motion artifacts in vivo, it is important to understand how the ensemble averaging affects the quality of our reconstructions.In figure 10, we show the observed data and reconstructed attenuation images for phantoms I-IV using only one experimental acquisition.Compared to figure 7, we notice a considerable increase in data noise, particularly at the margins where the focusing quality is the lowest.This noise is random and granular, meaning that the averaging allows us to reduce the imprint of random medium scatterers from the log-amplitude observations.This data noise increases the overall variability of reconstructed attenuation values, which translates to generally lower CNR values, although we still recover the circular inclusions relatively well.Therefore, averaging over repeated acquisitions is not critical to obtain meaningful attenuation images, although it improves their contrast resolution.

Discussion and conclusions
This work presents the first fully 2-D ultrasound attenuation imaging technique for pulseecho systems.Unlike state-of-the-art techniques, which estimate tissue attenuation from the axial variations of recorded echoes (Coila & Lavarello 2018, Labyed & Bigelow 2011), our approach uses observations of echo-amplitude changes between emissions at different steering angles.We thus constrain the attenuation values at each tissue location using multiple crossing wave paths, essential to resolve the spatial variations of this tissue property (Rawlinson & Spakman 2016).Furthermore, the proposed approach appears robust against echogenicity contrasts, which are an important source of artifacts for state-of-the-art amplitude-based techniques.These methods neglect the amplitude contributions arising from variations in tissue echogenicity, leading to inaccurate attenuation estimates at the interfaces (Kim & Varghese 2008).In contrast, our approach leverages beamformed images obtained for different emissions to effectively remove their imprint from amplitude observations.The images reconstructed with our technique show significantly better lateral than axial resolution, which is hampered by the limited angular aperture offered by pulseecho systems.Current ultrasound probes have an inter-element pitch close to the wavelength of emitted pulses, and grating lobes limit the maximum usable steering angle.The acquisition sequence chosen in this study has been optimized to minimize their influence while maximizing the angular aperture to improve the axial resolution.The small angular step allows us, through the coherent compounding process, to reduce the intensity of echoes arising from grating lobes relative to echoes of the main lobe.The optimal acquisition and processing parameters, however, depend on the respective ultrasound probe in use.High-quality probes featuring a larger element count at a smaller inter-element pitch could allow us to increase the usable angular aperture and thus improve the axial resolution.Regardless of the probe, we could explore other approaches to improve the axial resolution, for instance, by incorporating into our technique observations of axial echo-amplitude changes used in state-of-the-art approaches.
We base the theoretical derivations of the forward model on two major assumptions.First, we consider tissue as a diffuse scattering medium consisting of randomly and uniformly distributed scatterers.In reality, however, tissue also contains specular reflectors.They reflect ultrasound waves predominantly in directions that depend on the incident steering angle, and the amplitude of echoes detected at these reflectors will vary between different emissions, independently of tissue attenuation.This can introduce outliers in log-amplitude measurements, as we observe at the top and bottom edges of the inclusions in figure 7.For our experimental examples, only very few measurements showed such values, so they did not significantly affect our results.However, large measurement errors could become important sources of biases for our least-squaresbased reconstruction algorithm.It may be interesting to study the performance of our technique with optimization strategies that are more robust to outliers, for instance, by minimizing the ℓ 1 -norm (Sanabria et al. 2018).
Second, we have neglected amplitude terms related to the (de)focusing, diffraction, and interference effects caused by tissue speed-of-sound heterogeneities (Dahlen & Baig 2002).If present, they will be mapped as attenuation variations in the reconstructions.Such inaccuracies are not exclusive to our technique but will affect any attenuation imaging method based on amplitude observations.The reconstruction of phantom IV in figure 7 serves as an example to illustrate speed-of-sound-related artifacts.Here, the inclusion on the left has a negligible contrast in attenuation, but cross-correlation phases reveal a significant speed-of-sound contrast with respect to the background (see supplementary figure 3).This is probably why we observe tails in log-amplitude measurements, leading to artifacts at the edges of the inclusion.While these are not concerning in this example, future work should focus on understanding the effects of larger speed-of-sound contrasts and design strategies to disentangle their amplitude contributions from those related to attenuation.
Our inverse problem is defined to retrieve the frequency-dependent attenuation of tissues.This is a combination of two parameters characterizing tissues: the attenuation coefficient α 0 and the frequency power law exponent y.In our experimental results, we have assumed a constant value y = 1 to compare the reconstructed α 0 values with those provided by the phantom manufacturer, who reports the averaged measured values in the frequency range 2.25 -5.5 MHz.The calibrated exponent values of the inclusions, however, vary between 0.86 to 1.17 according to the manufacturer.This may explain part of the mismatch between certified and estimated values observed in figure 7. To characterize tissue attenuation accurately, both α 0 and y should be reconstructed in a multi-parametric fashion.For this, we could adapt our technique, for instance, following the approach suggested by Rau et al. (2021) and Chintada et al. (2022), where we first reconstruct an image per frequency band and then fit a power-law function at every spatial location.In this case, a reduced frequency bandwidth will enlarge the pointspread functions and increase the data noise level, as discussed in section 2.2.The feasibility of this approach should therefore be carefully studied in future works.
Similar to other attenuation imaging techniques, calibration is critical to obtain meaningful results with our method.It allows us to minimize amplitude contributions arising from transducer-and system-dependent effects, such as the directivity of transducer elements or time-gain compensation.This approach makes the quantitative accuracy of our images very dependent on a well-calibrated reference phantom and its acoustic properties, which should be similar to the properties of the tissue under study (e.g., see supplementary figure 4).It would be convenient to reduce this dependency by exploring calibration strategies that use, for example, adjacent frequencies to cancel the unwanted effects (Gong et al. 2019).
The attenuation imaging technique presented in this work relies on similar principles as the speed-of-sound tomography approach CUTE, with comparable computational and technical requirements.One of the main differences between both modalities is the type of observations extracted from the backscattered ultrasound signals: the phases of cross-correlations between beamformed images are related to the speed of sound, whereas their amplitudes are sensitive to attenuation in tissue.Therefore, both imaging modalities are complementary and can be easily combined into a single framework to provide simultaneous images of tissue attenuation and speed of sound.This work thus represents an important step towards exploiting the full waveform content of ultrasound recordings to characterize tissue acoustic properties comprehensively.

Figure 1 :
Figure 1: Workflow for pulse-echo ultrasound attenuation tomography.(a) We insonify the tissue by sequentially emitting a set of steered plane waves and reconstruct a complexvalued image per predefined emission angle using delay-and-sum beamforming and coherent compounding.For illustrative purposes, we show B-mode images and increase the angular step to 10 • .(b) We extract echo log-amplitude changes between emissions by cross-correlating every successive pair of images following equation (12).Nonilluminated areas, shown in white at the margins of the maps, are excluded from the observations.(c) We solve the regularized linear inverse problem in (15) to retrieve the spatial distribution of attenuation in tissue.The example shown here uses numerically computed ultrasound signals in a medium that contains a circular inclusion with +0.5 dB/cm/MHz attenuation contrast with respect to the background.

Figure 2 :
Figure 2: Log-amplitude differences between emissions steered at the angles indicated at the top.The medium is homogeneous with a circular inclusion (dashed lines) of positive attenuation contrast.(a) Ground truth.(b) Normalized cross-correlation log-amplitude measurements.Only three exemplary maps out of a total of 20 maps are shown in each case.

Figure 3 :
Figure 3: Reconstructed attenuation image using (a) the ground-truth (GT) data in figure 2(a) and (b) the cross-correlation (CC) data in figure 2(b).Both results use the same regularization parameter value λ = 2 × 10 −6 , which is optimized for the second reconstruction.The contrast-to-noise ratio (CNR) is computed by taking the values within and outside the true location of the inclusion, covering the whole image area.(c)-(d) Lateral and axial profiles of reconstructed images through the center of the inclusion.Dashed lines show the true attenuation coefficient of the medium.(e) Tomographic point-spread function (PSF) located at the center of the inclusion, together with its lateral and axial profiles.(f) Normalized variance of reconstructed attenuation values, illustrating the relative uncertainties of attenuation estimates at each spatial location.

Figure 6 :
Figure 6: (a) B-mode images, (b) an example of normalized cross-correlation log-amplitude measurements (steering angles: 12.5 • -15 • ), and (c) reconstructed attenuation images for a homogeneous medium (I) and the same medium containing a circular inclusion with +6 dB echogenicity contrast (II).Reconstructions use λ = 2 × 10 −6 .(d) Log-amplitudes of crosscorrelations and auto-correlations used to compute the measurements in (b) for phantom II.The values are normalized with respect to the maximum log-amplitude in each case.

Figure 7 :
Figure 7: Experimental results in tissue-mimicking phantoms containing cylindrical inclusions with contrast in attenuation and echogenicity.The B-mode images in the first row show the shape and location of the inclusions and their echogenicity contrast with respect to the background medium.In the second row, we show an example of normalized cross-correlation log-amplitude measurements (steering angles: 12.5 • -15 • ).Dashed lines indicate the location of the inclusions.The last two rows display the reconstructed attenuation images, where we indicate the average contrast-to-noise ratio (CNR) of both inclusions, and their lateral profiles through the center of the inclusions.The profiles also show the calibrated attenuation values provided by the phantom manufacturer (dashed gray line).All reconstructions use λ = 3•10 −5 or λ = 8 • 10 −5 depending on the L-curve optimization results (see supplementary figure 1).

Figure 8 :
Figure 8: Experimental results in tissue-mimicking phantoms containing a horizontal inclusion with contrast in either attenuation (phantom V) or echogenicity (phantom VI).B-mode images in the left column indicate the value of the echogenicity contrast in each case.Reconstructed attenuation coefficient images and their axial profiles are shown in the middle and right columns, respectively.The profiles compare the calibrated attenuation coefficient values provided by the phantom manufacturer (dashed gray line) with the estimated values at each lateral position (solid black line) and position x = 0 (solid red line).Optimal regularization parameter values are λ z = 7 • 10 −6 and λ z = 1 • 10 −5 , respectively.CNR: Contrast-to-noise ratio.

Figure 9 :
Figure 9: Sensitivity of experimental results to the speed-of-sound value c 0 used for delayand-sum beamforming.The first row shows the normalized cross-correlation log-amplitude measurements for steering angles 12.5 • -15 • .The middle and bottom rows display the reconstructed attenuation images and their lateral profiles, respectively.Considered c 0 values in each case are indicated at the top.All reconstructions use λ = 3 • 10 −5 .CNR: Contrast-tonoise ratio.

Figure 10 :
Figure 10: Normalized cross-correlation log-amplitude measurements (steering angles: 12.5 • -15 • ) and reconstructed attenuation images for phantoms in figure 7 with only one experimental realization.Dashed lines indicate the location of the inclusions.Reconstructions use the same regularization as in figure 7. CNR: Contrast-to-noise ratio.