2.1.

Acoustic Attenuation in Tissue

The acoustic attenuation within soft biological tissue over the frequency range of interest in biomedical photoacoustics has been experimentally shown to follow a frequency power law of the form¹¹

2.2.

Attenuation Compensation in Photoacoustic Tomography

A simple approach for attenuation compensation is to apply a frequency domain filter of the form $e^{{\alpha }_0{|\omega |}^{y}d}$ to each of the recorded time series before image reconstruction, where the distance $d$ is fixed and chosen in advance (removing the dependency of the filter on both $t$ and $\omega$). This was first discussed by Tan et al.,²⁵ who demonstrated a general improvement in image resolution after attenuation compensation using experimental measurements of a turbid gel phantom with three small hairs located in the central region. A similar approach has recently been used for attenuation compensation when photoacoustic waves propagate through an absorbing layer of known thickness.^15,26,27 However, since $d$ is fixed, this procedure is not suitable when the image features are positioned in an absorbing medium at varying distances from the receiver (as is most often the case in practice). In this case, signals arising from features shallower than the chosen $d$ will be overcompensated, and signals from deeper features will be undercompensated.

A more general approach for attenuation compensation was proposed by La Rivière et al., in which the time-varying pressure signals $p_{\mathrm{att}}$ measured at the tissue surface are related to the equivalent unattenuated (or corrected) signals $p$ via an integral relation of the form $p_{\mathrm{att}}=Lp$.^28,29 The original form of the operator $L$ was derived by relating $p$ and $p_{\mathrm{att}}$ in the temporal frequency domain using a lossy wave equation proposed by Szabo that accounts for power law absorption.^23,30 (This derivation is repeated in general form in the Appendix.) Attenuation compensation can then be achieved by forming a matrix equation and computing the inverse of $L$ to calculate an estimate of $p$ from $p_{\mathrm{att}}$. However, in general, the integral equation $p_{\mathrm{att}}=Lp$ is not well conditioned, which makes this inversion ill-posed.^31,32 In the original, a regularized pseudoinverse of $L$ was computed via singular value decomposition (SVD) by thresholding the eigenvalues.²⁸ To improve the stability of the inversion, Ammari et al. later suggested an asymptotic form of $L$ for small ${\alpha }_0$ that can be more easily inverted using a recursive scheme.³³ Kowar and Scherzer discussed a similar formulation to La Rivière for a range of other lossy wave equations and provided an analysis of the ill-conditioning of the integral equations via SVD.³²

In a slightly different approach, Modgil et al. derived an integral equation relating the measured pressured signals $p_{\mathrm{att}}$ directly to the reconstructed initial pressure distribution $p_0$ (including attenuation compensation) assuming a planar detection geometry.^34–36 The corresponding matrix equation was again inverted using truncated-SVD. In this case, the smaller eigenvalues were shown to correspond to wavelet-like eigenvectors with their energy concentrated deeper into the tissue. This is consistent with the expectation that high frequencies generated by deeper features in the tissue are more difficult to recover. Kowar provided similar integral relations relating $p_{\mathrm{att}}$ to $p_0$ for point, line, and planar detectors in simple geometries based on a variant of Szabo’s equation satisfying more stringent causality requirements.³⁷ These approaches have not yet been applied to experimental photoacoustic data.

An alternate scheme for attenuation compensation in photoacoustic tomography was proposed by Burgholzer et al. based on time reversal image reconstruction with a lossy wave equation.^38,39 In short, time reversal reconstruction works by running a numerical model of the acoustic forward problem backward in time. At each time step, the measured time-varying pressure signals are enforced (in time-reversed order) as a Dirichlet boundary condition at their recorded positions.^40,41 The initial pressure distribution is then recovered when $t$ reaches 0. To compensate for acoustic attenuation, Burgholzer et al. used a numerical implementation of Stokes’ lossy wave equation with the absorption term reversed in sign.³⁸ (Stokes’ equation accounts for power law absorption of the form $\alpha ={\alpha }_0{\omega }^2$.) Reversing the attenuation term causes the acoustic waves in the numerical model to grow, rather than decay, as they propagate back through the tissue. To keep noisy high-frequency signals from growing too quickly, the reconstruction was regularized by filtering the signals above a chosen cutoff frequency.

This approach was later extended by Treeby et al. to account for general power law absorption behavior using a lossy wave equation in which the absorption and dispersion are encapsulated by different terms.^6,19 The separation of the two allows the sign of the absorption term to be reversed to compensate for acoustic absorption, while leaving the sound speed dispersion term unchanged. The latter is important in time reversal attenuation compensation for recovering the correct shape of the propagating waves in dispersive media (i.e., when $y\ne 2$).⁶ This approach has recently been used to compensate for acoustic attenuation in a number of pre-clinical imaging studies.^42–44 The compensation improves the resolution of the reconstructed images and allows deeper features to be visualized.⁴² Huang et al. has also illustrated the applicability of this scheme to media with a heterogeneous distribution of acoustic absorption values.⁴⁵ This was demonstrated through a series of phantom measurements with pencil leads embedded in an agar phantom enclosed by an acrylic shell. The restriction of this approach to attenuation compensation is that time reversal must be used to reconstruct the image. For many canonical detector geometries, this may be less efficient than using other fast reconstruction algorithms.

2.3.

Attenuation Compensation in Seismology

A remarkably similar problem to attenuation compensation in photoacoustic tomography exists in seismology, though it appears that this connection has not been previously explored. In seismology, the propagating seismic waves (generated by an explosive charge, for example) are attenuated due to scattering and absorption in the subsurface of the earth. This is accompanied by a velocity dispersion in which the different frequency components of the waves travel at different speeds. The strength and characteristics of the attenuation are normally quantified by a dimensionless parameter called the earth quality factor $Q$. This is related to the attenuation and dispersion in the medium via the equation⁴⁶

Techniques for attenuation compensation in seismology (often referred to as inverse $Q$ filtering) have a long history, and there is a considerable volume of associated literature (the monograph by Wang provides an excellent starting point for examining the seismic literature on attenuation compensation⁴⁶). Most of these techniques fit into the framework of time-variant filtering. This allows the attenuated seismic traces (analogous to the recorded pressure signals in photoacoustics) to be corrected using a filter dependent on both time and frequency. The different compensation methods vary both in how the time-variant filter is constructed and in how it is applied. One of the simplest approaches is called time-variant spectral whitening.⁴⁷ In this technique, the seismic trace is filtered using a series of narrow bandpass filters as shown in Fig. 1(a). Attenuation compensation is then performed separately for each frequency band, assuming the attenuation is frequency independent (where $\omega$ has a fixed value for each band). The final corrected signal is then obtained by recombining the filtered signals. If the features in the seismic trace can be assumed to be of a similar magnitude, this approach can also be used to estimate the amount of attenuation within each frequency band by fitting an envelope to the traces.⁴⁷

Fig. 1

Schematic of attenuation compensation techniques used in seismology: (a) time-variant spectral whitening (after Yilmaz⁴⁷), (b) Gabor filtering, and (c) wave-field downward continuation.

Time-variant spectral whitening can be considered as a coarse method to decompose the seismic trace into a time-frequency distribution, which is then modified using a time-variant filter. A large number of extensions to this technique are therefore possible depending on how the time-frequency distribution and its inverse are constructed. For example, the Gabor transform (a particular type of short-time Fourier transform) has also been used for attenuation compensation in seismology.^48,49 The Gabor transform works by performing a series of localized Fourier transforms on the input signal using a moving time window. The resulting two-dimensional time-frequency representation of the signal can then be modified using a time-variant filter, as shown in Fig. 1(b). The filtered signal (in this case, the corrected seismic trace) is then recovered by performing an inverse Gabor transform.

An alternate way of applying the time-variant filter without explicitly decomposing the input signal into a time-frequency distribution is via wave-field downward continuation.⁴⁶ In this approach, the frequency spectrum of the seismic trace is multiplied by a frequency domain filter that corrects for acoustic attenuation assuming a particular discrete value of time or travel distance. In this case, the time is chosen to be $t=n\mathrm{\Delta }t$, where $n=0,1,2,\dots ,N-1$, and $\mathrm{\Delta }t$ and $N$ are chosen to match the acquisition settings. The resulting Fourier spectrum is then transformed back to the time-domain, and the $n$’th value of the signal is saved as the $n$’th value of the corrected time series. This procedure is repeated for each value of $n$, as shown in Fig. 1(c).

3.1.

Linear Time-Variant Filtering

While many of the approaches for attenuation compensation in seismology may warrant further investigation into their applicability to photoacoustics, here we focus on a general method for time-variant filtering using a particular form of nonstationary convolution. This is formulated as an extension of stationary convolution where the filter is also allowed to vary as a function of output time (Margrave refers to this as nonstationary combination⁵⁰). This can be expressed mathematically as

Fig. 2

Construction of a nonstationary convolution matrix used for time-variant filtering. (a) A frequency domain filter for each value of discrete time is constructed in the matrix rows. (b) The 1D inverse Fourier transform is performed along each matrix row to give the filter impulse response. (c) The rows are circular shifted (translated) based on their index, starting from 0. (d) The lower triangular part of the lower left quadrant and the upper triangular part of the upper right quadrant (shown in gray) are forced to be zero to make the filter acyclic.

For attenuation compensation in photoacoustics, the time-variant filter $F$ is constructed to counteract the effects of acoustic attenuation following a frequency power law as described by Eq. (5). This yields a filter of the form

The function $W(t,\omega )$ in Eq. (9) is a time-variant window used to regularize the filter. This is required to keep any high-frequency noise in the input signal from being amplified. Here the time-variant window is chosen to be a frequency domain Tukey or tapered cosine window for each discrete value of $t$. In other words, it consists of a series of scaled one-dimensional Tukey windows along each row of the convolution matrix. This choice allows the correct frequency dependence of the filter to be maintained within the flat component of the window (equivalent to a filter passband) while smoothly rolling off to zero at the cutoff frequency. An example of a time-variant filter constructed using Eq. (9) is shown in Fig. 3. In this case, the acoustic parameters are set to those of breast tissue, where $c_0=1510{\mathrm{ms}}^{-1}$, ${\alpha }_0=5.5\times {10}^{-10}\mathrm{Np}{(\mathrm{rad}/\mathrm{s})}^{-y}{\mathrm{m}}^{-1}$ (equivalent to $0.75\mathrm{dB}{\mathrm{MHz}}^{-y}{\mathrm{cm}}^{-1}$), and $y=1.5$. The filter is regularized using a series of scaled Tukey windows with a fixed taper ratio of 0.25 and a cutoff frequency varying from 10 MHz at $t=0$ to 4 MHz at $t=10\mu \mathrm{s}$.

Fig. 3

Example of a time-variant filter $F(t,\omega )$ designed using Eq. (9) to compensate for acoustic attenuation in biological tissue where $c_0=1510{\mathrm{ms}}^{-1}$, ${\alpha }_0=5.5\times {10}^{-10}\mathrm{Np}{(\mathrm{rad}/\mathrm{s})}^{-y}{\mathrm{m}}^{-1}$, and $y=1.5$ (here $\omega =2\pi f$). The filter is regularized using a scaled Tukey window with a fixed taper ratio of 0.25 and a cutoff frequency varying from 10 MHz at $t=0$ to 4 MHz at $t=10\mu \mathrm{s}$.

The computational complexity of applying the time-variant filter to a single time-domain signal is $O(N^2)$, where $N$ is the number of time samples. The complexity of sequentially applying attenuation compensation to multiple signals is thus $O({MN}^2)$, where $M$ is the number of signals or detectors. This means the approach is computationally efficient, as the number of time samples in each recorded signal is typically small, even for high-resolution images (e.g., $N=500$ in Laufer et al.⁴³). This restriction is because the required sampling frequency is ultimately limited by the bandwidth of the detectors, and the required acquisition time (or imaging depth) is ultimately limited by the sensitivity of the detectors and the penetration depth of the excitation laser light. In addition, the computational complexity of using time-variant filtering only grows linearly with the number of measurements or detectors, and this can be easily parallelized.

A simple example of using time-variant filtering to correct for acoustic attenuation and dispersion in a recorded time-domain signal is shown in Fig. 4. The forward signal was simulated using the second-order $k$-space model in the k-Wave MATLAB toolbox.^51,52 The initial photoacoustic pressure distribution was defined as two small discs and recorded by a point detector, as shown in Fig. 4(a). The recorded time signals for simulations in lossless (solid line) and lossy (dashed line) media without noise are shown in Fig. 4(b). The lossless signal is repeated in each panel for reference. Figure 4(c) shows the recorded signal after correction using a conventional time-invariant frequency domain filter of the form $e^{{\alpha }_0{|\omega |}^{y}d}$, equivalent to that used by Tan et al.²⁵ In this case, the propagation distance was fixed to $d=c_{0}t$ where $t=4\mu \mathrm{s}$. Since the filter characteristics do not depend on time, components of the signal arising from features shallower than the chosen distance have been overcorrected, and additional oscillations have been introduced. Figure 4(d) shows the recorded signal after correction using a time-variant filter based on Eq. (9) using only the absorption term (i.e., neglecting sound speed dispersion). In this case, the magnitude of the signal has been corrected, but components of the signal appear shifted due to the effects of sound speed dispersion. Figure 4(e) shows the time signal after correction using a time-variant filter based on Eq. (9) using both the absorption and dispersion terms. In this case, the corrected time signal closely matches the reference signal as desired. This illustrates the effectiveness of using time-variant filtering for attenuation compensation.

Fig. 4

Numerical example of using time-variant filtering to correct for acoustic attenuation and dispersion in a recorded time-domain signal. (a) Schematic of the simulation showing the initial pressure distribution and sensor. (b) The recorded time series in lossless (solid line) and lossy (dashed line) media. The lossless signal is repeated in each panel for reference. (c) Signal after correction using a time-invariant filter assuming a fixed distance (dashed line). (d) Signal after correction using a time-variant filter not accounting for sound speed dispersion (dashed line). (e) Signal after correction using a time-variant filter (dashed line).

3.2.

Adaptive Regularization Using Time-Frequency Analysis

In practice, the recorded photoacoustic time signals will always contain some level of noise. This means that any approach to attenuation compensation must include regularization to prevent noise at frequencies where the signal-to-noise ratio is low from being excessively amplified. In the case of time-variant filtering, the compensation is regularized by windowing the time-variant filter, as shown in Eq. (9). The choice of window shape and cutoff frequency is thus important for maintaining a balance between the fidelity of the compensation and the amount of noise introduced. In previous work on attenuation compensation using time reversal image reconstruction, regularization was achieved by filtering the absorption and dispersion terms in the governing equations using a Tukey window with a fixed cutoff frequency.⁶ The cutoff frequency was chosen by examining where the average power spectrum of the recorded time signals reached the noise floor. A similar approach is possible using time-variant filtering by using a time-invariant window $W(\omega )$ with a fixed cutoff frequency. However, the use of time-variant filtering also introduces the possibility of selecting a time-variant window to regularize the compensation based on the local frequency content of the signal. This is significant, because high frequencies generated by deeper features in the tissue will undergo more attenuation before reaching the detector compared to shallow features (the waves have travelled further). This means that high-frequency components from shallow features (which appear early in the time signal) may be above the noise floor, while the same frequency components from deeper features (which appear later in the time signal) may not. The use of a time-variant window to regularize the attenuation compensation means the filter can adapt to the local frequency content of the recorded signal.

One approach to select the time-varying cutoff frequency for the window is to decompose each recorded time signal into a two-dimensional time-frequency distribution. There are many different ways to decompose a time-domain signal into such a distribution, each with different characteristics.⁵³ In this work, the time-frequency distribution is used only to select the cutoff frequency for the window, not to apply the time-variant filter like the Gabor filtering method discussed in Sec. 2.3, so the exact choice is not critical. A representative example is given in Fig. 5. The Shepp-Logan phantom⁵⁴ was used to define the initial photoacoustic pressure distribution analogous to Fig. 6 in Treeby et al.⁶ The forward simulations were again computed in MATLAB using the k-Wave Toolbox. The computational domain was $1024\times 1024$ grid points ($52\times 52\mathrm{mm}$, including a perfectly matched layer) supporting frequencies up to 14 MHz. The acoustic parameters were set to $c_0=1510{\mathrm{ms}}^{-1}$, ${\alpha }_0=2.2\times {10}^{-9}\mathrm{Np}{(\mathrm{rad}/\mathrm{s})}^{-y}{\mathrm{m}}^{-1}$, and $y=1.5$. The signals were detected using a circle of 200 evenly distributed point detectors defined in Cartesian space with a radius of 25 mm. After the forward simulation, random Gaussian noise was added to the recorded signals to give a signal-to-noise ratio of 40 dB. The time signals recorded at one of the detector positions in both lossless (solid line) and lossy (dashed line) media are shown in Fig. 5(a). Three different time-frequency distributions of the lossy signal—the spectrogram, Wigner-Ville distribution, and Rihaczek distribution, respectively—are shown in Fig. 5(b)–5(d). Each distribution shows similar information on the local frequency content of the signal, but with different types and levels of artifacts (a detailed discussion of the merits of different time-frequency distributions is given by Boashash⁵³).

Fig. 5

Attenuation compensation using time-variant filtering. (a) Example of a time series recorded in a lossy medium (dashed line) against the equivalent signal recorded in a lossless medium (solid line) for an initial photoacoustic pressure distribution defined by the Shepp-Logan phantom. (b) Spectrogram time-frequency distribution of the recorded time signal. (c) Wigner-Ville time-frequency distribution of the recorded signal. (d) Rihaczek time-frequency distribution of the recorded signal. The extracted time-variant cutoff frequency is shown with a dashed white line. (e) Corrected time signal (dashed line) after time-variant filtering.

Fig. 6

Effect of using time-variant filtering to correct for acoustic attenuation on reconstructed images of the Shepp-Logan phantom. (a) Reconstruction with no compensation for acoustic attenuation. (b) Profile through $x=0$. The equivalent profile through the initial pressure distribution is shown with a light line for reference. (c) Reconstruction after time-variant filtering using a cutoff frequency chosen separately for each recorded signal based on the individual time-frequency distributions. (d) Profile through $x=0$. (e) Reconstruction after time-variant filtering using a cutoff frequency derived from the average time-frequency distribution of all the recorded signals. (f) Profile through $x=0$. (g) Reconstruction after time-variant filtering using a fixed cutoff frequency of 3 MHz. (h) Profile through $x=0$.

For the simulations presented here, the Rihaczek distribution was used to derive the time-variant cutoff frequency (other distributions could similarly be used if desired). The two-dimensional Rihaczek distribution $R_z(t,f)$ of a time-varying signal $s(t)$ is calculated by $R_z(t,f)=s(t)S^{*}(f)e^{-i2\pi ft}$, where $S^{*}(f)$ is the complex conjugate of the Fourier transform of $s(t)$.⁵³ Given this distribution, the cutoff frequency for each discrete time-point was selected as the frequency at which a particular percentage of the total energy was reached, in this case 98% (this was chosen heuristically). The variation in the cutoff frequency over time was then smoothed using a piecewise constant cubic spline. The resulting time-varying cutoff frequency for the signal shown in Fig. 5(a) is given as the white dashed line in Fig. 5(d). The corrected time signal after applying the filter is shown in Fig. 5(e).

The reconstructed photoacoustic images of the Shepp-Logan phantom are shown in Fig. 6. The recorded signals were first interpolated onto a continuous circular sensor surface using linear interpolation⁵² and then reconstructed using time reversal (any applicable reconstruction algorithm could similarly be used). The reconstructed image without attenuation compensation is shown in Fig. 6(a). A one-dimensional profile through $x=0$ is shown in Fig. 6(b), with the corresponding profile through the initial photoacoustic pressure distribution shown with a light line. The acoustic attenuation has blurred the edges of the reconstructed image, and the overall magnitude of the image features has been reduced. The image reconstructed after attenuation compensation using time-variant filtering is shown in Fig. 6(c) and 6(d). The filter matrix was constructed using Eq. (9), with the time-varying cutoff frequency for each recorded signal chosen separately based on the Rihaczek time-frequency distribution of each signal. The sharpness of the reconstructed image compared to the uncompensated case is considerably improved, and the reconstructed pressure magnitudes are much closer to quantitatively accurate values. Due to the band-limited frequency response of the recorded signals (some high frequencies are below the noise floor when the acoustic waves reach the detector), the compensation cannot completely recover the frequency content of the initial pressure distribution. This results in a small amount of overshoot (Gibbs phenomenon) near the sharp edges in the image.

The corresponding reconstruction using the average time-frequency distribution of all the recorded signals to select the cutoff frequency for the window is shown in Fig. 6(e) and 6(f). Using the average time-frequency distribution to construct the filter and then applying the same filter to all of the time signals is considerably faster than constructing and applying different filters for each signal. The difference is due to the computational cost of creating each nonstationary convolution matrix as described in Sec. 3.1. Again, the sharpness and magnitude of the reconstructed image are considerably improved relative to the uncompensated case. Using the individual filters gives a marginally sharper image as seen in Fig. 6(c) and 6(d). However, this comes at the expense of an increase in compensation time as mentioned above. For comparison, the reconstruction after attenuation compensation using time-variant filtering with a fixed cutoff frequency of 3 MHz is shown in Fig. 6(g) and 6(h). Using a fixed cutoff frequency is similar to using time reversal reconstruction with attenuation compensation, which is also regularized by a single cutoff frequency (see Fig. 6(e) of Treeby et al.⁶ for comparison). In this case, the edges of the reconstructed image appear sharp, but significantly more noise has been introduced. This is because the same cutoff frequency is used regardless of the local frequency content of the recorded time-domain signals.