CN104199067A - 多径环境下全球卫星导航系统接收机无模糊处理方法 - Google Patents

Abstract

本发明公开了多径环境下全球卫星导航系统接收机无模糊处理方法，其步骤为：根据BOC调制信号的统计特性，设计旁瓣抑制相关器，使得接收机跟踪到主瓣而不是旁瓣、输出具有单峰特性并且相关函数随时延的增加单调递减；建立多径环境下GNSS信号传播模型，接收机首先根据相关器输出判断接收信号是否存在视距信号，若存在视距信号，采用最大似然估计得到实时位置信息，运用牛顿迭代法更新幅度、时延、相位等参数，若不存在视距信号，则利用最先到达多径分量的统计特性修正视距传播时延。本发明设计的多径环境下GNSS无模糊处理接收机有效克服了多径衰落和BOC信号自相关函数多峰特性的影响，定位精度可以逼近克拉美罗界。

Description

多径环境下全球卫星导航系统接收机无模糊处理方法

技术领域

本发明涉及全球卫星导航系统定位技术领域，具体涉及多径环境下全球卫星导航系统接收机无模糊处理方法。

背景技术

全球卫星导航系统(global navigation satellite system，GNSS)可以全天候、广覆盖、实时地为海陆空各领域提供通信和导航服务。随着无线通信技术的发展，卫星导航系统凭借其在信息技术领域的独特优势，正逐步取代传统的地面无线电导航、天文测量和大地测绘技术，成为人类社会生产生活中普遍采用的导航定位手段。对于下一代GNSS，二进制偏移载波(binary offset carrier，BOC)作为一种新型的波形信号广泛应用于现代全球定位系统(global positioning system，GPS)，伽利略定位系统(Galileo positioning system)和北斗卫星导航系统(Beidou satellite navigation system，BDS)。在BOC调制中，信号由伪随机序列与方波亚载波相乘得到。与当前导航频段采用的二进制相移键控(binaryphase shift keying，BPSK)相比，BOC信号实现频谱分离的同时，提高了接收机的码跟踪精度和抗多径干扰能力。

对于GNSS接收机，通常认为BOC调制比BPSK调制具有更强的抗多径衰落性能。早期简单的码跟踪方法采用时延锁相环(delay lock loop，DLL)预测每个相关采样的时延。多径估计延迟锁相环(multipath estimation delay lock loop，MEDLL)技术比较具有不同幅度、时延、相位等参数的自相关函数，估计得出最符合接收序列的本地信号，但该方法在短时延和密集多径的情况下效果并不理想。

然而，BOC信号体制的主要缺陷在于自相关函数的多峰特性，这将在码跟踪环节出现潜在的模糊性问题。常见的无模糊捕获和跟踪方法包括BPSK-like技术，Bump-jumping技术和旁瓣抑制技术。BPSK-like技术和Bump-jumping技术通过比较自相关函数主峰和副峰幅度，确定本地信号是否与接收序列同步，但该方法不适用于信噪比较低的情况。旁瓣抑制技术通过生成不同与接受序列的本地信号，利用多个相关器和滤波器有效去除旁瓣对码跟踪精度的影响，但该方法降低了相关函数主峰的尖锐程度，运算复杂度也显著提高。

发明内容

发明目的：为了提高全球卫星导航系统接收机的定位精度，本发明提供多径环境下全球卫星导航系统接收机无模糊处理方法。

技术方案：为解决上述技术问题，本发明提供的多径环境下全球卫星导航系统接收机无模糊处理方法，其实现包括以下步骤：

(一)根据BOC调制信号的统计特性，设计旁瓣抑制相关器，使得接收机跟踪到主瓣而不是旁瓣、输出具有单峰特性，并且相关函数随时延的增加单调递减；

(二)建立多径环境下GNSS信号传播模型，定义接收到的复GNSS信号和本地载波的无模糊相关函数|R_rs(τ)|；

(三)接收机根据相关器输出判断接收信号是否存在视距信号，若存在视距信号，接收机采用最大似然估计得到实时位置信息，运用牛顿迭代法更新幅度、时延、相位参数，若不存在视距信号，则利用最先到达多径分量的统计特性修正理想状态下视距传播时延；

(四)设计多径环境下GNSS无模糊处理接收机，采用牛顿迭代计算多径时延，定位精度可以逼近克拉美罗界CRB。

所述步骤(一)中设计旁瓣抑制相关器采用波形为正弦信号的亚载波代替传统的方波亚载波，包括如下步骤：

(1)将接收到的BOC调制信号r(t)与本地同相正弦亚载波信号s_I(t)进行相关运算得到R_I(τ)；

(2)将接收到的BOC调制信号r(t)与本地正交余弦亚载波信号s_Q(t)进行相关运算得到R_Q(τ)；

(3)将同相相关函数R_I(τ)和正交相关函数R_Q(τ)进行平方和运算得到无模糊相关函数R(τ)。

所述步骤(二)多径环境下解调后接收复GNSS信号可以表示为

式中，A(t)表示视距分量的幅度，τ表示信号传播时间，表示载波相位偏移，P表示多径信号的数目，a_p(t)表示第p条多径分量的幅度，τ_p表示第p条多径分量的时延，表示第p条多径分量的相位偏移，n(t)表示复噪声。

假设信号的幅度和传播时延在单位观测时间内保持不变，在只有视距传播信号的情况下，传输过程可以看作加性高斯白噪声AWGN卫星信道，在存在多径反射信号的情况下，接收信号可以简化为视距分量和最先到达的非视距分量，根据中心极限定理，密集多径条件下接收信号模型可以简化为单路径传播，即接收到的复GNSS信号可以简化为

r(t)＝As（t-τ)+as(t-ε)+n(t)

式中，A表示视距相量，τ信号传播时间，a表示最先到达的非视距幅度相量，ε表示最先到达的非视距幅度分量；

定义接收到的复GNSS信号和本地载波的无模糊相关函数为

$R_{\mathrm{rs}}\left(\mathrm{\τ}\right)={\∫}_{t_0}^{t_0+T}r\left(t\right)s^{*}(t-\mathrm{\τ})\mathrm{dt}.$

所述步骤(三)包括如下步骤：

若无模糊相关函数|R_rs(τ)|大于或等于门限δ，则存在视距传播信号，此时接收信号r(t)＝As(t-τ)+as(t-ε)+n(t)，接收机采用最大似然估计得到实时位置信息，运用牛顿迭代法更新视距分量和最先到达的多径分量的幅度、时延、相位；

若无模糊相关函数|R_rs(τ)|小于门限δ，则不存在视距传播信号，此时接收信号r(t)＝as(t-ε)+n(t)，接收机首先估计得到最先到达的多径分量的幅度、时延、相位等参数，根据多径分量到达时间间隔服从指数分布的特性，采用最大似然估计修正理想状态下视距分量传播时延。

所述步骤(四)包括如下步骤：

多径环境下传播时延的克拉美罗界CRB表示为

${\mathrm{CRB}}_2\left(\mathrm{\τ}\right)=\frac{{(A-{\mathrm{BD}}^{-1}C)}_{\mathrm{2,2}}}{{(A-{\mathrm{BD}}^{-1}C)}_{\mathrm{1,1}}{(A-{\mathrm{BD}}^{-1}C)}_{\mathrm{2,2}}-{(A-{\mathrm{BD}}^{-1}C)}_{\mathrm{1,2}}^2}$

式中，

${(A-{\mathrm{BD}}^{-1}C)}_{\mathrm{1,1}}=4{\mathrm{\π}}^2{\left|A\right|}^2{f}_{\mathrm{Gabor}}^2-\frac{R{\left|A\right|}^2{R}^{\′\′2}\left(\mathrm{\δ}\right)}{R^2-R^2\left(\mathrm{\δ}\right)}$

${(A-{\mathrm{BD}}^{-1}C)}_{\mathrm{1,2}}=-\mathrm{Re}\left\{A{a}^{*}\right\}[R^{\′\′}\left(\mathrm{\δ}\right)+\frac{R\left(\mathrm{\Δ}\right)R^{\′2}\left(\mathrm{\δ}\right)}{R^2-R^2\left(\mathrm{\δ}\right)}]$

${(A-{\mathrm{BD}}^{-1}C)}_{\mathrm{2,2}}=4{\mathrm{\π}}^2{\left|a\right|}^2{f}_{\mathrm{Gabor}}^2-\frac{R{\left|a\right|}^2{R}^{\′\′2}\left(\mathrm{\δ}\right)}{R^2-R^2\left(\mathrm{\δ}\right)}$

当接收信号不存在视距分量时，接收信号可以看作最先到达的多径分量和噪声的叠加，此时，多径环境下概率密度函数为

$p(r;\mathrm{\θ})=\frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\σ}}\exp \{-\frac{1}{2{\mathrm{\σ}}^2}{\∫}_{t_0}^{t_0+T}{|r\left(t\right)-\mathrm{as}(t-\mathrm{\ϵ})|}^2\mathrm{dt}\}$

同理可得多径分量的幅度为

$\hat{a}=\frac{R_{\mathrm{rs}}\left(\mathrm{\ϵ}\right)}{R}$

多径时延的迭代公式为

${\widehat{\mathrm{\ϵ}}}_{i+1}={\widehat{\mathrm{\ϵ}}}_i-\frac{\mathrm{Re}\{R_{\mathrm{rs}}^{*}\left({\widehat{\mathrm{\ϵ}}}_i\right)\frac{\∂R_{\mathrm{rs}}\left({\widehat{\mathrm{\ϵ}}}_i\right)}{\∂\mathrm{\ϵ}}}{\mathrm{Re}\left\{R_{\mathrm{rs}}^{*}\left({\widehat{\mathrm{\ϵ}}}_i\right)\frac{{\∂}^2{R}_{\mathrm{ra}}\left({\widehat{\mathrm{\ϵ}}}_i\right)}{\∂{\mathrm{\ϵ}}^2}\right\}+{\left|\frac{\∂R_{\mathrm{rs}}\left({\widehat{\mathrm{\ϵ}}}_i\right)}{\∂\mathrm{\ϵ}}\right|}^2}$

根据早迟门同步法，上式可以改写为

${\widehat{\mathrm{\ϵ}}}_{i+1}={\widehat{\mathrm{\ϵ}}}_i-\frac{\mathrm{hRe}\left\{R_{\mathrm{rs}}^{*}\left({\widehat{\mathrm{\ϵ}}}_i\right)\right[R_{\mathrm{rs}}({\widehat{\mathrm{\ϵ}}}_i+h)-R_{\mathrm{rs}}({\widehat{\mathrm{\ϵ}}}_i-h)\left]\right\}}{2\mathrm{Re}\left\{R_{\mathrm{rs}}^{*}\left({\widehat{\mathrm{\ϵ}}}_i\right)\right[R_{\mathrm{rs}}({\widehat{\mathrm{\ϵ}}}_i+h)+R_{\mathrm{rs}}({\widehat{\mathrm{\ϵ}}}_i-h)-2R_{\mathrm{rs}}\left({\widehat{\mathrm{\ϵ}}}_i\right)\left]\right\}}$

多径分量到达时间间隔服从指数分布，将最先到达的多径时延ε₁，ε₂，...，ε_n作为样本，采用最大似然函数

$L({\mathrm{\ϵ}}_1,{\mathrm{\ϵ}}_2,...,{\mathrm{\ϵ}}_n;\widehat{\mathrm{\τ}})=\underset{\mathrm{\τ}\∈T}{\max }\underset{i=1}{\overset{n}{\mathrm{\Π}}}f({\mathrm{\ϵ}}_i;\mathrm{\τ})$

式中，f(ε₁；τ)表示指数分布函数，利用多径时延得到τ的最大似然估计量为

τ＝min{ε₁，ε₂，…，ε_n}。

当接收信号存在视距分量时，多径环境下接收信号的概率密度函数为

$p(r;\mathrm{\θ})=\frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\σ}}\exp \{-\frac{1}{2{\mathrm{\σ}}^2}{\∫}_{t_0}^{t_0+T}{|r\left(t\right)-\mathrm{As}(t-\mathrm{\τ})-\mathrm{as}(t-\mathrm{\ϵ})|}^2\mathrm{dt}\}$

式中，θ＝(τ，ε，A，a)∈Θ是最大似然估计的四个参数，σ²是AWGN的平均功率。假设积分时间T是GNSS信号周期的整数倍，最大似然函数可以表示为

$\begin{array}{c}\mathrm{\Γ}\left(\mathrm{\θ}\right)=-{\∫}_{t_0}^{t_0+T}{\left|r\left(t\right)\right|}^2\mathrm{dt}+2\mathrm{Re}\left\{A^{*}{\∫}_{t_0}^{t_0+T}r\left(t\right)s^{*}(t-\mathrm{\τ})\mathrm{dt}\right\}-{\left|A\right|}^2{\∫}_{t_0}^{t_0+T}{\left|s(t-\mathrm{\τ})\right|}^2\mathrm{dt}\\ +2\mathrm{Re}\left\{a^{*}{\∫}_{t_0}^{t_0+T}r\left(t\right)s^{*}(t-\mathrm{\ϵ})\mathrm{dt}\right\}-{\left|a\right|}^2{\∫}_{t_0}^{t_0+T}{\left|s(t-\mathrm{\ϵ})\right|}^2\mathrm{dt}-2\mathrm{Re}\left\{{\mathrm{Aa}}^{*}\right\}{\∫}_{t_0}^{t_0+T}s(t-\mathrm{\τ})s^{*}(t-\mathrm{\ϵ})\mathrm{dt}\end{array}$

式中，R_ss(ε-τ)表示s(t-τ)和s(t-ε)的自相关函数，它是一个仅取决于ε和τ的实数

$R_{\mathrm{ss}}(\mathrm{\ϵ}-\mathrm{\τ})={\∫}_{t_0}^{t_0+T}s(t-\mathrm{\τ})s^{*}(t-\mathrm{\ϵ})\mathrm{dt}$

幅度相量为

$\left\{\begin{array}{c}\hat{A}=\frac{R{R}_{\mathrm{rs}}\left(\mathrm{\τ}\right)-R_{\mathrm{ss}}(\mathrm{\ϵ}-\mathrm{\τ})R_{\mathrm{rs}}\left(\mathrm{\ϵ}\right)}{R^2-R_{\mathrm{ss}}^2(\mathrm{\ϵ}-\mathrm{\τ})}\\ \hat{a}=\frac{R{R}_{\mathrm{rs}}\left(\mathrm{\ϵ}\right)-R_{\mathrm{ss}}(\mathrm{\ϵ}-\mathrm{\τ})R_{\mathrm{rs}}\left(\mathrm{\τ}\right)}{R^2-R_{\mathrm{ss}}^2(\mathrm{\ϵ}-\mathrm{\τ})}\end{array}\right.$

采用牛顿迭代估计传播时延τ和ε

$\left\{\begin{array}{c}{\widehat{\mathrm{\τ}}}_{i+1}={\widehat{\mathrm{\τ}}}_i-\frac{\mathrm{hRe}\left\{R_{\mathrm{rA}}^{*}\left({\widehat{\mathrm{\τ}}}_i\right)\right[R_{\mathrm{rA}}({\widehat{\mathrm{\τ}}}_i+h)-R_{\mathrm{rA}}({\widehat{\mathrm{\τ}}}_i-h)\left]\right\}}{2\mathrm{Re}\left\{R_{\mathrm{rA}}^{*}\left({\widehat{\mathrm{\τ}}}_i\right)\right[R_{\mathrm{rA}}({\widehat{\mathrm{\τ}}}_i+h)+R_{\mathrm{rA}}({\widehat{\mathrm{\τ}}}_i-h)-2R_{\mathrm{rA}}\left({\widehat{\mathrm{\τ}}}_i\right)\left]\right\}}\\ {\widehat{\mathrm{\ϵ}}}_{i+1}={\widehat{\mathrm{\ϵ}}}_i-\frac{\mathrm{hRe}\left\{R_{\mathrm{ra}}^{*}\left({\widehat{\mathrm{\ϵ}}}_i\right)\right[R_{\mathrm{ra}}({\widehat{\mathrm{\ϵ}}}_i+h)-R_{\mathrm{ra}}({\widehat{\mathrm{\ϵ}}}_i-h)]\}}{2\mathrm{Re}\left\{R_{\mathrm{ra}}^{*}\left({\widehat{\mathrm{\ϵ}}}_i\right)\right[R_{\mathrm{ra}}({\widehat{\mathrm{\ϵ}}}_i+h)+R_{\mathrm{ra}}({\widehat{\mathrm{\ϵ}}}_i-h)-2R_{\mathrm{ra}}\left({\widehat{\mathrm{\ϵ}}}_i\right)\left]\right\}}\end{array}\right.$

有益效果：与现有技术相比，本发明所述的多径环境下全球卫星导航系统接收机无模糊处理方法定位精度显著提高；根据BOC调制信号的统计特性，本发明设计了一种简单有效的旁瓣抑制相关器，使得接收机的输出具有单峰特性，相关函数随时延的增加单调递减，有效克服由于BOC信号自相关函数的多峰特性而导致的潜在模糊性。针对多径环境下GNSS信号传播特点，本发明根据相关器输出判断接收信号是否存在视距信号，若存在视距信号，接收机采用最大似然估计得到实时位置信息，同时，运用牛顿迭代法更新幅度、时延、相位等参数，以适应外界环境的变化，若不存在视距信号，则利用最先到达的多径时延修正理想状态下视距传播时延。因此，本发明提供的技术方案实现了多径环境下GNSS接收信号的无模糊处理，定位精度可以逼近克拉美罗界。

附图说明

图1是本发明中的无模糊处理接收机原理框图。

图2是本发明中的归一化相关函数随时延变化仿真图。

图3是本发明中的S曲线随估计误差变化仿真图。

图4是本发明中的多径环境下GNSS信号传播模型示意图。

具体实施方式

下面结合附图对本发明作更进一步的说明。

如图1所示，本发明的多径环境下全球卫星导航系统接收机无模糊处理方法，对于下一代GNSS，现代全球定位系统，伽利略定位系统和北斗卫星导航系统已经采用不同参数的BOC调制信号。采用不同参数的BOC调制信号的思想主要是为了实现卫星信道不同GNSS信号的频谱分离。BOC信号可以描述为BPSK信号和方波亚载波的乘积，BOC发送信号的复包络可以表示为

$s\left(t\right)=e^{-\mathrm{j\θ}}\underset{k}{\mathrm{\Σ}}{a}_k{\mathrm{\μ}}_{n{T}_S}(t-\mathrm{kn}{T}_s-t_0)C_{T_s}(t-t_0)---\left(1\right)$

式中，a_k表示调制信号的扩频码，a_k∈{+1，-1}，C_Ts(t)表示周期为2T_s，μ_nTs(t)表示持续时间为nTs的矩形脉冲扩频符号，调制阶数n表示扩频码保持不变时半周期亚载波个数，θ和t₀分别表示复包络相对于参考点的相位偏移和时延。通常情况下，BOC信号可以表示为BOC(f_s，f_c)，亚载波频率f_s＝1/2T_s，传码率f_c＝1/nT_s，单位均是1.023MHz的倍数。

假设BOC扩频符号独立同分布。当n为偶数时，BOC调制的归一化基带功率谱密度可以表示为

$S\left(f\right)=\frac{f_c}{{\mathrm{\π}}^2{f}^2}{\tan }^2\left(\frac{\mathrm{\πf}}{{2f}_s}\right){\sin }^2\left(\frac{\mathrm{\πf}}{f_c}\right)---\left(2\right)$

当n为奇数时，归一化基带功率谱密度可以表示为

$S\left(f\right)=\frac{f_c}{{\mathrm{\π}}^2{f}^2}{\tan }^2\left(\frac{\mathrm{\πf}}{{2f}_s}\right){\cos }^2\left(\frac{\mathrm{\πf}}{f_c}\right)---\left(3\right)$

对于无限带宽的理想BOC信号，自相关函数由分段线性的主瓣和多个旁瓣组成。自相关函数幅度正负峰值的总数为2n-1。L表示自相关函数的指数，l＝0，1，...，n-1，峰值间隔为T_s。对于时延|τ|∈[lT_s，(l+1)T_s]，归一化自相关函数可以表示为

$R\left(\mathrm{\τ}\right)={(-1)}^l(2l+1-\frac{{2l}^2+2l}{n}-\frac{2n-2l-1}{{\mathrm{nT}}_s}|\mathrm{\τ}\left|\right)---\left(4\right)$

Gabor带宽f_Gabor定义为带限信号的均方根带宽，可以表示为

$f_{\mathrm{Gabor}}=\sqrt{{\∫}_{-B}^B{f}^{2}S\left(f\right)\mathrm{df}}---\left(5\right)$

式中，B是接收机前端的带宽。F_Gabor决定了码跟踪精度的下界。F_Gabor值越大，接收机获得的码跟踪精度越高。

由于BOC调制信号的自相关函数具有多峰特性，必须保证接收机跟踪到主瓣而不是旁瓣。本发明将接收到的BOC调制信号r(t)与本地正弦波信号进行相关运算，提供了一种简单有效的旁瓣去除技术，从而可以通过软件的方法去除码跟踪多峰效应。如图2所示，对于理想BOC信号，R(τ)具有单峰特性，并且随着时延|τ|的增加单调递减，确保最大似然估计通过牛顿迭代得到全局最优解。

具体说来，首先生成本地同相正弦亚载波信号

$s_I\left(t\right)=e^{-\mathrm{j\θ}}\underset{k}{\mathrm{\Σ}}{a}_k{\mathrm{\μ}}_{n{T}_s}(t-{\mathrm{knT}}_s-t_0)\text{sin}\frac{\mathrm{\π}}{T_s}(t-t_0)---\left(6\right)$

r(t)和本地同相正弦亚载波信号的相关函数定义为

$R_I\left(\mathrm{\τ}\right)={\∫}_{t_0}^{t_0+T}r\left(t\right)s_I^{*}(t-\mathrm{\τ})\mathrm{dt}---\left(7\right)$

式中，T表示相关运算的积分时间，(·)^*表示共轭运算。

同理，生成本地正交余弦亚载波信号

$s_Q\left(t\right)=e^{-\mathrm{j\θ}}\underset{k}{\mathrm{\Σ}}{a}_k{\mathrm{\μ}}_{n{T}_s}(t-{\mathrm{knT}}_s-t_0)\text{cos}\frac{\mathrm{\π}}{T_s}(t-t_0)---\left(8\right)$

r(t)和本地正交余弦亚载波信号的相关函数定义为

$R_Q\left(\mathrm{\τ}\right)={\∫}_{t_0}^{t_0+T}r\left(t\right)s_Q^{*}(t-\mathrm{\τ})\mathrm{dt}---\left(9\right)$

无模糊相关函数可以通过求解式(7)和式(9)的平方和得到

$R\left(\mathrm{\τ}\right)=R_I^2\left(\mathrm{\τ}\right)+R_Q^2\left(\mathrm{\τ}\right)---\left(10\right)$

式中，归一化同相相关函数R_I(τ)可以表示为

$R_I\left(\mathrm{\&tau;}\right)=\left\{\begin{array}{cc}\frac{1}{2n}+\frac{2n-1}{2n}\cos \frac{\mathrm{\&pi;}}{T_s}\mathrm{\&tau;}& 0\&le;\left|\mathrm{\&tau;}\right|<T_s\\ -\frac{1}{2n}+\frac{2n-3}{2n}\cos \frac{\mathrm{\&pi;}}{T_s}\mathrm{\&tau;}& T_s\&le;\left|\mathrm{\&tau;}\right|<{2T}_s\\ \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;\\ \frac{{(-1)}^l}{2n}+\frac{2n-2l-1}{2n}\cos \frac{\mathrm{\&pi;}}{T_s}\mathrm{\&tau;}& {\mathrm{lT}}_s\&le;\left|\mathrm{\&tau;}\right|<(l+1)T_s\\ \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;\\ \frac{{(-1)}^{n-1}}{2n}+\frac{1}{2n}\cos \frac{\mathrm{\&pi;}}{T_s}\mathrm{\&tau;}& (n-1)T_s\&le;\left|\mathrm{\&tau;}\right|<{\mathrm{nT}}_s\\ 0& \left|\mathrm{\&tau;}\right|\&GreaterEqual;{\mathrm{nT}}_s\end{array}\right.---\left(11\right)$

归一化正交相关函数R_Q(τ)可以表示为

$R_Q\left(\mathrm{\&tau;}\right)=\left\{\begin{array}{cc}\frac{2n-1}{2n}\sin \frac{\mathrm{\&pi;}}{T_s}\mathrm{\&tau;}& 0\&le;\left|\mathrm{\&tau;}\right|<T_s\\ \frac{2n-3}{2n}\sin \frac{\mathrm{\&pi;}}{T_s}\mathrm{\&tau;}& T_s\&le;\left|\mathrm{\&tau;}\right|<2T_s\\ \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;\\ \frac{2n-2l-1}{2n}\sin \frac{\mathrm{\&pi;}}{T_s}\mathrm{\&tau;}& l{T}_s\&le;\left|\mathrm{\&tau;}\right|<(l+1)T_s\\ \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;\\ \frac{1}{2n}\sin \frac{\mathrm{\&pi;}}{T_s}\mathrm{\&tau;}& (n-1)T_s\&le;\left|\mathrm{\&tau;}\right|<n{T}_s\\ 0& \left|\mathrm{\&tau;}\right|\&GreaterEqual;n{T}_s\end{array}\right.---\left(12\right)$

于是，本发明提供的无模糊归一化相关函数R(τ)为

$R\left(\mathrm{\&tau;}\right)=\left\{\begin{array}{cc}\frac{2n^2-2n+1}{2n^2}+\frac{2n-1}{2n^2}\cos \frac{\mathrm{\&pi;}}{T_s}\mathrm{\&tau;}& 0\&le;\left|\mathrm{\&tau;}\right|<T_s\\ \frac{2n^2-6n+5}{2n^2}-\frac{2n-3}{2n^2}\cos \frac{\mathrm{\&pi;}}{T_s}\mathrm{\&tau;}& T_s\&le;\left|\mathrm{\&tau;}\right|<2T_s\\ \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;\\ \frac{2{(n-l)}^2-2(n-l)+1}{2n^2}+\frac{{(-1)}^l(2n-2l-1)}{2n^2}\cos \frac{\mathrm{\&pi;}}{T_s}\mathrm{\&tau;}& l{T}_s\&le;\left|\mathrm{\&tau;}\right|<(l+1)T_s\\ \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;\\ \frac{1}{2n^2}+\frac{{(-1)}^{n-1}}{2n^2}\cos \frac{\mathrm{\&pi;}}{T_s}\mathrm{\&tau;}& (n-1)T_s\&le;\left|\mathrm{\&tau;}\right|<n{T}_s\\ 0& \left|\mathrm{\&tau;}\right|\&GreaterEqual;n{T}_s\end{array}\right.---\left(13\right)$

S曲线通常用来评价接收机的码跟踪性能，其表达式为

S(Δ)＝|R(Δ-h)|²-|R(Δ+h)|²   (14)

式中，Δ表示传播时延的估计误差，h是一个很小的正数。无模糊处理接收机的S曲线性能如图3所示，仿真结果表明与传统的BPSK、BOC调制相比，无模糊处理方法具有较好的码跟踪性能。

如图4所示，对于多径传播，接收信号包括视距信号，多径反射信号和噪声。在卫星移动通信中，忽略多普勒频移的影响，解调后接收复GNSS信号可以表示为

式中，A(t)表示视距分量的幅度，τ表示视距信号传播时间，表示载波相位偏移，P表示多径信号的数目，a_p(t)表示第p条多径分量的幅度，τ_p表示第p条多径分量的时延，表示第p条多径分量的相位偏移，n(t)表示复噪声。假设视距信号的幅度和传播时延在单位观测时间内保持不变。在只有视距传播信号的情况下，传输过程可以看作加性高斯白噪声(additive white Gaussian noise，AWGN)卫星信道。在存在多径反射信号的情况下，接收信号可以简化为视距分量和最先到达的非视距分量。根据中心极限定理，密集多径条件下接收信号模型可以简化为单路径传播。综上所述，接收到的复GNSS信号可以简化为

r(t)＝As(t-τ)+as(t-ε)+n(t)     (16)

式中，A表示视距相量，τ信号传播时间，a表示最先到达的非视距幅度相量，ε表示最先到达的非视距幅度分量。为了表述方便，定义接收到的复GNSS信号和本地载波的相关函数为

$R_{\mathrm{rs}}\left(\mathrm{\&tau;}\right)={\&Integral;}_{t_0}^{t_0+T}r\left(t\right)s^{*}(t-\mathrm{\&tau;})\mathrm{dt}---\left(17\right)$

对于接收到的复GNSS信号，接收机首先需要根据相关器的输出判断是否存在视距传播信号。当存在视距传播信号时，视距分量的到达时间应连续变化，多径分量的到达时间间隔服从指数分布。基于此，本发明采用如下判决准则：

若|R_rs(τ)|≥δ，则存在视距传播信号，此时接收信号r(t)＝As(t-τ)+as(t-ε)+n(t)；

若|R_rs(τ)|＜δ，则不存在视距传播信号，此时接收信号r(t)＝as(t-ε)+n(t)。

当接收信号存在视距分量时，多径环境下接收信号的概率密度函数为

$p(r;\mathrm{\&theta;})=\frac{1}{\sqrt{2\mathrm{\&pi;}}\mathrm{\&sigma;}}\exp \{-\frac{1}{{2\mathrm{\&sigma;}}^2}{\&Integral;}_{t_0}^{t_0+T}{|r\left(t\right)-\mathrm{As}(t-\mathrm{\&tau;})-\mathrm{as}(t-\mathrm{\&epsiv;})|}^2\mathrm{dt}\}---\left(18\right)$

式中，θ＝(τ，ε，A，a)∈Θ是最大似然估计的四个参数，σ²是AWGN的平均功率。假设积分时间T是GNSS信号周期的整数倍。因此，最大似然函数可以表示为

$\begin{array}{c}\mathrm{\&Gamma;}\left(\mathrm{\&theta;}\right)=-{\&Integral;}_{t_0}^{t_0+T}{\left|r\left(t\right)\right|}^2\mathrm{dt}+2\mathrm{Re}\left\{A^{*}{\&Integral;}_{t_0}^{t_0+T}r\left(t\right)s^{*}(t-\mathrm{\&tau;})\mathrm{dt}\right\}-{\left|A\right|}^2{\&Integral;}_{t_0}^{t_0+T}{\left|s(t-\mathrm{\&tau;})\right|}^2\mathrm{dt}\\ +2\mathrm{Re}\left\{a^{*}{\&Integral;}_{t_0}^{t_0+T}r\left(t\right)s^{*}(t-\mathrm{\&epsiv;})\mathrm{dt}\right\}-{\left|a\right|}^2{\&Integral;}_{t_0}^{t_0+T}{\left|s(t-\mathrm{\&epsiv;})\right|}^2\mathrm{dt}-2\mathrm{Re}\left\{A{a}^{*}\right\}{\&Integral;}_{t_0}^{t_0+T}s(t-\mathrm{\&tau;})s^{*}(t-\mathrm{\&epsiv;})\mathrm{dt}\end{array}---\left(19\right)$

式中，Re{·}表示取实部运算。R_ss(ε-τ)表示s(t-τ)和s(t-ε)的自相关函数，它是一个仅取决于ε和τ的实数

$R_{\mathrm{ss}}(\mathrm{\&epsiv;}-\mathrm{\&tau;})={\&Integral;}_{t_0}^{t_0+T}s(t-\mathrm{\&tau;})s^{*}(t-\mathrm{\&epsiv;})\mathrm{dt}---\left(20\right)$

注意到r(t)，s(t-τ)和s(t-ε)的自相关函数均并不依赖于θ(τ，ε，A，a)，为了方便起见，用R表示s(t)的自相关函数。最大似然估计可以表示为

$\widehat{\mathrm{\&theta;}}(\mathrm{\&tau;},\mathrm{\&epsiv;},A,a)=\underset{\mathrm{\&theta;}\&Element;\mathrm{\&Theta;}}{\arg \max }\mathrm{\&Gamma;}\left(\mathrm{\&theta;}\right)---\left(21\right)$

估计参数可以通过令偏导数为零计算得到

$\frac{\&PartialD;\mathrm{\&Gamma;}\left(\mathrm{\&theta;}\right)}{\&PartialD;\mathrm{\&theta;}}=0---\left(22\right)$

分别对幅度相量A和a求导得

$\left\{\begin{array}{c}\frac{\&PartialD;\mathrm{\&Gamma;}\left(\mathrm{\&theta;}\right)}{\&PartialD;A}=R_{\mathrm{rs}}^{*}\left(\mathrm{\&tau;}\right)-A^{*}R-a^{*}{R}_{\mathrm{ss}}(\mathrm{\&tau;}-\mathrm{\&epsiv;})\\ \frac{\&PartialD;\mathrm{\&Gamma;}\left(\mathrm{\&theta;}\right)}{\&PartialD;a}=R_{\mathrm{rs}}^{*}\left(\mathrm{\&epsiv;}\right)-a^{*}R-A^{*}{R}_{\mathrm{ss}}(\mathrm{\&tau;}-\mathrm{\&epsiv;})\end{array}\right.---\left(23\right)$

令偏导数为零，可以得到幅度相量为

$\left\{\begin{array}{c}\hat{A}=\frac{{\mathrm{RR}}_{\mathrm{rs}}\left(\mathrm{\&tau;}\right)-R_{\mathrm{ss}}(\mathrm{\&epsiv;}-\mathrm{\&tau;})R_{\mathrm{rs}}\left(\mathrm{\&epsiv;}\right)}{R^2-R_{\mathrm{ss}}^2(\mathrm{\&epsiv;}-\mathrm{\&tau;})}\\ \hat{a}=\frac{{\mathrm{RR}}_{\mathrm{rs}}\left(\mathrm{\&epsiv;}\right)-R_{\mathrm{ss}}(\mathrm{\&epsiv;}-\mathrm{\&tau;})R_{\mathrm{rs}}\left(\mathrm{\&tau;}\right)}{R^2-R_{\mathrm{ss}}^2(\mathrm{\&epsiv;}-\mathrm{\&tau;})}\end{array}\right.---\left(24\right)$

同理，计算传播时延τ和ε的偏导数

$\left\{\begin{array}{c}\frac{\&PartialD;\mathrm{\&Gamma;}\left(\mathrm{\&theta;}\right)}{\&PartialD;\mathrm{\&tau;}}=2\mathrm{Re}\left\{A^{*}\frac{\&PartialD;}{\&PartialD;\mathrm{\&tau;}}\right[R_{\mathrm{rs}}\left(\mathrm{\&tau;}\right)-{\mathrm{aR}}_{\mathrm{ss}}(\mathrm{\&epsiv;}-\mathrm{\&tau;})\left]\right\}\\ \frac{\&PartialD;\mathrm{\&Gamma;}\left(\mathrm{\&theta;}\right)}{\&PartialD;\mathrm{\&epsiv;}}=2\mathrm{Re}\left\{a^{*}\frac{\&PartialD;}{\&PartialD;\mathrm{\&epsiv;}}\right[R_{\mathrm{rs}}\left(\mathrm{\&epsiv;}\right)-{\mathrm{AR}}_{\mathrm{ss}}(\mathrm{\&epsiv;}-\mathrm{\&tau;})\left]\right\}\end{array}\right.---\left(25\right)$

定义R_rA(τ)＝R_rs(τ)-aR_ss(ε-τ)，R_ra(τ)＝R_rs(τ)-AR_ss(ε-τ)。与幅度相量A和a不同的求解，传播时延估计没有解析表达式。因此，本发明采用牛顿迭代估计τ和ε

$\left\{\begin{array}{c}{\widehat{\mathrm{\&tau;}}}_{i+1}={\widehat{\mathrm{\&tau;}}}_i-\frac{\mathrm{Re}\left\{R_{\mathrm{rA}}^{*}\left({\widehat{\mathrm{\&tau;}}}_i\right)\frac{{\&PartialD;R}_{\mathrm{rA}}\left({\widehat{\mathrm{\&tau;}}}_i\right)}{\&PartialD;\mathrm{\&tau;}}\right\}}{\mathrm{Re}\left\{R_{\mathrm{rA}}^{*}\left({\widehat{\mathrm{\&tau;}}}_i\right)\frac{{\&PartialD;}^2{R}_{\mathrm{rA}}\left({\widehat{\mathrm{\&tau;}}}_i\right)}{{\&PartialD;\mathrm{\&tau;}}^2}\right\}+{\left|\frac{{\&PartialD;R}_{\mathrm{rA}}\left({\widehat{\mathrm{\&tau;}}}_i\right)}{\&PartialD;\mathrm{\&tau;}}\right|}^2}\\ {\widehat{\mathrm{\&epsiv;}}}_{i+1}={\widehat{\mathrm{\&epsiv;}}}_i-\frac{\mathrm{Re}\left\{R_{\mathrm{ra}}^{*}\left({\widehat{\mathrm{\&epsiv;}}}_i\right)\frac{{\&PartialD;R}_{\mathrm{ra}}\left({\widehat{\mathrm{\&epsiv;}}}_i\right)}{\&PartialD;\mathrm{\&epsiv;}}\right\}}{\mathrm{Re}\left\{R_{\mathrm{ra}}^{*}\left({\widehat{\mathrm{\&epsiv;}}}_i\right)\frac{{\&PartialD;}^2{R}_{\mathrm{ra}}\left({\widehat{\mathrm{\&epsiv;}}}_i\right)}{{\&PartialD;\mathrm{\&epsiv;}}^2}\right\}+{\left|\frac{{\&PartialD;R}_{\mathrm{ra}}\left({\widehat{\mathrm{\&epsiv;}}}_i\right)}{\&PartialD;\mathrm{\&epsiv;}}\right|}^2}\end{array}\right.---\left(26\right)$

式中，i表示迭代次数，分母第二项可以忽略。在实际中，我们采用早迟门方法计算一阶和二阶偏导数

$\frac{{\&PartialD;R}_{\mathrm{rA}}\left(\mathrm{\&tau;}\right)}{\&PartialD;\mathrm{\&tau;}}=\frac{R_{\mathrm{rA}}(\mathrm{\&tau;}+h)-R_{\mathrm{rA}}(\mathrm{\&tau;}-h)}{2h}---\left(27\right)$

$\frac{{\&PartialD;}^2{R}_{\mathrm{rA}}\left(\mathrm{\&tau;}\right)}{\&PartialD;{\mathrm{\&tau;}}^2}=\frac{R_{\mathrm{rA}}(\mathrm{\&tau;}+h)+R_{\mathrm{rA}}(\mathrm{\&tau;}-h)-2R_{\mathrm{rA}}\left(\mathrm{\&tau;}\right)}{h^2}---\left(28\right)$

式中，h是一个较小的正数。因此，(26)可以改写成

$\left\{\begin{array}{c}{\widehat{\mathrm{\&tau;}}}_{i+1}={\widehat{\mathrm{\&tau;}}}_i-\frac{\mathrm{hRe}\left\{R_{\mathrm{rA}}^{*}\left({\widehat{\mathrm{\&tau;}}}_i\right)\right[R_{\mathrm{rA}}({\widehat{\mathrm{\&tau;}}}_i+h)-R_{\mathrm{rA}}({\widehat{\mathrm{\&tau;}}}_i-h)\left]\right\}}{2\mathrm{Re}\left\{R_{\mathrm{rA}}^{*}\left({\widehat{\mathrm{\&tau;}}}_i\right)\right[R_{\mathrm{rA}}({\widehat{\mathrm{\&tau;}}}_i+h)+R_{\mathrm{rA}}({\widehat{\mathrm{\&tau;}}}_i-h)-2R_{\mathrm{rA}}\left({\widehat{\mathrm{\&tau;}}}_i\right)\left]\right\}}\\ {\widehat{\mathrm{\&epsiv;}}}_{i+1}={\widehat{\mathrm{\&epsiv;}}}_i-\frac{\mathrm{hRe}\left\{R_{\mathrm{ra}}^{*}\left({\widehat{\mathrm{\&epsiv;}}}_i\right)\right[R_{\mathrm{ra}}({\widehat{\mathrm{\&epsiv;}}}_i+h)-R_{\mathrm{ra}}({\widehat{\mathrm{\&epsiv;}}}_i-h)\left]\right\}}{2\mathrm{Re}\left\{R_{\mathrm{ra}}^{*}\left({\widehat{\mathrm{\&epsiv;}}}_i\right)\right[R_{\mathrm{ra}}({\widehat{\mathrm{\&epsiv;}}}_i+h){+R}_{\mathrm{ra}}({\widehat{\mathrm{\&epsiv;}}}_i-h)-2R_{\mathrm{ra}}\left({\widehat{\mathrm{\&epsiv;}}}_i\right)\left]\right\}}\end{array}\right.---\left(29\right)$

对于每次迭代，采用牛顿迭代法传播时延τ和ε沿着Γ(θ)的梯度方向更新。根据式(24)可以直接计算幅度相量A和a。经过多次迭代，估计值逐渐收敛到实际值。然而，由于BOC调制信号的自相关函数具有多峰特性，牛顿法容易收敛到局部最优解。为了找到全局最优解，本发明采用无模糊相关函数，R_rs(τ)通过对接收到的方波BOC信号r(t)和生成的本地正弦波信号s_I(t-τ)作相关运算，R_ss(ε-τ)通过对接收到的本地方波BOC信号s(t-τ)和生成的本地正弦波信号s_I(t-ε)作相关运算。

当信号存在非视距传播时，多径抑制技术可以区分多径分量。通过对接收序列建立随机模型，可以采用最大似然(maximum posteriori，ML)估计信号传播时延。特别是在先验信息已知时，接收机可以显著提高定位精度。牛顿迭代作为一种有效的多径抑制计算方法，定位精度可以逼近克拉美罗界(Cramer-Rao bound，CRB)。对于无偏估计量，克拉美罗界(Cramer-Rao bound，CRB)是估计精度的理论下界。对于单径传播，传播时延的CRB为

${\mathrm{CRB}}_1\left(\mathrm{\&tau;}\right)=\frac{{\mathrm{\&sigma;}}^2}{4{\mathrm{\&pi;}}^2{A}^2{\&Integral;}_{-B}^B{f}^{2}S\left(f\right)\mathrm{df}}---\left(30\right)$

在多径环境下，对于实值估计向量θ，估计C(θ)的协方差矩阵受限于

C(θ)≥J^-1(θ)   (31)

式中，J(θ)表示Fisher信息矩阵(Fisher information matrix，FIM)，其逆矩阵是CRB矩阵。为了保证估计向量的所有元素均为示数，最大似然估计参数向量可以改写为

$J\left(\mathrm{\&theta;}\right)=-E\left[\frac{{\&PartialD;}^2\ln p(r;\mathrm{\&theta;})}{\&PartialD;\mathrm{\&theta;}\&PartialD;{\mathrm{\&theta;}}^T}\right]---\left(32\right)$

式中，E(·)表示数学期望。矩阵J(θ)中的每个元素可以表示为

$J{\left(\mathrm{\&theta;}\right)}_{i,j}=-E\left[\frac{{\&PartialD;}^2\ln p(r;\mathrm{\&theta;})}{{\&PartialD;\mathrm{\&theta;}}_i{\&PartialD;\mathrm{\&theta;}}_j}\right]---\left(33\right)$

首先对Fisher信息矩阵进行分块

$J\left(\mathrm{\&theta;}\right)=\left[\begin{array}{cc}{A}_{2\&times;2}& B_{2\&times;4}\\ C_{4\&times;2}& D_{4\&times;4}\end{array}\right]---\left(34\right)$

若D_4×4可逆，则J^-1(θ)可逆等价于(A-BD^-1C)可逆。基于此，分块矩阵的逆矩阵可以表示为

$J^{-1}\left(\mathrm{\&theta;}\right)=\left[\begin{array}{cc}{(A-{\mathrm{BD}}^{-1}C)}^{-1}& -(A-{\mathrm{BD}}^{-1}C){\mathrm{BD}}^{-1}\\ {-D}^{-1}C{(A-{\mathrm{BD}}^{-1}C)}^{-1}& D^{-1}+D^{-1}C{(A-{\mathrm{BD}}^{-1}C)}^{-1}{\mathrm{BD}}^{-1}\end{array}\right]---\left(35\right)$

式中，J^-1(θ)_1，1代表多径环境下传播时延的CRB，需要分别计算矩阵A_2×2，B_2×4，C_4×2，D_4×4。F_Gabor表示带限信号的均方根带宽。根据帕斯瓦尔定理，A_2×2可以表示为

$A_{2\&times;2}=\frac{1}{{\mathrm{\&sigma;}}^2}\left[\begin{array}{cc}{4\mathrm{\&pi;}}^2{\left|A\right|}^2{f}_{\mathrm{Gabor}}^2& -\mathrm{Re}\left\{{\mathrm{Aa}}^{*}\right\}{R}^{\&prime;\&prime;}\left(\mathrm{\&delta;}\right)\\ -\mathrm{Re}\left\{{\mathrm{Aa}}^{*}\right\}{R}^{\&prime;\&prime;}\left(\mathrm{\&delta;}\right)& {4\mathrm{\&pi;}}^2{\left|a\right|}^2{f}_{\mathrm{Gabor}}^2\end{array}\right]---\left(36\right)$

式中，δ＝ε-τ，R(δ)是本地方波s(t-τ)和生成的本地正弦波子载波s_I(t-ε)的相关函数

$R\left(\mathrm{\&delta;}\right)={\&Integral;}_{-B}^{B}F\left(f\right)e^{j2\mathrm{\&pi;f\&delta;}}\mathrm{df}---\left(37\right)$

式中，F(f)是相关函数R(δ)的傅里叶变换。R″(δ)是R(δ)对于δ的二阶偏导数

$R^{\&prime;\&prime;}\left(\mathrm{\&delta;}\right)=-4{\mathrm{\&pi;}}^2{\&Integral;}_{-B}^B{f}^{2}F\left(f\right)e^{j2\mathrm{\&pi;f\&delta;}}\mathrm{df}---\left(38\right)$

矩阵C_4×2是矩阵B_2×4的转置。因此，B_2×4和C_4×2可以表示为

$C_{4\&times;2}=B_{4\&times;2}^T=\frac{1}{{\mathrm{\&sigma;}}^2}\left[\begin{array}{cc}0& \mathrm{Re}\left\{a\right\}{R}^{\&prime;}\left(\mathrm{\&delta;}\right)\\ -\mathrm{Re}\left\{A\right\}{R}^{\&prime;}\left(\mathrm{\&delta;}\right)& 0\\ 0& \mathrm{Im}\left\{a\right\}{R}^{\&prime;}\left(\mathrm{\&delta;}\right)\\ -\mathrm{Im}\left\{A\right\}{R}^{\&prime;}\left(\mathrm{\&delta;}\right)& 0\end{array}\right]---\left(39\right)$

式中，R′(δ)是R(δ)对于δ的一阶偏导数

$R^{\&prime;}\left(\mathrm{\&delta;}\right)=-j2\mathrm{\&pi;}{\&Integral;}_{-B}^B\mathrm{fF}\left(f\right)e^{j2\mathrm{\&pi;f\&delta;}}\mathrm{df}---\left(40\right)$

R表示s(t)和生成的本地正弦亚载波s_I(t)的相关函数。因此，D_4×4可以表示为

$D_{4\&times;4}=\frac{1}{{\mathrm{\&sigma;}}^2}\left[\begin{array}{cccc}R& R\left(\mathrm{\&delta;}\right)& 0& 0\\ R\left(\mathrm{\&delta;}\right)& R& 0& 0\\ 0& 0& R& R\left(\mathrm{\&delta;}\right)\\ 0& 0& R\left(\mathrm{\&delta;}\right)& R\end{array}\right]---\left(41\right)$

综上所述，多径环境下传播时延的CRB可以表示为

${\mathrm{CRB}}_2\left(\mathrm{\&tau;}\right)=\frac{{(A-{\mathrm{BD}}^{-1}C)}_{\mathrm{2,2}}}{{(A-{\mathrm{BD}}^{-1}C)}_{\mathrm{1,1}}{(A-{\mathrm{BD}}^{-1}C)}_{\mathrm{2,2}}-{(A-{\mathrm{BD}}^{-1}C)}_{\mathrm{1,2}}^2}---\left(42\right)$

式中，

${(A-{\mathrm{BD}}^{-1}C)}_{\mathrm{1,1}}=4{\mathrm{\&pi;}}^2{\left|A\right|}^2{f}_{\mathrm{Gabor}}^2-\frac{R{\left|A\right|}^2{R}^{\&prime;\&prime;2}\left(\mathrm{\&delta;}\right)}{R^2-R^2\left(\mathrm{\&delta;}\right)}---\left(43\right)$

${(A-{\mathrm{BD}}^{-1}C)}_{\mathrm{1,2}}=-\mathrm{Re}\left\{{\mathrm{Aa}}^{*}\right\}[R^{\&prime;\&prime;}\left(\mathrm{\&delta;}\right)+\frac{R\left(\mathrm{\&Delta;}\right)R^{\&prime;2}\left(\mathrm{\&delta;}\right)}{R^2-R^2\left(\mathrm{\&delta;}\right)}]---\left(44\right)$

${(A-{\mathrm{BD}}^{-1}C)}_{\mathrm{2,2}}=4{\mathrm{\&pi;}}^2{\left|a\right|}^2{f}_{\mathrm{Gabor}}^2-\frac{R{\left|a\right|}^2{R}^{\&prime;\&prime;2}\left(\mathrm{\&delta;}\right)}{R^2-R^2\left(\mathrm{\&delta;}\right)}---\left(45\right)$

当接收信号不存在视距分量时，接收信号可以看作最先到达的多径分量和噪声的叠加。此时，多径环境下概率密度函数为

$p(r;\mathrm{\&theta;})=\frac{1}{\sqrt{2\mathrm{\&pi;}}\mathrm{\&sigma;}}\exp \{-\frac{1}{2{\mathrm{\&sigma;}}^2}{\&Integral;}_{t_0}^{t_0+T}{|r\left(t\right)-\mathrm{as}(t-\mathrm{\&epsiv;})|}^2\mathrm{dt}\}---\left(46\right)$

同理可得多径分量的幅度为

$\hat{a}=\frac{R_{\mathrm{rs}}\left(\mathrm{\&epsiv;}\right)}{R}---\left(47\right)$

多径时延的迭代公式为

${\widehat{\mathrm{\&epsiv;}}}_{i+1}={\widehat{\mathrm{\&epsiv;}}}_i-\frac{\mathrm{Re}\left\{R_{\mathrm{rs}}^{*}\left({\widehat{\mathrm{\&epsiv;}}}_i\right)\frac{\&PartialD;R_{\mathrm{rs}}\left({\widehat{\mathrm{\&epsiv;}}}_i\right)}{\&PartialD;\mathrm{\&epsiv;}}\right\}}{\mathrm{Re}\left\{R_{\mathrm{rs}}^{*}\left({\widehat{\mathrm{\&epsiv;}}}_i\right)\frac{{\&PartialD;}^2{R}_{\mathrm{ra}}\left({\widehat{\mathrm{\&epsiv;}}}_i\right)}{\&PartialD;{\mathrm{\&epsiv;}}^2}\right\}+{\left|\frac{\&PartialD;R_{\mathrm{rs}}\left({\widehat{\mathrm{\&epsiv;}}}_i\right)}{\&PartialD;\mathrm{\&epsiv;}}\right|}^2}---\left(48\right)$

根据早迟门同步法，公式(48)可以改写为

${\widehat{\mathrm{\&epsiv;}}}_{i+1}={\widehat{\mathrm{\&epsiv;}}}_i-\frac{\mathrm{hRe}\left\{R_{\mathrm{rs}}^{*}\left({\widehat{\mathrm{\&epsiv;}}}_i\right)\right[R_{\mathrm{rs}}({\widehat{\mathrm{\&epsiv;}}}_i+h)-R_{\mathrm{rs}}({\widehat{\mathrm{\&epsiv;}}}_i-h)\left]\right\}}{2\mathrm{Re}\left\{R_{\mathrm{rs}}^{*}\left({\widehat{\mathrm{\&epsiv;}}}_i\right)\right[R_{\mathrm{rs}}({\widehat{\mathrm{\&epsiv;}}}_i+h)+R_{\mathrm{rs}}({\widehat{\mathrm{\&epsiv;}}}_i-h)-2R_{\mathrm{rs}}\left({\widehat{\mathrm{\&epsiv;}}}_i\right)\left]\right\}}---\left(49\right)$

由于接收信号不存在视距分量，运用上述方法直接估计视距传播时间会产生很大误差。大量统计规律表明，多径分量到达时间间隔服从指数分布，将最先到达的多径时延ε₁，ε₂，...，ε_n作为样本，采用最大似然函数

$L({\mathrm{\&epsiv;}}_1,{\mathrm{\&epsiv;}}_2,\&CenterDot;\&CenterDot;\&CenterDot;,{\mathrm{\&epsiv;}}_n;\widehat{\mathrm{\&tau;}})=\underset{\mathrm{\&tau;}\&Element;T}{\max }\underset{i=1}{\overset{n}{\mathrm{\&Pi;}}}f({\mathrm{\&epsiv;}}_i;\mathrm{\&tau;})---\left(50\right)$

式中，f(ε₁；τ)表示指数分布函数。利用多径时延得到τ的最大似然估计量为

τ＝min{ε₁，ε₂，…，ε_n}   (51)

如上所述，尽管参照特定的优选实施例已经表示和表述了本发明，但其不得解释为对本发明自身的限制。在不脱离所附权利要求定义的本发明的精神和范围前提下，可对其在形式上和细节上作出各种变化。

Claims (6)

1.多径环境下全球卫星导航系统接收机无模糊处理方法，其特征在于：包括以下步骤：

(一)根据BOC调制信号的统计特性，设计旁瓣抑制相关器，使得接收机跟踪到主瓣而不是旁瓣、输出具有单峰特性，并且相关函数随时延的增加单调递减；

(二)建立多径环境下GNSS信号传播模型，定义接收到的复GNSS信号和本地载波的无模糊相关函数|R_rs(τ)|；

(三)接收机根据相关器输出判断接收信号是否存在视距信号，若存在视距信号，接收机采用最大似然估计得到实时位置信息，运用牛顿迭代法更新幅度、时延、相位参数，若不存在视距信号，则利用最先到达多径分量的统计特性修正理想状态下视距传播时延；

(四)设计多径环境下GNSS无模糊处理接收机，采用牛顿迭代计算多径时延，定位精度逼近克拉美罗界CRB。

2.根据权利要求1所述的多径环境下全球卫星导航系统接收机无模糊处理方法，其特征在于：所述步骤(一)中设计旁瓣抑制相关器采用波形为正弦信号的亚载波代替传统的方波亚载波，包括如下步骤：

(1)将接收到的BOC调制信号r(t)与本地同相正弦亚载波信号s_I(t)进行相关运算得到R_I(τ)；

(2)将接收到的BOC调制信号r(t)与本地正交余弦亚载波信号s_Q(t)进行相关运算得到R_Q(τ)；

(3)将同相相关函数R_I(τ)和正交相关函数R_Q(τ)进行平方和运算得到无模糊相关函数R(τ)。

3.根据权利要求1所述的多径环境下全球卫星导航系统接收机无模糊处理方法，其特征在于：所述步骤(二)多径环境下解调后接收复GNSS信号可以表示为

式中，A(t)表示视距分量的幅度，τ表示视距信号传播时间，表示载波相位偏移，P表示多径信号的数目，a_p(t)表示第p条多径分量的幅度，τ_p表示第p条多径分量的时延，表示第p条多径分量的相位偏移，n(t)表示复噪声。

假设信号的幅度和传播时延在单位观测时间内保持不变，在只有视距传播信号的情况下，传输过程可以看作加性高斯白噪声AWGN卫星信道，在存在多径反射信号的情况下，接收信号可以简化为视距分量和最先到达的非视距分量，根据中心极限定理，密集多径条件下接收信号模型可以简化为单路径传播，即接收到的复GNSS信号可以简化为

r(t)＝As(t-τ)+as(t-ε)+n(t)

式中，A表示视距相量，τ信号传播时间，a表示最先到达的非视距幅度相量，ε表示最先到达的非视距幅度分量；

定义接收到的复GNSS信号和本地载波的无模糊相关函数为

$R_{\mathrm{rs}}\left(\mathrm{\&tau;}\right)={\&Integral;}_{t_0}^{t_0+T}r\left(t\right)s^{*}(t-\mathrm{\&tau;})\mathrm{dt}.$

4.根据权利要求1所述的多径环境下全球卫星导航系统接收机无模糊处理方法，其特征在于，所述步骤(三)包括如下步骤：

若无模糊相关函数|R_rs(τ)|大于或等于门限δ，则存在视距传播信号，此时接收信号r(t)＝As(t-τ)+as(t-ε)+n(t)，接收机采用最大似然估计得到实时位置信息，运用牛顿迭代法更新视距分量和最先到达的多径分量的幅度、时延、相位；

若无模糊相关函数|R_rs(τ)|小于门限δ，则不存在视距传播信号，此时接收信号r(t)＝as(t-ε)+n(t)，接收机首先估计得到最先到达的多径分量的幅度、时延、相位等参数，根据多径分量到达时间间隔服从指数分布的特性，采用最大似然估计修正理想状态下视距分量传播时延。

5.根据权利要求1所述的多径环境下全球卫星导航系统接收机无模糊处理方法，其特征在于，所述步骤(四)包括如下步骤：

多径环境下传播时延的克拉美罗界CRB表示为

${\mathrm{CRB}}_2\left(\mathrm{\&tau;}\right)\frac{{(A-{\mathrm{BD}}^{-1}C)}_{\mathrm{2,2}}}{{(A-{\mathrm{BD}}^{-1}C)}_{\mathrm{1,1}}{(A-{\mathrm{BD}}^{-1}C)}_{\mathrm{2,2}}-{(A-{\mathrm{BD}}^{-1}C)}_{\mathrm{1,2}}^2}$

式中，

${(A-{\mathrm{BD}}^{-1}C)}_{\mathrm{1,1}}=4{\mathrm{\&pi;}}^2{\left|A\right|}^2{f}_{\mathrm{Gabor}}^2-\frac{R{\left|A\right|}^2{R}^{\&prime;\&prime;2}\left(\mathrm{\&delta;}\right)}{R^2-R^2\left(\mathrm{\&delta;}\right)}$

${(A-{\mathrm{BD}}^{-1}C)}_{\mathrm{1,2}}=-\mathrm{Re}\left\{{\mathrm{Aa}}^{*}\right\}[R^{\&prime;\&prime;}\left(\mathrm{\&delta;}\right)+\frac{R\left(\mathrm{\&Delta;}\right)R^{\&prime;2}\left(\mathrm{\&delta;}\right)}{R^2-R^2\left(\mathrm{\&delta;}\right)}]$

${(A-{\mathrm{BD}}^{-1}C)}_{\mathrm{2,2}}=4{\mathrm{\&pi;}}^2{\left|a\right|}^2{f}_{\mathrm{Gabor}}^2-\frac{R{\left|a\right|}^2{R}^{\&prime;\&prime;2}\left(\mathrm{\&delta;}\right)}{R^2-R^2\left(\mathrm{\&delta;}\right)}$

当接收信号不存在视距分量时，接收信号可以看作最先到达的多径分量和噪声的叠加，此时，多径环境下概率密度函数为

$p(r;\mathrm{\&theta;})=\frac{1}{\sqrt{2\mathrm{\&pi;}}\mathrm{\&sigma;}}\exp \{-\frac{1}{{2\mathrm{\&sigma;}}^2}{\&Integral;}_{t_0}^{t_0+T}{|r\left(t\right)-\mathrm{as}(t-\mathrm{\&epsiv;})|}^2\mathrm{dt}\}$

同理可得多径分量的幅度为

$\hat{a}=\frac{R_{\mathrm{rs}}\left(\mathrm{\&epsiv;}\right)}{R}$

多径时延的迭代公式为

${\widehat{\mathrm{\&epsiv;}}}_{i+1}={\widehat{\mathrm{\&epsiv;}}}_i\frac{\mathrm{Re}\left\{R_{\mathrm{rs}}^{*}\left({\widehat{\mathrm{\&epsiv;}}}_i\right)\frac{{\&PartialD;R}_{\mathrm{rs}}\left({\widehat{\mathrm{\&epsiv;}}}_i\right)}{{\&PartialD;}_3}\right\}}{R_e\left\{R_{\mathrm{rs}}^{*}\left({\widehat{\mathrm{\&epsiv;}}}_i\right)\frac{{\&PartialD;}^2{R}_{\mathrm{ra}}\left({\widehat{\mathrm{\&epsiv;}}}_i\right)}{{\&PartialD;\mathrm{\&epsiv;}}^2}\right\}+{\left|\frac{{\&PartialD;R}_{\mathrm{rs}}\left({\widehat{\mathrm{\&epsiv;}}}_i\right)}{\&PartialD;\mathrm{\&epsiv;}}\right|}^2}$

根据早迟门同步法，上式可以改写为

${\widehat{\mathrm{\&epsiv;}}}_{i+1}={\widehat{\mathrm{\&epsiv;}}}_i-\frac{\mathrm{hRe}\left\{R_{\mathrm{rs}}^{*}\left({\widehat{\mathrm{\&epsiv;}}}_i\right)\right[R_{\mathrm{rs}}({\widehat{\mathrm{\&epsiv;}}}_i+h)-R_{\mathrm{rs}}({\widehat{\mathrm{\&epsiv;}}}_i-h)\left]\right\}}{2\mathrm{Re}\left\{R_{\mathrm{rs}}^{*}\left({\widehat{\mathrm{\&epsiv;}}}_i\right)\right[R_{\mathrm{rs}}({\widehat{\mathrm{\&epsiv;}}}_i+h)+R_{\mathrm{rs}}({\widehat{\mathrm{\&epsiv;}}}_i-h)-2R_{\mathrm{rs}}\left({\widehat{\mathrm{\&epsiv;}}}_i\right)\text{]}}}$

多径分量到达时间间隔服从指数分布，将最先到达的多径时延ε₁，ε₂，...，ε_n作为样本，采用最大似然函数

$L({\mathrm{\&epsiv;}}_1,{\mathrm{\&epsiv;}}_2,\&CenterDot;\&CenterDot;\&CenterDot;,{\mathrm{\&epsiv;}}_n;\widehat{\mathrm{\&tau;}})=\underset{\mathrm{\&tau;}\&Element;T}{\max }\underset{i=1}{\overset{n}{\mathrm{\&Pi;}}}f({\mathrm{\&epsiv;}}_i;\mathrm{\&tau;})$

式中，f(ε₁；τ)表示指数分布函数，利用多径时延得到τ的最大似然估计量为

τ＝min{ε₁，ε₂，…，ε_n}。

6.根据权利要求4所述的多径环境下全球卫星导航系统接收机无模糊处理方法，其特征在于：

当接收信号存在视距分量时，多径环境下接收信号的概率密度函数为

$p(r;\mathrm{\&theta;})=\frac{1}{\sqrt{2\mathrm{\&pi;}}\mathrm{\&sigma;}}\exp \{-\frac{1}{2{\mathrm{\&sigma;}}^2}{\&Integral;}_{t_0}^{t_0+T}|r\left(t\right)-\mathrm{As}(t-\mathrm{\&tau;})-\mathrm{as}(t-\mathrm{\&epsiv;}){|}^2\mathrm{dt}\}$

式中，θ＝(τ，ε，A，a)∈Θ是最大似然估计的四个参数，σ²是AWGN的平均功率。假设积分时间T是GNSS信号周期的整数倍，最大似然函数可以表示为

$\begin{array}{c}\mathrm{\&Gamma;}\left(\mathrm{\&theta;}\right)=-{\&Integral;}_{t_0}^{t_0+T}{\left|r\left(t\right)\right|}^2\mathrm{dt}+2\mathrm{Re}\left\{A^{*}{\&Integral;}_{t_0}^{t_0+T}r\left(t\right)s^{*}(t-\mathrm{\&tau;})\mathrm{dt}\right\}-{\left|A\right|}^2{\&Integral;}_{t_0}^{t_0+T}{\left|s(t-\mathrm{\&tau;})\right|}^2\mathrm{dt}\\ +2\mathrm{Re}\left\{a^{*}{\&Integral;}_{t_0}^{t_0+T}r\left(t\right)s^{*}(t-\mathrm{\&epsiv;})\mathrm{dt}\right\}-{\left|a\right|}^2{\&Integral;}_{t_0}^{t_0+T}{\left|s(t-\mathrm{\&epsiv;})\right|}^2\mathrm{dt}-2\mathrm{Re}\left\{A{a}^{*}\right\}{\&Integral;}_{t_0}^{t_0+T}s(t-\mathrm{\&tau;})s^{*}(t-\mathrm{\&epsiv;})\mathrm{dt}\end{array}$

式中，R_ss(ε-τ)表示s(t-τ)和s(t-ε)的自相关函数，它是一个仅取决于ε和τ的实数

$R_{\mathrm{ss}}(\mathrm{\&epsiv;}-\mathrm{\&tau;})={\&Integral;}_{t_0}^{t_0+T}s(t-\mathrm{\&tau;})s^{*}(t-\mathrm{\&epsiv;})\mathrm{dt}$

求偏导，幅度相量为

$\left\{\begin{array}{c}\hat{A}=\frac{{\mathrm{RR}}_{\mathrm{rs}}\left(\mathrm{\&tau;}\right)-R_{\mathrm{ss}}(\mathrm{\&epsiv;}-\mathrm{\&tau;})R_{\mathrm{rs}}\left(\mathrm{\&epsiv;}\right)}{R^2-R_{\mathrm{ss}}^2(\mathrm{\&epsiv;}-\mathrm{\&tau;})}\\ \hat{a}=\frac{{\mathrm{RR}}_{\mathrm{rs}}\left(\mathrm{\&epsiv;}\right)-R_{\mathrm{ss}}(\mathrm{\&epsiv;}-\mathrm{\&tau;})R_{\mathrm{rs}}\left(\mathrm{\&tau;}\right)}{R^2-R_{\mathrm{ss}}^2(\mathrm{\&epsiv;}-\mathrm{\&tau;})}\end{array}\right.$

采用牛顿迭代估计传播时延τ和ε

$\left\{\begin{array}{c}{\widehat{\mathrm{\&tau;}}}_{i+1}={\widehat{\mathrm{\&tau;}}}_i-\frac{\mathrm{h\; Re}\left\{R_{\mathrm{rA}}^{*}\left({\widehat{\mathrm{\&tau;}}}_i\right)\right[R_{\mathrm{rA}}({\widehat{\mathrm{\&tau;}}}_i+h)-R_{\mathrm{rA}}({\widehat{\mathrm{\&tau;}}}_i-h)\left]\right\}}{2\mathrm{Re}\left\{R_{\mathrm{rA}}^{*}\left({\widehat{\mathrm{\&tau;}}}_i\right[R_{\mathrm{rA}}({\widehat{\mathrm{\&tau;}}}_i+h)+R_{\mathrm{rA}}({\widehat{\mathrm{\&tau;}}}_i-h)-2R_{\mathrm{rA}}\left({\widehat{\mathrm{\&tau;}}}_i\right)\left]\right)\right\}}\\ {\widehat{\mathrm{\&epsiv;}}}_{i+1}={\widehat{\mathrm{\&epsiv;}}}_i-\frac{\mathrm{h\; Re}\left\{R_{\mathrm{ra}}^{*}\left({\widehat{\mathrm{\&epsiv;}}}_i\right)\right[R_{\mathrm{ra}}({\widehat{\mathrm{\&epsiv;}}}_i+h)-R_{\mathrm{ra}}({\widehat{\mathrm{\&epsiv;}}}_i-h)\left]\right\}}{2\mathrm{Re}\left\{R_{\mathrm{ra}}^{*}\left({\widehat{\mathrm{\&epsiv;}}}_i\right)\right[R_{\mathrm{ra}}({\widehat{\mathrm{\&epsiv;}}}_i+h)+R_{\mathrm{ra}}({\widehat{\mathrm{\&epsiv;}}}_i-h)-{2R}_{\mathrm{ra}}\left({\widehat{\mathrm{\&epsiv;}}}_i\right)\left]{}_{}\right\}}\end{array}\right..$