Abstract
Playing an important role in electromechanical systems, hydraulic servo system is crucial to mechanical systems like engineering machinery, metallurgical machinery, ships, and other equipment. Fault diagnosis based on monitoring and sensory signals plays an important role in avoiding catastrophic accidents and enormous economic losses. This study presents a fault diagnosis scheme for hydraulic servo system using compressed random subspace based ReliefF (CRSR) method. From the point of view of feature selection, the scheme utilizes CRSR method to determine the most stable feature combination that contains the most adequate information simultaneously. Based on the feature selection structure of ReliefF, CRSR employs feature integration rules in the compressed domain. Meanwhile, CRSR substitutes information entropy and fuzzy membership for traditional distance measurement index. The proposed CRSR method is able to enhance the robustness of the feature information against interference while selecting the feature combination with balanced information expressing ability. To demonstrate the effectiveness of the proposed CRSR method, a hydraulic servo system joint simulation model is constructed by HyPneu and Simulink, and three fault modes are injected to generate the validation data.
1. Introduction
Hydraulic servo system plays a crucial role in electromechanical systems, like engineering machinery, metallurgical machinery, ships, and other equipment. Failures of hydraulic servo system caused by severe and complex conditions may lead to catastrophic accidents and enormous economic losses. Fault diagnosis based on monitoring and sensory signals is able to classify the current state of complex systems, which plays a key role in performance evaluation [1]. Feature set extracted from signals is an important index to reflect the fault mechanism and performance evolution laws. The quality of feature set plays a key role in improving the generalization ability of fault identification [2]. The common feature extraction methods are time-frequency index extraction, wavelet analysis, Hilbert transform, Duffing oscillator, and so on. Despite their respective applicable conditions and limitations, those methods are able to mine the health characteristics of the system from multiaspect [3, 4]. Same as machine learning, features extracted from images, speeches, and other signals often have certain correlations and hidden mutual influences. Information expressed by a single feature is usually inadequate, which can be greatly improved when the single feature is aggregated with others [5]. Similarly, due to the nonlinearity, instability, and nonconformity of complex electromechanical systems, the expression of the information on individual feature is often one-sided. Thus, a new challenge is how to utilize those features more effectively and efficiently, in other words, how to obtain the feature set that expresses the information sufficiently by eliminating the redundant and negatively correlated features [6–9].
To tackle the challenge mentioned above, on the premise of existing feature extraction techniques, feature processing techniques including feature selection and dimension reduction have gradually become an important research focus. Both feature selection and dimension reduction can reduce the scale of feature set by obtaining a set of principal variables. Such techniques often use a variety of feature extraction methods to integrate the features into a comprehensive representation of the signal [10]. For the purpose of enhancing the expressing ability of core information on multiclass feature sets, spatial transformation or importance measurement methods are used [11]. Such methods are able to reduce redundancy existing in features and improve learning efficiency while retaining the performance advantages. The data transformation may be linear, such as principal component analysis (PCA). But many nonlinear dimension reduction techniques also exist. The common feature processing techniques include linear dimensionality reduction methods (FDA, LPP, etc.), kernel function based dimension reduction methods (KPCA, KFDA, etc.), manifold based dimension reduction methods (Isomap, LLE, MDS, etc.), filter method based feature selection methods (Relief, Focus, information gain, CBFS, etc.), wrapper method based feature selection methods (genetic algorithm, distribution estimation, differential algorithm, etc.), and embedded method based feature selection methods (SVN-RFE, RF, etc.) [12–14].
For the nonlinear signals of complex electromechanical systems, although the dimension reduction methods can reduce the scale of input features for fault diagnosis, they change the basic attributes of the feature set. Such situation makes it difficult to give a clear understanding of the obtained feature subset. Meanwhile, existing feature selection methods sort the importance degree of the features according to the independent feature evaluation result. They ignore the interaction among features, which would lead to information loss in processing the data of electromechanical systems [15]. Aiming at the shortcomings of feature selection and dimension reduction techniques, such as low expandability, unclear evaluation indexes, and strong tendency, this study proposes compressed random space based ReliefF (CRSR) method. Based on ReliefF method, CRSR introduces ensemble strategy on feature level based on compressed random space. Furthermore, CRSR optimizes the objective function using information entropy and fuzzy theory. The main contributions are as follows:(i)This study analyzes the feasibility of the ReliefF based feature selection architecture for the fault diagnosis of complex electromechanical system. Meanwhile, the basic mechanism of measuring the contribution of the features based on ReliefF is also demonstrated.(ii)By converting the assessment process of ReliefF, which takes the entire feature set as object, into the construction and ensemble process of feature subspace, CRSR can improve the global optimization ability of ReliefF.(iii)Considering ReliefF based feature selection as a problem of maximizing the distance, CRSR replaces the traditional spatial distance with fuzzy membership degree, which is able to obtain a robust and steady objective function.
This paper is structured as follows: In Section 2, ReliefF method based feature selection structure is introduced. Then feature integration method based on compressed random subspace is described. Objective function optimization method based on information entropy and fuzzy theory is also introduced. Section 3 presents the overall diagnosis procedure for hydraulic servo systems, the details of construction and fault injection of the hydraulic servo system, and analysis and comparison on feature selection results using the proposed CRSR method which are discusses numerically. Section 4 concludes the paper.
2. Related Theories
2.1. ReliefF Method Based Feature Selection Structure
ReliefF is the extension of Relief method by estimating probabilities more reliably, which is able to handle incomplete and multiclass data sets while the complexity remains the same [16–19]. By calculating the distances between the sample distributions, ReliefF can obtain the correlation weight coefficients of the features which is similar to Relief.
For a specific feature from the feature set, if its difference in same class is much smaller than that in different class, it is considered that this feature contributes to class discrimination [20]. Given a sample set with instances, each sample, , has m-dimensional features. Meanwhile, the samples in only belong to two classes which are tagged as . The difference between each two samples ( and ) on feature is defined aswhere the attribute of feature is discrete value. If the attribute of feature is continuous value,
Features extracted from condition monitoring data of electromechanical system mostly are continuous data. Meanwhile and represent the maximum value and minimum value of the entire sample, respectively. The closest same-class instance of sample is called “near-hit (NH),” and the closest different-class instance of sample is called “near-miss (NM).” Meanwhile, the weight of feature is denoted as , and is updated bywhere the initial value of is 0, and is the last value of .
For reducing the randomness in feature evaluation, the whole process needs to repeat times to obtain the average value being the final weight. Although Relief method is very efficient in estimating the quality of attributes, it cannot deal with incomplete data and is limited to two-class problem. Thus ReliefF method is utilized in the paper to deal with the multiclass classification and regression problems for continuous data.
For a multiclass classification problem, assuming that the samples in belong to multiple classes and the tags for are , ReliefF updates on sample by taking near hits (NHS) and near misses (NMS) into consideration, which is different from Relief. Similarly, the weight of feature , which is , can be updated throughwhere represents the ratio of the entire samples in class to all the heterogeneous samples in . Furthermore, ReliefF method equalizes the differentiation of NHs and calculates the average differences between and other classes on feature to evaluate the classification ability of the samples nearby.
2.2. Feature Integration Method Based on Compressed Random Subspace
The purpose of CRSR method is finding the balance between the differences and the correlations of features. Specifically, in the premise of fully mining the correlations of features using ReliefF method, CRSR method is applied to make each feature subset keep a certain degree of difference. Based on random subspace and feature sorting strategy, the feature sets with higher contribution can be obtained in various feature combinations [21, 22]. and are denoted as two random subspaces, so the difference, denoted as , can be calculated aswhere symbol denotes the dimension of random subspaces.
The right side of (5) obtains the noncoincident features of and . is plus one when there is an unrepeated feature . The average difference between two random subspaces on all the features is defined aswhere is the dimension after feature ordering compression. It can be simplified as . Concretely, according to sorting strategy, RS12, the probability of the first strongly related features being selected is , and the relatively poorly related features being selected are . The average difference of feature evaluation based on compressed random subspace method, which is , can be denoted as follows:
Equation (7) shows that the ranking result can balance the difference of ReliefF to determine the dominant features that are crucial for classification, which would improve the feature selection efficiency.
Based on compressed random subspace method, this study proposes redundancy analysis method from statistics to reduce the redundancy of feature. The features are checked in pairs using redundancy analysis method [23, 24]. Firstly, two sets of feature vectors, and , are selected randomly from the feature set obtained by ReliefF. Then, the selected feature vectors are regarded as independent variable and dependent variable separately. The covariance matrixes, denoted as and , can be calculated, respectively, aswhere represents mathematical expectation of vector . Then, the correlation coefficient of and , denoted as , can be formulated aswhere denotes the covariance of and . If the correlation coefficient is greater than a presetting threshold, only the one with larger weight from and will be added to the final selected feature set. It is noticed that redundancy analysis based on matrix transformation focuses on the correlation between features instead of the similarity of data. Thus, CRSR can reduce the influence of numerical confusion existing in feature selection. Furthermore, compared with traditional methods based on data similarity, CRSR is able to obtain higher confidence level.
2.3. Objective Function Optimization Method Based on Information Entropy and Fuzzy Theory
From the aspect of maximizing the distance, Relief method can be seen as a distance optimization algorithm using feature weighting method [25]. Under this condition, the optimization objective function, denoted as , can be described as [26]where denotes the distance of the th sample. Based on (4), can be converted as
For complex electromechanical system signals, two problems occur when using (10) as the optimization objective function of CRSR. One problem is that the objective function concentrates the weight on one or some of the features, which leads to the result that assessment value of the remaining features tends to 0. Meanwhile, (10) regards the samples with stochastic volatility and noise similar to the normal sample. Another problem is the lack of consideration on the influence from the quality of samples on the feature selection process.
Aiming at the first problem, the information entropy theory is proposed based on compressed random subspace method, which combines the maximization of entropy together with maximization of distance to reduce the over tendency problem of the existing ReliefF method. After adding a sample estimation factor , the optimization objective function is denoted as
Supposing that and follow probability distribution, Shannon entropy is used to adjust the sample distribution, as shown below:
Aiming at the second problem, the fuzzy membership degree is chosen to replace the traditional nearest neighbor distance. The fuzzy membership degree has the advantage of being insensitive to sample fluctuations and noise and the ability of updating adaptively while changing the feature weights [27]. In a sample space of same class, the fuzzy membership degree of the feature , denoted as , can be calculated as [28, 29]where is the sample set of same class, and is the fuzzy difference between feature and feature , as shown below:
Then the same fuzzy distance, denoted as , is calculated as
Similarly, the fuzzy difference and the fuzzy membership degree of the heterogeneous sample sets of , which are denoted as and separately, can be calculated as follows:
Therefore, the heterogeneous fuzzy distance, denoted as , is calculated as
Based on the fuzzy distances obtained by (16) and (18), updated (11) can be formulated as follows:
Based on the formulas above, the objective function, which is , can be denoted aswhere the first item on the right side is the maximum distance of ReliefF, which is meant to determine the feature set that contributes most to classification. Meanwhile, from the aspect of entropy maximization, the second and the third items denote the sample evaluation operator and Shannon entropy of feature weight, respectively, which are used to avoid the over tendency problem of objective function. and are the balance coefficients for adjusting the differences between features. When the maximum point of objective function is achieved, the constraint condition of the sample evaluation operator is defined below [30]:where denotes the partial derivative of .
In the process of feature selection for mechanical system, the information entropy and fuzzy theory based optimization objective function ensures that the evaluation process of each subspace using ReliefF is adaptive and robust. Such advantage provides a new thought for feature processing of complex monitoring signal, in other words, under the premise of maintaining the high calculating efficiency of ReliefF method, reducing the bias and redundancy caused by methodological defects and external disturbances [31].
3. Method for Hydraulic Servo System Fault Diagnosis Based on CRSR
The CRSR based fault diagnosis method for hydraulic servo system consists of following successive steps: first, the average, standard deviation, skewness, and wavelet singular entropy features are extracted to form a feature matrix as the input of ReliefF model. Second, the initial contribution of the features sampled randomly is measured by calculating the inner-class distance and between-class distance, and the features of high contribution and features of low contribution are determined for fault classification. The th iteration operation is based on the result of the previous iteration, and the iteration stops as long as reaches the preset threshold. Third, based on the sorting result of the second step, the compressed evaluation of the features is realized using supervised sampling method in current iteration. Finally, keeping the iteration running until certain terms is satisfied to acquire the difference value as the criterion of feature selection. The detailed process of feature selection method using CRSR for hydraulic servo pump is shown in Figure 1.
Compared with the traditional feature processing methods [32–35], the CRSR method is designed to meet the robustness and accuracy requirements of fault diagnosis with smaller feature set and lower resource consumption. This study demonstrates the advantages of the proposed method by extracting a variety of categories of features from the simulation data of hydraulic servo system and optimizes the feature set extracted by CRSR method to verify the feasibility in feature selection technique.
3.1. Description of the Simulation Environment of Hydraulic Servo System
The data used in this case is generated in HyPneu and Simulink joint simulation environment. The joint simulation process of hydraulic servo system can be divided into mechanical part and control part considering the characteristics of the hydraulic components [36–38]. The mechanical hydraulic physical part, including hydraulic pumps, servo valves, and actuators, is modeled using component library in HyPneu. The failure is realized by adding the modules of fault injection for the relevant components. The dynamic control part, including feedback sensors and electronic amplifiers, is established using the relevant model in Simulink. Meanwhile, the control of the hydraulic servo system and the simulation of fault injection can be realized by transmitting the signal data through relevant interface files. Thus, the model architectures of hydraulic servo system constructed using HyPneu and Simulink are shown in Figures 2 and 3, respectively.
As shown in Figure 3, the Simulink model includes input signals, comparison elements, control elements, amplifying elements, and HyPneu module. The gain parameter of the electronic amplifying element is 80, and the control parameter of PID is set as proportional parameter of 1500, integral parameter of 0, and differential parameter of 5. To simulate the actual operating environment of hydraulic system, random noises in the range of −0.01 to 0.01 are added on the output of the HyPneu module.
The process of joint simulation based on HyPneu and Simulink is as follows. First, the actuator receives the feedback signal converted by servo valve and right after that drives the load to do reciprocated motion. Then, the real-time displacement information of the load is collected by the sensor and transferred to the control circuit. Finally, the input signal and the amplified displacement signal are compared and inputted into the servo valve to build the negative feedback control logic of closed loop. The joint simulation model has a good ability to match the hydraulic system and to realize the fault injection of the hydraulic system effectively.
3.2. Fault Injection for Hydraulic Servo System and Multidimensional Feature Extraction
This study simulates four kinds of state of the hydraulic system including normal state, electronic amplifier fault, sensor constant deviation fault, and hydraulic pump wear fault. The detailed fault injection scheme is shown in Table 1.
The input of all fault modes is sinusoidal signal with amplitude being 1 and frequency being . As the expression of the input signal is , the sampling frequency is 100 Hz and the sampling time is 70 s. The input-output relationships of the hydraulic servo system under circumstances of normal state and three fault states are collected, and every signal of the corresponding state contains 7000 points. The original output signals are shown in Figure 4.
(a) Normal state
(b) Amplifier fault
(c) Sensor fault
(d) Hydraulic pump wear fault
As mentioned above, the input of the proposed CRSR method is a feature matrix, and each column of the matrix represents a feature vector extracted from the original signal. As the prerequisite of feature selection, multidimensional features should be extracted for all of the signals. It should be noted that the cycle of hydraulic servo system signal collected for the feature selection should be an integer multiple of the input waveform cycle. The cycle of the input signal is 0.5 seconds and one input-output cycle contains 50 points, which means the object points of the feature extraction should be a geometric multiple of 50. In addition, in consideration of the fact that pump wear fault is a gradual degradation process and the degree of fault in the early phase may be weak, this study takes the last 4000 points of the original signal to extract the feature matrix. Moreover, a sampling window of 500-point length is used (with 5-point interval) to obtain 700 feature vectors.
Furthermore, this study adds 15 dB Gaussian white noise to all the original signals to validate the performance of CRSR method. The details of the extracted features are shown in Table 2.
As can be seen from Table 2, the quantitative range of different fault features fluctuates much, which could lead to the bias problem in the feature selection process of ReliefF. Therefore, this study normalizes the features to the range of zero to one. The visualization of the feature vectors of different states is shown in Figure 5.
As is shown in Figure 5, data confusion exists in the features of different fault modes, which makes it difficult to select satisfied feature subsets through artificial observation. Thus, CRSR method is designed to select the feature combination with higher contribution to the fault diagnosis adaptively.
3.3. Analysis and Comparison on Feature Selection Result
The input of CRSR feature selection model is a matrix acquired by regularizing the 9 sets of the obtained feature vectors. One of the advantages of CRSR is the ability to calculate the distances of either same class or different classes adaptively and iteratively. Based on the ReliefF feature selection method, the corresponding distribution vectors of feature weight can be obtained.
During the iteration process of ReliefF, feature selection constraint is established using compressed stochastic subspace method, which optimizes the performance and efficiency of CRSR. The parameters are as follows: number of iterations: 30; feature base of the subspace: 6. And the threshold condition of feature selection is that the difference between a feature being selected as “one with high contribution” and “one with low contribution” is greater than 10. With – representing the feature parameters in Figure 5, the statistical result of the features being selected as “one with high contribution” and “one with lo contribution” is shown in Table 3.
It can be seen from Table 3 that the feature combination selected by CRSR method contains average (), standard deviation (), mean square root (), crest factor (), and the maximum amplitude of FFT (). Although wavelet singular entropy () was selected as “one with high contribution” for 14 times during the iteration, it was selected as “one with low contribution” for 6 times. Such result reveals that wavelet singular entropy is able to contribute to classification sometimes, but the ensemble performance is not stable. Even in some feature subspace matrixes, it exerts a negative impact on the classification. The final statistical result of feature weights assessed by CRSR is shown in Figure 6.
As is shown in Figure 6, if the weight basis for ReliefF is set to 0.75, the feature combination selected by CRSR method will append features (total energy of wavelet) and (wavelet singular entropy). However, according to the statistical results from the CRSR model, the contribution of these two features in the integration process is unstable. In particular, as for the total energy of wavelet, the reason why its average weight coefficient is high is that the weight coefficient reaches 0.92 and 0.88 in the 8th and 22nd iterations, respectively, while in the remaining subspaces its contributions are lower than the average slightly. Such situation reveals that the total energy of wavelet could not meet the requirement of stability in feature selection. The result shows that the CRSR method is more reasonable than the traditional ReliefF, which reflects the advantages of ensemble learning in generalization.
Based on the feature selection result using CRSR, redundancy analysis for the feature sequence is carried out in this study. The threshold of correlation coefficient is set to 0.8, and the correlation matrix of feature is shown in Table 4.
It can be seen from Table 4 that the correlation coefficient of peak factor () and the maximum amplitude of FFT transform () are greater than the preset threshold, which indicates that redundancy exists in the fault diagnosis information provided by them. Therefore, only the maximum amplitude of FFT transform () is retained. In summary, the features selected by CRSR are average (), standard deviation (), mean square root (), and the maximum amplitude of FFT ().
The purpose of introducing CRSR method is improving the performance of fault diagnosis. In other words, for a classifier, the diagnostic performance using selected feature set should not be less than that using the original feature set. The comparative models used in this study contain the classical ReliefF algorithm, the Mean Impact Value (MIV) algorithm, the Locally Linear Embedding (LLE) algorithm, and Kernel Principal Component Analysis (KPCA). In this study, the feature sets selected by different feature selection methods are used as inputs for the classifier based on Radial Basis Function (RBF) neural network. The ratio of training samples to test samples is 50%. The fault diagnosis accuracy of the hydraulic servo system is obtained by using the 10-fold cross validation method as shown in Table 5.
In Table 5, the first row () represents the final feature set selected by CRSR; the second row () represents the feature set selected by CRSR before the redundancy analysis; the third row () represents the feature set selected by the classical ReliefF method; the fourth line () represents the feature set selected by MIV; the fifth row () and the sixth row () are the high-dimensional feature sequences obtained by dimension reduction using KPCA and LLE, respectively, where the dimension reduction target is set to be same as the feature numbers determined by CRSR. The last row () is the collection of original features. Taking the feature combinations mentioned above as input, the average and variance of fault diagnosis accuracy of hydraulic servo system calculated are shown in Figure 7, respectively.
According to the RBF network fault diagnostic results from Table 5 and Figure 7, the following can be obtained.
Compared with using the original feature set directly, the feature selection methods are able to improve the precision of fault diagnosis based on RBF classifier, as the averages of the 10-fold cross validation from to are 94.405%, 94.217%, 93.022%, 93.045%, 93.077%, 94.069%, and 92.983%, respectively. The result indicates that the feature selection process plays a positive role in the classification task using high-dimensional data.
Compared with the original ReliefF method, the introduction of compressed random subspace method eliminates those features with unstable contribution in the process of subspace integration and improves the performance and efficiency of feature selection. Eventually, higher fault diagnosis accuracy can be achieved with fewer input features.
The variances calculated for (with redundancy analysis) and (without redundancy analysis) are 0.29 and 0.37, respectively. It can be seen that although the average diagnostic accuracy of them is close, the redundancy analysis can optimize the information repetition among the features and reduce the computational resource consumption.
Compared with the MIV-based feature selection method, CRSR has a great advantage in the diagnosis performance, which indicates that it has higher confidence in feature selection. Compared with KPCA and LLE, although the latter two (especially the LLE reduction algorithm) achieve high diagnostic accuracy as well, the variance of the 10-fold cross validation results for KPCA and LLE is higher than that of CRSR (0.29 and 0.31, resp.), which indicates that a certain degree of volatility exists in the results obtained by dimension reduction algorithms. In addition, due to the lack of clear interpretability and the existence of ambiguity in optimal target parameter setting, it is difficult to determine the best reduction method among the dimension reduction algorithms. Thus adaptive distance metrics based CRSR method has better applicability.
To further illustrate the rationality of CRSR method, the diagnostic results were tested by rank-based nonparametric Kruskal-Wallis test, which is based on fault diagnosis accuracy as a basic indicator in Table 5. The Kruskal-Wallis test was used to determine whether the different features sets have significant influence on the diagnostic performance. The significance threshold is set to 0.05, and the test results returned by Kruskal-Wallis function are shown in Figure 8. The value returned by Kruskal-Wallis function is 0.00002 which is far less than 0.05. It indicates that the different combinations of features make significant impact on the fault diagnosis performance for hydraulic servo system.
Moreover, another improvement of CRSR is that the target optimization function is constructed based on the information entropy and the fuzzy theory. Compared with the distance measurement method used in the classical ReliefF model, the robustness of the feature selection process is improved in a complex environment. To illustrate the improvement of CRSR mentioned above, this study analyzes the original output signal and the output signal with 15 dB noise. Meanwhile, PCA is used to project two sets of the selected feature sequences into three-dimensional space. Furthermore, the K-means method is used to cluster the feature sets of different fault modes. The results are shown in Figure 9.
(a) CRSR cluster result without noise
(b) CRSR cluster result with 15 dB noise
(c) ReliefF cluster result without noise
(d) ReliefF cluster result with 15 dB noise
It can be seen from Figure 9 that the feature qualities selected by CRSR and ReliefF are higher before the noise is added, and the feature sets corresponding to each failure mode can be clearly distinguished by K-means method. After adding 15 dB of noise, due to the influence of external disturbances, the data points of the amplifier gain fault and the sensor constant deviation fault are obviously confused in the clustering result of ReliefF, which indicates that the feature differentiation declines. However, in the result of CRSR method, spatial distribution boundaries between the four fault modes are still identifiable, which proves that the proposed objective function optimization method has good robustness against noise. Thus, CRSR method is a promising technique for feature selection and subsequent fault diagnosis of hydraulic servo systems.
4. Conclusion
This study presents a fault diagnosis scheme for hydraulic servo system using compressed random subspace based ReliefF (CRSR) method. Based on the feature selection structure of ReliefF, the proposed CRSR method employs feature integration rules in the compressed domain and substitutes information entropy and fuzzy membership for traditional distance measure index. The advantage of the proposed method lies in the ability of determining the feature set with the better generalization performance and the less resource consumption. As a data-driven method, CRSR could be practical and flexible in engineering.
To demonstrate the effectiveness of the proposed CRSR method, validation data of three fault modes is generated through a hydraulic servo system joint simulation model. Comparing with existing feature reduction and feature selection methods, the result indicates that the feature selection process plays a positive role in the classification task using high-dimensional data, and CRSR based fault diagnosis method has higher average accuracy and smaller variance. Meanwhile, the compressed random subspace method can eliminate those features with unstable contribution in the process of subspace integration and improve the performance and efficiency of feature selection. Besides, due to the robustness and stability of the information entropy and fuzzy theory based objective function optimization, the result shows that CRSR method is more suitable for fault diagnosis problem under noisy conditions.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was supported by the National Natural Science Foundation of China [Grant nos. 51605014, 51105019, and 51575021], the Technology Foundation Program of National Defense [Grant no. Z132013B002], and the Fundamental Research Funds for the Central Universities [Grant no. YWF-16-BJ-J-18].