CN101609140B - 一种兼容导航接收机定位系统及其定位方法 - Google Patents

Abstract

本发明公开了一种兼容导航接收机定位系统及其定位方法。本系统包括参数配置模块、数据预处理器、卫星参数解算模块、伪距修正模块、坐标转换模块、定位解算模块、卡尔曼滤波器和报文形成模块。首先读取接收机的配置参数，然后提取星历数据、历书数据和观测数据，接着获取有效卫星的参数信息，再获取接收机位置和接收机速度，最后将导航信息生成报文。本发明实现了卫星导航系统的选择使用，北斗二代导航系统的独立导航摆脱了非常时期对其他导航系统的依赖，多卫星导航系统的组合导航增加了有效卫星的数量，保证了卫星定位服务的连续性、可用性和安全性，增强了观测卫星的几何图形强度，提高了整个卫星定位系统的可靠性。

Description

一种兼容导航接收机定位系统及其定位方法

技术领域

本发明属于卫星导航领域，具体涉及一种兼容导航接收机定位系统及其定位方法。

背景技术

卫星导航能够向各类用户和运动平台实时提供准确、连续的位置、速度和时间信息。全球卫星导航定位系统(Global Navigation Satellite System，简称为GNSS)是第二代卫星导航定位系统，具有陆海空全能性、全球性、全天候、连续性和实时性的特点。

目前世界上已经存在的两大全球卫星导航定位系统是美国的全球定位系统(GlobalPositioning System，简称为GPS)系统和俄罗斯的全球导航卫星系统(Global NavigationSatellite System，简称为GLONASS)，正在设计建设阶段的有欧洲的“伽利略”卫星定位系统(Galileo GNSS，简称为Galileo)和中国的北斗二代导航系统(BeiDou-2，简称为BD2)。

GPS是目前发展最为成熟的导航系统，但始终受美军方控制，在发生危机时，随时可能被关闭；BD2目前为区域导航系统，存在12颗定位卫星，单BD2在高遮蔽角的情况下定位精度较差，难以满足高精度定位要求，但其作为我国具有独立知识产权的全球卫星导航定位系统，它的应用在非常时期是必要的。

各系统的卫星星座规模有限，独立应用时还存在一定问题。首先，定位精度低，并且位置精度随用户接收机所处的位置和时间产生较大变化。其次，可靠性差，特别是在干扰较强的环境中，可能出现短暂信号中断或捕获卫星数少的情况，定位精度将受到很大影响。组合导航的应用得到广泛的研究，但迄今为止，尚未有包含BD2的多卫星导航系统兼容导航接收机定位系统，尚未提出包含BD2的多卫星导航系统兼容导航接收机定位方法。

GPS、GLONASS、Galileo、BD2系统定义的参考坐标系基本是一致的，但由于不可避免的存在测轨跟踪站站址坐标误差和测量误差，定义的坐标系与实际使用的坐标系存在一定的差异。因此，实际上各系统所采用的坐标系两两之间都存在差异。因此，在不同卫星导航系统组合应用时，需要对坐标系进行转换。确定两个坐标系转换参数的数学模型有很多。具有代表性的模型为布尔莎模型(Bursa Model，简称Bursa模型)和莫洛金斯基模型(Molodensky Model，简称Molodensky模型)，Bursa模型较为常用。

7参数Bursa模型：

其中，ΔX为地球中心偏移量在B系统参考坐标系X轴方向上的分量，ΔY为地球中心偏移量在B系统参考坐标系Y轴方向上的分量，ΔZ为地球中心偏移量在B系统参考坐标系Z轴方向上的分量，θ_x为B系统参考坐标系X轴到A系统参考坐标系X轴的旋转角，θ_y为B系统参考坐标系Y轴到A系统参考坐标系Y轴的旋转角，θ_z为B系统参考坐标系Z轴到A系统参考坐标系Z轴的旋转角， $m=\frac{l_a-l_b}{l_b}$ 为比例因子，l_a和l_b分别是同一条基线在A系统和B系统两坐标系中的基线长度。

GPS、GLONASS、Galileo、BD2系统皆为测时-测距导航系统，以测码伪距观测方式为例，用户接收机通过观测由卫星发射的测距码到用户接收机的传播时间，乘以光速c后获得接收机与卫星间通视距离的观测量。基于这一测量原理，四系统对时间测量以及对时间同步的要求都比较高，使得不同导航系统为了确保各自系统运行的独立性、可靠性以及稳定性，都独立维持其系统时间参考系统的运行。即使采用相同精度的设备，采用同样的技术也很难实现严格同步，导致多星座组合导航应用中需要进行时间系统统一，通常的方法是增加观测量，将时间同步参数作为未知数求解，便可以实现不同系统间的时间系统统一。

目前，普通接收机通用的单机定位方法是基于线性化的迭代求解算法，因其在大部分情况下都能取得良好的定位效果，并且算法简单，计算量相对较小，便于计算机实现。

基于线性化的迭代求解算法基于接收机与卫星S^j间的伪距公式构造导航定位方法：

${\widetilde{\mathrm{\ρ}}}^j\left(t\right)=R^j\left(t\right)+\mathrm{c\δ}\left(t\right)-\mathrm{c\δ}{t}^j\left(t\right)+{\mathrm{\Δ}}_I^j\left(t\right)+{\mathrm{\Δ}}_T^j\left(t\right)$

$=\sqrt{{(x\left(t\right)-x^j\left(t\right))}^2+{(y\left(t\right)-y^j\left(t\right))}^2+{(z\left(t\right)-z^j\left(t\right))}^2}+\mathrm{c\δ}\left(t\right)-\mathrm{c\δ}{t}^j\left(t\right)+{\mathrm{\Δ}}_I^j\left(t\right)+{\mathrm{\Δ}}_T^j\left(t\right)$

其中，j为观测卫星S^j的编号，

为t时刻卫星S^j的观测伪距，R^j(t)为用户接收机与卫星S^j间的几何距离，δ(t)为t时刻用户接收机时钟与卫星导航系统的系统时之间的钟差，δt^j(t)为t时刻卫星S^j的星钟与卫星导航系统的系统时之间的钟差，cδt^j(t)称为星钟修正误差，Δ_I ^j(t)为t时刻电离层折射对测码伪距的影响，称为电离层修正误差，Δ_T ^j(t)为t时刻大气对流层折射对测码伪距的影响，称为对流层修正误差，(x^j(t)，y^j(t)，z^j(t))^T为卫星S^j的位置坐标，(x(t)，y(t)，z(t))^T为用户接收机的位置坐标，c为光速；其中x(t)、y(t)、z(t)、δ(t)为未知量。

设进行星钟误差、电离层误差、对流层误差修正后的观测伪距为ρ^j(t)，则有：

${\mathrm{\ρ}}^j\left(t\right)={\widetilde{\mathrm{\ρ}}}^j\left(t\right)+\mathrm{c\δ}{t}^j\left(t\right)-{\mathrm{\Δ}}_I^j\left(t\right)-{\mathrm{\Δ}}_T^j\left(t\right)$

即： ${\mathrm{\ρ}}^j\left(t\right)=\sqrt{{(x\left(t\right)-x^j\left(t\right))}^2+{(y\left(t\right)-y^j\left(t\right))}^2+{(z\left(t\right)-z^j\left(t\right))}^2}+\mathrm{c\δ}\left(t\right)$

由于上述方程是非线性的，需运用泰勒一阶展开来线性化。

设用户接收机的近似位置的位置坐标为

真实位置和近似位置之间的偏离为(δx(t)，δy(t)，δz(t))^T，将R^j(t)按泰勒级数绕近似位置展开，将位置偏离(δx(t)，δy(t)，δz(t))^T表示为已知卫星坐标和伪距测量值的线性函数：

${\mathrm{\ρ}}^j\left(t\right)={\hat{R}}^j\left(t\right)-\frac{x^j\left(t\right)-\hat{x}\left(t\right)}{{\hat{R}}^j\left(t\right)}\mathrm{\δx}\left(t\right)-\frac{y^j\left(t\right)-\hat{y}\left(t\right)}{{\hat{R}}^j\left(t\right)}\mathrm{\δy}\left(t\right)-\frac{z^j\left(t\right)-\hat{z}\left(t\right)}{{\hat{R}}^j\left(t\right)}\mathrm{\δz}\left(t\right)+\mathrm{c\δt}$

即： ${\hat{R}}^j\left(t\right)-{\mathrm{\ρ}}^j\left(t\right)=l^j\left(t\right)\mathrm{\δx}\left(t\right)+m^j\left(t\right)\mathrm{\δy}\left(t\right)+n^j\left(t\right)\mathrm{\δz}\left(t\right)-\mathrm{c\δt}$

其中， ${\hat{R}}^j\left(t\right)=\sqrt{{(x^j\left(t\right)-\hat{x}\left(t\right))}^2+{(y^j\left(t\right)-\hat{y}\left(t\right))}^2+{(z^j\left(t\right)-\hat{z}\left(t\right))}^2}$ 为接收机近似位置到卫星S^j的几何距离；

$l^j\left(t\right)=\frac{x^j\left(t\right)-\hat{x}\left(t\right)}{{\hat{R}}^j\left(t\right)},$ $m^j\left(t\right)=\frac{y^j\left(t\right)-\hat{y}\left(t\right)}{{\hat{R}}^j\left(t\right)},$ $n^j\left(t\right)=\frac{z^j\left(t\right)-\hat{z}\left(t\right)}{{\hat{R}}^j\left(t\right)}$ 表示由近似用户位置指向卫星S^j的单位矢量的方向余弦；

当有效卫星数目满足系统定位要求时，可以联立线性方程生成矩阵，然后用最小二乘法迭代求解。

但是，用一阶泰勒展开不能完全表征信号的真实特性，在理论上导致明显的信号误差，以至于导致使用该方法进行定位时存在定位精度不高、恶劣定位条件下鲁棒性不强的不足。因此，为满足对定位精度要求较高的用户，考虑对线性迭代求解方法的结果进行滤波。

滤波的目的是估计当前载体的运动状态，对于导航接收机来说就是估计天线相位中心的运动状态，包括其在导航解算坐标系内的位置、速度等。目前，最基本的数据滤波方法包括线性回归滤波，两点外推滤波、为纳滤波、加权最小二乘滤波、α-β与α-β-γ滤波、卡尔曼滤波和简化的卡尔曼滤波等。从跟踪精度来看，以卡尔曼滤波为最好，两点外推滤波为最差；从计算量来看，以卡尔曼滤波为最大，两点外推滤波最小。虽然，随着现代微处理技术突飞猛进的发展，计算要求已不再成为应用卡尔曼滤波的主要障碍，但在多星座组合定位应用中，随着观测向量或状态向量维数的扩充，卡尔曼滤波的计算量将会有大量增长。因此，卡尔曼滤波的使用，需根据用户对跟踪精度和计算量的要求进行折衷考虑。

常规卡尔曼滤波需要准确地知道系统噪声和观测噪声的统计特性，而实际应用中观测噪声的随机特性取决于测量设备以及系统的工作环境，很难准确描述，通常采用自适应卡尔曼滤波解决这一问题。为合理的应用卡尔曼滤波技术进行动态定位，最重要的是要建立载体的运动模型。常速度模型和常加速度模型构造简单易于实现，但适用匀速或匀加速等无明显机动的场合；一阶时间相关模型对一般机动的载体具有很好地描述效果，但如突发剧烈机动效果迅速恶化；机动目标“当前”统计模型本质上是非零均值时间相关模型，其机动加速度的“当前”概率密度用修正的瑞利分布描述，均值为“当前”加速度预测值，相比之下能更为真实地反映目标机动范围和强度的变化。

常速度模型：

$\left[\begin{array}{c}\stackrel{\·}{x}\\ \stackrel{\·\·}{x}\end{array}\right]=\left[\begin{array}{cc}0& 1\\ 0& 0\end{array}\right]\left[\begin{array}{c}x\\ \stackrel{\·}{x}\end{array}\right]+\left[\begin{array}{c}0\\ 1\end{array}\right]\mathrm{\ω}\left(t\right)$

常加速度模型：

$\left[\begin{array}{c}\stackrel{\·}{x}\\ \stackrel{\·\·}{x}\\ \stackrel{\·\·\·}{x}\end{array}\right]=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}x\\ \stackrel{\·}{x}\\ \stackrel{\·\·}{x}\end{array}\right]+\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]\mathrm{\ω}\left(t\right)$

一阶时间相关模型：

$\left[\begin{array}{c}\stackrel{\·}{x}\\ \stackrel{\·\·}{x}\\ \stackrel{\·\·\·}{x}\end{array}\right]=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 0& 0& -\mathrm{\α}\end{array}\right]\left[\begin{array}{c}x\\ \stackrel{\·}{x}\\ \stackrel{\·\·}{x}\end{array}\right]+\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]\mathrm{\ω}\left(t\right)$

机动目标“当前”统计模型：

$\left[\begin{array}{c}\stackrel{\·}{x}\\ \stackrel{\·\·}{x}\\ \stackrel{\·\·\·}{x}\end{array}\right]=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 0& 0& -\mathrm{\α}\end{array}\right]\left[\begin{array}{c}x\\ \stackrel{\·}{x}\\ \stackrel{\·\·}{x}\end{array}\right]+\left[\begin{array}{c}0\\ 0\\ \mathrm{\α}\end{array}\right]\stackrel{\‾}{a}\left(t\right)+\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]\mathrm{\ω}\left(t\right)$

其中，x，

分别为运动载体的位置、速度和加速度分量；ω(t)是均值为零，方差为σ²的高斯白噪声；α为机动时间常数的倒数，即机动频率，通常其取值范围为：转弯机动α＝1/60，逃避机动α＝1/20，大气扰动α＝1，确切值只有通过实验测量确定；a(t)为机动加速度“当前”均值，在每一采用周期内为常数。

发明内容

本发明的目的是提供一种包含BD2的多卫星导航系统兼容导航接收机定位系统及其定位方法，可以选择用于进行定位的卫星导航系统，单BD2定位可以在非常时期摆脱对其他卫星导航系统的依赖，多卫星导航系统的组合定位可以增加有效卫星的数量、保证卫星定位服务的连续性和可用性、增强观测卫星的几何图形强度、提高整个卫星定位系统的可靠性。

本发明提供的一种兼容导航接收机定位系统，包括参数配置模块、数据预处理器、卫星参数解算模块、伪距修正模块、坐标转换模块、定位解算模块、卡尔曼滤波器和报文形成模块。其中，参数配置模块的第一个输出端与数据预处理器的第一个输入端相连，参数配置模块的第二个输出端与卫星参数解算模块的第一个输入端相连，参数配置模块的第三个输出端与伪距修正模块的第一个输入端相连，参数配置模块的第四个输出端与坐标转换模块的第一个输入端相连，参数配置模块的第五个输出端与定位解算模块的第一个输入端相连，参数配置模块的第六个输出端与卡尔曼滤波器的第一个输入端相连，参数配置模块的第七个输出端与报文形成模块的第一个输入端相连；导航电文从数据预处理器的第二个输入端输入，观测数据从数据预处理器的第三个输入端输入，数据预处理器的第一个输出端与卫星参数解算模块的第二个输入端相连，数据预处理器的第二个输出端与伪距修正模块的第二个输入端相连；卫星参数解算模块的输出端与坐标转换模块的第二个输入端相连，坐标转换模块的输出端与定位解算模块的第二个输入端相连；伪距修正模块的输出端与定位解算模块的第三个输入端相连；定位解算模块的输出端与卡尔曼滤波器的第二个输入端相连；卡尔曼滤波器的第一个输出端与伪距修正模块的第三个输入端相连，卡尔曼滤波器的第二个输出端与报文形成模块的第二个输入端相连。其中，参数配置模块读取兼容导航接收机的配置参数，并将配置参数输出给定位系统的各模块以控制其工作状态；数据预处理器采集有效卫星的导航电文和观测数据帧，并从中提取星历数据、历书数据和观测信息；卫星参数解算模块获取有效卫星的卫星位置和卫星速度；伪距修正模块获取有效卫星的观测伪距和伪距误差修正信息，并对观测伪距进行修正；坐标转换模块将来自不同卫星导航系统的有效卫星的坐标统一到同一参考坐标系中；定位解算模块获取接收机位置和接收机速度信息；卡尔曼滤波器对接收机位置和接收机速度信息进行滤波；报文形成模块形成导航报文。

本发明提出的一种兼容导航接收机定位方法，具体包含以下步骤：

步骤一：读取接收机的配置参数；

参数配置模块读取接收机的配置参数，

首先，根据用户需求，确定接收机的工作状态和使用的卫星导航系统；工作状态包括：单卫星导航系统定位、双卫星导航系统组合定位、三卫星导航系统组合定位和四卫星导航系统组合定位，使用的卫星导航系统包括：全球定位系统GPS、全球导航卫星系统GLONASS、北斗二代导航系统BD2和“伽利略”卫星定位系统Galileo；

然后，根据用户需求，确定卡尔曼滤波器的工作状态，确定导航信息的数据更新频率，确定导航报文的格式；

最后，将接收机的配置参数输出给数据预处理器、卫星参数解算模块、伪距修正模块、坐标转换模块、定位解算模块、卡尔曼滤波器和报文形成模块，以控制各模块的工作状态；

其中，接收机的配置参数包括：接收机的工作状态、使用的卫星导航系统、卡尔曼滤波器的工作状态、导航信息的数据更新频率以及导航报文的格式；

步骤二：提取星历数据、历书数据和观测数据；

根据参数配置模块输出的接收机的工作状态和使用的卫星导航系统，数据预处理器依据有效卫星所属的卫星导航系统的说明书对采集的有效卫星的导航电文进行拼接、校验、纠错，从中提取有效卫星的星历数据和历书数据；

根据参数解算模块输出的接收机工作状态和使用的卫星导航系统，数据预处理器依据捕获跟踪系统提供的观测数据帧格式从采集的有效卫星观测数据帧中提取毫秒数、码片偏移、相位偏移和多普勒频率；

步骤三：获取有效卫星的参数信息；

首先，根据参数配置模块输出的接收机的工作状态和使用的卫星导航系统，卫星参数解算模块依据有效卫星所属的卫星导航系统的说明书利用有效卫星的星历数据获取其位置坐标和速度坐标；

然后，根据参数配置模块输出的接收机的工作状态和使用的卫星导航系统，伪距修正模块依据有效卫星所属卫星导航系统的说明书利用有效卫星的星历数据和历书数据获取其星钟修正误差、电离层修正误差和对流层修正误差；

之后，伪距修正模块通过将有效卫星的周内秒计数、比特数、毫秒数、码片偏移、相位偏移进行累加获取其观测伪距；

最后，伪距修正模块利用获取的有效卫星的星钟修正误差、电离层修正误差和对流层修正误差对该卫星的观测伪距进行误差修正，获取修正后的观测伪距；

步骤四：获取接收机位置和接收机速度；

首先，根据参数配置模块输出的接收机的工作状态和使用的卫星导航系统，在单卫星导航系统定位时坐标转换模块不参与系统工作，在多卫星导航系统组合定位时坐标转换模块依据各卫星导航系统参考坐标系的关系利用7参数Bursa模型将来自不同卫星导航系统的有效卫星的坐标转换到同一参考坐标系，获取各有效卫星在同一参考坐标系中的三维坐标；

然后，根据参数配置模块输出的接收机的工作状态和使用的卫星导航系统，定位解算模块利用有效卫星的参数信息获取接收机位置和接收机速度；

利用有效卫星的位置坐标和修正后的观测伪距通过线性迭代方法获取接收机的位置，利用有效卫星的速度坐标和多普勒频率通过线性迭代方法获取接收机的速度；

最后，根据参数配置模块输出的接收机的工作状态、使用的卫星导航系统和卡尔曼滤波器的工作状态，确定卡尔曼滤波器是否对定位解算模块输出的接收机位置和接收机速度信息进行滤波处理；

如果卡尔曼滤波器的工作状态是工作，基于机动目标“当前”统计模型建立系统模型和测量模型，并通过自适应卡尔曼滤波获取接收机位置和接收机速度的滤波结果；

步骤五：将导航信息生成报文；

根据参数配置模块输出的用户对导航报文的要求，报文形成模块将导航信息以用户需要的格式形成报文。

本发明的优点在于：

(1)本发明采用了多卫星导航系统兼容定位技术，单BD2独立定位可以在非常时期摆脱对其他导航系统的依赖；

(2)本发明采用了多卫星导航系统兼容定位技术，多卫星导航系统组合定位可以大大增加有效卫星的数量，提高整个卫星定位系统的完善性监测性能，保证卫星定位服务的连续性、可用性，提高卫星导航定位应用的安全性。

(3)本发明采用了多卫星导航系统兼容定位技术，多卫星导航系统组合定位可以增强观测卫星的几何图形强度，提高整个卫星定位系统的可靠性。

附图说明

图1是本发明的系统框图；

图2是本发明的方法流程图；

图3是本发明组合定位与单GPS定位仰角限为10度时的有效卫星数目曲线图；

图4是本发明组合定位与单GPS定位仰角限为30度时的有效卫星数目曲线图；

图5是本发明组合定位与单GPS定位仰角限为10度时的DOP曲线图；

图6是本发明组合定位与单GPS定位仰角限为30度时的DOP曲线图；

图7是本发明组合定位与单GPS定位时的水平定位误差曲线图；

图8是本发明组合定位与单GPS定位时的垂直定位误差曲线图。

图中：1.参数配置模块  2.数据预处理器  3.卫星参数解算模块  4.伪距修正模块5.坐标转换模块  6.定位解算模块  7.卡尔曼滤波器  8.报文形成模块

具体实施方式

下面将以选择BD2和GPS进行双卫星导航系统组合定位为例，结合附图和实施例对本发明作进一步的详细说明。

本发明是一种兼容导航接收机定位系统及其方法，系统框图如图1所示，它包括参数配置模块1、数据预处理器2、卫星参数解算模块3、伪距修正模块4、坐标转换模块5、定位解算模块6、卡尔曼滤波器7和报文形成模块8。其中，参数配置模块1的第一个输出端与数据预处理器2的第一个输入端相连，参数配置模块1的第二个输出端与卫星参数解算模块3的第一个输入端相连，参数配置模块1的第三个输出端与伪距修正模块4的第一个输入端相连，参数配置模块1的第四个输出端与坐标转换模块5的第一个输入端相连，参数配置模块1的第五个输出端与定位解算模块6的第一个输入端相连，参数配置模块1的第六个输出端与卡尔曼滤波器7的第一个输入端相连，参数配置模块1的第七个输出端与报文形成模块8的第一个输入端相连；导航电文从数据预处理器2的第二个输入端输入，观测数据从数据预处理器2的第三个输入端输入，数据预处理器2的第一个输出端与卫星参数解算模块3的第二个输入端相连，数据预处理器2的第二个输出端与伪距修正模块4的第二个输入端相连；卫星参数解算模块3的输出端与坐标转换模块5的第二个输入端相连，坐标转换模块5的输出端与定位解算模块6的第二个输入端相连；伪距修正模块4的输出端与定位解算模块6的第三个输入端相连；定位解算模块6的输出端与卡尔曼滤波器7的第二个输入端相连；卡尔曼滤波器7的第一个输出端与伪距修正模块4的第三个输入端相连，卡尔曼滤波器7的第二个输出端与报文形成模块8的第二个输入端相连；

其中，参数配置模块1读取兼容导航接收机的配置参数，并将配置参数输出给定位系统的各模块以控制其工作状态；

数据预处理器2采集有效卫星的导航电文和观测数据帧，并从中提取星历数据、历书数据和观测信息；

卫星参数解算模块3获取有效卫星的卫星位置和卫星速度；

伪距修正模块4获取有效卫星的观测伪距和伪距误差修正信息，并对观测伪距进行修正；

坐标转换模块5将来自不同卫星导航系统的有效卫星的坐标统一到同一参考坐标系中；

定位解算模块6获取接收机位置和接收机速度信息；

卡尔曼滤波7对接收机的位置和速度信息进行滤波；

报文形成模块8形成导航报文。

本发明所述兼容导航接收机定位方法流程图如图2所示。

接收机完成了捕获跟踪，获得了一个完整的导航电文子帧并且观测数据已经更新后，需要根据用户的需要输出导航信息。本发明所述兼容导航接收机定位方法的具体实施方式如下：

步骤一：读取接收机的配置参数；

参数配置模块1读取接收机的配置参数，

首先，根据用户需求，确定接收机的工作状态和使用的卫星导航系统；工作状态为：双卫星导航系统组合定位，使用的卫星导航系统包括：北斗二代导航系统BD2和全球定位系统GPS；

然后，根据用户需求，确定卡尔曼滤波器7的工作状态为工作，确定导航信息的数据更新频率为10Hz，确定导航报文的格式为RINEX格式；

最后，将接收机的配置参数输出给数据预处理器2、卫星参数解算模块3、伪距修正模块4、坐标转换模块5、定位解算模块6、卡尔曼滤波器7和报文形成模块8，以控制各模块的工作状态；

其中，接收机的配置参数包括：接收机的工作状态、使用的卫星导航系统、卡尔曼滤波器7的工作状态、导航信息的数据更新频率以及导航报文的格式；

步骤二：提取星历数据、历书数据和观测数据；

根据参数配置模块1输出的接收机的工作状态和使用的卫星导航系统，接收机使用BD2和GPS进行双卫星导航系统组合定位，数据预处理器1依据BD2说明书，对采集的BD2有效卫星的导航电文子帧进行电文拼接和纠错，并从纠错后的电文中提取星历数据和历书数据；数据预处理器1依据GPS说明书，对采集的GPS有效卫星的导航电文子帧进行电文拼接和校验，并从校验后的电文中提取星历数据和历书数据；

其中，星历数据包括轨道长半轴方根

轨道偏心率e，参考时刻轨道倾角i₀，近地点角距ω，参考时刻升交点赤经Ω₀，参考时刻的平近点角M₀，平均运行速度差Δn，轨道倾角变率

升交点赤经变率

升交距角的调和改正项中的正弦项振幅C_us，升交距角的调和改正项中的余弦项振幅C_uc，卫星地心距的调和改正项中的正弦项振幅C_rs，卫星地心距的调和改正项中的余弦项振幅C_rc，轨道倾角的调和改正项中的正弦项振幅C_is，轨道倾角的调和改正项中的余弦项振幅C_ic，星历参数的参考历元t_oe，星历数据的龄期AODE；

其总，历书数据包括时钟数据基准时间t_oc，卫星时间修正参数a₀，卫星时间修正参数a₁、卫星时间修正参数a₂，卫星钟数据有效龄期AODC，延迟改正参数(GPS单频接收机延迟改正参数T_GD；BD2星上设备时延差T_GD1，BD2星上设备时延差T_GD2，T_GD1为BD2卫星发射的B1频点与B3频点信号之间存在的设备时延之差，T_GD2为BD2卫星发射的B2频点与B3频点信号之间存在的设备时延之差；电离层修正参数GPS：α_n、β_n，其中n＝1，2，3；BD2D1码：α_n、β_n，其中n＝1，2，3，D2码：A₁、B、α_n、β_n、γ_n，其中n＝1，2，3；

根据参数配置模块1输出的接收机的工作状态和使用的卫星导航系统，接收机使用BD2和GPS进行双卫星导航系统组合定位，数据预处理器2依据捕获跟踪系统提供的观测数据帧格式，从采集的观测数据帧中提取BD2和GPS双系统有效卫星的毫秒数、码片偏移和相位偏移以及观测信息多普勒频率；

步骤三：获取有效卫星的参数信息；

首先，根据参数配置模块1输出的接收机的工作状态和使用的卫星导航系统，接收机使用BD2和GPS进行双卫星导航系统组合定位，卫星参数解算模块3分别依据BD2说明书和GPS说明书，利用双系统有效卫星的星历数据获取其位置坐标和速度坐标；

然后，根据参数配置模块1输出的接收机的工作状态和使用的卫星导航系统，接收机使用BD2和GPS进行双卫星导航系统组合定位，伪距修正模块4分别依据BD2说明书和GPS说明书，利用双系统有效卫星的星历数据和历书数据获取其星钟修正误差、电离层修正误差和对流层修正误差；

之后，伪距修正模块4通过将双系统有效卫星的周内秒计数、比特数、毫秒数、码片偏移、相位偏移进行累加获取其观测伪距；

最后，伪距修正模块4利用获取的双系统有效卫星的星钟修正误差、电离层修正误差和对流层修正误差对其观测伪距进行误差修正，获取修正后的观测伪距；

步骤四：获取接收机位置和接收机速度；

首先，根据参数配置模块1输出的接收机的工作状态和使用的卫星导航系统，接收机使用BD2和GPS进行双卫星导航系统组合定位，BD2参考坐标系CGS2000与GPS参考坐标系WGS-84间的差异很小，在10cm以内，认为是等同的。故使用GPS与BD2双系统进行定位时，坐标转换模块5不需要参与工作，直接将卫星参数解算模块3输出的有效卫星坐标输入给定位解算模块6即可；

然后，根据参数配置模块1输出的接收机的工作状态和使用的卫星导航系统，接收机使用BD2和GPS进行双卫星导航系统组合定位，定位解算模块6利用BD2和GPS双系统的有效卫星的参数信息获取接收机位置和接收机速度；

A、获取接收机位置

大气折射对伪距观测量的影响已改正，且卫星钟差可由导航电文给出的参数加以修正，则站星伪距观测为：

BD2和GPS双系统工作模式下， ${\mathrm{\ρ}}_u^j=R_u^j+\mathrm{c\Δ}{t}_1+c{\mathrm{\Δt}}_2,$

其中，j＝1，2，…，n，j为观测卫星S^j的编号，n为观测卫星总数；ρ_u ^j为卫星S_j的观测伪距进行星历误差、电离层误差和对流层误差修正后的伪距，R_u ^j为接收机与卫星S_j间的几何距离，Δt₁为接收机时钟与BD2系统时钟之间的钟差，Δt₂为接收机时钟与GPS系统时钟之间的钟差，c为光速。

在一段时间内，接收机连续同步观测不少于5颗卫星时，利用有效卫星的位置坐标和修正后的观测伪距获取接收机的位置坐标偏移和钟差模型：

$\left[\begin{array}{c}\mathrm{\Δ}{x}_u\\ \mathrm{\Δ}{y}_u\\ \mathrm{\Δ}{z}_u\\ \mathrm{c\Δ}{t}_1\\ \mathrm{c\Δ}{t}_2\end{array}\right]=-{\left(H^{T}H\right)}^{-1}\left(H^T\mathrm{\Δ\ρ}\right)$

其中，上述位置坐标偏移和钟差模型等号右侧部分的 $H=\left[\begin{array}{ccccc}{a}_x^1& a_y^1& a_z^1& -1& 0\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ a_x^q& a_y^q& a_z^q& -1& 0\\ a_x^{q+1}& a_y^{q+1}& a_z^{q+1}& 0& -1\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ a_x^n& a_y^n& a_z^n& 0& -1\end{array}\right],$

为一个n×5的矩阵模块， $\mathrm{\Δ\ρ}=\left[\begin{array}{c}{\mathrm{\ρ}}_u^1\\ .\\ .\\ .\\ {\mathrm{\ρ}}_u^q\\ {\mathrm{\ρ}}_u^{q+1}\\ .\\ .\\ .\\ {\mathrm{\ρ}}_u^n\end{array}\right]-\left[\begin{array}{c}{\hat{R}}_u^1\\ .\\ .\\ .\\ {\hat{R}}_u^q\\ {\hat{R}}_u^{q+1}\\ .\\ .\\ .\\ {\hat{R}}_u^n\end{array}\right]$ 为一个n×1的矩阵模块，

其中，j为观测卫星S^j的编号，j＝1，2，…，n，n为观测卫星总数，S¹～S^q为BD2的卫星，S^q+1～Sⁿ为GPS的卫星；

(Δx_u，Δy_u，Δz_u)^T为用户接收机的位置坐标偏移；

Δt₁为用户接收机时钟与BD2系统时钟之间的钟差，Δt₂为用户接收机时钟与GPS系统时钟之间的钟差；

$a_x^j=\frac{x^j-{\hat{x}}_u}{{\hat{R}}_u^j},$ $a_y^j=\frac{y^j-{\hat{y}}_u}{{\hat{R}}_u^j},$ $a_z^j=\frac{z^j-{\hat{z}}_u}{{\hat{R}}_u^j}$ 各项分别表示从用户接收机近似位置指向卫星S^j的单位向量的方向余弦，j＝1，…，q，q+1，…，n，a_x ¹，a_y ¹，a_z ¹分别表示从用户接收机近似位置指向BD2的卫星S¹的单位向量的方向余弦，a_x ^q，a_y ^q，a_z ^q分别表示从用户接收机近似位置指向BD2的卫星S^q的单位向量的方向余弦，a_x ^q+1，a_y ^q+1，a_z ^q+1分别表示从用户接收机近似位置指向GPS的卫星S^q+1的单位向量的方向余弦，a_x ⁿ，a_y ⁿ，a_z ⁿ分别表示从用户接收机近似位置指向GPS的卫星Sⁿ的单位向量的方向余弦；

为用户接收机近似位置的位置坐标；

(x^j，y^j，z^j)^T为卫星S^j的位置坐标，j＝1，…q，q+1，…，n；

为用户接收机的近似位置与卫星S^j间的几何距离，j＝1，…q，q+1，…，n，

为用户接收机的近似位置与BD2的卫星S¹间的几何距离，

为用户接收机的近似位置与BD2的卫星S^q间的几何距离，为用户接收机的近似位置与GPS的卫星S^q+1间的几何距离，

为用户接收机的近似位置与GPS的卫星Sⁿ间的几何距离；

ρ_u ^j为对卫星S^j的观测伪距进行星钟误差、电离层误差和对流层误差修正后的伪距，j＝1，…q，q+1，…，n，ρ_u ¹为对BD2的卫星S¹的观测伪距进行星钟误差、电离层误差和对流层误差修正后的伪距，ρ_u ^q为对BD2的卫星S^q的观测伪距进行星钟误差、电离层误差和对流层误差修正后的伪距，ρ_u ^q+1为对GPS的卫星S^q+1的观测伪距进行星钟误差、电离层误差和对流层误差修正后的伪距，ρ_u ⁿ为对GPS的卫星Sⁿ的观测伪距进行星钟误差、电离层误差和对流层误差修正后的伪距；

c为光速；

用户接收机的位置坐标为：

$\left[\begin{array}{c}{x}_u\\ y_u\\ z_u\end{array}\right]=\left[\begin{array}{c}{\hat{x}}_u\\ {\hat{y}}_u\\ {\hat{z}}_u\end{array}\right]+\left[\begin{array}{c}\mathrm{\Δ}{x}_u\\ \mathrm{\Δ}{y}_u\\ \mathrm{\Δ}{z}_u\end{array}\right];$

B、获取接收机速度

由于用户接收机载体和卫星之间的相对运动，接收机接收到的载波信号与卫星发射的载波信号频率不同，产生多普勒频移。频移的大小与接收机与卫星之间距离的变率有关。改正大气折射对伪距观测量的影响，并由导航电文给出的参数修正卫星钟差，则可得站星观测伪距率为：

BD2和GPS双系统工作模式下， $d_u^j=\left(\begin{array}{ccc}{a}_x^j& a_y^j& a_z^j\end{array}\right)[\left[\begin{array}{c}{v}_x^j\\ v_y^j\\ v_z^j\end{array}\right]-\left[\begin{array}{c}{v}_{\mathrm{ux}}\\ v_{\mathrm{uy}}\\ v_{\mathrm{uz}}\end{array}\right]]+c{\stackrel{\·}{t}}_1+c{\stackrel{\·}{t}}_2,$

其中，j为观测卫星S^j的编号，j＝1，2，…，n，n为观测卫星总数，d_u ^j为卫星S_j的观测伪距率，a_x ^j，a_y ^j，a_z ^j各项分别表示从用户接收机位置指向卫星S^j的单位向量的方向余弦，(v_x ^j，v_y ^j，v_z ^j)^T为卫星S^j的速度坐标，(v_ux，v_uy，v_uz)^T为用户接收机的速度坐标，

为接收机时钟与BD2系统时钟之间的钟漂，

为接收机时钟与GPS系统时钟之间的钟漂，c为光速。

在一段时间内，接收机连续同步观测不少于5颗卫星时，利用接收机的位置坐标和有效卫星的多普勒频率获取接收机的速度坐标和钟漂模型：

$\left[\begin{array}{c}{v}_{\mathrm{ux}}\\ v_{\mathrm{uy}}\\ v_{\mathrm{uz}}\\ c{\stackrel{\·}{t}}_1\\ c{\stackrel{\·}{t}}_2\end{array}\right]=-{\left(H^{T}H\right)}^{-1}\left(H^{T}d\right),$

其中，上述速度坐标和钟漂模型等号右侧部分的 $H=\left[\begin{array}{ccccc}{a}_x^1& a_y^1& a_z^1& -1& 0\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ a_x^q& a_y^q& a_z^q& -1& 0\\ a_x^{q+1}& a_y^{q+1}& a_z^{q+1}& 0& -1\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ a_x^n& a_y^n& a_z^n& 0& -1\end{array}\right]$ 为一个n×5的矩阵模块，d＝[d¹…d^q d^q+1…dⁿ]^T为一个n×1的矩阵模块，

其中 $d^j=d_u^j-\left[\begin{array}{ccc}{a}_x^j& a_y^j& a_z^j\end{array}\right]\left[\begin{array}{c}{v}_x^j\\ v_y^j\\ v_z^j\end{array}\right],$

j为观测卫星S^j的编号，j＝1，…q，q+1，…，n，n为观测卫星总数，S¹～S^q为BD2的卫星，S^q+1～Sⁿ为GPS的卫星；

(v_ux，v_uy，v_uz)^T为用户接收机的速度坐标；

为用户接收机时钟与BD2系统时钟之间的钟漂，

为用户接收机时钟与GPS系统时钟之间的钟漂；

$a_x^j=\frac{x^j-x_u}{R_u^j},$ $a_y^j=\frac{y^j-y_u}{R_u^j},$ $a_z^j=\frac{z^j-z_u}{R_u^j}$ 各项分别表示从用户接收机指向卫星S^j的单位向量的方向余弦，j＝1，…q，q+1，…，n，a_x ¹，a_y ¹，a_z ¹分别表示从用户接收机指向BD2GPS卫星S¹的单位向量的方向余弦，a_x ^q，a_y ^q，a_z ^q分别表示从用户接收机指向BD2中的卫星S^q的单位向量的方向余弦，a_x ^q+1，a_y ^q+1，a_z ^q+1分别表示从用户接收机指向GPS的卫星S^q+1的单位向量的方向余弦，a_x ⁿ，a_y ⁿ，a_z ⁿ分别表示从用户接收机指向GPS的卫星Sⁿ的单位向量的方向余弦；

(x_u，y_u，z_u)^T为用户接收机的位置坐标；

(x^j，y^j，z^j)^T为卫星S^j的位置坐标，j＝1，…q，q+1，…，n；

R_u ^j为用户接收机与卫星S^j间的几何距离，j＝1，…q，q+1，…，n；

d_u ^j为卫星S^j的观测伪距率，j＝1，…q，q+1，…，n，d_u ¹为BD2的卫星S¹的观测伪距率，d_u ^q为BD2的卫星S^q的观测伪距率，j＝q+1时d_u ^q+1为GPS的卫星S^q+1的观测伪距率，j＝n时d_u ⁿ为GPS的卫星Sⁿ的观测伪距率；

(v_x ^j，v_y ^j，v_z ^j)^T为卫星S^j的速度坐标，j＝1，…q，q+1，…，n；

c为光速；

(3)根据参数配置模块1输出的接收机工作状态和卡尔曼滤波器7的工作状态，卡尔曼滤波器7基于机动目标“当前”统计模型建立系统模型和测量模型，并通过自适应卡尔曼滤波获取接收机位置和接收机速度的滤波结果；

采用递推方法给出滤波结果：

状态估计模型

估计均方误差模型P(k)＝[I-K(k)H(k)]P(k|k-1)；

其中，状态一步预测模型 $\hat{X}\left(k\right|k-1)={\mathrm{\Φ}}_1\left(k\right|k-1)\hat{X}(k-1),$

一步预测均方误差模型P(k|k-1)＝Φ(k|k-1)P(k-1)Φ^T(k|k-1)+Q(k-1)，

滤波增益模型K(k)＝P(k|k-1)H^T(k)[H(k)P(k|k-1)H^T(k)+R(k)]^-1；

其中，系统状态模型 $\hat{X}(k-1)=(x_u,{\stackrel{\·}{x}}_u,{\stackrel{\·\·}{x}}_u,y_u,{\stackrel{\·}{y}}_u,{\stackrel{\·\·}{y}}_u,z_u,{\stackrel{\·}{z}}_u,{\stackrel{\·\·}{z}}_u,\mathrm{c\Δ}{t}_1,c{\stackrel{\·}{t}}_1,\mathrm{c\Δ}{t}_2,c{\stackrel{\·}{t}}_2),$

分别为接收机在参考坐标系x轴方向的位置、速度和加速度，

分别为接收机在参考坐标系y轴方向的位置、速度和加速度，

分别为接收机在参考坐标系z轴方向的位置、速度和加速度，Δt₁，

分别为接收机与BD2系统时钟之间的钟差和钟漂，Δt₂，

分别为接收机与GPS系统时钟之间的钟差和钟漂，其初值本发明取为X₀＝[x_u，v_ux，0，y_u，v_uy，0，z_u，v_uz，0，0，0，0，0]^T，(x_u，y_u，z_u)^T，(v_ux，v_uy，v_uz)^T分别为定位解算模块6输出的接收机位置坐标和速度坐标；

其中，第一系统状态转移模型

${\mathrm{\Φ}}_1\left(k\right|k-1)=\mathrm{diag}[{\mathrm{\φ}}_{1x}\left(k\right|k-1),{\mathrm{\φ}}_{1y}\left(k\right|k-1),{\mathrm{\φ}}_{1z}\left(k\right|k-1),{\mathrm{\φ}}_{1t_1}\left(k\right|k-1),{\mathrm{\φ}}_{1t_2}(k|k-1)],$

${\mathrm{\φ}}_{1x}\left(k\right|k-1)={\mathrm{\φ}}_{1y}\left(k\right|k-1)={\mathrm{\φ}}_{1z}\left(k\right|k-1)=\left[\begin{array}{ccc}1& T& T^2/2\\ 0& 1& T\\ 0& 0& 1\end{array}\right],$

${\mathrm{\φ}}_{1t_1}\left(k\right|k-1)={\mathrm{\φ}}_{1t_2}\left(k\right|k-1)=\left[\begin{array}{cc}1& T\\ 0& 1\end{array}\right],$

T为采样周期，本发明取为T＝0.1s，

P(k-1)为状态协方差，其初值本发明取为

其中，第二系统状态转移模型 $\mathrm{\Φ}\left(k\right|k-1)=\mathrm{diag}[{\mathrm{\φ}}_x\left(k\right|k-1),{\mathrm{\φ}}_y\left(k\right|k-1),{\mathrm{\φ}}_z\left(k\right|k-1),{\mathrm{\φ}}_{t_1}\left(k\right|k-1),{\mathrm{\φ}}_{t_2}(k|k-1)],$

${\mathrm{\φ}}_x\left(k\right|k-1)=\left[\begin{array}{ccc}1& T& {\mathrm{\α}}_x^{-2}(-1+{\mathrm{\α}}_{x}T+e^{-{\mathrm{\α}}_{x}T})\\ 0& 1& (1-e^{-{\mathrm{\α}}_{x}T}){\mathrm{\α}}_x^{-1}\\ 0& 0& e^{-{\mathrm{\α}}_{x}T}\end{array}\right],$

${\mathrm{\φ}}_y\left(k\right|k-1)=\left[\begin{array}{ccc}1& T& {\mathrm{\α}}_y^{-2}(-1+{\mathrm{\α}}_{y}T+e^{-{\mathrm{\α}}_{y}T})\\ 0& 1& (1-e^{-{\mathrm{\α}}_{y}T}){\mathrm{\α}}_y^{-1}\\ 0& 0& e^{-{\mathrm{\α}}_{y}T}\end{array}\right],$

${\mathrm{\φ}}_z\left(k\right|k-1)=\left[\begin{array}{ccc}1& T& {\mathrm{\α}}_z^{-2}(-1+{\mathrm{\α}}_{z}T+e^{-{\mathrm{\α}}_{z}T})\\ 0& 1& (1-e^{-{\mathrm{\α}}_{z}T}){\mathrm{\α}}_z^{-1}\\ 0& 0& e^{-{\mathrm{\α}}_{z}T}\end{array}\right],$

${\mathrm{\φ}}_{1t_1}\left(k\right|k-1)={\mathrm{\φ}}_{1t_2}\left(k\right|k-1)=\left[\begin{array}{cc}1& T\\ 0& 1\end{array}\right],$

α_x，α_y，α_z为机动时间常数的倒数，本发明取为α_x＝α_y＝α_z＝0.2Hz；

其中，系统噪声方差模型

$Q_x\left(k\right)=\left[\begin{array}{ccc}{q}_{x11}& q_{x12}& q_{x13}\\ q_{x12}& q_{x22}& q_{x23}\\ q_{x13}& q_{x23}& q_{x33}\end{array}\right],$ $Q_{t_1}\left(k\right)=Q_{t_2}\left(k\right)=\left[\begin{array}{cc}{\mathrm{\σ}}_{\mathrm{\Δt}}^2& 0\\ 0& {\mathrm{\σ}}_{\stackrel{\·}{t}}^2\end{array}\right],$

$q_{x11}=0.5{\mathrm{\α}}_x^{-5}(1-e^{2{\mathrm{\α}}_{x}T}+2{\mathrm{\α}}_{x}T-4{\mathrm{\α}}_{x}T{e}^{-{\mathrm{\α}}_{x}T}-2{\mathrm{\α}}_x^2{T}^2+\frac{2}{3}{\mathrm{\α}}_x^3{T}^3),$

$q_{x12}=0.5{\mathrm{\α}}_x^{-4}(1-{2e}^{-{\mathrm{\α}}_{x}T}+e^{-2{\mathrm{\α}}_{x}T}-2{\mathrm{\α}}_{x}T+2{\mathrm{\α}}_{x}T{e}^{-{\mathrm{\α}}_{x}T}+{\mathrm{\α}}_x^2{T}^2),$

$q_{x13}=0.5{\mathrm{\α}}_x^{-3}(1-e^{2{\mathrm{\α}}_{x}T}-2{\mathrm{\α}}_{x}T{e}^{-{\mathrm{\α}}_{x}T}),$

$q_{x23}=0.5{\mathrm{\α}}_x^{-2}(1-2e^{-{\mathrm{\α}}_{x}T}+e^{-2{\mathrm{\α}}_{x}T}),$

$q_{x22}=0.5{\mathrm{\α}}_x^{-3}(-3+4e^{-{\mathrm{\α}}_{x}T}-e^{-2{\mathrm{\α}}_{x}T}+2{\mathrm{\α}}_{x}T),$

$q_{x33}=0.5{\mathrm{\α}}_x^{-1}(1-e^{-2{\mathrm{\α}}_{x}T}),$

Q_y(k)，Q_z(k)中的元素模型与Q_x(k)中的元素模型相同，将Q_x(k)中各元素模型中的α_x分别用α_y，α_z代替，即可得到Q_y(k)，Q_z(k)，

分别为与接收机在参考坐标系x，y，z坐标方向上的加速度分量对应的系统噪声方差，本发明取为 ${\mathrm{\σ}}_{a_x}^2={\mathrm{\σ}}_{a_y}^2={\mathrm{\σ}}_{a_z}^2={(7m/s^2)}^2,$ α_x，α_y，α_z为机动时间常数的倒数，本发明取为α_x＝α_y＝α_z＝0.2Hz，σ_Δt ²，

分别为与接收机钟差和钟漂对应的系统噪声方差，本发明取为 ${\mathrm{\σ}}_{\mathrm{\Δt}}^2={\left(0.3s\right)}^2,$ ${\mathrm{\σ}}_{\stackrel{\·}{t}}^2={(0.3s/s)}^2;$

其中，第一系统量测模型H(k)＝[H¹…H^j…Hⁿ]^T，

$H^j=\left[\begin{array}{ccccc}{A}_x^j& A_y^j& A_z^j& B_1& B_2\end{array}\right],$

$A_x^j=\left[\begin{array}{ccc}-a_x^j& 0& 0\\ 0& a_x^j& 0\end{array}\right],$ $A_y^j=\left[\begin{array}{ccc}-a_y^j& 0& 0\\ 0& a_y^j& 0\end{array}\right],$ $A_z^j=\left[\begin{array}{ccc}-a_z^j& 0& 0\\ 0& a_z^j& 0\end{array}\right],$

当卫星S^j属于BD2， $B_1=\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right],$ $B_2=\left[\begin{array}{cc}0& 0\\ 0& 0\end{array}\right],$

当卫星S^j属于GPS， $B_1=\left[\begin{array}{cc}0& 0\\ 0& 0\end{array}\right],$ $B_2=\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right],$

j为接收机定位时使用的有效卫星S^j的编号，j＝1，2，…，n，n为接收机定位使用的有效卫星的数目，

$a_x^j=\frac{x^j-x_u}{R_u^j},$ $a_y^j=\frac{y^j-y_u}{R_u^j},$ $a_z^j=\frac{z^j-z_u}{R_u^j}$ 各项分别表示从用户接收机指向卫星S^j的单位向量的方向余弦，(x_u，y_u，z_u)^T为用户接收机的位置坐标，(x^j，y^j，z^j)^T为卫星S^j的位置坐标，R_u ^j为用户接收机与卫星S^j间的几何距离，j＝1，2，…，n；

其中，量测噪声模型

分别为卫星S^j的伪距观测噪声方差和伪距率观测噪声方差，

分别为卫星S¹的伪距观测噪声方差和伪距率观测噪声方差，分别为卫星Sⁿ的伪距观测噪声方差和伪距率观测噪声方差，本发明取为 ${\mathrm{\σ}}_{{\mathrm{\ρ}}_u^j}^2={\left(12.5m\right)}^2,$ ${\mathrm{\σ}}_{d_u^j}^2={(0.8m/s)}^2,$ j＝1，2，…，n；

其中，观测模型 $Z\left(k\right)={\left[\begin{array}{cccccccc}{\mathrm{\ρ}}_u^1& d_u^1& \·\·\·& {\mathrm{\ρ}}_u^j& d_u^j& \·\·\·& {\mathrm{\ρ}}_u^n& d_u^n\end{array}\right]}^T,$ ρ_u ^j，d_u ^j分别为对卫星S^j的观测伪距进行星钟误差、电离层误差和对流层误差修正后的伪距和观测伪距率，j＝1，2，…，n，ρ_u ¹，d_u ¹分别为对卫星S¹的观测伪距进行星钟误差、电离层误差和对流层误差修正后的伪距和观测伪距率，ρ_u ⁿ，d_u ⁿ分别为对卫星Sⁿ的观测伪距进行星钟误差、电离层误差和对流层误差修正后的伪距和观测伪距率；

其中，第二系统量测模型h(X)＝[h¹…h^j…hⁿ]^T，

$h^j=\left[\begin{array}{c}\sqrt{{(x^j-x_u)}^2+{(y^j-y_u)}^2+{(z^j-z_u)}^2}+\underset{i=1}{\overset{2}{\mathrm{\Σ}}}c{k}_i^j\mathrm{\Δ}{t}_i\\ a_x^j(v_{\mathrm{ux}}-v_x^j)+a_y^j(v_{\mathrm{uy}}-v_y^j)+a_z^j(v_{\mathrm{uz}}-v_z^j)-\underset{i=1}{\overset{2}{\mathrm{\Σ}}}c{k}_i^j{\stackrel{\·}{t}}_i\end{array}\right],$

当卫星S^j属于BD2时， $k_1^j=1,$ $k_2^j=0,$

当卫星S^j属于GPS时， $k_1^j=0,$ $k_2^j=1;$

步骤五：将导航信息生成报文；

根据参数配置模块1输出的用户对导航报文的要求，报文形成模块8将数据更新频率为10Hz的导航信息以RINEX格式形成报文。

本发明的兼容导航接收机定位方法，单BD2系统独立工作可以在非常时期摆脱对其他卫星导航系统的依赖。多卫星导航系统组合定位可以大大增加有效卫星的数量，保证卫星定位服务的连续性、可用性。由于可见卫星数量大幅度增加，多卫星导航系统组合定位大大提高整个卫星定位系统的完善性监测性能，提高卫星导航定位应用的安全性。同时，多卫星导航系统组合定位可以增强观测卫星的几何图形强度，提高整个卫星定位系统的可靠性。

图3示出了BD2-GPS组合导航与单GPS定位仰角限为10度时的有效卫星数目曲线。图3中，实线表示应用本发明提供的兼容导航接收机定位方法进行组合定位仰角限为10度时的有效卫星数目，虚线表示应用单GPS进行定位仰角限为10度时的有效卫星数目。由图3可见，仰角限为10度时，单GPS全天有效卫星数目达5颗及以上，BD2-GPS组合系统全天有效卫星数目达13颗及以上，BD2-GPS组合星座为GPS星座有效卫星数目的两倍左右。

图4示出了BD2-GPS组合导航与单GPS定位仰角限为30度时的有效卫星数目曲线。图4中，实线表示应用本发明提供的兼容导航接收机定位方法进行组合定位仰角限为30度时的有效卫星数目，虚线表示应用单GPS进行定位仰角限为30度时的有效卫星数目。由图4可见，仰角限为30度时，单GPS全天有效卫星达4颗及以上的时刻占全天24h的90％，BD2-GPS组合系统全天有效卫星达7颗及以上，与单GPS相比，BD2-GPS组合系统把高度角的限制条件扩展了，即增强了卫星导航系统的可用性。

图5示出了BD2-GPS组合定位与单GPS定位仰角限为10度时的DOP曲线。图5中，实线表示应用本发明提供的兼容导航接收机定位方法进行组合定位仰角限为10度时的DOP，虚线表示应用单GPS进行定位仰角限为10度时的DOP。由图5可见，仰角限为10度时，单GPS的DOP值分布在1～6之间，波动较大，BD2-GPS组合系统DOP值明显减小，波动更加平缓。

图6示出了BD2-GPS组合定位与单GPS定位仰角限为30度时的DOP曲线。图6中，实线表示应用本发明提供的兼容导航接收机定位方法进行组合定位仰角限为30度时的DOP，虚线表示应用单GPS进行定位仰角限为30度时的DOP。由图6可见，仰角限为30度时，单GPS系统全天有10％时段观测不到4颗卫星，其他时段DOP较大，定位精度不能达到要求，BD2-GPS组合系统除部分超过4.0外，多介于2.0～4.0之间，波动平缓，精度明显优于单GPS，适合高精度的全天候导航定位。

图7示出了BD2-GPS组合定位与单GPS定位时的水平定位误差曲线。图中，实线表示应用本发明提供的兼容导航接收机定位方法进行组合定位时的水平定位误差，虚线表示应用单GPS进行定位时的水平定位误差。由图7可见，单GPS的水平定位精度低于BD2-GPS组合系统。

图8示出了BD2-GPS组合定位与单GPS定位时的垂直定位误差曲线。图中，实线表示应用本发明提供的兼容导航接收机定位方法进行组合定位时的垂直定位误差，虚线表示应用单GPS进行定位时的垂直定位误差。由图8可见，单GPS的垂直定位精度低于BD2-GPS组合系统。

Claims (3)

1.一种兼容导航接收机定位系统，其特征在于，该系统包括：参数配置模块(1)、数据预处理器(2)、卫星参数解算模块(3)、伪距修正模块(4)、坐标转换模块(5)、定位解算模块(6)、卡尔曼滤波器(7)和报文形成模块(8)；

其中，参数配置模块(1)的第一个输出端与数据预处理器(2)的第一个输入端相连，参数配置模块(1)的第二个输出端与卫星参数解算模块(3)的第一个输入端相连，参数配置模块(1)的第三个输出端与伪距修正模块(4)的第一个输入端相连，参数配置模块(1)的第四个输出端与坐标转换模块(5)的第一个输入端相连，参数配置模块(1)的第五个输出端与定位解算模块(6)的第一个输入端相连，参数配置模块(1)的第六个输出端与卡尔曼滤波器(7)的第一个输入端相连，参数配置模块(1)的第七个输出端与报文形成模块(8)的第一个输入端相连；导航电文从数据预处理器(2)的第二个输入端输入，观测数据帧从数据预处理器(2)的第三个输入端输入，数据预处理器(2)的第一个输出端与卫星参数解算模块(3)的第二个输入端相连，数据预处理器(2)的第二个输出端与伪距修正模块(4)的第二个输入端相连；卫星参数解算模块(3)的输出端与坐标转换模块(5)的第二个输入端相连，坐标转换模块(5)的输出端与定位解算模块(6)的第二个输入端相连；伪距修正模块(4)的输出端与定位解算模块(6)的第三个输入端相连；定位解算模块(6)的输出端与卡尔曼滤波器(7)的第二个输入端相连；卡尔曼滤波器(7)的第一个输出端与伪距修正模块(4)的第三个输入端相连，卡尔曼滤波器(7)的第二个输出端与报文形成模块(8)的第二个输入端相连；

其中，参数配置模块(1)读取兼容导航接收机的配置参数，并将配置参数输出给定位系统的各模块以控制其工作状态；

数据预处理器(2)采集有效卫星的导航电文和观测数据帧，并从中提取星历数据、历书数据和观测信息；

卫星参数解算模块(3)获取有效卫星的卫星位置和卫星速度；

伪距修正模块(4)获取有效卫星的观测伪距和伪距误差修正信息，并对观测伪距进行修正；

坐标转换模块(5)将来自不同卫星导航系统的有效卫星的坐标统一到同一参考坐标系中；

定位解算模块(6)获取接收机位置和接收机速度信息；

卡尔曼滤波器(7)对接收机的位置和速度信息进行滤波；

报文形成模块(8)形成导航报文。

2.一种利用权利要求1所述的兼容导航接收机定位系统的定位方法，其特征在于，该方法包括如下步骤：

步骤一：读取接收机的配置参数；

参数配置模块(1)读取接收机的配置参数，

首先，根据用户需求，确定接收机的工作状态和使用的卫星导航系统；工作状态包括：单卫星导航系统定位、双卫星导航系统组合定位、三卫星导航系统组合定位和四卫星导航系统组合定位，使用的卫星导航系统包括：全球定位系统GPS、全球导航卫星系统GLONASS、北斗二代导航系统BD2和“伽利略”卫星定位系统Galileo；

然后，根据用户需求，确定卡尔曼滤波器(7)的工作状态，确定导航信息的数据更新频率，确定导航报文的格式；

最后，将接收机的配置参数输出给数据预处理器(2)、卫星参数解算模块(3)、伪距修正模块(4)、坐标转换模块(5)、定位解算模块(6)、卡尔曼滤波器(7)和报文形成模块(8)，以控制各模块的工作状态；

其中，接收机的配置参数包括：接收机的工作状态、使用的卫星导航系统、卡尔曼滤波器(7)的工作状态、导航信息的数据更新频率以及导航报文的格式；

步骤二：提取星历数据、历书数据和观测数据；

根据参数配置模块(1)输出的接收机的工作状态和使用的卫星导航系统，数据预处理器(2)依据有效卫星所属的卫星导航系统的说明书对采集的有效卫星的导航电文进行拼接、校验、纠错，从中提取有效卫星的星历数据和历书数据；

根据参数配置模块(1)输出的接收机的工作状态和使用的卫星导航系统，数据预处理器(2)依据捕获跟踪系统提供的观测数据帧格式从采集的有效卫星观测数据帧中提取毫秒数、码片偏移、相位偏移和多普勒频率；

步骤三：获取有效卫星的参数信息；

首先，根据参数配置模块(1)输出的接收机的工作状态和使用的卫星导航系统，卫星参数解算模块(3)依据有效卫星所属的卫星导航系统的说明书利用有效卫星的星历数据获取其位置坐标和速度坐标；

然后，根据参数配置模块(1)输出的接收机的工作状态和使用的卫星导航系统，伪距修正模块(4)依据有效卫星所属卫星导航系统的说明书利用有效卫星的星历数据和历书数据获取其星钟修正误差、电离层修正误差和对流层修正误差；

之后，伪距修正模块(4)通过将有效卫星的周内秒计数、比特数、毫秒数、码片偏移、相位偏移进行累加获取其观测伪距；

最后，伪距修正模块(4)利用获取的有效卫星的星钟修正误差、电离层修正误差和对流层修正误差对该卫星的观测伪距进行误差修正，获取修正后的观测伪距；

步骤四：获取接收机位置和接收机速度；

首先，根据参数配置模块(1)输出的接收机的工作状态和使用的卫星导航系统，在单卫星导航系统定位时坐标转换模块(5)不参与系统工作，在多卫星导航系统组合定位时坐标转换模块(5)依据各卫星导航系统参考坐标系的关系利用7参数Bursa模型将来自不同卫星导航系统的有效卫星的坐标转换到同一参考坐标系，获取各有效卫星在同一参考坐标系中的三维坐标；

然后，根据参数配置模块(1)输出的接收机的工作状态和使用的卫星导航系统，定位解算模块(6)利用有效卫星的参数信息获取接收机位置和接收机速度；

A.获取接收机位置：

当有效卫星的数目不小于N+3时，利用有效卫星的位置坐标和修正后的观测伪距获取接收机的位置坐标偏移和钟差模型：

$\left[\begin{array}{c}\mathrm{\Δ}{x}_u\\ \mathrm{\Δ}{y}_u\\ \mathrm{\Δ}{z}_u\\ \mathrm{c\Δ}{t}_1\\ .\\ .\\ .\\ \mathrm{c\Δ}{t}_i\\ .\\ .\\ .\\ \mathrm{c\Δ}{t}_N\end{array}\right]=-{\left({\hat{H}}^T\hat{H}\right)}^{-1}\left({\hat{H}}^T\mathrm{\Δ\ρ}\right)$

其中，上述位置坐标偏移和钟差模型等号左侧部分为一个(N+3)×1的矩阵模块；上述位置坐标偏移和钟差模型等号右侧部分的

为一个n×(N+3)的矩阵模块，

为一个n×1的矩阵模块；

其中，i为卫星S^j所属卫星导航系统的编号，i＝1，2，…，N，N为用户接收机进行定位时使用的卫星导航系统的数目；j为用户接收机进行定位时使用的有效卫星S^j的编号，j＝1，…，o，o+1，…，p，…，q+1，…，n，n为用户接收机进行定位时使用的有效卫星总数，S¹～S^o为第一卫星导航系统的卫星，S^o+1～S^p为第二卫星导航系统的卫星，S^q+1～Sⁿ为第N卫星导航系统的卫星；

(Δx_u，Δy_u，Δz_u)^T为用户接收机的位置坐标偏移；

Δt_i为用户接收机时钟与第i卫星导航系统系统时钟之间的钟差，i＝1，2，…，N，Δt₁为用户接收机时钟与第一卫星导航系统系统时钟之间的钟差，……，Δt_N为用户接收机时钟与第N卫星导航系统系统时钟之间的钟差；

各项分别表示从用户接收机近似位置指向卫星S^j的单位向量的方向余弦，j＝1，…，o，o+1，…，p，…，q+1，…，n，

分别表示从用户接收机近似位置指向第一卫星导航系统中的卫星S¹的单位向量的方向余弦，

分别表示从用户接收机近似位置指向第一卫星导航系统中的卫星S^o的单位向量的方向余弦，

分别表示从用户接收机近似位置指向第二卫星导航系统中的卫星S^o+1的单位向量的方向余弦，分别表示从用户接收机近似位置指向第二卫星导航系统中的卫星S^p的单位向量的方向余弦，

分别表示从用户接收机近似位置指向第N卫星导航系统中的卫星S^q+1的单位向量的方向余弦，

分别表示从用户接收机近似位置指向第N卫星导航系统中的卫星Sⁿ的单位向量的方向余弦；

其中，

为用户接收机近似位置的位置坐标；

(x^j，y^j，z^j)^T为卫星S^j的位置坐标，j＝1，…，o，o+1，…，p，…，q+1，…，n；

为对卫星S^j的观测伪距进行星钟误差、电离层误差和对流层误差修正后的伪距，j＝1，…，o，o+1，…，p，…，q+1，…，n，为对第一卫星导航系统中的卫星S¹的观测伪距进行星钟误差、电离层误差和对流层误差修正后的伪距，

为对第一卫星导航系统中的卫星S^o的观测伪距进行星钟误差、电离层误差和对流层误差修正后的伪距，为对第二卫星导航系统中的卫星S^o+1的观测伪距进行星钟误差、电离层误差和对流层误差修正后的伪距，

为对第二卫星导航系统中的卫星S^p的观测伪距进行星钟误差、电离层误差和对流层误差修正后的伪距，

为对第N卫星导航系统中的卫星S^q+1的观测伪距进行星钟误差、电离层误差和对流层误差修正后的伪距，

为对第N卫星导航系统中的卫星Sⁿ的观测伪距进行星钟误差、电离层误差和对流层误差修正后的伪距；

为用户接收机的近似位置与卫星S^j间的几何距离，j＝1，…，o，o+1，…，p，…，q+1，…，n，

为用户接收机的近似位置与第一卫星导航系统中的卫星S¹间的几何距离，

为用户接收机的近似位置与第一卫星导航系统中的卫星S^o间的几何距离，为用户接收机的近似位置与第二卫星导航系统中的卫星S^o+1间的几何距离，

为用户接收机的近似位置与第二卫星导航系统中的卫星S^p间的几何距离，

为用户接收机的近似位置与第N卫星导航系统中的卫星S^q+1间的几何距离，

为用户接收机的近似位置与第N卫星导航系统中的卫星Sⁿ间的几何距离；

c为光速；

用户接收机的位置坐标为：

$\left[\begin{array}{c}{x}_u\\ y_u\\ z_u\end{array}\right]=\left[\begin{array}{c}{\hat{x}}_u\\ {\hat{y}}_u\\ {\hat{z}}_u\end{array}\right]+\left[\begin{array}{c}\mathrm{\Δ}{x}_u\\ \mathrm{\Δ}{y}_u\\ \mathrm{\Δ}{z}_u\end{array}\right];$

B.获取接收机速度：

当有效卫星的数目不小于N+3时，利用接收机的位置坐标和有效卫星的多普勒频率获取接收机的速度坐标和钟漂模型：

$\left[\begin{array}{c}{v}_{\mathrm{ux}}\\ v_{\mathrm{uy}}\\ v_{\mathrm{uz}}\\ c{\stackrel{\·}{t}}_1\\ .\\ .\\ .\\ c{\stackrel{\·}{t}}_i\\ .\\ .\\ .\\ c{\stackrel{\·}{t}}_N\end{array}\right]=-{\left(H^{T}H\right)}^{-1}\left(H^{T}d\right)$

其中，上述速度坐标和钟漂模型等号左侧部分为一个(N+3)×1的矩阵模块；

上述速度坐标和钟漂模型等号右侧部分的

为一个n×(N+3)的矩阵模块，

d＝[d¹…d^o d^o+1…d^p…d^q+1…dⁿ]^T，为一个n×1的矩阵模块，

其中 $d^j=d_u^j-\left[\begin{array}{ccc}{a}_x^j& a_y^j& a_z^j\end{array}\right]\left[\begin{array}{c}{v}_x^j\\ v_y^j\\ v_z^j\end{array}\right],$

j为接收机定位时使用的有效卫星S^j的编号，j＝1，…，o，o+1，…，p，…，q+1，…，n，其中n为接收机定位使用的有效卫星的数目，S¹～S^o为第一卫星导航系统的卫星，S^o+1～S^p为第二卫星导航系统的卫星，S^q+1～Sⁿ为第N卫星导航系统的卫星；

i为卫星S^j所属卫星导航系统的编号，i＝1，2，…，N，其中N为用户接收机进行定位时使用的卫星导航系统的数目；

各项分别表示从用户接收机指向卫星S^j的单位向量的方向余弦，j＝1，…，o，o+1，…，p，…，q+1，…，n，

分别表示从用户接收机指向第一卫星导航系统中的卫星S¹的单位向量的方向余弦，

分别表示从用户接收机指向第一卫星导航系统中的卫星S^o的单位向量的方向余弦，

分别表示从用户接收机指向第二卫星导航系统中的卫星S^o+1的单位向量的方向余弦，

分别表示从用户接收机指向第二卫星导航系统中的卫星S^p的单位向量的方向余弦，分别表示从用户接收机指向第N卫星导航系统中的卫星S^q+1的单位向量的方向余弦，

分别表示从用户接收机指向第N卫星导航系统N中的卫星Sⁿ的单位向量的方向余弦；

(x_u，y_u，z_u)^T为用户接收机的位置坐标；

(x^j，y^j，z^j)^T为卫星S^j的位置坐标，j＝1，…，o，o+1，…，p，…，q+1，…，n；

为用户接收机与卫星S^j间的几何距离，j＝1，…，o，o+1，…，p，…，q+1，…，n；(v_ux，v_uy，v_uz)^T为用户接收机的速度坐标；

为用户接收机时钟与第i卫星导航系统系统时钟之间的钟差漂移，简称钟漂，i＝1，2，…，N，

为用户接收机时钟与第一卫星导航系统系统时钟之间的钟漂，……，为用户接收机时钟与第N卫星导航系统系统时钟之间的钟漂；

为卫星S^j的观测伪距率，j＝1，…，o，o+1，…，p，…，q+1，…，n，

为第一卫星导航系统中的卫星S¹的观测伪距率，

为第一卫星导航系统中的卫星S^o的观测伪距率，

为第二卫星导航系统中的卫星S^o+1的观测伪距率，

为第二卫星导航系统中的卫星S^p的观测伪距率，为第N卫星导航系统中的卫星S^q+1的观测伪距率，

为第N卫星导航系统中的卫星Sⁿ的观测伪距率；

为卫星S^j的速度坐标，j＝1，…，o，o+1，…，p，…，q+1，…，n；

c为光速；

最后，根据参数配置模块(1)输出的接收机的工作状态、使用的卫星导航系统和卡尔曼滤波器(7)的工作状态，确定卡尔曼滤波器(7)是否对定位解算模块(6)输出的接收机位置和接收机速度信息进行滤波处理；如果卡尔曼滤波器(7)的工作状态是工作，基于机动目标“当前”统计模型建立系统模型和测量模型，并通过自适应卡尔曼滤波获取接收机位置和接收机速度的滤波结果；

步骤五：将导航信息生成报文；

根据参数配置模块(1)输出的用户对导航报文的要求，报文形成模块(8)将导航信息以用户需要的格式形成报文。

3.根据权利要求2所述的一种兼容导航接收机定位方法，其特征在于，步骤四中所述的滤波结果采用递推方法给出：

状态估计模型

为第二系统量测模型；

估计均方误差模型P(k)＝[I-K(k)H(k)]P(k|k-1)；

其中，状态一步预测模型

一步预测均方误差模型

P(k|k-1)＝Φ(k|k-1)P(k-1)Φ^T(k|k-1)+Q(k-1)，Q(k-1)为系统噪声方差模型；

滤波增益模型K(k)＝P(k|k-1)H^T(k)[H(k)P(k|k-1)H^T(k)+R(k)]^-1；

其中，系统状态模型

x_u，

分别为接收机在参考坐标系x轴方向的位置、速度和加速度，y_u，分别为接收机在参考坐标系y轴方向的位置、速度和加速度，z_u，

分别为接收机在参考坐标系z轴方向的位置、速度和加速度，Δt_i，

分别为接收机时钟与第i卫星导航系统系统时钟之间的钟差和钟漂，i为接收机定位使用的卫星导航系统的编号，i＝1，2，…，N，N为接收机定位时使用卫星导航系统的数目，Δt₁，

分别为接收机时钟与第一卫星导航系统系统时钟之间的钟差和钟漂，……，Δt_N，

分别为接收机时钟与第N卫星导航系统系统时钟之间的钟差和钟漂，初值为

X₀＝[x_u，v_ux，0，y_u，v_uy，0，z_u，v_uz，0，0，0，…，0，0]^T，(x_u，y_u，z_u)^T和(v_ux，v_uy，v_uz)^T分别为定位解算模块(6)输出的接收机位置坐标和速度坐标；

其中，第一系统状态转移模型

${\mathrm{\φ}}_{1x}\left(k\right|k-1)={\mathrm{\φ}}_{1y}\left(k\right|k-1)={\mathrm{\φ}}_{1z}\left(k\right|k-1)=\left[\begin{array}{ccc}1& T& T^2/2\\ 0& 1& T\\ 0& 0& 1\end{array}\right],$

${\mathrm{\φ}}_{{1t}_1}\left(k\right|k-1)=...={\mathrm{\φ}}_{{1t}_N}\left(k\right|k-1)=\left[\begin{array}{cc}1& T\\ 0& 1\end{array}\right];$

T为采样周期；

P(k-1)为状态协方差，初值P₀取为对角阵；

其中，第二系统状态转移模型

$\mathrm{\Φ}\left(k\right|k-1)=\mathrm{diag}[{\mathrm{\φ}}_x\left(k\right|k-1),{\mathrm{\φ}}_y\left(k\right|k-1),{\mathrm{\φ}}_z\left(k\right|k-1),{\mathrm{\φ}}_{t_1}\left(k\right|k-1),...,{\mathrm{\φ}}_{t_N}\left(k\right|k-1)],$

${\mathrm{\φ}}_x\left(k\right|k-1)=\left[\begin{array}{ccc}1& T& {\mathrm{\α}}_x^{-2}(-1+{\mathrm{\α}}_{x}T+e^{-{\mathrm{\α}}_{x}T})\\ 0& 1& (1-e^{-{\mathrm{\α}}_{x}T}){\mathrm{\α}}_x^{-1}\\ 0& 0& e^{-{\mathrm{\α}}_{x}T}\end{array}\right],$

${\mathrm{\φ}}_y\left(k\right|k-1)=\left[\begin{array}{ccc}1& T& {\mathrm{\α}}_y^{-2}(-1+{\mathrm{\α}}_{y}T+e^{-{\mathrm{\α}}_{y}T})\\ 0& 1& (1-e^{-{\mathrm{\α}}_{y}T}){\mathrm{\α}}_y^{-1}\\ 0& 0& e^{-{\mathrm{\α}}_{y}T}\end{array}\right],$

${\mathrm{\φ}}_z\left(k\right|k-1)=\left[\begin{array}{ccc}1& T& {\mathrm{\α}}_z^{-2}(-1+{\mathrm{\α}}_{z}T+e^{-{\mathrm{\α}}_{z}T})\\ 0& 1& (1-e^{-{\mathrm{\α}}_{z}T}){\mathrm{\α}}_z^{-1}\\ 0& 0& e^{-{\mathrm{\α}}_{z}T}\end{array}\right],$

${\mathrm{\φ}}_{t_1}\left(k\right|k-1)=...={\mathrm{\φ}}_{t_N}\left(k\right|k-1)=\left[\begin{array}{cc}1& T\\ 0& 1\end{array}\right],$

α_x，α_y，α_z为机动时间常数的倒数；

其中，系统噪声方差模型

$Q\left(k\right)=\mathrm{diag}[2{\mathrm{\σ}}_{a_x}^2{\mathrm{\α}}_x{Q}_x\left(k\right),2{\mathrm{\σ}}_{a_y}^2{\mathrm{\α}}_y{Q}_y\left(k\right),2{\mathrm{\σ}}_{a_z}^2{\mathrm{\α}}_z{Q}_z\left(k\right),Q_{t_1}\left(k\right),...,Q_{t_N}\left(k\right)],$

$Q_x\left(k\right)=\left[\begin{array}{ccc}{q}_{x11}& q_{x12}& q_{x13}\\ q_{x12}& q_{x22}& q_{x23}\\ q_{x13}& q_{x23}& q_{x33}\end{array}\right],$ $Q_{t_1}\left(k\right)=...=Q_{t_N}\left(k\right)=\left[\begin{array}{cc}{\mathrm{\σ}}_{\mathrm{\Δt}}^2& 0\\ 0& {\mathrm{\σ}}_{\stackrel{\·}{t}}^2\end{array}\right],$

$q_{x11}=0.5{\mathrm{\α}}_x^{-5}(1-e^{-2{\mathrm{\α}}_{x}T}+2{\mathrm{\α}}_{x}T-4{\mathrm{\α}}_{x}T{e}^{-{\mathrm{\α}}_{x}T}-2{\mathrm{\α}}_x^2{T}^2+\frac{2}{3}{\mathrm{\α}}_x^3{T}^3),$

$q_{x12}=0.5{\mathrm{\α}}_x^{-4}(1-2e^{-{\mathrm{\α}}_{x}T}+e^{-{2\mathrm{\α}}_{x}T}-2{\mathrm{\α}}_{x}T+2{\mathrm{\α}}_{x}T{e}^{-{\mathrm{\α}}_{x}T}+{\mathrm{\α}}_x^2{T}^2),$

$q_{x13}=0.5{\mathrm{\α}}_x^{-3}(1-e^{-2{\mathrm{\α}}_{x}T}-2{\mathrm{\α}}_{x}T{e}^{-{\mathrm{\α}}_{x}T}),$

$q_{x23}=0.5{\mathrm{\α}}_x^{-2}(1-2e^{-{\mathrm{\α}}_{x}T}+e^{-{2\mathrm{\α}}_{x}T}),$

$q_{x22}=0.5{\mathrm{\α}}_x^{-3}(-3+4e^{-{\mathrm{\α}}_{x}T}-e^{-2{\mathrm{\α}}_{x}T}+2{\mathrm{\α}}_{x}T),$

$q_{x33}=0.5{\mathrm{\α}}_x^{-1}(1-e^{-2{\mathrm{\α}}_{x}T}),$

Q_y(k)，Q_z(k)中的元素模型与Q_x(k)中的元素模型相同，将Q_x(k)中各元素模型中的α_x分别用α_y，α_z代替，即得到Q_y(k)，Q_z(k)；

分别为与接收机在参考坐标系x，y，z坐标方向上的加速度分量对应的系统噪声方差，α_x，α_y，α_z为机动时间常数的倒数，

分别为与接收机钟差和钟漂对应的系统噪声方差；

其中，第一系统量测模型H(k)＝[H¹…H^j…Hⁿ]^T，

$H^j=\left[\begin{array}{cccccccc}{A}_x^j& A_y^j& A_z^j& B_1& ...& B_i& ...& B_N\end{array}\right],$

$A_x^j=\left[\begin{array}{ccc}-a_x^j& 0& 0\\ 0& a_x^j& 0\end{array}\right],$ $A_y^j=\left[\begin{array}{ccc}-a_y^j& 0& 0\\ 0& a_y^j& 0\end{array}\right],$ $A_z^j=\left[\begin{array}{ccc}-a_z^j& 0& 0\\ 0& a_z^j& 0\end{array}\right],$

当卫星S^j属于第i卫星导航系统时，

当卫星S^j不属于第i卫星导航系统时，

j为用户接收机进行定位时使用的有效卫星S^j的编号，j＝1，2，…，n，n为用户接收机进行定位时使用的有效卫星的数目，

i为接收机定位时使用的卫星导航系统的编号，i＝1，2，…，N，N为接收机定位使用的卫星导航系统的数目，

各项分别表示从用户接收机指向卫星S^j的单位向量的方向余弦，(x_u，y_u，z_u)^T为用户接收机的位置坐标，(x^j，y^j，z^j)^T为卫星S^j的位置坐标，

为用户接收机与卫星S^j间的几何距离，j＝1，2，…，n；

其中，量测噪声模型

分别为卫星S^j的伪距观测噪声方差和伪距率观测噪声方差，j＝1，2，…，n，

分别为卫星S¹的伪距观测噪声方差和伪距率观测噪声方差，

分别为卫星Sⁿ的伪距观测噪声方差和伪距率观测噪声方差；

其中，观测模型

分别为对卫星S^j的观测伪距进行星钟误差、电离层误差和对流层误差修正后的伪距和观测伪距率，j＝1，2，…，n，

分别为对卫星S¹的观测伪距进行星钟误差、电离层误差和对流层误差修正后的伪距和观测伪距率，

分别为对卫星Sⁿ的观测伪距进行星钟误差、电离层误差和对流层误差修正后的伪距和观测伪距率；

其中，第二系统量测模型h(X)＝[h¹…h^j…hⁿ]^T

$h^j=\left[\begin{array}{c}\sqrt{{(x^j-x_u)}^2+{(y^j-y_u)}^2+{(z^j-z_u)}^2}+\underset{i=1}{\overset{N}{\mathrm{\Σ}}}c{k}_i^j\mathrm{\Δ}{t}_i\\ a_x^j(v_{\mathrm{ux}}-v_x^j)+a_y^j(v_{\mathrm{uy}}-v_y^j)+a_z^j(v_{\mathrm{uz}}-v_z^j)-\underset{i=1}{\overset{N}{\mathrm{\Σ}}}c{k}_i^j{\stackrel{\·}{t}}_i\end{array}\right]$

当卫星S^j属于第i卫星导航系统时，

当卫星S^j不属于第i卫星导航系统时，

j为用户接收机进行定位时使用的有效卫星S^j的编号，j＝1，2，…，n，其中n为接收机定位使用的有效卫星的数目，

i为接收机定位时使用的卫星导航系统的编号，i＝1，2，…，N，N为接收机定位使用的卫星导航系统的数目。

Images