Well-posedness for Photoacoustic Tomography with Fabry-Perot Sensors

TL;DR

This paper models photoacoustic tomography measurements using Fabry-Perot sensors, accounting for their directional response, and demonstrates the problem's well-posedness with numerical reconstructions via Landweber iterations.

Contribution

It introduces a realistic mathematical model for Fabry-Perot sensors in photoacoustic imaging and analyzes the problem's well-posedness under this model.

Findings

The photoacoustic problem is well-posed with Fabry-Perot sensors.

Numerical reconstructions are successfully performed using Landweber iterations.

The model incorporates the directional response of sensors, improving realism.

Abstract

In the mathematical analysis of photoacoustic imaging, it is usually assumed that the acoustic pressure (Dirichlet data) is measured on a detection surface. However, actual ultrasound detectors gather data of a different type. In this paper, we propose a more realistic mathematical model of ultrasound measurements acquired by the Fabry--Perot sensor. This modeling incorporates directional response of such sensors. We study the solvability of the resulting photoacoustic tomography problem, concluding that the problem is well-posed under certain assumptions. Numerical reconstructions are implemented using the Landweber iterations, after discretization of the governing equations using the finite element method.

Equations

$$p=p_{\mathrm{s}}\text{and}\frac{1}{\rho }\frac{\partial p}{\partial n}=\frac{1}{{\rho }_{\mathrm{s}}}\frac{\partial p_{\mathrm{s}}}{\partial n}\text{on Γ,}$$

$$p_{\mathrm{s}}=p_{\mathrm{b}}\text{and}\frac{1}{{\rho }_{\mathrm{s}}}\frac{\partial p_{\mathrm{s}}}{\partial n}=\frac{1}{{\rho }_{\mathrm{b}}}\frac{\partial p_{\mathrm{b}}}{\partial n}\text{on Γh,}$$

$$c_{\mathrm{s}}^{-2}\frac{{\partial }^2{p}_{\mathrm{s}}}{\partial t^2}=\frac{{\partial }^2{p}_{\mathrm{s}}}{\partial n^2}+2{\mathcal{H}}_r\frac{\partial p_{\mathrm{s}}}{\partial n}+{\Delta }_{{\Gamma }_r}{p}_{\mathrm{s}},\text{in (0,T)×Ωs.}$$

$$\frac{\partial p_{\mathrm{b}}}{\partial n}={\Lambda }_{\mathrm{b},h}{p}_{\mathrm{b}}\text{on Γh}$$

$${\Lambda }_{\mathrm{b},h}=-c_{\mathrm{b}}^{-1}{\partial }_t-{\mathcal{H}}_h$$

$${\mathcal{H}}_h=\frac{\mathcal{H}+h\mathcal{K}}{1+2h\mathcal{H}+h^2\mathcal{K}}=\mathcal{H}+h\left(\mathcal{K}-2{\mathcal{H}}^2\right)+\mathcal{O}(h^2)$$

$${\Lambda }_{\mathrm{b}}=-c_{\mathrm{b}}^{-1}{\partial }_t-\mathcal{H}$$

$$\frac{1}{\rho }{\frac{\partial p}{\partial n}|}_{\Gamma }$$

$$=\frac{1}{{\rho }_{\mathrm{s}}}\left[\frac{{\rho }_{\mathrm{s}}}{{\rho }_{\mathrm{b}}}{\Lambda }_{\mathrm{b},h}{\left(p_{\mathrm{s}}+h\frac{\partial p_{\mathrm{s}}}{\partial n}\right)|}_{\Gamma }-h{\frac{{\partial }^2{p}_{\mathrm{s}}}{\partial n^2}|}_{\Gamma }\right]+\mathcal{O}(h^2)$$

$$=\frac{1}{{\rho }_{\mathrm{s}}}\left[\frac{{\rho }_{\mathrm{s}}}{{\rho }_{\mathrm{b}}}{\Lambda }_{\mathrm{b},h}{\left(p+h\frac{{\rho }_{\mathrm{s}}}{\rho }\frac{\partial p}{\partial n}\right)|}_{\Gamma }-h{\left(c_{\mathrm{s}}^{-2}{\partial }_t^{2}p-2\mathcal{H}\frac{{\rho }_{\mathrm{s}}}{\rho }\frac{\partial p}{\partial n}-{\Delta }_{\Gamma }p\right)|}_{\Gamma }\right]+\mathcal{O}(h^2)$$

$$\frac{1}{\rho }\frac{\partial p}{\partial n}-\frac{1}{{\rho }_{\mathrm{b}}}{\Lambda }_{\mathrm{b}}p=h\left[\frac{{\rho }_{\mathrm{s}}}{\rho {\rho }_{\mathrm{b}}}{\Lambda }_{\mathrm{b}}\frac{\partial p}{\partial n}+\frac{1}{{\rho }_{\mathrm{b}}}\left(\mathcal{K}-2{\mathcal{H}}^2\right)p-\frac{1}{{\rho }_{\mathrm{s}}}\left(c_{\mathrm{s}}^{-2}{\partial }_t^{2}p-2\mathcal{H}\frac{{\rho }_{\mathrm{s}}}{\rho }\frac{\partial p}{\partial n}-{\Delta }_{\Gamma }p\right)\right]$$

$$\frac{1}{\rho }\frac{\partial p}{\partial n}=\frac{1}{{\rho }_{\mathrm{b}}}{\Lambda }_{\mathrm{b}}p\text{on Γ}.$$

$${\partial }_t^2\mathfrak{m}\sim \left({{\partial }_t^2{\text{u}}_{\mathrm{s}}|}_{\Gamma }-{{\partial }_t^2{\text{u}}_{\mathrm{s}}|}_{{\Gamma }_h}\right)\cdot \text{n}=\frac{1}{{\rho }_{\mathrm{s}}}\left({\frac{\partial p_{\mathrm{s}}}{\partial n}|}_{{\Gamma }_h}-{\frac{\partial p_{\mathrm{s}}}{\partial n}|}_{\Gamma }\right),$$

$$\frac{\partial p_{\mathrm{s}}}{\partial n}$$

$${p_{\mathrm{s}}|}_{{\Gamma }_h}={\left[p_{\mathrm{s}}+h\frac{\partial p_{\mathrm{s}}}{\partial n}\right]|}_{\Gamma }+\mathcal{O}(h^2)={\left[p+h\frac{{\rho }_{\mathrm{s}}}{\rho }\frac{\partial p}{\partial n}\right]|}_{\Gamma }+\mathcal{O}(h^2).$$

$${\partial }_t^2\mathfrak{m}\sim \left(\frac{1}{{\rho }_{\mathrm{b}}}{\Lambda }_{\mathrm{b}}p-\frac{1}{\rho }\frac{\partial p}{\partial n}\right)+h\left[\frac{{\rho }_{\mathrm{s}}}{\rho {\rho }_{\mathrm{b}}}{\Lambda }_{\mathrm{b}}\frac{\partial p}{\partial n}+\left(\mathcal{K}-2{\mathcal{H}}^2\right)p\right]+\mathcal{O}(h^2)$$

$${\partial }_t^2\mathfrak{m}\sim {\partial }_t^{2}p+2\mathcal{H}\frac{{\rho }_{\mathrm{s}}}{{\rho }_{\mathrm{b}}}\frac{c_{\mathrm{s}}^2}{c_{\mathrm{b}}}{\partial }_{t}p+2{\mathcal{H}}^2{c}_{\mathrm{s}}^2\frac{{\rho }_{\mathrm{s}}}{{\rho }_{\mathrm{b}}}p-c_{\mathrm{s}}^2{\Delta }_{\Gamma }p\text{on (0,T)×Γ}$$

$$p(\text{x},t)=e^{i\left(\text{x}\cdot \text{k}-\omega t\right)}+R{e}^{i\left(\text{x}\cdot {\text{k}}_{\mathrm{r}}-\omega t\right)}$$

$$R=\frac{\cos \theta -\alpha }{\cos \theta +\alpha },\text{where}\alpha =\frac{\rho c}{{\rho }_{\mathrm{b}}{c}_{\mathrm{b}}}.$$

$$\frac{\mathfrak{m}}{p_{\mathrm{inc}}}\sim \left(1+\frac{\cos \theta -\alpha }{\cos \theta +\alpha }\right)\left(1-\frac{c_{\mathrm{s}}^2}{c^2}{\sin }^2\theta \right).$$

$$\mathfrak{m}\sim p\text{on (0,T)×Γ,}\text{and ρs=ρb=ρ and cs=cb=c.}$$

$$\mathcal{F}:\left(p_0,p_1\right)\mapsto \mathfrak{m}$$

$${\partial }_t^{2}p$$

$${\rho }_{\mathrm{b}}{\partial }_{n}p$$

$$p=p_0$$

$$\Vert \left(p_0,p_1\right){\Vert }_{H^0(\Omega )\times H^{-1}(\Omega )}\le C\Vert \mathfrak{m}{\Vert }_{H^0((0,T);H^1(\Gamma ))}$$

$${\mathcal{F}}^{\ast }:\psi \mapsto \left({\partial }_t\xi (0),\xi (0)\right)$$

$${\partial }_t^2\xi$$

$${\rho }_{\mathrm{b}}{\partial }_n\xi -\rho {\Lambda }_{\mathrm{b}}^{\ast }\xi$$

$$\xi =0$$

$$\left({\partial }_t^{-2}v\right)(t)={\int }_0^t{\int }_0^{\tau }v(s)dsd\tau .$$

$${\left({\partial }_t^{2}v,\varphi \right)}_{H^0((0,T)\times \Gamma )}={\left(v,{\partial }_t^2\varphi \right)}_{H^0((0,T)\times \Gamma )}$$

$${\mathcal{G}}_1$$

Full text

{AMS}

35R30, 35L05, 35R01, 92C55

1 Introduction

Photoacoustic tomography (PAT) is a hybrid imaging technique based on the photoacoustic effect, which is the transformation of absorbed electromagnetic energy into pressure waves. This technique takes advantage of the fact that absorption exhibits high-contrast in soft biological tissues and that acoustic waves can be measured with broadband transducers leading to imaging with high-resolution. Therefore, high-contrast and high-resolution can be achieved simultaneously [14, 22, 23, 58, 65, 69, 70, 71, 72].

For qualitative photoacoustic tomography, the goal is to form an image of the initial state of the pressure field using boundary measurements of the transient pressure waves. Most of the reconstruction methods assume that the actual pressure field (Dirichlet data) can be measured at the boundary [2, 3, 4, 5, 11, 22, 32, 35, 36, 37, 38, 42, 44, 51, 52, 55, 59, 60, 61, 63]. In reality, ultrasound sensors are not able to measure the pressure field directly. Instead, they measure certain combinations of the field and its derivatives. This challenge has been noted in [77] and investigated by Finch [31] and by Zangerl, Moon and Haltmeier [79].

Fabry–Perot transducers offer an alternative to piezoelectric sensors for ultrasound-based imaging applications [10, 16, 24, 34, 62, 78, 80]. The design consists of a sensing film (10-50 $\mu$m thick) sandwiched between extremely thin optically reflective coatings ($\approx$ 50 nm thick) lying on an optically transparent backing substrate ($\approx$ 2 cm thick). An illustration is shown in Fig. 1. An interrogating laser beam is employed to generate reflections from both optically reflective coatings. When an incident pressure wave modulates the thickness of the sensing material, the change in the interference pattern from the reflected laser beam is used to estimate the distance between the reflective coatings. The deformation of the sensing material can then be related to pressure measurements. Cox and Beard provide a description of the Fabry–Perot design, and an excellent study of its frequency and directional responses [24].

In this paper, we model and investigate the mathematical solvability of the PAT problem for measurements acquired by sensors based on the Fabry–Perot design. We model idealized point-like ultrasound transducers and the physical variable being measured by these sensors. Our goal is to determine whether such measurements lead to the mathematical solvability of the PAT problem. We shall not account for resolution limitations of the Fabry–Perot design due to finite-size sensing elements. We refer the reader to [24, 66, 74, 75] for investigations concerning this issue. Our modeling is further simplified by ignoring shear waves that can travel in the sensing material and its backing substrate. In other words, we develop an analysis based entirely on the scalar wave equation.

The acoustic domain $\Omega$ contains soft tissue with density $\rho$ and wave speed $c$. The sensing film $\Omega_{\rm s}$ of thickness $h>0$ has density $\rho_{\rm s}$ and wave speed $c_{\rm s}$. The backing substrate $\Omega_{\rm b}$ has density $\rho_{\rm b}$ and wave speed $c_{\rm b}$. We assume that $\rho_{\rm s}$, $\rho_{\rm b}$, $c_{\rm s}$ and $c_{\rm b}$ are positive constants. However, $\rho$ and $c$ may vary within $\Omega$. The interface between the acoustic domain and the sensing film is denoted $\Gamma$. The interface between the sensing film and the backing substrate is denoted $\Gamma_{h}$. Typically, the sensing film and the backing substrate are acoustically more rigid than the biological soft tissue of interest. Hence, the presence of the sensors induces partial reflections of the waves. Other researchers have investigated PAT with reflecting boundaries assuming that the actual pressure can be measured [2, 21, 28, 39, 45, 50]. Here we seek to incorporate in our model the influence that the sensor exerts on the pressure waves, as well as the nature of the acoustic measurements for the Fabry–Perot design. The interplay between the sensors and the pressure field is modeled by the following transmission conditions at the interface $\Gamma$,

[TABLE]

where $p$ and $p_{\rm s}$ are the pressure in the acoustic medium and sensing material, respectively. The first condition in (1), known as dynamic transmission, ensures the continuity of the pressure field. The second condition in (1), known as kinematic transmission, ensures the continuity of particle motion in the normal direction. Similar transmission conditions hold at the interface $\Gamma_{h}$,

[TABLE]

where $p_{\rm b}$ is the pressure in the backing substrate.

In order to simply the analysis, in Section 2 we derive an effective boundary condition (valid for small $h>0$) to replace the transmission conditions (1)-(2). In Section 3 we mathematically model the measurements acquired by ultrasound sensors based on the Fabry–Perot design. For the effective boundary condition and modeled boundary measurements, in Section 4 we state and prove the solvability of the photoacoustic tomography problem. A reconstruction algorithm is proposed in Section 5 where some numerical experiments are presented as well. The conclusions follow in Section 6.

2 Effective boundary condition

For analytical and numerical purposes, it is convenient to replace the transmission conditions (1)-(2) for an asymptotically equivalent boundary condition for the acoustic pressure field at the boundary $\Gamma$. This condition is meant to account for the transmission into the sensing film $\Omega_{\rm s}$ and into the backing substrate $\Omega_{\rm b}$ without having to explicitly solve for the wave fields in those domains. See [7, 8, 18, 19, 41, 57] where similar problems are treated. This effective boundary condition also simplifies the model for the measurements as shown in Section 3.

We proceed by making some geometric assumptions about the domain $\Omega_{\rm s}$ occupied by the sensing film. We use the concept of parallel surfaces to define the shape of this extremely thin layer of material. These surfaces are parametrized by $0<r<h$ and defined by $\Gamma_{r}=\left\{\textbf{y}=\textbf{x}+r\textbf{n}(\textbf{x}):\textbf{x}\in\Gamma\right\}$. For smooth $\Gamma$ and sufficiently small $h$, each surface $\Gamma_{r}$ is well-defined and smooth. Moreover, the normal vector $\textbf{n}(\textbf{x}+r\textbf{n}(\textbf{x}))$ of the parallel surface $\Gamma_{r}$ coincides with the normal vector $\textbf{n}(\textbf{x})$ of $\Gamma$ for each $\textbf{x}\in\Gamma$. We let $\Omega_{\rm s}$ be the union of this family of parallel surfaces, where $h$ being sufficiently small ensures that each point $\textbf{y}\in\Omega_{\rm s}$ can be uniquely represented in the form $\textbf{y}=\textbf{x}+r\textbf{n}(\textbf{x})$ for $\textbf{x}\in\Gamma$ and $0<r<h$. See details in [43, §6.2] and [26, Probl. 11 §3.5].

The pressure field $p_{\rm s}$ in the sensing material satisfies the wave equation,

[TABLE]

For convenience, we have expressed the Laplacian in $\Omega_{\rm s}$ using the normal derivative $\partial/\partial n$ (which makes sense at any point in $\Omega_{s}$ given its definition in terms of parallel surfaces), the mean curvature $\mathcal{H}_{r}$ of $\Gamma_{r}$ and the Laplace–Beltrami operator $\Delta_{\Gamma_{r}}$ associated with $\Gamma_{r}$. See details in [9]. As in [24], we assume that the pressure field $p_{\rm b}$ in the backing substrate is outgoing. Therefore, the pressure field $p_{\rm b}$ satisfied the following radiation condition

[TABLE]

where $\Lambda_{{\rm b},h}$ is a nonreflecting boundary operator. The subject of nonreflecting or absorbing boundary conditions is beyond the scope of this paper. We refer to [9, 13, 6, 20, 1] for some relevant articles on that topic. We consider

[TABLE]

which is derived in [9] as a first order nonreflecting condition that takes into account the mean curvature $\mathcal{H}_{h}$ of the boundary $\Gamma_{h}$. As shown in [43, §6.2] or [26, Probl. 11 §3.5], the mean curvatures $\mathcal{H}$ and $\mathcal{H}_{h}$ of the surfaces $\Gamma$ and $\Gamma_{h}$, respectively, are related by

[TABLE]

where $\mathcal{K}$ is the Gaussian curvature of $\Gamma$. Notice that we are using the mean curvature sign convention from [9], not from [43] or [26]. Therefore, we have that at the surface $\Gamma$, the associated nonreflecting operator $\Lambda_{\rm b}$ given by

[TABLE]

satisfies $\Lambda_{{\rm b},h}=\Lambda_{\rm b}+h\left(\mathcal{K}-2\mathcal{H}^{2}\right)+\mathcal{O}(h^{2})$ where $\Lambda_{{\rm b},h}$ is defined in (5).

We proceed with a Taylor expansion for the normal derivative of the pressure field,

[TABLE]

where we have employed the transmission conditions (1)-(2), the wave equation (3) and the radiation condition (4). Using (6), re-grouping terms and neglecting $\mathcal{O}(h^{2})$ terms, we obtain a first order boundary condition

[TABLE]

for the acoustic wave field $p$ on the surface $\Gamma$. This conditions implies that most of the influence (zeroth order terms) that the sensor exerts on the pressure field at the boundary is provided by the thick backing substrate which reflects and refracts the wave field according to the mismatch in densities $\rho$ and $\rho_{\rm b}$, and in wave speeds $c$ and $c_{\rm b}$. The presence of the thin sensing film is accounted for by the first order terms in (7). If these latter terms are neglected, the pressure field satisfies the following zeroth order effective boundary condition

[TABLE]

However, the terms of order $\mathcal{O}(h)$ in (7) become important in Section 3 where we model the ultrasound measurements which are of order $\mathcal{O}(h)$.

3 Modeling ultrasound measurements

For an ultrasound sensor based on the Fabry–Perot design, the quantity being measured is proportional to the difference in the normal projection of the particle displacement on both sides of the sensing film [24]. Hence, up to a constant of proportionality, the measurements $\mathfrak{m}$ acquired by the ultrasound transducer satisfy

[TABLE]

where the symbol $\sim$ means equality up to a multiplicative constant, and the pressure–displacement formulation is valid in the absence of shear stress. We seek to express the measurement in terms of the pressure in the acoustic medium only. Using the transmission conditions (1)-(2) at both sides of the sensing film, and (4)-(6), we obtain

[TABLE]

We make the following Taylor approximation for the pressure field $p_{\rm s}$ within the sensing film and combine it with (1)-(2) and (10) to obtain,

[TABLE]

Therefore, combining (9)-(11) we obtain an expression for the measurements in terms of the acoustic pressure field (Dirichlet data) and its normal derivative (Neumann data) at the boundary $\Gamma$ as follows,

[TABLE]

Neglecting the $\mathcal{O}(h^{2})$ terms on the right-hand side of (12) and using the effective boundary conditions (7)-(8), we obtain a simplified or first order model for the measurements,

[TABLE]

where we have used the definition of the operator $\Lambda_{\rm b}$ given by (6). In order to fully determine $\mathfrak{m}$, initial conditions must be provided. In consistency with the PAT scenario, where the pressure field has vanishing initial Cauchy data in the exterior of $\Omega$, we let $\mathfrak{m}=\partial_{t}\mathfrak{m}=0$ on $\{t=0\}\times\Gamma$.

To illustrate the response associated with this sensor design, we briefly analyze its behavior for plane waves. For convenience, we momentarily assume that $\Gamma$ is a plane through the origin. Both, the effective boundary condition (8) and the form of the measurements (13) play a role in this analysis. From (6) we have $\Lambda_{\rm b}=-c_{\rm b}^{-1}\partial_{t}$ because for a flat surface $\Gamma$ the mean curvature is $\mathcal{H}=0$. A plane wave $p_{\rm inc}=e^{i\left(\textbf{x}\cdot\textbf{k}-\omega t\right)}$ with incidence angle $\theta$, induces a reflection governed by the effective boundary condition (8). The superposition of the incident and reflected waves has the form

[TABLE]

where $R$ is the reflection coefficient, $\textbf{k}_{\rm r}$ is the reflection wavenumber, such that $|\textbf{k}|=|\textbf{k}_{\rm r}|=\omega/c$ and $\textbf{n}\cdot\textbf{k}_{\rm r}=-\textbf{n}\cdot\textbf{k}$ where n is the outward normal vector. We also have $\textbf{n}\cdot\textbf{k}=|\textbf{k}|\cos\theta$. Once (14) is plugged into (8), the reflection coefficient is shown to satisfy

[TABLE]

Plugging (14)-(15) into (13) and evaluating at the origin $\textbf{x}=\textbf{0}$, we find that the measurements satisfy

[TABLE]

Figure 2 displays the directional response (16) for plane waves as a function of the incidence angle $\theta$. The parameters were taken from [24] for a Parylene sensing film and polycarbonate backing substrate. Figure 2 shows that incidence at approximately $42.99^{\circ}$ corresponds to a critical angle where the waves cause no particle motion in the normal direction. This occurs when the tangential phase speed of the acoustic wave equals the wave speed in the sensing film. The Fabry–Perot sensor does not capture such acoustic waves. We also note that for incidence angles greater than this critical angle, the pressure wave and the measurement have opposite signs. The response also approaches zero as the incidence angle approaches $90^{\circ}$. Therefore, for incidence at tangential angles, the sensor design is not able to measure the pressure waves adequately.

The above are some of the features not taken into account when it is assumed that the actual pressure field (Dirichlet data) is acquired at the measurement boundary. As explained in the Introduction, the overly idealized assumption for ultrasound sensors is that they measure the pressure field at the boundary and that the sensors themselves do not perturb the pressure waves. This overly idealized model can be expressed as follows,

[TABLE]

4 Main mathematical results

Here we define the photoacoustic tomography problem in terms of the wave equation, the effective boundary condition (8) and the ultrasound measurements modeled by (13). We also prove the solvability of this problem under the following geometric condition for the wave speed $c$ and the domain $\Omega$.

{assumption}

[Non–trapping Condition] Let $\Omega$ be a simply-connected bounded domain with smooth boundary $\Gamma$. Assume there exists $T_{\rm o}<\infty$ such that any geodesic ray of the manifold $(\Omega,c^{-2}dx^{2})$, originating from any point in $\Omega$ at time $t=0$ reached the boundary $\Gamma$ at a non-diffractive point before $t=T_{\rm o}$.

The forward mapping, which we seek to invert, is given by

[TABLE]

where the measurement $\mathfrak{m}$ satisfies (13) with vanishing initial Cauchy data. The initial velocity $p_{1}$ is known to be zero in the context of PAT. However, the mathematical theory allows us to recover it as well. The pressure field $p$ solves the following initial boundary value problem,

[TABLE]

The well-posedness of this problem for $\left(p_{0},p_{1}\right)\in H^{1}_{0}(\Omega)\times H^{0}(\Omega)$ has been established. See for instance [30, 46, 47]. The unique solution satisfies $p\in C([0,T];H^{1}(\Omega))$, $\partial_{t}p\in C([0,T];H^{0}(\Omega))$, $p|_{\Gamma}\in H^{1}((0,T)\times\Gamma)$ and $\partial_{n}p\in H^{0}((0,T)\times\Gamma)$. We work with the standard Sobolev spaces based on square-integrable functions over $\Omega$ or $(0,T)\times\Gamma$. The associated inner-product extends naturally as the duality pairing between functionals and functions. We should interpret the Hilbert space $H^{0}(\Omega)$ with the inner-product appropriately weighted by $c^{-2}$ so that $c^{2}\Delta$ is formally self-adjoint with respect to the duality pairing of $H^{0}(\Omega)$.

Now we state our main result in the form of a theorem.

Theorem 4.1** (Main Result).**

Under the non–trapping Section 4 for the manifold $(\Omega,c^{-2}dx^{2})$ and time $T>T_{\rm o}$, the forward mapping $\mathcal{F}:H^{1}_{0}(\Omega)\times H^{0}(\Omega)\to H^{0}((0,T);H^{1}(\Gamma))$ is injective, that is, the photoacoustic tomography problem is uniquely solvable. Moreover, the following stability estimate holds

[TABLE]

*for some constant $C>0$. *

We proceed to prove this theorem by showing that the adjoint of the forward mapping is surjective. This adjoint mapping $\mathcal{F}^{*}$ is given by

[TABLE]

where $\xi$ solves the following backwards–in–time boundary value problem,

[TABLE]

where $a=2\mathcal{H}c_{\rm s}^{2}\rho_{\rm s}/(c_{\rm b}\rho_{\rm b})$ and $b=2\mathcal{H}^{2}c_{\rm s}^{2}\rho_{\rm s}/\rho_{\rm b}$ are constants. The operator $\partial_{t}^{-2}$ can be defined as

[TABLE]

Notice that $\partial_{t}^{-2}\partial_{t}^{2}v=\partial_{t}^{2}\partial_{t}^{-2}v=v$ for all sufficiently smooth $v$ such that $v|_{t=0}=\partial_{t}v|_{t=0}=0$. Also notice that $\partial_{t}^{2}$ is formally self–adjoint. In particular,

[TABLE]

for all sufficiently regular $v$ and $\phi$ such that $v|_{t=0}=\partial_{t}v|_{t=0}=0$ and $\phi|_{t=T}=\partial_{t}\phi|_{t=T}=0$. The Laplace–Beltrami operator $\Delta_{\Gamma}$ is also self–adjoint since the manifold $\Gamma$ has no boundary. The nonreflecting operator $\Lambda_{\rm b}$ defined in (6) has an adjoint given by $\Lambda_{\rm b}^{*}=c_{\rm b}^{-1}\partial_{t}-\mathcal{H}$. Therefore, it stays as a differential operator with first and zeroth order terms.

The statement of Theorem 4.1 is a direct consequence of the following lemma.

Lemma 4.2**.**

*Under the non–trapping Section 4 for the manifold $(\Omega,c^{-2}dx^{2})$ and time $T>T_{\rm o}$, the operator $\mathcal{F}^{*}:H^{0}((0,T);H^{-1}(\Gamma))\to H^{0}(\Omega)\times H^{1}(\Omega)$ is surjective. *

Proof 4.3**.**

The mapping $\mathcal{F}^{*}$ can be composed as $\mathcal{F}^{*}=\mathcal{G}_{2}\circ\mathcal{G}_{1}$ where

[TABLE]

where the mapping $\mathcal{G}_{1}$ is given by

[TABLE]

where $\Psi$ has vanishing Cauchy data at $t=T$ and solves $\partial_{t}^{2}\Psi=\rho_{\rm b}\psi$. The mapping $\mathcal{G}_{2}$ is defined via $\xi$, the solution of

[TABLE]

Under the non–trapping assumption, the mapping $\mathcal{G}_{2}:\varphi\mapsto\left(\partial_{t}\xi(0),\xi(0)\right)$ defined by (27), is well–known to be surjective from $H^{0}((0,T);H^{0}(\Gamma))$ onto $H^{0}(\Omega)\times H^{1}(\Omega)$. That is the central theme of exact boundary controllability for the wave equation. See [33, Ch 6] and [12, Corollary 4.10].

*Hence, it only remains to show that the mapping $\mathcal{G}_{1}$ is surjective. This can be accomplished by proving that equation (26) is solvable for any forcing term $\varphi\in H^{0}((0,T);H^{0}(\Gamma))$ such that $\Psi$ has vanishing Cauchy data at $t=T$ and $\partial_{t}^{2}\Psi\in H^{0}((0,T);H^{-1}(\Gamma))$. This solvability is a well–established result. See [30, Theorems 3-5 in §7.2], [46] and [47, Ch 3 §8, Thm 8.1]. As a consequence, the operator $\mathcal{F}^{*}$ is surjective from $H^{0}((0,T);H^{-1}(\Gamma))$ onto $H^{0}(\Omega)\times H^{1}(\Omega)$. *

This Lemma 4.2 renders the proof of Theorem 4.1. Indeed, we first notice that since $\mathcal{F}^{*}:H^{0}((0,T);H^{-1}(\Gamma))\to H^{0}(\Omega)\times H^{1}(\Omega)$ is well-defined and surjective, then $\mathcal{F}:H^{0}(\Omega)\times H^{-1}(\Omega)\to H^{0}((0,T);H^{1}(\Gamma))$ is well-defined, injective, and has a closed range. The stability estimate (20) then follows from the Open Mapping Theorem (see [48, Ch 2] or [27, Ch 2]).

5 Numerical results

In this section we implement reconstruction algorithms to the solve the PAT problem at the discrete level. The reconstructions are based on the Landweber iterative method [29, Ch. 6]. In the context of PAT, the Landweber iteration has been employed previously because of its simplicity and compatibility with regularization methods [17, 25, 35, 36, 37, 52, 64]. Other approaches, such as the conjugate gradient method, may also be employed to solve PAT problems [2, 3, 37, 49, 52, 66, 67, 68, 76]. The Landweber iteration is defined in Algorithm 1.

The discretizations of the forward map (18) and adjoint map (21) are based on a piecewise linear finite element method (FEM) and second order finite difference for the time derivatives in the initial boundary value problems (19) and (22), respectively. The discretization parameters were chosen to satisfy the CFL stability condition. The FEM was implemented on triangular meshes of the unit disk and the physical parameters, described in Fig. 2, were non-dimensionalized accordingly. Figure 3 shows a coarse triangular mesh and the Shepp–Logan phantom to be reconstructed. Measured data was synthetically generated by discretizing the forward operator $\mathcal{F}$ using the FEM method. In all simulations, the mesh employed to generate the measurements had mesh size about half of the mesh size for the mesh employed to reconstruct the phantom.

Figure 4 displays the error $\|\phi_{k}-p_{0}\|$ (where $p_{0}$ is the exact solution) as a function of the iteration number $k\geq 0$ of the Landweber method. The figure shows results for three FEM meshes that were consecutively refined by halving the mesh size. We notice that initially, the error decays exponentially in $k$ (as the theory of this method predicts) but then it stagnates. The stagnation level decreases with mesh refinement. This phenomenon may be attributed to the fact that the discrete version of $\mathcal{F}^{*}$ is not the actual adjoint of the discrete version of $\mathcal{F}$. Thus, the discretization of the normal operator $\mathcal{F}^{*}\mathcal{F}$ is not symmetric positive definite as this method requires. However, as the mesh is refined, we expect this error to reduce. Implementation of an exact numerical adjoint, as done by Huang et al. [40], could remedy this issue.

Lastly, we compare the reconstruction of the initial pressure profile obtained by accounting for the structure of the Fabry–Perot measurement model (13) against the reconstruction obtained from the overly idealized (but commonly assumed) model (17). The latter reconstruction is obtained by synthetically producing the measurements following the model (13), but then incorrectly assuming that the measurements satisfy (17). Figure 5 displays the reconstruction results for both measurement models using 60 iterations of the Landweber method. For the reconstruction following the proposed model (13), the relative error is $0.73\%$. For the reconstruction following the overly idealized model (17), the relative error is $23.02\%$.

6 Conclusions

We have developed a model for the type of measurements acquired by sensors based on the Fabry–Perot design. This model takes the form shown in (13). The validity of this expression is limited to small values for the thickness $h$ of the sensing film with respect to the wavelength of the pressure fields. This means that $h\ll c_{\rm s}/f$ where $f$ is the oscillatory frequency. For instance, [24] considered a Fabry–Perot polymer film of thickness $h=40$ $\mu$m with compressional wave speed $c_{\rm s}=2200$ m/s and a frequency range $1-15$ MHz. At the higher end of this range, the film thickness is about one fourth of the wavelength. Therefore, the proposed model would be valid for most of this frequency range.

Our mathematical model of the measurements captures the directional response observed experimentally [24, 34, 62, 78, 80]. For instance, notice that for a pressure wave $p$ impinging the boundary $\Gamma$ in the normal direction, the sensor design fully captures the pressure field. However, for pressure waves traveling at other incidence angles, the sensor response may exhibit non-ideal behavior, such as vanishing response at critical angles, as shown in Figure 2. This is the mathematical description of the directivity associated with these ultrasound sensors. The incorporation of these features into reconstruction algorithms has been recognized as one of the challenges associated with improving photoacoustic inversion [22, 28, 66, 74].

Using the model (13) for the measurements, we studied the solvability of the PAT problem and concluded that the problem is well-posed in the appropriate spaces and norms. See the precise statements in Theorem 4.1. Following the analysis, a reconstruction algorithm was implemented based on the Landweber iteration. We carried out proof-of-concept numerical simulations to illustrate the reconstructions obtained from this method for synthetic data after discretization using FEM. For the chosen Shepp–Logan phantom, Figure 5 displays the results obtained by incorporating (Panel A) and by ignoring (Panel B) the model for the Fabry–Perot measurements. The respective error profiles are shown in Panels C and D of the same figure. Approximately, a $22\%$ relative error is added when the proposed model for the Fabry–Perot measurement is not incorporated in the reconstruction algorithm. We also highlight the qualitative difference between the error profiles from Panels C and D of Figure 5. By ignoring the Fabry–Perot model, the error profile exhibits prominent artifacts over the entire image, especially near the detection boundary. These artifacts may be explained by the directivity response shown in Figure 2. By design, the proposed reconstruction algorithm accounts for the directivity response of these sensors leading to the removal of those artifacts.

Finally, we propose a couple of directions for future research that may improve or extend this work. It remains to study the well-posedness of the PAT problem for a Fabry–Perot measurement model that includes both the p-waves and s-waves in the elastic sensing film and backing substrate of the sensor. Such a model would incorporate the influence of the shear waves on the measurements as studied by Cox and Beard for plane waves [24]. Also, non-trivial directivity responses are not only induced by the Fabry–Perot sensor design, but also by piezoelectric detectors [15, 53, 54, 56, 73]. Therefore, analysis of the well-posedness for the PAT problem using piezoelectric measurements is also needed.

Acknowledgments

The author would like to thank Texas Children’s Hospital for its support and for the research-oriented environment provided by the Predictive Analytics Lab.