Photoacoustic or thermoacoustic tomography is a non-ionizing imaging modality designed to advantageously combine the high contrast of optical absorption with the high resolution from broadband ultrasound waves. The imaging of optical absorption reveals important functional and pathological information about biological tissues Wang2003; Beard2011; Wang-Anastasio-2011; Wang2012.
One of the open challenges concerning photoacoustic inversion is the incorporation of realistic models for acoustic measurements. This need for modeling the physics of ultrasound sensors has been recognized in Johansson2003; Finch2005p; Wang2007; Acosta-Montalto-2015. It has been claimed that ultrasound measurements can be described as a linear combination of the pressure field and its normal derivative at the boundary. With that motivation, Dreier and Haltmeier Dreier2019 recently established explicit formulas for the inversion of the two-dimensional wave equation from Neumann boundary data for circular and elliptical domains. In a related effort, Zangerl, Moon and Haltmeier Zangerl2019 derived Fourier-based reconstruction formulas for the spherical detection geometry from knowledge of Robin boundary data.
Most other reconstruction algorithms assume that the pressure field (Dirichlet data) is measured on the detection boundary. This assumption is not satisfied in practice because acoustic transducers do not measure pressure directly. They measure a surrogate for pressure that depends on the actual transducer mechanism. The most common mechanisms for ultrasound applications are based on the piezoelectric effect Wilkens2007; Wear2018b; Wear2019a, on Fabry–Perot interferometry Beard1999b; Cox2007c; Zhang2008; Guggenheim2017 or on fiber-optic refractometry Wild2008; Wear2018a; Wissmeyer2018. In Acosta2019a we formulated a model specifically tailored to the Fabry–Perot sensor design. In this paper, we derive a similar model for piezoelectric sensors and analyze the well-posedness of photoacoustic imaging with such measurements. In order to attain a balance between accuracy and simplicity, the model we develop here is based on the following underlying idealizations:
The sensors are treated as point-like detectors. Hence, we do not account for resolution limitations due to the finite size of sensing elements. See Xu2003c; Haltmeier2010; Wang2011a; Roitner2014; Wear2018b; Wear2019a for investigations concerning this issue. We also assume that the domain to be imaged is fully enclosed by the detection surface.
We assume that in each detector, the piezoelectric film is flat and its thickness is small in comparison to the wavelengths under consideration. In practice, industrial processes can manufacture piezoelectric films with thickness 30 – 100 m approximately Fay1994; Sirohi2000; Su2005; Wilkens2007; Oreilly2010.
Although the sensors contain elastic materials that may support shear waves, our analysis is entirely based on compressional waves governed by the scalar wave equation.
For the piezoelectric film, the poling direction is along its thickness, and the piezoelectric properties are transversely isotropic in the plane perpendicular to the poling direction. The sensing film is mechanically isotropic.
Sensors may have complex structures, including a casing for structural integrity, electrodes, bonding layers and multiple paddings designed to match the mechanical impedance of the acoustic mediumBrown2000; Tichy2010; Szabo2014. However, we assume a simple design consisting of the piezoelectric film, sandwiched by electrode foils of negligible thickness, mounted on a much thicker backing layer. This follows models described in Fay1994; Sirohi2000; Szabo2014.
An illustration of the idealized setup is shown in Figure 1. The acoustic domain, denoted by , contains soft tissue with variable density and variable wave speed . The piezoelectric film has a small uniform thickness , constant density and constant wave speed . The thick backing layer has constant density and constant wave speed . The interface between the acoustic domain and the piezoelectric film is denoted . The interface between the piezoelectric film and the backing layer is denoted .
In section 2 we derive a model for the transduction from pressure to electrical voltage which is the physical quantity acquired by the piezoelectric sensors. In section 3 we derive an effective boundary condition for the transmission of waves from the acoustic medium of interest into the piezoelectric sensor. This effective transmission condition accounts for the influence that the sensor exerts on the acoustic waves. Using the coupled models for piezoelectric measurements and wave propagation, in section 4 we define the forward problem associated with photoacoustic imaging. Then in section 5 we state and prove the solvability of the imaging problem with piezoelectric measurements. A reconstruction algorithm is proposed in section 6 where some numerical simulations are presented. The conclusions follow in section 7.
2 Model for Piezoelectric Measurements
We start the modeling of the piezoelectric measurements from the basic constitutive relations for both piezoelectric and mechanical variables. Since the sensing material is mechanically isotropic, the stress tensoris related to the strain tensor as follows,
where the direction along the thickness of the piezoelectric film is denoted as the 3-axis, and the 1-axis and 2-axis are the transverse plane. Here is the Kronecker delta, and and are the first and second Lamé coefficients. The equation of mechanical motion is
is the material displacement vector. For irrotational deformations, i.e. in the absence of shear stress, the above equation can be simplified in order to relate the particle displacementu to the pressure in the piezoelectric film,
where the pressure is defined as
where the wave speed is defined by .
The piezoelectric transducer measures the electrical voltage across the piezoelectric film generated by the mechanical deformation due to the transmitted acoustic waves. We proceed to derive the mathematical relationship between the voltage and the pressure in the piezoelectric material. Our guiding references are (Tichy2010, Ch. 5), (Szabo2014, Ch. 5) and Sirohi2000. Under small perturbations of field conditions, the linearized constitutive relation for the piezoelectric effect is the following
where D is the electric displacement (electric charge per area), E is an externally applied electric field (voltage per length), is the stress tensor (force per area). The piezoelectric properties are the dielectric permittivity tensor (capacitance per length) and the piezoelectric tensor (electric charge per force). In the absence of shear stress and of an external electric field, the normal electric displacement is given by
As assumed above, the piezoelectric tensor in transversely isotropic, which allows us to simplify the notation as follows , , and . Combining the constitutive relations (7) and (1) we obtain the electric displacement in terms of the strain,
where and . As a consequence, using the definition of strain s in terms of the displacement u, we obtain
We have expressed where represents the surface Laplacian on the transverse plane.
The voltage generated across the film by the electric displacement , which is obtained from Gauss law
where is the thickness of the piezoelectric film, is its dielectric permittivity, and is the differential across the film. Combining (11) and (12), we obtain our model for the piezoelectric measurements
with vanishing initial state. Here the symbol
denotes equality up to a multiplicative constant (which is typically estimated through experimental calibration). In (13), it is assumed that the pressure field is constant across the piezoelectric film. This assumption is rigorously justified in the next section. The model (13) for the piezoelectric sensor is qualitatively similar to the model for the Fabry–Perot transducer proposed in Acosta2019a in spite of the completely different physical principles from which they are derived.
The symbol appearing in the model (13) is a unitless coefficient defined by the elastic and piezoelectric properties of the sensing film
where we have expressed and in terms of Poisson’s ratio to obtain the last equality. Common values for all these physical parameters are shown in Table 1 for polyvinylidene fluoride (PVDF) piezoelectric sensors.
We note from (13) that a theoretically perfect transduction from pressure to voltage would be attained if the coefficient . However, due to the nature of the poling processes employed to manufacture these piezoelectric materials, the coefficients and have opposite signs and generally . This implies that . Hence, in order for , the Poisson’s ratio would have to be which requires the piezoelectric material to be incompressible. In practice, Poisson’s ratio for PVDF films ranges from 0.2 to 0.4 approximately. We note that ranges from 0.3 to 1.5, for the realistic range of values for the Poisson’s ratio and the piezoelectric ratio displayed in Table 1.
|PVDF thickness||10 – 60||m|
|PVDF density||1780 – 1950||kg m|
|PVDF wave speed||1300 – 2300||m s|
|PVDF Poisson’s ratio||0.2 – 0.4|
|Piezoelectric coeff.||-(30 – 35)||pC/N|
|Piezoelectric coeff.||3 – 15||pC/N|
|Coefficient||0.3 – 1.5|
|Backing density||1900 – 2500||kg m|
|Backing wave speed||1000 – 4000||m s|
3 Effective Model for Wave Propagation
Typically, the piezoelectric film and the backing layer are acoustically more rigid and heavier than the biological medium of interest. Therefore, the presence of the sensors induces partial reflections of the waves. Here we seek to model how the sensors exert influence on the pressure waves. This model takes the form of an effective impedance boundary condition that replaces the more involved transmission process for waves traveling from the acoustic domain , through the piezoelectric film and into the backing layer . We assume that the pressure field in the backing layer is outgoing which translates into satisfying a radiation condition of the form,
where is the mean curvature of the surface . See Antoine1999; Acosta2017c for a derivation.
As in Acosta2019a, we make some geometric assumptions about the domain occupied by the piezoelectric film. We let . For sufficiently small , the domain can be expressed as a family of parallel surfaces parametrized by and defined by where is the normal vector at . For smooth and sufficiently small , the surfaces are well-defined, smooth and mutually disjoint. Each point can be uniquely represented in the form for and . In addition, the normal vector at coincides with the normal vector at . See details concerning parallel surfaces in (Kress-Book-1999, Sect. 6.2).
The transmission of the pressure field from the acoustic domain into the piezoelectric film is governed by the following transmission conditions at the interface ,
where and are the pressure in the acoustic medium and piezoelectric film, respectively. The first condition in (16) ensures the continuity of the pressure field. The second condition in (16) ensures the continuity of particle motion in the normal direction. Similar transmission conditions hold at the interface ,
where is the pressure in the backing layer. The pressures , and satisfy the wave equation with respective wave speeds , and .
Now, we proceed to use the method of matched asymptotic expansions to derive an effective model for the interplay between the pressure fields and the piezoelectric sensor. First we consider the formal asymptotic expansions for the pressure fields,
and introduce a change of variable in order to extract the effect of the piezoelectric film thickness ,
The boundary value problem for the pressure field in the piezoelectric film, is governed by the wave equation (5) and the transmission conditions (16)-(17) can be recast in terms of and terms with same powers of are gathered to obtain the following cases.
which imply that is constant as a function of and that the first effective transmission condition is that
which imply that is constant as a function of , and that the second effective transmission condition is
Similar models for photoacoustics are studied in Acosta-Montalto-2015; Acosta-Montalto-2016; Acosta2019a; Acosta2019b
As an example for the response of the piezoelectric sensor design, we can analyze its behavior for plane waves and for a flat boundary . Both the boundary value problem (22)-(23) and the model for the measurements (13) play an important role in this analysis. A plane wave of the form propagating in the direction of k, induces a reflection governed by (23). The total pressure field is the superposition of the incident and reflected wave,
where is the reflection coefficient, is the reflection wavenumber satisfying and , where n is the outward normal on . We can write where is the angle of incidence. Plugging (24) into (23) and neglecting the terms, we find that the reflection coefficient satisfies
where is given by (14). Figure 2 displays the directional response (26) in dB as a function of the incidence angle and piezoelectric coefficient over a realistic range of values shown in Table 1. We observe that for values of , a critical angle appears. For incidence at this critical angle, vanishing measurements are obtained by the piezoelectric sensor design. This critical angle is given by .
4 The Forward Problem
Now we proceed to define the inverse problem of photoacoustic imaging in terms of the wave propagation model (22)-(23) and the model for piezoelectric measurements (13). We neglect higher order terms studied in the previous section, so that the pressure field is assumed to satisfy the following initial value problem,
where is the measurement time to be determined later. Recall that the underlying assumption concerning media properties are that is bounded from below and above, and is smooth in , and that , , , , and are constants. The forward mapping, which we seek to invert for photoacoustic imaging, is given by
where, according to the piezoelectric model (13), the measured electric voltage satisfies
for the pressure field evolving according to (27) from the initial condition .
We work with the standard Sobolev spaces based on square-integrable functions defined on the domain or the boundary . The associated inner product extends as the duality pairing between functionals and functions. For the Sobolev space , the inner product is weighted by so that the differential operator is formally self-adjoint. The well-posedness in Sobolev spaces of the initial value problem (27) is a well-established result EvansPDE; Lio-Mag-Book-1972.
5 The Inverse Problem
The inverse problem associated with photoacoustic imaging is the following: Given the voltage measurements modeled by (29) on , induced by the pressure field satisfying (27), find the unknown initial condition . The solvability of this inverse problem depends on the geometry of the domain , the profile of the wave speed and the time . These conditions are made precise in the following assumption, known as the geometric control condition or a nontrapping condition for the geodesic flow Bardos1992; GlowinskiLionsHe2008.
Assumption 1 (Nontrapping condition)
Let be a simply connected bounded domain with smooth boundary . Assume there exists such that any geodesic ray of the manifold , originating from any point in at time , reaches the boundary at a nondiffractive point before .
With this assumption in place, we can state the main result of the paper in the form of a theorem.
Under the Assumption 1 for the manifold and time , the forward mapping is injective, that is, the photoacoustic imaging problem is uniquely solvable. Moreover, the following stability estimate,
holds for some constant .
We wish to make some comments before we proceed with the proof. Notice in (30) that we are only able to dominate in the norm of (rather than in the norm of its stated space ) with the measured data in the norm of . By contrast, when the Dirichlet data is measured on , then the imaging operator (left inverse of ) enjoys stability estimates as a mapping from to , or from to . See Bardos1992; Stefanov2009; Acosta-Montalto-2015 for details. Hence, there is an apparent loss of stability due to the nature of the piezoelectric measurement model (29). The double time-integration needed to invert the left-hand side of (29a) does not fully restore the regularity lost by the application of the hyperbolic differential operator on the right-hand side of (29a). This is a well-known property concerning regularity of hyperbolic equations EvansPDE. We also note that Theorem 5.1 is slightly different from what is presented in Acosta2019a where the imaging operator was shown to satisfy a stability estimate as a mapping from to . Hence, formally, there is a mild loss of stability of one degree either in space or in time, but not both.
Now, it is convenient to define the following operation
so that for any sufficiently smooth such that at . Now let
where and satisfy (29) and . Then it follows that solves the following initial value problem,
The following lemma is a well-established result. See (EvansPDE, §7.2, Thms. 3-5) or (Lio-Mag-Book-1972, Ch. 3, §8, Thm. 8.1) for details.
Let solve (33). If , then the field for . Moreover, the following stability estimate
holds for some constant .
Using this lemma, we proceed to prove the main theoretical result of the paper. In what follows, the generic constant changes from inequality to inequality, but it does not depend on , or .
Proof (Theorem 5.1)
Under Assumption 1 for the manifold and time , observability of waves from the boundary Bardos1992; GlowinskiLionsHe2008 yields that
6 Numerical Simulations
Now we propose and numerically implement a reconstruction algorithm to solve the PAT problem at the discrete level. The reconstructions presented here are based on the Landweber iterative method (EnglBook2000, Ch. 6) accelerated using Nesterov’s approach. The iterative process is defined in Algorithm 1. In this algorithms, the parameter is known as the relaxation factor which guarantees the stability and convergence of the iterations. The parameter , known as the momentum factor, allows for the acceleration of the convergence. In the presence of noise, the maximum number of iterations is chosen according to some stopping or regularization rule such as a discrepancy principle.
For the Landweber method, it is necessary to evaluate the adjoint of the forward operator . This evaluation amounts to solve the following final boundary value problem,
in order to define the adjoint mapping as
Both the forward map and its adjoint are discretized using a piecewise linear finite element method (FEM) in space and the explicit Newmark method for the time stepping. The discretization parameters are chosen to satisfy the CFL stability condition. The FEM is implemented on triangulations of the domain . For the numerical simulations, we have non-dimensionalized the physical parameters displayed in Table 1 in order to have , and . The final time . The non-dimensional parameters of the piezoelectric film have been chosen as follows , . The parameters of the backing layer are and . The piezoelectric-elastic coupling coefficient . These non-dimensional parameters are consistent with the ranges of their dimensional counterparts listed in Table 1.
Figure 3 displays a coarse mesh used for the FEM, the exact pressure profile to be reconstructed and the boundary measurements. These measurements were synthetically generated by applying the discrete version of the forward operator
using the aforementioned numerical method for the wave equation. The FEM for the reconstruction procedure has 109,762 degrees of freedom and 6,400 time steps covered the time window for
. The mesh employed to generate the measurements was more refined, with mesh size approximately half of the mesh size employed in the reconstruction steps, and the data was down-sampled to the reconstruction mesh using linear interpolation.
: White noise (constant spectral power density).B: Pink noise (spectral power decays as ). C: Red noise (spectral power decays as ). Relative error versus iteration number for various levels of noise. D: White noise. E: Pink noise. F: Red noise.
The performance of the Algorithm 1 for various values of the momentum factor is displayed in Figure 4. Significant improvements are observed for increasing values of . For instance, in order to reach below relative error, the original Landweber method () takes 36 iterations, whereas the accelerated method with takes 12 iterations. For values , some instability is observed. The error profiles for the initial guess (back-projection) and for the 12 iterate are shown in Figure 5. We observe that after a few iterations, the algorithm is able to recover most of the sharp features of the image.
The effect of noise added to the piezoelectric measurements is also explored here. We have considered three types of noise: white, pink and red. These are characterized by increasing autocorrelation distances or equivalently faster decay of their spectral power. White or Gaussian noise has equal power spectral density (PSD) at different frequencies. Pink noise has a PSD decaying as as . Red or Brownian noise has a PSD decaying as . Figure 6 displays samples from these three types of noise. Figure 6 also displays the relative error as a function of the iteration number and increasing levels of noise for the three types of noise. We observe that the implemented reconstruction method can handle best the white noise. The error for the pink noise case is greater. And the reconstruction error for the red noise case is the greatest of the three types of noise.
We have proposed a mathematical model for the acoustic measurements transduced by piezoelectric sensors. This model, taking the form (13), is highly idealized as described in the Introduction. This idealization allows us to attain a balance between correctness and simplicity in order to analyze the properties of the model (see section 3) and the solvability of the associated photoacoustic problem (see sections 4-5). The validity of the model is limited to small values for the thickness of the piezoelectric film with respect to the wavelength of the acoustic waves. Advanced industrial processes are able to manufacture piezoelectric films with thickness in the range 30 – 100 m approximately. As shown in Table 1, the wave speed in PVDF materials ranges from 1300 – 2300 m/s. Hence, we expect our model to be valid for frequencies 13 MHz.
The directional response for plane waves derived in section 3 and displayed in Figure 2 matches the experimental measurements carried out by other researchers Wilkens2007; Wear2018b; Wear2019a; Brown2000. Specifically, design considerations concerning the mechanical and electrical properties of the piezoelectric film and backing layer play a role in the appearance of critical angles where the sensitivity of the sensor vanishes. See Figure 2. The incorporation of sensor response into reconstruction algorithms has been highlighted as one of the challenges associated with photoacoustic imaging Wang2007; Cox2009a; Ellwood2014; Wang2011a; Xia2014 and partially investigated in Finch2005p; Acosta-Montalto-2015; Zangerl2019; Dreier2019; Acosta2019a. Our present paper represents a novel contribution, tailored to piezoelectric sensors, to this research effort.
From the theoretical perspective as stated in Theorem 5.1, we conclude that the photoacoustic imaging problem is solvable for piezoelectric measurements. However, some stability is lost compared to measuring directly the Dirichlet data. Based on these theoretical properties of the inverse problem, we have also proposed a reconstruction iterative Algorithm 1 based on an accelerated Landweber method. We implemented some proof-of-concept synthetic simulations with and without noise. Robustness in the presence of noise of different spectral characteristics is observed.
One limitation of the proposed method relates to the numerical characteristics of the FEM and Newmark time-stepping. Notice in Figure 4 that the error initially decays exponentially as the theory for the Landweber method predicts. However, as the iterations continue, the error eventually stagnates. This stagnation may be attributed to the numerical error introduced by the scheme employed to approximate the action of the forward operator and its adjoint . Among other issues, the discrete version of may not be an exact adjoint for the discrete version of . Also, the numerical scheme to solve the wave equation suffers from numerical dispersion. Different frequency components of the initial pressure profile travel at different group velocity towards the detection boundary. As a consequence, the discrete version of loses coercivity (becomes ill-conditioned) and the reconstructed images suffer from aberration. These complications are not mitigated by refinement of the FEM mesh. Remedies for this phenomenon have been investigated, including regularization, two-grid methods and numerical schemes for the wave equation with lower dispersion. See (GlowinskiLionsHe2008, Sect. 6.8-6.10), Zuazua2005 and references therein. However, these improvements fall outside of the scope of this paper.
Other possible areas for future research include efforts to relax the idealizations listed in the Introduction. From the mathematical point of view, it seems plausible to model the resolution limits due to the finite size of detectors, the effect of partial measurements, and the contribution from shear waves supported by the piezoelectric film and backing layer.
Acknowledgements.The author would like to thank Texas Children’s Hospital for its support and for the research oriented environment provided by the Predictive Analytics Laboratory.
Conflict of interest
The author declares no conflict of interest.