在火星的两颗天然卫星中，火卫一体积较大，距离火星距离较近。它的物理特性、分布、形成和演化是了解火星及其自身内部结构、物质组成、早期演化和生命形成的关键^[1]。火卫一微弱的引力场使得探测器在火卫一上降落和返回地球相比其他小行星更为容易^[2]，在未来，火卫一可以作为人类踏入深空的桥梁与补给站。

火卫一最早由Hall于1877年发现。美国国家航空航天局(National Aeronautics and Space Administration, NASA)的Mariner 9探测器对火卫一进行了近距离拍摄；Viking探测器任务执行过程中利用影像数据得到了火卫一的形状模型^[3]；火星全球勘探者(Mars Global Surveyor, MGS)以太阳同步轨道环绕火星，通过星载相机、热辐射光谱仪对火卫一进行了观测^[4]；火星探路者(Mars Pathfinder)使用自身携带的相机在火星表面上对火卫一进行了拍摄^[5]；欧洲航天局的火星快车任务(Mars Express，MEX)提供地球与其他国家部署的登陆车之间的通信中转服务，并由此成为国际火星探索工作的枢纽部分，同时MEX在飞掠火卫一期间对火卫一开展了形状和重力场方面的大量观测^[6-7]。

截至目前，由于没有专门的探测器对火卫一进行伴飞或绕飞探测，所以短时间内无法利用轨道跟踪数据反演得到重力场模型。近年来，国外针对火卫探测的MMX^[8](Martian Moons Exploration)任务和DePhine^[9](The Deimos and Phobos Interior Explorer)计划相继推出。其中，MMX任务计划当航天器到达火星后，轨道机动将航天器送到火卫一周围的准卫星轨道，绕火卫一进行一定时间的伴飞并进行各种观测。MMX任务的目的之一即是解算火卫一高阶重力场，并开展火卫一采样返回研究。在这些探测任务中，探测器的轨道设计和在星体表面着陆等任务的成功与否很大程度上取决于重力场模型的精确度，因此，火卫一重力场解算任务的作用至关重要。目前关于火卫一重力场模型的研究大多局限在低阶次。Davis等^[10]于1981年计算了均值密度假设下基于三轴椭球的表面格网形式的火卫一的重力场；Duxbury^[3]通过形状模型球谐系数和重力场球谐系数间的关系给出火卫一6阶次的球谐系数及主惯性矩，得到的3轴主惯性矩A、B、C分别为0.336 2、0.387 1、0.477 3；Martinec等^[11]根据Duxbury的球谐形式的形状模型，在均值密度假设下计算了6阶次的火卫一重力场模型；Thomas^[12]假设均值密度，根据表面控制点数据计算得到表面格网形式的火卫一重力场；Le Maistre等^[13]利用均值密度假设通过立方体叠加的方式逼近火卫一形状，计算每个立方体的球谐系数并相加得到火卫一的球谐系数，得到10阶次的火卫一重力场，并计算得到主惯性矩A、B、C分别为0.353 657、0.418 499、0.490 239。史弦等^[14]通过解析法获得17阶次的火卫一重力场，利用数值法获得火卫一表面重力的格网化模型。

关于火卫一重力场的研究，阶次最高的为史弦^[1]发表的外部14 km球面、阶次至20阶的重力场模型，其所使用的形状模型为参考半径为10.955 8 km的球谐展开形式，且是基于球谐系数的形状模型计算所得的引力位，形状模型所在坐标系为行星体固坐标系，其坐标系原点为行星历表积分中心，由于历表存在误差，因此这一坐标系统的原点既不与火卫一质心重合，也不与形状中心重合，并且坐标轴与行星主惯性轴不重合。这导致在进行探测器着陆、定轨等研究时，重力场模型需要转换到质心为原点、坐标轴和主惯性轴重合的主惯性轴坐标系中，较为不便。

针对当前研究中的不足，本文介绍了火卫一质心和惯性张量的计算思路，并给出火卫一从此行星体固坐标系转换到主惯性轴坐标系下的转换方法。采用64 800个特征点、49 152个三角形面组成的火卫一形状模型为基础数据，在计算中将参考系由此行星体固坐标系转换为主惯性坐标系，可以提供一个相对通用的重力场模型。假设其内部密度均匀，通过基于形状模型正演的方法，利用边长、面积、体积3种方式进行外部引力位的求解。3种方法实质上均为将天体划分成微元，此时天体即为微元的集合。以此为基础计算了火卫一外部14 km球面上1°为间隔的格网点的引力位，并反演求解了阶次至80阶的火卫一重力场模型。通过计算不同阶次下加速度的差值，以验证高阶项重力场阶次对加速度的影响。

1 计算模型与方法 1.1 转动惯量求取

因为火卫一的形状模型所使用的参考系统是行星体固坐标系，由国际天文学联合会（International Astronomical Union，IAU）给出的定义为x轴指向火星，z轴和平均旋转轴重合，经度是从本初子午线向西测量的。此坐标系为右手坐标系，利用这一准则确定y轴指向。因为低阶球谐系数决定了坐标系的原点和坐标轴的定向问题，所以为将形状模型所用参考坐标系和惯性系统一，需对其进行旋转平移，其中涉及到低阶球谐系数的计算。

0阶次至2阶次的球谐系数的求取公式如下：

$ \left\{ \begin{array}{l} {C_{lm}} = \smallint \smallint \smallint {R^n}{P_{{\rm{ l0}}}}(cos\theta ){\rm{d}}M,m = {\rm{ }}0\\ {C_{lm}} = 2\frac{{(n - m)!}}{{(n + m)!}}\smallint \smallint \smallint {R^n}{A_{lm}}(\theta ',\lambda '){\rm{d}}M,m \ne {\rm{ }}0\\ {S_{lm}} = 2\frac{{(n - m)!}}{{(n + m)!}}\smallint \smallint \smallint {R^n}{B_{lm}}(\theta ',\lambda '){\rm{d}}M,m \ne {\rm{ }}0 \end{array} \right. $

(1)

式中，A_lm=P_lm(cosθ)·cos(mλ)；B_lm=P_lm(cosθ)·sin(mλ)。

阶数为0时，球谐系数为：A₀₀=1，B₀₀=0，代入式(1)中可得C₀₀=∭dM=M，即零阶项系数为火卫一的质量，且无论参考系如何变换，都不会影响此项系数。

阶数为1时，球谐系数为：C₁₀=∭z'dM；C₁₁=∭x'dM；S₁₁=∭y'dM；且火卫一质心坐标为：${x_{{\text{origin}}}} = {\text{ }}\frac{1}{M}\smallint \smallint \smallint x'{\text{d}}M$；${y_{{\text{origin}}}} = {\text{ }}\frac{1}{M}\smallint \smallint \smallint y'{\text{d}}M$；${z_{{\text{origin}}}} = {\text{ }}\frac{1}{M}\smallint \smallint \smallint z'{\text{d}}M$。由此可得：C₁₀=M·x_origin；C₁₁=M·y_origin；S₁₁=M·z_origin。

阶数为2时，球谐系数为：${C_{{\text{20}}}} = {\text{ }}\frac{1}{2}\smallint \smallint \smallint \left( { - x{'^2} - y{'^2} + 2z{'^2}} \right){\text{d}}M$；${C_{{\text{21}}}} = {\text{ }}\frac{1}{2}\smallint \smallint \smallint \left( {x'z'} \right){\text{d}}M$；${S_{{\text{21}}}} = {\text{ }}\frac{1}{2}\smallint \smallint \smallint \left( {y'z'} \right){\text{d}}M$；${C_{{\text{22}}}} = {\text{ }}\frac{1}{4}\smallint \smallint \smallint \left( {x{'^2} - y{'^2}} \right){\text{d}}M$；${S_{{\text{22}}}} = {\text{ }}\frac{1}{2}\smallint \smallint \smallint \left( {x'y'} \right){\text{d}}M$。

对轴的惯性矩为：I_xx=∭(y'²+z'²) dM；I_yy=∭(x'²+z'²) dM；I_zz=∭(x'²+y'²) dM；I_xy=∭(x'⋅y') dM；I_xz=∭(x'⋅z') dM；I_yz=∭(y'⋅z') dM。

当坐标轴和惯性轴重合时，C₂₁作为关于oxz平面的惯量积，S₂₁作为关于oyz平面的惯量积，S₂₂作为关于oxy平面的惯量积则为零，所以可得：${C_{{\text{20}}}} = {\text{ }}\frac{{{I_{xx}} + {I_{yy}}}}{2} - {I_{zz}}$；S₂₀=0；${C_{{\text{21}}}} = {\text{ }}\frac{1}{2}{I_{xz}} = 0$；${S_{{\text{21}}}} = {\text{ }}\frac{1}{2}{I_{yz}} = 0$；${C_{{\text{22}}}} = {\text{ }}\frac{1}{4}\left( {{I_{yy}} - {I_{xx}}} \right)$；${S_{{\text{22}}}} = {\text{ }}\frac{1}{2}{I_{xy}} = 0$。

为更明确地排列，主惯性矩按照I_xx＜I_yy＜I_zz的方式进行排列，之后再进行引力位的计算。

1.2 质心的求取

质心坐标的求取公式如下：

$ {\text{Coo}}{{\text{r}}_{{\text{mass}}}} = \mathop \sum \limits_{{\text{faces}}} \Delta M \cdot {\text{ }}\Delta R/M = \mathop \sum \limits_{{\text{faces}}} \Delta V \cdot {\text{ }}\Delta R/V $

(2)

式中，Coor_mass为质心的坐标; ΔR=(D+E+F)/4；ΔV为以质心为顶点、三角形为底面的三棱锥的体积；V为所有的三棱锥的体积；D、E、F分别为三角形的3个顶点的坐标(逆时针方向取点)。

校正形状模型到质心坐标系的计算公式为：

式中，Coor为形状模型的坐标。

1.3 惯性张量的求取

火卫一坐标系间旋转矩阵的解算涉及到惯性张量的计算。转动惯量通过惯性张量矩阵来表述，且此矩阵为对称矩阵。惯性张量的表达形式为：

$ \begin{array}{*{20}{l}} {I = } \end{array}\left[ {\begin{array}{*{20}{c}} {{I_{xx}}}&{{I_{xy}}}&{{I_{xz}}}\\ {{I_{xy}}}&{{I_{yy}}}&{{I_{yz}}}\\ {{I_{xz}}}&{{I_{yz}}}&{{I_{zz}}} \end{array}} \right] $

(3)

式中，I_xx=P_yy+P_zz, I_yz=−P_yz，I_yy=P_xx+P_zz, I_xz=−P_xz，I_zz=P_xx+P_yy, I_xy=−P_xy。则${{P_{jk}} = \sum\limits_{{\rm{faces}}} \Delta {P_{jk}},\Delta {P_{jk}} = \frac{{\rho \Delta V}}{{20}}[2{D_j}{D_k} + 2{E_j}{E_k} + 2{F_j}{F_k} + {D_j}{E_k} + {D_k}{E_j} + {D_j}{F_k} + {D_k}{F_j} + {E_j}{F_k} + {E_k}{F_j}]} $。

得到惯性张量矩阵I后, 将I对角化，求取矩阵I的3个特征值和相对应的3个特征向量，则对角元素即为主惯性矩，将主惯性矩按I_xx＜I_yy＜I_zz的顺序排列，将其对应的特征向量同样按照此顺序排列，可得到旋转矩阵M_rotate：

$ {\text{Coo}}{{\text{r}}_{校正}}{\text{ = Coo}}{{\text{r}}_{质心校正}}{\text{·}}{M_{{\text{rotate}}}} $

(4)

经过校正以后，可得到以质心为原点、惯性轴与坐标轴重合的火卫一形状模型。

1.4 多面体模型引力位

天体的重力场可根据射电科学数据，通过分析航天器在场内的运动进行测量。迄今为止，这种重力场模型只适用于有限数量的太阳系天体。对于火卫一来说，由于近距离飞行的短弧数据和不利的观测几何条件，利用MEX射电跟踪数据仅能计算低阶重力场系数^[15]。对后续火卫一探测任务而言，基于形状模型所计算的高阶重力场，可以作为任务仿真分析中的先验模型，同时对导航也能起到一定的支撑作用。基于形状模型研究类似火卫一这样不规则的小天体重力场时，通常密度假设为均质密度^[16]。

多面体是将小行星表面分割成三角形去逼近其形状，火卫一表面形状由大量三角形组成，每个三角形由边组成。通过对其进行积分变换可以进行引力场的计算。Werner等^[17]于1996年提出将不规则形状小行星表面剖分为三角形，以多面体来逼近小行星的实际形状进而求取引力位的方法。本文通过假设火卫一的密度为均质密度，由地形点剖分所得多面体逼近实际天体形状，从而求取其外部固定半径球面上的引力位。引力位可以通过边长法、面积法、体积法3种方式进行积分计算得到。其中，边长法积分单元为三角形边长，面积法积分单元为三角形面积，体积法积分单元仍为三角形面积，但通过高斯散度定理将引力位的三维积分转换成二维积分时，新增的${{\hat r}_f} \cdot {{\hat n}_f}$项含义为所计算火卫一外部一点到面积所在平面的垂直距离。其中，${{\hat r}_f}$为场点到多面体面上的单位矢量；${{\hat n}_f}$为面上单位法向量。

3种方法中密度也不相同，需要单独计算。利用边长法所计算引力位为：

$ \left\{ {_{{\rho _{{\text{edge}}}} = M/(\mathop \sum \limits_{f \in {\text{ faces}}} \mathop \sum \limits_{e \in {\text{ edges}}} {\text{length}}_e^f)}^{{U_{{\text{edge}}}} = \frac{1}{2}G{\rho _{{\text{edge}}}}\mathop \sum \limits_{f \in {\text{ faces}}} (\mathop \sum \limits_{e \in {\text{ edges}}} L_e^f)}} \right\} $

(5)

利用面积法所计算引力位为：

$ \left\{ {\begin{array}{*{20}{l}} {{U_{{\rm{surface}}}} = G{\rho _{{\rm{surface}}}}\sum\limits_{f \in {\rm{faces}}} {\sum\limits_{e \in {\rm{edges}}} {\hat n_e^f} \cdot L_e^f} - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} G{\rho _{{\rm{surface}}}}\sum\limits_{f \in {\rm{faces}}} {{n_f} \cdot {r_f} \cdot {\omega _f}} }\\ {{\rho _{{\rm{surface}}}} = M/(\sum\limits_{f \in {\rm{faces}}} {{\rm{are}}} {{\rm{a}}_f})} \end{array}} \right. $

(6)

利用体积法所计算引力位为：

$ \left\{ {\begin{array}{*{20}{l}} {{U_{{\rm{polyhedron}}}}{\rm{ = }}\frac{1}{2}G{\rho _{{\rm{polyhedron}}}} \cdot {{\hat r}_f} \cdot {{\hat n}_f} \cdot }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} (\sum\limits_{f \in {\rm{faces}}} {(\sum\limits_{e \in {\rm{edges}}} {\hat n_e^f} \cdot L_e^f)} - \sum\limits_{f \in {\rm{faces}}} {{n_f} \cdot {r_f} \cdot {\omega _f}} )}\\ {{\rho _{{\rm{polyhedron}}}} = M/{\rm{volum}}{{\rm{e}}_{{\rm{polyhedron}}}}} \end{array}} \right. $

(7)

式（5）~（7）中，G为引力常量；ρ为火卫一体密度；${\hat n_e^f}$为组成面的边长的法向量；n_f为面的法向量；r_f为场点到多面体面上三角形3个顶点的矢量；${L_e^f}$为组成面的边计算所得引力位，ω_f为立体角，两者的计算式如下：

$ {L_e^f = {\rm{ln }}\frac{{{r_i} + {r_j} + {e_{ij}}}}{{{r_i} + {r_j} - {e_{ij}}}},{\omega _f} = {\rm{ }}2{\rm{ arctan }}\left( {\frac{{{r_1} \cdot {r_2} \times {r_3}}}{{{r_1}{r_2}{r_3} + {r_1}({r_2} \cdot {r_3}) + {r_2}({r_3} \cdot {r_1}) + {r_3}({r_1} \cdot {r_2})}}} \right)} $

(8)

1.5 球谐系数求取

由于引力场模型在布里渊球内收敛情况十分复杂，所以本文建立一个半径为14 km的完全将火卫一包附在其中的布里渊球，以1°为间隔划分格网，计算格网点上的引力位。对于特定的坐标系，惯性张量是唯一确定的值，而惯性张量的求取和多面体体积相关，所以体积法计算的引力位应当更为精准。将由上述体积法计算所得的引力位作为观测量，利用最小二乘法进行高阶球谐系数的求解，计算式为：

$ V = \frac{{GM}}{r}\mathop \Sigma \limits_{l = 0}^\infty \mathop \Sigma \limits_{m = 0}^l {\left( {\frac{R}{r}} \right)^n}{P_{lm}}\left( {\cos \theta } \right)\left[ {{C_{lm}} \cdot \cos \left( {m\lambda } \right) + {S_{lm}} \cdot \sin \left( {m\lambda } \right)} \right] $

(9)

式中，R为火卫一外部点距火卫一质量中心的距离；(λ, θ, r)为场点坐标(θ为纬度)；P_lm(cosθ)为勒让德多项式；C_lm、S_lm为球谐系数；l为阶数；m为次数。则有：

$ \left\{ {\begin{array}{*{20}{l}} {{P_{00}}{C_{00}} + \frac{R}{r}{P_{10}}{C_{10}} + \frac{R}{r}{P_{11}}\left( {{C_{11}}\cos \left( {m{\lambda _1}} \right) + {S_{11}}\sin \left( {m{\lambda _1}} \right)} \right) + \cdots + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{\left( {\frac{R}{r}} \right)}^n}{P_{nn}}\left( {{C_{nn}}\cos \left( {m{\lambda _1}} \right) + {S_{nn}}\sin \left( {m{\lambda _1}} \right)} \right) = \frac{{{V_1} \cdot r}}{{G \cdot M}}}\\ {{P_{00}}{C_{00}} + \frac{R}{r}{P_{10}}{C_{10}} + \frac{R}{r}{P_{11}}\left( {{C_{11}}\cos \left( {m{\lambda _2}} \right) + {S_{11}}\sin \left( {m{\lambda _2}} \right)} \right) + \cdots + }\\ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\left( {\frac{R}{r}} \right)^n}{P_{nn}}\left( {{C_{nn}}\cos \left( {m{\lambda _2}} \right) + {S_{nn}}\sin \left( {m{\lambda _2}} \right)} \right) = \frac{{{V_2} \cdot r}}{{G \cdot M}}\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \vdots \end{array}\\ {{P_{00}}{C_{00}} + \frac{R}{r}{P_{10}}{C_{10}} + \frac{R}{r}{P_{11}}\left( {{C_{11}}\cos \left( {m{\lambda _i}} \right) + {S_{11}}\sin \left( {m{\lambda _i}} \right)} \right) + \cdots + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{\left( {\frac{R}{r}} \right)}^n}{P_{nn}}\left( {{C_{nn}}\cos \left( {m{\lambda _i}} \right) + {S_{nn}}\sin \left( {m{\lambda _i}} \right)} \right) = \frac{{{V_i} \cdot r}}{{G \cdot M}}} \end{array}} \right. $

(10)

由于地图投影到平面上会发生变形，每条经线上纬度变化1°对应实际距离不相同。同一经度不同纬度处对应实际距离为：

$ \Delta l = {\text{ }}\frac{{180}}{\pi } \cdot {R_{{\text{equator}}}} \cdot {\text{ cos}}({\text{latitude}}) $

(11)

所以在利用最小二乘法进行计算时，将式(9)写为Ax=b的形式，需要给定权重W=cos(latitude)，即：

$ {A^{\text{T}}}AWx = {A^{\text{T}}}Wb \Rightarrow x = {({A^{\text{T}}}AW)^{ - 1}}{A^{\text{T}}}Wb $

(12)

由此可解得球谐系数C_lm和S_lm。

2 计算与分析 2.1 数据说明

本文中所使用的形状模型来自PDS(Planetary Data System)网站（https://nssdc.gsfc.nasa.gov/planetary/pds.html）, 所选取火卫一形状模型包含49 152个三角形面数，共计64 800个数据点。形状模型所在坐标系已由火卫一体固坐标系转换为火卫一质心为坐标原点、主惯性轴为坐标轴的坐标系。图 1给出了利用多面体模型逼近的火卫一形状。火卫一的物理与几何参数为：密度为1.876×10¹² kg/km^3[18], GM=0.712 7×10^-3 km³/s^2[18]，尺寸大小为13.00 km×11.39 km×9.07 km^[7]。

2.2 火卫一重力场计算与分析 2.2.1 引力位结果

根据上文所述，本节以火卫一作为算例，给出了火卫一外部距质心14 km处球面1°间隔格网点上的引力位。图 2(a)~2(c)分别为边长法、面积法、体积法计算所得引力位的结果，图中纵坐标为余纬度。

通过图 2可以发现，3种方法计算的引力位分布情况大致相同，均呈现条带状分布，其引力位在45~56 m²/s²范围间变化。火卫一外部引力位分布存在3处极大值，其中（1ºN, 49°W）附近，即Stickney陨石坑处所处位置存在最大值。

2.2.2 球谐系数结果

本节使用多面体体积积分得到的引力位计算阶次至80阶的球谐系数，并与已发表的重力场模型低阶球谐系数进行对比。表 1为零阶、一阶、二阶的归一化球谐系数。分析表 1数据可得，当坐标系原点为质心，坐标轴和主惯性轴重合时，二阶的球谐系数仅有C₂₀和C₂₂项不为零，这是主惯性轴作为坐标轴时的特征，也是和史弦^[1]及Le Maistre等^13]的结果的不同之处。

由多普勒数据计算所得二阶项系数C₂₀为-0.037 95 ±0.009 63，C₂₂为0.015 87±0.014 73^[19]，与本文所计算的二阶系数相比可发现，C₂₀和C₂₂项均在其误差范围内。因此可认为，在均质密度假设下，基于形状模型所计算重力位有一定参考价值。

为方便进行球谐系数的比较，本文通过功率谱来进行结果对比。功率谱的计算公式如下：

$ {\varepsilon _1} = \sqrt {\frac{{\mathop \Sigma \limits_{m = 0}^1 \left( {C_{lm}^2 + S_{lm}^2} \right)}}{{2l + 1}}} $

(13)

不同阶数功率谱结果对比见图 3。由表 1和图 3可以发现，和已发表的结果对比，低阶项尤其是一阶项和二阶项相差较大，其原因为一阶项与二阶项和坐标系的原点和定向有关。坐标系原点发生变动时，火卫一固坐标系和惯性系转换所导致的差异较大。二阶之后的功率谱吻合较好。虽然坐标系的不同对高阶的球谐系数也会产生影响，单独比较高阶的C_lm和S_lm仍不相同。但对于功率谱来说，相同阶次仅反映重力场频谱信号的强度，这也是高阶项功率谱相差不大的原因。

2.2.3 引力位结果比较

将边长模型、面积模型和体积模型计算所得引力场模型分别作差对比，以比较不同模型之间的差异。由图 4中3种模型所计算的均方根（root mean square，RMS）值可以看出，以体积法计算引力位作为检验值，边长法和面积法的最大误差都集中在Stickney陨石坑附近，地形显著变化处差异较大。其中，由于边长法和面积法仅积分微元不同，因此二者引力位相差较小，但最大误差仍集中在Stickney陨石坑处。而体积法计算引力位时，火卫一外部一点到面积所在平面的垂直距离是导致和面积模型，边长模型所计算引力位相差较大的原因。

2.2.4 加速度结果比较

为量化阶次对加速度产生的影响，本文比较了50阶次和20阶次重力场的加速度差值及80阶次和50阶次重力场的加速度差值（图 5）。火卫一外部14 km球面上的加速度在1~136 mGal间变化，通过图 5发现，50阶和20阶之间，80阶和50阶之间的差异最高分别达到9.55 mGal和0.96 mGal。且在（1°N, 49°W）附近出现最大差值，陨石坑附近加速度差值变化也较为明显。Stickney陨石坑所处位置地形变化较为剧烈，对加速度的影响较大。由图 5还发现，50阶和20阶重力场加速度差值相较80阶和50阶重力场加速度差值更为明显。这一结果表明火卫一高阶次重力场仍会产生显著的重力加速度。

3 结语

本文介绍了基于形状模型计算火卫一重力场的边长法、面积法、体积法3种方法。将火卫一形状模型转换到质心为原点坐标、惯性轴为坐标轴的坐标系中，进而计算引力位及阶数至80阶的高阶球谐系数，其中体积法所计算球谐系数结果更为精确。在地形变化较大处，误差也会相对较大。且通过比较20阶、50阶、80阶重力场所计算的加速度，表明阶次对加速度有明显影响。本文使用引力位计算方法，因为重力场仅在全部包裹天体的布里渊球面外收敛，球面内收敛情况十分复杂，所以选取的计算引力位所在球面需完全包覆天体。布里渊球面内火卫一外的部分则无法计算，这也是之后的研究中需要改进之处。

火卫一重力场信息对火卫一的科学研究和探测有重要意义^[20-21]。本文的计算结果可以提供一个高阶重力场作为任务仿真分析的先验模型，评估MMX任务等的可行性，并针对MMX任务的轨道特性，确认反演模型的有效阶次，为分析火卫一邻域内探测器的轨道特性、轨道设计仿真等提供理论依据；同时也对导航起一定的支撑作用，为探测器着陆定轨提供一个参考重力场模型。