卫星重力测量

卫星重力测量－基础、模型化方法与数据处理算法作者简介：张传定，男，1966年04月出生，1996年09月师从于解放军信息工程大学陆仲连教授，于2000年12月获博士学位。

摘要论文的中心内容是卫星重力测量中如何由星载传感器获得的观测数据恢复地球重力场这一过程的模型化问题。

旨在吸取前人的研究成果，提出更加合理的数据处理模型。

论文最突出的贡献是，改造并完善了大地重力学、空间大地测量、卫星轨道力学等学科模型化的理论与方法以适应卫星重力测量这一新型观测技术。

作者的主要工作和创新点有：1.在综合卫星重力测量有关最新研究成果的基础上，系统地论述了动态加速度测量、卫星重力梯度测量的基本原理；论证了它们的测量精度与姿态角加速度的关系以及卫星重力测量系统最终恢复地球重力场能力的判定准则；深入理解并掌握了现行SST、SGG卫星CHAMP、GRACE、GOCE各项指标及恢复地球重力场各频段的精度指标。

2.简要介绍了卫星重力测量中所涉及到的曲线坐标系下矢量、张量与曲线坐标之间的微分关系、坐标系之间的变换关系以及它们的矩阵表示。

详细研究了在地球重力场确定中常用的关于研究点P和流动点Q相互关联的球极坐标系，给出了球极坐标系下地球引力位V关于P点和关于Q点的微分公式以及它们与球坐标系下局部微分算子的关系。

深入研究了关于P和Q两点局部导数算子的相互作用问题，得到了扰动场元之间核函数和协方差函数的解析与级数展开式，首次给出了较为实用的明晰表达式。

此结果是对物理大地测量学关于这一论题的补充和完善。

这项工作是本文的一个创新点。

3.详细推导了地球、卫星、加速度传感器检验荷载这一特殊限定性三体问题的运动方程；指出星载加速度传感器的输出就是卫星所受非引力加速度和检验荷载相对于卫星中心地球引力的潮汐力之差；进而得到了由星载加速度传感器的比力测量和GPS跟踪测量数据直接恢复地球引力矢量的理论公式。

4.通过对扭秤、旋转梯度仪工作原理的考察和Molodensky关于垂线偏差推求高程异常的论述以及目前业已发现水平梯度分量的某种组合是球面正交函数系的事实，作者明确指出，在地球重力场的研究中，水平方向观测量的组合应作为复数使用。

扰动场元观测量的复数表达是本文立论和各种模型化(建模)工作的思想基础，也是本文最为突出的创新点。

5.在§2.7中，直接由体球谐函数水平梯度的复表示定义并证明了描述地球引力位直到二阶水平梯度所需的球面正交函数系。

它们关于纬度的函数是Legender函数及其导数的拟线性组合，可由目前熟知有关Legender函数及其导数的递推公式给予赋值。

连同球谐函数构成了描述引力矢量、引力梯度张量所需的正交函数系。

因而，利用它们可将引力矢量、引力梯度张量的复分量表达成一致的形式。

6.利用卫星重力测量数据恢复地球重力场，若从边值问题理论上可将其归结为平均轨道面上卫星重力测量超定边值问题。

通常又将利用单个边值条件确定扰动位问题称为单定问题。

在§3中，先以重力异常为例，类比依次给出直到二阶梯度球域单定连续边值问题恢复地球引力位系数的理论公式及其外部解析解和向下延拓截断核函数解；接着导出离散网格平均重力异常对应的简单调和分析公式和最小二乘调和分析公式；然后推广得到广义梯度调和分析公式和超定边值问题的最小方差解、最小二乘解。

并证明了最小方差解等价于单定边值问题调和分析解的频域加权平均；最小二乘解等价于单定离散边值问题最小二乘调和分析解法方程相加所得的解。

广义调和分析方法所需的有关勒让德函数及其导数的积分递推公式在§3.6中给出。

7.首次定义并推导出了水平一阶和二阶梯度平滑因子、。

在概念上澄清了它们与熟知的面球谐函数平滑因子是不同的。

尽管、与相差不大，但在实践上应严格区分它们，这样逻辑上才是严格的。

与观测量的对应关系是使用时应按格网均值数据类型，采用相应的平滑因子。

8.由于水平扰动场元之间的协方差并非各向同性，导致协方差矩阵结构复杂(子矩阵不是Toeplitz循环阵)，不能利用变换矩阵将其降阶，无法付诸实践,迄今尚无最小二乘配置理论应用于水平扰动场元观测量的模型化公式和数据处理方法。

作者研究发现，利用水平梯度的复组合，即复数表示后，扰动重力场元复组合之间的协方差函数尽管还是各向异性，但它们对应的协方差矩阵却具有分块Toeplitz循环阵的结构，因而水平分量复组合的配置问题与重力异常的配置问题相似，可以利用傅立叶变换矩阵进行降阶处理。

这表明，必须将最小二乘配置理论拓展，以适应复数信号的配置问题，本文将其称为最小二乘复配置。

作者将最小二乘配置理论拓展为既能处理复信号又能处理实信号的配置模型，得到了最小二乘复配置解所需的公式。

结合卫星重力测量观测量，详细研究了重力场元复分量之间协方差函数的级数展开式、扰动引力位系数与复分量间的协方差关系。

然后，利用最小二乘复配置理论和重力场复分量之间协方差函数表示，得到了与重力异常配置解接近于一致的各类单定离散边值问题的最小二乘复配置广义调和分析解和超定离散边值问题的最小二乘复配置广义调和分析解。

利用复数表示，解决了最小二乘配置理论难以用于重力场水平观测量这一瓶颈问题。

9.基于摄动力的S、T、W分解，首次给出了卫星受摄运动方程的反解公式。

它只需对轨道根数求一次时间微分，便可求得摄动力的S、T、W分解(不含中心引力)。

与由卫星的地心位置矢量求二次导数得卫星所受力(含中心引力)的直接方法相比，该方法因自然扣除中心引力且只需求时间一次导数，理论上精度要高一些，它具有Hill方程类似的作用，可用S、T、W的时间序列与非心引力矢量的S、T、W分解建立联系，得到时域观测方程。

10.成功地将高斯型求积公式用于常微分方程的数值解中，得到了高斯型隐式龙格库塔方法(GRK)，4级7阶GRK法的精度已与目前较为常用的10级8阶RK法的精度相当。

任意级GRK法的权系数恒为正值，这就是GRK法精度高的本质所在。

利用Legender多项式的零点可以得到任意阶GRK法的权系数，而传统高阶RK法的权系数只能手工推算，而且同牛顿柯斯特积分一样，高阶公式是不稳定的。

11.通过地面跟踪动力法观测方程的建立，得到了SST观测的动力法观测方程。

指出尽管SST测量是地面跟踪的空间拓展，但SST测量因增加了非引力摄动比力加速度观测量，其观测方程与地面跟踪观测方程相比，具有质的差异，这正是SST能以较高精度恢复地球重力场的优势所在。

12.利用最小二乘调和分析公式推导了确定恢复重力场最高阶次的理论准则，其水平敏感度准则系作者给出,利用累积代表误差和截断误差曲线的交点，即可求定各类观测量敏感地球重力场的最高阶次。

13.简要论述了卫星重力梯度测量数据向下延拓的截断核谱组合解及其截断阶次和频域加权准则。

假定梯度测量数据各分量间以及相邻采样点间独立不相关，则频域加权与阶次无关，只与独立观测量的个数有关。

14.配置方法的缺陷是，在数据范围和分辨率一定时，增加观测量类型，意味着协方差矩阵的维数升高。

若是全张量梯度测量数据，则协方差矩阵的维数将是单个分量协方差矩阵的5倍。

为解决矩阵降阶问题，在§7.3中，阐述了最小二乘矢量、张量配置的思想及其算法公式,为利用卫星重力测量数据逼近局部重力场提供了可用于实践的数据处理方法。

关键词：地球重力场，卫星重力测量，Toeplitz矩阵，正交函数系，复数表示，边值问题，调和分析，平滑因子，最小二乘配置，最小二乘复配置，最小二乘矢量配置，最小二乘张量配置，摄动运动方程ABSTRACTThis dissertation is mainly focused on the modeling process of the recovery of the earth's gravity field from observations acquired by sensors onboard satellites, with the aim of summarizing the works done before the author and establishing still more appropriate data processing models. The most important contribution this dissertation presents is that the modeling theories and methods of geodetic gravimetry, space geodesy and satellite orbit dynamics have been reconstructed and perfected so as to adapt to the new surveying techniques of satellite-borne gravimetry. The main work and innovations done by the author are presented as follows:1. On the basis of the generalization of the latest studies of satellite-borne gravimetry, the author discusses the principles of dynamic acceleration measurement and satellite-borne gradiometry systematically, and demonstrates both the relationship between accuracy and angle acceleration and the determinant rules of the ability of the earth's gravity field recovery from satellite-borne gravimetry. Then the author has a deep inquiry into the targets of the current SST and SGG satellites such as CHAMP, GRACE and GOCE as well as the accuracy targets of different frequencies of the earth's gravity field recovery.2. The author then introduces the differential relationship between vectors as well as tensors in a curve coordinate system and curve coordinates, and the transformation relationship between coordinate systems and their matrix expressions. After studying the frequently used spherical polar coordinate system concerning the calculating point P and the mobile point Q in detail, the author presents the differential formulae of the earth's gravitational potential V with respect to the points P and Q and their relationship with partial differential operators in a spherical coordinate system. Then after the author researches the interaction problem of the partial derivative operators with respect to the two points P and Q, analytical expressions as well as series expansions of the kernel function within disturbing gravity field elements are determined, and their explicit expressions which can be applied practically are given for the first time. This achievement, being an innovation of the dissertation, is a supplement and perfection of the topic in physical geodesy.3. The author deduces the motion equation of the specially restricted three-body problem that consists of the earth, satellite and acceleration sensor test payload in detail, and points out that the output of the acceleration sensor onboard the satellite is the very difference of the non-gravitation acceleration acting on the satellite and the tidal forces of the test payload relative to the gravitation at the center of the satellite, and thus obtains the theoretical formula of the earth's gravity vector recovery from measurement of the acceleration sensor onboard satellite and the GPS tracking observations.4. The author points out that the combination of horizontal observable should be used in the form of the complex number in the study of gravity field through an inspection of theprinciples of torsion balance and rotational gradiometer and the discussion of the acquisition of height anomaly from vertical deflections according to Molodensky's, as well as the fact that some combination of horizontal gradient components are orthogonal functions on a spherical surface. The complex expression of the observations of disturbing gravity elements is the basis of the argument of this dissertation and all kinds of modeling processes and it is also the most prominent innovation.5. In section 2.7, according to the complex expression of the horizontal gradient of spherical harmonic, the spherical orthogonal function family needed for the description of the earth's gravity potential up to 2 degree horizontal gradients is defined and proved directly. The functions, latitude being parameter, are the quasi linear combination of the Legender functions and their derivatives, which can be evaluated by the recursive formulae of the Legender function and its derivatives. Therefore, with these functions, the complex components of the gravity vector and the gravity gradient tensor can be expressed in a consistent form.6. When it comes to the gravity field recovery from satellite-borne gravimetry observations through the theory of the boundary value problem, this problem can be reduced to the over-determined boundary value problem of satellite-borne gravimetry on an average orbital plane. Usually the determination of the disturbing gravity potential with a single boundary value condition is called a single-determined problem. In section 3, taking the gravity anomaly for example firstly, the author presents theoretical formulae of the recovery of the earth's gravity potential coefficients, and its external analytical solution as well as the truncated kernel solution of downward continuation from up to two degree gradient by analogy, as to the single-determined continuous boundary value problem in sphere region. The author also presents simple harmonic analysis formulae and least squares harmonic analysis formulae corresponding to discrete grid mean gravity anomaly. Then these formulae are generalized to obtain what is called the generalized gradient harmonic analysis method in the dissertation. This step is followed by the minimum variance solution and least squares solution to the over-determined boundary value problem. It is also proved that the minimum variance solution is equivalent to the weighted average in frequency domain of the harmonic analysis solution to the single determined boundary value problem, and that the least squares solution is equivalent to the summation of the normal equations of the least squares harmonic analysis solution to the single determined discrete boundary value problem. The recursive formulae of the Legender and its derivatives requisite for the generalized harmonic analysis method are given in section 3.6.7. It is for the first time that the horizontal one- and two- degree gradient smoothing factors , are defined and deduced. And the conceptional difference between , and the well-known surface spherical harmonic is clarified. Although there is small difference between , and , yet they should be strictly distinguished in practice so as to keep strict in logic. The corresponding relation of and observable is shown as follows:In application, the selection of the smoothing factors should be based on the type of the grid mean values.8. On the ground that the covariance between the disturbing gravity field elements is notisotropic, the covariance matrix has a complicated structure (its submatrix being not repetitive Toeplitz matrix) and thus its degree cannot be degraded with transformation matrix method, so it cannot be put into practice. So far, the application of least squares collocation to the modeling and data processing of the horizontal disturbing gravity field observations has not been found in literature. After studying this phenomenon, the author finds that, the covariance matrix of the complex expressions of horizontal gradients of the disturbing gravity field elements has a structure of blockwise repetitive Toeplitz matrix, though the covariance within the disturbing gravity field elements is still anisotropic. Therefore, the collocation problem of the complex combination of horizontal components is similar to that of gravity anomaly, which makes the degradation be able to be completed by FFT transformation matrix. This indicates that the least squares collocation theory must be developed so as to adapt to the collocation of complex signals. And this is what is called the complex least squares collocation in the dissertation. The author develops the least squares collocation into the collocation model which can deal with both complex and real signals, thus obtaining formulae for complex least squares collocation. In reference to the satellite-borne gravimetry observations, the author studies the series expansion of covariance between the gravity field elements' complex components and the covariance relation between the disturbing gravity potential coefficients, and the complex components in detail. Then, through applying the complex least squares collocation theory and the covariance expressions between the complex components of the gravity field, the author acquires the generalized harmonic analysis solution of the complex least squares collocation of all kinds of single determined discrete boundary value problems and such solution to over-determined discrete boundary value problems. The solution is nearly consistent with the result of gravity anomaly collocation. By means of the complex expression, the author overcomes the bottleneck problem that it is hard to apply the least squares collocation theory to the horizontal observations of the gravity field.9. Based on the S, T and W decomposition of perturbing forces, the formulae of the reverse solution to the satellite's perturbed motion equation is given for the first time. As long as the first differential of the orbital elements with respect to time is made, the S, T, and W decomposition of perturbing forces (excluding the earth's central gravitation) can be obtained. Compared with the direct method of acquiring the forces exerted on satellite(containing the earth's central gravitation) from the second derivative of the satellite's geocentric position vector, this method is still more accurate theoretically in that it excludes the earth's central gravitation and just needs the first derivative with respect to time. And this method works like the Hill's equation, so the observation equation in time domain can be obtained by relating the time sequence of S,T and W with the S,T and W decomposition of non-central gravitation.10. The author successfully applies the Gauss-type integral formula to the numerical solution to the constant differential equation and creates the Gauss-type implicit Runge-Kutta Method(GRK) of which the accuracy of order 4 and degree 7 is already equivalent to that of the RK method with order 10 and degree 8 that is in common use. The GRK method of any order has a positive weight coefficient, which leads to its higher accuracy. The weight coefficients of GRK of any order can be acquired through the zero points of the Legender polynomial, whereas the weight coefficients of the traditional RK of higher orders can only be solved manually and，like the Newton-Coast Integral, the formula of higher order is unstable.11. Through the establishment of observation equation of ground tracking dynamic method,the dynamic observation equation of SST is obtained. The author points out that there is one more observable, that is to say, the acceleration of non-gravitational perturbing forces in the SST measurement, though the SST measurement is the extension of ground tracking. Therefore, the observation equation of the SST measurement is essentially different from that of the ground tracking, which is the advantage of the SST's recovery of the earth's gravity field with higher accuracy.12. By means of the least squares harmonic analysis formulae, the theoretical rule in the determination of the highest order and degree of the earth's gravity field recovery is deduced, in which the horizontal sensitivity rule is the author's. The highest order and degree of the gravity field sensed by all kinds of observations can be determined via the intersections of the accumulated representative error curve and the truncation error curve.13. The author discusses the solution of the truncated kernel spectral combination of the downward continuation of the satellite gravity gradiometry observations and the determinant rule of its truncation order and degree as well as weighting in frequency domain. Assuming that the components of gradient observation and the neighboring sampling points are independent of each other, the weighting in frequency domain has nothing to do with order and degree, only depending on the amount of independent observable.14. The limitation of the collocation method is that, the increment of observable types may cause the increase of the dimension of covariance matrix when the range and resolution of data are fixed. If it happens to be the full-tensor gradient observation data, then the dimension of covariance matrix will be five times that of single component covariance matrix. In order to solve this problem, the author discusses the practical computing method in applying the least squares collocation of vector and tensor, thus providing a data processing method that can be put into practice of the approximation of the local gravity field from satellite-borne gravimetric observations.Keywords: earth's gravity field, satellite-borne gravimetry, Toeplitz matrix, orthogonal function family, complex expression, boundary value problem, spherical harmonic analysis, smoothing factor, least squares collocation, least squares complex collocation, least squares vector collocation, least squares tensor collocation, perturbed motion equation。