A method for extracting effective Raman spectra of natural gas based on optimal wavelet decomposition layer number

Abstract

The application discloses a method for extracting natural gas effective Raman spectrum based on a wavelet optimal decomposition layer number, and comprises the following steps: carrying out wavelet transform on a Raman original spectrum signal f(t), carrying out j-layer and j+1-layer wavelet decomposition, combining a threshold function, judging whether ∑(D j+1 2 is greater than or equal to ∑cD j 2 , determining the current j value as the optimal decomposition layer number, and restoring the natural gas effective Raman spectrum. The application directly determines the optimal decomposition layer number by comparing the energy of each layer wavelet detail coefficient, and does not need to obtain a denoising signal through signal reconstruction, and then compare performance indexes such as a signal-to-noise ratio and a mean square error of multiple denoising results to determine the optimal decomposition layer number. The method has an accurate process result, a more simple process and a shorter time.

Description

Technical Field

[0001] This invention belongs to the field of natural gas spectral analysis, specifically relating to a method for extracting effective Raman spectra of natural gas based on the optimal number of wavelet decomposition layers. Background Technology

[0002] Natural gas spectral signals obtained by optical instruments contain a large amount of noise. Removing noise through effective preprocessing methods is a necessary prerequisite for ensuring the accuracy of the final analysis results, especially under low signal-to-noise ratio conditions.

[0003] Traditional signal denoising mainly relies on Fourier transform and smoothing window methods:

[0004] 1. The moving smoothing window method, based on the difference between the statistical characteristics of signal and noise, assumes that the noise is zero-mean noise. It aims to improve the signal-to-noise ratio by averaging the original signal. By selecting a smoothing window with an odd width of 2m+1, and using the center wavelength point k as a reference, the window is moved from left to right. The average value of all measurements within the window's coverage area replaces the measurement value corresponding to the center wavelength point until all points are smoothed. The width m of the smoothing window significantly affects the smoothing result; too small a width results in poor smoothing, while too large a width removes characteristic peak information, causing spectral distortion. Furthermore, there is a boundary problem: for a smoothing window with a width of 2m+1, m points at each end of the spectrum cannot be processed.

[0005] 2. The Fourier transform method involves performing a Fourier transform on the noisy signal to convert it from the time domain to the frequency domain. Then, certain frequency bands in the frequency domain are filtered out. Finally, the denoised signal is obtained through an inverse Fourier transform. For stationary signals, Fourier transform-based denoising methods have excellent denoising effects, but their effectiveness is poor for non-stationary signals.

[0006] Traditional denoising methods increase signal entropy and fail to describe the non-stationary characteristics and correlations of signals. Wavelet thresholding denoising offers advantages such as speed, accuracy, and flexible basis selection. Furthermore, wavelet transform exhibits excellent time-frequency characteristics and performs well in analyzing non-stationary signals and removing signal correlations, making it widely used in signal denoising.

[0007] When using wavelet analysis for denoising, insufficient decomposition levels result in a significant amount of noise remaining in the denoised signal; conversely, excessive decomposition levels lead to the loss of valuable information. To remove as much noise as possible from the original spectrum while preserving the complete effective information, existing methods determine the optimal number of decomposition levels through empirical judgment and comparative experiments. However, this method is heavily influenced by human factors and requires highly experienced operators. Especially when the optimal number of decomposition levels is large, this method becomes inefficient, cumbersome, and the determined decomposition level may not even be optimal. Summary of the Invention

[0008] In view of this, the present invention addresses the lack of a good method for determining the optimal number of decomposition layers in the current wavelet thresholding denoising process, and provides a method for extracting the effective Raman spectrum of natural gas based on the optimal number of wavelet decomposition layers, so as to solve the problems of low efficiency and cumbersome process in determining the optimal number of decomposition layers in the existing technology.

[0009] The technical solution is as follows: a method for extracting the effective Raman spectrum of natural gas based on the optimal wavelet decomposition layer number, the key of which includes the following steps:

[0010] S1. Acquire the raw Raman spectral signal f(t) of natural gas using optical instruments;

[0011] S2. After performing wavelet transform on the original Raman spectral signal f(t), substitute its wavelet energy coefficients into the expression for detail coefficient energy and total noise energy:

[0012]

[0013] E j For the wavelet detail coefficient energy under the j-level decomposition, cD j，i Let be the wavelet detail coefficients under the j-level decomposition, and n be the number of wavelet detail coefficients under the j-level decomposition.

[0014] After performing j-level wavelet decomposition on the raw Raman spectral signal f(t), the energy ∑cD of the j-level wavelet detail coefficients is obtained. j 2 In the formula, E represents the total noise energy in the original Raman spectral signal f(t). j For the wavelet detail coefficient energy under the j-level decomposition, cD j，i Let be the wavelet detail coefficients under level j decomposition, n be the number of wavelet detail coefficients under level j decomposition, and j be the number of wavelet decomposition levels. The initial value is equal to 1.

[0015] S3. Perform wavelet decomposition on the raw Raman spectral signal f(t) at level j+1 to obtain the energy ∑cD of the wavelet detail coefficients at level j+1. j+12 ;

[0016] S4. Determine ∑cD j+1 2 Is it greater than or equal to ∑cD? j 2 If not, then j = j + 1, return to S2 and S3, and repeat the loop until ∑cD j+1 2 ≥∑cD j 2 ;

[0017] By observing the wavelet detail coefficients at different decomposition levels, it can be found that at low decomposition levels, the wavelet detail coefficients only contain noise information; as the decomposition level increases, the wavelet detail coefficients begin to contain effective information; and finally, the wavelet detail coefficients contain a large amount of effective information. Observing the approximation coefficient transformation also reveals that as the decomposition level increases, the effective information in the signal gradually decreases.

[0018] From the perspective of signal energy, at low decomposition levels, the energy of noise signals in wavelet detail coefficients decreases as the decomposition level increases. However, as the level continues to increase, the detail coefficients begin to contain the energy of effective signals, which is much greater than that of noise signals. At this point, the energy contained in the detail coefficients increases more and more randomly with the increase of the decomposition level. The change in the energy trend of wavelet coefficients is also a sign that wavelet threshold denoising has reached the optimal decomposition scale.

[0019] S5. Determine the current j value as the optimal number of layers, perform j-layer wavelet decomposition on the original Raman spectrum signal to obtain the j-layer wavelet coefficients, and combine the threshold function to perform threshold processing on the wavelet coefficients. Finally, perform inverse wavelet transform to convert the current frequency domain signal into a time domain signal to restore the effective Raman spectrum of natural gas.

[0020] In steps S2 and S3, the original Raman spectral signal f(t) is transformed by layer-by-layer wavelet transform to obtain wavelet coefficients w. j，i Wavelet coefficients w layer by layer j，i It can be divided into wavelet detail coefficients cD j，i and wavelet approximation coefficient cA j，i Until step S5, a threshold function is used to perform thresholding on both to obtain the thresholded wavelet coefficients.

[0021] The wavelet transform formula is as follows:

[0022]

[0023] Here, is the wavelet basis function, a is the scaling parameter controlling the scaling, b is the position delay (translation) parameter controlling the position delay, and t is the Raman shift of the Raman spectrum abscissa.

[0024]

[0025] In the formula, WT f (m, n) represents the discrete wavelet transform of the original signal f(t), signifying a process; the wavelet coefficients obtained after the discrete wavelet transform of the original signal f(t) are denoted as w. j，i j represents the number of decomposition levels j, and i represents the i-th original element;

[0026] m and n are discrete values. The scale parameter a is discretized using the first discrete value m, and the scale parameter b is discretized using the second discrete value n. Here, only the discrete wavelet transform is emphasized because the discrete wavelet transform is actually used.

[0027] w j，i Divided into wavelet detail coefficients cD j，i and wavelet approximation coefficient cA j，i .

[0028] In step S5, the threshold function is a fixed threshold sqtwolog, and the expression for its threshold λ is:

[0029]

[0030]

[0031] Where σ is the standard deviation of the noise, and N is the wavelet coefficient w. j，i The quantity.

[0032] In step S5, the threshold function is the unbiased likelihood estimation criterion rigrsure threshold, used to estimate the wavelet coefficients w of the threshold. j，i Take the absolute values, sort them in ascending order, then square each element to obtain a new vector S to be estimated. For w... j，i For each element, calculate the risk vector using the following formula:

[0033]

[0034] Risk k Let N represent the risk vector of the k-th element, and let N be the wavelet coefficients w. j，i Quantity;

[0035] Find the minimum Risk in the risk vector. k As the minimum risk value, let k be the index corresponding to the minimum point of the risk vector Risk, thus obtaining the threshold:

[0036] In step S5, the threshold function is a heuristic threshold, which is a combination of the heuristic threshold and the sqtwolog threshold. Specifically, this involves comparing two variables, β and γ.

[0037]

[0038]

[0039] In the formula, w j，i For wavelet coefficients, if β is less than γ, the threshold T is the Sqtwolog threshold; if β is greater than or equal to γ, the threshold λ is the smaller of the Rigrsure threshold and the Sqtwolog threshold.

[0040] In step S5, the threshold function is a maximum-minimum threshold (minmax), and its expression is:

[0041]

[0042] Where N is w j，i The number of wavelet coefficients. This principle is used in statistics for designing estimators. Since the denoised signal can be assumed to be an estimator of an unknown regression function, the minimax estimator is the quantity that minimizes the mean square error under the worst-case condition, i.e., its selected threshold is the one that produces the minimum maximum variance.

[0043] Compared with existing technologies, this invention directly determines the optimal number of decomposition layers by comparing the energy of wavelet detail coefficients at each layer. It does not require obtaining a denoised signal through signal reconstruction and then comparing the signal-to-noise ratio and mean square error of multiple denoising results to determine the optimal number of decomposition layers. The method in this paper is more accurate, simpler, and faster. Attached Figure Description

[0044] Figure 1 This is a flowchart of the present invention;

[0045] Figure 2 The waveform of the original Raman spectrum signal;

[0046] Figure 3 A hierarchical diagram of wavelet detail coefficient decomposition;

[0047] Figure 4 This is a graph showing the energy variation of the detail coefficients;

[0048] Figure 5 The wavelet decomposition of the original signal yields a graph of detail coefficients and approximation coefficients.

[0049] Figure 6The waveform diagram shows the result of noise reduction processing for the noise signal. Detailed Implementation

[0050] The present invention will be further described below with reference to the embodiments and accompanying drawings.

[0051] Example 1, such as Figure 1 As shown, the method for extracting the effective Raman spectrum of natural gas based on the optimal wavelet decomposition layer number includes the following steps:

[0052] like Figure 2 As shown, S1, the raw Raman spectral signal f(t) of natural gas is acquired by optical instruments;

[0053] S2. After performing wavelet transform on the original Raman spectral signal f(t), substitute its wavelet energy coefficients into the expression for detail coefficient energy and total noise energy.

[0054] like Figure 3 , 4 As shown in Figure 5, the power spectral density of white noise is constant, meaning the energy is equal per unit frequency. In other words, the total energy of the white noise signal is the same constant at all frequencies. Combining this with the principle that the signal is equally divided into high and low frequencies during multi-scale wavelet decomposition, and considering that the frequency range of the wavelet detail coefficients gradually decreases by 50% as the decomposition level increases, it can be seen that when white noise undergoes multi-level wavelet decomposition, the energy in its detail coefficients decreases by 50% with each increase in the decomposition scale. Therefore, the relationship between the energy of the detail coefficients at each scale and the total noise energy is:

[0055]

[0056] E j For the wavelet detail coefficient energy under the j-level decomposition, cD j，i Let be the wavelet detail coefficients under the j-level decomposition, and n be the number of wavelet detail coefficients under the j-level decomposition.

[0057] After performing j-level wavelet decomposition on the raw Raman spectral signal f(t), the energy ∑cD of the j-level wavelet detail coefficients is obtained. j 2 In the formula, E represents the total noise energy in the original Raman spectral signal f(t). j For the wavelet detail coefficient energy under the j-level decomposition, cD j，i Let be the wavelet detail coefficients under level j decomposition, n be the number of wavelet detail coefficients under level j decomposition, and j be the number of wavelet decomposition levels. The initial value is equal to 1.

[0058] S3. Perform wavelet decomposition on the raw Raman spectral signal f(t) at level j+1 to obtain the energy ∑cD of the wavelet detail coefficients at level j+1. j+1 2 ;

[0059] S4. Determine ∑cD j+1 2 Is it greater than or equal to ∑cD? j 2 If not, then j = j + 1, return to S2 and S3, and repeat the loop until ∑cD j+1 2 ≥∑cD j 2 ;

[0060] like Figure 6 As shown in step S5, determine the current j value as the optimal number of layers, perform j-layer wavelet decomposition on the original Raman spectrum signal to obtain the j-layer wavelet coefficients, and combine the threshold function to perform threshold processing on the wavelet coefficients. Finally, perform inverse wavelet transform to convert the current frequency domain signal into a time domain signal to restore the effective Raman spectrum of natural gas.

[0061] In steps S2 and S3, the original Raman spectral signal f(t) is transformed by layer-by-layer wavelet transform to obtain wavelet coefficients w. j，i Layer-by-layer wavelet coefficients w j，i It can be divided into wavelet detail coefficients cD j，i and wavelet approximation coefficient cA j，i Until step S5, a threshold function is used to perform thresholding on both to obtain the thresholded wavelet coefficients.

[0062] The wavelet transform formula is as follows:

[0063]

[0064] Here, is the wavelet basis function, a is the scaling parameter controlling the scaling, b is the position delay (translation) parameter controlling the position delay, and t is the Raman shift of the Raman spectrum abscissa.

[0065]

[0066] In the formula, WT f (m, n) denotes the discrete wavelet transform of the original signal f(t); the wavelet coefficients obtained after the discrete wavelet transform of the original signal f(t) are denoted as w. j，i j represents the number of decomposition levels j, and i represents the i-th original element;

[0067] m and n are discrete values. The scale parameter a is discretized using the first discrete value m, and the scale parameter b is discretized using the second discrete value n.

[0068] w j，i Divided into wavelet detail coefficients cD j，i and wavelet approximation coefficient cAj，i .

[0069] In Example 2, in step S5, the threshold function is a fixed threshold sqtwolog, and the expression for its threshold λ is:

[0070]

[0071]

[0072] Where σ is the standard deviation of the noise, and N is the wavelet coefficient w. j，i The quantity.

[0073] In step S5, the threshold function is the unbiased likelihood estimation criterion rigrsure threshold, used to estimate the wavelet coefficients w of the threshold. j，i Take the absolute values, sort them in ascending order, then square each element to obtain a new vector S to be estimated. For w... j，i For each element, calculate the risk vector using the following formula:

[0074]

[0075] Risk k Let N represent the risk vector of the k-th element, and let N be the wavelet coefficients w. j，i Quantity;

[0076] Find the minimum Risk in the risk vector. k As the minimum risk value, let k be the index corresponding to the minimum point of the risk vector Risk, thus obtaining the threshold:

[0077] In Example 3, in step S5, the threshold function is a heuristic threshold, which is a combination of the rigrsure threshold and the sqtwolog threshold. Specifically, this involves comparing two variables, β and γ.

[0078]

[0079]

[0080] In the formula, w j，i For wavelet coefficients, if β is less than γ, the threshold T is the Sqtwolog threshold; if β is greater than or equal to γ, the threshold λ is the smaller of the Rigrsure threshold and the Sqtwolog threshold.

[0081] In Example 4, in step S5, the threshold function is a minima threshold minmax, and its expression is:

[0082]

[0083] Where N is w j，i The number of wavelet coefficients.

[0084] Compared to existing methods, this method directly determines the optimal number of decomposition layers by comparing the energy of wavelet detail coefficients at each layer. It eliminates the need for signal reconstruction to obtain a denoised signal, and then compares the signal-to-noise ratio and mean square error of multiple denoising results to determine the optimal number of decomposition layers. This method is more accurate, simpler, and faster. The accuracy of this method is shown in Tables 1-4.

[0085] Table 1. Noise ratio of block signals at different scales

[0086]

[0087] Table 2. Noise ratio of bump signals at different scales

[0088]

[0089] Table 3. Noise ratio of heavy sine signal at different scales

[0090]

[0091] Table 4. Noise ratio of Doppler signal at different scales

[0092]

[0093] Four signals—blocks, bumps, heavy sine, and doppler—were selected from Matlab and white noise was added as the original signals. Thresholds were selected using heuristic thresholding, sqtwolog thresholding, and minimaxi thresholding, respectively. A soft thresholding method was used with the sym8 wavelet basis. For the four original signals, the optimal decomposition level K calculated by the proposed method was 3, 3, 6, and 4, respectively. When using heursure and minimaxi thresholding, the proposed method accurately found the optimal decomposition level. When using sqtwolog thresholding, the optimal level calculated by the proposed method deviated from the actual optimal level by one level for the blocks, bumps, and heavy sine signals. However, the signal-to-noise ratio difference between the optimal scale denoising result determined by the proposed method and the denoising result with the actual optimal level was very small.

[0094] The noise reduction effect of four simulated signals under the Minimaxi threshold selection rule, using the decomposition scale determined by this method, is as follows: Figure 6As shown in the figure, the waveform is well recovered while improving the signal-to-noise ratio. Finally, it should be noted that the above description is merely a preferred embodiment of the present invention. Those skilled in the art, under the guidance of the present invention, can make various similar representations without departing from the spirit and claims of the present invention, and all such modifications fall within the protection scope of the present invention.

Claims

1. A method for extracting effective Raman spectrum of natural gas based on the optimal decomposition level of wavelet, characterized in that, The method comprises the following steps: S1, collecting a Raman original spectrum signal f(t) of natural gas by optical instruments; S2, after wavelet transformation of the Raman original spectrum signal f(t), substituting wavelet energy coefficients into an expression of detail coefficient energy and total noise energy: E j To j cDj-1is the energy of the wavelet detail coefficients at layer j-1, n is the number of wavelet detail coefficients at layer j-1. j，i cDjis the energy of the wavelet detail coefficients at layer j, n is the number of wavelet detail coefficients at layer j. After performing j-level wavelet decomposition on the raw Raman spectral signal f(t), the energy ∑cD of the j-level wavelet detail coefficients is obtained. j 2 In the formula, E represents the total noise energy in the original Raman spectral signal f(t). j To j Energy of wavelet detail coefficients under layer decomposition, cD j，i Let be the wavelet detail coefficients under the j-level decomposition, n be the number of wavelet detail coefficients under the j-level decomposition, j be the number of wavelet decomposition levels, the initial value is equal to 1, and i represent the i-th wavelet coefficient; S3, performing j+1 layer wavelet decomposition on the Raman original spectrum signal f(t), and obtaining j+1 layer wavelet detail coefficient energy ∑cD j+1 2 ; S4, judge ∑cD j+1 2 whether greater than or equal to ∑cD j 2 , no, then j = j + 1, return to S2 and S3, loop until ∑cD j+1 2 ≥∑cD j 2 ; S5, determining a current j value as an optimal decomposition layer number, performing j-layer wavelet decomposition on the original Raman spectrum signal to obtain j-layer wavelet coefficients, performing threshold processing on the wavelet coefficients in combination with a threshold function, and finally inversely transforming a current frequency domain signal into a time domain signal by wavelet inverse transformation to restore an effective Raman spectrum of natural gas.

2. The method for extracting effective Raman spectrum of natural gas based on the optimal decomposition level of wavelet according to claim 1, characterized in that, In the step S2 and the step S3, the Raman original spectrum signal f(t) is changed by the wavelet to obtain the wavelet coefficient w j，i The wavelet coefficient w j，i The wavelet coefficient w is decomposed layer by layer to obtain the wavelet detail coefficient cD of the layer j，i And the wavelet approximation coefficient cA j，i After that, until the step S5, the threshold function is used to threshold process the two to obtain the threshold processed wavelet coefficient 3. The method for extracting effective Raman spectrum of natural gas based on the optimal decomposition level of wavelet according to claim 2, characterized in that, The wavelet transformation formula is: where w is a wavelet basis function, a is a scale parameter that controls the stretching, b is a position delay (translation) parameter that controls the position, and t is the Raman spectrum abscissa Raman shift. wherein WT f (m, n) denotes the discrete wavelet transform of the original signal f(t); the original signal f(t) is transformed by the discrete wavelet transform to obtain a series of wavelet coefficients, denoted as w j，i j denotes the decomposition j level. m and n are discrete values, a scale parameter a is discretized by the first discrete value m, and a scale parameter b is discretized by the second discrete value n; w j，i cAand cDare the wavelet approximation and detail coefficients, respectively. j，i j，i cAand cDare the wavelet approximation and detail coefficients, respectively.​ 4. The method for extracting effective Raman spectrum of natural gas based on the optimal decomposition level of wavelet according to claim 2, characterized in that, In the step S5, the threshold function is a fixed threshold sqtwolog, and an expression of a threshold value λ of the threshold function is: where σ is the standard deviation of the noise, N is the number of wavelet coefficients w j，i the number of the quantity.

5. The method for extracting effective Raman spectrum of natural gas based on the optimal decomposition level of wavelet according to claim 2, characterized in that, In the step S5, the threshold function is a rigrsure threshold of unbiased likelihood estimation criterion, which is used to estimate the threshold of the wavelet coefficient w j，i Taking absolute value, sorting from small to large, then squaring each element, a new to-be-estimated vector S is obtained, and the risk vector is calculated according to the following formula for each element of w j，i ​ Risk k Let N represent the risk vector of the k-th element, and let N be the wavelet coefficients w. j，i Quantity; Find the minimum Risk in the risk vector k As the lowest risk value, record the index value corresponding to the minimum point of the risk vector Risk, so that the threshold is:

6. The method for extracting effective Raman spectrum of natural gas based on optimal decomposition level of wavelet according to claim 2 or 4 or 5, characterized in that, In the step S5, the threshold function is a heuristic threshold heursure, which is a comprehensive form of a rigrsure threshold and a sqtwolog threshold, and a specific form is to compare two variables β and γ first, where w j，i The threshold T takes the Sqtwolog threshold if β is less than γ, and takes the smaller of the Rigrsure threshold and the Sqtwolog threshold if β is greater than or equal to γ.

7. The method for extracting effective Raman spectrum of natural gas based on optimal decomposition level of wavelet according to claim 1, characterized in that, In the step S5, the threshold function is a max-min threshold minmax, and an expression of the threshold function is: where N is w j，i The number of wavelet coefficients.

Images