一种应用于地磁导航中的地磁场模型误差补偿方法

摘要

本发明提供了一种应用于地磁导航中的地磁场模型误差补偿方法，该方法是通过分析地磁导航滤波算法进入收敛期后的地磁场测量数据与地磁模型输出数据之间的误差，找到该误差随地理纬度的变化规律；利用傅里叶变换对误差数据进行频率辨识，然后基于主要频率进行数据拟合，获得地磁模型误差的数学模型；以此数学模型为基础，结合地磁误差数据与地理纬度之间的变化规律，在滤波过程中对地磁场模型误差进行实时预测与补偿，最终实现提高地磁导航精度的目的。

说明书

技术领域

本发明涉及卫星自主导航技术领域，具体涉及一种应用于地磁导航中的地磁场模型误差补偿方法。

背景技术

地磁导航通过星上磁强计的测量值与地磁模型的输出值(预测值)进行比较而得到轨道修正信息，利用该修正信息可以实现轨道的确定，进而达到航天器自主导航的目的。地磁导航是近地航天器的一种重要的自主导航方式，可实现航天器位置、速度以及姿态信息的自主确定。与传统的GPS卫星导航相比，具有抗干扰能力强、隐蔽性强的优势；与惯性导航相比，地磁导航误差不随时间累积，并且导航系统体积小、功耗低。因此在近地轨道卫星、水下航行器等方面具有很好的应用前景。

利用地磁场为近地航天器进行导航时，需要利用地磁模型来估计预测值，所以地磁场模型是实现地磁导航不可或缺的基础，地磁模型与实际地磁场之间的差别，即地磁模型精度是制约导航定位精度的核心因素，所以提高地磁模型的精度是实现高精度地磁导航的最重要的途径。目前，常用的地磁模型，如国际地磁场参考模型(InternationalGeomagnetic Reference Field，IGRF)、世界地磁场模型(World Magnetic Model，WMM)是采用高斯球谐函数，建立了地磁场与地理位置(经度、纬度、高度)数学关系表达式，这种方法的优点是能够得到一个合理的、全面反映磁场的宏观分布的模型，但是忽略了局部的磁异常信息。另外由于地磁场复杂的空间结构和时间演化特性，尽管人们在研究地磁场方面做了很多努力，地磁模型每五年更新一次，仍然难以得到全球的高精度磁场模型。受到观测资料、分析方法以及模型阶数的限制，国际地磁场模型IGRF的误差可能达到250nT，显然不能满足高精度卫星导航的需求。

基于以上分析可知，研究地磁场模型误差补偿方法是地磁导航研究的重要方向。而国内在这个方向的研究较少，国外研究由于军事应用的背景，不会公开其研究成果。从已公开的资料来看，常用的模型补偿方法是通过建立高斯球谐系数以及系数时间变化量(其中表示高斯球谐系数，表示球谐系数的变化量，下标1表示阶数，上标1表示次数；)的状态方程，并将这些系数作为状态变量，在导航过程中实时修正，从而达到提高导航精度的目的。但是，由于模型的球谐系数多达N*(N+2)个(N为高斯球谐系数的阶数)，这么多状态量的实时估计对星载计算机的计算能力是一个很大的挑战，并且此方法也不适用于工程上的应用。

发明内容

本发明的目的在于，考虑到高精度地磁导航和其工程应用两方面因素，本发明提出一种简单有效的地磁模型误差补偿方法。该方法利用傅里叶变换进行量测残差频率的在线辨识，然后建立残差与地理特征之间的函数表达式，最终达到补偿地磁模型误差，提高地磁导航精度的目的。

为实现上述目的，本发明提供了一种应用于地磁导航中的地磁场模型误差补偿方法，该方法是通过分析地磁导航滤波算法进入收敛期后的地磁场测量数据与地磁模型输出数据之间的误差，找到该误差随地理纬度的变化规律；利用傅里叶变换对误差数据进行频率辨识，然后基于主要频率进行数据拟合，获得地磁模型误差的数学模型；以此数学模型为基础，结合地磁误差数据与地理纬度之间的变化规律，在滤波过程中对地磁场模型误差进行实时预测与补偿，实现提高地磁导航精度的目的。

所述的应用于地磁导航中的地磁场模型误差补偿方法，具体包括：

步骤1)当地磁导航滤波算法稳定后，计算当前卫星轨道周期内地磁场测量数据与地磁场模型输出数据之间的量测残差；

步骤2)利用傅里叶变换辨识量测残差的频率信息；

步骤3)依据量测残差在频率下的表达式拟合量测残差与地理特征之间的关系，建立残差数学模型；

步骤4)将下一卫星轨道周期内地理纬度数据输入到残差的数学模型中，计算获得地磁场模型在下一卫星轨道周期内的量测残差的补偿信息。

作为上述技术方案的进一步改进，所述步骤2)中傅里叶变换的表达式为：

其中：A_i,ω_i,分别表示各余弦分量的幅值、频率和相位；t表示时间；C为常值项；ΔB_x/z表示北东下坐标系的x和z两个方向的量测残差；所述地磁场测量数据包括x、y和z三个矢量方向的测量数据。

作为上述技术方案的进一步改进，所述的步骤2)与步骤3)之间还包括：设定幅值的阈值为4nT，利用阈值与步骤2)获得的频率信息中各余弦分量的幅值进行比较，选择幅值大于该阈值的余弦分量，作为步骤3)中参与拟合的频率信息。

作为上述技术方案的进一步改进，所述步骤3)中的残差数学模型表示为：

其中，t＝a*θ+b，η_i＝aω_i；θ为纬度值；a为斜率，b为常值。

作为上述技术方案的进一步改进，所述步骤4)中量测残差的补偿计算公式表示为：

其中，u,d分别表示在卫星运行周期内轨道上升和下降两个阶段，ΔB₂表示后一轨道周期内补偿后的量测残差，B_测,B_模分别表示后一轨道周期内地磁场强度的实测值与地磁场模型的输出值，分别表示轨道上升和下降两个阶段中三个矢量方向的量测残差补偿值，θ₂表示后一轨道周期内的纬度值，f₁(θ₂)表示将后一轨道周期内的纬度值θ₂代入前一轨道周期所建立的残差数学模型。

本发明的一种应用于地磁导航中的地磁场模型误差补偿方法优点在于：

1、导航精度高

本发明能够在模型预测的基础上，对地磁场模型的固有误差进行补偿。基于卫星的地磁场实测数据进行的仿真表明，应用本发明的补偿方法后，使得单一地磁导航的位置与速度精度分别可以提高63％和58％。(经仿真测试，在相同条件下，基于地磁模型误差在线补偿地磁导航精度为2.72km和3.26m/s，传统地磁导航精度为7.4km和7.75m/s)

2、工程应用前景广阔

本发明所提出的补偿方法不需要增加其他的测量设备，仅利用滤波过程中的量测残差就能够对模型误差进行预测，因此具有较强的工程应用价值。

附图说明

图1a为地磁场模型在BX方向上的误差周期性变化曲线。

图1b为地磁场模型在BY方向上的误差周期性变化曲线。

图1c为地磁场模型在BZ方向上的误差周期性变化曲线。

图2为本发明提供的一种应用于地磁导航中的地磁场模型误差补偿方法流程图。

图3为地磁场模型强度F的全球分布图。

图4a为地磁场模型在BX方向上的卫星轨道上升阶段与下降阶段的误差变化曲线。

图4b为地磁场模型在BZ方向上的卫星轨道上升阶段与下降阶段的误差变化曲线。

图5为本发明的方法中建立残差数学模型的操作流程图。

图6为本发明的方法中地磁场模型误差在线预测与补偿的流程图。

图7为利用本发明的方法进行地磁场模型误差补偿的地磁导航流程图。

图8为地磁导航位置误差变化曲线。

图9a为地磁场模型在BX方向上的残差数据区间。

图9b为地磁场模型在BZ方向上的残差数据区间。

图10a为无地磁场模型误差补偿的位置误差曲线。

图10b为有地磁场模型误差补偿的位置误差曲线。

具体实施方式

下面结合附图和实施例对本发明所述的一种应用于地磁导航中的地磁场模型误差补偿方法进行详细说明。

本发明是针对地磁导航中地磁场模型固有误差提出的在线补偿方法。首先基于卫星的实测数据获得地磁场模型的三个矢量方向固有误差。从图1a、图1b和图1c中可以看出，BX、BY、BZ三个方向的误差都具有周期性变化的特点。在地磁导航中，量测残差信息中包含全部的模型误差。如图2所示，本发明提供的一种应用于地磁导航中的地磁场模型误差补偿方法，包括以下步骤：

步骤1)当地磁导航滤波算法稳定后，计算当前卫星轨道周期内地磁场测量数据与地磁场模型输出数据之间的量测残差。

步骤2)利用傅里叶变换辨识量测残差的频率信息。

步骤3)依据量测残差在频率下的表达式拟合量测残差与地理特征(纬度)之间的关系，建立残差数学模型；该残差数学模型能够以纬度为自变量预测下一周期内的模型误差，同时可以补偿到地磁导航算法中。

步骤4)将下一卫星轨道周期内地理纬度数据输入到残差的数学模型中，输出值作为地磁场模型在后一卫星轨道周期内的量测残差的补偿信息，完成地磁场模型的误差预测与补偿。

上述误差补偿方法的具体实施步骤如下：

1、利用傅里叶变换进行量测残差的频率辨识

任何周期函数都可以表示为不同频率的正弦或者余弦函数之和的形式，每个正弦函数或者余弦函数分别乘以不同的系数，然后进行线性组合便可以构成原来的周期函数，甚至非周期函数(该函数下的面积是有限的)也可以用正弦或余弦乘以加权函数的积分来表示。不同频率的正弦函数或者余弦函数的线性组合就是傅里叶级数。假设f(x)为任意周期函数，则其用傅里叶级数表示为：

其中：n表示序列；a₀,a_n,b_n都是实常数；x为自变量。傅里叶变换是将时域信号转换为频域信号，不同性质的信号如常值信号、周期信号在频域下会有不同的特性。通过傅里叶变换可以得到一组信号的幅频特性曲线和相频特性曲线，而幅频特性曲线和相频特性曲线包含了本发明所需要的所有信息。傅里叶变换公式表示为：

而且傅里叶变换满足线性叠加原理，若f₁(x),f₂(x)的傅里叶变换分别为F₁(ξ),F₂(ξ)，则有：

af₁(x)+bf₂(x)＝aF₁(ξ)+bF₂(ξ)(3)

a,b分别表示系数，ξ表示频率，i表示复数，x为自变量。傅里叶变换是信号处理领域普遍采用的一种方法，能够将一组信号中不同频率的部分区分出来。

为此，基于上述傅里叶变换原理，在本发明中通过傅里叶变换对地磁场模型量测残差进行辨识，获得其频率信息。如图5所示，以滤波算法稳定后某周期内的轨道上升阶段为例，对BX、BZ两个方向的量测残差进行傅里叶变换，可以得到由不同频率余弦分量和常值项组成的表达式。考虑到噪声的影响，可近似认为噪声是由较小幅度的分量组成，这里设定幅值的阈值为4nT进行判断。利用阈值与步骤2)获得的频率信息中各余弦分量的幅值进行比较，选择幅值大于该阈值的余弦分量组成数学表达式，作为步骤3)中参与拟合的频率信息，即：

其中：A_i,ω_i,分别表示各余弦分量的幅值、频率和相位；t表示时间；C为常值项；ΔB_x/z表示北东下坐标系的x和z两个方向的量测残差，所述地磁场测量数据包括x、y和z三个矢量方向的测量数据。

北东下坐标系是指坐标原点在卫星的质心，X轴和Y轴位于当地铅垂面内，其中X轴指向北，Y轴指向东，Z轴垂直向下指向地心，并与X/Y轴构成右手直角坐标系。通过对数据分析发现，量测残差在By方向的变化量较小，且周期性不明显。因此，如果对Bx/By/Bz三个方向进行模型误差的补偿，那么仿真效果较差，对导航精度的提高没有效果。为此，在本发明的方法中，可只针对残差变化量较大、且周期性明显的两个方向Bx/Bz进行预测和误差补偿，而不再对By方向上的误差进行计算。

由于BX/BZ方向量测残差随纬度的变化由余弦曲线的趋势，所以一般情况下，只需要3-4项余弦分量就可以表示残差的频域特性。

2、地磁场模型误差拟合

地磁场具有复杂的空间结构和时变特性，虽然它以地心为中心，但是地磁场并不是关于地轴对称的模型。从图3地磁场模型强度F的全球分布图也可以明显看出其分布的不规律性。同样对于地磁场模型固有误差来讲，它的变化规律也不是对称的，如图4a、4b所示，地磁场模型在BX、BZ方向上，卫星轨道上升阶段与下降阶段两个阶段相比，误差变化并不相同。因此，考虑到曲线拟合的准确性，将每个周期分为两个阶段进行表达式的拟合。

如图5所示，卫星在轨道运行过程中，每个周期内纬度θ与时间t之间是线性变化的关系，即：

t＝a*θ+b(5)

其中：a为斜率，b为常值。且相邻两个周期内，纬度基本保持不变，利用这个关系，将量测残差表示为关于纬度的函数，这个表达式可以用来预测相邻周期内上升或者下降阶段的模型误差。

将上述公式(5)代入公式(4)，即获得残差数学模型表示为：

其中：η_i＝aω_i；

利用上述曲线拟合方法可以得到量测残差与时间的关系，而在卫星运动过程中，每个周期内纬度与时间基本保持线性关系，因此最终能够获得滤波算法稳定后任意周期内量测残差与纬度之间的函数表达式。

3、地磁场模型误差在线预测与补偿

如图6所示，地磁场模型误差在线预测与补偿的具体过程包括：(1)选择滤波算法稳定后某个周期为初始1T周期，利用上述曲线拟合的步骤，分别对该周期内轨道上升和下降阶段内的BX、BZ两个方向的量测残差进行拟合，得到如下四个函数表达式：

其中，u,d分别表示轨道上升和下降两个阶段；θ₁表示该周期内卫星的纬度；ΔB_1x,ΔB_1z分别表示该周期内，北东下坐标系下x和z两个方向的量测残差。

(2)2T周期内地磁场模型量测残差的补偿及拟合

以上式(7)中四个拟合表达式为基础，对下一个周期内的滤波进行地磁场模型的预测与补偿，则补偿后的量测残差如下式(8)所示：

其中，ΔB₂表示后一轨道周期内补偿后的量测残差，B_测,B_模分别表示2T周期内地磁场强度的实测值与地磁模型的输出值；分别表示轨道上升和下降两个阶段中三个矢量方向的量测残差补偿值，θ₂为2T周期内的纬度信息；f₁(θ₂)表示将后一轨道周期内的纬度值θ₂代入前一轨道周期所建立的残差数学模型。通过对图1中误差曲线的分析，已知模型误差具有周期性变化的规律，所以以基于1T的拟合表达式和2T的纬度，能够预测该周期内的模型误差。然后将其作为模型的固有误差对地磁场模型进行补偿。

对2T周期滤波后的量测残差进行曲线拟合，可以预测并补偿3T周期内的模型误差。但是这里需要注意的是，由于在2T周期内已对模型固有误差进行了补偿，所以补偿后的量测残差ΔB_2T与正常滤波时的量测残差相比，数值上有了很大的变化。如果仍然以量测ΔB_2T为基础进行曲线拟合，肯定无法得到正确的误差预测表达式。所以这里为了保证量测残差拟合的准确，当2T周期误差补偿后，需以为基础进行曲线拟合，即将2T周期的原始量测残差数据(未补偿的数据)作为3T周期内模型误差预测的表达式。

(3)3T周期内补偿及拟合

利用2T周期内补偿预测、补偿以及残差拟合方法，对3T周期模型误差进行补偿；以此类推，利用上述拟合方法对4T、5T…NT周期进行相同的误差补偿操作。

实施例一

在本实施例中，利用本发明的方法进行地磁场模型误差补偿的地磁导航流程如图7所示：首先建立系统的测量方程和以地球固联坐标系下的轨道动力学方程为基础的状态方程，然后利用卫星的姿态信息，将磁强计的测量数据进行坐标转换，与地磁场模型的输出作差后得到量测残差，并作为导航修正信息同状态方程、量测方程输入到所设计的扩展卡尔曼滤波器中，最终得到精度较高的导航信息。本发明提出的补偿方法是用于稳定后的滤波器中，如图7所示，具体实施包括获取量测残差数据，提取频率信息并建立残差数学模型，以及在线补偿模型误差三个步骤。

1、获取滤波算法稳定后的残差基准数据

分别建立系统状态方程和量测方程，以卫星上磁强计实测数据为测量信息，设计扩展卡尔曼滤波器，然后设定仿真初始条件(如表1所示)。得到如图8所示的地磁导航位置误差变化曲线，选取滤波算法稳定后的某个周期作为残差数据基准区间，如图8中标注的区间所示。图9a和图9b分别示出了对应基准区间内BX和BZ方向上的量测残差。

表1滤波初值

Table 1 Initial parameters of filter

扩展卡尔曼滤波器可以分为时间序列更新和量测序列更新，如下所示：

时间序列更新表示为：

量测序列更新表示为：

上式中：为状态量预测值，为状态量的滤波值，分别为状态预测值和状态滤波值的误差协方差矩阵，Z(k)为磁强计的测量信息，K_k为增益矩阵，I为单位矩阵；Φ(k,k-1)表示状态转移矩阵；R(k),Q(k)分别为测量噪声协方差矩阵和系统噪声协方差矩阵。

2、辨识量测残差的频率信息并建立残差数学模型

选取基准周期内的BX和BZ方向上的量测残差，利用傅里叶变换的方法提出残差的主要频率信息。由于地磁场的分布关于地球自转轴不是对称的，卫星在升轨和降轨两个阶段内，地磁场模型误差的变化规律也是不相同。所以这里将分为升轨和降轨两个阶段进行残差频率信息的提取。另外，因为残差信息中包含小部分噪声的干扰，设定判断条件，舍去幅值小于4nT的函数项，从而保留能够反映残差变化规律的主要余弦项。记升轨阶段为u，降轨阶段为d，则基准区间内BX和BZ方向上的残差数据进行傅里叶变换的结果如下式(11)-(14)所示：

由以上各式可知，经过傅里叶变换和剔除小幅值噪声两个处理步骤后，基准周期内的残差的表达式都可以由若干项余弦函数分量表示，例如：在5.68*cos(0.00224*t+1.67)表达式中，5.68表示该项余弦信号的幅值，括号内0.00224表示余弦信号的频率，t表示时间，1.67表示余弦信号的相位。每个表达式右端最后一项，如-63.54表示信号的常值项。

根据卫星的运行特点，在升轨u和降轨d两个阶段，地理纬度θ与时间t都是线性变化的关系，如式(15)和(16)所示：

t(u)＝898.5θ(u)+1406(15)

t(d)＝-898.5θ(d)+1403(16)

将式(15)代入式(11)和(12)，式(16)代入式(13)和(14)，可以得到Bx/Bz两个方向上的量测残差与地理纬度θ之间的关系表达式，也叫做残差的数学模型，如下式(17)-(20)所示：

上式中：θ_u,θ_d分别表示升轨和降轨阶段的地理纬度。

3、在线预测并补偿地磁场模型误差

根据BX/BZ方向上的量测残差具有周期性变化的规律，可以利用上述步骤得到的基准周期残差数学模型，在滤波过程中实时预测相邻周期内的残差信息。假定相邻周期内的地理纬度分别为将其代入式(17)-(20)中量测残差与纬度的关系式，可以表示如下：

将预测的残差补偿到地磁场模型中，则扩展卡尔曼滤波算法的量测更新方程式(10)变换为：

升轨阶段：

降轨阶段：

通过实施上述步骤，可以实现地磁导航过程中地磁场模型误差补偿的目的。使用本发明所提出的方法后，导航误差曲线如图10a、10b所示，图10a、10b给出了补偿区间为2T-6T时的位置误差对比曲线图。从图中可以得出，相比于无地磁场模型误差补偿的位置误差曲线，经过误差补偿后的位置精度可以由补偿前的7.4km提高到2.7km，位置精度由补偿前的7.35m/s提高到3.26m/s。由此可知，利用本发明提供的地磁场模型误差补偿方法，能够有效改善地磁导航的精度。

最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。