根据骆遥等^[1]计算的地磁参考球面上的正常地磁场总强度的全球垂直梯度，其值范围大约为6~34 nT/km，而全球大地水准面起伏大约为-107~85 m^[2]，若忽略大地水准面起伏，粗略估算将致使正常地磁场总强度的计算误差约为-3.6~0.5 nT。但是，以往研究均认为，忽略大地水准面起伏引起的主磁场计算误差应该低于1 nT^[3-4]。尽管该误差远低于IGRF模型本身的误差^[5-6]，但是认识该误差的全球真实分布具有理论意义，对未来高精度的IGRF模型构建也具有参考价值。

1 国际地磁参考场的计算方法

在球面局部直角参考系(X轴指北、Y轴指东、Z轴指向地心)中，地球表面及其以上任一点的主磁场3个分量可以表达为^[7]：

$ {B_x}\left( {r,\theta ,\lambda ,t} \right) = \sum\limits_{n = 1}^{{n_{\max }}} {\sum\limits_{m = 0}^n {{{\left( {\frac{R}{r}} \right)}^{n + 2}}\left[ {g_n^m\left( t \right)\cos m\lambda + h_n^m\left( t \right)\sin m\lambda } \right]\frac{{\rm{d}}}{{{\rm{d}}\theta }}\tilde P_n^m\left( {\cos \theta } \right)} } $

(1)

$ {B_y}\left( {r,\theta ,\lambda ,t} \right) = \sum\limits_{n = 1}^{{n_{\max }}} {\sum\limits_{m = 0}^n {{{\left( {\frac{R}{r}} \right)}^{n + 2}}m\left[ {g_n^m\left( t \right)\sin m\lambda - h_n^m\left( t \right)\cos m\lambda } \right]\frac{{\tilde P_n^m\left( {\cos \theta } \right)}}{{{\rm{sin}}\theta }}} } $

(2)

$ {B_z}\left( {r,\theta ,\lambda ,t} \right) = - \sum\limits_{n = 1}^{{n_{\max }}} {\sum\limits_{m = 0}^n {\left( {n + 1} \right){{\left( {\frac{R}{r}} \right)}^{n + 2}}\left[ {g_n^m\left( t \right)\cos m\lambda + h_n^m\left( t \right)\sin m\lambda } \right]\tilde P_n^m\left( {\cos \theta } \right)} } $

(3)

式中，r、θ与λ分别为地心球坐标系中的地心距、余纬度与经度，R为地磁场参考球面半径(6 371.2 km)，n与m分别为球谐函数的阶与次，n_max为球谐模型的终止阶数，$ \mathit{\tilde P}_\mathit{n}^\mathit{m}$(cosθ)为施密特准归一化的缔合勒让德函数，g_n^m(t)与h_n^m(t)为t时刻(单位a)的球谐系数，d表示微分算子。

对于IGRF模型，由于每5 a更新1次，可以利用线性公式计算t时刻的高斯系数^[7]：

$ g_n^m\left( t \right) = g_n^m\left( {{T_0}} \right) + \left( {t - {T_0}} \right)\dot g_n^m\left( {{T_0}} \right) $

(4)

$ h_n^m\left( t \right) = h_n^m\left( {{T_0}} \right) + \left( {t - {T_0}} \right)\dot h_n^m\left( {{T_0}} \right) $

(5)

式中，T₀为初始年份(对于IGRF-12即为2015-01-01 00：00)，$ \mathit{\dot g}_\mathit{n}^\mathit{m}$(T₀)与$ \mathit{\dot h}_\mathit{n}^\mathit{m}$(T₀)为初始年份给出的未来5 a的长期变化系数(单位nT/a)，即t≤T₀+5。

由于式(1)含有施密特准归一化的缔合勒让德函数在余纬度方向的一阶导数、式(2)中含有奇异项(即1/sinθ)，因此采用无奇异性球谐表达式^[8]进行计算。其中，由于IGRF模型的球谐展开阶数较低，因此施密特准归一化的缔合勒让德函数计算采用Belikov递推法^[9]。

在实际情况中，待计算点的坐标往往采用大地坐标(即大地纬度φ′或余纬度θ′、大地经度λ′与大地高H)，因此需要事先将大地坐标转换为地心球坐标(即地心纬度φ或余纬度θ、地心经度λ与地心距r)，转换公式如下^[10]：

$ \begin{array}{*{20}{c}} {r = \left\{ {H\left[ {H + 2{{\left( {{a^2}{{\sin }^2}\theta ' + {b^2}{{\cos }^2}\theta '} \right)}^{1/2}}} \right] + } \right.}\\ {{{\left. {\frac{{{a^4}{{\sin }^2}\theta ' + {b^4}{{\cos }^2}\theta '}}{{{a^2}{{\sin }^2}\theta ' + {b^2}{{\cos }^2}\theta '}}} \right\}}^{1/2}}} \end{array} $

(6)

$ \theta = \theta ' + \alpha $

(7)

$ \cos \alpha = \frac{{H + {{\left( {{a^2}{{\sin }^2}\theta ' + {b^2}{{\cos }^2}\theta '} \right)}^{1/2}}}}{r} $

(8)

式中，a与b分别为参考椭球的长半轴与短半轴(此处采用WGS-84参考椭球，即a=6 378.137 km，b=6 356.752 km)，α为大地纬度与地心纬度之差。在式(6)~(8)中，若采用海拔高h，则需事先计算大地高H，若忽略垂线偏差则H=h+N，N为大地水准面起伏。

采用式(1)计算得到球面局部直角参考系中的3个分量(B_x、B_y与B_z)后，若忽略垂线偏差，则可以利用参考椭球面局部直角参考系表达地表待计算点处的3个分量(B_X、B_Y与B_Z)，转换公式如下^[11]：

$ {B_X} = {B_x}\cos \alpha + {B_z}\sin \alpha $

(9)

$ {B_Y} = {B_y} $

(10)

$ {B_Z} = {B_z}\cos \alpha - {B_x}\sin \alpha $

(11)

进一步地，其他地磁要素(水平分量B_H、总模量B_F、地磁偏角D与地磁倾角I)可以利用3个分量进行计算，即

$ {B_H} = {\left( {B_X^2 + B_Y^2} \right)^{1/2}} $

(12)

$ {B_F} = {\left( {B_X^2 + B_Y^2 + B_Z^2} \right)^{1/2}} $

(13)

$ D = \arccos \left( {{B_X}/{B_H}} \right) $

(14)

$ I = \arctan \left( {{B_Z}/{B_H}} \right) $

(15)

2 大地水准面起伏对国际地磁参考场计算的间接影响

首先，采用全球陆地表面与海底起伏模拟主磁场解算点的空间位置，所用数据来源于ETOPO1(ice surface version)模型^[12]，点位水平位置采用大地坐标(即大地纬度与大地经度)，分辨率为1′×1′，但是高程与水深相对于平均海平面，若忽略稳态海面地形(在全球范围不会超过±2 m)^[13]，实际上ETOPO1的高程与水深误差约10 m，已经远大于该起伏幅度，因此可以利用大地水准面作为ETOPO1的高程基准面(为了减少计算量，利用球冠滑动平均法^[14]将ETOPO1数据的分辨率降低为7.5′×7.5′)。然后，基于球谐展开阶数达360阶的全球重力场模型EGM96^[2]，采用Rapp^[15]提出的计算方法得到全球大地水准面相对于WGS84参考椭球面的起伏，解算数据间隔为7.5′×7.5′。最后，计算得到考虑与不考虑大地水准面起伏时全球陆地与海底表面主磁场分布之间的差异(图 1)，统计结果见表 1。

由图 1与表 1可以看出，大地水准面起伏对国际地磁参考场计算的间接影响与大地水准面起伏表现出对应性。对于垂向分量、水平分量与总模量，若大地水准面上凸，则致使计算值偏大而呈现正值误差；相反，若大地水准面下凹，则致使计算值偏小而呈现负值误差。

3 结语

本文通过实际计算，得到大地水准面起伏对国际地磁参考场计算间接影响的全球分布。计算结果表明，若忽略大地水准面起伏，可引起地表主磁场的北向、东向、垂向与水平分量、总模量的误差范围分别为-2.37~1.59 nT、-0.51~0.33 nT、-1.95~2.09 nT、-2.37~1.59 nT、-2.39~1.82 nT，误差最大值超出了以往粗略估计的最大误差(例如总模量最大误差的绝对值应该低于1 nT)，而对于磁场倾角和偏角误差影响非常微小。因此，在今后的全球与区域地磁参考场建模以及主磁场解算时，建议尽量采用大地高。