Depto. de Análisis Matemático y Matemática AplicadaFac. de Ciencias MatemáticasFALSEunpu

Carpio, Ana

Rapún, ML

Carpio Rodríguez, Ana María

Rapún Banzo, María Luisa

Docta Complutense

Domain and parameter reconstruction inphotothermal imagingAna Carpio and Marı́a–Luisa RapúnAbstract In this work we address the inverse problem of reconstructing inclusionsand their thermal parameters given temperature measurements at the accessible sideof a material. We describe an iterative descent method that combines topologicalderivative computations to reconstruct the geometry of the defects with gradientiterations to approximate the material parameters. A numerical experiment showingthe ability of the method to obtain reasonable reconstructions in a few iterations ispresented.1 Statement of the problemPhotothermal imaging techniques are suitable means of inspecting composite mate-rials with nondestructive tests. In this work we develop techniques to detect defectsΩ buried in a medium by surface thermal measurements. We are interested in a pho-tothermal technique that consists in heating the surface of a semi–infinite mediumby a laser beam and recording the temperature at several receptors located on thesame surface during a time interval, see Fig. 1. Recent physical experiments usingthis kind of technique can be found in [9, 18].The forward problem is modelled by a heat diffusion equation in the half planeR2− := {(x,y) ∈ R2, y < 0}. The surface of the sample Π := {(x,0), x ∈ R} isthermically excited with a delta–pulse located at a source point x0 ∈ Π , generatinga thermal wave of the formAna CarpioDep. Matemática Aplicada, Fac. Matemáticas, Universidad Complutense de Madrid, 28040Madrid, Spain e-mail: ana carpio@mat.ucm.esMarı́a–Luisa RapúnDep. Fundamentos Matemáticos, ETSI Aeronáuticos, Universidad Politécnica de Madrid, 28040Madrid, Spain e-mail: marialuisa.rapun@upm.es12 Ana Carpio and Marı́a–Luisa Rapúnsource point ΠΩ ΩFig. 1 Geometrical configuration. The source and observation points are located on the boundaryΠ of the medium. The objects Ω and their physical parameters are the unknowns of the phothermalimaging problem.Uinc(x, t) = (1/t) exp(−ρe|x−x0|2/(4κet)), x ∈ R2, t > 0. (1)Here κe is the thermal conductivity of the exterior medium Ωe := R2− \Ω and ρeis the density multiplied by its specific heat. The corresponding thermal parametersinside the inclusions Ω are κi and ρi. The temperature distributionU(x, t) :={U+(x, t), in Ωe × (0,∞),U−(x, t), in Ω × (0,∞),satisfies the heat equationsρe∂tU+ = κe∆U+, in Ωe × (0,∞), ρi∂tU− = κi∆U−, in Ω × (0,∞). (2)In the exterior domain Ωe, the total temperature Utotal = U+ +Uinc is the super-position of U+ and the incident wave defined in (1). The temperature satisfies thefollowing transmission conditions at the common interface:U−−U+ =Uinc, κi∂nU−−κe∂nU+ = κe∂nUinc, on ∂Ω × (0,∞). (3)The forward problem is completed imposing an adiabatic boundary condition on theupper boundary Π , and homogeneous initial conditions:∂nU+ = 0, on Π × (0,∞), U+(x,0) =U−(x,0) = 0, ∀x ∈ R2−. (4)The solution U of the forward problem can be numerically approximated usingthe following strategy [12, 14]: if we consider the Laplace transform of U and Uinc,u(x,s) =∫ ∞0 e−stU(x, t)dt and uinc,s(x) =∫ ∞0 e−stUinc(x, t)dt, then, for each value ofs the function us(x) := u(x,s) is a radiating solution of the stationary problemDomain and parameter reconstruction in photothermal imaging 3κe∆us − sρeus = 0, in Ωe, κi∆us − sρius = 0, in Ω ,u−s −u+s = uinc,s, on ∂Ω κi∂nu−s −κe∂nu+s = κe∂nuinc,s, on ∂Ω ,∂nus = 0, on Π . (5)To invert the Laplace transform we choose the hyperbolic paths of the form [17]:γ(θ) := µ(1− sin(π/4+ ıθ)), θ ∈ R, where µ > 0. Then, the solution of (2–4) isU(x, t) =12πı∫ ∞−∞etγ(θ)u(x,γ(θ))γ ′(θ)dθ .Numerical approximations of U can be calculated using a truncated trapezoidal ruleU(x, t)≈L∑ℓ=−Lcℓetsℓu(x,sℓ),sℓ = γ(log(L)L ℓ)and weights cℓ =log(L)2πıL γ′(log(L)L ℓ).The inverse problem consists in finding the objects Ω and the parameters κi, ρisuch that the solution of the forward transmission problem (2–4) equals the mea-sured values of the total wave Umeas(xi, t j) at the detector locations x1, . . . ,xM ∈ Πat the time instants t1, . . . , tN . Since this problem is ill–posed, we consider a weakervariational reformulation: find Ω ,κi,ρi minimizing the functionalJ(R2− \Ω ,κi,ρi) =12M∑i=1N∑j=1f (t j)(Utotal(xi, t j)−Umeas(xi, t j))2 , (6)where Utotal is the solution of the forward problem (2–4) when the object is Ω andthe interior thermal parameters are κi and ρi. The weight function f (t) normalizesthe time decay of the solutions of the heat equation. For our numerical experimentin section 3 we select f (t) = maxx∈{x1,...,xM} |Umeas(x, t)|−1. Other possibilities wereexplored in [4, 6].Based on the Laplace–transform strategy described above, we proposed in [4, 6]to substitute the cost functional (6) by the approximated funcionalJ(R2− \Ω ,κi,ρi) =12M∑i=1N∑j=1f (t j)(L∑ℓ=−Lcℓet jsℓu(xi,sℓ)−Umeas(xi, t j))2, (7)having now 2L+1 stationary constraints.2 Iterative method to reconstruct inclusions and parametersTo solve the optimization problem we combine gradient and topological derivative(TD in the sequel) methods to generate sequences of parameters and objects in such4 Ana Carpio and Marı́a–Luisa Rapúna way that the cost functional decreases throughout the iterative procedure. Ourchoice of a TD strategy is based on the following advantageous features:• Without a priori information, the TD provides a good first guess of the num-ber, size and location of the inclusions. This has been tested in a wide range ofphysical settings, including acoustics, electromagnetism, elastodynamics, electri-cal impedance tomography, fluorescence optical tomography, and photothermalimaging [1, 5, 7, 8, 10, 15, 20].• Iterative TD methods allow for topological changes during the iterations, in con-trast to classical shape deformation strategies [11,13,19] that require knowledgeof the number of objects from the start. Using iterative TD based methods, newobjects may be created in the course of the iterations, existing contours maymerge and holes inside existing objects may be detected, see [2,3]. Furthermore,even if the number of inclusions is known (assumption that in most practicalapplications is not realistic), TD–iterative methods are a powerful alternative tothese classical methods, providing accurate reconstructions at a low computa-tional cost, as extensively checked by the authors in different contexts (see [2–6]and references therein).• In comparison with other strategies allowing for topological changes (as i.e. levelset methods [16, 21]), the number of iterations with respect to the domain isusually much smaller.In our previous papers [4, 6] we used a non–standard formulation of the pho-tothermal problem (2–4), involving two interior parameters related with κe,κi,ρe,and ρi with no physical meaning. In this paper we adapt the results in [4, 6] to dealwith the reconstruction of defects and of their physical parameters κi and ρi.The TD of a shape functional J (R) is a pointwise function defined as [22]:DT (x,R) = limε→0J (R \Bε(x))−J (R)πε2, x ∈ R, (8)where Bε(x) is a ball centered at x with radius ε . Then, it follows the expansion:J (R \Bε(x)) = J (R)+DT (x,R)πε2 +o(ε2), as ε → 0.This motivates the key idea for the reconstruction technique: if we locate smallobjects Bε(x) at the points x∈R where DT (x,R) is negative, then J (R \Bε(x))<J (R), that is, the value of the functional decreases. Hence we will identify thepoints where the TD attains the larger negative values with the regions where it ismore likely to have an object.The next result can be proved following Theorem 3.2 in [4].Theorem 1. The TD of the functional J(R2− \Ω ,κi,ρi) defined in (7) isDT (x) = Re( L∑ℓ=−L2κe(κe −κi)κe +κi∇utotal,sℓ(x)∇psℓ(x)+(ρe −ρi)sℓutotal,sℓ(x)psℓ(x))(9)Domain and parameter reconstruction in photothermal imaging 5for x ∈ R2− \Ω , where utotal,sℓ = uinc,sℓ +usℓ and usℓ is the solution of (5) for s = sℓ.The adjoint fields psℓ are solutions of:κe∆ psℓ − sℓρe psℓ = gsℓ , in Ωe, κi∆ psℓ − sℓρi psℓ = 0, in Ω ,p−sℓ − p+sℓ = 0, on ∂Ω κi∂n p−sℓ −κe∂n p+sℓ = 0, on ∂Ω ,∂n psℓ = 0, on Π ,(10)with gsℓ(x) := ∑Mi=1 ∑Nj=1 f (t j)csℓet jsℓ(Umeas(xi, t j)−∑Lk=−L cket jsk usk(xi))δxi(x).By iteratively applying Theorem 1, we construct a monotone sequence of ap-proximate domains Ωd ⊂ Ωd+1 adding to the current approximation Ωd the pointswhere the TD attains pronounced negative values. To be able to remove points fromΩd , we need to compute the TD inside the inclusions. The definition (8) can be ex-tended to the points inside Ω [3], and an analogous expression to (9) can be foundfor x ∈ Ω [5, 6]. This extension is the basis to develop iterative strategies able tocorrect an approximation of Ω by removing the points where the TD attains pro-nounced positive values.To determine the thermal parameters we proceed as follows. If κ̃i, ρ̃i are ap-proximate values of κi and ρi, and Ω̃ is an approximation of Ω , then we correctthe values κ̃i, ρ̃i by a gradient method. The idea is to define κi = κ̃i +ηϕ , ρi =ρ̃i +ηψ , where η > 0 is small and ϕ ,ψ are selected calculating the derivative ofJ(η) := J(Ω̃ , κ̃i +ηϕ , ρ̃i +ηψ) with respect to η to ensure that J′(0)< 0. A proce-dure to obtain explicit formulae in terms of forward and adjoint fields for this kindof functionals is explained in [3, 6]. In our case, it can be proven that the choiceϕ = Re(∫Ω̃L∑ℓ=−L∇usℓ∇psℓ), ψ = Re(∫Ω̃L∑ℓ=−Lsℓ usℓ psℓ), (11)makes J′(0)< 0. Here usℓ , psℓ are solutions of (5) and (10), respectively, with Ω =Ω̃ , κi = κ̃i and ρi = ρ̃i.Finally, our procedure is as follows. In a first step we consider initial guesses ofthe parameters κi = κ0i , ρi = ρ0i and compute the TD in R2− for these parameters,that is, the TD of J(R2−,κ0i ,ρ0i ). We find then a first approximation Ω1 of Ω as theunion of all the points where the TD is smaller than a negative constant (see [3,4] forguidelines of the selection of such constant). Once Ω1 is set, we update the valuesof the parameters performing Q iterations of the gradient method (Q = 8 in ournumerical example in Section 3) as explained above:κqi = κq−1i +ηϕq, ρqi = ρq−1i +ηψq, q = 1, . . . ,Q,with ϕ q,ψq defined as in (11) with Ω̃ = Ω1, κ̃i = κq−1i and ρ̃i = ρq−1i . Once theparameters are corrected, we compute the TD of J(R2− \Ω1,κQi ,ρQi ) to update thedomain Ω1 by adding to it the points x ∈ R2− \Ω1 where the TD attains the largernegative values, and removing from Ω1 the points inside it, if any, where the TDattains the larger positive values. Once the approximation of the domains is im-6 Ana Carpio and Marı́a–Luisa Rapún(a) (b)(c) (d)Fig. 2 (a) Topological derivative when κ0i = 3/4, ρ0i = 1/3 and Ω = /0. (b)–(d) Approximateddomains Ωd , d = 1,2,3, superimposed to the TD computed for Ω = Ωd and the correspondingupdated values of the thermal parameters κi and ρi.proved we perform further gradient iterations to update the parameters and so on.The algorithm stops if any of the following stopping criteria is satisfied:• meas(Ωd \Ωd−1) is small,• |κqi −κq−1i |+ |ρqi −ρq−1i | is small and ∥Umeas −Utotal∥ is small,• J(R2− \Ωd ,κqi ,ρqi ) is small.3 A numerical exampleIn this section we present a numerical example to illustrate the feasibility of ourreconstruction algorithm. We consider a simple geometry where Ω is the ellipseΩ = {(x,y) ∈ R2, x2/0.552 +(y+ 1)2/0.352 < 1}, with thermal parameters κi =1/2 and ρi = 1. In the exterior medium the values of the parameters are κe = 1 andρe = 1/5.Synthetic data are created solving (2–4) by means of the Laplace transform withrespect to time and a boundary element formulation in space (see [4, 12] for de-tails). A relative 1% Gaussian error was added at each observation point to bothavoid inverse crimes and to simulate measurement errors. We have considered sixincident waves of the form (1) generated at the uniformly distributed source pointsrepresented in all the plots in figure 2 by ’•’ marks. Measurements of the temper-ature were taken at the seven observation points represented by ’×’ marks at 10uniformly distributed times in the time interval [0.05,0.5].We started the algorithm by choosing the initial values κ0i = 3/4 and ρ0i = 1/3.To obtain an initial guess for the domain, we computed the TD in the samplingregion [−3,3]× [−2.5,0]. This yields the colormap represented in figure 2(a). DarkDomain and parameter reconstruction in photothermal imaging 7(a) (b)−3 −2 −1 0 1 2 3−2.5−2−1.5−1−0.500 10 20 30 40 50 60 70 80 900.50.60.7κ iiteration0 10 20 30 40 50 60 70 80 900.40.60.81ρ iiterationFig. 3 Final reconstrucion after 10 iterations with respect to the domain. (a) Initial (green dashedline), predicted (red dashed line) and true (blue solid line) objects. (b) Values of the thermal pa-rameters κi and ρi versus the number of iterations.blue colors indicate the regions where the TD takes large negative values, at whichthe objects should be located. The boundary of the true defect is represented by asolid white line in all plots in figure 2. Our initial guess Ω1 is represented in figure2(b). We set now Ω = Ω1 and perform eight iterations with the gradient methodto correct the values of κi and ρi. In the next step, the values of the parameters arefixed and the TD is again computed. In figure 2(c) we represent the updated objectΩ2. The hybrid algorithm alternates eight iterations with the gradient method withone TD computation. It stopped at the tenth iteration with respect to the domain.The first three approximated domains are represented in figure 2. In figure 3(a) weshow the true object (solid blue line), the initial guess Ω1 (dashed green line), andthe final reconstruction Ω10 (dashed red line). The values of κi and ρi throughoutthe iterations are given in figure 3(b). Two identical values mark a TD computationto update the domain. The final approximations were κ f inali = 0.5703 (recall that thetrue value is 0.5), and ρ f inali = 0.8266 (while the true value is 1). We have obtained asatisfactory reconstruction taking into account that no a priori knowledge about thenumber, size or location of the objects is assumed, and that few data were available.In our example, we found a sequence of enlarging sets, i.e., satisfying Ωd ⊂Ωd+1. An example where the TD provides a sequence of defects where Ωd * Ωd+1for some values of d ∈ N can be found in [6]. The interested reader may find somereconstructions with other geometries, multiple objects and different weight func-tions f (t) in [4, 6]. Furthermore, a gallery of comparisons varying the different pa-rameters of the problem, namely, the number of source points and/or observationpoints, the number of time observations, etc, can be found in [4] for a simplifiedsituation where the interior parameters are assumed to be known and in [6] for arelated problem with unknown domains and parameters.Acknowledgements The authors are partially supported by the Spanish Government researchproject TRA2010–18054 and the Spanish Ministerio de Economı́a y Competitividad Grants No.FIS2011–222888–C02–02 and FIS2010–22438–E.8 Ana Carpio and Marı́a–Luisa RapúnReferences1. Bonnet, M., Constantinescu, A.: Inverse problems in elasticity. Inverse Problems, vol. 21, pp.1–50 (2005)2. Carpio, A., Rapún, M.–L.: Topological derivatives for shape reconstruction, In: Inverse Prob-lems and Imaging, Lecture Notes in Mathematics, pp. 85–134. Springer (2008)3. Carpio, A., Rapún, M.–L.: Solving inhomogenous inverse problems by topological derivativemethods. Inverse Problems vol. 24, art. num. 045014 (2008)4. Carpio, A., Rapún, M.–L.: Domain reconstruction using photothermal techniques. J. Comput.Phys. vol. 15, pp. 8083–8106 (2008)5. Carpio, A. Rapún, M.–L.: Hybrid topological derivative and gradient based methods for non–destructive testing. Abstract and Applied Analysis 2013, art. num. 816134 (2013)6. Carpio, A., Rapún, M.–L.: Parameter identification in photothermal imaging. J. Math. ImagingVis. vol. 49, pp. 273–288 (2014)7. Feijoo, G.R.: A new method in inverse scattering based on the topological derivative. InverseProblems, vol. 20, pp. 1819–1840 (2004)8. Garreau, S., Guillaume, P., Masmoudi, M.: The topological asymptotic for PDE sistems: theelasticity case. SIAM J. Control Optim., vol. 39, pp. 1756–1778 (2001)9. Garrido, F., Salazar, A.: Thermal wave scattering by spheres. J. Appl. Phys. vol. 95, pp. 140–149(2004)10. Guzina, B.B., Bonnet, M.: Small–inclusion asymptotics of misfit functionals for inverse prob-lems in acoustics. Inverse Problems, vol. 22, pp. 1761–1785 (2006)11. Hettlich, F.: Fréchet derivatives in inverse scattering problems. Inverse Problems vol. 11, pp.371–382 (1995)12. Hohage, T., Sayas, F.–J.: Numerical approximation of a heat diffusion problem by boundaryelement methods using the Laplace transform. Numer. Math. vol. 102, pp. 67–92 (2005)13. Kirsch, A.: The domain derivative and two applications in inverse scattering theory. InverseProblems, vol. 9, pp. 81–96 (1993)14. Laliena, A., Sayas, F.–J.: LDBEM in diffusion problems. In: Proceedings of XIX CEDYA/IXCMA (electronic version) (2005)15. Laurain, A., Hintermüller, M., Freiberger, M., Scharfetter, H.: Topological sensitivity analysisin fluorescence optical tomography. Inverse Problems, vol. 29, art. num. 025003 (2013)16. Litman, A., Lesselier, D., Santosa, F.: Reconstruction of a two–dimensional binary obstacleby controlled evolution of a level–set. Inverse Problems, vol. 14, pp. 685–706 (1998)17. López-Fernández, M., Palencia, C.: On the numerical inversion of the Laplace transform ofertain holomorphic mappings. Appl. Numer. Math. vol. 51, pp. 289–303 (2004)18. Mendioroz, A., Castelo, A., Celorrio, R., Salazar, A.: Characterization of vertical buried de-fects using lock–in vibrothermography: I Direct problem. Meas. Sci. Technol. vol. 24, art. num.065601 (2013)19. Potthast, R.: Fréchet differentiability of the solution to the acoustic Neumann scattering prob-lem with respect to the domain. Journal of Inverse and Ill–Possed Problems, vol. 4, pp. 67–84(1996)20. Samet, B., Amstutz, S., Masmoudi, M.: The topological asymptotic for the Helmholtz equa-tion. SIAM J. Control Optim. vol. 42, pp. 1523–1544 (2003)21. Santosa, F.: A level–set approach for inverse problems involving obstacles. Optimisation etCalcul des Variations, vol. 1, pp. 17–33 (1996)22. Sokolowski, J., Zolésio, J.P.: Introduction to shape optimization. Shape Sensitivity Analysis.Springer, Heidelberg (1992).

Spanish

Domain and parameter reconstruction in photothermal imaging

http://eprints.ucm.es/40463/1/CarpioRapunECMI2014.pdf

Domain and parameter reconstruction in photothermal imaging

Abstract

Similar works

Full text

Available Versions

Docta Complutense