地磁导航技术具有自主性、全天时、全天候、高精度等优点，是现有导航方式的重要补充，而高精度的局部地磁基准图是实现精确地磁导航的重要基础。地磁基准图的制备以实测地磁数据为基础，并可分为2种方式，一是局部地磁场建模方法^{[1, 2, 3, 4]}，二是基于插值估计的数据加密方法。前者在描述小尺度的磁场特征，尤其是磁异常信息时，需要很高的模型阶数，计算量大，而且很难充分利用高密度的局部磁场测量数据，造成浪费；后者主要是根据未知位置周围的已知数据对其进行估计，能够充分利用所有测量数据资源，而且计算效率高，因而在精确地磁导航研究中，得到了更为广泛的应用。

黄学功^[5]、王哲^[6]、张晓明^[7]、杨功流^[8]等对制备地磁图的几种网格插值方法进行评估，实验结果显示径向基函数、克里金法具有较高的插值精度。但是对奇异性较强的局部区域的刻画误差较大。Spector^[9]等在对航磁数据谱进行统计分析的过程中证实了地磁异常具有分形特征。大多研究^{[10, 11, 12]}还局限于对整体分形特征的描述，在小尺度上，例如对文中涉及的局部地磁异常场分形特征的研究较少。文献[13]利用逐步插值校正法刻画了磁异常场在小尺度上的奇异特征，但该方法在实测数据较丰富的情况下插值效果较好。

针对上述问题，文中提出一种基于多重分形克里金的逐步插值校正法，实验表明，该方法构建的局部地磁基准图精度较高，为后续将其应用于地磁匹配导航作保障。

克里金法是插值精度较高的空间插值方法之一。它是在平均残差接近于零、估计误差方差最小的前提下，将已有数据加权线性结合得到估计值，是一种最优、线性且无偏的估计。它通过以距离为自变量的变异函数计算权值，相对于其他空间插值方法而言，更加符合实际。克里金插值法在估计待插值点属性值的过程中，不仅考虑了待插值点与邻近已知点的空间位置关系，还考虑了各相邻已知点之间的位置关系，充分利用了已知数据点的空间分布结构特征。但由于克里金插值的推导过程是建立在空间相关系数矩阵极小意义下的，这决定了它是一种滑动加权平均或低通滤波过程，因此不可避免地削弱了高频、局部及弱信号，无法重现局部区域隆起或凹陷的奇异特征。

普通克里金插值法的估值公式为

式中：Z^*为待插值点x₀处的属性值；Z_i为估计区域内1~n个采样点的实测地磁异常场强度值；λ_i为待求权重系数。

多重分形研究的关键问题之一就是测度与尺度的关系。对于曲面而言，δ尺度下的测度是以s点为中心、δ为边长的正方形邻域的质量。表示为

根据分维定义，拟合最佳逼近时，尺度与测度的关系可表示为

式中：a、b为常数。a为奇异系数，表征了待插值点小邻域内的凹凸特性。在测度—尺度的双对数坐标中，a是在最小方差意义下所拟合直线的斜率。

多重分形克里金插值法的原理为利用多重分形的理论对待插值点小邻域内克里金插值结果进行奇异性校正，即用较大邻域内地磁异常场强度的均值去估计中心小邻域内的均值，作为待插值点的插值结果。

据式(1)、(2)，将尺度为N·δ邻域及尺度为δ邻域内的测度表示为

式(3)与式(4)联立得

若克里金法的插值邻域半径为R，将圆域等效为多重分形插值的正方形区域时，为保证测度不变，取正方形边长N·δ=πR。由于插值区域内采样点是按照等步长均匀分布的，由盒计数法知，可以待插值点为中心，相邻数据点间最小距离的奇数倍为边长建立若干正方形盒子，则待插值点的面密度等价于盒子内所有采样点的地磁异常场强度平均值，再依据式(1)用面密度乘以盒子面积可求测度值m，如图 1所示。尺度δ和盒子数N是在盒计数法中设定的，即是已知值，则可依据式(2)，可求奇异系数a和校正系数N^2－α。

图 1 二维空间测度计算示意

图选项

本文选取3个实验区域对插值算法进行仿真实验，并选取4个有效性评价指标对插值后的网格数据进行评估，如表 1所示。由于在研究区域的边界缺少参考点，以至于边界处插值精度较低，因此文中在构建地磁基准图的仿真过程中均是在考虑边界效应的情况下进行的。经实验验证，第1周和第2周的参考点对边界区域的构图精度有着重要的影响，不可缺少。为了方便，将待插值数据的最外围5周作为边界区域。

表 1 有效性评价指标

指标	公式	作用
平均估计误差百分比(PAEE)		PAEE趋近于0，说明估计是无偏的
相对均方差(RMSE)		反映了随机变量之间的离散程度
均方根预测误差(RMSPE)		RMSPE越小，预测值就越接近真实值
平均绝对离差(ADD)	实测数据与插值估算数据之间的差值的绝对值加权后再取平均	避免了残差计算中正负值抵消的问题

表选项

研究区域的编号分别为4 145、4 053、4 109B，它们的地磁异常场强度数据来源于美国国家海洋和大气管理局(national oceanic and atmospheric administration，NOAA)发布的航空磁测数据。表 2列出了原始数据、基准数据及中心区域数据的经纬度范围。基准数据由原始数据网格化得到，共257×257个点，经、纬度分辨率均为0.001°。从基准数据中每隔3个点提取1个点作为待插值点，则一共可提取出65×65个点组成待插值区域(即中心区域)，其经、纬度分辨率均为0.004°。

对编号为4 053、4 145、4 109B数据的中心区域进行多重分形克里金插值，流程如图 2所示，比较多重分形克里金插值法与克里金插值法的构图精度。

图 2 多重分形克里金插值法示意

图选项

改变盒子数N的值，观察其对插值精度的影响。以4 053数据为例，实验结果见表 3。

由表 3的数据可以看出，多重分形克里金插值法的构图精度不高，原因可能有2点：

1)现有数据较稀疏，在用局部去估计整体的过程中，不可避免地导致部分分形特征丢失。

2)在某些点上克里金插值的结果已经与基准数据很接近，再用校正系数去校正，就会增大结果的误差。

为解决上述不足，提出基于多重分形克里金的逐步插值校正法。

为解决分形特征丢失的问题，现从增加采样点密度的思路出发，采用逐步插值校正的方法完成多重分形克里金插值过程，提出基于多重分形克里金的逐步插值校正法。其原理为先对分辨率为0.004°×0.004°，65×65个已知待插值数据点进行克里金插值，插值成分辨率为0.002°×0.002°的数据，共129×129个点，完成加密重构过程。再利用插值后的数据计算校正系数，完成多重分形克里金插值，如图 3所示，提取中心区域数据进行精度验证。

图 3 基于多重分形克里金的逐步插值校正法示意

图选项

观察表 3可知，对于4 053这组数据，多重分形克里金法在N=2时插值精度最高，因此仿真过程中盒子数N均取为2。同理4 145、4 109B这两组数据的盒子数分别取为2和3。盒子数是通过反复实验测得的，文中假设已知数据两点之间的距离为单位1，得到盒子数的最佳取值。仿真结果见表 4和图 4。

分析表 4的数据可知，文中方法的构图精度优于克里金插值法的构图精度。由图 4也可看出，基于逐步插值校正的多重分形克里金插值法对局部位置的刻画较细致，能够表现出更多的隆起或下陷等奇异特征，与基准数据的三维图形态更加相近，具有良好的插值效果。

图 4 区域4053三维比较图

图选项

为了更清楚地观察文中方法所构地磁基准图的准确性，将观测角度细化到每一条剖线上。图 5为4 053数据中x=82的剖线，图 6为4145数据中x=166的剖线，图 7为4 019B数据中x=106的剖线。

图 5 区域4 053中的剖线x=82

图选项

图 6 区域4 145中的剖线x=166

图选项

图 7 区域4 109B中的剖线x=106

图选项

从剖线图中更能清晰地看出，文中方法所构地磁基准图在细节上与基准数据曲线更加相似。文中方法是在克里金插值法描绘整体趋势的基础上，采用多重分形插值对局部位置逐步进行奇异性校正。

实验表明，在数据非常稀疏的情况下，文中方法相对于常用插值算法而言仍具有较高的插值精度，如表 5所示。表中数据的分辨率为0.016°，中心区域共7×7个采样点。

综上所述，基于多重分形克里金的逐步插值校正法的构图精度较高，局部位置奇异特征较明显，在局部地磁基准图构建方面具有显著的优势。

文中为改善克里金插值法的低通滤波性，引入多重分形理论，增加对局部位置的奇异性描述，还原输入信号中更多的高频信息，研究了多重分形克里金插值方法。由于稀疏的已知采样点直接插值会导致部分分形特征丢失，插值精度降低，因此文中从增加参考数据密度的角度出发，在考虑边界效应问题的条件下，提出基于多重分形克里金的逐步插值校正法。实验证明该方法相对于克里金插值法得到了较大改善，适合于构建高精度的局部地磁基准图。