一种基于地磁匹配的地磁日变修正方法

摘要

本发明公开了一种基于地磁匹配的地磁日变修正方法。针对地磁匹配导航定位中的地磁日变效应，本发明提出基于FMI方法与多维特征量匹配对地磁日变效应进行修正，包括如下步骤：1)选取地磁匹配特征量；2)基于FMI方法从地磁匹配特征量的实测数据中提取各地磁匹配特征量的日变量；3)通过地磁匹配特征量的实测数据减去与其对应的日变量，建立地磁实时图；4)基于多维特征量匹配对地磁实时图和地磁基准图进行匹配；5)输出匹配结果。本发明能较为有效地消除地磁匹配导航定位中的地磁日变效应，能显著提高地磁匹配导航定位精度。

说明书

技术领域

本发明涉及一种基于地磁匹配的地磁日变修正方法，属于高精度地磁匹配导航的关键 技术领域。

背景技术

地磁导航是当前国际导航领域的一大研究热点，具有高度自主性、隐蔽性、无辐射、 能耗低等诸多优点，能弥补惯导系统位置误差长期累积的不足，也能克服GPS等卫星有源 导航系统易被侦测和受干扰等缺陷，已成为隐蔽性运动载体无源导航的重要手段，广泛应 用于潜艇、导弹和航天器辅助导航等领域，展现出巨大的军事价值和应用前景。然而，地 磁匹配导航定位精度不高严重制约了其在相关领域的工程应用。

地磁日变是影响地磁匹配导航精度的重要因素之一。在地磁匹配导航研究中，许多学 者指出了地磁日变的重要性，但没有开展深入的研究。目前的研究表明地磁日变对地磁导 航的影响是不容忽视的，但存在以下问题：

①多数研究对地磁日变影响只进行了简单的探讨；

②部分研究直接将日变以噪声的形式加入到仿真中；

③部分研究完全忽略地磁日变效应。

传统日变场的分析方法有日均值法和子夜均值法。其中日均值法模糊了日变白天变化 大、夜间变化小的特点。子夜均值法以子夜均值为基线值，但基线值的选择具有很大的人 为影响因素，容易给匹配定位带来新的误差源。FMI（Finish Meteorological Institute，芬兰 气象局）方法是一种基于三天的地磁场X、Y分量分钟值观测数据为数据源计算中间一天K 指数的计算方法，此方法由该气象局提出。与前述两者相比，FMI方法是对观测数据的直 接拟合，不需选取基线值，减少了人为误差，所提取的日变更为客观、真实，具有较大的 适用性。

发明内容

本发明所要解决的技术问题是提供一种能用于地磁匹配导航定位中的地磁日变修正 方法，基于FMI方法提取地磁日变成分，在多维地磁特征量匹配中消除地磁日变效应，从 而提高地磁匹配导航定位精度。

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

一种基于地磁匹配的地磁日变修正方法，包括以下步骤：

步骤1，选取地磁匹配特征量，所述地磁匹配特征量为三个地磁分量，包括北向分量 X、东向分量Y、垂直分量Z；

步骤2，基于FMI方法从地磁匹配特征量的实测数据中提取各地磁匹配特征量的日变 量；

步骤3，通过地磁匹配特征量的实测数据减去与其对应的日变量，建立地磁实时图；

步骤4，基于多维特征量匹配对地磁实时图和地磁基准图进行匹配；

步骤5，输出匹配结果。

作为本发明的基于地磁匹配的地磁日变修正方法的进一步优化方案，所述步骤2涉及 的实测数据是由磁传感器采集的；所述基于FMI方法从实测数据中提取日变量包括以下步 骤：

步骤2.1，根据地磁匹配特征量的扰动幅度a分别确定与各地磁匹配特征量相对应的K 指数的初始值K_X0、K_Y0、K_Z0，取K₀=max(K_X0,K_Y0,K_Z0)，所述扰动幅度a为每3小时时段 内地磁匹配特征量的实测数据的最大值和最小值间的差值；

步骤2.2，地磁匹配特征量的时均值由该小时及其前后各(m+n)分钟的数据平均值来计 算，然后利用5次谐波函数对时均值进行拟合，生成S_R曲线；所述n＝K³³；m与地方时有 关：地方时3点到6点和18点到21点，m取60；地方时7点到17点，m取0；其余时间， m取90；

步骤2.3，分别计算地磁匹配特征量的实测数据与S_R曲线的残差序列，由残差序列的 3小时扰动幅度a分别确定与各地磁匹配特征量相对应的K指数的中间值K_Xm、K_Ym、K_Zm， 取K_m=max(K_Xm,K_Ym,K_Zm)；

步骤2.4，重复步骤2.2和步骤2.3的计算一次，再次得到K指数的中间值

K_m=max(K_Xm,K_Ym,K_Zm)，此K_m即为K指数的最终值；

步骤2.5，由K指数的最终值，代入步骤2.2计算，生成最终S_R曲线；

步骤2.6，计算最终S_R曲线与地球主磁场和磁异常间的差值，得到地磁匹配特征量的 日变量。

作为本发明的基于地磁匹配的地磁日变修正方法的进一步优化方案，所述步骤4涉及 的地磁基准图是指基于地磁模型建立的基准数据库，或基于测量资料绘制的地磁图。

作为本发明的基于地磁匹配的地磁日变修正方法的进一步优化方案，步骤4涉及的多 维特征量匹配是指通过3个地磁匹配特征量按照先后顺序分别进行匹配，即以北向分量X、 东向分量Y、垂直分量Z这三个独立特征量进行匹配，具体包含以下步骤：

步骤4.1，以北向分量X作为第一匹配特征量，根据子夜均值附近北向分量X的日变 大小设定阈值λ_X，以北向分量X的实时数据B_X在地磁基准数据库R_x，y，z={B_x,B_y,B_z}中搜 索满足给定地磁匹配算法条件的匹配点，设该匹配结果为点集R_y,z={B′_y,B′_z}；其中，B_x， B_y，B_z分别表示三个地磁分量上的地磁基准数据，B′_y,B′_z分别表示在基准数据库中搜索出 的东向分量、垂直分量上的地磁基准数据；

步骤4.2，以东向分量Y作为第二匹配特征量，根据子夜均值附近东向分量Y的日变 大小设定阈值λ_Y，以东向分量Y的实时数据B_Y在集合R_y,z中搜索满足给定地磁匹配算法条 件的匹配点，设该匹配结果为点集R_z={B″_z}，其中，B″_z表示在集合R_y，z中进一步搜索出的 垂直分量上的地磁基准数据；

步骤4.3，以垂直分量Z作为第三匹配特征量，以垂直分量Z的实时数据B_Z在集合R_z中搜索满足给定地磁匹配算法条件的匹配点，该点即为最佳匹配位置。

作为本发明的基于地磁匹配的地磁日变修正方法的进一步优化方案，所述地磁匹配算 法为常用的相关分析匹配算法，包括平均平方差算法MSD、平均绝对差算法MAD、平方 差算法SD、绝对差算法AD、积相关算法PROD。

作为本发明的基于地磁匹配的地磁日变修正方法的进一步优化方案，所述地磁匹配算 法为平均绝对差算法MAD，具体如下：

$\mathrm{MAD}(P,Q)=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}|p_i-q_i|,$

其中，P为实时数据序列，Q为基准数据序列，实时数据序列P和基准数据序列Q的 长度相同，N表示实时数据序列P或基准数据序列Q的数据序列长度，p_i表示实时数据序 列P中第i个值，q_i表示基准数据序列Q中第i个值，i、N均为正整数。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：

本发明将FMI方法应用到地磁匹配导航中，基于该方法提取地磁静日期间地球变化磁 场的最主要成分——地磁日变，并建立地磁实时图，通过多维地磁特征量进行匹配，消除 地磁匹配中的地磁日变效应，从而提高地磁匹配导航定位精度。

同时，本发明能有效地消除地磁匹配中的日变效应，显著地减小了匹配误差。

附图说明

图1是日变修正方法流程图。

图2是基于FMI方法提取日变量的流程图。

图3是多维特征量匹配流程图。

图4是2005年7月25日广州台站X、Y、Z分量的观测值曲线、拟合曲线和日变曲线。

图5是2005年7月25日兰州台站X、Y、Z分量的观测值曲线、拟合曲线和日变曲线。

图6是2005年7月25日北京十三陵台站X、Y、Z分量的观测值曲线、拟合曲线和日 变曲线。

图7是2009年7月2日广州台站X、Y、Z分量的观测值曲线、拟合曲线和日变曲线。

图8是2009年7月2日兰州台站X、Y、Z分量的观测值曲线、拟合曲线和日变曲线。

图9是2009年7月2日北京十三陵台站X、Y、Z分量的观测值曲线、拟合曲线和日 变曲线。

图10是2005年7月25日广州、兰州和北京十三陵三台站日变修正前后的匹配误差 图。

图11是2009年7月2日广州、兰州和北京十三陵三台站日变修正前后的匹配误差图。

具体实施方式

下面结合附图对本发明的技术方案做进一步的详细说明：

本实施例是基于地磁台站观测数据、FMI方法和多维特征量匹配进行的。以分布在不 同纬度上的三个地磁台站（广州GZH、兰州LZH和北京十三陵BMT）在太阳活动高年（2005 年）和太阳活动低年（2009年）的地磁观测数据分钟值模拟地磁实测数据。三个台站的地 理位置见下表：

地磁台站 纬度(°) 经度(°) 海拔（m） 广州 23.093N 113.343E 11 兰州 36.087N 103.845E 1560 北京十三陵 40.3N 116.2E 183

对实测数据进行预处理：所选用的地磁台站数据已对环境干扰进行了处理，但需要进 行地方时校正，并剔除坏值。

参照图1，基于地磁匹配的地磁日变修正方法的分析步骤如下：

步骤1，选取地磁匹配特征量：地磁场有七个分量，分别为北向分量X、东向分量Y、 垂直分量Z、磁偏角D、水平强度H、磁倾角I和总强度F。其中，X、Y、Z三个分量是独 立的。本实施例选取X、Y、Z三分量作为地磁匹配特征量。

步骤2，基于FMI方法从地磁匹配特征量的实测数据中提取各地磁匹配特征量的日变 量：参照图2所示提取流程，具体步骤为：

步骤2.1，分别从各地磁台站提取X、Y、Z三分量的分钟值观测数据；根据观测数据 的3小时扰动幅度a（3小时时段内的最大值和最小值间的差值）分别确定各匹配特征量 对应的K指数（描述单个地磁台站3小时时段内地磁扰动强度的指数）的初始值K_X0、K_Y0、 K_Z0，取K₀=max(K_X0,K_Y0,K_Z0)；所述K指数从0到9共分10级，由3小时扰动幅度a确 定：a<3nT，K取0；3≤a<6nT，K取1；6≤a<12nT，K取2；12≤a<24nT，K取3； 24≤a<40nT，K取4；40≤a<70nT，K取5；70≤a<120nT，K取6；120≤a<200nT， K取7；200≤a<300nT，K取8；a≥300nT，K取9。

步骤2.2，地磁X、Y、Z三分量的时均值由该小时及其前后各(m+n)分钟的数据平均值 来计算，然后利用5次谐波函数(对时间序列x=[x₁,x₂，...,X_N]，x包含N个数据，其5次 谐波函数为$T_5\left(t\right)=\frac{a_0}{2}+\underset{j=1}{\overset{5}{\Sigma }}{a}_j\cos \left(2\pi f_{j}t\right)+b_j\sin \left(2\pi f_{j}t\right),$其中，$a_0=\frac{x_1+x_N}{2},$$a_j=\frac{2}{N}\underset{t=1}{\overset{N}{\Sigma }}{x}_t\cos \left(2\pi f_{j}t\right),$$b_j=\frac{2}{N}\underset{t=1}{\overset{N}{\Sigma }}{x}_t\sin \left(2\pi f_{j}t\right),$$f_i=\frac{j}{N},$t＝1,2,....,N，j=1,2,...,5)对时 均值进行拟合，生成S_R曲线（Solar Regular variations，规则日变化）；所述n=K³³，m与 地方时有关：地方时3点到6点和18点到21点，m取60；地方时7点到17点，m取0； 其余时间，m取90。

步骤2.3，分别计算地磁X、Y、Z三分量的观测值与S_R曲线的残差序列，由残差序列 的3小时扰动幅度a分别确定与各地磁匹配特征量相对应的K指数的中间值K_Xm、K_Ym、 K_Zm，取K_m=max(K_Xm,K_Ym,K_Zm)；

步骤2.4，重复步骤2.2和步骤2.3的计算一次，再次得到K指数的中间值

K_m=max(K_Xm,K_Ym,K_Zm)，此时的K_m即为K指数的最终值；

步骤2.5，由K指数的最终值，代入步骤2.2计算，生成最终S_R曲线；

步骤2.6，计算最终S_R曲线与地球主磁场和磁异常间的差值，得到地磁三分量X、Y、 Z的日变量。

参照图4至图9，分别为2005年7月25日和2009年7月2日广州、兰州和北京十三 陵三个地磁台站X、Y、Z三分量的日变量处理结果。

步骤3，通过地磁X、Y、Z三分量的观测数据分别减去对应的日变量建立地磁实时图。

步骤4，基于多维特征量匹配对地磁实时图和地磁基准图进行匹配。

本实施例基于多维特征量匹配对地磁实时图和基准图进行匹配，具体如下：

采用IGRF11（the11^th generation International Geomagnetic Reference Field，第11代国际参考地磁 场）+磁异常，分别以广州、兰州和北京十三陵三个台站所在位置为中心建立2005年7月 25日和2009年7月2日的地磁基准图，位置参数采用经纬度，取网格大小为1000×1000， 网格精度为0.001°×0.001°。

参照图3，本实施例所涉及的多维特征量匹配具体是指：

(1)选择MAD（平均绝对差）作为匹配算法：

其中，P为实时数据序列，Q为基准数据序列，实时数据 序列P和基准数据序列Q的长度相同，N表示实时数据序列P或基准数据序列Q的数据 序列长度，p_i表示实时数据序列P中第i个值，q_i表示基准数据序列Q中第i个值，i、N 均为正整数。

本发明适用的相关分析匹配算法还有很多，常规的比如：MSD（Mean Square Deviation， 平均平方差）、MAD（Mean Absolute Deviation，平均绝对差）、SD（Square Deviation，平 方差）、AD（Absolute Deviation，绝对差）、PROD（Product Correlation，积相关）等，都 属于常用的匹配算法，在这里就不一一赘述。

(2)匹配过程：

①以北向分量X作为第一匹配特征量，设定阈值λ_X（地磁日变具有在子夜附近几乎 不变的特点，因此，阈值的设定参考子夜均值附近北向分量X的日变大小），以北向分量X 的实时数据B_X在基准地磁数据库R_x，y，z={B_x,B_y,B_z}中搜索使B_X满足MAD(B_X,B_x)<λ_X的匹配点，设该匹配结果为点集R_y，z={B′_y,B′_z}；其中，B_x，B_y，B_z分别表示X、Y、Z三 个地磁分量上的地磁基准数据，B′_y,B′_z分别表示在基准数据库中搜索出的东向分量、垂直 分量上的地磁基准数据；

②以东向分量Y作为第二匹配特征量，设定阈值λ_Y（地磁日变具有在子夜附近几乎 不变的特点，因此，阈值的设定参考子夜均值附近东向分量Y的日变大小），以东向分量Y 的实时数据B_Y在集合R_y，z中搜索使B_Y满足MAD(B_Y,B′_y)<λ_Y的匹配点，设该匹配结果为点 集R_z={B″_z}；其中，B″_z表示在集合R_y，z中进一步搜索出的垂直分量上的地磁基准数据；

③以垂直分量Z作为第三匹配特征量，以垂直分量Z的实时数据B_Z在集合R_z中搜索 使B_Z满足min(MAD(B_Z,B″_z))的点，该点即为最佳匹配位置。图10和图11分别是2005年7 月25日和2009年7月2日广州、兰州和北京十三陵三个地磁台站日变修正前后的匹配误 差图。参阅图10和图11可知，本发明在地方时6点至18点之间修正效果尤为明显，匹 配定位精度能提高一个数量级。

(1)综合来看，地磁匹配中的日变效应显著，随太阳活动性、地方时、纬度而分布。日 变修正前后的匹配误差曲线具有明显差异：修正前，匹配误差波动大，变化剧烈；修正后， 匹配误差曲线波动小，变化较平缓。尤其在地方时6点到18点之间，日变修正后匹配定 位误差明显减小，匹配定位精度得到显著提高。在图10中，广州台站白天14时，未进行 日变修正前经度误差达到0.4°，而日变修正后经度误差仅为0.05°左右，匹配定位精度 提高了约一个数量级。

(2)从图10和图11对比来看，日变修正效果在太阳活动低年要优于太阳活动高年。这 是因为，在太阳活动高年，地磁活动剧烈，除地磁日变外，还存在一些较大的扰动，如地 磁脉动、地磁钩扰等。图10中可以看出，在太阳活动高年，在低纬度上的日变效应尤为 显著。经日变修正后，广州台站纬度方向上误差仍可达到0.1°，经度方向上达0.2°。而 图11中，在太阳活动低年，经日变修正后，三个台站纬度方向上误差控制在0.005°左右， 经度方向上则控制在0.05°左右。

(3)三个台站地理位置不同，但修正效果都比较明显，说明该方法具有较大的适用性。