New Frontiers in Computational Mathematics Inverse Problems in Electromagnetics Bill Lionheart Department of Mathematics UMIST, Manchester, UK 11 January 2004
Electromagnetic Inverse Problems In Inverse Boundary Value Problems we attempt to infer material properties in the interior of a domain from measurements at the boundary. In electromagnetics the material properties include conductivity, permittivity, and permeability which could be scalar or tensor functions of space. Applications include • Electrical Impedance Tomography for medical diagnosis • Permittivity of conductivity imaging of industrial processes • Geophysical exploration and monitoring • Non-destructive testing – eddy current imaging/testing, microwave tomography, photoelas- tic tomography, director field measurement in LCDs 1
Some challenges of Electromagnetic Inverse Problems Typically we are considering inverse boundary value problems . We have a (system of) PDEs with unknown spatially varying coefficients but with multiple over-determined boundary data . Challenges include • They are ill-posed specifically we are trying to solve F ( x ) = d for x ∈ R n the unknown pa- rameters, d ∈ R m the data where, although F is smooth and injective, the inverse of F is unstable with respect to inversion. In the continuum problem F is a compact operator in any reasonable spaces. • These problems are inherently non linear. For example in the conductivity equation ∇ · ( σ∇φ ) = 0 for given boundary data (eg Dirichlet), φ also depends on σ so the over-specified data (eg Neumann) is a non linear function of σ . For problems with a large contrast non-linearity is important. • Iterative solution methods must be employed. This means the forward problem, typically a FEM calculation must be done repeatedly with for several fixed sets of boundary data but with varying coefficients . 2
Maxwell’s equations Time harmonic Maxwell’s equations ∇ × E = − i ω H (1) ∇ × H = ( σ + i ωε ) E (2) Together with ∇ · µ H = 0 , ∇ · ε E = 0 . Any of the material properties σ , ε and µ might be anisotropic (symmetric matrix valued functions). Typical boundary data would be all (theoretically) or some (practically) pairs of tangential com- ponents of E and H at the boundary. 3
Typical scalar inverse problem Electrical Impedance Tomography used in geophysics, process monitoring, medical diagnosis. There are uniqueness results for the idealized problem ∇ · σ∇φ = 0 For σ scalar or symmetric tensor valued function. The boundary data is pairs ( φ | ∂Ω , σ∂φ / ∂ n ) or equivalently the Dirichlet-to-Neumann map Λ σ : φ | ∂Ω �→ σ∂φ / ∂ n We seek to find σ from Λ σ . Typically σ is permittivity or conductivity. Under various smoothness hypotheses σ �→ Λ σ is known to be injective (scalar σ ). Of course inverse Λ σ �→ σ not continuous in any reasonable topology. 4
More realistic boundary conditions Current is applied via a discrete set of metallic electrodes E l with a thin resistive layer caused by electro-chemistry. On the l -the electrode E l φ + z l σ∂φ ∂ n = V l where V l is constant on E l , and z l is the contact impedance. The total current on E l is σ∂φ Z I l = ∂ n dS E l and on the parts of the boundary ∂Ω of the domain Ω with no elecrodes, we have no current σ∂φ / ∂ n = 0 . These non-local boundary conditions are non-standard, and one reason we have to make our own custom FEM code to solve the forward problem. 5
Solving the Inverse Problem For the EIT case x is a discrete parameterization of σ and d is the set of voltages on each electrode, when a set of a set current patterns are applied. To get a stable solution we must add a prori information about x . Typical method is to use a penalty term and minimize � F ( x ) − d � 2 + G ( x ) Were a common choice for the penalty term is G ( x ) = � Lx � 2 for a discretized differential operator L . This has the advantage of being a smooth optimization problem to which we can apply Gauss- Newton or non-linear conjugate gradient for example. The resulting solution will be smooth by assumption. Note the Bayesian interpretation is a prior probability density e − G ( x ) / 2 and uniform Gaussian noise, in which case our minimum is the Maxi- mum a-posteriori estimate. (details later) 6
Some smooth regularization, Gauss-Newton solutions on experimental data. reconstruction by Nick Polydorides and experimental results in collaboration with Robert West at the University of Leeds. 7
Probabilistic interpretation of regularization The statistical approach to regularization gives an alternative justifi cation of generalized Tikhonov regularization. Bayes’ theorem relates conditional probabilities of random variables. For simplicity consider the linear inverse problem of solving Ax = b P ( x | b ) = P ( b | x ) P ( x ) . P ( b ) The probability of x given b is the probability of b given x times P ( x ) / P ( b ) . We now want the most likely x , so we maximize P ( x | b ) , the so called Maximum A-Posteriori (MAP) estimate. This is easy to do if we assume x is multivariate Gaussian with mean x 0 and covariance Cov [ x ] = P − 1 , and e has mean zero and Cov [ e ] = Q − 1 � � � � 1 − 1 − 1 2 � Ax − b � 2 2 � x − x 0 � 2 P ( x | b ) = · exp P ( b ) exp Q P where we have used that x and e are independent so P b ( b | x ) = P e ( b − Ax ) . We notice that P ( x | b ) is maximized by minimizing � Ax − b � 2 Q + � x − x 0 � 2 P . where � x � P = � Px � . 8
Beyond simple Regularization The assumption that the noise in the data is an independent identically distributed Gaussian variable is reasonable in some situations – for example a carefully calibrated optics experiment averaged over many measurements. In that case 2-norm minimization is sensible. In other situations, such as geophysical data gathered in the field, or clinical EIT data there are many sources of correlated error, some systematic. It is often more realistic to use a norm which rejects outliers. Often the MAP estimate is not so useful – for example when the posterior distribution is bimodal. In any case we don’t usually just want a picture but error bars at least. One way to explore the entire posterior distribution is to use Markov Chain Monte Carlo Method (MCMC). This involves a very large number of forward problem solutions and is currently only computationally feasible for small problems. 9
Forward solution issues • Rapid assembly of FEM system matrix, geometric information precalculated and system matrix updated as conductivity parameters change. Another reason for custom solver. • For some applications σ is complex, system matrix is complex symmetric • Have multiple right-hand-side vectors. Although iterative methods can still be used to good effect as previous potentials φ can be used to start iteration with updated σ . • Both multi-grid and adaptive meshing can improve accuracy and speed. 10
Less smooth conductivity In medical applications organs boundaries give rises to jump discontinuities, as do phase bound- aries in industrial applications. One approach is to replace the smooth regularization term a Total Variation term Z | ∇σ | this works well but is difficult to implement efficiently, see Nick Polydorides’ talk next. Other approaches are possible including shape reconstruction (including level set methods), Linear Sampling method and monotonicity method. All have advantages and disadvantages. 11
Anisotropic problems Materials may be anisotropic because of • Layered or fibrous structure (rocks, muscle, composite materials) • Deformation – permittivity in photoelastic effect, compression of soil or rock. • Crystaline structure , eg calcite, or liquid crystals • Alignment of particles in fluid flow. For the low frequency case, a gauge condition, most easily seen by interpreting the conduc- tivity as a Riemannian metric, means that the solution is not unique and extra information is needed. However there are high frequency problems where the permeability is isotropic, but the permittivity is anisotropic, and there is no problem with gauge conditions. 12
Photoelastic Tomography One such application is Photoelastic Tomography. The permittivity is a known linear function of strain, one attempts to measure the permittivity tensor using polarized light. Currently, the technique is restricted to two dimensional objects, and models must be stress frozen and then sliced, which is expensive and time consuming. In a collaboration with Rachel Tomlinson and Joab Winkler at Sheffield, Hanno Hammer and I are developing a 3D photoelastic tomography system to overcome these problems. For small stresses one can assume ray propagation of light, and the linearized problem can be reduced to scalar Radon transforms, and inverted with conventional X-Ray CT algorithms. For higher stresses we hope to use a Newton-Kantarovich method. In contrast to many inverse problems where the regularization is ad hoc in this problem we can use the divergence of the stress as a penalty term, as this should vanish in the interior where there are no applied forces. 13
Photoelastic measurement system. Courtesy Rachel Tomlinson, University of Sheffi eld 14
Two dimensional photoelasticity. Turbine, thread on nut and bolt, loaded crane hook. Courtesy Rachel Tomlinson, University of Sheffi eld 15
Recommend
More recommend
