高精度高分辨率地球重力场观测信息在大地测量学、地球物理学、地球动力学、海洋学和空间科学等研究领域中具有重要支撑作用，在地质矿产资源勘查、远程空间飞行器轨道控制和相关军事活动保障等工程领域也具有重要的应用价值[1-4]．精化大地水准面是地球重力场观测信息最重要的应用领域之一，也始终是现代大地测量科学研究的核心问题．特别是，随着全球导航卫星系统(GNSS)高精度测高技术的发展，GNSS大地高“+”大地水准面模型的测高新模式，已经成为当今高程测量现代化的重点发展方向[5]．GNSS高精度定位技术实现了以参考椭球面为基准的大地高测量现代化以后，建立高精度高分辨率的大地水准面数字模型，就成为当前突破海拔高程测量现代化发展瓶颈的核心关键．对大地水准面模型不断精化的现实需求，正在持续推动地球重力场新型探测装备和技术手段的发展，因为除了计算理论方法和手段方面的因素外，大地水准面模型的最终精度主要取决于地球重力场观测数据资源的覆盖率和精细化程度．在海洋重力场信息观探测领域，除了具有全球海域探测能力的卫星测高重力反演手段外，海面船载和航空重力测量仍是目前获取中高频局部重力场信息最有效的两种技术手段，二者各具特点，具有较好的互补性．海面船载重力测量技术已经有近百年的发展历史，其理论方法体系较为成熟[3,6-7]．航空重力测量技术取得突破主要得益于GNSS动态精密定位技术的支持，解决了飞机载体运动加速度的高精度测定难题，同时得益于重力传感器和稳定平台两大核心技术的发展进步[8]．但与陆地静态重力测量模式相比，海洋和航空(以下简称海空)重力测量具有明显的动态效应特性，由此带来了观测数据分析与处理上的复杂性．具体表现在：一方面，受海上风、流、压及海水运动、海浪起伏、近地空间大气湍流等各类环境因素及飞机发动机引起机身震动的干扰，海空重力观测记录会同时受到由载体运动引起的水平和垂直方向干扰加速度、厄特弗斯(Eötvös)、交叉耦合(Cross-Coupling，CC)、杆臂(也称动态偏心)等动态效应的综合影响；另一方面，受重力传感器自身结构设计和制造工艺上的制约，海空重力观测记录还会受到来自仪器固有零点漂移、参数标定误差及元器件不稳定或老化引起的仪器掉格及其他异常突变等多种因素干扰[6-7,9]，因此海空重力测量具有比陆地测量更多更难以掌控的观测噪声干扰源．各类误差源的综合影响不仅会给海空重力观测带来一定大小的随机性误差，还会带来不可忽视的系统性偏差和粗差，这些误差干扰的变化幅度往往比待探测的地球重力加速度大百倍甚至千倍．尽管通过数字滤波和测线网平差手段，可在一定程度上减弱各类误差源的影响，但其剩余误差影响对海空重力测量成果的干扰仍是不可忽略的[10]．特殊的动态测量作业模式也给海空重力测量数据质量控制和评估带来较大的困难，一方面是海上和空中均缺少高精度的重力控制点信息，另一方面是海上和空中都不易实施严格意义上的重复点观测．实践中，一般采用测线交叉或重复线测量方法进行海空重力测量数据内部符合精度评估，采用同船同机观测或同区域不同航次观测数据比对方式进行外部符合精度评估[3,6]．显然，内符合精度检核方法很难发现海空重力测量中的系统偏差和粗差，外部质量检核方法的可行性和有效性主要取决于外部数据源的可用性和可靠性，但不管是同船同机观测，还是同区域不同航次观测数据比对，都存在质量检核成本和实施难度过高的问题，因此，外部质量检核方法一般只用于新型测量仪器技术性能指标的验证和评估，不宜作为常态化的海上作业观测数据质量检核和评估手段．为了突破这方面的限制，基于对国际上最新发布的卫星测高重力模型的分析和评估结果，提出将其中最具代表性的高分辨率测高重力数值模型作为外部数据源，引入海空重力测量数据质量控制领域，发挥其作为局部重力场趋势性变化的控制作用，通过点、线、面多种比对方式直接检测海空重力测量数据中的粗差和系统偏差．1 卫星测高重力模型分析评估卫星测高技术诞生于20世纪70年代，将该项技术获取的几何观测量转换为物理场参量可得到不同类别的测高反演重力模型，这是此项新技术的典型应用．经过几十年的发展和积累，卫星测高技术已经取得了一系列重要突破和显著进展，新型测高卫星陆续发射，有越来越多的卫星测高数据资源可供利用，国际上有诸多研究机构和学者依据不同时期获取的多代卫星测高数据集，先后反演得到了多个版本不同分辨率的卫星测高重力模型[3,11-14]．在这些众多的国际模型中，最具代表和影响力的模型当属由美国加州大学圣地亚哥分校(UCSD)斯克里普斯海洋研究所(SIO)Sandwell和Smith研究团队推出的SIO系列模型及由丹麦科技大学(DTU)Andersen和Knudsen研究团队研发的DTU系列模型．由于受卫星轨迹分布密度和测高数据质量的制约，早期发布的卫星测高重力模型精度和分辨率都较低，模型内插值与海面船测点值的互差均方根值一般在10～30 mGal(1 mGal=10-5 m/s2)[3]．随着新型测高数据的补充和积累，以及数据处理新方法的采用，卫星测高重力模型精度和分辨率均得到显著提升，测量点值比对精度很快提升到几个毫伽水平，模型标称分辨率达到2′×2′甚至1′×1′(实际分辨率不超过2′×2′)．文献[12]曾使用加利福尼亚湾海域的海面船测重力资料，对文献[11]、文献[13]和文献[14]等各自发布的卫星测高重力模型做过精度检核评估，3个模型与船测数据的比对精度依次为12.1，10.9和11.4 mGal；文献[14]发布的测高重力模型在全球12个不同海区，与船测数据的比对精度在5～14 mGal之间；文献[15]发布的DTU10模型在部分海区与船测重力的比对精度已经达到4 mGal．最近一个时期，卫星测高重力模型技术得到迅猛发展，不同版本的模型更新周期越来越短，模型精度越来越高，SIO系列模型已经推出2019年版的V28模型[16]，DTU系列模型也已推出2018年版的DTU17模型[17]．自2010年以后，卫星测高重力模型精度的改善主要得益于包括Jason-1/2，CryoSat-2和SARAL/AltiKa等在内的新型卫星测高资料的支持，特别是SARAL/AltiKa卫星观测数据的贡献．根据文献[16]，加入新型卫星测高资料以后，新版本的V28模型与墨西哥湾海区的船测重力数据(精度为0.5 mGal)作比较，在12 km滤波长度下，V28模型的比对精度达到了1.23 mGal，而前面两个版本V23和V18模型的比对精度分别为1.52 mGal和2.05 mGal．由此可见：新版本相对旧版本模型的精度改善幅度还是比较可观的．基于上述比对结果，同时依据每个5′×3′区块内所有测高卫星获取的沿轨海面倾斜偏差的中位平均值及对应区块内的测点数量，文献[16]给出了V28模型在全球海域的不确定度分布图，其量值变化幅度在0～4 mGal之间．文献[17]分别在大西洋西北部和北冰洋海域开展了DTU系列模型的精度分析和评估，其中，DTU10，DTU13，DTU15和DTU17等4个模型与大西洋西北部超过1.4×106个船测重力点值(由美国国家地理空间情报局(NGA)提供，精度为1.5～2.0 mGal)的互比均方根差值分别为3.16，2.83，2.51和2.51 mGal，与北冰洋海域超过5.5×104个航空重力测量点值(来自国际航空重力测量项目LomGrav，精度优于2 mGal)的互比均方根差值分别为8.81，5.91，5.45和3.78 mGal．从这些有限的分析评估结果可以看出，国际上近期推出的卫星测高重力模型精度在大部分海区已经优于4 mGal．为了进一步了解掌握卫星测高重力模型的可靠性和可用性，分别在中国近海、南海和西北太平洋3个海域各挑选一个区块的海面船测重力数据作为基准，对V28模型进行比对分析和精度评估，3个海域对应的区块简称为区块1，2和3，3个区块船测重力数据的观测精度在1～2 mGal之间，模型值与观测重力点值的互比统计结果如表1所示．作为比较，表1同时列出了地球重力位模型EGM2008与观测点值的互比统计结果．10.13245/j.hust.220918.T001表1模型值与观测重力点值的互比统计结果参量区块1区块2区块3EGM2008V28EGM2008V28EGM2008V28最大差/mGal13.4411.1252.6327.7349.1344.52最小差/mGal-14.07-9.50-23.19-17.91-36.66-39.21平均差/mGal-0.03-0.28-0.34-0.33-1.60-0.54均方根值/mGal4.323.224.213.014.473.07比对点数5.335 11×1055.335 11×1054.249 622×1064.249 622×1063.431 619×1063.431 619×106从表1比对结果可以看出：V28模型在3个海区与实测数据的比对互差均方根值均在3 mGal左右，如果考虑船测数据可能有1～2 mGal的观测误差，那么可以认为，V28模型的外符合精度已经达到2～3 mGal的水平，与国外学者和研究机构的同类比对结果基本一致[16-17]．相比较而言，EGM2008模型的比对精度要略低一些，互差均方根值均超过4 mGal，相对应的外符合精度应当在3～4 mGal水平．表1同时显示，除了EGM2008模型在区块3的比对结果存在1.6 mGal大小的系统差外，其他比对结果都不存在明显的系统性偏差．上述比对检核结果再次说明：最新发布的国际卫星测高重力模型在可靠性和精准度两方面都得到了很大提升，因此具有较高的参考使用价值．2 海空重力测量误差检测中的应用如前所述，近期发布的卫星测高重力模型在精准度和可靠性方面都已经达到了较高的水平，在全球海底地形反演、大尺度重力场逼近计算和空间基准建立等领域得到了广泛应用[16-17]．尽管与高精度高分辨率的海面船载重力测量相比，卫星测高重力模型在表示海洋重力场细部特征方面还有一定差距，但作为海洋重力场中长波趋势性变化的基准是值得信赖的[18-20]．据此，研究探讨了将卫星测高重力模型应用于海空重力测量误差检测和校正的可行性及有效性，主要聚焦粗差检测和系统偏差补偿两个方面．2.1　在粗差检测中的应用根据定义，粗差是指在一定的测量条件下，明显超出统计规律预期值的异常误差[21]．海空重力测量中的粗差是指在作业过程中，偶尔出现测量环境条件的剧烈变化迫使测量载体突然改变航向和航速，或测量仪器偶尔受到突然冲击、强烈震动而出现影响工作稳定性的故障，或作业人员操作设备时偶尔出现差错等，使得相应的观测值超出了正常的预期值，这样的观测值又称为异常值．显然，受粗差污染的测量成果不仅会扭曲海空重力场的变化形态，还会影响测量成果质量评估的可靠性，因此必须预先将粗差从观测数据系列中剔除．测量中的粗差一般有大粗差和小粗差之分，大粗差指量值远大于由正常随机干扰因素引起的误差，小粗差指量值比较接近于正常随机干扰因素引起的误差．可以采用不同的处理策略来检测和剔除两类不同性质的粗差，对于大粗差，可分别采用阈值排除法和基准比对法进行处理；对于小粗差，可采用现代稳健估计方法进行定位．阈值排除法的依据是：理论上，全球重力异常变化幅度最大不会超过正负几百个mGal，故海空重力测量获取的重力异常值不可能超过某个阈值．就某个测量区域而言，可以借助前面介绍的卫星测高重力模型预先设定一个有一定余量的阈值．考虑到卫星测高重力模型已经具有较高的可靠性，对于重力场变化特征比较剧烈的区域，可将排除观测粗差的阈值(Ymax/min)设定为Ymax/min=[Δgmin-50,Δgmax+50]，(1)式中Δgmin和Δgmax分别表示由卫星测高重力模型计算得到的重力异常数据序列中的最小值和最大值．对于重力场变化特征相对平缓的区域，可将排除观测粗差的阈值设定为Ymax/min=[Δgmin-30,gmax+30]．(2)当海空重力测量的某个重力异常观测量Δg超出式(1)或式(2)设定的阈值时，有理由认为该观测值为大粗差．基准比对法的理论依据来自数理统计原理，具体做法与阈值排除法有一定的相似性，同样须要首先确定一个阈值，然后判断比对对象是否超出预设阈值．这里的比对对象选择为海空重力测量观测值与测高重力模型值的互差，也就是以测高重力模型值作为比对基准，比对阈值(ΔYmax/min)可依据海空重力测量精度估值(mc)和测高重力模型精度估值(mw)进行确定，具体为ΔYmax/min=[-3(3mc2+mw2),3(3mc2+mw2)], (3)式中：右端第一个3倍系数的意义在于确认可能出现小粗差的置信限；第二个3倍系数则是为了确认可能出现大粗差的置信限．式(3)的意义等价于数理统计理论中的极差限定条件，它规定了一个正常样本序列值的极端离散范围[22]．当海空重力测量精度估值和测高重力模型精度估值分别取为：mc=2 mGal，mw=4 mGal时，可求得ΔYmax/min=[-40，40] mGal．这样的量值与前面阈值排除法确定的[Δgmin-50,Δgmax+50]和[Δgmin-30,Δgmax+30]基本一致．当海空重力测量某个观测值与测高重力模型计算值的互差超出式(3)设定的阈值时，也有理由认为该观测值为大粗差．为了提高粗差检测的准确率，实践中可结合前面的数值计算和三维立体图形或测线剖面曲线变化形态显示情况及误差形成机理分析，对可疑观测点做出综合判断．海空重力测量大粗差大多出现在测线起始和结束、航向和航速突变及仪器掉格线段，数据处理中应给予重点关注．在前面基于基准比对法检测大粗差基础上，可进一步采用现代稳健估计方法对剩余的小粗差进行定位．具体原理是：以海空重力测量观测值与测高重力模型计算值的互差作为新的观测值系列，将每个观测点作为被检测点，利用每个检测点周围的新观测值(不含检测点处的新观测值)，通过抗差M估计的选权迭代方法，对每个检测点进行推估计算，然后将估算值与检测点上的新观测值作比较，最后根据比较结果大小判断该检测点是否为粗差点．选用抗差估计方法做推估计算，是为了减弱局域内出现多个粗差时引起的交叉干扰．这里选用相对简单的加权平均数模型作为检测点的推估模型，计算公式为[23]Lpk+1=∑i=1np¯ikLi∑i=1np¯ik，(4)式中：Li为被检测点周围邻域内参加计算的新观测值；n为参加计算的新观测值个数；Lpk+1为检测点处的第k+1次推估值；k为迭代次数；p¯i为第i个新观测值对应的等价权因子．其中p¯i一般采用中国学者提出的IGGⅢ方案进行计算，具体为[24]p¯i=pi                      (vi'≤k0);pik0[(k1-vi')/(k1-k0)]2vi'        (k0vi'≤k1);0                  (vi'k1), (5)式中：vi=Lpk-Li；vi'=vi/σp，σp=[(pipivi)/(n-1)]1/2为由残差vi计算得到的单位权中误差，pi=1/(di+ε)2为第i个新观测值对应的原始权因子，依据观测点至被检测点的距离计算，di为第i个观测点至被检测点的距离，ε为一小正数，其作用是防止权函数的分母趋于零；系数k0和k1可分别取为1和3．为了提高迭代计算式(4)的稳定性，可选用Lp0=medi{Li}确定迭代计算的初值，同时改用σpk=medi{vik}/0.674 5计算单位权中误差[22]，med{⋅}表示取该数列的中位数．按照式(4)进行迭代计算得到稳定的推估结果(L̂P)后，可将其与被检测点上的新观测值(LP)做比较，求两者的互差值ΔLP=L̂P-LP．根据ΔLP的绝对值大小和预先设定的限差标准，即可对检测点观测值是否包含小粗差做出判断，这里推荐采用的判别计算式为ΔLPk2σP2+mc2+mw2．(6)当取σp=3 mGal，mc=2 mGal，mw=4 mGal，k2=2时，式(6)右端等于11 mGal．应用实践表明，对于小粗差，设定这样的量化判别标准较为合适．2.2　在系统误差检测中的应用根据定义，广义的系统误差是指：对同一量的多次观测过程中，保持恒定或以可预知方式变化的观测误差分量[21]．测量领域的系统误差一般指大小和符号都表现为系统性和规律性变化的观测误差．系统误差有定值和变值系统误差之分，变值系统误差又有线性变化、周期性变化和复杂规律变化系统误差之分．海上作业环境复杂多变，干扰海空重力测量的外部因素多种多样，因此海空重力测量数据存在系统性偏差也属正常现象．小的系统偏差可以忽略不计，当系统偏差达到一定量值时，比如超过仪器正常观测精度的1/4时，必须对其进行分析判定并做适当补偿．从误差形成机理角度分析，海空重力测量全过程都可能出现上述各类性质的系统误差，仪器参数标定比如格值测定误差、重力基点传递误差和仪器出现突然掉格等都可引起测量成果产生恒定(定值)的系统偏差，仪器出现缓慢掉格可引起观测数据产生线性变化的系统偏差，海上规律和非规律的风、流、浪、涌干扰可引起海空重力测线数据产生周期性和非周期性变化系统偏差．海空重力测量系统误差一般可通过建模方法加以判别和修正，比如可利用东西正反向重复测线数据检测和修正海空重力仪格值测定误差[25]，可通过始点和终点闭合重复测量方法，对仪器缓慢掉格引起的线性化误差进行改正[3]，可分别采用测线网平差和相关分析方法对各类变值系统误差进行补偿[6-7,26-28]．对于表现为测区或测线维度上的海空重力测量定值系统误差，可考虑采用前面介绍过的卫星测高重力模型作为比对基准，对此类误差进行检测和补偿．因测高重力模型对应的是全球坐标和重力基准系统，在表示中大尺度重力场变化特征方面具有较高的可靠性，故特别适合作为大区域大跨度海空重力测量的质量控制基准．在完成观测粗差剔除以后，可以整个测区或单一测线(最好跨度在100 km以上)为单元，计算每个测点的观测值与测高重力模型值的互差，进而计算单元内所有互差的平均值，当这样的平均值超过海空重力测量正常观测精度的1/4时，比如超过±0.5 mGal时，可认为单元内的观测量存在比较明显的定值系统偏差，应对其产生的可能原因进行机理分析和技术确认，并做出相应的补偿，比较简单的补偿方法是从计算单元的每个观测值中减去对应的互差平均值．上述方法可简称为单元平均值判别法．理论分析和应用实践都表明，计算单元划分也就是测量数据覆盖区域越大，测高重力模型作为比对基准的可靠性就越高，在此基础上完成定值系统偏差补偿的有效性也就越高．但计算单元划分过大，可能不利于局部系统偏差或半系统误差(比如测线系统差)的发现和补偿．如何在两者之间取得较好的平衡，须要在作业实践中结合数据应用需求，不断总结和积累误差分析处理的相关经验并加以改进．3 应用实例分析为了验证前述检测粗差和系统误差方法的有效性，这里特别选用两组比较典型的海上测量航次重力观测数据作为试验样例，分别对大粗差、小粗差及系统误差的检测效果进行计算和分析．其中，第1组为一个独立航次的海面船载观测网数据，包含161条有效测线，交叉点870个，测点总数接近1.8×105．因该测区位于岛礁密集区，测量过程中改变航向航速的频率较高，加上使用的测量载体吨位相对较小、受海面风浪影响较大，故该航次获取的重力观测数据质量相对较差，原始观测数据交叉点重力不符值均方根值为3.82 mGal；第2组为涵盖一个大测区、包含来自两个部门于不同时间段执行的5个航次海面船测数据，测点总数超过5.55×106，尽管各个航次的观测海况条件都比较正常，但不能排除这些航次的测量数据中仍存在少量的粗差和一定大小的系统偏差．首先采用测高重力模型V28作为比对基准和前面介绍的粗差检测方法，对第1组观测数据依次进行大粗差和小粗差检测，发现大粗差点423个、小粗差点2 710个，其中大粗差采用基准比对检测法．剔除大粗差点后，重力观测数据均方根值从原来的9.85 mGal减小到9.30 mGal，进一步剔除小粗差点后，其均方根值减小为9.15 mGal．在此基础上，进一步开展系统误差检测，分两种检测单元划分方式进行：一种以该航次覆盖的整个测区作为检测单元，另一种以每一条测线作为检测单元．统一对系统偏差超过±0.5 mGal的单元数据作补偿处理，具体检测结果为：测区单元检测发现存在-0.82 mGal的系统偏差；测线单元检测发现有104条测线存在系统偏差，最大偏差值为1.87 mGal，最小偏差值为-3.24 mGal．从以上检测结果可以看出：受特殊海况条件的影响，该航次原始观测数据质量远低于测量规范规定的2 mGal精度要求，不仅存在超过3 000个测点的粗差(占测点总数的1.76%)，还存在比较明显的测区和测线系统性偏差，存在测线系统偏差的测线数占有效测线总数的12%，偏差最大值超过3 mGal．经采用前述方法对大小粗差进行剔除和系统偏差补偿后，该航次的测线交叉点重力不符值均方根差值从原来的3.82 mGal减小到1.18 mGal，充分说明了本文处理方法的合理性和有效性，也从另一个侧面说明测高重力模型V28具有较高的可靠性．这里须要补充说明的是，测区系统差对所有测点都是常值，因此其大小并不改变交叉点不符值的统计结果；测线系统差在单一测线内部是常值，但在测线之间不存在相关性，因此其大小必然在交叉点重力不符值中得到一定程度的体现，对测线系统差作补偿后，交叉点重力不符值均方根差值显著减小，正好反映了卫星测高重力模型作为系统差检测和补偿基准的有效性．从检测到的粗差分布图(这里从略)也可以看出：粗差测点大多都出现在测线的头部和尾部，以及航向航速突变较大的线段，这样的结果与先前所作的理论分析相符．采用同样的处理步骤对第2组观测数据进行粗差和系统误差检测，表2列出了各航次的粗差检测统计结果，其中的互差均方根是指测点观测值与模型计算值的互差均方根值；表3给出了分别以航次和测线为检测单元且系统偏差超过±0.5 mGal的对应航次及作系统误差补偿处理前后的交叉点不符值统计结果．10.13245/j.hust.220918.T002表2第2组观测数据粗差检测统计结果航次测点数/106粗差个数最大粗差/mGal最小粗差/mGal剔除前互差均方根/mGal剔除后互差均方根/mGalQ011.143 5225 02444.52-27.7519.5419.49Q021.314 1061 30034.25-11.8926.0226.02Q031.161 1059312.90-16.909.379.37Q041.126 9921 33838.61-39.2110.3910.39Q050.805 4354731.4110.0211.3311.3310.13245/j.hust.220918.T003表3第2组观测数据系统误差检测统计结果航次原始数据统计测区单元检测测线单元检测测线总数交叉点个数交叉点不符值均方根/mGal系统偏差/mGal存在系统差测线数最大偏差/mGal最小偏差/mGal补偿后交叉点不符值均方根/mGalQ011981 0061.11-0.341233.00-2.640.35Q021271 5630.55-0.34360.55-1.410.15Q032141 1280.76-0.171061.96-2.560.22Q042299880.91-1.141685.21-3.590.21Q051139311.050.60763.05-3.240.16从表2和表3的统计结果可以看出：即使是一些常态化的海上测量航次，也难免出现少量由各种异常现象引起的粗差数据和由规律性变化干扰因素引起的系统偏差．算例中，虽然各航次的测区系统偏差较小，但都出现了比较明显的测线系统偏差，这是海空重力测量采用测线作业方式带来的固有缺陷，也是评估海空重力测量数据质量须要重点关注的问题之一．经粗差剔除和系统偏差补偿后，各航次的测线交叉点重力不符值均方根差均大幅减小，下降幅度都在3倍以上，再次验证了本文处理方法的可行性和有效性．4 结语误差分析处理是海空重力测量数据精细化处理的一个重要环节，粗差和系统误差检测是提升海空重力测量数据质量的决定性因素，这一步骤的关键点是确立一个稳定可靠的数据比对基准。为此，开展了卫星测高重力模型精度分析与评估研究，使用多个区域的海上实测重力数据做比对，结果表明：国际上最新发布的卫星测高重力模型与实测数据的符合度在绝大多数海域已经达到±3 mGal的水平，基本满足作为海空重力测量数据质量控制基准的精度要求．在此基础上提出了基于卫星测高重力模型基准检测大粗差的阈值排除法和基准比对法，检测小粗差的抗差选权迭代互比法，以及分别以测区和测线为计算单元，检测补偿系统误差的平均值判别法．多个航次的数据处理应用结果表明：本研究推荐的上述粗差和系统误差检测方法是可行和有效的，具有较高的实用价值．