3D Non-destructive Imaging of Punctures in Polyethylene Composite Armor by THz Time Domain Spectroscopy

An ultra-high molecular weight polyethylene composite sample totally punctured by a projectile was examined by THz TDS raster scanning method in reflection configuration. The scanning results correctly match the distribution of delaminations inside the sample, which was proven with cross-sectional and frontal views after waterjet cutting. For further analysis, a signal-processing algorithm based on the deconvolution method was developed and the modified reference signal was used to reduce disturbances. The complex refractive index of the sample was determined by transmission TDS technique and was later used for the simulation of pulse propagation by the finite difference time domain method. These simulations verified the correctness of the proposed method and showed its constraints. Using the proposed algorithm, the ambiguous raw THz image was converted into a binary 3D image of the sample, which consists only of two areas: sample—polyethylene and delamination—air. As a result, a clear image of the distribution of delaminations with their spatial extent was obtained which can be used for further comparative analysis. The limitation of the proposed method is that parts of the central area of the puncture cannot be analyzed because tilted layers deflect the incident signal.


Introduction
Nowadays, newly developed composite materials introduce significant improvements in many domains of security including the protection of law enforcement and soldiers. This group of materials contains mainly composites made of the following fibers: carbon, glass, aramid, ultra-high molecular weight polyethylene (UHMWPE), and others [1][2][3][4]. In comparison with classic materials, these composites have high strength, rigidity, and big resistance to environmental conditions, and can be formed in a whole variety of different shapes [2][3][4][5].
The UHMWPE composite, being the topic of this study, consists of densely packed perpendicular layers of polyethylene fibers glued together with polyurethane resin [1][2][3][4][5]. The whole sample forms very hard and lightweight plates about 1-2 cm thick, which can be used as ballistic shields in vehicles, bulletproof jackets, and helmets [2][3][4].
Here, we consider a UHMWPE composite sample manufactured from HB26 material (Dyneema [2]) which was totally punctured by a 7.62-mm projectile from an AK-47 assault rifle (Kalashnikov). The first phase of the impact results in total destruction of the fibers and resin due to shearing of the material [5,6]. Simultaneously, the tip of the projectile is flattened, which causes a noticeable decrease of its effectiveness. Later, the resistance of UHMWPE is totally revealed and the stretching of the composite material can be clearly seen. The analyzed sample acts like a system of nets linked together, which can eventually stop the projectile. Such progressive moderation of the velocity of the projectile results in bulges formed on the surface of the plate and causes internal delaminations occurring not only in the area close to the puncture. The spatial distribution of these delaminations and their position, range, direction, and width can be compared with the results of numerical simulations of the sample-projectile interaction which leads to better understanding of the underlying physical phenomena. All these parameters allow for thorough investigation of this interaction and enable the determination of more details connected with the destruction of individual layers of the composite material and delaminations occurring between them. The detailed study of the underlying physical phenomena, in particular the process of destruction, is crucial for designing materials assuring higher ballistic protection.
The classical method of determining the internal structure of UHMWPE composites is cutting, mostly using a waterjet cutting machine [6]. Usually, the shearing line is crossing the middle of the puncture which facilitates the examination of the lateral structure of the sample. However, cutting is a destructive technique and some features of the structure of the puncture can be displaced or even destroyed. Moreover, the process of cutting heats up the composite material causing melting of fibers and resin and sometimes even clogging of some delaminations, making it difficult or even impossible for visual identification. In reality, only a limited number of such cuts may be carried out; thus, gathering the full 3D information is almost impossible.
X-ray computed tomography is widely used for non-destructive 3D imaging of punctures in composites and offers high resolution and large penetration depths as well as clear visualization of the internal features [7,8]. Its main drawbacks are the high price of the scanners and the use of the ionizing radiation. Composite materials can also be examined using THz radiation in the range between 0.1 and 3 THz which is non-ionizing and can provide contactless and nondestructive testing (NDT). The internal structure of many materials like thicknesses of internal layers, hidden defects, and non-uniformities inside those materials can also be revealed [9][10][11][12][13]. For many applications, THz techniques can give competitive results in comparison with the established methods like X-ray, microwave, ultrasonic, acoustics, or thermographic techniques [12][13][14][15]. The investigated UHMWPE composites have a very small absorption coefficient (around 0.3-1.5 cm −1 over the whole THz spectral range) and low dispersion [16] which makes them perfect test pieces for THz examination techniques. Despite these advantageous properties, research concerning UHMWPE composites is until now only superficially described [4,[17][18][19].
The TDS reflection method was selected from many available THz scanning techniques because it offers high-resolution 3D imaging [9-13, 15, 17-19] through the time-of-flight analysis [9]. It has been noticed that scanning of the skewed sample (by 45°in relation to the scanner axis) improves the contrast ratio of pulses reflected from delaminations. The scanning results presented here correctly match the distribution of delaminations inside the sample, which was shown in the cross-sectional and frontal views-so-called B-and C-scans. A thorough comparison of positions and widths of five delaminations for one side of the sample was carried out. The results were consistent with the real condition of the sample. The only disadvantage is the fact that a part of the central area of the puncture cannot be analyzed by this method.
The raw experimental data can hardly be used for further comparative analysis, because positions and widths of delaminations are not unequivocally determined. Hence, the signal processing algorithm, which bases on the deconvolution, was developed and positively verified by the numerical simulation. Thanks to the algorithm, the raw THz image was converted into a binary 3D image of the sample, with a clear distribution of all delaminations with their spatial extent. To best of our knowledge, it is the first time that such a complicated internal structure of the composite plate hit by a projectile was determined in such detailed manner by means of terahertz radiation.

Characterization of the Sample
The 250×250-mm 2 and 16-mm-thick UHMWPE sample was manufactured from a HB26 laminated composite from Dyneema [2]. The laminate is composed of 17-μm-diameter SK 76 crystalline fibers with the extended molecular chains highly aligned along the fiber axis. Particular fibers are formed into yarns, which are composing a single layer of the material (a ply). The thickness of a single ply is about 60 μm, which means that the sample is composed of around 266 plies. Plies are then positioned in the [0°/90°] arrangement and joined by polyurethane resin forming so-called pre-pregs (Bpre-impregnated^composite fibers), each having the thickness 0.27 mm. Each ply consists of 83 % of fibers [5]. Then, many cut and arranged pre-preg layers were pressed at an elevated temperature, forming the investigated sample.
The sample was placed in a holder mounted perpendicularly to the trajectory of the projectile, and then, it was shot by the 7.62-mm projectile with a steel core using a ballistic barrel. The bullet weighed 7.9 g, and the impact took place with a velocity of 720 m/s according to the Polish Standard PN-V-87000:2011.
The impact point was chosen not in the middle of the sample to enable observing the asymmetry of the delaminations' distribution. The projectile totally punctured the sample. A small bulge is visible around the inlet (Fig. 1a, e), whereas the outlet bulge is significantly bigger (Fig. 1b, e). Metallic tape was used in the THz measurements as markers to indicate the back surface of the sample. Both inlet and outlet surfaces of the sample are slightly curved after being hit by the projectile, and thus, a perfect alignment in the THz scanner is impossible.
On the lateral surfaces of the sample, some delaminations (up to 100 mm long) are clearly visible (Fig. 1c).
For further comparative analysis, the area around the point B, marked in Fig. 1, was chosen and precisely characterized. The lateral surface of the sample (Fig. 1c) indicates widths of five hardly visible delaminations. The magnified area corresponding to the point B (Fig. 1d) shows the same five delaminations marked on the left with roman numerals (I-V). Numbers on the right indicate estimated distances between those delaminations.
After performing the THz scanning, the sample was cut using a waterjet cutting machine along the axis marked with the broken line in Fig. 1a. In the photograph illustrating the lateral surface of the cutting (Fig. 1e), we can distinguish an 8-mm-diameter bullet entry channel (for the first 5 mm of the puncture) and a 10-14-mm-diameter bullet exit channel with a characteristic narrowing shape. The edges of the puncture are jagged and covered with grey matter. The puncture channel is surrounded by an area having a thickness of a few millimeters, where the sample material is soft and fibers are torn, melted, and convoluted. The inhomogeneous and asymmetrical inlet bulge has a diameter of 40 mm and a height of 3 mm. The projectile leaving the sample forms an oval, uniform and symmetric outlet bulge having a diameter of 120 mm and a height of 10 mm.
Delaminations are radially extending from the puncture and have different widths and ranges. Figure 1f illustrates the estimation of widths of chosen delaminations determined using a microscope. Layers from the outlet side were also Bcontracted,^causing their pullback of about 1-3 mm, which can be noticed in Fig. 1e as steps (area around point B). Only the steps indicated as delaminations II and IV are visible, and all others cannot be seen due to the fact that waterjet cutting slightly warms up the sample. This leads to a melting of the resin and sticks together the delaminations so that they are not visible in the cut region. This fact was verified by carefully opening the delamination behind the melted resin with a sharp and thin knife. The OT cutting line corresponds to the notation from Fig. 2

Experimental Setup
The sample was measured by the TPI Imaga unit from TeraView with an XY stage for raster scanning in reflection configuration (Fig. 2). The TDS-based setup uses fiber-coupled emitter and receiver heads mounted in arms. A THz beam with Gaussian shape has linear polarization E parallel to the X axis of the scanner. The beam is focused on the sample by means of a 150mm focal length Tsurupica lens to a spot size of approximately 3.6 mm (at 1 THz). The incident angle and the reflection angle are equal to 7°. The XY stage uses servomotors equipped with an encoder for accurate real-time positioning of the sample. The linear spatial resolution for all motor stage components is 20 μm with the repeatability of approximately 20 μm.
The registered 196-ps-long waveform had 4096 data points and the averaging was 100 times, which lead to a total measurement time of approximately 100 h. Such a high averaging was selected to suppress short time fluctuations of delamination-related pulses (standard deviations of pulse amplitudes) below about 3 %. It provides the appropriate smoothness of the waveforms, which simplifies further signal processing. A signal with FWHM = 0.33 ps reflected from a metallic mirror was used as a reference ( Fig. 3(a)). Small oscillations behind the main pulse (plotted in the inset in Fig. 3(a)) are caused by two main factors: the propagation of the THz radiation through a partially dried atmosphere having a humidity around few percent and small reflections coming from the inside elements of the TDS system. Such a pulse can be observed around 27 ps after the main pulse, but its amplitude is very small in comparison to the main pulse (which is around 15 times larger) and, therefore, it has no significance in further measurements and analysis. The power spectrum of the reference pulse is in the 0.2-3.5-THz range with the maximum dynamic range equal to about 65 dB.
Two signals reflected from the sample in the corner point A (Fig. 2) are illustrated in Fig. 3(b). The sample was mounted at angles of 0°and 45°between the selected axis of the sample and the direction of the electric field vector E of the THz radiation. In both cases, it can be noticed that waveforms start with the high front surface echo at 20 ps and Fig. 2 The experimental setup end with the small back surface echo at about 180 ps. The difference between front and back surface echoes is equal to 163 ps. Then, knowing the thickness of the sample in point A (equal to 16 mm), the group refractive index was determined and is equal to n g =1.529.
UHMWPE fibers show a distinct refractive index anisotropy [16,17], and thus, oscillations in the area of 20-100 ps can be noticed. If the THz field is polarized along the direction of the fibers in a layer (0°orientation), the refractive index is higher by about 0.05 than for the perpendicular case (90°orientation). Therefore, if the radiation is parallel to the direction of fibers of one layer (0°), the modulation of the refractive index is observed along the propagation axis of the THz radiation. Reflections from each interface between successive layers result in oscillations plotted in Fig. 3(b). Theoretically, for an angle of 45°, oscillations should disappear because there is no modulation of the refractive index. In the real case ( Fig. 3(b)), small oscillations are present due to the imperfect production process and the fact that the refractive index of particular layers of fibers may slightly differ and interfaces between layers may not be totally smooth.
In the following measurement, the sample was mounted at an angle of 45°to the polarization of the E field of the illuminating THz radiation, thus to the X axis of the scanner (Fig. 2). Such solution enables minimizing reflections from composite layers and therefore improving the visibility of small defects. Moreover, the biggest delaminations of the composite material ( Fig. 1c) extend perpendicularly to the edge of the sample and, therefore, the 45°orientation (marked with the broken line square in Fig. 2) allows for optimum measurements of the area of interest. The PR and QT diagonals form an axis of the sample parallel to the corresponding edges of the sample and crossing the center of the puncture indicated with the letter O. After the scanning, the sample was cut along the QT section.
During the measurement, the 200-mm-side-length square area was scanned line after line and point after point with a 2-mm step width, which resulted in the 3D matrix of points having the dimensions 101×101×4096. The sample was scanned from the inlet side, because the inlet bulge is causing smaller disturbances compared to the bigger outlet bulge. The THz radiation was focused on the surface of the sample in such a way as to obtain the maximal value of the reflected signal in point A. Fluctuations of the height of the sample outside the inlet bulge area were around ±1 mm which does not influence the value of the reflected signal because of the relatively large depth of field of the used THz setup.

Results
As a result of the scanning, we obtained a 3D matrix with X and Y axes representing position and Z axis representing the optical delay time. For further analysis, the unprocessed waveform was chosen corresponding to point B (Fig. 4). The pulse reflected from the front surface of the sample is clearly seen, as well as the pulse corresponding to the back surface covered with the aluminum foil. The presence of the foil causes the change of the phase of the pulse in comparison with the pulse reflected from the back surface without the foil (Fig. 3).
The first 5 mm of the puncture from the inlet side is called the shearing area. There, delaminations occur only in the area closely surrounding the puncture [5,6]. In the 0-50-ps time range (corresponding to the shearing area), there are many small oscillations being the result of successive composite layers and residual influence of the water vapor from atmosphere. Therefore, it was assumed that the 0-50-ps time range for the whole sample (excluding the 50-mm-diameter central area) will be treated as the area without delaminations, which simplifies further analysis.
In the middle of the waveform, there are five pulses (I-V) corresponding to delaminations illustrated in Fig. 1d. Assuming that widths of these delaminations are negligible in relation to the thickness of the sample (equal to 16 mm) and knowing the time difference between the front and back surface echoes (equal to 163 ps), the distances between delaminations were determined (Fig. 4). The results stay in good agreement with estimations from direct measurements (Fig. 1d).
Between distinct pulses corresponding to delaminations (the inset in Fig. 4), there are many small oscillations. Such oscillations may result from several reasons like, for example, internal reflections inside one delamination and reflections between two delaminations or between the delamination and the front surface of the sample; also, they can be caused by a not ideal incident pulse shape and by the residual influence of the water vapor. Another reason for such oscillations may be the existence of thinner delaminations having the width of few micrometers. The THz setup response for such structure has the form of pulses with small amplitudes; however, macroscopic verification of this fact is difficult. It is assumed that such oscillations are not noise due to the fact that noise was minimized by averaging. Moreover, two cross-sectional views present measurement results in 2D space: perpendicular to the surface of the sample (B-scan) and parallel (C-scan). In the discussed case, more interesting are phenomena proceeding in the sample along axes oriented at an angle of 45°to the X and Y axes. Therefore, B-scans (Figs. 5 and 6) are plotted as a cutting line along P-O-R and Q-O-T lines. The illustrated images are artificially saturated to better visualize pulses reflected from the inside details of the structure. Moreover, images are extended along the time axis in relation to the real size of the sample (Fig. 1e).
Lines of the front surface echo for 15 ps are not horizontal but curved because the sample was deformed by the interaction with the projectile. In the area of the inlet bulge, the surface of the sample is not perpendicular to the incident radiation which explains the low signal level and signal deformations. The back surface echo is quite significant in both B-scans (in areas outside the outlet bulge), because it is the reflection from the marking aluminum foil.
The THz technique is perfectly detecting delaminations being out of the central area of the puncture. The position and range of delaminations may be compared to the real sample defects. In Fig. 5, the area corresponding to point B is marked with a broken line. Ranges (lengths) of all delaminations (I-V) can be clearly defined. In the lateral cross section (Fig. 1e), only two of the delaminations can be observed (II and IV) because the others are not visible after waterjet cutting and their existence may be proven only by mechanically removing the melted resin with a knife. In the same way, delaminations localized on the other lateral edges of the sample (P, Q, R) were positively verified.
The approximate positions of borders of the puncture are marked in Fig. 5 with a broken line and were determined on the basis of Fig. 1e. The lack of the reflected signal defines the 8mm-diameter inlet channel area for the first 5 mm of the puncture (as in reality). Next, the puncture is broadening and forms a 10-14-mm-diameter hole. Delaminations determined by THz scanning in the central area of the sample are deflected as an effect of the projectile impact, but they do not approach directly the real border of the puncture (broken line). After analyzing both B-scans (Figs. 5 and 6), it can be concluded that delaminations are not arranged symmetrically around the center of the sample; their amount and depth vary for different directions. It may be noticed that more visible delaminations are oriented along the P direction, which is the closest one to the sample edge. For the R direction, which is corresponding to the biggest distance between the puncture and the edge, delaminations are less visible.
Complementary information about the horizontal range of delaminations is taken from Cscans plotted in Fig. 7 in a peak-to-peak slice mode. Time intervals (slices creating 2D images) corresponding to the best visibility of delaminations (I-V) along the T direction were chosen from the entire measurement range.    (Fig. 1a), which is not uniform and is inclined. Scratches along the line O-P corresponds to the waviness of the surface. The last Cscan (Fig. 7h) plots radial stripes determining the reflection from the aluminum foil and the back surface of the sample. Sizes of delaminations and puncture in horizontal direction can be investigated on C-scans illustrating inner layers (Fig. 7b-g). The shearing area consists of small delaminations diverging from the puncture area and can be seen for the time interval 27-50 ps (Fig. 7b). Figure 7c, d plots symmetrically oriented delaminations radially diverging from the center along the fibers' direction and having the shape of regular stripes 15 to 40 mm wide. Figure 7e, f illustrates broadening of striped delaminations which start to cover a larger area of the sample. Finally, striped delaminations become very weak and delaminations that in other intervals were circular now become square-like (Fig. 7g). Generally, it can be noticed that Fig. 7 C-scans in peak-to-peak mode for eight time intervals corresponding to eight slices. All dimensions are given in millimeters. The image is saturated to emphasize delaminations delaminations in the P-Q-T area are much more visible due to the asymmetrical location of the puncture in the sample. C-scans ( Fig. 7c-g) contain estimated widths of delaminations on the edge of the sample, which stays in good agreement with measurements carried out by the use of a ruler (Fig. 1c). Widths obtained from THz measurements are slightly bigger because determining the precise position of the end of each delamination is not so accurate.
Interaction of the high-energy projectile with UHMWPE with a relatively low melting temperature (150°C [5]) can locally change its propagation properties, especially the refractive index. In the 3-5-mm-wide area surrounding the puncture, composite fibers are melted or frayed after the interaction with the projectile, which scatters THz radiation and significantly decreases the reflection. Moreover, in this area, delaminations are not perpendicular to the incident radiation, so the reflected pulse does not reach the detector. Therefore, the discussed technique is not capable of proper determination and investigation of the area limited by the dotted line in the region of the outlet bulge. However, the horizontal lines in Figs. 5 and 6 outside the puncture region, representing delaminations, are rather flat, which suggests that the refractive index is not changed by the projectile and the optical depth is converted into thickness with good accuracy.
Summing up, after precise analysis of the selected area of the sample, it was proven that the proposed method may be applied to qualitative and rough determination of positions and ranges of delaminations in the HB26 sample. The investigation of ply separation is possible out of the dead zone corresponding to the outlet bulge.

Simulations and Signal Processing
The analysis presented so far qualitatively depicts the clear image of the distribution of all delaminations inside the sample. However, such a result cannot be used for further analysis, like the comparison of numerical simulations with experimental results. The main problem is the fact that the registered sample is not binary, which would mean that it consists only of two areas: polyethylene B1^and air (delamination) B0.^Therefore, based on the deconvolution method, a very simple and robust time signal processing technique was elaborated. As a result, the transformation of the registered THz image to the 3D binary image is possible.

Transmission Characteristic of the Sample
A standard time domain spectroscopy system Spectra 3000 from TeraView was used to investigate a 3.8-mm-thick delamination-free HB26 sample prepared in the same manner as the sample described above. The determination of the complex refractive index N 1 = n + iκ, where n is the refractive index and κ is the extinction coefficient, using the procedure described in [20], was very important for further simulations. Both parameters are plotted in Fig. 8. Due to many weak internal reflections inside the multilayer sample, the registered characteristics have many small oscillations and a feature at about 0.85 THz. The absorption peak of the polyethylene at about 2.2 THz can also be clearly identified. One can notice (Fig. 8a) that n is decreasing almost linearly in the 1.533-1.526 range and the group refractive index n g = 1.529 (determined in Section 3) is consistent with these calculations. The extinction coefficient κ is considered constant and is equal to 2.5 10 −3 .

Finite Difference Time Domain Simulation
The 16-mm-thick planar-parallel dielectric structure has the complex refractive index (N 1 ) determined in Section 5.1 and contains five air delaminations (with the refractive index N 0 ). These delaminations (indicated with numbers I-V) have the widths 300, 75, 25, 75, and 300 μm, starting at distances from the beginning of the sample: 5, 7, 10, 11, and 13 mm, respectively. The positions and widths of the delaminations were adjusted to be similar with the presented measured waveform (Figs. 11 and 12).
The propagation of the THz pulse in the structure was simulated in 2D by means of the software OptiFDTD 5.0 from OptiWave. It uses the finite difference time domain (FDTD) method with a perfect magnetic conductor and perfectly matched layers boundary conditions. Proper discretization of time (2.2 fs) and space (1 μm) provides good convergence. The simulation lasted about 10 min using a tabletop PC. The structure was illuminated perpendicularly by the pulse plotted in Fig. 9a which was determined by replacing the structure with a perfect conductor. The shape of the pulse is very similar to the shape of the pulse from the experiment (Fig. 2), and it contains two side lobes. The FWHM is equal to 0.33 ps. After carrying out numerical simulations, the number of waveform samples was decreased to a value of 4096, which corresponds to the number of samples in the experiment. A small noise was added to the resulting waveform which cannot be seen in Fig. 9b.
The reflected waveform (Fig. 9b) reveals a series of pulses corresponding to the front and back surfaces of the structure and to inside delaminations. Dispersion can be observed as  Bstretching^pulses in time, and therefore, the FWHM of the pulse reflected from the back surface of the structure is about 1.1 ps.
Each delamination forms a double interface between the media-PE/air and air/PE-and therefore, in accordance with Fresnel's equations, we expect two time-shifted and opposite phase reflections from each delamination. Nevertheless, the width of the incident pulse is quite large in comparison with the propagation time through examined delaminations and, therefore, both pulses overlap in different manners. For delamination I (300 μm), two separated pulses reflected from both interfaces may be resolved. In the case of delamination II (75 μm), these pulses cannot be separated, but amplitudes of constituent sub-pulses are relatively large. Delamination III (25 μm) is similar to the second one, but sub-pulse amplitudes are significantly smaller due to the fact that for such width of the delamination, consecutive pulses add destructively. Dispersion causes pulses 4 and 5 to become Bstretched^copies of pulse 2 and 1, respectively.
Moreover, many additional internal reflections can be observed in this structure, which is called the Fabry-Perot effect (FP). The reflection inside delaminations (FP1) occurs directly after pulses coming from corresponding delaminations. Reflections between interfaces of different delaminations and between these surfaces and the beginning of the sample (FP2) are separated in time. It may happen that the reflection FP2 will overlap with the pulse reflected from the delamination, which will complicate further analysis.

Deconvolution Method
The analysis carried out here shows that the resulting waveform has a quite complicated form as time pulses partially overlap and it is difficult to unambiguously determine the beginning and the end of the delamination. Knowing the reference waveform illuminating the sample (h) and the one reflected from the sample, namely the sample signal (g), the sample pulse response (f) may be determined by means of the 1D deconvolution method and then interfaces between media [21]. After testing various algorithms [21][22][23][24], we decided to use inverse filtering coupled with a band pass double Gaussian filter [22]. The sample pulse response function (f) can be calculated from the formula [22]: where FFT is the fast Fourier transform. A common approach is to use a band pass filter to suppress the amplified noise effects caused by the division operation in Eq. (1). The filter is formed by combining two Gaussian functions with pulse width determined by high-frequency (HF) and low-frequency (LF) coefficients [22]: where t is the time. Since all signals consist of N = 4096 discrete points, e.g., g = g (i), i = 1, 2,… N, for further calculations, the time domain can be substituted with equally spaced points (1, 2,… N) and, therefore, values of HF and LF are now dimensionless numbers. Values of HF and LF can be set manually. Later, it is assumed that LF = 4096. Figure 9b (blue waveform) plots the signal which is deconvolved and filtered by the band pass double Gaussian filter for HF=5. The transition from the medium with a lower refractive index (N 0 ) to the medium with a higher refractive index (N 1 ), which corresponds to the air/PE interface (either beginning of the sample or second interface in a delamination), results in the Gaussian-like positive reflected pulse. The inverse situation results in the negative pulse, which corresponds to the end of the sample or the PE/air interface.
Pulse peaks are determining positions of corresponding borders of delamination, and the time difference (τ) between these pulses defines the width of the delamination (d) according to the formula: where c is the speed of light. Thus, for delamination I, the time difference is equal to τ 1 =2 ps which corresponds exactly to the 300-μm width. But for delamination II, the time difference τ 2 =0.95 ps corresponds to the width 142.5 μm, which is two times bigger than the real width of this delamination. For delaminations III, IV, and V, time differences are 0.8, 1, and 2 ps, respectively, and corresponding widths are 120, 150, and 300 μm, respectively, while real widths are 25, 75, and 300 μm, respectively. In other simulations, it was estimated that the limit of the proper determination of the width of the delamination is about 150 μm. Obviously, this is a significant limitation of the used method; however, taking into account the practical point of view, it is not so important, because the widths of delaminations may fluctuate depending on local material parameters of the sample and defining the width as being in the 50-150-μm range seems to be sufficient.
Simulations carried out here indicate that pulses reflected from about 25-μm-wide delaminations will have relatively small amplitudes which may be difficult to distinguish from other disturbances existing in real measurements.
Deconvolution carried out for the measured signals (Section 6) was strongly disturbed by slowly changing oscillations. A few signal-filtering techniques were tested, but the best results can be achieved by using the deconvolution method for Bcut^truncate reference waveform (h ′). Such waveform has all values after the main pulse ( Fig. 9) equal to 0 (marked with an arrow). The signal after deconvolution (f′) using the h′ reference is plotted with the blue line in Fig. 9b and has a very similar shape to the f signal. The f′ signal has a negligible lower amplitude than f, and the time difference between peaks is slightly smaller. Nevertheless, taking into account earlier described limitations, these differences seem to be insignificant.

Derivatives Analysis
On the basis of the waveform shape after applying the deconvolution method, the last step in the analysis is related to the automatic determination of minima and maxima corresponding to delaminations. Taking into account the waveform registered in real conditions, the difficulty lies in the differentiation between slowly changing disturbance signals and fast-changing components corresponding to delaminations. Therefore, the method using the time derivative of the deconvolved signal was proposed and denoted with the letter q. Figure 10 plots signals after deconvolution (f 1 ′ and f 2 ′) and their derivatives (q 1 and q 2 ) for delaminations I-III computed for HF 1 =5 and HF 2 =20. Amplitudes of signals f and q were adjusted to the scale of the image.
The way of selecting the optimum value of the HF coefficient is simply choosing the smallest possible value which does not introduce additional distortions to the signal f [22]. Moreover, the smaller the value of HF is, the larger is the visibility of pulses (Fig. 10) and the determination of widths of delaminations is more consistent with reality. Therefore, instead of looking for one value for the entire spectrum of waveforms (~10,000), only two filter values will be applied: the smaller one, HF 1 =5, which assures large contrast value of pulses but may cause distortions, and the bigger one, HF 2 =20, which allows for compensation of distortions caused by HF 1 .
The interesting part of the waveform f 1 ′ between the peaks p 2 and r 2 corresponds to the positive part of the waveform q 1 (Fig. 10b). The conventional zero value for the q 1 and q 2 curves is marked with the broken line in Fig. 10. Due to the existing noise and Fabry-Perot reflections, the constant threshold level th (marked in Fig. 10 with a dashed line) was introduced. In the simplest case, the proposed algorithm at first is seeking for areas of the q 1 signal having a value bigger than the threshold value and on this basis is determining the preliminary duration time corresponding to the delamination (t 2 ). Then it seeks the zero values of the q 1 signal, determining the final duration T 2 .
Due to the low value of HF 1 , waveform q 1 for delamination I has a central dip (Fig. 10a) and, therefore, the algorithm determines two preliminary durations (t 11 and t 12 ). If the time difference between t 11 and t 12 is lower than the threshold value (usually 1-2 ps), it is suspected that this can be a single delamination and the non-distorted q 2 waveform is analyzed. It can be shown (see also Fig. 10a) that a positive value of q 2 in the range between t 11 and t 12 is a sufficient condition that indicates a single delamination. The proper preliminary duration t 1 contains t 11 , t 12 , and the interval between them, so the final duration T 1 is determined referring to zero values of q 1 as mentioned above. Figure 10b, c illustrates that a properly defined threshold rejects small oscillations (FP1 and FP2) and passes the pulse coming from delamination III. The duration T 3 was determined in the same manner as that of T 2 . The widths of these delaminations are the same as determined in Section 5.3.

Signal Processing and Results
This section contains test results of the developed algorithm for experimental data processing and 2D and 3D visualizations. The wavy signal after deconvolution (f in Fig. 11b) was (1)-(2) to the measured signal reflected from point B (Fig. 5), when HF 1 =5. The original reference signal (h) is plotted in Fig. 3, and its initial part is also shown in Fig. 11a. It is worth noticing that pulses corresponding to delaminations I, II, and IV are clearly visible, while those connected with delaminations III and V could be difficult to detect, especially automatically.
Some oscillations (e.g., at 70, 120, 170 ps) having high amplitude values may be also difficult to distinguish from comparable pulses corresponding to delaminations (I, II, and IV). The oscillations result from the fact that the deconvolution method requires the uniform illumination (exactly the same pulse shape) of the reference mirror and all points of the sample. All these points should be equally distant from the THz emitter, which means they must be placed exactly in the same place along the propagation axis (to have the same value of time delay) [23]. These requirements are not satisfied because the sample is not ideally flat and time fluctuations of the incident pulse occur during the measurement. We must be aware of the fact that the pulse generated in the TDS setup is not time-invariant during long measurements.
After testing few signal-filtering methods, it occurred that the best results give the deconvolution method applied to the truncated reference waveform for which all values after the main pulse are equal to 0 (marked in blue in Fig. 11a). The signal after the deconvolution f′ is now much more flat with clear peaks, and it resembles the waveforms obtained in the simulations (Fig. 9). The visibility and differentiation of pulses have significantly increased, and therefore, the waveform h′ was chosen for further analysis. Fig. 11 The reference signal (a) and the signal after deconvolution (b) Fig. 12 Signals after deconvolution and their derivatives Figure 12 plots two waveforms after applying the deconvolution method (f 1 ′ and f 2 ′) obtained for the reference signal h′ for two filter values HF 1 =5 and HF 2 =20. Based on the waveforms f 1 ′ and f 2 ′ the automatic detection of pulses corresponding to all delaminations is still difficult due to different amplitudes of these pulses. Using two threshold values (for example th 1 and th 2 ) enables finding minima and maxima corresponding to real delaminations (like p 4 and r 4 ) and those coming from different features of the analyzed medium (like p 3 and r 3 ).
Derivatives (q 1 and q 2 ) of signals f 1 ′ and f 2 ′ considerably simplify the process of unambiguous detection of delaminations under inspection. One threshold value (TH) enables separation of pulses coming from delaminations (I-V) from other oscillations (noise, Fabry-Perot reflections, and other disturbances).
The algorithm (derivative analysis) described in Section 5.4 enables the estimation of preliminary time durations (t 11 , t 12 , t 2 ), and on this basis, the final durations of pulses for particular delaminations (T 1 -T 5 ) are determined. Then, by using Eq. 3, we can calculate the widths of all these delaminations (D 1 -D 5 ) which are presented in Table 1.
The developed algorithm was verified for many measuring points and then applied for all points of the sample, which enabled transforming the analog image into the 3D binary digital one. For 10,000 measuring points, the filtering by deconvolution, derivative analysis, and transformation to the binary image lasted about 50 min using a tabletop PC. Processing of the raw data into the binary image is crucial for further automatic detection of delaminations and becomes possible by properly adjusted threshold values. Figure 13 shows the binary B-scan (bB-scan) corresponding to the one in Fig. 5, carried out along the Q-T axis. The black color indicates delaminations and the sample surrounding (air), the white color corresponds to the UHMWPE material of the sample, while the area of the puncture and the Bdead^zone under the outlet bulge are marked with grey colors. Easy-to-determine front and back surface echoes enable determining the beginning and the end of the sample. In Fig. 13, the Y axis was scaled and now it denotes the depth.
Broad delaminations with ambiguously indicated borders illustrated in the B-scan in Fig. 5 are now replaced with the borders having sharp-edge interfaces with delaminations' widths close to their real values in the bB-scan. The comparison of the photograph of the sample from Fig. 1e with the bB-scan from Fig. 13 shows their good agreement, taking into account both the position and the width of each delamination.
Here, in comparison with the B-scan, all ambiguous areas with small amplitudes (below the TH threshold value) were removed. These areas correspond to Fabry-Perot reflections, noise, and disturbances coming from the deconvolution. Some delaminations are marked as discontinuous lines, like for example 50-120 mm width distance for delamination V. Probably it results from the fact that this delamination is too narrow in these regions; thus, the amplitude of the signal decreases and the q signal has the value below the threshold level TH. This conversion process may also remove the information about small-width delaminations (<20 μm) from the image, but their significance for further analysis seems to be not essential.
The fact that the determination of the position and width of delaminations in the area under the outlet bulge is impossible seems to be a real constraint of this method. Therefore, in  Fig. 13, borders of the outlet bulge are marked with grey color and then delaminations that bend towards the back of the sample were extended and extrapolated (red lines), which is some approximation of their real distribution.
The main advantage of THz scanning methods is the possibility of 3D visualization of delaminations. Figure 14 contains the distribution of delaminations I, II, IV, and V (for clarity, delamination III has been removed). In comparison to corresponding C-scans from Fig. 7, additional data about the width of each delamination and its real shape can be obtained.

Discussion and Conclusions
THz TDS raster scanning in reflection configuration was used for the investigation of an UHMWPE composite sample totally punctured by a projectile. The sample, the experimental setup, and the measurement procedure were precisely described. The internal structure of the sample was determined experimentally. It was found that the contrast ratio of pulses reflected from delaminations can be improved by scanning the sample with an angle of 45°of the electric field vector in relation to the scanner axis. The qualitative analysis revealed that the scanning results using the raw measurement data correctly match the distribution of delaminations inside the sample, which was shown in the B-and C-scans. Five delaminations on a Fig. 13 The binary B-scan along the Q-T axis, where the white color corresponds to the composite material, the black color to the air (delamination or surrounding), and the grey color to the Bdead^zone under the outlet bulge A comparison of measured and simulated waveforms cannot be carried out for raw data due to the fact that positions and widths of delaminations are not precisely known. To suppress the influence of additional defects on the waveform corresponding to delaminations of the analyzed structure, a deconvolution method with an additional derivative analysis was used. Using a modified reference signal, a reduction of disturbances becomes possible. The numerical simulations confirm the correctness of the proposed algorithm and show the limits of its applicability. The complex refractive index of the sample was determined by THz TDS technique in transmission mode and applied for the determination of pulse propagation using the finite difference time domain (FDTD) method. Finally, the binary 3D image of the sample clearly resolves all delaminations.
The elaborated signal processing method is simple and robust and should work properly for all 10,000 waveforms. On the other hand, the proposed algorithm does not always return proper results, because the spectra of individual signals may be quite different and difficult to predict caused by noise, time instability of the setup, internal reflections, and other perturbations. However, the use of the binary image does not require very precise determination of widths of delaminations, as their features are related to local material properties and may be characterized with some inaccuracy. The described signal processing method gives quite accurate positions of particular delaminations and a little bit less accurate, but still sufficient for further analysis, thicknesses of those delaminations. The small error in determining the position may result from local refractive index fluctuations caused either by sample-projectile interaction (especially in the area of the puncture) or by manufacturing inhomogeneity of the material. Moreover, the group refractive index used for the processing was calculated for the thinner sample at a single point and also can fluctuate.
Averaging over 100 individual measurements increased the time of the total measurement of the sample to around 100 h but provided small noise. This time can be significantly reduced using fast TDS scanners with a better S/N ratio and possible filtering. The 3.6-mm spot diameter together with the method of oversampling (the distance between measuring points is 2 mm) ensures a good lateral resolution. An improvement of the resolution can be expected applying a collinear configuration with 0°incidence angle, in which the focal length can be shortened leading to a decrease in spot diameter. The 0°i ncidence angle is optimum for the discussed multilayer structures; however, in such a setup, the dynamic range is inherently worse by 6 dB. The biggest drawback of TDS setups working in reflection mode is the limited possibility to scan surfaces oriented at an angle other than 90°to the incident radiation (Bdead^zone), so that the reflected radiation cannot be detected. A determination of the incident angle with respect to the surface and an automatic alignment system may solve this problem.
In summary, the results obtained in this investigation closely reflect the real conditions of our sample, as far as the real 3D distribution of all delaminations is known. The TDS THz technique seems to be a good choice for 3D non-destructive, non-ionizing, and contactless scanning of delaminations in UHMWPE composites. Clear images of the distribution of all delaminations including their extension and widths were obtained which can be used for further analysis of particular properties of the composite material under investigation.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.