卫星导航系统的信号精度和定位服务可用性的评估方法

摘要

本发明提供一种卫星导航系统的信号精度和定位服务可用性的评估方法，该方法包括：从IGS组织的网站下载广播星历文件的数据和精密星历文件的数据；参照GNSS接口文件的描述，使用所述广播星历文件计算卫星概略坐标和卫星概略钟差；使用所述精密星历文件计算卫星精密坐标和卫星精密钟差；计算每一分钟的卫星坐标误差和卫星钟差误差；计算t时刻网格点用户测距误差；计算网格点用户测距误差精度，并针对全球，统计卫星导航系统的用户测距误差精度；计算t时刻网格点的水平定位精度因子和垂直定位精度因子；计算水平定位精度和水平定位服务可用性性能，计算垂直定位精度和垂直定位服务可用性性能。本发明可以适用于各种不同的卫星导航系统。

说明书

技术领域

本发明涉及一种卫星导航系统的信号精度和定位服务可用性的评估方 法，属于卫星导航技术领域。

背景技术

随着卫星导航系统技术的逐步发展，未来不仅有美国的全球定位系统 (Global Positioning System，简称：GPS)导航系统，还有俄罗斯的格洛 纳斯(GLONASS)导航系统、欧洲的伽利略(Galileo)导航系统、以及中国 的北斗导航系统。不同的卫星导航系统提供的卫星导航信号精度和定位服务 可用性不尽相同。卫星导航系统包括空间星座系统、地面控制系统和用户终 端系统。对用户定位而言，空间星座系统提供的卫星导航信号性能与传播环 境误差、以及用户接收机性能一道决定了用户的基本导航定位性能。其中， 卫星导航信号精度是卫星导航系统服务提供者唯一能够控制和保证的部分。 对用户定位可用性而言，不同的服务可用性决定了卫星导航是否能用于特定 的应用领域，并且是开发相应卫星增强系统的基础。为了评估不同卫星导航 系统的信号精度和服务可用性，有必要发展一种统一的评估方法。

发明内容

本发明提供一种卫星导航系统的信号精度和定位服务可用性的评估方 法，以实现针对不同卫星导航系统的统一评估。

本发明一方面提供一种卫星导航系统的信号精度和定位服务可用性的评 估方法，其中包括：

从IGS组织的网站下载广播星历文件的数据和精密星历文件的数据；

参照GNSS接口文件的描述，使用所述广播星历文件计算卫星概略坐标 (X_k Y_k Z_k)和卫星概略钟差Δt_sv；

使用所述精密星历文件计算卫星精密坐标(X_p，Y_p，Z_p)和卫星精密钟差 Δt_sv，p；

计算每一分钟的卫星坐标误差(ΔX＝X_k-X_p，ΔY＝Y_k-Y_p，ΔZ＝Z_k-Z_p)和卫星 钟差误差Δt＝Δt_sv-Δt_sv，p；

计算t时刻网格点用户测距误差 ${\mathrm{URE}}_{\mathrm{grid},t}=({\Delta }_x^{\mathrm{rs}}·\mathrm{\Delta X}+{\Delta }_y^{\mathrm{rs}}·\mathrm{\Delta Y}+{\Delta }_z^{\mathrm{rs}}·\mathrm{\Delta Z})-\mathrm{c\Delta t};$

计算网格点用户测距误差精度并针对全球，统计卫 星导航系统的用户测距误差精度

计算t时刻网格点的水平定位精度因子和垂直定位精度 因子$\mathrm{VDOP}=\sqrt{g_{33}^s};$

计算t时刻网格点处的水平定位精度HACCU_grid，t＝1.96HDOP*RMS_grid，并计 算网格点卫星导航系统的水平定位服务可用性性能 ${\mathrm{Avail}}_{\mathrm{Hor},\mathrm{grid}}=\frac{\underset{t}{\Sigma }\mathrm{Count}({\mathrm{HACCU}}_{\mathrm{grid},t}<17)}{\underset{t}{\Sigma }\mathrm{Count}\left(H_{\mathrm{grid},t}\right)}\times 100\%;$

计算t时刻网格点处的垂直定位精度VACCU_grid，t＝1.96VDOP*RMS_grid，并计 算网格点卫星导航系统的垂直定位服务可用性性能 ${\mathrm{Avail}}_{\mathrm{ver},\mathrm{grid}}=\frac{\underset{t}{\Sigma }\mathrm{Count}({\mathrm{VACCU}}_{\mathrm{grid},t}<37)}{\underset{t}{\Sigma }\mathrm{Count}\left({\mathrm{VACCU}}_{\mathrm{grid},t}\right)}\times 100\%.$

本发明另一方面提供一种卫星导航系统的信号精度和定位服务可用性 的评估装置，其中包括：

数据下载模块，用于从IGS组织的网站下载广播星历文件数据和精密星 历文件数据；

概略运算模块，用于参照GNSS接口文件的描述，使用所述广播星历文件 计算卫星概略坐标(X_k Y_k Z_k)和卫星概略钟差Δt_sv；

精密运算模块，用于使用所述精密星历文件计算卫星精密坐标(X_p，Y_p，Z_p) 和卫星精密钟差Δt_sv，p；

误差运算模块，用于计算每一分钟的卫星坐标误差 (ΔX＝X_k-X_p，ΔY＝Y_k-Y_p，ΔZ＝Z_k-Z_p)和卫星钟差误差Δt＝Δt_sv-Δt_sv，p，并计算t 时刻网格点用户测距误差${\mathrm{URE}}_{\mathrm{grid},t}=({\Delta }_x^{\mathrm{rs}}·\mathrm{\Delta X}+{\Delta }_y^{\mathrm{rs}}·\mathrm{\Delta Y}+{\Delta }_z^{\mathrm{rs}}·\mathrm{\Delta Z})-\mathrm{c\Delta t};$

精度运算模块，用于计算网格点用户测距误差精度并针对全球，统计卫星导航系统的用户测距误差精度

因子运算模块，用于计算t时刻网格点的水平定位精度因子和垂直定位精度因子

水平性能模块，用于计算t时刻网格点处的水平定位精度 HACCU_grid，t＝1.96HDOP*RMS_grid，并计算网格点卫星导航系统的水平定位服务可 用性性能${\mathrm{Avail}}_{\mathrm{Hor},\mathrm{grid}}=\frac{\underset{t}{\Sigma }\mathrm{Count}({\mathrm{HACCU}}_{\mathrm{grid},t}<17)}{\underset{t}{\Sigma }\mathrm{Count}\left(H_{\mathrm{grid},t}\right)}\times 100\%;$

垂直性能模块，用于计算t时刻网格点处的垂直定位精度 VACCU_grid，t＝1.96VDOP*RMS_grid，并计算网格点卫星导航系统的垂直定位服务 可用性性能${\mathrm{Avail}}_{\mathrm{ver},\mathrm{grid}}=\frac{\underset{t}{\Sigma }\mathrm{Count}({\mathrm{VACCU}}_{\mathrm{grid},t}<37)}{\underset{t}{\Sigma }\mathrm{Count}\left({\mathrm{VACCU}}_{\mathrm{grid},t}\right)}\times 100\%.$

本发明按照用户定位需求，根据卫星导航系统提供的实际导航数据和卫 星星座的实际布局，提出了评估导航数据的准确程度和定位服务可用性的可 行方案，可以适用于各种不同的卫星导航系统。

附图说明

图1为本发明所述卫星导航系统的信号精度和定位服务可用性的评估方 法实施例的流程图；

图2为图1所示步骤120的具体流程图；

图3为本发明所述卫星导航系统的信号精度和定位服务可用性的评估装 置实施例的结构示意图；

图4为图3所示精密运算模块30的具体结构示意图。

具体实施方式

图1为本发明所述卫星导航系统的信号精度的评估方法实施例的流程 图，如图所示，包括如下步骤：

步骤110，从IGS组织的网站下载广播星历文件数据和精密星历文件数 据。

其中，所述IGS组织是指国际GPS服务组织，是指导各国卫星导航技术 发展的国际组织机构；具体地，可以通过文件传输协议(File Transfer  Protocol，简称：FTP)实现下载。

下载后，还可以删除所述广播星历文件数据中的不健康数据，例如，IODC (Issue of Data，Clock，简称：时钟数据龄期)与IODE(Issue of Data， Ephemeris，简称：星历数据龄期)不一致的数据、卫星健康字不为0的数据、 卫星URA(User Range Accuracy，简称：测距精度指数)大于等于48米的数 据等。

步骤120，参照GNSS接口文件的描述，使用所述广播星历文件计算卫星 概略坐标(X_k Y_k Z_k)和卫星概略钟差Δt_sv；

其中，所述GNSS接口是指全球导航卫星系统(Global Navigation  Satellite System)接口；所述广播星历文件包含表1所示的卫星星历参数：

表1

步骤130，使用所述精密星历文件计算卫星精密坐标(X_p，Y_p，Z_p)和卫星精 密钟差Δt_sv，p。

具体地，IGS组织的网站公布的SP3(Standard Product 3 Orbit Format， 简称：标准轨道产品格式)精密星历文件按照一定时间间隔(例如15分钟) 给出了相应GPS时刻的卫星精密坐标与卫星精密钟差。性能评估系统可自定 义时间间隔(例如1分钟)计算相应时刻的用户测距误差。如果性能评估系 统自定义时间间隔与SP3精密星历文件给定的时间间隔不一致，可以通过拉 格朗日插值法，由精密星历文件计算选择的时刻的卫星精密坐标和卫星精密 钟差。一般为了降低时间边界上的插值误差，可以选择三天的精密星历文件 计算当中一天的卫星精密星历。拉格朗日插值的阶数可以自定义。以七阶拉 格朗日插值为例，某一时刻拉格朗日插值公式形式为：

$L_7\left(x\right)=\underset{i=0}{\overset{7}{\Sigma }}\underset{j\ne i}{\overset{7}{\underset{j=0}{\Pi }}}\frac{x-x_j}{x_i-x_j}{y}_i$

精密星历文件给出的卫星坐标是卫星质心坐标，而广播星历给出的卫星 坐标是卫星相位中心坐标。使用IGS定期公布的相位中心偏差数据，将拉格 朗日插值计算的卫星质心坐标在星固坐标系下转化为卫星相位中心坐标，该 卫星相位中心坐标即为卫星精密坐标(X_p，Y_p，Z_p)；根据精密星历文件数据按照 拉格朗日插值公式计算得到任意时刻的卫星精密钟差Δt_sv，p。

步骤140，计算每一分钟的卫星坐标误差 (ΔX＝X_k-X_p，ΔY＝Y_k-Y_p，ΔZ＝Z_k-Z_p)和卫星钟差误差Δt＝Δt_sv-Δt_sv，p。

步骤150，计算t时刻网格点用户测距误差。

具体地，将全球划分为5度x5度的网格，按照每个网格点的坐标，计算 卫星坐标误差到相应网格点的误差投影。对每个网格点，在计算过程中只考 虑仰角在5度以上的卫星。按照1分钟时间间隔对每个网格每个时刻的用户 测距误差进行计算，持续24小时。

假设网格点坐标为X_u，Y_u，Z_u，计算网格点至卫星的坐标向量为：

${\Delta }_y^{\mathrm{rs}}=(X_k-X_u)/\sqrt{{\Delta }_x^{\mathrm{rs}}·{\Delta }_x^{\mathrm{rs}}+{\Delta }_y^{\mathrm{rs}}·{\Delta }_y^{\mathrm{rs}}+{\Delta }_z^{\mathrm{rs}}·{\Delta }_z^{\mathrm{rs}}}$

${\Delta }_x^{\mathrm{rs}}=(Y_k-Y_u)/\sqrt{{\Delta }_x^{\mathrm{rs}}·{\Delta }_x^{\mathrm{rs}}+{\Delta }_y^{\mathrm{rs}}·{\Delta }_y^{\mathrm{rs}}+{\Delta }_z^{\mathrm{rs}}·{\Delta }_z^{\mathrm{rs}}}$

${\Delta }_z^{\mathrm{rs}}=(Z_k-Z_u)/\sqrt{{\Delta }_x^{\mathrm{rs}}·{\Delta }_x^{\mathrm{rs}}+{\Delta }_y^{\mathrm{rs}}·{\Delta }_y^{\mathrm{rs}}+{\Delta }_z^{\mathrm{rs}}·{\Delta }_z^{\mathrm{rs}}}$

则该网格点处时刻t时的用户测距误差为：

${\mathrm{URE}}_{\mathrm{grid},t}=({\Delta }_x^{\mathrm{rs}}·\mathrm{\Delta X}+{\Delta }_y^{\mathrm{rs}}·\mathrm{\Delta Y}+{\Delta }_z^{\mathrm{rs}}·\mathrm{\Delta Z})-\mathrm{c\Delta t}.$

步骤160，计算网格点用户测距误差精度，并针对全球，统计卫星导航 系统的用户测距误差精度。

具体地，首先要删除不符合GPS信号规范的异常值。删除URE大于某一 门限的样本，门限根据GPS信号规范选择，取为25米。然后针对每个网格点， 根据网格点观测到的所有卫星24小时内用户测距误差样本，统计网格点的用 户测距误差精度为针对全球，统计卫星导航系统的用户 测距误差精度为：

步骤170，计算t时刻网格点的水平定位精度因子和垂 直定位精度因子

具体地，首先，计算由网格点到所有可视卫星的观测矩阵：

$H=\left(\begin{array}{cccc}{l}_1& m_1& n_1& 1\\ l_2& m_2& n_2& 1\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ l_n& m_n& n_n& 1\end{array}\right)$

其中，l_j、m_j和n_j表示由近似位置指向第j颗卫星的单位矢量的方向余 弦；然后计算G＝(H^TH)^-1，并假设G的分量表示为，

$G=\left(\begin{array}{cccc}{g}_{11}& g_{12}& g_{13}& g_{14}\\ g_{12}& g_{22}& g_{23}& g_{24}\\ g_{13}& g_{23}& g_{33}& g_{34}\\ g_{14}& g_{24}& g_{34}& g_{44}\end{array}\right)$

上面计算G矩阵是基于ECEF(Earth-Centered Earth-Fixed，简称：地 心地固)坐标系的，需要转换为当地ENU(East North Upper，简称：东北 高)坐标系，转换到当地ENU坐标系后，G矩阵的各分量变为，

$g_{11}^s=g_{11}{\sin }^2\phi {\cos }^2\lambda +g_{22}{\sin }^2\phi {\sin }^2\lambda +g_{33}{\cos }^2\phi +$

$2g_{12}{\sin }^2\phi \sin \lambda \cos \lambda -2g_{13}\sin \phi \cos \phi \cos \lambda -2g_{23}\sin \phi \cos \phi \sin \lambda$

$g_{22}^s=g_{11}{\sin }^2\lambda -2g_{12}\sin \lambda \cos \lambda +g_{22}{\cos }^2\lambda$

$g_{33}^s=g_{11}{\cos }^2\phi {\cos }^2\lambda +g_{22}{\cos }^2\phi {\sin }^2\lambda +g_{33}{\sin }^{22}\phi +$

${2g}_{12}{\cos }^2\phi \sin \lambda \cos \lambda +{2g}_{13}\sin \phi \cos \phi \cos \lambda +{2g}_{23}\sin \phi \cos \phi \sin \lambda$

$g_{44}^s=g_{44}$

基于ENU坐标系的水平定位精度因子为：垂直定位精 度因子为：$\mathrm{VDOP}=\sqrt{g_{33}^s}.$

步骤180，计算t时刻网格点处的水平定位精度 HACCU_grid，t＝1.96HDOP*RMS_grid，并计算网格点卫星导航系统的水平定位服务 可用性性能${\mathrm{Avail}}_{\mathrm{Hor},\mathrm{grid}}=\frac{\underset{t}{\Sigma }\mathrm{Count}({\mathrm{HACCU}}_{\mathrm{grid},t}<17)}{\underset{t}{\Sigma }\mathrm{Count}\left(H_{\mathrm{grid},t}\right)}\times 100\%.$

其中，Count(*)是计数函数。

步骤190，计算t时刻网格点处的垂直定位精度 VACCU_grid，t＝1.96VDOP*RMS_grid，并计算网格点卫星导航系统的垂直定位服务可 用性性能${\mathrm{Avail}}_{\mathrm{ver},\mathrm{grid}}=\frac{\underset{t}{\Sigma }\mathrm{Count}({\mathrm{VACCU}}_{\mathrm{grid},t}<37)}{\underset{t}{\Sigma }\mathrm{Count}\left({\mathrm{VACCU}}_{\mathrm{grid},t}\right)}\times 100\%.$

以下参照图2，介绍上述步骤120的具体计算过程，如图所示，包括如 下步骤：

步骤121，计算长半轴A、平角速度n₀、从需要时刻到参考时刻的时间差 t_k、改正平角速度n和平近点角M_k。

其中，GM＝398600.5(km³/s²)为地球引力常数； t_k＝t-t_oe，时刻t必须是以GPST计的卫星信号发射时刻，它由信号接收时刻 减去传播时间延迟得到，当t_k＝t-t_oe大于302400秒时，应在t_k中减去一个GPS 星期相应的秒数604800秒，当t_k＝t-t_oe小于-302400秒时，应在t_k中加上604800 秒；n＝n₀+Δn；M_k＝M₀+n×t_k。

步骤122，根据等式M_k＝E_k-e×sinE_k迭代计算偏近点角E_k。

步骤123，根据如下方程组计算真近点角v_k，并计算升交点角距参数 φ_k＝v_k+ω，

$\left(\begin{array}{c}\cos v_k=\frac{(\cos E_k-e)}{(1-e\times \cos E_k)}\\ \sin v_k=\frac{\sqrt{1-e^2}\sin E_k}{(1-e\times \cos E_k)}\end{array}\right)$

步骤124，根据如下方程组计算周期改正项。

$\left(\begin{array}{c}\delta u_k=C_{\mathrm{us}}\sin \left(2{\phi }_k\right)+C_{\mathrm{uc}}\cos \left(2{\phi }_k\right)\\ \delta r_k=C_{\mathrm{rs}}\sin \left(2{\phi }_k\right)+C_{\mathrm{rc}}\cos \left(2{\phi }_k\right)\\ \delta i_k=C_{\mathrm{is}}\sin \left(2{\phi }_k\right)+C_{\mathrm{ic}}\cos \left(2{\phi }_k\right)\end{array}\right)$

步骤125，计算改正后的升交点角距参数u_k＝φ_k+δu_k，以及改正后的向经 r_k＝A(1-e cos E_k)+δr_k。

步骤126，计算卫星在轨道平面内的坐标x_k＝r_kcos(u_k)，y_k＝r_k sin(u_k)。

步骤127，计算改正后的倾角i_k＝i_o+δi_k+(IDOT)_tk，并计算升交点赤经 ${\Omega }_k={\Omega }_0+(\stackrel{·}{\Omega }-{\stackrel{·}{\Omega }}_e)t_k-{\stackrel{·}{\Omega }}_e{t}_{\mathrm{oe}}.$

其中，是地球自转的角速度(7.2921151467x1E-5rad/sec)。

步骤128，根据如下方程组计算所述卫星概略坐标(X_k Y_k Z_k)，并计算 所述卫星概略钟差$\Delta t_{\mathrm{sv}}=a_{f0}+a_{f1}(t-t_{\mathrm{oc}})+a_{f2}{(t-t_{\mathrm{oc}})}^2+\mathrm{Fe}\sqrt{A}\sin E_k-T_{\mathrm{GD}},$

$\left(\begin{array}{c}{X}_k=x_k\times \cos \left({\Omega }_k\right)-y_k\times \cos \left(i_k\right)\times \sin \left({\Omega }_k\right)\\ Y_k=x_k\times \sin \left({\Omega }_k\right)+y_k\times \cos \left(i_k\right)\times \cos \left({\Omega }_k\right)\\ Z_k=y_k\times \sin \left(i_k\right)\end{array}\right)$

其中，F为常数(-4.442807633x1E-10sec/sqrt(meter))。

本实施例所述方法按照用户定位需求，根据卫星导航系统提供的实际导 航数据和卫星星座的实际布局，提出了评估导航数据的准确程度和定位服务 可用性的可行方案，可以适用于各种不同的卫星导航系统。

图3为本发明所述卫星导航系统的信号精度和定位服务可用性的评估 装置实施例的结构示意图，用以实现上述方法，如图所示，该装置包括： 数据下载模块10、概略运算模块20、精密运算模块30、误差运算模块40、 精度运算模块50、因子运算模块60、水平性能模块70及垂直性能模块80， 图中的粗直线表示总线连接，采用其他连接方式也可是可以的。该装置的工 作原理如下：

数据下载模块10从IGS组织的网站下载广播星历文件数据和精密星历文 件数据，下载后，还可以通过数据删除模块90删除由数据下载模块10下载 的所述广播星历文件数据中的不健康数据，有关不健康数据可参见上述方法 实施例；概略运算模块20参照GNSS接口文件的描述，使用所述广播星历文 件计算卫星概略坐标(X_k Y_k Z_k)和卫星概略钟差Δt_sv，所述广播星历文件中的 卫星星历参数可参见表1。

精密运算模块30使用所述精密星历文件计算卫星精密坐标(X_p，Y_p，Z_p)和 卫星精密钟差Δt_sv，p。具体地，如图4所示，通过精密运算模块30中的插值 单元31对所述精密星历文件数据中的卫星质心坐标进行拉格朗日插值计算； 然后由转化单元32将经所述拉格朗日插值计算后的卫星质心坐标在星固坐 标系下转化为卫星相位中心坐标作为所述卫星精密坐标(X_p，Y_p，Z_p)；进而由钟 差单元33根据所述精密星历文件数据进行拉格朗日插值计算任意时刻的卫 星精密钟差Δt_sv，p。具体的插值公式可参见上述方法实施例的相关内容。

此后，由误差运算模块40计算每一分钟的卫星坐标误差 (ΔX＝X_k-X_p，ΔY＝Y_k-Y_p，ΔZ＝Z_k-Z_p)和卫星钟差误差Δt＝Δt_sv-Δt_sv，p，并计算t 时刻网格点用户测距误差${\mathrm{URE}}_{\mathrm{grid},t}=({\Delta }_x^{\mathrm{rs}}·\mathrm{\Delta X}+{\Delta }_y^{\mathrm{rs}}·\mathrm{\Delta Y}+{\Delta }_z^{\mathrm{rs}}·\mathrm{\Delta Z})-\mathrm{c\Delta t};$此 后，由精度运算模块50计算网格点用户测距误差精度并针对全球，统计卫星导航系统的用户测距误差精度此后， 由因子运算模块60计算t时刻网格点的水平定位精度因子和垂直定位精度因子

基于上述运算结果，由水平性能模块70计算t时刻网格点处的水平定位 精度HACCU_grid，t＝1.96HDOP*RMS_grid，并计算网格点卫星导航系统的水平定位服 务可用性性能${\mathrm{Avail}}_{\mathrm{Hor},\mathrm{grid}}=\frac{\underset{t}{\Sigma }\mathrm{Count}({\mathrm{HACCU}}_{\mathrm{grid},t}<17)}{\underset{t}{\Sigma }\mathrm{Count}\left(H_{\mathrm{grid},t}\right)}\times 100\%,$并由垂直性能模 块80计算t时刻网格点处的垂直定位精度VACCU_grid，t＝1.96VDOP*RMS_grid，并 计算网格点卫星导航系统的垂直定位服务可用性性能

${\mathrm{Avail}}_{\mathrm{ver},\mathrm{grid}}=\frac{\underset{t}{\Sigma }\mathrm{Count}({\mathrm{VACCU}}_{\mathrm{grid},t}<37)}{\underset{t}{\Sigma }\mathrm{Count}\left({\mathrm{VACCU}}_{\mathrm{grid},t}\right)}\times 100\%.$

本实施例所述装置按照用户定位需求，根据卫星导航系统提供的实际导 航数据和卫星星座的实际布局，提出了评估导航数据的准确程度和定位服务 可用性的可行方案，可以适用于各种不同的卫星导航系统。

本领域普通技术人员可以理解：实现上述方法实施例的全部或部分步骤 可以通过程序指令相关的硬件来完成，前述的程序可以存储于一计算机可读 取存储介质中，该程序在执行时，执行包括上述方法实施例的步骤；而前述 的存储介质包括：ROM、RAM、磁碟或者光盘等各种可以存储程序代码的介质。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其 限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术 人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或 者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技 术方案的本质脱离本发明各实施例技术方案的精神和范围。