k-Wave Toolbox Previous   Next

Modelling Wave Propagation

Modelling ultrasound propagation

When an acoustic (eg. ultrasonic) wave passes through a medium there are fluctuations in pressure, density, temperature, particle velocity, etc. As most microphones and hydrophones are sensitive to some combination of acoustic pressure and particle velocity, it is these fields that are most interesting and which k-Wave gives as outputs for the user.

Another reason for focussing on pressure and particle velocity is that most acoustic sources can be characterised as either pressure or velocity sources (sometimes categorised equivalently as 'mass' and 'force' sources) and k-Wave allows either to be defined. The key difference between them is that velocity is a vector so a velocity source has an inherent direction associated with it, whereas a pressure source will radiate in all directions (although the shape of a pressure transducer may focus the field more strongly in one direction than another). An example of a velocity source is the surface of a piezoelectric transducer. If the velocity of the surface can be measured then that could be used as an input to k-Wave to model the acoustic field from that transducer.

k-Wave treats the medium as a sound-absorbing fluid in which the absorption has the same frequency-dependence as in real tissue. This is accurate for situations in which the shear modulus is negligible, such as in soft biological tissue. The ability to model, efficiently, ultrasound propagation in which the absorption depends on an arbitrarily specifiable power of the frequency is one of k-Wave's strengths. It can also handle heterogeneities in the background density and sound speed, and model cummulative nonlinear effects. Combined with its computational efficiency, these features make it a useful tool for those working in acoustics, ultrasound, or related fields.

 Back to Top

Initial value problem for modelling photoacoustics

The forward problem in photoacoustics is to calculate the acoustic (ultrasonic) field as a function of time following the absorption of an optical pulse by an elastic fluid with heterogeneous optical absorption. (When RF radiation is used in place of optical, the term thermoacoustic is often preferred.) If, by absorbing the optical pulse, a region of fluid is heated sufficiently quickly to generate a local pressure increase, then an acoustic wave will result. When the heating occurs on a time-scale much shorter than the time necessary for the tissue to relax (for the density to change), the forward problem reduces to an initial value problem (IVP), with the particle velocity zero and an initial acoustic pressure distribution proportional to the absorbed optical energy density. Exact solutions to this IVP are only known in simplified cases, so in general it is necessary to seek a numerical solution. The forward model included in this toolbox uses a computationally efficient k-space approach, a technique closely allied to pseudo-spectral methods, and it can be used to model wave propagation in one, two, or three dimensional space.

 Back to Top

Photoacoustic image reconstruction

The image reconstruction problem in photoacoustic imaging is an acoustic inverse initial value problem: estimate the initial pressure distribution (photoacoustic source region) given a set of acoustic pressure time histories recorded at every point on a measurement surface. A number of algorithms have been proposed to solve this problem. Two are included here, as they are closely related to the k-space forward models described above:

In time-reversal image reconstruction, the k-space forward model (which is symmetrical with respect to time) is used with zero initial conditions, but with the measured acoustic pressure time histories imposed as a time-varying Dirichlet boundary condition at the position of the detectors on the measurement surface. When the measurement surface completely encloses the initial pressure distribution, this model will rebuild the acoustic pressure field inside the measurement surface in time-reversed order. When this is continued back to time t=0, the final pressure field will be the initial pressure distribution - the required photoacoustic image. Examples are given for two or three dimensional mediums.

When the measurement surface is a plane and the medium is acoustically homogeneous and non-absorbing, a much quicker algorithm is available that calculates the initial pressure distribution in a single step. In essence, the algorithm maps the measured data, recorded as a function of position on the plane (2D) and time, to the image (3D); it maps the time-domain information into a spatial dimension. (It relates the temporal frequency to spatial frequency in the depth direction, which are linked by a dispersion relation.) The advantage of this algorithm is that it can be implemented using FFTs, and is thus very fast. Examples are given for two or three dimensional mediums.

 Back to Top

The model equations

The model of linear sound propagation in soft biological tissue used in k-Wave assumes that the medium is isotropic (its material properties do not depend on the direction the wave is travelling), quiescent (there is no net flow) and that shear waves can be neglected. Absorption is included as a phenomenological energy loss term (ie. it does not model specific absorption mechanisms - such as thermal conduction - but instead models the observed frequency-dependent behaviour of absorption in real tissue). The propagation medium is characterised by four parameters, sound speed and density, which can vary spatially, and the absorption coefficient prefactor and exponent.

The three equations from which the linear acoustic equations are derived embody three physical principles: conservation of momentum, conservation of mass, and conservation of energy. The respective equations are:

conservation of momentum

conservation of mass

conservation of energy

where p is the pressure, v is the particle velocity vector, rho density, and s the specific entropy. A is the rate of energy absorption per unit volume and will depend on the acoustic field (it will be frequency dependent in general). The variables p, v, rho and s consist of the ambient values plus the acoustic fluctuation. By considering the acoustic parts (indicated with a prime ' ) as small compared to the ambient values (indicated with a zero subscript), the momentum and mass equations can be linearised to give the following linear equations in the acoustic variables:

linearised conservation of momentum

linearised conservation of mass

The thermodynamic identity (which holds for every fluid particle even when moving)

can then be used to transform the energy equation into a linearised form (here c is the sound speed, beta the volume thermal expansivity, and Cp the specific heat capacity)

linearised conservation of energy

The absorption term A is assumed to depend on the rate of change of density

where, to simplify the notation slightly and for computational convenience, a new quantity has been defined as

Note that in a medium with uniform background density rho'' = rho'. In this notation the energy equation becomes

By integrating this equation with respect to time, and by writing the mass equation in terms of the new variable, the three equations used in k-Wave can be written as:

linearised conservation of momentum

linearised conservation of mass

acoustic equation of state

k-Wave solves these coupled first order acoustic equations using a k-space approach, described below. When medium.BonA is specified, additional nonlinear terms are also included. The equations can be solved either in conjunction with space- and time-varying velocity or pressure source terms, or with an initial condition for the velocity or pressure. The first will be useful for modelling ultrasound from conventional transducers, the last for modelling photoacoustic wave generation.

A unique and useful feature of k-Wave is that the absorption operator has been designed to model the effect of a frequency-dependent absorption coefficient of the form

absorption coefficient

which is often a good model for absorption in biological tissues. Classical absorption can only explain a frequency squared (y=2) dependence, whereas a dependence between y=1 and y=2 is much more typical of soft biological tissue. This absorption prefactor and exponent (alpha0 and y) are the two material properties, along with density and sound speed, that are used to characterise the tissue. To model this power law dependence accurately, efficiently, and in a way that satisfies causality (Kramers-Kronig) an operator based on the fractional Laplacian is used (further details can be found in Treeby & Cox, J. Acoust. Soc. Am., 127 (5), 2741-2748, 2010.

 Back to Top

Pseudo-spectral and k-space methods

The most commonly-used numerical methods for solving partial differential equations are finite-difference (FD) and finite-element methods. Although excellent for many applications, for time-domain modelling of high-frequency wave applications they can become cumbersome and slow. This is due to the requirements for many mesh points per wavelength and small time steps to minimize unwanted numerical dispersion. The pseudospectral method can help reduce the first of these problems, and the k-space method can help to overcome the second.

In a simple FD scheme, the gradient of the field is estimated using linear interpolation between its values at the grid points (i.e., the mesh nodes). A better estimate of the gradient can be obtained by fitting a higher-order polynomial to a greater number of nodes and calculating the derivative of the polynomial. The more points used, the higher the degree of polynomial required, and the more accurate the estimate of the derivative. The pseudospectral method method takes this idea further and fits a Fourier series to all of the data; it is therefore sometimes referred to as a global, rather than local, method. There are two significant advantages to using Fourier series: first, the amplitudes of the Fourier components can be calculated efficiently using the fast Fourier transform (FFT), and second, the basis functions are sinusoidal, so only two nodes per wavelength are required, rather than the six to ten required in other methods.

While the pseudospectral method improves efficiency in the spatial domain, conventional finite-difference schemes are still necessary to calculate the gradients (rates of change) in the time domain. The finite-difference approximation introduces instability into the numerical simulation that can only be controlled by limiting the size of the time-step. The techniques broadly classed as k-space methods attempt to relax this limitation in order to allow larger time-steps to be used without compromising accuracy. By comparing a simple pseudospectral time domain model for acoustically homogeneous media to an exact solution to the corresponding homogeneous wave equation, it is possible to find replacement expressions for either the temporal or the spatial derivative such that the numerical solutions are exact for arbitrarily large time-steps. In effect, this substitution incorporates a priori information about the form of the derivative specific to the governing wave equation. These k-space adjustments also lead to improved numerical stability in the case of acoustically heterogeneous media. For the range of heterogeneity evident in soft biological tissue, this allows the use of much larger time steps for the same degree of accuracy.

For further information see the reference list or post a question to the online k-Wave User Forum.

 Back to Top


© 2009-2012 Bradley Treeby and Ben Cox.