Compensation of Ultrasound Attenuation in Photoacoustic Imaging

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.


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

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 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.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/cm 2 in the visible spectral range and 100 mJ/cm 2 in the near infrared light.
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 .The initial pressure 0 p at time 0 t = 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 cylindrical surface with the radius ct ⋅ 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.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 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.

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).
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 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 0 () (Szabo, 1994) has suggested adding the loss term () ( ) , Lt p t * r to the wave equation in order to account for such attenuation behavior: 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 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  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) 2) is seen to be the special case of a single relaxation mechanism 1 N = , 1 κ χ = .
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( ) ) . 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.

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 , where τ is related to attenuation and η to dispersion (Treeby & Cox, 2010)  For the inverse problem, the same algorithm as for the direct problem is applied, but with zero initial conditions (, 0 ) 0 p = r , (, 0 )= ur 0 and time-varying Dirichlet boundary conditions on detector surface points S r in the form of the time-reversed sensor data (, ) (, ) Sm e a s S p tp T t =− rr .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., 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 1.5 y = 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 1 r = 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 0 0 (,) (, ) e x p ( ()) () This relation yields the temporal Fourier transform of the attenuated pressure at a given point S = rr (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 (,) . Since this solution is very sensitive to measurement noise, i.e. noise in the spectrum of (,) , 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.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 s f will be constant and equal 2 /12 s qf (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.

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 thermodynam ics b y a n e n t r o p y 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
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: The components of the vector α are the random variables (1 , 2 , . . . , )i α in = , 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: 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 0 α , the mean value at a later time t will be: 0 () 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: www.intechopen.com Acoustic Waves -From Microdevices to Helioseismology 204 with the Boltzmann constant B k .For a constant covariance matrix the above integration results in: 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 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 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 0 v .Following eq. ( 7) the mean value () vt shows an exponential decay: The variance of the velocity () Var v is 2 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: 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: 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: Eq. ( 14) has been derived previously by the equipartition theorem: the equilibrium energy associated with fluctuations in each degree of freedom is k B T/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 2 σ , and the strength of the dissipation γ .This is the fluctuation-dissipation theorem its simplest form for uncorrelated white noise.() vV a r v can be interpreted as signal-to-noise-ration (SNR) of () vt .

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).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 ): These equations can be combined using a two dimensional random vector (,) xp = α and were solved already by (Chandrasekhar, 1943) for definite initial conditions (0) x and (0) p .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 0 3 m γω = .Using the fluctuation-dissipation theorem 2 2 kT σγ = one obtains for the distribution function where () xt and () p t are the solutions of the ordinary (non-stochastic) damped harmonic oscillator.Then one gets for the information loss from eq. ( 9) which is equal to the entropy from thermodynamics, where E pot +E kin 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 2 2 B kT σγ = can be derived from the equity of entropy production and information loss.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 p 0 =10 has been added to the scaled momentum.After some time the information about the value of p 0 gets more and more lost due to the fluctuations For the mean value of the momentum () p t one gets from eq. ( 7): 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 0 /(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.

Attenuated acoustic wave as a stochastic process
In photoacoustic imaging, the laser pulse at a time t = 0 generates an initial pressure distribution () 0 p r (see section 2).For numerical calculations we use a discrete space j = rr (j=1,..,N), where j r are N points on a cubic lattice with a spacing of r ∆ within the sample volume V.At a time t the pressure distribution ( ) j p t can be represented by a Fourier series (Barret et al. 1995), including the time t = 0 with the initial pressure distribution: () D r is a support function which is one within the sample volume V and zero outside.k ρ are integer points on an infinite 3D lattices.From eq. ( 1) we get (if acoustic attenuation can be neglected): 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 ˆ() k p t not only an oscillating term but also an exponential decay: where the phase factors k A and k B can be derived from the initial conditions and k ω and k λ are functions of k ρ .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.
Fourier space "k-space" Real space

Time t
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 and2010b).As in 1D the reconstructed image is just the shifted 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 ) ( ˆt p k not only an oscillating term but also an exponential decay: Putting eq. ( 24) into the Stokes equation (2) one gets:  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 cm 2 and the SVD regularization from section 4 the compensated ˆ() k 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).

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.

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.

Fig. 1 .
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

Fig. 2 .
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

Fig. 4
Fourier transformed attenuation results for subcutaneous fat of pig, human blood and olive oil are shown.A power law to those substances, where f denotes the frequency (Bauer-Marschallinger et al., 2011).

Fig. 4 .
Fig. 4. Attenuation as a function of frequency for three biological substances (double logarithmic scale) 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 (liquids as well as in soft biological tissues of density ρ .Stokes' equation (

=
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 and the inverse Fourier transform (IFT) in d spatial dimensions

Fig. 5 .
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 any material volume V at temperature T pressu re fluctuations of variance (1980) will occur where χ is the adiabatic compressibility and B 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 − .Second, the fluctuation-dissipation theorem states that sound absorption processes in media create pressure fluctuations.If the absorptio n i s d e s c r i b e d b y S t o k e s ' e q u a t i o n ( 2 ) , a straightforward calculation yields a power spectral density of these pressure fluctuations

Fig. 6 .(
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 velocity coordinate.At the time t=0 a value of v 0 =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 r v () ()

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

Fig. 8 .
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 p 0 =10 has been added to the scaled momentum.After some time the information about the value of p 0 gets more and more lost due to the fluctuations

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 γ 2 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: 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 .
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-