CN104808249A - 地磁场高精度建模方法 - Google Patents

Abstract

本发明涉及地磁场建模领域，公开了一种地磁场高精度建模方法，通过获取测试点的地磁场异常值，选取研究区域周边预设范围的测试点参与地磁场模型的构建，确定矩形区域，将测试点的地理坐标转换为矩形区域内的矩谐坐标，将地理坐标系下所述测试点的地磁场分量值转换为所述矩形区域内的矩谐坐标系下的地磁场分量值，建立地磁场模型，对地磁场模型的精度进行验证，当精度符合要求时，地磁场模型构建完成。有效减小了地磁场模型的边界效应，使用奇异值分解对地磁场模型的系数矩阵进行求解，减小了运算量，使得地磁场模型的构建更加精确和快速。

Description

地磁场高精度建模方法

技术领域

本发明涉及地磁场建模领域，尤其涉及一种地磁场高精度建模方法。

背景技术

地磁信息在航空航天、环境监测、地球科学、国防建设、资源探测、地震预报等诸多领域具有十分重要的实际应用价值，成为了各国宝贵的战略性资源，但地磁应用依赖于高精度、高分辨率的地磁场建模，因此，地磁场建模是决定地磁应用广度和深度的关键性因素。

地磁建模是通过对已有的离散的地磁测量点进行插值计算，进而得到区域三维空间的磁场结构的过程。通常的建模方法如多项式、样条曲面和多面函数等数值方法缺乏对位场特征的约束，这些方法的误差估计是模型对观测数据拟合近似的程度，往往会忽略其在物理上的合理性，其对地磁场单要素的拟合误差比较小，但各地磁分量之间不能体现磁场的物理意义。而现有技术中的矩谐分析类似于球谐分析，是从满足位场的拉普拉斯方程出发，通过构建研究区域大小的矩形平面来近似球面，在矩形区域内进行协和分析，求解拉普拉斯方程，并写成相应的级数形式的一种方法，其不同阶的系数反映了不同波长的空间磁场分布情况，进而得到区域空间的全波长磁场分布。现有技术中的矩谐分析方法存在很明显的边界效应，在研究区域的边界产生龙格现象，进而影响地磁建模的精度。

发明内容

本发明提供一种地磁场高精度建模方法，解决现有技术中的地磁场建模存在明显边界效应的技术问题。

本发明的目的是通过以下技术方案实现的：

一种地磁场高精度建模方法，包括：

获取测试点的地磁场测试值，并计算所述测试点的地磁场值，将所述地磁场测试值与所述地磁场值做差，以获取所述测试点的地磁场异常值；

根据所述测试点的分布情况和有限元法中的三角剖分方法，选取研究区域周边预设范围的测试点参与地磁场模型的构建，计算研究区域周边预设范围的测试点的地磁场值，确定矩形区域的大小、矩形坐标的原点以及地磁场模型的截止阶数；

将所述测试点的地理坐标转换为所述矩形区域内的矩谐坐标；

将地理坐标系下所述测试点的地磁场分量值转换为所述矩形区域内的矩谐坐标系下的地磁场分量值；

选取部分测试点作为地磁场模型的精度检查点，根据剩余测试点在所述矩形区域内的矩谐坐标以及矩谐坐标系下的地磁场分量值，构建地磁场模型，并使用奇异值分解算法对模型系数矩阵进行求解计算；

通过地磁场模型的精度检查点对建立的地磁场模型的精度进行检查，判断地磁场模型的精度是否符合要求，当精度符合要求时，地磁场模型构建完成。

通过本发明提供的一种地磁场高精度建模方法，获取测试点的地磁场异常值，选取研究区域周边预设范围的测试点参与地磁场模型的构建，确定矩形区域，将测试点的地理坐标转换为矩形区域内的矩谐坐标，将地理坐标系下所述测试点的地磁场分量值转换为所述矩形区域内的矩谐坐标系下的地磁场分量值，建立地磁场模型，对地磁场模型的精度进行验证，当精度符合要求时，地磁场模型构建完成。有效减小了地磁场模型的边界效应，使用奇异值分解对地磁场模型的系数矩阵进行求解，减小了运算量，使得地磁场模型的构建更加精确和快速。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可根据这些附图获得其他的附图。

图1为本发明实施例提供的一种地磁场高精度建模方法的流程图。

具体实施方式

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

如图1为本发明实施例提供的一种地磁场高精度建模方法的流程图，包括：

步骤101、获取测试点的地磁场异常值；

其中，主地磁场在地表处约为3万纳特-7万纳特，占地磁场的95％左右，地磁异常场只占地磁场的5％左右，因此，需要获取测试点的地磁场测试值，并计算所述测试点的地磁场值，将所述地磁场测试值与所述地磁场值做差，以获取所述测试点的地磁场异常值。

利用IGRF(国际地磁参考场)模型，来获得测点的主地磁场值，并取n<14的低阶项来描述地球的主磁场，计算公式为

$V=a\underset{n=1}{\overset{13}{\mathrm{\&Sigma;}}}\underset{m=0}{\overset{n}{\mathrm{\&Sigma;}}}{\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(\cos \mathrm{\&theta;}\right)---\left(1\right)$

$X=\underset{n=1}{\overset{13}{\mathrm{\&Sigma;}}}\underset{m=0}{\overset{n}{\mathrm{\&Sigma;}}}{\left[\frac{a}{r}\right]}^{n+2}(g_n^m\cos \mathrm{m\&lambda;}+h_n^m\sin \mathrm{m\&lambda;})\frac{d}{\mathrm{d\&theta;}}{P}_n^m\left(\cos \mathrm{\&theta;}\right)---\left(2\right)$

$Y=\underset{n=1}{\overset{13}{\mathrm{\&Sigma;}}}\underset{m=0}{\overset{n}{\mathrm{\&Sigma;}}}{\left[\frac{a}{r}\right]}^{n+2}\frac{m}{\sin \mathrm{\&theta;}}(g_n^m\cos \mathrm{m\&lambda;}-h_n^m\sin \mathrm{m\&lambda;})P_n^m\left(\cos \mathrm{\&theta;}\right)---\left(3\right)$

$Z=\underset{n=1}{\overset{13}{\mathrm{\&Sigma;}}}\underset{m=0}{\overset{n}{\mathrm{\&Sigma;}}}-(n+1){\left[\frac{a}{r}\right]}^{n+2}(g_n^m\cos \mathrm{m\&lambda;}+h_n^m\sin \mathrm{m\&lambda;})P_n^m\left(\cos \mathrm{\&theta;}\right)---\left(4\right)$

其中，V为地磁场总场强；X,Y,Z分别为地磁场的三分量值；a为地球半径；θ是地理余纬；λ是地理经度；为高斯球谐系数(该系数可利用均匀分布于全球的足够多的测点值通过最小二乘法拟合求得，也可以通过查询官方公布的高斯系数表获得)；n,m分别为球谐系数的阶和次，这里取n＝12,m＝10；为n阶m次缔合勒让德函数；r是测点到地心的距离，其具体计算方法如下：

其中为地理纬度；λ为地理经度；N为卯酉圈曲率半径；a为椭球长半径，a＝6378137m；e为椭球第一偏心率，e＝0.00669437999013。

将测试点的地理坐标代入公式(1)-(4)，分别得到该测试点的三分量场值x1，y1，z1，将其与实际的地磁场测试值的三分量场值x0，y0，z0做差Δx＝x₀-x₁，Δy＝y₀-y₁，Δz＝z₀-z₁，以获取所述测试点的地磁场异常值Δx、Δy、Δz。

步骤102、选取研究区域周边预设范围的测试点参与地磁场模型的构建，确定矩形区域及模型参数；

其中，步骤102具体可以包括：

步骤102-1、根据研究区域整体面积及平均点距，确定研究区域周边的所述预设范围；

步骤102-2、根据有限元法中的三角剖分方法，选取所述预设范围内的测试点参与地磁场模型的构建；

其中，选取研究区域周边预设范围的测试点的原则为1、距离研究区域越近剖分的越精细，测试点数越密集；距离研究区越远剖分的越粗糙，测试点数越稀疏。2、对地磁数据梯度变化较大的区域进行加密，增加测试点数。

步骤102-3、计算研究区域周边预设范围的测试点的地磁场值，确定矩形区域的大小、矩形坐标的原点以及地磁场模型的截止阶数。

其中，根据所述测试点的分布情况和有限元法中的三角剖分方法，选取研究区域周边预设范围的测试点参与地磁场模型的构建，以减小边界效应的影响，并通过公式(1)至(4)计算研究区域周边预设范围的测试点的地磁场值。矩形区域的大小以能完全覆盖测试区所有点位为准，南北边长Lx，东西边长Ly，一般取矩形的中心为矩形坐标的原点，视数据质量确定地磁场模型的截止阶数，数据质量好选取较高的阶数，反之选取较低的阶数。

步骤103、将测试点的地理坐标转换为矩形区域内的矩谐坐标；

其中，地理坐标包括余纬θ和经度λ，矩谐坐标

θ₀,λ₀分别为矩形坐标的原点的地理坐标系的余纬和经度，R_e为地球半径。

步骤104、将地理坐标系下所述测试点的地磁场分量值转换为所述矩形区域内的矩谐坐标系下的地磁场分量值；

其中，所述测试点的地磁场分量值为Bx、By、Bz，矩谐坐标系下的地磁场分量值B_x'＝-X'cosθ₀cosλ₀-Y'cosθ₀sinλ₀+Z^/sinθ₀，

B_y'＝-X'sinλ₀+Y'cosλ₀，B_z'＝-X'sinθ₀cosλ₀-Y'sinθ₀sinλ₀-Z^/cosθ₀，其中，

X'＝-B_xcosθcosλ-B_y sinλ-B_zsinθcosλ，Y'＝-B_x cosθsinλ+B_y cosλ-B_zsinθsinλ，

Z'＝B_xsinθ-B_zcosθ，θ₀,λ₀分别为矩形坐标的原点的地理坐标系的余纬和经度。

步骤105、建立地磁场模型；

其中，选取部分测试点作为地磁场模型的精度检查点，根据剩余测试点在所述矩形区域内的矩谐坐标以及矩谐坐标系下的地磁场分量值，构建地磁场模型，并使用奇异值分解算法对模型系数矩阵进行求解计算；

根据矩谐坐标系下的拉普拉斯方程方程的解即为地磁场模型的构建公式，可以表达为

$\begin{array}{c}V=\mathrm{Ax}+\mathrm{By}+\mathrm{Cz}+\underset{q=0}{\overset{N_{\max }}{\mathrm{\&Sigma;}}}\underset{m=0}{\overset{q}{\mathrm{\&Sigma;}}}[E_{\mathrm{mn}}\cos \left(\mathrm{mvx}\right)\cos \left(\mathrm{nwy}\right)+F_{\mathrm{mn}}\cos \left(\mathrm{mvx}\right)\sin \left(\mathrm{nwy}\right)+G_{\mathrm{mn}}\sin \left(\mathrm{mvx}\right)\cos \left(\mathrm{nwy}\right)\\ +H_{\mathrm{mn}}\sin \left(\mathrm{mvx}\right)\sin \left(\mathrm{nwy}\right)]e^{\mathrm{uz}}\end{array}$

公式中，v＝2π/L_x,w＝2π/L_y,n＝q-m,N_max为模型的阶数，L_x和L_y分别为矩形区域的南北和东西长度，A、B、C、E_ij、F_ij、G_ij、H_ij为待定系数，共有2N_max(N_max+1)+3个。

地磁场模型的构建公式的矩阵形式可以表达为

其中，d为地磁场测试值，A为系数矩阵，M为地磁场模型矩阵，使用奇异值分解算法对模型系数矩阵A进行求解。

步骤106、对地磁场模型的精度进行验证，判断精度是否符合要求；

其中，通过地磁场模型的精度检查点对建立的地磁场模型的精度进行检查，判断地磁场模型的精度是否符合要求，当地磁场模型的精度符合要求时，跳转至步骤107，地磁场模型构建完成。当精度不符合要求时，跳转至步骤102，调整参数(研究区域周边的预设范围的大小、三角剖分的疏密程度、地磁场模型的截止阶数等)，重新进行模型构建。

步骤106具体可以包括如下步骤：

步骤106-1、将地磁场模型的精度检查点的矩谐坐标带入到地磁场模型中，计算出地磁场模型的精度检查点的磁场各分量场值d'；

步骤106-2、计算地磁场模型的精度检查点的磁场各分量场值d'与检查点的实测磁场各分量值d的方差；

其中，方差等于 $\mathrm{RMS}=\frac{\mathrm{\&Sigma;}[{(d_1-{d_1}^{\&prime;})}^2+\&CenterDot;\&CenterDot;\&CenterDot;{(d_m-d_m^{\&prime;})}^2]}{m}.$

步骤106-3、判断所述方差是否小于预设误差阈值。

步骤107、地磁场模型构建完成；

通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到本发明可借助软件加必需的硬件平台的方式来实现，当然也可以全部通过硬件来实施，但很多情况下前者是更佳的实施方式。基于这样的理解，本发明的技术方案对背景技术做出贡献的全部或者部分可以以软件产品的形式体现出来，该计算机软件产品可以存储在存储介质中，如ROM/RAM、磁碟、光盘等，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例或者实施例的某些部分所述的方法。

以上对本发明进行了详细介绍，本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

Claims (8)

1.一种地磁场高精度建模方法，其特征在于，包括：

获取测试点的地磁场测试值，并计算所述测试点的地磁场值，将所述地磁场测试值与所述地磁场值做差，以获取所述测试点的地磁场异常值；

根据所述测试点的分布情况和有限元法中的三角剖分方法，选取研究区域周边预设范围的测试点参与地磁场模型的构建，计算研究区域周边预设范围的测试点的地磁场值，确定矩形区域的大小、矩形坐标的原点以及地磁场模型的截止阶数；

将所述测试点的地理坐标转换为所述矩形区域内的矩谐坐标；

将地理坐标系下所述测试点的地磁场分量值转换为所述矩形区域内的矩谐坐标系下的地磁场分量值；

选取部分测试点作为地磁场模型的精度检查点，根据剩余测试点在所述矩形区域内的矩谐坐标以及矩谐坐标系下的地磁场分量值，构建地磁场模型，并使用奇异值分解算法对模型系数矩阵进行求解计算；

通过地磁场模型的精度检查点对建立的地磁场模型的精度进行检查，判断地磁场模型的精度是否符合要求，当精度符合要求时，地磁场模型构建完成。

2.根据权利要求1所述的方法，其特征在于，当精度符合要求时，调整相应参数，根据所述测试点的分布情况和有限元法中的三角剖分方法，选取研究区域周边预设范围的测试点重新参与地磁场模型的构建，所述相应参数包括研究区域周边的预设范围的大小、三角剖分的疏密程度、地磁场模型的截止阶数。

3.根据权利要求1所述的方法，其特征在于，所述获取测试点的地磁场测试值，并计算所述测试点的地磁场值，将所述地磁场测试值与所述地磁场值做差，以获取所述测试点的地磁场异常值，包括：

获取测试点的地磁场三分量的测试值；

利用IGRF国际地磁参考场模型和所述测试点的地理坐标，计算所述测试点的地磁场的三分量值，其中，所述IGRF模型取小于14的低阶项；

将地磁场三分量的测试值与地磁场的三分量值做差，以获取所述测试点的地磁场三分量的异常值。

4.根据权利要求1所述的方法，其特征在于，所述根据所述测试点的分布情况和有限元法中的三角剖分方法，选取研究区域周边预设范围的测试点参与地磁场模型的构建，包括：

根据研究区域整体面积及平均点距，确定研究区域周边的所述预设范围；

根据有限元法中的三角剖分方法，选取所述预设范围内的测试点参与地磁场模型的构建。

5.根据权利要求1所述的方法，其特征在于，所述根据剩余测试点在所述矩形区域内的矩谐坐标以及矩谐坐标系下的地磁场分量值，构建地磁场模型，并使用奇异值分解算法对模型系数矩阵进行求解计算，包括：

根据矩谐坐标系下的拉普拉斯方程获得地磁场模型的构建公式的矩阵形式为

其中，其中，d为地磁场测试值，A为系数矩阵，M为地磁场模型矩阵；

使用奇异值分解算法对模型系数矩阵A进行求解。

6.根据权利要求1所述的方法，其特征在于，通过地磁场模型的精度检查点对建立的地磁场模型的精度进行检查，判断地磁场模型的精度是否符合要求，包括：

将地磁场模型的精度检查点的矩谐坐标带入到地磁场模型中，计算出地磁场模型的精度检查点的磁场各分量场值d'；

计算地磁场模型的精度检查点的磁场各分量场值d'与检查点的实测磁场各分量值d的方差；

判断所述方差是否小于预设误差阈值。

7.根据权利要求1所述的方法，其特征在于，所述将所述测试点的地理坐标转换为所述矩形区域内的矩谐坐标，包括：

将所述测试点的地理坐标通过公式转换为所述矩形区域内的矩谐坐标，其中，地理坐标包括余纬θ和经度λ，矩谐坐标 θ₀,λ₀分别为矩形坐标的原点的地理坐标系的余纬和经度，R_e为地球半径。

8.根据权利要求1所述的方法，其特征在于，所述将地理坐标系下所述测试点的地磁场分量值转换为所述矩形区域内的矩谐坐标系下的地磁场分量值，包括：

将所述测试点的地理坐标系下所述测试点的地磁场分量值通过公式转换为所述矩形区域内的矩谐坐标系下的地磁场分量值，所述测试点的地磁场分量值为Bx、By、Bz，矩谐坐标系下的地磁场分量值B_x'＝-X'cosθ₀cosλ₀-Y'cosθ₀sinλ₀+Z^/sinθ₀，B_y'＝-X'sinλ₀+Y'cosλ₀，B_z'＝-X'sinθ₀cosλ₀-Y'sinθ₀sinλ₀-Z^/cosθ₀，其中，X'＝-B_xcosθcosλ-B_ysinλ-B_zsinθcosλ，Y'＝-B_xcosθsinλ+B_ycosλ-B_zsinθsinλ，Z'＝B_xsinθ-B_zcosθ，θ₀,λ₀分别为矩形坐标的原点的地理坐标系的余纬和经度。