GPS Sparse Multipath Signal Estimation Based on Compressive Sensing

Abstract

A GPS sparse multipath signal estimation method based on compressive sensing is proposed. A new 0 norm approximation function is designed, and the parameter of the approximate function is gradually reduced to realize the approximation of 0 norm. The sparse signal is reconstructed by a modified Newton method. The reconstruction performance of the proposed algorithm is better than several commonly reconstruction algorithms at different sparse numbers and noise intensities. The GPS sparse multipath signal model is established, and the sparse multipath signal is estimated by the proposed reconstruction algorithm in this paper. Compared with several commonly used estimation methods, the estimation error of the proposed method is lower.

This work is licensed under http://creativecommons.org/licenses/by/4.0/ (the “License”). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.

1. Introduction

In wireless communication system, due to the influence of scattering or reflection, radio wave may reach the receiver through multiple different paths, so that the received signal at the receiver is the superposition of multiple signals. The multipath interference is inevitable, which is one of the main reasons that affect the performance of GNSS receiver. It will lead to the distortion of the autocorrelation waveform of the pseudocode and affect the detection of the coherent peak, thus reducing the positioning accuracy of the receiver. It is necessary to detect the multipath signals [1, 2].

Candès et al. [3] and Donoho [4] proposed the theory of compressive sensing (CS) in 2006, which can sample the sparse signal at low sampling rate. Currently, the research and application fields of CS include radar signal arrival angle estimation [5, 6], satellite navigation signal processing [7–9], and sparse radar design [10–13]. The CS signal reconstruction needs to recover the original high-dimensional signal from a small amount of measurement data. There are several kinds of reconstruction algorithms. The relaxation algorithm transforms 0 norm into 1 norm. It converts the problem into a convex optimization problem which is easy to be solved [14, 15]. The greedy algorithm has the advantages of simple calculation and easy implementation. However, the reconstruction accuracy is not high and the sparse number is needed [16, 17]. In recent years, some 0 norm approximate solution algorithms have been developed, which transform the 0 norm into mathematical analyzable approximation functions. The sparse signal is reconstructed by a gradient method, Newton method, or other optimization methods [18, 19]. This method can reduce the difficulty of signal reconstruction and obtain better results.

A new 0 norm approximation function is proposed. The sparse reconstruction problem of CS is translated into a Lagrange multiplier problem. It gradually reduces the parameters and reconstructs the original signal by the modified Newton method. A GPS sparse multipath signal model is established, and the estimation of GPS sparse multipath signal is realized by using the approximate function which is proposed in this paper. The rest of this paper is organized as follows. In Section 2, the GPS sparse multipath model is introduced. Section 3 analyzes the proposed 0 norm approximate function and reconstruction method. Sections 4 and 5 give the reconstruction and estimation results of the new proposed method. Finally, the paper is concluded in Section 6.

2. GPS Sparse Multipath Model

Assuming that the GPS receiver has accurately tracked the carrier frequency of the receiver signal and stripped the carrier, the received multipath signal in one navigation message can be expressed as$$\begin{array}{}\text{(1)}& rt=\sum _{k=1}^L{A}_{k}st-{\tau }_k{e}^{j{\varphi }_k}+nt.\end{array}$$

In the formula, $A_k$, ${\tau }_k$, and ${\varphi }_k$ are the amplitude, delay time, and phase of the multipath signal, respectively. $st$ is the spread spectrum code, and $nt$ is the noise [20]. The corresponding signals of all sampling points can be obtained by cyclic shift, so the upper formula can be represented as a matrix form:$$\begin{array}{}\text{(2)}& \mathbf{r}=\mathbf{S}\mathbf{a}+\mathbf{n},\end{array}$$where $\mathbf{S}={\mathbf{s}}_{{\tau }_1},{\mathbf{s}}_{{\tau }_2},\cdots ,{\mathbf{s}}_{{\tau }_N}$ and $a_k=A_k{e}^{j{\varphi }_k}$. When there are multipath signals, the vector values equal to the amplitude value of the multipath signal, and the other values are 0. The vector $a$ can be represented as$$\begin{array}{}\text{(3)}& ai=\begin{array}{cc}{a}_k,\hfill & \text{if}{\tau }_i={\tau }_k,\hfill \\ 0,\hfill & \text{elsewhere}.\hfill \end{array}\end{array}$$

The length of a GPS spread spectrum chip is 1/1023 ms. Suppose the sampling frequency of the GPS signal is 51.15 MHz. So there are 50 sampling points in a chip. Assuming that the number of multipath in a chip is 3; that is to say, three elements are not 0 and the others are 0. $\mathbf{r}$ is a sparse vector. The GPS multipath signal represented by formula (3) is a sparse model.

There are many GPS multipath estimation methods. In reference [21], a TK operator method is easy to implement but the estimation effect is not ideal in strong noise environment. The method based on maximum likelihood (ML) estimation is proposed in reference [22]. The ML method is improved, and the Generalized Likelihood Ratio (GLRT) is used to improve the performance of the estimation [23]. The MEDLL is used to analyze the signal, and the effect of MEDLL is good, but the method is complex [24]. There are other wavelet-based analysis methods and so on [25, 26]. These methods estimate GPS multipath signals without considering the sparse characteristics of GPS multipath. If sparse constraints are added to the estimation, the performance of multipath estimation can be improved theoretically.

3. New 0 Norm Approximation Function

3.1. Design of Approximate 0 Norm

Donoho et al. proposed a CS-based sparse reconstruction method for magnetic resonance imaging (MRI) images [27]. Instead of 0 norm, 1 norm is used to reconstruct the MRI images. The 1 norm is changed to the following formula:$$\begin{array}{}\text{(4)}& {\mathbf{x}}_1\approx \sum _{i=1}^N{x_i^2+p}^{1/2}.\end{array}$$

In the formula, $p$ is a small positive constant. The approximate formula is used to transform the nondifferential 1 norm into a differentiable function. Then, the MRI image is reconstructed by using the optimization algorithm. A CS reconstruction algorithm based on 0 norm approximation function is proposed [18]. The approximation function is a Gaussian function:$$\begin{array}{}\text{(5)}& {\mathbf{x}}_0=\underset{\sigma \to 0}{\lim }\sum _{i=0}^N\exp -\frac{x_i^2}{2{\sigma }^2}.\end{array}$$

The function realizes the approximate expression of 0 norm by decreasing the parameter $\sigma$ gradually. Then, the sparse signal is reconstructed by the gradient method. Both methods convert a nonconvex 0 norm problem to a convex function which has some deviation in principle. Furthermore, the parameter $p$ of formula (4) remains unchanged, and the reconstruction may have some deviation.

Based on the approximation of 0 norm given by Donoho and Gaussian function, a new approximate function is proposed. The coefficient $1/2$ is changed to $p/2$ by modifying the approximation of 1 norm given by Donoho. An approximate 0 norm can be obtained when the coefficient $p$ decreases gradually to 0. The designed 0 norm approximation function is$$\begin{array}{}\text{(6)}& {\mathbf{x}}_0=\underset{p\to 0}{\lim }\sum _{i=0}^N{x_i^2+p^2}^{p/2}.\end{array}$$

The model has only one parameter, and the approximate function is a smooth differentiable function. So the gradient and Hessian matrix of the function can be obtained, which facilitates the solution of the model.

3.2. Sparse Signal Reconstruction Based on Approximate 0 Norm

Based on the proposed approximation function, the CS reconstruction problem can be described as the following mathematical constraint expression:$$\begin{array}{c}\text{(7)}& \underset{x}{\min }\underset{p\to 0}{\lim }\sum _{i=0}^N{x_i^2+p^2}^{p/2}\text{s}.\text{t}.\hspace{1em}\mathbf{y}=\mathbf{\Phi }\mathbf{x}.\end{array}$$

$\mathbf{\Phi }\in R^{m\times n}$ is a measurement matrix. In reference [18], the gradient method is used to reconstruct the sparse signal. The descending direction is “zigzag,” which affects the efficiency of reconstruction. In this paper, the Newton method is used to reconstruct the sparse signal for improving the efficiency [28, 29]. Newton’s method requires positive definite Hessian matrix, so the Hessian matrix needs to be modified in iteration. The CS reconstruction problem is transformed into a Lagrange solution model:$$\begin{array}{}\text{(8)}& \underset{x}{\min }f=\frac{1}{2}{\mathbf{\Phi }\mathbf{x}-\mathbf{y}}^2+\lambda \sum _{i=0}^N{x_i^2+p^2}^{p/2}.\end{array}$$

$y$ is Lagrange multiplier. The gradient of formula (8) is$$\begin{array}{}\text{(9)}& \mathbf{g}={\mathbf{\Phi }}^{\text{T}}\mathbf{\Phi }\mathbf{x}-\mathbf{y}+\lambda \begin{array}{c}p{x_1^2+p^2}^{p/2-1}{x}_1\\ \begin{array}{c}.\\ .\end{array}\\ p{x_N^2+p^2}^{p/2-1}{x}_N\end{array}.\end{array}$$

The Hessian matrix of formula (8) is obtained by derivation of expression (9):$$\begin{array}{}\text{(10)}& \mathbf{H}={\mathbf{\Phi }}^{\text{T}}\mathbf{\Phi }+\lambda \ast \begin{array}{c}{u}_{1}0\cdots 0\\ \begin{array}{cccc}0& u_2& \cdots & 0\end{array}\\ \begin{array}{cccc}.& .& .& .\end{array}\\ \begin{array}{cccc}0& 0& \cdots & u_N\end{array}\end{array}.\end{array}$$

$u_i=p{x_i^2+p^2}^{p/2-2}p-1x_i^2+p^2$. When $p=1$, $u_i$ is a positive value, so the Hessian matrix is a positive definite matrix. That is to say, the function is a convex function. The value $p=1$ can be used as the starting point of iteration. According to Newton’s method, the iterative formula can be expressed as$$\begin{array}{}\text{(11)}& {\mathbf{x}}_{k+1}={\mathbf{x}}_k+{\alpha }_k{\mathbf{d}}_k,\text{(12)}& {\mathbf{d}}_k=-{\mathbf{H}}^{-1}{\mathbf{g}}_{k-1}.\end{array}$$

In the formula ((11) and (12)), ${\alpha }_k$ is the scale factor and $d_k$ is the direction of gradient descent. $g_{k-1}$ is the gradient of the $k-1$ iteration which can be obtained by formula (9). In the iteration, $u_i$ may be a negative value in formula (10). Formula (12) requires a positive definite Hessian matrix. It needs to modify the value of$u_i$ to guarantee a positive definite matrix in the iteration. The modification formula is$$\begin{array}{}\text{(13)}& u_i=\begin{array}{c}\begin{array}{cc}{u}_i,\hfill & \text{if}\hspace{1em}{u}_i>\delta >0,\hfill \end{array}\hfill \\ \begin{array}{cc}\delta ,\hfill & \text{else}.\hfill \end{array}\hfill \end{array}\end{array}$$

$\delta$ is a very small positive value which ensures the Hessian matrix is a positive definite matrix. The empirical value is ${10}^{‐5}$. The initial value of the parameter $p$ in the approximate function is 1. The value is decreased in iteration according to the following formula:$$\begin{array}{}\text{(14)}& \begin{array}{cc}{p}_k={cp}_{k-1},& k=1,\cdots ,J\end{array}.\end{array}$$

It can obtain better reconstruction effect when $c$ belong to $0.5,1$.

The general idea of sparse reconstruction algorithm in this paper is given as below. The reconstruction has two search processes. For each given value $p_k$ in the external loop, an optimal value is obtained by the modified Newton method. An upgrading value is gotten by decreasing $p$ gradually. Then, $p_k$ in the outer loop is reduced to the value $p_{k+1}$, and the new optimal value is searched again according to the previous optimal solution. The sparse solution is found iteratively. The implementation steps of this reconstruction Algorithm 1 are as follows:

Algorithm 1: GPS sparse multipath signal estimation based on compressive sensing.

Input: random measurement matrix $\mathbf{\Phi }$, measurement value $\mathit{y}$;

Output: reconstruction signal $\mathit{x}\mathit{R}$;

(1) Initialization: Set the initial parameter: $p_0$, $\lambda$, $c$, ${\alpha }_k$, Total iteration times $J$, inner iteration $L_k$, reconstruction error threshold $\epsilon$, the signal initial value ${\mathit{x}}_0=\mathbf{0}$, external iteration count $k$, internal iteration count $t$;

(2) The external iteration: if $k>J$, turn to step 5. Otherwise, $x_k$ is calculated by Newton method;

(3) The internal iteration: $t=t+1$, if $t>L_k$, the internal iteration is over and turn to the next step. Otherwise, calculate ${\mathit{x}}_k^t$ by ${\mathit{x}}_k$ and $p_k^t$ which decrease gradually. If $E_t^k=\mathrm{norm}{\mathbf{H}}_t^k<\epsilon$, turn to the next step;

(4) $p_{k+1}=p_k$, $k=k+1$, and turn to step 2, searching the optimal value based on ${\mathit{x}}_k^t$;

(5) All the iteration is over and get the sparse solution $\mathit{x}\mathit{R}={\mathit{x}}_k^t$;

Output2: There is no satellite in the line of sight.

3.3. Convergence Analysis of Approximate Function

In this section, the convergence of the algorithm is analyzed. In order to facilitate the analysis, a lemma is first introduced.

Lemma 1.

Suppose that $\mathbf{A}\in R^{m\times n}$ is a Unique Representation Property (URP) matrix. All $m\times m$ submatrixes have inverse matrix [30]. A vector $\mathbf{s}={s_1,s_2,\cdots ,s_n}^T$ belongs to zero space of $\mathbf{A}$. If more than $n-m$ elements converge to 0 in the vector, then the vector converges to 0.

Proof.

Suppose that $\mathbf{A}$ has been normalized by a column vector. For any $\gamma >0$, $I_{\gamma }$ is a subscript set that all the corresponding $s_i$ are larger than $\gamma$. $\widehat{\mathbf{A}}$ is a submatrix which is composed of the columns corresponding to the subscript set.

Because $\mathbf{A}$ is a URP matrix, each column satisfies the linear independence. The matrix $\widehat{\mathbf{A}}$ has a left inverse matrix. $M=\max {\stackrel{\wedge }{\mathbf{A}}}^{-1}$; $I_{\gamma }$ is defined as the number of elements in the subscript set. From the known conditions,$$\begin{array}{}\text{(15)}& \mathbf{A}\mathbf{s}=\sum _{i=1}^m{s}_i{a}_i=\sum _{i\in I_{\gamma }}{s}_i{a}_i+\sum _{i\notin I_{\gamma }}{s}_i{a}_i.\end{array}$$

So$$\begin{array}{}\text{(16)}& \sum _{i\in I_{\gamma }}{s}_i{a}_i=-\sum _{i\in I_{\gamma }}{s}_i{a}_i.\end{array}$$

There are 2 norms for the upper formulas:$$\begin{array}{}\text{(17)}& X=x_1,x_2,\cdots ,x_L\in R^{N\times L}=\sum _{i\notin I_{\gamma }}\gamma =n-I_{\gamma }\gamma \le n\gamma .\end{array}$$

$\widehat{\mathbf{A}}=a_i,i\in I_{\gamma }$; rewrite the upper formula to a matrix form:$$\begin{array}{}\text{(18)}& \sum _{i\in I_{\gamma }}{s}_i{a}_i={\mathbf{A}\mathbf{s}}_{i\in I_{\gamma }}\Rightarrow {{\mathbf{s}}_{i\in I_{\gamma }}}_2={{\stackrel{\wedge }{\mathbf{A}}}^{-1}\sum _{i\in I_{\gamma }}{s}_i{a}_i}_2\le {{\stackrel{\wedge }{\mathbf{A}}}^{-1}}_2{\sum _{i\in I_{\gamma }}{s}_i{a}_i}_2\le {{\stackrel{\wedge }{\mathbf{A}}}^{-1}}_{2}n\gamma .\end{array}$$

And ${\mathbf{s}}_2\le {s_{i\in I_{\gamma }}}_2+{s_{i\notin I_{\gamma }}}_2\le {{\stackrel{\wedge }{\mathbf{A}}}^{-1}}_2+1n\gamma$, so$$\begin{array}{}\text{(19)}& {\mathbf{s}}_2\le {{\stackrel{\wedge }{\mathbf{A}}}^{-1}}_2+1n\gamma \le M+1n\gamma .\end{array}$$

Proof of lemma 1 is over.

Corollary 2.

Suppose that $\mathbf{A}\in R^{m\times n}$ is a URP matrix and the equation $\mathbf{y}=\mathbf{A}\mathbf{s}$ has a sparse solution. If there is another solution $\tilde{\mathbf{s}}={{\tilde{\mathbf{s}}}_1,{\tilde{\mathbf{s}}}_2,\cdots ,{\tilde{\mathbf{s}}}_n}^T$ and ${\tilde{\mathbf{s}}}_0=k<m/2$, the vector converges to 0.

The proof of corollary can be drawn from the lemma. Because $0\le {\tilde{\mathbf{s}}-{\mathbf{s}}^0}_0\le 2k<m$, more than $n-m$ elements are 0, and $\tilde{\mathbf{s}}-{\mathbf{s}}^0$ belongs to the zero space of $\mathbf{A}$. The matrix is a URP matrix, so the difference vector $\tilde{\mathbf{s}}-{\mathbf{s}}^0$ converges to 0. That is to say, the sparse solution is unique.

For the approximate function proposed in this paper, suppose that $\zeta$ is a very small positive number. If $\gamma =\pm {1-\zeta }^{1/p}$, for a given $\gamma \le s_i\le 1$, the following equations hold:$$\begin{array}{}\text{(20)}& \underset{p\to 0}{\lim }{f}_p{s}_i=\underset{p\to 0}{\lim }{s_i^2+p^2}^{p/2}\ge {{1-\zeta }^{2/p}}^{p/2}\approx 1.\end{array}$$

Because the range of $p$ is $0,1$, so ${1-\zeta }^{1/p}<1-\zeta <1$. For a very small coefficient $s_i$ in the vector, the value of the proposed function is approximately 1.

For any given $\gamma =\pm {1+\zeta }^{1/p}$, the following equations hold:$$\begin{array}{}\text{(21)}& \underset{p\to 0}{\lim }{f}_p{s}_i=\underset{p\to 0}{\lim }{s_i^2+p^2}^{p/2}=\underset{p\to 0}{\lim }{{1+\zeta }^{2/p}+p^2}^{p/2}=\underset{p\to 0}{\lim }1+\zeta {1+\frac{p^2}{{1+\zeta }^{2/p}}}^{p/2}\approx 1.\end{array}$$

For all coefficients which are greater than 1, the function values are approximately 1. Therefore, the proposed approximate function in this paper realizes the function of 0 norm.

The above analysis theoretically ensures that the approximate function can find the sparse solution of the signal. The constraint condition $k<m/2$ is a sufficient rather than a necessary condition, and the constraint condition is strong. If this condition is not satisfied in practical application, the algorithm can reconstruct the sparse signal also, but the probability of reconstruction is very low.

4. Simulation Experiment of Reconstruction Algorithm

Three experiments are designed to verify the feasibility and effectiveness of the algorithm. In experiment 1, the measurement matrix is fixed. The reconstruction performance is analyzed at different sparse numbers. Antinoise performance is analyzed in experiment 2. The running time of the algorithm is also an important standard to test the algorithm. Experiment 3 analyzes the running time of algorithm.

4.1. Different Sparse Numbers

According to the proposed approximate function, a CS reconstruction model is defined as$$\begin{array}{c}\text{(22)}& \underset{s}{\min }\underset{p\to 0}{\lim }\sum _{i=0}^N{s_i^2+p^2}^{p/2}\text{s}.\text{t}.\hspace{1em}{\mathbf{y}-\mathbf{\Phi }\mathbf{s}}_2<\epsilon .\end{array}$$

In order to compare the reconstruction effect, the proposed algorithm is compared with OMP, Donoho method (1 norm) [27], and the reference [18] algorithm (SL0). The noise type is zero mean Gaussian white noise, and the noise variance is 0.01 in experiment 1. The measurement matrix is $\Phi \in R^{200\times 500}$. The length of measurement value and sparse signal is fixed at 200 and 500, respectively. The simulation is run 1000 times at each sparse number, and the total reconstruction times are counted if the signal-to-noise ratio (SNR) of the reconstructed signal is greater than 25 dB. The experimental results are shown in Table 1.

Table 1

Successful times at different sparse numbers.

Sparse number	50	60	70	80	90	100	110	120	130
OMP	1000	998	981	944	729	228	25	0	0
1 norm	1000	1000	995	978	864	403	96	0	0
SL0	1000	1000	1000	996	947	429	185	3	0
This paper	1000	1000	1000	1000	981	506	227	7	0

As can be seen from Table 1, each algorithm can reconstruct the signal completely when the sparse number is less than 50. When the sparse number is greater than 130, the reconstruction SNR of all algorithms are less than 25 dB. When the sparse number increases, the successful reconstruction times of the proposed algorithm in this paper are greater than other algorithms. For example, when the number of nonzero data reaches 120, there are still 7 times reconstructed signals with SNR above 25 dB, but SL0 algorithm only has 3 times and 1 norm and the OMP algorithm can no longer reconstruct the signal above 25 dB. Comparing the data in Table 1, the reconstruction effect algorithm proposed in this paper is better.

4.2. Robust Analysis of the Proposed Algorithm

The robust analysis of the proposed algorithm is analyzed at different noise intensities in experiment 2. The noise variance is [0.01, 0.02, 0.03, 0.04]. The other parameters are consistent with experiment 1. The average value of the reconstructed SNR is obtained from 1000 experiments. The experimental results are shown in Table 2.

Table 2

Robust analysis of the proposed algorithm.

Noise variance	Sparse number
	40	50	60	70	80	90	100	110
0.01	33.108	32.689	32.179	31.574	31.129	29.956	25.272	14.594
0.02	26.785	26.777	26.510	26.280	26.139	25.003	19.559	11.744
0.03	20.602	20.693	20.678	20.652	20.411	19.488	15.945	9.897
0.04	16.952	16.920	16.827	16.747	16.570	14.910	11.887	7.680

From Table 2, it can be seen that the reconstruction SNR is above 25 dB at low sparse number when the noise variance is 0.01 and 0.02. The proposed algorithm is effective. The reconstruction SNR can reach 25.272 dB when the sparse number is 100 and noise variance is 0.01. The proposed algorithm shows strong robustness. However, when the noise variance is 0.04, the reconstruction effect of this algorithm is not ideal.

4.3. Average CPU Time of Algorithm

The reconstruction time is an important index to estimate the algorithm. The experimental parameter is designed as follows: the measurement matrix $\mathbf{\Phi }\in R^{m\times n}$, $m/n=0.4$, and sparse number $m=4k$. According to the previous experiments, this numerical setting can effectively reconstruct sparse signals in a noiseless environment. The length of the signal is constantly adjusted, and the CPU times for each of the algorithms averaged over 100 trials are tabulated in Table 3.

Table 3

Reconstruction CPU time of different algorithms (second).

Algorithm	Signal length
	500	1000	2000	3000	4000
OMP	0.160	0.262	1.816	9.095	21.696
1 norm	0.245	0.258	1.365	6.220	15.461
SL0	0.261	0.322	1.313	5.998	14.669
This paper	0.204	0.270	1.299	5.881	11.577

From Table 3, it can be found that the running time of several algorithms is close when the total length of the signal is less than 2000 sampling points. However, when the signal length reaches 4000, the reconstruction CPU time of OMP is almost twice that of the proposed algorithm. For high-dimensional matrix, the pseudoinverse calculation of the matrix in the OMP algorithm takes too much time. SL0 is solved by the optimal gradient method which appears “zigzag” when finding the optimal value. This reduces the efficiency of finding the optimal solution. In this paper, the modified Newton method is used to improve the reconstruction speed. According to Table 3, the algorithm is faster than the other algorithms, so the algorithm is more suitable for high-dimensional signal processing.

5. GPS Sparse Multipath Signal Estimation

The GPS sparse multipath model is given in formula (3). The sparse multipath signal is estimated by CS theory in this section. First, the GPS signal is measured by measurement matrix $\mathbf{\Phi }\in R^{m\times n}$:$$\begin{array}{}\text{(23)}& z=\mathbf{\Phi }r=\mathbf{\Phi }\mathbf{S}\mathbf{a}+\mathbf{\Phi }\mathbf{n}.\end{array}$$

Then, sparse multipath signal is estimated by using the 0 norm approximation function which is proposed in this paper:$$\begin{array}{}\text{(24)}& \underset{a}{\min }fa=\frac{1}{2}{\mathbf{\Phi }\mathbf{S}\mathbf{a}-\mathit{z}}^2+\lambda \sum _{i=0}^N{a_i^2+p^2}^{p/2}.\end{array}$$

The multipath model discussed in this paper is based on the following assumptions [31, 32]:

(1) The transmission delay of the nondirect signal is greater than that of the direct signal, which is in line with the actual model. The nondirect signal is often formed by refraction or reflection which takes more time

(2) The amplitude of nondirect signal is lower than that of the direct signal. In the shelter or indoor environment, this assumption is generally not true which needs to estimate multipath with an additional method. But this situation is beyond the scope of this paper which is not analyzed here

(3) In this paper, the multipath delay in on chip is only considered. That is to say, the multipath delay is limited to one chip

(4) Multiple sparse values may be estimated due to noise interference. Assuming that the number of sparse multipath is known, the maximum sparse value is selected as the multipath signal. The small sparse values is considered as noise and excluded

5.1. Comparison of Different Sparse Multipath Estimation Methods

The TK operator [21] and ML estimate method [22] are selected to estimate the multipath signal. Assuming that the GPS signal has two nondirect signals. The direct signal delay is 0, and the normalized amplitude is 1. The delay time of the nondirect signal changes randomly in a chip. The ranges of absolute amplitude vary randomly between 0.2 and 0.8. The GPS signal length is 1 ms and sampling rate is 51.15 MHz. The SNR of correlation result is 35 dB. The compression rate is $m/n=0.4$. The estimation results of different methods are shown in Figure 1. The amplitude of two random nondirect signals is (0.6, 0.4), and the position is (0.2, 0.58) chip.

[figure omitted; refer to PDF]

From Figure 1, three methods can correctly estimate the multipath signal delay with some deviation in magnitude. The TK operator and ML estimate method can correctly find multipath signals. Due to noise interference, the estimation amplitude deviation of multipath signals is large, and the estimated values of nonsparse positions fluctuate greatly which result in more fake multipath signals. The CS estimation method proposed in this paper adds sparse prior information, the estimated multipath signal is consistent with the simulation signal position, and the estimated multipath signal amplitudes are close to the real multipath signal.

ML and TK operators are not effective because of noise interference. After the multipath estimation of the signal, some detection algorithm is used to detect it, and the false multipath signal is excluded in practical application. For example, the GLRT criterion is used to analyze the estimated multipath signal in reference [23]. The proposed CS method adds sparse constraints and selects the largest signal as a sparse multipath signal. It can be seen from Figure 1 that the performance is better than that of ML and TK operators. It is unfair to compare the CS estimation method with the common methods, so it only compares the method of this paper with the method based on the CS theory in the subsequent performance analysis.

5.2. Multipath Signal Estimation Based on CS

It is mainly affected by three factors: noise, sparse number (number of multipath), and measurement matrix dimension $m$ to estimating multipath signal based on CS. This section analyzes the estimated effect from these aspects. First, the sparse estimation results based on CS are analyzed by an experiment. 1 ms GPS signal is selected with sampling rate 40.92 MHz. The sparse number is set to 4, and the delay of nondirect signal varies randomly in a chip. The absolute amplitude varies randomly between 0.2 and 0.8. The delay and amplitude are uniformly distributed. SNR is set to 0 dB. Compression ratio $m/n$ is 0.5. In the simulation, assuming that the number of sparse multipath is known, the estimated results of the proposed method are shown in Figure 2. In a strong noise environment, the amplitude of the multipath signal and the position of the multipath signal may be deviated. Multipath signal estimation error is used to evaluate the estimation results, which is defined as$$\begin{array}{}\text{(25)}& \text{Error}=\frac{1}{K}\sum _{i=1}^N{a}_i-{\widehat{a}}_i.\end{array}$$

[figure omitted; refer to PDF]

In the formula, $N$ is the total number of sampling point of a chip. $K$ is a sparse number. $a_i$ and ${\widehat{a}}_i$ represent real and estimated multipath signals, respectively. All multipath signal amplitudes are normalized. It can be seen from Figure 2 that if the position and phase of the estimated signal are consistent with the real multipath signal, the estimation error is small. If the position or phase of estimated signal is different from the real signal, the real value and the estimated value are counted into error. Formula (25) can objectively evaluate the estimated performance.

The multipath estimation effect is analyzed at different multipath numbers and noise intensities in experiment 1. The measurement matrix is fixed and $m/n=0.3$. The SNR and sparse number vary from -20 dB to 20 dB and 2-8, respectively. The other parameters are set as above. The experiment repeats 100 times at each sparse number. The average estimation error is shown in Figure 3. From the figure, it can be found that the estimation error increases when the sparse number increase. It is consistent with the CS theory that the reconstruction effect decreases with the increase of sparsity at the same measurement matrix. The best estimation effect is obtained when there is only one nondirect multipath signal. The reconstruction becomes unstable when $K=8$.

[figure omitted; refer to PDF]

When the multipath number is fixed, the influence of measurement matrix is analyzed in experiment 2. The multipath number is set to 6. The column of the matrix remains unchanged, and $m/n$ range from 0.2 to 0.8. The other related parameters are set in the same experiment 1. The experiment repeats 100 times at each $m/n$. The average estimation error is shown in Figure 4. The estimation error will be relatively reduced when the number of random measurement is increased, which can retain more information of sparse signal. The estimation effect is not good when the SNR is less than 0 dB.

[figure omitted; refer to PDF]

The reconstruction algorithm has some effect on sparse signal reconstruction in CS, which also has some influence on sparse multipath estimation. The OMP, 1 norm, and SL0 are selected to estimate the multipath signal. The estimation effect of different algorithms is analyzed in experiment 3. The sparse number of signal is set to 5. The $m/n$ of the measurement matrix is 0.5. The GPS signal length is 1 ms with 40.92 MHz sampling rate. The experiment repeats 100 times at each SNR, and the average estimation error is shown in Figure 5. The estimation error of all methods is large at low SNR. When the SNR is higher than 0 dB, the estimation effect is good. As the SNR increases, three multipath estimation errors based on optimization algorithms are slightly lower than OMP methods.

[figure omitted; refer to PDF]

From the three experiments, it can be seen that the sparse multipath signal can be estimated by the CS method effectively. When the SNR is low, the estimation effect is poor. For example, the CS method cannot estimate the multipath signal at -20 dB strong noise interference. As the SNR increases, the estimation error decreases. When the SNR is higher than 5 dB, the estimated multipath signal is basically consistent with the added multipath signal. In this paper, the length of GPS signal is 1 ms. In practical application, the SNR can be improved by increasing the GPS signal length because the time of navigation message is 20 ms, which can obtain lower estimation error than 1 ms.

6. Conclusions

In this paper, we proposed a new GPS multipath signal estimation method based on CS. A new 0 norm approximate function is proposed. The sparse signal can be reconstructed by using the modified Newton method. A large number of simulation results show that the new method can reconstruct the sparse signal effectively at different sparse numbers and SNR. The reconstruction performance of the proposed algorithm is better than OMP and other optimization algorithms. The GPS signal sparse multipath model is established. The GPS multipath signal can be estimated by the proposed method in this paper. The estimation result is better than the TK and ML method and other CS reconstruction methods. Because there are 20 periods of spread spectrum code in a navigation message, the SNR can be improved by accumulating for a long time, and then, the estimation effect can be improved. Our research team will continue to delve into how to improve the estimation effect of this method.

Acknowledgments

This research was partly funded by Special Program of National Major Scientific Instruments Equipments Development of China, grant number 2013YQ200607; Major Projects of Natural Science Foundation of Anhui Province universities, grant number KJ2019ZD35; and Ministry of Education Key Lab. of Cognitive Radio and Information Processing, grant number CRKL200202.

References

[1] E. D. Kaplan, Understanding GPS: Principles and Application, 2005.

[2] L. Cheng, K. Wang, M. Ren, G. Yan, "Adaptive Filter Approach for GPS Multipath Estimation Under Correntropy Criterion in Dynamic Multipath Environment," IEEE Transactions on Signal Processing, vol. 67 no. 22, pp. 5798-5810, DOI: 10.1109/TSP.2019.2946028, 2019.

[3] E. J. Candès, J. Romberg, T. Tao, "Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information," IEEE Transactions on Information Theory, vol. 52 no. 2, pp. 489-509, DOI: 10.1109/TIT.2005.862083, 2006.

[4] D. L. Donoho, "Compressed sensing," IEEE Transactions on Information Theory, vol. 52 no. 4, pp. 1289-1306, DOI: 10.1109/TIT.2006.871582, 2006.

[5] D. Meng, X. Wang, M. Huang, L. Wan, B. Zhang, "Robust weighted subspace fitting for DOA estimation via block sparse recovery," IEEE Communications Letters, vol. 24 no. 3, pp. 563-567, DOI: 10.1109/LCOMM.2019.2958913, 2020.

[6] X. Wang, M. Huang, L. Wan, "Joint 2D-DOD and 2D-DOA estimation for coprime EMVS–MIMO radar," Circuits Systems and Signal Processing,DOI: 10.1007/s00034-020-01605-5, 2021.

[7] G. He, M. Song, X. He, Y. Hu, "GPS signal acquisition based on compressive sensing and modified greedy acquisition algorithm," IEEE Access, vol. 7, pp. 40445-40453, DOI: 10.1109/ACCESS.2019.2906682, 2019.

[8] G. He, M. Song, S. Zhang, P. Song, X. Shu, "Sparse GLONASS signal acquisition based on compressive sensing and multiple measurement vectors," Mathematical Problems in Engineering, vol. 2020,DOI: 10.1155/2020/9654120, 2020.

[9] A. Osman, M. M. E. Moussa, M. Tamazin, M. J. Korenberg, A. Noureldin, "DOA elevation and azimuth angles estimation of GPS jamming signals using fast orthogonal search," IEEE Transactions on Aerospace and Electronic Systems, vol. 56 no. 5, pp. 3812-3821, DOI: 10.1109/TAES.2020.2988424, 2020.

[10] F. Wen, J. Shi, "Fast direction finding for bistatic EMVS-MIMO radar without pairing," Signal Processing, vol. 173, article 107512,DOI: 10.1016/j.sigpro.2020.107512, 2020.

[11] Q. Cheng, A. A. Ihalage, Y. Liu, Y. Hao, A. Noureldin, "Compressive Sensing Radar Imaging With Convolutional Neural Networks," IEEE Access, vol. 8 no. 5, pp. 212917-212926, DOI: 10.1109/access.2020.3040498, 2020.

[12] J. Shi, F. Wen, T. Liu, "Nested MIMO radar: coarrays, tensor modeling and angle estimation," IEEE Transactions on Aerospace and Electronic Systems, vol. 57 no. 1, pp. 573-585, DOI: 10.1109/TAES.2020.3034012, 2021.

[13] J. He, L. Li, T. Shu, "Sparse nested arrays with spatially spread square acoustic vector sensors for high accuracy underdetermined direction finding," IEEE Transactions on Aerospace and Electronic Systems,DOI: 10.1109/TAES.2021.3057682, 2021.

[14] X. Wang, L. Wang, M. Li, G. Bi, "Nuclear norm minimization framework for DOA estimation in MIMO radar," Signal Processing, vol. 153, pp. 147-152, DOI: 10.1016/j.sigpro.2016.12.031, 2017.

[15] E. Michael, Sparse and Redundant Representations: From Theory to Applications in Signal and Image Processing, 2010.

[16] J. Tropp, A. C. Gilbert, "Signal recovery from random measurements via orthogonal matching pursuit," IEEE Transactions on Information Theory, vol. 53 no. 12, pp. 4655-4666, DOI: 10.1109/TIT.2007.909108, 2007.

[17] J. A. Tropp, A. C. Gilbert, M. J. Strauss, "Algorithms for simultaneous sparse approximation. Part I: greedy pursuit," Signal Processing, vol. 86 no. 3, pp. 572-588, DOI: 10.1016/j.sigpro.2005.05.030, 2006.

[18] H. Mohimani, M. Babie-Zadeh, C. Jutten, "A fast approach for overcomplete sparse decomposition based on smoothed ℓ 0 norm," IEEE Transactions on Signal Processing, vol. 57 no. 1, pp. 289-301, DOI: 10.1109/TSP.2008.2007606, 2009.

[19] F.-Y. Wu, K. Yang, R. Duan, "Compressed sensing of underwater acoustic signals via structured approximation $l_0$-norm," IEEE Transactions on Vehicular Technology, vol. 67 no. 9, pp. 8504-8513, DOI: 10.1109/TVT.2018.2850305, 2018.

[20] B. Kai, D. M. Akos, N. Bertelsen, A Software-Defined GPS and Galileo Receiver, 2007.

[21] R. Hamila, E. S. Lohan, M. Renfors, "Subchip multipath delay estimation for downlink WCDMA system based on Teager-Kaiser operator," IEEE Communications Letters, vol. 7 no. 1,DOI: 10.1109/LCOMM.2002.807439, 2003.

[22] E. Ertin, U. Mitra, S. Siwamogsatham, "Maximum-likelihood-based multipath channel estimation for code-division multiple-access systems," IEEE Transactions on Communications, vol. 49 no. 2, pp. 290-302, DOI: 10.1109/26.905887, 2001.

[23] N. Sokhandan, J. T. Curran, A. Broumandan, G. Lachapelle, "An advanced GNSS code multipath detection and estimation algorithm," GPS Solutions, vol. 20 no. 4, pp. 627-640, DOI: 10.1007/s10291-015-0475-z, 2016.

[24] B. R. Townsend, P. C. Fenton, K. J. Van Dierendonck, D. J. R. Van Nee, "Performance evaluation of the multipath estimating delay lock loop," Navigation, vol. 42 no. 3, pp. 502-514, DOI: 10.1002/j.2161-4296.1995.tb01903.x, 1995.

[25] L. Cheng, K. Wang, M. Ren, G. Yan, "Adaptive filter approach for GPS multipath estimation under correntropy criterion in dynamic multipath environment," IEEE Transactions on Signal Processing, vol. 67 no. 22, pp. 5798-5810, DOI: 10.1109/TSP.2019.2946028, 2019.

[26] P. Zhong, X. L. Ding, D. W. Zheng, W. Chen, D. F. Huang, "Adaptive wavelet transform based on cross-validation method and its application to GPS multipath mitigation," GPS Solutions, vol. 12 no. 2, pp. 109-117, DOI: 10.1007/s10291-007-0071-y, 2008.

[27] M. Lustig, D. Donoho, J. M. Pauly, M. R. I. Sparse, "Sparse MRI: the application of compressed sensing for rapid MR imaging," Magnetic Resonance in Medicine, vol. 58 no. 6, pp. 1182-1195, DOI: 10.1002/mrm.21391, 2007.

[28] L. Wan, X. Kong, F. Xia, "Joint range-Doppler-angle estimation for intelligent tracking of moving aerial targets," IEEE Internet of Things Journal, vol. 5 no. 3, pp. 1625-1636, DOI: 10.1109/JIOT.2017.2787785, 2018.

[29] A. Antoniou, W. S. Lu, Practical Optimization: Algorithms and Engineering Applications, 2010.

[30] I. F. Gorodnitsky, B. D. Rao, "Sparse signal reconstruction from limited data using FOCUSS: a re-weighted minimum norm algorithm," IEEE Transactions on Signal Processing, vol. 45 no. 3, pp. 600-616, DOI: 10.1109/78.558475, 1997.

[31] P. Misra, P. Enge, Global Positioning System: Signals, Measurements, and Performance, 2011.

[32] Y. Hu, M. Z. Song, X. Y. Dang, H. L. Yan, "Interference mitigation for the GPS receiver utilizing the cyclic spectral analysis and RR-MSWF algorithm," Radioengineering, vol. 26 no. 3, pp. 798-807, DOI: 10.13164/re.2017.0798, 2017.