高精度磁测梯度改正的 Matlab程序实现

陈林 李暖 陈伟刚

(青海有色地质矿产勘查局地质矿产勘查院，青海西宁 810007)

地球磁场磁化作用于岩(矿)石所产生的磁性，是研究磁异常的理论基础，通过研究现代样品的原生剩磁推测地质时期的磁性特征，使得研究地质构造的演化提供了可能(罗效宽等，1991)［1］。当进行大面积高精度磁测工作时，需要进行正常梯度改正。此时要用查全国地磁图的办法将不能满足正常场梯度改正的精度要求，而球谐分析方法则可满足这一要求。目前，由于球谐分析方法公式复杂，参数较多，需要花大量的人力物力才能进行正常场的梯度改正，同时，由于市面上出现的同类软件价格昂贵，且购买一个软件只能安装在一台电脑上。因此，作者利用Matlab软件编程的简易性、函数的多样性、数学计算的快速性等特点，编写了一套高精度磁测梯度改正的简易程序，解决了这一难题。

1 地球磁场球谐分析的原理

1838年由高斯首先提出球谐分析的方法，该方法是表示全球范围地磁场的分布及其长期变化的数学方法。假设地球是均匀磁化球体，球体半径为R，N为地理北极。这里设地球旋转轴与地磁轴重合，故N也表示地磁北极。若采用球坐标系，如图1所示。

图1 球极坐标系示意图

坐标原点为球心，球外任一点的地心距为r，余纬度为θ(θ= 90°－ψ，ψ为纬度)，经度为λ。则在地磁场源区之外空间域坐标系(r，θ，ψ)中，磁位u的拉普拉斯方程可以写成:

对上式采用分离变量法，即令U(r，θ，λ)=R(r)·H(θ)·Φ (λ)，则可解得拉普拉斯方程的一般解，从而可分别获得其内源场和外源场的磁位球谐表达式。若设外源场磁位为零，则内源场的磁位球谐一般表达式为:

其中，Pn(cosθ)为勒让德多项式。

由此，便可得到相应三个轴向磁场强度的三分量的球谐表达式。地磁场感应强度的三个分量即北向水平分量X、东向水平分量Y、垂直分量Z(注意这里定义X轴指北为正，Z轴向下为正):

2 Matlab编程及应用

2.1 技术路线

由上述球谐分析方法的原理，需要进行n阶(n+1)次的迭代计算，因此，地磁场磁感应强度的三个分量X，Y，Z值在Matlab软件中将分别进行两次套合循环计算。程序运行路线见图2。

2.2 Matlab软件编程源代码

根据上述编程技术路线，其程序源代码如下:

图2 Matlab软件程序运行路线图

其中，需要注意的是:

2.3 实际应用

为验证程序计算结果，实际选取4个点进行计算，计算结果见表1。

表1 实际计算结果表

3 结语

1)本程序经过在高精度磁测实际操作中，快速解决了正常场梯度改正的问题，并且有效精度能达到4位小数，从而快速解决了利用人力所需大量时间的问题。2)Matlab软件由于含有大量的数学计算函数，即便是勒让德多项式这一复杂度计算，在编程过程中，也只需一个函数代码即能解决。3)由于本程序只是为解决人力大量计算的辅助简易程序，并且Matlab软件是以C语言及数理方程为基础，其编程有相应的数学及语言格式，人机窗口操作上则为不便，因此，使得程序存在变量较多，操作繁琐这一局限性。

［1］ 罗孝宽，郭绍雍.应用地球物理教重力：磁法篇［M］.北京：地质出版社，1991.