CN109668562A - 一种考虑偏差时引入伪测量的重力梯度运动学导航方法 - Google Patents

Abstract

一种考虑偏差时引入伪测量的重力梯度运动学导航方法，属于自主导航领域。其步骤如下：一：准备工作：给出重力梯度相关的计算公式，以及重力梯度在不同坐标系下的表示；二：设置解算初始参数，准备数据，例如迭代初值、测量方差的选取等；三：利用EGM2008模型计算重力梯度值和重力梯度对位置的偏导阵；四：引入伪测量，构建运动学导航观测模型；五：利用非线性最小二乘方法进行迭代求解。通过以上步骤，引入伪测量，增加了观测信息，然后构建考虑测量偏差的重力梯度序贯运动学导航观测方程。摆脱了传统卫星定轨应用中对卫星动力学约束的依赖，将重力梯度导航的应用扩展到一般物体，从而大大提高了重力梯度定位导航的应用潜力和应用价值。

Description

一种考虑偏差时引入伪测量的重力梯度运动学导航方法

技术领域

本发明提供一种考虑偏差时引入伪测量的重力梯度运动学导航方法，它涉及一种基于重力梯度测量的序贯运动学定位方法和偏差固定方法，即一种对无动力学约束的静止或者运动载体利用重力梯度测量进行位置估计和测量偏差估计的方法，属于导航技术领域。

背景技术

全张量重力梯度匹配导航技术是一种新兴的自主导航方式。其利用地球本身的重力梯度信息，通过建模、匹配等方法对载体位置进行估计。重力梯度匹配导航技术具有实时、全自主、抗干扰等诸多优势，通过在载体平台上安装一个重力梯度仪作为测量仪器，得到六个输出量，如果不考虑仪器的测量偏差，则一个时刻的测量数据可以直接解算出位置信息。由于重力梯度匹配导航技术具有完全自主、抗干扰、无误差积累等优点，在未来有望成为新一代的主流导航系统。

重力梯度是重力势对位置的二阶导数，是一个3阶对称方阵，其中只有五个元素是相互独立的。在通常的重力梯度辅助导航系统或者作为单独的导航系统中，一般不考虑重力梯度仪的测量偏差。在重力梯度导航方式作为单独工作的导航系统中，目前已经有融合卫星动力学信息的导航方式出现。重力梯度定位精度依赖于重力梯度仪的测量精度和重力场模型的精度，同时依赖于是否有其它信息的约束，比如动力学信息。依据是否利用卫星动力学信息，可以将重力梯度导航方式分为动力学导航和运动学导航。其中运动学导航算法因为不需要动力学信息的约束，适用范围更广，可以为车辆、舰船、飞机、火箭和卫星等提供导航信息。在不考虑仪器测量偏差的情况下，运动学导航算法可以利用六个重力梯度测量量来估计三个位置信息，根据仪器测量数据和相对应的测量方差，利用迭代非线性最小二乘解算得到位置信息。此种方法相当于对于每一时刻进行单点定位解算。在考虑仪器测量偏差时，由于仪器测量偏差未知，因此需要将仪器测量偏差当作未知量进行估计。这样就增加了未知量的个数，有九个未知量，而观测方程还是只有六个，因此观测方程秩亏，无法解算。由于考虑重力梯度仪测量偏差时运动学解算方程秩亏的问题，目前重力梯度导航的研究主要有两方面，第一是作为辅助导航方式应用在舰船、潜艇等载体上，第二是结合卫星动力学信息对卫星进行导航。

本发明在运动学导航方法的基础上提出了一种考虑测量偏差的序贯运动学导航方法，该方法引入了伪测量，从而增加了每个历元的观测信息，解决了普通运动学导航观测方程秩亏的问题，而且可以同时对载体位置和仪器测量偏差进行估计。由于是运动学导航方法，因此该方法可以摆脱对卫星动力学信息的依赖，能够用于普通载体如汽车、火箭、飞机等的导航。

发明内容

(一)发明目的：本发明以序贯运动学导航算法基础，针对考虑带有偏差的运动学导航中秩亏的问题，提出了一种引入伪测量的序贯运动学导航方法。

(二)技术方案

本发明提出一种考虑偏差时引入伪测量的重力梯度运动学导航方法，其步骤如下：

步骤一：准备工作

首先，重力梯度定义为重力势函数对位置的二阶导数，以U表示地球重力势，其表达形式如下：

式中GM表示地球的引力常数，表示地球半径，C_nm和S_nm表示重力场模型中的球谐项系数，其值由重力场模型给出。m和n表示重力场模型的阶次数(order，degree)。V_nm和W_nm的计算可参见Montenbruck所著《Satellite Orbits》一书。现总结如下：

式中，r表示位置[x，y，z]处与地球球心之间的距离，根据上面式(2)和式(3)，可以分别迭代计算各个阶次的V_nm和W_nm，给出零阶零次的V₀₀和W₀₀，

其递推计算过程如下：

重力梯度张量是重力势对位置的二阶导数，单位为E(1E＝10^-9s^-2)，其张量矩阵在坐标系ECEF下的表达式如下：

为V_nm，W_nm相对于位置的二阶偏导，式(7)为V_nm对位置的二阶导数，式(8)为W_nm对位置的二阶导数。GGT^E为ECEF坐标系下的重力梯度张量矩阵。

由上面(2)式到(8)式即可计算出在地球固连坐标系ECEF下某一位置处的重力梯度。重力梯度仪的的测量值是在重力梯度仪坐标系(GRF)下表示的，在实际定位中涉及到重力梯度在不同坐标系下的表示，下面给出重力梯度在地球固连坐标系(ECEF)和重力梯度仪坐标系(GRF)之间的转换关系：

其中T_ecef2grf表示ECEF到GRF坐标系的转换矩阵，在实际中可以由其他手段如星敏感器获得。GGT^E表示ECEF坐标系下的重力梯度张量矩阵。T_ecef2grf的表达形式如下：

c_ij表示矩阵第i行第j列的元素。考虑到重力梯度张量矩阵是一个对称矩阵，重力梯度仪实际输出是张量矩阵GGT中上三角的6个测量结果，组成一个包含6个元素的向量。其表达式为：

式中V^G表示重力梯度在重力梯度仪坐标系(GRF)下的分量。

实际重力梯度仪测量中，包含较大偏差，给出考虑测量偏差时的重力梯度原始观测方程如下：

式中表示GRF坐标系下的观测值，HT表示重力梯度从ECEF到GRF的转换矩阵，是一个6×6的矩阵，矩阵中的元素c_ij即是T_ecef2grf中第i行，第j列的元素。表示重力梯度在ECEF坐标系位置r处的真实值，b表示重力梯度仪的测量偏差，w表示重力梯度仪的测量噪声。重力梯度与位置之间是强非线性关系，在利用重力梯度求解位置信息时需将原始观测方程线性化，为统一计算将测量方程在[r₀，b₀]处线性化，过程如下：

将上式整理得到考虑测量偏差的重力梯度线性化观测方程如下：

其中表示在位置r＝r₀处重力梯度向量对位置向量的偏导阵。

步骤二：初始参数设置

给出历元k时的位置初值r_k0和测量偏差初值b_k0，可将上一历元的估计结果作为本历元的位置初值r_k0和偏差初值b_k0：

以及实际观测测量方差阵R_rk和偏差观测测量方差阵R_bk。实际观测测量方差取值为GOCE星载重力梯度仪测量方差：

R_rk＝diag([15 15 15 350 15 500])² (17)

且实际测量方差不随时间变化。由于测量偏差是伪测量，不存在实际的测量，也就不存在测量方差，利用测量偏差不随时间变化的特性，由于将k-1历元的估计值作为本历元的测量值，故可将k-1历元的偏差估计的协方差P_b(k-1)作为本历元的测量方差R_bk：

R_bk＝P_b(k-1) (18)

上面的位置初值、偏差初值和偏差测量方差阵都牵涉到上一历元的信息，对于第一个历元，由于没有上一历元信息，需要给出先验信息，假定位置初值r₁₀、偏差初值b₁₀，以及偏差的测量方差R_b1。第一历元的信息是任意给出的，存在很大不确定性，因此第一时刻的测量方差阵需要给出一个稍大一些的值，例如第一时刻六个偏差的初值误差都是10⁴mE，则取测量噪声都为10⁴mE量级，即：

R_b1＝diag([10⁴ 10⁴ 10⁴ 10⁴ 10⁴ 10⁴]) (19)

然后，将历元k时ECEF坐标系与GRF坐标系之间的转换矩阵T_ecef2grf提取出来，并根据式(13)构造出相应的HT_k，在实际测量中，T_ecef2grf可以由其它手段获得，如星敏等设备。

步骤三：计算位置r_k0处的重力梯度值GGV(r_k0)和偏导阵Hg(r_k0)

由式(2)到式(8)，可以求出在ECEF坐标系下位置r_k0处的重力梯度张量矩阵，由于矩阵是对称的，因此，只需取下(上)三角元素即可，得到GGV(r_k0)如下：

GGV(r_k0)＝GGT^E[1，5，9，2，3，6]＝[GGV_xx，GGV_yy，GGV_zz，GGV_xy，GGV_xx，GGV_yz] (20)

同时，还需要求出在位置r_k0处的重力梯度GGV(r_k0)对位置的偏导阵Hg(r_k0)，经推导，GGV第一个元素对位置的偏导表达形式如下：

其中，

同理，对于重力梯度向量的第二个元素GGV_yy，将其对位置求偏导，我们得到：

其中，

对于重力梯度向量的第三个元素，将其对位置求偏导，我们得到：

其中，

对于重力梯度的其余元素，其对位置的偏导有以下关系：

其中，

至此，我们已经求得重力梯度对位置的偏导阵Hg，其表达形式如下：

式中所有的元素都已在上面求出。

步骤四：引入伪测量，构建运动学导航观测模型

重力梯度只是位置的函数，对于考虑偏差时的运动学导航，其原始观测方程与其线性化形式都已经在步骤一中给出，见式(12)及式(15)。

由于方程只有六个观测量，分别为重力梯度仪的六个输出量，而待估量有9个，分别为三个位置待估量x，y，z，六个偏差待估量b_xx，b_yy，b_zz，b_xy，b_xz，b_yz，因此，观测量的个数少于待估量的个数，观测方程秩亏，无法解算。为解决方程秩亏的问题，本方法在重力梯度运动学导航方法中引入了伪测量的概念，考虑到每个历元的测量偏差都相等，因此可将上一历元的估计结果作为本历元的观测值，其具体表达式如下：

其中b_k为本历元k的偏差真值，w为测量噪声。上式本就已经是线性方程，为计算上的统一性，将其考虑成非线性方程，对其在b_k＝b_k0处进行线性化，可以得到伪观测方程线性化后的形式如下：

在步骤二中已经说明，由于测量偏差b_k不随时间变化，因此上一历元的估计协方差P_b(k-1)即可作为本历元的观测方差R_bk：

R_bk＝P_b(k-1) (30)

由于每个历元偏差的测量方差都是由上一历元的估计协方差传递而来，因此，每个历元偏差的测量方差都是不断变化的。结合式(15)和式(29)，可以得到考虑测量偏差的重力梯度序贯运动学导航的观测方程如下：

其中，令

得到线性化表达式为：

z_k-Z_k＝H_k·dX+w (32)

步骤五：非线性最小二乘方法迭代求解

根据线性化之后的观测方程式(32)，由线性最小二乘理论可以直接得到每次迭代修正量的大小dX，

得到修正量之后，更新待估量，表达式如下：

X_k＝X_k0+dX (34)

此时需要判断修正量的模是否超过预定阈值，此处我们设定阈值为1，即：

norm(dX)＞1？ (35)

如果判断为是，即修正量大于阈值，则更新后的位置和偏差作为初值r_k0和b_k0，即：

然后重新进入步骤三。

如果判断为否，即修正量足够小，则可以结束本历元的解算，然后保存本历元数据，作为下一历元的初值，进入步骤二进行下一历元的解算。解算结果的协方差矩阵为：

此处需要指出的是，步骤二中对于偏差伪测量的测量方差，需要取前一历元的估计方差，此处所选协方差方法如下：

R_k＝P_k-1(4：9，4：9) (38)

通过以上步骤，引入伪测量，增加了观测信息，然后构建考虑测量偏差的重力梯度序贯运动学导航观测方程。摆脱了传统卫星定轨应用中对卫星动力学约束的依赖，将重力梯度导航的应用扩展到一般物体，从而大大提高了重力梯度定位导航的应用潜力和应用价值。

(三)优点

本发明提供的一种考虑测量偏差时引入伪测量的重力梯度序贯运动学导航方法的优点在于：

①本发明对于不需要卫星动力学信息，仅仅依靠重力梯度仪的测量输出就能进行定位，尤其适用于普通物体的定位导航。

②本发明在解决了传统运动学定位无法解决测量偏差和位置同时估计的问题。

附图说明

说明书附图说明了本发明一种考虑测量偏差时引入伪测量的重力梯度序贯运动学导航方法的具体实施过程。

具体实施方式

本发明一种考虑测量偏差时引入伪测量的重力梯度序贯运动学导航方法，见附图中的图1所示，其具体实施步骤如下：

步骤一：初始参数设置

给出历元k时的位置初值r_k0和测量偏差初值b_k0，可将上一历元的估计结果作为本历元的位置初值r_k0和偏差初值b_k0：

以及实际观测测量方差阵R_rk和偏差观测测量方差阵R_bk。实际观测测量方差取值为GOCE星载重力梯度仪测量方差：

R_rk＝diag([15 15 15 350 15 500])² (40)

且实际测量方差不随时间变化。由于测量偏差是伪测量，不存在实际的测量，也就不存在测量方差，利用测量偏差不随时间变化的特性，由于将k-1历元的估计值作为本历元的测量值，故可将k-1历元的偏差估计的协方差P_b(k-1)作为本历元的测量方差R_bk：

R_bk＝P_b(k-1) (41)

上面的位置初值、偏差初值和偏差测量方差阵都牵涉到上一历元的信息，对于第一个历元，由于没有上一历元信息，需要给出先验信息，假定位置初值r₁₀、偏差初值b₁₀，以及偏差的测量方差R_b1。第一历元的信息是任意给出的，存在很大不确定性，因此第一时刻的测量方差阵需要给出一个稍大一些的值，例如第一时刻六个偏差的初值误差都是10⁴mE，则取测量噪声都为10⁴mE量级，即：

R_b1＝diag([10⁴ 10⁴ 10⁴ 10⁴ 10⁴ 10⁴]) (42)

然后，将历元k时ECEF坐标系与GRF坐标系之间的转换矩阵T_ecef2grf提取出来，并根据式(51)构造出相应的HT_k，在实际测量中，T_ecef2grf可以由其它手段获得，如星敏等设备。此部分对应于说明书附图中的1～2.

步骤二：计算位置r_k0处的重力梯度值GGV(r_k0)和偏导阵Hg(r_k0)

由式(2)到式(8)，可以求出在ECEF坐标系下位置r_k0处的重力梯度张量矩阵，由于矩阵是对称的，因此，只需取下(上)三角元素即可，得到GGV(r_k0)如下：

GGV(r_k0)＝GGT^E[1，5，9，2，3，6]＝[GGV_xx，GGV_yy，GGV_zz，GGV_xy，GGV_xz，GGV_yz] (43)

同时，还需要求出在位置r_k0处的重力梯度GGV(r_k0)对位置的偏导阵Hg(r_k0)，经推导，GGV第一个元素对位置的偏导表达形式如下：

其中，

同理，对于重力梯度向量的第二个元素GGV_yy，将其对位置求偏导，我们得到：

其中，

对于重力梯度向量的第三个元素，将其对位置求偏导，我们得到：

其中，

对于重力梯度的其余元素，其对位置的偏导有以下关系：

其中，

至此，我们已经求得重力梯度对位置的偏导阵Hg，其表达形式如下：

式中所有的元素都已在上面求出。此部分对应于说明书附图中的框图3.

步骤三：引入伪测量，构建运动学导航观测模型

重力梯度只是位置的函数，对于考虑偏差时的运动学导航，其原始观测方程与其线性化形式都已经在步骤一中给出，见式(12)及式(15)。

由于方程只有六个观测量，分别为重力梯度仪的六个输出量，而待估量有9个，分别为三个位置待估量x，y，z，六个偏差待估量b_xx，b_yy，b_zz，b_xy，b_xz，b_yz，因此，观测量的个数少于待估量的个数，观测方程秩亏，无法解算。为解决方程秩亏的问题，本方法在重力梯度运动学导航方法中引入了伪测量的概念，考虑到每个历元的测量偏差都相等，因此可将上一历元的估计结果作为本历元的观测值，其具体表达式如下：

其中b_k为本历元k的偏差真值，w为测量噪声。上式本就已经是线性方程，为计算上的统一性，将其考虑成非线性方程，对其在b_k＝b_k0处进行线性化，可以得到伪观测方程线性化后的形式如下：

在步骤二中已经说明，由于测量偏差b_k不随时间变化，因此上一历元的估计协方差P_b(k-1)即可作为本历元的观测方差R_bk：

R_bk＝P_b(k-1) (53)

由于每个历元偏差的测量方差都是由上一历元的估计协方差传递而来，因此，每个历元偏差的测量方差都是不断变化的。结合式(15)和式(29)，可以得到考虑测量偏差的重力梯度序贯运动学导航的观测方程如下：

其中，令

得到线性化表达式为：

z_k-Z_k＝H_k·dX+w (55)

此部分对应于说明书附图中的框图4.

步骤五：非线性最小二乘方法迭代求解

根据线性化之后的观测方程式(55)，由线性最小二乘理论可以直接得到每次迭代修正量的大小dX，

得到修正量之后，更新待估量，表达式如下：

X_k＝X_k0+dX (57)

此时需要判断修正量的模是否超过预定阈值，即：

norm(dX)＞Threshold？ (58)

如果判断为是，即修正量大于阈值，则更新后的位置和偏差作为初值r_k0和b_k0，即：

然后重新进入步骤三。

如果判断为否，即修正量足够小，则可以结束本历元的解算，然后保存本历元数据，作为下一历元的初值，进入步骤二进行下一历元的解算。解算结果的协方差矩阵为：

此处需要指出的是，步骤二中对于偏差伪测量的测量方差，需要取前一历元的估计方差，此处所选协方差方法如下：

R_k＝P_k-1(4：9，4：9) (61)

此部分对应于说明书附图中的5～8.

通过上述步骤，提出了一种考虑偏差时引入伪测量的重力梯度运动学导航方法。该方法基于重力梯度，考虑重力梯度仪的测量偏差，通过引入伪观测，解决了普通运动学导航方法中的方程秩亏问题，能够同时估计观测偏差和位置。该方法为全自主运动学导航方法，不受外界干扰，可对地面、海面以及空中静止或者运动的载体进行定位导航。

Claims (4)

1.一种考虑偏差时引入伪测量的重力梯度运动学导航方法，其特征在于：首先以重力梯度作为单独的导航系统，并将其应用到航天器导航上；其次，仅仅将重力梯度值为测量值，结合高精度重力场模型，通过模型解算重力梯度；并考虑重力梯度仪所带有的测量偏差，利用运动学序贯方法估计出载体所处的位置信息及速度信息；另外，为弥补考虑测量偏差时重力梯度运动学导航观测方程秩亏的问题，方法引入了伪测量的概念，解决了观测方程秩亏的问题，可以实现卫星位置与偏差的同时估计。本方法主要实施步骤为：设置初始参数；采用高精度的EGM2008模型计算地球固联坐标系(ECEF)相应位置处的重力梯度值和重力梯度对位置的偏导阵Hg；为解决考虑测量偏差时同时估计偏差和位置造成的观测方程秩亏的问题，首先在重力梯度运动学实时定位导航中引入了伪测量的概念，利用测量偏差几乎不随时间变化这一特点，将前一历元的测量偏差的估计量作为本历元测量偏差的虚拟观测量，并构建引入伪测量的运动学导航观测方程。最后，对测量方程进行线性化，并通过迭代最小二乘方法进行求解。得到卫星位置和重力梯度仪的测量偏差。

通过以上步骤，引入伪测量，增加了观测信息，然后构建考虑测量偏差的重力梯度序贯运动学导航观测方程。摆脱了传统卫星定轨应用中对卫星动力学约束的依赖，将重力梯度导航的应用扩展到一般物体。

2.权利要求1中所述的一种考虑偏差时引入伪测量的重力梯度运动学导航方法，其特征在于：将运动学导航的方式引入重力梯度导航体系中，不使用卫星所受到的动力学约束信息。考虑测量偏差时重力梯度运动学导航的线性化观测方程可以写成如下形式：其中表示重力梯度仪坐标系(GRF)下测量得到的重力梯度向量，表示在ECEF坐标系下位置r₀处的重力梯度向量，HT表示将重力梯度向量从ECEF坐标系到GRF坐标系下的转换矩阵，大小为6×6，在实际中由其它方式得到。Hg(r₀)表示ECEF坐标系下的重力梯度向量对位置r₀的偏导阵，表达形式为r和r₀表示分别表示在本历元的实际位置和线性化位置，即位置初值。b和b₀分别表示本历元重力梯度仪的测量偏差真值和测量偏差的线性化点，即偏差初值。w表示重力梯度仪真实测量的测量偏差，考虑为满足高斯分布的均值为零的白噪声。作者所采用的GOCE卫星星载重力梯度仪的测量方差量级为R_g＝diag([15 15 15 350 15 500]mE)²，且假设此测量方差对于每一历元都是相同。其中mE表示重力梯度的单位，1mE＝10^-12s^-2。

3.权利要求1中所述的一种考虑偏差时引入伪测量的重力梯度运动学导航方法，其特征在于：在重力梯度运动学导航定位方法中引入了伪测量的概念，并用其构建运动学观测方程。

所述的引入伪测量的方法，是将前一个历元k-1，(k≥2)的偏差估计结果作为本历元i的观测量b_k，其伪测量的观测表达式为同时，引入伪测量之后，还须给出相应的测量方差，以便后面求解。方差按照相同的思想给出，由于把前一历元的偏差估计值作为本历元的观测量，因此将前一历元的偏差的估计协方差阵作为本历元的测量方差。具体形式为：R_bk＝P_b(k-1)，其中R_bk表示本历元k时刻的测量方差，P_b(k-1)表示上一历元k-1时刻偏差的估计方差。这里对于伪测量，每一个历元的测量数据都可以由上一个历元的估计结果得到，只有第一个历元除外。由于第一个历元没有上一历元信息，因此对于第一个历元的伪测量，需要先给出先验信息，假定第一历元位置初值为r₁₀、偏差初值b₁₀，以及偏差的测量方差R_b1。第一历元的信息是任意给出的，存在很大不确定性，因此第一时刻的测量方差阵需要给出一个稍大一些的值，例如第一时刻六个偏差的观测量与真实的偏差相比，误差都是10⁴mE，则取测量噪声都为10⁴mE量级，即：R_b1＝diag([10⁴ 10⁴ 10⁴ 10⁴ 10⁴ 10⁴])，然后将这些测量传递给解算程序即可。

4.权利要求1中所述的一种考虑偏差时引入伪测量的重力梯度运动学导航方法，其特征在于：利用序贯最小二乘方法对运动学导航问题进行求解。通过上面引入伪观测量，每个历元时刻的观测量和观测方差都已经得到，线性化后的观测方程为：

其中，b_k和b_k0分别表示在k历元时的偏差迭代初值，也即线性化点，r_k和r_k0表示历元k时的位置真值和位置初值，即线性化位置，HT_k表示历元k时重力梯度向量从ECEF到GRF坐标系的转换矩阵，0_6×3表示一个6×3的零矩阵，I₆表示一个6阶的单位阵，这里待估量为r_k和b_k。令 则待估量为X_k。得到线性化表达式为：

z_k-Z_k＝H_k·dX+w

dX＝X_k-X_k0

根据最小二乘理论可以计算得到每次迭代修正量其中，然后更新待估量X_k＝X_k0+dX。将此过程进行迭代直至精度到预设阈值，最终得到本历元的位置和偏差的解算结果。