-
月球磁场是月球重要的物理场[1-2]，是研究太阳风与月球相互作用的基础资料[3-4]，也是探索月球内部结构、月球起源与演化历史的重要渠道[5-6]，对月球探测器的轨道控制和着陆器的防护设计等具有重要意义[7-8]。
月球没有类似地球的全球性偶极子磁场，但是在阿波罗计划期间，着陆点附近的磁力仪原位测量数据、子卫星测得的磁场和电子反射率数据以及月岩样品磁学测量数据，均显示月球存在剩余磁化强度[9-10]。为了获得月球全球内源磁场分布，国外学者相继利用月球探测者号(Lunar prospector，LP)获得的电子反射率与三轴磁通门磁力仪测量数据，以及日本Kaguya卫星获得的三轴磁通门磁力仪测量数据进行月球内源磁场建模[11-13]。在所有的磁场建模方法中，球谐分析方法由于其适用于全球性问题，且考虑了磁场的物理性质(例如在场源以外空间满足拉普拉斯方程以及在无穷远处磁位为零等)，而且能够快速实现任意位置磁场预测、分量计算以及频谱分析等处理，已经成为全球磁场建模的标准方法[14-15]。但是，观测点与(或)待解算点往往是不规则分布的，而传统的球谐解算方法计算效率较低并且存在计算不稳定性。为了更高效地进行磁场球谐建模，本文借鉴重力学领域相关方法[16-17]开展月球起伏表面的磁场解算工作。这对于研究月表磁场环境及其与太阳风的相互作用、月球内部磁性结构以及月球起源与演化历史探究也具有重要意义，也为将来进行更高阶/次的月球内源磁场球谐建模提供基础与参考。
-
在月心球坐标系下，令r为观测点至月心的距离，φ与θ分别为纬度与余纬度，λ为经度；在场源即月球物理表面之外，内源磁场的磁位满足拉普拉斯方程，根据边界条件求解该方程，可得到磁位的球谐展开表达式。具体计算如下：
$$ V\left(r, \theta , \lambda \right)=a\sum\limits_{l=0}^{L}\sum\limits_{m=0}^{l}{\left(\frac{a}{r}\right)}^{l+1}\left[{g}_{l}^{m}\mathrm{c}\mathrm{o}\mathrm{s}\right(m\lambda )+ $$ $$ {h}_{l}^{m}\mathrm{s}\mathrm{i}\mathrm{n}\left(m\lambda \right)]{\tilde{P}}_{l}^{m}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right) $$ (1) 式中，a为月球内源磁场参考球面半径；$ {g}_{l}^{m}\mathrm{、}{h}_{l}^{m} $为l阶m次高斯球谐系数；L为截断阶数或最大展开阶数；$ {\tilde{P}}_{l}^{m}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right) $为施密特准归一化的l阶m次缔合勒让德函数。对式(1)求负梯度得到磁场三分量，将其表达在解算点所在的球面局部指北直角坐标系(X轴、Y轴与Z轴分别指向月球北极、东向与径向向内)中，得到：
$$ \left\{\begin{array}{l}{B}_{x}\left(r, \theta , \lambda \right)=\sum\limits_{l=0}^{L}\sum\limits_{m=0}^{l}{\left(\frac{a}{r}\right)}^{l+2}\left[{g}_{l}^{m}\mathrm{c}\mathrm{o}\mathrm{s}\left(m\lambda \right)+{h}_{l}^{m}\mathrm{s}\mathrm{i}\mathrm{n}\left(m\lambda \right)\right]\frac{\mathrm{d}}{\mathrm{d}\theta }{\tilde{P}}_{l}^{m}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right)\\ {B}_{y}\left(r, \theta , \lambda \right)=\sum\limits_{l=0}^{L}\sum\limits_{m=0}^{l}{\left(\frac{a}{r}\right)}^{l+2}m\left[{g}_{l}^{m}\mathrm{s}\mathrm{i}\mathrm{n}\left(m\lambda \right)-{h}_{l}^{m}\mathrm{c}\mathrm{o}\mathrm{s}\left(m\lambda \right)\right]\frac{{\tilde{P}}_{l}^{m}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right)}{\mathrm{s}\mathrm{i}\mathrm{n}\theta }\\ {B}_{z}\left(r, \theta , \lambda \right)=-\sum\limits_{l=0}^{L}\sum\limits_{m=0}^{l}\left(l+1\right){\left(\frac{a}{r}\right)}^{l+2}\left[{g}_{l}^{m}\mathrm{c}\mathrm{o}\mathrm{s}\left(m\lambda \right)+{h}_{l}^{m}\mathrm{s}\mathrm{i}\mathrm{n}\left(m\lambda \right)\right]{\tilde{P}}_{l}^{m}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right)\end{array}\right. $$ (2) 若所有待解算点的月心距离r相等、余纬度θ与经度λ等间隔分布，为了提高解算效率，可将式(2)改写为：
$$ \left\{\begin{array}{l}{B}_{x}\left(r, \theta , \lambda \right)=\sum\limits_{m=0}^{L}\left[{A}_{m}\left(\theta \right)\mathrm{c}\mathrm{o}\mathrm{s}\left(m\lambda \right)+\right.\\ \left.\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }{B}_{m}\left(\theta \right)\mathrm{s}\mathrm{i}\mathrm{n}\left(m\lambda \right)\right]\\ {B}_{y}\left(r, \theta , \lambda \right)=\sum\limits_{m=0}^{L}m\left[{C}_{m}\left(\theta \right)\mathrm{s}\mathrm{i}\mathrm{n}\left(m\lambda \right)-\right.\\ \left.\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }{D}_{m}\left(\theta \right)\mathrm{c}\mathrm{o}\mathrm{s}\left(m\lambda \right)\right]\\ {B}_{z}\left(r, \theta , \lambda \right)=\sum\limits_{m=0}^{L}\left[{E}_{m}\left(\theta \right)\mathrm{c}\mathrm{o}\mathrm{s}\left(m\lambda \right)+\right.\\ \left.\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }{F}_{m}\left(\theta \right)\mathrm{s}\mathrm{i}\mathrm{n}\left(m\lambda \right)\right]\end{array}\right. $$ (3) 式中，
$$ \left\{\begin{array}{c}{A}_{m}\left(\theta \right)=\sum\limits_{l=m}^{L}{\left(\frac{a}{r}\right)}^{l+2}{g}_{l}^{m}\frac{\mathrm{d}}{\mathrm{d}\theta }{\tilde{P}}_{l}^{m}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right)\\ {B}_{m}\left(\theta \right)=\sum\limits_{l=m}^{L}{\left(\frac{a}{r}\right)}^{l+2}{h}_{l}^{m}\frac{\mathrm{d}}{\mathrm{d}\theta }{\tilde{P}}_{l}^{m}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right)\end{array}\right. $$ (4) $$ \left\{\begin{array}{c}{C}_{m}\left(\theta \right)=\sum\limits_{l=m}^{L}{\left(\frac{a}{r}\right)}^{l+2}{g}_{l}^{m}\frac{{\tilde{P}}_{l}^{m}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right)}{\mathrm{s}\mathrm{i}\mathrm{n}\theta }\\ {D}_{m}\left(\theta \right)=\sum\limits_{l=m}^{L}{\left(\frac{a}{r}\right)}^{l+2}{h}_{l}^{m}\frac{{\tilde{P}}_{l}^{m}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right)}{\mathrm{s}\mathrm{i}\mathrm{n}\theta }\end{array}\right. $$ (5) $$ \left\{\begin{array}{c}{E}_{m}\left(\theta \right)=-\sum\limits_{l=m}^{L}\left(l+1\right){\left(\frac{a}{r}\right)}^{l+2}{g}_{l}^{m}{\tilde{P}}_{l}^{m}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right)\\ {F}_{m}\left(\theta \right)=-\sum\limits_{l=m}^{L}\left(l+1\right){\left(\frac{a}{r}\right)}^{l+2}{h}_{l}^{m}{\tilde{P}}_{l}^{m}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right)\end{array}\right. $$ (6) 上述方法称之为集总系数法(lumped coefficients approach，LCA)[18-19]。
在解算起伏月表上的磁场分布时，所有待解算点按照余纬度θ与经度λ等间隔分布，但月心距离r均不相等。
因此，本文选择平均半径r0在月心距离方向进行泰勒级数展开，从而求得相对平均半径起伏观测点上的磁场三分量，即：
$$ \left\{\begin{array}{l}{B}_{x}\left(r, \theta , \lambda \right)=\sum\limits_{k=0}^{K}\frac{1}{k!}{\left(1-\frac{r}{{r}_{0}}\right)}^{k}\sum\limits_{m=0}^{L}\left[{C}_{m}^{k}\left(\theta \right)\mathrm{c}\mathrm{o}\mathrm{s}\left(m\lambda \right)+{D}_{m}^{k}\left(\theta \right)\mathrm{s}\mathrm{i}\mathrm{n}\left(m\lambda \right)\right]\\ {B}_{y}\left(r, \theta , \lambda \right)=\sum\limits_{k=0}^{K}\frac{1}{k!}{\left(1-\frac{r}{{r}_{0}}\right)}^{k}\sum\limits_{m=0}^{L}m\left[{E}_{m}^{k}\left(\theta \right)\mathrm{s}\mathrm{i}\mathrm{n}\left(m\lambda \right)-{F}_{m}^{k}\left(\theta \right)\mathrm{c}\mathrm{o}\mathrm{s}\left(m\lambda \right)\right]\\ {B}_{z}\left(r, \theta , \lambda \right)=\sum\limits_{k=0}^{K}\frac{1}{k!}{\left(1-\frac{r}{{r}_{0}}\right)}^{k}\sum\limits_{m=0}^{L}\left[{G}_{m}^{k}\left(\theta \right)\mathrm{c}\mathrm{o}\mathrm{s}\left(m\lambda \right)+{H}_{m}^{k}\left(\theta \right)\mathrm{s}\mathrm{i}\mathrm{n}\left(m\lambda \right)\right]\end{array}\right. $$ (7) 式中，k表示第k阶径向导数；K表示最大泰勒级数展开次数。此时，令：
$$ {t}_{l}^{k}=\left\{\begin{array}{l}\prod\limits_{i=1}^{k}\left(l+2+i\right), k>0\\ 1, k=0\end{array}\right. $$ 则：
$$ \left\{\begin{array}{c}{A}_{m}^{k}\left(\theta \right)=\sum\limits_{l=m}^{L}{\left(\frac{a}{{r}_{0}}\right)}^{l+2}{g}_{l}^{m}{t}_{l}^{k}\frac{\mathrm{d}}{\mathrm{d}\theta }{\tilde{P}}_{l}^{m}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right)\\ {B}_{m}^{k}\left(\theta \right)=\sum\limits_{l=m}^{L}{\left(\frac{a}{{r}_{0}}\right)}^{l+2}{h}_{l}^{m}{t}_{l}^{k}\frac{\mathrm{d}}{\mathrm{d}\theta }{\tilde{P}}_{l}^{m}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right)\end{array}\right. $$ (8) $$ \left\{\begin{array}{c}{C}_{m}^{k}\left(\theta \right)=\sum\limits_{l=m}^{L}{\left(\frac{a}{{r}_{0}}\right)}^{l+2}{g}_{l}^{m}{t}_{l}^{k}\frac{{\tilde{P}}_{l}^{m}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right)}{\mathrm{s}\mathrm{i}\mathrm{n}\theta }\\ {D}_{m}^{k}\left(\theta \right)=\sum\limits_{l=m}^{L}{\left(\frac{a}{{r}_{0}}\right)}^{l+2}{h}_{l}^{m}{t}_{l}^{k}\frac{{\tilde{P}}_{l}^{m}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right)}{\mathrm{s}\mathrm{i}\mathrm{n}\theta }\end{array}\right. $$ (9) $$ \left\{\begin{array}{c}{E}_{m}^{k}\left(\theta \right)=-\sum\limits_{l=m}^{L}\left(l+1\right){\left(\frac{a}{{r}_{0}}\right)}^{l+2}{g}_{l}^{m}{t}_{l}^{k}{\tilde{P}}_{l}^{m}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right)\\ {F}_{m}^{k}\left(\theta \right)=-\sum\limits_{l=m}^{L}\left(l+1\right){\left(\frac{a}{{r}_{0}}\right)}^{l+2}{h}_{l}^{m}{t}_{l}^{k}{\tilde{P}}_{l}^{m}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right)\end{array}\right. $$ (10) 针对上述表达式中的不稳定项$ \frac{\mathrm{d}}{\mathrm{d}\theta }{\tilde{P}}_{l}^{m}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right) $与$ {\tilde{P}}_{l}^{m}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right)/\mathrm{s}\mathrm{i}\mathrm{n}\theta $，采用无奇异性公式计算[20]：
$$ \left\{\begin{array}{l}\frac{\mathrm{d}}{\mathrm{d}\theta }{\tilde{P}}_{l}^{m}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right)={a}_{l, m}{\tilde{P}}_{l}^{m-1}-{b}_{l, m}{\tilde{P}}_{l}^{m+1}\\ \frac{{\tilde{P}}_{l}^{m}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right)}{\mathrm{s}\mathrm{i}\mathrm{n}\theta }={c}_{l, m}{\tilde{P}}_{l-1}^{m-1}+{d}_{l, m}{\tilde{P}}_{l-1}^{m+1}, m>0\end{array}\right. $$ (11) 式中，
$$ \left\{\begin{array}{l}{a}_{l, m}=0.5\sqrt{l+m}\sqrt{l-m+1}\sqrt{\frac{{C}_{m}}{{C}_{m-1}}}\\ {b}_{l, m}=0.5\sqrt{l+m+1}\sqrt{l-m}\sqrt{\frac{{C}_{m}}{{C}_{m+1}}}\\ {c}_{l, m}=0.5\sqrt{l+m}\sqrt{l+m-1}\sqrt{\frac{{C}_{m}}{{C}_{m-1}}}/m, \mathrm{ }\\ \mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }m>0\\ {d}_{l, m}=0.5\sqrt{l-m}\sqrt{l-m-1}\sqrt{\frac{{C}_{m}}{{C}_{m+1}}}/m, \\ \mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }m>0\end{array}\right. $$ (12) 其中，$ {C}_{m}=\left\{\begin{array}{c}1, m=0\\ 2, m\ne 0\end{array}\right. $，$ {\tilde{P}}_{l}^{-m}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right)=-{\tilde{P}}_{l}^{m}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right) $。
-
根据泰勒级数展开理论，影响计算精度的主要因素在于观测点相对平均半径r0的起伏程度、泰勒级数展开次数。在进行月球起伏表面的磁场解算时，r0应该选取为所有观测点月心距离的平均值，而且观测点偏离平均观测面的程度越剧烈，需要的泰勒级数展开次数就越高。此外，磁场球谐模型的分辨率越高(即最大球谐展开阶数越大)，短波长成分的磁场在径向的变化就越快，因此想要描述短波长磁场的径向快速变化，需要的泰勒级数展开次数则越高。
为了验证上述理论分析的正确性，选取月球地形起伏最为剧烈的背面高地和部分艾特肯盆地作为本文分析的研究区域，如图 1所示。采用美国圣杯号探测器月球轨道激光测高数据构建的720阶次月球全球地形起伏球谐模型[21]，解算了研究区域内的地形起伏，为了测试本文方法的计算效率，设置解算分辨率为0.05°×0.05°。在研究区域内，面积加权的平均月心距离为1 738.244 km，相对该球面，地形起伏的最大值与最小值分别为10 812 m与-7 552 m，可见地形起伏非常剧烈。
图 1 研究区域的月表高程
Figure 1. Lunar Topography of Study Aera
月球内源磁场球谐展开模型最高阶次为450[13]，该模型主要基于日本Kaguya号和美国LP号卫星的三轴磁通门磁力仪测量数据，采用边值法构建获得，其空间分辨率为0.4°，在赤道地区，空间分辨率约为12 km。首先利用式(2)计算得到研究区域内实际月表的磁场，将其作为理论值；然后利用式(7)中的径向泰勒级数展开进行起伏月表的磁场解算，其中r0为月表月心距离的面积加权平均值1 738.244 km。采用不同泰勒级数展开次数范围计算的磁场与理论值之间的差异参见表 1，其中B为磁场总强度。由表 1可知，随着泰勒级数展开次数的增高，计算误差逐渐递减；但在7次展开之后，误差趋于稳定。
表 1 不同最大泰勒级数展开次数计算的磁场与理论值之间差异的统计参数表
Table 1. Statistical Parameters of Differences Between the Theoretical Magnetic Fields and Those Obtained Using Radial Taylor's Series Expansion up to Different Maximum Orders
最大泰勒级数展开次数 磁场分量 误差统计参数/nT 计算时间/s 最小值 最大值 算数平均值 均方根 0次 Bx -167.461 119.944 -0.004 388 7.177 392 312 By -95.129 139.648 -0.002 834 6.365 829 Bz -180.157 195.700 0.000 873 9.564 835 B -195.376 158.133 -2.940 556 12.980 010 1次 Bx -103.487 115.474 -0.001 936 3.293 964 340 By -49.593 75.785 0.003 619 2.702 255 Bz -116.609 144.618 0.005 074 4.255 824 B -79.270 113.639 -1.702 738 4.587 907 2次 Bx -68.192 75.703 0.000 199 1.419 042 378 By -49.300 33.681 -0.000 070 1.045 869 Bz -86.556 74.740 0.000 044 1.763 508 B -22.054 85.126 -0.131 610 2.134 567 3次 Bx -39.906 34.905 0.000 021 0.545 539 420 By -16.926 23.968 0.000 479 0.367 488 Bz -36.923 45.339 0.000 595 0.658 072 B -36.819 44.777 -0.131 752 0.657 029 4次 Bx -13.992 16.604 0.000 047 0.178 899 453 By -9.144 6.558 -0.000 031 0.112 901 Bz -18.476 14.627 0.000 146 0.211 320 B -9.768 13.479 0.017 096 0.204 595 5次 Bx -5.732 4.609 0.000 128 0.061 734 493 By -2.176 3.087 0.000 171 0.042 037 Bz -5.144 6.118 0.000 083 0.074 518 B -4.403 6.046 0.011 798 0.083 024 6次 Bx -1.262 1.756 0.000 089 0.034 876 532 By -0.783 0.557 0.000 102 0.029 811 Bz -1.804 1.381 0.000 065 0.045 669 B -1.811 1.085 0.020 331 0.062 416 7次 Bx -0.658 0.943 0.000 102 0.033 742 575 By -0.697 0.463 0.000 123 0.029 884 Bz -1.037 0.870 0.000 049 0.044 864 B -0.526 1.093 0.019 979 0.063 464 8次 Bx -0.667 0.943 0.000 098 0.033 597 615 By -0.697 0.463 0.000 118 0.029 882 Bz -1.038 0.872 0.000 049 0.044 748 B -0.296 1.094 0.020 294 0.063 427 图 2为10次泰勒级数展开法计算结果与理论值之间的差异分布。由图 2可以看出，误差较大值主要集中于地形下凹的艾特肯盆地区域，而在地形较高的高地区域误差趋于零值。造成该误差的原因在于理论值的计算存在误差，由式(2)可知，当计算点月心距离r小于月球磁场参考球面半径a时，(a/r)l+2是一个不稳定项，且随着球谐阶数l的增大，传统算法可能存在不稳定性。尽管式(7)所示的泰勒级数展开法也存在不稳定项(a/r0)l+2与(1-r/r0)k，但是r0往往选取为平均月心距且与a相当，最大泰勒级数展开次数一般较小(如表 1所示的8次)，因此相比传统算法可以降低计算不稳定性。
图 2 10次泰勒级数展开法计算结果与理论值之间的差异分布
Figure 2. Differences Between the Theoretical Magnetic Fields and Those Obtained Using Radial Taylor's Series Expansion up to 10th Order
分别采用实验区域的平均月心距离1 738.244 km与月球磁场参考球面半径1 737.400 km作为r0，计算磁场误差的均方根，结果如图 3所示。
图 3 磁场计算误差的均方根随泰勒级数展开最大次数的变化
Figure 3. Root Mean Square of Calculating Error of Magnetic Field Intensity Changing with the Maximum Orders of Taylor's Series Expansion
由图 3可以看出，前者的计算收敛速度要高于后者；对于0次与1次泰勒级数展开，两者的计算误差比较一致，其原因在于0次与1次径向梯度主要描述长波长磁场的径向变化，对于r0的敏感度较低，计算误差主要起因于中短波长磁场；当展开次数增高时，逐渐描述中短波长磁场的径向变化，此时对于r0的敏感度较高。
对于本算例，若采用传统算法计算，耗时约为6 h，而采用泰勒级数展开法15次展开耗时也仅约为960 s，可见计算效率得到大幅度提升。根据表 1可以看出，随着泰勒级数展开次数的增大，计算耗时呈线性增加。本文实验在笔记本电脑上完成，其处理器为IntelⓇ XeonⓇ CPU E3-1505M v6 @ 3.00 GHz，内存为32.0 GB。由于执行本文实验采用了式(11)、式(12)的无奇异性算法，若不考虑极区奇异性而采用普通算法，计算耗时还将大幅度降低。为了了解月球全球物理表面与参考球面磁场之间的差异，本文采用720阶次月球全球地形起伏球谐模型[21]，解算得到月球全球地形起伏，如图 4所示，地形最高点(10 768 m)与最低点(-8 424 m)分别位于月球背面的高地与艾特肯盆地。
图 4 月球全球表面高程
Figure 4. Global Topography of the Moon
采用450阶次的月球内源磁场球谐模型，先利用式(7)20次泰勒级数展开法计算得到月球全球实际地形起伏面上的磁场分布，再利用式(4)计算得到月球地形参考球面上的磁场分布，得到两个观测面上磁场之间的差异。相关统计参数见表 2，差异的空间分布如图 5所示。
表 2 月球物理表面与参考球面上的磁场及两者之间差异的统计参数表
Table 2. Statistical Parameters of Magnetic Fields on Lunar Physical Surface and on the Reference Sphere, and Their Differences
解算面与磁场差异 磁场类型 统计参数/nT 最小值 最大值 算数平均值 均方根 月球物理表面 Bx -487.314 370.278 0.020 435 7.035 171 By -284.449 236.321 0.000 489 6.088 395 Bz -424.985 514.243 -0.005 698 9.307 018 B 0.005 590.282 6.413 270 13.159 893 月球参考球面(1 737.152 km) Bx -384.399 317.358 0.018 055 7.265 267 By -332.734 197.771 0.000 000 6.052 672 Bz -388.501 567.346 -0.002 315 9.479 470 B 0.004 575.511 6.368 949 13.389 521 月球物理表面与参考球面的磁场差异 δBx -146.160 163.098 0.002 380 3.145 313 δBy -114.539 154.202 0.000 489 2.552 891 δBz -247.566 183.685 -0.003 383 4.046 276 δB -251.560 155.363 0.044 322 5.500 285 
图 5 月球全球地形起伏面上的磁场分布及其与参考球面上磁场之间的差异分布
Figure 5. Distribution of Global Magnetic Field on Lunar Relief Surface and Difference Distribution Between Magnetic Fields on Lunar Relief Surface and on the Reference Sphere
由表 2可以看出，地形起伏面与参考球面上磁场差异较大，差异幅度达到了参考球面上磁场幅度的一半左右，由于磁场随观测点与场源之间的距离l按照l-3衰减，因此说明月球内源磁场磁性载体深度位置较浅，致使内源磁场在月表附近径向变化剧烈。
由图 5可知，磁场总强度B在实际地形起伏面与参考球面上的差异与地形起伏呈现反相关性，地形越低越靠近场源，磁场越强，反之亦然。但是，磁场径向分量Bz与地形不存在相关性，这是因为月球磁性岩石圈磁性载体的磁化方向并非径向所致。
-
在基于球谐模型解算起伏面上的空间磁场时，传统算法计算效率较低且存在不稳定性。针对该问题，本文提出了径向泰勒级数展开法，通过实验分析，发现本文方法具有计算的高效性与较高的稳定性，具体结论如下：
1) 相比传统算法，本文实验表明径向泰勒级数展开法计算效率提升近20倍，计算时间随级数展开次数增大呈线性增加，而且随着解算点数增多、球谐展开阶次增大，计算效率的提升效果会更加明显。
2) 计算精度主要取决于泰勒级数展开次数，观测面起伏越剧烈、磁场球谐模型展开阶次越高，所需要的泰勒级数展开次数也就越高；泰勒级数展开法中的平均向径建议选取为所有待解算点月心距离的面积加权平均值；计算稳定性主要取决于泰勒级数展开的次数，相比传统算法，本文方法的计算不稳定性得到了大幅度降低。
3) 月表地形起伏面与参考球面上磁场差异较大，说明月球内源磁场磁性载体位置较浅，而且磁场总强度在实际地形起伏面与参考球面上的差异与地形起伏呈现反相关性，而磁场径向分量与地形不存在相关性，说明月球内源磁场场源的磁化方向并非径向。
Fast Spherical Harmonic Synthesis of Magnetic Field at Lunar Surface Based on Radial Taylor's Series Expansion
-
摘要: 在利用球谐模型解算月球起伏表面的磁场时, 传统算法计算效率较低且存在不稳定性, 因此提出一种基于球谐分析的径向泰勒级数展开法。首先采用450阶次的月球内源磁场球谐模型在月表地形起伏剧烈的背面区域进行不同泰勒级数展开次数范围以及不同平均半径的解算精度实验, 验证了所提方法的有效性, 并且通过与传统算法的对比证明了所提方法具有更高的计算效率和计算稳定性; 然后将计算得到的月球全球起伏表面磁场分布与月球参考球面上的磁场分布进行了对比分析, 结果表明, 月表地形起伏面与参考球面上的磁场差异较大, 磁场总强度在实际地形起伏面与参考球面上的差异和地形起伏呈现反相关性, 磁场径向分量与地形起伏不存在相关性, 这些说明月球内源磁场磁性载体埋深较浅且磁化方向并非径向。Abstract:
Objectives For spherical harmonic synthesis of magnetic field at lunar irregular surface, traditional calculating algorithm has a low computing efficiency and the possible instability. Therefore, the radial Taylor's series expansion method is proposed. Methods Firstly, the formulae of the traditional and newly presented methods are given. Then, we adopt a spherical harmonic model of the lunar magnetic field with maximum degree and order of 450, the numerical tests are performed in the lunar highland region with a drastic topographic relief. Finally, the global magnetic fields on lunar relief surface and on the reference sphere are calculated, respectively. Results The results of numerical tests and practical application show that the calculating accuracy of our method mainly depends on the maximum order of the Taylor's series expansion, and if the observing surface changes more dramatically and the spherical harmonic model of the magnetic field has higher maximum degree and order, a higher maximum order of the Taylor's series expansion is required. Moreover, the average radius of the algorithm is suggested to be the area-weighted mean of the radius values of all calculating points. The computing time varies linearly with the maximum order of the Taylor's series expansion. If more calculating points and higher maximum degree and order of spherical harmonic model of the magnetic field, the computational efficiency will be improved more significantly by our proposed method. The differences between magnetic fields on lunar relief surface and on the reference sphere are very large, which suggests that the source depth of the internal magnetic field of the Moon is relatively shallow. Especially, the differences of the intensity of the magnetic field have a negative correlation with the topography variations, but there is no correlation between the radial component of the magnetic field and the topography variations. These indicate that the sources' magnetized directions of the internal magnetic field of the Moon are not radial. Conclusions In brief, this study significantly certifies that our proposed method has a high computing efficiency and a very weakly instability. -
图 2 10次泰勒级数展开法计算结果与理论值之间的差异分布
Figure 2. Differences Between the Theoretical Magnetic Fields and Those Obtained Using Radial Taylor's Series Expansion up to 10th Order
图 3 磁场计算误差的均方根随泰勒级数展开最大次数的变化
Figure 3. Root Mean Square of Calculating Error of Magnetic Field Intensity Changing with the Maximum Orders of Taylor's Series Expansion
图 5 月球全球地形起伏面上的磁场分布及其与参考球面上磁场之间的差异分布
Figure 5. Distribution of Global Magnetic Field on Lunar Relief Surface and Difference Distribution Between Magnetic Fields on Lunar Relief Surface and on the Reference Sphere
表 1 不同最大泰勒级数展开次数计算的磁场与理论值之间差异的统计参数表
Table 1. Statistical Parameters of Differences Between the Theoretical Magnetic Fields and Those Obtained Using Radial Taylor's Series Expansion up to Different Maximum Orders
最大泰勒级数展开次数 磁场分量 误差统计参数/nT 计算时间/s 最小值 最大值 算数平均值 均方根 0次 Bx -167.461 119.944 -0.004 388 7.177 392 312 By -95.129 139.648 -0.002 834 6.365 829 Bz -180.157 195.700 0.000 873 9.564 835 B -195.376 158.133 -2.940 556 12.980 010 1次 Bx -103.487 115.474 -0.001 936 3.293 964 340 By -49.593 75.785 0.003 619 2.702 255 Bz -116.609 144.618 0.005 074 4.255 824 B -79.270 113.639 -1.702 738 4.587 907 2次 Bx -68.192 75.703 0.000 199 1.419 042 378 By -49.300 33.681 -0.000 070 1.045 869 Bz -86.556 74.740 0.000 044 1.763 508 B -22.054 85.126 -0.131 610 2.134 567 3次 Bx -39.906 34.905 0.000 021 0.545 539 420 By -16.926 23.968 0.000 479 0.367 488 Bz -36.923 45.339 0.000 595 0.658 072 B -36.819 44.777 -0.131 752 0.657 029 4次 Bx -13.992 16.604 0.000 047 0.178 899 453 By -9.144 6.558 -0.000 031 0.112 901 Bz -18.476 14.627 0.000 146 0.211 320 B -9.768 13.479 0.017 096 0.204 595 5次 Bx -5.732 4.609 0.000 128 0.061 734 493 By -2.176 3.087 0.000 171 0.042 037 Bz -5.144 6.118 0.000 083 0.074 518 B -4.403 6.046 0.011 798 0.083 024 6次 Bx -1.262 1.756 0.000 089 0.034 876 532 By -0.783 0.557 0.000 102 0.029 811 Bz -1.804 1.381 0.000 065 0.045 669 B -1.811 1.085 0.020 331 0.062 416 7次 Bx -0.658 0.943 0.000 102 0.033 742 575 By -0.697 0.463 0.000 123 0.029 884 Bz -1.037 0.870 0.000 049 0.044 864 B -0.526 1.093 0.019 979 0.063 464 8次 Bx -0.667 0.943 0.000 098 0.033 597 615 By -0.697 0.463 0.000 118 0.029 882 Bz -1.038 0.872 0.000 049 0.044 748 B -0.296 1.094 0.020 294 0.063 427 
下载: 导出CSV
表 2 月球物理表面与参考球面上的磁场及两者之间差异的统计参数表
Table 2. Statistical Parameters of Magnetic Fields on Lunar Physical Surface and on the Reference Sphere, and Their Differences
解算面与磁场差异 磁场类型 统计参数/nT 最小值 最大值 算数平均值 均方根 月球物理表面 Bx -487.314 370.278 0.020 435 7.035 171 By -284.449 236.321 0.000 489 6.088 395 Bz -424.985 514.243 -0.005 698 9.307 018 B 0.005 590.282 6.413 270 13.159 893 月球参考球面(1 737.152 km) Bx -384.399 317.358 0.018 055 7.265 267 By -332.734 197.771 0.000 000 6.052 672 Bz -388.501 567.346 -0.002 315 9.479 470 B 0.004 575.511 6.368 949 13.389 521 月球物理表面与参考球面的磁场差异 δBx -146.160 163.098 0.002 380 3.145 313 δBy -114.539 154.202 0.000 489 2.552 891 δBz -247.566 183.685 -0.003 383 4.046 276 δB -251.560 155.363 0.044 322 5.500 285
下载: 导出CSV
-
[1] 李斐, 鄢建国. 月球重力场的确定及构建我国自主月球重力场模型的方案研究[J]. 武汉大学学报·信息科学版, 2007, 32(1): 6-10 http://ch.whu.edu.cn/article/id/1811 Li Fei, Yan Jianguo. Principle and Method of Lunar Gravity Field Fetermination and Project on Self-determinational Lunar Gravity Field[J]. Geomatics and Information Science of Wuhan University, 2007, 32(1): 6-10 http://ch.whu.edu.cn/article/id/1811 [2] 钟振, 李斐, 鄢建国, 等. 新近月球重力场模型的比较与分析[J]. 武汉大学学报·信息科学版, 2013, 38(4): 390-393 http://ch.whu.edu.cn/article/id/755 Zhong Zhen, Li Fei, Yan Jianguo, et al. Comparison and Analysis on Main and Newly Lunar Gravity Field Models[J]. Geomatics and Information Science of Wuhan University, 2013, 38(4): 390-393 http://ch.whu.edu.cn/article/id/755 [3] Wang X Q, Cui J, Wang X D, et al. The Solar Wind Interactions with Lunar Magnetic Anomalies: A Case Study of the Chang'E-2 Plasma Data near the Serenitatis Antipode[J]. Advances in Space Research, 2012, 50(12): 1 600-1 606 doi: 10.1016/j.asr.2011.12.003 [4] Deca J, Divin A, Lue C, et al. Reiner Gamma Albedo Features Reproduced by Modeling Solar Wind Standoff[J]. Communications Physics, 2018, 12(1): 1-8 http://www.nature.com/articles/s42005-018-0012-9 [5] Wieczorek M A, Weiss B P, Stewart S T. An Impactor Origin for Lunar Magnetic Anomalies[J]. Science, 2012, 335(6 073): 1 212-1 215 http://search.ebscohost.com/login.aspx?direct=true&db=aph&AN=74262035&site=ehost-live [6] Takahashi F, Tsunakawa H, Shimizu H, et al. Reorientation of the Early Lunar Pole[J]. Nature Geoscience, 2014, 7(6): 409-412 doi: 10.1038/ngeo2150 [7] 潘永信, 纪新林, 朱日祥. 月球磁学观测与研究进展[J]. 地球化学, 2010, 39(1): 32-36 https://www.cnki.com.cn/Article/CJFDTOTAL-DQHX201001008.htm Pan Yongxin, Ji Xinlin, Zhu Rixiang. A Review of Lunar Magnetism[J]. Geochimica, 2010, 39(1): 32-36 https://www.cnki.com.cn/Article/CJFDTOTAL-DQHX201001008.htm [8] 李泳泉, 刘建忠, 欧阳自远, 等. 月球磁场与月球演化[J]. 地球物理学进展, 2005, 20(4): 1 003-1 008 https://www.cnki.com.cn/Article/CJFDTOTAL-DQWJ200504019.htm Li Yongquan, Liu Jianzhong, Ouyang Ziyuan, et al. Lunar Magnetism and Its Evolution[J]. Progress in Geophysics, 2005, 20(4): 1 003-1 008 https://www.cnki.com.cn/Article/CJFDTOTAL-DQWJ200504019.htm [9] 张昌达, 陈超. 月球引力场和磁场探测新进展[J]. 地质科技情报, 2004, 23(4): 6-11 doi: 10.3969/j.issn.1000-7849.2004.04.002 Zhang Changda, Chen Chao. Progress in Prospecting for Lunar Gravity and Magnetic Fields[J]. Geological Science and Technology Information, 2004, 23(4): 6-11 doi: 10.3969/j.issn.1000-7849.2004.04.002 [10] Fuller M, Cisowski S M. Lunar Paleomagnetism[M]. San Diego: Academic Press, 1987 [11] Hakekas J S, Mitchell D L, Lin R P, et al. Mapping of Crustal Magnetic Anomalies on the Lunar near Side by the Lunar Prospector Electron Reflectometer [J]. Journal of Geophysical Research: Planets, 2001, 106(E11): 27 841-27 852 doi: 10.1029/2000JE001380 [12] Purucker M E, Nicholas J B. Global Spherical Harmonic Models of the Internal Magnetic Field of the Moon Based on Sequential and Coestimation Approaches[J]. Journal of Geophysical Research: Planets, 2010, 115(E12): 27 932-27 943 doi: 10.1029/2010JE003650/full [13] Tsunakawa H, Takahashi F, Schimizu H, et al. Surface Vector Mapping of Magnetic Anomalies over the Moon Using Kaguya and Lunar Prospector Observations [J]. Journal of Geophysical Research: Planets, 2015, 120(6): 1 160-1 185 doi: 10.1002/2014JE004785 [14] 徐文耀, 区加明, 杜爱民. 地磁场全球建模与区域建模[J]. 地球物理学进展, 2011, 26(2): 398-415 doi: 10.3969/j.issn.1004-2903.2011.02.002 Xu Wenyao, Ou Jiaming, Du Aimin. Geomagnetic Field Modelling for the Globe and a Limited Region[J]. Progress in Geophysics, 2011, 26(2): 398-415 doi: 10.3969/j.issn.1004-2903.2011.02.002 [15] 黄朝艳, 刘浩, 苏斌, 等. 月球磁场建模研究[J]. 测绘工程, 2013, 22(6): 77-80 doi: 10.3969/j.issn.1006-7949.2013.06.020 Huang Chaoyan, Liu Hao, Su Bin, et al. A Study of Modeling of Lunar Magnetic Field[J]. Engineering of Surveying and Mapping, 2013, 22(6): 77-80 doi: 10.3969/j.issn.1006-7949.2013.06.020 [16] Hirt C. Efficient and Accurate High-degree Spherical Harmonic Synthesis of Gravity Field Functionals at the Earth's Surface Using the Gradient Approach[J]. Journal of Geodesy, 2012, 86(9): 729-744 doi: 10.1007/s00190-012-0550-y [17] Bucha B, Janák J. A MATLAB-Based Graphical User Interface Program for Computing Functionals of the Geopotential up to Ultra-high Degrees and Orders: Efficient Computation at Irregular Surfaces[J]. Computers & Geosciences, 2014, 66: 219–227 http://www.sciencedirect.com/science/article/pii/S0098300414000302 [18] Colombo O L. Numerical Methods for Harmonic Analysis on the Sphere[R]. Department of Geodetic Science and Surveying, The Ohio State University, Columbus, Ohio, 1981 [19] 钟波. 基于GOCE卫星重力测量技术确定地球重力场的研究[D]. 武汉: 武汉大学, 2010 Zhong Bo. Study on the Determination of the Earth's Gravity Field from Satellite Gravimetry Mission GOCE[D]. Wuhan: Wuhan University, 2010 [20] Du J, Chen C, Lesur V, et al. Non-singular Spherical Harmonic Expressions of Geomagnetic Vector and Gradient Tensor Fields in the Local North-Oriented Reference Frame[J]. Geoscientific Model Development, 2015, 8(7): 1 979–1 990 doi: 10.5194/gmd-8-1979-2015 [21] Smith D E, Zuber M T, Neumann G A, et al. Initial Observations from the Lunar Orbiter Laser Altimeter (LOLA) [J]. Geophysical Research Letters, 2010, 37: L18204 -
点击查看大图
图(5) / 表(2)
计量
- 文章访问数: 476
- HTML全文浏览量: 116
- PDF下载量: 26
- 被引次数: 0
