坐标计算相对误差（武汉大学测绘学院赵闯）

坐标计算相对误差（武汉大学测绘学院赵闯）4. 广州海洋地质调查局资料处理研究所 广东 广州 5100003. 地球空间环境与大地测量教育部重点实验室 湖北 武汉 430079; 秦朋波4 杨连俊11. 武汉大学测绘学院 湖北 武汉 430079;2. 湖北珞珈实验室 湖北 武汉 430079;

本文内容来源于《测绘学报》2023年第

顾及残差约束的多面函数法融合卫星测高和船载重力

赵闯¹

金涛勇^{1 2 3}

秦朋波⁴ 杨连俊¹

1. 武汉大学测绘学院 湖北 武汉 430079;

2. 湖北珞珈实验室 湖北 武汉 430079;

3. 地球空间环境与大地测量教育部重点实验室 湖北 武汉 430079;

4. 广州海洋地质调查局资料处理研究所 广东 广州 510000

基金项目：国家自然科学基金(41974020;42192533;41721003);湖北珞珈实验室专项基金

摘要：融合多源重力数据是精化局部区域重力场的必要手段。本文通过引入残差约束因子 改进了基于多面函数的多源海洋重力数据解析融合方法。以日本群岛附近海域为例 使用卫星测高反演重力场模型和船载重力数据 分析了分块大小、融合船载重力数据量及其分布对该方法的影响 并与改进前多面函数法及最小二乘配置法的结果进行比较 验证了改进方法的可靠性和有效性。与验证船载重力数据相比 顾及残差约束的多面函数法精度最好 极值和标准差均显著下降。与多面函数法相比 不仅较好地改善了控制船载测点附近的不符值 而且将其较好地推估至整个区域的船载观测数据之间 残差分布合理 且无明显测线状误差。

关键词：卫星测高 船载重力 残差约束 多面函数拟合 最小二乘配置

引文格式：赵闯 金涛勇 秦朋波 等. 顾及残差约束的多面函数法融合卫星测高和船载重力[J]. 测绘学报，2023，52(4)：605-613. DOI: 10.11947/j.AGCS.2023.20210444

ZHAO Chuang JIN Taoyong QIN Pengbo et al. An improved multi-surface function method with residual constraint for the fusion of shipborne and satellite altimetry derived gravity data[J]. Acta Geodaetica et Cartographica Sinica 2023 52(4): 605-613. DOI: 10.11947/j.AGCS.2023.20210444

阅读全文：http://xb.chinasmp.com/article/2023/1001-1595/20230409.htm

引 言

地球重力场结构反映了地球物质质量分布特征，是地球内部结构探测、矿产资源勘探、地球物理和地球动力学研究等的重要基础^[1-3]。约71%的地球表面被海洋覆盖，海洋重力场信息是获取地球重力场精细结构的重要组成部分。当前，获取海洋重力场信息的手段包括船载重力测量、航空重力测量、卫星重力测量和卫星测高反演等。其中，船载重力测量是获得高精度海洋重力场信息的主要方式，在沿测线方向上可得到高频率采样重力数据，但由于测线分布稀疏、不均匀，会存在大量空白，导致整体空间分辨率较低，且测量周期长、成本高，单纯依靠船载测量得到高精度高分辨率的全球海洋重力场难度较大^[4-5]。卫星测高技术的发展很大程度上改变了海洋重力测量的模式，解决了广阔海洋区域的重力数据空缺问题^[6-7]。但卫星测高反演得到的海洋重力场高频信息不如船载数据，在近岸浅水等区域，因受到陆地的干扰，数据精度还有所下降^[8]。

针对卫星测高数据空间覆盖较高但在近岸、极地等部分海域精度较差，同时，船载重力数据整体空间分辨率低、测线间距大，存在很多空白区，但在近岸浅水等区域精度较高的情况，综合利用船载重力测量和卫星测高技术，进行多源重力数据融合，是精化局部海洋重力场的必要途径^[9-12]。多源重力数据融合的方法主要有统计法和解析法两大类型。最小二乘配置法是统计法的典型代表，但由于必须以足够分辨率的观测数据来构建经验协方差函数，并且协方差矩阵求逆还存在不稳定性问题，在一定程度上限制了其融合效果^[13-14]。点质量法^[15]因计算模型简单、使用方便成为解析法的主要代表，但求解过程涉及向下延拓，并且矩阵求逆常受到边界面和边界值粗糙化或奇异性的影响，有一定不稳定性^[16-17]。

如果船载与卫星测高重力异常均已归算到平均海平面上，两种数据的融合可以看作网格化重力异常数字模型构建问题，即数据的拟合插值过程^[18]。但对船载与卫星测高重力异常进行网格化拟合插值(如最小二乘配置法)的过程中，往往会在船测线中间地带出现数据空缺的情况，导致不能有效修正卫星测高重力异常，使得融合后局部重力异常的精度不均匀。针对这一问题，文献[19]提出了一步融合和分步融合的解析方法，利用融合多面函数的趋势面函数对多源重力数据的系统偏差进行平差改正，再考虑数据的空间相关性和精度水平差异，采用双权因子进行调配，有效修正了船载重力测线间的阶梯状变化。在此基础上，本文在多面函数模型中进一步引入残差约束因子，对其进行改进，并以日本群岛附近海域为例，利用船载和卫星测高反演的重力数据进行融合试算和验证。

1 数据与方法

为对多源重力数据融合方法进行比较和验证，本文考虑高精度船载重力数据收集的可行性，选取日本群岛附近海域(137°E—142°E，30°N—35°N)为研究区域，收集卫星测高反演海洋重力场和船载重力观测数据，设计对比试验，分析了方法的可行性和可靠性。

1.1 海洋重力场模型

卫星测高反演海洋重力场模型采用了丹麦科技大学研发的DTU13海洋重力场模型。该模型联合Geosat、ERS-1、T/P、GFO、ERS-2、Jason-1、Envisat和Cryosat-2测高数据，并使用ICESat激光测高数据和ArcGP重力测量数据填补极区空白，通过波形重跟踪、地球物理和环境误差改正等系列数据处理后，采用DOT07A模型考虑稳态海面地形的影响，最后使用最小二乘配置法获取规则网格化的大地水准面高，利用逆Stokes公式、FFT快速解算、Wiener滤波和EGM2008参考模型恢复得出全球1′×1′的海洋重力异常，与部分区域船载重力数据相比，精度可达3.15 mGal^[20-22]。

1.2 船载重力数据及处理

船载重力测量数据采用日本海洋地球科学技术厅发布的自由空气重力异常数据(http://www.godac.jamstec.go.jp/darwin/e)。研究区域内卫星测高重力异常模型及船载重力数据分布如图 1(a)所示，包括371条船载重力测线，总计3 304 907个数据点，测线在整个区域均有分布，部分海域较为密集。卫星测高重力异常在反演中采用移去-恢复技术，并经过滤波处理，整体较为平滑，在近岸浅水区域精度可能较低。船载重力测量是在运动状态和复杂环境下进行，易受到海浪、航速、海风等扰动因素影响^{[2 23]}，并且仪器误差和导航误差也不能忽视^[24]，数据一般会存在粗差，应进行海面大气压改正、漂移校正、厄特弗斯效应校正及粗差剔除等一系列预处理^[25-26]。本文使用的数据已进行相应改正，并已通过质量控制删除了以下3种情况下的观测数据：自由空气异常瞬时变化超过10 mGal/km、厄特弗斯修正值变化超过3 mGal/min、船速低于3 km^[27]。

图 1 研究区域内测高重力异常模型及船载数据分布Fig. 1 The satellite altimetry derived marine gravity field model and the distribution of shipborne gravity data in the study area

图选项

为进一步增强两种重力数据的一致性，对船载重力数据按如下流程进行了筛选和编辑：①将DTU13模型内插至船载测点处，计算模型值与船载测量值之间的残差；②在船载测量值中减去每条测线残差的平均值以改正两种观测方式之间的系统偏差；③以卫星测高重力异常数值作为参考，利用窗口移动粗差探测法^[28]对各测线进行粗差剔除，最终剩余3 254 212个数据点，约占原始总数据的98.5%。随后，统计各测线实测数据与卫星测高反演数据的残差中误差，结果显示每条测线的残差中误差均小于7.32 mGal，与其他学者对卫星测高与船载重力数据比较所得的差异基本一致^{[23 29-30]}。剔除的数据分布情况如图 1(b)所示，可以看出这些数据点大多分布在近岸半开阔海域或地形起伏较大的区域。

1.3 顾及残差约束的多面函数拟合法

若卫星测高反演海洋重力和船载测量重力都已归算到海平面，则将两种重力数据的融合看作是数据拟合内插过程。将卫星测高重力异常内插到船载重力数据每个测点上，并将其看作公共重复测点。由于数据融合的目的在于用较高质量的船载数据校正卫星测高反演的重力数据，故假设船载测量重力为真值，在公共重复测点处建立以下数学模型^[19]

(1)

式中，Δg_sa为卫星测高重力异常；Δg_sh为船载重力异常；F(φ λ)为卫星测高重力异常与船载重力异常之间的误差趋势面函数；Δ为观测白噪声。

对于误差趋势面函数F(φ λ)，文献[19]设计了以下形式

(2)

(3)

式中，φ、λ为测点的纬度和经度；采用(φ_min~φ_max λ_min~λ_max)表示两种数据重叠区域大小，Δφ=φ－φ_min，Δλ=λ－λ_min；a_i、b_i为待定系数；Q_i(φ λ)为多面函数，(φ_i λ_i)为数据节点坐标，一般选测点局部特征点(即局部最大最小值点)作为数据节点，为了使选取的所有数据节点分布尽可能均匀，通常将研究区域均匀分为若干小区域，在小区域内分别选取局部特征点；c为拓展因子，一般选取与测点最小间距大小相当。

设式(1)中Δg_sa的改正数为V，将式(1)写成矩阵形式，则可在卫星测高和船载测量重力异常的公共重复测点上，建立误差方程，并利用最小二乘进行求解

(4)

式中，L代表卫星测高重力异常和船载重力异常的不符值向量；X为式(2)中误差趋势面函数的待求参数向量；A为其系数矩阵；P为L的权矩阵，计算时可视为等精度观测。

文献[19]以上述趋势面函数为基础，对卫星测高和船载重力数据间的残差进行改正后，再顾及数据的空间相关性和精度水平差异进行综合定权，实现了两种重力数据的融合，有效改善了船载测线间卫星测高重力数据的修正效果，避免了局部区域融合重力异常的阶梯式变化。但所采用的多面函数仅考虑了经纬度位置及测点间距，在此基础上，本文保持误差趋势面函数的形式不变，在式(3)的基础上，引入残差约束因子，得到

(5)

(6)

式中，d为残差约束因子；ω为测点处卫星测高重力异常和船载重力异常观测值的残差；ω_i为数据节点处卫星测高重力异常和船载重力异常观测值的残差。该约束因子的主要作用为：在数据节点残差与测点残差相差过大时，可以抑制误差的传递；当数据节点残差与测点残差相差不大时，能够较好地保留两种重力异常的残差特征。通过上述过程，使用船载重力数据对卫星测高重力异常进行修正，将得到的修正值恢复到卫星测高重力异常上，即得到计算区域内最终的格网重力异常模型。

具体解算步骤如下：

(1) 利用卫星测高重力异常格网模型，内插得到船载测点上的模型重力异常，并计算模型重力异常与船载测量重力异常之间的残差。

(2) 选取试验区域内一定比例船载测点的重力异常测量值与卫星测高重力异常的残差作为融合数据，剩余船载测点用于精度验证。

(3) 确定特征点，为使得选取的所有数据节点在研究区域内较均匀分布，在1°×1°范围内均匀划分为N×N的区块，在每个区块内检索两种重力异常数据残差的最大值和最小值点作为局部特征点。

(4) 利用研究区域内所有数据节点和参与融合的船测点重力异常残差数据计算系数矩阵A，进一步求得参数向量X，并利用参数向量X计算所有点的残差改正。

(5) 将计算得到的残差改正恢复到卫星测高重力异常格网，得到融合后格网重力异常模型，再内插到用于检核的船载测点上，与船载测量重力数据比对，进行融合效果的精度评估。

2 结果与分析

根据趋势面函数和多面函数的形式可知，在计算过程中数据节点的选取对融合效果会产生影响，而分块数量N、参与融合的船载重力数据量及其空间分布情况会影响最终所有数据节点的数量和分布。因此，下面分别分析分块数量、融合船载重力数据的比例和空间分布对融合结果的影响，并与最小二乘配置法及改进前多面函数拟合法的融合结果进行对比。

2.1 分块数量选择

为评估改进后的多面函数解析法在改善卫星测高重力异常精度的效果，利用移去-恢复技术，随机选取试验区域内70%的船载测点的重力异常测量值与卫星测高重力异常的残差作为融合数据，剩余30%的船载测点用于精度验证，在研究区域内选择不同分块数量N，比较验证点处船载重力与DTU13重力异常模型及改进多面函数法融合后结果的差值。此外，为了得到相对稳定的结果，不同分块数量的情况下，保证融合数据和验证数据的比例不变，均进行了10次随机数据选取，并取其平均作为计算结果(表 1)。

表 1 不同分块数量下改进多面函数法融合结果与船载重力差值统计Tab. 1 The statistics of the difference between fusion results of improved multi-surface function method and shipborne validation data using different block numbers

mGal

表选项

由表 1可以看出，随着分块数量N从1增加到4，特征点数量增加，融合结果精度逐渐提高，而当增加到5时，结果明显变差。同时，在10次随机选取的计算中发现：对于不同融合数据，分块数量不超过4时结果较为稳定，并且每次结果相差不大；而当分块数量为5时，10次计算中有9次均出现较大的极值(表 2)。这可能是因为，当分块数量过大时，选取的数据节点过多，导致多面函数解析模型的系数矩阵变得病态，求逆过程出现奇异，进一步造成改正结果不合理。因此，针对不同区域或不同分布特征的数据，最优分块数量随着船载重力数据分布情况的改变会有一定差异。此外，随着分块数量的增加，计算效率迅速降低。综合考虑模型的适定型、计算效率和精度，针对该区域的融合数据，本文随后均选取N为4进行计算。并且由于分块数量为4时结果稳定，后续均以一次计算为准，不采用10次计算取平均的方式。

表 2 分块数量为5时改进多面函数法的10次随机计算结果与船载重力差值统计Tab. 2 The statistics of the difference between 10 random calculations of the improved multi-surface function method with 5 blocks and the shipborne validation data

mGal

表选项

2.2 融合船载重力数据分布的影响

在确定试验区域多面函数方法分块数量后，因参与融合的船载重力数据分布对融合效果有一定影响，而船载数据一般以测线为单位，因此，在保证参与船载重力数据比例为70%的前提下，分别采用以船载测点为单位和船载测线为单位的随机取样方式，仍以剩余30%船载重力数据作为验证点，分析了船载重力数据分布对融合结果的影响。以船载测点为单位和测线为单位的70%随机船载数据分布分别如图 2(a)、(b)所示。

图 2 测点取样和测线取样选取的70%随机船载重力数据分布Fig. 2 The distribution of 70% random shipborne gravity data selected by observation points and tracks

图选项

利用两种取样方式得到船载重力数据，采用改进多面函数法融合后结果与剩余30%船载重力数据比较，其差值统计见表 3。可以看出，融合数据的不同分布情况对结果精度有较大影响。基于船载测点的随机数据比基于船载测线的随机数据分布较为均匀，故融合后结果标准差更小，精度更优。主要原因在于，以船载测线为单位的随机取样，会不可避免地造成一些测线稀疏区域的融合数据空白，进一步导致局部特征点的选取缺失问题。因此，当融合数据分布更均匀时，融合结果的精度更好。

表 3 70%比例下不同分布船载重力数据融合结果与剩余船载重力的差值统计Tab. 3 The statistics of the difference between the fusion results and the validation data under different distributions of shipborne gravity data at 70% ratio

mGal

表选项

2.3 融合船载重力数据比例的影响

基于上述两个原则，选择分块数量N为4，以船载重力测点为单位的取样方式选择参与融合的重力数据，下面分析参与融合的船载重力数据比例对融合结果的影响。分别选取随机取样比例为90%、70%、50%、30% 4种情况，剩余的10%、30%、50%、70%作为其对应的验证数据。不同比例下融合后结果与验证船载重力的差值统计见表 4。可以看出，选用不同比例的数据进行融合，融合结果的验证精度变化不大，可能原因是，在总数据量较大并且不分测线随机选取融合数据的情况下，当选取的数据达到一定比例，如50%时，已能保证有分布较为均匀的船载重力数据，足够对整个区域进行较好的数据质量控制，因此融合后结果精度较好，相差较小。而当参与融合的数据比例较小，如30%时，作为控制数据的分布并不理想，使得融合后结果精度下降较大。

表 4 不同比例船载重力数据融合结果与剩余船载重力数据的差值统计Tab. 4 The statistics of the difference between the fusion results and validation data under different proportions of shipborne gravity data

mGal

表选项

2.4 与其他方法的精度比较

基于上述分析，选定分块数量N为4，以船载测点为单位随机选取70%的数据参与融合，剩余30%船载数据作为验证，比较了改进多面函数法与改进前多面函数法及最小二乘配置法的融合效果(表 5)，同时也给出了DTU13卫星测高重力场模型与验证船载重力数据的比较结果。其中，最小二乘配置法选择二阶Markov协方差函数，相关距离选择70 km。

表 5 不同方法融合重力场及DTU13模型与验证船载重力数据的差值统计

Tab. 5 The statistics of the difference between the fusion results by different methods DTU13 model and the shipborne validation data

mGal

表选项

由表 5可以看出，与验证船载重力数据相比，利用3种方法，融合船载重力数据后结果较未融合的DTU13测高重力场模型的精度均有提升，反映了船载重力数据对精化区域重力异常场的贡献。其中，从平均值和标准差来看，最小二乘配置法和多面函数法有相当的效果，融合后精度提升均达到约60%。但从最大值和最小值来看，最小二乘配置法出现了较原模型更大的极值，表明融合后结果精度不均匀，在部分区域引入了更大的极值误差，而多面函数法极值与原模型相当，仅将两种重力数据的不符值进行了分配。总体而言，改进后的多面函数法在所有统计量上均表现最优，平均值和标准差均接近于零，表明将卫星测高与船载重力数据之间的不符值进行了较为均匀合理的分配，避免了不符值极值误差对其他区域融合结果的影响。但需要说明的是，基于该统计结果，可以看出，改进多面函数法将两种重力数据的残差进行了全区域较为均匀分配，这就对参与融合的船载重力数据质量提出了较高要求，如船载重力数据中存在粗差，将不可避免地对融合后结果产生影响。

为了直观体现不同方法融合效果，图 3给出了最小二乘配置法、多面函数法和改进多面函数法相对于DTU13重力场模型的残差图，以及改进多面函数法与多面函数法融合结果的残差图。可以看出，最小二乘配置法融合结果有很明显的测线状残差存在，仅在船载测线附近对原测高重力进行了改善，测线间改善极小，进而导致在部分复杂区域出现较大的极值；多面函数法融合结果比最小二乘配置法融合结果有一定改善，融合前后残差分布较为平滑，但在船载测线附近仍出现了较明显的误差；改进后多面函数法融合结果没有明显的船载测线误差影响，不仅在船载测点附近改正了卫星测高重力数据，还将公共测点不符值推估至船载测量线间的区域，相较于原模型，残差分布较为合理地体现了重力异常的显著变化，如海底地形的突变边界等，表明引入残差约束因子后，改善了残差的分配状况，提升了整体融合精度水平。

图 3 不同方法融合海洋重力场相对DTU13模型及多面函数法改进前后融合结果的残差

Fig. 3 The differences between the fused marine gravity by different methods and DTU13 model and by the multi-surface function method improved before and after

图选项

3 结论

卫星测高反演海洋重力场与航空重力、船载重力等多源实测数据的融合是提升海洋重力场精度和分辨率的重要途径。本文利用基于多面函数的多源海洋重力数据解析融合方法，引入残差约束因子，更为合理地分配不同重力数据之间的不符值，并选择日本群岛附近海域，基于DTU13卫星测高重力异常模型和船载重力观测数据，考虑到研究区域内数据节点数量与分布对多面函数法的影响，分析了分块数量、参与融合的船载重力数据量及其分布对融合结果的影响，并与改进前多面函数法和最小二乘配置法进行了比较，验证了改进方法的可靠性和有效性。结果表明，引入残差约束因子后，与验证船载重力数据相比，改正前后极值和标准差改善效果提升显著，不仅较好地改善了控制船载测线附近的不符值，而且将其较好地推估至整个区域的船测测线之间，残差分布合理，且无明显测线状阶梯误差，提高了整体融合精度，可利用有限船载重力数据提升全区域海洋重力场精度。

随着海洋重力测量仪器和技术的不断发展，船载重力数据的质量越来越好。在保证其精度的前提下，应将所有的船载重力数据融合到测高重力异常模型中，但每个区域的实际观测数据量、分布和精度有较大差异，本文针对参与融合的船载重力数据量及其分布的研究可为其他区域数据融合提供一定参考。然而，本文方法也存在一定局限性，其核心在于数据节点的数量、分布及其残差大小，使得针对不同船载重力数据情况，最优分块数量存在差异，且需要首先剔除船载重力数据粗差，以保证数据节点的残差不受影响。

作者简介

第一作者简介：赵闯(1997—) 男 硕士生 主要从事多源重力数据融合算法研究。E-mail: 2015301140010@whu.edu.cn

通信作者：金涛勇 E-mail：tyjin@sgg.whu.edu.cn

初审：张艳玲

复审：宋启凡

终审：金 君

资讯