CN113048976A

CN113048976A - 一种双磁参量坐标反演定位方法及装置

Info

Publication number: CN113048976A
Application number: CN202110180232.XA
Authority: CN
Inventors: 李子圆; 于化鹏; 沈同圣
Original assignee: National Defense Technology Innovation Institute PLA Academy of Military Science
Current assignee: National Defense Technology Innovation Institute PLA Academy of Military Science
Priority date: 2021-02-08
Filing date: 2021-02-08
Publication date: 2021-06-29
Anticipated expiration: 2041-02-08
Also published as: CN113048976B

Abstract

本申请公开了一种双磁坐标反演定位方法及装置，该方法包括：获取当前位置的地磁场总强度F和磁倾角I；利用非线性关系模型解算与所述地磁场总强度F和磁倾角I对应的当前位置的地理经纬度；以及根据当前位置与目标位置的地理经纬度之间的差值，判断是否到达目标位置。本申请通过建立地磁坐标与经纬度的非线性映射关系，只需测量所在地的总磁场强度和磁倾角即可得到该位置的地理经纬度，具有很好的实用性以及产品便携性。

Description

一种双磁参量坐标反演定位方法及装置

技术领域

本申请属于定位导航技术领域，具体涉及利用地磁场进行地球物理导航技术，尤其是涉及一种双磁坐标反演定位方法及装置。

背景技术

利用地磁场导航具有隐蔽性高，自主性强的特点，是主动导航的重要分支。目前，针对磁场导航的研究主要集中在基于磁图的匹配算法和基于仿生进化的磁参量搜索方法。

然而，基于磁图匹配算法需要提前获取目标区域的地球磁场图，包括磁参数的7个分量，并且导航精度与磁图分辨率呈正相关，计算负荷也会随磁图分辨率提高而迅速增加。

尽管基于磁参数的目标搜索方法不需要先验磁图，但这种方法导航实时性差，具有边定位边探索的特点，如果需要多个参数在目标点同时收敛也是相当困难，使得这种方法即便在理想环境下，通常也存在一定的导航偏差。

发明内容

根据本申请的第1方面，公开了一种双磁坐标反演定位方法，包括：

获取当前位置的地磁场总强度F和磁倾角I；

利用非线性关系模型解算与所述地磁场总强度F和磁倾角I对应的当前位置的地理经纬度；以及

根据当前位置与目标位置的地理经纬度之间的差值，判断是否到达目标位置；

其中，所述非线性关系模型将地理经纬度坐标Y表示为泰勒展开式系数矩阵A与双磁坐标矩阵X的乘积，所述双磁坐标由地磁场总强度F、磁倾角I表示。

在其他的一些示例中，通过以下方法确定所述泰勒展开式系数矩阵A：

构建训练模型数据库；

将所述区域内的各提取点的磁坐标及地理经纬度坐标带入所述关系模型，得到n阶泰勒展开式系数矩阵的初始值；

利用最小二乘法解算迭代，得到泰勒展开式系数矩阵的优化值，进而确定最终的非线性关系模型。

在其他的一些示例中，将当前位置的地磁场总强度F和磁倾角I代入所述最终的非线性关系模型，得到当前位置的地理经纬度。

在其他的一些示例中，所述地理经纬度坐标表示为λ+φi，λ、

分别为地理经纬度，i为复数因子。这里，通过将地理经纬度坐标表示为复数形式，将二维数据转换为一维数据，大大方便了非线性模型的解算过程。

在其他的一些示例中，如果当前位置的地理经纬度

与目标位置的地理经纬度

之间的差值分别小于等于设定阈值，则认为已到达目标位置。

在其他的一些示例中，还包括对利用非线性关系模型解算出的当前位置的经度进行卡尔曼滤波处理。

在其他的一些示例中，所述非线性关系模型为：

式中，左边的i为复数因子，当前位置的地理经纬度为

右边的n代表泰勒展开阶数，A_i,j代表泰勒展开式系数矩阵，F₀代表总磁场强度展开点，I₀代表磁倾角展开点，F代表当前位置磁场强度，I代表当前位置磁倾角。

根据本申请的第2方面，公开了一种双磁坐标反演定位装置，包括：

磁场数据采集单元，用于获取当前位置的地磁场总强度F和磁倾角I；

计算单元，用于利用非线性关系模型解算与所述地磁场总强度F和磁倾角I对应的当前位置的地理经纬度；

判断单元，用于根据当前位置与目标位置的地理经纬度之间的差值，判断是否到达目标位置；

在其他的一些示例中，所述计算单元还包括滤波器，该滤波器用于对利用所述非线性关系模型解算得到的当前位置的经度进行卡尔曼滤波处理。

在其他的一些示例中，所述非线性关系模型为：

式中，左边的i为复数因子，当前位置的地理经纬度为

本申请的方法主要利用地磁场总强度F和磁倾角I反演经纬度信息，通过建立(F，I)组合坐标与经纬度的非线性映射关系，然后借助最小二乘法求解得到非线性关系表达式系数，最后实现只需测量所在地的总磁场强度和磁倾角即可得到该位置的地理经纬度。整个过程不涉及磁图匹配或磁场数据库存储，实时测量当地磁场参量便可立即得到当地经纬度信息，具有很好的实用性以及产品便携性。

通过以下参照附图对本申请的示例性实施例的详细描述，本申请的其它特征及其优点将会变得清楚。

附图说明

此处所说明的附图用来提供对本申请的进一步理解，构成本申请的一部分，本申请的示意性实施例及其说明用于解释本申请，并不构成对本申请的不当限定。在附图中：

图1为根据本申请实施例的双磁坐标反演定位方法示意性工作流程；

图2为根据本申请实施例的确定系数矩阵的示意性工作流程；

图3-5分别为地磁场日变化、磁场风暴以及磁异常三种环境下对双磁坐标反演定位方法的精度影响情况；

图6为根据本申请实施例的双磁坐标反演定位装置结构示意图。

具体实施方式

下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。以下对至少一个示例性实施例的描述实际上仅仅是说明性的，决不作为对本申请及其应用或使用的任何限制。基于本申请中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本申请保护的范围。

除非另外具体说明，否则在这些实施例中阐述的部件和步骤的相对布置、数字表达式和数值不限制本申请的范围。对于相关领域普通技术人员已知的技术、方法和设备可能不作详细讨论，但在适当情况下，所述技术、方法和设备应当被视为授权说明书的一部分。在这里示出和讨论的所有示例中，任何具体值应被解释为仅仅是示例性的，而不是作为限制。因此，示例性实施例的其它示例可以具有不同的值。

为了解决现有利用地磁导航方法中对先验数据库过度依赖,以及多磁参量搜索进化算法收敛困难等问题，本申请通过建立磁场总强度F、磁倾角I与地理经纬度之间的非线性映射关系，进而通过实时测量所在地磁参量F和I得到双磁参量定位坐标，借助所述非线性关系将双磁参量坐标定位系统转化为经纬度定位系统。

图1示出了根据本申请第一实施方式的双磁参量坐标反演定位方法的示意性流程。如图1所示，该方法具体包括如下步骤：

步骤102，获取当前位置的地磁场总强度F和磁倾角I；

步骤104，利用非线性关系模型解算与所述地磁场总强度F和磁倾角I对应的当前位置的地理经纬度；

所述非线性关系模型用于表示双磁坐标与地理经纬度坐标(地理坐标)之间的非线性映射关系，所述双磁坐标由地磁场总强度F、磁倾角I表示；

本申请中，双磁坐标F、I与地理经纬度坐标λ、

之间的关系模型为：

其中，式(1)左边的i为复数因子；式(1)右边的n代表泰勒展开阶数，A_i,j代表泰勒展开式系数矩阵，F₀代表总磁场强度展开点，I₀代表磁倾角展开点，F代表当前位置磁场强度，I代表当前位置磁倾角。

式(1)右边是对当前位置磁场强度F和磁倾角I的n阶泰勒展开式，具体为：

令F＝F-F₀，I＝I-I₀,则上式为：

表达成矩阵形式为：

这里，Α＝[A_0,0,A_1,0,...,A_n,n]，X＝[F⁰I⁰,F¹I⁰,...,FⁿIⁿ]^T。

设：

Y＝λ+φi

则式(1)可简化为：

Y＝Α·X (2)

即，将地理经纬度坐标Y表示为泰勒展开式系数矩阵A与双磁坐标矩阵X的乘积。

通过以下方法确定所述泰勒展开式系数矩阵A：

步骤1042，构建训练模型数据库；

在地球表面选择预定区域，在该预定区域中选择一定数量的提取点，从地磁模型中提取该预定区域的总磁场强度F和磁倾角I，以及各提取点的地理经纬度，作为训练模型数据库。可以理解，该训练模型数据库包含足够多的区域信息。

可选地，设置各提取点在训练模型数据库中的磁场网格分辨率为0.05°，空间距离为5km。

例如，选择一条美国纽约(41.04N,73.07W)至法国巴黎(48.79N,2.43E)的跨大西洋航线，如图3所示，并从世界地磁模型WMM2020中提取(40.64N,73.07W)，(41.44N,73.07W),(49.19N,2.43E),(48.39N,2.43E)四点围成的矩形区域的总磁场强度F和磁倾角I以及每个提取点的地理经纬度，构建训练模型数据库。

此外，在导航系统中，载体安装有用于提供航向和姿态信息的设备，例如陀螺仪等。

步骤1044，将所述区域内的各提取点的磁坐标及地理经纬度坐标带入所述关系模型，得到n阶泰勒展开式系数矩阵的初始值；

例如，采用4阶泰勒展开，则所述关系模型为：

则4阶泰勒展开式系数矩阵A的初始值如下表所示：

A<sub>0,0</sub>	A<sub>1,0</sub>	A<sub>1,1</sub>	A<sub>2,0</sub>	A<sub>2,1</sub>	A<sub>2,2</sub>	A<sub>3,0</sub>	A<sub>3,1</sub>
								6.294	8.116	9.412	-7.461	9.143	6.983	8.277	-0.293
A<sub>3,2</sub>	A<sub>3,3</sub>	A<sub>4,0</sub>	A<sub>4,1</sub>	A<sub>4,2</sub>	A<sub>4,3</sub>	A<sub>4,4</sub>
								8.680	-9.077	2.647	6.006	3.575	-8.058	-6.262

步骤1046，利用最小二乘法解算迭代，得到泰勒展开式系数矩阵的优化值，进而确定最终的非线性关系模型。

例如，设置最小二乘法解算迭代次数为100，经过迭代优化，得到最优的系数矩阵A，进而确定最终的关系模型。

将当前位置的地磁场总强度F和磁倾角I代入所述最终的非线性关系模型，得到当前位置的地理经纬度。

例如，离开出发位置第一次测得总磁场强度为53711.56nT，磁倾角为70.03°，将其代入所述关系模型中(其中，系数矩阵A为经过优化的系数矩阵)中，得到当前位置的经纬度分别为-72.4641W，41.0919N。

步骤106，根据当前位置与目标位置的地理经纬度之间的差值，判断是否到达目标位置。

设出发点位置与目标位置的地理经纬度坐标

其中λ_o、λ_t分别代表出发位置和目标位置的经度，

分别代表出发点位置和目标位置的纬度。

如果当前位置的地理经纬度

与目标位置的地理经纬度

之间的差值分别小于等于设定阈值，即：

|λ-λ_t|≤ε_λ

其中，ε_λ为经度误差阈值，

为纬度误差阈值，则认为已到达目标位置；否则，实时获取当前位置的地磁场总强度F和磁倾角I，利用所述最终的非线性关系模型计算当前位置的地理经纬度，再次进行上述判断。

在一些示例中，经度误差阈值ε_λ为0.5°，纬度误差阈值

为0.1°。

由此可见，本申请提出的方法更像是一个带有磁力计的“盲盒”，“盲盒”测量地球磁场参数，输出的是所在地的经纬度坐标。与全球定位系统(GPS)不同的是，GPS需要在有信号的区域接收卫星信号，而本申请提出的“盲盒”只要能测到地磁场就能给出定位信息，幸运的是，在地球上任何位置都有磁场分布。

在本申请的第二实施方式中，针对经度对磁场变化敏感，易受磁场噪声干扰的情况，引入扩展卡尔曼滤波算法，提高经度反演的鲁棒性。以下描述中，对于未明确说明的内容，认为沿用第一实施方式的方法。

假设载体在k点测得当前位置的总磁场强度F和磁倾角I，利用所述最终的非线性关系模型计算得到当前位置的经度为λ_k，作为观测量可以写成：

Z_k＝λ_k

在卡尔曼算法中状态向量可以被写成：

R(k)＝[x_k,v_k]^T

其中，x_k，v_k分别代表经度和经度方向速度。

本申请中，离散时间过程的状态方程为：

R(k+1)＝AR(k)+Bu(k)+Gw(k)

观测方程为：

Z(k)＝HR(k)+v(k)

式中，A为状态转移矩阵，G为噪声驱动矩阵，H为观测矩阵，w(k)，v(k)分别为过程噪声和观测噪声。

为降低模型复杂度，本申请将载体设为匀速运动，即u(k)＝0。这里，w(k)，v(k)均为相互独立的白噪声，分布满足如下关系：

p(w)～N(0,Q)

p(v)～N(0,M)

其中，Q，M分别是过程噪声和观测噪声的方差。

可选地，设置Q＝0.0001*diag([1,1])，M＝0.0785。

然后，基于k时的观测值更新k+1时的状态向量R(k+1)，实现基于卡尔曼滤波的双磁坐标反演定位方法。

在本实施例中，所选仿真预定区域为长方形，长边沿东西走向，短边沿南北走向，所选区域经度跨度为75.50°，纬度跨度为7.75°，在构建的双磁坐标反演定位模型中经度的变化幅度更大，表现在泰勒多项式的系数A对经度反演的精确控制相比于纬度难度更大。此外，在所选预定区域总磁场F和磁倾角I的等值线大致沿东西走向，因此在仿真路线上相邻两位置点的磁场元素变化不大情况下，经度会有较大变动，表现在受噪声干扰时反演的经度误差比纬度误差要大。

以下分析地球磁场噪声对所述方法的影响，例如，地磁场日变化、磁场风暴以及磁异常三种环境下对双磁坐标反演定位方法的精度产生影响。

为了分析地磁场日变化对双磁坐标反演定位的影响，在局部区域(46.59-47.79N，19.07-7.57W)加入连续均匀分布噪声，作为实时测量到总磁场强度F和磁倾角I日变干扰。

例如，设置大磁场日变干扰环境中，总磁场强度变化幅度60nT，磁倾角I变化幅度0.15°。小磁场日变干扰环境中，总磁场强度变化幅度30nT，磁倾角I变化幅度0.10°。

在所选局部区域测量当地磁场数据，并利用所述第二实施方式的方法计算出经纬度。仿真结果如图3所示，图中灰色直线是磁场测量点的真实经纬度，黑色圆代表小磁场日变干扰下用本申请所述方法反演得到的位置，灰色圆代表小磁场日变干扰下用本申请所述方法反演得到的位置。

为了分析磁风暴引起的地磁场变化对双磁坐标反演定位的影响，在局部区域(23.57W-2.43E，46.14-48.79N)加入噪声，作为实时测量总磁场强度F和磁倾角I的日变干扰。

可以设置两种不同强度的磁风暴干扰，第一种是总磁场强度变化幅值70nT，第二种是总磁场强度变化幅值150nT。两种磁风暴随时间的变化特征相同，不同的仅是总磁场强度F变化量的最大值。

下式描绘了第一种磁风暴环境下噪声随时间分布情况，式中k代表时间。

在所选局部区域测量当地磁场数据，并利用本申请第二实施方式的方法计算出经纬度。仿真结果如图4所示，图中灰色直线是磁场测量点的真实轨迹，黑色圆代表小磁风暴干扰下用本申请所述方法反演轨迹，灰色圆代表磁风暴干扰下用本申请所述方法反演轨迹。

为了分析地磁异常对双磁坐标反演定位的影响，在区域(20.37-14.87W，45-48N)加入噪声，作为实时测量到总磁场强度F和磁倾角I日变干扰。

可选地，从NOAA发布的2020年WMM模型中，提取设计航线的磁异常值作为实测量总磁场强度F和磁倾角I的数据噪声。

在所选局部区域测量当地磁场数据，并利用本申请第二实施方式的方法计算出经纬度。仿真结果如图5所示，图中灰色直线是磁场测量点的真实轨迹，黑色圆代表磁异常干扰下用本申请所述方法反演轨迹。

从图中可以看出，经过卡尔曼滤波算法处理，本申请提供的方法能够有效降低地球磁场噪声的影响，从而得到更加准确的地理经纬度坐标。

根据本申请的第三实施方式，还提供了一种双磁坐标反演定位装置200，如图6所示，该装置包括：

磁场数据采集单元201，用于获取当前位置的地磁场总强度F和磁倾角I；

计算单元202，用于利用非线性关系模型解算与所述地磁场总强度F和磁倾角I对应的当前位置的地理经纬度；

判断单元203，用于根据当前位置与目标位置的地理经纬度之间的差值，判断是否到达目标位置。

其中，所述非线性关系模型为：

式中，左边的i为复数因子；右边的n代表泰勒展开阶数，A_i,j代表泰勒展开式系数矩阵，F₀代表总磁场强度展开点，I₀代表磁倾角展开点，F代表当前位置磁场强度，I代表当前位置磁倾角。

其中，所述计算单元202还包括滤波器2021，该滤波器用于对利用所述非线性关系模型解算得到的经度数据进行卡尔曼滤波处理，以提高地理经纬度的计算精度。

最后应当说明的是，以上实施例仅用以说明本申请的技术方案而非对其限制。尽管参照较佳实施例对本申请进行了详细的说明，所属领域的普通技术人员应当理解，依然可以对本申请的具体实施方式进行修改或者对部分技术特征进行等同替换。而不脱离本申请技术方案的精神，其均应涵盖在本申请请求保护的技术方案范围当中。

Claims

1.一种双磁坐标反演定位方法，其特征在于，包括：

获取当前位置的地磁场总强度F和磁倾角I；

2.根据权利要求1所述的双磁坐标反演定位方法，其特征在于，通过以下方法确定所述泰勒展开式系数矩阵A：

构建训练模型数据库；

3.根据权利要求2所述的双磁坐标反演定位方法，其特征在于，将当前位置的地磁场总强度F和磁倾角I的代入所述最终的非线性关系模型，得到当前位置的地理经纬度。

4.根据权利要求2所述的双磁坐标反演定位方法，其特征在于，所述地理经纬度坐标表示为λ+φi，λ、

分别为地理经纬度，i为复数因子。

5.根据权利要求1所述的双磁坐标反演定位方法，其特征在于，如果当前位置的地理经纬度

与目标位置的地理经纬度

6.根据权利要求1-5任一项所述的双磁坐标反演定位方法，其特征在于，还包括对利用非线性关系模型解算出的当前位置的经度进行卡尔曼滤波处理。

7.根据权利要求1-6任一项所述的双磁坐标反演定位方法，其特征在于，所述非线性关系模型为：

式中，左边的i为复数因子，当前位置的地理经纬度为

8.一种双磁坐标反演定位装置，其特征在于，包括：

9.根据权利要求8所述的双磁坐标反演定位装置，其特征在于，所述计算单元还包括滤波器，该滤波器用于对利用所述非线性关系模型解算得到的当前位置的经度进行卡尔曼滤波处理。

10.根据权利要求8或9所述的双磁坐标反演定位装置，其特征在于，所述非线性关系模型为：

式中，左边的i为复数因子，当前位置的地理经纬度为