**Compensation of Ultrasound Attenuation in Photoacoustic Imaging**

P. Burgholzer1,2, H. Roitner1,2, J. Bauer-Marschallinger1,2, H. Grün2, T. Berer1,2 and G. Paltauf3 *1Christian Doppler Laboratory for Photoacoustic Imaging and Laser Ultrasonics, 2Research Center for Non Destructive Testing (RECENDT), 3Institute of Physics, Karl-Franzens-University Graz Austria* 

## **1. Introduction**

190 Acoustic Waves – From Microdevices to Helioseismology

Singh, A., Houser, D. R. & Vijayakar, S. (1996). Early detection of gear pitting, *American Society of Mechanical Engineers*, Vol. 88, (1996) pp.673-678, ISSN 15214613. Siores, E. & Negro, A.A. (1997). Condition monitoring of a gear box using acoustic emission

Tan, C.K. & Mba, D. (2005). Limitation of acoustic emission for identifying seeded defects in

Tan, C.K., Irving, P. & Mba, D. (2007). A comparative experimental study on the diagnostic

Wu, J.D., Hsu, C.C. & Wu, G.Z. (2009). Fault gear identification and classification using

Yang, Y., Yu, D. & Cheng, J. (2007). A fault diagnosis approach for roller bearing based on

*Applications*, Vol. 36, No. 3, (April 2009), pp.6244-6255, ISSN 0957-4174. Wu, J.D. & Chen, J.C. (2006). Continuous wavelet transform technique for fault signal

00255327

0263-2231.

28, ISSN 0195-9298.

(January 2007), pp.208-233, ISSN 0888-3270

(June 2006), pp.304-311, ISSN 0963-8695.

testing, *Material Evaluation*, Vol. 55, No. 2, (February 1997), pp.183-187, ISSN

gearboxes, *Journal of Non-Destructive Evaluation*, Vol. 24, No. 1, (March 2005), pp.11-

and prognostic capabilities of acoustics emission, vibration and spectrometric oil analysis for spur gears, *Mechanical Systems and Signal Processing*, Vol. 21, No. 1,

discrete wavelet transform and adaptive neuro-fuzzy inference, *Expert Systems with* 

diagnosis of internal combustions engines, *NDT&E International*, Vol. 39, No. 4,

IMF envelope spectrum and SVM, *Measurement: Journal of the International Measurement Confederation*, Vol. 40, No. 9-10, (November 2007), pp.943-950, ISSN Photoacoustic imaging is a non-destructive method to obtain information about the distribution of optically absorbing structures inside a semitransparent medium. It is based on thermoelastic generation of ultrasonic waves by the absorption of a short laser pulse inside the sample. From the ultrasonic waves measured outside the object, the interior distribution of absorbed energy is reconstructed. The ultrasonic waves, which transport information from the interior to the surface of the sample, are scattered or absorbed to a certain extent by dissipative processes. The scope of this work is to quantify the information loss which is equal to the entropy production during these dissipative processes and thereby to give a principle limit for the spatial resolution which can be gained in photoacoustic imaging. This theoretical limit is compared to experimental data. In this book chapter stateof-the-art methods for modeling ultrasonic wave propagation in the case of attenuating media are described. From these models strategies for compensating ultrasound attenuation are derived which may be combined with well-known reconstruction algorithms from the non-attenuating case for photoacoustic imaging.

Section 2 gives a short description of photoacoustic imaging, especially photoacoustic tomography, and the available image reconstruction algorithms to reconstruct the interior structure from the detected ultrasound signal at the sample surface. Beside small point-like detectors also large detectors, so called integrating detectors are used for photoacoustic tomography. The latter ones require different image reconstruction algorithms. Spatial resolution is an essential issue for any imaging method. Therefore we describe the influencing factors of the resolution in photoacoustic tomography.

Section 3 is dedicated to acoustic attenuation. The spatial resolution in photoacoustic imaging is limited by the acoustic bandwidth. To resolve small objects shorter wavelengths with higher frequencies are necessary. For such high frequencies, however, the acoustic attenuation increases. This effect is usually ignored in photoacoustic image reconstruction but as small objects or structures generate high frequency components it limits the minimum detectable size, hence the resolution. Several models for acoustic attenuation, especially used for ultrasound propagation in biological tissue, are compared with experimental data.

Compensation of Ultrasound Attenuation in Photoacoustic Imaging 193

ultrasound wave arriving from all different locations of the source. A map or "image" of the photo-generated pressure distribution in the sample can be made by collecting the ultrasound at many different locations and processing it using a suitable algorithm e.g. by a

Only if the pulse is short enough, thermal expansion causes a pressure rise proportional to the locally absorbed energy density. Short enough means that the pressure wave does not "run out" of the smallest structure which should be resolved in the photoacoustic image during the pulse time. This so called "stress confinement" is therefore fulfilled if the sound velocity multiplied by the pulse time of the laser is small compared to the spatial resolution one wants to achieve in imaging. Another constraint is the "thermal confinement" which is fulfilled if the heat induced in a structure by the absorbed laser pulse does not diffuse out of this structure during the time of the laser pulse. As heat diffusion is usually slower than the propagation of sound the thermal confinement is fulfilled if stress confinement is fulfilled.

Fig. 1. Photoacoustic Imaging – the spatial resolution is determined by excitation,

resolution

infrared light.

propagation, and detection of the acoustic wave. (a) Thermoelastic generation of acoustic wave by laser light (arrows indicate scattered photons): excited pressure is proportional to the absorbed optical energy density, if laser pulse is short enough to satisfy thermal and stress confinement. (b) Propagation of ultrasonic wave to sample surface: frequency dependent acoustic attenuation causes entropy production and therefore a loss of

information. (c) Detection of ultrasonic wave: bandwidth and size of detector limits spatial

Any photons, either unscattered or scattered (see arrows in Fig. 1), contribute to the absorbed energy as long as the photon excitation is relaxed thermally. Therefore PAT visualizes the product of the optical absorption distribution and the local light fluence. Using a Nd:YAG laser and an optical parametric oscillator (OPO) light pulses from the infrared to the visible regime can be selected with a repetition rate from 10 Hz up to 100 Hz. Some high speed PAT systems can even go up to 1000 Hz. The pulse duration in the nanosecond range enables a theoretical resolution of several microns in tissue (sound velocity similar to water at approx. 1500 m/s). For biomedical applications the light energy should not exceed 20 mJ/cm2 in the visible spectral range and 100 mJ/cm2 in the near

filtered backprojection algorithm or by a time reversal algorithm.

Section 4 describes two different attempts to compensate this acoustic attenuation: either to include the compensation directly in the image reconstruction, e. g. in a modified time reversal method, or to calculate first the acoustic signal without attenuation from the measured attenuated signal and then perform the conventional photoacoustic image reconstruction. As any compensation of acoustic attenuation is mathematically an ill-posed problem both methods need regularization to prevent small measurement fluctuations from growing infinitely high in the reconstructed image. The possible degree of compensation depends on the size of these fluctuations. On the other hand acoustic attenuation is a dissipative process that causes entropy production equal to a loss of information, which cannot be compensated by any compensation algorithm. Therefore one can use the entropy production caused by acoustic attenuation to determine the minimal fluctuations in the measurement data, which turn out to be equal to thermal fluctuations. In statistical physics this fact is well known as the fluctuation-dissipation theorem, but the information theoretical background as a starting point to derive this theorem was not mentioned before in the literature.

Section 5 uses stochastic processes to understand theoretically how information can be lost and its connection to entropy production. Therefore, the measured pressure signal is treated as a random variable with a certain mean value as a function of time and certain fluctuations around that value. First for the simple model of a damped harmonic oscillator, it is shown how information is lost during a dissipative process and to what extent we can reconstruct the original information after some time. Then attenuated acoustic waves can be treated in a similar way: the spatial Fourier transform of the pressure wave can be described by a similar stochastic process as the damped harmonic oscillator – only in a higher dimension. Each wave vector is represented by a damped oscillator of different frequency.

Thinking about acoustic attenuation as a stochastic process helps to understand how entropy production and loss of information "work" on a microscopic scale. Beside a better theoretical insight the stochastic view on the acoustic wave answers a very important question: which is the best compensation method and the corresponding practicable spatial resolution in photoacoustic imaging? This question can be answered without taking fluctuations on a microscopic scale into account: the entropy production, which can be calculated from macroscopic mean values, is set equal to the information loss.

## **2. Photoacoustic imaging**

In 1880, Alexander Graham Bell discovered that pulsed light striking a solid substrate can produce a sound wave, a phenomenon called the photoacoustic effect (Bell, 1880). Practical imaging methods based on this effect have been developed and reported the last decade (Xu & Wang, 2006). Today, photoacoustic imaging, which is also referred to as optoacoustic imaging or, when using microwaves instead of light for excitation, as thermoacoustic imaging, is attracting intense interest for cross-sectional or three-dimensional imaging in biomedicine.

In photoacoustic imaging, short laser pulses are fired at a sample and the absorbed energy causes local heating (Fig. 1). This heating causes thermoelastic expansion and thereby generation of broadband elastic pressure waves (ultrasound) which can be detected outside the sample, for example by a piezoelectric device or by an optical detector. Two methods are used for photoacoustic imaging: photoacoustic microscopy uses focused ultrasonic detectors and the sample is imaged by scanning the focus through the sample. In photoacoustic tomography (PAT) an unfocused detector is used which detects the pressure from the

Section 4 describes two different attempts to compensate this acoustic attenuation: either to include the compensation directly in the image reconstruction, e. g. in a modified time reversal method, or to calculate first the acoustic signal without attenuation from the measured attenuated signal and then perform the conventional photoacoustic image reconstruction. As any compensation of acoustic attenuation is mathematically an ill-posed problem both methods need regularization to prevent small measurement fluctuations from growing infinitely high in the reconstructed image. The possible degree of compensation depends on the size of these fluctuations. On the other hand acoustic attenuation is a dissipative process that causes entropy production equal to a loss of information, which cannot be compensated by any compensation algorithm. Therefore one can use the entropy production caused by acoustic attenuation to determine the minimal fluctuations in the measurement data, which turn out to be equal to thermal fluctuations. In statistical physics this fact is well known as the fluctuation-dissipation theorem, but the information theoretical background as a starting point to derive this theorem was not mentioned before

Section 5 uses stochastic processes to understand theoretically how information can be lost and its connection to entropy production. Therefore, the measured pressure signal is treated as a random variable with a certain mean value as a function of time and certain fluctuations around that value. First for the simple model of a damped harmonic oscillator, it is shown how information is lost during a dissipative process and to what extent we can reconstruct the original information after some time. Then attenuated acoustic waves can be treated in a similar way: the spatial Fourier transform of the pressure wave can be described by a similar stochastic process as the damped harmonic oscillator – only in a higher dimension. Each

Thinking about acoustic attenuation as a stochastic process helps to understand how entropy production and loss of information "work" on a microscopic scale. Beside a better theoretical insight the stochastic view on the acoustic wave answers a very important question: which is the best compensation method and the corresponding practicable spatial resolution in photoacoustic imaging? This question can be answered without taking fluctuations on a microscopic scale into account: the entropy production, which can be

In 1880, Alexander Graham Bell discovered that pulsed light striking a solid substrate can produce a sound wave, a phenomenon called the photoacoustic effect (Bell, 1880). Practical imaging methods based on this effect have been developed and reported the last decade (Xu & Wang, 2006). Today, photoacoustic imaging, which is also referred to as optoacoustic imaging or, when using microwaves instead of light for excitation, as thermoacoustic imaging, is attracting intense interest for cross-sectional or three-dimensional imaging in biomedicine. In photoacoustic imaging, short laser pulses are fired at a sample and the absorbed energy causes local heating (Fig. 1). This heating causes thermoelastic expansion and thereby generation of broadband elastic pressure waves (ultrasound) which can be detected outside the sample, for example by a piezoelectric device or by an optical detector. Two methods are used for photoacoustic imaging: photoacoustic microscopy uses focused ultrasonic detectors and the sample is imaged by scanning the focus through the sample. In photoacoustic tomography (PAT) an unfocused detector is used which detects the pressure from the

wave vector is represented by a damped oscillator of different frequency.

calculated from macroscopic mean values, is set equal to the information loss.

in the literature.

**2. Photoacoustic imaging** 

ultrasound wave arriving from all different locations of the source. A map or "image" of the photo-generated pressure distribution in the sample can be made by collecting the ultrasound at many different locations and processing it using a suitable algorithm e.g. by a filtered backprojection algorithm or by a time reversal algorithm.

Only if the pulse is short enough, thermal expansion causes a pressure rise proportional to the locally absorbed energy density. Short enough means that the pressure wave does not "run out" of the smallest structure which should be resolved in the photoacoustic image during the pulse time. This so called "stress confinement" is therefore fulfilled if the sound velocity multiplied by the pulse time of the laser is small compared to the spatial resolution one wants to achieve in imaging. Another constraint is the "thermal confinement" which is fulfilled if the heat induced in a structure by the absorbed laser pulse does not diffuse out of this structure during the time of the laser pulse. As heat diffusion is usually slower than the propagation of sound the thermal confinement is fulfilled if stress confinement is fulfilled.

Fig. 1. Photoacoustic Imaging – the spatial resolution is determined by excitation, propagation, and detection of the acoustic wave. (a) Thermoelastic generation of acoustic wave by laser light (arrows indicate scattered photons): excited pressure is proportional to the absorbed optical energy density, if laser pulse is short enough to satisfy thermal and stress confinement. (b) Propagation of ultrasonic wave to sample surface: frequency dependent acoustic attenuation causes entropy production and therefore a loss of information. (c) Detection of ultrasonic wave: bandwidth and size of detector limits spatial resolution

Any photons, either unscattered or scattered (see arrows in Fig. 1), contribute to the absorbed energy as long as the photon excitation is relaxed thermally. Therefore PAT visualizes the product of the optical absorption distribution and the local light fluence. Using a Nd:YAG laser and an optical parametric oscillator (OPO) light pulses from the infrared to the visible regime can be selected with a repetition rate from 10 Hz up to 100 Hz. Some high speed PAT systems can even go up to 1000 Hz. The pulse duration in the nanosecond range enables a theoretical resolution of several microns in tissue (sound velocity similar to water at approx. 1500 m/s). For biomedical applications the light energy should not exceed 20 mJ/cm2 in the visible spectral range and 100 mJ/cm2 in the near infrared light.

Compensation of Ultrasound Attenuation in Photoacoustic Imaging 195

cylindrical surface with the radius *c t*⋅ where *c* is the speed of sound in the medium and *t* the time. Thus, integrating line detectors arranged in an array around the sample, e.g. in a circle, measure projection images of the object in a first measurement and reconstruction step. By rotating the sample and measuring such projection images from different angles it is possible to reconstruct a 3D image (Fig. 2). Three dimensional image reconstruction from a set of projection images requires only the application of the inverse Radon transform. Therefore 3D imaging using integrating line detectors is not computationally intensive compared with other algorithms which reconstruct a 3D image from a set of signals

Fig. 2. Line detectors around a sample rotated on one axis. Either one line detector is

Since the first measurement results from (Burgholzer et al, 2005) the integrating line detector was further improved in sensitivity and spatial resolution. Several types of such line detectors were implemented. One premature approach was a line made of PVDF (piezo foil) which provides high sensitivity but with the drawback of directivity. The next consequential step was an optical line detector realized by an interferometer. A laser beam that is part of an interferometer measures variations of the refractive index induced by the acoustic pressure (elasto-optic effect) (Paltauf et al., 2006). Optical detectors offer a broadband characteristic and due to the circular shape of a laser beam or light guiding fiber there is no such directivity like when using a film of piezo foil. Two main approaches can be distinguished: free-beam implementations and the use of fiber-based interferometers. Independent of the realization different types of the interferometer can be used, e.g. a Mach-Zehnder or a Fabry-Perot interferometer which is in general more sensitive than the first one. (Paltauf et al., 2006) presented measurements using a free-beam Mach-Zehnder interferometer. They used a focused laser beam as line detector. When placing the object next to the beam waist of the focused laser beam the best spatial resolution due to the

(Grün et al, 2010) implemented fiber-based line detectors. The advantages of fiber-based line detectors are the easy handling and the small and constant beam diameter in the fiber. A small diameter of the laser beam is necessary for a good spatial resolution. The smaller the

scanning around the object or a detector array is used

smallest beam diameter could be achieved (Paltauf et al., 2008).

acquired from point like detectors.

If acoustic attenuation and shear waves (in liquid and soft tissue) are neglected, the acoustic pressure *p* as a function of time and space obeys the equation

$$
\Delta p(\mathbf{r}, t) - \frac{1}{c^2} \frac{\partial^2}{\partial t^2} p(\mathbf{r}, t) = -\frac{\beta}{C\_p} \frac{\partial}{\partial t} \mathcal{H}(\mathbf{r}, t) \tag{1}
$$

where Δ is the three-dimensional Laplace operator, *c* the sound velocity, β the thermal expansion coefficient, *Cp* the specific heat capacity and *H t* (,) **r** the deposited energy per time and volume ("heating function") caused by the absorption of the electromagnetic radiation in the sample.

For short electromagnetic pulses *Ht A t* (,) () () **r r** = ⋅δ , where *A*( )**r** is the energy density of the absorbed electromagnetic radiation and δ ( )*t* the Dirac delta function. Then the acoustic pressure *p*(,) **r** *t* solves the homogeneous scalar wave equation with the initial conditions ( ) () () <sup>2</sup> <sup>0</sup> ,0 ( ) *<sup>p</sup> p p cC A A* **rr r r** = = ⋅ ≡Γ⋅ β and ∂∂ = / ( ,0) 0 *t p* **r** . The initial pressure 0 *p* at time *t* = 0 is therefore directly proportional to the absorbed energy density *A* with the dimensionless constant Γ , the Grüneisen coefficient.

As shown in Fig. 1 (c) bandwidth and size of the detectors for collecting the ultrasound signals are important for the resolution of this imaging modality. A photoacoustically generated ultrasound signal is a broadband signal and contains frequencies in the range from kilohertz up to a few megahertz. Conventional point like piezo elements such as known from arrays for medical ultrasonic imaging have their maximum sensitivity close to their center frequency and can detect frequencies only within a certain bandwidth around this frequency. Therefore high frequencies are not detected which correspond to small structures and are necessary for image reconstruction with a high spatial resolution. Other approaches are necessary for high resolution photoacoustic imaging. A hydrophone could be one solution (Wang, 2008), or the utilization of optical point like detection as demonstrated e.g. by (Zhang et al., 2009) or (Berer et al., 2010).

Point like detectors show a limit in achievable resolution by their size. The smaller the point detector the better is the spatial resolution. Unfortunately thermal and other fluctuations increase for a smaller detector, which results in a reduction of resolution. A totally different approach is the use of so called integrating detectors which are at least in one dimension larger than the object. This way the drawback of finite dimensions of point like detectors can be overcome. Such an integrating detector for photoacoustic imaging was introduced by (Haltmeier et al, 2004). They showed the mathematical proof of integrating area and line detectors and introduced new reconstruction methods which are necessary when using such a detector. First measurements using an integrating detector were shown by (Burgholzer et al., 2005). The first integrating detector was an area detector which was bigger than the object in two dimensions. For sufficient data for 3D image reconstruction the area detector had to be scanned around the object tangential to the surface of a sphere. This detector movement is difficult to realize. Hence the idea of the integrating line detector was developed. A fragmentation of the area detector into an array of line detectors results in an easier setup with only one rotation axis for the object and a linear motion of the integrating line detector (Fig. 2).

An integrating line detector is a line which has at least a length 8 \* *D* , where *D* is the diameter of a circle enclosing the sample and tangentially touching the line detector (Haltmeier et al., 2004). The line detector integrates the pressure along the line on a

If acoustic attenuation and shear waves (in liquid and soft tissue) are neglected, the acoustic

( ) ( ) ( ) <sup>2</sup>

∂ ∂

*p*

**rr r** (1)

β

, where *A*( )**r** is the energy density of the

( )*t* the Dirac delta function. Then the acoustic

and ∂∂ = / ( ,0) 0 *t p* **r** . The initial pressure 0 *p* at time

the thermal

β

<sup>1</sup> ,, ,

*p t pt t ct Ct* ∂ ∂

expansion coefficient, *Cp* the specific heat capacity and *H t* (,) **r** the deposited energy per time and volume ("heating function") caused by the absorption of the electromagnetic

δ

pressure *p*(,) **r** *t* solves the homogeneous scalar wave equation with the initial conditions

*t* = 0 is therefore directly proportional to the absorbed energy density *A* with the

As shown in Fig. 1 (c) bandwidth and size of the detectors for collecting the ultrasound signals are important for the resolution of this imaging modality. A photoacoustically generated ultrasound signal is a broadband signal and contains frequencies in the range from kilohertz up to a few megahertz. Conventional point like piezo elements such as known from arrays for medical ultrasonic imaging have their maximum sensitivity close to their center frequency and can detect frequencies only within a certain bandwidth around this frequency. Therefore high frequencies are not detected which correspond to small structures and are necessary for image reconstruction with a high spatial resolution. Other approaches are necessary for high resolution photoacoustic imaging. A hydrophone could be one solution (Wang, 2008), or the utilization of optical point like detection as

Point like detectors show a limit in achievable resolution by their size. The smaller the point detector the better is the spatial resolution. Unfortunately thermal and other fluctuations increase for a smaller detector, which results in a reduction of resolution. A totally different approach is the use of so called integrating detectors which are at least in one dimension larger than the object. This way the drawback of finite dimensions of point like detectors can be overcome. Such an integrating detector for photoacoustic imaging was introduced by (Haltmeier et al, 2004). They showed the mathematical proof of integrating area and line detectors and introduced new reconstruction methods which are necessary when using such a detector. First measurements using an integrating detector were shown by (Burgholzer et al., 2005). The first integrating detector was an area detector which was bigger than the object in two dimensions. For sufficient data for 3D image reconstruction the area detector had to be scanned around the object tangential to the surface of a sphere. This detector movement is difficult to realize. Hence the idea of the integrating line detector was developed. A fragmentation of the area detector into an array of line detectors results in an easier setup with only one rotation axis for the object and a linear motion of the integrating

An integrating line detector is a line which has at least a length 8 \* *D* , where *D* is the diameter of a circle enclosing the sample and tangentially touching the line detector (Haltmeier et al., 2004). The line detector integrates the pressure along the line on a

Δ − =− Η

δ

2 2

where Δ is the three-dimensional Laplace operator, *c* the sound velocity,

pressure *p* as a function of time and space obeys the equation

For short electromagnetic pulses *Ht A t* (,) () () **r r** = ⋅

dimensionless constant Γ , the Grüneisen coefficient.

demonstrated e.g. by (Zhang et al., 2009) or (Berer et al., 2010).

absorbed electromagnetic radiation and

( ) () () <sup>2</sup> <sup>0</sup> ,0 ( ) *<sup>p</sup> p p cC A A* **rr r r** = = ⋅ ≡Γ⋅ β

radiation in the sample.

line detector (Fig. 2).

cylindrical surface with the radius *c t*⋅ where *c* is the speed of sound in the medium and *t* the time. Thus, integrating line detectors arranged in an array around the sample, e.g. in a circle, measure projection images of the object in a first measurement and reconstruction step. By rotating the sample and measuring such projection images from different angles it is possible to reconstruct a 3D image (Fig. 2). Three dimensional image reconstruction from a set of projection images requires only the application of the inverse Radon transform. Therefore 3D imaging using integrating line detectors is not computationally intensive compared with other algorithms which reconstruct a 3D image from a set of signals acquired from point like detectors.

Fig. 2. Line detectors around a sample rotated on one axis. Either one line detector is scanning around the object or a detector array is used

Since the first measurement results from (Burgholzer et al, 2005) the integrating line detector was further improved in sensitivity and spatial resolution. Several types of such line detectors were implemented. One premature approach was a line made of PVDF (piezo foil) which provides high sensitivity but with the drawback of directivity. The next consequential step was an optical line detector realized by an interferometer. A laser beam that is part of an interferometer measures variations of the refractive index induced by the acoustic pressure (elasto-optic effect) (Paltauf et al., 2006). Optical detectors offer a broadband characteristic and due to the circular shape of a laser beam or light guiding fiber there is no such directivity like when using a film of piezo foil. Two main approaches can be distinguished: free-beam implementations and the use of fiber-based interferometers. Independent of the realization different types of the interferometer can be used, e.g. a Mach-Zehnder or a Fabry-Perot interferometer which is in general more sensitive than the first one.

(Paltauf et al., 2006) presented measurements using a free-beam Mach-Zehnder interferometer. They used a focused laser beam as line detector. When placing the object next to the beam waist of the focused laser beam the best spatial resolution due to the smallest beam diameter could be achieved (Paltauf et al., 2008).

(Grün et al, 2010) implemented fiber-based line detectors. The advantages of fiber-based line detectors are the easy handling and the small and constant beam diameter in the fiber. A small diameter of the laser beam is necessary for a good spatial resolution. The smaller the

Compensation of Ultrasound Attenuation in Photoacoustic Imaging 197

Acoustic attenuation is an irreversible process and therefore the wave equation is not invariant to time reversal. Several important reasons for acoustic attenuation have been reported: viscosity, heat conduction, relaxation processes and chemical reactions. Stokes

> 2 2 2 <sup>1</sup> <sup>0</sup> *p p <sup>p</sup> ct t* τ

∂ ∂ Δ− + Δ =

under the assumption of adiabatic conditions (thus neglecting the loss due to heat conduction) (Stokes, 1845). This equation can also be found in (Shutilov, 1988) and can generally be derived from a relaxation behavior of pressure and density (Royer & Dieulesaint, 2000), where the density change follows the pressure change after a relaxation

known as the thermoviscous equation. Eq. (2) describes acoustic attenuation which is approximately proportional to the square of the frequency. Other wave equations can describe a more general power law frequency dependence of the attenuation of the form

( ) ( ) () ( ) <sup>2</sup>

<sup>∂</sup> Δ − +∗ =

<sup>1</sup> *p t p t Lt p t* , , ,0

Other models for acoustic wave propagation in acoustic media have been proposed also by

Fig. 3. Experimental set-up to measure broadband high frequency acoustic attenuation in tissue with ultrasound generated by short laser pulses (a) and by piezoelectric transducers

(b). In both cases a piezoelectric transducer receives the ultrasonic signals

2 2 0

∂

*c t*

the wave equation in order to account for such attenuation behavior:

**Measurement and simulation of broadband acoustic attenuation:** 

(Nachman et al., 1990) and (Treeby & Cox, 2010).

= ⋅ (0 < y < 3). (Szabo, 1994) has suggested adding the loss term *L t*() ( ) ∗ *p* **r**,*t* to

is further expressed in terms of viscosity and specific heat, this equation is also

**rrr** (3)

∂ ∂ (2)

derived the scalar wave equation

time τ . If τ

αω

<sup>0</sup> ( ) *<sup>y</sup>*

 α ω

diameter of the detecting part the better is the spatial resolution. A typical single mode fiber for near infrared has a core diameter of 9 microns; single mode fibers for the visible range of detection wavelength have typically about 6 microns. Due to the constant diameter along the whole line detector this type of integrating line detector is dedicated for the imaging of big samples. After first implementations of a Mach-Zehnder and a Fabry-Perot interferometer in glass fibers now polymer fibers are used. Due to the much better impedance matching of polymer fibers to the surrounding water, their sensitivity is higher than for glass fibers, where approximately 2/3 of the incoming signal is reflected before reaching the core (Grün et al., 2010). Furthermore the Young's modulus is much lower in polymer fibers than in glass fibers for which reason the deformation of a polymer fiber is bigger than of a glass fiber applying the same pressure wave. As the strain optic coefficients are in the same order, this results in an enhanced change of refractive index and thus to higher signal amplitudes in the polymer fiber (Kiesel et al., 2007).

(Nuster et al, 2009) did a comparison of the different implementations of an integrating line detector. At this stage of development the free-beam Mach-Zehnder interferometer was the most sensitive integrating line detector. But these measurements showed some new approaches how the fiber-based line detector could be made more sensitive, e.g. by building up a Fabry-Perot interferometer in a single mode polymer fiber.

The next step after developing a sensitive line detector, no matter of which approach is the most sensitive one, is the creation of an array of many integrating line detectors, e.g. 200 detectors arranged in a curve around the sample. This way one could acquire all data for a projection image within one excitation pulse of the laser.

### **3. Acoustic attenuation**

The imaging resolution in photoacoustic imaging is limited by the acoustic bandwidth and therefore by the laser pulse duration as mentioned above, but also by the attenuation of the acoustic wave on its way to the sample surface, and finally by the bandwidth and size of the ultrasonic transducer. The acoustic attenuation can be substantial for high frequencies. This effect is usually ignored in reconstruction algorithms but can have a strong impact on the resolution of small objects or structures within objects. Stokes could already show in 1845 that for liquids with low viscosity, such as water, the acoustical absorption increases by the square of the frequency (Stokes, 1845). If we do not take acoustic attenuation for photoacoustic image reconstruction into account, especially small structures (corresponding to shorter wavelengths and therefore higher frequencies) appear blurred (La Riviere et al., 2006). To what extent this blurring can be compensated by regularization methods (as performed by (La Riviere et al., 2006)) and how much information is lost due the irreversibility of attenuation is investigated in this chapter.

For thin layers (1D), small cylinders (2D), and small spherical inclusions (3D) the effect of attenuation is simulated and experimental results for several types of tissue are given. For photoacoustic tomography a new description of attenuation seems to be useful: like for a standing wave in a resonator the wave number is real but the frequency is complex. The complex part of the frequency is the damping in time. The resulting pressure wave as the solution of the wave equation is of course the same as by decomposing into plane waves with complex wave number. But with the complex frequency description acoustic attenuation can be included in all "k-space" methods well known in photoacoustic tomography just by introducing a factor describing the exponential decay in time (Roitner & Burgholzer, 2011).

diameter of the detecting part the better is the spatial resolution. A typical single mode fiber for near infrared has a core diameter of 9 microns; single mode fibers for the visible range of detection wavelength have typically about 6 microns. Due to the constant diameter along the whole line detector this type of integrating line detector is dedicated for the imaging of big samples. After first implementations of a Mach-Zehnder and a Fabry-Perot interferometer in glass fibers now polymer fibers are used. Due to the much better impedance matching of polymer fibers to the surrounding water, their sensitivity is higher than for glass fibers, where approximately 2/3 of the incoming signal is reflected before reaching the core (Grün et al., 2010). Furthermore the Young's modulus is much lower in polymer fibers than in glass fibers for which reason the deformation of a polymer fiber is bigger than of a glass fiber applying the same pressure wave. As the strain optic coefficients are in the same order, this results in an enhanced change of refractive index and thus to

(Nuster et al, 2009) did a comparison of the different implementations of an integrating line detector. At this stage of development the free-beam Mach-Zehnder interferometer was the most sensitive integrating line detector. But these measurements showed some new approaches how the fiber-based line detector could be made more sensitive, e.g. by building

The next step after developing a sensitive line detector, no matter of which approach is the most sensitive one, is the creation of an array of many integrating line detectors, e.g. 200 detectors arranged in a curve around the sample. This way one could acquire all data for a

The imaging resolution in photoacoustic imaging is limited by the acoustic bandwidth and therefore by the laser pulse duration as mentioned above, but also by the attenuation of the acoustic wave on its way to the sample surface, and finally by the bandwidth and size of the ultrasonic transducer. The acoustic attenuation can be substantial for high frequencies. This effect is usually ignored in reconstruction algorithms but can have a strong impact on the resolution of small objects or structures within objects. Stokes could already show in 1845 that for liquids with low viscosity, such as water, the acoustical absorption increases by the square of the frequency (Stokes, 1845). If we do not take acoustic attenuation for photoacoustic image reconstruction into account, especially small structures (corresponding to shorter wavelengths and therefore higher frequencies) appear blurred (La Riviere et al., 2006). To what extent this blurring can be compensated by regularization methods (as performed by (La Riviere et al., 2006)) and how much information is lost due the

For thin layers (1D), small cylinders (2D), and small spherical inclusions (3D) the effect of attenuation is simulated and experimental results for several types of tissue are given. For photoacoustic tomography a new description of attenuation seems to be useful: like for a standing wave in a resonator the wave number is real but the frequency is complex. The complex part of the frequency is the damping in time. The resulting pressure wave as the solution of the wave equation is of course the same as by decomposing into plane waves with complex wave number. But with the complex frequency description acoustic attenuation can be included in all "k-space" methods well known in photoacoustic tomography just by introducing a factor describing the exponential decay in time (Roitner & Burgholzer, 2011).

higher signal amplitudes in the polymer fiber (Kiesel et al., 2007).

up a Fabry-Perot interferometer in a single mode polymer fiber.

projection image within one excitation pulse of the laser.

irreversibility of attenuation is investigated in this chapter.

**3. Acoustic attenuation** 

Acoustic attenuation is an irreversible process and therefore the wave equation is not invariant to time reversal. Several important reasons for acoustic attenuation have been reported: viscosity, heat conduction, relaxation processes and chemical reactions. Stokes derived the scalar wave equation

$$
\Delta p - \frac{1}{c^2} \frac{\partial^2 p}{\partial t^2} + \pi \Delta \frac{\partial p}{\partial t} = 0 \tag{2}
$$

under the assumption of adiabatic conditions (thus neglecting the loss due to heat conduction) (Stokes, 1845). This equation can also be found in (Shutilov, 1988) and can generally be derived from a relaxation behavior of pressure and density (Royer & Dieulesaint, 2000), where the density change follows the pressure change after a relaxation time τ . If τ is further expressed in terms of viscosity and specific heat, this equation is also known as the thermoviscous equation. Eq. (2) describes acoustic attenuation which is approximately proportional to the square of the frequency. Other wave equations can describe a more general power law frequency dependence of the attenuation of the form <sup>0</sup> ( ) *<sup>y</sup>* αω α ω = ⋅ (0 < y < 3). (Szabo, 1994) has suggested adding the loss term *L t*() ( ) ∗ *p* **r**,*t* to the wave equation in order to account for such attenuation behavior:

$$
\Delta p(\mathbf{r}, t) - \frac{1}{c\_0^2} \frac{\partial^2}{\partial t^2} p(\mathbf{r}, t) + L(t) \ast p(\mathbf{r}, t) = 0 \tag{3}
$$

Other models for acoustic wave propagation in acoustic media have been proposed also by (Nachman et al., 1990) and (Treeby & Cox, 2010).

**Measurement and simulation of broadband acoustic attenuation:** 

Fig. 3. Experimental set-up to measure broadband high frequency acoustic attenuation in tissue with ultrasound generated by short laser pulses (a) and by piezoelectric transducers (b). In both cases a piezoelectric transducer receives the ultrasonic signals

Compensation of Ultrasound Attenuation in Photoacoustic Imaging 199

energies. With these parameters a complex wave number is obtained involving a

*n n t t t Ht* κ

 <sup>∞</sup> =

χ

ω

over the fat thickness and then applying fat attenuation over the same fat thickness.

( ) ( ) exp( / ) ( ) *N n*

τ

. Stokes' equation (2) is seen to be the special case of a single relaxation mechanism 1 *N* = ,

All three simulation methods produce waveforms in good agreement to the measured waveform. Of course, the approximation will become less accurate with longer propagation

Simulation results may also be validated if we recall that the frequency domain transfer function of the pressure signal propagating through a layer with complex wavenumber

may be obtained from the measured 'water-only' waveform by inverting water attenuation

When taking acoustic attenuation into account, the wave equation, e.g. (2) or (3), is not invariant to time reversal any more. First order terms in time t or higher order odd terms change sign if time is reversed. The equations then do not behave "well" any more, noise is amplified exponentially and regularization methods have to be used for solving these time reversed equations. Using such regularized methods the spatial resolution, that is limited by the frequency-dependent damping, is improved (Burgholzer et al., 2007). To prevent highfrequency noise from growing exponentially, Fourier spectral methods (Trefethen, 2000) in space have been used. They utilize the spatial Fourier transform to calculate the Laplacian and therefore allow for incorporating a damping of higher frequencies when calculating the time reversal. A similar algorithm to compensate acoustic attenuation step by step was

In the second strategy, attenuation is compensated in one step where an approximation to the 'un-attenuated' measured signal is calculated from the attenuated measured signal by solving an appropriate integral equation (La Riviere et al., 2006); then an ordinary reconstruction algorithm for photoacoustic tomography is applied as in the absence of

The Matlab Toolbox by Treeby et al. implements a state-of-the-art algorithm to simulate ultrasound wave propagation ('the direct problem') and image reconstruction ('the inverse problem') for PAI. For the direct and inverse problems arbitrary source distributions and detector geometries, variable medium densities and sound speeds and an arbitrary constant medium absorption are supported in 1D-3D. The algorithm is first-order finite difference in time and pseudospectral in space (working in spatial Fourier space = 'k-space'). The

> *<sup>p</sup> <sup>t</sup>* ρ∗ <sup>∂</sup> =− ∇

**<sup>u</sup>** , *<sup>t</sup>*

ρ

∂

. In the case of absorbing media the equation of state is extended to

ρ

**u** together with an

<sup>∂</sup> <sup>∗</sup> =− ∇⋅

are solved up to time *T* for the sound velocity vector **u** ,

∂

*n*

=+ − where 1

*N n n*

= = − and

κ

κ ∞ χ

denotes the usual compressibility calculated with

*d* . So the simulated 'water-fat-water' waveform

τ

( ) *c* <sup>−</sup> = which holds in liquids as well as in soft biological tissues of density

generalized compressibility 1

and thickness *d* equals exp(iK( ) )

**4. Compensation of acoustic attenuation** 

linearized Euler equations in a fluid <sup>1</sup>

ρ

<sup>0</sup> *p* = *c* ρ

*H t*( ) is the Heaviside step function.

 ρ

distance in the absorbing medium.

proposed by (Treeby et al., 2010c).

adiabatic equation of state <sup>2</sup>

pressure *p* and density

attenuation.

the formula 2 1 χ

ρ

κ <sup>1</sup> = χ.

*K*( ) ω

κ

 κδ

Fig. 3 shows the set up to determine frequency dependent acoustic attenuation in tissue by a single transmission experiment. High frequency ultrasonic pulses are either generated by a pulse laser (pulse duration 6 ns) heating up the surface of a silicon wafer (Fig. 3a) or by a piezoelectric transducer (Fig. 3b). The resulting planar-like ultrasound waves propagate through a fat tissue with varying thickness and through distilled water as a coupling medium to a piezoelectric transducer. The ultrasonic attenuation α is determined by comparing transmission measurements of the investigated samples and distilled water. In Fig. 4 Fourier transformed attenuation results for subcutaneous fat of pig, human blood and olive oil are shown. A power law 0 ( ) *<sup>y</sup>* α α *f f* = ⋅ can be applied to those substances, where f denotes the frequency (Bauer-Marschallinger et al., 2011).

Fig. 4. Attenuation as a function of frequency for three biological substances (double logarithmic scale)

The attenuated planar ultrasound waveform can be calculated by a number of available simulation methods. In a comprehensive study (Roitner et al., 2011) we compared measured waveforms for two fat thicknesses (3.2 mm and 6.2 mm) to simulated waveforms obtained with three current simulation methods. The methods using the Matlab toolbox by (Treeby & Cox, 2010b) or the relation (4) by (LaRiviere et al., 2006) rely on a frequency power law absorption and are described in detail below. A third method by (Nachman et al., 1990) assumes multiple molecular relaxation processes as the cause of attenuation. These processes are characterized by their contributions (1 ) κ *<sup>n</sup> n N* = … to the isothermal compressibility of the tissue and the relaxation times (1 ) τ*<sup>n</sup> n N* = … of their vibrational

Fig. 3 shows the set up to determine frequency dependent acoustic attenuation in tissue by a single transmission experiment. High frequency ultrasonic pulses are either generated by a pulse laser (pulse duration 6 ns) heating up the surface of a silicon wafer (Fig. 3a) or by a piezoelectric transducer (Fig. 3b). The resulting planar-like ultrasound waves propagate through a fat tissue with varying thickness and through distilled water as a coupling

comparing transmission measurements of the investigated samples and distilled water. In Fig. 4 Fourier transformed attenuation results for subcutaneous fat of pig, human blood and

α

*f f* = ⋅ can be applied to those substances, where f

is determined by

medium to a piezoelectric transducer. The ultrasonic attenuation

α

 α

Fig. 4. Attenuation as a function of frequency for three biological substances (double

processes are characterized by their contributions (1 )

compressibility of the tissue and the relaxation times (1 )

The attenuated planar ultrasound waveform can be calculated by a number of available simulation methods. In a comprehensive study (Roitner et al., 2011) we compared measured waveforms for two fat thicknesses (3.2 mm and 6.2 mm) to simulated waveforms obtained with three current simulation methods. The methods using the Matlab toolbox by (Treeby & Cox, 2010b) or the relation (4) by (LaRiviere et al., 2006) rely on a frequency power law absorption and are described in detail below. A third method by (Nachman et al., 1990) assumes multiple molecular relaxation processes as the cause of attenuation. These

κ

τ

*<sup>n</sup> n N* = … to the isothermal

*<sup>n</sup> n N* = … of their vibrational

olive oil are shown. A power law 0 ( ) *<sup>y</sup>*

logarithmic scale)

denotes the frequency (Bauer-Marschallinger et al., 2011).

energies. With these parameters a complex wave number is obtained involving a generalized compressibility 1 ( ) ( ) exp( / ) ( ) *N n n n n t t t Ht* κ κ κδ τ τ <sup>∞</sup> = =+ − where 1 *N n n* κ ∞ χ κ = = − and *H t*( ) is the Heaviside step function. χ denotes the usual compressibility calculated with the formula 2 1 χ ρ ( ) *c* <sup>−</sup> = which holds in liquids as well as in soft biological tissues of density ρ . Stokes' equation (2) is seen to be the special case of a single relaxation mechanism 1 *N* = , κ <sup>1</sup> = χ.

All three simulation methods produce waveforms in good agreement to the measured waveform. Of course, the approximation will become less accurate with longer propagation distance in the absorbing medium.

Simulation results may also be validated if we recall that the frequency domain transfer function of the pressure signal propagating through a layer with complex wavenumber *K*( ) ω and thickness *d* equals exp(iK( ) ) ω *d* . So the simulated 'water-fat-water' waveform may be obtained from the measured 'water-only' waveform by inverting water attenuation over the fat thickness and then applying fat attenuation over the same fat thickness.
