Application of the wavelet transform to sediment grain sizes analysis with an impact plate for bedload monitoring in sediment bypass tunnels

An impact plate (IP) is a bedload transport monitoring device developed for a part of the sediment bypass tunnels management. In the measurement, the impact produced by bedload is recorded as the number of impulses (Ip) which is widely used in Japan. Ip, however, has several shortcomings attributed to the insufficient raw signal denoising. In this study, the discrete wavelet transform (DWT), an advanced signal processing technique especially for noisy, non-periodical, and transient signals, was introduced to devise an improved Ip count system solving the problems in the original signal denoising process. The presented results revealed that the DWT is useful for water noise reduction, signal overlap reduction, and mitigating Ip saturation at grain sizes Ds = 50 and 100 mm compared to the original Ip counting system.


Introduction
A crucial issue in the long term use of dams is reservoir sedimentation caused by incoming sediment from the upstream reaches, and the sediment bypass tunnel (SBT) is an effective measure to mitigate the problem.Meanwhile floods occur, the SBTs route sediment coming from the upstream to downstream river hence the sedimentation is mitigated [1].Although SBTs are working effectively in Japan and Switzerland, the long-time use of SBTs has revealed that the SBTs inevitably suffer from severe invert abrasion which requires regular rehabilitation costs.Accordingly, the accurate abrasion estimation is desired for the better management and operation of SBTs.For the use of existing abrasion prediction models, the amount and grain sizes of the impacted sediment onto the invert are essential [2].
An impact plate (IP) is one of the surrogate bedload monitoring systems developed for sediment monitoring in SBTs.In general, the surrogate bedload monitoring systems observe the oscillation or acoustic energy caused by sediment impact onto the device placed in the river bed to estimate bedload movement with this.The IP was developed based on widely used two systems, namely Swiss plate geophone and Japanese pipe microphone (JPM) [3] with aiming high sensitivity to detect gravels with grain sizes Ds > 2 mm and adequate robustness against the large gravels impact with high flow velocities in SBTs V > 10 m/s [4].In the bedload measurement with JPM in Japan, sediment impact is recorded as a summary value called the number of impulses (Ip), and the IP is also following the recording system [5].Ip was introduced for estimating the sediment amount and grain size information with a small volume of data.Ip, however, still has several shortcomings as listed below, and presumably, these are associated with insufficient signal denoise in the original Ip count:  Ip produced by water noise hinders that by small particle impact  Underestimated record of Ip due to signal overlapping (long reverberation)  Limited detectable upper Ds due to signal saturation Recent technology has allowed us applying more advanced signal processing techniques such as the wavelet transform (WT).WT was developed for overcoming disadvantages of the short-time Fourier transform (STFT).The use of WT gives profits to analyze signals with non-periodical, noisy, and transient features [6].Given the IP signals containing the aforementioned features, WT is worth applying to improve the original Ip count system.
This paper attempts to apply WT for denoising sediment impact signals recorded by IP.Then, we evaluate the availability in reducing aforementioned shortcomings for the improved Ip and raw signals computed with WT.

Wavelet transform
In this section the principle of WT is briefly introduced.Firstly, the continuous wavelet transform (CWT) is presented with contrasting STFT, then the discrete wavelet transform (DWT), used in the present study, is explained.
STFT is a technique to overcome a disadvantage of the Fourier transform (FT) that the time information is lost in the frequency domain (frequency spectrum).STFT of function x(t) at x = x0 is defined as : where ω is a frequency and w(x-x0) is a window function.Eq. ( 1) means that the signal x(t) is split by the window function defined by t0 and the splitted signal is Fourier transformed.Therefore, STFT can give the frequency spectrum of x(t) around at t = t0.
Although STFT makes it possible to analyze time variation of frequency spectrum, it does not allow to obtain both good time and frequency resolutions.This causes, for example, setting the window short for observing time variation in detail inevitably invites the reduction of the frequency resolution.This trade-off relationship between time and frequency resolutions is known as Heisenberg uncertainty principle.
WT is a technique developed based on the idea to expand and contract the window function for avoiding the problem in STFT.The CWT is defined as: where continuous parameters a and b are dilation and translation factors, respectively, and ψa,b defined as Eq. ( 3) is a function called the mother wavelet.Comparing with Eq. ( 1), STFT can be interpreted as a method to express a signal with the superposition of sine and cosine curves, whereas CWT expresses that with the superposition of ψa,b.The mother wavelet, dilated by a and shifted by b, is corresponding to e -iωt w(t-t0) in STFT.The dilation and shift of ψa,b make it possible to apply long time resolution to low frequency components of the signal and vice versa, hence CWT realizes both good time and frequency resolutions.
For the preparation of ψa,b, many functions have been suggested and all of them have different pros and cons [7].In this study, Daubechies wavelet, one of the most widely used mother wavelets, was selected to analyze IP signals.
In DWT, the dilation and translation factors (a and b) have discrete values.The discrete wavelet is defined by discretizing a = 2 j and b = 2 j k in Eq. (3) as : /2 , ( ) 2 (2 ) here, j is a parameter called level which scales the wavelet by every 2 times and k plays a role of shifting.A remarkable benefit of DWT is the availability of multiresolution analysis when ψj,k is an orthogonal wavelet [6].In this instance, a signal can be expanded in terms of discrete wavelets.This expansion is equivalent to the decomposition of a signal as : here, dj is called detail at level j and aj is approximation at level j.Eq. ( 5) illustrates a signal is decomposed into d1+a1, then the a1 is decomposed into d2+a2, and this process is repeated until j = J.Literally, aj provides an approximative component of x(t) and the approximation becomes coarser as j increases.Therefore, respective dj and aj are corresponding to high-pass and low-pass filters (see Section 3.4).The detailed information of the multiresolution analysis is introduced in e.g., [6] and [8].

The Impact plate
The impact plate (IP, manufactured by Hydrotech Co., Ltd., Japan, Fig. 1), consists of a microphone and an accelerometer mounted underneath a steel plate (49.2 cm × 35.8 cm × 1.5 cm), records impact caused by sediment hitting on the plate.The plate is embedded in an SBT invert and isolated with a rubber material in order not to detect impacts transmitted from the plate vicinity.Particle impact on the plate is registered with a sampling rate of 50 kHz.Figure 2a exemplifies an acoustic signal produced by several gravels of Ds = 10 mm.

Flume experiment setup
A flume experiment was done to determine the behavior of the IP for various water velocities V and grain size fractions Ds at the Laboratory of Hydraulics, Hydrology and Glaciology (VAW) of ETH Zürich, Switzerland [4].For that, we used a flat horizontal flume with 7.56 m length, 0.5 m width, and the flow depth was fixed to 0.08 m.The IP was mounted at the flume end and the test sediment was fed from the upstream reach.As in Table 1, uniform grain size sediment with five Ds were chosen (mean Ds = 2, 5, 10, 50 and 100 mm), with two flow discharges Q = 0.10 and 0.18 m 3 /s resulting in respective V = 2.5 and 4.5 m/s.Each case was repeated around 50 times, hence approximately 500 runs were conducted in total.The IP started recording concurrently with the feeding of sediment into the flume and recorded sediment impact until all the sediment completely flowed out of the flume.
All gravels used in the experiment were sampled from a natural river, and, thus, had various shapes, including round, irregular, and angular shapes.Test sediment in cases with Ds = 2 and 5 mm weighted 50 g of sediment, and cases with Ds = 10, 50, and 100 mm numbered 20 particles in a run.

The original impulse count
Since the recording raw signals requires a huge volume of storage and electricity, a data compression method called the number of impulses (Ip) has been widely used for sediment monitoring with JPM [5].Ip is counted simultaneously with the raw signal recording through analog signal processing in a data logger (Fig. 2).Firstly, the raw signal (Fig. 2a) is sent through a band-path filter to extract a frequency of around 4.6 kHz, which was previously computed as the most effective frequency by the manufacturer (Fig. 2b).Subsequently, the filtered raw signal is transformed into absolute values hence all the negative values are changed to positive ones (Fig. 2c).Then, an envelope curve is generated and the enveloped signal is exported to 6 different channels in which the wave is amplified by factors of 2, 4, 16, 64, 256, and 1024.Finally, the Ip value for each amplification factor (Amp) is counted for each envelope, exceeding a pre-selected threshold of 2 V (Fig. 2d).In this original Ip count process, the band-pass filter plays a role to denoise signals.

The improved impulse count with DWT
For finding well denoised raw signals, each raw signal was once decomposed through DWT with Daubechies 6 wavelet at level 5.The original signal was decomposed into five levels of details and one approximation.Then, each detail coefficient was shrinked (so called softthresholding [9]).The threshold was the maximum value of detail coefficients at each level, that appeared in a signal produced by water flow of V = 4.5 m/s without any sediment.This procedure meant to reduce noise caused by water flow.Finally, each coefficient was reconstructed to the time domain signal.2cd).The threshold was defined as a value being double of the maximum value appeared in each generated signal in Case 5 (Ds = 100 mm) of the first run, since the value might be enough large for all the experiments.12 levels of amplification factors (Amp = 2 0 ,…,2 11 ) were employed.

Result and Discussion
In this section, Ip computed from denoised signals with DWT is compared with the original Ip.In addition, the denoised signals are also shown to visually confirm the improvement in reducing noise.Hereafter, let Ipi be Ip computed with the amplification factor of j.

The original number of Impulse
Relationship between the original Ip and Amp for all Ds are shown in Figure 4.Each line provides the mean Ip with an error bar (σ = 1) versus six base-2 logarithm of Amp (= 1, 2, 4, 6, 8, and 10).The detailed interpretation of the original Ip result is described in [4].
These results indicate several shortcomings on estimating the volume and grain sizes of the transported sediment.Firstly, Ip both of Ds = 50 and 100 mm counted almost the same Ip with Amp = 1, additionally, both of V = 4.5 m/s provided almost the same shape.This result explains the signals by particles of Ds = 50 and 100 mm saturated even after the denoising.Therefore, the upper limit of the detectable grain size might be limited to Ds = 50 mm.Secondary, this procedure counted Ip in water noises under the condition of V = 4.5 m/s.According to the measurement of water flow without sediment revealed that water noise generates Ip1024 = 1.8 and Ip256 = 0.4 per a second in average.Since single experimental run spent 10 to 20 seconds of recording, Ip1024 and Ip256 in Figure 4b are not valid data.Furthermore, it has been often reported that high flux sediment flow affects the linear relationship between Ip and sediment volume because the condition gives frequent signal overlapping resulting in the underestimation of Ip.This also implicates the insufficient denoising in getting rid of unnecessary reverberant of signals.

The number of Impulse with DWT
The three results illustrating relationship between Ip with DWT denoising and Amp are in Figure 5. Ip was computed from signals consist of d1+d5+a5 (Fig. 5a), a5 (Fig. 5b), and d2+d3+d4+d5+a5 (Fig. 5c), respectively.Overall, Ip increases as Amp similar to Figure 4.The variance of Ip seems equal or less than the original Ip.
In Figure 5a, a clear monotonic increase can be seen except for Ds = 2 mm.Particularly in Amp = 1 and 2, particles of Ds = 100 mm produced Ip but that of Ds = 50 didn't, as the saturation of signals were alleviated.A deficit in the result is difficulty in detecting Ds = 2 mm, thus we attempted to use different reconstructions for different Ds.
Figure 5b seems to have an advantage for analyzing large grain sizes Ds = 50 and 100 mm.Clear monotonic increase was observed with a small variance, moreover the Ip of small Ds particles were cancelled by the denoising.According the fact that the result was computed from a5 only, signals produced by small Ds particles contained few relatively low frequency components.This tendency has been seen in the measurement with other surrogate monitoring systems and also be in line with Hertz's contact theory [10].This finding might be useful in the point of monitoring a large Ds portion of sediment in SBTs because it is known that the saltating coarse gravels are the major cause of SBT invert abrasion [2].
For relatively small Ds (2, 5, and 10 mm), Figure 5c provides clear count of Ip.As for Ds = 50 and 100 mm, Ip once decreased seemingly because of the signal overlapping.Although the variance of them are slightly larger than other results, small Ds particles generated relatively higher number of Ip to that of Ds = 50 and 100 mm.Therefore, the incorporative use of a5, and d2+d3+d4+d5+a5 possibly increases the accuracy for grain size analysis.Finally, the benefit of DWT denoising is checked for raw signals (Fig. 6). Figure 6a shows one of the d2+d3+d4+d5+a5 signal where Ds = 2 mm, V = 4.5 m/s (a base of Fig. 5c, positive portion) with the original raw signal and the 4.6 kHz band-passed signal.Apparently, DWT effectively denoised the signal and accentuated impact caused by the particles.It is also favorably seen that the threshold corresponding to the highest Amp was enough larger than the water noise, hence few water flow-derived Ip were counted.
While, Figure 6b shows d1+d5+a5 signal where Ds = 100 mm, V = 2.5 m/s (a base of Fig. 5a).This indicates the DWT well weaken the reverberant time of signals, hence the problem of signal overlapping would be mitigated.

Conclusions
We applied DWT for denoising sediment impact signals recorded by IP to develop improved impulse count system.Results revealed that the DWT denoising was useful for water noise reduction, signal overlap reduction, and mitigating Ip saturation at Ds = 50 and 100 mm compared to the original Ip count system.We also suggested the combined use of the different DWT denoising for different grain size fractions to improve the accuracy.
The next challenge is the estimation of grain size distribution from mixed grain size test results.Ideally, the grain size fraction (solution space) can be estimated from the observed Ip (observation space) by solving a discrete linear model with using the pre obtained Ip-Amp relation [4].Beside the improvement of Ip count discussed in this study, still some mathematical conditions are required for successfully solving the model.For selecting the best DWT with considering the point, we need to solve the model with various DWT based Ip-Amp relations and compare their accuracies.

Fig. 1 .
Fig. 1.(a) IP installed in an SBT invert (Koshibu dam, Japan), (b) the back side of the IP with a microphone and an accelerometer.

Fig. 2 .
Fig. 2. Process to count the number of impulses Ip, (a) a recorded raw signal, (b) the band-pass filtered signal, (c) transformation into absolute value and enveloping, (d) amplification and Ip count.

Figure 3
exemplifies the DWT decomposing procedure.The original signal (s) is decomposed into five levels of high-passed signals dj after soft thresholding, here dj illustrates level j reconstructed signals (j = 1,...,5).a5 signal is a residual of the decomposition and corresponding to the most approximative (low-passed) composition of the original signal.In contrast, d1 shows the most detailed (high frequency) component of the original signal.From one original signal, totally 64 patterns of signals were generated as a combination of six signals (d1-d5 and a5).Ip was counted based on the generated signals in the same manner as the original Ip count (Fig.

Fig. 3 .
Fig. 3.A schematic chart illustrating the signal decomposition process, the original signal (s) is decomposed with DWT into high-pass filtered signals (d1-d5) and the residual (a5).

Fig. 4 .
Fig. 4. Relation between the original number of impulses Ip and base-2 logarithm of amplification factors Amp for different grain size fractions Ds, (a) flow velocity V = 2.5 m/s and (b) V = 4.5 m/s.

Fig. 5
Fig. 5 Relation between the number of impulses Ip and base-2 logarithm of amplification factors Amp for different grain size fractions Ds, each row presents flow velocity V = 2.5 m/s (left) and V = 4.5 m/s (right).

Table 1 .
Experimental conditions