CN101655546B - 卫星导航直接信号参数与多径信号数目和参数的估计方法 - Google Patents

Abstract

本发明涉及一种卫星导航直接信号参数与多径信号数目和参数的估计方法，技术特征在于：采用Dirichlet过程DP(α，G₀)作为多径信号数目M的先验分布，则M的分布符合Dirichlet过程：

根据此先验分布采样得到多径信号可能的数目；根据多径信号可能的数目，构建以直接信号和多径信号的参数为状态量的状态方程x(k+1)＝Fx(k)+w(k)和观测向量的观测方程y＝aR(τ)+n；根据多径信号可能的数目和状态方程和观测方程，采用粒子滤波方法估计出直接信号和多径信号的参数，同时可以确定多径的数目。本发明可以在不了解多径信号数目的情况下估计出多径信号的数目和参数，不但适用于直接信号大于多径信号的情况，还适用于直接信号小于多径信号的情况。

Description

卫星导航直接信号参数与多径信号数目和参数的估计方法

技术领域

本发明涉及一种卫星导航直接信号参数与多径信号数目和参数的估计方法，属于卫星导航信号处理方法。 

背景技术

卫星导航多径信号估计方法通过将多径信号估计出来以消除卫星导航***(包括GPS，GLONASS，GALILEO，北斗***等)多径误差对定位结果的影响。 

对卫星导航多径信号的估计方法的研究分为基于最大似然的估计方法和基于贝叶斯估计的方法。基于最大似然估计的方法通过最大化似然函数来找到参数的最优估计值，对于这种方法，可参考“K.Kaindl，N.Niklasch，Improving Path Separation in a Multipath Environment by Applying Nonlinear Estimation Theory[C]，ION GPS2002，Portland，OR，24-27September 2002”。在这类技术中，多径估计延迟锁定环(Multipath Estimating Delay Lock Loop，MEDLL)技术是其典型代表，这种技术将估计器引入跟踪环路，可以实现对直接信号和多径信号的估计。 

MEDLL技术出现较早，其精度很高，但效率不高，为了提高计算效率，又出现了MMT(multipath mitigation technology)技术，见“L.R.Weill，Multipath Mitigation Using Modernized GPS Signals：How Good Can it Get？[C]ION GPS 2002，Portland，OR，24-27 September 2002”，MMT采用非线性变换将参数空间的维数减少，提高计算效率。MMT技术在实际中已经获得应用，NovAtel公司的VC(vision correlator)技术就是以MMT为基础的，VC技术见“M.Sahmoudi and R.Jr.Landry，Multipath Mitigation Techniques Using Maximum-Likelihood Principle[J]，Inside GNSSnovember～december 2008”。 

基于最大似然的多径信号估计技术的精度虽然很高，但其缺点是只能针对静态情况，这限制了这一类估计技术的应用范围。 

多径信号的贝叶斯估计技术采用状态空间模型，将运动信息引入估计过程，可以克服最大似然估计的缺点，适用于动态情况的估计。文献“G.I.Jee，H.S.Kim，Y.J.Lee，C.G.Park，A GPS C/A Code Tracking Loop Based on Extended Kalman Filter with Multipath Mitigation[C]，ION GPS 2002”提出了基于扩展卡尔曼滤波的码跟踪环路来估计多径信号，直接信号与多径信号的幅度、延迟和相位都可由扩展卡尔曼滤波估计出。文献“P.Closas，C.Fernandez-Prades，A.Juan，Bayesian DLL for Multipath Mitigation In Navigation Systems Using Particle Filters[C]，IEEE International Conference on Acoustics，Speech，and Signal Processing，2006”，“M.Lentmaier，B.Krach and P.Robertson，Bayesian Time Delay Estimation of GNSS Signals in Dynamic Multipath Environments[J]，International Journal of Navigation and Observation Volume 2008”，“P.Closas，C.Fernandez-Prades and J.A.Fernandez-Rubio，A Particle Filtering Tracking Algorithm for GNSSSynchronization Using Laplace’S Method[C]，ICASSP’081347，2008”都是采用基于粒子滤波估计器的技术对多径信号进行估计。其中，文献“P.Closas，C.Fernandez-Prades，A.Juan，Bayesian DLL for Multipath Mitigation In Navigation Systems Using Particle Filters[C]，IEEE International Conference on Acoustics，Speech，and Signal Processing，2006”采用最小二乘技术对信号幅度进行估计，减少了状态维数，提高了计算效率。文献“P.Closas，C.Fernandez-Prades and J.A.Fernandez-Rubio，A Particle Filtering Tracking Algorithm for GNSS Synchronization Using Laplace’S Method[C]，ICASSP’08 1347，2008”采用 Rao-Blackwellization粒子滤波来减少状态维数，重要性函数的选择采用Laplace方法，效率有一定的提高。 

各种已有的多径信号估计方法都有一个共同的特点，就是在应用它们估计多径信号之前，必须首先对多径信号的数目加以假设，然后才能构造相应的模型进行多径信号参数的估计，如果事先不了解多径信号的数目，则无法确定模型的结构，当对多径信号数目的假设与实际不符时，误差会很大，估计的精度就会大大降低。 

在实际环境中，比如市区这种多径影响大的环境，当携带接收机的载体运动时，来自建筑物和车辆等反射体的多径反射在不断变化，因此很难确定在某一时刻多径信号的数目，在这种情况下，应用已有方法就会导致很大的误差。 

发明内容

要解决的技术问题 

为了避免现有技术的不足之处，本发明提出一种卫星导航直接信号参数与多径信号数目和参数的估计方法，可以在不需要多径信号数目假设的条件下估计出多径信号的数目和多径信号的参数，克服了已有方法的缺陷，能够在多径变化大的条件下更好的估计出多径。 

本发明的思想在于：卫星导航直接信号指从导航卫星直接传送到导航卫星接收机的信号，多径信号(又称为多路径信号)指导航卫星发射的信号经过反射、散射过程间接传送到导航卫星接收机的信号。直接信号只有一个，而多径信号却可能有很多个。本发明提出的方法可以估计出直接信号的参数、多径信号的数目和多径信号的参数。采用Dirichlet过程这种概率模型对多径信号的数目建模，可以很好的根据实际情况估计出多径信号的数目，然后根据多径信号的数目构造信号参数的估计模型，以此对直接信号和多径信号的参数进行估计。 

技术方案 

一种卫星导航直接信号参数与多径信号数目和参数的估计方法，其特征在于步骤如下： 

步骤1：采用Dirichlet过程DP(α，G₀)作为多径信号数目M的先验分布，则M的分布符合Dirichlet过程： $\begin{array}{c}G~\mathrm{DP}(\mathrm{\α},G_0)\\ M\left(k\right)~G\end{array},$ 根据此先验分布采样得到多径信号可能的数目M，其中：α为参数，G₀为参数为p的几何分布Geo(p)，p的参考范围0＜p＜1； 

步骤2：根据多径信号可能的数目，构建以直接信号和多径信号的参数为状态量的状态方程x(k+1)＝Fx(k)+w(k)和观测向量的观测方程y＝aR(τ)+n； 

其中：k为离散时间变量； 

F＝diag{F_a，F_τ}： $F_a=\left[\begin{array}{cccc}{\mathrm{\λ}}_0& 0& ...& 0\\ 0& {\mathrm{\λ}}_1& ...& 0\\ ...& ...& ...& ...\\ 0& 0& ...& {\mathrm{\λ}}_M\end{array}\right],$ $F_{\mathrm{\τ}}=\left[\begin{array}{cccc}{\mathrm{\γ}}_0& 0& ...& 0\\ 0& {\mathrm{\γ}}_1& ...& 0\\ ...& ...& ...& ...\\ 0& 0& ...& {\mathrm{\γ}}_M\end{array}\right],$ 此两矩阵对角线元素参考取值为0.9～1； 

w(k)为高斯白噪声：其协方差阵为∑_k，x＝diag{∑_k，a，∑_k，τ}，∑_k，a，∑_k，τ为对角阵，对角线元素参考取值为0.0001～0.01； 

y为观测向量，其中的一个分量为： 

$y_{{\mathrm{\δ}}_j}\left(t_k\right)={\mathrm{\α}}_0\left(t_k\right)R(\mathrm{\τ}-{\mathrm{\τ}}_0+{\mathrm{\δ}}_j)e^{j{\mathrm{\Δ\φ}}_0}+\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\α}}_i\left(t_k\right)R(\mathrm{\τ}-{\mathrm{\τ}}_i+{\mathrm{\δ}}_j)e^{j{\mathrm{\Δ\φ}}_i}+n\left(t_k\right),$ N维观测向量为 

δ_j为产生第j个观测分量的卫星导航接收机延迟锁定环相关累加支路的延迟； 

a为复向量， $a=[a_0,a_1,...,a_M]=[{\mathrm{\α}}_0{e}^{j{\mathrm{\Δ\φ}}_0},{\mathrm{\α}}_1{e}^{j{\mathrm{\Δ\φ}}_1},...,{\mathrm{\α}}_M{e}^{j{\mathrm{\Δ\φ}}_M}],$ 其中Δφ₁，Δφ₂，…Δφ_M为相位差； 

R(τ)为包含直接信号和多径信号延迟的相关函数矩阵： 

$R\left(\mathrm{\τ}\right)=\left[\begin{array}{cccc}R(\mathrm{\τ}-{\mathrm{\τ}}_0+{\mathrm{\δ}}_1)& R(\mathrm{\τ}-{\mathrm{\τ}}_0+{\mathrm{\δ}}_2)& ...& R(\mathrm{\τ}-{\mathrm{\τ}}_0+{\mathrm{\δ}}_N)\\ R(\mathrm{\τ}-{\mathrm{\τ}}_1+{\mathrm{\δ}}_1)& R(\mathrm{\τ}-{\mathrm{\τ}}_1+{\mathrm{\δ}}_2)& ...& R(\mathrm{\τ}-{\mathrm{\τ}}_1+{\mathrm{\δ}}_N)\\ ...& ...& ...& ...\\ R(\mathrm{\τ}-{\mathrm{\τ}}_M+{\mathrm{\δ}}_1)& R(\mathrm{\τ}-{\mathrm{\τ}}_M+{\mathrm{\δ}}_2)& ...& R(\mathrm{\τ}-{\mathrm{\τ}}_M+{\mathrm{\δ}}_N)\end{array}\right],$ 这里R(·)为卫星信号和接收机自身参考信号的相关函数，n为噪声，大小由信噪比确定； 

所述的直接信号和多径信号的参数为直接信号和多径信号的幅度向量和延迟向量；其中 

所述幅度向量为 所述延迟向量为 

其中：M为多径信号的数目，a₀，τ₀为直接信号的幅度和延迟，a₁...a_M为多径信号1到多径信号M的幅度，τ₁...τ_M为多径信号1到多径信号M的延迟； 

步骤3：根据步骤1中的多径信号可能的数目，和步骤2中的状态方程和观测方程，采用粒子滤波方法估计出直接信号和多径信号的参数，同时可以确定多径的数目。 

所述的参数α为1～100。 

有益效果 

本发明提出一种卫星导航直接信号参数与多径信号数目和参数的估计方法，可以在不了解多径信号数目的情况下估计出多径信号的数目和参数，不但适用于直接信号大于多径信号的情况，还适用于直接信号小于多径信号的情况(经常发生在室内定位情况)。 

以GPS***为例通过仿真验证方法的效果，仿真分两种情况：室外定位情况和室内定位情况，GPS室外定位时，信号强度较大，多径信号较少，多径信号强度小于直接信号，即|α_l|＜|α₀|(τ_l＞τ₀， 

)，这里α₀，τ₀为直接信号的幅度和延迟，α_l，τ_l为多径信号的幅度和延迟；室内定位时，GPS信号强度较小，而且多径信号强度可能大于直接信号。 

附图说明

图1：给出了室外情况的信号延迟估计误差，，单位为米。图中(a)代表直接信号延迟的估计误差，(b)代表第一路多径信号延迟的估计误差，(c)代表第二路多径信号延迟的估计误差。 

图2：给出了室外情况的信号幅度估计误差，图中(a)代表直接信号幅度的估计误差，(b)代表第一路多径信号幅度的估计误差，(c)代表第二路多径信号幅度的估计误差。 

图3：给出了室外情况的信号延迟的真实值及其估计值，单位为米。图中(a)代表直接信号的延迟及其估计，(b)代表第一路多径信号的延迟及其估计，(c)代表第二路多径信号的延迟及其估计。实线代表真实延迟，虚线代表估计值。 

图4：给出了室外情况的信号幅度真实值及其估计值，图中(a)代表直接信号的幅度及其估计，(b)代表第一路多径信号的幅度及其估计，(c)代表第二路多径信号的幅度及其估计。实线代表真实幅度，虚线代表估计值。 

图5：给出了室外情况的信号数目估计，其中“×”为真实信号数目，包括直接信号和多径信号，“o”为估计的信号数目。 

图6：给出了室外情况20次实验的信号延迟估计RMS误差，单位为米。图中(a)代表直接信号延迟估计的RMS误差，(b)代表第一路多径信号延迟估计的RMS误差，(c)代表第二路多径信号延迟估计的RMS误差。 

图7：给出了室外情况20次实验的信号幅度估计RMS误差，图中(a)代表直接信号幅度估计的RMS误差，(b)代表第一路多径信号幅度估计的RMS误差，(c)代表第二路多径信号幅度估计的RMS误差。 

图8：给出了室内情况的信号延迟估计误差，单位为米。图中(a)代表直接信号 延迟的估计误差，(b)代表第一路多径信号延迟的估计误差，(c)代表第二路多径信号延迟的估计误差，(d)代表第三路多径信号延迟的估计误差。 

图9：给出了室内情况的信号幅度估计误差，图中(a)代表直接信号幅度的估计误差，(b)代表第一路多径信号幅度的估计误差，(c)代表第二路多径信号幅度的估计误差，(d)代表第三路多径信号幅度的估计误差。 

图10：给出了室内情况的信号延迟的真实值及其估计值，单位为米。图中(a)代表直接信号的延迟及其估计，(b)代表第一路多径信号的延迟及其估计，(c)代表第二路多径信号的延迟及其估计，(d)代表第三路多径信号的延迟及其估计。实线代表真实延迟，虚线代表估计值。 

图11：给出了室内情况的信号幅度真实值及其估计值，图中(a)代表直接信号的幅度及其估计，(b)代表第一路多径信号的幅度及其估计，(c)代表第二路多径信号的幅度及其估计，(d)代表第三路多径信号的幅度及其估计。实线代表真实幅度，虚线代表估计值。 

图12：给出了室内情况的信号数目估计，其中“×”为真实信号数目，包括直接信号和多径信号，“o”为估计的信号数目。 

图13：给出了室内情况10次实验的信号延迟估计RMS误差，单位为米。图中(a)代表直接信号延迟估计的RMS误差，(b)代表第一路多径信号延迟估计的RMS误差，(c)代表第二路多径信号延迟估计的RMS误差，(d)代表第三路多径信号延迟估计的RMS误差。 

图14：给出了室内情况10次实验的信号幅度估计RMS误差，图中(a)代表直接信号幅度估计的RMS误差，(b)代表第一路多径信号幅度估计的RMS误差，(c)代表第二路多径信号幅度估计的RMS误差，(d)代表第三路多径信号幅度估计的RMS误差。

具体实施方式

现结合实施例、附图对本发明作进一步描述： 

本发明提出的方法采用Dirichlet过程对多径信号的数目建模，然后根据多径信号的数目构造不同的模型，来对多径信号的参数进行估计。 

采用Dirichlet过程对多径信号的数目进行建模的好处在于能够更好的估计出变化的多径信号数目。 

本发明提出的新方法以基于状态空间的贝叶斯估计技术为基础，采用粒子滤波方法进行后验概率的计算。其数学模型可表示为： 

$\left\{\begin{array}{c}{x}_k=f(x_{k-1},r_k,v_k)\\ z_k=h(x_k,r_k,w_k)\end{array}\right.---\left(1\right)$

其中，x_k为***状态，表示了待估计量，包括直接信号的延迟、幅度和多径信号的延迟、幅度，f_k为状态方程，表示了状态的变化，v_k为***噪声，z_k为观测量，h_k为观测方程，w_k为观测噪声。当多径信号的数目不确定时，状态x_k的维数也不确定，因此***模型结构维数也不确定，状态方程和观测方程的结构维数不确定性可以由模态变量r_1:k表示，r_1:k也代表多径数目。 

由于在多径环境中，***模型不断变化，在新的情况下可能需要新的模型来描述***，r_1:k所代表的多径数目不可能是固定的，根据新的情况需要不断的调整变化，所以采用Dirichlet过程作为r_1:k的先验分布： 

G～DP(α，G₀) 

r_k～G                                     (2) 

$r_{k+1}|r_1,r_2,...,r_k~\frac{1}{\mathrm{\α}+k}(\mathrm{\α}{G}_0+\underset{i=1}{\overset{k}{\mathrm{\Σ}}}\mathrm{\δ}{r}_i)$

采用Dirichlet过程作为先验分布可以不需要事先假设多径信号的数目，在估计过程中，多径信号的数目与多径信号参数可以同时被估计出，当多径信号数目变化时，其变化也可以被估计出来。在(2)式中，DP为Dirichlet过程，G₀为Dirichlet过程的基准分布，α为参数，～表示从...采样，δ为Diracδ函数。 

以(1)为基础，r_1:k采用(2)的模型，由状态量x_k表示的直接信号的延迟、幅度和多径信号的延迟、幅度就可以采用粒子滤波等贝叶斯估计方法估计出。 

下面我们以GPS多径信号估计为例，提出具体的方法。当不考虑多径时，接收到的GPS信号可以表示为： 

$s\left(t\right)=\sqrt{2P\left(t\right)}d(t-\mathrm{\τ})c(t-\mathrm{\τ})\sin (({\mathrm{\ω}}_c+\mathrm{\Δ\ω})t+{\mathrm{\φ}}_0)+n_T\left(t\right)---\left(3\right)$

上式中，d(t-τ)为导航数据，c(t-τ)为伪随机码，P(t)为信号功率，τ为码相位延迟，ω_c为载波频率，Δω为多普勒频移，φ₀为未知初始相位，n_T(t)为热噪声。GPS信号在接收机中，经过前端处理，下变频后，就成为中频信号，将中频信号采样后，得到： 

$s\left(t_k\right)=\sqrt{2P\left(t_k\right)}d\left(t_k\right)c\left(t_k\right)\sin (({\mathrm{\ω}}_{\mathrm{IF}}+\mathrm{\Δ\ω})t_k+{\mathrm{\φ}}_0)+n_T\left(t_k\right)---\left(4\right)$

上式中，ω_IF为中频信号频率，t_k为采样时刻。采样后的中频信号和接收机本地同相信号 

Ac(t_k+δ)sin(ω_IFt_k+φ_r)    (5) 

及正交信号 

Ac(t_k+δ)cos(ω_IFt_k+φ_r)    (6) 

相关，得： 

$\sqrt{2P\left(t_k\right)}\mathrm{Ac}\left(t_k\right)c(t_k+\mathrm{\δ})d\left(t_k\right)\sin (({\mathrm{\ω}}_{\mathrm{IF}}+\mathrm{\Δ\ω})t_k+{\mathrm{\φ}}_0)\sin ({\mathrm{\ω}}_{\mathrm{IF}}{t}_k+{\mathrm{\φ}}_r)+n_I^{\′}$

$=\frac{\sqrt{2P\left(t_k\right)}}{2}\mathrm{Ac}\left(t_k\right)c(t_k+\mathrm{\δ})d\left(t_k\right)\cos ({\mathrm{\Δ\ωt}}_k+\mathrm{\Δ\φ})+n_I^{\′}---\left(7\right)$

和 

$\sqrt{2P\left(t_k\right)}\mathrm{Ac}\left(t_k\right)c(t_k+\mathrm{\δ})d\left(t_k\right)\sin (({\mathrm{\ω}}_{\mathrm{IF}}+\mathrm{\Δ\ω})t_k+{\mathrm{\φ}}_0)\cos ({\mathrm{\ω}}_{\mathrm{IF}}{t}_k+{\mathrm{\φ}}_r)+n_I^{\′}$

$=\frac{\sqrt{2P\left(t_k\right)}}{2}\mathrm{Ac}\left(t_k\right)c(t_k+\mathrm{\δ})d\left(t_k\right)\sin ({\mathrm{\Δ\ωt}}_k+\mathrm{\Δ\φ})+n_I^{\′}---\left(8\right)$

上式中，δ为本地信号与实际信号的相对延迟，A为本地信号幅度，φ_r为本地信号初始相位，n′_I为噪声，Δφ为相位差，设A＝2，则上述相关值的累加输出为： 

$I_{\mathrm{\δ}}\left(t_k\right)=\sqrt{2P\left(t_k\right)}{M}^{\′}R(\mathrm{\τ}+\mathrm{\δ})d\left(t_k\right)\sin (\mathrm{\Δ\ωT}+\mathrm{\Δ\φ})\sin c\left(\mathrm{\ΔfT}\right)+n_I---\left(9\right)$

$Q_{\mathrm{\δ}}\left(t_k\right)=\sqrt{2P\left(t_k\right)}{M}^{\′}R(\mathrm{\τ}+\mathrm{\δ})d\left(t_k\right)\cos (\mathrm{\Δ\ωT}+\mathrm{\Δ\φ})\sin c\left(\mathrm{\ΔfT}\right)+n_Q---\left(10\right)$

上式中，I_δ为同相信号输出，Q_δ为正交信号输出，M′为累加采样数，R(·)为自相关函数，n_I，n_Q为噪声，T为累加时间，Δf＝Δφ/2π，sin c(·)为sinc函数，其定义为： 

$\sin c\left(\mathrm{\ΔfT}\right)=\frac{\sin \left(\mathrm{\π\ΔfT}\right)}{\mathrm{\π\ΔfT}}---\left(11\right)$

在GPS接收机中，由(9)和(10)给出的相关累加输出将被输入到鉴别器，并最终获得信号延迟的测量值。 

当存在多径时，情况将变得比较复杂。由于在多径环境下，信道可建模为时变***： 

$h(t,\mathrm{\τ})=\underset{n=0}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\α}}_n\left(t\right)\mathrm{\δ}(\mathrm{\τ}-{\mathrm{\τ}}_n)---\left(12\right)$

其中M代表多径的数目，α_n为多径系数，τ_n为多径延迟，这里δ为冲激函数。因此， 

当考虑到多径信号后，接收到的中频信号变为： 

$s\left(t\right)=\sqrt{{2P}_0\left(t\right)}d(t-{\mathrm{\τ}}_0)c(t-{\mathrm{\τ}}_0)\sin (({\mathrm{\ω}}_{\mathrm{IF}}+{\mathrm{\Δ\ω}}_0)t+{\mathrm{\φ}}_0)$

$+\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\sqrt{{2P}_i\left(t\right)}d(t-{\mathrm{\τ}}_i)c(t-{\mathrm{\τ}}_i)\sin (({\mathrm{\ω}}_{\mathrm{IF}}+{\mathrm{\Δ\ω}}_i)t+{\mathrm{\φ}}_i)+n_T\left(t\right)---\left(13\right)$

上式中，下标i代表第i路多径信号。当本地信号延迟为δ时，上述带有多径的接收信号的相关累加输出为： 

$I_{\mathrm{\δ}}\left(t_k\right)=\sqrt{{2P}_0\left(t_k\right)}{M}^{\′}R(\mathrm{\τ}-{\mathrm{\τ}}_0+\mathrm{\δ})d(t_k-{\mathrm{\τ}}_0)\sin ({\mathrm{\Δ\ω}}_{0}T+{\mathrm{\Δ\φ}}_0)\sin c\left(\mathrm{\ΔfT}\right)$

$+\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\sqrt{{2P}_i\left(t_k\right)}{M}^{\′}R(\mathrm{\τ}-{\mathrm{\τ}}_i+\mathrm{\δ})d(t_k-{\mathrm{\τ}}_i)\sin ({\mathrm{\Δ\ω}}_{i}T+{\mathrm{\Δ\φ}}_i)\sin c\left(\mathrm{\ΔfT}\right)+n_I---\left(14\right)$

$Q_{\mathrm{\δ}}\left(t_k\right)=\sqrt{{2P}_0\left(t_k\right)}{M}^{\′}R(\mathrm{\τ}-{\mathrm{\τ}}_0+\mathrm{\δ})d(t_k-{\mathrm{\τ}}_0)\cos ({\mathrm{\Δ\ω}}_{0}T+{\mathrm{\Δ\φ}}_0)\sin c\left(\mathrm{\ΔfT}\right)$

$+\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\sqrt{{2P}_i\left(t_k\right)}{M}^{\′}R(\mathrm{\τ}-{\mathrm{\τ}}_i+\mathrm{\δ})d(t_k-{\mathrm{\τ}}_i)\cos ({\mathrm{\Δ\ω}}_{i}T+{\mathrm{\Δ\φ}}_i)\sin c\left(\mathrm{\ΔfT}\right)+n_I---\left(15\right)$

为了简便起见，我们合并同相和正交分量，并写成指数形式，则(14)和(15)可合并为： 

$y_{\mathrm{\δ}}\left(t_k\right)={\mathrm{\α}}_0\left(t_k\right)R(\mathrm{\τ}-{\mathrm{\τ}}_0+\mathrm{\δ})e^{{\mathrm{j\Δ\φ}}_0}+\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\α}}_i\left(t_k\right)R(\mathrm{\τ}-{\mathrm{\τ}}_i+\mathrm{\δ})e^{{\mathrm{j\Δ\φ}}_i}+n\left(t_k\right)---\left(16\right)$

在上式中， 

n为噪声，同时认为锁相环能够完全跟踪载波频率，这样就可以忽略频率误差，并假设在积分累加周期T内，数据位保持不变，即d＝1。 

方程(16)表出了观测量的基础，通过获得不同延迟δ的多路相关累加支路的输出信息，估计器得到所需的观测向量，若设有N路相关累加支路，它们的延迟分别是δ₁，δ₂，...，δ_N，则在t_k时刻，N维观测向量的一个分量为： 

$y_{{\mathrm{\δ}}_j}\left(t_k\right)={\mathrm{\α}}_0\left(t_k\right)R(\mathrm{\τ}-{\mathrm{\τ}}_0+{\mathrm{\δ}}_j)e^{{\mathrm{j\Δ\φ}}_0}+\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\α}}_i\left(t_k\right)R(\mathrm{\τ}-{\mathrm{\τ}}_i+{\mathrm{\δ}}_j)e^{{\mathrm{j\Δ\φ}}_i}+n\left(t_k\right)---\left(17\right)$

其中，δ_j∈{δ₁，δ₂，...，δ_N}，N维观测向量为 

方程(17)给出了观测向量，其中α₀，τ₀，Δφ₀代表直接信号的大小，延迟和相位， α_i，τ_i，Δφ_i，i＝1，2，...，M代表多径信号的大小，延迟和相位，注意这些量都是时变的，α₀，τ₀，Δφ₀和α_i，τ_i，Δφ_i，i＝1，2，...，M为待估计量，构成状态向量x。如果将信号幅度写为复数形式，即令 

和 i＝1，2，...，M，则状态量变为a₀，τ₀和a_i，τ_i，i＝1，2，...，M，由此列出状态方程： 

$a_0(k+1)={\mathrm{\λ}}_0{a}_l\left(k\right)+w_{a_0}\left(k\right)$

${\mathrm{\τ}}_0(k+1)={\mathrm{\γ}}_0{\mathrm{\τ}}_0\left(k\right)+w_{{\mathrm{\τ}}_0}\left(k\right)$ (18) 

$a_i(k+1)={\mathrm{\λ}}_i{a}_i\left(k\right)+w_{a_i}\left(k\right)$

${\mathrm{\τ}}_i(k+1)={\mathrm{\γ}}_i{\mathrm{\τ}}_i\left(k\right)+w_{{\mathrm{\τ}}_i}\left(k\right)$

其中，i＝1，2，...，M，括号中的k+1代表t_k+1时刻，k代表t_k时刻，令 

且 $w\left(k\right)={[w_{a_0}\left(k\right),w_{a_1}\left(k\right),...,w_{a_M}\left(k\right),w_{{\mathrm{\τ}}_0}\left(k\right),w_{{\mathrm{\τ}}_1}\left(k\right),...,w_{{\mathrm{\τ}}_M}\left(k\right)]}^T,$ 则***状态空间模型可写为： 

x(k+1)＝Fx(k)+w(k) 

                                    (19) 

y(k)＝H(x(k))+n(k) 

在上式中，F＝diag{F_a，F_τ}，其中 

$F_a=\left[\begin{array}{cccc}{\mathrm{\λ}}_0& 0& ...& 0\\ 0& {\mathrm{\λ}}_1& ...& 0\\ ...& ...& ...& ...\\ 0& 0& ...& {\mathrm{\λ}}_M\end{array}\right],$ $F_{\mathrm{\τ}}=\left[\begin{array}{cccc}{\mathrm{\γ}}_0& 0& ...& 0\\ 0& {\mathrm{\γ}}_1& ...& 0\\ ...& ...& ...& ...\\ 0& 0& ...& {\mathrm{\γ}}_M\end{array}\right]$

上两矩阵对角线元素参考取值为0.9～1。 

w(k)为高斯白噪声，其协方差阵为∑_k，x＝diag{∑_k，a，∑_k，τ}，其中∑_k，a，∑_k，τ为对角阵，对角线元素参考取值为0.0001～0.01，n(k)也为高斯白噪声，协方差阵为∑_k，R，取值由信号信噪比确定，F为线性矩阵，H为非线性函数，形式为(17)式。 

 如果我们将t_k时刻的状态向量和观测向量分别写为x_k和y_k，则GPS多径信号的估计问题就可表述为对条件概率p(x_k|y_k)的估计，在这里我们需要注意，当假设的多径数目不同时，即M不同时，状态的维数不同，这将使状态方程的结构不同，因此，在对p(x_k|y_k)进行估计时，必须首先确定多径信号的个数M，如果不能事先确定M，就不能正确的估计出状态向量。 

问题的模型由(19)和(17)给出，为了减少运算量，我们可以采用一种基于最小二乘的降维技巧，为了应用这种技术，我们首先将(17)写为矩阵形式： 

y＝aR(τ)+n    (20) 

其中，y为观测向量，a为复向量， 

R(τ)为包含直接信号和多径信号延迟的相关函数矩阵，n为噪声。 

根据式(20)，可以利用最小二乘法降低状态维数，因为当τ给定时，a的最小二乘估计值为： 

${\hat{a}}_{\mathrm{LM}}=\mathrm{yR}{\left(\mathrm{\τ}\right)}^H{\left(R\left(\mathrm{\τ}\right)R{\left(\mathrm{\τ}\right)}^H\right)}^{-1}---\left(21\right)$

因此，在得到τ的估计值 

时，可以利用(21)得到a的估计值 

$\hat{a}=\mathrm{yR}{\left(\widehat{\mathrm{\τ}}\right)}^H{\left(R\left(\widehat{\mathrm{\τ}}\right)R{\left(\widehat{\mathrm{\τ}}\right)}^H\right)}^{-1}---\left(22\right)$

通过这样的处理，我们可以将状态维数减少一半，这样可以大大降低运算量。 

这里给出多径估计方法： 

方法利用粒子滤波来完成统计计算，在方法中，N为粒子数，τ₀为直接信号延迟，τ_mp为多径信号延迟，U为均匀分布，T_u为码元长度，方法首先进行初始化，从均匀分布中得到第一批直接信号延迟和多径信号延迟的粒子采样，并将权值设为1/N， 

表示多径信号延迟大于直接信号。G为Dirichlet过程先验，模态变量r_k表示为M(k)，令M(k)的先验分布符合Dirichlet过程： 

G～DP(α，G₀) 

                                            (23) 

M(k)～G 

令G₀取参数为p的几何分布Geo(p)，p的参考范围0＜p＜1，α的参考值1～100，则M(k)的取值范围为{1，2，...，∞}，其值代表多径的数目，即模型的结构，可以看到，多径数目的取值范围为1到∞，这就意味着该模型可以适用于多径数目未知的任何情况，由于采用几何分布作为基准分布，在对多径数目进行估计时，这种方法更倾向于选择多径数目较少的模型，这样做计算更简单，而且符合实际情况。 

在初始化后，就进入迭代计算阶段，利用(23)式，给出多径数目的重要性采样，并由状态方程给出延迟的采样，计算出a⁽ⁱ⁾(k)，然后更新权值。为了避免退化现象，方法还进行了重采样，并将权值重新置为1/N。 

首先给出室外定位情况的仿真结果。令信号C/N₀为43dBHz，这符合一般情况，除直接信号外，多径信号数目在0～2之间变化，事先并不知道多径信号数目。由于当相位差为0度或90度时，多径误差的影响最大，因此我们在仿真中假设相位差为0度或90度，这里假设相位差为0度。直接信号幅度a₀的初值设为1，两路多径信号的初值分别设为0.5，在仿真过程中，直接信号和多径信号的延迟是不断变化的，这表示我们的估计是在动态情况下的估计。为了符合实际情况，我们将直接信号幅度和多径信号幅度也设为可变的，这表示接收到的信号功率不断变化，相当于接收机载体穿过林地、街区等环境对GPS信号的衰减影响较大的地区。仿真中延迟的初值随机给定。23路相关累加支路，粒子数60000。 

给出了在室外环境下的多径估计结果，在这种情况下，直接信号功率较大，多径信号小于直接信号。给出了信号数目的估计，从图中可以看到，在仿真过程中，信号的数目是在不断变化的，采用本发明提出的方法，我们可以跟踪这样的变化，基本给出正确的信号数目估计，这样就确定了模型的结构。显示了各路信号的估计误差和跟踪结果，图中(b)，(c)的间断性与对信号实际数目的估计有关，当这一路信号不存 在时，各图中的(b)，(c)应为0，如果对信号数目的估计有误，则可能不一定为0。延迟和幅度估计误差中的尖峰是由于多径信号数目估计不准确(即模型不匹配)引起的，当估计出真实多径信号数目(即估计出真实的模型结构)后，估计误差将趋于平稳。给出了延迟和幅度估计的RMS误差，对误差求平均，可以得到直接信号、第一路多径信号、第二路多径信号的平均延迟分别为：2.6028米、15.1419米和40.7180米，考虑到延迟和幅度估计误差中的尖峰或野值对RMS的影响，这样的结果充分证明了方法的有效性，采用我们的方法可以在事先不确定多径信号数目的情况下很好的对直接信号和多径信号进行估计和跟踪。 

现在我们给出室内定位情况的仿真结果。令信号C/N₀为23dBHz，这样，信号功率比室外情况小20dB，这符合室内定位弱信号情况，除直接信号外，多径信号数目在0～3之间变化，事先也不知道多径数目。多径信号与直接信号之间信号的相位差假设为90度。直接信号幅度a₀的初值设为1，多径信号的初值分别设为2、2、1，这表明在室内环境下，多径信号的幅度大于直接信号，在仿真过程中，直接信号的幅度和延迟与多径信号的幅度和延迟也都是可变的。仿真中延迟的初值随机给定，采用23路相关累加支路，粒子数60000。 

给出了在室内环境下的多径估计结果，在这种情况下，直接信号功率微弱，多径信号大于直接信号。给出了信号数目的估计，从图中可以看到，采用本发明提出的方法，我们可以跟踪信号数目的变化，在仿真过程中实时估计模型的结构。显示了各路信号的估计误差和跟踪结果。同样，延迟和幅度估计误差中的尖峰是由于多径数目估计不准确(即模型不匹配)引起的，当估计出真实的模型结构后，估计误差将趋于平稳。给出了延迟和幅度估计的RMS误差，对误差求平均，可以得到直接信号、第一路多径信号、第二路多径信号和第三路多径信号的平均延迟分别为：2.6148米、13.8008米、32.3000米和42.8040米，由此可见，采用我们的方法可以在事先不确定多径数目的情况下很好的对直接信号和多径信号进行估计和跟踪。 

Claims (2)

1.一种卫星导航直接信号参数与多径信号数目和参数的估计方法，其特征在于步骤如下：

步骤1：采用Dirichlet过程DP(α，G₀)作为多径信号数目M的先验分布，则M的分布符合Dirichlet过程：

根据此先验分布采样得到多径信号可能的数目M，其中：α为参数，G₀为参数为p的几何分布Geo(p)，p的参考范围0＜p＜1；

步骤2：根据多径信号可能的数目，构建以直接信号和多径信号的参数为状态量的状态方程x(k+1)＝Fx(k)+w(k)和观测向量的观测方程y＝aR(τ)+n；

其中：k为离散时间变量；

F＝diag{F_a，F_τ}：

此两矩阵对角线元素参考取值为0.9～1；

w(k)为高斯白噪声：其协方差阵为∑_k，x＝diag{∑_k，a，∑_k，τ}，∑_k，a，∑_k，τ为对角阵，对角线元素参考取值为0.0001～0.01；

y为观测向量，其中的一个分量为：

N维观测向量为 

δ_j为产生第j个观测分量的卫星导航接收机延迟锁定环相关累加支路的延迟；

a为复向量，其中Δφ₁，Δφ₂，…Δφ_M为相位差；n(t_k)为t_k时刻的噪声；

R(τ)为包含直接信号和多径信号延迟的相关函数矩阵： 

这里R(·)为卫星信号和接收机自身参考信号的相关函数，n为噪声，大小由信噪比确定；

所述的直接信号和多径信号的参数为直接信号和多径信号的幅度向量和延迟向量；其中 

所述幅度向量为 

所述延迟向量为 

其中：M为多径信号的数目，a₀，τ₀为直接信号的幅度和延迟，a₁...a_M为多径信号1到多径信号M的幅度，τ₁...τ_M为多径信号1到多径信号M的延迟；

步骤3：根据步骤1中的多径信号可能的数目，和步骤2中的状态方程和观测方程，采用粒子滤波方法估计出直接信号和多径信号的参数，同时可以确定多径的数目。

2.根据权利要求1所述的卫星导航直接信号与多径信号数目和参数的联合估计方法，其特征在于：所述的参数α为1～100。 

Images