by

L
. ,5... 1.

it..0u.§f§11§a&£$§ww . .1 . 1 . _ .. gwmwnzhsé .fiafifiwmgg imfimfifla .,
..1. . 1 ..1 :91... ..1 1... . ... .2“. ..1 , . , J . .. . ., ..!#Q

..I 1 .11..

1 ‘ I‘

 

WESls
l
J4
215

n .5‘
593v 3

This is to certify that the
thesis entitled

Independent Component Analysis for Enhanced Feature
Extraction in NDE Applications

presented by
Byung Hyuk Shin

has been accepted towards fulfillment
of the requirements for the

MS. degree in Electrical Engineering

 

 

(La/MC «Z Lat/wt,

 

Major Professor’sgignature

L

Date

MSU is an Affirmative Action/Equal Opportunity Institution

 

LIBRARY
Michigan State
University

 

 

 

| PLACE IN RETURN Box to remove this checkout from your record.
To AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/01 cJClRC/DateDuopes-p. 1 5

 

Independent Component Analysis for Enhanced Feature Extraction
from Non-destructive Evaluation

By

Byung Hyuk Shin

A THESIS

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

MASTER OF SCIENCE

Department of Electrical and Computer Engineering

2004

ABSTRACT

Independent Component Analysis for Enhanced Feature Extraction
in NDE Applications

By

Byung Hyuk Shin

In this thesis. independent component analysis (ICA) is proposed for enhancing flaw
information in eddy cun‘ent nondestructive evaluation (NDE). ICA provides a method for
representing data as weighted combination of independent components and higher—order
statistics of the data is used to minimize the dependence between the components of the
system output. Mum-frequency eddy current testing is a widely used NDE method in
situations where the defect signal is con'upted by noise and contributions from extemal
supports that makes the analysis challenging. ICA, along with an affine transfomiation as a
preprocessing stage. is shown to extract the defect signal from a combination of defect.

support and noise signals. while improving the SNR.

ACKNOWLEDGEMENTS

I would like to express my deep gratitude to major professor Dr. Lalita Udpa for
her kind guidance, support, and encouragement during my graduate study. I also thank
her for taking time to revise my thesis. I would sincerely like to thank Dr. Satish Udpa
and Dr. Pradeep Ramuhalli for their technical support and invaluable advice through the
project. I wish to express my special thanks to the members of NDEL at Michigan State
University for sharing their knowledge and time.

My deep gratitude goes to my parents their warm love, support, and
encouragement during my study in East Lansing. I am also thankful to all my friends in

Korea and MSU for their suggestions and kind help during my study.

iii

TABLE OF CONTENTS

ABSTRACT ..................................................................................................................... ii
ACKNOWLEDGEMENTS ............................................................................................ iii
LIST OF FIGURES ........................................................................................................ vi
LIST OF TABLES ........................................................................................................ viii
1. INTRODUCTION ......... 1
2. EDDY CURRENT TECHNIQUE ................................................................................ 5
2.1. Principles of Eddy Current Technique (ECT) ............................................... 5

2.2. The Eddy Current Phenomenon ..................................................................... 8

2.3. Skin Effect ................................................................................................... 11

2.4. Multi-Frequency Eddy Current Testing ....................................................... 13

2.5. Analysis of Multi-Frequency Eddy Current Testing (MFECT) .................. 17

2.5.1. Algebraic Method ......................................................................... 18

2.5.2. Phase Rotation and Subtraction Method ....................................... 21

3. INDEPENDENT COMPONENT ANALYSIS (ICA) ................................................ 28
3.1. Principles of ICA ......................................................................................... 28

3.2. Preprocessing ............................................................................................... 32

3.3. Independent Component Analysis Based on Gaussianity ........................... 33

3.4. Independent Component Analysis Based on Mutual Information ............... 44

4. APPLICATION OF ICA TO NDT ............................................................................. 47
4.1. Feature extraction algorithm ........................................................................ 47

4.2. Finite Element Method (FEM) ..................................................................... 51

5. SIMULATION RESULTS AND DISCUSSIONS ..................................................... 59

5.1. Simulation of FEM Data .............................................................................. 60
5.1.1. Source Signals Separation (Defect and TSP) ................................ 62
5.1.2. Source Signals Separation (Defect, TSP and Noise) .................... 66
5.2. Performance on Experimental Data ............................................................. 69

5.2.1. Extraction of the Defect Signal from Mixture of Defect

and Noise ................................................................................................ 69

5.2.2. Extraction of the Source Signal from Defect and TSP ................. 73

5.2.3. Extraction of the Source Signals from Defect and Noise ............. 77

5.2.4. Use of ICA for defect detector ...................................................... 82

6. CONCLUSION AND FUTURE WORK ................................................................... 84
6.1. Conclusion ................................................................................................... 84

6.2. Future Work ................................................................................................. 84
References ....................................................................................................................... 86

LIST OF FIGURES

Figure1.1 Schematic of a Nuclear Power Plant ........................................................................ 2
Figure 2.1 Eddy Current Principle ............................................................................................ 6
Figure 2.2 Simplified Model of ECT Probe ............................................................................. 7
Figure 2.3 Steam generator tube with support plate ............................................................... 14
Figure 2.4 Cutaway View of steam generator tube with support plate and defects ............... 15
Figure 2.5 Simulated Eddy Current Impedance Plane Trajectories due to: ........................... 16
Figure 2.6 Impedance Plane Trajectories of a defect close to a support plate ....................... 16
Figure 2.7 Architecture of RBF neural network ..................................................................... 24

Figure 2.8 Normalized Simulated Impedance Plane Trajectories of Defect at excitation

frequency: (a) 50 kHz, (b) 75 kHz, (c) 100 kHz, (d) 125 kllz, (e) 150 kllz, and Support Plate at

excitation frequency (f) 50 ltHZ, (g) 75 Ht, (h) 100 kHz, (i) 125 kllz, (i) 150 kllz ......................... 27
Figure 3.1 Schematic overview of ICA Process ..................................................................... 31
Figure 3.2 Flowchart of Fast—ICA Algorithm ........................................................................ 42
Figure 4.1 Application of ICA to NDT ................................................................................... 47
Figure 4.2 Schematic diagram of ICA analysis in NDE ......................................................... 49
Figure 4.3 Impedance trajectory plots of ECT data after RST ............................................... 51
Figure 4.4 Geometry for eddy current NDT [34] ................................................................... 53
Figure 4.5 A finite element mesh for eddy current geometry ................................................ 53

vi

Figure 5.1 ECT data from FEM at frequency of 35 kill (a) Impedance trajectory (b) Real
component of eddy current signal (upper) and imaginary component of eddy current signal
(lower) ..................................................................................................................................... 61
Figure 5.2 ECT data from FEM with excitation frequency of 200 kill (a) Impedance
trajectory (b) Real component of eddy current signal (upper) and imaginary component of
eddy current signal (lower) ..................................................................................................... 61
Figure 5.3 ECT data from FEM with excitation frequency of 400 kill (a) Impedance
trajectory (b) Real component of eddy current signal (upper) and imaginary component of
eddy current signal (lower) ..................................................................................................... 62
Figure 5.4 Impedance trajectories using FEM at (a) 35 kllz (b) 400 ldiz (c) separated defect
signal ((1) separated support plate signal ................................................................................. 64
Figure 5.5 Real component (Upper) and imaginary component (Lower) of (a) FEM—ECT
signals at 35 kllz (b) FEM—ECT signal at 400 kHz (0) extracted defect signal (d) extracted
support plate signal .................................................................................................................. 65
Figure 5.6 Source signal separation with noise filtering (simulated trajectories): (a) Noisy
ECT data at 35 kHz, (b) Noisy ECT data at 200 kHz, (0) Noisy ECT data at 400 kHz, ((1)
Separated Defect Signal, (e) Separated Support Plate Signal, (f) Separated Noise ............... 67
Figure 5.7 Source signal separation with noise filtering (simulated data): (a) noisy horizontal

ECT data at 35 kllz (top), 200 kl’a(middle) and 400 kllz(bottom), (b) noisy vertical ECT data at
35 kllz (top). 200 kllz(middle) and 400 kllz(bottom), (c) separated signals (horizontal

component), ((1) separated signal (vertical component) ......................................................... 68

vii

Figure 5.8 Test 1 —ICA implementation with ECT field data: (a) Horizontal ECT data, (b)
Horizontal Separated Signal, (c) Vertical ECT data, (d) Vertical Separated Signal .............. 71
Figure 5.9 Test 1 — Impedance trajectories of ECT field data at frequency: (a) 35 kHz, (b)
200 kllz, (c) 400 kHz, ((1) 600 kl‘b ................................................................................................. 72
Figure 5.10 Test 1 — Impedance trajectories of the reconstructed signal: (a) Defect, (b) Noise73
Figure 5.11 Test 2 — ICA implementation with ECT field data: (a) Horizontal ECT data, (b)
Vertical ECT data, (c) Separated horizontal signals, (d) Separated vertical signals .............. 75
Figure 5.12 Test 2 — Impedance trajectories of ECT field data at frequency: (a) 35 kllz, (b)
200 kHz, (c) 400 kHz, ((1) 600 kllz ................................................................................................. 76

Figure 5.13 Test 2 — Impedance trajectories of the extracted signals: (a) Defect, (b) TSP....77
Figure 5.14 Test 3 — ICA implementation with ECT field data: (a) Horizontal ECT data, (b)

Horizontal Separated Signal, (0) Vertical ECT data, ((1) Vertical Separated Signal .............. 79

Figure 5.15 Test 3 — Impedance trajectories of ECT field data at frequency: (a) 35 kHz, (b)

200 kill, (0) 400 kHz, ((1) 600 kllz ................................................................................................. 80

Figure 5.16 Test 3 — Impedance trajectories of the reconstructed signal: (a) Defect, (b) Noise81

viii

LIST OF TABLES

Table 5.1 Parameters for Finite Element Model ..................................................................... 60
Table 5.2 Statistics on the ICA performance .......................................................................... 82
Table 5.3 Statistics on the ICA performance .......................................................................... 83
Table 5.4 Statistics on the mixing algorithm performance ..................................................... 83

1. INTRODUCTION

Steam generators (SG) are heat exchange units in nuclear power plants (NPP) [1]. The
function of the steam generator is to transfer the heat generated in the primary loop to a
mixture of water and steam circulating outside the tube. The water around the tubes absorbs
the heat and vaporizes into steam which is used to run the turbines. Ferromagnetic support
plates in the steam generator located at periodic intervals are used to anchor the tubes. The
steam generator tubes and support plates when exposed to extreme environments such as
high temperature and pressure can result in corrosion products which are deposited in the
crevice between tubes and support plates. This, in turn, leads to denting, thinning and
cracking of the tube. When cracking occurs in the tubes, radioactive primary water can leak
and contarrrinate the water and steam used to drive the turbines. Hence, it is very critical to
find small flaws in the SG tubes before they develop into a propagating crack.

It is critical that there is no leakage of fluid from the primary loop into the secondary
loop. However, the harsh environmental conditions cause the tubes in the SG unit to be
subjected to various types of degradation mechanisms such as mechanical wear between
tube and support plate, outer diameter stress corrosion cracking (ODSCC), pitting,
volumetric changes, primary water stress corrosion cracking (PWSCC) and inter granular
attack (IGA), all of which result in tube thinning and cracking. Since the outage of a NPP
costs almost $500,000 a day, there is a strong motivation to develop a reliable method to
assess the safety of the unit. The steam generator tubes are inspected periodically for

cracks, and degradation.

Shield

Steam
Generator Seconda ry LOOP Electricity
Turbine Generator

rillllllfl'mm
..iriim rile.

mu-mum-miu

        
 
   
  

Pressurize

 

rriruullrun-muummr

.illll It.

 

Condensing Chamber

  
      
   

Primary Lccp

Ocean

Figurel.ll Schematic of a Nuclear Power Plant

Eddy current technique [2,3] is a widely used practical tool for inspecting SG tubes.
Although the measured eddy current data contain information related to defects, it is
sometimes corrupted by noise and contribution from external support plates that make the
analysis challenging. This problem is usually solved by means of a rrrixing algorithm [4]
implemented using multi-frequency eddy current data where a signal at a low excitation
frequency is transformed to the corresponding signal at high frequency and the transformed
signal is subtracted from the measured high frequency signal to suppress the tube support

plate (TSP) contribution.

This thesis investigates the use of Independent Component Analysis (ICA) to separate
out the signals due to defects and tube supports from noise. The measured data is modeled
by a combination (linear or non-linear) of the sources that is represented using a mixing
matrix. The objective of ICA is to invert the unknown mixing system and to estimate the
individual source signals. In general, the mixing and the separation problem can be
formulated as

§=Wx where x=As (1.1)

where s is source signal, x is measurement, § represents the estimated source signal.
A is mixing matrix and W is separation matrix that is the inverse of A.

Under the assumptions that (i) all sources are mutually independent, (ii) the signals are
linearly mixed and (iii) no individual source signal is Gaussian, the separation system can be
obtained by maxirrrizing the mutual independence between components.

ICA was originally developed to solve the cocktail party problem of extracting speech
signals of individuals from data obtained by recording sounds from a roomful of people. It
was originally addressed as the problem of blind source separation by Christian Jutten and
Jeanny Herault in the early of 1980‘s [5, 6]. In the context of blind source separation, the
inverse problem has been well studied and many algorithms have been developed depending
on the nature of the mixing model. The problem of blind source separation is more difficult
since the mixing matrix is unknown. The only assumption made by Jutten and Herault was
the independence of the source signals. A fast and efficient algorithm was required and Bell
and Sejnowski obtained the result based on infomax approach to ICA [7]. S. Amari soon
realized that the infomax ICA algorithm could be improved by using the natural gradient,

which multiplies the gradient of the feed forward weight matrix by a positive definite matrix

and makes it converge faster by eliminating the matrix inversion [8]. This allowed infomax
ICA algorithm to be practical for various problems in the real world. In 1998, Te-Won Lee,
Mark Girolami and Sejnowski developed the extended version of infomax ICA algorithm
that is suitable for general non-Gaussian signals [9]. Finnish researchers, A. Hyv'arinen and E.
Oja, developed the fast and robust metric known as fast ICA algorithm for fixed-point data
[10, 11].

This thesis is organized into 6 chapters. Chapter 2 describes the physical principles of
eddy current technique and multi-frequency eddy current (MFEC) data while Chapter 3
presents the primary approach proposed in this thesis, namely independent component
analysis (ICA). Chapter 4 introduces the application of ICA to ECT data and finite element
method (FEM) for NDE. The simulation results are presented and discussed in Chapter 5.

Conclusions and directions for future work are given in Chapters 6.

2. EDDY CURRENT TECHNIQUE

2.1. Principles of Eddy Current Technique (ECT)

Eddy Current Method is a widely used nondestructive testing technique in nuclear,
aerospace, power, transport, petroleum and other industries to inspect conducting
samples for detection and characterization of defects, corrosion and other variations. ECT
is based on electromagnetic induction and hence direct electrical contact with the
material is not required. ECT is widely used largely, because it is a non-contact method
and can detect surface and sub-surface defects (fatigue cracks, stress corrosion cracks)
with hi gh-speed, accuracy and reliability.

ECT is principally based on the Maxwell-Ampere law and Faraday’s law of
electromagnetic induction. When a coil is excited by alternating currents, a time-varying

magnetic field is set up in accordance with the Maxwell-Ampere law [12].
jCHodfzjjS J-dE (2.1)

When the coil is brought close to a conducting material, the primary field associated

with the coil induces EMF in the medium according to the Maxwell-Faraday law.
§C E -di =— ”8:3 . at (2.2)
t

which produces eddy currents that flows in closed paths and it is shown in Figure 2.1.
According to Lenz’s law, the EMF and induced currents are directed so as to oppose
the change that produces them. The magnetic field set up by the induced eddy current
opposes the primary magnetic field associated with the coil. When the test specimen is
nonferromagnetic the net flux linkages of the coil decreases which in turn decreases the

inductance of the coil. Accompanying the decrease in the inductance is an increase in

resistance of the coil since the eddy current losses incurred in the specimen has to be met
by the source of primary excitation. The presence of discontinuity in the material results
in the different magnetic flux linkages. Therefore, the changes in coil impedance due to
the presence of flaws can be analyzed to estimate the surface properties of the specimen

[13].

 

Alternating Current

 

Material

 

 

 

Figure 2.1 Eddy Current Principle

Eddy current signals are displayed as horizontal and vertical channel data of the
complex impedance in the impedance plane and trajectory. The eddy current probe data

can be measured using an AC bridge as shown in Figure 2.2.

 

I l
VR .

v, T

 

Oscillator I

 

 

Figure 2.2 Simplified Model of ECT Probe

In a direct current (DC) circuit, current and voltage can be represented by the magnitude.
However, the analysis of alternating current (AC) circuit is more complicated. Since the
magnitude of current and voltage vary with time, not only the magnitude information but
the phase information must be taken into account. Assume that a resistor and inductor are
connected in series with an oscillator (Figure 2.2). This is a typical AC circuit and can be
considered to be a simplified model of probe. The inductor models the reactive part of the
coil while the resistor models both coil wire and cable resistance. For this circuit, the voltage
across the inductor (VL) leads the current (I) by 90° while voltage across the resistor (VR)
has the same phase with the current. Therefore, the current can serve as a point of reference.
As a result, the voltage across the inductor leads the voltage across the resistor by 90°. The
voltage across the resistor and inductor (VL) leads the current (or the voltage across the
resistor) by an angle less than 90°.

To evaluate the total voltage (VT), we add the voltage across the resistor (V R) and the

inductor (VL).

VT = VR + VL
= I(R + ij) (2-3)
= IRsin(wt + 0) +ijsin(wt + rr/2)

Hence, the impedance is represented as

V
Z 2 —Il = Rsin(wt) + ijsin(wt + 1r/2) (2.4)

The two terms in the last equation (2.3) contain the amplitude and phase shift. Hence,
they can be represented by phasors. The amplitude of the first term is R and the phase

shift is 0 while the amplitude of the second term is wL and the phase shift is 1r/2.

2.2. The Eddy Current Phenomenon

The governing equation for eddy current phenomenon can be derived from Maxwell‘s

equations [12, 14]. In differential form, the Maxwell‘s equations can be written as

_33

V x E = 3’— (2.5)
V x H =J + 8—D (2.6)
a:

V o B = 0 (2.7)
V o D = p (2.8)

and the constitutive relations for isotropic, linear and homogeneous medium are
B = p H (2.9)
D=6 E (2.10)
J = a E (2.11)

where
5 is the electric permittivity (farads/m)
,u is the magnetic permeability (Henry/m)

a is the electric conductivity (mhos/m)

Since V- B=O , B can be expressed as the curl of the vector magnetic potential A
given by
B = V x A (2.12)

Substituting for B in equation (2.5)

Vsz—Vxéi
at
or Vx£E+a—A—)=0 (2.13)
at
Using the vector identities we express
E+%A=-V(p (2.14)
t

where (o is the scalar electric potential.

Substituting equation (2.11) in equation (2. 14) we have

8A 8A
J “(tyywwi 41—87111, (2.15)

where J S =—0V¢ is the source current density and 0%]: is the induced eddy current
t

density.

. . . . . 8D .
At the excrtatron frequencres for eddy current testing the displacement current a— in
r

equation (2.6) is negligibly small incomparison with the conduction current density J
and equation (2.6) reduces to

VXH=J (2.16)

Substituting equations (2.9) and (2.15) into (2.16)

[Vx-Ejals—O'EA— (2.17)
[1 3:
Using equation (2.12)
1 8A
[Vx—(VXA)]=JS —0'— (2.18)
,u 3:

Assuming a homogeneous medium and using the vector identity

Vx(VxA)=V(VoA)—V2A (2.19)
we arrive
iv2A=a§i‘——JS (2.20)
,u at

where V o A = 0 choosing the Coulomb gauge.

Assuming that the fields vary harmonically in steady state we can express A as
A = Aoe’j‘" (2.21)
where w is the angular frequency.

Substituting (2.21) into (2.20), we obtain

iv2A= jwaA-JS (2.22)
,1:

The solution of equation (2.22) can be obtained analytically only for very simple

geometries such as infinite half plane media.

10

2.3. Skin Effect

Eddy currents are closed loops of induced currents circulating in planes perpendicular

to the magnetic flux. They are oriented parallel to the coil winding and are limited to the

area of the inducing magnetic field. Excitation frequency determines the depth of

penetration into the specimen; as frequency increased, the depth of penetration decreases

and the eddy current distribution becomes denser near the specimen’s surface. This is

known as the skin effect. This implies that test frequency also affects the sensitivity to

changes in material properties and defects.

The depth of penetration depends on the operating frequency and properties of material

such as electrical conductivity and magnetic permeability. The value of the skin depth can

be determined from Maxwell’s equation as follows [13].

Using equation (2.11), equation (2.16) can be written as

VxH=aE
Taking the curl on both sides
Vx (V x H)=Vx (0E)
Using the vector identities in equation (2.19), equation (2.24) can be written as
—V2H+Vx(VoH)=Vx(0E)
From equation (2.7), V o H = 0 and hence equation (2.25) reduces to
V 2 H = — V x (0E)
From equation (2.5) and (2.9), we have

VxE=—-a£=—,u-a—H-
a: 3:

Substituting equation (2.27) into equation (2.26)

11

(2.23)

(2.24)

(2.25)

(2.26)

(2.27)

vzrr = poi"; (2.28)

Assuming the field varies harmonically in the steady state with time H can be

expressed as
H =H0e’” (2.29)
where w is the angular frequency.

Substituting equation (2.29) into (2.28), we have

VZH =(jwp0')H (2'30)
=k2H
where k2 = jwtra.
Foran infinite sheet of current in the Y direction on the YZ plane the magnetic field
intensity is in the Z direction with no components in X or Y directions. Hence,

equation (2.30) reduces to

82H
‘aTzz'JKZHZ (2.31)

And the solution to equation (2.31) is given by
Hz = Hoe—l“ (2.32)

where k is defined by equation (2.30), called the propagation constant, is given by

k = Gwyn)“
1/2 1/2 (2.33)
0’ . 0'
=14) + (4)
Hence, equation (2.32) can be written as
Hz = Hoe—We—jw (2.34)

12

 

where 5 = that is called as the skin depth.

1
Jam
The standard depth of penetration is defined as the depth at which eddy current density

has decreased to l/e or 36.8% of the surface value. The skin depth (d) is often used as a

guideline to select the excitation frequency for a given test specimen. It is computed as

 

(5': (2.35)

where f is excitation frequency, ,u is magnetic permeability of material and o is electrical
conductivity of the target material. Although the induced fields and currents can reach
deeper than one standard depth of penetration, the eddy current density decreases very

rapidly. At the two standard depths of penetration (26), the eddy current density decreases to

(He) 2or 13.5% of the surface density. At the three standard depths of penetration (3(5), it

decreases to 5% of the surface density.

2.4. Multi-Frequency Eddy Current Testing

Multi-Frequency eddy current testing is usually used in non-destructive evaluation
where the defect signal is distorted, due to superposition of contribution from other
structures in the test geometry. Artifacts, such as specimen geometric boundaries or attached
conductive objects, generate their own impedance plane trajectories during eddy current
testing. A major problem in the analysis of SG tube inspection is the corruption of defect
signal by a large contribution from tube support plate (TSP).

As an example, consider the tube and support plate shown in Figure 2.3. This

configuration is typical of coolant tubes found in pressurized steam generators, such as

13

those used in nuclear power plants. The tube projects through support plates arranged within
the steam generator [15]. Eddy current testing of these tubes is conducted to determine their
structural condition. Specifically, the tubes are examined for corrosion pitting, dents and

build-up of sludge [16].

Support Plate
(Carbon Steel)

 
  

Steam generator tube

(Inconel) \

Figure 2.3 Steam generator tube with support plate

   
   
    

Support Plate
“x

\

\\

Defect B Defect A
\

\
\

Eddy Current
Probe

Figure 2.4 Cutaway View of steam generator tube with support plate and defects

The problem of distortion of defect signals in single-frequency eddy current testing is
illustrated in Figure 2.4 that shows a cut-away view of Figure 2.3. There are two defects,
A and B , in the tube. As the Eddy current probe passes through the tube, the impedance
plane trajectory of defect B will be visible because there is no structure in the vicinity.
However, the trajectory generated by defect A will be distorted, since the support plate is
a conducting material and will also generate an eddy current signal trajectory as the probe

passes it. The superposition of the support plate’s signal with that of defect A makes the

resulting trajectory more complicated. This distortion makes defect detection difficult.

15

 

 

 

Figure 2.5 Simulated Eddy Current Impedance Plane Trajectories due to:
(a) defect (b) support plate

0.5 r

 

.1
q
--—4
—1
—1
—(

 

 

 

 

Figure 2.6 Impedance Plane Trajectories of a defect close to a support plate

16

-‘1'"'"h.-_ _"'" I. u— -

 

Trig-4w"-
‘5’,“

Figure 2.5 shows simulated data representing eddy current trajectories generated by
defect in tube and support plate and Figure 2.6 shows the impedance plane trajectory when
the defect is in close proximity to the support plate. In addition to this distortion shown in
Figure 2.6, the presence of noise also deteriorates the impedance plane trajectory of a defect,
making data interpretation difficult. Filtering of the eddy current signal may eliminate much
of the problems due to noise, but does not make the task of extracting defect information
any easier due to the problem of superposition.

Multi-frequency Eddy current technique (MFECT) is commonly used for solving this
problem. In MFECT, two or more sinusoidal signals of different frequencies are fed into a
single eddy current probe and the gain and the phase signal from each frequency is collected
separately. As it is shown in the previous section, the depth of penetration and phase lag are
a function of test frequency; increasing test frequency reduces the depth of penetration and
increases phase lag. Hence, it is possible to change the eddy current response by changing

the test frequency and use this information to remove the unwanted signals.

2.5. Analysis of Multi-Frequency Eddy Current Testing (MFECT)

The aim of multi-frequency eddy current testing is to suppress artefacts affecting the
defect impedance trajectory. This method of defect discrimination begins with identification
of desired and undesired features. Desired features are typically defects and undesired
features are artefacts that make eddy current characterization of defect difficult. When a
high frequency is used as excitation frequency, inner diameter (ID) defect are more
predominant in the signal. At intermediate frequency, all features in the wall are detectable

and there is phase discrimination between internal and external defect signals. At low

17

frequency, there is a predominant signal from TSP with little phase separation between
internal and external defect signals. In MFECT, signals are collected separately at a low and
high frequencies. The high frequency ECT signal is rotated to match the low frequency ECT
signal. The amplitude of high frequency ECT data is also scaled to match the low frequency
ECT signal. On subtracting the transformed high frequency signal from the low frequency
signal, the signal from TSP is largely reduced and only defect signal will remain as output
signal. This process is also referred to as mixing.

Signal processing methods employed in analyzing MFECT data can be classified into
three classes. Specifically, these are the algebraic, coordinate transform and phase rotation
and subtraction method. All the methods assume that the signal contribution from the defect
and the artefacts superpose linearly to form the eddy current impedance trajectory. This is
often an erroneous assumption. Even though it is reasonable to assume linearity with respect
to the excitation function at low excitation level, the system response is not linear with
respect to the material property distribution in the domain of interest. Nevertheless, under

some restrictive conditions, linearity may be assumed [17].

2.5.1. Algebraic Method

The algebraic method assumes that the multi-frequency data are linear function of n
parameters. This assumption has been shown to be valid experimentally, under small signal
conditions [18]. Furthermore, these parameters represent n inhomogeneities or attributes
within the specimen [19]. To maintain the validity of this model, small signal conditions are
assumed such that there is a linear relationship between variations in the parameters and the

probe coil impedance. The complex impedance phasor of the test coil is represented by its

18

real (x-axis) and imaginary (y-axis) values. In case of two frequencies testing of four
parameters, the coil impedance is expressed as

x1 =a11P1+012P2+a13P3+ar4P4
Yr =a21Pr +“221,241531334112413:
x2 =031Pr +032P2 +“33133 +034P4
Y2 :a4rPr +042P2 +“431,3 +0441”),

where
. . (2.36)
x, : ohmrc resrstance at frequencyl
y , :inductive reactance at frequencyl
x2: ohmic resistance at frequency2
y2 :inductive reactance at frequency2
Pi : Observed specimen parameters
a”. :Coefficient of jth parameter
In matrix notation, it is expressed as
AP = C
where (2.37)
l- " f" " F -'
an an 013 an Pr x1
(1 a,, a a, P y
A = 21 ..1. 23 .4 , P = 2 and C = l
0‘41 032 an “34 P3 x2
_a4r “42 “43 “44‘ _P4‘ _y2_

 

 

 

 

 

 

In this model the matrix A is unknown and only the elements of C are observed
directly. Elements of P are measured prior to eddy current testing, during a calibration
stage. Equation (2.37) shows the modulation of the parameter vector P by a modulation
matrix A , generating a signal vector C. The modulation equation is assumed to represent
the interaction between the system current flow within the test coil and the specimen
characteristics [18].

The matrix A is determined in the calibration stage using calibration standards with

19

[€31

COB

ltlit

llllfi.

eltn

mail

known elements of P, and observing the elements of C. Once the matrix A has been
determined, it is applied to the signal vector C during eddy current testing of a specimen.
The goal of the multiparameter method is to detect and isolate the actual parameters by
separating or decoupling the system of equations represented by equation (2.37). A complete

solution has the form
P=A"C (2.38)
The signal vector elements represent the in-phase (x) and quadrature (y) of eddy
current data collected at different frequencies. Specifically, for the case of two frequencies
testing, x1 and )21 are equivalent to the in-phase and quadrature coil impedance
components at frequency 1, respectively. In the case of four parameters, a similar
relationship holds for elements x2 and y2 at frequency 2.

Since decoupling the system of equations represented by equation (2.37) is of primary
interest, a complete solution to equation (2.37) is not necessary. The decoupling of the

elements of P is accomplished by the adjoint of matrix A , contained within the inverse
matrix A‘I [21]. It reduces equation (2.38) to
P=adjiA“JC (2.39)

Therefore,

20

 

Pr =b11xr+b12YI+br3x2 +brtryz

P2 :bzrxr +b22y1+b23x2 +b24)’2

P3 =b31x1+b32y1+b33x2 +b34Y2

P4 =b4lxI +b42y, +1943x2 44944))2

where (2.40)
x, , y , :frequency 1 data elemenets of C

x2 , y2 :frequency 2 data elements of C

Pi : Observed specimen parameters

by. :Coefficient of j‘h parameter of adj[A]
The coefficients by are determined using equation (2.39), or by using iterative

numerical methods [19].

2.5.2. Phase Rotation and Subtraction Method

The objective of the phasor rotation and subtraction method is to discriminate the defect
by subtracting out the effect of the artifacts. For suppression of unwanted artefact, an
auxiliary frequency is utilized. This frequency is selected such that the eddy current probe is
particularly sensitive to the presence of the undesired artefact. Due to the skin effect, low
excitation frequencies are used to obtain signals that are sensitive to support plate and high
excitation frequencies are used to obtain signals that are sensitive to probe wobble.

The subtraction of these undesired effects occurs in the mixing stage, after detection of
the phase and magnitude of the probe coil signal. At this stage of signal processing, the eddy
current data is resolved to the in-phase and quadrature components of the impedance plane.
By transforming the auxiliary frequency trajectory to resemble its primary frequency
trajectory, the data required for subtraction from the defect signal is obtained. The

transformation of the auxiliary frequency artifact trajectory is accomplished through rotation,

21

 

translation and scaling. This procedure is typically executed manually, where an NDE
analyst observes the stored complex impedance plane primary and auxiliary frequency
trajectories due to the tube support plate.

The input data, xa and ya are applied to the scaling stage, which applies the
transformation parameters SX and S y. The output of this step is applied to 6 rotation

combined with the translation parameters Tx and Ty . Then, the process is complete when

the transformed auxiliary artifact signal is subtracted in the mixing stage from the distorted
primary frequency defect trajectory. The output of typical two-parameter transformation is
x, =xp —l(Sxxa cost9+Syya sin Qj—ij

y, =yp -[(—Sxxa sint9+Syya cosB)—Ty]

where

x p , y p :Input primary frequency data ordered pair (2 41)
xa , y 0 :Input auxiliary frequency data ordered pair ‘

x, , y, :Mixed output data ordered pair
S x,S y:Scaling parameters

Tx , T y :Translation parameters

6: Rotation parameter, degrees

Least Square Estimation
With properly chosen transformation parameters in equation (2.41), the effect of the
artifact is minimized and the defect signal is more clearly defined. The optimal affine

transform matrix parameters are estimated using a least squares estimation scheme [22]. By

defining 5,, as the primary frequency signal vector and Ca as the auxiliary frequency
signal vector, the transformed Ca is expressed as
E; :6“ .A(Sx,sy,rx,ry,e) (2.42)

22

The cost function E ,5 is expressed as

 

 

(2.43)

 

 

where H“2 is the Euclidean norm.

The optimal values of the affine transform parameters are estimated by minimizing the

cost function. This involves the solution of a simultaneous system of five equations

Egg—i=0 (2.44)

:i—tw (2.45)

%=o (2.46)

95%.:0 (2.47)
y

95%:0 (2.48)

Solutions of these equations yield the optimal affine transformation parameters. One of
the iterative approaches to minimizing the cost function is the method of steepest descent.
By choosing initial values for all five affine matrix parameters, the values are stepped
iteratively so as to minimize the cost function. At each step, the gradient of the cost function

is evaluated at the current parameter values.

Radial Basis Function Method

23

lll‘c‘
Wei

lleg

toil

In mixing, a signal at a low excitation frequency is transformed to the corresponding
signal at high frequency and the transformed signal is subtracted from the measured high
frequency signal [4, 22]. In this approach, a radial basis function (RBF) neural network is
applied to obtain transformation parameters. The architecture of RBF neural network is

shown in Figure 2.7.

Hidden layer of radial

basis functions Output layer

Input layer

 

 

 

Figure 2.7 Architecture of RBF neural network

The RBF network output is a function F that has the following form.
N
th.) = >3w.¢l|x-x.-Il) (2.49)
i=1

where xi is the input at the ith input node, ¢(0) is a set of arbitrary functions known as
the radial basis functions, "0" denotes Euclidean norm and w,- constitutes a set of linear

weights. During calibration, the set of weights are computed using known signals at

frequencies f1 and f2.
Using one of these methods described above, the signal at frequency fl is transformed

to the corresponding signal at frequency f2. The residual is obtained by subtracting

24

transformed signal from the signal at frequency f2 that contains only the defect signal with
the contribution from artifacts suppressed to a large degree.

However, such approaches have significant errors when the data is not similar to data
used to generate the mixing parameters. Consequently the parameters of mixing algorithm
depend strongly on the data in the standard calibration file. When the parameters derived
using calibration data is not optimal for a test signal the complementation of ICA results in a
large residual which is difficult to interpret. In this thesis, an alternate approach to linear
mixing algorithms based on the use of independent component analysis (ICA) is investigated
to separate out signals due to defects and external support plates from noise. ICA is a method
that depends on only the observed signal. In this thesis it is has also been used for noise
filtering. The theoretical basis of independent component analysis is discussed in the next

chapter.

25

III'IIIIIIIIIIIIIIIIIIIIIIIII'I'I.

 

 

 

 

I'liIl-IlllIIIIII'IIIIIIIIIIIJ

 

 

 

 

lilllltiilll'llllll'tlrtilllilll.

 

 

 

 

 

(h)

(C)

26

 

 

 

 

 

 

(d)

 

 

 

 

(l)

 

frequency: (a) 50 kHz, (b) 75 kllz, (c) 100 kHz, ((1) 125 kllz, (e) 150 kHz, and Support Plate at excitation

Figure 2.8 Normalized Simulated Impedance Plane Trajectories of Defect at excitation
frequency (0 50 kHz (g) 75 kHz, (h) 100 kHz, (i) 125 kllz, (j) 150 kllz

27

3. INDEPENDENT COMPONENT ANALYSIS (ICA)

3.1. Principles of ICA

Independent Component Analysis (ICA) is a statistical method for extracting the source
signals from a set of observed signals that can be expressed as mixtures of multiple sources.
This problem is analogous to the human ability to distinguish a certain sound in the presence
of different kinds of sounds. For instance, a human being can distinguish the sound of a
violin from the recording of music by an orchestra. Hence, it can be said that ICA attempts
to mimic this ability.

The data from a set of two source signals 3, and 52 can be represented by a set of

simple linear equations.

x1=ausl+a1252 (31)
x2 = 02131+02232

where a,,,a,2,a2, and a22 are mixing parameters that depend on the experimental
conditions, x, and x2 are the measured signals from a mixture of 2 source signals, s,
and s2.

When the coefficients a 1 ,,a , 2,a 2 1 and c222 are known, it is straightforward to find the

source signals from the above equations using classical inversion methods. In real world
problems, however, these coefficients are seldom known which makes the solution of these
equations very complex. The Independent Component Analysis method based on the central
limit theorem provides a technique for achieving this goal [23]. The central limit theorem
states that, as the number of independent, identically distributed random variables increases,

the cumulative density function (CDF) of the sum will approach the CDF of a Gaussian

28

random variable. It is assumed that the source signals are independent and identically
distributed in the rest of this thesis. Hence, measuring the gaussianity of a signal and
minimizing the gaussianity provides a way of finding individual components. A major issue
with this approach is selection of a metric that measures the gaussianity of a signal. The
details of the method for measuring gaussianity will be discussed in this chapter.

The general statistical model is defined below. Suppose that there are n linear mixtures
and n independent sources. Each of the mixtures is assumed to be a linear weighted sum
of the original signals from n sources.

x]: 01131 +012S2 +~-+alnsn

x2 = 02151 + 02252 +‘“+ (121131:

(3.2)
Xj =ajlsl +aj252 +"'+ajnS

x" = an1s1+anzsz +---+a,msn
The time index, I is omitted because x and s are regarded as random variables and
x(t) is obtained as samples of the random variable. The observed signal )1 and original
signal s are assumed to have zero mean without loss of generality. Otherwise, a pre-
processing step can be applied to subtract the sample mean to make them zero-mean
random variables.

The equations (3.2) can be represented in matrix form as

x 2 As (3.3)
where x=[x1 x2 xnjr and s=[s, s2 snjr.
The above ICA model is a generative model, which means that the observed data are

generated by a process of mixing the components 8 using the mixing matrix A which is

29

unknown. The aim of independent component analysis is to estimate the unrrrixing matrix
W using statistical principles and computations. Eventually, the original source signals can

be found as

m)
ll

2

X

(3.4)
where s is the estimated source signal.

The ICA method estimates the source signals up to a constant. This is because both the
mixing matrix and original signals are unknown and any scalar multiplier can be cancelled
by dividing the corresponding column by the scalar. The sign of the signal is ambiguous for
the same reason. Another characteristic of ICA is that the order or position of the original
signal cannot be determined. The order of source signals can be changed freely during the
process and any component can be the first source because both the mixing matrix and
original signals are unknown. Figure 3.1 shows the schematic overview of the ICA
algorithm. The input data collected consists of a number of recordings which is larger than
the number of expected source signals. The next two steps consist of preprocessing and ICA
algorithm which are investigated at chapter 3.2 and 3.3, respectively. At the output stage, we
obtain the estimated source signals using the unmixing matrix found using the ICA

algorithm.

30

 

INPUT

Observed Random vector

 

 

 

 

 

Preprocessing

 

 

 

 

 

Fast ICA algorithm

 

 

 

 

Output

 

 

Unmixing Matirx

 

 

Figure 3.1 Schematic overview of ICA Process

31

3.2. Preprocessin g

A preprocessing stage is required before applying ICA for better conditioning of the
input data. Generally it consists of centering and whitening the data. Centering is the
simple process of subtracting mean value from original data to obtain a zero mean signal.
After estimating the original signals, it is necessary to add the subtracted mean value to the
estimated signals. The second step is whitening of the observed variables. The purpose of
whitening is to obtain a new vector in which the components are uncorrelated and variance

of each signal is unity. The most popular method for whitening is using eigenvalue

decomposition (EVD) of covariance matrix E{xxT} =EDET where E is the
orthogonal matrix of eigenvectors of E{xxT} and D is the diagonal matrix of

eigenvalues of E{xxT} [24]. The whitening transform is implemented as

-1
x'= ED /2ETx (3.5)

Using the whitening process, the orthogonal mixing matrix A' can be computed as

-1 _1
x'=ED AETx=ED /2ETAs=A's (3 6)
A’ = ED’%ETA
By using the above transform, the orthogonal mixing matrix A' is obtained. When the

variance of independent component is assumed to be unity, it can be seen that
1:.{11'11'T } = [scramT }A'T = A'A'T =1 (3.7)

Because A' is orthogonal matrix, A'T is equal to the inverse of A'. Hence, the
transformation reduces the number of parameters that need to be calculated to find the

unmixing matrix. When the dimension of matrix is large, the whitening step reduces the

32

computational time and storage considerably. As a result of whitening, the measured data is
also transformed to be uncorrelated and have a unit variance. In this thesis, it is assumed

that centering and whitening is performed on the measured signals.

3.3. Independent Component Analysis Based on Gaussianity

The main concept of ICA is inspired by the central limit theorem which states that, as
the number of independent, identically distributed (i.i.d.) random variables increases, the
cumulative density function (CDF) of the sum approaches the CDF of Gaussian random
variable. Under the assumption that the different source signals are independent to each
other, the original signal can be reconstructed when minimum gaussianity of the random

variable is observed. Let x ,,x2,...,x,, be a sequence of i.i.d. variables and let s" be the

sum of n random variables.
sn=xl+x2+...+xn (3.8)
No matter what the distributions of random variables x 1.x2,. . . ,x,, are, the CDF of sn
approaches the Gaussian distribution as it increases. Therefore, if gaussianity is measured
using some method, it will be a monotonically increasing function of the number of
independent components. On the other hand, gaussianity will be minimized when 3,,
contains only one component. When minimum gaussianity is observed, the corresponding
signal can be regarded as an estimate of one component of the original signals. Hence, it is
very critical to find a metric to measure the gaussianity of random variables.

Measuring Gaussianity

The classical method for measuring gaussianity is based on using kurtosis, the fourth

33

order moment [25]. Kurtosis is widely used to measure gaussianity because of its
computational simplicity. Another method, which uses negentropy, is based on information

theory and differential entropy. These two methods are described next.
For a random variable s, kurtosis is defined as
Kurtosis(§) = E{§4 I — 3(E{§2 })2 (3.9)
where E is expectation operator.
Since s is assumed to be zero-mean, unit-variance random variable, kurtosis is
simplified as
Kurtosis(§) = E{§4 } - 3 (3.10)
This is the normalized version of fourth moment of a random variable. For a Gaussian

random variable, the fourth moment is equal to 3(E{§}2)2 . Therefore, kurtosis is zero for

Gaussian random variable and is positive or negative for non—Gaussian variables. Positive
kurtosis is referred to as supergaussian property for which the distribution is more similar to
that of a Gaussian distribution. A spike at the rrriddle and approaching zero at the tails
characterize positive kurtosis. Negative kurtosis refers to subgaussian property and is
characterized by a relatively flat distribution.

To illustrate how the original signal is reconstructed using kurtosis, we need to know

the properties of kurtosis. When x, and x2 are independent to each other, kurtosis has

following properties.

Kurtosis(x1 + x2 ) = Kurtosis(x1)+ Kurtosis(x2) (3 ll)
Kurtosis( ax 1 ) = a4Kurtosis(x1 ) .

where a isascalar.

34

Let us assume that E{(wa)2} =1 and let the two-dimensional mixing model be
given by x = As where random variable x is the measured data, 8 is the original source

signals and A is a 2 by 2 mixing matrix. Define another vector 2 = ATw. It can be

shown that

E{(wa)2}=E{wTAATw}=||z“2=1 where E{ssT}=I (3.12)

. 2 . . . .
Under the constraint “2“ =1, kurtosrs functron has a number of local rrnnrma and

maxima. To make the argument clearer, let us assume that one component has negative
kurtosis and the other component has positive kurtosis.

Let § bethe weighted vector of x.

=wa=wTAs=sz= zls, +2:st (3.13)

M)

From properties of kurtosis, we have

Kurtosis(§) = Kurtosis(z 131 )+ Kurtosis(z2s2 ) (3 14
- )
= 2,4Kurtosis(s, )+ z24Kurtosis(s2 )

Since s and s are assumed to have unit-variance as a result of whitening in the

preprocessing stage, we have unit variance of s and s.
E{§2} = z,2E{512}+222E{522}= 2,2 + z22 =1 (3.15)

It is clear that (21,22) is located on the unit circle that satisfies above equation. Now,

the problem is to find one of the maximum values of 2,4Kurtosis(s, )+ z 24Kurtosis(s 2 )j

with a constraint of z ,2 + 222 = l . Due to the assumption that kurtosis of s, and 32 have

different sign, it can be easily shown that the maximum value occurs when one element is

zero and the other is +1 or —1. Therefore, there are 4 maximum values of

35

214Kurtosis(sl )+ z24Kurtosis(sz)| with a constraint of Z12 + 222 :1.

In practice, the strongest direction of initial weighting vector is computed (increasing
direction for positive kurtosis, decreasing direction for negative kurtosis) based on the
available samples of x. Then, the gradient method is applied to find a new weighting
vector w. Because of this computational simplicity, kurtosis has been widely used for
estimating ICA model. However, the disadvantage in using kurtosis as a metric of
measuring gaussianity is that the accuracy of the estimate of kurtosis depends on the
availability of adequate data particularly in the tail of distribution. This implies that kurtosis
is not a robust method for measuring gaussianity. To compensate for this weakness, an
alternate technique based on negentropy is used and is described next.

Negentropy is an alternate method to measure gaussianity of random variables that is
based on information theory and provides a quantitative measurement of differential entropy.
Entropy is an essential concept of information theory and it measures the amount of
information contained in an observation. Entropy increases as the random variable becomes
more unpredictable and unstructured [26].

Entropy, H , of a discrete random variable is defined as

H(§discrete) : _Z P(gdiscrete : k1 )log P(gdiscrete _—._ kt) (3'16)
1

where sdmm is a discrete random variable and ki is a possible value of sdiscm with
probability density function (PDF), P(§ discrete) .
For a continuous random variable, it is defined as
H(§ continuous) : —.i f(§ continuous )log f6 continuous )d§continuous (3'17)

where s is the continuous random variable with PDF f(§commuous).

confinuous

36

From the fundamental proposition of information theory, a Gaussian random variable
has larger entropy than any other random variable with the same variance. This implies that
entropy can be used as a quantitative measure of gaussianity. To obtain a measure of
gaussianity, the modified definition of differential entropy is used. It is called as negentropy

and is defined as
J(§) = H(§Gauss ) — H(§) (3.18)
where scam is a Gaussian random variable which has the same covariance matrix with 11.

By definition, negentropy is always positive and zero if and only if s is Gaussian
random variable. When the statistical property is considered, negentropy is a well-defined
metric to measure non-gaussianity [26]. However, in general, negentropy has the
disadvantage that it is very difficult to compute. Therefore, a modified approximation of
negentropy that combines the strengths of kurtosis and negentropy is often used to measure
the gaussianity.

A commonly used approximation of negentropy is given by

A 1 .3 2 l . A 2
J s z—E s +—Kurtosrs s 3.19
() 12 { } 48 () ( )

where s is zero mean, unit variance random variable. The above approximation includes
the kurtosis function, which leads to issues of robustness of the metric. Also, it has been
argued that cumulant-based approximations of negentropy are inaccurate in many cases
[27,28]. Therefore, a new approximation method based on maximum-entropy principles is

often used and is expressed as

J(§) z Zki [Etci (§)} —E{Gi (§GMSS)]2 (3.20)

37

where ki is a positive constant and scans, is a standardized Gaussian random variable
with zero—mean, the same variance as s and G (s) is a function that does not grow faster
than quadratically as a function of Isl [28]. According to equation (3.20), negentropy is

always positive and zero when s is Gaussian. The choice of function G1 is very

important. In the practical choice of the functions Gi(s), the following criteria must be

emphasized: First, E{Gi (1‘)} should be estimated without statistical difficulty. Second, the

maximum entropy method assumes that the density function f (s) is integrable, because
the expectation is obtained as

E{Gi (5)} = —jf(s)Gi(s)ds (3.21)

Therefore, Gi must not grow faster than quadratic function, because a function

growing faster might lead to non-integrability [28]. Usually the following choices show
robust and reliable results and are called contrast functions [10, 11].

G](u) = ilogcosh(alu)
3' (where 1 s a, s 2) (3.22)

2
(4

G2 (14) = —e 7
Approximations with the above G function have combined the advantages of two
metrics, kurtosis and negentropy, used to measure gaussianity. It contains all the statistical
information and is much easier to compute than negentropy. A practical algorithm, which is
based on this approximation, is developed and is known as the Fast-ICA algorithm.
F ast-ICA Algorithm
Fast-ICA algorithm is one of the most powerful algorithms for estimating individual

non-gaussian components in a mixture of source signals. The procedure of Fast-ICA

38

algorithm starts with choosing an initial weight vector w , which is updated by the learning
rule. With the chosen weight vector, gaussianity is measured by J(§) given in the equation

(3.20). The purpose of Fast-ICA algorithm is to find a vector w such that its projection on

T1:, minimizes gaussianity. There is a constraint that the variance of wa is

x , s = w
unity. When it is preprocessed, the norm of weight vector can be calculated from this

constraint. Recall that the norm of preprocessed data is unity.
r 2 T r T
“w x“ =w xx w=w w=1 (3.23)
The Fast-ICA algorithm is based on a fixed points iteration scheme for finding a weight

vector that minimizes the gaussianity of wa. The algorithm is derived from derivatives of
the non—quadratic function G in (3.21). If the two functions in (3.21) are chosen, the

derivatives are given by

 

 

3:1“) = gl (u) = tanh(alu)
u
"2 (where 1 S a1 S 2) (3.24)
3G2(u) = g2(u) = “6‘?
Bu

The details of Fast-ICA algorithm is summarized as follows:
1. Preprocessing: The data is first preprocessed with centering and whitening.

Centering subtracts the mean to obtain a zero-mean signal and whitening produces

components that are uncorrelated with unit variance, giving ||§|| = "wan =1 when

"w“ = l . Whitening also transforms the mixing matrix A to be orthogonal.

2. Select an initial, random weight vector w with the constraint "w“ = l.

39

3. From (3.20), the negentropy is maximized when E{G(§)} is maximum or

minimum since E{G(§Gauss )} is a constant. Thus, the optimization process seeks to
find optima of E{G(§)} under the constraint of “éllzllexllzl. It can be

expressed as

F(w) = E{G(wa)}— 8(j|w|| — 1) (3.25)

4. Using the Kuhn-Tucker theorem [29], the solution is obtained as a solution of

following equation.
F'(w) = E{xg(wa)}— 13w,- 2 0 (3.26)
where B is a constant called a Lagrange multiplier.

5. The updated weight vector w+ is obtained using Newton iteration [30] and is given

by

w+ = E{xg(wT )}- E{g’(wa)}w (3.27)

6. Normalize the updated vector

+

w"""' = w+/ w (3.28)

 

 

 

 

7. If w""”’ 2 wow, §= wa is an estimate of the source. Go to Step 2 to compute
the next source signal. If not, go to Step 5.

When solving for multiple sources, wi,i =1,2,...,n are obtained from the above

. T T T
process. The estimated sources w] x,w2 x,...,w,, x are then uncorrelated after each

component is obtained. Decorrelation is achieved by Gram—Schmidt decorrelation [31].

From this, the independent component is calculated one by one. When one obtains a new

40

weighting vector w HT, the projection w p+1ij, j=1,2,..., p on each of the

p

previously estimated vectors is subtracted from w p +1 , which is then normalized. Therefore,

the separated independent components are unique up to a multiplicative sign and scale

factor. The schematic overview of the process is illustrated in Figure 3.2.

41

 

I Choose an initial weight vector : W I

 

 

 

 

% Update w by learning rule : WWW I

 

 

 

I normalize w I

 

 
  

wnew has the same direction
with w°1cl

 

    

 

Decori‘elate Wne‘” from obtained w I

 

 

 

 

normalize wnew

 

 

 

 

WWW is a colurrrrr of the inverse matrix (W) I

 

 

Figure 3.2 Flowchart of Fast—ICA Algorithm

42

Noisy ICA Model (Residue Theory)

The ICA model cab be used to extract source signals successfully when noise is
usually present in the measurement. In general, methods taking noise explicitly into
account assume that the noise has Gaussian distribution. As discussed before, the Fast ICA
algorithm measures gaussianity of the mixture and minirrrizes it to extract a source signal.
Therefore, Gaussian source signals cannot be extracted using independent component
analysis.

However, when only one Gaussian noise signal is present in the mixture, ICA can find

noise using residual theory. Residual theory states that a function f (z) has a simple pole at
20 , we can write

a0

 

f (z) = + g(z) (3.29)

Z—ZO
where g(z) is analytic at 20, so that (z—z0)g(z)—>0 as z——> 20. It follows that the

residue is given by

“0 = lim (Z — Zo)f(Z) (3.30)

z—>z0
It is assumed that fj(z) is analytic in C except at 2].. For j=1,2,---,k, let
f j(z) denote the principal part of f (2) at 2].. Using this residue theory, Gaussian

noise g(z) can be found as
k
g(z) = f(z)- Zlfjiz) (331)
j:

where g(z) is analytic in C, so that Ig(z)=0.
C

43

3.4. Independent Component Analysis Based on Mutual Information

Mutual information is another approach for ICA which measures the mutual dependence
of the components of a transformed vector, 5 and then optimizes the unmixing matrix, W
to minimize the measured dependence of each components. There are many methods to
measure the mutual dependence and Shannon’s mutual information is considered to be one
of the popular choices.

Shannon‘s mutual information between 11 random variables, ii,i=1,2,~-,n , is

defined as
I(§) = item—He) (3.32)
i=1

The mutual information, 1(8), measures the amount of shared information in the

components of 8. It can also be considered as a measure of the independence of random
variables. Naturally, it is also non-negative and is zero if and only if the components of 8
are mutually independent with respect to each other. Hence, the primary purpose is to find a
transform that minimizes the mutual information between measurements of original
components [32]. One of the most important properties of mutual information is that it is
invariant under continuous and monotonic transformations. Therefore, mutual information
does not change unless the components are mixed. For invertible linear transformation

§ = Wx, mutual information is defined as:

1(§,,§2,...,§,,) 2 211(3) — H(x) —logIdet WI (3.33)

If 5,- is uncorrelated and has unit variance, then the covariance matrix of 8 is an

identity matrix. Hence,

44

E ԤT}= WEIXXTIVVT =1 (3.34)
det EI§§T I: det WEIxxT IWT = det Wdet EIxxT Idet WT =1 (3.35)

which implies that det W always has a constant value. In this case, the equation (3.33)

reduces to

I(§]$§29”'9§n):C—ZJ(§1°) (3.36)

From equation (3.36), we see that negentropy is in inverse proportion to mutual
information. Mutual information and negentropy differs by a constant and sign difference
for variables with unit variance. When the source signals are uncorrelated, minimizing
mutual information of the estimated signals is equal to maximizing the gaussianity.

Therefore, maximizing negentropy is equivalent to minimizing the mutual information of

the separated components of §.

The use of mutual information is also motivated by the Kullback_Leibler divergence

which is defined as:

5(fhf2): If] (9108115245 (337)
his)

where f1 and f2 are probability density functions (PDF). The Kullback_Leibler
divergence can be considered as a distance between two probability functions. If 8,. in
(3.32) are independent to each other, the joint probability function could be factorized as
f (8182 ...8, ) = f1 (81 ) f2 (82) f" (8") . Thus, the independence of 5,- can be measured by
the Kullback-Leibler divergence between the real PDF and the factorized density function
_/’E = fl(§l)f2(§2)... f" (8") , where f,-(-) are the marginal density function of 3“,- .
Theoretically, the measured independence of i,- is equal to the mutual information of 5} .

45

Kullback—Leibler divergence also shows that minimizing mutual information is closely
related with maximizing likelihood. The likelihood is represented as a Kullback—Leibler
distance between the observed density and a factorized density. The problem with mutual
information is that the estimate of PDF is needed to obtain the entropy that restricts the use
of mutual information in ICA implementation. This problem can be solved by using an
approximation of mutual information based on polynomial density expansion [32], which is
related to the Taylor expansion that gives an approximation of PDF. For example, if § is a
random variable having zero mean and unit variance, the first term of the Edgeworth

expansion gives [33]:
ftri=¢trI1+ks<§>hs<%+k4<§>h4<%4+...I 1....)

where (o is the density function of a standardized Gaussian random variable, the k,- (8) is

the cumulants of § (zero-mean, unit variance) and [1,-(1) is a polynomial function. Using

the above equations, the following approximation is obtained for mutual information.
. 1 m ,. . . . .
1(s) z C + 48 2 [41.3 (s,- )2 + k4(s,. )2 + 71.46,.)4 — 6k3 (s, )2k4(s,. )I (3.39)
i=1

where C is constant and 8,- are uncorrelated to each other.

46

4. APPLICATION OF ICA TO NDT

4.1. Feature extraction algorithm

Multi-frequency eddy current testing (MFECT) is one of the most widely used methods
for inspecting steam generator (SG) tubes in nuclear power plant. However, signal
processing is often required to extract the relevant information from measurements because
there are many factors including noise that affect the eddy current measurement. This thesis
addresses the problem of separating the defect signal that is distorted by contributions from
benign sources such as tube support plate (TSP) in steam generator inspection. Novel signal
processing techniques have been proposed for mitigating difficulties encountered in solving
the underlying inverse problem. Feature extraction techniques for eddy current signals are
used to identify and eliminate irrelevant parts of the data generated by edge effects, probe
lift-off noise and white noise. It is possible to look at the noise as derived signal from an
additional source and develop a method, whose rationale is to separate the desired signal

from the other interfering signals. Figure 4.1 shows the basic scheme of using ICA to

 

 

 

MFECT data analysis.
MFECT at f1 —+ ——+Defect
MFECT at f2——~ ICA HS?
MFECT at f3 ——b fiJ‘Nojge

 

 

 

Figure 4.1 Application of ICA to NDT

47

In order to apply ICA to steam generator tube inspection problem, there is one
essential requirement in that the number of ECT measurements should be equal to or
more than the number of expected source components. In this study, the measured signal
is assumed to be made up of 3 components namely, defect, TSP and noise. Since there are
3 sources to be estimated, it is necessary to have at least 3 different measurements. In this
application, Multi—frequency Eddy Current Technique (MFECT) is used to obtain the
measurements where the probe is excited at 3 different frequencies. Due to skin effect
phenomenon, signal at low frequency carries information about defects that are at a larger
depth than the high frequency signal that carries only surface information due to the shallow
skin depth. The details of this phenomenon were described in chapter 2. Hence, each set of
ECT data includes different combinations of the source signals depending on the excitation
frequency that enables ICA to separate out the source signals. Figure 4.2 illustrates the

schemes proposed method in this thesis.

48

 

MFEC Test Data

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 
 

i
Segmentation
*
RST
9 4
Real Imaginary
t L
ICA ICA

 

Combine Signal

 
 

 

Support Plate

 

 

 

Defect

 

 

Noise

 

 

 

Figure 4.2 Schematic diagram of ICA analysis in NDE

In Figure 4.2, the overall scheme consists of 4 stages. Prior to using MFECT data in the

ICA algorithm, they are tuned to a common frequency by means of the rotating, scaling

and translation (RST) transform described in section 2.5.2.

1. Segmentation:

In steam generator tube inspection, the probe is pulled through the entire tube
resulting in a large amount of data that typically includes multiple TSPs and defects. The
ECT signals are calibrated and preprocessed to extract segments of signals with single
indication typically at the TSP locations. Consequently, each of these segments can
contain contribution from TSP alone or TSP and defects or TSP, defect and noise. These

signals are then applied to the ICA algorithm. The MFECT data collected are segmented

49

 

 

 

to include contribution of a defect and tube support plate (TSP).
2. Signal Transformation:

Rotating, scaling and translation (RST) transformation is applied to the segmented
MFECT data. RST transforms a low frequency eddy current data to corresponding high
frequency signal. The RST parameters, such as rotating angles, scales and translation
parameters, are obtained from the calibrated standard ECT data. The result of RST

transformation is given below. Figure 4.3 shows the result of RST transform applied to

ECT data at 200klb. The transformed signal is compared to the signal predicted when the

excitation frequency is 400kllz.
3. ICA algorithm:

After MFECT data sets are transformed using RST, ICA is applied separately to the
real and the imaginary parts of MFECT data. As a result of ICA, the real and
imaginary parts of component source signals are estimated. The real part and the
imaginary part of the same component — defect, TSP or noise — can then be combined
to estimate the impedance trajectory of the defect and TSP.

4. Combine the extracted signals that are a similar distribution to produce the impedance

plane trajectory of the estimated signal.

50

Impedance of RST data

 

 

 

 

 

500 I I I I’ r
CV41“ 5 E I I -+- imp200 after RST
0% 5 i i -*> ""9400
400 ~\‘\\& ..... ............. -*~ 'mpzm p
5\ 4.3.: ' . 5 1* E
5 \e‘+.‘~.\ ‘ 5
\~. 4

----------- \‘X ---------- a" ------------- .... ....... .

 

 

 

 

: EA 2 2 i
'200“ """" («3’ +H\G\\\s """""" -
5 it“ E “‘34 E
: : . \\0‘~.\'
-400 ~------------; ------------- ;. ------------- ------ 1.3.8.33. --------- -
: : : g \‘~‘;;%:G\
- - = 5 “4‘9
-500 1 1 1 1 1
-600 400 -200 0 200 400 600

Figure 4.3 Impedance trajectory plots of ECT data after RST

4.2. Finite Element Method (FEM)

For evaluating the performance of ICA algorithm, simulated MFECT data were
obtained from a theoretical model that was used to generate ECT signal with TSP and
defect parameters. Finite element method (FEM) has been widely used for solving the
governing partial differential equations underlying physical process. Analytical,
numerical and hybrid methods have been in existence for solving the governing equations.
Finite element method is based on the principles of variational calculus. The solution to

the governing differential equation involves the incorporation of the equation in an

51

integral form using energy functional. The energy functional, which represents the energy
of the system whose stationary value is a minimum, is minimized resulting in the solution
to the governing equation.

The geometry consists of an infinitely long tube within which moves a probe as seen
in Figure 4.4. The axisymmetric nature of the current source and the vector potential
allows the reduction of the problem to the two-dimensional plane and a coarse finite
element mesh spanning the region is given in Figure 4.5. The finite element method
consists of following steps.

1. The region of interest is discretized with a suitable mesh consisting of a number of

elements connected at the common nodal points as shown in Figure 4.5.

2. The node and the elements of the different materials in the regions are identified and
numbered.

3. An interpolating function, which approximates the continuous field over each element
in terms of the nodal point value is defined in such a way that the field is continuous
across the element boundaries. The interpolating function can be linear or nonlinear
and depends on the variations of the field in the test geometry.

4. Minimization of the energy functional with respect to the unknown nodal point values
results in a matrix equation.

5. The solution of the matrix equation yields the field values in the region of interest.

52

Defect

 

Figure 4.4 Geometry for eddy current NDT [34]

Eddy current
Probe

Tube support ..........
Plate """ >

Defect

 

Figure 4.5 A finite element mesh for eddy current geometry

53

The finite element formulation in terms of magnetic vector potential A for the two
dimensional axisymmetric eddy current phenomena was developed by Lord and Palanisamy
[35]. Once the magnetic vector potential values at all the nodes in the mesh region are
determined, the probe impedance that is our parameter of interest can be computed. TWO
different approaches commonly used to estimate this value, are the direct and energy
methods and are described in more detail next.

Direct Method

The impedance of a single turn coil of radius r carrying an alternating current of Is

amperes is given by

_Y_ (4.1)

l:
I

where V is the RMS phasor voltage induced in the coil, expressed in terms of the

electric field intensity E as

V=—pE-E (4»
From equation (2.13), we have
E=_§§_v¢ (4a
3.

Assuming that the field varies harmonically with frequency w given by the

equafion

A = Aoe'l'w’ (4.4)
we have

E = -—ij — V¢ (4.5)

54

Since the induced voltage is independent of the gradient of the scalar potential V¢,
substituting equation (4.5) into (4.2)

V==mpAw me)

From equation (4.1), the impedance of the coil is

2:43pAw (4n

S

which for a single turn coil of radius r is

 

(4.8)

The real and imaginary part of which can be interpreted as the resistive and reactive
components of the impedance.

The coil impedance however is calculated in an approximate manner using the finite
element method. Consider the cross-section of the coil, discretized by triangular elements.
If the dimensions of the elements are small than the vector magnetic potential of all the

turns covered by the element i can be approximated by the centroidal value Aci and

similarly the radius of all the turns in the element can be approximated by the centroidal

value rci. From equation (4.8), the impedance for each turn within the element i is
given by

Z : jZflwrciAci

13 (4.9)

If N3 is the total number of turns in the coil cross-section, Nt tums/mz, the turn

density of the element and ai the area of the element, the total impedance of all the

turns in the element 1' is given by

55

27M ~A . ~N
Z=j rClI Clal I (4.10)
S

 

If the number of elements in the cross-section of the coil is NC , the total impedance

2 -A - N.
z: jfluzairciAc, (4.11)
l

I s i:

The coil impedance in a two dimensional or axisymmetric problem can be computed
using the method explained here. This method assumes that the magnetic vector potential
values are constant along the source, i.e. along the circumferential direction of the coil.
However, this is not true in the case of 3-dimensional problems. An alternate method to
compute the impedance based on the calculation of the stored and dissipated energies is
explained next.

Energy Method

The impedance of a coil can be calculated from the energy of the system since the

inductance and the resistance are associated with the stored and dissipated energies, in

the system, respectively. The stored energy E is given by

E: IVE-Haw (4.12)

From equation (4.12) assuming constant reluctivity in each direction the energy

stored in a finite element of volume Vi can be written in terms of B as
1
Bi =3(va,(2 +vyBy2+szZZIVi (4.13)

where vx , v and vZ are the reluctivities in the corresponding directions.

Y

From the relation between the magnetic flux density and vector magnetic potential

56

 

g

B we have

 

 

B. _8A_z_iA_
By 82
8A 8A.
B = X _ l 4.14
y dz 8x ( )
Bx : dAy _an
8x 8y

Substituting this in equation (4.13) and adding all the elements in the mesh region,

the total stored energy in the system is given by

N . 8A. 2 . . 2 3A. . 2
EZZ vx[aAZI_ 1’1] +vy[an1_aA21)+vz[ Yl_an1] i

 

 

 

 

i=1 dyi dzi 821 3x1 Bxi dyi
(4.15)
The inductance of the coil can then be calculated using
LI 2
E: 25 or L=12—I;: (4.16)

S
where IS is the current in the source coil.
The resistance of the coil is associated with the dissipated energy in the system. The

dissipated energy in a finite element of volume Vi is given by

2
V- l
P. = ——'| ed”! (4.17)
0'
where chdy is the eddy current density and was derived in chapter 2 as
Jeddy : _jwaAci (4.18)

where Ad is the centroidal magnetic vector potential value for the element i.

Substituting equation (4.18) into (4.17) and summing all the elements in the mesh

57

region, the total dissipated energy is given by

N N 2 2
P: 2P1= [View IAciI (4.19)
i=1 i=1
Using P = 12R , the resistance of the coil is given by
P
R = 1—2 (4.20)
and the coil impedance is
Z=R+ij (4.21)

The above formulation is applicable to two dimensional and axisymmetric problems.

By definition the total flux linkages in a coil with N5 turns is equal to NS¢ where (I)

is the flux linkages due to a coil with 1 turn. Multiplying the corresponding values of the
magnetic flux density and the vector potential value by N, we have the impedance of the

coil with N3 turns to be

zN = N522 (4.22)
where Z is the impedance of the coil with 1 turn.
Flux Density and Eddy Current Density

Once the vector magnetic potentials are computed, the formulations of the flux and

eddy current densities for the axisymmetric two-dimensional case are obtained as

B=VxA

. (4.23)
Jeddy : ’JWUA

where B is the flux density and J eddy is the eddy current density.

58

5. SIIVIULATION RESULTS AND DISCUSSIONS

The proposed ICA algorithm was first applied to synthetic data obtained from the finite
element model (FEM) and then on experimental data collected from the steam generator
tubes in a nuclear power plant. The function of the source signals separation from the
observed multi-frequency eddy current (MFEC) data is focused as the essential feature of
ICA. However, application of ICA for noise filtering is also tested. Fast-ICA algorithm and
preprocessing are employed to extract the original signals from MFEC data. For the FEM
simulations, only defects, support plates and noise are considered as source components.
For 3 distinct sources, we need at least three sets of MFEC data at different excitation

frequencies. Hence, FEM simulation of MFEC inspection process was conducted at three
excitation frequencies, namely, 35kllz, 200kllz and 400kllz. Due to the property of ICA, the
sign or scale is not guaranteed for the extracted source signal. Hence, the sign or scale has
little meaning at the result. Instead, the distribution of signal gives information about the
extracted source signals.

Two sets of results were obtained from implementation of ICA on FEM data. The first
set of results considered a defect in a support plate region at 35kllz and 400lde and the
second set of results evaluates the performance of ICA in filtering noise as well as
separating source signals. For evaluating noise filtering, random white noise (2 dB SNR
with 50 dB of signal power) is added to ECT data. The noise component is assumed to be

the same at each excitation frequency. ECT data at 35kllz, 200kllz and 400klt are used as

input to ICA and three source signals due to defect, support plate and noise, are extracted as

output. In addition, the signal to noise ratio of the extracted flaw signal is seen to be

59

significantly higher, indicating that ICA can be used to denoise data. These results are
shown in section 5.1.

The ICA algorithm was also applied to MFEC field data collected at excitation

frequencies of 35kllz, 200kll2, 400kllz and 600le. The MFEC field data and the corresponding

results are shown in section 5.2.

5.1. Simulation of FEM Data

The finite element model allows us to obtain signal under controlled conditions of defect
geometry, excitation frequency and material properties. The geometry including support
plate and defect is shown in Figure 4.5. Table 1 shows the value of the parameters used in
the model. The real and imaginary components of the probe impedance for the geometry in

Figure 4.5 at frequencies 35kll2, 200kllZ and 400kllz are shown in Figure 5.1, 5.2 and 5.3,

 

 

 

 

 

respectively.
Electiiical Relative Current Densit
Category Conductivity Permeability (Almz) y
(ohm-meter)’l (it)
Tube 1x106 1 —
Support Platel 10x 106 10x106 -
Coil 0 1 106

 

 

 

 

Table 5.1 Parameters for Finite Element Model

60

 

 

4n \1,‘ r T r T i r
'\\|_ i i i i l .
\(K~ ; I I 1 I l
. ‘-. : : : : : :
30 )..........|.>~. -\. ..... ... ............................................ 1. ........ 1.4
2 \ “-1. z r : : :
: \ : : : : : :
20 . \ . \\ . . . . .
IN; ........ ‘ ...: ...... 1.... ...... ..;........;. ...-..- ........ .4
: x r... : : 1 :
10 : 1 : ~ : .
r- ....... , .......... ...-\.; "“1"; .......... .... ....... ..,
‘<,_ \; : :
U ‘5: It) 2 . .
.. ......................... ;... T‘: ’ - ...; ' , ...
' 1: \‘-~=:a~. : :
- \-. : : :
'10 \t ‘3. : '
p. ................... i ................... I". ....... 3.3“- ..... l'- ....... - ......... .1
I \ - V, 1 1
1 \ ' \xi l
-20 1 “ '
. ‘ _'\‘
p- y l ..g. 1 ‘1" . .3.. \1..-‘ .4
. \ . 4,‘
3 \ I .
-30 \I ‘1‘:
r q \,i
.. ....\.. 11‘ ....... a
1 1 1 1 ..1- _- _ 1...- -_ -1-.- .J

 

 

(a)

 

 

 

 

 

 

 

 

 

I I I I I I l 5 TI
1 1 i 1 1 i i i 1
AU 50 3:] III] 120 140 150 1m 2]] 220

 

 

 

 

 

 

l l l l I. I l l
i 1 i i i 1 1
40 50 33 III] 120 140 160

 

(b)

Figure 5.1 ECT data from FEM at frequency of 35 kllz (a) Impedance trajectory (b) Real
component of eddy current signal (upper) and imaginary component of eddy current signal

(lower)

 

 

 

 

(a)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(b)

Figure 5.2 ECT data from FEM with excitation frequency of 200 ldlz (a) Impedance trajectory
(b) Real component of eddy current signal (upper) and imaginary component of eddy current

signal (lower)

61

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(a) (b)

Figure 5.3 ECT data from FEM with excitation frequency of 400 kill (a) Impedance trajectory
(b) Real component of eddy current signal (upper) and imaginary component of eddy current
signal (lower)

5.1.1. Source Signals Separation (Defect and TSP)

The main objective of ICA is to extract defect signal from multi-frequency eddy current
(MFEC) data in the presence of signals from external structures such as TSP, tube sheets
(TS) and conductive deposits in addition to the flaws in the tube. The purpose of this
simulation is to evaluate the ability of ICA to separate out the defect information from the
irrelevant signals from other components using ECT data generated by finite element

model (FEM), which is noise-free. MFEC data at 35kllz and 400kl’c are shown in Figure

5.1 and 5.3, respectively. The real part (horizontal component) and the imaginary part
(vertical component) of ECT data are considered separately, processed by ICA and

reassembled to construct impedance trajectories of source components. The real part and

imaginary part of MFEC signals at 35kllz and 400kllz are shown in Figure 5.5 where the two

62

indications at the end correspond to support plates and the indication in the middle is
generated from the defect. Reconstructed defect signal and support plate signal are shown in
(c) and (d) along with the corresponding impedance plot. The impedance trajectories of the
defect and the support plate signal extracted from mixed MFEC data are shown in (c) and
(d) in Figure 5.4.

In the mixing algorithm described in chapter 2, the “mixed” signal still contains some
residual signal at the origin indicating that the effect of support plate is not completely
removed. However, the relevant source information due to the defect extracted from the
mixed signal using ICA allows distinguishing the defect information without any
contribution from the TSP. The real and imaginary components of data are shown in Figure

5.5.

63

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

x102

Figure 5.4 Impedance trajectories using FEM at (a) 35 kllz (b) 400 kHz (c) separated defect
signal ((1) separated support plate signal

64

 

 

 

 

 

   

 

 

 

 

 

 

501111111.
0204060831m1201401601m2m220

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

     

 

 

 

 

 

 

 

. l I l l i l l i

50 ..................................................... -

100 i i r 1 1 1 r r 1 1

20 40 50 80 1CD 120 140 160 180 ID 220
(a)
5m
4m
2m

0.
cm 0
cm ' 0 : 0 C if
an i 1 1 1 i 1 i 4

0 m 40 a] a] 140 1m 1m m 220 0 2D 40 50 3] 1m 120 I40 160 180 2111 2m
6C0 I I I I I r I I l i i i f T. r r r 7
4m»- ....................... .4- 5. ................... : ............. ...... ........................ _
m —....... .........+ ----------------------- L-----¢------c---o--E-o---—i ' i g I g .

o 1 ,- Saw-i 3 I
.ms.----.:.-..-..-----+-----.--u...---...-“...-.na-n-non ”i ..... —- E I E E 5
we -------------------------- f ------ ------ ------ i ------ : ----- -
In, 1 1 1 1 ' 1 r 1 1 1 i 1 1 1 i 1

0 20 40 50 a] 40 160 180 200 220 80 III] 13] 140 160 1m 10 220

(b) (d)

Figure 5.5 Real component (Upper) and imaginary component (Lower) of (a) FEM—ECT
signals at 35 kHz (b) FEM-ECT signal at 400 kHz (c) extracted defect signal ((1) extracted
support plate signal

65

 

5.1.2. Source Signals Separation (Defect, TSP and Noise)

The purpose of this simulation is to evaluate the ability of ICA to function as a noise
filter as well as source separator. ECT data is often corrupted by the presence of noise or
interference signals depending on inspection conditions. Sometimes, defect signal is buried
in noise signal which makes it difficult to extract meaningful defect information from the
measured ECT data. This section describes the use of ICA to remove noise. Usually noise is
modeled by a Gaussian distribution which violates one of conditions for ICA. However,
when there is only one component that has Gaussian distribution, residue theory [36] can be
employed to filter out Gaussian noise from observed ECT signals.

The ECT signal generated by finite element model (FEM) at excitation frequencies of
35 kHz , 200 kHz and 400kHz are corrupted by adding random white noise signal (2 dB
SNR with 50 dB of signal power) to each data set. The noisy signal is input to ICA which
produces 3 separated source signals. These results are shown in Figure 5.6. Figure 5.7
shows noisy version of the horizontal (real) and vertical (imaginary) components of FEM-
ECT data in (a) and (b), respectively. The resulting source signals due to — defect, support
plate and noise — are extracted and shown as impedance trajectories in Figure 5.6. The result
shows that source components — defect, support plate and noise are successfully extracted

from noisy ECT measurements.

66

 

(C) (1)

Figure 5.6 Source signal separation with noise filtering (simulated trajectories): (a) Noisy ECT
data at 35 kHz. (19) Noisy ECT data at 200 kHz, (c) Noisy ECT data at 400 kilz, (d) Separated Defect
Signal, (e) Separated Support Plate Signal, (1) Separated Noise

67

 

 

 

 

 

 

 

 

 

 

Emlmlmufllfifllwmm

IEUIBOZIDZQU

140

 

 

 

 

 

 

 

 

 

 

40

1

80101120

50

 

 

180200220

103 120 140 160

20406080

 

 

 

60801001201401601832C02200

40

 

 

 

(c)

(a)

 

 

 

 

 

 

 

 

608010312014015018010220

 

 

 

 

 

 

 

 

220

10112014016013] 2m

80

204050

 

 

 

 

 

 

 

 

1802113220

1CD 120 140 160

020406080

180200220

50 80101120140160

40

(d)

(b)

Figure 5.7 Source signal separation with noise filtering (simulated data): (a) noisy horizontal ECT

(b) noisy vertical ECT data at 35 kl’a (top),

200 kllz(middle) and 400 kilz(bottom), (c) separated signals (horizontal component), ((1) separated

data at 35 kill (mp). 200 kilz(middle) and 400 kilz(bottom),
signal (vertical component)

68

5.2. Performance on Experimental Data

While the synthetic ECT data is obtained from finite element model (FEM), the
experimental data are collected during inspection of steam generating tubes. In fact, a variety
of operational parameters and material conditions such as lift—off, local variations in
permeability, conductivity and material composition, surface roughness, etc. give rise to a
noisy signal, which makes the detection of the true defect signals a formidable task.

The impedance signals from the EC probe are measured at 4 excitation frequencies —
35kllz, 200kllz, 400le and 600kli2. The maximum number of components that can be
estimated is equal to the number of available data sets. Hence, at most 4 source components
can be estimated. In this test, 3 source components are assumed, namely defect, support
plate and noise. In the preprocessing stage, rotation, scaling and translation (RST)

parameters are obtained from the calibration standard as described in 2.5.2 and 4.1.

5.2. 1. Extraction of the Defect Signal from Mixture of Defect and Noise

ECT field data obtained from Nuclear Power Plant (NPP) SG tube inspector was used in
this test. Figure 5.8 (a) presents a set of 4 horizontal channel EC signals obtained at different
excitation frequencies — 35kiiz, 200kilz, 400kiiz and 600kilz. Preprocessing using RST
transformation is first applied to the signals to transform all signals to match the
corresponding signal at 600kilz. The parameters for RST are obtained from the standard

calibration. The RST processed signals are applied as input to the ICA algorithm. Figure 5.8
(b) shows 4 independent components extracted as a result of ICA. The first separated

horizontal signal has a big pulse at the corresponding defect location (77~85) in the top

69

figure as indicated. The signal at bottom of Figure 5.8 (b) is regarded as horizontal
component of noise.

The result of independent component analysis of the corresponding vertical channel
signals is shown in Figure 5.8 (b) and (d). In the separated components in Figure 5.8 (d), the
defect signal is extracted as the second independent component and the defect location is
indicated. Extracted horizontal and vertical channel noise components are the last signals in
the bottom of Figure 5.8 (b) and ((1), respectively. Figure 5.9 shows the impedance plane

trajectories of defect, TSP signal at 35kHz, 200kil2, 400kil2 and 600kilz. Figure 5.10 shows the

impedance plane trajectories of defect and TSP extracted by ICA algorithm. The result
indicates that it is possible to extract defect information from MFECT data by using ICA

algorithm.

70

 

 

20134050

10

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

.1 m m em
.i... ... F----
f fi\._u
..-.-4..-1m milii
-\\" M. ”
11M---Mr-r-lm m. r J---l
T -..Win.-lmw v. wT-M.lmiiirl

1W 01m wI-Vuu.uwurinl
..... WIIIM-fl 1 “fl .-M-am----l
.- -.-m-.. ... .-.-m
.- 4-1--.... m. 1
m a m m
.b 0 U I».
«is, 0 .11 al 0 4|
0 0 0, 0

339310]

203040506070

10

 

marmsosuroaosorm

10

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

m m
m w
m a
.m- .nr.
m m
)
Au. ..w m
w w
w m
m m
m. m.
00 U
m m
m m a
w a w
m M m
m mw
) m
m. w w
w . w
n m
m M mm
m ma...-
0 .. 00

 

 

 

 

 

 

 

 

 

((1)

Figure 5.8 Test 1 —ICA implementation with ECT field data: (a) Horizontal ECT data, (b)

Horizontal Separated Signal, (c) Vertical ECI‘ data, (d) Vertical Separated Signal

)

(c

71

 

 

015 I i i 1 0% i i i T r If i

 

 

 

 

 

 

 

 

I i i j i 00 i i i i i r i
«ll-IS 0.1 0.05 0 0 05 0.1 43.02 0 0.02 0.04 0.03 0.00 0.1 0.12 0.14
(a) (b)

 

““8 i l I I I 0‘
0.075 .
008 : ........ : ________ :
007 : .
0.085
0.08
0.06
0055
005
0.045
004

0.035

 

 

 

0 ' I ' ' '
03.04 0.05 0.08 0.07 0.08 0.09 0.1 0.11 0.12
(d)

 

003
0

 

Figure 5.9 Test 1 - Impedance trajectories of ECT field data at frequency: (a) 35ltHz, (b) 200kilz, (c)
400kilz, (d) 600kllz

72

 

 

 

 

 

 

 

 

 

 

Figure 5.10 Test 1 - Impedance trajectories of the reconstructed signal: (a) Defect, (b) Noise

5.2.2. Extraction of the Source Signal from Defect and TSP

Figures 5.11 (a) and (b) show another example of field data from steam generator (SG)

inspection MFECT data at 35kilz, 200kllz, 400kilz and 600kilz. The horizontal channel MFEC

data in Figure 5.11 includes a support plate characterized by two pulses at the beginning and
at the end of data. The presence of support plate makes it difficult to distinguish defect
component from the impedance plane trajectories of MFEC data shown in Figure 5.12.
After preprocessing using RST transformation, the signals are applied to the ICA
algorithm. The horizontal components of extracted source signals are shown in Figure 5.11
(b). The first extracted signal is interpreted as the horizontal component of TSP since it
exhibits large indication at the location of TSP. The third extracted signal is considered to be
the defect signal since it has a big indication at the location of defect. The horizontal or
vertical signal of a true defect signal is typically characterized by two spatially separated

peaks arising from the two edge points of the defect. The same procedure is implemented

73

with the vertical components of MFEC data. The vertical channel signals of MFEC data are
shown in Figure 5.11 (c) and extracted signals are presented in Figure 5.11 (d). In the
signals presented in Figure 5.11 (d), we see that the second signal has a distribution that
corresponds to the defect signal. The fourth signal in Figure 5.11 (d) is seen to be that of the
vertical component of TSP signal.

The horizontal and vertical signals of defect component yield the impedance plane
trajectory shown in Figure 5.13 (a). The impedance trajectory of TSP is presented in Figure
5.13 (b). From the extracted impedance trajectory of defect, we can distinguish that this
defect is categorized into wear because the impedance trajectory of defect has only one side
loop. Support plate is extracted as another independent component and the impedance

trajectory is shown in Figure 5.13 (b).

74

 

 

 

 

90

III
I

70

 

)

 

(c

203040505070809]
2030405080

10
10

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2030405050
(

10

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

((1)

Figure 5.11 Test 2 — ICA implementation with ECT field data: (a) Horizontal ECT data, (b) Vertical

ECT data, (c) Separated horizontal signals, ((1) Separated vertical signals
7 5

(b)

 

 

 

-----'-----—

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0.0% 1 1 1 i i 1 1 1 1 0 DI) 1 1 i 1 1 1
.l -0 [B 0&5 004 002 0 0.02 0.04 0.00 0.08 0] .2 015 -0.1 0 05 0 0.05 0.l 0.15
(b) ((1)

Figure 5.12 Test 2 - Impedance trajectories of ECT field data at frequency: (a) 35kiiz, (b) 200kiiz,
(c) 400kiiz, (d) 600kilz

76

 

 

 

 

 

 

 

 

 

 

Figure 5.13 Test 2 - Impedance trajectories of the extracted signals: (a) Defect, (b) TSP

5.2.3. Extraction of the Source Signals from Defect and Noise

A third experimental signal from an ‘Outside Diameter Inter Granular Attack (ODIGA)’
is analyzed next. Figure 5.14 (a) presents the horizontal channel data and resulting source
signals are shown in Figure 5.14 (b). The signal at the top of Figure 5.14 (b) is recognized as
the defect signal. The third signal is considered to be noise signal. The same procedure is
implemented on the vertical channel MFEC data. MFEC data and the extracted signals are
shown in Figure 5.14 (c) and ((1), respectively. Figure 5.14 ((1) shows the defect signal as the
second source signal. The third extracted signal in Figure 5.14 (d) is considered as noise
signal.

The impedance plane trajectory of MFEC data is shown in Figure 5.15. 4 impedance
plane trajectories obtained at different excitation frequencies are presented and it is seen that

the defect signal is barely visible. Notice that in contrast to the ideal noise free condition, the

77

 

background signal (pertaining to the defect region) has significant amplitude. However, the
impedance plane trajectories synthesized from the extracted horizontal and vertical defect

signals clearly show the presence of defect in Figure 5.16.

78

 

 

.....

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1C0 ill 140

 

 

      
 

 

 

 

 

 

 

 

(C) ((0

Figure 5.14 Test 3 — ICA implementation with ECT field data: (a) Horizontal ECT data, (b)
Horizontal Separated Signal, (c) Vertical ECT data, ((1) Vertical Separated Signal

79

 

0 015 01 015 02 025 0.3 035 04 0.45 0E0 009 01 011 012 0.13 014 015 015 UV

(c) ((1)

Figure 5.15 Test 3 — Impedance trajectories of ECT field data at frequency: (a) 35kilz, (b) 200kilz,
(c) 400kllz, (d) 600kilz

80

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 5.16 Test 3 — Impedance trajectories of the reconstructed signal: (a) Defect, (b) Noise

81

5.2.4. Use of ICA for defect detector

The ICA algorithm was applied to a database of steam generator (SG) tube inspection
signals for the purpose of defect detector. In this experiment, the separated signals are

tested for presence of a defect using a criterion function depicted as
= —— 5.1
A1 S [x] ( )

where s[] is standard deviation, x is the entire source signal (100 points) and d is the
segment of x at the defect location (15~20 points).

Using ,u a decision is made to classify a signal as “defect” or “no defect”. The
database consists of a total of 110 signals with only TSP and 27 signals comprising TSP
and defect. The classification performance is summarized in Table 5.2. The performance

showed a detection probability of about 90% (24 out of 27) and probability of false call

of about 12% (14 out of 110).

 

 

 

 

Cate or Predicted Predicted Total
g y Detect No Defect
True
Defect 24 3 27
True
No Defect 96 14 “0

 

 

 

 

Table 5.2 Statistics on the ICA performance

82

 

In order to compare the performance of ICA with that of conventional mixing
algorithm the second database was used for which results obtained using mixing
algorithm was available. The performance of ICA and mixing algorithm are summarized
in Table 5.3 and Table 5.4, respectively. It can be seen from theses results that ICA

performs better than mixing algorithm both in terms of detection probability and number

of false calls.

 

 

 

 

 

 

 

 

 

 

 

 

P3222?“ £26325;

Digit 30 4 34

NOTIgleEfect 8 30 38
Table 5.3 Statistics on the ICA performance

Pgszzztd £3332:

1322:; 25 9 34

NoTlguetfect 10 28 38

 

 

 

 

Table 5.4 Statistics on the mixing algorithm performance

83

 

 

6. CONCLUSION AND FUTURE WORK

6.1 . Conclusion

The application of ICA in enhancing differential eddy current probe signals from
flaws in the vicinity of tube support plate (TSP) in steam generator tubing is evaluated in
this thesis. ICA is a general-purpose statistical technique in which observed random data
from mixture of sources are decomposed into components that are maximally
independent from each other. The proposed Fast ICA algorithm gives a computationally
efficient method performing the feature extraction from multi-frequency eddy current
(MFEC) data. The algorithm has been applied to both, simulated data using finite
element method (FEM) and experimental MFEC field data. The results indicate that
independent component analysis (ICA) is a feasible approach for extracting signals from
flaws in steam generator tubes, even when the flaw is close to support plates. The

algorithm can also be used to enhance the flaw signal in the presence of noise.

6.2. Future Work

Future work will focus on a study of the following issues.

0 Test the algorithm on more exhaustive data set.

0 More robust and automated method to pick up the defect signal from the output of
ICA algorithm is required.

0 The performance of ICA algorithm depends largely on the linearity of MFECT
data. It is assumed that the observed MFECT data is a linear mixture of the source

signals. The validity of this assumption at different frequencies and its defect on

84

the result should be further studied. An alternate approach is to develop ICA
model to solve the nonlinear mixture problem.

0 ICA is a statistical method to extract source signal using a method that minimizes
gaussianity. Due to this property, ICA cannot be used to extract a signal that has
Gaussian distribution. In this paper, Gaussian noise signal is extracted using
residual theory. However, this is useful when only one Gaussian noise component
is present in the mixture. In this thesis, only one Gaussian measurement noise is
assumed to be present in the MFECT data. However, in order to model variation
in lift-off and surface roughness, which usually have Gaussian distributions
additional work need to be done. A simple approach is to eliminate noise signal

before applying ICA algorithm.

85

 

References

[1] Kenneth. C. Weston, “Energy Conversion: Chapter 10. Nuclear Power Plants”, e-book,
available at http://www.personal.utulsa.edu/~kenneth-weston/”

[2] S. Udpa, L. Udpa, “Eddy Current Nondestructive Evaluation” in Wiley Encyclopedia
of Electrical and Electronics Engineering, edited by John G. Webster, Wiley-
Interscience Publication, New York, 1999,Vol. 6.

[3] SP. Sullivan, V.S. Cecco, L.S. Obrutsky, J .R. Lakhan and AH. Park, “Validating Eddy
Current Array Probes for Inspecting Steam Generator Tubes”, European Commission
JRC, Institute for Advanced Materials, 1998 January, Vol.3, No.1

[4] K. Arunchalam, P. Ramuhali, L. Udpa, S. Udpa, “Nonlinear Mixing Algorithm For
Suppression of TSP Signals In Bobbin Coil Eddy Current Data,” Review of Progress of
Progress in QNDE, Vol. 21, D. 0. Thompson and D. E. Chimenti (Eds), American
Institute of Physics, 2002, pp.631-638.

[5] C. Jutten and J. Herault, “Blind separation of sources, part I : An adaptive algorithm
based on neurornimetic architecture,” Signal Processing 24, 1~10, 1991.

[6] C. Jutten, J. Herault and P. Comon, “Blind separation of sources, Part II : Problem
statement,” Signal Processing, Vol. 24, 11~20, 1991.

[7] A. J. Bell and T. J. Sejnowski, “An Information-Maximization Approach to Blind
Separation and Blind Deconvolution,” Neural Computation 7, 1129~1159, 1995.

[8] S. Amati, “Narural Gradient works efficiently in learning,” Neural Computation, 10,
251~276, 1997.

[9] Te-Won Lee, M. Girolami and T. J. Sejnowski, “Independent Component Analysis
using an Extended Infomax Algorithm for Mixed Sub-Gaussian and Super Gaussian
Sources,” 4‘h Joint Synposium on Neural Computation, 7, 132-140, Institute for
Neural Computation, 1997.

[10] A. Hyv‘arinen, “Fast and Robust Fixed-Point Algorithms for Independent Component
Analysis,” IEEE Trans. on Neural Computation, 10(3):626-634, 1995.

[11] A. Hyv'arinen and E. Oja, “A fast fixed-point algorithm for independent component
analysis,” IEEE Trans. on Neural Computation, 9(7):l483-l492, 1997.

[12] Udpa, L. “Imaging of Electromagnetic NDT Phenomena.” Ph.D. Dissertation,
Colorado State University, Ft. Collins, 1983.

[13] SN. Rajesh, “Probability of detection models for eddy current NDE method.”

86

Master Thesis, Ames, 1993.

[14] Cecco V.S. and Van Drunen, "Recognizing the Scope of Eddy Current Testing", in
Research Techniaues in Nondestructive Testing, Vol. 8, 1985, pp 269-301, ed. by RS.
Sharpe, Academic Press.

[15] Mundis, J. A. and DeYoung, G W. “Solution to NDE Problems in PWR Steam
Generators — An Overview.” Quantitative NDE in the Nuclear Industry. Edited by RB.
Clough. Metals Park, Ohio: American Society for Metals, 1983. P.93-98.

[16] Brown, S. D. “Multifrequency/Multiparameter Eddy Current Steam Generator
NDE.” Quantitative NDE in the Nuclear Industry. Edited by RB. Clough, Metals
Park, Ohio: American Society for Metals, 1983. p.99-105.

[17] Smith, J ., Dodd, C., and Chitwood, L. “Multifrequency Eddy Current Examination of
Seam Weld in Stell Sheath.” Materials Evaluation Vol.43. November, 1985. p. 1566-
1572.

[18] Libby, H. L. Introduction to Electromagnetic Nondestructive Test Models. New
York: John Wiley and Sons, Inc., 1971.

[19] Saglio, R. and Pigeon, M. “Concepts of Multifrequency Eddy Current Testing.”
Electromagnetic Testing. Vol. 4, 2nd ed. Edited by P. McIntire, R. McMaster.
Columbus, Ohio: American Society for Nondestructive testing, 1986. p.592-605.

[20] Solete, J ., Udpa, L., and Lord. W. “Multifrequency Eddy Current Testing of Steam
Generator Tubes Using Optimal Affine Transform.” Review of Progress in
Quantitative Nondestructive Evaluation, Vol. 7A. Edited by DD. Thompson and DE.
Chimenti. New York: Plenum Press, 1988. p.821-830.

[21] Libby, H. L. “Eddy Current Test for Tubing Flaws in Support Regions.” Research
Techniques in Nondestructive Testing, Vol. H. Edited by RS. Sharpe. London:
Academic Press, 1973. p151-184.

[22] CL. Owsley, “A Multifrequency Eddy Current NDE System.” Master Thesis, Ames,
1996.

[23] A. Leon-Garcia, “Probability and Random Processes for Electrical Engineering,”
Addison Wesley, 2nd edition,

[24] J .H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford University Press, Oxford,
1965.

[25] M. Kendall and A. Stuart, The advanced theory of statistics, Mcgraw Hill, 1969.

[26] P. Comon, “Independent component analysis - a new concept”, Signal Processing,

87

36:287-314, 1994.

[27] MC. Jones and R. Sibson, “What is projection pursuit?”, J. of the Royal Statistical
Society, set: A, 15011-36, 1987.

[28] A. Hyv'arinen. New approximations of differential entropy for independent
component analysis and projection pursuit. In advances in Neural Information
Processing Systems, Vol. 10, p 273-279, MIT Press, 1998.

[29] G. Arfken, Mathematical Method for Physicists, 3rd edition, American Press, 1985.

[30] L. V. Kantorovich and GP. Akilov, Functional Analysis in Normed Spaces,
Pergamon Press, Elmsford, New York, 1964. '

[31] J. Karhunen, E. Oja, L. Wang, R. Vigario, and J. Joutsensalo. “A class of neural
networks for independent component analysis”, IEEE Trans. on Neural Networks,
8(3):486-504, 1997.

[32] A. Hyv'arinen, “One-unit contrast functions for independent component analysis: A
statistical analysis”, In Neural Networks for Signal Processing VII (Proc. IEEE
Workshop on Neural Networks for Signal Processing), pages 388-397, Amelia Island,
Florida, 1997.

[33] M. Kendall and A. Stuart, “The Advanced Theory of Statistics”, Charles Griffin &
Company, 1958.

[34] The Eddy Current Simulation with Finite Element Method in MATLAB code.
Available at “http://www.egr.msu.edu/~ zengzhiwl”.

[35] Lord, W. and Palanisamy, R. “Magnetic Probe Inspection of steam Generator
Tubing.” Materials Evaluation, Vol. 38, No. 5, May 1980, p478-485

[36] J. Freund, Mathematical statistics, Prentice Hall, Upper Saddle River, New Jersey,
six edition, 1999.

88

 

[[11]];111311][[11

 

 

flag——