根据张量不变理论，对称二阶引力梯度张量$\mathit{\Gamma }=\left[V_{ij}\right](i,j=\mathrm{1,\; 2},\mathrm{3})$具有3个独立的标量不变量，即一、二、三阶主子式之和，它们组成了一个完备的不变系统。该不变系统可表达为^[9]：

式中，V_ij（i，j=1，2，3）表示引力梯度张量的6个分量；I₁、I₂、I₃为引力梯度张量的3个不变量。

在局部指北坐标系中，引力位可表达为：

式中，$G\mathrm{、}M$为地心引力常数；$R$为地球平均半径；$r\mathrm{、}\theta \mathrm{、}\lambda$分别为地心向径、地心余纬和地心经度；${\overline{P}}_{nm}(\mathrm{c}\mathrm{o}\mathrm{s}\theta )$为完全规格化缔合勒让德函数，简记为${\overline{P}}_{nm}$；${\overline{C}}_{nm}$、${\overline{S}}_{nm}$为完全规格化的引力位球谐系数；n、m为球谐展开的阶和次。将引力位截断至一定的阶次N，代入式（1），可得到3个不变量的表达式^[9]。

该不变量是引力梯度张量系数矩阵3个主对角元素的和，且引力位系数与不变量观测值间的函数模型为线性模型，即：

式中，$W_{n_{\mathrm{1}}{m}_{\mathrm{1}}}(\theta )$为余纬$\theta$和阶次$n_{\mathrm{1}}$、$m_{\mathrm{1}}$的函数：

式（4）即勒让德差分方程。

I₁不变量为引力梯度张量的迹。由于$\mathrm{t}\mathrm{r}\mathit{\Gamma }=\mathrm{0}$，因此无法利用其进行地球重力场的恢复。但是它可以作为最小二乘方法的先验约束条件，从而提高求解的精度。

该不变量为引力梯度张量6个分量的组合形式，是引力梯度分量平方的线性组合形式。但相对于球谐系数，I₂为一次非线性模型。

$\begin{array}{c}{I}_{\mathrm{2}}=-\frac{\mathrm{1}}{\mathrm{2}}(V_{\mathrm{11}}^{\mathrm{2}}+V_{\mathrm{22}}^{\mathrm{2}}+V_{\mathrm{33}}^{\mathrm{2}})-V_{\mathrm{12}}^{\mathrm{2}}-V_{\mathrm{13}}^{\mathrm{2}}-V_{\mathrm{23}}^{\mathrm{2}}=\\ {\left(\frac{GM}{R^{\mathrm{3}}}\right)}^{\mathrm{2}}{\displaystyle \sum _{n_{\mathrm{1}}=\mathrm{0}}^N}{\displaystyle \sum _{m_{\mathrm{1}}=\mathrm{0}}^{n_{\mathrm{1}}}}{\displaystyle \sum _{n_{\mathrm{2}}=\mathrm{0}}^N}{\displaystyle \sum _{m_{\mathrm{2}}=\mathrm{0}}^{n_{\mathrm{2}}}}{\left(\frac{R}{r}\right)}^{n_{\mathrm{1}}+n_{\mathrm{2}}+\mathrm{6}}\times \\ [K_{n_{\mathrm{1}}{m}_{\mathrm{1}}{n}_{\mathrm{2}}{m}_{\mathrm{2}}}^{\mathrm{1}}(\theta )(-{\overline{C}}_{n_{\mathrm{1}}{m}_{\mathrm{1}}}\mathrm{s}\mathrm{i}\mathrm{n}{m}_{\mathrm{1}}\lambda +{\overline{S}}_{n_{\mathrm{1}}{m}_{\mathrm{1}}}\mathrm{c}\mathrm{o}\mathrm{s}{m}_{\mathrm{1}}\lambda )\times \\ (-{\overline{C}}_{n_{\mathrm{2}}{m}_{\mathrm{2}}}\mathrm{s}\mathrm{i}\mathrm{n}{m}_{\mathrm{2}}\lambda +{\overline{S}}_{n_{\mathrm{2}}{m}_{\mathrm{2}}}\mathrm{c}\mathrm{o}\mathrm{s}{m}_{\mathrm{2}}\lambda )+K_{n_{\mathrm{1}}{m}_{\mathrm{1}}{n}_{\mathrm{2}}{m}_{\mathrm{2}}}^{\mathrm{2}}(\theta )\times \\ ({\overline{C}}_{n_{\mathrm{1}}{m}_{\mathrm{1}}}\mathrm{c}\mathrm{o}\mathrm{s}{m}_{\mathrm{1}}\lambda +{\overline{S}}_{n_{\mathrm{1}}{m}_{\mathrm{1}}}\mathrm{s}\mathrm{i}\mathrm{n}{m}_{\mathrm{1}}\lambda )({\overline{C}}_{n_{\mathrm{2}}{m}_{\mathrm{2}}}\mathrm{c}\mathrm{o}\mathrm{s}{m}_{\mathrm{2}}\lambda +\\ {\overline{S}}_{n_{\mathrm{2}}{m}_{\mathrm{2}}}\mathrm{s}\mathrm{i}\mathrm{n}{m}_{\mathrm{2}}\lambda )](\mathrm{5})\end{array}$

球谐系数的乘积形式，例如${\overline{C}}_{n_{\mathrm{1}}{m}_{\mathrm{1}}}{\overline{C}}_{n_{\mathrm{2}}{m}_{\mathrm{2}}}$、${\overline{C}}_{n_{\mathrm{1}}{m}_{\mathrm{1}}}{\overline{S}}_{n_{\mathrm{2}}{m}_{\mathrm{2}}}$等，导致了不变量$I_{\mathrm{2}}$的非线性特征。其中，n₂、m₂分别为球谐展开的阶和次，且：

式中，$\frac{\mathrm{d}{\overline{P}}_{n_{\mathrm{1}}{m}_{\mathrm{1}}}}{\mathrm{d}\theta }$表示完全规格化缔合勒让德函数${\overline{P}}_{n_{\mathrm{1}}{m}_{\mathrm{1}}}$相对于地心余纬$\theta$的导数。

不变量I₃亦为一非线性函数模型：

相对于不变量I₂，I₃的非线性特征由球谐系数的乘积形式$C_{n_{\mathrm{1}}{m}_{\mathrm{1}}}{C}_{n_{\mathrm{2}}{m}_{\mathrm{2}}}{C}_{n_{\mathrm{3}}{m}_{\mathrm{3}}}$等引起，为二次非线性模型。其中，$n_{\mathrm{3}}\mathrm{、}{m}_{\mathrm{3}}$为球谐展开的阶和次，且：

需要指出的是，梯度仪相对于惯性空间的旋转特性所带来的旋转效应会对张量不变量的计算带来比较大的影响。因此，在利用卫星重力梯度数据恢复地球重力场时，需要首先从观测张量中消去欧拉和离心张量^{[5，9，17-19]}。为了消去这种旋转效应，在张量不变方法中，角速度和角加速度信息同样需要获取，但不需要精确的定向参数信息，这与经典的重力梯度反演方法有别。下文的研究中，假定已经对引力梯度张量和旋转张量进行了分离处理。

为了在最小二乘准则下利用张量不变量进行地球重力场的恢复研究，寻找有效的线性化方法，实现函数模型式（5）和式（8）的线性化处理至关重要。用于评价线性化效果的标准一般有3个：（1）线性化误差的量级；（2）收敛特性的优劣；（3）数值运算的效率。对于大型的最小二乘问题，第3个标准特别重要。Brute-force线性化方法，即利用泰勒级数展开实现模型的线性化，并不适用于大型的最小二乘问题，其带来的嵌套求和过程将会增加大量的计算量，即便对于中低阶的重力位模型的求解亦会存在这种问题。

针对张量不变量观测值恢复地球重力场中的实际问题，本文采用文献[9]中的线性化方法对张量不变量进行了线性化处理^[9，20]，最大限度地减少了海量数据带来的计算负担，即利用一先验的参考重力位模型计算张量不变量的扰动量，同时达到线性化处理的目的。从第二次迭代开始，可以利用上次迭代的求解结果来建立线性函数模型。张量不变量扰动量的表达式为^[9]：

式中，$I_{\mathrm{2}}^{\mathrm{r}\mathrm{e}\mathrm{f}}\mathrm{、}{I}_{\mathrm{3}}^{\mathrm{r}\mathrm{e}\mathrm{f}}$分别为利用参考模型进行球谐综合得到的张量不变量；U_ij（i，j=1，2，3）为参考模型的引力梯度张量分量，且$V_{ij}=U_{ij}+\mathrm{\delta }{V}_{ij}$，$\mathrm{\delta }{V}_{ij}$为对应于参考量U_ij的改正量。

该线性化处理可以在任意的参考系下进行，即线性化结果与引力梯度张量的方向无关。当参考梯度值的频谱分辨率与参数估计的频谱分辨率相等时，即阶数N^ref=N时，线性化效果最优。该线性化方法增加的计算量主要集中在每一次迭代前参考模型值U_ij的球谐综合计算。

本文利用EIGEN-6C4的前300阶计算得到了30 d、5 s采样的卫星引力梯度观测值，卫星轨道仿真计算时的初始轨道根数为：轨道长半径6 628.000 km，偏心率0.001，轨道倾角96.6°，升交点赤经0.0°，近地点角距0.0°，平近点角0.0°。图 1给出了全张量梯度观测值下，利用先验重力场模型进行模型线性化后的求解结果。为有效测试解的收敛特性，初始线性化位模型选为EGM96模型的前200阶。张量不变量选为I₂，解算阶数为200阶。从图 1中可以看出，在第2次迭代后，精度不再有明显提高，即第2次迭代已满足收敛条件。此外，第1次迭代后的模型阶误差明显小于单独利用V₃₃分量的求解精度，说明了线性化方法的误差量级非常小，可有效实现张量不变模型的线性化处理。

图  1  张量不变方法的线性化分析

Figure 1.  Tailored Linearization Analysis of Tensor Invariant Method

图 2显示了分别利用EGM96模型和GGM05C模型作为初始线性化模型的第一次迭代后的阶误差。从图 2中可以看出，以GGM05C作为线性化模型第一次迭代后的精度较EGM96模型要好，其差异主要集中在长波部分，且GGM05C模型第一次迭代阶误差非常接近于线性化后的最终收敛值。这是由于GGM05C与EIGEN-6C4模型中均使用了GRACE观测数据参与解算，相对于EGM96模型，GGM05C初始线性化模型更为精确。初始线性化模型的优劣与初始迭代精度是密切相关的，但对最终收敛结果的影响非常小，不会导致收敛特性的改变。在实测数据处理时，为了提高收敛速度，保证收敛特性不变，应选取精度更高的重力场模型作为初始线性化模型。

图  2  参考模型选择对线性化的影响

Figure 2.  Influences of the Reference Model on the Tailored Linearization

从图 1和图 2中还可以看出，张量不变方法的重力场模型解算精度要优于V₃₃分量的解算精度。一方面，本文梯度观测值仿真计算时采用了300阶模型，而解算只解算到了200阶，这主要是考虑到卫星在轨的实际情况，这种频谱混淆在实际重力场模型解算中也是存在的。另一方面，由于张量不变法采用了先验模型进行线性化，观测数据的数学统计特性同时得到了改善和优化，加之线性化误差较为微小，因此其精度结果优于V₃₃分量。

函数模型（5）和（8）的计算需要全张量梯度仪的观测值。不同梯度张量分量的误差特性会统一反映在张量不变量中。因此，所有的梯度张量分量应具有同等量级的测量精度。对于GOCE重力任务，张量分量V₁₂和V₂₃的精度较之其他分量要低，需要寻找一定的方法对这种情况进行处理。本文中，在进行不变量的计算时，首先假定V₁₂和V₂₃观测值不可知，利用先验模型，结合球谐综合方法得到或直接将其赋零，然后进行重力场模型的解算。

图 3给出了将V₁₂和V₂₃置为零后，利用张量不变方法解算得到的重力场模型阶误差。从图 3中可以看出，第一次迭代的精度较之全张量梯度观测值的求解结果要低一个数量级，但第二次迭代即满足收敛条件，精度仅仅稍逊于全张量的求解结果。图 4给出了利用EGM96位模型的球谐综合值代替V₁₂和V₂₃分量后的张量不变方法的解算结果。从图 4中可以看出，第一次迭代结果较之赋零策略的结果精度要高，且第二次迭代收敛值与全张量求解结果吻合很好。这说明，对于非全张量梯度观测值，在利用张量不变方法求解时，借助先验位模型的球谐综合值代替低精度分量的数据处理策略是可行和有效的。

图  3  V₁₂和V₂₃分量置零后的张量不变解

Figure 3.  Tensor Invariant Solution After V₁₂ and V₂₃ Are Set to Zero

图  4  V₁₂和V₂₃分量球谐综合计算后的张量不变解

Figure 4.  Tensor Invariant Solution After Spherical Harmonic Synthesis of V₁₂ and V₂₃

GOCE卫星搭载的梯度仪具有天底指向特性，使得引力梯度张量非对角元素的变化范围在几个Eötvös左右。与之相比，主对角线分量的量级要大3个量级左右。图 5给出了将低精度观测值V₁₂和V₂₃分别置为零以及利用EGM96位模型球谐综合值进行代替时的张量不变量的相对误差。从图 5中可以看出，将V₁₂和V₂₃置为零或者利用先验位模型的球谐综合值进行代替对张量不变量的量级影响很微小，这种良好的数学特性也是选择不变系统进行重力场研究的初衷之一。因此，在利用张量不变方法求解卫星引力梯度重力场时，建议借助先验位模型的球谐综合值代替低精度分量进行处理。

图  5  非全张量不变量的相对误差

Figure 5.  Relative Errors of Non-full Tensor Invariants

对梯度张量模拟值的V₁₁、V₂₂、V₃₃和V₁₃分量分别加入标准差为3.2 mE的白噪声，利用张量不变方法进行重力位模型的解算。其中，V₁₂和V₂₃分量采用EGM96模型的球谐综合值代替，正则化方法选用Kaula正则化，结果见图 6。从图 6中可以看出，白噪声情况下，张量不变方法的模型阶误差要优于单分量的解算结果，这与无噪声情况下的结果一致。但迭代次数的增加并不能改进模型阶误差的精度。这说明模型的线性化误差以及采用先验模型的球谐综合值代替低精度分量所引入的误差较之白噪声水平要低，利用先验重力场模型进行线性化，以及球谐综合值代替低精度分量的数据处理策略仍然是有效的。

图  6  白噪声情况下张量不变解的模型阶误差

Figure 6.  Model Degree Errors of Tensor Invariant Solution in White Noise

本文对张量不变理论在卫星重力梯度数据恢复地球重力场中的应用进行了数值分析。针对观测模型的非线性特征以及GOCE任务的非全张量观测特征，对基于先验重力场模型的线性化方法以及球谐综合方法的数据处理策略进行了探讨，结果表明：

1）张量不变方法不依赖于梯度仪的空间定向信息的获取，较之传统的地球重力场解算方法具有特别的优越性。该方法实现了不同观测值的联合求解，降低了计算负担，特别是对于海量的卫星重力梯度观测值尤为可贵，且张量不变解的精度明显高于单梯度分量的求解结果。

2）基于先验重力场模型的线性化方法能够在保证计算效率的前提下，有效实现观测模型的线性化，解的收敛非常快，精度较高，线性化误差较小，且初始线性化模型的选择对解的收敛特性影响较小。

3）利用先验位模型的球谐综合值代替低精度梯度分量实现非全张量梯度观测值的处理是有效的，不会改变解的收敛速度和最终精度。这主要源于梯度仪的天底指向特性。

4）白噪声情况下，基于先验重力场模型的线性化方法以及采用先验模型的球谐综合值代替低精度分量的数据处理策略仍然是有效的，其引入的误差较之白噪声水平要低，从而迭代次数的增加并不能带来模型阶误差精度的显著提高。