A wavelet thresholding method for vibration signals denoising of high-piled wharf structure based on a modified artificial bee colony algorithm

2.1. Stationary wavelet transform denoising

Denoising by discrete wavelet transform (DWT) thresholding is widely utilized in the noise reduction [14-16]. However, DWT is a non-redundant decomposition analysis, which leads to the translation non-invariance. The stationary wavelet transform (SWT) [17] was put forward to make the wavelet decomposition invariance as it can remove the down samplers and up samplers in the WT. This property ensures the same size of the subbands at different wavelet decomposition levels, which makes the study of optimal subband thresholds easier. Besides, the redundant decomposition of SWT can also provide a more accurate estimation of the variances at each wavelet decomposition level for recognizing noise [18]. Given the above advantages of SWT over DWT, the stationary transform is utilized in this paper to eliminate noise. As the construction of the wavelet function in the wavelet domain is difficult, Mallat [19] introduced a computationally efficient algorithm known as multi-resolution analysis. On the basis of this method, the signal denoising using the SWT can be divided into the following three steps.

1) With an appropriate wavelet basis function, decompose the noisy signal through SWT, and extract the corresponding approximation coefficients and detail coefficients. Suppose $x\left(k\right)$ is the noisy signal, the following is the stationary wavelet decomposition equation [18]:

where $a_{j,k}$ and $d_{j,k}$ are the wavelet approximation coefficients and detail coefficients, respectively, $h$ and $g$ are impulse responses of low-pass and high-pass filters, respectively, $j$ is the wavelet decomposition scale, $k=\text{0, 1,}\text{…}\text{,}n-1$, and $n$ is the scattered sampling point.

Fig. 1 displays an example of the SWT decomposition process of a signal $f$.

Fig. 1An example of SWT decomposition of signal f

2) The obtained wavelet coefficients are shrunk by a thresholding function, which has an effect on the quality of the denoised signal. The commonly used thresholding functions are hard and soft thresholding functions. As the reconstructed signal is smoother using the soft thresholding function, it is adopted in this paper to modify the wavelet decomposition coefficients. The following is the expression of the soft thresholding function:

where ${W(d}_{j,k})$ denotes the wavelet detail coefficients after thresholding, $\mathrm{s}\mathrm{g}\mathrm{n}(\bullet )$ is the Signum function, $\lambda$ is the selected threshold value.

As the pre-estimation of the noisy signals is always required in the existing rigrsure, sqtwolog, heursure and minimaxi threshold selection strategies, the GCV [20] criterion-based threshold selection method is employed in this paper. This method depends on the input and output data only [21] and is proven to be asymptotically optimal in the Euclidean norm. Therefore, it does not require any prior knowledge of the noisy signals. The corresponding GCV function can be defined by [20]:

where $\mathbf{w}$ is the wavelet coefficient of the real measured noisy signal, ${\mathbf{w}}_{\lambda }$ is the modified wavelet coefficient after using the threshold $\lambda$, $M$ is the total number of wavelet coefficients and $M_0$ is the number of these wavelet coefficients that were set to zero. Under the spline smoothing conditions [22], this threshold choice is asymptotically optimal estimation for the minimizer of the $GCV\left(\lambda \right)$.

3) Wavelet reconstruction of the signal using approximation and modified detail coefficients by means of inverse SWT (ISWT) to estimate the original signals. The wavelet reconstruction equation can be presented as:

where $Wd$ are the modified wavelet coefficients by the thresholding function. From Eq. (4), we can obtain $Wx\left(n\right)=a_{0,n}\text{,}$ which is the estimated signal after wavelet denoising. The reconstruction process of signal $f$ can be seen in Fig. 2.

Fig. 2Reconstruction process of signal f

2.2. Wavelet-based optimization technique thresholding

According to the above denoising steps, we can obtain the goal of SWT denoising is to select the optimal thresholds. With the optimal thresholds, we can shrink the wavelet detail coefficients through thresholding function. Finally, the optimal estimation of the original signal can be reached employing the approximation and thresholded wavelet coefficients. Hence, the signal denoising can be presented as an optimization problem. Different optimization techniques can be employed to solve this problem. If we initialize the thresholds as the possible solutions of the optimization techniques, the optimal thresholds can be obtained by minimizing the objective function within iterative optimization procedure. Fig. 3 shows the basic theory of wavelet-based optimization technique thresholding [14-16].

Fig. 3Wavelet-based optimization technique thresholding [14-16]

3.1. Basic ABC algorithm

In basic ABC (BABC) algorithm, the swarm of artificial bees is composed of three groups of bees: employed bees, onlookers and scouts [23, 24]. The employed bees are occupied to exploit a specific food source and share information about this particular source with the onlookers through waggle dancing. The shared information includes the position and profitability of the food source. An onlooker chooses a food source with a probability related to its nectar information that taken from the employed bees. The chosen probability of a food source increases with the increasing nectar amount. When a food source is abandoned by the employed bees, the corresponding employed bee becomes a scout, and a new food source is randomly generated to replace the abandoned one. It is important to remark that the number of the employed bees or onlookers is equal to the number of the food sources, and both of them account for half of the whole bee colony. The position of a food source corresponds to a possible solution in the search space and the nectar amount represents the quality of the related solution. The implementation steps of the BABC algorithm are presented as follows.

1) Initialize population.

Similar to other evolutionary algorithms, the BABC algorithm searches for the optimal fitness through a randomly initialized bee colony. Suppose $N$ is the number of food sources, where the food source ${\mathbf{X}}_i=({\mathbf{X}}_{i,1},{\mathbf{X}}_{i,2},\mathrm{}\dots ,\mathrm{}{\mathbf{X}}_{i,D})$ represents a candidate solution, $D$ is the dimension of the associated problem space, $i=\text{1, 2,}\text{…}\text{,}\mathrm{}N$. The initialization of the food source ${\mathbf{X}}_i$ can be expressed as:

where ${\mathbf{X}}_{i,d}\left(0\right)$ is the $d$th dimensional component of the $i$th food source, $R$ is a random constant in the range of [0, 1], ${\mathbf{X}}_{\mathrm{m}\mathrm{i}\mathrm{n},d}$ is the minimum value of the $d$th dimensional component and ${\mathbf{X}}_{\mathrm{m}\mathrm{a}\mathrm{x},d}$ is the maximum value of the $d$th dimensional component.

2) Employed bees.

At this stage, each employed bee searches the neighbor of its corresponding source to generate a new food source ${\mathbf{V}}_i=({\mathbf{V}}_{i,1},{\mathbf{V}}_{i,2},\dots ,{\mathbf{V}}_{i,D})$ according to the search rule. If the nectar of the new food source is better than the old one, then retain the new food and abandon the old one. The new food source determination equation can be presented as follows:

where $\varphi$ is a random number in the range of [–1, 1], $j=\text{1, 2,}\text{…}\text{,}N$ is a randomly selected food source.

3) Onlooker bees.

After the searching of all employed bees, the onlookers choose a food source to exploit further according to the probability associated to the nectar amount, which is shared by the employed bees on the dance area. Generally, this choice process is conducted based on the fitness of the possible solution. The corresponding fitness function is computed as:

where $f_i\left(t\right)$ is the objective function of the $i$th food source, $F_i\left(t\right)$ is the corresponding fitness value.

4) Scout bees.

If the fitness of a food source is not improved for a predefined number of trials called “limit”, then that food source is assumed to be exhausted and should be abandoned. The corresponding employed bee becomes a scout, and a new food source is randomly determined by Eq. (5).

3.2. Proposed MABC algorithm

3.2.1. Kent chaotic sequence

The characteristics of chaos technique, including the ergodicity, randomicity and regularity, can be utilized to improve the performance of the BABC algorithm. Attempts have been made in literature [25, 26] to modify the dynamic property of the BABC algorithm by introducing the Logistic chaos, and the so called Logistic chaos-based ABC (LABC) algorithm is proposed herein. However, the Logistic map is a nonuniform ergodicity, which may affect the search speed and efficiency of the LABC algorithm. As compared with the Logistic chaos, the Kent chaos mapping has better uniform ergodicity. Besides, the Kent and Logistic chaos have the same characteristic of topological conjugacy. Therefore, the Kent chaos can also be employed to modify the BABC algorithm. Considering the better uniform ergodicity of the Kent chaos, the chaotic sequences in this paper are produced using Kent map, which can be expressed as:

8

$z_{n+1}=\left\{\begin{array}{c}\frac{z_n}{\beta },0<z_n\le \beta ,\\ \frac{(1-z_n)}{(1-\beta )},\beta <z_n\le 1,\end{array}\right.$

where $z_n$ is a random number in the range of (0, 1), $\beta$ is a control parameter in the range of [0, 1]. The Kent chaotic sequence $z$ can be obtained through gradually iteration. When $\beta$ is taken as 0.4, the sequence of Kent map is chaotic, and the corresponding probability density function satisfies uniform distribution in the rage of [0, 1], namely:

9

$\rho \left(z\right)=1.$

Fig. 4 shows the probability distribution histograms of Kent and Logistic chaos (its control parameter is set to 4) for 40, 000 independent iterations. As clearly presented in this figure, the probability of Logistic map is higher in the ranges of [0, 0.1] and [0.9, 1] compared with the other ranges. However, the probability of Kent map is nearly uniformly distributed in the whole range of [0, 1]. This indicates that the ergodicity of Kent map is better than Logistic map. According to the following steps, the Kent chaotic sequences can be generated.

1) Produce a random initial value $z_0$ in the feasible region.

2) Generate a chaotic sequence $z$ through gradually iteration using Eq. (8).

3) Repeat the above iterative process until the maximum number of allowable iterations is achieved.

Fig. 4Probability distribution histograms of chaotic maps

a) Kent map

b) Logistic map

3.3. Signal denoising using the MABC algorithm

The proposed MABC algorithm is designed to denoise the real measured signals. The flow chart of the designed approach is displayed in Fig. 5. As shown in this figure, the specific steps of this method are presented as follows.

Step 1. Initialize related parameters such as swarm size $S$, number of food sources $N$, search dimension $D$ , maximum number of allowable evolutionary iterations $T$, maximum number of allowable chaotic iterations $t_{\mathrm{m}\mathrm{a}\mathrm{x}}$, exploiting trials of the same food $trial\left(i\right)$, maximum number of allowable trails “limit”, and so forth.

Step 2. Produce D-dimensional chaotic sequence $z_{id}\left(t\right)$ according to the steps of Kent map presented in Section 3.2.1, and generate the chaotic individual $\mathbf{K}\mathbf{X}\left(S\right)$ by the following carrier equation:

Step 3. Compute the opposite solution $\mathbf{O}\mathbf{X}\left(S\right)$ corresponds to the chaotic individuals according to the GOBL strategy described in Section 3.2.2.

Step 4. Sort the solutions by the corresponding fitness values, and take the first $N$ solutions associated with better fitness values as the initial food source $\mathbf{X}\left(N\right)$.

Step 5. Each employed bee searches the neighbor food for the new position ${\mathbf{V}}_i$, and compute the corresponding fitness value $F\left({\mathbf{V}}_i\right)$. If $F\left({\mathbf{V}}_i\right)>F\left({\mathbf{X}}_i\right)$, then ${{\mathbf{X}}_i=\mathbf{V}}_i$, $trial\left(i\right)=0$; else ${\mathbf{X}}_i$ does not change, $trial\left(i\right)=trial\left(i\right)+1$.

Step 6. On the basis of tournament selection probability presented in Eq. (14), the onlooker bees choose a food source to conduct neighborhood search.

Step 7. Memorize the optimal food position if all onlookers are distributed onto food source otherwise return to Step 6.

Step 8. If the value of the counter $trial\left(i\right)$ that recorded the exploit time of the same food is larger than the predefined “limit”, then conduct the Kent chaotic search to produce new food source ${\mathbf{V}}_i$; otherwise turn to next step.

Step 9. Repeat Step 5-8 till maximum number of evolutionary iterations is achieved.

Step 10. Output obtained optimal thresholds and reconstruct the signals using ISWT method.

Fig. 5Flow chart of measured signal denoising using the proposed wavelet-based MABC thresholding

4.1. Simulation setting

The compared four benchmark signals are Blocks, Bumps, Heavy sine, and Doppler. Two degrees of noise are added to these signals and the corresponding SNR of the contaminated noisy signals is 5.0000 db and 10.0000 db for each benchmark, respectively. The classical soft and hard thresholding [9] methods and the existing optimization algorithms-based (i.e., PSO [14], BABC [24], and LABC [25]) thresholding approaches are performed to verify the performance of the proposed thresholding denoising technique. For a fair comparison, all of the denoising algorithms are executed on the same computer configuration (ACPI Multiprocessor PC with Intel-Pentium® 2.13 GHz CPU and 1.92 Gbyte RAM) with MATLAB 2009a. The total number of samples is 2048. The Toolbox for wavelets 4.4 is utilized for the library of wavelet filters. The wavelet transform basis function is db4, and the decomposition level is 3. In the classical thresholding methods, the Universal threshold is adopted for the denoising procedure. For the optimization algorithms-based thresholding approaches, the level-dependent threshold values are adopted so that adapt to the signals and noise features. For different decomposition levels, different threshold values are used. Hence, the number of threshold values should be set to 3, and the possible solution of the optimization algorithms-based thresholding methods can be presented as a three-dimensional vector. The parameters of the four algorithms are set as follows. In PSO algorithm, the swarm size is taken as twice the dimensional search space. The maximum number of allowable evolutionary iterations is set to 50, and the other parameter sets recommended in the literature is used. In BABC, LABC and KABC algorithms, the following basic settings of parameters are the same. The swarm size is set to 6, which includes 3 employed bees and 3 onlookers. The maximum number of allowable evolutionary iterations is the same as PSO algorithm. The maximum number of allowable trails “limit” is given by $(S\times D)/2=9$. Particularly, the maximum number of chaotic iteration $t_{\mathrm{m}\mathrm{a}\mathrm{x}}$ of LABC and KABC algorithms is set to 30. All denoising results of the optimization algorithms-based thresholdings are averaged over 20 independent trails to eliminate the random discrepancy.

4.2. Results and analysis

Table 1 presents the SNR results of the classical wavelet thresholding and mean values of the 20 trails using the optimization algorithms-based thresholding methods. The best results of these methods are indicated by boldface. As shown in this table, except for the same denoising results for the noisy Heavy sine signal with SNR is 5 db, all the optimization algorithms-based thresholding techniques obtain better results than the classical soft and hard thresholding methods. This result reveals the advantages of the adopted optimization algorithms-based wavelet thresholding methods. Besides, the MABC algorithm-based wavelet thresholding obtains the best results among the noise reduction methods in Blocks signal and Doppler signal denoising. For Bumps signal, the MABC algorithm-based thresholding obtains the same results as BABC and LABC algorithms-based thresholding methods, which are the best results among the denoising methods. When the SNR of noisy Heavy sine is 10 db, the MABC algorithm gets the same result as LABC algorithms-based thresholding method, which is also the best one. Hence, we can find the proposed wavelet thresholding method can effectively reduce the noise.

5.1. Experimental process

In order to verify the performance of the proposed wavelet thresholding for actual vibration signal denoising, the test of a physical model for high-piled wharf was carried out for this study. This experiment was conducted at the Key Laboratory of Harbor & Marine Structure Safety of Ministry of Transport in Tianjin Research Institute Water Transport Engineering, China [33]. The experimental physical model is structural segment of a high-piled wharf in the Tianjin Harbor, which is designed with a scale of 1:10. As the soil is involved in this physical model, the experiment cannot fully satisfy the similarity criterion. The wharf is constructed out of reinforced concrete, and is established on the soft clay foundation, which is composed of two layers of clay. The width, length, and breadth of the soil bin are 6.5, 14.5, and 3 m, respectively. The width and length of the wharf model are 2 and 6 m, respectively. The framed bents are 0.7 m apart, and the total number of bents is 9. The sectional drawing of this model [33] can be seen in Fig. 6. All of the components and corresponding sizes of the bents are displayed in this figure.

Fig. 6Sectional drawing of the high-piled wharf model [33]

Forced vibration testing was performed on the high-piled wharf. The excitations were applied by use of hammering method. The instrumented accelerometers number of each vertical piles and oblique piles were 4 and 3, respectively, and the corresponding locations are shown in Fig. 7(a). The accelerometers used in the test were INV 9828 (Piezotronics). The measured vibration signals were collected using a high speed data acquisition system model number INV 3020C. The tests of vertical piles and oblique piles are conducted separately, and each set of test includes two directions (i.e., alongshore and vertical alongshore directions). The total sets of the tests are 4, which are carried out under the same humidity and temperature and by the same operator. The data was acquired at the rate of 256 Hz for 10 seconds. The accelerometers, acquisition instrument and physical model connections are given in Fig. 7(b). As the frequency range of interest is less than 20 Hz, the vibration data is low-pass filtered using a Chebyshev filter with passband cut-off frequency 20 Hz. After pre-processing, the measured data was down sampled to 51.2 Hz, and the corresponding number of sampling points is 512.

Fig. 7Sensor locations (alongshore direction) and experimental connections

a) Schematic diagram of accelerometer locations

b) Experimental connections for model tests

5.2. Real data denoising results and analysis

As the original noise-free vibration signals are unknown, the SNR of the denoised signals, which can be used to evaluate the classical thresholding methods, cannot be calculated. Besides, take the advantage of the optimization algorithms over the classical thresholding methods into consideration, we only use the optimization algorithms-based thresholding methods to denoise the actual measured vibration signals. The measured acceleration response with respect to the three vertical piles and one oblique pile in alongshore direction is observed. To simplify this problem, the response corresponds to the first accelerometer of the four piles is considered. For illustrative purposes, the investigated four sets of signals are recorded as $f_{P1}$, $f_{P2}$, $f_{P3}$ and $f_{P4}$, respectively. The above mentioned parameter sets of the four optimization algorithms are used.

All densoing results of the compared optimization algorithms-based thresholding methods are also computed 20 independent trails. The computation objective values of the four different optimization algorithms, including the best, worst, mean and standard deviation (std) values, are listed in Table 2. The best results of these methods are indicated by boldface. As shown in this table, the proposed MABC algorithm-based thresholding can obtain the optimal denoising results except for the best objective values of $f_{P3}$ and $f_{P4}$. This comparison result indicates the superiority of the MABC algorithm in terms of the global optimization.

Fig. 8 presents the box-and-whisker diagrams relevant to the vibration signals denoising results belonging to the four different algorithms for 20 independent trials. Normally, the bottom and top of the boxes in these diagrams represent the lower and upper quartiles; the bands within the boxes give the median; the whiskers extend to the maximum and minimum of the obtained objective values; the symbols ‘+’ stand for the outliers. From the box-and-whisker displayed in Fig. 8, one may easily notice that the MABC algorithm provides denoising results that read very low amplitude of maximum objective values. The MABC method obtains denoising performance characterized by narrower boxes and lower median values. This result indicates the proposed MABC algorithm can obtain the minimum objective value.

Fig. 8Statistics of the denoising results of the four optimization algorithms-based thresholdings for different measured vibration signals

a)$f_{P1}$ signal

b)$f_{P2}$ signal

c)$f_{P3}$ signal

d)$f_{P4}$ signal

In Fig. 9, the mean best objective values of the above four optimization algorithms for 20 independent trials are compared. From Fig. 9, it is clear that the MABC algorithm converges faster and finds better solutions than all of the other methods for $f_{P1}$ and $f_{P2}$ signals. For $f_{P3}$ and $f_{P4}$ signals, the MABC method has similar convergence speed with LABC method. However, the MABC method has shown to have improved the solutions throughout the denoising process. This result indicates the MABC-based wavelet thresholding outperforms all of the other denoising strategies.

In view of the advantages of the proposed MABC approach, we use this algorithm to denoise the measured signals. Fig. 10 presents the measured signals and the corresponding power spectrum density (PSD) at the first pile of the physical model before and after denoising. As illustrated in Fig. 10(a), the measured signals are contaminated by the interference of the environment because of the relative wide bandwidth of the receiver. After filtering by the proposed thresholding, the signals are clearer than the measured signals. Besides, the denoised signals show lower amplitude of maximum values and have a smaller spread of coordinates. This result indicates the proposed thresholding method is effective in signal denoising. The efficiency of the proposed denoising method can also be verified by the PSD of the signals before and after filtering given in Fig. 10(b). The proposed method can clearly reduce the noise spectrum compared with the unprocessed signals.

Fig. 9Evolution process curves of four optimization thresholdings densoing for different measured vibration signals

a)$f_{P1}$ signal

b)$f_{P2}$ signal

c)$f_{P3}$ signal

d)$f_{P4}$ signal

Fig. 10Acceleration response and corresponding PSD before and after filtering

a) Measured and denoised response

b) PSD of the measured and denoised signals

5.3. Discussion

The denoising experiment is conducted on the basis of physical model of a high-piled wharf structure, which is designed with a certain reduced-scale. As the physical model cannot fully satisfy the similarity criterion, the denoising experiment cannot completely reflect the actual engineering. Therefore, the model errors may impact the filtering results that are not considered in the denoising experiments. However, this problem can be solved if we employ the measured signals from the actual engineering. Besides, we only take one wavelet basis function and fixed wavelet decomposition level so that investigate the impact of threshold determination on the thresholding denoising methods. The validity of the proposed wavelet thresholding method could be verified by the experiment results. However, different wavelet basis functions and decomposition levels may also affect the performance of the denoising approaches. We could consider the influences of these factors in future research. What’s more, the control parameters of the filtering algorithm, namely, threshold values in this paper, are initialized as the position of the optimization algorithm. The optimal threshold values are reached by minimizing the associated objective function. For the random discrepancy of the optimization algorithm, each computation may vary the denoising results together with the threshold values. Further, the optimization algorithm also has sensitive control parameters and their variation may seriously change the denoising results, which may lead to different threshold values as well. We also more concerned about the denoising results than the specific threshold values. Hence, this paper makes the parameters of the filtering algorithm controlled by the optimization algorithm and does not present the specific set of the threshold values.