CN102590830A - 一种用于估计gps信号源的多普勒频率的方法 - Google Patents

Abstract

本发明属于卫星导航信号源的参数估计领域，提出了一种用于估计GPS信号源的多普勒频率的方法，包括根据当前时刻的用户速度、位置信息和导航历书计算当前时刻卫星位置信息、下一时刻用户位置信息和下一时刻卫星位置信息，利用这些位置信息分别计算当前时刻和下一时刻卫星和用户的相对距离，计算两时刻间卫星和用户的相对距离变化量，计算卫星和用户相对传输时延变化量，估算从当前时刻到下一时刻多普勒频率的平均值等步骤。该方法降低了运算量，更适合对多普勒频率进行连续的估计。

Description

一种用于估计GPS信号源的多普勒频率的方法

技术领域

本发明涉及卫星导航信号源的参数估计领域，具体涉及一种用于估计GPS信号源的多普勒频率的方法。

背景技术

GPS信号源可以为卫星导航领域的各项研究、仿真和测试提供可重复的模拟信号，有助于研究工作的推进。同时，在不便于架设天线的环境中，如偏远山区和气候恶劣的地区，GPS信号源也可以保证导航终端的正常工作，大大降低了架设天线的建设成本和维护费用。

在GPS信号源模拟产生GPS导航信号时，准确估计出多普勒频率是其中一个重要的步骤，是正确模拟GPS信号的基础。目前常用的方法是，利用用户和所选卫星的速度信息，计算卫星与用户的相对速度，并基于此计算信号的多普勒频率参数f_d，即

$f_d=f_T\×\frac{v_{u}g{v}_s}{c}$

其中，f_T为GPS信号的载波频率，c为光速，v_u和v_s分别为当前时刻t的用户和卫星的速度矢量，v_u·v_s计算的是当前时刻卫星与用户的相对速度。

在实际计算时，以上方法主要包括如下步骤：

步骤1：输入用户状态信息和所选卫星的导航历书，并设定多普勒频率更新间隔T_u。

用户状态信息包括：当前时刻时间t、以及对应的ECEF(地心地固)用户位置信息p_u＝[x_u，y_u，z_u]^T和速度信息v_u＝[v_ux，v_uy，v_uz]^T。如果输入的用户位置和速度信息不是ECEF(地心地固)坐标，则需要对其进行坐标转换，转换方法为本领域普通技术人员所熟知的内容。

导航历书可以是GPS导航技术规范上记载的历书，也可以是通过GPS导航终端接收的实际历书。常用的导航历书参数包括：卫星轨道离心率e、历书生成时刻t_oe、卫星轨道和赤道的倾角i₀、升交点赤经的变化率′Ω、卫星轨道长半轴a、每周历元Ω_e、近地点弧角ω、参考时刻的平近点角M₀、轨道倾角的变化率di/dt、平均角速度的校正值Δn、对纬度幅角余弦的校正值C_uc、对纬度幅角正弦的校正值C_us、对轨道半径余弦的校正值C_rc、对轨道半径正弦的校正值C_rs、对轨道倾角余弦的校正值C_ic、对轨道倾角正弦的校正值C_is等。

步骤2：利用当前时刻t和导航历书，计算当前时刻的卫星位置，具体包括如下步骤：

步骤2-1：计算归一化时间：

t_k＝t-t_oe

步骤2-2：计算卫星运行的平均角速度n_k

$n_k=n_0+\mathrm{\Δn}=\sqrt{G/a^3}+\mathrm{\Δn}$

其中，G＝3.986005×10¹⁴是地心引力常数，a是卫星椭圆轨道长半轴，Δn为平均角速度的校正值，均来自导航历书。

步骤2-3：计算卫星在t_k时刻的平近点角M_k

M_k＝n_k×t_k+M₀

其中，M₀为历书中提供的参考时刻t_oe的平近点角。

步骤2-4：迭代计算卫星在t_k时刻的偏近点角F_k，

(1)置初值，令E_k0＝M_k；

(2)迭代计算，令i＝1，进行如下计算：

E_ki＝M_k+esin E_k(i-1)，ΔE_k＝|E_ki-E_k(i-1)|

判断，若ΔE_k小于事先设定的门限，则停止迭代，E_k＝E_ki；

否则，令i＝i+1，返回第(2)步。

步骤2-5：计算卫星在t_k时刻的真近点角f_k，

$f_k={\cos }^{-2}\left(\frac{{\cos E}_k-e}{1-{e\cos E}_k}\right)\×\mathrm{sign}\left(\frac{\sqrt{1-e^2}\sin E_k}{1-{e\cos E}_k}\right)$

其中，卫星轨道离心率e由导航历书提供。

步骤2-6：计算未校准的升交点角距Φ_k

Φ_k＝f_k+ω

其中，ω为历书中提供的近地点弧角。

步骤2-7：计算摄动修正项δ_μ，δ_r和δ_i。

δμ，δr和δi分别为升交点角距Φ_k、卫星失径r_k和轨道倾角i_k的摄动量，其具体的计算方法如下式所示：

δ_u＝C_uc cos 2Φ_k+C_us sin 2Φ_k

δ_r＝C_rc cos 2Φ_k+C_rs sin 2Φ_k

δ_i＝C_ic cos 2Φ_k+C_is sin 2Φ_k

在此基础上，计算经过摄动改正的升交点角距Φ_k，卫星失径r_k，和轨道倾角i_d，具体计算公式为：

Φ_k＝Φ_k+δ_u，r_k＝a(1-ecos E_k)+δ_r，i_k＝i₀+(di/dt)×t_k+δ_i

其中，i₀为卫星轨道倾角，di/dt为轨道倾角的变化率，均由历书提供。

步骤2-8：计算卫星在t_k时刻的升交点赤经Ω_k

Ω_k＝Ω_e+′Ω×t_k

其中，Ω_e为历书中提供的参考时刻的升交点赤经，′Ω为历书提供的升交点赤经的变化率。

如考虑地球自转，还需要对其进行地球自转校正

Ω_k＝Ω_e+(′Ω-ω_e)t_k-ω_et_oe＝Ω_e+′Ωt_k-ω_et_in

其中，ω_e为地球自转速度，取值为ω_e＝7.2921151467×10^-5(rad/s)。

步骤2-9：计算t_k时刻卫星在ECEF(地心地固)坐标系的坐标

步骤2-9-1：计算卫星在椭圆轨道直角坐标系中的位置坐标

_xk＝r_kcosΦ_k，y_k＝r_ksinΦ_k

步骤2-9-2：计算卫星在ECEF坐标系中的坐标p_s＝[x_s，y_s，z_s]^T：

$\left[\begin{array}{c}{x}_s\\ y_s\\ z_s\end{array}\right]=\left[\begin{array}{c}{x}_k{\cos \mathrm{\Ω}}_k-y_k{\cos i}_k\×{\sin \mathrm{\Ω}}_k\\ x_k{\sin \mathrm{\Ω}}_k+y_k{\cos i}_k\×{\cos \mathrm{\Ω}}_k\\ y_k{\sin i}_k\end{array}\right]$

步骤3：利用当前时刻卫星的位置信息p_s，对其进行求导，计算当前时刻卫星的速度信息v_s

$v_s=\left[\begin{array}{c}{v}_{\mathrm{sx}}\\ v_{\mathrm{sy}}\\ v_{\mathrm{sz}}\end{array}\right]=\left[\begin{array}{c}{-y}_s{\mathrm{\Ω}}_k^{\′}-(y_k^{\′}{\cos i}_k-z_s{i}_k^{\′})\·{\sin \mathrm{\Ω}}_k+x_k^{\′}{\cos \mathrm{\Ω}}_k\\ x_s{\mathrm{\Ω}}_k^{\′}+(y_k^{\′}{\cos i}_k-z_s{i}_k^{\′})\·{\cos \mathrm{\Ω}}_k+x_k^{\′}{\sin \mathrm{\Ω}}_k\\ y_k^{\′}{\sin i}_k+y_k{i}_k^{\′}{\cos i}_k\end{array}\right]$

上式中的x_k’和y_k’是由x_k和y_k对时间求导得到，即

x′_k＝r′_k cos u_k-r_ku′_k sin u_k，y′_k＝r′_k sin u_k+r_ku′_k cos u_k

其中u_k’，r_k’，i_k’和Ω_k’可分别由u_k，r_k，i_k和Ω_k对时间求导得到，即

u′_k＝Φ′_k+δ′_u，r′_k＝aeE′_ksin E_k+δ′_r

i′_k＝di/dt+δ′_i，Ω′_k＝Ω′-ω_e

其中δ_μ’，δ_r’和δ_i’可分别由δ_μ，δ_r和δ_i对时间求导得到，即

δ′_u＝2Φ′_k[C_us cos(2Φ_k)-C_uc sin(2Φ_k)]

δ′_u＝2Φ′_k[C_rs cos(2Φ_k)-C_rc sin(2Φ_k)]

δ′_u＝2Φ′_k[C_is cos(2Φ_k)-C_ic sin(2Φ_k)]

其中Φ_k′可由Φ_k对时间求导得到，即

${\mathrm{\Φ}}_k^{\′}=f_k^{\′}=\frac{\sqrt{1-e^2}{E}_k^{\′}}{1-e\cos E_k}$

其中E_k可由E_k对时间求导得到，即

$E_k^{\′}=\frac{M_k^{\′}}{1-e\cos E_k}$

其中M_k′可由M_k对时间求导得到，即

M′_k＝n_k

步骤4：利用当前时刻用户和卫星的速度信息，计算卫星信号的多普勒频率。

$f_d=f_T\×\frac{v_{u}g{v}_s}{c}$

步骤5：如果到达多普勒频率更新间隔T_u，则更新当前时刻为t＝t+T_u，以及当前时刻的用户位置信息p_u和速度信息v_u，返回步骤2。

上述方法，在电子工业出版社2007年出版的由JAMES BAO-YEN TSUI著，陈军等译.GPS软件接收机基础(第2版)中有详细记载，在电子工业出版社2009年出版的由谢钢著的GPS原理与接收机设计中也有详细记载。

我们经过分析发现，以上方法在具体实现时存在如下难点：

真实卫星导航信号的多普勒频率都是连续变化的，而在使用以上方法对导航信号进行模拟时，只能以多普勒频率的离散值来逼近其连续变化的值。这就存在一个参数的更新率问题。理论上，参数的更新率越高，越能逼近其真实的变化情况。但是，参数的更新率越高，对实现系统的性能和成本要求也就越高，因为这要求系统能够在更短的时间内完成参数的计算，即将上述多普勒频率的估计方法重新计算一次。在以上方法中，多普勒频率的计算需要卫星和用户的速度信息，计算卫星的速度需要在求得其位置信息后，进一步对位置求导来获得，运算复杂度较高。因此，多普勒频率参数的更新较为复杂，降低该参数更新的运算量有助于简化系统的实现难度和降低其硬件平台的成本。

发明内容

本发明提供一种用于估计GPS信号源的多普勒频率的方法。和现有技术相比，该方法更为简单，仅使用用户和卫星的位置信息来对多普勒频率进行估计，不再需要计算卫星的速度信息。

本发明提供的一种用于估计GPS信号源的多普勒频率的方法，包括利用卫星的导航历书参数计算该卫星在当前时刻t₁和下一时刻t₂的位置信息，分别记为p_s1＝[x_s1，y_s1，z_s1]^T和p_s2＝[x_s2，y_s2，z_s2]^T，并利用用户在当前时刻t₁的速度信息v_u1＝[v_xu1，v_yu1，v_zu1]^T和位置信息p_u1＝[x_u1，y_u1，z_u1]^T，估计用户在下一时刻t₂的位置信息p_u2＝[x_u2，y_u2，z_u2]^T＝[x_u1+v_xu1×(t₂-t₁)，y_u1+v_yu1×(t₂-t₁)，z_u1+y_zu1×(t₂-t₁)]^T，其特征是还依次包括如下步骤：

步骤1：根据用户的位置信息和卫星的位置信息，分别计算出在当前时刻t₁和下一时刻t₂时，用户和卫星的相对距离

${\mathrm{\ρ}}_1={\left|\right|p_{u1}-p_{s1}\left|\right|}_2=\sqrt{{(x_{u1}-x_{s1})}^2+{(y_{u1}-y_{s1})}^2+{(z_{u1}-z_{s1})}^2}$

${\mathrm{\ρ}}_2={\left|\right|p_{u2}-p_{s2}\left|\right|}_2=\sqrt{{(x_{u2}-x_{s2})}^2+{(y_{u2}-y_{s2})}^2+{(z_{u2}-z_{s2})}^2}$

步骤2：计算从t₁时刻到t₂时刻，用户与卫星相对距离的变化量Δρ＝ρ₂-ρ₁；

步骤3：计算从t₁时刻到t₂时刻，用户与卫星相对传输时延的变化量Δτ＝Δρ/c，其中，c为光速；

步骤4：估计从t₁时刻到t₂时刻，多普勒频率的平均值f_d＝f_T×Δτ/T_u，其中，f_T为GPS信号的载波频率，多普勒频率更新的时间间隔T_u＝t₂-t₁。

所述步骤4完成后，更新当前时刻的值t₁＝t₂，更新当前时刻的用户位置信息p_u1＝p_u2，更新速度信息v_u1＝v_u2，更新用户和卫星的相对距离ρ₁＝ρ₂，并计算出下一时刻t₂的用户位置信息p_u2和卫星位置信息p_s2，然后返回步骤1。

通过以上对技术方案的描述可知，本发明在计算卫星的位置信息采用的是和现有技术相同的方法，但是后续的计算步骤采用了更新时间间隔内的平均多普勒频率来代替更新时刻点上的瞬时多普勒频率，因此不再需要计算更新时刻的卫星瞬时速度，仅利用更新时刻的卫星位置信息就可完成多普勒频率的估计，从而大大降低了多普勒频率的计算量。而且，连续两次多普勒频率估计之间，第一次估计计算出的“下一时刻”相关参数，正好是第二次估计的“当前时刻”相关参数，这样，连续两次估计得很多参数就可以共用，运算量能进一步降低，便于流水实现。

和现有技术相比，在相同的硬件平台上，本发明的技术方案可以支持更高的参数更新率，从而更准确地逼近这些参数连续变化的情况。或者说，在相同的参数更新率要求下，本发明的技术方案需要的运算量更少，可以在一个成本更低的平台上进行实现。

附图说明

附图1为本发明提出的一种用于估计GPS信号源的多普勒频率的方法的总体流程图。

具体实施方式

下面具体说明本发明技术方案的实施方式。

步骤1：输入用户状态信息和导航历书。用户当前时刻为t₁＝0s，用户ECEF位置坐标为p_u＝[x_u，y_u，z_u]^T，用户ECEF速度坐标为v_u＝[0，0，0]^T，即用户处于静止状态；导航历书为网站上下载的简易历书。同时设定其它参数，包括光速c＝3×10⁸m/s，GPS信号载波频率f_T＝1575.42MHz，多普勒频率更新间隔T_u＝1ms。

步骤2：利用用户当前时间信息t₁，结合导航历书参数，计算当前时刻所选GPS卫星的位置信息p_s＝[x_s，y_s，z_s]^T。具体的计算方法采用背景技术中提及的现有的计算方法。

步骤3：用户和卫星的位置信息p_u和p_s，计算当前时刻t₁用户和卫星的相对距离ρ₁：

${\mathrm{\ρ}}_1={\left|\right|p_u-p_s\left|\right|}_2=\sqrt{{(x_u-x_s)}^2+{(y_u-y_s)}^2+{(z_u-z_s)}^2}$

步骤4：计算多普勒频率的更新时刻t₂

t₂＝t₁+T_u＝0.001s

步骤5：利用导航电文，计算t₂时刻的卫星位置，仍记为p_s三[x_s，y_s，z_s]^T。具体的计算方法采用背景技术中提及的现有的计算方法。

步骤6：利用用户的位置和速度信息，计算t₂＝0.001s时刻的用户位置，仍记为p_u，由于用户静止，故位置保持不变。

p_u＝p_u+v_u×T_u＝[x_u y_u z_u]^T

步骤7：用户和卫星的位置信息p_u和p_s，计算时刻t₂＝0.001s用户和卫星的相对距离ρ₂：

${\mathrm{\ρ}}_2={\left|\right|p_u-p_s\left|\right|}_2=\sqrt{{(x_u-x_s)}^2+{(y_u-y_s)}^2+{(z_u-z_s)}^2}$

步骤8：计算从t时刻到t₁时刻，用户与卫星相对距离的变化量

Δρ＝ρ₂-ρ₁

步骤9：计算从t时刻到t₁时刻，用户与卫星相对传输时延的变化量

Δτ＝Δρ/c＝Δρ/(3×10⁸)

步骤10：估计从t₁时刻到t₂时刻，多普勒频率的平均值

f_d＝f_T×Δτ/T_u＝1575.42MHz×Δτ/1ms

步骤11：如果到达多普勒频率更新间隔T_u，则更新当前时刻为t₂＝t₁+T_u，以及当前时刻的用户速度信息v_u＝[0，0，0]^T，同时，令ρ₁＝ρ₂，返回步骤4，即可对多普勒频率进行连续估计。

由上述描述可以看出，连续两次多普勒频率估计之间，很多参数可以共用，这样运算量能进一步降低，同时便于流水实现。

Claims (2)

1.一种用于估计GPS信号源的多普勒频率的方法，包括利用卫星的导航历书参数计算该卫星在当前时刻t₁和下一时刻t₂的位置信息，分别记为p_s1＝[x_s1，y_s1，z_s1]^T和p_s2＝[x_s2，y_s2，z_s2]^T，并利用用户在当前时刻t₁的速度信息v_u1＝[v_x1，v_yu1，v_zu1]^T和位置信息p_u1＝[x_u1，y_u1，z_u1]^T，估计用户在下一时刻t₂的位置信息p_u2＝[x_u2，y_u2，z_u2]^T＝[x_u1+v_xu1×(t₂-t₁)，y_u1+v_yu1×(t₂-t₁)，z_u1+v_zu1×(t₂-t₁)}^T，其特征是还依次包括如下步骤：

步骤1：根据用户的位置信息和卫星的位置信息，分别计算出在当前时刻t₁和下一时刻t₂时，用户和卫星的相对距离

${\mathrm{\ρ}}_1={\left|\right|p_{u1}-p_{s1}\left|\right|}_2=\sqrt{{(x_{u1}-x_{s1})}^2+{(y_{u1}-y_{s1})}^2+{(z_{u1}-z_{s1})}^2}$

${\mathrm{\ρ}}_2={\left|\right|p_{u2}-p_{s2}\left|\right|}_2=\sqrt{{(x_{u2}-x_{s2})}^2+{(y_{u2}-y_{s2})}^2+{(z_{u2}-z_{s2})}^2}$

步骤2：计算从t₁时刻到t₂时刻，用户与卫星相对距离的变化量Δρ＝ρ₂-ρ₁；

步骤3：计算从t₁时刻到t₂时刻，用户与卫星相对传输时延的变化量Δτ＝Δρ/c，其中，c为光速；

步骤4：估计从t₁时刻到t₂时刻，多普勒频率的平均值f_d＝f_T×Δτ/T_u，其中，f_T为GPS信号的载波频率，多普勒频率更新的时间间隔T_u＝t₂-t₁。

2.根据权利要求1所述的一种用于估计GPS信号源的多普勒频率的方法，其特征是：所述步骤4完成后，更新当前时刻的值t₁＝t₂，更新当前时刻的用户位置信息p_u1＝p_u2，更新速度信息v_u1＝v_u2，更新用户和卫星的相对距离ρ₁＝ρ₂，并计算出下一时刻t₂的用户位置信息p_u2和卫星位置信息p_s2，然后返回步骤1。

Images