Citation

## Material Information

Title:
Low-frequency imaging algorithm for inverting arbitrary electromagnetic systems
Creator:
Todorovski, Dalibor J.
Place of Publication:
Denver, CO
Publisher:
Publication Date:
Language:
English

## Thesis/Dissertation Information

Degree:
Master's ( Master of science)
Degree Grantor:
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:
Holding Location:
Auraria Library
Rights Management:
Copyright Dalibor J. Todorovski. Permission granted to University of Colorado Denver to digitize and display this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.

Full Text
LOW-FREQUENCY 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)
Low-Frequency 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 low-frequency 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 highly-conductive half-space 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 computationally-efficient 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

CHAPTER
I. MOTIVATION........................................................... 1
II. THEORY OF ELECTROMAGNETIC SCATTERING................................ 4
II. 1 Equivalent Dipole Formulation................................. 4
11.2 Scattering in Free-Space....................................... 6
11.3 Scattering by a PEC Sphere Buried in a Homogeneous Half-Space . 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 L-Curve Method......................................... 24
111.4.2 Cross-Validation Method ............................... 25
111.5 Imaging Resolution Limits ................................... 26
IV. ALGORITHM DESCRIPTION............................................. 28
IV. 1 High-Level 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.2 Inversion results for Homogeneous Half-Space 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 free-space illuminated by a TE-polarized plane wave. In the long-wavelength limit, the fields scattered by this can be approximated by the fields of a â€”z-oriented ideal magnetic dipole.......................................................................... 7
11.2 Comparison of analytical scattered magnetic fields due to ideal dipole
against MoM-SIE 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 MoM-SIE solution of Js of PEC sphere buried in homogeneous half-space
(a = 1 S/m) illuminated by TE-polarized 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 half-space. 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 free-space 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 log-log plot of L-Curve 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 5-fold 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 Unite-conductivity 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 high-level 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 free-space handled in FEKO. The unknown PEC scatterers are placed in this model with the known problem boundaries. The scattered near-held data are requested at Nr points..................................................... 31
IV. 5 2D Cross-section 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 near-held 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 near-held 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 near-held requests below bottom face................................ 48
V.14 Synthetic forward data of near-held 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 man-made underground structures [26], non-destructive 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 highly-conductive media typically can not be located with the aid of traditional electromagnetic imaging methods. For instance, ground-penetrating (GPR) systems (10MHz < / < 1GHz) [4] typically image in soils at depths of up to 100 meters [28], but highly-conductive 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, high-frequency (> 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 low-frequency (ELF/VLF) electromagnetic signals (/ <30 kHz) with larger skin-depths 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 trade-off with imaging resolution. Conventional imaging techniques such as GPR, X-ray tomography, or CT scans rely on the far-held (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 ray-like 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 so-called near-held region and the corresponding quasi-static 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 quasi-static magnetic fields in the near-held 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 ill-posedness 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 EEG-like techniques is to infer the magnitude and location of these currents in a non-invasive manner. Thus, the
(II-1)
(11*2)
4

EEG imaging methodology is mathematically identical to the general near-held 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.
(II-6)
5

fm = 6{\r-rg\) (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 over-determined and inversion of (II.9) must be solved in a least-squares 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 Free-Space
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 TE-polarized 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 free-space illuminated by a TE-polarized plane wave. In the long-wavelength limit, the fields scattered by this can be approximated by the fields of a â€”^-oriented ideal magnetic dipole.
In the long-wavelength 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 â€” 0-directed. Computation of the scattering due these surface currents can be simplified by considering the fields radiated by an equivalent â€” i-directed 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 near-held as (11-13).
m=â€”2ir a3#' (11.13)
These analytical expressions for Hsc were compared to method of moments -surface integral equation (MoM-SIE) solutions obtained in the commercial electromagnetic solver FEKO along a cut-line 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 MoM-SIE 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 near-held 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 Half-Space
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 half-space with parameters (Â£,p., Figure II.3: Configuration of PEC sphere of radius a embedded in homogeneous half-space (Region 2) with conductivity a illuminated by an incident plane wave.
For simplification of the analysis, the plane wave is assumed to be a TE-polarized 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 half-space (Region 2) and free-space (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
(II-14)
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
(H-17)
rpTM
2 tk0
tko + tok2
(II-18)
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 half-space. 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 non-uniform induced Js over the â€”0-direction. The attenuation of H-2y(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: MoM-SIE solution of Js of PEC sphere buried in homogeneous halfspace ( Since Js is not uniform across the â€”0-direction, 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 half-space 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 p-direction (p = xcos(f> + ysliuf>) and a plane-wave propagating in the z-
of the Sommerfeld approach is that each cylindrical wave can be reflected from the air-ground 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 k-space (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\ _ ftTEe-jk2,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 half-space 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 half-space.
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 thin-shell 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 ill-posed. 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 well-posed, 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. A-1 is defined for all B and is continuous.
The electromagnetic scattering problem typically fails both of these conditions and is Hadamard ill-posed.
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(A|6), where (A\b) is the augmented matrix. By
Rouche-Capelli 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 non-singular 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 free-space.
17

Figure III. 1: Figure from [19]. A collection of N ideal dipoles in free-space 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 Free-Space 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 well-posedness in this frequency regime. While the inverse problem is ill-posed 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||:r|E (HI-4)
Equation (HI.4) represents the general regularized minimization problem, where A is the regularization parameter. The second term in (HI.4) represents the general p-norm 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 p-norm in the literature are the p = 2 Tikhonov regularization [14] andp = 1 LASSO regularization (or regression) [30]. These p-norm 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 L2-norm minimization problem (Tikhonov Regularization problem).
min ||6 â€” Ax\\\ + A||rc||| (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 â€” Ax||2 + A||x||2 is differentiable for all x in the problem space.
Â£ = || b â€” Ax \\l + A||x||2 = xt(AtA + A I)x â€” bT Ax + bTb (HI-7)
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 Li-norm 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 Li-norm 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 closed-form solution for the orthonormal design case (XTX = I) can be obtained by considering sub-domain 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 sub-domain 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 sub-domains 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Â°)[2|ftÂ«|-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 Kuhn-Tucker conditions [30]. The inequality constraint JT \(3j\ < A is expressed as 8f/3 < A, where 8i for i = 1,2,..., 2â„¢ are the n-tuples 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 Kuhn-Tucker 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 L-Curve method and Cross-Validation (CV) method.
23

III.4.1 L-Curve Method
The L-Curve method, named for the characteristic shape of the log-log plot as in Fig. III.4, attempts to minimize the trade-off 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 over-smoothing. 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 log-log plot of L-Curve 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 L-curve 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 Cross-Validation Method
The cross-validation method is a statistical technique of estimating the mean squared error (MSE) of the solution without explicit knowledge of the true solution. Only k-fold 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 k-fold CV involves dividing the data set b into fc-equally 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 5-fold 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 far-held (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 finite-conductivity 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 far-held approximation [11], However, this near-held 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 High-Level Overview of Algorithm
The algorithms described attempt to solve the inverse problem for the low-frequency near-held 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 L-Curve method, or solved with the aid of the CVX MATLAB toolbox [13] for the Li-norm 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 high-level 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 half-space 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 free-space as well as the case of a conductive half-space (half-space 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 (SIE-MoM) 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 near-held 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 free-space handled in FEKO. The unknown PEC scatterers are placed in this model with the known problem boundaries. The scattered near-held 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 half-space 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 (half-space 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 regularized-inversion, or L2 Regularization, attempts to minimize the 2-norm 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 L-Curve 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 L-curve 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 1-norm 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.
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 1-5% 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 Cross-section 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.
This inversion algorithm was validated with the scattering by a PEC sphere of radius a = 20cm located in free-space 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 xy-plane (similar to the incident fields described in Sec. II. 3). A factor of ej27r^, with / = 25kHz, is assumed. The near-held 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 free-space limited to the sub-domain 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

Hy-Field Magnitude [uA/m]: Free-Space 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 near-held scattered data computed by FEKO (domain limited to (x, y) G [â€”2m, 2m]).
Tikhonov Inversion w/ Regularization: A = 5.4163e-08
^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 L-Curve 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 over-regularized 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 = 1e-05
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 cross-validation 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 (10-5 < 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 free-space problem.
-1.4 -1.6--1.8--2-
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.216e-08
â€”^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 sub-domain.
L1 Inversion w/ Regularization: A = 1e-05
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 sub-domain 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 free-space problem are likely attributable to the choice of the measurement plane. The lack of variation in the z-direction 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 Half-Space Inversion
The forward problem was constructed identically to the free-space 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 free-space problem, though the measured fields are twice as large due to reflection from the interface between the free-space region and the finite-conductivity half-space 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 free-space problem, it is not clear what the geometry and depth of the PEC scatterer are from direct examination of the forward data.
42

Hy-Field 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 near-held scattered data computed by FEKO (domain limited to (x,y) E [â€”2m, 2m]).
Â£ -1
-2;
2
Tikhonov Inversion w/ Regularization: A = 4.3545e-08
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 free-space 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 free-space and conductive halfspace problem are not clear.
L1 Inversion w/ Regularization: A = 1e-05
^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 free-space problem. This again suggests that the L\ penalty is simply not suitable to the problem of scattering of a PEC sphere in free-space or the scattering of a PEC sphere buried in a homogeneous half-space.
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 free-space 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 ill-conditioned and a least squares solution serves as a poor ht for the data set.
45

Tikhonov Inversion w/ Regularization: A = 1.0194e-08
> 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 iry-plane 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 rry-plane could be attained through the inclusion of physical constraints in the inversion problem.
46

L1 Inversion w/ Regularization: A = 1e-05
| 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 iry-plane is poor as in the results of the Tikhonov-regularized 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 near-held requests below bottom face.
The near-held 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 sub-domain is given by (xi,yi) E [0.4m, 0.6m] x [0.4m, 0.6m] and (:r-2, t/2) G [0.4m, 0.6m] x [2mm, 0.202m],
48

Hz-Field 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 near-held 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 rry-plane. 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 = 1e-05
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.3518e-05
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.6102e-05
â€”Â»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 sub-domains corresponding to the contours in that plot. These cells correspond to the two sub-domains (Vi,yi) G [0.4, 0.6] and x\ G [0.4, 0.6] U y-2 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 free-space region inside the imaging domain. A similar number of cells was used to construct each sub-domain which achieved a better resolution compared to the unconstrained problem. The uniform imaging grid was generated for each of the two sub-domains 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
16--------------1-------------1-------------1------------
^,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 10-4 in the domain considered (10-5 < 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'------------1-----------1-----------1-----------1------------1-----------
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 near-held imaging problems. Direct applications of these classical inversion techniques are insufficient in solving the ill-posed near-held imaging problem. The results shown in Figs. V.19-V.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 ill-posedness 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 sub-domain 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 Monte-Carlo-like 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 near-field 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 low-frequency near-held 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 ill-posedness 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 free-space, a PEC sphere in a conductive half-space, and PEC cubes shielded by a highly-conductive 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 ill-posed inverse problems using iterative deep neural networks. Inverse Problems, 33(12): 1 24, 2017.
[2] Ahmed I. Al-Shammaâ€™a, Andrew Shaw, and Saher Saman. Propagation of electromagnetic waves at MHz frequencies through seawater. IEEE Transactions on Antennas and Propagation, 52(ll):2843-2849, 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 7-8. Elsevier B.V., 2009.
[5] Sylvain Baillet and Line Garnero. A Bayesian Approach to Introducing Anatomo-Functional Priors in the EEG/MEG Inverse Problem. IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, 44(5):374-385, 1997.
[6] Glenn D. Boreman. Modulation Transfer Function in Optical and Electro-Optical 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 near-held 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 trans-ionospheric attenuation estimates of VLF signals. Journal of Geophysical Research: Space Physics, 118(5):2708-2720, 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 Van-rumste. 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. Prentice-Hall, Inc., 3rd edition, 1999.
[16] Robert E. Grimm, Barry Berdanier, Robert Warden, James Harrer, Raymond Demara, James Pfeiffer, and Richard Blohm. A time-domain 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 L-Curve in the Regularization of Discrete Ill-Posed Problems. SIAM J. Sci. Comput., 14(6): 1487â€”1504, 1993.
[19] Vijay Harid. Resolution limits of near-field electromagnetic imaging. In 2018 International Applied Computational Electromagnetics Society Symposium in Denver, ACES-Denver 2018, number 2, pages 1-2. Applied Computational Electromagnetics Society (ACES), 2018.
[20] Vijay Harid, Mark Golkowski, Stephen Gedney, Ronald Rorrer, Poorya Hos-seini, Morris Cohen, Nathan Opalinski, and Sarah Patch. Modeling Near-Field 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(1-3):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 III-Posed 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):301-320, 2005.
62

Full Text

PAGE 1

LOW-FREQUENCYIMAGINGALGORITHMFORINVERTINGARBITRARY 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 Low-FrequencyImagingAlgorithmforInvertingArbitraryElectromagneticSystems ThesisdirectedbyAssistantProfessorVijayHarid ABSTRACT MethodsandsoftwareweredevelopedinMATLABforthepurposesofinverting arbitrarylow-frequencyelectromagneticscatteringproblemsandlocatingtheequivalentcurrentdensitiespresentintheforwardproblem.Whiletheclassesofproblems consideredforthismethodarethosewhereforwardscatterersareconcealedwithin ahighly-conductivehalf-spaceorametallicshell,themethodscanbegeneralizedto otherproblemsduetotherobustnessoftheelectromagneticsolverused.Synthetic electromagneticelddatawascreatedfortheconsideredcaseswiththeuseofthe commercialelectromagneticsolverFEKO.Similarly,theGreen'sfunctionresponse wasconstructednumericallyinthissolverusingtheknownphysicalboundariesof theproblemanditsmaterialproperties.Withthesecomputedquantities,alinear relationshipwasformedbetweentheforwardelectromagneticdataandthenumerical Green'sfunctionNGFthatcouldbeinvertedwiththeaidofvariousregularization techniques.ResultsarepresentedforinversionsusingbothTikhonovregularization andL 1 regularizationonafullvolumetricreconstructionoftheimagingspaceaswell asamorecomputationally-ecientiterativereconstruction. Theformandcontentofthisabstractareapproved.Irecommenditspublication. Approved:VijayHarid iii

PAGE 4

ACKNOWLEDGMENTS ThankyoutoalloftheprofessorsattheUniversityofColorado-Denver.Youhave alltaughtmeinvaluableknowledgeandskillsasanengineerandresearcherthatI willcarryforwardwithme.AndIwouldliketoextendmygratitudetoDr.Vijay Harid.Yourguidanceandencouragementoverthesepasttwoyearshavebeencritical inhelpingmendmypathforward. iv

PAGE 5

TABLEOFCONTENTS CHAPTER I.MOTIVATION.................................1 II.THEORYOFELECTROMAGNETICSCATTERING...........4 II.1EquivalentDipoleFormulation.....................4 II.2ScatteringinFree-Space.........................6 II.3ScatteringbyaPECSphereBuriedinaHomogeneousHalf-Space.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.1L-CurveMethod.........................24 III.4.2Cross-ValidationMethod....................25 III.5ImagingResolutionLimits.......................26 IV.ALGORITHMDESCRIPTION........................28 IV.1High-LevelOverviewofAlgorithm...................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.1InversionresultsforFree-SpaceInversion...............36 V.2InversionresultsforHomogeneousHalf-SpaceInversion.......42 V.3InversionresultsforConductiveCubicShellInversion........47 VI.FUTUREWORK................................57 VII.CONCLUSION.................................59 REFERENCES...................................60 vi

PAGE 7

FIGURES II.1ScatteringofPECsphereofradius a infree-spaceilluminatedbyaTEpolarizedplanewave.Inthelong-wavelengthlimit,theeldsscattered bythiscanbeapproximatedbytheeldsofa )]TJ/F15 11.9552 Tf 10.007 0 Td [(^ z -orientedidealmagnetic dipole......................................7 II.2Comparisonofanalyticalscatteredmagneticeldsduetoidealdipole againstMoM-SIEsolutionsofspherescatteringandequivalentmagnetic dipole......................................8 II.3CongurationofPECsphereofradius a embeddedinhomogeneoushalfspaceRegion 2 withconductivity illuminatedbyanincidentplane wave......................................9 II.4MoM-SIEsolutionof ~ J s ofPECsphereburiedinhomogeneoushalf-space =1S/milluminatedbyTE-polarizedplanewave f =1kHz.A larger ~ J s isinducedatthetopofsphererelativetothebottomdueto attenuationoverthelengthofthesphereinsidethehalf-space.......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 idealdipolesinfree-spaceuniformly 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.4Characteristiclog-logplotofL-CurveMethodforTikhonovRegularization.Thecornerorclosestpointtotheoriginrepresentstheoptimal parameter ..................................24 III.5Equipartitionofdatasetfor5-foldCV.Thisprovidesanestimateofthe MSEofthesolutionvector x thatcanbeusedtooptimizetheregularizationparameter ................................25 III.6MinimumSpatialResolutionforsignalpropagatinginvarioushomogeneousnite-conductivitygrounds.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.1Flowchartdescribinghigh-leveldesignofsoftwareforbothvolumetric 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.4GenerationofForwardProblemofPECSpherelocatedinfree-spacehandledinFEKO.TheunknownPECscatterersareplacedinthismodel withtheknownproblemboundaries.Thescatterednear-elddataare requestedat N r points.............................31 IV.52DCross-sectionofinitialcoarsegridleft.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.1Syntheticforwarddataofnear-eldscattereddatacomputedbyFEKO 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.7Syntheticforwarddataofnear-eldscattereddatacomputedbyFEKO 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 andnear-eldrequestsbelowbottomface..................48 V.14Syntheticforwarddataofnear-eldscattereddatacomputedbyFEKO.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],detectionofnaturalandman-made undergroundstructures[26],non-destructiveevaluationforvehiclesandcontainerinspections[25],underwaterdetection[2,16],andseveralotheractiveareasofresearch. Thediculty,however,isthatobjectsofinterestconcealedbehindhighly-conductive mediatypicallycannotbelocatedwiththeaidoftraditionalelectromagneticimaging methods.Forinstance,ground-penetratingGPRsystemsMHz 30kHzelectromagneticradiationcannotpenetratethrough metalliccontainersofeven1 )]TJ/F15 11.9552 Tf 12.968 0 Td [(2mmthicknesswithoutundergoingextremelosses > 300dB.Fortunately,thedrawbackofhighlossescanbeovercomewiththeuseof extremelyandverylow-frequencyELF/VLFelectromagneticsignals f< 30kHz withlargerskin-depths 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,thereisaninherenttrade-owithimaging resolution.ConventionalimagingtechniquessuchasGPR,X-raytomography,orCT scansrelyonthefar-eldraypropertiesofelectromagneticradiationtoreconstruct 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 whenusingimagingalgorithmsthatrelyontheray-likefeaturesofelectromagnetic waves.Tocircumventthislimit,itisnecessarytodevelopnovelstrategiesthatdonot relyontypicalelectromagneticwavefeaturesforimaging.Imagingsmallobjectsin theELF/VLFrangewilloccurintheso-callednear-eldregionandthecorresponding quasi-staticnatureofthesesignalsmustbeusedtoreconstructanimage.Physically, incidentELF/VLFsignalscaninduceloopsofelectriccurrentonanobjectofinterest whichwillinturngeneratequasi-staticmagneticeldsinthenear-eldoftheobject. Thesemagneticeldswillhaveaspatialamplitudeandpolarizationprolethatare specictotheobject'sshapeandsize.Thus,anappropriateELF/VLFimagingalgo2

PAGE 14

rithmcanbedesignedbyleveragingthemappingfromtheobject'sgeometrytothe measuredmagneticeldprole.Evenso,inferringthelocationofhiddenobjectsfrom elddataisinitselfaconsiderablechallengeduetotheill-posednessoftheproblem 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].ThepurposeofEEG-liketechniquesistoinfer themagnitudeandlocationofthesecurrentsinanon-invasivemanner.Thus,the 4

PAGE 16

EEGimagingmethodologyismathematicallyidenticaltothegeneralnear-eldimage 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 ,thislinearsystemisover-determinedandinversionof II.9mustbesolvedinaleast-squaressense.Accuratereconstructionoftheforward problemfromthismeasurementdata E sc ;B sc requiresregularizationandapriori constraintswhicharediscussedinChapterIII. II.2ScatteringinFree-Space Considertheproblemofelectromagneticscatteringbyaperfectlyelectrically conductingPECsphereofradius a infreespace.Thisscattererisilluminatedby aTE-polarizedplanewavewithwavenumber k propagatinginthe+^ y directionas showninFig.II.1. 6

PAGE 18

FigureII.1: ScatteringofPECsphereofradius a infree-spaceilluminatedbya TE-polarizedplanewave.Inthelong-wavelengthlimit,theeldsscatteredbythis canbeapproximatedbytheeldsofa )]TJ/F15 11.9552 Tf 10.007 0 Td [(^ z -orientedidealmagneticdipole. Inthelong-wavelengthlimit 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 inthenear-eldasII.13. ~ m = )]TJ/F15 11.9552 Tf 9.299 0 Td [(2 a 3 ~ H i II.13 Theseanalyticalexpressionsfor H sc werecomparedtomethodofmomentssurfaceintegralequationMoM-SIEsolutionsobtainedinthecommercialelectromagneticsolverFEKOalongacut-line2metersabovethePECsphere z =2mas inFig.II.2. FigureII.2: Comparisonofanalyticalscatteredmagneticeldsduetoidealdipole againstMoM-SIEsolutionsofspherescatteringandequivalentmagneticdipole. Theseanalyticalresultsfor H sc agreeconsiderablywiththeMoMsolutionsof thisscatteringproblem.Thissimpleapproachofrepresentingmagneticscattering byconductiveobjectswithequivalentmagneticdipolesisavalidapproximationin thenear-eldfortheunknownelectromagneticscatterersintheimagingproblem. Sinceamagneticdipoleisidenticaltoaloopofelectricalcurrent,thisapproachis 8

PAGE 20

identicaltosearchingforloopsofcurrentasopposedtoelectricaldipolesthatare inducedonanobject.Thisisaconvenientabstractionthathelpsalleviatesomeof thecomputationalburdenassociatedwithimagereconstruction. II.3ScatteringbyaPECSphereBuriedinaHomogeneousHalf-Space ConsidertheelectromagneticscatteringofaPECsphereilluminatedbyanincidentplanewaveasshowninFig.II.3.ThePECsphereisembeddedwithina homogeneoushalf-spacewithparameters ";; atadepthof d meters.Theregion aboveisfree-space. FigureII.3: CongurationofPECsphereofradius a embeddedinhomogeneous half-spaceRegion 2 withconductivity illuminatedbyanincidentplanewave. Forsimplicationoftheanalysis,theplanewaveisassumedtobeaTE-polarized 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 half-spaceRegion 2 andfree-spaceRegion 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 theimpedancemismatchbetweenthetwomediaattheinterfaceandthenfromattenuationasitpropagatesdeeperintothehalf-space.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 anon-uniforminduced ~ 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: MoM-SIEsolutionof ~ J s ofPECsphereburiedinhomogeneoushalfspace =1S/milluminatedbyTE-polarizedplanewave f =1kHz.Alarger ~ J s isinducedatthetopofsphererelativetothebottomduetoattenuationoverthe lengthofthesphereinsidethehalf-space. 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 .Includingtheeectofaniteconductivityhalf-spacerequiresreformulatingtheGreen'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 andaplane-wavepropagatinginthe^ z directionforallwavenumbersink-space k 2 = k 2 + k 2 z [9].Theprimaryadvantage oftheSommerfeldapproachisthateachcylindricalwavecanbereectedfromthe air-groundinterfaceviathetraditionalFresnelcoecients.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 thehomogeneoushalf-spaceproblemisconsideredinthiswork. 13

PAGE 25

ComputationoftheSommerfeldIntegralsinII.21areingeneralanalytically intractableandthusrequirenumericaltechniques.Typicalasymptoticmethodsdetailedin[9]arehighlyproblemdependentandpotentiallyposecomplicationsfor scatteringproblemswithlittletonoaprioriinformationaboutthescatterer.Inthis work,thisproblemofevaluatingtheSommerfeldIntegralsissimpliedbyutilizingthe commercialelectromagneticsolverFEKOtoquicklycomputeboththeforwardscatteringproblemaswellastogenerateanumericalGreen'sFunctionresponseNGF foranequivalentmagneticmagneticdipoleburiedinaconductivehalf-space. 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. Forthecaseofasphericalthin-shellshield,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],ishighlyill-posed.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.1canbesolvedexactlyifitiswell-posed,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 isHadamardill-posed. 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 idealdipolesinfree-spaceuniformly 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 conditionsforwell-posednessinthisfrequencyregime.Whiletheinverseproblem isill-posedanddiculttosolve,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.Twocommonchoicesofp-normintheliteraturearethe p =2Tikhonov regularization[14]and p =1LASSOregularizationorregression[30].Thesep-norm 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. Aclosed-formsolutionfortheorthonormaldesigncase X T X = I canbeobtainedbyconsideringsub-domainwheretheabsolutevaluefunction 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 Solutionstothemodiedobjectivecanbefoundforeachsub-domain 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.15fortheunionofthesesub-domains 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 theKuhn-Tuckerconditions[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 Kuhn-Tuckerconditions.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.TwotechniquescommonlyappliedintheliteraturearetheL-Curvemethodand Cross-ValidationCVmethod. 23

PAGE 35

III.4.1L-CurveMethod TheL-Curvemethod,namedforthecharacteristicshapeofthelog-logplotasin Fig.III.4,attemptstominimizethetrade-obetweenthegeneralizedregularization 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 dominatedbyregularizationerrorssuchasover-smoothing.Otherwise,if istoo small,thesolutionwillbedominatedbythemodelingerrorsandmeasurementnoise [18]. FigureIII.4: Characteristiclog-logplotofL-CurveMethodforTikhonovRegularization.Thecornerorclosestpointtotheoriginrepresentstheoptimalparameter . Therearetwomethodsoflocatingthiscorner[18]. 1.Findingthecornerisequivalenttondingtheclosestpointonthecurvetothe origin.Thismetricdiersbasedonthechoiceofregularization.ForTikhonov regularization,thismetricis d = + 2 . 2.Thecornerisalsothepointofmaximumcurvature .Thisisindependentof thechoiceofregularizationsince isapurelygeometricquantity.Assuming 24

PAGE 36

theL-curveisasmoothcurveandthat^ =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.2Cross-ValidationMethod Thecross-validationmethodisastatisticaltechniqueofestimatingthemean squarederrorMSEofthesolutionwithoutexplicitknowledgeofthetruesolution. Onlyk-foldCVisconsideredfortheseproblems[22].Thisestimateconvergestothe trueMSEsincethesesamplessharethestatisticalpropertiesofthedataset.The k-foldCVinvolvesdividingthedataset b into k -equallysizedpartitions.Oneofthese partitionsiskeptasaholdoutvalidationsetandtheremaining k )]TJ/F15 11.9552 Tf 12.264 0 Td [(1areusedas trainingsets,where k istypicallytakentobe5or10. FigureIII.5: Equipartitionofdatasetfor5-foldCV.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 operateinthefar-eld ka 1,thisbestresolutionisdeterminedbytheRayleigh limit,givenbyIII.19[11]. R ff =1 : 22 2 III.19 Spatialfeaturesthatarelargerthan R ff areresolvablebythesystem. FigureIII.6: MinimumSpatialResolutionforsignalpropagatinginvarioushomogeneousnite-conductivitygrounds.Thisassumesnoattenuationofthesignal. ItisevidentfromFig.III.6thattheobjectsofinterestwhicharetypicallysizes ontheorderof1marenotresolvableatallbyanimagingsystemoperatinginthefareldassuminganoiselessandlosslesssignal.Theinformationofthesesmallspatial 26

PAGE 38

featuresarecontainedintheevanescentwaveswhichareneglectedinthefar-eld approximation[11].However,thisnear-eldresolutionlimitwilldiergreatlyfor 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. TheresolutionlimitsattainablebyaTikhonov-or L 1 -regularizedsolutionarenot examinedbutareassumedtobeequal,ifnotbetter,thanthoseforanunregularized solutionasin[19].Intheproblemswhereafullvolumetricinversionisconsidered, theresolutionoftheuniformgridisgreaterthanthelimitposedfortheordinary leastsquaressolutionandareassumednottobelimitedbythisresolutionlimit. 27

PAGE 39

CHAPTERIV ALGORITHMDESCRIPTION IV.1High-LevelOverviewofAlgorithm Thealgorithmsdescribedattempttosolvetheinverseproblemforthelowfrequencynear-eldelectromagneticscatteringproblem.Itreliesonthecreationof syntheticforwarddata B data andtheNGF G B inthecommercialelectromagnetic solverFEKO[3].Thisinversionprcessisagnostictothechoiceofelectromagnetic solver.Thissolverwaschosenduetotheavailabilityandadvantagesofthemoments methodMoMelectromagneticsolveroverothercommercialsolvers.TheMoM solverhassimpleandecienttreatmentsofidealdipolesourcesandonlyrequires meshingofthesurfaceelementstosolvethescatteringproblem.Theinverseproblem issolvedeitheranalyticallyfortheTikhonovregularizedproblem,whichisoptimized bytheL-Curvemethod,orsolvedwiththeaidoftheCVXMATLABtoolbox[13] forthe L 1 -normregularizedproblem,optimizedbytheCVmethod.Twotechniques forthegenerationoftheimaginggridarepresented:asimplestructuredvolumetric generationoftheimaginggridandaniterativeandlesscomputationallyintensive generationoftheimaginggridtermedadaptivemeshrenement"AMR. 28

PAGE 40

FigureIV.1: Flowchartdescribinghigh-leveldesignofsoftwareforbothvolumetric meshandadaptivemeshrenementprocedures. IV.2GenerationofSyntheticData TheforwardproblemisfullymodelledinFEKO.Onlythreecasesareconsidered fortheforwardproblem:ithescatteringproblemposedinSec.II.2,iiscattering fromaPECsphereburiedinaconductivehalf-spaceasinSec.II.3,andiiithe scatteringfromPECobjectslocatedinsideinametalliccubicshield.FigureIV.2 showsanexampleforwardsimulationsetupoftwoPECcubeslocatedinsidea1m 1m 1mconductivecubicshellof2mmthicknessandconductivity =10 7 S/m. FigureIV.4showsthesetupofaPECsphereinfree-spaceaswellasthecaseofa conductivehalf-spacehalf-spacenotvisibleinsimulationCADle. 29

PAGE 41

FigureIV.2: ForwardProblemGeometryforMetallicBoxwithtwoPECCubes. Themetallicshellhas =10 7 S/mandthickness=2mm.Receiversareplacedon planesparalleltoeachfaceosetby5mm.Anidealmagneticdipolesourceisplaced overeachfacethatoperatesintheLFregime. TheprocessofgeneratingsyntheticforwarddataintheFEKOsolverisgeneral andsotheinversionprocesscanpotentiallybeextendedtoanyclassofelectromagneticinversescatteringproblems.Thegeometryoftheforwardproblem,which includesboththeunknownscatterersandtheknownproblemboundaries,ismodeled withtheCADtoolsprovidedintheFEKOsolver.ThissolverthensolvestheelectromagneticscatteringproblemusingtheSurfaceIntegralEquationMoMSIE-MoM approach[3].Thesurfacesofthemodeledgeometryareapproximatedbyameshof triangularelementsbyFEKO. 30

PAGE 42

FigureIV.3: Surfacesaremeshedwithmanytriangularelements.ThebasisfunctionsusedtosolvetheSIEaremappedtotheedgesoftheseelementsandusedto computesurfacecurrentdensities. Whileincreasingthenumberofthesemeshelementsimprovestheaccuracyofthe model,itdoessoataincreasedcomputationalcost.Materialpropertiessuchas conductivityandthicknessareappliedtotheapplicablesurfacesofthemodel.Last, thenear-elddataofthescatteredmagneticeldsarerequestedat N r arbitrary pointsplacedoutsidetheproblemspace. FigureIV.4: GenerationofForwardProblemofPECSpherelocatedinfree-space handledinFEKO.TheunknownPECscatterersareplacedinthismodelwiththe knownproblemboundaries.Thescatterednear-elddataarerequestedat N r points. 31

PAGE 43

IV.3ComputationofNumericalGreen'sFunctionNGF ClosedformexpressionsfortheGreen'sfunctionresponseofthehalf-spaceand 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 isausefulmethodtoincludetheeectsofthebackgroundshieldhalf-spaceorshell withouthavingtoexplicitlymodelthebackgroundintheinversionalgorithms. IV.4InversionAlgorithm Threedierenttypesofinversionmethodsareattemptedtosolvefortheunknown modelparameters M inv .Eachattemptstoimposecertainfavorablefeaturesinthe modelspacetoaccuratelyinverttheprobleminquestion. IV.4.1TikhonovRegularization Tikhonovregularized-inversion,or L 2 Regularization,attemptstominimizethe 2-normofthesolutionvector 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 foundthroughtheL-CurvemethodasdetailedinSecIII.4.1.TheRegularization ToolspackageforMATLABisutilized[17]toaidincomputationof .Thistreatsthe discreteL-curveasadierentiablecurvetocomputethepointofmaximumcurvature . 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 Sincethe1-normpenaltyisnotdierentiable,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,thisthresholdwastypicallychosentobe1-5%of thedimensionsoftheglobalmodelspace. Inthisway,eachiterationeectivelyadds7gridcellstothemodelsincetheoriginal gridcellisremovedduringtherenementprocess.Thealgorithmforasingleiteration isdemonstratedinFig.IV.5. FigureIV.5: 2DCross-sectionofinitialcoarsegridleft.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.1InversionresultsforFree-SpaceInversion ThisinversionalgorithmwasvalidatedwiththescatteringbyaPECsphereof radius a =20cmlocatedinfree-spaceat 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.Thenear-elddataissampledalong 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 forwarddataofthescatteringprobleminfree-spacelimitedtothesub-domain 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: Syntheticforwarddataofnear-eldscattereddatacomputedbyFEKO 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 theL-Curvemethodhadselectedsuchasmallvaluefortheregularizationparameter =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 cellsclosesttothemeasurementplaneintheover-regularizedinversionisduetothe 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,sincetheestimatedcross-validationcurvewasmonotonicallyincreasingoverthetestedregularizationparametersforthisdataset.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 techniqueisunsuitableforavolumetricinversionofthefree-spaceproblem. 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 anadequatemethodofselectingtherenementsub-domain. 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 convergenceoftheseresultsislikelyduetotheselectioncriteriaforsub-domainrenement,whichlikelyresultsinabiastowardsthetopoftheimagingdomain.Another likelyexplanationisthatduetotheunsuitabilityofthe L 1 regularizationtothis problemduetotheabsenceofalocalminimumintheCVcurve,theinversion isnobetterthandoinganordinaryleastsquaresinversioni.e.minimizingtheunregularizedminimizationproblemIII.3duetotheabsenceoflocalminimainthe domain. Aproposedtesttoevaluatethesuitabilityoftheregularizationtechniquetoa datasetwouldbetoexaminethebehavioroftheCVcurveas ! 0.IftheCV 41

PAGE 53

curvedecreasesasymptoticallyinthislimit,thentheregularizationmethodisno betterthananordinaryleastsquaressolutionduetotheabsenceoflocalminimain thecurve. Theseinconclusiveinversionresultsforthefree-spaceproblemarelikelyattributabletothechoiceofthemeasurementplane.Thelackofvariationinthe ^ z -directioncoulddirectlyberesponsiblefortheseissues.Theaddeddiversityin measurementscouldbecrucialinresolvingdepthinformation,thoughnoresultsare availableforthisscatteringproblemtoverifytheeectsofdierentchoicesofmeasurementgeometry. V.2InversionresultsforHomogeneousHalf-SpaceInversion Theforwardproblemwasconstructedidenticallytothefree-spaceprobleminSec. V.1.AgroundplanesimulatedwiththeexactSommerfeldintegralswasincludedat z =0with =10mS/m.Theforwarddatasharesthecharacteristicshapeasinthe free-spaceproblem,thoughthemeasuredeldsaretwiceaslargeduetoreection fromtheinterfacebetweenthefree-spaceregionandthenite-conductivityhalf-space at z =0.Similarly,thelocationofthePECscatterercanbedeterminedtobenear theorigin x =0 ;y =0duetothelocationofthemaximuminthemagneticeld data.Asinthefree-spaceproblem,itisnotclearwhatthegeometryanddepthof thePECscattererarefromdirectexaminationoftheforwarddata. 42

PAGE 54

FigureV.7: Syntheticforwarddataofnear-eldscattereddatacomputedbyFEKO 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.Unliketheresultsofthefree-spaceproblem,thiswasbiased towardsthebottomoftheimaginggridwithasimilarregularizationparameter . Thebehaviorofbiasingtowardstheboundariesoftheimaginggridcouldbehighly problemdependent.Testingforoverttinginthisproblembiasesthelocationof theoptimalinversionresultstothebottomoftheimaginggrid.Thereasonforthe disparityinthislocalizationbehaviorbetweenthefree-spaceandconductivehalfspaceproblemarenotclear. 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

thisproblemsharedthesamebehaviorasthefree-spaceproblem.Thisagainsuggests thatthe L 1 penaltyissimplynotsuitabletotheproblemofscatteringofaPECsphere infree-spaceorthescatteringofaPECsphereburiedinahomogeneoushalf-space. 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.Asinthefree-spaceproblem,thecomputedMSEofthe ordinaryleastsquaressolutionforeachgridcellwasidentical,givingtherstelement inthevector M inv .Attemptstotthisproblemwithindividualdipolesarelikely unsuitableforeitherithecellsofuniformgridarenotcorrectlychosenastheydo notcoincidewiththeinteriorofthePECsphereoriitheforwardoperatorfora singledipoleisill-conditionedandaleastsquaressolutionservesasapoortforthe 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 -planeispoorasintheresultsoftheTikhonov-regularizedinversion. 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 topfaceandnear-eldrequestsbelowbottomface. Thenear-eldscattereddataisrequested1mmbelowthebottomface.The outlinesofthePECcubesincontactwiththisfaceclearlyappearinthemeasured elddataasinFig.V.14.Suchinformationisvaluableinestablishingphysical constraintstotheinversionproblemforobjectsnearorincontactwiththemetallic shell.Thiswasusedtoconstraintheimagingdomainsothat M inv isnonzeroonly forgridcellscorrespondingtothisoutline:thissub-domainisgivenby 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: Syntheticforwarddataofnear-eldscattereddatacomputedby 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 cellslocatedonthesub-domainscorrespondingtothecontoursinthatplot.These cellscorrespondtothetwosub-domains 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 ratherthanafree-spaceregioninsidetheimagingdomain.Asimilarnumberofcells wasusedtoconstructeachsub-domainwhichachievedabetterresolutioncompared totheunconstrainedproblem.Theuniformimaginggridwasgeneratedforeachof thetwosub-domainswith6 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 andaccuratereconstructionsoftheconsideredclassesofnear-eldimagingproblems. Directapplicationsoftheseclassicalinversiontechniquesareinsucientinsolving theill-posednear-eldimagingproblem.TheresultsshowninFigs.V.19-V.22indicatethatthesetechniquesbenetgreatlyfromtheinclusionofphysicalconstraints. Furtheranalysisoftheforwardproblemsarenecessarytondessentialphysicalconstraintstooptimallysolvingtheinverseproblem. TheadaptivemeshrenementAMRalgorithmdoesindeedreducethecomputationalcosts,sinceitnecessitatesthecomputationoffewerdipolescomparedtothe volumetricinversion,butitsresultsinaidingtheinversionprocessareinconclusive. Itisdiculttodeterminewhethertheinvertedsolutiondivergesfromtheforward solutionafterafewiterationsduetotheselectioncriteriainthegridrenementprocessorduetothegeneralill-posednessoftheinverseproblem.Aproposedsolution involvesturningthisstructuredmeshrenementprocessintoastochasticone.The imaginggridwouldbeinitializedwithanitequantityofgridcellstobedetermined fromastatisticaldistributionintheglobaldomain.Subsequentiterationsrenethe selectedsub-domainusingthisstatisticaldistributionuntiltheconvergencecriteria aremet.ThisstochasticAMRprocedureisrepeatedmanytimesandtheacquired resultsareaveragedtocomputetheinvertedcurrentdensity.Whilethisremains untested,themotivationistoapplyMonte-Carlo-liketechniquestoaidinsolving thisdeterministicinverseproblem. Last,theinaccuraciesinthereconstructioncouldbeattributedtothechoiceof regularizationintheinverseproblem.Simply,the L 1 penaltyortheTikhonovpenalty maynotbesuitabletotheseclassesofproblems.Onealterationwouldbeanelastic 57

PAGE 69

netpenalty[33]whichcombinestheeectsofthesetwopenaltyterms.However,the problemofrequiringthecorrectchoiceofthepenaltytermscanlikelybeavoided entirelyusingmachinelearningmethodstosolvetheinverseproblem[1]. Whileallthesepotentialpathstoimprovingtheaccuracyofthereconstruction ofthenear-eldelectromagneticimagingproblemarepromisingcandidates,these arenotguaranteedtoaccuratelyandmeaningfullysolvethesespecicinverseproblems.Eithermoresophisticatedtechniquesarenecessaryorthattheseclassesof electromagneticforwardproblemsmustbetreatedindividually. 58

PAGE 70

CHAPTERVII CONCLUSION Presentedhereisaninversionalgorithmforsolvingthegenerallow-frequency near-eldelectromagneticinverseproblem.Theelectromagneticscatteringproblem istreatedsimplyasalinearsystemwithaforwardoperatordeterminedbytheGreen's functionresponseoftheproblem G B andthemeasuredelddatainregionsexternal totheimagingregion B sc .TheforwardproblemandthenumericalGreen'sfunction responsesarecomputedwiththeaidofthecommercialelectromagneticmethodof momentssolverFEKO.Theunknownscatterersinthisregionareidenticallyrepresentedbyequivalentvolumetriccurrentdensities J eq withunknownweights.Dueto theill-posednessoftheinverseproblemintheLFregime,itisexceedinglydicult tosolvefortheseunknownmodelparameters.Thiscanbealleviatedwiththeinclusionofregularizationandphysicalconstraintsintotheinverseproblem.Various regularizedinversiontechniquesweretestedwiththisalgorithmonthethreecases considered:aPECsphereinfree-space,aPECsphereinaconductivehalf-space, andPECcubesshieldedbyahighly-conductivethinshell.Directapplicationsofthe regularizedinverseproblemonthesethreeconsideredcasesdonotproduceconclusive ormeaningfulresults.However,theintroductionofphysicalconstraintsintotheinverseproblemconsiderablyimprovedthereconstructionsincestrongresponseswere invertedtothesurfacesoftheinternalPECcubes,thoughitdoesnotconclusively reconstructthegeometryofthescatterers.Furtherresearchmustbedoneinorder tooptimizethisinversionalgorithm. 59

PAGE 71

REFERENCES [1]JonasAdlerandOzaOktem.Solvingill-posedinverseproblemsusingiterative deepneuralnetworks. InverseProblems ,33:1{24,2017. [2]AhmedI.Al-Shamma'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 Anatomo-FunctionalPriorsintheEEG/MEGInverseProblem. IEEETRANSACTIONSONBIOMEDICALENGINEERING ,44:374{385,1997. [6]GlennD.Boreman. ModulationTransferFunctioninOpticalandElectro-Optical 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 anddenoisinginnear-eldimaging. InverseProblems ,22:1437{1456,2006. [12]K.L.Graf,N.G.Lehtinen,M.Spasojevic,M.B.Cohen,R.A.Marshall,and U.S.Inan.Analysisofexperimentallyvalidatedtrans-ionosphericattenuation 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 .Prentice-Hall,Inc.,3rd edition,1999. [16]RobertE.Grimm,BarryBerdanier,RobertWarden,JamesHarrer,Raymond Demara,JamesPfeier,andRichardBlohm.Atime-domainelectromagnetic sounderfordetectionandcharacterizationofgroundwateronMars. Planetary andSpaceScience ,57:1268{1281,2009. [17]P.C.Hansen.RegularizationToolsVersion4.0forMatlab7.3,2007. [18]PerChristianHansenandDianneProstO'Leary.TheUseoftheL-Curve intheRegularizationofDiscreteIll-PosedProblems. SIAMJ.Sci.Comput. , 14:1487{1504,1993. [19]VijayHarid.Resolutionlimitsofnear-eldelectromagneticimaging.In 2018InternationalAppliedComputationalElectromagneticsSocietySymposiuminDenver,ACES-Denver2018 ,number2,pages1{2.AppliedComputationalElectromagneticsSocietyACES,2018. [20]VijayHarid,MarkGolkowski,StephenGedney,RonaldRorrer,PooryaHosseini,MorrisCohen,NathanOpalinski,andSarahPatch.ModelingNear-Field 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 ,33-3: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 MethodsfortheSolutionofIll-PosedProblems .SpringerNetherlands,1edition, 1995. [32]StuartM.Wentworth. FundamentalsofElectromagneticswithEngineeringApplications .Wiley,2006. [33]HuiZouandTrevorHastie.RegularizationandVariableSelectionviatheElastic Net. JournaloftheRoyalStatisticalSociety.SeriesBStatisticalMethodology , 67:301{320,2005. 62