留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

基于联邦卡尔曼滤波的月球车联合定位方法

魏二虎 汤深权 刘建栋 董翠军

魏二虎, 汤深权, 刘建栋, 董翠军. 基于联邦卡尔曼滤波的月球车联合定位方法[J]. 武汉大学学报 ● 信息科学版, 2017, 42(10): 1431-1437. doi: 10.13203/j.whugis20160137
引用本文: 魏二虎, 汤深权, 刘建栋, 董翠军. 基于联邦卡尔曼滤波的月球车联合定位方法[J]. 武汉大学学报 ● 信息科学版, 2017, 42(10): 1431-1437. doi: 10.13203/j.whugis20160137
WEI Erhu, TANG Shenquan, LIU Jiandong, DONG Cuijun. A Method for Joint Positioning of Lunar Roverbased on Federated Kalman Filter[J]. Geomatics and Information Science of Wuhan University, 2017, 42(10): 1431-1437. doi: 10.13203/j.whugis20160137
Citation: WEI Erhu, TANG Shenquan, LIU Jiandong, DONG Cuijun. A Method for Joint Positioning of Lunar Roverbased on Federated Kalman Filter[J]. Geomatics and Information Science of Wuhan University, 2017, 42(10): 1431-1437. doi: 10.13203/j.whugis20160137

基于联邦卡尔曼滤波的月球车联合定位方法

doi: 10.13203/j.whugis20160137
基金项目: 

国家自然科学基金 41374012

详细信息

A Method for Joint Positioning of Lunar Roverbased on Federated Kalman Filter

Funds: 

The National Natural Science Foundation of China 41374012

More Information
图(3) / 表(1)
计量
  • 文章访问数:  1936
  • HTML全文浏览量:  57
  • PDF下载量:  317
  • 被引次数: 0
出版历程
  • 收稿日期:  2016-08-16
  • 刊出日期:  2017-10-05

基于联邦卡尔曼滤波的月球车联合定位方法

doi: 10.13203/j.whugis20160137
    基金项目:

    国家自然科学基金 41374012

    作者简介:

    魏二虎, 博士, 教授, 主要从事空间大地测量和地球动力学研究.ehwei@sgg.whu.edu.cn

    通讯作者: 汤深权, 硕士生.tangshenquan@whu.edu.cn
  • 中图分类号: P228

摘要: 针对月球车定位,设计和推导了一种基于甚长基线干涉测量(very long baseline interferometry,VLBI)和天文导航系统(celestial navigation system,CNS)相结合的月球车联合定位方法,并采用联邦卡尔曼滤波来实现位置信息的最优估计以增强系统的可靠性和容错能力。最后,通过嫦娥-3号(CE-3)实测数据的解算证明了该方法较使用最小二乘法进行联合解算和单独采用VLBI方法定位,可以获得更高的月球车定位精度。同时,也有效地保障了月球车定位的可靠性和稳定性。

English Abstract

魏二虎, 汤深权, 刘建栋, 董翠军. 基于联邦卡尔曼滤波的月球车联合定位方法[J]. 武汉大学学报 ● 信息科学版, 2017, 42(10): 1431-1437. doi: 10.13203/j.whugis20160137
引用本文: 魏二虎, 汤深权, 刘建栋, 董翠军. 基于联邦卡尔曼滤波的月球车联合定位方法[J]. 武汉大学学报 ● 信息科学版, 2017, 42(10): 1431-1437. doi: 10.13203/j.whugis20160137
WEI Erhu, TANG Shenquan, LIU Jiandong, DONG Cuijun. A Method for Joint Positioning of Lunar Roverbased on Federated Kalman Filter[J]. Geomatics and Information Science of Wuhan University, 2017, 42(10): 1431-1437. doi: 10.13203/j.whugis20160137
Citation: WEI Erhu, TANG Shenquan, LIU Jiandong, DONG Cuijun. A Method for Joint Positioning of Lunar Roverbased on Federated Kalman Filter[J]. Geomatics and Information Science of Wuhan University, 2017, 42(10): 1431-1437. doi: 10.13203/j.whugis20160137
  • 嫦娥三号(CE-3) 月球探测器于2013年12月14日21时11分在月球正面的虹湾地区软着陆成功,随后着陆器和巡视器进行分离[1]。巡视器月面移动期间,中国甚长基线干涉测量网(vevy long baseline interferometry, VLBI)对整个过程进行了跟踪观测。月球车在此期间要提高自身的生存能力以及顺利地开展月面相应的科学任务,其关键是要获取月球车的位置[2-3]。VLBI技术历经50余年的发展,其已经具备高精度、高分辨率等特点,NASA、ESA等机构也早已将VLBI技术成功应用于深空探测工程中。天文导航系统(celestial navigation system, CNS)作为传统导航方法,其具有自主性强、稳定性高、导航误差不随时间积累等特点[4-5],也早已被应用于卫星、军事武器等设备的导航中。

    针对月球车定位,国内外学者进行了大量的研究,提出了同波束VLBI相对定位、航位推算、天文导航和视觉导航等方法。文献[4-5]详细介绍了月球车自主天文导航的原理与数学模型。文献[6]利用天文导航和惯性导航相结合的方式进行了月球车自主导航的研究。文献[7]给出了月球车视觉定位的方法并成功将其应用于我国的嫦娥工程中。文献[8]推导了SBI数据和USB数据联合定位的算法,但其存在一定的局限性:在不具备测距、测速条件时无法使用。文献[9-10]利用统计定位方法确定着陆器和巡视器的坐标,其结果和NASA探测器LRO的航拍图像结果差异小于50 m。文献[11-12]详细论述了VLBI差分相时延在深空探测器定轨、定位中的应用,其得到的事后相对定位结果能够达到米级的精度。

    本文针对月球车定位方法中天文导航方法不受时间、距离长短的影响,能够稳定地提供位置、姿态信息,但是其短时定位精度较差[13-14];VLBI定位精度较高,但是由于外部不可控因素的影响,VLBI信号有时无法接收,非自主性易受干扰的固有缺陷对月球车定位的稳定性产生一定的影响。针对这一问题,本文设计和推导了VLBI联合天文导航的月球车定位算法,利用嫦娥三号VLBI实测数据和天文导航仿真数据实例解算了月球车单点位置。采用联邦卡尔曼滤波实现了VLBI和CNS两个子滤波系统的最优估计和信息融合。通过解算结果可以看出联合定位方法较单一方法的定位精度有所提升;并且采用联邦卡尔曼滤波方法进行联合解算较采用最小二乘方法进行解算,精度也有所提升。

    • AB两点表示的是地面上两个VLBI站点,其测站ITRF2000坐标转换到月固坐标系下的坐标分别为(x1, y1, z1)、(x2, y2, z2)。设月球车信号到达A, B两个站的时间分别为t1t2,时间差为τ12。月球车与VLBI站之间的距离分别为r1r2。VLBI观测定位模型如图 1所示。

      图  1  VLBI定位示意图

      Figure 1.  Schematic of VLBI Positioning

      $$ \begin{array}{l} c{\tau _{12}} = c({t_2}-{t_1}) = {r_2}-{r_1} = \\ \sqrt {{{({x_2}-{x_s})}^2} + {{({y_2} - {y_s})}^2} + {{({z_2} - {z_s})}^2}} - \\ \sqrt {{{({x_2} - {x_s})}^2} + {{({y_1} - {y_s})}^2} + {{({z_1} - {z_s})}^2}} \end{array} $$ (1)

      式中, c表示光速;(xs, ys, zs)表示月球车的月固坐标并将其设为系统状态量$ {\hat X} $。上式为非线性观测方程,将上述公式围绕滤波值$ {\hat X(k + 1/k)}$处展成泰勒级数并略去二阶及以上项线性化得:

      $$ \begin{array}{l} c{\tau _{12}} = ({r_2}-{r_1}){|_{\hat X(k + 1/k)}} + \frac{{\partial (c{\tau _{12}})}}{{\partial {x_s}}}{|_{\hat X(k + 1/k)}}\\ ({{\hat x}_s}\left( {k + 1} \right)-{{\hat x}_s}\left( {k + 1/k} \right)) + \\ \frac{{\partial (c{\tau _{12}})}}{{\partial {y_s}}}{|_{\hat X(k + 1/k)}}({{\hat y}_s}\left( {k + 1} \right)-{{\hat y}_s}\left( {k + 1/k} \right)) + \\ \frac{{\partial (c{\tau _{12}})}}{{\partial {z_s}}}{|_{\hat X(k + 1/k)}}({{\hat z}_s}\left( {k + 1} \right) - {{\hat z}_s}\left( {k + 1/k} \right)) \end{array} $$ (2)

      式中:

      $$ \left\{ \begin{array}{l} \frac{{\partial (c{\tau _{12}})}}{{\partial {x_s}}} = \frac{{{x_1}-{x_s}}}{{{r_1}}}-\frac{{{x_2}-{x_s}}}{{{r_2}}}\\ \frac{{\partial (c{\tau _{12}})}}{{\partial {y_s}}} = \frac{{{y_1} - {y_s}}}{{{r_1}}} - \frac{{{y_2} - {y_s}}}{{{r_2}}}\\ \frac{{\partial (c{\tau _{12}})}}{{\partial {z_s}}} = \frac{{{z_1} - {z_s}}}{{{r_1}}} - \frac{{{z_2} - {z_s}}}{{{r_2}}} \end{array} \right. $$

      利用中国VLBI网的4个站点,可以组成6条观测基线。但是仅仅使用VLBI得到的时延信息来转换成角位置信息时,会有较大的误差。因此,在实际解算月球车位置时往往需要加入月球车的距离约束条件。可以建立如下的条件约束方程:

      $$ \sqrt {x_{_s}^{^2} + y_{_s}^{^2} + z_{_s}^{^2}}-\sqrt {x_{_{s0}}^{^2} + y_{_{s0}}^{^2} + z_{_{s0}}^{^2}} = 0 $$ (3)

      (xs0, ys0, zs0)表示为t0时刻的月球车的坐标初始值,同理其线性化公式为:

      $$ \begin{array}{l} \frac{{{x_{S0}}}}{{\sqrt {x_{_{S0}}^{^2} + y_{_{S0}}^{^2} + z_{_{S0}}^{^2}} }}({{\hat x}_s}\left( {k + 1} \right)-\\ {{\hat x}_s}\left( {k + 1/k} \right)) + \frac{{{y_{S0}}}}{{\sqrt {x_{_{S0}}^{^2} + y_{_{S0}}^{^2} + z_{_{S0}}^{^2}} }}\cdot\\ ({{\hat y}_s}\left( {k + 1} \right)-{{\hat y}_s}\left( {k + 1/k} \right)) + \frac{{{z_{S0}}}}{{\sqrt {x_{_{S0}}^{^2} + y_{_{S0}}^{^2} + z_{_{S0}}^{^2}} }}\\ ({{\hat z}_s}\left( {k + 1} \right)-{{\hat z}_s}\left( {k + 1/k} \right)) = 0 \end{array} $$ (4)

      以时延值τ为观测量,则6条基线的观测方程和式(4) 联合起来就组成了VLBI定位的观测方程:

      $$ {Z_1}\left( t \right) = \left[{\tau \left( t \right)} \right] = {h_1}\left( t \right)d\hat X + {V_1}\left( t \right) $$ (5)

      式中:

      $$ \begin{array}{l} {h_1}\left( t \right) = \left[ {\begin{array}{*{20}{c}} {\frac{{\partial (c{\tau _{12}})}}{{\partial {x_s}}}}&{\frac{{\partial (c{\tau _{12}})}}{{\partial {y_s}}}}&{\frac{{\partial (c{\tau _{12}})}}{{\partial {z_s}}}}\\ {\frac{{\partial (c{\tau _{13}})}}{{\partial {x_s}}}}&{\frac{{\partial (c{\tau _{13}})}}{{\partial {y_s}}}}&{\frac{{\partial (c{\tau _{13}})}}{{\partial {z_s}}}}\\ {\frac{{\partial (c{\tau _{14}})}}{{\partial {x_s}}}}&{\frac{{\partial (c{\tau _{14}})}}{{\partial {y_s}}}}&{\frac{{\partial (c{\tau _{14}})}}{{\partial {z_s}}}}\\ {\frac{{\partial (c{\tau _{23}})}}{{\partial {x_s}}}}&{\frac{{\partial (c{\tau _{23}})}}{{\partial {y_s}}}}&{\frac{{\partial (c{\tau _{23}})}}{{\partial {z_s}}}}\\ {\frac{{\partial (c{\tau _{24}})}}{{\partial {x_s}}}}&{\frac{{\partial (c{\tau _{24}})}}{{\partial {y_s}}}}&{\frac{{\partial (c{\tau _{24}})}}{{\partial {z_s}}}}\\ {\frac{{\partial (c{\tau _{34}})}}{{\partial {x_s}}}}&{\frac{{\partial (c{\tau _{34}})}}{{\partial {y_s}}}}&{\frac{{\partial (c{\tau _{34}})}}{{\partial {z_s}}}}\\ {{h_{71}}}&{{h_{72}}}&{{h_{73}}} \end{array}} \right],\\ \begin{array}{*{20}{c}} \begin{array}{l} {h_{71}} = \sqrt {{x_{S0}}x_{_{S0}}^{^2} + y_{_{S0}}^{^2} + z_{_{S0}}^{^2}} \\ {h_{72}} = \sqrt {{y_{S0}}x_{_{S0}}^{^2} + y_{_{S0}}^{^2} + z_{_{S0}}^{^2}} \\ {h_{73}} = \sqrt {{z_{S0}}x_{_{S0}}^{^2} + y_{_{S0}}^{^2} + z_{_{S0}}^{^2}} \end{array}\\ {} \end{array}, \end{array} $$

      V1(t)为VLBI系统的测量噪声。

    • 设月球车位于J2000.0历元时刻的月固坐标系下赤经为λ,赤纬为φ。天文导航的观测方程[4]为:

      $$ \left\{ \begin{array}{l} h = {\rm{arcsin}}\left[{{\rm{sin}}\varphi {\rm{sin}}\delta + {\rm{cos}}\varphi {\rm{cos}}\delta {\rm{cos}}\left( {\lambda + GHA-\alpha } \right)} \right]\\ A = {\rm{arctan}}\left[{\frac{{{\rm{sin}}\left( {\lambda + GHA-\alpha } \right)}}{{{\rm{cos}}\left( {\lambda-\alpha } \right){\rm{sin}}\varphi-{\rm{cos}}\varphi {\rm{tan}}\delta }}} \right] \end{array} \right. $$ (6)

      式中, hA分别表示的是目标天体在月球车所在地的水平坐标系中的高度角和方位角;αδ分别表示的是目标天体的赤经、赤纬;GHA表示的是春分点的格林尼治时角,这3个参数可以从星表或给定的计算公式中获得。为方便计算,本文将sinh、tanA当作“虚拟观测值”并对上述观测方程进行线性化,得到如下结果:

      $$ \left\{ \begin{array}{l} {\rm{sin}}h = {\left( {{\rm{sin}}\;h} \right)_0} + \frac{{\partial {\rm{sin}}h}}{{\partial \varphi }}d\varphi + \frac{{\partial {\rm{sin}}h}}{{\partial \lambda }}d\lambda + 0\;{\rm{d}}H\\ {\rm{tan}}A = {\left( {{\rm{tan}}\;A} \right)_0} + \frac{{\partial {\rm{tan}}A}}{{\partial \varphi }}d\varphi + \frac{{\partial {\rm{tan}}A}}{{\partial \lambda }}d\lambda + 0\;{\rm{d}}H \end{array} \right. $$ (7)

      式中,(sinh)0、(tanA)0分别表示的是观测值的近似值,其他的系数分别为:

      $$ \left\{ \begin{array}{l} \frac{{\partial {\rm{sin}}h}}{{\partial \varphi }} = {\rm{cos}}\varphi \cdot{\rm{sin}}\delta-{\rm{sin}}\varphi \cdot{\rm{cos}}\delta \cdot{\rm{cos}}\left( {LHA} \right)\\ \frac{{\partial {\rm{sin}}h}}{{\partial \lambda }} =-{\rm{cos}}\varphi \cdot{\rm{cos}}\delta {\rm{sin}}\left( {LHA} \right)\\ \frac{{\partial {\rm{tan}}A}}{{\partial \varphi }} = {\rm{sin}}\left( {LHA} \right)\cdot-\left( {\frac{{{\rm{cos}}\varphi \cdot{\rm{cos}}\left( {LHA} \right) + {\rm{sin}}\varphi {\rm{{\cdot}tan}}\delta }}{{{{({\rm{sin}}\varphi \cdot{\rm{cos}}\left( {LHA} \right) - {\rm{cos}}\varphi \cdot{\rm{tan}}\delta )}^2}}}} \right)\\ \frac{{\partial {\rm{tan}}A}}{{\partial \lambda }} = \frac{{{\rm{cos}}\left( {LHA} \right)\cdot({\rm{sin}}\varphi \cdot{\rm{cos}}\left( {LHA} \right) - {\rm{cos}}\varphi {\rm{{\cdot}tan}}\delta )}}{{{{({\rm{sin}}\varphi {\rm{{\cdot}cos}}\left( {LHA} \right) - {\rm{cos}}\varphi \cdot{\rm{tan}}\delta )}^2}}} + \\ \frac{{{\rm{sin}}\left( {LHA} \right)\cdot{\rm{sin}}\varphi \cdot{\rm{sin}}\left( {LHA} \right)}}{{{{({\rm{sin}}\varphi \cdot{\rm{cos}}\left( {LHA} \right) - {\rm{cos}}\varphi {\rm{{\cdot}tan}}\delta )}^2}}}\\ LHA = GHA + \lambda - \alpha \end{array} \right. $$ (8)

      式中,LHA表示的是目标天体的地方时角、H表示的是月面高程。在天文导航观测值仿真的过程中,以太阳和地球为目标观测天体。则观测方程的系数矩阵BCNS为:

      $$ {\mathit{\boldsymbol{B}}_{{\rm{CNS}}}} = \left[{\begin{array}{*{20}{c}} {\frac{{\partial {\rm{sin}}{h_E}}}{{\partial \varphi }}}&{\frac{{\partial {\rm{sin}}{h_E}}}{{\partial \lambda }}}&0\\ {\frac{{\partial {\rm{tan}}{A_E}}}{{\partial \varphi }}}&{\frac{{\partial {\rm{tan}}{A_E}}}{{\partial \lambda }}}&0\\ {\frac{{\partial {\rm{sin}}{h_S}}}{{\partial \varphi }}}&{\frac{{\partial {\rm{sin}}{h_S}}}{{\partial \lambda }}}&0\\ {\frac{{\partial {\rm{tan}}{A_S}}}{{\partial \varphi }}}&{\frac{{\partial {\rm{tan}}{A_S}}}{{\partial \lambda }}}&0 \end{array}} \right] $$

      式中,hEAE分别表示的是地球对应的高度角和方位角;hSAS分别表示的是太阳对应的高度角和方位角。

      分别取太阳敏感器测量精度为6″,地球敏感器测量精度为36″,采样间隔为5 s (同VLBI采样间隔一致)的数据进行天文导航仿真实验。观测值仿真的具体流程如下:

      1) 通过JPL星历(本文采用DE423) 获取目标天体相对于月球的三维相对坐标,通过坐标转换,转换成观测星体的赤纬、赤经;

      2) 将通过VLBI解算得到的月球车的月固直角坐标通过坐标转换,转换成月固系下的赤经、赤纬;

      3) 将步骤1)、2) 得到的值代入式(6) 中模拟出天文导航的观测值;

      4) 加入观测误差,如星敏感器误差、星历误差等。

      由于月心直角坐标和月心大地坐标存在如下的函数关系:

      $$ \left\{ \begin{array}{l} {x_s} = \left( {N + H} \right){\rm{cos}}\varphi {\rm{cos}}\lambda \\ {y_s} = \left( {N + H} \right){\rm{cos}}\varphi {\rm{sin}}\lambda \\ {z_s} = \left[{N{{\left( {1-e} \right)}^2} + H} \right]{\rm{sin}}\varphi \end{array} \right. $$ (9)

      式中,N表示的是卯酉圈半径,$ N = \frac{a}{{\sqrt {1-{e^2}{\rm{si}}{{\rm{n}}^2}\varphi } }} $,a为椭球长半轴,e为椭球第一离心率。先对式(9) 求全微分,设B2矩阵为:

      $$ \begin{array}{l} \mathit{\boldsymbol{B}}2 = \left[{\begin{array}{*{20}{c}} {CN{\cdot}{\rm{cos}}\lambda {\rm{{\cdot}cos}}\varphi-{\rm{cos}}\lambda {\rm{sin}}\varphi {\cdot}K}&{-{\rm{sin}}\lambda {\cdot}K{\cdot}{\rm{cos}}\varphi }&{{\rm{cos}}\lambda {\cdot}{\rm{cos}}\varphi }\\ {CN{\cdot}{\rm{cos}}\varphi {\cdot}{\rm{sin}}\lambda-{\rm{sin}}\lambda {\rm{sin}}\varphi {\cdot}K}&{{\rm{cos}}\lambda {\cdot}K{\cdot}{\rm{cos}}\varphi }&{{\rm{cos}}\varphi {\cdot}{\rm{sin}}\lambda }\\ {{\rm{cos}}\varphi N(1 - {e^2}) + H + {\rm{sin}}\varphi (1 - {e^2}){\cdot}CN}&0&{{\rm{sin}}\varphi } \end{array}} \right]\\ CN = a{e^2}{(1 -{e^2}si{n^2}\varphi )^{ -\frac{3}{2}}}{\rm{sin}}\varphi {\rm{cos}}\varphi, K = N + H \end{array} $$ (10)

      (B2)3×3为两组参数的转换矩阵,转换关系为:

      $$ \left( {\begin{array}{*{20}{c}} {d{x_s}}\\ {d{y_s}}\\ {d{z_s}} \end{array}} \right) = {\left( {\mathit{\boldsymbol{B}}2} \right)_{3 \times 3}}\cdot\left( {\begin{array}{*{20}{c}} {d\varphi }\\ {d\lambda }\\ {dH} \end{array}} \right) $$ (11)

      通过B2矩阵,将天文导航的(, )参数转换为(dx, dy, dz)参数,这样VLBI系统观测方程和天文导航观测方程的参数就统一起来了。此后同理将天文导航观测方程围绕滤波值$ \hat X\left( {k + 1/k} \right) $进行展开。以sinh,tanA为观测值,则天文导航的观测方程为:

      $$ {\mathit{\boldsymbol{Z}}_2}\left( t \right) = \left[{\begin{array}{*{20}{c}} {{\rm{sin}}\;h}\\ {{\rm{tan}}A} \end{array}} \right] = {h_2}\left( t \right)d\hat X + {V_2}\left( t \right) $$ (12)

      式中,h2(t)=BCNS·B2;V2(t)为天文导航系统的测量噪声。

    • 在月球车月固坐标求解的过程中,需要将某一观测历元某测站的协议地球坐标转换成相应的月固坐标,该转换过程最后一步是从月心天球坐标系转换到月固坐标系。J2000惯性坐标系到月固坐标系的转换矩阵Mt(t)的计算式为:

      $$ {M_t}\left( t \right) = {R_Z}\left( \mu \right){R_X}\left( i \right){R_Z}\left( \mathit{\Omega} \right) $$ (13)

      式中,RXRZ分别表示的是绕X轴、Z轴旋转的坐标转换矩阵;μiΩ可分别由JPL星历查得。对于单点解算,由于状态量$ {\hat X} $从J2000惯性坐标系转换到月固坐标系后随时间变化保持不变[15],所以其状态转移矩阵可表示为:

      $$ \begin{array}{l} \hat X\left( {k + 1} \right) = {M_t}({t_{k + 1}}){({M_t}({t_k}))^{-1}}\hat X\left( k \right)\\ \;\;\;\;\;\;\;\;\; = \mathit{\Phi} \left( {k + 1, k} \right)\hat X\left( k \right) \end{array} $$ (14)

      本文将VLBI和CNS有机结合起来进行信息融合。在子滤波器中,采用状态方程式(14) 和观测方程式(5) 构成VLBI子滤波系统,采用状态方程式(14) 和观测方程式(12) 构成CNS子滤波系统。并采用联邦卡尔曼滤波技术完成系统状态量的最优估计和信息融合。各子滤波器卡尔曼滤波算法[16-17]为:

      1) 状态一步预测方程:

      $$ {X_i}\left( {k + 1/k} \right) = \mathit{\Phi} \left( {k + 1/k} \right){X_i}\left( k \right) $$ (15)

      2) 状态估值计算方程:

      $$ \begin{array}{l} {X_i}\left( {k + 1} \right) = {X_i}\left( {k + 1/k} \right) + {K_i}\left( {k + 1} \right)\cdot\\ ({Z_i}\left( {k + 1} \right)-{H_i}\left( {k + 1} \right){X_i}\left( {k + 1/k} \right)) \end{array} $$ (16)

      3) 滤波增益方程:

      $$ \begin{array}{l} {K_i}\left( {k + 1} \right) = {P_i}\left( {k + 1/k} \right){H_i}^T\left( {k + 1} \right)\cdot[{H_i}\left( {k + 1} \right)\\ {P_i}\left( {k + 1/k} \right){H^T}_i\left( {k + 1} \right) + {R_i}\left( {k + 1} \right){]^{ -1}} \end{array} $$ (17)

      4) 一步预测均方差方程:

      $$ \begin{array}{l} {P_i}\left( {k + 1/k} \right) = \ \mathit{\Phi} \left( {k + 1/k} \right){P_i}\left( k \right){\mathit{\Phi} ^T}\left( {k + 1/k} \right) + {Q_i}\left( k \right) \end{array} $$ (18)

      5) 估计均方差方程:

      $$ \begin{array}{l} {P_i}\left( {k + 1} \right) = \\ [I-{K_i}\left( {k + 1} \right){H_i}\left( {k + 1} \right)]{P_i}\left( {k + 1/k} \right) \end{array} $$ (19)

      式中,i=1、2,1表示的是VLBI系统子滤波器,2表示的是CNS子滤波器;Ki为增益矩阵;Pi为状态量估计误差方差矩阵;Qi为系统过程噪声协方差矩阵;Ri为系统观测协方差矩阵。

      若令主系统不占任何全局信息,仅对子系统估计进行综合运算。此时系统的全局最优估计状态为:

      $$ \left\{ \begin{array}{l} {X_g}\left( k \right) = {P_g}\left( k \right)[{P^{-1}}_1\left( k \right){X_1}\left( k \right) + {P^{-1}}_2\left( k \right){X_2}\left( k \right)]\\ {P_g}\left( k \right) = {[{P^{-1}}_1\left( k \right) + {P^{-1}}_2\left( k \right)]^{ -1}} \end{array} \right. $$ (20)

      之后将全局最优估计值传给两个子滤波器,作为下一时刻两个子滤波器的估计值(i=1, 2)。

      $$ \left\{ \begin{array}{l} {X_i}\left( k \right) = {X_g}\left( k \right)\\ P_{_i}^{^{-1}}\left( k \right) = P_{_g}^{^{-1}}{\beta _i}\\ Q_{_i}^{^{-1}}\left( k \right) = Q_{_g}^{^{ - 1}}{\beta _i} \end{array} \right. $$ (21)

      式中,β1为VLBI子系统所对应的信息分配因子;β2为天文导航子系统所对应的信息分配因子。信息分配因子的选择满足信息守恒式β1+β2=1(0≤β1, β2≤1) 且其数值大小跟子滤波器的滤波精度成正比。为了使联合系统具有更好的自适应能力,使用基于估计误差矩阵P的范数动态分配因子,其计算公式为[17]

      $$ {\beta _i}\left( k \right) = \frac{{{{({{\left\| {{P_i}\left( {k-1} \right)} \right\|}_F})}^{-1}}}}{{\sum\limits_{i = 1}^2 {{{({{\left\| {{P_i}\left( {k-1} \right)} \right\|}_F})}^{ - 1}}} }} $$ (22)

      式中,Pi(k-1) 为各子滤波器状态量估计误差方差矩阵;‖·‖F表示的是Frobenius范数。对于任意一个矩阵A来说,其Frobenius范数可以表示为$ {\left\| A \right\|_F} = \sqrt {\sum {\rm{diag}}({A^T}{\cdot}A)} $。

    • 本文采用CE-3月球探测器2013年12月20日地面北京(BJ)、昆明(KM)、乌鲁木齐(UR)和上海(TM)4个VLBI站19:41:57.439 125~20:48:32.439 156时段的实测观测值来进行月球车位置的解算。在运用最小二乘算法(LS)进行VLBI和CNS联合解算时,刚开始给定VLBI和CNS两类观测值的权阵是不合理的,给定后再采用赫尔默特验后定权方法[18]实现最终的权值确定。之后通过误差方程求解获得改正数,再将初值和改正数之和作为待求参数,进一步构造误差方程进行解算,直到小于一定的阈值时解算得到最终值。在运用联邦卡尔曼滤波(FKF)解算的过程中,以最小二乘方法解算得到的月球车月固坐标作为状态量的初始值(单位:m)为:$ {{\hat X}_0} $=(1 172 330.9,-416 020.8, 1 208 219.9),Q矩阵可根据JPL星历的精度得到,Q矩阵过大或过小都会对滤波的结果产生较大的影响,本文取(单位:km2):

      Q=diag(10-8,10-8,10-8)

      R1R2由两个观测方程的噪声可以得出。VLBI子滤波系统中(单位:km2):

      P1=diag[1, 1, 1]

      CNS子滤波系统中:

      P2=((0.01°) 2, (0.01°) 2)

      图 2所表示的是采用单一解算和分别采用最小二乘方法、联邦卡尔曼滤波方法联合VLBI和CNS进行解算,所得到的月球车在XYZ 3个方向的坐标中误差。图 2中最下面一幅图表示的是联合定位解算和单一解算三维坐标中误差绝对值之和的差值。从图中可以看出,联合解算较单一解算、三维坐标方向的精度都有所提高。在使用最小二乘方法进行联合解算时,X方向精度提高的最大值为2.0 m,最小值为0.6 m,平均值为1.0 m;Z方向精度提高的最大值为2.1 m,最小值为0.7 m,平均值为1.1 m;Y方向提升的最大值也就只有0.06 m。在使用联邦卡尔曼滤波算法进行联合解算时,XYZ三维方向精度最高提高分别为29.4 m、13.4 m、28.1 m;平均值分别为10.5 m、2. 2 m、9.3 m。

      图  2  联合解算和VLBI单独解算的对比

      Figure 2.  Comparison of Joint Solution and Single Solution of VLBI

      可以看出,采用FKF法进行联合解算的精度要比采用LS法解算高。对于三维坐标中误差绝对值之和的差值,LS法精度提高的平均值为2.1 m,而FKF法精度提高的平均值为22.03 m。这是因为在各个子系统运用卡尔曼滤波进行解算时有初始值约束,在初始值给定较合理时,结果收敛较快,同时FKF法利用信息分配原则来消除各子状态估计的相关性,实现两个子系统的信息综合和最优估计[19]。另外,不管是采用哪种方法进行联合解算,Y方向的精度提升都相对较小。这是因为观测信号的投影大都集中在Y方向,Y方向则反映了较真实的观测精度,故联合解算较单一解算,改变较小。

      在解算经纬度数值时,CNS单独解算和分别采用LS法、FKF法进行联合解算时所得到的经纬度方向的误差如图 3所示。

      图  3  经纬度方向误差对比

      Figure 3.  Comparison of Error in Longitude Direction & Latitudinal Direction

      CNS短时测量精度受测量仪器精度的限制,导致其精度较低,在单独采用CNS进行解算时,经纬度方向的误差较大并且波动也较大。而当采用联合解算的时候,因为加入了VLBI观测量和距离约束条件,较单独采用角度信息进行求解时,使得经纬度方向的误差的波动小很多。这反映了在加入VLBI观测量之后,CNS的解算更加的稳定且解算精度有较大的提升。而且使用FKF法解算比LS法解算精度高,这是因为FKF法能得到更优的局部估计和信息综合。

      在进行实例解算时,同样使用了地面BJ、KM、UR和TM 4个VLBI测站的数据。分别使用单一解算和联邦卡尔曼滤波联合算法计算了地面UR-KM-TM、UR-KM-BJ、UR-BJ-TM和TM-BJ-KM 4条基线的残差平均值和残差闭合值。以UR-KM-TM基线为例,将KM-UR基线残差值减去KM-TM基线残差值,再加上UR-TM基线残差值得到这3台站的残差闭合值。表 1列出了相应的解算结果,可见联合解算较单独解算,4条基线的残差平均值和残差闭合值的精度都有所提高。

      表 1  基线残差平均值和闭合值/m

      Table 1.  The Mean Value and the Closure Value of Baseline Residual/m

      残差 UR-KM-TM UR-KM-BJ UR-BJ-TM TM-BJ-KM
      单一残差平均值 -2.56 0.31 0.67 -2.61
      单一残差闭合差 3.44 0.16 -0.15 -3.42
      联合残差平均值 -2.51 0.32 0.66 -2.57
      联合残差闭合差 3.40 0.16 -0.14 -3.37
    • 本文根据月球车导航的实际应用要求,设计并实现了一种VLBI联合CNS的月球车定位方法。分别介绍了VLBI观测定位模型、CNS观测模型,之后根据联邦卡尔曼滤波原理联合了VLBI和CNS系统。通过实测数据实际解算结果可以得出以下结论。

      1) 联合解算较单一解算,月球车定位的精度得到提升。

      2) 采用联邦卡尔曼滤波进行联合解算较最小二乘方法,精度更高。且组合系统的可靠性更高,系统容错能力更强。

      但该组合导航模型也存在一些不足之处,VLBI定位的优势主要是横向精度较高,但其径向精度较差。如果能将VLBI数据和测距、测速数据结合,则横向、径向精度会更高,二者结合将会有更好的互补性。

参考文献 (19)

目录

    /

    返回文章
    返回