Research Article | Open Access
Jianlong Zhao, Yongchao Zhang, Qingguang Chen, "Rolling Bearing Fault Feature Extraction Based on Adaptive Tunable Q-Factor Wavelet Transform and Spectral Kurtosis", Shock and Vibration, vol. 2020, Article ID 8875179, 19 pages, 2020. https://doi.org/10.1155/2020/8875179
Rolling Bearing Fault Feature Extraction Based on Adaptive Tunable Q-Factor Wavelet Transform and Spectral Kurtosis
The fault feature of the rolling bearing is difficult to extract when weak fault occurs and interference exists. The tunable Q-factor wavelet transform (TQWT) can effectively extract the weak fault characteristic of the rolling bearing, but the manual selection of the Q-factor affects the decomposition result and only using TQWT presents interference. Aiming at the above problems, an adaptive tunable Q-factor wavelet transform (ATQWT) and spectral kurtosis (SK) method is proposed in this paper. Firstly, the method applies particle swarm optimization (PSO) to seek the optimized Q-factor for avoiding manual selection, which uses the kurtosis value of the transient impact component as the particle fitness function. The rolling bearing fault signal is decomposed into continuous oscillation component and transient impact component containing fault feature by the optimized Q-factor. Then, due to the presence of interference in the decomposition result of ATQWT, the SK analysis of the transient impact component is used to determine the frequency band of periodic impact component characterizing fault feature by fast kurtogram. Finally, the band-pass filter established through the above frequency band is employed to filter the interference in the transient impact component. Simulation and experimental results indicate that the ATQWT can highlight the periodic impact component characterizing rolling bearing fault feature, and the SK can filter interference in the transient impact component, which improves feature extraction effect and has great significance to enhance fault diagnosis accuracy of the rolling bearing. Compared with EEMD-TQWT and TQWT-SK, the fault feature extracted by the proposed method is prominent and effective.
Rolling bearing, a harsh working environment and high incidence of failure, is a critical part of rotating machinery. When fault occurs, the measured vibration signal is the combination effect from its structure, machining and assembly errors, movement and force of other parts, and operating failure, etc. [1–5]. It contains abundant feature information and the excitation with different fault leads to distinct response. Therefore, it is the most important task in the rolling bearing fault diagnosis, extracting the fault feature from the measured vibration signal .
At the moment of rolling bearing failure, the continuous contact between the failure surface and other surfaces generates a series of impact pulse force with broadband . The resonance is excited on account of impact pulse force covered high-frequency natural vibration of the rolling bearing, and its attenuation time is much smaller than the interval between the impulse forces, thus generating a series of periodic impact components . In the early stage of rolling bearing fault, the periodic impact component that characterizes the rolling bearing fault is weak, has interference, and is difficult to extract, so a lot of research has been conducted on the fault feature extraction method of rolling bearing in recent years. Time-frequency analysis methods are used to extract fault feature of rolling bearing such as short-time Fourier transform (STFT), Winger-Ville distribution (WVD), wavelet transform (WT), and empirical mode decomposition (EMD) . The size of the time-frequency window in STFT does not change and lacks adaptability. When the WVD is used for multi-component signals, it is easy to generate cross term and difficult to extract feature frequency . The WT needs to select the basis function and decomposition layer in advance, and the selection has a greater impact on the final decomposition result. EMD has a certain mode confusion .
On account of the problems existing in the above feature extraction methods, Selesnick proposed the tunable Q-factor wavelet transform (TQWT), which can decompose the vibration signal into the continuous oscillation component and the transient impact component according to the difference of pre-set Q-factor [12, 13]. Wang used the TQWT to decompose the IMF of the maximum kurtosis value obtained by EEMD, and the low resonance component was carried out by envelope spectrum analysis to realize the rolling bearing fault diagnosis . Li decomposed the vibration signal of rolling bearing through ICD and extracted the fault feature of rolling bearing using TQWT . Bharath utilized TQWT to extract rolling bearing fault feature and realized the classification of rolling bearing fault via classification algorithm . Consequently, the TQWT can extract the rolling bearing fault feature. Spectral kurtosis (SK) can detect the non-Gaussian component in the signal and provide its statistical indicator in the frequency position [17–19]. Abbas founded the resonance band of the defective rolling bearing through SK and achieved the rolling bearing fault identification . Wan applied SK to select the resonance frequency band and extract the fault feature using the envelope demodulation method . After the above analysis, we can find that the periodic impact component characterizing the rolling bearing fault feature is not obviously directly obtained by the TQWT, has a lot of interferences, and is influenced by the selected Q-factor. The SK can locate the frequency band of the periodic impact component owing to its non-Gaussian property, but the filtered rolling bearing fault signal has many interferences only using the above frequency band. Therefore, the fault feature extraction method needs to be further improved.
In view of the advantages of TQWT and SK in extracting rolling bearing fault feature, this paper proposes a fault feature extraction method based on adaptive tunable Q-factor wavelet transform (ATQWT) and SK. The high and low factors are optimized by particle swarm optimization (PSO), and the optimized Q-factor is used to perform tunable Q-factor wavelet transform to highlight the periodic impact component in the transient impact component; then the SK analysis of the transient impact component determines the frequency range of the periodic impact component. The filter established by the above frequency band can filter out the interference component in the transient impact component, highlight the fault feature, and improve the accuracy of fault diagnosis of rolling bearing.
2. Fault Feature Extraction Method of Rolling Bearing
2.1. Tunable Q-Factor Wavelet Transform
To adapt to the feature of the analyzed signal, the Q-factor of the wavelet transform should be selected according to its property. It should have a high Q-factor in analyzing continuous oscillation signal; otherwise it should own low Q-factor in the transient impact signal, but its Q-factor is difficult to adjust. The TQWT proposed by Selesnick can select the Q-factor according to the characteristic of the analyzed signal and possess perfect reconstruction and moderate overcompleteness and is completely discrete .
The decomposition and reconstruction filter bank of the TQWT is shown in Figure 1 and a total of 4 sub-band signals are obtained, where is the sampling frequency, is the low-pass scaling factor, and is the high-pass scaling factor. is the decomposition layer, and the high-frequency and low-frequency sub-band signal obtained by the j-level filter bank are and , respectively, and the sampling rates are and , respectively. is the input of the j + 1 level filter bank. and are low-pass and high-pass filter response function. and are the complex conjugate of and .
In order to ensure that the wavelet transform is not excessively redundant and has perfect reconstruction property, the scale factor should satisfywhere is the low-pass scaling factor and is the high-pass scaling factor.
For perfect reconstruction, the low-pass filter and the high-pass filter should be constructed bywhere and are low-pass and high-pass filter response function, respectively, and is the -periodic power complementary function.
The center frequency and bandwidth of the decomposition layer j can be calculated by equation (5), and the scale factor can be calculated by equation (6):where is the center frequency, BW is the bandwidth, Q is the quality factor, r is the oversampling rate, and is the sampling frequency of the input signal.
The specified Q-factor should be chosen subject to . Setting leads to a wavelet transform for which the wavelet resembles the second derivative of a Gaussian. Higher values of Q lead to more oscillatory wavelets. The specified oversampling rate r must be strictly greater than 1. If r is close to unity, then the transition bands of and will be relatively narrow and the time-domain response will not be well localized. For , the wavelet will resemble the sinc wavelet. In order to avoid this issue, it is sufficient to select .
Different correspond to different decomposition layer J. It is determined by the maximum decomposition layer. The maximum decomposition layer can be obtained bywhere N is the length of the input signal and is the maximum decomposition layer .
The rolling bearing fault signal is composed of continuous oscillation component and transient impact component, and it is difficult to separate them owing to their overlap frequency band or existing interference. Morphological component analysis can effectively separate the continuous oscillation component and the transient impact component according to their oscillation property; the signal x is supposed to be a linear combination of high resonance component and low resonance component with morphological difference; then, x is expressed as
In order to achieve the separation of and employing morphological component analysis, and are assumed as overcomplete dictionaries for signals and , respectively, and can only be represented sparsely by and can only be represented sparsely by , and and are the coefficient vectors of sparse representation; then, the solution to the sparse decomposition coefficient of the signal is transformed into an optimization problem as in
The solution of the resonance sparse decomposition coefficient is transformed into optimization problem of equation (10) through morphological component analysis:
By iteratively solving the decomposition coefficient and corresponding to the minimum value of equation (10), the high resonance component and low resonance component obtained by the resonance sparse decomposition are, respectively, expressed as
2.2. Adaptive Tunable Q-Factor Wavelet Transform
The Q-factor of the TQWT should be selected according to the oscillation property of the decomposed signal. According to the oscillation property of the signal, the high Q-factor of TQWT is greater than or equal to 3, and the low Q-factor is between 1 and 3. The greater high Q-quality does not mean the better decomposition effect. The search space is relatively small, and long search time will cause casualty for fault diagnosis. Moreover, the PSO is simple to implement and has fast convergence speed and strong practicability. Considering the above factors, the PSO was selected among several algorithms such as genetic algorithm and ant colony algorithm, etc. The PSO is used to optimize the Q-factor of the tunable Q-factor wavelet transform. The kurtosis value of the low resonance component is used as a particle fitness function. After continuous iteration, the Q-factor adapted to the decomposed signal is finally obtained. The PSO is an algorithm for solving practical optimization problem obtained by studying the migration and foraging behavior of birds. To randomly initialize a group of particles without mass and volume, each particle is regarded as a potential solution of the optimization problem and continuously flying in the search space at a certain speed, and the particle can find the optimal solution through iteratively updating fitness .
Suppose there are m particles in the n-dimensional target search space, the position and velocity of the i-th particle are and , and the i-th particle individual extremum is , and the global extremum of the entire population is . Due to the iterative update of the current particle velocity and position according to its velocity and position, individual extremum, and global extremum, the k + 1 iteration of the i-th particle is carried out according towhere is the inertial weight, m is the total number of particles, n is the dimension of the solution space, and are the learning factors, and and are the random numbers between (0, 1).
2.3. Spectral Kurtosis
Spectral kurtosis can detect Gaussian components contained in the signal  and can provide statistical indicators of its position in the frequency domain. It calculates the kurtosis of time domain data corresponding to each spectral line to determine the position of the transient impact components . The SK is defined as follows.
Spectral moment is a very valuable statistical indicator for measuring nonstationary processes, which is defined aswhere is the spectral moment and is the time-varying function.
The calculation of the spectral moment needs to be repeated and the alternative method of spectral moment averages the instantaneous moment along the time axis, as shown inwhere represents the time average operator and T is time.
Under the conditions of stationarity and ergodicity of the complex envelope , equation (16) can be proved:
The spectral cumulant is very sensitive to nonstationary processes. The non-Gaussian process’s order spectral cumulant is nonzero, and the fourth-order spectral cumulant is defined aswhere is the fourth-order spectral cumulant.
The fourth-order spectral accumulation is larger when the signal much more deviates the Gaussian. The energy-normalized fourth-order spectral accumulation will give the peak measure of the probability density function of the signal at frequency f, also the so-called SK. It is shown inwhere is the spectral kurtosis.
2.4. Fault Feature Extraction Based on Adapted TQWT and SK
Using kurtosis as the objective function to optimize the Q-factor through PSO, the tunable Q-factor wavelet transform can highlight the periodic impact component that characterizes the fault of the rolling bearing in the low resonance component. The SK can locate the frequency band of the periodic impact component. The filter established by the frequency band can filter out interference component and highlight fault feature. The rolling bearing fault type can be determined by the envelope spectrum analysis. The algorithm flow chart is shown in Figure 2. The specific algorithm is as follows:(1)Initializing , m, Loopcount, , the fitness function is established through the kurtosis of the low resonance component of the tunable quality wavelet transform(2)Setting the redundancy factor and and the decomposition layer and , the Q-factor of the tunable factor wavelet transform is optimized by PSO(3)The rolling bearing fault signal is decomposed by tunable quality wavelet transform using optimized Q-factor, and the high and low resonance components are acquired(4)The fast kurtogram gained by SK of the low resonance component finds out central frequency and bandwidth BW of the periodic impact component characterizing rolling bearing fault feature(5)The finite impulse response filter established by the central frequency and bandwidth BW filter the low resonance component(6)The envelope spectrum of the filtered signal is gained by Hilbert envelope spectrum analysis, and the rolling bearing fault type is identified by the amplitude distribution in envelope spectrum
3. Simulation Analysis
In order to verify the effectiveness of the ATQWT and SK in rolling bearing fault feature extraction, a simulation signal x is simulated to the rolling bearing fault vibration signal, where with a frequency of 100 Hz simulates the periodic impact component characterizing rolling bearing fault feature, and the sinusoidal signals with a frequency of 30 Hz and 2358 Hz express the rotational frequency and harmonic, and the Gaussian white noise is used as the background noise with the intensity of −9db:
Setting the sampling frequency to 10 kHz and the sampling points to 2048, and x are shown in Figures 3(a) and 3(b), respectively. The periodic impact component representing the fault feature of rolling bearing was submerged, and it was impossible to judge whether the rolling bearings failed in Figure 3(b). As shown in Figure 4, the envelope spectrum of is most prominent in 100 Hz and its frequency doubling in Figure 4(a), which is consistent with the rolling bearing fault characteristic frequency, but the envelope spectrum of x has no noteworthy spectral lines in Figure 4(b), so it is not possible to directly apply Hilbert envelope to extract rolling bearing fault feature frequency.
The Q-factor optimized by PSO is , and the rolling bearing fault signal is carried out ATQWT at parameters , and the result is shown in Figure 5. In Figure 5(b), the periodic impact component in the low resonance component is not prominent and has interference, and its Hilbert envelope spectrum is presented in Figure 6. The 100 Hz and its frequency doubling emerge in Figure 6, but its amplitude is small and possesses interference in the frequency domain. Therefore, the periodic impact component in the low resonance component acquired by ATQWT is not prominent, has interference, and affects the accuracy of the fault diagnosis of the rolling bearing.
The low resonance component is performed by the SK analysis, and the fast kurtogram is shown in Figure 7. Its peak value presents in , and the above frequency band includes the 2000 Hz high-frequency natural vibration characterizing faults rolling bearing fault, so the SK can accurately locate the frequency band of the periodic impact component characterizing rolling bearing fault. Building a finite impulse response filter with order = 30, the low resonance component is filtered by the above , and the time domain and envelope spectrum of filter signal is shown in Figure 8. The periodic impact component in Figure 8(a) is prominent and the fault characteristic frequency of 100 Hz and its multiple frequency in Figure 8(b) have same outcome. Hence, the ATQWT and SK can extract the fault feature of rolling bearing, filter out the interference, and highlight the fault feature.
In order to highlight the advantages of the proposed method in rolling bearing fault feature extraction, the decomposition results of the TQWT-SK and EEMD-TQWT methods are shown in Figures 9 and 10; it can be found that the 100 Hz and its multiple frequency are not prominent and have very small amplitude, so the effect of the rolling bearing fault feature extraction proposed in this paper is significantly better than TWQT-SK and EEMD-TQWT.
4. Experimental Verification
To prove the effectiveness and feasibility of ATQWT and SK method in rolling bearing fault feature extraction, the bearing-rotor testbed is established and is presented in Figure 11. The fault rolling bearing is made by electrical discharge machining and then using it to drill 0.2 mm diameter holes simulating the inner and outer race fault, respectively. The bearing designation is 6308, and the specific parameters are shown in Table 1. The inner and outer race fault signals of the rolling bearing are, respectively, collected in the rotational frequency , the sampling frequency of 10 kHz, and the sampling points of 2048. According to the calculation formula of the fault characteristic frequency and the data in the Table 1, the theoretical fault characteristic frequency of the inner and outer race of the rolling bearing is and .
4.1. Rolling Bearing Inner Race Fault
Figure 12 shows the fault signal and the envelope spectrum of the inner race fault; there is no periodic impact component in Figure 12(a) and no prominent spectral line of in Figure 12(b); it does not judge the rolling bearing fault type. The Q-factors and optimized by PSO are used to conduct ATQWT of the inner race fault signal and other parameters are set to and . The decomposition result is shown in Figure 13, and there exist impact component and a lot of interference in the low resonance component. Figure 14 is not prominent in and its multiple frequency; therefore, directly performing ATQWT on the signal cannot accurately extract the periodic impact component of the inner race fault.
Figure 15 shows the fast kurtogram of low resonance component, its peak present in the central frequency of 3750 Hz and bandwidth of 2500 Hz; it indicates that the periodic impact component is contained in the frequency band. The low resonance component is filtered by a finite impulse response filter with order = 30, and its result is presented in Figure 16, the periodic impact component in Figure 16(a) and the spectral lines of in Figure 16(b) are prominent, and the fault feature of inner race was extracted successfully.
The decomposition results of the TQWT-SK and EEMD-TQWT in inner race fault are shown in Figures 17 and 18; it can be found that the time domain waveform has many periodic impact components, but and its multiple frequency are not prominent and have small amplitude and present a series of modulation frequency. So the feature extraction effect of the inner race fault proposed in this paper is significantly better than TWQT-SK and EEMD-TQWT.
4.2. Rolling Bearing Outer Race Fault
The time domain waveform and envelope spectrum of the outer race fault are shown in Figure 19, the periodic impact component in the time domain waveform is obliterated, and no prominent spectral lines are found in the envelope spectrum, and the failure state cannot be determined. Setting the parameters of performs the ATQWT to obtain high and low resonance components as shown in Figure 20, where Q-factor is obtained by the PSO; there presents periodic impact component and exists a series of interference. The envelope spectrum of the low resonance component in Figure 21 is prominent in ; it affects the accuracy of the rolling bearing fault diagnosis.
The peak value in Figure 22 presents in the central frequency of 3125 Hz and bandwidth of 1250 Hz; it indicates that the periodic impact component is contained in the frequency band. It is presented in Figure 23 that time domain waveform and envelope spectrum of the low resonance component filtered by with order = 30 and the periodic impact component in Figure 23(a) and the spectral lines of in Figure 16(b) are prominent; the fault feature of outer race was extracted successfully.
As present in Figures 24 and 25, the outer race fault decomposition result of TQWT-SK and EEMD-TQWT presents a series of harmonic interference, and the periodic impact periodic component is not prominent and submerged. The frequency domain has frequency modulation and has a small amplitude in outer fault feature frequency and its multiple frequency. Therefore, the feature extraction effect of the outer race fault proposed in this paper is significantly better than TWQT-SK and EEMD-TQWT.
In order to overcome the difficulty of fault feature extraction of the rolling bearing, the selection of Q-quality, and existing interference in TQWT, the ATQWT and SK extraction method is proposed. The method applies PSO to seek the optimized Q-factor of TQWT and the SK analysis is used to determine the frequency band of periodic impact component characterizing fault feature by fast kurtogram. Its effectiveness is verified through simulation and experiment, and the following conclusions are obtained:(1)The Q-factor optimization through PSO can be adapted to the decomposed signal, highlighting the periodic impact component characterizing the rolling bearing fault feature in the low resonance component(2)The SK analysis of the low resonance component can accurately locate the frequency band of the periodic impact component, and the filter established by this band can filter out the interference in the rolling bearing fault signal(3)The method based on ATQWT and SK can extract weak fault feature, improve the feature extraction effect, and enhance the fault diagnosis accuracy of rolling bearing
The data used to support the findings of this study are available from the corresponding author.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
This work was supported by the Shandong Natural Science Foundation (ZR2018MEE036), the National Natural Science Foundation of China (5197050945), and Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD).
- C. Rajeswari, B. Sathiyabhama, S. Devendiran, and K. Manivannan, “Bearing fault diagnosis using wavelet packet transform, hybrid PSO and support vector machine,” Procedia Engineering, vol. 97, pp. 1772–1783, 2014.
- H. Zhao, H. Liu, J. Xu, and W. Deng, “Performance prediction using high-order differential mathematical morphology gradient spectrum entropy and extreme learning machine,” IEEE Transactions on Instrumentation and Measurement, vol. 69, no. 7, pp. 4165–4172, 2020.
- W. Deng, H. Liu, J. Xu, H. Zhao, and Y. Song, “An improved quantum-inspired differential evolution algorithm for deep belief network,” IEEE Transactions on Instrumentation and Measurement, vol. 2020, Article ID 2983233, 1 page, 2020.
- H. Zhao, J. Zheng, W. Deng, and Y. Song, “Semi-supervised broad learning system based on manifold regularization and broad network,” IEEE Transactions on Circuits and Systems I: Regular Papers, vol. 67, no. 3, pp. 983–994, 2020.
- W. Deng, J. Xu, and H. Zhao, “An improved ant colony optimization algorithm based on hybrid strategies for scheduling problem,” IEEE Access, vol. 7, Article ID 2897580, 20292 pages, 2019.
- X. W. Qin, J. L. Guo, X. G. Dong, and Y. Guo, “The fault diagnosis of rolling bearing based on variational mode decomposition and random forest,” Shock and Vibration, vol. 2020, Article ID 1576150, 11 pages, 2020.
- H. B. Lin, F. T. Wu, and G. L He, “Rolling bearing fault diagnosis using impulse feature enhancement and nonconvex regularization,” Mechanical Systems and Signal Processing, vol. 142, Article ID 106790, 2020.
- C. Mishra, A. K. Samantaray, and G. Chakraborty, “Rolling element bearing fault diagnosis under slow speed operation using wavelet de-noising,” Measurement, vol. 103, pp. 77–86, 2017.
- H. M. Liu, L. F. Li, and J. Ma, “Rolling bearing fault diagnosis based on STFT-Deep learning and sound signals,” Shock and Vibration, vol. 2016, Article ID 6127479, 2016.
- B. Vishwash, P. S. Pai, N. S. Sriram, R. Ahmed, H. S. Kumar, and G. S. Vijay, “Multiscale slope feature extraction for gear and bearing fault diagnosis using wavelet transform,” Procedia Materials Science, vol. 5, pp. 1650–1659, 2014.
- X. W. Qin, Q. L. Li, X. G. Dong, and S. Lv, “The fault diagnosis of rolling bearing based on ensemble empirical mode decomposition and random forest,” Shock and Vibration, vol. 2017, Article ID 2623081, 9 pages, 2017.
- I. W. Selesnick, “Resonance-based signal decomposition: a new sparsity-enabled signal analysis method,” Signal Processing, vol. 91, no. 12, pp. 2793–2809, 2011.
- I. W. Selesnick, “Wavelet transform with tunable Q-factor,” IEEE Transactions on Signal Processing, vol. 59, no. 8, pp. 3560–3575, 2011.
- H. C. Wang, J. Chen, and G. M. Dong, “Feature extraction of rolling bearing’s early weak fault based on EEMD and tunable Q-factor wavelet transform,” Mechanical Systems and Signal Process, vol. 48, no. 1-2, pp. 103–119, 2014.
- Y. Li, X. Liang, M. Xu, and W. Huang, “Early fault feature extraction of rolling bearing based on ICD and tunable Q-factor wavelet transform,” Mechanical Systems and Signal Processing, vol. 86, pp. 204–223, 2017.
- I. Bharath, S. Devendiran, D. M. reddy, and A. T. Mathew, “Bearing condition monitoring using tunable Q-factor wavelet transform, spectral features and classification algorithm,” Materials Today: Proceedings, vol. 5, no. 5, pp. 11476–11490, 2018.
- J. Antoni, “The spectral kurtosis: a useful tool for characterising non-stationary signals,” Mechanical Systems and Signal Processing, vol. 20, no. 2, pp. 282–307, 2006.
- Y. Hu, W. Bao, X. Tu, F. Li, and K. Li, “An adaptive spectral kurtosis method and its application to fault detection of rolling element bearings,” IEEE Transactions on Instrumentation and Measurement, vol. 69, no. 3, pp. 739–750, 2020.
- W. Bao, X. Tu, Y. Hu, and F. Li, “Envelope spectrum L-Kurtosis and its application for fault detection of rolling element bearings,” IEEE Transactions on Instrumentation and Measurement, vol. 69, no. 5, pp. 1993–2002, 2020.
- A. R. Bastami and A. Bashari, “Rolling element bearing diagnosis using spectral kurtosis based on optimized impulse response wavelet,” Journal of Vibration & Control, vol. 26, no. 3-4, pp. 175–185, 2019.
- S. Wan, X. Zhang, and L. Dou, “Compound fault diagnosis of bearings using an improved spectral kurtosis by MCDK,” Mathematical Problems in Engineering, vol. 2018, pp. 1–12, 2018.
- J. Kennedy, R. C. Eberhart, and Y. Shi, “The particle swarm,” Swarm Intelligence, Springer, Berlin, Germany, 2001.
- J. Antoni, “Fast computation of the kurtogram for the detection of transient faults,” Mechanical Systems and Signal Processing, vol. 21, no. 1, pp. 108–124, 2007.
Copyright © 2020 Jianlong Zhao et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.