Bayesian EM approach for GNSS parameters of interest estimation under constant modulus interference

originally


Introduction
Global navigation satellite systems (GNSSs) [1] have a wide range of applications, extending beyond navigation and timing to fields like Earth observation, attitude estimation, and space weather characterization.As a result, the accuracy of position, navigation, and timing information is crucial, especially for critical applications like intelligent transportation systems and autonomous unmanned ground/air vehicles.While GNSS has become the primary source of positioning, it was originally designed for optimal performance in clear-sky conditions, making its reliability susceptible to degradation in challenging environments.For instance, phenomenon such as multipath (reflections) [2], spoofing and interferences (intentional or unintentional) are the most challenging ones, being a key issue in safety-critical scenarios [3], such as civilian aviation [4].These effects have been reported in the state of the art, and several mitigation countermeasures have already been proposed [5].In the field of intentional interference, real-world scenarios have identified jammers broadcasting interference characterized by a constant modulus (CM).Initially, these devices Page 2 of 24 Lesouple and Ortega EURASIP Journal on Advances in Signal Processing (2024) 2024:32 utilized constant amplitude tones to disrupt receiver functionality.With this type of signal (even with their straightforward structure), they were able to prevent the receiver from functioning.As a countermeasure, in the time domain, two common methods are employed • pulse blanking [6], which involves zeroing-out samples of the incoming signal exceeding a predefined power threshold to mitigate the impact of pulsed interference • adaptive notch filtering [7], where the jamming signal's instantaneous frequency is continuously estimated using a recurrence equation in the time domain, and the corresponding frequency components are filtered out from the incoming signal.This approach avoids the need for frequency-domain transformations.
However, they proved ineffective when the tone's frequency varied, paving the way for chirp interference, which remains a significant issue today.Notably, notch filters fail to attenuate chirp interference adequately, particularly with nonlinear frequency variations.Even with linear variations, these countermeasures encounter challenges.
Alternative approaches involve signal processing techniques such as the discrete Fourier transform (DFT), which project the signal into the frequency domain, allowing the application of a threshold to eliminate suspicious elements.Another transformation with potential for interference mitigation is the Karhunen-Loeve transform (KLT) [8,9], which relies on the eigenvalues and eigenvectors of the incoming signal's autocorrelation.However, these methods often overly degrade the signal, especially with wide chirp bandwidths.
In this article, we introduce a novel approach to mitigate interference characterized by a CM.This interference category is among the most prominent forms of intentional interference reported in the literature and includes signals like pure tones and chirped signals with time-varying tones.The constant modulus property results in a complex circular search space when attempting to identify interference at the receiver.To characterize these circles, latent variables are introduced.The primary contribution of this article is the computation of the maximum likelihood estimator (MLE) for key parameters, specifically the time delay and Doppler shift, in the presence of these latent variables.To calculate the MLE, we choose independent von Mises distributions with unknown parameters for the interference phases and we employ the expectation-maximization (EM) algorithm, which has demonstrated asymptotic efficiency in similar scenarios and has proven effective for N-hypersphere estimation [10].To evaluate the performance of our proposed algorithm, we compare it against the theoretical limits of time-delay and Doppler shift estimation under the following particular cases: i) we consider the scenario where no interference corrupts the GNSS signal.Note that this is the best possible scenario and the theoretical limits are provided by the Cramér-Rao bound (CRB) derived in [11], ii) the misspecified conditional model [12,13].In this scenario, the signal is corrupted by an interference but the receiver estimates the parameters of interest without considering it.Then, the time-delay and Doppler estimates are biased.Note that this is the worst possible case and the performance limits are characterized by the misspecified CRB (MCRB) derived in [14].
The fair comparison of our algorithm should be with respect to the CRB, which takes into account the parameters describing the structure of the interference.However, the corresponding derivation is intractable and we resort to the so-called modified CRB (MoCRB) [15], which is a looser bound, normally used in problems involving missing variables.
The article's structure is organized as follows: In Sect.2, we present the GNSS received signal model in the presence of a CM bandlimited interfering signal.In Sect.3, we derive a closed-form expression of the MoCRB for the parameters of interest, considering that the signal is bandlimited.This expression only depends on the baseband samples and the parameters that define the structure of the interference.Section 4 provides an in-depth description of the proposed EM algorithm, which is employed for mitigating interference under the CM hypothesis.Section 5 presents the simulation results that confirm the effectiveness of the proposed approach for two synthetic signal scenarios.Finally, Sect.6 offers the concluding remarks.

Signal model
In this article, we consider a bandlimited signal s(t), with bandwidth B, transmitted over a carrier frequency f c and traveling at the speed of light c, from a GNSS satellite to a receiver.The transmitter and receiver are assumed to be in uniform linear motion such as the distance can be modeled by a first-order d − v distance-velocity model [16].At the receiver, a narrow-band signal model is assumed and the received signal x(t) at the output of the receiver's Hilbert filter can be approximated by [11,17] where ρ and φ are the amplitude and phase of the complex coefficient α = ρe jφ ∈ C induced by the propagation characteristics, τ = d/c is the unknown propagation delay, b = v/c is the unknown Doppler shift and n(t) is a zero-mean white complex circu- lar Gaussian noise.An interfering signal I(t), unknown and bandlimited within the frequency band of interest, is also arriving at the receiver.Then, the received signal x becomes: Considering the acquisition of N = N 2 − N 1 + 1 samples at the sampling frequency F s = B = 1/T s , and assuming that the observation window [N 1 T s , N 2 T s ] is short enough to consider constant amplitude, delay and Doppler shift, the discrete signal model yields to: where µ(η) = . . ., s(kTs − τ )e −j2π f c b(kT s −τ ) , . . . . .T ∈ C N and n = [. . . ,n(kT s ), . . .T ∼ CN (0, σ 2 I N ) .Under constant modulus interference, all the components of the vector I have the same modulus A. We therefore propose the parametrization (1) In other words, each component of the vector Ĩ belongs to the complex unit circle, meaning that the vector Ĩ belongs to the complex hyper-torus of dimension N. Hence, there exists θ = θ 1 . . .θ N T ∈ (−π, π ] N

Complete likelihood
The resulting problem has the following likelihood where ε = η T , ρ, φ, A, σ 2 is the vector gathering the parameters of interest.We con- sider the angles θ k as latent (random) variables with independent prior distributions p(θ k ) .We can therefore form the joint likelihood of observed random variables x and unobserved random variables θ as p(x, θ|ε) = p(x|θ, ε)p(θ ) using ( 6) and the chosen prior p(θ) = N k=1 p(θ k ) .One way to overcome the fact that θ is unobserved is to mar- ginalize p(x, θ |ε) w.r.t.θ and to maximize this marginalized likelihood.For the particular case where p(θ k ) follows a uniform distribution over [0, 2π ] , the marginalized likelihood expression is: where B ν is the modified Bessel function of the first kind and order ν [18, Chap.3.5.4].This expression cannot be optimized wrt.ε , and one cannot derive closed-form expres- sions for the maximum likelihood estimators of the parameters in ε .Moreover, when the chosen prior over θ k is more complicated than the uniform one, the likelihood expres- sion would be even more complicated and closed-form expressions for the maximum likelihood estimators cannot be derived.One way to bypass this limit is to resort to the EM algorithm.The EM algorithm [19] can be handy to evaluate the maximum likelihood estimator of parameters when missing variables appear in the estimation framework.
The complete likelihood of the parameters ε given the observations x and missing variables θ can be expressed as: where p(x|θ, ε) is given in (6).For the prior p(θ ) , in this article we choose independent Von Mises distributions with parameter γ and κ for the interference phases θ with (4) where γ is the mean direction of the von Mises distribution, κ is the concentration parameter.In the following, we use the notation to describe the interference phases.One can note that when κ is set to 0, the uniform distribution is recovered, and hence, the results should be the same as presented in [20].This prior distribution depends on a set of two hyperparameters ϕ = {γ , κ} .These hyperparameters might be set by the user or unknown.Therefore, in the following, we will consider the general case where they are unknown and will estimate them jointly with the vector of parameters ε .Using prior ( 9) and the conditional likelihood (6), we can rewrite the complete likelihood as where �z� 2 = z H z for any z ∈ C N and where it is assumed all θ k are in (−π , π ] to avoid the indicator functions.

Modified Cramér-Rao bound
To assess the performance of the proposed method, it would be appealing to derive the Cramér-Rao Bound (CRB) of the likelihood (7).However, the corresponding derivations are intractable, mainly due to the expectation of the B 0 function, having no closed-form expression.To alleviate this problem, we propose to resort to the so-called modified CRB (MoCRB) [21], which has been designed for such problems involving missing variables.The MoCRB results in a looser bound of the asymptotic estimation performance of the parameter of interest, but with a closed-form formulation.The MoCRB of parameters of interest ε is then defined in its vector form as [15] where the matrix to be inverted is the so-called modified Fisher information matrix (MoFIM) of the vector ε , denoted F M (ε) in the following.Moreover, we have Hence, of a complex Gaussian model, and we can resort to the Slepian-Bangs formula [22, (8.34)] to find its expression (10) Recalling ε = η T , ρ, φ, A, σ 2 and given the previous formula, the FIM is where F θ (η T , ρ, φ) = F (η T , ρ, φ) is the FIM for the case without interference, F θ (σ 2 ) = N σ 4 and both are derived in [11] and independent of θ .On the other hand, To derive the MoFIM, one has to take the expected value of ( 16) wrt.θ .Given the previous results, one only requires the expected value ĨH (θ ) , which yields to [18, (3.5.25),(3.5.26)] with 1 N the vector of ones with size 1 × N .Then, after some calculations that can be found in "Appendix 1", (18) yields to with w c = 2π f c , s = . . ., s(kT s ), . . ., e φ−γ (f ) = . . ., e j(φ−γ −2π fkT s ) , . . .and D = diag(N 1 , . . ., N 2 ) .Finally, the MoFIM is Finally, from the matrix inversion lemma [23, 14.11-(a)] we have the MoCRB expression for the parameters of interest (15) Note that when κ = 0 , i.e., when the prior distribution is uniform, the second term van- ishes due to the term B 1 (κ) B 0 (κ) which is 0 when κ = 0 , and the MoCRB becomes the CRB in the absence of interference.

EM approach for interference mitigation under the CM hypothesis
In this section, we introduce the proposed EM algorithm.The EM algorithm iterates between the expectation (E) and the maximization (M) steps to obtain a maximum of the likelihood function: • E-step: the derivation of the function where t represents the iteration index.

Derivation of the Q function
At the t + 1 iteration, the E-step approximates the loglikelihood (which is to be maxi- mized) around the parameters ε (t) and ϕ (t) .This approximation is given by: where ε (t) , reps.ϕ (t) , represent the current value for the set of parameters, resp.hyper- parameters.From (12), we have (22) , where arg(.): C → (−π , π ] the argument of a complex number.The only terms depending on θ in (26) are in the sum of cosines, leading to

Conditional distribution
This expected value considers the distribution of θ |x, ε (t) , ϕ (t) .The distribution of θ |x, ε (t) , ϕ (t) can be expressed as (see "Appendix 2" for details) and To express the terms in this distribution, we define x k as the k-th component of the vec- tor x and µ k (η) the k-th component of the vector µ(η) , and then,

M-step
The second step of the EM algorithm is to maximize (37) w.r.t.ε and ϕ .Note that this equation is the sum of two independent terms in ε and ϕ , so the optimization can be done independently.

Optimization with respect to ρ, ϕ, η
The second term in function (41) is independent of ρ, ϕ and η , and then, we can recast the equation as where K ′′ represents the terms independent of ρ, ϕ and η .Let us denote S = span(B) , with B a matrix, which is the linear span of the set of its column vectors.
Note that maximizing the previous expression wrt.ρ, ϕ yields to the following update equations: , where the Pythagora's theorem has been applied.Then, maximizing Q(ε|ε (t) , ϕ (t) ) w.r.t.η yields to the fol- lowing update equation

M-step: optimization with respect to A
To carry out the optimization, we first remind that the objective function (41) can be expressed as: where K ′′′ represents the terms independent of A. Differentiating this expression w.r.t.A yields to We can note that can be injected in (51) to give the update equation Note that (53) can be injected into (49) to provide the update equation ( 46) (47) . (49) 4.4 M-step: optimization with respect to κ and γ In the case where the hyperparameters are to be estimated jointly with the model parameters, we can update their values maximizing Q(ϕ|ε (t) , ϕ (t) ) , previously defined in (42).

Let us denote
Then, C t = R t cos γ t and S t = R t sin γ t .Moreover, yielding From the previous equation, we can recognize the loglikelihood of a von Mises distribution as in [18, (5.3.1)], which leads the update equations where the function A(κ) = B 1 (κ)/B 0 (κ) has no closed form, but can be approximated using an iterative algorithm [24].

Implemented algorithm
For t = 0 , i.e., for initialization, we might compute η (0) T , ρ (0) , φ (0) , σ 2 (0) thanks to the standard MLE [11].In other words, (54) Once the interference is detected, due to the fact that the interference is much more powerful than the GNSS signal masked in the noise, we keep the MLE estimators for η (1) , ρ (1) ϕ (1) and (σ 2 ) (1) , and we initialize the interference power as If we do not know the parameters for the prior distribution, we can set γ (1) = 0 and κ (1) = 0 .Those values indicate that the prior distribution is uniform.Then, we compute the Von Mises parameters, γ (t)  k and κ (t) k , and trigonometric moments w (t) k and a t , using Eqs.( 32), (33), and (40), respectively.Finally, we update the parameters of interest in the order η, A, ρ, ϕ and σ 2 with Eqs. (43), (47), (48), (53),and (54), respectively.Note that Eqs. ( 54) and (63) do not have closed-form expressions.However, a grid search approach can be used to maximize the corresponding functions.Also, a little simplification can be made to optimize (49), using (54) where A (t+1) is replaced by A (t) to update η .With this approach, the same function can be used to solve both (54) and (63).

Simplification of the algorithm based on an uniform prior
The algorithm presented in this article is a generalization of the algorithm presented in [20].In that article, the prior distribution is considered to be a uniform distribution, i.e. κ (t) = 0 and γ (t) can be undefined.Therefore, as demonstrated in "Appendix 2", k and M-step presented in Sect.4.4 are not required.Once these simplifications have been made, the algorithm presented in [20] is the same as the one proposed in Sect.4.5.

Experiments
In order to evaluate the performance of the algorithms presented in Sects.4.5 and 4.6, we consider two possible scenarios.

• Scenario 1:
Let us consider the case where a GPS L1 C/A signal [1] is attacked by a jammer that is generating a linear frequency modulation (LFM) signal [25], which is defined as (63) (67) where β c is the chirp rate and A is the amplitude.For this particular scenario, we set the waveform period as T = N • T s , i.e., equal to the integration time.The instantane- ous frequency is f (t) = 1 2π d dt πβ c t 2 = β c t , and therefore, the waveform bandwidth is B c = β c T .We consider the case where, after the Hilbert filter, the chirp is located at the baseband frequency, i.e., the central frequency of the chirp is f i = 0 .Then, the waveform can be rewritten as: We set the chirp bandwidth B c = 1 MHz, with initial phase φ = 0 , amplitude A = 40 .In Fig. 1, we illustrate the power spectral density (PSD) of the linear chirp considered in scenario 1.We would also like to point out that the phase distribution for this case is a uniform distribution.

• Scenario 2:
We consider a GPS L1 C/A signal attacked by a jammer generating a nonlinear frequency modulation [26], which is defined as where ϕ(β c ; t) is a nonlinear function with β c a parameter that controls the bandwidth of the chirp.For this particular scenario, we set ϕ(β c ; T ; t) = sin π T πβ c t .Note that the chirp is located at the baseband frequency, i.e., the central frequency of the chirp is f i = 0 .For this scenario, we set β c = 51 , T = 1 ms (which it is the duration of the (68)  The root-mean-squared error (RMSE) for the parameters of interest η T is displayed in Figs. 4 and 5 for scenario 1 and in Figs. 6 and 7 for scenario 2. These figures depict how the RMSE varies with respect to the signal-to-noise ratio (SNR) at the output of the matched filter, labeled as SNR OUT .The evaluation is conducted under the following conditions: a GNSS receiver with a sampling frequency of F s = 4 MHz and integration times of T = 1 ms.
The EM algorithm incorporates a stopping criterion that assesses the change in noise variance at each iteration.The maximum iteration limit is set to 15, and 1000 Monte • the √ CRB (as referred to in [11]), which signifies the asymptotic estimation perfor- mance of the parameters in the absence of interference • the MCRB + Bias 2 , representing the asymptotic estimation performance of the parameters when the receiver is unaware of the presence of interference (as discussed in [13,27]).This includes the root MSE of the misspecified maximum likelihood esti- effectiveness of this EM algorithm by evaluating its RMSE for time-delay and Doppler parameters.The evaluation encompasses scenarios involving chirp interference jamming and a GPS L1 C/A signal.The results clearly illustrate the strong performance of the proposed algorithm.Finally, we would like to point out that numerous array processing solutions leverage Riemannian optimization to mitigate interference while adhering to constant modulus constraints [28,29].These methods could represent a promising advancement in the search for low-complexity interference mitigation algorithms for GNSS.

Appendix 1: MoFIM derivation for bandlimited signals
In this appendix, we focus on the derivation of the MoFIM.In order to do that, we need to derive a closed-form expression of the elements within the vector F M A, η T , ρ, φ .
For the element F M (A, ρ) , it is required to compute: and applying the Nyquist-Shannon theorem for bandlimited signals, we have where and Thus, the closed-form expression of the element yields to Similarly, to compute the element F M (A, φ) , it is required to compute: (71) Then, it is simple to check that closed form yields to Now, we compute the elements within the vector F M A, η T .Then, we need to compute where s (1) = ds(t) dt .Then, and Applying the Nyquist-Shannon theorem for bandlimited signals, we have where (77)   (101 T ∈ C N with η = (τ , b) and k ∈ {N 1 , . . ., N 2 } , I = [. . . ,I(kT s ), .
B = B(B H B) −1 B H and ⊥ B = I − B are the orthogonal projectors over S and S ⊥ , respectively.Then, and

Fig. 1
Fig. 1 PSD of centered linear chirp signal of bandwidth B c = 1 MHz and amplitude A = 40 .The sampling frequency is set to F s = 4 MHz

Fig. 2 Fig. 3
Fig. 2 PSD of centered nonlinear chirp signal of bandwidth B c = 0.16 MHz and amplitude A = 40 .The sampling frequency is set to F s = 4 MHz

1 Fig. 4 4 1 Fig. 5
Fig. 4 RMSE for time-delay estimation of the GPS L1 C/A signal received along with a centered LFM chirp signal of bandwidth B c = 1 MHz and amplitude A = 40 .The sampling frequency is set to F s = 4 MHz

2 Fig. 6 4 2 Fig. 7
Fig. 6 RMSE for time-delay estimation of the GPS L1 C/A signal received along with a centered nonlinear chirp signal with β c = 51 , T = 1 ms and amplitude A = 40 .The sampling frequency is set to F s = 4 MHz