Wavelet Decomposition Layer Selection for the φ-OTDR Signal

Abstract

The choice of wavelet decomposition layer (DL) not only affects the denoising quality of wavelet denoising (WD), but also limits the denoising efficiency, especially when dealing with real phase-sensitive optical time-domain reflectometry ($\phi$-OTDR) signals with complex signal characteristics and different noise distributions. In this paper, a straightforward adaptive DL selection method is introduced, which dose not require known noise and clean signals, but relies on the similarity between the probability density function (PDF) of method noise (MN) and the PDF of Gaussian white noise. Validation is carried out using hypothetical noise signals and measured $\phi$-OTDR vibration signals by comparison with conventional metrics, such as peak signal-to-noise ratio (PSNR) and structural similarity (SSIM). The proposed wavelet DL selection method contributes to the fast processing of distributed fiber optic sensing signals and further improves the system performance.

1. Introduction

Wavelet transform (WT) stands as a widely employed technique for denoising, filtering, and compressing signals in distributed fiber optic sensing systems, contributing to the enhancement of various sensing technologies, such as Brillouin optical time-domain analyzer (BOTDA), Brillouin optical time-domain reflectometry (BOTDR), Raman optical time-domain reflectometry (ROTDR), phase-sensitive optical time-domain reflection ($\phi$-OTDR), and optical frequency-domain reflectometer (OFDR), among others technologies [1,2,3,4,5,6,7,8,9,10,11,12,13,14]. Its synergy with artificial intelligence methods also helps in time-frequency feature analysis and facilitates feature extraction [15,16,17,18,19,20,21,22]. In addition, signal processing methods such as Empirical Mode Decomposition (EMD), Moving Average (MA), Moving Differential (MD), and Frequency Domain Dynamic Averaging (FDDA) can also eliminate the frequency deviation of $\phi$-OTDR signals, improve their signal-to-noise ratio (SNR), and enhance the accuracy of vibration mode recognition [23,24,25]. Although image processing techniques, such as non-local mean (NLM) [26,27,28], Block-Matching and three-dimensional filtering (BM3D) [29,30], and the neural network-based filtering method [31,32], show superior denoising effects compared to two-dimensional WT, two-dimensional WT has higher processing efficiency with good denoising effects [33]. Additionally, the WD method can be combined with other methods to obtain better denoising results, but essentially, the parameters of WD play a decisive role in the denoising results. Therefore, current research focuses on optimizing wavelet denoising (WD) parameters, including wavelet threshold, threshold function, decomposition level (DL), wavelet basis, etc. [34]. Various strategies have been explored, such as employing an adaptive wavelet threshold based on a simulated annealing algorithm (SAA) to enhance SNR in challenging environments [35]. Furthermore, an improved WT called Maximum Overlapping Discrete Wavelet Transform (MODWT) is introduced to process the vibration signals, achieving remarkable sensing lengths, positioning accuracy, and processing times [36]. Notably, recent efforts in our research have enhanced the wavelet threshold and threshold function to improve the SNR of $\phi$-OTDR and electrocardiogram (ECG) signals [33,37,38]. Despite these advancements, limited attention has been given to the optimization of wavelet DL in recent years. In previous work, we have attempted to select the DL of WD using a signal-based method [33], namely first cycle lag of the autocorrelation function (FCL-ACF). However, the optimal DL cannot be determined because FCL-ACF increased monotonically with DL. In addition, it is difficult to determine the optimal DL by conventional methods. The visual method is a non-objective metric and the optimal DL varies greatly depending on each individual’s subjective vision [39]. The co-occurrence matrices for selecting the optimal DL are limited to image texture detection and periodic filtering [40]. The dominant entropy-based DL optimization methods are monotonically varying exponentials, making it difficult to directly and accurately determine the optimal DL [41,42]. The noise-energy-based methods do not accurately estimate the noise level in practical scenarios. Consequently, there remains a lack of an effective method for determining the optimal wavelet DL [43].

This paper aims to address this gap by proposing a simple and effective wavelet DL optimization method based on the similarity between the probability density functions (PDF) of method noise (MN) obtained from different DLs. This method adaptively selects the optimal DL for the WD, thus effectively avoiding over-filtering caused by unnecessary signal decomposition. The principles of WD and MN are introduced, followed by the validation of the MN-based DL selection method against conventional evaluation indexes, such as the peak signal-to-noise ratio (PSNR) and structural similarity (SSIM), using hypothetical and measured $\phi$-OTDR signals, respectively. The study concludes by investigating the fidelity of the filtered $\phi$-OTDR signal, considering aspects of signal distortion, signal fading, and spatial resolution.

2. Principle

2.1. Wavelet Denoising

The fundamental principle of Wavelet Denoising (WD) involves a threshold shrinkage process applied to wavelet coefficients decomposed into various scales. The sequential steps are outlined as follows:

• Decompose the noisy signal using a designated wavelet basis and Decomposition Level (DL).
• Eliminate the noise component from the wavelet coefficients of each DL using a specified threshold and threshold function.
• Reconstruct the wavelet coefficients subjected to threshold shrinkage through wavelet inversion, resulting in the denoised signal.

By comparing different wavelet families (e.g., daubechies and symlets), we can choose wavelet bases that are similar to the signal or its properties, yielding better signal and noise separation as well as sparsity. Thus, dmey and bior2.2 wavelet bases were chosen to process the hypothetical signals and the measured $\phi$-OTDR signals, respectively. The approximation of these two wavelet bases are shown in Figure 1a and Figure 1b, respectively. The threshold of each wavelet decomposition layer is optimized and selected using an autocorrelation function (ACF)-based thresholding method [28]. The FCL-ACF can be expressed as:

$$G\left(n\right)=\sum _x^M\sum _y^{N}f(x,y)\ast f(x,y-n),$$

(1)

where $M\times N$ is the image size, $f(x,y)$ is the image intensity at position $(x,y)$, n represents the lag from the corresponding position ($x,y$) along the direction of trace, and the symbol ‘∗’ represents the convolution. Since the FCL-ACF monotonically increases with SNR, it can be used to assess the quality of the denoised image [38]. The higher the FCL-ACF is, the better the denoising quality is. While the precise wavelet threshold for each DL can be determined using the image periodic ACF algorithm, there lacks an effective method for selecting the optimal DL. Inadequate denoising may occur with too few DLs, failing to capture signal details, while an excess of DLs can alter original signal components and increase computational costs.

2.2. Method Noise

Method Noise (MN) is defined as the disparity between the noise-containing signal and the denoised signal [44]. It is commonly employed to identify and mitigate potential biases introduced by data analysis methods, thus enhancing the accuracy and reliability of results [28,45,46]:

$$n(i,j)=x(i,j)-y(i,j),$$

(2)

where $n(i,j)$, $x(i,j)$, and $y(i,j)$ represent MN, the noisy signal, and the denoised signal, respectively, and $(i,j)$ denotes the pixel coordinates. If the MN solely comprises Gaussian white noise components, its Probability Density Function (PDF) distribution mirrors a Gaussian distribution. Conversely, if it includes additional signal components, the PDF may deviate from the Gaussian norm. The PDF is visualized through histograms, where MN data are segmented into 1000 equally spaced intervals along the x-axis, with each histogram bar’s height reflecting the data amount within each partition [47]. Due to the high frequency of noise, the MN from the first layer of wavelet decomposition predominantly consists of noise with a PDF distribution closest to the Gaussian normal distribution. Consequently, the SSIM between this PDF and the PDF distributions obtained from other wavelet DLs serves as a determinant for the presence of only a noise component in the MN. This aids in identifying whether the loss of signal components results from an excessive number of DLs. The implementation steps are as follows:

• Obtain the MN and its PDF distribution by subtracting the wavelet denoised signal obtained using the first DL from the noisy signal.
• Similarly, obtain the corresponding MN and PDF distributions for various DLs by subtracting denoised signals obtained from different DLs from the noisy signal.
• Calculate the similarity between the PDF distributions of different DLs and the PDF distribution of the first DL.

The algorithm flowchart of optimal wavelet DL selection method is shown in Figure 2.

• Initialize the parameter of WD as “Dmey” wavelet base, soft threshold function, DL = 1. Then, the wavelet coefficient of each DL are obtained using Matlab command “wavedec2”, and the threshold of each DL is optimized by FCL-ACF based thresholding method [33].
• Denoise the noisy image by WD with the initial parameters.
• Calculate the MN₀ of raw image $Im{g}_0$ and denoised image $Im{g}_1$.
• Calculate the PDF of MN₀ that is used as the reference.
• Increase the DL from 1 to 8 with step of 1.
• Calculate the MN of different DLs, and the SSIM of the denoised and clean images for each decomposition layer. Meanwhile, calculate the FCL-ACF and PSNR of the denoised image according to Equations (1) and (3).
• Repeat the above steps to obtain the SSIM of MN of each DL until DL reaches 8.
• Identify the DL that has the largest rate of change in SSIM of the PDF of MN, i.e., the best layer.

Usually, we use two traditional evaluation metrics, i.e., PSNR and SSIM, to evaluate the denoising results. The conventional PSNR calculated with Equation (3):

$$\mathrm{PSNR}=10{log}_{10}\left(\frac{{255}^2}{\mathrm{MSE}}\right),$$

(3)

$$\mathrm{MSE}=\frac{1}{M\times N}\sum _{x=1}^M\sum _{y=1}^N[Im{g}_n(x,y)-Im{g}_0(x,y)],$$

(4)

where MSE is the mean square error between the clean reference image $Im{g}_0(x,y)$ and the denoised image $Im{g}_n(x,y)$. The SSIM between two non-negative images $Im{g}_x$ and $Im{g}_y$ can be expressed as:

$$SSIM\left(Im{g}_x,Im{g}_y\right)=\frac{\left(2{\mu }_x{\mu }_y+C_1\right)\left(2{\sigma }_{xy}+C_2\right)}{\left({\mu }_x^2+{\mu }_y^2+C_1\right)\left({\sigma }_x^2+{\sigma }_y^2+C_2\right)}$$

(5)

where ${\mu }_x$ and ${\mu }_y$ are the mean intensity of images $Im{g}_x$ and $Im{g}_y$, respectively, ${\sigma }_x$ and ${\sigma }_y$ are the standard deviation and of the image $Im{g}_x$ and $Im{g}_y$, respectively, ${\sigma }_{xy}$ is the covariance of image $Im{g}_x$ and $Im{g}_y$, and $C_1$ and $C_2$ are constants included to avoid instability. It is well known that the larger the PSNR and SSIM values of an image are, the more similar the denoised signal is to the clean reference signals and the better the denoising effect is. For the measured vibration signals, it is not possible to accurately estimate the noise level without a clean signal as a reference, so we choose to use the FCL-ACF metrics as its evaluation standard.

2.3. $\phi$-OTDR

A narrow linewidth laser (NLL) pulse of $\phi$-OTDR forms a distributed coherent speckle pattern along with the fiber that changes with the local light phase. Disturbances such as vibration and temperature will affect the refractive index of the sensing fiber, and thus, the intensity of the speckle pattern is changed. As shown in Figure 3, a 3 kHz NLL emits a continuous laser that passes through a 50:50 optical coupler (OC1) 50% of the light intensity goes through a variable optical attenuator (VOA) as a reference light, and the other 50% is modulated by an acoustic-optic modulator (AOM) and becomes pulse light with a repetition frequency of 10 kHz. Then, it is amplified by an erbium-doped optical fiber amplifier (EDFA) and injected into an 8 km sensing fiber through a circulator (port2). The Rayleigh scattering light returns from port2 to port3 of the circulator. A 1550.12 nm light is filtered out by a dense wavelength division multiplexer (DWDM) and interfered with the reference light through a 50:50 OC2. In the end, the interfered light is detected by a avalanche photodiode (APD) and discretized by a data acquisition (DAQ) card with a sample rate of 80 MS/s. A 100 m sensing fiber at the end of 8 km was used to detect the vibration.

3. Decomposition Layer Optimization Based on Method Noise

3.1. Hypothetical Signal

To validate the efficacy of the MN-based DL optimization algorithm, we conduct a comparison with two traditional indicators: Peak Signal-to-Noise Ratio (PSNR) and Structural Similarity Index (SSIM), the latter representing the similarity between denoised and clean images. Because any signal can be decomposed into a linear superposition of a series of sinusoidal (or cosinusoidal) components at different frequencies, we use a sinusoidal clean image $f_0$ to simulate the hypothetical theoretical $\phi$-OTDR signal:

$$f_0=0.3sin\left(\frac{4\pi }{63}x\right)sin\left(\frac{4\pi }{63}y\right),$$

(6)

where $x,y\in [1,126]$ shows the pixel position of the image. The hypothetical clean image $f_0$ with an amplitude of 0.3 is the same as the one in our previous work [33], one dimension of which is distance and the other is trace, as shown in Figure 4a. The pixels in the image with different colors represent different signal intensities. A noisy image $f_x$ is obtained by adding additional noise X to the clean image $f_0$, as illustrated in Figure 4b:

$$f_x=f_0+X.$$

(7)

where noise X with the amplitude of 0.2 is generated by the function of ‘randn’ in the MATLAB software and it follows a Gaussian distribution. As a result, a noisy image with a SNR of 3.5 dB is obtained.

The hypothetical noisy image undergoes denoising with wavelet Decomposition Levels (DLs) ranging from 1 to 6 and using Dmey wavelet base, resulting in the corresponding Method Noise (MN) and its Probability Density Function (PDF), as illustrated in Figure 5. As depicted in Figure 5a, the denoised image is notably blurred when DL is 1, and its MN (Figure 5b) solely contains noise components, with a PDF distribution closely resembling a Gaussian distribution (Figure 5c). With an increase in DL, clarity is restored in the images (Figure 5d,g), and their corresponding MNs (Figure 5e,h) exhibit no evident periodic signal components, maintaining PDF distributions akin to Gaussian norms. However, beyond DL 3, the denoised images progressively blur compared to Figure 5g (Figure 5j,m,p). Concurrently, the corresponding MNs reveal clear periodic signal components (Figure 5k,n,q), and their PDF distributions distinctly deviate from the Gaussian norm (Figure 5l,o,r). Consequently, the optimal wavelet DL for the denoised image in this scenario is determined to be 3.

To quantitatively assess the optimal DL, the Structural Similarity Index (SSIM) between the PDFs of MN obtained from each DL and the one obtained from the first DL is computed. Specifically, the SSIM for each PDF of Figure 5c,f,i,l,o,r compared to the PDF of Figure 5c is calculated, as depicted in Figure 6a. Using the PDF from the first DL as a reference—given its exclusive inclusion of noise components—reveals a significant decrease in the difference between the third and fourth DL. This signifies that the PDF distribution of the fourth DL deviates significantly from the Gaussian distribution, with the MN evidently containing signal components Figure 5k. Simplifying the determination of the optimal DL, the absolute difference in the SSIMs of adjacent DLs is plotted against DL in Figure 6b. The peak value (red asterisk) in Figure 6b indicates that, as DL increases from 3 to 4, MN contains the most signal components, implying a notable loss of signal components in the denoised image. Thus, the optimal DL is identified as 3. Furthermore, changes in Peak SNR (PSNR) and SSIM with DL are graphed in Figure 6c and Figure 6d, respectively, revealing that the optimal DL determined by our proposed algorithm aligns with those determined by traditional methods (red asterisk in Figure 6c,d) and the visual judgment mentioned earlier. Notably, the calculation of PSNR and SSIM typically necessitates a clean signal or known noise for reference, which is often impractical. In contrast, our proposed algorithm surpasses traditional methods by eliminating the requirement for clean signals or known noise, enhancing its practical utility.

3.2. Measured $\phi$-OTDR Signal

To validate the efficacy of the proposed method, it is implemented on an actual $\phi$-OTDR vibration signal for wavelet noise reduction. The $\phi$-OTDR signal was measured over 100 m at the termination point of a $\phi$-OTDR sensing fiber, located approximately 8 km away, as documented in our previous work [28]. Some backscattering light traces were collected to verify the parameter optimization method based on ACF. The 100 adjacent traces were subtracted to obtain the vibration signals, as depicted in Figure 7a, where the horizontal coordinate is the distance, the vertical coordinate is the trace, and the different colors of the pixels represent different vibration intensities.

In practice, SSIM and PSNR cannot be used to optimize DL due to the lack of clean signals to use as a reference and the difficulty of accurately estimating noise levels. As shown in Figure 7b, FCL-ACF, the signal-based method proposed in our previous work [33] also fails to adaptively optimize DL due to the lack of significant peaks. It is known that a higher FCL-ACF indicates a stronger periodicity between texture signals, which usually means a better denoising effect. However, signal components of the denoised image may be removed as DL increases, while FCL-ACF still increases.

WD is applied to the measured signal using the bior2.2 wavelet base, employing an autocorrelation-based thresholding method [33], across DLs 1 to 6. The resulting denoised images for each DL are presented in the initial column of Figure 8, while their MNs and respective PDFs are displayed in the subsequent two columns of Figure 8. Visual judgment based on denoised images or MNs alone proves challenging for determining the optimal DL. Moreover, the unavailability of PSNR and SSIM in practical scenarios renders traditional evaluation methods impractical. However, a notable deviation from the Gaussian distribution in the PDF of the fifth DL suggests that the optimal DL may be 4. To quantitatively assess this, SSIM between each DL’s PDF and the PDF from the first DL is computed to gauge the extent of deviation from the Gaussian distribution. As DL increases, the PDF distribution of the MN gradually deviates from the Gaussian distribution, which means that the component of the MN is not the pure noise, as depicted in Figure 9a. Since the MN of some neighboring DLs may be close to pure noise, their PDF distribution is more similar to a Gaussian distribution. This is why it may increase slightly, as shown in point 7 of the horizontal coordinate in Figure 9a. However, the selection of the optimal DL depends on the amount of change in SSIM, so this slight upward trend does not affect the selection. Plotting the absolute difference between adjacent SSIMs in Figure 9b reveals the fourth DL as the optimal DL (indicated by a red asterisk). Additionally, the denoised image of the fourth DL exhibits superior clarity compared to other DLs, and its corresponding MN displays a more uniformly distributed noise without conspicuous signal components. Consequently, our proposed algorithm identifies the optimal DL as 4.

4. Discussion

Traditional methods for determining the optimal wavelet DL exhibit shortcomings, including the subjective nature of visual inspection, limitations of co-occurrence matrices for image texture detection, and the inability of noise-energy-based methods to accurately measure noise changes in practical scenarios [39,40,43]. As reported in previous studies, the entropy-based evaluation metrics are similar to FCL-ACF metrics in that they vary monotonically with DL [41,42], and the change curve usually tends to stabilize and continue to grow after reaching the optimal value. However, our proposed method can easily find the maximum value, and thus, adaptively select the optimal DL for WD. More importantly, different from the traditional methods, our method does not require precise knowledge of the noise characteristics or the clean signals, and thus, it is generally applicable to DL optimization for WD denoising of various types of noise-containing signals. According to the evaluation results of FCL-ACF, listed in

Table 1. Comparison of time required for WD with different DLs.

DL	1	2	3	4	5	6	7
FCL-ACF	0.3470	0.4142	0.4163	0.4359	0.4452	0.4755	0.5094
Time (s)	1.0824	2.2738	3.6594	5.0467	6.3418	7.8427	9.3678

, when the fourth layer is selected as the most suitable DL, the FCL-ACF reaches 0.4359, which is nearly 25.62% higher than the 0.3470 obtained when the first layer is less effective in denoising. The noise is not correlated, while the signal is. Therefore, by choosing the appropriate DL, the autocorrelation of the denoised signal is significantly improved, indicating that the noise is effectively filtered. In addition, the processing time for the image denoising with WD in different DLs are also listed in Table 1, which is almost linearly related to DL. The higher the DL, the longer the time required. Therefore, finding the optimal DL quickly and accurately can save processing time. In this case, the computation time with an optimal DL of 4 is nearly 46.13% less than that with DL of 7, which greatly improves the denoising efficiency of WD. MATLAB R2023a software and a computer with a 13th generation Intel(R) Core(TM) i9-13905H 2.6 GHz processor and 32 GB of RAM were used for this calculation.

To assess the impact of Wavelet Denoising (WD) on spatial resolution and signal fading, changes in spatial trace and the time series signal are examined, as illustrated in Figure 10a. The 100 m trace featuring the maximum signal strength before and after undergoing four-layer wavelet decomposition is depicted in Figure 10b. Notably, the denoised trace (red dashed line) closely aligns with the raw trace (blue solid line), exhibiting only a slight decrease in peak value. The signal components within the initial 40-m sensing fiber remain intact, while noise components in the subsequent 60 m are noticeably attenuated. The full width at half maximum of the trace waveform with the highest intensity remains at 7 m before and after WD, affirming that spatial resolution remains unaltered. Examining the difference in time series signal (green dotted line) between the raw signal (blue solid line) and the denoised signal (red dashed line) in Figure 10c reveals that peak values are generally consistent, indicating that WD may not exacerbate signal fading. For a comprehensive analysis of signal distortion, the standard deviation (STD) of Method Noise (MN) at various sensing points along the optical fiber is computed, as shown in Figure 10d. The STD of the initial 40-m fiber optic sensing signal during vibration and the STD of the subsequent 60 m at rest both fall below 0.121. This suggests that the fluctuation range of MN is minimal, and WD does not introduce significant signal distortion.

5. Conclusions

We propose an MN-based signal analysis method that addresses the difficulty of wavelet DL selection, relying on the resemblance between the PDF distribution of MN and Gaussian noise distribution. Validation against traditional PSNR, SSIM methods, and signal-based FCL-ACF methods, using both an assumed noisy image and a measured $\phi$-OTDR signal, confirmed the superiority of our proposed method. The results show that using the MN-based method, significant signal components can be observed in the MN of a signal when the wavelet DL of the hypothetical signal exceeds the optimal value of 3. Similarly, when the optimal wavelet DL exceeds 4, this phenomenon is clearly observed by analysing the MN of the measured $\phi$-OTDR vibration signal. Moreover, not limited to periodic signals, the optimal wavelet DLs for all types of signals can be found efficiently by this method. By choosing the appropriate DL, the time required to process the signal is reduced by almost 46.13%. This research contributes to enhancing measurement accuracy and spatial resolution in various distributed fiber optic sensing systems, augmenting their signal pattern recognition capabilities.