邢琮琮，徐行

中国地质调查局广州海洋地质调查局，广州 510760

地磁场模型的计算与研究是地磁学的重要研究内容之一，它在地球科学、空间科学、地球物理勘探、岩石层物理学和地球深部研究等许多领域都有重要的学术意义和广泛的实际实用价值[1]。在海洋磁测中，地磁日变改正是资料处理中的一个重要环节，决定了后期地磁测量数据处理成果的质量，影响着地球物理勘探成果的学术水平。通常，磁法勘探需在测区附近布设地磁日变观测站，同步收集地磁日变数据[2]。而在深远海域中解决此类问题一直受技术方法的制约。针对这种情形，前人先后采用过“动态差分”和“静态差分”法来消除地磁日变影响。其中，所谓的“动态差分”法就是用双磁探头模式的地磁梯度测量方法来消除测量过程中的地磁日变影响；“静态差分”法是用定点的地磁日变观测数据去改正动态的海洋地磁测量值[3-4]。前者在小比例尺海洋磁测中的实践应用效果不好；而后者的效果虽好，但存在潜标系统设备收放复杂、采集成功率不稳定以及海底观测数据受海水中的感应地磁场干扰等因素的影响[5]。对于地磁场模型的使用，有国内学者通过研究认为采用综合模型一并进行日变改正、正常场及其长期变化改正,可以更好地解决日变改正中遇到的问题[6]；海外学者在海洋磁测中使用CM4和IGRF地磁场模型计算值来改正海洋磁测中的日变影响，在大塔木海山调查中取得了成功[7]。但由于相关成果中没有做与定点地磁日变站观测数据改正的对比分析，因而未能对其方法及精度做系统的评价。

随着人们对全球地磁场模型研究程度的加深和研究资料的不断积累，以及参与地磁场模型计算的卫星磁测数据质量的提高，将地磁场模型的计算成果应用到无限海域地磁测量的时机已成熟。然而，将地磁场模型计算值用于海洋磁测中的日变改正与其他方法的差异如何？测量成果的精度怎样？本文作者通过对地磁场模型的计算与应用研究，结合海上地磁测量和资料处理的具体实践，来回答这些问题。

1 数据和方法

1.1 数据来源

本文共收集了3类数据：①西太平洋某工区的实测数据（海洋磁测的测线数据及地磁日变观测站的观测日变数据）；②广东省肇庆地磁台K指数数据；③地磁场模型计算日变数据（包括更新的CM4模型卫星数据和最新的IGRF模型数据）。

地磁日变观测站（下文简称日变站）是由一条4 000 m长的主绳，6个正浮力玻璃浮球，2台加拿大生产的Sentinel磁力仪，3台声学释放器和1个自制的218 kg（海水中的重量）的水泥块组成。日变站的结构图见图1。其中，Sentinel磁力仪主要技术指标为：直接测程：18 000～120 000 nT，灵敏度：0.015 nT，分辨率：0.001 nT，绝对精度：0.2 nT；该设备的工作温度范围为−25～+60 ℃，工作水深可达 6 000 m，数据存储量为1×107个读数，数据采样率可达10 Hz，无方向性误差，无温度漂移，无测量盲区。

西太平洋某工区的日变站布放位置为16°23′N、129°35′E，投放点水深 4 280 m，图2 为日变站的布放位置图。其中，Sentinel磁力仪沉放深度400 m，采集的数据时间跨度为24.7 d。海洋磁测的测线共15条，包括8条东西向主测线和7条南北向联络测线。海洋磁测中使用了Seaspy海洋磁力仪，与Sentinel磁力仪同一厂家，其技术指标也相同。以上实测数据从广州海洋地质调查局的调查项目中收集，野外资料的技术要求和质量均满足海洋地质调查规范的要求。

图1 地磁日变观测站结构图Fig.1 The structure of geomagnetic diurnal observatory

图2 地磁日变观测站布放位置图Fig.2 Location of the geomagnetic diurnal observatory

K指数是单个地磁台用来描述每日每3小时内的地磁扰动强度的指数，称为3小时指数或磁情指数。它是一种定量的分级指数，从0到9共分10级，数字越大表示地磁扰动越强。每日分为8个时段（世界时），00～03时为第1时段，03～06时为第 2 时段等，···，21～24 时为第 8 时段，每个时段确定一个K值[8-9]。本文使用了广东省地震局肇庆基准地磁台（22°58.20′N、112°27.00′E）的 K 指数数据，相应时段与本文中的日变站记录时段一致。K指数的波动与日变站观测数据所反映出的磁情变化基本一致，数据质量可靠。

CM4模型的卫星数据使用了1960—2002年间POGO卫星和CHAMP卫星的标量数据，Magsat卫星和Oersted卫星的标量与矢量数据、更新的磁层环电流指数（Dst指数）和太阳通量指数（F10.7指数）；台站数据使用了地磁台站每月最静时时均值，分为两种不同采样率[10-11]：①1960—1985年最接近01:00地方时的时均值（OHM-1AM）；②在POGO和MAGSAT卫星任务期间，每月磁静日中以2小时为间隔的台站时均值（OHM-MUL）。IGRF模型使用了国际上每5年发表一次的球谐系数，由一组高斯球谐系数、）和年变率系数、）组成[12-13]。

1.2 模型计算方法

CM4模型的适用时间为1960年1月—2002年6月，它的计算分为内源场与外源场两部分。对于近地表的观测，内源场是由地核磁场、地幔感应电流和地壳磁场组成；外源场是由磁层场及其感应场、电离层场及其感应场和环形电流组成[14-15]。为突破CM4模型的适用时间限制，根据内源场的成分组成及其随时间变化缓慢的性质，本文使用全球地磁场模型国际地磁参考场（IGRF）替代；外源场则是更新CM4模型中相应的卫星数据之后进行计算[7,10]。由于IGRF模型计算内源场和CM4模型计算外源场的结果均为三分量数据，最终将它们合成为地磁总场值。具体的计算公式和方法见作者的前期成果[16]。

1.3 磁测数据的处理和精度评价

海洋磁测的原始数据需要经预处理以及一系列的改正才可得出磁异常值。根据海洋地质规范中的处理方法[17]，其数据处理的流程图见图3所示。

海洋磁测精度评价的重要信息来源是主测线和联络测线的交点差。一般采用交叉点的内符合精度来衡量数据质量，其统计方法[18-19]如公式1所列：

2 数据处理

2.1 整体日变曲线对比

对比的地磁日变曲线共两种：一是日变站观测的记录曲线（下文简称，实测日变曲线），二是用地磁场模型计算获得的计算曲线（下文简称，计算日变曲线）。其中，日变站数据的时段为世界时2018年12月4日06:00至2018年12月28日23:00，经去除跳点、插值及滑动窗口平均之后，得到实测日变曲线。在CM4模型和IGRF模型计算中，所设参数均为本文日变站布放位置和采集时段，高程为0 m；其中，IGRF模型在采集时间内每天计算一次。最终将CM4模型和IGRF模型合成地磁总场值得到计算日变曲线。

图4为实测日变曲线与计算日变曲线的对比图。从K指数的波动程度分析，2018年12月4日—11日和16日—19日，K值在1～3之间波动；12日—15日，K值多为0～2，偶有3的跳变；20日—26日，K值最大为2，是磁静日时段；27日—28日的K值为3～4，为磁扰日时段。

比较实测日变曲线和计算日变曲线（图4）：在形态方面，两者变化趋势相似；其中的低频部分拟合较好，高频细节部分有差异。在幅值方面，两者在曲线的极大值和极小值位置有一定的幅值差，计算日变曲线的幅值变化幅度大多数情况下小于实测日变曲线；特别是在出现磁扰（K指数较大）的位置上，其幅值变化的差异更加明显。在相位方面，两者基本相同，计算日变曲线略有延后。

2.2 典型时段日变曲线及其改正效果对比

在图4中分别选择K指数在不同数值区间的时段作图分析，a、b、c区域分别代表K指数在0～2、2～3和3～4之间的细节图。

图3 海洋磁测数据处理流程图Fig.3 The processing flow chart of marine magnetic data

图4 地磁日变曲线对比图 （红色线为计算日变曲线，蓝色线为实测日变曲线，黑色线为K指数；a、b、c虚线框为典型细节图范围）Fig.4 Comparison of geomagnetic diurnal variation curves （The red line is the calculated curve,the blue line is the measured curve,and the black line is the K index; The dotted lines of a,b,and c are those showing typical details）

图5 日变曲线细节对比图 （红色线为计算日变曲线，蓝色线为实测日变曲线，黑色线为K指数）Fig.5 Comparison of diurnal variation curve details （The red line is the calculated curve,the blue line is the measured curve,and the black line is the curve of K index）

图5a的时段为12月22日06:00至24日18:00。K指数有3个时段最大为2，4个时段最小为0，其他时段为1，为磁静时段。两曲线形态基本相同，拟合程度较好，曲线的细节波动有差异，幅值差小，在两个波峰位置处，两曲线的极大值差值分别为2.22和3.56 nT，曲线的后半段可看出模型计算日变曲线的相位有延后。

图5b的时段为12月7日00:00至8日18:00。K指数有5个时段为3，9个时段为2，有轻度磁扰。虽两曲线形态相似，但存在一定的幅值差。其中的计算日变曲线较为圆滑，缺失高频成分。

图5c的时段为12月27日21:00至28日22:00。K指数有2个时段最大为4，1个时段最小为2，其他时段为3，是典型的磁扰时段。虽然两曲线总体变化趋势类似，但不仅形态相差大，而且幅值差也较大。在本时段中的计算日变曲线幅值变化幅度为 32.14 nT，小于日变站的 89.14 nT。

在西太平洋某工区的15条海洋磁测测线的数据处理过程中，分别使用了模型计算数据和日变站实测数据进行地磁日变改正，其他所用数据以及处理方法均相同。使用模型计算数据得到磁异常简称为“模型磁异常”，使用日变站实测数据处理得到的磁异常简称为“实测磁异常”。取其典型海洋磁测的测线进行分析（图6）。

图6中a、b、c 3条测线分别对应图5中 a、b、c 3个时段，图中纵轴幅值比例相同。图6a两条磁异常曲线重合度高，相关性较好，形态和相位相同，幅值在K指数为2的时段存在微小差异。图6b两条磁异常曲线形态和相位基本相同，相关性一般，幅值在K指数由3降至2前后存在一定差异。图6c两条磁异常曲线形态和相位相似，相关性较差，两者的幅值重合度较低，有较大差异。

图6 测线磁异常对比图 （红色线为模型磁异常曲线，蓝色线为实测磁异常曲线，黑色线为K指数）Fig.6 Magnetic anomaly comparison of survey line （The red line is the magnetic anomaly curve of the model,the blue line is the measured magnetic anomaly curve,and the black line is the K index）

3 讨论

3.1 模型计算和日变站曲线分析

根据图4中的计算日变曲线与实测日变曲线对比分析，两者形态相似，低频拟合好，高频有缺失，幅值变化相对平缓，相位略延后。为评价两曲线的拟合程度，计算两者的相关系数，如公式2所列：

式中r为相关系数，Cov(X,Y)为X与Y的协方差，Var[X]为X的方差，Var[Y]为Y的方差，其中X、Y分别代表计算日变曲线与实测日变曲线。经计算，两者的相关系数为0.85，拟合程度较高。

分析两曲线差异的原因，应为计算外源场时使用Dst指数时间分辨率低，导致计算结果高频不足；而且Dst指数是由4个近赤道中低纬地磁台H变化均值归算到赤道后求平均得到的。当有地磁扰动时，不能精确反映计算点实际磁扰情况，导致曲线幅值变化相对平缓。

图5显示了K指数在不同区间时计算日变曲线与实测日变曲线的细节对比。计算a，b，c细节图的曲线相关系数，分别为0.95，0.93，0.87。说明当磁扰变大时，K指数相应增大，变化地磁场复杂，幅值差异逐渐增大，导致了两种曲线的相关程度逐渐降低，代表3≤K≤4的c图的相关系数明显低于a、b两图。

表1为在不同K指数下两类曲线的幅值差绝对值的均值与最大值统计表，表中差值绝对值与K指数为非线性正相关关系，K=4时差值绝对值的均值明显增大，结合图5中c图的相关系数明显减小，证明K=4时的模型计算精度较差。

表1 K指数与两种日变曲线差值绝对值统计表Table 1 Statistical of K index and absolute value of difference between two kinds of diurnal variation curves

3.2 两种磁异常的幅值差分析

综合图6中a、b、c 3个图的分析，随着K指数的增大，模型磁异常与实测磁异常的差异随之增大。图7反映了图6中3条典型测线的模型磁异常曲线和实测磁异常曲线之间的差值变化曲线。表2为两种磁异常差值绝对值的平均值和最大值统计表。显然，两种磁异常的差异主要体现在幅值方面。

在磁测数据处理时，尽管使用了两种不同类型的地磁日变数据，但计算方法和过程是相同的，所以得到的磁异常差异归咎于地磁日变数据类型的差异。在图7中，测线a的统计值均小于8 nT，磁异常差异较小；测线b的统计值多数小于10 nT，少数为10～20 nT，磁异常差异略有增大；测线c的统计值多数为 10～30 nT，少数小于 10 nT，磁异常差异较大。分析表2的相关信息，当K指数从0～2变化至2～3时，磁异常差值绝对值的平均值由3.08 nT增大至5.23 nT，增长量不大；而当K指数从2～3变化至3～4时，其平均值是阶跃式增长，从5.23 nT增大至 12.98 nT。

表2 两种磁异常差值绝对值统计表Table 2 Statistical of absolute value of difference between two kinds of magnetic anomalies

图7 两种磁异常差值绝对值分析图 （红色线为磁异常差值绝对值曲线，黑色线为K指数）Fig.7 Absolute value analysis of difference between two magnetic anomalies （The red line is the absolute value curve of magnetic anomaly difference,and the black line is K index）

William W.Sager等在对大塔木海山的海洋磁测中同样使用了CM4模型计算外源场，筛选数据的标准为Kp指数（地磁行星指数）大于5。由于Kp指数反映的是全球磁场扰动水平，而K指数反映的是单个地磁台的地磁扰动水平[20-21]，本文认为调查区附近的地磁台站K指数更适于反映调查区的地磁扰动。

系统地对比两种不同类型的地磁日变数据以及典型测线的数据处理结果，可认为：在0≤K≤2时，用地磁场模型计算值进行日变改正的成果精度较高；在3≤K≤4时，其成果精度要低一些；在K＞4时，需谨慎使用，甚至建议不使用。

3.3 精度评价

交点差及其内符合精度是评价数据精度的重要指标。在15条海洋磁测测线的数据处理中，有8条东西向主测线，7条南北向联络测线，测线交点共计56个。图8为处理得到的实测磁异常交叉点和模型磁异常交叉点差值空间分布图。图8中两者的交点差均为调差前交点差。

图8 交点差对比图 （交点右上红色值为模型交点差，左上蓝色值为实测交点差）Fig.8 Comparison of intersection difference （The red value on the upper right of the intersection point is the model intersection point difference,and the blue value on the upper left is the measured intersection point difference）

在56个交点中，有37个交点的模型交点差绝对值大于实测交点差绝对值，占总数的三分之二。模型交点差绝对值的均值为5.74 nT，最大值为20.29 nT，内符合精度为 5.04 nT；实测交点差绝对值的均值为 2.93 nT，最大值为 9.12 nT，内符合精度为2.55 nT。由对测线磁异常的讨论可知，在K＞3时模型计算数据的处理结果可信度低，去除处于12.27—28日磁扰时段的6个交点再次统计，模型内符合精度为4.46 nT，实测内符合精度为2.64 nT。两者的精度变化表明：在3≤K≤4时，实测数据的改正效果比模型计算数据好，所以去除磁扰时段交点后实测内符合精度变化小而模型内符合精度明显提升，从中验证了模型计算数据在海洋磁测资料处理中的改正能力。

根据海洋磁力测量规范中的精度要求为调差后交叉点内符合精度小于4 nT，对本文中的56个交点进行了调差。调差后，用模型计算和实测数据做海洋磁测处理的内符合精度分别为3.74和1.68 nT，两者均符合要求，但前者的精度低于后者。

4 结论

（1）用计算数据和实测数据做海洋磁测的处理时，在调差前，两者的内符合精度分别为5.04和2.55 nT，基本满足海洋地质规范的要求；在调差后，两者的内符合精度分别为3.74和1.68 nT，均符合海洋地质规范的要求，但前者的精度低于后者。

（2）磁静时，计算数据和实测数据拟合程度高；磁扰时，两者幅值差异较大，磁扰越强相关系数越小。在0≤K≤2时，用计算数据进行日变改正的成果精度较高；在3≤K≤4时，成果精度低；在K＞4时，建议谨慎使用或不使用。

致谢：感谢广东省地震局肇庆地磁台提供的论文K指数数据，云南省地震局预报中心姚休义博士在计算方法上给予的重要帮助，海军测绘研究所黄漠涛研究员在论文修改中提出的宝贵意见。