Finite element approximation for quantitative photoacoustic tomography in a diffusive regime

TL;DR

This work tackles the reconstruction of the diffusion coefficient $D$ and absorption coefficient $\sigma$ in quantitative photoacoustic tomography (QPAT) from internal data $H(x)=\sigma(x)u(x)$ by decoupling into an inverse diffusivity problem (IDP) and a direct diffusion problem. It uses randomly chosen boundary illuminations to guarantee a local nonzero gradient condition with high probability, enabling a conditional Hölder stability and a weighted-energy error analysis for a finite-element regularized least-squares scheme. An $L^2(\Omega')$ error bound for the IDP is derived and extended to QPAT, yielding a two-stage reconstruction where the intermediate coefficient $q$ is recovered first and then $(D,\sigma)$ from the recovered $u^{(1)}$, with recommended parameter scaling $h\sim\sqrt{\delta}$ and $\alpha\sim\delta^2$ to achieve $\delta^{1/4}$-type accuracy. Numerical experiments in 2D with both smooth and nonsmooth coefficients validate the theory, show that the nonzero-gradient region expands with more random illuminations, and demonstrate robustness to measurement noise.

Abstract

In this paper, we focus on the numerical analysis of quantitative photoacoustic tomography. Our goal is to reconstruct the optical coefficients, i.e., the diffusion and absorption coefficients, using multiple internal observational data. The foundation of our numerical algorithm lies in solving an inverse diffusivity problem and a direct problem associated with elliptic equations. The stability of the inverse problem depends critically on a non-zero condition in the internal observations, a condition that can be met using randomly chosen boundary excitation data. Utilizing these randomly generated boundary data, we implement an output least squares formulation combined with finite element discretization to solve the inverse problem. In this scenario, we provide a rigorous error estimate in $L^2(Ω)$ norm for the numerical reconstruction using a weighted energy estimate, inspired by the analysis of a newly proposed conditional stability result. The resulting error estimate serves as a valuable guide for selecting appropriate regularization parameters and discretization mesh sizes according to the noise levels present in the data. Several numerical experiments are presented to support our theoretical results and illustrate the effectiveness of our numerical scheme.

Figures (10)

• Figure 1: Boundary illuminations and the non-zero region of Example \ref{['ex:1']}. Top left: plot of boundary data $f^{(\ell)}=g^{(\ell+1)}$. Top middle to bottom right: region where the non-zero condition is satisfied as the number of boundary inputs increases.
• Figure 2: Example \ref{['ex:1']}. First row: reconstructions of $D^\dag$. Second row: reconstructions of $\sigma^\dag$.
• Figure 3: Boundary illuminations and the non-zero region of Example \ref{['ex:2']}. Top left: plot of boundary data $f^{(\ell)}=g^{(\ell+1)}$. Top middle to bottom right: region which satisfying the non-zero condition as number of boundary input increasing.
• Figure 4: Example \ref{['ex:2']}. First row: reconstructions of $D^\dag$. Second row: reconstructions of $\sigma^\dag$.
• Figure 5: Boundary illuminations and the non-zero region of Example \ref{['ex:3']}. Top left: plot of boundary data $f^{(\ell)}=g^{(\ell+1)}$. Top middle to bottom right: region which satisfying the non-zero condition as number of boundary input increasing.

Theorems & Definitions (26)

• Remark 2.1
• Proposition 2.1
• Proof 1
• Theorem 2.1
• Proof 2
• Remark 2.2
• Remark 2.3
• Lemma 2.1
• Proof 3
• Corollary 2.1