A GNSS multipath signal suppression method based on TVF-EMD-SVD

Abstract

A GNSS multipath signal suppression method based on TVF-EMD-SVD, firstly, a GNSS signal containing multipath effect is subjected to time-varying filtering empirical mode decomposition to obtain m intrinsic mode function components, and multipath signals are often concentrated in low-frequency intrinsic mode function components. The m intrinsic mode function components obtained by decomposition are used to construct a time-frequency matrix, and the time-frequency matrix constructed by IMF components is subjected to singular value decomposition to decompose the frequency components of the multipath signal into different time-frequency subspaces, and then the correlation coefficient between the autocorrelation function of each time-frequency subspace and the autocorrelation function of the original signal is calculated, and the time-frequency subspace with a correlation coefficient greater than or equal to 0.5 is selected for signal reconstruction, thereby achieving suppression of multipath signals. The present invention can avoid the mode aliasing that occurs when EMD is decomposed, and at the same time, for the problem that the reconstruction order cannot be accurately selected when SVD is used for noise reduction, the autocorrelation function and the autocorrelation coefficient are introduced to help SVD to accurately select the reconstruction order.

Description

Technical Field

[0001] The present invention belongs to the field of satellite navigation and relates to a GNSS multipath signal suppression method, in particular to a GNSS multipath signal suppression method based on TVF-EMD-SVD. Background Art

[0002] With the construction and application of the Global Navigation Satellite System (GNSS), satellite navigation has ushered in a new era of human navigation, and people have increasingly higher requirements for its positioning accuracy. However, satellite signals will be reflected or refracted by obstacles during propagation, resulting in the receiver antenna receiving not only an electromagnetic wave signal transmitted from the satellite and propagating in a straight line, but also one or more signals reflected by the electromagnetic wave from surrounding objects. Each reflected signal may reach the antenna after one or more reflections, thus causing multipath errors. Therefore, the multipath effect has become one of the main error sources that restrict positioning accuracy.

[0003] At present, there are three main methods for suppressing GNSS multipath signals: one is based on the design of the receiver antenna; the second is based on the improvement of the receiver tracking loop; and the third is multipath suppression based on data post-processing. The first two have certain difficulties and disadvantages in implementation: based on the design of the antenna end, such as the use of chokes, photonic crystal antennas, etc., although multipath signals can be suppressed from the source of the receiver, they are large in size and high in cost; multipath suppression technology based on the receiver tracking loop, such as multipath estimation delay lock loop (Multipath Estimation Delay Lock Loop, MEDLL), Strobe correlator, PAC technology, etc., experiments with multipath signal suppression by increasing the number of correlators, but it increases the amount of calculation inside the receiver, making the cost of the receiver higher. In recent years, with the continuous development of signal processing technology, various filtering and noise reduction algorithms based on data processing are becoming the main methods for multipath suppression. Empirical Mode Decomposition (EMD) is a data-driven adaptive time domain decomposition method with strong local adaptive decomposition characteristics. It is suitable for processing nonlinear and non-stationary signals. However, this method will cause changes in amplitude during the decomposition process, and there are usually errors in reconstruction. In view of the low-frequency characteristics of GNSS multipath signals, the intrinsic mode function (IMF) obtained after traditional EMD decomposition has a small amount of modal aliasing in the low-frequency band, resulting in the reconstructed signal still containing a small amount of multipath signals. Summary of the invention

[0004] In order to solve the above technical problems, the present invention proposes a GNSS multipath signal suppression method based on TVF-EMD-SVD.

[0005] A GNSS multipath signal suppression method based on TVF-EMD-SVD, the specific steps are as follows:

[0006] Step 1: Select the GNSS signal with multipath effect as the input signal s(t);

[0007] Step 2: Calculate the local cutoff frequency;

[0008] First, the instantaneous amplitude A(t) and instantaneous frequency of the signal s(t) are obtained by Hilbert transform Determine the local maximum of A(t) max}) and the local minimum A({t min}),as well as The local maximum With local minimum Then perform the difference operation to get η 1 (t) and η 2 (t), and at the same time, the input signal is filtered using a time-varying filter to obtain the instantaneous mean a 1 (t) and instantaneous envelope a 2 (t), and then calculate the local cutoff frequency

[0009] Step 3: Reconstruct the signal into h(t) according to the adjusted cutoff frequency, let

[0010]

[0011] Step 4: Determine whether the stopping criteria are met and calculate the cutoff criteria θ(t)

[0012]

[0013] in is the weighted mean instantaneous frequency; B L (t) is the Loughlin instantaneous bandwidth; is the instantaneous frequency of the i-th order component;

[0014] Determine whether the stopping criterion is met. If θ(t)≤ξ, where ξ is a given bandwidth threshold, then s(t) is determined to be an IMF. Otherwise, the extreme points of h(t) are taken as nodes, h(t) is divided into n segments, each with a step length of m. B-spline interpolation is used to approximate the signal s(t), and the approximation result is m(t), that is, let s 1(t) = s(t) - m(t), repeat the above steps (2) to (4); after performing TVF-EMD decomposition on the signal, a set of IMF components imf is obtained i (t), that is Among them, IMF refers to the intrinsic mode function, i is the order of the modal component, and imf i (t) is the i-th order modal component sequence, N is the total number of IMFs obtained by decomposition;

[0015] Step 5: Use the obtained IMF to construct the time-frequency matrix P;

[0016] Step 6: Perform singular value decomposition on the time-frequency matrix;

[0017] Step 7: Calculate the correlation function and correlation coefficient between the time-frequency subspace and the original signal;

[0018] Step 8: Select the time-frequency subspace with a correlation coefficient greater than or equal to 0.5 for reconstruction, and perform inverse singular value decomposition transform on the selected subspace;

[0019] Step 9: Convert the reconstructed signal into a one-dimensional signal, that is, the obtained signal is a signal after multipath suppression.

[0020] As a further improvement of the present invention, the multipath signal s(t) in step 1 should satisfy:

[0021] Where A is the signal amplitude, p(t) is the XOR sum of the data code and the pseudo code with a value of ±1, f is the carrier frequency after considering the Doppler effect, and α i is the attenuation coefficient of the reflected wave, τ i is the propagation delay of the reflected wave i relative to the direct wave, It is the phase change of the signal before and after reflection on each reflection surface.

[0022] As a further improvement of the present invention, the bandwidth threshold ξ in step 4 is 0.08, and the B-spline order is 20.

[0023] As a further improvement of the present invention, the specific operation steps of the singular value decomposition in step 6 are as follows:

[0024] (1) Perform singular value decomposition on the time-frequency matrix P and obtain the decomposition form shown in the following formula

[0025] P=UQV T

[0026] Among them, U is an m×m left orthogonal matrix, V is an n×n right orthogonal matrix, and Q is a diagonal matrix, which can be specifically expressed as follows:

[0027] Ω=diag(σ1 ,σ 2 ,σ 3 ,…,σ q )

[0028] In the above formula, σ i (i=1,2,3,…,q) are the singular values ​​of the time-frequency matrix P, and σ 1 >σ 2 >σ 3 >...>σ q > 0, q is the total number of singular values ​​obtained by decomposition

[0029] (2): According to the correlation coefficient calculated by the time-frequency subspace, the time-frequency subspace with a correlation coefficient greater than or equal to 0.5 is selected, and the number of valid singular values ​​is determined to be k.

[0030] (3): Keep the first k valid singular values ​​of the diagonal matrix Ω, and then set the remaining smaller singular values ​​to 0 to obtain a new singular matrix Q′. Substitute the new singular matrix Q′, the left orthogonal matrix U m×m and the right orthogonal matrix V n×n Substitute it in and calculate to get the new matrix P′.

[0031] Compared with the prior art, the present invention adopts the above technical solution and has the following technical effects:

[0032] 1. TVF-EMD is fully adaptive and suitable for the analysis of linear and non-stationary signals. Compared with EMD, the proposed method is able to improve the frequency separation performance, as well as the stability at low sampling rates. The local cutoff frequency is adaptively designed by making full use of the instantaneous amplitude and frequency information. Then, a non-uniform B-spline approximation is used as a time-varying filter. In order to solve the intermittent problem, a cutoff frequency rearrangement algorithm is also introduced. In order to improve the performance at low sampling rates, the bandwidth criterion of the intrinsic mode function is proposed. In addition, the proposed method is robust to noise interference.

[0033] 2. In the traditional SVD denoising algorithm, some of the retained singular values ​​are mainly noise contributions, and the number of singular values ​​cannot be accurately selected, so the pure signal matrix estimate obtained contains more noise. In this method, the time-frequency characteristic information of multipath signals is decomposed into a series of singular vector time-frequency subspaces composed of U and V, and the time domain information is expressed by V, and the frequency domain information is expressed by U. Therefore, the effective time domain components of the left and right singular quantity signals are selected to reconstruct the order of the matrix. Compared with the traditional SVD denoising, this method can retain more useful signals. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] Figure 1 Algorithm flow chart of the present invention;

[0035] Figure 2 Multipath signal waveform diagram;

[0036] Figure 3 Decomposition results of multipath signal after TVF-EMD decomposition;

[0037] Figure 4 Time-frequency subspace autocorrelation function;

[0038] Figure 5 The signal after suppression. DETAILED DESCRIPTION

[0039] The present invention is further described in detail below in conjunction with the accompanying drawings and specific embodiments:

[0040] like Figure 1 As shown, the present invention discloses a GNSS multipath signal suppression method based on TVF-EMD-SVD. First, a GNSS signal containing multipath effect is decomposed by TVF-EMD to obtain m intrinsic mode function components. The multipath signal is often concentrated in the low-frequency intrinsic mode function components. The m intrinsic mode function components obtained by decomposition are used to construct a time-frequency matrix. The time-frequency matrix constructed by the IMF component is subjected to singular value decomposition to decompose the frequency information of the multipath signal into different time-frequency subspaces. Then, the correlation coefficient between the autocorrelation function of each time-frequency subspace and the autocorrelation function of the original signal is calculated. According to the correlation coefficient, a suitable time-frequency subspace is selected to reconstruct the signal, thereby realizing the suppression of the multipath signal.

[0041] Step 1: Input the GNSS signal s(t) containing low-frequency multipath signals.

[0042]

[0043] Where A is the signal amplitude, p(t) is the XOR sum of the data code and the pseudo code with a value of ±1, f is the carrier frequency after considering the Doppler effect, and α i is the attenuation coefficient of the reflected wave, τ i is the propagation delay of the reflected wave i relative to the direct wave, It is the phase change of the signal before and after reflection on each reflection surface.

[0044] Step 2: Use Hilbert transform to calculate the instantaneous amplitude A(t) and instantaneous frequency of the multipath signal s(t) As shown in formula (2) and (3):

[0045]

[0046] In the formula, is the Hilbert transform of the signal.

[0047] Step 3: Determine the local maximum of A(t) max}) and the local minimum A({t min}), and determine The local maximum With local minimum

[0048] Step 4: Interpolate the extreme points of A(t) to obtain curves β 1 (t) and β 2 (t).

[0049] Step 5: Calculate the instantaneous mean a 1 (t) and the instantaneous envelope a 2 (t), where:

[0050] a 1 (t) = [β 1 (t)+β 2 (t)] / 2(4)

[0051] a 2 (t) = [β 1 (t)-β 2 (t)] / 2 (5)

[0052] Step 6: Right and Interpolate respectively and get η 1 (t) and η 2 (t), and calculated according to equations (6) and (7) and

[0053]

[0054] Step 7: Calculate the local cutoff frequency

[0055]

[0056] Step 8: Reconstruct the signal according to the adjusted cutoff frequency, as shown in equation (9):

[0057]

[0058] The extreme points of h(t) are taken as nodes, and h(t) is divided into n segments, each with a step length of m, where n is called the order of the B-spline.

[0059] Formula (10) is used for B-spline interpolation approximation, and the approximation result is recorded as m(t), which represents the mean of the curve.

[0060]

[0061] Where: [·] is the downsampling operation, is a prefilter and

[0062]

[0063] β n (t) is the B-spline function; n is the B-spline order, m is the node, t is the time, and * represents the convolution operation.

[0064] Step 9: Calculate the cutoff frequency θ(t) and the given bandwidth threshold ξ. If θ(t)≤ξ, the signal can be considered as a narrowband signal, that is, s(t) is determined to be an IMF. Otherwise, let s 1 (t) = s(t) - m(t), repeat steps 2 to 9. Weighted average instantaneous frequency With B L Instantaneous bandwidth B L The calculation of is as shown in formula (13) to (15);

[0065]

[0066]

[0067] After the above steps, the original signal can be expressed as follows:

[0068] s(t)=∑s i (t) (16)

[0069] In the formula, s i (t) is the decomposed subsequence.

[0070] Assume that the number of sampling points of the one-dimensional multipath signal s(t) is n, and the signal is decomposed by TVF-EMD to obtain m IMF components. The data points of each IMF component are a ik (i=1…m;k=1…n), then the time-frequency matrix P can be expressed as:

[0071]

[0072] Step 10: Perform singular value decomposition on the time-frequency matrix P to obtain the decomposition form shown in the following equation (18):

[0073] P=UQV T (18)

[0074] Among them, U is an m×m left orthogonal matrix, V is an n×n right orthogonal matrix, and Q is a diagonal matrix, which can be specifically expressed as follows:

[0075]

[0076] In the above formula, σ i (i=1,2,3,…,q) are the singular values ​​of the time-frequency matrix P, and σ 1 >σ 2 >σ 3 >...>σ q > 0, q is the total number of singular values ​​obtained by decomposition.

[0077] Since the matrix P contains the time-frequency characteristic information of the signal, the series of sub-matrices P obtained by singular value decomposition of the matrix P are i It also contains time-frequency feature information.

[0078] Step 11: Calculate the autocorrelation functions R of the time-frequency subspace s1 ,R s2 ,...,R sj and the autocorrelation function R of the original signal x , the autocorrelation function is shown in formula (20):

[0079]

[0080] In the formula, s(i) is the state of the signal at a certain moment, and N is the number of sampling points of the signal. Normalize the autocorrelation function and find R s1 ,R s2 ,...,R sj The correlation coefficient ρ(j) is calculated as shown in formula (21):

[0081]

[0082] Where j is the jth component.

[0083] The larger the absolute value of the correlation coefficient is, the stronger the correlation is; the closer the correlation coefficient is to 1 or -1, the stronger the correlation is; the closer the correlation coefficient is to 0, the weaker the correlation is. The threshold value of the correlation coefficient in this method is 0.5.

[0084] Step 12: According to the correlation coefficient calculated in the time-frequency subspace, the number of valid singular values ​​is determined to be k.

[0085] Step 13: retain the first k valid singular values ​​of the diagonal matrix Ω, and then set the remaining smaller singular values ​​to 0 to obtain a new singular matrix Q′. Substitute the new singular matrix Q′, the left orthogonal matrix U and the right orthogonal matrix V into the above formula (19) to obtain the new matrix P′.

[0086] Step 14: Convert the matrix A′ into a one-dimensional signal s′(t), which is the signal after noise reduction.

[0087] Example 1

[0088] The present invention provides a GNSS multipath signal suppression method based on TVF-EMD-SVD, which uses simulated GNSS multipath signals for suppression. The relevant parameter settings and steps are as follows:

[0089] The sampling rate of the multipath signal is 1000 Hz, and the sampling time is 2 s. In this embodiment, the amplitude of the direct signal is 1, the frequency is 50 Hz, and there is no phase delay; the amplitude of the first multipath signal is 0.8, the frequency is 3.5 Hz, and the phase delay is The amplitude of the second multipath signal is 0.62, the frequency is 15 Hz, and the phase delay is At the same time, we also add a mean of 0 and a variance of σ 2 =1 Gaussian white noise, the waveform of the multipath signal is as follows Figure 2 shown.

[0090] After steps 2-9, the above signal is decomposed into 9 IMF components, where the bandwidth threshold ξ=0.08 and the B-spline order is 20. Figure 3 shown.

[0091] The autocorrelation function is calculated in step 10 as follows: Figure 4 The autocorrelation coefficients are shown in Table 1, where the components with autocorrelation coefficients greater than 0.5 are S7 and S9.

[0092] Table 1 Correlation coefficients between S1-S9 and the original signal

[0093]

[0094] After steps 11-14, the signals before and after inhibition are as follows Figure 5 As shown. By calculating the signal-to-noise ratio and root mean square error before and after suppression, it is found that the signal-to-noise ratio before suppression is -0.4939dB and the root mean square error is 0.7466; the signal-to-noise ratio after suppression is 1.7610 and the root mean square error is 0.5775. After comparison, it is found that this method can effectively suppress multipath signals and improve the signal-to-noise ratio of signals.

[0095] It will be understood by those skilled in the art that, unless otherwise defined, all terms (including technical and scientific terms) used herein have the same meaning as those generally understood by those skilled in the art in the art to which the present invention belongs. It should also be understood that terms such as those defined in common dictionaries should be understood to have meanings consistent with the meanings in the context of the prior art, and will not be interpreted with idealized or overly formal meanings unless defined as herein.

[0096] The above description is only a preferred embodiment of the present invention and does not constitute any other form of limitation to the present invention. Any modification or equivalent change made based on the technical essence of the present invention still falls within the scope of protection required by the present invention.

Claims

1. A GNSS multipath signal suppression method based on TVF-EMD-SVD, characterized by: The specific steps are as follows: Step 1: Select the GNSS signal containing multipath effect as the input signal s(t); Step 2: Calculate the local cutoff frequency; First, the instantaneous amplitude A(t) and instantaneous frequency of the signal s(t) are obtained by Hilbert transform Determine the local maximum of A(t) max }) and the local minimum A({t min }),as well as The local maximum With local minimum Then, the difference operation is performed to obtain η1(t) and η2(t). At the same time, the input signal is filtered using a time-varying filter to obtain the instantaneous mean a1(t) and instantaneous envelope a2(t), and then the local cutoff frequency is calculated. Step 3: Reconstruct the signal into h(t) according to the adjusted cutoff frequency, let Step 4: Determine whether the stopping criteria are met and calculate the cutoff criteria θ(t) in is the weighted mean instantaneous frequency; B L (t) is the Loughlin instantaneous bandwidth; is the instantaneous frequency of the i-th order component; Determine whether the stopping criterion is met. If θ(t)≤ξ, where ξ is a given bandwidth threshold, then s(t) is determined to be an IMF. Otherwise, the extreme point of h(t) is used as a node, h(t) is divided into n segments, each with a step length of m, and the signal s(t) is approximated using B-spline interpolation. The approximation result is m(t), that is, s1(t)=s(t)-m(t), and the above steps 2 to 4 are repeated; After performing TVF-EMD decomposition on the signal, a set of IMF components imf is obtained i (t), that is Among them, IMF refers to the intrinsic mode function, i is the order of the modal component, and imf i (t) is the i-th order modal component sequence, N is the total number of IMFs obtained by decomposition; Step 5: Use the obtained IMF to construct the time-frequency matrix P; Step 6: Perform singular value decomposition on the time-frequency matrix; Step 7: Calculate the correlation function and correlation coefficient between the time-frequency subspace and the original signal; Step 8: Select the time-frequency subspace with a correlation coefficient greater than or equal to 0.5 for reconstruction, and perform inverse singular value decomposition transform on the selected subspace; Step 9: Convert the reconstructed signal into a one-dimensional signal, that is, the obtained signal is a signal after multipath suppression.

2. The GNSS multipath signal suppression method based on TVF-EMD-SVD according to claim 1, characterized in that: The multipath signal s(t) in step 1 should satisfy: Where A is the signal amplitude, p(t) is the XOR sum of the data code and the pseudo code with a value of ±1, f is the carrier frequency after considering the Doppler effect, and α i is the attenuation coefficient of the reflected wave, τ i is the propagation delay of the reflected wave i relative to the direct wave, It is the phase change of the signal before and after reflection on each reflection surface.

3. The GNSS multipath signal suppression method based on TVF-EMD-SVD according to claim 1, characterized in that: The bandwidth threshold ξ in step 4 is 0.08, and the B-spline order is 20.

4. The GNSS multipath signal suppression method based on TVF-EMD-SVD according to claim 1, characterized in that: The specific operation steps of the singular value decomposition in step 6 are as follows: (1) Perform singular value decomposition on the time-frequency matrix P and obtain the decomposition form shown in the following formula P=UQV T Among them, U is an m×m left orthogonal matrix, V is an n×n right orthogonal matrix, and Q is a diagonal matrix, which can be specifically expressed as follows: In the above formula, σ i (i=1,2,3,…,q) are the singular values ​​of the time-frequency matrix P, and σ1>σ2>σ3>…>σ q > 0, q is the total number of singular values ​​obtained by decomposition; (2): According to the correlation coefficient calculated by the time-frequency subspace, the time-frequency subspace with a correlation coefficient greater than or equal to 0.5 is selected, and the number of valid singular values ​​is determined to be k; (3): Keep the first k valid singular values ​​of the diagonal matrix Ω, and then set the remaining smaller singular values ​​to 0 to obtain a new singular matrix Q′. Substitute the new singular matrix Q′, the left orthogonal matrix U m×m and the right orthogonal matrix V n×n Substitute it in and calculate to get the new matrix P′.

Images