DETERMINATION OF BURN DEPTH ON THE BASIS OF SKIN SURFACE TEMPERATURE – SOLUTION OF INVERSE PROBLEM USING THE GRADIENT METHOD

The burned and healthy layers of skin tissue are considered. The temperature distribution in the domains is described by the system of two Pennes equations with the different thermophysical parameters. In the healthy layer the metabolic and perfusion heat sources are taken into account, while the burned layer is dead and the blood perfusion and metabolic do not occur in this region. At the surface between burned and healthy layers the ideal contact is assumed (continuity of heat flux and temperature field), at the internal surface limiting the system the body temperature is known. Heat transfer between skin surface and environment is described by the well known Robin boundary condition (ambient temperature and heat transfer coefficient are given). It is assumed that the shape of surface between burned and healthy tissue is unknown. Additional information necessary to solve the inverse problem formulated results from a knowledge of skin surface temperature distribution. At the stage of direct problem solution the multiple reciprocity boundary element method is used. This variant allows one to avoid the discretization of domain interior (inside of the healthy tissue sub-domain the volumetric internal heat sources must be taken into account). To solve the inverse problem above formulated the gradient method is used and the shape sensitivity analysis is applied to determine the coefficients appearing in the least square criterion. Here the implicit variant of shape sensitivity analysis is used. In the final part of the paper the results of computations are shown.


INTRODUCTION
A very important problem in the burns therapy is a correct evaluation of skin burn depth [1,2,3,4].One of the possibilities depends on the measurements of the temperature distribution at the skin surface (thermographs) and the solution of the problem using an inverse approach.So, the burned and healthy layers of skin tissue should be considered.The temperature distribution in the domains is described by the system of two Pennes equations supplemented by the boundary conditions.The position of surface between burned and healthy tissues is un-known.Additional information necessary to solve the identification problem results from a knowledge of skin surface temperature distribution.At the stage of direct problem solution the classical boundary element method for burned subdomain and multiple reciprocity boundary element method for healthy tissue sub-domain are used.These methods are coupled by the boundary condition given at the contact surface between sub-domains considered.This approach allows one to avoid the discretization of domain interior (inside the healthy tissue sub-domain the volumetric internal heat sources must be taken into account).To solve the inverse problem formulated the gradient method is used.The least square criterion is applied in which the sensitivity coefficients appear.To determine these coefficients the implicit variant of shape sensitivity analysis is introduced.In the final part of the paper the results of computations are shown.

GOVERNING EQUATIONS
Let us consider the domain formed by two sub-domains (Figure 1).: λ ( ) ( ) 0 where T 2 (x) is the tissue temperature, λ 2 is the tissue thermal conductivity, W B is the blood perfusion rate, c B is the specific heat of blood, T B is the arterial blood temperature, Q met is the metabolic heat source, x are the spatial co-ordinates, x={x 1 , x 2 } for 2D problem, x={x 1 , x 2 , x 3 } for 3D problem.
For burned tissue, blood perfusion and metabolic heat generation are equal to zero, because the tissue is dead.So where λ 1 is the thermal conductivity of burned tissue.
At the surface between sub-domains the continuity of heat flux and temperature field is assumed where T e /n , e=1, 2 denotes the normal derivative, n is the normal outward vector, for 2D problem: n= [cos where α is the heat transfer coefficient, T a is the ambient temperature.At the internal surface the body core temperature is known For all others boundaries the zero heat flux can be assumed.

Boundary integral equation for burned and healthy tissue
The direct problem has been solved using the boundary element method [5,6].The boundary integral equation corresponding to the equation ( 2) is the following where ξ is the observation point, the coefficient B(ξ) is dependent on the location of source point ξ, T * ( ξ, x) is the fundamental solution, and then where r is the distance between the points ξ and x, while For the healthy tissue sub-domain the multiple reciprocity boundary element method [7,8] has been applied.This variant of the BEM allows one to avoid the discretization of domain interior.So, the following integral equation corresponding to the equation ( 1) is considered where where The heat fluxes for 3D problem

Numerical realization of the BEM
To solve the equations ( 6) and ( 10) the boundary Γ is divided into N elements Γ j =1,2,…,N.
Next, the integrals in the equations ( 6), ( 10) can be replaced by the sums of integrals over these elements, leading to the following equations where N 1 is the number of elements distinguished at the boundary limiting domain Ω 1 .
It should be pointed out that different types of boundary elements can be used, namely constant elements, linear or parabolic ones [5,6].Here the linear boundary elements are applied [5. 6].
After the mathematical manipulations one obtains the following system of algebraic equations corresponding to the burned tissue and the system of equations corresponding to the healthy tissue where K 1 is the number of boundary nodes located at the boundary limiting sub-domain Ω 1 , KK 1 is the number of boundary nodes located at the boundary limiting sub-domain Ω 2 .For example, for 2D problem presented in Figure 2: K 1 =48, K=108.It should be pointed out that taking into account the boundary conditions, at the contact surface between sub-domains Ω 1 and Ω 2 (nodes 2869 in Figure 2) the double boundary nodes should be introduced.The remaining double nodes (shown in Figure 2) for example 3-4, 24-25 correspond to different boundary conditions, namely nodes 4, 25 are connected with the Robin condition, while nodes 3, 25 are connected with the Neumann condition.
The systems of equations ( 17), ( 18) can be written in the matrix form  G q H T P . (20) The way of calculation of matrix elements G 1 , H 1 , G 2 , H 2 , P is described in details in [7].
For the needs of further considerations concerning the temperature field calculation the following denotations are introduced (c.f.Figures 1 and 2) T T T q q q are the vectors of functions T and q at the boundary T T q q are the vectors of functions T and q on the contact surface Γ c between sub-domains Ω 1 and Ω 2 , T T q q q are the vectors of functions T and q at the boundary Using above notations, one obtains the following systems of equations  for the burned region -for the healthy tissue domain The condition (3) written in the form should be introduced to the equations ( 21), ( 22).Next, coupling these system of equations one has The remaining boundary conditions should be also taken into account.So where A is the main matrix of the system of equations ( 25), Y is the vector of unknowns and B is the free terms vector.
The system of equations ( 26) allows one to find the ,missing'' boundary values.Knowledge of boundary temperatures and heat fluxes at all nodes constitutes a basis for determination of internal temperatures at the optional points selected from the domain considered [7,8].

SHAPE SENSITIVITY ANALYSIS
In literature one can find two basic approaches to sensitivity analysis using the boundary element method: the continuous approach and the discretized one [9].In the case of continuous approach (explicit differentiation method) the mathematical model of sensitivity is formulated and next the solution is found numerically using the BEM.The implicit differentiation method, which belongs to the discretized approach, basis on the differentiation of algebraic boundary element matrix equations.
In the paper the implicit differentiation method is applied.Let us assume that b is the shape parameter [9,10].The system of equations ( 25) should be differentiated with respect to parameter b and then It should be pointed out that the derivatives of the boundary element matrices are calculated analytically [11].The additional problem defined by the system of equations (28) has been solved using the BEM.The main matrix of the systems of equations ( 26) and ( 28) is the same.

INVERSE PROBLEM
The rectangular domain of dimensions 2LL with internal boundary Γ c shown in Figure 1 is considered.The inverse problem formulated here bases on the assumption that the temperature distribution at the skin surface Γ ex is known (e.g.thermographs), while the position of Γ c is unknown.The 'measured' temperatures at the skin surface are denoted by T d i , i = 1,2,…M, where M is the number of sensors.
The surface Γ c is defined by the parabolic function where (L, y p ) is the parabolic vertex.The inverse problem consists in the identification of value b = y p which determines maximum burn depth in the domain considered.
The criterion which should be minimized is of the form [12] 2 1 1 ( ) ( ) where T i are the calculated temperatures.These temperatures are obtained from the direct problem solution with an estimate for unknown value of b.
Using the necessary condition of minimum, one obtains 1 ( ) 0 where are the sensitivity coefficients, k is the number of iteration, for k = 0 b 0 is the arbitrary assumed value of parameter b, while b k for k > 0 results from the previous iteration.It should be pointed out that in order to determine the coefficients (32) the shape sensitivity analysis described in chapter 4 is used.Function T i is expanded in the Taylor series about known value of b k , this means 1 () Introducing (33) into (31) one obtains where K is the assumed number of iterations.

RESULTS OF COMPUTATIONS
The domain of dimensions 0.04m 0.02m has been considered.At first, the direct problem described in Chapter 1 has been solved.The following input data have been assumed: thermal conductivity of burned tissue λ 1 =0.1 W/(mK), thermal conductivity of healthy tissue λ 2 =0.2 W/(mK), blood perfusion rate W B =0.5 kg/((m 3 s), specific heat of blood c B =4200 J/(kgK), arterial blood temperature T B =37 o C, metabolic heat source Q met =200 W/m 3 , heat transfer coefficient α=10 W/ (m 2 K), ambient temperature T a =20 o C. The discretization of boundaries of burned and healthy tissue sub-domains using the linear boundary elements is shown in Figure 2.
In Figure 3 the temperature distribution at the skin surface for three different values of parabolic vertex b (c.f. Figure 2) this means 0.012 m, 0.014 m and 0.016 m is shown.The last curve corresponds to the constant burn depth, of course.Next, the inverse problem has been solved using the values of temperatures at the skin surface (nodes 4-24 shown in Figure 2) obtained from the solution of direct problemc.f.

CONCLUSIONS
The possibilities of burn shape estimation on the basis of knowledge of skin surface temperature are shown.To this end the mathematical model basing on the system of two Pennes equations for burned and healthy tissue has been formulated and least square criterion has been applied.The inverse problem has been solved by gradient method coupled with the BEM.In future, real values of skin temperatures obtained from experiments will be applied.
from the fundamental solutions (11) can be calculated analytically and then

Figure 3 .
The parameter b has been estimated in iterative way (equation (35)).In Figures4, 5, 6 the results of inverse problems solution are shown.These Figures illustrate the values of parameter b for successive iterations and the curves correspond to the different initial values of parameter b 0 .In the all cases presented the iteration process was convergent and the exact value of b has been obtained after the several iterations.

Figure 4 .
Figure 4. Results of identification for different initial values of b 0 (b real =0.012m).

Figure 5 .
Figure 5. Results of identification for different initial values of b 0 (b real =0.014m).

Figure 6 .
Figure 6.Results of identification for different initial values of b 0 (b real =0.016m).