MONOTONICITY BASED IMAGING METHODS FOR PULSED EDDY

CURRENT TESTING

By

Zhiyi Su

A DISSERTATION

Submitted to

Michigan State University

in partial fulfillment of the requirements

for the degree of

Electrical Engineering—Doctor of Philosophy

2018

MONOTONICITY BASED IMAGING METHODS FOR PULSED EDDY CURRENT

ABSTRACT

TESTING

By

Zhiyi Su

Eddy current testing (ECT) is used extensively as a nondestructive evaluation (NDE)

method to detect cracks and corrosion in critical structures, particularly in the aerospace

and nuclear industries. The major advantages of ECT is its non-contact nature allowing

high-speed inspection, simplicity of operation and robustness in hazardous environments.

While the detection of defects is relatively simple and straightforward in eddy current test-

ing, the imaging problem, i.e. reconstruction of the defect profile is a challenging inverse

problem.

Iterative methods generally applied to solve the inverse problem for the defect

profile are computationally intensive and are prone to be trapped in local minima. Non-

iterative methods are required for real-time imaging capability. Monotonicity based imaging

method is one class of non-iterative methods for reconstructing defect profiles from pulsed

eddy current measurements.

This dissertation presents two major contributions to the fields of NDE and eddy current

imaging:

1. Formulated two important properties, namely monotonicity of time constants and

monotonicity of transfer function, in time-domain eddy current measurements.

2. Developed non-iterative, real-time imaging algorithms based on the monotonicity.

The transient response of the eddy current system can be characterized by a discrete set

(countable) of real, non-negative time constants. Time constants are global properties of the
α ≥ τ k
specimen under test and they are monotonic with the anomaly profile: Dα ⊆ Dβ ⇒ τ k
β ,

where Dα and Dβ are two anomalies and τ k

α and τ k

β are the associated time constants

arranged in decreasing order.

This property is exploited in this thesis to implement a non-iterative imaging method

which determines if a voxel/test element in the sample is part of the anomaly or not by

comparing the time constants. However, the implementation of this method has two major

challenges, namely non-uniqueness of solutions due to the isometry and extraction of time

constants. In this thesis, additional conductive blocks are introduced as part of the inspection

system to break the isometry. Different algorithms are implemented to extract the time

constants from transient waveforms. A Laplace-Pad`e approximant approach is shown to be

able to extract a few significant time constants.

An alternate monotonicity property found in transient eddy current measurements is the
monotonicity of the transfer function: Dα ⊆ Dβ ⇒ Zα − Zβ is negative semidefinite, where

Z is the transfer function evaluated on the negative real axis. An imaging method based

on this property is also studied and compared to the earlier method using time constants.

Validation results using simulated and experimental data are presented.

Copyright by
ZHIYI SU
2018

ACKNOWLEDGMENTS

I would like to express my most sincere gratitude to my two major advisors: Dr. Antonello

Tamburrino and Dr. Lalita Udpa. This thesis would not be possible without their support,

without their patience to my day-by-day discussion requests, and I would not have been

exposed to this very interesting topic at first place. Looking back to the five years working

with them, I truly felt it changes my belief to mathematics, physics and engineering.

I would like to thank all my committee members: Dr. Satish Udpa, Dr. Yiming Deng

and Dr. Guowei Wei. They never hesitated to provide the best help when I came across any

problems in my research.

Special thanks to Dr. Salvatore Ventre for providing the numerical software to solve

the 3D eddy current problem, to Dr. Carmine Abbate for providing the voltage-controlled

current source circuit, to Dr. Zhengfang Zhou for providing a lot of help in the theory of

partial differential equations.

I am so grateful to the NDE lab. This is the most diverse lab I’ve ever seen. There are

experts in electrical engineering, mechanical engineering, material science, etc. One can learn

and get inspired so much by just talking to them every day. They are the best colleagues

and friends. They raps and rocks!

I would like to thank Yijie. Thanks for the daily companion. Thanks for all your

suggestions, they are truly very helpful.

And finally, last but by no means least, to my parents, for the care, for the listening, for

the support, for always being there for me.

Thanks for all your encouragement!

v

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

ix

Chapter 1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1 Nondestructive Testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1.1 Overview of Other Nondestructive Testing Methods . . . . . . . . . .
1.2 Electromagnetic Testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 Eddy Current Testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.1 Advantages of Eddy Current Testing . . . . . . . . . . . . . . . . . .
1.3.2 Limitations of Eddy Current testing . . . . . . . . . . . . . . . . . . .
1.4 Pulsed Eddy Current Testing . . . . . . . . . . . . . . . . . . . . . . . . . .

Chapter 2

Inverse Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1 General Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2
Inverse Problem in NDT . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 Solution to Inverse Problems . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.1
Iterative Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.2 Non-iterative Methods . . . . . . . . . . . . . . . . . . . . . . . . . .

Chapter 3 Monotonicity Principle Method in PECT problems . . . . . . .
3.1 Mathematical Model of Pulsed Eddy Current Problems . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.1 Mathematical Model
3.1.2 Numerical Model
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.3 Natural Modes for Eddy Current Problems . . . . . . . . . . . . . . .
3.2 Monotonicity of Time Constants . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.1 Monotonicity Principle . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.2
Isometry and Time Constants . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . .
Interaction between Natural Modes and Anomalies
3.2.3
3.3 Application of MPM to Imaging . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.1
Imaging Method based on Monotonicity Principle . . . . . . . . . . .
3.3.2 Consideration of Noise . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Chapter 4 Estimation of Time Constants . . . . . . . . . . . . . . . . . . . .
4.1 Previous Works on Exponential Analysis . . . . . . . . . . . . . . . . . . . .
4.2 A Modified Method Based on the Laplace-Pad´e Approximant Approach . . .
4.3 Numerical Examples on Estimating Time Constants . . . . . . . . . . . . . .
4.4 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1
1
1
4
5
6
8
8

11
11
12
13
14
16

19
19
19
23
25
26
28
28
31
35
35
37
38

39
39
43
50
53

vi

Chapter 5 Numerical Examples . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1 Validation of Monotonicity . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2 Defect Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3 Application of MPM in Defect Sizing . . . . . . . . . . . . . . . . . . . . . .
5.4 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Chapter 6 Monotonicity of Transfer Function and Imaging Method . . .
6.1 Definition of Transfer Function . . . . . . . . . . . . . . . . . . . . . . . . .
6.2 Monotonicity of Transfer Function and Imaging Scheme . . . . . . . . . . . .
6.3 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.3.1 Measurement of Transfer Function . . . . . . . . . . . . . . . . . . .
6.3.2 Noise Free Reconstruction . . . . . . . . . . . . . . . . . . . . . . . .
6.3.3 Noisy Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.4 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Chapter 7 Experiment Setup and Results . . . . . . . . . . . . . . . . . . . .
7.1 Experiment Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.2
Imaging Scheme and Results . . . . . . . . . . . . . . . . . . . . . . . . . . .

55
55
56
62
65

66
66
69
71
71
73
77
84

86
87
91

97
Chapter 8 Conclusions and Future Work . . . . . . . . . . . . . . . . . . . .
8.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
97
8.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

APPENDIX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104

vii

LIST OF TABLES

Table 3.1: The perturbation of time constants when introducing a surface breaking
anomaly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33

Table 3.2: The perturbation of time constants when introducing a buried anomaly .

34

Table 4.1: Estimation of time constants from synthesized signals in Example 1.

. . .

52

Table 4.2: Estimation of time constants from synthesized signals in Example 2.

. . .

54

Table 5.1: Classification of anomalies via time constants . . . . . . . . . . . . . . . .

65

viii

LIST OF FIGURES

Figure 1.1: Principle of eddy current testing method. . . . . . . . . . . . . . . . . . .

Figure 1.2: A typical PECT signal and common time-domain features. . . . . . . . .

6

9

Figure 2.1: Flowchart of the iterative methods to solve inverse problems . . . . . . .

14

Figure 2.2: Flowchart of the non-iterative methods to solve inverse problems . . . . .

17

Figure 3.1: Schematics of a typical ECT system.

. . . . . . . . . . . . . . . . . . . .

20

Figure 3.2:

Isometry in Euclidean space: translation, reflection and rotation.

. . . .

29

Figure 3.3:

If two configurations can be transformed into each other by an isome-
try, they have identical time constants. (a) Schematic of two isometric
configurations, and (b) the related time constants. . . . . . . . . . . . . .

Figure 3.4:

Introducing an additional conductive patch can help distinguish two iso-
metric configurations. (a) Schematic of the test specimen together with
. . . . . . . . . . . . . . .
the patch, and (b) the related time constants.

Figure 3.5: The natural modes of a conductive cube. (a) The cube, and (b)-(l) the
natural modes associated with the k−th time constant. . . . . . . . . . .

Figure 3.6: The imaging method exploiting the union of test elements. (a) If the test

element Tj is part of the anomaly V , then τ

i , for all i. (b) If there
i , then Tj is not part of V . . . . . . . . . . .

Tj
i < τ V

exist an i such that τ

i ≥ τ V
Tj

Figure 3.7: The imaging method exploiting the intersection of test elements. (a) If

the anomaly V contains the test element Tj, then τ V
If there exists an i such that τ V

Tj
i < τ
i

, then V does not contain Tj. . . .

i ≥ τ

Tj
i

, for all i. (b)

Figure 4.1: The example of Lanczos. Twenty-four data points (circles) are fitted
by a double-exponential function f1(t) = 2.202e−t/4.45 + 0.305e−t/1.58
(blue solid line) and by a triple-exponential function f2(t) = 0.0951e−t +
0.8607e−t/3 + 1.5576e−t/5 (red dash line).
. . . . . . . . . . . . . . . . .

30

31

32

35

36

41

Figure 5.1: 3D mesh of the test structure: a cylinder pipe section . . . . . . . . . . .

56

Figure 5.2: Validation of the monotonicity of time constants . . . . . . . . . . . . . .

57

ix

Figure 5.3: A homogeneous slab as the test specimen. (a) Finite element mesh. (b)
. . . . . . . . . . . . . . . . . . . . . . . . . . . .

The region of interest.

Figure 5.4: The inspection system consists of (a) a 2D coil array, and (b) a large
conductive patch above the specimen. . . . . . . . . . . . . . . . . . . . .

Figure 5.5:

Imaging example using union of test anomalies, assuming noise level ∆ =
0. The red dashed frames represent the target anomaly. . . . . . . . . . .

Figure 5.6:

Imaging example using intersection of test anomalies, assuming noise level
∆ = 0. The black bold frames represent individual test anomalies that
are found through the monotonicity test (3.47). The red dashed frames
represent the target anomaly.
. . . . . . . . . . . . . . . . . . . . . . . .

Figure 5.7:

Imaging example with various noise levels. The red dashed frames repre-
. . . . . . . . . . . . . . . . . . . . . . . . . . .
sent the target anomaly.

Figure 5.8:

Imaging example with various numbers of available time constants, as-
suming noise level ∆ = 0.1%. The red dashed frames represent the target
anomaly. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 5.9:

Illustration of the numerical examples. The coil array is above the test
. . . . . . .
specimen. The dark colored voxels represent the air domain.

Figure 5.10: Illustration of the anomalies of various sizes. (a) small anomaly (S), (b)
medium anomaly (M) and, (c) large anomaly (L). . . . . . . . . . . . . .

58

58

59

60

61

62

63

63

Figure 5.11: The time constants associated with anomalies of different sizes group into

three clusters. The measured signals have (a) 0.1% noise and (b) 1% noise. 64

Figure 6.1: An example of the excitation waveform . . . . . . . . . . . . . . . . . . .

67

Figure 6.2: An example problem: the mesh of specimen and the transmitter/receiver
coil setup. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 6.3: Transient waveform and transfer function. (a) input current in coil 1 (b)
pick-up voltages (c) impedance, the asymptotic values are the transfer
function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

72

74

Figure 6.4: Sensor array configurations.

. . . . . . . . . . . . . . . . . . . . . . . . .

75

Figure 6.5: Example defect profiles.

. . . . . . . . . . . . . . . . . . . . . . . . . . .

75

x

Figure 6.6: Reconstruction via monotonicity of transfer function using noise free data.
Top figure: real defect profiles (in red color). Bottom figure: reconstructed
defect profiles (in light blue color) . . . . . . . . . . . . . . . . . . . . . .

Figure 6.7: Statistical reconstruction results for defect #1, τ = 5τ1 . . . . . . . . . .

Figure 6.8: Statistical reconstruction results for defect #2, τ = 5τ1 . . . . . . . . . .

Figure 6.9: Statistical reconstruction results for defect #3, τ = 5τ1 . . . . . . . . . .

Figure 6.10: Statistical reconstruction results for defect #4, τ = 5τ1 . . . . . . . . . .

Figure 6.11: Statistical reconstruction results for defect #1, τ = 2τ1 . . . . . . . . . .

Figure 6.12: Statistical reconstruction results for defect #2, τ = 2τ1 . . . . . . . . . .

Figure 6.13: Statistical reconstruction results for defect #3, τ = 2τ1 . . . . . . . . . .

Figure 6.14: Statistical reconstruction results for defect #4, τ = 2τ1 . . . . . . . . . .

76

78

79

80

81

82

83

84

85

Figure 7.1: Schematic of experiment setup . . . . . . . . . . . . . . . . . . . . . . .

86

Figure 7.2: Design of the voltage controlled current source.

. . . . . . . . . . . . . .

88

Figure 7.3: The voltage controlled current source is connected to a pair of transmitter-
receiver coils for testing performance. . . . . . . . . . . . . . . . . . . . .

89

Figure 7.4: Frequency response of the current source circuit. . . . . . . . . . . . . . .

90

Figure 7.5: The excitation current agrees very well with the input signal.

. . . . . .

91

Figure 7.6: The two coils consisting the array. The smaller coil is inserted into the
. . . . . . . . . . . . . . . . . . . . . . . . . . . .

larger one when used.

91

Figure 7.7: Three artificial defects are fabricated on an aluminum test sample.

. . .

92

Figure 7.8: A test element and the coil positions where the measurements are relevant. 94

Figure 7.9: The measurements of experimental current and voltage waveforms and
their fitting to exponential curves. . . . . . . . . . . . . . . . . . . . . . .

Figure 7.10: Experimental results using the imaging approach based on monotonicity
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

of transfer function.

95

96

xi

Chapter 1

Introduction

1.1 Nondestructive Testing

Nondestructive testing (NDT) is the process of inspecting, testing, or evaluating materials,

components or assemblies for discontinuities or anomalies without destroying the service-

ability of the part or system [1]. A variety of nondestructive techniques are used in industry

to detect variations in structure, changes in surface finish, the presence of cracks or other

physical discontinuities, and assess the integrity of critical structures and components [2].

Modern nondestructive tests are widely used in almost every industry, civil engineering

[3], automotive [4], aircraft [5], nuclear [6], renewable energy [7], etc. NDT is also used

as a control mechanism to ensure that manufacturing processes meet design performance

requirements and in turn reliability.

1.1.1 Overview of Other Nondestructive Testing Methods

In order to use a nondestructive testing effectively, it is necessary first to understand the

principles of the method. The six most frequently used test methods according to the Amer-

ican Society for destructive Testing (ASNT) are Visual Testing, Magnetic Particle Testing,

Liquid Penetrant Testing, Ultrasonic Testing, Electromagnetic Testing and Radiographic

Testing. The following section briefly describes these methods and their applications.

1

• Visual Testing. As the name implies, visual testing is the visual observation of a test

object to evaluate the presence of surface discontinuities. Visual testing could be per-

formed directly using eyes, or indirectly using optical instruments such as magnifying

glasses, mirrors, borescopes or computer-assisted viewing systems. The test usually re-

quires proper cleaning of the surface and adequate illumination. It is the most widely

used method for detecting and examining surface discontinuities. The following sur-

face discontinuities may be detected by a visual test: cracks, misalignment, warping,

corrosion, wear and physical damage.

• Magnetic Particle Testing. Magnetic Particle Testing makes use of leakage mag-

netic field to locate surface and near-surface discontinuities in ferromagnetic materials.

When the test object is magnetized, discontinuities that lie in a direction generally

transverse to the direction of the magnetic field will cause a leakage field at and above

the discontinuities. Very fine colored ferromagnetic particles (”magnetic particles”)

are then applied to the surface of the test object. The magnetic particles are attracted

to the leakage field around discontinuities and hence, form visible indication of the

presence of this leakage field and therefore the discontinuity itself. The magnetic field

can be applied using a permanent magnet or an electromagnet such as yokes, prods,

electric coils or central conductor. Some of the typically detected discontinuities are

surface discontinuities, seams, cracks and laps.

• Liquid Penetrant Testing. Liquid penetrant testing reveals fissures and voids on the

surfaces of solid and nonporous materials. The liquid penetrant is a liquid of very low

viscosity, when applied to the surface of the test object it seeps into the open-surface

fissures and voids. Once the excess penetrant is removed, the trapped penetrant forms

2

a visible indication and the test object is inspected visually by eyes, or with aid of

an ultraviolet light in case of fluorescent penetrants. The test requires the surface

of interest to be clean and free of any material or liquids that prevent the penetrant

from entering the voids and fissures. The following surface discontinuities are typically

detected by a liquid penetrant inspection: surface discontinuities, seams, cracks, laps,

porosity and leak paths.

• Ultrasonic Testing. Ultrasound refers to acoustic waves of frequencies higher than

the upper audible limit of human hearing and it does not travel through air. Ultrasonic

testing is a nondestructive method in which beams of acoustic waves are introduced

into test object for the detection of discontinuities in the material. These acoustic

waves travel through the material with some attenuation and are reflected at disconti-

nuities, or more precisely, at the interfaces between two materials of different acoustic

impedance. The reflected waves are then collected and analyzed to detect and locate

the discontinuities. Sound waves of lower frequencies have better penetration but are

less sensitive to small defects. The excitation sound is usually injected by means of

an ultrasonic transducer which converts electrical impulses into acoustic waves. Since

ultrasound does not propagate in air, usually a liquid or gel called ”couplant” is used

between the surface of the transducer and the test object. Ultrasonic testing is very

widely used to detect internal discontinuities in most engineering metals and alloys, as

well bonds produced by welding, brazing, soldering and adhesives.

• Radiographic Testing. Radiographic testing involves exposing the test object to

penetrating radiation – either electromagnetic, e.g. x-ray and gamma ray or particulate

radiation, e.g. positron, neutron – and a recording medium on the other side of the

3

test object. Different portions of the test object absorb different amount of penetrating

radiation depending on the differences in density, thickness, absorption rate, etc. Thus,

by monitoring the unabsorbed radiation that passes through the object and arrives at

the recording media, one can discover variations in the absorption of radiation and

therefore the variation in material properties. Typically detected discontinuities and

conditions in radiographic testing include inclusions, lack of fusion, cracks, corrosion,

porosity, leak paths, missing or incomplete components and debris.

There are many other nondestructive testing methods such as Acoustic Emission Testing,

Ground Penetrating Radar, Thermal/Infrared Testing and Vibration Analysis.

1.2 Electromagnetic Testing

Electromagnetic methods are an important and very widely used method in nondestructive

testing. The electromagnetic test method involves exposing the test object to an electro-

magnetic field and measuring the response of field/flaw interaction, in order to evaluate the

property of the test object. Electromagnetic testing methods cover a wide range of the fre-

quency spectrum and include magnetic flux leakage testing, eddy current testing (ECT) and

microwave testing [2]. Magnetic flux leakage testing typically uses excitation frequencies

near 0 Hz. Eddy current testing uses excitation frequencies from about 100 Hz to about

10 MHz, where the electromagnetic field is said to be quasi-static. Microwave testing uses

excitation sources usually beyond 10 MHz, where the energy propagates in the form of waves

into the test object.

4

1.3 Eddy Current Testing

Eddy Current Testing is an electromagnetic NDT technique operating at quasi-static regime,

namely from 100 Hz to about 10 MHz, where the displacement current is negligible. ECT

relies on the principles of electromagnetic induction to interrogate test object. When a

coil excited by alternating source current is brought close to the specimen, the primary

field Bs generated by the coil induces eddy currents Js within the electrically conductive

specimen. According to the Lenz’s law, the secondary field Beddy produced by these eddy

currents opposes the change in the primary field (Figure 1.1) thereby altering the terminal

impedance of the coil.

Conventional ECT measures the terminal impedance of the coil. If the test specimen is

non-ferromagnetic, the magnetic flux linked with the coil decreases since the secondary field

Beddy opposes the primary field Bs. Because the self-inductance of the coil is linear with

its flux linkage, the inductance of the coil decreases. At the same time, the resistance of the

coil increases due to the fact that eddy current losses incurred within the specimen have to

be met by the source of primary excitation.

Eddy current is sensitive to changes of the electromagnetic properties, namely electrical

resistivity η and permeability µ, within the specimen. For instance, the presence of a discon-

tinuity in the specimen causes a reduction as well as a redistribution of the eddy currents.

Consequently the changes in the inductance and resistivity of the coil are reduced.

In case of ferromagnetic materials, when the coil is brought close to a ferromagnetic

specimen, its inductance increases due to the higher permeability of the material and there

is an increase in its resistance due to the eddy current and hysteresis losses.

The secondary field generated by eddy currents can also be measured using separate

5

pick-up coils or field sensors. The terminal voltage of a pick-up coil is linearly dependent on

the rate of change of the magnetic flux linkage associated with it:

v(t) = −N

dφ(t)

dt

,

(1.1)

where v(t) is the induced voltage, φ(t) is the magnetic flux and N is the number of turns of

the coil. Hence, pick-up coils are used to measure only varying magnetic field, usually of a

moderate frequency which is not too low. Magnetic field sensors based on different sensing

technologies such as Hall Effect, Anisotropic Magnetoresistance (AMR), Giant Magnetoresis-

tance (GMR), and others such as Tunneling Magnetoresistance (TMR) and Superconducting

Quantum Interface Device (SQUID) are also used to measure the induced magnetic field di-

rectly..

Figure 1.1: Principle of eddy current testing method.

1.3.1 Advantages of Eddy Current Testing

Application of eddy current tests are numerous, particularly in automotive, nuclear [8, 9],

aircraft [10, 5] and renewable energy [7, 11] industries. Modern eddy current techniques offer

low cost means for inspecting a variety of materials including metals, alloy and composites.

For instance, carbon-fiber-reinforced plastic (CFRP) materials, which are commonly seen

in components of vehicles or wind turbines, can be inspected using eddy current method

6

[11, 12]. In particular, the main advantages of ECT are:

• Versatile applications.

In the proper circumstances, eddy currents can be used for

crack detection, material thickness measurement, coating thickness measurement, and

conductivity measurement.

• Low cost. Minimum part preparation, wire-would coil(s), is required.

• Simplicity of operation. Very little pre-cleaning and preparation time is required. On

relatively uniform parts, eddy current method can inspect very quickly and be easily

automated.

• Non-contact. This feature allows eddy current detect through non-conductive surface

coatings, or work under harsh operation conditions that other NDT methods are not

capable of[13].

• Rapid inspection. Eddy current testing gives immediate results regarding crack detec-

tion applications.

• Works with single side access. The nature of magnetic field penetration allows eddy

current to detect defects in multi-layer structures without interference from the planar

interfaces. This is very valuable particularly for in-service inspections where only one-

side access is possible.

Consequently, the ECT method is widely used in a production environment for detection

of anomalies in materials.

7

1.3.2 Limitations of Eddy Current testing

The Limitations of eddy current testing are direct consequences of electromagnetic fields.

Eddy current tests provide maximum sensitivity to surface or near-surface layers of the test

object due to the fact that eddy current density decreases exponentially from its value at the

surface according to the depth from the surface, which is commonly recognized as the skin

effect. In general, eddy current testing are only applicable to test materials with significant

electrical conductivity, where eddy currents are induced. Also, eddy current signals are not

sensitive to discontinuities that lie parallel to the induced eddy currents. The EC methods

produces a strong signal in the case of discontinuities that lie in the direction generally

transverse to the flow of eddy currents, as they significantly modify the current flow paths.

At last, eddy current signals are sensitive to the coil lift-off and tilt. They are usually

considered as noise sources and are undesired in defect detection.

1.4 Pulsed Eddy Current Testing

Pulsed eddy current testing (PECT) employs excitation waveforms of repetitive pulses, such

as square, triangular or saw tooth waveforms. Compared to conventional eddy current

method which exploits a single frequency waveform, PECT is advantageous since it contains

a large spectrum of frequencies and hence, provides more information that can be used in

the detection and characterization of defects.

Figure 1.2 shows a typical PECT signal and its time-domain features. Analysis of PEC

data is usually carried out in the time domain. Some commonly used features are the peak

value, the arrival time of the positive peak, the arrival time of zero-crossing, the arrival time

of the rising point, and the arrival time of the descending point [14]. Other signal process

8

techniques exploiting the statistical features such as principle component analysis (PCA)

[15] have also been used, and showed promising results in defect detection and classification.

However, these features are sensitive to probe lift-off and tilt, as well the driving system and

excitation waveforms. In some cases, the variation due to lift-off can be compensated by

introducing extra reference signals [16].

Figure 1.2: A typical PECT signal and common time-domain features.

In PECT, the induced eddy current density J, after the excitation is switched off, can be

represented as sum

+∞(cid:88)

i=1

J(r, t) =

ciJi(r)e−t/τi

(1.2)

where Ji(r) are the natural modes, τi the time constants and ci the linear coefficients. In this

report, we present time constants as a new feature of pulsed eddy current testing for defect

detection and identification. Time constants characterize the source-free response in PECT

and they depend only on intrinsic properties of the test specimen.

In other words, time

constants are invariant w.r.t. probe lift-off, tilt, driving system and excitation waveforms.

Pulsed eddy current method has been applied extensively by several researchers. Yang et

9

al. [5] exploited pulsed eddy current excitation and giant magnetoresistive (GMR) sensors

to detect cracks as small as 1 mm in the 2nd and 3rd layer of a three-layer structure, which

demonstrates the potential of this method in detecting flaws at large depths. Multi-modality

is another research area in PEC techniques. Cheng et al. [12] combined PEC and infrared

imaging to successfully detect defects over a relatively large area within a short time.

This report consists of four parts. In the first part (Chapter 1 and Chapter 2), a brief

introduction to the background of this study and important concepts of inverse problems

in NDT are given. The second part (Chapter 3 to Chapter 5) starts off with a mathe-

matical model of the eddy current problem, introduces the concepts of natural modes and

time constants, and proves the monotonicity property of time constants. A few important

properties of time constants are also explored, including the global interaction and isometry.

Finally, a non-iterative imaging algorithm is presented and validated with numerical models,

along with addressing the practical problem of extracting time constants. The third part

(Chapter 6 to Chapter 7) focuses on another monotonicity property of eddy current problem:

the monotonicity of transfer function. While time constants are related to the source-free

response of the problem, transfer function is derived from forced response. Monotonicity

of time constants and monotonicity of transfer function are two complementary properties

that can be applied to the problem of eddy current imaging. Concluding remarks and some

considerations on the extension of this study are given at the end of this report (Chapter 8).

10

Chapter 2

Inverse Problem

The inverse problem in nondestructive evaluation is defined as estimation of material prop-

erty distribution in a structure, given the NDE sensor measurements. In this chapter we

introduce the basic concept of inverse problem and its definition in nondestructive testing.

We discuss the general ideas of solving an inverse problem, including iterative methods and

non-iterative methods.

2.1 General Remarks

Consider the equation

A(x) = y

(2.1)

where A is a continuous operator mapping a subset X of a Banach space to a subset Y of
another Banach space, and the target is to find x ∈ X given y ∈ Y . This problem is said to

be well-posed if the following conditions hold [17]:

1. for any y ∈ Y there exists a solution x ∈ X (existence);
2. for any y ∈ Y there is only one solution x ∈ X satisfying (2.1) (uniqueness);
3. (cid:107)x − x∗(cid:107)X → 0 when (cid:107)y − y∗(cid:107)Y → 0 (stability).

(2.2)

(2.3)

(2.4)

11

In (2.4), x∗ ∈ X, y∗ ∈ Y and A(x∗) = (y∗). (cid:107)·(cid:107)X and (cid:107)·(cid:107)Y are proper norms defined in X

and Y , respectively. They compute the distance between two elements in the corresponding

space.

If any of the conditions from (2.2) to (2.4) is not satisfied, the problem (2.1) is said to be

ill-posed. Inverse problems are typically ill-posed. Of the three conditions for a well-posed

problem, the condition of stability (2.4) is most often violated. In other words, small errors

in the measurement of y are greatly amplified in the solution x.

A most common approach for solving (2.1) is by means of regularization, which essentially

replaces the original problem with an alternative equation that includes a regularization term

[17]. The alternative equation can be solved in a stable way and its solution is ”close” to that

of the original problem. The nature of regularization is to introduce additional information,

i.e. the smoothness or bounds on the solution x, to prevent over-fitting.

2.2

Inverse Problem in NDT

Inverse problem in NDT involves the characterization of the specimen parameters given an

NDT probe signal [18]. This usually involves the evaluation of material properties, which is

equivalent to detection/imaging of anomalies. Adapting the idea of (2.1), X is the parameter

space of the specimen and Y is NDT probe signal in this context. Ideally, one wishes to
develop an analytical model of A that gives a closed form representation of the probe signal

given the specimen geometry. However, this is usually not possible in NDT problems due to

the material nonlinearity and arbitrary shapes of the specimen and anomalies. This results

in the use of numerical methods such as finite element modeling to describe the model of
forward system A. However, numerical models do not yield a closed form representation of

12

the problem making inverse problem solution more difficult. A typical model-based approach

for solving the inverse problem follows an iterative procedure as described in section 2.3.1.

This report focuses on the problem of imaging an anomaly given the pulsed eddy current

signals. Let A be the operator for mapping a prescribed resistivity distribution η(r) onto
the measured quantity m, i.e. A : η → m. The inverse problem in PECT problems is then

defined as estimation of η(r) given the measured quantity m. Typically, the inverse problem

is formulated in terms of minimization of an appropriate cost function plus a regularization

term [19]:

(cid:107)A (η) − ˜m(cid:107)2 + αΩ(η)

min
η∈C

(2.5)

where ˜m is the experimental, noisy measurement, Ω a penalty term which takes into account

a priori information about the unknown resistivity distribution η, α the corresponding regu-
larization parameter, C the constraint set and (cid:107)·(cid:107), an appropriate norm. η is then estimated

by minimizing the cost function over C. For instance, C can be the set of non-negative

electrical resistivity distributions.

2.3 Solution to Inverse Problems

This section discusses model-based methods to solve the inverse problem. There are other

possible non model-based methods, for instance using neural networks to map the measure-

ment data to the parameter space [20, 21]. However, it is not effective to train networks to

cover a broad class of defects. The trained networks only provide good performance over a

narrow range of signals compatible with the training dataset [22]. Model based solutions to

inverse problems can be generally classified as iterative methods and non-iterative methods.

13

2.3.1 Iterative Methods

The iterative methods start from an initial guess of the resistivity distribution η0. The

numerical model is solved to predict the measurement signal. The error between simulation

and measurement forms the cost function. The solution η is then iteratively updated to

reduce the value of objective function (2.5), until a stopping criterion is met. This approach

is presented schematically in Figure 2.1.

Figure 2.1: Flowchart of the iterative methods to solve inverse problems

Several iterative methods have been applied to solve the minimization problem in (2.5)

and they have been proved to be quite successful in solving NDT inverse problems [23,

14

24]. However, iterative algorithms suffer from: (i) false solutions since the algorithm may

be trapped in local minima, and (ii) high computational cost. Generally speaking, high

computational cost is the bottleneck for iterative methods. All these algorithms require to

solve the forward problem (sometimes the gradient of the forward problem as well), which is

usually a 3-D FEM model in ECT problems, at least once in each iteration. Further, with

Gauss-Newton algorithms including the steepest descent algorithm, the conjugate gradient

algorithm, and the Levenberg-Marquardt algorithm, the computational cost increases also

due to the need for computing the gradient of the objective function. This typically requires

the solutions of N forward problems at every iteration, N being the number of degrees of

freedom of the electrical resistivity in discrete domain. Norton and Bowler [25] reduced the

computation of N forward problems to two by deriving an explicit expression for the gradient

in terms of solutions of an ”ordinary” forward problem corresponding to the current estimate

of the resistivity, and an ”adjoint” problem. We also mention Lionheart and Soleimani [26]

who derived a formula for sensitivity of the induced voltage in a measurement coil due to a

small perturbation in the resistivity.

Also, it is possible to avoid the risk of local minima by using meta-heuristic algorithms to

update the estimation of solution. Commonly seen meta-heuristic algorithms used in NDT

studies include simulated annealing, tabu search or genetic algorithms [27, 28]. Formally,

a meta-heuristic is a high-level design to find, generate or select a subordinate heuristic

for intelligently exploring and exploiting the search space, in order to achieve global or

near-global optimal solutions [29]. For instance, genetic algorithms use mechanisms inspired

by biological evolution, such as reproduction, mutation, recombination and selection.

It

is called a population-based meta-heuristic algorithm. These algorithms again requires a

large number of forward problem solutions, and is hence computational intensive. With

15

development in large-scale parallel computation techniques, Yusa et al. [30] conducted a tabu

search algorithm on a supercomputer using up to 128 CPUs to reconstruct the natural stress

corrosion crack profile with piecewise constant resistivity values in a reasonable computation

time (600 seconds).

2.3.2 Non-iterative Methods

There are practical situations where it is sufficient to retrieve only the support of an anomaly

in an otherwise homogeneous background. An example is nondestructive testing of conduct-

ing materials where the target is to detect, locate and image fatigue and corrosion cracks

in homogeneous structures.

In these cases, non-iterative imaging methods are appealing

because of their real-time, or near real-time, capabilities. The common concept in this class

of methods is the computation of an indicator, i.e. a function of the space, assuming differ-

ent values if evaluated inside or outside the anomalous region. This approach is presented

schematically in Figure 2.2.

A few of the non-iterative imaging methods are mentioned here. Colton and Kirsch in-

troduced the first non-iterative approach named Linear Sampling Method (LSM) [31] for

wave-propagation phenomena. The Factorization Method (FM) was proposed by Kirsch

[32] for wave-propagation phenomena and extended by Hanke and Br¨uhl to Electrical Re-

sistance Tomography (ERT) [33]. MUSIC (MUltiple SIgnal Classification) algorithm is well

known in signal processing as a method to estimate the constant parameters upon which

the measurement depends such as frequency estimation. Devaney [34] pointed out that it

could also be used in imaging problems. Ammari and Lesselier [35] developed a MUSIC

algorithm for locating small inclusions buried in a half-space, given the scattered amplitude

data. This method is able to detect multiple inclusions at the same time but limited to the

16

Figure 2.2: Flowchart of the non-iterative methods to solve inverse problems

case that inclusions are small and well separated. Eventually, the Monotonicity Principle

Method (MPM) was introduced by Tamburrino and Rubinacci [36] and applied to electrical

resistance tomography (ERT). In this method a simple test determines if a given small test

anomaly is contained in the unknown anomaly or not. MPM was first introduced in the area

of ERT [36] (elliptic PDE problems) and then extended to parabolic PDEs such as those

governing eddy current testing in small skin-depth regime [37] and large skin-depth regime

[38] with experimental validations in [39]. MPM was also extended to hyperbolic PDEs such

as those governing wave propagation inverse scattering problems [40]. MPM provides upper

17

and lower bounds for the unknown object in the case of a finite number of measurements

(limited aperture data) [36]. In [41] it was proven that MPM provides the exact shape of

the anomalies in the ideal case of an infinite number of measurements (full aperture data),

under the hypothesis that each connected component of the anomalies is contractible.

Among non-iterative methods, MPM is a topic of active research and has been shown to be

promising by many investigators. In [42] the authors have formulated a monotonicity-based

shape reconstruction scheme and regularized against noise and modeling error. Further, two

reconstruction algorithms are formulated and validated by numerical examples with both

simulated complete electrode model (CEM) data and experimental measurement data. In

medical imaging area, a novel bimodal imaging technique combining both breast microwave

radar (BMR) and electrical impedance tomography (EIT) methods are presented to form

a resistivity distribution map of a breast region that can be used to assess the presence of

malignant lesions [43, 44]. This technique uses a priori information obtained from BMR

images to estimate the location of the dense breast regions and then, the monotonicity of

the impedance matrix is used to reconstruct a profile of the tissue distribution in the breast

region.

Monotonicity of time constants in pulsed eddy current problems was first introduced and

numerically validated in [45]. The imaging method based on this property was developed in

subsequent publications [46, 47]. These studies demonstrated that MPM of time constants

is capable of fast reconstruction of an anomaly within practical 3D structures.

18

Chapter 3

Monotonicity Principle Method in

PECT problems

3.1 Mathematical Model of Pulsed Eddy Current

Problems

This section first describes the specific eddy current system setup for ECT problems, based

on which we deduce the mathematical model describing the relationship between the eddy

current density and given excitation current source. The domain is discretized for numerical

implementation, and the concepts of natural modes and their related time constants are

introduced.

3.1.1 Mathematical Model

The system setup of ECT is showed in Figure 3.1. It is assumed that: (i) the probe comprises

an array of Nc coils and (ii) the measurement consists of voltages across this set of coils for

prescribed currents injected in the same set of coils. Hereafter ik (t) (k = 1, ..., Nc) stands
for the current impressed at the k−th coil and vk (t) stands for measured voltages across the
k−th coil. The system is operated in time-domain.

The forward problem, that is the computation of the induced voltages, can be conve-

19

Figure 3.1: Schematics of a typical ECT system.

niently solved by describing the underlying eddy current problem via an integral equation in

terms of induced eddy current density J. It is assumed that (i) the conducting domain D is
surrounded by an insulating material which implies that J · ˆn=0 on ∂D, where ˆn is the unit

normal vector of the boundary ∂D, (ii) the measurement system operates in a contactless

manner and (iii) the current circulating in the excitation coils is impressed and controlled

by the measurement system.

Assuming that the conductive material is linear, non-dispersive, isotropic, non-magnetic

and neglecting the displacement current, the eddy current density J satisfies the following

equations (rotation, gradient and divergence are meant in weak sense) [48, 49]:

ηJ= − ∂tAi [J] − ∇ϕ − ∂tAc, in D

∇ · J = 0, in D
J · ˆn=0 on ∂D

(3.1)

(3.2)

(3.3)

20

where D is the conducting domain, η (·) is the electrical resistivity of the conductor, ϕ is

the electric scalar potential, Ac is the magnetic vector potential due to the assigned driving
currents and the integral operator Ai, acting on spatial coordinates, is defined as

(cid:90)

D

J(cid:0)r(cid:48), t(cid:1)
(cid:107)r − r(cid:48)(cid:107)dV (cid:48).

Ai : J (r, t) → µ0
4π

(3.4)

(3.5)

Following [38], the vector potential Ac can be expressed as:

(cid:88)

Ac (r, t) =

Ac,k (r) ik (t) ,

where Ac,k is the magnetic vector potential produced by a unit current flowing into the k-th

k

coil:

(cid:90)

γk

µ0
4π

dl(cid:48)

(cid:107)r − r(cid:48)(cid:107),

Ac,k (r) =

where γk is the curve representing the path of the wire of k−th coil.

From the Faraday-Neumann-Lenz law, the voltage across the k−th coil is given by

vk (t) = ηkik (t) +

dφc
k (t)
dt

+

dφeddy
dt

k

(t)

,

where

φc
k (t) =

φeddy
k

(t) =

NC(cid:88)
(cid:90)

k=1
µ0
4π

L0,knin (t) ,

(cid:90)

dl(cid:48) · J(cid:0)r(cid:48)(cid:48), t(cid:1)
(cid:107)r(cid:48) − r(cid:48)(cid:48)(cid:107) dV (cid:48)(cid:48),

γk

D

(3.6)

(3.7)

(3.8)

(3.9)

ηk is the resistance of the k−th coil, φc

k is the free-space magnetic flux linked with the k−th

21

is the magnetic flux linked with the k−th coil due to the eddy current density

coil and φeddy

k

J.

The variational formulation (w.r.t. the spatial coordinates) for the eddy current problem

(Eq. (3.1)-(3.3)) isi

(cid:90)

D

J(cid:48)(r)·ηJ(r, t)dV = − ∂t

(cid:90)

D

J(cid:48)(r)·Ai [J] dV − ∂t

(cid:90)

D

J(cid:48)(r)·Ac(r, t)dV , ∀J(cid:48) ∈ ˆH (3.10)

where ˆH = {v ∈ H (div; D)|∇ · v =0, v · ˆn=0 on ∂D}, ˆn is the unit normal vector of the
boundary ∂D and J (·, t) ∈ ˆH for any given t. This implies that there is no source or sink

of eddy currents within the conducting domain D and the eddy current does not flow out of

D, which is imposed by the assumption that D is surrounded by an insulating material.

(cid:2)i1(t), . . . , iNc(t)(cid:3)T

and

(3.11)

(3.12)

Let

i(t)

and φeddy(t)

defined

as

i(t)

=

φeddy(t) =

φeddy
1

(t), . . . , φeddy
Nc

. It is worth noting that the operators

(cid:104)

be

(cid:105)T

(t)

M : i (t) → Ac (r, t)|D
N : J (r, t) → φeddy (t)

are one the adjoint of the other, i.e. MT = N [50]. Indeed,
i(cid:82)
D J(cid:48) · ∇ϕ dV = 0 because ∇ · J(cid:48) = 0 in D and J(cid:48) · ˆn = 0 on ∂D.

22

(cid:88)
(cid:88)
(cid:88)

k

k

(J,Mi)L2(D) =

=

=

(cid:90)

D

ik(t)

J · Ac,k(r)ik(t)dV
(cid:90)

(cid:20) µ0

(cid:90)

4π

D

γk

J(r, t)
(cid:107)r − r(cid:48)(cid:107)dV

φeddy
k

(t)ik(t)

(cid:21)

· dl

(3.13)

k

= (N J, i)RNC

.

In the following, we assume that η ∈ L∞ (D) and that 0 < ¯η0 ≤ η (r) ≤ ¯η1 < +∞ where

¯η0 and ¯η1 are two constants.

3.1.2 Numerical Model

The numerical model (see [48, 49]) is obtained by introducing the electric vector potential
T (J =∇ × T) to enforce (3.2), and expanding T by means of edge-element shape functions

Ni’s, i.e. T(r, t) =(cid:80)N

i=1 Ti(t)Ni(r) and

T (t) = s (t) ,

(3.15)

23

J (r, t) =

Ti (t)∇ × Ni (r) .

(3.14)

i=1

Here N is the number of Degrees-of-Freedom (DOF) of the discrete representation. The tree-

cotree decomposition technique is used to impose both the uniqueness (gauge condition) of
the electric vector potential and the boundary condition J · ˆn=0 on ∂D [49, 51].

The Galerkin method applied to the weak form (3.10) yields the linear system

N(cid:88)

(cid:18)

(cid:19)

R+

d
dt

L

where

(cid:90)

D

Rij (cid:44)

(cid:90)
(cid:90)
∇ × Ni · η∇ × NjdV,
(cid:90)

(cid:107)r − r(cid:48)(cid:107)

∇ × Ni (r) · ∇ × Nj

D

D

∇ × Ni (r) · Ac (r, t) dV.

D

Lij (cid:44) µ0
4π
si (t) (cid:44) − d
dt

(cid:0)r(cid:48)(cid:1)

dV dV (cid:48),

It is worth noting that [38]

(cid:90)

D

∇ × Ni(r) · Ac(r, t)dV =

φeddy
k

(t) =

(cid:88)
(cid:88)

k

i

Mikik(t),

MikTi(t),

or in matrix form

where

φeddy(t) = MT T,

(cid:90)

D

Mik =

∇ × Ni (r) · Ac,k (r) dV.

Thus, (3.15) can also be written as

(cid:18)

(cid:19)

R+

L

d
dt

T (t) = −M

d
dt

i(t),

(3.16)

(3.17)

(3.18)

(3.19)

(3.20)

(3.21)

(3.22)

(3.23)

where i is the column vector containing the unknown coefficient ik’s of impressed current in

the coil array.

Substitute (3.20) into (3.7), the output voltage across the k-th coil when the excitation

24

current is switched off is

or in matrix form

(cid:80)

i MikdTi(t)

dt

,

veddy
k

(t) =

veddy(t) = MT dT(t)
dt

.

3.1.3 Natural Modes for Eddy Current Problems

(3.24)

(3.25)

Following [50], the natural modes are the homogeneous solutions of (3.15):

(cid:18)

(cid:19)

d
dt

R+

L

T (t) = 0,

(3.26)

and, as well known, are of the type ue−t/τ where u satisfies:

(cid:19)

(cid:18)

− 1
τ

L + R

u = 0.

(3.27)

τ and u are, therefore, solutions of the following generalized eigenvalue problem:ii

Lu = τ Ru.

(3.28)

The τ(cid:48)s are the time constants of the first order differential system (3.15). The eigenvectors
{u1, ..., uN} form a basis that is orthogonal w.r.t. the scalar product induced by L and R

i Luj=0 for i (cid:54)= j,
uT
i Ruj=0 for i (cid:54)= j,
uT

(3.29)

(3.30)

iiWe highlight that both L and R are positive definite if uniqueness of the electric vector potential has

been enforced.

25

i.e. the matrix U = [u1, ..., uN ] diagonalizes L and R:

UT LU=DL,

UT RU=DR,

(3.31)

(3.32)

where DL = diag (l1, ..., lN ), DR = diag (r1, ..., rN ) and li > 0 and ri > 0, i = 1, ...N . The

generalized eigenvalues τi’s can be expressed as τi = li/ri.

Let us assume that the generalized eigenvalues are ordered in a decreasing sequence
τ1 ≥ τ2 ≥ ··· ≥ τN > 0, we have the following min-max characterization (See Appendix A):

(cid:32)

(cid:33)

τi = max

dim(U )=i

min
x∈U

xT Lx
xT Rx

.

(3.33)

3.2 Monotonicity of Time Constants

The source free response, i.e. the response when the source current is switched off, of pulsed

eddy current problem can be represented as sum of exponential terms:

∞(cid:88)

J(r, t) =

ciJi(r)e−t/τi,

(3.34)

i=1

where J(r, t) is the eddy current density, as a function of spatial coordinate r and time t

and τi is the time constant as obtained by solving (3.28). τi’s are real, positive, bounded

and approach to zero. According to the Biot-Savart law, the corresponding magnetic flux

26

(cid:90)

density is

Beddy(r, t) =

=

=

µ
4π

∞(cid:88)
∞(cid:88)

i=1

J(r, t) × (r − r(cid:48))dV (cid:48)

(cid:90)

D

D

ci

µ
4π

(cid:107)r − r(cid:48)(cid:107)3
Ji(r) × (r − r(cid:48))dV (cid:48)

(cid:107)r − r(cid:48)(cid:107)3

ciBeddy

i

(r)e−t/τi.

i=1

Adopting (3.7) and (3.9), the voltage across the k-th coil is

e−t/τi

(3.35)

vk(t) =

=

=

=

(cid:90)
(cid:90)

(t)

(cid:90)

d
dt

dφeddy
k
dt
µ0
4π

(cid:32)
∞(cid:88)
∞(cid:88)

(cid:32)

i=1

ci

γk
µ0
4π

γk
civk,ie−t/τi.

(cid:33)

D

(cid:90)

dl(cid:48) · J(r, t)
(cid:107)r(cid:48) − r(cid:107) dV
dl(cid:48) · Ji(r)
(cid:107)r(cid:48) − r(cid:107) dV

D

(cid:33)

· −1

τi

e−t/τi

(3.36)

i=1

It is worth notice that, these different quantities share the same set of time constants.

This is straightforward according to Biot-Savart law and Faraday’s law of induction. From a

practical point of view, it states that the time constants can be measured exploiting different

sensors and techniques. It can be achieved with any standard pulsed eddy current testing

system. The only requirement is that the response should be measured when the excitation

current is switched off, which stresses the idea that the excitation current is impressed and

controlled by the measurement system.

27

3.2.1 Monotonicity Principle

Given a conductive domain D and two different anomalies V1 and V2, we assume that the

resistivity of the anomalies is larger than that of the background. The time constants related

to anomaly V1 are denoted as τ

V1
i

V2
and those related to anomaly V2 as τ
i

, i = 1, 2, 3, ..., N .

The time constants are assumed to be ordered in a decreasing sequence. The monotonicity

of time constants is stated as [50]:

V1 ⊆ V2 ⇒ τ

i ≥ τ
V1

V2
i

, ∀i.

(3.37)

To prove (3.37), we notice that V1 ⊆ V2 implies η1 (r) ≤ η2 (r), and therefore R1 − R2 ≤

0iii. Then from xT R1x ≤ xT R2x we have

xT Lx
xT R1x

≥ min
x∈U

xT Lx
xT R2x

min
x∈U

(cid:32)

(cid:33)

(cid:32)

(3.38)

(cid:33)

= τ

V2
i

.

(3.39)

and, from the min-max characterization (3.33), it turns out that:

τ

V1
i = max

dim(U )=k

xT Lx
xT R1x

min
x∈U

≥ max

dim(U )=k

xT Lx
xT R2x

min
x∈U

3.2.2 Isometry and Time Constants

Isometry is a distance-preserving bijective mapping between two metric spaces. Formally, let

X and Y be two spaces with metrics, i.e. distance function between two elements of a set,
denoted as dX and dY , A map f : X → Y is called an isometry if the following condition is
(cid:12)(cid:12)2dV = xT R2x, ∀x
iiiWe notice that xT R1x =(cid:82)

(cid:12)(cid:12)2dV ≥(cid:82)

k xk∇ × Nk

(cid:12)(cid:12)(cid:80)

D η1

k xk∇ × Nk

(cid:12)(cid:12)(cid:80)

D η2

28

true:

dY (f (a), f (b)) = dX (a, b) for any a, b ∈ X.

(3.40)

Isometry is by definition a bijective, i.e. one-to-one, mapping. By contradiction, if a
map g : X → Y maps two different elements a, b ∈ X to the same element c ∈ Y , then

dY (g(a), g(b)) = dY (c, c) = 0 < dX (a, b). Thus, g is not an isometry.

Any translation, reflection and rotation and any combinations of them are isometries on

Euclidean spaces. Figure 3.2 shows an example of isometries in a 2D Euclidean space.

Figure 3.2: Isometry in Euclidean space: translation, reflection and rotation.

Time constants depend on properties of the conductive specimen only. Two configura-

tions that can be transformed into each other by an isometry are said to be congruent. The

time constants of two congruent configurations are identical. For instance, let us consider a
20 mm×20 mm × 20 mm cube having a resistivity of 1× 10−7 Ω m and two possible anoma-

lies as shown in Figure 3.3(a). These two configurations are identical up to a reflection

about the x = 0 plane, which is an isometry and, their related time constants are identical

as shown in Figure 3.3(b). This means that, in principle, two congruent configurations are

indistinguishable from the knowledge of time constants.

However, it is possible to overcome this issue by introducing additional conductive patches

external to the specimen, as part of the inspection system. It is required that the original

isometry which relates two congruent configurations is not equivalent to an identity map for

29

(a)

(b)

Figure 3.3: If two configurations can be transformed into each other by an isometry, they
have identical time constants. (a) Schematic of two isometric configurations, and (b) the
related time constants.

these additional conductive patches. When source currents are impressed in the excitation

coils, eddy currents are induced not only in the specimen but also in the external patches

and hence, the eddy current distributions are different even for two congruent anomalies.

This idea exploits the freedom in the design of inspection system as we have control not only

on the probes and excitation waveforms, but also on the problem domain setup.

For example, let us introduce a conductive patch at 1 mm from the specimen. As shown

in Figure 3.4(a), these two configurations are no longer congruent related by reflection about

the x = 0 plane. The thickness of the patch is 3 mm and the resistivity of the patch is
1.25 × 10−7 Ω m. The resistivity of the additional patch should be chosen smartly. As a

general rule, the resistivity of the patch and of the specimen should be of the same order of

magnitude, but the optimal value is yet to be determined depending on the size of the patch,

on the distance between them, etc. If the resistivity is chosen too large, the eddy currents

mainly circulate in the specimen and the patch does not introduce enough interference in the

eddy current path and, eventually little disturbance in the time constants. If the resistivity

is chosen too small, the eddy currents mainly circulate in the patch and the time constants

30

reveal the properties more of the patch rather than of the specimen, the useful information

in time constants is little.

(a)

(b)

Figure 3.4: Introducing an additional conductive patch can help distinguish two isometric
configurations. (a) Schematic of the test specimen together with the patch, and (b) the
related time constants.

The time constants related to the two anomalies after the introduction of additional patch

are given in Figure 3.4(b). From this plot it is evident that the two anomalies can be clearly

distinguished.

Summing up, when introducing the additional conductive patches: (i) the isometry which

relates two congruent configurations can not be equivalent to an identity map for the patches,

(ii) the shape, position and resistivity of the patches should be chosen smartly in order to

interfere properly with the eddy current in the specimen.

3.2.3 Interaction between Natural Modes and Anomalies

In this section we investigate how anomalies may affect the time constants. Let us consider
a simple example consisting of a conductive cube of 20 mm×20 mm×20 mm. The resistivity
of the cube is 1× 10−7 Ω m. By solving the generalized eigenvalue problem (3.28), we obtain

the time constants τi, the generalized eigenvector ui and the associated natural modes in

31

terms of eddy current density Ji(r) =

N(cid:88)

ui[k]∇ × Nk (r), where ui[k] is the k-th element

of ui. Figure 3.5 plots the natural modes in terms of Ji for some selected i.

k=1

(b) i = 1

(c) i = 2

(d) i = 3

(e) i = 4

(f) i = 5

(a) Conductive cube

(g) i = 6

(h) i = 7

(i) i = 8

(j) i = 9

(k) i = 10

(l) i = 11

Figure 3.5: The natural modes of a conductive cube. (a) The cube, and (b)-(l) the natural
modes associated with the k−th time constant.

As depicted in Figure 3.5, the eddy currents related to the 1st, 2nd and 3rd time con-

stants circulate mainly on the surfaces of the cube. As one can expect, these three modes

correspond to time constants of identical value since they are equivalent up to a proper ro-

tation. Similarly the eddy current modes i = 4, 5, i = 6, 7, 8 and i = 9, 10, 11 are three other

32

groups of equivalent modes. The natural modes in the same group correspond to identical

time constants.

When there is an anomaly in the cube, the eddy current is forced to redistribute around

the anomaly. The anomaly affects a mode if and only if it interacts significantly with the

associated eddy current. For instance, if an anomaly is located in the center of the cube,

then it will not interacts significantly with modes like the 1st, 2nd and 3rd which circulate

mainly close to the surface. On the other hand, the same anomaly will interact significantly

with the 6th, 7th and 8th modes modifying the corresponding time constants.

As a comparison, a surface breaking anomaly on the top surface and a buried anomaly

in the center of the specimen were considered and the time constants were computed with
and without the anomaly. The dimensions of the anomaly are 4 mm×4 mm×4 mm. The

surface breaking anomaly interacts mainly with modes i = 2, 3 and i = 10, 11; the buried

anomaly interacts mainly with modes i = 6, 7, 8. The corresponding results are summarized

in Table 3.1 and Table 3.2.

Table 3.1: The perturbation of time constants when introducing a surface breaking anomaly

Anomaly No. No anomaly (µs) With anomaly (µs) Difference (µs)

1
2
3
4
5
6
7
8
9
10
11

178
178
178
96.2
96.2
89.8
89.8
89.8
78.2
78.2
78.2

178
172
172
96.0
96.0
89.6
89.6
88.4
77.9
74.2
74.2

(a)

0
6
6
0.2
0.2
0.2
0.2
1.4
0.3
4
4

In Table 3.1 a surface breaking anomaly is introduced at the center of top surface. Modes

33

Table 3.2: The perturbation of time constants when introducing a buried anomaly

Anomaly No. No anomaly (µs) With anomaly (µs) Difference (µs)

1
2
3
4
5
6
7
8
9
10
11

178
178
178
96.2
96.2
89.8
89.8
89.8
78.2
78.2
78.2

177
177
177
96.2
96.2
81.6
81.6
81.6
78.2
78.2
78.2

(b)

1
1
1
0
0
8.2
8.2
8.2
0
0
0

i = 2, 3 interact significantly with this anomaly and, therefore, the associated time constants

are strongly affected, whereas the time constant for mode i = 1 is barely changed. Similarly

for another group of modes i = 9, 10, 11 which circulate close to the surface. Modes i = 4, 5

form a loop on the surface and hence barely interact with the anomaly in the center of this

loop. Modes i = 6, 7, 8 circulate mainly in the center of the cube and the associated time

constant are only slightly affected by the presence of the surface breaking anomaly.

In Table 3.2 a buried anomaly is introduced at the center of the cube. Modes i = 6, 7, 8

interact significantly with this anomaly and, therefore, the corresponding time constants are

strongly reduced. Modes i = 1, 2, 3, i = 4, 5 and i = 9, 10, 11 circulate mainly on the surface

of the cube and the associated time constant are only slightly changed by the presence of

buried anomaly.

Concluded from these simple examples, time constants associated with modes circulating

mainly in the innermost region of the conductive specimen are sensitive to buried anomalies.

Therefore, the analysis of time constants is a valuable tool for inspecting deep regions of the

specimen which is usually a difficult task with conventional eddy current methods.

34

3.3 Application of MPM to Imaging

Monotonicity of time constants yields a fast and simple imaging method, as presented in

[45, 46, 47].

3.3.1 Imaging Method based on Monotonicity Principle

Let Tj, j = 1, 2, ..., P , be the test elements. That is, Tj is an a priori prescribed volume
inside the specimen. And let V be the unknown anomaly. From (3.37), if Tj ⊆ V , then
i ≥ τ V
Tj
i , then Tj (cid:42) V , i.e.
Tj is completely or partially external to V (Figure 3.6). By repeating this simple test with

i , ∀i. Thus, if there exists at least one i such that τ

Tj
i < τ V

τ

different Tj we may estimate V as follows:

(cid:91){Tj|τ

VU =

i ≥ τ V
Tj

i ,∀i}.

(3.41)

(a)

(b)

Figure 3.6: The imaging method exploiting the union of test elements.

(a) If the test

element Tj is part of the anomaly V , then τ

i , for all i. (b) If there exist an i such

i ≥ τ V
Tj

that τ

Tj
i < τ V

i , then Tj is not part of V .

Similarly, if V ⊆ Tj then τ V

i ≥ τ

Tj
i

, ∀i. Thus, if there exists at least one i such that

35

Tj
τ V
i < τ
i

then V (cid:42) Tj, i.e. V is completely or partially external to Tj (Figure 3.7). By

repeating this test with different Tj we may estimate V as follows:

(cid:92){Tj|τ V

VI =

i ≥ τ

,∀i}.

Tj
i

(3.42)

(a)

(b)

Figure 3.7: The imaging method exploiting the intersection of test elements.
anomaly V contains the test element Tj, then τ V
such that τ V
, then V does not contain Tj.

i ≥ τ

Tj
i

Tj
i < τ
i

, for all i. (b) If there exists an i

(a) If the

VU and VI are the anomaly profiles estimated via the union/intersection of the test

elements, respectively.

We also point out that, when the unknown anomaly is precisely represented by the union

of the test elements:

(cid:91)

Tj2

(cid:91)···(cid:91)

TjM

V = Tj1

(3.43)

where {Tjm|m = 1, 2, ..., M} is a subset of all test elements, the reconstruction with the

union method is equal or larger than the unknown anomaly. According to (3.41):

36

Tjm ⊂ V ⇒ τ

Tjm
i

≥ τ V

test elements V = Tj1

V ⊂ Tjm ⇒ τ V

i ≥ τ

Tjm
i

i ,∀i ⇒ Tjm ⊆ VU ⇒(cid:91){Tjm} ⊆ VU ⇒ V ⊆ VU .
(cid:84)···(cid:84) TjM
,∀i ⇒ VI ⊆ Tjm ⇒ VI ⊆(cid:92){Tjm} ⇒ VI ⊆ V.

, we have

(cid:84) Tj2

(3.44)

(3.45)

Similarly, when the unknown anomaly is precisely represented by the intersection of the

In other words, combining the union and the intersection method, we can reconstruct

the upper and lower bounds of the unknown anomaly.

3.3.2 Consideration of Noise

In practice, the measured waveforms are noisy, which in turn affects the estimated time

constants. Here, as a preliminary step, we assume the noise affecting the time constants to

be a multiplicative one.

Let us denote the noise-free time constants for V as τ V

version as ˜τ V
i . We assume a multiplicative noise, i.e. ˜τ V
uniformly distributed random variable in the interval [−∆, ∆].

i and the corresponding noisy
i = τ V

i (1 + ξi) where ξi is a

Adapting the idea of Harrach [52], if a test element Tj is part of V , we have τ

which corresponds to τ

i ≥ ˜τ V
Tj
1+ξi

i

,∀i. Therefore, ˜τ V

i ≤ τ

Tj
i

(1 + ξi) ≤ τ

Tj
i

i ≥ τ V
Tj

i ,∀i
(1 + ∆),∀i and the

reconstruction algorithm based on union becomes:

(cid:91){Tj|τ

Tj
i

VU =

(1 + ∆) ≥ ˜τ V

i ,∀i}.

(3.46)

37

Similarly, in the presence of noise, the reconstruction algorithm based on intersection

becomes:

VI =

(cid:92){Tj|˜τ V

i ≥ τ

Tj
i

(1 − ∆),∀i}.

(3.47)

3.4 Concluding Remarks

As discussed through this chapter, time constants is a very powerful tool in eddy current

testing. Given the comprehensive knowledge of natural modes and time constants, we are

able to reveal the defect profiles even at the innermost region of test sample, which is not

possible with traditional eddy current methods. However, time constants is not a directly

measurable quantity.

It is of great significance to correctly and accurately estimate time

constants from the time-domain eddy current response. Next chapter will be dedicated to

this problem, exploring different methods to estimate time constants, and discussing the

merits and limitations of different approaches.

38

Chapter 4

Estimation of Time Constants

In the following section, we assume that the measurements consist of the terminal voltages

of a set of Nc pick-up coils due to eddy currents only. Adopting the expression of (3.36), we

denote the voltage across the k-th coil as

vk(t) =

civk,ie−t/τi =

N(cid:88)

i=1

N(cid:88)

i=1

βk,ie−t/τi, k = 1, 2, ..., Nc.

(4.1)

The problem of estimating time constants is to estimate the number of exponential terms

N , the time constants τi and the linear coefficients βk,i given the knowledge of vk(t) at a

finite number of discrete samples at tj = j∆t, assuming uniform sampling.

4.1 Previous Works on Exponential Analysis

edge of v(t) = (cid:80)N

The estimation of time constants from transient waveform, i.e. finding τi’s given the knowl-
i=1 βie−t/τi, is challenging. A more general framework of this problem

is recognized as the exponential analysis, which is proved to be highly ill-posed when the

number of exponential terms is large or when two consecutive time constants are close to

each other [53]. This problem has been tackled from the mathematical, physics and engi-

neering perspectives by previous researchers, but none of these methods yielded a general

solution to it. The major numerical techniques of exponential analysis include: the nonlinear

39

least square fit, the Prony method [54, 55], the method of modulating functions [56, 57], the

method of moments [58], the Laplace-Pad´e approximation [59], the Tikhonov regularization

method [60], the Gardner transformation [61, 62], the method of maximum entropy [63, 64]

and others.

Most of these methods, similar to the Fourier transform approach, rely on generating

peaks in some spectral function for the detection of the exponential components with real

exponents. The general scope of these methods corresponds to numerically performing the

inverse Laplace transform of v(t). Hence, they suffer from difficulties of unresolved peaks in

the noisy data [59]. It can be proved that the closest exponential components that can be

resolved satisfy [53]:

δ =

τi+1
τi

= eπ/ωmax

(4.2)

where ωmax is determined by the SNR (signal-to-noise ratio) in the transient waveform:

cosh(πωmax) = π(SN R)2

(4.3)

The ill-posed nature of the exponential analysis problem was demonstrated by the exam-

ple of Lanczos [65]. He showed that a sum of two exponentials could be reproduced to within

two decimal places by a sum of three exponentials with entirely different time constants and

amplitudes. This example of Lanczos is reproduced in Figure 4.1, where the experiment data

(circles) and the double-exponential and triple-exponential fits (lines) are plotted. Though

in Figure 4.1 there are two separate lines, they are not distinguishable because the difference

between them is smaller than the line width.

Of course, it is important to note that a priori knowledge can partly compensate for the

40

Figure 4.1: The example of Lanczos. Twenty-four data points (circles) are fitted by a
double-exponential function f1(t) = 2.202e−t/4.45 + 0.305e−t/1.58 (blue solid line) and by a
triple-exponential function f2(t) = 0.0951e−t + 0.8607e−t/3 + 1.5576e−t/5 (red dash line).

information lost in noisy data, and can be used to improve the accuracy. For instance, if

it is known that all linear coefficients βi are positive, specific algorithms could provide fits

orders of magnitude more accurate than other generic algorithms [66, 67]. Other useful a

priori information includes the range of τ , the number of exponential terms, etc.

Istratova and Vyvenko [53] emphasized three cases of exponential analysis. In the simplest

case, further referred to as ”monoexponential analysis”, the transient is assumed to be a

single exponential, which is characterized by the linear coefficient A and exponential λ:

f (t) = Ae−λt

(4.4)

If the decay consists of a sum of N exponentials of the form (4.4), it is within the scope

of ”multiexponential analysis”:

N(cid:88)

i=1

Aie−λit

f (t) =

(4.5)

The goal of the multiexponential analysis is to determine the number of exponential

41

components N , their amplitudes Ai, and decay rate λi. Finally, in the general case when

the decay is described by a continuous spectral function g(λ) rather than by a sum of

discrete exponential transients, it is discussed as the analysis of ”nonexponential transients”

or spectroscopic methods of exponential analysis:

(cid:90) ∞

f (t) =

g(λ)e−λtdλ

(4.6)

0

Analysis of nonexponential transients is aimed at determining the spectral function g(λ).

Equation (4.6) reduces to (4.5) if the spectral function g(λ) can be represented as a sum of

N delta functions:

N(cid:88)

i=1

g(λ) =

Aiδ(λ − λi)

(4.7)

It is thus straightforward that the problem of estimating time constants can be tackled

with those techniques to solve the problem of ”multiexponential analysis” and ”nonexponen-

tial analysis”. An extensive review of these techniques from various branches of engineering,

physics and mathematics is given in [53]. However, when N is a large integer, most of these

techniques become unstable and the solution is not unique.

Despite all challenges as discussed above, it is important to mention that previous studies

on this problem focused on one single waveform only. Whereas in our study of estimating

time constants, there are multiple waveforms sharing the same set of time constants. It is

potential that by exploiting this extra information, an properly-designed algorithm could

enhance the signal to noise raito (SNR) and eventually, increase the resolution in time

constants space, i.e. distinguish two ”close” time constants.

This study is currently still under active research. The following sections will present a

modified method based on the Laplace-Pad´e approximant approach for accurately estimating

42

the time constants and preserving the monotonicity property.

4.2 A Modified Method Based on the Laplace-Pad´e

Approximant Approach

This report introduces a modified Laplace-Pad´e approximant approach to estimate the time

constants.

The Laplace-Pad´e approximant approach was first introduced by Yeramian and Claverie

in [59]. This approach is based upon the combined use of Laplace transform and of Pad´e

approximants. It does not require a hypothesis of the number of exponential terms, which is

an output of the analysis. Further, this is an iterative method that allows a priori knowledge

to be integrated into the stopping criteria.

The proposed method is based on the Laplace-Pad´e approximant approach, tailored for

the pulsed eddy current problem.

Let us assume the voltage measured across the k−th coil,

N(cid:88)

vk(t) =

βk,ie−t/τi + ξk(t),

(4.8)

i=1

where ξk(t) is the noise signal.

In practice, we have at our disposal only a finite number of samples of vk(t) at tj =

j∆t, j = 0, ..., M , where ∆t is the time step and tM should be large enough such that all the

waveforms decays to almost zero. We assume that all individual waveforms vk(t) have been

normalized to their maximum values.

43

By definition, the Laplace transform of vk(t) is

(cid:90) ∞

0

Vk(s) =

vk(t)e−stdt.

(4.9)

Specifically, the Laplace transform of a exponential decay function f (t) = e−αt, α > 0 is

F (s) = 1

s+α , whose region of convergence (ROG) is Re(s) > −α. Thus,

N(cid:88)

Vk(s) =

s + 1/τi

i=1

βk,i

+ Ξk(s), Re(s) > −1/τ1.

(4.10)

where Ξk(s) is the Laplace transform of ξk(t) and τ1 > τ2 > ... > τN .

The Pad´e approximant is the approximation of a function as a ratio of two power series.
Given a function f and two integers L ≥ 0 and M ≥ 1, the Pad´e approximant of order

[L/M ] is the rational function:

(cid:80)L
1 +(cid:80)M

RL/M (x) =

PL(x)
QM (x)

=

l=0 alxl
m=1 bmxm

=

a0 + a1x + a2x2 + ... + aLxL
1 + b1x + b2x2 + ... + bM xM

.

(4.11)

Notice that the first term of (4.10) is exactly its own Pad´e approximant of order [N−1/N ].

As a general remark, the Pad´e approximant is unique for given L and M . Especially for

a power series

The coefficients are found by setting

∞(cid:88)

j=0

djxj

A(x) =

A(x) − PL(x)
QM (x)

= 0,

44

(4.12)

(4.13)

which gives the set of equations

d0 = a0

d1 + d0b1 = a1

d2 + d1b1 + d0b2 = a2

...

dL + dL−1b1 + ... + d0bL = aL

dL+1 + dLb1 + ... + dL−M +1bM = 0

...

dL+M + dL+M−1b1 + ... + dLbM = 0.

(4.14)

The Laplace-Pad´e approximant algorithm attempts to equate the Taylor expansion of
Vk(s) at point s0 up to order 2n − 1 and the Pad´e approximant of Vk(s) of order [n − 1/n],

where n is estimate number of exponential terms. n changes during the execution of the

v(t) =(cid:80)N

algorithm and is essentially an output of the algorithm. Assuming only one single waveform

i=1 βie−t/τi + ξ(t), the original algorithm for estimating the time-constants based

on the Laplace-Pad´e approximation [59] is summarized as follows:

Step 1 Set the estimated number of time constants n to 1;

Step 2 Calculate the coefficients of Taylor expansion for the Laplace transform of v(t) at

point s0 and its successive derivatives up to order 2n − 1;

Step 3 Equate the Taylor expansion and the Pad´e approximant of order [n − 1/n] at point

s0:

2n−1(cid:88)

i=0

di(s − s0)i =

a0 + a1x + a2x2 + ... + an−1xn−1

1 + b1x + b2x2 + ... + bnxn

;

(4.15)

45

Step 4 Find the poles and residues of the Pad´e approximant, which correspond to the

estimated time constants ˆτi and linear coefficients ˆβi, respectively;

Step 5 If a linear coefficient ˆβi is significantly smaller than the remaining ones, the algo-

rithm stops;

Step 6 Set n = n + 1 and go to step 2.

We notice that the Laplace transform referred to in Step 2 may be numerically computed

from the samples of the waveform as follows:

1

2

(cid:16)

V (s) ≈ ∆t

e−st1v(t1) + e−stM v(tM )

 .

−stj v(tj)

e

(cid:17)

+

M−1(cid:88)

j=2

The Taylor expansion of V (s) at s = s0 is

+∞(cid:88)

V (s) =

di(s − s0)i,

i=0

. Adopting the expression of (4.16), the derivatives of V (s) w.r.t

(4.16)

(4.17)

(cid:18)

d(i)V
ds(i)

where di = 1
i!

(cid:19)

s=s0

s is numerically calculated as

(cid:33)

(cid:32)
1

2

di =

1
i!

d(i)V
ds(i)

≈ ∆t

(cid:16)

s=s0

(cid:17)
(−t1)ie−st1v(t1) + (−tM )ie−stM v(tM )

+

M−1(cid:88)

j=2

 .

(4.18)

(−tj)ie

−stj v(tj)

Also, adopting (4.14) in Step 3, the coefficients of the denominator bi can be obtained

46

by solving an set of n equations with n unknowns:



dn−1

. . . d1

dn
...

. . . d2
...
. . .

d0

d1
...

d2(n−1)

. . . dn dn−1





b1

b2
...

bn

 =



−dn
−dn+1

...

−d2n−1

 .

(4.19)

Then the coefficients of the numerator ai are

ai = di + di−1b1 + ··· + d0bi(i = 0, 1, ..., n − 1).

(4.20)

The algorithm proposed in this study is modified in the following manner. First, multiple

waveforms sharing the same set of time constants are considered. Multiple waveforms are

simultaneously integrated in the solution to enhance the robustness w.r.t. noise. When
trying to equate the Taylor expansion of order 2n − 1 and the Pad´e approximant of order
[n − 1/n] at point s0, in the original algorithm, there are n independent equations and

n unknowns; in the algorithm proposed here, there are Ncn independent equations and n

unknowns, which gives an over-determined equation set.

[b1, b2, ..., bn]T are determined in

the least-squares sense.

Second, the selection of s0, which is the only parameter to be set in reconstruction

procedure, was originally chosen by experience. Specifically, [68] states that a reasonable

choice for s0 is the inverse of the time the waveform takes to decay to half of its initial value,

hereafter denoted as s1/2. In the proposed algorithm, s0 is determined by minimizing the

47

following function:

(cid:118)(cid:117)(cid:117)(cid:117)(cid:116)(cid:88)

k,j

[vk(tj) − N(cid:88)

i=1

ϕ(s0) =

ˆβk,i(s0)e

−tj /ˆτi(s0)]2

(4.21)

where ˆβk,i(s0) and ˆτi(s0) are the estimations of βk,i and τi, as obtained in Step 4 of the

original algorithm. ˆβk,i(s0) and ˆτi(s0) depend on the point of expansion s0.

In numerical fits, residuals, i.e. the difference between the experimental data and the

fit, are often used to evaluate the quality of the fit. Equation (4.21) minimizes the sum

of the squares of the residuals at every time step. By implementing (4.21) it eliminates

the arbitrariness in the choice of s0. Some preliminary tests about s0 were carried out in

this study and we assumed that the s0 which minimize the l2-norm of the residual is the

best choice. Test results showed that s1/2 is usually close to the optimal choice but not the

optimum. Therefore, in this study the choice of s0 is determined by minimizing the objective

function ϕ(s0) in the range [0.25s1/2, 4s1/2].

Third, a proper stopping criteria is proposed considering the physics of pulsed eddy

current problem. The algorithm stops when at least one of the following conditions is met:

• A negative time constant appears;

• A pair of conjugate complex time constants appear;

• The error reaches a minimum;

• The error is smaller than the prescribed noise level.

A step-by-step procedure of the proposed algorithm is given as follows:

Step 1 Set the estimated number of time constants n to 1;

48

Step 2 Set the point of Taylor expansion s0 to s1/2;

Step 3 Calculate the coefficients of Taylor expansions for the Laplace transform of vk(t), k =

1, 2, ..., Nc at point s0 and their successive derivatives up to order 2n − 1;

Step 4 Equate the Taylor expansions and the Pad´e approximant of order [n− 1/n] at point

s0;

Step 5 Find the poles and residues of the Pad´e approximant, which correspond to the

estimated time constants ˆτi(s0) and linear coefficients ˆβk,i(s0), respectively;

Step 6 Update s0 to minimize the error function (4.21) in the interval [0.25s1/2, 4s1/2].

Assuming the error function ϕ(s0) is a unimodal function, this is implemented as a

golden-section search [69].

Step 7 If ϕ(s0) is the minimum, let s∗ = s0 and go to Step 8, otherwise go to Step 3;

Step 8 If one of the following stopping criteria is met, the algorithm stops:

• one of the estimated time constants (at s∗) is either negative or complex;

• the error, as a function of n, increases;

• the error is smaller than a prescribed threshold related to the noise level.

Step 9 Set n = n + 1 and go to step 2.

49

4.3 Numerical Examples on Estimating Time

Constants

This section provides numerical examples of application of the aforementioned algorithm on

synthesized signals consisting of the sum of exponential terms. In the framework of MPM,

it is important not only to accurately estimate the time constants, but also to preserve

the inherent monotonicity of time constants. In other words, if anomaly A is contained in

anomaly B, all the time constants related to A, when ordered in a decreasing sequence,

have to be greater than the corresponding time constants related to B. The time constants

estimated via the proposed algorithm from measurement of A must be greater than those

from the measurement of B. If there exists τ A

i > τ B
i

for some i and τ A

k < τ B

k for some k,

we say that the time constants related to A and B have cross-overs, which means neither A

is part of B nor B is part of A. For a good algorithm which estimates the time constants,

if the true time constants related to anomaly A and B have cross-overs, then the estimated

time constants have cross-overs as well.

To demonstrate the validity of the proposed algorithm, two numerical examples are pre-

sented in this section.

Example 1:

Group 1

v1(t) = e−t/10 + e−t/5 + e−t/2 + e−t/1 + e−t/0.5 + e−t/0.2 + e−t/0.1
v2(t) = e−t/10 + 2e−t/5 + 3e−t/2 + 4e−t/1 + 5e−t/0.5 + 6e−t/0.2 + 7e−t/0.1
v3(t) = 7e−t/10 + 6e−t/5 + 5e−t/2 + 4e−t/1 + 3e−t/0.5 + 2e−t/0.2 + e−t/0.1

50

Group 2

v1(t) = e−t/9 + e−t/4.5 + e−t/1.8 + e−t/0.9 + e−t/0.45 + e−t/0.18 + e−t/0.09
v2(t) = e−t/9 + 2e−t/4.5 + 3e−t/1.8 + 4e−t/0.9 + 5e−t/0.45 + 6e−t/0.18 + 7e−t/0.09
v3(t) = 7e−t/9 + 6e−t/4.5 + 5e−t/1.8 + 4e−t/0.9 + 3e−t/0.45 + 2e−t/0.18 + e−t/0.09

The first example comprises two groups of signals. There are three time-domain signals

in each group. Each of these signals is sum of 7 exponential terms. In other words, these two

groups of signals are considered as in (4.8) with k = 1, 2, 3 and N = 7. The time constants

in group 2 are 90% of those in group 1 and the linear coefficients βk,i are the same in both

groups.

In pulsed eddy current testing, the signals measured at multiple pick-up coils share the

same set of time constants. In order to imitate the experimental measurements, each group

consists of three waveforms representing the voltages measured at three different pick-up

coils.

All waveforms are sampled from 0 to 50s, which is five times of the largest time constant,

and the sampling rate is 1000 samples per second. An additive noise proportional to the

magnitude of each individual waveform is imposed in this study. Specifically, the noise

ξk(t) added to the k-th waveform, is a uniformly distributed random variable onto interval
[−∆ · max{|vk(t)|}, +∆ · max{|vk(t)|}], where ∆ = 0.001 (0.1% noise level). The energy of
k,j ξk(tj)2, is 0.0647 for Group 1 and 0.0649 for

the noise signal, which is defined as

(cid:113)(cid:80)

Group 2.

The estimation results are summarized in Table 4.1. ”n” in the first column represents

the number of exponential terms to be estimated and ”Error” gives the value of the error

51

Table 4.1: Estimation of time constants from synthesized signals in Example 1.

n Group

Estimated Time Constant(s)

1

2

3

4

5

6

7

1
2
1
2
1
2
1
2
1
2
1
2
1
2

3.8819
3.4824
6.2542
5.6453
7.4942
6.7439
8.4740
7.6238
9.0699
8.1189
10.0454
8.7065
14.4586
10.1987

0.7429
0.6723
1.4266
1.2837
2.3258
2.0866
3.2387
2.8636
5.0594
3.9017
11.5321
6.4859

0.2840
0.2556
0.7369
0.6605
1.1743
1.0544
1.9430
1.4457
5.6226
2.9230

0.1785
0.1606
0.4482
0.4070
0.8373
0.5987
2.0509
1.2800

0.1383
0.1254
0.3209
0.2321
0.8739
0.5395

0.1143
0.1016
0.3362
0.1898

0.1171
0.0847

Error
10.155
9.6351
3.0323
2.8787
0.9604
0.9143
0.3173
0.3047
0.1383
0.1406
0.0661
0.0675
0.0687
0.0681

function as defined in (4.21). Starting from n = 1, the proposed algorithm attempts to

estimate n time constants from the transient signals. The value of error function decreases

as n increases, which means the reconstructed signals from estimated parameters, namely ˆτi

and ˆβk,i, become better fits to the original signals. The error function is minimized at n = 6.

Its values, 0.0661 and 0.0675, are very comparable to the energy of the noise signal, 0.0647

and 0.0649. The estimated time constants follow the same monotonicity as in the case of the

true time constants: all time constants in group 1 are larger than the corresponding ones in

group 2.

Example 2:

Group 1

v1(t) = e−t/10 + e−t/5 + e−t/2 + e−t/1 + e−t/0.5 + e−t/0.2 + e−t/0.1
v2(t) = e−t/10 + 2e−t/5 + 3e−t/2 + 4e−t/1 + 5e−t/0.5 + 6e−t/0.2 + 7e−t/0.1
v3(t) = 7e−t/10 + 6e−t/5 + 5e−t/2 + 4e−t/1 + 3e−t/0.5 + 2e−t/0.2 + e−t/0.1

52

Group 2

v1(t) = e−t/9 + e−t/5.5 + e−t/1.8 + e−t/1.1 + e−t/0.45 + e−t/0.22 + e−t/0.09
v2(t) = e−t/9 + 2e−t/5.5 + 3e−t/1.8 + 4e−t/1.1 + 5e−t/0.45 + 6e−t/0.22 + 7e−t/0.09
v3(t) = 7e−t/9 + 6e−t/5.5 + 5e−t/1.8 + 4e−t/1.1 + 3e−t/0.45 + 2e−t/0.22 + e−t/0.09

The second example also comprises two groups of signals. In fact they are the same as

the signals in example 1 except that: the 1st, 3rd, 5th and 7th time constants in group 2

are 90% of those in group 1 while the rest are 110% of those in group 1. In other words,

the true time constants have cross-overs. The energy of the noise signal is 0.0647 for Group

1 and 0.0648 for Group 2. The estimation results of this example are summarized in Table

4.2. The error function is minimized at n = 6, which means the reconstruction signals give

a best fit to the original signals. The estimated time constants of group 1 and group 2 have

cross-overs, as in the case of the true time constants.

It is worth noting that at n = 7, there are a pair of complex values in the estimated

time constants. Imaginary time constants represent sinusoidal components in the waveform,

which can be understood as the attempt to fit the random, bounded noise signal.

4.4 Concluding Remarks

This chapter focused on the problem of estimating time constants from the time-domain

eddy current response. This problem is broadly seen across different physics and engineering

disciplines, and is more generally recognized as the problem of exponential analysis. A lot

of numerical approaches have been suggested to deal with this problem. Among them the

Laplace-Pad´e approximant approach is most suitable for the specific problem of estimat-

53

Table 4.2: Estimation of time constants from synthesized signals in Example 2.

n Group

Estimated Time Constant(s)

1

2

3

4

5

6

7

1
2
1
2
1
2
1
2
1
2
1
2
1
2

3.8819
3.7243
6.2542
5.9852
7.4942
7.0223
8.4740
7.5744
9.0699
7.9509
10.0454
9.2359
14.4586
8.3945

0.7429
0.7036
1.4266
1.2793
2.3258
1.7254
3.2387
2.7509
5.0594
5.8665
11.5321
4.1456

0.2840
0.2508
0.7369
0.5645
1.1743
1.1882
1.9430
1.8783
5.6226
1.3912

0.1785
0.1637
0.4482
0.3724
0.8373
1.0074
2.0509
0.4-0.2i

0.1383
0.1199
0.3209
0.3364
0.8739
0.4+0.2i

0.1143
0.0995
0.3362
0.3299

0.1171
0.0998

Error
10.155
9.8610
3.0323
2.8069
0.9604
0.2686
0.3173
0.3047
0.1383
0.1389
0.0661
0.0707
0.0687
0.0746

ing time constants, because it does not require the a priori knowledge of number of time

constants as an input, instead, the proper number of exponential terms is an output of the

analysis.

In this study, the Laplace-Pad´e approximant approach was further modified to fit the

problem of estimating time constants. To be more specific, the following modifications were

made: (1) multiple waveforms sharing the same set of time constants are considered; (2) the

value of s0 is determined by minimizing the error function; (3) the estimated time constants

are limited to real numbers only. The modified algorithm has been applied to synthesized

waveforms and was capable of estimating up to 7 time constants with acceptable errors.

54

Chapter 5

Numerical Examples

This chapter demonstrates the main features of the monotonicity principle method by means

of numerical examples. We first provide numerical evidence of the monotonicity of time

constants in a realistic 3D model. Next, examples of imaging via monotonicity of time

constants are provided. A parametric study w.r.t. the noise level and different number of

available time constants is also presented.

Monotonicity of time constants implies that when the size of an anomaly increase, all

related time constants decrease. This property can be applied to the problem of defect

detection, imaging, sizing, etc.

5.1 Validation of Monotonicity

In this section we validate numerically the monotonicity of time constants. A cylindrical

conductive pipe section of length 10 mm, inner radius 5 mm and outer radius 9 mm is

considered. In the finite element model, this conductor is discretized into 30 elements in the

azimuthal direction, 10 elements in the radial direction and 10 elements in the axial direction,

as illustrated in Figure 5.1(a). To break the isometry, three additional conductive patches

are placed close to the inner surface of the pipe as part of the inspection system, see Figure

5.1(b). An array of coils, acting as both transmitters and receivers, are also placed close

to the inner surface of the pipe. The electrical resistivity of the pipe and of the conductive

55

patches is 10−7 Ω m and that of the anomaly is 10−5 Ω m, which is 100 times larger.

(a) Mesh of the discretized model

(b) Additional patches and array of coils

Figure 5.1: 3D mesh of the test structure: a cylinder pipe section

In the first example we assume two outer surface anomalies A and B, where A is contained

in B. Figure 5.2(a) shows the differences between the largest 200 time constants related to A

and B, arranged in decreasing order. As expected, the differences are all positive according

to monotonicity principle.

In the second example we assume two outer surface anomalies B and C, where B and

C are of the same size and are partially overlapped with each other. Figure 5.2(b) shows

the differences between the largest 200 time constant related to B and C, arranged in the

decreasing order. As expected, some of the differences are positive and some are negative

since the monotonicity condition is not satisfied. Similar results are observed in the third

example of an outer surface anomaly B and an inner surface anomaly D, see Figure 5.2(c).

Anomaly B is completely external to anomaly D.

5.2 Defect Reconstruction

In this section we consider the problem of defect imaging in a conductive slab via mono-

tonicity of time constants.

56

Anomaly A

Anomaly B

(a) Anomaly A is part of anomaly B ⇒ τ A

related to anomalies A and B

i ≥ τ B

i , for all i.

Differences between time constants

Anomaly B

Anomaly C

(b) Anomaly B is partially overlapped with anomaly C ⇒

related to anomalies B and C

i ≥ τ C
τ B

i , for some i and τ B

i < τ C

i , for the rest.

Differences between time constants

Anomaly B

Anomaly D

Differences between time constants

(c) Anomaly B is completely external to anomaly D ⇒

related to anomalies B and D

i ≥ τ D
τ B

i

, for some i and τ B

i < τ D
i

, for the rest.

Figure 5.2: Validation of the monotonicity of time constants

The conducting domain is a homogeneous slab of 20 mm×20 mm×3 mm. The region of

interest (ROI) is a thin layer in the middle of the slab: 0< x <0.4 mm, -3 mm< y <3 mm,

-3 mm< z <0, as highlighted in Figure 5.3(a) and enlarged in Figure 5.3(b). We assume
that all potential anomalies are inside the ROI. The resistivity of the plate is 10−7 Ω m and
that of the anomaly is 10−5 Ω m. The ROI is discretized into 9 × 18 voxels in the forward

problem, as shown in Figure 5.3(b). A typical test anomaly in the inverse problem consists

of combinations of these voxels.

The excitation and measurement system consists of a 2D rectangular array of coils as

57

(a)

(b)

Figure 5.3: A homogeneous slab as the test specimen. (a) Finite element mesh. (b) The
region of interest.

seen in Figure 5.4(a). This array is placed at 0.1 mm above the specimen and covers a square
region of 6 mm×6 mm. The dimension of each coil is 1 mm×1 mm and the thickness is

0.1 mm. In order to overcome the issue of isometry, an additional conductive patch, whose
resistivity is 1.25 × 10−7 Ω m, is placed above the specimen, as shown in Figure 5.4(b).

(a)

(b)

Figure 5.4: The inspection system consists of (a) a 2D coil array, and (b) a large conductive
patch above the specimen.

In the following examples all test anomalies consist of the voxels in the ROI, which is

uniformly divided into 1 element in x direction, 18 elements in y direction and 9 elements
in z direction. We use ”a × b” to describe a test anomaly that consists of a consecutive

elements in the y direction and b consecutive elements in the z direction. For simplicity only

test anomalies of square shape are considered. Figure 5.3(b) gives an example of 2 partially

58

overlapped 3 × 3 test anomalies. The target is to reconstruct an anomaly of 5 × 5 voxels. In

this first example we assume that the time constants are noise-free and the first 200 largest

time constants are available.

According to (3.46), we can determine if a test anomaly is part of the unknown target

anomaly by comparing their time constants. Figure 5.5 represents the reconstructed anomaly

using union of test anomalies which satisfy the monotonicity condition. Notice that we are
not presenting the reconstruction results using the union of 6× 6, 7× 7, 8× 8 and 9× 9 test

anomalies since they are empty. These test anomalies are larger than the target anomaly

and thus, the corresponding estimation VU is empty.

Union of 2×2 test anomalies

Union of 3×3 test anomalies

Union of 4×4 test anomalies

Union of 5×5 test anomalies

Figure 5.5: Imaging example using union of test anomalies, assuming noise level ∆ = 0. The
red dashed frames represent the target anomaly.

Similarly, the profile of the target anomaly can also be estimated using the intersection of

test anomalies according (3.47) as shown in Figure 5.6. The results presented in this example

are related to noise-free measurements and the first 200 largest time constants are considered.

Dual to the union method, we have no reconstruction when using test anomalies smaller than

the target anomaly because the set at the r.h.s. of (3.47) is void. The reconstructions using
8 × 8 test anomalies and using 9 × 9 test anomalies have a few extra voxels compared to the

59

Intersection of 5×5 test anomalies

Intersection of 6×6 test anomalies

Intersection of 7×7 test anomalies

Intersection of 8×8 test anomalies

Intersection of 9×9 test anomalies

Figure 5.6: Imaging example using intersection of test anomalies, assuming noise level ∆ =
0. The black bold frames represent individual test anomalies that are found through the
monotonicity test (3.47). The red dashed frames represent the target anomaly.

target anomaly. This is due to the fact that these test anomalies are to large to move freely

in the ROI, especially in z direction, otherwise the intersection area would be reduced.

Also notice that, as a general rule, the reconstruction result is better when the sizes of

test elements are closer to that of the unknown anomaly. In all subsequent discussions, we

provide unions VU using the largest test elements and intersections VI using the smallest

test elements before the reconstruction result vanishes.

In the second example, for a prescribed noise level, we repeated the numerical experiment

for N = 1000 times by generating the time constants using multiplicative noise model (Eq.

(3.46)-(3.47)). The probability of the voxel i being identified as part of the anomaly is defined

as pi =

Ni
N , where Ni is the number of times the voxel is identified as part of the anomaly.

200 largest time constants are used in these examples. Figure 5.7 gives the plots of pi as

60

the average reconstruction results for noise level ∆ = 1%, ∆ = 0.1% and ∆ = 0.01%. The

voxel of a darker color is more probable part of the anomaly. As expected, the reconstruction

result is better with lower noise level. In this experiment setup a proper reconstruction result

is obtained when the noise level is less than ∆ = 0.1%.

Union, ∆ = 1%

Intersection, ∆ = 1%

Union, ∆ = 0.1%

Intersection, ∆ = 0.1%

Union, ∆ = 0.01%

Intersection, ∆ = 0.01%

Figure 5.7: Imaging example with various noise levels. The red dashed frames represent the
target anomaly.

Finally, we considered the impact of the number of available time constants on the re-

constructions. In the previous results, it was assumed 200 time constants are available. In

practice, this number strongly depends on the signal condition as well as the method to

estimate the time constants. In this examples, we reduced the number of available time con-

stants to 150, 100 and 50, and revisit the reconstruction results. The noise level is prescribed

as ∆ = 0.1% and imaging is averaged over N measurements. Results are summarized in

Figure 5.8.

Generally speaking, reconstruction using the union strategy gives an estimate of the

61

Union, 150 time constants

Intersection, 150 time constants

Union, 100 time constants

Intersection, 100 time constants

Union, 50 time constants

Intersection, 50 time constants

Figure 5.8: Imaging example with various numbers of available time constants, assuming
noise level ∆ = 0.1%. The red dashed frames represent the target anomaly.

upper bound of the anomaly while reconstruction using the intersection strategy gives the

estimate of the lower bound. When more time constants are available, the reconstructed

boundaries are more close to the real profile of the target anomaly, and further the difference

between the reconstructions of two strategies is smaller. It is worth noting that, even in the

case where only 50 time constants are available, the reconstructions correctly estimate the

lower and upper bound, which is very reasonable in an experimental realization.

5.3 Application of MPM in Defect Sizing

This section presents the application of MPM to the problem of defect sizing. Time constants

decrease as the size of anomaly increases. Moreover, time constants depend on the position of

anomaly since the way an anomaly interacts with natural modes is related to its position. As

62

a consequence, anomalies of the same size can be clustered via time constants and clusters

related to different sizes are well separated. By investigating if the time constants of an

anomaly fall in a cluster, one can determine its size.

In this examples, the probe comprises an array of 14 coils above an L-shape conductive

block, which is the test specimen, as shown in Figure 5.9. Anomalies of various sizes and

locations are assumed buried in this specimen. More specifically the anomalies are considered
to be cubes of size 2 × 2 × 2 voxels denoted as S (for small), 3 × 3 × 3 voxels denoted as M
(for medium) and 4 × 4 × 4 voxels denoted as L (for large). The excitation coil was driven

with a square wave current of 50% duty cycle. The excitation coils were excited one at a

time and the transient signal was measured at the same coil when the injected current was

switch-off.

Figure 5.9: Illustration of the numerical examples. The coil array is above the test specimen.
The dark colored voxels represent the air domain.

(a)

(b)

(c)

Figure 5.10: Illustration of the anomalies of various sizes. (a) small anomaly (S), (b) medium
anomaly (M) and, (c) large anomaly (L).

The Laplace-Pad´e approximant algorithm as described in section 4.2 was used to ana-

63

lyze the time domain waveforms measured from the coil array. Different noise levels were

superimposed on the measured waveforms. Two time constants can be estimated for each

configuration associated with a specific anomaly profile. Each point in Figure 5.11 represents

an anomaly in the 2D time constant space. The anomalies of three different sizes naturally

group into three clusters and can be separated by linear decision boundaries. It is clearly

illustrated in the figure that configurations associated with larger anomalies have smaller

time constants. Table 5.1 summarizes the probability of anomalies being correctly classified.

For instance, the first row shows that in the case the signals had 0.1% noise, 94.81% of the

small anomalies were correctly classified as small, 4.44% and 0.74% of them were classified

as medium and large, respectively. This result is very encouraging as it shows the proposed
method has a high probability (≥ 90%) of correctly classifying anomalies of different sizes,

even when the noise level is as high as 1%.

(a)

(b)

Figure 5.11: The time constants associated with anomalies of different sizes group into three
clusters. The measured signals have (a) 0.1% noise and (b) 1% noise.

64

Table 5.1: Classification of anomalies via time constants

waveforms have 0.1% noise waveforms have 1% noise

real
size
S
M
L

identified as size

identified as size

S

0.9481

0
0

M

0.0444
0.9643

0

L

0.0074
0.0357
1.0000

S

0.9067
0.0786

0

M

0.0896
0.8929
0.0577

L

0.0037
0.0286
0.9423

5.4 Concluding Remarks

In this chapter, algorithms based on monotonicity of time constants were applied to the

problem of defect reconstruction and defect sizing. Numerical examples were provided first

to validate the monotonicity property in a pipe section, which is a commonly seen structure

in eddy current problems. Next, we applied monotonicity of time constants to the problem

of imaging cracks in a metal plate structure. Numerical examples showed promising results.

However, in all examples with good imaging results, we used at least 50 most significant time

constants, while in practice we can at best estimate less than 10 time constants accurately,

from the analysis of Chapter 4. This raises concerns about applying monotonicity of time

constant in real-life inspections. Luckily, we discovered another nice property in pulsed eddy

current testing, namely the monotonicity of transfer function. Besides, transfer function is

directly measurable in experiments. This topic will be covered in the next two chapters.

The monotonicity of time constants was also applied to the problem of defect sizing. In

this example, only two most significant estimated time constants were used. The classifi-

cation result is very promising as it showed the proposed algorithm can correctly classify

defects of different sizes, with a correct rate of more than 90% even when suffering 1% noise

in the time domain waveforms.

65

Chapter 6

Monotonicity of Transfer Function

and Imaging Method

A transfer function defines the input-output relation of an eddy current system. Specifically

in this chapter, we evaluate the transfer function on the negative real axis and prove that
this transfer function is monotonic with the anomaly profile, i.e. Dα ⊆ Dβ ⇒ Zα − Zβ

is negative semidefinite, where Z is the transfer function. This property is exploited to

implement a non-iterative imaging algorithm. This chapter discusses the resolution of this

algorithm, sensitivity to noise and its implementation with experimental considerations.

6.1 Definition of Transfer Function

As described in Chapter 3, the linear equations that defines an eddy current system can be

written as:

(cid:18)

(cid:19)

R+

L

d
dt

T (t) = −M

d
dt

i(t),

v(t) = MT dT(t)
dt

.

(6.1)

(6.2)

where v(t) represents the voltage induced by the magnetic flux due to eddy current only.
Let us assume the input current is an exponentially decaying one, i(t) = i01(t)e−t/τ +
i01(−t), where i0 is a vector describing the magnitude of currents injected in the coil array,

66

τ is a prescribed positive constant, and 1(t) is the heaviside step function:

0,

1,

1(t) =

t < 0

t ≥ 0

(6.3)

An illustration of the current waveform is given in figure 6.1.

Figure 6.1: An example of the excitation waveform

The solution to (6.1) consists of two parts, a source-free response and a forced response:

T(t) = Tfree(t) + Tforce(t)

(6.4)

The source free response is the solution to the homogeneous problem and can be expressed

as sum of the natural modes:

Tfree(t) =

The forced solution satisfies

(cid:16)

R+ d

dtL

(cid:17)

(cid:88)

k

ckuke−t/τk

(6.5)

T (t) = −M d

dti0e−t/τ and can be expressed as

67

T(t) = Ae−t/τ , where A can be determined by the following deduction:

e−t/τ

Mi0

τ

e−t/τ =
RAe−t/τ − LA
τ
Mi0

⇒RA − LA
τ
⇒τ RA − LA = Mi0
⇒A = (τ R − L)−1Mi0

=

τ

Substituting this into (6.2) gives,

v(t) = MT d
dt

[(τ R − L)−1Mi0e−t/τ ]
MT (τ R − L)−1Mi0e−t/τ

= − 1
τ

(6.6)

(6.7)

A matrix Z, defined as Z = − 1

τ MT (τ R − L)−1M, satisfies the input-output relation
v(t) = Zi(t) = Zi0e−t/τ when the input current is an exponential decaying one. Notice

that Z is not a function of time t. By definition of the Laplace transform, Z is the transfer
function evaluated on the negative real axis s = −1/τ .

The measurement of the transfer function Z is straightforward. The elements of Z are

the self/mutual impedance relates the forced response and the input current. To estimate
Z, one chooses τ >> τ1 ≥ τ2 ≥ ··· ≥ τN and measures the pick-up coil voltages when the

source free response vanishes (falls below the noise floor). To be more specific,

Zmn =

vm(t)
in(t)

, t >> (

)−1

1
τ1

− 1
τ

(6.8)

68

6.2 Monotonicity of Transfer Function and Imaging

Scheme

Given a conductive domain D and two anomalies Vα and Vβ, we assume that the resistivity

of the anomalies is larger than that of the background. The transfer function related to

anomaly Vα is denoted as Zα and that related to anomaly Vβ as Zβ. The monotonicity of

transfer functions is:

Vα ⊆ Vβ ⇒ Zα ≤ Zβ

(6.9)

where Zα ≤ Zβ means that Zα−Zβ is a negative semidefinite matrix, i.e. all of its eigenvalues

are non-positive.

To prove (6.9), we notice that Vα ⊆ Vβ implies ηα(r) ≤ ηβ(r), and therefore Rα ≤ Rβ,

Dα ⊆ Dβ ⇒ R1 ≤ R2

⇒ (τ R1 − L) ≤ (τ R2 − L)
⇒ (τ R1 − L)−1 ≥ (τ R2 − L)−1
⇒ MT (τ R1 − L)−1M ≥ MT (τ R2 − L)−1M[70]
⇒ − 1
τ
⇒ Z1 ≤ Z2

MT (τ R1 − L)−1M ≤ − 1

MT (τ R2 − L)−1M

τ

(6.10)

Notice that the third line of (6.10),(τ R1−L)−1 ≥ (τ R2−L)−1, is true only when matrices
τ R1 − L and τ R2 − L are invertible. By definition of time constants, (τkR − L)uk = 0
for non-trivial u(cid:48)
ks, i.e. (τkR − L) is not invertible for k = 1, 2,··· , N . Because the time
constants τk’s are bounded and approach zero, it is necessary that the decay constant τ of

69

the injected current is chosen to be τ >> τ1 ≥ τ2 ··· .

Similar to monotonicity of time constants, monotonicity of the transfer function can be

exploited to formulate an non-iterative imaging algorithm. Let Tj, j = 1, 2, ..., P , be the test
elements, V be the unknown anomaly. From (6.9), if Tj ⊆ V , matrix ZTj
− ZV is negative
semidefinite. By repeating this simple test with different Tj we may estimate V as follows:

(cid:91){Tj|ZTj

VU =

≤ ZV }.

(6.11)

The proposed method has low computational cost as it involves only a very simply test
− ZV , which is usually small in size: the

on the eigenvalues of the different matrix ZTj

dimension of the matrix is equal to the number of coils in the excitation array. To retrieve

the full profile (VU ), one needs to repeat this simple test for all test elements. The number of

test elements is predefined by the inspector and should be adjusted according to the target

resolution. In practical cases, this number is of the order of hundreds.

When taking into account the noise, only a polluted version ˜ZV of the genuine matrix ZV

is available. In line with [52], we assume that the noise level is known, i.e.

(cid:13)(cid:13)(cid:13)˜ZV − ZV

(cid:13)(cid:13)(cid:13) ≤ δ,

where δ is a measure of the noise affecting the measured data. It is trivial to see the following

inequality holds:

˜ZV − δI ≤ ZV ≤ ˜ZV + δI

In this case, the unknown anomaly profile can be estimated as:

(cid:91){Tj|ZTj

VU =

≤ ˜ZV + δI}.

70

(6.12)

(6.13)

6.3 Numerical Results

This section is devoted to defect reconstruction exploiting the monotonicity of transfer func-

tions. The first part of this section demonstrates the measurement of transfer function from

transient waveforms. The second part presents reconstruction examples of four typical defect

configurations under noise-free and noisy assumptions.

6.3.1 Measurement of Transfer Function

The transfer function can be accurately measured from the transient waveforms of pick-up

coil voltages and input currents. Given a system described by (6.1) and (6.2), the θ-method

evaluates the transient response as:

T(tn+1) = T(tn) + h[θ

T(tn+1) + (1 − θ)

d
dt

d
dt

T(tn)]

(6.14)

where the variables are discretized in time domain, tn = nh, n = 0, 1,··· and h is the time

step. θ is a number between 0 and 1. When θ = 0, (6.14) reduces to the Euler Method.

2, (6.14) is the trapezoidal method. Notice that for θ = 1

When θ = 1
2, θ-method is a second-
order method, which means the local error at each step is of order O(h3), giving a global
error of order O(h2). For all other θ values, θ-method is a first-order method.

From (6.1),

T(tn) = L−1[−RT(tn) − M

d
dt

d
dt

i(tn)],

T(tn+1) = L−1[−RT(tn+1) − M

d
dt

d
dt

i(tn+1)].

(6.15)

(6.16)

71

Substituting (6.15) and (6.16) into (6.14), and applying θ = 1

2, we have

T(tn+1) = T(tn) − 1
2

hL−1[RT(tn+1) + M

d
dt

i(tn+1) + RT(tn) + M

d
dt

i(tn)]

(6.17)

This gives an explicit computation scheme of T(t) given the knowledge of i(t) and T(0) =

0,

(I +

1
2

hL−1R)T(tn+1) = T(tn) − 1

2

hL−1[M

d
dt

i(tn+1) + RT(tn) + M

d
dt

i(tn)]

(6.18)

In practice, (I + 1

2 hL−1R) is a large matrix and computing its inverse is time consuming.
LU decomposition and back substitution technique is used to update T(tn+1) at each step.

To validate the measurement of transfer function from transient waveforms, let us assume
an example problem. The test sample is a 3” × 3” × 1/4” metal plate with a electrical
conductivity of 5 × 106 S/m. The sample is discretized into 20 × 20 × 6 elements in x, y, z

axes respectively. A pair of transmitter/receiver coil is placed on top of the sample with

0.825 mm lift-off. The coil thickness is 1 mm, the inner and outer diameter are 4 mm and 5

mm, respectively (figure 6.2).

Figure 6.2: An example problem: the mesh of specimen and the transmitter/receiver coil
setup.

72

The true transfer function computed with Z = − 1

τ MT (τ R − L)−1M is:

−8.2605 × 10−7 −3.3128 × 10−7

−3.3128 × 10−7 −8.2605 × 10−7



Z =

(6.19)

On the other hand, the elements of ”measured” transfer function is calculated as the

ratio between pick-up coil voltage and drive current. Figure 6.3 illustrates these transient

waveforms. It is clear that this ratio converges to a constant value at large t. In practice,

the asymptotic values of these curves are estimated as the transfer function. These values,
as read in figure 6.3(c), are −8.2605 × 10−7 and −3.3128 × 10−7, which are identical to the

elements of true transfer function Z in (6.19).

6.3.2 Noise Free Reconstruction

This part of the results consider the ideal condition where the transfer function is obtained

with zero noise. To compare with the results obtained using monotonicity of time constants,

the test sample as well as the region of interest (ROI) remain the same as in section 5.2. The
length, width and depth of the ROI are 6 mm, 0.33 mm and 3 mm, respectively. A 2 × 6

coil array is arranged above the ROI (see figure 6.4) and hence the dimension of impedance
matrix is 12 × 12. The coils are placed at 0.1 mm above the sample surface, the inner and

outer diameter of the coils are 0.8 mm and 1 mm, respectively. The number of turns is 100.

Four defect profiles are considered: (a) L-shape top surface opening (b) square shape buried

defect (c) T-shape bottom surface opening (d) two separate square defects. An illustration

of all defects can be found in figure 6.5.

The reconstruction results using noise free transfer function are presented in the figure 6.6,

73

(a)

(b)

(c)

Figure 6.3: Transient waveform and transfer function. (a) input current in coil 1 (b) pick-up
voltages (c) impedance, the asymptotic values are the transfer function

74

Figure 6.4: Sensor array configurations.

defect a

defect b

defect c

defect d

Figure 6.5: Example defect profiles.

together with a comparison with the real defect profiles. The reconstruction is successful for

all configurations, including multiple defects. When compared with time constants, transfer

function is a localized property, and hence is not insensitive to sensor liftoff and tilting,

but it has the following advantages: (i) it is easy to measure, (ii) it does not require extra

conductive patches to break the isometry.

75

defect a

defect b

defect c

defect d

Figure 6.6: Reconstruction via monotonicity of transfer function using noise free data. Top
figure: real defect profiles (in red color). Bottom figure: reconstructed defect profiles (in
light blue color)

76

6.3.3 Noisy Reconstruction

The specimen and sensor system geometry are identical to the one described in section 6.3.2.

An additive noise is assumed on the transfer function

˜ZV = ZV + N,

(6.20)

where ZV is the true transfer function of unknown defect V , of which ˜ZV is the noisy version,

and N is the noise matrix. Equivalently, in terms of matrix elements,

( ˜ZV )i,j = (ZV )i,j + ξi,j · ∆Zmax

(6.21)

where ∆Zmax = maxi,j |(ZV )i,j − (ZBG)i,j| and BG denotes the background (defect-free)

configuration. ξi,j is a Gaussian random variable with zero mean and standard deviation σ,
− x2
2σ2 ,∀i, j. The estimate of noise

1√
2πσ2

i.e. the probability density function f (ξi,j = x) =
level δ, as introduced in section 6.2, is selected to be δ = (cid:107)N(cid:107)

e

Different values of the prescribed decay constant τ were selected. The maximum natural
time constant of the background specimen is τ1 = 2.9405× 10−5, and τ is selected to be 5τ1

and 2τ1, respectively.

77

Case 1: τ = 5τ1, the L-2 norm of matrix Zdef ect is 0.0027, the maximum eigenvalue of

matrix Zdef ect is −1.5572 × 10−4 and the minimum eigenvalue is -0.0027.

For a prescribed noise level σ, as introduced in (6.21), we repeated the numerical ex-

periment for N = 1000 times by generating the time constants using additive noise model

(6.21). The probability of a voxel i being identified as part of the anomaly is defined as

pi =

Ni
N , where Ni is the number of times the voxel i is identified as part of the anomaly.

The imaging results at different noise levels are presented in figure 6.7 to figure 6.10. The

voxel of a darker color is more probable a part of the anomaly.

(a) true anomaly profile

(b) σ = 0.0001, SNR = 46

(c) σ = 0.001, SNR = 35

(d) σ = 0.01, SNR = 25

Figure 6.7: Statistical reconstruction results for defect #1, τ = 5τ1

78

(a) true anomaly profile

(b) σ = 0.0001, SNR = 56

(c) σ = 0.001, SNR = 46

(d) σ = 0.01, SNR = 36

Figure 6.8: Statistical reconstruction results for defect #2, τ = 5τ1

79

(a) true anomaly profile

(b) σ = 0.0001, SNR = 65

(c) σ = 0.001, SNR = 56

(d) σ = 0.01, SNR = 46

Figure 6.9: Statistical reconstruction results for defect #3, τ = 5τ1

80

(a) true anomaly profile

(b) σ = 0.0001, SNR = 59

(c) σ = 0.001, SNR = 49

(d) σ = 0.01, SNR = 38

Figure 6.10: Statistical reconstruction results for defect #4, τ = 5τ1

81

Case 2: τ = 2τ1, the L-2 norm of matrix Zdef ect is 0.0214, the maximum eigenvalue of
matrix Zdef ect is −9.7582×10−4 and the minimum eigenvalue is -0.0214. Similarly, different

values of σ are selected and the reconstruction results averaged over 1000 implementations

are presented in figure 6.11 to figure 6.14.

(a) true anomaly profile

(b) σ = 0.0001, SNR = 47

(c) σ = 0.001, SNR = 37

(b) σ = 0.01, SNR = 27

Figure 6.11: Statistical reconstruction results for defect #1, τ = 2τ1

A few conclusions can be drawn from these two examples:

• All defect profiles, including shape, location, and numbers, can still be properly recon-

structed under noisy conditions.

• When the noise level increases, the performance of the reconstruction deteriorates.

When the SNR is smaller than 40, the reconstruction error increases at larger depths..

This is reasonable because these pixels are far away from the excitation coil array and

82

(a) true anomaly profile

(b) σ = 0.0001, SNR = 56

(c) σ = 0.001, SNR = 46

(b) σ = 0.01, SNR = 36

Figure 6.12: Statistical reconstruction results for defect #2, τ = 2τ1

hence, they are less coupled and produce less disturbance to the measurement.

• The reconstruction using smaller τ value (in this case, τ = 2τ1) is better than larger

τ value (τ = 5τ1) under same noise levels. However, there is a lower limit of τ = τ1 to

satisfy the monotonicity (6.9). Besides, in practice, it is more difficult to measure the

transfer function with smaller τ values since measurements have to be made after the

natural response vanishes and before the forced response falls below the noise floor. It

is important that there is a significant separation between τ and τ1.

83

(a) true anomaly profile

(b) σ = 0.0001, SNR = 65

(c) σ = 0.001, SNR = 55

(b) σ = 0.01, SNR = 45

Figure 6.13: Statistical reconstruction results for defect #3, τ = 2τ1

6.4 Concluding Remarks

This chapter has been focused on the concept, validation and application examples of the

monotonicity of transfer function, with support of numerical evidences. Monotonicity of

time constants and monotonicity of transfer function are two complementary properties.

While time constants are related to the source-free response of the eddy current problem, the

transfer function is obtained from the forced response only. Compared to time constants, the

transfer function is a localized property and it is not insensitive to probe lift-off and tilting.

However, for the same exact reason, monotonicity of transfer function does not suffer from

the problem of isometry when applied to defect imaging. Moreover, the transfer function

is directly measurable from eddy current responses and hence, is much more practical in

84

(a) true anomaly profile

(b) σ = 0.0001, SNR = 58

(c) σ = 0.001, SNR = 48

(b) σ = 0.01, SNR = 38

Figure 6.14: Statistical reconstruction results for defect #4, τ = 2τ1

real-life inspections.

The estimation of transfer function requires exponentially decaying drive current, and the

decay constants should be significantly greater than the time constants of natural modes.

Next chapter will present the efforts in building an experiment system which satisfies such

conditions. Experimental results of defect imaging via monotonicity-based approach will

also be provided.

85

Chapter 7

Experiment Setup and Results

The experiment setup consists of three major parts: excitation system, coil array, and mea-

surement system, which are all driven and controlled by a computer/control unit (Figure

7.1).

Figure 7.1: Schematic of experiment setup

86

7.1 Experiment Setup

A exponential decaying current is applied to properly estimate the transfer function. While

such a signal can be generated easily with an arbitrary waveform generator, a waveform

generator does not always drive current of identical waveform. The drive current depends

on the load connected to it. Specifically, a typical waveform generator has a internal output

impedance of 50 Ω and the excitation coil features a very small resistance and significant

inductance. In practice, the actual drive current can be monitored by a small resistor which

is connected serially to the excitation coil. It is heavily distorted when the coil is directly

connected to the waveform generator.

A voltage controlled current source (VCCS) drives current that strictly follows the input

voltage waveform. The schematic of the specific implementation of VCCS in this study is

shown in figure 7.2. When an operational amplifier is operating in negative feedback mode,

the two input terminals (node #9 and node #10) are virtually short. The output current is

solely dependent on the input voltage and the feedback resistor R8, I = Vin/R8, regardless

of the load connected between node J5 and J8.

This circuit is built using the audio power amplifier LM3886. LM3886 is a high perfor-

mance amplifier capable of delivering 68W of continuous average power to a 4Ω load. It has

a large gain-bandwidth product of 8 MHz, which is also highly desired in this experiment.

A considerable amount of drive current of the order of a few Amps is generated, so resistor

R8 is chosen to be a power resistor and proper heat dissipation is ensured. A photograph of

the excitation system is presented in figure 7.3.

The frequency response of this circuit under different loads are measured. Specifically, a

5.6 Ω resistor and a 100-turn coil (with its DC resistance measured as 6.8 Ω) are used. The

87

Figure 7.2: Design of the voltage controlled current source.

frequency response, in terms of magnitude and phase, are given in figure 7.4. The magnitude

is measured as the ratio between the voltage VR8 across the feedback resistor R8 and the

input signal Vin. It is clear that the magnitude is almost unity and the phase shift is very

close to zero up to 100 kHz. Below this frequency, the drive current follows the input signal

very well.

An exponentially decaying signal was programmed via an arbitrary waveform generator

(Agilent 33500B) and fed to this circuit. The decay constant is 50 µs. The drive current,

represented by the voltage across the feedback voltage, was plotted in figure 7.5. Please notice

that the current waveform is measured when connecting to a 100-turn coil. An exponentially

decaying current was achieved.

A sensor comprising two coaxial pancake coils (figure 7.6 shows the two coils separately,

the smaller coil is inserted into the larger one when used) is employed as the probe in this

experiment. The exterior/larger coil has an inner diameter of 7 mm, outer diameter of 11.5

88

Figure 7.3: The voltage controlled current source is connected to a pair of transmitter-receiver
coils for testing performance.

mm, height of 5.4 mm and its number of turns is 60. The interior/smaller coil has an inner

diameter of 3 mm, outer diameter of 6 mm, height of 2.8 mm and its number of turns is 20.

Copper wire of AWG=30 is used to prepare these coils.

The specimen under test is a plate of aluminum alloy 6061 (length = 12 inch, width
= 6 inch, depth = 1/4 inch, electrical conductivity = 2.5 × 107 S/m). Three defects were

artificially fabricated with dimensions as shown in Figure 7.7. Defect 1 is of dimension 10
mm × 10 mm and serves as the test element in the inversion procedure. Defects 2 and 3

serve as the unknown defects to be reconstructed. The depth of all these defects is 3 mm.

The time domain waveforms including both the drive current and the induced coil voltage,

89

(a)

(b)

Figure 7.4: Frequency response of the current source circuit.

are collected via the oscilloscope (model InfiniiVision DSO-X 4052A). The sampling rate is
2 × 106 samples per second.

90

Figure 7.5: The excitation current agrees very well with the input signal.

Figure 7.6: The two coils consisting the array. The smaller coil is inserted into the larger
one when used.

7.2

Imaging Scheme and Results

When conducting an inspection, the sensor array was mounted on an automated 3D scanner.

A 2D area is scanned while keeping a constant 0.5 mm lift-off from the specimen. The scan

91

Figure 7.7: Three artificial defects are fabricated on an aluminum test sample.

92

area was subdivided into pixels of dimensions 10 mm × 10 mm, which is identical to the

dimensions of a test element. The transfer function was measured when the barycenter of

the sensor was on the center of each pixel.

In practice, due to the localized nature of eddy current phenomenon, for a given test

element Tk, only the measurements from positions j in the proximity of this test element is

meaningful. Adapting the idea from [39], the unknown defect V can be estimated as

(cid:91){Tk|Zcoil,j

Tk

˜V =

≤ Zcoil,j

V

,∀j ∈ Ik}.

(7.1)

where the superscript coil, j means this is the transfer function measured at coil position j,

and Ik is the set of proper measurement positions for test element Tk. In this experiment,

the 9 positions centering the test element were considered the set of proper positions, as

numbered in figure 7.8 (the test element is at the center of these positions, which is position

5). The transfer function for this single pixel (test element) was experimentally measured

and processed in advance. Thanks to the translation invariance of the problem (assuming

the region of interest is far away from any edge or other defects), the measurement of a single

test element also comprises the transfer function of any other arbitrary pixel in the grid.

In practice, the measurements are corrupted by noise. To acquire a constant transfer

function for each position, an exponential curve was fitted to both the drive current and

pick-up coil voltage. Since we prescribed the decay constant in the drive current (τ = 50µs),

curve fitting problem reduced to a simple linear regression. With this experimental setup,

the noise level in the voltage measurement was of the order of 10 mV. Figure 7.9 shows an

example of experimental waveforms and the corresponding fitted exponential curves.

For each ”unknown” defect, a region consisting of 25 pixels was scanned and the transfer

93

Figure 7.8: A test element and the coil positions where the measurements are relevant.

function was measured at each pixel. These measurements are used to determined whether

any of the 9 elements in the center is part of the defect, as formulated in eq. (7.1). The

reconstruction results are illustrated in figure 7.10. With this relatively simple surface-

breaking problem, the imaging results are error-free even under experimental noise.

94

Figure 7.9: The measurements of experimental current and voltage waveforms and their
fitting to exponential curves.

95

Figure 7.10: Experimental results using the imaging approach based on monotonicity of
transfer function.

96

Chapter 8

Conclusions and Future Work

8.1 Conclusions

This dissertation developed a fast non-iterative method for imaging defects in electrically

conductive materials using time-domain eddy current test data and the monotonicity prop-

erties.Two monotonicity properties and related imaging algorithms were formulated, and

results of implementation are also presented to validate the approach.

This study started by formulating the governing equation of the eddy current problem in

terms of electric vector potentials. The result is a parabolic partial differential equation. The

solution to this governing equation consists of two parts: the source-free response and the

forced response. The source-free response, i.e. the solution to the homogeneous problem, can

be expressed as sum of natural modes: ufree(r, t) =(cid:80)

k ckuk(r)e−t/τk . The time constants

τk are real, positive, bounded and approach to zero when arranged in decreasing order.

It was proved that if an anomaly in a conductive domain expands, all the time constants

associated with this anomaly decrease. This property is called the monotonicity of time

constants.

It leads to a fast, non-iterative imaging method, which is referred to as the

Monotonicity Principle Method (MPM): the region of interest (ROI) is sampled with test

elements, if a test element is part of an unknown defect, all the time constants associated

with the test element are greater than or equal to the time constants associated with the

97

unknown defect. If the test element is not part of the unknown defect, their time constants are

not monotonic. In practice, the time constants of test elements are numerically computed

or experimentally measured in advance. Hence, the imaging procedure requires only the

comparison of time constants, which takes negligible computation cost. Moreover, different

test elements could be processed in parallel to further reduce processing time.

From the study of natural modes and time constants, it was noticed that natural modes

are global. Time constants alters significantly when its related natural mode interacts

strongly with a defect, i.e. the natural mode will be changed notably if this defect is intro-

duced. Since natural modes are eigenvectors of the homogeneous problem, they are global

and independent of the excitation system. Some natural modes circulates mainly in the

center part of the problem domain, in other words, their time constants may help to find

deeply buried defects at the center of the sample. However, for the same reason that natural

modes are independent of excitation, time constants can not distinguish two configurations

if they are isometric.

In this study, we solved this problem by using an extra conductor

patch to break the isometry.

Imaging of an anomaly requires the knowledge of time constants estimated from the tran-

sient signals of a pulsed eddy current testing. This report reviewed numerical techniques for

multiexponential analysis and nonexponential analysis and, proposed a modified algorithm

based on the Laplace-Pad´e approximant approach. This algorithm takes into account that

multiple measurements share the same set of time constants. Numerical examples showed it

is robust with respect to noise and preserve the monotonicity property in the estimated time

constants. We exploited this algorithm to estimate time constants from realistic transient

signals, the results are promising in the problem of classifying anomalies of different sizes.

MPM of time constants is advantageous over other imaging methods in PECT for two

98

major reasons: (a) this is a non-iterative method requiring extremely low computation cost

and hence is suitable for real-time imaging and, (b) the time constants are intrinsic to test

samples and independent of the inspection system, and hence not sensitive to probe lift-off,

tilt and excitation waveforms. The inspection system, however, may affect the observability

of the time constants.

The second monotonicity property is derived from the forced response of the eddy cur-

rent problem. Specifically, when the system is excited with an exponentially decaying cur-
rent i(t) = i0e−t/τ , the input-output relationship can be expressed as a constant matrix:
uforced(t) = Zi(t). Notice that in practice, this requires τ >> τk and the transfer function

should be measured after the natural response vanishes and before the forced response falls

below the noise floor.

The monotonicity of transfer function states that if an anomaly in a conductive domain

expands, the transfer function of the system increases. This dissertation also developed an

non-iterative imaging method exploiting the monotonicity of transfer function and validated

it with numerical experiments. When compared to MPM of time constants, transfer function

does not have the advantages of being insensitive to probe lift-off and tilting, but it is much

more practical and easier to implement.

An experimental system was built to measure transfer functions. This system consisted

of a voltage-controlled current source, a coil array, a 3D scanner, and a digital oscilloscope.

Data was collected using an automated data collection software developed in C#. Experi-

mental results validated the approach based on monotonicity of transfer function. Error-free

reconstruction examples using the monotonicity-based imaging algorithm show significant

promise.

99

8.2 Future Work

Monotonicity of time constants and the imaging method based on it are very appealing, as

it is not sensitive to probe perturbation and excitation waveforms. However, the extraction

of time constants from transient waveform remains a very challenging problem. Researchers

from different areas have tried to address this issue. Different methods, e.g.

statistical

methods, direct inversion methods and optimization methods, have been applied to it but

none yields an accurate enough result yet. A physics-based a priori condition may be

required to regularize this problem.

Transfer function is a local property. Good measurement of transfer function requires a

compact and careful design of coil array. Coils in the array should be nicely coupled with

each other in order to achieve high SNR. The example in Chapter 7 (two coaxial coil pair) is

just a preliminary attempt, a systematic approach to design the coil array is yet to develop.

100

APPENDIX

101

Appendix

Min-Max Characterization of Time

Constants

Following [50], let Uk and r (x) be defined as

Uk = span{u1, ..., uk}

r (x) =

xT Lx
xT Rx

.

r (x) is the well know Rayleigh quotient.

Here we first prove that

and, finally, that

τk = min
x∈Uk

(cid:32)

xT Lx
xT Rx

τk = max

dim(U )=k

min
x∈U

xT Lx
xT Rx

(cid:33)

.

(.1)

(.2)

First, we notice that (see [71]) the stationary points of the Rayleigh quotient provide the gen-

eralized eigenvectors and the corresponding values of r (x) give the generalized eigenvalues.

Indeed,

∇r (x) =

2Lx − 2r (x) Rx

xT Rx

and, therefore, ∇r (x) = 0 ⇔ Lx − r (x) Rx = 0. (.1) follows easily because the stationary

102

points of r (x) in Uk are u1, ..., uk and the decreasing ordering of the τk’s.

To prove (.2) we notice that for a given k dimensional subspace U , we have a non vanishing

v ∈U ∩ span{uk, ..., uN}. Then

min
x∈U

and

r (x) ≤ r (v) ≤

max

x∈span{uk,...,uN} r (x) = τk
(cid:18)

(cid:19)

max

dim(U )=k

min
x∈U

r (x)

≤ τk.

(.3)

On the other hand, from (.1) it follows that

(cid:32)

(cid:33)

≥ min
x∈Uk

xT Lx
xT Rx

= τk

xT Lx
xT Rx

max

dim(U )=k

min
x∈U

which, together with (.3) proves (.2).

103

BIBLIOGRAPHY

104

BIBLIOGRAPHY

[1] “Introduction to Nondestructive Testing,” https://www.asnt.org/MinorSiteSections/

AboutASNT/Intro-to-NDT/, [Online; accessed 24-March-2017].

[2] S. S. Udpa and P. O. Moore, Nondestructive Testing Handbook, volume 5–, 3rd ed., ser.

Nondestructive Testing Handbook. ASNT, 2004, vol. 5.

[3] D. McCann and M. Forde, “Review of ndt methods in the assessment of concrete and

masonry structures,” Ndt & E International, vol. 34, no. 2, pp. 71–84, 2001.

[4] M. Jolly, A. Prabhakar, B. Sturzu, K. Hollstein, R. Singh, S. Thomas, P. Foote, and
A. Shaw, “Review of non-destructive testing (ndt) techniques and their applicability to
thick walled composites,” Procedia CIRP, vol. 38, pp. 129–136, 2015.

[5] G. Yang, A. Tamburrino, L. Udpa, S. S. Udpa, Z. Zeng, Y. Deng, and P. Que, “Pulsed
eddy-current based giant magnetoresistive system for the inspection of aircraft struc-
tures,” IEEE transactions on magnetics, vol. 46, no. 3, pp. 910–917, 2010.

[6] J. Garc´ıa-Mart´ın, J. G´omez-Gil, and E. V´azquez-S´anchez, “Non-destructive techniques

based on eddy current testing,” Sensors, vol. 11, no. 3, pp. 2525–2565, 2011.

[7] F. Rachidi, M. Rubinstein, J. Montanya, J.-L. Bermudez, R. R. Sola, G. Sola, and
N. Korovkin, “A review of current issues in lightning protection of new-generation wind-
turbine blades,” IEEE Transactions on Industrial Electronics, vol. 55, no. 6, pp. 2489–
2496, 2008.

[8] G. Chen, A. Yamaguchi, and K. Miya, “A novel signal processing technique for eddy-
current testing of steam generator tubes,” IEEE Transactions on magnetics, vol. 34,
no. 3, pp. 642–648, 1998.

[9] Z. Chen, N. Yusa, and K. Miya, “Enhancements of eddy current testing techniques
for quantitative nondestructive testing of key structural components of nuclear power
plants,” Nuclear Engineering and Design, vol. 238, no. 7, pp. 1651–1656, 2008.

[10] G. L. Fitzpatrick, D. K. Thome, R. L. Skaugset, E. Y. Shih, and W. C. Shih, “Magneto-
optic/eddy current imaging of aging aircraft: A new ndi technique,” Materials Evalua-
tion;(United States), vol. 51, no. 12, 1993.

[11] R. Lange and G. Mook, “Structural analysis of cfrp using eddy current methods,” NDT

& E International, vol. 27, no. 5, pp. 241–248, 1994.

105

[12] L. Cheng and G. Y. Tian, “Surface crack detection for carbon fiber reinforced plas-
tic (cfrp) materials using pulsed eddy current thermography,” IEEE Sensors Journal,
vol. 11, no. 12, pp. 3261–3268, 2011.

[13] J.-S. Bae and S.-Y. Kim, “Hot wire inspection using eddy current,” in Instrumentation
and Measurement Technology Conference, 2001. IMTC 2001. Proceedings of the 18th
IEEE, vol. 2.

IEEE, 2001, pp. 962–965.

[14] T. Chen, G. Y. Tian, A. Sophian, and P. W. Que, “Feature extraction and selection for
defect classification of pulsed eddy current ndt,” Ndt & E International, vol. 41, no. 6,
pp. 467–476, 2008.

[15] A. Sophian, G. Y. Tian, D. Taylor, and J. Rudlin, “A feature extraction technique based
on principal component analysis for pulsed eddy current ndt,” NDT & e International,
vol. 36, no. 1, pp. 37–41, 2003.

[16] G. Y. Tian and A. Sophian, “Reduction of lift-off effects for pulsed eddy current ndt,”

NDT & E International, vol. 38, no. 4, pp. 319–324, 2005.

[17] V. Isakov, Inverse problems for partial differential equations. Springer Science & Busi-

ness Media, 2006, vol. 127.

[18] L. Udpa and W. Lord, “A discussion of the inverse problem in electromagnetic ndt,”
Springer, 1986, pp.

in Review of progress in quantitative nondestructive evaluation.
375–382.

[19] A. Tikhonov and V. Y. Arsenin, Methods for solving ill-posed problems.

John Wiley

and Sons, Inc, 1977.

[20] E. Coccorese, R. Martone, and F. C. Morabito, “A neural network approach for the
solution of electric and magnetic inverse problems,” IEEE Transactions on Magnetics,
vol. 30, no. 5, pp. 2829–2839, 1994.

[21] A. Ayad, F. Benhamida, A. Bendaoud, Y. Bihan, and M. Bensetti, “Solution of inverse
problems in electromagnetic ndt using neural networks,” Przeglad Elektrotechniczny,
vol. 2011, no. 9a, pp. 330–333, 2011.

[22] S. R. H. Hoole, “Artificial neural networks in the solution of inverse electromagnetic
field problems,” IEEE transactions on Magnetics, vol. 29, no. 2, pp. 1931–1934, 1993.

[23] F. Hettlich, “The landweber iteration applied to inverse conductive scattering prob-

lems,” Inverse Problems, vol. 14, no. 4, p. 931, 1998.

106

[24] P. Pinheiro, W. Loh, and F. Dickin, “Three-dimensional reconstruction algorithm for
electrical resistance tomography,” IEE Proceedings-Science, Measurement and Technol-
ogy, vol. 145, no. 3, pp. 85–93, 1998.

[25] S. J. Norton and J. R. Bowler, “Theory of eddy current inversion,” Journal of Applied

Physics, vol. 73, no. 2, pp. 501–512, 1993.

[26] W. Lionheart, M. Soleimani, and A. Peyton, “Sensitivity analysis in 3d magnetic in-
duction tomography (mit),” in Proc. 3rd World Congr. Industrial Process Tomography,
2003, pp. 239–244.

[27] S. Caorsi, A. Massa, and M. Pastorino, “A crack identification microwave procedure
based on a genetic algorithm for nondestructive testing,” IEEE Transactions on Anten-
nas and Propagation, vol. 49, no. 12, pp. 1812–1820, 2001.

[28] M. Benedetti, G. Franceschini, R. Azaro, and A. Massa, “A numerical assessment of
the reconstruction effectiveness of the integrated ga-based multicrack strategy,” IEEE
Antennas and Wireless Propagation Letters, vol. 6, pp. 271–274, 2007.

[29] I. H. Osman and G. Laporte, “Metaheuristics: A bibliography,” 1996.

[30] N. Yusa, Z. Chen, K. Miya, T. Uchimoto, and T. Takagi, “Large-scale parallel compu-
tation for the reconstruction of natural stress corrosion cracks from eddy current testing
signals,” NDT & E International, vol. 36, no. 7, pp. 449–459, 2003.

[31] D. Colton and A. Kirsch, “A simple method for solving inverse scattering problems in

the resonance region,” Inverse problems, vol. 12, no. 4, p. 383, 1996.

[32] A. Kirsch, “Characterization of the shape of a scattering obstacle using the spectral

data of the far field operator,” Inverse problems, vol. 14, no. 6, p. 1489, 1998.

[33] M. Hanke and M. Br¨uhl, “Recent progress in electrical impedance tomography,” Inverse

Problems, vol. 19, no. 6, p. S65, 2003.

[34] A. J. Devaney, “Super-resolution processing of multi-static data using time reversal and

music,” 2000.

[35] H. Ammari, E. Iakovleva, and D. Lesselier, “A music algorithm for locating small in-
clusions buried in a half-space from the scattering amplitude at a fixed frequency,”
Multiscale Modeling & Simulation, vol. 3, no. 3, pp. 597–628, 2005.

[36] A. Tamburrino and G. Rubinacci, “A new non-iterative inversion method for electrical

resistance tomography,” Inverse Problems, vol. 18, no. 6, p. 1809, 2002.

107

[37] A. Tamburrino, S. Ventre, and G. Rubinacci, “Recent developments of a monotonicity
imaging method for magnetic induction tomography in the small skin-depth regime,”
Inverse Problems, vol. 26, no. 7, p. 074016, 2010.

[38] A. Tamburrino and G. Rubinacci, “Fast methods for quantitative eddy-current tomog-
raphy of conductive materials,” IEEE transactions on magnetics, vol. 42, no. 8, pp.
2017–2028, 2006.

[39] A. Tamburrino, F. Calvano, S. Ventre, and G. Rubinacci, “Non-iterative imaging
method for experimental data inversion in eddy current tomography,” NDT & E Inter-
national, vol. 47, pp. 26–34, 2012.

[40] A. Tamburrino, L. Barbato, D. Colton, and P. Monk, “Imaging of dielectric objects via
monotonicity of the transmission eigenvalues,” Proc. of The 12th International Confer-
ence on Mathematical and Numerical Aspects of Wave Propagation, 2015.

[41] B. Harrach and M. Ullrich, “Monotonicity-based shape reconstruction in electrical
impedance tomography,” SIAM Journal on Mathematical Analysis, vol. 45, no. 6, pp.
3382–3403, 2013.

[42] H. Garde and S. Staboulis, “Convergence and regularization for monotonicity-based
shape reconstruction in electrical impedance tomography,” Numerische Mathematik,
pp. 1–31, 2016.

[43] D. Flores-Tapia and S. Pistorius, “Electrical impedance tomography reconstruction us-
ing a monotonicity approach based on a priori knowledge,” in Engineering in Medicine
and Biology Society (EMBC), 2010 Annual International Conference of the IEEE.
IEEE, 2010, pp. 4996–4999.

[44] D. Flores-Tapia, M. O’Halloran, and S. Pistorius, “A bimodal reconstruction method for
breast cancer imaging,” Progress In Electromagnetics Research, vol. 118, pp. 461–486,
2011.

[45] A. Tamburrino, Z. Su, N. Lei, L. Udpa, and S. Udpa, “The monotonicity imaging
method for pect,” Electromagnetic Nondestructive Evaluation (XVIII), vol. 40, p. 159,
2015.

[46] Z. Su, A. Tamburrino, S. Ventre, L. Udpa, and S. Udpa, “Time domain monotonicity
based inversion method for eddy current tomography,” in Applied Computational Elec-
tromagnetics (ACES), 2015 31st International Review of Progress in.
IEEE, 2015, pp.
1–2.

[47] A. Tamburrino, Z. Su, S. Ventre, L. Udpa, and S. S. Udpa, “Monotonicity based imang
method in time domain eddy current testing,” Electromagnetic Nondestructive Evalua-
tion (XIX), vol. 41, p. 1, 2016.

108

[48] R. Albanese and G. Rubinacci, “Integral formulation for 3d eddy-current computation
using edge elements,” IEE Proceedings A-Physical Science, Measurement and Instru-
mentation, Management and Education-Reviews, vol. 135, no. 7, pp. 457–462, 1988.

[49] ——, “Finite element methods for the solution of 3d eddy current problems,” Advances

in Imaging and Electron Physics, vol. 102, pp. 1–86, 1997.

[50] A. Tamburrino, “Monotonicity for parabolic pdes: Eddy current testing,” personal

communication, pp. 1–12, 2014.

[51] R. Albanese and G. Rubinacci, “Treatment of multiply connected regions in two-
component electric vector potentials formulations,” IEEE Transactions on Magnetics,
vol. 26, no. 2, pp. 650–653, 1990.

[52] B. Harrach and M. Ullrich, “Resolution guarantees in electrical impedance tomography,”

IEEE transactions on medical imaging, vol. 34, no. 7, pp. 1513–1521, 2015.

[53] A. A. Istratov and O. F. Vyvenko, “Exponential analysis in physical phenomena,”

Review of Scientific Instruments, vol. 70, no. 2, pp. 1233–1257, 1999.

[54] R. Prony, “Essai experimental–,-,” J. de lEcole Polytechnique, 1795.

[55] R. G. Cornell, “A method for fitting linear combinations of exponentials,” Biometrics,

vol. 18, no. 1, pp. 104–113, 1962.

[56] J. Loeb and G. Cahen, “Extraction `a partir des enregistrements de mesures des

param`etres dynamiques d’un syst`eme,” Automatisme, vol. 8, pp. 479–486, 1963.

[57] ——, “More about process identification,” IEEE Transactions on Automatic Control,

vol. 10, no. 3, pp. 359–361, 1965.

[58] Z. Bay, “Calculation of decay times from coincidence experiments,” Physical Review,

vol. 77, no. 3, p. 419, 1950.

[59] E. Yeramian and P. Claverie, “Analysis of multiexponential functions without a hy-
pothesis as to the number of components,” Nature, vol. 326, no. 6109, pp. 169–174,
1987.

[60] A. Tikhonov and V. Y. Arsenin, Methods for solving ill-posed problems.

John Wiley

and Sons, Inc, 1977.

[61] D. G. Gardner, J. C. Gardner, G. Laush, and W. W. Meinke, “Method for the analysis
of multicomponent exponential decay curves,” The journal of chemical physics, vol. 31,
no. 4, pp. 978–986, 1959.

109

[62] D. G. Gardner, “Resolution of multi-component exponential decay curves using fourier
transforms,” Annals of the New York Academy of Sciences, vol. 108, no. 1, pp. 195–203,
1963.

[63] C. E. Shannon and W. Weaver, The mathematical theory of communication. University

of Illinois press, 2002.

[64] S. Meshkov and D. Berkov, “Fast maximum entropy algorithm for ill-posed problems,”

International Journal of Modern Physics C, vol. 5, no. 06, pp. 987–995, 1994.

[65] C. Lanczos, Applied analysis. Courier Corporation, 1988.

[66] W. Wiscombe and J. Evans, “Exponential-sum fitting of radiative transmission func-

tions,” Journal of Computational Physics, vol. 24, no. 4, pp. 416–444, 1977.

[67] K. Holmstr¨om and J. Petersson, “A review of the parameter estimation problem of
fitting positive exponential sums to empirical data,” Applied Mathematics and Compu-
tation, vol. 126, no. 1, pp. 31–61, 2002.

[68] E. H. Hellen, “Pad´e–laplace analysis of signal averaged voltage decays obtained from a

simple circuit,” American journal of physics, vol. 73, no. 9, pp. 871–875, 2005.

[69] J. Kiefer, “Sequential minimax search for a maximum,” Proceedings of the American

mathematical society, vol. 4, no. 3, pp. 502–506, 1953.

[70] C. R. J. Roger A. Horn, “Definitions and properties,” in Matrix Analysis (Second Edi-

tion). Cambridge University Press, 2012, ch. 7.1, p. 431.

[71] T. Melzer, “Svd and its application to generalized eigenvalue problems,” Vienna Uni-

versity of Technology, 2004.

110