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

摘要

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

主权项

一种基于多重分形克里金法的局部地磁图构建方法，其特征在于，包括以下几个步骤：步骤一：去除实测数据中的主磁场成分；判断实测地磁数据中是否包含地球主磁场成分，如果包含，利用地磁场模型计算各测点主磁场强度，并在实测数据中将其去除；所述的步骤一具体为：(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{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{\left(\frac{a}{r}\right)}^{n+2}(g_n^m\cos \mathrm{m\lambda }+h_n^m\sin \mathrm{m\lambda })\frac{{\partial P}_n^m\left(\theta \right)}{\partial \theta }---\left(1\right)$$Y=\underset{n=1}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{\left(\frac{a}{r}\right)}^{n+2}(g_n^m\sin \mathrm{m\lambda }-h_n^m\cos \mathrm{m\lambda })\frac{{\mathrm{mP}}_n^m\left(\theta \right)}{\sin \theta }---\left(2\right)$$Z=-\underset{n=1}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}(n+1){\left(\frac{a}{r}\right)}^{n+1}(g_n^m\cos \mathrm{m\lambda }+h_n^m\sin \mathrm{m\lambda })P_n^m\left(\theta \right)---\left(3\right)$其中，a为地球半径；r是地球表面或表面以上空间中的计算点与地心的距离；θ是地理余纬度，设为地理纬度，则有λ是地理经度；为高斯球谐系数；n、m分别为球谐系数的阶和次；( )为n阶m次缔合勒让德函数；在公式(1)、(2)、(3)中将n＞13的高阶项截掉，得公式(4)、(5)、(6)，剩余低阶项部分能够准确的描述地球主磁场；$X=\underset{n=1}{\overset{13}{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{\left(\frac{a}{r}\right)}^{n+2}(g_n^m\cos \mathrm{m\lambda }+h_n^m\sin \mathrm{m\lambda })\frac{{\partial P}_n^m\left(\theta \right)}{\partial \theta }---\left(4\right)$$Y=\underset{n=1}{\overset{13}{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{\left(\frac{a}{r}\right)}^{n+2}(g_n^m\sin \mathrm{m\lambda }-h_n^m\cos \mathrm{m\lambda })\frac{{\mathrm{mP}}_n^m\left(\theta \right)}{\sin \theta }---\left(5\right)$$Z=-\underset{n=1}{\overset{13}{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}(n+1){\left(\frac{a}{r}\right)}^{n+1}(g_n^m\cos \mathrm{m\lambda }+h_n^m\sin \mathrm{m\lambda })P_n^m\left(\theta \right)---\left(6\right)$将实测地磁数据中包含主磁场成分的测点依次带入式(4)、式(5)、式(6)中，分别得到各个测点的地球主磁场三分量X_Ri,Y_Ri,Z_Ri，即得此测点的地球主磁场矢量$\overrightarrow{g_{\mathrm{Ri}}}=(X_{\mathrm{Ri}},Y_{\mathrm{Ri}},Z_{\mathrm{Ri}});$将g_i与作矢量差，得出第i个测点的地磁异常矢量g_Ui,相应的地磁异常强度w_i＝|g_Ui|；步骤二：选定估值区域；根据实测地磁数据的分辨率以及地磁异常场的空间特性，以待插值点为中心，选定估值区域；步骤三：选定高斯模型描述地磁异常场的变异特性，拟合高斯模型中的未知参数，并利用此模型计算估值区域内各实测点之间以及各实测点和待插值点之间的半方差；在估值区域内，利用实测地磁数据，计算不同的分离距离对应的变异函数值，从而得到分离距离与变异函数值序列对，根据此序列对，利用拟合法推算出高斯模型中的未知参数，并利用此模型计算估值区域内各实测点之间以及实测点和待插值点之间的半方差；所述的步骤三具体为：(1)选定高斯模型描述地磁异常场的变异特性，如(7)式：$\gamma \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)得到估值区域内的变异函数值：${\gamma }^{*}\left(h^{\prime }\right)=\frac{1}{{2N}_{n^{\prime }}}\underset{j=1}{\overset{N_{h^{\prime }}}{\Sigma }}{[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}{\gamma }_{11}& ...& {\gamma }_{1P}& 1\\ .& & .& .\\ .& & .& .\\ .& & .& .\\ {\gamma }_{P1}& ...& {\gamma }_{\mathrm{PP}}& 1\\ 1& ...& 1& 0\end{array}\right]\left[\begin{array}{c}{\lambda }_1\\ .\\ .\\ .\\ {\lambda }_P\\ \mu \end{array}\right]=\left[\begin{array}{c}{\gamma }_{01}\\ .\\ .\\ .\\ {\gamma }_{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}{\Sigma }}({\lambda }_i*w_{\mathrm{ai}})---\left(15\right)$再根据式(14)可得式(15)得出地磁异常场在待插值点邻域内的多重分形克里金插值公式为：$m_L(x_0,y_0)=k^{2-a}\underset{i=1}{\overset{P}{\Sigma }}({\lambda }_i*w_{\mathrm{ai}})---\left(16\right)$即$w(x_0,y_0)=k^{2-a}\underset{i=1}{\overset{P}{\Sigma }}({\lambda }_i*w_{\mathrm{ai}})---\left(17\right)$其中，w(x₀,y₀)即为待插值点的地磁异常估计值；步骤七：采用交叉验证法验证插值的精度并构建地磁图；设X_i为测区内一实测点，将此实测点暂时去除，利用步骤二至步骤六所述方法对此点进行插值，得到这一点的地磁异常估计值按步骤二至步骤六得到测区内所有实测点的地磁异常估计值，并采用均方根预测误差、平均估计误差百分比、相对均方差三个指标对插值结果进行评估；设定均方根预测误差、平均估计误差百分比、相对均方差三个指标的数值，如果插值结果符合指标要求，则根据插值后的地磁数据构建地磁基准图；否则，返回步骤二，重新进行插值；所述的步骤七中均方根预测误差、平均估计误差百分比、相对均方差具体为：(1)平均估计误差百分比PAEE；$\mathrm{PAEE}=\frac{1}{m*\bar{w}}\underset{i=1}{\overset{m}{\Sigma }}{[\hat{w}\left(X_i\right)-w\left(X_i\right)]}^2\times 100\%$式中：为位置X_i处的地磁异常随机变量的估计值，w(X_i)表示位置X_i处的地磁异常实测值，为测区内所有测点实测地磁要素平均值，m表示估值区域内实测点数；(2)相对均方差RMSE；$\mathrm{RMSE}=\frac{1}{m*s^2}\underset{i=1}{\overset{m}{\Sigma }}{[\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}{\Sigma }}{[\hat{w}\left(X_i\right)-w\left(X_i\right)]}^2}.$