Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

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 field 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 [11]–[13], on Fabry–Perot the case for piezoelectric sensors that are commonly employed for photoacoustic imaging. In this paper, using the method interferometry [14]–[17] or on fiber-optic refractometry [18]– of matched asymptotic expansions and the basic constitutive [20]. In [21] we formulated a model specifically 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 finite size of sensing elements. See [12], [13], [22]–[25] 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 film I. I NTRODUCTION is flat 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 films with thick- modality designed to advantageously combine the high ness 30 – 100 m approximately [11], [26]–[29]. 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 [1]–[4]. (d) For the piezoelectric film, 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 film is mechanically isotropic. sensors has been recognized in [5]–[8]. 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 field and its normal derivative at multiple paddings designed to match the mechanical the boundary. With that motivation, Dreier and Haltmeier [9] impedance of the acoustic medium [30]–[32]. 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 film, 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 [10] derived Fourier-based reconstruction follows models described in [26], [27], [32]. 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 field (Dirichlet data) is measured exactly. These film has a small uniform thickness  > 0, constant density assumptions are not satisfied 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, film is denoted . The interface between the piezoelectric film USA. (sebastian.acosta@bcm.edu) 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 film, @ u = rp (3) p p where the pressure p is defined as p = ( + 2) div u: (4) Combining (3) and (4), we find that the pressure field p satisfies the wave equation, 2 2 @ p = c p (5) p p t p where the wave speed c is defined by c = ( + 2) = . p p The piezoelectric transducer measures the electrical voltage V across the piezoelectric film 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 [27]. Under small and wave speed c . The interface between the acoustic domain b perturbations of field conditions, the linearized constitutive and the piezoelectric film is denoted by . The interface relation for the piezoelectric effect is the following between the piezoelectric film 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 field 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 [27], 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 influence 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 define 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 field 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 definition of strain s in terms of film 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 first and second Lame ´ coefficients. 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 simplified 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 [11], [26], [27], [30], [32]–[35]. 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 definition of voltage as an electric potential, the voltage Piezoelectric coeff. d 3 – 15 pC/N V generated across the piezoelectric film by the 3-component Coefficient  0.3 – 1.5 of the electric displacement D satisfies 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 film 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 film and into the backing layer . We assume that p b @ V / @ p c  p (13) p ? p t t p the pressure field 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 field is constant across the piezoelectric film. This where H is the mean curvature of the surface . See [36], [37] assumption is rigorously justified in the next section. for a derivation. The symbol  appearing in the model (13) is a unitless As in [21], we make some geometric assumptions about the coefficient defined by the elastic and piezoelectric properties domain occupied by the piezoelectric film. We let p p of the sensing film fy 2 : 0 < dist(y; ) < g. For sufficiently small , the e e 2 (d d ) domain can be expressed as a family of parallel surfaces ? ? p = = parametrized by 0 < z <  and defined 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 sufficiently small , the surfaces are well- 1  (1 2d =d) defined, 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 fluoride concerning parallel surfaces in [38, Sect. 6.2]. (PVDF) piezoelectric sensors. The transmission of the pressure field from the acoustic We note from (13) that a theoretically perfect transduction domain into the piezoelectric film is governed by the following from pressure to voltage would be attained if the coefficient 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 coefficients 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 film, respectively. The first condition in order for  = 0, the Poisson’s ratio would have to be  = 0:5 (16) ensures the continuity of the pressure field. 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 films 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 film 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 reflections of the waves. Here we seek to model how between the pressure fields and the piezoelectric sensor. For the sensors exert influence on the pressure waves. This model an introduction to this method, see [39]. First we consider the 4 formal asymptotic expansions for the pressure fields, coefficient satisfies 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 find that the piezoelectric measurements satisfy the following directivity pattern and introduce a change of variable in order to extract the effect of the piezoelectric film thickness , V cos = 1 + 1  sin  (26) z =   for  2 [0; 1]: (19) p cos  + c inc The boundary value problem for the pressure field p in the p where  is given by (14). Figure 2 displays the directional piezoelectric film, 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 coefficient  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 first 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 field in the domain a function of the incidence angle  and coefficient . The parameters correspond to a PVDF film 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 [8], [21], [40], corresponds to water with wave speed c = 1500 m/s and [41]. 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 flat 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 define the forward problem of pho- inc a reflection governed by (23). The total pressure field p is the toacoustic imaging in terms of the wave propagation model superposition of the incident and reflected 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 field p is assumed to satisfy the following initial value problem, where R is the reflection coefficient, k is the reflection 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 find that the reflection 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 [8], [44], [46] 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 satisfies 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 [42]. We also for the pressure field 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 [21] 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 profile 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 define the following operation We work with the standard Sobolev spaces based on square- integrable functions defined 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 sufficiently 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 [42], [43]. 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 field p satisfying The following lemma is a well-established result. See [42, (27), find 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 profile of the wave speed c and the time T < 1. These k 1k the field 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 flow. We work with the manifold kuk k 1k  CkVk 1 0 (34) C ([0;T ];H ()) H ([0;T ];H ()) 2 2 endowed with the Riemannian metric c dx . See [44], [45] holds for some constant C > 0. for details. 1 0 The definition of the Bochner spaces H ([0; T ]; H ()) and Assumption 1 (Nontrapping condition): Let be a simply k 1k C ([0; T ]; H ()) can be found in [42], [43]. 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 [44], [45] 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 definition 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 [48] for an excellent analysis of the Landweber method for a similar problem. The iterative process is defined 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 definite 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 final boundary value problem, 2 2 @ ' c ' = 0 (0; T ) (37a) @ ' c @ ' + H' =   (0; T ) (37b) b n t b ' = 0 and @ ' = 0 ft = Tg (37c) where  solves 2 2 2 @  = @ c  in (0; T ) (38a) t t p = 0 and @  = 0 on ft = Tg (38b) in order to define the adjoint mapping as F : 7! @ 'j : (39) t t=0 It can be shown, using straight-forward integration by parts, that indeed F is the adjoint of F or equivalently that hf; @ 'j i = hp; i t t=0 [0;T ] 2 2 = hp; I c  @ i ? [0;T ] p t = hV; i [0;T ] for all sufficiently smooth f and , where V = F (f ) according to (28) and @ 'j = F ( ) according to (39). t t=0 The brackets h;i denote the inner product of the Hilbert space H (X ) which extends as the duality pairing between the Sobolev functionals and functions. Algorithm 1 Accelerated Landweber iteration Set 0 < K , 0 < < 2kFk and 0   < 1. Initial guesses u = F V and v = u 0 0 0 for k = 1; 2; :::; K do v = u (F Fu u ) =ku k k k1 k1 0 0 Fig. 3: A: Coarse mesh for the FEM method. B: Exact profile u = v +  (v v ) k k k k1 to be reconstructed showing the brain vasculature imaged with MRI technology [49]. C: Synthetic piezoelectric measurements return u as modeled by (28). 7 Both the forward map F and its adjoint F 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 I in order to have c = 1, diam( ) = 2 and  = 1. The final time T = 2. The non-dimensional parameters of the piezoelectric film have been chosen as follows  = 1:5, c = 1:0. The p p parameters of the backing layer are  = 2:0 and c = 1:0. b b The piezoelectric-elastic coupling coefficient  = 0:9. These non-dimensional parameters are consistent with the ranges of their dimensional counterparts listed in Table I. Figure 3 displays a coarse mesh used for the FEM, the exact pressure profile to be reconstructed and the boundary measure- ments. These measurements were synthetically generated by applying the discrete version of the forward operator F 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 T = 2. 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. 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 1% relative error, the original Landweber method ( = 0:0) takes 36 iterations, whereas the accelerated method with  = 0:6 takes 12 iterations. For values  > 0:7, some instability is observed. =0.0 =0.2 10% =0.4 =0.6 1% Fig. 5: Error profiles for synthetic measurements obtained using the piezoelectric model (28). A: Reconstruction using an interpretation of the measurements as piezoelectric data. The relative error is 0:54%. B: Reconstruction using a naive interpretation of the measurements as Dirichlet data. The relative error is 6:43%. 0.1% 0 10 20 30 40 50 Iteration the same synthetic measurements are used. For the first Fig. 4: Relative error versus iteration number, for various reconstruction, we properly interpret the given measurements values of the momentum factor  to accelerate the Landweber as generated by the piezoelectric model (28) and carry out iterations. The relaxation factor = 5 10 in all cases. the reconstruction using Algorithm 1. After 50 iterations we In order to visualize the impact of improperly modeling obtain a relative error of 0:54%. For the second reconstruction, the piezoelectric measurements, we have implemented two we improperly interpret the measurements as the Dirichlet reconstructions of the initial acoustic profile. In both case, data of the pressure field. This is the naive model commonly Rel Error 8 K k XX employed by others but inconsistent with the piezoelectric = F c (FF ) V (40) kj transduction. The reconstruction is carried out using Algorithm k=0 j=0 1 modified by setting  = 0 which is equivalent to assuming for some coefficients c . This last expression motivates the kj that the measurements are Dirichlet data. After 50 iterations study of how the discrete normal operator (FF ) responds to we obtain a relative error of 6:43%. The error profiles for both the noise contained in the piezoelectric measurements. Figure reconstructions are displayed in Figure 5. 7 displays the average power spectral response to time- and space-tracings of white noise as the input to (FF ). We observe that the operator (FF ) suppresses the high frequency white noise components. This explains why the reconstruction algorithm pink noise can handle white noise better than pink or red noise. The latter red noise noise types have a larger portion of their power residing over 10% low frequencies. Therefore, the reconstruction algorithm does not suppress those types of noise. -1 1% 0 5 10 15 20 25 -2 Iteration Fig. 6: Relative error versus iteration number for white, pink -3 and red noise. In all cases the noise level is at 10%. -4 The effect of noise added to the piezoelectric measurements space is also explored. We have considered three types of noise: time white, pink and red. These are characterized by increasing -1 0 1 2 3 10 10 10 10 10 autocorrelation distances or equivalently faster decay of their Normalized Frequency spectral power. White or Gaussian noise has equal power spectral density (PSD) at different frequencies. Pink noise has Fig. 7: Frequency response of the normal operator (FF ) to a PSD decaying like ! as ! ! 1. Red or Brownian noise white noise as the input. has a PSD decaying like ! . Figure 6 displays the relative error as a function of the iteration number for the three types of noise at a 10% level. We observe that the reconstruction VII. D ISCUSSION AND CONCLUSIONS method can handle best the white noise. The error for the We have proposed a mathematical model for the acous- pink noise is greater. And the reconstruction error for the red tic measurements transduced by piezoelectric sensors. This noise is the greatest of the three types of noise. model, taking the form (28), is highly idealized as described in This phenomenon could be explained by studying the spec- the Introduction. This idealization allows us to attain a balance tral characteristics of the discrete normal operators (F F ) between correctness and simplicity in order to analyze the or (FF ). The iterations from the Landweber algorithm properties of the model (see section III) and the solvability of correspond to truncated Neumann series [48], the associated photoacoustic problem (see sections IV-V). The directional response for plane waves derived in sec- tion III and displayed in Figure 2 matches the experimental u = (I F F ) (F V ); measurements carried out by other researchers [11]–[13], [30]. k=0 Design considerations for the mechanical and electrical prop- where binomial formula yields erties of the piezoelectric film (condensed in the parameter defined in (14)) play a role in the appearance of critical angles k j (I F F ) = ( F F ) : where the sensitivity of the sensor vanishes. The incorporation j=0 of this type of sensor response into reconstruction algorithms Therefore, the Landweber iterate u could be expressed as has been highlighted as one of the challenges associated with photoacoustic imaging [7], [24], [50]–[52] and partially K k XX investigated in [6], [8]–[10], [21]. Our present paper represents u = c (F F ) (F V ) K kj a novel contribution to this research effort. k=0 j=0 Rel Error Spectral Power 9 We note that our mathematical model for the piezoelectric an impedance Z such that Z < Z < Z . Z so that the ml ml p b sensor is qualitatively similar to the model for the Fabry– wave field experiences a more gradual change of media across Perot transducer proposed in [21] in spite of the completely the fluid-sensor interface. different physical principles from which they are derived. The validity of the proposed asymptotic model is limited to Hence, a common mathematical framework can be used to small values for the thickness  of the piezoelectric film with analyze both types of sensors. From the theoretical perspective, respect to the wavelength of the acoustic waves. Advanced as stated in Theorem 1, we conclude that the photoacoustic industrial processes are able to manufacture piezoelectric films with thickness in the range 30 – 100 m approximately. As imaging problem is solvable for piezoelectric measurements. shown in Table I, the wave speed in PVDF materials ranges However, some stability is lost compared to measuring directly from 1300 – 2300 m/s. Hence, we expect our model to be the Dirichlet data. valid for frequencies . 12 MHz. For higher frequencies, the We have also proposed an iterative reconstruction Algorithm wave field is affected by the presence of the piezoelectric 1 based on an accelerated Landweber method. We imple- film and the thin matching layer [30]. In such a scenario, mented proof-of-concept synthetic simulations to highlight the our transmission model would be inappropriate. Hence, there importance of incorporating the proper modeling of the piezo- remains a need for a more complete model that can accurately electric transducer. We carried out reconstructions with and handle multiple frequency scales. without the proper piezoelectric model. See the error profiles shown in Figure 5. We observe that certain features, whose ACKNOWLEDGMENT wavefronts may have reached the measurement boundary at a The author would like to thank Texas Children’s Hospital non-normal incidence, are missed by the naive reconstruction. for its support and for the research oriented environment We also investigated the effect of noise of different spectral provided by the Predictive Analytics Laboratory. The author characteristics. The reconstruction method can handle white gratefully acknowledges support by the NSF grant DMS- noise better than pink or red noise. See Figures 6. This 1712725. The author also thanks the anonymous reviewers is due to suppression of high-frequency components in the whose observations were most constructive. measurements. See Figure 7. One limitation of the proposed method relates to the numer- REFERENCES ical characteristics of the FEM and Newmark time-stepping. [1] X. Wang, Y. Pang, G. Ku, X. Xie, G. Stoica, and L. V. Wang, Notice in Figure 4 that the error initially decays exponentially “Noninvasive laser-induced photoacoustic tomography for structural and as the theory for the Landweber method predicts. However, functional in vivo imaging of the brain,” Nature Biotechnology, vol. 21, as the iterations continue, the error eventually stagnates. This no. 7, pp. 803–806, 2003. [2] P. Beard, “Biomedical photoacoustic imaging,” Interface Focus, vol. 1, stagnation may be attributed to the fact that the measurements pp. 602–31, 2011. were synthetically generated on a mesh different from the [3] K. Wang and M. Anastasio, “Photoacoustic and Thermoacoustic To- reconstruction mesh. Among other issues, the discrete version mography: Image Formation Principles,” in Handbook of Mathematical Methods in Imaging, O. Scherzer, Ed. Springer New York, 2011, pp. ofF may not be an exact adjoint for the discrete version ofF . 781–815. Also, the numerical scheme to solve the wave equation suffers [4] L. V. Wang and S. Hu, “Photoacoustic tomography: In vivo imaging from organelles to organs,” Science, vol. 335, pp. 1458–1462, 2012. from numerical dispersion. Different frequency components [5] G. Johansson and A. J. Niklasson, “Approximate dynamic boundary of the initial pressure profile travel at different group velocity conditions for a thin piezoelectric layer,” International Journal of Solids towards the detection boundary. As a consequence, the discrete and Structures, vol. 40, no. 13-14, pp. 3477–3492, 2003. [6] D. Finch, “On a thermoacoustic transform,” in Proc. 8th Int. Meeting version of (F F ) loses coercivity (becomes ill-conditioned) on Fully 3D Image Reconstruction in Radiology and Nuclear Medicine, and the reconstructed images suffer from aberration. Reme- 2005, pp. 150–151. dies for this phenomenon have been investigated, including [7] L. V. Wang and X. Yang, “Boundary conditions in photoacoustic tomography and image reconstruction,” Journal of Biomedical Optics, regularization, two-grid methods and numerical schemes with vol. 12, no. 1, p. 014027, 2007. lower dispersion. See [45, Sect. 6.8-6.10], [53] and references [8] S. Acosta and C. Montalto, “Multiwave imaging in an enclosure with therein. However, these improvements fall outside of the scope variable wave speed,” Inverse Problems, vol. 31, no. 6, p. 065009, 2015. of this paper. [9] F. Dreier and M. Haltmeier, “Explicit inversion formulas for the two- dimensional wave equation from Neumann traces,” arXiv:1905.03460, As future research, it would be very important to explicitly [10] G. Zangerl, S. Moon, and M. Haltmeier, “Photoacoustic tomography model the influence of matching layers commonly employed with direction dependent data : An exact series reconstruction approach,” in the design of piezoelectric sensors [30], [32]. Matching Inverse Problems, vol. 35, p. 114005, 2019. layers play an important role in optimizing the transmission of [11] V. Wilkens and W. Molkenstruck, “Broadband PVDF membrane hy- high-frequency acoustic energy into the sensor. In this paper, drophone for comparisons of hydrophone calibration methods up to 140 MHz,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control, vol. 54, no. 9, this transmission is modeled by (27b) which for a flat surface pp. 1784–1791, 2007. simplifies to [12] K. A. Wear, C. Baker, and P. Miloro, “Directivity and Frequency- 1 @p 1 @p Dependent Effective Sensitive Element Size of Needle Hydrophones: + = 0: Predictions From Four Theoretical Forms Compared With Measure- @n Z @t ments,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Fre- Here Z =  c is the acoustic impedance of the thick backing quency Control, vol. 65, no. 10, pp. 1781–1788, 2018. b b b [13] ——, “Directivity and Frequency-Dependent Effective Sensitive Ele- substrate which usually does not match the acoustic impedance ment Size of Membrane Hydrophones: Theory versus Experiment,” Z = c of the fluid or the acoustic impedance Z =  c of p p p IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Con- the piezoelectric film. Matching layers are designed to have trol, vol. 66, no. 11, pp. 1723–1730, 2019. 10 [14] P. Beard, F. Perennes, and T. Mills, “Transduction mechanisms of pyroelectric effect,” Acoustical Physics, vol. 64, no. 6, pp. 789–795, the Fabry-Perot polymer film sensing concept for wideband ultrasound 2018. detection,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Fre- [36] X. Antoine, H. Barucq, and A. Bendali, “Bayliss-Turkel-like radiation quency Control, vol. 46, no. 6, pp. 1575–1582, 1999. conditions on surfaces of arbitrary shape,” Journal of Mathematical [15] B. T. Cox and P. C. Beard, “The frequency-dependent directivity of Analysis and Applications, vol. 229, pp. 184–211, 1999. a planar Fabry-Perot polymer film ultrasound sensor,” IEEE Trans. [37] S. Acosta, “High order surface radiation conditions for time-harmonic Ultrason. Ferroelectr. Freq. Control, vol. 54, no. 2, pp. 394–404, 2007. waves in exterior domains,” Computer Methods in Applied Mechanics and Engineering, vol. 322, pp. 296–310, 2017. [16] E. Zhang, J. Laufer, and P. Beard, “Backward-mode multiwavelength [38] R. Kress, Linear Integral Equations, 2nd ed., ser. Applied mathematical photoacoustic scanner using a planar Fabry-Perot polymer film ultra- sound sensor for high-resolution three-dimensional imaging of biological sciences. New York : Springer-Verlag, 1999, vol. 82. tissues,” Applied Optics, vol. 47, no. 4, pp. 561–577, 2008. [39] E. Hinch, Perturbation Methods. Cambridge Univ. Press, 2012. [40] S. Acosta and C. Montalto, “Photoacoustic imaging taking into account [17] J. A. Guggenheim, E. Z. Zhang, and P. C. Beard, “A method for thermodynamic attenuation,” Inverse Problems, vol. 32, no. 11, p. measuring the directional response of ultrasound receivers in the range 115001, 2016. 0.3-80 Mhz using a laser-generated ultrasound source,” IEEE Trans. [41] S. Acosta, “Recovery of pressure and wave speed for photoacoustic Ultrason. Ferroelectr. Freq. Control, vol. 64, no. 12, pp. 1857–1863, imaging under a condition of relative uncertainty,” Inverse Problems, vol. 35, pp. 46–49, 2019. [18] G. Wild and S. Hinckley, “Acousto-ultrasonic optical fiber sensors: [42] L. C. Evans, Partial Differential Equations, ser. Graduate Studies in Overview and state-of-the-art,” IEEE Sensors Journal, vol. 8, no. 7, Mathematics. Providence, RI: American Mathematical Society, 1998, pp. 1184–1193, 2008. vol. 19. [19] K. A. Wear and S. M. Howard, “Directivity and frequency-dependent [43] J.-L. Lions and E. Magenes, Non-homogeneous boundary value prob- effective sensitive element size of a reflectance-based fiber-optic hy- lems and applications, ser. Grundlehren der mathematischen Wis- drophone: Predictions from theoretical models compared with mea- senschaften in Einzeldarstellungen, Bd. 181-183. Berlin, New York, surements,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Fre- Springer-Verlag, 1972. quency Control, vol. 65, no. 12, pp. 2343–2348, 2018. [44] C. Bardos, G. Lebeau, and J. Rauch, “Sharp sufficient conditions for [20] G. Wissmeyer, M. A. Pleitez, A. Rosenthal, and V. Ntziachristos, the observation, control, and stabilization of waves from the boundary,” “Looking at sound: optoacoustics with all-optical ultrasound detection,” SIAM Journal on Control and Optimization, vol. 30, no. 5, pp. 1024– Light: Science and Applications, vol. 7, p. 53, 2018. 1065, 1992. [21] S. Acosta, “Well-posedness for photoacoustic tomography with Fabry- [45] R. Glowinski, J.-L. Lions, and J. He, Exact and Approximate Control- Perot sensors,” SIAM Journal on Imaging Sciences, vol. 12, no. 4, pp. lability for Distributed Parameter Systems : A Numerical Approach, 1669–1685, 2019. ser. Encyclopedia of Mathematics and its Applications. Cambridge [22] M. Xu and L. V. Wang, “Analytic explanation of spatial resolution University Press, 2008, vol. 117. related to bandwidth and detector aperture size in thermoacoustic or [46] P. Stefanov and G. Uhlmann, “Thermoacoustic tomography with variable photoacoustic reconstruction,” Physical Review E, vol. 67, p. 056605, sound speed,” Inverse Problems, vol. 25, p. 075011, 2009. [47] H. Engl, M. Hanke, and A. Neubauer, Regularization of Inverse Prob- [23] M. Haltmeier and G. Zangerl, “Spatial resolution in photoacoustic lems, ser. Mathematics and Its Applications. Springer, 2000, vol. 375. tomography: Effects of detector size and detector bandwidth,” Inverse [48] P. Stefanov and Y. Yang, “Multiwave tomography with reflectors: Problems, vol. 26, no. 12, p. 125002, 2010. Landweber’s iteration,” Inverse Problems & Imaging, vol. 11, no. 2, [24] K. Wang, S. Ermilov, R. Su, H. Brecht, A. Oraevsky, and M. Anastasio, pp. 373–401, 2017. “Imaging model incorporating ultrasonic transducer properties for three- [49] W. Koroshetz, J. Gordon, A. Adams, A. Beckel-Mitchener, J. Churchill, dimensional optoacoustic tomography,” IEEE Transactions on Medical G. Farber, M. Freund, J. Gnadt, N. S. Hsu, N. Langhals, S. Lisanby, Imaging, vol. 30, no. 2, pp. 203–214, 2011. G. Liu, G. C. Peng, K. Ramos, M. Steinmetz, E. Talley, and S. White, [25] H. Roitner, M. Haltmeier, R. Nuster, D. P. O’Leary, T. Berer, G. Paltauf, “The state of the NIH brain initiative,” Journal of Neuroscience, vol. 38, H. Grun, ¨ and P. Burgholzer, “Deblurring algorithms accounting for the no. 29, pp. 6427–6438, 2018. finite detector size in photoacoustic tomography,” Journal of Biomedical [50] B. Cox, J. Laufer, and P. Beard, “The challenges for quantitative Optics, vol. 19, no. 5, p. 056011, 2014. photoacoustic imaging,” Proc. SPIE, vol. 7177, pp. 1–9, 2009. [26] B. Fay, G. Ludwig, C. Lankjaer, and P. A. Lewin, “Frequency response [51] R. Ellwood, E. Zhang, P. Beard, and B. Cox, “Photoacoustic imaging of PVDF needle-type hydrophones,” Ultrasound in Medicine and Biol- using reflectors to enhance planar arrays,” Journal of Biomedical Optics, ogy, vol. 20, no. 4, pp. 361–366, 1994. vol. 19, p. 126012, 2014. [27] J. Sirohi and I. Chopra, “Fundamental understanding of piezoelectric [52] J. Xia, J. Yao, and L. V. Wang, “Photoacoustic tomography: Principles strain sensors,” Journal of Intelligent Material Systems and Structures, and advances,” Progress In Electromagnetics Research, vol. 147, pp. vol. 11, no. 4, pp. 246–257, 2000. 1–22, 2014. [28] Y. Su, F. Zhang, K. Xu, J. Yao, and R. Wang, “A photoacoustic [53] E. Zuazua, “Propagation, observation, and control of waves approxi- tomography system for imaging of biological tissues,” Journal of Physics mated by finite difference methods,” SIAM Review, vol. 47, no. 2, pp. D: Applied Physics, vol. 38, pp. 2640–2644, 2005. 197–243, 2005. [29] M. A. Oreilly and K. Hynynen, “A PVDF receiver for ultrasound mon- itoring of transcranial focused ultrasound therapy,” IEEE Transactions on Biomedical Engineering, vol. 57, no. 9, pp. 2286–2294, 2010. [30] L. F. Brown, “Design considerations for piezoelectric polymer ultra- sound transducers,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control, vol. 47, no. 6, pp. 1377–1396, 2000. [31] J. Tichy, ´ J. Erhart, E. Kittinger, and J. P´ ıvratska, ´ Fundamentals of Sebastian Acosta received a B.Sc. degree in me- piezoelectric sensorics: Mechanical, dielectric, and thermodynamical chanical engineering and a M.Sc. degree in math- properties of piezoelectric materials. Springer-Verlag Berlin Heidel- ematics from Brigham Young University, Provo, berg, 2010. Utah, USA in 2009 and 2011 respectively, and a [32] T. L. Szabo, “Transducers,” in Diagnostic Ultrasound Imaging: Inside Ph.D. degree in computational and applied mathe- Out, 2nd ed. Academic Press, 2014, ch. 5, pp. 121–165. matics from William Marsh Rice University, Hous- [33] B. Zeqiri, P. N. Gelat, ´ J. Barrie, and C. J. Bickley, “A novel pyroelectric ton, Texas, USA in 2014. He is currently a re- method of determining ultrasonic transducer output power: Device search faculty member in the Department of Pedi- concept, modeling, and preliminary studies,” IEEE Transactions on atrics at Baylor College Medicine and a member Ultrasonics, Ferroelectrics, and Frequency Control, vol. 54, no. 11, pp. of the Predictive Analytics Laboratory at Texas 2318–2330, 2007. Children’s Hospital. His current research interests [34] M. Bedard ´ and A. Berry, “Development of a directivity-controlled include biomedical imaging and advanced machine learning for medical piezoelectric transducer for sound reproduction,” Journal of Sound and monitoring of critically ill patients. Vibration, vol. 311, no. 3-5, pp. 1271–1285, 2008. [35] Y. Cao, Q. Chen, H. Zheng, L. Lu, Y. Wang, and J. Zhu, “Study on the mechanism of ultrasonic power measurement sensor based on http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Computing Research Repository arXiv (Cornell University)

Solvability for Photoacoustic Imaging with Idealized Piezoelectric Sensors

Acosta, Sebastian
Computing Research Repository , Volume 2021 (2002) – Feb 23, 2020
Download PDF
10 pages

Loading next page...
 
/lp/arxiv-cornell-university/solvability-for-photoacoustic-imaging-with-idealized-piezoelectric-c9ufdZbvpX

References (54)

ISSN
0885-3010
eISSN
ARCH-3344
DOI
10.1109/TUFFC.2020.3005037
Publisher site
See Article on Publisher Site

Abstract

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 field 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 [11]–[13], on Fabry–Perot the case for piezoelectric sensors that are commonly employed for photoacoustic imaging. In this paper, using the method interferometry [14]–[17] or on fiber-optic refractometry [18]– of matched asymptotic expansions and the basic constitutive [20]. In [21] we formulated a model specifically 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 finite size of sensing elements. See [12], [13], [22]–[25] 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 film I. I NTRODUCTION is flat 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 films with thick- modality designed to advantageously combine the high ness 30 – 100 m approximately [11], [26]–[29]. 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 [1]–[4]. (d) For the piezoelectric film, 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 film is mechanically isotropic. sensors has been recognized in [5]–[8]. 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 field and its normal derivative at multiple paddings designed to match the mechanical the boundary. With that motivation, Dreier and Haltmeier [9] impedance of the acoustic medium [30]–[32]. 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 film, 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 [10] derived Fourier-based reconstruction follows models described in [26], [27], [32]. 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 field (Dirichlet data) is measured exactly. These film has a small uniform thickness  > 0, constant density assumptions are not satisfied 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, film is denoted . The interface between the piezoelectric film USA. (sebastian.acosta@bcm.edu) 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 film, @ u = rp (3) p p where the pressure p is defined as p = ( + 2) div u: (4) Combining (3) and (4), we find that the pressure field p satisfies the wave equation, 2 2 @ p = c p (5) p p t p where the wave speed c is defined by c = ( + 2) = . p p The piezoelectric transducer measures the electrical voltage V across the piezoelectric film 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 [27]. Under small and wave speed c . The interface between the acoustic domain b perturbations of field conditions, the linearized constitutive and the piezoelectric film is denoted by . The interface relation for the piezoelectric effect is the following between the piezoelectric film 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 field 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 [27], 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 influence 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 define 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 field 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 definition of strain s in terms of film 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 first and second Lame ´ coefficients. 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 simplified 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 [11], [26], [27], [30], [32]–[35]. 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 definition of voltage as an electric potential, the voltage Piezoelectric coeff. d 3 – 15 pC/N V generated across the piezoelectric film by the 3-component Coefficient  0.3 – 1.5 of the electric displacement D satisfies 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 film 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 film and into the backing layer . We assume that p b @ V / @ p c  p (13) p ? p t t p the pressure field 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 field is constant across the piezoelectric film. This where H is the mean curvature of the surface . See [36], [37] assumption is rigorously justified in the next section. for a derivation. The symbol  appearing in the model (13) is a unitless As in [21], we make some geometric assumptions about the coefficient defined by the elastic and piezoelectric properties domain occupied by the piezoelectric film. We let p p of the sensing film fy 2 : 0 < dist(y; ) < g. For sufficiently small , the e e 2 (d d ) domain can be expressed as a family of parallel surfaces ? ? p = = parametrized by 0 < z <  and defined 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 sufficiently small , the surfaces are well- 1  (1 2d =d) defined, 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 fluoride concerning parallel surfaces in [38, Sect. 6.2]. (PVDF) piezoelectric sensors. The transmission of the pressure field from the acoustic We note from (13) that a theoretically perfect transduction domain into the piezoelectric film is governed by the following from pressure to voltage would be attained if the coefficient 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 coefficients 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 film, respectively. The first condition in order for  = 0, the Poisson’s ratio would have to be  = 0:5 (16) ensures the continuity of the pressure field. 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 films 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 film 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 reflections of the waves. Here we seek to model how between the pressure fields and the piezoelectric sensor. For the sensors exert influence on the pressure waves. This model an introduction to this method, see [39]. First we consider the 4 formal asymptotic expansions for the pressure fields, coefficient satisfies 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 find that the piezoelectric measurements satisfy the following directivity pattern and introduce a change of variable in order to extract the effect of the piezoelectric film thickness , V cos = 1 + 1  sin  (26) z =   for  2 [0; 1]: (19) p cos  + c inc The boundary value problem for the pressure field p in the p where  is given by (14). Figure 2 displays the directional piezoelectric film, 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 coefficient  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 first 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 field in the domain a function of the incidence angle  and coefficient . The parameters correspond to a PVDF film 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 [8], [21], [40], corresponds to water with wave speed c = 1500 m/s and [41]. 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 flat 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 define the forward problem of pho- inc a reflection governed by (23). The total pressure field p is the toacoustic imaging in terms of the wave propagation model superposition of the incident and reflected 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 field p is assumed to satisfy the following initial value problem, where R is the reflection coefficient, k is the reflection 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 find that the reflection 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 [8], [44], [46] 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 satisfies 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 [42]. We also for the pressure field 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 [21] 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 profile 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 define the following operation We work with the standard Sobolev spaces based on square- integrable functions defined 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 sufficiently 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 [42], [43]. 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 field p satisfying The following lemma is a well-established result. See [42, (27), find 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 profile of the wave speed c and the time T < 1. These k 1k the field 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 flow. We work with the manifold kuk k 1k  CkVk 1 0 (34) C ([0;T ];H ()) H ([0;T ];H ()) 2 2 endowed with the Riemannian metric c dx . See [44], [45] holds for some constant C > 0. for details. 1 0 The definition of the Bochner spaces H ([0; T ]; H ()) and Assumption 1 (Nontrapping condition): Let be a simply k 1k C ([0; T ]; H ()) can be found in [42], [43]. 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 [44], [45] 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 definition 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 [48] for an excellent analysis of the Landweber method for a similar problem. The iterative process is defined 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 definite 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 final boundary value problem, 2 2 @ ' c ' = 0 (0; T ) (37a) @ ' c @ ' + H' =   (0; T ) (37b) b n t b ' = 0 and @ ' = 0 ft = Tg (37c) where  solves 2 2 2 @  = @ c  in (0; T ) (38a) t t p = 0 and @  = 0 on ft = Tg (38b) in order to define the adjoint mapping as F : 7! @ 'j : (39) t t=0 It can be shown, using straight-forward integration by parts, that indeed F is the adjoint of F or equivalently that hf; @ 'j i = hp; i t t=0 [0;T ] 2 2 = hp; I c  @ i ? [0;T ] p t = hV; i [0;T ] for all sufficiently smooth f and , where V = F (f ) according to (28) and @ 'j = F ( ) according to (39). t t=0 The brackets h;i denote the inner product of the Hilbert space H (X ) which extends as the duality pairing between the Sobolev functionals and functions. Algorithm 1 Accelerated Landweber iteration Set 0 < K , 0 < < 2kFk and 0   < 1. Initial guesses u = F V and v = u 0 0 0 for k = 1; 2; :::; K do v = u (F Fu u ) =ku k k k1 k1 0 0 Fig. 3: A: Coarse mesh for the FEM method. B: Exact profile u = v +  (v v ) k k k k1 to be reconstructed showing the brain vasculature imaged with MRI technology [49]. C: Synthetic piezoelectric measurements return u as modeled by (28). 7 Both the forward map F and its adjoint F 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 I in order to have c = 1, diam( ) = 2 and  = 1. The final time T = 2. The non-dimensional parameters of the piezoelectric film have been chosen as follows  = 1:5, c = 1:0. The p p parameters of the backing layer are  = 2:0 and c = 1:0. b b The piezoelectric-elastic coupling coefficient  = 0:9. These non-dimensional parameters are consistent with the ranges of their dimensional counterparts listed in Table I. Figure 3 displays a coarse mesh used for the FEM, the exact pressure profile to be reconstructed and the boundary measure- ments. These measurements were synthetically generated by applying the discrete version of the forward operator F 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 T = 2. 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. 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 1% relative error, the original Landweber method ( = 0:0) takes 36 iterations, whereas the accelerated method with  = 0:6 takes 12 iterations. For values  > 0:7, some instability is observed. =0.0 =0.2 10% =0.4 =0.6 1% Fig. 5: Error profiles for synthetic measurements obtained using the piezoelectric model (28). A: Reconstruction using an interpretation of the measurements as piezoelectric data. The relative error is 0:54%. B: Reconstruction using a naive interpretation of the measurements as Dirichlet data. The relative error is 6:43%. 0.1% 0 10 20 30 40 50 Iteration the same synthetic measurements are used. For the first Fig. 4: Relative error versus iteration number, for various reconstruction, we properly interpret the given measurements values of the momentum factor  to accelerate the Landweber as generated by the piezoelectric model (28) and carry out iterations. The relaxation factor = 5 10 in all cases. the reconstruction using Algorithm 1. After 50 iterations we In order to visualize the impact of improperly modeling obtain a relative error of 0:54%. For the second reconstruction, the piezoelectric measurements, we have implemented two we improperly interpret the measurements as the Dirichlet reconstructions of the initial acoustic profile. In both case, data of the pressure field. This is the naive model commonly Rel Error 8 K k XX employed by others but inconsistent with the piezoelectric = F c (FF ) V (40) kj transduction. The reconstruction is carried out using Algorithm k=0 j=0 1 modified by setting  = 0 which is equivalent to assuming for some coefficients c . This last expression motivates the kj that the measurements are Dirichlet data. After 50 iterations study of how the discrete normal operator (FF ) responds to we obtain a relative error of 6:43%. The error profiles for both the noise contained in the piezoelectric measurements. Figure reconstructions are displayed in Figure 5. 7 displays the average power spectral response to time- and space-tracings of white noise as the input to (FF ). We observe that the operator (FF ) suppresses the high frequency white noise components. This explains why the reconstruction algorithm pink noise can handle white noise better than pink or red noise. The latter red noise noise types have a larger portion of their power residing over 10% low frequencies. Therefore, the reconstruction algorithm does not suppress those types of noise. -1 1% 0 5 10 15 20 25 -2 Iteration Fig. 6: Relative error versus iteration number for white, pink -3 and red noise. In all cases the noise level is at 10%. -4 The effect of noise added to the piezoelectric measurements space is also explored. We have considered three types of noise: time white, pink and red. These are characterized by increasing -1 0 1 2 3 10 10 10 10 10 autocorrelation distances or equivalently faster decay of their Normalized Frequency spectral power. White or Gaussian noise has equal power spectral density (PSD) at different frequencies. Pink noise has Fig. 7: Frequency response of the normal operator (FF ) to a PSD decaying like ! as ! ! 1. Red or Brownian noise white noise as the input. has a PSD decaying like ! . Figure 6 displays the relative error as a function of the iteration number for the three types of noise at a 10% level. We observe that the reconstruction VII. D ISCUSSION AND CONCLUSIONS method can handle best the white noise. The error for the We have proposed a mathematical model for the acous- pink noise is greater. And the reconstruction error for the red tic measurements transduced by piezoelectric sensors. This noise is the greatest of the three types of noise. model, taking the form (28), is highly idealized as described in This phenomenon could be explained by studying the spec- the Introduction. This idealization allows us to attain a balance tral characteristics of the discrete normal operators (F F ) between correctness and simplicity in order to analyze the or (FF ). The iterations from the Landweber algorithm properties of the model (see section III) and the solvability of correspond to truncated Neumann series [48], the associated photoacoustic problem (see sections IV-V). The directional response for plane waves derived in sec- tion III and displayed in Figure 2 matches the experimental u = (I F F ) (F V ); measurements carried out by other researchers [11]–[13], [30]. k=0 Design considerations for the mechanical and electrical prop- where binomial formula yields erties of the piezoelectric film (condensed in the parameter defined in (14)) play a role in the appearance of critical angles k j (I F F ) = ( F F ) : where the sensitivity of the sensor vanishes. The incorporation j=0 of this type of sensor response into reconstruction algorithms Therefore, the Landweber iterate u could be expressed as has been highlighted as one of the challenges associated with photoacoustic imaging [7], [24], [50]–[52] and partially K k XX investigated in [6], [8]–[10], [21]. Our present paper represents u = c (F F ) (F V ) K kj a novel contribution to this research effort. k=0 j=0 Rel Error Spectral Power 9 We note that our mathematical model for the piezoelectric an impedance Z such that Z < Z < Z . Z so that the ml ml p b sensor is qualitatively similar to the model for the Fabry– wave field experiences a more gradual change of media across Perot transducer proposed in [21] in spite of the completely the fluid-sensor interface. different physical principles from which they are derived. The validity of the proposed asymptotic model is limited to Hence, a common mathematical framework can be used to small values for the thickness  of the piezoelectric film with analyze both types of sensors. From the theoretical perspective, respect to the wavelength of the acoustic waves. Advanced as stated in Theorem 1, we conclude that the photoacoustic industrial processes are able to manufacture piezoelectric films with thickness in the range 30 – 100 m approximately. As imaging problem is solvable for piezoelectric measurements. shown in Table I, the wave speed in PVDF materials ranges However, some stability is lost compared to measuring directly from 1300 – 2300 m/s. Hence, we expect our model to be the Dirichlet data. valid for frequencies . 12 MHz. For higher frequencies, the We have also proposed an iterative reconstruction Algorithm wave field is affected by the presence of the piezoelectric 1 based on an accelerated Landweber method. We imple- film and the thin matching layer [30]. In such a scenario, mented proof-of-concept synthetic simulations to highlight the our transmission model would be inappropriate. Hence, there importance of incorporating the proper modeling of the piezo- remains a need for a more complete model that can accurately electric transducer. We carried out reconstructions with and handle multiple frequency scales. without the proper piezoelectric model. See the error profiles shown in Figure 5. We observe that certain features, whose ACKNOWLEDGMENT wavefronts may have reached the measurement boundary at a The author would like to thank Texas Children’s Hospital non-normal incidence, are missed by the naive reconstruction. for its support and for the research oriented environment We also investigated the effect of noise of different spectral provided by the Predictive Analytics Laboratory. The author characteristics. The reconstruction method can handle white gratefully acknowledges support by the NSF grant DMS- noise better than pink or red noise. See Figures 6. This 1712725. The author also thanks the anonymous reviewers is due to suppression of high-frequency components in the whose observations were most constructive. measurements. See Figure 7. One limitation of the proposed method relates to the numer- REFERENCES ical characteristics of the FEM and Newmark time-stepping. [1] X. Wang, Y. Pang, G. Ku, X. Xie, G. Stoica, and L. V. Wang, Notice in Figure 4 that the error initially decays exponentially “Noninvasive laser-induced photoacoustic tomography for structural and as the theory for the Landweber method predicts. However, functional in vivo imaging of the brain,” Nature Biotechnology, vol. 21, as the iterations continue, the error eventually stagnates. This no. 7, pp. 803–806, 2003. [2] P. Beard, “Biomedical photoacoustic imaging,” Interface Focus, vol. 1, stagnation may be attributed to the fact that the measurements pp. 602–31, 2011. were synthetically generated on a mesh different from the [3] K. Wang and M. Anastasio, “Photoacoustic and Thermoacoustic To- reconstruction mesh. Among other issues, the discrete version mography: Image Formation Principles,” in Handbook of Mathematical Methods in Imaging, O. Scherzer, Ed. Springer New York, 2011, pp. ofF may not be an exact adjoint for the discrete version ofF . 781–815. Also, the numerical scheme to solve the wave equation suffers [4] L. V. Wang and S. Hu, “Photoacoustic tomography: In vivo imaging from organelles to organs,” Science, vol. 335, pp. 1458–1462, 2012. from numerical dispersion. Different frequency components [5] G. Johansson and A. J. Niklasson, “Approximate dynamic boundary of the initial pressure profile travel at different group velocity conditions for a thin piezoelectric layer,” International Journal of Solids towards the detection boundary. As a consequence, the discrete and Structures, vol. 40, no. 13-14, pp. 3477–3492, 2003. [6] D. Finch, “On a thermoacoustic transform,” in Proc. 8th Int. Meeting version of (F F ) loses coercivity (becomes ill-conditioned) on Fully 3D Image Reconstruction in Radiology and Nuclear Medicine, and the reconstructed images suffer from aberration. Reme- 2005, pp. 150–151. dies for this phenomenon have been investigated, including [7] L. V. Wang and X. Yang, “Boundary conditions in photoacoustic tomography and image reconstruction,” Journal of Biomedical Optics, regularization, two-grid methods and numerical schemes with vol. 12, no. 1, p. 014027, 2007. lower dispersion. See [45, Sect. 6.8-6.10], [53] and references [8] S. Acosta and C. Montalto, “Multiwave imaging in an enclosure with therein. However, these improvements fall outside of the scope variable wave speed,” Inverse Problems, vol. 31, no. 6, p. 065009, 2015. of this paper. [9] F. Dreier and M. Haltmeier, “Explicit inversion formulas for the two- dimensional wave equation from Neumann traces,” arXiv:1905.03460, As future research, it would be very important to explicitly [10] G. Zangerl, S. Moon, and M. Haltmeier, “Photoacoustic tomography model the influence of matching layers commonly employed with direction dependent data : An exact series reconstruction approach,” in the design of piezoelectric sensors [30], [32]. Matching Inverse Problems, vol. 35, p. 114005, 2019. layers play an important role in optimizing the transmission of [11] V. Wilkens and W. Molkenstruck, “Broadband PVDF membrane hy- high-frequency acoustic energy into the sensor. In this paper, drophone for comparisons of hydrophone calibration methods up to 140 MHz,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control, vol. 54, no. 9, this transmission is modeled by (27b) which for a flat surface pp. 1784–1791, 2007. simplifies to [12] K. A. Wear, C. Baker, and P. Miloro, “Directivity and Frequency- 1 @p 1 @p Dependent Effective Sensitive Element Size of Needle Hydrophones: + = 0: Predictions From Four Theoretical Forms Compared With Measure- @n Z @t ments,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Fre- Here Z =  c is the acoustic impedance of the thick backing quency Control, vol. 65, no. 10, pp. 1781–1788, 2018. b b b [13] ——, “Directivity and Frequency-Dependent Effective Sensitive Ele- substrate which usually does not match the acoustic impedance ment Size of Membrane Hydrophones: Theory versus Experiment,” Z = c of the fluid or the acoustic impedance Z =  c of p p p IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Con- the piezoelectric film. Matching layers are designed to have trol, vol. 66, no. 11, pp. 1723–1730, 2019. 10 [14] P. Beard, F. Perennes, and T. Mills, “Transduction mechanisms of pyroelectric effect,” Acoustical Physics, vol. 64, no. 6, pp. 789–795, the Fabry-Perot polymer film sensing concept for wideband ultrasound 2018. detection,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Fre- [36] X. Antoine, H. Barucq, and A. Bendali, “Bayliss-Turkel-like radiation quency Control, vol. 46, no. 6, pp. 1575–1582, 1999. conditions on surfaces of arbitrary shape,” Journal of Mathematical [15] B. T. Cox and P. C. Beard, “The frequency-dependent directivity of Analysis and Applications, vol. 229, pp. 184–211, 1999. a planar Fabry-Perot polymer film ultrasound sensor,” IEEE Trans. [37] S. Acosta, “High order surface radiation conditions for time-harmonic Ultrason. Ferroelectr. Freq. Control, vol. 54, no. 2, pp. 394–404, 2007. waves in exterior domains,” Computer Methods in Applied Mechanics and Engineering, vol. 322, pp. 296–310, 2017. [16] E. Zhang, J. Laufer, and P. Beard, “Backward-mode multiwavelength [38] R. Kress, Linear Integral Equations, 2nd ed., ser. Applied mathematical photoacoustic scanner using a planar Fabry-Perot polymer film ultra- sound sensor for high-resolution three-dimensional imaging of biological sciences. New York : Springer-Verlag, 1999, vol. 82. tissues,” Applied Optics, vol. 47, no. 4, pp. 561–577, 2008. [39] E. Hinch, Perturbation Methods. Cambridge Univ. Press, 2012. [40] S. Acosta and C. Montalto, “Photoacoustic imaging taking into account [17] J. A. Guggenheim, E. Z. Zhang, and P. C. Beard, “A method for thermodynamic attenuation,” Inverse Problems, vol. 32, no. 11, p. measuring the directional response of ultrasound receivers in the range 115001, 2016. 0.3-80 Mhz using a laser-generated ultrasound source,” IEEE Trans. [41] S. Acosta, “Recovery of pressure and wave speed for photoacoustic Ultrason. Ferroelectr. Freq. Control, vol. 64, no. 12, pp. 1857–1863, imaging under a condition of relative uncertainty,” Inverse Problems, vol. 35, pp. 46–49, 2019. [18] G. Wild and S. Hinckley, “Acousto-ultrasonic optical fiber sensors: [42] L. C. Evans, Partial Differential Equations, ser. Graduate Studies in Overview and state-of-the-art,” IEEE Sensors Journal, vol. 8, no. 7, Mathematics. Providence, RI: American Mathematical Society, 1998, pp. 1184–1193, 2008. vol. 19. [19] K. A. Wear and S. M. Howard, “Directivity and frequency-dependent [43] J.-L. Lions and E. Magenes, Non-homogeneous boundary value prob- effective sensitive element size of a reflectance-based fiber-optic hy- lems and applications, ser. Grundlehren der mathematischen Wis- drophone: Predictions from theoretical models compared with mea- senschaften in Einzeldarstellungen, Bd. 181-183. Berlin, New York, surements,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Fre- Springer-Verlag, 1972. quency Control, vol. 65, no. 12, pp. 2343–2348, 2018. [44] C. Bardos, G. Lebeau, and J. Rauch, “Sharp sufficient conditions for [20] G. Wissmeyer, M. A. Pleitez, A. Rosenthal, and V. Ntziachristos, the observation, control, and stabilization of waves from the boundary,” “Looking at sound: optoacoustics with all-optical ultrasound detection,” SIAM Journal on Control and Optimization, vol. 30, no. 5, pp. 1024– Light: Science and Applications, vol. 7, p. 53, 2018. 1065, 1992. [21] S. Acosta, “Well-posedness for photoacoustic tomography with Fabry- [45] R. Glowinski, J.-L. Lions, and J. He, Exact and Approximate Control- Perot sensors,” SIAM Journal on Imaging Sciences, vol. 12, no. 4, pp. lability for Distributed Parameter Systems : A Numerical Approach, 1669–1685, 2019. ser. Encyclopedia of Mathematics and its Applications. Cambridge [22] M. Xu and L. V. Wang, “Analytic explanation of spatial resolution University Press, 2008, vol. 117. related to bandwidth and detector aperture size in thermoacoustic or [46] P. Stefanov and G. Uhlmann, “Thermoacoustic tomography with variable photoacoustic reconstruction,” Physical Review E, vol. 67, p. 056605, sound speed,” Inverse Problems, vol. 25, p. 075011, 2009. [47] H. Engl, M. Hanke, and A. Neubauer, Regularization of Inverse Prob- [23] M. Haltmeier and G. Zangerl, “Spatial resolution in photoacoustic lems, ser. Mathematics and Its Applications. Springer, 2000, vol. 375. tomography: Effects of detector size and detector bandwidth,” Inverse [48] P. Stefanov and Y. Yang, “Multiwave tomography with reflectors: Problems, vol. 26, no. 12, p. 125002, 2010. Landweber’s iteration,” Inverse Problems & Imaging, vol. 11, no. 2, [24] K. Wang, S. Ermilov, R. Su, H. Brecht, A. Oraevsky, and M. Anastasio, pp. 373–401, 2017. “Imaging model incorporating ultrasonic transducer properties for three- [49] W. Koroshetz, J. Gordon, A. Adams, A. Beckel-Mitchener, J. Churchill, dimensional optoacoustic tomography,” IEEE Transactions on Medical G. Farber, M. Freund, J. Gnadt, N. S. Hsu, N. Langhals, S. Lisanby, Imaging, vol. 30, no. 2, pp. 203–214, 2011. G. Liu, G. C. Peng, K. Ramos, M. Steinmetz, E. Talley, and S. White, [25] H. Roitner, M. Haltmeier, R. Nuster, D. P. O’Leary, T. Berer, G. Paltauf, “The state of the NIH brain initiative,” Journal of Neuroscience, vol. 38, H. Grun, ¨ and P. Burgholzer, “Deblurring algorithms accounting for the no. 29, pp. 6427–6438, 2018. finite detector size in photoacoustic tomography,” Journal of Biomedical [50] B. Cox, J. Laufer, and P. Beard, “The challenges for quantitative Optics, vol. 19, no. 5, p. 056011, 2014. photoacoustic imaging,” Proc. SPIE, vol. 7177, pp. 1–9, 2009. [26] B. Fay, G. Ludwig, C. Lankjaer, and P. A. Lewin, “Frequency response [51] R. Ellwood, E. Zhang, P. Beard, and B. Cox, “Photoacoustic imaging of PVDF needle-type hydrophones,” Ultrasound in Medicine and Biol- using reflectors to enhance planar arrays,” Journal of Biomedical Optics, ogy, vol. 20, no. 4, pp. 361–366, 1994. vol. 19, p. 126012, 2014. [27] J. Sirohi and I. Chopra, “Fundamental understanding of piezoelectric [52] J. Xia, J. Yao, and L. V. Wang, “Photoacoustic tomography: Principles strain sensors,” Journal of Intelligent Material Systems and Structures, and advances,” Progress In Electromagnetics Research, vol. 147, pp. vol. 11, no. 4, pp. 246–257, 2000. 1–22, 2014. [28] Y. Su, F. Zhang, K. Xu, J. Yao, and R. Wang, “A photoacoustic [53] E. Zuazua, “Propagation, observation, and control of waves approxi- tomography system for imaging of biological tissues,” Journal of Physics mated by finite difference methods,” SIAM Review, vol. 47, no. 2, pp. D: Applied Physics, vol. 38, pp. 2640–2644, 2005. 197–243, 2005. [29] M. A. Oreilly and K. Hynynen, “A PVDF receiver for ultrasound mon- itoring of transcranial focused ultrasound therapy,” IEEE Transactions on Biomedical Engineering, vol. 57, no. 9, pp. 2286–2294, 2010. [30] L. F. Brown, “Design considerations for piezoelectric polymer ultra- sound transducers,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control, vol. 47, no. 6, pp. 1377–1396, 2000. [31] J. Tichy, ´ J. Erhart, E. Kittinger, and J. P´ ıvratska, ´ Fundamentals of Sebastian Acosta received a B.Sc. degree in me- piezoelectric sensorics: Mechanical, dielectric, and thermodynamical chanical engineering and a M.Sc. degree in math- properties of piezoelectric materials. Springer-Verlag Berlin Heidel- ematics from Brigham Young University, Provo, berg, 2010. Utah, USA in 2009 and 2011 respectively, and a [32] T. L. Szabo, “Transducers,” in Diagnostic Ultrasound Imaging: Inside Ph.D. degree in computational and applied mathe- Out, 2nd ed. Academic Press, 2014, ch. 5, pp. 121–165. matics from William Marsh Rice University, Hous- [33] B. Zeqiri, P. N. Gelat, ´ J. Barrie, and C. J. Bickley, “A novel pyroelectric ton, Texas, USA in 2014. He is currently a re- method of determining ultrasonic transducer output power: Device search faculty member in the Department of Pedi- concept, modeling, and preliminary studies,” IEEE Transactions on atrics at Baylor College Medicine and a member Ultrasonics, Ferroelectrics, and Frequency Control, vol. 54, no. 11, pp. of the Predictive Analytics Laboratory at Texas 2318–2330, 2007. Children’s Hospital. His current research interests [34] M. Bedard ´ and A. Berry, “Development of a directivity-controlled include biomedical imaging and advanced machine learning for medical piezoelectric transducer for sound reproduction,” Journal of Sound and monitoring of critically ill patients. Vibration, vol. 311, no. 3-5, pp. 1271–1285, 2008. [35] Y. Cao, Q. Chen, H. Zheng, L. Lu, Y. Wang, and J. Zhu, “Study on the mechanism of ultrasonic power measurement sensor based on

Journal

Computing Research Repository – arXiv (Cornell University)

Published: Feb 23, 2020

There are no references for this article.