CN103364841B - 一种用于航空重力测量中星座跳变误差的平滑消除方法 - Google Patents

Abstract

本发明公开了一种用于航空重力测量中星座跳变误差的平滑消除方法，首先对每条测线出现导航星座跳变的时刻进行分类和标记，并按照跳变的类型进行误差估计和补偿；然后，对跳变前后历元的计算位置的位置进行加权平滑补偿，以消除星座跳变对载***置解算的影响，以获取平滑的位置序列。本发明具有原理简单、操作简便、可完全消除跳变误差、有利于提高解算精度等优点。

Description

一种用于航空重力测量中星座跳变误差的平滑消除方法

技术领域

本发明主要涉及到航空重力测量技术领域，特指一种用于航空重力测量中星座跳变误差的平滑消除方法。

背景技术

航空重力测量是目前获取近地信息的重要手段之一，确定航空载体自身的加速度等状态参数是实现航空重力测量数据处理的前提和关键。目前，普遍利用卫星导航***的高精度定位能力进行实载体加速度的精密测量，因此测量结果受到各项GNSS误差源的影响，其中星座跳变误差是对结果影响较大的重要误差源。

航空重力测线长度可达数百公里，由于仰角约束、信号遮挡等原因，导航卫星的可见性发生改变是不可避免的现象，进而导致不同历元解算加速度时将采用不同的空间卫星星座。前后历元解算的空间基准发生改变时，解算的载***置会引起跳变，该跳变误差在求取加速度的差分过程中将被逐次放大，对提取的重力信号产生严重影响。

跳变误差的量级与站的可见卫星总数、卫星仰角、卫星轨道类型、星座几何构型等因素有关。研究发现，当载波相位观测值的测量误差为厘米级时，星座构型跳变虽然对定位结果的影响为毫米级，可忽略不计，但对差分运算得到的载体加速度的影响可达数百mGal（1mGal=10^-5m/s²）甚至上千mGal的影响。所提取的重力异常信号的量级一般在100mGal量级，因此该项误差对微弱重力信号的提取将会产生严重影响，低通滤波时会产生边缘效应将影响测量的多个重力异常值，目前在数据处理中一般采取以下几种手段来减弱该项误差的影响：

（1）卫星剔除法。在测线数据处理前，首先剔除该测线上发生可见性变化的卫星。该方法存在的弊端是可利用的卫星资源减少，解算结果的精度降低；数据处理程序复杂；仅适用于较短测线。

（2）协方差加权法。通过设计加权矩阵，低仰角卫星将被赋予较小的权重以降低对星座跳变误差的影响。该方法存在的弊端是加权矩阵的设计缺乏理论性，只能在一定程度上减小跳变误差的影响而无法完全消除，并且对非低仰角原因导致的星座跳变误差修正没有效果，例如接收机由于被遮挡等原因导致的对高仰角卫星的失锁等。

发明内容

本发明要解决的技术问题就在于：针对现有技术存在的技术问题，本发明提供一种原理简单、操作简便、可完全消除跳变误差、有利于提高解算精度的用于航空重力测量中星座跳变误差的平滑消除方法。

为解决上述技术问题，本发明采用以下技术方案：

一种用于航空重力测量中星座跳变误差的平滑消除方法，首先对每条测线出现导航星座跳变的时刻进行分类和标记，并按照跳变的类型进行误差估计和补偿；然后，对跳变前后历元的计算位置的位置进行加权平滑补偿，以消除星座跳变对载***置解算的影响，以获取平滑的位置序列。

作为本发明的进一步改进：

所述加权平滑补偿的过程为：

（1）某计算时刻t，星座跳变前的位置x^a可通过最小二乘算法求解，并对跳变后的解x^b进行估计：

x^b=x^a+△x

（2）在解x^a和x^b之间进行加权平均，即：

$\stackrel{\‾}{x}\left(t\right)=\mathrm{\ω}\left(t\right)x^b+[1-\mathrm{\ω}\left(t\right)]x^a$

上式中，0≤ω(t)≤1称为平滑因子，且满足：T为位置序列跳变时刻，T_L为数据平滑的位置序列长度，称为平滑窗口长度；

（3）由以上步骤（1）和步骤（2）得到加权平滑方案：

$\stackrel{\&OverBar;}{x}\left(t\right)=\left\{\begin{array}{cc}{x}^a& t<T\\ \mathrm{\&omega;}\left(t\right)x_t^b+[1-\mathrm{\&omega;}\left(t\right)]x_t^a& T\&le;t\&le;T+T_L\\ x^b& t>T+T_L\end{array}\right..$

作为本发明的进一步改进：

如果当前历元存在星座跳变，则置位跳变标记，判断星座跳变误差的类型：

如果为I类跳变，按照下式计算跳变的大小：△x_p=x₁-x₂；

如果为II类跳变，按照下式计算跳变的大小：△x^p=σ₀·△GDOP^p，△GDOP即为星座构型改变前后GDOP值的改变量，σ_x=σ₀·GDOP；然后，按照步骤（3）所示的加权平均方案对位置解算结果进行平滑处理。

作为本发明的进一步改进：

所述I类跳变为导航卫星增加；由于仰角变高或者卫星重新锁定原因，导航卫星出现在测线的后半段，导航卫星可见后导致定位精度正向跳变；

所述II类跳变为导航卫星减少；由于仰角变低或者卫星被遮挡原因，导航卫星出现在测线的前半段，导航卫星不可见后导致定位精度负向跳变。

与现有技术相比，本发明为基于卫星导航星座的优化技术，首先对星座跳变引起的位置跳变进行估计，进而通过平滑滤波器对跳变误差进行补偿。相比于传统的跳变补偿方法，该方法具有以下优点：

1、本发明充分利用了空间卫星资源。导航卫星作为载体加速度解算的空间基准，应尽可能多的采用，有利于提高解算精度，本发明由于不需剔除可见性发生变化的卫星，因此导航接收机可见的卫星均可参与计算，有利于获得载体加速度的最优解。

2、本发明处理流程简单，自动化程序容易实现。在数据处理中利用基于本方法设计的自动化算法，不需人为干预，自动化程度高。

3、本发明应用后跳变误差可完全消除。

附图说明

图1是位置序列跳变误差对差分计算的载体加速度影响的示意图。

图2是本发明用于星座跳变误差的平滑消除方法在具体应用实例中的流程示意图。

具体实施方式

以下将结合说明书附图和具体实施例对本发明做进一步详细说明。

航空重力测量中的载体加速度的求解一般事后处理。首先求解出载体的位置序列，并对位置序列进行数值差分以获取载体的速度、加速度值。因此，位置序列中的跳变会在差分过程中引入额外的误差。另外，所求取的载体加速度由于包含大量的随机性测量噪声，这就需要通过低通滤波器抑制高频噪声，提取重力信号。上述跳变误差也会在低通滤波中引起边缘效应，进而对所求取的加速度序列产生影响。

在进行载体加速度解算时，首先，通过采用高精度载波差分GNSS***（DGNSS）确定载体的位置序列，再通过两次数据差分，获得载体的动态速度和加速度序列。

设GNSS基准站为m，航空测量飞机机载接收机为k，则双差观测方程可写为：

$\mathrm{\&lambda;}(\mathrm{\&Delta;}{\&dtri;\mathrm{\&Phi;}}_{m,k}^{p.q}+\mathrm{\&Delta;}\&dtri;N_{m,k}^{p.q})=\mathrm{\&Delta;}\&dtri;R_{m,k}^{p.q}+\mathrm{\&epsiv;}---\left(1\right)$

式（1）中，上标p,q表示卫星，下标m,k表示卫星导航接收机，“△▽”为双差运算符，并且定义λ为载波波长；Φ为载波相位观测值；N为在载波相位模糊度；R为星地几何距离；ε为双差测量误差。

线性化方程（1），可得双差解算的法方程为：

A·x=y（2）

式（2）中：x为待求的载***置矢量，A为观测矩阵，反映导航卫星星座的几何构型，y为线性化之后自由项，并且有：

$A=\left[\begin{array}{c}{\left(h_k^{p,1}\right)}^T\\ {\left(h_k^{p,2}\right)}^T\\ .\\ .\\ .\\ {\left(h_k^{p,N}\right)}^T\end{array}\right]---\left(3\right)$

式（3）中，上标表示导航卫星；N为基准站和机载接收机共视的卫星总数；p为选定的参考卫星；称为接收机k对卫星q的观测矢量，并且：

$h_k^{p,q}=\left[\begin{array}{c}{e}_{\mathrm{kx}}^q-e_{\mathrm{kx}}^p\\ e_{\mathrm{kq}}^q-e_{\mathrm{ky}}^p\\ e_{\mathrm{kz}}^q-e_{\mathrm{kz}}^p\end{array}\right]---\left(4\right)$

式（4）中，e为站星观测视线方向的单位矢量，并且有：

$e_{\mathrm{kx}}^i=\frac{x^i-x_k}{R_k^i},e_{\mathrm{ky}}^i=\frac{y^i-y_k}{R_k^i},e_{\mathrm{kz}}^i=\frac{z^i-z_k}{R_k^i}---\left(5\right)$

式（5）中，xⁱ、yⁱ、zⁱ分别为导航卫星i在地心地固坐标系中的位置分量，可由导航星历进行计算。根据最小二乘估计原理，载体加速度可解算为：

$\hat{x}={\left(A^{T}A\right)}^{-1}{A}^{T}y---\left(6\right)$

载***置解算结果的精度可估计为：

${\mathrm{\&sigma;}}_x^2={\mathrm{\&sigma;}}_0^2\&CenterDot;{\left(A^{T}A\right)}^{-1}---\left(7\right)$

式（7）中，为相位观测值的方差，则双差定位精度衰减因子GDOP为：

$\mathrm{GDOP}=\sqrt{\mathrm{trace}\left[{\left(A^{T}A\right)}^{-1}\right]}---\left(8\right)$

式（8）中，trace(·)为矩阵求迹运算。

解算出载体的位置序列后，可采用一阶中心差分滤波器求解载体的速度、加速度，即：

$\left\{\begin{array}{c}\stackrel{\&CenterDot;}{x}=\frac{x(t+\mathrm{\&Delta;t})-x(t-\mathrm{\&Delta;t})}{2\mathrm{\&Delta;t}}\\ \stackrel{\&CenterDot;\&CenterDot;}{x}=\frac{\stackrel{\&CenterDot;}{x}(t+\mathrm{\&Delta;t})-\stackrel{\&CenterDot;}{x}(t-\mathrm{\&Delta;t})}{2\mathrm{\&Delta;t}}\end{array}\right.---\left(9\right)$

式（9）中，△t为解算的位置序列的时间间隔。

根据以上分析，基于GNSS确定的载***置发生跳变的根本原因是：星座构型（以GDOP值的形式体现）发生跳变，以至对差分得到的载体速度和加速度序列产生影响。

如图1所示，用于加速度确定的导航星座构型在时刻T发生了跳变，进而引起确定的载***置序列发生跳变，该跳变在微分求取载体速度、加速度的过程中将被逐级放大而产生重大影响，例如当△x=5cm时，加速度的跳变影响将达到5000mGal。

本发明将星座跳变分成以下三种情况：

I类跳变：导航卫星增加。由于仰角变高或者卫星重新锁定等原因，导航卫星出现在测线的后半段，导航卫星可见后导致定位精度正向跳变；

II类跳变：导航卫星减少。由于仰角变低或者卫星被遮挡等原因，导航卫星出现在测线的前半段，导航卫星不可见后导致定位精度负向跳变；

III类跳变：导航卫星短时段内连续发生可见性变化。由于卫星短时段被遮挡等原因，在测线内的某个时段内，导航卫星连续发生可见性变化，每次变化时均会导致导航精度发生跳变。

通过式（2）可知，当采用不同的星座进行载体加速度解算时，会导致解算结果存在差异。假定在某一时刻，分别选取两个不同的导航星座，分别记为S1、S2，导航卫星数目分别为N1、N2，其中前N颗可见星一致，则观测矩阵可分别表示为：

$\left\{\begin{array}{c}{A}_{S1}={(h_0,h_1,...,h_N,...{h_N}_1)}^T\stackrel{\mathrm{\&Delta;}}{=}\left(A\right|A_1)\\ A_{S2}={(h_0,h_1,...,h_N,...{h_N}_2)}^T\stackrel{\mathrm{\&Delta;}}{=}\left(A\right|A_2)\end{array}\right.---\left(10\right)$

星座S1、S2对应的自由项可分别记为：

$\left\{\begin{array}{c}{L}_{S1}={\left[L\right|L_1]}^T\\ L_{S2}={\left[L\right|L_2]}^T\end{array}\right.---\left(11\right)$

根据式（6）可得到采用S1、S2的解算结果：

$\left\{\begin{array}{c}{x}_1={\left({A_{S1}}^T{A}_{S1}\right)}^{-1}{A_{S1}}^T{L}_{S1}={(A^{T}A+{A_1}^T{A}_1)}^{-1}(A^{T}L+{A_1}^T{L}_1)\\ x_2={\left({A_{S2}}^T{A}_{S2}\right)}^{-1}{A_{S2}}^T{L}_{S2}={(A^{T}A+{A_2}^T{A}_2)}^{-1}(A^{T}L+{A_2}^T{L}_2)\end{array}\right.---\left(12\right)$

根据矩阵和求逆公式，可以得到：

$\left\{\begin{array}{c}{(A^{T}A+{A_1}^T{A}_1)}^{-1}={\left(A^{T}A\right)}^{-1}-{\left(A^{T}A\right)}^{-1}\left({A_1}^T{A}_1\right){\left(A^{T}A\right)}^{-1}{(I+{A_1}^T{A}_1{\left(A^{T}A\right)}^{-1})}^{-1}\\ {(A^{T}A+{A_2}^T{A}_2)}^{-1}={\left(A^{T}A\right)}^{-1}-{\left(A^{T}A\right)}^{-1}\left({A_2}^T{A}_2\right){\left(A^{T}A\right)}^{-1}{(I+{A_2}^T{A}_2{\left(A^{T}A\right)}^{-1})}^{-1}\end{array}\right.---\left(13\right)$

由此得到解算结果差异为：

$x_1-x_2={\left(A^{T}A\right)}^{-1}(A_1^T{L}_1-A_2^T{L}_2)-[{\left(A^{T}A\right)}^{-1}{A}_1^T{A}_1{x}_1-{\left(A^{T}A\right)}^{-1}{A}_2^T{A}_2{x}_2]---\left(14\right)$

化简，立得：

$x_1={(I+{\left(A^{T}A\right)}^{-1}{A_1}^T{A}_1)}^{-1}(I+{\left(A^{T}A\right)}^{-1}{A_2}^T{A}_2)x_2+{(A^{T}A+{A_1}^T{A}_1)}^{-1}({A_1}^T{L}_1-{A_2}^T{L}_2)---\left(15\right)$

式（15）反应了不同导航星座S1、S2的解算结果x₁,x₂之间的关系。在解算出x₂后，可计算出基于星座S1的解x₁，进而可计算出星座由S1改变为S2后解算结果的跳变值x₁-x₂。

对于I类跳变误差的估计：

当单颗卫星的可见性发生变化。假设S1比S2多一颗导航星p，即A_S1=[A_S2|h_p]，令 $G_{s1}={\left(A_{s1}^T{A}_{s1}\right)}^{-1};G_{s2}={\left(A_{s2}^T{A}_{s2}\right)}^{-1},B_{s1}=A_{s1}^T{L}_{s1};B_{s2}=A_{s2}^T{L}_{s2},$ L_s1=[(L_s2)^Tσ_p]T，有：

$\left\{\begin{array}{c}{x}_1=G_{S1}{B}_{S1}\\ x_2=G_{S2}{B}_{S2}\end{array}\right.---\left(16\right)$

定义卫星p对解算结果的影响△x_p为：

△x_p=x₁-x₂（17）

对II型跳变误差的估计：

II型跳变误差按照式（17）解算时，x₁无法求得。可按照以下程序进行估计，假设用于载体加速度计算的各卫星载波行为观测值的测量精度相同，则载***置确定精度为：

σ_x=σ₀·GDOP（18）

当卫星星座构型发生改变时，描述其几何构型的量GDOP会发生跳变，进而引起所确定的位置发生跳变，卫星p对解算结果的影响为：

△x^p=σ₀·△GDOP^p（19）

式（19）中，△GDOP即为星座构型改变前后GDOP值的改变量。星座中单颗卫星p对GDOP的影响可计算为：

$\mathrm{\&Delta;GDO}{P}^p=\mathrm{trace}\left(\frac{A_n{h}_p^T{h}_p{A}_n}{1-h_p{A}_n{h}_p^T}\right)---\left(20\right)$

式（20）中，A_n为星座n颗卫星对应的观测矩阵；h_p为卫星p的观测矢量；trace为求迹函数。

对III型跳变误差的估计：

III型跳变误差可分解为I型跳变误差和II型跳变误差，然后分别按照式（17）和式（19）估计。

基于上述原理和分析，本发明的用于航空重力测量中星座跳变误差的平滑消除方法为：首先对每条测线出现导航星座跳变的时刻进行分类和标记，并按照跳变的类型进行误差估计和补偿；然后，对跳变前后历元的计算位置的位置进行加权平滑补偿，以消除星座跳变对载***置解算的影响，以获取平滑的位置序列，进而求得加速度的最优解算结果。

具体而言，在本发明中采用以下平滑消除方法：

某计算时刻t，星座跳变前的位置x^a可通过上式（6）所示的最小二乘算法求解；而跳变后的解x^b可联合式（19）进行估计：

x^b=x^a+△x（21）

本发明的平滑方法就是在解x^a和x^b之间进行加权平均，即：

$\stackrel{\&OverBar;}{x}\left(t\right)=\mathrm{\&omega;}\left(t\right)x^b+[1-\mathrm{\&omega;}\left(t\right)]x^a---\left(22\right)$

式（22）中，0≤ω(t)≤1称为平滑因子，本发明采用以下方案计算该平滑因子：

$\mathrm{\&omega;}\left(t\right)=\frac{t-T}{T_L},(T\&le;t\&le;T+T_L)---\left(23\right)$

式（23）中，T为位置序列跳变时刻，T_L为数据平滑的位置序列长度，称为平滑窗口长度。由此知：

$\mathrm{\&omega;}\left(t\right)=\left\{\begin{array}{cc}0& \mathrm{if\; t}=T\\ 1& \mathrm{if\; t}=T+T_L\end{array}\right.---\left(24\right)$

加权平滑方案：根据式（22）和式（23），可得以下加权平滑方案：

$\stackrel{\&OverBar;}{x}\left(t\right)=\left\{\begin{array}{cc}{x}^a& t<T\\ \mathrm{\&omega;}\left(t\right)x_t^b+[1-\mathrm{\&omega;}\left(t\right)]x_t^a& T\&le;t\&le;T+T_L\\ x^b& t>T+T_L\end{array}\right.---\left(25\right)$

进一步，对本发明方法的平滑效果进行评估：平滑算法应消除位置序列跳变误差的影响，而同时对其他时刻的加速度提取不产生或产生可忽略的影响。通过上式（25）可知，平滑算法对t<T或t>T+T_L时段内的位置序列值没有影响，而当T≤t≤T+T_L时，二阶微分可得：

$\stackrel{\&CenterDot;\&CenterDot;}{\stackrel{\&OverBar;}{x}}\left(t\right)={\stackrel{\&CenterDot;\&CenterDot;}{x}}_t^a+\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&omega;}}\left(t\right)\&CenterDot;(x_t^b-x_t^a)---\left(26\right)$

根据式（23），可计算得：

$\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&omega;}}\left(t\right)=0,(T\&le;t\&le;T+T_L)---\left(27\right)$

由上可知，平滑算子对载体加速度的提取没有影响，并且平滑前的位置序列在t=T处发生跳变，而平滑后该处的跳变被完全消除，从而证明了本发明的有效性。

如图2所示，为本发明方法的一个具体应用实例。其中，航空重力星座跳变误差的数据处理具体包含以下步骤：

（1）GNSS接收机观测数据下载和预处理；

在航空重力测量结束后，应从基准站导航接收设备、航空机载导航接收设备中导出GNSS卫星观测数据，进行必要的数据预处理，如剔除野值、整周模糊度解算等。

（2）航空重力测量测线分段；

按照实际航空重力测量过程对测线进行分段，一般每条测线中应不包含测量飞机的转弯过程。

（3）开始数据处理，设定初始参数；

从测线起始时刻开始数据处理，并设定相应的初始参数，重置跳变补偿标记。如，设定初始时刻和跳变标记，t=0；Flag=0。

（4）利用定位算法确定载***置（t时刻）；按照式（6）所示的最小二乘算法进行载***置解算。

（5）采用本发明的方法进行跳变误差平滑补偿。

如果当前历元存在星座跳变，则置位跳变标记，判断星座跳变误差的类型，如果为I类跳变误差，按照式（17）计算跳变的大小；如果为II类跳变则按照式（19）计算跳变的大小，然后按照式（25）所示的加权平均方案对位置解算结果进行平滑处理。在测线计算结束后采用数值差分求取的载体加速度可避免跳变型误差的影响。

以上仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，应视为本发明的保护范围。

Claims (3)

1.一种用于航空重力测量中星座跳变误差的平滑消除方法，其特征在于，首先对每条测线出现导航星座跳变的时刻进行分类和标记，并按照跳变的类型进行误差估计和补偿；然后，对跳变前后历元的计算位置进行加权平滑补偿，以消除星座跳变对载***置解算的影响，以获取平滑的位置序列；

所述加权平滑补偿的过程为：

(1)某计算时刻t，星座跳变前的位置x^a可通过最小二乘算法求解，并对跳变后的解x^b进行估计：

x^b＝x^a+△x

(2)在解x^a和x^b之间进行加权平均，即：

$\stackrel{\&OverBar;}{x}\left(t\right)=\mathrm{\&omega;}\left(t\right)x^b+\&lsqb;1-\mathrm{\&omega;}\left(t\right)\&rsqb;x^a$

上式中，0≤ω(t)≤1称为平滑因子，且满足：其中T≤t≤T+T_L，T为位置序列跳变时刻，T_L为数据平滑的位置序列长度，称为平滑窗口长度；

(3)由以上步骤(1)和步骤(2)得到加权平滑方案：

$\stackrel{\&OverBar;}{x}\left(t\right)=\left\{\begin{array}{cc}{x}^a& t<T\\ \mathrm{\&omega;}\left(t\right)x^b+\&lsqb;1-\mathrm{\&omega;}\left(t\right)\&rsqb;x^a& T\&le;t\&le;T+T_L\\ x^b& t>T+T_L\end{array}\right..$

2.根据权利要求1所述的用于航空重力测量中星座跳变误差的平滑消除方法，其特征在于，如果当前历元存在星座跳变，则置位跳变标记，判断星座跳变误差的类型：

如果为I类跳变，按照下式计算跳变的大小：△x_p＝x₁-x₂；x₁为基于星座S1的解，x₂为基于星座S2的解；

如果为II类跳变，按照下式计算跳变的大小：△x^p＝σ₀·△GDOP^p，GDOP为双差定位精度衰减因子，△GDOP即为星座构型改变前后GDOP值的改变量，△GDOP^p为星座中单颗卫星p对GDOP的影响，σ₀为相位观测值，载***置确定精度σ_x＝σ₀·GDOP；然后，按照步骤(3)所示的加权平均方案对位置解算结果进行平滑处理。

3.根据权利要求2所述的用于航空重力测量中星座跳变误差的平滑消除方法，其特征在于：

所述I类跳变为导航卫星增加；由于仰角变高或者卫星重新锁定原因，导航卫星出现在测线的后半段，导航卫星可见后导致定位精度正向跳变；

所述II类跳变为导航卫星减少；由于仰角变低或者卫星被遮挡原因，导航卫星出现在测线的前半段，导航卫星不可见后导致定位精度负向跳变。