Solvability for Photoacoustic Imaging with Idealized Piezoelectric Sensors
Solvability for Photoacoustic Imaging with Idealized Piezoelectric Sensors
Solvability for Photoacoustic Imaging with Idealized Piezoelectric Sensors Sebastian ´ Acosta Abstract—Most reconstruction algorithms for photoacoustic affected by the presence of the sensors and the acoustic sensors imaging assume that the pressure ﬁeld is measured by ultrasound do not measure pressure directly. They measure a surrogate sensors placed on a detection surface. However, such sensors do for pressure that depends on the actual transducer mechanism. not measure pressure exactly due to their non-uniform directional The most common mechanisms for ultrasound applications are and frequency responses, and resolution limitations. This is based on the piezoelectric effect –, on Fabry–Perot the case for piezoelectric sensors that are commonly employed for photoacoustic imaging. In this paper, using the method interferometry – or on ﬁber-optic refractometry – of matched asymptotic expansions and the basic constitutive . In  we formulated a model speciﬁcally tailored to the relations for piezoelectricity, we propose a simple mathematical Fabry–Perot sensor design. In this paper, we derive a similar model for piezoelectric transducers. The approach simultane- model for piezoelectric sensors and analyze the well-posedness ously models how the pressure waves induce the piezoelectric of photoacoustic imaging with such measurements. In order measurements and how the presence of the sensors affects the pressure waves. Using this model, we analyze whether the data to attain a balance between accuracy and simplicity, the gathered by piezoelectric sensors leads to the mathematical model we develop here is based on the following underlying solvability of the photoacoustic imaging problem. We conclude idealizations: that this imaging problem is well-posed in certain normed spaces and under a geometric assumption. We also propose an (a) The sensors are treated as point-like detectors. Hence, iterative reconstruction algorithm that incorporates the model we do not account for resolution limitations due to the for piezoelectric measurements. A numerical implementation of ﬁnite size of sensing elements. See , , – the reconstruction algorithm is presented. for investigations concerning this issue. We also assume Index Terms—Inverse problems, thermoacoustics, optoacous- that the domain to be imaged is fully enclosed by the tics, tomography, ultrasound detection surface. (b) We assume that in each detector, the piezoelectric ﬁlm I. I NTRODUCTION is ﬂat and its thickness is small in comparison to the wavelengths under consideration. In practice, industrial HOTOACOUSTIC tomography is a non-ionizing imaging processes can manufacture piezoelectric ﬁlms with thick- modality designed to advantageously combine the high ness 30 – 100 m approximately , –. contrast of optical absorption with the high resolution from (c) Although the sensors contain elastic materials that may broadband ultrasound waves. The imaging of optical absorp- support shear waves, our analysis is valid for compres- tion reveals important functional and pathological information sional waves governed by the scalar wave equation. about biological tissues –. (d) For the piezoelectric ﬁlm, the poling direction is along One of the open challenges concerning photoacoustic inver- its thickness, and the piezoelectric properties are trans- sion is the incorporation of realistic models for acoustic mea- versely isotropic in the plane perpendicular to the poling surements. This need for modeling the physics of ultrasound direction. The sensing ﬁlm is mechanically isotropic. sensors has been recognized in –. It has been claimed (e) Sensors may have complex structures, including a casing that ultrasound measurements can be described as a linear for structural integrity, electrodes, bonding layers and combination of the pressure ﬁeld and its normal derivative at multiple paddings designed to match the mechanical the boundary. With that motivation, Dreier and Haltmeier  impedance of the acoustic medium –. However, recently established explicit formulas for the inversion of the we assume a simple design consisting of the piezoelectric two-dimensional wave equation from Neumann boundary data ﬁlm, sandwiched by electrode foils of negligible thick- for circular and elliptical domains. In a related effort, Zangerl, ness, mounted on a much thicker backing layer. This Moon and Haltmeier  derived Fourier-based reconstruction follows models described in , , . formulas for the spherical detection geometry from knowledge of Robin boundary data. An illustration of the idealized setup is shown in Figure 1. Most other reconstruction algorithm assume that waves The acoustic domain, denoted by , contains soft tissue with propagate freely across the detection boundary and that the variable density and variable wave speed c. The piezoelectric pressure ﬁeld (Dirichlet data) is measured exactly. These ﬁlm has a small uniform thickness > 0, constant density assumptions are not satisﬁed in practice. The pressure wave are and constant wave speed c . The thick backing layer p p b has constant density and constant wave speed c . The b b S. Acosta, Department of Pediatrics, Baylor College of Medicine and interface between the acoustic domain and the piezoelectric Predictive Analytics Laboratory, Texas Children’s Hospital, Houston TX, ﬁlm is denoted . The interface between the piezoelectric ﬁlm USA. (firstname.lastname@example.org) and the backing layer is denoted . arXiv:2002.09929v2 [math.AP] 30 Jun 2020 2 particle displacement u to the pressure p in the piezoelectric ﬁlm, @ u = rp (3) p p where the pressure p is deﬁned as p = ( + 2) div u: (4) Combining (3) and (4), we ﬁnd that the pressure ﬁeld p satisﬁes the wave equation, 2 2 @ p = c p (5) p p t p where the wave speed c is deﬁned by c = ( + 2) = . p p The piezoelectric transducer measures the electrical voltage V across the piezoelectric ﬁlm generated by the mechanical deformation due to the transmitted acoustic waves. We proceed Fig. 1: Acoustic domain with density and wave speed c. to derive the mathematical relationship between the voltage V Piezoelectric material of uniform thickness , density p p and the pressure p in the piezoelectric material. Our guiding and wave speed c . The thick backing layer has density p b b references are [31, Ch. 5], [32, Ch. 5] and . Under small and wave speed c . The interface between the acoustic domain b perturbations of ﬁeld conditions, the linearized constitutive and the piezoelectric ﬁlm is denoted by . The interface relation for the piezoelectric effect is the following between the piezoelectric ﬁlm and the backing layer is denoted D = " E + d (6) by . where D is the 3 1 electric displacement vector (electric charge per area), E is an externally applied 3 1 electric ﬁeld In section II we derive a model for the transduction from (voltage per length) and " is the 3 3 dielectric permittivity pressure to electrical voltage which is the physical quantity tensor (capacitance per length). Following , it is convenient acquired by the piezoelectric sensors. In section III we derive to express the symmetric stress tensor (force per area) as an effective boundary condition for the transmission of waves a 6 1 vector and the piezoelectric tensor d (electric charge from the acoustic medium of interest into the piezoelectric per force) as a 3 6 matrix, sensor. This effective transmission condition accounts for the 2 3 inﬂuence that the sensor exerts on the acoustic waves. Using 6 7 2 3 the coupled models for piezoelectric measurements and wave 6 7 0 0 0 0 d 0 6 7 propagation, in section IV we deﬁne the forward problem 6 7 4 5 d = 0 0 0 d 0 0 ; = : 6 7 associated with photoacoustic imaging. Then in section V 6 7 d d d 0 0 0 31 32 33 4 5 we state and prove the solvability of the imaging problem with piezoelectric measurements. A reconstruction algorithm is proposed in section VI where some numerical simulations In the absence of an external electric ﬁeld in (6), the normal are presented. The conclusions follow in section VII. electric displacement D is given by D = d + d + d : (7) 3 31 11 32 22 33 33 II. M ODEL FOR P IEZOELECTRIC M EASUREMENTS As assumed above, the piezoelectric tensor is transversely We start the modeling of the piezoelectric measurements isotropic, which allows us to simplify the notation as follows from the basic constitutive relations for both piezoelectric and d = d = d and d = d . Combining the constitutive ? 31 32 33 mechanical variables. Since the sensing material is mechani- relations (1) and (7) we obtain the electric displacement in cally isotropic, the 3 3 symmetric stress tensor is related terms of the strain, to the 3 3 symmetric strain tensor s as follows, D = e s + e s + es (8) 3 ? 11 ? 22 33 = (s + s + s ) + 2 s (1) ij ij 11 22 33 ij where e = 2d ( + ) + d and e = d ( + 2) + 2d . ? ? ? where the direction along the thickness of the piezoelectric As a consequence, using the deﬁnition of strain s in terms of ﬁlm is denoted as the 3-axis, and the 1-axis and 2-axis are the the displacement u, we obtain transverse plane. Here is the Kronecker delta, and and ij are the ﬁrst and second Lame ´ coefﬁcients. The equation of D = e div u + (e e ) @ (n u) (9) 3 ? ? n mechanical motion is where n is the normal vector on and @ represents the @ u = r (2) derivative along the normal direction or 3-axis. Now we take two time-derivatives of (9) and combine with (3)-(4) to obtain where u is the 3 1 material displacement vector. For irrotational deformations, i.e. in the absence of shear stress, e (e e ) ? ? 2 2 2 @ D = @ p + @ p : (10) 3 p p t t n the above equation can be simpliﬁed in order to relate the + 2 p 3 2 TABLE I: Estimates for the physical parameters of PVDF Express = @ + where represents the surface ? ? 2 piezoelectric sensors , , , , –. Laplacian on the transverse plane, recall that c = (+2)= and solve for @ p in (5) to plug it into (10) to obtain Parameter Value Units e (e e ) ? ? 2 2 2 2 2 @ D = c @ p + c @ p p PVDF thickness 10 – 60 m 3 p p ? p t p t p t p p 3 PVDF density 1780 – 1950 kg m ec (e e ) p ? 2 2 PVDF wave speed c 1300 – 2300 m s = @ p c p : (11) p ? p t p PVDF Poisson’s ratio 0.2 – 0.4 Piezoelectric coeff. d -(30 – 35) pC/N By deﬁnition of voltage as an electric potential, the voltage Piezoelectric coeff. d 3 – 15 pC/N V generated across the piezoelectric ﬁlm by the 3-component Coefﬁcient 0.3 – 1.5 of the electric displacement D satisﬁes Backing density 1900 – 2500 kg m Backing wave speed c 1000 – 4000 m s @ V = (12) where " is the dielectric permittivity along the 3-axis. Integrating (12) across the piezoelectric ﬁlm and combining takes the form of an effective impedance boundary condition the result with (11), we obtain our model for the piezoelectric that replaces the more involved transmission process for waves measurements traveling from the acoustic domain , through the piezoelec- 2 2 2 tric ﬁlm and into the backing layer . We assume that p b @ V / @ p c p (13) p ? p t t p the pressure ﬁeld p in the backing layer is outgoing which with vanishing initial state. Here the symbol / denotes equal- translates into satisfying a radiation condition of the form, ity up to a multiplicative constant (which is typically estimated @ p + c @ p +Hp = 0 on (15) n b t b b through experimental calibration). In (13), it is assumed that b the pressure ﬁeld is constant across the piezoelectric ﬁlm. This where H is the mean curvature of the surface . See ,  assumption is rigorously justiﬁed in the next section. for a derivation. The symbol appearing in the model (13) is a unitless As in , we make some geometric assumptions about the coefﬁcient deﬁned by the elastic and piezoelectric properties domain occupied by the piezoelectric ﬁlm. We let p p of the sensing ﬁlm fy 2 : 0 < dist(y; ) < g. For sufﬁciently small , the e e 2 (d d ) domain can be expressed as a family of parallel surfaces ? ? p = = parametrized by 0 < z < and deﬁned by = fy = x + e d ( + 2) + 2d z zn(x) : x 2 g where n(x) is the normal vector at x 2 . For (1 2) (1 d =d) = (14) smooth and sufﬁciently small , the surfaces are well- 1 (1 2d =d) deﬁned, smooth and mutually disjoint. Each point y 2 where we have expressed = c =(1 ) and = p can be uniquely represented in the form y = x + zn(x) for c (1 2)=(2 2) in terms of Poisson’s ratio to p x 2 and 0 < z < . In addition, the normal vector at obtain the last equality. Common values for all these physical y 2 coincides with the normal vector at x 2 . See details parameters are shown in Table I for polyvinylidene ﬂuoride concerning parallel surfaces in [38, Sect. 6.2]. (PVDF) piezoelectric sensors. The transmission of the pressure ﬁeld from the acoustic We note from (13) that a theoretically perfect transduction domain into the piezoelectric ﬁlm is governed by the following from pressure to voltage would be attained if the coefﬁcient transmission conditions at the interface , = 0. However, due to the nature of the poling processes 1 1 p = p and @ p = @ p on , (16) p n n p employed to manufacture these piezoelectric materials, the p coefﬁcients d and d have opposite signs and generally where p and p are the pressure in the acoustic medium jdj> jd j. This implies that 1 < (1 d =d) < 2. Hence, in ? ? and piezoelectric ﬁlm, respectively. The ﬁrst condition in order for = 0, the Poisson’s ratio would have to be = 0:5 (16) ensures the continuity of the pressure ﬁeld. The second which requires the piezoelectric material to be incompressible. condition in (16) ensures the continuity of particle motion in In practice, Poisson’s ratio for PVDF ﬁlms ranges from 0.2 to the normal direction. Similar transmission conditions hold at 0.4 approximately. We note that ranges from 0.3 to 1.5, for the interface , the realistic range of values for the Poisson’s ratio and the 1 1 piezoelectric ratio d =d displayed in Table I. p = p and @ p = @ p on , (17) p b n p n b ? p where p is the pressure in the backing layer. The pressures III. EFFECTIVE MODEL FOR WAVE PROPAGATION p, p and p satisfy the wave equation with respective wave p b speeds c, c and c . Typically, the piezoelectric ﬁlm and the backing layer are p b acoustically more rigid and heavier than the biological medium Now we proceed to use the method of matched asymptotic of interest. Therefore, the presence of the sensors induces expansions to derive an effective model for the interplay partial reﬂections of the waves. Here we seek to model how between the pressure ﬁelds and the piezoelectric sensor. For the sensors exert inﬂuence on the pressure waves. This model an introduction to this method, see . First we consider the 4 formal asymptotic expansions for the pressure ﬁelds, coefﬁcient satisﬁes cos c 0 1 2 p(t; x) = p (t; x) + p (t; x) +O( ) (18a) R = ; where = : (25) cos + c 0 1 2 b b p (t; x; z) = p (t; x; z) + p (t; x; z) +O( ) (18b) p p After plugging (24)-(25) into the model (13) and evaluating at 0 1 2 p (t; x) = p (t; x) + p (t; x) +O( ) (18c) b b the origin x = 0, we ﬁnd that the piezoelectric measurements satisfy the following directivity pattern and introduce a change of variable in order to extract the effect of the piezoelectric ﬁlm thickness , V cos = 1 + 1 sin (26) z = for 2 [0; 1]: (19) p cos + c inc The boundary value problem for the pressure ﬁeld p in the p where is given by (14). Figure 2 displays the directional piezoelectric ﬁlm, governed by the wave equation (5) and the response (26) in decibels as a function of the incidence transmission conditions (16)-(17), can be recast in terms of angle and piezoelectric coefﬁcient over a realistic range and terms with same powers of are gathered to obtain the of values shown in Table I. We observe that for values of 2 2 following cases. > c =c , a critical angle appears. For incidence at this a) O( )-terms: critical angle, vanishing measurements are obtained by the piezoelectric sensor design. This critical angle is given by 2 0 @ p = 0; for 2 (0; 1), 1 1=2 = arcsin c c . cr 0 0 0 p = p and @ p = 0; at = 0, p p 0 0 0 p = p and @ p = 0; at = 1, p b p which imply that p is constant as a function of and that the ﬁrst effective transmission condition is that 0 0 0 p = p = p ; on : (20) p b b) O( )-terms: 2 1 @ p = 0; for 2 (0; 1), 1 1 1 0 1 1 p = p and @ p = @ p ; at = 0, p p p 1 1 1 1 1 0 p = p and @ p = @ p ; at = 1, p b p p b b which imply that @ p is constant as a function of and that the second effective transmission condition is 1 0 1 1 1 0 @ p = @ p = @ p on : (21) n n p p b Combining (15) and (20)-(21), we obtain closed-form ef- fective governing equations for the leading order term p of Fig. 2: Directivity (26) in decibels for piezoelectric sensor as the acoustic ﬁeld in the domain a function of the incidence angle and coefﬁcient . The parameters correspond to a PVDF ﬁlm with (compressional) 2 0 2 0 @ p c p = 0 in ft > 0g ; (22) wave speed c = 2000 m/s and density = 1800 kg/m p p 0 1 0 0 @ p + c @ p + Hp = 0 on ft > 0g : (23) b n t and a backing layer with (compressional) wave speed c = 1000 m/s and density = 2000 kg/m . The acoustic medium Similar models for photoacoustics are studied in , , , corresponds to water with wave speed c = 1500 m/s and . density = 1000 kg/m . As an example for the response of the piezoelectric sensor design, we can analyze its behavior for plane waves and for a ﬂat boundary . Both the boundary value problem (22)-(23) and the model for the measurements (13) play an IV. T HE FORWARD PROBLEM important role in this analysis. A plane wave of the form i(xk !t) p = e propagating in the direction of k, induces Now we proceed to deﬁne the forward problem of pho- inc a reﬂection governed by (23). The total pressure ﬁeld p is the toacoustic imaging in terms of the wave propagation model superposition of the incident and reﬂected wave, (22)-(23) and the model for piezoelectric measurements (13). We neglect higher order terms O() studied in the previous i(xk !t) i(xk !t) p(x; t) = e + Re +O() (24) section, so that the pressure ﬁeld p is assumed to satisfy the following initial value problem, where R is the reﬂection coefﬁcient, k is the reﬂection wavenumber satisfying jkj= jk j= !=c and n k = n k, r r 2 2 @ p c p = 0 in (0; T ) (27a) where n is the outward normal on . We can write n k = @ p + c @ p + Hp = 0 on (0; T ) (27b) b n t jkjcos where is the angle of incidence. Plugging (24) into b (23) and neglecting the O() terms, we ﬁnd that the reﬂection p = f and @ p = 0 on ft = 0g (27c) t 5 where 0 < T < 1 is the measurement time to be determined f in the norm of H ( ) (rather than in the norm of its later. Recall that the underlying assumption concerning media stated space H ( )) with the measured data V in the norm 1 0 properties are that c is bounded from below and above, and is of H ([0; T ]; H ( )). By contrast, when the Dirichlet data smooth in , and that c , c , , , and are constants. is measured on [0; T ] , then the imaging operator (left p b p b The forward mapping is given by inverse of F ) enjoys stability estimates as a mapping from 0 0 1 1 H ([0; T ] ) to H ( ), or from H ([0; T ] ) to H ( ). F : f 7! V (28) See , ,  for details. Hence, there is an apparent loss of stability due to the nature of the piezoelectric measurement where, according to the piezoelectric model (13), the measured model (29). The double time-integration needed to invert the electric voltage V satisﬁes left-hand side of (29a) does not fully restore the regularity 2 2 2 @ V = @ p c p on (0; T ) (29a) t t p lost by the application of the hyperbolic differential operator V = @ V = 0 on ft = 0g (29b) on the right-hand side of (29a). This is a well-known property concerning regularity of hyperbolic equations . We also for the pressure ﬁeld p evolving according to (27) from the note that Theorem 1 is slightly different from what is presented initial condition f . The mapping (28), which we seek to invert in  where the imaging operator was shown to satisfy a for photoacoustic imaging, encodes the physical principles 0 1 stability estimate as a mapping from H ([0; T ]; H ( )) to for acoustic wave propagation from the unknown pressure H ( ). Hence, formally, there is a mild loss of stability of proﬁle f to the measured electrical voltage V generated by one degree either in space or in time, but not both. the piezoelectric sensors. Now, it is convenient to deﬁne the following operation We work with the standard Sobolev spaces based on square- integrable functions deﬁned on the domain or the boundary (@ v)(t) = v(s) ds (31) (0; T ) . The associated inner product extends as the duality pairing between functionals and functions. For the Sobolev 1 1 so that @ @ v = @ @ v = v for any sufﬁciently smooth v 0 2 t t t t space H ( ), the inner product is weighted by c so that such that v = 0 at t = 0. Now let the differential operator c is formally self-adjoint. The well- posedness in Sobolev spaces of the initial value problem (27) u = @ p (32) is a well-established result , . where p and V satisfy (29) and p(0) = f 2 H ( ). Then it follows that u solves the following initial value problem, V. THE I NVERSE P ROBLEM 2 2 @ u c u = @ V on (0; T ) (33a) The inverse problem associated with photoacoustic imaging t p ? t is the following: Given the voltage measurements V modeled u = @ u = 0 on ft = 0g : (33b) by (29) on (0; T ), induced by the pressure ﬁeld p satisfying The following lemma is a well-established result. See [42, (27), ﬁnd the unknown initial condition f . The solvability of §7.2, Thms. 3-5] or [43, Ch. 3, §8, Thm. 8.1] for details. this inverse problem depends on the geometry of the domain 1 0 Lemma 1: Let u solve (33). If V 2 H ([0; T ]; H ( )), then , the proﬁle of the wave speed c and the time T < 1. These k 1 k the ﬁeld u 2 C ([0; T ]; H ( )) for k = 0; 1. Moreover, the conditions are made precise in the following assumption, following stability estimate known as the geometric control condition or a nontrapping condition for the geodesic ﬂow. We work with the manifold kuk k 1 k CkVk 1 0 (34) C ([0;T ];H ( )) H ([0;T ];H ( )) 2 2 endowed with the Riemannian metric c dx . See ,  holds for some constant C > 0. for details. 1 0 The deﬁnition of the Bochner spaces H ([0; T ]; H ( )) and Assumption 1 (Nontrapping condition): Let be a simply k 1 k C ([0; T ]; H ( )) can be found in , . Using Lemma connected bounded domain with smooth boundary . Assume 1, we proceed to prove the main theoretical result of the paper. there exists T < 1 such that any (unit speed) geodesic ray 2 2 In what follows, the generic constant C > 0 changes from of the manifold ( ; c dx ), originating from any point in inequality to inequality, but it does not depend on f , p or V . at time t = 0, reaches the boundary at a nondiffractive point Proof of Theorem 1: Under Assumption 1 for the before t = T . 2 2 manifold ( ; c dx ) and time T > T , observability of With this assumption in place, we can state the main result o waves from the boundary ,  yields that of the paper in the form of a theorem. Theorem 1: Under the Assumption 1 for the manifold 0 0 kfk Ckpk (35) H ( ) H ([0;T ] ) 2 2 ; c dx ) and time T > T , the forward mapping F : 1 1 0 H ( ) ! H ([0; T ]; H ( )) is injective, that is, the pho- for some constant C = C ( ; c; T ). Now, from the deﬁnition toacoustic imaging problem is uniquely solvable. Moreover, (32) of u and the stability estimate in Lemma 1 for k = 1, we the following stability estimate, obtain that kfk 0 CkVk 1 0 (30) kpk 1 0 CkVk 1 0 : (36) H ( ) H ([0;T ];H ( )) C ([0;T ];H ( )) H ([0;T ];H ( )) 1 0 holds for some constant C > 0. Since the norm of C ([0; T ]; H ( )) dominates the norm of We wish to make some comments before we proceed with H ([0; T ] ), combining (35) and (36) we obtain the desired the proof. Notice in (30) that we are only able to dominate result (30). 6 VI. NUMERICAL SIMULATIONS Now we propose and numerically implement a reconstruc- tion algorithm to solve the PAT problem at the discrete level. The reconstructions presented here are based on the Landwe- ber iterative method [47, Ch. 6] with Nesterov’s acceleration. See Stefanov and Yang  for an excellent analysis of the Landweber method for a similar problem. The iterative process is deﬁned in Algorithm 1. In brief, the Landweber method is based on inverting (28) by solving (I (I
F F )) f =
F V where the normal operator (F F ) is positive deﬁnite and > 0 is chosen small enough so that the spectrum of (I
F F ) is contained in ( 1; 1) making it a contraction. The parameter
is known as the relaxation factor. The parameter , known as the momentum factor, allows for the acceleration of the convergence. Unfortunately, it is hard to choose
and optimally. We resort to trial and error to set them satisfactorily. In the presence of noise, the number K of iterations is chosen according to a regularization rule. For the Landweber method, it is necessary to evaluate the adjoint F of the forward operator F . This evaluation amounts to solve the following ﬁnal boundary value problem, 2 2 @ '