Integration of Wavelet Denoising and HHT Applied to the Analysis of Bridge Dynamic Characteristics

When the dynamic characteristics of a bridge structure are analyzed though Hilbert–Huang transform (HHT), the noise contained in the bridge dynamic monitoring data may seriously affect the performance of the first natural frequency identification. A time-frequency analysis method that integrates wavelet threshold denoising and HHT is proposed to overcome this deficiency. The denoising effect of the experimental analysis on the simulated noisy signals proves the effectiveness of the proposed method. This method is used to perform denoising pre-processing on the dynamic monitoring data of Sutong Bridge, and the denoised results of different methods are compared and analyzed. Then, the best denoising data are selected as the input data of Hilbert spectrum analysis to identify the structural first natural frequency of the bridge. The results indicate that the wavelet-empirical mode decomposition (EMD) method effectively reduces the interference of random noise and eliminates useless intrinsic modal function (IMF) components, and the excellent properties of the signal evaluation index after denoising make the method suitable for processing non-stationary signals with noise. When Hilbert spectrum analysis is applied to the denoised data, the first natural frequency of the bridge structure can be identified clearly and is consistent with the theoretical calculation. The proposed method can effectively determine the natural vibration characteristics of the bridge structure.


Introduction
Affected by ambient excitation and traffic loads, the complex structure of long-span bridges may have internal and non-linear responses that influence the health of bridge structures. The analysis of the dynamic characteristics of bridge structures under various conditions through the health monitoring of large bridges is essential to the construction and safe operation of bridges. As an important parameter reflecting the dynamic characteristics of bridge structures, the main vibration frequency has become the focus of attention in bridge monitoring. Many studies have showed that the modal frequencies of different structures and their changes reflect the health state of a bridge structure [1][2][3][4][5][6][7][8]. The accelerometers, optical fibers, piezoelectric sensors, and intelligent materials attached to important components of a bridge structure can monitor their local natural frequency characteristics [1][2][3][4][5][6][7][8][9]. However, the displacement deformation characteristics are difficult to determine [9][10][11]. Although total station (even robotic total station), 3D laser scanners, close-range photogrammetry equipment, and the Global Navigation Satellite System (GNSS) have been adopted to monitor the deformation information of bridges, their respective shortcomings limit their wider application [10][11][12][13]. Nakamura [14] used GNSS technology in 1998 to monitor the dynamic deformation of a suspension bridge with a main span of 720 m and verified that the vibration displacement and main frequency of the main beam under wind load were consistent with the results of the wind tunnel experiments and finite element calculations. Given that GNSS monitoring can provide long-term, continuous, high-frequency, and dynamic displacement information and extract dynamic frequencies from it, GNSS technology has become an important means of bridge health monitoring and is widely used in the safety monitoring of bridge structures [15][16][17][18][19].
GNSS dynamic observation data often contain considerable noise due to the interference of external environmental factors, such as multipath effects, and noise-containing GNSS signals usually exhibit nonlinearity and non-stationarity. Thus, exploring the abundant structural health information hidden in GNSS monitoring signals is crucial. Domestic and foreign scholars have conducted studies on non-stationary signal processing methods, such as short-time Fourier transform, wavelet transform, and the Wigner-Ville distribution [20][21][22][23][24]. Although short-time Fourier and wavelet transform can analyze the time-frequency characteristics of non-stationary noise-containing signals, these methods on non-stationary signals may lead to the phenomenon of false components due to the limitation imposed by the uncertainty principle. In addition, the selection of the Fourier window and wavelet basis functions entails high subjectivity, which limits the utilization scope of the two methods to some extent [25]. Hilbert-Huang transform (HHT) is an adaptive local time-frequency analysis method proposed by Huang et al. in 1998; HHT is composed of empirical mode decomposition (EMD) and the Hilbert transform (HT) [26]. The method carries out blind adaptive decomposition in accordance with the signal itself and separates the non-stationary signal into several intrinsic mode functions (IMF) according to the frequency content (from high to low frequency). On this basis, HT is implemented on each IMF component to explore its spectral characteristics [26]. However, due to the influence of noise and signal discontinuity, the mode-mixing problem often occurs when EMD is applied to the decomposition of non-stationary noise-containing signals. This may result in interference on the EMD decomposition results and the subsequent results of Hilbert spectrum analysis; consequently, effective physical information may not be accurately determined [27,28]. Xu et al. [29] improved the HHT method, applied it to the processing of the GNSS monitoring data of Baishazhou Yangtze River Bridge in Wuhan, and verified the effectiveness and rationality of this method when used for modal parameters' identification from non-stationary vibration signals.
With regard to the development of EMD-related algorithms, Wu et al. [30] proposed a noiseassisted data analysis method called ensemble empirical mode decomposition (EEMD), which suppresses the mode mixing problem of EMD by adding white noise to the initial data many times. However, the performance of the EEMD algorithm depends largely on the noise amplitude and the number of trials. The IMFs generated by EEMD become highly polluted and even yield pseudocomponents when noises with inappropriate amplitudes are added and the number of trials is changed inappropriately [31]. Furthermore, the time required by EEMD-related algorithms may increase when the number of trails is increased [32].
On the basis of the remarkable effect of noise when the HHT method is applied to non-stationary data and in consideration of the excellent time-frequency localization properties in the field of filtering and denoising, a time-frequency analysis method based on the combination of wavelet threshold denoising and HHT is proposed in this study. The GNSS dynamic monitoring data of Sutong Bridge are adopted as the study object, to identify the first natural frequency of the bridge structure. The implementation process is presented.
The method implements noise reduction preprocessing on the dynamic monitoring data to obtain reduced decomposition layers in the EMD decomposition process and decrease the marginal effects on the quality of useful signal decomposition. The vibration component is reconstructed according to the correlation of each component and subsequently serves as input data for Hilbert spectral analysis to identify the first natural frequency of the bridge structure. The denoising effect of the experimental analysis on simulated noisy signals proves the effectiveness of the proposed method. The application analysis results of real bridge monitoring data indicate that the proposed method can highlight the dynamic characteristics of bridge structure vibration signals in noisy environments and effectively identify the first natural frequency, which is conducive to the scientific evaluation of the safety status of the bridge structure.

Basic Principle of HHT
Hilbert spectrum analysis is performed using HHT from the non-stationary signal. The HHTbased time-frequency analysis method consists of two main steps: (1) multi-scale decomposition, which is implemented on non-stationary signals using the EMD method and whose result is a series of IMFs, and (2) HT, which is performed on each IMF component represented in the joint timefrequency domain. Then, the Hilbert spectrum of the signal is acquired. The steps of EMD method [26][27][28] are as follows: 1) The extreme points of non-stationary signal x(t) are fitted using the cubic spline function mentioned in [26] to obtain the upper and lower envelopes of the signal. The mean is calculated as: 2) The remainder is calculated as: Afterwards, whether 1 ( ) h t satisfies the definition of IMF is determined [26]. If the definition is satisfied, then 1 ( ) h t is the first IMF component and denoted as: where 1 ( ) k h t is the value of 1 ( ) h t after the k th iteration. Otherwise, iterative calculation continues until the condition mentioned in [26] is met.
3) 1 ( ) c t is separated from   x t . Then, a differential signal with the high-frequency components removed is obtained as follows: 4) The preceding steps are repeated many times by taking 1 ( ) x t (the latter counterpart is recorded as ( ) n x t ) as a new signal [26], and the individual IMF components are decomposed sequentially. The decomposition process terminates when residual signal x is a monotonic function. Then,   x t is decomposed into n IMF components and one remainder through the above- The decomposition process shown above indicates that the EMD model separates the nonstationary signal into several intrinsic mode functions, which are all stationary or stationary-trended signals (i.e., it is an adaptive signal decomposition method). The obtained ( ) The analytic signal is constructed as: The amplitude function is obtained as: Furthermore, the instantaneous frequency is described as: Therefore, where Re is the real part of the original signal, and residual function ( ) n r t can be ignored.
Equation (10) is the Hilbert spectrum, which is recorded as: The Hilbert marginal spectrum is defined as: where T represents the total length of the signal,

Time-Frequency Analysis Method Based on Wavelet Threshold Denoising and HHT
In real engineering applications, the dynamic observation data collected from the GNSS receiver may contain many components aside from the vibration signals from the bridge structure itself, and the data may also contain interference from environmental noise. Therefore, the signal-to-noise ratio (SNR) is relatively low, which greatly affects the recognition accuracy of the main vibration frequency. Only by filtering out noise effectively can useful information and reasonable and reliable conclusions be obtained. For the EMD method, the modal aliasing phenomenon may occur due to intermittent events and the interaction between signals, resulting in incomplete decomposition. The number of decomposition layers from the EMD decomposition process and the margin effect on signal analysis must be reduced via denoising pre-processing to improve the accuracy and timeliness of signal feature extraction.

Wavelet Threshold Denoising
The wavelet threshold method is popular in signal denoising, whose basic idea is to denoise the non-stationary signal as follows: Firstly, multi-scale decomposition is performed on the noisecontaining signal. Secondly, an appropriate threshold function is selected to threshold the wavelet coefficients. Lastly, the inverse wavelet transform is used to reconstruct each signal and achieve denoising. If signal f(t) is a square integrable function, then the wavelet transform of f(t) is the inner product of the signal and wavelet function where a represents the scale factor, b represents the translation or displacement factor, , a b R  , and the conjugate function of  is denoted as  . Wavelet transform analyzes signal by scaling and shifting the position of ( ) t  . Thus, the wavelet transform results possess good time-and frequency-domain localities when the wavelet function is selected properly. Selecting the appropriate threshold function in the practical application of wavelet threshold denoising is extremely important because it minimizes the noise in noisy data and preserves the local characteristics of the effective signal. Although hard and soft threshold denoising methods are widely used in data filtering, they still have shortcomings. The works in [25][26][27][28][29] proved that the semithreshold function shown in Equation (14) may obtain good denoising results.
where B and C is a constant that usually has a value greater than one. The threshold function has continuity and high-order conductivity. It overcomes the discontinuity and oscillation problems of the signal reconstructed from the hard threshold function and addresses the deficiency that although the soft threshold function is continuous, the denoising results need to be subtracted from the threshold; systematic deviation exists when the coefficient exceeds the threshold.

Time-Frequency Analysis Method Combining Wavelet Threshold Denoising and HHT
Although the wavelet threshold method can eliminate white noise in the original signal, the noise reduction preprocessing may not sufficiently eliminate the interference of noise in practical applications due to the noise interference of many different properties. Therefore, the IMF correlation must be adopted for appropriate post-processing, and it may improve the recognition accuracy of the main vibration frequency. The denoising method based on wavelet threshold denoising and EMD is called the wavelet-EMD method, and its specific steps are as follows: 1) Wavelet threshold denoising: To eliminate the influence of noise, the original noise-containing non-stationary signal ( ) x t is decomposed by the wavelet threshold, and the signal is reconstructed according to the energy of the decomposed frequency band. The denoising signal that eliminates interference noise is then obtained. This denoising process via the wavelet threshold can decrease the influence of high-frequency noise and the decomposition layers of EMD. Hence, it may provide a relatively "clean" input signal for EMD.
2) Multi-scale EMD is carried out on the denoised signal ( ) x t  processed by the wavelet threshold denoising method to obtain the IMF components of the denoised signal, 1 2 , , , n c c c  .
3) The correlation value of each IMF component is calculated to identify the vibration component of the bridge structure.
For the decomposition results of large-scale building dynamic observation data, the highfrequency parts of the IMF component are generally dominated by noise, and the low-frequency part is dominated by the vibration signal. The investigation of the decomposition results of EMD indicates that several components have little energy after decomposition and cannot represent the original signal, which is not conducive to analyzing the spectral characteristics of the signal. To identify the effective vibration component easily, the component correlation i  according to the correlation between each component and the original signal is defined to distinguish the effective IMF component from the unwanted IMF component.
where ( ) where ( ) The quality of denoising performance is the key to evaluating the advantages and disadvantages of a denoising algorithm in practical applications. Therefore, an objective denoising performance evaluation index must be selected. Similar to the evaluation index of image denoising [33], the root mean squared error (RMSE) and SNR are usually chosen for such an evaluation. RMSE and SNR are expressed as follows: where ( ) y i is the denoised signal, ( ) x i is the original signal that serves as the standard signal, and N is the length of the signal. RMSE reflects the proximity of the denoised signal to the original signal.
The smaller the value is, the more obvious the obtained filtering effect is. Meanwhile, SNR reflects the ratio between the denoised and the noisy signal, and its value symbolism is contrary to that of RMSE. With regard to the denoising the bridge GNSS dynamic observation data, the pre-processing of dynamic data must retain the useful signal as much as possible whilst removing the noise to the greatest extent. Thus, the two evaluation indicators must be considered together when evaluating the denoising effect of non-stationary signals. That is, the RSME of the denoised and original signals and the SNR of the filtered signal should meet the requirements.

Comparative Analysis of the Denoising Effect on Simulation Signals
Several commonly used standard analogue signals, such as the blocks, bumps, Doppler, and heavy sine signal, were selected for a simulation analysis, which aimed to verify the performance of the proposed denoising method and examine the effective suppression of the modal aliasing phenomenon in the traditional EMD method. The denoising performance of the proposed method was compared with that of wavelet denoising and EMD methods. The two evaluation indicators were used to evaluate the denoising performance of the four methods quantitatively.
Classic simulation signals were generated by the Wnoise function in MATLAB 2014b. These signals had the same sampling length and different SNR. Their waveforms (SNR were all 4 dB) are shown in Figure 3.       Table 1 indicates that each method performed well to a certain degree in terms of the denoising effect on the simulation signals. Overall, the denoising effect of the wavelet denoising method was not as good as that of EMD and was inferior to that of wavelet-EMD. The wavelet-EMD method performed well not only in the evaluation indicator SNR, but also in RMSE, thereby satisfying the demands mentioned in Section 3.2.
To investigate this issue further, noisy signals with different SNRs (i.e., 6, 2 and −2dB) were analyzed in the same manner as above. The wavelet functions for the comparative analysis used the db8 function and the semi-soft threshold method. The waveforms of the denoising results from the different methods were omitted to save space; only the denoising effect indicators after processing are listed in Table 2.  Table 2 together with Table 1 show that the denoising effect of the different methods was verified again. Several conditions were inapplicable to the blocks signal. The denoising effect varied with the SNR of the noisy signal. The wavelet-EMD method exhibited good adaptability, and its denoising effect had little correlation with the noise intensity of the noisy signal. In particular, even when the SNR of the noisy signal was relatively low, the denoising result from the wavelet-EMD method was smooth, and the method could still denoise effectively. The superiority of the proposed algorithm was proven again.

Application Analysis
Sutong Bridge, which is in the lower reaches of Yangtze River, is a steel box girder cable-stayed bridge with double towers and double cables. Its main span is 1088 m long, and its main tower is 300.4 m high. Sutong Bridge ranked second amongst similar-type bridges in the world at the time of its completion. Evaluating the stage line shape, force, and safety of the bridge structure and determining the safety status of the bridge structure under environmental excitation conditions are vital. An all-weather dynamic geometric monitoring system based on GNSS technology had been established to monitor the geometry and structural state of towers and beams continuously and in real-time, to analyze the construction status, and to avoid unfavorable construction conditions [16]. The GNSS monitoring points were set synchronously according to the needs of construction control and monitoring during the installation process of the steel box girder. The overall layout of the monitoring points is shown in Figure 7. The experimental data were obtained from the monitoring of a point on Sutong Bridge under a normal construction environment on May 1, 2007; the weather condition was good; the altitude cutoff angle of the satellites was set to 15°; and the continuous observation in the dynamic observation mode lasted for approximately 3 h. Given that the natural frequency of long-span bridges usually falls within the range of 0-2 Hz, the monitoring sequence with the frequency range of 0-5 Hz could be identified when the sampling frequency of the GPS receiver was set to 10 Hz according to the Nyquist sampling theorem of signal processing. This identified range satisfied the requirements of first natural frequency identification. The GNSS dynamic observation data included the vibration signal of the bridge structure under environmental excitation and accidental errors, such as the multipath effect, because of the influence of the construction and the natural environment. Therefore, the rich bridge structure health information embedded in GNSS monitoring signals must be mined. The observation sequences of the monitoring point in the x (lateral) and y (longitudinal) directions lasting for 1 h were selected to identify the first natural frequency of the bridge structure using the proposed method and to analyze the dynamic characteristics of the bridge structure. The obtained coordinate sequence was transformed from a geodetic coordinate system into a bridge axis coordinate system by using the same method as that in [11]. The time history curve of the averaged data from the GNSS monitoring point is shown in Figure 8. The accuracy analysis of the observation sequence in [16]    processing on the original observation sequence was necessary to improve the accuracy of identifying the first natural frequency. In the experiment, the wavelet-EMD denoising method based on wavelet threshold denoising and the EMD method were used to denoise the observation sequence in the horizontal and vertical directions in Figure 8. The results were compared with those of hard, soft, and semi-threshold wavelet denoising (the wavelet function selects the Sym wavelet, and the decomposition layer selected four layers; the wavelet-EMD method used the semi-threshold function). The results of the four denoising methods are shown in Figure 9, where H-Wavelet represents hard threshold denoising, S-Wavelet represents soft threshold denoising, I-Wavelet represents half threshold denoising, and WE represents wavelet-EMD denoising. Figure 9a shows the denoising results in the horizontal direction obtained using the four methods, and Figure 9b shows the denoising results in the vertical direction derived with the four methods.  Figure 9 shows that the wavelet-EMD method obtained a large reduction in Gaussian white noise in the horizontal and vertical directions compared with the three other methods. The signal denoised by the wavelet-EMD method also filtered out most of the invalid IMF components, which could highlight the original signal information. The SNR and RMSE of the denoised signal and the linear correlation coefficient between the denoised signal and the original observation sequence were calculated to evaluate the denoising performance of the four methods quantitatively. Table 3 shows a comparison of the denoising performance of the four methods in the horizontal and vertical directions. A comparative analysis of Figure 9 and Table 3 indicates that all four denoising methods played a certain denoising role, but the denoising effects differed slightly in different directions. In the comparison of the H-Wavelet and S-Wavelet methods, the three evaluation indexes indicated that the H-Wavelet method had a better denoising effect in the horizontal direction, whereas the S-Wavelet denoising method had a better denoising effect in the vertical direction. The wavelet-EMD method exerted a stable denoising effect in both directions. The method had a lower RMSE of the denoised signal and a higher SNR than the three methods, thus suppressing the noise effectively. This result proved that wavelet-EMD achieved the best denoising effect amongst the four methods. In addition, the linear correlation coefficient between the denoised and original observed signals in Table 3 indicated that the correlation between the signal denoised by the H-Wavelet method and the original signal was the weakest; the correlation between the denoised signal using the wavelet-EMD method and the original signal was the strongest, which also reflected the good denoising performance of the wavelet-EMD method. The spectral values of the signals in the horizontal and vertical observation sequences denoised by the four methods were calculated to compare the methods' denoising performance further. The results are shown in Figure 10. The comparative analysis of the spectral values of the denoised signals from the four methods in Figure 10 (a) and (b) indicated that the main frequencies of the observed sequences filtered by the four methods were prominent, which meant that these methods played a certain role in noise inhibition. However, the noise spectrum power of the denoised signal via the H-Wavelet and S-Wavelet methods was relatively large, and the frequency of the noise still filled the entire frequency band. The noise spectrum power of the denoised signal via the I-Wavelet method was reduced to a certain extent, and the noise spectrum power was the smallest in the case of the wavelet-EMD method. Noise was suppressed well in the high-frequency band greater than 2 Hz, indicating that the method effectively removed the interference of noise. Thus, the denoised data obtained by the wavelet-EMD method served as the input data of the Hilbert spectrum analysis in this study to extract the first natural frequency of the bridge structure.

First Natural Frequency Identification from the Denoised Results
Given that the HT method was sensitive to errors during spectral analysis, the accuracy of Hilbert spectrum analysis could be improved to some extent if the quality of the original observation data were improved by appropriate data preprocessing methods. Therefore, the "clean" monitoring sequence after denoising that served as the input data was vital for the identification of the main vibration frequency of the bridge structure, which was why noise reduction pre-processing of the dynamic observation data was performed in Section 4.1. In this section, the Hilbert spectrum analysis is performed on the observed data sequence after denoising, and the vibration frequencies of the bridge structure in the horizontal and vertical directions are identified. Meanwhile, the identification value in this study and the empirical calculation value must be compared to determine the health and structural safety status of the bridge, which would verify the reliability and accuracy of the first natural frequency identification result.
The EMD process was implemented on the original noise-containing observation sequence in the horizontal direction of the bridge, then on the denoised signal via the wavelet-EMD method. The purpose was to verify that the time-frequency analysis method based on wavelet threshold denoising and HHT could effectively extract the first natural frequency of the bridge structure under a noisy    A comparative analysis of Figures 11 and 13 indicates that the frequencies of IMF1 and IMF2 components obtained from the EMD method before denoising were distributed in a disorderly manner, and the correlation value of the component was low. IMF1 and IMF2 were high-frequency noise components according to the global and high distribution of the accidental error. Figure 12 shows that the high-frequency noise in the original observation sequence was suppressed effectively via wavelet-EMD denoising, and the number of decomposition layers of EMD was decreased. Given that the natural frequency of a long-span bridge structure is low, which is generally in the range of 0.1-1 Hz [29], the correlation values of the first three IMF components of the denoised signal were high, and those of the last two IMF components were low ( Figure 14). Therefore, invalid trend components IMF4 and IMF5 may be removed from the denoised signal, and effective components IMF1-IMF3 were selected for reconstruction to obtain the vibration component of the bridge structure.
The Hilbert time spectrum and the marginal spectrum of the bridge structure were obtained via Hilbert spectral transformation from the reconstructed vibration components and are shown in Figures 15 and 16, respectively. The results proved that the Hilbert marginal spectrum could accurately reflect the main frequency with the actual physical component. Similarly, multi-scale decomposition was performed on the observation data of the bridge in the vertical direction denoised by the wavelet-EMD method, and the correlation value was calculated. Then, the effective vibration components were selected for reconstruction. Finally, the Hilbert time spectrum and marginal spectrum of the bridge structure in the vertical direction were calculated and shown in Figures 17  and 18     However, the reliability and precision of the first natural frequency identification results based on time-frequency analysis methods are often difficult to ensure in practical applications. The parameter values identified in this study must be compared with those of other different methods. The theoretical calculation value of the bending natural vibration frequency of Sutong Bridge in the horizontal direction was 0.145 Hz according to the theoretical and empirical calculation methods of the dynamic characteristics of long-span bridges [34][35][36], and the value of the vertical direction was 0.133 Hz. A comparison of the identification results with the theoretical calculations showed that the relative errors of the bridge structure in horizontal and vertical directions were 5.52% and 4.67%, respectively, and the statistical results are shown in Table 4.  Table 4 shows that the amplitude of the main frequency of the bridge in the horizontal direction (4.50 mm) was significantly larger than that in the longitudinal direction (1.02 mm), and the main frequency in the vertical direction of the bridge was not as prominent as the horizontal one. This result indicated that under environment excitation conditions, the loads, such as the construction environment and wind, exerted a larger influence on the longitudinal direction of the bridge. Similarly, analysis of the observation data from the monitoring points in other periods indicated that the first natural frequency was the same in horizontal and vertical directions, except for a slight difference in amplitude. This result indicated that the wavelet-EMD method was used to denoise the original observation data firstly. Then, Hilbert spectrum analysis was performed on the denoised data to obtain the accurate first natural frequency of the bridge structure. This scheme had high identification accuracy and could be used for dynamic characteristic extraction of bridge structures.

Conclusions
With the GNSS dynamic monitoring data of Sutong Bridge as the study object and in consideration of the low precision of first natural frequency identification due to noise interference caused by HHT processing of non-stationary signals, this study proposed a time-frequency analysis method based on the combination of wavelet threshold denoising and HHT. This method was applied to the analysis of the dynamic characteristics of a bridge structure. The conclusions were as follows: 1) The proposed wavelet-EMD method suppressed high-frequency noise effectively, decreased the decomposition layers of EMD, and reduced the margin effects on the quality of effective signal decomposition. The SNR and linear correlation coefficient of the denoised signal were the largest, and RMSE was the smallest. The method had a stable denoising effect in the horizontal and vertical directions.
2) The Hilbert spectrum analysis of the denoised data clearly reflected the spectral value of the bridge structure, and the numerical results agreed well with the theoretical calculations. The relative errors of the natural frequency identification in the horizontal and vertical directions were 5.52% and 4.67%, respectively, which meant that the natural vibration characteristics of the bridge structure were identified effectively.

HHT
Hilbert-Huang transform EMD Empirical mode decomposition wavelet-EMD Wavelet-empirical mode decomposition IMF Intrinsic modal component GNSS Global Navigation Satellite System HT Hilbert transform Author Contributions: X.W. and S.H. participated in the design of this study, and they both performed the statistical analysis. X.W. carried out the study and collected important background information. All authors read and approved the final manuscript. X.W. and S.H. carried out the concepts, the design of data acquisition, the definition of the new algorithm, the literature search, the data analysis, and the manuscript preparation. X.W., C.K., G.L., and C.L. carried out the data acquisition, the literature search, and the manuscript editing. S.H., X.W., C.K., and C.L. performed the manuscript review. All authors read and approved the content of the manuscript.