
Citation 
 Permanent Link:
 http://digital.auraria.edu/AA00007161/00001
Material Information
 Title:
 Lowfrequency imaging algorithm for inverting arbitrary electromagnetic systems
 Creator:
 Todorovski, Dalibor J.
 Place of Publication:
 Denver, CO
 Publisher:
 University of Colorado Denver
 Publication Date:
 2019
 Language:
 English
Thesis/Dissertation Information
 Degree:
 Master's ( Master of science)
 Degree Grantor:
 University of Colorado Denver
 Degree Divisions:
 Department of Electrical Engineering, CU Denver
 Degree Disciplines:
 Electrical engineering
 Committee Chair:
 Harid, Vijay
 Committee Members:
 Golkowski, Mark
Gedney, Stephen
Record Information
 Source Institution:
 University of Colorado Denver
 Holding Location:
 Auraria Library
 Rights Management:
 Copyright Dalibor J. Todorovski. Permission granted to University of Colorado Denver to digitize and display this item for nonprofit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.

Downloads 
This item has the following downloads:

Full Text 
LOWFREQUENCY IMAGING ALGORITHM FOR INVERTING ARBITRARY
ELECTROMAGNETIC SYSTEMS
by
DALIBOR J. TODOROVSKI B.S., University of Central Florida, 2016
A thesis submitted to the Faculty of the Graduate School of the University of Colorado in partial fulfillment of the requirements for the degree of Master of Science Electrical Engineering
2019
This thesis for the Master of Science degree by Dalibor J. Todorovski has been approved for the Electrical Engineering Program by
Vijay Harid, Chair Vijay Harid, Advisor Mark Golkowski Stephen Gedney
Date: May 9, 2019
Todorovski, Dalibor J. (M.S., Electrical Engineering)
LowFrequency Imaging Algorithm for Inverting Arbitrary Electromagnetic Systems Thesis directed by Assistant Professor Vijay Harid
ABSTRACT
Methods and software were developed in MATLAB for the purposes of inverting arbitrary lowfrequency electromagnetic scattering problems and locating the equivalent current densities present in the forward problem. While the classes of problems considered for this method are those where forward scatterers are concealed within a highlyconductive halfspace or a metallic shell, the methods can be generalized to other problems due to the robustness of the electromagnetic solver used. Synthetic electromagnetic held data was created for the considered cases with the use of the commercial electromagnetic solver FEKO. Similarly, the Greenâ€™s function response was constructed numerically in this solver using the known physical boundaries of the problem and its material properties. With these computed quantities, a linear relationship was formed between the forward electromagnetic data and the numerical Greenâ€™s function (NGF) that could be inverted with the aid of various regularization techniques. Results are presented for inversions using both Tikhonov regularization and Li regularization on a full volumetric reconstruction of the imaging space as well as a more computationallyefficient iterative reconstruction.
The form and content of this abstract are approved. I recommend its publication.
Approved: Vijay Harid
m
ACKNOWLEDGMENTS
Thank you to all of the professors at the University of Colorado  Denver. You have all taught me invaluable knowledge and skills as an engineer and researcher that I will carry forward with me. And I would like to extend my gratitude to Dr. Vijay Harid. Your guidance and encouragement over these past two years have been critical in helping me find my path forward.
IV
TABLE OF CONTENTS
CHAPTER
I. MOTIVATION........................................................... 1
II. THEORY OF ELECTROMAGNETIC SCATTERING................................ 4
II. 1 Equivalent Dipole Formulation................................. 4
11.2 Scattering in FreeSpace....................................... 6
11.3 Scattering by a PEC Sphere Buried in a Homogeneous HalfSpace . 9
11.4 PEC Object Inside a Thin Conductive Shell..................... 14
III. INVERSE PROBLEM FORMALISM......................................... 16
111.1 Necessity for Regularization................................. 16
111.2 L2 Tikhonov Regularization Theory.......................... 20
111.3 Li Regularization Theory..................................... 20
111.3.1 Tibshiraniâ€™s Algorithm for Solving LASSO............. 22
111.4 Optimal Choice of Regularization Parameter................. 23
111.4.1 LCurve Method......................................... 24
111.4.2 CrossValidation Method ............................... 25
111.5 Imaging Resolution Limits ................................... 26
IV. ALGORITHM DESCRIPTION............................................. 28
IV. 1 HighLevel Overview of Algorithm............................ 28
IV.2 Generation of Synthetic Data.................................. 29
IV.3 Computation of Numerical Greenâ€™s Function (NGF)............... 32
IV.4 Inversion Algorithm........................................... 32
IV.4.1 Tikhonov Regularization................................. 32
IV.4.2 Li Regularization....................................... 33
IV.4.3 Iterative Least Squares Solution........................ 33
IV.5 Generation of Imaging Grid.................................. 34
IV.5.1 Adaptive Mesh Refinement Algorithm...................... 34
v
V. RESULTS..................................................... 36
V.l Inversion results for FreeSpace Inversion............. 36
V.2 Inversion results for Homogeneous HalfSpace Inversion. 42
V.3 Inversion results for Conductive Cubic Shell Inversion. 47
VI. FUTURE WORK................................................ 57
VII. CONCLUSION................................................ 59
REFERENCES..................................................... 60
vi
FIGURES
II. 1 Scattering of PEC sphere of radius a in freespace illuminated by a TEpolarized plane wave. In the longwavelength limit, the fields scattered by this can be approximated by the fields of a â€”zoriented ideal magnetic dipole.......................................................................... 7
11.2 Comparison of analytical scattered magnetic fields due to ideal dipole
against MoMSIE solutions of sphere scattering and equivalent magnetic dipole............................................................... 8
11.3 Configuration of PEC sphere of radius a embedded in homogeneous half
space (Region 2) with conductivity a illuminated by an incident plane wave................................................................. 9
11.4 MoMSIE solution of Js of PEC sphere buried in homogeneous halfspace
(a = 1 S/m) illuminated by TEpolarized plane wave (/ = 1kHz). A larger Js is induced at the top of sphere relative to the bottom due to attenuation over the length of the sphere inside the halfspace. 11
11.5 HMD (a = â€”y) serves as a rough approximation for the electromagnetic
scattering by a PEC sphere (centered at z = 0)........................... 13
11.6 A spherical shield with thickness A = d and outer radius a illuminated by a uniform magnetic held as in [20]. A single PEC sphere is located inside
the shield............................................................... 15
III. 1 Figure from [19]. A collection of N ideal dipoles in freespace uniformly spaced (with spacing Az) along a line of length L are used to compute the condition number of the forward operator Gb to determines its stability. 18
III. 2 Condition number of Gb is computed for a collection of dipoles N = L/Az bound on line of length L. For L = 1 m, the normalized A correspond to waves of / = 300 GHz, 300 MHz, and 300 kHz, respectively.................... 18
vii
111.3 Pseudocode describing Tibshiraniâ€™s algorithm for solving LASSO. This
minimizes the objective of the primal problem with the constraints in the equality set if........................................................ 23
111.4 Characteristic loglog plot of LCurve Method for Tikhonov Regulariza
tion. The corner or closest point to the origin represents the optimal parameter A*........................................................... 24
111.5 Equipartition of data set for 5fold CV. This provides an estimate of the
MSE of the solution vector x that can be used to optimize the regularization parameter A....................................................... 25
111.6 Minimum Spatial Resolution for signal propagating in various homoge
neous Uniteconductivity grounds. This assumes no attenuation of the signal................................................................. 26
111.7 From [19]. This suggests for an imaging region of size L m, the minimum
resolution obtainable for an ordinary least squares inversion (when the forward operator A is full rank and invertible) is ~ 30 â€” 40L mm....... 27
IV. 1 Flowchart describing highlevel design of software for both volumetric
mesh and adaptive mesh refinement procedures................................. 29
IV.2 Forward Problem Geometry for Metallic Box with two PEC Cubes. The metallic shell has a = 107 S/m and thickness A = 2mm. Receivers are placed on planes parallel to each face offset by 5mm. An ideal magnetic dipole source is placed over each face that operates in the LF regime. . . 30
IV.3 Surfaces are meshed with many triangular elements. The basis functions used to solve the SIE are mapped to the edges of these elements and used to compute surface current densities.............................................. 31
viii
IV.4 Generation of Forward Problem of PEC Sphere located in freespace handled in FEKO. The unknown PEC scatterers are placed in this model with the known problem boundaries. The scattered nearheld data are requested at Nr points..................................................... 31
IV. 5 2D Crosssection of initial coarse grid (left). Distance between grid centers
dr = \rf â€” rf'\, for j ^ i, is taken to be the size of the grid cell. At the next iteration (right), the cell with the strongest response is removed and replaced with a subgrid with 8 cells of size dr/2.................... 35
V. I Synthetic forward data of nearheld scattered data computed by FEKO
(domain limited to (x,y) E [â€”2m, 2m])................................ 37
V.2 Volumetric Inversion with Tikhonov Regularization. Optimal Inverted
Current Density \Mmv\ found at (x = â€”2.5m,y = â€”2.5m, z = â€”0.1m) . . 37
V.3 Volumetric Inversion with L\ Regularization. Optimal Inverted Current
Density \Mmv\ found at (x = â€” lm,y = â€” lm, z = â€”2m) .............. 38
V.4 Volumetric Inversion with Iterative Least Squares. Optimal Inverted Current Density \Mmv\ found at (x = â€” lm, y = â€”lm, z = â€”2m)................ 39
V.5 AMR Inversion with Tikhonov Regularization after 5 Iterations. Optimal Inverted Current Density MCTM' found at (x = â€”1.875m,y =
â€” 1.875m, z = â€”0.15937m)............................................. 40
V.6 AMR Inversion with Li Regularization after 5 Iterations. Optimal Inverted Current Density \Mmv\ found at (x = â€” 1.875m,y = â€” 1.875m, z =
0.15937m)........................................................... 41
V.7 Synthetic forward data of nearheld scattered data computed by FEKO
(domain limited to (x,y) E [â€”2m, 2m])................................ 43
V.8 Volumetric Inversion with Tikhonov Regularization. Optimal Inverted
Current Density \Mmv\ found at {x = â€” 2m, y = â€” 2m, z = â€”0.91429m) . 43
IX
V.9 Volumetric Inversion with L\ Regularization. Optimal Inverted Current
Density \Mmv\ found at (x = 2m,y = 2m,z = â€”2m) ................. 44
V.10 Volumetric Inversion with Iterative Least Squares. Optimal Inverted Current Density MCTM' found at {x = â€” 2m, y = â€” 2m, z = â€”2m).............. 45
V.ll AMR Inversion with Tikhonov Regularization after 5 Iterations. Optimal Inverted Current Density MCTM' found at {x = â€”1.875m,y =
â€” 1.875m, z = â€”0.87187m)................................................. 46
V.12 AMR Inversion with Li Regularization after 5 Iterations. Optimal Inverted Current Density \Mmv\ found at (x = â€” 1.875m,y = â€” 1.875m, z =
0.63438m)............................................................... 47
V.13 Configuration of Forward Problem with planar array of dipoles on top face
and nearheld requests below bottom face................................ 48
V.14 Synthetic forward data of nearheld scattered data computed by FEKO. 49 V.15 Volumetric Inversion with Tikhonov Regularization. Optimal Inverted
Current Density MCTM' found at {x = 2mm, y = 2mm, z = 0.14429m). . . 50
V.16 Volumetric Inversion with L\ Regularization. Optimal Inverted Current
Density MCTM' found at {x = 2mm, y = 2mm, z = 2mm).................. 51
V.17 AMR Inversion with Tikhonov Regularization. Optimal Inverted Current
Density \Mmv\ found at {x = 34.0643mm, y = 34.0643mm, z = 34.0643mm). 52 V.18 AMR Inversion with Li Regularization. Optimal Inverted Current Density
\Mmv\ found at (x = 34.0643mm, y = 34.0643mm, z = 34.0643mm). . . 53
V.19 Constrained Volumetric Inversion with Tikhonov Regularization. Optimal Inverted Current Density MCTM' found at {x = 0.4m, y = 3mm, z =
0.145m). This point lies along one of the edges of Cube 2............. 54
V.20 Current Density magnitude of constrained volumetric inversion with Tikhonov Regularization. Grid Index of maximum corresponds to optimal inverted current density MCTM'.................................................. 55
x
V.21 Constrained Volumetric Inversion with L\ Regularization. Optimal In
verted Current Density found at {x = 0.4m, y = 0.4m, z = 3mm).
This point lies along one of the edges of Cube 1........................ 55
V.22 Current Density magnitude of constrained volumetric inversion with L\ Regularization. Grid Index of maximum corresponds to optimal inverted current density \Mmv\. The effects of the Li penalty are clearly demonstrated in these sparse inversion results..................................... 56
xi
CHAPTER I
MOTIVATION
The problem of detecting features and objects of interest that are obstructed by
conductive barriers is pervasive in many engineering applications. These include subsurface imaging for geophysical exploration [8], detection of natural and manmade underground structures [26], nondestructive evaluation for vehicles and container inspections [25], underwater detection [2, 16], and several other active areas of research. The difficulty, however, is that objects of interest concealed behind highlyconductive media typically can not be located with the aid of traditional electromagnetic imaging methods. For instance, groundpenetrating (GPR) systems (10MHz < / < 1GHz) [4] typically image in soils at depths of up to 100 meters [28], but highlyconductive subsurface sediments can limit this maximum depth to a few centimeters beneath the surface.
The utility of GPR and other similar systems is greatly diminished when considering imaging through highly conductive barriers such as shipping containers. In such cases, highfrequency (> 30kHz) electromagnetic radiation cannot penetrate through metallic containers of even 1 â€” 2mm thickness without undergoing extreme losses (> 300dB). Fortunately, the drawback of high losses can be overcome with the use of extremely and very lowfrequency (ELF/VLF) electromagnetic signals (/ <30 kHz) with larger skindepths 5 in that conductive medium [32],
Here uj is the wave angular frequency, a is the conductivity, and is the permeability of free space. These electromagnetic signals will attenuate while passing through a
(
(1.1)
depth d in the conductive media as e d/s. For instance, the skin depth in aluminum
1
signals would allow for detection of the scattered fields produced by objects hidden behind these conductive barriers while minimizing these attenuation losses.
ELF/VLF signals serve as a promising alternative for the detection of objects hidden by conductive barriers. Flowever, there is an inherent tradeoff with imaging resolution. Conventional imaging techniques such as GPR, Xray tomography, or CT scans rely on the farheld (ray) properties of electromagnetic radiation to reconstruct an image. For these aforementioned imaging modalities, the imaging resolution is limited by the Rayleigh Criterion (also known as the diffraction limit) [11]. Specifically, the smallest resolvable feature is approximately half a wavelength in size, which determines the highest achievable resolution. The Rayleigh criterion is the spatial analog of the Nyquist frequency [6] and is an inherent feature of traditional image reconstruction. For example, GPR operating at 1GHz has an associated resolution limit of 15cm according to the Rayleight Criterion. The problem with ELF/VLF imaging is that the wavelengths are tens to hundreds of kilometers. The Rayleigh Criterion suggests that only objects that are several kilometers in size are resolvable using ELF/VLF radiation. This limit is reasonable when studying large scale geological features of the Earth [29]. On the other hand, imaging underground bunkers, tunnels, or objects inside shipping containers is not possible with ELF/VLF signals when using imaging algorithms that rely on the raylike features of electromagnetic waves. To circumvent this limit, it is necessary to develop novel strategies that do not rely on typical electromagnetic wave features for imaging. Imaging small objects in the ELF/VLF range will occur in the socalled nearheld region and the corresponding quasistatic nature of these signals must be used to reconstruct an image. Physically, incident ELF/VLF signals can induce loops of electric current on an object of interest which will in turn generate quasistatic magnetic fields in the nearheld of the object. These magnetic helds will have a spatial amplitude and polarization prohle that are specihc to the objectâ€™s shape and size. Thus, an appropriate ELF/VLF imaging algo
2
rithm can be designed by leveraging the mapping from the objectâ€™s geometry to the measured magnetic held profile. Even so, inferring the location of hidden objects from held data is in itself a considerable challenge due to the illposedness of the problem in this frequency regime. In this work, strategies for detecting underground objects as well as objects inside metallic containers using ELF/VLF signals are considered.
3
CHAPTER II
Theory of Electromagnetic Scattering
II. 1 Equivalent Dipole Formulation
In general, the electric field that is scattered by a conductive object can be formulated by the volume integral equation (VIE) given by (II. 1).
Here, the quantity J is the equivalent current density in the material and is given by J = aE, where E is the total and unknown electric held in the material. Once the electric held in the material (and corresponding currents) is determined, the scattered magnetic held can be computed using (II.2).
The problem of determining the conductivity prohle is in general nonlinear since it is a product of the unknown electric held and unknown conductivity. However, instead of inverting the conductivity distribution of an object, inferring the equivalent current prohle from magnetic held measurements is a good proxy for recovering the geometry of the object in question. This approach has the added beneht of linearizing the image reconstruction problem.
With this equivalent current reconstruction method in mind, formulation of the LF electromagnetic inversion problem takes inspiration from the source localization techniques of electroencephalography (EEG) and other related medical imaging techniques [14]. Specifically, electric currents inside the human body typically operate at frequencies that less than 1kHz [5]. The purpose of EEGlike techniques is to infer the magnitude and location of these currents in a noninvasive manner. Thus, the
(II1)
(11*2)
4
EEG imaging methodology is mathematically identical to the general nearheld image reconstruction problem with ELF/VLF signals.
As suggested by [10], neural activity in the brain can be modeled as a series of equivalent current dipoles with dipole moment di = c^ei that is located at rdiP;i. Each current dipole contributes to the voltage potential measured at the electrode location r. As such, the full description of each equivalent dipole requires six parameters in R3. By superposition, the potential measured at an electrode due to these equivalent sources can be expressed as (II. 3) [14].
m(rj) = X^(rlâ€™rdip,i>di) = X^(rlâ€™rdip,i) d* (H.3)
i i
Assuming N observation locations and p equivalent sources, (II.3) can be expressed as a linear system relating the electrode measurements to the equivalent source dipole moments (II.4) or similarly in its condensed form as (II.5).
m(rx)
g(ri,rdiP;i) â€¢ â€¢ â€¢ g(n, rdiP;P)
m(rN)
gij^N> Iâ€™dipq) ' ' ' gipTSii rdiP)P)
(IF 4)
M = GD (II.5)
This formulation (II.5) takes the same form as the integral equation formulation of the scattered electromagnetic fields described in (II. 2). To formulate the electromagnetic integral equation in the same manner as in EEG imaging, equations (II.2) and (II. 1) can be discretized using pulse basis functions (II.6) and point matching test functions (II.7) on the discretized imaging region.
A?) = 5^a*P*(r), pi(*)
i
1, ? e Si
0, o.w.
(II6)
5
fm = 6{\rrg\) (II.7)
The dipole moment in M3 is described by a weighted vector cq over the arbitrarily defined cell Si. The point matching testing functions assures sampling at the center of the ith cell [23]. With this choice, the VIE of Esc (II. 1) and Bsc (II.2) transforms into a linear system resembling (II.4).
Bsc{ ri)
GB(r 1,r1') â– â– â€¢ Gfl(ri,rp')
Ji
Bsc{ rN)
Gb{fNTiO â€¢ â€¢ â€¢ Gb{rN,rp')
(IL8)
Bsc = GbJ (II.9)
The problem of equivalent current reconstruction using magnetic held measurements is now recast in an identical mathematical form to the EEG problem (II.5). However, since generally N > p, this linear system is overdetermined and inversion of (II.9) must be solved in a leastsquares sense. Accurate reconstruction of the forward problem from this measurement data (ESC,BSC) requires regularization and a priori constraints which are discussed in Chapter III.
II.2 Scattering in FreeSpace
Consider the problem of electromagnetic scattering by a perfectly electrically conducting (PEC) sphere of radius a in free space. This scatterer is illuminated by a TEpolarized plane wave with wavenumber k propagating in the +y direction as shown in Fig. II. 1.
6
Oo^o)
Figure II. 1: Scattering of PEC sphere of radius a in freespace illuminated by a TEpolarized plane wave. In the longwavelength limit, the fields scattered by this can be approximated by the fields of a â€”^oriented ideal magnetic dipole.
In the longwavelength limit ka 1, the spatial variation of the plane wave is small and can be approximated with a uniform held of ZP(y) = H0e~jkyz ~ H0z = Ho (cos Or â€” sin0 0). The boundary conditions on the PEC sphere require the normal component of the fields to be continuous at the interface and a surface current Js to be induced tangentially to the surface [23].
n â€¢ H = 0
h x H = Js
(II. 10)
Application of the normal boundary condition (II. 10), in a spherical coordinate frame, requires (II. 11).
r â– Hsc\r=a =cosOH0 (11.11)
Application of the tangential boundary condition at the interface requires the induced surface current Js to be â€” 0directed. Computation of the scattering due these surface currents can be simplified by considering the fields radiated by an equivalent â€” idirected magnetic dipole located at the center of the PEC sphere [15].
7
(11.12)
> 1 ml
Hsc(r) = â€”[3(m â€¢ r)r â€” m] = â€” â€”[2cos0 r + sind 9}
Airr6 47t rA
The combination of these results give exact expressions of the magnetic dipole moment
in the nearheld as (1113).
m=â€”2ir a3#' (11.13)
These analytical expressions for Hsc were compared to method of moments surface integral equation (MoMSIE) solutions obtained in the commercial electromagnetic solver FEKO along a cutline 2 meters above the PEC sphere (c = 2m) as in Fig. II.2.
y(m)
Figure II.2: Comparison of analytical scattered magnetic fields due to ideal dipole against MoMSIE solutions of sphere scattering and equivalent magnetic dipole.
These analytical results for Hsc agree considerably with the MoM solutions of this scattering problem. This simple approach of representing magnetic scattering by conductive objects with equivalent magnetic dipoles is a valid approximation in the nearheld for the unknown electromagnetic scatterers in the imaging problem. Since a magnetic dipole is identical to a loop of electrical current, this approach is
8
identical to searching for loops of current (as opposed to electrical dipoles) that are induced on an object. This is a convenient abstraction that helps alleviate some of the computational burden associated with image reconstruction.
II.3 Scattering by a PEC Sphere Buried in a Homogeneous HalfSpace
Consider the electromagnetic scattering of a PEC sphere illuminated by an incident plane wave as shown in Fig. II.3. The PEC sphere is embedded within a homogeneous halfspace with parameters (Â£,p.,
Figure II.3: Configuration of PEC sphere of radius a embedded in homogeneous halfspace (Region 2) with conductivity a illuminated by an incident plane wave.
For simplification of the analysis, the plane wave is assumed to be a TEpolarized plane wave propagating with normal incidence to the interface such that ko = â€”koz, where k0 is the wavenumber in free space. It is worth nothing that signals from existing VLF navy communications transmitters are typically TM polarized and at grazing incidence [12], Nevertheless, the TE normal incidence case is a conservative estimate for analytical calculations. Due to an impedance mismatch between the halfspace (Region 2) and freespace (Region 1), reflections occur at the interface
0 = o).
9
Vo
V2
â€” Â« 377Q
to
a
Vo
eo  J
JUJ/lO
a
r a
for â€” > e0
OJ
(II14)
The magnetic field in Region 1 is given by,
Hly{z) = H0 e?koZ  RteH0 e~jkÂ°z (11.15)
where the quantity RTE is the Fresnel reflection coefficient [9]. A fraction of the magnetic held gets transmitted into Region 2 with Fresnel transmission coefficient
rpTE
H2y{z) = TTEH0elk2Z (11.16)
For normal incidence, the Fresnel coefficients for the magnetic held are related by 1 + Rte = Tte. Specifically, the Fresnel rehection and transmission coefficients are given by (IF 17) and (11.18) [9].
Râ„¢
tko â€” to k2 tko + tok2
(H17)
rpTM
2 tk0
tko + tok2
(II18)
The signal transmitted into Region 2 suffers considerable losses initially due to the impedance mismatch between the two media at the interface and then from attenuation as it propagates deeper into the halfspace. The attenuation losses result from the complex nature of the wavenumber k2 = Re^) + jlm(k2).
Scattering by the buried PEC sphere now will only be due to the helds in Region 2 (11.16). However, the approaches detailed in Section II.2 are not entirely valid for this problem since the helds can potentially vary across the length of the scatterer
10
due to the background medium. Satisfying the boundary conditions (11.10) results in a nonuniform induced Js over the â€”0direction. The attenuation of H2y(z) over the length of the scatterer results in a stronger response in the illuminated region of the PEC sphere (c > â€”cl) compared to the shadow region (c < â€”cl) by the simulation results depicted in Fig. II.4.
Surface current [mA/m]
â–
6.600
6.140
5.680
5.220
4.760
4.300
3.840
3.380
2.920
2.460
2.000
Figure II.4: MoMSIE solution of Js of PEC sphere buried in homogeneous halfspace (
Since Js is not uniform across the â€”0direction, a single horizontal magnetic dipole (HMD) (â€”y) is not sufficient to describe the radiated scattered fields. Instead, this must be approximated with a series of HMD with decreasing magnitude over the length of the sphere. For the purposes of this analysis, this will be roughly approximated by a single HMD with dipole moment m.
11
The fields radiated by a HMD dipole with current moment IÂ£ are given by [9] as
(11.19).
g jkr
47t r
(11.19)
E(r) = V x
e~jkr
For the case considered, a is taken to be â€”y. Including the effect of a finite conductivity halfspace requires reformulating the Greenâ€™s functions via a spectral expansion [9]. Specifically, the Sommerfeld identity can be applied for a spectral representation of the point source and is given by (11.20).
This represents a spherical wave expanded as a series of cylindrical waves propagating in the pdirection (p = xcos(f> + ysliuf>) and a planewave propagating in the z
of the Sommerfeld approach is that each cylindrical wave can be reflected from the airground interface via the traditional Fresnel coefficients. The net field is thus a continuous superposition of each reflected and transmitted wave. This allows the total field to be written as a single integral expression that can then be approximated with analytical or numerical calculations.
(11.20)
direction for all wavenumbers in kspace (k2 = k2 + k2) [9]. The primary advantage
12
(% P o)
1
z = d
2
Figure II.5: HMD (ct = â€”y) serves as a rough approximation for the electromagnetic scattering by a PEC sphere (centered at c = 0).
H 2,z(z)
Hu(:)
g^cos0 j dkpkpH[1] (kpp)
_ m POO ^
â€”^cos (f) / dkpkpH[1] (kpp)
^ e~jk2,z\z\ _ ftTEejk2,t(z+2d)
rpTE jk0,z(z+d)+jk2 td 12,1 e
(11.21)
R
â– TE
2,1
Pokâ€˜2Z â€” pko^ Pokâ€˜2^ + pkoâ€ž
T,
TE
2,1
2p,0h
2,z
(11.22)
Pok^'Z + pkO'Z
In (11.21), (â€”) is chosen for the primary held for waves propagating in Region 2 at c > 0.
This methodology can be extended to account for multilayered planar media by taking a similar approach as in (11.21), that is by considering reflections and transmissions in each layer of interest and the source region [9]. For simplicity, only the homogeneous halfspace problem is considered in this work.
13
Computation of the Sommerfeld Integrals in (11.21) are in general analytically intractable and thus require numerical techniques. Typical asymptotic methods detailed in [9] are highly problem dependent and potentially pose complications for scattering problems with little to no a priori information about the scatterer. In this work, this problem of evaluating the Sommerfeld Integrals is simplified by utilizing the commercial electromagnetic solver FEKO to quickly compute both the forward scattering problem as well as to generate a numerical Greenâ€™s Function response (NGF) for an equivalent magnetic magnetic dipole buried in a conductive halfspace.
II.4 PEC Object Inside a Thin Conductive Shell
Of particular interest is the detection of metallic objects that are located inside a metal shield (such as a shipping container). Consider a modification to the electromagnetic scattering problem posed in Sec. II.2. The PEC sphere of radius a is now shielded by an arbitrary thin, conductive shell with thickness A = d and is illuminated by a uniform magnetic held H = Hoy. The problem is now considerably more complex as the shield attenuates the signal that illuminates the PEC object, and the scattered fields have to once again traverse the shield which results in additional attenuation. The magnetic held that is scattered by the PEC object can still be approximated by that of a magnetic dipole, however, the free space response must be scaled by 1/SE#. The quantity SEh represents the shielding ehectiveness of the outer shield. The reason for squaring SEH is that the incident signal is hrst shielded before illuminating the PEC sphere, and is again shielded when exiting the shield. For the case of a spherical thinshell shield, there exists a simple analytical solution for SEH and is given by (11.23).
14
feo^o)
Figure II.6: A spherical shield with thickness A = d and outer radius a illuminated by a uniform magnetic held as in [20]. A single PEC sphere is located inside the shield.
H 11
SEh = â€”â€” = cosh kmd + ~(kma + ) sinh kmd
11 s O rCmCL
(11.23)
The quantity km corresponds to the wavenumber inside the shell and is given by km = (1 + j)/8, where 8 is skin depth from equation (1.1). The fields scattered by the PEC sphere inside the shield are thus given by (11.23), where H0 are the fields in the absence of the shield while Hs correspond to the shielded fields [20]. Equation (11.23) is an oversimplified approximation that is used to estimate the strength of the scattered fields. However, the problem of magnetic scattering from a PEC object that is inside a metal shield requires a numerical approach.
15
CHAPTER III
INVERSE PROBLEM FORMALISM III. 1 Necessity for Regularization
Consider the general linear electromagnetic scattering problem posed in (II.9). Generally, this problem, like the EEG inverse problem [14], is highly illposed. Let (II.9) be recast as (III. 1), with A being the forward operator, b G B being the problem data set, and x G X being the unknown model parameters. These parameters are generally complex valued and so can be assumed to be vectors in B G Cm and Xe Câ„¢, where Câ„¢ = {a + bj \ (a, b) G Eâ„¢}.
Ax = b (III.l)
This problem (III.l) can be solved exactly if it is wellposed, i.e. satisfies both Hadamard conditions [31]:
1. For each b G B, a unique solution exists.
2. The solution is stable under perturbation, i.e. A1 is defined for all B and is continuous.
The electromagnetic scattering problem typically fails both of these conditions and is Hadamard illposed.
With the equivalent dipole technique, infinitely many dipoles are required to accurately model the essentially continuous imaging region. Discretizing the imaging region still would require many more equivalent dipoles than there are measurement locations (n > m) for a sufficient model. Eq. (III.l) will only be consistent if b G 7Z(A), since by definition the range 71(A) = {Ax \ x G Câ„¢} [7].This is equivalent to requiring rank(A) = rank(A6), where (A\b) is the augmented matrix. By
RoucheCapelli theorem [27], this solution is unique if and only if the rank(A) = n.
16
This cannot occur for the underdetermined system since rank(A) < n and A is full rank only if rank(A) = min{m,n} [7]. Thus, if the underdetermined system is consistent, the solution is not unique and violates the uniqueness Hadamard condition. Additionally, the computational costs of building a sufficient and accurate model for each inverse problem is large due to the number of equivalent dipoles required. It is instead realistic to require n < m. However, the conditions for uniqueness and existence of the model parameters that satisfy (III. 1) are not guaranteed since A is not guaranteed to be full rank.
For all possible data b, the solution is also unstable  it is highly sensitive to noise or other small perturbations in the data [14], This instability can be quantified by the condition number k(A) [21],
k(A) = A A 11,for nonsingular A
k(A) â€˜7â€â€™"r
(HI.2)
(7n
where omax, amin are the maximum and minimum singular values of A [24].
Consider the problem posed in [19] of a collection of N ideal dipoles along a line of length L that is located in freespace.
17
Figure III. 1: Figure from [19]. A collection of N ideal dipoles in freespace uniformly spaced (with spacing Ac) along a line of length L are used to compute the condition number of the forward operator Gb to determines its stability.
The dipole array is surrounded by Nr receivers that are bound to the surface of a sphere of radius L/2. The Greenâ€™s function of an ideal magnetic dipole G#(r, r') is computed for this configuration. The condition number of this forward operator is computed using (III.2). The results are shown below in Fig. III.2.
Condition Number of FreeSpace Gg
L/Az
Figure III.2: Condition number of Gb is computed for a collection of dipoles N = L/Ac bound on line of length L. For L = 1 m, the normalized A correspond to waves of / = 300 GHz, 300 MHz, and 300 kHz, respectively.
18
It is evident that the solution becomes quickly unstable at low frequencies. Solutions at higher frequency will remain relatively stable while the minimum resolution Az > Rff, with Rff = 1.22A being the Rayleigh limit [If]. The higher frequency considered in Fig. III.2 is likely to become unstable at N zz 1700 dipoles. For small N = L/Az, this condition number is relatively small due to the sufficient spacing between dipole pairs. The fields of any dipole is distinguishable in this limit.
Therefore, it is clear that (III. 1) fails at least one, if not both, of the Hadamard conditions for wellposedness in this frequency regime. While the inverse problem is illposed and difficult to solve, an optimal estimator x exists that both minimizes the error between the data set b and the estimate Ax and also represents the most likely â€œimageâ€ or distribution of equivalent dipoles expected in the forward problem. A typical metric for measuring error between the data and observations is the unconstrained ordinary least squares minimization problem is given by (III.3).
min 6 â€” Ax\\\ (III.3)
X
However, the least squares solution doesnâ€™t enforce certain desirable features in the image such as sparse or minimal model parameters. This is accomplished through the addition of a regularization term to (III. 3) that acts as a mathematical constraint to the minimization problem.
min 6 â€” Ax\\\ + A:rE (HI4)
Equation (HI.4) represents the general regularized minimization problem, where A is the regularization parameter. The second term in (HI.4) represents the general pnorm and is given by (III.5).
IXI
(III.5)
19
The choice of A corresponds to a certain family of solutions to the optimization problem. Two common choices of pnorm in the literature are the p = 2 Tikhonov regularization [14] andp = 1 LASSO regularization (or regression) [30]. These pnorm regularized problems are guaranteed to have unique and globally optimal solutions in the problem space since these functions are convex for p > 1 [7].
III.2 L2 Tikhonov Regularization Theory
Consider the L2norm minimization problem (Tikhonov Regularization problem).
min 6 â€” Ax\\\ + Arc (III.6)
X
It is a fundamental property of convex optimization problems that the locally optimal solution is also globally optimal [7]. The globally optimal solution to (III.6) has an analytical form since the objective function Â£ = 6 â€” Ax2 + Ax2 is differentiable for all x in the problem space.
Â£ =  b â€” Ax \\l + Ax2 = xt(AtA + A I)x â€” bT Ax + bTb (HI7)
Terms independent of x can be neglected in the minimization problem. To aid in the computation of its derivative, the quadratic term (ATA + XI) is a symmetric matrix.
dÂ£
dx
2xt(AtA + XI)  bTA = 0
x={A1A + XI)~lA1b
(III.8)
This regularization problem attempts to minimize all predicted components of the estimator x uniformly.
III.3 Li Regularization Theory
Consider the Linorm minimization problem (III.9).
min 6 â€” Ax\\\ + X\\x
i
(III.9) 20
While the objective function is convex, it is not differentiable due to the Linorm penalty. This is because the absolute value function
C = xtAtAx â€” bTAx + A ^ \xi\ (III. 10)
i
For simplicity, (III.9) can be recast as its dual [7].
min 6 â€” Ax\\\ subject to y \xi\ â€” ^
i
The LASSO (least absolute shrinkage and selection operator) [30] takes an identical form to the primal problem (III. 11).
N
(aj) = arg mini ^(y* a ^PjXijf ^ i= 1 j
subject to EIAI Â£ A
(III.12)
For all regularization parameters, the optimal estimate of a is taken to be the mean of the data y [30]. Without loss of generality, this constant term is often taken to be a = 0.
A closedform solution for the orthonormal design case (XTX = I) can be obtained by considering subdomain where the absolute value function f(x) = x is differentiable for x ^ 0.
dx
f(x)
1, x < 0
(III.13)
1, x > 0
Solutions to the modified objective can be found for each subdomain B+ and O, where B+ is defined for f3j > 0 and O for f3j < 0. Let f3Â° = yTX be the ordinary least squares solution for the orthonormal design case.
21
(III. 14)
ft = i [2/3?  A], ft > 0
k = \[ 2i3â€ + A],ft<0
This simply becomes (III. 15) for the union of these subdomains O = O UO+, where sgn(x) is the sign function. The constrained domain of this problem requires the regularization parameter to be \/3j\ > A/2.
ft = ygn(ftÂ°)[2ftÂ«A] (111,15)
III.3.1 Tibshiraniâ€™s Algorithm for Solving LASSO
Tibshirani had proposed an algorithm for a quadratic program to solve the LASSO optimization problem. For a fixed regularization parameter, (III. 12) is expressed as a least squares problem with m = 2n inequality constraints. These inequality constraints are introduced sequentially until a feasible solution is obtained that satisfies the KuhnTucker conditions [30]. The inequality constraint JT \(3j\ < A is expressed as 8f/3 < A, where 8i for i = 1,2,..., 2â„¢ are the ntuples of the form (Â±1,..., Â±1). Let Ge be the matrix whose rows contain 5i for i e E, where E = {i : 5f (3 = A} is the equality set and S = {i : 5f (3 < A} is the slack set, and 1 be the vector of ones of the same length as Gg.
for n
Gf
1 On 1 1 1
s2 1 1
8s 1 1
1 On l 1 1 Iâ€”1 1 Iâ€”1 1
(III.16)
Tibshiraniâ€™s algorithm is initialized with E = i0 with 8i0 = sgn(/5Â°). If the original inequality constraint JT \(3j\ < A is satisfied, the algorithm completes. Otherwise, the violated constraint is removed from the slack set and added to the equality set.
22
Initialize ftÂ°
Initialize E = {i0\ 8foft = X] do while ^
add {i} Â» E with = sgn(/?) ft = min$(/?) s. t. GEp < A1
Figure III.3: Pseudocode describing Tibshiraniâ€™s algorithm for solving LASSO. This minimizes the objective of the primal problem with the constraints in the equality set E.
Tibshiraniâ€™s algorithm is guaranteed to converge in at most 2â€ steps. This is because the final iterate is the solution to the original problem which does satisfy the KuhnTucker conditions. With the application of standard quadratic programming techniques, this can reduce the number of steps required for convergence to 2n + 1 steps [30]. However, it is clear that this is impractical for large n due to the length of the tuples required to model the constraints. Instead, different techniques for solving the quadratic program are desirable.
III.4 Optimal Choice of Regularization Parameter
Each choice of regularization parameter A > 0 corresponds to a different family of solutions of the regularized optimization problem. The optimal choice of A is critical for solving for the true forward model parameters corresponding to the imaging problem in question. Choosing this optimal A* is equivalent to minimizing the error between the model and the data. However, the difficulty in finding this A* arises from the lack of explicit knowledge of the forward problem. Since these true forward model parameters are typically unknown in the imaging problem, statistical and imaging techniques must be used to provide an estimate of the true model error. Two techniques commonly applied in the literature are the LCurve method and CrossValidation (CV) method.
23
III.4.1 LCurve Method
The LCurve method, named for the characteristic shape of the loglog plot as in Fig. III.4, attempts to minimize the tradeoff between the generalized regularization bias q = 11L(x) \  and the residual p = \\bâ€”Ax\(2 [14], Qualitatively, the optimal choice of A* is given by the corner of the curve. If A is too large, the inversion solution is dominated by regularization errors such as oversmoothing. Otherwise, if A is too small, the solution will be dominated by the modeling errors and measurement noise
[18].
log B. , G J III:
3 11 data inv"2
Figure III.4: Characteristic loglog plot of LCurve Method for Tikhonov Regularization. The corner or closest point to the origin represents the optimal parameter A*.
There are two methods of locating this corner [18].
1. Finding the corner is equivalent to finding the closest point on the curve to the origin. This metric differs based on the choice of regularization. For Tikhonov regularization, this metric is cl(A) = p + A2q.
2. The corner is also the point of maximum curvature /i(A). This is independent of
the choice of regularization since n{X) is a purely geometric quantity. Assuming
24
the Lcurve is a smooth curve and that p = log 116 â€” Ax\[2 and i) = log 11112 for the Tikhonov Regularization, then the curvature is given by (III. 17), where differentiation is taken with respect to A.
k(X)
p'if â€” p"if
m2 + wÂ»3/2
(III.17)
III.4.2 CrossValidation Method
The crossvalidation method is a statistical technique of estimating the mean squared error (MSE) of the solution without explicit knowledge of the true solution. Only kfold CV is considered for these problems [22], This estimate converges to the true MSE since these samples share the statistical properties of the data set. The kfold CV involves dividing the data set b into fcequally sized partitions. One of these partitions is kept as a holdout (validation) set and the remaining k â€” 1 are used as training sets, where k is typically taken to be 5 or 10.
Figure III.5: Equipartition of data set for 5fold CV. This provides an estimate of the MSE of the solution vector x that can be used to optimize the regularization parameter A.
For a fixed regularization parameter, the regularized optimization problem is then solved using the k â€” 1 training sets. The optimal solution x is projected onto the holdout set by b' = A!x, where A' is the forward operator corresponding to the holdout
25
set. This procedure is repeated for the remaining sets until each of the partitions is used as a holdout set [22],
cv(A) = ipÂ«A&f> = iÂ£ MSBÂ® (III.1S)
i= 1 i= 1
III.5 Imaging Resolution Limits
The smallest details able to be resolved by the imaging system is highly dependent on the source and the regularization techniques considered. In imaging systems that operate in the farheld (ka 1), this best resolution is determined by the Rayleigh limit, given by (III. 19) [11].
Rff = 1.22^ (III.19)
Spatial features that are larger than are resolvable by the system.
Far Field Resolution Limit (Ravleiqh Limit)
0 5 10 15 20 25 30 35 40 45 50
Frequency [kHz]
Figure III.6: Minimum Spatial Resolution for signal propagating in various homogeneous finiteconductivity grounds. This assumes no attenuation of the signal.
It is evident from Fig. III.6 that the objects of interest (which are typically sizes
on the order of lm) are not resolvable at all by an imaging system operating in the far
held assuming a noiseless and lossless signal. The information of these small spatial
26
features are contained in the evanescent waves which are neglected in the farheld approximation [11], However, this nearheld resolution limit will differ greatly for each regularization technique applied. Harid [19] gives the minimum resolution limits for an imaging system solvable by an ordinary least squares inversion.
Figure III. 7: From [19]. This suggests for an imaging region of size L m, the minimum resolution obtainable for an ordinary least squares inversion (when the forward operator A is full rank and invertible) is ~ 30 â€” 40L mm.
The resolution limits attainable by a Tikhonov or L\ regularized solution are not examined but are assumed to be equal, if not better, than those for an unregularized solution as in [19]. In the problems where a full volumetric inversion is considered, the resolution of the uniform grid is greater than the limit posed for the ordinary least squares solution and are assumed not to be limited by this resolution limit.
27
CHAPTER IV
ALGORITHM DESCRIPTION
IV. 1 HighLevel Overview of Algorithm
The algorithms described attempt to solve the inverse problem for the lowfrequency nearheld electromagnetic scattering problem. It relies on the creation of synthetic forward data (Bdata) and the NGF (Gb) in the commercial electromagnetic solver FEKO [3]. This inversion prcess is agnostic to the choice of electromagnetic solver. This solver was chosen due to the availability and advantages of the moments method (MoM) electromagnetic solver over other commercial solvers. The MoM solver has simple and efficient treatments of ideal dipole sources and only requires meshing of the surface elements to solve the scattering problem. The inverse problem is solved either analytically for the Tikhonov regularized problem, which is optimized by the LCurve method, or solved with the aid of the CVX MATLAB toolbox [13] for the Linorm regularized problem, optimized by the CV method. Two techniques for the generation of the imaging grid are presented: a simple structured volumetric generation of the imaging grid and an iterative and less computationally intensive generation of the imaging grid termed â€œadaptive mesh refinementâ€ (AMR).
28
Figure IV. 1: Flowchart describing highlevel design of software for both volumetric mesh and adaptive mesh refinement procedures.
IV.2 Generation of Synthetic Data
The forward problem is fully modelled in FEKO. Only three cases are considered for the forward problem: (i) the scattering problem posed in Sec. II.2, (ii) scattering from a PEC sphere buried in a conductive halfspace as in Sec. II.3, and (iii) the scattering from PEC objects located inside in a metallic cubic shield. Figure IV.2 shows an example forward simulation setup of two PEC cubes located inside a lm x lm x lm conductive cubic shell of 2mm thickness and conductivity a = 10' S/m. Figure IV.4 shows the setup of a PEC sphere in freespace as well as the case of a conductive halfspace (halfspace not visible in simulation CAD file).
29
Figure IV.2: Forward Problem Geometry for Metallic Box with two PEC Cubes. The metallic shell has a = 10' S/m and thickness A = 2mm. Receivers are placed on planes parallel to each face offset by 5mm. An ideal magnetic dipole source is placed over each face that operates in the LF regime.
The process of generating synthetic forward data in the FEKO solver is general and so the inversion process can potentially be extended to any class of electromagnetic inverse scattering problems. The geometry of the forward problem, which includes both the unknown scatterers and the known problem boundaries, is modeled with the CAD tools provided in the FEKO solver. This solver then solves the electromagnetic scattering problem using the Surface Integral Equation MoM (SIEMoM) approach [3]. The surfaces of the modeled geometry are approximated by a mesh of triangular elements by FEKO.
30
Figure IV.3: Surfaces are meshed with many triangular elements. The basis functions used to solve the SIE are mapped to the edges of these elements and used to compute surface current densities.
While increasing the number of these mesh elements improves the accuracy of the model, it does so at a increased computational cost. Material properties such as conductivity and thickness are applied to the applicable surfaces of the model. Last, the nearheld data of the scattered magnetic fields are requested at Nr arbitrary points placed outside the problem space.
Generation of synthetic forward data is handled from FEKO
Figure IV.4: Generation of Forward Problem of PEC Sphere located in freespace handled in FEKO. The unknown PEC scatterers are placed in this model with the known problem boundaries. The scattered nearheld data are requested at Nr points.
31
IV.3 Computation of Numerical Greenâ€™s Function (NGF)
Closed form expressions for the Greenâ€™s function response of the halfspace and
analytically. Instead, computation of these are handled numerically with the aid of the LUA scripting tools available in FEKO. Assume a single ideal magnetic dipole with unit dipole moment is placed at r and its magnetic held is observed at n. Then, the Greenâ€™s function response Gb for that collection of source and observation locations (r',fi) are computed from (IV.1) by iterating through all possible unit magnetic dipoles along a and solving for the respective matrix entry.
This computation is computed for all observation locations Nr and source locations N, resulting in Gb, a matrix of size 3Nr x 3N. The use of numerical Greenâ€™s function is a useful method to include the effects of the background shield (halfspace or shell) without having to explicitly model the background in the inversion algorithms.
IV.4 Inversion Algorithm
Three different types of inversion methods are attempted to solve for the unknown model parameters MmG Each attempts to impose certain favorable features in the model space to accurately invert the problem in question.
IV.4.1 Tikhonov Regularization
Tikhonov regularizedinversion, or L2 Regularization, attempts to minimize the 2norm of the solution vector Mâ€œ uniformly (IV.2).
metallic box problems are unavailable due to the infeasiblility of solving for these
(IV. 1)
inv 112
+ X\\Mmv\\22
(IV.2)
32
X
Solution to (IV. 2) is straightforward since the minimization problem can be solved analytically as described in Section III.2. The optimal regularization parameter is found through the LCurve method as detailed in Sec III.4.1. The Regularization Tools package for MATLAB is utilized [17] to aid in computation of A*. This treats the discrete Lcurve as a differentiable curve to compute the point of maximum curvature Â«(**)â€¢
IV.4.2 Li Regularization
The Li regularized inversion problem attempts to enforce sparsity in the solution vector Mâ„¢ (IV.3).
min Bdata  GBMinv\\22 + AHMHIi (IV.3)
X
Since the 1norm penalty is not differentiable, it cannot be computed analytical solution to (IV.3). This is further described in Sec III.3. Instead, the CVX toolbox for MATLAB is implemented in the algorithm to handle both constrained and unconstrained forms of (IV.3) to solve Mmv [13]. The regularization parameter is optimized through CV method as in Sec III.4.2.
IV.4.3 Iterative Least Squares Solution
This method greatly differs from the previous two in that it applies no traditional regularization to the minimization problem. Instead, it takes Ndip dipoles, where Ndip < N, and iterates them through all possible locations in the imaging grid. There will be M = ) total combinations. The ordinary least squares solution is
taken for each combination. The combination that minimizes the mean squared error (MSE) is taken as the solution (IV.4).
MSE(j)
ji)data __ q
B,jMJnv
II Bd>
lata 112
j
1,2 ,...,M
(IV.4)
33
IV. 5 Generation of Imaging Grid
The imaging grid is a discretization of the continuous forward model space or problem space. Each cell in the imaging grid represents a voxel (volumetric pixel) in the final image. Increasing the number of cells N improves the resolution of the image as well as increasing the computational costs of constructing the NGF and solving the inverse problem. The volumetric reconstruction is the simplest method of constructing this imaging grid. This involves the generation of a uniform rectangular grid inside the bounds of the physical model space. If N = NxNyNz, then the resolution of this grid is given by (IV.5).
Ax Ay A z
\%ub %lb\
Nx 1 \yub ~ yib\ Ny~ 1
\Zub %lb\
Nz 1
(IV.5)
It is necessary to require these upper and lower bounds rÂ«, < r < ruy be offset from the physical problem boundaries. Intersecting grid cells with the physical boundaries often results in singular values of the Greenâ€™s function response. However, the use of this volumetric grid requires large computational costs to achieve an adequate resolution. Additional modeling errors to the inversion could arise from the inappropriately chosen grid cells. To counteract these adverse consequences from the volumetric inversion, an efficient technique called â€œadaptive mesh refinementâ€ (AMR) is implemented.
IV.5.1 Adaptive Mesh Refinement Algorithm
The motivation of the AMR algorithm is to use the outputs of a very coarse mesh (typically N = 8) to aid in local refinements of the grid. The algorithm of a structured AMR is as follows:
34
1. Initialize with N = 8 grid cells with each grid cell located at the center of each octant of the full model space.
2. Solve the regularized inversion problem for this configuration.
3. Refine the imaging grid with a refined subgrid with N = 8 at the location of the strongest response. The size of new grid cells is half of the size of the original grid cell.
4. Repeat until the size of these new grid cells becomes negligible compared to the model space. In experiments, this threshold was typically chosen to be 15% of the dimensions of the global model space.
In this way, each iteration effectively adds 7 grid cells to the model since the original grid cell is removed during the refinement process. The algorithm for a single iteration is demonstrated in Fig. IV.5.
y C r2
Figure IV.5: 2D Crosssection of initial coarse grid (left). Distance between grid centers dr = \rf â€” r/f\, for j % i, is taken to be the size of the grid cell. At the next iteration (right), the cell with the strongest response is removed and replaced with a subgrid with 8 cells of size dr/2.
The selection criterion typically considered for the grid refinement is the greatest  Minv .
35
CHAPTER V
RESULTS
Summarized below are the inversion results from the application of the algorithm and techniques described to the three cases considered. It is worth noting that these results do not demonstrate accurate reconstructions of the forward problem, thus indicating the necessity for further optimization of the inversion process.
V.l Inversion results for FreeSpace Inversion
This inversion algorithm was validated with the scattering by a PEC sphere of radius a = 20cm located in freespace at (x,y,z) = (0,0,â€”0.5m). The scatterer is illuminated by a plane wave source with El = â€”1 V/m x propagating with normal incidence to the xyplane (similar to the incident fields described in Sec. II. 3). A factor of ej27r^, with / = 25kHz, is assumed. The nearheld data is sampled along points on the plane z = +10 cm on x G [â€”10m, 10m] and y G [â€”10m, 10m], The forward data of the scattering problem in freespace limited to the subdomain x G [â€”2m, 2m] and y G [â€”2m, 2m] is shown in Fig. V.l. The characteristic shape of this held data indicates the strongest response in the center (x = 0,y = 0). This suggests the presence of the PEC scatterer near the origin, though the depth and geometry cannot be conclusively determined from direct examination of this forward data.
36
HyField Magnitude [uA/m]: FreeSpace Illuminated by Plane Wave
X position [m]
Near field bof structure (Frequency = 25 kHz; Z position = 100 mm)  empty_forward_model
Figure V.l: Synthetic forward data of nearheld scattered data computed by FEKO (domain limited to (x, y) G [â€”2m, 2m]).
Tikhonov Inversion w/ Regularization: A = 5.4163e08
^Inverted Current Density Mmv Â° Forward Scatterer * Maximum M'nv___________
0 ^ Max Minv found at x = 2.5, y = 2.5, z = 0.1
ylm] 3 3 ' x [m]
Figure V.2: Volumetric Inversion with Tikhonov Regularization. Optimal Inverted
Current Density Mtraâ€™ found at (x = â€”2.5m, y = â€”2.5m, c = â€”0.1m)
37
A volumetric inversion regularized with the Tikhonov penalty is performed shown in Fig. V.2. A uniform grid of 8 x 8 x 8 cells is constructed on the domain (x,y) G [â€”2.5m, â€”2.5m] and c G [â€”2m, 0.1m]. The largest magnitude inverted current density was selected as the optimal value. However, the optimal response was located at the top of the discretized imaging grid at (x = â€”2.5m, y = â€”2.5m, c = â€”0.1m). Since the LCurve method had selected such a small value for the regularization parameter A* = 5.4163e â€” 08, the results are likely dominated by modeling errors, since all in the imaging grid are all relatively large and uniform. Choosing a larger regularization parameter biases the inversion results of this data set toward the top of the imaging region (closest to the receiver plane at z = 10cm). The bias towards cells closest to the measurement plane in the overregularized inversion is due to the minimization of Tikhonov penalty. Since the magnitude of the equivalent current density is minimized, the sources must be inverted closest to the measurement plane to produce the model the forward held data accurately. The Tikhonov penalty is not suitable for the configuration of this volumetric inversion.
L1 Inversion w/ Regularization: A = 1e05
0 ^ Max Minv found atx=1,y = 1,z = 2
^ Inverted Current Density Mmv â€¢ Maximum Minv
Â° Forward Scatterer
1.5
Figure V.3: Volumetric Inversion with Li Regularization. Optimal Inverted Current
Density Mtraâ€™ found at (x = â€” lm, y = â€” lm, c = â€”2m)
38
Constructed similarly to Fig. V.2, a volumetric inversion regularized with L\ penalty was performed as depicted in Fig. V.3. The CV method had selected the smallest A* tested, since the estimated crossvalidation curve was monotonically increasing over the tested regularization parameters for this data set. This behavior is likely due to either that the L1 penalty is not applicable to this problem or that a different choice of regularization parameters is necessary for this configuration. It could be the case that the optimal regularization parameters lied outside the range of A tested (105 < A < 1). This choice of small regularization parameter had not enforced sparsity in the solutions since the current densities were relatively uniform across all grid cells and this is likely to not occur for smaller regularization parameters. Unlike the Tikhonov inversion, the L\ regularized problem inverted the optimal current density to a location at the bottom of the imaging grid (x = â€” lm, y = â€” lm, c = â€”2m). It is inconclusive whether this is biased towards the bottom of the region or that this technique is unsuitable for a volumetric inversion of the freespace problem.
1.4 1.61.82
Iterative Least Squares Inversion
Optimal Minv found at x = 1, y = 1, z = 2
Inverted Volumetric Current Density Forward Scatterer Location
y[m]
0.5
x [m]
Figure V.4: Volumetric Inversion with Iterative Least Squares. Optimal Inverted
Current Density Mtraâ€™ found at (x = â€” lm, y = â€” lm, c = â€”2m)
39
The volumetric inversion was performed using the iterative least squares method as shown in Fig. V.4. This solution is highly unsuited for this data set. The normalized MSE (IV.4) was computed across all grid cells in the uniform 8x8x8 cell grid and each were found to be identical. Due to identical MSE values, the minimum index chosen corresponded to the first index in the vector of Mmv.
Tikhonov Inversion w/ Regularization: A = 1.216e08
â€”^Inverted Volumetric Current Densities 0 Forward Scatterer Location â€¢ Strongest Inversion Result_____
Max Minv found at x = 1.875, y = 1.875, z = 0.15937
Figure V.5: AMR Inversion with Tikhonov Regularization after 5 Iterations. Optimal Inverted Current Density Mtraâ€™ found at (;r = â€” 1.875m, y = â€”1.875m, c = 0.15937m)
The Tikhonov regularized inversion was performed with the AMR algorithm as depicted in Fig. V.5. After 5 iterations of the mesh refinement algorithm, the optimal inverted current density was located at (x = â€” 1.875m, y = â€”1.875m, c = â€”0.15937m). This was not tested for more than 5 iterations since the optimal values quickly diverged away from the location of the scatterer in the forward problem. Further iterations appeared to bias towards the top of the imaging domain. This bias is likely due to the same reasons that the volumetric inversion had posed inadequate for this problem: locations closest to measurement plane are those that simultaneously
minimize the Tikhonov penalty and the residual between the data and the prediction.
40
This could also indicate that the current metric applied to the AMR algorithm is not an adequate method of selecting the refinement subdomain.
L1 Inversion w/ Regularization: A = 1e05
Max Minv found at x = 1.875, y = 1.875, z = 0.15937
â€”> Inverted Volumetric Current Densities Â° Forward Scatterer Location â€¢ Strongest Inversion Result
Figure V.6: AMR Inversion with Li Regularization after 5 Iterations. Optimal Inverted Current Density Mtraâ€™ found at (x = â€” 1.875m, y = â€” 1.875m, c = 0.15937m)
The Li regularized inversion was performed with the AMR algorithm as depicted in Fig. V.6. After 5 iterations, the results indicated that the same optimal choice of inverted current density was selected as in the Tikhonov regularized problem. The convergence of these results is likely due to the selection criteria for subdomain refinement, which likely results in a bias towards the top of the imaging domain. Another likely explanation is that due to the unsuitability of the L1 regularization to this problem (due to the absence of a local minimum in the CV curve), the inversion is no better than doing an ordinary least squares inversion i.e. minimizing the unregularized minimization problem (III.3) due to the absence of local minima in the domain.
A proposed test to evaluate the suitability of the regularization technique to a
data set would be to examine the behavior of the CV curve as A â€”> 0. If the CV
41
curve decreases asymptotically in this limit, then the regularization method is no better than an ordinary least squares solution due to the absence of local minima in the curve.
These inconclusive inversion results for the freespace problem are likely attributable to the choice of the measurement plane. The lack of variation in the zdirection could directly be responsible for these issues. The added diversity in measurements could be crucial in resolving depth information, though no results are available for this scattering problem to verify the effects of different choices of measurement geometry.
V.2 Inversion results for Homogeneous HalfSpace Inversion
The forward problem was constructed identically to the freespace problem in Sec.
V.l. A ground plane simulated with the exact Sommerfeld integrals was included at z = 0 with a = 10 mS/m. The forward data shares the characteristic shape as in the freespace problem, though the measured fields are twice as large due to reflection from the interface between the freespace region and the finiteconductivity halfspace at z = 0. Similarly, the location of the PEC scatterer can be determined to be near the origin {x = 0,y = 0) due to the location of the maximum in the magnetic held data. As in the freespace problem, it is not clear what the geometry and depth of the PEC scatterer are from direct examination of the forward data.
42
HyField Magnitude [uA/m]: Conductive Ground Illuminated by Plane Wave
X position [m]
Near field bof structure (Frequency = 25 kHz; Z position = 100 mm)  ground_forward_model
Figure V.7: Synthetic forward data of nearheld scattered data computed by FEKO (domain limited to (x,y) E [â€”2m, 2m]).
Â£ 1
2;
2
Tikhonov Inversion w/ Regularization: A = 4.3545e08
Max Minv found at x = 2, y = 2, z = 0.91429
'' V
' \N
Vlfc
V.
\,
'' \ \ '' 
' â– H.,
â– v.
v.
o
W.., 'V "W.
\s
N\,
\
 \ ,
\ ,
^ Inverted Current Density M Maximum Minv
Forward Scatterer
y [m]
2 2
x [m]
Figure V.8: Volumetric Inversion with Tikhonov Regularization. Optimal Inverted
Current Density \Mmv\ found at (x = â€”2m, y = â€”2m, c = â€”0.91429m)
43
The volumetric inversion regularized with the Tikhonov penalty is performed as depicted in Fig. V.8. The optimal current density was located at (x = â€” 2m, y = â€” 2m, c = â€”0.91429m). Unlike the results of the freespace problem, this was biased towards the bottom of the imaging grid with a similar regularization parameter A*. The behavior of biasing towards the boundaries of the imaging grid could be highly problem dependent. Testing for overfitting in this problem biases the location of the optimal inversion results to the bottom of the imaging grid. The reason for the disparity in this localization behavior between the freespace and conductive halfspace problem are not clear.
L1 Inversion w/ Regularization: A = 1e05
^Inverted Current Density Mmv Â° Forward Scatterer * Maximum M'nv________________
0 Max Minv found at x = 2, y = 2, z = 2
Figure V.9: Volumetric Inversion with L\ Regularization. Optimal Inverted Current Density \Mmv\ found at (x = 2m, y = 2m, c = â€”2m)
The volumetric inversion regularized with the Li penalty is performed as depicted in Fig. V.9. The optimal inverted current density computed is identical to the one found in the Tikhonov regularized problem. The results are again inconclusive as to whether this primarily due to modeling errors of the forward operator Gb or due to the choice of cell locations in the uniform grid. However, the CV curve for optimizing
44
this problem shared the same behavior as the freespace problem. This again suggests that the L\ penalty is simply not suitable to the problem of scattering of a PEC sphere in freespace or the scattering of a PEC sphere buried in a homogeneous halfspace.
Iterative Least Squares Inversion
Inverted Volumetric Current Density Forward Scatterer Location
0.2 0.4 0.6 0.8 
I 1"
N 1.2
1.4 
Optimal M  found at x = 2, y = 2, z = 2
y N
x N
Figure V.10: Volumetric Inversion with Iterative Least Squares. Optimal Inverted Current Density Mtraâ€™ found at (V = â€”2m, y = â€”2m, c = â€”2m)
The volumetric inversion was performed using the iterative least squares method as shown in Fig. V.10. As in the freespace problem, the computed MSE of the ordinary least squares solution for each grid cell was identical, giving the first element in the vector Mmv. Attempts to ht this problem with individual dipoles are likely unsuitable for either (i) the cells of uniform grid are not correctly chosen as they do not coincide with the interior of the PEC sphere or (ii) the forward operator for a single dipole is illconditioned and a least squares solution serves as a poor ht for the data set.
45
Tikhonov Inversion w/ Regularization: A = 1.0194e08
> Inverted Volumetric Current Densities Forward Scatterer Location Strongest Inversion Result
Max Minv found at x = 1.875, y = 1.875, z = 0.87187
0.5 ^ A\\ 
\\Wy
E, "1 " \%vA' \
N W\
1.5 2 : , . 4 \ ' v, \
2 0 y[mI 4 A x [m]
Figure V.ll: AMR Inversion with Tikhonov Regularization after 5 Iterations. Optimal Inverted Current Density Mtraâ€™ found at (A = â€”1.875m, y = â€”1.875m, c = 0.87187m)
The Tikhonov regularized inversion was performed with the AMR algorithm as depicted in Fig. V.ll. After 5 iterations, the optimal inverted current density was located at (x = â€” 1.875m, y = â€”1.875m, c = â€”0.87187m). While the localization in the iryplane was poor, the inversion results demonstrate significant improvement from the volumetric problem in the localization of the depth of the PEC sphere. The error in the depth from the bottom of the sphere is Az = 0.17187m. Improvement in the localization rryplane could be attained through the inclusion of physical constraints in the inversion problem.
46
L1 Inversion w/ Regularization: A = 1e05
 Max Minv found at x = 1.875, y = 
0.5 ^ ft
 v
A \
N \ 'v'' \
1.5 A
2 : , . 4
> Inverted Volumetric Current Densities Forward Scatterer Location Strongest Inversion Result
y [m]
4 4
x [m]
Figure V.12: AMR Inversion with L\ Regularization after 5 Iterations. Optimal Inverted Current Density Mtraâ€™ found at {x = â€” 1.875m, y = â€” 1.875m, c = 0.63438m)
The Li regularized inversion was performed with the AMR algorithm as depicted in Fig. V.12. fter 5 iterations, the optimal inverted current density was located at (x = â€” 1.875m, y = â€” 1.875m, c = â€”0.63438m). The inversion localized the depth of the optimal current density to the interior of the PEC sphere, though localization in the iryplane is poor as in the results of the Tikhonovregularized inversion.
The improvement in localization for these two cases could be thought of as the inclusion of an additional constraint into the general optimization problem (III.4). Effectively, the AMR algorithm serves as an additional regularization penalty by constraining the domain of the problem to coarse grid cells. The technique has potential but requires optimization to handle issues such as divergence in the refinement process and determining the optimal amount of iterations required to accurately invert the forward model parameters.
V.3 Inversion results for Conductive Cubic Shell Inversion
The final problem considered is imaging through a thin, highly conductive cu
47
bic shell of size lm3 (thickness A = 2mm) with a = 10' S/m. Two smaller PEC cubes with side length Â£ = 0.2m are placed in the bottom of the shell at (xi, yi) = (0.5m, 0.5m) and (^2,1/2) = (0.5m, 0.1m). This is so that the Cube 1 (center cube) only makes contact with the bottom face while the Cube 2 makes contact with the bottom and y = 0 faces. The configuration is illuminated by a planar array of nine dipoles 5mm above the top face as in Fig. V.13.
Figure V.13: Configuration of Forward Problem with planar array of dipoles on top face and nearheld requests below bottom face.
The nearheld scattered data is requested 1mm below the bottom face. The outlines of the PEC cubes in contact with this face clearly appear in the measured held data as in Fig. V.14. Such information is valuable in establishing physical constraints to the inversion problem for objects near or in contact with the metallic shell. This was used to constrain the imaging domain so that Mmv is nonzero only for grid cells corresponding to this outline: this subdomain is given by (xi,yi) E [0.4m, 0.6m] x [0.4m, 0.6m] and (:r2, t/2) G [0.4m, 0.6m] x [2mm, 0.202m],
48
HzField Magnitude [dBA/m]
0.0
0.8
1.6
2.4
3.2
4.0
4.8
5.6
6.4
7.2
8.0
X
Figure V.14: Synthetic forward data of nearheld scattered data computed by FEKO.
49
Tikhonov Inversion w/ Regularization: A = 0.12026
Inverted Current Density Mmv
â€¢ Maximum Minv
Forward Scatterer
1 . Max Minv found at x = 0.002, y = 0.002, z = 0.14429
y[ml o'"o x [m]
Figure V.15: Volumetric Inversion with Tikhonov Regularization. Optimal Inverted Current Density Mtraâ€™ found at (x = 2mm, y = 2mm, c = 0.14429m).
The volumetric inversion regularized with the Tikhonov penalty is performed as depicted in Fig. V.15. The optimal current density was located at (x = 2mm, y = 2mm, c = 0.14429m). The inversion performed as well as the L\ regularized inversion in Fig. V.12. The inversion localizes the optimal Mmv to a cell at the edge of the imaging, demonstrating poor localization in the rryplane. However, the ^component of the optimal Mmv was localized to a height corresponding to the two PEC cubes of 0.002m < c < 0.202m. Improvements to this xy localization are attained through the inclusion of the aforementioned physical constraints to the inversion problem depicted in Fig. V.19.
50
L1 Inversion w/ Regularization: A = 1e05
Figure V.16: Volumetric Inversion with L\ Regularization. Optimal Inverted Current Density found at {x = 2mm, y = 2mm, c = 2mm).
The volumetric inversion regularized with the L\ penalty is performed as depicted in Fig. V.16. The optimal current density was located at (x = 2mm, y = 2mm,; = 2mm). This localized to the edge of the imaging domain. The CV curve that optimizes this inversion problem share the same behavior as the previous volumetric Li regularized results. The absence of local minima suggest that this regularization method is not suitable to this problem.
51
Tikhonov Inversion w/ Regularization: A = 4.3518e05
Figure V.17: AMR Inversion with Tikhonov Regularization. Optimal Inverted Current Density \Mmv\ found at (x = 34.0643mm, y = 34.0643mm, c = 34.0643mm).
The Tikhonov regularized inversion was performed with the AMR algorithm as depicted in Fig. V.17. After 5 iterations, the optimal inverted current density was located at (x = 34.0643mm, y = 34.0643mm, c = 34.0643mm). These results again demonstrate a poor localization of the optimal Mmv towards the bottom of the imaging domain.
52
L1 Inversion w/ Regularization: A = 2.6102e05
â€”Â»Inverted Volumetric Current Densities
â€¢ Strongest Inversion Result
Forward Scatterer Location
Max Minv found at x = 0.034063, y = 0.034063, z = 0.034063
Figure V.18: AMR Inversion with L\ Regularization. Optimal Inverted Current Density \Mmv\ found at (x = 34.0643mm, y = 34.0643mm, z = 34.0643mm).
The Li regularized inversion was performed with the AMR algorithm as depicted in Fig. V.18. After 5 iterations, the optimal inverted current density was located at (x = 34.0643mm, y = 34.0643mm, c = 34.0643mm). These results again demonstrate a poor localization of the optimal Mmv towards the bottom of the imaging domain.
53
Tikhonov Inversion w/ Regularization: A = 0.0075354
0.9 0.8 0.7 0.6 Â£0.5 N 0.4 0.3 0.2 0.1
Figure V.19: Constrained Volumetric Inversion with Tikhonov Regularization. Optimal Inverted Current Density Mtraâ€™ found at (x = 0.4m, y = 3mm, c = 0.145m). This point lies along one of the edges of Cube 2.
The volumetric inversion regularized with the Tikhonov penalty is performed as depicted in Fig. V.19. The volumetric mesh is constructed with the inclusion of the physical constraints determined from Fig. V.14 such that Mmv is nonzero for grid cells located on the subdomains corresponding to the contours in that plot. These cells correspond to the two subdomains (Vi,yi) G [0.4, 0.6] and x\ G [0.4, 0.6] U y2 G [0.002,0.202], The optimal Mt1w was localized to a cell on the edge of one of the PEC cubes at (x = 0.4m, y = 3nun,; = 0.145m). This technique appears to be promising in that (i) the constrained domain effectively improves the resolution of the volumetric grid and (ii) the Mmv was localized to the surface the PEC cubes rather than a freespace region inside the imaging domain. A similar number of cells was used to construct each subdomain which achieved a better resolution compared to the unconstrained problem. The uniform imaging grid was generated for each of the two subdomains with 6x6x8 cells.
The magnitude of the inverted current density Mmv across all grid cells is depicted
â€”Â»Inverted Volumetric Current Densities
â€¢ Strongest Inversion Result
Forward Scatterer Location
Max Minv found at x = 0.4, y = 0.003, z = 0.145
54
in Fig. V.20. The uniform features imposed by the Tikhonov penalty are apparent
in this plot of \Mmv\. Maxima of correspond to cells on the edge of the PEC
cubes.
Current Density Magnitude
16111
^,14 
â€¢n
E
>,12 
a>
T3
0 100 200 300 400 500 600
Grid Index
Figure V.20: Current Density magnitude of constrained volumetric inversion with Tikhonov Regularization. Grid Index of maximum corresponds to optimal inverted current density \Mmv\.
L1 Inversion w/ Regularization: A = 0.00074989
â€”^Inverted Volumetric Current Densities â€¢ Strongest Inversion Result IM Forward Scatterer Location______
Figure V.21: Constrained Volumetric Inversion with Li Regularization. Optimal Inverted Current Density Mtraâ€™ found at (x = 0.4m, y = 0.4m, c = 3mm). This point lies along one of the edges of Cube 1.
55
The volumetric inversion regularized with the L\ penalty is performed as depicted in Fig. V.21. The volumetric mesh is constructed identically to the mesh of Fig. V.19 with the inclusion of the physical constraints determined from Fig. V.14. The optimal Mmv was localized to a cell on the edge of one of the PEC cubes at (x = 0.4m, y = 0.4m, c = 3mm), inverting to one of the corners of the central PEC cube. Again, inverting an optimal Mmv to the surface of the PEC cube does suggest that this method of constraining the imaging domain is crucial to imaging the interior of the shell.
The magnitude of the inverted current density Mmv across all grid cells is depicted in Fig. V.22. The sparse features imposed by the L1 penalty are apparent in this plot of \Mmv\. Only a single strong response is evident in this plot while all other responses are exceedingly small. The L\ penalty succeeded in inverting this data set. Unlike the previous applications of this penalty, the CV curve did possess a local minimum A* = 7.4989 x 104 in the domain considered (105 < A < 1), suggesting that these constraints also improve the behavior of the regularization process.
Current Density Magnitude
80 i i i
â€” 70 "
co
E
>60 
CD
"O
Â£ 50 c
CD
CD
^ 40 >>
W
Â§j 30 Q
2 20 Zj
O'11111
0 100 200 300 400 500 600
Grid Index
Figure V.22: Current Density magnitude of constrained volumetric inversion with Li Regularization. Grid Index of maximum corresponds to optimal inverted current density \Mmv\. The effects of the L\ penalty are clearly demonstrated in these sparse inversion results.
56
CHAPTER VI
FUTURE WORK
The inversion algorithm must be further refined to in order to produce meaningful and accurate reconstructions of the considered classes of nearheld imaging problems. Direct applications of these classical inversion techniques are insufficient in solving the illposed nearheld imaging problem. The results shown in Figs. V.19V.22 indicate that these techniques beneht greatly from the inclusion of physical constraints. Further analysis of the forward problems are necessary to hnd essential physical constraints to optimally solving the inverse problem.
The adaptive mesh rehnement (AMR) algorithm does indeed reduce the computational costs, since it necessitates the computation of fewer dipoles compared to the volumetric inversion, but its results in aiding the inversion process are inconclusive. It is difficult to determine whether the inverted solution diverges from the forward solution after a few iterations due to the selection criteria in the grid rehnement process or due to the general illposedness of the inverse problem. A proposed solution involves turning this structured mesh rehnement process into a stochastic one. The imaging grid would be initialized with a hnite quantity of grid cells to be determined from a statistical distribution in the global domain. Subsequent iterations rehne the selected subdomain using this statistical distribution until the convergence criteria are met. This stochastic AMR procedure is repeated many times and the acquired results are averaged to compute the inverted current density. While this remains untested, the motivation is to apply MonteCarlolike techniques to aid in solving this deterministic inverse problem.
Last, the inaccuracies in the reconstruction could be attributed to the choice of regularization in the inverse problem. Simply, the L\ penalty or the Tikhonov penalty may not be suitable to these classes of problems. One alteration would be an elastic
57
net penalty [33] which combines the effects of these two penalty terms. However, the problem of requiring the correct choice of the penalty terms can likely be avoided entirely using machine learning methods to solve the inverse problem [1],
While all these potential paths to improving the accuracy of the reconstruction of the nearfield electromagnetic imaging problem are promising candidates, these are not guaranteed to accurately and meaningfully solve these specific inverse problems. Either more sophisticated techniques are necessary or that these classes of electromagnetic forward problems must be treated individually.
58
CHAPTER VII
CONCLUSION
Presented here is an inversion algorithm for solving the general lowfrequency nearheld electromagnetic inverse problem. The electromagnetic scattering problem is treated simply as a linear system with a forward operator determined by the Greenâ€™s function response of the problem Gb and the measured held data in regions external to the imaging region Bsc. The forward problem and the numerical Greenâ€™s function responses are computed with the aid of the commercial electromagnetic method of moments solver FEKO. The unknown scatterers in this region are identically represented by equivalent volumetric current densities Jeq with unknown weights. Due to the illposedness of the inverse problem in the LF regime, it is exceedingly difficult to solve for these unknown model parameters. This can be alleviated with the inclusion of regularization and physical constraints into the inverse problem. Various regularized inversion techniques were tested with this algorithm on the three cases considered: a PEC sphere in freespace, a PEC sphere in a conductive halfspace, and PEC cubes shielded by a highlyconductive thin shell. Direct applications of the regularized inverse problem on these three considered cases do not produce conclusive or meaningful results. However, the introduction of physical constraints into the inverse problem considerably improved the reconstruction since strong responses were inverted to the surfaces of the internal PEC cubes, though it does not conclusively reconstruct the geometry of the scatterers. Further research must be done in order to optimize this inversion algorithm.
59
REFERENCES
[1] Jonas Adler and Oza Oktem. Solving illposed inverse problems using iterative deep neural networks. Inverse Problems, 33(12): 1 24, 2017.
[2] Ahmed I. AlShammaâ€™a, Andrew Shaw, and Saher Saman. Propagation of electromagnetic waves at MHz frequencies through seawater. IEEE Transactions on Antennas and Propagation, 52(ll):28432849, 2004.
[3] Altair Engineering Inc. Altair FEKO, 2018.
[4] A. P. Annan. ELECTROMAGNETIC PRINCIPLES OF GROUND PENETRATING RADAR. In Harry M. Jol, editor, Ground Penetrating Radar Theory and Applications, pages 78. Elsevier B.V., 2009.
[5] Sylvain Baillet and Line Garnero. A Bayesian Approach to Introducing AnatomoFunctional Priors in the EEG/MEG Inverse Problem. IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, 44(5):374385, 1997.
[6] Glenn D. Boreman. Modulation Transfer Function in Optical and ElectroOptical Systems. SPIE Press, 2001.
[7] Stephen Boyd and Lieven Vandenberghe. Convex Optimization. Cambridge University Press, 1st edition, 2004.
[8] Leonardo Carrer and Lorenzo Bruzzone. Solving for ambiguities in radar geophysical exploration of planetary bodies by mimicking bats echolocation. Nature Communications, 8(1), 2017.
[9] Weng Cho Chew. Waves and Fields in Inhomogeneous Media. IEEE Press, 1995.
[10] Jan C. de Munck, Bob W. van Dijk, and Henk Spekreijse. Mathematical Dipoles are Adequate to Describe Realistic Generators of Human Brain Activity. IEEE Transactions on Biomedical Engineering, 35(11):960 966, 1988.
[11] Gregoire Derveaux, George Papanicolaou, and Chrysoula Tsogka. Resolution and denoising in nearheld imaging. Inverse Problems, 22(4): 1437â€”1456, 2006.
[12] K. L. Graf, N. G. Lehtinen, M. Spasojevic, M. B. Cohen, R. A. Marshall, and U. S. Inan. Analysis of experimentally validated transionospheric attenuation estimates of VLF signals. Journal of Geophysical Research: Space Physics, 118(5):27082720, 2013.
[13] Michael Grant and Stephen Boyd. CVX: Matlab software for disciplined convex programming, 2013.
[14] Roberta Grech, Tracey Cassar, Joseph Muscat, Kenneth P Camilleri, Simon G Fabri, Michalis Zervakis, Petros Xanthopoulos, Vangelis Sakkalis, and Bart Vanrumste. Review on solving the inverse problem in EEG source analysis. Journal of N euro Engineering and Rehabilitation, 5(25): 1â€”33, 2008.
60
[15] David J. Griffiths. Introduction to Electrodynamics. PrenticeHall, Inc., 3rd edition, 1999.
[16] Robert E. Grimm, Barry Berdanier, Robert Warden, James Harrer, Raymond Demara, James Pfeiffer, and Richard Blohm. A timedomain electromagnetic sounder for detection and characterization of groundwater on Mars. Planetary and Space Science, 57(11): 1268â€”1281, 2009.
[17] P.C. Hansen. Regularization Tools Version 4.0 for Matlab 7.3, 2007.
[18] Per Christian Hansen and Dianne Prost Oâ€™Leary. The Use of the LCurve in the Regularization of Discrete IllPosed Problems. SIAM J. Sci. Comput., 14(6): 1487â€”1504, 1993.
[19] Vijay Harid. Resolution limits of nearfield electromagnetic imaging. In 2018 International Applied Computational Electromagnetics Society Symposium in Denver, ACESDenver 2018, number 2, pages 12. Applied Computational Electromagnetics Society (ACES), 2018.
[20] Vijay Harid, Mark Golkowski, Stephen Gedney, Ronald Rorrer, Poorya Hosseini, Morris Cohen, Nathan Opalinski, and Sarah Patch. Modeling NearField Magnetic Shielding by Conductive Enclosures at Very Low Frequencies. 2019.
[21] Nicholas J. Higham. Accuracy and Stability of Numerical Algorithms. Society for Industrial and Applied Mathematics, 2nd edition, 2002.
[22] Gareth James, Daniela Witten, Trevor Hastie, and Robert Tibshirani. An Introduction to Statistical Learning: with Applications in R. 2013.
[23] Jian Ming Jin. Theory and Computation of Electromagnetic Fields. John Wiley & Sons, Incorporated, 2nd edition, 2015.
[24] Daniel Lichtblau and Eric W. Weisstein. Condition Number.
[25] Yuri Alvarez Lopez and Jose Angel Martinez Lorenzo. Compressed sensing techniques applied to ultrasonic imaging of cargo containers. Sensors (Switzerland), 17(1): 1â€”13, 2017.
[26] Lorenzo Lo Monte, Danilo Erricolo, Francesco Soldovieri, and Michael C. Wicks. Radio frequency tomography for tunnel detection. IEEE Transactions on Geoscience and Remote Sensing, 48(3 PART 1): 1128â€”1137, 2010.
[27] Igor R. Shafarevich and Alexy O. Remizov. Linear Algebra and Geometry. Springer, 2013.
[28] Deraid G Smith and Harry M Job Ground penetrating radar: antenna frequencies and maximum probable depths of penetration in Quaternary sediments. Journal of Applied Geophysics, 33(13):93â€”100, 1995.
61
[29] W. M. Telford, W. F. King, and A. Becker. VLF MAPPING OF GEOLOGICAL STRUCTURE. Technical report, Energy, Mines and Resources Canada, 1977.
[30] Robert Tibshirani. Regression Shrinkage and Selection via the Lasso. Journal of the Royal Statistical Society, 58(1):267â€”288, 1996.
[31] A.N. Tikhonov, A. Goncharsky, V.V. Stepanov, and A.G. Yagola. Numerical Methods for the Solution of IIIPosed Problems. Springer Netherlands, 1 edition, 1995.
[32] Stuart M. Wentworth. Fundamentals of Electromagnetics with Engineering Applications. Wiley, 2006.
[33] Hui Zou and Trevor Hastie. Regularization and Variable Selection via the Elastic Net. Journal of the Royal Statistical Society. Series B (Statistical Methodology), 67(2):301320, 2005.
62

Full Text 
PAGE 1
LOWFREQUENCYIMAGINGALGORITHMFORINVERTINGARBITRARY ELECTROMAGNETICSYSTEMS by DALIBORJ.TODOROVSKI B.S.,UniversityofCentralFlorida,2016 Athesissubmittedtothe FacultyoftheGraduateSchoolofthe UniversityofColoradoinpartialfulllment oftherequirementsforthedegreeof MasterofScience ElectricalEngineering 2019
PAGE 2
ThisthesisfortheMasterofSciencedegreeby DaliborJ.Todorovski hasbeenapprovedforthe ElectricalEngineeringProgram by VijayHarid,Chair VijayHarid,Advisor MarkGolkowski StephenGedney Date:May9,2019
PAGE 3
Todorovski,DaliborJ.M.S.,ElectricalEngineering LowFrequencyImagingAlgorithmforInvertingArbitraryElectromagneticSystems ThesisdirectedbyAssistantProfessorVijayHarid ABSTRACT MethodsandsoftwareweredevelopedinMATLABforthepurposesofinverting arbitrarylowfrequencyelectromagneticscatteringproblemsandlocatingtheequivalentcurrentdensitiespresentintheforwardproblem.Whiletheclassesofproblems consideredforthismethodarethosewhereforwardscatterersareconcealedwithin ahighlyconductivehalfspaceorametallicshell,themethodscanbegeneralizedto otherproblemsduetotherobustnessoftheelectromagneticsolverused.Synthetic electromagneticelddatawascreatedfortheconsideredcaseswiththeuseofthe commercialelectromagneticsolverFEKO.Similarly,theGreen'sfunctionresponse wasconstructednumericallyinthissolverusingtheknownphysicalboundariesof theproblemanditsmaterialproperties.Withthesecomputedquantities,alinear relationshipwasformedbetweentheforwardelectromagneticdataandthenumerical Green'sfunctionNGFthatcouldbeinvertedwiththeaidofvariousregularization techniques.ResultsarepresentedforinversionsusingbothTikhonovregularization andL 1 regularizationonafullvolumetricreconstructionoftheimagingspaceaswell asamorecomputationallyecientiterativereconstruction. Theformandcontentofthisabstractareapproved.Irecommenditspublication. Approved:VijayHarid iii
PAGE 4
ACKNOWLEDGMENTS ThankyoutoalloftheprofessorsattheUniversityofColoradoDenver.Youhave alltaughtmeinvaluableknowledgeandskillsasanengineerandresearcherthatI willcarryforwardwithme.AndIwouldliketoextendmygratitudetoDr.Vijay Harid.Yourguidanceandencouragementoverthesepasttwoyearshavebeencritical inhelpingmendmypathforward. iv
PAGE 5
TABLEOFCONTENTS CHAPTER I.MOTIVATION.................................1 II.THEORYOFELECTROMAGNETICSCATTERING...........4 II.1EquivalentDipoleFormulation.....................4 II.2ScatteringinFreeSpace.........................6 II.3ScatteringbyaPECSphereBuriedinaHomogeneousHalfSpace.9 II.4PECObjectInsideaThinConductiveShell..............14 III.INVERSEPROBLEMFORMALISM.....................16 III.1NecessityforRegularization.......................16 III.2 L 2 TikhonovRegularizationTheory..................20 III.3 L 1 RegularizationTheory........................20 III.3.1Tibshirani'sAlgorithmforSolvingLASSO...........22 III.4OptimalChoiceofRegularizationParameter.............23 III.4.1LCurveMethod.........................24 III.4.2CrossValidationMethod....................25 III.5ImagingResolutionLimits.......................26 IV.ALGORITHMDESCRIPTION........................28 IV.1HighLevelOverviewofAlgorithm...................28 IV.2GenerationofSyntheticData......................29 IV.3ComputationofNumericalGreen'sFunctionNGF.........32 IV.4InversionAlgorithm...........................32 IV.4.1TikhonovRegularization.....................32 IV.4.2 L 1 Regularization.........................33 IV.4.3IterativeLeastSquaresSolution.................33 IV.5GenerationofImagingGrid.......................34 IV.5.1AdaptiveMeshRenementAlgorithm.............34 v
PAGE 6
V.RESULTS....................................36 V.1InversionresultsforFreeSpaceInversion...............36 V.2InversionresultsforHomogeneousHalfSpaceInversion.......42 V.3InversionresultsforConductiveCubicShellInversion........47 VI.FUTUREWORK................................57 VII.CONCLUSION.................................59 REFERENCES...................................60 vi
PAGE 7
FIGURES II.1ScatteringofPECsphereofradius a infreespaceilluminatedbyaTEpolarizedplanewave.Inthelongwavelengthlimit,theeldsscattered bythiscanbeapproximatedbytheeldsofa )]TJ/F15 11.9552 Tf 10.007 0 Td [(^ z orientedidealmagnetic dipole......................................7 II.2Comparisonofanalyticalscatteredmagneticeldsduetoidealdipole againstMoMSIEsolutionsofspherescatteringandequivalentmagnetic dipole......................................8 II.3CongurationofPECsphereofradius a embeddedinhomogeneoushalfspaceRegion 2 withconductivity illuminatedbyanincidentplane wave......................................9 II.4MoMSIEsolutionof ~ J s ofPECsphereburiedinhomogeneoushalfspace =1S/milluminatedbyTEpolarizedplanewave f =1kHz.A larger ~ J s isinducedatthetopofsphererelativetothebottomdueto attenuationoverthelengthofthesphereinsidethehalfspace.......11 II.5HMD^ = )]TJ/F15 11.9552 Tf 10.09 0 Td [(^ y servesasaroughapproximationfortheelectromagnetic scatteringbyaPECspherecenteredat z =0...............13 II.6Asphericalshieldwiththickness= d andouterradius a illuminatedby auniformmagneticeldasin[20].AsinglePECsphereislocatedinside theshield....................................15 III.1Figurefrom[19].Acollectionof N idealdipolesinfreespaceuniformly spacedwithspacing z alongalineoflength L areusedtocomputethe conditionnumberoftheforwardoperator G B todeterminesitsstability.18 III.2Conditionnumberof G B iscomputedforacollectionofdipoles N = L= z boundonlineoflength L .For L =1m,thenormalized correspondto wavesof f =300GHz,300MHz,and300kHz,respectively........18 vii
PAGE 8
III.3PseudocodedescribingTibshirani'salgorithmforsolvingLASSO.This minimizestheobjectiveoftheprimalproblemwiththeconstraintsinthe equalityset E .................................23 III.4CharacteristicloglogplotofLCurveMethodforTikhonovRegularization.Thecornerorclosestpointtotheoriginrepresentstheoptimal parameter ..................................24 III.5Equipartitionofdatasetfor5foldCV.Thisprovidesanestimateofthe MSEofthesolutionvector x thatcanbeusedtooptimizetheregularizationparameter ................................25 III.6MinimumSpatialResolutionforsignalpropagatinginvarioushomogeneousniteconductivitygrounds.Thisassumesnoattenuationofthe signal......................................26 III.7From[19].Thissuggestsforanimagingregionofsize L m,theminimum resolutionobtainableforanordinaryleastsquaresinversionwhenthe forwardoperator A isfullrankandinvertibleis 30 )]TJ/F15 11.9552 Tf 11.955 0 Td [(40 L mm.....27 IV.1Flowchartdescribinghighleveldesignofsoftwareforbothvolumetric meshandadaptivemeshrenementprocedures...............29 IV.2ForwardProblemGeometryforMetallicBoxwithtwoPECCubes.The metallicshellhas =10 7 S/mandthickness=2mm.Receiversare placedonplanesparalleltoeachfaceosetby5mm.Anidealmagnetic dipolesourceisplacedovereachfacethatoperatesintheLFregime...30 IV.3Surfacesaremeshedwithmanytriangularelements.Thebasisfunctions usedtosolvetheSIEaremappedtotheedgesoftheseelementsandused tocomputesurfacecurrentdensities.....................31 viii
PAGE 9
IV.4GenerationofForwardProblemofPECSpherelocatedinfreespacehandledinFEKO.TheunknownPECscatterersareplacedinthismodel withtheknownproblemboundaries.Thescatterednearelddataare requestedat N r points.............................31 IV.52DCrosssectionofinitialcoarsegridleft.Distancebetweengridcenters dr = j r C i )]TJ/F19 11.9552 Tf 12.241 0 Td [(r C j j ; for j 6 = i ,istakentobethesizeofthegridcell.Atthe nextiterationright,thecellwiththestrongestresponseisremovedand replacedwithasubgridwith8cellsofsize dr= 2...............35 V.1SyntheticforwarddataofneareldscattereddatacomputedbyFEKO domainlimitedto x;y 2 [ )]TJ/F15 11.9552 Tf 9.298 0 Td [(2m ; 2m]....................37 V.2VolumetricInversionwithTikhonovRegularization.OptimalInverted CurrentDensity j M inv j foundat x = )]TJ/F15 11.9552 Tf 9.299 0 Td [(2 : 5m ;y = )]TJ/F15 11.9552 Tf 9.299 0 Td [(2 : 5m ;z = )]TJ/F15 11.9552 Tf 9.298 0 Td [(0 : 1m..37 V.3VolumetricInversionwith L 1 Regularization.OptimalInvertedCurrent Density j M inv j foundat x = )]TJ/F15 11.9552 Tf 9.298 0 Td [(1m ;y = )]TJ/F15 11.9552 Tf 9.299 0 Td [(1m ;z = )]TJ/F15 11.9552 Tf 9.298 0 Td [(2m.........38 V.4VolumetricInversionwithIterativeLeastSquares.OptimalInvertedCurrentDensity j M inv j foundat x = )]TJ/F15 11.9552 Tf 9.299 0 Td [(1m ;y = )]TJ/F15 11.9552 Tf 9.298 0 Td [(1m ;z = )]TJ/F15 11.9552 Tf 9.299 0 Td [(2m.......39 V.5AMRInversionwithTikhonovRegularizationafter5Iterations.OptimalInvertedCurrentDensity j M inv j foundat x = )]TJ/F15 11.9552 Tf 9.298 0 Td [(1 : 875m ;y = )]TJ/F15 11.9552 Tf 9.299 0 Td [(1 : 875m ;z = )]TJ/F15 11.9552 Tf 9.298 0 Td [(0 : 15937m..........................40 V.6AMRInversionwith L 1 Regularizationafter5Iterations.OptimalInvertedCurrentDensity j M inv j foundat x = )]TJ/F15 11.9552 Tf 9.299 0 Td [(1 : 875m ;y = )]TJ/F15 11.9552 Tf 9.299 0 Td [(1 : 875m ;z = )]TJ/F15 11.9552 Tf 9.299 0 Td [(0 : 15937m..................................41 V.7SyntheticforwarddataofneareldscattereddatacomputedbyFEKO domainlimitedto x;y 2 [ )]TJ/F15 11.9552 Tf 9.298 0 Td [(2m ; 2m]....................43 V.8VolumetricInversionwithTikhonovRegularization.OptimalInverted CurrentDensity j M inv j foundat x = )]TJ/F15 11.9552 Tf 9.299 0 Td [(2m ;y = )]TJ/F15 11.9552 Tf 9.298 0 Td [(2m ;z = )]TJ/F15 11.9552 Tf 9.299 0 Td [(0 : 91429m.43 ix
PAGE 10
V.9VolumetricInversionwith L 1 Regularization.OptimalInvertedCurrent Density j M inv j foundat x =2m ;y =2m ;z = )]TJ/F15 11.9552 Tf 9.298 0 Td [(2m...........44 V.10VolumetricInversionwithIterativeLeastSquares.OptimalInvertedCurrentDensity j M inv j foundat x = )]TJ/F15 11.9552 Tf 9.299 0 Td [(2m ;y = )]TJ/F15 11.9552 Tf 9.298 0 Td [(2m ;z = )]TJ/F15 11.9552 Tf 9.299 0 Td [(2m.......45 V.11AMRInversionwithTikhonovRegularizationafter5Iterations.OptimalInvertedCurrentDensity j M inv j foundat x = )]TJ/F15 11.9552 Tf 9.299 0 Td [(1 : 875m ;y = )]TJ/F15 11.9552 Tf 9.299 0 Td [(1 : 875m ;z = )]TJ/F15 11.9552 Tf 9.299 0 Td [(0 : 87187m..........................46 V.12AMRInversionwith L 1 Regularizationafter5Iterations.OptimalInvertedCurrentDensity j M inv j foundat x = )]TJ/F15 11.9552 Tf 9.299 0 Td [(1 : 875m ;y = )]TJ/F15 11.9552 Tf 9.299 0 Td [(1 : 875m ;z = )]TJ/F15 11.9552 Tf 9.299 0 Td [(0 : 63438m..................................47 V.13CongurationofForwardProblemwithplanararrayofdipolesontopface andneareldrequestsbelowbottomface..................48 V.14SyntheticforwarddataofneareldscattereddatacomputedbyFEKO.49 V.15VolumetricInversionwithTikhonovRegularization.OptimalInverted CurrentDensity j M inv j foundat x =2mm ;y =2mm ;z =0 : 14429m...50 V.16VolumetricInversionwith L 1 Regularization.OptimalInvertedCurrent Density j M inv j foundat x =2mm ;y =2mm ;z =2mm..........51 V.17AMRInversionwithTikhonovRegularization.OptimalInvertedCurrent Density j M inv j foundat x =34 : 0643mm ;y =34 : 0643mm ;z =34 : 0643mm.52 V.18AMRInversionwith L 1 Regularization.OptimalInvertedCurrentDensity j M inv j foundatx=34.0643mm,y=34.0643mm,z=34.0643mm...53 V.19ConstrainedVolumetricInversionwithTikhonovRegularization.OptimalInvertedCurrentDensity j M inv j foundat x =0 : 4m ;y =3mm ;z = 0 : 145m.ThispointliesalongoneoftheedgesofCube2.........54 V.20CurrentDensitymagnitudeofconstrainedvolumetricinversionwithTikhonov Regularization.GridIndexofmaximumcorrespondstooptimalinverted currentdensity j M inv j .............................55 x
PAGE 11
V.21ConstrainedVolumetricInversionwith L 1 Regularization.OptimalInvertedCurrentDensity j M inv j foundat x =0 : 4m ;y =0 : 4m ;z =3mm. ThispointliesalongoneoftheedgesofCube1...............55 V.22CurrentDensitymagnitudeofconstrainedvolumetricinversionwith L 1 Regularization.GridIndexofmaximumcorrespondstooptimalinverted currentdensity j M inv j .Theeectsofthe L 1 penaltyareclearlydemonstratedinthesesparseinversionresults...................56 xi
PAGE 12
CHAPTERI MOTIVATION Theproblemofdetectingfeaturesandobjectsofinterestthatareobstructedby conductivebarriersispervasiveinmanyengineeringapplications.Theseincludesubsurfaceimagingforgeophysicalexploration[8],detectionofnaturalandmanmade undergroundstructures[26],nondestructiveevaluationforvehiclesandcontainerinspections[25],underwaterdetection[2,16],andseveralotheractiveareasofresearch. Thediculty,however,isthatobjectsofinterestconcealedbehindhighlyconductive mediatypicallycannotbelocatedwiththeaidoftraditionalelectromagneticimaging methods.Forinstance,groundpenetratingGPRsystemsMHz 30kHzelectromagneticradiationcannotpenetratethrough metalliccontainersofeven1 )]TJ/F15 11.9552 Tf 12.968 0 Td [(2mmthicknesswithoutundergoingextremelosses > 300dB.Fortunately,thedrawbackofhighlossescanbeovercomewiththeuseof extremelyandverylowfrequencyELF/VLFelectromagneticsignals f< 30kHz withlargerskindepths inthatconductivemedium[32]. r 2 ! 0 I.1 Here ! isthewaveangularfrequency, istheconductivity,and 0 isthepermeability offreespace.Theseelectromagneticsignalswillattenuatewhilepassingthrougha depth d intheconductivemediaas e )]TJ/F20 7.9701 Tf 6.586 0 Td [(d= .Forinstance,theskindepthinaluminum 10 7 S/mis16 mat100MHz,butisonly16mmat100Hz.UsingtransmittedLF 1
PAGE 13
signalswouldallowfordetectionofthescatteredeldsproducedbyobjectshidden behindtheseconductivebarrierswhileminimizingtheseattenuationlosses. ELF/VLFsignalsserveasapromisingalternativeforthedetectionofobjects hiddenbyconductivebarriers.However,thereisaninherenttradeowithimaging resolution.ConventionalimagingtechniquessuchasGPR,Xraytomography,orCT scansrelyonthefareldraypropertiesofelectromagneticradiationtoreconstruct animage.Fortheseaforementionedimagingmodalities,theimagingresolutionis limitedbytheRayleighCriterionalsoknownasthediractionlimit[11].Specically,thesmallestresolvablefeatureisapproximatelyhalfawavelengthinsize,which determinesthehighestachievableresolution.TheRayleighcriterionisthespatial analogoftheNyquistfrequency[6]andisaninherentfeatureoftraditionalimage reconstruction.Forexample,GPRoperatingat1GHzhasanassociatedresolution limitof15cmaccordingtotheRayleightCriterion.TheproblemwithELF/VLF imagingisthatthewavelengthsaretenstohundredsofkilometers.TheRayleigh Criterionsuggeststhatonlyobjectsthatareseveralkilometersinsizeareresolvable usingELF/VLFradiation.ThislimitisreasonablewhenstudyinglargescalegeologicalfeaturesoftheEarth[29].Ontheotherhand,imagingundergroundbunkers, tunnels,orobjectsinsideshippingcontainersisnotpossiblewithELF/VLFsignals whenusingimagingalgorithmsthatrelyontheraylikefeaturesofelectromagnetic waves.Tocircumventthislimit,itisnecessarytodevelopnovelstrategiesthatdonot relyontypicalelectromagneticwavefeaturesforimaging.Imagingsmallobjectsin theELF/VLFrangewilloccurinthesocalledneareldregionandthecorresponding quasistaticnatureofthesesignalsmustbeusedtoreconstructanimage.Physically, incidentELF/VLFsignalscaninduceloopsofelectriccurrentonanobjectofinterest whichwillinturngeneratequasistaticmagneticeldsintheneareldoftheobject. Thesemagneticeldswillhaveaspatialamplitudeandpolarizationprolethatare specictotheobject'sshapeandsize.Thus,anappropriateELF/VLFimagingalgo2
PAGE 14
rithmcanbedesignedbyleveragingthemappingfromtheobject'sgeometrytothe measuredmagneticeldprole.Evenso,inferringthelocationofhiddenobjectsfrom elddataisinitselfaconsiderablechallengeduetotheillposednessoftheproblem inthisfrequencyregime.Inthiswork,strategiesfordetectingundergroundobjects aswellasobjectsinsidemetalliccontainersusingELF/VLFsignalsareconsidered. 3
PAGE 15
CHAPTERII TheoryofElectromagneticScattering II.1EquivalentDipoleFormulation Ingeneral,theelectriceldthatisscatteredbyaconductiveobjectcanbeformulatedbythevolumeintegralequationVIEgivenbyII.1. ~ E sc = G r ; r 0 J r 0 d r 0 II.1 Here,thequantity J istheequivalentcurrentdensityinthematerialandisgivenby J = E ,where E isthetotalandunknownelectriceldinthematerial.Oncethe electriceldinthematerialandcorrespondingcurrentsisdetermined,thescattered magneticeldcanbecomputedusingII.2. ~ B sc = G B r ; r 0 J r 0 d r 0 II.2 Theproblemofdeterminingtheconductivityproleisingeneralnonlinearsince itisaproductoftheunknownelectriceldandunknownconductivity.However, insteadofinvertingtheconductivitydistributionofanobject,inferringtheequivalent currentprolefrommagneticeldmeasurementsisagoodproxyforrecoveringthe geometryoftheobjectinquestion.Thisapproachhastheaddedbenetoflinearizing theimagereconstructionproblem. Withthisequivalentcurrentreconstructionmethodinmind,formulationofthe LFelectromagneticinversionproblemtakesinspirationfromthesourcelocalization techniquesofelectroencephalographyEEGandotherrelatedmedicalimagingtechniques[14].Specically,electriccurrentsinsidethehumanbodytypicallyoperateat frequenciesthatlessthan1kHz[5].ThepurposeofEEGliketechniquesistoinfer themagnitudeandlocationofthesecurrentsinanoninvasivemanner.Thus,the 4
PAGE 16
EEGimagingmethodologyismathematicallyidenticaltothegeneralneareldimage reconstructionproblemwithELF/VLFsignals. Assuggestedby[10],neuralactivityinthebraincanbemodeledasaseriesof equivalentcurrentdipoleswithdipolemoment ~ d i = d i ^ e i thatislocatedat ~ r dip ; i .Each currentdipolecontributestothevoltagepotentialmeasuredattheelectrodelocation ~ r .Assuch,thefulldescriptionofeachequivalentdipolerequiressixparametersin R 3 .Bysuperposition,thepotentialmeasuredatanelectrodeduetotheseequivalent sourcescanbeexpressedasII.3[14]. m r j = X i g r j ; r dip ; i ; d i = X i g r j ; r dip ; i d i II.3 Assuming N observationlocationsand p equivalentsources,II.3canbeexpressed asalinearsystemrelatingtheelectrodemeasurementstotheequivalentsourcedipole momentsII.4orsimilarlyinitscondensedformasII.5. 2 6 6 6 6 4 m r 1 . . . m r N 3 7 7 7 7 5 = 2 6 6 6 6 4 g r 1 ; r dip ; 1 g r 1 ; r dip ; p . . . . . . . . . g r N ; r dip ; 1 g r N ; r dip ; p 3 7 7 7 7 5 2 6 6 6 6 4 d 1 . . . d p 3 7 7 7 7 5 II.4 M = GD II.5 ThisformulationII.5takesthesameformastheintegralequationformulation ofthescatteredelectromagneticeldsdescribedinII.2.ToformulatetheelectromagneticintegralequationinthesamemannerasinEEGimaging,equationsII.2 andII.1canbediscretizedusingpulsebasisfunctionsII.6andpointmatching testfunctionsII.7onthediscretizedimagingregion. ~ J ~ r = X i ~ i P i ~ r ;P i ~ r = 8 > > < > > : 1 ; ~ r 2 S i 0 ;o:w: II.6 5
PAGE 17
~ T m = j ~ r )]TJ/F19 11.9552 Tf 10.856 0.166 Td [(~ r C m j II.7 Thedipolemomentin R 3 isdescribedbyaweightedvector ~ i overthearbitrarily denedcell S i .Thepointmatchingtestingfunctionsassuressamplingatthecenter ofthe i th cell[23].Withthischoice,theVIEof E sc II.1and B sc II.2transforms intoalinearsystemresemblingII.4. 2 6 6 6 6 4 B sc r 1 . . . B sc r N 3 7 7 7 7 5 = 2 6 6 6 6 4 G B r 1 ; r 1 0 G B r 1 ; r p 0 . . . . . . . . . G B r N ; r 1 0 G B r N ; r p 0 3 7 7 7 7 5 2 6 6 6 6 4 J 1 . . . J p 3 7 7 7 7 5 II.8 B sc = G B J II.9 TheproblemofequivalentcurrentreconstructionusingmagneticeldmeasurementsisnowrecastinanidenticalmathematicalformtotheEEGproblemII.5. However,sincegenerally N>p ,thislinearsystemisoverdeterminedandinversionof II.9mustbesolvedinaleastsquaressense.Accuratereconstructionoftheforward problemfromthismeasurementdata E sc ;B sc requiresregularizationandapriori constraintswhicharediscussedinChapterIII. II.2ScatteringinFreeSpace Considertheproblemofelectromagneticscatteringbyaperfectlyelectrically conductingPECsphereofradius a infreespace.Thisscattererisilluminatedby aTEpolarizedplanewavewithwavenumber k propagatinginthe+^ y directionas showninFig.II.1. 6
PAGE 18
FigureII.1: ScatteringofPECsphereofradius a infreespaceilluminatedbya TEpolarizedplanewave.Inthelongwavelengthlimit,theeldsscatteredbythis canbeapproximatedbytheeldsofa )]TJ/F15 11.9552 Tf 10.007 0 Td [(^ z orientedidealmagneticdipole. Inthelongwavelengthlimit ka 1,thespatialvariationoftheplanewaveis smallandcanbeapproximatedwithauniformeldof H i y = H 0 e )]TJ/F20 7.9701 Tf 6.587 0 Td [(jky ^ z H 0 ^ z = H 0 cos ^ r )]TJ/F15 11.9552 Tf 11.246 0 Td [(sin ^ .TheboundaryconditionsonthePECsphererequirethenormal componentoftheeldstobecontinuousattheinterfaceandasurfacecurrent J s to beinducedtangentiallytothesurface[23]. ^ n ~ H =0 ^ n ~ H = ~ J s II.10 ApplicationofthenormalboundaryconditionII.10,inasphericalcoordinateframe, requiresII.11. ^ r ~ H sc j r = a = )]TJ/F15 11.9552 Tf 11.291 0 Td [(cos H 0 II.11 Applicationofthetangentialboundaryconditionattheinterfacerequirestheinduced surfacecurrent J s tobe )]TJ/F15 11.9552 Tf 10.812 3.155 Td [(^ directed.Computationofthescatteringduethesesurface currentscanbesimpliedbyconsideringtheeldsradiatedbyanequivalent )]TJ/F15 11.9552 Tf 10.008 0 Td [(^ z directedmagneticdipolelocatedatthecenterofthePECsphere[15]. 7
PAGE 19
~ H sc ~ r = 1 4 r 3 [3 m ^ r ^ r )]TJ/F16 11.9552 Tf 11.955 0 Td [(m ]= )]TJ 13.197 8.087 Td [(j m j 4 r 3 [2cos ^ r +sin ^ ]II.12 Thecombinationoftheseresultsgiveexactexpressionsofthemagneticdipolemoment intheneareldasII.13. ~ m = )]TJ/F15 11.9552 Tf 9.299 0 Td [(2 a 3 ~ H i II.13 Theseanalyticalexpressionsfor H sc werecomparedtomethodofmomentssurfaceintegralequationMoMSIEsolutionsobtainedinthecommercialelectromagneticsolverFEKOalongacutline2metersabovethePECsphere z =2mas inFig.II.2. FigureII.2: Comparisonofanalyticalscatteredmagneticeldsduetoidealdipole againstMoMSIEsolutionsofspherescatteringandequivalentmagneticdipole. Theseanalyticalresultsfor H sc agreeconsiderablywiththeMoMsolutionsof thisscatteringproblem.Thissimpleapproachofrepresentingmagneticscattering byconductiveobjectswithequivalentmagneticdipolesisavalidapproximationin theneareldfortheunknownelectromagneticscatterersintheimagingproblem. Sinceamagneticdipoleisidenticaltoaloopofelectricalcurrent,thisapproachis 8
PAGE 20
identicaltosearchingforloopsofcurrentasopposedtoelectricaldipolesthatare inducedonanobject.Thisisaconvenientabstractionthathelpsalleviatesomeof thecomputationalburdenassociatedwithimagereconstruction. II.3ScatteringbyaPECSphereBuriedinaHomogeneousHalfSpace ConsidertheelectromagneticscatteringofaPECsphereilluminatedbyanincidentplanewaveasshowninFig.II.3.ThePECsphereisembeddedwithina homogeneoushalfspacewithparameters ";; atadepthof d meters.Theregion aboveisfreespace. FigureII.3: CongurationofPECsphereofradius a embeddedinhomogeneous halfspaceRegion 2 withconductivity illuminatedbyanincidentplanewave. Forsimplicationoftheanalysis,theplanewaveisassumedtobeaTEpolarized planewavepropagatingwithnormalincidencetotheinterfacesuchthat ^ k 0 = )]TJ/F19 11.9552 Tf 9.298 0 Td [(k 0 ^ z , where k 0 isthewavenumberinfreespace.Itisworthnothingthatsignalsfrom existingVLFnavycommunicationstransmittersaretypicallyTMpolarizedandat grazingincidence[12].Nevertheless,theTEnormalincidencecaseisaconservative estimateforanalyticalcalculations.Duetoanimpedancemismatchbetweenthe halfspaceRegion 2 andfreespaceRegion 1 ,reectionsoccurattheinterface z =0. 9
PAGE 21
0 = r 0 0 377 2 = r = r 0 0 )]TJ/F19 11.9552 Tf 11.955 0 Td [(j ! r j! 0 for ! 0 II.14 ThemagneticeldinRegion 1 isgivenby, H 1 y z = H 0 e jk 0 z )]TJ/F19 11.9552 Tf 11.955 0 Td [(R TE H 0 e )]TJ/F20 7.9701 Tf 6.587 0 Td [(jk 0 z II.15 wherethequantity R TE istheFresnelreectioncoecient[9].Afractionofthe magneticeldgetstransmittedintoRegion 2 withFresneltransmissioncoecient T TE . H 2 y z = T TE H 0 e jk 2 z II.16 Fornormalincidence,theFresnelcoecientsforthemagneticeldarerelatedby 1+ R TE = T TE .Specically,theFresnelreectionandtransmissioncoecientsare givenbyII.17andII.18[9]. R TM = k 0 )]TJ/F19 11.9552 Tf 11.955 0 Td [( 0 k 2 k 0 + 0 k 2 II.17 T TM = 2 k 0 k 0 + 0 k 2 II.18 ThesignaltransmittedintoRegion 2 suersconsiderablelossesinitiallydueto theimpedancemismatchbetweenthetwomediaattheinterfaceandthenfromattenuationasitpropagatesdeeperintothehalfspace.Theattenuationlossesresult fromthecomplexnatureofthewavenumber k 2 =Re k 2 + j Im k 2 . ScatteringbytheburiedPECspherenowwillonlybeduetotheeldsinRegion 2 II.16.However,theapproachesdetailedinSectionII.2arenotentirelyvalidfor thisproblemsincetheeldscanpotentiallyvaryacrossthelengthofthescatterer 10
PAGE 22
duetothebackgroundmedium.SatisfyingtheboundaryconditionsII.10resultsin anonuniforminduced ~ J s overthe )]TJ/F15 11.9552 Tf 10.238 3.155 Td [(^ direction.Theattenuationof H 2 y z overthe lengthofthescattererresultsinastrongerresponseintheilluminatedregionofthe PECsphere z> )]TJ/F19 11.9552 Tf 9.299 0 Td [(d comparedtotheshadowregion z< )]TJ/F19 11.9552 Tf 9.298 0 Td [(d bythesimulation resultsdepictedinFig.II.4. FigureII.4: MoMSIEsolutionof ~ J s ofPECsphereburiedinhomogeneoushalfspace =1S/milluminatedbyTEpolarizedplanewave f =1kHz.Alarger ~ J s isinducedatthetopofsphererelativetothebottomduetoattenuationoverthe lengthofthesphereinsidethehalfspace. Since J s isnotuniformacrossthe )]TJ/F15 11.9552 Tf 10.238 3.155 Td [(^ direction,asinglehorizontalmagnetic dipoleHMD )]TJ/F15 11.9552 Tf 10.09 0 Td [(^ y isnotsucienttodescribetheradiatedscatteredelds.Instead, thismustbeapproximatedwithaseriesofHMDwithdecreasingmagnitudeover thelengthofthesphere.Forthepurposesofthisanalysis,thiswillberoughly approximatedbyasingleHMDwithdipolemoment m . 11
PAGE 23
TheeldsradiatedbyaHMDdipolewithcurrentmoment I ~ ` aregivenby[9]as II.19. H r = )]TJ/F19 11.9552 Tf 9.299 0 Td [(k 2 I + rr k 2 ^ j m j e )]TJ/F20 7.9701 Tf 6.587 0 Td [(jkr 4 r E r = r ^ j! j m j e )]TJ/F20 7.9701 Tf 6.587 0 Td [(jkr 4 r II.19 Forthecaseconsidered,^ istakentobe )]TJ/F15 11.9552 Tf 10.09 0 Td [(^ y .IncludingtheeectofaniteconductivityhalfspacerequiresreformulatingtheGreen'sfunctionsviaaspectralexpansion [9].Specically,theSommerfeldidentitycanbeappliedforaspectralrepresentation ofthepointsourceandisgivenbyII.20. e )]TJ/F20 7.9701 Tf 6.586 0 Td [(jkr r = )]TJ/F19 11.9552 Tf 10.663 8.087 Td [(j 2 1 dk k k z H 0 k e )]TJ/F20 7.9701 Tf 6.587 0 Td [(jk z j z j II.20 Thisrepresentsasphericalwaveexpandedasaseriesofcylindricalwavespropagating inthe^ direction =^ x cos +^ y sin andaplanewavepropagatinginthe^ z directionforallwavenumbersinkspace k 2 = k 2 + k 2 z [9].Theprimaryadvantage oftheSommerfeldapproachisthateachcylindricalwavecanbereectedfromthe airgroundinterfaceviathetraditionalFresnelcoecients.Theneteldisthusa continuoussuperpositionofeachreectedandtransmittedwave.Thisallowsthe totaleldtobewrittenasasingleintegralexpressionthatcanthenbeapproximated withanalyticalornumericalcalculations. 12
PAGE 24
FigureII.5: HMD^ = )]TJ/F15 11.9552 Tf 10.091 0 Td [(^ y servesasaroughapproximationfortheelectromagnetic scatteringbyaPECspherecenteredat z =0. H 2 ;z z = j m j 8 cos 1 dk k 2 H 1 k e )]TJ/F20 7.9701 Tf 6.587 0 Td [(jk 2 ;z j z j )]TJ/F15 11.9552 Tf 14.508 3.022 Td [(~ R TE 2 ; 1 e )]TJ/F20 7.9701 Tf 6.587 0 Td [(jk 2 ;z z +2 d H 1 ;z z = j m j 8 cos 1 dk k 2 H 1 k ~ T TE 2 ; 1 e )]TJ/F20 7.9701 Tf 6.586 0 Td [(jk 0 ;z z + d + jk 2 ;z d II.21 ~ R TE 2 ; 1 = 0 k 2 ;z )]TJ/F19 11.9552 Tf 11.955 0 Td [(k 0 ;z 0 k 2 ;z + k 0 ;z ~ T TE 2 ; 1 = 2 0 k 2 ;z 0 k 2 ;z + k 0 ;z II.22 InII.21, )]TJ/F15 11.9552 Tf 9.298 0 Td [(ischosenfortheprimaryeldforwavespropagatinginRegion 2 at z> 0. Thismethodologycanbeextendedtoaccountformultilayeredplanarmedia bytakingasimilarapproachasinII.21,thatisbyconsideringreectionsand transmissionsineachlayerofinterestandthesourceregion[9].Forsimplicity,only thehomogeneoushalfspaceproblemisconsideredinthiswork. 13
PAGE 25
ComputationoftheSommerfeldIntegralsinII.21areingeneralanalytically intractableandthusrequirenumericaltechniques.Typicalasymptoticmethodsdetailedin[9]arehighlyproblemdependentandpotentiallyposecomplicationsfor scatteringproblemswithlittletonoaprioriinformationaboutthescatterer.Inthis work,thisproblemofevaluatingtheSommerfeldIntegralsissimpliedbyutilizingthe commercialelectromagneticsolverFEKOtoquicklycomputeboththeforwardscatteringproblemaswellastogenerateanumericalGreen'sFunctionresponseNGF foranequivalentmagneticmagneticdipoleburiedinaconductivehalfspace. II.4PECObjectInsideaThinConductiveShell Ofparticularinterestisthedetectionofmetallicobjectsthatarelocatedinside ametalshieldsuchasashippingcontainer.ConsideramodicationtotheelectromagneticscatteringproblemposedinSec.II.2.ThePECsphereofradius a is nowshieldedbyanarbitrarythin,conductiveshellwiththickness= d andisilluminatedbyauniformmagneticeld ~ H = H 0 ^ y .Theproblemisnowconsiderably morecomplexastheshieldattenuatesthesignalthatilluminatesthePECobject, andthescatteredeldshavetoonceagaintraversetheshieldwhichresultsinadditionalattenuation.ThemagneticeldthatisscatteredbythePECobjectcanstill beapproximatedbythatofamagneticdipole,however,thefreespaceresponsemust bescaledby1 = SE 2 H .ThequantitySE H representstheshieldingeectivenessofthe outershield.ThereasonforsquaringSE H isthattheincidentsignalisrstshielded beforeilluminatingthePECsphere,andisagainshieldedwhenexitingtheshield. Forthecaseofasphericalthinshellshield,thereexistsasimpleanalyticalsolution forSE H andisgivenbyII.23. 14
PAGE 26
FigureII.6: Asphericalshieldwiththickness= d andouterradius a illuminated byauniformmagneticeldasin[20].AsinglePECsphereislocatedinsidethe shield. SE H = H 0 H S =cosh k m d + 1 3 k m a + 1 k m a sinh k m d II.23 Thequantity k m correspondstothewavenumberinsidetheshellandisgivenby k m =+ j = ,where isskindepthfromequationI.1.Theeldsscatteredby thePECsphereinsidetheshieldarethusgivenbyII.23,where H 0 aretheeldsin theabsenceoftheshieldwhile H S correspondtotheshieldedelds[20].Equation II.23isanoversimpliedapproximationthatisusedtoestimatethestrengthof thescatteredelds.However,theproblemofmagneticscatteringfromaPECobject thatisinsideametalshieldrequiresanumericalapproach. 15
PAGE 27
CHAPTERIII INVERSEPROBLEMFORMALISM III.1NecessityforRegularization ConsiderthegenerallinearelectromagneticscatteringproblemposedinII.9. Generally,thisproblem,liketheEEGinverseproblem[14],ishighlyillposed.Let II.9berecastasIII.1,with A beingtheforwardoperator, b 2 B beingtheproblem dataset,and x 2 X beingtheunknownmodelparameters.Theseparametersare generallycomplexvaluedandsocanbeassumedtobevectorsin B 2 C m and X 2 C n , where C n = f a + bj j a;b 2 R n g . Ax = b III.1 ThisproblemIII.1canbesolvedexactlyifitiswellposed,i.e.satisesboth Hadamardconditions[31]: 1.Foreach b 2 B ,auniquesolutionexists. 2.Thesolutionisstableunderperturbation,i.e. A )]TJ/F17 7.9701 Tf 6.586 0 Td [(1 isdenedforall B andis continuous. Theelectromagneticscatteringproblemtypicallyfailsbothoftheseconditionsand isHadamardillposed. Withtheequivalentdipoletechnique,innitelymanydipolesarerequiredtoaccuratelymodeltheessentiallycontinuousimagingregion.Discretizingtheimaging regionstillwouldrequiremanymoreequivalentdipolesthantherearemeasurement locations n>m forasucientmodel.Eq.III.1willonlybeconsistentif b 2R A ,sincebydenitiontherange R A = f Ax j x 2 C n g [7].Thisisequivalenttorequiringrank A =rank A j b ,where A j b istheaugmentedmatrix.By Rouche{Capellitheorem[27],thissolutionisuniqueifandonlyiftherank A = n . 16
PAGE 28
Thiscannotoccurfortheunderdeterminedsystemsincerank A
PAGE 29
FigureIII.1: Figurefrom[19].Acollectionof N idealdipolesinfreespaceuniformly spacedwithspacing z alongalineoflength L areusedtocomputethecondition numberoftheforwardoperator G B todeterminesitsstability. Thedipolearrayissurroundedby N r receiversthatareboundtothesurfaceof asphereofradius L= 2.TheGreen'sfunctionofanidealmagneticdipole G B r;r 0 iscomputedforthisconguration.Theconditionnumberofthisforwardoperatoris computedusingIII.2.TheresultsareshownbelowinFig.III.2. FigureIII.2: Conditionnumberof G B iscomputedforacollectionofdipoles N = L= z boundonlineoflength L .For L =1m,thenormalized correspondto wavesof f =300GHz,300MHz,and300kHz,respectively. 18
PAGE 30
Itisevidentthatthesolutionbecomesquicklyunstableatlowfrequencies.Solutionsathigherfrequencywillremainrelativelystablewhiletheminimumresolution z> R FF ,withR FF =1 : 22 2 beingtheRayleighlimit[11].Thehigherfrequency consideredinFig.III.2islikelytobecomeunstableat N 1700dipoles.Forsmall N = L= z ,thisconditionnumberisrelativelysmallduetothesucientspacing betweendipolepairs.Theeldsofanydipoleisdistinguishableinthislimit. Therefore,itisclearthatIII.1failsatleastone,ifnotboth,oftheHadamard conditionsforwellposednessinthisfrequencyregime.Whiletheinverseproblem isillposedanddiculttosolve,anoptimalestimator^ x existsthatbothminimizes theerrorbetweenthedataset b andtheestimate A ^ x andalsorepresentsthemost likelyimage"ordistributionofequivalentdipolesexpectedintheforwardproblem.Atypicalmetricformeasuringerrorbetweenthedataandobservationsisthe unconstrainedordinaryleastsquaresminimizationproblemisgivenbyIII.3. min x jj b )]TJ/F19 11.9552 Tf 11.955 0 Td [(Ax jj 2 2 III.3 However,theleastsquaressolutiondoesn'tenforcecertaindesirablefeaturesin theimagesuchassparseorminimalmodelparameters.Thisisaccomplishedthrough theadditionofaregularizationtermtoIII.3thatactsasamathematicalconstraint totheminimizationproblem. min x jj b )]TJ/F19 11.9552 Tf 11.955 0 Td [(Ax jj 2 2 + jj x jj p p III.4 EquationIII.4representsthegeneralregularizedminimizationproblem,where istheregularizationparameter.ThesecondterminIII.4representsthegeneral p normandisgivenbyIII.5. jj X jj p = X i j x i j p 1 =p III.5 19
PAGE 31
Thechoiceof correspondstoacertainfamilyofsolutionstotheoptimization problem.Twocommonchoicesofpnormintheliteraturearethe p =2Tikhonov regularization[14]and p =1LASSOregularizationorregression[30].Thesepnorm regularizedproblemsareguaranteedtohaveuniqueandgloballyoptimalsolutionsin theproblemspacesincethesefunctionsareconvexfor p 1[7]. III.2 L 2 TikhonovRegularizationTheory Considerthe L 2 normminimizationproblemTikhonovRegularizationproblem. min x jj b )]TJ/F19 11.9552 Tf 11.956 0 Td [(Ax jj 2 2 + jj x jj 2 2 III.6 Itisafundamentalpropertyofconvexoptimizationproblemsthatthelocallyoptimal solutionisalsogloballyoptimal[7].ThegloballyoptimalsolutiontoIII.6hasan analyticalformsincetheobjectivefunction L = jj b )]TJ/F19 11.9552 Tf 12.019 0 Td [(Ax jj 2 2 + jj x jj 2 2 isdierentiable forall x intheproblemspace. L = jj b )]TJ/F19 11.9552 Tf 11.956 0 Td [(Ax jj 2 2 + jj x jj 2 2 = x T A T A + I x )]TJ/F19 11.9552 Tf 11.955 0 Td [(b T Ax + b T b III.7 Termsindependentof x canbeneglectedintheminimizationproblem.Toaidinthe computationofitsderivative,thequadraticterm A T A + I isasymmetricmatrix. d L dx =2 x T A T A + I )]TJ/F19 11.9552 Tf 11.955 0 Td [(b T A =0 = ^ x = A T A + I )]TJ/F17 7.9701 Tf 6.586 0 Td [(1 A T b III.8 Thisregularizationproblemattemptstominimizeallpredictedcomponentsofthe estimator^ x uniformly. III.3 L 1 RegularizationTheory Considerthe L 1 normminimizationproblemIII.9. min x jj b )]TJ/F19 11.9552 Tf 11.956 0 Td [(Ax jj 2 2 + jj x jj 1 III.9 20
PAGE 32
Whiletheobjectivefunctionisconvex,itisnotdierentiableduetothe L 1 norm penalty.Thisisbecausetheabsolutevaluefunction L = x T A T Ax )]TJ/F19 11.9552 Tf 11.955 0 Td [(b T Ax + X i j x i j III.10 Forsimplicity,III.9canberecastasitsdual[7]. min x jj b )]TJ/F19 11.9552 Tf 11.955 0 Td [(Ax jj 2 2 subjectto X i j x i j III.11 TheLASSOleastabsoluteshrinkageandselectionoperator[30]takesanidentical formtotheprimalproblemIII.11. ^ ; ^ =argmin N X i =1 y i )]TJ/F19 11.9552 Tf 11.956 0 Td [( )]TJ/F25 11.9552 Tf 11.955 11.358 Td [(X j j x ij 2 subjectto X j j j j III.12 Forallregularizationparameters,theoptimalestimateof^ istakentobethemean ofthedata y [30].Withoutlossofgenerality,thisconstanttermisoftentakentobe ^ =0. Aclosedformsolutionfortheorthonormaldesigncase X T X = I canbeobtainedbyconsideringsubdomainwheretheabsolutevaluefunction f x = j x j is dierentiablefor x 6 =0. d dx f x = 8 > > < > > : )]TJ/F15 11.9552 Tf 9.298 0 Td [(1 ;x< 0 1 ;x> 0 III.13 Solutionstothemodiedobjectivecanbefoundforeachsubdomain D + and D )]TJ/F15 11.9552 Tf 7.084 4.339 Td [(, where D + isdenedfor j > 0and D )]TJ/F15 11.9552 Tf 11.575 4.339 Td [(for j < 0.Let ^ 0 = y T X betheordinary leastsquaressolutionfortheorthonormaldesigncase. 21
PAGE 33
^ j = 1 2 2 ^ 0 j )]TJ/F19 11.9552 Tf 11.955 0 Td [( ; j > 0 ^ j = 1 2 2 ^ 0 j + ; j < 0 III.14 ThissimplybecomesIII.15fortheunionofthesesubdomains D = D )]TJ/F22 11.9552 Tf 8.665 4.338 Td [([ D + ,where sgn x isthesignfunction.Theconstraineddomainofthisproblemrequiresthe regularizationparametertobe j ^ j j >= 2. ^ j = 1 2 sgn ^ 0 j 2 j ^ 0 j j)]TJ/F19 11.9552 Tf 17.932 0 Td [( III.15 III.3.1Tibshirani'sAlgorithmforSolvingLASSO TibshiranihadproposedanalgorithmforaquadraticprogramtosolvetheLASSO optimizationproblem.Foraxedregularizationparameter,III.12isexpressedas aleastsquaresproblemwith m =2 n inequalityconstraints.Theseinequalityconstraintsareintroducedsequentiallyuntilafeasiblesolutionisobtainedthatsatises theKuhnTuckerconditions[30].Theinequalityconstraint P j j j j isexpressed as T i ,where i for i =1 ; 2 ;:::; 2 n arethe n tuplesoftheform 1 ;:::; 1. Let G E bethematrixwhoserowscontain i for i 2 E ,where E = f i : T i = g is theequalitysetand S = f i : T i < g istheslackset,and ^ 1 bethevectorofonesof thesamelengthas G E . forn=2, G E = 2 6 6 6 6 6 6 6 4 1 2 3 4 3 7 7 7 7 7 7 7 5 = 2 6 6 6 6 6 6 6 4 11 1 )]TJ/F15 11.9552 Tf 9.298 0 Td [(1 )]TJ/F15 11.9552 Tf 9.298 0 Td [(11 )]TJ/F15 11.9552 Tf 9.298 0 Td [(1 )]TJ/F15 11.9552 Tf 9.298 0 Td [(1 3 7 7 7 7 7 7 7 5 III.16 Tibshirani'salgorithmisinitializedwith E = i 0 with i 0 =sgn ^ 0 .Iftheoriginal inequalityconstraint P j j j j issatised,thealgorithmcompletes.Otherwise, theviolatedconstraintisremovedfromtheslacksetandaddedtotheequalityset. 22
PAGE 34
FigureIII.3: PseudocodedescribingTibshirani'salgorithmforsolvingLASSO.This minimizestheobjectiveoftheprimalproblemwiththeconstraintsintheequality set E . Tibshirani'salgorithmisguaranteedtoconvergeinatmost2 n steps.Thisis becausethenaliterateisthesolutiontotheoriginalproblemwhichdoessatisfythe KuhnTuckerconditions.Withtheapplicationofstandardquadraticprogramming techniques,thiscanreducethenumberofstepsrequiredforconvergenceto2 n +1 steps[30].However,itisclearthatthisisimpracticalforlarge n duetothelengthof thetuplesrequiredtomodeltheconstraints.Instead,dierenttechniquesforsolving thequadraticprogramaredesirable. III.4OptimalChoiceofRegularizationParameter Eachchoiceofregularizationparameter > 0correspondstoadierentfamilyofsolutionsoftheregularizedoptimizationproblem.Theoptimalchoiceof iscriticalforsolvingforthetrueforwardmodelparameterscorrespondingtothe imagingprobleminquestion.Choosingthisoptimal isequivalenttominimizing theerrorbetweenthemodelandthedata.However,thedicultyinndingthis arisesfromthelackofexplicitknowledgeoftheforwardproblem.Sincethesetrue forwardmodelparametersaretypicallyunknownintheimagingproblem,statistical andimagingtechniquesmustbeusedtoprovideanestimateofthetruemodelerror.TwotechniquescommonlyappliedintheliteraturearetheLCurvemethodand CrossValidationCVmethod. 23
PAGE 35
III.4.1LCurveMethod TheLCurvemethod,namedforthecharacteristicshapeoftheloglogplotasin Fig.III.4,attemptstominimizethetradeobetweenthegeneralizedregularization bias = jj L x jj andtheresidual = jj b )]TJ/F19 11.9552 Tf 9.833 0 Td [(Ax jj 2 [14].Qualitatively,theoptimalchoice of isgivenbythecornerofthecurve.If istoolarge,theinversionsolutionis dominatedbyregularizationerrorssuchasoversmoothing.Otherwise,if istoo small,thesolutionwillbedominatedbythemodelingerrorsandmeasurementnoise [18]. FigureIII.4: CharacteristicloglogplotofLCurveMethodforTikhonovRegularization.Thecornerorclosestpointtotheoriginrepresentstheoptimalparameter . Therearetwomethodsoflocatingthiscorner[18]. 1.Findingthecornerisequivalenttondingtheclosestpointonthecurvetothe origin.Thismetricdiersbasedonthechoiceofregularization.ForTikhonov regularization,thismetricis d = + 2 . 2.Thecornerisalsothepointofmaximumcurvature .Thisisindependentof thechoiceofregularizationsince isapurelygeometricquantity.Assuming 24
PAGE 36
theLcurveisasmoothcurveandthat^ =log jj b )]TJ/F19 11.9552 Tf 12.539 0 Td [(Ax jj 2 and^ =log jj x jj 2 fortheTikhonovRegularization,thenthecurvatureisgivenbyIII.17,where dierentiationistakenwithrespectto . = ^ 0 ^ 00 )]TJ/F15 11.9552 Tf 13.023 0 Td [(^ 00 ^ 0 ^ 0 2 +^ 0 3 = 2 III.17 III.4.2CrossValidationMethod Thecrossvalidationmethodisastatisticaltechniqueofestimatingthemean squarederrorMSEofthesolutionwithoutexplicitknowledgeofthetruesolution. OnlykfoldCVisconsideredfortheseproblems[22].Thisestimateconvergestothe trueMSEsincethesesamplessharethestatisticalpropertiesofthedataset.The kfoldCVinvolvesdividingthedataset b into k equallysizedpartitions.Oneofthese partitionsiskeptasaholdoutvalidationsetandtheremaining k )]TJ/F15 11.9552 Tf 12.264 0 Td [(1areusedas trainingsets,where k istypicallytakentobe5or10. FigureIII.5: Equipartitionofdatasetfor5foldCV.Thisprovidesanestimate oftheMSEofthesolutionvector x thatcanbeusedtooptimizetheregularization parameter . Foraxedregularizationparameter,theregularizedoptimizationproblemisthen solvedusingthe k )]TJ/F15 11.9552 Tf 12.605 0 Td [(1trainingsets.Theoptimalsolution^ x isprojectedontothe holdoutsetby b 0 = A 0 ^ x ,where A 0 istheforwardoperatorcorrespondingtotheholdout 25
PAGE 37
set.Thisprocedureisrepeatedfortheremainingsetsuntileachofthepartitionsis usedasaholdoutset[22]. ~ CV = 1 k k X i =1 [ b 0 i )]TJ/F19 11.9552 Tf 11.955 0 Td [(A 0 i ^ x i ] 2 = 1 k k X i =1 MSE i III.18 III.5ImagingResolutionLimits Thesmallestdetailsabletoberesolvedbytheimagingsystemishighlydependent onthesourceandtheregularizationtechniquesconsidered.Inimagingsystemsthat operateinthefareld ka 1,thisbestresolutionisdeterminedbytheRayleigh limit,givenbyIII.19[11]. R ff =1 : 22 2 III.19 Spatialfeaturesthatarelargerthan R ff areresolvablebythesystem. FigureIII.6: MinimumSpatialResolutionforsignalpropagatinginvarioushomogeneousniteconductivitygrounds.Thisassumesnoattenuationofthesignal. ItisevidentfromFig.III.6thattheobjectsofinterestwhicharetypicallysizes ontheorderof1marenotresolvableatallbyanimagingsystemoperatinginthefareldassuminganoiselessandlosslesssignal.Theinformationofthesesmallspatial 26
PAGE 38
featuresarecontainedintheevanescentwaveswhichareneglectedinthefareld approximation[11].However,thisneareldresolutionlimitwilldiergreatlyfor eachregularizationtechniqueapplied.Harid[19]givestheminimumresolutionlimits foranimagingsystemsolvablebyanordinaryleastsquaresinversion. FigureIII.7: From[19].Thissuggestsforanimagingregionofsize L m,the minimumresolutionobtainableforanordinaryleastsquaresinversionwhenthe forwardoperator A isfullrankandinvertibleis 30 )]TJ/F15 11.9552 Tf 11.955 0 Td [(40 L mm. TheresolutionlimitsattainablebyaTikhonovor L 1 regularizedsolutionarenot examinedbutareassumedtobeequal,ifnotbetter,thanthoseforanunregularized solutionasin[19].Intheproblemswhereafullvolumetricinversionisconsidered, theresolutionoftheuniformgridisgreaterthanthelimitposedfortheordinary leastsquaressolutionandareassumednottobelimitedbythisresolutionlimit. 27
PAGE 39
CHAPTERIV ALGORITHMDESCRIPTION IV.1HighLevelOverviewofAlgorithm Thealgorithmsdescribedattempttosolvetheinverseproblemforthelowfrequencyneareldelectromagneticscatteringproblem.Itreliesonthecreationof syntheticforwarddata B data andtheNGF G B inthecommercialelectromagnetic solverFEKO[3].Thisinversionprcessisagnostictothechoiceofelectromagnetic solver.Thissolverwaschosenduetotheavailabilityandadvantagesofthemoments methodMoMelectromagneticsolveroverothercommercialsolvers.TheMoM solverhassimpleandecienttreatmentsofidealdipolesourcesandonlyrequires meshingofthesurfaceelementstosolvethescatteringproblem.Theinverseproblem issolvedeitheranalyticallyfortheTikhonovregularizedproblem,whichisoptimized bytheLCurvemethod,orsolvedwiththeaidoftheCVXMATLABtoolbox[13] forthe L 1 normregularizedproblem,optimizedbytheCVmethod.Twotechniques forthegenerationoftheimaginggridarepresented:asimplestructuredvolumetric generationoftheimaginggridandaniterativeandlesscomputationallyintensive generationoftheimaginggridtermedadaptivemeshrenement"AMR. 28
PAGE 40
FigureIV.1: Flowchartdescribinghighleveldesignofsoftwareforbothvolumetric meshandadaptivemeshrenementprocedures. IV.2GenerationofSyntheticData TheforwardproblemisfullymodelledinFEKO.Onlythreecasesareconsidered fortheforwardproblem:ithescatteringproblemposedinSec.II.2,iiscattering fromaPECsphereburiedinaconductivehalfspaceasinSec.II.3,andiiithe scatteringfromPECobjectslocatedinsideinametalliccubicshield.FigureIV.2 showsanexampleforwardsimulationsetupoftwoPECcubeslocatedinsidea1m 1m 1mconductivecubicshellof2mmthicknessandconductivity =10 7 S/m. FigureIV.4showsthesetupofaPECsphereinfreespaceaswellasthecaseofa conductivehalfspacehalfspacenotvisibleinsimulationCADle. 29
PAGE 41
FigureIV.2: ForwardProblemGeometryforMetallicBoxwithtwoPECCubes. Themetallicshellhas =10 7 S/mandthickness=2mm.Receiversareplacedon planesparalleltoeachfaceosetby5mm.Anidealmagneticdipolesourceisplaced overeachfacethatoperatesintheLFregime. TheprocessofgeneratingsyntheticforwarddataintheFEKOsolverisgeneral andsotheinversionprocesscanpotentiallybeextendedtoanyclassofelectromagneticinversescatteringproblems.Thegeometryoftheforwardproblem,which includesboththeunknownscatterersandtheknownproblemboundaries,ismodeled withtheCADtoolsprovidedintheFEKOsolver.ThissolverthensolvestheelectromagneticscatteringproblemusingtheSurfaceIntegralEquationMoMSIEMoM approach[3].Thesurfacesofthemodeledgeometryareapproximatedbyameshof triangularelementsbyFEKO. 30
PAGE 42
FigureIV.3: Surfacesaremeshedwithmanytriangularelements.ThebasisfunctionsusedtosolvetheSIEaremappedtotheedgesoftheseelementsandusedto computesurfacecurrentdensities. Whileincreasingthenumberofthesemeshelementsimprovestheaccuracyofthe model,itdoessoataincreasedcomputationalcost.Materialpropertiessuchas conductivityandthicknessareappliedtotheapplicablesurfacesofthemodel.Last, thenearelddataofthescatteredmagneticeldsarerequestedat N r arbitrary pointsplacedoutsidetheproblemspace. FigureIV.4: GenerationofForwardProblemofPECSpherelocatedinfreespace handledinFEKO.TheunknownPECscatterersareplacedinthismodelwiththe knownproblemboundaries.Thescatterednearelddataarerequestedat N r points. 31
PAGE 43
IV.3ComputationofNumericalGreen'sFunctionNGF ClosedformexpressionsfortheGreen'sfunctionresponseofthehalfspaceand metallicboxproblemsareunavailableduetotheinfeasiblilityofsolvingforthese analytically.Instead,computationofthesearehandlednumericallywiththeaidofthe LUAscriptingtoolsavailableinFEKO.Assumeasingleidealmagneticdipolewith unitdipolemomentisplacedat r 0 i anditsmagneticeldisobservedat r i .Then,the Green'sfunctionresponse G B forthatcollectionofsourceandobservationlocations r 0 i ;r i arecomputedfromIV.1byiteratingthroughallpossibleunitmagnetic dipolesalong^ andsolvingfortherespectivematrixentry. 2 6 6 6 6 4 B sc x B sc y B sc z 3 7 7 7 7 5 = 2 6 6 6 6 4 G B xx G B xy G B xz G B yx G B yy G B yz G B zx G B zy G B zz 3 7 7 7 7 5 2 6 6 6 6 4 m x m y m z 3 7 7 7 7 5 IV.1 Thiscomputationiscomputedforallobservationlocations N r andsourcelocations N ,resultingin G B ,amatrixofsize3 N r 3 N .TheuseofnumericalGreen'sfunction isausefulmethodtoincludetheeectsofthebackgroundshieldhalfspaceorshell withouthavingtoexplicitlymodelthebackgroundintheinversionalgorithms. IV.4InversionAlgorithm Threedierenttypesofinversionmethodsareattemptedtosolvefortheunknown modelparameters M inv .Eachattemptstoimposecertainfavorablefeaturesinthe modelspacetoaccuratelyinverttheprobleminquestion. IV.4.1TikhonovRegularization Tikhonovregularizedinversion,or L 2 Regularization,attemptstominimizethe 2normofthesolutionvector M inv uniformlyIV.2. min x jj B data )]TJ/F19 11.9552 Tf 11.956 0 Td [(G B M inv jj 2 2 + jj M inv jj 2 2 IV.2 32
PAGE 44
SolutiontoIV.2isstraightforwardsincetheminimizationproblemcanbesolved analyticallyasdescribedinSectionIII.2.Theoptimalregularizationparameteris foundthroughtheLCurvemethodasdetailedinSecIII.4.1.TheRegularization ToolspackageforMATLABisutilized[17]toaidincomputationof .Thistreatsthe discreteLcurveasadierentiablecurvetocomputethepointofmaximumcurvature . IV.4.2 L 1 Regularization The L 1 regularizedinversionproblemattemptstoenforcesparsityinthesolution vector M inv IV.3. min x jj B data )]TJ/F19 11.9552 Tf 11.956 0 Td [(G B M inv jj 2 2 + jj M inv jj 1 IV.3 Sincethe1normpenaltyisnotdierentiable,itcannotbecomputedanalyticalsolutiontoIV.3.ThisisfurtherdescribedinSecIII.3.Instead,theCVXtoolboxfor MATLABisimplementedinthealgorithmtohandlebothconstrainedandunconstrainedformsofIV.3tosolve M inv [13].Theregularizationparameterisoptimized throughCVmethodasinSecIII.4.2. IV.4.3IterativeLeastSquaresSolution Thismethodgreatlydiersfromtheprevioustwointhatitappliesnotraditional regularizationtotheminimizationproblem.Instead,ittakes N dip dipoles,where N dip
PAGE 45
IV.5GenerationofImagingGrid Theimaginggridisadiscretizationofthecontinuousforwardmodelspaceor problemspace.Eachcellintheimaginggridrepresentsavoxelvolumetricpixel inthenalimage.Increasingthenumberofcells N improvestheresolutionofthe imageaswellasincreasingthecomputationalcostsofconstructingtheNGFand solvingtheinverseproblem.Thevolumetricreconstructionisthesimplestmethodof constructingthisimaginggrid.Thisinvolvesthegenerationofauniformrectangular gridinsidetheboundsofthephysicalmodelspace.If N = N x N y N z ,thenthe resolutionofthisgridisgivenbyIV.5. x = j x ub )]TJ/F19 11.9552 Tf 11.955 0 Td [(x lb j N x )]TJ/F15 11.9552 Tf 11.956 0 Td [(1 y = j y ub )]TJ/F19 11.9552 Tf 11.955 0 Td [(y lb j N y )]TJ/F15 11.9552 Tf 11.955 0 Td [(1 z = j z ub )]TJ/F19 11.9552 Tf 11.955 0 Td [(z lb j N z )]TJ/F15 11.9552 Tf 11.955 0 Td [(1 IV.5 Itisnecessarytorequiretheseupperandlowerbounds r lb r r ub beosetfromthephysicalproblemboundaries.Intersectinggridcellswiththephysical boundariesoftenresultsinsingularvaluesoftheGreen'sfunctionresponse.However,theuseofthisvolumetricgridrequireslargecomputationalcoststoachievean adequateresolution.Additionalmodelingerrorstotheinversioncouldarisefromthe inappropriatelychosengridcells.Tocounteracttheseadverseconsequencesfromthe volumetricinversion,anecienttechniquecalledadaptivemeshrenement"AMR isimplemented. IV.5.1AdaptiveMeshRenementAlgorithm ThemotivationoftheAMRalgorithmistousetheoutputsofaverycoarse meshtypically N =8toaidinlocalrenementsofthegrid.Thealgorithmofa structuredAMRisasfollows: 34
PAGE 46
1.Initializewith N =8gridcellswitheachgridcelllocatedatthecenterofeach octantofthefullmodelspace. 2.Solvetheregularizedinversionproblemforthisconguration. 3.Renetheimaginggridwitharenedsubgridwith N =8atthelocationofthe strongestresponse.Thesizeofnewgridcellsishalfofthesizeoftheoriginal gridcell. 4.Repeatuntilthesizeofthesenewgridcellsbecomesnegligiblecomparedtothe modelspace.Inexperiments,thisthresholdwastypicallychosentobe15%of thedimensionsoftheglobalmodelspace. Inthisway,eachiterationeectivelyadds7gridcellstothemodelsincetheoriginal gridcellisremovedduringtherenementprocess.Thealgorithmforasingleiteration isdemonstratedinFig.IV.5. FigureIV.5: 2DCrosssectionofinitialcoarsegridleft.Distancebetweengrid centers dr = j r C i )]TJ/F19 11.9552 Tf 12.017 0 Td [(r C j j ; for j 6 = i ,istakentobethesizeofthegridcell.Atthenext iterationright,thecellwiththestrongestresponseisremovedandreplacedwitha subgridwith8cellsofsize dr= 2. Theselectioncriteriontypicallyconsideredforthegridrenementisthegreatest j M inv j . 35
PAGE 47
CHAPTERV RESULTS Summarizedbelowaretheinversionresultsfromtheapplicationofthealgorithm andtechniquesdescribedtothethreecasesconsidered.Itisworthnotingthatthese resultsdonotdemonstrateaccuratereconstructionsoftheforwardproblem,thus indicatingthenecessityforfurtheroptimizationoftheinversionprocess. V.1InversionresultsforFreeSpaceInversion ThisinversionalgorithmwasvalidatedwiththescatteringbyaPECsphereof radius a =20cmlocatedinfreespaceat x;y;z = ; 0 ; )]TJ/F15 11.9552 Tf 9.298 0 Td [(0 : 5m.Thescattereris illuminatedbyaplanewavesourcewith ~ E i = )]TJ/F15 11.9552 Tf 9.298 0 Td [(1V/m^ x propagatingwithnormal incidencetothe xy planesimilartotheincidenteldsdescribedinSec.II.3.A factorof e j 2 ft ,with f =25kHz,isassumed.Thenearelddataissampledalong pointsontheplane z =+10cmon x 2 [ )]TJ/F15 11.9552 Tf 9.298 0 Td [(10m ; 10m]and y 2 [ )]TJ/F15 11.9552 Tf 9.298 0 Td [(10m ; 10m].The forwarddataofthescatteringprobleminfreespacelimitedtothesubdomain x 2 [ )]TJ/F15 11.9552 Tf 9.298 0 Td [(2m ; 2m]and y 2 [ )]TJ/F15 11.9552 Tf 9.299 0 Td [(2m ; 2m]isshowninFig.V.1.Thecharacteristicshapeofthis elddataindicatesthestrongestresponseinthecenter x =0 ;y =0.Thissuggests thepresenceofthePECscattererneartheorigin,thoughthedepthandgeometry cannotbeconclusivelydeterminedfromdirectexaminationofthisforwarddata. 36
PAGE 48
FigureV.1: SyntheticforwarddataofneareldscattereddatacomputedbyFEKO domainlimitedto x;y 2 [ )]TJ/F15 11.9552 Tf 9.299 0 Td [(2m ; 2m]. FigureV.2: VolumetricInversionwithTikhonovRegularization.OptimalInverted CurrentDensity j M inv j foundat x = )]TJ/F15 11.9552 Tf 9.298 0 Td [(2 : 5m ;y = )]TJ/F15 11.9552 Tf 9.299 0 Td [(2 : 5m ;z = )]TJ/F15 11.9552 Tf 9.299 0 Td [(0 : 1m 37
PAGE 49
AvolumetricinversionregularizedwiththeTikhonovpenaltyisperformedshown inFig.V.2.Auniformgridof8 8 8cellsisconstructedonthedomain x;y 2 [ )]TJ/F15 11.9552 Tf 9.298 0 Td [(2 : 5m ; )]TJ/F15 11.9552 Tf 9.298 0 Td [(2 : 5m]and z 2 [ )]TJ/F15 11.9552 Tf 9.299 0 Td [(2m ; 0 : 1m].Thelargestmagnitudeinvertedcurrentdensity wasselectedastheoptimalvalue.However,theoptimalresponsewaslocatedatthe topofthediscretizedimaginggridat x = )]TJ/F15 11.9552 Tf 9.298 0 Td [(2 : 5m ;y = )]TJ/F15 11.9552 Tf 9.299 0 Td [(2 : 5m ;z = )]TJ/F15 11.9552 Tf 9.298 0 Td [(0 : 1m.Since theLCurvemethodhadselectedsuchasmallvaluefortheregularizationparameter =5 : 4163 e )]TJ/F15 11.9552 Tf 13.001 0 Td [(08,theresultsarelikelydominatedbymodelingerrors,sinceall j M inv j intheimaginggridareallrelativelylargeanduniform.Choosingalarger regularizationparameterbiasestheinversionresultsofthisdatasettowardthetop oftheimagingregionclosesttothereceiverplaneatz=10cm.Thebiastowards cellsclosesttothemeasurementplaneintheoverregularizedinversionisduetothe minimizationofTikhonovpenalty.Sincethemagnitudeoftheequivalentcurrent densityisminimized,thesourcesmustbeinvertedclosesttothemeasurementplane toproducethemodeltheforwardelddataaccurately.TheTikhonovpenaltyisnot suitableforthecongurationofthisvolumetricinversion. FigureV.3: VolumetricInversionwith L 1 Regularization.OptimalInvertedCurrent Density j M inv j foundat x = )]TJ/F15 11.9552 Tf 9.299 0 Td [(1m ;y = )]TJ/F15 11.9552 Tf 9.299 0 Td [(1m ;z = )]TJ/F15 11.9552 Tf 9.298 0 Td [(2m 38
PAGE 50
ConstructedsimilarlytoFig.V.2,avolumetricinversionregularizedwith L 1 penaltywasperformedasdepictedinFig.V.3.TheCVmethodhadselectedthe smallest tested,sincetheestimatedcrossvalidationcurvewasmonotonicallyincreasingoverthetestedregularizationparametersforthisdataset.Thisbehavior islikelyduetoeitherthatthe L 1 penaltyisnotapplicabletothisproblemorthat adierentchoiceofregularizationparametersisnecessaryforthisconguration.It couldbethecasethattheoptimalregularizationparametersliedoutsidetherangeof tested )]TJ/F17 7.9701 Tf 6.587 0 Td [(5 1.Thischoiceofsmallregularizationparameterhadnotenforced sparsityinthesolutionssincethecurrentdensitieswererelativelyuniformacrossall gridcellsandthisislikelytonotoccurforsmallerregularizationparameters.Unlike theTikhonovinversion,the L 1 regularizedprobleminvertedtheoptimalcurrentdensitytoalocationatthebottomoftheimaginggrid x = )]TJ/F15 11.9552 Tf 9.298 0 Td [(1m ;y = )]TJ/F15 11.9552 Tf 9.299 0 Td [(1m ;z = )]TJ/F15 11.9552 Tf 9.298 0 Td [(2m. Itisinconclusivewhetherthisisbiasedtowardsthebottomoftheregionorthatthis techniqueisunsuitableforavolumetricinversionofthefreespaceproblem. FigureV.4: VolumetricInversionwithIterativeLeastSquares.OptimalInverted CurrentDensity j M inv j foundat x = )]TJ/F15 11.9552 Tf 9.298 0 Td [(1m ;y = )]TJ/F15 11.9552 Tf 9.299 0 Td [(1m ;z = )]TJ/F15 11.9552 Tf 9.298 0 Td [(2m 39
PAGE 51
Thevolumetricinversionwasperformedusingtheiterativeleastsquaresmethod asshowninFig.V.4.Thissolutionishighlyunsuitedforthisdataset.ThenormalizedMSEIV.4wascomputedacrossallgridcellsintheuniform8 8 8cell gridandeachwerefoundtobeidentical.DuetoidenticalMSEvalues,theminimum indexchosencorrespondedtotherstindexinthevectorof M inv . FigureV.5: AMRInversionwithTikhonovRegularizationafter5Iterations.OptimalInvertedCurrentDensity j M inv j foundat x = )]TJ/F15 11.9552 Tf 9.298 0 Td [(1 : 875m ;y = )]TJ/F15 11.9552 Tf 9.299 0 Td [(1 : 875m ;z = )]TJ/F15 11.9552 Tf 9.298 0 Td [(0 : 15937m TheTikhonovregularizedinversionwasperformedwiththeAMRalgorithmas depictedinFig.V.5.After5iterationsofthemeshrenementalgorithm,theoptimalinvertedcurrentdensitywaslocatedat x = )]TJ/F15 11.9552 Tf 9.299 0 Td [(1 : 875m ;y = )]TJ/F15 11.9552 Tf 9.299 0 Td [(1 : 875m ;z = )]TJ/F15 11.9552 Tf 9.299 0 Td [(0 : 15937m.Thiswasnottestedformorethan5iterationssincetheoptimalvalues quicklydivergedawayfromthelocationofthescattererintheforwardproblem.Furtheriterationsappearedtobiastowardsthetopoftheimagingdomain.Thisbiasis likelyduetothesamereasonsthatthevolumetricinversionhadposedinadequatefor thisproblem:locationsclosesttomeasurementplanearethosethatsimultaneously minimizetheTikhonovpenaltyandtheresidualbetweenthedataandtheprediction. 40
PAGE 52
ThiscouldalsoindicatethatthecurrentmetricappliedtotheAMRalgorithmisnot anadequatemethodofselectingtherenementsubdomain. FigureV.6: AMRInversionwith L 1 Regularizationafter5Iterations.OptimalInvertedCurrentDensity j M inv j foundat x = )]TJ/F15 11.9552 Tf 9.298 0 Td [(1 : 875m ;y = )]TJ/F15 11.9552 Tf 9.298 0 Td [(1 : 875m ;z = )]TJ/F15 11.9552 Tf 9.299 0 Td [(0 : 15937m The L 1 regularizedinversionwasperformedwiththeAMRalgorithmasdepicted inFig.V.6.After5iterations,theresultsindicatedthatthesameoptimalchoiceof invertedcurrentdensitywasselectedasintheTikhonovregularizedproblem.The convergenceoftheseresultsislikelyduetotheselectioncriteriaforsubdomainrenement,whichlikelyresultsinabiastowardsthetopoftheimagingdomain.Another likelyexplanationisthatduetotheunsuitabilityofthe L 1 regularizationtothis problemduetotheabsenceofalocalminimumintheCVcurve,theinversion isnobetterthandoinganordinaryleastsquaresinversioni.e.minimizingtheunregularizedminimizationproblemIII.3duetotheabsenceoflocalminimainthe domain. Aproposedtesttoevaluatethesuitabilityoftheregularizationtechniquetoa datasetwouldbetoexaminethebehavioroftheCVcurveas ! 0.IftheCV 41
PAGE 53
curvedecreasesasymptoticallyinthislimit,thentheregularizationmethodisno betterthananordinaryleastsquaressolutionduetotheabsenceoflocalminimain thecurve. Theseinconclusiveinversionresultsforthefreespaceproblemarelikelyattributabletothechoiceofthemeasurementplane.Thelackofvariationinthe ^ z directioncoulddirectlyberesponsiblefortheseissues.Theaddeddiversityin measurementscouldbecrucialinresolvingdepthinformation,thoughnoresultsare availableforthisscatteringproblemtoverifytheeectsofdierentchoicesofmeasurementgeometry. V.2InversionresultsforHomogeneousHalfSpaceInversion TheforwardproblemwasconstructedidenticallytothefreespaceprobleminSec. V.1.AgroundplanesimulatedwiththeexactSommerfeldintegralswasincludedat z =0with =10mS/m.Theforwarddatasharesthecharacteristicshapeasinthe freespaceproblem,thoughthemeasuredeldsaretwiceaslargeduetoreection fromtheinterfacebetweenthefreespaceregionandtheniteconductivityhalfspace at z =0.Similarly,thelocationofthePECscatterercanbedeterminedtobenear theorigin x =0 ;y =0duetothelocationofthemaximuminthemagneticeld data.Asinthefreespaceproblem,itisnotclearwhatthegeometryanddepthof thePECscattererarefromdirectexaminationoftheforwarddata. 42
PAGE 54
FigureV.7: SyntheticforwarddataofneareldscattereddatacomputedbyFEKO domainlimitedto x;y 2 [ )]TJ/F15 11.9552 Tf 9.299 0 Td [(2m ; 2m]. FigureV.8: VolumetricInversionwithTikhonovRegularization.OptimalInverted CurrentDensity j M inv j foundat x = )]TJ/F15 11.9552 Tf 9.298 0 Td [(2m ;y = )]TJ/F15 11.9552 Tf 9.299 0 Td [(2m ;z = )]TJ/F15 11.9552 Tf 9.298 0 Td [(0 : 91429m 43
PAGE 55
ThevolumetricinversionregularizedwiththeTikhonovpenaltyisperformedas depictedinFig.V.8.Theoptimalcurrentdensitywaslocatedat x = )]TJ/F15 11.9552 Tf 9.299 0 Td [(2m ;y = )]TJ/F15 11.9552 Tf 9.298 0 Td [(2m ;z = )]TJ/F15 11.9552 Tf 9.298 0 Td [(0 : 91429m.Unliketheresultsofthefreespaceproblem,thiswasbiased towardsthebottomoftheimaginggridwithasimilarregularizationparameter . Thebehaviorofbiasingtowardstheboundariesoftheimaginggridcouldbehighly problemdependent.Testingforoverttinginthisproblembiasesthelocationof theoptimalinversionresultstothebottomoftheimaginggrid.Thereasonforthe disparityinthislocalizationbehaviorbetweenthefreespaceandconductivehalfspaceproblemarenotclear. FigureV.9: VolumetricInversionwith L 1 Regularization.OptimalInvertedCurrent Density j M inv j foundat x =2m ;y =2m ;z = )]TJ/F15 11.9552 Tf 9.298 0 Td [(2m Thevolumetricinversionregularizedwiththe L 1 penaltyisperformedasdepicted inFig.V.9.Theoptimalinvertedcurrentdensitycomputedisidenticaltotheone foundintheTikhonovregularizedproblem.Theresultsareagaininconclusiveasto whetherthisprimarilyduetomodelingerrorsoftheforwardoperator G B ordueto thechoiceofcelllocationsintheuniformgrid.However,theCVcurveforoptimizing 44
PAGE 56
thisproblemsharedthesamebehaviorasthefreespaceproblem.Thisagainsuggests thatthe L 1 penaltyissimplynotsuitabletotheproblemofscatteringofaPECsphere infreespaceorthescatteringofaPECsphereburiedinahomogeneoushalfspace. FigureV.10: VolumetricInversionwithIterativeLeastSquares.OptimalInverted CurrentDensity j M inv j foundat x = )]TJ/F15 11.9552 Tf 9.298 0 Td [(2m ;y = )]TJ/F15 11.9552 Tf 9.299 0 Td [(2m ;z = )]TJ/F15 11.9552 Tf 9.298 0 Td [(2m Thevolumetricinversionwasperformedusingtheiterativeleastsquaresmethod asshowninFig.V.10.Asinthefreespaceproblem,thecomputedMSEofthe ordinaryleastsquaressolutionforeachgridcellwasidentical,givingtherstelement inthevector M inv .Attemptstotthisproblemwithindividualdipolesarelikely unsuitableforeitherithecellsofuniformgridarenotcorrectlychosenastheydo notcoincidewiththeinteriorofthePECsphereoriitheforwardoperatorfora singledipoleisillconditionedandaleastsquaressolutionservesasapoortforthe dataset. 45
PAGE 57
FigureV.11: AMRInversionwithTikhonovRegularizationafter5Iterations. OptimalInvertedCurrentDensity j M inv j foundat x = )]TJ/F15 11.9552 Tf 9.299 0 Td [(1 : 875m ;y = )]TJ/F15 11.9552 Tf 9.298 0 Td [(1 : 875m ;z = )]TJ/F15 11.9552 Tf 9.298 0 Td [(0 : 87187m TheTikhonovregularizedinversionwasperformedwiththeAMRalgorithmas depictedinFig.V.11.After5iterations,theoptimalinvertedcurrentdensitywas locatedat x = )]TJ/F15 11.9552 Tf 9.298 0 Td [(1 : 875m ;y = )]TJ/F15 11.9552 Tf 9.299 0 Td [(1 : 875m ;z = )]TJ/F15 11.9552 Tf 9.298 0 Td [(0 : 87187m.Whilethelocalizationin the xy planewaspoor,theinversionresultsdemonstratesignicantimprovementfrom thevolumetricprobleminthelocalizationofthedepthofthePECsphere.Theerror inthedepthfromthebottomofthesphereis z =0 : 17187m.Improvementinthe localization xy planecouldbeattainedthroughtheinclusionofphysicalconstraints intheinversionproblem. 46
PAGE 58
FigureV.12: AMRInversionwith L 1 Regularizationafter5Iterations.OptimalInvertedCurrentDensity j M inv j foundat x = )]TJ/F15 11.9552 Tf 9.298 0 Td [(1 : 875m ;y = )]TJ/F15 11.9552 Tf 9.298 0 Td [(1 : 875m ;z = )]TJ/F15 11.9552 Tf 9.298 0 Td [(0 : 63438m The L 1 regularizedinversionwasperformedwiththeAMRalgorithmasdepicted inFig.V.12.fter5iterations,theoptimalinvertedcurrentdensitywaslocatedat x = )]TJ/F15 11.9552 Tf 9.298 0 Td [(1 : 875m ;y = )]TJ/F15 11.9552 Tf 9.299 0 Td [(1 : 875m ;z = )]TJ/F15 11.9552 Tf 9.299 0 Td [(0 : 63438m.Theinversionlocalizedthedepthof theoptimalcurrentdensitytotheinteriorofthePECsphere,thoughlocalizationin the xy planeispoorasintheresultsoftheTikhonovregularizedinversion. Theimprovementinlocalizationforthesetwocasescouldbethoughtofasthe inclusionofanadditionalconstraintintothegeneraloptimizationproblemIII.4. Eectively,theAMRalgorithmservesasanadditionalregularizationpenaltybyconstrainingthedomainoftheproblemtocoarsegridcells.Thetechniquehaspotential butrequiresoptimizationtohandleissuessuchasdivergenceintherenementprocessanddeterminingtheoptimalamountofiterationsrequiredtoaccuratelyinvert theforwardmodelparameters. V.3InversionresultsforConductiveCubicShellInversion Thenalproblemconsideredisimagingthroughathin,highlyconductivecu47
PAGE 59
bicshellofsize1m 3 thickness=2mmwith =10 7 S/m.Twosmaller PECcubeswithsidelength ` =0 : 2mareplacedinthebottomoftheshellat x 1 ;y 1 = : 5m ; 0 : 5mand x 2 ;y 2 = : 5m ; 0 : 1m.ThisissothattheCube1 centercubeonlymakescontactwiththebottomfacewhiletheCube2makescontactwiththebottomand y =0faces.Thecongurationisilluminatedbyaplanar arrayofninedipoles5mmabovethetopfaceasinFig.V.13. FigureV.13: CongurationofForwardProblemwithplanararrayofdipoleson topfaceandneareldrequestsbelowbottomface. Theneareldscattereddataisrequested1mmbelowthebottomface.The outlinesofthePECcubesincontactwiththisfaceclearlyappearinthemeasured elddataasinFig.V.14.Suchinformationisvaluableinestablishingphysical constraintstotheinversionproblemforobjectsnearorincontactwiththemetallic shell.Thiswasusedtoconstraintheimagingdomainsothat M inv isnonzeroonly forgridcellscorrespondingtothisoutline:thissubdomainisgivenby x 1 ;y 1 2 [0 : 4m ; 0 : 6m] [0 : 4m ; 0 : 6m]and x 2 ;y 2 2 [0 : 4m ; 0 : 6m] [2mm ; 0 : 202m]. 48
PAGE 60
FigureV.14: Syntheticforwarddataofneareldscattereddatacomputedby FEKO. 49
PAGE 61
FigureV.15: VolumetricInversionwithTikhonovRegularization.OptimalInverted CurrentDensity j M inv j foundat x =2mm ;y =2mm ;z =0 : 14429m. ThevolumetricinversionregularizedwiththeTikhonovpenaltyisperformedas depictedinFig.V.15.Theoptimalcurrentdensitywaslocatedat x =2mm ;y = 2mm ;z =0 : 14429m.Theinversionperformedaswellasthe L 1 regularizedinversion inFig.V.12.Theinversionlocalizestheoptimal M inv toacellattheedgeofthe imaging,demonstratingpoorlocalizationinthe xy plane.However,the z component oftheoptimal M inv waslocalizedtoaheightcorrespondingtothetwoPECcubesof 0 : 002m
PAGE 62
FigureV.16: VolumetricInversionwith L 1 Regularization.OptimalInvertedCurrentDensity j M inv j foundat x =2mm ;y =2mm ;z =2mm. Thevolumetricinversionregularizedwiththe L 1 penaltyisperformedasdepicted inFig.V.16.Theoptimalcurrentdensitywaslocatedat x =2mm ;y =2mm ;z = 2mm.Thislocalizedtotheedgeoftheimagingdomain.TheCVcurvethatoptimizesthisinversionproblemsharethesamebehaviorasthepreviousvolumetric L 1 regularizedresults.Theabsenceoflocalminimasuggestthatthisregularization methodisnotsuitabletothisproblem. 51
PAGE 63
FigureV.17: AMRInversionwithTikhonovRegularization.OptimalInverted CurrentDensity j M inv j foundat x =34 : 0643mm ;y =34 : 0643mm ;z =34 : 0643mm. TheTikhonovregularizedinversionwasperformedwiththeAMRalgorithmas depictedinFig.V.17.After5iterations,theoptimalinvertedcurrentdensitywas locatedat x =34 : 0643mm ;y =34 : 0643mm ;z =34 : 0643mm.Theseresultsagain demonstrateapoorlocalizationoftheoptimal M inv towardsthebottomoftheimagingdomain. 52
PAGE 64
FigureV.18: AMRInversionwith L 1 Regularization.OptimalInvertedCurrent Density j M inv j foundatx=34.0643mm,y=34.0643mm,z=34.0643mm. The L 1 regularizedinversionwasperformedwiththeAMRalgorithmasdepicted inFig.V.18.After5iterations,theoptimalinvertedcurrentdensitywaslocatedat x =34 : 0643mm ;y =34 : 0643mm ;z =34 : 0643mm.Theseresultsagaindemonstrate apoorlocalizationoftheoptimal M inv towardsthebottomoftheimagingdomain. 53
PAGE 65
FigureV.19: ConstrainedVolumetricInversionwithTikhonovRegularization.OptimalInvertedCurrentDensity j M inv j foundat x =0 : 4m ;y =3mm ;z =0 : 145m. ThispointliesalongoneoftheedgesofCube2. ThevolumetricinversionregularizedwiththeTikhonovpenaltyisperformedas depictedinFig.V.19.Thevolumetricmeshisconstructedwiththeinclusionofthe physicalconstraintsdeterminedfromFig.V.14suchthat M inv isnonzeroforgrid cellslocatedonthesubdomainscorrespondingtothecontoursinthatplot.These cellscorrespondtothetwosubdomains x 1 ;y 1 2 [0 : 4 ; 0 : 6]and x 1 2 [0 : 4 ; 0 : 6] [ y 2 2 [0 : 002 ; 0 : 202].Theoptimal M inv waslocalizedtoacellontheedgeofoneofthe PECcubesat x =0 : 4m ;y =3mm ;z =0 : 145m.Thistechniqueappearstobe promisinginthatitheconstraineddomaineectivelyimprovestheresolutionof thevolumetricgridandiithe M inv waslocalizedtothesurfacethePECcubes ratherthanafreespaceregioninsidetheimagingdomain.Asimilarnumberofcells wasusedtoconstructeachsubdomainwhichachievedabetterresolutioncompared totheunconstrainedproblem.Theuniformimaginggridwasgeneratedforeachof thetwosubdomainswith6 6 8cells. Themagnitudeoftheinvertedcurrentdensity M inv acrossallgridcellsisdepicted 54
PAGE 66
inFig.V.20.TheuniformfeaturesimposedbytheTikhonovpenaltyareapparent inthisplotof j M inv j .Maximaof j M inv j correspondtocellsontheedgeofthePEC cubes. FigureV.20: CurrentDensitymagnitudeofconstrainedvolumetricinversionwith TikhonovRegularization.GridIndexofmaximumcorrespondstooptimalinverted currentdensity j M inv j . FigureV.21: ConstrainedVolumetricInversionwith L 1 Regularization.Optimal InvertedCurrentDensity j M inv j foundat x =0 : 4m ;y =0 : 4m ;z =3mm.This pointliesalongoneoftheedgesofCube1. 55
PAGE 67
Thevolumetricinversionregularizedwiththe L 1 penaltyisperformedasdepicted inFig.V.21.ThevolumetricmeshisconstructedidenticallytothemeshofFig.V.19 withtheinclusionofthephysicalconstraintsdeterminedfromFig.V.14.Theoptimal M inv waslocalizedtoacellontheedgeofoneofthePECcubesat x =0 : 4m ;y = 0 : 4m ;z =3mm,invertingtooneofthecornersofthecentralPECcube.Again, invertinganoptimal M inv tothesurfaceofthePECcubedoessuggestthatthis methodofconstrainingtheimagingdomainiscrucialtoimagingtheinteriorofthe shell. Themagnitudeoftheinvertedcurrentdensity M inv acrossallgridcellsisdepicted inFig.V.22.Thesparsefeaturesimposedbythe L 1 penaltyareapparentinthis plotof j M inv j .Onlyasinglestrongresponseisevidentinthisplotwhileallother responsesareexceedinglysmall.The L 1 penaltysucceededininvertingthisdataset. Unlikethepreviousapplicationsofthispenalty,theCVcurvedidpossessalocal minimum =7 : 4989 10 )]TJ/F17 7.9701 Tf 6.586 0 Td [(4 inthedomainconsidered )]TJ/F17 7.9701 Tf 6.586 0 Td [(5 1,suggesting thattheseconstraintsalsoimprovethebehavioroftheregularizationprocess. FigureV.22: CurrentDensitymagnitudeofconstrainedvolumetricinversionwith L 1 Regularization.GridIndexofmaximumcorrespondstooptimalinvertedcurrent density j M inv j .Theeectsofthe L 1 penaltyareclearlydemonstratedinthesesparse inversionresults. 56
PAGE 68
CHAPTERVI FUTUREWORK Theinversionalgorithmmustbefurtherrenedtoinordertoproducemeaningful andaccuratereconstructionsoftheconsideredclassesofneareldimagingproblems. Directapplicationsoftheseclassicalinversiontechniquesareinsucientinsolving theillposedneareldimagingproblem.TheresultsshowninFigs.V.19V.22indicatethatthesetechniquesbenetgreatlyfromtheinclusionofphysicalconstraints. Furtheranalysisoftheforwardproblemsarenecessarytondessentialphysicalconstraintstooptimallysolvingtheinverseproblem. TheadaptivemeshrenementAMRalgorithmdoesindeedreducethecomputationalcosts,sinceitnecessitatesthecomputationoffewerdipolescomparedtothe volumetricinversion,butitsresultsinaidingtheinversionprocessareinconclusive. Itisdiculttodeterminewhethertheinvertedsolutiondivergesfromtheforward solutionafterafewiterationsduetotheselectioncriteriainthegridrenementprocessorduetothegeneralillposednessoftheinverseproblem.Aproposedsolution involvesturningthisstructuredmeshrenementprocessintoastochasticone.The imaginggridwouldbeinitializedwithanitequantityofgridcellstobedetermined fromastatisticaldistributionintheglobaldomain.Subsequentiterationsrenethe selectedsubdomainusingthisstatisticaldistributionuntiltheconvergencecriteria aremet.ThisstochasticAMRprocedureisrepeatedmanytimesandtheacquired resultsareaveragedtocomputetheinvertedcurrentdensity.Whilethisremains untested,themotivationistoapplyMonteCarloliketechniquestoaidinsolving thisdeterministicinverseproblem. Last,theinaccuraciesinthereconstructioncouldbeattributedtothechoiceof regularizationintheinverseproblem.Simply,the L 1 penaltyortheTikhonovpenalty maynotbesuitabletotheseclassesofproblems.Onealterationwouldbeanelastic 57
PAGE 69
netpenalty[33]whichcombinestheeectsofthesetwopenaltyterms.However,the problemofrequiringthecorrectchoiceofthepenaltytermscanlikelybeavoided entirelyusingmachinelearningmethodstosolvetheinverseproblem[1]. Whileallthesepotentialpathstoimprovingtheaccuracyofthereconstruction oftheneareldelectromagneticimagingproblemarepromisingcandidates,these arenotguaranteedtoaccuratelyandmeaningfullysolvethesespecicinverseproblems.Eithermoresophisticatedtechniquesarenecessaryorthattheseclassesof electromagneticforwardproblemsmustbetreatedindividually. 58
PAGE 70
CHAPTERVII CONCLUSION Presentedhereisaninversionalgorithmforsolvingthegenerallowfrequency neareldelectromagneticinverseproblem.Theelectromagneticscatteringproblem istreatedsimplyasalinearsystemwithaforwardoperatordeterminedbytheGreen's functionresponseoftheproblem G B andthemeasuredelddatainregionsexternal totheimagingregion B sc .TheforwardproblemandthenumericalGreen'sfunction responsesarecomputedwiththeaidofthecommercialelectromagneticmethodof momentssolverFEKO.Theunknownscatterersinthisregionareidenticallyrepresentedbyequivalentvolumetriccurrentdensities J eq withunknownweights.Dueto theillposednessoftheinverseproblemintheLFregime,itisexceedinglydicult tosolvefortheseunknownmodelparameters.Thiscanbealleviatedwiththeinclusionofregularizationandphysicalconstraintsintotheinverseproblem.Various regularizedinversiontechniquesweretestedwiththisalgorithmonthethreecases considered:aPECsphereinfreespace,aPECsphereinaconductivehalfspace, andPECcubesshieldedbyahighlyconductivethinshell.Directapplicationsofthe regularizedinverseproblemonthesethreeconsideredcasesdonotproduceconclusive ormeaningfulresults.However,theintroductionofphysicalconstraintsintotheinverseproblemconsiderablyimprovedthereconstructionsincestrongresponseswere invertedtothesurfacesoftheinternalPECcubes,thoughitdoesnotconclusively reconstructthegeometryofthescatterers.Furtherresearchmustbedoneinorder tooptimizethisinversionalgorithm. 59
PAGE 71
REFERENCES [1]JonasAdlerandOzaOktem.Solvingillposedinverseproblemsusingiterative deepneuralnetworks. InverseProblems ,33:1{24,2017. [2]AhmedI.AlShamma'a,AndrewShaw,andSaherSaman.PropagationofelectromagneticwavesatMHzfrequenciesthroughseawater. IEEETransactionson AntennasandPropagation ,52:2843{2849,2004. [3]AltairEngineeringInc.AltairFEKO,2018. [4]A.P.Annan.ELECTROMAGNETICPRINCIPLESOFGROUNDPENETRATINGRADAR.InHarryM.Jol,editor, GroundPenetratingRadarTheory andApplications ,pages7{8.ElsevierB.V.,2009. [5]SylvainBailletandLineGarnero.ABayesianApproachtoIntroducing AnatomoFunctionalPriorsintheEEG/MEGInverseProblem. IEEETRANSACTIONSONBIOMEDICALENGINEERING ,44:374{385,1997. [6]GlennD.Boreman. ModulationTransferFunctioninOpticalandElectroOptical Systems .SPIEPress,2001. [7]StephenBoydandLievenVandenberghe. ConvexOptimization .Cambridge UniversityPress,1stedition,2004. [8]LeonardoCarrerandLorenzoBruzzone.Solvingforambiguitiesinradargeophysicalexplorationofplanetarybodiesbymimickingbatsecholocation. Nature Communications ,8,2017. [9]WengChoChew. WavesandFieldsinInhomogeneousMedia .IEEEPress,1995. [10]JanC.deMunck,BobW.vanDijk,andHenkSpekreijse.MathematicalDipoles areAdequatetoDescribeRealisticGeneratorsofHumanBrainActivity. IEEE TransactionsonBiomedicalEngineering ,35:960{966,1988. [11]GregoireDerveaux,GeorgePapanicolaou,andChrysoulaTsogka.Resolution anddenoisinginneareldimaging. InverseProblems ,22:1437{1456,2006. [12]K.L.Graf,N.G.Lehtinen,M.Spasojevic,M.B.Cohen,R.A.Marshall,and U.S.Inan.Analysisofexperimentallyvalidatedtransionosphericattenuation estimatesofVLFsignals. JournalofGeophysicalResearch:SpacePhysics , 118:2708{2720,2013. [13]MichaelGrantandStephenBoyd.CVX:Matlabsoftwarefordisciplinedconvex programming,2013. [14]RobertaGrech,TraceyCassar,JosephMuscat,KennethPCamilleri,SimonG Fabri,MichalisZervakis,PetrosXanthopoulos,VangelisSakkalis,andBartVanrumste.ReviewonsolvingtheinverseprobleminEEGsourceanalysis. Journal ofNeuroEngineeringandRehabilitation ,5:1{33,2008. 60
PAGE 72
[15]DavidJ.Griths. IntroductiontoElectrodynamics .PrenticeHall,Inc.,3rd edition,1999. [16]RobertE.Grimm,BarryBerdanier,RobertWarden,JamesHarrer,Raymond Demara,JamesPfeier,andRichardBlohm.Atimedomainelectromagnetic sounderfordetectionandcharacterizationofgroundwateronMars. Planetary andSpaceScience ,57:1268{1281,2009. [17]P.C.Hansen.RegularizationToolsVersion4.0forMatlab7.3,2007. [18]PerChristianHansenandDianneProstO'Leary.TheUseoftheLCurve intheRegularizationofDiscreteIllPosedProblems. SIAMJ.Sci.Comput. , 14:1487{1504,1993. [19]VijayHarid.Resolutionlimitsofneareldelectromagneticimaging.In 2018InternationalAppliedComputationalElectromagneticsSocietySymposiuminDenver,ACESDenver2018 ,number2,pages1{2.AppliedComputationalElectromagneticsSocietyACES,2018. [20]VijayHarid,MarkGolkowski,StephenGedney,RonaldRorrer,PooryaHosseini,MorrisCohen,NathanOpalinski,andSarahPatch.ModelingNearField MagneticShieldingbyConductiveEnclosuresatVeryLowFrequencies.2019. [21]NicholasJ.Higham. AccuracyandStabilityofNumericalAlgorithms .Society forIndustrialandAppliedMathematics,2ndedition,2002. [22]GarethJames,DanielaWitten,TrevorHastie,andRobertTibshirani. AnIntroductiontoStatisticalLearning:withApplicationsinR .2013. [23]JianMingJin. TheoryandComputationofElectromagneticFields .JohnWiley &Sons,Incorporated,2ndedition,2015. [24]DanielLichtblauandEricW.Weisstein.ConditionNumber. [25]Yuri AlvarezLopezandJose AngelMartnezLorenzo.Compressedsensingtechniquesappliedtoultrasonicimagingofcargocontainers. SensorsSwitzerland , 17:1{13,2017. [26]LorenzoLoMonte,DaniloErricolo,FrancescoSoldovieri,andMichaelC.Wicks. Radiofrequencytomographyfortunneldetection. IEEETransactionsonGeoscienceandRemoteSensing ,48PART1:1128{1137,2010. [27]IgorR.ShafarevichandAlexyO.Remizov. LinearAlgebraandGeometry . Springer,2013. [28]DeraldGSmithandHarryMJol.Groundpenetratingradar:antennafrequenciesandmaximumprobabledepthsofpenetrationinQuaternarysediments. JournalofAppliedGeophysics ,333:93{100,1995. 61
PAGE 73
[29]W.M.Telford,W.F.King,andA.Becker.VLFMAPPINGOFGEOLOGICAL STRUCTURE.Technicalreport,Energy,MinesandResourcesCanada,1977. [30]RobertTibshirani.RegressionShrinkageandSelectionviatheLasso. Journal oftheRoyalStatisticalSociety ,58:267{288,1996. [31]A.N.Tikhonov,A.Goncharsky,V.V.Stepanov,andA.G.Yagola. Numerical MethodsfortheSolutionofIllPosedProblems .SpringerNetherlands,1edition, 1995. [32]StuartM.Wentworth. FundamentalsofElectromagneticswithEngineeringApplications .Wiley,2006. [33]HuiZouandTrevorHastie.RegularizationandVariableSelectionviatheElastic Net. JournaloftheRoyalStatisticalSociety.SeriesBStatisticalMethodology , 67:301{320,2005. 62

