第1章卫星重力场恢复的参数估计
于尔根 库什和安妮 斯普林格
摘要：简单来说，这里的参数估计是从不确定、错误、冲突，甚至可能同时不完整的数据中，提取地球物理系统或测量系统中具有明确定义的相关值最佳猜测的过程。我们可以通过直接观测的方式，从数据中得到有关所需参数的信息，但是在很多情况下数据与所需参数之间的关系是间接的。实际上，我们永远无法找到“绝对真实”的参数，但可以寻求在某方面最佳的估计参数，比如对于给定类型的数据该估计参数传播误差最小。
所有的自然科学以及与“经验性”相关的大部分社会科学，都会面临参数估计这一问题。在卫星重力场测量领域，地球物理场是所要求解的相关量，它包含无限阶次的待定位系数，并且没有*佳的截断阶次，这使得标准方法和流传数百年的方法[如*小二乘(least square，LS)法]的应用更加复杂，并造成了来自不同学术背景的数据分析师之间工作的混乱。基于泛函分析的构造逼近理论提供了一种将问题转化为求解有限数量参数的现代方法，并且不会遗漏掉重要细节。但是，逼近理论的误差界限取决于一些非常抽象的概念，比如描述空间光滑性的核函数等。
由于上述原因，在卫星重力场测量领域中*小二乘法及其他算法都会受到质疑。然而，最小二乘法也存在一定优势，它不仅易于应用，而且可对实际测量系统引入的不确定性进行误差估计并提供解决方法。
最小二乘法可以解决当前重力卫星数据分析中面临的一些挑战：①高效处理大量数据，获取所需的多项参数；②全面量化数据误差和球谐分析中的背景模型误差带来的影响，并定量评价由此衍生的科学成果；③在存在不一致性的情况下，将不同仪器、卫星和处理链的数据集不间断地组合起来。这种不一致性是已知存在的，但很难被精准确定，并且从基本原则的角度来看其移除成本过高。
最小二乘法的发展源于对拟合方法的需求，这种拟合方法在古代已经实现，在最早的天文/大地观测中就涉及地球半径以及月球、地球轨道半径的拟合[109]。众所周知，在大地测量领域中，为确定椭球形地球的扁率，人们进行了子午线测量。在这些测量过程中，我们无法直接从一个单*的测量中导出所需参数，而是对测量数据进行了平均处理。拉普拉斯(Laplace)提出了几种拟合方法，可以用于减小误差，例如早期用于减小拟合椭球面和数据点之间的最大误差，后来用于减小误差的平均绝对值，使全部误差的绝对值总和为零。这些拟合工作受到了意大利天文学家和数学家博斯科维奇(Boscovitch)的影响，这位并不为人所熟知的学者致力于研究拟合问题，*次提出了“平差理论”的这一基本原理。勒让德(Legendre)提出了可用于弧度测量分析的*小二乘标准方程。尽管当时并未采用矩阵表示法，高斯(Gauss)提出了我们今天所知的加权*小二乘法，他将最小二乘法用于行星及小行星轨道确定，也用于大规模三角测量校准中。
这些讲义是根据“基于星星跟踪测量的全球重力场恢复”国际秋季学校讲座编写的，该秋季学校于2015年10月4日～9日在德国巴特洪内夫举办，主办方为DFG下属的协同研究中心项目团队“geo-Q”。国际秋季学校的主旨是向来自不同学术背景的学生们介绍重力卫星数据分析中有关参数估计的一些常见而有用的概念，例如求解一组球谐系数。讲座的重点在于介绍概念，因此省略了技术性证明。
1.1符号说明
矢量由a、x等小写字母表示，矩阵则由A、X等大写字母表示。这样可以使符合表示更加统一，但是更重要的目的是区分随机变量与它们的实现。
表示期望算子，表示方差算子，表示参数的估计值，注意和不同。表示球谐函数，当时它等于；当时它等于。这些函数默认是完全规格化的，因而省略了上横线符号。
1.2相关文献
平差问题在不同的科学领域得到了研究，因而存在大量的关于平差问题的参考文献。这里，我们对大地测量领域文献进行有选择性的概述。
文献[109]详细描述了*小二乘法平差理论的发展历史，同时列举了几个**的大地测量问题，比如与地球形状相关的问题。为深入了解平差理论及统计学领域问题，我们**文献[85]。文献[13]从一位应用数学家的视角来探讨这个问题，主要集中于离散型最小二乘问题，其中包括有效求解器的设计问题。对于卫星重力场测量领域的文献，我们重点**文献[127]，它讨论了以去相关和迭代解为核心的大型平差问题数值策略。文献[22]进一步拓展了上述工作，解决了大规模并行编程环境下利用*小二乘法求解重力场参数的若干问题。文献[19]综述了可应用于无约束和有约束条件下的全球重力场分析/求解的质量度量问题。
卫星重力场恢复是一种不适定问题，解决该问题的关键在于正则化。文献[33]针对GOCE重力场测量中的正则化问题开展了深入研究。方差分量估计(variance-component estimation，VCE)是正则化参数选择的有效方法，1979年发表的开创性论文[40]提出了VCE方法的迭代过程。如果想了解*大似然框架下方差分量估计的理论推导与证明，文献[84]是很好的参考材料。文献[87]对方差分量估计进行了更深入的研究，并与Lerch*优加权法进行了比较，此外还介绍了一种可用于方差分量估计的蒙特卡罗方法。文献[145]对几种适用于地球重力场恢复的方差分量估计方法进行了全面讨论，其中包括对辅助参数的假设检验。
1.3高斯-马尔可夫模型
在对测量数据yi.1，2，.n的统计分析中，最重要的模型之一就是高斯-马尔可夫模型(Gauss-Mark，ov model，，GMM)。现代大地测量领域的很多问题无法通过高斯-马尔可夫模型来表示，需要应用参数估计理论，包括整数-浮点数混合估计[例如，在全球导航卫星系统(global navigation satellite system，GNSS)模糊度解算问题中]、乘性噪声估计(例如，在雷达测高和InSAR中)、分类问题(例如，在测高和遥感中)等。
1.3.1高斯-马尔可夫模型的基本假设
在高斯-马尔可夫模型使用中，我们假设：①数据的数学期望可以与一组未知但固定的参数，线性相关，不严格地讲，这里的数据要求是无错误的；②使用者熟练掌握数据误差的方差和协方差知识。
高斯-马尔可夫模型提供的分析框架非常通用，例如它可以应用于时间序列分析。与待估参数相关的信息可以同时在不同的地点采集，也可以利用同一设备在不同的时刻采集，或者是在这些情况简单组合的条件下采集。
1.3.1.1高斯-马尔可夫模型公式
采用矢量形式来表示数据和参数，高斯-马尔可夫模型*基本的公式如下：
(1.1)
其中，y是大小为的数据矢量；x是大小为的参数矢量；A是满秩的.观测方程矩阵或设计矩阵.为方差因子。假设所有的观测量具有相同的误差方差，且互不相关。
期望和方差可表示为
不严格地来讲，期望E.y.可以对无限数量的数据矢量取平均值得到，其中
每个观测量出现的频率用概率密度函数p.y.表示。方差可以通过计算无限数量的数据矢量与其期望的偏差，进而计算每个偏差与其转置乘积的期望得到。在对为误差方差；在非对角线上，角线上，为误差协方差。这里假设协方差
为0。对于未知的观测误差，存在
(1.2)
也就是说，对于第i个观测量.可以通过对取平均得到。
可以通过对ee取平均得到。
1.3.1.2高斯-马尔可夫模型的性质
在式(1.1)中，y为观测值矢量。矩阵A是已知的，它反映了观测矢量与待估参数x之间的数学或物理联系，是将可观测变量与隐藏参数相关联的模型。这里，仅考虑nm.这种*常见的情况，即观测数据比待估参数的个数多。显然，矩阵A具有不同的数学形式，这与待估参数x的定义有关。
例1.1
利用一组观测量拟合多项式，**种模型为
第二种模型为
很明显，
在**种模型中，设计矩阵的第
行对应第i个观测值。
待估参数x是未知的，但也是确定不变的、非随机的。这一点与贝叶斯分析方法是不同的，因为在后者中待估参数被认为是随机的，需要指定待估参数的概率分布函数。在贝叶斯分析中，上述模型可表示为
其中，分别表示在待估参数x和.2给定的条件下。
当进行重力场参数估计时，需要建立均匀采样的测量序列(如重力卫星测量中的测距数据)与地球重力场中的球谐系数Cnm、Snm之间的联系。另一方面，我们可以从时间序列的角度看待测量数据流，将p个连续样本在特定的时间窗口内拟合到与重力场无关的简单模型中，如简单多项式、切比雪夫多项式等。在不引入新的误差的条件下，上述公式改写为，该方程不唯一且无解，至少对于一个给定的测量数据矢量是无解的。暂时假设A矩阵满秩，后面可以把这一条件放松。
按照定义，观测误差e是未知的。另外，在公式中无须假设误差呈正态分布。在高斯-马尔可夫模型中，数据误差被认为是随机的，所做的唯一假设是它们的期望值为零，并且其协方差矩阵中某些因子已知。特别需要说明的是，我们无须假设误差服从正态分布。这意味着有
假设方差因子.2是未知的，需要进行估计。事实上，在高斯-马尔可夫模型的标准解中不需要用到.2值，因此可以根据从观测数据到“预测数据”的拟合结果来估计。为保证观测值的先验精度，通常使用一个整体比例因子来确定方差因子。如果数据误差方差的先验假设成立，那么可以估计出一个接近1的方差因子。同时，这就意味着如果事先已知.2值，对它进行估计也会是一个很好的检验。上述分析中假设误差协方差矩阵为对角矩阵，后面可以将这一假设放松。
1.3.1.3高斯-马尔可夫模型的一般化
在高斯-马尔可夫模型中，可以将所有的单个观测结果yi与未知参数相关联。但是实际上，在一个方程中可能会有多个观测结果。这样，就会出现更多通用模型，即更多地将原始数据与未知数据相关联的线性函数，比如，这里不考虑这些情况。其中的一个例子是雷达测量中用于表示散斑噪声的乘性噪声模型，比如，其中。
在实际情况中，矩阵A中经常包含观测量。例如，在卫星重力场测量中，需要利用卫星轨道数据来构造矩阵A中的元素。矩阵A中的观测误差可能非常关键，但事实上在标准高斯-马尔可夫模型中并没有考虑这一误差。有一些矩阵考虑了随机误差更通用的统计模型，比如高斯-赫尔默特模型(Gauss-Helmert model)、变量误差和整体*小二乘模型等，但是这些模型通常难以应用。在重力场测量领域，更常见的做法是将这些误差参数化，即将它们与确定性参数相关联，这样也会增加参数矢量的长度。
1.3.2高斯-马尔可夫模型中的参数估计
在上述模型中，数据及其误差的“真实值”是未知的，因此只能推断出参数矢量x的统计估计值x.。统计估计值的精度取决于未知的数据误差，因此估计值相对期望的偏差也是未知的。因此，我们希望能够定理确定估计x.时的统计不确定性。通常，这意味着需要推导协方差矩阵.
的估计值。
为实现这一目的，需要引入统计理论，基于某种原理将一些可以合理解释的量最小化，从而实现估计过程。换句话说，我们需要将合理性引入数学框架中。对于上述模型而言，最佳线性无偏估计(best linear unbiased estimation，BLUE)、最小二乘法和最大似然(maximum likelihood，ML)法假设误差服从高斯概率密度函数(probability density function，PDF)分布，可以得到相同的参数估计结果：
(1.3)
当时，定义该估计过程为无偏估计。这意味着，如果能够采集到的同类观测数据越多，那么估计值就越接近真实值。在这种情况下有
这说明估计是无偏的。在高斯-马尔可夫模型中，无偏估计通常称为线性回归。通过对观测数据的加权组合得到估计值。如果观测数据误差服从正态分布，那么估计误差也服从正态分布。显然，计算x.时不需要已知方差因子的值。