# Nondestructive Inspection of Thin Basalt Fiber Reinforced Composites Using Combined Terahertz Imaging and Infrared Thermography

Authors: Przemyslaw Lopato; Grzegorz Psuj; Barbara Szymanik
Journal: Advances in Materials Science and Engineering
 (2016)

---

## Abstract

The inspection of thin basalt fiber reinforced composite materials was carried out using two nondestructive methods: terahertz time domain imaging and infrared thermography. In order to combine the information about the defects arising in examined materials the inspection results were parametrized. In order to acquire more information content, new approximation based features are proposed. Then, a knowledge extraction based multivariate analysis of preselected features’ vector was carried out. Finally, in order to integrate features distributions of representing different dynamic level of information, a multiresolution wavelet based data fusion algorithm was applied. The results are presented and discussed.


---

## Body

## 1. Introduction

The increasing application of composite materials in modern structures of the aerospace, construction, and automotive industries puts new demands on the quality control units. Mostly, nondestructive testing (NDT) methods relate to metals inspection. Therefore, there is a need for development of new or adaptation of existing methods to composite materials properties [1]. The selection of the most appropriate method for testing of composite materials depends on the nature of the fibers and the polymer matrix [2]. However, no single NDT technique allows a full assessment of the material under test structure’s integrity. In essence, each method presents some limits of detection and determination of the nonhomogeneity of the material [3]. In order to increase the probability of defects’ detection and a proper evaluation of the material’s structure, a need of fusion of data obtained from various testing methods is arising [4, 5].In this paper the inspection of thin basalt fiber reinforced composite (BFRC) materials (shortly described in Section2) is carried out using two electromagnetic NDT methods. The selected techniques allow assessing different aspects of the defects arising in specimen. First, the acquired data were processed and parameterized. The analysis of each method’s results allowed evaluating the usability of each applied testing technique for the inspection of basalt composite. Then, in order to enhance the proper evaluation of the material stage and combine information supplied by each method, a knowledge extraction based multivariate analysis of features vector was performed followed by multiresolution decomposition based data fusion. Finally, the performance of the data fusion was assessed and compared with the results obtained for each method separately.


## 2. Basalt Fiber Reinforced Composite

Basalt fiber is a perspective polymer-reinforcing material and can be applied in polymer matrix composites instead of glass fiber. Fabrics of varying surface densities are made depending upon the application type and are in the range from 160 g/m2 to 1100 g/m2 [6]. The comparison of selected properties of basalt and glass fibers is shown in Table 1. One can observe an advantage of basalt over glass fibers in case of mechanical and thermal properties.Table 1
Selected properties of basalt and glass fibers [7, 8].
Properties
Unit
Basalt
E-Glass
S-Glass
Density
g/cm3
2.7
2.57
2.48
Thermal linear expansion coefficient
ppm/°C
8.0
5.4
2.9
Tensile strength
GPa
4.84
3.45
4.71
Elastic modulus
GPa
89
77
89
Elongation at break
%
3.15
4.7
5.6
Compression strength
MPa
3792
3033
3033
Maximum application temperature
°C
982°
650°
650°
Sustained operating temperature
°C
820°
480°
480°
Minimum operating temperature
°C
−260°
−60°
−60°
Melting temperature
°C
1450°
1120°
1120°
Absorption of humidity (65% RAH)
%
<0.1
<0.1
<0.1
Emissivity at 20°C
—
0.92
0.65
0.65
Thermal conductivity
W/m K
0.031–0.038
0.034–0.04
0.034–0.04The photo of the utilized basalt fibers reinforced composite specimen is presented in Figure1. The material consists of six layers of basalt fiber fabric within a polyester resin (Polimal 1094 AWTP-1). The thickness of the samples was 1–1,2 mm; thus thickness of single layer was less than 200 μm. The evaluated basalt test specimens were exposed to mechanical impacts of 2 and 5J energy (sample s
I
M
P
_
2
J and s
I
M
P
_
5
J). During samples production process, artificial inclusions (mica film of 145 μm thickness) were also introduced under selected layers (samples s
I
N
C
_
A and s
I
N
C
_
B).Figure 1
Photo of the examined specimens
I
M
C
_
5
J (basalt reinforced composite material exposed to impact energy): (a) general view and (b) enlargement of measured region.
(a)
(b)


## 3. Multisource Inspection Procedure

In the process of inspection of thin basalt fiber reinforced composite materials, two nondestructive methods were utilized. The measurements were carried out using terahertz (THz) inspection technique and active infrared thermography (IRT). This approach enabled more efficient detection of defects based on changes of both electrical and thermal parameters.Terahertz electromagnetic radiation enables noninvasive, nonionizing, and noncontact examination of dielectric materials such as plastics, dry wood, explosives ceramics, foams, and composites [9, 10]. Any defect which noticeably disturbs refractive index can be detected, for example, void, delamination, inclusion, material inhomogeneity (fiber/matrix distribution), surface roughness, fiber waviness, and internal interfaces between layers (in layered structures). Principle of operation is similar to radar: electromagnetic wave in terahertz frequency range is emitted, focused on surface of examined material, and acquired after interaction (transmission or reflection) with material. In most cases, defects are detected by reflection and transmission imaging based on pulsed terahertz TDS (Time Domain Spectroscopy) [11]. In this case electromagnetic excitation is in form of very short (order of picoseconds) pulses. The method is well suited for evaluation of layered materials [12]. Due to the different intrinsic impedance values, presence of defect or interface between separate layers causes reflection of the incident THz pulse and attenuation of the transmitted one. Differences in delays of the propagated pulses and their echo (delayed layer reflections) enable characterization of the inner structure state. Very short pulses contain wide frequency bandwidth (0.05–3 THz) and therefore it is possible to carry one single point broadband measurements [13]. Photo of utilized terahertz measuring setup is presented in Figure 2(a).Figure 2
The measurement setups: (a) terahertz inspection system and (b) IRT system.
(a)
(b)Active infrared technique (IRT) finds its application in inspection of variety of materials not only composite. In this method thermal nonequilibrium of the system is obtained using external energy source [2, 3]. The temperature difference in examined samples can be induced using several excitation methods [3]. The choice of proper method depends on the tested materials’ properties. For conductive materials the induction heating can be used, while in case of nonconductive composites materials other energy sources have to be utilized, for example, halogen lamps, convective heating, and microwaves. The temperature distribution at the specimen’s surface can be observed during both heating and cooling stage using thermovision camera (in this study, the FLIR A325 camera was used). Defects are detected as (depending on the damage type) under- or overheated spots. In this study, the halogen lamp heating system was used to induce the temperature differences within the examined composite samples. In this method the specimen is placed in front of two halogen lamps of maximum power 2000 W each, positioned in the manner which ensures uniform illumination of the sample. Photo of IRT measuring setup is presented in Figure 2(b).


## 4. Results of Measurements and Parameters Calculation

### 4.1. Terahertz Inspection Results

Terahertz inspection was performed using imaging system based on TRay4000 pulsed spectroscope of Picometrix. Measuring head was scanned over the surface of a sample, and for each(
x
,
y
) position time domain signal s
(
t
) was acquired. An exemplary result of pulsed terahertz inspection of thin basalt fiber reinforced composite without defects in given (
x
,
y
) position is shown in Figure 3. One can distinguish seven peaks caused by step change of wave impedance Z
_
w on layers boundaries. First peak caused by front surface reflection (FSR) is negative, and the other—including the last one caused by back surface reflection (BSR)—is positive. An existence of defects changes refractive index distribution and/or introduces additional boundaries. Those effects influence obtained signal s
(
t
). Results of line scan (B-scan) over defected area are presented in Figure 4.Figure 3
Exemplary A-scan signals
(
t
) acquired over thin basalt fiber sample without defects: FSR: pulse caused by front surface reflection; BSR: pulse caused by back surface reflection.Figure 4
Selected B-scan results of pulsed terahertz inspection: (a) 5J impact caused defect and (b) thin inclusion under 1th layer.
(a)
(b)In order to detect hidden defects inside of examined structure, parametrization procedure was utilized. The following parameters ofs
(
t
) signal were calculated for each (
x
,
y
) position:(1)
M
a
x
v
a
l
x
,
y
=
max
⁡
s
t
x
,
y
,
M
i
n
v
a
l
x
,
y
=
min
⁡
s
t
x
,
y
,
S
u
m
a
b
s
x
,
y
=
∫
t
1
t
2
s
t
x
,
y
d
t
,
S
u
m
d
e
r
i
v
x
,
y
=
∫
t
1
t
2
d
2
s
t
x
,
y
d
t
2
d
t
,
M
a
x
d
e
r
i
v
x
,
y
=
max
⁡
d
2
s
t
x
,
y
d
t
2
,where s
t
(
x
,
y
) is the A-scan signal acquired in given position (
x
,
y
), t is a time, and t
1, t
2 are start and stop times for integration (time gating of A-scan signal).Additionally, approximation procedure was utilized. The material response signals
(
t
) is approximated using the following function: (2)
s
a
p
p
t
=
∑
k
=
0
2
p
k
t
k
+
∑
l
=
1
7
a
l
e
-
t
-
b
l
/
c
l
2
,where t is a time, p
k is a polynomial component coefficient, and a
l, b
l, c
l are Gaussian component coefficients.The approximation function is dedicated to thin material consisting of 6 layers. It consists of combination of polynomial and Gaussian components. The resulting number of coefficients is high in case of proposed approximation function, but it enables very precise evaluation of materials state at different depths (layers). Spatial distributions of selected parameters in case of various defects are presented in further sections. Such distributions will be utilized by data fusion algorithm.


### 4.2. Infrared Thermography Inspection

During active thermographic inspection the halogen lamps heating was utilized. Taking into consideration the possibility of heating phase observation, the Pulsed Phase Thermography (PPT) was applied in order to extract the information about defects. This method combines the experimental procedure used in Pulsed Thermography (PT), with signal analysis used in Modulated Thermography (MT). First, the thermograms’ sequence is recorded while the heat pulse is applied to examined specimen and for the certain time after the heating phase, to observe the cooling phase as well. In this study the sample was heated for 30 seconds and then natural process of cooling was observed for additional 30 seconds. The recorded sequence consisted of 60 thermograms (recording frequency was set to 1 image per second). Then, the analysis of obtained sequence is based on Discrete Fourier Transform (DFT), which allows evaluating the output as the combination of phase (P
h
a
I
R
T
n
(
x
,
y
)) and amplitude (A
m
p
I
R
T
n
(
x
,
y
)) images of nth harmonic of IRT image sequence [14–16]. The selected results of the thermographic inspection of the basalt fiber reinforced composite material achieved for halogen lamps heating method are presented in Figure 5.Figure 5
Selected results of thermographic inspection: (a) amplitude imageA
m
p
I
R
T
4
(
x
,
y
) and (b) phase image P
h
a
I
R
T
4
(
x
,
y
) obtained using halogen lamps heating system.Obtained amplitudeA
m
p
I
R
T
(
x
,
y
) and phase P
h
a
I
R
T
(
x
,
y
) distributions can be utilized as parameters for further analysis—data fusion and defects detection procedure. Another interesting parameter that can be calculated from frequency transformed time sequence of images is total harmonic distortion, defined as(3)
T
H
D
x
,
y
=
∑
n
=
2
N
A
m
p
I
R
T
n
x
,
y
2
A
m
p
I
R
T
1
x
,
y
,where n is the harmonics number of image sequences.Temperature changes in case of selected(
x
,
y
) position are shown in Figure 6. During first 31 seconds temperature is growing semilinearly (heating phase). After this heat source is turned off and the sample temperature is decreasing (cooling phase). Based on heating phase observation the following linear approximation model is proposed:(4)
T
x
,
y
,
t
=
α
I
R
T
x
,
y
·
t
+
β
I
R
T
x
,
y
,where T
(
x
,
y
) is the temperature at given position (
x
,
y
), α
I
R
T
(
x
,
y
) is the multiplicative coefficient (slope of the heating phase line), and β
I
R
T
(
x
,
y
) is the additive coefficient.Figure 6
Measured temperature changes of selected(
x
,
y
) position.Selected parameters will be utilized by data fusion procedure presented in next sections.


## 4.1. Terahertz Inspection Results

Terahertz inspection was performed using imaging system based on TRay4000 pulsed spectroscope of Picometrix. Measuring head was scanned over the surface of a sample, and for each(
x
,
y
) position time domain signal s
(
t
) was acquired. An exemplary result of pulsed terahertz inspection of thin basalt fiber reinforced composite without defects in given (
x
,
y
) position is shown in Figure 3. One can distinguish seven peaks caused by step change of wave impedance Z
_
w on layers boundaries. First peak caused by front surface reflection (FSR) is negative, and the other—including the last one caused by back surface reflection (BSR)—is positive. An existence of defects changes refractive index distribution and/or introduces additional boundaries. Those effects influence obtained signal s
(
t
). Results of line scan (B-scan) over defected area are presented in Figure 4.Figure 3
Exemplary A-scan signals
(
t
) acquired over thin basalt fiber sample without defects: FSR: pulse caused by front surface reflection; BSR: pulse caused by back surface reflection.Figure 4
Selected B-scan results of pulsed terahertz inspection: (a) 5J impact caused defect and (b) thin inclusion under 1th layer.
(a)
(b)In order to detect hidden defects inside of examined structure, parametrization procedure was utilized. The following parameters ofs
(
t
) signal were calculated for each (
x
,
y
) position:(1)
M
a
x
v
a
l
x
,
y
=
max
⁡
s
t
x
,
y
,
M
i
n
v
a
l
x
,
y
=
min
⁡
s
t
x
,
y
,
S
u
m
a
b
s
x
,
y
=
∫
t
1
t
2
s
t
x
,
y
d
t
,
S
u
m
d
e
r
i
v
x
,
y
=
∫
t
1
t
2
d
2
s
t
x
,
y
d
t
2
d
t
,
M
a
x
d
e
r
i
v
x
,
y
=
max
⁡
d
2
s
t
x
,
y
d
t
2
,where s
t
(
x
,
y
) is the A-scan signal acquired in given position (
x
,
y
), t is a time, and t
1, t
2 are start and stop times for integration (time gating of A-scan signal).Additionally, approximation procedure was utilized. The material response signals
(
t
) is approximated using the following function: (2)
s
a
p
p
t
=
∑
k
=
0
2
p
k
t
k
+
∑
l
=
1
7
a
l
e
-
t
-
b
l
/
c
l
2
,where t is a time, p
k is a polynomial component coefficient, and a
l, b
l, c
l are Gaussian component coefficients.The approximation function is dedicated to thin material consisting of 6 layers. It consists of combination of polynomial and Gaussian components. The resulting number of coefficients is high in case of proposed approximation function, but it enables very precise evaluation of materials state at different depths (layers). Spatial distributions of selected parameters in case of various defects are presented in further sections. Such distributions will be utilized by data fusion algorithm.


## 4.2. Infrared Thermography Inspection

During active thermographic inspection the halogen lamps heating was utilized. Taking into consideration the possibility of heating phase observation, the Pulsed Phase Thermography (PPT) was applied in order to extract the information about defects. This method combines the experimental procedure used in Pulsed Thermography (PT), with signal analysis used in Modulated Thermography (MT). First, the thermograms’ sequence is recorded while the heat pulse is applied to examined specimen and for the certain time after the heating phase, to observe the cooling phase as well. In this study the sample was heated for 30 seconds and then natural process of cooling was observed for additional 30 seconds. The recorded sequence consisted of 60 thermograms (recording frequency was set to 1 image per second). Then, the analysis of obtained sequence is based on Discrete Fourier Transform (DFT), which allows evaluating the output as the combination of phase (P
h
a
I
R
T
n
(
x
,
y
)) and amplitude (A
m
p
I
R
T
n
(
x
,
y
)) images of nth harmonic of IRT image sequence [14–16]. The selected results of the thermographic inspection of the basalt fiber reinforced composite material achieved for halogen lamps heating method are presented in Figure 5.Figure 5
Selected results of thermographic inspection: (a) amplitude imageA
m
p
I
R
T
4
(
x
,
y
) and (b) phase image P
h
a
I
R
T
4
(
x
,
y
) obtained using halogen lamps heating system.Obtained amplitudeA
m
p
I
R
T
(
x
,
y
) and phase P
h
a
I
R
T
(
x
,
y
) distributions can be utilized as parameters for further analysis—data fusion and defects detection procedure. Another interesting parameter that can be calculated from frequency transformed time sequence of images is total harmonic distortion, defined as(3)
T
H
D
x
,
y
=
∑
n
=
2
N
A
m
p
I
R
T
n
x
,
y
2
A
m
p
I
R
T
1
x
,
y
,where n is the harmonics number of image sequences.Temperature changes in case of selected(
x
,
y
) position are shown in Figure 6. During first 31 seconds temperature is growing semilinearly (heating phase). After this heat source is turned off and the sample temperature is decreasing (cooling phase). Based on heating phase observation the following linear approximation model is proposed:(4)
T
x
,
y
,
t
=
α
I
R
T
x
,
y
·
t
+
β
I
R
T
x
,
y
,where T
(
x
,
y
) is the temperature at given position (
x
,
y
), α
I
R
T
(
x
,
y
) is the multiplicative coefficient (slope of the heating phase line), and β
I
R
T
(
x
,
y
) is the additive coefficient.Figure 6
Measured temperature changes of selected(
x
,
y
) position.Selected parameters will be utilized by data fusion procedure presented in next sections.


## 5. Evaluation Using Multiple Features Data Fusion

None of the methods allow us to fully assess the integration stage of the structure. Both inspection techniques bring crucial information in the process of the nondestructive structural integrity imaging. Even within a single method different features of acquired data can provide one with unique information pertaining to current stage of the material’s structure; that is, time response allows monitoring material’s condition at different depth. Therefore, in order to fully conduct the process of nondestructive imaging of the structure state data fusion of multiple features extracted from both methods results was carried out. The multiple sources’ inspection using the methods coming from different physical origins and utilizing various phenomena to visualize the structure state make the evaluation system more robust for the unwanted disturbances and increase the overall performance of the accession process [4, 5, 17].The block diagram of the algorithm utilized in this paper is presented in Figure7. First, features extracted from both inspection methods are transformed into common representation format undergoing the spatial and resolution registration process. Then, a representative set of all transformed features allowing monitoring of different aspects of material structure were chosen for final definition of database feature vector. Finally a knowledge extraction for the purpose of vector dimension reduction was carried out and a multiresolution data fusion was proceeded for final material’s structure imaging.Figure 7
The functional block diagram of the data fusion process.


### 5.1. Spatial Registration of Extracted Features

In order to transform the features of both methods the data registration process must be carried out. Therefore several aspects have to be considered to carry out the registration process. The most important factors that should be taken under consideration are each method’s sensing element geometrical distortions and position (alignment and rotation) with respect to the evaluated material or measuring resolution. The block diagram of the applied algorithm is presented in Figure8. High resolution photo of the sample was used as the reference image. First, dedicated metallic markers (control points CP) detectable by both methods were used. Then the CP matching process was proceeded. Taking into consideration that during the experiments different position of the sensing element in respect to the examined composite was applied, different transformation model had to be considered. Visualization of the measuring procedure and setup of the sensing devices with respect to sample surface were presented in Figure 9. In case of terahertz inspection both wave source and detector were placed parallel to the composite material (parallel projection), while in case of the thermographic inspection the infrared camera was observing a sample from some perspective (perspective diametric projection). Therefore, considering the projection type, the similarity and the projective transformation of data were applied, respectively, in first and second case.Figure 8
The functional block diagram of the data registration process.Figure 9
Schematic view of the experiments setup and spatial relationship between both methods sensing units: (a) top view and (b) side view.
(a)
(b)Both conversions are a standard geometrical transformation utilized in image processing algorithms and are defined as follows [18, 19]:(i)
Similarity transformation:(5)
x
′
=
s
x
cos
⁡
α
-
y
sin
⁡
α
+
t
x
,
y
′
=
s
x
sin
⁡
α
+
y
cos
⁡
α
+
t
y
,
where(
x
,
y
) are original coordinates, (
x
′
,
y
′
) are new coordinates, (
t
x
,
t
y
) are translation coefficients which specify the movement of the system’s center, s is a scaling factor, and α corresponds to a rotation angle.
(ii)
Perspective transformation:(6)
x
′
=
a
11
x
+
a
12
y
+
t
x
b
1
x
+
b
2
y
+
1
,
y
′
=
a
21
x
+
a
22
y
+
t
y
b
1
x
+
b
2
y
+
1
,
wherea
11, a
12, a
21, a
22 are coefficients responsible for rotation and scaling and b
1 and b
2 are coefficients defining the projection.The similarity is shape and angles preserving transformation, allowing scaling, rotation, and reflection operation. The projective one is used to transform an image perspective. It preserves collinearity and incidence; however it affects parallelism, length, and angle.After the procedure of matching of the control points (Figure10 presents the projection of each method’s results on CCD image of the sample), the images geometrical distortions were eliminated. Selected results of the registration process are shown in Figure 11.Figure 10
Results of data to inspection area matching: (a) photo, (b) photo with depicted exemplary THz result, and (c) photo with depicted exemplary IRT result.
(a)
(b)
(c)Figure 11
Results of the data registration process: IRT exemplary results (a) before and (b) after the registration process, THz exemplary results (c) before and (d) after the registration process; dpi before registration process, 40 (IRT) and 25 (THz), and after, 500.
(a)
(b)
(c)
(d)


### 5.2. Final Feature Vector Definition

After the data registration process the features distributions obtained for each single method were analyzed and preselected. During the preselection process two aspects were taken into account: first to minimize the amount of similar information constancy and second to preserve the possibility of observation of different aspects of structure state assessment.Finally the set of 20 features representing both methods were used to define the features’ vectorF for processing of the data fusion algorithm. The chosen feature list is presented in Table 2 while the features distributions obtained for impacted (s
I
M
P
_
2
J
,
s
I
M
P
_
5
J) and with inclusions (s
I
N
C
_
A
,
s
I
N
C
_
B) samples are shown in Figures 12, 13, 14, and 15, respectively.Table 2
Features’ vectorF of the database.
Feature number
Name
Measuring method
1
M
a
x
v
a
l
THz
2
M
i
n
v
a
l
THz
3
S
u
m
a
b
s
THz
4
S
u
m
d
e
r
i
v
THz
5
M
a
x
d
e
r
i
v
THz
6
a
1
THz
7
b
2
THz
8
b
4
THz
9
s
a
p
p(b
1)
THz
10
s
a
p
p(b
3)
THz
11
s
a
p
p(b
7)
THz
12
AmpIRT1
IRT
13
AmpIRT2
IRT
14
AmpIRT3
IRT
15
PhaIRT2
IRT
16
PhaIRT3
IRT
17
PhaIRT4
IRT
18
THD
IRT
19
α
I
R
T
IRT
20
β
I
R
T
IRTFigure 12
Results of the data registration process carried out for features vector obtained for sample impacted by 2J.Figure 13
Results of the data registration process carried out for features vector obtained for sample impacted by 5J.Figure 14
Results of the data registration process carried out for features vector obtained for sample with inclusion A.Figure 15
Results of the data registration process carried out for features vector obtained for sample with inclusion B.The THz signals features allow observing rapid changes in the structure of the examined composite materials. In particular inclusions are very well visible. It is possible to precisely localize rectangular shaped inclusions in case of both samples in distributions of features 4, 6, 7, 8, and 10. In case of samples
I
N
C
_
A during the production process additional defects were unintentionally introduced beside the inclusion. These defected areas can be visible as straight lines ranging from top left corner to bottom right in the distributions of features 1, 3, 8, and 11. For impacted samples the THz features rather allow indicating the localization where the highest damage occurred rather than the range of the defected area.Contrary to THz, the IRT imaging features are characterized by the low dynamics of the heating system. In case of composites, which generally have low value of thermal diffusivity factor, the heat conduction process within material is relatively slow. Therefore, the indications of defects observed at the sample’s surface are, in most cases, blurred. This effect can be noticed for samples with inclusions (it is not possible to identify the shape of the inclusions). On the other hand, the obtained features distributions make it possible to fully assess the range of the impacted area even in the sample impacted with lower energy. Despite the lower ratio between the response to defects and the background values than in case of the THz features, it is possible also to observe the additional defects in included samples
I
N
C
_
A (see feature 16, 18, and 19 distributions).


### 5.3. Multiresolution Decomposition Based Data Fusion Algorithm

The results of both methods represent different spatial dynamics of the feature changes. In case of the THz the inspecting wave is focused on local point of the material, so the response is also from that specific point, while in case of IRT testing, obviously, the heat generated by the halogen lamps cannot be focused at one point but is radiated to the sample’s surface and then is conducted evenly throughout the material volume. Therefore the thermal response in the specific point is disturbed by the surround. This results in lower spatial dynamics of the IRT imaging in comparison to THz one. Therefore, in order to preserve the information content of both rapid and slowly changing responses, a multiresolution decomposition (MRD) based data fusion algorithm was applied to the features vector.The block diagram of the procedure is presented in Figure16. First, in order to reduce the dimensionality of the problem, the knowledge extraction was carried out for the features’ vector. The purpose of this stage is to transform the original features’ vector into a new coordinate system of reduced dimensions level by extracting the most useful information from the original data. In order to carry out the multivariate analysis two algorithms were used: Independent Component Analysis (ICA) and Principal Component Analysis (PCA).Figure 16
Multiresolution decomposition based data fusion algorithm diagram; DWT and DWT−1: wavelet multiresolution decomposition and its inverse.The general idea of the ICA is to extract from multidimensional set of data the independent components under the assumption that the sources of data are statistically independent [19]. This means that each feature value gathered in the vector F is a combination of the original sources data. Then, the ICA analysis allows estimating the original sources from the collected features vector.The objective of PCA is to convert the database of original possibly correlated features into a database having new, statistically independent variables using orthogonal linear transformation [19, 20]. The new parameters called principal components pc are obtained from eigenvectors of covariance matrix of the original database and ordered in accordance with quantity of their variance (described by eigenvalues) in the database. The greater the variance of the variable is, the more significant it is. Typically the first few pc variables express the majority of database variance. Therefore they can be used to represent the whole database without losing a significant amount of information.In both cases the analysis of the two new components (combined featuresc
f
1 and c
f
2) was taken for further processing: two independent and first two principal components representing the greatest variance of the whole set.After the new components extraction MRD was carried out. The general idea of MRD is to represent the source of data with a collection of hierarchical representations of basis function at different resolutions (frequency bands). The original data is decomposed into a set of spatial frequency bandpass representatives, obtained by convolving and subsampling operations. The data can be decomposed using various methods [5]. In this paper the wavelet decomposition WT was utilized. WT is one of the most effective decomposition methods allowing us to obtain a good resolution in both time and frequency domain [21, 22]. After the decomposition the two distributions inputs at given sublevel are fused using the specified fusion rule. Then, the inverse transform DWT−1 is computed and the fused distribution is reconstructed. The whole process can be described by(7)
F
D
F
x
,
y
=
D
W
T
-
1
ϕ
D
W
T
F
C
1
x
,
y
,
D
W
T
F
C
2
x
,
y
,where ϕ is the fusion rule, F
C
1 and F
C
2 are components 1 and 2 obtained after extraction process, and F
D
F is the fused distribution.The result achieved for the proposed algorithm with the utilization of Daubechiesψ
D
10 wavelet function for the decomposition of c
f
1 and c
f
2 and the maximum selection as the fusion rule [22] are presented in Figures 17 and 18, respectively, for ICA and PCA knowledge extraction procedure. The general advantage of the used data fusion algorithm is that no matter which feature would provide the information about the defect it would be indicated in the fused distribution. Analyzing the obtained results one can notice that the fused distributions carry crucial information about the defects and the structure delivered by all features. All defects were indicated in single fused distributions, which can be especially seen in case of samples with inclusions s
I
N
C
_
A and s
I
N
C
_
B. Besides clearly visible rectangular shaped inclusions it is possible to observe additional defect of the structure in the neighborhood of the inclusions (Figures 17, 18(c), and 18(d)). Also in case of impacted samples the achieved results allow distinguishing the areas of different damage level. The applied knowledge extraction procedures allowed compressing information in much less data representations, and the multiresolution decomposition fused both rapid and slowly changing signals. This was useful once again especially in case of samples with inclusions. Inclusions are characterized mainly by the high dynamic signals’ components while additional anomalies in the composite structure being results of faulty production process are indicated in the low dynamic components.Figure 17
Results of the ICA-DWT data fusion algorithm obtained for sample: (a)s
I
M
P
_
2
J, (b) s
I
M
P
_
5
J, (c) s
I
N
C
_
A, and (d) s
I
N
C
_
B.
(a)
(b)
(c)
(d)Figure 18
Results of the PCA-DWT data fusion algorithm obtained for sample: (a)s
I
M
P
_
2
J, (b) s
I
M
P
_
5
J, (c) s
I
N
C
_
A, and (d) s
I
N
C
_
B.
(a)
(b)
(c)
(d)


## 5.1. Spatial Registration of Extracted Features

In order to transform the features of both methods the data registration process must be carried out. Therefore several aspects have to be considered to carry out the registration process. The most important factors that should be taken under consideration are each method’s sensing element geometrical distortions and position (alignment and rotation) with respect to the evaluated material or measuring resolution. The block diagram of the applied algorithm is presented in Figure8. High resolution photo of the sample was used as the reference image. First, dedicated metallic markers (control points CP) detectable by both methods were used. Then the CP matching process was proceeded. Taking into consideration that during the experiments different position of the sensing element in respect to the examined composite was applied, different transformation model had to be considered. Visualization of the measuring procedure and setup of the sensing devices with respect to sample surface were presented in Figure 9. In case of terahertz inspection both wave source and detector were placed parallel to the composite material (parallel projection), while in case of the thermographic inspection the infrared camera was observing a sample from some perspective (perspective diametric projection). Therefore, considering the projection type, the similarity and the projective transformation of data were applied, respectively, in first and second case.Figure 8
The functional block diagram of the data registration process.Figure 9
Schematic view of the experiments setup and spatial relationship between both methods sensing units: (a) top view and (b) side view.
(a)
(b)Both conversions are a standard geometrical transformation utilized in image processing algorithms and are defined as follows [18, 19]:(i)
Similarity transformation:(5)
x
′
=
s
x
cos
⁡
α
-
y
sin
⁡
α
+
t
x
,
y
′
=
s
x
sin
⁡
α
+
y
cos
⁡
α
+
t
y
,
where(
x
,
y
) are original coordinates, (
x
′
,
y
′
) are new coordinates, (
t
x
,
t
y
) are translation coefficients which specify the movement of the system’s center, s is a scaling factor, and α corresponds to a rotation angle.
(ii)
Perspective transformation:(6)
x
′
=
a
11
x
+
a
12
y
+
t
x
b
1
x
+
b
2
y
+
1
,
y
′
=
a
21
x
+
a
22
y
+
t
y
b
1
x
+
b
2
y
+
1
,
wherea
11, a
12, a
21, a
22 are coefficients responsible for rotation and scaling and b
1 and b
2 are coefficients defining the projection.The similarity is shape and angles preserving transformation, allowing scaling, rotation, and reflection operation. The projective one is used to transform an image perspective. It preserves collinearity and incidence; however it affects parallelism, length, and angle.After the procedure of matching of the control points (Figure10 presents the projection of each method’s results on CCD image of the sample), the images geometrical distortions were eliminated. Selected results of the registration process are shown in Figure 11.Figure 10
Results of data to inspection area matching: (a) photo, (b) photo with depicted exemplary THz result, and (c) photo with depicted exemplary IRT result.
(a)
(b)
(c)Figure 11
Results of the data registration process: IRT exemplary results (a) before and (b) after the registration process, THz exemplary results (c) before and (d) after the registration process; dpi before registration process, 40 (IRT) and 25 (THz), and after, 500.
(a)
(b)
(c)
(d)


## 5.2. Final Feature Vector Definition

After the data registration process the features distributions obtained for each single method were analyzed and preselected. During the preselection process two aspects were taken into account: first to minimize the amount of similar information constancy and second to preserve the possibility of observation of different aspects of structure state assessment.Finally the set of 20 features representing both methods were used to define the features’ vectorF for processing of the data fusion algorithm. The chosen feature list is presented in Table 2 while the features distributions obtained for impacted (s
I
M
P
_
2
J
,
s
I
M
P
_
5
J) and with inclusions (s
I
N
C
_
A
,
s
I
N
C
_
B) samples are shown in Figures 12, 13, 14, and 15, respectively.Table 2
Features’ vectorF of the database.
Feature number
Name
Measuring method
1
M
a
x
v
a
l
THz
2
M
i
n
v
a
l
THz
3
S
u
m
a
b
s
THz
4
S
u
m
d
e
r
i
v
THz
5
M
a
x
d
e
r
i
v
THz
6
a
1
THz
7
b
2
THz
8
b
4
THz
9
s
a
p
p(b
1)
THz
10
s
a
p
p(b
3)
THz
11
s
a
p
p(b
7)
THz
12
AmpIRT1
IRT
13
AmpIRT2
IRT
14
AmpIRT3
IRT
15
PhaIRT2
IRT
16
PhaIRT3
IRT
17
PhaIRT4
IRT
18
THD
IRT
19
α
I
R
T
IRT
20
β
I
R
T
IRTFigure 12
Results of the data registration process carried out for features vector obtained for sample impacted by 2J.Figure 13
Results of the data registration process carried out for features vector obtained for sample impacted by 5J.Figure 14
Results of the data registration process carried out for features vector obtained for sample with inclusion A.Figure 15
Results of the data registration process carried out for features vector obtained for sample with inclusion B.The THz signals features allow observing rapid changes in the structure of the examined composite materials. In particular inclusions are very well visible. It is possible to precisely localize rectangular shaped inclusions in case of both samples in distributions of features 4, 6, 7, 8, and 10. In case of samples
I
N
C
_
A during the production process additional defects were unintentionally introduced beside the inclusion. These defected areas can be visible as straight lines ranging from top left corner to bottom right in the distributions of features 1, 3, 8, and 11. For impacted samples the THz features rather allow indicating the localization where the highest damage occurred rather than the range of the defected area.Contrary to THz, the IRT imaging features are characterized by the low dynamics of the heating system. In case of composites, which generally have low value of thermal diffusivity factor, the heat conduction process within material is relatively slow. Therefore, the indications of defects observed at the sample’s surface are, in most cases, blurred. This effect can be noticed for samples with inclusions (it is not possible to identify the shape of the inclusions). On the other hand, the obtained features distributions make it possible to fully assess the range of the impacted area even in the sample impacted with lower energy. Despite the lower ratio between the response to defects and the background values than in case of the THz features, it is possible also to observe the additional defects in included samples
I
N
C
_
A (see feature 16, 18, and 19 distributions).


## 5.3. Multiresolution Decomposition Based Data Fusion Algorithm

The results of both methods represent different spatial dynamics of the feature changes. In case of the THz the inspecting wave is focused on local point of the material, so the response is also from that specific point, while in case of IRT testing, obviously, the heat generated by the halogen lamps cannot be focused at one point but is radiated to the sample’s surface and then is conducted evenly throughout the material volume. Therefore the thermal response in the specific point is disturbed by the surround. This results in lower spatial dynamics of the IRT imaging in comparison to THz one. Therefore, in order to preserve the information content of both rapid and slowly changing responses, a multiresolution decomposition (MRD) based data fusion algorithm was applied to the features vector.The block diagram of the procedure is presented in Figure16. First, in order to reduce the dimensionality of the problem, the knowledge extraction was carried out for the features’ vector. The purpose of this stage is to transform the original features’ vector into a new coordinate system of reduced dimensions level by extracting the most useful information from the original data. In order to carry out the multivariate analysis two algorithms were used: Independent Component Analysis (ICA) and Principal Component Analysis (PCA).Figure 16
Multiresolution decomposition based data fusion algorithm diagram; DWT and DWT−1: wavelet multiresolution decomposition and its inverse.The general idea of the ICA is to extract from multidimensional set of data the independent components under the assumption that the sources of data are statistically independent [19]. This means that each feature value gathered in the vector F is a combination of the original sources data. Then, the ICA analysis allows estimating the original sources from the collected features vector.The objective of PCA is to convert the database of original possibly correlated features into a database having new, statistically independent variables using orthogonal linear transformation [19, 20]. The new parameters called principal components pc are obtained from eigenvectors of covariance matrix of the original database and ordered in accordance with quantity of their variance (described by eigenvalues) in the database. The greater the variance of the variable is, the more significant it is. Typically the first few pc variables express the majority of database variance. Therefore they can be used to represent the whole database without losing a significant amount of information.In both cases the analysis of the two new components (combined featuresc
f
1 and c
f
2) was taken for further processing: two independent and first two principal components representing the greatest variance of the whole set.After the new components extraction MRD was carried out. The general idea of MRD is to represent the source of data with a collection of hierarchical representations of basis function at different resolutions (frequency bands). The original data is decomposed into a set of spatial frequency bandpass representatives, obtained by convolving and subsampling operations. The data can be decomposed using various methods [5]. In this paper the wavelet decomposition WT was utilized. WT is one of the most effective decomposition methods allowing us to obtain a good resolution in both time and frequency domain [21, 22]. After the decomposition the two distributions inputs at given sublevel are fused using the specified fusion rule. Then, the inverse transform DWT−1 is computed and the fused distribution is reconstructed. The whole process can be described by(7)
F
D
F
x
,
y
=
D
W
T
-
1
ϕ
D
W
T
F
C
1
x
,
y
,
D
W
T
F
C
2
x
,
y
,where ϕ is the fusion rule, F
C
1 and F
C
2 are components 1 and 2 obtained after extraction process, and F
D
F is the fused distribution.The result achieved for the proposed algorithm with the utilization of Daubechiesψ
D
10 wavelet function for the decomposition of c
f
1 and c
f
2 and the maximum selection as the fusion rule [22] are presented in Figures 17 and 18, respectively, for ICA and PCA knowledge extraction procedure. The general advantage of the used data fusion algorithm is that no matter which feature would provide the information about the defect it would be indicated in the fused distribution. Analyzing the obtained results one can notice that the fused distributions carry crucial information about the defects and the structure delivered by all features. All defects were indicated in single fused distributions, which can be especially seen in case of samples with inclusions s
I
N
C
_
A and s
I
N
C
_
B. Besides clearly visible rectangular shaped inclusions it is possible to observe additional defect of the structure in the neighborhood of the inclusions (Figures 17, 18(c), and 18(d)). Also in case of impacted samples the achieved results allow distinguishing the areas of different damage level. The applied knowledge extraction procedures allowed compressing information in much less data representations, and the multiresolution decomposition fused both rapid and slowly changing signals. This was useful once again especially in case of samples with inclusions. Inclusions are characterized mainly by the high dynamic signals’ components while additional anomalies in the composite structure being results of faulty production process are indicated in the low dynamic components.Figure 17
Results of the ICA-DWT data fusion algorithm obtained for sample: (a)s
I
M
P
_
2
J, (b) s
I
M
P
_
5
J, (c) s
I
N
C
_
A, and (d) s
I
N
C
_
B.
(a)
(b)
(c)
(d)Figure 18
Results of the PCA-DWT data fusion algorithm obtained for sample: (a)s
I
M
P
_
2
J, (b) s
I
M
P
_
5
J, (c) s
I
N
C
_
A, and (d) s
I
N
C
_
B.
(a)
(b)
(c)
(d)


## 6. Conclusions

Methods of inspection of the composite materials still need to be developed. In comparison to the solid materials such as steel, the structure of composite materials results in complicated response even in no defect cases. Therefore, a need of implementation of the various techniques arises. Two different nondestructive techniques, pulsed THz inspection and active IR thermography, were utilized. Both of these methods are contactless; thus their common implementation is simpler. Proposed techniques have various contrast mechanisms. THz method, as a pure electromagnetic one, is sensitive to electromagnetic parameters like permittivity (or refractive index) changes. Due to the very small power of picosecond pulses utilized as excitation, temperature changes in the evaluated material can be omitted. In case of IRT, from the other hand, the heat transfer within the examined structure is crucial, and observed temperature values are connected directly to the changes of different physical properties, like thermal conductivity or material density. It was shown that different features of THz and IRT inspection methods signals can provide strong indication of different defects that arose in thin basalt fiber reinforced composites. However no single feature allows observing all aspects of material’s structure state. Therefore, the next step in inspection procedures, which was presented in this paper, is to combine the information from various techniques to present the results in common format. It was shown that utilization of methods based on both electromagnetic and thermal parameters enables better detection of various types of defects in basalt fiber reinforced composites. The applied multiresolution data fusion algorithm allowed preserving information representing diverse dynamic ranges of signals. Different defects types reflect different spatial frequency of the obtained parameters distribution. Moreover both inspection techniques are also representing different defects response projections on the acquired signals. Therefore the utilized MRD with support of knowledge extraction procedures can result in an effective evaluation of the examined materials’ structure state.The performance of data fusion algorithms can be assessed both qualitatively and quantitatively [23]. The subjective estimation of quality (qualitative assessment) of achieved results confirmed that the fusion of few inspection methods can be effective in detection of different defects such as resin and fibers degradation or inclusions. The complementary information about the defects details gathered by single inspection methods is fused into the common representation. Nevertheless, quantitative evaluation of the data fusion results is a complex process, because there exists no ground truth reference distribution contacting all information about the structure state.



---

*Source: 1249625-2016-07-03.xml*