GNSS tunnel portal slope deformation time sequence noise reduction method and system combined with CEEMD and adaptive wavelet packet threshold function, and medium

Abstract

The invention relates to the cross technical field of geological disaster monitoring and signal processing, and provides a GNSS tunnel portal slope deformation time sequence noise reduction method, system and medium combining CEEMD and an adaptive wavelet packet threshold function, the method integrates the advantages of CEEMD and wavelet packet adaptation, processing of mixed entropy demarcation, adaptive threshold, dynamic smoothing and multi-path compensation is used, and the noise reduction efficiency is improved. Useful deformation features can be reserved while more comprehensive and finer noise reduction is performed, smooth compression can be performed in a large signal part through improved threshold function processing, excessive suppression can be avoided, smooth transition is performed in a small signal part through exponential attenuation, features of the large signal part are reserved, noise suppression is realized, signal details are reserved to the greatest extent, and the noise suppression effect is improved. The noise reduction requirement of tunnel portal side slope monitoring data is met, and particularly in side slope monitoring of a tunnel portal in a complex environment, the precision and reliability of GNSS deformation monitoring are improved.

Description

Technical Field

[0001] The present invention relates to the technical field of intersection of geological disaster monitoring and signal processing, and in particular to a GNSS deformation monitoring time series noise reduction technology for tunnel vibration and complex multipath noise scenarios. A GNSS tunnel entrance slope deformation time series noise reduction method, system and medium combining CEEMD and an adaptive wavelet packet threshold function are proposed. Background Art

[0002] Frequent geological disasters pose a serious threat to human society. Landslides, a typical example, are crucial for monitoring and early warning, ensuring they are crucial for reducing losses. The application of Global Navigation Satellite System (GNSS) technology in landslide monitoring has become a research hotspot. GNSS technology, with its high-precision, all-weather, continuous three-dimensional positioning capabilities, and the absence of line-of-sight, provides strong support for accurately acquiring landslide deformation information.

[0003] In practical landslide monitoring applications, coordinate time series acquired by GNSS are often subject to various interferences and noises, which can obscure the true landslide signal and reduce the accuracy of the monitoring data. Among these interferences, multipath is the primary factor affecting deformation monitoring. Multipath is caused by signal reflection or diffraction due to the influence of the receiver's surrounding environment.

[0004] Current research on multipath effects primarily focuses on high-precision GNSS processing algorithms and improvements to receiver hardware. Receiver hardware improvements include the use of antenna methods, multipath estimation techniques (MET), and multipath estimation delay-locked loops (MEDLLs). While hardware improvements can mitigate some multipath errors, residual errors remain a major source of error in deformation monitoring. GNSS data post-processing algorithms include multi-antenna array technology and stellar day filtering. Research on stellar filtering methods, based on the spatiotemporal repeatability of satellites, primarily encompasses both the coordinate and observation domains, and has been applied to studies that ignore satellite effects. For example, Genrich and Bock first proposed the stellar filtering (SF) algorithm to reduce multipath errors at observatory stations. Since then, research on SF algorithms has deepened, including the development of SF algorithms in the coordinate domain, analysis of multipath effects for landslide monitoring in complex environments, and improved monitoring accuracy by correcting monitoring sequences in the coordinate domain. Although coordinate-domain methods can correct for the combined multipath errors of all satellites, differences in the orbital repetition periods of different satellites can lead to time-shift errors in the calculated coordinate sequences.

[0005] To address this issue, existing technologies have proposed observation-domain multipath error correction models. These models utilize the repetition period of each satellite to extract multipath errors from single- or double-difference residuals and apply corrections in the next period. Their corrections are superior to coordinate-domain methods. However, these methods do not perform any decomposition and only process GNSS data in a single dimension. Due to the inherent multi-scale nature of GNSS data, these methods are unable to decompose the signal into different time scales.

[0006] In addition, in existing designs or research, most of the GNSS data for multipath effect analysis are collected in field environments or simulated environments. However, there is no relevant existing technical research on landslide monitoring applications in complex environments such as undulating terrain and vegetation cover. Summary of the Invention

[0007] In view of the defects and shortcomings of the existing technology, the purpose of the present invention is to provide a GNSS tunnel entrance slope deformation time series denoising method that combines CEEMD with an adaptive wavelet packet threshold function. The adaptive wavelet packet threshold function is used to identify and remove the noise of the tunnel entrance slope monitoring deformation sequence, thereby realizing the correction capability of the multipath effect and reducing the adverse effects of the multipath effect on the monitoring deformation sequence in the tunnel entrance slope environment.

[0008] According to a first aspect of the present invention, a method for denoising GNSS tunnel entrance slope deformation time series by combining CEEMD with an adaptive wavelet packet threshold function is proposed, comprising the following steps:

[0009] A GNSS tunnel entrance slope deformation time series denoising method combining CEEMD and adaptive wavelet packet threshold function is characterized by comprising the following steps:

[0010] Step 1: Obtain the original GNSS signal for GNSS tunnel entrance slope deformation monitoring;

[0011] Step 2: Add a Gaussian white noise signal with opposite positive and negative values ​​to the original GNSS signal and use the CEEMD algorithm to decompose the signal into multiple IMF components.

[0012] Step 3: Construct the mixed entropy of Pearson correlation coefficient and fuzzy entropy, detect the local minimum based on the mixed entropy, and determine the component IMF corresponding to the first mixed entropy that takes the local minimum. k As the boundary between noise and signal, IMF is determined based on this k The components after sequence number k are dominated by noise, and the noisy signals are filtered out;

[0013] Step 4: For the filtered noisy signal, noise weights are constructed by performing sub-band frequency band energy-entropy joint analysis based on wavelet packet decomposition, and a threshold is adaptively determined based on the sub-band noise level and the noise weight to perform segmented threshold function processing to obtain sub-bands after threshold processing;

[0014] Step 5: Perform inverse wavelet packet transform on the sub-band after threshold processing to obtain the denoised noise IMF component as the reconstructed signal;

[0015] Step 6: dynamically determine the Gaussian kernel width according to the local fluctuation level of the reconstructed signal, perform convolution processing on the reconstructed signal using the Gaussian kernel function, and smooth the denoising result through adaptive Gaussian convolution; and

[0016] Step 7: extract the multipath effect component from the residual obtained by decomposition based on the CEEMD algorithm and design a notch filter accordingly, perform multipath effect filtering compensation on the denoising result after adaptive Gaussian convolution smoothing, and eliminate periodic multipath interference.

[0017] According to a second aspect of the present invention, a computer system is provided, comprising:

[0018] one or more processors; and

[0019] Memory, which stores instructions that can be operated;

[0020] Among them, when the instruction is executed by one or more processors, the one or more processors mentioned above perform operations, and the operations include executing the process of the GNSS tunnel entrance slope deformation time series denoising method combining CEEMD and adaptive wavelet packet threshold function.

[0021] According to a third aspect of the present invention, a computer-readable storage medium is provided for storing one or more programs, wherein the one or more programs include instructions or instruction sets that can be executed by one or more processors;

[0022] When the instructions or instruction set are executed by one or more processors, the process of the aforementioned GNSS tunnel entrance slope deformation time series denoising method combining CEEMD with an adaptive wavelet packet threshold function is performed.

[0023] The GNSS tunnel entrance slope deformation time series denoising method combining CEEMD and adaptive wavelet packet threshold function of the technical solution of the present invention combines the advantages of CEEMD and wavelet packet adaptation, uses mixed entropy demarcation, adaptive threshold, dynamic smoothing and multipath compensation processing to more comprehensively and finely reduce noise while retaining useful signals, especially in the slope monitoring of tunnel entrances in complex environments, which have the characteristics of more noise and non-stationary signals, thereby improving the accuracy and reliability of GNSS deformation monitoring.

[0024] Compared with the existing technology, the GNSS tunnel entrance slope deformation time series denoising method of the present invention has the following significant advantages:

[0025] (1) In modal decomposition, Gaussian white noise with opposite positive and negative noises is introduced, and the noise residue is offset by complementary integrated averaging. In particular, compared with traditional EMD decomposition (which is sensitive to noise and easily produces false modes (modal aliasing), resulting in impure signal decomposition) and EEMD decomposition (although modal aliasing is suppressed by noise injection, residual noise still affects the decomposition accuracy), CEEMD significantly improves the purity of IMF components by canceling positive and negative noises, effectively reducing the modal aliasing problem.

[0026] (2) By constructing a hybrid entropy index of the Pearson correlation coefficient and fuzzy entropy, the correlation (signal-dominated) and complexity (noise-dominated) of the IMF components are comprehensively evaluated, and the boundary between noise and signal is determined by the local minimum. The dual-index joint analysis proposed in this invention enhances the robustness of the boundary discrimination, which is particularly suitable for separating non-stationary noise and weak signals in the tunnel entrance slope deformation signal; by combining CEEMD with hybrid entropy to discriminate the boundary point, the decomposition accuracy of non-stationary noise is improved, and high robustness is achieved;

[0027] (3) In the complex scenario of high-frequency vibration and low-frequency deformation mixed in GNSS monitoring at the tunnel entrance, noise weights are constructed based on the joint analysis of sub-band frequency band energy-entropy, and the threshold is dynamically adjusted in combination with the sub-band noise level to achieve segmented threshold processing, thereby avoiding the loss of high-frequency weak signals or low-frequency noise residue caused by the indistinguishable noise energy differences between frequency bands; through the sub-band energy-entropy weight distribution of the present invention, a stricter threshold is applied to the high-noise sub-band, and the low-noise sub-band retains more effective signals; thereby achieving dynamic adjustment from threshold selection to smoothing kernel width, avoiding the loss of details and noise residue caused by manual intervention and fixed thresholds, adapting to complex scenarios, and improving adaptive capabilities;

[0028] (4) In the process of signal reconstruction, the Gaussian kernel width is dynamically adjusted according to the local fluctuation of the reconstructed signal to achieve adaptive smoothing, thereby avoiding the problem of detail distortion caused by signal mutation areas (such as sudden deformation of the slope). Dynamic smoothing can enhance the denoising effect in the stable area, retain details in the mutation area, and adapt to the local non-stationary characteristics of the slope deformation signal; and after reconstruction, the residual analysis is combined with the notch filter to solve the dynamic interference problem of the multipath effect in GNSS tunnel entrance monitoring. By extracting the multipath effect components of the residual after CEEMD decomposition, a notch filter is designed for targeted compensation to cope with the time-varying characteristics of the multipath effect at the tunnel entrance. The multipath period is extracted in real time through residual analysis, and the notch parameters are dynamically adjusted. Based on the residual adaptive filtering, the suppression capability of periodic interference is significantly improved.

[0029] It should be understood that all combinations of the foregoing concepts and the additional concepts described in more detail below, as long as such concepts are not mutually inconsistent, can be considered part of the inventive subject matter of this disclosure. In addition, all combinations of the claimed subject matter are considered part of the inventive subject matter of this disclosure.

[0030] The foregoing and other aspects, embodiments, and features of the present invention will be more fully understood from the following description in conjunction with the accompanying drawings. Other additional aspects of the present invention, such as features and / or beneficial effects of the exemplary embodiments, will become apparent from the following description or through practice of specific embodiments according to the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0031] The accompanying drawings are not intended to be drawn to scale. In the accompanying drawings, each identical or nearly identical component shown in various figures may be represented by the same reference numeral. For clarity, not every component is labeled in every figure. Embodiments of various aspects of the present invention will now be described by way of example and with reference to the accompanying drawings.

[0032] Figure 1 3 is a flow chart of a GNSS tunnel entrance slope deformation time series denoising method combining CEEMD with an adaptive wavelet packet threshold function according to an embodiment of the present invention.

[0033] Figure 2 3 is a comparison diagram of threshold functions of a GNSS tunnel entrance slope deformation time series denoising method combining CEEMD with an adaptive wavelet packet threshold function according to an embodiment of the present invention.

[0034] Figure 3 1 is a schematic diagram of four analog signal waveforms according to an embodiment of the present invention.

[0035] Figure 4 3 is a comparison diagram of the results of four analog signals processed by different denoising methods according to an embodiment of the present invention.

[0036] Figure 5 3 is a schematic diagram comparing the RMSE accuracies of four analog signals processed by different denoising methods according to an embodiment of the present invention.

[0037] Figure 6 3 is a schematic diagram comparing the SNR accuracies of four analog signals processed by different denoising methods according to an embodiment of the present invention. DETAILED DESCRIPTION

[0038] In order to better understand the technical content of the present invention, specific embodiments are given and described below with reference to the accompanying drawings.

[0039] Various aspects of the present invention are described in this disclosure with reference to the accompanying drawings, in which a number of illustrative embodiments are shown. The embodiments of the present disclosure are not necessarily intended to include all aspects of the present invention. It should be understood that the various concepts and embodiments introduced above, as well as those described in more detail below, can be implemented in any of many ways, because the concepts and embodiments disclosed herein are not limited to any embodiment. In addition, some aspects of the present disclosure may be used alone or in any appropriate combination with other aspects disclosed herein.

[0040] {Example 1}

[0041] Combined with attachment Figure 1 As shown, according to an embodiment of the present invention, a GNSS tunnel entrance slope deformation time series denoising method combining CEEMD and an adaptive wavelet packet threshold function includes the following steps:

[0042] Step 1: Obtain the original GNSS signal for GNSS tunnel entrance slope deformation monitoring;

[0043] Step 2: Add a Gaussian white noise signal with opposite positive and negative values ​​to the original GNSS signal and use the CEEMD algorithm to decompose the signal into multiple IMF components.

[0044] Step 3: Construct the mixed entropy of Pearson correlation coefficient and fuzzy entropy, detect the local minimum based on the mixed entropy, and determine the component IMF corresponding to the first mixed entropy that takes the local minimum. k As the boundary between noise and signal, IMF is determined based on this k The components after sequence number k are dominated by noise, and the noisy signals are filtered out;

[0045] Step 4: For the filtered noisy signal, noise weights are constructed by performing sub-band frequency band energy-entropy joint analysis based on wavelet packet decomposition, and a threshold is adaptively determined based on the sub-band noise level and the noise weight to perform segmented threshold function processing to obtain sub-bands after threshold processing;

[0046] Step 5: Perform inverse wavelet packet transform on the sub-band after threshold processing to obtain the denoised noise IMF component as the reconstructed signal;

[0047] Step 6: Dynamically determine the Gaussian kernel width according to the local fluctuation level of the reconstructed signal, perform convolution processing on the reconstructed signal using the Gaussian kernel function, and smooth the denoising result through adaptive Gaussian convolution;

[0048] Step 7: extract the multipath effect component from the residual obtained by decomposition based on the CEEMD algorithm and design a notch filter accordingly, perform multipath effect filtering compensation on the denoising result after adaptive Gaussian convolution smoothing, and eliminate periodic multipath interference.

[0049] Therefore, in the GNSS tunnel entrance slope deformation time series denoising method of an embodiment of the present invention, the advantages of CEEMD and wavelet packet adaptive processing are integrated, and the processing of mixed entropy demarcation, adaptive threshold, dynamic smoothing and multipath compensation are combined to more comprehensively and finely reduce noise while retaining useful signals, especially in the slope monitoring of tunnel entrances in complex environments, which have the characteristics of more noise and non-stationary signals, thereby improving the accuracy and reliability of GNSS deformation monitoring.

[0050] It should be understood that CEEMD is an improved method based on the EEMD (Ensemble Empirical Mode Decomposition) algorithm. It overcomes the problems of mode aliasing and endpoint effects that are common in traditional EMD (Empirical Mode Decomposition) decomposition. The CEEMD decomposition algorithm improves the accuracy and robustness of signal decomposition by introducing white noise perturbations and integrating the results of multiple decompositions.

[0051] In the traditional EMD method, the signal is decomposed into several intrinsic mode functions, each of which represents the local characteristics of the signal. However, a major drawback of the EMD method is that it is easily affected by endpoint effects, which causes artifacts in the decomposed IMFs at the beginning and end of the signal. To alleviate this problem, the EEMD method introduces white noise. By adding noise multiple times and performing EMD decomposition, the multiple IMF results obtained are averaged, thereby reducing the impact of noise on the decomposition. CEEMD further improves EEMD by reducing noise interference and improving the stability of the decomposition through a fully integrated process, ensuring the representativeness and accuracy of each mode function.

[0052] As an optional implementation, in step 2, a Gaussian white noise signal with opposite positive and negative values ​​is added to the original GNSS signal, and the CEEMD algorithm is used to decompose the signal into multiple IMF components, including:

[0053] Step 2-1: Add n pairs of positive and negative Gaussian white noise signals to the original GNSS signal x(t) to generate n pairs of positive and negative Gaussian white noise signals. Each noise pair satisfies And the amplitude standard deviation σ e Set to 10% to 20% of the standard deviation of the original signal, thus obtaining two groups of noisy signals (i.e., noisy signal groups), expressed as:

[0054]

[0055] Where, are i noisy signals with positive noise added; is the i-th noisy signal with negative noise added; and are Gaussian white noises with opposite sign and same amplitude respectively; by adding noise in pairs, the subsequent averaging can cancel out the noise residue;

[0056] Step 2-2: Use EMD algorithm to add Gaussian white noise to the noisy signal and Decompose and obtain two sets of positive and negative IMF components in:

[0057]

[0058] Where, and are the jth IMF components of the i-th addition of positive and negative noise, respectively, j = 1, 2, 3, ..., m, where m represents the number of EMD decomposition layers; and are the positive and negative residual components respectively;

[0059] As an optional approach, it should be understood that in the EMD decomposition algorithm, the number of EMD decomposition layers is adaptively terminated by the signal characteristics (stopped when the residual is a monotonic function);

[0060] Step 2-3: Add the two sets of positive and negative IMF components obtained from the k-th decomposition result and calculate the average value to obtain the j-th IMF component, which is expressed as:

[0061]

[0062] Step 2-4: Output the decomposition results to obtain m IMF components {IMF1,IMF2,...,IMF m} and the residual term r(t), satisfying:

[0063]

[0064] Therefore, CEEMD is used to reduce modal aliasing and preliminarily separate high-frequency noise from low-frequency signals.

[0065] As an optional implementation, in step 3, a mixed entropy of the Pearson correlation coefficient and the fuzzy entropy is constructed, and the local minimum is detected based on the mixed entropy to determine the component IMF corresponding to the first mixed entropy that takes the local minimum. k As the boundary between noise and signal, determine IMF k The components after sequence number k are dominated by noise, and the noisy signals are filtered out, including:

[0066] Step 3-1: For each IMF component j(t), calculate its Pearson correlation coefficient ρ with the original GNSS signal x(t) j ;

[0067] Step 3-2: IMF j (t) is converted into an m-dimensional vector sequence X(i) and the distance d between all vector pairs is calculated. ij And fuzzy membership similarity, based on the fuzzy membership similarity, the fuzzy entropy FE is calculated j ;

[0068] Step 3-3, based on fuzzy entropy FE j And the correlation coefficient ρ j , construct the mixed entropy α j As a joint indicator:

[0069] α j＝ ρ j / (1+FE j );

[0070] Step 3-4: For each IMF component, traverse and calculate the mixing entropy α j Sequence, and use the sliding window method to detect the local minimum, select the first component IMF corresponding to the mixed entropy where the local minimum occurs k As the boundary between noise and signal, determine IMF k The components after sequence number k are dominated by noise, and the noisy signals are filtered out.

[0071] As an alternative, for each IMF component IMF j (t), calculate its Pearson correlation coefficient ρ with the original GNSS signal x(t) j The existing method can be used, and the calculation formula is:

[0072]

[0073] It should be understood that the closer the correlation coefficient ρ is to 1, the j It indicates that the stronger the linear correlation between IMF and the original signal (the more signal components).

[0074] As an optional method, in step 3-2, IMF j (t) is converted into an m-dimensional vector sequence X(i) and the distance d between all vector pairs is calculated. ij And fuzzy membership similarity, based on the fuzzy membership similarity, the fuzzy entropy FE is calculated j , including the following processes:

[0075] Step a. Sequence reconstruction: IMF j(t) Convert to an m-dimensional (where m represents the embedding dimension) vector sequence X(i), {X(i) = [IMF j (i),IMFj(i+1),...,IMF j (i+m-1)]};

[0076] Step b. Distance calculation: Calculate the distance d between all vector pairs ij :

[0077] d ij =max k＝0,...,m-1 |X(i+k)-X(j+k)|;

[0078] Step c. Fuzzy membership calculation: For each distance d ij , calculate the similarity D ij =exp(-(d ij / r) n ), where r represents the similarity tolerance, which is set to r=0.2*std(IMF j );as well as

[0079] Step d. Fuzzy entropy calculation:

[0080] FE j =ln φ m (t)-lnφ m+1 (r);

[0081] in,

[0082] In this embodiment, the embedding dimension is 2. Fuzzy entropy FE j The higher the value, the higher the IMF complexity (noise-dominated), and the lower the value, the stronger the regularity (signal-dominated).

[0083] Further, on this basis, as mentioned above, based on fuzzy entropy FE j And the correlation coefficient ρ j , construct the mixed entropy α j as a joint indicator.

[0084] In this step, a joint index is constructed through normalization to balance correlation and complexity and avoid high amplitude noise causing ρ j Inflated problem.

[0085] Further in step 3-4, the mixing entropy α is calculated by traversing j Sequence, and use the sliding window method to detect the local minimum, select the first component IMF corresponding to the mixed entropy where the local minimum occurs k As the boundary between noise and signal, determine IMF kThe components after sequence number k are dominated by noise, and the noisy signals are filtered out. As an optional implementation method, the specific implementation process includes:

[0086] First, for each IMF component (j=1,2,...,m) obtained by CEEMD decomposition, the mixing entropy α is calculated j , generate the sequence {α1,α2,...,α m};

[0087] Then, use the sliding window method (window length is set to 3) to determine α j Whether it meets:

[0088] α j <α j-1 And α j <α j+1

[0089] If it is satisfied, then it is judged that j corresponds to a local minimum point;

[0090] Therefore, window sliding is used to avoid misjudgment caused by noise fluctuations, avoid interference from high-frequency fluctuations, and ensure boundary stability;

[0091] Then, the IMF number k corresponding to the first local minimum is selected and the IMF is considered k And the subsequent components are dominated by noise; for example, if α3 is the first minimum, then k = 3, and thus:

[0092] IMF3,IMF4,...,IMF m is the noise IMF.

[0093] Furthermore, we can conduct further verification and adjustment, for example: if all α j Monotonically decreasing, the default k = argmin j α j ; If there is no minimum value, k is set to the empirical value m / 2, where m is the number of layers of CEEMD decomposition.

[0094] Therefore, the joint index constructed by the embodiment of the present invention can improve the robustness of noise boundary and avoid misjudgment of traditional methods when the noise amplitude suddenly changes (such as blasting vibration interference).

[0095] As an optional implementation, in step 4, for the filtered noisy signal, noise weights are constructed by performing a sub-band frequency band energy-entropy joint analysis based on wavelet packet decomposition, and a threshold is adaptively determined based on the sub-band noise level and the noise weight to perform piecewise threshold function processing to obtain a sub-band after threshold processing, including:

[0096] Step 4-1: Perform wavelet packet decomposition on the noisy signal layer by layer according to the selected wavelet basis and decomposition level, and further decompose each subband of the previous layer into low-frequency and high-frequency parts layer by layer until the set decomposition level is reached;

[0097] Step 4-2: For each sub-band, calculate its energy proportion E b , and calculate the permutation entropy PE of each permutation pattern based on subband reconstruction b ;

[0098] Step 4-3: Sub-band-based energy ratio E b And the permutation entropy PE of each permutation pattern b , construct the noise weight w b ;

[0099] Step 4-4: Based on the sub-band noise level and noise weight w b , calculate the adaptive threshold T b ;

[0100] Step 4-5: Using the adaptive threshold T b As a benchmark, perform segmented threshold function processing, including:

[0101] 1) For |W l,b |≥T b The coefficients of retain the signal components:

[0102]

[0103] 2) For |W l,b |<T b The coefficient of , exponential decay:

[0104]

[0105] Therefore, through the design of an improved threshold function, the threshold is adaptively selected according to the noise level at each scale, and the advantages of hard threshold and soft threshold are combined to retain signal features and suppress noise, overcoming the limitations of traditional wavelet hard threshold, soft threshold and adaptive threshold methods in complex noise backgrounds, especially in high-frequency deformation feature extraction and signal fidelity. By combining the advantages of hard threshold and soft threshold, efficient denoising is achieved and signal details are retained to the maximum extent, which significantly reduces the generation of artifacts and improves the authenticity and integrity of the signal.

[0106] In the embodiment of the present invention, the wavelet packet decomposition process can use the Db4 wavelet basis (Daubechies 4th order), which is suitable for capturing the local mutation characteristics of GNSS signals by taking advantage of its compact support and moderate smoothness. Set the decomposition level L (for example, L=5) and decompose the noise IMF into 2 L sub-band {Wl,b} (i.e. sub-band), as mentioned above, is decomposed into 2 5 = 32 sub-bands {W l,b}, where l = 1, 2, ..., L, L is the number of wavelet basis decomposition layers, b = 1, 2, ..., 2 l , which is represented by the subband number (i.e. index) of each layer.

[0107] It should be understood that the number of layers of the aforementioned wavelet basis decomposition can be configured to other numbers, and a higher number of layers (such as L=5) can refine the frequency band division. In this embodiment, the computational efficiency is weighed and the value L=5 is taken as an example for illustration.

[0108] In an embodiment of the present invention, the process of the wavelet packet decomposition algorithm is:

[0109] Perform wavelet packet decomposition on the noise IMF layer by layer:

[0110] Layer 1: Decompose the signal into low-frequency approximations (W 1,1 ) and high frequency details (W 1,2 );

[0111] Layer l: Each subband of the previous layer is further decomposed into low-frequency and high-frequency parts until the L layer is reached.

[0112] Therefore, by layer-by-layer wavelet packet decomposition, the high-frequency sub-band (such as W 5,32 ) concentrated random noise, low frequency sub-band (such as W 5,1 ) preserves signal details.

[0113] As an optional implementation, in step 4-2, for each sub-band, calculate its energy proportion E b , and calculate the permutation entropy PE of each permutation pattern based on subband reconstruction b ,include:

[0114] First, based on each sub-band obtained by wavelet packet decomposition, calculate its energy proportion E b :

[0115]

[0116] Where N represents the number of sample points in the subband, that is, the number of sampling points of the signal in each subband; L represents the number of wavelet packet decomposition layers, which determines the number of subbands. The total number of subbands is 2 L ; b represents the index of the target subband, which is used to identify the subband whose energy ratio is currently being calculated; b' represents the index of all subbands to be traversed, b'=1,2,3,…,2 L , used to calculate the total energy; W l，b(t) represents the wavelet packet coefficient of the lth layer and the bth subband. The square of its absolute value represents the energy of the sample point. The total energy of the target subband b is obtained by summing the energies of the N sample points in the subband.

[0117] It should be understood that the energy ratio is low (for example, E b <0.05) are usually noise-dominated;

[0118] Then, each subband W is transformed into l，b (t) Reconstructed into an m-dimensional phase space vector to capture the dynamic structure of the sequence: X(i) = [W l，b (i), W l，b (i+τ),...,W l，b (i+(m-1)τ)]; for each vector X(i), calculate its arrangement mode and count the probability P of each mode (such as ascending, descending, etc.) j , based on which we get the permutation entropy PE b :

[0119]

[0120] Where m! represents the total number of permutations and combinations of m-dimensional phase space vectors; P j Represents the probability of the jth arrangement appearing, satisfying

[0121] Among them, the embedding dimension, that is, the dimension of phase space reconstruction, determines the length of the phase space vector; the delay time is the time interval used to construct the phase space vector, which determines the time span of adjacent elements in the phase space.

[0122] In an embodiment of the present invention, PE b It represents the permutation entropy of sub-band b, reflecting the complexity of the sub-band signal. The negative value comes from the information entropy formula H = -∑P j ln P j , where dividing by ln(m!) is to normalize the result to [0,1] to facilitate comparison of the results under the embedding dimension.

[0123] It should be understood that in the embodiment of the present invention, permutation entropy is used to measure the complexity and regularity of time series. For each phase space vector, the order of the size of its elements is statistically defined as a "permutation pattern", and the probability of all patterns appearing is calculated as P. j The uncertainty is calculated based on this; the more uniform the probability distribution (i.e., the probability of occurrence of multiple permutation patterns is similar), the higher the sequence complexity and the greater the permutation entropy; if only a few patterns dominate, the smaller the entropy value, the stronger the sequence regularity. The formula quantifies this uncertainty by operating on probability. Permutation entropy PEb The closer the value of is to 1, the more random the subband is (noise dominated).

[0124] As an optional implementation, in step 4-3, the sub-band-based energy proportion E b And the permutation entropy PE of each permutation pattern b , construct the noise weight w b ,include:

[0125]

[0126] Where ∈ is the adjustment factor, which is 10 -6 ~10 -5 , used to prevent the denominator from being zero and to adjust the weight ratio of sub-bands with extremely low energy.

[0127] Among them, high PE b (strong randomness) and low E b (weak energy) subband weight w b The larger the value, the more aggressive the denoising is required. The configuration of ∈ is designed to prevent the denominator from being zero, and for very low energy sub-bands (such as E b ≈0) for numerical stabilization.

[0128] As an optional embodiment, in step 4-4, based on the sub-band noise level and the noise weight w b , calculate the adaptive threshold T b ,include:

[0129] First, the sub-band noise level σ is determined using the median absolute error estimate b :

[0130]

[0131] Among them, 0.67450.6745 is the conversion factor between the standard deviation of the Gaussian distribution and the absolute deviation of the median;

[0132] Then, based on the subband noise level σ b and the noise weight w b Calculate the adaptive threshold T b :

[0133] Represents a general threshold benchmark based on Gaussian white noise.

[0134] in, The larger the scale l (low-frequency subband), the smaller the threshold is to preserve details; the noise weight w b The higher it is, the larger the adaptive threshold.

[0135] As an optional implementation, in step 6, the Gaussian kernel width is dynamically determined according to the local fluctuation level of the reconstructed signal, the reconstructed signal is convolved with the Gaussian kernel function, and the denoising result is smoothed by adaptive Gaussian convolution, including:

[0136] Step 6-1. For the reconstructed signal s'(t), calculate the local variance Var(t) using a sliding window [tL, t+L]:

[0137]

[0138] Where t represents the window time length, L represents the sliding step size; the larger the local variance Var(t), the more severe the regional noise fluctuation;

[0139] Step 6-2: Dynamically adjust the Gaussian kernel width σ(t) based on the local variance Var(t), that is, the width of the Gaussian kernel function after dynamic adjustment at time t:

[0140]

[0141] Among them, σ0 represents the initial kernel width, which serves as the benchmark value for dynamic adjustment and determines the basic magnitude of the kernel width; Var(t) represents the local variance corresponding to time t, which characterizes the degree of fluctuation of the data at that time; max(Var) represents the maximum value of the local variance, which is used to normalize Var(t) to constrain the adjustment coefficient

[0142] Step 6-3: Perform Gaussian convolution on the reconstructed signal s'(t) to achieve smoothing and denoising. The convolution process is as follows:

[0143]

[0144] The width of the convolution process is set to 6σ(t)+1, which is used to cover more than 99.5% of the Gaussian distribution energy.

[0145] In combination with the above examples, the local variance Var(t) is used to measure the degree of data fluctuation at time t. When the data fluctuates greatly, such as in areas of sudden noise changes (large fluctuations, large Var(t)), the Gaussian kernel width σ(t) is increased to enhance the smoothing strength, strengthen the smoothing effect of the kernel function, and expand the range of influence of the data points. Conversely, when the data fluctuates less, the kernel width is relatively reduced, and the stable areas are weakly smoothed to preserve data details.

[0146] And, through Normalize the current variance to ensure that the adjustment range is within a reasonable range and avoid excessive changes in the kernel width; use the initial kernel width as a benchmark and dynamically scale the normalized variance to adapt the kernel width to the real-time characteristics of the data.

[0147] Therefore, the denoising results are smoothed by adaptive Gaussian convolution, high-frequency noise is removed and low-frequency structure is retained, further improving the denoising effect. This method shows obvious advantages in noise suppression and signal recovery, and is especially suitable for complex signal processing tasks under high-noise backgrounds. It strikes a balance between noise suppression and signal detail preservation, which not only reduces the impact of residual noise, but also enhances the main features of the signal.

[0148] In an embodiment of the present invention, the above-mentioned improved threshold function is combined with a segmented processing strategy and adaptive Gaussian convolution processing. On the one hand, different denoising strategies are adopted for different signals: when the signal amplitude is greater than the threshold, smooth compression is used to avoid excessive weakening of the signal amplitude, eliminate and suppress strong noise, and retain significant signal features; when the signal amplitude is less than the threshold, exponential decay is used to retain weak signal features, smooth the transition area, avoid pseudo-Gibbs oscillation caused by hard thresholding, and avoid information loss caused by direct zeroing.

[0149] Compared with traditional hard and soft thresholding methods, the improved threshold function proposed in this paper demonstrates superior performance in denoising, particularly in preserving signal detail. While hard thresholding effectively suppresses noise, it can lead to sudden changes at the threshold point, which can easily introduce artifacts. Soft thresholding mitigates this problem through smooth transitions, but can also result in the loss of high-frequency detail.

[0150] In contrast, the improved threshold function can smoothly compress large signals to avoid over-suppression, and smoothly transition through exponential decay in small signal parts to retain their characteristics, thus achieving both noise suppression and maximally retaining signal details.

[0151] Furthermore, while the improved threshold function lacks the dynamic adjustment capabilities of the adaptive threshold function, its processing maintains greater smoothness and continuity at the threshold point, avoiding the discontinuities that can occur in adaptive methods. While the adaptive threshold can dynamically adjust the threshold based on signal characteristics, it has a high computational complexity and is not suitable for applications requiring high real-time performance. In summary, the improved threshold function proposed in this invention achieves a good balance between noise suppression, signal recovery, and computational efficiency, providing an efficient and stable denoising solution.

[0152] After completing the threshold processing, in order to further eliminate the high-frequency noise in the wavelet basis coefficients, the present invention introduces an adaptive Gaussian smoothing mechanism to smooth the denoising result and further reduce the impact of residual noise on signal recovery.

[0153] Furthermore, the noisy signal is decomposed by wavelet packets, and the wavelet packet coefficients after denoising based on improved threshold processing and adaptive Gaussian smoothing are used to reconstruct the GNSS signal.

[0154] In this process, not only the interference of noise is suppressed, but also the important information contained in the original signal is retained, thereby achieving the effect of signal denoising. Through the method proposed in the present invention, the stability and clarity of the signal are enhanced by reducing the high-frequency noise in the signal, which significantly improves the quality and practicality of the signal, laying a solid foundation for subsequent signal analysis and processing.

[0155] As an optional implementation, in step 7, extracting multipath effect components from the residual obtained by decomposition based on the CEEMD algorithm and designing a notch filter based on the components, performing multipath effect filtering compensation on the denoising result after adaptive Gaussian convolution smoothing, and eliminating periodic multipath interference, including:

[0156] Step 7-1: Extract the multipath effect component r from the residual r(t) obtained by decomposing the original GNSS signal using the CEEMD algorithm. mp (t), perform residual signal extraction:

[0157] r mp (t)=r(t)-trend top

[0158] Among them, the aforementioned trend item is set to a layer that can be removed by polynomial fitting;

[0159] Step 7-2: Perform Fourier spectrum analysis on the extracted residual signal to obtain the spectrum R(f), and identify the significant peak spectrum f mp :

[0160]

[0161] Among them, the multipath effect in the typical frequency band of the tunnel entrance environment is concentrated in 0.1-1Hz;

[0162] Step 7-3: Design a notch filter to perform notch filtering on the denoised signal after the Gaussian convolution operation to suppress the energy of the frequency band corresponding to the significant peak frequency, wherein the transfer function of the notch filter is configured as follows:

[0163]

[0164] Among them, f s represents the GNSS sampling frequency; ρ represents the pole radius, which is used to suppress the stopband width Δf:

[0165] Δf=(1-ρ)f mp ;

[0166] Therefore, the denoised signal after the Gaussian convolution operation is passed through a notch filter to suppress f mp The energy in the ±Δf / 2 frequency band is used to eliminate periodic multipath interference and retain the main components of the deformation signal.

[0167] Based on the above specific implementation process, the GNSS tunnel entrance slope deformation time series denoising method proposed in the present invention, which combines CEEMD with the adaptive wavelet packet threshold function, has better signal retention and noise suppression effects. Through simulation experiments and comparisons of four groups of simulated signals, the SNR and RMSE of the denoising method proposed in the present invention are the best.

[0168] In an application to real-world tunnel entrance slope monitoring data, the proposed GNSS tunnel entrance slope deformation time-series denoising method effectively identifies and retains high-frequency deformation information, while also providing a more stable denoising process. Compared to traditional wavelet hard thresholding, soft thresholding, and adaptive thresholding methods, the improved method significantly improves denoising by 46.96%, 42.10%, and 47.14%, respectively. This demonstrates that the proposed method can more accurately remove noise while preserving important deformation features when processing real-world GNSS monitoring data, meeting the denoising requirements of tunnel entrance slope monitoring data.

[0169] The proposed GNSS tunnel entrance slope deformation time-series denoising method has demonstrated promising results in monitoring the impact of blasting during tunnel mining. By monitoring the deformation response caused by tunnel blasting vibration, the improved threshold GNSS tunnel entrance slope deformation time-series denoising method can more accurately identify the displacement effects of blasting on the tunnel entrance slope. In particular, it demonstrates greater sensitivity and stability in detecting abnormal high-frequency displacement changes caused by blasting vibration. This method holds great promise for application in tunnel mining projects, providing effective technical support for blasting vibration monitoring and further improving the reliability and accuracy of engineering monitoring.

[0170] {Example 2}

[0171] In order to verify the denoising effect of the GNSS tunnel entrance slope deformation time series denoising method, four groups of simulation signals were constructed in this embodiment, such as Figure 3 shown.

[0172] (1) Signal 1 (f1)—combination of periodic term and trend term

[0173] The signal f1(t) consists of a periodic term and a trend term. The periodic term contains sine waves and cosine waves of multiple frequencies, and the trend term is a linear growth function. Its mathematical expression is:

[0174] f1(t)=3sin(2π·6·t)+2cos(2π·6·t)+sin·(4π·6·t)+2cos(4π·6·t)+0.005·t·2048

[0175] Where t∈[0,1] is the time vector in seconds. The periodic term of this signal consists of a sine wave and a cosine wave with a frequency of 6 Hz. In addition, it includes a linear growth term to simulate the signal's trend. This signal can simulate scenarios with periodic fluctuations and long-term trends, making it suitable for verifying the performance of denoising algorithms on such signals.

[0176] (2) Signal 2 (f2)—Blocks signal

[0177] The signal f2(t) uses the Blocks signal model, which is a signal with step changes and is often used to simulate sudden or drastic changes in system states. Its mathematical expression is:

[0178]

[0179] Among them, h i is the amplitude of the step function, t i is the starting time of the step function, K(t) is the step function, defined as:

[0180]

[0181] This signal consists of multiple step functions of varying amplitudes and time positions, simulating a sudden change or jump in a signal at a specific moment. The Blocks signal model is often used to simulate interference from nonstationary signals when testing signal recovery and denoising algorithms.

[0182] (3) Signal 3 (f3)—HeaviSine signal

[0183] The signal f3(t) uses the HeaviSine signal model, which is a signal that combines a periodic sine wave and a step function. Its mathematical expression is:

[0184] f3(t)=4sin(4πt)-sgn(t-0.3)-sgn(0.72-t)

[0185] Here, 4sin(4πt) is a sine wave with a frequency of 2 Hz, and sgn(t-0.3) and sgn(0.72-t) are step functions with sudden changes at t = 0.3 and t = 0.72, respectively. The HeaviSine signal combines periodic fluctuations with sudden changes and is widely used to test the responsiveness of denoising algorithms to periodic components and transient changes.

[0186] (4) Signal 4 (f4)—the sum of Signal 2 and Signal 3

[0187] Signal f4(t) is the superposition of signal f2(t) and signal f3(t), and its expression is:

[0188] f4(t)=f2(t)+f3(t)

[0189] Signal 4 combines the characteristics of step changes and periodic fluctuations, simulating the system's response to the superposition of different types of signals. Due to its complex structure, f4(t) can be used to test the performance of denoising algorithms when faced with the overlap of multiple different signal sources.

[0190] In order to simulate the noise interference in the actual signal, in this embodiment, colored noise e is added to the above simulation signal. t The noise is generated by the Gaussian white noise model wgn(t), specifically:

[0191] e t =wgn(t)-0.5·wgn(t-1)

[0192] Where wgn(t) is a Gaussian white noise with a mean of zero and a variance of 10. The noisy signal is obtained by replacing the noise e t Adding this to the original signal gives:

[0193] noisySignal(t)=f(t)+e t

[0194] The noisy signal contains the interference of white noise on the original signal and has typical noise pollution characteristics, which can fully verify the recovery ability of the denoising algorithm under different noise interference conditions.

[0195] Combined with attachment Figure 4 、 5 As shown in Figures 6 and 7, the results of the method proposed in the present invention are compared with the denoising methods using soft threshold, hard threshold and adaptive threshold function. It can be seen that the GNSS tunnel entrance slope deformation time series denoising method proposed in the present invention is significantly higher than the results of the other three traditional algorithms in terms of RMSE and SNR accuracy.

[0196] The GNSS denoising method proposed in the present invention can accurately remove noise while retaining important deformation features, meeting the denoising needs of tunnel entrance slope monitoring data. It not only suppresses noise interference, but also retains the important information contained in the original signal, achieving the effect of signal denoising. At the same time, by reducing the high-frequency noise in the signal, the stability and clarity of the signal are enhanced, significantly improving the quality and practicality of the signal.

[0197] At the same time, the GNSS tunnel entrance slope deformation time series denoising method with improved threshold in the present invention can more accurately identify the displacement impact of blasting on the tunnel entrance slope, especially in the detection of abnormal high-frequency displacement changes caused by blasting vibration. It has higher sensitivity and stability, which makes this method have good application prospects in tunnel mining projects, and can provide effective technical support for blasting vibration monitoring, further improving the reliability and accuracy of engineering monitoring.

[0198] {Example 3}

[0199] According to an embodiment disclosed in the present invention, a computer system is further provided, comprising:

[0200] one or more processors; and

[0201] Memory, which stores instructions that can be operated;

[0202] Among them, when the instruction is executed by one or more processors, it causes the aforementioned one or more processors to perform operations, which include the process of executing the GNSS tunnel entrance slope deformation time series denoising method of the aforementioned embodiment, so as to accurately remove noise while retaining important deformation features, thereby meeting the noise reduction requirements of tunnel entrance slope monitoring data.

[0203] {Example 4}

[0204] According to an embodiment disclosed in the present invention, a computer-readable storage medium is further provided for storing one or more programs, the one or more programs including instructions or instruction sets that can be executed by one or more processors;

[0205] Among them, when the aforementioned instructions or instruction sets are executed by one or more processors, the process of the GNSS tunnel entrance slope deformation time series denoising method of the aforementioned embodiment is executed to accurately remove noise while retaining important deformation features, thereby meeting the noise reduction requirements of tunnel entrance slope monitoring data.

[0206] While the present invention has been disclosed above with reference to preferred embodiments, this is not intended to limit the present invention. Persons skilled in the art will readily appreciate that various modifications and variations can be made without departing from the spirit and scope of the present invention. Therefore, the scope of protection of the present invention shall be determined by the claims.

Claims

1. A GNSS tunnel entrance slope deformation time series denoising method combining CEEMD and adaptive wavelet packet threshold function, characterized in that: The following steps are involved: Step 1: Obtain the original GNSS signal for GNSS tunnel entrance slope deformation monitoring; Step 2: Add a Gaussian white noise signal with opposite positive and negative values ​​to the original GNSS signal and use the CEEMD algorithm to decompose the signal into multiple IMF components. Step 3: Construct the mixed entropy of Pearson correlation coefficient and fuzzy entropy, detect the local minimum based on the mixed entropy, and determine the component IMF corresponding to the first mixed entropy that takes the local minimum. k As the boundary between noise and signal, IMF is determined based on this k The components after sequence number k are dominated by noise, and the noisy signals are filtered out; Step 4: For the filtered noisy signal, noise weights are constructed by performing sub-band frequency band energy-entropy joint analysis based on wavelet packet decomposition, and a threshold is adaptively determined based on the sub-band noise level and the noise weight to perform segmented threshold function processing to obtain sub-bands after threshold processing; Step 5: Perform inverse wavelet packet transform on the sub-band after threshold processing to obtain the denoised noise IMF component as the reconstructed signal; Step 6: Dynamically determine the Gaussian kernel width according to the local fluctuation level of the reconstructed signal, perform convolution processing on the reconstructed signal using the Gaussian kernel function, and smooth the denoising result through adaptive Gaussian convolution; Step 7: extract the multipath effect component from the residual obtained by decomposition based on the CEEMD algorithm and design a notch filter accordingly, perform multipath effect filtering compensation on the denoising result after adaptive Gaussian convolution smoothing, and eliminate periodic multipath interference.

2. The GNSS tunnel entrance slope deformation time series denoising method combining CEEMD and adaptive wavelet packet threshold function according to claim 1 is characterized in that: In step 2, a Gaussian white noise signal with opposite positive and negative values ​​is added to the original GNSS signal, and the signal is decomposed into multiple IMF components using the CEEMD algorithm, including: Step 2-1: Add n pairs of Gaussian white noise signals with opposite positive and negative values ​​to the original GNSS signal x(t) to obtain two sets of noisy signals, expressed as: Where, are i noisy signals with positive noise added; is the i-th noisy signal with negative noise added; and are Gaussian white noises with opposite positive and negative values ​​and the same amplitude, and the amplitude standard deviation σ e Set to 10% to 20% of the standard deviation of the original signal; Step 2-2: Use EMD algorithm to add Gaussian white noise to the noisy signal and Decompose and obtain two sets of positive and negative IMF components in: Where, and are the jth IMF components of the i-th addition of positive and negative noise, respectively, j = 1, 2, 3, ..., m, where m represents the number of EMD decomposition layers; and are the positive and negative residual components respectively; Step 2-3: Add the two sets of positive and negative IMF components obtained from the k-th decomposition result and calculate the average value to obtain the j-th IMF component, which is expressed as: Step 2-4: Output the decomposition results to obtain m IMF components {IMF1,IMF2,...,IMF m } and the residual term r(t), satisfying:

3. The GNSS tunnel entrance slope deformation time series denoising method combining CEEMD and adaptive wavelet packet threshold function according to claim 1 is characterized in that: In step 3, the mixed entropy of the Pearson correlation coefficient and the fuzzy entropy is constructed, and the local minimum is detected based on the mixed entropy to determine the component IMF corresponding to the first mixed entropy with the local minimum value. k As the boundary between noise and signal, determine IMF k The components after sequence number k are dominated by noise, and the noisy signals are filtered out, including: Step 3-1: For each IMF component j (t), calculate its Pearson correlation coefficient ρ with the original GNSS signal x(t) j ; Step 3-2: IMF j (t) is converted into an m-dimensional vector sequence X(i) and the distance d between all vector pairs is calculated. ij And fuzzy membership similarity, based on the fuzzy membership similarity, the fuzzy entropy FE is calculated j ; Step 3-3, based on fuzzy entropy FE j And the correlation coefficient ρ j , construct the mixed entropy α j As a joint indicator: a j＝ r j / (1+FE j ); Step 3-4: For each IMF component, traverse and calculate the mixing entropy α j Sequence, and use the sliding window method to detect the local minimum, select the first component IMF corresponding to the mixed entropy where the local minimum occurs k As the boundary between noise and signal, determine IMF k The components after sequence number k are dominated by noise, and the noisy signals are filtered out.

4. The GNSS tunnel entrance slope deformation time series denoising method combining CEEMD and adaptive wavelet packet threshold function according to claim 1 is characterized in that: In step 4, for the filtered noisy signal, noise weights are constructed by performing sub-band frequency band energy-entropy joint analysis based on wavelet packet decomposition, and a threshold is adaptively determined based on the sub-band noise level and the noise weight to perform segmented threshold function processing to obtain the sub-band after threshold processing, including: Step 4-1: Perform wavelet packet decomposition on the noisy signal layer by layer according to the selected wavelet basis and decomposition level, and further decompose each subband of the previous layer into low-frequency and high-frequency parts layer by layer until the set decomposition level is reached; Step 4-2: For each sub-band, calculate its energy proportion E b , and calculate the permutation entropy PE of each permutation pattern based on subband reconstruction b ; Step 4-3: Sub-band-based energy ratio E b And the permutation entropy PE of each permutation pattern b , construct the noise weight w b ; Step 4-4: Based on the sub-band noise level and noise weight w b , calculate the adaptive threshold T b ; Step 4-5: Using the adaptive threshold T b As a benchmark, perform segmented threshold function processing, including: 1) For |W l,b |≥T b The coefficients of retain the signal components: 2) For |W l,b |<T b The coefficient of , exponential decay:

5. The GNSS tunnel entrance slope deformation time series denoising method combining CEEMD and adaptive wavelet packet threshold function according to claim 4 is characterized in that: In step 4-2, for each sub-band, calculate its energy proportion E b , and calculate the permutation entropy PE of each permutation pattern based on subband reconstruction b ,include: First, based on each sub-band obtained by wavelet packet decomposition, calculate its energy proportion E b : Where N represents the number of sample points in the subband, that is, the number of sampling points of the signal in each subband; L represents the number of wavelet packet decomposition layers, which determines the number of subbands. The total number of subbands is 2 L ; b represents the index of the target subband, which is used to identify the subband whose energy ratio is currently calculated; b ’ Indicates the index of traversing all sub-bands, b ’ =1,2,3,…,2 L , used to calculate the total energy; W l，b (t) represents the wavelet packet coefficient of the lth layer and the bth subband. The square of its absolute value represents the energy of the sample point. The total energy of the target subband b is obtained by summing the energies of the N sample points in the subband. Then, each subband W is embedded according to the configuration parameters of embedding dimension m and delay time τ. l，b (t) Reconstructed into an m-dimensional phase space vector: X(i) = [W l,b (i),W l,b (i+τ),...,W l，b (i+(m-1)τ)]; for each vector X(i), calculate its arrangement pattern and count the probability P of each pattern appearing j , based on which we get the permutation entropy PE b : Where m! represents the total number of permutations and combinations of m-dimensional phase space vectors; P j Represents the probability of the jth arrangement appearing, satisfying 6. The GNSS tunnel entrance slope deformation time series denoising method combining CEEMD and adaptive wavelet packet threshold function according to claim 4 is characterized in that: In step 4-3, the sub-band-based energy proportion E b And the permutation entropy PE of each permutation pattern b , construct the noise weight w b , include: Where ∈ is the adjustment factor, which is 10 -6 ~10 -5 , used to prevent the denominator from being zero and to adjust the weight ratio of sub-bands with extremely low energy.

7. The GNSS tunnel entrance slope deformation time series denoising method combining CEEMD and adaptive wavelet packet threshold function according to claim 4 is characterized in that: In step 4-4, based on the sub-band noise level and the noise weight w b , calculate the adaptive threshold T b ,include: First, the sub-band noise level σ is determined using the median absolute error estimate b : Among them, 0.67450.6745 is the conversion factor between the standard deviation of the Gaussian distribution and the absolute deviation of the median; Then, based on the subband noise level σ b and the noise weight w b Calculate the adaptive threshold T b : in, Represents a general threshold benchmark based on Gaussian white noise.

8. The GNSS tunnel entrance slope deformation time series denoising method combining CEEMD and adaptive wavelet packet threshold function according to any one of claims 1 to 7, characterized in that: In step 6, the Gaussian kernel width is dynamically determined according to the local fluctuation level of the reconstructed signal, the reconstructed signal is convolved using the Gaussian kernel function, and the denoising result is smoothed by adaptive Gaussian convolution, including: Step 6-1. For the reconstructed signal s'(t), calculate the local variance Var(t) using a sliding window [tL, t+L]: Where t represents the window time length and L represents the sliding step size; Step 6-2: Dynamically adjust the Gaussian kernel width σ(t) based on the local variance Var(t), that is, the width of the Gaussian kernel function after dynamic adjustment at time t: Among them, σ0 represents the initial kernel width, which serves as the benchmark value for dynamic adjustment and determines the basic magnitude of the kernel width; Var(t) represents the local variance corresponding to time t, which characterizes the degree of fluctuation of the data at that time; max(Var) represents the maximum value of the local variance, which is used to normalize Var(t) to constrain the adjustment coefficient Step 6-3: Perform Gaussian convolution on the reconstructed signal s'(t) to achieve smoothing and denoising. The convolution process is as follows: The width of the convolution process is set to 6σ(t)+1, which is used to cover more than 99.5% of the Gaussian distribution energy.

9. The GNSS tunnel entrance slope deformation time series denoising method combining CEEMD and adaptive wavelet packet threshold function according to any one of claims 1 to 7, characterized in that: In step 7, the multipath effect component is extracted based on the residual obtained by decomposition using the CEEMD algorithm, and a notch filter is designed accordingly. The multipath effect filtering compensation is performed on the denoising result after adaptive Gaussian convolution smoothing to eliminate periodic multipath interference, including: Step 7-1: Extract the multipath effect component r from the residual r(t) obtained by decomposing the original GNSS signal using the CEEMD algorithm. mp (t), perform residual signal extraction: r mp (t) = r(t) - trend term Among them, the aforementioned trend item is set to a layer that can be removed by polynomial fitting; Step 7-2: Perform Fourier spectrum analysis on the extracted residual signal to obtain the spectrum R(f), and identify the significant peak spectrum f mp : Among them, the multipath effect in the typical frequency band of the tunnel entrance environment is concentrated in 0.1-1Hz; Step 7-3: Design a notch filter to perform notch filtering on the denoised signal after the Gaussian convolution operation to suppress the energy of the frequency band corresponding to the significant peak frequency, wherein the transfer function of the notch filter is configured as follows: Among them, f s represents the GNSS sampling frequency; ρ represents the pole radius, which is used to suppress the stopband width Δf: Δf=(1-ρ)f mp ; Therefore, the denoised signal after the Gaussian convolution operation is passed through a notch filter to suppress f mp The energy in the ±Δf / 2 frequency band is used to eliminate periodic multipath interference and retain the main components of the deformation signal.

10. A computer system, characterized in that: include: one or more processors; as well as Memory, which stores instructions that can be operated; Wherein, when the instruction is executed by one or more processors, the aforementioned one or more processors perform an operation, and the operation includes the process of executing the GNSS tunnel entrance slope deformation time series denoising method combining CEEMD and adaptive wavelet packet threshold function as described in any one of claims 1-9.

11. A computer-readable storage medium for storing one or more programs, characterized in that: The one or more programs include instructions or sets of instructions that can be executed by one or more processors; Wherein, when the instruction or instruction set is executed by one or more processors, the process of the GNSS tunnel entrance slope deformation time series denoising method combining CEEMD and adaptive wavelet packet threshold function described in any one of claims 1 to 9 is executed.

Images