Studies On The Numerical Methods Of The Inverse Problems In EIT And MREIT

Electrical impedance tomography (EIT) is a technique to recover spatial prop-erties of the interior of a conducting object for biologic tissue imaging from electro-static measurements taken on its boundary. Magnetic resonance electrical impedance tomography (MREIT) is to recover the tissue conductivity using the internal cur-rent density and the induced magnetic flux density information. Compared with EIT, MREIT has some improvement in the data acquisition and the reconstruction method. At present, EIT and MREIT have attracted extensive attention worldwide in medical imaging, geophysics and so on. In this paper, we propose two numeri-cal algorithms for the cavity reconstruction problem in EIT, and the conductivity reconstruction in Bz based MREIT. We specify the results as follows:I. Cavity reconstruction within homogeneous medium in EITWe consider a disc BR (?) R2 with its radius R. Assuming that the object is homogeneous and conducting except for a insulating cavityΩ, which is convex with smooth boundary (?)Ω, and we assume that there exists a constant r0>0, satisfying Bro={(x,y)\x2+y2≤r02}(?)Ω(?)BRThere is a prescribed boundary current f∈L02((?)BR), the corresponding boundary potential g∈L02((?)BR) can be measured without physical damage to the object, which is satisfied the NtD map associated with f:ΛN-D f=g((?)BR). Then the potential u satisfies the Laplace equation with boundary conditions n is the outer (relative to BR\Ω) unit normal vector.Cauchy problem (40),(41) has a unique solution u∈H01(BR\Ω) if it exists[128]. The idea of the algorithm is for given boundary data (41), to solve the Laplace equation (40), then to recover (?)Ωdepending on the cavity boundary condition (42).Ⅰ.1. Cauchy problem in circular ring domain We transform (40),(41) into the Cauchy problem under the polar coordinates as: and n is written as Then, the cavity boundary condition (42) is transformed as follows: Thus Then the cavity problem can be characterized as:solving r0(θ)∈C2π1, satisfies where u is the solution to (43). By the assumption ofΩ, there exists a constant r0, satisfying 0<r0<r0(θ)< R, then there exists two domains satisfying: D=(r0(θ),R]×[0,2π](?) D0=[r0,R]×[0,2π] We consider the Cauchy problem in Do: From the analytic continuation, the solutions to the Cauchy problem (43) and (46) are equal to each other in D[148], i.e. u=v, in D. we solve the Cauchy problem (46), and hope to find the point set where the Neumann data of the solution vanish, that is the boundary of the cavityΩ.By the variables separation, we have a general expansion of the solution to the Cauchy problem (46) as follows: and the boundary conditions can be given as where a0, a'0, an, a'n, bn, b'n are given. Comparing the superposition coefficients in (47), (48) and (49), we have Then the solution (47) to the Cauchy problem (46) can be written as:where Pn,Qn,Sn,Wn are given in (50). Moreover, we haveFinally, our problem is characterized as:Solving r0=r0(θ)∈C2π1, satisfyingⅠ.2. Optimization procedureWe use the functional minimization method to obtain an approximate solution. Forwe construct a functional asDefine represented the finite dimensional space of trigonometric polynomial with order of less than or equal to N-1,cn is the complex coefficients.From the compound trapezoid formula by the uniform subdivision to [0,2π] with the stepΔθ, we yield the discrete expression to the functional as: Then the aim is to determine a set of expansion coefficient{c0,c1,…,cN-1),such that the functional can reach the minimum.Let then the nonlinear least-squares problem can be described as: where G:RN→RM is a nonlinear function on{cn}.Denote S({cn})to be the Jacobian matrix of G on{cn} Then the gradient of Jh({cn})is and the Hessian matrix of Jh({cn})is The difficulty is that the second-order term T({cn})in the Hessian matrix H({cn}) is hard to calculate.To avoid this,we adopt Gauss-Newton(G-N)algorithm by ignoring the second-order term to solve the nohlinear least-squares problem(57): If we choose appropriate parameterλ,and make Then the iteration procedure can be ensured decreasing(not increasing at least).The experiments demonstrate that G-N method is always valid if the initial value is not far from the local point of minimum to the nonlinear least-squares problem,with the requirement that the matrix S({cn*})is not a singular matrix.But unfortunately,S({cn*})is always singular in fact,such that the algorithm converge to other points,while we should re-choose the decreasing direction.We now consider the following trust region model: where hk is called trust region radius. The above model can be characterized by solving the following equation: Then,we have If letμk=0;otherwiseμk>0. Since the matrix(S({cn(k)})T S({cn(k)})+μkI)is positive defined,the direction solved by(62)is decreasing.This algorithm is called Levenberg-Marquardt(L-M)method,with the L-M parameterμk.Ⅱ.Successive linearization method for nonlinear Bz-based MREITWe propose a numerical algorithm for Bz-based MREIT in view of the basic idea of Newton iteration method:successive linearization. We transform the non-linear operator equation into a series of linear problems,we also deduce the Frechet derivative,and give the equivalent expression of the solution to the linear problems. Ⅱ.1. Abstract operator equationWe consider a smooth and simply connected domainΩ(?)R3, with the isotropic conductivity distributionσ. The impressed electrical current I through the elec-trodes E±attached on the boundary (?)Ωresults an internal current density J= (Jx,Jy,Jz).Consider the perfect electrodes, then there is no tangential current through the electrodes, i.e.J×n|'E±=0. From the Ohm's law J=-σ▽u, the resulting voltage u satisfies the following boundary value problem: From the Biot-Savart's law. the induced magnetic flux density B=(Bx,By,Bz) and the current density J satisfy whereμ0 is the vacuum permeability. Then Bz-based MREIT is to recover the tissue conductivity a using z-component of the internal magnetic flux density information.The basic idea of MREIT forΩ(?)R3 is to recoverσin every two dimensional section perpendicular to z-direction by assuming a and the prescribed boundary cur-rent to be invariant along z-direction, then the resulting voltage u is invariant along z-direction either. Therefore, we consider the two dimensional MREIT problem only.LetΩz0=Ω∩(x,y,z):z=z0}(?)R2 be the cross section perpendicular to z-direction, consider the boundary value problem inΩz0: where n is the outer unit normal vector of boundary (?)Ωz0.On the other hand, due to the prescribed boundary current being invariant along z-direction, the internal electrical current density J can be regarded as a plane current density in every two dimensional cross section perpendicular to z-direction, the induced internal magnetic flux density is perpendicular to z-direction too, therefore, the z-component of the magnetic flux density can be expressed as:Thus, the two dimensional MREIT problem can be characterized as:solving the scalar binary functionσ(x,y), to satisfy (65) and (66).Consider the function space C(Ωz0) consisting of continuous function inΩz0. Forσ∈C(Ωz0), define a nonlinear operator F:C(Ωz0)(?)C(Ωz0): where u(x, y)=u(x, y;σ) is the solution to the equation (65).Therefore, we transform the inverse MREIT problem into the nonlinear oper-ator equation in C(Ωz0) where B=Bz(x,y) is a known function. II.2. Successive linearizationWe will transform the nonlinear operator equation (68) into a series of linear problems in view of the basic idea of Newton iteration method:successive lineariza-tion. Letσ0∈(Ωz00) is the given initial approximation, make iteration until the sequence{σn} is convergence.Δσn is the solution to the linear problem where DF(σn) is the Frechet derivative of F onσn.For givenσn∈C(Ωz0), solvingΔσn is the key of the iterative algorithm, so we consider the following problem:for givenσ, g∈C(Ωz0), solving p, to satisfy For anyσ,u∈C(Ωz0), we define an operator in two variables from the definition of the operator F(a), we have where u(σ) is the solution to the equation (65).Theorem 1 For given a, g∈C(Ωz0),ρbeing the the solution to the equation (71) is equavalent toρand v satisfying the following equation: and the operator DσG(σ,u) and DuG(σ,u) are defined as: where u(x,y)=u(x,y;σ) is the solution to the equation (65). theorem 1 has transformed a nonlinear equation into the coupling of a linear partial differential equation and a linear integral equation, it's the most important part of the iterative method. We give the specific procedure as follows:Step1:For a given initial approximative conductivityσ0,we solve the forward problem(65), and obtain the corresponding voltage u0;Step2:Giving an initial increment p0,0 of the conductivity for givenσ0 and u0, solve the equation(74), and obtain an intermediate inviable v0,0, then solve the integral equation (75) using v0,0,σ0 and u0, and we can obtain p1,0;Step3:Fixing n=0,1,2,…, we do iteration for k=1,2,…as follows:Step(1):Solve the boundary value problem (74) usingσn, un andρk,n, and obtain an intermediate inviable vk,n;Step (2):Solve the integral equation (75) using vk,n,σn and un, and obtainρk+1,n; Step (3):Setting an error∈1,if the iteration stops, otherwise, return to Step(1) by k←k+1.Step4:Letσn+1=σn+ρk+1,n, setting an error∈2,∈2, the iteration stops, otherwise, return to Step3 by n←n+1.We present numerical experiments of the unit disk with radially symmetric inclusions centered at the origin. This is a purely numerical setting with no actual measured data. For the measured data with noise, the regularization strategy can be considered. On the other hand, improving the stable algorithm for solving the coupling problem of the partial differential equation and the integral equation and proving the convergency of the algorithm are our further research.