CN103674030A - 基于天文姿态基准保持的垂线偏差动态测量装置和方法 - Google Patents

Abstract

本发明公开了一种基于天文姿态基准保持的垂线偏差动态测量装置及方法，构建了INS/GPS姿态测量子系统和LGU/GPS姿态测量子系统，其中LGU/GPS姿态测量子系统进行姿态更新采用的初始值由星敏感器输出的姿态信息提供。将LGU/GPS姿态测量子系统输出的姿态和INS/GPS姿态测量子系统输出的姿态求差，进而计算垂线偏差。最后修正垂线偏差测量值中的跳变误差，并利用全球重力模型数据修正垂线偏差测量值中的低频误差。本发明系统简单；本发明的方法鲁棒性强，测量结果稳定；降低了GPS精度的需求，相比传统的矢量航空重力仪依赖于高精度差分GPS，本发明只需一般的GPS单点定位即可满足精度要求，因此拓展了测量方法的应用范围。

Description

基于天文姿态基准保持的垂线偏差动态测量装置和方法

技术领域

本发明涉及测绘领域,特别是一种基于天文姿态基准保持的垂线偏差动态测量方法。

背景技术

由于地球形状的不规则、地形的起伏、密度分布不均等因素，实际的重力场与正常重力模型之间存在差异，称之为重力扰动，其垂直分量称为重力异常，水平分量表现为垂线偏差。垂线偏差是大地测量学和空间科学中非常重要的数据。通过垂线偏差可以计算高程异常，确定大地水准面，还可应用于天文大地测量观测数据的归算、弹道计算、空间技术及高精度惯性导航领域。此外，通过研究垂线偏差的时变特性与地震、火山活动之间的关系，可以为地震和火山监测提供新的手段。

Journal of Surveying Engineering（《测绘工程杂志》）2010年第一期中“ModernDetermination of Vertical Deflections Using Digital Zenith Cameras”（利用数字天顶相机的垂线偏差现代测量方法）论文公开了一种利用天顶数字相机测量垂线偏差的方法，其原理是利用高精度的天顶数字相机获取载体相对于参考椭球面的水平姿态角，同时利用高精度电子水平仪测量载体相对于真实大地水准面的水平姿态角，将两者求差可得到垂线偏差。该方法的测量精度极高，可以达到0.1″，但由于只能进行静态测量，测量效率很低，因而对人力、物力及时间的消耗巨大。

美国专利US005112126公开了一种用于矿井内垂线偏差测量的方法。其实施方法为：利用天文经纬仪测量得到本地参考椭圆面的法线方向，即正常重力的方向，采用三个高精度激光陀螺进行姿态解算，用于该法线方向的姿态保持，并将激光陀螺组合体放入矿井中，同时将一个摆仪与激光陀螺组合体固联安装，利用摆仪测量当地真实垂线的方向，该发明设计了一种光学测量装置用于测量真实垂线的方向与正常重力方向之间偏差角，即可得到垂线偏差。该专利的缺点在于不能实现动态测量，且测量精度容易受陀螺测量误差的影响，所设计的用于测量垂线偏差的光学测量装置结构较为复杂，通用性不强。

美国专利US3803916公开了一种重力扰动矢量测量的方案。其实施方法为：利用LORAN（long range electronic navigation system远距电子导航系统）进行定位，测量载体在天文坐标系中的位置；同时利用惯性导航系统进行定位解算，将解算结果和LORAN定位结果做差，该差值反映的是陀螺零偏与垂线偏差对惯性导航系统的影响。在扣除陀螺零偏引起的位置误差后，通过捷联惯性导航的反向滤波技术，利用位置误差反算出垂线偏差。这种专利主要缺陷有两点：一是重力扰动模型建模（包括模型的形式和模型参数辨识）较为困难，当模型存在误差时，所估计得到的重力扰动矢量误差较大；二是导航的水平姿态误差角与垂线偏差存在严重的耦合，垂线偏差将引起导航的水平姿态角误差，同时水平姿态角的误差也会引起严重的水平重力扰动估计误差，这两种状态量在只有位置辅助的条件下是不可分离的。

美国专利US5924056和Journal of Geodesy（《大地测量杂志》）2001年的论文“A newapproach for airborne vector gravimetry using GPS-INS”（利用GPS/INS实现航空矢量重力测量的新方法）中公开了一类矢量重力测量方法。它的基本原理是：利用INS（InertialNavigation System惯性导航系统）和GPS（Global Position System全球定位系统）进行组合导航实现比力测量，利用高精度动态差分GPS测量载体的运动加速度，将比力输出和运动加速度做差即可得到矢量重力扰动，这是目前国外研究中普遍采用的矢量重力测量方案。该方案的主要缺陷在于比力的水平分量测量精度不高，垂线偏差的存在将导致水平姿态角测量误差，进而将天向重力加速度耦合到比力的水平分量上，国内外采用多次迭代求解的方法抑制这一影响，但效果欠佳，没有从本质上解决耦合问题，此外这种方案对GPS精度要求极高，需要建立差分GPS基站，在远洋条件下无法使用。根据目前公开的文献，该类矢量重力测量技术尚未有成熟的产品，其精度和分辨率难以令人满意。

重力梯度仪也是垂线偏差动态测量的有力工具，但目前只有少数国家发展出成熟的产品，如Bell Geospace公司的Air-FTG重力梯度仪、美国专利US5357802公开的一种旋转加速度计重力梯度仪等。由于涉及军事领域，重力梯度仪严格限制出口，并对关键技术严格保密。国内的研究进展较为缓慢，离实用还有较大的差距。

鉴于以上原因，全球只有少数国家和地区具有较为丰富的高精度垂线偏差数据，且这些垂线偏差数据大多数都未公开。近年来，卫星重力测量技术的发展使全球重力模型的精度和分辨率不断提高，促使全球垂线偏差数据精度有了较大提高，并在大地测量、地球物理、惯性导航等领域中得到广泛应用。然而，对于某些特殊领域的应用，如地震和火山的监测、矿藏勘探及超高精度导航等，仍难以满足要求。

针对高精度垂线偏差数据获取困难的问题，迫切需要发展一种能够在动态条件下实现垂线偏差测量的新方法，以更高效地获取高精度、高分辨率的垂线偏差数据，降低测量所消耗的成本和时间，推动其他相关学科领域的发展。

发明内容

本发明所要解决的技术问题是，针对现有技术不足，提供一种基于天文姿态基准保持的垂线偏差动态测量装置及方法，在满足测量精度和分辨率需求的基础上，降低垂线偏差测量所消耗的成本和时间，提高测量效率，同时使测量精度稳定，无测量实施的地域限制。

为解决上述技术问题，本发明所采用的技术方案是：一种基于天文姿态基准保持的垂线偏差动态测量装置，包括由惯性导航系统（INS）、GPS天线、GPS接收机组成的INS/GPS姿态测量子系统、星敏感器，所述惯性导航系统包括三个正交安装的激光陀螺，所述三个正交安装的激光陀螺与GPS接收机组成LGU/GPS姿态测量子系统，三个激光陀螺均与GPS接收机通信；所述惯性导航系统、GPS天线、GPS接收机、星敏感器均安装于测量载体上；所述惯性导航系统、GPS接收机、星敏感器均与数据处理计算机连接；所述GPS天线与所述GPS接收机通信。

INS、GPS天线以及星敏感器固联安装于测量载体上，测量载体可以是测量船、测量车等运载工具。INS、GPS接收机和星敏感器的测量数据通过数据线传输到数据处理计算机中，在数据处理计算机中完成垂线偏差的解算。

所采用的星敏感器安装于测量载体时，应使其光轴尽可能朝向天顶。由于星敏感器的光轴指向测量精度优于横滚角测量精度，当光轴指向天顶时可以提高水平姿态角的测量精度，此外还可以抑制大气折射误差的影响。INS采用龙兴武等在2010年《中国惯性技术学报》第2期“激光陀螺单轴旋转惯性导航系统”论文中公开的单轴旋转式结构，通过单轴旋转结构对惯性器件进行周期性调制，可以有效抑制惯性器件的零偏误差，提高INS的姿态解算精度。

本发明还提供了一种利用上述装置动态测量垂线偏差的方法，该方法为：

1）启动INS/GPS姿态测量子系统，进行8小时以上对准。对准过程中INS/GPS姿态测量子系统即开始进行姿态测量，但无垂线偏差数据输出。

INS/GPS组合姿态测量的方法可参照全伟等著、国防工业出版社2011年出版的《惯性/天文/卫星组合导航技术》中120～123页介绍的方法建立INS/GPS组合导航的状态方程和观测方程，再根据第50～52页介绍的Kalman滤波算法，对INS的姿态误差、速度误差、位置误差、陀螺零偏、加速度计零偏进行估计，并利用估计得到的东向、北向和天向姿态误差

对INS的姿态矩阵输出进行修正，得到修正后的姿态信息为INS坐标系（b系）相对于计算导航坐标系n'的姿态矩阵

修正方法参考秦永元等著、西北工业大学出版社2012年出版的《卡尔曼滤波与组合导航原理》中337页中的公式（8.5.37）～（8.5.39）。Kalman滤波估计得到的X、Y、Z三个陀螺的零偏分别记为ε_x、ε_y、ε_z。在整个测量过程中INS/GPS组合姿态测量子系统连续输出姿态信息

和三个陀螺零偏的估计值ε_x、ε_y、ε_z。

2）完成INS/GPS姿态测量子系统的对准后，启动LGU/GPS测量子系统，开始进行垂线偏差测量，

3）令i=1；

4）利用星敏感器观测恒星，计算t_i,1时b系相对于n系的姿态矩阵

其中b系是指惯性导航坐标系O_b-x_by_bz_b；n系是指真实导航坐标系，即东-北-天坐标系O-xyz；

5）利用上述姿态矩阵对LGU/GPS姿态测量子系统进行姿态初始化，令LGU/GPS姿态测量子系统的姿态矩阵为

6）令j=2；

7）利用GPS输出的位置信息计算t_i,j-1到t_i,j时刻n系的变化矩阵

将t_i,j时刻测量载体所在的经度和纬度分别记为λ_i,j、L_i,j，则t_i,j-1和t_i,j时刻n系相对于地球坐标系e系的姿态矩阵分别为和

的计算公式如下：

$C_e^{n(i,j)}=\left[\begin{array}{ccc}-{\sin \mathrm{\λ}}_{i,j}& \cos {\mathrm{\λ}}_{i,j}& 0\\ -\sin L_{i,j}{\cos \mathrm{\λ}}_{i,j}& {-\sin L}_{i,j}{\sin \mathrm{\λ}}_{i,j}& {\cos L}_{i,j}\\ {\cos L}_{i,j}{\cos \mathrm{\λ}}_{i,j}& {\cos L}_{i,j}{\sin \mathrm{\λ}}_{i,j}& {\sin L}_{i,j}\end{array}\right];$

其中所述t_i,j时刻是指第i次利用星敏感器观测恒星后，第j次LGU/GPS姿态测量子系统姿态更新的时刻；t_i,j=t_i,j-1+δt，δt为激光陀螺的采样周期；

则 $C_{n(i,j-1)}^{n(i,j)}=C_e^{n(i,j)}\·{\left[C_e^{n(i,j-1)}\right]}^T,$ 其中，[]^T表示矩阵的转置；

8）将t_i,j-1到t_i,j时刻采样周期δt内三个激光陀螺的输出角增量分别记为并利用t_i,j时刻INS/GPS姿态测量子系统中Kalman滤波估计得到的三个激光陀螺的零偏ε_x、ε_y、ε_z修正三个激光陀螺的输出角增量，得到修正后的三个激光陀螺角增量分别为： ${\mathrm{\θ}}_x={\widehat{\mathrm{\θ}}}_x-{\mathrm{\ϵ}}_x,$ ${\mathrm{\θ}}_y={\widehat{\mathrm{\θ}}}_y-{\mathrm{\ϵ}}_y,$ ${\mathrm{\θ}}_z={\widehat{\mathrm{\θ}}}_z-{\mathrm{\ϵ}}_z;$

9）利用下式构造四元数q：

$q=\left[\begin{array}{c}\cos \frac{\mathrm{\θ}}{2}\\ \frac{{\mathrm{\θ}}_x}{\mathrm{\θ}}\sin \frac{\mathrm{\θ}}{2}\\ \frac{{\mathrm{\θ}}_y}{\mathrm{\θ}}\sin \frac{\mathrm{\θ}}{2}\\ \frac{{\mathrm{\θ}}_z}{\mathrm{\θ}}\sin \frac{\mathrm{\θ}}{2}\end{array}\right];$

其中 $\mathrm{\θ}=\sqrt{{\mathrm{\θ}}_x^2+{\mathrm{\θ}}_y^2+{\mathrm{\θ}}_z^2};$

10）计算t_i,j-1到t_i,j时刻测量INS坐标系的变化矩阵

$C_{b(i,j-1)}^{b(i,j)}=\left[\begin{array}{ccc}{q}_0^2+q_1^2-q_2^2-q_3^2& 2(q_1{q}_2-q_0{q}_3)& 2(q_1{q}_3+q_0{q}_2)\\ 2(q_1{q}_2+q_0{q}_3)& q_0^2-q_1^2+q_2^2-q_3^2& 2(q_2{q}_3-q_0{q}_1)\\ 2(q_1{q}_3-q_0{q}_2)& 2(q_2{q}_3+q_0{q}_1)& q_0^2-q_1^2-q_2^2+q_3^2\end{array}\right];$

其中，q₀、q₁、q₂、q₃分别为q的第1～4个元素；

11）计算t_i,j时刻b系相对于n系的姿态矩阵

$C_{b(i,j)}^{n(i,j)}=C_{n(i,j-1)}^{n(i,j)}\·C_{b(i,j-1)}^{n(i,j-1)}\·{\left[C_{b(i,j-1)}^{b(i,j)}\right]}^T;$

12）计算t_i,j时刻n'系到n系之间的姿态旋转矩阵

其中n'系是指O'-x'y'z'，即计算导航坐标系，

为t_i,j时刻INS/GPS姿态测量子系统输出的姿态矩阵；；

13）利用下式计算t_i,j时刻测量位置东向垂线偏差η_i,j和北向垂线偏差ξ_i,j： $\left\{\begin{array}{c}{\mathrm{\ξ}}_{i,j}=C_{n^{\′}(i,j)}^{n(i,j)}\left(\mathrm{3,1}\right)\\ {\mathrm{\η}}_{i,j}=C_{n^{\′}(i,j)}^{n(i,j)}\left(\mathrm{2,3}\right)\end{array}\right.;$ 其中，

表示

矩阵的第3行，第1列的元素；

表示

矩阵的第2行，第3列的元素； $C_{n^{\′}(i,j)}^{n(i,j)}\≈\left[\begin{array}{ccc}1& 0& -{\mathrm{\ξ}}_{i,j}\\ 0& 1& {\mathrm{\η}}_{i,j}\\ {\mathrm{\ξ}}_{i,j}& -{\mathrm{\η}}_{i,j}& 1\end{array}\right];$

14）判断是否接收到“结束测量”的命令，如果是，跳转到步骤16）；否则，执行步骤15）；

15）当j·δt>T_i时，令i=i+1，并跳转到步骤4）；否则，令j=j+1，并跳转到步骤7）；其中，T_i为星敏感器观测恒星时刻t_i,1到下一次星敏感器观测恒星时刻t_i+1,1所经历的时间，T_i的取值范围为6～8小时；

16）修正上述步骤13）中的东向垂线偏差η_i,j和北向垂线偏差ξ_i,j的低频误差和跳变误差；

17）结束测量。

与现有技术相比，本发明所具有的有益效果为：本发明构建了LGU/GPS姿态测量子系统进行姿态测量，将LGU/GPS的姿态输出与INS/GPS姿态测量子系统的姿态输出求差，该装置结构简单，可以直接测量n系与n′系坐标系偏差，不需要测量实施区域内垂线偏差的先验信息，因而本发明方法鲁棒性强，测量结果稳定；降低了GPS精度的需求，相比传统的矢量航空重力仪依赖于高精度差分GPS，本发明只需一般的GPS单点定位即可满足精度要求，因此拓展了测量方法的应用范围；惯性导航系统（INS）采用的是单轴旋转式的惯性导航系统，通过单轴旋转结构对惯性器件进行周期性调制，可以有效抑制惯性器件的零偏误差，提高了INS的导航解算精度；本发明可利用某测量点的垂线偏差数据替代星敏感器对LGU/GPS姿态测量子系统进行初始化，从而减小了对星敏感器等天文姿态基准的依赖。

附图说明

图1是本发明垂线偏差动态测量装置结构示意图；

图2是导航坐标系与垂线偏差的关系；

图3是IMU的结构示意图；

图4是本发明的垂线偏差动态测量方法流程图；

图5是利用LGU/GPS姿态测量子系统解算姿态矩阵

结合INS/GPS姿态测量子系统的姿态输出

计算东向和北向垂线偏差η、ξ的方法流程图。

具体实施方式

本发明采用的垂线偏差动态测量装置如图1所示，由INS1、GPS天线2、GPS接收机3、星敏感器4、测量载体5和数据处理计算机6构成。其中INS1、GPS天线2及星敏感器4固联安装于测量载体5上，测量载体可以是测量船、测量车等运载工具。INS1、GPS接收机3和星敏感器4的测量数据通过数据线传输到数据处理计算机6中，在数据处理计算机6中完成垂线偏差的解算。

所采用的星敏感器安装于测量载体时，应使其光轴尽可能朝向天顶。由于星敏感器的光轴指向测量精度优于横滚角测量精度，当光轴指向天顶时可以提高水平姿态角的测量精度，此外还可以抑制大气折射误差的影响。INS采用单轴旋转式结构，通过单轴旋转结构对惯性器件进行周期性调制，可以有效抑制惯性器件的零偏误差，提高INS的导航解算精度。

图4为本发明测量方法总体流程图，主要包含以下七个步骤：

步骤1：构建INS/GPS姿态测量子系统。

步骤2：构建LGU/GPS姿态测量子系统。

步骤3：启动INS/GPS姿态测量子系统，进行8小时以上对准。

步骤4：启动星敏感器和LGU/GPS姿态测量子系统，开始测量。

步骤5：利用LGU/GPS姿态测量子系统解算姿态矩阵

结合INS/GPS姿态测量子系统的姿态输出

计算东向和北向垂线偏差η、ξ。

步骤6：修正垂线偏差中的低频误差和跳变误差。

步骤7：结束测量。

如图1所示，定义INS的坐标系为b系O_b-x_by_bz_b，星敏感器坐标系为O_ss-x_ssy_ssz_ss。如图2所示，定义东-北-天坐标系O-xyz为真实导航坐标系，即n系，其中O-xy平面与地球的参考椭球面平行。由于垂线偏差的存在，惯性导航解算得到的计算导航坐标系与真实的导航坐标系存在偏差，定义计算导航坐标系为O'-x'y'z'，即n'系，计算导航坐标系的定义可参考高钟毓著、清华大学出版社2012年出版的《惯性导航系统技术》第225-226页的描述。

本发明具体实施步骤如下：

步骤1：构建INS/GPS姿态测量子系统。

由图1装置中的INS1、GPS天线2和GPS接收机3构建INS/GPS姿态测量子系统，INS/GPS组合姿态测量的方法可参照全伟等著、国防工业出版社2011年出版的《惯性/天文/卫星组合导航技术》中120～123页介绍的方法建立INS/GPS组合导航的状态方程和观测方程，再根据第50～52页介绍的Kalman滤波算法，对INS的姿态误差、速度误差、位置误差、陀螺零偏、加速度计零偏进行估计，并利用估计得到的东向、北向和天向姿态误差

对INS的姿态矩阵输出进行修正，得到修正后的姿态信息为INS坐标系（b系）相对于计算导航坐标系n'的姿态矩阵修正方法参考秦永元等著、西北工业大学出版社2012年出版的《卡尔曼滤波与组合导航原理》中337页中的公式（8.5.37）～（8.5.39）。Kalman滤波估计得到的X、Y、Z三个陀螺的零偏分别记为ε_x、ε_y、ε_z。在整个测量过程中INS/GPS组合姿态测量子系统连续输出姿态信息

和三个陀螺零偏的估计值ε_x、ε_y、ε_z。

步骤2：构建LGU/GPS姿态测量子系统。图1中INS1的核心测量部件为IMU（InertialMeasure Unit惯性测量单元），如图3所示，IMU中包含3个正交安装的激光陀螺111、112、113和3个正交安装的加速度计12。利用IMU中的3个激光陀螺构建新的测量组件，称为LGU，再利用LGU和GPS可构建另一个组合姿态测量子系统，称为LGU/GPS姿态测量子系统。LGU/GPS实现姿态测量算法将在后续步骤中详细说明。

步骤3：启动INS/GPS姿态测量子系统，进行8小时以上对准。对准过程中无垂线偏差数据输出。

步骤4：启动星敏感器和LGU/GPS姿态测量子系统，开始测量。

步骤5：利用LGU/GPS姿态测量子系统解算姿态矩阵

结合INS/GPS姿态测量子系统输出的姿态矩阵

计算东向和北向垂线偏差η、ξ，如图2所示。具体实施方法如下：

记第i次利用星敏感器观测恒星的时刻记为t_i,1。

步骤5.1：令i=1。

步骤5.2：利用星敏感器观测恒星，获取b系相对于真实导航坐标系n系的姿态矩阵

并对LGU/GPS姿态计算进行初始化，具体实施方法为：

步骤5.2.1：利用星敏感器观测恒星，获取星敏感器坐标系相对于n系的姿态矩阵

利用星敏感器与INS之间的安装矩阵

通过式（1）计算b系相对于n系的姿态矩阵

$C_b^n=C_{\mathrm{ss}}^n\·C_b^{\mathrm{ss}}---\left(1\right)$

其中“·”表示矩阵相乘，

采用王岩在2009年《战术导弹控制技术》期刊第3期“一种星敏感器与捷联惯导高精度安装误差标定方法”论文中公开的方法事先标定得到。

t_i,1时刻由式（1）计算得到的b系相对于n系的姿态矩阵记为

还可以采用以下替代方案计算

矩阵：

当已知t_i,1时刻测量载体所在位置上的垂线偏差η_i,1,ξ_i,1时，可利用η_i,1,ξ_i,1代替星敏感器数据计算b系相对于n系的姿态矩阵计算方法如下：

利用η_i,1,ξ_i,1由式（2）计算n系和n′系之间的旋转矩阵

$C_{n^{\′}(i,1)}^{n(i,1)}\≈\left[\begin{array}{ccc}1& 0& -{\mathrm{\ξ}}_{i,1}\\ 0& 1& {\mathrm{\η}}_{i,1}\\ {\mathrm{\ξ}}_{i,1}& -{\mathrm{\η}}_{i,1}& 1\end{array}\right]---\left(2\right)$

利用INS/GPS姿态测量子系统在t_i,1时刻输出的姿态矩阵

由式（3）计算

$C_{b(i,1)}^{n(i,1)}=C_{n^{\′}(i,1)}^{n(i,1)}\·C_{b(i,1)}^{n^{\′}(i,1)}---\left(3\right)$

使用该替代方案可以脱离对星敏感器等天文姿态基准的依赖。

步骤5.2.2：用矩阵对LGU/GPS姿态测量子系统进行姿态初始化，令LGU/GPS的姿态矩阵为

步骤5.3：以

为初值，利用激光陀螺和GPS的输出进行LGU/GPS姿态更新，计算b系相对于n系的姿态矩阵

更新算法如下：

记第i次利用星敏感器观测恒星后，第j次LGU/GPS姿态更新的时刻记为t_i,j，t_i,j=t_i,j-1+δt，δt为激光陀螺的采样周期。记t_i,j时刻b系相对于n系的姿态矩阵为

步骤5.3.1：令j=2。

步骤5.3.2：利用GPS输出的位置信息计算t_i,j-1到t_i,j时刻导航坐标系的变化矩阵。

t_i,j时刻载体所在的经度和纬度分别记为λ_i,j,L_i,j。则t_i,j-1和t_i,j时刻导航坐标系相对于地球坐标系e系的姿态矩阵分别记为

和

可由式(4)计算：

$C_e^{n(i,j)}=\left[\begin{array}{ccc}-{\sin \mathrm{\λ}}_{i,j}& \cos {\mathrm{\λ}}_{i,j}& 0\\ -\sin L_{i,j}{\cos \mathrm{\λ}}_{i,j}& {-\sin L}_{i,j}{\sin \mathrm{\λ}}_{i,j}& {\cos L}_{i,j}\\ {\cos L}_{i,j}{\cos \mathrm{\λ}}_{i,j}& {\cos L}_{i,j}{\sin \mathrm{\λ}}_{i,j}& {\sin L}_{i,j}\end{array}\right]---\left(4\right)$

将式（4）中j替换为j-1可计算得到

由式(5)计算

$C_{n(i,j-1)}^{n(i,j)}=C_e^{n(i,j)}\·{\left[C_e^{n(i,j-1)}\right]}^T---\left(5\right)$

步骤5.3.3：利用激光陀螺输出的角增量计算t_i,j-1到t_i,j时刻INS坐标系的变化矩阵 $C_{b(i,j-1)}^{b(i,j)}.$

t_i,j-1到t_i,j时刻采样周期δt内，三个激光陀螺111、112、113输出的角增量分别记为

采用步骤1由Kalman滤波估计得到的三个陀螺的零偏ε_x、ε_y、ε_z对陀螺角增量进行修正，修正后的陀螺角增量为： ${\mathrm{\θ}}_x={\widehat{\mathrm{\θ}}}_x-{\mathrm{\ϵ}}_x,$ ${\mathrm{\θ}}_y={\widehat{\mathrm{\θ}}}_y-{\mathrm{\ϵ}}_y,$ ${\mathrm{\θ}}_z={\widehat{\mathrm{\θ}}}_z-{\mathrm{\ϵ}}_z;$

由式（6）构造旋转四元数：

$q=\left[\begin{array}{c}\cos \frac{\mathrm{\θ}}{2}\\ \frac{{\mathrm{\θ}}_x}{\mathrm{\θ}}\sin \frac{\mathrm{\θ}}{2}\\ \frac{{\mathrm{\θ}}_y}{\mathrm{\θ}}\sin \frac{\mathrm{\θ}}{2}\\ \frac{{\mathrm{\θ}}_z}{\mathrm{\θ}}\sin \frac{\mathrm{\θ}}{2}\end{array}\right]---\left(6\right)$

其中

q₀、q₁、q₂、q₃分别为q的第1～4个元素。

可由式(7)计算：

$C_{b(i,j-1)}^{b(i,j)}=\left[\begin{array}{ccc}{q}_0^2+q_1^2-q_2^2-q_3^2& 2(q_1{q}_2-q_0{q}_3)& 2(q_1{q}_3+q_0{q}_2)\\ 2(q_1{q}_2+q_0{q}_3)& q_0^2-q_1^2+q_2^2-q_3^2& 2(q_2{q}_3-q_0{q}_1)\\ 2(q_1{q}_3-q_0{q}_2)& 2(q_2{q}_3+q_0{q}_1)& q_0^2-q_1^2-q_2^2+q_3^2\end{array}\right]---\left(7\right)$

步骤5.3.4：t_i,j时刻b系相对于n系的姿态矩阵可由式（8）计算：

$C_{b(i,j)}^{n(i,j)}=C_{n(i,j-1)}^{n(i,j)}\·C_{b(i,j-1)}^{n(i,j-1)}\·{\left[C_{b(i,j-1)}^{b(i,j)}\right]}^T---\left(8\right)$

其中[]^T表示矩阵的转置。

步骤5.4：计算t_i,j时刻n'系到n系之间的姿态旋转矩阵

t_i,j时刻INS/GPS姿态测量子系统输出的姿态矩阵记为由式（9）可计算矩阵 $C_{n^{\′}(i,j)}^{n(i,j)}.$

$C_{n^{\′}(i,j)}^{n(i,j)}=C_{b(i,j)}^{n(i,j)}\·{\left[C_{b(i,j)}^{n^{\′}(i,j)}\right]}^T---\left(9\right)$

步骤5.5：由

矩阵计算t_i,j时刻测量位置上的垂线偏差。

记东向垂线偏差为η_i,j，北向垂线偏差为ξ_i,j。η_i,j、ξ_i,j与

矩阵近似满足式（10）的关系：

$C_{n^{\′}(i,j)}^{n(i,j)}\≈\left[\begin{array}{ccc}1& 0& -{\mathrm{\ξ}}_{i,j}\\ 0& 1& {\mathrm{\η}}_{i,j}\\ {\mathrm{\ξ}}_{i,j}& -{\mathrm{\η}}_{i,j}& 1\end{array}\right]---\left(10\right)$

可以通过式（11）求得垂线偏差：

$\left\{\begin{array}{c}{\mathrm{\ξ}}_{i,j}=C_{n^{\′}(i,j)}^{n(i,j)}\left(\mathrm{3,1}\right)\\ {\mathrm{\η}}_{i,j}=C_{n^{\′}(i,j)}^{n(i,j)}\left(\mathrm{2,3}\right)\end{array}\right.---\left(11\right)$

其中

表示

矩阵的第3行，第1列的元素，

表示

矩阵的第2行，第3列的元素。

步骤5.6：判断是否收到“结束测量”的命令，如果接受到“结束测量”的命令则跳转到步骤6。否则，继续执行步骤5.7。

步骤5.7：本次星敏感器观测恒星时刻t_i,1到下一次星敏感器观测恒星时刻t_i+1,1所经历的时间记为T_i。当j·δt>T_i时，令i=i+1，并跳转到步骤5.2。否则，令j=j+1，并跳转到步骤5.3.2。

T_i可设置的范围为：6小时<T_i<8小时，可以根据测量的条件灵活设定。一般选择星敏感器观测到的恒星成像清晰，无云雾遮挡，测量载体运动状态平稳的时刻观测恒星。

步骤6：修正垂线偏差的低频误差和跳变误差，其具体修正方法如下：

步骤6.1：修正垂线偏差数据的跳变误差。

由步骤5得到的垂线偏差测量序列为η_i,j、ξ_i,j（i=1…M，j=1…N_i），M为全部测量过程利用星敏感器对LGU/GPS姿态初始化的次数，N_i为第i次观测恒星后LGU/GPS连续进行步骤5.3姿态更新的次数。

由于星敏感器具有一定的姿态测量误差，因此每次在步骤5.2对LGU/GPS进行姿态初始化时，t_i,1时刻的垂线偏差η_i,1、ξ_i,1相对于其前一时刻

的垂线偏差存在跳变误差。记第i次对LGU/GPS进行姿态初始化时产生的东向垂线偏差和北向垂线偏差的跳变误差分别为

和

的计算方法如下：

步骤6.1.1：令 ${\&dtri;}_{\mathrm{\&eta;},1}=0,{\&dtri;}_{\mathrm{\&xi;},1}=0$

步骤6.1.2：令i=i+1

步骤6.1.3： $\left\{\begin{array}{c}{\&dtri;}_{\mathrm{\&eta;},i}={\&dtri;}_{\mathrm{\&eta;},i-1}+({\mathrm{\&eta;}}_{i,1}-{\mathrm{\&eta;}}_{i-1,N_{i-1}})\\ {\&dtri;}_{\mathrm{\&xi;},i}={\&dtri;}_{\mathrm{\&xi;},i-1}+({\mathrm{\&xi;}}_{i,1}-{\mathrm{\&xi;}}_{i-1,N_{i-1}})\end{array}\right.---\left(12\right)$ 步骤6.1.4： $\left\{\begin{array}{c}{\mathrm{\&eta;}}_{i,j}={\mathrm{\&eta;}}_{i,j}-{\&dtri;}_{\mathrm{\&eta;},i}\\ {\mathrm{\&xi;}}_{i,j}={\mathrm{\&xi;}}_{i,j}-{\&dtri;}_{\mathrm{\&xi;},i}\end{array}\right.---\left(13\right)$

其中对全部1≤j≤N_i执行式（13）。

步骤6.1.5：若i<M，，跳转到步骤6.1.2，直到i=M。

步骤6.2：修正垂线偏差数据的低频误差。

为了便于描述，将全部测量时刻t_i,j按时间顺序排列，并记为单下标的形式t_k（1≤k≤N），k为第k个垂线偏差测量数据点的序号，N为垂线偏差测量数据点的总个数。将t_i,j时刻对应的垂线偏差数据η_i,j、ξ_i,j记为η_k、ξ_k，相应的EGM2008垂线偏差数据

记为 ${\widehat{\mathrm{\&eta;}}}_k.{\widehat{\mathrm{\&xi;}}}_k..$

步骤6.2.1：根据t_k时刻载体所在的位置读取全球重力数据库EGM2008的垂线偏差数据，记为

主要包含垂线偏差数据的低频信息。

步骤6.2.2：将测量得到的垂线偏差数据η_k、ξ_k中的低频部分利用式（14）移除，得到剩余量δη_k、δξ_k包含垂线偏差的高频信息。

$\left\{\begin{array}{c}{\mathrm{\&delta;\&eta;}}_k={\mathrm{\&eta;}}_k-{\widehat{\mathrm{\&eta;}}}_k\\ {\mathrm{\&delta;\&xi;}}_k={\mathrm{\&xi;}}_k-{\widehat{\mathrm{\&xi;}}}_k\end{array}\right.---\left(14\right)$

对全部1≤k≤N执行式（14）的操作。

步骤6.2.3：利用多项式模型拟合δη_k、δξ_k中存在的缓变误差，并将该缓变误差扣除，得到垂线偏差的高频信息

具体的修正方法如下：

以k（1≤k≤N）为自变量，对数据序列δη_k、δξ_k分别进行多项式拟合，采用的拟合方法为最小二乘法，拟合的模型为：

$\left\{\begin{array}{c}{y}_{\mathrm{\&eta;}}\left(k\right)=a_{\mathrm{\&eta;},0}+a_{\mathrm{\&eta;},1}\&CenterDot;k+a_{\mathrm{\&eta;},2}\&CenterDot;k^2\\ y_{\mathrm{\&xi;}}\left(k\right)=a_{\mathrm{\&xi;},0}+a_{\mathrm{\&xi;},1}\&CenterDot;k+a_{\mathrm{\&xi;},2}\&CenterDot;k^2\end{array}\right.---\left(15\right)$

其中a_η,0、a_η,1、a_η,2为由自变量k和应变量δη_k拟合得到的模型系数，y_η(k)为由拟合模型计算得到的第k个垂线偏差测量数据点的δη_k缓变误差。a_ξ,0、a_ξ,1、a_ξ,2为由自变量k和应变量δξ_k拟合得到的模型系数，y_ξ(k)为由拟合模型计算得到的第k个数据点的δξ_k缓变误差。

对全部1≤k≤N执行式（15）的操作。

利用式（16）扣除δη_k和δξ_k中的缓变误差：

$\left\{\begin{array}{c}\mathrm{\&delta;}{\widehat{\mathrm{\&eta;}}}_k={\mathrm{\&delta;\&eta;}}_k-y_{\mathrm{\&eta;}}\left(k\right)\\ \mathrm{\&delta;}{\widehat{\mathrm{\&xi;}}}_k={\mathrm{\&delta;\&xi;}}_k-y_{\mathrm{\&xi;}}\left(k\right)\end{array}\right.---\left(16\right)$

对全部1≤k≤N执行式（16）的操作。

步骤6.2.4：将EGM2008的垂线偏差数据序列

与式（17）求得的垂线偏差高频量δη_k、δξ_k相加即可得到垂线偏差。

$\left\{\begin{array}{c}{\mathrm{\&eta;}}_k={\widehat{\mathrm{\&eta;}}}_k+\mathrm{\&delta;}{\widehat{\mathrm{\&eta;}}}_k\\ {\mathrm{\&xi;}}_k={\widehat{\mathrm{\&xi;}}}_k+\mathrm{\&delta;}{\widehat{\mathrm{\&xi;}}}_k\end{array}\right.---\left(17\right)$

对全部1≤k≤N执行式（17）的操作。

步骤7：结束测量。

Claims (5)

1.一种基于天文姿态基准保持的垂线偏差动态测量装置，包括由惯性导航系统（1）、GPS天线（2）、GPS接收机（3）组成的INS/GPS姿态测量子系统、星敏感器（4），所述惯性导航系统（1）包括三个正交安装的激光陀螺（11），所述三个正交安装的激光陀螺（11）与GPS接收机（3）组成LGU/GPS姿态测量子系统，三个激光陀螺（11）均与GPS接收机（3）通信；其特征在于，所述惯性导航系统（1）、GPS天线（2）、GPS接收机（3）、星敏感器（4）均安装于测量载体（5）上；所述惯性导航系统（1）、GPS接收机（3）、星敏感器（4）均与数据处理计算机（6）连接；所述GPS天线（2）与所述GPS接收机（3）通信。

2.一种利用权利要求1所述装置动态测量垂线偏差的方法，其特征在于，该方法为：

1）启动INS/GPS姿态测量子系统，进行8小时以上对准；

2）启动LGU/GPS测量子系统；

3）令i=1；

4）利用星敏感器观测恒星，计算t_i,1时b系相对于n系的姿态矩阵

其中b系是指惯性导航系统坐标系O_b-x_by_bz_b；n系是指真实导航坐标系，即东-北-天坐标系O-xyz；

5）利用上述姿态矩阵

对LGU/GPS姿态测量子系统进行姿态初始化，令LGU/GPS姿态测量子系统的姿态矩阵为

6）令j=2；

7）利用GPS输出的位置信息计算t_i,j-1到t_i,j时刻n系的变化矩阵 $C_{n(i,j-1)}^{n(i,j)}=C_e^{n(i,j)}\&CenterDot;{\left[C_e^{n(i,j-1)}\right]}^T,$ 其中，[]^T表示矩阵的转置；

和

分别为t_i,j-1和t_i,j时刻n系相对于地球坐标系e系的姿态矩阵，

的计算公式如下：

$C_e^{n(i,j)}=\left[\begin{array}{ccc}-{\sin \mathrm{\&lambda;}}_{i,j}& \cos {\mathrm{\&lambda;}}_{i,j}& 0\\ -\sin L_{i,j}{\cos \mathrm{\&lambda;}}_{i,j}& {-\sin L}_{i,j}{\sin \mathrm{\&lambda;}}_{i,j}& {\cos L}_{i,j}\\ {\cos L}_{i,j}{\cos \mathrm{\&lambda;}}_{i,j}& {\cos L}_{i,j}{\sin \mathrm{\&lambda;}}_{i,j}& {\sin L}_{i,j}\end{array}\right];$

其中t_i,j时刻是指第i次利用星敏感器观测恒星后，第j次LGU/GPS姿态测量子系统姿态更新的时刻；t_i,j=t_i,j-1+δt，δt为激光陀螺的采样周期；λ_i,j、L_i,j分别为t_i,j时刻测量载体所在的经度和纬度；

8）将t_i,j-1到t_i,j时刻采样周期δt内三个激光陀螺的输出角增量分别记为

并利用三个激光陀螺的零偏ε_x、ε_y、ε_z修正三个激光陀螺的输出角增量，得到修正后的三个激光陀螺角增量分别为： ${\mathrm{\&theta;}}_x={\widehat{\mathrm{\&theta;}}}_x-{\mathrm{\&epsiv;}}_x,$ ${\mathrm{\&theta;}}_y={\widehat{\mathrm{\&theta;}}}_y-{\mathrm{\&epsiv;}}_y,$ ${\mathrm{\&theta;}}_z={\widehat{\mathrm{\&theta;}}}_z-{\mathrm{\&epsiv;}}_z;$

9）利用下式构造四元数q：

$q=\left[\begin{array}{c}\cos \frac{\mathrm{\&theta;}}{2}\\ \frac{{\mathrm{\&theta;}}_x}{\mathrm{\&theta;}}\sin \frac{\mathrm{\&theta;}}{2}\\ \frac{{\mathrm{\&theta;}}_y}{\mathrm{\&theta;}}\sin \frac{\mathrm{\&theta;}}{2}\\ \frac{{\mathrm{\&theta;}}_z}{\mathrm{\&theta;}}\sin \frac{\mathrm{\&theta;}}{2}\end{array}\right];$

其中 $\mathrm{\&theta;}=\sqrt{{\mathrm{\&theta;}}_x^2+{\mathrm{\&theta;}}_y^2+{\mathrm{\&theta;}}_z^2};$

10）计算t_i,j-1到t_i,j时刻测量惯性导航系统坐标系的变化矩阵

$C_{b(i,j-1)}^{b(i,j)}=\left[\begin{array}{ccc}{q}_0^2+q_1^2-q_2^2-q_3^2& 2(q_1{q}_2-q_0{q}_3)& 2(q_1{q}_3+q_0{q}_2)\\ 2(q_1{q}_2+q_0{q}_3)& q_0^2-q_1^2+q_2^2-q_3^2& 2(q_2{q}_3-q_0{q}_1)\\ 2(q_1{q}_3-q_0{q}_2)& 2(q_2{q}_3+q_0{q}_1)& q_0^2-q_1^2-q_2^2+q_3^2\end{array}\right];$

其中，q₀、q₁、q₂、q₃分别为q的第1～4个元素；

11）计算t_i,j时刻b系相对于n系的姿态矩阵

$C_{b(i,j)}^{n(i,j)}=C_{n(i,j-1)}^{n(i,j)}\&CenterDot;C_{b(i,j-1)}^{n(i,j-1)}\&CenterDot;{\left[C_{b(i,j-1)}^{b(i,j)}\right]}^T;$

12）计算t_i,j时刻n'系到n系之间的姿态旋转矩阵

其中n'系是指O'-x'y'z'，即计算导航坐标系，

为t_i,j时刻INS/GPS姿态测量子系统输出的姿态矩阵；

13）利用下式计算t_i,j时刻测量位置东向垂线偏差η_i,j和北向垂线偏差 ${\mathrm{\&xi;}}_{i,j}:\left\{\begin{array}{c}{\mathrm{\&xi;}}_{i,j}=C_{n^{\&prime;}(i,j)}^{n(i,j)}\left(\mathrm{3,1}\right)\\ {\mathrm{\&eta;}}_{i,j}=C_{n^{\&prime;}(i,j)}^{n(i,j)}\left(\mathrm{2,3}\right)\end{array}\right.;$ 其中，

表示

矩阵的第3行，第1列的元素；

表示矩阵的第2行，第3列的元素； $C_{n^{\&prime;}(i,j)}^{n(i,j)}\&ap;\left[\begin{array}{ccc}1& 0& -{\mathrm{\&xi;}}_{i,j}\\ 0& 1& {\mathrm{\&eta;}}_{i,j}\\ {\mathrm{\&xi;}}_{i,j}& -{\mathrm{\&eta;}}_{i,j}& 1\end{array}\right];$

14）判断是否接收到“结束测量”的命令，如果是，跳转到步骤16）；否则，执行步骤15）；

15）当j·δt>T_i时，令i=i+1，并跳转到步骤4）；否则，令j=j+1，并跳转到步骤7）；其中，T_i为星敏感器观测恒星时刻t_i,1到下一次星敏感器观测恒星时刻t_i+1,1所经历的时间，T_i的取值范围为6～8小时；

16）修正上述步骤15）中的东向垂线偏差η_i,j和北向垂线偏差ξ_i,j的低频误差和跳变误差，得到修正后的东向垂线偏差和北向垂线偏差；

17）结束测量。

3.根据权利要求2所述的方法，其特征在于，所述步骤2）中，的计算公式如下：

其中，

为星敏感器与惯性导航系统之间的安装矩阵；

为星敏感器坐标系相对于n系的姿态矩阵。

4.根据权利要求2或3所述的方法，其特征在于，所述步骤14）中，东向垂线偏差η_i,j和北向垂线偏差ξ_i,j的跳变误差修正方法如下：

1）初始化东向垂线偏差的跳变误差令

2）令i=i+1；

3）利用下式计算东向垂线偏差的跳变误差

${\&dtri;}_{\mathrm{\&xi;},i}:\left\{\begin{array}{c}{\&dtri;}_{\mathrm{\&eta;},i}={\&dtri;}_{\mathrm{\&eta;},i-1}+({\mathrm{\&eta;}}_{i,1}-{\mathrm{\&eta;}}_{i-1,N_{i-1}})\\ {\&dtri;}_{\mathrm{\&xi;},i}={\&dtri;}_{\mathrm{\&xi;},i-1}+({\mathrm{\&xi;}}_{i,1}-{\mathrm{\&xi;}}_{i-1,N_{i-1}})\end{array}\right.;$ 其中N_i为第i次观测恒星后LGU/GPS姿态测量子系统连续进行更新的次数；

4）利用下式修正东向垂线偏差η_i,j和北向垂线偏差ξ_i,j，得到修正后的东向垂线偏差η_i,j'和北向垂线偏差 ${{\mathrm{\&xi;}}_{i,j}}^{\&prime;}:\left\{\begin{array}{c}{{\mathrm{\&eta;}}_{i,j}}^{\&prime;}={\mathrm{\&eta;}}_{i,j}-{\&dtri;}_{\mathrm{\&eta;},i}\\ {{\mathrm{\&xi;}}_{i,j}}^{\&prime;}={\mathrm{\&xi;}}_{i,j}-{\&dtri;}_{\mathrm{\&xi;},i}\end{array}\right.;$

5）重复上述步骤2）～步骤4），直到i=M；其中M为利用姿态矩阵

对LGU/GPS姿态测量子系统进行姿态初始化的次数。

5.根据权利要求2或3所述的方法，其特征在于，东向垂线偏差η_i,j和北向垂线偏差ξ_i,j的低频误差修正方法如下：

1）根据t_i,j时刻载体所在的位置读取全球重力数据库EGM2008的东向和北向垂线偏差数据，记为

2）利用下式移除东向垂线偏差η_i,j和北向垂线偏差ξ_i,j的低频部分，得到包含东向垂线偏差η_i,j和北向垂线偏差ξ_i,j高频信息的东向垂线偏差余量和北向垂线偏差余量δη_i,j、δξ_i,j：

$\left\{\begin{array}{c}{\mathrm{\&delta;\&eta;}}_{i.j}={\mathrm{\&eta;}}_{i.j}-{\widehat{\mathrm{\&eta;}}}_{i.j}\\ {\mathrm{\&delta;\&xi;}}_{i.j}={\mathrm{\&xi;}}_{i.j}-{\widehat{\mathrm{\&xi;}}}_{i.j}\end{array}\right.;$

3）利用多项式模型拟合δη_i,j、δξ_i,j中存在的缓变误差，并去除该缓变误差，得到东向垂线偏差和北向垂线偏差的高频信息所述多项式模型为：

$\left\{\begin{array}{c}{y}_{\mathrm{\&eta;}}\left(k\right)=a_{\mathrm{\&eta;},0}+a_{\mathrm{\&eta;},1}\&CenterDot;k+a_{\mathrm{\&eta;},2}\&CenterDot;k^2\\ y_{\mathrm{\&xi;}}\left(k\right)=a_{\mathrm{\&xi;},0}+a_{\mathrm{\&xi;},1}\&CenterDot;k+a_{\mathrm{\&xi;},2}\&CenterDot;k^2\end{array}\right.;$

其中a_η,0、a_η,1、a_η,2为由自变量k和应变量

拟合得到的模型系数，y_η(k)为由多项式模型计算得到的第k个垂线偏差测量数据点的的缓变误差；a_ξ,0、a_ξ,1、a_ξ,2为由自变量k和应变量拟合得到的模型系数，y_ξ(k)为由多项式模型计算得到的第k个垂线偏差测量数据点的

的缓变误差；1≤k≤N，N为垂线偏差测量测量数据点的总个数；

4）将

分别与

相加，得到修正后的东向垂线偏差η_i,j'和北向垂线偏差ξ_i,j'：

$\left\{\begin{array}{c}{\mathrm{\&eta;}}_{{i,j}^{\&prime;}}={\widehat{\mathrm{\&eta;}}}_{i,j}+\mathrm{\&delta;}{\widehat{\mathrm{\&eta;}}}_{i,j}\\ {\mathrm{\&xi;}}_{{i,j}^{\&prime;}}={\widehat{\mathrm{\&xi;}}}_{i,j}+\mathrm{\&delta;}{\widehat{\mathrm{\&xi;}}}_{i,j}\end{array}\right..$

Images