**An Intelligent Fault Diagnosis Method for Bearings with Variable Rotating Speed Based on Pythagorean Spatial Pyramid Pooling CNN**

### **Sheng Guo , Tao Yang \* , Wei Gao, Chen Zhang and Yanping Zhang**

School of Energy and Power Engineering, Huazhong University of Science & Technology, Wuhan 430074, China; levykwok@hust.edu.cn (S.G.); gw@hust.edu.cn (W.G.); zhangchen710@yeah.net (C.Z.); zyp2817@hust.edu.cn (Y.Z.)

**\*** Correspondence: hust\_yt@hust.edu.cn; Tel.: +86-027-87542817

Received: 16 October 2018; Accepted: 7 November 2018; Published: 9 November 2018

**Abstract:** Deep learning methods have been introduced for fault diagnosis of rotating machinery. Most methods have good performance when processing bearing data at a certain rotating speed. However, most rotating machinery in industrial practice has variable working speed. When processing the bearing data with variable rotating speed, the existing methods have low accuracies, or need complex parameter adjustments. To solve this problem, a fault diagnosis method based on continuous wavelet transform scalogram (CWTS) and Pythagorean spatial pyramid pooling convolutional neural network (PSPP-CNN) is proposed in this paper. In this method, continuous wavelet transform is used to decompose vibration signals into CWTSs with different scale ranges according to the rotating speed. By adding a PSPP layer, CNN can process CWTSs in different sizes. Then the fault diagnosis of variable rotating speed bearing can be carried out by a single CNN model without complex parameter adjustment. Compared with a spatial pyramid pooling (SPP) layer that has been used in CNN, a PSPP layer locates as front layer of CNN. Thus, the features obtained by PSPP layer can be delivered to convolutional layers for further feature extraction. According to experiment results, this method has higher diagnosis accuracy for variable rotating speed bearing than other methods. In addition, the PSPP-CNN model trained by data at some rotating speeds can be used to diagnose bearing fault at full working speed.

**Keywords:** convolutional neural network; spatial pyramid pooling; fault diagnosis; bearing; wavelet transform

### **1. Introduction**

As most rotating machinery are the key equipment in production and work in a high-speed rotating environment, their failures will cause major economic losses and safety accidents. Therefore, it is important to detect equipment failure as soon as possible. An intelligent fault diagnosis method has the motivation that it can detect the fault of operating machinery in real time without manual operation. With the popularization of online vibration monitoring systems, manufacturers have accumulated a large amount of data that can support intelligent fault diagnosis methods. Many intelligent fault diagnosis methods have been proposed to diagnose bearing faults [1–5]. However, most of them use a simple classifier and focus on fault feature extraction algorithms. When these methods are used to diagnose bearings with complex operating conditions, the simple classifier cannot process large amounts of monitoring data and can easily cause over fitting.

As the advanced representation of intelligent algorithms, deep learning methods have greatly changed our daily life. They are successfully applied in many different areas such as computer vision, object detection, natural language processing, and even disease diagnosis [6–10]. Deep learning

methods can recognize high-dimensional complex input, get rid of the reliance on signal processing techniques and hand-engineered feature extraction algorithms. Unlike most previous artificial intelligent methods that can only process one-dimensional hand-engineered features [11], they can process the two-dimensional results of some basic signal processing methods directly. With these advantages, deep learning methods have been introduced to the fault diagnosis of rotating machinery, such as convolutional neural network (CNN), deep belief nets (DBN), and recurrent neural networks (RNN) [12–17]. Most methods perform well when dealing with bearing data at a certain rotating speed. However, in practical use, most machinery, such as wind turbines, pumps, and fans, have varying rotating speed. The performance of the methods has not been verified when processing data with varying rotating speed. Liu et al. [18] propose a dislocated time series CNN method for bearing diagnosis and apply it to varying rotating speed data by testing different network parameters to achieve the best result. However, the number of parameters that can be tested is limited, and the method is difficult to realize in practical applications. Meanwhile, in practical use, faults may not happen fully in all working speeds. Then there will not be enough fault data that covers full working speeds. Therefore, there needs to be a new deep learning fault diagnosis framework which can deal with data with varying rotating speed without complex parameters attempt. In addition, it needs to be able to realize diagnosis in all working speed based on limited fault data.

The continuous wavelet transform (CWT) has been proved to be a useful method to analyze vibration signals [19–21]. As a time-frequency analysis method, the result contains the complete time-frequency domain information of the vibration signal and avoids information loss of the original signal. In addition, it is suitable for detecting bearing faults which is usually presented as shock signals. When used in fault diagnosis of bearings, it has advantages compared with some other time-frequency methods, such as Short-time Fourier Transform (STFT), Discrete Wavelet Transform (DWT), Wavelet Package Decomposition (WPD) and Empirical Mode Decomposition (EMD).

Short-time Fourier Transform is the time-frequency transform based on Fourier transform. Because the window size is fixed, it only applies to stationary signals with small frequency fluctuations. In addition, the result is susceptible to noise interference. DWT is the discretization of the scale and displacement of CWT according to the power of 2. It retains less time-frequency information than CWT and may lose the key information near fault characteristic frequencies of bearing. WPD is an improved method of CWT. It provides a more detailed decomposition of high-frequency components. The fault characteristic frequencies of a bearing are usually less than 12 multiples of the rotating frequency. With the high sampling frequency, all the information near fault characteristic frequencies can be got by CWT. There is no need to use WPD, which is a more time-consuming method than CWT. Empirical Mode Decomposition decomposes the signal based on the time scale characteristics of the signal itself. Without a preset basis function, the location of the fault feature is uncertain in EMD when it applies to signals of different sensors. Therefore, when used for deep learning method, CWT, which uses a fixed wavelet basis function to decompose all the signals, is a better choice.

However, in most cases, a feature exaction method is used to exact a one-dimensional vector from the two-dimensional continuous wavelet transform coefficients. It may result in the loss of key fault information. Continuous wavelet transform scalogram (CWTS) contains all the continuous wavelet transform coefficients. It is a two-dimensional matrix which contains the complete time-frequency domain information of the vibration signal. With its powerful image recognition ability, CNN is the most suitable deep learning method to deal with CWTS. When applied to data with varying rotating speed, the CWTSs will have different size if they have the same frequency multiplication range of the rotating frequency. Without a cropping or warping operation, ordinary CNN can only process the input in the same size. Therefore, a CNN with new structure is needed to process CWTSs in different sizes.

To overcome the problems and challenges above, this paper proposes a fault diagnosis method based on continuous wavelet transform and Pythagorean spatial pyramid pooling (PSPP) CNN. This method uses a continuous wavelet transform to decompose vibration signals into CWTS in

different scales according to the rotating speed. Using the PSPP strategy, CNN could then process the different size scalograms. Therefore, the fault diagnosis of data at variable rotating speed can be carried out by a single CNN model. The PSPP strategy is an improvement on spatial pyramid pooling (SPP). A PSPP layer can locate as front layer of CNN. Thus, the features obtained by the PSPP layer can be delivered to the convolutional layers for further feature extraction. Experiments are carried out on data from two different testbeds, constant rotating speed data and variable rotating speed data, respectively. The results demonstrate the effectiveness of the proposed approach. The contributions of the proposed approach are as follows:


The paper is organized as follows. Section 2 presents the fault diagnosis method that combines CWTS and PSPP-CNN for fault diagnosis, with a detailed procedure of the proposed method and SPP, the proposed PSPP layer, and the structure of PSPP-CNN used in this paper. Experimental verification of the method which includes constant rotating speed data and variable rotating speed data is described in Section 3. Finally, the concluding remarks are given in Section 4.

### **2. Proposed Method**

As described above, CNN has been successfully applied to fault diagnosis. However, the proposed fault diagnosis models lack the diagnosis of variable rotating speed data, which has practical engineering value for online diagnosis of variable speed equipment such as the wind turbine. Therefore, this paper proposes a fault diagnosis framework based on CWTS and PSPP-CNN. Figure 1 illustrates the procedure of the proposed method. First, accelerators are used to collect the vibration signals of bearing. Second, continuous wavelet transform is used to decompose vibration signals into CWTSs. Next, as fault characteristic frequencies of bearings are related to rotating speed, the CWTSs are cropped into different sizes according to rotating speed. Then using PSPP strategy, CNN can deal with the input of different sizes. Therefore, the CWTSs in different sizes can be trained using a single PSPP-CNN. Finally, the test signals of bearing with variable rotating speed need to be decomposed into CWTSs with different scale ranges according to rotating speed. Using the CWTSs as the input of the trained PSPP-CNN, fault diagnosis of the signals can be achieved. Details of the main steps of the proposed method are described as follows:

**Figure 1.** Flow chart of proposed fault diagnosis method. CWTS represents Continuous Wavelet Transform Scalogram, and PSPP-CNN represents Pythagorean Spatial Pyramid Pooling Convolutional Neural Network.

### *2.1. Continuous Wavelet Transform Scalogram*

The continuous wavelet transform decomposes a signal in the time-frequency domain by using a family of wavelet functions to obtain feature values. Next, by analyzing the continuous wavelet coefficients or using the classification algorithm, we gain insight about the fault condition of the equipment. The process of continuous wavelet transformation can be described as:

$$\Psi\_{a,b}(t) = |a|^{-\frac{1}{2}} \Psi\left(\frac{t-b}{a}\right) \text{ a, } b \in \mathbb{R} \text{ a } \neq 0 \tag{1}$$

$$\mathbb{C}\_{a}(k) = \int \mathbf{x}(t) \overline{\mathbf{P}}\_{a,b}(t) dt \tag{2}$$

where Ψ*a,b*(*t*) is a wavelet function whose shape and displacement are determined by *a*, the scale parameter, and *b*, the translation parameter. *x* is a signal with *m* data points. The wavelet coefficient of *x*(*t*) at the *a*-th scale is *Ca* (*a* = 1,2,3,··· ,*l*). *k* is time order (*k* = 1,2,3,··· ,*m*). Ψ*a*,*b*(*t*) is the complex conjugate of the wavelet function at scale *a* and translation *b*.

To show the change of wavelet coefficients intuitively, a CWTS is proposed [22]. The CWTS expresses continuous wavelet coefficients by a two-dimensional image in the time-frequency domain. Put all wavelet coefficients in a matrix *<sup>P</sup>* = [*C*1,*C*2,··· ,*Cl*]. The graph of wavelet coefficients matrix *<sup>P</sup>* is called a CWTS.

Figure 2 shows the CWTS of a ball fault bearing signal with a 2400 rpm rotating speed sampled at 12 kHz. It is decomposed by the Morlet wavelet from 1 to 300-scale and has 300 data points in time series. The horizontal axis represents the position along the time direction, and the vertical axis represents the scale. Morlet wavelet is chosen as the wavelet used in this paper. Because it has the similar shape with the shock signal caused by bearing faults [23]. In addition, the signal extracted by the Morlet wavelet has the higher energy-to-Shannon entropy ratio than the other common wavelet types. Energy-to-Shannon entropy ratio is an important indicator to measure the fitness of wavelet functions [24].

**Figure 2.** Continuous Wavelet Transform Scalogram (CWTS) of a ball fault bearing signal. The darker pixels correspond to larger wavelet coefficients.

### *2.2. Continuous Wavelet Transform Scalogram Cropping*

The object recognition using CNN is concerned with the shape of the object. If the shape appears in the image, the existence of the object is detected. However, in signal processing area, the axes of images constructed usually have clearly defined meanings. The appearance of the same shape in different location may indicate different fault modes. Therefore, the location of the features should also be paid attention to.

As the vertical axis of CWTS represents the scale, different positions on the vertical axis relate to different frequencies. As we know, the fault characteristic frequencies of bearings are related to rotating speed. If the input CWTSs of PSPP-CNN can be ensued to have the same time domain range and frequency multiplication of the rotating frequency, the fault characteristic of the same fault could appear at the similar position in CWTSs. Thus, the classification of CWTSs will achieve a comparatively accurate result. Therefore, a CWTS cropping step is proposed in this paper.

Suppose a vibration signal *x*(*i*) (*i* = 1,2, ... *m*) is collected at a sampling frequency *f* (Hz) with *m* sampling data points. The rotating speed is *n* (rpm), corresponding to a machine rotating frequency *fm* = *n*/60. The integer multiple *f* to *fm* is

$$q = \frac{f}{f\_m} + \frac{1}{2} = \frac{60f}{n} + \frac{1}{2} \tag{3}$$

From Equation (1), we can calculate that the central frequency of a wavelet function is inversely proportional to scale *a*. Suppose *fj* = *k*/*j*, where *fj* is the central frequency at scale *j*, and *k* is the proportionality coefficient. According to the Morlet wavelet function,

$$k = f\_0 \times f \tag{4}$$

where *f* <sup>0</sup> is the center frequency of the wavelet function, and the range is from 0.796 to 0.955. In this paper, 0.955 is chosen as *f* 0. Therefore, we choose scale 1 as the starting scale of cropping which corresponds to a high frequency. To make the end scale the same multiple of the rotating frequency, we choose *q* as the end scale. Scale *q* corresponds to the frequency:

$$f\_{\emptyset} = k/q = k/f \times f\_{m} = f\_{0} \times f\_{m} \tag{5}$$

Thus, *fq* is the *f* <sup>0</sup> multiplication of the rotating frequency. As the fault characteristic frequencies of bearings are larger than the rotating frequency, the cropping from scale 1 to *q* at the scale axis is sufficient for bearing fault diagnosis. The cropped CWTS will contain all the fault characteristic frequencies of bearings needed for analysis.

For the time domain axis, the time of *q* length data points is *t* = *q*/*f* ≈ 60/*n*, which is about the time duration of a rotor rotation cycle.

Therefore, in this paper, the original CWTS is cropped from scale 1 to *q*, and *q* length in the time domain axis. Thus, the cropped CWTSs with different rotating speeds have the same time domain range and frequency multiplication relative to the rotating speed.

### *2.3. Pythagorean Spatial Pyramid Pooling Convolutional Neural Network Training*

### 2.3.1. Pythagorean Spatial Pyramid Pooling Convolutional Neural Network

A CNN comprises convolutional layers, pooling layers, and fully connected layers. Most CNNs require a fixed input size. So before being sent into the first CNN layer, images need a cropping or warping operation. Convolutional layers and pooling layers do not need a fixed input size. However, the fully connected layers require a fixed input and output size to maintain constant number of the full connections. SPP can pool the mixed-size images into fixed-length outputs, thus meeting the need for fixed inputs in the fully connected layers.

Spatial pyramid pooling (or spatial pyramid matching) was first used in computer vision. Used together with feature extraction and classification algorithms, it has shown good results in image classification [25–27], object recognition [28–30], semantic concept detection [31], and image memorability [32]. Next, the SPP layer was introduced to CNN to remove the fixed-size input constraint of CNN [33]. The SPP-CNN method has been used in remote sensing hyperspectral image classification [34], handwritten word image categorization [35], and action recognition [36]. According to these applications, SPP is useful in CNN. It can reduce the cropping and warping operations used to fit a fixed-size CNN input, and avoid information loss in the operations.

The ordinary pooling layer in CNN is used to compress the input feature maps. It helps to reduce the feature maps and simplify the computational complexity of the network. It also extracts the main features from the original maps. There are two general types of pooling operations: average pooling and max pooling. Figure 3 shows an example of max pooling process in CNN, where *filter* is the filter size that indicates the range of pooling operation. *stride* is the space between pooling operations. It is clear that if the input image size changes more than the stride, the output size will change. This will make the classification algorithm impossible to continue.

**Figure 3.** Max pooling process. Get the maximal value in the range of each pooling operation.

To resolve the requirement for a fixed input size, SPP is introduced to CNN as the last pooling layer. As shown in Figure 4, the feature images are pooled to different levels in the SPP layer. We will get an *l* × *l* size image in level *l*. Thereafter, a fixed-length output can be obtained by *n* level pooling. Feature values *<sup>n</sup>* ∑ *i*=1 *i* <sup>2</sup> will be sent into the fully connected layer for classification.

Spatial pyramid pooling layer

**Figure 4.** Convolutional neural network (CNN) with a spatial pyramid layer. Each input image is pooled to several levels. The results are transformed into one-dimension vectors to form the spatial pyramid pooling (SPP) output.

To get an *l* × *l* size image in level *l*, the filter size and stride should change by level. The *filter* and *stride* can be computed by

$$filter = \left[ m/l \right] \tag{6}$$

$$stride = \begin{bmatrix} m/l \end{bmatrix} \tag{7}$$

where *m* × *m* is the size of the feature maps from the last layer. Table 1 shows an example of 4-level SPP. Two input images with different size 15 × 15 and 20 × 20 get the same 30 length output by a 4-level SPP layer. If the input size changes, the pooling parameters will change to ensure the outputs have the same length.


**Table 1.** Parameters of a 4-level SPP.

Although SPP-CNN has shown good performance in image classification, there are still some problems. As the SPP layer lies at the last, before the fully connected layer in CNN, the outputs of the SPP layer are sent directly into the fully connected layer for classification. However, with no convolutional layer, the features obtained by the SPP layer may not be fully used. In addition, some of the location information will be lost. Meanwhile, the fully connected layer will have a large input matrix. It will greatly increase the connections in fully connected layer. Then there will be much more parameters to be trained. Therefore, a PSPP layer is proposed to make full use of the features and reduce the parameters.

The structure of the PSPP layer is shown in Figure 5. SPP is used to pool the input images into two different levels. Next, the output of the two levels will be used to compose new feature maps rather than a feature vector in the SPP layer. Thus, the output feature maps can be delivered to the convolutional layer for another round of feature extraction.

Pythagorean spatial pyramid pooling layer

**Figure 5.** Structure of Pythagorean spatial pyramid pooling (PSPP) layer. It constructs the output using the results of two pooling levels. The position information of the higher-level pooling result *A* are retained.

To facilitate the composition of two SPP outputs, the pooling levels *a, b* (*a* > *b*) are chosen from the smaller two numbers of the Pythagorean triple. Hence, the output feature maps will have the size of the largest number *c* in the Pythagorean triple. To retain some position information of the feature maps, the composition is processed in the following way. The output matrix of the higher pooling level *<sup>A</sup>* will be retained. Next, the smaller output matrix *<sup>B</sup>* is reshaped as (*<sup>c</sup>* <sup>+</sup> *<sup>a</sup>*) <sup>×</sup> (*<sup>c</sup>* <sup>−</sup> *<sup>a</sup>*). The reshaped matrix is used to expand *A* to *C* on the right side and down side.

Using the PSPP layer in CNN, the fixed input problem can be solved. In addition, the output of the layer are square matrices which can be extracted in the following steps. The structure of the PSPP-CNN used in multi-size training of this paper is shown in Figure 6. Two convolutional layers are added after the PSPP layer for further feature extraction. The convolutional layers will also reduce the size of feature maps. Then the connections in fully connected layer is reduced.

**Figure 6.** Structure of PSPP-CNN used in this paper. It is the PSPP-CNN used in all the following experiments, as PSPP-CNN has ability to process multi-size input.

It is recommended that the PSPP layer be in the middle layers of PSPP-CNN. As the PSPP layer has a larger feature reduction than normal 2 × 2 or 3 × 3 max pooling, the ahead position PSPP layer will lead to the premature loss of features. The PSPP layer at the back position is more like an SPP layer without enough convolutional layers to use the features obtained.

The size of convolutional and pooling kernels is changeable according to the input image size. However, big kernel size may result in information loss and increase computational complexity. Therefore, we choose to add more convolutional layers when processing large input images.

### 2.3.2. Pythagorean Spatial Pyramid Pooling Convolutional Neural Network Training Method

According to the previous description, the forward process is easy to realize. The filter size and stride can be pre-computed before the pooling. However, the back-propagation process in PSPP-CNN training requires some strategy.

When a back-propagation result is received from the last layer, the result is first divided to levels in the same order as the forward process. Next, the result in each PSPP level is restructured as a square matrix. The square matrices apply back-propagation separately. Thus, 2 back-propagation matrices are obtained. Hence, there are two ways to calculate the back-propagation matrix of an PSPP layer: (1) the mean of 2 back-propagation matrices, and (2) the weighted mean of 2 back-propagation matrices according to the level. The calculations can be presented by (8) and (9)

$$d(i) = \frac{1}{2} \sum\_{k=1}^{2} d\_k(i+1) \tag{8}$$

$$d(i) = \frac{1}{\sum\_{k=1}^{2} k^2} \sum\_{k=1}^{2} k^2 d\_k(i+1) \tag{9}$$

where *d*(*i*) is the back-propagation matrix of layer *i* in CNN, *dk*(*i* + 1) is the level *k* back-propagation matrix of layer *i* + 1 which is an PSPP layer. The two methods are tested using the same CNN structure, samples, and learning rate. The samples are part of the data used in Section 3.1. The results are listed in Table 2. Training error rate less than 0.05% is consider as achieving convergence.


**Table 2.** Convergence time and accuracy using two back-propagation methods.

As shown in Table 2, the two methods have similar accuracy and during time of each training step. However, the back-propagation using the second method has a faster convergence rate. This is because that the high-level pooling in PSPP reserves more features from feature maps. Therefore, the high weight of high-level pooling will lead the training to the right direction.

When PSPP-CNN is applied to multi-size images training, an important problem is the training order of the multi-size samples. To prevent the network from fitting a single image size, the multi-size samples in our work will be trained by turns. After all the samples of one size are trained, we will switch to another size. When the training error rates of samples in each size are less than 0.1%, the PSPP-CNN is considered to be achieving convergence.

### **3. Experiment**

To verify the validity of the proposed method, two series of experiments are presented in this paper. Fault diagnosis of constant and variable rotating speed data are conducted using the proposed method.

### *3.1. Constant Rotating Speed Data*

The bearing fault data from the Case Western Reserve University (CWRU) Bearing Data Center [37] is selected to verify the validity of the method in a constant rotating speed environment. The bearing test stand used in the experiment is shown in Figure 7. There are four bearing states: normal, ball fault, inner race fault, and outer race fault. In each bearing fault state, the bearings have fault diameters of 0.007 inches, 0.014 inches, and 0.021 inches. There are also 0.028 inches fault data of ball fault and inner race fault. Thus, there are 12 conditions in total. The fault bearings are installed on the drive end. Three accelerometers are installed on the fan end, drive end, and the base, respectively. The rotating speed of the shaft is about 1800 rpm with the motor load ranging from 0 to 3 HP. All the data selected were sampled at a frequency of 12 kHz.

**Figure 7.** Bearing test stand used by Case Western Reserve University (CWRU).

In our fault diagnosis experiment, all the data are divided into 12 conditions according to the fault states and fault diameters. The influences of fault bearing location, accelerometer location, and motor load are ignored. Because the data were stored as a long array with more than 250,000 data points, the data were divided into several samples. Each sample contains 1024 data points. The size of each condition used for training set and test set are listed in Table 3. The selection of training samples is random. The ratio of training samples to test samples is 2 to 1.

**Table 3.** Sizes of training and test sets in 12 conditions.


Because the sampling frequency of all data is the same, and the change in rotating speed is very small, it can be considered that the sampling frequency of all data is equal to the same multiple of the rotating frequency, namely a 400 multiple. At the same time, because the characteristic frequency of these bearing faults is higher than a two multiple of the rotating frequency, in this experiment, the continuous wavelet transform of the bearing data is carried out from 1 to 200 scales. Hence, 200 × 200 CWTSs are obtained as the CNN input.

To compare the diagnosis effectiveness of the three networks, CNN, SPP-CNN, and PSPP-CNN, on constant rotating speed data, three models are built to diagnose the samples. The structures of the CNNs are listed in Table 4:

As shown in Table 4, the first five layers of the CNNs are set as the same. Because the SPP layer and PSPP layer can reduce more image size than the max pooling layer, the original CNN has more layers and more connections. PSPP-CNN has two more convolutional layers than SPP-CNN for further feature exaction. The selection of CNNs parameters is problem dependent and obtained by trial and error. The selection of the parameters of the first 5 layers was based on the principles proposed in [38]. Then a validation set was built to optimize the parameters of the layers after layer 5 in three networks. The parameters that have best performance on validation set were chosen as the final CNN parameters. The training rate of all three CNNs is set to 0.002 and changed to 0.0005 when the training error is reduced to 1%.



MATLAB is used to implement the training on the computer with two E5-2667 v3 CPUs, a GTX1080Ti GPU, 32 GB memory, and a 1 TB driver. Based on Matconvnet toolkit [39], the SPP-CNN and PSPP-CNN layer are implemented by adding new layer types to it. After the convergence of the CNNs, the test samples are sent into CNNs for fault diagnosis. The convergence and accuracy of CNNs are listed in Table 5:


**Table 5.** Convergence and accuracy of Convolutional neural networks.

From Table 5, it is clear that the accuracy of SPP-CNN is a bit lower than that of CNN in Constant rotating speed data classification. SPP-CNN may lose some important information during SPP layer which is a big feature reduction. However, the convolutional layer following a PSPP layer can re-extract fault features. PSPP-CNN has a diagnosis accuracy similar with the CNN method, and better than SPP-CNN. It has a much shorter training time of each steps, for it has much less parameters to be trained. This shows that the PSPP layer retains the main fault features while reducing the total number of features. The fault diagnosis accuracy of PSPP-CNN in 12 conditions is shown in Table 6. All the conditions have an accuracy greater than 91.03%. Compared with diagnosis accuracies listed in [37], PSPP-CNN has equivalent accuracy, lower proportion of training samples and more conditions. This shows that the method we proposed has good performance in a constant rotating speed data diagnosis.

**Table 6.** Accuracy of test samples in 12 conditions using Pythagorean Spatial Pyramid Pooling Convolutional Neural Network (PSPP-CNN).


### *3.2. Variable Rotating Speed Data*

The fault diagnosis of equipment with variable rotating speed is an important issue that has not been solved satisfactorily. The proposed method is capable of handling variable rotating speed data. Therefore, a variable rotating speed experiment is carried out to show its advantages.

The test bed used in this experiment, the Machinery Fault Simulator-Rotor Dynamics Simulator (MFS-RDS), is shown in Figure 8. Bearing fault experiments were conducted on this test bed. The bearing used is ER-16K LINK-BELT (LBX Company LLC, Lexington, KY, USA). There are four bearing fault modes: normal (NO), ball fault (BA), inner race fault (IR), and outer race fault (OR). The fault bearings can be installed at the drive or non-drive end, and each experiment has at most one fault bearing. Two accelerators are installed on the vertical direction of the two bearing housings. The load is constant in the experiment. Under each bearing fault condition, three sets of data are collected at the rotating speed of 1800 rpm, 2400 rpm, and 2900 rpm, respectively, and the sampling time of each set is approximately 10 min. The sampling frequency of all data is 12 kHz.

**Figure 8.** Machinery Fault Simulator-Rotor Dynamics Simulator (MFS-RDS) test bed. It can be used to simulate shaft, motor and bearing faults. The eddy current sensors are used to monitor the state of shaft. The data of them are not used for bearing diagnosis experiment. (1) speed controller, (2) rigid coupling, (3) accelerometer, (4) electromotor, (5) bearing base, (6) bearing, (7) shaft, (8) eddy current sensor, (9) rotary table.

In the fault diagnosis experiment, the data are divided into 12 cases according to the fault status and rotating speeds. The influence of fault bearing position and sensor position is neglected. Because the data are stored as continuous time series, the data are divided into several samples. Each sample contains 1200 data points. Next, 5376 samples are obtained in each case; half of them, 2688 samples, are used for training, and the remaining 2688 samples are for test. The ratio of training samples to test samples is 1:1. Therefore, we have 32,256 training samples and 32,256 test samples totally. The fault diagnosis aims to classify the data into four categories based on the fault status.

Because all data samples have the same sampling frequency of 12 kHz, this corresponds to 400, 300, and 248 multiples of the rotating speed at 1800 rpm, 2400 rpm, and 2900 rpm, respectively. Thus, the continuous wavelet transform of three sets of bearing data is carried out from 1 to 400 scales, 300 scales, and 248 scales, respectively. In the time axis, the middle 400, 300, and 248 coordinates are chosen, because they have data points of one rotating period and can neglect the first and last few points of each sample. Hence, the CWTSs are chopped to square CWTSs that have a size of 400 × 400, 300 × 300, and 248 × 248, respectively.

The PSPP-CNN structure used in this experiment is shown in Figure 4 and Table 4. There are five convolutional layers, three max pooling layers, one PSPP layer, and one fully connected layer. The training rate is initially set to 0.002 and changed to 0.0005 when the training error is reduced to 1%. The training environment is the same as the constant rotating speed data training. It takes 317 min to achieve convergence after 44 training steps which means the error of training samples is less than 0.1%.

The confusion matrix of fault diagnosis result is shown in Table 7. The first row represents the rotating speed and labels of the test data. The first shows the diagnosis result labels. The method has a high diagnosis accuracy of 99.11%, and an accuracy of more than 97.61% is obtained for the data of each fault condition. This indicates that the PSPP-CNN method is suitable for bearing fault diagnosis with variable rotating speed.

To compare the diagnosis effect with CNN and SPP-CNN, two other models are built using the CNN and SPP-CNN structures, as listed in Table 4. As CNN can only accept fixed-size images, all the data to 400 × 400, 300 × 300, and 248 × 248 CWTSs are translated to train the CNN separately. Accordingly, the CNN structure changes by adding a convolutional layer of a different size before the fully connected layer. SPP-CNN uses the same input images as those used by PSPP-CNN. 400 × 400, 300 × 300, and 248 × 248 CWTSs are also used to train SPP-CNN and PSPP-CNN. The diagnosis accuracy of seven different CNN models are listed in Table 8 and Figure 9.


**Table 7.** Accuracy of test samples at different rotating speeds using PSPP-CNN.

**Table 8.** Accuracy of test samples using different CNN models.


**Figure 9.** Accuracy of test samples using different input size and CNN models.

As shown in Table 8, PSPP-CNN has a better diagnosis accuracy than other CNN and SPP-CNN models with single-size or multi-size input. Figure 9 shows that the accuracy increases along with the size and diversity of input samples. The method we propose using PSPP-CNN and multi-size input has the best diagnosis accuracy 96.58%. It shows a nearly two-percent improvement over SPP-CNN. The SPP-CNN model with multi-size does not get a big accuracy improvement. It is because that different from image recognition tasks, fault diagnosis needs precisely feature positioning in CWTS. A SPP layer before fully connect layer will lose more position information in CWTS than a PSPP layer which can locate as front layer of PSPP-CNN. In addition, using the normal CNN, the diagnosis accuracy increases with the increase of input image size. According to our analysis, this is mainly owing to the increase of 1800 rpm data accuracy which are 92.05%, 94.22%, 95.53% separately. The reason may be that the 300 × 300 and 248 × 248 cropping of CWTSs will lose some of the fault features of 1800 rpm data. An input size larger than 400 × 400 will not increase the accuracy.

To compare this PSPP-CNN method with other CNN-based fault diagnosis methods, several proposed methods are applied using the same dataset. Deep Convolution Neural Network with Wide first-layer kernels (WDCNN) [12] is a bearing diagnosis method that uses raw vibration signals as input, and wide kernels in first layer for feature extraction and high-frequency noise suppression. Dislocated Time Series Convolutional Neural Network (DTS-CNN) [18] is proposed for

mechanical signals by adding a dislocate layer to CNN. DTS-CNN can extract the relationship between signals with different intervals in periodic mechanical signals. Resample-CNN [38] uses a resample method to normalize the data to make sampling frequencies are the same multiples of the rotating frequency. The three proposed methods are used to compare the fault diagnosis accuracy using the same dataset, for they are all proposed for the diagnosis of machinery using vibration signals based on CNN. Because the data used in the papers have similar testbed structure, sampling frequency and rotating speed, WDCNN and DTS-CNN use the same network parameters as those in the papers. Resample-CNN uses the 400 × 400 CNN structure for they have the same input size. As all the three methods use the ordinary CNN and have the similar CNN structures, the changes of CNN parameters will not change the diagnosis result greatly. What we need to focus on is the construction of deep learning input. As the accuracies shown in Table 9, the method proposed in this paper has the best performance in the diagnosis of variable rotating speed data.

**Table 9.** Fault diagnosis result using other proposed CNN-based methods.


The reason for the comparison made on this set of data is that the accuracies on this case can reflect the effectives of the methods both on constant speed and variable speed data. In addition, in practical use, most variable speed machines work like this case at a certain speed range or some preset optimal speed points.

To further study the effectiveness of the method at full working speed range, an experiment to diagnosis the full working speed data is carried on using the PSPP-CNN trained above. The data of four bearing conditions were collected with 12 kHz on MFS-RDS. The fault bearing was installed on drive end. The rotation speed ranges from 300 to 3000 r/min. The rotation frequency increases at 0.15 Hz/s. Figure 10 shows the vibration signals of drive end accelerator in four fault conditions.

**Figure 10.** Vibration signals of four fault conditions. (**a**) normal (**b**) ball (**c**) inner race (**d**) outer race.

To test the PSPP-CNN trained, the data of drive end accelerator in each condition are divided into 440 samples. Each sample contains 8192 data points. We get the cropped CWTS of each sample according to rotating speed. Therefore, the size of cropped CWTSs are ranges from 240 × 240 to 2400 × 2400. Then the cropped CWTSs are sent into PSPP-CNN for fault diagnosis. The accuracy of each fault condition is listed in Table 10.


**Table 10.** Accuracy of four fault conditions at full rotating speeds using PSPP-CNN.

As shown in Table 10, PSPP-CNN has diagnosis accuracies more than 90% for each fault conditions. It means that the PSPP-CNN trained by data at some rotating speed can be used to diagnosis bearing fault in full working speed. Through analysis, the accuracy of data under 1200 rpm is a little lower. Adding an incremental training using low speed data will increase the accuracy. It shows that the PSPP-CNN trained using data of few certain rotating speeds can be used to diagnose bearing fault in full working speed.

Through the fault diagnosis experiments of constant and variable rotating speed data, we can know that the PSPP-CNN method proposed is an effective solution for fault diagnosis of bearing. When applied to intelligent diagnosis system, the method has some advantages. First, the PSPP-CNN proposed in this paper can be easily implemented by adding a PSPP layer to the ordinary CNN code based on max pooling layer. There are some mature CNN frameworks based on MATLAB, Net Framework or Python. All of them can be easily built. Second, the fault diagnosis process has high power efficiency. As shown in Table 5, PSPP-CNN has less parameters than an ordinary CNN with the same front layers. It reduces the computation of each training and test epoch. Although the training of PSPP-CNN is still time-consuming in the computer without GPU accelerated computing, the diagnosis process using trained PSPP-CNN model can be completed quickly even at laptop computer. Third, with the de-noising ability of wavelet transform, the diagnosis system has good robustness. Diagnosis result will not be affected by the background noise in signal of practical equipment.

In practical applications, we can gather the experiment data or online monitoring data from the varying working speed bearing and use the data to train PSPP-CNN for fault diagnosis. The intelligent fault diagnosis method proposed in this paper has been used for online fault diagnosis of wind turbine bearings in a wind farm. The fault diagnosis software is exploited by C#& MATLAB combined programming. The signal is transmitted to MATLAB for CWTS calculation and classification. The diagnosis result can be obtained in 5 s. The training data was collected from the online vibration monitoring system installed on wind turbines. The fault data was picked out by referencing fault records.

### **4. Conclusions**

In this paper, we propose an intelligent fault diagnosis method for a variable rotating speed bearing. The proposed approach is built upon CWTS and PSPP-CNN. This method decomposes vibration signals of bearing into CWTSs of different scales according to the rotating speed. In addition, the different size CWTSs are sent into PSPP-CNN for fault diagnosis. The PSPP-CNN that we proposed is an improvement of the SPP-CNN. The PSPP layer can fully use the pooling result for further feature extraction than SPP layer as they all can pool the input of different sizes to a fixed size. A series of experiments are carried out using constant rotating speed data and variable rotating speed data. The results show that the proposed approach is an effective solution. The method has been used in practical applications.

Although many fault diagnosis methods based on deep learning have been proposed, most methods are totally data-driven and focus on the small improvement of deep learning algorithm. The domain knowledge that has been used for fault diagnosis in recent decades is barely used. In addition, the working condition information is not considered in the input of deep learning algorithm. This paper takes the fault characteristics frequency of bearing and working speed into consideration. However, more domain knowledge and working condition information, such as load and output power, can be combined with deep learning. It may improve the accuracy and robustness, and the features extracted will be more interpretable.

**Author Contributions:** T.Y. and W.G. conceived this study. S.G. and T.Y. designed the experiments. S.G. and C.Z. performed the experiments. S.G. analyzed the data and wrote the paper. T.Y., W.G. and Y. Z. reviewed the paper.

**Funding:** This research was supported by the National Program of International Science and Technology Cooperation of China (Project No. 2014DFA60990).

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

### **References**


© 2018 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 (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Collision Detection and Identification on Robot Manipulators Based on Vibration Analysis**

**Feiyan Min 1,2, Gao Wang 1,2 and Ning Liu 1,2,\***


Received: 25 December 2018; Accepted: 27 February 2019; Published: 3 March 2019

**Abstract:** Robot manipulators should be able to quickly detect collisions to limit damage due to physical contact. Traditional model-based detection methods in robotics are mainly concentrated on the difference between the estimated and actual applied torque. In this paper, a model independent collision detection method is presented, based on the vibration features generated by collisions. Firstly, the natural frequencies and vibration modal features of the manipulator under collisions are extracted with illustrative examples. Then, a peak frequency based method is developed for the estimation of the vibration modal along the manipulator structure. The vibration modal features are utilized for the construction and training of the artificial neural network for the collision detection task. Furthermore, the proposed networks also generate the location and direction information about contact. The experimental results show the validity of the collision detection and identification scheme, and that it can achieve considerable accuracy.

**Keywords:** manipulator; model independent method; collision detection; collision identification; vibration analysis; artificial neural network

### **1. Introduction**

Industrial robots play an important role in the modern manufacturing industry, and humanfriendly robots will soon become flexible and versatile coworkers in the industrial setting [1]. One of the core problems in human–robot interaction is the detection of collisions between robots and the industrial environment, including humans and other manufacturing structures. Indeed, industrial robots should be able to operate in very dynamic, unstructured, and partially unknown environments, sharing the workspace with the human user, and preventing upcoming and undesired collisions [2].

Furthermore, researchers are getting more interested in gathering the maximum amount of information from the impact event, such as the contact position, direction and intensity, in order to let the robot react in the most appropriate fashion. Related concepts include collision avoidance [3], and collision isolation, identification, classification and reaction [2]. Particularly, the concept collision isolation aims at localizing the contact point, or at least which link out of the n-body robot collided. Furthermore, collision identification is to determine the directional information of the generalized collision force [2].

Different approaches for detection of robot collisions have been presented in the literature. A first intuitive approach is to monitor the current transient in robot electric drives, looking for shock changes within currents caused by collisions [4–7]. A second approach is based on the tactile sensors laying inside robot skins. The more common approach is the so-called model-based method (the state observer or Kalman filter method), and the detection algorithms are mainly based on the evaluation of monitoring signals (motor currents, difference between actual and predicted torques, etc.), which should be below some setting values, otherwise collision alarms are generated.

A major practical problem in these methods is the selection of thresholds for the monitoring signals, since the modeling error and sensor noise affect the monitoring signal in the same way as collision disturbance. A good detection algorithm therefore must distinguish the effect of modeling error and sensor noise on monitoring signal from that of a real collision. For this reason, it usually leads to a tradeoff between sensitivity and false alarm rate, with a risk of excessively conservative threshold [2,5]. To overcome this problem, some different methods are proposed. A dynamic threshold is defined in [4] to represent the residual dependence on the state of the robot (position, velocity, acceleration) using fuzzy logic rules. Furthermore, in [6], authors propose an adaptive detection algorithm based on a state-dependent dynamic threshold. In recent years, some extended state observer [8] and sliding mode observer [9] methods are proposed to obtain more effective detection performance.

Most of these evaluation algorithms mainly concentrate on the time domain information of motor torque deviation (The motor currents can be considered as another form of torque). A considerable alternative is to use frequency domain features for detection purposes.

This paper describes a novel detection method based on vibration modal feature generated by collision. The natural frequency and vibration modal features of collision experiments are extracted for the constructing and training of several structures of Back Propagation (BP) neural network. The research result shows that this method can be used for the detection of collision with considerable accuracy, not only for the detection of collision occurrence, but also for the positioning and direction determination.

A remarkable contribution of this paper is to introduce frequency and vibration information for the detection of robot collision. The vibration information is independent of the dynamic model and can be easily acquired with acceleration sensor, which has greatly developed in recent years. Since the occurrence of collision can happen anywhere along the robot structure, the acceleration sensor is easy to stall and suitable for the detection use. Furthermore, the industrial robots are not usually equipped with torque sensor because of cost and structural constraints, and this research provides an opportunity for MEMS accelerometer.

Accurate analysis of natural frequency and modal shape is fundamental for mechanical design, dynamic identification and control of high-speed manipulator [10,11]. In recent years, the acceleration and vibration analysis method has been gradually used for the status monitoring of manipulator robot. In [12], the authors propose a fault detection method for industrial welding robot. In their study, joint acceleration of robot is considered as evaluation criteria and their evaluation algorithm is based on neural network. The article [13] presents a method for processing and analyzing the measurement signals used in the problem of diagnosing the state of a manipulator's tool. The analysis algorithm is performed within the time and frequency domain. The signal utilized in the research is the mechanical vibrations and the rotation speed of the tool. In the work [14], frequency domain analysis method is researched for the event classification in robot assisted deburring. Power spectrum density (PSD) of sensor data acquired between certain sample rates is calculated, and it is then used for classifying vibration signal generated by the spindle from the vibration signal acquired during the deburring process. The paper [15] presents a signal fusion method based on accelerometer and encoder in serial robots. Besides, a number of studies have presented the data-based method for mechanical system monitoring based on vibration signals [16–20]. In most of this research, specific forms of artificial neural network are proposed and implemented [21–25].

The rest of this paper is organized as follows. Section 2 describes the frequency domain features and parameters of vibration of manipulator in case of collisions. Section 3 presents the architecture of our detection neural network, together with vibration modal analysis method; Section 4 presents some results and experiments obtained with the proposed method; and finally, Section 5 addresses the main conclusion and future work.

### **2. Vibration Modeling and Feature Extraction**

This section address the vibration presentation and its characteristics. First, the modeling method of vibration response is proposed with dynamic equations and transfer functions. Then, the typical features are analyzed based on the mathematical model, together with illustrative examples.

### *2.1. Vibration Response under Collision and Its Mathematical Modeling*

In this paper, our research focuses on the vibration response of multiple test points along manipulator structure under several experiments. An effective method for dealing with vibration of kinematic chain mechanism is the elastodynamic modeling and analysis [10,11]. In this method, each critical mechanical structure and component is simplified with stiffness, viscous and mass parameter. As for a *n*−dof manipulator, it is composed of a series of motors, gears, links and joints, some of which will generate deformation and vibration under collision. For this reason, the dominant vibration structures can be considered as elastic bodies with certain stiffness and viscous coefficients.

We consider a *n*−dof manipulator with *m* dominant vibration structures. Its axis displacement vector can be defined as:

$$q = \left[ \begin{array}{c} q\_D \\ q\_M \end{array} \right] \in \mathbb{R}^N \tag{1}$$

where *qD* = [*q*1, *q*2, ··· , *qm*] *<sup>T</sup>* is the vibration deviation, and *qM* <sup>=</sup> [*qm*+1, *qm*+2, ··· , *qm*+*n*] *<sup>T</sup>* denotes the joint displacements. Furthermore, we assume the equilibrium points of *qD* is *qD* = [*q*1, ··· , *qm*] *T*. The dynamic equation of manipulator can be written as:

$$M(q)\ddot{q} + \mathbb{C}(q, \dot{q})\dot{q} + \mathbb{G}(q) = \begin{bmatrix} \tau\_D\\ \tau\_M + \tau\_f \end{bmatrix} \tag{2}$$

the variable *τM*, *τ<sup>f</sup>* and *τ<sup>D</sup>* denotes the joint torque generated by motor, friction and structure deformation respectively. We denote *τ<sup>D</sup>* = [*τ*1, ··· , *τm*] *<sup>T</sup>*, and *<sup>τ</sup><sup>M</sup>* <sup>=</sup> [*τm*+1, ··· , *<sup>τ</sup>m*+*n*] *<sup>T</sup>*. The subscript *D* denotes the dominant vibration structures, and *M* denotes the drive motor of manipulator.

On the other hand, the torque generated by the deformation of vibration structures can be given by:

$$
\pi\_D = K\_p(\overline{q\_D} - q\_D) - K\_v q\_D \tag{3}
$$

where *Kp* = *diag*(*kp*1, *kp*2, ··· , *kpm*), *Kv* = *diag*(*kv*1, *kv*2, ··· , *kvm*) is the vibration dynamic coefficient matrix. Furthermore, *kpi*, *kvi* denotes the stiffness and vicious coefficient of *ith* dominant vibration structure respectively.

We assume that the gravity and friction is compensated by feedback control loop, and the Coriolis and centrifugal effect generated by structure deformation is relatively small. In most cases, the displacement of *qD* is much smaller than *qM*, and the inertia matrix *M*(*q*) is mainly determined by joint displacement. Furthermore, the torque vector generated by collision is assumed as *fext*. Then we get the dynamic function of manipulator as follows

$$
\begin{bmatrix} 0 \\ \tau\_M \end{bmatrix} = M(q\_M) \begin{bmatrix} q\_D \\ q\_M \end{bmatrix} + \begin{bmatrix} K\_{\overline{\nu}} & \mathbb{C}\_D(q\_M, q\_M) \\ 0 & \mathbb{C}\_M(q\_M, q\_M) \end{bmatrix} \begin{bmatrix} q\_D \\ q\_M \end{bmatrix} + \begin{bmatrix} K\_p \\ 0 \end{bmatrix} (q\_D - \overline{q\_D}) + \begin{bmatrix} J\_D^T \\ J\_M^T \end{bmatrix} f\_{\text{ext}} \tag{4}
$$

where *JD<sup>T</sup>* and *JM<sup>T</sup>* is the associated geometric contact Jacobian matrix to *m* vibration structures and *n* robot motors, respectively. However, the collision torque vector *fext* and contact Jacobian *JD<sup>T</sup>*, *JM<sup>T</sup>* is typically unknown.

Denoting *x*<sup>1</sup> = *q*˙*D qD* − *qD* , *y*<sup>1</sup> = *qD* − *qD*, and *x*<sup>2</sup> = *q*˙*<sup>M</sup> qM* , *y*<sup>2</sup> = *q*˙*M*, an *n* inputs and *m* + *n* outputs state-space equation is obtained

$$\begin{cases} \dot{x}\_1 = A\_1 x\_1 + B\_{11} y\_2 + B\_{12} f\_{ext} \\ \dot{x}\_2 = A\_2 x\_2 + B\_{21} u + B\_{22} f\_{ext} \\ y\_1 = C\_1 x\_1 \\ y\_2 = C\_2 x\_2 \end{cases} \tag{5}$$

where *u* is the active motor torque, and *fext* is torque vector generated by collision, and the related parameter matrix is given by:

$$\begin{aligned} A\_{1} &= \begin{bmatrix} -\mathcal{M}(q\_{M})^{-1}K\_{\upsilon} & -\mathcal{M}(q\_{M})^{-1}K\_{p} \\ I & 0 \end{bmatrix}, \\ B\_{11} &= \begin{bmatrix} -\mathcal{M}(q\_{M})^{-1}\mathcal{C}\_{D}(q\_{M}, \dot{q}\_{M}) \\ 0 \end{bmatrix}, \quad B\_{12} = \begin{bmatrix} I\_{D}^{T} \\ 0 \end{bmatrix}, \quad \mathcal{C}\_{1} = \begin{bmatrix} 0 \\ I \end{bmatrix}^{T}, \\ A\_{2} &= \begin{bmatrix} -\mathcal{M}(q\_{M})^{-1}\mathcal{C}\_{M}(q\_{M}, \dot{q}\_{M}) & 0 \\ I & 0 \end{bmatrix}, \\ B\_{21} &= \begin{bmatrix} -\mathcal{M}(q\_{M})^{-1} \\ 0 \end{bmatrix}, \quad B\_{22} = \begin{bmatrix} I\_{M}^{T} \\ 0 \end{bmatrix}, \quad \mathcal{C}\_{2} = \begin{bmatrix} I \\ 0 \end{bmatrix}^{T} \end{aligned}$$

Then we get the modal analyzed transfer function from collision torque *fext* to vibration deformation *y*<sup>1</sup> as follows

$$P(\mathbf{s}) = \frac{Y\_1(\mathbf{s})}{F\_{\text{ext}}(\mathbf{s})} = \mathbb{C}\_1(\mathbf{s}I - A\_1)^{-1} \left[B\_{12} + B\_{11}\mathbb{C}\_2(\mathbf{s}I - A\_2)^{-1}B\_{22}\right] \tag{6}$$

We considered the structure of system matrix *A*<sup>1</sup> and *A*2, and we got that *rank*(*A*1) = 2*m* and *rank*(*A*2) = *n*. Let *λDi*, *λDi*(*i* = 1, 2, ··· , *m*) be 2*m* eigenvalues of *A*1, and *λ*0(= 0), *λMi*(*i* = 1, 2, ··· , *n*) be the *n* + 1 eigenvalues of *A*2.

From the partial fraction expansion, the modal analyzed transfer function *P*(*s*) becomes

$$P(s) = \sum\_{k=1}^{m} \frac{\Phi\_k}{s^2 + 2\xi\_k \omega\_{Dk} s + \omega\_{Dk}^2} + P\_1(s) \tag{7}$$

By modal analysis, *P*(*s*) can be expressed as a linear sum of *m* + *n* vibration modes, *m* modes of which are generated from the *m* dominant vibration structures. Furthermore, the other *n* poles which come from *A*2, correspond to the *n*-dof active dynamic of manipulator.

Then the natural frequency of vibration along manipulator can be obtained by *ωDi* = |*λDi*|(*i* = 1, 2, ··· , *m*) and *ωMj* = |*λMj*|(*j* = 1, 2, ··· , *n*).

The matrix Φ*<sup>k</sup>* corresponds to the rigid body vibration mode and is positive semidefinite, and it has the following structure

$$
\Phi\_{\cdots i} = \begin{bmatrix}
\Phi\_{11i} & \Phi\_{12i} & \cdots & \Phi\_{1ni} \\
\Phi\_{21i} & \Phi\_{22i} & \cdots & \Phi\_{2ni} \\
\cdots & \cdots & \cdots & \cdots & \cdots \\
\Phi\_{l1i} & \Phi\_{l2i} & \cdots & \Phi\_{lni}
\end{bmatrix} \tag{8}
$$

For any *φkji* ∈ Φ··*i*, the subscript *k* ∈ (1, 2, ··· , *l*) denotes the serial number of vibration test position, *j* ∈ (1, 2, ··· , *n*) denotes the serial number of collision torque in vector *fext*, and *i* ∈ (1, 2, ··· , *m*) is the number of vibration modes.

We considered the structure of system matrix *A*<sup>1</sup> = −*M*(*qM*) −1 *Kv* −*M*(*qM*) −1 *Kp I* 0 , it is mainly determined by inertia matrix *M*(*qM*), vibration dynamic coefficient matrix *Kp* and *Kv*, of which the latter two remain nearly unchanged. That means the natural frequencies *ωDi*(*i* = 1, 2, ··· , *m*) will vary with joint displacement *qM*.

As the stiffness coefficients *kpi*(*i* = 1, 2, ··· , *m*) of manipulator structure are usually considerably large, there exist several eigenvalues of *A*<sup>1</sup> relatively larger than that of *A*2, which comes from the active dynamic of robot and reflects the normal working frequency of manipulator. For this reason, there exist several natural frequencies of vibration that are obviously bigger than the frequency band of active dynamic of robot, i.e., ∃*ωDi*(*i* = 1, 2, ··· , *m*), and *ωDi* > *ωMj*(*j* = 1, 2, ··· , *n*) hold.

So we get an important conclusion: the dominant natural frequency of vibration under collisions is independent of the dynamic property of robot. When the inertia matrix is given, the natural frequency under collision should remain the same during different dynamic processes or static statuses.

Furthermore, the modal matrix Φ*<sup>k</sup>* depends on the value of geometric contact Jacobian matrix *JD<sup>T</sup>* and *JM<sup>T</sup>*. For this reason, the contact position information can be induced by the value of vibration mode matrix Φ*k*.

### *2.2. Vibration Features Analysis and Illustrative Example*

This section illustrates the characteristics of vibration modal shape of manipulator under collisions with example. The test data were collected on the STR6-05 robot arm (see Figure 1), a 6-DOF heavy-load manipulator. Four consecutive collision experiments were conducted while the end-effector moved in line; the collision conditions are listed in Table 1.


**Table 1.** Experiment conditions for vibration modal test.

**Figure 1.** Experiment setup.

The joint displacement and motor current signals of joint 1, 2 and 3 are shown in Figure 2. We label the contact time of corresponding collision events on the figure. We find that there are some shock changes in the motor currents but little change in joint displacements. However, it is hard to detect collision events from current signals directly because of dynamic change and the noise involved in signals, particular for heavy-load manipulator. As shown in the figure, we cannot find some obvious features for collision 1.

**Figure 2.** Joint displacement and current during the experiment.

For vibration signal measuring, the 1A113E and 1A114A industrial accelerometer (made by Donghua Testing Technology) with NI-9232 data acquisition module is used. The 1A113E accelerometer (uniaxial) is mounted beside joint 2, perpendicular to the direction of joint 1 and joint 2. The 1A114A accelerometer (triaxial) is equipped beside the end-effector of manipulator.

The vibration acceleration and vibration modal is shown in Figure 3. The first signal comes from 1A113A located beside joint 2, and the latter comes from 1A114A beside end-effector, corresponding to two perpendicular directions.

**Figure 3.** Vibration modes in experiment.

The frequency characteristic of each accelerometer within sliding windows is shown in Figure 3. There are four obvious frequency charts corresponding to four collisions in each vibration signal. The peak frequencies, i.e., the natural frequencies of three accelerations corresponding to one collision

are approximately the same, with different energy densities. The magnitudes of all the peak frequencies constitute the modal matrix Φ in Equation (8).

Figure 3 shows that the vibration modes of different parts along robot structure have the following characteristics:


Limiting the range of magnitude below 1 dB, we got the frequency characteristic of normal dynamic comparative to collisions in Figure 4. It shows that the natural frequency of active dynamic appears at low frequency segment (always below 50 Hz), while the vibration frequency by collision appears at intermediate segment. That means the eigenvalues of *A*<sup>2</sup> is much smaller than that of *A*<sup>1</sup> in Equation (5), and the dynamics of normal operation do not affect the natural frequencies of collisions.

**Figure 4.** Low magnitude segment of vibration modes in experiment.

Obviously, the natural frequencies and modal is mainly dependent on the inertia matrix, and independent of the robot dynamic. This means that the detection algorithm can be designed without considering the dynamic property of the robot, and the training and test samples for the model independent method can be collected in some simple or static scenarios.

Since collision will generally cause shock vibration, and its natural frequencies are mainly contained in the high-frequency domain, several symptom parameters in the frequency domain can be selected to represent the collision event. Extracting collision information from frequency domain signal requires proper understanding of the process. For frequency features, such as natural frequencies, the vibration modal shown in the spectrum often has direct or indirect connection to certain dynamic events.

### **3. Collision Detection and Identification Method Based on Vibration Features**

Based on the vibration features discussed above, we proposed a learning-based algorithm for the detection, isolation and identification of collisions.

The proposed method mainly contains three parts: vibration mode analysis, collision detection, and collision identification. As shown in Figure 5, the vibration signals are first recorded by accelerometer sensors. The vibration mode related features are then extracted from the vibration signals. Any change in the vibration mode from normal condition can indicate the occurrence of collision. Three different Back Propagation neural networks are developed for the detection of collisions (BP1), the isolation of contact part (BP2), and identification of direction (BP3). BP2 and BP3 should be activated only if some collisions are detected by BP1, and the input of BP3 varies with the output of BP2.

**Figure 5.** Vibration signal based detection framework.

### *3.1. Vibration Mode Estimation*

The main objective of vibration mode estimation is to determine the main natural frequencies *ωDi*(*i* = 1, 2, ··· , *m*), and the values of vibration mode corresponding to each natural frequency *φkji* ∈ Φ··*i*, where *k* denotes the serial number of sensors, *j* denotes the serial number of collision force, and *i* ∈ (1, 2, ··· , *m*) denotes the number of vibration modes.

In the collision experiment, the result of measurement is strings of acceleration values in discrete moments. The acceleration signal should be properly processed in order to build the collision classifier. For the discrete string *a*1(*k*), *a*2(*k*), and *a*3(*k*), fast Fourier transform (FFT) is used to determine vibration spectrum within a sliding window which has lower computational cost and better accuracy performance. In our research, the sampling rate is 3.2 k Hz, and the width of sliding window is 320 samples with 50% overlap. Furthermore, the cycle time of detection algorithm is 0.1 s.

For the spectrum of each sampling window, we proposed a peak frequency-based method to determine the natural frequency *<sup>ω</sup>Di*(*<sup>i</sup>* <sup>=</sup> 1, 2, ··· , *<sup>m</sup>*), the estimation of *<sup>ω</sup>Di*. For each, the proposed approach consists of the following steps:

*Step 1: For acceleration signal k* ∈ (1, 2, ··· , *l*)*, set a tolerance error* Δ *and δ for spectrum analysis;*

*Step 2: Find all the local maximum power and local minimum power from start frequency f*<sup>0</sup> *to cut-off frequency fe, such that there is no other power value larger than current maximum value between the two adjacent local minimum powers, and there is no other power value less than current minimum value between the two adjacent local minimum power densities; The difference between adjacent maximum and minimum power density should be larger than* Δ*;*

*Step 3: Collect all the local maximum power values and corresponding frequencies;*

*Step 4: Calculate Num, the number of local maximum power density;*

*Step 5: If Num* > *m, then* Δ = Δ − *δ, repeat step 2 to step 4;*

*Step 6: Record current local maximum power values Pki*(*i* = 1, 2, ··· , *m*) *and corresponding frequencies fki*(*i* = 1, 2, ··· , *m*)*;*

*Step 7: For other acceleration signal k* ∈ (1, 2, ··· , *l*)*, repeat step 1 to step 6*;

*Step 8: Rank all the frequency values fki*(*k* = 1, 2, ··· , *l*; *i* = 1, 2, ··· , *m*)*, and divide them into m groups with maximum intervals.*

*Step 9: Calculate the average frequencies within each group, these values are the estimation of modal frequency, <sup>ω</sup>Di*(*<sup>i</sup>* <sup>=</sup> 1, 2, ··· , *<sup>m</sup>*)*.*

With the estimation of modal frequencies, the vibration modal matrix Φ can be determined by extracting the magnitude of corresponding frequency in the spectrum chart, as it shows in Figure 6. In this figure, we get 6 natural frequencies *<sup>ω</sup><sup>D</sup>*1, *<sup>ω</sup><sup>D</sup>*2, ··· , *<sup>ω</sup><sup>D</sup>*<sup>6</sup> from the measurements of three acceleration signals in collision experiment *j*. The magnitude values at each characteristic frequency are the estimation of vibration modal. Obviously, there is little error of estimation, e.g., the estimation *φ*<sup>2</sup>*j*5. However, this error has limited influence on the final detection result because this detection is based on the synthesis of multiple vibration modal. In this way, we get the estimation of vibration modal matrix of experiment *j*:

$$
\hat{\Phi}\_{\cdot j\cdot} = \begin{bmatrix}
\hat{\Phi}\_{1/1} & \hat{\Phi}\_{1j2} & \cdots & \hat{\Phi}\_{1j6} \\
\hat{\Phi}\_{2j1} & \hat{\Phi}\_{2j2} & \cdots & \hat{\Phi}\_{2j6} \\
\hat{\Phi}\_{3j1} & \hat{\Phi}\_{3j2} & \cdots & \hat{\Phi}\_{3j3}
\end{bmatrix} \tag{9}
$$

**Figure 6.** Joint displacement and current during the experiment.

### *3.2. Proposed Artificial Neural Network*

In our method, 3 BP networks are implemented for the collision detection, positioning and direction identification respectively. The collision detection artificial neural network together with the modal analysis algorithm should be executed within each sliding window. Once a collision event is detected, the collision positioning artificial neural network is launched with the current vibration modal data. In addition, the input of the 3rd network is dependent on the output of collision position information. The operation process of relevant algorithms is displayed in Figure 7.

**Figure 7.** Procedure of detection algorithm.

B-P ANN consists of an input layer, hidden layers, and an output layer of neurons. A neuron serves as a processing unit in which output is a linear or nonlinear transformation of its inputs. The neurons, as a group, serve to map the input vibration modal features to the desired collision patterns. The structure of BP network is shown in Figure 8.

The output signal of hidden layer and output layer can be described in the following equations:

$$m\_{\dot{j}}(t) = f(\sum\_{j=1}^{10} w\_{\dot{i}\dot{j}}(t)t + b\_{\dot{j}}) \tag{10}$$

where *mj*(*t*) is the output of current neurons, *wij* is the weight of the connection between current layer neurons and its input layer neurons, *bj* is the bias of the jth neuron of current layer. The activation function of output layer is linear function, while the hidden layer uses sigmoid function

$$f(t) = \frac{1}{1 + e^{-t}} \tag{11}$$

**Figure 8.** Back Propagation (BP) artificial neural network structure.

Considering the characteristic of vibration modes discussed in Section 2, we selected the most important features for each kind of detection task, as shown in Tables 2–4.

As the natural frequencies are mainly dependent on inertia matrix *MqM*, we selected the displacement of joint 2 and joint 3 as input features of detection. The displacements of other joints have little influence on inertia matrix.


**Table 2.** The features used for the detection of collisions.

The contact position and direction of collision can affect the comparative vibration magnitude of sensor on different positions and directions. We select relative modal between test position for the positioning of contact, and relative modal between test direction beside contact position for direction identification.

**Table 3.** The features used for the isolation of contact position.



### **4. Experiment and Discussion**

### *4.1. Experiment Dataset Preparation and Training Procedure*

The experiment dataset is collected on the platform in Figure 1. The robot is controlled with PD control law during the experiment. An uniaxial accelerometer is mounted beside joint 2, perpendicular to the direction of joint 1 and joint 2. In addition, a triaxial accelerometer is equipped beside the end-effector of manipulator. All acceleration data is recorded by NI-9232 data acquisition module. Our experiments are conducted with an aluminum impact hammer, which is also connected to NI-9232 data acquisition module. In this way, the collision time and force data is recorded, as shown in Figure 9.

**Figure 9.** Data acquisition module for collision experiment.

The collision experiments are performed with eight selected contact points, with different direction and force intensity. As the natural frequencies of vibration are mainly dependent on inertial matrix and joint displacement (see Section 2.1), we select five typical working patterns (see Figure 10) as testing standard. Our collision experiments are performed while the robot is transforming randomly from one pattern to another.

**Figure 10.** Five testing patterns of manipulator.

Hundreds of collision experiments are performed on the platform. The vibration modal data of acceleration signals during these experiments are analyzed with the sliding-window method introduced above. Typical vibration modal data sets with corresponding experiment conditions are collected for classification research. Meanwhile, we also collect the vibration modal data during collision-free operations. We get in total 800 experiment samples, of which 50 percent are collision samples performed on 8 different contact positions along robot structure. Eighty-five percent of these samples are selected randomly for training and the remaining are used for testing.

### *4.2. Detection and Identification Results*

In this subsection, several test divisions of experiment data are utilized to evaluate the efficiency of the proposed method. Our BP-ANN algorithms are taken from the Neural Network toolbox in Matlab. At the training stage, we optimize the weights and bias parameters by minimizing mean squared error, according to Levenberg-Marquardt optimization.

For the detection of collisions, BP networks with different architecture are constructed and trained. Table 5 lists the accuracy of detection ANN of different architectures (The structure *i* − *j* − *k* − *o* means a network with *i* input neurons and *o* output neurons, and *j* and *k* stands for the number of neurons in hidden layers ). It is clear to see from the table that the proposed method has considerable accuracy for collision detection, and an average accuracy of 0.95 can be obtained with 1 or 2 hidden layers.


All the vibration modal data of actual collision are used for the training and testing of positioning neural networks. Considering the geometric layout of the STR6-05 robot, the moveable structure is divided as two parts, i.e., one part is link 3, and the other contains link 4, link 5 and link 6, as shown in Figure 10. Table 6 lists the accuracies of different layers of positioning networks. In addition, the BP networks with 2 hidden layers are suitable for the positioning task.


**Table 6.** Accuracy of collision positioning.

As the vibration modal features of collision direction are mainly reflected by the acceleration signals nearby the collision point, the input of direction identification network is determined by the output of positioning network. Table 7 lists the accuracy of network of 3 architectures. The training and testing samples come from the vibration data of collisions on link 4, link 5 and link 6, which is near to the acceleration sensors on end-effector. It is obvious that the BP network with 2 hidden layers can obtain comparative stabilization accuracy.

**Table 7.** Accuracy of direction identification.


Considering all the results comprehensively, the proposed method can be utilised for the detection, positioning and identification of collisions with considerable accuracy. By analyzing, we find that the positioning error and identification error is mainly derived from boundary samples, that is the collisions near joint 3. One more accelerometer located beside joint 3 may be used for the enhancement of accuracy.

By analysing the training and test procedure, we can find that an experiment of 300 collision samples is enough for the artificial network training for any 6 − *do f* manipulator with different sensor placement scheme, and it can be accomplished in half or one hour with well designed scenarios and test procedures.

### *4.3. Rapid Prototyping System Design*

In order to realize the online test of the proposed method, the computation complexity and detection time of the related algorithms should be estimated. As discussed above, the main calculation consumption includes three parts, namely, the fast Fourier transform of vibration acceleration data, the estimation of natural frequency and modal, and the node outputs update of neural network.

It can be seen from Section 2.2 that the collision vibration frequency of the robot is mainly between 50 Hz and 1500 Hz, and we set the sampling rate to 3200 samples/second. On the other hand, in order to ensure the spectral characteristics have a sufficient resolution, we select a sliding window of 0.1 s for each cycle, and the number of sampling points *N* = 320 for Fourier transform, as shown in Figure 11. The overlap amount of the sliding window is 50 percent, that is, the calculation cycle of the proposed algorithm is 0.05 s. Within each computing period, Butterfly fast Fourier transform is adopted, and the total number of real multiplications required by FFT of *N* points is 2*N* ∗ *log*2(*N*), and the total number of real numbers added is 2*N* ∗ *log*2(*N*) for each vibration signal.

**Figure 11.** Sliding window for online vibration detection test.

For the natural frequency and modal estimation algorithm, the main computational cost of the algorithm is the sorting algorithm of spectrum amplitude, and the computational complexity of the algorithm is *O*(*N* ∗ *log*2(*N*)) [26], that means, the computational cost of the frequency and modal estimation algorithm can be ignored relative to the fast Fourier transform.

For the status update of the neural network, each neuron includes multiple addition operations, multiplication operations and exponential operations, and the exponential operation can be converted into a number of addition and multiplication operations. For the collision detection neural network with typical structure of 27-6-4-1, the addition calculation times of a single update is 278, and the total number of multiplication is 98. For the positioning and direction identification neural network with one or two hidden layers, the calculation cost is of the same magnitude.

Therefore, the main computational cost of this method is derived from the fast Fourier transform. In addition, within a single update cycle (50 ms), the addition and multiplication times of the algorithm is about 22,000 and 21,000 respectively. The total computation consumption of this method is about the same size as the traditional state observer method with Newton-Euler Function, as discussed in [27].

Based on the above analysis, we use Simulink/Data Acquisition Toolbox to develop a rapid prototyping system, as shown in Figure 12. The vibration data within each sliding window is buffered, and then used for FFT, modal estimated, and neural network status updates. The system is a slower-than-real-time system due to data caching and operating system. It shows that the collision detection time is about 0.1 s, and the isolation and identification time will be a little longer.

**Figure 12.** A Rapid Prototyping System.

### **5. Conclusions and Discussion**

In this work, we present a model independent method for robot collision detection, positioning and identification. The vibration signals are analyzed and used for the construction and training of BP neural network. The test results of the experiment confirm that it is possible to build a monitoring algorithm with considerable accuracy. The conclusions may be summarized as follows:


**Author Contributions:** F.M. conceived the original ideas, designed and conducted all the experiments, and subsequently drafted the manuscript. G.W. contributed to the construction of the experiment platform and writing-review. N.L. provided supervision to the project.

**Funding:** This research was funded by the Fundamental Research Fund for the Central University (No. 11618403) and the Guangdong Natural Science Foundation (No. 2018030310482).

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

### **References**


c 2019 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 (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Integration of Terrestrial Laser Scanning and NURBS Modeling for the Deformation Monitoring of an Earth-Rock Dam**

### **Hao Xu 1, Haibo Li 2, Xingguo Yang 2, Shunchao Qi <sup>1</sup> and Jiawen Zhou 1,\***


Received: 9 November 2018; Accepted: 19 December 2018; Published: 21 December 2018

**Abstract:** A complete picture of the deformation characteristics (distribution and evolution) of the geotechnical infrastructures serves as superior information for understanding their potential instability mechanism. How to monitor more completely and accurately the deformation of these infrastructures (either artificial or natural) in the field expediently and roundly remains a scientific topic. The conventional deformation monitoring methods are mostly carried out at a limited number of discrete points and cannot acquire the deformation data of the whole structure. In this paper, a new monitoring methodology of dam deformation and associated results interpretation is presented by taking the advantages of the terrestrial laser scanning (TLS), which, in contrast with most of the conventional methods, is capable of capturing the geometric information at a huge amount of points over an object in a relatively fast manner. By employing the non-uniform rational B-splines (NURBS) technology, the high spatial resolution models of the monitored geotechnical objects can be created with sufficient accuracy based on these point cloud data obtained from application of the TLS. Finally, the characteristics of deformation, to which the geotechnical infrastructures have been subjected, are interpreted more completely according to the models created based on a series of consecutive monitoring exercises at different times. The present methodology is applied to the Changheba earth-rock dam, which allows the visualization of deformation over the entire dam during different periods. Results from analysis of the surface deformation distribution show that the surface deformations in the middle are generally larger than those on both sides near the bank, and the deformations increase with the increase of the elevations. The results from the present application highlight that the adhibition of the TLS and NURBS technology permits a better understanding of deformation behavior of geotechnical objects of large size in the field.

**Keywords:** earth-rock dam; 3D visualization; deformation monitoring; terrestrial laser scanning (TLS); NURBS

### **1. Introduction**

The deformation distribution and evolution are present as important indications of the instability of large artificial and natural structures such as tunnels, bridges, and landslide [1–4]. Characterizing the in-situ deformation behavior of huge/important structures helps to understand the underlying mechanism of sliding and to predict the possibilities of catastrophic collapse. Thus, periodic monitoring of deformation of the structure has been an important task over the whole lifecycle of infrastructures of great importance, e.g., for ensuring their safety during construction as well as for their post-construction maintenances. Take the example of a dam: many monitoring instruments have been adopted with

proper methodologies to appraise the condition and safety of the dam over the past decades [5–7]. Using traditional methods are possible for monitoring over a relatively large area. However, one of the common features of all these approaches e.g. general geodesy method and GPS static method, is that they can only provide point-wise information. For the dams with large scale, the number of the monitoring points is usually rather limited in terms of sufficiently characterizing the deformation characteristics. Synthetic Aperture Radar (SAR) is a relatively new technology, whose variant, GB-SAR, is used for the deformation measurement and monitoring of dam [8–10]. It needs to be recognized that SAR offers high sensitivity to small displacements. The most advanced GB-SAR is even capable of providing millimeter precision. Conversely, the main shortcoming of SAR is that it can only allow for deformation detection along the sensors-target line of sight. Furthermore, the data processing and analysis tools of SAR are quite complex, which makes it difficult to be applied into practice.

In contrast, the terrestrial laser scanning has a definite advantage of providing point clouds composed of millions of points with high accuracy and high spatial resolution, and the point clouds can be used to detect vertical and horizontal deformation through relatively simple process. Nevertheless, the point cloud data acquired from the laser scanner cannot be processed in the same way as that for the data from the traditional methods, since the laser pulse emitted from the machine is not necessarily aimed at the same place of the target object in different scanning operations [11]. Additionally, the point density of the target object varies in different measurement campaigns, which making direct comparison between points cloud obtained sequentially is not favorable.

Hence, how to extract information of the dam deformation from the terrestrial laser scanning (TLS) data is a challenging and important task, which is the focus of many recent works carried out on the deformation detection by the application of TSL [12–16]. A simple means of extracting the deformation is to compare the point-cloud gathered in different epochs [17–19]. This, however, can only show the general deformation tendency exhibited by the target object, which implies that it only allows for a qualitative analysis of the deformation rather than a quantitative analysis. The traditional method for surface reconstruction is point cloud gridding. A typical type of gridding algorithm, called the Delaunay triangulation algorithm, is widely used in modern programs, like Rapid Form and Riscan Pro. Due to its flexibility and adaptability, the Delaunay triangulation algorithm is widely accepted for reconstructing the surface of the objects with irregular shape, such as land relief and mechanical parts [20–22]. Unfortunately, the accuracy of the triangulated irregular network model is not satisfactory, making the subsequent deformation detection difficult to accomplish. Besides, the triangulation algorithm is rather disappointing in the aspect of noise reduction.

In this paper, a detailed deformation measurement method, based on combination of TLS and NURBS technologies, is presented. The method is established on the basis of the quadrangular surface domain parameters for deformation monitoring. By specifying the quadrangular surface domain parameters, the four corner points of the NURBS fitting surface are exactly positioned at the four vertexes of the control mesh. Thus, the precision of the NURBS fitting surface reaches up to a couple of millimeters and is higher than the traditional NURBS surface. The performance of this present method in the field is illustrated by an application to a selected hydraulic structure, the Changheba Dam, located in Southwest China. To acquire the detailed and accurate 3D data rapidly and efficiently, the TLS technology is first applied to take inventories. To fully take advantage of the point cloud data, the digital surface model of the Changheba Dam with high accuracy and high spatial resolution is established by using the NURBS surface modeling technology. Through the comparison between multi-temporal models, the deformation characteristics of the Changheba Dam in different stages are analyzed. The results presented in this study illustrate the applicability of the present methodology to the precise deformation monitoring over large regions.

### **2. Background**

### *2.1. The Changheba Dam*

Located about 360 km southwest of Chengdu, the capital city of Sichuan province in Southwest China, the Changheba Hydropower Station retains the Dadu River over a basin area of 56,648 km2 (Figure 1a,b). The hydro-technical system, with the primary purpose of exploiting the hydropower potential of the Dadu River, consists of an earth-rock dam, a spillway system and a water diversion and power generation facilities. The construction of the whole project was completed in April, 2018. With the maximum height of 240.0 m, the Changheba dam is the highest earth-rock dam built with gravel and soil core wall in China at present (Figure 1c). The length and width of the dam crest are 502.8 m and 16.6 m, respectively. Its gross filling volume is about 3.42 × 107 m3, the largest volume of core wall earth-rock dam constructed in China. The dam can be divided into eight different zones according to the filling materials. The distinct mechanical properties of the filling materials make the deformation of the dam very intricate.

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

(**c**)

**Figure 1.** Location and layout of the dam at Changheba Hydropower Station: (**a**) and (**b**) location of the Changheba Hydropower Station, (**c**) layout of the dam at Changheba Hydropower Station.

### *2.2. Geodetic Network*

A geodetic network in the nearby of the earth-rock dam has been constructed since the commencement of the project in the year of 2011, which consists of nine control points materialized by observation monuments with a forced centering device. Of the nine control points, five are located on the left bank of the Dadu River, with other four on the right bank, and Figure 2 shows the layout of the geodetic network. The maximum and minimal side lengths of the network are about 1200 m and 210 m, respectively. Datum points TN03, TN04 and TN06 are equipped with inverted plumb line system. Thanks to the good conservation of these observation monuments, the geodetic network has been working well since the beginning. Two measurements of the geodetic network have been performed in October 2016, April 2017, respectively, by means of a Leica TM30 GeoRobot (provided by Leica Geosystems AG, Heerbrugg, Switzeland). The observation results of the TN04 and TN06 based on the inverted plumb line system shows the two points of the control net are stable. Thus, they are used as starting point for the classic free network adjustment calculation. The mean square error of points is less than +2 mm using the classic free network adjustment. Of the two measurements, the mean square error of the weakest point is 1.61 mm and 1.76 mm, respectively. And the average relative mean square error of the weakest side is 1/283,000 and 1/258,000, respectively. The results of the measurements and adjustment calculation indicate that the geodetic network is fine enough for the coordinate measurement of the scanner.

**Figure 2.** Layout of the geodetic network.

### **3. Methods**

### *3.1. Terrestrial Laser Scanning*

Terrestrial laser scanning is a quite new surveying method in geodesy. The scanner illuminates the target object with pulsed laser and records the returned pulse of the laser. It calculates the distance from itself to the target object surface automatically by timing the round-trip time of one pulse of the laser. Utilizing the time-of-flight (TOF) technology, the TLS is able to record dense point-clouds over an extremely short period of time. The modern laser scanner provided by the manufacturer shows a significant higher speed of data acquisition, compared to the conventional surveying instruments like total station [23]. Moreover, the TLS can also record the intensity of the reflected pulsed laser and RGB color data of the target object.

Let the origin of the Cartesian (rectangular) coordinate system coincide with the center of the laser scanner, where X and Y are two axes perpendicular to each other lying in the horizontal plane, and the Z axis is oriented upwards and perpendicular to the horizontal plane (Figure 3). Then, the coordinates of the laser point in the there-dimensional (3-D) space can be calculated from Equations (1) and (2) [24].

$$S = c \times \left(\frac{TOF}{2}\right) \tag{1}$$

$$\begin{cases} \ x = S \cdot \cos \theta \cdot \cos \alpha \\\ y = S \cdot \cos \theta \cdot \sin \alpha \\\ z = S \cdot \sin \theta \end{cases} \tag{2}$$

where *c* is the speed of light; *TOF* is the time of flight of the laser pulse; *S* is the distance from the scanner to the reflecting surface; *θ* is the angle between the line *OP* and the *XY* plane; and α is the angle between *X* axis and the orthogonal projection of the *OP* onto the plane *XY*.

**Figure 3.** Measuring principle of the terrestrial laser scanning (TLS).

In this paper, a pulse-based scanner, Riegl VZ-400 (provided by RIEGL Laser Measurement Systems, Horn, Austria), is applied to survey the study area. The scanner is capable of measuring the distance ranging from 1.5 m to 600 m. Due to the echo digitization and online waveform processing technique, the highest angular resolution that this instrument can achieve is 0.0005◦, and its efficient measurement rate is up to 120,000 points per second. Besides, it offers a wide field of view up to 100◦ vertical and 360◦ horizontal. As such, this equipment is suitable for the data acquisition in the vested condition.

### *3.2. Data Acquisitions and Preprocessing*

Measurement campaigns were performed after the filling of the dam. Due to the great size of the dam, it is impossible to acquire the point cloud data of the entire dam surface with only one scan. In addition, a relatively larger number of scans can increase the density of target point clouds, making sure that the overlap between different scans is enough for the alignment of different data sets.

Here, both the particular geologic of this case and the requirement of data processing, nine scans were performed at nine different stations distributed in the surrounding of the dam for the first time in October 2016, when the project began to impound. For each scan, the scanner was placed on the observation station. Then the coordinates of the scanner were measured based on the geodetic network. In this way, it was possible to get more precise coordinates of the scanner and reduce alignment error. It took approximately nine minutes to complete each 360◦ scan. The whole data set acquired in the first measurement campaign comprises one hundred and thirty-six million points in total, which is treated as the reference point cloud. The second measurement campaign was performed in April 2017 in a similar way as that in the first campaign, and both measurements were referred to the geodetic coordinate system.

These two sets of TLS data were preprocessed using the RISCAN PRO software. First, the geodetic coordinates of each scanner were extracted from the total station based on the geodetic network. By simply adding the instrument height, the geodetic coordinates of the scanner center were acquired. Then the geodetic coordinate data was input into the RISCAN PRO software for the raw scans by the backsighting orientation. In this way, the location of each scan center was determined. After the input of the coordinates, the orientation of each scan remained random. Therefore, a manual modification of the orientations was carried out to make the orientations in space correct. The second process was multi station adjustment. This process was performed using the iterative closet point (ICP) algorithm proposed by Besl and McKay [25]. In the first place scans were aligned pair by pair by means of lowering the "search distance" parameter from meters to some centimeters step by step. After a few adjustments, the alignment of the pair of scans led to an optimal rota-translation alignment matrix. Following this procedure a global alignment was employed for the whole scans to obtain a best fit alignment [26]. Both steps were based on the ICP algorithm and the latter step was executed for the purpose of distributing the residual registration error more homogeneously across the scans [25]. The global alignment of the nine scans resulted in the overall standard deviation of 0.0017 m, which meant the average distance between two closest points. The alignment error is mainly owing to the point spacing of the datasets in the overlap area and the point measurement error in practice. By increasing the point density in the overlap area and reducing the distance between the scanner and the overlap area, the standard deviation can be significantly lowered. In consideration of the high point density of the point clouds, this alignment error is primarily caused by point measurement error which can be reduced by fitting technique. Thus, the alignment error is acceptable for the deformation monitoring of the dam. Meanwhile, the four circular reflector targets were used as tie point for the calibration of the alignment. Similarly, the geodetic coordinates of its center points were measured by total station. Moreover, the geodetic coordinates of the center points were extracted from the scanning data. Then there were two sets of geodetic coordinates of the center points. By taking the geodetic coordinates acquired from the total station as datum, the mean square error of the four center points was calculated out. The mean square error of the points was no more than 0.0013 m, which proved that the alignment of the scans was successful. The alignment of the partial scans leaded to a single point cloud data set consisting of all scanned points.

After alignment of all nine scans, there were some unneeded objects such as dust and vegetation in the point cloud. Via automatic and manual operation, the data was run through terrain filter to remove such points. Eventually, a single point cloud colored with RGB information obtained from the calibrated camera on top of the terrestrial laser scanner came into being, whose partial view is shown in Figure 4 below.

**Figure 4.** Aligned 3D model of the Changhheba dam (rendered view).

After the alignment of the point cloud, it comes to the stage of the triangulation. The basic input data that can be processed for NURBS surface reconstruction is the triangulated irregular network (TIN) rather than the point cloud. Hence, the TIN should be constructed based on the point cloud first, i.e., triangulation. The triangulation of points is also called "tessellation". By using the Delaunay triangulation algorithm, each set of three closest points in the point cloud are connected to form a triangle, resulting in a non-overlapping triangulation as a whole. However, the initial triangulation generated automatically by the Delaunay triangulation algorithm usually has deficiencies, like acute angled triangles and holes. With automatic analysis, inconsistencies are detected and then repaired through automatic or manual operation. The final operation carried out to refine the TIN is the smoothing. The purpose of smoothing operation is to reduce the noise and therefore to ensure both the quality and accuracy of NURBS surface that will be constructed later. The downstream face of the Changheba dam is made of masonry; the dam surface is flat macroscopically. Thus, the smoothing operation is of central significance to the surface reconstruction.

The point clouds for the Changheba dam are converted into a polygon object with 3.48 billion triangles in total (Figure 5). Only points representing the dam surface are employed and converted into TIN, and other points beyond the survey area are not processed in this section.

**Figure 5.** Triangulated irregular network (TIN) model of the Changheba dam.

The previously-constructed TIN is exported and then processed into a NURBS fitting surface model of the Changheba dam. With high accuracy and high spatial resolution, this NURBS surface model can be applied to detect the deformation of the dam surface for further study.

### *3.3. NURBS Modeling*

NURBS is the abbreviation of non-uniform rational B-splines, in which, the non-uniform means that the spacing of the knots is uneven, the Rational implies that the control point can be weighted, and the B-spline represents that B-spline is used as basis function. NURBS is a piecewise rational vector polynomial function advocated by Versprille in 1975 [27]. It can generate and represent arbitrary curves and surfaces better than other functions like radial basis function [28], which can be in either standard shapes or free-form shapes. The NURBS is widely used in computer graphics and the CAD/CAM industry, due to its great flexibility and precision. Besides, fitting surface is a useful means of reducing noise. Thus, NURBS is a useful tool to build surface model of natural land relief [29].

*Sensors* **2019**, *19*, 22

A NURBS surface *P*(*u, v*)={*x*(*u, v*), *y*(*u, v*), *z*(*u, v*)} is a piecewise rational surface defined by Equation (3):

$$P(u,v) = \frac{\sum\_{i=0}^{m} \sum\_{j=0}^{n} w\_{i,j} N\_{i,k}(u) N\_{j,l}(v) d\_{i,j}}{\sum\_{i=0}^{m} \sum\_{j=0}^{n} w\_{i,j} N\_{i,k}(u) N\_{j,l}(v)} \tag{3}$$

where, the *di,j* (*i* = 0, 1, ... , m; *j* = 0, 1, ... , n) are the control points representing a topological mesh, *wi,j* is the so-called weights, and the *Ni,k*(*u*), *Nj,l*(*v*) are the normalized B-spline basis functions defined on the non-periodic knot. For instance, the mathematical expression of the *Ni,k*(*u*) can be defined recursively as follow,

$$N\_{i,0}(u) = \begin{cases} \quad 1, \quad u\_i \le u \le u\_{i+1} \\ \quad 0, \quad otherwise \end{cases} \tag{4}$$

$$N\_{i,k}(\boldsymbol{u}) = \frac{\boldsymbol{u} - \boldsymbol{u}\_{i}}{\boldsymbol{u}\_{i+k} - \boldsymbol{u}\_{i}} N\_{i,k-1}(\boldsymbol{u}) + \frac{\boldsymbol{u}\_{i+k+1} - \boldsymbol{u}}{\boldsymbol{u}\_{i+k+1} - \boldsymbol{u}\_{i+1}} N\_{i+1,k-1}(\boldsymbol{u}) \tag{5}$$

To generate a NURBS surface, three groups of parameters, control points *di,j*, weight factors *wi,j* and the knot vectors *U* and *V* must be determined. On the grounds of the NURBS theory, the dam surface is constructed from the TIN following the four procedures: panel demarcation, surface patch insertion, grid generation, and NURBS surface construction. The four procedures are executed in sequence to construct the NURBS surface, as explained in the following.

Data segmentation must be carried out to precisely construct the NURBS model of dam surface primarily. In this step, the characteristic lines are extracted from the TIN model. The surface curvature detection technique is adhibited for data segmentation. After setting the proper curvature level, the curvature detecting is performed and the contour lines are placed in the areas of curvature. In the meantime, Markers are highlighted on the contour lines. There are two different kinds of markers. A red marker indicates a corner, while a yellow marker represents a non-corner point of inflection. The contour lines and the boundary lines are used as panel demarcation lines dividing the TIN model into forty-three panels (Figure 6a). If a high curvature level was specified, fewer panels would be acquired. Then the precision the NURBS fitting surface would be affected. But more panels require more processing power and lower the operability. In this study, both the precision requirement and operability, the forty-three panels are proper for the NURBS fitting surface. The demarcated panels will be the containers for the surface patches in the next step.

On the basis of the result of the panel demarcation, the patch boundary structure is generated, which means that four-sided patches are inserted into each panel (Figure 6b). A surface patch is a four-sided subdivision of a panel that is approximately equilateral. Based on the auto estimate technique, the target patch count is automatically calculated depending on the size and the smoothness of the panels. As the panels are very irregular and the patches are relatively regular, the adaptive insertion technique is applied in order to acquire fairly uniform patches, which contributes to the precision of the NURBSB model. Using the adaptive surface patch insertion technique, the shape of the surface patches is determined. The model consists of eight hundred and thirty-nine surface patches. From Figure 6b, it can be seen that the patches that comprise a panel differ in size, and the shapes of most patches are close to the rectangle. The inserted patches will be the containers for the grids in the following procedure.

After the surface patches are inserted, a further process is grid generation. A grid is a quadrangular mesh constructed in every patch. The mesh is made up of a forty-by-forty set of rectangles, which means that each panel consists of one thousand and six hundred grids. The quadrangular mesh density can be specified by the user. Needless to say, a finer grid produces greater precision in the eventual NURBS surface. The grid generation process creates an ordered u-v grid in every patch on the model (Figure 6c).

**Figure 6.** Non-uniform rational B-splines (NURBS) surface construction procedures: (**a**) panel demarcation; (**b**) surface patch insertion; (**c**) grid generation; (**d**) NURBS surface construction.

For the free-form surface of the dam, the exact NURBS surface can be ultimately generated on the TIN using the surface fitting technique. The NURBS surface of the dam is shown in Figure 6d. The consequent high accuracy and high spatial resolution surface model captures very well the morphology and geometry features of downstream face of the Changheba dam, which can be then inquired into the deformation characteristics of the dam [30,31].

### *3.4. Deformation Measurement by Shortest Distance (SD) Comparison*

Deformation can be detected by making geometrical comparison between multi-temporal surface models. In this study, the shortest distance algorithm is applied to the deformation measurement. It is noted that the normal directions of the two surfaces should be generally consistent before the comparison. The algorithm can still work even when the surface normal vector is biased. For each point *i* (*xi.tes*, *yi.tes*, *zi.tes*) *<sup>T</sup>* in the test model, the algorithm searches for its nearest corresponding point *j* (*xj.ref*, *yj.ref*, *zj.ref*) *<sup>T</sup>* in the reference model and computes the SD vector, *Vi*, that starts in point *i* and ends in point *j* (Equation (6)) [32].

$$V\_i = \begin{pmatrix} \Delta \mathbf{x}\_i, \Delta \mathbf{y}\_i, \Delta \mathbf{z}\_i \end{pmatrix}^T = \begin{pmatrix} \mathbf{x}\_{i.ts\star} y\_{i.ts\star} z\_{i.ts\star} \end{pmatrix}^T - \begin{pmatrix} \mathbf{x}\_{j.ref} \ y\_{j.ref} \ z\_{j.ref} \end{pmatrix}^T \tag{6}$$

For the NURBS fitting surface, the central points of grid are used to calculate the shortest distance. The calculated SD vectors do not necessarily represent the real displacements of the object. They are the shortest distances from the point in the test model to the nearest neighbor point in the reference model

mathematically (Figure 7). Yet, the shortest distances (SDs) are useful for deformation measurement, since they allow for the detection of the vertical, horizontal, and oblique distances. In this paper, the sign conventions are defined as follow. Positive SDs indicated that the points in the test model are in front of or above those in the reference model. They may not appear to be continuous and can be interpreted as material stack owing to human activity. Negative SDs signified that the detecting part is behind or below the reference dataset, which are related to the vertical settlement and subsidence.

**Figure 7.** Conceptual comparison of the profiles.

To analyze the deformation of the structure, it is necessary to choose the proper coordinate system [33]. Herein, the main deformation of the dam surface taking place during the operation period is the settlement along the vertical direction. Hence, the Z direction is set as the vertical direction. And the Y direction is set in the direction parallel to the river. The X direction is set in the direction perpendicular to the river. Then, the shortest distance calculated from multi-temporal models can represent the deformation. The calculated data is useful for the dam deformation monitoring [32], with which the deformation of the downstream face can be analyzed in detail.

### **4. Results**

### *4.1. Deformation Analysis*

It is well known that the earth-rock dams are usually subjected to deformation over a relatively long time after the filling. Obviously, the deformation distribution is an essential indicator of the potential instability of the dam. Thus, periodic monitoring of the dam deformation during operation period is of great importance. The deformation of the surface is always larger than the interior due to the effect of displacement cumulation. Therefore, monitoring the surface deformation is an effective way of determining the serviceability of the dam. By introducing the TLS and NURBS technology, making comparison between different multi-temporal scans provides the plentiful and vital information for analysis of the deformation distribution of the dam.

For the Changheba dam studied herein, the NURBS model in October 2016 was treated as reference model, while the NURBS model in April 2017 was regarded as the test model. The reference model was subtracted from the test model. In this way, more than eight hundred thousand central points of grid in the test model were calculated and color-coded mappings of the differences were generated, showing the deformation distribution over the monitored time interval during the operation period.

The deformation distribution of the Changheba dam from October 2016 to April 2017 is shown in Figure 8. On account of the high water level in April 2017, the point cloud data of the upstream face was not obtained, at which the deformation cannot be calculated with only point cloud data in October 2016. As a result, only the deformation of the downstream face is presented. Negative changes representing

the dam settlement are shown in cold colors. As shown, the dam has experienced settlements continuously after the filling, which can be interpreted as consolidation. The differential deformation is relatively significant. On top of the dam, the maximum deformation value is −0.0976 m, while the deformation values near the dam toe are close to zero. Along the stream direction, the deformation values of the dam surface get gradually smaller from the top to the toe. In the direction perpendicular to the river, the envelope of deformation exhibits a counter-arch shape (Figure 8). That is to say, the deformation in the middle of the dam is larger than that at the sides of the dam at the same elevation. The deformation of the zigzag road is much larger than its near regions due to the intense human activities.

**Figure 8.** Deformation distribution of the Changheba dam. SDs: shortest distances.

There are three main regions with positive changes in Figure 8, marked as R1, R2, and R3. R3 was in the construct platform in elevation 1551 m in virtue of the construction activities. The other two parts R1 and R2 are located near two abutments. The point cloud data here is missing in the second measurement campaign in April 2017. The hole was filled manually in the process of NURBS model construction. Thus, the calculated result does not indicate the real changes happening in the site.

### *4.2. Deformation Mechanism*

Figure 9 shows some deformation distribution along the cross sections A-A, B-B (identified in Figure 8). As shown in Figures 8 and 9, this dam continued to undergo deformation after the completion of filling of the dam, with the maximum deformation value up to virtually 90 mm in the two cross sections. In the meanwhile, the deformations in the middle at both cross-sections are normally larger than those near the banks, which caused by the effects of constraint from the banks on both sides. The envelopes of deformations of both cross sections are asymmetric, reflecting the fact that the earth-rock dam has been constructed in an asymmetric canyon.

**Figure 9.** The deformation distribution along the two cross-sections.

It is worth mentioning that the deformations in cross section A-A are much smaller compared to those in cross section B-B. One of the reasons is that the lower part of the dam has been subjected to settlement over a relatively longer time than the upper part of the dam. The lower part of dam started to deform right after its filling, at this time the upper part has not been filled. Another reason is that, after filling of the upper part, additional settlement of the lower part will be induced due to the addition of gravity load on the top. Thus, the relative smaller deformation of lower part is the consequence of combination of a longer period of consolidation and higher stress conditions.

Figure 10 show some deformation distribution along the longitudinal sections I-I, II-II, III-III, and IV-IV, as marked in Figure 8. The magnitudes of deformations increase nonlinearly with the increasing elevation. The relatively larger deformations of the upper part demonstrate that the inner part of the dam has experienced a rapid consolidation during the investigation period. It can also been seen that the deformations of longitudinal sections II-II, III-III are always larger than those at the other two sections I-I, IV-IV, at the same elevation, which means the deformations in the middle are generally bigger than those on the sides at the same elevation. This corresponds with the conclusion drawn from Figure 9. In the meanwhile, there exist several abnormal regions circled in red in Figure 10, which represent the positive changes taking place on the dam surface. This phenomenon results from the human activities and the data missing.

**Figure 10.** The deformation distribution along longitudinal sections.

### **5. Discussions**

In this paper, a new deformation monitoring method for large structure is formulated. The field experiments are carried out to test its feasibility. In practical works, reducing error is of vital importance for deformation monitoring. According to the different stages of the methodology, the possible errors can be divided into three parts: instrument error, alignment error and modeling error. Many researches have been done on the instrument error and alignment error for the past decades [34,35]. Here, the modeling error of the NURBS technology is analyzed. The NURBS surface model constructed is a kind of fitting surface. There exist some errors in the constructed surfaces, and the error distribution is shown in Figure 11 by making comparison between the original point cloud and the NURBS fitting surface. The overall error of the fitting surface is relatively small, and the distribution of the errors is quite uniform. Errors in most parts are within ±0.002 m. The relatively larger errors are concentrated mainly in the places with large curvatures. One example is in the area of the construction platform including D1 and D2. Since the construction platform is not the focus of this current study, where the NURBS fitting surface is constructed roughly in order to reduce the size of the NURBS surface file. Nevertheless, the fitting surface may not be subtle enough to represent surface of objects with the large curvatures where large errors are likely to be induced. Another area with relatively large errors is at the two sides of the zigzag road form dam crest to dam toe. The NURBS surface is about 0.006 m higher than the original points inside the zigzag road. On the contrary, the NURBS surface is approximately 0.006 m lower than the original points outside the zigzag road. This is also due to the large curvatures of the steps and the drainage ditch. As the larger errors don't exist in the main study reaches, the constructed NURBS surface can be applied to the analysis of the deformation distribution. Error analysis above proves that the NURBS fitting surface is useful way of surface reconstruction.

**Figure 11.** Error distribution in the NURBS modeling process.

### **6. Conclusions**

The deformation distribution has a significant influence on the stability and safety of large artificial and natural structures. By taking full advantage of the TLS and NURBS technology, a new methodology of deformation monitoring and analysis is presented. Using the earth-rock dam as an example, the point cloud gained by TLS allows the detection of small deformation for its high accuracy and capability of target acquisition. The NURBS model based on the point cloud is characterized by high precision and high spatial resolution. Eventually, the holistic deformation distribution of the downstream face is shown in the cloud chart. The deformation monitoring achieved great success with the proposed methodology. The TLS has proved to be a valid solution for the deformation acquisition in three-dimensional space. The NURBS modeling technology is capable of dealing with a huge number of points and making use of them. In comparison with the traditional monitoring, the methodology that integrates TLS and NURBS technologies permits a better grasp of the deformation distribution of the large structures.

The millimeter level measurement requires that the data acquisition is performed with great patience and plenty of time. In the future works, data acquisition in the field should be optimized to reduce the field working time. Besides, the data processing is also time-consuming and complicated. Several procedures require personal experience and expertise. These are several potential areas in terms of improving the performance of the presented methodology, in which further research will be conducted.

**Author Contributions:** H.X. and H.L. conceived the original ideas, designed the experiments, processed the data and drafted the manuscript. S.Q. and J.Z. helped to write and edit the manuscript. J.Z. and X.Y. supervised the project.

**Funding:** This research was funded by the Key Project of the Power Construction Corporation of China (ZDZX-5) and the Youth Science and Technology Fund of Sichuan Province (2016JQ0011).

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

### **References**


© 2018 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 (http://creativecommons.org/licenses/by/4.0/).
