CN102841385B - 一种基于多重分形克里金法的局部地磁图构建方法 - Google Patents

Abstract

本发明公开了一种基于多重分形克里金法的局部地磁图构建方法，包括步骤一：去除实测数据中的主磁场成分；步骤二：选定估值区域；步骤三：拟合高斯模型，计算估值区域内各控制点间以及各控制点与待插值点之间的半方差；步骤四：计算估值区域内各测点的权重系数；步骤五：利用多重分形理论推算待插值点邻域内的测度表达式；步骤六：建立待插值点在其小邻域内的多重分形克里金法插值方程；步骤七、采用交叉验证法验证插值的精度并构建地磁图。本发明充分考虑地磁异常场的空间相关性及其在待插值点邻域内的奇异特征，根据实测数据进行精确的插值，使得插值结果更加符合实际地磁异常场的空间变化趋势，为地磁导航构建高精度、高分辨率的局部地磁图。

Description

一种基于多重分形克里金法的局部地磁图构建方法

技术领域

本发明属于地磁场数据处理领域，具体涉及一种基于多重分形克里金法的局部地磁图构建方法。

背景技术

高精度的地磁基准图的获取是实现精确地磁导航的关键。地磁异常场是地磁场的重要组成部分，由地壳表面分布的矿产、岩石、人造磁场等产生，分布很稳定，几乎不随时间变化，通常只在近地空间存在；它的空间结构复杂，波形尺度小，具有局部、随机、异常的特征，非常适合作为地磁导航的参照。

由于我国地域广阔，开展广泛的陆地及海洋磁测工作需要耗费巨大的成本，因此依靠直接测量获取细致的地磁异常数据难度很大。在已有数据的基础上通过插值来获得高精度、高分辨率的地磁异常网格数据是对实测工作的一项重要补充。

陈丽等研究了双线性插值方法在构建高精度地磁图中的应用，实验证明该方法计算量小，实时性较好，但它只适合应用在与位置关系密切且在大范围内具有明显变化规律的地球主磁场的处理中，而对于局部、随机特性比较明显的地磁异常场，该方法并不适用；王金玲等讨论了反距离加权法在空间数据插值中的应用，它根据搜索范围内实测点和待插值点的距离作为加权依据，忽视了实测点之间的位置关系，如果实测点的位置相对待插值点具有明显的重复性，那么插值结果将出现较大的误差；张晓明等单纯使用克里金法进行了局部地磁图的构建，并通过实验证明该方法能够比较准确的反映局部地区地磁信息，但是这种方法在充分描述地磁数据空间相关性的同时，没有考虑地磁异常数据的小尺度奇异性，即忽略了地磁异常场的高频信息。

克里金法是一种用于空间插值的地理统计方法，充分吸收了地理统计思想，认为任何在空间连续变化的属性是非常不规则的，不能简单用平滑的数学函数进行模拟，而采用随机表面可以恰当的描述，它是一种无偏估计，是空间相关意义下的最优（最小方差）估计。克里金法假设某种属性的空间变化既不是完全随机也不是完全确定的，可以表示为3种主要成分的和：与恒定均值和趋势有关的结构成分；与空间变化有关的随机变量，即区域性变量；与空间无关的随机噪声。地磁场中的地球主磁场、地磁异常场及扰动磁场3个主要成分与克里金法中假定的结构性成分、随机变量以及随机噪声具有一一对应关系。因此克里金法较为准确的反应地磁异常场的空间相关特性规律，可以用于地磁异常的空间插值。

但是克里金法建立在空间相关系数矩阵极小意义下的求解方式，决定了它是一个空间滑动平均或低通滤波过程，高频、局部与弱信号在相关系数矩阵的方差中占得比例很小，难以反映地磁异常场在小尺度上的奇异特征。而多重分形理论可以弥补这个缺点，它是由Mandelbrot在研究湍流时提出来的，专门应用于研究物理量和其他量在几何支撑上的奇异性分布，它是定义在分形上，由多个标度指数的奇异测度所组成的集合。Spector等人通过对航磁数据谱的统计分析，表明磁异常具有分形特征。

发明内容

本发明的目的是为了解决上述问题，提出一种基于多重分形克里金法的局部地磁图构建方法，本发明在已有较稀疏的实测地磁异常数据的基础上，针对地壳表面空间结构复杂且非常稳定的地磁异常场，采用多重分形克里金法对实测地磁数据进行插值，建立高精度、高分辨率的局部地磁异常网格数据，进而为实现地磁导航构建局部地磁图。本发明不仅能充分的表征磁异常的区域性变化，同时能保留更多的高频信息。

一种基于多重分形克里金法的局部地磁图构建方法，包括以下几个步骤：

步骤一、去除实测数据中的主磁场成分；

判断实测地磁数据中是否包含地球主磁场成分，如果包含，利用地磁场模型计算各测点主磁场强度，并在实测数据中将其去除；

步骤二、选定估值区域；

根据实测地磁数据的分辨率以及地磁异常场的空间特性，以待插值点为中心，选定估值区域；

步骤三：选定高斯模型描述地磁异常场的变异特性，拟合高斯模型中的未知参数，并利用此模型计算估值区域内各控制点间以及各控制点与待插值点之间的半方差；

在估值区域内，利用实测地磁数据，计算不同的分离距离对应的变异函数值，从而得到分离距离与变异函数值序列对，根据此序列对，利用拟合法推算出高斯模型中的未知参数，并利用此模型计算估值区域内各实测点之间以及实测点和待插值点之间的半方差；

步骤四、计算估值区域内各测点的权重系数；

利用步骤三中得出的半方差值，通过解空间相关系数矩阵方程求出各实测点的权重系数；

步骤五、利用多重分形理论推算待插值点邻域内的测度表达式；

在估值区域内，利用测度与尺度的在双对数坐标系中的线性关系，用最小二乘拟合法得出插值点邻域内的分形维数，进而推算出以带插值点为中心的小正方形（邻域）的测度与整个估值区域测度的关系式；

步骤六、建立待插值点在其小邻域内的多重分形克里金法插值方程；

利用步骤四中求出的估值区域内各实测点的权重系数，求出整个估值区域测度表达式，并根据步骤五中得到的关系式，推算出待插值点在其小邻域内的多重分形克里金法插值方程；

步骤七、采用交叉验证法验证插值的精度并构建地磁图；

在测区内每一个有观测值的地方，将此观测值暂时去除，用以此点为中心建立的估值区域内的其他观测值以及本专利发明的插值方法，估计这一观察点的值，然后将暂时去除的观察值放回，重复以上步骤，对测区所有的观察点进行估值。并采用均方根预测误差、平均估计误差百分比、相对均方差三个指标对插值结果进行评估。设定评价指标，如果插值结果符合指标要求，则根据插值后的地磁数据构建地磁基准图；否则，返回步骤二，重新进行插值。

本发明的优点在于：

本发明的优点在于在插值过程中，采用多重分形理论精确的反映出地磁异常场在待插值点邻域内的奇异特征，并利用克里金方法充分考虑了地磁异常场的空间相关特性，能够在已有实测数据的基础上，通过插值得到高精度高、分辨率的地磁异常网格数据。

附图说明

图1是本发明的方法流程图；

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

本发明是一种基于多重分形克里金法的局部地磁图构建方法，流程如图1所示，包括以下几个步骤：

步骤一：去除实测数据中的主磁场成分；

判断实测地磁数据中是否包含地球主磁场成分，如果包含，利用地磁场模型计算各测点主磁场强度，并在实测数据中将其去除，具体为：

（1）判断实测地磁数据中是否包含地球主磁场成分；

实测地磁数据为N个测点序列，测点序列为(λ_i,g_i),其中，1≤i≤N，λ_i,g_i分别为测点的地理经度、地理纬度、地磁测量值，g_i为矢量；

因为主磁场在地表处的强度为30000—70000nT，占地磁场总量的95%以上，地磁异常场强度占地磁场强度的4％—5%，在1000nT-4000nT之间，所以，通过观察实测数据的大小，判断实测地磁数据中是否包含有主磁场成分，当实测地磁强度大于30000nT时，实测地磁数据中包含主磁场成分，当其位于1000nT-4000nT之间时，实测地磁数据中不包含主磁场成分。

（2）如果实测地磁数据中包含主磁场成分，利用地磁场模型计算各测点主磁场强度，并在实测数据中将其去除；

根据国际地磁参考场模型（IGRF），利用截断阶数法获取测点在北向、东向、竖直（向下）方向上的地磁场分量X、Y、Z，X、Y、Z通过式（1）、式（2）、式（3）获取，具体为：

$X=\underset{n=1}{\overset{\∞}{\mathrm{\Σ}}}\underset{m=0}{\overset{n}{\mathrm{\Σ}}}{\left(\frac{a}{r}\right)}^{n+2}(g_n^m\cos \mathrm{m\λ}+h_n^m\sin \mathrm{m\λ})\frac{{\∂P}_n^m\left(\mathrm{\θ}\right)}{\∂\mathrm{\θ}}---\left(1\right)$

$Y=\underset{n=1}{\overset{\∞}{\mathrm{\Σ}}}\underset{m=0}{\overset{n}{\mathrm{\Σ}}}{\left(\frac{a}{r}\right)}^{n+2}(g_n^m\sin \mathrm{m\λ}-h_n^m\cos \mathrm{m\λ})\frac{{\mathrm{mP}}_n^m\left(\mathrm{\θ}\right)}{\sin \mathrm{\θ}}---\left(2\right)$

$Z=-\underset{n=1}{\overset{\∞}{\mathrm{\Σ}}}\underset{m=0}{\overset{n}{\mathrm{\Σ}}}(n+1){\left(\frac{a}{r}\right)}^{n+1}(g_n^m\cos \mathrm{m\λ}+h_n^m\sin \mathrm{m\λ})P_n^m\left(\mathrm{\θ}\right)---\left(3\right)$

其中，a为地球半径；r是地球表面或表面以上空间中的计算点与地心的距离；θ是地理余纬度（设为地理纬度，则有）；λ是地理经度； 为高斯球谐系数；n、m分别为球谐系数的阶和次；为n阶m次缔合勒让德函数。

IGRF模型反映的主要是地球主磁场的变化趋势，在公式（1）、（2）、（3）中将n＞13的高阶项截掉，得公式（4）、（5）、（6），剩余低阶项部分可以准确的描述地球主磁场。

$X=\underset{n=1}{\overset{13}{\mathrm{\Σ}}}\underset{m=0}{\overset{n}{\mathrm{\Σ}}}{\left(\frac{a}{r}\right)}^{n+2}(g_n^m\cos \mathrm{m\λ}+h_n^m\sin \mathrm{m\λ})\frac{{\∂P}_n^m\left(\mathrm{\θ}\right)}{\∂\mathrm{\θ}}---\left(4\right)$

$Y=\underset{n=1}{\overset{13}{\mathrm{\Σ}}}\underset{m=0}{\overset{n}{\mathrm{\Σ}}}{\left(\frac{a}{r}\right)}^{n+2}(g_n^m\sin \mathrm{m\λ}-h_n^m\cos \mathrm{m\λ})\frac{{\mathrm{mP}}_n^m\left(\mathrm{\θ}\right)}{\sin \mathrm{\θ}}---\left(5\right)$

$Z=-\underset{n=1}{\overset{13}{\mathrm{\Σ}}}\underset{m=0}{\overset{n}{\mathrm{\Σ}}}(n+1){\left(\frac{a}{r}\right)}^{n+1}(g_n^m\cos \mathrm{m\λ}+h_n^m\sin \mathrm{m\λ})P_n^m\left(\mathrm{\θ}\right)---\left(6\right)$

将实测地磁数据中包含主磁场成分的测点(λ_i,g_i)依次带入式（4）、式（5）、式（6）中，分别得到各个测点的地球主磁场三分量X_Ri,Y_Ri,Z_Ri，即得此测点的地球主磁场矢量

${\stackrel{\→}{g}}_{\mathrm{Ri}}=(X_{\mathrm{Ri}},Y_{\mathrm{Ri}},Z_{\mathrm{Ri}}).$

将g_i与g_Ri作矢量差，得出第i个测点的地磁异常矢量g_Ui,相应的地磁异常强度w_i＝|g_Ui|。

步骤二：选定估值区域；

根据实测地磁数据的分辨率以及地磁异常场的空间特性，以待插值点为中心，选定估值区域，具体为：

地磁异常场具有明显的空间相关特性，估值区域内必须包含充足的测点才能对待插值点的地磁异常值进行准确的估计。在测区内，以待插值点为中心，建立以D为尺度（边长）的正方形区域作为估值区域。D不宜过小也不宜过大，过小会给估值结果带来较大误差，过大会带来不必要的计算量，影响实时效果。

设测区内实测点的平均间距为d（以航空磁测为例，在经、纬度方向上4″≤d≤6″），取D＝10d~15d。

步骤三：选定高斯模型描述地磁异常场的变异特性，拟合高斯模型中的未知参数，并利用此模型计算估值区域内各实测点之间以及各实测点和待插值点之间的半方差；

在估值区域内，利用实测地磁数据，计算不同的分离距离对应的变异函数值，从而得到分离距离与变异函数值序列对，根据此序列对，利用拟合法推算出高斯模型中的未知参数，并利用此模型计算估值区域内各实测点之间以及实测点和待插值点之间的半方差，具体为：

（1）选定高斯模型描述地磁异常场的变异特性，如（7）式：

$\mathrm{\γ}\left(h\right)=\left\{\begin{array}{c}0,h=0\\ C_0+C_1(1-e^{-{(h/a)}^2}),h>0\end{array}\right.---\left(7\right)$

其中：C₀为块金值，C₀+C₁为基台值，即h为无穷大时的变异函数值，h为样本间的距离,变程为C₀、C₁、a均为未知量。

（2）利用实测地磁数据，计算不同的分离距离对应的变异函数值，从而得到分离距离与变异函数值序列对，根据此序列对，利用拟合法推算出高斯模型中的未知参数；

根据实测地磁数据，获取估值区域内实测点的实测地磁异常数据序列为w_a1、w_a1、...、w_aP，（w_i表示的是整个实测区域内地磁异常数据序列中的一个值，而w_a1、w_a2、...w_aP是估值区域内的地此异常数据序列，知识全部实测地磁异常数据的一部分）P为估值区域内实测点的个数，根据式（8）得到估值区域内的变异函数值：

${\mathrm{\γ}}^{*}\left(h^{\′}\right)=\frac{1}{{2N}_{h^{\′}}}\underset{i=1}{\overset{N_{h^{\′}}}{\mathrm{\Σ}}}{[w_i^{*}-w_i]}^2---\left(8\right)$

其中，h′为分离距离，N_h′是估值区域内距离小于h′实测点对数，和w_j为距离小于h′的一对实测数据。

（3）利用拟合法推算出高斯模型中的未知参数，并利用此模型计算估值区域内各实测点之间以及实测点和待插值点之间的半方差；

通过式（8），根据估值区域内各实测点的之间的距离，通过调整h的值，可以得到变异函数值与分离距离序列对，这一序列对同时满足式（7），由此可以拟合出高斯模型中的三个未知参数C₀、C₁、a，并可求得变程

利用已知的高斯模型计算出估值区域内各实测点之间的半方差γ_ij(1≤i,j≤P)以及各实测点与待插值点之间的半方差γ_0i(1≤i≤P)。

步骤四：计算估值区域内各测点的权重系数；

利用步骤三中得出的半方差值，通过解空间相关系数矩阵方程求出各实测点的权重系数，具体为：

设估值区域内各实测点的权重系数为λ₁,λ₂,...,λ_P，空间相关系数矩阵方程如下：

$\left[\begin{array}{cccc}{\mathrm{\γ}}_{11}& ...& {\mathrm{\γ}}_{1P}& 1\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ {\mathrm{\γ}}_{P1}& ...& {\mathrm{\γ}}_{\mathrm{PP}}& 1\\ 1& ...& 1& 0\end{array}\right]\left[\begin{array}{c}{\mathrm{\λ}}_1\\ .\\ .\\ .\\ {\mathrm{\λ}}_P\\ \mathrm{\μ}\end{array}\right]=\left[\begin{array}{c}{\mathrm{\γ}}_{01}\\ .\\ .\\ .\\ {\mathrm{\γ}}_{0P}\\ 1\end{array}\right]---\left(9\right)$

其中：μ为拉格朗日乘数，γ_jk为式（5）中求出的各实测点之间的半方差，γ_0j为式（5）求出的实测点与待插值点之间的半方差，1≤j,k≤P，通过解此方程（9）可以得到各测点的权重系数。

步骤五：利用多重分形理论推算待插值点邻域内的测度表达式；

在估值区域内，利用测度与尺度的在双对数坐标系中的线性关系，用最小二乘拟合法得出插值点邻域内的分形维数，进而推算出以待插值点为中心的小正方形（邻域）的测度与整个估值区域测度的关系式，具体为：

根据分维定义，设F为估值区域内实测数据集，M(F)为尺度l下F的测度，当l→0时，有：

$M\left(F\right)\≅{\mathrm{cl}}^a---\left(10\right)$

其中，c为常数，a为数据集的分维。测度M(F)与尺度l服从指数定律，在双对数坐标系中表现为一条直线，即

$\ln M\left(F\right)\≅\ln c+a\ln l---\left(11\right)$

采用最小二乘方法拟合式（11），通过双对数图中直线的斜率即可求出a。

在估值区域内，以待插值点为中心，建立尺度为L的小正方形（小邻域），并令D＝kL（D估值区域的边长，k＞1），则该邻域的测度M_L和估值区域的测度M_kL分别如（12）式和（13）式所示：

M_L＝m_L*L²＝c*L^a    (12)

M_kL＝m_kL*(kL)²＝c*(kL)^a    (13)

式中，m_L和m_kL分别为小正方形和大正方形的平均测度。

综合式（12）和式（13）可得待插值点邻域内的平均测度表达式，如式（14）所示：

m_L＝k^2-a*m_kL    (14)

步骤六：根据步骤四中利用克里金法求出的各实测点的权重系数，推算出待插值点邻域内的多重分形克里金插值方程；

利用步骤四中求出的估值区域内各实测点的权重系数，求出整个估值区域测度表达式，并根据步骤五中得到的关系式，推算出待插值点在其小邻域内的多重分形克里金法插值方程，具体为：

根据利用克里金法求出的估值区域内各实测点的权重系数λ₁,λ₂,...,λ_P,可得尺度为kL的大正方形（估值区域）的平均测度m_kL，如式（15），

$m_{\mathrm{kL}}=\underset{i=1}{\overset{P}{\mathrm{\Σ}}}({\mathrm{\λ}}_i/w_{\mathrm{ai}})---\left(15\right)$

再根据式（14）可得式（15）得出地磁异常场在待插值点邻域内的多重分形克里金插值公式为：

$m_L(x_0,y_0)=k^{2-a}\underset{i=1}{\overset{P}{\mathrm{\Σ}}}({\mathrm{\λ}}_i*w_{\mathrm{ai}})---\left(16\right)$

即

$w(x_0,y_0)=k^{2-a}\underset{i=1}{\overset{P}{\mathrm{\Σ}}}({\mathrm{\λ}}_i*w_{\mathrm{ai}})---\left(17\right)$

其中，w(x₀,y₀)即为待插值点的地磁异常估计值。

步骤七：采用交叉验证法验证插值的精度并构建地磁图

设Xi为测区内一实测点，将此实测点暂时去除，利用步骤二至步骤六所述方法对此点进行插值，得到这一点的地磁异常估计值按此方法得到测区内所有实测点的地磁异常估计值，并采用均方根预测误差、平均估计误差百分比、相对均方差三个指标对插值结果进行评估。

（1）平均估计误差百分比PAEE；

$\mathrm{PAEE}=\frac{1}{m*\stackrel{\‾}{w}}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{[\hat{w}\left(X_i\right)-w\left(X_i\right)]}^2\×100\%$

式中：为位置X_i处的地磁异常随机变量的估计值，w(X_i)表示位置X_i处的地磁异常实测值，为测区内所有测点实测地磁要素平均值，m表示估值区域内实测点数。

（2）相对均方差RMSE；

$\mathrm{RMSE}=\frac{1}{m*s^2}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{[\hat{w}\left(X_i\right)-w\left(X_i\right)]}^2$

式中：s²为测区内所有测点地磁要素的样本方差。

（3）均方根预测误差RMSPE；

$\mathrm{RMSPE}=\sqrt{\frac{1}{n}\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{[\hat{w}\left(X_i\right)-w\left(X_i\right)]}^2}$

设定均方根预测误差、平均估计误差百分比、相对均方差三个指标的数值，如果插值结果符合指标要求，则根据插值后的地磁数据构建地磁基准图；否则，返回步骤二，重新进行插值。

Claims (3)

1.一种基于多重分形克里金法的局部地磁图构建方法，其特征在于，包括以下几个步骤：

步骤一：去除实测数据中的主磁场成分；

判断实测地磁数据中是否包含地球主磁场成分，如果包含，利用地磁场模型计算各测点主磁场强度，并在实测数据中将其去除；

所述的步骤一具体为：

(1)判断实测地磁数据中是否包含地球主磁场成分；

实测地磁数据为N个测点序列，测点序列为其中，1≤i≤N，分别为测点的地理经度、地理纬度、地磁测量值，g_i为矢量；

通过实测数据的大小，判断实测地磁数据中是否包含有主磁场成分，当实测地磁强度大于30000nT时，实测地磁数据中包含主磁场成分，当其位于1000nT—4000nT之间时，实测地磁数据中不包含主磁场成分；

(2)如果实测地磁数据中包含主磁场成分，利用地磁场模型计算各测点主磁场强度，并在实测数据中将其去除；

根据国际地磁参考场模型，利用截断阶数法获取测点在北向、东向、竖直向下方向上的地磁场分量X、Y、Z，X、Y、Z通过式(1)、式(2)、式(3)获取，具体为：

$X=\underset{n=1}{\overset{\∞}{\mathrm{\Σ}}}\underset{m=0}{\overset{n}{\mathrm{\Σ}}}{\left(\frac{a}{r}\right)}^{n+2}(g_n^m\cos \mathrm{m\λ}+h_n^m\sin \mathrm{m\λ})\frac{{\∂P}_n^m\left(\mathrm{\θ}\right)}{\∂\mathrm{\θ}}---\left(1\right)$

$Y=\underset{n=1}{\overset{\∞}{\mathrm{\Σ}}}\underset{m=0}{\overset{n}{\mathrm{\Σ}}}{\left(\frac{a}{r}\right)}^{n+2}(g_n^m\sin \mathrm{m\λ}-h_n^m\cos \mathrm{m\λ})\frac{{\mathrm{mP}}_n^m\left(\mathrm{\θ}\right)}{\sin \mathrm{\θ}}---\left(2\right)$

$Z=-\underset{n=1}{\overset{\∞}{\mathrm{\Σ}}}\underset{m=0}{\overset{n}{\mathrm{\Σ}}}(n+1){\left(\frac{a}{r}\right)}^{n+1}(g_n^m\cos \mathrm{m\λ}+h_n^m\sin \mathrm{m\λ})P_n^m\left(\mathrm{\θ}\right)---\left(3\right)$

其中，a为地球半径；r是地球表面或表面以上空间中的计算点与地心的距离；θ是地理余纬度，设为地理纬度，则有λ是地理经度；为高斯球谐系数；n、m分别为球谐系数的阶和次；( )为n阶m次缔合勒让德函数；

在公式(1)、(2)、(3)中将n＞13的高阶项截掉，得公式(4)、(5)、(6)，剩余低阶项部分能够准确的描述地球主磁场；

$X=\underset{n=1}{\overset{13}{\mathrm{\Σ}}}\underset{m=0}{\overset{n}{\mathrm{\Σ}}}{\left(\frac{a}{r}\right)}^{n+2}(g_n^m\cos \mathrm{m\λ}+h_n^m\sin \mathrm{m\λ})\frac{{\∂P}_n^m\left(\mathrm{\θ}\right)}{\∂\mathrm{\θ}}---\left(4\right)$

$Y=\underset{n=1}{\overset{13}{\mathrm{\Σ}}}\underset{m=0}{\overset{n}{\mathrm{\Σ}}}{\left(\frac{a}{r}\right)}^{n+2}(g_n^m\sin \mathrm{m\λ}-h_n^m\cos \mathrm{m\λ})\frac{{\mathrm{mP}}_n^m\left(\mathrm{\θ}\right)}{\sin \mathrm{\θ}}---\left(5\right)$

$Z=-\underset{n=1}{\overset{13}{\mathrm{\Σ}}}\underset{m=0}{\overset{n}{\mathrm{\Σ}}}(n+1){\left(\frac{a}{r}\right)}^{n+1}(g_n^m\cos \mathrm{m\λ}+h_n^m\sin \mathrm{m\λ})P_n^m\left(\mathrm{\θ}\right)---\left(6\right)$

将实测地磁数据中包含主磁场成分的测点依次带入式(4)、式(5)、式(6)中，分别得到各个测点的地球主磁场三分量X_Ri,Y_Ri,Z_Ri，即得此测点的地球主磁场矢量 $\stackrel{\→}{g_{\mathrm{Ri}}}=(X_{\mathrm{Ri}},Y_{\mathrm{Ri}},Z_{\mathrm{Ri}});$

将g_i与作矢量差，得出第i个测点的地磁异常矢量g_Ui,相应的地磁异常强度w_i＝|g_Ui|；

步骤二：选定估值区域；

根据实测地磁数据的分辨率以及地磁异常场的空间特性，以待插值点为中心，选定估值区域；

步骤三：选定高斯模型描述地磁异常场的变异特性，拟合高斯模型中的未知参数，并利用此模型计算估值区域内各实测点之间以及各实测点和待插值点之间的半方差；

在估值区域内，利用实测地磁数据，计算不同的分离距离对应的变异函数值，从而得到分离距离与变异函数值序列对，根据此序列对，利用拟合法推算出高斯模型中的未知参数，并利用此模型计算估值区域内各实测点之间以及实测点和待插值点之间的半方差；

所述的步骤三具体为：

(1)选定高斯模型描述地磁异常场的变异特性，如(7)式：

$\mathrm{\γ}\left(h\right)=\left\{\begin{array}{c}0,h=0\\ C_0+C_1(1-e^{-{(h/a)}^2}),h>0\end{array}\right.---\left(7\right)$

其中：h为样本间的距离,变程为C₀为块金值，C₀+C₁为基台值，即h为无穷大时的变异函数值，C₀、C₁、a均为未知量；

(2)利用实测地磁数据，计算不同的分离距离对应的变异函数值，从而得到分离距离与变异函数值序列对，根据此序列对，利用拟合法推算出高斯模型中的未知参数；

根据实测地磁数据，获取估值区域内实测地磁异常数据序列为w_a1、w_a1、...、w_aP，P为估值区域内实测点的个数，根据式(8)得到估值区域内的变异函数值：

${\mathrm{\γ}}^{*}\left(h^{\′}\right)=\frac{1}{{2N}_{n^{\′}}}\underset{j=1}{\overset{N_{h^{\′}}}{\mathrm{\Σ}}}{[w_j^{*}-w_j]}^2---\left(8\right)$

其中，h′为分离距离，N_h′是估值区域内距离小于h′实测点对数，和w_j为距离小于h′的一对实测数据；

(3)利用拟合法推算出高斯模型中的未知参数，并利用此模型计算估值区域内各实测点之间以及实测点和待插值点之间的半方差；

通过式(8)，根据估值区域内各实测点的之间的距离，通过调整h′的值，得到变异函数值与分离距离序列对，序列对同时满足式(7)，由此拟合出高斯模型中的三个未知参数C₀、C₁、a，并可求得变程

利用已知的高斯模型计算出估值区域内各实测点之间的半方差γ_jk以及各实测点与待插值点之间的半方差γ_0i，其中1≤j,k≤P；

步骤四：计算估值区域内各测点的权重系数；

利用步骤三中得出的半方差值，通过解空间相关系数矩阵方程求出各实测点的权重系数；

所述的步骤四具体为：

设估值区域内各实测点的权重系数为λ₁,λ₂,...,λ_P，空间相关系数矩阵方程如下：

$\left[\begin{array}{cccc}{\mathrm{\γ}}_{11}& ...& {\mathrm{\γ}}_{1P}& 1\\ .& & .& .\\ .& & .& .\\ .& & .& .\\ {\mathrm{\γ}}_{P1}& ...& {\mathrm{\γ}}_{\mathrm{PP}}& 1\\ 1& ...& 1& 0\end{array}\right]\left[\begin{array}{c}{\mathrm{\λ}}_1\\ .\\ .\\ .\\ {\mathrm{\λ}}_P\\ \mathrm{\μ}\end{array}\right]=\left[\begin{array}{c}{\mathrm{\γ}}_{01}\\ .\\ .\\ .\\ {\mathrm{\γ}}_{0P}\\ 1\end{array}\right]---\left(9\right)$

其中：μ为拉格朗日乘数，γ_jk为步骤三中(3)求出的各实测点之间的半方差，γ_0j为步骤三中(3)求出的实测点与待插值点之间的半方差，1≤j,k≤P，通过解此方程(9)得到各测点的权重系数；

步骤五：利用多重分形理论推算待插值点邻域内的测度表达式；

在估值区域内，利用测度与尺度的在双对数坐标系中的线性关系，用最小二乘拟合法得出插值点邻域内的分形维数，进而推算出以待插值点为中心的小正方形的测度与整个估值区域测度的关系式；

所述的步骤五具体为：

根据分维定义，设F为估值区域内实测数据集，M(F)为尺度l下F的测度，当l→0时，有：

M(F)≌cl^a    (10)

其中，c为常数，a为数据集的分维；测度M(F)与尺度l服从指数定律，在双对数坐标系中表现为一条直线，即

lnM(F)≌lnc+alnl    (11)

采用最小二乘方法拟合式(11)，通过双对数图中直线的斜率求出a；

在估值区域内，以待插值点为中心，建立尺度为L的小正方形，并令D＝kL，D估值区域的边长，k＞1，则该邻域的测度M_L和估值区域的测度M_kL分别如(12)式和(13)式所示：

M_L＝m_L*L²＝c*L^a    (12)

M_kL＝m_kL*(kL)²＝c*(kL)^a    (13)

式中，m_L和m_kL分别为小正方形和大正方形的平均测度；

综合式(12)和式(13)得到待插值点邻域内的平均测度表达式，如式(14)所示：

m_L＝k^2-a*m_kL    (14)

步骤六：根据步骤四中求出的各实测点的权重系数，推算出待插值点邻域内的多重分形克里金插值方程；

利用步骤四中求出的估值区域内各实测点的权重系数，求出整个估值区域测度表达式，并根据步骤五中得到的关系式，推算出待插值点在其小邻域内的多重分形克里金法插值方程；

所述的步骤六具体为：

根据利用克里金法求出的估值区域内各实测点的权重系数λ₁,λ₂,...,λ_P,得到尺度为kL的大正方形的平均测度m_kL，如式(15)，

$m_{\mathrm{kL}}=\underset{i=1}{\overset{P}{\mathrm{\Σ}}}({\mathrm{\λ}}_i*w_{\mathrm{ai}})---\left(15\right)$

再根据式(14)可得式(15)得出地磁异常场在待插值点邻域内的多重分形克里金插值公式为：

$m_L(x_0,y_0)=k^{2-a}\underset{i=1}{\overset{P}{\mathrm{\Σ}}}({\mathrm{\λ}}_i*w_{\mathrm{ai}})---\left(16\right)$

即

$w(x_0,y_0)=k^{2-a}\underset{i=1}{\overset{P}{\mathrm{\Σ}}}({\mathrm{\λ}}_i*w_{\mathrm{ai}})---\left(17\right)$

其中，w(x₀,y₀)即为待插值点的地磁异常估计值；

步骤七：采用交叉验证法验证插值的精度并构建地磁图；

设X_i为测区内一实测点，将此实测点暂时去除，利用步骤二至步骤六所述方法对此点进行插值，得到这一点的地磁异常估计值按步骤二至步骤六得到测区内所有实测点的地磁异常估计值，并采用均方根预测误差、平均估计误差百分比、相对均方差三个指标对插值结果进行评估；

设定均方根预测误差、平均估计误差百分比、相对均方差三个指标的数值，如果插值结果符合指标要求，则根据插值后的地磁数据构建地磁基准图；否则，返回步骤二，重新进行插值；

所述的步骤七中均方根预测误差、平均估计误差百分比、相对均方差具体为：

(1)平均估计误差百分比PAEE；

$\mathrm{PAEE}=\frac{1}{m*\stackrel{\‾}{w}}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{[\hat{w}\left(X_i\right)-w\left(X_i\right)]}^2\×100\%$

式中：为位置X_i处的地磁异常随机变量的估计值，w(X_i)表示位置X_i处的地磁异常实测值，为测区内所有测点实测地磁要素平均值，m表示估值区域内实测点数；

(2)相对均方差RMSE；

$\mathrm{RMSE}=\frac{1}{m*s^2}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{[\hat{w}\left(X_i\right)-w\left(X_i\right)]}^2$

式中：s²为测区内所有测点地磁要素的样本方差；

(3)均方根预测误差RMSPE；

$\mathrm{RMSPE}=\sqrt{\frac{1}{m}\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{[\hat{w}\left(X_i\right)-w\left(X_i\right)]}^2}.$

2.根据权利要求1所述的一种基于多重分形克里金法的局部地磁图构建方法，其特征在于，所述的步骤二具体为：

在测区内，以待插值点为中心，建立以D为尺度的正方形区域作为估值区域；

设测区内实测点的平均间距为d，取D＝10d～15d。

3.根据权利要求2所述的一种基于多重分形克里金法的局部地磁图构建方法，其特征在于，当航空磁测时，在经、纬度方向上4″≤d≤6″。