Design and Comparison of Filters of Gravity Field Inversion from GOCE Satellite Gravity Gradient Data
-
摘要: 地球重力场和海洋环流探测（gravity field and steady-state ocean circulation explorer，GOCE）卫星重力梯度数据有色噪声和低频系统误差的滤波处理是反演高精度地球重力场的一个关键问题。针对GOCE卫星重力梯度数据的滤波处理，基于移动平均（moving average，MA）方法和CPR（circle per revolution）经验参数方法设计了两类低频系统误差滤波器，并分别将这两类滤波器与基于自回归移动平均（auto-regressive and moving average，ARMA）模型设计的有色噪声滤波器组合起来形成级联滤波器。为了分析滤波器处理的实际效果，基于空域最小二乘法采用70 d的GOCE观测数据，并联合重力恢复与气候实验（gravity recovery and climate experiment，GRACE）数据分别反演了224阶次的重力场模型GOGR-MA（MA+ARMA级联滤波）和GOGR-CPR（CPR+ARMA级联滤波）。将反演模型与采用同期数据求解的第一代GOCE系列模型及GOCE和GRACE联合模型GOCO01S进行比较，并通过高精度的全球定位系统水准数据和稳态海面地形模型进行检核，结果表明：构建的MA和CPR经验参数滤波器均能削弱低频系统误差的影响，并且后者效果更为明显，而ARMA去相关滤波器能够有效地对重力梯度数据有色噪声进行白化处理；反演的联合模型GOGR-MA和GOGR-CPR的精度接近，并且它们都优于欧空局采用同期数据研制的GOCO01S模型。反演结果验证了设计的级联滤波器的正确性和有效性，可为GOCE卫星重力梯度数据的滤波处理及高精度重力场反演提供参考。
-
关键词:
- 卫星重力梯度 /
- 有色噪声 /
- 低频系统误差 /
- 滤波器设计 /
- 地球重力场和海洋环流探测重力场反演
Abstract:Objectives The filtering of colored noise and low-frequency systematic errors of gravity gradient data from gravity field and steady-state ocean circulation explorer (GOCE) satellite is a key issue in the inversion of high-precision Earth gravity fields.Methods For the filtering of GOCE satellite gravity gradient data, we design two kinds of low-frequency systematic error filters based on moving average (MA) method and circle per revolution (CPR) empirical parameter method, and then combine them with the colored noise filter based on auto-regressive moving average (ARMA) model to form cascaded filters. To analyze the actual effect of these filters, the gravity field models GOGR-MA (MA+ARMA cascade filtering) and GOGR-CPR (CPR+ARMA cascade filtering) gravity field models of order 224 are inversed by using 70 d of GOCE data based on the space-wise least square method, respectively. The inversion models are compared with the first generation of GOCE series models, and GOCE and GRACE combined model GOCO01S solved by contemporaneous data. Then the model is checked by high-precision global positioning system leveling data and mean dynamic topography model.Results The results show that both the MA filter and CPR empirical parameter filter can weaken the influence of low-frequency system errors, and the latter is better. Moreover, ARMA decorrelation filter can effectively whiten the colored noise of gravity gradient data. The combined models GOGR-MA and GOGR-CPR have similar accuracy and they perform better than GOCO01S model developed by European Space Agency (ESA) using the same period data.Conclusions The inversion results verify the correctness and effectiveness of the proposed cascaded filters, which can provide a reference for filtering of GOCE gravity gradient data and inversion of high-precision gravity field. -
数据融合最初仅仅局限于军事领域，是针对异质传感器产生的海量数据进行一致性估计的方法[1-3]。20世纪90年代后，数据融合的概念被引入中国[4]。文献[5]将数据融合的结构分为3层，即数据层、特征层和决策层，在不同层进行数据融合会得到不同的结果，其差异体现在融合精度和处理效率等方面，此结构也作为基本的数据融合结构应用于其他领域。
在大地测量领域，文献[6]从函数模型误差和随机模型误差两方面详细讨论了不同的融合结果，具有较大的理论意义，相关的数据融合应用也较为广泛，体现在GPS与INS组合导航[7]、GPS与InSAR[8]、GPS不同系统间的融合[9]、遥感图像及地图数据融合[9, 10]以及3S集成[11]等方面。在工业测量中，相关的数据融合理论方法并不常见，文献[12]介绍了在盾构自动导向中综合利用全站仪和倾斜仪数据联合平差解算盾构姿态的方法，体现了采用多源异质传感器进行数据融合的优势。
在mm级精度的工业测量中，是采用全站仪进行特定点的三维坐标的采集，全站仪是测距仪(获取斜距)和经纬仪(获取角度)的集成，其获取的坐标实际是测距数据和测角数据的融合结果；而在船舶、飞机制造等工业测量领域，由于一些大型构件无法一站式获取全部待测点位三维信息，需采用移站测量(moving station method, MSM)的方法[13]，通过在不同测站上重复测量公共点，将不同测站的坐标系统一起来。但是，随着移站次数的增加，测量误差会不断积累，影响整体测量的精度。
本文以工业测量整体平差(overall adjustment, OA)为例，讨论了两种不同数据融合方法，即基于原始观测值(水平角、竖直角、斜距)的平差模型和基于坐标观测值(直角坐标系中的三维坐标)的平差模型，从解算精度、多余观测值的敏感性以及随机模型偏差对结果的影响等3方面，详细对比了两者的区别并和传统的移站测量结果进行比较，得到在工业测量中融合策略选取的相关结论。
1 工业测量整体平差函数模型
1.1 基于原始观测值的模型
在传统的工业测量中，采用全站仪获取数据是使用最多的方法之一。在每站测量前经过整平后，全站仪的原始观测值为测站点到待测点的水平方向、竖直角和斜距，这3个观测值互相独立，并在全站仪出厂时会给出角度和距离测量的标称精度。通常所有的待测点无法在一个测站全部观测，故采用移站测量的方法将不同测站的坐标系通过重复测量公共点进行统一，一般都以起始测站的坐标系作为全局坐标系。
在基于原始观测值的整体平差函数模型中，设第j测站在全局坐标系中的坐标为$ \left( {{{\hat{X}}}_{\text{stat}-j}}, \ \ {{{\hat{Y}}}_{\text{stat}-j}}, \ \ {{{\hat{Z}}}_{\text{stat}-j}} \right) $，测量了n个全局坐标系中的待测点$ \left( {{{\hat{X}}}_{i}}, \ \ {{{\hat{Y}}}_{i}}, \ \ {{{\hat{Z}}}_{i}} \right) $，其中i=1, 2, …, n，则有：
$$ \left\{ \begin{array}{l} {L_{\rm{s}}} + {V_{\rm{s}}} = \sqrt {{{\left( {{{\hat X}_i} - {{\hat X}_{{\rm{stat}} - j}}} \right)}^2} + {{\left( {{{\hat Y}_i} - {{\hat Y}_{{\rm{stat}} - j}}} \right)}^2} + {{\left( {{{\hat Z}_i} - {{\hat Z}_{{\rm{stat}} - j}}} \right)}^2}} \\ {L_{\rm{H}}} + {V_{\rm{H}}} = \arctan \frac{{{{\hat Y}_i} - {{\hat Y}_{{\rm{stat}} - j}}}}{{{{\hat X}_i} - {{\hat X}_{{\rm{stat}} - j}}}} - {{\hat A}_{{\rm{azi}} - j}}\\ {L_{\rm{V}}} + {V_{\rm{V}}} = \arctan \frac{{{{\hat Z}_i} - {{\hat Z}_{{\rm{stat}} - j}}}}{{\sqrt {{{\left( {{{\hat X}_i} - {{\hat X}_{{\rm{stat}} - j}}} \right)}^2} + {{\left( {{{\hat Y}_i} - {{\hat Y}_{{\rm{stat}} - j}}} \right)}^2}} }} \end{array} \right. $$ (1) 式中，Ls、LH、LV分别为斜距、水平角、竖直角观测值；对应Vs、VH、VV为观测值的改正数；$ {{\hat{A}}_{\text{azi}-j}} $为第j测站在全局坐标系中水平角初始零方向参数。将式(1)线性化后有：
$$ \left\{ \begin{array}{l} {V_{\rm{s}}} = {p_1}\delta {{\hat X}_{{\rm{stat}} - j}} + {p_2}\delta {{\hat Y}_{{\rm{stat}} - j}} + {p_3}\delta {{\hat Z}_{{\rm{stat}} - j}} - {p_1}\delta {{\hat X}_i} + {p_2}\delta {{\hat Y}_i} + {p_3}\delta {{\hat Z}_i} - {l_s}\\ {V_H} = - {p_4}\delta {{\hat X}_{{\rm{stat}} - j}} + {p_5}\delta {{\hat Y}_{{\rm{stat}} - j}} + {p_4}\delta {{\hat X}_i} - {p_5}\delta {{\hat Y}_i} - \delta {{\hat A}_{{\rm{az}} - j}} - {l_H}\\ {V_V} = {p_6}\delta {{\hat X}_{{\rm{stat}} - j}} + {p_7}\delta {{\hat Y}_{{\rm{stat}} - j}} + {p_8}\delta {{\hat Z}_{{\rm{stat}} - j}} - {p_6}\delta {{\hat X}_i} + {p_7}\delta {{\hat Y}_i} + {p_8}\delta {{\hat Z}_i} - {l_V} \end{array} \right. $$ (2) $$ \left\{ \begin{array}{l} {p_1} = \frac{{{X_{{\rm{stat}} - {j_0}}} + {X_{{i_0}}}}}{{{S_{{i_0}}}}}\\ {p_2} = \frac{{{Y_{{\rm{stat}} - {j_0}}} + {Y_{{i_0}}}}}{{{S_{{i_0}}}}}\\ {p_3} = \frac{{{Z_{{\rm{stat}} - {j_0}}} + {Z_{{i_0}}}}}{{{S_{{i_0}}}}}\\ {p_4} = \frac{{{Y_{{\rm{stat}} - {j_0}}} + {Y_{{i_0}}}}}{{D_{{i_0}}^2}}\\ {p_5} = \frac{{{X_{{\rm{stat}} - {j_0}}} + {X_{{i_0}}}}}{{D_{{i_0}}^2}}\\ {p_6} = \frac{{\left( {{Z_{{\rm{stat}} - {j_0}}} - {Z_{{i_0}}}} \right) \cdot \left( {{X_{{\rm{stat}} - {j_0}}} - {X_{{i_0}}}} \right)}}{{{D_{{i_0}}} \cdot \left( {D_{{i_0}}^2 + {{\left( {{Z_{{\rm{stat}} - {j_0}}} - {Z_{{i_0}}}} \right)}^2}} \right)}}\\ {p_7} = = \frac{{\left( {{Z_{{\rm{stat}} - {j_0}}} - {Z_{{i_0}}}} \right) \cdot \left( {{Y_{{\rm{stat}} - {j_0}}} - {Y_{{i_0}}}} \right)}}{{{D_{{i_0}}} \cdot \left( {D_{{i_0}}^2 + {{\left( {{Z_{{\rm{stat}} - {j_0}}} - {Z_{{i_0}}}} \right)}^2}} \right)}}\\ {p_8} = \frac{{ - D_{{i_0}}^2}}{{{D_{{i_0}}} \cdot \left( {D_{{i_0}}^2 + {{\left( {{Z_{{\rm{stat}} - {j_0}}} - {Z_{{i_0}}}} \right)}^2}} \right)}}\\ {l_{\rm{s}}} = {L_{\rm{s}}} - {S_{{i_0}}};\\ {l_H} = {L_H} + {A_0} - {A_1}\\ {l_V} = {L_V} - \arctan \frac{{{Z_i} - {Z_{{\rm{stat}} - {j_0}}}}}{{{D_{{i_0}}}}} \end{array} \right. $$ (3) 式中，Xstat-j0、Ystat-j0、Zstat-j0为在全局坐标系中第j测站坐标的近似值；Xi0、Yi0、Zi0为在全局坐标系中第i目标点的坐标近似值；A0为初始零方向近似值；A1为通过测站和目标点近似值解算的初始零方向，两者的范围都为0°~360°；Si0和Di0分别为通过测站点和目标点近似值求得的点间距及其水平投影。由于式(1)为非线性函数，求解通过最小二乘法迭代，这些近似值在不断迭代的过程中经求解的改正数调整，直到求解的改正数小于某一给定的阈值δ。
在某次测量任务中，假设共测了s个测站，目标点一共有t个，以第一站的坐标系作为全局坐标系，则需要解得参数共有(s－1)×3+t×3+(s－1)个。需要注意的是，在t个目标点中某些点是可以在多个测站中重复观测的，也就是在移站测量中的公共点。参数解为：
$$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{\hat X}} = }\\ {\left( {\underbrace {\begin{array}{*{20}{c}} {\delta {{\hat X}_{{\rm{stat}} - 2}}}&{\delta {{\hat Y}_{{\rm{stat}} - 2}}}&{\delta {{\hat Z}_{{\rm{stat}} - 2}}}& \cdots &{\delta {{\hat X}_{{\rm{stat}} - s}}}&{\delta {{\hat Y}_{{\rm{stat}} - s}}}&{\delta {{\hat Z}_{{\rm{stat}} - s}}} \end{array}}_{s - 1个测站坐标}} \right.}\\ {\underbrace {\begin{array}{*{20}{c}} {\delta {{\hat X}_1}}&{\delta {{\hat Y}_1}}&{\delta {{\hat Z}_1}}& \cdots &{\delta {{\hat X}_t}}&{\delta {{\hat Y}_t}}&{\delta {{\hat Z}_t}} \end{array}}_{t\;个目标点坐标}}\\ {{{\left. {\underbrace {\begin{array}{*{20}{c}} {\delta {{\hat a}_{{\rm{azi}} - 2}}}& \cdots &{\delta {{\hat a}_{{\rm{azi}} - s}}} \end{array}}_{s - 1个测站初始零方向参数}} \right)}^{\rm{T}}}} \end{array} $$ (4) 由于第一个测站的测站坐标和初始零方向都为0，因此该测站列立的误差方程仅有目标点的坐标参数。
观测方程的随机模型由仪器标定精度给定，测距精度σs=(a+b×Ls), a为固定误差，b为比例误差, Ls为测站点到目标点的斜距；由于一般测量目标点的角度值都是半测回测一次的结果，而标称精度一般为一测回方向的精度，所以实际的测角精度σangle为给定测角精度的$ \sqrt{2} $倍，在空间上带来的误差影响为σangle×Ls。此外，距离观测与角度观测之间是独立的，因此基于原始观测值的整体平差权阵为对角阵。对于一次观测目标点的3个原始观测值权阵为：
$$ {\mathit{\boldsymbol{P}}_1} = {\rm{diag}}\left[ {\begin{array}{*{20}{c}} {\frac{{\sigma _0^2}}{{\sigma _s^2}}}&{\frac{{\sigma _0^2}}{{{{\left( {{\sigma _{{\rm{angle}}}} \times {L_s}} \right)}^2}}}}&{\frac{{\sigma _0^2}}{{{{\left( {{\sigma _{{\rm{angle}}}} \times {L_s}} \right)}^2}}}} \end{array}} \right] $$ (5) 式中，σ0为给定的验前单位权中误差，一般设定为1 mm。验后单位权中误差由式(6)求得：
$$ \hat \sigma _0^2 = \frac{{{\mathit{\boldsymbol{V}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_1}\mathit{\boldsymbol{V}}}}{{n - \left( {s - 1} \right) \times 4 - t \times 3}} $$ (6) 式中，n为总的观测值个数; V为式(2)中的残差向量，则求解参数$ \mathit{\boldsymbol{\hat{X}}} $的精度：
$$ {\mathit{\boldsymbol{D}}_{\hat X\hat X}} = \hat \sigma _0^2{\mathit{\boldsymbol{Q}}_{\hat X\hat X}} = \hat \sigma _0^2{\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_1}\mathit{\boldsymbol{B}}} \right)^{ - 1}} $$ (7) 式中，矩阵B为式(2)中的系数阵。
1.2 基于坐标观测值的模型
基于坐标观测值的整体平差中，在第j测站，观测i号点的函数模型为：
$$ \left( {\begin{array}{*{20}{c}} {{L_X}}\\ {{L_Y}}\\ {{L_Z}} \end{array}} \right) + \left( {\begin{array}{*{20}{c}} {{V_X}}\\ {{V_Y}}\\ {{V_Z}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {\Delta {{\hat X}_j}}\\ {\Delta {{\hat Y}_j}}\\ {\Delta {{\hat Z}_j}} \end{array}} \right) + \left( {\begin{array}{*{20}{c}} {\cos {{\hat \gamma }_j}}&{ - \sin {{\hat \gamma }_j}}&0\\ {\sin {{\hat \gamma }_j}}&{\cos {{\hat \gamma }_j}}&0\\ 0&0&1 \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{{\hat X}_i}}\\ {{{\hat Y}_i}}\\ {{{\hat Z}_i}} \end{array}} \right) $$ (8) 式中，$ \Delta {{\hat{X}}_{j}} $、$ \Delta {{\hat{Y}}_{j}} $、$ \Delta {{\hat{Z}}_{j}} $、$ {\hat{\gamma }}_{j} $为第j测站相对于全局坐标系的平移参数和绕Z轴的旋转角；$ {{\hat{X}}_{i}} $、$ {{\hat{Y}}_{i}} $、$ {{\hat{Z}}_{i}} $为i号待测点在全局坐标系中的三维坐标值。
式(8)中，求解站与站之间转换参数采用的观测值为坐标，而坐标值是由笛卡尔坐标系中原始观测的斜距、水平角和竖直角推算求得的：
$$ \left\{ \begin{array}{l} {L_X} = {L_{\rm{s}}} \cdot \cos {L_H} \cdot \sin {L_V}\\ {L_Y} = {L_{\rm{s}}} \cdot \sin {L_H} \cdot \sin {L_V}\\ {L_Z} = {L_{\rm{s}}} \cdot \cos {L_V} \end{array} \right. $$ (9) 可以看出，虽然原始观测值之间是独立的，但是坐标观测值是相关的，由误差传播可知，坐标观测值的方差阵为：
$$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{D}}_{XYZ}} = \mathit{\boldsymbol{A}} \cdot {\rm{diag}}}\\ {\left[ {\begin{array}{*{20}{c}} {\sigma _s^2}&{{{\left( {{\sigma _{{\rm{angle}}}} \times {L_s}} \right)}^2}}&{{{\left( {{\sigma _{{\rm{angle}}}} \times {L_s}} \right)}^2}} \end{array}} \right] \cdot {\mathit{\boldsymbol{A}}^{\rm{T}}}} \end{array} $$ (10) 式中，
$$ \mathit{\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}} {\cos {L_H}\sin {L_V}}&{\sin {L_H}\sin {L_V}}&{\cos {L_H}\cos {L_V}}\\ {\sin {L_H}\sin {L_V}}&{\cos {L_H}\sin {L_V}}&{\sin {L_H}\cos {L_V}}\\ {\cos {L_V}}&0&{\sin {L_V}} \end{array}} \right] $$ (11) 仍设定验前单位权中误差为σ0，则对应权阵为：
$$ {\mathit{\boldsymbol{P}}_2} = {\left( {\frac{{{\mathit{\boldsymbol{D}}_{XYZ}}}}{{\sigma _0^2}}} \right)^{ - 1}} $$ (12) 与式(1)一样，式(8)也是非线性的函数，需进行线性化，可以参考大角度七参数三维基准转换模型[14]。需要注意的是，由于以第一测站的坐标系作为全局坐标系，其误差方程更为简单：
$$ \left( {\begin{array}{*{20}{c}} {{L_X}}\\ {{L_Y}}\\ {{L_Z}} \end{array}} \right) + \left( {\begin{array}{*{20}{c}} {{V_X}}\\ {{V_Y}}\\ {{V_Z}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {{{\hat X}_i}}\\ {{{\hat Y}_i}}\\ {{{\hat Z}_i}} \end{array}} \right) $$ (13) 传统移站测量方法中，站与站之间通过观测公共点将坐标系统一并传递下去。实际公共点在两套坐标系的观测值中都含有误差，而实际上一般仅采用最小二乘方法处理，忽略了系数阵的误差，因此随着移站的增多，误差累积越来越大。
对比§1.1和§1.2可以看出，两者采用不同的函数模型，但最终都能通过整体平差求解出待测点在第一测站坐标系下的坐标，而且没有传统移站测量中出现的误差累积的问题；此外，两者的系数阵中都不存在含观测误差的观测值，因此通过传统最小二乘算法和准确的随机模型即可获得准确的参数估计，而避免采用整体最小二乘等较复杂的方法来提高解算精度。
2 算例与分析
为了对比研究两种数据融合模型的优劣，笔者利用索佳Net05高精度全站仪，围绕某旗杆墩架设了4个测站，并在墩上的4个面观测了24个贴有反射片的待测点，在每个测站上都可以观测到两个墩面共12个目标点，测站1观测1~12号点，测站2观测7~18号点，测站3观测13~24号点，测站4观测19~24号点和1~6号点，如图 1所示。
对于索佳Net05，当反射物为反射片时，测距标称精度为σs = 0.5+1× ppm mm，测角精度为0.5″，则半测回某方向的精度为σangle=0.5×$ \sqrt{2} $≈0.7″，采用§1.1中的融合方法及式(7)的精度评定公式，以第一测站坐标系作为全局坐标系，得到该坐标系下除测站1的其他3个测站的测站坐标及初始零方向参数和24个待测点的坐标，解算结果见表 1、表 2。同理可以得到按照§1.2的融合方法的解算结果，见表 3、表 4。
表 1 基于原始观测值的整体平差结果Table 1. Result of Overall Adjustment Based on Original Data类型 编号 X/m Y/m Z/m $ \hat{\sigma } $X/mm $ \hat{\sigma } $Y/mm $ \hat{\sigma } $Z/mm 待测点 1 -1.847 0 -11.186 5 0.764 1 0.2 0.2 0.1 2 -1.483 9 -10.830 7 0.779 2 0.2 0.1 0.1 3 -1.064 8 -10.431 0 0.796 8 0.1 0.2 0.1 4 -1.872 2 -11.153 5 0.230 3 0.2 0.2 0.1 5 -1.521 4 -10.806 9 0.196 2 0.2 0.1 0.1 6 -1.098 7 -10.396 3 0.173 7 0.1 0.2 0.1 7 -0.786 9 -10.499 8 0.823 1 0.1 0.2 0.1 8 -0.454 7 -10.842 2 0.828 7 0.2 0.2 0.2 9 -0.108 5 -11.190 8 0.835 5 0.2 0.2 0.2 10 -0.768 9 -10.449 1 0.260 0 0.1 0.2 0.1 11 -0.417 1 -10.811 3 0.272 4 0.2 0.2 0.1 12 -0.082 2 -11.157 0 0.318 2 0.2 0.2 0.2 13 -0.093 4 -11.504 2 0.759 5 0.2 0.2 0.2 14 -0.439 7 -11.837 9 0.779 4 0.3 0.2 0.2 15 -0.826 0 -12.211 0 0.767 9 0.3 0.2 0.2 16 -0.054 3 -11.531 9 0.179 0 0.2 0.2 0.2 17 -0.442 7 -11.904 3 0.191 9 0.3 0.2 0.2 18 -0.798 3 -12.253 4 0.168 8 0.3 0.2 0.2 19 -1.116 3 -12.211 5 0.733 2 0.3 0.2 0.1 20 -1.468 5 -11.854 8 0.745 2 0.3 0.2 0.1 21 -1.845 4 -11.474 2 0.735 1 0.2 0.2 0.1 22 -1.150 9 -12.234 0 0.235 4 0.3 0.2 0.1 23 -1.493 8 -11.888 1 0.222 7 0.3 0.2 0.1 24 -1.877 1 -11.502 8 0.192 6 0.2 0.2 0.1 2 11.670 7 -12.423 4 -0.018 8 0.3 2.5 0.1 测站 3 -1.140 5 -23.007 8 -0.013 3 2.1 0.1 0.1 4 -11.095 8 -11.317 5 0.039 3 0.1 1.8 0.1 表 2 测站其他参数平差结果Table 2. Result of Other Parament编号 初始零方向
/rad解算精度/(″) 验后单位权中误差/mm 2 -2.922 588 41.6 0.5 3 0.741 673 37.4 4 5.704 257 37.9 表 3 基于坐标观测值的整体平差结果Table 3. Result of Overall Adjustment Based on Coordinate Data类型 编号 X/m Y/m Z/m $ \hat{\sigma } $X/mm $ \hat{\sigma } $Y/mm $ \hat{\sigma } $Z/mm 待测点 1 -1.847 1 -11.186 5 0.764 1 0.1 0.4 0.0 2 -1.483 8 -10.831 0 0.779 3 0.1 0.3 0.1 3 -1.064 9 -10.431 4 0.796 8 0.1 0.3 0.1 4 -1.872 2 -11.153 6 0.230 3 0.1 0.4 0.0 5 -1.521 4 -10.807 0 0.196 2 0.1 0.3 0.0 6 -1.098 7 -10.396 7 0.173 7 0.1 0.3 0.0 7 -0.786 8 -10.499 2 0.823 0 0.1 0.3 0.1 8 -0.455 0 -10.841 7 0.828 7 0.1 0.3 0.1 9 -0.108 4 -11.190 4 0.835 5 0.1 0.4 0.1 10 -0.768 7 -10.448 7 0.260 0 0.1 0.3 0.1 11 -0.417 4 -10.811 0 0.272 4 0.1 0.3 0.1 12 -0.082 0 -11.156 6 0.318 2 0.1 0.4 0.1 13 -0.093 3 -11.504 0 0.759 5 0.5 0.4 0.1 14 -0.439 8 -11.837 5 0.779 5 0.5 0.3 0.1 15 -0.825 9 -12.210 5 0.767 9 0.6 0.3 0.1 16 -0.054 5 -11.531 7 0.179 1 0.5 0.4 0.1 17 -0.442 8 -11.903 8 0.191 9 0.5 0.3 0.1 18 -0.798 1 -12.252 8 0.168 8 0.6 0.3 0.1 19 -1.116 1 -12.211 9 0.733 2 0.6 0.3 0.1 20 -1.468 5 -11.855 1 0.745 2 0.5 0.3 0.1 21 -1.845 6 -11.474 2 0.735 1 0.5 0.4 0.1 22 -1.150 8 -12.234 4 0.235 4 0.6 0.3 0.1 23 -1.494 0 -11.888 4 0.222 7 0.5 0.3 0.1 24 -1.877 4 -11.502 9 0.192 6 0.5 0.4 0.1 2 8.697 1 -14.663 0 0.018 8 5.2 1.5 0.0 测站 3 16.385 9 16.191 1 0.013 3 3.7 4.0 0.0 4 3.090 5 15.541 9 -0.039 3 4.4 2.2 0.0 表 4 测站其他参数平差结果Table 4. Result of Other Parament编号 绕Z轴旋转角/rad 解算精度/(″) 验后单位权中误差/mm 2 2.922 953 100.6 3 -0.742 009 98.3 1.7 4 0.579 411 91.0 从表 1、表 2可以看出，待测点的点位解算精度可达0.3 mm左右，且精度分布均匀。从表 3、表 4可以看出，点位解算精度基本在0.5 mm左右，且测站的定位参数解算精度也不如表 1、表 2的结果高。虽然从数学角度看，§1.2的融合方法中，观测值是由§1.1的原始观测值等价推导得到的，随机模型也是从§1.1的误差模型通过误差传播得来的，但是融合结果却不是完全等价的。
按照传统移站测量的思路，采用§1.1的函数模型进行分站平差，即第一测站观测的1~12点都为全局坐标系的坐标；第二测站观测7~12点来解算测站参数，根据解算得到的测站参数和13~18点的原始观测值解算出13~18点的全局坐标系的坐标；同理，第三测站也可求解出19~24点的全局坐标系坐标。此方法称为基于原始观测值的移站测量法。同理采用§1.2的函数模型进行传统移站测量的方法称为基于坐标观测值的移站测量法。
通过不同的方法得到24个待测点的全局坐标，有如下4种方案：①基于原始观测值的整体平差法；②基于坐标观测值的整体平差法；③基于原始观测值的移站测量法；④基于坐标观测值的移站测量法。分别对比表 1、表 2和表 3、表 4，可以看出方案①的参数解算精度较高。因此，以方案①解算的待测点坐标为基准，可以得到方案①与方案②、方案①与方案③以及方案②与方案④的坐标点位偏差，如图 2所示。
可以看出，方案①与方案②的各点点位偏差较均匀，范围在0.1~0.6 mm之间；而在传统的MSM中，无论是基于原始观测值或是坐标观测值，其解算的点位偏差随着搬站的增加而递增。
图 2中MSM采用6个公共点进行搬站，在经过2次搬站后测得的19~24号点与整体平差的点位偏差可达3~4 mm左右，而实际测量中一般不会采用这么多公共点进行移站测量。如若在24个待测点中选取12个点，如选取1、3、5、…、19、21、23号待测点，每两个测站之间只有3个公共点，即减少了多余观测数，可以得到两种OA结果与对应传统MSM的点位偏差，如图 3所示。
![]()
图 3 3个公共点时整体平差与移站测量的点位偏差Figure 3. Coordinate Differences Between OA and MSM with 3 Common Points从图 2和图 3可以看出，基于原始观测值的OA和MSM受多余观测值数量的影响较小，而基于坐标观测值的OA和MSM受多余观测数的影响非常明显。当搬站公共点仅有3个的时候，经两次搬站后测量的目标点与整体平差点位偏差可达将近3 cm，因此基于坐标观测值的MSM的点位偏差更加难以控制。
为了进一步比较当随机模型与实测数据存在偏差，即存在随机模型误差时，方案①和方案②的解算稳定性，在所有斜距原始观测值分别添加中误差为1 mm和5 mm的白噪声(标称精度仅为0.5 mm)，而随机模型仍然采用仪器的标称精度计算，可以得到两种方案在相同额外白噪声影响下解算的坐标偏差，如图 4所示。
![]()
图 4 添加不同白噪声对两种解算点位偏差影响Figure 4. Influence on the Coordinate Deviation when Adding Different White Noise方案①在外加1 mm斜距误差后解算的24个待测点位与不加误差时偏差在0.2~1.5 mm，在外加5 mm斜距误差后解算的24个待测点位与不加误差时偏差在1.1~5.9 mm；同样的条件，方案②的偏差分别在0~0.8 mm和1.1~2.9 mm。从以上偏差的范围和稳定性都可以看出，当随机模型含有误差时，方案②参数解算稳定性更高。
3 结语
1) 当随机模型的误差不明显，即仪器实际测量精度和标称精度相差不大时，采用基于原始观测值的整体平差融合算法来解算待定点坐标的方法精度最高。在传统的移站测量解算中，采用基于原始观测值的方法也比基于坐标观测值的结果要稳定。
2) 当多余观测数较少，即公共点个数较少时，若采用传统移站测量法解算点位坐标，不宜采用基于坐标观测值的融合方法，采用基于原始观测值的融合方法会更可靠。
3) 当随机模型含有明显误差时，基于坐标观测值的整体平差法的融合结果要比基于原始观测值的整体平差融合结果稳定可靠。
-
![]()
图 1 GOCE卫星重力梯度$ {V}_{zz} $分量的误差PSD
Figure 1. Error PSD of Gravity Gradient Component $ {V}_{zz} $ of GOCE Satellite
![]()
图 2 MA滤波器处理后的$ {V}_{zz} $分量误差PSD
Figure 2. Error PSD of $ {V}_{zz} $ Processed by MA Filter
![]()
图 3 CPR滤波器处理后的$ {V}_{zz} $误差PSD
Figure 3. Error PSD of $ {V}_{zz} $ Processed by CPR Filter
![]()
图 4 MA+ARMA级联滤波器处理后4个重力梯度各分量的误差PSD
Figure 4. Error PSD of 4 Gravity Gradient Components Processed by MA+ARMA Cascaded Filter
![]()
图 5 CPR+ARMA级联滤波器处理后4个重力梯度分量的误差PSD
Figure 5. Error PSD of 4 Gravity Gradient Components Processed by CPR+ARMA Cascaded Filter
![]()
图 6 GOCE-MA和GOCE-CPR相比GOCO06S模型的阶误差RMS
Figure 6. Degree-Error RMS of GOCE-MA and GOCE-CPR Compared with GOCO06S
![]()
图 7 6种模型相比GOCO06S模型的阶误差RMS
Figure 7. Degree-Error RMS of 6 Models Compared with GOCO06S
![]()
图 8 不同模型计算的黑潮表层地转流
Figure 8. Surface Geostrophic Current of Kuroshio Calculated by Different Models
表 1 全球$ 1°\times 1° $格网大地水准面和重力异常误差
Table 1 Statistical Results of Global $ 1°\times 1° $ Grid Geoid Errors and Gravity Anomaly Errors
阶次 模型 大地水准面误差/m 重力异常误差/mGal 平均值 RMS 平均值 RMS 224 GOCO01S 0.003 0.208 0.046 6.510 224 GOGR-MA 0.001 0.155 0.001 4.664 224 GOGR-CPR 0.001 0.168 -0.003 4.980 224 DIR_R1 0.000 0.134 0.013 3.988 224 TIM_R1 -0.053 0.595 -0.359 8.436 210 SPW_R1 0.003 0.219 0.032 6.134 
下载: 导出CSV
表 2 GPS水准外部检核统计结果
Table 2 External Checking Results of GPS Leveling
阶次 模型 美国 欧洲 中国西南地区 平均值/m 标准差/m 平均值/m 标准差/m 平均值/m 标准差/m 224 GOCO01S 1.039 0.525 -0.030 0.296 0.622 0.783 224 GOGR-MA 1.049 0.515 -0.028 0.278 0.613 0.769 224 GOGR-CPR 1.042 0.510 -0.028 0.291 0.670 0.777 224 DIR_R1 1.036 0.509 0.000 0.274 0.633 0.691 224 TIM_R1 1.038 0.529 -0.015 0.299 0.601 0.778 210 SPW_R1 1.034 0.544 -0.028 0.294 0.589 0.749
下载: 导出CSV
-
[1] Floberghagen R, Fehringer M, Lamarre D, et al. Mission Design, Operation and Exploitation of the Gravity Field and Steady-State Ocean Circulation Explorer Mission[J]. Journal of Geodesy, 2011, 85(11): 749-758. doi: 10.1007/s00190-011-0498-3
[2] Pail R, Bruinsma S, Migliaccio F, et al. First GOCE Gravity Field Models Derived by Three Different Approaches[J]. Journal of Geodesy, 2011, 85(11): 819-843. doi: 10.1007/s00190-011-0467-x
[3] Yu J H, Wan X Y. Recovery of the Gravity Field from GOCE Data by Using the Invariants of Gradient Tensor[J]. Science China Earth Sciences, 2013, 56(7): 1193-1199. doi: 10.1007/s11430-012-4427-y
[4] Lu B, Luo Z C, Zhong B, et al. The Gravity Field Model IGGT_R1 Based on the Second Invariant of the GOCE Gravitational Gradient Tensor[J]. Journal of Geodesy, 2018, 92(5): 561-572. doi: 10.1007/s00190-017-1089-8
[5] 徐新禹, 李建成, 姜卫平, 等. 基于空域最小二乘法求解GOCE卫星重力场的模拟研究[J]. 测绘学报, 2011, 40(6): 697-702. https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201106006.htm Xu Xinyu, Li Jiancheng, Jiang Weiping, et al. Simulation Study for Recovering GOCE Satellite Gravity Model Based on Space-Wise LS Method[J]. Acta Geodaetica et Cartographica Sinica, 2011, 40(6): 697-702. https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201106006.htm
[6] 钟波, 宁津生, 罗志才, 等. 联合GOCE卫星轨道和重力梯度数据严密求解重力场的模拟研究[J]. 武汉大学学报(信息科学版), 2012, 37(10): 1215-1220. http://ch.whu.edu.cn/article/id/344 Zhong Bo, Ning Jinsheng, Luo Zhicai, et al. Simulation Study of Rigorous Gravity Field Recovery by Combining GOCE Satellite Orbit and Gravity Gradient Data[J]. Geomatics and Information Science of Wuhan University, 2012, 37(10): 1215-1220. http://ch.whu.edu.cn/article/id/344
[7] 钟波, 罗志才, 李建成, 等. 联合高低卫-卫跟踪和卫星重力梯度数据恢复地球重力场的谱组合法[J]. 测绘学报, 2012, 41(5): 735-742. https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201205020.htm Zhong Bo, Luo Zhicai, Li Jiancheng, et al. Spectral Combination Method for Recovering the Earth’s Gravity Field from High-Low SST and SGG Data[J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(5): 735-742. https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201205020.htm
[8] Rummel R, Yi W Y, Stummer C. GOCE Gravitational Gradiometry[J]. Journal of Geodesy, 2011, 85(11): 777. doi: 10.1007/s00190-011-0500-0
[9] 万晓云, 于锦海, 曾艳艳. GOCE引力梯度的频谱分析及滤波[J]. 地球物理学报, 2012, 55(9): 2909-2916. https://www.cnki.com.cn/Article/CJFDTOTAL-DQWX201209011.htm Wan Xiaoyun, Yu Jinhai, Zeng Yanyan. Frequency Analysis and Filtering Processing of Gravity Gradients Data from GOCE[J]. Chinese Journal of Geophysics, 2012, 55(9): 2909-2916. https://www.cnki.com.cn/Article/CJFDTOTAL-DQWX201209011.htm
[10] 苏勇, 范东明, 蒲星钢, 等. 联合GOCE卫星数据和GRACE法方程确定SWJTU-GOGR01S全球重力场模型[J]. 武汉大学学报(信息科学版), 2018, 43(3): 457-463. doi: 10.13203/j.whugis20150100 Su Yong, Fan Dongming, Pu Xinggang, et al. New Static Gravity Field Model SWJTU-GOGR01S Derived from GOCE Data and GRACE Normal Equation[J]. Geomatics and Information Science of Wuhan University, 2018, 43(3): 457-463. doi: 10.13203/j.whugis20150100
[11] Pail R, Goiginger H, Schuh W D, et al. Combined Satellite Gravity Field Model GOCO01S Derived from GOCE and GRACE[J]. Geophysical Research Letters, 2010, 37(20), DOI: 10.1029/2010GL044906.
[12] 罗志才, 钟波, 宁津生. 卫星重力梯度测量确定地球重力场的理论与方法[M]. 武汉: 武汉大学出版社, 2015. Luo Zhicai, Zhong Bo, Ning Jinsheng. Theory and Method for Determining the Earth's Gravity Field from Satellite Gravity Gradiometry[M]. Wuhan: Wuhan University Press, 2015.
[13] Koch K, Kusche J. Regularization of Geopotential Determination from Satellite Data by Variance Components[J]. Journal of Geodesy, 2002, 76(5): 259-268.
[14] Hashemi F H, Ditmar P, Klees R, et al. Validation of Static Gravity Field Models Using GRACE K-Band Ranging and GOCE Gradiometry Data[J]. Geophysical Journal International, 2013, 194(2): 751-771.
[15] Schuh W D. The Processing of Band-Limited Measurements; Filtering Techniques in the Least Squares Context and in the Presence of Data Gaps[J]. Space Science Reviews, 2003, 108(1/2): 67-78.
[16] Mayer-Gürr T, Kurtenbach E, Eicker A. The Satellite-Only Gravity Field Model ITG-GRACE 2010s[EB/OL]. [2020-10-20]. http://www.igg.uni-bonn.de/apmg/index.php?id=itg-grace2010, 2010.
[17] Kvas A, Mayer G T, Krauss S, et al. The Satellite-Only Gravity Field Model GOCO06S[EB/OL]. [2020-06-30]. https://doi.org/10.5880/ICGEM.2019.002.
[18] Schaeffer P, Faugére Y, Legeais J F, et al. The CNES_CLS11 Global Mean Sea Surface Computed from 16 Years of Satellite Altimeter Data[J]. Marine Geodesy, 2012, 35(sup1): 3-19.
[19] Dinardo S. GUT Tutorial Version 6.4[R]. ESA, 2012.
-
期刊类型引用(4)
1. 朱杰，郑加柱，陈红华，杨静，胡平昌，陆敏燕. 结合POI数据的南京市商业中心识别与集聚特征研究. 现代测绘. 2022(06): 34-39 . 
百度学术
2. 金澄，安晓亚，陈占龙，马啸川. 矢量居民地多边形多级图划分聚类方法. 武汉大学学报(信息科学版). 2021(01): 19-29 . 百度学术
3. 张铭龙，何贞铭. 基于因子分析法的城市商业中心抽取研究. 地理空间信息. 2021(08): 58-60+64+5 . 百度学术
4. 李卫东，张铭龙，段金龙. 基于POI数据的南京市空间格局定量研究. 世界地理研究. 2020(02): 317-326 . 百度学术
其他类型引用(3)
计量
- 文章访问数: 293
- HTML全文浏览量: 81
- PDF下载量: 61
- 被引次数: 7
