Journal of Medical Ultrasonics (2019) 46:27–33
https://doi.org/10.1007/s10396-018-0915-9

ORIGINAL ARTICLE

Estimation method for sound velocity distribution for high‑resolution
ultrasonic tomographic imaging
Keiichiro Abe1 · Mototaka Arakawa1,2 · Hiroshi Kanai1,2
Received: 20 March 2018 / Accepted: 3 October 2018 / Published online: 15 November 2018
© The Japan Society of Ultrasonics in Medicine 2018

Abstract
Purpose With commercial ultrasonic equipment, the sound velocity is fixed to a constant value of 1530 or 1540 m/s, which
is used for beam formation. However, the assumption of a constant sound velocity is not optimal, as the sound velocity in aliving body is heterogeneous. In this study, a novel method was proposed to estimate the distribution of the sound velocity
in a region of interest.
Methods The sound velocity distribution was estimated by fitting the theoretical propagation time of the ultrasonic wave
from the scatterer to each of the probe elements with measured values.
Results In a phantom experiment, the sound velocity distribution was estimated by the proposed method with a maximum
estimation error of 0.6%, and the resultant local sound velocity values successfully improved the quality of the ultrasonicimage.
Conclusion The proposed method has the potential to improve ultrasonic image quality in in vivo experiments by estimating
the sound velocity distribution.
Keywords Ultrasound imaging · Ultrasound velocity

Introduction
Ultrasound imaging can be employed to characterize soft
tissues. It can be repeatedly applied owing to the high levelof safety and low cost [1–5]. Therefore, it is very useful inthe diagnosis of various diseases, including cardiovasculardiseases [6].
Linear and sector array probes are often used in ultrasound imaging. In conventional in vivo ultrasound imaging
methods, the sound velocity in a tissue is typically assumed
to be constant in the delay-and-sum beam-forming (in manycases, it is assumed to be 1540 m/s) [7]. However, when thesound velocity of a tissue is inhomogeneous in a clinicalcondition, the spatial resolution of a conventional beamformer deteriorates [8–10]. This issue is particularly conspicuous in the diagnosis of obese patients [11, 12], as the* Keiichiro Abe
abekei@ecei.tohoku.ac.jp
1

Graduate School of Engineering, Tohoku University,
Sendai 980‑8579, Japan

2

Graduate School of Biomedical Engineering, Tohoku
University, Sendai 980‑8579, Japanaverage sound velocity of fat (1450 m/s) is different fromthat of other soft tissues (1540 m/s). This problem is alsopresent in the diagnosis of breast diseases [13, 14].
Determination of the sound velocity distribution in a tissue not only improves the image quality but can also be
used for discrimination between normal and abnormal tissues for a diagnosis of various diseases. For example, thesound velocity of the normal liver decreases when the liver
contains fat [15, 16]. Therefore, it is possible to diagnose
fatty liver by estimating the sound velocity distribution.
Many studies have been performed to correct the timedelay caused by tissue inhomogeneity. The cross-correlation
technique is one of the most common methods for phaseaberration correction [17, 18]. This method is highly accurate in the correction of the time delay caused by inhomogeneity [19]; however, it is susceptible to noise and its computational load is high [12, 20]. In addition, studies havebeen performed to estimate the optimal sound velocity forbeam-forming to improve image quality. Methods have beendeveloped to estimate the optimal sound velocity for imaging from the focusing quality [12, 21, 22]. The crossed beam
method [23] has also been proposed as a method for estimating the local sound velocity in the living body. However, it is13Vol.:(0123456789)

28

necessary to know the sound velocity between the ultrasonicprobe and the region of interest (ROI) when estimating the
local sound velocity. There is a large variation in the estimated sound velocity.
In this study, a new method was developed to estimate
the sound velocity distribution from the time delay values
in the radio-frequency (RF) signals received at 96 elementsof an ultrasonic probe. This method utilizes the curvature
set of the time delay, which depends on the average soundvelocity distribution along the propagation path, yielding
an estimated sound velocity distribution in the tissue. The
performance of the proposed method was evaluated throughexperiments.

Methods
Estimation of sound velocity using time delay
The sound velocity distribution is estimated from the timedelay set of the RF signals received at multiple elements ofan ultrasonic probe, as follows. Let us assume that there are
N (n = 1, 2,…, N) regions with different sound velocities {cn}in the beam direction. Rn is defined as the sound-velocityestimation region where the depth is ln ; cn in Rn is assumed
to be constant.
An array probe transmits an ultrasound beam with anincidence angle of 𝜃1 to a scatterer in the shallowest region
R1, as shown in Fig. 1, where r1 is the distance from the central element to the strong target, and 𝜃1 is the direction angleof the target with respect to the central element ( xk = 0 ).
First, for the RF signal received at the central element, thearrival time of the echo returning from the strong target
is determined. Second, the time delay of the echo at eachelement, with respect to the arrival time of the echo at the

Journal of Medical Ultrasonics (2019) 46:27–33

central element, is obtained using the cross-correlation function. The time delay T(xk ) at the k th element is

T1 (xk ) ={}√1
r1 + xk2 − 2r1 sin 𝜃1 × xk + r12 ,c1()−KΔx ≤ xk ≤ KΔx ,(1)

where Δx is the element pitch, xk = kΔx is the lateral position of the k th element, and (2K + 1) is the number of receiving elements (k = −K, … , K). As the forward path from thecentral element to the target coincides with its backward path(from the target to the central element), i.e., T1 (0) = 2r1 ∕c1,
the time delay TR1 (xk ) at the k th element with respect to thecentral element is( )( ) T(0)TR1 xk = T xk −
2
√
1
2
=
xk − 2r1 sin 𝜃1 × xk + r12 .
c

(2)

( )
Similarly, the time delay TRn xk from the scatterer in the
nth region is
�
� n−1
∑
� li
rn cos 𝜃n − n−1
l
� �
i=1 i
+
TRn xk =
c
cn
i=1 i
(3)
�
1
2
2
×
xk − 2rn sin 𝜃n × xk + rn ,rn cos 𝜃nwhere {rn } and {𝜃n } are the distance and beam angle fromthe central element to the target point in the region {Rn},respectively. Here, the influence of refraction when the signal is incident on regions with different sound velocities
is ignored. In addition, the influence of the heterogeneous
sound velocity distribution in the lateral direction is also
ignored. For example, in the liver, there is uniform soundvelocity in the lateral direction in almost every case. It isalso possible to suppress the influence by confirming anddiscriminating the inhomogeneity region using the ultrasound B-mode image.
( )The square of Eq. (3) becomes aquadratic function fn xk
� �� �
fn xk = TR2 xk
n
�2
� n−1
∑
� li
rn cos 𝜃n − n−1
l
i=1 i
+
=
c
cn
i=1 i
(4)
� 2
�
1
x − 2rn sin 𝜃n × xk + rn2
× 2
rn cos2 𝜃n k

= a2 xk2 + a1 xk + a0 .
The( )coefficients a2 , a1 , and a0 are determined by fittingTR2 xk measured for the scatterer in each region to the quadn

Fig. 1  Schematic of the system setup with an ultrasonic probe andtarget scatterer

13

ratic function of Eq. (4) using the weighted least square method;{cn}, {rn}, and {θn} can be obtained using a2, a1, and a0.

Journal of Medical Ultrasonics (2019) 46:27–33

29

P is located at a distance r from the center of the ultrasonicprobe, in the sound-velocity-estimation region Rn ; the local
sound velocity estimated in that region is denoted as cn.
Figure 3 shows a schematic of the ultrasonic propagation
path. The propagation time tk (P) of the ultrasonic wave transmitted at x = 0, scattered at the point P, and received at theelement k with a lateral position xk is calculated as follows.
The length in the region ­Rn of the forward path from the∑central element to the point P is r − n−1
l , while that of thei=1 i
return
path
from
the
point
P
to
the
element k is��
∑n−1 � 2
r − i=1 li
r + xk2 ∕r.
In the sound-velocity-estimation{ } regions R1 to Rn−1, theli , while those of the return
lengths of the forward
path
are
( )√
2
2
r + xk . The ultrasonic propagation time inpath are li ∕rFig. 2  Measurement grids aligned vertically in an ROI to estimate thesound velocity distribution

This estimation process is applied to selected strong targets in each ROI. The propagation times in Eq. (4) include
information not only about the sound velocity in the nth
region but also about those up to the (n − 1)th region. Therefore, the sound velocities from the first to the nth region canbe simultaneously estimated by fitting the propagation times
from the scatterers in regions one to n to their propagationtime theoretical equations. The amplitude of the received
signal is large at the element above the target, and decreases
with the increase of the distance from the element, owing
to the different distance and direction from the scatterer toeach element of the ultrasonic probe. A weighting function
determined by the amplitude of the received signal at eachelement was applied to the above fitting. Therefore, it is possible to estimate the sound velocity distribution along thebeam direction, as shown in Fig. 2.

Beam‑forming of the received signal using
the estimated sound velocity
The delay time calculated using the velocity distributions
obtained in “Estimation of sound velocity using time delay”
is added to the received signal at each element. The received
beams are then compounded. When the amplitude A(P) of
the RF signal is estimated by forming the received beamat the point of interest P, the delay time to be applied to
the received signal at each element is calculated as follows.
Assuming that a linear probe is used, the ultrasonic wave istransmitted vertically to the probe surface, and the centerof the transmitted beam (x = 0) is set as the origin in the
direction parallel to the probe surface. The point of interest

each region can be obtained by dividing the propagation pathby the local sound velocity of the region. tk (P) can beobtained by adding up the round-trip propagation times fromthe probe surface to the region Rn:��li � 22l +tk (P) =r + xk
c i r
i=1 i
��
�
�
∑
n−1
�
r − n−1
l�
1
i=1 i
r−
li +
+
r2 + xk2
cn
r
(5)
i=1
�
�⎛
� n−1
⎞
∑
r2 + xk2 ⎟� lir − n−1li=1 i ⎜
+
=
⎜1 +
⎟.
c
cn
r
i=1 i
⎜
⎟
⎝
⎠
n−1
�
1

If the sound velocities are equal in{ all}regions, i.e., ci = c0,
Eq. (5) corresponds to Eq. (1) with 𝜃n = 0. Assuming theultrasound is transmitted at t = 0, and the received RF signal
at the element k is gk (t), the amplitude corresponding to thepoint P of the RF signal after the formation of the reception
beam is gk {tk (P)}. By adding them for all elements of theprobe, the reception beam RF signal can be formed. Therefore, A(P) can be expressed as

A(P) =

K∑{}
gk tk (P) .

(6)

k=−K

Results and discussion
We demonstrated the usefulness of the proposed method
through a basic experiment using a phantom including
strong targets. Figure 4 shows the experimental system. An
ultrasonic diagnostic unit (ProSound α10, Aloka) with a linear probe (UST-5412), whose center frequency is 10 MHz,was used with a transmission frequency of 7.5 MHz. The
probe contained 192 elements, with a pitch of 0.2 mm. A

13

30

Journal of Medical Ultrasonics (2019) 46:27–33

Fig. 4  Experimental setup of the phantom experiment and B-mode
image

Fig. 3  Schematic of the ultrasonic propagation path

total of 96 elements around the center were used to transmit
and receive signals. The sampling frequency to obtain the
RF signal in each element was 40 MHz. The temperature of
the water was 16 °C, and the sound velocity was 1469 m/s
in the water [24]. The sound velocity was 1540 m/s in thesilicone phantom.
Two molybdenum wires with a diameter of 50 µm were
installed in the water. Two nylon wires with a diameter of
100 µm were embedded in the silicone phantom to serve as
target points. Four estimation regions were provided; thelengths l1, l2, l3, and l4 in the depth direction of the gridwere 16.0, 17.0, 10.0, and 10.0 mm, respectively. The region
size was assumed to be judged from the B-mode image
for the target organ. Therefore, the length of each region
was selected as above. As a linear probe was used in thisexperiment, the beam angles {θn}, introduced above, wereassumed to be 0◦. Four focused beams with different focus
points set at the targets were transmitted. The measurements
were repeated eight times.
The delay-time distributions from the target at eachregion and parabolic approximations are shown with blue
dots and red curves in Fig. 5, respectively. The maximum
differences of the propagation times from the target in the
deep region are smaller than those in the shallow region. The
propagation times could not be obtained as a smooth curve
for a part of the results from the target in the deep region.
The estimated sound velocity and distance to the target
in each region are compared with the true values in Fig. 6.
Table 1 summarizes the local sound velocity in each region,
its estimation error, and standard deviation. Table 2 shows
the depth to the target. The estimated local sound velocity
in each region was close to the true value. The estimatedsound velocities of healthy persons and patients with hepatic

13

steatosis were significantly different by about 10–40 m/s
in the liver [25]. Therefore, we aimed for an accuracy of± 10 m/s in sound velocity estimation by this proposed
method. The error of 10 m/s did not significantly affect theresolution of the received beam. The standard deviations in
the deeper region were larger than those in the shallower
region, owing to the following three possible factors. The
first factor is attributed to the maximum difference in propagation times, which is small in the deep region, as shown in
Fig. 5, leading to a worse estimation accuracy. The second
factor is attributed to the decreases in the amplitudes of thereflected signals and resultant reductions in the signal-tonoise ratio owing to the large attenuation along the large

Fig. 5  Delay-time distributions from the target at each region (bluedots) and parabolic approximations (red curves)

Journal of Medical Ultrasonics (2019) 46:27–33

31

Fig. 7  Lateral beam profiles obtained using the a assumed soundvelocity of 1540 m/s, b assumed sound velocity of 1469.0 m/s (soundvelocity in water), and c estimated sound velocity

Fig. 6  Estimation results of the sound velocity distribution
Table 1  Estimated local sound velocities in each region, estimation
errors, and standard deviations in the phantom experiment
RegionEstimated (m/s)

True (m/s)

Estimation
error (%)

Standard
deviation
(m/s)

R1
R2
R3
R4

1469.91470.21531.41547.51469.41469.41540.01540.0

0.034
0.052
0.556
0.485

0.22
0.31
4.76
5.76

Table 2  Estimated depths in the phantom experiment
RegionEstimated (mm)True (mm)

R1
R2
R3
R4

16.7
25.9
39.7
53.4

16.0
25.0
39.0
53.0

propagation distance from the deep region. The third factor is attributed to the parabolic approximation in the deep
region, which is performed using the local sound velocities obtained in the shallow region; a slight difference fromthe true value in the shallow region would affect the soundvelocity estimation in the deeper region.
Furthermore, the received beam was formed using theestimated sound velocities. Figure 7 shows the lateral profile of the received beam composed of the signals from thewires in the third region. A sharp peak was obtained, and thepeak amplitude was improved by 3.1 dB and 1.3 dB, compared with the profiles beam-formed with the assumed soundvelocities of 1540 m/s (widely used in ultrasonic diagnostic
units) and 1469.0 m/s (sound velocity in water), respectively.

The nominal sound velocity of the silicone phantom (CIRS
054 GS) from the product information was treated as thetrue value. The full-width-at-half-maximum of the amplitude was also improved to 0.61 mm, from 0.92 mm with
1540.0 m/s and 0.76 mm with 1469.0 m/s. Regarding theimaging, tomograms obtained with the assumed velocitiesof 1540 m/s and 1469.0 m/s, and those with the estimatedvelocities, are shown in Fig. 8. Figure 9 shows the envelope
amplitude of the beam direction formed using the assumed
and estimated sound velocities. The results showed that thepeak was slightly sharp in beam-forming using the estimated sound velocities. From these results, we confirmed
improvements in the resolution in both the lateral and beam
directions on the received beam using the sound velocities
estimated by the proposed method. The blurring of the wire
image became smaller using the estimated sound velocity,and the improvement in the image quality of the tomogram
was confirmed by the proposed method.
In the conventional method, the average sound velocity
from the ultrasonic probe to the ROI is estimated to obtaina high-resolution ultrasonic tomographic image in the ROI.
On the other hand, the proposed method estimates the soundvelocity in each region along the beam direction from theprobe. As a result, it is possible to realize a high-resolution
ultrasonic tomographic image in the entire region by beamforming of the received signal. For example, if the soundvelocity distribution inside the organ is obtained, it is possible to estimate the local fat content and others, and early
detection of the lesion can be expected.
In the phantom experiment, the sound velocities estimated by the proposed method were close to the true values
as well as those estimated by the conventional method. In
addition, since the procedure of this method is composedof only (1) calculation of the cross-correlation function
between the elements of the received signals and (2) parabolic approximation by the least squares method, the computational effort is very small. On the other hand, in the

13

32

Journal of Medical Ultrasonics (2019) 46:27–33

phantom and wires were used as objects in the experiment.
Clutter suppression [26] would be necessary in in vivo application and should be studied in the future. In future studies,
we aim to investigate estimation methods using the planeboundary such as artery and estimate the sound speed usingreflection from strong scatterers present in living bodies.
There are high-intensity scatterers such as microvessels in
the B-mode image of the liver [27]. We confirmed the propagation time distribution from a strong target in the liver bymeasuring with the same experimental setting in the present
study. In the future, estimation of the sound velocity distribution in the liver by this method will be examined. Afterconfirming the usefulness of this method in the liver, weplan to expand the application to the other organs including
the breast.
Fig. 8  Images obtained by the conventional method with velocities ofa 1540 m/s and b 1496.0 m/s, and c estimated sound velocity

Compliance with ethical standards
Ethical statement This article does not contain studies with human or
animal subjects performed by the authors.
Conflict of interest The authors have no conflicts of interest with regard to the presented research.

References

Fig. 9  In the beam direction, beam profiles obtained using the aassumed sound velocity of 1540 m/s, b assumed sound velocity of1469.0 m/s (sound velocity in water), and c estimated sound velocity

conventional method, for example, the method [12, 21, 22]for estimating the optimum sound velocity using the focusing quality requires more computational cost because thefocusing quality has to be calculated for the beam-formed
data while changing the sound velocity from 1400 m/s to1600 m/s in steps of 5 m/s.

Conclusions
We proposed a method to estimate the sound velocity anddepth of a target in a ROI using the time delays of the signals
received by each element of the array probe. In the phantom
experiment, the local sound velocity in each region could be
estimated, which was close to the true value. The usefulnessof the proposed method was confirmed by improvement in
the blurring of the wire image. In this study, a stationary

13

1. Nagai Y, Hasegawa H, Kanai H. Improvement of accuracy
in ultrasonic measurement of luminal surface roughness ofcarotid arterial wall by deconvolution filtering. Jpn J Appl Phys.
2014;53:07KF19-1–9.
2. Miyachi Y, Hasegawa H, Kanai H. Automated detection ofarterial wall boundaries based on correlation between adjacent receive scan lines for elasticity imaging. Jpn J Appl Phys.
2015;54:07HF18-1–-11.
3. Kurokawa Y, Taki H, Yashiro S, et al. Estimation of size of red
blood cell aggregates using backscattering property of highfrequency ultrasound: in vivo evaluation. Jpn J Appl Phys.
2016;55:07KF12-1–8.
4. Sakai Y, Taki H, Kanai H. Accurate evaluation of viscoelasticity
of radial artery wall during flow-mediated dilation in ultrasoundmeasurement. Jpn J Appl Phys. 2016;55:07KF11-1–6.
5. Mochizuki Y, Taki H, Kanai H. Three-dimensional visualizationof shear wave propagation generated by dual acoustic radiationpressure. Jpn J Appl Phys. 2016;55:07KF13-1–5.
6. Tanaka M, Sakamoto T, Sugawara S, et al. Blood flow structureand dynamics, and ejection mechanism in the left ventricle: analysis using echo-dynamography. J Cardiol. 2008;52:86–101.
7. Chang JH, Raphael DT, Zhang YP, et al. Proof of concept:
in vitro measurement of correlation between radiodensity andultrasound echo response of ovine vertebral bodies. Ultrasonics.
2010;51:253–7.
8. Anderson ME, McKeag MS, Trahey GE. The impact of soundspeed errors on medical ultrasound imaging. J Acoust Soc Am.
2000;107:3540–8.
9. Goss SA, Johnston RL, Dunn F. Comprehensive compilation of
empirical ultrasonic properties of mammalian tissues. J Acoust
Soc Am. 1978;64:423–57.

Journal of Medical Ultrasonics (2019) 46:27–33
10. Goss SA, Johnston RL, Dunn F. Compilation of empirical ultrasonic properties of mammalian tissues. II. J Acoust Soc Am.
1980;68:93–108.
11. Tranquart F, Grenier N, Eder V, et al. Clinical use of ultrasoundtissue harmonic imaging. Ultrasound Med Biol. 1999;25:889–94.
12. Yoon C, Seo H, Lee Y, et al. Optimal sound speed estimation
using modified nonlinear anisotropic diffusion to improve spatial
resolution in ultrasound imaging. IEEE Trans Ultrason Ferroelectr
Freq Control. 2012;59:905–14.
13. Nakashima K, Endo T, Ikedo Y, et al. Importance of changes
in images by controlling the speed of the sound received during
B mode ultrasonography examination of the breast. J Jpn Assoc
Breast Cancer Screen. 2008;16:179–89.
14. Yang JN, Murphy AAD, Madsen EL, et al. A method for in vitro
mapping of ultrasonic speed and density in breast tissue. Ultrason
Imaging. 1991;13:91–109.
15. Bamber JC, Hill CR. Acoustic properties of normal and cancerous
human liver—I. Dependence on pathological condition. Ultrasound Med Biol. 1981;7:121–33.
16. Bamber JC, Hill CR, King JA. Acoustic properties of normal andcancerous human liver—II. Dependence on tissue structure. Ultrasound Med Biol. 1981;7:135–44.
17. Flax SW, O’Donnell M. Phase-aberration correction using signals
from point reflectors and diffuse scatterers: basic principles. IEEE
Trans Ultrason Ferroelectr Freq Control. 1988;35:758–67.
18. O’Donnell M, Flax SW. Phase-aberration correction using signals
from point reflectors and diffuse scatterers: measurements. IEEE
Trans Ultrason Ferroelectr Freq Control. 1988;35:768–74.

33
19. Nock L, Trahey GE, Smith SW. Phase aberration correction in
medical ultrasound using speckle brightness as a quality factor. J
Acoust Soc Am. 1989;85:1819–33.
20. Schneider FK, Yoo YM, Agarwal A, et al. New demodulation filter in digital phase rotation beamforming. Ultrasonics.
2006;44:265–71.
21. Napolitano D, Chou CH, McLaughlin G, et al. Sound speed correction in ultrasound imaging. Ultrasonics. 2006;44:e43–6.
22. Yoon C, Lee Y, Chang JH, Song TK. In vitro estimation of meansound speed based on minimum average phase variance in medical
ultrasound imaging. Ultrasonics. 2011;51:795–802.
23. Kondo M, Takamizawa K, Hirata M, et al. An evalution of anin vivo local sound speed estimation technique by the crossd beam
method. Ultrasound Med Biol. 1990;16:65–72.
24. Kroebel W, Mahrt KH. Recent results of absolute sound velocity
measurements in pure water and sea water at atmospheric pressure. Acustica. 1976;35:154–64.
25. Imbault M, Faccinetto A, Osmanski BF, et al. Robust sound speed
estimation for ultrasound based hepatic steatosis assessment. Phys
Med Biol. 2017;62:3582–98.
26. Jaeger M, Held G, Peeters S, et al. Computed ultrasound tomography in echo mode for imaging speed of sound using pulseecho sonography: proof of principle. Ultrasound Med Biol.
2015;41:235–41.
27. Mori S, Hirata S, Yamaguchi T, et al. Probability image oftissue characteristics for liver fibrosis using multi-Rayleighmodel with removal of nonspeckle signals. Jpn J Appl Phys.
2015;54:07HF20-11–8.

13