As a broadband ultrasound wave passes through a medium, its waveform changes according to the phase velocity and attenuation of the medium. Ultrasound waves travelling through a viscoelastic medium can be defined in the frequency domain as: U d , ω = G ω e − i k ̂ ω d (1a) w i t h k ̂ ω = ω c ω − i α ω (1b)

where U (d, ω) is the FT of the ultrasonic wave u (d, t), ω is the angular frequency, d is the propagation distance, k ̂ ( ω ) is the complex wave number, c(ω) is the phase velocity, and α(ω) is attenuation at ω, and G(ω) accounts for the transmit and receive transducer responses at frequency ω.

As seen in Eq. 1a, ultrasonic wave amplitude decreases as a function of propagation distance d due to absorption, relaxation, scattering, etc. Indeed using this wave equation, ultrasound phase velocity and attenuation can be derived using amplitude and phase spectrums by taking the natural logarithm of amplitude |U (d, ω)| and separating the phase terms to equate to phase angle ∠U (d, ω), which yield respectively: ∠ U d , ω = ω c ω d + ∠ G ω , (2a) ln | U d , ω | = ln G ω − α ω d . (2b)

Phase and amplitude spectra for the reflector echo corresponding to Tx element t and Rx element r relates, respectively, to slowness s _p (ω) = 1/c(ω), i.e., the inverse of phase velocity, and attenuation α(ω) along the traversed ray path p as follows: ∠ U t , r ω = ω ∫ p s ω d l + ∠ G t , r ω , θ t , r , (3a) | U t , r ω | = exp − ∫ p α ω d l + | G t , r ω , θ t , r | . (3b)

where p represents the acoustic ray path and l is the traversed distance. ∠G _t,r (ω, θ _t,r ) and |G _t,r (ω, θ _t,r )| are the confounding factors attributable to transducer Tx-Rx phase and amplitude transfer characteristics at different angles (e.g., the aperture opening and side-lobes) as well as the acoustic reflector reflection characteristics, i.e., the incidence angle dependent specular and diffuse reflection and transmission characteristics of the reflector surface. These have to be compensated for in order to be able to measure the actual ultrasound phase velocity and attenuation effects. For this purpose, we calibrate these measurements by normalizing with measurements in water ( U t , r w a t e r ( ω ) ) at a known water temperature and using the same ultrasound pulses (bandwidth) as in the targeted tissue imaging setup. From the given water temperature, we compute the ground-truth water SoS using the relationship in [42], which was 1,483.1 m/s in our setup. Water UA is known from [43] to be α(f) = 2.17 × 10⁻¹⁵ f ².

We used the procedure described in [34] to compensate for the confounding effects in |G _t,r (ω, θ _t,r )|. The reflector does not change the phase of the reflected ultrasound signals, since its acoustic impedance is higher than water and soft tissues. Thus, for factoring out the confounding effects in ∠G _t,r (ω), we used simply the phase angles of reflector echo profiles from water calibration experiments. We denote these calibrated amplitude and phase spectrum values as | U t , r ′ ( ω ) | and ∠ U t , r ′ ( ω ) , respectively, which are given by ∠ U t , r ′ ω = ∠ U t , r ω − ∠ U t , r water ω (4a) | U t , r ′ ω | = U t , r ω R water θ U t , r water ω R θ . (4b)

where R _water(θ) and R(θ) are the ultrasound reflection coefficients of the reflector in water and in targeted imaging medium, respectively, for incident angle θ. It can be estimated using Snell’s law with the SoS in water (c _water) and an approximate SoS of the targeted medium, as in [34].

For acoustic propagation, refractions are herein ignored and the travel of acoustic wave from a transmitting to receiving element is modeled using a straight ray approximation as depicted in Figure 1A. For a straight acoustic path, the first wavefront from a Tx element, reflecting from the reflector and arriving at an Rx element can be assumed to follow the shortest path, thus being reflected at the mid-point between these Tx and Rx element locations projected on the reflector, cf. Figure 1A. To cast this as a tomographic reconstruction problem, the relation in Eq. 3 above can be discretized on a Cartesian grid as ∠ U t , r ′ ω = ω ∑ n = 1 N i s k ω l i k , (5a) | U t , r ′ ω | = exp − ∑ n = 1 N i α k ω d l . (5b)

where s _k (ω) and α _k (ω) represents the slowness and attenuation, respectively, in pixel k along ray i spanning N _i pixels, where l _ik represents the partial acoustic path length of ray i within this pixel k.

Collocating Eq. 5 for each ray path along all Tx-Rx combinations and taking logarithm on both sides of Eq. 5b, the forward-problem of ultrasound propagation can then be represented with the following systems of linear equations ω L s ω = ϕ ω , (6a) L α ω = b ω , (6b)

where ϕ(ω) and b(ω) are, respectively, the column vectors of phase ∠ U t , r ′ ( ω ) and amplitude ln | U t , r ′ ( r , ω ) | values along all ray paths at frequency ω; s(ω) and α(ω) are the column vectors of, respectively, the slowness and attenuation values for all reconstruction grid pixels at frequency ω; and the system matrix L encodes the discretized ray integrals for each ray onto the pixel grid as in [25]. Note that L is constant for given transducer geometry and reflector position.

Given calibrated phase spectrum ϕ(ω) and amplitude spectrum b(ω) measurements, one can then reconstruct the local slowness and UA maps at frequency ω by solving the following inverse problems: s ̂ ω = arg min s ω ‖ ω L s ω − ϕ ω ‖ 1 + λ ‖ D s ω ‖ 1 (7a) α ̂ ω = arg min α ω ‖ L α ω − b ω ‖ 1 + λ ‖ D α ω ‖ 1 (7b)

where λ is the weight parameter of regularization required to help robustly solve these ill-posed problems. Similarly to [22,25,32,34], we herein use l ₁-norm for both the data and regularization terms for robustness to outliers in, respectively, the measurements and the reconstructions [44]. For regularization, in order to suppress the typical streaking artifacts associated with LA-CT due to missing orthogonal projections as well as the uneven distribution of measurement paths, we use anisotropic weighting of directional gradients proposed in [45], with a weighting of k = 0.9 as in [25]. Accordingly, the regularization matrix D is composed of a mix of Sobel (axis-aligned) and Robert (diagonal) kernels with corresponding weights and locations, as in [23]. In this paper we empirically set λ = 0.6 for all experiments of phase velocity and UA reconstructions. For the numerical solution of the optimization problem Eq. 7, a limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm [46–49] is used from the unconstrained optimization package minFunc ¹ .

Ultrasound phase velocity c(ω) and attenuation α(ω) reconstructions are computed for a range of frequencies [ω _c − nΔω, … , ω _c , … , ω _c + nΔω] around the center frequency ω _c of the ultrasound pulse with Δω being the frequency resolution used for computing the FT of the reflector profiles. Sample reconstructions at three frequency values are shown in Figures 3A,B. Phase velocities for wide range of frequencies may exhibit complex profiles as in [50–52]. Nevertheless, for small ranges of ultrasound imaging frequencies, such frequency dependence can be approximated with a linear disperse model, similarly to the shear-wave phase velocity studies in [16,53]. Given the relatively small bandwidth of the utilized ultrasound pulse, the phase velocity is expected to vary linearly within this bandwidth, as illustrated by two sample profiles in Figure 3C. For frequency characterization of attenuation, it is well known in the literature that attenuation as a function of frequency obeys the power-law, which corroborates the observation in Figure 3D. One could fit these models pixel-wise for the given range frequencies and solve for the model parameters in the least square sense, as in [34]. However, this may be error prone as it ignores the spatio-spectral continuity in the medium. Herein, we incorporate frequency-dependent models of speed and attenuation into the reconstruction process, such that the inverse problem solution is informed by and can leverage the measurements at each frequency and path concurrently, while estimating the frequency-dependent medium parameters.

As models of the frequency-dependent nature of ultrasound phase velocity and attenuation, we use herein the following linear and power-law relations, respectively: c f = c 0 + c y f (8a) α f = α 0 f y ; (8b)

where c(f) and α(f) are the frequency-dependent SoS and attenuation, c ₀ and α ₀ the medium-specific SoS and attenuation coefficients, c _y the medium-specific SoS dispersion coefficient, and y the medium-specific attenuation exponent.

By colocating SoS reconstructions for all measurement frequencies from Eq. 7a in Eq. 8a, a comprehensive system of equations that encodes both SoS model parameters for all image locations is arrived as follows: θ ̂ = arg min θ ‖ A θ − v ‖ 1 + λ ‖ D ′ θ ‖ 1 (9) where θ = c 0 c y T , v = c f 1 c f 2 ⋯ c f n T , D′ = [D D], and A = I n p f 1 ∗ I n p I n p f 2 ∗ I n p ⋮ ⋮ I n p f n ∗ I n p with bold variables denoting size-n _p row vectors of the sought parameters in the imaged field-of-view, I n p is the identity matrix of size n _p and n _p is the total number of reconstruction grid pixels.

For estimating UA parameters, we colocate UA reconstructions for all measurement frequencies from Eq. 7b in Eq. 8b. Taking the logarithm of both sides (in order to drop down the exponent) a similar linear problem as in Eq. 9 is arrived, but this time with the variables θ = log α 0 , y T , v = log α f 1 , log α f 2 , … , log α f n T , and A = I n p log f 1 ∗ I n p I n p log f 2 ∗ I n p ⋮ ⋮ I n p log f n ∗ I n p .

Solving these inverse problem formulations we estimate the spatial maps of SoS model parameters c ₀ and c _y , and of UA model parameters α ₀ and y. Sample reconstructions from the earlier simulated example can be seen in Figures 3E–H.

Given ultrasound phase velocity c(ω) and attenuation α(ω) maps reconstructed above, complex bulk modulus K ̂ can be determined in a model-independent way to characterize the viscoelastic nature of the tissue through its relationship to complex wave number as follows: K ̂ ω = K ′ ω + i K ″ ω = ρ ω 2 k ̂ 2 ω , (10)

where ρ represents tissue density. Substituting Eq. 1b in Eq. 10 allows to derive the bulk storage modulus K′ and bulk loss modulus K ^″ as follows: K ′ ω = ρ ω 2 ω c ω 2 − α ω 2 ω c ω 2 + α ω 2 2 , (11a) K ″ ω = 2 ρ ω 2 ω c ω α ω ω c ω 2 + α ω 2 2 . (11b)

In this study, the following metrics are used for a quantitative analysis of the simulation results:

• Root-mean-squared-error (RMSE)

      R M S E = 1 n p x ⋆ − x ̂ 2 , where x represents a reconstructed parameter map with the number of pixels n _p . The ⋅^⋆ and ⋅ ̂ indicate, respectively, the ground-truth and a reconstruction.

• Contrast-ratio fraction (CRF)

      C R F = C ̂ C ⋆ , where C = | μ inc − μ bg | | μ inc | + | μ bg | with mean inclusion and background values μ _inc and μ _bg, respectively, where the inclusion is delineated using the ground-truth map.

• Contrast-to-noise ratio (CNR)

      C N R = | μ inc | − | μ bg | σ inc 2 + σ bg 2 , where σ ² represents the variance.

Numerical simulations were conducted to study the reconstruction accuracy of the proposed methods. These simulations were conducted by varying local distribution of SoS and UA patterns and varying their frequency dependency characteristics. An open-source acoustics toolbox k-Wave [54] was used for this purpose. k-Wave uses Kramers–Kronig relations [50–52] to simulate the dispersive phase speed based on given spatial maps of UA coefficient α ₀ (x, y) and SoS c f c ( x , y ) at the center frequency, as well as the power-law dependency parameter y. Note that in k-Wave, y is a single constant for the entire domain/simulation and cannot be controlled spatially.

The transducer was simulated as a linear array of 128 elements with a pitch of 0.3 mm. Full-matrix multistatic data was simulated at a center frequency of 5 MHz with a 5-cycle input pulse. For the simulation, we used a temporal resolution of 160 MHz and a spatial grid resolution of 37.5 μm. In the physical experimental setup, the acoustic mirror is made of plexiglass, so it was simulated using plexiglass’s acoustic properties (density ρ = 1180 kg/m³, SoS c = 2700 m/s), as positioned at a distance of 42 mm from the transducer surface in all three simulations. A single circular inclusion with 5 mm radius is placed at the center of the imaging field-of-view. SoS at center frequency ( c f c ) in the background and inclusion were set, respectively, 1,540 m/s and 1,601.6 m/s (for 4% contrast). The inclusion and background attenuation coefficient (α ₀) were set, respectively, 0.2 dB/cm⋅MHz ^y and 0.1 dB/cm⋅MHz ^y .

We repeated the simulations three times by varying the power exponent y = {1.1, 1.5, 1.9}, which led to three numerical phantoms referred hereafter as sim{#1, #2, #3}, respectively. Note that given constraints in k-Wave, SoS and UA cannot be controlled independently and arbitrarily, and the variation of such exponent y changes both SoS and UA with a known relation.

These were conducted using the data from [34]. In the first phantom #A we used, the background has 10% gelatin 1% Sigmacell Cellulose Type 50 (Sigma Aldrich, St. Louis, MO, United States), into which a bovine skeletal muscle sample was inserted as an inclusion. For the second phantom #B, the background and inclusion compositions were interchanged, i.e., a gelatin phantom piece was inserted as the inclusion inside a bulk muscle sample. Data acquisitions were conducted on a research ultrasound machine (Verasonics, Kirkland, WA, United States) with a 128-element linear-array transducer (Philips, ATL L7-4). Similarly to the simulation setup, a wideband Tx pulse with a center frequency of 5.2 MHz and a pulse length of five cycles was used. For the calibration procedure, multi-static data sets were acquired in distilled water at room temperature by placing the reflector at multiple depths {31, 35, 39, 43, 47} mm. For any intermediate reflector depths used in the phantom measurements, calibration data was interpolated from the available water measurements above, as described in [32].

SoS and UA reconstructions for sim#2 are presented in Figure 2. Reconstructed c ₀ (i.e., SoS at the center frequency c (f _c )) and the linear dispersion coefficient c _y for the three simulations are depicted in Figure 4. From the ground-truth and reconstructed SoS maps at the center frequency in Figures 4A,B, it can be seen that the background SoS at the center frequency is reconstructed with some success in all cases. However, the inclusion SoS values are seen to be shifted towards background values, which corroborates the findings in [25]. This is due to the limited-angle tomographic nature of the problem, where the lack of lateral projections combined with the regularization needed for robust solutions cause the smoothing and axial elongation of the inclusions, hence spreading them over a larger area. Since the cumulative SoS effect shall stay the same, per-pixel inclusion values are effectively averaged with the background, inversely proportional to such artifactual area increase.

Note that as y varies across the simulations, so does c _y , according to Kramers–Kronig relationship [51]. Therefore, the c _y maps in Figures 4D are expected to present contrast between the three simulations with different y values. An analytical computation of expected SoS dispersion in the utilized frequency range of 4–6 MHz according to the Kramers–Kronig relationship is illustrated in Figure 5. Accordingly, in this frequency range a linear SoS dispersion (c _y ) of {6.5, 9.7, 5.3} cm/s⋅MHz for the background and {14.1, 21.0, 11.4} cm/s⋅MHz for the inclusion are expected for each respective simulation. These are used as ground-truth c _y values for the following evaluation.

For a quantitative evaluation, we computed the metrics RMSE, CNR, and CRF for each simulation, as reported in Table 1. We computed the mean and standard deviation for the background and inclusion using their respective masks from simulations. The average RMSE of c (f _c ) is 8.50 ± 0.24 m/s and the average mean values of c (f _c ) in the background region is 1,542.3 ± 0.1 m/s, demonstrating that an accurate reconstruction of the speed of sound is possible. The contrast metrics CNR and CRF for c (f _c ) and for c ₀ indicate that the inclusions can be successfully distinguished in SoS reconstructions, with a contrast similar to their prescribed ground-truth value. Note that c _y reconstructions for all the simulations show large RMSE errors of 39–53 cm/s⋅MHz. This is mainly due to the relatively minute dispersion values of [5.3–9.7] cm/s⋅MHz and hence their minimal effect in SoS change within the given frequency interval (see the changes in y-axis in Figure 5 being < 0.4 m/s for inclusion < 0.2 m/s for the background). Note that this is partly due to our earlier-mentioned experimental limitations with k-Wave parameterizations, where SoS frequency dispersion can only be generated by varying y, which then has a much more prominent impact on UA. To capture such variation, very high SNR and precision would be required in SoS reconstructions. Despite the accuracy being low for these, nevertheless, the relative differences and contrasts are recovered successfully: For instance, c _y for the inclusion of sim #2 is higher than that of sim #1 and #3, which corroborates the analytical expectation from Figure 5. Similarly, corroborating the analytical expectation, the inclusions in sim #1 and #2 show roughly twice the dispersion that their backgrounds. The average ratio of μ i n c c y to μ b g c y is ≈2.4, which is close to the ground truth c _y contrast of 2 between inclusion and background. Therefore, despite not being accurate reconstructed c _y values are still correlated with ground truth counter-parts and hence may potentially act as an imaging biomarker.

Reconstructed UA at the center frequency α _fc , UA coefficient α ₀, and exponent y are plotted in Figure 6. Given the groundtruth and reconstructed α _fc maps in Figures 6A,B, it can be seen that the UA at center frequency is successfully reconstructed in all cases. In Figures 6B,C it is observed that the reconstructed α _fc and α ₀ values in the inclusions are shifted towards the background values, similarly to SoS results in Figure 4B and also in line with previous observations in [34]. The UA exponent y maps in Figure 6C should not present any contrast since the power exponent y are constant across each domain. Although this is true for sim #1, the other two simulations exhibit a slight deviation of y in the inclusion region. This is due to the relatively more pronounced artifactual elongations of the inclusion in UA reconstructions with increased frequency dispersion, which then slightly reduce the per-pixel reconstructed UA value (through averaging as reasoned above), which in turn is erroneously attributed by the model fitting to an increased y in the inclusion.

To quantify UA reconstructions, we report RMSE, CNR, and CRF in Table 1. Average RMSE of α(f _c ) is 0.16 ± 0.09 dB/cm at center frequency 5 MHz. Considering the range of ground-truth SoS values, this demonstrates that an accurate reconstruction of UA is possible. With an RMSE of {0.04, 0.05, 0.06} respectively for prescribed y values {1.1, 1.5, 1.9}, the relative estimation error becomes below 3.3%, indicating a robust estimation of frequency exponent y. The contrast metrics CNR and CRF for α(f _c ), α ₀, and y demonstrate that the inclusions can be successfully distinguished in UA reconstructions, with a contrast similar to their prescribed ground-truth values.

The SoS reconstruction results of the gelatin phantom and ex-vivo experiments are shown in Figure 7 with their corresponding B-Mode images in Figure 7A. Reconstructions of SoS maps at the center frequency are shown in Figure 7B, while reconstructed c ₀ and c _y are shown in Figures 7C,D. The reconstructed results in Figure 7B were overlaid with inclusion markings annotated from the B-mode images. A clear contrast is visible in Figure 7 for all our proposed reconstructions. The reconstructed maps c (f _c ) and c _y for gelatin and bovine tissue regions are consistent in phantoms #A and #B; i.e., the reconstructions show inverted images between these two phantoms, indicating that these quantities can be reproducibly used for tissue differentiation. From Figure 7C, it is observed that the tissue samples show much higher dispersion than in the simulations. This indicates the selection of frequency for SoS quantification being an important factor in practice. This also signifies the need for imaging C ₀ and C _y to image the SoS of the targeted medium comprehensively. In other words, UA dispersion estimation can help both as additional parameters to characterize tissues and also for disambiguation of SoS measurement, which are seen here to be highly frequency dependent; also corroborating tissue observations in the literature, e.g. [33,39].

The UA results of the gelatin phantom and ex-vivo experiments are shown in Figure 8. A clear contrast is visible for all the UA reconstructions using our proposed methods. Furthermore, all three reconstructed parameters between phantoms #A and #B are consistent, i.e., inclusion and background values inverted, showing reproducibility under different experimental settings. Axial edge artifacts, especially in y and α ₀ maps were much reduced compared to results in [34] thanks to solving model parameters collectively in a closed form as in Eq. 9. The average and standard deviation values of SoS and UA obtained in gelatin and bovine tissue were tabulated along with CNR of the reconstructions in Table 2.

Using the relationship in Eq. 11, the complex bulk modulus of the gelatin phantom and ex-vivo tissues at ultrasound center frequency can be computed to be 2.31 + i (2.42 × 10^–3) GPa and 2.33 + i (10.01 × 10^–3) GPa, respectively. These values indicate 1) that the bulk loss modulus in both media is three orders of magnitude smaller than the bulk storage modulus; and 2) that although both media have similar bulk storage modulus, the bulk loss modulus in ex-vivo bovine tissue is ≈4 times higher than that of gelatin. This striking difference indicates that the complex bulk modulus may act as a new imaging biomarker, where the loss modulus component potentially having superior tissue differentiation compared to its storage component.