专利名称：一种基于多重分形克里金法的局部地磁图构建方法
技术领域：
本发明属于地磁场数据处理领域，具体涉及一种基于多重分形克里金法的局部地磁图构建方法。
背景技术：
高精度的地磁基准图的获取是实现精确地磁导航的关键。地磁异常场是地磁场的重要组成部分，由地壳表面分布的矿产、岩石、人造磁场等产生，分布很稳定，几乎不随时间变化，通常只在近地空间存在；它的空间结构复杂，波形尺度小，具有局部、随机、异常的特征，非常适合作为地磁导航的参照。由于我国地域广阔，开展广泛的陆地及海洋磁测工作需要耗费巨 大的成本，因此依靠直接测量获取细致的地磁异常数据难度很大。在已有数据的基础上通过插值来获得高精度、高分辨率的地磁异常网格数据是对实测工作的一项重要补充。陈丽等研究了双线性插值方法在构建高精度地磁图中的应用，实验证明该方法计算量小，实时性较好，但它只适合应用在与位置关系密切且在大范围内具有明显变化规律的地球主磁场的处理中，而对于局部、随机特性比较明显的地磁异常场，该方法并不适用；王金玲等讨论了反距离加权法在空间数据插值中的应用，它根据搜索范围内实测点和待插值点的距离作为加权依据，忽视了实测点之间的位置关系，如果实测点的位置相对待插值点具有明显的重复性，那么插值结果将出现较大的误差；张晓明等单纯使用克里金法进行了局部地磁图的构建，并通过实验证明该方法能够比较准确的反映局部地区地磁信息，但是这种方法在充分描述地磁数据空间相关性的同时，没有考虑地磁异常数据的小尺度奇异性，即忽略了地磁异常场的高频信息。克里金法是一种用于空间插值的地理统计方法，充分吸收了地理统计思想，认为任何在空间连续变化的属性是非常不规则的，不能简单用平滑的数学函数进行模拟，而采用随机表面可以恰当的描述，它是一种无偏估计，是空间相关意义下的最优(最小方差)估计。克里金法假设某种属性的空间变化既不是完全随机也不是完全确定的，可以表示为3种主要成分的和与恒定均值和趋势有关的结构成分；与空间变化有关的随机变量，即区域性变量；与空间无关的随机噪声。地磁场中的地球主磁场、地磁异常场及扰动磁场3个主要成分与克里金法中假定的结构性成分、随机变量以及随机噪声具有一一对应关系。因此克里金法较为准确的反应地磁异常场的空间相关特性规律，可以用于地磁异常的空间插值。但是克里金法建立在空间相关系数矩阵极小意义下的求解方式，决定了它是一个空间滑动平均或低通滤波过程，高频、局部与弱信号在相关系数矩阵的方差中占得比例很小，难以反映地磁异常场在小尺度上的奇异特征。而多重分形理论可以弥补这个缺点，它是由Mandelbrot在研究湍流时提出来的，专门应用于研究物理量和其他量在几何支撑上的奇异性分布，它是定义在分形上，由多个标度指数的奇异测度所组成的集合。Spector等人通过对航磁数据谱的统计分析，表明磁异常具有分形特征。

发明内容
本发明的目的是为了解决上述问题，提出一种基于多重分形克里金法的局部地磁图构建方法，本发明在已有较稀疏的实测地磁异常数据的基础上，针对地壳表面空间结构复杂且非常稳定的地磁异常场，采用多重分形克里金法对实测地磁数据进行插值，建立高精度、高分辨率的局部地磁异常网格数据，进而为实现地磁导航构建局部地磁图。本发明不仅能充分的表征磁异常的区域性变化，同时能保留更多的高频信息。一种基于多重分形克里金法的局部地磁图构建方法，包括以下几个步骤步骤一、去除实测数据中的主磁场成分；判断实测地磁数据中是否包含地球主磁场成分，如果包含，利用地磁场模型计算各测点主磁场强度,并在实测数据中将其去除；
步骤二、选定估值区域；根据实测地磁数据的分辨率以及地磁异常场的空间特性，以待插值点为中心，选定估值区域；步骤三选定高斯模型描述地磁异常场的变异特性，拟合高斯模型中的未知参数，并利用此模型计算估值区域内各控制点间以及各控制点与待插值点之间的半方差；在估值区域内，利用实测地磁数据，计算不同的分离距离对应的变异函数值，从而得到分离距离与变异函数值序列对，根据此序列对，利用拟合法推算出高斯模型中的未知参数，并利用此模型计算估值区域内各实测点之间以及实测点和待插值点之间的半方差；步骤四、计算估值区域内各测点的权重系数；利用步骤三中得出的半方差值，通过解空间相关系数矩阵方程求出各实测点的权重系数；步骤五、利用多重分形理论推算待插值点邻域内的测度表达式；在估值区域内，利用测度与尺度的在双对数坐标系中的线性关系，用最小二乘拟合法得出插值点邻域内的分形维数，进而推算出以带插值点为中心的小正方形(邻域)的测度与整个估值区域测度的关系式；步骤六、建立待插值点在其小邻域内的多重分形克里金法插值方程；利用步骤四中求出的估值区域内各实测点的权重系数，求出整个估值区域测度表达式，并根据步骤五中得到的关系式，推算出待插值点在其小邻域内的多重分形克里金法插值方程；步骤七、采用交叉验证法验证插值的精度并构建地磁图；在测区内每一个有观测值的地方，将此观测值暂时去除，用以此点为中心建立的估值区域内的其他观测值以及本专利发明的插值方法，估计这一观察点的值，然后将暂时去除的观察值放回，重复以上步骤，对测区所有的观察点进行估值。并采用均方根预测误差、平均估计误差百分比、相对均方差三个指标对插值结果进行评估。设定评价指标，如果插值结果符合指标要求，则根据插值后的地磁数据构建地磁基准图；否则，返回步骤二，重新进行插值。本发明的优点在于本发明的优点在于在插值过程中，采用多重分形理论精确的反映出地磁异常场在待插值点邻域内的奇异特征，并利用克里金方法充分考虑了地磁异常场的空间相关特性，能够在已有实测数据的基础上，通过插值得到高精度高、分辨率的地磁异常网格数据。


图1是本发明的方法流程具体实施例方式下面将结合附图和实施例对本发明作进一步的详细说明。本发明是一种基于多重分形克里金法的局部地磁图构建方法，流程如图I所示，包括以下几个步骤步骤一去除实测数据中的主磁场成分；判断实测地磁数据中是否包含地球主磁场成分，如果包含，利用地磁场模型计算各测点主磁场强度，并在实测数据中将其去除，具体为(1)判断实测地磁数据中是否包含地球主磁场成分；实测地磁数据为N个测点序列，测点序列为(λ φ gi),其中，I彡i ^N, Ai, ^gi分别为测点的地理经度、地理纬度、地磁测量值，gi为矢量；因为主磁场在地表处的强度为30000— 70000nT，占地磁场总量的95%以上，地磁异常场强度占地磁场强度的4 % — 5%，在1000ηΤ-4000ηΤ之间，所以，通过观察实测数据的大小，判断实测地磁数据中是否包含有主磁场成分，当实测地磁强度大于30000ηΤ时，实测地磁数据中包含主磁场成分，当其位于1000ηΤ-4000ηΤ之间时，实测地磁数据中不包含主磁场成分。(2)如果实测地磁数据中包含主磁场成分，利用地磁场模型计算各测点主磁场强度，并在实测数据中将其去除；根据国际地磁参考场模型(IGRF)，利用截断阶数法获取测点在北向、东向、竖直(向下)方向上的地磁场分量乂、￥、2，乂、￥、2通过式(1)、式(2)、式(3)获取，具体为
权利要求
1.一种基于多重分形克里金法的局部地磁图构建方法，其特征在于，包括以下几个步骤 步骤ー去除实测数据中的主磁场成分； 判断实测地磁数据中是否包含地球主磁场成分，如果包含，利用地磁场模型计算各测点主磁场强度，并在实测数据中将其去除； 步骤ニ 选定估值区域； 根据实测地磁数据的分辨率以及地磁异常场的空间特性，以待插值点为中心，选定估值区域； 步骤三选定高斯模型描述地磁异常场的变异特性，拟合高斯模型中的未知參数，井利用此模型计算估值区域内各实测点之间以及各实测点和待插值点之间的半方差； 在估值区域内，利用实测地磁数据，计算不同的分离距离对应的变异函数值，从而得到分离距离与变异函数值序列对，根据此序列对，利用拟合法推算出高斯模型中的未知參数，并利用此模型计算估值区域内各实测点之间以及实测点和待插值点之间的半方差； 步骤四计算估值区域内各测点的权重系数； 利用步骤三中得出的半方差值，通过解空间相关系数矩阵方程求出各实测点的权重系数； 步骤五利用多重分形理论推算待插值点邻域内的測度表达式； 在估值区域内，利用測度与尺度的在双对数坐标系中的线性关系，用最小ニ乘拟合法得出插值点邻域内的分形维数，进而推算出以待插值点为中心的小正方形的測度与整个估值区域測度的关系式； 步骤六根据步骤四中利用克里金法求出的各实测点的权重系数，推算出待插值点邻域内的多重分形克里金插值方程； 利用步骤四中求出的估值区域内各实测点的权重系数，求出整个估值区域測度表达式，井根据步骤五中得到的关系式，推算出待插值点在其小邻域内的多重分形克里金法插值方程； 步骤七采用交叉验证法验证插值的精度并构建地磁图； 设Xi为测区内ー实测点，将此实测点暂时去除，利用步骤ニ至步骤六所述方法对此点进行插值，得到这一点的地磁异常估计值ル(孓)，按步骤ニ至步骤六得到测区内所有实测点的地磁异常估计值，并采用均方根预测误差、平均估计误差百分比、相对均方差三个指标对插值结果进行评估； 设定均方根预测误差、平均估计误差百分比、相对均方差三个指标的数值，如果插值结果符合指标要求，则根据插值后的地磁数据构建地磁基准图；否则，返回步骤ニ，重新进行插值。
2.根据权利要求I所述的ー种基于多重分形克里金法的局部地磁图构建方法，其特征在于，所述的步骤一具体为 (I)判断实测地磁数据中是否包含地球主磁场成分； 实测地磁数据为N个测点序列，测点序列为(Xi, φ{, gi),其中，I ≤ i≤ N, λ i； φi, gi分别为测点的地理经度、地理纬度、地磁测量值，gi为矢量； 通过实测数据的大小，判断实测地磁数据中是否包含有主磁场成分，当实测地磁强度大于30000nT时，实测地磁数据中包含主磁场成分，当其位于IOOOnT—4000ηΤ之间时，实测地磁数据中不包含主磁场成分； (2)如果实测地磁数据中包含主磁场成分，利用地磁场模型计算各测点主磁场强度，并在实测数据中将其去除； 根据国际地磁參考场模型，利用截断阶数法获取测点在北向、东向、竖直向下方向上的 地磁场分量乂3、2，乂、￥、2通过式(1)、式(2)、式(3)获取，具体为
3.根据权利要求I所述的ー种基于多重分形克里金法的局部地磁图构建方法，其特征在于，所述的步骤ニ具体为 在测区内，以待插值点为中心，建立以D为尺度的正方形区域作为估值区域； 设测区内实测点的平均间距为山取D = 10cTl5d。
4.根据权利要求3所述的ー种基于多重分形克里金法的局部地磁图构建方法，其特征在于，当航空磁测时，在经、纬度方向上4"彡d<6"。
5.根据权利要求I所述的ー种基于多重分形克里金法的局部地磁图构建方法，其特征在于，所述的步骤三具体为 (I)选定高斯模型描述地磁异常场的变异特性，如(7)式
6.根据权利要求I所述的ー种基于多重分形克里金法的局部地磁图构建方法，其特征在于，所述的步骤四具体为设估值区域内各实测点的权重系数为A1, λ 2，...，λ ρ，空间相关系数矩阵方程如下
7.根据权利要求I所述的ー种基于多重分形克里金法的局部地磁图构建方法，其特征在于，所述的步骤五具体为根据分维定义，设F为估值区域内实测数据集，M(F)为尺度I下F的測度，当I — O吋，有
8.根据权利要求I所述的ー种基于多重分形克里金法的局部地磁图构建方法，其特征在于，所述的步骤六具体为 根据利用克里金法求出的估值区域内各实测点的权重系数λ ” λ 2，. . .，λ p，得到尺度为kL的大正方形的平均测度1 ，如式(15)，
9.根据权利要求I所述的ー种基于多重分形克里金法的局部地磁图构建方法，其特征在于，所述的步骤七中均方根预测误差、平均估计误差百分比、相对均方差具体为 (1)平均估计误差百分比PAEE；
全文摘要
本发明公开了一种基于多重分形克里金法的局部地磁图构建方法，包括步骤一去除实测数据中的主磁场成分；步骤二选定估值区域；步骤三拟合高斯模型，计算估值区域内各控制点间以及各控制点与待插值点之间的半方差；步骤四计算估值区域内各测点的权重系数；步骤五利用多重分形理论推算待插值点邻域内的测度表达式；步骤六建立待插值点在其小邻域内的多重分形克里金法插值方程；步骤七、采用交叉验证法验证插值的精度并构建地磁图。本发明充分考虑地磁异常场的空间相关性及其在待插值点邻域内的奇异特征，根据实测数据进行精确的插值，使得插值结果更加符合实际地磁异常场的空间变化趋势，为地磁导航构建高精度、高分辨率的局部地磁图。
文档编号G01V3/40GK102841385SQ201210237318
公开日2012年12月26日 申请日期2012年7月10日 优先权日2012年7月10日
发明者赵玉新, 常帅, 刘厂, 张振兴, 沈志峰 申请人:哈尔滨工程大学