# **Advanced Sensing, Fault Diagnostics, and Structural Health Management**

Edited by

Yongbo Li, Bing Li, Jinchen Ji and Hamed Kalhori Printed Edition of the Special Issue Published in *Sensors*

www.mdpi.com/journal/sensors

## **Advanced Sensing, Fault Diagnostics, and Structural Health Management**

## **Advanced Sensing, Fault Diagnostics, and Structural Health Management**

Editors

**Yongbo Li Bing Li Jinchen Ji Hamed Kalhori**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editors* Yongbo Li Northwestern Polytechnical University China

Jinchen Ji University of Technology Australia

Bing Li Northwestern Polytechnical University China

Hamed Kalhori The University of Technology Sydney Australia

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Sensors* (ISSN 1424-8220) (available at: https://www.mdpi.com/journal/sensors/special issues/ ASFDSHM).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-6181-3 (Hbk) ISBN 978-3-0365-6182-0 (PDF)**

Cover image courtesy of unsplash

© 2023 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**


## *Editorial* **Advanced Sensing, Fault Diagnostics, and Structural Health Management**

**Yongbo Li 1,\*, Bing Li 1, Jinchen Ji <sup>2</sup> and Hamed Kalhori <sup>3</sup>**


Advanced sensing, fault diagnosis, and structural health management are important parts of the maintenance strategy of modern industries. With the advancement of science and technology, modern structural and mechanical systems are becoming more and more complex. Due to the continuous nature of operation and utilization, modern systems are heavily susceptible to fault. Hence, the operational reliability and safety of the systems can be greatly enhanced by using the multifaced strategy of designing novel sensing technologies and advanced intelligent algorithms, and constructing modern data acquisition systems and structural health monitoring techniques. As a result, this research domain is receiving a significant amount of attention from researchers in recent years. Furthermore, the research findings have been successfully applied in a wide range of fields such as aerospace, manufacturing, transportation and processes.

This Special Issue of *Sensors* aims to collect the latest research results and developments encompassing all the areas of advanced sensor design, fault diagnosis, and structural health management. This collection contains 10 papers that represent state-of-the-art technology developed for reliability engineering.

Addressing the challenges of maintenance at high altitudes, Guo et al. [1] designed a novel intelligent technique for detecting the rust on transmission fitting lines with the help of unmanned aerial vehicles (UAVs). A novel convolutional neural network (CNN), namely, the R-CNN algorithm, was proposed to extract rich information about the transmission line from the UAV images. Secondly, a feature enhancement technique was added after the pooling layer of the region of interest (ROI) to enhance the feature representations of the regions that have real fittings.

Aiming to improve the interpretability aeromagnetic data by compensating for the onboard electronic interference (OBE), a data-driven OBE interference compensation method was proposed by Wang et al. [2]. Unlike existing linear algorithms, the proposed method can reduce OBE interference without relying on any reference sensors. A long short-term memory (LSTM) network was combined with wavelet decomposition for detecting and predicting OBE interference and, subsequently, local variations of the magnetic field were estimated to remove the drift of the interference.

A novel crack monitoring technique for engineered cementitious composite (ECC) beams reinforced with hybrid bars using piezoceramic-based smart aggregates was designed by Qian et al. [3]. The designed ECC bars were fabricated and tested under the influence of cyclic loading. Two smart aggregate (SA)-based active sensing methods pasted onto both ends of the beams were used as the actuator and sensor to monitor the growth of the damage. Furthermore, a novel self-repairing index for monitoring the self-repairing capacity was proposed.

Wang et al. [4] proposed an evaluation method using hierarchical analysis based on a combination of expert industry opinions for XLPE cables utilized in coal mines. After

**Citation:** Li, Y.; Li, B.; Ji, J.; Kalhori, H. Advanced Sensing, Fault Diagnostics, and Structural Health Management. *Sensors* **2022**, *22*, 9087. https:// doi.org/10.3390/s22239087

Received: 15 November 2022 Accepted: 16 November 2022 Published: 23 November 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

<sup>1</sup> School of Aeronautics, Northwestern Polytechnical University, Xi'an 710072, China

<sup>2</sup> School of Mechanical and Mechatronic Engineering, Faculty of Engineering and IT, University of Technology, Sydney, P.O. Box 123, Broadway, NSW 2007, Australia

classifying the mining cable health status according to the degree of severity, a cloud model theory was utilized to transform the standard status level into a visualized status space. The membership degrees of each quantitative and qualitative index was calculated. An improved analytic hierarchy process (AHP) was utilized to calculate the weight of each indicator in the indicator layer. The cable health status was judged after fusing the indicator membership and the weight via the D-S evidence theory.

A novel technique for impact calibration of a wide-range triaxial force transducer was proposed by Wang et al. [5] using the Hopkinson bar technique. In the proposed technique, the reference input force for the transducer was generated by the Hopkinson bar. Furthermore, different from the existing methods, the transverse sensitiveness of the triaxial transducer was given importance in the proposed method. The calibration results expressed by sensitivity coefficients were linearly fit by using the least square method (LSM) in a sensitivity matrix.

Targeting the limitations of the conventional velocity stress-dependent acoustoelastic effect for short bolts, a stress-dependent attenuation estimation model was developed by Fu et al. [6]. The effect of axial stress on ultrasonic scattering attenuation was investigated by calculating the change in the energy attenuation coefficient of ultrasonic echoes after the application of an axial preload. Additionally, to overcome the challenge of obtaining a frequency-dependent attenuation coefficient, the bandwidth of the measured echoes was divided into several frequency bands to select the frequency band sensitive to the axial stress changes. Under 20-step axial preloads, the final estimation model between the axial stress and energy attenuation coefficient in the frequency band was established.

Xie et al. [7] studied the combined effect of corrosion and crack defects on the growth of cracks in pipelines. The interaction effect of corrosion on the fatigue crack was obtained by studying the stress intensity factor (SIF) interaction impact ratio after calculating the SIF by using finite element models. One direct and one indirect approach based on extreme gradient boosting (XGboost) had been used to predict the SIF interaction impact ratio. Finally, a crack propagation model was designed based on the XGboost models, and the Paris law and corrosion growth model was proposed for pipelines with interacting crack and corrosion defects.

Aiming to solve the limitations of prognosis operation for Francis turbine units (FTUs) in practical engineering environments, an ensemble prognostic method under variable operating conditions was proposed by Duan et al. [8]. After constructing the running data set with the help of values of the water head, active power, and vibration amplitude of the top cover, a density-based spatial clustering of applications with noise (DBSCAN) was introduced for filtering the raw data from outliers and singularities. From the cleaned raw data, a healthy state model was constructed from a Gaussian mixture model, and subsequently, a performance degradation indicator was calculated by using the negative log-likelihood. Finally, based on the designed indicator, a multiobjective prediction model was proposed based on the non-dominated sorting genetic algorithm and Gaussian process regression.

Zhang et al. [9] analysed the nonlinear dynamic behaviour of tethered satellite formation systems in tethered space systems (TSSs) based on a simplified rigid-rod tether model. Two new stability control laws were proposed based on the tether release rate and tether tension for controlling the variation of tether length. In addition, based on the Floquet theory proposed in 1868, the periodic stability of the post-deployment time-varying control system was analysed.

Focusing on the improvement of the sampling performance of a wireless sensor network (WSN) in real-time situations, Ha et al. [10] proposed an optimal remote monitoring system platform for SHM, which is based on the pulsed eddy current (PEC). The proposed method was utilized for measuring the corrosion of a steel-framed structure. A new circuit was designed to delay the PEC response to tackle the fast-varying output signal in an actual scenario for an efficient sampling performance.

All the Editors would like to extend their gratitude to the authors for their significant and noteworthy contributions. We also give our sincere thanks to all the anonymous

reviewers who took their time and provided feedback to the authors for the improvement of their research. Additionally, we want to thank the publishers and all the members of the *Sensors* board for providing this opportunity to present these works.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **A Robust Faster R-CNN Model with Feature Enhancement for Rust Detection of Transmission Line Fitting**

**Zhimin Guo, Yangyang Tian \* and Wandeng Mao**

Electric Power Research Institute, State Grid Henan Electric Power Company, Zhengzhou 450007, China **\*** Correspondence: yytian.mail@gmail.com; Tel.: +86-150-3730-1821

**Abstract:** Rust of transmission line fittings is a major hidden risk to transmission safety. Since the fittings located at high altitude are inconvenient to detect and maintain, machine vision techniques have been introduced to realize the intelligent rust detection with the help of unmanned aerial vehicles (UAV). Due to the small size of fittings and disturbance of complex environmental background, however, there are often cases of missing detection and false detection. To improve the detection reliability and robustness, this paper proposes a new robust Faster R-CNN model with feature enhancement mechanism for the rust detection of transmission line fitting. Different from current methods that improve feature representation in front end, this paper adopts an idea of back-end feature enhancement. First, the residual network ResNet-101 is introduced as the backbone network to extract rich discriminative information from the UAV images. Second, a new feature enhancement mechanism is added after the region of interest (ROI) pooling layer. Through calculating the similarity between each region proposal and the others, the feature weights of the region proposals containing target object can be enhanced via the overlaying of the object's representation. The weight of the disturbance terms can then be relatively reduced. Empirical evaluation is conducted on some real-world UAV monitoring images. The comparative results demonstrate the effectiveness of the proposed model in terms of detection precision and recall rate, with the average precision of rust detection 97.07%, indicating that the proposed method can provide an reliable and robust solution for the rust detection.

**Keywords:** rust detection; transmission lines fitting; object recognition; faster R-CNN; transmission safety

#### **1. Introduction**

Due to long-term exposure to the wild environment, transmission line fittings are prone to defects such as aging, damage and rust, resulting in heavy risk to the transmission safety. It is significantly important to detect and deal with the rust of transmission line fitting in a timely manner. Presently, unmanned aerial vehicle (UAV) inspection has replaced labor routing inspection in many scenarios due to some merits such as no terrain limitation, fast speed, high efficiency, low labor costs, strong safety and so on. In the UAV inspection mode, however, UAVs generally collect monitoring data for the artificial check, which is with low efficiency. Machine vision with artificial intelligence techniques is currently becoming a promising tool to analyze the UAV monitoring data, and has shown prevailing performance compared artificial check. It is of great significance to develop a robust and accurate rust detection method for transmission line fittings.

From the theoretical perspective of machine vision, the rust detection problem can be viewed as the problem of object detection [1]. With the rapid development of convolutional neural network (CNN), deep learning techniques have become a promising tool in object detection [2]. In summary, these techniques can be divided into two strategies: one-stage detection and two-stage detection. The one-stage algorithm, such as YOLO [3], SSD [4,5], RetinaNet [6], uses a unified deep neural network (e.g., CNN) for feature extraction, target

**Citation:** Guo, Z.; Tian, Y.; Mao, W. A Robust Faster R-CNN Model with Feature Enhancement for Rust Detection of Transmission Line Fitting. *Sensors* **2022**, *22*, 7961. https://doi.org/10.3390/s22207961

Academic Editor: Manuel Pineda-Sanchez

Received: 13 September 2022 Accepted: 10 October 2022 Published: 19 October 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

classification and bounding box regression, achieving end-to-end object detection. It has a faster detection speed and relatively lower detection accuracy. The two-stage algorithms, mainly the variants of R-CNN, i.e., R-CNN [7], Fast R-CNN [8], Faster R-CNN [9] and Mask R-CNN [10], adopt a classical sliding window mechanism to extract interested region and then carry out classification with the features of the regions. In these algorithms, Faster R-CNN is on par with, or even outperforms, the other algorithms in terms of detection accuracy. Nevertheless, the classical Faster R-CNN still has some limitations in the detection of small-size objects, especially under complex background. Many studies have been devoted to overcoming the limitations. For instance, Cui et al. [11] adopted a feature pyramid network in Faster R-CNN with attention module. By highlighting the saliency of object's features, the detection accuracy can be improved. Lim et al. [12] introduced a residual attention mechanism to obtain rich information of small-size objects. Aside from considering feature representation, Xue et al. [13] also introduced coordinate attention mechanism into Faster R-CNN for incorporating the location information that is believed helpful to the channel information. Hong et al. [14] designed a quartile attention mechanism that uses four branches to capture internal and cross-latitude interactions between channels and spatial locations, making better use of contextual information.

These studies can improve the detection robustness under challenging environments by extracting rich semantic information. According to our empirical study, however, these methods do not work well in the rust detection of transmission line fittings. The reason is that the rust detection of transmission line fittings has some special challenges. In most actual applications, the transmission line is long and widely distributed, leading to complex background for the detection. Too many disturbance items such as tree, car, village, house, etc., exist and raise false detection. Moreover, UAV graph usually contains several fittings, each of which has relatively small size, also raising missing detection. Figure 1 shows some real-world examples for each challenge. It is clear that the small size of fittings, as well as various kinds of disturbance items, brings heavy obstacle for the rust detection. The current methods all work to improve feature representation in front end, e.g., using attention mechanism and pyramid architecture. However, for the rust detection, these front-end improvements cannot guarantee the valid detection for small-size fittings and effectively eliminate the disturbance from the background environment. According to our literature survey, there have been some studies for solving similar problems. For instance, Zhai et al. [15] proposed a new cascade reasoning graph network for multi-fitting detection on transmission lines. This network incorporates three kinds of domain knowledge, i.e., cooccurrence knowledge, semantic knowledge and spatial knowledge, to represent the corelation of different mini-size fittings. With these knowledge reasoned by graph attention network, more discriminative features can be extracted based on the original visual features to recognize and position the fittings. However, this method still works in front end and is devoted to feature enhancement before generating accurate proposals. It aims to develop the detection accuracy and does not consider the disturbance of complex background which will reduce the detection robustness. As shown in Figure 2, missing detection, as well as false detection, has occurred many times in our experiment when running the methods discussed above. For an actual applications, missing detection and false detection should be significantly avoided from the rust detection, especially in online scenarios. For online tour-inspection, UAVs, which are equipped the detection algorithms, need to provide more reliable and robust detection results. It is necessary to enhance the feature representation of fittings based on the current front-end techniques to improve the detection accuracy and robustness.

(**a**) (**b**) (**c**) (**d**)

**Figure 1.** Rust detection examples of transmission line fitting with different challenges: (**a**,**b**) are of small object, while (**c**,**d**) are with complex background.

**Figure 2.** Examples of missing detection and false detection using Lim's method [12] that is an improved version of Faster R-CNN for small-size object detection, where (**a**,**b**) are the two examples regarding of missing detection and false detection.

Based on the analysis mentioned above, the main challenge for the rust detection of fittings in complex environment is developing feature representations of the fittings against the background disturbance. We observe an interesting phenomenon from our empirical evaluations. Despite of many cases of missing detection and false detection, Faster R-CNN can still obtain the interesting regions, most of which have a certain degree of feature representation of the fittings. In other words, most of the region proposals in Faster R-CNN actually are related to the fitting object. Then it motivates us a new idea: enhance the feature representation from these regions themselves. Following this idea, we build a new Faster R-CNN model for the robust rust detection of transmission line fitting in this paper. The backbone network, VGG16 network, is replaced by a deeper network ResNet-101 for extracting more rich information about the fitting object from UAV images. More importantly, a new feature enhancement mechanism is built after the region of interest (RoI) pooling layer to improve the feature representation of the regions that have real fittings. The weight of the disturbance terms can then be relatively reduced. Comparative results on some real-world UAV monitoring images verify that the proposed model can significantly increase the detection accuracy and robustness.

The main contributions of this paper can be summarized as follows: (1) From the application perspective, this paper proposes an lightweight but effective solution for the rust detection of transmission line fittings. The proposed method is simple and of high accuracy as well as robustness. To our best knowledge, the study of the rust detection for transmission line fittings is still at its infancy. (2) From the theoretical perspective, this paper constructs a new feature enhancement mechanism in the back end of classical object detection algorithms. Different from most of current methods, this mechanism can further enhance the feature representation based on the generated features. This mechanism can apply for the current two-stage detection methods without too much modification on the algorithmic architecture. We believe this mechanism can provide a different aspect to improve the detection reliability and robustness.

The remaining part of this paper is as follows. Section 2 is dedicated to the implementation of the classical Faster R-CNN. Section 3 describes the proposed model in detail. Section 4 carries out a set of comparative experiments, followed by a conclusion in Section 5.

#### **2. Background of Faster R-CNN**

Faster R-CNN was developed from R-CNN and Fast R-CNN. R-CNN is the first algorithm to apply CNN to an object detection task. It uses a selective search algorithm to obtain region proposals with fine-tuning the CNN, and trains a support vector machine (SVM) classifier that also performs border regression. This method does not work end-toend. Based on the spatial pyramid pooling network (SPP-Net [16]), Fast R-CNN inputs the whole image instead of each candidate region into R-CNN for feature extraction, also with the region proposals generated through selective search. The biggest improvement of Faster R-CNN is the use of region proposal network (RPN) to generate regions of interest (ROI), no longer using the selective search strategy again. Another interesting point is that the whole training process can run under GPU environment, indicating computationally inexpensive.

The classical Faster R-CNN algorithm is composed of an RPN network and a Fast R-CNN network. The whole architecture includes four parts: convolution layer, RPN layer, ROI pooling layer and classification regression layer, as shown in Figure 3. To improve readability, here we take the rust detection problem as an example to describe the algorithmic details.

**Figure 3.** Architecture of Faster R-CNN. To improve readability, we take the rust detection problem as the background.

Faster R-CNN first scales each UAV image of size *P* × *Q* to the size of *M* × *N*, then inputs the image to a CNN network (e.g., the commonly used VGG16) to obtain a feature map. The feature map is then fed into the RPN that generates region proposals on the feature map. The object category, as well as its position, in the region proposals can be

obtained through the classification regression layer. Specifically, the RPN distinguishes between the foreground and background of region proposals, and outputs the region proposal in the foreground region. The ROI pooling layer reshape the region proposal in foreground region to a fixed size (7 × 7) by combining the CNN features and RPN information. The region is connected to a detection network for judging the object category and fine-tuning its position as well.

In Figure 3, the foreground classification is the key. It requires to compare the region proposals with the ground-truth box manually annotated by experts, and further calculates the intersection ratio of the two boxes, defined as Intersection over Union (IoU): *IoU* = *<sup>A</sup>*∩*<sup>B</sup> <sup>A</sup>*∪*<sup>B</sup>* , just as shown in Figure 4. When the IoU of one region proposal is greater than 0.7, the region is set as positive sample, i.e., the foreground. If the IoU < 0.3, the region is set as negative sample, i.e., the background. The region proposal with IoU value of 0.3–0.7 is not involved in the training. The positive and negative samples are then used to train RPN.

**Figure 4.** Schematic drawing of IoU calculation, where A and B are the ground-truth box and region proposal respectively.

#### **3. The Proposed Faster R-CNN Model with Feature Enhancement**

In this section, a new Faster R-CNN model is proposed, in which two developments are made: updating the backbone network using the residual network ResNet-101, and designing a new feature enhancement mechanism after the ROI pooling layer. The structure of the proposed model is shown in Figure 5. The details will be elaborated as follows.

**Figure 5.** Structure of the proposed new Faster R-CNN model with feature enhancement.

#### *3.1. Feature Extraction Network with ResNet-101*

It has been proven that deeper network is capable of extracting more robust feature representations from images [17]. However too deep a network would raise gradient disappearance and gradient explosion. Residual connection is an effective trick to extend the depth of a deep convolutional network [18], as shown in Figure 6. Here we adopt a version of a well-known residual network, named ResNet-101 network, as the backbone network of Faster R-CNN. This network is believed to obtain richer feature information, thereby improving the feature representation for the detection. The structure is listed in Table 1. Since ResNet-101 has been widely studied, we would not analyze it in detail. Please find the reference [18] for the implementation details.

**Figure 6.** Structure of residual connection.

**Table 1.** Structure of the ResNet-101 network used in this paper.


#### *3.2. Feature Enhancement Mechanism*

Motivated by the self-attention mechanism that can extract richer information by learning the similarity between the target object and the other ones, this section builds a new back-end feature enhancement mechanism after the ROI pooling layer rather than in the feature extraction network. Through calculating the similarity between each region proposal and the others, the feature weights of the region proposals containing target object can be enhanced via the overlaying of the object's representation. This operation is based on the observation: the majority of the obtained region proposals contain fitting objects, and their feature representations are essentially similar. Then we can put weights on the features of the obtained region proposals. The weight on the region with fittings is pushed to be greater, while the weight on the region with disturbance item is reduced. Then the feature representation of the actual fittings can be enhanced to reach a more robust detection. This operation is called *feature enhancement mechanism*.

Specifically, denote the input as the feature map *Ri* that is the output of ROI pooling layer. The calculation process is as follows:

(1) Without loss of generality, calculate the similarity between the 1st region proposal and the other ones, and get the weight *α*1,*i*:

$$\alpha\_{1,i} = \frac{\det(R\_1 \, R\_i)}{\sqrt{d\_k}} \tag{1}$$

where *R*<sup>1</sup> and *Ri* are the feature map of the 1st and *i*th region proposal respectively, *dot*(·, ·) means dot product that is chosen as the similarity measure, *dk* is the input feature dimension. Certainly, different similarity measures can also be adopted.

(2) Normalize the weight *α*1,*<sup>i</sup>* via the Softmax layer to obtain the final weights *α* 1,*i* , as shown in Equation (2). This operation can also be visualized in Figure 7. Obviously, *α* 1,*i* indicates the influence of the *i*-th region proposal on the 1st region proposal.

**Figure 7.** Sketch of calculating the similarity between region proposals.

(3) Multiply *α* 1,*<sup>i</sup>* by *Ri* and sum up all region proposals to obtain an output *α* <sup>1</sup> that has the same dimension as the input data, as shown in Equation (3). The operation can be visualized in Figure 8.

**Figure 8.** Sketch of calculating the output feature.

From the analysis mentioned above, the feature enhancement mechanism can increase the weight of the target region proposals by calculating the similarity between the obtained regions, which helps to eliminate missing detection. Obviously, for the regions proposal containing disturbance item, the weight will be relatively decreased, which helps to lessen false detection.

#### *3.3. Loss Function*

The loss of the whole network mainly consist of classification loss *Lcls* and regression loss *Lreg*, as follows:

$$L(\{p\_{\dot{l}}\}, \{t\_{\dot{l}}\}) = \frac{1}{N\_{\rm cls}} \sum\_{\dot{l}} L\_{\rm cls}(p\_{\dot{l}\prime}p\_{\dot{l}}^{\*}) + \lambda \frac{1}{N\_{\rm reg}} \sum\_{\dot{l}} p\_{\dot{l}}^{\*} L\_{\rm reg}(t\_{\dot{l}\prime}t\_{\dot{l}}^{\*}) \tag{4}$$

where *i* is the anchor index, *pi* is the probability of the *i*-th anchor to be predicted as the ground-truth label, *p*∗ *<sup>i</sup>* = 1 if it is a positive sample, otherwise *p*<sup>∗</sup> *<sup>i</sup>* = 0. *ti* is a vector representing the four parameterized coordinates of the predicted bounding box, and *t* ∗ *i* is that of the ground-truth box associated with a positive anchor. *λ* is the regularization parameter to tradeoff *Lcls* and *Lreg*. The calculation of *Lcls* and *Lreg* are as follows:

$$L\_{\rm cls}(p\_{i\prime}p\_{i\prime}^\*) = -\log[p\_i p\_i^\* + (1 - p\_i)(1 - p\_i^\*)]\tag{5}$$

$$L\_{\text{reg}}(t\_{i\prime}t\_{i\prime}^\*) = \sum\_{i \in \{\ge, y, w, h\}} smooth\_{L\_1}(t\_i - t\_i^\*) \tag{6}$$

where {*x*, *y*, *w*, *h*} denotes the two coordinates of the box center, width and height, *smoothL*<sup>1</sup> (*x*) = 0.5*x*<sup>2</sup> if|*x*|<<sup>1</sup> <sup>|</sup>*x*|−0.5 otherwise .

To improve the readability of the proposed method, we provide the flowchart of the methodology in Figure 9. The key of the methodology is the proposed feature enhancement mechanism. Please note that this mechanism can also be applied to the other two-stage object detection architectures.

**Figure 9.** Flowchart of the whole methodology.

#### **4. Experimental Results**

In this section, the effectiveness of the proposed model is verified. The programming environment is Linux Mint 19.2, PyTorch 1.0 and CUDA10.2, configured with GeForce RTX1080 graphics card.

#### *4.1. Dataset Preprocessing*

The UAV images used in this experiment come from our real-world application, as shown in Figure 10. The dataset consists of 245 images for training and 105 images for test, which were collected from southern China. We use the LabelImg software to mark the images containing transmission line fittings as 'Fittings'. To increase the amount of training data, the images are preprocessed through data enhancement, such as adjusting brightness, adding noise, mirror processing, translation and rotation processing. The effect after data enhancement is shown in Figure 11.

**Figure 10.** Examples of UAV images used in this experiment. For better illustrative effect, we divide the examples into the two groups, as shown in the subfigures (**a**,**b**).

**Figure 11.** UAV images of transmission line fittings after data enhancement. The columns (**a**–**d**) are the four examples to show the effect of data enhancement..

#### *4.2. Ablation Validation*

In the proposed model, we integrate ResNet-101 and feature enhancement mechanism in Faster R-CNN. To evaluate the effect of each component, we set up two ablation experiments for the evaluation.

#### 4.2.1. Change of Backbone Network

To evaluate the effect of feature extraction network, we replace the ResNet-101 by another CNN network, i.e., VGG16. The comparative results are shown in Figure 12. It is clear that ResNet-101 can better recognize the rusted fittings with lower missing detection rate. It can demonstrate that deeper network can enhance the feature learning of the target object by extracting abundant feature information, which enables the model to obtain high accuracy in the following detection. The stronger the learning ability of the model is, the more robust the detection and the higher the detection accuracy will be. However, we also find that the phenomenon of missing detection has not been completely eliminated, also shown in Figure 12.

**Figure 12.** Comparative results of different backbone networks, where (**a**,**c**) are with VGG16, (**b**,**d**) are with ResNet-101.

#### 4.2.2. Change of Feature Enhancement Mechanism

Figure 13 clearly shows the detection results adding the feature enhancement mechanism. The corresponding feature heat maps are shown in Figure 14. The results demonstrate that the false detection in the red box can be better improved by employing the feature enhancement mechanism. Obviously, the feature enhancement mechanism can fully utilize the correlation of region proposals to enhance the feature information of the target region. Meanwhile, it improves the ability of accurate recognition of the target objects, which can well solve the problem of false detection in complex background environments.

(**c**) (**d**)

**Figure 13.** Comparative results for evaluating the feature enhancement mechanism, where (**a**,**c**,**e**) are the results without the mechanism, (**b**,**d**,**f**) are the results of using the mechanism.

(**a**) (**b**)

(**e**) (**f**)

**Figure 14.** Feature heat maps of the results in Figure 13, where (**a**,**c**,**e**,**g**) are the results without the mechanism, (**b**,**d**,**f**,**h**) are the results of using the mechanism.

#### *4.3. Comparative Results*

We also employ the two indexes, Recall and Precision, to numerically evaluate the detection performance, as listed in Table 2. To provide a straightforward comparison, we further plot the P-R curve of the methods, as shown in Figure 15. The formulations of the two indexes are as follows:

$$Recall = \frac{TP}{TP + FN} \tag{7}$$

$$Precision = \frac{TP}{TP + FP} \tag{8}$$

where *TP* (True Positive) represents positive samples that are correctly classified; *FP* (False Positive) is negative samples wrongly categorized as positive ones; *TN* (True Negative) is the negative examples that are correctly classified; *FN* (False Negative) represents negative samples wrongly categorized as positive ones.

From Table 2 and Figure 15, the detection accuracy of the proposed model is much higher than the other methods. This further demonstrates that the proposed model can effectively improve the feature representation and enhance the feature information of the region of interest as well. Consequently, a more accurate and reliable detection of rusted fittings can be achieved.

**Table 2.** Numerical comparison of the Faster R-CNN model based on the backbone network.


**Figure 15.** P-R curve of the models for comparison. The closer the curve is to the right-hand top corner, the better the detection performance will be.

Furthermore, we compare the proposed model with two typical object detection algorithms SSD and YOLOv3. We also introduce a state-of-the-art small-size object objection algorithm, called Lim's method [12], for comparison. The detection results of the four methods are shown in Figure 16. No surprisingly, the proposed model gets the best detection performance, which proves again the effectiveness of the feature enhancement mechanism.

**Figure 16.** P-R curve of the four object detection methods.

#### **5. Conclusions**

In this paper, a new robust Faster R-CNN model is proposed for the rust detection of transmission line fitting. This model aims at solving the two challenges of the rust detection: disturbance of complex environment and small size of fitting object. The proposed model focuses on the feature enhancement based on the obtained region proposals. With the proposed feature enhancement mechanism, the feature representation of the rusted fittings can be improved in an targeted mode. Moreover, the mechanism is of good application universality, since it can work on different kinds of two-stage detection architectures. With self-learning the rich information about the fitting object, the detection robustness as well as accuracy can then be developed with much lower missing detection rate and false detection rate. Then the reliability of the detection results can be much improved. The proposed model is easy to implement and has better deployment capacity for real-world applications, especially for online scenarios.

In future works, we plan to exploit the structured information about fittings. It can be observed that the appearance of fittings must be accompanied with transmission lines, which indicates sort of structured information. This information is believed beneficial for the rust detection. Moreover, for an actual engineering, the trustworthy decision is more preferable. Interpretability analysis will be applied to the rust detection. How to understand the detection results is another interesting problem. Online rust detection should be also paid more attention since online tour-inspection is an actual demand for UAV applications. In our current engineering, the online detection task is made by loading the offline-trained detection algorithm into the UAV, which motivated our study in this paper, i.e., enhancing the robustness of Faster R-CNN. We think another feasible solution is updating the detection algorithm online with the sequentially-collected images, i.e., in an incremental mode. For example, if the UAV tours around some special terrains, such as forest, villages, rivers, etc., the images with such terrain characteristics should bring more kinds of feature representation for the detection algorithm. The detection model is then required to be updated automatically and incrementally. How to incrementally update the online detection model is interesting, of course, not easy to realize, for the online tour-inspection. We will study this problem in the future work.

**Author Contributions:** Conceptualization, Z.G. and Y.T.; Methodology, Y.T.; Software, Y.T.; Validation, Y.T. and W.M.; Formal analysis, W.M.; Investigation, W.M.; Resources, Z.G.; Data curation, W.M.; Writing–original draft preparation, Y.T.; Writing–review and editing, Z.G.; Visualization, Y.T.; Supervision, Z.G.; Funding acquisition, Z.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Science and Technology Project of SGCC (Research and application of audiovisual active perception and collaborative cognitive technology for smart grid operation and maintenance scenarios) (No. 5600-202046347A-0-0-00).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


### *Article* **A Data-Driven OBE Magnetic Interference Compensation Method**

**Yizhen Wang, Qi Han \*, Dechen Zhan and Qiong Li**

Faculty of Computing, Harbin Institute of Technology, Harbin 150001, China

**\*** Correspondence: qi.han@hit.edu.cn; Tel.: +86-1393-662-2926

**Abstract:** Aeromagnetic compensation is a technology used to reduce aircraft magnetic interference, which plays an important role in aeromagnetic surveys. In addition to maneuvering interferences, the onboard electronic (OBE) interference has been proven to be a significant part of aircraft interference, which must be reduced before further interpretation of aeromagnetic data. In the past, most researchers have focused on establishing linear models to compensate for OBE magnetic interference. However, such methods can only work using accurate reference sensors. In this paper, we propose a data-driven OBE interference compensation method, which can reduce OBE interference without relying on any other reference sensor. This network-based method can integrally detect and repair the OBE magnetic interference. The proposed method builds a prediction model by combining wavelet decomposition with a long short-term memory (LSTM) network to detect and predict OBE interference, and then estimates the local variation of the magnetic field to remove the drift of the interference. In our tests, we construct 10 semi-real datasets to quantitatively evaluate the performance of the proposed method. The F1 score of the proposed method for OBE interference detection is over 0.79, and the RMSE of the compensated signal is less than 0.009 nT. Moreover, we also test our method on real signals, and the results show that our method can detect all interference and significantly reduce the standard deviation of the magnetic field.

**Keywords:** aeromagnetic survey; OBE interference compensation; LSTM network

**Citation:** Wang, Y.; Han, Q.; Zhan, D.; Li, Q. A Data-Driven OBE Magnetic Interference Compensation Method. *Sensors* **2022**, *22*, 7732. https:// doi.org/10.3390/s22207732

Academic Editors: Yongbo Li, Bing Li, Jinchen Ji and Hamed Kalhori

Received: 14 August 2022 Accepted: 7 October 2022 Published: 12 October 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

Aeromagnetic surveys are a way of measuring the Earth's magnetic field using a magnetometer mounted on an aircraft. Because of their low cost and high efficiency, aeromagnetic surveys are widely used in archaeological surveys, geological research, mineral exploration, and so on [1]. Aeromagnetic compensation is a technique for eliminating the interference of aircraft and is a key part of aeromagnetic surveys. The most widely used compensation model, called the T-L model, was proposed by Tolles and Lawson. They divided the aircraft interference field into three parts: constant magnetic field, induced magnetic field, and eddy current field, and established a linear regression model to estimate these interferences [2]. However, in practical applications, some onboard electronic (OBE) systems will generate magnetic interferences (called OBE interferences in this paper), which is not described in the T-L model. With the improvement of the accuracy of magnetic field measurement and aeromagnetic compensation, the influence of OBE interference becomes more and more obvious. Therefore, it is very important to continuously monitor and compensate for the OBE interference in aeromagnetic exploration.

At present, there are not many studies on OBE interference compensation. These methods always rely on reference sensors, such as current/voltage sensors and reference magnetometers. In [3–5], the current are voltage are measured by corresponding sensors and used to estimate the OBE magnetic interference according to the assumption that the OBE magnetic interference is proportional to the current or voltage. However, there are some problems with this kind of method: (1) The ON/OFF of OBE devices may not be an instantaneous operation [6], during the switching process, the interference is very hard to be well compensated without a very good synchronization between the current/voltage sensor and the magnetometer. (2) On the ground, the ambient magnetic field is too complex and noisy, and the calibration parameters calculated on the ground are hard to be made accurate. The method proposed in [7] is different from the above methods. The authors use a reference magnetometer to adaptively estimate the OBE magnetic interference. This method is particularly effective for time-varying interference with unknown signal characteristics. However, there is also a problem: this method can only handle the case that there is only one interference source, but there are more interference sources on aircraft. Moreover, there is a common problem with these two kinds of methods relying on reference sensors: the installation of reference sensors is difficult, and they must meet Civil Aviation requirements.

In [8], the authors propose a method without reference sensors (called the COOE method in this paper). In this work, the OBE magnetic interference is detected according to the difference between the variances in the magnetic field in the adjacent windows. Then the interference is roughly compensated using linear interpolation. The flaw of this method is that some parameters, such as the sliding window length and the detection threshold, should be manually adjusted carefully to avoid lots of missing alarms. Moreover, the linear interpolation compensation misses magnetic details.

In our opinion, the method of time series processing can be used to detect and compensate for the OBE interference. In recent years, methods based on neural networks, especially long short-term memory (LSTM), have been widely used in context anomaly detection. Since aeromagnetic data are stored in chronological order [9], it can be regarded as a time series. Meanwhile, the OBE interferences can be regarded as context-dependent anomalies. Hundman et al. [10] used LSTM networks to detect anomalies in spacecraft data and proposed a dynamic threshold segmentation method based on past data. Markus Thill et al. [11] proposed an unsupervised ECG anomaly detection method based on stacked LSTM. The anomalies in ECG sequences can be detected by predicting normal sequence behavior and establishing a statistical model of normal behavior prediction error. In [12], the LSTM and Bi-LSTM models are generated to identify rice crops using a Sentinel-1 time series. Ding et al. [13] used an LSTM model to detect errors in industrial manipulator systems. In the remote sensing field, Sun et al. applied LSTM to crop yield prediction [14], and Wang et al. proposed a land cover classification and supervision framework based on LSTM [15].

To avoid the limitations of OBE interference compensation methods with reference sensors, we proposed a data-driven method without any reference sensors. This method is more accurate and reliable than other data-driven methods. We use the LSTM network to identify OBE interference and design a pipeline to repair it. Compared with the COOE method, our proposed compensation method can reduce OBE interference while preserving the details of magnetic field data. To quantitatively evaluate the proposed method, 10 semi-real datasets are constructed using real measured magnetic fields and simulated interferences, each of which contains about 10% OBE interferences. We also test the proposed method in real data and compare it with the COOE method.

In this work, we propose an integrated method to detect and repair the OBE magnetic interference without relying on any reference sensor. To detect the OBE magnetic interference, we propose an LSTM-based network to predict the normal magnetic field and calculate an adaptive threshold of the error between the prediction result and the measured magnetic field. Before that, we use the maximum overlap discrete wavelet transform (MODWT) to decompose the magnetic field into multi-resolution terms, which makes the prediction more accurate. After the detection, we analyze two typical OBE interference types and propose an algorithm to repair them using the mentioned prediction result and the local signal variation. In addition, we also utilize a Gaussian kernel convolution to remove the trend term, which can be embedded into a network and improve the model generalization.

The organization of this paper is as follows. The proposed method is described in Section 2. Section 3 discusses the experimental details, including dataset preparation, model parameters, training configuration, and evaluation metrics. The results and discussion are provided in Section 4. The conclusions are summarized in Section 5.

#### **2. Methodology**

#### *2.1. Background*

An aeromagnetic survey usually refers to the collection of geomagnetic signals using aircraft installed with a high-precision magnetometer. The aircraft is always equipped with various OBE devices, such as air conditioning, beacon light, radio, and so on [6,16]. When these devices are working, they produce DC currents, which can generate a magnetic field [6]. These magnetic fields are projected onto the geomagnetic field and are captured by the scalar magnetometer. In this case, the OBE interference can be described by Biot– Savart law and can be viewed as proportional to the current [3,4]. When the device is operating for a long time, OBE interference usually presents as a long-term magnetic field shift. When the device is switched frequently, such as the beacon light continuously sends a pulsed current [17], OBE interference also appears as a pulse-like signal. All these magnetic interferences are named OBE interferences. Without corresponding compensation, such interferences will disturb magnetic field analysis. Therefore, it is very important to compensate for the OBE interference in aeromagnetic surveys.

There are two typical forms of OBE interferences. One type is short-time, peaky, and usually periodic. This kind of interference is called short-term interference, such as interference caused by the beacon light or strobe light. The other type of interference changes quickly at the beginning and end, like short-term interference, but lasts longer. This kind of interference is called long-term interference, such as interference from the radio or rudder motor.

In Figure 1, two typical forms of magnetic interference are given. Figure 1a is the short-term interference with peaks. Figure 1b is the long-term interference with data drift.

**Figure 1.** Typical OBE generated by the software interferences: (**a**) the short-term interference, (**b**) the long-term interference, and (**c**) the description of interferences.

For the two types of OBE interferences, we give a brief description in Figure 1c, and use light yellow boxes to mark the start and end positions of OBE interferences. In our work, the interferences at the start and end positions are named anomalies, and the part between the start and end positions is named the drift. As can be seen from the figure, these two types of interferences are in the same form at the start and end positions; the difference is there is a data drift in the long-term type. Therefore, we can adopt a unified approach to deal with the two kinds of interference and then deal with the data drift separately.

Based on the above analysis, we propose an unsupervised automatic OBE interference compensation framework, as shown in Figure 2. First, the magnetic field data is detrended, normalized, and decomposed using a wavelet. These procedures are collectively called preprocessing. Then we use the LSTM network to model and predict the magnetic field data and detect interferences by comparing the errors between the prediction and the input with an adaptive threshold based on extreme value theory (EVT). With the detection results, we repair the interference anomalies with the predictions and then repair the data drift with calculated bias. Finally, the repaired results are superimposed with the magnetic field trend to obtain the final compensated results.

**Figure 2.** Overview of method structure.

#### *2.2. Data Preprocessing*

#### 2.2.1. Trend Separation

The Earth's magnetic field has a large magnitude of about 50,000 nT [9] and wide range of variation. The variation comes from the sum of diurnal variations, micropulsations, magnetic storms, and long-term variations [18]. Due to the above influence, the variation in the Earth's magnetic field is complex, and its variation range can reach several hundred nanoteslas in a period of time [9,19]. However, the magnitude of OBE interference is very small. For example, a beacon light can only cause magnetic interference of about 0.25 nT, which is far less than the variation of the Earth's magnetic field. Therefore, the large variations of the Earth's magnetic field should be reduced.

As a low-pass filter, a Gaussian filter can extract the low-frequency part of the signal [20]. We take the extracted low-frequency part as the extraction trend and reduce the influence of variations of the Earth's magnetic field and enhance the interference characteristics by subtracting this part. In addition, since the filter is actually implemented by convolution, this part can be directly embedded into the network.

The Gaussian kernel is expressed as [21]

*<sup>g</sup>*(*s*, *<sup>σ</sup>g*) = <sup>1</sup> <sup>√</sup>2*πσ<sup>g</sup> e* <sup>−</sup> *<sup>s</sup>*<sup>2</sup> 2*σ*2 *<sup>g</sup>* , (1)

where *σ<sup>g</sup>* is the standard deviation of the Gaussian function, and *s* is the magnetic field value of a sampling point.

For the measured magnetic field **S***Total*, the trend obtained after filtering is **S***Trend* = *conv*(**S***Total*, **g**), where *conv* represents the one-dimensional convolution operation. The smoothness of the trend can be adjusted by changing the *σ<sup>g</sup>* value. Empirically, the *σ<sup>g</sup>* is set as 50. The magnetic field after removing the trend is expressed as **S** = **S***Total* − **S***Trend*.

#### 2.2.2. Normalization

The data should be normalized into the network to obtain better performance. Because it is hard to give a suitable minimum or maximum value, we use z-score normalization as follows:

$$\mathbf{X} = \frac{\mathbf{S} - \mu\_s}{\lambda \sigma\_s},\tag{2}$$

where the standardized value is **X**, the input data are **S**, the average of the input data is *μs*, and the standard deviation of the input data is *σs*.

To keep the value range as close to [−1, 1] as possible, we bring *λ* in to scale *σs*. In our work, *λ* is experimentally set as 4.

#### 2.2.3. Wavelet Decomposition

Wavelet transform often plays a positive role in time series prediction models [22–24]. As stated in [25], wavelet transform is an adaptive time-frequency domain analysis method that can effectively deal with non-stationary time series and non-Gaussian noise and is an effective data enhancement method. For complex sequences, combining wavelet transform with a prediction model has been a widely used data enhancement method, which can effectively improve the accuracy of prediction results [26]. In addition, wavelet transform has a better processing effect for non-stationary signals such as magnetic field measured by a magnetometer [27]. Therefore, based on the good time-frequency localization ability of the wavelet [22], we use wavelet decomposition to decompose the magnetic field data into different scales and obtain the near-periodic expression. We feed signals of different scales into the LSTM network to obtain more accurate prediction results.

In this paper, the MODWT is adopted to process a magnetic field with the trend removed. The MODWT is a time-shift invariant, redundant, and non-orthogonal transformation method. Compared with other wavelet transform methods, the MODWT has the following advantages [25]: (1) the ability to handle any sample size; (2) increased coarse scale resolution; (3) a more asymptotically efficient estimation of wavelet square difference than DWT; (4) it can deal with non-stationary time series and non-Gaussian noise more effectively. More information about MODWT can be found in [22,26,28].

Here, MODWT [29,30] is used to analyze the time series. When the MODWT is applied to the time series **X**, the wavelet and the scaling coefficients of level *j* are

$$\tilde{\mathcal{W}}\_{j,t} = \sum\_{l=0}^{L\_j - 1} \tilde{h}\_{j,l} \mathbf{X}\_{t - l \text{mod}\mathcal{N}'} \tag{3}$$

$$
\tilde{V}\_{j,t} = \sum\_{l=0}^{L\_j - 1} \tilde{\mathfrak{F}}\_{j,l} X\_{t - l \text{mod}N'} \tag{4}
$$

where {˜ *hj*,*l*} *Lj*−1 *<sup>l</sup>*=<sup>0</sup> and {*g*˜*j*,*l*} *Lj*−1 *<sup>l</sup>*=<sup>0</sup> are the wavelet and scaling filters of level *j*, respectively. The filter width is expressed as

$$L\_{\hat{j}} = (2^{\hat{j}} - 1)(L\_1 - 1) + 1,\tag{5}$$

where *L*<sup>1</sup> is the width of the unit-level Daubechies wavelet coefficient. In order to avoid the influence caused by the boundary conditions of the wavelet transform, we first remove the *Lj* wavelet and the scale coefficients (determined by Equation (5)), and obtain the "boundary correction" wavelet and scaling coefficients [26].

#### *2.3. OBE Interference Detection*

In this section, we design an unsupervised OBE interference detection algorithm. The algorithm is trained on normal magnetic field data without OBE interference and outputs the prediction error. The Extreme Value Theory (EVT) method is used to set the threshold automatically to detect OBE interference.

#### 2.3.1. Magnetic Field Predictor Based on LSTM

The LSTM network was first proposed by Hochreiter and Schmidhuber as a variant of RNN [15]. The LSTM can make full use of the historical information and time dependence of modeling signals [31]. It allows subsequent states at different time intervals to be stored by regularly connecting hidden layer nodes, where parameters are shared between different parts of the model. Many studies have proven that LSTM is very suitable for processing time series data [10,32–34]. The basic structure of the LSTM includes an input gate, forgetting gate, output gate, and internal storage unit. In this paper, we adopted the LSTM formula proposed by Graves et al. [35], and the key formulas are expressed as

$$\begin{aligned} f\_t &= \sigma(\mathcal{W}\_f[h\_{t-1}, \mathbf{x}\_t] + b\_f), \\ i\_t &= \sigma(\mathcal{W}\_l[h\_{t-1}, \mathbf{x}\_t] + b\_l), \\ \tilde{\mathsf{C}}\_t &= \tanh(\mathsf{W}\_{\mathbb{C}}[h\_{t-1}, \mathbf{x}\_t] + b\_{\mathsf{C}}), \\ \mathsf{C}\_t &= f\_t \* \mathsf{C}\_{t-1} + i\_t \* \tilde{\mathsf{C}}\_t, \\ o\_t &= \sigma(\mathcal{W}\_o[h\_{t-1}, \mathbf{x}\_t] + b\_o), \\ h\_l &= o\_l \* \tanh(\mathsf{C}\_t), \end{aligned} \tag{6}$$

where the subscript *t* − 1 is the previous time, and *t* is the current time. The operator ∗ represents the Hadamard product. Two activation functions are used, the hyperbolic tangent function *tanh*(·) and the sigmoid function *σ*(·). The different weight matrices *W* and deviation *b* are the parameters to be trained. The forgetting gate *f* determines how much information is forgotten in the output *ht*−<sup>1</sup> and the current time input *xt*. Similarly, the input gate *i* determines which values will be updated. *C* is the updated state for the current cell. The output gate *o* determines which parts of the cell state will be output. Finally, the updated hidden state is *h*.

Suppose that the sequence of the preprocessed magnetic field is **X**<sup>∗</sup> = [**x**1, **x**2, ... , **x***n*], where **<sup>x</sup>***<sup>i</sup>* <sup>∈</sup> **<sup>R</sup>***m*(*<sup>i</sup>* <sup>=</sup> 1, ... , *<sup>n</sup>*) is the *<sup>m</sup>*-dimensional component of magnetic field data of the *i*th sampling point obtained after wavelet decomposition, and *n* is the number of points. Assuming the embedding time step of LSTM is *s* and the predicting step length *l* = 1. For the magnetic field of the *t* sampling point, the sequence **X**∗ *<sup>t</sup>* = [**x***t*−*s*, **x***t*−*s*+1, ... , **x***t*−1] with length *s* is used as the input of the LSTM network to obtain the predicted value **x**ˆ*t*. At time *t* − 1, the structure of the multidimensional LSTM is shown in Figure 3 below.

The predicted sequence is **X**ˆ <sup>∗</sup> , **X**ˆ obtained using wavelet reconstruction to **X**ˆ <sup>∗</sup> .

**Figure 3.** Structure of the LSTM.

2.3.2. Adaptive Threshold Based on EVT

Predictors are trained on normal magnetic data without OBE interference. The output prediction error *Et* is defined as the square of the difference between the observed value and the corresponding predicted value at time *t* and is expressed as

$$E\_t = (X\_t - \hat{X}\_t)^2. \tag{7}$$

All prediction errors in the training dataset are expressed as a vector **E** = [*E*1, *E*2, ... , *En*]. We use the Peak Over Threshold (POT) method based on EVT to adaptively adjust the detection threshold. The EVT is a statistical theory aiming to find patterns of extreme values, which are usually at the tail of probability distributions [36]. The advantage of EVT

is that there is no need to make assumptions about data distributions when looking for extreme values [37,38]. The POT is the second theorem in EVT, and the basic idea of that is to use the generalized Pareto distribution (GPD) with parameters to fit the tail of the probability distribution. The GPD function is as follows:

$$\bar{F}\_{th}(e) = P(E - th > e \mid E > th) \sim (1 + \frac{\gamma e}{\sigma})^{-\frac{1}{\gamma}},\tag{8}$$

where *th* is the initial threshold, and *γ* and *σ* are the shape and scale parameters of GPD. *E* represents the data point in **E**, and *E* − *th* represents the portion that exceeds the threshold *th*. The values of parameters *γ* and *σ* are estimated by maximum likelihood estimation (MLE). Then the final threshold *Th* is calculated [37].

In the detection process, the prediction error **E** of the test dataset is used to detect whether there is OBE interference in the test dataset. If a value in **E** exceeds threshold *Th*, the corresponding magnetic data are considered anomalies. The detection result can be expressed as **E***<sup>a</sup>* = {*Ei* ∈ **E**|*Ei* > *Th*}, *i* = 1, 2, ... , *n*. The detected anomalies **D** are represented as continuous sequences of *Ea* ∈ **E***a*; the start and end point of *j*th detected anomaly *Dj* are denoted as *ano<sup>j</sup> begin* and *ano<sup>j</sup> end*, respectively.

#### *2.4. OBE Interference Repair*

The mode of OBE interference includes not only spike-like short-term interference but also drift-like long-term interference. The detection results of these two types of OBE interference are shown in Figure 4. The parts where the OBE interference appears and disappears are detected and marked with red backgrounds. These parts (named anomalies) should be corrected. In addition, notice that if the OBE interference is long-term, there is a data drift in the magnetic field, which should also be corrected.

**Figure 4.** Detection results of two kinds of OBE interferences: (**a**) A detected short-term interference. (**b**) A detected long-term interference. The anomaly segments are marked with a red background.

For each detected anomaly, the repair method can be decomposed into two steps: (1) repair anomaly segment with predicted points; (2) if a data drift exists, repair it with the evaluated drift magnitude.

#### 2.4.1. Anomaly Repair

The anomaly segments are detected according to the method in Section 2.3. For each detected anomaly segment, we use the LSTM-based predictor to obtain the predicted value and repair the segment point by point. For each point in this anomaly segment, we obtain the predicted value and update the corresponding point using the LSTM network. The operation is looped until this anomaly segment is processed.

Algorithm 1 shows the detailed steps of anomaly repair.

**Data:** detection result **D***a*, test dataset **X**


**<sup>2</sup> for** *k = 1:* (*ano<sup>j</sup> end* <sup>−</sup> *ano<sup>j</sup> begin*) **do**


```
6 end
```
**<sup>7</sup>** do wavelet reconstruction and get the repaired value **X**ˆ ;

#### 2.4.2. Drift Repair

As can be seen from Figure 4, when there is a drift between the data in front of the anomaly and the data behind the anomaly, the mean value of the data will change significantly. Therefore, the existence of drift can be judged according to the change in the mean value.

Drift repair can be expressed in the following steps:


Algorithm 2 shows the detailed steps of drift repair.

#### **Algorithm 2:** Drift repair algorithm

**Data:** an detected anomaly *Dj*, test dataset **X** after anomaly repair **Result:** drift repaired value**X <sup>1</sup>** get the data in front of the anomaly segment **<sup>X</sup>***pre* <sup>=</sup> {*Xanoj*−<sup>1</sup> *end* ,..., *Xanoj begin*−1 } and the data behind the anomaly segment **<sup>X</sup>***a f ter* <sup>=</sup> {*Xanoj end*+1 ,..., *Xanoj*<sup>+</sup><sup>1</sup> *begin* } ; **<sup>2</sup>** get the intercepts **<sup>X</sup>***pre*\_*cut* <sup>=</sup> {*Xanoj begin*−*npre* ,..., *Xanoj begin*−1 } and **<sup>X</sup>***a f ter*\_*cut* <sup>=</sup> {*Xanoj end*+1 ,..., *Xanoj begin*+*na f ter* }, where *npre* <sup>=</sup> *min*{(*ano<sup>j</sup> begin* <sup>−</sup> <sup>1</sup>) <sup>−</sup> *anoj*−<sup>1</sup> *end* , *npre*\_*num*}, *na f ter* <sup>=</sup> *min*{(*anoj*+<sup>1</sup> *begin* <sup>−</sup> *ano<sup>j</sup> end* + 1), *na f ter*\_*num*} ; **<sup>3</sup>** calculate the change in mean <sup>Δ</sup>*avg* <sup>=</sup><sup>|</sup> **<sup>X</sup>**¯ *pre*\_*cut* <sup>−</sup> **<sup>X</sup>**¯ *a f ter*\_*cut* <sup>|</sup> ; **<sup>4</sup>** for **X***pre*\_*cut*, calculate the threshold value *thavg* using Neyman – Pearson (N – P) criterion ; **<sup>5</sup> if** Δ*avg* > *thavg* **then <sup>6</sup>** calculate *di f f*(**X***a f ter*) = *Xt* <sup>−</sup> *Xt*−1, *<sup>t</sup>* <sup>∈</sup> (*ano<sup>j</sup> end* <sup>+</sup> 1, *anoj*+<sup>1</sup> *begin*] ; **<sup>7</sup>** repair drift using **X** *a f ter* <sup>=</sup> *Xanoj end* + ∑ *di f f*(**X***a f ter*) ; **<sup>8</sup> else <sup>9</sup>** process the next drift ;

```
10 end
```
#### **3. Experimental Details**

*3.1. Data Preparation*

3.1.1. Background Data

The background data were obtained outdoors using an optically pumped cesium vapor magnetometer and a self-made collector. The setup of the collection experiment is the same as our previous work [39]. We set the sampling rate to 10 Hz, as is often used in aeromagnetic surveys [40]. The range of the measured magnetic field is about 50,104 to 50,130 nT. We collected magnetic field data for two days, denoted as Day1\_Data and Day2\_Data, and used them for training and validation, respectively.

#### 3.1.2. Semi-Real Data

In aeromagnetic surveys, it is difficult to obtain a pure magnetic field that can be used as ground truth when OBE interference exists. Therefore, we use semi-real data to quantitatively evaluate the performance of the proposed method. We generate 10 semi-real datasets using the background field and simulated OBE interferences. We divide the background magnetic field Day2\_Data into 10 segments of similar length and denote them as C-1 to C-10. For each segment, we randomly select 10% to add simulation OBE interferences. The OBE interference is generated according to the model in [4]. The interference length is 1–40 sampling points, and the intensity is between 0.04 nT and 0.2 nT. The 10 semi-real datasets are called S-1 to S-10. Figure 5 shows a segment of the semi-real datasets.

**Figure 5.** A segment of the semi-real datasets.

#### 3.1.3. Magnetic Field with Real OBE Interferences

To validate the method, we also collected real magnetic field data with OBE interference, as shown in Figure 1. Figure 1a,b are interferences caused by a beacon light and a radio, respectively. The collection equipment consists of a scalar magnetometer and a three-axis vector magnetometer, and the sampling rate is consistent with Section 3.1.1. The interference magnetic field of the beacon light is collected by the scalar magnetometer mounted at the end of the elongated tail. The interference magnetic field of the radio is composed of the vector magnetometer data.

#### *3.2. Model Training*

#### 3.2.1. Model Parameters

The structure and parameters of the LSTM model are shown in Table 1.


**Table 1.** Parameters of LSTM model.

The model consists of two hidden layers, each with 128 cells. We find that this structure provides sufficient capacity to predict the magnetic field well. Increasing the model size provides few benefits. Moreover, the sequence length *s* is set as 32, it provides a good balance between performance and time cost. The model was trained for 200 iterations.

In order to train the LSTM model correctly, the background dataset Day1\_Data is divided into an 80% training set and a 20% verification set. The Adam optimizer uses the mean absolute error (MAE) as a loss function. The batch size and learning rate are set as 64 and 0.0001, respectively. For the training step, we use a computer with an Intel Core i7-8700K CPU and an NVIDIA GeForce 1080Ti GPU.

#### 3.2.2. Evaluation Indicators

Since the real datasets lack the ground truth value, we use different evaluation indicators on semi-real datasets and real datasets to evaluate the performance of the proposed method.

For the semi-real datasets, the real label of OBE interference and the pure magnetic field is known, so the detection label and the compensated magnetic field are quantitatively evaluated with the ground truth.

We use the modified range-based precision, recall, F1 score [41,42], and the ROC (Receiver Operating Characteristic) curve to evaluate the interference detection performance. For a dataset, the set of real anomaly ranges is denoted as *R*, *Ri* represents the *i*th real anomaly range, and *Nr* is the number of real anomalies. The set of detected anomaly ranges is denoted as *P*, *Pi* is the *i*th detected anomaly range, and *Np* is the number of detected anomalies. The *Rewardexistence* represents the number of intersections between the detected anomaly ranges and the real anomaly ranges, and the *Rewardoverlap* is a function that describes the overlap situation between the detected anomaly ranges and the real anomaly ranges (see [41] for more details). The recall is defined as Equation (9), which colligates the *Rewardexistence* and *Rewardoverlap* using factor *α*. In this paper, we set *α* as 0.5 to treat *Rewardexistence* and *Rewardoverlap* fairly. The precision is defined as Equation (10), which mainly evaluates the overlap between real anomaly ranges and detects anomaly ranges. The F1 score is defined as Equation (11), which is the harmonic mean of the above two metrics and reflects the robustness of the detection algorithm. The ROC represents the detection ability of the model when the discriminative threshold changes, and the area under the ROC curve (often referred to as AUC) is used to measure the probability that a model can be classified correctly [43].

$$\text{Recall} = \frac{\sum\_{i=1}^{N\_r} \left( a \times \text{Reward}\_{\text{existence}}(\mathbf{R}\_{i\prime}P) + (1 - a) \times \text{Reward}\_{\text{overlap}}(\mathbf{R}\_{i\prime}P) \right)}{N\_r},\tag{9}$$

$$\text{Precision} = \frac{\sum\_{i=1}^{N\_p} \text{Reward}\_{overlap}(R\_\prime P\_i)}{N\_p} \,\tag{10}$$

$$\text{F1 score} = \frac{2 \times \text{Precision} \times \text{Recall}}{\text{Precision} + \text{Recall}}.\tag{11}$$

For the OBE interference repair results, we use the root mean squared errors (RMSE) to evaluate the performance [44,45], which is defined as Equation (12). The RMSE represents the deviation between the repaired results and the clean reference data of all OBE interferences. Therefore, the calculation of RMSE is dependent on the detection results and only includes the repair results of the interference part rather than the complete magnetic field series.

$$\text{RMSE} = \sqrt{\frac{1}{n} \sum\_{i=1}^{n} (\mathbf{X}\_i^{\prime} - \mathbf{X}\_i^{pure})^2} \text{ \tag{12}$$

where **X** *<sup>i</sup>* is the *<sup>i</sup>*th repaired OBE interference, **<sup>X</sup>***pure <sup>i</sup>* is the corresponding pure magnetic field without interference, and *n* is the number of anomalies.

For real datasets, we adopted the standard deviation (STD) and improvement ratio (IR) after bandwidth filtering to evaluate the repair results, which are common performance evaluation indexes in magnetic exploration [46,47].

#### **4. Results and Discussion**

#### *4.1. Selection of the Mother Wavelet and Wavelet Decomposition Level*

The selection of the mother wavelet and the level of wavelet decomposition will affect the performance of the predictor and then affect the detection and repair effect. In this paper, several orthogonal mother wavelets of the wavelet family are compared, including haar (db1), Daubechies (db2, db3, db4, db5), Symlets (sym3, sym5), and Coiflets (coif1, coif2) [48]. For each mother wavelet, we explore 1 to 5 levels of decomposition, denoted as L1 to L5, respectively. At the same time, we also test the case without using wavelet decomposition, denoted as L0. The selection of mother wavelet and decomposition levels is based on the prediction results in the pure magnetic field. We use RMSE to measure the prediction error, and the lower RMSE value represents that the predictor performs better. We use Day\_1 Data for training and test each combination of mother wavelet and decomposition level on C-1 to C-10 and calculate the average RMSE for 10 datasets, as shown in Table 2.


**Table 2.** The RMSEs using different mother wavelets and decomposition levels. **Bold** indicates the best performance.

The test results show that using the haar wavelet with four decomposition levels can achieve the best prediction performance. The haar wavelet has three advantages: (1) it can capture the fluctuations between adjacent data, (2) it does not have an aliasing problem, and (3) it can express the low-frequency features well [49]. The main components of the geomagnetic field are also mostly distributed in low frequency [50], so using the haar wavelet decomposition is conducive to the LSTM network for learning geomagnetic field characteristics.

#### *4.2. Semi-Real Dataset Test Results*

The detection results on the semi-real dataset are shown in Table 3. We compare the proposed detection method with the COOE method. Note that the parameters of the COOE method are set to the maximum of the F1 score of each dataset.


**Table 3.** Comparisons of detection results between the proposed method and the COOE method in semi-real datasets.

As shown in Table 3, our method performs better in three indicators. The recalls of our method are larger than 0.92 in all datasets and equal to 1 in 7 datasets, which indicates that our method can detect almost all anomalies. However, the recalls of the fine-tuned COOE algorithm are below 0.7 in all datasets, mostly in the range of 0.5 to 0.64, which means that it misses a lot of anomalies. The precisions of our method are significantly higher than the COOE method, which indicates that our method mislabels fewer normal cases. The precisions of the COOE method are very small because the anomaly length defined in this paper is small, and the detection results of the COOE method often have large biases from real anomalies. Moreover, the higher F1 scores of our method indicate that our method has better comprehensive capacity.

We also compare the ROC curves and calculate the AUCs of the two methods, as shown in Figure 6. The ROCs and AUCs were calculated by counting all detection results from S-1 to S-10. Our method achieves a higher AUC value, which is close to 1, indicating that our method can detect anomalies more accurately.

Figure 7 shows the detection results of a segment of semi-real data. We can find that our method can detect the OBE interference well, almost consistent with the ground truth. The COOE method can also detect all the OBE interferences, but there are some offset points for some OBE interferences. This phenomenon is caused by the fixed-length and non-overlapping sliding window used in COOE. As a note, the fixed window length is hard to pick.

Based on the detection results above, we compare the OBE interference repair performance of the proposed method with that of the interpolation method in COOE on the same datasets. The comparisons of the repair results are shown in Table 4.

**Table 4.** Comparisons of repair results between the proposed method and the COOE method in semi-real datasets.


As shown in Table 4, the proposed method obtained lower RMSEs in all datasets, indicating that the error between the repaired results of our method and the pure magnetic field without OBE interferences is smaller.

Figure 8 shows the repair results of the same semi-real data segment. It can be seen that the results of our method are closer to the ground truth. However, due to the deviation of the detection results, the interpolation results obtained by the COOE method are not satisfactory, and many OBE interferences can not be completely removed.

**Figure 8.** The comparison of repair results in the semi-real data segment.

#### *4.3. Cross-Validation of Semi-Real Datasets*

To verify the generalization of the model, we conduct cross-validation on Day\_2 Data. We use the leave-one-out method to carry out 10-fold cross-validation and calculate the detection and repair criteria of 10 validations. For each fold, we select 9 segments from C-1 to C-10 datasets for model training and use the remaining dataset with OBE interferences for verification. For example, for the first fold, we choose C-1 to C-9 as the training set and use S-10 as the test dataset. The final indicators of cross-validation are shown in Table 5. These results are similar to those in Section 4.2 using Day\_1 Data as the training set, which proves the good generalization of our model.


**Table 5.** The average criterions of 10-fold cross-validation.

#### *4.4. Real Datasets Test Results*

We verified the effectiveness of the proposed method on the real collected magnetic fields with beacon light and radio interferences; these two datasets are named R-1 and R-2. Since a clean magnetic field cannot be obtained for comparison in real aeromagnetic measurement, we compared the detection results with manual marks and evaluated the magnetic series before and after repair using standard deviation (Std) and improvement ratio (IR).

Figure 9a,b are the detection and repair results of beacon light interferences, respectively. From this figure, we find that our method can detect all the interferences and repair them well. In contrast, the COOE method misses some interferences in the detection procedure, which leads to the repaired signal having obvious mutations.

**Figure 9.** Detection and repair results of the aircraft's beacon light interferences.

Figure 10a shows the detection results of the aircraft's radio interference, both methods can detect the interference. Figure 10b shows the repair results of the radio interference. Because the COOE method uses the linear interpolation method to repair, its result completely loses the data characteristics of the magnetic field. Our method can retain the characteristics of the magnetic field data.

Table 6 shows the STDs and IRs of repaired magnetic fields after band-pass filtering. It is a conventional operation in aeromagnetic compensation to perform band-pass filtering before calculating and evaluating STD and IR, and the filter bandwidth is set to 0.1 to 0.6 Hz, according to [51]. Our method achieves smaller STD and larger IR in R-1. The IR of COOE is even less than 1 because the data mutation caused by linear interpolation will produce a large oscillation after filtering. Although the COOE method has a better STD and IR in R-2, the data characteristics of the magnetic field are totally removed. In conclusion, our method can effectively eliminate the current magnetic interference without affecting the characteristics of the magnetic field.

**Table 6.** Comparisons of repair results between the proposed method and the COOE method on real datasets.


**Figure 10.** Detection and repair results of the aircraft's radio interference.

#### **5. Conclusions**

In aeromagnetic surveys, the interference of OBE magnetic fields is non-negligible. The previous methods often use reference current or magnetic sensors to estimate and remove the OBE interference. In this paper, we propose an unsupervised and unreferenced method to integrally detect and repair them. The detection is determined by an LSTM-based predictor. The threshold of the error between the prediction result and the measured magnetic field is adaptively calculated by the POT algorithm. Moreover, wavelet decomposition is also utilized to improve the prediction accuracy. After detection, based on the prediction result, we design an algorithm to repair the OBE magnetic interference. In contrast with the methods based on interpolation, our method can retain the detailed signal characteristics. In addition, we embed a Gaussian kernel convolution layer into the network, which can detrend the signal and improve the model generalization.

We compare the proposed method with a previous work relying on no reference sensor. On semi-real datasets, it is shown that our proposed method is better than the COOE method in the range-based recall, precision, F1 score, AUC, and RMSE. On real datasets, the results also show that our method can effectively compensate for OBE interferences and increases the improvement ratio. In addition, our method can retain the normal magnetic field characteristics in long-term interference, which ensures the validity of the magnetic field data.

**Author Contributions:** Conceptualization , Y.W.; methodology, Y.W.; software, Y.W.; validation, Y.W.; formal analysis, Y.W. and Q.H.; investigation, Y.W.; resources, Q.H., D.Z. and Q.L.; data curation, Y.W. and Q.H.; writing—original draft preparation, Y.W.; writing—review and editing, Y.W.; visualization, Y.W.; supervision, Q.H.; project administration, Y.W.; funding acquisition, Q.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Natural Science Foundation of China (Grant No. 61771168).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** This work was supported by the National Natural Science Foundation of China (Grant No. 61771168). The authors would like to thank the Institute of Information Countermeasures Technology for providing deep learning servers.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Damage Monitoring of Engineered Cementitious Composite Beams Reinforced with Hybrid Bars Using Piezoceramic-Based Smart Aggregates**

**Hui Qian, Yuqing Zhang, Yuechang Li, Jundong Gao \* and Jianxue Song**

School of Civil Engineering, Zhengzhou University, Zhengzhou 450001, China

**\*** Correspondence: gaojundong@zzu.edu.cn

**Abstract:** In order to explore the crack development mechanism and damage self-repairing capacity of ECC beams reinforced with hybrid bars, the smart aggregate-based active sensing approach were herein adopted to conduct damage monitoring of ECC beams under cyclic loading. A total of six beams, including five engineered cementitious composite (ECC) beams reinforced with different bars and one reinforcement concrete counterpart, were fabricated and tested under cyclic loading. The ultimate failure modes and hysteresis curves were obtained and discussed herein, demonstrating the multiple crack behavior and excellent ductility of ECC material. The damage of the tested beams was monitored by smart aggregate-based (SA) active sensing method, in which two SAs pasted on both beam ends were used as actuator and sensor, respectively. The time domain analysis, wavelet packet-based energy analysis and wavelet packet-based damage index analysis were performed to quantitatively evaluate the crack development. To evaluate the self-repairing capacity of the beams, a self-repairing index defined by the difference of damage index at loading and unloading peak points was proposed. The results in time domain and wavelet packed analysis were in close agreement with the observed crack development, revealing the feasibility of smart aggregate-based active sensing approach in damage detection for ECC beams. Especially, the proposed damage self-repairing index can describe the same structural re-centering phenomena with the test results, showing the proposed index can be used to evaluate the damage self-repairing capacity.

**Keywords:** PZT; smart aggregate; ECC beam; damage monitoring

#### **1. Introduction**

Engineered cementitious composites (ECCs) are fiber-reinforced cementitious composite materials [1,2], where an appropriate amount of polyvinyl alcohol (PVA) fibers are randomly distributed to form a three-dimensional space supporting system. Hence, the tensile strain capacity, the toughness, the durability and the impact fatigue of concrete members are significantly enhanced [3–5]. Specifically, two of the most important mechanical characteristics of ECC are the quasi-strain hardening properties and the multiple micro-cracking behavior with self-controlled crack widths [3,6]. The ultimate tensile strain attained by ECC is 200–600 times greater than that of regular concrete, and the multiple cracking behavior of ECC is distinguished from that of regular concrete. However, the damage monitoring of ECC is not well studied in the literature despite its great importance.

There have been several techniques for damage monitoring or health monitoring of structures in the past few decades [7]. Piezoceramic-based smart aggregates are multi-functional and can perform various tasks [8]. The PZT (Lead Zirconate Titanate) is the most popular piezoceramic material due to its strong piezoelectric effect [9], high bandwidth [10,11], fast response [12,13], and availability in different forms [14,15]. The applicability of PZT smart aggregates to health monitoring and damage detection of various civil structures has been demonstrated by experimental results. Gu et al. [16] conducted the early-age strength monitoring of concrete cylinder specimens and predicted the concrete strength development based

**Citation:** Qian, H.; Zhang, Y.; Li, Y.; Gao, J.; Song, J. Damage Monitoring of Engineered Cementitious Composite Beams Reinforced with Hybrid Bars Using Piezoceramic-Based Smart Aggregates. *Sensors* **2022**, *22*, 7184. https://doi.org/10.3390/s22197184

Academic Editors: Hamed Kalhori, Yongbo Li, Bing Li and Jinchen Ji

Received: 29 August 2022 Accepted: 20 September 2022 Published: 22 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

on the output voltage of the sensors embedded into the concrete structures. Jiang et al. [17] presented a stress wave-based active sensing method to detect the crack in FRP-reinforced concrete beams, and the results show that the developed piezoceramic-based active sensing method can monitor the crack-induced damage and estimate the process of damage degree in real-time. Song et al. [18] developed a smart aggregate-based impact detection and evaluation system, which had been used to detect the impacts on a concrete beam. Song et al. [19] also performed the structural health monitoring of a specially designed concrete bent-cap, indicating that the smart aggregate is able to capture the moment of concrete cracking. More recently, smart aggregates have been increasing employed in the monitoring field, involving the strength monitoring of early age concrete [20], impact detection and evaluation [21], health monitoring [22–27] and damage detection [28,29].

On the other hand, PZT smart aggregates were experimentally validated to be applicable to other fields, such as bolt looseness monitoring [30–33], soil freeze–thaw process monitoring [34], monitoring of water content in sandy soil [35,36], soil compaction monitoring [37], degree of water permeability [38,39], and damage diagnosis of hydraulic structure [40]. At the same time, to meet the monitoring environment under different conditions, the forms of piezoelectric intelligent aggregates are becoming more and more diverse. Gao et al. [41] designed, fabricated, and tested a novel embeddable tubular smart aggregate (TSA) based on a piezoceramic tube for use in two dimensional (2D) structures, and through test results showed that the proposed TSA is suitable for monitoring the health condition of a 2D concrete structure. Lu et al. [42] developed a novel piezoceramic stack-based smart aggregate (PiSSA) with piezoceramic wafers in series or parallel connection to increase the efficiency and output performance over the conventional smart aggregate. Moreover, the research study on wireless smart aggregates (WSAs) was conducted by Yan et al. [43], and the efficiency of the WSA health monitoring system was experimentally validated in a bridge health monitoring system [44]. Voutetaki et al. [45] detected and evaluated the damage severity of shear critical concrete beams with a new portable real-time wireless impedance monitoring system.

In order to explore the crack development in ECC beams, the smart aggregate-based active sensing approach are herein adopted to conduct damage monitoring of ECC beams under cyclic loading. Five beams made of ECC materials reinforced with different types of hybrid reinforcement materials and one reinforcement concrete counterpart were designed and fabricated. The one-way cyclic loading test results of six beams are presented to discuss the reinforcement effect of composited materials on crack development, together with the self-centering effect. In addition, the crack development and self-centering effect are investigated by the time domain analysis and wavelet packet analysis.

#### **2. Principle of Damage Monitoring**

#### *2.1. Smart Aggregate-Based Active Sensing Approach*

In general, piezoceramic materials cannot be directly used in structural health monitoring, owing to the inherent fragility. A smart aggregate is thus designed by sandwiching the PZT patch into two marble blocks with epoxy, as illustrated in Figure 1a. Meanwhile, the cable with a Bayonet Neill–Concelman (BNC) connector is soldered to the PZT patch of the smart aggregate, as shown in Figure 1b. The smart aggregates are employed in the active sensing approach. Specifically, one smart aggregate connected to the waveform generator is used as an actuator to send excitation waves, and other distributed smart aggregates are regarded as sensors to simultaneously detect the propagated signals [34]. The values of wave amplitude and transmission energy will decrease with the occurrence of the cracks or damage inside the concrete, and the dropped values are associated with the degree of damage inside.

**Figure 1.** Piezoceramic-based smart aggregate.

In this paper, the damage development of ECC beams under cyclic loading is monitored by the active sensing approach, wherein the cracks or damage inside the tested beams can be reflected by the signals recorded by the sensors. During the damage monitoring of the ECC beam, the stress wave will change with the generation and closure of cracks. The signal received by the collector can be analyzed by the wavelet packet algorithm in terms of energy and damage index.

#### *2.2. Wavelet Packet Analysis*

The principle of wavelet packet analysis can be explained by the Figure 2, wherein the sensor signal S is decomposed by an *n*-level wavelet packet decomposition into 2*<sup>n</sup>* signal sets {*X*1, *X*2,... ,X2 *<sup>n</sup>*}. *Xj* is given by

$$X\_{\mathbf{j}} = [\mathbf{x}\_{\mathbf{j},1} + \mathbf{x}\_{\mathbf{j},2} + \dots + \mathbf{x}\_{\mathbf{j},m}] \tag{1}$$

where *m* is the number of sampling data and *j* is the frequency band (*j* = 1, ... , 2*n*). The energy (*Ei*, *<sup>j</sup>*) of the band signal *j* at time *i* is defined as,

$$E\_{i,j} = \|X\_j\|\_2^2 = \mathbf{x}\_{j,1}^2 + \mathbf{x}\_{j,2}^2 + \dots + \mathbf{x}\_{j,m}^2 (j = 1, 2, 3, \dots, 2^n) \tag{2}$$

where *i* is the time index.

**Figure 2.** Wavelet packet decomposition signal.

The difference between the signatures of health and damaged states have been commonly compared by root mean-square deviation (RMSD), which is the widely used damage index for health monitoring of concrete structures. The energy vector at health states *E*0,*<sup>j</sup>* = [*E*0,1, *E*0,2, ... , *E*0,2*n*], whilst the energy vector for the damaged data at time index *i* is marked as,

$$E\_i = \begin{bmatrix} E\_{i,1}, E\_{i,2}, \dots, E\_{i,2^n} \end{bmatrix} \tag{3}$$

Therefore, the damage index at time *i* can be defined as,

$$D = \sqrt{\frac{\sum\_{j=1}^{2^n} \left(E\_{i,j} - E\_{0,j}\right)^2}{\sum\_{j=1}^{2^n} E\_{0,j}^2}} \tag{4}$$

The transmission energy loss portion caused by structural damage can be quantitatively evaluated by the damage index *D*. When the concrete structure stays in a healthy state, the damage index is 0. While the damage index is larger than the initial value, indicating that damage appears in the concrete structures. Greater index indicates more serious damage. The value of damage index at the complete failure of a concrete structure is close to 1.

#### *2.3. Damage Self-Repairing Index*

For the re-centering structures, in order to evaluate the damage self-repairing capacity of the structure member, a damage self-repairing index is proposed in this paper. The self-repairing effect of structural members can be evaluated by proposed index, which is defined by the ratio of the difference of damage index at loading peak point, and unloading and the damage index at loading peak point. Therefore, the damage self-repairing index can be defined as,

$$R\_m = \frac{D\_{m,l} - D\_{m,n}}{D\_{m,l}} \tag{5}$$

where *Rm* represents the self-repairing ratio of damage during the cyclic loading; *m* stands for a certain cycle of structure loading; *l* stands for the loading peak point during the cycle; and *u* stands for the unloading phase during the cycle. If *Rm* = 1, it indicates the structure has a good self-repairing effect. While if *Rm* = 0, it indicates that the structure has no self-repairing effect.

By comparing the damage self-repairing index of each loading cycle, whether the structure has self-repairing effect can be obtained. By comparing the damage self-repairing index of different components with the same loading cycle, we can obtain the self-repairing condition between different components.

#### **3. Test of ECC Beams**

#### *3.1. Test Specimens*

A total of six concrete beams were designed with the same cross-sectional area of 100 mm × 100 mm, and the length was set to 1100 mm, as shown in Table 1. The specimen RC is a reinforced concrete beam, while another five specimens are made of ECC material, and are reinforced with steel rebars, steel strands, glass fiber reinforced plastics (GFRP) rods, shape memory alloy (SMA) rods, and both GFRP and SMA rods, respectively, as shown in Figure 3. The material of both reinforcements and stirrups is HRB 400, and the arrangements of stirrups are the same for six specimens. The crack changes of six beams can be schematically illustrated by Figure 4.

**Table 1.** Specimens' geometric figures.


(**a**) General view

(**b**) SMA anchorage (**c**) Steel strand anchorage (**d**) GFRP anchorage

**Figure 3.** Photo of six beam specimens.

**Figure 4.** Schematic drawing of crack changes.

#### *3.2. Materials*

The compressive strength of concrete was measured by the uniaxial compression test of three cubes with the same size of 150 mm × 150 mm × 150 mm, and the average value is equal to 44.02 MPa. The ECC material employed in the specimens was mixed with 2% of polyvinyl alcohol (PVA) fiber. The values of both tensile strength and compressive strength were determined, as summarized in Table 2.

**Table 2.** Material properties of ECC.


The SMA material selected in this study is a Ni-Ti alloy (56.35% Ni) with a diameter of 8 mm, and the diameter of GFRP and steel rebar is also 8 mm, while the nominal diameter of steel strand is 4.5 mm. The uniaxial tensile test at room temperature was performed to four kinds of reinforcement material, and the stress–strain relationships are presented in Figure 5. The brittle failure occurs for both steel strand and GFRP when the values of the tensile strength reach 1592.15 MPa and 881.35 MPa, respectively. The steel rebar and SMA experience four stages, involving the elastic stage, the elastoplastic stage, the plastic stage and the failure stage, the maximum tensile strengths of which reach 635.31 MPa and 802.35 MPa, respectively.

**Figure 5.** Stress–strain relationship of four different materials under uniaxial tensile loading.

#### *3.3. Test Setup and Loading Procedure*

The setup for the cyclic loading tests on six beam specimens is displayed in Figure 6a, wherein the span was designed as 1000 mm. Each beam specimen was tested under four-point bending condition. The load applied at mid-span was transmitted to the beam specimen by the spreader beam, and both beam ends were simply placed at two steel bearings including the pin support and the roller support. Two smart aggregates marked as SA1 and SA2, were attached to both ends of the beam with epoxy to detect the damage development and self-centering effect, wherein SA2 was responsible for recording the wave signal excited by SA1. Figure 6b is a picture of test setup.

(**b**) Photo of test setup

**Figure 6.** Test setup.

All the specimens were loaded according to the same loading protocol, as displayed in Figure 7. Displacement control was employed, and the displacement applied at mid-span gradually increased in equal increment of 2 mm. Only one cycle at each step was required. The loading process was terminated until the load decreased to 85% of the peak load.

**Figure 7.** Loading protocol.

The employed test system for damage monitoring is schematically shown in Figure 8. During the loading procedure, the repeated swept sine wave at a frequency range of 100 Hz to 130 kHz was generated by the function generator and then amplified by power

amplifier, and the PZT smart aggregate SA1 was stimulated. Afterwards, the stress wave propagated through the beam and was received by the sensor SA2. The detected wave signals were recorded and analyzed by the collector and the laptop with supporting software, respectively. It should be mentioned that the initial healthy state of each specimen monitored before applying load was supposed as the benchmark, and the subsequent damage development monitoring for the beam under cyclic loading was carried out. At the end of loading and unloading for each cycle, the signal propagated through the specimen was recorded by the damage monitoring system.

**Figure 8.** Schematic of test system for damage monitoring.

#### **4. Experiment Results and Analysis**

#### *4.1. Experiment Results*

The obtained failure modes and the crack widths of six tested beams are illustrated in Figure 9. Noticeable cracks from the tension zone at beam bottom until the neutral axis can be observed for each specimen. The termination of the loading process is associated with the occurrence of a major crack, except for specimen SS-ECC. This can be explained by the fact that the rupture of pretensioned steel strands occurred for specimen SS-ECC. It was also noted that only the concrete in the compression zone of the specimen RC was crushed, another five beam specimens exhibit multiple cracking behavior.

(**a**) RC 20 mm (**b**) R-ECC 26 mm

(**c**) SS-ECC 26 mm (**d**) GFRP-ECC 26 mm

(**e**) SMA-ECC 26 mm (**f**) GFRP/SMA-ECC 26 mm

**Figure 9.** Failure modes and the crack widths of tested beams.

Hysteresis curves of six specimens are illustrated in Figure 10, wherein the vertical coordinate is the applied load, and the abscissa is the displacement at the mid-span of the spreader beam. It is clear that larger displacement amplitude can be achieved by the beam specimens with ECC material, compared with the concrete beam specimen RC. Specimen GFRP-ECC experiences the largest deformation, and the applied load is also maximum among the tested beams, indicating the excellent load-carrying capacity. The load applied to specimen SS-ECC drops sharply owing to the fracture of steel strands. It should be mentioned that the load of beam specimen SMA-ECC is the least.

The residual deformation of the beam with increasing loading step is plotted in Figure 11. Significant residual deformation of six beams can be observed, implying the occurrence of a poor self-centering effect. At the same loading level, the residual deformation of the reinforced concrete beam specimen RC is higher than the others, while which of specimen GFRP-ECC is the least. It is noted that the values of specimens R-ECC and SMA-ECC in terms of residual deformation are the close to each other, and the cumulative residual deformation to these two beams are the most significant.

#### *4.2. Damage Monitoring Results*

Time-domain analysis is employed herein to reflect the development of cracks. The received sensor voltage signals of six beam specimens are shown in Figure 12, where the signal amplitudes collected at specified displacement level are presented. It can be clearly found that the signal strength in the case of no crack is greater than that at damage status. Furthermore, the amplitudes of sensor voltage decrease with increasing displacement, and the decline amount is gradually narrowing, reflecting that cracks develop significantly until the failure of the beam.

Figure 13 displays the energy indices with increasing loading cycle. It is clear that the energy propagated through the beam length gradually decreases with the development of cracks. The detected energy reduced to 25.6% of that at the health status, when the displacement applied to beam RC is equal to 2 mm, indicating the occurrence of crack, as illustrated in Figure 13a, wherein the final remaining energy is less than 1% of that in the initial state. Compared with detected energy at the initial status, the values of observed energy loss portion for five beam specimens with ECC material after the first loading cycle are equal to 56%, 6.6%, 63.3%, 18.9% and 23.4%, respectively, and the final energy ratios range from 1% to 10%. Therefore, it can be found that the energy loss of specimen R-ECC is less serious than its concrete counterpart RC, revealing that the development of cracks can be effectively suppressed by ECC materials. With the completion of the first loading step, the energy loss portion of beam SS-ECC among five ECC beams is the least, indicating the good reinforcement effect of steel strands. On the other hand, a poor reinforcement effect of GFRP is observed, since the beam GFRP-ECC displays the most apparent energy loss.

**Figure 10.** Hysteresis curves of tested beams.

The values of damage indices calculated by wavelet packet analysis are shown in Figure 14, including the values corresponding to both loading and unloading peak points. The initial value of damage index is equal to zero since each beam specimen is in health status. The increment of the damage index correlates with the increment of the loading step. Specifically, the damage index increases greatly with the application of the first-level loading, and the increment of beam RC is largest among all specimens. After the fifth-level loading, the values of damage index for six tested beams are very close to one, indicating extreme structural damage.

**Figure 11.** Residual deformation with increasing loading step.

**Figure 12.** Time-domain analysis.

The development of self-repairing index is herein shown in Figure 15, including both histogram and fitting curves. It is evident that the beam GFRP-ECC experiences mild self-repairing owing to the small values of self-repairing index throughout the loading process. Meanwhile, another five specimens exhibit noticeable self-repairing effect in the initial five levels of cyclic loading, and then drops sharply to near zero, indicating that no self-repairing phenomenon can be observed in the later stage. The unexpected self-repairing phenomenon may be explained by the fact that the unexpected adhesion effect between ECC concrete and hybrid bars was noticed owing to the smooth surfaces of hybrid bars. In addition, tested beam specimens without enough reinforcements yielded in the initial five loading levels, and experienced plastic deformation with growing loading cycles. Hence, the residual displacement was remarkable in the later loading stage, and little self-repairing phenomenon was observed.

(**b**) Damage index at each unloading peak point

**Figure 14.** Damage index.

(**a**) Self-repairing index

**Figure 15.** *Cont*.

(**b**) Fitting curves of self-repairing index

**Figure 15.** Self-repairing index.

#### **5. Conclusions**

The smart aggregate-based active sensing approach is employed to monitor the crack development of six concrete beams under cyclic loading. The results of failure modes, hysteresis cures and residual deformation were analyzed in detail. The time domain analysis and wavelet packet analysis in terms of energy indices and damage index were conducted, and the self-centering effect of tested beams were evaluated. The following conclusions can be drawn:

(1) Noticeable multiple crack behavior and better ductility was observed in beams with ECC. The ECC beam strengthened with GFRP exhibited favorable performance in terms of the load-carrying capacity, ductility, residual deformation and self-centering effect.

(2) The damage monitoring results were consistent with the observed crack development, indicating the feasibility of damage detection for ECC beams using smart aggregatebased active sensing approach. The experimental results provided the basis for the application of PZT smart aggregates in the ECC structures.

(3) The proposed damage self-repairing index can describe the same structural recentering phenomenon with the test results, showing the proposed index can be used to evaluate the damage self-repairing capacity.

**Author Contributions:** Conceptualization, H.Q.; methodology, H.Q.; software, Y.L. and Y.Z.; validation, Y.L.; formal analysis, Y.L.; investigation, Y.Z.; resources, Y.Z.; data curation, J.G.; writing original draft preparation, J.G. and Y.Z.; writing—review and editing, J.S.; visualization, J.G.; supervision, H.Q. and J.S.; project administration, H.Q.; funding acquisition, H.Q. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by National Natural Science Foundation of China, grant number 51978631.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** The research presented in this paper was supported by the National Natural Science Foundation of China (No. 51978631). The authors hereby would like to express their thanks for the support.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **A Comprehensive Operation Status Evaluation Method for Mining XLPE Cables**

**Yanwen Wang 1, Peng Chen 1,\*, Yanying Sun <sup>1</sup> and Chen Feng <sup>2</sup>**


**Abstract:** At present, the online insulation monitoring and fault diagnosis of mining cables are extensively discussed, while their operation status assessment has not been deeply studied. Considering that mining cables are closely related to the safe and stable operation of coal mine power supply systems, a comprehensive evaluation method including the Analytic Hierarchy Process (AHP), the membership cloud theory, and the D-S evidence theory is proposed in this paper in order to accurately assess the operation status of the mining XLPE cable. Firstly, the membership cloud is introduced to solve the index membership degree and the weights are calculated by an improved weight vector calculation method. Secondly, the conversion from the base layer indicator membership degree to the target layer trust degree is realized based on the D-S evidence theory. Then, the cable operation status is judged via the trust degree maximum and the distribution of conflict coefficients is further analyzed to warn the indicators with a bad status in the base layer. Finally, the feasibility of the proposed evaluation method is verified by a sufficient and detailed case analysis.

**Keywords:** status evaluation; mining XLPE cables; membership cloud; D-S evidence theory

**Citation:** Wang, Y.; Chen, P.; Sun, Y.; Feng, C. A Comprehensive Operation Status Evaluation Method for Mining XLPE Cables. *Sensors* **2022**, *22*, 7174. https://doi.org/10.3390/s22197174

Academic Editors: Yongbo Li, Bing Li, Jinchen Ji and Hamed Kalhori

Received: 30 August 2022 Accepted: 19 September 2022 Published: 21 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

As an important component of a coal mine's power grid, the mining XLPE cable (hereafter referred to as mining cable) is the core component of a coal mine's power supply system, and its operational safety status directly affects the stable operation of the coal mine's power grid and even concerns the production safety of the coal mine itself. The operating safety status of these mining cables is mainly affected by their operating environment and operating conditions [1]. Unlike the operating environment of cables in ordinary power grids, the air humidity in underground coal mines is high and the temperature varies greatly in different areas. Because of this, the insulation in mining cables is easily aged which leads to insulation degradation, and the space in underground coal mines is narrow, making mining cables susceptible to smashing, touching, and dragging, which cause the cable insulation to be damaged, wherein grounding or leakage faults can occur. The grounding method of the neutral point through to the arc extinguishing coil is usually used in a coal mine's power grid, although this grounding method can, in principle, allow the power supply system to operate with faults for 2–3 h. However, the underground environment of coal mines is different from that on the ground. The closed underground environment is filled with a large amount of gas and coal dust. When the mine cable discharges due to insulation failure or single-phase grounding, the generated sparks can easily lead to the environment catching "fire" and/or "explosive" conditions which could cause serious coal mine safety accidents such as electromechanical accidents, cable "release", underground fires, and explosions [2–5]. According to the national and provincial coal mine accident analysis reports released in 2021, there were 122 coal mine accidents and 225 fatalities in 2020 of which electromechanical, gas, and cable discharge

accidents accounted for 19% of the total number of accidents and 27% of the total number of fatalities, with some provinces accounting for a higher percentage of electromechanical, gas, and cable discharge accidents than the national accident rate. In some provinces, the proportion of electromechanical, gas, and cable discharge accidents was higher than that of average national accidents, and, among these, cable failure in coal mines was a direct cause of cable discharge accidents, an important factor that led to a number of electromechanical accidents and was a major external source of ignition-causing gas accidents.

Mining cable faults pose a long-term threat to the economic, safe, and reliable operation of coal mine power grids. Therefore, it is crucial to diagnose the operation status and fault monitoring of mine cables within coal mines. If we can intelligently sense and evaluate the operational safety status of cables in a coal mine power grid in real-time and accurately detect abnormalities and warnings of faults before they occur in order to prevent accidents before they happen, this will break through the current technical bottleneck and, with a small amount of investment, solve the actual demand problem and effectively reduce the personal and property losses caused by cable faults. This is a necessary precondition to ensure the safety of coal mine production and can bring about strong economic and social benefits and broad application prospects with significant research necessity and urgency.

In order to achieve an accurate assessment of the operating status of mining cables, this paper establishes a mining cable evaluation index system using hierarchical analysis based on a combination of expert industry opinions, relevant literature, protocols, and research. We first divide the mining cable status into severe, abnormal, attention, and normal, and then transform the standard status level into a visualized status space according to the cloud model theory; then, we calculate the membership degrees of quantitative and qualitative indicators, respectively; then, we use the improved AHP weight calculation algorithm to calculate the weights of each indicator in the indicator layer; and, finally, in order to reduce the uncertainty within the evaluation process, the fusion of the indicator membership and weights is realized step-by-step based on the D-S evidence theory and the current cable status is judged based on the fusion results.

The structure of this paper is as follows: Section 2 provides a literature review and summarizes previous research results; Section 3 introduces the AHP algorithm, the improved AHP weight calculation algorithm, and the D-S evidence theory; Section 4 describes the mining cable condition evaluation system and the calculation of the membership degrees of quantitative and qualitative indicators; Section 5 presents the detailed numerical work in the evaluation model; and, finally, the conclusions of this paper are given in Section 6.

#### **2. Literature Review**

Currently, cross-linked polyethylene (XLPE) power cables are widely used because of their excellent insulation and heat resistance properties [6–8]. Due to uncertainties in the design and production process, the frequency of faults has gradually increased, reducing the safety of the power grid. Cable insulation or fault monitoring methods mainly include temperature, DC components, dielectric loss, partial discharge, and traveling wave detection methods [9–11]. DC component and dielectric loss methods can only be used to perform overall insulation condition assessment with low accuracy. The local discharge method is influenced by the surrounding environment and signals propagation distance [12]. High-frequency signals decay rapidly making long-distance measurements difficult. Temperature monitoring methods are not sensitive enough for fault identification and the data generated by them is of limited use. For the traveling wave method, it is difficult to detect cable faults because the amplitude of the wave head is significantly attenuated after the reflected wave propagates through the long cable, and it is easily affected by the interference signal. To make up for the shortcomings of the previous methods, many novel methods for online cable monitoring have been proposed by many experts and scholars in recent years. Guangya Zhu [13] proposed a new online monitoring method of power cable insulation conditions based on low-frequency signal injection. For this method, a low-frequency signal is injected into the power system via the potential transformer's (PT)

open delta configuration. The cable conductor voltage and leakage current are detected. The interpolating windowed fast Fourier transform (FFT) algorithm is applied to calculate the dielectric loss angle. Then, the cable tangent delta (tanδ) can be deduced and the cable condition can be assessed. Wei Zhao [14] drew a two-dimensional trajectory map by simultaneously measuring two circulating currents in a coaxial cable. Fault criterion and databases are established to detect faults by analyzing the changes in trajectory characteristic parameters. Wenxia Pan [15] proposes a distributed online monitoring method for cable PD based on the phase-sensitive fiber-optic time domain reflection (ϕ-OTDR) principle. When the cable has PD, the backscattered Rayleigh light intensity change of the PD position is higher than the intact position. Yang Wu [16] proposed a monitoring scheme based on CM leakage current measurements at selected monitoring frequencies and developed an aging feature extraction method based on principal component analysis (PCA) which provides an estimate of the insulation's aging severity.

However, the online monitoring of the cable is only for online monitoring of a certain index of the cable and does not evaluate the overall running status of the cable in many aspects. In response to such problems, more and more scholars put their research focus on cable condition assessment. Heqian Liu [17] introduced the theory of dielectric response such as the isothermal relaxation current. Combined with the cable-aging equivalent model, the isothermal relaxation current peak-split fitting method, to represent the different processes of relaxation according to the attenuation characteristics of isothermal relaxation current, is provided. Lulu Li [18] proposed a non-invasive aging assessment method using transient disturbances originating from the system. The relative dielectric constant of the cable is extracted from the response of transient disturbances instead of the conventional dielectric loss angle in order to characterize the aging state more sensitively. Yanqun Liao [19] proposed a novel holistic approach in order to facilitate the implementation of risk-based maintenance strategies for cable conduits, cable terminations, joints, bodies, and grounding systems for each cable loop. Based on the polymer trap theory and the extended Debye model, the shape of the PDC curve, depolarization charge, parameters of the extended Debye model, aging factor (*A*), elongation at break retention rate (EB%), and their relationships under different thermal aging degrees were analyzed by Yiyi Zhang [20]. It is worth noting that there are many factors that affect the running state of the cable and no relevant scholars have carried out an in-depth analysis of the factors affecting the running state of the cable in all aspects and angles. Therefore, an index system that can scientifically evaluate the running state of the cable has not yet been established. At the same time, the comprehensive evaluation method of cable operation status also needs to be studied more deeply.

#### **3. Evaluation Models and Principles**

#### *3.1. AHP Method*

AHP is a multi-criteria decision-making method (MCDM) developed by T.L. Satty in order to evaluate and select alternatives based on a set of selected criteria [21]. This process can combine judgments from intangible qualitative criteria with tangible quantitative criteria. The specific steps of the AHP are as follows:

**STEP 1:** Establish a hierarchical structure.

We need to stratify the problem to be analyzed and establish a three-level structure model including the target layer, factor layer, and base layer. **STEP 2:** Construct a comparison judgment matrix.

$$A = \begin{bmatrix} a\_{11} & a\_{12} & \cdots & a\_{1n} \\ a\_{21} & a\_{22} & \cdots & a\_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a\_{n1} & a\_{n2} & \cdots & a\_{nn} \end{bmatrix} \tag{1}$$

We need to construct a matrix *A* according to the relative importance of each indicator. The basis for judging relative importance is shown in Table 1. *aij* represents the importance of the *i*-th indicator relative to the *j*-th indicator. When *i* = *j*, there is *aij* = 1; otherwise, there are *aij* > 0 and *aji* = 1/*aij*.

**Table 1.** Meaning of the scale.


**STEP 3:** Calculate the weight vector.

$$\omega\_{i} = \frac{\left(\prod\_{j=1}^{n} a\_{ij}\right)^{\frac{1}{n}}}{\sum\_{k=1}^{n} \left(\prod\_{j=1}^{n} a\_{kj}\right)^{\frac{1}{n}}} \tag{2}$$

According to Formula (2), *W* = [*ω*1, *ω*2,... , *ωn*] <sup>T</sup> can be calculated. **STEP 4:** Consistency test.

In order to verify whether the weight vector is reasonable, we need to check its consistency. Formulas for calculating the random consistency ratio *CR* are as follows:

$$CR = \frac{CI}{RI} \tag{3}$$

$$CI = \frac{\lambda\_{\text{max}} - n}{n - 1} \tag{4}$$

$$
\lambda\_{\text{max}} = \frac{1}{n} \sum\_{i=1}^{n} \frac{(\mathbf{A}\mathbf{W})\_i}{\omega\_i} \tag{5}
$$

where (*AW*)*<sup>i</sup>* is the *i*-th element of the product of *A* and *W*.

The average random consistency index RI of the multilevel matrix can be obtained from Table 2.

**Table 2.** Random consistency index.


If *CR* < 0.1, the comparison judgment matrix has a satisfactory consistency index; otherwise, the comparison judgment matrix needs to be readjusted.

In order to solve the problem that AHP is too subjective, we propose an improved weight vector calculation method. The method replaces the comparison judgment matrix with the interval judgment matrix and searches for the matrix with the highest consistency ratio in the interval to calculate the weight vector. This can reduce the subjectivity of experts and improve the objectivity of weights. The specific weight vector calculation has been introduced in detail in our previous study, and, therefore, in this paper, we will not describe too much due to the limitation of space and the details can be referred to in the literature [22].

#### *3.2. Membership Cloud Theory*

The membership cloud is a model proposed by academician Deyi Li in 1995 for converting qualitative evaluation and quantitative values, which represents the conversion relationship between numerical and linguistic values and can take into account the randomness and fuzziness of linguistic evaluation. The fuzziness is described by the width of the cloud and the randomness is described by the thickness of the cloud [23].

The definition of the membership cloud is: let *U* be an exact numerical representation of the universe of discourse. In the corresponding qualitative concept *A* on *U*, for any element x in the universe, there is a random number *y* ∈ [0, 1] with a tendency to be stable, which is called the membership degree of *x* to *A*, and the distribution of the membership degree in the universe is called the membership cloud, referred to as a cloud. Clouds are composed of several cloud droplets. Cloud droplets are quantitative descriptions of qualitative concepts. The generation process of cloud droplets expresses the uncertainty mapping relationship between qualitative concepts and quantitative values. According to the dimension of the universe *U*, a cloud can be divided into a one-dimensional cloud, a two-dimensional cloud, a multi-dimensional cloud, and so on.

The portrayal of the cloud model relies on three parameters, which are expectation *Ex*, entropy *En*, and super entropy *He*. Among these, *Ex* reflects the central value of a concept corresponding to a theoretical domain, *En* reflects the ambiguity of the concept, and *He* reflects the discrete degree of the cloud drops. The forward cloud generator forms a number of random numbers with stable tendencies based on the numerical characteristics of the cloud model to form an evaluation cloud; the inverse cloud generator calculates the numerical characteristics of the cloud model based on finite expert evaluation.

#### *3.3. D-S Evidence Theory*

In the 1960s, Dempster proposed the concept of evidence theory and his student Shafer redefined it and created the "mathematical theory of evidence", which was later called D-S evidence theory. Because D-S evidence theory has the ability of uncertainty reasoning and can represent, fuse, and decide uncertain information, it has been widely used in the fields of decision analysis, pattern recognition, and information fusion. The basic principles of D-S evidence theory are as follows [24,25]:

Assuming that *U* is the identification frame, it is a finite and complete universe of discourse, and *A* is a subset of *U*. If there is a set function m:P(*U*)→[0, 1] that satisfies the condition of Formula (6):

$$\begin{cases} \sum \, m(A) = 1\\ A \subseteq \mathcal{U} \\ \, m(\mathcal{Q}) = 0 \end{cases} \tag{6}$$

where *m* is the probability distribution function on the identification frame *U*. When *m*(*A*) > 0, *A* is a focal element; *m*(*A*) is the function value of the probability distribution function corresponding to event *A*. When the identification frame *U* is incomplete, *m*(∅) = 0. This paper only discusses the case when the identification frame *U* is complete and the elements are limited.

Dempster's rule for fusing N pieces of evidence is shown in Formula (7):

$$m(A) = \frac{1}{1 - K} \sum\_{A\_1 \cap \ldots \cap A\_M = A} \prod\_{i=1}^N m\_i(A\_i) \tag{7}$$

In Formula (7):

$$K = \sum\_{A\_1 \cap \dots \cap A\_M = \mathcal{D}} \prod\_{i=1}^N m\_i(A\_i) \tag{8}$$

where *K* represents the conflict size between the evidence bodies. When *K* = 1, the combination rule is invalid and the evidence bodies completely conflict. When *K*→1, the fusion

decision result may be contrary to common sense, so effectively resolving evidence conflicts is an important part of obtaining reliable fusion results.

#### **4. Empirical Application of the Evaluation Model**

In accordance with the requirements of AHP, we analyzed the main factors affecting the operation status of cables for coal mines, determined the interval judgment matrix of each layer, and obtained the optimal weight vector through the improved AHP algorithm. The membership matrix of each layer to the target layer was established through the membership cloud theory and the optimal weight vector was synthesized with the membership matrix by evidence synthesis through the D-S evidence theory in order to obtain a comprehensive judgment on the operating state of the mine cable.

#### *4.1. Section of the Voltage Situation Evaluation*

The operating condition of mining cables is influenced by a variety of factors. The selection of evaluation indexes plays a crucial role in the accuracy of evaluation results. In this paper, when constructing the evaluation index system, we strictly follow the five basic principles of systematicity, objectivity, measurability, scientificity, and hierarchy. Combining many references and expert opinions, a total of 11 individual indicators are selected which finally form a progressive evaluation index system as shown in Figure 1.

In Figure 1, x11 is the index of the insulation resistance test, x12 is the index of the pressure-tight test, and x13 is the index of the pulse current. x21 is the index of the leakage current, x22 is the index of the dielectric loss angle, x23 is the index of the core temperature of the cables, x24 is the index of the partial discharge, x31 is the index of the operating life, x32 is the index of the operating environment, x33 is the index of the load condition, and x34 is the index of the history of the fault.

**Figure 1.** Operation status evaluation model for mining XLPE cables.

#### *4.2. State Space*

In this paper, the operating status of the mining cables is classified as severe, abnormal, attention, and normal, denoted as *sk*(*k* = 1, 2, 3, 4). Subsequently, the boundary values *c1*, *c2*, and *c3* of adjacent states are determined using the Weibull distribution model based on real-time and historical data to determine the boundary interval (*d*min, *d*max)k of the kth state level as shown in Table 1. Where *d*min and *d*max are the left and right values of the kth bounding interval, respectively.

In view of the inherent vagueness of the division of state levels and the randomness of the appearance of each state, for this reason, this paper uses a cloud model to portray each state level, i.e., the state cloud. The grade boundary interval and the maximum possible interval of the kth state cloud (*Exk* − 3*Enk*, *Exk* + 3*Enk*) form an equation relationship, as shown in Equation (9), from which the expectation (*Exk*) and the entropy (*Enk*) of each state cloud can be obtained. At the same time, the distribution range of each state cloud is constrained by the limit condition of cloud image fogging, as shown in Formula (10), to obtain the super entropy (*Hek*) of each state cloud:

$$\begin{cases} Ex\_k - 3En\_k = d\_{\text{min}} \\ Ex\_k + 3En\_k = d\_{\text{max}} \end{cases} \Rightarrow \begin{cases} Ex\_k = \frac{d\_{\text{min}} + d\_{\text{max}}}{2} \\ En\_k = \frac{d\_{\text{max}} - d\_{\text{min}}}{6} \end{cases} \tag{9}$$

$$He\_k = \frac{En\_k}{18} \tag{10}$$

where *Exk* represents the point value that best reflects the *k*-th state level; *Enk* represents the measured range of the *k*-th state level; *Hek* represents the degree of cohesion of the data in the *k*-th state.

Send *Exk*, *Enk*, and *Hek* into the forward cloud generator to randomly generate the N cloud droplets (*g*, *μk(g*)), which can be visualized in the form of clouds for each operating state. The steps for generating cloud drops in the forward cloud generator are as follows:

Step 1: Enn = Randn(*Enk, Hek*). That is, *Enk* is the expectation and *Hek* is the standard deviation to generate a normally distributed random number *Enn*.

Step 2: g = Randn(*Exk*, *Enn*). That is, *Exk* is the expectation and *Enn* is the standard deviation to generate a normally distributed random number g.

Step 3: *<sup>μ</sup>k*(*g*) = exp <sup>−</sup>(*g*−*Exk* ) 2(*Enn*) 2 . The membership is calculated by this equation, and the number pair (*g*, *μk(g*)) represents a cloud droplet distributed over the theoretical domain.

Step 4: Repeat steps 1 to 3 until enough cloud droplets are generated (generally, N is 5000) to restore different operating states in the form of cloud models.

In addition, the subsequent solution of the quantitative and qualitative index membership in this paper is not the same, so two types of state spaces are formed as shown in Figures 2 and 3. In order to take into account the deterioration process of the quantitative index and the tolerance of the equipment to potential adverse factors, the adjacent state clouds in the quantitative space have a certain degree of transition trend; however, the qualitative space is only used as the limit measure, so the adjacent state clouds in the qualitative space are independent of each other.

**Figure 2.** State space of quantitative indicator.

**Figure 3.** State space of qualitative indicator.

#### *4.3. Quantitative Indicators*

In this paper, before determining the membership degree of quantitative indicators, the degradation degree function is introduced in order to normalize them and transform them uniformly to the range of [0, 1]. The quantitative indicators within the evaluation system can be divided into three categories: cost type (the smaller the measured value is, the better), benefit type (the larger the measured value is, the better), and interval type (the more centered the measured value is, the better), which are synthetically represented in Figure 4 and Equation (11) in this paper:

$$d = \begin{cases} d\_0 + \frac{1 - d\_0}{\mathbf{x}\_1 - \mathbf{x}\_{\min}} (\mathbf{x} - \mathbf{x}\_{\min}) & \mathbf{x}\_{\min} \le \mathbf{x} \le \mathbf{x}\_1 \\ 1 & \mathbf{x}\_1 \le \mathbf{x} \le \mathbf{x}\_2 \\ d\_0 + \frac{d\_0 - 1}{\mathbf{x}\_{\max} - \mathbf{x}\_2} (\mathbf{x} - \mathbf{x}\_2) & \mathbf{x}\_2 \le \mathbf{x} \le \mathbf{x}\_{\max} \\ 0 & \mathbf{x} \notin [\mathbf{x}\_{\min}, \mathbf{x}\_{\max}] \end{cases} \tag{11}$$

where *x* is the measured value of the quantitative index; *x*min and *x*max are the left and right values of the warning range; *x*<sup>1</sup> and *x*<sup>2</sup> are the left and right values of the allowed range; *d* is the degradation degree; and *d*<sup>0</sup> is the lower limit of the warning range.

**Figure 4.** Degradation calculation of different types of indicators.

After the degradation degree is calculated, the membership degree can be obtained by substituting the result into the expected curve Formula (12) of each state cloud in the quantitative space:

$$r\_k^{(1)} = \exp\left[-\frac{\left(d - E\boldsymbol{x}\_k\right)}{2\left(E\boldsymbol{n}\_k\right)^2}\right] \tag{12}$$

where *rk* (1) is the membership of the quantitative index in the *k*-th state.

#### *4.4. Qualitative Indicators*

Unlike quantitative indexes, the qualitative index membership is further determined only after the experts give their scores empirically by combining the field inspection situation with the test results. In this paper, to weaken the influence of subjectivity in a single expert, h (h = 10 in this paper) experts are invited to score; then numerical feature values are extracted from the discrete scoring results by an inverse cloud generator (Equations (13)–(15)) and then combined with a forward cloud generator to present them visually. The results are called floating clouds, as shown in Figure 5. The more discrete the cloud droplets of the floating cloud are, the greater the degree of disagreement among the experts.

**Figure 5.** Floating cloud.

*Exf* <sup>=</sup> <sup>1</sup> *h h* ∑ *j*=1 *pj* (13)

$$E n\_f = \sqrt{\frac{\pi}{2}} \times \frac{1}{h} \sum\_{j=1}^{h} \left| p\_j - E x\_f \right| \tag{14}$$

$$He\_f = \sqrt{\frac{1}{\hbar} \sum\_{j=1}^{\hbar} (p\_j - Ex\_f)^2 - (En\_f)^2} \tag{15}$$

where *pj* is the rating given by the *j*-th expert (1-point scale); *Exf*, *Enf*, and *Hef* are the expectation, entropy, and hyperentropy of the floating cloud, respectively.

At the same time, in view of the atomization nature of the cloud model (Equation (16)), it is suitable for the consensus judgment of group cognition and can be used as a critical condition in order to judge whether the scoring result of group experts is reasonable; and, if it is not satisfied, it will be re-scored:

$$\text{He}\_f \lhd \text{En}\_f / 3 \tag{16}$$

The greater the degree of overlap between the floating cloud and the *k*-th state cloud in the space (as shown in Figure 6), the stronger the correlation is between the two. Therefore, this paper uses the Formula (17) to calculate the overlapping area O*<sup>k</sup>* of the floating cloud and the *k*-th state cloud. In Figure 6, O2 and O3 are the overlapping areas of the floating cloud, the abnormal cloud, the floating cloud, and the attention cloud, respectively; the standardized O*<sup>k</sup>* is used as the membership degree of the qualitative index, as shown in Formula (18):

$$O\_k = \begin{cases} \int\_{-\infty}^{p\_0} \mu\_f dt + \int\_{p\_0}^{+\infty} \mu\_k dt & p\_0 \in \text{Ex}\_f\\ \int\_{-\infty}^{p\_0} \mu\_k dt + \int\_{p\_0}^{+\infty} \mu\_f dt & p\_0 \ge \text{Ex}\_f \end{cases} \tag{17}$$

$$r\_k^{(2)} = \frac{2O\_k}{\sqrt{2\pi}(En\_f + En\_k)}\tag{18}$$

where *rk* (2) is the membership of the qualitative indicator in the *k*th state, however, it should be noted that *rk* (1) and *rk* (2) are only the distinctions between the membership of quantitative and qualitative indicators and all subsequent use of *rk* to indicate the membership of an indicator in the *k*-th state level; *uf* and *uk* are the expectation curves of the *k*-th state cloud in the floating cloud and the qualitative space, respectively; *p*<sup>0</sup> is the intersection value of the two curves.

**Figure 6.** The degree of overlap between the floating cloud and the state cloud.

#### *4.5. Evaluation Algorithm Process*

In this paper, based on the principle of AHP algorithm, the evaluation system of mining cable operation status is established based on a combination of expert opinions and literature references, and the weight vectors of indicator layer and criterion layer are calculated according to the improved weight calculation method. According to the respective characteristics of quantitative and qualitative indicators, different methods are adopted to calculate the membership degree of the two types of indicators; finally, the indicator weights are realized according to the D-S evidence theory and membership fusion according to the D-S evidence theory. The specific process is as shown in Figure 7:

**Figure 7.** Evaluation algorithm flowchart.

Before conducting evidence fusion, it is necessary to standardize the identification framework of the identification object *θ* = {serious state, abnormal state, attention state, normal state, uncertain state} = {*At*}, *t* = 1, 2, 3, 4, 5; then, treat each indicator as evidence and construct the basic confidence assignment function with a generalized fuzzy number for it [26]. Considering that there are differences in the importance degree for each evidence, again, in this paper, we also add the importance factor to its correction, as shown in Equation (19):

$$m\_i(A\_t) = \begin{cases} \beta\_l \frac{r\_{ik}}{\sum\_{k=1}^4 r\_{ik} + (1 - \max\{r\_{ik}\})} & t \in [1, 4] \\\ t & 1 - \sum\_{t=1}^4 m\_i(A\_t) \\\ \end{cases} \tag{19}$$

$$
\beta\_{\dot{l}} = \lambda\_{\dot{l}} \frac{\omega\_{\dot{l}}}{\omega\_{\text{max}}} \tag{20}
$$

where *rik* is the membership degree of the *i*-th indicator in the *k*-th state; *mi*(*At*) is the trust degree of the ith indicator within *θ*; *β<sup>i</sup>* is the importance factor of the i-th evidence; *λ<sup>i</sup>* is the priority trustworthiness coefficient, usually taken as 0.9; and *ω*max is the maximum value of the combination weight.

Finally, according to the synthesis rules of evidence theory (Formulae (21) and (22)), the fusion is carried out step-by-step. In order to solve the problem of the poor status of the base layer indicators that cannot be shown in the final fusion results, this paper proposes the concept of a deviation coefficient *ζa*, as shown in Formula (23):

$$m(B\_t) = \frac{\sum\_{\substack{\cap A\_t = B \ t \le i \le 5}} m\_i(A\_t)}{1 - K} \tag{21}$$

$$K = \sum\_{\cap A\_l = \mathcal{D}} \prod\_{t \le i \le 5} m\_l(A\_t) \tag{22}$$

$$\mathfrak{f}\_a = |K\_{Ri} - K\_R| \tag{23}$$

where *KRi* and *KR* are the factor layer conflict coefficient and target layer conflict coefficient, respectively, which can be obtained from Equation (22); *ζ<sup>a</sup>* is the deviation coefficient of the *a*-th indicator in the factor layer, and if *ζ<sup>a</sup>* > 0.05, the state of the *a*-th indicator is taken as the final state of the cable, otherwise, the maximum value within the fusion result of the target layer is taken as the final state of the terminal.

#### **5. Numerical Work**

In order to verify the feasibility of the method proposed in this paper, we used the operational data of a 6 kV cable of a coal mine grid as the basis to determine the operational status of this cable and analyzed the superiority of the method proposed in this paper with the actual calculation results.

#### *5.1. Weight Vector Calculation*

Due to the limitation of space, the detailed numerical calculation of the weight vector calculation method is not presented in this paper (the detailed numerical calculation has been introduced in the previously published literature [22]). In order to save time in the calculation, we have written an arithmetic program for the weight calculation algorithm using MATLAB (see Supplementary Materials), into which we only need to input the interval judgment matrix reflecting the relative importance among the indicators provided by the experts to obtain the weight vector that we are seeking. The weight vectors of the indicators in each layer are calculated as follows:


*W* is the weight vector of each indicator in the criterion layer, and *W*1, *W*2, *W*<sup>3</sup> are the weight vectors of each indicator in the base layer, respectively.

#### *5.2. Fuzzy Evalution Matrix*

The quantitative and qualitative index membership of the base layer can be calculated according to Formulas (12) and (18), respectively, as shown in Table 3. In this paper, in order to show the feasibility and superiority of the requested qualitative index membership, the results of the membership of x31, x32, x33, and x34 under the operating condition information are compared using different methods, as shown in Table 4. Through Table 4, it is easy to find that the most proximate state levels of x31, x32, x33, and x34 are identical under the three methods, but the method in this paper is more convenient for us to visually determine the proximate state levels of the four indicators, as shown in Figures 8–11. In addition, the fuzzy statistics and the gray theory assign the possibility of x31, x32, x33, and x34 to four pre-set state levels completely, i.e., the sum of the membership degree of x31, x32, x33, and x34 is homogeneously normalized. Since the existing state levels are not carefully divided, there are deviations between the state space model and the actual model. Experts will express it with an accurate numerical value within the maximum possible range of a certain state and ignore the occurrence of the remaining states. In view of the above considerations, this paper uses the overlap between the floating cloud and the qualitative space in order to reduce the interference of the above factors and make the obtained results more conservative.


**Table 3.** Index membership of base layer.

**Table 4.** Membership solution results of x31, x32, x33, and x34 under different methods.


**Figure 8.** Membership degree of x31.

**Figure 9.** Membership degree of x32.

**Figure 10.** Membership degree of x33.

**Figure 11.** Membership degree of x34.

#### *5.3. Voltage Safety Level Judgement*

According to the relevant principles of D-S evidence theory, the relevant parameters of the factor layer indicators can be calculated by combining Equations (19) and (20), as shown in Table 5. Combined with Table 5, the target layer *m*(*B*) = {0, 0.0059, 0.1261, 0.6805, 0.1875} can be calculated by Formula (21). According to the principle of the maximum trust degree, it can be determined that the operating state of the cable is in a normal state and there is also a weak degree of trust in the attention state and abnormal state. Therefore, it is necessary to continue to analyze the conflict coefficient between the factor layer and the target layer. Taking the conflict coefficient of the target layer as the reference value and by comparing the conflict coefficients of x1, x2, and x3 with the reference value, the deviation coefficient of each index of the factor layer can be calculated by Formula (23), and it is found that the deviation coefficients are all less than 0.05, so there is no need to correct the judgment result.

**Table 5.** Factor layer related parameters.


#### **6. Conclusions**

The safe and stable operation of mining cables directly affects the safe production of coal mines. If we can accurately assess the operating status of mining cables, we can prevent problems before they occur and respond in time before accidents happen. In order to realize the accurate evaluation of the operating state of the mining cable, this paper proposes a comprehensive evaluation method based on AHP, membership cloud theory, and D-S evidence theory. Considering the ambiguity and randomness of the state level, this paper introduces the membership cloud theory to visualize the cable running state in the form of a cloud. When calculating the membership degree of the qualitative index, the overlapping degree of the floating cloud and the state space is used as the membership degree of the qualitative index which can more intuitively reflect the relationship between the expert score and the state level. In the fusion of index weight and membership degree, this paper replaces the traditional fuzzy synthesis algorithm with D-S evidence theory. The main conclusions are as follows:


**Supplementary Materials:** The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/s22197174/s1, Weight Calculation Program of AHP.

**Author Contributions:** Conceptualization, Y.W. and P.C.; methodology, P.C.; software, P.C.; validation, Y.S. and C.F.; writing—original draft preparation, P.C.; writing—review and editing, Y.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **New Technique for Impact Calibration of Wide-Range Triaxial Force Transducer Using Hopkinson Bar**

**Qinghua Wang 1, Feng Xu 2,\*, Weiguo Guo <sup>1</sup> and Meng Gao <sup>1</sup>**


**\*** Correspondence: xufeng@nwpu.edu.cn

**Abstract:** At the current stage, there is an urgent need for new techniques to dynamically calibrate advanced wide-range (up to 10<sup>4</sup> N~10<sup>5</sup> N) triaxial force transducers. Based on this background, this paper proposes a novel impact calibration method, specifically for the triaxial force transducer, with a wide range and high-frequency response. In this method, the Hopkinson bar, which is typically used to test the dynamic mechanical properties of materials, was used as a generator to generate reference input force for the transducer. In addition, unlike conventional methods, the transverse sensitivities of the transducer were given necessary importance in the proposed method. The calibration result of the triaxial force transducer was expressed in a sensitivity matrix, containing three main sensitivity coefficients and six transverse sensitivity coefficients. The least squares method (LSM) was used to fit the sensitivity matrix linearly. Calibration experiments were performed on a typical triaxial force transducer. Several key issues, involving the validity and the test range, of the method were further investigated numerically. The feasibility and validity of the method were eventually confirmed. The test range of the method can be up to 106 N.

**Keywords:** triaxial force transducer; wide-range; dynamic calibration; Hopkinson bar; sensitivity matrix

#### **1. Introduction**

Triaxial force transducer is a kind of sensing system which can detect and measure force quantity in three-dimensional space [1]. There are significant demands for a triaxial force transducer, having wide range, in many important applications, such as load testing for aircraft landing, crash testing of automobiles, etc. The development of the advanced wide-range triaxial force transducer is receiving increasing attention from researchers in sensors and other related fields [2]. In engineering, the triaxial force transducer needs to be calibrated before it can be put into service. Moreover, the force transducer used for dynamic force measuring ought to be calibrated using a dynamic calibration method, because a transducer calibrated by the static method will have larger measurement error in measuring dynamic forces [3]. However, the dynamic calibration of the wide-range force transducer is currently facing a major challenge in how to excite reference input forces with high amplitudes (10<sup>5</sup> N level and above) and narrow pulse widths (10<sup>2</sup> μs level and below).

Currently, the dynamic calibration methods for force transducers can be divided into three types: vibration calibration, impact calibration and step calibration. In the general vibration calibration method [4–7], the mounting end of the force transducer being calibrated is rigidly connected to a shaker, while the sensitive end of the force transducer is rigidly connected to an external mass. The external mass has the role of generating dynamic (sinusoidal) reference input force for the transducer by oscillating vertically with the shaker. Up to now, the upper limit of the calibration range of vibration calibration facilities can reach 10 kN [8]. In the general impact calibration method [3,9–14], an object of known mass is usually made to collide with the force transducer being calibrated. An impulse

**Citation:** Wang, Q.; Xu, F.; Guo, W.; Gao, M. New Technique for Impact Calibration of Wide-Range Triaxial Force Transducer Using Hopkinson Bar. *Sensors* **2022**, *22*, 4885. https:// doi.org/10.3390/s22134885

Academic Editor: Francesc Pozo

Received: 11 June 2022 Accepted: 27 June 2022 Published: 28 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

will be generated during the collision, and it will be input to the transducer as a reference force. Currently available impact calibration facilities can test up to 20 kN, but the pulse width of the generated reference force can only reach the millisecond level [10,11]. In the step calibration method [15], an object of known mass is first suspended above the force transducer being calibrated. Then, by releasing the object, the gravitational force of the object will be suddenly applied on the transducer as a step reference force. The test range of the step calibration device developed in literature [15] can only reach 102 N.

According to the above, it is difficult to generate a reference input force with an amplitude higher than 20 kN using the dynamic calibration methods available at present. Moreover, it is difficult to generate a reference force with pulse width shorter than the millisecond order using the existing methods. Therefore, the currently available methods cannot cope with the need for dynamic calibration of the advanced force transducer with a range up to 10<sup>5</sup> N and frequency response up to 10 kHz. On the other hand, in the calibration of multi-dimensional force transducers, the previous studies usually focused only on the main axis sensitivities of the transducer [4,5]. The transverse sensitivities, which may have an influence on measurement accuracy, have not received the attention they deserve.

In this paper, an impact calibration method specially for the triaxial force transducer, with wide range and high-frequency response, was proposed, based on the Hopkinson bar technique. Calibration experiments were then designed and conducted on a typical wide-range, high-frequency response triaxial force using the newly developed method. The sensitivity of the transducer was expressed in a matrix. The transverse sensitivities of the transducer were evaluated by the off-diagonal elements of this matrix. Finally, several key issues, involving the validity and the test range of the method, were investigated and discussed using the finite element method.

#### **2. Theory and Method**

#### *2.1. Least Squares Method*

Similar to the physical quantity of acceleration [16], the force in motion space can be considered as a vector quantity with both magnitude and direction. The physical process of detecting a three-dimensional force in motion space with a triaxial force transducer can be considered mathematically as the process of projecting a three-dimensional vector in force space to a three-dimensional vector in signal space. When the nonlinearities between the inputs and outputs of the transducer are not considered, the projection can be described as:

$$\begin{cases} \mathcal{U}\_X = \mathcal{S}\_{XX}F\_X + \mathcal{S}\_{XY}F\_Y + \mathcal{S}\_{XZ}F\_Z \\ \mathcal{U}\_Y = \mathcal{S}\_{YX}F\_X + \mathcal{S}\_{YY}F\_Y + \mathcal{S}\_{YZ}F\_Z \\ \mathcal{U}\_Z = \mathcal{S}\_{ZX}F\_X + \mathcal{S}\_{ZY}F\_Y + \mathcal{S}\_{ZZ}F\_Z \end{cases} \tag{1}$$

or, in vector form:

$$
\begin{bmatrix} \mathcal{U}\_X \\ \mathcal{U}\_Y \\ \mathcal{U}\_Z \end{bmatrix} = \begin{bmatrix} \mathcal{S}\_{XX} & \mathcal{S}\_{XY} & \mathcal{S}\_{XZ} \\ \mathcal{S}\_{YX} & \mathcal{S}\_{YY} & \mathcal{S}\_{YZ} \\ \mathcal{S}\_{ZX} & \mathcal{S}\_{ZY} & \mathcal{S}\_{ZZ} \end{bmatrix} \begin{bmatrix} \mathcal{F}\_X \\ \mathcal{F}\_Y \\ \mathcal{F}\_Z \end{bmatrix} \tag{2}
$$

The subscripts *X*, *Y* and *X* represent the three sensitive axes of the transducer. The value *Ui*(*i* = *X*, *Y*, *Z*) represents the output voltage of the *i*-axis, *Fi*(*i* = *X*, *Y*, *Z*) represents the reference input force of the *i*-axis and *Sij*(*i*, *j* = *X*, *Y*, *Z*) represents the sensitivity coefficient of the transducer, where the subscript *i* denotes the output axis, and the subscript *j* denotes the input axis. In particular, when *i* = *j*, *Sii* is referred to as a main axis sensitivity; besides, when *i* = *j*, *Sij* is referred to as a transverse sensitivity [16].

Considering *n* sets of linearly independent inputs and outputs of the triaxial force transducer, then Equation (2) becomes:

$$
\begin{bmatrix}
\mathcal{U}\_{X1} & \cdots & \mathcal{U}\_{Xn} \\
\mathcal{U}\_{Y1} & \cdots & \mathcal{U}\_{Yn} \\
\mathcal{U}\_{Z1} & \cdots & \mathcal{U}\_{Zn}
\end{bmatrix} = \begin{bmatrix}
\mathcal{S}\_{XX} & \mathcal{S}\_{XY} & \mathcal{S}\_{XZ} \\
\mathcal{S}\_{YX} & \mathcal{S}\_{YY} & \mathcal{S}\_{YZ} \\
\mathcal{S}\_{ZX} & \mathcal{S}\_{ZY} & \mathcal{S}\_{ZZ}
\end{bmatrix} \begin{bmatrix}
F\_{X1} & \cdots & F\_{Xn} \\
F\_{Y1} & \cdots & F\_{Yn} \\
F\_{Z1} & \cdots & F\_{Zn}
\end{bmatrix} \tag{3}
$$

or, in general matrix form:

$$\mathbf{U} = \mathbf{S} \mathbf{F} \tag{4}$$

where **F** (3 × *n*) is the matrix for the reference input force of the transducer. **U** (3 × *n*) is the matrix for the output voltage of the transducer. **S** (3 × 3) is the sensitivity matrix of the transducer. The diagonal elements in **S** are the main sensitivity coefficients of transducer. The off-diagonal elements in **S** are the transverse sensitivity coefficients of transducer.

In fact, the relationship from the input to the output of a multi-dimensional force transducer is generally not purely linear [17,18], which leads to a certain deviation from the actual output to its linear regression. Let

$$\mathbf{e} = \begin{bmatrix} \mathfrak{e}\_{X1} & \cdots & \mathfrak{e}\_{Xn} \\ \mathfrak{e}\_{Y1} & \cdots & \mathfrak{e}\_{Yn} \\ \mathfrak{e}\_{Z1} & \cdots & \mathfrak{e}\_{Zn} \end{bmatrix} \tag{5}$$

be the deviation matrix of the triaxial force transducer, then Equation (3) will become:

$$
\begin{bmatrix}
\mathcal{U}\_{\mathcal{X}1} & \cdots & \mathcal{U}\_{\mathcal{X}n} \\
\mathcal{U}\_{\mathcal{Y}1} & \cdots & \mathcal{U}\_{\mathcal{Y}n} \\
\mathcal{U}\_{\mathcal{Z}1} & \cdots & \mathcal{U}\_{\mathcal{Z}n}
\end{bmatrix} = \begin{bmatrix}
\mathbb{S}\_{\mathcal{X}X} & \mathbb{S}\_{\mathcal{X}Y} & \mathbb{S}\_{\mathcal{X}Z} \\
\mathbb{S}\_{\mathcal{Y}X} & \mathbb{S}\_{\mathcal{Y}Y} & \mathbb{S}\_{\mathcal{Y}Z} \\
\mathbb{S}\_{\mathcal{Z}X} & \mathbb{S}\_{\mathcal{Z}Y} & \mathbb{S}\_{\mathcal{Z}Z}
\end{bmatrix} \begin{bmatrix}
\mathcal{F}\_{\mathcal{X}1} & \cdots & \mathcal{F}\_{\mathcal{X}l} \\
\mathcal{F}\_{\mathcal{Y}1} & \cdots & \mathcal{F}\_{\mathcal{Y}l} \\
\mathbb{G}\_{\mathcal{Z}1} & \cdots & \mathcal{F}\_{\mathcal{Z}l}
\end{bmatrix} + \begin{bmatrix}
\varepsilon\_{\mathcal{X}1} & \cdots & \varepsilon\_{\mathcal{X}l} \\
\varepsilon\_{\mathcal{Y}1} & \cdots & \varepsilon\_{\mathcal{Y}l} \\
\varepsilon\_{\mathcal{Z}1} & \cdots & \varepsilon\_{\mathcal{Z}l}
\end{bmatrix} \tag{6}
$$

Equation (4) will become:

$$\mathbf{U} = \mathbf{S}\mathbf{F} + \mathbf{e} \tag{7}$$

Based on Equation (7), the sensitivity matrix **S** satisfying the least square principle is:

$$\mathbf{S}^\* = \underset{\mathbf{S}}{\operatorname{argmin}} ||\mathbf{S}||\_2^2 = \underset{\mathbf{S}}{\operatorname{argmin}} \mathbf{e}^T = \underset{\mathbf{S}}{\operatorname{argmin}} (\mathbf{U} - \mathbf{S} \mathbf{F}) (\mathbf{U} - \mathbf{S} \mathbf{F})^T \tag{8}$$

**S**∗ can be solved according to Equation (9).

$$\frac{\partial \left[ \left( \mathbf{U} - \mathbf{S} \mathbf{F} \right) \left( \mathbf{U} - \mathbf{S} \mathbf{F} \right)^{\mathrm{T}} \right]}{\partial \mathbf{S}} = 0 \tag{9}$$

The solution of **S**∗ is given in Equation (10).

$$\mathbf{S}^\* = \mathbf{U} \mathbf{F}^T (\mathbf{F} \mathbf{F}^T)^{-1} \tag{10}$$

Actually, what is really required in engineering is the inverse of sensitivity matrix **S**∗. As described in Equation (11), the force to be measured is obtained by multiplying the output of the transducer by the inverse of the sensitivity matrix **S**∗.

$$\mathbf{F} = \mathbf{S}^{\*\; -1} \mathbf{U} \tag{11}$$

#### *2.2. Principle of the Hopkinson Bar Method*

Figure 1 depicts the typical structure and the usage of the triaxial force transducers commonly used today. As can be seen, the structure of the transducer can be divided into three parts: the mounting end, the sensitive end and the main part. In practical use, the transducer is held on a platform by its mounting end. A specially designed adapter is usually required to capture and transfer the force to be measured to the sensitive end of the transducer. The sensitive axes of the transducer will then detect the force transferred to the sensitive end and output it as a three-channel voltage signal. The sensing components of transducers are encapsulated in the main part. In this paper, the sensitive axes of the triaxial transducer are represented by X, Y and Z, as shown in Figure 1a.

**Figure 1.** (**a**) Schematic diagram of the typical structure of a triaxial force transducer; (**b**) Schematic diagram of the assembly and the usage of a triaxial force transducer.

The principle of using a Hopkinson bar to calibrate the force transducer is schematically illustrated in Figure 2. Bullets with a taper at its front end are generally adopted to strike the Hopkinson bar. A shaper is attached to the front end of the bar for pulse shaping and to prevent the bar from being damaged by the impact of the bullet. The force transducer being calibrated is mounted on a block. The block constrains the displacement of the transducer in axial direction. The sensitive end (usually requiring an adapter) of the transducer is in good contact with the back end of the bar. The bullet strikes the shaper and excites a compressive wave (the incident wave) in the bar. The wave propagates along the bar to the right and can be detected by the strain gauge glued on the surface of the bar at approximately the middle length. According to the theory of one-dimensional elastic wave propagation, the incident wave will reflect at the interface between the bar and force transducer. The reflected wave can also be detected by the strain gauge. Assuming that the waves propagate along the bar with no dispersion and attenuation, then the strain history on the interface between the Hopkinson bar and transducer can be determined by superimposing the incident and reflected waves as:

$$
\varepsilon(t) = \varepsilon\_\mathrm{I}(t) + \varepsilon\_\mathrm{R}(t) \tag{12}
$$

where *ε*(*t*) is the strain history on the interface, *ε*I(*t*) and *ε*R(*t*) are the incident and reflected waves detected by the strain gauge, respectively. Thus, the reference input force of the transducer can be calculated as:

$$F(t) = EA[\varepsilon\_\mathrm{I}(t) + \varepsilon\_\mathrm{R}(t)]\tag{13}$$

where, *E* is the elastic modulus of the bar, *A* is the cross-sectional area of the bar.

**Figure 2.** Schematic diagram of the principle of using Hopkinson bar to calibrate force transducer.

For the triaxial force transducer, the principle of its calibration using the Hopkinson bar is shown in Figure 3. The transducer being calibrated was fixed on a mounting block. The block constrained the degrees of freedom of the transducer. A cubic adapter was adopted to capture and transfer the input force for the transducer. The sensitive axes of the triaxial force transducer were tested sequentially by the Hopkinson bar. The attitude of the transducer, as well as the position of the mounting block, had to be changed so that the direction of the force could always be aligned with the direction of the sensitive axis. It is worth noting that in the calibration, shown in Figure 3, only one of the three sensitive axes of the transducer had force input in each test. However, all the three sensitive axes

would have voltages output at the same time, due to the coupling between the sensitive axes. Therefore, in this particular case, there would be only one non-zero element in the column vector of the force matrix **F** in Equation (10), and the other two elements would both be zero. However, the elements in voltage matrix **U** would all be non-zero elements.

**Figure 3.** Schematic diagram of the principle of calibrating triaxial force transducer using Hopkinson bar.

#### **3. Experiment**

#### *3.1. The Triaxial Force Transducer*

The triaxial force transducer used in this paper was a B25B piezoresistive transducer, as shown in Figure 4, provided by AVIC (Aviation Industry Corporation of China, Ltd., Beijing, China). Table 1 lists the characteristic parameters of the B25B transducer. The main axis sensitivity coefficients of the transducer given in Table 1 were obtained using a static method. The transverse sensitivity coefficients of the transducer were unknown before this paper.

**Figure 4.** B25B triaxial force transducer from AVIC.


#### *3.2. Experiment Set-Ups*

The calibration system, established based on the Hopkinson bar, is shown in Figure 5. An air gun was adopted to launch the bullet. When valve 2 was closed and valve 1 was open, the high-pressure air chamber would be inflated by the air compressor. A barometer provided the real-time display of the air pressure in the chamber. When the pressure in the chamber reached an expected value, valve 1 should be closed to stop inflation. At this point, the bullet would be launched by the high-pressure air once valve 2 was opened. The bullet used was made of a 45# high-strength steel. Its geometry and dimensions are given in Figure 5. The velocity of the bullet could be adjusted by changing the pressure in the chamber at the moment of launching. The shapers used were made of a 2024 aluminum alloy. The dimension of the shaper was Φ15 × 5. The Hopkinson bar used was made of a 7075-aluminum alloy. The bar was a square bar with a section of 20 mm × 20 mm and a length of 1800 mm. The adapter used was a cube with a side length of 20 mm. It was connected rigidly to the transducer with a thread of M22× 15. The adapter, as well as the mounting block and its supporting, were all made of 45# high-strength steel.

**Figure 5.** The experiment set-ups based on Hopkinson bar.

The strain gauges used were BE120-3AA strain gauges from AVIC. Two strain gauges were symmetrically glued on the upper and lower surfaces of the bar. These two strain gauges were then connected to a Wheatstone bridge as its two opposing arms. In this way, the stress in the bar would then be converted to a voltage signal by the bridge. The voltage output from the bridge was amplified by a super dynamic voltage amplifier (SDY2107A from BDHSD Co., Ltd., Beidaihe, China; frequency response could be up to 2.5 MHz). The excitation voltages of the transducer were supplied by the power module of the amplifier through the Wheatstone bridge. The voltages output from the transducer were reversely input to the amplifier through the bridge and amplified by the amplifier. Each sensitive axis of the transducer required a signal channel. The voltages output from the amplifier would be captured and recorded by a high-speed digitizer (USB8502, a USB-driven card

from ART Technology, Co., Ltd., Beijing, China; sampling frequency up to 40 M/s). The digitizer was driven by an industrial personal computer (IPC). The sampling frequency in experiments was set to 10 M/s. In the experiments, the sensitive axes of the transducer were tested in the order of Z − X − Y.

#### **4. Results and Discussion**

#### *4.1. The Calibration Results*

In the experiments, the *Z*-axis of the transducer was tested by the bullet at pressures of 0.05, 0.1, 0.15, 0.2, 0.25, 0.3 and 0.35 MPa. The *X*-axis and *Y*-axis of the transducer were tested at pressures of 0.05, 0.075, 0.1, 0.125, 0.15, 0.175 and 0.2 MPa. The typical input forces and output voltages of the transducer are shown in Figure 6. Specifically, the forces and voltages in Figure 6a were recorded in the case where the *X*-axis was tested at the pressure of 0.15 MPa. The forces and voltages in Figure 6b were recorded in the case where the *Y*-axis was tested at the pressure of 0.15 MPa. The forces and voltages in Figure 6c were recorded in the case where the *Z*-axis was tested at the pressure of 0.3 MPa. It can be seen that the duration of the force pulses excited by the bullet was about 140 μs. In addition, the amplitude of the force pulse increased with the pressure, while the duration of the pulse seemed to be independent of the pressure. On the other hand, it can be seen from Figure 6 that, although only one of the three sensitive axes of transducer had force input, the other two axes would also have voltages output at the same time. This indicated that coupling effects existed between the sensitive axes of the transducer. Moreover, the coupling between *X*-axis and *Y*-axis was positive, and the coupling between *Z*-axis and *X*-axis, as well as the coupling between *Z*-axis and *Y*-axis, were both negative, as illustrated in Figure 6.

**Figure 6.** Typical input forces and output voltages of the transducer in the cases where (**a**) the *X*-axis was tested at pressure of 0.15 MPa; (**b**) the *Y*-axis was tested at pressure of 0.15 MPa; (**c**) the *Z*-axis was tested at pressure of 0.3 MPa.

The input force curves and the output voltage curves of the sensitive axes of the transducer had roughly the same trend. This demonstrated that the input forces applied to the transducer could be effectively detected by the sensitive axes. The data used to calculate the sensitivity matrix were the peak values of the force pulses and voltage pulses. The calculation result of the sensitivity matrix of transducer is given in Equation (14). It can be seen that the main axis sensitivity coefficients of the transducer obtained with the dynamic (impact) calibration method were different from those obtained with the static method (Table 1). More specifically, the sensitivity coefficients obtained with the dynamic method were smaller than those obtained with the static method. This coincided with the

result of a dynamic evaluation of a single-axis piezoresistive force transducer conducted previously [3].

$$\mathbf{S} = \begin{bmatrix} 0.440 & 0.042 & -0.027 \\ 0.037 & 0.460 & -0.032 \\ -0.039 & -0.063 & 0.211 \end{bmatrix} \tag{14}$$

#### *4.2. Discussion*

In this section, the numerical model of the calibration system was built using the commercially available finite element software ADAQUS (Version 2017, Dassault Aviation, Paris, France). A number of simulations were then conducted to investigate the validity and the calibration range of the newly proposed method. The typical model built is shown in Figure 7. Specifically, Figure 7a shows the numerical model, wherein the *Z*-axis was tested. Figure 7b shows the numerical model wherein the *X*-axis was tested. Figure 7c shows the numerical model wherein the *Y*-axis was tested. In addition, a typical cross section of the Hopkinson bar is shown in Figure 7d. The points P, Q, and R were the sampling points on the cross-sections. The surface of the adapter in contact with the bar is shown in Figure 7e. The center point C was the sampling point on this surface.

**Figure 7.** Typical numerical model of the calibration system when (**a**) *Z*-axis; (**b**) *X*-axis; (**c**) *Y*-axis was tested; (**d**) typical cross section of the Hopkinson bar; (**e**) the surface of the adapter which is in contact with the bar.

The materials of the components in numerical simulations were consistent with the settings in experiments. The bullet was assumed to be made of 45# steel. The shaper was assumed to be made of 2024 aluminum alloy. The bar was assumed to be made of 7075-T6 aluminum alloy. The material property constants used in simulations are all listed in Table 2. Since the shaper would deform plastically, due to the impact of bullet, a plastic model should be considered in addition to the elastic model in the material model of the shaper. The plastic model selected for the 2024 aluminum alloy was a Johnson-Cook model developed in literature [19]. The dimensional settings of the components in numerical simulations were also consistent with the settings in experiments.

The contact properties in the numerical models were all set to frictional contacts with a coefficient of 0.1. The threaded holes at the sensitive end and mounting end of the transducer were simplified into conventional holes with smooth inner surfaces. The main part of the transducer was simplified into a solid elastomer. The threaded bolt of the adapter was simplified into a cylinder with smooth outer surfaces. In the numerical models, the adapter was assembled into the transducer by using a tie constraint. The inner surface of the hole at the mounting end of the transducer was constrained in all the six degrees of freedom to simulate the configuration wherein the transducer was fixed on a mounting block. The type of the tetrahedral elements in the numerical models was set to C3D4 and the type of the hexahedral elements was set to C3D8R. The total number of elements of C3D4 and C3D8R in the numerical model was 127854.

**Table 2.** The material property constants used in the numerical simulations.


In order to verify the numerical models, we conducted experimental tests and numerical simulations using the same bullet and at the same impact velocities. Figure 8 shows the comparisons of the input forces of the transducer obtained from the experiments and the numerical simulations. Specifically, the force pulses in Figure 8a were obtained from the experimental test and numerical simulation wherein the *Z*-axis of the transducer was tested at the velocity of 20 m/s. The force pulses in Figure 8b were obtained from the experimental test and numerical simulation wherein the *X*-axis of the transducer was tested at the velocity of 12 m/s. It can be seen that the force pulses obtained from experimental tests and numerical simulations were in good agreement. So, the numerical models built could be considered valid and accurate to simulate the experimental tests. The slight differences between the amplitude and the duration of the force pulses obtained from experiments and simulations might be due to damping and defects of material [20], which were not taken into account in the numerical simulations.

**Figure 8.** Comparisons of the input forces of the transducer obtained from the experimental tests and the numerical simulations wherein (**a**) the *Z*-axis of the transducer was tested by the bullet at velocity of 20 m/s; (**b**) the *X*-axis of the transducer was tested by the bullet at velocity of 12 m/s.

#### 4.2.1. Validation of the Method

According to the principle described in Section 2.2, the input force of the transducer was obtained indirectly from the signal of the strain gauge. Therefore, the validity of the proposed method was predicated on the assumption that the impact force applied to the adapter could be accurately measured by the strain gauge. Let ෨ be the input force of the transducer calculated from the strain signals at the midpoint of bar; F be the actual force on the contact surface between bar and adapter, i.e., the actual input force of the transducer. According to the simulations, the stress distribution on the contact surface between bar and adapter could be considered uniform. So, F could be obtained by the product of the stress at point C (see Figure 7) and the cross-sectional area of the bar. The ෨ and F were compared in several cases, as shown in Figure 9. It could be seen that the ෨ curves were all in very good agreement with the F curves in all the cases. In fact, the relative differences between the durations of the pulses in Figure 9 were less than 2%. The relative differences between

the amplitudes of the pulses in Figure 9 were less than 1.5%. It could be considered that the actual input force of the transducer could be accurately measured by the strain signal at the midpoint of the bar. Therefore, the method could achieve the impact calibration of the transducer effectively by generating a measurable reference input force for it. It should be noted that if the structure and usage of the transducer being calibrated are different from Figure 1, the mounting form of the transducer should be changed to ensure the generated reference force can be effectively applied to the sensitive end of the transducer. The method will be valid and applicable, as long as the generated reference load can be effectively applied to the transducer sensitive end and be accurately measured.

**Figure 9.** The comparisons of ෨ and F in the cases wherein (**a**) the *Z*-axis was tested at the velocity of 20 m/s; (**b**) the *Z*-axis was tested at the velocity of 10 m/s; (**c**) the *X*-axis was tested at the velocity of 12 m/s; (**d**) the *Y*-axis was tested at the velocity of 10 m/s.

#### 4.2.2. Influence of Bullet Geometry

Figure 10 shows the waveforms of the input forces of *Z*-axis excited by the bullets with various geometries. The bullets I–V in Figure 10 were all 40 mm in length and all had a maximum diameter of 28 mm. The impact velocity of the bullets was 20 m/s. It can be seen that the waveform of the input force was related to the geometry of the bullet. To be more exact, it was related to the taper at the impact end of the bullet. The input force pulse excited by the bullet with a smaller taper at its impact end had a lower amplitude and a wider duration. Yet, the input force pulse excited by the bullet with a larger taper had a higher amplitude and a shorter duration. For an extreme case, the input force pulse excited by bullet I, a cylindrical bullet with a maximum taper at its impact end, had the highest amplitude and the shortest duration.

**Figure 10.** The waveforms of the input forces of *Z*-axis excited by the bullets with various geometries.

#### 4.2.3. Influence of Wave Propagation

In order to study the influence of wave propagation, a series of cross sections were selected at 10 cm intervals along the bar, as shown in Figure 7. The P, Q, R points were the sampling points on each section (Figure 7). Let *σp*, *σq*, *σ<sup>r</sup>* be the stress amplitude at the points P, Q, R, respectively. The values of *σp*, *σ<sup>q</sup>* and *σ<sup>r</sup>* on each of the selected sections were recorded. Figure 11 shows the typical trends of the *σp*, *σ<sup>q</sup>* and *σ<sup>r</sup>* along the bar. It can be seen the stress distribution on the impact end face of the bar had a significant dispersion. A certain distance was required for the stress to reach a uniform distribution across the section of bar. The distance in Figure 11 was about 160 mm, 8 times the diameter of the bar. Fortunately, the stress wave would propagate with almost a constant amplitude once the stress distribution reached a uniform state. Therefore, the attenuation of the stress amplitude caused by the wave propagation was negligible.

**Figure 11.** The trends of the *σp*, *σq* and *σr* along the bar in the case wherein bullet IV (also the bullet in Figure 5) was used.

#### 4.2.4. Calibration Range

Building on Equation (13), the reference force input to the transducer could also be described as:

$$F = A\sigma \tag{15}$$

where *A* is the sectional area of the bar, *σ* is the stress on the bar end face in contact with the adapter. According to the principle described in Section 2.2, the deformation of the bar during calibration should be elastic. So, the calibration range of the method was limited by the elastic limit of the bar. Let *σ<sup>e</sup>* be the elastic limit of the bar, then the input force of the transducer will satisfy the following condition.

$$F < A\sigma\_{\ell} \tag{16}$$

In the present study, the sectional area of the bar was 4×10−<sup>4</sup> <sup>m</sup>2, the elastic limit of the bar was about 330 MPa. Therefore, the upper limit of the calibration method would be up to 106 N when *σ* in Equation (15) reached 75.76% of the elastic limit of the bar. By adjusting the velocity and geometry of the bullet, a wide range from 104 N to 106 N could be achieved by the proposed method.

#### **5. Conclusions and Outlook**

In this study, a novel impact calibration method applicable for the triaxial force transducer, with wide range and a high-frequency response, was proposed. The Hopkinson bar was used as a force generator in this method. Calibration experiments were conducted on a typical wide-range triaxial force transducer. The transducer sensitivity was expressed in a matrix form to attend simultaneously to main axis sensitivities and transverse sensitivities. The main axis sensitivity coefficients were obtained using the proposed method and were then compared with the coefficients obtained with the static method. Finally, the validity and the test range of the method were investigated and discussed using a numerical method. The conclusions of this study can be summarized as follows:


This study provides an effective solution for the calibration of the wide-range triaxial force transducer at this stage. However, the greatest challenge for the calibration of the triaxial force transducer with wide range and high-frequency response is the generating of reference input forces synchronously along the three sensitive axes. This really is a challenge since it is hard to generate measurable high-amplitude force pulses and keep them synchronized within a short duration of tens of or hundreds of microseconds. How to achieve synchronous calibration of the triaxial force transducer deserves great attention from researchers in related fields. This will also be the focus of the authors' future work.

**Author Contributions:** Methodology, W.G.; Experiments, Q.W.; Data Processing, M.G.; writing original draft preparation, Q.W.; writing—review and editing, F.X.; supervision, F.X. and W.G.; funding acquisition, F.X. and W.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by National Natural Science Foundation of China (NSFC), grant number 11872051.

**Data Availability Statement:** The data that support the findings of this study are available from the corresponding authors upon reasonable request.

**Acknowledgments:** The authors would like to thank the Aviation Industry Corporation of China (AVIC) for providing the triaxial force transducer for free.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **A New Axial Stress Measurement Method for High-Strength Short Bolts Based on Stress-Dependent Scattering Effect and Energy Attenuation Coefficient**

**Tong Fu, Ping Chen \* and Aijun Yin**

College of Mechanical Engineering, Chongqing University, Chongqing 400044, China; fifaft@163.com (T.F.); aijun.yin@cqu.edu.cn (A.Y.)

**\*** Correspondence: chempion@126.com

**Abstract:** The accurate estimation of axial stresses is a major problem for high-strength bolted connections that needs to be overcome to improve the assembly quality and safety of aviation structures. However, the conventional acoustoelastic effect based on velocity-stress dependence is very weak for short bolts, which leads to large estimation errors. In this article, the effect of axial stress on ultrasonic scattering attenuation is investigated by calculating the change in the energy attenuation coefficient of ultrasonic echoes after applying axial preload. Based on this effect, a stressdependent attenuation estimation model is developed to measure the bolt axial stress. In addition, the spectrum of the first and second round-trip echoes is divided into several frequency bands to calculate the energy attenuation coefficients, which are used to select the frequency band sensitive to the axial stress changes. Finally, the estimation model between axial stress and energy attenuation coefficients in the sensitive frequency band is established under 20 steps of axial preloads. The experimental results show that the energy attenuation coefficient in the sensitive band corresponds well with axial stress. The average relative error of the predicted axial stress is 6.28%, which is better than that of the conventional acoustoelastic effect method. Therefore, the proposed approach can be used as an effective method to measure the axial stress of short bolts in the assembly of high-strength connections.

**Keywords:** bolt axial stress; acoustoelastic effect; scattering attenuation; stress-dependent attenuation coefficient; sensitive frequency band

#### **1. Introduction**

The aviation structure contains many precise bolt connection structures, and a loose bolt may lead to abnormal mechanical vibration or even local disintegration, thus causing serious aviation safety accidents. Controlling the magnitude of the axial force of the bolt is an important measure to improve the connection quality of the aviation structure and reduce the incidence of aviation accidents. As the aviation structure has high speed and high vibration characteristics, it is necessary to strengthen the accurate measurement of the bolt axial stress in the aviation structure to make the connection stable without breaking the ring and then improve the overall connection level. Therefore, accurate prediction of axial stresses and control of bolt preload is important to ensure machine functioning and structural stability. There are several methods for estimating axial stress and controlling bolt preload, including the torque method [1,2], the strain gauge method [3,4], and the ultrasonic method [5]. The torque method generates the expected tightening torque through manual torque wrenches or pneumatic, hydraulic wrenches to control the bolt preload, which depends on the fact that the torque applied to the bolt can be effectively converted into the preload of the bolt. However, only about 10–15% of the torque can be converted into axial preload due to the different friction coefficients between the bearing surfaces of the thread and the nut. Some reports have shown that the error in axial stress measurements based on

**Citation:** Fu, T.; Chen, P.; Yin, A. A New Axial Stress Measurement Method for High-Strength Short Bolts Based on Stress-Dependent Scattering Effect and Energy Attenuation Coefficient. *Sensors* **2022**, *22*, 4692. https://doi.org/10.3390/s22134692

Academic Editors: Yongbo Li, Bing Li, Jinchen Ji and Hamed Kalhori

Received: 10 May 2022 Accepted: 9 June 2022 Published: 22 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the torque method exceeds 30%. The strain gauge method uses the strain in the surface of the bolt to obtain its axial stress. This method is limited by the size of the measurement conditions, which makes it difficult to install the strain gauge on the stressed part of the bolt. The ultrasonic method has the advantages of high adaptability, accuracy, and stability, and this is the main development direction of axial stress measurement in the future.

Ultrasonic non-destructive methods, such as the piezoelectric ultrasonic method [6–8], the EMAT ultrasonic method [9–11], and the acoustoelastic effect method [12–17], have been reported to measure bolt axial stress. The piezoelectric ultrasonic method is primarily used for loosening the monitoring of bolted connections. Recently, the smart piezoelectric bolt is developed by embedding a piezoceramic transducer in the bolt head to improve measurement accuracy [18,19]. However, the devices used in this method are costly and complicated to install. As far as the EMAT method is concerned, a high-frequency generator is required. Furthermore, the response signal of the EMAT is easily disturbed by noise [20]. Among these methods, the acoustoelastic effect method is the most widely used. By adopting the acoustoelastic effect, ultrasonic longitudinal and transversal velocities can be used to characterize axial stress [15]. Nohyu Kim [12] used the velocity ratio of the mode conversion wave propagated in the bolt to estimate the axial stress of high-strength bolts. To improve the measurement accuracy, Yongmeng Liu et al. [13] proposed a method of measuring the fastening force with dry coupling. Dry coupling is used instead of liquid coupling to improve coupling conditions affected by the bolt surface. Qinxue Pan et al. [17] focused on the non-uniform distribution of axial stress in the effective stressed region of the bolt. The shape factor is proposed to eliminate the impact of the stress distribution on the propagation path of ultrasonic waves on the measurement. YASUI et al. [16] developed a calibration method based on the longitudinal and transverse waves for the short bolt. However, it needs to establish a large number of calibration curves to compensate for the difference in stress conditions in the actual stress measurement. Enxiao Liu et al. [14] investigated the influence of coupling layer thickness on the bolt axial stress. The result shows that the influence of the thickness change of the coupling layer on the accuracy of axial stress measurement is greater for short bolts than for long bolts. Apparently, the results of these methods are highly dependent on the accuracy of the time of flight (ToF) measurements, which can be significantly influenced by the bolt length. The principle of the conventional acoustoelastic effect method is to determine the axial stress by the difference in wave propagation time. Both the ultrasonic signal before and after bolt tightening should be measured precisely to obtain an accurate ToF. For the short bolts, the absolute axial elongation under the same load will be reduced accordingly, so that the difference in propagation time before and after applying axial preload will become difficult to distinguish, leading to invalid or large error measurement results. In this case, the acquisition of ToF is a great challenge for the employed AD digitizer. Therefore, the conventional acoustoelastic effect method is not ideal for the axial stress evaluation of short bolts.

In recent decades, ultrasonic scattering theory based on the interaction of microstructure with waves has been greatly developed, making attenuation-based methods an effective nondestructive evaluation tool. It has been widely used to extract microstructural parameters such as grain size, boundary, and dislocation [21–23]. In research on the ultrasonic attenuation method, the mechanism of scattering attenuation is highlighted [24–28]. Recently, Turner and Ghoshal [24] presented a theoretical basis to extract stress information from polycrystalline microstructure based on ultrasonic scattering attenuation. It was later extended and generalized by Kube [26]. The covariance tensor of elastic modulus variations was included as part of the previously developed ultrasonic grain scattering models and is proportional to the attenuation and backscatter coefficients [25,27,29]. Kube et al. [28] confirmed the stress dependence of the covariance tensor by investigating the change of the spatial variance amplitude under an applied uniaxial load. The results show that attenuation-stress variations are two to four times larger than the wave velocity-stress variations for metallic and molecular bonding materials. Therefore, the stress-dependent

attenuation effect is more suitable for high-strength short bolts that are not sensitive to the conventional acoustoelastic effect.

The stress-dependent scattering effect is based on the change in the attenuation of acoustic energy when ultrasound waves propagate through materials subjected to different uniaxial stresses. The relationship between stress and attenuation coefficient can be modeled with the help of the attenuation coefficient spectrum. Furthermore, the attenuation coefficient spectrum provides a detailed analysis of attenuation in the frequency domain, which can comprehensively show the relationship between stress, attenuation coefficient, and frequency. However, the main challenge is to obtain reliable frequency-dependent attenuation coefficients. In the traditional attenuation coefficient measurement method, the center frequency of the ultrasonic echo signal may shift downward due to the scattering attenuation effect of the material grain boundary, making it difficult to select a reasonable effective bandwidth range for further research [30], and the peak frequency and peak value only reflect the fundamental frequency components of the ultrasonic echo signal; the corresponding energy amount is not sufficient to reflect all the characteristics of the ultrasonic echo. Xiongbing Li et al. [31] established a multi-scale ultrasonic attenuation evaluation model based on wavelet transform, which uses attenuation coefficients at different scales to effectively control the errors arising from the fitting of the model itself. In addition, Min Li et al. [32] discovered that the frequency bands are sensitive to grain size variations based on the energy attenuation coefficient spectrum. The results show that the evaluation modes built using the energy attenuation coefficients of the sensitive bands have higher accuracy compared to other bands. Gao et al. [33] proposed a frequency bandwidth selection rule for amplitude ratio linear regression to determine the attenuation spectrum and they proposed a method to further improve the stability and accuracy. Based on the current research introduced above, the extraction of the spectral characteristic parameters confirms the fact that the ultrasonic scattering attenuation is sensitive to the frequency band, which means that the properties of the material make it attenuate differently in different frequency bands.

In this study, a novel bolt axial stress measurement method based on the stressdependent attenuation effect and energy attenuation coefficient is developed. The bandwidth of the measured echoes is divided into several frequency bands to select the frequency band sensitive to the axial stress changes. This article is organized as follows. First, the important concepts of stress-dependent elastic wave scattering theory are introduced. Second, the estimation model of axial stress on ultrasonic scattering attenuation is derived. The process of the calculation of the energy attenuation coefficient and the selection of the sensitive frequency band is presented. Following that, a series of experiments are then implemented to validate the proposed method. The results of axial stress measurements with different frequency bands and clamping lengths are analyzed. Furthermore, a comparison with the conventional acoustoelastic effect method was conducted to verify the accuracy and stability of the proposed method.

#### **2. Ultrasonic Measurement Model of Bolt Axial Stress Based on the Energy Attenuation Coefficient**

#### *2.1. Stress-Dependent Attenuation Theory*

In the models of Turner and Ghoshal [24] and Kube [26], the acoustoelastic of single crystals was derived by considering polycrystals as an ensemble of crystals and deriving the associated tensor of the effective stress-dependent elastic modulus. The stress-independent covariance tensor is defined as:

$$\Xi\_{ijkl}^{a\mathcal{B}\gamma\delta} = \left\langle \mathbf{C}\_{ijkl}^{\*} \mathbf{C}\_{a\mathcal{B}\gamma\delta}^{\*} \right\rangle - \left\langle \mathbf{C}\_{ijkl}^{\*} \right\rangle \left\langle \mathbf{C}\_{a\mathcal{B}\gamma\delta}^{\*} \right\rangle \tag{1}$$

where *C*∗ *ijkl* is the stress-dependent effective elastic modulus tensor of the polycrystalline materials. denote the ensemble average moduli. *C*∗ *ijklC*<sup>∗</sup> *αβγδ*! is performed as an average over all possible crystal orientations. For a homogeneous single crystal, it can be written as *Gijkl* = *Cijkl* + " <sup>−</sup>*CijklSrrpq* <sup>+</sup> *CijklSrlpq* <sup>+</sup> *CijrlSrkpq* <sup>+</sup> *CirklSrjpq* <sup>+</sup> *CrjklSripq* <sup>+</sup> *CijklmnSnmpq*# *σpq*, where *Cijkl*, *Cijklmn* are the second-order elastic modulus tensor and third order elastic modulus tensor, *Sijkl* = *C*−<sup>1</sup> *ijkl* is the elastic compliance tensor, and *σpq* is the stress tensor. For crystallites of cubic symmetry, the third-order elastic modulus tensor can be written as [24]

$$\mathbb{C}\_{ijkl} = \mathbb{C}\_{ijkl}^{I} + \upsilon \sum\_{n=1}^{3} a\_{in} a\_{jn} a\_{kn} a\_{ln} \tag{2}$$

where *<sup>ν</sup>* = *<sup>c</sup>*<sup>11</sup> − *<sup>c</sup>*<sup>12</sup> − <sup>2</sup>*c*<sup>44</sup> is the anisotropy coefficient, and *<sup>C</sup><sup>I</sup> ijkl* is the isotropic fourth-rank tensor. *aij* is rotation matrices with respect to the rotating coordinate system.

According to Equations (1) and (2), Kube derived a relationship between the attenuation coefficients and the anisotropy third-order elastic constants. The stress-dependent attenuation coefficient is defined as the scattered energy lost from the wave with polarization I scatter into polarization p [28]:

$$\mathfrak{a}\_{\mathbf{I}\rightarrow P} = \frac{\left(wl\right)^{4}}{2\rho^{2}V\_{\mathbf{I}}^{3}V\_{P}^{5}} \times \int\_{-1}^{1} \frac{\mathbf{h} \cdot \mathbf{I}^{\mathbf{I}\rightarrow S} \cdot \mathbf{g}^{\mathbf{T}}}{\left[1 + \left(wl\right)^{2}\left(V\_{\mathbf{I}}^{-2} + V\_{P}^{-2} - 2V\_{\mathbf{I}}^{-1}V\_{P}^{-1}\cos\theta\right)\right]^{2}} \times d(\cos\theta) \tag{3}$$

which is a general expression for the longitudinal and shear wave modes. *w* is the wave angular frequency, *l* is the mean diameter of the grains, *ρ* is the density of the material, and *VI* and *VP* are the velocities of the polarization I and the polarization *P*, respectively. *θ* is the scattering angle between the incident wave vector and I→*P* denotes the wave with polarization I scatters into polarization *P*. **h** = [1 T T2] and **g** = [1 cos2*θ* cos4*θ*] are the row-vectors of the stress magnitude and the scattering angle. The term **Γ**I→<sup>S</sup> in Equation (3) is a 3 × 3 matrix for the different possible scattering modes. For cubic symmetry crystallites, the total longitudinal matrices are:

$$\begin{aligned} \Gamma^{L \to L} &= \begin{bmatrix} 9 & 6 & 1 \\ & -9\zeta & -6\zeta & -\zeta \\ -9\zeta^2/4 & 3\zeta^2/2 & \zeta^2/4 \end{bmatrix} \\ \Gamma^{L \to SH} &= \begin{bmatrix} 5 & 5 & 0 \\ -5\zeta & -\zeta & 0 \\ 5\zeta^2/2 & \zeta^2/4 & -\zeta^2/4 \\ 10 & 1 & -1 \\ -10\zeta & -\zeta & \zeta \\ 5\zeta^2/2 & \zeta^2/4 & -\zeta^2/4 \end{bmatrix} \end{aligned} \tag{4}$$

where *ξ* = 2(*ν* + *η*)(3ᢖ*ν*) is a parameter related to the anisotropy constants and bulk modulus. ᢖ = (*c*<sup>11</sup> + *c*12)/3 is the bulk modulus of the cubic crystallite and *ν* = *c*<sup>11</sup> − *c*<sup>12</sup> − 2*c*<sup>44</sup> and *η* = *c*<sup>111</sup> − *c*<sup>123</sup> − 2*c*<sup>144</sup> − 4*c*<sup>155</sup> are second and third-order elastic anisotropy constants for cubic crystals, respectively. This assumes that the spatial correlation function *η*(*r*) = exp(−*r*/*l*) that described the two random points separated by the distance *r* is generally linear in the same grain. The total longitudinal attenuation constants are then:

$$
\mathfrak{a}\_L = \mathfrak{a}\_{L \to L} + \mathfrak{a}\_{L \to SH} + \mathfrak{a}\_{L \to SV} \tag{5}
$$

The subscripts *αL*→*L*, *αL*→*SH* and *αL*→*SV* represent the three scattering conversion modes of the incident longitudinal waves, respectively. Finally, in the Rayleigh scattering attenuation zone (*wl*2/*V*<sup>2</sup> *<sup>L</sup>* )1, the integration in Equation (3) for the longitudinal wave attenuation coefficient can be expressed as

$$a\_T^L = \frac{(2\pi f \ell)^4 \nu^2}{375 \rho^2 V\_L^8} \left( 8 + 12 \frac{V\_L^5}{V\_S^5} \right) \left( 1 - \frac{\zeta}{2} T \right)^2 \tag{6}$$

where *VL*, and *Vs* are the velocities of the longitudinal wave and the shear wave, respectively. Equation (6) builds a connection between the magnitude of the stress *T* and the strength of the ultrasonic attenuation. It should be noted that the formula is valid only under the condition that stress is not causing the polycrystal to yield.

#### *2.2. The Axial Stress Measurement Model Based on the Stress-Dependent Attenuation*

When an ultrasonic wave propagates through the polycrystalline material, its energy attenuates due to wave diffraction, scattering, and absorption by the interaction between the waves and the grains. In ultrasonic testing, the scattering-induced attenuation is the loss due to the perturbation of ultrasonic waves when the acoustic impedance between heterogeneous media does not match the interface (grain boundaries). Because polycrystalline materials are usually not viscoelastic, the absorption attenuation can be negligible in our study. In addition, the sidewalls of the rod component can suppress the diffusion of the acoustic beam by reflecting the dissipated energy. Therefore, in the bolt specimen, the attenuation of ultrasonic waves is essentially produced by the scattering process. To illustrate the dependency between stress and ultrasonic attenuation, a stress-free polycrystalline cylindrical specimen (with cubic crystal symmetry) subjected to tensile preload δ in the axial direction is considered, as shown in Figure 1. Based on the unified scattering theory [21], the attenuation of ultrasound waves in the MHz range propagating in polycrystals strongly depends on the scattering occurring at the grain boundaries. Due to the stress dependence of crystallites in their effective elastic characteristics, the attenuation coefficient, which results in the scaling of the wave amplitude during propagation, carries information about the stress and frequency properties.

**Figure 1.** Ultrasonic wave propagation in bolt under the axial stress.

We consider a one-dimensional model to calculate the response of ultrasonic waves to axial preload, and the axial stress distribution over the length of the bolt is shown in Figure 2. The stresses in the middle of the bolt are approximately uniform, whereas the stresses near the head and threads show a gradient pattern along the axial direction [17]. To simplify the problem, we assume that the full bolt length is the sum of the equivalent stressed length and the unstressed portion. The entire stressed part of the bolt is equivalent to a cylinder with uniformly distributed stresses, the length of which is defined as the equivalent stressed length. In Figure 2, *L*<sup>0</sup> is the full length of the bolt, *Le* is the equivalent stressed length, and *Lp* is the clamping length. Therefore, the attenuation coefficient of the bolt under axial preload can be expressed as:

$$
\overline{a}\_T = \frac{L\_\ell}{L\_0} a(\delta) + \frac{L\_0 - L\_\ell}{L\_0} a(0) \tag{7}
$$

where *α*(0) and *α*(*δ*) are the attenuation coefficients of the stress-free and equivalent stressed sections of the bolt. It should be noted that the effective stressed length *Le* is dependent on the axial preload, which means that *Le* will vary as the axial stress changes. Because the clamping length *Lp* is found to be linearly related to *Le* [17], the effective length ratio *β* can be used as a substitute for *Le*/*L*0. Substituting (6) into (7):

#2

4 *v*2 %

*L*

&

(8)

**Figure 2.** The one-dimensional axial stress model.

4 *v*2 %

*L*

&"

The effective length ratio *β* is expressed as

$$
\beta = \frac{kL\_p + b}{L\_0} \tag{9}
$$

where *k* and *b* are the coefficients with respect to the clamping length and axial preload. Equation (8) can be further simplified as follows:

$$\overline{\alpha}\_{T} = \beta A \left[ \left( 1 - \frac{\mathfrak{T}}{2} T \right)^{2} - 1 \right] + A = \beta A \frac{\mathfrak{T}^{2}}{4} T^{2} - \beta A \xi T + A \tag{10}$$

$$A = \frac{\left(2\pi f l\right)^4 v^2}{375\rho^2 V\_L^8} \left(8 + 12\frac{V\_L^5}{V\_S^5}\right) \tag{11}$$

where *A* is a coefficient for the model associated with material parameters *ρ*, *ν*, *VL*, *VS*, and *f* (the frequency in the employed frequency band). It can be seen that *A*, *β*, and ζ are stressdependent constants related to the applied stress preload and frequency. Equation (10) builds a connection between the magnitude of the bolt axial stress *T* and the attenuation coefficient.

*ζ* is negative for almost all metals [28]. Several observations can be derived from Equation (10). First, the attenuation coefficient *α<sup>T</sup>* will increase with the tensile preload. Secondly, *β* is positively correlated with the attenuation coefficient *α<sup>T</sup>* under the same

stress. Finally, because the ratio of *ζ* to *T* is large (*ζ*/*T* > 10<sup>3</sup> MPa), the scattering response is expected to be nearly linear for low stress (<500 MPa).

#### *2.3. Energy Attenuation Coefficients*

As discussed in the introduction, the main problem in establishing an estimation model using Equation (10) is to obtain a reliable frequency-dependent attenuation coefficient. On the one hand, it is desired to obtain the attenuation coefficient that can effectively reflect the relationship between attenuation and stress. On the other hand, it is desired to know whether the attenuation coefficients can better distinguish the change of axial stress in a specific frequency range. As shown in Figure 3, to model the relationship between the stress, attenuation, and frequency, the frequency spectrum energy is proposed to calculate the attenuation coefficient, and the calculation is performed in multiple frequency bands to find what is sensitive to axial stress changes. The detailed process is described as follows.

**Figure 3.** Flow chart of the calculation of energy attenuation coefficient.

During the propagation process, the ultrasonic waves will occur in a series of reflections and mode conversions at the bolt boundary. When the ultrasonic longitudinal transducer is located at the surface of the bolt head, part of the signal in the main flap region is unaffected by the boundary and reaches the end of the bolt directly. Therefore, the normally incident longitudinal-to-longitudinal ultrasound (travels one and two round-trips in the bolt of longitudinal wave polarization direction) is used to characterize scattering attenuation under uniaxial stress, as shown in Figure 3. A rectangular window is employed to intercept ultrasonic echoes traveling the first round-trip *v*1(*t*) and second *v*2(*t*)

round-trip of the longitudinal polarization direction and the echoes of the time domain are converted to the frequency domain, which can be expressed as:

$$S\_i(f) = \int\_{t\_1}^{t\_2} v\_i(t)e^{-j2\pi ft}dt\tag{12}$$

where *i* represents *i*th round-trip echo, *i* = 1, 2. Then, the frequency spectrum is divided into n bands within the effective bandwidth of the transducer, and each band has an equal bandwidth. The frequency spectrum of each band is Δ *Si*,*j*(*f*) , *<sup>j</sup>* represents the *<sup>j</sup>*th frequency band, *j* = 1,2 ... , n. It should be noted that n should be determined by the corresponding effective frequency bandwidth, and a high value of n does not necessarily achieve good results.

The spectrum energy *Ei;j*, by integrating each band Δ *Si*,*j*(*f*) is expressed as:

$$E\_{i,j} = \int\_{f\_j^b}^{f\_j^a} (\Delta |S\_{i,j}(f)|)^2 df \tag{13}$$

where *f <sup>t</sup> <sup>j</sup>* and *<sup>f</sup> <sup>b</sup> <sup>j</sup>* represent the top and bottom limits of each frequency band, respectively.

The spectrum energies of the *j*th frequency band of the first round-trip and second round-trip echoes, respectively, are obtained as *E*1,*<sup>j</sup>* and *E*2,*j*. Then, the energy attenuation coefficient of the *j*th frequency band can be expressed as:

$$\alpha\_{j} = \frac{10}{L\_{0}} \ln \left( E\_{1,j} / E\_{2,j} \right) \tag{14}$$

where *L*<sup>0</sup> is the length of the bolt. The energy attenuation coefficient in the n frequency band, *α*1, ... *αj*, ... *α<sup>n</sup>* can be obtained, which can reflect the relationship between attenuation coefficient *α* and frequency *f*.

In this study, 20 steps of axial preload were applied to the bolt and the signals collected under each axial preload step were converted into the frequency spectrum. Each frequency spectrum was then divided into n equal parts in the effective frequency band to obtain the energy attenuation coefficient matrix *α*20×*n*. One row of the matrix corresponds to the data under one axial preload step. Each row contains energy attenuation coefficients for n frequency bands. One column of the matrix contains the energy attenuation coefficients for 20 steps of axial preload in the same frequency band, denoted by *αj*. If *α<sup>j</sup>* exhibits obvious fluctuations compared to other frequency bands, then the energy attenuation coefficient in the *j*th frequency band is more sensitive to axial stress change. This frequency band is called the sensitive band and is denoted by *αsensitive*. To find the sensitive band, the coefficient of variation (CV) was employed to analyze the fluctuations of the energy attenuation coefficient matrix in n frequency bands [32].

$$\mathbf{CV} = \sigma\_{\mathfrak{u}} / \mu\_{\mathfrak{u}} \tag{15}$$

where *σα* and *μα* are the standard deviation and mean value of the energy attenuation coefficient in the *j*th frequency band, respectively. A high CV value implies a large fluctuation in the energy attenuation coefficient. Consequently, the estimation model between the axial stress and energy attenuation coefficient in the sensitive frequency band can be established. It should be noted that the dependence of the stress in the attenuation model comes from the quadratic form with the coefficient *ξ* and also from the dependence of the attenuation on the frequency itself; thus, *ξ* is not an indicator of the sensitivity of the material to the stress. The choice of the sensitive band is based on the attenuation of the material in the effective band of the transducer.

#### **3. Experimental Setup**

#### *3.1. Experimental Equipment*

The axial stress measurement experiment system is built as shown in Figure 4. The ultrasonic transmitting and receiving card JSR PRC50 (produced by JSR Ultrasonics) with a gain range from 14 dB to 60 dB is used to generate a high-voltage spike to excite the ultrasonic transducers. The AD-link PCIe-9852 oscilloscope (produced by AD-link Technologies) with a sampling rate of up to 200 MHz is selected as the data acquisition hardware. They are integrated into a portable industrial computer to form a DAQ system. To ensure stable coupling, the magnetic transducer PT7 (produced by the Dakota Ultrasonics of USA) with a 6 mm piezoelectric disc is used, and its parameters are shown in Table 1. The magnetic force brings the transducer in close contact with the end face of the bolt. The couplant B2 (produced by the Olympus Corporation) is applied on the bolt head surface to reduce the wave reflection in the bolt–transducer interface; the coupling layer thickness remains uniform and thin. Two coarse thread hex head bolts are used, and their detailed parameters are provided in Tables 2 and 3. To avoid echo distortion and improve the coupling condition, the end faces of each bolt were finished with a milling machine before the measurement experiments. Figure 5 shows the image of the polished surface from bolt specimens. The bolt is installed on the CTM2200s tensile testing machine, which is produced by Xieqiang Co. Ltd. of China, through a specially designed jig. The device will ensure that the bolt is applied with a uniaxial preload; the clamping fixture can provide a fixed force to ensure stable contact between the transducer and the bolt at each step.

**Figure 4.** Schematic of the experimental configuration.

**Table 1.** Parameters of the used ultrasonic transducer.


**Table 2.** Parameters of the bolt specimens.


**Table 3.** Specimens' material properties: density, ρ; bulk modulus, κ; second-and third-order singlecrystal anisotropy coefficients, ν and η.


**Figure 5.** The polished surface of the bolt specimens.

In addition, because the temperature is an important factor affecting the accuracy of bolt stress measurements [28], a heater and a thermometer with a resolution of 0.1 ◦C (JXB312, Berrcom Ltd., China) were used to control and monitor the measurement environment temperature. To avoid the influence of temperature changes on material attenuation and acoustic velocity, the entire experimental setup is located in a laboratory with a stable ambient temperature within ±1 ◦C by central air conditioning during the stress measurement tests.

#### *3.2. Experimental Processes*

According to Equation (10), the model parameters need to be calibrated before axial stress measurement. The calibrated process based on the energy attenuation coefficient is as follows:

(a) Install the bolt and attached the magnetic transducer on the surface of the bolt head and tighten the fixture device to ensure stable coupling between the transducer and the bolt before applying axial preload.

(b) Apply preload to the bolt to gradually increase its axial stress, and collect the ultrasonic echo signal data under the corresponding axial stress.

(c) Calculate the energy attenuation coefficient matrix *αm*×*<sup>n</sup>* based on the first roundtrip and second round-trip echo signals corresponding to each axial preload. Furthermore, obtain the energy attenuation coefficient in the sensitive frequency band *αsensitive* by selecting the highest CV value in n frequency bands.

(d) Obtain the α–T curve from the axial stress and *αsensitive* using a nonlinear least-squares fitting technique with Equation (10). Then, the calibrated parameters are determined.

After the calibration parameters are obtained, the axial stress can be measured. The measurement process is as follows:

(a) Tighten the bolt to produce axial stress *Tmea* and collect ultrasonic echo signal data under this axial stress.

(b) Obtain the energy attenuation coefficient in the sensitive frequency band *αmea*.

(c) Substitute the energy attenuation coefficient *αmea* into the estimation model based on Equation (10) to obtain the axial stress.

A flow chart of these steps is shown in Figure 6.

**Figure 6.** Modeling and predicting process.

#### **4. Analysis of Experimental Results**

To verify the validity of the proposed method, the following experiments were conducted. First, the bolt axial preload process was divided into 20 steps; each step of the tensile testing machine is set as 10 MPa, the range of loading is 0–200 MPa, and ultrasonic signals after each step were measured. Next, the energy attenuation coefficient matrix *α*20×*<sup>n</sup>* of the 20 steps was obtained to select the sensitive frequency band. The estimation model was built using the calibrated parameters. Thereafter, the results of other frequency bands and different clamping lengths were analyzed. Finally, experiments based on the conventional acoustic elastic method were carried out to compare them with the proposed method.

#### *4.1. Analysis of the Ultrasonic Frequency Spectrum*

Figure 7 shows the changes of the first and second round-trip echo amplitude under 20 steps of axial preload for the clamping length of 35.70 mm. It can be clearly seen that the

echo amplitude decreases with the increase in axial preload. This phenomenon indicates that the axial preload enhances the scattering effect and leads to more severe attenuation.

**Figure 7.** The echo signals under 20 steps of axial preload for the clamping length of 35.70 mm.

Rectangular windows were used to extract the first and second round-trip echo signals, and then the fast Fourier transform (FFT) was used to analyze their frequency-domain characteristics. Figure 8a,b show the corresponding spectrum obtained from the FFT results. It can be seen that the center frequencies of the first and second round-trip echo decrease as the axial preload increases. The center frequency of the first round-trip echo was reduced from 8.2 to 7.9 MHz, and the second round-trip echo was reduced from 7.8 to 7.6 MHz, which is lower than the nominal value. This result is due to the fact that the attenuation of the high-frequency content is more severe than that of the low-frequency content, and the scattering attenuation effect is enhanced with the increase in the stress, which makes this phenomenon more obvious. In addition, it can also be found that the center frequency of the second round-trip echo is further reduced and moved to the left along the frequency spectrum compared to the first round-trip echo. The maximum center frequency of the first round-trip echo signal is close to 8.2 MHz. The center frequency of the second round-trip echo signal is further reduced with a maximum value of 7.8 MHz. The reason is that the second round-trip echo travels a longer distance and encounters more grain boundaries than the first round-trip echo.

It should be emphasized that the spectrum behavior of both the first and second round-trip echoes should be considered. From the frequency domain characteristics shown in Figure 8, the frequency range within the full width of the half-maximum of these first round-trip echoes is approximately 5.5–11 MHz, whereas the range of the second roundtrip echo signal is approximately 5–10.5 MHz. Based on overall considerations, the range between 5.5 and 10.5 MHz was chosen as the effective bandwidth for calculating the energy attenuation coefficient.

**Figure 8.** The frequency spectrum of the echo signals under 20 steps axial preload for clamping length of 35.70 mm: (**a**) frequency spectrum of the first round-trip echoes; (**b**) frequency spectrum of the second round-trip echoes.

#### *4.2. Selection of Sensitive Frequency Bands*

To improve the prediction accuracy of axial stress, it is necessary to select the sensitive frequency band energy attenuation coefficient to build the estimation model. The effective spectrum bandwidth in the range of 5.5–10.5 MHz is divided into five parts: 5.5–6.5 MHz, 6.5–7.5 MHz, 7.5–8.5 MHz, 8.5–9.5 MHz, and 9.5–10.5 MHz. The energy attenuation coefficients are calculated for each axial preload step in each frequency band. Thus, the energy attenuation coefficient matrix *α*20×<sup>5</sup> can be obtained. One matrix row contains the energy attenuation coefficients of five frequency bands, and one matrix column contains the energy attenuation coefficients of 20 axial preload steps in the same frequency band.

The CV values for each band were calculated using Equation (15). The results are shown in Figure 9. It can be seen that the CV value of the energy attenuation coefficient in the 8.5–9.5 MHz band is the highest. This indicates that the fluctuation of the energy attenuation coefficient is the largest in this band and is highly sensitive to the stressdependent attenuation changes. It is worth noting that the energy attenuation coefficient is calculated using the spectrum energy of the first and second round-trip echoes, which increases the CV value in the frequency band compared to the frequency peak-based calculation method. In this study, the frequency band of 8.5–9.5 MHz was chosen as the sensitive band to build the estimation model based on Equation (10).

The polynomial least squares fitting is used to model the measurement data. The correlation index *R*<sup>2</sup> was employed to evaluate the validity of the model in the modeling process. The range of *R*<sup>2</sup> is [0, 1]; a larger *R*<sup>2</sup> indicates a better approximation capability of the proposed model. Figure 10a shows the calibrated curve in the sensitive frequency band 8.5–9.5 MHz. The correlation coefficient *R*<sup>2</sup> of the estimation model is 0.9951. It is worth noting that the energy attenuation coefficient is positively correlated with the increase in axial stress.

**Figure 9.** Value of CV in each frequency band.

**Figure 10.** Fitting curves: (**a**) the calibrated curve of L6 = 35.70 in 8.5–9.5 MHz; (**b**) the calibrated curves of L1-L6 in 8.5–9.5 MHz.

In the proposed method, the calculation of the axial stress depends mainly on the energy attenuation coefficient of the ultrasonic echo in the stress region. According to Equation (10), the attenuation of ultrasonic echoes is also influenced by the clamping length, through the 20 steps of the axial preload calibrated process for six different clamping lengths. The relationship between clamping length and echo energy attenuation is illustrated in Figure 10b. It can be seen that the energy attenuation coefficient increases with the increase in the clamping length under the same axial preload. This is because the increasing clamping length will correspondingly increase the propagation distance of the ultrasonic echoes in the effective stress region, which enhances the scattering attenuation.

#### *4.3. Measurement Results in Other Bands*

In the process of axial stress estimation, the relative error is used to evaluate the accuracy of the prediction. A small relative error indicates the high accuracy of the prediction results. We further validate the effectiveness of the proposed method. Table 4 shows the prediction results for the sensitive frequency band and entire frequency band, and Table 5 shows the prediction results for the other four frequency bands.


**Table 4.** Prediction results in the sensitive and entire frequency bands.

**Table 5.** Prediction results in the other four bands.


The comparison between Tables 4 and 5 shows that the minimum relative error in the 8.5–9.5 MHz band is only 2.35%. The relative errors for the other four frequency bands in Table 5 show large fluctuations. For example, in the results of the 8.5–9.5 MHz band, the relative error of 115.01 MPa is 2.70%. However, the relative errors of the 6.5–7.5 MHz band and 7.5–8.5 MHz band are both above 4.5%. Similarly, the relative errors of the prediction results based on all frequency bands are all larger than that of the sensitive frequency band. Furthermore, the correlation index *R*<sup>2</sup> of 8.5–9.5 MHz is the closest to 1 of all the other bands.

Three clamping lengths of 34.43 mm, 33.42 mm, and 32.93 mm were selected to verify the sensitive frequency band. The actual axial preloads were re-selected, and eight axial stress measurements were carried out under each clamping length. The average relative errors for the five frequency bands and the entire band are shown in Figure 11. The results show that the average relative error in the 8.5~9.5 MHz band is also the smallest, which is approximately 7%. The average relative errors of other frequency bands are more than 9%. Therefore, selecting the energy attenuation coefficient in the sensitive band to build the estimation model is effective and necessary, which can better reflect the effect of axial stress on ultrasonic scattering attenuation.

**Figure 11.** Average relative errors of three different clamping lengths.

#### *4.4. Comparison of the Proposed Method and Conventional Acoustoelastic Method*

Next, we tested the accuracy and repeatability of the proposed method in practical application. First, 20 steps of axial preload were performed with bolt A1 to build the estimation model. Second, 10 steps of axial preload were applied to bolt A2 and the energy attenuation coefficient in the sensitive band was calculated for each step. Then, the axial stress of the same type of bolt A2 was measured by substituting the calculated energy attenuation into the estimation model. Repeating the measurement experiment under three different clamping lengths, the measurement results are shown in Table 6.


**Table 6.** Measurement results of axial stress based on the energy attenuation coefficient method.

Figure 12a shows the predicted relative error and the average relative error in the three clamping lengths. The average relative error is 6.28% and the variation trend of the predicted relative error in the three clamping lengths is relatively stable. It should be noted that the predicted results for the low-stress region (<80 MPa) have a larger error relative to the high-stress region (>80 MPa) in all clamping lengths. The relative errors of the high-stress region (>80 MPa) are less than 6%. However, the minimum relative error of the low region (<80 MPa) is 9.57%. On the one hand, the stress-dependent attenuation effect becomes weaker as the axial stress decreases. On the other hand, the nonlinearity of the stress distribution is enhanced in the low-stress region, which reduces the accuracy of the estimated model. For high-strength short bolts (working stress > 100 MPa), this method has a good measurement accuracy.

**Figure 12.** Relative error: (**a**) relative error of the proposed method; (**b**) relative error of the conventional acoustoelastic effect method.

Axial stress measurement experiments based on the conventional acoustoelastic method were also conducted for the same three clamping lengths. Thirty sets of axial stresses were measured and the results are shown in Table 7. Compared with the three clamping length evaluation results of the two methods, the error band of the two methods can be observed in Figure 12. The average relative errors of the axial stress measured by the stress-dependent attenuation method and the conventional acoustoelastic effect method were 6.28% and 16.68%, respectively. Furthermore, the error bands of the predicted values were reduced significantly, from 8.8–5.9% of the attenuation method to 30.78–3.04% of the conventional method, which means that the measurement stability of the attenuation method is better.



Therefore, the attenuation energy ultrasonic method is more suitable for the axial stress measurement of high-strength short bolts, which makes the prediction results more accurate and stable. Furthermore, the following advantages of the stress-dependent attenuation method in practical applications give it better application prospects. First, compared with the conventional acoustoelastic effect, the stress-dependent attenuation effect is more sensitive to the stress changes, and the energy attenuation coefficient in the sensitive frequency band benefits from distinguishing the subtle changes of stress and improving the accuracy of prediction. Second, compared with the time-domain characterization, the spectrum analysis has better noise immunity. Third, the stress-dependent attenuation method does not need to measure the signal characteristics under zero stress, which improves the measurement efficiency and availability.

#### **5. Conclusions**

This study developed an axial stress measurement method for short bolts based on the stress-dependent scattering effect and energy attenuation coefficient. The estimation model of axial stress on the ultrasonic scattering attenuation is established based on the energy attenuation coefficient in the sensitive band. The energy attenuation coefficient in different frequency bands can show a more comprehensive analysis of attenuation in the frequency domain. Compared with other frequency bands, the axial stress estimation model based on the energy attenuation coefficient in the sensitive frequency band is more accurate. Ultrasonic experiments were performed with a magnetic transducer for the 45steel short bolt M10\*54. The experimental results show that the energy attenuation coefficient in the sensitive band is well correlated with the axial stress of the bolt. Compared with the conventional acoustoelastic effect method based on the change of propagation time before and after bolt tightening, the proposed method has higher prediction accuracy and better stability. In addition, the method proposed in this article has not only been experimentally proven to be feasible for evaluating axial stresses in bolts but it is also applicable to validate some other similar shaft-like bars, and thus is well suited for structural health monitoring systems.

It is worth mentioning that the sensitive frequency band in this study is based on testing data acquired from the 45steel short bolt. The effective sensitive frequency band needs to be determined by the selected bolt material properties and transducer characteristics. Furthermore, the scattering assumption in this study is based on the equiaxial grain shape. When the ultrasonic wave passes through the material with irregular grain shapes, the model may produce large relative errors. In addition, the effect of non-uniform stress distribution on the estimation and the uncertainties in the model and experiment will be further investigated in future work.

**Author Contributions:** T.F. conceived this article, perfomed the investigation, and wrote the original draft; P.C. data curation, and validation; and A.Y. revised the article and provided some valuable suggestions. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (grant number 51374264) and the Science and Technology Major Project of Chongqing (grant number cstc2018jszx-cyztzxX0032).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **A Crack Propagation Method for Pipelines with Interacting Corrosion and Crack Defects**

**Mingjiang Xie 1, Yifei Wang 1, Weinan Xiong 2, Jianli Zhao <sup>1</sup> and Xianjun Pei 1,\***


**Abstract:** Corrosion and crack defects often exist at the same time in pipelines. The interaction impact between these defects could potentially affect the growth of the fatigue crack. In this paper, a crack propagation method is proposed for pipelines with interacting corrosion and crack defects. The finite element models are built to obtain the Stress Intensity Factors (SIFs) for fatigue crack. SIF interaction impact ratio is introduced to describe the interaction effect of corrosion on fatigue crack. Two approaches based on extreme gradient boosting (XGBoost) are proposed in this paper to predict the SIF interaction impact ratio at the deepest point of the crack defect for pipelines with interacting corrosion and crack defects. Crack size, corrosion size and the axial distance between these two defects are the factors that have an impact on the growth of the fatigue crack, and so they are considered as the input of XGBoost models. Based on the synthetic samples from finite element modeling, it has been proved that the proposed approaches can effectively predict the SIF interaction impact ratio with relatively high accuracy. The crack propagation models are built based on the proposed XGBoost models, Paris' law and corrosion growth model. Sensitivity analyses regarding corrosion initial depth and axial distance between defects are performed. The proposed method can support pipeline integrity management by linking the crack propagation model with corrosion size, crack size and the axial distance. The problem of how the interaction between corrosion and crack defects impacts crack defect growth is investigated.

**Keywords:** pipeline; fatigue crack; corrosion; stress intensity factor; finite element; XGBoost

#### **1. Introduction**

Pipelines are widely used to transport oil and gas products over long distances. Ensuring pipeline safety is a prerequisite for the transportation of fuels such as oil and natural gas. Researchers are committed to constructing more accurate and effective health management models and improving the integrity management system of pipelines. Researchers [1–3] summarized the existing models in the field of pipeline integrity management and pointed out that, although the current models consider the accuracy of inline inspection tools, they are still too ideal and challenging to accurately reflect the proper working conditions of the pipeline. Metal-loss corrosion defects are significant threats to pipeline integrity. Some researchers use stochastic processes to describe uncertainties associated with the degradation of wall thickness incurred by corrosion defects. Wang et al. [4] proposed a stochastic corrosion growth model using the geometric Brownian bridge process. Ossai et al. [5] used a non-homogeneous linear growth pure birth Markov model to predict the degradation of internal corrosion defects in oil and gas pipelines. Bazan and Beck [6] employed a Poisson square wave process to describe the corrosion growth rate and compared the proposed non-linear stochastic model with the linear corrosion growth model. Qin et al. [7] proposed a corrosion growth model based on Inverse Gaussian process and Markov Chain Monte Carlo simulation method. Pan et al. [8] also used Inverse Gaussian process to characterize

**Citation:** Xie, M.; Wang, Y.; Xiong, W.; Zhao, J.; Pei, X. A Crack Propagation Method for Pipelines with Interacting Corrosion and Crack Defects. *Sensors* **2022**, *22*, 986. https://doi.org/ 10.3390/s22030986

Academic Editors: Hamed Kalhori, Yongbo Li, Bing Li and Jinchen Ji

Received: 10 January 2022 Accepted: 26 January 2022 Published: 27 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the degradation process of defects. Peng et al. [9] proposed a Bayesian framework of Inverse Gaussian process models. Remaining useful life of pipelines with multiple defects was predicted in refs. [10–12]. Although these corrosion growth models take multiple corrosion defects into account, they hardly consider the interacting effects among these defects, let alone the interactions between different types of defects.

There are a number of papers investigating pipelines with interacting corrosion defects. Benjamin et al. [13,14] presented a detailed literature review of pipelines with interacting corrosion defects and a database of corroded pipe tests. Amandi et al. [15] proposed a finite element model combined with a curve fitting method to estimate the remaining strength of pipelines with interacting corrosion defects. Sun and Cheng [16] also implemented a 3D finite element model to investigate mechano-electrochemical interaction of multiple longitudinal corrosion defects. Soares et al. [17] presented a model to analyze the integrity of pipelines with interacting corrosion defects under internal pressure and thermal stresses. Chen et al. [18] used a nonlinear finite element model to study the failure pressure of X80 pipelines with interacting corrosion defects. Kuppusamy et al. [19] investigated the effect of interaction of corrosion defects on the buckling strength of pipelines. Assessing and managing crack defects is also a vital part of pipeline integrity management. The remaining useful life prediction for pipelines with a single crack defect was conducted in refs. [20–22]. As for pipelines with interacting crack defects, Zhang et al. [23] presented a numerical model and fatigue simulations to analyze the fatigue behaviors. The corrosive environment will affect the growth of the crack, which is called Stress Corrosion Cracking (SCC). Hu et al. [24] applied the Monte Carlo method to predict and evaluate SCC. Lu et al. [25] established an SCC crack growth model in a high pH environment and verified it through experiments. Sekhar [26] summarized the effects of various crack interactions. This study shows that it is necessary to include the analysis of the interaction coupling between crack defects. These studies are all about the interaction between different defects of the same type. However, the exploration of interaction impact between different types of defects is still lacking in the existing literature.

In pipelines, common pipeline defects, such as crack and corrosion, exist at the same time. Specifically, there is an interacting effect between the fatigue crack and corrosion defect in the same pipeline segment. Pipeline corrosion will change the strength of the pipeline in the surrounding area. If the corrosion and crack defects are adjacent, a certain interaction coupling will occur and impact the Stress Intensity Factor (SIF) of the crack surface, thereby affecting the propagation of fatigue crack. Therefore, the crack propagation model that considers the interaction between these two types of defects is conducive to formulating more accurate detection and maintenance strategies. Motivated by this need, this paper plans to study the interacting effects of corrosion and crack defects on pipeline crack propagation.

In this paper, a method was developed for predicting the propagation of fatigue crack for pipelines with interacting corrosion and crack defects. Crack length, crack depth, corrosion length, corrosion depth and the axial distance between the crack and corrosion defects are all considered when developing this method. The finite element models are built to obtain the SIF values with and without considering the interaction impact between these defects. The powerful regression model, XGBoost [27], is applied in this paper to predict the SIF interaction impact ratio. With synthetic data from finite element analysis modeling, two approaches are provided to fit and predict the SIF interaction impact ratio at the deepest point of the crack defect, considering the interaction between corrosion and crack defects. The first one uses the data samples to directly fit and predict the SIF interaction impact ratio with a XGBoost model. As for the second one, it is an indirect prediction approach. It fits the SIF with and without considering the interaction impact, respectively. Therefore, two XGBoost models are acquired in this approach. The prediction results from these two models are utilized to calculate SIF interaction impact ratio. SIF interaction impact ratio is defined as the ratio of the SIF considering the interaction impact divided by the SIF without considering the interaction impact. With the proposed XGBoost

models and traditional crack and corrosion growth models, a crack propagation model is proposed for pipelines with interacting crack and corrosion defects, and simulation results are obtained for sensitivity analysis.

The novelty of this paper is three-fold: (1) it studies the interaction impact between different types of defects in pipelines, viz. crack and corrosion defects, depending upon the crack size, corrosion size and the axial distance between them; (2) it introduces SIF interaction impact ratio to describe the degree of the interaction impact and employs an advanced machine learning algorithm XGBoost to fit and predict the SIF interaction impact ratio; and (3) it proposes a method to predict the propagation of fatigue crack considering the interaction impact.

The rest of the paper is organized as follows. Section 2 presents the finite element analysis model for a pipeline with interacting corrosion and crack defects. Section 3 presents the proposed crack propagation model based on XGBoost. In Section 4, experimental results are obtained to analyze the interaction impact. Conclusions are presented in Section 5.

#### **2. The Pipeline Finite Element Analysis Model**

#### *2.1. The FEA Model*

In this section, the finite element software ANSYS*®* is used to model the pipeline with interacting fatigue crack and external corrosion defects. In the modeling process, the pipeline models with and without corrosion defects are established, respectively, to analyze the interaction impact of corrosion defect on crack propagation. The material of the modelled pipeline is API 5L X70. The outside diameter of the pipeline is set as 914.4 mm, and the wall thickness is 15.875 mm. The internal pressure is assumed to be 1 MPa for modeling. The fatigue crack is modeled as a semi-elliptical shape with a length of 15.2 mm and a crack depth in the range of 2 mm–12 mm. The SIF values corresponding to the deepest point and edge point can be obtained through stress analysis. The internal pressure of the pipeline is 1 MPa. At the same time, there are cuboid corrosion defects on the outer surface of the pipeline, and the axial distance from the crack center to the corrosion center moves from 150 mm to 500 mm. The depth of the corrosion defect is from 2 mm to 14 mm, with an increment of 1 mm each time. The geometric modeling of a corroded pipeline is shown in Figure 1, and the finite element model built in this paper is shown in Figure 2.

**Figure 1.** Geometric modeling of corroded pipeline.

**Figure 2.** The grid division of pipeline. (**a**) Partial refinement of corrosion defect grid. (**b**) Partial refinement of grid at crack.

#### *2.2. Validation of FEA*

Generally, fracture in engineering structures can be classified into three types: opening mode (I), sliding mode (II) and tearing mode (III), and SIF is used to reflect these modes. Compared with mode II and III, SIF corresponding to mode I is much larger, so the mode I SIF dominates the propagation of fatigue crack. In this paper, mode I SIF was only considered in the pipeline remaining useful life prediction. The method based on API579 for the partial verification of FEA model was employed. According to API579 criterion, the SIF of mode I of the pipeline is calculated as follows:

$$K = \frac{pR\_i^2}{R\_o^2 - R\_i^2} \left[ 2G\_0 + 2G\_1 \left(\frac{a}{R\_o}\right) + 3G\_2 \left(\frac{a}{R\_o}\right)^2 + 4G\_3 \left(\frac{a}{R\_o}\right)^3 + 5G\_4 \left(\frac{a}{R\_o}\right)^4 \right] \sqrt{\frac{\pi a}{Q}} \tag{1}$$

$$Q = 1.0 + 1.464 \left( \frac{a}{c} \right)^{1.65}, \text{ for } \frac{a}{c} \le 1 \tag{2}$$

$$G\_0 = A\_{0,0} + A\_{1,0}\boldsymbol{\beta} + A\_{2,0}\boldsymbol{\beta}^2 + A\_{3,0}\boldsymbol{\beta}^3 + A\_{4,0}\boldsymbol{\beta}^4 + A\_{5,0}\boldsymbol{\beta}^5 + A\_{6,0}\boldsymbol{\beta}^6\tag{3}$$

$$G\_1 = A\_{0,1} + A\_{1,1}\beta + A\_{2,1}\beta^2 + A\_{3,1}\beta^3 + A\_{4,1}\beta^4 + A\_{5,1}\beta^5 + A\_{6,1}\beta^6 \tag{4}$$

$$
\beta = \frac{2\phi}{\pi} \tag{5}
$$

$$G\_2 = \frac{\sqrt{2Q}}{\pi} \left( \frac{16}{15} + \frac{1}{3}M\_1 + \frac{16}{105}M\_2 + \frac{1}{12}M\_3 \right) \tag{6}$$

$$G\_3 = \frac{\sqrt{2\mathbb{Q}}}{\pi} \left( \frac{32}{35} + \frac{1}{4}M\_1 + \frac{32}{315}M\_2 + \frac{1}{20}M\_3 \right) \tag{7}$$

$$G\_4 = \frac{\sqrt{2Q}}{\pi} \left( \frac{256}{315} + \frac{1}{5}M\_1 + \frac{256}{3465}M\_2 + \frac{1}{30}M\_3 \right) \tag{8}$$

$$M\_1 = \frac{2\pi}{\sqrt{2Q}}(\Im G\_1 - G\_0) - \frac{24}{5} \tag{9}$$

$$M\_2 = 3\tag{10}$$

$$M\_3 = \frac{6\pi}{\sqrt{2Q}}(G\_0 - G\_1) + \frac{8}{5} \tag{11}$$

where *p* is the internal pressure; *Ri* is the internal radius; *Ro* is the outer radius; *a* is the crack depth; *Q* is a parameter based on crack geometry; *G*0, *G*1, *G*2, *G*3, *G*4, *M*1, *M*2, *M*3, *Ai,j* (*i* ∈ {0,1,2,3,4,5,6}, {*j* ∈ 0,1}), *β* are influence coefficients; *φ* is the included angle; *c* is the half crack length; and *K* is the mode I SIF.

The finite element simulation results are compared with the SIF results calculated according to API 579 criterion. The results are shown in Figure 3. It can be found from the figure that for the pipeline without corrosion defects, the SIF obtained by finite element simulation is very close to the results of theoretical calculation for a large portion of crack depth range, and the maximum error is less than 5%. The accuracy of finite element simulation is proved. Then, the SIF of pipeline with corrosion defects is studied. As is obtained from Figure 3, for the same crack depth, the SIF of the pipeline with corrosion defects is greater than that without corrosion defects. With the increase in crack depth, SIF also increases gradually. The comparison results demonstrate that there is an interacting impact of corrosion and crack defects on SIF values. Therefore, it is necessary to study the interacting impact between corrosion and crack defects.

**Figure 3.** Comparison of pipeline SIF results.

#### **3. The Proposed Crack Propagation Method Based on Extreme Gradient-Boosting Algorithm** *3.1. The Extreme Gradient-Boosting Model*

Extreme Gradient Boosting (XGBoost) is an ensemble machine learning algorithm based on Decision Tree and uses Gradient Boosting as the framework. It is developed from Gradient-Boosting Decision Tree (GBDT). GBDT is an additive model based on boosting, which is a general ensemble method. It employs a forward stagewise algorithm for greedy learning in the training process. In each iteration, GBDT learns a Classification and Regression Tree (CART), where Figure 4 is an example of CARTs, to fit the residual error between the prediction result from previous CARTs and the actual value of the training dataset. In other words, it is there to build a model from the training dataset and create a second model to correct the residual error from the first model. Then, the models are added until the training dataset is predicted relatively accurately, or a maximum number of models is added.

**Figure 4.** An example of CARTs.

Several optimization strategies are added into XGBoost model. Firstly, in order to improve computational accuracy, XGBoost uses the second-order derivative to optimize the objective function. Conversely, GBDT only uses the first-order derivative for optimization. In addition, the objective function of XGBoost utilizes regularization term to simplify the model and avoid overfitting. On the contrary, the GBDT does not have any regularization term in the objective function. XGBoost is able to automatically process default values and compute in parallel through a block storage structure, which cannot be implemented in GBDT. Since XGBoost has a high precision on the second-order derivative and fast parallel computation speed, it is very efficient in data processing and data modeling. In addition, XGBoost is relatively flexible, as it supports classification and regression, and it is able to provide customized objective function. XGBoost can be used with multiple programing languages and platforms. Therefore, XGBoost is widely used in the areas of data mining, recommender system and so on.

The objective function of XGBoost in the training process consists of two parts: loss function and regularization term:

$$Obj(\Theta) = L(\Theta) + \Omega(\Theta) \tag{12}$$

where Θ is the parameters obtained from the training processing; *L*(Θ) is the training error, which denotes the matching degree of the model to the training dataset; Ω(Θ) is the regularization term, which represents the complexity of the model. Assuming that the training dataset is *S* = {(*x*1, *y*1), (*x*2, *y*2), ... ,(*x <sup>n</sup>*, *yn*)}, the training error *L* can be expressed as the following equation:

$$L = \sum\_{i=1}^{n} l(y\_i, \mathfrak{Y}\_i) \tag{13}$$

where *yi* and *<sup>y</sup>*ˆ*<sup>i</sup>* are the target output and the predicted output of the *<sup>i</sup>*-th sample *xi* (*xi* ∈ *Rz*, *<sup>z</sup>* is the number of features of the dataset), respectively, and *n* is the number of samples in the training dataset. For the proposed gradient-boosted machine, *l*(*yi*, *y*ˆ*i*) = (*yi* − *y*ˆ*i*) 2 . The objective is to minimize *Obj*(Θ), which means *L*(Θ) and Ω(Θ) should be relatively small. During the training process, it is required to balance the tradeoffs between bias and variance. Bias is controlled by *L*(Θ) and variance is controlled by Ω(Θ). *L*(Θ) and Ω(Θ) would be relatively large if underfitting. If overfitting, Ω(Θ) would also be relatively large, since the model is weak on scalability and stability. Assuming there are *V* CARTs in the model, then

$$\mathcal{Y}\_{l} = \sum\_{\upsilon=1}^{V} f\_{\upsilon}(\mathbf{x}\_{l})\_{\prime} \; f\_{\upsilon} \in F \tag{14}$$

where *F* is the function space of all the CARTs in the model. *fv*(*xi* ) represents the weight of the *i*-th sample falling on the leaf in the *v*-th tree. For the example in Figure 4, *f*1(sample2) = *w*1-1, *f*2(sample2) = *w*2-2, *f* (sample2) = *w*1-1 + *w*2-2. Then, the model parameters that will be optimized from the training process are Θ = {*f* <sup>1</sup>, *f* <sup>2</sup>, ..., *f <sup>V</sup>*}, where *fv* denotes the weight distribution of the samples falling on the leaf in the *v*-th tree. The objective function is shown in Equation (15):

$$Obj = \sum\_{i=1}^{n} l(y\_i, \hat{y}\_i) + \sum\_{v=1}^{V} \Omega(f\_v) \tag{15}$$

Next, the objective function will be optimized in three steps. The first step is to use the second-order Taylor series expansion to optimize the loss function. The predicted values can also be expressed as

$$\mathfrak{Y}\_{i}^{(u)} = \mathfrak{Y}\_{i}^{(u-1)} + f\_{\mathfrak{u}}(\mathfrak{x}\_{i}) \tag{16}$$

which is the same as the expression of the GBDT. *y*ˆ (*u*) *<sup>i</sup>* is the predicted value of *xi* in tree *u* after the *i*-th iteration. Then, the objective function after the *i*-th iteration can be represented using Equation (17):

$$\Omega \mathcal{O} \mathfrak{j}^{(u)} = \sum\_{i=1}^{n} l(\mathfrak{y}\_{i\prime} \mathfrak{j}\_{i}^{(u-1)} + f\_{\mathfrak{u}}(\mathfrak{x}\_{i})) + \sum\_{v=1}^{\mathfrak{u}} \Omega(f\_{v}) \tag{17}$$

Using the second-order Taylor series expansion, the loss function becomes

$$\sum\_{i=1}^{n} l(y\_i, \boldsymbol{\hat{y}}\_i^{(u-1)} + f\_u(\mathbf{x}\_i)) \approx \sum\_{i=1}^{n} \left[ l(y\_i, \boldsymbol{\hat{y}}\_i^{(u-1)}) + \mathbf{g}\_i f\_u(\mathbf{x}\_i) + \frac{1}{2} h\_i f\_u^{\;2}(\mathbf{x}\_i) \right] \tag{18}$$

where

$$\log\_i = \frac{d(l(y\_{i\prime}\hat{y}\_i^{(\mu-1)}))}{d(\hat{y}^{(\mu-1)})} \tag{19}$$

$$h\_l = \frac{\partial^2(l(y\_{l'} \, \hat{y}\_i^{(u-1)}))}{\partial(\hat{y}^{(u-1)})^2} \tag{20}$$

and the objective function is expressed in Equation (21):

$$\{Ob\}^{(u)} \approx \sum\_{i=1}^{n} \left[ l(y\_i, \mathfrak{f}\_i^{(u-1)}) + g\_i f\_{\mu}(\mathbf{x}\_i) + \frac{1}{2} h\_i f\_{\mu}{}^2(\mathbf{x}\_i) \right] + \sum\_{\upsilon=1}^{\mu} \Omega(f\_{\upsilon}) \tag{21}$$

The second step is to optimize the regularization term by expanding the regularization term and removing the constant term. Since forward calculation is adopted in XGBoost, then the structure of the (*u* − 1)-th tree has been confirmed:

$$l(y\_i, \hat{y}\_i^{(u-1)}) = \text{constant} \tag{22}$$

$$\begin{aligned} \mathop{\rm v}\limits\_{v=1}^{\mu} \Omega(f\_{\upsilon}) &= \Omega(f\_{\mu}) + \mathop{\rm v}\limits\_{\upsilon=1}^{\mu-1} \Omega(f\_{\upsilon}) \\ &= \Omega(f\_{\mu}) + constant \end{aligned} \tag{23}$$

Then the objective function is expressed as follows:

$$\begin{array}{lcl}Obj^{(u)} \approx \sum\_{i=1}^{n} \left[ l(y\_i, \hat{y}\_i^{(u-1)}) + \mathfrak{g}\_i f\_u(\mathbf{x}\_i) + \frac{1}{2} h\_i f\_u^{\ 2}(\mathbf{x}\_i) \right] + \Omega(f\_u) + constant \\ = \sum\_{i=1}^{n} \left[ \mathfrak{g}\_i f\_u(\mathbf{x}\_i) + \frac{1}{2} h\_i f\_u^{\ 2}(\mathbf{x}\_i) \right] + \Omega(f\_u) + \left[ \sum\_{i=1}^{n} l(y\_i, \hat{y}\_i^{(u-1)}) + constant \right] \\ = \sum\_{i=1}^{n} \left[ \mathfrak{g}\_i f\_u(\mathbf{x}\_i) + \frac{1}{2} h\_i f\_u^{\ 2}(\mathbf{x}\_i) \right] + \Omega(f\_u) + constant \end{array} \tag{24}$$

After removing the constant term, the simplified objective function is

$$\left[Obj^{(u)}\right] \approx \sum\_{i=1}^{n} \left[g\_if\_{\mu}(\mathbf{x}\_i) + \frac{1}{2}h\_if\_{\mu}{}^2(\mathbf{x}\_i)\right] + \Omega(f\_{\mu})\tag{25}$$

The last step of the optimization process is to merge the coefficients of the first-degree term and the quadratic term. Regarding the definition of a tree, the weight vector of leaves is set as *<sup>w</sup>* ∈ *<sup>R</sup>T*and the mapping relationship between the leaves (viz. the structure of the tree) is defined as *<sup>q</sup>* : *<sup>R</sup><sup>Z</sup>* <sup>→</sup> {1, 2, 3, . . . , *<sup>T</sup>*} where *<sup>T</sup>* is the number of leaves in the tree. Then, *q*(*x*) denotes the location of the leaf, for sample *x*. For the example in Figure 4, *q*(sample2) = 1 in tree 1, *q*(sample2) = 2 in tree 2. *ft*(*x*) can be represented by

$$f\_{\boldsymbol{\mu}}(\boldsymbol{x}) = w\_{\boldsymbol{\eta}(\boldsymbol{x})} \tag{26}$$

Here, the number of leaves *T* and smoothness of leaf weight (viz. L2 norm of leaf weights) are used to describe the complexity of the tree, so

$$
\Omega(f\_{\mu}) = \gamma T + \frac{1}{2}\lambda \sum\_{j=1}^{T} w\_j^2 \tag{27}
$$

For the example in Figure 4, Ω(*f*1) = *γ*3 + <sup>1</sup> 2*λ* ' *w*1−<sup>1</sup> <sup>2</sup> <sup>+</sup> *<sup>w</sup>*1−<sup>2</sup> <sup>2</sup> <sup>+</sup> *<sup>w</sup>*1−<sup>3</sup> 2 ( , Ω(*f*2) = *γ*2 + <sup>1</sup> 2*λ* ' *w*2−<sup>1</sup> <sup>2</sup> <sup>+</sup> *<sup>w</sup>*2−<sup>2</sup> 2 ( . *Ij* = {*i* | *q*(*xi*) = *j*} is the instance set in leaf *j*, *j* = 1, 2, ... , *T*. Grouping all the training samples based on leaves and utilizing Equations (26) and (27), then, the objective function is

$$\left[Obj^{(u)} \approx \sum\_{j=1}^{T} \left[ (\sum\_{i \in I\_{j}} g\_{i}) w\_{j} + \frac{1}{2} (\sum\_{i \in I\_{j}} h\_{i} + \lambda) w\_{j}^{2} \right] + \gamma T \tag{28}$$

At last, merging the first-degree term and the quadratic term, then

$$\langle Obj^{(u)} \approx \sum\_{j=1}^{T} \left[ G\_j w\_j + \frac{1}{2} (H\_j + \lambda) w\_j^2 \right] + \gamma T \tag{29}$$

where

$$G\_j = \sum\_{i \in I\_j} \mathcal{g}\_{i\prime} H\_j = \sum\_{i \in I\_j} h\_i \tag{30}$$

For each leaf *j*, the objective function is expressed as follows:

$$f(w\_{\dot{l}}) = G\_{\dot{l}}w\_{\dot{l}} + \frac{1}{2}(H\_{\dot{l}} + \lambda)w\_{\dot{j}}^2\tag{31}$$

As the objective function of each leaf in the overall objective function is independent, then the overall objective function will achieve the minimum value when each leaf's objective function is minimized. The optimal solution of the quadratic function of one variable is

$$w\_j^\* = -\frac{G\_j}{H\_{\bar{l}} + \lambda} \tag{32}$$

At this point, each leaf weight is optimized, and the overall objective function achieves its optimal value, viz. the minimum value:

$$Obj^\* = -\frac{1}{2} \sum\_{j=1}^{T} \frac{G\_j}{H\_j + \lambda} + \gamma T \tag{33}$$

The structure of the tree is also the best at this time. The optimal objective functions of Figure 4 are shown in Figure 5. The fewer objective functions there are, the better the tree structures are.

**Figure 5.** Objective function of the example in Figure 4.

In the actual training process, finding the optimal split point is a key problem. The applicable methods include greedy algorithm, approximate algorithm, weighted quantile sketch and sparsity-aware split finding. The greedy algorithm is the most commonly used.

#### *3.2. The Proposed Model Based on XGBoost*

In the proposed method, the scikit-learn wrapper interface for XGBoost was utilized to construct models to predict the SIF interaction impact ratio at the deepest point of the crack defect for pipelines with corrosion and crack defects. Based on the observations from finite element modeling, the size of crack and corrosion defects, and the axial distance between them, can affect SIF results. Therefore, the input variables of the proposed model are the length and depth of the crack defect, the length and depth of the corrosion defect, and the axial distance between the crack and the corrosion defects. The output variable is the interaction impact ratio *α*. The input and output variables are shown in Table 1. The crack length is assumed in the range of 15.2 mm–76.0 mm, and the crack depth is in the range of 2 mm–12 mm. The axial distance between the corrosion and crack defects is from 150 mm to 500 mm. The depth of the corrosion defect is from 2 mm to 14 mm.

**Table 1.** Input and output variables of the proposed XGBoost models.


In this paper, two approaches are provided to fit *α*. The first one is to directly construct a XGBoost model to predict *α*. The second one is to construct two XGBoost models to fit SIF values with and without considering the interaction impact, which are *K*\* and *K*, respectively. Then, the interaction impact ratio can be calculated with the formula *α* = *K*\*/*K*. It is worth noting that only crack length and crack depth have an impact on SIF without considering interaction impact. Then, in the process of fitting *K*, there are only two input variables, viz. crack depth and crack length.

The samples used for modeling are synthetic data from finite element modeling. In total, 385 pieces of data are generated. A share of 80% of these data was randomly selected as the training set. The remaining data are the testing set. The scikit-learn API for XGBoost regression has a lot of parameters to set. In this paper, five parameters are selected for parameter tuning to get the best model structure and parameters: the number of gradientboosted trees, the maximum depth of a tree, the minimum sum of instance weight needed in a child, *L*1 and *L*2 regularization terms on weights. The adjusting ranges for these five parameters are shown in Table 2. Increasing the maximum depth of a tree will make the model more complex, and it will be more likely to overfit, so the maximum value for this parameter is set to 10 in this paper. If the sum of instance weight in a leaf node is less than the minimum sum of instance weight needed in a child, the building process will stop further partitioning. Regarding the *L*1 and *L*2 regularization terms on weights, increasing their values will make the model more conservative. The learning rate is set at 0.1, which updates the weights to prevent overfitting and makes the boosting process more conservative. For the other parameters, such as the initial prediction score of all instances (global bias), minimum loss function required to make a further partition on a leaf node of the tree, etc., the default values in the scikit-learn API are applied.

**Table 2.** Parameter tuning for XGBoost models.


In the training process, a grid search method with 5-fold cross-validation was applied to select the best combination of the tuning parameters based on the determination coefficient *R*2, which describes the goodness of fit of the current trained model. In other words, the original training set was re-segmented into the training set and validation set with the ratio of 4:1 five times, as shown in Figure 6. For each combination of the tuning parameters, the training set was used to train the model, and the validation set was used to evaluate the model's performance five times and compute the average performance, viz. average *R*2, with these five times' results. This method can reduce training bias and improve the model's stability. After all the combinations' results are obtained, the model with the highest *R*<sup>2</sup> has the best combination of the tuning parameters. The values of *R*<sup>2</sup> are between 0 and 1. A value much closer to 1 indicates the regression model has a higher fitting degree.

In the actual training process, a pipeline of transforms with a final estimator (viz. model to be fitted) is utilized. This method is to sequentially apply a list of transforms and a final model. Intermediate steps of the pipeline must implement fit and transform methods, while the final model only needs to implement the fit method. In this paper, a pipeline consisting of a standard scaler and an XGBoost model is applied. The standard scaler is to normalize data to make its features have zero mean and unit variance. The standard scaler fits to the training set and transforms the training set and validation set.

The overall process of the first approach for constructing the XGBoost model is as follows:

Step 1. Randomly split the samples (output variable is SIF interaction impact ratio *α*) into training set and testing set with the ratio 8:2.

Step 2. Employ a pipeline consisting of a standard scaler and an XBoost model to the original training set for training. In detail, the 5-fold cross-validated grid search method is applied to the original training set to select the best model structure and parameters among all the combinations of the tuning parameters. The model with the highest *R*<sup>2</sup> is the best

model. The best model is then saved and can be directly applied to new data to acquire prediction values.

Step 3. Feed the testing set to the trained model to obtain the value of *R*2, which indicates the ability fitting to new data with the trained model. The closer that *R*<sup>2</sup> is to 1, the better structured the model is. If the value is close to 1, then the trained model can be used to directly predict interaction impact ratio *α*.

**Figure 6.** The procedure of the proposed algorithm.

Similarly, the overall process of the second approach for constructing the two XGBoost models is:

Step 1. Randomly split the samples (output variables are SIF values with and without considering interaction impact, viz. *K\** and *K*) into training set and testing set with the ratio 8:2.

Step 2. When the output variable is *K*, the input's variables are crack depth and crack length, and employ a pipeline of a standard scaler and an XGBoost model to the original training set for training. In the same way, the 5-fold cross-validated grid search method is applied to the original training set to select the best model. The best model is saved to predict *K*.

Step 3. When the output variable is *K\**, the input includes all the five input variables and employs a pipeline of a standard scaler and an XGBoost model to the original training set for training. In the same way, the 5-fold cross-validated grid search method is applied to the original training set to select the best model. The best model is then saved to predict *K\**.

Step 4. Respectively, feed the two testing sets to the two trained models to obtain the values of *R*2. If the two values of *R*2are close to 1, then the two trained models can be used to predict SIF with and without considering interaction impact, respectively.

Step 5. Calculate predicted interaction impact ratio *α* on the testing set with predicted *K* and *K\** and compare the predicted values with the target ratio values by calculating *R*2.

#### *3.3. The Pipeline Corrosion and Fatigue Crack Growth Models*

In the proposed model, corrosion defect is assumed to grow linearly. The growth of the corrosion depth is characterized by

$$d(t) = d\_0 + \mathcal{g}\_d t \tag{34}$$

where *d*<sup>0</sup> represents the corrosion initial depth, *gd* is the growth rate of corrosion depth, and *t* is the propagation time. The corrosion depth is used as the input variable in the XGBoost model to calculate the SIF interaction impact ratio. In this paper, the corrosion depth growth rate is assumed to be 0.3 mm/year [11].

Pipeline fatigue crack growth is predicted using the physics-based methods governed by Paris' law, which was employed in [28–30]. Based on Paris' law and the proposed model for evaluating the SIF interaction impact between corrosion and crack defects, the fatigue crack growth model is introduced in the following equation:

$$da/dN = \mathbb{C}(\Delta \mathcal{K} \alpha)^m \tag{35}$$

where *da*/*dN* is crack growth rate; *a* is crack depth; *N* is the number of loading cycles; α is the SIF interaction impact ratio; and Δ*K* is the range of SIF. *C* and *m* are material-related model parameters, which can be estimated via experiments. In this paper, it is assumed that model parameters *<sup>C</sup>* = 5 × <sup>10</sup>−12, *<sup>m</sup>* =3[21]. Methods based on FE and XGBoost models are employed to calculate the SIF and SIF interaction impact ratio at the deepest point of the fatigue crack. In this study, this paper focuses on the crack depth growth, since the length is mostly unchanged.

#### **4. Results**

When directly fitting the SIF interaction impact ratio, the average determination coefficient on the validation sets during cross validation is 0.9935, and the standard deviation is 0.0041. Thus, it can be seen that the trained model has a relatively high stability. The prediction result on the testing set is as Figure 7 shows. On the testing set, the determination coefficient *R*<sup>2</sup> is 0.9876, which means the developed model can accurately predict the SIF interaction impact ratio. At this point, the number of gradient-boosted trees is 110, the maximum depth of a tree is 6, the minimum sum of instance weight needed in a child is 1, and the *L*1 and *L*2 regularization terms on weights are 0.05 and 0.1, respectively.

The prediction results of the SIF interaction impact ratio are shown in Figure 8. As observed in Figure 8, it can be found that the interaction impact ratio decreases as the crack depth *a* increases. From Figure 8a to Figure 8c, as corrosion depth increases from 4 mm to 10 mm, the SIF interaction impact ratio increases a lot when the axial distance between two corrosion and crack defects remains the same. Thus, the corrosion depth does affect the SIF interaction impact ratio a lot. The highest SIF interaction impact ratio in Figure 8c is 1.1848, which means it is necessary to consider the interaction impact between these two defects in the crack propagation process. The comparison results for different corrosion depths when crack depth is equal to 6 mm are shown in Figure 8d. From all these four figures, it can be found that the SIF interaction impact ratio is overall decreasing as the axial distance increases. These ratios first decrease quickly when the axial distance is smaller than around 175 mm and then decrease relatively slowly when the axial distance is in the range of 175 mm and 240 mm. When the axial distance is bigger than 240 mm, the decreasing speed is getting even smaller. This is because the corrosion defect moves away from the stress concentration zone of the crack defect.

**Figure 7.** Prediction results of SIF interaction impact ratio based on approach 1.

**Figure 8.** Prediction results of SIF interaction impact ratio based on approach 1. Corrosion depth *d*: (**a**) 4 mm; (**b**) 6 mm; (**c**) 10 mm; (**d**) Comparison results for different corrosion depths when crack depth *a* is 6 mm.

For the second approach, the average determination coefficient on the validation sets is 1.0000 and the standard deviation is 0.0000 when predicting SIF values without considering interaction impact, which means the trained model is relatively stable. Here, the number of gradient-boosted trees is 110, the maximum depth of a tree is 3, the minimum sum of instance weight needed in a child is 1, and the *L*1 and *L*2 regularization terms on weights are both 0.05. The prediction result on the testing set is shown in Figure 9, where the

determination coefficient *R*<sup>2</sup> is 1.0000. When considering interaction impact, the average determination coefficient on the validation sets is still 0.9998, and the standard deviation is still 0.0001. However, the selected structure and parameters of the model are different. The number of gradient-boosted trees is 110, the maximum depth of a tree is 4, the minimum sum of instance needed in a child is 2, and the *L*1 and *L*2 regularization terms on weights are 0.1 and 1, respectively. On the testing set, the determination coefficient *R*<sup>2</sup> is 0.9992, and the prediction result is as displayed in Figure 10. It can be seen that for these two predictive models, the performance is quite stable on the validation sets and very accurate on the testing set. Therefore, it can be concluded that these two XGBoost models can predict SIF values with and without considering interaction impact efficiently and accurately. Furthermore, this indicates these two models are able to predict interaction impact ratio efficiently and accurately, since the ratio is calculated from the predicted results of these two models. After the predictive results are obtained from these two models, the SIF interaction impact ratio *α* can be calculated with the equation *K\**/*K*. For the testing set in this paper, the result is shown in Figure 11. At this time, the determination coefficient *R*<sup>2</sup> of the SIF interaction impact ratio on the testing set is 0.9852. From the experimental result, it can also be concluded that the two trained XGBoost models can predict SIF interaction impact ratio at the deepest point of the crack defect, considering the interaction between corrosion and crack defects accurately and efficiently.

**Figure 9.** Prediction results of SIF values without considering interaction impact based on approach 2.

The comparison results of SIF values with and without considering interaction impact are shown in Figure 12. When increasing the crack depth *a* from 2 mm to 9 mm, *K*\* and *K* values both gradually increase as expected. From Figure 12a to Figure 12c, *K*\* values gradually increase, while *K* values remain the same as the corrosion depths increases. From the observations of these in Figure 12d, *K*\* and *K* values have relatively big differences when the axial distance between two defects is smaller than around 240 mm.

**Figure 10.** Prediction results of SIF values considering interaction impact based on approach 2.

**Figure 11.** Prediction results of SIF interaction impact ratio based on approach 2.

The crack propagation models are built based on the proposed XGBoost model and Paris' law. The crack initial depth is set at 2 mm, since SIF interaction impact ratio is relatively large when crack depth is small. The corrosion initial depth is assumed at 6 mm. Figure 13a–d show the comparison results of crack depth growth models for different axial distances between two defects using approach 1 and 2, respectively. The red dash lines represent the crack critical depth, which is approximately 80% of the wall thickness. When the crack depth exceeds the crack critical depth, it is considered a failure. The comparison results shown in Figure 13 indicate that the crack depth predicted by approach 2 reaches the threshold more quickly than approach 1. The comparison results of the crack depth propagation curves for different axial distances based on approach 2 are shown in Figure 14. If the interaction impact between two defects is not considered, it takes about 14.8 years to fail. Meanwhile, when considering the interaction impact between two defects by

implementing the proposed model, the failure time changes from 9.2 to 13.8, 14.0 and 14.5 years, as the axial distance changes from 150 to 200, 300 and 350 mm. There is a big difference between the two crack propagation curves in Figure 13a, since the corrosion defect is in the stress concentration zone (axial distance smaller than 175 mm).

**Figure 12.** The comparison of SIF results based on approach 2. Corrosion depth *d*: (**a**) 4 mm; (**b**) 6 mm; (**c**) 10 mm; (**d**) Comparison results for different corrosion depths when crack depth *a* is 6 mm.

To perform sensitivity analysis regarding corrosion initial depth, the aggressive case was studied, in which the axial distance was set at 150 mm. The corrosion initial depth varies from 2 mm to 8 mm. Figure 15a–d show the results of crack depth growth models for different corrosion initial depths. These figures indicate that the crack depth grows more quickly using approach 2 than approach 1. The comparison results of the crack depth propagation curves for different corrosion initial depths based on approach 2 are shown in Figure 16. If not considering the interaction impact between two defects, it also takes about 14.8 years to fail. From the comparison results in Figures 15 and 16, the time to reach critical crack depth is 10.8, 9.9, 9.2 and 8.8 years, respectively, as the corrosion initial depth changes from 2 to 3, 6 and 8 mm. From the experimental results obtained from Figures 13–16, it can be concluded that the interaction impact between corrosion and crack defects affects the propagation of fatigue crack a lot. Thus, it is necessary to consider the SIF interaction impact ratio in the remaining useful life prediction, especially when the corrosion defect is in the stress concentration zone.

**Figure 13.** Investigations of the interaction impact on crack depth growth for different axial distances. Axial distance: (**a**) 150 mm; (**b**) 200 mm; (**c**) 300 mm; (**d**) 500 mm.

**Figure 14.** Crack depth growth curves for different distances based on approach 2.

**Figure 15.** Investigations of the interaction impact on crack depth growth. Corrosion initial depth: (**a**) 2 mm; (**b**) 4 mm; (**c**) 6 mm; (**d**) 8 mm.

**Figure 16.** Crack depth growth curves for different corrosion initial depths based on approach 2.

#### **5. Conclusions**

The existing reported work only focuses on pipeline life prediction with single or multiple defects of the same type. The interaction impacts between different types of defects are not considered. In this work, the interaction impacts between crack and corrosion defects were studied, and a fatigue crack propagation method considering these impacts was proposed based on XGBoost models and Paris' law. Crack size, corrosion size, and the axial distance between these two defects were all considered in the proposed method. In addition, this paper introduced SIF interaction impact ratio to describe how the corrosion defect affects the stress concentration zone of the fatigue crack. Two approaches were implemented for SIF interaction impact ratio prediction. The first one directly fitted and predicted SIF interaction impact ratio with the synthetic samples from finite element modeling. The second one fitted and predicted the SIF with and without considering interaction impacts, respectively, and then calculated the SIF interaction impact ratio. Examples were used to demonstrate the proposed method. The determination coefficients of these two approaches on the testing sets were 0.9876 and 0.9852, respectively, which was quite close to 1. Therefore, it can be concluded that the developed method can predict fatigue crack growth accurately. Several key findings are listed below:

The SIF interaction impact ratio decreases as the crack depth increases. It increases as the corrosion depth increases.

The SIF interaction impact ratio is gradually decreasing as the axial distance increases. This ratio is relatively large when the axial distance is smaller than 240 mm.

The time to reach critical crack depth decreases as the corrosion initial depth increases or the axial distance decreases.

The method developed in this paper can support the decision making in pipeline integrity planning, especially when the corrosion defect is relatively close to the crack defect. However, the proposed method only considered the interacting impact between two defects. More efficient crack and corrosion propagation models considering more than two defects are desired in future research. Another research topic is to develop crack propagation models for different types of crack shapes instead of semi-elliptical shapes.

**Author Contributions:** Conceptualization, Y.W., W.X. and J.Z.; methodology, M.X.; software, Y.W. and M.X.; validation, Y.W., J.Z. and M.X.; formal analysis, M.X., W.X. and X.P.; investigation, M.X.; resources, M.X.; data curation, M.X.; writing—original draft preparation, Y.W. and M.X.; writing review and editing, M.X. and W.X.; visualization, M.X. and X.P.; supervision, M.X. and X.P.; project administration, M.X. and X.P.; funding acquisition, M.X. and X.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China, grant number 72001039, 71671035, 12102090.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**



#### **References**


## *Article* **An Ensemble Prognostic Method of Francis Turbine Units Using Low-Quality Data under Variable Operating Conditions**

**Ran Duan, Jie Liu, Jianzhong Zhou \*, Pei Wang and Wei Liu**

School of Civil and Hydraulic Engineering, Huazhong University of Science and Technology, Wuhan 430074, China; Duan\_Ran@hust.edu.cn (R.D.); jie\_liu@hust.edu.cn (J.L.); w\_sara@hust.edu.cn (P.W.); l\_wei@hust.edu.cn (W.L.)

**\*** Correspondence: JianzhongZhou@hust.edu.cn

**Abstract:** The prognostic is the key to the state-based maintenance of Francis turbine units (FTUs), which consists of performance state evaluation and degradation trend prediction. In practical engineering environments, there are three significant difficulties: low data quality, complex variable operation conditions, and prediction model parameter optimization. In order to effectively solve the above three problems, an ensemble prognostic method of FTUs using low-quality data under variable operation conditions is proposed in this study. Firstly, to consider the operation condition parameters, the running data set of the FTU is constructed by the water head, active power, and vibration amplitude of the top cover. Then, to improve the robustness of the proposed model against anomaly data, the density-based spatial clustering of applications with noise (DBSCAN) is introduced to clean outliers and singularities in the raw running data set. Next, considering the randomness of the monitoring data, the healthy state model based on the Gaussian mixture model is constructed, and the negative log-likelihood probability is calculated as the performance degradation indicator (PDI). Furthermore, to predict the trend of PDIs with confidence interval and automatically optimize the prediction model on both accuracy and certainty, the multiobjective prediction model is proposed based on the non-dominated sorting genetic algorithm and Gaussian process regression. Finally, monitoring data from an actual large FTU was used for effectiveness verification. The stability and smoothness of the PDI curve are improved by 3.2 times and 1.9 times, respectively, by DBSCAN compared with 3-sigma. The root-mean-squared error, the prediction interval normalized average, the prediction interval coverage probability, the mean absolute percentage error, and the R2 score of the proposed method achieved 0.223, 0.289, 1.000, 0.641%, and 0.974, respectively. The comparison experiments demonstrate that the proposed method is more robust to low-quality data and has better accuracy, certainty, and reliability for the prognostic of the FTU under complex operating conditions.

**Keywords:** Francis turbine unit; prognostic; performance state evaluation; degradation trend prediction; DBSCAN; Gaussian mixture model; NSGA-II; Gaussian process regression

#### **1. Introduction**

With the optimization of global energy structure, hydropower has been vigorously developed [1]. As critical equipment of hydropower utilization, it is important to ensure the safe and stable operation of Francis turbine units (FTUs) [2]. In order to transform the service mode of the FTUs from traditional schedule-based maintenance to state-based maintenance, the prognostic of FTUs has received more and more attention [3]. There are two main parts in the prognostic of FTUs, including performance state evaluation and degradation trend prediction [4]. Based on this, the current and future states of FTUs can be determined so as to formulate a targeted maintenance strategy. However, in practical engineering environments, there are three significant difficulties, including low data quality, complex variable operation conditions, and prediction model parameter optimization.

**Citation:** Duan, R.; Liu, J.; Zhou, J.; Wang, P.; Liu, W. An Ensemble Prognostic Method of Francis Turbine Units Using Low-Quality Data under Variable Operating Conditions. *Sensors* **2022**, *22*, 525. https:// doi.org/10.3390/s22020525

Academic Editors: Yongbo Li, Bing Li, Jinchen Ji and Hamed Kalhori

Received: 12 December 2021 Accepted: 5 January 2022 Published: 11 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The state monitoring systems of FTUs work in the harsh environment of high humidity, high vibration, and strong electromagnetic. Due to sensor failure, electromagnetic interference, and missing communication packets, the data quality of the actual measured signal is usually low, which is mainly manifested as data anomaly and data loss [5,6]. Most signal denoising methods, including spectrum analysis and signal reconstruction, are effective for signals with a high sampling rate and consistent sampling frequency [7–9]. However, data loss has a great influence on their practical application. Suppression of abnormal data in practical low-quality data sets of FTUs is rarely discussed in existing studies. As an unsupervised learning algorithm, clustering can adaptively detect potential patterns among multi-dimensional data [10]. The density-based spatial clustering of applications with noise (DBSCAN) has the ability to identify the isolated noise from data sets with arbitrary shapes, and the DBSCAN is confirmed to be efficient and robust for low-quality data sets [11,12]. Therefore the DBSCAN is adopted to clean the actual monitoring data set of the FTU in this study.

Generally, FTUs need to participate in the load and frequency regulation of the power grid, and the incoming water levels fluctuate significantly with seasonal variation. Therefore, the operating condition parameters of FTU change frequently [13]. Since the monitoring vibration signals are highly correlated to operation condition parameters, the traditional performance evaluation method with a fixed threshold is difficult to reflect the actual state of FTUs accurately [14]. Machine learning has been widely used in equipment fault diagnosis and performance evaluation due to its good pattern recognition capability [15]. To establish the performance state evaluation model under variable operation conditions, An et al. and Shan et al. adopted the radial basis function and backpropagation neural network (BPNN), respectively, to fit the mapping relationship between water head, active power, and vibration amplitude, as the healthy state model (HSM) of the FTU [4,16]. These studies established the definite functional relation between operation condition parameters and monitoring signal, but the randomness of the signal may still affect the accuracy and stability of these value-to-value mapping models. To deal with it, Rai et al. adopted the Gaussian mixture model (GMM) to fit the probability density distribution (PDD) of vibration features as the HSM of a rolling bearing. The results showed the HSM based on GMM was more accurate and monotonous [17]. However, in this work, the operation condition parameters were constant. Thus, the variable operation condition parameters and the signal randomness are sufficiently considered in this paper, and the GMM is adopted to fit the probability distribution of the running data set as the HSM of the FTU under complex operating conditions.

After constructing the proposed HSM, the real-time state of the FTU can be quantified as performance degradation indicator (PDI) values. Thus, time-series prediction methods can be used to solve the degradation trend prediction problems. For example, Li et al. combined the convolutional filters and the gated recurrent unit to construct a degradation trend prediction model for a turbofan based on data [18]. Jin et al. presented a novel adaptive residual long short-term memory network to predict cutter head torque across domains. The prediction performance is effectively improved by using the knowledge of the source domain dataset [19]. However, most deep learning models only output point prediction values, and it is difficult to quantify the uncertainty of results directly [20]. As a probability prediction model, the Gaussian process regression (GPR) has been widely proven to perform well in uncertain prediction problems [21–23]. However, its property is greatly affected by model parameters. Manual parameter tuning is inefficient and relies on prior knowledge. Therefore, some research adopted intelligent optimization algorithms to optimize prediction model parameters automatically [24–26]. However, the accuracy of the prediction model was taken as the only optimization objective in these researches. The confidence interval (CI) width represents the certainty of probability prediction, which is also an important objective to consider. As one of the most popular multiobjective optimization algorithms, the non-dominated sorting genetic algorithm (NSGA-II) reduces the complexity of genetic algorithms with fast calculation and good convergence [27–29]. Thus the NSGA-II is adopted to optimize the GPR on the two objectives of accuracy and certainty to construct the multiobjective GPR (MOGPR) for degradation trend prediction of FTUs.

According to the above discussion, in the relative field of the prognostic of FTUs, there are few targeted processing methods for low-quality data obtained in practical engineering environments. Meanwhile, existing evaluation methods of FTU under complex conditions rarely consider the randomness of monitoring data. Moreover, in the interval prediction of PDIs, it is still difficult to optimize model parameters automatically while considering both the accuracy and certainty of results.

In this paper, an ensemble prognostic method of FTUs using low-quality data under variable operation conditions is proposed. The major contributions are outlined as follows:


The remainder of this paper is organized as follows: In Section 2, the proposed approach framework and related theories are expounded. In Section 3, an engineering application of the proposed method is presented. In Section 4, different data cleaning methods, HSM construction, and trend prediction are compared and discussed. Finally, a summary is presented in Section 5.

#### **2. Proposed Approach**

To evaluate and predict the performance state of FTUs more accurately under the circumstance of low-quality data and variable operation conditions, an ensemble prognostic method of FTUs using low-quality data under variable operation conditions is proposed in this paper, which mainly includes four steps: data acquisition, data cleaning, performance state evaluation, and degradation trend prediction. The proposed framework flow chart is shown in Figure 1. First, the monitoring data of the water head, active power, and vibration amplitude of the top cover are integrated to construct the running data set of the FTU. Second, aiming at data loss and data anomaly, the data cleaning operation based on DBSCAN is implemented. Third, to solve the problem of monitoring data fluctuating with operation conditions, the HSM is established based on GMM. Then the negative log-likelihood probability (NLLP) is calculated as the PDI of the FTU. The relative trend of PDIs over time reflects the process of performance degradation. Fourth, the MOGPR model is constructed to predict the performance degradation trend of the FTU and takes both prediction accuracy and confidence interval into consideration. Finally, the validity of the performance evaluation model and degradation trend prediction model is evaluated by multi-criterions.

**Figure 1.** Flowchart of the proposed approach.

#### *2.1. Data Acquisition*

Due to their task of regulating the power grid, the operation condition parameters of FTUs change more frequently than other kinds of rotating machines. Therefore, unlike most traditional methods, operation condition parameters are taken into full consideration in the performance state evaluation in this study. The major operation condition parameters of FTUs include rotation speed (*n*), water head (*H*), active power (*P*), flow rate (*Q*), guide vanes opening degree (*α*), etc. Since the duration of transient operation conditions is particularly short compared with the total working time of FTUs, only steady operation conditions are considered in this research. So n can be considered equal to the rated value of the FTU. For a specific FTU, there is a certain relationship between *H*, *P*, *Q*, and *α*. If two of them are identified, others can be inquired from the comprehensive operation curve of the FTU [30]. Thus H and P are chosen as the studied operation condition parameters. As a critical component of the FTU, the top cover is used to seal the runner and support the main shaft. Its vibration amplitude (*V*) can reflect the performance state of the FTU. In conclusion, the running sample set of the FTU is formed by (*H*, *P*, *V*), including both operation condition parameters and monitoring data.

#### *2.2. Data Cleaning Based on DBSCAN*

To solve the data anomaly in the raw running sample set, a data cleaning approach based on DBSCAN is proposed. Traditional statistics-based methods such as the 3-sigma principle are effective in detecting singular points, but they can hardly identify outliers whose value is within the normal range [31]. DBSCAN clustering algorithm has the advantages of high adaptability, extensible dimension, and applicability to arbitrary data shape [32]. It can automatically identify the singularities and outliers in the multidimensional data set. Its schematic diagram is shown in Figure 2, and the main steps are as follows [33]:

**Figure 2.** Schematic diagram of the DBSCAN.

**Step 1**: For the sample set Φ = {*ϕ*1, *ϕ*2, ··· , *ϕN*}, the region whose Euclidean distance from the sample point *ϕ<sup>i</sup>* is less than *ε* is defined as the *ε* neighborhood of *ϕi*. If the sample number in the neighborhood of *ϕ<sup>i</sup>* is greater than the threshold M, *ϕ<sup>i</sup>* is defined as the core point. Samples that are not core points but are in the *ε* neighborhood of a core point are defined as boundary points. Samples that are not core points and are not boundary points are defined as noise points. The *ε* neighborhood of a core point is defined as a temporary cluster *C j*.

**Step 2**: Traversal the sample set Φ, if the sample point *ϕ<sup>i</sup>* in the temporary cluster *C m* is also a core point in another temporary cluster *C n*, the union set *C <sup>m</sup>* ∪ *C n* is defined as a new temporary cluster.

**Step 3**: Repeat Step 2 until the sample points in each temporary cluster are all core points or boundary points, then each temporary cluster is determined as a cluster *C k*.

#### *2.3. Performance State Evaluation Based on GMM*

The performance degradation process of the FTU is reflected in the PDD variation of the running sample set (*H*, *P*, *V*). Thus the GMM is adopted to fit the three-dimensional PDD function of the healthy data as the HSM. In GMM, the population distribution of the sample set is assumed to be a combination of a series of Gaussian distributions:

$$P(X|\theta) = \sum\_{k=1}^{K} w\_k N\_k(\mathbf{x}|\theta\_k) \tag{1}$$

where *X* represents the healthy data set, *wk* represents the weight coefficient. *K* ∑ *k*=1 *wk* = 1. *<sup>θ</sup><sup>k</sup>* <sup>=</sup> (*μk*, *<sup>σ</sup>k*) is the distribution parameters. *Nk*(*x*|*θk*) is the PDD function of the *<sup>k</sup>*th Gaussian component given by:

$$N\_k(\boldsymbol{x}|\theta\_k) = \frac{1}{\sqrt{2\pi}\sigma\_k} \exp\left(-\frac{(\boldsymbol{x}-\boldsymbol{\mu}\_k)^2}{2\sigma\_k^2}\right) \tag{2}$$

The expectation-maximization algorithm is used to estimate *wk*, and *θ<sup>k</sup>* of the GMM [34]:

$$\hat{\theta} = \operatorname{argmax} P(\mathbf{x}|\theta) = \operatorname{argmax} \sum\_{k=1}^{K} w\_k N\_k(\mathbf{x}|\theta\_k) \tag{3}$$

where ˆ *θ* is the maximum likelihood estimate value of *θ*.

The NLLP represents the probability that current data is observed based on the prior given by the GMM. It represents the difference between the PDD of running data *X* and the constructed HSM. So the NLLP is calculated as the PDI of the FTU.

$$NLPP = -\log P(X'|\theta) \tag{4}$$

#### *2.4. Degradation Trend Prediction Based on MOGPR*

Through the above procedures, the performance state of the FTU is quantified as PDIs. The performance degradation prediction of the FTU is transformed to a time series prediction task. In this section, the GPR and NSGA-II are combined to construct the MOGPR model for the degradation trend prediction.

#### 2.4.1. GPR Algorithm

As a nonparametric Bayesian inference model, GPR is widely used in probability interval prediction [35,36]. In GPR, the distribution of possible values at each time point is assumed to obey the Gaussian distribution, which can be expressed in terms of the expectation function *μ<sup>f</sup>* and the covariance function *κ*(·) as:

$$Y = f(X) \sim N\left(\mu\_f, \kappa\_{ff}\right) \tag{5}$$

where *X* is the independent variables, *κ*(·) is also called the kernel function, *κ f f* = *κ*(*X*, *X*). The joint distribution of the actual observed values *Y*∗ and *Y* also follows the Gaussian distribution.

$$N\left[\begin{array}{cc} Y \\ Y^\* \end{array}\right] \sim N\left(\left[\begin{array}{cc} \mu\_f \\ \mu\_y \end{array}\right] , \left[\begin{array}{cc} \kappa\_{ff} & \kappa\_{fy} \\ \kappa\_{fy}^T & \kappa\_{yy} \end{array}\right] \right) \tag{6}$$

where *κ f y* = *κ*(*X*, *X*∗), *κyy* = *κ*(*X*∗, *X*∗).

Finally, the GPR model can be expressed as:

$$Y \sim N\left(\kappa\_{fy}^T \kappa\_{ff}^{-1} Y^\* + \mu\_{f'} \kappa\_{yy} - \kappa\_{fy}^T \kappa\_{ff}^{-1} \kappa\_{fy}\right) \tag{7}$$

The kernel functions *κ*(·) are used to enhance the representation of relationships between input samples. Various kernel functions have different properties and characteristics. The commonly used kernel functions include radial basis function (RBF), matern (MA), rational quadratic (RQ). The kernel function adopted in this study is composed of them.

$$\kappa\_{RBF} \left( X\_{i\prime} X\_{j} \right) = \exp \left( -\frac{||X\_{i} - X\_{j}||^{2}}{2l^{2}} \right) \tag{8}$$

$$\begin{cases} \kappa\_{MA} \left( X\_i, X\_j \right) = \frac{1}{\sqrt{2\pi}} \rho^{1.5} \left( 1 + \frac{\sqrt{3}\rho}{l} \right) \exp\left( -\frac{\sqrt{3}\rho}{l} \right) \\\ \rho = \frac{\sqrt{3}}{l} ||X\_i - X\_j|| \end{cases} \tag{9}$$

$$\kappa\_{RQ}(X\_{i\prime}X\_j) = 1 + \frac{\left\|X\_i - X\_j\right\|^2}{2l^2} \tag{10}$$

These kernel functions all contain a length scale *l*. To improve the performance of the model, the common approach is to adjust *l* manually, which is less efficient. Aiming at this problem, the NSGA-II algorithm is introduced to automatically optimize the GPR model. Since GPR outputs the PDD of predicted values, quantiles and confidence intervals (CIs) can be calculated directly. Interval prediction results are usually evaluated by multiple

criteria, including the root-mean-square error (RMSE), which reflects the accuracy, and the prediction interval normalized average (PINAW), which reflects certainty, defined as:

$$RMSE = \sqrt{\frac{1}{N} \sum\_{i=1}^{N} (\hat{y}\_i - y\_i^\*)^2} \tag{11}$$

$$PINAW = \frac{1}{N\Delta y^\*} \sum\_{i=1}^{N} (u\_i - l\_i) \tag{12}$$

where *N* is the number of series data, *y*ˆ and *y*∗ represents the prediction value and the actual value, respectively. Δ*y*∗ is the difference between the maximum and minimum values of the actual series, *ui* and *li* represents the upper and lower boundary of 95% CI. Lower RMSEs and PINAWs indicate better accuracy and certainty of the prediction model.

#### 2.4.2. NSGA-II Algorithm

NSGA-II makes excellent improvements on NSGA, and it has faster computational efficiency and population diversity [37]. There are two basic concepts in NSGA-II: nondominated sorting and crowding distance. The procedure of non-dominated sorting begins with the identification of non-dominated solutions. As shown in Figure 3, for two members *mi* and *mj*, if all the objectives of *mi* are better than *mj*, *mi* is defined to dominate *mj*. Then, the members which are not dominated by others constitute the current front. Next, the members of the current front are removed, and the sorting is performed on the remained population. The procedure is repeated until all the members are distributed to different fronts.

**Figure 3.** Non-dominated sorting and crowding distance.

The crowding distance is used to measure the density of members. It is defined as the sum of the side length of the cuboid shown in Figure 3.

$$cd(m\_i) = \sum\_{k=1}^{K} \left| o\_k^{i+1} - o\_k^{i-1} \right| \tag{13}$$

where *K* is the dimension number of the objectives. Selecting members with a high crowding distance can improve the diversity of the population.

The brief schematic of NSGA-II is illustrated in Figure 4, and the main steps are as follows:

**Step 1**: The population is initialized. Then, the offspring population *Qt* is generated from the parental population *Pt* through crossover and mutation operations. The population sizes of *Pt* and *Qt* are both *N*.

**Step 2**: The *Pt* and *Qt* are merged to form the *Rt*. The non-dominated sorting is performed on *Rt*, and a series of fronts *Fi* are obtained.

**Step 3**: The fronts are selected in sort order to form the *Pt*+1 until the population size of *Pt*+1 exceeds *N*. The members of the last front *Fl* are sorted by the crowding distance. The members are selected in order until the population size of *Pt*+1 equals *N*.

**Step 4**: The above steps are repeated until the maximum number of evolution is reached. Then the first front is selected as the Pareto front.

**Figure 4.** Schematic diagram of the NSGA-II.

#### 2.4.3. MOGPR Model

The GPR and the NSGA-II are coupled to construct the MOGPR model. The threedimensional decision vector is formed by the length scales of RBF, MA, and RQ in the GPR model. The multiple objectives include both RMSE and PINAW. The schematic diagram of the MOGPR is illustrated in Figure 5. The GPR model sets parameter values according to the decision vector generated by NSGA-II and calculates the prediction results. Then the results are compared with the actual values to obtain multi objectives. NSGA-II updates population location according to the multi objectives. Through several epochs of evolution, the Pareto front of the optimal result is finally output.

**Figure 5.** Schematic diagram of the MOGPR.

#### **3. Engineering Application**

In this section, the long-term monitoring data of a large actual FTU is adopted to conduct the experiments. The basic information of the FTU and the data source are described first. Then the data cleaning based on DBSCAN is implemented on the raw running data set. Next, the HSM of the FTU under variable operation conditions is constructed based on GMM. Finally, the performance degradation trend of the FTU is predicted by the proposed MOGPR.

#### *3.1. Data Description*

The studied FTU is located in the Dadu River basin, Sichuan, China. Its basic parameters are listed in Table 1. A set of PSTA-2100 state monitoring systems, including on-site monitoring cabinet and upper computer system, is configured on the FTU. The

on-site monitoring cabinet is formed by a sensor power module, data acquisition module, synchronous clock module, and industry cabinet. It is located near the FTU. The top cover vibration (*V*) is monitored by the acceleration sensor (PCB 352A60). The water head (*H*) and the active power (*P*) are obtained through communication with the supervisory control system of the plant following the modbus 485 protocol. The collected data are transmitted to the upper computer through the communication link following TCP/IP protocol. The actual length of the communication link, which consists of optical cables, switches, routers, repeaters, and transceivers, is above 1000 m. The analyzed sample set (*H*, *P*, *V*) is constructed by the data records exported from the upper computer. The overall structure of the FTU and the data acquisition system are illustrated in Figure 6. To sum up, the data acquisition process involves a series of devices and processes, which are prone to data loss and data abnormality.

**Table 1.** Basic parameters of the FTU.


**Figure 6.** Structure of the FTU and the data acquisition system.

The exported data includes 104,454 historical records from 20 January 2019 to 11 October 2019 for a total of 264 days, as shown in Figure 7. The on-site dataset is quite different from the high-quality data obtained in the ideal experimental environment, such as the bearing life cycle dataset published by the Center for Intelligent Maintenance Systems,

University of Cincinnati [38]. Most experimental data sets are less disturbed by external interference and have stable sampling frequency. However, it can be seen that the data loss and data anomaly are very obvious in the raw on-site data set. Data loss leads to unequal time stamps of samples. As shown in Figure 8, the time intervals between adjacent samples ranged from 10 s to 10<sup>3</sup> s, which makes it difficult to analyze the frequency domain features of signals. In addition, because of the task of grid regulation and seasonal changes in the water head, the operation condition parameters of the FTU change frequently and drastically. Low data quality and variable operation conditions greatly influence the performance state evaluation and degradation trend prediction of the FTU.

**Figure 7.** On-site measured data: (**a**) Water head; (**b**) Active power; (**c**) Vibration amplitude of top cover.

**Figure 8.** Frequency distribution histogram of logarithms of time intervals.

#### *3.2. Data Cleaning*

The raw running sample set (*H*, *P*, *V*) is exhibited in Figure 9. Due to the characteristic of the FTU, there is a specific limited operation region within the operating conditions. In the limited area, the stability and efficiency of the FTU decrease. In practice, the FTU is avoided from working in the limited area, so the sample points are concentrated in two regions. Because of various interference factors, the data anomaly is very obvious in the raw data set. The data anomaly mainly includes singular points whose values deviate significantly from the normal level and outliers whose values are within the normal scope, but their distribution deviates from the valid samples. To reduce the impact of data anomaly, DBSCAN was adopted to identify the singulars and outliers in the original sample set. The *ε* and *M* were set as 7.5 and 200, respectively. The processing result is shown in Figure 10. The two sample concentrated regions are marked as cluster 1 and cluster 2. Singulars and outliers are both marked as anomaly data automatically. Since the DBSCAN considers the distribution density of samples rather than the value of data, it performs well in identifying both outliers and singularities in the low-quality data set.

**Figure 9.** Raw running sample set.

**Figure 10.** Processing result of DBSCAN.

#### *3.3. Performance State Evaluation of the FTU*

Routine maintenance was implemented before 20 January 2019 on the studied FTU, and the restart test run performed normally. Moreover, as shown in Figure 7, the period from 20 January to 6 May 2019 includes most of the possible operating conditions of the FTU. Therefore, a total of 58,175 valid data points during this period were selected to construct the healthy sample set.

The GMM was adopted to fit the three-dimensional PDD function of the healthy sample set for the HSM construction of the FTU, as illustrated in Figure 11a, where different colors represent different Gaussian components of the HSM. On this basis, the NLLP was calculated as the PDI of the FTU. The distribution of PDIs on the *H* = 125 m section and *P* = 110 MW section of the HSM is shown in Figure 11b,c. It shows that there is a complex mapping relationship between PDIs and operation condition parameters. The PDIs are relatively lower in regions with dense distribution of healthy samples.

**Figure 11.** HSM based on GMM (**a**) Result of GMM clustering (**b**) PDI distribution on *H* = 125 m section; (**c**) PDI distribution on *P* = 110 MW section.

The valid data after 6 May 2019 were selected to construct the evaluation sample set. All 300 evaluation samples were input into the HSM as a group to calculate the PDIs. The timestamp corresponding to the last evaluation sample in each group was taken as the evaluation time. The obtained performance degradation curve, including 135 PDIs, is demonstrated in Figure 12. Although the PDI curve fluctuates locally due to low data quality, the curve has an obvious upward trend reflecting the performance degradation process of the FTU.

**Figure 12.** Performance degradation curve of the FTU.

#### *3.4. Degradation Trend Prediction of the FTU*

To forecast the performance degradation trend of the FTU, a rolling prediction model based on MOGPR was established. Every seven PDIs, as the time window, were input into the prediction model, and the next PDI was predicted. The sliding step was set as one PDI point. As the time window slides, the prediction model is constantly updated.

The RBF, MA, and RQ were chosen to construct the kernel function of the GPR model. The length scales *lRBF*, *lMA* and *lRQ* were optimized with the NSGA-II according to two objectives of RMSE and PINAW. The search scopes of the above three length scales were all set as [10−5, 1]. The parameter configuration of the NSGA-II is listed in Table 2. The prediction result that is closest to the origin in the Pareto front of multiobjective optimization is exhibited in Figure 13. The proposed MOGPR takes both accuracy and certainty of interval prediction into account. The maximum error does not exceed 0.47, and the maximum width of CI does not exceed 3.41.

**Table 2.** Parameter configuration of the NSGA-II.


**Figure 13.** Prediction result of MOGPR: (**a**) Prediction result; (**b**) Width of CI; (**c**) Error distribution.

#### **4. Comparison Analyses**

To validate the effectiveness of the proposed approach, different methods of data cleaning, HSM construction, and trend prediction were compared. The experiments were conducted by Python 3.7.9 in a calculation station with an Intel Core i9-10900K CPU, an NVIDIA GeForce RTX 2080 Super GPU, and 64 GB RAM.

#### *4.1. Comparison of Different Data Cleaning Methods*

To validate the effectiveness of the proposed data cleaning method, the DBSCAN and 3-sigma principle are compared in this section. Different data sets were adopted to establish HSMs, including the raw data set Φ0, the 3-sigma-processed data set Φ1, and the DBSCAN-processed data set Φ2. The processing result of the 3-sigma principle is shown in Figure 14. Compared with Figure 10, it can be discovered that the 3-sigma principle can recognize partial singularities but can hardly identify outliers.

Then the PDI curves were calculated in the same way as Section 3.3. The *STD* and *S* were adopted as the criteria of stability and smoothness, defined as follows:

$$STD = \sqrt{\frac{1}{N} \sum\_{i=1}^{N} \left(I\_i - \overline{I}\right)^2} \tag{14}$$

$$S = \frac{\sum\_{i=1}^{N-1} |I\_{i+1} - I\_i|}{N - 1} \tag{15}$$

where *I* represents the PDI and *N* is the number of PDIs in the curve. Smaller *S* and *STD* indicate that the curve is smoother and the represented performance degradation trend is clearer.

**Figure 14.** Processing result of 3-sigma principle.

The comparison results are shown in Table 3 and Figure 15. It can be seen that:


**Table 3.** Comparison of different data cleaning methods.


**Figure 15.** Comparison of different data cleaning methods.

#### *4.2. Comparison of Different Methods for HSM Construction*

To compare the robustness of HSMs towards data loss, different ratios (*r*) of randomly selected healthy samples were adopted to build HSMs based on GMM, SVM, and BPNN, respectively, *r* ∈ [0.2, 0.4, 0.6, 0.8, 1]. The major parameters of the models mentioned above are listed in Table 4. Then the entire healthy sample set was input into the HSMs to calculate PDIs. Each experiment was repeated 100 times, and the STDs of results were calculated. In order to make the PDI results of different HSMs comparable, the maximum-minimum normalization was implemented to transform the value range of PDIs to [0, 1]. Unlike GMM, SVM and BPNN establish the mapping relationship between operation condition parameters and vibration amplitude values. Thus the PDI of HSMs based on SVM and BPNN was defined as:

$$PDI = \frac{1}{N} \sum\_{i=1}^{N} \left| V^\* - \hat{V}(H^\*, P^\*) \right| \tag{16}$$

where *V*∗, *H*∗, and *P*∗ represent the real values of vibration amplitude, water head, and active power. *V*ˆ denotes the mapping value outputted by the HSM.


**Table 4.** Parameter setting of compared models.

The comparison results are illustrated in Figure 16. As the sample quantity decreases, the results' STD based on GMM increases slowly when *r* is greater than 0.6. The results' STDs based on SVM and BPNN are relatively higher and have no significant relationship with the data quantity of HSM construction. Since the HSM based on GMM looks at the PDD of samples rather than the specific values of individual samples, the actual distribution of the population can be estimated from the sample distribution as long as the total quantity of samples is enough. Therefore, the proposed HSM is less affected by data loss.

**Figure 16.** Comparison of different HSM.

*4.3. Comparison of Different Prediction Methods*

To verify the effectiveness of the proposed MOGPR, different prediction models are compared, including two multiobjective optimization machine learning models named MOQRLSTM and MOQRNN, two machine learning models without parameter optimization named QRLSTM and QRNN [39], three nonparametric statistical regression models named locally weighted linear regression (LWLR) [40] and kernel regression (KR) [41]. The loss functions of QRLSTM and QRNN are defined as quantile regression functions [42], which makes them predict conditional quantiles. According to the two objectives of RMSE and PINAW, the node numbers of each hidden layer of QRLSTM and QRNN are optimized by NSGA-II to construct the MOQRLSTM and MOQRNN. The NSGA-II is set as the same as Section 3.4. LWLR and KR are widely used to predict time series data with trend terms of different scales because their algorithms are simple, and the time costs are low. The main parameters of the mentioned models are listed in Table 5.

**Table 5.** Parameter setting of compared models.


<sup>1</sup> For the three multiobjective optimization prediction models, some parameters are optimized in a range. For the five prediction models without parameter optimization, the parameters are certain values.

For a more complete performance comparison of the various models, in addition to the two objectives of RMSE and PINAW, the other fourcriteria named prediction interval coverage probability (PICP), mean absolute percentage error (MAPE), R2\_score, and probability integral transform (PIT) are introduced [43]. The RMSE represents the accuracy of prediction results; the PINAW indicates the certainty of the prediction interval. Their

definitions are given in Equations (11) and (12), respectively. The MAPE is complementary to RMSE and reflects the relative accuracy of the prediction, defined as:

$$MAPE = \frac{1}{N} \sum\_{i=1}^{N} \left| \frac{\mathcal{Y}\_{i} - \mathcal{Y}\_{i}^{\*}}{\mathcal{Y}\_{i}^{\*}} \right| \times 100\text{\%} \tag{17}$$

The R2\_score represents the proportion that can be explained by the fitted regression relationship in the overall change trend of the dependent variable, defined as:

$$R2\\_score = \frac{\sum\_{i=1}^{N} (\mathcal{Y}\_i - \overline{\mathcal{Y}})^2}{\sum\_{i=1}^{N} (\mathcal{Y}^\*\_i - \overline{\mathcal{Y}})} \tag{18}$$

where *y* is the mean of the actual value *y*∗. A R2\_score close to 1 indicates that the prediction model is more reliable.

The PICP represents the probability that the predicted CI contains the observed values, expressed as:

$$\begin{array}{c} PICP = \frac{1}{N} \sum\_{i=1}^{N} f\_i\\ f\_{\bar{i}} = \left\{ \begin{array}{c} 0, y\_{\bar{i}}^\* \notin [l\_{\bar{i}}, u\_{\bar{i}}]\\ 1, y\_{\bar{i}}^\* \in [l\_{\bar{i}}, u\_{\bar{i}}] \end{array} \right. \end{array} \tag{19}$$

where *fi* indicates whether the observed value is within the predicted CI. A higher PICP indicates the interval prediction result is more accurate.

The PIT values reflect the correlation between the prediction distribution and the actual distribution, defined as:

$$PIT = \int\_{-\infty}^{y^\*} \hat{P}(y) \, \mathrm{d}y \tag{20}$$

where *P*ˆ is the prediction distribution and *y*<sup>∗</sup> is the observed value. If the distribution of PIT values is similar to the uniform distribution, the prediction reliability is better.

The Pareto fronts of optimization results of MOGPR, MOQRLSTM, and MOQRNN are illustrated in Figure 17. The results of MOGPR are generally closer to the origin. The lowest PINAW values reached by the 3 models are similar (about 0.243~0.271), while the lowest RMSE value obtained by MOGPR (about 0.199) is smaller than MOQRLSTM's (about 0.262) and MOQRNN's (about 0.273). As two independent optimization objectives, RMSE and PINAW have the same dimension of quantity and similar value range. In this study, the importance of the two objectives is considered equal. So the result closest to the origin in the Pareto front is selected as the final result of three multiobjective optimization prediction models.

The multi metrics of different models are listed in Table 6, including RMSE, PINAW, PICP, MAPE, R2\_score, and cost time. The prediction results are shown in Figure 18, and the quantile-quantile (QQ) plots of PITs are illustrated in Figure 19. The areas between the purple dotted lines in QQ plots represent PITs that follow the uniform distribution at a 5% significance level. The analysis of the results is as follows:


**Figure 17.** Pareto fronts of optimization results.

**Table 6.** Results of comparison experiments.


**Figure 18.** Results of interval prediction: (**a**) MOGPR; (**b**) MOQRLSTM; (**c**) MOQRNN; (**d**) GPR; (**e**) QRLSTM; (**f**) QRNN; (**g**) LWLR; (**h**) KR.

**Figure 19.** *Cont*.

**Figure 19.** QQ plot of PITs: (**a**) MOGPR; (**b**) MOQRLSTM; (**c**) MOQRNN; (**d**) GPR; (**e**) QRLSTM; (**f**) QRNN; (**g**) LWLR; (**h**) KR.

In summary, among the eight compared models, MOGPR has better accuracy, certainty, and reliability of probability prediction. Besides, the calculation effectiveness of MOGPR is higher than other multiobjective optimization prediction models.

#### **5. Conclusions**

In this study, aiming at the three major practical engineering difficulties of low data quality, complex variable operation conditions, and prediction model parameter optimization, an ensemble prognostic method of FTUs using low-quality data under variable operation conditions is proposed. Firstly, to comprehensively reflect the performance of the FTU under complex operation conditions, the running data set is constructed by combining operation condition parameters and monitoring data. Secondly, to reduce the impact of anomaly data, the DBSCAN is adopted to clean both outliers and singulars in the raw running data set. Thirdly, based on the GMM and the probability theory, the HSM is established, which improves the robustness of the evaluation model against data missing and data randomness. Fourthly, the MOGPR is proposed to predict the performance degradation trend with a confidence interval and to automatically optimize model parameters on both accuracy and certainty. Finally, a series of comparison experiments were implemented on practical data set from an actual large FTU.

The experimental results demonstrate that: (1) The data cleaning approach based on DBSCAN performs better in identifying both outliers and singularities. The stability and smoothness of PDI curves are improved by 3.2 times and 1.9 times, respectively, by DBSCAN compared with 3-sigma. (2) Compared with SVM and BPNN, the HSM based on GMM has better robustness against data loss. (3) The proposed MOGPR has better accuracy, certainty, and reliability of probability prediction. The root-mean-squared error, the prediction interval normalized average, the prediction interval coverage probability, the mean absolute percentage error, and the R2 score of the proposed method achieved 0.223, 0.289, 1.000, 0.641%, and 0.974, respectively.

Thus, the proposed method can be applied to the performance evaluation and degradation trend prediction of FTUs in a practical engineering environment. In further research, if the running data across multiple maintenance periods are available, the corresponding relationship between PDI values and actual state can be established according to maintenance records. On this basis, the multi-level degradation state of the FTU can be divided, and corresponding maintenance strategies can be formulated to provide technical support for state-based maintenance.

**Author Contributions:** Conceptualization, R.D.; Data curation, P.W.; Formal analysis, J.L.; Methodology, R.D. and J.Z.; Project administration, R.D.; Validation, J.L.; Visualization, W.L.; Writing—original draft, R.D.; Writing—review & editing, J.Z., P.W. and W.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Natural Science Foundation of China (No. 72171096) and the Fundamental Research Funds for the Central Universities (HUST: 2020JYCXJJ046).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**


#### **References**


## *Article* **Dynamic Behavior Analysis and Stability Control of Tethered Satellite Formation Deployment**

**Kangyu Zhang 1, Kuan Lu 1,\*, Xiaohui Gu 2, Chao Fu <sup>1</sup> and Shibo Zhao <sup>1</sup>**


**Abstract:** In recent years, Tethered Space Systems (TSSs) have received significant attention in aerospace research as a result of their significant advantages: dexterousness, long life cycles and fuel-less engines. However, configurational conversion processes of tethered satellite formation systems in a complex space environment are essentially unstable. Due to their structural peculiarities and the special environment in outer space, TSS vibrations are easily produced. These types of vibrations are extremely harmful to spacecraft. Hence, the nonlinear dynamic behavior of systems based on a simplified rigid-rod tether model is analyzed in this paper. Two stability control laws for tether release rate and tether tension are proposed in order to control tether length variation. In addition, periodic stability of time-varying control systems after deployment is analyzed by using Floquet theory, and small parameter domains of systems in asymptotically stable states are obtained. Numerical simulations show that proposed tether tension controls can suppress in-plane and out-ofplane librations of rigid tethered satellites, while spacecraft and tether stability control goals can be achieved. Most importantly, this paper provides tether release rate and tether tension control laws for suppressing wide-ranging TSS vibrations that are valuable for improving TSS attitude control accuracy and performance, specifically for TSSs that are operating in low-eccentricity orbits.

**Keywords:** tethered satellite formation; dynamic behavior; control; stable deployment; Floquet theory

#### **1. Introduction**

In recent years, satellite development has rapidly increased worldwide [1,2]. In particular, Tethered Space Systems (TSSs) have received significant attention in aerospace research as a result of their significant advantages: dexterousness, long life cycles and fuel-less engines [3–5]. TSSs are a new class of space vehicle that join two or more spacecraft together into a single structure by using soft tethers [6,7]. TSSs are utilized in man-made microgravity environments [8,9], spacecraft orbit transfers [10,11] and space debris cleanup [12,13], and they exhibit stronger reliability, higher stability and more diversified functions [14,15] when compared to traditional satellites. Tethered satellites are in unstable states during deployment without effective control as a result of disturbances produced by the space environment [16,17]. Due to their structural peculiarities, gravitational forces, aerodynamic drag, solar radiation pressure and other disturbances produced by the special environment in outer space, TSS vibrations are easily produced. These types of vibrations are extremely harmful to spacecraft.

TSS dynamics and control aspects have received considerable attention in recent decades [18–20]. Rigid-rod tether models provide analytical solutions of TSS and are widely used in basic research. For example, Williams [21] proposed a new feedback control scheme in which electrodynamic tether vibrations were suppressed, and tether stability was effectively controlled during deployment by using only electric current modulations.

**Citation:** Zhang, K.; Lu, K.; Gu, X.; Fu, C.; Zhao, S. Dynamic Behavior Analysis and Stability Control of Tethered Satellite Formation Deployment. *Sensors* **2022**, *22*, 62. https://doi.org/10.3390/s22010062

Academic Editor: Mohammad N Noori

Received: 12 October 2021 Accepted: 16 December 2021 Published: 23 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

In another study, interorbital rendezvous with small relative inclination was also analyzed, and a nonlinear receding horizon controller was considered for tracking highly nonlinear systems by producing disturbances in system mass distributions and perturbations to initial system conditions [22]. Pradeep and Kumar [23] proposed nonlinear feedback tension control laws based on Liapunov's method, used a linear state variable feedback control and affirmed the desired length of extended tethers in a reasonable amount of time.

Stability analysis is a core research focus in mechanism studies of dynamic systems. The Floquet theory is a stability theory of solutions of linear ordinary differential equations with periodic variable coefficients [24] that was proposed by G. Floquet in 1868. Few researchers used the Floquent theory to study TSS stability and dynamic behavior. Yu et al. [25] analyzed the spinning stability of a three-body Tethered Satellite Formation (TSF) by using Floquet theory, and stability analysis indicated that unstable motion occurs if its spinning angular rate is less than the critical value |−2.8| or 0.65 times its orbital angular rate. In another study, an analytical tether length rate control was designed, and parameter regions for stable deployment in order to maintain a tensile tether state were obtained [26]. Ellis and Hall [27] analyzed the stability of out-of-plane vibrations of a spinning TSS, and two satellites as point masses were connected by a rigid rod, constraining the system's mass center to a circular orbit.

The aforementioned references indicate that little attention has been given to numerical studies on the accuracy of simplified TSS models, control stability and the influence of orbital eccentricity. However, tethers produce in-plane and out-of-plane swings and longitudinal and transverse vibrations as a result of complex perturbations [28,29], and a slight change in orbital eccentricity can significantly affect the original system. Previous studies show that analytical solutions for complex nonlinear models of TSS are difficult to obtain. This paper evaluates the stability of periodic TSS motions by using Floquet theory and provides control laws and small parameter domains of stability based on a simplified rigid-rod tether model.

In this paper, nonlinear dynamic behavior and stability of TSS during deployment are analyzed. A simplified rigid-rod model of a two-body tethered satellite is described in Section 2. Two simplified models of three-DOF equations are discussed in Section 3. Two control laws of tether release rate and tether tension are proposed in Section 4. The periodic stability of time-varying control systems is analyzed by using Floquet theory in Section 5. Conclusions are discussed in Section 6.

#### **2. TSS Equations of Motion Using Lagrangian Method**

Since the 1990s, TSS theory has developed rapidly from rigid-rod models to bead models and from continuous models to discrete models, and model accuracy is constantly improving [30,31]. However, other models' dynamics equations are more complex and the quantity of computation is larger when compared with rigid-rod models [32]. The classical rigid-rod model is widely used. As shown in Figure 1, a rigid tethered satellite system consists of a mother satellite, *m*1, and a subsatellite, *m*2, (mass points), respectively. Both satellites are connected by a rigid tether. In this model, regardless of tether flexibility and tether elasticity, tethers released into outer space are considered as straight rods of infinite stiffness that do not bend or twist. Tether wounds on a spool of a deployment device in the mother satellite and tether length can be effectively controlled by the device.

The inertial geocentric frame, *O* − *XYZ*, the orbital frame, *o* − *xyz*, and the tether body frame, *o* − *x*b*y*b*z*b, are all radial–transversal out-of-plane frames, and these are established in order to describe TSS position and attitude in Figure 1. Orbital radius and tether length are expressed as *R*, *l*, respectively; right ascension and declination to the center of mass are expressed as *α*, *δ*, respectively; and the tether in-plane angle and tether out-of-plane angle are expressed as *θ*, *φ*, respectively.

**Figure 1.** Rigid tether dynamical model.

The radius vector with respect to the center of mass is written in inertial coordinates as follows.

$$\mathbf{R} = R\cos\delta\cos a\mathbf{i} + R\cos\delta\sin a\mathbf{j} + R\sin\delta\mathbf{k}.\tag{1}$$

Total kinetic energy, *Tk*, consists of the translation of the center of mass, *Tt*, system rotation, *Tr*, and tether deployment, *Te*:

$$T\_l = \frac{1}{2}m\left(\dot{\vec{R}}^2 + R^2\dot{\delta}^2 + R^2\dot{\alpha}^2\cos^2\delta\right) \tag{2}$$

where *m* = *m*<sup>1</sup> + *m*<sup>2</sup> + *mt* is the total system mass, *mt* = *ρl* = *m*<sup>0</sup> <sup>1</sup> − *m*<sup>0</sup> is the tether mass released into the external environment, *ρ* is the tether linear density and *m*<sup>0</sup> <sup>1</sup> is the mass of the mother satellite before tether deployment that includes the tether mass:

$$T\_r = \frac{1}{2} \{\boldsymbol{\omega}\}^T [\boldsymbol{I}] \{\boldsymbol{\omega}\} \tag{3}$$

where *I* is the tensor matrix of moment of inertia of the tethered satellite system, and *ω* is the inertial angular velocity of the tether in the inertial frame [33], which can be expressed as follows:

$$\begin{aligned} \boldsymbol{\omega} &= \left( \dot{a} \sin \delta \cos \theta \cos \phi + \dot{\theta} \sin \phi - \dot{\delta} \sin \theta \cos \phi + \dot{a} \cos \delta \sin \phi \right) \mathbf{i} \\ &- \left( \dot{\phi} + \dot{a} \sin \delta \sin \theta + \dot{\delta} \cos \theta \right) \mathbf{j} \\ &+ (\dot{\theta} \cos \phi - \dot{a} \sin \delta \cos \theta \sin \phi + \dot{\delta} \sin \theta \sin \phi + \dot{a} \cos \delta \cos \phi) \mathbf{k} \end{aligned} \tag{4}$$
 
$$T\_{\mathbf{r}} = \frac{1}{2} m^{\star \mathbf{l}} l^2 \left[ (\dot{\phi} + \dot{a} \sin \delta \sin \theta + \dot{\delta} \cos \theta) \right]^2 \tag{5}$$

$$\begin{aligned} T\_r &= \frac{1}{2} m^\* l^2 \Big[ \left( \dot{\phi} + \dot{a} \sin \delta \sin \theta + \dot{\delta} \cos \theta \right)^\ast \\ &+ \left( \dot{\delta} \sin \theta \sin \phi - \dot{a} \sin \delta \cos \theta \sin \phi + \dot{a} \cos \delta \cos \phi + \dot{\theta} \cos \phi \right)^2 \end{aligned} \tag{5}$$

where *m*\* = (*m*1+*mt* /2)(*m*2+*mt* /2)/*m* − *mt*/6 is the reduced mass of the system.

$$T\_{\varepsilon} = \frac{1}{2} \frac{m\_1 (m\_2 + m\_l)}{m} \dot{l}^2. \tag{6}$$

The tether is assumed to be stationary relative to the mother satellite within the winch control mechanism, and its speed is provided during deployment.

When TSS systems are active in space, they are still within Earth's gravitational field, and their potential energy is caused by Earth's attraction. Potential energy is obtained from the mother satellite, the subsatellite and the tether, which is simplified by taking the

first term of Maclaurin's series expansion. We assume that a spherical earth is considered as follows:

$$V = -\frac{\mu\_\varepsilon m}{R} + \frac{\mu\_\varepsilon m^\ast l^2}{2R^3} (1 - 3\cos^2\phi \cos^2\theta) \tag{7}$$

where *μ<sup>e</sup>* = 398, 600 km3/s2 is Earth's gravitational coefficient. The Lagrange function can be formed as follows.

$$L = T\_l + T\_r + T\_\varepsilon - V.\tag{8}$$

By substituting Equation (8) into Lagrange's equations and by assuming that the system's center of mass is running in a constant orbital plane (*δ* = 0), the system's equations of motion can be obtained as follows:

$$m\ddot{R} - mR\dot{\alpha}^2 + \frac{\mu m}{R^2} - \frac{3\mu m^\*l^2}{2R^4}(1 - 3\cos^2\theta\cos^2\phi) = Q\_R\tag{9}$$

$$2m\dot{R}\dot{R}\dot{a} + \dot{m}^\* l^2 (\dot{a} + \dot{\theta}) \cos^2 \phi + 2m^\* l \dot{l} [(\dot{a} + \dot{\theta}) \cos^2 \phi] \tag{10}$$
 
$$\dots \gamma\_{\epsilon \times \cdots} \quad \dots \quad \gamma\_{\epsilon \times \cdots} \quad \dots \quad \vdots \quad \dots \quad \dots \quad \dots \tag{10}$$

$$\begin{aligned} &+m^{\*}l^{2}[(\ddot{\alpha}+\ddot{\theta})\cos^{2}\phi-2(\dot{\alpha}+\dot{\theta})\dot{\phi}\sin\phi\cos\phi]+mR^{2}\ddot{\alpha}=Q\_{d} \end{aligned} \tag{10}$$
 
$$\begin{aligned} \left[\begin{array}{cccc} \dots \gamma\epsilon & \ddot{\gamma} & \ddots & \ddots & \ddots & \cdots & \dot{\gamma} & \dot{\alpha} \end{array}\right] \end{aligned} \tag{10}$$

$$\begin{aligned} m^\*l^2[(\ddot{\alpha} + \ddot{\theta})\cos^2\phi - (\dot{\alpha} + \dot{\theta})\dot{\phi}\sin 2\phi] + 2m^\*l\dot{l}(\dot{\alpha} + \dot{\theta})\cos^2\phi\\ + \dot{m}^\*l^2(\dot{\alpha} + \dot{\theta})\cos^2\phi + \frac{3\mu m^\*l^2}{2R^3}\sin 2\theta\cos^2\phi = Q\_\theta \end{aligned} \tag{11}$$

$$\dot{m}^{\ast\ast}l^{2}\dot{\phi} + 2m^{\ast}l\dot{l}\dot{\phi} + m^{\ast}l^{2}\ddot{\phi} + \frac{1}{2}m^{\ast}l^{2}(\dot{a} + \dot{\theta})^{2}\sin 2\phi + \frac{3\mu m^{\ast}l^{2}}{2R^{3}}\sin 2\phi\cos^{2}\theta = Q\_{\phi} \tag{12}$$

$$\begin{aligned} \dot{m}^\# \dot{l} + m^\# \dot{l} - \frac{1}{2} (m^\star)' l^2 \left[ \dot{\phi}^2 + (\dot{a} + \dot{\theta})^2 \cos^2 \phi \right] - m^\star l \left[ \dot{\phi}^2 + (\dot{a} + \dot{\theta})^2 \cos^2 \phi \right] \\\ \dot{m}^\star \ddot{l} - \frac{1}{2} m^\star \left[ \dot{a} + \dot{a} + \dot{a}^\star \dot{a} + \dot{a}^\star \dot{a}^\star \right] \end{aligned} \tag{13}$$

$$\left(-\frac{1}{2}(m^{\#})\dot{l}^2 + \frac{\mu(m^{\star})l^2}{2\hbar^3}(1 - 3\cos^2\theta\cos^2\phi) + \frac{\mu m^{\star}l}{\hbar^3}(1 - 3\cos^2\theta\cos^2\phi) = Q\_l$$

where ( ) = <sup>d</sup>( )/d*l*, *<sup>m</sup>*\* = *mt*(3*m*<sup>1</sup> − <sup>3</sup>*m*<sup>2</sup> − *<sup>m</sup>*)/(6*m*), *<sup>m</sup>*# = *mt*(2*m*<sup>1</sup> − *<sup>m</sup>*)/*m*, *Ql* = −*<sup>T</sup>* is the tension control and the generalized forces, *Q<sup>θ</sup>* and *Qφ*, are typically assumed to be negligible as a result of distributed forces along the tether. It should be noted that tether length can be controlled by deployment/retrieval of the winch control mechanism in the mother satellite; therefore, *m*1, *mt* in Equations (9)–(13) are functions of tether length.

Based on the premise of a Keplerian reference orbit for the center of mass as the independent variable, the orbit's true anomaly, *ν*, is used to replace the generalized coordinate, *α*, in order to withdraw the premise, which can be expressed as follows:

$$\dot{\nu} = \sqrt{\frac{\mu\_{\varepsilon}}{a^3 (1 - e^2)}} \kappa^2 \text{ , } R = \frac{a(1 - e^2)}{\kappa} \tag{14}$$

where *e* is the orbital eccentricity, and *a* is the semi-major axis of the orbit, *κ* = 1 + *e* cos *ν*. d(*i*) <sup>d</sup>*<sup>t</sup>* <sup>=</sup> <sup>d</sup>(*i*) <sup>d</sup>*<sup>ν</sup>* <sup>×</sup> <sup>d</sup>*<sup>ν</sup>* <sup>d</sup>*<sup>t</sup>* <sup>⇒</sup> . *i* = *i* . *ν*, .. *i* = *i* . *ν* <sup>2</sup> <sup>+</sup> .. *νi* . *<sup>ν</sup>* , *<sup>i</sup>* = *<sup>θ</sup>*, *<sup>φ</sup>*, *<sup>l</sup>* is utilized by Equations (9)–(13). Nondimensional equations of motion can be written as follows:

$$\theta^{\prime\prime} = 2(\theta^{\prime} + 1)[\frac{e\sin\nu}{\kappa} + \phi^{\prime}\tan\phi - \Phi\_1 \frac{\Lambda^{\prime}}{\Lambda}] - \frac{\mathfrak{Z}}{2\kappa}\sin 2\theta \tag{15}$$

$$\phi'' = \frac{2e\sin\upsilon}{\kappa}\phi' - 2\Phi\_1\frac{\Lambda'}{\Lambda}\phi' - \frac{1}{2}[\left(\theta' + 1\right)^2 + \frac{3}{\kappa}\cos^2\theta]\sin 2\phi \tag{16}$$

$$\begin{split} \Lambda^{\prime\prime} &= \frac{2\varepsilon \sin \upsilon}{\kappa} \Lambda^{\prime} - \Phi\_2 \frac{\Lambda^{\prime^2}}{\Lambda} + \Phi\_3 \Lambda \left[ \phi^{\prime^2} + (\theta^{\prime} + 1)^2 \cos^2 \phi \right. \\ &\left. + \frac{1}{\kappa} (3 \cos^2 \theta \cos^2 \phi - 1) \right] - \frac{mT}{m\_1 \dot{\upsilon}^2 L (m\_2 + m\_1)} \end{split} \tag{17}$$

where ( ) = d( )/d*ν*, Λ = *l*/*L* is the nondimensional tether length, and *L* is the reference tether length. Φ*i*, *i* = 1, 2, 3 is the nondimensional coefficient.

$$\Phi\_1 = \frac{m\_1(m\_2 + m\_l/2)}{mm^\*}, \ \Phi\_2 = \frac{(2m\_1 - m)m\_l}{2m\_1(m\_2 + m\_l)}, \ \Phi\_3 = \frac{m\_2 + m\_l/2}{m\_2 + m\_l} \tag{18}$$

#### **3. Dynamic Analysis of Simplified Models of Single-DOF and Two-DOFs**

In order to clarify dynamic behavior mechanisms of TSS and to explore the influence of various parameters on dynamic responses, Equations (15)–(17) need to be simplified, assuming that the tether length is fixed when system configurations remain fixed.

#### *3.1. Single-DOF (θ)*

Numerous studies show that tether in-plane angles are much larger than out-of-plane angles; when TSSs are operating in orbital planes, *φ* = 0. In this case, Equations (15)–(17) can be rewritten as follows:

$$\theta'' = 2(\theta' + 1)\frac{e\sin\nu}{\kappa} - \frac{3}{\kappa}\sin\theta\cos\theta \tag{19}$$

assuming that a spherical earth is considered. Orbital eccentricity, *e*, is a small quantity; thus, the perturbation method was selected in order to calculate approximate analytical solutions for Equation (19), where *e* is regarded as a tiny perturbation that is substituted into Equation (19).

$$
\theta'' - 2(\theta' + 1)\frac{e\sin\nu}{1 + e\cos\nu} + \frac{3\sin 2\theta}{2(1 + e\cos\nu)} = 0\tag{20}
$$

The power series form of the periodic solution can be written as follows.

$$\theta\_p(\nu, e) = e \cdot \theta\_1(\nu) + e^2 \cdot \theta\_2(\nu) + e^3 \cdot \theta\_3(\nu) + e^4 \cdot \theta\_4(\nu) + e^5 \cdot \theta\_5(\nu) \tag{21}$$

The linear ordinary differential equation is written as follows.

$$\begin{aligned} \theta'' &= 3\sin\upsilon\\ \theta''\_2 + 3\theta\_2 &= 2\theta'\_1\sin\upsilon - \theta''\_1\cos\upsilon,\ \theta''\_3 + 3\theta\_3 = 2\theta'\_2\sin\upsilon - \theta''\_2\cos\upsilon\\ \theta''\_4 + 3\theta\_4 &= 2\theta'\_3\sin\upsilon - \theta''\_3\cos\upsilon,\ \theta''\_5 + 3\theta\_5 = 2\theta'\_4\sin\upsilon - \theta''\_4\cos\upsilon \end{aligned} \tag{22}$$

Equation (22) can be executed by the periodic initial condition, *θi*(0, *θi*0) = *θi*(2*π*, *θi*0)*i* = 1, 2, . . . , 5, and the analytical solution of Equation (19) can be expressed as follows.

$$\begin{aligned} \theta\_p &= e \sin \nu - \frac{3}{2} e^2 \sin 2\nu + e^3 \sin 3\nu - \frac{3}{25} e^4 (5 \sin 4\nu + 13 \sin 2\nu) \\ &+ \frac{3}{143} e^5 (143 \sin \nu + 55 \sin 3\nu - 5 \sin 5\nu) \end{aligned} \tag{23}$$

Figure 2a shows the tether in-plane vibration angle versus the orbit's true anomaly expressed in radians with various *e*. The system moves periodically and repeatedly with a period of 2*π* in the direction of *θ*, and *e* = 0.1, *θ*max = 5.93◦ appears at 1/4 and 3/4 of the period, respectively. Figure 2b illustrates that *θ*max, *θ* max increase as orbital eccentricity increases *e* → 1 (the elliptical orbit is flatter), which means that the system tends to move towards an unstable equilibrium state.

**Figure 2.** (**a**) Tether in-plane angle, *θ*, versus *ν*; (**b**) angular velocity of the in-plane angle, *θ* , versus *θ*. *3.2. Two-DOFs (θ*, *φ)*

The tether out-of-plane vibration angle, *φ*, is considered based on a simplified single-DOF model. In this case, Equations (15)–(17) can be rewritten as follows.

$$\theta'' = 2(\theta' + 1)[\frac{e\sin\nu}{\kappa} + \phi'\tan\phi] - \frac{3}{2\kappa}\sin 2\theta \tag{24}$$

$$\phi'' = \frac{2e\sin\nu}{\kappa}\phi' - \left[\left(\theta' + 1\right)^2 + \frac{3}{\kappa}\cos^2\theta\right]\frac{1}{2}\sin 2\phi \tag{25}$$

Figure 3 shows tether in-plane and out-of-plane vibration angles, *θ*, *φ*, versus *ν*.

**Figure 3.** Tether vibration angles, *e* = 0.1: (**a**) tether in-plane angles versus *ν*; (**b**) tether out-plane vibration angles versus *ν*.

In Equations (24) and (25), orbital eccentricity is 0.1, and initial values of the single-DOF motion solution *θ*0, *θ* <sup>0</sup> are adopted by Equation (24). As shown in Figure 3, numerical simulations show that the out-of-plane vibration angle is relatively small, and the maximum of *φ* is 0.0247◦, far less than the in-plane angle, which has a slight effect on TSS dynamic response. Hence, the effect of the out-of-plane angle is negligible. However, coupling errors caused by out-of-plane vibrations to in-plane vibrations still require further numerical verification.

In order to further verify the accuracy of the single-DOF simplified model, different orbital eccentricity values are substituted into Equations (24) and (25). Both curves almost coincide in Figure 4a, and the maximum error of the in-plane vibration angle is 0.0945◦ with *e* = 0.1 in Figure 4b, which illustrates that coupling effects of the out-of-plane vibration angle are negligible. In particular, the solution of the single-DOF with the first five orders demonstrates sufficient accuracy in Equation (23), which proves that the error of the perturbation method is negligible. However, as shown in Figure 5, it can be easily observed that the error of the single-DOF simplified model increases as orbital eccentricity increases. In Figure 5, the error of the single-DOF simplified model was significantly smaller when value *e* decreased from 0.42 to 0.1, which means that a single DOF-simplified model can be applied to orbits with low orbital eccentricity so that accuracy can be guaranteed.

**Figure 4.** Coupling effects of *φ*, *e* = 0.1: (**a**) numerical comparison for *θ*; (**b**) absolute error for *θ*.

**Figure 5.** Coupling effects of *φ*: (**a**) *e* = 0.30; (**b**) *e* = 0.40; (**c**) *e* = 0.41; (**d**) *e* = 0.42.

**Remark 1.** *The accuracy of the simplified model of TSS is strongly influenced by orbital eccentricity. For low-eccentricity orbits, a simplified model of TSS can significantly reduce calculation time.*

#### **4. Stable Deployment Laws of Tether Release Rate and Tether Tension Control**

Entire TSS configurations require transformation according to different mission requirements. Changing the tether length is the most direct control method of system conception transformation, which includes tether release rate control and tether tension control. During TSS transformation, tether release rate and tether tension control parameters are controlled by a deployment device in the mother satellite.

#### *4.1. Tether Release Rate Control*

The tether release rate is directly controlled by a winch control mechanism in the mother satellite, and the influences of tether tension and *φ* are ignored. Hence, Equations (15)–(17) can be rewritten as follows:

$$\theta'' = 2(\theta' + 1)[\frac{\mathcal{e}\sin\nu}{\kappa} - \frac{l'}{l}] - \frac{\mathfrak{Z}}{2\kappa}\sin 2\theta \tag{26}$$

where *l* /*l* is the pseudo damping term, which makes *θ*, *φ* convergent.

#### 4.1.1. Fixed Angle *θ*

The system's in-plane angle is assumed to be a fixed angle of nonrotating motion. Corresponding to actual conditions, a Global Positioning System (GPS) rotates around the earth at a fixed angle in order to produce a stable state. *θ* = *θ*0, *θ* = *θ* = 0 are substituted into Equation (26), which can be written as follows.

$$\frac{l'(\nu)}{l(\nu)} = \frac{e\sin\nu}{1 + e\cos\nu} - \frac{3\sin 2\theta\_0}{4(1 + e\cos\nu)}\tag{27}$$

The log ratio of tether length, ln[*l*(*ν*)/*l*0(*ν*)], can be obtained by integrating Equation (27). Figure 6 shows that Equation (27) achieves a unified analytic solution when *θ*<sup>0</sup> = *kπ*/2 (*k* is an integer). Tether length is positively related to orbital eccentricity, and the abscissa corresponding to the highest point is *ν* = (2*k* + 1)*π*.

**Figure 6.** Log ratio of tether length versus *ν*, *θ* = *const*.

#### 4.1.2. Fixed Angular Velocity *θ*

The system is assumed to rotate steadily and uniformly in the direction of the in-plane angle, *θ* = *ων* + *θ*0, where *ω* is the angular velocity, and *θ*<sup>0</sup> is the initial value. *θ* = *ω* and *θ* = 0 are substituted into Equation (26), which can be written as follows.

$$\frac{l'(\nu)}{l(\nu)} = \frac{e\sin\nu}{1 + e\cos\nu} - \frac{3\sin[2(\omega\nu + \theta\_0)]}{4(\omega + 1)(1 + e\cos\nu)}\tag{28}$$

The log ratio of tether length can be obtained by applying integration.

As shown in Figure 7, tether length amplitude is positively related to *e*, and the system requires longer tethers to achieve stability control. The results above can be used to guide TSS system attitude control. When tether release conditions satisfy Equation (28), TSS systems can operate at a fixed angular velocity, *θ* .

**Figure 7.** Log ratio of tether length versus *ν*, *θ* = *const*: (**a**) *ω* = 1, *θ*<sup>0</sup> = 0 and *π*/2; (**b**) *ω* = 3/4, *θ*<sup>0</sup> = 0 and *π*/2.

#### *4.2. Tether Tension Control*

In order to enhance efficiency, applicability and stability, tether tension control is facilitated by TSS [34]. It is assumed that tether tension is the same at every point of the tether and is equal to the tension at the point of deployment/retrieval. The tether braking mechanism is modeled after the SEDS deployer, which uses a friction brake in order to control tether deployment speed [35]. Tether tension is expressed as follows:

$$T = \left[T\_0 + I\rho\dot{l}^2 \left(1 - A\_{\rm sol}l/L\_{\rm ref}\right)^{-E}\right] \exp(f\_\theta|\theta - \theta\_0| + 2\pi f\_\text{fl} n^\*)\tag{29}$$

where . *l* is the tether release rate. Tether tension control parameters are listed in Table 1.

A numerical simulation was performed in order to demonstrate control law performance, and simulation parameters are listed in Table 2, where *μ<sup>e</sup>* = 398, 600 km3/s2 is Earth's gravitational coefficient.



**Table 2.** TSS parameter values.


Figure 8 shows the dynamic response of the TSS deployment process. Figure 8a,b show variations of in-plane and out-of-plane pitch angles and roll angles versus the true anomaly, *ν*. It can be concluded that the system approaches the expected angle, 0 rad, after swinging under an initial perturbation. This result illustrates that a controlled deployment process is asymptotically stable and demonstrates the validity of the tether tension control equation (Equation (29)). Figure 8c shows tether length during deployment, which exhibits a smooth deployment curve, and the tether eventually reaches a stable length.

**Figure 8.** Deployment process of TSS: (**a**) tether in-plane angles versus *ν*; (**b**) tether out-plane vibration angles versus *ν*; (**c**) tether length versus *ν*.

**Remark 2.** *Analysis results of tether release rate control and tether tension control laws can provide effective feedback for TSS position and attitude.*

#### **5. Stability Analysis of TSS Deployment Using Floquet Theory**

Floquet theory is used to analyze the stability of solutions of linear ordinary differential equations with periodic variable coefficients. Local stability of deployment along preassigned pitch angles and roll angles can be analyzed by using Floquet theory. Tether length remains unchanged once the non-dimensional tether length equals one after accomplishing deployment in Section 4.2. Steps of Floquet theory applied to TSS systems are shown as follows.

In the case of *p*<sup>1</sup> = *θ*, *p*<sup>2</sup> = *θ* , *p*<sup>3</sup> = *φ*, *p*<sup>4</sup> = *φ* , the matrix form of Equations (15)–(17) can be summarized as follows:

$$\mathbf{P} = \begin{bmatrix} p'\_1 \\ p'\_2 \\ p'\_3 \\ p'\_4 \end{bmatrix} = \begin{bmatrix} \theta' \\ 2(\theta' + 1)(\frac{\epsilon \sin \upsilon}{\kappa} + \phi' \tan \phi) - \frac{3}{2\kappa} \sin(2\theta) \\ \phi' \\ \frac{2\epsilon \phi' \sin \upsilon}{\kappa} - \frac{1}{2} [(\theta' + 1)^2 + \frac{3}{\kappa} \cos^2 \theta] \sin(2\phi) \end{bmatrix} \tag{30}$$

where **p** = (*θ*, *θ* , *φ*, *φ* ) <sup>T</sup> are state-space vectors, and **<sup>p</sup>***<sup>s</sup>* = (*θs*, *<sup>θ</sup> <sup>s</sup>*, *φs*, *φ s*) <sup>T</sup> are equilibrium points, *θ <sup>s</sup>* = <sup>d</sup>*θ<sup>s</sup>* <sup>d</sup>*<sup>ν</sup>* . Equation (30) can be expressed as follows:

$$
\Phi' = \mathbf{A}(\mathbf{p}\_s)\Phi \tag{31}
$$

where **A**(**p***s*) is the Jacobian matrix of vector function, **P**, in a small neighborhood near the equilibrium point, **p***s*.

$$\mathbf{A(p\_s)} = \begin{bmatrix} 0 & 1 & 0 & 0 \\ -\frac{3}{k}\cos(2p\_{1s}) & 2(\frac{\varepsilon\sin\nu}{k} + p\_{4s}\tan p\_{3s}) & \frac{2p\_{4s}(p\_{2s}+1)}{\cos^2 p\_{3s}} & 2(p\_{2s}+1)\tan p\_{3s} \\ 0 & 0 & 0 & 1 \\ \frac{3}{2k}\sin(2p\_{1s})\sin(2p\_{3s}) & -(p\_{2s}+1)\sin(2p\_{3s}) & -[(p\_{2s}+1)^2 + \frac{3}{k}\cos^2 p\_{1s}]\cos(2p\_{3s}) & \frac{2\varepsilon\sin\nu}{k} \\ \end{bmatrix} \tag{32}$$

The period is 2*π*, and it can be expressed as follows.

$$\mathbf{A}(\nu, e) = \mathbf{A}(\nu + 2\pi, e) \tag{33}$$

The monodromy matrix can be obtained by integrating Equation (32) for one period from initial time *ν* = 0, which combines with the initial condition **Φ**(*ν*0, *e*) = **I**, where **I**4×<sup>4</sup> is the identity matrix.

$$\mathbf{M} = \Phi(2\pi, e) = e^{\int\_0^{2\pi} \mathbf{A}(\nu) d\nu} \tag{34}$$

According to Floquet theory, the stability of the zero solution of Equations (15)–(17) can be assessed by a Floquet multiplier:

$$\begin{cases} |\lambda\_i|\_{\max} < 1, \ (i = 1, \ 2, \ 3, \ 4) & \text{Aymptotially stable} \\ \quad |\lambda\_i|\_{\max} = 1, \ (i = 1, \ 2, \ 3, \ 4) & \text{Underermined} \\ \quad |\lambda\_i|\_{\max} > 1, \ (i = 1, \ 2, \ 3, \ 4) & \text{Unstable} \end{cases} \tag{35}$$

where |*λi*|max is a Floquet multiplier that is the maximum of the absolute value of *λi*, and *λ* is an eigenvalue of the monodromy matrix, **M**, in Equation (34).

Figure 9a shows the relationship between Floquet multipliers and stability of expected inplane angles for *e* = 0, where the case of *φ<sup>s</sup>* = 0 is discussed. As shown in Figure 9b, this result shows that Floquet multipliers are less than one for *θ<sup>s</sup>* ∈ (−1.584, −1.563) and (1.563, 1.584) by symmetry, which illustrates that deployment in a short interval of *θ<sup>s</sup>* is asymptotically stable. However, once *θ<sup>s</sup>* lies outside the specified range, Floquet multipliers are always greater than one, which shows that the deployment process is undetermined since the

expected in-plane angle lies outside the domain of stability. A similar result is achieved with *e* ∈ (0, 0.5).

**Figure 9.** (**a**) Floquet multipliers versus *θs*, *e* = 0; (**b**) local enlargement of Floquet multipliers versus *θs*.

#### **6. Conclusions**

In this paper, nonlinear dynamic characteristics of TSS during a configuration conversion process were analyzed based on a simplified rigid-rod model. Tether tension control was proposed, and numerical simulations show that the proposed law can suppress in-plane and out-of-plane librations of rigid tethered satellites during deployment, and spacecraft and tether stability control goals can be achieved. The periodic stability of time-varying control systems was analyzed by using Floquet theory, and small parameter regions of TSS in asymptotically stable states were expressed.

In summary, this paper has provided tether release rate and tether tension control laws for suppressing wide-ranging TSS vibrations that are valuable for improving TSS attitude control accuracy and performance, specifically for TSSs that are operating in loweccentricity orbits. Additionally, future studies based on existing research can be conducted with respect to two aspects: (1) A more accurate model can be established since, in the current study, the tether was discretized into a series of lumped masses connected by springs and dampers with mass. (2) Applications of accurate models can generate more dimensions, and Floquet theory used to analyze the stability of high-dimensional dynamic systems requires further verification.

**Author Contributions:** Conceptualization, K.Z. and K.L.; methodology, K.Z.; software, K.Z.; validation, K.Z., K.L. and X.G.; formal analysis, K.Z., K.L., C.F.; investigation, K.Z., K.L.; resources, K.L.; data curation, K.L.; writing—original draft preparation, K.Z.; writing—review and editing, C.F., S.Z.; visualization, K.L.; supervision, C.F., S.Z.; project administration, X.G.; funding acquisition, K.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Key R&D program of China (Grant No. 2020YFA0711700); the National Natural Science Foundation of China (Grant Nos. 12072263, 11802235); Natural Science Foundation of Shaanxi Province (Grant No. 2020JQ-129); State Key Laboratory of Mechanical Behavior and System Safety of Traffic Engineering Structures (Grant No. KF2020-26); and China Postdoctoral Science Foundation (2021T140033, 2021M690274).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** All necessary and relevant data are included in this paper.

**Acknowledgments:** The authors thank the editor and reviewers for their valuable suggestions.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Development of a Wireless Corrosion Detection System for Steel-Framed Structures Using Pulsed Eddy Currents**

**Namhoon Ha 1, Han-Seung Lee <sup>2</sup> and Songjun Lee 1,\***


**Abstract:** Structural health monitoring (SHM) can be more efficient with the application of a wireless sensor network (WSN). However, the hardware that makes up this system should have sufficient performance to sample the data collected from the sensor in real-time situations. High-performance hardware can be used for this purpose, but is not suitable in this application because of its relatively high power consumption, high cost, large size, and so on. In this paper, an optimal remote monitoring system platform for SHM is proposed based on pulsed eddy current (PEC) that is utilized for measuring the corrosion of a steel-framed construction. A circuit to delay the PEC response based on the resistance–inductance–capacitance (RLC) combination was designed for data sampling to utilize the conventional hardware of WSN for SHM, and this approach was verified by simulations and experiments. Especially, the importance of configuring sensing modules and the WSN for remote monitoring were studied, and the PEC responses caused by the corrosion of a specimen made with steel were able to be sampled remotely using the proposed system. Therefore, we present a remote SHM system platform for diagnosing the corrosion condition of a building with a steel structure, and proving its viability with experiments.

**Keywords:** structural health monitoring; wireless sensor network; steel-framed construction; corrosion; pulsed eddy current

#### **1. Introduction**

SHM, which evaluates the durability of building structures, diagnoses points with damage and finds their location by collecting data using a sensor system in real time [1]. Additionally, advanced methods were introduced to reconstruct the lost data for the precise SHM. [2,3]. However, the conventional wired system used in SHM is uneconomical because a large amount of wire is necessary and requires substantial labor during the installation and maintenance periods. If a WSN is applied to the SHM system, building structures can be conveniently maintained with a low cost [4,5]. When considering the total economic costs of WSN SHM, the operating time should be taken into account, because the WSN is generally powered by a battery.

SHM includes measuring temperature, humidity, wind speed, earthquake incidence, and corrosion. Corrosion, which affects the durability of buildings, occurs in all steel structural materials. Various steel sections, such as wide flanges, I-beams, and channels, are used for structures such as buildings, roads, and bridges, which can be easily seen around us, and these structures are classified by the following construction method. Steelframed construction (SC) [6] is constructed quickly and has low local environmental pollution and can be applied for the construction of building parking lots, but, in this case, severe corrosion can occur due to direct exposure to weather conditions [7]. Steel-framed reinforced concrete (SRC) has a light self-weight, high strength, and high stiffness by combining steel beams with reinforced concrete (RC) [8] and is used in large construction projects, such as skyscrapers, bridges, and tunnels. However, bridges are easily corroded

**Citation:** Ha, N.; Lee, H.-S.; Lee, S. Development of a Wireless Corrosion Detection System for Steel-Framed Structures Using Pulsed Eddy Currents. *Sensors* **2021**, *21*, 8199. https://doi.org/10.3390/s21248199

Academic Editors: Yongbo Li, Bing Li, Jinchen Ji and Hamed Kalhori

Received: 4 November 2021 Accepted: 6 December 2021 Published: 8 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

by ionic deicing chemicals used in winter [9], and subsea tunnels are also corroded due to electrochemical action caused by chloride ion invasion [10]. The corrosion of steel-framed construction, which is caused by various factors in the diverse environments, should be measured to provide a warning before breakdown, because it can cause cracks, the spalling of concrete, and structural collapse [11].

There are several technologies used to measure corrosion, including the following methods: eddy current [12,13], ground penetrating radar (GPR) [14], galvanostatic pulse method (GPM) [15], fiber Bragg grating (FBG) [16], ultrasonic pulse velocity (UPV) [17], and infrared thermography (IRT) [18]. Eddy current testing (ECT) is the method used in this research. The mechanism of the eddy current method is to measure the conductivity and permeability changes caused by corrosion by inducing an eddy current that is generated by a sine wave or pulse [19–23]. Although the method using a single-wavelength sine wave successfully measures corrosion, the detectable depth is limited by the skin effect. Since corrosion occurs not only on the surface of a steel frame but also inside of it, it is appropriate to use multiple-frequency waves as input signals to evaluate the durability of structures. However, the method of generating multiple frequencies requires much more complex and expensive electronics than a single-frequency system, involving generating input signals and measuring output signals, and it is not suitable for minimizing cost or power consumption. On the other hand, the method of PEC using a pulse signal as an input has the possibility to minimize the power consumption [24,25]. Further, since it covers multiple frequencies which can detect various depths of objects without actually changing the frequency [12,26], PEC can be used as an alternative to explicit multiplefrequency methods for inspecting corrosion deep inside a steel frame. Thus, to detect corrosion, the pulse input is applied for the steel structure, and the output from the PEC response should be measured. A typical PEC response appears as an exponential decay for several milliseconds [27], shown in Figure 1. Data acquisition (DAQ) equipment, such as an analog-to-digital converter (ADC) board, is used to measure the response [28,29], but this configuration is unsuitable for application in an actual situation that requires a low-power system, such as the typical WSN environment. For a WSN, the low cost, small size, and low power consumption are usually taken into account. Therefore, measuring the PEC response is difficult when applying a conventional system for use in an actual SHM application.

**Figure 1.** A typical PEC response.

In this paper, we propose a method for detecting the corrosion of a steel-framed construction with a convenient monitoring system using WSN. A circuit designed to delay the PEC response makes it possible to easily deploy in an actual construction environment. After a PEC is induced to detect corrosion, in order to measure the response, a delay circuit for the response signal should be provided, because the output signals vary too fast for sampling them using a conventional measurement system. By applying the proposed method using the WSN, the remote monitoring system can be utilized for a more convenient real-time analysis of the corrosion state. In order to achieve this, tiny sensor modules, without large-sized, general-purpose measuring equipment, should be developed and installed in various locations of the building for a more efficient and precise SHM. Additionally, the more collected data from the several parts of a building, the more efficient SHM is. Additionally, various analysis environments with a convenient user interface are provided in the hardware configurations.

#### **2. Pulsed Eddy Current Response**

PEC response must be measured to evaluate the amount of corrosion, but it is difficult to measure this with the hardware configuration generally used for WSNs. (Figure 1) In this section, a delay circuit designed to sample PEC responses, in order to compare the results from the hardware configuration of a conventional WSN system, is described.

Figure 2 shows the configuration of the proposed delay circuit of the sensor module. It is based on an RLC circuit and constitutes a loop circuit by connecting all components in series. *C*<sup>1</sup> and *C*<sup>2</sup> are capacitors, and a sensor coil is placed between them. The sensor coil has an inductance (*Ls*) and resistance (*Rs*). *VPulse* is the input source, and *ILoop* is the current flowing in the loop circuit. According to Kirchhoff's voltage law (KVL), the circuit is described as follows:

**Figure 2.** The delay circuit of a sensor module for sampling PEC response.

However, Equation (1) is based on the circuit in Figure 2 without any specimen. Therefore, when measuring corrosion, the effective resistance of the coil sensor (Δ*RV*) is varied by the specimen. It is defined as follows:

$$R\_{VS} = R\_S + \Delta R\_V \tag{2}$$

The effective inductance of the sensor (Δ*LV*) coil is also altered by the specimen, which could be a conductive material, and is defined as follows:

$$L\_{VS} = L\_S + \Delta L\_V \tag{3}$$

The substitution of Equations (2) and (3) into (1) is expressed as follows:

$$V\_{\rm Pulse} + \frac{\dot{j}}{\omega \mathbf{C}\_1} I\_{\rmLoop} - (R\_{VS} + j\omega L\_{VS}) I\_{\rmLoop} + \frac{\dot{j}}{\omega \mathbf{C}\_2} I\_{\rmLoop} = 0 \tag{4}$$

*ILoop* is calculated as Equation (5):

$$I\_{Loop} = -\frac{V\_{Pulse}}{j\left(\frac{1}{\omega \cdot \mathbf{C}\_1} + \frac{1}{\omega \cdot \mathbf{C}\_2} - \omega L\_{VS}\right) - (R\_{VS})}\tag{5}$$

The voltage of *C*<sup>2</sup> is expressed as the following Equation (6):

$$V\_{\mathbb{C2}} = -j \frac{I\_{Loop}}{\omega \mathbb{C2}} \tag{6}$$

Using Equations (5) and (6) can be described as follows:

$$V\_{C2} = j \frac{V\_{Pulse}}{\omega \mathcal{C}\_2 \left[ j \left( \frac{1}{\omega \mathcal{C}\_1} + \frac{1}{\omega \mathcal{C}\_2} - \omega L\_{VS} \right) - \left( R\_{VS} \right) \right]} \tag{7}$$

Therefore, it can be seen that a change in the specimen's electrical properties produces a variation in the PEC.

In order to operate the delay circuit for a simulation, 4.7 nF ceramic capacitors are selected for both *C*<sup>1</sup> and *C*2, and the inductance and resistance of the sensor coil are set to 195 uH and 2.5 Ω, respectively. The *VPulse* square wave with 1 Hz, is an input to the RLC circuit. Figure 3 shows the result achieved by an electronic circuit simulator, EveryCircuit, confirming that the PEC response is delayed by the proposed circuit. The exponentially decaying signal obtained by simulation appeared for about 500 ms, and demonstrated that it took approximately 100 times longer than the general PEC response shown in Figure 1.

**Figure 3.** A simulation of the delay circuit for monitoring the PEC response.

#### **3. Corrosion Remote Monitoring System**

In this section, we propose a system for implementing the design described in the previous section, which is able to sample the delayed PEC response. The sensor node is made by combining the proposed circuit with the hardware configuration of a conventional Zigbee-based WSN for remote monitoring. Additionally, a user-friendly interface is implemented for convenience.

#### *3.1. Sensor Node Design*

For the detection and remote monitoring of corrosion, the proposed system requires several functions, including data collection and communication. Figure 4 shows the hardware configuration of the sensor node. An Arduino-based system, which has many commercial modules for expansion, sufficient open library codes, and a high compatibility, provides a complimentary integrated development environment (IDE) for developers. Thus, the Arduino Pro Mini was chosen to develop the hardware for the proposed system. In addition, it is utilized as a pulse generator, because a digital input signal can be generated by controlling its GPIO. On the other hand, in the case of a conventional ECT, a magnetic sensor, such as a Hall sensor, is used to detect the intensity of the eddy current, but, in this research, the only variation of the voltage of the coil caused by the PEC is measured without any additional sensor, while the PEC is induced on a specimen through the same sensor coil. If the number of turns of the coil for the sensor is increased, the intensity of PEC response becomes high. Thus, the sensor coil with a larger size has a higher sensitivity. In order to determine the size of the sensor, it is designed according to the size and shape of the part that is measured, optimized by a simulation or experiment, and then applied to the actual building structures. In this research, the sensor coil has a planar square shape. It can be placed close to a steel component, and is suitable for measuring the corrosion of a large area. In order to make the planar coil sensor, a foam board, double-sided adhesive tape, and AWG26 wire were used, and its size is arbitrarily chosen to be 150 mm by 150 mm for an experiment. For sampling the PEC response signals detected by the sensor coil, the ADS1015 analog-to-digital converter is used, which has a 12-bit resolution and a sampling rate of 3300 samples per second. The sampled data are stored in an SD card for every sampling as backup data, which can prevent data loss caused by communication errors. In this research, a sampling rate of 500 samples per second is chosen for the experiment, and a smaller sampling rate is possible to detect corrosion for SHM. The proposed system configures the WSN to remotely transmit the measured data to the master node of a Zigbee communication device (xBee s2c, DIGI, Hopkins, MN, USA). The Zigbee module based on the IEEE 802.15.4 standard has various advantages such as a low power consumption, low cost, and high compatibility with various network topologies. Therefore, it is suitable for an actual application of a WSN for structural health monitoring.

**Figure 4.** Hardware configuration of the sensor node.

#### *3.2. Networking and Monitoring*

Figure 5 shows the architecture of the system for remotely monitoring the corrosion data collected from the sensor node. The master node in the gateway layer is implemented using a Raspberry Pi 3 B+, which is a single-board computer the size of a credit card and has a built-in operating system. It is used as both a master node and server. In order to use a relational database management system, MariaDB is installed to collect data from the sensor node and for storage in the database. The data are visualized for analysis, and Grafana, a web-based interactive application, provides convenient features to visualize the data from a database. Additionally, the Raspberry Pi has both an ethernet port and a wireless LAN, and it can be connected to the Internet without an additional device. As a result, users can access the master node through the Internet and monitor the data measured from the construction site anytime and anywhere.

**Figure 5.** The proposed architecture of the monitoring system.

#### **4. Experimental Section**

Experiments were conducted to confirm the delay in the PEC response and verify the detectability of corrosion of the steel plate. Figure 6 shows the experimental setup to examine corrosion using the proposed system. The lift-off effect is caused by variations in the distance between the sensor and the specimen. The thickness of the rust on the specimen is regarded as the variation of distance, and it becomes a factor that interferes with the signal that is collected. Several studies aimed to reduce this effect [30–32]. However, when the sensor coil is buried or fixed, the effect is generated by the corrosion of the specimen, and the same principle of a thickness measurement is applied to the evaluation of corrosion [13]. Accordingly, in the experimental setup, the thickness reduction due to corrosion is replaced by the change in the thickness of the steel plates, and the consequent

lift-off effect was simulated by placing non-conductive acrylic plates between the steel plate and the coil. The size of the steel plates and the acrylic plates used in the experiment were each 150 mm by 150 mm, and had the same dimensions as the sensor coil with thicknesses of 2 mm, 4 mm, 6 mm, 8 mm, and 10 mm. The collected experimental data were recorded and monitored by the proposed system. This experiment measured the change in the response of PEC according to the thickness of corrosion. When installed in an actual building, the sensor coil should not prevent the corrosion of the steel frame, so it should be installed at a suitable distance.

**Figure 6.** The experimental setup.

#### **5. Result**

In order to verify the delay of the PEC response using the sensor circuit explained in Section 2, the signal was measured using an oscilloscope. Figure 7 shows the output signal when the pulse is applied to the sensor circuit. The signal induced by the proposed sensor circuit showed a peak value of approximately 2.5 V and an attenuation of a signal of about 500 ms, which were very similar to the simulation result in Figure 3. It shows that the sampling is possible by the ADC, which is generally used in the actual situation.

Figure 8 shows the experimental result of sampling the delayed PEC response using the proposed system and examines if corrosion is detected from the steel-framed construction sample. Each line in Figure 8 indicates the voltage variation of *VC*<sup>2</sup> according to the level of corrosion, which represents the PEC response that is sampled 500 times for about one second by the sensor node. The thicknesses of the steel plates (S) and the acrylic plates (A) used to simulate the depth of corrosion were 10 mm for the steel plate and 0 mm for the acrylic plate for Line 1; 8 mm for the steel plate and 2 mm for the acrylic plate for Line 2; 6 mm for the steel plate and 4 mm for the acrylic plate for Line 3; 4 mm for the steel plate and 6 mm for the acrylic plate for Line 4; and 2 mm for the steel plate and 8 mm for the acrylic plate for Line 5, respectively. According to the results, the delayed PEC response was successfully sampled by the sensor node with the proposed circuit, and the

corrosion was identified by considering the fact that the peak value of *VC*<sup>2</sup> consistently decreased as the corrosion progressed. The significant thickness of corrosion depends on the kinds of building structures and their environments [33]. The amount of corrosion can be evaluated without a high-end data acquisition device nor a signal conditioner with a computer system, and the cost of the sensor module can be reduced because the coil is utilized for both inducing eddy current and detecting the effect of the eddy current without an additional device. Thus, the developed low-cost wireless SHM system can be utilized to conveniently measure the corrosion of steel-framed construction with convenience at an actual construction site. In addition, for the precise SHM, the simulation or experiment for the specific condition should be carried out before installing the system for an actual building structure. In case of the long-term stability of the system, the accumulated measurement error of sensor inside the building structure, due to various causes, should be corrected by periodic calibration so that it is not regarded as corrosion.

**Figure 7.** The delayed PEC response measured by the sensor circuit.

**Figure 8.** The variation of the PEC response by corrosion level.

Additionally, the variations of the inductance and the resistance of the sensor coil were examined using an LCR meter to verify the relationship with the corrosion level. (Figure 9) The red line indicates the inductance of the sensor coil, and the blue line indicates the resistance of sensor coil. As a result, the inductance and resistance varied by up to 26%, and 35%, respectively, due to the corrosion level. The PEC induced in the steel plate generates a magnetic field whose direction is opposite to that of the sensor coil. The inductance of the sensor coil decreases due to the reduction in magnetic flux in the coil. Additionally, the resistance of the sensor coil increases because of the energy dissipation caused by the PEC [34]. This was verified by experiments and analyses in the study.

**Figure 9.** Resistance and inductance of the sensor coil by corrosion level.

Figure 10 shows the detectable range of the sensor, which is experimental result measured using a steel plate and acrylic plates (0 mm, 25 mm, 50 mm). According to the blue line (S10 A50), the maximum distance for the detection of steel plate is about 50 mm because it is almost similar to the pink line (S10 A50+), which shows the voltage values in the case that the steel plate is placed further than 50mm from the sensor coil.

**Figure 10.** Detectable range of the proposed sensor.

In the experiments, a viability of applying the proposed method to a steel-framed structure was confirmed by the fact that the corrosion level was detected by variation of PEC. According to the previous work [13], detecting the corrosion of building structures is possible using PEC methodology. Tiny sensor modules developed by the proposed method can be installed in various parts of the real building structures more easily than conventional hardware equipment for SHM.

Figure 11 is the screen that shows the results, according to corrosion status, using the Grafana web-based monitoring application, based on the hardware configuration explained in Section 3. Grafana provides various visualization tools, such as graphs, tables, and bar charts, which use data from the sensor nodes. Moreover, the server computer, implemented with a Raspberry Pi, is suitable for the proposed system because of its low cost and tiny size. Therefore, users can easily access the monitoring system without the restrictions of time and space, provided that they have a smartphone or tablet PC connected to the Internet.

**Figure 11.** The dashboard for monitoring corrosion.

#### **6. Conclusions**

In this paper, we proposed a platform that could perform the remote monitoring of the corrosion of steel-framed construction in real time using the PEC method. For the SHM of the steel-framed construction, the PEC response was sampled with the designed hardware configuration containing a delay circuit, and the performance of the sensor was confirmed by measuring the corrosion of the test sample in experiments. Furthermore, by applying WSN and IoT technology, a real-time remote monitoring system was implemented that was easily accessible for user and could efficiently analyze the status of corrosion with database and visualization software. By using the proposed method, the tiny, low-cost hardware module for SHM can be manufactured without a function generator, oscilloscope, and other general-purpose measuring equipment. Additionally, real-time and remote monitoring become possible if they are applied to many locations in the building structure with wireless networks. Thus, the efficiency becomes higher than the conventional SHM method. Moreover, the system also provides convenience, efficiency, and portability by using a Raspberry Pi for the server computer.

Furthermore, the proposed system can be applied in other fields, such as non-destructive testing related to cracks or the thickness of paint for aircraft and ships by designing sensor coils with various shapes depending on the particular desired application.

The limitation of the proposed system is the power consumption for the long-term measurements. The power consumption of sensor modules can be reduced by data compression and more efficient communication algorithms, etc. Additionally, the operation time can be improved by applying energy harvesting techniques such as wind and solar power, etc.

**Author Contributions:** Conceptualization, N.H., H.-S.L. and S.L.; methodology, N.H., H.-S.L. and S.L.; software, N.H.; validation, N.H.; writing—original draft preparation, N.H.; writing—review and editing, N.H. and S.L.; supervision, S.L.; funding acquisition, H.-S.L. and S.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by a National Research Foundation of Korea (NRF) grant funded by the Korean government (MSIT) (No. 2020R1F1A107534412 and 2015R1A5A1037548).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Sensors* Editorial Office E-mail: sensors@mdpi.com www.mdpi.com/journal/sensors

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34

www.mdpi.com

ISBN 978-3-0365-6182-0