2012, Vol. 55
2. 中国科学院空间天气学国家重点实验室, 北京 100080;
3. 中国科学院研究生院, 北京 100049
, DU Ai-Min1,2
, XU Wen-Yao1, HONG Ming-Hua1, CHEN Geng-Xiong1, LUO Hao1, ZHAO Xu-Dong1, BAI Chun-Hua1, WANG Yuan1
2. State Key Laboratory of Space Weather, Chinese Academy of Sciences, Beijing 100080, China;
3. Graduate University of Chinese Academy of Sciences, Beijing 100049, China
地磁场的组成主要包括主磁场部分和岩石圈磁场部分.其中主磁场占据主要的比重.排除外源场的影响，来源于岩石圈磁性物质的剩余磁场也能对地表磁场观测数据带来重要的影响.岩石圈磁异常特性通常表现在1000km 以上的水平空间尺度上，因此在小尺度范围内仍然保持着比较好的一致性.但浅层地壳(深度≤20km)岩石的分布变得相当复杂，所以在矿产资源密集分布的小尺度区域(≤10km)，该剩余磁场能引起的地表磁场强度变化往往能达到1000nT 量级，而用现有的全球或区域模型却无法表征如此小尺度的地磁场空间分布特性.
使用全球或区域的地磁场建模方法建立小尺度磁异常模型存在很多的局限.首先，球谐分析方法作用在全球区域，其水平空间分辨率理论上依赖于模型的阶次数.NGDC720模型的球谐函数最高达720阶次，能反映的磁场最小空间尺度为55km，远大于上述的小尺度范围[1].第二，高阶次的模型需要大量的实际观测数据拟合，而大部分的全球或区域模型都要借助卫星数据实现高的空间分辨率[2-4].小尺度的磁场模型一般需要反映近地面的磁场分布，数据的不匹配将影响模型上下延拓的精度.第三，许多区域模型例如矩谐分析和球冠谐分析方法[5-6]，由于理论上的缺陷，运用在较小范围的区域时，模型的精度和空间分辨率会明显降低，不能满足小尺度地磁场模型高分辨率和高精度的要求[7-8].
多项式方法早在1969年就应用在地磁场建模当中[9]．当时使用的是泰勒多项式方法，泰勒多项式虽然能用指数形式的基函数表征物理量，但由于不是正交多项式，在反演过程使用矩阵计算会大大增加计算量，不适用于高阶运算．多项式模型没有从磁场的位场理论出发，而是单独拟合磁场的一个分量，虽然仅局限于二维平面，但由于模型精度较高，而且可以拟合磁场的磁偏角、磁倾角等非线性分量，因此在矿产资源开发，石油勘探等实际应用中具有较高的小尺度磁场建模能力．勒让德多项式作为一类具有正交特性的多项式形式，除了能大大简化模型运算，还可以进一步提高模型的截止水平，反映更小尺度的磁场分布．赵建虎等人运用勒让德多项式方法建立海洋局域的地磁模型[10]，但截止水平仅达到8阶．本文通过实验数据建模，介绍勒让德多项式在小尺度范围内建立高截止水平磁场模型的情况，提出一些数据预处理的过程，并通过误差分析，讨论勒让德多项式建模方法的使用问题.
2 方法介绍 2.1 勒让德多项式介绍勒让德多项式是在球坐标下利用分离变量法解拉普拉斯方程得到
|
(1) |
勒让德多项式在[-1,1]之间具有正交性，即
|
(2) |
式(1)中x为自变量，k为勒让德多项式的阶数，m为求和各项的序数.用磁场测量值确定多项式系数是一个多元回归的问题[11].当组成多项式的各项不是正交的时候，回归系数存在一定的相关性，以致剔除一个自变量后，还必须重新计算.但因为勒让德多项式具有正交性，能使系数矩阵变成对角阵，所以不仅可以大大简化计算，而且能消去回归系数之间的相关性.
用多项式拟合磁异常分布时，常常遇到这样的问题:低次多项式虽然简单、稳定，但只能表示变化平缓、结构简单的磁场；为了很好地表示一个复杂的磁异常分布，特别是当计算区域内包含一些小尺度磁异常时，必须使用高次项，而高次多项式的拟合和插值过程有数值不稳定的缺点，在没有数据点控制的地方，甚至会给出毫无意义的结果.但如果能解决上述问题，运用高次项建立模型，就可以把区域当中小尺度的磁异常区域表示出来.
2.2 建模的基本方法采用实验数据作为方程组的输入，把每一个数据用同一组参数的勒让德多项式表示，建立方程组.由于勒让德多项式只在[-1,1]区间内具有正交性，所以首先把每个数据的坐标值进行归一化的处理，使其都落在[-1,1]的区间里，然后利用最小二乘原则对系数矩阵进行求解，得到最优的参数解.
若地磁的总强度为B，则基于勒让德多项式所建立的区域磁异常基本模型方程为[10, 12]
|
(3) |
其中B为观测数据，aij为要求的参数，N是勒让德多项式的截止阶数，Pi(Δφ)为i阶的勒让德函数，Δφ 和Δλ是归一化以后观测数据所在的纬度和经度值.公式(3)是由两套互相正交的勒让德多项式组构成，分别以经度和纬度作为位置变量，代表在各自的一维方向上的基函数波形，因此不同阶数的勒让德多项式，可用于表征磁场不同频率成分的基本要素，而系数aij可认为是不同频率成分所占的能量大小.由于勒让德多项式只在[-1,1]区间内具有正交性，所以需要把观测数据的经纬度坐标作归一化处理:
|
(4) |
然后把一系列的观测数据作为输入建立方程组，利用已知的Pi(Δφ)或Pi-j(Δλ)，简化方程组为B=AX.其中B为观测数据序列，B= [B1，…，Bk] T；A为勒让德函数组成的矩阵，
|
(5) |
X为要求参数的矩阵，X=[x00，x01，…，xN0].
利用最小二乘原则可求得参数矩阵:
|
(6) |
其中C-1为协方差矩阵，这里定义为单位矩阵.
3 利用实际观测数据建立模型 3.1 实验区域选取和观测数据来源实验数据是来自内蒙古某区域，范围为20km×20km.其中图 1 给出的是本次实验观测数据点在二维平面上的分布图，总共包含143个数据点.本次实验采用的数据均为磁场总强度.
|
图 1 原始观测数据二维分布图，经纬度坐标均经过处理，维持原来的尺度大小，数据点的不同颜色代表不同的磁场总强度值 Fig. 1 Two dimensional distribution of original survey data with different colors for different values of the total geomagnetic field intensity |
该实验区域的数据点空间分布具有两个特点:首先是点分布区域范围比较小，经纬度的最大值和最小值差值都约为0.2°，空间尺度为20km；其次数据点磁场强度的空间分布结构:该区域大致可分为上下两个部分.上部的数据点比较集中，同时存在几个磁场强度异常较大的区域.而且分布非常不均匀，磁场结构复杂；而下部的数据则相反，数据的分布集中且均匀，其磁场的位形大致平缓，但存在几个数据空白的区域.这些典型的特征都经常能在实际的地磁数据测量中碰到.
实验数据测量仪器包括:质子磁力仪、DI观测仪、高精度差分GPS.
GSM-86 是一种标量磁力计，具有很强的准确度(0.1nT)和很好的稳定性(漂移量为0.05nT/a).这一仪器的核心优势功能在于具有长时间的稳定性和高精度以及高分辨率和低噪声.
DI观测仪能够满足实验和地磁测量中对精度和多功能性的各种要求.它可以与微型的轴向或横向磁通门探针或安装在经纬仪上的探针配合使用来对地磁场进行高精度测量.其范围从1nT 至2 mT(依赖于探针)，极好的线性，良好的温度稳定性和较低的噪声输出.
差分GPS 则采用Ashtecha的ProMark2TM.它是一种将后处理系统、厘米级精度的静态测量系统以及3~5 m 精度的实时导航系统有机结合的GPS接收机产品，能够为各种工程应用提供高质量的测量数据.
观测数据均经过临近固定台站的日变校正.
3.2 初步的模型拟合若只采用观测数据的总共143个点作为模型参数计算的输入，由于实验观测数据比较少，从理论上来说，勒让德多项式模型最大的截止阶数是15 阶，才能保证未知参数的个数少于参与方程运算的观测数据个数.随着阶数的增加模型拟合观测数据的效果会越来越好，由于15阶模型要求的未知参数过于接近方程组输入观测数据的个数，而14阶以下的误差水平比较接近，所以不妨分别选用5阶和10阶的模型在区域内产生均匀网格数据，并制作等值线图.如图 2所示.
|
图 2 原始观测数据（a) 5阶模型和（b) 10阶模型磁场总强度等值线图，经纬度坐标均经过处理，维持原来的尺度大小，等值线间隔为200 nT Fig. 2 Contours of geomagnetic total intensity with intervals of 200 nT by (a) the 5-degree model and (b) the 10-degree model derived from original survey data |
比较图 2与图 1，可以发现5阶模型里面，原始观测数据中一些大磁异常点在低阶模型中得不到反映，模型所拟合出来的整体数据误差很大程度上来自这些大磁异常点的贡献；与此同时，实验区域中一些缺乏数据的地方在10阶的模型中失去控制，得到毫无意义的数据，产生严重的边界效应.如图 2b中所标的区域所反映的情况所示:区域1 和3 是因为周边的数据点从外到内有变小的趋势，导致模型在该区域中拟合出非常低的值；而区域2刚好相反，由于周边的数据点从外到内有变大的趋势，导致模型在该区域中拟合出非常大的值.通过初步的模型拟合结果可以看出，仅用原始的观测数据建立模型会带来两个问题，第一是精度:在低阶的时候拟合出来的数据与原始数据之间存在较大的误差；第二是边界效应:由于观测数据是非均匀的，拟合曲线简单地从受数据点控制的区域按照该区域的变化趋势延伸到没有数据点控制的区域，由于这些区域通常处在边界，所以容易导致在边界上出现非常异常的与实际不符的数据，特别当阶数比较高的时候这种情况会更加显现出来.所以有必要对原始的观测数据作合理的预处理.
3.3 进一步的数据处理和模型改进首先对原始的观测数据的边界进行预处理，在比较缺乏数据的边界上人为添加数据，添加的值与最接近边界的同纬度或同经度点的值保持一致，这样是避免数据点边界的梯度趋势在实验区域的边界上被放大，产生无意义的拟合数据.同时在个别异常的数据点附近添加数据点来控制该区域异常，降低对其附近空白的区域的影响.处理后的数据区域分布图如图 3所示.
|
图 3 加入控制点后的数据点分布图 彩色的实心圆圈为原始观测数据，红色和蓝色的 空心圆圈均为人为添加的控制点数据. Fig. 3 Data distribution map Red and blue hollow circles are the manually added boundary controlling points，solid circles are original data. |
图中可以看到，在数据区域的左边边界增加了一些数据点，如红色空心圆圈所示，而旁边一些小尺度的异常点用与其附近最近点的值相同的点去限制该区域梯度向空白区域发散，避免影响梯度平缓的区域，如蓝色空心圆圈内所示.
在保证异常点不会在数据空白区域以及边界发散的前提下，对区域的数据进行插值的预处理一方面可以在没有数据点控制的区域内添加合适的、不影响该区域磁场变化特征的数据点，另一方面使数据点的分布更为合理，有利于模型的数值拟合.本次实验采用的插值方法是最小曲率法[13].
最小曲率法(Minimum Curvature)是一种在地球科学中广泛应用的插值方法.其插值函数的几何图像是一个尽可能与原始数据点吻合的最平滑的曲面，相当于一个具有线性弹性的薄板.最小曲率法试图在准确还原原始数据的同时，生成尽可能圆滑的曲面[14].
假设在平面(xi，yi)上有n个数据点，最小曲率法的基础函数是
|
(7) |
其中d表示(x0，y0)和(x1，y1)之间的距离；当k=1时表示最小曲率法，k=2时表示薄片样条法.从而得到的插值结果如图 4所示.
|
图 4 原始数据预处理后的磁场总强度等值线图，等值线间隔为200nT Fig. 4 Contours of geomagnetic total intensity with interval of 200 nT after data preprocessing |
经过处理后的观测数据再经过插值以后的结果比较好地控制了边界.最后把插值以后的数据矩阵加上原始数据作为方程组的输入用于求解模型参数.此时可根据数据样本数的大小进一步提高模型的截止阶数，反映实验区域中的小尺度磁异常情况.
4 实验结果最后通过实验数据的预处理和提高模型阶数，得到截止水平为50阶的勒让德多项式模型，表 1为模型对观测数据进行不同阶数拟合估计得到的误差值.图 5为模型拟合原始观测数据的误差分布图.图 6为模型拟合的磁场总强度等值线图.
|
表 1 模型误差 Table 1 Modeling errors |
|
图 5 用插值后的数据建立50阶模型 拟合结果的误差分布图 Fig. 5 Distribution of misfits of the 50-degree model compared to the original data after interpolation |
|
图 6 50阶模型拟合结果的磁场总强度等值线图，等值线间隔为200 nT Fig. 6 Contours of geomagnetic total intensity with interval of 200 nT predicted by the 50-degree model |
我们利用勒让德多项式对小尺度(20km×20km)地磁场进行建模.由于该实验区域测量数据点分布不均匀，同时存在若干强异常点，所以对建立模型的方法的测试很有帮助.从实验结果可以看出，经过改进的模型其精度有很大的提高.只用原始的观测数据点建立的5阶模型只能描述4km 范围尺度的磁异常信息，而平均误差达到400nT 以上，所反映的区域分布不显著，而10阶的模型虽然能描述大概2km范围尺度的磁异常区域信息，但平均误差也超过了300nT，虽然能反映一些局部异常，但异常区域的拟合误差较大，而且存在严重的边界效应.
通过一定的边界处理、异常点处理、以及数据插值可以提高模型阶数并减低边界效应带来的影响.由于勒让德函数阶数与区间内的波数存在线性关系，因此可以推算50阶的模型能描述0.4km 尺度范围的磁异常，而模型拟合实验数据的平均误差小于30nT，异常点也得到很好地控制.拟合出的等值线图与模型输入数据的等值线图基本保持一致.由于模型中n=0阶代表着磁场的平均值，所以进一步利用模型结果计算区域内原始数据点的磁能密度平均值在n=1~50阶的分布情况，如图 7所示.可以看到模型n=10阶附近包含磁场数值的主要部分，n=30~40阶反映着低阶模型的主要误差来源，而n=50阶时所代表的磁能密度数值基本为零.因此选取50阶的截止水平已经可以很好地拟合原始数据，而再高阶的模型则容易出现不稳定的情况，即在没有控制点的区域产生不合理的结果.若数据分布均匀合理，同时基本处于一个平面内，则勒让德多项式方法均可适用，包括海洋测量和航空测量数据.
|
图 7 模型拟合原始数据磁能密度随模型阶数的分布图 Fig. 7 Distribution of predicted magnetic power density with the model degree |
个别的误差来源于小尺度的磁异常区域:从表 1所反映的误差情况可以看到，模型在对观测数据拟合过程中存在的个别大误差点也能对平均误差造成影响.通过误差分布图 5可以看到，误差较大的点均存在于数据点分布比较密集且差值较大的区域.在数据点增多而且截止阶数有限的情况下，这些梯度大的异常是难以准确地通过模型描述出来.边界效应控制与模型的截止阶数之间存在着矛盾:模型的截止阶数较低则反映的区域异常情况较差，而相反提高模型的截止阶数则容易引起边界的不稳定，而过多的边界数据预处理就会降低实验数据的可信度.因此观测数据预处理的方法只能对控制边界起作用，不能减小此类误差.这为模型阶数的选择带来不便，因为不同的区域出现磁异常的情况均不一样，没必要都采用很高阶数的模型去描述[15-16].虽然高阶的模型能显示出高精度的结构，但并不代表实际中的区域磁场就如模型所反映的那样复杂，所以如何合理地选取截止阶数也是非常值得考虑的问题.另外，模型误差还有可能来源于外部磁场.在磁场观测过程中，受外部电流体系影响而产生的感应磁场也能对本底固有的地磁场造成干扰(如Sq电流体系)，尤其是在磁扰日发生磁暴和亚暴期间.因此，需要对测量数据进行前期的校对和修正[17-24].
边界效应虽然能通过边界数据的预处理以及异常点的控制得到一定的抑制，但过多人为控制点的引入同样会改变磁场的分布位形.人为控制点的磁场值选取同样重要，最简单的做法是选取最临近边界的实际数据点作为参考值，或者通过临近的几个实际数据点所反映的磁场梯度作为参考.边界的数据处理通常不会要求特别严谨，实际测量过程中人们通常扩大测量范围，而实际关心的区域比实际测量的小，从而能通过建立模型取得比较可信的拟合结果.不同的区域的边界只能根据不同的情况作特殊的处理，本文并没有总结出一种比较通用的处理边界的方法.
不同插值的方法的使用同样影响模型拟合的情况，采用什么样的插值方法取决于该区域磁场的分布情况以及实验区域的大小、研究的磁场特征尺度.本文并没有估计不同的插值方法对反映磁场分布位形的影响.在经过数据插值的预处理后，勒让德多项式建模方法的拟合效果是本文关心的问题，而本文所采取的最小曲率法插值方法，则倾向于匹配原始数据点附近的磁场分布情况.
| [1] | Maus S, Sazonova T, Hemant K, et al. National geophysical data center candidate for the world digital magnetic anomaly map. Geochem. Geophys. Geosyst. , 2007, 8(6): Q06017. DOI:10.1029/2007GC001643 |
| [2] | Sabaka T J, Olsen N, Purucker M. Extending comprehensive models of the Earth’s magnetic field with rsted and CHAMP data. Geophys. J. Int. , 2004, 159(2): 521-547. DOI:10.1111/gji.2004.159.issue-2 |
| [3] | Maus S, Rother M, Hemant K, et al. Earth's lithospheric magnetic field determined to spherical harmonic degree 90 from CHAMP satellite measurements. Geophys. J. Int. , 2006, 164(2): 319-330. DOI:10.1111/gji.2006.164.issue-2 |
| [4] | Olsen N, Leuhr H, Sabaka T J, et al. CHAOS—a model of the Earth’s magnetic field derived from CHAMP, Orsted, and SAC-C magnetic satellite data. Geophys. J. Int. , 2006, 166(1): 67-75. DOI:10.1111/gji.2006.166.issue-1 |
| [5] | Alldredge L R. Rectangular harmonic analysis applied to the geomagnetic field. J. Geophys. Res. , 1981, 86(B4): 3021-3026. DOI:10.1029/JB086iB04p03021 |
| [6] | Haines G V. Spherical cap harmonic analysis. J. Geophys. Res. , 1985, 90(B3): 2583-2591. DOI:10.1029/JB090iB03p02583 |
| [7] | Thebault E, Gaya-Pique L. Applied comparisons between SCHA and R-SCHA regional modeling techniques. Geochem. Geophys. Geosyst. , 2008, 9(1): Q07005. DOI:10.1029/2008GC001953 |
| [8] | Ou J M, Du A M, Thebault E. A high resolution lithospheric magnetic field model over China. Sci. China. Ser. D-Earth Sci., 2012, in press. |
| [9] | Le Mouël J L. Sur la distribution des elements magnétiques en France. Paris: University de Paris, 1969 . |
| [10] | 赵建虎, 王胜平, 刘辉, 等. 海洋局域地磁正常场勒让德多项式模型的建立. 地球物理学进展 , 2008, 23(6): 1802–1808. Zhao J H, Wang S P, Liu H, et al. Construction of normal geomagnetic field model with Legendre polynomial model in local marine area. Progress in Geophysics (in Chinese) , 2008, 23(6): 1802-1808. |
| [11] | 徐文耀. 地磁学. 北京: 地震出版社, 2003 . Xu W Y. Geomagnetism (in Chinese). Beijing: Seismological Press, 2003 . |
| [12] | 安振昌, 谭东海. 1980.0年东亚地区地磁场的勒让德多项式模型. 地球物理学 , 1995, 38(2): 227–233. An Z C, Tan D H. Legendre polynomial models of the geomagnetic field over the eastern Asia at epoch 1980.0. Chinese J. Geophys. (in Chinese) , 1995, 38(2): 227-233. |
| [13] | 张珺, 李星. 利用IDL 进行地学数据处理的多种插值法. 工程地球物理学报 , 2006, 3(1): 49–53. Zhang J, Li X. Several interpolation methods to process geo-science data by IDL. Chinese Journal of Engineering Geophysics (in Chinese) , 2006, 3(1): 49-53. |
| [14] | Briggs I C. Machine contouring using minimum curvature. Geophysics , 1974, 22(3): 48-49. |
| [15] | 徐文耀, 区家明, 杜爱民. 地磁场全球建模与局域建模. 地球物理学进展 , 2011, 26(2): 398–415. Xu W Y, Qu J M, Du A M. Geomagnetic field modelling for the globe and a limited region. Progress in Geophysics (in Chinese) , 2011, 26(2): 398-415. |
| [16] | 徐文耀, 区家明, 杜爱民. 地磁场模型误差分析中的几个问题. 地球物理学进展 , 2011, 26(5): 1485–1509. Xu W Y, Qu J M, Du A M. An analysis of errors in geomagnetic field models. Progress in Geophysics (in Chinese) , 2011, 26(5): 1485-1509. |
| [17] | Xu W Y. Energy budget in the coupling processes of the solar wind, magnetosphere and ionosphere. Chin. J. Space. Sci. , 2011, 31(1): 1-14. |
| [18] | Du A M, Huang Q H, Yang S F. Epicenter location by abnormal ULF electromagnetic emissions. Geophys. Res. Lett. , 2002, 29(10). DOI:10.1029/2001GL013616 |
| [19] | Du A M, Tsurutani B T, Sun W. Anomalous geomagnetic storm of 21-22 January 2005: A storm main phase during northward IMFs. J. Geophys. Res. , 2008, 113(A10): A10214. DOI:10.1029/2008JA013284 |
| [20] | Du A M, Xu W Y, Sun W. Experimental evidence of direct penetration of upstream ULF waves from the solar wind into the magnetosphere during the strong magnetic storm of November 9, 2004. Planet Space Sci. , 2010, 58(7-8): 1040-1044. DOI:10.1016/j.pss.2010.04.001 |
| [21] | Du A M, Tsurutani B T, Sun W. Solar wind energy input during prolonged, intense northward interplanetary magnetic fields: A new coupling function. J. Geophys. Res. , 2011, 116(A12): A12215. DOI:10.1029/2011JA016718 |
| [22] | Du A M, Nakamura R, Zhang T L, et al. Fast tailward flows in the plasma sheet boundary layer during a substorm on 9 March 2008: THEMIS observations. J. Geophys. Res. , 2011, 116(A3): A03216. DOI:10.1029/2010JA015969 |
| [23] | Du A M, Sun W, Tsurutani B T, et al. Observations of dawn-dusk aligned polar cap aurora during the substorms of January 21, 2005. Planet Space Sci. , 2011, 59(13): 1551-1558. DOI:10.1016/j.pss.2011.06.021 |
| [24] | Luo H, Chen G X, Du A M, et al. THEMIS multipoint observations of Pi2 pulsations inside and outside the plasmasphere. J. Geophys. Res , 2011, 116(A12): A12206. DOI:10.1029/2011JA016746 |