**4. Compensation of acoustic 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 proposed by (Treeby et al., 2010c).

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 attenuation.

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 linearized Euler equations in a fluid <sup>1</sup> *<sup>p</sup> <sup>t</sup>* ρ∗ <sup>∂</sup> =− ∇ ∂ **<sup>u</sup>** , *<sup>t</sup>* ρ ρ <sup>∂</sup> <sup>∗</sup> =− ∇⋅ ∂ **u** together with an

adiabatic equation of state <sup>2</sup> <sup>0</sup> *p* = *c* ρ are solved up to time *T* for the sound velocity vector **u** , pressure *p* and density ρ. In the case of absorbing media the equation of state is extended to

Compensation of Ultrasound Attenuation in Photoacoustic Imaging 201

A Tukey (tapered cosine) window is chosen for this filter. This window is parameterized by

space) and its taper ratio *r* . The ratio controls the filter roll-off: raising *r* suppresses more high-frequency content in the measurement noise. In Fig. 5a we display a complex phantom consisting of nine circles of different sizes and source strengths, surrounded by a circular array of detectors (which could represent line detectors in 3D). This phantom was reconstructed first with taper ratio 0.5 *r* = (Fig. 5b) and then with taper ratio 0.75 *r* = (Fig. 5c). The absorption parameters were an exponent of *y* = 1.5 and 1.5

Note that in Fig. 5b ring-like artifacts are produced during reconstruction, which are much reduced in Fig. 5c because there the regularization filter is more restrictive. With taper ratio *r* = 1 they disappear altogether. If a simpler phantom is used, e.g. just a single circle, taper ratio 0.5 *r* = does not produce any artifacts. To sum up, complex, fine-structured objects exhibit more signal content and hence more measurement noise at higher frequencies. Then the regularization filter must be parameterized more restrictively than for compact objects. The formula used for calculating and compensating attenuation in one step was first

(,) ( , )exp( ( ) ) ( ) *att S ideal S <sup>p</sup> p t ic K t dt c K*

This relation yields the temporal Fourier transform of the attenuated pressure at a given point = *<sup>S</sup>* **r r** (and, via an IFT, the attenuated pressure itself) if the ideal (i.e. un-attenuated) pressure is known over time at that point. This can be exploited for the direct problem where the evaluation involves only two nested numerical integrations. For the inverse problem, the left-hand side of (4) is known and the relation represents an integral equation to be solved for ( ,) *ideal S p* **r** *t* . Since this solution is very sensitive to measurement noise, i.e. noise in the spectrum of ( ,) *att S p* **r** *t* , regularization is necessary, for instance in the form of

If the detector geometry is regular, fast reconstruction algorithms using series expansions in eigenfunctions of the Laplacian on the regular domain are available for the lossless case making the whole reconstruction with compensation much faster than reconstruction with the mentioned toolbox algorithm. In Fig. 5d we display a reconstruction of the nine-circle phantom using this procedure. To apply the SVD, a certain noise level in temporal Fourier space must be assumed. It is desirable to base this assumption on physical reasons - this will be discussed below. Here we assume a white Gaussian noise with a standard deviation of

Fig. 5d shows that the reconstruction with the one-step compensation is almost of the same quality as the reconstruction with the Matlab toolbox algorithm. However, computation

We already observed that the quality of reconstructions depends crucially on a good analysis of the amplitude and spectral distribution of the measurement noise. On the source side, advance knowledge of the type of object to be imaged is helpful (more bulky and compact, or more fine-structured). On the detector side, noise is created by fundamental

⋅ *rad m* , corresponding to 2

0

= ⋅ <sup>⋅</sup> **r r** (4)

 ω π

0 α

⋅ 15*MHz* in

= ⋅ 3 /( ) *dB cm MHz* ,

ω-

π

its cutoff-frequency in *k* -space (here <sup>4</sup> 2 10 /

presented by (LaRiviere et al., 2006) and reads

0,1% of the largest amplitude in the spectrum.

physical processes as well as by technological limitations.

and the 'measurement noise' consisted only of numerical noise.

0

ω

ω

ω

truncated singular value decomposition (SVD, see e.g. La Riviere et al., 2006).

time is reduced by a factor of 10 in comparison to the toolbox algorithm.

+∞

−∞

2 1 2 1 0 ˆ( ,) (,) (,) ( ˆ( , )) *y y <sup>t</sup> p tc t k k t t* ρ ρ τ ηρ <sup>−</sup> − − <sup>∂</sup> =+ + ∂ **k r r** F **k** , where τ is related to attenuation and ηto

dispersion (Treeby & Cox, 2010) and the inverse Fourier transform (IFT) in *d* spatial dimensions is defined by <sup>1</sup> ˆ ˆ ( ) ( ( )) (2 ) ( )exp( ) *d d R f f f id* π− − == ⋅ **r k** <sup>F</sup> **k kr k** . With this extension of

the equation of state frequency-domain ultrasound absorption power laws of the form <sup>0</sup> ( ) *<sup>y</sup>* αω α ω= ⋅ can be supported.

For the inverse problem, the same algorithm as for the direct problem is applied, but with zero initial conditions *p*( ,0) 0 **r** = , **ur 0** ( ,0) = and time-varying Dirichlet boundary conditions on detector surface points *S***r** in the form of the time-reversed sensor data ( ,) ( , ) *S meas S p* **r r** *t p Tt* = − . Absorption is compensated by inverting the sign of τ but leaving η unchanged. Additionally, a regularization filter is applied in *k*-space that suppresses the higher frequencies.

Fig. 5. Reconstruction of nine-circle phantom: (a) initial pressure distribution, (b) reconstruction Matlab toolbox with taper ratio r = 0.5, (c) reconstruction Matlab toolbox with taper ratio r = 0.75, (d) reconstruction one-step compensation plus fast series algorithm

dispersion (Treeby & Cox, 2010) and the inverse Fourier transform (IFT) in *d* spatial

the equation of state frequency-domain ultrasound absorption power laws of the form

For the inverse problem, the same algorithm as for the direct problem is applied, but with zero initial conditions *p*( ,0) 0 **r** = , **ur 0** ( ,0) = and time-varying Dirichlet boundary conditions on detector surface points *S***r** in the form of the time-reversed sensor data

unchanged. Additionally, a regularization filter is applied in *k*-space that suppresses the

pixel Nr.

200

300

400

500

600 **d**

pixel Nr.

reconstruction Matlab toolbox with taper ratio r = 0.5, (c) reconstruction Matlab toolbox with taper ratio r = 0.75, (d) reconstruction one-step compensation plus fast series algorithm

Fig. 5. Reconstruction of nine-circle phantom: (a) initial pressure distribution, (b)

200

300

400

500

<sup>600</sup> **<sup>b</sup>**

( ,) ( , ) *S meas S p* **r r** *t p Tt* = − . Absorption is compensated by inverting the sign of

*d d R f f f id* π

τ

− − == ⋅ **r k** <sup>F</sup> **k kr k** . With this extension of

is related to attenuation and

τ

pixel Nr.

pixel Nr.

200 300 400 500 600

200 300 400 500 600

but leaving

ηto

2 1 2 1

*t* ρ

dimensions is defined by <sup>1</sup> ˆ ˆ ( ) ( ( )) (2 ) ( )exp( )

<sup>−</sup> − − <sup>∂</sup>

∂ **k r r** F **k** , where

> pixel Nr. 200 300 400 500 600

> > pixel Nr.

200 300 400 500 600

 ηρ

ˆ( ,) (,) (,) ( ˆ( , )) *y y <sup>t</sup> p tc t k k t*

 τ

=+ +

= ⋅ can be supported.

0

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

higher frequencies.

200

300

400

500

600 **a**

pixel Nr.

pixel Nr.

200

300

400

500

600 **c**

 α ω

αω

η

ρ

A Tukey (tapered cosine) window is chosen for this filter. This window is parameterized by its cutoff-frequency in *k* -space (here <sup>4</sup> 2 10 / π ⋅ *rad m* , corresponding to 2π ⋅ 15*MHz* in ω space) and its taper ratio *r* . The ratio controls the filter roll-off: raising *r* suppresses more high-frequency content in the measurement noise. In Fig. 5a we display a complex phantom consisting of nine circles of different sizes and source strengths, surrounded by a circular array of detectors (which could represent line detectors in 3D). This phantom was reconstructed first with taper ratio 0.5 *r* = (Fig. 5b) and then with taper ratio 0.75 *r* = (Fig. 5c). The absorption parameters were an exponent of *y* = 1.5 and 1.5 0 α = ⋅ 3 /( ) *dB cm MHz* , and the 'measurement noise' consisted only of numerical noise.

Note that in Fig. 5b ring-like artifacts are produced during reconstruction, which are much reduced in Fig. 5c because there the regularization filter is more restrictive. With taper ratio *r* = 1 they disappear altogether. If a simpler phantom is used, e.g. just a single circle, taper ratio 0.5 *r* = does not produce any artifacts. To sum up, complex, fine-structured objects exhibit more signal content and hence more measurement noise at higher frequencies. Then the regularization filter must be parameterized more restrictively than for compact objects.

The formula used for calculating and compensating attenuation in one step was first presented by (LaRiviere et al., 2006) and reads

$$
\tilde{p}\_{\rm att}(\mathbf{r}\_{\rm s}, o) = \frac{o\nu}{c\_0 \cdot K(o)} \stackrel{\leftarrow}{\underset{-}{\longleftrightarrow}} \tilde{p}\_{\rm ideal}(\mathbf{r}\_{\rm s}, t) \exp(i\varepsilon\_0 K(o) \cdot t) \, dt \tag{4}
$$

This relation yields the temporal Fourier transform of the attenuated pressure at a given point = *<sup>S</sup>* **r r** (and, via an IFT, the attenuated pressure itself) if the ideal (i.e. un-attenuated) pressure is known over time at that point. This can be exploited for the direct problem where the evaluation involves only two nested numerical integrations. For the inverse problem, the left-hand side of (4) is known and the relation represents an integral equation to be solved for ( ,) *ideal S p* **r** *t* . Since this solution is very sensitive to measurement noise, i.e. noise in the spectrum of ( ,) *att S p* **r** *t* , regularization is necessary, for instance in the form of truncated singular value decomposition (SVD, see e.g. La Riviere et al., 2006).

If the detector geometry is regular, fast reconstruction algorithms using series expansions in eigenfunctions of the Laplacian on the regular domain are available for the lossless case making the whole reconstruction with compensation much faster than reconstruction with the mentioned toolbox algorithm. In Fig. 5d we display a reconstruction of the nine-circle phantom using this procedure. To apply the SVD, a certain noise level in temporal Fourier space must be assumed. It is desirable to base this assumption on physical reasons - this will be discussed below. Here we assume a white Gaussian noise with a standard deviation of 0,1% of the largest amplitude in the spectrum.

Fig. 5d shows that the reconstruction with the one-step compensation is almost of the same quality as the reconstruction with the Matlab toolbox algorithm. However, computation time is reduced by a factor of 10 in comparison to the toolbox algorithm.

We already observed that the quality of reconstructions depends crucially on a good analysis of the amplitude and spectral distribution of the measurement noise. On the source side, advance knowledge of the type of object to be imaged is helpful (more bulky and compact, or more fine-structured). On the detector side, noise is created by fundamental physical processes as well as by technological limitations.

Compensation of Ultrasound Attenuation in Photoacoustic Imaging 203

harmonic oscillator – called kicked damped oscillator in this chapter). Reconstruction of the size of this perturbation at time t=0 from the measurement after a time t shows how the information about the size of this kick at t=0 gets lost with increasing time if dissipative processes occur. This change of the initial conditions has a significant advantage: it turns out that the variance stays constant in time, while the mean value is a function of time. This facilitates the calculations of the entropy production and of the information loss due to the

In this section we shall assume that the time varying stochastic processes will have Gauss-Markov character. In doing so we do not wish to assume that all dissipative macroscopic processes considered belong to this specific class of Gauss-Markov processes. It may, however, be surmised that a number of real phenomena may, with a certain approximation, be adequately described by such processes (De Groot & Mazur). The advantage of specifying more precisely the nature of the processes considered is that it enables us to discuss, on the level of the theory of random processes, the behavior of entropy production

Following the theory of random fluctuations given e.g. by (De Groot & Mazur), we take as a starting point equations analogous to the Langevin equation used to describe the velocity of

( ) *<sup>d</sup> <sup>t</sup>*

The components of the vector **α** are the random variables ( 1,2,..., ) *<sup>i</sup> α i n* = , having zero mean value at equilibrium. The matrix **M** of real phenomenological coefficients is independent of time. The vector **ε**( )*t* represents white noise, which is uncorrelated at different times. The

Σ

with the mean value **α**( )*t* and the covariance matrix Σ , which is usually also a function of time but for the initial conditions used later on the covariance matrix will be constant in time. Σ is the determinant of the covariance matrix. If the initial value of **α** is given by a

By an adequate coordinate transformation of **α** the matrix **M** can be diagonalized: the

According to the second law of thermodynamics the entropy of an adiabatically insulated system must increase monotonously until thermodynamic equilibrium is established within

> ( ,) ( ) ( , )ln (, ) *<sup>B</sup> f t St k <sup>f</sup> t d f t* = − → ∞ **<sup>α</sup>**

**α α α**

<sup>1</sup> <sup>1</sup> ( ( )) ( ( ))

**αα αα**

<sup>−</sup> −− Σ −

*<sup>T</sup> t t*

**α** (6)

<sup>0</sup> ( ) *<sup>t</sup> t e*<sup>−</sup> = ⋅ **<sup>M</sup> α α** (7)

(8)

**<sup>α</sup> <sup>M</sup> α ε** (5)

=− ⋅ +

*dt*

distribution density of **α** turns out to be an n-dimensional Gaussian distribution:

<sup>2</sup> <sup>1</sup> ( ,) (2 )

*<sup>n</sup> ft e* π

the system. Then the entropy will be set to zero and the entropy at a time *t* is:

=

given **α**<sup>0</sup> , the mean value at a later time *t* will be:

eigenvalues of **M** are the elements of the diagonal matrix.

stochastic process.

and of information loss.

a Brownian particle:

**Gauss-Markov processes and entropy** 

First, for any material volume *V* at temperature *T* pressure fluctuations of variance *Var p k T V* () / = *<sup>B</sup>* χ (Landau&Lifshitz, 1980) will occur where χ is the adiabatic compressibility and *<sup>B</sup> k* is the Boltzmann constant. So if detectors are e.g. thin foils with small volume, this may play a role since noise amplitudes grow like 1/2 *V* <sup>−</sup> .

Second, the fluctuation-dissipation theorem states that sound absorption processes in media create pressure fluctuations. If the absorption is described by Stokes' equation (2), a straightforward calculation yields a power spectral density of these pressure fluctuations

$$S\_P(\text{co}) = \frac{2k\_B T}{\pi} \cdot \mathcal{X}^{-2} \frac{\pi}{(1 + \cos^2 \pi^2)^2} \dots$$

Third, noise is created by the analog-to-digital conversion process during electronic acquisition of the pressure signal. If this quantization noise is assumed to be white, its power spectral density at sampling frequency *sf* will be constant and equal <sup>2</sup> /12 *<sup>s</sup> q f* (Widrow & Kollar, 2008) where *q* is the quantization interval.

We have given only an outline of the considerations that have to be taken into account for a good estimation of the measurement noise. In a concrete application the relative importance of the mentioned contributions will depend on sizes and shapes of objects and detectors, the physics of the ultrasound propagation medium (compressibility, temperature) and the properties of the measurement devices.

### **5. Stochastic processes for modeling acoustic attenuation**

As for any other dissipative process the energy of the attenuated acoustic wave is not lost but is transferred to heat, which can be described in thermodynamics by an entropy increase. This increase in entropy is equal to a loss of information, as defined by Shannon, and no compensation algorithm can compensate this loss of information. This is a limit given by the second law of thermodynamics. But can this limit be found in the algorithms for compensation of acoustic attenuation given in the previous section? Fluctuations of the measured pressure are "amplified" exponentially during the compensation and therefore we suspect a relation between the entropy production caused by acoustic attenuation and the fluctuations. Indeed such a relation is known in statistical physics as the fluctuationdissipation theorem and is due to (Callen & Welton, 1951), (Callen & Green, 1952), and (Greene & Callen, 1952). It represents in fact a generalization of the famous (Johnson, 1928) (Nyquist, 1928) formula in the theory of electric noise.

In this section we use the theory of non-equilibrium thermodynamics presented by S.R. de Groot and P. Mazur (De Groot & Mazur). More about random variables and stochastic processes can be found e.g. in the book about Statistical Physics from (Honerkamp, 1998). An introduction to stochastic processes on an elementary level has been published by. (Lemons, 2002), also containing "On the Theory of Brownian Motion" by (Langevin, 1908). An introduction to Markov Processes is given by (Gillespie, 1992).

In this section the measured pressure signal is treated as a time-dependent random variable with a mean value and a variance as a function of time. To be able to use the results of some "model" stochastic processes given in literature (Ornstein-Uhlenbeck process or damped harmonic oscillator) for a model of photoacoustic reconstruction we have changed the initial conditions: instead of a defined initial value (with zero variance) we have taken the stochastic process at equilibrium before time zero and at a time zero a certain perturbation has been applied to the process (e. g. a rapid change in momentum for the damped harmonic oscillator – called kicked damped oscillator in this chapter). Reconstruction of the size of this perturbation at time t=0 from the measurement after a time t shows how the information about the size of this kick at t=0 gets lost with increasing time if dissipative processes occur. This change of the initial conditions has a significant advantage: it turns out that the variance stays constant in time, while the mean value is a function of time. This facilitates the calculations of the entropy production and of the information loss due to the stochastic process.

### **Gauss-Markov processes and entropy**

202 Acoustic Waves – From Microdevices to Helioseismology

First, for any material volume *V* at temperature *T* pressure fluctuations of variance

compressibility and *<sup>B</sup> k* is the Boltzmann constant. So if detectors are e.g. thin foils with

Second, the fluctuation-dissipation theorem states that sound absorption processes in media create pressure fluctuations. If the absorption is described by Stokes' equation (2), a straightforward calculation yields a power spectral density of these pressure fluctuations

Third, noise is created by the analog-to-digital conversion process during electronic acquisition of the pressure signal. If this quantization noise is assumed to be white, its power spectral density at sampling frequency *sf* will be constant and equal <sup>2</sup> /12 *<sup>s</sup> q f*

We have given only an outline of the considerations that have to be taken into account for a good estimation of the measurement noise. In a concrete application the relative importance of the mentioned contributions will depend on sizes and shapes of objects and detectors, the physics of the ultrasound propagation medium (compressibility, temperature) and the

As for any other dissipative process the energy of the attenuated acoustic wave is not lost but is transferred to heat, which can be described in thermodynamics by an entropy increase. This increase in entropy is equal to a loss of information, as defined by Shannon, and no compensation algorithm can compensate this loss of information. This is a limit given by the second law of thermodynamics. But can this limit be found in the algorithms for compensation of acoustic attenuation given in the previous section? Fluctuations of the measured pressure are "amplified" exponentially during the compensation and therefore we suspect a relation between the entropy production caused by acoustic attenuation and the fluctuations. Indeed such a relation is known in statistical physics as the fluctuationdissipation theorem and is due to (Callen & Welton, 1951), (Callen & Green, 1952), and (Greene & Callen, 1952). It represents in fact a generalization of the famous (Johnson, 1928)

In this section we use the theory of non-equilibrium thermodynamics presented by S.R. de Groot and P. Mazur (De Groot & Mazur). More about random variables and stochastic processes can be found e.g. in the book about Statistical Physics from (Honerkamp, 1998). An introduction to stochastic processes on an elementary level has been published by. (Lemons, 2002), also containing "On the Theory of Brownian Motion" by (Langevin, 1908).

In this section the measured pressure signal is treated as a time-dependent random variable with a mean value and a variance as a function of time. To be able to use the results of some "model" stochastic processes given in literature (Ornstein-Uhlenbeck process or damped harmonic oscillator) for a model of photoacoustic reconstruction we have changed the initial conditions: instead of a defined initial value (with zero variance) we have taken the stochastic process at equilibrium before time zero and at a time zero a certain perturbation has been applied to the process (e. g. a rapid change in momentum for the damped

χ

is the adiabatic

(Landau&Lifshitz, 1980) will occur where

small volume, this may play a role since noise amplitudes grow like 1/2 *V* <sup>−</sup> .

*Var p k T V* () / = *<sup>B</sup>*

*P*

*k T <sup>S</sup>*

π

<sup>−</sup> ω= ⋅

χ

2

properties of the measurement devices.

<sup>2</sup> ( ) ) *B*

χ

2 22

(Widrow & Kollar, 2008) where *q* is the quantization interval.

**5. Stochastic processes for modeling acoustic attenuation** 

(Nyquist, 1928) formula in the theory of electric noise.

An introduction to Markov Processes is given by (Gillespie, 1992).

 τ

(1 + ω .

τ

In this section we shall assume that the time varying stochastic processes will have Gauss-Markov character. In doing so we do not wish to assume that all dissipative macroscopic processes considered belong to this specific class of Gauss-Markov processes. It may, however, be surmised that a number of real phenomena may, with a certain approximation, be adequately described by such processes (De Groot & Mazur). The advantage of specifying more precisely the nature of the processes considered is that it enables us to discuss, on the level of the theory of random processes, the behavior of entropy production and of information loss.

Following the theory of random fluctuations given e.g. by (De Groot & Mazur), we take as a starting point equations analogous to the Langevin equation used to describe the velocity of a Brownian particle:

$$\frac{d\mathbf{a}}{dt} = -\mathbf{M} \cdot \mathbf{a} + \varepsilon(t) \tag{5}$$

The components of the vector **α** are the random variables ( 1,2,..., ) *<sup>i</sup> α i n* = , having zero mean value at equilibrium. The matrix **M** of real phenomenological coefficients is independent of time. The vector **ε**( )*t* represents white noise, which is uncorrelated at different times. The distribution density of **α** turns out to be an n-dimensional Gaussian distribution:

$$f(\mathbf{a}, t) = \frac{1}{\sqrt{(2\pi)^{n} |\Sigma|}} e^{-\frac{1}{2} (\mathbf{a} - \mathfrak{A}(t))^{T} \Sigma^{-1} (\mathbf{a} - \mathfrak{A}(t))} \tag{6}$$

with the mean value **α**( )*t* and the covariance matrix Σ , which is usually also a function of time but for the initial conditions used later on the covariance matrix will be constant in time. Σ is the determinant of the covariance matrix. If the initial value of **α** is given by a given **α**<sup>0</sup> , the mean value at a later time *t* will be:

$$
\overline{\mathbf{a}}(t) = e^{-\mathbf{M}t} \cdot \overline{\mathbf{a}}\_0 \tag{7}
$$

By an adequate coordinate transformation of **α** the matrix **M** can be diagonalized: the eigenvalues of **M** are the elements of the diagonal matrix.

According to the second law of thermodynamics the entropy of an adiabatically insulated system must increase monotonously until thermodynamic equilibrium is established within the system. Then the entropy will be set to zero and the entropy at a time *t* is:

$$S(t) = -k\_{\beta} \int f(\mathbf{a}, t) \ln \frac{f(\mathbf{a}, t)}{f(\mathbf{a}, t \to \ast)} \, d\mathbf{a} \tag{8}$$

Compensation of Ultrasound Attenuation in Photoacoustic Imaging 205

Eq. (14) has been derived previously by the equipartition theorem: the equilibrium energy associated with fluctuations in each degree of freedom is *kBT/2*. We have used the equity of entropy production and information loss. Eq. (14) states a connection between the strength

fluctuation-dissipation theorem in its simplest form for uncorrelated white noise.

Fig. 6. Points on a sample path of the normalized kicked Ornstein-Uhlenbeck process defined by the Langevin eq. (10). The solid lines represent the mean, and mean ± standard deviation of the scaled velocity coordinate. At the time t=0 a value of v0=10 has been added to the scaled velocity. After some time the information of the amplitude gets more and more

It is instructive to determine the least square estimator (Honerkamp, 1998) for the initial

<sup>0</sup> ( )*<sup>r</sup> v v* − :

0

22 2

σ γ

*e v* γ

−

*t*

( ) <sup>2</sup>

γ σ γ

This gives the Tikhonov regularization with 22 2

*t <sup>e</sup> R t*

parameter. The inverse square root of the regularization parameter 0 *v Var v*( ) can be

For modeling of acoustic waves one needs in addition to the dissipation also an oscillating

term. For pure oscillation without damping we have no loss of information (Fig. 7).

() () *<sup>r</sup> v Rt vt* = ⋅ (15)

<sup>−</sup> <sup>=</sup> <sup>+</sup> (16)

2 () *v Var v v* = as regularization

0 0

, and the strength of the dissipation

γ

. This is the

σ

of the fluctuations, given by <sup>2</sup>

lost due to the fluctuations

velocity 0 *v* . If we write for the estimated initial velocity *<sup>r</sup> v*

we calculate *R t*( ) by minimizing the mean error <sup>2</sup>

interpreted as signal-to-noise-ration (SNR) of *v t*( ) .

**Example: kicked harmonic osciallator** 

with the Boltzmann constant *<sup>B</sup> k* . For a constant covariance matrix the above integration results in:

$$S(t) = -\frac{1}{2}k\_{\underline{\mu}}\overline{\mathbf{a}}(t)^{\mathrm{T}}\Sigma^{-1}\overline{\mathbf{a}}(t)\tag{9}$$

Before modeling the attenuated acoustic wave as a Gauss-Markov process we give two simple examples: the Orstein-Uhlenbeck process with only one component of **α** as a model for the velocity of a Brownian particle and the damped harmonic oscillator with two components of **α** .

### **Example: kicked Ornstein-Uhlenbeck process**

If the random vector **α** in eq. (5) has only one component we get the Langevin equation

$$\frac{dv(t)}{dt} = -\gamma \cdot v(t) + \sigma \cdot \eta(t) \tag{10}$$

which was used to describe Brownian motion of a particle. The random variable *v* is the particle velocity, − ⋅ γ *v* is the viscous drag, and σ is the amplitude of the random fluctuations. The Langevin equation governs an Ornstein-Uhlenbeck process, after (Uhlenbeck & Ornstein, 1930), who formalized the properties of this continuous Markov process. Now we assume that initially we have thermal equilibrium with zero mean velocity. At time zero the particle is kicked which causes an immediate change in velocity of <sup>0</sup> *v* . Following eq. (7) the mean value *v t*( ) shows an exponential decay:

$$
\overline{\boldsymbol{\upsilon}}(t) = \boldsymbol{e}^{-\boldsymbol{\chi}\boldsymbol{t}} \cdot \overline{\boldsymbol{\upsilon}}\_0 \tag{11}
$$

The variance of the velocity *Var v*( )is <sup>2</sup> σ 2γ and is constant in time. In Fig. 6 time and velocity are scaled to be dimensionless and the standard deviation (square root of the variance) of the velocity is normalized.

From eq. (9) the information loss equal to the entropy production till the time t after the kick is:

$$
\Delta S(t) = k\_{\text{g}} \frac{\mathcal{Y}}{\sigma^2} \overline{v}(t)^2 \tag{12}
$$

On the other hand the entropy production known from thermodynamics is the dissipated energy Δ*Q* , which is the kinetic energy of the Brownian particle of mass m, divided by the temperature T:

$$
\Delta S(t) = \frac{\Delta Q}{T} = \frac{m\overline{v}(t)^2}{2T} \tag{13}
$$

The thermodynamic entropy production in eq. (12) has to be equal to the loss of information in eq. (13), and therefore we get for the variance of the velocity:

$$\frac{\sigma^2}{2\gamma} = \frac{kT}{m} \tag{14}$$

with the Boltzmann constant *<sup>B</sup> k* . For a constant covariance matrix the above integration

<sup>1</sup> <sup>1</sup> () () () <sup>2</sup> *<sup>T</sup> St k t t <sup>B</sup>*

Before modeling the attenuated acoustic wave as a Gauss-Markov process we give two simple examples: the Orstein-Uhlenbeck process with only one component of **α** as a model for the velocity of a Brownian particle and the damped harmonic oscillator with two

If the random vector **α** in eq. (5) has only one component we get the Langevin equation

*dt*

*v* is the viscous drag, and

<sup>0</sup> *v* . Following eq. (7) the mean value *v t*( ) shows an exponential decay:

( ) () () *dv t vt t*

which was used to describe Brownian motion of a particle. The random variable *v* is the

fluctuations. The Langevin equation governs an Ornstein-Uhlenbeck process, after (Uhlenbeck & Ornstein, 1930), who formalized the properties of this continuous Markov process. Now we assume that initially we have thermal equilibrium with zero mean velocity. At time zero the particle is kicked which causes an immediate change in velocity of

> <sup>0</sup> ( ) *<sup>t</sup> vt e v* <sup>−</sup>γ

velocity are scaled to be dimensionless and the standard deviation (square root of the

From eq. (9) the information loss equal to the entropy production till the time t after the kick

<sup>2</sup> () () *St k vt <sup>B</sup>* γ

On the other hand the entropy production known from thermodynamics is the dissipated energy Δ*Q* , which is the kinetic energy of the Brownian particle of mass m, divided by the

> <sup>2</sup> ( ) ( ) <sup>2</sup> *Q mv t S t T T*

The thermodynamic entropy production in eq. (12) has to be equal to the loss of information

*kT m*

2 2

γ

σ

in eq. (13), and therefore we get for the variance of the velocity:

σ

2

σ 2γ σ η

σ

=− ⋅ + γ

<sup>−</sup> =− Σ **α α** (9)

(10)

is the amplitude of the random

= ⋅ (11)

Δ = (12)

<sup>Δ</sup> Δ= = (13)

= (14)

and is constant in time. In Fig. 6 time and

results in:

components of **α** .

particle velocity, − ⋅

is:

temperature T:

**Example: kicked Ornstein-Uhlenbeck process** 

γ

The variance of the velocity *Var v*( )is <sup>2</sup>

variance) of the velocity is normalized.

Eq. (14) has been derived previously by the equipartition theorem: the equilibrium energy associated with fluctuations in each degree of freedom is *kBT/2*. We have used the equity of entropy production and information loss. Eq. (14) states a connection between the strength of the fluctuations, given by <sup>2</sup> σ , and the strength of the dissipation γ . This is the fluctuation-dissipation theorem in its simplest form for uncorrelated white noise.

Fig. 6. Points on a sample path of the normalized kicked Ornstein-Uhlenbeck process defined by the Langevin eq. (10). The solid lines represent the mean, and mean ± standard deviation of the scaled velocity coordinate. At the time t=0 a value of v0=10 has been added to the scaled velocity. After some time the information of the amplitude gets more and more lost due to the fluctuations

It is instructive to determine the least square estimator (Honerkamp, 1998) for the initial velocity 0 *v* . If we write for the estimated initial velocity *<sup>r</sup> v*

$$
\upsilon\_r = R(t) \cdot \upsilon(t) \tag{15}
$$

we calculate *R t*( ) by minimizing the mean error <sup>2</sup> <sup>0</sup> ( )*<sup>r</sup> v v* − :

$$R(t) = \frac{e^{-\gamma t}}{e^{-2\gamma t} + \sigma^2 \sqrt{2\gamma} v\_0^{-2}}\tag{16}$$

This gives the Tikhonov regularization with 22 2 0 0 σ γ2 () *v Var v v* = as regularization parameter. The inverse square root of the regularization parameter 0 *v Var v*( ) can be interpreted as signal-to-noise-ration (SNR) of *v t*( ) .

### **Example: kicked harmonic osciallator**

For modeling of acoustic waves one needs in addition to the dissipation also an oscillating term. For pure oscillation without damping we have no loss of information (Fig. 7).

Compensation of Ultrasound Attenuation in Photoacoustic Imaging 207

scaled momentum (p sqrt(2 γ / m)/σ )

<sup>2</sup> <sup>0</sup> sin( ) ( ) [cos( ) ] <sup>2</sup>

γ

*m*

 ω

ω

*p tp t e*

Compared to the Ornstein-Uhlenbeck process (eq. (11)) the mean value has not only an exponential decay in time but also an oscillation with a frequency 22 2

mentioned above (Fig. 7) the oscillating term does not change the entropy and no information is lost. In the average only the exponential decay causes production of entropy and the information about the value of 0 *p* gets more and more lost due to the fluctuations.


scaled time (t γ /m )

/*m* (dashed line)

Fig. 9. Total energy of the damped harmonic oscillator (solid line) shows in the average an

In photoacoustic imaging, the laser pulse at a time *t* = 0 generates an initial pressure distribution *p*<sup>0</sup> ( )**r** (see section 2). For numerical calculations we use a discrete space = *<sup>j</sup>* **r r** *(j=1,..,N)*, where *<sup>j</sup>* **r** are *N* points on a cubic lattice with a spacing of Δ*r* within the sample

γ

ω

The entropy production which is proportional to the total energy is shown in Fig. 9.

Fig. 8. Points on a sample path of the normalized kicked damped harmonic oscillator (eq. (17) and (18)). The solid lines represent the mean, and mean ± standard deviation of the scaled oscillation (left) and momentum (right). At the time t=0 a value of p0=10 has been added to the scaled momentum. After some time the information about the value of p0 gets


*t <sup>m</sup> t*

γ

− = − (21)

scaled time (t γ /m )

ωωγ

= − <sup>0</sup> /(4 ) *m* . As


scaled time (t γ /m )

more and more lost due to the fluctuations

10

exponential decay with the time constant

**Attenuated acoustic wave as a stochastic process** 

20

30

energy

40

50

For the mean value of the momentum *p*( )*t* one gets from eq. (7):

scaled osciallation (x ω sqrt(2 m γ)/ σ)

Fig. 7. Points on a sample path of the kicked harmonic oscillator without damping. The momentum used to kick the oscillator can be reconstructed without loss of information at any time after the kick if ω is known

The stochastically damped harmonic oscillator combines the oscillatory and the diffusive behavior and therefore it is a good starting point to model attenuated acoustic waves. The equations of motion are (using the momentum *p* instead of the velocity *v* ):

$$\frac{d\chi(t)}{dt} = \frac{1}{m}\rho(t) \tag{17}$$

$$\frac{dp(t)}{dt} = -\frac{\mathcal{Y}}{m} \cdot p(t) - m\alpha\_0^2 \ge +\sigma \,\,\eta(t) \tag{18}$$

These equations can be combined using a two dimensional random vector **α** = (,) *x p* and were solved already by (Chandrasekhar, 1943) for definite initial conditions *x*(0) and *p*(0). Again we have changed the initial conditions to an oscillator with zero mean values kicked by an initial momentum 0 *p* . In Fig. 8 the damping is chosen to be <sup>0</sup> γ ω= *m* 3 .

Using the fluctuation-dissipation theorem <sup>2</sup> σ γ2 = *kT* one obtains for the distribution function

$$f(\mathbf{x}, p, t) = \frac{1}{2\pi \frac{kT}{a\_0}} e^{-\frac{1}{2mkT} \left(p - \overline{p}(t)\right)^2 - \frac{1}{2kT} ma\_0^2 \left(\mathbf{x} - \overline{\mathbf{x}}(t)\right)^2} \tag{19}$$

where *x t*( ) and *p*( )*t* are the solutions of the ordinary (non-stochastic) damped harmonic oscillator. Then one gets for the information loss from eq. (9)

$$S(t) = -\frac{1}{T}(\frac{1}{2}m\alpha\rho\_b^2\overline{\mathbf{x}}(t)^2 + \frac{\overline{p}(t)^2}{2m}) = -\frac{1}{T}(E\_{\text{pot}} + E\_{\text{kin}})\,\tag{20}$$

which is equal to the entropy from thermodynamics, where *Epot+Ekin* is the total energy of the harmonic oscillator (sum of the potential and kinetic energy). This fact confirms again that the fluctuation-dissipation theorem for the damped harmonic oscillator <sup>2</sup> 2 *<sup>B</sup>* σ γ = *k T* can be derived from the equity of entropy production and information loss.

Fig. 7. Points on a sample path of the kicked harmonic oscillator without damping. The momentum used to kick the oscillator can be reconstructed without loss of information at

equations of motion are (using the momentum *p* instead of the velocity *v* ):

*dt m*

by an initial momentum 0 *p* . In Fig. 8 the damping is chosen to be <sup>0</sup>

2

π

ω

can be derived from the equity of entropy production and information loss.

<sup>1</sup> ( , ,)

oscillator. Then one gets for the information loss from eq. (9)

Using the fluctuation-dissipation theorem <sup>2</sup>

function

γ

0

2 2 0

ω

The stochastically damped harmonic oscillator combines the oscillatory and the diffusive behavior and therefore it is a good starting point to model attenuated acoustic waves. The

> () 1 ( ) *dx t p t*

( ) ( ) ( ) *dp t <sup>p</sup> tmx t*

These equations can be combined using a two dimensional random vector **α** = (,) *x p* and were solved already by (Chandrasekhar, 1943) for definite initial conditions *x*(0) and *p*(0). Again we have changed the initial conditions to an oscillator with zero mean values kicked

> σ γ

*p pt m x xt mkT kT f xpt e kT*

where *x t*( ) and *p*( )*t* are the solutions of the ordinary (non-stochastic) damped harmonic

1 1 ( ) <sup>1</sup> () ( () ) ( ) 2 2 *pot kin p t St m xt E E T mT* = −

which is equal to the entropy from thermodynamics, where *Epot+Ekin* is the total energy of the harmonic oscillator (sum of the potential and kinetic energy). This fact confirms again that the fluctuation-dissipation theorem for the damped harmonic oscillator <sup>2</sup> 2 *<sup>B</sup>*

=− ⋅ − +

2 0

 ση

2 22 <sup>0</sup> 1 1 ( ( )) ( ( )) 2 2

= , (19)

+ =− + , (20)

ω

− −− −

2

ω

*dt m*= (17)

γ ω= *m* 3 .

2 = *kT* one obtains for the distribution

(18)

σ γ= *k T*

any time after the kick if ω is known

Fig. 8. Points on a sample path of the normalized kicked damped harmonic oscillator (eq. (17) and (18)). The solid lines represent the mean, and mean ± standard deviation of the scaled oscillation (left) and momentum (right). At the time t=0 a value of p0=10 has been added to the scaled momentum. After some time the information about the value of p0 gets more and more lost due to the fluctuations

For the mean value of the momentum *p*( )*t* one gets from eq. (7):

$$\overline{p}(t) = p\_o[\cos(\alpha t) - \frac{\mathcal{Y}}{2m}\frac{\sin(\alpha t)}{\alpha}]e^{-\frac{\mathcal{Y}t}{2m}}\tag{21}$$

Compared to the Ornstein-Uhlenbeck process (eq. (11)) the mean value has not only an exponential decay in time but also an oscillation with a frequency 22 2 ωωγ = − <sup>0</sup> /(4 ) *m* . As mentioned above (Fig. 7) the oscillating term does not change the entropy and no information is lost. In the average only the exponential decay causes production of entropy and the information about the value of 0 *p* gets more and more lost due to the fluctuations. The entropy production which is proportional to the total energy is shown in Fig. 9.

Fig. 9. Total energy of the damped harmonic oscillator (solid line) shows in the average an exponential decay with the time constant γ/*m* (dashed line)

### **Attenuated acoustic wave as a stochastic process**

In photoacoustic imaging, the laser pulse at a time *t* = 0 generates an initial pressure distribution *p*<sup>0</sup> ( )**r** (see section 2). For numerical calculations we use a discrete space = *<sup>j</sup>* **r r** *(j=1,..,N)*, where *<sup>j</sup>* **r** are *N* points on a cubic lattice with a spacing of Δ*r* within the sample

Compensation of Ultrasound Attenuation in Photoacoustic Imaging 209

measured signal, the effect of attenuation and of its compensation can be directly seen. The photoacoustic signal of a 0.2 mm thin absorbing layer in glycerin is calculated by using the scheme of Fig. 10 after a time of 4.5 microseconds. In glycerin the acoustic pressure can be described well by the Stokes' equation with a relaxation time of 244 picoseconds (Shutilov,

For an attenuated acoustic wave instead of the wave equation (1) we have used the Stokes equation (2) or the wave equation (3), giving for the spatial Fourier components ˆ *tp* )( *<sup>k</sup>* not

> τ= − *c c* **ρ ρ** and 2 2

The initial pressure distribution is a one dimensional square pulse corresponding to an

The calculated signal after a time of 4.5 microseconds is shown in Fig. 11 (dashed line). The dashed dotted line shows the signal without attenuation (no relaxation time) and therefore

<sup>1</sup> <sup>ˆ</sup> (0) sin( ) 2 *ikr k k p e ka with k <sup>k</sup>*

Fig. 11. Simulation result as an example for Stoke's equation: it is shown how an initial square pulse would look after a time of 4.5 microseconds in glycerine (τ=244ps) (dashed line). The dashed dotted line shows the pulse without attenuation (τ=0), which correspond

For compensation of the acoustic attenuation we use the time reversed process of the Gauss-Markov process from Fig. 10. Similar to the stochastically damped oscillator the entropy production for the acoustic wave is set equal to the loss of information from the Gauss-

to the reconstruction of the initial pressure pulse, and the solid line is the SVD

reconstruction of the initial pressure pulse (see text)

*k k* λπ

πρ

<sup>−</sup> = = (26)

 τ

= *c* **ρ** (25)

2 22 2 2 22 <sup>1</sup> 4 (1 ) <sup>4</sup> *kkk*

absorbing layer of thickness 2a with a = 0.1 mm. In Fourier space we get:

π

only an oscillating term but also an exponential decay: Putting eq. (24) into the Stokes equation (2) one gets:

ωπ

represents the ideal reconstructed image.

1988).

volume *V*. At a time *t* the pressure distribution ( ) *<sup>j</sup> p t* can be represented by a Fourier series (Barret et al. 1995), including the time *t* = 0 with the initial pressure distribution:

$$p\_{\rangle}(t) = p(\mathbf{r}\_{\rangle}, t) = \sum\_{k=1}^{N} \hat{p}\_{k}(t) \, \varphi\_{k}(\mathbf{r}\_{\rangle}) \quad \text{with} \quad \varphi\_{k}(\mathbf{r}) = e^{2\pi i \mathbf{p}\_{k} \cdot \mathbf{r}} \cdot D(\mathbf{r}) \tag{22}$$

*D*( )**r** is a support function which is one within the sample volume *V* and zero outside. **ρ***<sup>k</sup>* are integer points on an infinite 3D lattices. From eq. (1) we get (if acoustic attenuation can be neglected):

$$
\hat{p}\_k(t) = \cos(\alpha\_k t)\hat{p}\_k(0) \qquad \text{with} \quad \alpha\_k^{-2} = 4\pi^2 c^2 \mathbf{p}\_k^{-2}.\tag{23}
$$

Therefore we have only an oscillating term as shown in Fig. 7, but in higher dimensions, and no information is lost.

For an attenuated acoustic wave instead of the wave equation (1) we have used the Stokes equation (2) or the wave equation (3), giving for the spatial Fourier components ˆ ( ) *<sup>k</sup> p t* not only an oscillating term but also an exponential decay:

$$
\hat{p}\_k(t) = \left[A\_k \sin(a\theta\_k t) + B\_k \cos(a\theta\_k t)\right] e^{-\hat{\lambda}\_k t} \tag{24}
$$

where the phase factors *Ak* and *Bk* can be derived from the initial conditions and ω*<sup>k</sup>* and λ*<sup>k</sup>* are functions of **ρ***<sup>k</sup>* . Like for the damped harmonic oscillator this exponential decay causes the loss of information and can be modeled by a Gauss-Markov process. The information content in Fourier space is the same as in real space (Fig. 10). Therefore we can describe the information loss in the attenuated acoustic wave by the same model as for the damped harmonic oscillator, only in higher dimensions, as for each wave-vector-index *k* an oscillator is needed.

Fig. 10. The initial pressure distribution just after the laser pulse is Fourier transformed (FT). The time evolution of the Fourier series coefficients can be described similar to the mean value of a stochatsically damped harmonic oscillator. The pressure distribution after a time t is then calculated by an inverse Fourier transform (IFT)

### **One dimensional (1D) example: photoacoustic signal of a 0.2 mm thin layer in glycerin:**

The effect of acoustic attenuation for the reconstructed image is similar in 1D, 2D, and 3D (Burgholzer et al., 2010a and 2010b). As in 1D the reconstructed image is just the shifted

volume *V*. At a time *t* the pressure distribution ( ) *<sup>j</sup> p t* can be represented by a Fourier series

() ( ,) () ( ) () () ˆ *<sup>k</sup> <sup>N</sup> <sup>i</sup>*

*D*( )**r** is a support function which is one within the sample volume *V* and zero outside. **ρ***<sup>k</sup>* are integer points on an infinite 3D lattices. From eq. (1) we get (if acoustic attenuation can

2 22 2 ˆ ˆ ( ) cos( ) (0) 4 *k k k kk p t t p with c* = =

Therefore we have only an oscillating term as shown in Fig. 7, but in higher dimensions, and

For an attenuated acoustic wave instead of the wave equation (1) we have used the Stokes equation (2) or the wave equation (3), giving for the spatial Fourier components ˆ ( ) *<sup>k</sup> p t* not

> ˆ ( ) [ sin( ) cos( )] *kt k k kk k pt A t B t e*

*<sup>k</sup>* are functions of **ρ***<sup>k</sup>* . Like for the damped harmonic oscillator this exponential decay causes the loss of information and can be modeled by a Gauss-Markov process. The information content in Fourier space is the same as in real space (Fig. 10). Therefore we can describe the information loss in the attenuated acoustic wave by the same model as for the damped harmonic oscillator, only in higher dimensions, as for each wave-vector-index *k* an

)( Real space <sup>0</sup> *rp*

*kkkkk <sup>k</sup> etBtAtp*

Fig. 10. The initial pressure distribution just after the laser pulse is Fourier transformed (FT). The time evolution of the Fourier series coefficients can be described similar to the mean value of a stochatsically damped harmonic oscillator. The pressure distribution after a time t

**One dimensional (1D) example: photoacoustic signal of a 0.2 mm thin layer in glycerin:**  The effect of acoustic attenuation for the reconstructed image is similar in 1D, 2D, and 3D (Burgholzer et al., 2010a and 2010b). As in 1D the reconstructed image is just the shifted

<sup>−</sup> ˆ = [)( sin( ) + cos( )]

ω

FT IFT

ω

where the phase factors *Ak* and *Bk* can be derived from the initial conditions and

 ϕ

 ωπ

> ω

*p t p t p t with e D*

ϕ

2

λ

<sup>−</sup> = + (24)

<sup>⋅</sup>

= = = ⋅ **<sup>ρ</sup> <sup>r</sup> r rr r** (22)

π

**ρ** . (23)

ω*<sup>k</sup>* and

Fourier space

*t*

λ

ω

*trp* ),(

"k-space"

(Barret et al. 1995), including the time *t* = 0 with the initial pressure distribution:

1

ω

only an oscillating term but also an exponential decay:

=

be neglected):

λ

no information is lost.

oscillator is needed.

)0( *<sup>k</sup> p* 

Time t

is then calculated by an inverse Fourier transform (IFT)

*j j k kj k k*

measured signal, the effect of attenuation and of its compensation can be directly seen. The photoacoustic signal of a 0.2 mm thin absorbing layer in glycerin is calculated by using the scheme of Fig. 10 after a time of 4.5 microseconds. In glycerin the acoustic pressure can be described well by the Stokes' equation with a relaxation time of 244 picoseconds (Shutilov, 1988).

For an attenuated acoustic wave instead of the wave equation (1) we have used the Stokes equation (2) or the wave equation (3), giving for the spatial Fourier components ˆ *tp* )( *<sup>k</sup>* not only an oscillating term but also an exponential decay:

Putting eq. (24) into the Stokes equation (2) one gets:

$$\left\|\boldsymbol{\alpha}\right\|^{2} = 4\pi^{2}c^{2}\boldsymbol{\mathfrak{p}}\_{k}{}^{2}\left(1 - \frac{1}{4}c^{2}\boldsymbol{\mathfrak{p}}\_{k}{}^{2}\boldsymbol{\pi}^{2}\right) \text{ and } \left.\mathcal{J}\_{k} = \pi c^{2}\boldsymbol{\mathfrak{p}}\_{k}{}^{2}\boldsymbol{\mathfrak{r}}\right.\tag{25}$$

The initial pressure distribution is a one dimensional square pulse corresponding to an absorbing layer of thickness 2a with a = 0.1 mm. In Fourier space we get:

$$
\hat{p}\_k(0) = \frac{1}{\pi k} e^{-ikr} \sin(ka) \,\text{with}\, k = 2\pi p\_k \tag{26}
$$

The calculated signal after a time of 4.5 microseconds is shown in Fig. 11 (dashed line). The dashed dotted line shows the signal without attenuation (no relaxation time) and therefore represents the ideal reconstructed image.

Fig. 11. Simulation result as an example for Stoke's equation: it is shown how an initial square pulse would look after a time of 4.5 microseconds in glycerine (τ=244ps) (dashed line). The dashed dotted line shows the pulse without attenuation (τ=0), which correspond to the reconstruction of the initial pressure pulse, and the solid line is the SVD reconstruction of the initial pressure pulse (see text)

For compensation of the acoustic attenuation we use the time reversed process of the Gauss-Markov process from Fig. 10. Similar to the stochastically damped oscillator the entropy production for the acoustic wave is set equal to the loss of information from the Gauss-

Compensation of Ultrasound Attenuation in Photoacoustic Imaging 211

Barret, H. H., Denny, J. L., Wagner, R. F. & Myers, K. J. (1995). Objective assessment of

Bauer-Marschallinger, J., Berer, T., Roitner, H., Grün, H., Reitinger, B. & Burgholzer, P.

Bell, A. G. (1880). On the production and reproduction of sound by light: the photophone,

Berer, T., Hochreiner, A., Zamiri, S. and Burgholzer, P. (2010) Remote photoacoustic imaging on solid material using a two-wave mixing interferometer, *Opt. Lett*. 35, 4151-4153 Burgholzer, P., Hofer, C., Paltauf, G., Haltmeier, M. & Scherzer, O. (2005). Thermoacoustic

Burgholzer, P., Grün, H., Haltmeier, M., Nuster, R. & Paltauf, G. (2007). Compensation of

Burgholzer, P., Roitner, H., Bauer-Marschallinger, J. & Paltauf, G. (2010a). Image

Chandrasekhar, S. (1943). Stochastic Problems in Physics and Astronomy, *Reviews of Modern* 

De Groot, S.R. & Mazur, P. (1984) *Non-Equilibrium Thermodynamics,* Dover Publications, Inc.,

Grün, H., Berer, T., Nuster, R., Paltauf, G. & Burgholzer, P. (2010). Three-dimensional

Haltmeier, M., Scherzer, O., Burgholzer, P. & Paltauf, G. (2004). Thermoacoustic computed

Hansen, P. C. (1987). *Rank-deficient and discrete ill-posed problems: Numerical aspects of linear* 

Kiesel, S., Peters, K., Hassan, T. & Kowalsky, M. (2007). Behavior of intrinsic polymer optical

Sciences (Paris) 146, 530-533. Translation by Anthony Gythiel, published in

tomography with large planar receivers, *InverseProbl.* 20, 1663–1673

Johnson, J. B. (1928). Thermal Agitation of Electricity in Conductors, *Phy. Rev*. 32, 97

Landau, L. D. & Lifshitz, E. M. (1980). *Statistical physics*, *Part 1*, Pergamon Press, Oxford Langevin, P. (1908). Sur la théorie du mouvement brownien, Compets rendus Académie des

fibre sensor for large-strain, *Meas. Sci.Technol.* 18, 3144–3154

American Journal of Physics 65, 1079-1081 (1997)

photoacoustic imaging using fiber-based line detectors, *Journal of Biomedical Optics*,

(2011). Ultrasonic attenuation of biomaterials for compensation in photoacoustic

tomography with integrating area and line detectors, *IEEE Trans. Ultrason.* 

acoustic attenuation for high resolution photoacoustic imaging with line detectors,

Reconstruction in Photoacoustic Tomography using Integrating Detectors accounting for Frequency-Dependent Attenuation, *Proc. of SPIE* Vol. 7564 756423-1 Burgholzer, P., Berer, T., Gruen, H., Roitner, H., Bauer-Marschallinger, J., Nuster, R. &

Paltauf, G. (2010b). Photoacoustic Tomography using Integrating Line Detectors

image quality II, *J. Opt. Soc. Am*. A 12(5) 834-852

imaging, *Proc. SPIE* 7899, 789931

*American Journal of Science* 20, 305–324

*Ferroelectr. Freq. Control* 52, 1577–1583

Invited, *Journal of Physics: Conference Series*, 214(1).

Gillespie, D. T. (1992). *Markov Processes*, Academic Press, New York

*Proc. of SPIE* Vol. 6437 643724-1

Callen, H. B. & Welton, T. A. (1951), *Phys. Rev*. 83, 34 Callen, H. B. & Greene, R. F. (1952). *Phys. Rev*. 86, 702

Greene, R. F. & Callen, H. B. (1952). *Phys. Rev*. 88, 1387

15(2), 021306-1 - 021306-8

*inversion*, SIAM, Philadelphia

Honerkamp, J. (1998). *Statistical Physics*, Springer-Verlag

*Physics* 15, 1-89

New York

**8. References** 

Markov process with standard deviation *sk* and is approximated by an exponential decay in time (Fig. 9):

$$
\Delta S(t) = \frac{1}{2} k\_B \sum\_k \frac{1}{s\_k^{\frac{2}{2}}} \exp(-2\lambda\_k t) \hat{p}\_k(0)^2 \tag{27}
$$

$$
\Delta S(t) = \frac{\Delta Q}{T} = \frac{1}{2c\_0^2 \rho\_0 T} \sum\_k \frac{1}{s\_k^2} \exp(-2\lambda\_k t) \hat{p}\_k(0)^2 \tag{28}
$$

which gives a variance of ˆ ( ) *<sup>k</sup> p t* in the spatial Fourier space ("k-space") independent from k:

$$\left|\mathbf{s}\_k\right|^2 = \left|\mathbf{k}\_B\mathbf{c}\_0\right|^2 \rho\_0 T \tag{29}$$

In real space this gives the thermodynamic fluctuations. Or the other way around, the thermodynamic fluctuations have a quantity such that the information loss of the Gauss-Markov process is equal to the dissipated heat divided by the temperature. Using the thermodynamic fluctuations for a detector size of 1 cm2 and the SVD regularization from section 4 the compensated ˆ ( ) *<sup>k</sup> p t* is calculated. By applying a subsequent inverse Fourier transform (IFT) one gets the compensated pressure profile in real space (solid line in Fig. 11).

### **6. Summary and conclusions**

Acoustic attenuation is modeled as a stochastic process: this helps to understand how thermodynamic entropy production and the decrease of information, which is "transported" in the acoustic wave, are closely connected on a microscopic scale. This theoretical insight enables to answer an important question in photoacoustic imaging: what is the highest possible compensation of attenuation and therefore the best spatial resolution one can achieve?

We could show that for thermal fluctuations the information loss of the reconstructed image is equal to the entropy production due to attenuation of the acoustic wave. Therefore it is sufficient to calculate the entropy production from the macroscopic mean values and it is not necessary to take the fluctuations of the pressure into account.

The size and locations of detectors in photoacoustic imaging should be optimized to get the best resolution and sensitivity. Up to now in such models it is not assumed that the pressure is a random variable, which favors small point like detectors. Taking thermal fluctuations and other noise into account will help to get a more realistic model for detectors and the reconstructed images from measured signals with these detectors.

### **7. Acknowledgments**

This work has been supported by the Christian Doppler Research Association, by the Federal Ministry of Economy, Family and Youth, by the Austrian Science Fund (FWF) project numbers S10503-N20 and TRP102-N20, by the European Regional Development Fund (EFRE) in the framework of the EU-program Regio 13, and the federal state Upper Austria.

### **8. References**

210 Acoustic Waves – From Microdevices to Helioseismology

Markov process with standard deviation *sk* and is approximated by an exponential decay in

2 1 1 ( ) exp( 2 ) (0) <sup>ˆ</sup> <sup>2</sup> *B k <sup>k</sup> k k St k t p s* Δ≈ −

> 2 2 0 0

which gives a variance of ˆ ( ) *<sup>k</sup> p t* in the spatial Fourier space ("k-space") independent from k:

In real space this gives the thermodynamic fluctuations. Or the other way around, the thermodynamic fluctuations have a quantity such that the information loss of the Gauss-Markov process is equal to the dissipated heat divided by the temperature. Using the thermodynamic fluctuations for a detector size of 1 cm2 and the SVD regularization from section 4 the compensated ˆ ( ) *<sup>k</sup> p t* is calculated. By applying a subsequent inverse Fourier transform (IFT) one gets the compensated pressure profile in real space (solid line in Fig.

Acoustic attenuation is modeled as a stochastic process: this helps to understand how thermodynamic entropy production and the decrease of information, which is "transported" in the acoustic wave, are closely connected on a microscopic scale. This theoretical insight enables to answer an important question in photoacoustic imaging: what is the highest possible compensation of attenuation and therefore the best spatial resolution one can

We could show that for thermal fluctuations the information loss of the reconstructed image is equal to the entropy production due to attenuation of the acoustic wave. Therefore it is sufficient to calculate the entropy production from the macroscopic mean values and it is

The size and locations of detectors in photoacoustic imaging should be optimized to get the best resolution and sensitivity. Up to now in such models it is not assumed that the pressure is a random variable, which favors small point like detectors. Taking thermal fluctuations and other noise into account will help to get a more realistic model for detectors and the

This work has been supported by the Christian Doppler Research Association, by the Federal Ministry of Economy, Family and Youth, by the Austrian Science Fund (FWF) project numbers S10503-N20 and TRP102-N20, by the European Regional Development Fund (EFRE) in the framework of the EU-program Regio 13, and the federal state Upper

not necessary to take the fluctuations of the pressure into account.

reconstructed images from measured signals with these detectors.

2 2 *k B* 0 0 *s kc T* = ρ

ρ

1 1 ( ) exp( 2 ) (0) <sup>ˆ</sup> <sup>2</sup> *k k k k <sup>Q</sup> S t t p T c Ts*

2

λ

<sup>Δ</sup> Δ= ≈ <sup>−</sup> (28)

2

(29)

(27)

λ

time (Fig. 9):

11).

achieve?

Austria.

**7. Acknowledgments** 

**6. Summary and conclusions** 


**10** 

*France* 

Georges Nassar

*Université de Lille –Nord de France* 

**Low Frequency Acoustic Devices for** 

**Viscoelastic Complex Media Characterization** 

The evolution in consumer expectations in terms of quality and safety in the agro-industry has led to the need to develop new methods of investigating product quality and the processes involved. Many fields of production still rely too much on the know-how of the operators who, with their experience acquired over time, have become key players in the company. In addition, the manufacturing of quality food products frequently relies on artisanal know-how that is difficult to industrialise and often synonymous of high production losses, therefore prohibitive costs. In contrast, so as to limit costs, the industrial production process is often associated with poorer quality. The objective evaluation of product quality involves the development of methods and sensors adapted to the product or

Indeed, beneath an apparent simplicity, agro-industry products have complex physical properties linked to elasticity, viscosity and plasticity. One of the major difficulties lies in the complexity of the processes which depend on numerous physical parameters. The matter is subjected to numerous mechanical, thermal or chemical treatments thus migrating towards

The originality of the approach adopted consists in the study and set up of an ultrasonic measuring device associated with its electronic environment in order to reply to a specific need due to the complexity of the physico-chemical phenomena involved. A global approach to this problem is very tricky as the physical properties of the media evolve

We thus focused on the development of sensors and methods of characterisation dedicated to different phases of the industrial processes. Two very closely linked aspects were

Work has been carried out to develop acoustic and ultrasonic instrumentation designed to monitor the change in state of the matter (liquid-gel transition and product cohesion), then to monitor the evolution of its elastic properties. The process control applications concern the development of a very low frequency, non-destructive monitoring method to reply to

In this document, we report the scientific approach highlighting the design of the ultrasonic sources which dispenses with classic design through the choice of specific resonance modes for the sensors. Their design aims at promoting low frequency resonance in a relatively small scale composite structure. This sensor technology was adapted according to the

viscoelastic or even plastic properties that are more difficult to quantify.

therefore studied targeting product characterisation and process control.

the specificities of the physical properties of the matter.

**1. Introduction** 

the manufacturing process.

significantly throughout processing.

