一种基于GNSS的单频高精度定位方法与流程

本发明提供一种基于GNSS(Global Navigation Satellite System)的高精度单频定位方法，它涉及一种对飞行器或地面，海面运动载体利用GNSS进行绝对位置解算的方法，属于导航技术领域。

背景技术：

全球定位系统(GPS，Global Positioning System)是美国建立的一种卫星导航系统,能够为用户提供实时的三维导航服务。伴随着GPS的广泛应用，许多科研机构开始注意其单个接收机在低轨卫星精密定轨，以及地面和海面，静止或运动载体的高精度定位方面的应用。目前，上述方面的研究可以分为两个方向。第一，采用双频GPS的精密定位；该方法充分利用GPS播发的多频信息，通过各种不同的线性组合构建出无电离层观测量，从而消除电离层延迟的影响，以提高导航精度。采用双频GPS进行的运动学后处理低轨卫星定轨精度在厘米级。但双频方法对接收机要求较高，数据处理复杂，计算量大。第二，采用单频GPS的精密定位；单频GPS定位只利用GPS卫星播发的一个频率上的信息，其只包含伪距信息和载波信息，采用电离层模型或GRAPHIC(Group and Phase Ionospheric Calibration)消电离层。相较于双频GPS精密定位，单频方法的精度要低一些，其定位精度一般在分米级，但由于单频接收机造价低廉，数据量小，处理简单，使得单个单频接收机的高精度定位方法研究仍具有巨大的实用意义。

如今，各国的卫星导航系统都在如火如荼的建设当中，除了已经成熟完善的美国的GPS系统，我国的北斗卫星导航系统(BeiDou Navigation Satellite System)，俄罗斯的格洛纳斯(GLONASS，GLOBAL NAVIGATION SATELLITE SYSTEM)以及欧洲的伽利略卫星导航系统(Galileo satellite navigation system)都在建设中。可以预想未来众多卫星导航系统的建立完善会使得基于GNSS的精密导航定位服务研究拥有更广阔的平台，而本专利所提出的一种基于GNSS的单频高精度定位方法便是其中之一。

在进行精密定位的过程中，电离层延迟误差是最大的误差来源之一。在双频GPS中利用电离层延迟一阶项与载波频率的关系，通过两个载波上信息的线性组合来消除电离层延迟的影响。在单频GPS中利用电离层延迟一阶项在伪距观测量与载波观测量中大小相等符号相反的性质，构建GRAPHIC(Group and Phase Ionospheric Calibration)观测量消除电离层延迟的影响。而由于伪距观测量的噪声为半米左右，载波观测量的噪声为毫米级，导致了构建的GRAPHIC组合观测量的噪声约为伪距噪声的一半，可达分米级。大噪声的观测量严重限制了定位结果的精度。

在不同的情况与目的下，对所得到的伪距观测量与载波相位观测量有许多处理方法。目前在导航定位方面主要有星间单差消除接收机卫星钟差、站间单差消除GPS卫星钟差以及历元单差消除载波相位模糊度，三种单差又可以进行组合从而得到更加复杂的双差观测量以及三差观测量。

综上，基于GNSS的精密定位具有巨大的研究与应用价值，其中的采用单接收机单频GPS测量量的定位算法，其设备简单，间隔低廉，数据处理容易，具有更大的普及与实用价值。但由于其定位精度相比双频GPS较差，从而限制了其在高精度情景下的应用。本发明提出了一种基于GNSS的单频高精度定位方法，提高了单频GPS定位精度，从而实现了低轨卫星的高精度定轨、地面、水面以及空中的静止或运动载体高精度定位与导航。

技术实现要素：

(一)发明目的：本发明采用单接收机的单频GPS接收机数据，以运动学GRAPHIC组合定位方法为基础，创新的在每个历元定位过程中引入该历元时刻位置信息虚拟观测量，而当前历元位置信息虚拟观测量是使用上一历元的绝对位置加上载波相位历元差分所求解的相对位置变化量得到。由于载波相位观测量的观测精度很高，其历元差分求解的相对位置变化量误差在厘米级，在传统的GRAPHIC组合运动学定轨中加入这些信息，增加了观测方程，而没有引入大的误差，可以有效的抑制传统GRAPHIC组合定位结果的噪声，因此对于定位精度有很大的提高。

(二)技术方案

本发明一种基于GNSS的单频高精度定位方法，其步骤如下：

步骤一：准备工作

首先给出伪距观测方程：

P＝ρ+c(δtr-δt^s)+δion+δrel+εp (1)

式(1)中，P伪距观测量，以米(m)为单位；ρ为GPS卫星到接收机的真实几何距离；c为真空中的光速；δtr为接收机钟差；δt^s为GPS卫星钟差；δion为电离层延迟；δrel为相对论效应引起的误差；εp为伪距观测噪声，其符合方差为的白噪声；

给出载波相位观测方程：

式(2)中，Φ为载波相位观测量(以米为单位)，为载波相位观测量(以周为单位)；λ为载波波长；ρ为GPS卫星到接收机的真实几何距离；c为真空中光速；δtr为接收机钟差；δt^s为GPS卫星钟差；N为模糊度；δion为电离层延迟；δrel为相对论效应引起的误差；为载波相位观测噪声，其符合方差为的白噪声；

为了消除电离层误差，利用伪距观测方程(1)和载波相位观测方程(2)中电离层延迟一阶项，大小相等，符号相反的性质，得到以下消除电离层误差项的GRAPHIC组合观测方程为：

上式中δt^s可以通过外部星历得到，为已知项；δrel可以通过建模得到并进行修正。合并已知项的观测方程为：

Gc＝ρ+cδtr+λNG+εG (4)

其中Gc为误差修正后的GRAPHIC组合观测量，NG为GRAPHIC组合模糊度，表达

式如下：

下文在不引起歧义的情况下，统一将NG简写为N，并称为模糊度；

εG为GRAPHIC组合观测噪声，表达式如下：

其噪声方差如下：

式(7)是GRAPHIC组合观测量的噪声方差，其中D(·)为求方差算子(下文出现，不再说明)，可见其是由伪距噪声方差与载波噪声方差组合得到；伪距的测量噪声可达米级，而载波相位的测量噪声一般为毫米级；由上式可以看到GRAPHIC组合的噪声水平处于伪距噪声和载波相位噪声之间，为分米级，使用这种组合进行传统的GRAPHIC组合运动学定轨、定位和导航，结果的误差在分米级，难以满足高精度定位要求；

给出载波相位历元差分观测方程：

上式(8)中，k表示当前历元，k-1表示上一历元；其他符号与上文相同。其中，GPS卫星钟差项δt^s可通过外部星历得到；相对论引起的误差项δrel可通过建模得到；而当历元时间较短时，相邻历元电离层变化不大，残留项(δion,k-δion,k-1)可以忽略；因此得到整理后的载波相位历元差分观测方程：

上式中ΔΦc为GPS接收机钟差，相对论效应修正后的载波相位历元差分观测量；其观测噪声方差为：

从(10)中可以看出，载波相位历元差分观测量的噪声水平与载波相位观测量的水平相当；

以上，给出了载波相位历元差分测量模型，该模型用于相邻历元的相对位置求解，其结果将传递到下一历元用于构建位置信息虚拟观测量；使用最小二乘估计对当前时刻接收机的位置和钟差信息进行估计；以下给出一种基于GNSS的单频高精度定位方法具体步骤；

步骤二：初始历元位置解算

由于初始历元是导航定位过程的第一个历元，其没有来自上一历元的先验信息，因此对于初始历元要进行单独处理，初始化它的位置信息与模糊度信息；

构建GRAPHIC组合观测量计算值Z；

由公式(4)，可知GRAPHIC组合观测方程是一个关于状态变量x＝(xr,yr,zr,δtr,N1…Np)^T的非线性方程，其中模糊度下标表示当前历元共观测到的GPS卫星数目。GRAPHIC组合观测量计算值Z，计算如下：

其中的关系式如下：

其中，上式符号“～”表示估计值，为接收机的三个位置估计分量，(x^s,y^s,z^s)^T为GPS卫星的三个位置分量；

采用卡尔曼滤波器一步估计状态变量x＝(xr,yr,zr,δtr,N1…Np)^T。因为GRAPHIC组合观测方程是关于状态变量x＝(xr,yr,zr,δtr,N1…Np)^T的非线性方程，因此要对其线性化；线性化过程如下：

将GRAPHIC组合观测方程(4)，写成如下形式：

Gc＝h(x)+εG (13)

在参考点处线性化：

其中便是上面所求的GRAPHIC组合观测量的计算值Z。Δx为状态变量x的实际值与估计值之差；H模型观测量对状态量x的偏导数雅克比矩阵，有：

认为每个GRAPHIC组合观测量是相互独立的，观测噪声阵有如下形式：

状态变量初始值的先验协方差阵P^-如下：

式(17)，是一个分块对角阵。其中为位置信息初始先验方差阵；为接收机钟差初始先验方差阵；为模糊度方差阵，取值如下：

其中，取值一般为伪距单点定位误差的一半，现取30米；

式(19)，表示初始历元钟差δtr的初始带入值是准确的，在后面估计出的各个历元接收机钟差项全都是相对于第一历元时刻的相对钟差，后边不再叙述；

式(20)，模糊度方差值为一个极大值，这是因为初始历元并没有模糊度信息的先验值，所以如此处理表示方差的初始带入值是完全不可信的；

给状态变量的估计值赋初值,其中位置状态变量的估计值赋值如下：

其中，(xspp,yspp,zspp)^T表示由伪距单点定位得出的接收机位置，单点定位在此不再叙述。

钟差状态变量的赋值如下：

上式与式(19)联立，其意义表示初始历元钟差为0，并且是准确的。以后历元的钟差估计值皆表示相对于初始历元的相对钟差；

模糊度状态变量估计值的赋值如下：

卡尔曼滤波，滤波系数如下：

K＝P^-H^T(HP^-H^T+R) (24)

则状态量的估计值更新如下：

上式中，上标“+”表示更新后的量，下文出现不再叙述；

以式(25)的一步估计值，作为初始历元的最终解算结果；

状态变量方差阵P更新如下：

P⁺＝(1-KH)P^- (26)

步骤三：构建载波相位历元差分观测模型，求解相对位置变化量

由式(9)给出的载波相位历元差分方程，带入k历元与k-1历元的坐标位置如下：

上式中，GPS卫星位置、钟差以及k-1历元接收机位置皆为已知量；

k历元与k-1历元接收机位置、钟差有如下关系：

上式中，Δx,Δy,Δz,Δδtr分别表示从k-1历元到k历元，接收机位置变换量与钟差变化量；

将式(28)带入式(27)，则ΔΦc成为一个关于状态变量y＝(Δx,Δy,Δz,cΔδtr)^T的非线性函数，记为下式：

ΔΦc＝f(Δx,Δy,Δz,Δδtr) (29)

给出状态变量估计值的初始值如下:

计算载波相位历元差分观测量的计算值ΔL:

非线性方程(29)在带入估计值泰勒展开，取一阶项，求H阵：

上式，各个量含义与(14)中各量含义类似，不再叙述；

由式(10)，给出载波相位历元差分观测量的噪声阵R，如下：

由最小二乘估计：

Δy＝(H^TR^-1H)^-1(H^TR^-1(ΔΦc-ΔL)) (34)

更新状态变量的估计值：

将作为新的初值，带入式(30)，计算式(31)～(35)；重复以上过程直到满足下式：

||Δy||＜10^-3m (36)

或超出最大迭代次数时终止；

状态变量y的方差阵P由下式给出：

上式中PdR为位置变化量信息方差阵；PtdR,PdRt为位置变化量与钟差变化量协方差阵；Pt为钟差变化量方差阵；I为单位阵；Q为电离层控制因子，用来控制不同历元间电离层残差影响；

步骤四：模糊度信息与位置信息传递

k-1历元的位置信息与模糊度信息已知，由步骤三也已解算出k-1历元到k历元的位置变化量；此步骤便是要把位置信息与模糊度信息传递到k历元；

k历元的位置信息传递如下：

上式中下标“0”，表示初值；

K历元的模糊度信息传递如下：

方差信息传递如下：

上式中下标“0”，表示初值；

因为不对钟差进行传递，所以，对式(41)中钟差的方差设为10⁶，对均设为0；

步骤五：模糊度信息与方差信息调整

若从k-1历元到k历元若发生换星，则进行该步骤。否则跳过步骤五，执行步骤六；

从k-1历元到k历元若发生换星，位置信息不发生变化，模糊度改变；要对步骤四中的式(40)与式(41)中和模糊度有关的项进行修改；

对式(40)的调整，保留k-1历元与k历元共同观测的GPS卫星模糊度，去掉只在k-1历元观测到的卫星模糊度，增加k历元新增的GPS卫星模糊度，将其模糊度设为0。得到模糊度信息为(N0)q×1，‘q’表示k历元所观测到的总的GPS卫星个数；

对式(41)的调整，保留k-1历元与k历元共同观测的GPS卫星模糊度方差与协方差，去掉只在k-1历元观测到的卫星的方差与协方差，增加k历元新增的GPS卫星方差与协方差，方差值取10⁶，协方差取0；得到的模糊度方差阵为

步骤六：构建模糊度与位置信息虚拟观测方程

使用步骤四传递，步骤五调整过的位置信息与模糊度信息作为虚拟测量的测量值，如下：

上式中x,y,z为接收机三个位置坐标，标“v”，表示虚拟测量；q表示k历元所观测到的GPS卫星总数；

虚拟观测方程如下：

上式中，(Nk)q×1为k历元代估模糊度向量；

观测方程矩阵形式如下：

上式中Hv表达式如下：

Hv＝I(3+q)×(3+q) (45)

虚拟观测量的噪声阵由下式给出：

若发生了换星，则式(46)中PN0由步骤五得到；若未发生换星，则

步骤七：GRAPHIC组合观测方程与虚拟观测方程共同进行最小二乘估计，解算最终定位结果

观测方程如下：

上式中，HG的计算步骤与步骤二中GRAPHIC组合H阵计算相同；Hv由步骤六中给出；υ为噪声向量；

观测量的噪声阵，由下式给出：

上式中RG由(16)给出；

最小二乘估计：

Δx＝(H^TR^-1H)^-1R^-1(z-Z) (49)

上式中，z为观测量的测量值，如下：

Z为观测量的计算值，如下：

上式中G'由式(11)计算，V'为下式：

更新状态变量x：

将更新后的状态量估计值作为新的初值带入式(47)计算H阵；重复(47)～(53)，直到满足下式：

或迭代次数超过最大迭代次数时停止；

得到状态量的最小二乘估计，其k历元的最终定位结果为(xk yk zk)^T，计算状态变量的协方差阵：

回到步骤三，重复以后步骤，进行下一个历元的相对定位，信息传递以及运动学定位解算；

通过上述步骤，提出了一种基于GNSS的单频高精度定位方法；该方法仅采用单个接收机的单频GPS数据，使用载波相位历元差分辅助的GRAPHIC组合运动学定位算法，提高了运动学导航定位精度；采用精密星历可以对低轨卫星进行高精度定轨，采用广播星历则可以实现实时导航定位。该方法，不受导航定位过程中换星影响，可靠度高。并且该方法的实现只需一个单频GPS接收机，具有成本低廉的优势，数据处理简单，利于普及与应用。

(三)优点

本发明提供的一种基于GNSS的单频高精度定位方法，即一种基于GNSS的单频高精度定位方法相较于传统的GNSS导航方法的优点在于：

①本发明中，创新的使用了载波相位历元差分求解相对位置变化量辅助GRAPHIC组合运动学导航，极大的抑制了单频接收机GRAPHIC运动学定位的结果噪声。

②本发明提出的方法，具有更加广泛的应用性。使用广播星历，可为船舶、汽车、航拍无人机等需实时定位的载体提供高精度定位结果。使用精密星历进行事后处理，可以得到更高精度的定位结果；

③本发明提出的方法，实现成本低廉。只需要一个单频GNSS接收机便可实现，对接收机的性能要求大幅度降低，成本控制方面具有突出优势，在民用高精度定位方面具有普及优势。

附图说明

图1是本发明所述方法流程图。

具体实施方式

下面将结合图1和技术方案对本发明的具体实施过程做进一步的详细说明。

本发明一种基于GNSS的单频高精度定位方法，见图1所示，其步骤如下：

步骤一：初始历元位置解算

构建GRAPHIC组合观测量理论值Z。

由公式(4)，可知GRAPHIC组合观测方程是一个关于状态变量x＝(xr,yr,zr,δtr,N1…Np)^T的非线性方程，其中模糊度下标表示当前历元共观测到的GPS卫星数目。GRAPHIC组合观测量计算值Z，计算如下：

其中ρ0的关系式如下：

其中，上式符号“～”表示估计值，为接收机的三个位置估计分量，(x^s,y^s,z^s)^T为GPS卫星的三个位置分量。

采用卡尔曼滤波器一步估计状态变量x＝(xr,yr,zr,δtr,N1…Np)^T。因为GRAPHIC组合观测方程是关于状态变量x＝(xr,yr,zr,δtr,N1…Np)^T的非线性方程，因此要对其线性化。线性化过程如下：

将GRAPHIC组合观测方程(4)，写成如下形式：

Gc＝h(x)+εG (58)

在估计值处线性化：

其中便是上面所求的GRAPHIC组合观测量的计算值Z。Δx为状态变量x的实际值与估计值之差。H模型观测量对状态量x的偏导数雅克比矩阵，有：

认为每个GRAPHIC组合观测量是相互独立的，观测噪声阵有如下形式：

状态变量的协方差阵P^-如下：

式(62)，是一个分块对角阵。其中为位置信息先验方差阵；为接收机钟差先验方差阵；为模糊度先验方差阵，取值如下：

其中，取值一般为伪距单点定位误差的一半，现取30米。

式(64)，表示初始历元钟差δtr的初始估计值是准确的，在后面估计出的各个历元接收机钟差项全都是相对于第一历元时刻的相对钟差，下文不再叙述。

式(65)，给模糊度了一个非常大的方差值，这是因为初始历元并没有模糊度信息的先验值，所以如此处理表示方差的初始带入值是完全不可信的。给状态变量的估计值赋初值,其中位置状态变量的估计值赋值如下：

其中，(xspp,yspp,zspp)^T表示由伪距单点定位得出的接收机位置，单点定位在此不再叙述。

钟差状态变量的赋值如下：

上式与式(64)联立，其意义表示初始历元钟差为0，并且是准确的。以后历元的钟差估计值皆表示相对于初始历元的相对钟差。

模糊度状态变量估计值的赋值如下：

卡尔曼滤波，滤波系数如下：

K＝P H^T(HP^-H^T+R) (69)

则状态量的估计值更新如下：

上式中，上标“+”表示更新后的量，下文出现不再叙述。

以式(70)的一步估计值，作为初始历元的最终解算结果。

状态变量方差阵P更新如下：

P⁺＝(1-KH)P^- (71)

该步骤对应图1中的第二个方框。

步骤三：构建载波相位历元差分观测模型，求解位置变化量

由式(9)给出的载波相位历元差分方程，带入k历元与k-1历元的坐标位置如下：

上式中，GPS卫星位置、钟差以及k-1历元接收机位置皆为已知量。

k历元与k-1历元接收机位置、钟差有如下关系：

上式中，Δx,Δy,Δz,Δδtr分别表示从k-1历元到k历元，接收机位置变换量与钟差变化量。

将式(73)带入式(72)，则ΔΦc成为一个关于状态变量y＝(Δx,Δy,Δz,cΔδtr)^T的非线性函数，记为下式：

ΔΦc＝f(Δx,Δy,Δz,Δδtr) (74)

给出状态变量估计值的初始值如下:

计算载波相位历元差分观测量的计算值ΔL:

非线性方程(74)在带入估计值处线性化，求H阵：

上式，各个量含义与(14)中各量含义类似，不再叙述。

由式(10)，给出载波相位历元差分观测量的噪声阵R，如下：

由最小二乘估计：

更新状态变量的估计值：

将作为新的初值，带入式(75)，计算式(76)～(80)。重复以上过程直到满足下式：

||Δy||＜10^-3m (81)

或超出最大迭代次数时终止。

状态变量y的方差阵P由下式给出：

上式中PdR为位置变化量信息协方差阵；PtdR,PdRt为位置变化量与钟差变化量协方差阵；Pt为钟差变化量协方差阵；I为单位阵；Q为电离层控制因子。

该步骤对应图1中的第三个方框。

步骤四：模糊度信息与位置信息传递

k-1历元的位置信息与模糊度信息已知，由步骤三也已解算出k-1历元到k历元的位置变化量。此步骤便是要把位置信息与模糊度信息传递到k历元。

k历元的位置信息传递如下：

上式中下标“0”，表示初值。

K历元的模糊度信息传递如下：

方差信息传递如下：

上式中下标“0”，表示初值。

不对钟差进行传递，所以，对式(86)中钟差的方差设为10⁶，对均设为0。

该步骤对应图1中的第四个方框。

步骤五：构建GRAPHIC组合观测模型

如步骤一中，构建k历元的GRAPHIC组合观测模型，得到GRAPHIC组合测量量的计算值Z。

该步骤对应图1中的第五个方框。

步骤六：判断是否发生换星

判断从k-1历元到k历元，是否发生换星。若是，则执行步骤七；否则，执行步骤八。

该步骤对应附图一种的第六个方框

步骤七：模糊度信息与方差信息的调整

从k-1历元到k历元若发生换星，位置信息不发生变化，模糊度改变。要对步骤四中的式(85)与式(86)中和模糊度有关的项进行修改。

对式(85)的调整，保留k-1历元与k历元共同观测的GPS卫星模糊度，去掉只在k-1历元观测到的卫星模糊度，增加k历元新增的GPS卫星模糊度，将其模糊度设为0。得到模糊度信息为(N0)q×1。

对式(86)的调整，保留k-1历元与k历元共同观测的GPS卫星模糊度方差与协方差，去掉只在k-1历元观测到的卫星的方差与协方差，增加k历元新增的GPS卫星方差与协方差，方差值取10⁶，协方差取0。得到的模糊度方差阵为

该步骤对应图1中的第七个方框。

步骤八：构建虚拟测量方程

使用步骤四传递，步骤七调整过的位置信息与模糊度信息作为虚拟测量量的测量值，如下：

上式中下标“v”，表示虚拟测量；q表示k历元所观测到的GPS卫星总数。

虚拟观测方程如下：

上式中，(Nk)q×1为k历元代估模糊度向量。

观测方程矩阵形式如下：

上式中Hv表达式如下：

Hv＝I(3+q)×(3+q) (90)

虚拟测量量的噪声阵由下式给出：

若发生了换星，则上式中PN0由步骤七得到；若未发生换星，则

该步骤对应图1中的第八个方框。

步骤九：GRAPHIC组合观测方程与虚拟观测方程共同进行最小二乘估计，解算最终定位结果

测量方程如下：

上式中，HG的计算步骤与步骤一中GRAPHIC组合H阵计算相同；Hv由步骤八中给出；υ为噪声向量。

观测量的噪声阵，由下式给出：

上式中RG由(16)给出。

最小二乘估计：

Δx＝(H^TR^-1H)^-1R^-1(z-Z) (94)

上式中，z为测量量的测量值，如下：

Z为测量量的计算值，如下：

上式中G'由式(11)计算，V'为下式：

更新状态变量x：

将更新后的状态量估计值作为新的初值带入式(92)计算H阵。重复(92)～(98)，直到满足下式：

或迭代次数超过最大迭代次数时停止。

得到状态量的最小二乘估计，其k历元的最终定位结果为(xk yk zk)^T，计算状态变量的协方差阵：

该步骤对应附图一中的第九个方框。

步骤十：判断导航是否结束

若导航结束，则终止导航流程；否则执行步骤三。

该步骤对应图1中的第十个方框。

通过上述步骤，提出了一种基于GNSS的单频高精度定位方法。该方法仅采用单个单频GPS接收机，使用载波相位历元差分辅助的GRAPHIC组合运动学导航算法，抑制了定位噪声误差，提高了运动学导航定位精度。采用精密星历可以对低轨卫星进行高精度定轨，采用广播星历则可以实现实时导航定位。该方法，不受导航定位过程中换星影响，可靠度高。可对地面、海面以及空中静止的或运动的载体进行的高精度的定位导航。单频GPS接收机，具有成本低廉的优势，利于在民用高精度定位导航方面的普及与应用。