4 aha-a 3.7. 1

1"“? ,

a ..

3;. .....%t.
.% anagfimwfi ‘.
«a

...%

.4

.1}.

11:91:
,9: ‘15

. 1, 4-2....
x .. . ..
a!

>~x§fl1¥a~

.-\’

 

 

us-m nvm-‘e-nw-«qus

El

’1
new

{If 5..., I

51.1...1? , . ‘ .. ‘ ‘ . LI... .
Mara 3 . , , . . . . .

7..

 

r

 

 

Tums LIBRARY
2. Michigan State
lb or? University

 

 

This is to certify that the
dissertation entitled

ULTRASOUND PHASED ARRAY SIMULATIONS FOR
HYPERTHERMIA

presented by

XIAOZHENG ZENG

has been accepted towards fulfillment
of the requirements for the

 

 

 

 

Ph.D. degree in Electrical Enflweering
QM 0. WCM
Mafir Professor’s Signétdre
Nov, 7! 200?
Date

MSU is an aflinnative-action, equal-opportunity employer

To AVOID FINES return on or before date due.

l
1 PLACE IN RETURN BOX to remove this checkout from your record.
|
' MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

5/08 KzlProjIAoclLPresICIRClDateDue.indd

ULTRASOUND PHASED ARRAY SIMULATIONS FOR
HYPERTHERMIA

By

Xiaozheng Zeng

A DISSERTATION

Submitted to
Michigan State University.
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Electrical Engineering

2008

ABSTRACT

ULTRASOUND PHASED ARRAY SIMULATIONS FOR
HYPERTHERMIA

By
Xiaozheng Zeng

Hyperthermia uses elevated temperature to achieve a therapeutic benefit in cancer
treatments. Ultrasound is one of the most promising modalities in hyperthermia.
The numerical simulations of the pressure, power, and temperature fields generated
by ultrasound phased arrays are explored in this thesis.

Angular spectrum approach rapidly and accurately simulates the pressure fields
for ultrasound phased arrays. Angular spectrum simulations with the spatial propa-
gator produce accurate results in the central portion of the computational grid while
significant errors are produced near the edge due to the finite extent of the window
applied to the spatial propagator. Angular spectrum simulations with the spectral
propagator are more erroneous due to the undersampling in the spatial frequency do-
main. However, the error is greatly reduced by angular restriction in non—attenuating
media. The spatial and spectral propagators achieve similar accuracy in attenuating
media or for apodized sources. Angular spectrum simulations with pressure sources
are more accurate than those with normal velocity sources. The pressure source
plane should be larger than the array aperture and at least one wavelength away
from the array. While the above criteria are established for homogeneous and linear
simulations, the angular spectrum approach can also simulate nonlinear harmonic
propagations in layered tissue media. Results show that the absorption of these non-
linear harmonics is important to high temperature therapy, but is negligible to mild
temperature hyperthermia. The pressure fields simulated with the angular spectrum

approach are used as inputs to the bio-heat transfer equation to model temperatures.

The power and temperature distributions are simulated in tumor-tissue models, where
the objective is a relatively uniform temperature between 41 — 43°C in tumor and
minimized temperature rises in normal tissue. The power fields are produced either
with single-focus spot scans optimized by a thermal inverse method or with multiple-
focus scans optimized by the waveform diversity method. Results show that the
single-focus spot scan with the thermal inverse optimization achieves good results
for small tumors but produces excessive intervening tissue heating for large tumors.
The multiple-focus scan with the waveform diversity optimization achieves superior
results in heating large deep-seated tumors without inducing excessive intervening

tissue heating.

© 2008

Xiaozheng Zeng
All Rights Reserved

To Kang

ACKNOWLEDGMENTS

First and foremost, I would like to express my deepest gratitude and respect to
my Ph.D. adviser Dr. Robert J. McGough - a man of great wisdom and generosity.
He guided me into the world of therapeutic ultrasound, and advised me with great
insight during my exploration. I have benefited tremendously from his inspiration,
encouragement, and many helpful discussions. Without his advice and help, this
thesis will not be possible.

I am also deeply grateful to my Ph.D. committee members Dr. Edward Rothwell,
Dr. Leo Kempel, and Dr. Mei Zhuang. I thank Dr. Rothwell in particular for
sharing his expertise in angular spectrum approach. I own my sincere gratitude
to Dr. Shanker Balasubramaniam, who taught me three classes on computational
methods applied to acoustics and electromagnetics.

I thank my colleagues at the Biomedical Ultrasonics and Electromagnetics Labo-
ratory at MSU, Liyong Wu, James F. Kelly, Duo Chen, Don Chorman, Jason Biel,
Khawar Khurshid, Ruihua Ding, Josh Wong, for their friendship and help during
the period my doctoral study. In particular, I would like to thank Liyong Wu for
his collaboration on the hybrid thermal therapy project, and James F. Kelly for his
inspiring discussions which help me understanding the physics and mathematics in
acoustics.

Last but not least, I would like to thank my family. I thank my parents Jianming
Zeng and Zheng Song for their unconditional love and support, my husband Kang,
for his help and support to my research and my life. Without them I would have

never completed this achievement.

vi

TABLE OF CONTENTS

LIST OF TABLES .............................. x
LIST OF FIGURES ............................. xi
1 Introduction ................................ 1
1.1 Clinical aspects of hyperthermia ..................... 1
1.2 Engineering aspects of hyperthermia .................. 2
1.2.1 Ultrasound hyperthermia ..................... 3

1.2.2 Numerical modeling ....................... 4

1.3 Content of thesis ............................. 5

2 Pressure Simulation Methods for Transducers with Canonical

Shapes ................................... 7
2.1 Introduction ................................ 7
2.2 Circular piston .............................. 9
2.2.1 Point source superposition method for a circular piston . . . . 9
2.2.2 Spatial impulse response method for a circular piston ..... 9
2.2.3 Fast nearfield method for a circular piston ........... 10
2.2.4 Farfield approximation method for a circular piston ...... 11
2.3 Rectangular piston ............................ 11
2.3.1 Rectangular radiator method for a rectangular piston ..... 11
2.3.2 Spatial impulse response method for a rectangular piston . . . 12
2.3.3 Fast nearfield method for a rectangular piston ......... 13
2.4 Spherically focused piston ........................ 14

2.4.1 Point source superposition method for a spherically focused piston 14
2.4.2 Spatial impulse response method for a spherically focused piston 15

2.4.3 Fast nearfield method for a spherically focused piston ..... 17

3 Angular Spectrum Approach ......... h ............. 18
3. 1 Introduction ................................ 18

3.2 Theory ‘of the angular spectrum approach ............... 20
3.3 Influence of the propagators ....................... 22

vii

3.3.1 Spatial propagator ........................ 22

3.3.2 Spectral propagator ........................ 25
3.3.3 Spectral propagator with angular restriction .......... 28
3.3.4 Comparison between the spatial and spectral propagators . . . 32
3.4 Angular spectrum approach in attenuating medium .......... 40
3.5 Angular spectrum approach for apodized source ............ 44
3.5.1 Velocity source distribution ................... 44
3.5.2 Pressure source distribution ................... 47
3.6 Backward propagation .......................... 50
Ultrasound Phased Array Simulations ................ 56
4.1 Introduction ................................ 56
4.2 Phased array designs ........................... 57
4.2.1 Curvilinear array ......................... 57
4.2.2 Two dimensional planar array .................. 58
4.2.3 Spherically focused array ..................... 58
4.2.4 Spherical-section array ...................... 59
4.3 Ultrasound phased array beamforming ................. 59
4.3.1 Single focus ............................ 59
4.3.2 Sector-vortex focus ........................ 62
4.3.3 Mode scanning .......................... 63
4.4 Pressure simulations for phased arrays using integral approaches . . . 69
4.4.1 Phased array models ....................... 69
4.4.2 Computational coordinate system ................ 69
4.4.3 Simulation results ......................... 70
4.4.4 Computation time and numerical error ............. 72

Angular Spectrum Approach for Ultrasound Phased Array Simula-

tions in Linear and Homogeneous Media ............... 76
5. 1 Introduction ................................ 76
5.2 Comparison between the pressure and velocity sources ........ 77
5.3 Optimal parameters for the source field ................. 81
5.4 Evaluation of input and output errors .................. 86

viii

5.5 Computation time and numerical error with the angular spectrum ap—
proach ................................... 90

6 Angular Spectrum Approach for Ultrasound Phased Array Simula-

tions in Layered and Nonlinear Media ................ 92
6.1 Introduction ................................ 92
6.2 Acoustic simulations in layered and nonlinear media .......... 95
6.2.1 Angular spectrum approach in layered media .......... 95
6.2.2 Angular spectrum approach in nonlinear media ........ 97
6.3 Power and temperature simulations in layered and nonlinear media . 103
6.3.1 Acoustic intensity and power deposition ............ 103
6.3.2 Bio-heat transfer equation .................... 106
6.3.3 Temperature simulation in a layered tissue model ....... 108
6.3.4 The effect of nonlinear ultrasound propagation on the simulated
temperature ............................ 1 13

7 Power and Temperature Optimization in Ultrasound Hyperthermia 122

7.1 Introduction ................................ 122
7.2 Single focus spot scanning ........................ 124
7.3 Multiple focus mode scanning ...................... 133
7.4 Power optimization with waveform diversity method .......... 135
7.4.1 Introduction ............................ 135

7.4.2 Formulation ............................ 138

7.4.3 Simulation for a deep seated spherical tumor .......... 140

7.4.4 Simulation for a breast model with an embedded tumor . . . . 146

8 Conclusion ................................. 152
8.1 Summary ................................. 152
BIBLIOGRAPHY .............................. 155

ix

2.1

4.1

5.1

LIST OF TABLES

Combinations of 1,, j based on the location of the observation point
(3:,y, z). ..................................

The computation time and the root mean squared error (RMSE) for the
simulated pressures generated by (a) a planar array with 1024 square el-
ements, (b) a spherical-section array with 1444 square elements, and (c)
a spherically focused array with 1441 circular elements. The pressures
are simulated with two methods: the fast nearfield method (FNM)
and the point source superposition of the Rayleigh-Sommerfeld inte-
gral (Rayleigh-Sommerfeld). The simulated pressures are compared
with the reference pressure computed with the spatial impulse response
method. ..................................

The computation time and the root mean squared error (RMSE) for the
simulated pressure generated by (a) a planar array with 1024 square e1-
ement, (b) a spherical-section array with 1444 square elements, and (c)
a spherically focused array with 1441 circular elements. The pressures
are simulated with the angular spectrum approach using the particle
velocity source (ASA with U0) and the pressure source (ASA with
P0). The simulated fields are compared with the reference pressure
computed with the spatial impulse response method. .........

13

75

91

6.1 The acoustic and thermal properties of water and some biological tissues. 108

2.1

3.1

3.2

3.3

3.4

3.5

LIST OF FIGURES

The geometry of a spherically focused piston evaluated in a cross sec-
tion passing the z axis ...........................

The D x D source plane consisting of a non-zero normal particle ve-
locity distribution in a 20. x 2a square area on the piston surface. The
remaining area is filled with zeros .....................

Two dimensional cross section of the three dimensional reference pres-
sure generated by a 3cm x 3cm square piston in non-attenuating media.
The excitation frequency is 1MHz, and the normal particle velocity dis-
tribution is uniform across the piston surface. The reference pressure is
computed with the fast nearfield method. The result is normalized to
the maximum pressure amplitude computed in the three dimensional
volume. ..................................

Simulated pressure generated by a 30m x 3cm square piston in non—
attenuating media computed by the angular spectrum approach, where
(a) the spatial propagator is truncated with a 9cm x 9cm window and
evaluated in a 301 x 301 grid for FFT calculations, and (b) the spatial
propagator is truncated with a 9cm x 9cm window and zero-padded to a
512 x 512 grid for FFT calculations. The sample spacing is (5 = 0.075cm
in both (a) and (b) .............................

Simulated pressure generated by a 3cm x 3cm square piston in non-
attenuating media computed by the angular spectrum approach. where
(a) the spectral propagator is evaluated in a 301 x 301 grid for FFT
calculations, and (b) the spectral propagator is evaluated in a 512 x 512
grid for FFT calculations. The spatial frequency extent for the spectral
propagator is 2km = 27r/6, and the sample spacing is 6 = 0.075cm in
both (a) and (b). .............................

Simulated pressure generated by a 3cm x 3cm square piston in non-
attenuating media computed by the angular spectrum approach with
angular restriction, where (a) the spectral propagator is evaluated in a
301 x 301 grid for FFT calculations, and (b) the spectral propagator
is evaluated in a 512 x 512 grid for FFT calculations. The spatial
frequency extent for the spectral propagator is 2km = 27r/6, and the
sample spacing is 6 = 0.075cm in both (a) and (b). ..........

xi

16

23

24

26

29

3.6

3.7

3.8

3.9

3.10

3.11

Axial and transverse plots of the simulated pressure generated by a
3cm x 30m square piston in non-attenuating medium. Results are
shown for the reference computed by the fast nearfield method, the
angular spectrum approach using the spectral propagator without and
with angular restriction, and the spatial propagator. Results are com-
puted using N = 301 and D = 9cm ....................

The real and imaginary parts of the spectral propagator (without an-
gular restriction) and the FFT of the spatial propagator in a non-
attenuating medium evaluated at z = 15cm. The spectrum is plotted
along the kg; direction for Icy = 0. The spatial propagator is computed
with D = 9cm, and N = 301 is used for FFT calculations. The radi-
ally symmetric window applied to the spectral propagator shows the
angular threshold is at kc = 0.3917 k. ..................

The real and imaginary parts of the spatial propagator and the inverse
FFT of the spectral propagator (with angular restriction) in a non-
attenuating medium evaluated at z = 15cm. The pressure field is
plotted along the :c direction for y = 0 and normalized to the maximum
pressure amplitude along this line. The spatial propagator is computed
with D = 9cm, and N = 301 is used for the inverse FFT calculation.
The angular threshold for the spectral propagator is kc = 0.391719. . .

Root mean squared errors for the pressure generated by a 3cm x 3cm
square piston excited with uniform normal particle velocity and eval-
uated in a non-attenuating medium. The errors are evaluated in three
dimensional volumes where the lateral dimensions are (a) D x D and

(b)LwaithL=D—2a. ......................

A cross section of the attenuation term S(k$, kg, 2) along the k3, direc-
tion for kg = 0. The results are shown for z = 0.6cm (solid line) and
z = 15cm (dash-dot line). ........................

The absolute value of the complex reference pressure generated by a
3cm x 3cm square piston in attenuating media (oz = 1dB/cm/MHz).
The excitation frequency is 1MHz. The reference pressure is computed
with the fast nearfield method. .....................

xii

33

35

36

39

41

3.12

3.13

3.14

3.15

3.16

3.17

3.18

Simulated pressure in an attenuating medium computed by the angular
spectrum approach using (a) the spectral propagator without angular
restriction, (b) the spectral propagator with angular restriction, and
(c) the spatial propagator. All fields are calculated in successive D x D
transverse planes (D = 9cm) and truncated to L x L (L = 6cm), where
L = D —- 2a .................................

Root mean squared errors evaluated for the pressure generated in at-
tenuating media (a = 1dB/cm/MHz) by a 3cm x 3cm uniformly excited
square piston. The errors are evaluated in three dimensional volumes
where the lateral dimensions are L X L with L = D — 2a. The markers
on the curves indicate the corresponding results shown in Figure 3.12.

A square piston (illustrated by the rectangle at the bottom) and the
particle velocity distribution on the piston surface for (a) uniformly
excitation and (b) quadratically apodization (indicated by the mesh
on the top). ................................

The reference pressure calculated with the Rayleigh—Sommerfeld inte—
gral in a non-attenuating medium, where the pressure is generated by
a 3cm x 3cm square piston with an apodized normal particle velocity
distribution. The excitation frequency is 1MHz. ............

Simulated pressure in a non-attenuating medium generated by a
3cm x 3cm square piston with an apodized normal particle velocity
distribution. The pressure is computed by the angular spectrum ap-
proach using (a) the spectral propagator without angular restriction,
(b) the spectral propagator with angular restriction, and (c) the spa-
tial propagator. All fields are calculated in successive D x D transverse

43

44

45

47

planes (D = 9cm) and truncated to L x L (L = 6cm), where L = D — 2a. 48

Root mean squared errors for the pressure generated in a non-
attenuating medium by a 3cm X 30m square piston with an apodized
normal particle velocity distribution. The errors are evaluated in three

dimensional volumes where the lateral dimensions are L x L with
L = D — 2a. ...............................

The spectrum of the uniform velocity source and the apodized velocity
source along along the kg: direction evaluated at ky = 0 .........

xiii

49

3.19 Root mean squared errors for the pressure generated in non-attenuating
media by a 30m x 3cm square piston calculated with velocity and pres-
sure sources. The errors are evaluated in successive planes in different
places in the z direction. The pressure is computed by the angular
spectrum approach using (a) the spectral propagator without angular
restriction, (b) the spectral propagator with angular restriction, and
(c) the spatial propagator. All fields are calculated in successive D x D
transverse planes (D = 9cm) and truncated to L x L (L = 60m), where
L = D — 2a .................................

3.20 The pressure field in my plane at z = 15cm generated by a 3cm x 3cm
square piston in non-attenuating media. The excitation frequency is
1MHz. The pressure field is computed with the fast nearfield method
and the result is normalized to the maximum pressure amplitude in
the z = 15cm plane. ...........................

3.21 A cross section of the simulated pressure generated by a 3cm x 3cm
square piston in non-attenuating media. The 3D pressure is calculated
from a windowed pressure distribution at z = 15cm and propagated
forward and backward with the angular spectrum approach using (a)
the spectral propagator with angular restriction evaluated in a 301x 301
grid, and (b) the spatial propagator truncated by a 90m x 9cm window
and evaluated in a 301 x 301 grid. ...................

3.22 Thesimulated pressures along the :1: direction evaluated at y = z = 0.
The solid line is the reference pressure, the dashed line is the pres-
sure computed with the spectral propagator without angular restric-
tion (POHp), and the dash-dot line is the pressure computed with the
spatial propagator truncated by a 90m x 9cm window (pg (8) hp). The
source plane is located in the plane at z = 15cm, and the pressure is
backward propagated with the angular spectrum approach .......

4.1 A curvilinear array comprised of 51 rectangular elements ........
4.2 A two dimensional planar array comprised of 1024 square elements.

4.3 A spherically focused phased array comprised of 1441 circular elements:
(a) top view and (b) side view. .....................

4.4 A spherical—section array comprised of 1444 square elements: (a) top
view and (b) side view. .........................

xiv

51

52

53

55

57

58

4.5

4.6

4.7

4.8

4.9

4.10

4.11

4.12

Contours describing the pressure field generated by a 661-element
spherically focused array excited by sector-vortex phases with mode
number M=3. The array has an 80m aperture diameter and an 8cm
radius of curvature. The focal spots are located in the my plane at
z = 8cm ...................................

Illustration of a group of four symmetric elements in a two dimensional
planar phased array. ...........................

Illustration of a spherically focused phased array. The elements labeled
from 1 to 6 are symmetric with respect to the hexagonal axes.

The contours of the pressure field generated by a 661 element spheri-
cally focused array focused with the mode scanning method. The array
aperture diameter is 8cm, and the radius of curvature is 80m. The fo-
cal spots are located in the my plane at z = 8cm, and each focal spot
is located 0.50m from the z axis. The pressure is simulated with the
fast nearfield method. ..........................

The contours of the pressure field generated by a 661 element spheri-
cally focused array focused with the mode scanning method. The array
aperture diameter is 8cm, and the radius of curvature is 8cm. The fo-
cal spots are located in the my plane at z = 8cm, and each focal spot
is located 1cm from the z axis. The pressure is simulated with the fast
nearfield method. . . . ., ........................

Contour plot describing the pressure field generated by a 661 element
spherically focused array focused with the mode scanning method. The
aperture diameter of the array is 8cm, and the radius of curvature is
8cm. The focal spots are steered off axis in the my plane at z = 8cm,

64

65

66

67

and the center of the focal pattern is located at (m, y, z) = (0, 0.80m, 80m). 68

The simulated pressure generated by a planar array with 1024 square
elements (a) in the yz plane at m = 0 and (b) in the my plane at
z = 10cm. The pressure is normalized by the peak pressure amplitude
in the three dimensional volume. ....... h .............

The simulated pressure generated by a spherical-section array with
1444 square elements (a) in the yz plane at m = 0 and (b) in the my
plane at z = 12cm. The pressure is normalized by the peak pressure
amplitude in the three dimensional volume. ..............

XV

4.13 The simulated pressure generated by a spherically focused array with

5.1

5.2

5.3

5.4

1441 circular elements (a) in the yz plane at m = 0 and (b) in the my
plane at z = 12cm. The pressure is normalized by the peak pressure
amplitude in the three dimensional volume ................

The simulated axial pressures generated by a planar phased array
located in the my plane at z = Ocm and electronically focused at
(0,0,100m). The reference pressure calculated by the spatial impulse
response method is represented by the solid line, the pressure com-
puted with the angular spectrum approach using the velocity source is
represented by the dash—dot line, and the pressure computed with the
angular spectrum approach using the pressure source is represented by
the dashed line. .............................

Root mean squared errors evaluated in transverse planes for z rang-
ing from 4cm to 16cm. The errors are obtained from pressure fields
generated by the 1024 element planar phased array, which is electron-
ically focused at (O, 0, 10)cm. The pressure is calculated with the an—
gular spectrum approach using a velocity source plane (dash-dot line)
and a pressure source plane (dashed line). The velocity and pres-
sure source planes are both located at 2 == 0cm and truncated with a
16.2cm X 16.2cm square window. The source field is sampled at an in-
terval of (a) 5 = 0.075cm and (b) 5 = 0.0375cm. The FFT calculations
are evaluated in N X N grids with N = 512 and N = 1024 in (a) and
(b), respectively. .............................

The angular spectra of the particle velocity source and the pressure
source for the 1024 element planar array in Figure 4.2 along the kg;
direction for log = 0. The angular spectra are obtained by applying a
2D FFT to the particle velocity or pressure source at z = 0, where the
source is sampled with 6 = A / 2 (solid line) and A / 4 (solid-dot line).

Root mean squared errors evaluated in 3D volumes as a function of
the pressure source plane location 20. The errors are obtained from
pressure fields generated by the 1024 element planar phased array,
which is electronically focused at (0,0,10)cm. The size of the source
plane is (a) 7.8cm X 7.8cm and (b) 16.2cm X 16.2cm, which is sampled
at A / 2 and zero-padded to a 512 X 512 grid. The location of the source

74

78

80

82

plane ranges between 20 = 0 and 20 = 3cm with an increment of 0.150m. 84

xvi

5.5

5.6

5.7

5.8

Axial pressures simulated with the angular spectrum approach using
pressure source planes that are truncated by square windows of dif-
ferent sizes. The pressures are generated by the 1024 element planar
phased array, which is electronically focused at (0, 0, 10)cm. The refer-
ence pressure is indicated by the solid line, the pressure computed with
the angular spectrum approach using a 6cm X 6cm pressure source plane
is represented by the dash-dot line, and the axial pressure computed
with the angular spectrum approach using a 7.80m x 7.80m pressure
source plane is represented by the dashed line. The solid line and the
dashed line are nearly coincident, which indicates that L = 7.8cm is
sufficiently large for a 7 .31cm X 7.31cm phased array. .........

Root mean squared errors plotted as a function of the source plane
extent L. The errors are obtained from pressure fields generated by
the 1024 element planar phased array, which is electronically focused
at (0,0,10)cm. The pressure source plane for angular spectrum sim-
ulations is located at 20 = 0.15cm and truncated by L X L square
windows, where L ranges from 6cm to 20.4cm with an increment of

0.15cm. The two markers indicate the values of L shown in Figure 5.5.

Root mean squared errors in the source plane plotted as a function
of the number of source plane abscissas used in each direction. The
errors are obtained from pressure fields generated by the 1024 element
planar phased array, which is electronically focused at (0, O, 10)cm. The
source plane is located at 20 = 0.15cm and truncated by a 7.8cm X
7.8cm window. The source pressure is computed with the Rayleigh-
Sommerfeld integral using 2 X 2 to 10 X 10 abscissas and with the fast
nearfield method using 2 to 10 abscissas. ...............

Root mean squared errors for 3D pressure field outputs plotted as
a function of the number of source plane abscissas used in each di-
rection. The errors are obtained from pressure fields generated by
the 1024 element planar phased array, which is electronically focused
at (0,0,10)cm. The source plane pressures are calculated with the
Rayleigh-Sommerfeld integral using 2 X 2 to 10 X 10 abscissas and with
the fast nearfield method using 2 to 10 abscissas. The errors obtained
from both methods approach the same limiting value, but the fast
nearfield method achieves convergence with far fewer abscissas. . . . .

xvii

85

86

88

89

6.1

6.2

6.3

6.4

6.5

6.6

6.7

The axial pressures generated by a spherically focused transducer in
a medium with water and fatty tissue layers. The planar interface
between the water and fatty tissue layers is located at 6cm, 7cm, 8cm,

and 9cm. The pressure is calculated with the angular spectrum approach. 96

Simulated fundamental, second, and third harmonics along the line at
m = y = 0 generated by a plane circular transducer in water. The
transducer is excited by a 2.25MHz continuous wave input with a peak
amplitude of 100kPa. The solid line with the largest peak amplitude
represents the fundamental harmonic, the lines with the second and
third peaks represent the second and the third harmonics. ......

Simulated fundamental, second, and third harmonics generated by a
plane circular transducer in water. The transducer is excited by a
2.25MHz continuous wave input with a peak amplitude of 100kPa. The
pressures are plotted along the :1: direction for y = 0 at (a) z = 27.5cm
and (b) 2 = 500m. The solid line with the largest peak amplitude
represents the fundamental harmonic, the lines with the second and
third peaks represent the second and the third harmonics. ......

Nonlinear pressure waveforms generated by a spherically focused trans-
ducer with a 10cm aperture diameter and a 20cm radius of curvature
at the focus (a) in fatty tissue and (b) in water. ............

Simulated fundamental, second, and third harmonics along the line
at m = y = 0 generated by a spherically focused transducer with a
3W/cm2 input intensity. The simulation model includes a water layer
and a fatty tissue layer, and the interface is at z = 9cm. The solid line
with the largest peak amplitude represents the fundamental harmonic,
the lines with the second and third peaks represent the second and the
third harmonics ...............................

A model of the human abdomen that includes four soft tissue layers. .

The normalized axial pressure generated by a 2575 element spherically
focused phased array evaluated along the z axis where m = y = 0. The
pressure is simulated in a homogeneous fatty tissue medium (solid line)
and a model containing layers of water, skin, fat, muscle and viscera

(dash-dot line). ..............................

xviii

102

104

105

109

111

6.8

6.9

6.10

6.11

6.12

6.13

7.1

The simulated normalized power deposition and temperature rise gen-
erated by the 2575 element spherically focused phased array in the mz
plane at y = 0. The ultrasound bio-heat simulations are performed in
a model containing layers of water, skin, fat, muscle, and viscera.

The simulated axial pressure along the z axis at m = y = 0 generated
by (a) a 16cm X 16cm planar phased array and (b) a 4cm X 4cm plar
nar phased array. The dashed line represents the pressure simulated
with the linear model and the rest of the curves are harmonics sim-
ulated with the nonlinear model. The simulations are evaluated in a
homogeneous fatty tissue medium .....................

The simulated peak intensity and temperature rise as functions of the
input intensity produced by the 4cm X 4cm planar phased array. The
solid lines represent the results simulated with the linear model, and the
dotted lines are the results simulated with the nonlinear model. Figures
(a) and (b) contain the results from the steady state simulations, and
Figures (c) and (d) contain the results from transient simulations with
10 second exposures. ...........................

The temperature rise generated by the 4cm X 4cm planar array with
a 2.5W/cm2 input intensity using (a) the fundamental and (b) the
second harmonic obtained from nonlinear simulations. ........

The temperature rise generated by the 4cm X 4cm planar array with
a 2.5W/cm2 input intensity using the linear model (top half of the
figure) and the nonlinear model (bottom half of the figure) .......

The normalized difference between the peak temperatures obtained
from linear and nonlinear simulations for (a) steady state and (b) 10
second transient exposures. The difference is normalized by the maxi-
mum temperature rise obtained from the linear model and plotted as
a function of the linear temperature rise. ................

The simulated temperature rise in the yz plane evaluated at m = 0 pro—
duced by a 1444 element spherical-section array that is focused within
a 1.2cm diameter spherical tumor. The heating pattern is created by
scanning six focal spots repeatedly around the circumference of the
tumor. ..................................

xix

112

115

118

119

120

7.2

7.3

7.4

7.5

7.6

7.7

7.8

The simulated temperature rise produced by a 1444 element spherical-
section array in a spherical tumor with a 3cm diameter. The heating
pattern is created by spot scanning 172 focal points that are evenly
distributed in the tumor with uniform weighting. ...........

The distribution of control points within the spherical tumor and in
the normal tissue region for direct thermal inverse optimization. The
3cm spherical tumor is represented by the sphere, and the 1444 element
spherical-section phased array is located at the bottom of the figure. .

Simulated power deposition in the my plane evaluated at z = 15cm
produced by the 1444 element spherical-section array in a spherical
tumor with a 3cm diameter. The heating pattern is created by spot
scanning 172 foci within the tumor, where the power weight for each
focus is optimized with the direct thermal inverse method. . . .....

The simulated temperature rise produced by the 1444 element
spherical-section array in a spherical tumor with a 3cm diameter. The
heating pattern is created by spot scanning 172 foci within the tu-
mor, where the power weight for each scan is optimized with the direct
thermal inverse method. .........................

A transparent mesh that indicates the size and location of the spherical
tumor model and an isothermal surface that indicates the region that
achieves a minimum of 4°C temperature rise. The heating pattern is
created by spot scanning 172 foci within the tumor, where the power

128

129

130

131

weight for each focus is optimized with the direct thermal inverse method. 132

The power deposition in the my plane evaluated at z = 15cm pro-
duced by the 1444 element spherical-section array in a spherical tumor
with 3cm diameter. The heating pattern is created by mode scanning
172 foci within the tumor, where the power weight for each focus is
optimized with the direct thermal inverse method ............

The simulated temperature rise produced by the 1444 element
spherical—section array in a spherical tumor with a 3cm diameter. The
heating pattern is created by mode scanning 172 foci in the tumor,
where the power weight for each focus is optimized with the direct
thermal inverse method. .........................

135

7.9

7.10

7.11

7.12

7.13

7.14

7.15

A transparent mesh that indicates the size and location of the spherical
tumor model and an isothermal surface that indicates the region that
achieves a minimum of 4°C temperature rise. The heating pattern is
created by mode scanning 172 foci within the tumor, where the power

weight for each focus is optimized with the direct thermal inverse method. 137

The 1444 element spherical-section array, the 3cm spherical tumor and
the control points used in the waveform diversity optimization. One
quarter of the array elements and control points that are located in the
same quadrant of the coordinates are used in the optimization, that is,
a total of 361 array elements (painted in red), 43 tumor control points,
and 67 tissue control points. ......................

Simulated power deposition in the my plane evaluated at z = 15cm
produced by the 1444 element spherical-section array in a spherical
tumor with a 3cm diameter. The heating pattern is created by the
waveform diversity method combined with mode scanning. A total
of 361 array elements, 43 tumor control points, and 67 tissue control
points are used in the optimization ....................

Simulated temperature rise produced by the 1444 element spherical-
section array in a spherical tumor with a 30m diameter. The heating
pattern is created by the waveform diversity method combined with
mode scanning. 361 array elements, 43 tumor control points, and 67
tissue control points are defined for optimization. ...........

The 3cm spherical tumor model and the 4°C isothermal surface. The
heating pattern is created by the waveform diversity method combined
with mode scanning. 361 array elements, 43 tumor control points, and
67 tissue control points are defined for optimization ...........

The power deposition in (a) the mz plane evaluated at y = 0 and (b)
the yz plane evaluated at m = 0 produced by a 256 element spherical-
section array in a breast model with an embedded 3cm diameter spher-
ical tumor. The heating pattern is optimized by the waveform diversity
method. ..................................

The temperature rise in (a) the 1132 plane evaluated at y = 0 and (b)
the yz plane evaluated at m = 0 produced by a 256 element spherical-
section array in a breast model with an embedded 3cm diameter spher-
ical tumor. The heating pattern is optimized by the waveform diversity
method. ..................................

141

142

144

145

149

7.16 3D illustration of the 256 element spherical-section array, the breast
model, the embedded spherical tumor, and the 4°C isothermal surface.
The heating pattern is optimized with the waveform diversity method. 151

CHAPTER 1

Introduction

1.1 Clinical aspects of hyperthermia

Surgery, chemotherapy, radiation therapy, and hyperthermia are available options for
cancer treatment. Current clinical cancer treatments often involve undesirable side—
effects [1]. Surgery may cause infection, hemorrhage, and scarring. Radiation therapy
and chemotherapy heavily damage both normal cells and cancerous cells. Common
side effects of chemotherapy include anemia, fatigue, and hair loss. Hyperthermia
is a noninvasive therapy that can sometimes causes burns, especially at muscle/ fat
interfaces. However, there can be a synergistic effect between heat and radiation that
improves treatment outcome.

Hyperthermia uses high temperatures to eliminate cancerous tissues. Biological
studies show that 41 — 45°C is sufficient to induce a direct cytotoxic effect on cells,
including destruction of the cell membrane and cytoskeleton and disruption of the
cell metabolism and cell reproductive cycles [2]. Hyperthermia is especially effective
in solid and poorly vascularized tumors that are often resistant to radiation ther-
apy. Heat also increases the blood flow and the tissue oxygenation, enhances the
toxicity of certain drugs and their uptake by tissues, and inhibits the repair of cells
after radiation damage [3]. Therefore, hyperthermia is a direct treatment as well as

an adjuvant to chemotherapy and radiotherapy. Moreover, malignant cells are more

thermally sensitive than normal cells, and tumor cells suffer from chronic hypoxia, nu-
tritional deprivation and low pH [4]. These facts support the possibility of destroying
malignant tissue while minimizing the harm to healthy tissue. Higher temperatures
above 48°C and below 100°C can induce irreversible damage to critical cellular pro-
teins, tissue structural components, and vasculature, and these effects lead to tissue
destruction in a short time. Tissue ablation at these temperatures therefore provides

a minimally invasive alternative to traditional surgery [2].

1.2 Engineering aspects of hyperthermia

Ultrasound (US) and electromagnetic (EM) heating, including mircrowave (MW)
and radio-frequency (RF), are the most popular modalities in hyperthermia. Radio-
frequency current (0.1MHz to 100MHz) heats tissue by coupling electromagnetic en-
ergy capacitively or inductively [2]. Microwave devices are driven under higher fre-
quencies from 100MHz to 3GHz. Typical frequencies for clinical MW hyperthermia
systems are 433MHz and 915MHz [2]. The mechanisms of EM for hyperthermia in-
clude the friction between molecular dipoles during dielectric relaxation and the ionic
conduction current induced by ion oscillations. The absorption of electromagnetic
waves is proportional to the permittivity and conductivity of biological tissue [3].
Tissue that contains a high water content, for example, muscle or skin, has higher
dielectric parameters than fat and bone; therefore, higher power deposition is induced
in these water-based tissues. Clinical ultrasound frequencies range from 500kHz to
5MHz [2]. The heating mechanism of ultrasound is the microscopic friction during
periodical movement of particles [3]. The absorption of ultrasound in water is pro-
portional to the square of the ultrasound frequency, and the absorption in biological
tissues is roughly proportional to the ultrasound frequency [5].

Duo to the small value of the penetration depth, microwaves are ineffective when

applied to deep seated tumors. Microwaves are good for heating superficial regions

but not for deep seated targets. Ultrasound penetrates human tissues better than
microwaves. EM hyperthermia tends to heat large regions uniformly, while US hy-

perthermia is better at heating smaller regions.

1.2.1 Ultrasound hyperthermia

Basic clinical ultrasound hyperthermia systems contain RF signal generators, power
amplifiers, electronic matching circuits, water cooling systems, ultrasound applicators,
and temperature monitoring systems that provide feedback control. Over the past
three decades, hundreds of hyperthermia devices have been investigated theoretically,
several of these have been turned into prototypes in research facilities, and a few have
been commercialized by industry.

Early researchers attempted to perform hyperthermia with single disc transduc-
ers [6] or spherically focused transducers [7, 8]. Since the emergence of phased array
technology, ultrasound phased arrays become the state of the art in ultrasound hyper-
thermia. For superficial cancer treatment, a 16 transducer planar array Sonotherm
1000 (Labthermics Technologies Inc, Champaign, IL) applicator was developed [9]. A
scanning ultrasound reflector linear arrays system (SURLAS) was also developed, and
this system has been applied to small animal experiments [10, 11, 12]. Several phased
arrays have also been designed for hyperthermia in deep seated tumors, including a
cylindrical-section array developed by Ebbini and Cain [13], a spherical-section ar-
ray for deep regional hyperthermia [14], a 200 element quasi-random sparse array
for trans-skull brain therapy by Pernot et al [15], and concentric-ring and sector-
vortex arrays by Cain and Umemura [16, 17]. Recent developments in hyperthermia
systems focus on integrating flexible or site-specific multi-elements applicators with
more sophisticated patient treatment planning. New methods and technologies have
been developed to provide real-time monitoring and feedback control with the aid of

microwaves, ultrasound, and magnetic resonance imaging.

1.2.2 Numerical modeling

The precise control of power deposition and temperature distribution in both tumor
and healthy tissue is essential to hyperthermia. Computer simulations provide valu-
able predictions of the treatment results before a complicated hyperthermia system
is constructed. As most hyperthermia devices are at the experimental or clinical
trial stage, comparisons between computer simulations and experimental results are
extremely helpful in optimizing the design and operation of future clinical systems.
Power deposition patterns, thermal responses, and tissue parameters can be predicted
from computer modeling. These simulation results are essential in the assessment of
the effectiveness of hyperthermia.

Patient models are needed for computer simulations. These models are usually
constructed with patient anatomical images obtained from computed-tomography
(CT) or magnetic—resonance—imaging (MRI). Normal and cancerous tissues are iden-
tified from these images. Then, the number of focal points and the scan paths for
focused ultrasound are initially determined based on the geometry and the underly-
ing physics. The power deposition and the temperature response or the thermal dose
are then computed. The simulations determine if the treatment goal is attainable.
If not, certain parameters such like the scan path and the input power weights are
altered until optimal results are obtained. The speed of these calculations and the
optimizations is extremely important, as these simulations are usually repeated many
times.

In acoustic simulations, a linear and homogeneous medium is traditionally as-
sumed, where the linear wave equation is the governing equation for ultrasound
wave propagations. Various integral methods have been developed for the numer-
ical solutions to the linear wave equation, including the point source superposition
of the Rayleigh—Sommerfeld integral [18, 19], the spatial impulse response method
[20, 21, 22], and the fast nearfield method [23, 24, 25, 26]. Heterogeneity can be

handled with the finite difference time domain method (FDTD), the finite-element
method (FEM), and the pseudo-spectral (PS) method. The finite-element method is
especially useful for modeling heterogeneous tissue with irregular geometries. How-
ever, these methods are both time and memory consuming and are therefore imprac-
tical for large scale three dimensional simulations of hyperthermia applications.

Within the last two decades, nonlinear wave propagation has become more im-
portant in simulation of medical ultrasound. This is in part motivated by the per-
formance of harmonics in imaging, and in part due to the emergence of high in-
tensity focused ultrasound (HIFU) as a therapeutic modality. The nonlinear wave
propagation can be modeled by various nonlinear wave equations under different ap—
proximations, for example, the Burger’s equation, the Westervelt equation, and the
Khokhlov-Zabolotshaya—Kuznetsov (KZK) equation [27].

Thermal modeling for hyperthermia is conventionally performed with the Penne’s
bio—heat transfer equation (BHTE). BHTE effectively predicts the temperature re—
sponse to local heating subject to certain boundary conditions. The power deposition
produced by ultrasound provides the input to this inhomogeneous partial difl’erential
equation, and this is the quantity that must be optimized in order to achieve the
desired temperature distribution. Available optimization methods for this purpose
include the pseudo-inverse pattern synthesis method [28] and the direct thermal in-

verse method [29].

1.3 Content of thesis

This dissertation investigates fast pressure simulation methods for ultrasound phased
arrays and robust power and temperature optimization methods for hyperthermia
applications. Chapter 2 reviews several integral approaches that model the pressure
fields produced by transducers with canonical shapes. Chapter 3 introduces the an-

gular spectrum approach for pressure field calculations in a homogeneous and linear

medium. The spatial and spectral propagators are compared, and the artifact associ-
ated with each propagator is identified. The influence of the attenuation of the media
and the source distribution are evaluated. Chapter 4 introduces some configurations
of large ultrasound phased arrays and beamforming methods. The pressure fields
generated by these ultrasound phased arrays are computed with integral approaches
and the computation time and numerical errors for these simulations are compared.
Chapter 5 discusses the optimal parameters used by the angular spectrum approach
to simulate large phased arrays. Chapter 6 introduces pressure simulations with the
angular spectrum approach in a layered and nonlinear medium. The pressure, power,
and temperature simulations in this complex model are demonstrated and compared
with simulation results in a homogeneous medium. Chapter 7 evaluates several single-
focus and multiple-focus scanning strategies for heating deep seated tumors, where
the power field is optimized with a direct thermal inverse method and a waveform
diversity method. The thermal inverse optimization produces excessive intervening
tissue heating, but the waveform diversity optimization produces result that closely

matches the required temperature distribution. Chapter 8 concludes the thesis.

CHAPTER 2

Pressure Simulation Methods for

Transducers with Canonical Shapes

2. 1 Introduction

The acoustic pressure generated by a planar transducer embedded in a rigid and
infinite baffle can be modeled with the Rayleigh-Sommerfeld integral [18, 30],

jpc
p(ra t) = —

/\ ds', (2.1)

 

ej(wt — k [r — r’l)
[2.

Ir - 1‘’l
where p is the pressure, 8’ is the surface area of the transducer, 8’ is the integration
variable, w is the angular frequency of the excitation signal, /\ is the wavelength, k
is the wavenumber, c is the speed of sound, p is the density of the medium, u is the
normal velocity on the transducer surface, t represents time, r is the spatial coordinate
of an observation point, and r’ is the spatial coordinate of a point on the transducer
surface. Equation 2.1 evaluates the Green’s function solution to the homogeneous
wave equation under a rigid boundary condition. This equation describes the pressure
produced by a planar radiator for a single frequency continuous wave excitation.
The pressure simulations for hyperthermia applications conventionally employ
point source superposition of the Rayleigh-Sommerfeld integral [18, 30] or the rect—
angular radiator method [31]. The accuracy of these methods is extremely limited in

the nearfield of the source due to the numerical singularity in the denominator. Both

methods evaluate double integrals; therefore, the computation time quadratically in-
creases with respect to the number of subdivisions in each dimension. In ultrasound
imaging, the impulse response method is widely used. The impulse response method
evaluates a 1D convolution integral between the time derivative of the excitation sig—
nal and the impulse response of the radiator. The popular ultrasound simulation
program Field II [32] uses the impulse response method to compute the scattered
pulse-echo field from a computer generated tissue target. However, since transient
waves are used in medical ultrasound imaging, the impulse response is usually evalu-
ated in the time domain and transformed into the temporal frequency domain through
the Fourier transform. A large amount of time samples is required to accurately sam—
ple the impulse response function, so this computation is very time consuming. For
continuous wave (CW) excitations, this Fourier transform can be expressed analyti-
cally to reduce the computation time. The CW fast nearfield method is specialized
for continuous wave simulations. The fast nearfield method replaces the time variable
with a distance variable in the integration, removes the numerical singularity in the
denominator of the impulse response integrand, and employs an uniform integral that
applies everywhere in the spatial domain. The CW fast nearfield method has been
derived for circular transducers [23], rectangular transducers [24], triangular trans-
ducers [26], and spherical transducers. The fast mearfield method is accurate for both
nearfield simulations of single transducers and simulations of phased arrays.

This chapter contains a comprehensive review of the numerical methods that sim-
ulate the continuous wave pressure generated by transducers with canonical shapes.
The point source superposition of the Rayleigh-Sommerfeld integral, the impulse re-
sponse method, the rectangular radiator method, the fast nearfield method, and the
farfield expressions for circular, rectangular and spherical transducers are summarized

in this chapter.

2.2 Circular piston

2.2.1 Point source superposition method for a circular piston

For a circular radiator with uniform particle velocity distribution on the surface, the

Rayleigh-Sommerfeld integral is formulated as [25]

 

—jk\/r2 + z2 + a2 - 2racosz/2

2J-PC'U’O 'wt fa /7r 8
p 7', z,t = ——e] adodzp, (2.2)
( ) A O 0 \/T2+Z2+0'2—2TUCOS’I/J

 

 

where a is the radius of the circular transducer, no is the uniform normal particle
velocity on the piston surface, (0, 1,0) is the polar coordinate of a source point on the
transducer surface, and (r, z) is the cylindrical coordinate of an observation point.
Circular symmetry is exploited in Equation 2.2 by using 7‘ = W.

The point source superposition method [19] divides the transducer surface into a
number of radiating point sources. The total pressure field produced by the transducer
is approximated by the superposition of the fields produced by all point sources.
The two-dimensional integral in Equation 2.2 is evaluated with Riemann’s mid—point
quadrature rule, which is equivalent to dividing the circular transducer into sectors
and concentric rings and approximating each resulting section with a point source.
The numerical convergence of point source superposition is typically very slow in
the nearfield due to the difficulty of sampling the function near the singularity at

r=z=0.

2.2.2 Spatial impulse response method for a circular piston

The pressure field generated by a transducer driven uniformly by a single frequency
excitation is proportional to the Fourier transform of the spatial impulse response
of the transducer evaluated at that particular frequency [21, 33]. The analytical

expression for the pressure produced by a circular transducer calculated using the

spatial impulse response method is [20, 34]
(e—jkz _ e—jlcx/z2 + a2) 7. = 0

p(7', z,t) = pcquJWt jk11 + (e—jkz _ e—jk\/z2 + (r -— (1)2) 0 < r < a . (2.3)

jkll r _>_ a.

 

Two of the three expressions in Equation 2.3 depend on the numerical evaluation of

the integral

1 ,/(.+..)2..2 g2_22+.~2_.2 .
1 =_ -1 "de 2.4
1 «ffi_.)2+zzc°s ( 2r\/a_2—_z2 e ‘3’ ( )

where 6 = ct is the distance variable and c is the speed of sound. The impulse re-

 

 

 

sponse method computes the pressure field faster than the point source superposition
method because the single integral in Equation 2.4 converges more rapidly than the
double integral in Equation 2.2. However, the impulse response method encounters
some numerical artifacts at r z a for all values of z. Therefore, slow numerical con-
vergence is expected near 1” = 0.. These artifacts are caused by the infinite slopes in
the integrand of Equation 2.4. The integral in Equation 2.4 can be evaluated with
Riemann’s mid-point quadrature or Gauss quadrature, where a faster convergence is

achieved with Gauss quadrature.

2.2.3 Fast nearfield method for a circular piston

An integral equation that combines Equation 2.3 and Equation 2.4 is derived by

McGough et al [23], where

 

 

7f
- a rcosw—a
PU, zit) = pcquJWt_/ 2 2
7r 0 r +a —2arcosz/)
(e—jkx/r?+ a2 — 2arcos¢ + :52 _ e-jkz) dzb. (25)

Equation 2.5 eliminates the numerical singularities that otherwise occur where r z

a. Unlike Equation 2.3, which uses different expressions in different regions, the

10

expression in Equation 2.5 applies everywhere. Equation 2.5 computes the nearfield
pressure more rapidly and more accurately than the spatial impulse response method
[23]. The single integral in Equation 2.5 converges faster by using Riemann’s mid—

point quadrature.

2.2.4 Farfield approximation method for a circular piston

When the distance from the acoustic source to the field point is much longer than
the dimension of the source, i.e., R = VT? + 22 >> a, the pressure can be calculated

by the farfield Bessel approximation [18],

t _ kR) Ira2 2J1(ka sin 0) (2.6)
2R ka. sinl9

 

 

190‘. z. t) = J'pcuoej(°"

where 0 = arctan (r/z) and J1 is the first order Bessel function of the first kind. Since
no integration is involved in Equation 2.6, pressures can be computed rapidly using
this formula. However, in order to achieve a reasonable accuracy with the farfield
approximation method, the field point need to be at least several farfield transition

distances (az/A) away from the source.

2.3 Rectangular piston

2.3.1 Rectangular radiator method for a rectangular piston

The rectangular radiator method divides a rectangular piston into sub—rectangles of
width Aw and height Ah, and the pressure field is then calculated by superposing

the contributions from those sub-elements [31] according to

 

 

 

 

e_ij’ _ I _ I
AwAhZe SianAw($ mnlsianA-My yn),

t = - jwt
p($ayazv ) 9901106 R, 2R1 2R!

(2.7)

 

where R’— — (22+ +(m — m’ 2+ (y - ya)? and m’,, and y], are the coordinates of the

center point of each sub—element. Ref [31] also provides a time-saving scheme by

11

increasing Aw and Ah and reducing the number of sub-elements as the field point

moves further away from the piston. This approximation required that Aw and Ah

Aw and Ah 5 V fig: (2.8)

where F is a control constant with typical values between 10 and 20. At points suffi-

satisfy the criteria

ciently far from the aperture in the farfield, no subdivision is needed. The summation
disappears in Equation 2.7, and the pressure is proportional to the product of two
sinc functions. The resulting expression is sinc(kam/ R)sinc(kby/ R), where a and b

are the half width and the half height of the rectangular piston.

2.3.2 Spatial impulse response method for a rectangular piston

The spatial impulse response method for calculating the pressure generated by a

rectangular piston is evaluated with [21]

2 2
P($,y,z,t) =J'Pcuoejwt Z Z fishy-(513.31%), (2.9)
i=1j=l

where

_ 1:1[23 (e—jkv 2:2 + s2 +l2 _ e—jkz)

3’1 - 27r 2k
\/22 + 32 +12 8 .
— cos‘l (—————) e-Jkfidfi
/ 22 + s2 V :32 - z

, Fz2"+'—32 —+ 12
[V 2:2 + 12

The values of s and l are related to the coordinates of the observation point and the

cos-1 (7%) e—jkfldfi]. (2.10)

size of the rectangular transducer so that 31 = [m — aI, $2 = [m + a], ll = [y — bl, and
l2 = [y + b]. In the plane that contains the transducer surface, four sub-rectangles
are formed by the lines passing through the four sides of the rectangular piston and

the horizontal and vertical lines across the projection of the observation point. The

12

Table 2.1. Combinations of Ii, j based on the location of the observation point (m, y, z).
IIEI Sa, lylsb I2,1+71,2+I1,1+172,2
m>a.|yl_<_5 I2L1+12,2-~’1,1-11,2
x<—a.lylsb 11,1+11,2-12,1-Igg

 

 

 

 

 

 

 

lmlSa,y>b I1,2+T2,2-11,1-12,1
[93] Sa, y< -b I1,1+I'2,1 -11,2-12,2
m>a, >b

x<_§iy<_b I2,2+11,1---11,2-12,1
m< —a,y>b

 

 

 

m>a,y< —b [1.2+12,1 -12.2-11,1

 

physical interpretation of the variables 3 and l are the width and the height of the four
sub—rectangles in the m and y directions. I 3,, l]. is the Fourier transform of the impulse
response from any of these sub-rectangles. Depending on the contribution from each
sub-rectangle at the observation point, a + or a — Sign is chosen for the summation
in Equation 2.9. According to the location of the observation point relative to the
rectangular piston, there are seven combinations of I 31, 1]., or Ii, j for simplicity, which
are listed in the table below. When using Equation 2.10, slow numerical convergence
is expected in the regions where [m] = a and [y] = I) because singularities occur in

these regions.

2.3.3 Fast nearfield method for a rectangular piston

The fast nearfield method for a rectangular piston [24] is formulated as

 

_' 2 2 2 .
12 e ]k\/Z +0 +81—eT-7kz

 

 

 

 

 

' 1
m, ,z,t = ' cu eJWt—(s/ (10
p( y ) Jp 0 27f 1 —11 02“?
32 e—jk‘/22 + 0'2 + l%_e_jkz
+l1/ 2 2 d0
)2 e—jk\/22+02+s§_e_jkz
+82] 2 2 d0
—ll 0 +32
82 €_jk‘/22+02+lg—€_jkz
+12] 2 2 do), (2.11)
-81 0 H2

13

where 31 = a — m, [1 = b — y, .92 = a + m, and Z2 = b + y. In the plane that contains
the transducer surface, four sub-triangles are formed by the lines passing through
the four sides of the rectangular piston and the lines connecting the projection of
the observation point to the four corners of the rectangular piston. The absolute
values of s and l are equal to the base lengths and the heights of these triangles.
By using sub-triangles instead of sub-rectangles, overlapping upper and lower limits
are combined, and the total number of integrals [is reduced from 8 to 4 in Equation
2.11. Both positive and negative values of s and I define the upper and lower limits
of the integrals so that Equation 2.11 holds everywhere in space. The singularities
in Equation 2.10 at [m] = a, [y] = b and z = 0 are eliminated in Equation 2.11 by

subtracting 6.37“ from the integrand of each expression.

2.4 Spherically focused piston

2.4.1 Point source superposition method for a spherically focused piston

The integration of the Rayleigh-Sommerfeld integral over a spherically focused piston

is expressed as

2- , 7r 7r ——jkd
17(7“, 2, t) = __Jp:uoe]wt/() / 0 e
- 0

 

R0 sin Odddw, (2.12)

where

 

d = \/R(2) sin2 0 — 2R0r sin0cos¢ + r2 + [R0 c030 — (z — 120)]2 (2.13)

and R0 is the radius of curvature of the spherically focused piston, 60 is the half
opening angle between the z axis and a line connecting the geometrical focus to the
edge of the piston, d is the distance from a source point on the piston to a field
point, (r, z) is the cylindrical coordinate of the observation point in the plane of zero
azimuthal angle, (120,9,1/2) describes the spherical coordinates of a source point on

the piston surface, 6 and w are the polar angle and the azimuthal angle, respectively.

14

The center of the piston is located at the origin at m = y = z = 0, and the z axis
passes through the geometrical focus.
While Equation 2.12 is straightforward to evaluate, the Rayleigh-Sommerfeld inte-

gral for a spherically focused piston can also be implemented analogously to Equation

 

 

 

 

2.2, where
2' '
mm) = image)...
_' 2 _ 2 2_ _ 2_ 2 2
a "e ]k((/r +(z Z1) +0 2racos¢+R0 \/R0 a +0 )c0320
adadw,
0 0 \/r2+(z—zl)2+a2—2racosr,bcos€0
(2.14)

where a represents the radius of aperture of the spherically focused piston, (0,211) is
the radial coordinate and the azimuthal angle of a source point on the piston, and 2:1
is the distance from the bottom of the piston to the plane that contains the piston
aperture. The geometrical description of the parameters for a spherically focused
piston are described in Figure 2.1. Equation 2.14 is derived by projecting the normal
particle velocity no on the piston surface to the aperture that contains the piston

aperture. The projected particle velocity 11.6 in the aperture plane is

2 0 ,
ab = EClS—uoe’flle, (2.15)

 

where c036 = MR?) —a2/\/R(2)—a2+02, c0390 = R3—a2/RO, and R1 = R0 —

\/R(2) — a2 + 02 is the distance from arbitrary point on the shell surface to a point on

 

the aperture plane along a line that passes through the geometrical focus. Equation

2.14 is obtained by plugging Equation 2.15 into Equation 2.2.

2.4.2 Spatial impulse response method for a spherically focused piston

The impulse response method for a spherically focused piston is derived by Arditi

[22]. The computational region is divided into two regions divided by a double-ended

15

 

 

 

Figure 2.1. The geometry of a spherically focused piston evaluated in a cross section
passing the z axis.

cone formed by lines passing through the geometrical focus and the edge of the piston.

The pressure outside of the double—ened cone is calculated by

 

 

 

 

 

' t r2 230 — 'kfi
p2(m, y, z, t) = jpcuoejw / —_—cos’lfle J (1,8, (2.16)
7., Ar
where
EcosB+? R3 — a2
Q = , (2.17)
r(/Rg — cos2 B
R2 + 22 — a?
B = 0 .
COS 2F 3 (2 18)
and
E = Z — R0,
F = \/m2 + y2 + '22, (2.19)
r = \/$2 + 112.
r1=\/(a—r)2+(3+ Rg—a2)2,
(2.20)

 

r2= \/(a+r)2+(3+ R3-a2)2.

16

The pressure inside the double-ended cone is the sum of p2(m,y, z,t) and an extra

term,
t t . jwthO e—jkrl — e—jkro, E < 0, 2 21
P1($,y, 21 )_ p2($aya Z) )‘l’JPC’UOe ? e_jk7~0 _ e_jkr2, .2. > 0, ( ' )
where
R0 — F, E < 0,
= 2.22
To { R0 +13 2 > 0. ( )

A singularity occurs at the geometrical focus of the spherically focused piston (E = 0),

and the pressure at the geometrical focus is given by

me, y. z. t) = jpcuoejwtkmo — R3 — a2)e—ijo (2.23)

2.4.3 Fast nearfield method for a spherically focused piston

The fast nearfield method for a spherically focused piston is expressed by a single

 

integral,
' a
19(23, :1. z, t) = J'pcuoejwt—RO-
7r
7f rcostl) R3—a2+2a
/ Xdi/J, (2.24)
0 1287‘2 + 2am cosw R3 - a2? —- a2r2 cos2 w + 32:12

where

 

_- 2 —2 _ — 2 _ 2 . . _ _
X = e ]k\/RO + r 2am COSt/J + 2z(/RO a _ e—]k[R0+SIgn(z)r], (2.25)
and the function sign is defined as

. 1, z > 0,
Slgn(z) = { _1 z < 0 (2.26)

The singularities are eliminated, and Equation 2.24 holds everywhere in space except

at the geometrical focus where the pressure is defined in by Equation 2.23.

17

CHAPTER 3

Angular Spectrum Approach

3. 1 Introduction

The angular spectrum approach describes the diffraction of acoustic waves from finite
apertures by superposing plane waves traveling in different directions [30] and prop—
agating these components in the spatial frequency domain. As opposed to integral
approaches that calculate the field at each observation point, the angular spectrum
approach computes the pressure field in successive planes with a 2D FFT, which
speeds up these calculations significantly. The angular spectrum approach uses ei-
ther the normal particle velocity or the pressure as the source, and then the spectral
propagator or the 2D Fourier transform of the spatial propagator is multiplied by the
2D Fourier transform of the source to simulate the propagation of acoustic waves. The
spectral propagator is described as an analytical function [35, 36, 37, 38, 39, 40, 41],
whereas the spatial propagator is calculated analytically and then transformed into
the spatial frequency domain with a 2D FFT [42, 43, 44].

The numerical accuracy of the angular spectrum approach has been explored with
great detail. Williams and Maynard [45, 46] discuss the difference between the an-
alytical Fourier transform and the discrete Fourier transform in angular spectrum
simulations. Williams and Maynard also analyze the effect of the windowing func-

tion, the aliasing errors due to the lack of spectral samples, and the singularity of

18

the spectral propagator. Orofino and Pedersen [40, 41] derive a relationship between
the angular range for the decomposed plane wave components and the sampling rate
of the spectra in the spatial frequency domain. These parameters are correlated to
the spatial sampling rate and the combination of these determine the accuracy of the
angular spectrum simulation. Christopher and Parker [42] simplify the angular spec-
trum calculation for an axisymmetric radiator by using the discrete Hankel transform
instead of the discrete Fourier transform. A spatial frequency truncation technique is
proposed in Ref. [42] to reduce the error caused by the undersampling of the spectral
propagator. Wu and Stepinski [37, 38, 39] investigate the spatial frequency truncation
technique and derive the optimal truncation frequency for accurate angular spectrum
simulations.

This chapter compares the performance of the spectral propagator and the spatial
propagator in terms of the numerical accuracy and the computation time. Section
3.2 presents the theory of the angular spectrum approach. In the remaining sections
of chapter 3, new simulation results are presented. Numerical implementations of
the spatial and spectral propagators are described in sections 3.3.1, 3.3.2, and 3.3.3.
Spatial truncation is found to be the main source of error for the spatial propagator.
The aliasing error caused by the undersampling of the spectra induces most of the
errors when calculating with the spectral propagator. Section 3.3.4 compares the
spatial and spectral propagators in both spatial and spatial frequency domains. The
performance of the two propagators is evaluated for a square piston in a homogeneous
and linear medium. Modeling attenuation with the angular spectrum approach is
introduced in section 3.4, and the frequency filtering effect of the attenuation term
is explored. The influence of the source distribution on the accuracy of the angular
spectrum simulation is discussed in section 3.5. Backward propagation of the pressure
field is described in section 3.6. The loss of high spatial frequency components during

backward propagation is found to be the limiting factor for predicting the initial

19

wavefront in these simulation.

3.2 Theory of the angular spectrum approach

The complex field excited by a monochromatic source can be described by the super-
position of an infinite number of plane waves [30]. Each of these plane waves travels
in one direction corresponding to an unique wavevector (km, kg, 16;) that obeys the

relationship

1.3+ 1.3+ k3 = k2, (3.1)

where k = 27r//\ is the spatial wavenumber. If some of the plane wave components
travel in the z direction, then the complex field in a plane orthogonal to the z direction
can be decomposed in the spatial frequency domain (k space) with the 2D Fourier

transform,
P(km,kyrzl= / / Maude—“km“kyy)dxdy. (32)
—oo —oo

where p(m, y, 2) represents the complex acoustic pressure in a transverse plane with
constant 2, and P(kx, Icy, z) is the corresponding angular spectrum. Likewise, the 2D

inverse Fourier transform recovers the pressure from the angular spectrum with

1 oo oo ,

In linear acoustics, the Helmholtz equation is satisfied by a monochromatic field

in the spatial frequency domain [36],

62P(kx, Icy, 2)
(3‘22

 

+ (k2 — kg. — k§)P(kx, ky, z) = 0. (3.4)
The solution of Equation 3.4 satisfies
P(k$,ky,2) = P0(k$,ky,ZO)Hp(kx,ky,AZ), (3.5)

and

 

_' 2_ 2_ 2
Hp(k$,ky,AZ) =8 JAZ\/k k3 kg,

20

where P0(k$, kg, 20) is the angular spectrum of the source in the transverse plane at
20, P(kx, ky, z) is the angular spectrum of the pressure field in a parallel plane at z,
and H (kx, ky, A2) is the spectral propagator function with A2 = z — 20.

The angular spectrum of the complex pressure can also be obtained from the

normal particle velocity distribution on the radiator surface,
P(kg;, ky, Z) = ijkU0(k$, ky, 0)Hu(k$, ky, A2), (3.7)

where U0(kx, kg, 0) is the angular spectrum of the normal particle velocity on the
radiator surface at z = 0, Hu(km, ky, A2) is the corresponding spectral propagator
function, and p and c are the density and the speed of sound, respectively. Hu can

be derived from the linear Euler’s equation, which is

-J'pckU(:v. y, 2) = Vp(rv. y, 2)- (3-8)

The linear Euler’s equation is transformed into the spatial frequency domain via

spatial Fourier transform,
kaU<kx,ky,Z) = kzP(k$,ky,Z). (3.9)

By dividing Equation 3.5 by Equation 3.7, and substituting P0 for U0 according to

Equation 3.9, the spectral propagator for a velocity source is obtained,

 

1 e—jAz\/k2 — kg —— k5,
j\/k2 — k3, — k5

assuming that the acoustic radiator is located in the plane of z = 0 and mounted in

Hu(k$, kg, AZ) =

 

(3.10)

 

an infinite rigid baffle.

The spatial equivalents of Equation 3.5 and Equation 3.7 are [46]

p($1y1 Z) =p0(x7ya ZO) ®hp(xaya AZ): (3'11)

19(56, :1. z) = jpckuohz. y. 0) ® hams/.132). (3-12)

21

where the symbol 8) denotes spatial convolution, 110(1), y,0) and p0(m,y, .20) are the
normal particle velocity and pressure distributions in the source plane, respectively,
and hu(m,y,A2) and hp(m,y,Az) are the spatial propagators associated with the
particle velocity and pressure sources, respectively. By applying a 2D inverse Fourier
transform to Hu(k$, ky, A2) in Equation 3.10, hu(m, y, A2) is obtained as [43]
e—jkr

27rr‘ ’

 

hu(m, y, Az) = (3.13)

 

where r = \/m2 + y2 + A22. This is the Green’s function for the Rayleigh-Sommerfeld
integral. Similarly, the spatial propagator hp(m,y, A2) for a pressure source is ob-

tained as

hp(m, y, A2) = £0 +jkr)e_jkr. (3.14)
2nr3

3.3 Influence of the propagators

The spatial and spectral propagators each have distinct advantages and disadvan-
tages, and these features determine the trade-offs between the numerical accuracy
and the computation time for angular spectrum calculations. The spatial and spec-
tral propagators for a velocity source, hu(m, y, Az) and Hu(kx, ky, Az), are compared
in section 3.3. The behaviors of hp(m,y, A2) and Hp(k$, ky, Az) are similar to the
corresponding terms for the particle velocity source, and therefore results for these

are not shown.

3.3.1 Spatial propagator

In computer implementations of the angular spectrum approach, the infinitely extend-
ing spatial propagator hu(m, y, A2) is truncated in both the m and the y directions so
that the extent of each computational plane is D X D. A 20. X 2a square piston with
2a < D is illustrated in Figure 3.1. With a sample spacing of 6, the D X D plane

is discretized to N X N grid points, where N = D/d + 1. The spatial propagator is

22

 

 

A m
V V

Figure 3.1. The D x D source plane consisting of a non-zero normal particle velocity
distribution in a 2a X 2a square area on the piston surface. The remaining area is filled with
zeros.

evaluated in this discretized grid, and the discrete coordinates are

m=m6, m=—-N/2+1+¢,...,N/2+¢>,

.I
y=n6, n=——N/2+1+¢,---,N/2+¢a (3 5)
where
—1/2, when N is odd,
qb _ { 0, when N is even. (316)

With this value of q), the grid is symmetric for odd N. For even N, the grid is biased
so that points on the central axis m = y = 0 are included in the three dimensional
grid. A two dimensional FFT is then applied to both hu(m, y, Az) and u(m, y, 0), and
the results are multiplied in the spatial frequency domain. The complex pressure is
the inverse two dimensional FFT of the product.

Figure 3.2 shows the normalized reference pressure generated by a 3cm X 3cm
square piston (a = b = 1.50m) in the mz plane evaluated at y = 0. The piston
is excited at a frequency of 1MHz, and the speed of sound is 1500m/s; therefore

the wavelength is A = 0.150m. The reference pressure is computed with the fast

23

nearfield method and normalized to the maximum pressure amplitude. The pressure

  
    
 

93
3
8’. 1
g 0.8

0.6
U 30
o 0.4 I,
E‘ 0.2 25
g 0 20
54.5 15
z

_3 z(Cm)

Figure 3.2. Two dimensional cross section of the three dimensional reference pressure
generated by a 3cm X 3cm square piston in non-attenuating media. The excitation frequency
is 1MHz, and the normal particle velocity distribution is uniform across the piston surface.
The reference pressure is computed with the fast nearfield method. The result is normalized
to the maximum pressure amplitude computed in the three dimensional volume.

is also simulated with the angular spectrum approach using the spatial propagator
in a three dimensional volume, and a cross section of the pressure field in the mz
plane at y = 0 is shown in Figure 3.3. The extent of the three dimensional volume is
defined by D = 9cm in both the m and the y directions, and the three dimensional
volume extends 30cm in the z direction. With a sample spacing of 6 = 0.03cm
(A/5), the transverse plane is discretized to a 301 X 301 grid (N = 301). In Figure
3.3(a), a 301 X 301 point FFT is used to transform the spatial propagator into the
spatial frequency domain, while in Figure 3.3(b), the extent of the spatial propagator
is enlarged with zero padding. A 512 X 512 point FFT is then computed, and the
resulting field is truncated to a 301 X 301 grid. The pressure fields in both Figure
3.3(a) and Figure 3.3(b) closely resemble the reference field in the central portion of
the computational grid; however, near the edge of the computational grid, a discrete

jump appears in the computed field. These numerical errors are caused by the finite

.24

extent of the D X D window that truncates the spatial propagator. In Figure 3.3(a),
the error is represented by an artifact of circular convolution between the spatial
propagator and the velocity source, where the spatial propagator is evaluated in a
D X D window, and the size of the non-zero velocity source is 2a X 20.. The circular
convolution performed over a D X D area causes the fields to wrap around at the
edge. In Figure 3.3(b), where adequate zero padding is used to prevent the circular
convolution error, the truncation of the spatial propagator causes errors in the same
region. With zero padding, as soon as the non-zero part of the velocity source is
shifted to a location that overlaps with a value outside of the D X D rectangular
window, the source is convolved with zero instead of the correct value of the spatial
propagator. As a result, the pressure fields computed with the spatial propagator
is erroneous starting from a distance a from the edge. With zero padding, the edge
errors in Figure 3.3(b) are lower than those in Figure 3.3(a); however, errors are
always associated with the spatial propagator within a peripheral band of width 0..
The L X L area in the center consistently avoids these problems, where L = D — 2a.
The spatial propagator hu(m, y, A2) is smooth everywhere, although a singularity
occurs on the piston surface at (0,0,0). Close to the radiator surface, hu(m,y, A2)
has a steep slope, and a high spatial sampling rate is required. Away from the piston
surface, the Nyquist sampling rate of the spatial propagator is easily satisfied when

6 S A / 2 is the sampling interval [40].

3.3.2 Spectral propagator

The properties of the spectral propagator Hu(kg;, Icy, A2) in Equation 3.10 change
depending on the coordinates in the spatial frequency domain. In the region where
kg + kg 3 k2, Hu(kx, kg, A2) propagates the field in the 2 direction by applying a
complex weight to each plane wave component. In the region where k3: + [(232] > k2,

Hu(kx, Icy, A2) exponentially attenuates evanescent waves. Therefore, Equation 3.10

25

099994
mem

Normalized pressure
4:.
'01

 

(3.) Simulated pressure obtained from the spatial propagator
with D = 9cm and without zero padding (N = 301).

\ \ \ .
,(x\ (.1 (n\

\\ \\\ N \\
\\\\\\(\\ \

t \ \
\“t\‘\\“\ \ \\\‘\\‘\ ‘
,. M“ (o ,\\\ \(\\\(\\\\\\\\\\‘(
\‘\ M \\ \ \‘\‘\‘\
\‘\\‘(\\\\\\\\\\\\(\‘\\\\ \\ \\ 1‘. .
\
\\\\\\\\\(\\\(\\‘,\\ \
M
,\

    

\\

\ > ‘ "if
. Ni \ \‘\ \ ‘
all " ,0 ill, ‘ “ 30

l] 1
“ill-"1"“! 1'4 1. 25

    

09999;
mama

P
0':

Normalized pressure

z(cm)

(b) Simulated pressure obtained from the spatial propagator
with D = 9cm and with zero-padding to N = 512.

Figure 3.3. Simulated pressure generated by a 3cm X 3cm square piston in non-attenuating
media computed by the angular spectrum approach, where (a) the spatial propagator is
truncated with a 9cm X 9cm window and evaluated in a 301 x 301 grid for FFT calculations,
and (b) the spatial propagator is truncated with a 9cm X 9cm window and zero-padded to
a 512 X 512 grid for FFT calculations. The sample spacing is 6 = 0.075cm in both (a) and

(b)-

26

is equivalently stated,

 

e—jA2\/k2 — kg — kg,

 

k§+k§sk2

j k2 _ k2 _ k2
—Az,/k,, + ky — k
6 k3, + 1:5 > k2

flag—1.2

When the spectral propagator is discretized to an N X N grid, the discretized trans-

 

 

verse wavenumbers are

k, = mAk, m = —N/2 + 1+ (1), N/2 + a,

.1
k, = nAk, n = —N/2 + 1 + as, ...,N/2 + 4;, (3 8)

where qfi is defined in Equation 3.16. The maximum value of kg; or Icy is 7r / 6 , and the
sample spacing in the spatial frequency domain is Ak = 21r / (N 6).
If the effects of discretization are included, the inverse Fourier transform of the

angular spectrum in Equation 3.7 is

pD(xa y, 2) = ijka-1{U0(k$a ky, O)H’U(k$: ky, AZ)C0mb(k$/Aka [Cy/AM}
= J'pckp($.y12) ® comb(55/ (N 5). y/ (N 5)),

(3.19)
where comb() represents the comb function
comb m(/a, y/b)= 219:2“ m —pa)6 (—y qb). (3.20)
The Fourier transform of the comb function is
a {comb(z/a y/b» = Z 2601: — P2—7:)5(y - q-ZI). (3.21)
’ p q a 5

According to Equation 3.19, discretization in the spatial frequency domain leads to
circular convolution in the spatial domain. The corresponding spatial field p(m, y, z)
is replicated and shifted to multiple locations centered at (pN 6, qN 6), where p and
q are arbitrary integers, and the results are added. Since p(m, y, 2) extends to +/—
infinity in the spatial domain, the tails from the replicas of the field leak into the

computational area and cause aliasing errors. These errors are illustrated in the

27

pressure meshes shown in Figure 3.4. Figure 3.4(a) shows a cross section of the
pressure field simulated with the spectral propagator using N = 301. Figure 3.4(b)
shows the pressure simulated with the spectral propagator using N = 512. The
pressure field in Figure 3.4(a) contains significant errors due to the aliasing of the
spectral propagator. These errors are somewhat reduced when a larger N is used in
Figure 3.4(b). Larger values of N correspond to greater extent of the spatial field
N6, which moves the spatial replica farther away from one another. Therefore, the
error induced by aliasing is reduced as N increases.

The spectral propagator Hu(k$, ky, A2) also has singularities at kg + kg, = k2,
which causes some difficulties in the numerical evaluation of the spectral propagator.
In an effort to address the problem with these singularities, different values were
substituted for the infinite values encountered in these locations, including 0, the
amplitude of an adjacent point, and an average value of the adjacent points around
the singularity [45]. Each of these strategies produces very similar errors. This
suggests that the numerical singularities at k3c + kg = k2 have little influence on the

overall error.

3.3.3 Spectral propagator with angular restriction

The partial derivatives of the spectral propagator with respect to kg; and kg of the
phase represent how rapidly the complex spectrum oscillates. According to Equation

3.10, the phase and the derivatives of the phase are

 

 

(0(kz, kg) = —jA2\/k2 — k2,. — kg, (3.22)
68—20— : z'Azkx
9’ k2 — k3, — 23’
y \/k2 — kg — 1.3

Since flap/8kg; and ago/[My approach infinity when kg; +k§ approaches k2, the real and

imaginary parts of the spectral propagator oscillate more rapidly near kg + k5 = k2.

28

   

  

l‘I’phhj‘i'

 
     
 
  
  

  
 
 
 
 

   
 

   
  

    
     

a.)

3

U)

(D

2 018 1’1,

Q ' 1' 1'11, r"!

a 0.6 ‘ 1.114111 "" ‘

ID 0.4 110111" 11 Id" ""

N ””4““ “3103.191“ “‘1‘“ A.

:= 0.2 at? 11".” lil'l'“ 'I'

g 0 “I,“ 'l""\l:1‘lrr' 11"” ‘ 20

‘54-5 3 934 15

z 10
z(cm)

(a) Simulated pressure obtained from the spectral propaga-
tor with N=301.

   

 
  
 
 

.’
III"
‘;5. . I":

   
 
  

         
 

H
(“1'7
‘ ",r/ .“v’ ’0"
"I-’ 10'1" 5,2,1, "“3 [5‘1” 1",?”
[1‘11 15,): 0H
”/31”;

  

\|'/,')'n'

  
    
  
 
 
 

.1, MM" I
”If," ”1’”
/l ‘U
', ’l‘flvl'

'.5."‘I

 
 

09999;
Namm

:‘i‘
0':

Normalized pressure

(b) Simulated pressure obtained from‘the spectral propaga—
tor with N=512.

Figure 3.4. Simulated pressure generated by a 3cm X 3cm square piston in non-attenuating
media computed by the angular spectrum approach, where (a) the spectral propagator is
evaluated in a 301 X 301 grid for FFT calculations, and (b) the spectral propagator is
evaluated in a 512 X 512 grid for FFT calculations. The spatial frequency extent for the
spectral propagator is 2km = 27r/ 6, and the sample spacing is 6 = 0.075cm in both (a) and
(b).

29

At kg + k; = k2, 61,0/ 8km and acp/aky encounter a singularity, and Hu(kx, ky, A2) is
inherently undersampled. This undersampling in the spatial frequency domain leads
to aliasing errors in the spatial domain [39]. The aliasing error is especially severe
for large 2 because the rate of change 690/ka or 61p/ 819;, is proportional to A2. The
undersampling can be reduced by increasing N at the expense of increased computa-
tional cost. Christopher [42] and Wu [38, 39, 37] use a spatial frequency truncation
technique as an alternative solution to this problem. This technique reduces the alias-
ing errors without increasing the size of the computational grid. Spatial frequency
truncation is achieved through angular restriction, which applies a radially symmetric
low-pass filter to the spectral propagator function Hu. The cut-off spatial frequency,

or angular threshold, is given by [37]

D2/2
k = 19‘1— 3.24
C D2/2+z2 ( )

where D = N 6. Equation 3.24 is obtained by relating the maximum spatial extent
D to the Nyquist frequency in the spatial frequency domain. Angular restriction
removes the under-sampled angular spectra and prevents the high spatial frequency
components from leaking into the propagating field.

Figure 3.5 shows the pressure simulated using the spectral propagator with angular
restriction. The results in Figure 3.5(a) and Figure 3.5(b) are computed with N = 301
and N = 512, respectively. The aliasing errors in Figures 3.5(a) and 3.5(b) are
greatly reduced relative to the corresponding Figures 3.4(a) and 3.4(b), where no
angular restriction is used. The results demonstrate that angular restriction effectively
reduces the aliasing errors caused by under-sampling the highly oscillatory spectral

propagator.

30

Normalized pressure
N «D- 03 co

    

x(cm)' —3 -4.5 0 z(cm)

(a) Simulated pressure obtained from the spectral propaga-
tor with N=301 with angular restriction.

 

    
 

       
     
  

          

g

(I, 00’

u) 4:,1'910- .

a, 1 1129516723,
58.g 60!,”

'0 ' Il’ ‘ 30
a”, 8'3 W“ 25
= . //’//

g 0 WI! 20
c,4.5

Z

z(cm)

(b) Simulated pressure obtained from the spectral propaga-
tor with N=512 with angular restriction.

Figure 3.5. Simulated pressure generated by a 3cm x 3cm square piston in non-attenuating
media computed by the angular spectrum approach with angular restriction, where (a) the
spectral propagator is evaluated in a 301 x 301 grid for FFT calculations, and (b) the spectral
propagator is evaluated in a 512 x 512 grid for FFT calculations. The spatial frequency
extent for the spectral propagator is 2km = 27r/6, and the sample spacing is 6 = 0.075cm
in both (a) and (b).

31

3.3.4 Comparison between the spatial and spectral propagators

Although the spatial and spectral propagators are analytically equivalent to each
other, the performance of the two propagators differs. Figure 3.6(a) and Figure
3.6(b) contain the axial and transverse pressure data extracted from Figures 3.2,
3.3(a), 3.4(a) and 3.5(a). The reference pressure is provided by the fast nearfield
method, and the remaining results are obtained from the angular spectrum approach
with the spatial propagator using D = 9cm and the spectral propagator with and
without angular restriction using N = 301. The axial pressure in Figure 3.6(a) shows
that the pressure computed with the spectral propagator without angular restriction
(dotted line) is corrupted by high spatial frequency errors. The field computed with
the spectral propagator using angular restriction (dash-dot line) contains some ripples
corresponding to intermediate frequencies. The field computed with the spatial prop—
agator (dashed line) closely tracks the reference (solid line). Figure 3.6(b) emphasizes
the circular convolution errors in the result from the spatial propagator (dashed line)
within an edge band of width a (a = 1.5cm) along each side. The errors from the spec—
tral propagator without angular restriction (dotted line) and with angular restriction
(dash-dot line) are not restricted to this edge band. Figure 3.6 demonstrates that, in
non—attenuating media, the spatial propagator computes pressures more accurately
than the spectral propagator in the central portion of the grid, while relatively large
errors are induced by the spatial propagator starting at a distance a from the edge of
the computational planes.

The FFTS of the spatial propagator and the highly oscillatory spectral propagator
are compared in Figure 3.7. In Figure 3.7, the spectral propagator Hu(k$,ky,z)
(dashed line) is computed with N = 301 at z = 15cm and shown along the kg;
direction at Icy = 0. The real and imaginary parts of Hu(kx, ky, A2) vary rapidly
at large kg; and Icy, and the amplitude of Hu(k;,;, Icy, Az) increases quickly with kg:
and ky until the singularities are encountered at ”kg + k5 = k. The rectangular

32

 

L

 

—Reference
------- Unrestricted spectral propagator
1.5. '---' Restricted spectral propagator i
- - - Spatial progator

 

 

 

:5 i

5. .
3. -‘ :g -‘-~:

Normalized pressure

 

 

 

 

0 1 0 20 30
z(cm)
(a) Simulated aixal pressure along 2 direction at a: = y = 0.

2

 

 

—Reference

------ Unrestricted spectral propagator
----- Restricted spectral propagator

- - - Spatial propagator

.3
U1

 

 

.0
or

Normalized pressure
A

 

 

 

94.5 -3 -1.5 0 1:5 31 4.5

 

(b) Simulated radial pressure along a: direction at y = 0 and
z = 15cm.

Figure 3.6. Axial and transverse plots of the simulated pressure generated by a 3cm x 3cm
square piston in non-attenuating medium. Results are shown for the reference computed
by the fast nearfield method, the angular spectrum approach using the spectral propagator
without and with angular restriction, and the spatial propagator. Results are computed
using N = 301 and D = 9cm.

33

window superposed on the spectral propagator plot illustrates the radially symmetric
low-pass filter for angular restriction. The width of the window kc = 0.39219 is
calculated with Equation 3.24. The spatial propagator hu(:r,y, A2) is evaluated in
a 9cm x 9cm window in the my plane at z = 15cm with a sample spacing of 6 =
/\/ 5 = 0.03cm, and a two dimensional FFT is then applied to hu(:r, y, A2). The
FFT of the spatial propagator overlaps with the spectral propagator in the central
portion corresponding to low spatial frequencies. However, a significant difference
appears in the higher spatial frequency region. The spectral propagator (dashed
line) has a full spatial bandwidth within —k and 19; however, the spectrum outside
of the truncation window is under-sampled, and the aliasing errors in Figure 3.6 is
induced. If the sampling rate in the spatial frequency domain is increased, the aliasing
error is reduced correspondingly. The spatial propagator is evaluated analytically and
truncated by a 9cm x 90m window. In the spatial frequency domain, this is equivalent
to convolving the spectral propagator with a two dimensional sinc function. The
high spatial frequency components are attenuated, however, the more important low
spatial frequency components are preserved.

The real and imaginary parts of the spatial propagator are plotted as solid lines
in Figure 3.8. The spatial propagator h(:2:, y, A2) is evaluated analytically from a: =
—4.5cm to a: = 4.5cm where y = 0 and z = 150m. The dashed line in Figure 3.8 is the
inverse FFT of the spectral propagator with angular restriction. A two dimensional
inverse FFT with N = 301 converts the spectral propagator to the spatial domain,
and the result is shown in a line along the :1: direction at y = 0 and z = 15cm. Even
with angular restriction, some under-sampled components with intermediate spatial
frequencies still appear in the spectral propagator. The ripples in the dashed lines
are caused by these aliased components.

To quantitively evaluate the errors involved with the spatial and spectral propa-

gators, a root mean squared error (RMSE) is introduced and the overall error in a

34

 

 

- - - Hu(kx,ky,z)

P a)
or

 

 

 

    

 

 

 

S _FFT of hu(x,y,z)
‘8 3 - . ~
3 l. :
“a 15 l
t l
N .I
2' 0 M d
an a r-
m a r-
m -1 5+ " ll
:. ii
I 'r
_3 1 1L

-3 -é —1 d
kx(normalized by k)

(a) The real part of the spectral propagator and the FFT of
the spatial propagator along the kx direction for kg = O.

 

 

 

 

 

 

 

3 . : ,
E ii
*3 15 l: . :l'iil
w .,.
a o .. . ,
“5 '. l E:
5 ‘1 5 . 1:3
Q l
.5? 4.5? - - - Hu(kx,ky,z) _
g _FFT of hu(x,y,z)

 

 

 

 

—3 -é -i d 1 2 3
kx(normalized by k)

(b) The imaginary part of the spectral propagator and the
FFT of the spatial propagator along the km direction for
Icy = 0.

Figure 3.7. The real and imaginary parts of the spectral propagator (without angular
restriction) and the FFT of the spatial propagator in a non-attenuating medium evaluated
at z = 15cm. The spectrum is plotted along the k3; direction for kg = 0. The spatial
prepagator is computed with D = 9cm, and N = 301 is used for FFT calculations. The
radially symmetric window applied to the spectral pr0pagator shows the angular threshold
is at kc = 0.3917k.

35

 

 

 

 

 

 

 

 

'5 3 - , ,

2; —h(x.y.2)

g 2 _ ---lnverse FFT of H(k ,k ,z).

o l- , x y .

E. :l r. Fl '0 I‘| :: I‘ I: I... f:

'U 1 ’ I ' ' I‘ . . ’I ' ' I '

G) H i ' H

g I: s :: o "\

g 0 I gr. : ‘I ' i... ‘ I, |‘ ‘.5 |

8 I ' .' I ' ‘ | " I II

3 -1; 3; z. '. : .: is '1

5,2. 'e -: .

a

‘3‘ 33 5 —3 -1 5 o 1.5 3 4.5
x(cm)

(a) The real part of the spatial propagator and the inverse
FFT of the spectral propagator along the km direction for

 

 

 

 

 

 

 

 

kg = 0.

a 3 - -

g —h(X.Y.Z)

g 2 _ ---Inverse FFT of H(kx,ky.Z)1
a ::

8 1( I I: . ' i
g ' ‘ ' ‘ I 3
cu h. g , ~ ‘ '. 1...;
E 0: :. , " i
8 ': \ i

u— ‘1 .: :l I, . ’ II ' l.

O ., - ‘0 V '.' I
t

a-z-

§

— :25 -3 -1.5 o 1.5 3 4.5

x(cm)
(b) The imaginary part of the spatial propagator and the in-

verse FFT of the spectral propagator along the k1; direction
for by = 0.

Figure 3.8. The real and imaginary parts of the spatial propagator and the inverse FFT of
the spectral propagator (with angular restriction) in a non-attenuating medium evaluated
at z = 15cm. The pressure field is plotted along the a: direction for y = O and normalized to
the maximum pressure amplitude along this line. The spatial propagator is computed with
D = 9cm, and N = 301 is used for the inverse FFT calculation. The angular threshold for
the spectral propagator is kc = 0.3917k.

36

three dimensional volume is calculated. This error metric defined for this evaluation

is given by

 

 

M E: 2: mak— ’3’ 2 3.25

$.71

where the superscripts (i, j, k) represent discrete field points in the computational
grid, and 72.13, ny, and n2 describe the number of points in each direction. The variable
pmf denotes the complex reference pressure, and the variable p is the complex pressure
computed with the angular spectrum approach.

In Figure 3.9(a), the root mean squared errors are calculated in three dimensional
volumes for different D values, where the parameter D represents the extent of the
computational volume in both the a: and the y directions. The results are plotted as a
function of D, where the value of D ranges from 6cm to 15cm. The root mean squared
error for the spatial propagator (dashed line) is between that of the spectral prop-
agator without angular restriction (dotted line) and with angular restriction (solid
line). Thus, when the circular convolution errors produced by the spatial propagator
are included, the spectral propagator with angular restriction achieves a smaller root
mean squared error.

In Figure 3.9(b), the volumes are truncated at the edges to exclude the windowing
errors produced by the spatial propagator, and then the root mean squared errors
are computed. The truncated volume has a lateral extent of L x L, so the results are
plotted as a function of L with L = D - 2a, where L ranges from 3cm to 12cm. In
Figure 3.9(b), the errors from the spatial propagator are 2 to 4 times smaller than
those in Figure 3.9(a). The errors computed for the spectral propagator with and
without angular restriction are influenced less by the truncation, where the errors
without angular restriction are about 1.5 to 2 times higher than those with angular
restriction in both Figure 3.9(a) and Figure 3.9(b). The difference between the dashed
curves in Figure 3.9(a) and Figure 3.9(b) reinforces the need to truncate the field

computed with the spatial propagator to exclude the errors produced at the edge of

37

the grid.

The spectral propagator function is analytical, and the non-zero portion of
Hu(k$, ky, Az) is confined within kg. + kg S k2 except for some evanescent waves at
small values of Az. However, Hu(k1;, Icy, z) encounters under-sampling problems with
rapidly oscillating real and imaginary components as «kg + 1632/ approaches It. This
problem is alleviated by removing the under-sampled spectra near k = (Hag; + 1:32] = k
via angular restriction. The spatial propagator is analytical, smooth, and slowly-
varying everywhere except at locations adjacent to the radiator. Numerical imple-
mentations require truncation of the infinite spatial propagator in both the a: and the
y directions, and a two dimensional FFT then converts the result into the spatial fre-
quency domain. The finite windowing of the spatial propagator produces errors near
the edge of each transverse plane, whereas the central portion of the computational
grid is unaffected by these errors.

Computations with the spatial propagator take longer than calculations with the
spectral propagator because an extra 2D FFT calculation is required in each plane to
convert the spatial propagator into the spatial frequency domain. Simulations of the
pressure generated by a 3cm x 3cm piston with N ranging from 301 to 1024 show that
the computation time for the spatial propagator is 1.1 to 1.9 times longer than that
for the spectral propagator. This ratio approaches the smaller value as N increases.
For example, after the pressure is calculated with N = 301 (D = 9cm), the field is
truncated to a 201 x 201 x 200 grid (6cm x 60m x 300m volume) for error evaluations.
On a 2.4GHz Pentium 4 PC with 1GB RAM, the spectral propagator computes the
complex pressure in 15.22 seconds with a root mean squared error of 0.044, and
the spatial propagator computes the result in 28.90 seconds with a normalized root
mean squared error of 0.025. When N = 512 (D = 15.330m), the angular spectrum
approach with the spectral propagator computes the complex pressures with a root

mean squared error of 0.027 in 1.23 minutes, whereas the spatial propagator takes

38

 

 

 

 

 

 

 

 

 

0.2 a a . .
------ Unrestricted spectral propagator

3 — Restricted spectral propagator
at: 0.15. - - - Spatial propagator
E .
(U
:3
O'
U)
C
(U
a:
E
‘5'
o
a:

0 l I A I

6 8 10 12 14 16
D(cm)

(a) Root mean squared errors computed in a three dimen-
sional volume with lateral dimensions D x D.

 

 

 

 

 

 

 

 

 

0.2 1 . . .
------- Unrestricted spectral propagator

B —Restricted spectral propagator
at: 0.1533”. ---Spatial propagator
(U
:3
D'
(I)
C
(U
0)
E
‘5
O
a:

0 3 5 7 a 11 13

(b) Root mean squared errors computed in a three dimen-
sional volume with lateral dimensions L x L.

Figure 3.9. Root mean squared errors for the pressure generated by a 3cm x 3cm square
piston excited with uniform normal particle velocity and evaluated in a non-attenuating
medium. The errors are evaluated in three dimensional volumes where the lateral dimensions
are(a)DxDand(b)LwaithL=D—2a. ,

39

1.99 minutes using the same parameters while achieving a root mean squared error of
0.025. As the computational grid becomes larger, the difference between the spatial

and spectral propagator becomes smaller in terms of both error and time.

3.4 Angular spectrum approach in attenuating medium

Simulations in non-attenuating media are informative, as many experiments are con-
ducted in water, which has a very small attenuation coefficient. However, most bio-
logical tissue and fluids have significantly larger attenuation coefficient. In biological
tissue media, the acoustic pressure is attenuated exponentially as 6“"ber — til, where
b indicates the power law dependence of the attenuation coefficient on the frequency.
The attenuation is incorporated into the spatial propagator by replacing the real-
valued wavenumber k in hu(:r, y, z) with the complex wavenumber k — ja f b. For the

spectral propagator Hu(ka;, ky, 2), an exponential term is multiplied by

 

_ b 2_ 2_ 2
S(kx,ky,z)=e “f “M“ k1” kg, (3.26)

which is a simplified expression that is mathematically equivalent to the trigonometric
attenuation term in Equation 12 of Ref [42]. S (km, kg, 2) accounts for the attenuation
of plane waves traveling in the (km, ky) direction. Within the region k3, + 11:32, 3 k2,
S(kx, Icy, z) is analogous to a spatial low-pass filter. Figure 3.10 contains two cross
sections of S (km, kg, 2) along the kg; direction for ky = 0. The solid line shows
the result for z = 0.6cm, and the dash-dot line shows the result for z = 15cm.
As 2 increases, the peak amplitude of 3(km, Icy, 2) decreases, and the higher spatial
frequency components are increasingly attenuated. As shown in Figure 3.10, the
3dB bandwidth for the attenuation term at z = 0.6cm is 1.98117, whereas the 3dB
bandwidth for the attenuation term at z = 15cm is reduced to 1.3919.

Figure 3.11 shows the reference pressure generated by the same 3cm x 3cm square

piston in attenuating media with a = 1dB/cm/MHz. The reference pressure is cal-

40

 

 

 

 

 

 

 

 

 

 

' ' —z=0.6cm
0 8- K N ----- z=15cm
E
2
c 0.6“
.2
E
2 0.4L
8
<
0.21 , ...... ~
0 . i’ '\

-3 -2 -1 d 1 .9. 3
kx(normalized by k)

Figure 3.10. A cross section of the attenuation term S(kx, log, z) along the km direction
for Icy = 0. The results are shown for z = 0.6cm (solid line) and z = 15cm (dash-dot line).

culated with the fast nearfield method [24] using 100 Gauss abscissas, which achieves
11 digits of accuracy compared to the result obtained with 2000 Gauss abscissas.
The pressures from the angular spectrum approach in Figure 3.12 are computed with
D = 90m (N = 301) and then truncated to L = 6cm so that the windowing errors are
excluded from the spatial propagator result. Figure 3.12(a) shows the result for the
spectral propagator without angular restriction. The field contains some high spa-
tial frequency ripples due to under-sampling in the spatial frequency domain. Figure
3.12(b) is calculated using the spectral propagator with angular restriction, and the
ripples are reduced somewhat. Figure 3.12(c), which contains the result obtained
with the spatial propagator, is relatively smooth. All three results are quite simi—
lar to the reference in Figure 3.11. Figure 3.12(a) contains smaller aliasing errors
relative to the simulation result in non—attenuating media in Figure 3.5(a) because
the attenuation term S (km, kg, 2) is a low-pass filter that effectively attenuates the
high spatial frequency components in Hu(kx, Icy, 2). Since attenuation reduces the

highly oscillatory spectrum near kg + kg, = 1:2, aliasing errors in attenuating media

41

...

......
u ._
c,’ .

‘9
...
'.
‘.

   

    

e N... HL'~ ................

a I. ...................

8 1 ......... :

Q 0 -'- i .
0.6 " .. . _.

8 0.4 .W .. '« ,‘Q‘nftjf " 3‘. ~ I 30

N l”"\\"‘, as 25

g 0 ......... 20

o 3 ..... 15

2

Figure 3.11. The absolute value of the complex reference pressure generated by a 3cm x
3cm square piston in attenuating media (a = 1dB/cm/MHz). The excitation frequency is
1MHz. The reference pressure is computed with the fast nearfield method.

are less severe than those in non-attenuating media, therefore angular restriction is
not required in attenuating media. The spectral propagator is preferable in attenu-
ating media, because the spectral propagator generates acceptable results in a D x D
region while the spatial propagator only produces accurate results in an L x L region.
Meanwhile, the spectral propagator uses less computation time because fewer FFT
calculations are required.

Figure 3.13 shows the root mean squared errors obtained using the spectral prop-
agator (with and without angular restriction) and the spatial propagator. The pres-
sures are calculated in D x D transverse planes and then truncated to L x L planes
with L = D — 2a. The root mean squared error is then computed in the truncated
volume. In Figure 3.13, L ranges from 3cm to 12cm, and the three curves approx-
imately converge to the same value for large L. The difference between the result
from the unrestricted spectral propagator (dotted line) and the restricted spectral
propagator (solid line) is smaller in Figure 3.13 than in Figure 3.9(b), which implies

that angular restriction is less important in attenuating media.

42

     
  

Normalized pressure

 
     

5 —1.5 5
0 z(cm) x(cm) _3 0 z(cm)
(a) Simulated pressure obtained from the spec— (b) Simulated pressure obtained from the spec-

tral propagator without angular restriction. tral propagator with angular restriction.

Nonnallzed pressure
1?.
(.0
O

   

-1.5

x(cm) -3 0 z(cm)
(c) Simulated pressure obtained from the spa-
tial propagator.

Figure 3.12. Simulated pressure in an attenuating medium computed by the angular
spectrum approach using (a) the spectral propagator without angular restriction, (b) the
spectral propagator with angular restriction, and (c) the spatial propagator. All fields
are calculated in successive D x D transverse planes (D = 9cm) and truncated to L x L
(L = 6cm), where L = D —— 2a.

43

 

_o
O
4:.

 

------- LUnrestri'cted spectral propagator
'2, —Restricted spectral propagator
---Spatial propagator ,

 

 

.0
o
o)

 

 

Root mean squared error
9 o
O O
-8 N

 

 

 

3 5 7‘ a 11 13
L(cm)

Figure 3.13. Root mean squared errors evaluated for the pressure generated in attenuating
media (a = 1dB/cm/MHz) by a 3cm x 3cm uniformly excited square piston. The errors
are evaluated in three dimensional volumes where the lateral dimensions are L x L with
L = D — 2a. The markers on the curves indicate the corresponding results shown in Figure
3.12.

3.5 Angular spectrum approach for apodized source

3.5.1 Velocity source distribution

The distribution of the normal particle velocity on the piston surface also influences
the numerical accuracy of the angular spectrum approach. To demonstrate the effect
of apodization, a tapered window is applied to the aperture. The apodization func—
tion evaluated here is a product of quadratic polynomials in both the a: and the y

directions, where

um) = { [1 — (as/am [1— (y/bfl 1x1 < a. lyl < b (3,7)

0 |$|201|y|2b'

This source velocity has a peak at the origin and smoothly decays to zero at the edges
of the piston. Two-dimensional illustrations of the uniform velocity distribution used
in previous calculations and the tapered velocity distribution in Equation 3.27 are

shown in Figure 3.14(a) and Figure 3.14(b), respectively.

44

l

l 1' *

1 llll

' l 1)”
ll

llllllll” lllllilllll In ”III: II" A

(a) Uniform velocity distribution.

 

 
        

  
 

"If"

    

m . .~ ,m 3

11,, he ‘\ “u
"n I 'I N t s‘ " o“
”’1' 4"}..a" .i‘\| “ \\\

 
    
 

-. . u
,1: HI I, \ ‘
m, '1, I, I I I O“ ‘\‘\\\\‘
’14, I, 1'." as .\ ‘\ ‘\\

   

II
Ila, 1‘. .
”I”! 0, l ’I 0 ‘ I‘.‘\\‘ \\\““\l {
, I "flu“ s \\ “m \
4.0.9; \ \ \\“\\\\\\\‘

    
  
 
 

       

\
“hf \lilh
\‘.\ ' .‘
““‘lliitxxixmv‘.
\

    
   
   

\ .

\ . \

, at“ -\\\.\\\\\u\m“w

‘6w‘.‘:o‘lx“l\\\“lll\l“‘“lim\\\\§l“§li§.ww'

0 ‘ \\ \\ \\\ “V!“ “M.“
I,o‘15“ll1\\\l“ll(i_\\\i\\l\)“

“‘lllliltfiml‘

   
 

    

(b) Quadratically apodized velocity distribution.

Figure 3.14. A square piston (illustrated by the rectangle at the bottom) and the particle
velocity distribution on the piston surface for (a) uniformly excitation and (b) quadratically
apodization (indicated by the mesh on the top).

45

The pressure from the apodized piston is calculated with the point source super-
position of the Rayleigh-Sommerfeld integral, where each point source is weighted
by the value of the normal velocity at the sub-element center. In the numerical
evaluation of the Rayleigh-Sommerfeld integral, the reference field is computed with
5, 000 x 5, 000 uniformly excited subdivisions on the piston surface. This reference is
accurate to 7 digits, as determined from a comparison with the same result calculated
with 10,000 X 10, 000 subdivisions. Figure 3.15 shows the reference field in the my
plane at y = 0 in non-attenuating media. The pressure fields computed with the an-
gular spectrum approach are shown in Figure 3.16. Figure 3.16(a) and Figure 3.16(b)
show the results from the spectral propagator without angular restriction and with
angular restriction, respectively. The result obtained with the spatial propagator is
shown in Figure 3.16(c). All pressure fields are computed with D = 90m (N = 301)
and then truncated to L = 6cm so that windowing errors are excluded from the spa-
tial propagator result. Figure 3.16(c) is almost identical to the reference in Figure
3.15, while Figures 3.16(a) and 3.16(b) also closely resemble the reference except for
some small ripples caused by under-sampling in the spatial frequency domain.

Figure 3.17 shows the root mean squared errors of the pressures computed with the
spectral propagator (with and without angular restriction) and the spatial propagator.
The pressures are calculated in three dimensional volumes with a lateral extent of
D x D, then the volume is truncated in the lateral directions to L x L with L = D—2a,
and the root mean squared error is evaluated. The errors in all three curves are much
lower than those computed with the uniform velocity source in non-attenuating media
(Figure 3.9(b)) or attenuating media (Figure 3.13).

The two dimensional Fourier transform of the uniform velocity distribution in
Figure 3.14(a) is represented by a two dimensional sinc function with large sidelobes.
The spatial frequencies in the sidelobes of the sine function are under-sampled, there-

fore, aliasing errors are caused in the pressure field. In contrast, the apodized piston

46

  
      

L
U) \ \ o \\ 1‘ \ v
m 1 :\\‘\:\‘\‘\t\“‘\\\\\‘\\\“\\\\“‘\\“
\ \ \ \
\‘ \“ \\ \‘ \\‘ \
a) \x“\\\\§\\\“\\‘ \\‘,\
h 0 8
\ \ \ \\
o 6
. \\\ ‘\ ‘1 ‘\ ,1 \\
\ \\\\‘\ ‘\ ‘5‘
a) 0 4 ‘ ‘ “
_ \ x“
_N 0 2 25
L 3 1 5

1.5
x(cm) -3 0

Figure 3.15. The reference pressure calculated with the Rayleigh-Sommerfeld integral in
a non-attenuating medium, where the pressure is generated by a 3cm X 3cm square piston
with an apodized normal particle velocity distribution. The excitation frequency is 1MHz.

has a tapered velocity distribution, the spectrum of which has lower sidelobes and
is therefore less prone to high spatial frequency errors. Figure 3.18 shows the FFTS
of the two velocity sources in a line along the km direction at ky = 0, where the
spectrum of the uniform piston is represented by a solid line and the spectrum of the
apodized piston is represented by a dashed line. The first sidelobes for the uniform
velocity distribution and the apodized velocity distributions are observed at -6.8dB
and -11.2dB, respectively. Near k1 = k and kg = 0, where the spectrum becomes
evanescent outside of this range, the sidelobes are -17.7dB and -30.8dB for the uniform
and the apodized velocity distribution, respectively. Thus, by reducing the sidelobe
levels in the source spectrum, apodization achieves a significant reduction in the error

and eliminates the need for angular restriction.

3.5.2 Pressure source distribution

For a pre—computed pressure distribution in one transverse plane, the three dimen-

sional pressure field in successive parallel planes is computed either with equation 3.5

47

  
  

Normalized pressure

   

5
-1.5
x(cm) -3 0

0 5 z(cm)

(a) Simulated pressure obtained from the spec— (b) Simulated pressure obtained from the spec-
tral propagator without angular restriction. tral propagator with angular restriction.

Normalized pressure

 
   

5
-1.5
x(cm) -3 0 z(cm)
(c) Simulated pressure obtained from the spa—
tial propagator.

Figure 3.16. Simulated pressure in a non-attenuating medium generated by a 3cm x
3cm square piston with an apodized normal particle velocity distribution. The pressure
is computed by the angular spectrum approach using (a) the spectral propagator without
angular restriction, (b) the spectral propagator with angular restriction, and (c) the spatial
propagator. All fields are calculated in successive D x D transverse planes (D = 90m) and
truncated to L x L (L = 6cm), where L = D — 2a.

48

 

 

0.02 L _ - . .
------- Unrestncted spectral propagator

—Restricted spectral prepagator
0.01 5. - - - Spatial propagator

 

 

 

Root mean squared error

 

 

 

 

7 9
L(cm)

Figure 3.17. Root mean squared errors for the pressure generated in a non-attenuating
medium by a 3cm x 3cm square piston with an apodized normal particle velocity distribution.
The errors are evaluated in three dimensional volumes where the lateral dimensions are L x L
with L = D -— 20..

N

 

_\
O

 

:-Uniform velocity
- --Apodized velocity

 

 

 

_s
O
o

I
N

I
38

Spectrum (Log scale)
3 ES

 

 

 

 

10’6 -
-3 -1.5 0 1 .5 3
kx(normalized by k)

Figure 3.18. The spectrum of the uniform velocity source and the apodized velocity source
along along the kg; direction evaluated at kg = 0.

49

or with equation 3.11, thus, either the spectral propagator in equation 3.6 or the spa-
tial propagator in equation 3.14 is used to propagate the source field. In Figure 3.19,
the root mean squared errors for the simulated pressure using the velocity and the
pressure sources are compared. Figure 3.19(a) and Figure 3.19(b) show the results
obtained from the spectral propagator without angular restriction and with angular
restriction, respectively. Figure 3.19(c) shows the results from the spatial propagator.
The pressure field is generated by a 3cm x 3cm square piston with uniform normal
particle velocity in non-attenuating media. The source pressure at z = 0 is calcu-
lated by the fast nearfield method using 50 abscissas. In Figures 3.19(a), 3.19(b),
and 3.19(c), the root mean squared errors from the pressure source (dashed line) are
slightly smaller than the errors from the velocity source (solid line). The error gener-
ated by the pressure source is consistently smaller because the pressure distribution
varies smoothly at the edges of the piston, while the normal velocity distribution has
discontinuities in these locations. Aliasing of the velocity source spectrum is more
severe than the pressure source spectrum. Therefore, sampling the pressure source

field is advantageous for angular spectrum calculations.

3.6 Backward propagation

With a known pressure distribution in a plane, acoustic wave propagations can be
simulated with the angular spectrum approach in both forward and backward di-
rections orthogonal to the source plane plane. The backward propagation algorithm
is discussed in terms of ultrasonic wavefront design [35, 43], pressure field charac-
terization [36], and wavefront distortion compensation during phased array imaging
[47].

The spectral propagator for backward propagation is given by Equation 3.6. The
spatial propagator, however, differs from Equation 3.14 in that the Sign of the jkr

term is opposite because the direction of propagation is reversed. The formula for

50

 

—U0Hu

 

 

 

 

Root mean squared error

 

 

 

13
z(cm)

(a) Root mean squared errors computed us-

ing the spectral propagator without angular re-

striction.

 

 

 

 

 

 

 

0.2
E _UOHU
0 0.15- "'POHp
E
(U
3
8 0.1
C
m
0
E o 05»
£3
0 3 A T 11 13

z(cm)

(b) Root mean squared errors computed using
the spectral propagator with angular restric-
tion.

 

 

 

 

 

 

 

0.2
E —uo<8>hu
0 0.15 ---Do®hp
E
(U
3
3 01
C
on
O
E
*5 0.05
O
m \
0 3 1‘1 13

2(cm)

(c) Root mean squared errors computed using
the spatial propagator.

Figure 3.19. Root mean squared errors for the pressure generated in non-attenuating
media by a 3cm x 3cm square piston calculated with velocity and pressure sources. The
errors are evaluated in successive planes in different places in the z direction. The pressure
is computed by the angular spectrum approach using (a) the spectral propagator without
angular restriction, (b) the spectral prepagator with angular restriction, and (c) the spatial
propagator. All fields are calculated in successive D x D transverse planes (D = 9cm) and
truncated to L x L (L = 60m), where L = D — 2a.

51

the backward spatial propagator is

Az
7”.3

 

hpps, y, A2) = — (1 — jkr)ejk", (3.28)

 

where A2 = z —- 20 and r = V22 +y2+Az2.

A
I

pressure

.0
‘1‘

 

” Owl“
inn"; I will“

Normalized

f]!
I | I"
4,11,72,46.“ “‘0”I’.’7I\'I.

”‘9

 

Figure 3.20. The pressure field in my plane at z = 15cm generated by a 3cm x 3cm square
piston in non-attenuating media. The excitation frequency is 1MHz. The pressure field
is computed with the fast nearfield method and the result is normalized to the maximum
pressure amplitude in the z = 15cm plane.

Figure 3.20 shows the pressure generated by a 3cm x 30m square piston in a non-
attenuating medium in the my plane at z = 15cm. The pressure field is truncated
by a 9cm x 9cm square window. Using this pressure distribution as the source,
the pressure field is propagated forward and backward to parallel planes, and the
results are shown in Figure 3.21. The pressure in Figure 3.21(a) is computed with
the spectral propagator without angular restriction using N = 301. Even without
angular restriction, the backward propagated field in the range from z = 0 to 15cm is
very smooth, and the forward propagated field from 2 = 15cm to 30cm demonstrates
fewer aliasing errors compared to Figure 3.4(a). This improvement occurs because the
pressure source at z = 15cm contains less high spatial frequency components than the

velocity source at z = 0, and most under-sampled spectral components in the spectral

52

  
     
   

099994
memo

Normalized pressure
.5
'01

z(cm)

(a) Simulated pressure obtained from the spectral propaga-
tor with N=301 without angular restriction.

     
  

CD
5
8 1
gas -
0.6 30
8 0.4
g 0.2
g o 20
54.5 15
Z
z(cm)

(b) Simulated pressure obtained from the spatial propagator
with N=301.

Figure 3.21. A cross section of the simulated pressure generated by a 3cm x 3cm square
piston in non-attenuating media. The 3D pressure is calculated from a windowed pressure
distribution at z = 15cm and propagated forward and backward with the angular spectrum
approach using (a) the spectral propagator with angular restriction evaluated in a 301 x 301
grid, and (b) the spatial propagator truncated by a 9cm x 9cm window and evaluated in a

53

propagator are filtered out after the spectral propagator and the source spectrum are
multiplied. In the region between 2 = O and 15cm, some of the sidelobes are absent
from Figure 3.21(a) and Figure 3.21(b) that are present in the reference field in Figure
3.2, and the pressure field adjacent to the transducer also appears to be smoother than
the reference. This is because of the absence of high spatial frequency components in
the source spectrum at z = 15cm, and the restricted bandwidth ( k3; + k; S k) of
the spectral propagator.

Figure 3.21(b) is computed from a forward and backward spatial propagator that
is truncated to D = 9cm in both the :c and the y directions and transformed with
a 301 x 301 point 2D FFT. In Figure 3.21(b), the forward propagated field from
2 = 15cm to 30cm is corrupted with oscillatory errors throughout the transverse
direction. Unlike Figure 3.3(a), where the errors are restricted in a band near the edge
of the computational grid, the errors in Figure 3.21(b) span the entire computational
grid. This error is caused by the circular convolution between the pressure source
and the spatial propagator, both of which are truncated by a D x D window and
both have non-zero values throughout the D x D area. Since the non-zero part of the
velocity source is confined within a 2a x 2a window, the L x L area in Figure 3.3(a) is
not affected by this circular convolution error, where L = D — 2a. In Figure 3.21(b),
where the pressure source in a D x D window is used, the circular convolution error
affects the entire D x D area. In Figure 3.21(b), the circular convolution error also
occurs in the range from 2 = 0 to 15cm, but the error in this region is less severe.

Figure 3.22 displays the reference pressure (solid line) and the backward propa-
gated pressures along a line in the :1: direction evaluated at y = 0 in the plane that
contains the transducer (2 = 0). The dashed line represents the pressure computed
with the spectral propagator without angular restriction, and the dash-dot line repre-
sents the pressure computed with the spatial propagator. Due to the limited spatial

frequency bandwidth of the source pressure spectrum, the results from the spectral

54

and spatial propagators are both smoother than the reference field. The high spatial
frequency components are not retained during backward propagation. Figure 3.22
also shows that the pressure computed from the spatial propagator has large errors

near the edge due to the circular convolution error.

 

 

 

 

 

 

 

 

1.5 . e
—Reference
---P H
91.2l 0 f1 ~
'''''
a0.9' ." ‘W‘c’ . ~
13
3
=3 0.61 l
E
O
2 0.3- l
0 ~ ’ ' ,
-4.5 -3 -1 5 0 1 5 3 4 5

Figure 3.22. The simulated pressures along the a: direction evaluated at y = z = 0.
The solid line is the reference pressure, the dashed line is the pressure computed with the
spectral prepagator without angular restriction (POHp), and the dash-dot line is the pressure
computed with the spatial propagator truncated by a 9cm x 90m window (po 8) hp). The
source plane is located in the plane at z = 150m, and the pressure is backward propagated
with the angular spectrum approach.

55

CHAPTER 4

Ultrasound Phased Array Simulations

4.1 Introduction

Ultrasound phased arrays often contain many piezoelectric transducers. Each trans-
ducer has an electrode attached, and a voltage signal that excites the transducer is
applied to the electrode. Applying signals with different time delays to the array
transducers can focus and steer the ultrasound beam. The abilities of phased ar-
rays to electronically focus, swiftly scan, and generate multiple focal spots enable
phased arrays to treat large tumors in deep sites, which bring enormous advantages
to thermal therapy.

Computer simulations provide valuable information for the design of phased arrays
and for patient treatment planning. The size of the array aperture, the operating fre-
quency, and the size and number of array elements heavily influence the performance
of phased arrays for therapeutic applications. A large aperture size provides a wide
acoustic window and thus reduces the amount of ultrasound energy that accumulates
in front of the focus. A large number of transducer elements across a broad aperture
ensures sufficient power delivery within the target. The size of the individual trans-
ducer elements and the center-to—center spacing between adjacent elements needs to
be sufficiently small so that unwanted grating lobes are reduced. The attenuation

of ultrasound is a function of frequency, which determines the penetration depth of

56

ultrasound.

The three dimensional pressure field is computed by superposing of the pressures
produced by all transducer elements combined with the corresponding weights and
phase shifts. The transducer sources in ultrasound arrays are modeled by the point
source superposition of the Rayleigh-Sommerfeld integral [11, 28, 48], the rectangular
radiator method [49, 50], the spatial impulse response method [51, 52], and other
integral approaches based on these methods.

In this chapter, several 1.5D and 2D ultrasound phased array configurations are
first introduced. A single focus beamforming method and a few multiple-foci beam-
forming methods are then described. The pressure field generated by a 2D pla—
nar phased array and two concave phased arrays are simulated with the Rayleigh-
Sommerfeld integral and the fast nearified method. The computation time and nu-

merical error are compared among these numerical results.

4.2 Phased array designs

4.2.1 Curvilinear array

A curvilinear array is a concave one dimensional array. The curvilinear array in Figure
4.1 comprises rectangular strip elements, and the normal of each element points to
the geometrical focus of the array. The array in Figure 4.1 contains 51 rectangular

elements that are 0.135cm wide and 5cm high. The curvature of the array is 5cm.

 

Figure 4.1. A curvilinear array comprised of 51 rectangular elements.

57

4.2.2 Two dimensional planar array

A two dimensional planar array is illustrated by Figure 4.2, where each small square
represents a transducer element. In Figure 4.2, 1024 square elements are packed into
a 7.31cm x 7.31cm area. The size of each element is 0.18cm x 0.18cm, and the kerf
between adjacent elements is 0.05cm. For a 1MHz excitation frequency, the center-

to—center spacing between adjacent elements is 1.53/\, where A is the wavelength.

Figure 4.2. A two dimensional planar array comprised of 1024 square elements.

4.2.3 Spherically focused array

A spherically focused array that consists of densely packed circular elements on a
spherical shell [53] is shown in Figure 4.3. Six longitudinal axes are filled with ele-
ments, and these axes divide the array into equal sectors that demonstrate hexagonal
symmetry. Each element is rotated such that the normal is pointed at the geometri-

cal focus of the array. The array in Figure 4.3 contains 1441 circular elements with

58

0.24cm diameter. In general, concave arrays have higher focal gain than planar arrays

and therefore are preferred in thermal therapy applications.

4.2.4 Spherical-section array

A spherical-section array is comprised of square elements distributed on a spherical
shell. The surface of this array is a spherical section truncated by a rectangular
window [14]. Figure 4.4 shows a spherical-section array comprised of 1444 square
elements, and the size of each square element is 0.24cm x 0.24cm. The geometrical
focus of this array is located at z = 12cm, the angle defined by the aperture with
respect to the geometrical focus is 7r / 3 in each lateral direction, and the aperture

extends about 12cm in each direction.

4.3 Ultrasound phased array beamforming

Beamforming technique is applied to ultrasound phased array to produce directional
beams. Electronic focusing and beam steering are realized by adjusting the time
delay of the excitation signals applied to array elements. As a result, constructive and
destructive interferences maximize the pressure field at one or more spatial locations

and minimize the pressure field at other locations.

4.3.1 Single focus

In homogeneous media, a single focus can be generated with an ultrasound phased ar-
ray by applying ejklrl) - rf I to each array element, where k = 27r f / c is the wavenum-
ber, and [r6 — r f] is the distance from the center of a transducer element r6 to the
focal point r f. This focusing scheme assumes that the contributions from all array
elements form a spherical wavefront. The pressure generated by each element has the
same phase at the focal point; therefore, constructive interference is realized at that

point.

59

00.00
assesses...
ooomoomooowoowooMo

OOOWOOOOOOOOOOWOOWOOO
OOOOOOOOOOOCOOOOOOOOOOOO
WoemmommomowmoWowowoowowwwooowoo
oooooooooeooooooooa
oooooooooo oooooooooo
coco o 0000 o o 000 o o 0000 0000
o o 0000 o 000000 0000 o o 0000 o
sooooo 000000 0 o o 0 000000 0 coco
5.00000 0 0 000000000 0 0 0000.00
oomooeemomomomooooeoooowowowoweooowoo
omoooooo eeooooooooeooooooooeoooooooowo
oooooooomooeeooooooeooooooeeoooooooooooo
o 0 000000 e 0 000000000 0 o oo o oo

wouowe 0

000000 0

000000 0

000000 0

000000 0

000000 0

00 000 o

 

(a) Top view

            

ooooo

mwmmmmmmm“

bonuses m u .....
eases . . ....

\

(b) Side view

Figure 4.3. A spherically focused phased array comprised of 1441 circular elements: (a)

top view and (b) side view.

60

a) Top view

(

 

 

(b) Side view

section array comprised of 1444 square elements: (a) top view and

Figure 4.4. A spherical-

(b) side view.

61

The phase conjugation method is another related focusing scheme [54]. In this
method, the phase angle of the excitation signal applied to a transducer element is
equal to the phase angle of the complex pressure produced by that element at the
focal point multiplied by negative one. The excitation signal is thus equivalent to
the conjugate of the normalized complex pressure produced by a transducer element
at the focal point. The phase conjugation method ensures that pressures generated
from all elements have zero phase in a specified location, which produces constructive

interference at the focal point.

4.3.2 Sector-vortex focus

A sector-vortex array [16, 17] is formed by dividing an annular transducer or a con-
centric ring array into sectors. Each of these sectors is multiplied with a rotational
signal in order to generate a ring-type focus. Similar schemes can be applied to planar
arrays or spherically focused arrays. As shown in Figure 4.3, a spherically focused
array is divided into six sectors with hexagonal symmetry. The rotational phase 679’”

is applied to all elements in the ith sector,
d,- = 27rMi/6, M = 1,2, or 3. (4.1)

In Equation 4.1, M is the number of phase rotations along the entire array aperture
over a period of 27r. The value of M is bounded by the number of sectors.

Figure 4.5 shows the contours of the pressure fields generated by a spherically
focused array using the sector-vortex excitation. The phased array is comprised of
661 circular elements and the diameter of each element is 0.24cm. The aperture of
the phased array is 8cm in diameter, and the radius of curvature is also 8cm. The
value of M = 3 is used, and the excitation signals applied to the six sectors are
[87”, 1, 677', 1, 677T, 1] in sequence. Since the phases of the excitation signals between
two adjacent sectors are opposite, the pressure is canceled in the plane of symmetry

between two sectors. Constructive interference is generated between adjacent planes

62

of symmetry, and focal spots are produced near the geometrical focus of the array.
Figure 4.5(a) shows the contour of the pressure field within a 2cm x 2cm area evaluated
in the focal plane at z = 8cm. Six separate focal spots are formed, and each focal
spot is 0.8cm from the z axis. Six hot spots produced by the grating lobes also
appear in Figure 4.5(a). Figure 4.5(b) and Figure 4.5(c) show the pressure contours
within a 6cm x 6cm area in the pre-focal plane at z = 6cm and the post-focal plane
at z = 10cm, respectively. Figures 4.5(b) and 4.5(c) demonstrate similar contours
in a broader area relative to Figure 4.5(a), and the maximum pressure amplitudes
in Figures 4.5(b) and 4.5(c) are 5dB and 10dB lower than that in Figure 4.5(a). In
Figure 4.5(d), the pressure contour in the yz plane at a: = 0 shows two focal spots
and an empty area between the focal spots where the pressure field is canceled.

The sector-vortex beamforming method produces multiple focal spots and covers
a larger area than a single focus beam pattern. However, the standard focal pattern
produced by the sector-vortex array formula lacks flexibility. Once the structure of
the phased array and the mode number M are chosen, the locations of the focal spots
are fixed. To change the focal depth or the spacing between the focal spots requires

a reconfiguration of the array or a more flexible focusing strategy.

4.3.3 Mode scanning

Mode scanning [55] produces multiple focal spots simultaneously according to the
symmetry of a phased array. McGough et al describes how this method generates
two focus or four focus modes using a two dimensional planar phased array [55].
Specifically, the elements in the planar array are divided into four groups by the
horizontal and the vertical lines across the center of the array. In Figure 4.6, a
planar array comprised of 100 square transducers is shown, and the four elements
belonging to four different groups are labeled with 1, 2, 3, and 4. Rotational signals

[1, 677‘, 1, 67“]e_j‘p are applied to elements 1, 2, 3, and 4 in sequence. The resulting

63

   

—1 —o.5 o 0.5 1 o
ylcm) ylcm)

(a) Pressure in the focal plane at z=8cm. (b) Pressure in the pre-focal plane at z=6cm.

3

   

10 12

133 —2 -1 8
z(cm)

o
ylcm)

(0) Pressure in the post-focal plane at z=10cm. (d) Pressure in the axial plane at x=Ocm.

Figure 4.5. Contours describing the pressure field generated by a 661-element spherically
focused array excited by sector-vortex phases with mode number M=3. The array has an
8cm aperture diameter and an 8cm radius of curvature. The focal spots are located in the

my plane at z = 8cm.

64

 

 

IIHIIIINI
IIEIIIIHI
Ila-Ell.-

Figure 4.6. Illustration of a group of four symmetric elements in a two dimensional planar
phased array.

fields are canceled along the horizontal and the vertical planes of symmetry, and four
focal spots are generated in the designated locations.

Similarly, a six focus mode can be generated with a spherically focused array using
the mode scanning technique. The mode pattern is symmetric with respect to the
projection of the hexagonal axes in the transverse planes, and the fields are canceled
in the planes of symmetry that are perpendicular to the transverse planes. This is
achieved by applying excitation signals with rotational phases to each group of array
elements that are determined by the symmetry of the array, as for the elements with
labels from 1 to 6 in Figure 4.7. The excitation signals applied to these elements
are [1,67“, 1, 67”, 1, 67”]6-790. Let ff,- represent the complex pressure generated by
element 1' at an arbitrary field point. In order to maximize the sum of the pressure

fields at that point from all six elements, (,0 satisfies

so=arg<7£—7z§+i§—E+EE—Féi (42)

Similar to the sector-vortex method, the mode scanning method cancels the fields

65

Figure 4.7. Illustration of a spherically focused phased array. The elements labeled from
1 to 6 are symmetric with respect to the hexagonal axes.

in six planes of symmetry, which reduces intervening tissue heating for deep regional
hyperthermia. However, the mode scanning method provides a more flexible focusing
scheme than the sector-vortex method. Figure 4.8 and Figure 4.9 demonstrate the
contours of the pressure field generated by a 661 element spherically focused array
with a six focus mode. Figure 4.8(a) and Figure 4.8(b) show the pressure contours
in the focal plane at z=8cm and the axial plane at :r=0cm, respectively, where the
distance from each focal spot to the z axis is 0.5cm. Figures 4.9(a) and 4.9(b) demon—
strate two pressure contours similar to Figure 4.8(a) and Figure 4.8(b), except that
the distance from each focal spot to the z axis is lcm. As shown in Figure 4.8 and
Figure 4.9, the spacing between focal spots produced by mode scanning can be ad-
justed arbitrarily. Scans with different modes can be designed fill in a large target
region.

Steering of focal spots generated by mode scanning can be achieved by shifting

66

   

0 8 12
31(le z(cm)
(a) The pressure in the focal plane evaluated (b) The pressure in the axial plane at x=0crn.

at z=8cm.

Figure 4.8. The contours of the pressure field generated by a 661 element spherically
focused array focused with the mode scanning method. The array aperture diameter is
8cm, and the radius of curvature is 8cm. The focal spots are located in the my plane at
z = 80m, and each focal spot is located 0.5cm from the z axis. The pressure is simulated
with the fast nearfield method.

   

4 6 8 10 12
y(cm) z(cm)

(a) The pressure in the focal plane evaluated (b) The pressure in the axial plane evaluated
at z=8cm. at x=0cm.

Figure 4.9. The contours of the pressure field generated by a 661 element spherically
focused array focused with the mode scanning method. The array aperture diameter is
8cm, and the radius of curvature is 8cm. The focal spots are located in the my plane at
z = 8cm, and each focal spot is located 1cm from the z axis. The pressure is simulated
with the fast nearfield method.

67

 

10 12

8
z(cm)

Figure 4.10. Contour plot describing the pressure field generated by a 661 element spher-
ically focused array focused with the mode scanning method. The aperture diameter of the
array is 8cm, and the radius of curvature is 8cm. The focal spots are steered off axis in the my
plane at z = 8cm, and the center of the focal pattern is located at (m, y, z) = (0, 0.80m, 8cm).

the phases that are applied to the excitation signals. According to Fraunhofer ap—
proximation, the pressure in the farfield is proportional to the Fourier transform of
the normal particle velocity distribution on the transducer or array aperture [30].
The forward Fourier transform maps variables (m’, y’) to (km’/z, ky’/z), and the in—
verse Fourier transform maps (km’/z,ky’/z) to (m’,y’), where (m’,y’,0) and (m,y,z)
are the coordinates of the transducer element and the observation point, respectively.
To steer the center of the six focal spots from (0,0, 2) to (ms,ys, z), the signal ap-
plied to each of the elements is multiplied by e—jk(m3mb + ysyb)/z in addition to the
phase term computed for the mode scanning, where (mb,y6) denotes the center of
each circular element. This is because a coordinate translation in the spatial domain
(m’ — 936, y’ — 90) is equivalent to a phase shift in the spatial frequency domain. Figure
4.10 shows a pressure contour similar to Figure 4.8(b) except that the focal spots are

moved 0.8cm in the positive y direction.

68

4.4 Pressure simulations for phased arrays using integral ap-

proaches

4.4.1 Phased array models

Three types of ultrasound phased arrays are simulated. The first array is a two di-
mensional planar array comprised of 1024 square transducers, where each transducer
is 0.18cm x 0.18cm, and the kerf between adjacent elements is 0.05cm. The sec-
ond array is a spherical-section array comprised of 1444 square elements, where each
transducer is 0.24cm x 0.24cm. The spherical-section array is geometrically focused at
(0, 0, 12cm), and the opening angle from the geometrical focus to the lateral extent of
the array is 7r / 3 in the m and the y direction. The third array is a spherically focused
array comprised of 1441 circular elements, where the diameter of each transducer is
0.24cm and the center-to—center spacing between adjacent elements is approximately
0.315cm. The aperture diameter of the spherically focused array is 12cm, and the
radius of curvature is 12cm. All three arrays are centered at the origin of the coordi-
nate system, and the direction normal at the center of each array is in the positive 2
direction. The excitation frequency is 1MHz, the speed of sound is 1500m/s, and the
attenuation coefficient is a = 1dB/cm/MHz in all simulations.

The reference pressure fields produced by the planar array, the spherical-section
array, and the spherically focused array are calculated by the spatial impulse response
method [20, 21, 22] using 1000 Gauss abscissas on each element. The pressure is also
computed with the fast nearfield method [23, 24, 26] and the point source superposi-
tion of the Rayleigh-Sommerfeld integral [18, 19].

4.4.2 Computational coordinate system

The computational domain is a 16.2cm x 16.20m X 12cm rectangular volume. With

a sampling interval of 6 = 0.075cm (A / 2), the computational volume is discretized to

69

a 217 x 217 x 161 grid in a Cartesian coordinate system.

To calculate the pressure produced by an individual element in the global coor-
dinate system, the computational grid is translated with respect to the center of the
element, then multiplied by a rotation matrix that describes the grid coordinate rela-
tive to the orientation of the element. The rotation matrix is formulated as a function

of the Euler angles (6,11), (b),

cos¢c0s6 —- sinqzisinipsinél sinqficosfi +cos¢sin1psin6 —cosr/)sin6
ROT = — sin 45 cos 2p cos 45 cos 1]) sin 1p ,
cos qb sin 6 + sin 45 sin 11) cos 6 sin ([5 sin 0 — cos g1) sin 1,!) cos 6 cos 1,!) cos 0
(4.3)

where 9, 11) and Q5 are the rotation angles of the element along the y, m, and z axes,
respectively. Evaluating the matrix-vector product for this rotation matrix and the
coordinate vector is equivalent to rotating this point by (,25 degrees with respect to
the z axis first, then It degrees with respect to the m axis, and finally 0 degrees with

respect to the y axis.

4.4.3 Simulation results

An example of the simulated pressure generated by the planar phased array is shown
in Figure 4.11, where the array is electronically focused at (0,0,10cm). The pressure
field is calculated in a 16.2cm x 16.2cm x 12cm volume. The computational volume
extends from —8.1cm to 8.1cm in the m and y directions and extends from 4cm to
16cm in the z direction. Figure 4.11(a) shows the pressure field in the yz plane at
m = 0 and Figure 4.11(b) shows the pressure field in the my plane at z = 10cm.

The pressure fields generated by the spherical-section array and the spherically
focused array are shown in Figure 4.12 and Figure 4.13, respectively. The elements
in these two phased arrays are excited in phase to produce a geometrical focus at
(0,0,12cm). For both arrays, the pressure field is computed in a 16.2cm x 16.2cm x

12cm volume that extends from —8.1cm to 8.1cm in the m and y directions and from

70

Normalized pressure

    

_ ~ 6
y(cm) 4 -8 4 z(cm)

(a) Simulated pressure in the yz plane at m = 0.

  

99.09..
thnoo

Normalized pressure

mo

-4

4 .
x(cm) '8 '8 y(cm)

(b) Simulated pressure in the my plane at z = 10cm.

Figure 4.11. The simulated pressure generated by a planar array with 1024 square elements
(a) in the yz plane at m = 0 and (b) in the my plane at z = 10cm. The pressure is normalized
by the peak pressure amplitude in the three dimensional volume.

71

6cm to 18cm in the z direction. Figure 4.12(a) and Figure 4.13(a) show the pressure
fields in the yz plane, whereas Figure 4.12(b) and Figure 4.13(b) show the pressure
fields in the my plane. Figures 4.12 and 4.13, when compared with Figure 4.11, shows
that the pressure fields generated by the concave arrays, i.e., the spherical-section
array and the spherically focused array, have lower grating lobes than that generated

by the planar array.

4.4.4 Computation time and numerical error

The pressure field generated by the planar array, the spherical-section array, and
the spherically focused array is computed with the fast nearfield method and the
point source superposition of the Rayleigh-Sommerfeld integral in a 217 x 217 x 161
grid. The computation times and the root mean squared error (RMSE) values for
the simulated pressure are listed in Table 4.1. All simulations were performed on a
2.4GHz Pentium 4 PC (1GB RAM) running Windows XP operating system. The
routines were written in C language, compiled by Microsoft Visual C/C++ version
7.0 and called by MATLAB 7.1 as MEX files.

The integrals in the fast nearfield method are evaluated with Gauss quadrature
for rectangular elements and with Riemann’s mid-point rule for circular elements. In
the point source superposition method, the Rayleigh-Sommerfeld integral is evalu-
ated with a 2D integration over the transducer surface using the mid-point rule. The
numerical accuracy of these integral methods primarily depends on the number of
abscissas used. In the first line of Table 4.1(a), 4.1(b) and 4.1(c), the n term in the
parenthesis denotes the number of abscissas for the FNM or the Rayleigh-Sommerfeld
integral. The values for n in Table 4.1 represent the smallest number of abscissas re-
quired for each method to achieve a root mean squared error below 0.01 in the 3D
volume. In Table 4.1, the computation time of the fast nearfield method is much

shorter than that of the point source superposition of the Rayleigh-Sommerfeld inte-

72

Normalized pressure

    

_ 8
y(cm) 4 -8 6 “W"

(a) Simulated pressure in the yz plane at m = 0.

      
 

9.0.0.0..
memos

 

Normalized pressure

we

4 _
x(cm) _8 '8 y(cm)

(b) Simulated pressure in the my plane at z = 12cm.

Figure 4.12. The simulated pressure generated by a spherical-section array with 1444
square elements (a) in the yz plane at m = 0 and (b) in the my plane at z = 12cm. The
pressure is normalized by the peak pressure amplitude in the three dimensional volume.

73

Normalized pressure

     

_ 8
y(cm) 4 —8 6 z(cm)

(a) Simulated pressure in the yz plane at m = 0.

Normalized pressure

 

x(cm) _8 "8 y(cm)

(b) Simulated pressure in the my plane at z = 12cm.

Figure 4.13. The simulated pressure generated by a spherically focused array with 1441
circular elements (a) in the yz plane at m = 0 and (b) in the my plane at z = 12cm. The
pressure is normalized by the peak pressure amplitude in the three dimensional volume.

74

Table 4.1. The computation time and the root mean squared error (RMSE) for the simu-
lated pressures generated by (a) a planar array with 1024 square elements, (b) a spherical-
section array with 1444 square elements, and (c) a spherically focused array with 1441 cir-
cular elements. The pressures are simulated with two methods: the fast nearfield method
(FNM) and the point source superposition of the Rayleigh-Sommerfeld integral (Rayleigh-
Sommerfeld). The simulated pressures are compared with the reference pressure computed
with the spatial impulse response method.

Method | FNM (n=2) | Rayleigh-Sommerfeld (n=4 x 4)
Time (min) | 346.3 I 646.1

RMSE | 0.0042 | 0.0087
(a) Planar array

Method | FNM (n=2) | Rayleigh-Sommerfeld (n=5 x 5)
Time (min) | 948.36 | 1650.3

RMSE ] 0.007 | 0.0064
(b) Spherical-section array

Method | FNM (n=3) | Rayleigh-Sommerfeld (n=5 x 5)
Time (min) | 256.5 | 1524.2

RMSE | 0.0086 | . 0.0079
(c) Spherically focused array

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

gral, because a 1D integration is evaluated in the former whereas a 2D integration is
used in the latter. This is consistently observed in the simulations evaluated for each

of the array structures.

75

CHAPTER 5

Angular Spectrum Approach for
Ultrasound Phased Array Simulations in

Linear and Homogeneous Media

5.1 Introduction

Pressure fields generated by ultrasound therapy arrays are typically calculated by su-
perposing the fields from individual transducer sources. Traditionally, these sources
are modeled with integral approaches at each grid point. Therefore, the simulation
time is proportional to the number of array elements multiplied by the size of the
computational grid. These simulations are relatively slow because the fields are re-
peatedly computed at each point for every transducer source. In contrast, the angular
spectrum approach [30] rapidly computes pressures in parallel planes with a signifi-
cant reduction of computation time. The field generated by each transducer source
is calculated at most once in a single transverse slice of the computational grid, and
the pressure field throughout the rest of the grid points is calculated in the spatial
frequency domain.

The angular spectrum approach is widely used for phased array simulations. For
example, the pressure fields generated by concentric ring arrays and sector vortex

arrays have been compared with experimental results [36], and the angular spec-

76

trum approach has also been applied to high intensity focused ultrasound (HIFU)
simulations [56]. Zemp et al [44] compares the sampling of the spatial and spectral
propagators for phased array simulations and then derives the maximum unaliased
spectral sampling rate for the spectral propagator. Other important issues that im-
pact the accuracy of thermal therapy simulations remain unsolved; therefore, more
thorough evaluations of the angular spectrum approach are needed.

In this chapter, the performance of the angular spectrum approach is first com-
pared for velocity and pressure source plane inputs. Second, the effect of the location
of the source plane and the size of the window that truncates the pressure source
plane are determined. Third, the impact of the numerical accuracy of the 2D pressure
source plane on the 3D output field computed with the angular spectrum approach
is evaluated. Finally, the significance of the angular spectrum approach in reducing

the computation time is demonstrated.

5.2 Comparison between the pressure and velocity sources

Although angular spectrum calculations with the pressure source (Equation 3.5) and
the particle velocity source (Equation 3.7) are analytically equivalent, the numerical
errors differ. To illustrate the difference between these two approaches, the pressure
generated by the planar array with 1024 square elements in Figure 4.2 is simulated
with the angular spectrum approach using pressure and particle velocity sources.
In these simulations, the particle velocity source distribution is uniform across each
element on the array aperture. The source pressure field ideally extends to infinity
in both lateral directions; however, for these computer simulations, the pressure field
is truncated by a 16.20m X 16.2cm square window. The source pressure is calculated
with the fast nearfield method using 20 abscissas for each integral.

Figure 5.1 shows the reference axial pressure at m = y = 0 (solid line) and the axial

pressure computed with the angular spectrum approach using the normal particle

77

velocity (dash-dot line) and the pressure source at z = 0 (dashed line). In each
lateral direction, the spatial sampling interval is 6 = 0.075cm, and N = 512. The
pressure fields are normalized by the maximum amplitude of the reference pressure.
Figure 5.1 shows that the fields obtained from both the pressure source and the normal
velocity source match the reference closely near the focus. The field computed with
the pressure source is slightly larger than the reference in front of the focus and
slightly smaller than the reference beyond the focus. The field computed with the
normal velocity source is smaller in amplitude compared to the reference in front of
the focus and slightly larger than the reference beyond the focus. Overall, the axial
pressure computed with the normal velocity source has a somewhat larger error than

that obtained with the pressure source.

 

 

 

 

 

 

 

 

1_5 . . .
—Reference
m _ -----Velocity source
5 1'2 ---Pressure source ‘
t
50.91
1:
£1
'7; 0.6-
E .
o 1‘
z 0.3' ‘ 3"
o “ I . . .
4 8 10 12 14 16
z(cm)

Figure 5.1. The simulated axial pressures generated by a planar phased array located in
the my plane at z = 0cm and electronically focused at (0,0,10cm). The reference pressure
calculated by the spatial impulse response method is represented by the solid line, the pres-
sure computed with the angular spectrum approach using the velocity source is represented
by the dash—dot line, and the pressure computed with the angular spectrum approach using
the pressure source is represented by the dashed line.

The root mean squared errors calculated in transverse planes are plotted as a
function of the z coordinate in Figure 5.2. In Figure 5.2(a), the spatial sampling

interval is 6 = 0.075cm (A/ 2 at 1MHz), and the spectral propagator is discretized in

78

an N x N points grid with N = 512. In Figure 5.2(b), the spatial sampling interval
is 6 = 0.0375cm (/\/4 at 1MHz), and the spectral propagator is discretized in an
N x N points grid with N = 1024. In both Figure 5.2(a) and Figure 5.2(b), the
errors decrease monotonically as 2 increases, and the errors produced by the velocity
source (dash-dot line) are higher than those produced by the pressure source (dashed
line) everywhere. By decreasing 6 from A / 2 to A / 4, the maximum root mean squared
error for the result obtained from the velocity source is reduced from 0.098 in Figure
5.2(a) to 0.041 in Figure 5.2(b), and the maximum root mean squared error for the
result obtained from the pressure source is reduced from 0.055 to 0.007.

The spatial sampling interval 6 is an important parameter in angular spectrum
simulations. Undersampling in the spatial domain produces aliased spectra in the
spatial frequency domain. In Figure 5.2, the large errors produced by the particle
velocity or pressure source plane located at 20 = 0 are largely due to aliasing. In
the phased array model evaluated here, the velocity distribution on each transducer
element is represented by a 2D rect function. The discontinuities at the edge of
each element and at the edge of the array aperture introduce high spatial frequency
components that are inherently aliased. Thus, the velocity source encounters some
sampling difficulties that cause aliasing. In contrast, the changes in the pressure
distribution are less abrupt. The angular spectrum of the pressure source contains
fewer high spatial frequency components; therefore, aliasing errors in the pressure
source plane are less severe. Figure 5.3(a) and Figure 5.3(b) show the angular spectra
of the particle velocity source and the pressure source located at z = 0cm, respectively.

The angular spectrum for either source is obtained by sampling the source with

6 = A/ 2 (solid line) and 6 = A/4 (dash-dot line), and then applying a 2D FFT.

 

The angular spectrum of the velocity source is multiplied by k/ (j \/k2 — k3; — kg) to
compensate for the amplitude difference with respect to the angular spectrum of the

pressure source. Figure 5.3 shows that the angular spectra for the velocity source

79

 

.9
_A
J

 

 

 

 

 

 

 

‘1, --'-'Velocity source
'5 "s. - - -Pressure source
5 0.08- ‘-,. 1
E A"!
g 0.06; 1
CT \ t
(D x t
c \s ‘v‘.
8 0.04-
8 002. “4“ ........... 4
a: ‘ ................... -
0 1 L
4 6 8 10 12 14 16
z(cm)

(a) 6 = A/2 and N = 512.

 

9
o
or

 

----- Velocity source
-- -Pressure source

 

 

 

9.0.0
00
row?
I
I
I
I
I

Root mean squared error
.0
O
_‘L

 

 

I
Q
0‘.
~l
V.
§,.
I
......
b .I
.............
---
--
--------
-----
------------

 

4 5 81101121416
z(cm)
(b)6=A/4 andN=1024.

Figure 5.2. Root mean squared errors evaluated in transverse planes for z ranging from
4cm to 160m. The errors are obtained from pressure fields generated by the 1024 element
planar phased array, which is electronically focused at (0,0,10)cm. The pressure is calcu-
lated with the angular spectrum approach using a velocity source plane (dash-dot line) and
a pressure source plane (dashed line). The velocity and pressure source planes are both
located at z = 0cm and truncated with a 16.2cm x 16.2cm square window. The source field
is sampled at an interval of (a) 6 = 0.075cm and (b) 6 = 0.0375cm. The FFT calculations
are evaluated in N x N grids with N = 512 and N = 1024 in (a) and (b), respectively.

80

plane sampled with 6 = A/ 2 and A/ 4 differ by a significant amount in Figure 5.3(a),
while the spectra for the pressure source plane sampled with 6 = 2V 2 and A/ 4 are
very similar in Figure 5.3(b). Therefore, the velocity source contains more aliasing
errors than the pressure source when sampled with 6 = /\ / 2. For either source plane,

the aliasing error is reduced by increasing the spatial sampling rate.

5.3 Optimal parameters for the source field

When the spectral propagator Hp(km, ky, A2) for the pressure source in equation
3.6 is used for angular spectrum simulations, the source pressure field is truncated
by a rectangular window in the m and y directions. The size of the window and
the location of the pressure source plane both influence the accuracy of the angular
spectrum approach. These parameters are evaluated for the planar phased array in
Figure 4.2, where the spatial sampling interval is 6 = A/ 2, and the N x N grid for
the 2D FFT is evaluated with N = 512.

To demonstrate the impact of source plane location on angular spectrum calcula-
tions, the root mean squared errors in a three dimensional volume are calculated as a
function of the source plane location 20. The pressure sources used for these calcula-
tions are located on a series of evenly spaced planes for 20 ranging from 0 to 3cm with
an interval of 0.15cm. The pressure in each source plane is calculated with the fast
nearfield method using 20 abscissas for each integral. In Figure 5.4(a), the pressure
in the source plane is truncated by a 7.8cm x 7.8cm square window, which is slightly
larger than the array aperture (7.31cm x 7.31cm). The resulting errors decrease for a
short distance, then oscillate between 0 and 0.02 as 20 increases. In Figure 5.4(b), the
pressure in the source plane is truncated by a 16.2cm x 16.2cm square window, and
the errors monotonically decrease as 20 increases. The error in Figure 5.4(b) drops
sharply when zo increases from 0 to 0.15cm (one wavelength at 1MHz), and then the

error remains small for 20 2 A. Figure 5.4(a) shows that, when a small truncation

81

9
01

9
4;

Normalized spectrum

9
co

9
N

9
J

O

 

 

 

 

 

 

:9}: 1
' .- _ Iii!
r .ffilfi':

 

 

 

 

0 1 2
kx (normalized by k)

(a) Angular spectrum for the particle velocity source.

0.5

9
A

.o .o
N 00

9
_s

 

 

 

 

 

Nonnallzed spectrum

 

 

 

 

0 1
k)( (normalized by k)
(b) Angular spectrum for the pressure source.

Figure 5.3. The angular spectra of the particle velocity source and the pressure source for
the 1024 element planar array in Figure 4.2 along the kx direction for kg = 0. The angular
spectra are obtained by applying a 2D FFT to the particle velocity or pressure source at
z = 0, where the source is sampled with 6 = A/2 (solid line) and A/4 (solid-dot line).

82

window is used, the errors oscillate as the location of the source plane changes, but
the error is approximately constant for 20 Z 0.15cm with a much larger source plane
window. In Figure 5.4(a), the minimum error occurs at 20 = 0.450m; however, the er-
ror is also small at 20 = 0.15cm. In Figure 5.4(b), the error is very small for all values
of 20 Z A. Figure 5.4 shows that the source plane should be at least one wavelength
from the array aperture to avoid sampling problems with evanescent waves near the
array aperture. Figure 5.4(a) also suggests that, when the source plane window is
slightly larger than the array aperture, the resulting error is sufficiently small for 20
equal to 0.150m (one wavelength). Thus, for either source plane size, 20 = A (0.15cm)
is an appropriate source plane location for the 1024 element array.

To demonstrate the influence of the size of the window that truncates the source
plane on the computed pressure, the pressure source plane for the 1024 element phased
array is located at 20 = 0.15cm, and the angular spectrum calculations are evaluated
for a 6cm X 6cm source plane and a 7.8cm X 7.8cm source plane. The resulting axial
pressures are shown in Figure 5.5. The source plane extent specified by L = 6cm is
smaller than the array, so the resulting pressure field (dash-dot line) deviates from the
reference by a significant amount, especially in the region around the focus. When
the extent of the source plane is specified by L = 7.80m, the axial pressure (dashed
line) closely matches the reference pressure (solid line) in Figure 5.5.

Figure 5.6 shows the root mean squared errors evaluated in a three dimensional
volume, where the pressure source planes are truncated by square windows with sizes
ranging from 6cm to 20.4cm. For all of the results shown in Figure 5.6, the source
plane is located at 2:0 = 0.15cm. The two markers in Figure 5.6 indicate the values
of L evaluated in Figure 5.5, where the circle denotes L = 6cm and the star denotes
L = 7.8cm. In Figure 5.6, the errors monotonically decrease as L increases. The
errors are larger when the window size is smaller than the size of the array aperture

(7.31cm), and the errors are smaller when the window size is larger than the array

83

 

 

 

 

 

0.025 . - . . -
‘5
50.02
'0
93
30.015
0'
U)
E 0.01~
Q)
E
30.0051
0:
0 I I I I L
0 0.51 1.5 2 2.5 3

zo(cm)

(a) Root mean square errors obtained when the pressure
source plane is truncated by a 7.8cm X 7 .8cm window.

0.025

 

Pg?
338

9
o
o
or

Root mean squared error

 

 

 

 

%

0 0511152 25 3
zo(cm)

 

(b) Root mean square errors obtained when the pressure
source plane is truncated by a 16.2cm X 16.2cm window.

Figure 5.4. Root mean squared errors evaluated in 3D volumes as a function of the
pressure source plane location 20. The errors are obtained from pressure fields generated by
the 1024 element planar phased array, which is electronically focused at (0,0,10)cm. The
size of the source plane is (a) 7.80m X 7.8cm and (b) 16.2cm X 16.2cm, which is sampled at
Al/ 2 and zero-padded to a 512 X 512 grid. The location of the source plane ranges between
20 = 0 and 20 = 3cm with an increment of 0.15cm.

84

 

.3
U1

 

 

 

 

 

 

 

1 —R;ference I
''''' Source plane L=60m ]
g 1-2' ---Source plane L=7.80m
g
a 0.9
'U
a.)
'7; 0.6
E
o
2 0.3
0 4

 

Figure 5.5. Axial pressures simulated with the angular spectrum approach using pres-
sure source planes that are truncated by square windows of different sizes. The pressures
are generated by the 1024 element planar phased array, which is electronically focused at
(0,0,10)cm. The reference pressure is indicated by the solid line, the pressure computed
with the angular spectrum approach using a 6cm X 6cm pressure source plane is represented
by the dash-dot line, and the axial pressure computed with the angular spectrum approach
using a 7.8cm X 7 .8cm pressure source plane is represented by the dashed line. The solid
line and the dashed line are nearly coincident, which indicates that L = 7.8cm is sufficiently
large for a 7.310m X 7.310m phased array.

aperture. Figure 5.6 suggests that only a moderate reduction in the error is achieved
for values of L > 7.8cm. Furthermore, in Figure 5.5, the results for L = 7.80m are
nearly coincident with the reference. Thus, L = 7.8cm is identified as the optimal size
of the source plane for the 1024 element array. If L < 7.8cm is used, the truncation
of the pressure wavefront causes an increase in the root mean square errors. A source
plane window with L > 7.80m consistently produces small errors. However, there is
a trade-off between the accuracy and the efficiency of these calculations for larger L.
On the one hand, accurate results are achieved when large values of L are used, and
larger source plane sizes are often necessary for array simulations that also include
grating lobes, as in Figure 4.11. On the other hand, L should be as small as possible

because the computation time and computer memory required are proportional to

85

the number of grid points, which is determined by the value of L when (5 is fixed.
Figure 5.4 indicates that the minimum error is achieved at 20 = 0.45cm, and Figure
5.6 suggests that L = 7.8cm is optimal for the 1024 element planar array. The results
in Figure 5.6 also show that the pressure source plane located at 20 = 0.15cm achieves

acceptably small errors.

 

0.12

Root mean squared error
0
O
O)

 

 

 

 

 

 

 

0.04
0.02- 1
0 I
6 10 14 18 22
L(cm)

Figure 5.6. Root mean squared errors plotted as a function of the source plane extent L.
The errors are obtained from pressure fields generated by the 1024 element planar phased
array, which is electronically focused at (0, 0,10)cm. The pressure source plane for angular
spectrum simulations is located at zo = 0.15cm and truncated by L X L square windows,
where L ranges from 6cm to 20.4cm with an increment of 0.15cm. The two markers indicate
the values of L shown in Figure 5.5.

5.4 Evaluation of input and output errors

In angular spectrum calculations that use the spectral propagator Hp(k$, ky,Az),
the source pressure in an initial plane is typically simulated with analytical integral
approaches, and this result provides the input to angular spectrum simulations that
evaluate the pressure in a three dimensional volume. For these simulations, fast and

accurate calculations of the source pressure are desirable. This motivates numerical

86

evaluations of the input error associated with the pressure source plane that determine
the impact of the input error on the computed 3D pressure output.

Using the spatial impulse response method as the reference, two analytical in-
tegral methods are compared for simulations of the source pressure: the Rayleigh-
Sommerfeld integral and the fast nearfield method. The accuracy of these methods is
determined by the number of abscissas used in evaluations of the integrals in equation
2.1 and equation 2.11, respectively. As demonstrated in Ref [24], the fast nearfield
method achieves higher accuracy with fewer abscissas relative to other single integral
approaches in the simulations of single transducers. This result also holds for phased
array simulations. To demonstrate the change in the root mean squared error in the
source plane as the number of abscissas increases, the source pressure is computed in
a 7.8cm X 7.8cm source plane at a depth of zo = 0.15cm, where 7.8cm is the optimal
value of L determined in the previous section for the 1024 element phased array. With
a spatial sampling interval of 0.075cm (A / 2), the computational plane is discretized
to 105 points in both the m and the y directions. The root mean squared errors
are evaluated for the source plane pressure computed with the Rayleigh-Sommerfeld
integral using 2 X 2 to 10 X 10 abscissas and with the fast nearfield method using
2 to 10 abscissas. The results are plotted in Figure 5.7, where the errors from both
methods decrease as the number of abscissas increases. The root mean squared error
obtained with the Rayleigh-Sommerfeld integral approach in the source plane is 0.216
when 2 X 2 abscissas are used, and the error in the source plane decreases to 0.006
for 10 X 10 abscissas. In contrast, the root mean squared error obtained with the fast
nearfield method in the source plane is 0.076 for 2 abscissas, and the error quickly
decreases to 0.0004 when only 4 abscissas are used. In Figure 5.7, errors less than or
equal to 0.0004 are coincident with the horizontal axis when the values in the range
shown are plotted on a linear scale.

Figure 5.8 demonstrates the influence of the number of abscissas used for source

87

 

 

 

 

 

0.25 . .
--Fast nearfield method
*- --Ra lei h-Sommerfeld
g 0.2 y g -
d)
3 .5
(U 0.
é a
E 2
‘3' .5
0:

 

 

 

 

Abscissas

Figure 5.7. Root mean squared errors in the source plane plotted as a function of the
number of source plane abscissas used in each direction. The errors are obtained from
pressure fields generated by the 1024 element planar phased array, which is electronically
focused at (0,0,10)cm. The source plane is located at 20 = 0.15cm and truncated by a
7.80m X 7.80m window. The source pressure is computed with the Rayleigh-Sommerfeld
integral using 2 X 2 to 10 X 10 abscissas and with the fast nearfield method using 2 to 10
abscissas.

88

0.015 a

 

 

--Fasf nearfield method
-- Rayleigh-Sommerfeld

 

 

 

 

 

 

 

 

5 *5
5.90.012

3
'o o
e e
g @0009
3 e
C Q.
8 0 0.006
E 0’)
4.. d) *
8 f 0.003-
D: o

0 2 4 5 a 10

Abscissas

Figure 5.8. Root mean squared errors for 3D pressure field outputs plotted as a function
of the number of source plane abscissas used in each direction. The errors are obtained
from pressure fields generated by the 1024 element planar phased array, which is electroni-
cally focused at (0, 0, 10)cm. The source plane pressures are calculated with the Rayleigh-
Sommerfeld integral using 2 X 2 to 10 X 10 abscissas and with the fast nearfield method
using 2 to 10 abscissas. The errors obtained from both methods approach the same limiting
value, but the fast nearfield method achieves convergence with far fewer abscissas.

plane calculations on the three dimensional angular spectrum results. In Figure
5.8, the horizontal axis contains the number of abscissas used in each direction for
source plane calculations with the Rayleigh-Sommerfeld integral and the fast nearfield
method, and the vertical axis contains the output root mean squared errors in the
three dimensional pressure field computed with the angular spectrum approach. The
output root mean squared error approaches a limiting value of 0.004 when the fast
nearfield method with 3 or more abscissas calculates the source pressure. When the
source pressure is calculated with the Rayleigh-Sommerfeld integral, the output error
asymptotically approaches the same value. The saturation of the root mean squared
error in Figure 5.8 indicates that there is a lower limit for the output error in angular
spectrum calculations that is determined by the grid size, the grid spacing, and the lo-

cation of the source plane. Moreover, Figure 5.8 shows that the fast nearfield method

89

requires far fewer abscissas to achieve the minimum output error. Therefore, the fast

nearfield method is a much more efficient approach for source plane calculations.

5.5 Computation time and numerical error with the angular

spectrum approach

The angular spectrum approach computes pressures in parallel planes by propagating
the fields in the spatial frequency domain. This method reduces the computation
time significantly compared to conventional integral approaches, which compute the
pressure at individual field points. Table 5.1 shows the computation time (minutes)
and the root mean squared error for the pressures generated by the planar array, the
spherical-section array and the spherically focused array. The pressures are calculated
with the angular spectrum approach (ASA) using the particle velocity source (U0)
and the pressure source (P0).

For the results shown in Table 5.1, the source plane for the planar array is located
at z = 40m, and the source planes for the spherical-section array and the spherically
focused array are both located at z = 60m. The source pressures are truncated by
a 16.2cm X 16.2cm square window in all three cases, which yields a 217 X 217 point
grid with 6 = 0.075cm = A/2. This grid is zero-padded to 512 X 512 before the
2D FFT is applied to the source field. The source pressures are computed by the
fast nearfield method with 2 abscissas for the planar array, with 2 abscissas for the
spherical-section array, and with 3 abscissas for the spherically focused array. Angular
spectrum calculations with the pressure source requires extra time to compute the
source. Table 5.1(a) shows that the angular spectrum approach with the velocity
source computes the 3D pressure field in less time but with lower numerical accuracy.
Due to the concave geometry of the spherical-section array and the spherically focused

array, only the planar array is computed with the normal particle velocity source. A

90

Table 5.1. The computation time and the root mean squared error (RMSE) for the
simulated pressure generated by (a) a planar array with 1024 square element, (b) a spherical—
section array with 1444 square elements, and (c) a spherically focused array with 1441
circular elements. The pressures are simulated with the angular spectrum approach using
the particle velocity source (ASA with U0) and the pressure source (ASA with P0). The
simulated fields are compared with the reference pressure computed with the spatial impulse
response method.

 

Method | ASA with U0 | ASA with P0
Time (min) | 2.1 I 4.35

RMSE | 0.0435 | 0.0028

(a) Planar array.

Method | ASA with U0 | ASA with P0
Time (min) | - ] 3.78

RMSE | - | 0.0068

(b) Spherical-section array.

Method | ASA with U0 | ASA with P0
Time (min) | - | 8.56

RMSE | - | 0.00041

(c) Spherically focused array.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Comparison of Table 4.1 and Table 5.1 shows that the total time for the angular
spectrum calculation and the source pressure calculation is significantly lower than
the time for each of the integral approaches, yet the numerical accuracies of these
methods are comparable. Thus, the angular spectrum approach is extremely efficient

for simulations of large phased arrays.

91

CHAPTER 6

Angular Spectrum Approach for
Ultrasound Phased Array Simulations in

Layered and Nonlinear Media

6. 1 Introduction

Biological tissues are conventionally modeled as homogeneous and linear elastic me—
dia in hyperthermia simulations. Homogeneous media have uniform characteristics
as impedance, speed of sound, and attenuation. However, living tissues have more
complex structures that lead to heterogeneity. When ultrasound propagates in het-
erogeneous tissues, the ultrasound waves reflect on tissue boundaries, scatter, and
diffract due to interactions with tissue structures of different scales. Furthermore,
the speed of sound and attenuation vary in different tissues. The heterogeneity of
tissues sometimes causes deformation in the propagating wavefront, which is called
aberration in ultrasound imaging. The pulsed wave propagation in heterogeneous tis-
sue media can be simulated with a finite difference time domain (FDTD) approach or
a full-wave k-space method [57]. Ultrasound propagations in abdominal wall [58, 59],
breast tissue [60, 61] and chest wall [62] have been simulated while incorporating
different properties of these anatomical structures. However, due to the memory and

time intensive nature of the FDTD and the k-space algorithms, these simulations

92

often employ two dimensional models. Large scale simulations in three dimensional
volumes can be realized with the ray theory or the angular spectrum approach, but
the inhomogeneous tissues are usually modeled as planar or curved layers in these
simulations. For example, the wave reflection and refraction at soft tissue boundaries
was simulated by Fan and Hynynen [63, 64] using a ray theory based method. The
calculation of the reflection of ultrasound waves on planar interfaces using the angu—
lar spectrum approach was proposed by Christoper and Parker [42] and extended by
Clement and Hynynen [65] to model the propagation of ultrasound waves in nonpar-
allel layers.

Nonlinearity can also appear as ultrasound waves propagate in tissue. Nonlinear
wave distortion is a progressive process that depends on the input power intensity,
the ultrasound frequency, the nonlinearity of the media, and the distance the wave
travels. The linear theory breaks down when the role of nonlinearity increases, and
the linear wave equation no longer describes ultrasound wave propagation. The non-
linear wave equation is instead adopted for more realistic simulations. Efforts have
focused on solutions of the Burger’s equation for plane waves [66, 67], the generalized
Burger’s equation for plane, cylindrical, and spherical waves [68], and the Khokhlov—
Zabolotshaya-Kuznetsov (KZK) equation [69] for general diffracted waves using a
parabolic approximation. The numerical solutions to these equations either use a
temporal frequency domain solution [69, 70, 71] that utilizes Fourier series expan-
sions, or use a time domain approach [72, 73, 74] for pulsed wave simulations. The
modeling of nonlinearity is incorporated into the angular spectrum approach mainly
through the contributions from Christopher [75, 76] and Vecchio [77, 78].

In biological tissues, harmonic waves with higher frequencies are absorbed more
rapidly than those with lower frequencies. Therefore, the power deposition and tem-
perature distribution induced by nonlinear propagation are different from those en-

countered in linear wave propagation. The increased heating caused by nonlinear

93

waves has been investigated in fluid media with low attenuation, e. g., water [79], agar
[80], urine and amniotic fluid [81], and in tissues with higher attenuation such as hu-
man calf muscle [82] and dog thigh [7]. Shock waves are more likely to form in fluid
media with low absorption, which produces intensive surface heating when the wave
enters a medium with higher absorption coefficient. In tissues with higher attenua-
tion, shock waves are rarely generated due to the rapid absorption of harmonics. The
investigation of the media with smaller attenuation values provides some insight into
the design of a water bolus, however, this chapter focuses on nonlinear modeling in
tissue where the harmonics induce additional heating.

This chapter discusses the influences of inhomogeneity and nonlinearity on ul—
trasound hyperthermia simulations. The angular spectrum approach for modeling
ultrasound wave propagation in parallel layers and in nonlinear media is reviewed
in sections 6.2.1 and 6.2.2, respectively. The bio-heat transfer equation that models
the local temperature variation in tissue is then introduced. The power and tem-
perature distributions in a layered tissue medium that mimics the human abdomen
are simulated in section 6.3.3. The power deposition produced by focused ultrasound
in the layered model is very similar to that in the homogeneous model, where the
homogeneous model uses the weighted average of different acoustic properties in the
layered model. Nonlinear propagation of ultrasound and the resulting temperature
in fatty tissue is simulated in section 6.3.4. The temperature elevation produced by
nonlinear harmonics is shown to be negligible in mild temperature hyperthermia, but

nonlinear contributions are important for high intensity ultrasound ablation.

94

6.2 Acoustic simulations in layered and nonlinear media

6.2.1 Angular spectrum approach in layered media

If the acoustic properties of a medium, e.g., density and speed of sound, only vary
a small amount as a function of distance evaluated on a wavelength scale, the wave
propagation in the medium can be approximately modeled with ray theory. Ac-
cording to linear ray theory, as an acoustic wave encounters the interface between
two fluids, part of the wave is reflected back into the first medium, and the rest is
transmitted into the second medium. The refracted wave potentially differs from
the incident wave in both phase and amplitude. However, the complex pressure and
the normal acoustic velocity are continuous across the interface. The reflection and
transmission coefficients at a plane wave on a planar interface are well-known [18].
The simulation of acoustic wave propagation through layered media by the angular
spectrum approach has also been established [42, 65]. In the spatial frequency do-
main, each of the plane wave components is associated with an angle of incidence,

and the transmission coefficient for each plane wave is

 

 

plcl cosflt
T k k = - = __ _ .
p( :c, y) Pt/Pz 2/(1+ p2C2 cosei)’ (61)
c c050-
mkx. Icy) = zit/14: 2/(1+ 9—2 z). (6.2)
plcl cosélt

where 6; and 9t are the incident and refraction angles, 1); and pt are the pressures
produced by the incident and refracted waves, u; and at are the corresponding particle
velocities, Tp and Tu are the transmission coefficients for the pressure and the particle
velocity, respectively, and p101 and p202 are the acoustic impedances in medium 1
and medium 2, respectively. Multiple reflections between boundaries are ignored in
Equation 6.2 [42]. Snell’s law holds for each plane wave component in the spatial

frequency domain; therefore, the incident angle and the refraction angle in the plane

95

of interface can be expressed as

 

 

 

 

\/k2 — kg — k5
cos 6.; = k . (6.3)
2 2
6 = — — . .4
cos t \/ (Cl) 1:2 (6 )

The angular spectra of the pressures at the interface between two different media
are related by
132(ka kya 3) = P1(km, ky, ZlTp(ka (9y), (6-5)

where P1 and P2 are the angular spectra in medium 1 and medium 2, respectively.

Similarly, the particle velocities at the interface between two media are related by

 

 

 

 

 

 

 

 

0.4 . .
—lnterface at z=60m
----- Interface at z=7cm
A 03 - --lnterface at z=8cm.
§ ------- Interface at z=90m
12',
e 0 2 K’I—‘Hsm. .,
5’, .
U)
9
‘1 0.1
O 0

 

Figure 6.1. The axial pressures generated by a spherically focused transducer in a medium
with water and fatty tissue layers. The planar interface between the water and fatty tissue
layers is located at 6cm, 7cm, 8cm, and 9cm. The pressure is calculated with the angular
spectrum approach.

Figure 6.1 shows the axial pressure generated by a spherically focused transducer

with a 2cm aperture diameter and a 10cm radius of curvature. The bottom of the

96

transducer is located at the origin, and the normal at the center of the transducer
coincides with the z axis. The pressure generated by this transducer propagates in
a medium composed of one layer of water and one layer of fatty tissue. The planar
interface between two layers is perpendicular to the normal of the transducer, and
the location of the interface is equal to z = 6cm, 7cm, 8cm, and 9cm. The intensity
on the transducer surface is 0.03W/cm2. Figure 6.1 shows that the fatty tissue has a
higher attenuation than the water layer. Figure 6.1 also indicates that the pressure

is continuous at the planar interface for all results shown.

6.2.2 Angular spectrum approach in nonlinear media

The relationship between the pressure and the density in fluid media can be described

by a Fourier series [18]

119-290: 00(g%)p=po:m:: +8233:ng P0(p _Bfl)2+"' (6.7)

where p is the total pressure, p0 and p0 are the pressure and the density at equi-
librium. Equation 6.7 is valid under adiabatic conditions, that is, the heat transfer
during density fluctuations is small, and the entropy of the medium remains nearly
constant. In a linear medium, the pressure and the density fluctuations follow a linear
relationship. Therefore, only the first term of Equation 6.7 is retained, with A = pocg,
where CO is the speed of sound at equilibrium. In a nonlinear medium, the higher
order terms in Equation 6.7 contribute to the right hand side. The third order terms
and above are generally ignored. The nonlinearity of a medium is characterized by

the parameter of nonlinearity

[3 = 1 + B/2A, (6.8)

where B and A are defined in Equation 6.7.
I As a sinusoidal plane wave propagates through a nonlinear lossless medium, the

particle velocity/ pressure waveform steepens until a sawtooth-like waveform (or N-

97

wave) is formed. The progressive distortion of the waveform corresponds to the gen-
eration of harmonics that are multiples of the fundamental frequency. The particle
velocity of the acoustic wave is usually much smaller than the speed of sound in fluid
media. Therefore, this distortion is called finite amplitude distortion, as opposed to
the small amplitude assumption used in linear theory. Convection and the nonlin-
earity of the medium are the two main reasons for finite amplitude distortion [27].
When the slope of the waveform becomes large, a shock wave is formed, and energy is
dissipated at the shock wavefront. The lossless model fails after the shock is formed.
The shock wave keeps propagating with a decaying amplitude. Eventually, the wave
returns to the original shape with a much lower amplitude, and this amplitude is
independent of the initial amplitude of the wave.

Nonlinear simulations can be incorporated into the angular spectrum approach [44,
75, 76, 78, 83] through the solution of Burger’s equation. Burger’s equation describes
the propagation and dissipation of a finite amplitude plane wave in themoviscous

lossy media. Burger’s equation in terms of pressure is

8;» [3 ap_ 6821)

a “ zaps: — 2729?? (69)

where T = t — z/c is retarded time, 6 is the parameter of nonlinearity, and 6 is the
diffusivity of sound. Under the weak attenuation condition 0 << k, the relationship

between the diffusivity 5 and the absorption coefficient a is
6w2/2c3 a: a, (6.10)

where w = 27r f is the radian frequency of the monochromatic excitation. By substi-

tuting Equation 6.10 into Equation 6.9, the Burger’s equation for pressure is rewritten,

0p 6 8p_ aazp

a; ‘ gape—T - @3577 (6“)

Through a Fourier series expansion, the pressure can be expressed as the summation

98

of an infinite number of harmonic waves,
1): pa + Z pn cos(nwr + no): :2 PnejWT, (6.12)
n = 12—n — —00
where Pn = pnemff, and Pn = Pin. By substituting Equation 6.12 into Equation

6.11, the following expression is derived for each harmonic pn,

3pm . fiw
.87 = 3%?(16 n2: kPkPn—k +k 23+ lnpkpk_ n) —om2 pn. (6.13)

The left hand side of Equation 6.13 can be written as a forward difference,

_(ZPE : pn(z + A2) — 1971(2)
82 Az °

 

(6.14)

Equation 6.13 and Equation 6.14 are then incorporated after each linear angular
spectrum step to model nonlinearity. This is done by taking a virtual step that

projects the linear field onto the nonlinear field with [75],

A
pn(z+Az)= p’n(z+Az)+ jZfiw pc3:;z(

(6.15)
where all terms with 13’ represent the pressures obtained from the linear angular spec-
trum calculation, and pn is the nth harmonic obtained from nonlinear propagation.
The term p’n(z+ A2) is computed with the spectral propagator Hp using p.’n(z) as the
source. The spatial propagator is not appropriate for this purpose, because the prop—
agation distance A2 is on the order of a wavelength at the fundamental frequency,
and the spatial propagator approaches a delta function when A2 approaches zero,
which causes significant sampling errors. The initial value of pn with n > 1 is zero
and is progressively filled with non-zero values during propagation. The value of pn
is computed at each point in each successive plane. The second series on the right
hand side of Equation 6.15 is usually truncated at N’, where N’ is the maximum
number of harmonics included in the simulation. Because of this truncation, the en—

ergy flow from lower harmonics to higher harmonics stops at order N’, and the higher

99

harmonics may reach abnormally large amplitudes and induce a numerical overflow.
An ad-hoc method that automatically updates N’ compares the maximum amplitude
of the present highest harmonic to the maximum amplitude of the fundamental har-
monic, and if the former exceeds 1% of the latter, N’ is increased by 1 [78]. Using
this scheme, the simulation can start with a relatively small N’, and N’ is increased
during propagation. Diffraction is not modeled by Equation 6.15, because Burger’s
equation only models plane wave propagation. A more accurate model is provided
by the KZK equation [69]. However, the diffraction term in the KZK equation cal-
culates the Laplacian of the pressure on the transverse coordinate. The nonlinear
dependency on the transverse grid is incompatible with angular spectrum approach,
which utilizes linear plane wave superposition.

In order to validate the nonlinear angular spectrum approach, nonlinear harmonics
generated by a plane circular transducer in water are simulated, and the simulation
results are compared with the experimental results in Ref [84]. The transducer is
a planar disc with 1.90m radius. The operating frequency (fundamental frequency)
is 2.25MHz, the speed of sound is 1486m/s, the parameter of nonlinearity fl is 3.5,
and the attenuation coefl‘icient is 2.5 x 10‘4Neper/cm/MH22. The peak pressure
on the transducer surface is 100kPa. For angular spectrum simulations, the spatial
sampling interval is 0.0165cm (A / 4), the source plane size is 15.85cm x15.85cm, which
corresponds to 961 x 961 grid points in one transverse plane, and this grid is zero-
padded to 1600 x 1600 for FFT calculations. A total of 10 harmonics are included
in the simulation. The results from the angular spectrum simulations are shown in
Figure 6.2 and Figure 6.3. Figure 6.2 shows the absolute axial pressures of the first
three harmonics along the z direction for a: = y = O. The fundamental, the second
harmonic, and the third harmonic have descending peak amplitudes in Figure 6.2,
and the amplitudes of the second and third harmonics grow as 2 increases. Except

for the ripples in the fundamental harmonic field due to undersampling in the spatial

100

 

 

Pressure (MPa)
0

 

 

 

 

 

o 20 40 6b 80
z(cm)

Figure 6.2. Simulated fundamental, second, and third harmonics along the line at a: =
y = 0 generated by a plane circular transducer in water. The transducer is excited by a
2.25MHz continuous wave input with a peak amplitude of 100kPa. The solid line with the
largest peak amplitude represents the fundamental harmonic, the lines with the second and
third peaks represent the second and the third harmonics.

domain, Figure 6.2 closely resembles Figure 1 in Ref [84]. The first three simulated
harmonics along the a: direction for y = 0 at z = 27.5cm and z = 50cm are shown in
Figures 6.3(a) and 6.3(b), respectively. These two radial pressure plots closely match
Figure 5 in Ref [84].

Figure 6.4 shows two snapshots of the time harmonic nonlinear pressure waveforms
generated by a spherically focused transducer with a 10cm aperture diameter and a
20cm radius of curvature. The excitation signal is a sinusoidal wave with f = 1.5MHz.
The pressure is evaluated at the geometrical focus of the transducer. The pressure
waveform is given by

p(z, t) = Z pcun(z)ej2’mf(t — z/C). (6.16)

17.

Figure 6.4(a) shows two cycles of the output pressure waveform in fatty tissue,
where a = 0.07Neper/cm/MHz, b = 1, and fl = 10.5. Due to the relatively high
absorption of fatty tissue, only 5 harmonics are required in this calculation. Figure
6.4(b) is the output waveform in water, where a = 2.5 x 10'4Neper/cm/MH22, b =
2, and fl = 3.5. A total of 50 harmonics are used for the calculation in water.

Figure 6.4(a) and Figure 6.4(b) are calculated with the same amount of input power;

101

 

 

 

 

 

 

 

 

 

0
’5
‘5
"" -20' ‘
9
CD
3
92
a -40
w
9
D.

"3% -25 0‘ 2.55

x(cm)
(a) Simulated radial pressures at z = 27.5cm.
0

75?
‘5'
‘- -20"
9
CD
3
9
c?) -40'
m
9
0..

-§% -2.5 o' 2.5 5

x(cm)

(b) Simulated radial pressures at z = 50cm.

Figure 6.3. Simulated fundamental, second, and third harmonics generated by a plane
circular transducer in water. The transducer is excited by a 2.25MHz continuous wave
input with a peak amplitude of 100kPa. The pressures are plotted along the :5 direction for
y = 0 at (a) z = 27.50m and (b) 2 = 50cm. The solid line with the largest peak amplitude
represents the fundamental harmonic, the lines with the second and third peaks represent
the second and the third harmonics.

102

however, the waveform in water has much higher output amplitude and stronger finite
amplitude distortion than that in fatty tissue. Note that unlike the initial waveform,
neither waveform in Figure 6.4 is symmetric with respect to the horizontal axis. This
is a result of superposing harmonics with different phases. Also, because of the phase
differences, each cycle of the waveform in Figure 6.4(a) has two steep peaks and one
slow-varying trough instead of one sawtooth-like peak and trough. The result in
Figure 6.4(b) agrees with the experimental result in Figure 3(b) of Ref [7].

The pressure generated by a spherically focused transducer with a 20m diameter
and a 10cm radius of curvature is simulated using a higher input intensity so that
a nonlinear waveform is produced. In Figure 6.1, the transducer is excited with an
input intensity of 0.03W/cm2, and linear wave propagation is assumed. When the
input intensity is raised to 3W/cm2, the nonlinear harmonics have a much greater
influence on the pressure field. Figure 6.5 shows the first three harmonics generated
by the spherically focused transducer with a 3W/cm2 input intensity. The simulation
model includes one water layer and one fatty tissue layer, and the interface between

the two layers is located at z = 9cm.

6.3 Power and temperature simulations in layered and non-

linear media

6.3.1 Acoustic intensity and power deposition

As acoustic waves propagate through fluid media, energy is transported through the
kinetic movement of particles and the compression of media. For a monochromatic
wave, the acoustic intensity (W/m2) is defined as the time averaged rate of energy

transmission through a unit area [18],

T
I = 511/0 p(r,t)u(r,t)dt, (6.17)

103

0.6

 

Pressure(MPa)

 

 

 

 

"0'40 025 1:0 1.5
Time(x1 0'6)
(a) Nonlinear waveform in fatty tissue.
15 - a

 

Pressure (MPa)

 

 

 

 

0 0:5 1T0 1.5
Time (x 10'6 s)

(b) Nonlinear waveform in water.

Figure 6.4. Nonlinear pressure waveforms generated by a spherically focused transducer

with a 10cm aperture diameter and a 200m radius of curvature at the focus (a) in fatty
tissue and (b) in water.

104

 

Pressure (MPa)
O .0 .0
h g) (I)

.0
N

 

 

 

 

 

0 i 5 5 9 12 15
z(cm)

Figure 6.5. Simulated fundamental, second, and third harmonics along the line at a: =
y = 0 generated by a spherically focused transducer with a 3W/cm2 input intensity. The
simulation model includes a water layer and a fatty tissue layer, and the interface is at z =
9cm. The solid line with the largast peak amplitude represents the fundamental harmonic,
the lines with the second and third peaks represent the second and the third harmonics.

where p and u are the instantaneous pressure and particle velocity, respectively, r is
the position vector, and T is one period of the monochromatic wave. Both 19 and u are
functions of time and space. For a plane harmonic wave, Equation 6.17 is equivalent
to

I = P2/(2pc), (6.18)

where P is the amplitude of the plane wave, p is the density, and c is the speed of
sound.
Similarly, the power density (W/m3) is defined as the time averaged rate of power

transmission in a unit volume,

1

T
Q = T [0 p(r,t)V-u(r,t)dt. (6-19)

In an attenuating medium, ultrasound power is dissipated during propagation and

converted into heat. For a plane wave, the power deposition is proportional to the

105

absorption coefficient and the square of the pressure amplitude [85],

Qp(rc,y, 2) = imam), (6.20)

where oz is the absorption coefficient, and "' denotes the complex conjugate of the
pressure. Equation 6.20 is derived from Equation 6.19 by using the equation of state
and the equation of continuity, and plugging in the complex wavenumber k = w / c— j a

[85]. Equation 6.20 also assumes homogeneous attenuation with zero shear velocity.

6.3.2 Bio-heat transfer equation

Heat transfer in tissue is mainly achieved via conduction and convection. Conduction
is the exchange of molecular kinetic energy due to collisions among particles, while
convection is the heat transfer due to fluid movement, i.e., the circulation of blood [1].
Currently, most analytical thermal modeling methods assume tissue is a continuum
that only consist of micro-circulatory vessels, and the heat convection in large vessels
is treated separately. The Pennes bio-heat transfer equation (BHTE) [86] is one
model for local heat transfer in biological tissue. In the bio-heat transfer equation,

the temperature field as a function of time and space is modeled by [87]

BT
par—5t— = KVZT — Wbe(T — Ta) + Qp, (6.21)

where p is the tissue density (kg/m3), K is the thermal conductivity of tissue
(W/m/OC), C and Cb are the specific heats (J /kg/°C) of tissue and blood, respec-
tively, Wb is the blood perfusion rate (kg/m3/s), Qp is the power deposition in a
unit volume (W/m3), T = T(:r, y, z, t) is the total temperature, and Ta is the arterial
blood temperature, which is 37°C. Equation 6.21 ignores the power deposition due to
metabolism, which is usually much smaller than the heat produced by ultrasound de-
vices. K is represented by a constant in Equation 6.21, and the term K VZT accounts
for heat dissipation by conduction. The term WbeT models the convection losses

due to blood circulation, where non-directional blood flow is assumed [3]. Equation

106

6.21 is the transient BHTE, where 8T/Bt denotes the time derivative of the temper-
ature. When a steady state temperature is achieved, the temperature is modeled by

the steady state BHTE [88]
KV2T — Wbe(T — Ta) + Qp = 0. (6.22)

Equation 6.22 is usually solved with a the finite difference (FD) approach. By
discretizing the computational volume in a rectilinear grid and using the central
difference to approximate the second order derivative, Equation 6.22 in a discretized

grid is expressed as

 

',',k_ K 52 mph ’+1,',k '-1,',k
TH _ 6K+Wbe62(KQP +T’ J “LT: J

6.23
+Ti,j+l,k+Ti,j—1,k+Ti,j,k+1+Ti,j,k—1), ( )

where (2', j, [6) denotes the indices of the grid points in the :r, y and 2 directions, and
6 is the uniform spacing between adjacent grid points in the x, y and 2 directions.
Equation 6.23 is typically solved with known temperatures at the boundaries of the
computational volume, i.e., in the planes where i = 1, N53, j = 1, Ny, and k = 1, N z.
The boundary conditions for hyperthermia are usually determined at the skin surface
by water bath and internally by the core temperature. The finite difference approach
is a highly localized numerical scheme, where calculations at each point only depend
on the values at the immediate neighbors. The solution of Equation 6.23 converges
after multiple iterations.

Equation 6.21 is implemented with a finite difference solution in the time domain.

Equation 6.21 is discretized in both spatial and temporal domains, and the result is

" AtWC " — AtK ' 1'1: —1
T2,],k,m= 1_6AtK_ 9 l2 T2,],k,m 1+ (Tz+ ,J, ,m
62pC P 62pC
+Ti—1,j,k,m—1+Ti,j+1,k,m— +Tz',j—1,k,m—1+Ti.j,k+1,m-1 (6-24)
+Tl,j,k—l,m—1)+%szijikym—1,
where At is the duration of each time step, and m denotes the index of the time step.

Equation 6.24 is iteratively advanced in both the spatial and temporal domains. To

107

avoid dispersion, At and 6 must satisfy the relationship

62pC

6.3.3 Temperature simulation in a layered tissue model

In a tissue model with parallel layers, wave propagation is simulated by the angular
spectrum approach where reflection and refraction are modeled in the spatial fre-
quency domain. A layered tissue model that mimics the human abdomen [16] is
shown in Figure 6.6. The water layer (2.8cm in depth) represnts the water bolus that
is inserted between the patient and the ultrasound applicator for coupling and cooling
purposes. The tissue layers consist of 0.4cm skin, 2cm fatty tissue, 3cm muscle, and
13.6cm viscera. The thermal and acoustic properties for these materials are listed in

Table 6.1.

Table 6.1. The acoustic and thermal prOperties of water and some biological tissues

 

Tissue I mks/ma) I C(m/S) I 0z(l‘lP/Cm/NHLIZ) I 3 I K(W/m/°C) I C(J/k8/°C)

 

 

 

 

Water | 1000 | 1500 | 2.5x10-4 | 3.5 | 0.615 | 4180
Fat | 921 | 1445 | 0.07 |10.5| 0.223 | 2325
Muscle] 1138 | 1569 | 0.09 (4.75) 0.4975 | 3720
Skin | 1200 | 1498 | 0.14 | - | 0.266 | 3430
Visceral 1060 | 1540 I 0.032 | - | 0.512 | 3600

 

The pressure generated by a spherically focused phased array is simulated in a
homogeneous medium and a layered tissue medium in Figure 6.6. The array contains
257 5 circular elements, where the diameter of each element is 0.24cm and the center-
to-center spacing is 0.2475cm. The normal evaluated at the center of each element
points towards the geometrical focus of the array. The diameter of the aperture is
16cm, and the radius of curvature is 16cm. The array is driven by a 1MHz continu-
ous wave excitation. The pressure field is first simulated in a homogeneous medium
by the fast nearfield method using Riemann’s mid—point rule evaluated with 20 ab-

scissas. The attenuation coefficient and the density of the homogeneous medium are

108

, i Irllltiateri(vary) ‘
Skin (0.4cm)
Fat (2 cm)

Muscle (3 cm)

Viscer‘a (13.6 cm)

 

Figure 6.6. A model of the human abdomen that includes four soft tissue layers.

109

the weighted averages of the parameters in the layered model, where the weight is
the thickness of each layer divided by the total thickness. The average attenuation
coefficient is 5 = 0.047Np/cm, and the average density is 'p' = 1087kg/m3. The speed
of sound in the homogeneous medium is the average speed of sound in the layered
model, E = 1513m/s. The pressure field in the layered tissue medium is computed
with the linear angular spectrum approach. In the angular spectrum simulations,
the spatial sampling interval is 6 = 0.075cm, the grid size is specified by N = 512,
and the pressure source plane is located at 20 = 2.25cm (within the water layer) and
truncated by a 16.2cm x 16.2cm square window.

Figure 6.7 compares the axial pressures generated by the spherically focused
phased array in the homogeneous model (solid line) and in the layered tissue model
(dash-dot line). The water-skin interface in the layered tissue model is located at
z = 50m. The pressure fields in Figure 6.7 are normalized by the peak amplitude
of the pressure in the homogeneous medium. The pressure amplitude in the layered
medium is very close to that in the homogeneous medium when average parameters
are used in the simulation. The focus of the field in the layered medium is shifted
about 0.3cm towards the array as opposed to that in the homogeneous medium. This
is due to the variation of the speed of sound in the layered medium. As the sound
beam travels from a medium with lower speed of sound (fat) to a medium with higher
speed of sound (muscle), the beam is bent away from the normal axis of the array.
Overall, the difference between the simulated pressure in homogeneous and layered
media is relatively small.

Figure 6.8(a) shows the normalized power deposition in the 3:2 plane at y = 0.
Figure 6.8(b) shows the temperature rise in the same plane with the peak tempera-
ture rise normalized to 6°C. The water layer is excluded from power and temperature
calculations because the water bolus temperature is usually remained at a constant

value with a circulating water bath. The temperature at the boundary of the compu-

110

1 a

 

    

 

 

 

 

 

 

 

 

—Homogeneous model g,
----- Layered tissue model - '3
9 0.8- ' s
3 I
a) !
8 a
50.6 a,
'0 3
f3 3.
=5 0.4- l,
E z
0 i

2 0.2 g,

0 1 ' | \.‘l ‘ ‘-
4 12 16 20 24

Figure 6.7. The normalized axial pressure generated by a 2575 element spherically focused
phased array evaluated along the z axis where a: = y = 0. The pressure is simulated in a

homogeneous fatty tissue medium (solid line) and a model containing layers of water, skin,
fat, muscle and viscera (dash-dot line).

tational grid is 37°C, which corresponds to the arterial blood temperature. Although
the power deposition in the region between the phased array and the focal spot is
relatively small in Figure 6.8(a), the temperature rise in this region is relatively large
due to the large absorption coefficients of fat and muscle. Thus, Figure 6.8 shows an

example of intervening tissue heating, which causes great difliculties in hyperthermia
[89].

111

4".

 

Nonnallzed power deposition

 

x(cm) "8 5 z(cm)

. _.-
- .
6" ‘
.-'

9’

  
  
 
   
    
    
 

l

_.
.'"
....

."
,-

I

.,. 6

_.
,6' .
. .
. ..1
.I'

.....
......
....
......
.....
_,.I

l

NNNNNNN
,,,,,,,,,,
27;;

l

 

Temperature rise
0 —t N E «D- (II

H

.;;; 11111 II/l/llll/I ”ll/I” ”III/1,, .

8 ”””//IIII// III/1, [”11” "If” "is.
4 Ill/Illl///”[[I[I[I” Ill/III .;..6 .6.66
Ill/III III/711]?)
‘4 " ' 1o

x(cm) '8 5 z(cm)

(b) Temperature rise in the 2:2 plane at y = 0.

Figure 6.8. The simulated normalized power deposition and temperature rise generated by
the 2575 element spherically focused phased array in the 3:2 plane at y = 0. The ultrasound
bio-heat simulations are performed in a model containing layers of water, skin, fat, muscle,
and viscera.

112

6.3.4 The effect of nonlinear ultrasound propagation on the simulated

temperature

The nonlinear distortion of therapeutic ultrasound and the influence of nonlinear
harmonics on temperature simulations for hyperthermia is investigated in this section.
Linear and nonlinear models of the pressure generated by three planar phased arrays
are evaluated, and the power depositions and the temperature fields produced by
these arrays under various input intensities are also investigated. The three phased
arrays contain 25 x 25, 50 x 50 and 100 x 100 square elements, respectively. In each
array, the element size is 0.15cm x 0.15cm, and the kerf between adjacent elements
is 0.0lcm. The sizes of these arrays are 4cm x 4cm, 8cm x 8cm, and 16cm x 160m,
respectively. The arrays are driven at 1MHz and electronically focused at (0,0,8cm).
The acoustic wave propagation is evaluated in fatty tissue medium.

The linear axial pressure and the fundamental, the second harmonic, and the third
harmonic from nonlinear simulations are shown in Figure 6.9. Figure 6.9(a) shows
the pressure field generated by the 16cm x 16cm array with an input intensity of
0.11W/cm2. The pressure of the fundamental obtained from the nonlinear simula-
tions is almost coincident with the pressure calculated with the linear model, and
the absolute value of the second harmonic is 9.7dB below that of the fundamental
harmonic at the peak. After 10 seconds of insonation under this intensity, the peak
temperature increases by 34°C if calculated with the linear model, and by 345°C if
calculated with the nonlinear model. Figure 6.9(b) shows the pressures generated
by the 4cm x 4cm array with an input intensity of 1.7W/cm2. Figure 6.9(b) shows
that the fundamental harmonic differs from the linear pressure field near the focal
region, while the second and the third harmonics are about 5.1dB and 8.7dB below
the fundamental at the peak. After 10 seconds of insonation, the simulated peak
temperature rise is 34°C using the linear model, but the peak temperature rise ob-

tained from the nonlinear model is 37°C. The power emitted by the 16cm x 16cm

113

array is 28.16W, and the power emitted by the 4cm x 4cm array is 27.2W. While both
arrays focus at the same depth, consume a similar amount of power, and achieve the
same peak temperature rise under the linear model, the pressure of the fundamental
produced by the larger array is higher at the focus, yet the amplitude of the second
and third harmonics are lower. Based on the comparison of the simulations obtained
from these two arrays, the nonlinear contributions appear to be strongly influenced
by the input intensity for array excitations with equal input power.

The effect of the nonlinear harmonics on the temperature rise is studied quanti-
tively and illustrated in Figure 6.10. The peak intensity and the temperature elevation
above 37°C generated by the 4cm X 4cm planar array are shown in Figure 6.10, where
Figures 6.10(a) and 6.10(b) are the results of steady state temperature simulations,
and Figures 6.10(c) and 6.10(d) are the results of transient temperature simulations
with 10 seconds exposures. The :1: axis in Figures 6.10(a) and 6.10(b) ranges from
0.1W/cm2 to 0.57W/cm2, and the corresponding temperature elevation is from 87°C '
to 50°C. To reach the same range of temperature elevations, the transient simula-
tions in Figures 6.10(c) and 6.10(d) require input intensities between 0.1W/cm2 and
2.48W/cm2. The steady state simulation requires a much lower input intensity than
the transient simulation to achieve the same maximum temperature increase. In
Figure 6.10, the results obtained from linear simulations (solid line) demonstrate
strict linear relationship between the input intensity and the peak output intensity
or the temperature rise. The results obtained from nonlinear simulations (dashed
line) in Figures 6.10(a) and 6.10(b) closely match those calculated from the linear
model. However, the results obtained from the nonlinear model in Figures 6.10(c)
and 6.10(d) diverge from those from the linear model in the high intensity range. The
nonlinear relationship between the input intensity and the output temperature rise is
particularly pronounced in Figure 6.10(d).

The temperature distribution generated by the 4cm x 4cm planar array with the

114

 

 

 

 

 

 

 

 

 

  

 

 

 

2.5 - .
— 1st harmonic
2 ----- 2nd harmonic
A 2 _ ....................... ....................... ....... 3rd harmonic .-
0“: ' ' ---linear
E15“ .......................
93
:3
a 1 .............................................................................................
93
a
05 ........................ ............................................ _,
0 4 6 ‘ 8" ”-1 12

z(cm)

(a) Pressure generated by the 16cm x 16cm planar array.

 

 

 

 

 

    

 

 

 

2 6 6 4
—1st harmonic
‘ i ----- 2nd harmonic
A1 5 .............. 3rd harmonic .
o. _ linear
é '
0.3 1 ........................................................
3 i
U) E
0’ E
g z
0.5, ...................
0 _______ _ .13": """""""""""" i ---»-...'..‘.'..'.:::.-..-:
10 12

z(cm)

(b) Pressure generated by the 4cm x 4cm planar array.

Figure 6.9. The simulated axial pressure along the z axis at a: = y = 0 generated by
(a) a 16cm x 16cm planar phased array and (b) a 4cm x 4cm planar phased array. The
dashed line represents the pressure simulated with the linear model and the rest of the
curves are harmonics simulated with the nonlinear model. The simulations are evaluated

in a homogeneous fatty tissue medium.

115

 

 

 

 

 
    

 
   

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 
 
   

 

 

 

 

 

 

 

 

 

 

 

 

 

3G . . . 50 . .
A —Linear simulations . ,0? —Linear simulations
NE 25 ---Nonlinear simulatigs * - 9 3 40_ ---Nonlinear simulations ,,,,,,,,,,,,,,,,,
g . . B, . . ;
$20 . g
E g 30’
a2, 15- -- E
E 10 .320 .
‘5 .
e a
10-
8 5 5
a 3 1 F 5 3 i
(0.1 04.2 033 0.4 0.5 0.6 CI11 0.2 03 0.4 0.5 0.6
input intensity (chmz) Input intensity (W/cmz)
(a) Output intensity as a function of input in— (b) Output temperature rise as a function of
tensity in steady state simulations. input intensity in steady state simulations.
120 . . . 1 5c; , . .
A —Linear simulations 5 o‘ 1; —Linear simulations .
N - - -Nonlinear simulations ' 3 - - -Nonlinear simulations 2
E .. . ~40. : . , ' 1
\ 90 I - .. - - ~ g: . -
5 2.
g. g 30*
U) ‘C
c 60
a: 2
3:: a 20* ' 1
‘5 30 g
0 a:
; 6 -i ~ '- 3 -: a ‘
0 A n A 1 0 A A i 1
0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5
Input intensity (chmz) Input intensity (W/cmz)
(c) Output intensity as a function of input in- ((1) Output temperature rise as a function of
tensity in transient simulations. input intensity in transient simulations.

Figure 6.10. The simulated peak intensity and temperature rise as functions of the input
intensity produced by the 4cm x 4cm planar phased array. The solid lines represent the
results simulated with the linear model, and the dotted lines are the results simulated
with the nonlinear model. Figures (a) and (b) contain the results from the steady state
simulations, and Figures (c) and (d) contain the results from transient simulations with 10
second exposures.

116

nonlinear pressure simulation model is illustrated in Figure 6.11. Figure 6.11 demon-
strates the simulated temperature rise after a 10 second transient insonation. The
input intensity in this case is 2.5W/cm2. The temperature rise in Figure 6.11(a) is
produced by the fundamental with a maximum value of 47.230 The temperature
rise in Figure 6.11(b) is produced by the second harmonic with a maximum value of
7.7400. The temperature distributions from the linear and the nonlinear models are
compared in Figure 6.12. The top half of Figure 6.12 is the simulated temperature
rise generated by the linear model in the 502 plane at y = 0, and the bottom half of
Figure 6.12 is the total temperature rise produced by the first three harmonics in the
nonlinear model. Figure 6.12 magnifies the difference near the focus. The tempera-
ture contours from the linear and the nonlinear models have similar shapes, but the
temperature from the nonlinear model is higher around the focal spot.

Figure 6.13 compares the difference in peak temperature elevations with the lin-
ear model and the nonlinear model for the three arrays. The :1: axis shows the peak
temperature rise obtained from the linear model, and the differences in tempera-
ture between the linear and nonlinear models are normalized by the maximum linear
temperature rise and plotted on the y axis. The temperatures simulated with the
nonlinear model are always higher than those simulated with the linear model due
to the extra heating from the nonlinear harmonics. Figure 6.13(a) shows the result
from steady state simulations, and Figure 6.13(b) shows the results from transient
simulations with a 10 second exposure time. The solid, dashed and dash-dot lines
represent the results from the 4cm x 4cm, 8cm x 80m, and 16cm x 16cm planar arrays,
respectively. Figure 6.13(a) and Figure 6.13(b) show that the nonlinear contributions
from the smaller arrays produces higher temperature increases than the larger ar-
rays for approximately equal power inputs. For example, at the point corresponding
to a 50°C temperature rise in the transient linear simulation (Figure 6.13(b)), the

temperature produced by the 4cm x 4cm array with the nonlinear model is 10.5%

117

 

10

8
z(cm)

(a) Temperature rise produced by the fundamental.

3

    

4 6 8 1 0
z(cm)

(b) Temperature rise produced by the second harmonic.

Figure 6.11. The temperature rise generated by the 4cm x 4cm planar array with a
2.5W/cm2 input intensity using (a) the fundamental and (b) the second harmonic obtained
from nonlinear simulations.

118

 

8
z(cm)

Figure 6.12. The temperature rise generated by the 4cm x 4cm planar array with a
2.5W/cm input intensity using the linear model (top half of the figure) and the nonlinear
model (bottom half of the figure).

higher than that with the linear model, whereas the temperature produced by the
160m x 16cm array with the nonlinear model is only 1.7% higher than that obtained
with the linear model. Figure 6.13(a) and Figure 6.13(b) also show that, when the
same temperature is achieved in the linear simulations, the additional temperature
rise produced by the nonlinear contributions in the 10 second transient exposure is
an order of magnitude higher than that in the steady state simulation.

The nonlinear simulation results for several planar phased arrays of different sizes
demonstrate that the nonlinear harmonics supplement the heating contribution pro-
vided by the fundamental during hyperthermia. A smaller array requires a higher
input intensity to achieve the same temperature rise that is achieved by a larger
array; therefore, smaller arrays can be associated with more contributions from non—
linear harmonics. The nonlinear contributions are negligible for the steady state
temperature simulations, but the higher harmonics can be significant for transient
simulations at high intensities. These simulation results suggest that nonlinear har—

monies should be considered for pulsed high intensity focused ultrasound (HIFU)

119

 

fi

 

0.012 .
—4x4cm2
0-01"---8x3cm2 _

..... 2%
0.008t’ 162(16601 ::

  
 
 

......................................................

 

 

 

{
’-
'I
‘0
'Q
'|
u
..............................................
.....
. “
.‘ .
"c .
‘3-
't
‘u
"

' :o 1 3.
p............... r..' ............ u .................. :...................s ..................
0.002 ., ,,,,,, ; ; :

 

 

 

Normalized temperature difference
0
O
O
O)

O

0 1‘0 2‘0 3‘0 4‘0 50
Linear temperature rise (degrees)

(a) Steady state simulation results.

 

 

   
  

 

 

 

 

 

 

8 0.12 . t 1 —,r
E; —4x4c:m2
g 01 .. "'8x80m2 .......................................................
U .-.-. 2
g 0.08) 16%16001 .; ........................................................
E :
8.0 06 ............................................................................................
g "‘
S 004. .v' ..........
a)
.5
To 0.02" 3 ”.7 _____________
E ”'3 ........... :._._,_,..-..-
g 0 ..:. ------ 1 """ =. i .=

0 10 2 30 40 50

Linear temperature rise (degrees)

(b) Transient simulation results.

Figure 6.13. The normalized difference between the peak temperaturas obtained from
linear and nonlinear simulations for (a) steady state and (b) 10 second transient exposures
The difference is normalized by the maximum temperature rise obtained from the linear
model and plotted as a function of the linear temperature rise.

120

treatments, but during mild hyperthermia, i.e., treatments with peak temperatures

between 41°C — 43°C, the influence of nonlinear fields may be ignored.

121

CHAPTER 7

Power and Temperature Optimization in

Ultrasound Hyperthermia

7. 1 Introduction

The goal of ultrasound hyperthermia is to increase the temperature in tumor tissue
above a certain threshold value while maintaining the temperature in normal tissue
at sub-therapeutic levels. High intensity focused ultrasound (HIFU) typically heats
tumor targets to above 48 — 50°C with pulsed ultrasound during a short exposure
time (on the order of tens of seconds); mild hyperthermia aims at increasing the tem-
perature in tumors to 41°C to 45°C uniformly for a longer period of time, typically
half an hour to one hour. The thermal effect for HIFU treatment is usually character-
ized with an equivalent exposure time (in minutes) at 43°C [90]. The temperatures
obtained from HIFU and mild hyperthermia are usually simulated with the Pennes
bio-heat transfer equation (BHTE) [86].

Various techniques have been developed to effectively deliver ultrasound to tu-
mor sites and to optimize the temperature distribution in tissue. Lin et al uses
the conjugate-gradient search and the golden section line search to find the opti-
mal set of parameters, including focal depth, scan radius and emitting power, for a

scanned focused ultrasound hyperthermia system that mechanically scans focal spots

122

throughout the target volume [91]. This optimization scheme requires multiple ini-
tial guesses, and often reaches non-unique local minima. A pseudo—inverse pattern
synthesis method [28] obtains a minimum norm least-squares solution to a series of
linear equations that satisfy the Rayleigh-Sommerfeld integral. These equations de—
scribe the relationship between the array excitation signals and the pressures at a
set of specified field points. An intensity gain maximization method [92] is combined
with the pseudo-inverse method to optimize the phase distribution of excitation sig—
nals and achieve maximum gain at desired focal points. Based on the pseudo-inverse
method, McGough [29] develops a method to optimize the temperature field directly
by evaluating the pseudo-inverse of a matrix that relates the power to the temper-
ature. This method can determine the power weights for a series of focal spots to
achieve uniform temperatures in tumors. A waveform diversity method based on the
multi-input multi-output (MIMO) radar system design has been applied to hyperther—
mia [93, 94]. This method achieves maximum and minimum power values in multiple
sites simultaneously through semidefinite programming. Preliminary simulations of
the waveform diversity method have demonstrated successful results in optimizing
the power deposition in an artificial 2D breast model [93].

Temperature optimization for mild hyperthermia is the focus of this chapter. The
simulation results from single focus spot scans are first demonstrated for hyperther-
mia treatment in a deep-seated tumor with a large phased array, where the power
weights for each scan are optimized with the inverse temperature synthesis method
[29]. Second, the mode scanning method generates four focus modes in the same
tumor. Third, the waveform diversity method combined with the mode scanning
technique generates heating patterns for this tumor. Finally, the waveform diversity

method optimizes the power and temperature fields in a breast-tumor model.

123

7 .2 Single focus spot scanning

The focal spots produced by ultrasound phased arrays are usually very small com-
pared to tumor size due to the short wavelength of therapeutic ultrasound. In order
to heat the entire tumor, focal spots are scanned across the target region. If the focal
spots are scanned swiftly and repetitively in space, the accumulated power deposition
is equal to the linear superposition of the power produced by each scan.

Figure 7.1 shows the temperature rise (above 37°C) produced by scanning the
ultrasound beams generated by a spherical-section phased array in a spherical tumor.
The array is comprised of 1444 square 0.24cm >< 0.24cm elements. The geometrical
focus of the array is located at z = 12cm, and the opening angle defined at the
geometrical focus relative to the lateral extent of the array in the a: and the y directions
is 7r / 3. The array is driven by a 1MHz continuous wave excitation. The diameter
of the spherical tumor model is 1.20m, and the tumor is centered at (0,0,15cm) on
the z axis. The surrounding medium is modeled as normal tissue. The attenuation
coefficients for the tumor and the normal tissue are both a = 0.5dB/cm/MHz. Six
focal spots are evenly distributed on the circumference of the tumor in the plane at
z = 15cm. The pressure generated by each scan is calculated with the fast nearfield
method in the my plane at z = 80m and then propagated with the angular spectrum
approach in a 15cm x 15cm x 120m volume. The spatial sampling rate is 6 = 0.0750m,
and the F FT calculations are performed on an N x N grid with N = 512. The power
deposition from each scan is calculated from the pressure field, and the results are
then superposed. The total power deposition is subsequentlly used as the input to the
bio-heat transfer equation. The simulated temperature rise in the yz plane at a: = 0
is shown in Figure 7.1, where the tumor cross section is denoted a circle, and the
4°C and 5°C isothermal contours are represented by the elliptical curves. The peak
temperature rise in the computational volume is 6°C. Figure 7.1 shows that the small

tumor is uniformly heated with single focus spot scans, but the normal tissue between

124

    

"7'58 10 12 14 16 18 20
z(cm)

Figure 7.1. The simulated temperature rise in the yz plane evaluated at :1: = 0 produced
by a 1444 element spherical-section array that is focused within a 1.2cm diameter spherical
tumor. The heating pattern is created by scanning six focal spots repeatedly around the
circumference of the tumor.

the tumor and the phased array applicator also suffers from intervening tissue heating
near z=13cm and z=14cm.

As the size of the tumor and the required number of focal spots increases, there
is a greater need for optimization to determine the location of the focal spots and
the power weight for each focus. The direct thermal inverse method [29] optimizes
the power weight at each focal spot by solving a group of linear equations that relate
the power to the temperature distribution. The total power deposition is modeled
as the linear superposition of the powers produced by the single focus scans, and

the relationship between the temperature and the peak power at each focal spot is

125

expressed in matrix form as,

 

”T17 'b1,1 b1,2 b1,m- “(11 -

T2 52,1 52,2 b2, (12

' = : m : = Bn x QO x 11 (7'1)
b Tn _ . bnql bn,2 ... bn’m _[ _ qm _]

 

 

 

 

 

where m is the number of focal spots and n is the number of selected measurement
points including focal spots and some additional points in normal tissue, qm x 1 is a
vector of power deposition values at all focal spots, T n x 1 is a vector of temperature
values at all measurement points. The matrix 3,, x m consists of entries b1, j that
relate the temperature produced by the jth scan evaluated at the ith measurement
point. The vector qm x 1 contains the power weights to be calculated. The vector
Tn x 1 contains the desired temperature rise at all measurement points, for instance,
6°C at focal points and 4°C at points in normal tissue. Computing all of the entries in
B is time-consuming because the 3D power deposition and the temperature response
need to be calculated for each scan. An interpolation based approximation scheme is
instead used to populate the matrix B. A reference temperature field with a single
on-axis focus is calculated in a 3D grid. The other temperature fields with steered
focal spots are obtained from rotating and interpolating the reference field. The
interpolation are calculated according to the locations of the measurement points
relative to the focus in each scan. This approximation is applicable if the variation
of the size of the focus is small throughout the target region and if the grating lobes
are small compared to the main lobe at the focus. The fields produced by ultrasound
therapy arrays usually satisfy these constraints. After the entries in the B matrix and
the T vector are determined, the unknown vector q is computed by solving Equation
7.1 in a least-square sense subject to a non-negative constraint.

To examine the heating strategy of single focus spot scans, the power deposition
generated by a 1444 elements spherical-section array is modeled in a 3cm diameter

spherical tumor. The center of the array is located at the origin, and the center of

126

the tumor is located at (0,0,15cm). A total of 172 focal spots are evenly distributed
within the spherical tumor with a 0.45cm spacing between adjacent focal spots. The
total power deposition is the superposition of the powers produced by individual
scans. The temperature field is calculated with the bio-heat transfer equation using
the total power deposition as the input. Figure 7.2 contains the resulting temperature
distribution evaluated in two orthogonal planes when all focal spots have the same
power weights. Figure 7.2(a) is the temperature distribution across the center of
the tumor in the my plane at z = 15cm, and Figure 7.2(b) shows the temperature
distribution in the yz plane at a: = 0. The circle in Figure 7.2 illustrates the boundary
of the spherical tumor. The isothermal contours associated with the 4°C and 5°C
temperature rise are also indicated. Figure 7.2(a) shows that only the center of the
tumor is heated to the desired temperature. Figure 7.2(b) shows that the region
close to the z axis contains the highest temperature increase due to the accumulation
of power along the scan paths, whereas the off-axis regions are below the threshold
temperature.

In order to increase heating near the tumor periphery, the power weight at each
focus is optimized by the direct thermal inverse method using Equation 7.1. The
distribution of control points is shown in Figure 7.3, along with illustrations of the
spherical-section array and the spherical tumor. Figure 7.3 contains 249 control points
that are densely distributed in a 1.6cm x 1.6cm x 4cm volume in the normal tissue
region proximal to the tumor, and the spacing between adjacent tissue control points
is 0.4cm. A total of 172 focal spots are evenly distributed within the spherical tumor
with a 0.450m spacing between adjacent focal spots. The optimization objective is to
increase the temperature at the tumor control points to 43°C while maintaining the
temperature at the tissue control points at 40°C to limit intervening tissue heating.
The power weights obtained from the direct thermal inverse method are multiplied

by the power deposition obtained for each focus before the individual contributions

127

    

"73.5 —5 -2.5 o 2.5 5 7.5
y(cm)

(a) Temperature rise in the xy plane evaluated at z = 15cm.

7.5 6

 

-7.
581012141618 20
z(cm)

(b) Temperature rise in the yz plane evaluated at a: = 0.

Figure 7 .2. The simulated temperature rise produced by a 1444 element spherical—section
array in a spherical tumor with a 3cm diameter. The heating pattern is created by spot
scanning 172 focal points that are evenly distributed in the tumor with uniform weighting.

128

 

 

 

 

Figure 7 .3. The distribution of control points within the spherical tumor and in the
normal tissue region for direct thermal inverse optimization. The 3cm spherical tumor is
represented by the sphere, and the 1444 element spherical-section phased array is located
at the bottom of the figure.

are superposed. The resulting power deposition across the center of the tumor in
the my plane evaluated at z = 15cm is shown in Figure 7.4, where the power field
is normalized by the peak value in the 3D volume. Although the focal spots are
distributed uniformly in the tumor, Figure 7.4 shows that the power weights on the
focal points near the periphery of the tumor are larger than those located closer to
the center of the tumor. The uneven power distribution within this cross section is a
result of the direct thermal inverse optimization.

Figure 7.5 shows the temperature rise (above 37°C) calculated by the BHTE using
the power deposition in Figure 7.4 as the input. Figure 7.5(a) shows the temperature

rise in the my plane evaluated at z = 150m, where the 4°C isothermal contour is

129

7.5

—5

    

73.5 -5 —2.5 o 2.5 5 7.5
y(cm)

Figure 7.4. Simulated power deposition in the my plane evaluated at z = 15cm produced
by the 1444 element spherical-section array in a spherical tumor with a 3cm diameter. The
heating pattern is created by spot scanning 172 foci within the tumor, where the power
weight for each focus is optimized with the direct thermal inverse method.

almost coincident with the tumor periphery, and the peak temperature in this plane
is about 45°C. Figure 7.5(b) shows the temperature rise in the yz plane evaluated at
a: = 0, where most of the tumor region achieves a 4°C temperature rise, and the peak
temperature increase of 6°C occurs between the tumor and the phased array. Figure
7.6 shows the 3D spherical tumor model that is illustrated as a sphere and the isother-
mal surface that corresponds to temperature rise above 4°C. A significant portion of
the tumor achieves the desired temperature elevation, whereas a substantial amount
of the normal tissue between the tumor and the array are inevitably overheated. The
results in Figure 7.5 and Figure 7.6 show that single focus spot scanning using opti-
mized power weights achieves a better result than uniform weights. However, due to
severe intervening tissue heating, the single focus spot scan appears to be suboptimal

for deep regional hyperthermia.

130

  
  

 

4:75 -5 —2.5 o 2.5 5 7.5
y(cm)

(a) Temperature rise evaluated in the xy plane at z = 15cm.

7.5 6

 

O

'7'58 10 12 14 16 18 20
z(cm)

(b) Temperature rise evaluated in the yz plane at a: = 0.

Figure 7.5. The simulated temperature rise produced by the 1444 element spherical-
section array in a spherical tumor with a 3cm diameter. The heating pattern is created by
spot scanning 172 foci within the tumor, where the power weight for each scan is optimized
with the direct thermal inverse method.

131

 

 

 

 

y(cm)

Figure 7 .6. A transparent mesh that indicates the size and location of the spherical tumor
model and an isothermal surface that indicates the region that achieves a minimum of 4°C
temperature rise. The heating pattern is created by spot scanning 172 foci within the tumor,
where the power weight for each focus is optimized with the direct thermal inverse method.

132

7 .3 Multiple focus mode scanning

In hyperthermia, multiple focus beamforming can be superior to single focus beam-
forming for at least two different reasons. First, the chance of over-heating one spot
and generating undesirable cavitation is reduced because the power is distributed
among several focal spots in a multiple focus scan. Second, less phase channels may
be needed because the tumor region is covered by fewer scans when targeting mul-
tiple locations simultaneously. Mode scanning method [55], one of the few multiple
focus beamforming techniques, utilizes the symmetry of the phased array to generate
symmetric focal spots while canceling the pressure along planes of symmetry. An
array that also contains multiple planes of symmetry is separated into M sections,
and the array elements in each section are indexed according to the planes of symme-
try so that the distances from one pair of elements to the dividing plane are equal.
Excitation signals with rotational phases are applied to symmetric elements in each
section so that M focal spots are generated simultaneously. The rotational signal is
defined so that if w1(t) is applied to array section 1, the excitation signal for the lth
section of the array is w1(t) = 101(1.‘)e~7'(l — I)“. The rotational signals with a period
of 1r cause pressure cancellation on the planes of symmetry.

The mode scanning method is applied to the 1444 element phased array to generate
multiple focal spots and heat the 3cm spherical tumor centered at z = 15cm. A total
of 172 focal spots are used according to the locations described in section 7 .2. A
group of four focal spots are produced in each scan, and these spots are symmetric
with respect to the planes at :1: = 0 and y = 0. The steady state power deposition
is equal to the superposition of the power fields produced by all scans. The direct
thermal inverse method also optimizes the power weight at each focus. Only one
quarter of the control points are needed for optimization due to the symmetry of the
control points.

Figure 7.7 shows the overall power deposition in the my plane evaluated at z =

133

 

15cm. As a result of the direct thermal inverse optimization, the focal spots near the
periphery of the tumor have higher power weights than those closer to the center.
Figure 7.7 also shows that the power is zero in the mz plane at y = 0 and in the yz
plane at m = 0 because of the cancellation of the pressure fields in these planes. Figure
7.8 shows the temperature rise (above 37°C) calculated by the BHTE using this power
deposition as the input. Figure 7.8(a) shows the temperature distribution in the my
plane evaluated at z = 15cm, and Figure 7.8(b) is the temperature distribution in
the yz plane evaluated at m = 0. The results in Figure 7.8 and Figure 7.5 are similar,
where a broad region in the proximal part of the tumor is heated above 4°C, but
the peak temperature occurs in the normal tissue region between the tumor and the
phased array. Figure 7.9 shows the 3D isothermal surface with 4°C temperature
rise. Relative to Figure 7.6, the result in Figure 7 .9 contains a similar amount of
intervening tissue heating. For this combination of parameters, multiple focus mode
scanning and single foCus spot scanning are roughly equivalent in the deep-seated

tumor target evaluated here.

134

7.5

2.5

x(cm)
O

    

'735 -5 -2.5 o 2.5 5 7.5
y(cm)

Figure 7.7. The power deposition in the my plane evaluated at z 2 15cm produced by the
1444 element spherical-section array in a spherical tumor with 3cm diameter. The heating
pattern is created by mode scanning 172 foci within the tumor, where the power weight for
each focus is optimized with the direct thermal inverse method.

7 .4 Power optimization with waveform diversity method

7 .4.1 Introduction

The waveform diversity method is originally used for the multi—input multi-output
(MIMO) design in radar signal processing. The waveform diversity method can gen-
erate diverse beam patterns to match the desired power distribution. The effect of
this method is equivalent to multiple focus beamforming in the field of ultrasound.
Preliminary results have demonstrated the superior performance of the waveform di-
versity method in hyperthermia through computer simulations of a two dimensional
breast model with sphere and hemisphere geometries [93, 94]. While other single
focus and multiple focus temperature synthesis methods solve the power or temper—
ature optimization through multiple steps, i.e., from the temperature distribution to

the power distribution, then to the pressure, and finally the excitation signals on

135

    

7:95 —5 —2.5 o 2.5 5 7.5
y(cm)

(a) Temperature rise in the my plane evaluated at z = 15cm.

7.5 6

-2.5

-5

 

'7'58 10 12 14 16 18 2o
z(cm)

(b) Temperature rise in the yz plane evaluated at m = 0.

Figure 7.8. The simulated temperature rise produced by the 1444 element spherical-
section array in a spherical tumor with a 3cm diameter. The heating pattern is created by
mode scanning 172 foci in the tumor, where the power weight for each focus is optimized
with the direct thermal inverse method.

136

 

 

 

    

0-1

-2 -2

y(cm) x(cm)
Figure 7.9. A transparent mesh that indicates the size and location of the spherical tumor
model and an isothermal surface that indicates the region that achieves a minimum of 4°C
temperature rise. The heating pattern is created by mode scanning 172 foci within the
tumor, where the power weight for each focus is optimized with the direct thermal inverse
method.

137

array elements, the waveform diversity method solves a convex optimization problem
that connects the power pattern directly to the covariance matrix of the excitation
signals. The excitation signals are then synthesized from this covariance matrix. The
optimization problem can be solved with public domain convex optimization software
like SeDuMi [95]. When the two dimensional models are extended to three dimen-
sional models that use larger phased arrays, the simulations sometimes suffer from
insufficient memory and long computation times. Therefore, the array elements and
the control points are divided into groups with the mode scanning technique, and the

solution produces symmetric results.

7.4.2 Formulation

The complex pressure generated by a phased array at an arbitrary field point r is

denoted by
M
190‘, t) = Z pm(r)wm(t), (7-2)
= 1

where M is the number of array elements, pm(r) is the complex pressure produced
by the mth array element at location r when the element is excited by a sinusoidal
signal with unit amplitude and zero phase, and wm(t) is the complex excitation

h array element. The waveform diversity method defines each

signal applied to the mt
excitation signal wm(t) as a sinusoidal signal exp(j27rf0t) multiplied by a complex
envelope signal. The complex envelope signal is assumed to be a random variable
over time, and the envelope signal varies more slowly than the sinusoidal signal with
frequency f0. If wm(t) is discretized with a sampling rate much lower than fo, the

excitation signal can be regarded as a constant amplitude continuous wave signal

between each time sample. After discretization, Equation 7.2 becomes

M
p<r.n>= Z pm<r>wm<n>=mr>w<n1 n=1,2.....N. (7.3)
= 1

138

where n is the index of the time samples that range from 1 to N, 13(r) is an 1 x N
vector comprised of all values of pm(r), and iii—(n) is an N X 1 vector comprised of
the excitation signals applied to all array elements at time sample n. The power at

location r is proportional to the auto-covariance of pressure waveforms,

620‘) = %E{p(r1n)p*(rm}, (7.4)

where a is the attenuation coefficient, p and c represent the medium density and the
Speed of sound, respectively, and E denotes the expectation. By plugging Equation

7.3 into Equation 7.4, the expression for the power becomes

620‘) = fi(r)E{'U7(n)W(nl}17*(r) = 1700131570, (7-5)

where R is the covariance matrix of the excitation signals “117(n). R is independent of
the location of the observation points. Although 117(7).) represents a vector of random
variables, the covariance matrix R uniquely defines a power distribution. One way to

synthesize {711111)} from R is to apply
w(n) = Rl/2z(n), (7.6)

where R”2 is the square root of R, and the sequences 3(7).) represent M 1 x N i.i.d.
random vectors with zero means and an M x M identity matrix as the covariance
matrix. Given that N is sufficiently large, the resulting 13(12) are M random variables
with zero means and a covariance matrix that approximates R.

The goal of hyperthermia is to distribute relatively uniform temperature within
tumor region and minimize the temperature increase in normal tissue. If the power

distribution is instead optimized, this problem is defined as,

min“; —t
5.1:. 13(ro)RI‘9*(ro) - flulRWu) > t, u E 0N
W!) 19*(11) Z 0-9§(1‘0)Rfi*(r0)1 V 6 QT (7 7)
13(V)RI‘9*(V) S 1.123(rolR1‘9*(r0)1 V 6 RT '
R 2 0,
Rmm=1/M, m=1,2,...,M

139

where )1 denotes an arbitrary point in the normal tissue region ON, V denotes an
arbitrary point in the tumor region QT, and r0 is a typical point within QT, i.e., the
center of the tumor. The cost function maximizes the difference between the power
at r0 and the power at all tissue control points )1. The first constraint (where t > 0)
guarantees that the power at r0 is always larger than the power at p. The second and
the third constraints ensure that the power at the tumor control points 11 is no more
than 1.1 times and no less than 0.9 times the power at r0. The fourth constraint
defines the matrix R as a positive semi-definite matrix. The fifth constraint ensures
that R has uniform diagonal entries, which guarantees that all of the array elements
emit an equal amount of power. The function p(r) = p(r)Rp*(r) is convex, there-
fore, the optimization problem in Equation 7.7 can be solved by convex optimization
softwares such as SeDuMi [95] and SDPT3 [96, 97] and the corresponding MATLAB
interface tools Yalmip [98] and CVX [99].

7 .4.3 Simulation for a deep seated spherical tumor

The waveform diversity method optimizes the power deposition generated by the 1444
element spherical-section array to heat a deep-seated spherical tumor. To reduce the
computational cost due to the large number of transduCer elements in the array, the
spherical-section array is divided into four section according to the array symmetry.
The control points are also accordingly divided into four symmetric groups. The exci-
tation signals for one quarter of the array elements are calculated with the waveform
diversity method, the excitation signals are then multiplied by the rotational phases
ejm7r (m = 2, 3, 4) determined by the mode scanning technique, and the results are
applied to the elements in the other three quarters of the array. Initial optimization
results show that intervening tissue heating tends to occur distal to tumor along the
four scan paths that connect each symmetric quadrant of the array to the tumor.

The tissue control points are distributed in a 0.3cm x 0.3cm X 1.5cm volume distal

140

20
16

z(cm)

(DO-km

        

6

3 3

o o
-3 -3
y(cm) '6 '6 x(cm)

Figure 7 .10. The 1444 element spherical-section array, the 3cm spherical tumor and
the control points used in the waveform diversity Optimization. One quarter of the array
elements and control points that are located in the same quadrant of the coordinates are
used in the optimization, that is, a total of 361 array elements (painted in red), 43 tumor
control points, and 67 tissue control points.

to the tumor and in a 0.6cm x 0.6cm x 3cm volume proximal to the tumor in one of
the quadrants. The spacing of the tissue control points distal to the tumor is 0.15cm,
and the spacing of the tissue control points proximal to the tumor is 0.60m. The
tumor control points are evenly distributed in the spherical tumor with a 0.45cm
spacing between adjacent points. The control point distribution is shown in Figure
7.10, where one quarter of the array (361 elements) is shown in red, and 43 tumor
control points and 67 tissue control points are indicated by dots. The active array
elements and the control points are shown in the quadrant that has negative :1: and y
coordinates.

Figure 7.11 shows the power deposition calculated with the waveform diversity

141

7.5

    

"735 -5 -2.5 o 2.5 5 7.5
y(cm)

Figure 7.11. Simulated power deposition in the my plane evaluated at z = 15cm produced
by the 1444 element spherical-section array in a spherical tumor with a 3cm diameter. The
heating pattern is created by the waveform diversity method combined with mode scanning.
A total of 361 array elements, 43 tumor control points, and 67 tissue control points are used
in the optimization.

method in the my plane located at z = 15cm. In Figure 7.11, the focal spots are
somewhat smeared, and the resulting tumor power deposition is relatively uniform
in the tumor in this cross section. This result is different from the single focus spot
scanning or the multiple focus mode scanning results, where the power gradient is
large between focal points and other regions. The power distribution in Figure 7.11
is symmetric with respect to to the m and y axes, and the power along these axes is
zero because of the destructive interference produced by mode scanning.

Figure 7.12 shows the temperature rise (above 37°C) calculated with the BHTE
using the power deposition shown in Figure 7.11 as the input. Figure 7.12(a) shows
the temperature distribution in the my plane evaluated at z = 150m, where the
temperature is increased by 4 —— 6°C in most parts of the tumor. Figure 7.12(b)

shows the temperature distribution in the yz plane evaluated at m = 0, where the

142

peak temperature is located inside of the tumor, and the intervening tissue heating
between the tumor and the array is eliminated in this plane. A significant amount
of normal tissue distal to the tumor is also overheated due to the intersection of the
beams produced by mode scanning.

Figure 7 .13 shows 4°C isothermal surface superposed on the spherical tumor
model. Figure 7.13 shows that the a large portion of the tumor reaches the tar-
get temperature. The hot spot behind the tumor is only moderate. The intervening
tissue heating between the array and the tumor is also largely reduced relative to the
results from the single focus spot scanning and multiple focus mode scanning. The
results in Figures 7 .12 and 7.13 show that the waveform diversity method has signif-

icant promise for deep regional hyperthermia applications that involve large arrays.

The computational cost of the waveform diversity optimization is an important
issue. For thermal therapy, a densely populated large array aperture filled with
small elements (less or equal to half wavelength) are desirable because these features
influence intervening heating and the formation of grating lobes. Therefore, phased
arrays designed to heat deep seated tumors often contain thousands of elements.
Furthermore, for the waveform diversity method, a dense distribution of both tumor
and tissue control points are needed so that the optimization goals can be achieved.
According to Equation 7 .7 , the number of constraints is equal to N N + 2NT + M, and
the number of unknowns is N N + 2NT + M 2, where M is the number of transducer
elements, and N N and NT are the number of control points in normal tissue and tumor,
respectively. An intermediate matrix of size (N N + 2NT + M) x (N N + 2NT + M 2)
is constructed during the optimization process, and the number of nonzero entires in
the matrix is (N N + 2NT + M) x (1 + M 2). There is a huge amount of data, which
prevents the waveform diversity method from directly optimizing the power deposition

generated by a large array comprised of many elements. For example, optimizing the

143

    

'7475 —5 —2.5 o 2.5 5 7.5
y(cm)

(a) Temperature rise in the my plane evaluated at z = 15cm.

7.5 6

 

"7'58 10 12 14 16 18 20
z(cm)

(b) Temperature rise in the yz plane evaluated at m = 0.

Figure 7.12. Simulated temperature rise produced by the 1444 element spherical-section
array in a spherical tumor with a 3cm diameter. The heating pattern is created by the
waveform diversity method combined with mode scanning. 361 array elements, 43 tumor
control points, and 67 tissue control points are defined for optimization.

144

 

 

 

 

-1
y(cm) '2 '2

Figure 7 .13. The 3cm spherical tumor model and the 4°C isothermal surface. The heat-
ing pattern is created by the waveform diversity method combined with mode scanning.
361 array elements, 43 tumor control points, and 67 tissue control points are defined for
optimization.

145

power generated by the 1444 element array in 172 tumor control points and 268 tissue
control points requires more than 4GB of RAM for storage purposes only. Storing
and processing this amount of data is very costly even if compression is performed on
the data. In contrast, with 361 array elements, 43 tumor control points, and 67 tissue
control point, the required RAM is reduced to 67MB. The memory and computational
time are largely reduced when the mode scanning technique is combined with the
waveform diversity method. The waveform diversity optimization combined with the
mode scanning method is performed on a 64 bit dual-core computer workstation with
2.4GHz CPU and 4GB RAM. The simulation time is roughly half an hour for the

problem described in the text.

7 .4.4 Simulation for a breast model with an embedded tumor

The power distribution in a breast model with an embedded spherical tumor is op-
timized with the waveform diversity method. The breast model is obtained from
segmented MRI images of a female patient. Breast tissue, water, and tumor tissue
are identified in the simulation model. A spherical tumor mode with 3cm diameter
is used.

A spherical-section array comprised of 256 square elements is defined for breast
tumor heating. The array is driven by 1MHz continuous wave. The size of each
square element is 0.24cm (1.6x\). The geometrical focus of this array is located at
z = 8cm, and the opening angle defined at the geometric focus is 7r /4. The center of
the array is located at the origin of the coordinates, and the normal passing through
the center of the array is coincident with the z axis. The computational grid extends
from —7cm to 7cm in the m and the y directions, and the grid extends from 2cm to
10cm in the z direction. The volume is sampled at an interval of 6 = 0.05cm, which
results in a 281 X 281 X 161 computational grid. The pressure field is simulated with

the angular spectrum approach, where the model is approximately by parallel layers

146

of water and fatty tissue. The source pressure in water is simulated with the fast
nearfield method in a 14cm X 14cm source plane located at z = 2cm. The angular
spectrum approach then propagates the pressure to a plane at z = 2.5cm, where the
breast tissue (including skin) starts to contact with water. Reflection and refraction
of ultrasound waves are calculated by the angular spectrum simulations, assuming
the breast / water interface is approximately planar. The angular spectrum approach
then propagates the pressure in breast tissue using parameters for fat. The tumor
model has the same attenuation as fatty tissue.

The power distribution result obtained from the waveform diversity optimization
is shown in Figure 7.14, where Figure 7.14(a) shows the power deposition in the
mz plane evaluated at y = 0, and Figure 7.14(b) shows the power deposition in
the yz plane at m = 0. 118 tumor control points and 73 tissue control points are
selected for waveform diversity optimization. The tumor control points are uniformly
distributed in the spherical tumor with a sample spacing of 0.45cm. The tissue control
points are uniformly distributed in a 6cm X 6cm X 1.5cm volume with a 0.8cm sample
spacing located in the breast tissue region between the tumor and the phased array
applicator. In Figure 7.14, the power deposition in water is extremely small due to
the small value of the absorption coefficient in water, and a significant amount of
power is delivered to the tumor target. The corresponding temperature distributions
in the mz plane and the yz plane are shown in Figure 7.15(a) and Figure 7.15(b),
respectively. The temperature field is calculated with the bio-heat transfer equation,
and the temperature increase in water is maintained at 0°C. A temperature increase
above 4°C is achieved in a large portion of the tumor shown in Figure 7.15. Figure
7.16 shows the heating effect in 3D, where the spherical-section array is located at the
bottom, the green surface represents the boundary of the breast, the purple sphere
represents the tumor, and the red surface is the 4°C isothermal surface. Figure 7.16

shows that the temperature in a large portion of the tumor is increased 4°C or more

147

without excessive intervening tissue heating.

148

x(cm)

 

z(cm)
(a) Power deposition in the mz plane evaluated at y = 0.

7

11(le

 

6
z(cm)

(b) Power deposition in the yz plane evaluated at m = 0.

Figure 7.14. The power deposition in (a) the m2 plane evaluated at y = 0 and (b) the
yz plane evaluated at m = 0 produced by a 256 element spherical-section array in a breast
model with an embedded 3cm diameter spherical tumor. The heating pattern is optimized
by the waveform diversity method.

149

x(cm)

 

6
z(cm)

(a Temperature rise in the .732 plane evaluated at
y = 0.

   

—7
2 4 6 8 10
z(cm)
(b) Temperature rise in the yz plane evaluated at

m=0.

Figure 7.15. The temperature rise in (a) the mz plane evaluated at y = 0 and (b) the
yz plane evaluated at m = 0 produced by a 256 element spherical-section array in a breast
model with an embedded 3cm diameter spherical tumor. The heating pattern is optimized
by the waveform diversity method.

150

0 RM
0 '0 '0:
4.2.3.2?
O
‘.o . .0..o.~:'
. o .9 o on.
x. O o 0‘...»

   

z(cm)
mo M AC) 003

3

  

0
3
y(cm) -6 -6 x(cm)

Figure 7 .16. 3D illustration of the 256 element spherical-section array, the breast model,
the embedded spherical tumor, and the 4°C isothermal surface. The heating pattern is
optimized with the waveform diversity method.

151

CHAPTER 8

Conclusion

8.1 Summary

This thesis contains analysis that show the angular spectrum approach is a rapid
and robust method for acoustic simulations. Implementation details of the angular
spectrum approach applied to phased array simulations are described. Methods to
improve thermal therapy efficiency are demonstrated through power optimizations
and temperature simulations. In chapter 2, several conventional integral methods
are presented for pressure field calculations in linear and homogeneous media. These
methods include the point source superposition method of the Rayleigh-Sommerfeld
integral, the spatial impulse response method, and the fast nearfield method. These
methods are relatively straightforward but are usually computationally inefficient in
simulating large ultrasound phased arrays in three dimensional domains. In chapter
3, the angular spectrum approach is investigated for nearfield pressure simulations
of single transducers. The angular spectrum approach is implemented for both the
spatial and spectral propagators. Results show that the spatial propagator computes
the pressure field more accurately, whereas the spectral propagator is more time
efficient. Furthermore, the results obtained from the spatial and spectral propagators
achieve similar accuracy in simulations that involve attenuation or source apodization.

In chapter 4, ultrasound phased array simulation in a linear and homogeneous medium

152

is investigated. Several ultrasound phased array configurations are demonstrated and
different beamforming methods are introduced. One planar 2D array and two concave
arrays are simulated with the Rayleigh-Sommerfeld integral approach and the fast
nearfield method. Results show that the fast nearfield method achieves the same
numerical accuracy in much less time than the Rayleigh-Sommerfeld integral. In
chapter 5, the angular spectrum approach is evaluated for simulations of ultrasound
phased arrays. Results show that the pressure source plane produces more accurate
simulation results than the velocity source plane in angular spectrum computations.
Small errors are produced when the pressure source plane is located one wavelength
away from the array aperture and when the source plane size is larger than the
array aperture. Furthermore, the output errors from angular spectrum computations
asymptotically approach a limiting value as the number of abscissas used for pressure
source plane calculations increases. To achieve this error limit in the computed 3D
pressure field, fewer abscissas are required by the fast nearfield method than the
Rayleigh-Sommerfeld integral approach for 2D pressure source plane calculations due
to the rapid convergence of the fast nearfield method. The computational efficiency
of the angular spectrum approach is also demonstrated. In chapter 6, acoustic and
thermal simulations are performed for a layered tissue model for nonlinear media.
Results show that the heating produced by nonlinear harmonics is small in mild
temperature hyperthermia, but nonlinear propagation is important for high intensity
focused ultrasound ablation. In chapter 7, the temperature fields produced by single
focus spot scans and multiple focus scans are simulated. The single focus beam
pattern is synthesized with the phase conjugation method, and the power weights
at the focal spots are optimized with the direct thermal inverse method. Results
show that spot scanning can cause severe intervening heating in normal tissue, which
limits the penetration depth and the size of treatable tumors. Multiple focus beam

patterns are synthesized first with the mode scanning method, then with the waveform

153

diversity method. The waveform diversity method is combined with mode scanning
to reduce the computational cost. Results Show that the waveform diversity method

combined with mode scanning show great promise in deep region hyperthermia.

154

BIBLIOGRAPHY

[1] L. J. Anghileri and J. Robert. Hyperthermia in cancer treatment. CRC Press,
Inc., 1st edition, 1986.

[2] R. B. Romer. Engineering aspect of hyperthermia therapy. Annu. Rev. Biomed.
Eng, 01:347—376, 1999.

[3] J. W. Hand and J. R. James. Physical techniques in clinical hyperthermia. Re-
search studies press Ltd., Herts, England, 1986.

[4] F. K. Storm. Hyperthermia in cancer therapy. G. K. Hall Medical Publishers,
Boston, 1983.

[5] T. L. Szabo. Diagnostic ultrasound imaging: inside out. Elsevier Academic
Press, 2004.

[6] P. M. Corry, W. J. Spanos, Tilchen E. J ., B. Barlogie, H. T. Barkley, and E. P.
Armour. Combined ultrasound and radiation therapy treatment of human su-
perficial tumors. Radiology, 145(1):165—169, 1982.

[7] K. Hynynen. The role of nonlinear ultrasound propagation during hyperthermia
treatments. Med Phys, 18(6):1156—63, 1991.

[8] C. Damianou and K. Hynynen. Focal spacing and near-field heating during
pulsed high temperature ultrasound therapy. Phy. Med. Biol, 19(9):777—787,
1993.

[9] T. V. Samulski, W. J. Grant, J. R. Oleson, K. A. Leopold, M. W. Dewhirst,
P. Vallario, and J. Blivin. Clinical experience with a multi-element ultrasonic

hyperthermia system: Analysis of treatment temperatures. Int. J. Hyperthermia,
6(5):909—922, 1990.

[10] E. G. Moros, W. L. Straube, and M .J. Myerson. A reflected-scanned ultrasound
system for external simultaneous thermoradiotherapy. IEEE Trans. Ultrason.
Ferroelect. Freq. Contr., 43(3):441-449, 1996.

155

 

[11] E. G. Moros, X. Fan, and W. L. Straube. An investigation of penetration depth
control using parallel opposed ultrasound arrays and a scanning reflector. J.
Acoust. Soc. Am., 101(3):1734—1741, 1997.

[12] P. Novak, E. G. Moros, W. L. Straube, and M. J. Myerson. SURLAS: A new
clinical grade ultrasound system for sequential or concomitant thermoradiother-
apy of superficial tumors: Applicator description. Med. Phys, 32(1):230—240,
2005.

[13] E. S. Ebbini, S. i. Umemura, M. Ibbini, and C. A. Cain. A cylindrical-section ul-
trasound phased-array applicator for hyperthermia cancer therapy. IEEE Trans.
Biomed. Eng, 35(5):561—572, 1988.

[14] E. S. Ebbini and C. A. Cain. A spherical-section ultrasound phased array applica-
tor for deep localized hyperthermia. IEEE Trans. Biomed. Eng, 38(7):634—643,
1991.

[15] M. Pernot, J. F. Aubry, M. Tanter, J. L. Thomas, and M. Fink. High power
transcranial beam steering for ultrasonic brain therapy. Phys. Med. Biol,
48(16):2577—2589, 2003.

[16] C. A. Cain and S. I. Umemura. Concentric-ring and sector-vortex phased-array
applicators for ultrasound hyperthermia. IEEE Trans. Microwave Theory Tech,
MTT-34(5):542—551, 1986.

[17] S. i.. Umemura and C. A. Cain. The sector-vortex array: Acoustic field synthesis
for hyperthermia. Trans. Ultrason. Ferroelect. Freq. Contr., 36(2):249—257, 1989.

[18] L. E. Kinsler, A. R. Frey, A. B. Coppens, and J. V. Sanders. Fundamentals of
Acoustics. John Wiley and Sons, New York, fourth edition, 2000.

[19] J. Zemanek. Beam behavior within the nearfield of a vibrating piston. J. Acoust.
Soc. Am., 49(1):181—191, 1971.

[20] F. Oberhettinger. On transient solutions of the “baffled piston” problem. Journal
of Research of the National Bureau of Standards-B. Mathematics and Mathemat-
ical Physics, 65B(1):1—6, 1961.

[21] J. C. Lockwood and J. G. Willette. High-speed method for computing the exact
solution for the pressure variations in the nearfield of a baffled piston. J. Acoust.
Soc. Am., 53(3):735—741, 1973.

[22] M. Arditi, F. S. Foster, and J. W. Hunt. Transient fields of concave annular
arrays. Ultrason. Imaging, 3(1):37—61, 1981.

156

 

[23] R. J. McGough, J. F. Kelly, and T. V. Samulski. An efficient grid sectoring
method for calculations of the nearfield pressure generated by a circular piston.
J. Acoust. Soc. Am., 115(5):1942—1954, 2004.

[24] R. J. McGough. Rapid calculations of time-harmonic nearfield pressures pro-
duced by rectangular piston. J. Acoust. Soc. Am., 115(5):1934—1941, 2004.

[25] J. F. Kelly and R. J. McGough. A time—space decomposition method for cal-
culating the nearfield pressure generated by a pulsed circular piston. J. Acoust.
Soc. Am., 53(6):1150—1159, 2006.

[26] D. Chen and R. J. McGough. A fast nearfield method for calculations of time-
harmonic and transient pressures produced by triangular pistons. J. Acoust. Soc.
Am., 120(5):2450—2459, 2006.

[27] M. F. Hamilton and D. T. Blackstock. Nonlinear acoustics. Academic Press,
San Diego, 1998.

[28] E. S. Ebbini and C. A. Cain. Multiple-focus ultrasound phased-array pattern
synthesis: Optimal driving signal distributions for hyperthermia. IEEE Trans
Ultrason Ferroelectr Freq Control, 36(5):540—548, 1989.

[29] R. J. McGough, E. S. Ebbini, and C. A. Cain. Direct computation of ultra-
sound phased—array driving signals from a specified temperature distribution for
hyperthermia. IEEE Trans. Biomed. Eng, 39(8):825-835, 1992.

[30] J. W. Goodman. Introduction to Fourier optics. McGraw-Hill, New York, second
edition, 1996.

[31] K. B. Ocheltree and L. A. Frizzell. Sound field calculation for rectangular sources.
IEEE Trans. Ultrason. Ferroelect. Freq. Contr., 36(2):242—248, 1989.

[32] J. A. Jensen and N. B. Svendsen. Calculation of pressure fields from arbitrarily
shaped, apodized, and excited ultrasound transducers. IEEE Trans. Ultrason.
Ferroelect. Freq. Contr., 39(1):262—267, 1992.

[33] P. R. Stepanishen. Transient radiation from pistons in an infinite planar baffle.
J. Acoust. Soc. Am., 49(5):1629—1638, 1971.

[34] D. A. Hutchins, H. D. Mair, P. A. Puhach, and A. J. Osei. Continuous-wave
pressure fields of ultrasonic transducers. J. Acoust. Soc. Am., 80(1):1—12, 1986.

[35] P. R. Stepanishen and K. C. Benjamin. Forward and backward projection of
acoustic fields using FFT methods. J. Acoust. Soc. Am., 71(4):803—812, 1982.

157

 

[36] G. T. Clement and K. Hynynen. Field characterization of therapeutic ultrasound
phased arrays through forward and backward planar projection. J. Acoust. Soc.
Am., 108(1):441—446, 2000.

[37] P. Wu, R. Kazys, and T. Stepinski. Optimal selection of parameters for the
angular spectrum approach to numerically evaluate acoustic fields. J. Acoust.
Soc. Am., 101(1):125—134, 1997.

[38] P. Wu, R. Kazys, and T. Stepinski. Analysis of the numerically implemented
angular spectrum approach based on the evaluation of two-dimensional acoustic
fields. part I. error due to the discrete fourier transform and discretization. J.
Acoust. Soc. Am., 99(3):1339—1348, 1996.

[39] P. Wu, R. Kazys, and T. Stepinski. Analysis of the numerically implemented
angular spectrum approach based on the evaluation of two-dimensional acoustic
fields. part II. characteristics as a function of angular range. J. Acoust. Soc. Am.,
99(3):1349—1359, 1996.

[40] D. P. Orofino and P. C. Pedersen. Efficient angular spectrum decomposition of
acoustic sources—part 1: Theory. IEEE Trans. Ultrason. Ferroelect. Freq. Contra,
40(3):238—249, 1993.

[41] D. P. Orofino and P. C. Pedersen. Efficient angular spectrum decomposition of
acoustic sources—part II: Results. IEEE Trans. Ultrason. Ferroelect. Freq. Contr.,
40(3):250—257, 1993.

[42] P. T. Christopher and K. J. Parker. New approaches to the linear propagation
of acoustic fields. J. Acoust. Soc. Am., 90(1):507—521, 1991.

[43] D. L. Liu and R. C. Wagg. Propagation and backpropagation for ultrasonic
wavefront design. IEEE Trans. Ultrason. Ferroelect. Freq. Contr., 44(1):1—13,
1997.

[44] R. J. Zemp and J. T. Tavakkoli. Modeling of nonlinear ultrasound propagation
in tissue from array transducers. J. Acoust. Soc. Am., 113(1):139—152, 2003.

[45] E. G. Williams and J. D. Maynard. Numerical evaluation of the Rayleigh integral
for planar radiators using FFT. J. Acoust. Soc. Am., 72(6):2020—2030, 1982.

[46] E. G. Williams. Fourier acoustics: Sound radiation and nearfield acoustical
holography. Academic press, London, 1999.

[47] D. L. Liu and R. C. Waag. Correction of ultrasonic wavefront distortion using
backpropagation and a reference waveform method for time-shift compensation.
J. Acoust. Soc. Am., 96(2 Pt 1):649—60, 1994.

158

 

[48] K. Y. Saleh and N. B. Smith. Two-dimensional ultrasound phased array de-
sign for tissue ablation for treatment of benign prostatic hyperplasia. Int. J.
Hyperthenn, 20(1):7—31, 2004.

[49] J. S. Tan, L. A. Frizzell, and N. Sanghvi. Ultrasound phased arrays for prostate
treatment. J. Acoust. Soc. Am., 109(6):3055—3064, 2001.

[50] R. J. McGough, D. Cindric, and T. V. Samulski. Shape calibration of a confor-
mal ultrasound therapy array. IEEE Trans. Ultrason. Ferroelect. Freq. Contr,
48(2):494—505, 2001. ‘

[51] B. Piwakowski and K. Sbai. A new approach to calculate the field radiated from
arbitrarily structured transducer arrays. IEEE Trans. Ultrason. Ferroelect. Freq.
Contr., 46(2):422-—440, 1999.

[52] B. Piwakowski and K. Sbai. A new calculation procedure for spatial impulse
responses in ultrasound. J. Acoust. Soc. Am., 105(6):3266—3274, 1999.

[53] X. Zeng and R. J. McGough. Fast pressure field calculations applied to large
spherical ultrasound phased arrays designed for thermal therapy. Proceedings of
SPIE, 5698, 2005.

[54] M. S. Ibbini and C. A. Cain. A field conjugation method for direct synthesis of
hyperthermia phased-array heating patterns. IEEE Trans. Ultrason. Ferroelect.
Freq. Contr, 36(2):3—9, 1989.

[55] R. J. McGough, H. Wang, E. S. Ebbini, and C. A. Cain. Mode scanning:
heating pattern synthesis with ultrasound phased arrays. Int. J. Hyperthermia,
10(3):433—442, 1994.

[56] P. Godden, G. ter Haar, and I. Rivens. Numerical modelling of high intensity
focused ultrasound phased arrays by an angular spectrum decomposition using
highly-oscillatory quadrature. International Symposium on Therapeutic Ultra-
sound, pages 36—42, 2007.

[57] T. D. Mast, L. P. Souriau, D.-L. D. Liu, M. Tabei, A. I. Nachman, and R. C.
Waag. A k-space method for large-scale models of wave propagation in tissue.
IEEE Trans. Ultrason. Ferroelect. Freq. Contr., 48(2):341—354, 2001.

[58] T. D. Mast, L. M. Hinkelman, M. J. Orr, V. W. Sparrow, and R. C. Waag. Sim-
ulation of ultrasonic pulse propagation through the abdominal wall. J. Acoust.
Soc. Am., 102(2):1177—1190, 1997.

[59] T. D. Mast, L. M. Hinkelman, M. J. Orr, and R. C. Waag. The effect of abdominal
wall morphology on ultrasonic pulse distortion. part ii. simulations. J. Acoust.
Soc. Am., 104(6):3651?664, 1998.

159

[60] T. D. Mast. Two— and three-dimensional simulations of ultrasonic propagation
through human breast tissue. Acoust. Res. Lett. Online, 3(2):53?8, 2002.

[61] M. Tabei, T. D. Mast, and R. C. Waag. Simulation of ultrasonic focus aberration
and correction through human tissue. J. Acoust. Soc. Am., 113(2):1166—1176,
2003.

[62] T. D. Mast, L. M. Hinkelman, L. A. Metlay, M. J. Orr, and R. C. Waag. Simu-
lation of ultrasonic pulse propagation, distortion, and attenuation in the human
chest wall. J. Acoust. Soc. Am., 106(6):3665—3677, 1999.

[63] X. Fan and K. Hynynen. The effect of wave reflection and refraction at soft

tissue interfaces during ultrasound hyperthermia treatments. J. Acoust. Soc.
Am., 91(3):1727—1736, 1992.

[64] X. Fan and K. Hynynen. The effects of curved tissue layers on the power de-
position patterns of therapeutic ultrasound beams. Med. Phys, 21(1):25—34,
1994.

[65] G. T. Clement and K. Hynynen. Forward planar projection through layered
media. IEEE Trans Ultrason Ferroelectr Freq Control, 50(12):1689-98, 2003.

[66] M. E. Haran. Distortion of finite amplitude ultrasound in lossy media. J. Acoust.
Soc. Am., 73(3):774—779, 1983.

[67] W. Swindell. A theoretical study of nonlinear effects with focused ultrasound in
tissues: an ”Acoustic Bragg peak”. Ultrasound Med. Biol, 11(1):121—30, 1985.

[68] D. H. Trivett and A. L. Van Buren. Propagation of plane, cylindrical and spher-
ical finite amplitude waves. J. Acoust. Soc. Am., 69(4):943—949, 1981.

[69] S. I. Aanonsen and T. Barkve. Distortion and harmonic generation in the
nearfield of a flute amplitude sound beam. J. Acoust. Soc. Am., 75(3):749—768,
1984.

[70] A. Korpel. Frequency approach to nonlinear dispersive waves. J. Acoust. Soc.
Am., 67(6):1954—1958, 1980.

[71] G. Du and M. A. Breazeale. Harmonic distortion of a finite amplitude gaussian
beam in a fluid. J. Acoust. Soc. Am., 80(1):212—216, 1986.

[72] T. S. Hart and M. F. Hamilton. Nonlinear effects in focused sound beams. J.
Acoust. Soc. Am., 84(4):1488—1496, 1988.

[73] Y. S. Lee and M. F. Hamilton. Time-domain modeling of pulsed finite-amplitude
sound beams. J. Acoust. Soc. Am., 97(2):906—917, 1995.

160

[74] G. Wojcik, J. Mould Jr., N. Abboud, M. Ostromogilsky, and D. Vaughan.
Nonlinear modeling of therapeutic ultrasound. In Proceedings of
the IEEE Ultrasonics Symposium, volume 2, pages 1617—1622, 1995.
http: //control.ee.ethz.ch/ joloef/yalmip.php.

[75] P. T. Christopher and K. J. Parker. New approaches to nonlinear diffractive field
propagation. J. Acoust. Soc. Am., 90(1):488—499, 1991.

[76] T. Christopher. Modeling acoustic field propagation for medical devices. PhD
thesis, University of Rochester, 1993.

[77] C. J. Vecchio and P. A. Lewin. Finite amplitude acoustic propagation modeling
using the extended angular spectrum method. J. Acoust. Soc. Am., 95(5):2399—
2408, 1994.

[78] C. J. Vecchio. Finite amplitude acoustic propagation modeling using the emtended
angular spectrum method. PhD thesis, Drexel University, 1992.

[79] E. L. Carstensen, W. K. Law, N. D. McKay, and T. G. Muir. Demonstration of
nonlinear acoustical effects at biomedical frequencies and intensities. Ultrasound
Med. Biol, 6(4):359—68, 1980.

[80] D. Dalecki, E. L. Carstensen, and K. J. Parker. Absorption of finite amplitude
focused ultrasound. J. Acoust. Soc. Am., 89(5):2435—2447, 1991.

[81] D. R. Bacon and E. L. Carstensen. Increased heating by diagnostic ultrasound
due to nonlinear propagation. J. Acoust. Soc. Am, 88(1):26—34, 1990.

[82] H. C. Starritt, M. A. Perkins, F. A. Duck, and V. F. Humphrey. Evidence
for ultrasonic finite—amplitude distortion in muscle using medical equipment. J.
Acoust. Soc. Am., 77(1):302—6, 1985.

[83] K. D. Wallace. Characterization of the nonlinear propagation of difiracting, finite
amplitude ultrasonic fields. PhD thesis, Washington University, St Louis, 2001.

[84] A. C. Baker, K. Anastasiadis, and V. F. Humphrey. The nonlinear pressure
field of a plane circular piston: Theory and experiment. J. Acoust. Soc. Am.,
84(4):1483—1487, 1988.

[85] T. J. Cavicchi and W. D. O’Brien Jr. Heat generated by ultrasound in an
absorbing medium. J. Acoust. Soc. Am., 76(4):1244—1245, 1984.

[86] H. A. Pennes. Analysis of tissue and arterial blood temperature in the resting
human forearm. J. Appl. Physiol, 1(2):93, 1948.

161

[87] K. B. Ocheltree and L. A. Frizzell. Determination of power deposition patterns
for localized hyperthermia: a transient analysis. Int. J. Hypertherm., 4(3):281—
296, 1988.

[88] K. B. Ocheltree and L. A. Frizzell. Determination of power deposition patterns for
localized hyperthermia: a steady state analysis. Int. J. Hypertherm, 3(3):269—
279, 1987.

[89] E. G. Moros, R. B. Roemer, and K. Hynynen. Pre-focal plane high-temperature
regions induced by scanning focused ultrasound beams. Int. J. Hypertherm,
6(2):351—366, 1990.

[90] S. A. Sapareto and W. C. Dewey. Thermal dose determination in cancer therapy.
Int. J. Radiation Oncology Biol. Phys, 10(6):787—796, 1984.

[91] W. L. Lin, R. B. Roemer, Moros E. G., and K. Hynynen. Optimization of
temperature distributions in scanned, focused ultrasound hyperthermia. Int. J.
Hyperthermia, 8(1):61—78, 1992.

[92] E. S. Ebbini and C. A. Cain. Optimization of the intensity gain of multiple-focus
phased-array heating patterns. Int. J. Hyperthermia, 7(6):953—73, 1991.

[93] B. Guo, L. Xu, and J. Li. Time reversal based microwave hyperthermia treatment
of breast cancer. Microwave and Optical Technology Letters, 47(4):335—338, 2005.

[94] B. Guo and J. Li. Waveform diversity based ultrasound system for hyperthermia
treatment of breast cancer. IEEE Trans. Biomed. Eng, 55(2):822—826, 2008.

[95] J. F. Sturm. Using SeDuMi 1.02, a matlab toolbox for optimization over sym—
metric cones. Optimization Methods and Software, 11-12:625?53, 1999.

[96] K. C. Toh, M. J. Todd, and R. H. 'I‘utuncu. SDPT3 - a matlab software package
for semidefinite programming. Optimization Methods and Software, 11:545—581,
1999.

[97] R. H. Tutuncu, K. C. Toh, and M. J. Todd. Solving semidefinite—quadratic—linear
programs using SDPT3. Mathematical Programming Ser. B, 95:189—217, 2003.

[98] J. Lofberg. Yalmip : A toolbox for modeling and optimization in MAT-

LAB. In Proceedings of the CACSD Conference, Taipei, Taiwan, 2004.
http: //control.ee.ethz.ch/ joloef/yalmip.php.

[99] M. Grant and S. Boyd. CVX: Matlab software for disciplined convex program—
ming, 2008. http://stanford.edu/ boyd/cvx.

162