论文专区▏海空重力测量及应用技术研究进展与展望（四）：数值模型构建与综合应用技术
一、引言
在文献[1]中，我们简要介绍了海空重力测量及应用技术的研究背景，详细论述了海空重力测量信息的应用价值及开展该项技术研究的目的意义和应用需求，分析了海空重力测量的技术特点，提出了海空重力测量及应用技术体系的基本架构，明确了该项技术的研究内容、作业流程及各个技术环节的相互关系。在文献[2]中，我们分三个阶段系统分析总结了海空重力仪的发展进程及研究现状，简要介绍了作为海空重力测量重要组成部分的定位系统的技术发展状况，分析展望了海空重力仪稳定性测试与评估技术的发展前景，简要论述了海洋重力场特征的分析方法及研究方向，分析讨论了海洋重力场信息应用需求与海空重力测量规划设计的内在联系及海洋重力测线布设技术的研究进展。在文献[3]中，我们分析总结了海空重力测量环境效应改正技术的研究现状及发展方向，简要论述了海空重力测量数据滤波与误差补偿技术的研究动态及发展前景，分析讨论了海空重力测量精度评估技术的研究进展及发展思路。本文将着重对海空重力测量技术体系中的数值模型构建和数据综合应用两个后端环节进行分析和总结，其中，数值模型构建环节涵盖重力向上延拓、向下延拓和多源数据融合处理三个部分，数据综合应用环节包含地球外部重力场赋值和大地水准面精化两个部分。
二、海空重力测量数值模型构建技术
数值模型构建是海空重力测量成果向应用领域转换，即从测量成果转化为数字产品必不可少的关键环节。由于海洋和航空重力测量成果涉及不同的数据归算基准面，不同数据源之间具有比较明显的异构性特征，故这个技术环节的研究内容不仅包括多源重力数据的融合处理，还包括前端的重力向上延拓和向下延拓两项关键技术。
⒈ 重力向上延拓技术
重力向上和向下延拓是地球外部重力场赋值和多源数据融合处理最常用的技术手段之一[4,5]。如前所述，最近一个时期，在矿产资源开发和航天技术保障需求的强劲牵引下，我国的航空重力测量技术取得了较大发展，重力向上和向下延拓技术也因此成为了这个领域的研究热点。由于航空重力测量成果是飞行高度面上的重力扰动或重力异常，因此在实际应用中，遇到更多的需求是将空中重力观测参量向下延拓到地面或大地水准面，以便联合其他类型观测数据进行大地测量应用计算[6-8]。但应用中也不可避免会遇到相反的地球外部重力场赋值需求[5]，一方面，外部空间航天器飞行重力场保障需要综合利用包括地面和海面在内的各类重力测量数据源；另一方面，开展航空重力测量系统技术性能评估与观测数据质量外部检核研究也需要重力向上技术的支持，即需要将已知的地面或海面重力向上延拓到飞行高度面，直接与航空重力测量成果作比对和分析，以便对其观测精度做出客观的评价[9,10]。
关于重力向上延拓问题，国内外学者已经提出了多种解决方案，包括直接数值积分法、梯度法、FFT方法、最小二乘配置法和等效源法等[11]，最传统最常用的解算模型是著名的Poisson积分方程[4]。Cruz等曾对Poisson积分计算模型作过深入研究，详细分析探讨了该模型在顾及地形改正、不顾及地形改正和采用移去-恢复技术等各种条件下的实际应用效果[12]；为了评估航空重力测量的精度，Argeseanu全面分析比较了包括Poisson积分公式在内的几种重力向上延拓模型的适用性和有效性[11]；石磐等提出通过Poisson积分公式计算空中重力异常，反过来评估联合应用航空重力测量和DEM数据确定地面重力场的内符合精度[13]；Kern在讨论卫星、航空和地面重力等多源数据融合处理问题时，也推荐使用Poisson积分计算模型[14]；王兴涛等在对航空重力测量数据向下延拓方法做比较分析时，使用Poisson积分公式作为地面重力向上延拓的计算模型[6]；翟振和等在讨论空中重力异常代表误差时，同样使用Poisson积分公式作为地面重力向上延拓的计算模型[15]。这里需要指出的是，由于Poisson积分模型是源于底律希勒（Dirichlet）边值问题的球面解，并未顾及地形高度起伏的影响，因此严格讲，该模型只适用于海域重力向上延拓解算[10]。而在实际应用中，国内外有不少作者往往忽略了该模型的适用条件，直接将其应用于陆部重力向上延拓解算，这样必然会给计算结果带来一定的偏差，其影响量值大小取决于计算区域重力场和地形变化的复杂程度。由于地形因素是客观存在的，因此在陆部要想达到一定的解算精度，必须考虑地形效应的作用。实际上，早在1966年Moritz就推出了根据地面重力异常和地形高数据计算外部空间重力异常的精密公式[5]。但因该模型计算过程过于复杂，就是现有的数据条件也难以保证计算结果精度，故在实践中一直未能得到推广应用。目前在实施重力向下延拓计算中，更多采用所谓“移去—恢复”技术，即首先移去地形质量对空中重力异常的影响，然后使用Poisson积分完成向下延拓计算，最后在计算结果中恢复地形质量的影响[7,16]。显然，在实施重力向上延拓计算时，也完全可以采用类似的思路来顾及地形效应的影响。但值得注意的是，移去地形引力对地面重力影响后，仍需要将地面重力异常残差向下延拓到海平面（可近似为球面）才能应用于Poisson积分计算，这一过程又涉及不规则地形面处理问题。可见，关于重力向上延拓精密计算，当前需要破解的技术瓶颈是，如何精确、有效地响应不规则地形质量对重力向上延拓计算结果的影响。由此引出的向上延拓不同解算模型的适用性和精度评估问题，应当是这个研究领域需要关注的重点。
⒉ 重力向下延拓技术
精化大地水准面一直是现代物理大地测量学的主要研究任务之一[17-18]。精密确定大地水准面需要联合利用卫星、航空、地面、海洋等多源重力和地形高信息，其计算过程不可避免涉及多源数据的融合处理问题。如前所述，向下延拓是多源重力数据融合处理必不可少的技术环节，在数据准备阶段，需要将卫星和航空重力测量成果延拓到地面作联合处理；在Stokes边值问题解算阶段，需要将地面重力数据延拓到大地水准面。由于物理场向下延拓在数学上属于不适定反问题，求解此类问题存在很大的不确定性[19]，一直以来国内外诸多学者为解决这一棘手问题付出了不懈的努力。目前解决向下延拓问题主要沿用三种途径，第一种是直接求逆Poisson积分方程，是应用比较广泛的主流途径[20]，国内外学者主要围绕求逆过程引起的不稳定性问题开展不同形式的正则化方法研究。Martinec曾从理论上全面分析研究了重力向下延拓过程的稳定性问题[21]；Xu提出使用截断奇异值分解（TSVD）方法求解不适定向下延拓问题[22]；Hansen等提出使用L-曲线法确定求解向下延拓问题的正则化参数[23]；沈云中等分析推导了不适定方程正则化算法的谱分解公式[24]；王兴涛等研究探讨了重力向下延拓的正则化算法及其谱分解式，同时分析比较了几种常用向下延拓方法的应用效果[6,25]；Kern、Mueller等、Alberts等学者也都曾对包括正则化和非正则化在内的不同类型向下延拓算法的适用性作过分析和比较[14,26,27]；王振杰给出了不适定问题正则化解法的统一表达式[28]；张辉等、刘东甲等提出了位场向下延拓的波数域算法及波数域迭代法[29,30]；顾勇为等分析研究了基于信噪比和参数分组修正的重力向下延拓正则化解法[31,32]；邓凯亮等提出了重力向下延拓的Tikhonov双参数正则化算法[33]；蒋涛等分析比较了Tikhonov、岭估计和广义岭估计三种正则化算法的技术特点及应用效果[34]；孙文等提出了基于波数域迭代的向下延拓Tikhonov正则化算法[35]；刘晓刚等研究探讨了重力和磁力测量数据向下延拓中最优正则化参数的确定问题[36]。当前解决向下延拓问题通常采用的第二种途径是间接求逆法，包括最小二乘配置方法[9,16,37]、虚拟点质量方法[17,18]和矩谐分析法[38]等，此类方法虽然避开了求逆Poisson方程，但仍涉及矩阵求逆过程，因此也不可避免存在不稳定性问题[39,40]。解决向下延拓问题的第三种途径可统称为非求逆方法，石磐等提出的利用航空重力测量和DEM数据确定地面重力场的直接代表法[13]，欧阳永忠、黄谟涛等提出的联合使用超高阶位模型和地形信息确定向下延拓改正数的方法(也称差分延拓法)等[41-43]都属于非求逆途径的范畴，此类方法不受传统求逆过程不稳定性的影响。Novak等基于带限（band-limited）航空重力测量数据固有的频谱特性，提出了向下延拓的直接积分公式(也称频谱截断法)[20]，并将其推广应用于大地水准面的直接解算[44]，其计算过程也避开了积分方程求逆难题。Kern将频谱截断方法和基于求逆过程的各类正则化方法作了全面的数值比较和分析[14]，得出的结论是：前者的计算精度和效率都明显优于其他方法。
如前所述，困扰向下延拓问题的关键是其本身固有的不适定性。Martinec曾对向下延拓的稳定性问题作过深入分析，同时提出了改善计算稳定性的地形效应补偿方法[21]；Kern、Alberts等、Novak等专门针对航空重力向下延拓问题，开展了大量卓有成效的数值计算和比对工作[14,201,27,44,45]，得出了一些极具参考价值的结论。由已有的研究成果得知，虽然向下延拓问题本身是不适定的，但其解算结果的不稳定程度取决于延拓计算高度和数据分辨率（即网格间距大小）两个方面。在一定条件下，向下延拓解算方程可以是良态的，而在超越一定界限以后，即使采用正则化处理技术，也无法取得有效的解算结果[14,45]。因此，对于实际应用来说，最要紧的是预先掌控前面所指的一定“条件”和“界限”，尽可能在规定的“条件”下开展航空重力测量作业，避免超越已知的“界限”，只有这样才能获得预期的重力测量归算成果。另一方面，地形效应对重力向下延拓解算结果具有重要影响，因为求逆过程的稳定性不仅取决于积分方程系数矩阵的结构，还取决于观测向量的频谱特性[21]，故可通过地形效应的“移去-恢复”运算来改变重力观测量的频谱特征，从而改善向下延拓解算的稳定性。刘敏等深入分析研究了当前国内外最具代表性的3种向下延拓计算模型的技术特点和适用条件，提出了应用超高阶位模型、局部地形改正和移去—恢复技术顾及地形效应，以及位场延拓结果球面化曲面的工程化方法，分析探讨了计算模型的稳定性及数据观测误差对延拓计算结果的影响，定量评估了不同向下延拓模型的解算精度及其可靠性，为下一步的工程化应用提供了具有可操作性的延拓计算方案[46]。
⒊ 多源数据融合技术
如前所述，进入21世纪，随着现代航空航天技术的飞速发展，地球重力场信息获取手段得到了全面拓展，目前已经形成了陆、海、空、天等全方位的地球重力场观测体系，重力测量数据种类日益丰富，观测精度逐步提高[17,18]。但随之而来我们也不得不面临这样一个问题，即如何有效地处理由不同测量手段、在不同界面获取的重力场观测数据。这些数据具有不同的频谱属性、分辨率、空间分布和误差特性。为了充分发挥各类数据资源的自身优势，准确刻画地球重力场的变化规律，必须采用现代数据处理手段，对存在异构的多源重力数据进行有效的融合处理，在融合过程中消除不同种类数据结构差异带来的矛盾，从而达到提高最终数据产品（也称数值模型）质量和可靠性的目的。
关于多源重力数据融合问题，国内外学者已经就此进行了广泛而深入的研究，提出了许多富有成效的处理方法，主要有统计法和解析法两大类型。最小二乘配置（简称配置法）是统计法的典型代表，由于该方法可以联合处理不同类型的重力数据，因此在多源重力数据融合处理中得到了广泛应用。配置法最早由Krarup提出[47]，后由Moritz等人发展和完善[37]，是20世纪60年代以来物理大地测量学在局部重力场逼近理论研究领域取得的重要进展。在此期间，Balmino、Schwarz、Rapp等学者都曾对配置法涉及的数学基础理论及其统计特性作过深入分析和研究[48-50]。在应用领域，Tscherning等研究探讨了采用配置法融合处理航空重力测量数据与地面重力的技术途径[51]；Strykowski等提出采用配置法解决丹麦格陵兰地区卫星、航空和地面重力数据之间的不协调性问题[52]；Tziavos等提出采用配置法融合处理卫星测高和海面船测重力数据[53]；Kern等分析研究了采用配置法融合处理卫星测高、航空和地面重力测量数据的实际效果[54]；邹贤才等研究探讨了配置法在局部重力场逼近计算中的适用性[55]；成怡开展了配置法在海洋多源重力数据融合处理中的应用研究[56]；翟振和将配置法应用于陆海交界区域多源重力数据的融合处理[57]；欧阳永忠、黄谟涛等提出了融合海域多源重力数据的Tikhonov正则化配置法[39,41]。从已有研究结果可以看出，协方差模型构建是配置法应用的前提和关键，尽管可以通过自适应调整模型参数的方式来改善协方差函数的特性[58]，但由于经验协方差模型的建立必须以足够分辨率的观测数据为基础，因此在实际应用中要想获得较高逼近度的协方差函数模型特别是三维空间协方差函数模型不是一件很容易的事，配置法的融合处理效果也因此受到了很大的制约。需要指出的是，协方差函数的拟合过程虽然可以近似的方式得以实现，但随之而来的问题是这种近似能否得到正确的解算结果，即拟合解的收敛问题，目前还无法得到证明。另外，基于随机过程的配置法与地球物理学事实相违背的是，地球重力场并不是一个严格意义的随机场[59]。尽管存在上述缺陷，配置法在多源数据融合处理中仍不失为一种可供选择的实用型方法。除配置法外，谱组合法和多输入/单输出法也属于统计法的范畴，它们在联合多种数据确定全球和局部重力场模型中得到了较好的应用[57,59-61]，Li等对采用配置法、频域最小二乘法和多输入/单输出法，联合卫星测高和船测重力数据恢复局部重力场的实际效果进行了全面分析和比较[62]。联合平差法和迭代法是多源数据融合处理解析方法的主要代表。刘晓刚分析研究了利用多种类型重力测量数据确定地球重力场的联合平差模型，提出了不同类型重力测量数据的最优权估计方法[59]；Kern、郝燕玲等依据残差重力异常修正重力位模型系数的思想，提出了融合卫星、航空、地面（海面）重力数据的迭代计算方法[14,63]；黄谟涛等、欧阳永忠提出的融合海域多源重力数据的正则化点质量方法也属于解析法的范畴[40,41]；欧阳永忠、黄谟涛等在深入分析海域多源重力数据误差特性基础上，提出了基于双权因子的多源数据网格化一步融合处理方法和基于分步平差、拟合、推估和内插相结合的多步融合处理方法[41,64]。这个研究领域下一步的关注重点，仍然是如何依据不同手段获取数据的异构性特点，构造出相匹配和最优化的数据融合处理模型。
三、海空重力测量数据综合应用技术
海空重力测量数据在国民经济建设、军事应用保障和地球科学研究等诸多领域都具有重要的应用价值。考虑到研究目标的指向性和研究内容的限定性，本文主要针对军事应用保障中的外部重力场赋值（即地球外部重力场逼近）和大地测量学研究中的大地水准面精化（即地球局部重力场逼近）两大应用主题作综合评述。
⒈ 地球外部重力场赋值技术
地球外部重力场赋值是海空重力测量数据应用于航天飞行器飞行轨迹控制保障的一个非常重要的方面。如文献[5]所述，外部扰动引力场计算是地球重力场逼近理论研究的主要内容之一，也是大地测量边值问题解算的最终目的。此项研究在空间科学研究领域具有重要的应用价值[4,18]。因为在地球外部空间飞行的人造卫星和各类飞行器始终受到地球扰动引力场的作用，其飞行轨道必然会偏离正常的椭圆轨道，即产生所谓的轨道摄动。为了精准确定飞行轨道摄动，必须精确计算飞行器在飞行过程中所受的地球扰动引力矢量，这正是开展外部重力场赋值研究的目的所在[18,65]。而这项工作的数据基础就是由包括海空重力测量资料在内的多源数据生成的不同分辨率数值模型。
当前应用于外部重力场赋值的计算模型主要有三类：一类是基于全球重力位模型的球谐展开式[4,18,65,66]；另一类是基于经典Stokes边值理论的直接积分模型[4,18,65,67]；第三类是基于现代Bjerhammar换置理论的点质量模型[18,68,69]。黄谟涛、吴晓平、赵东明等、张嗥、郑伟等、江东等、马彪等人的研究结果表明[18,69-74]，三类模型各具不同的技术特征和适用条件，球谐展开模型主要用于扰动重力的长波段计算，但受计算时间的限制，球谐展开式的阶次不宜取得太高；积分模型不需要对地面重力数据作进一步的转换处理，但其在超低空存在积分奇异性和无法顾及地形效应的缺陷，在一定程度上限制了该模型的应用；点质量模型对超低空扰动重力场具有较强的恢复能力，但赋值之前需要完成地面重力数据到点质量的转换，其中涉及较大规模线性方程组的解算问题。就军事应用保障而言，远程飞行器对重力场的保障需求主要体现为扰动引力矢量计算精度和速度两个方面，在确保精度指标要求前提下，如何突破扰动引力的计算速度和稳定性指标，一直是这个研究领域的关注点。为了满足扰动引力快速赋值的工程化要求，人们进而研究提出了扰动引力数值逼近方法，即通过前期的数值建模，构建扰动引力与空间位置一一对应的简单化关系，来进一步缩短后期扰动引力的计算时间。赵东明等分析研究了扰动引力快速赋值的有限元逼近法[75]；张嗥、王建强等提出了扰动引力多项式拟合法[71,76]；郑伟等、周世昌等进一步提出了广域多项式拟合法[72,77]；王继平等研究探讨了扰动引力赋值的神经网络逼近算法[78]。这个研究领域的最新进展是，针对直接积分模型在计算低空扰动引力时出现的数值不稳定性和积分奇异性问题，刘长弘提出了四种改进型Stokes积分算法：即去中央奇异点算法、中央格网加密算法、奇异点积分值修正算法和改进积分式算法[79]；根据扰动引力场变化特征和广域多项式计算特点，常岺提出了构建覆盖全球的空间分层扰动引力赋值模型[80]；翟曜研究分析了利用附加参数最小二乘配置方法进行扰动引力逼近计算的实际效果[81]。
这里需要指出的是，对于海上军事应用重力场保障需求，由于海域飞行器发射阵位具有高度机动性和覆盖范围广阔的显著特点，因此，与陆地固定发射阵位相比较，远程飞行器对海域外部扰动引力场赋值的技术要求更高，特别是数据准备阶段的实时性要求更高。针对这种特殊的需求，必须在深入分析现有赋值模式技术特点的基础上，探讨联合应用多种模式构建一种能够兼备各自模式优良特性的综合模型的可能性，同时从工程化应用角度出发，继续对传统赋值模式进行改进和优化，为建立完善的远程飞行器重力场保障体系提供必要的技术支撑。
⒉ 大地水准面精化技术
不断精化大地水准面始终是现代物理大地测量学永恒的研究主题之一[17,18]，同时也是海空重力测量数据应用的重点方向[42,82,83]。大地水准面，一方面作为最贴近地球形状的一个封闭重力等位面，在地球科学相关学科研究中具有重要的科学意义；另一方面作为高程系统的起算面，在统一全球空间地理信息基准框架方面具有不可替代的实用价值[17,18,84]。随着高精度GNSS测高技术的发展，“GNSS+大地水准面模型”技术已经从根本上改变了传统高程基准的维持模式和高程测定的作业模式。正如李建成院士所述[85]：新的维持模式是一种无需建立地面标石的“绿色模式”，新的作业模式是一种地表“无障碍模式”，也是一种相对“独立测高模式”。精密确定大地水准面模型因此已经成为当前全球高程基准现代化基础设施建设的核心任务之一。
由物理大地测量学知[4,17]，确定重力大地水准面，理论上归结为求解相对于一个正常重力场的扰动位函数，即求解人们熟知的大地测量边值问题。经过100多年的研究和发展，边值问题经历了以三类边值理论为代表的三个主要发展阶段，即以Stokes理论为代表的大地水准面边值理论发展阶段[4]，以Molodensky理论为代表的地球自然表面边值理论发展阶段[38]，以及以Bjerhammar理论为代表的虚拟球面边值理论发展阶段[68]。Stokes理论因涉及大地水准面外部质量调整和地形质量密度假设问题，故也被称为调整大地水准面边值问题。与上述三类理论相对应的边值问题解式分别为：Stokes积分、Molodensky级数解和Bjerhammar-广义Stokes积分。从理论上讲，只要地形重力归算及其间接影响考虑周密，各类不同的地形质量调整方案都能给出比较一致的Stokes积分解[4]。目前在实践中应用比较广泛的地形质量调整方案是Helmert第二类凝集法，其相对应的边值问题解被称为Stokes-Helmert方法[85]。Stokes边值问题解存在的主要问题是：地形质量密度假设偏差带来的不严密性。Molodensky边值理论不要求调整地形质量，因此Molodensky级数解在理论上是严密的，但由于高阶项计算过程的复杂性和不稳定性，实践中一般只考虑到级数的一阶项[17,86]，至多考虑到级数的二阶项[87]。考虑一阶项时通常使用局部地形改正代替重力改正，即相当于使用Faye异常的Stokes积分解，这是基于空间重力异常与高程存在线性关系的基本假设下所做的近似处理[17,37,88]。然而在地形起伏较大的山区，重力异常与地形高度之间未必严格满足简单的线性关系，因此这种近似处理必然会带来一定的偏差[89]。Bjerhammar边值理论也不存在地形质量调整问题，不需要做地形效应改正和补偿处理。实际上，Bjerhammar理论可以看作是一类广义的“等效场源”方法，Bjerhammar提出的虚拟重力异常法[90]，Weightman提出的虚拟点质量法[91,92]，许厚泽等提出的虚拟单层密度法[93]，申文斌提出的虚拟压缩恢复法[94]，甚至包括Moritz提出的解析延拓法[37]等，都是建立在一种所谓的位场“等效原理”之上的，因此它们都可以归类为上述的“等效场源”方法。这类方法不需要调整地形质量，但在地面观测与等效场源的转换过程中已经巧妙地顾及了地形效应影响，因此具有理论上的严密性。这类方法在实用上的不足是，求解等效场源都包含一种“逆”过程，即涉及地面观测向下延拓的不适定性问题，这是此类方法在数学结构上的“先天性”共同弱点[17]。如何有效克服这一弱点，是推广应用“等效场源”方法的关键所在[40,95]。
这里需要补充说明的是，对于传统的Stokes边值问题解，考虑到无论采用哪一种类型的地形质量调整方案，最终都需要将经过“改正”后的地面重力异常向下延拓到大地水准面[4]，因此其转换过程也必然会遇到与Bjerhammar“等效场源”解法相类似的不适定反问题。尽管与“原始”的地面重力异常相比，经过“改正”后的地面重力异常可能要变得平缓许多，此时的延拓空间也不再存在地形质量（其前提条件是地形质量密度假设没有偏差），即能够保证这个阶段的向下延拓是正则的，但向下延拓问题本身仍然是一种“逆”过程，并不改变固有的不适定特性，对地面重力观测误差仍会有放大作用，因此，与Moritz解析延拓法相类似，其解算过程的不稳定性问题仍然需要关注。另外，关于“移去-恢复”过程中参考场的使用问题，不同的边值解算方法应当使用不同的边界面计算方案。对于未经调整的大地水准面Stokes边值问题解，既可以移去地形面上的参考场重力异常，也可以移去大地水准面上的参考场重力异常，前者意味着已经使用位模型将地面重力异常向下延拓到大地水准面，后者则意味着把地面重力异常等同于大地水准面重力异常，但两者都必须在大地水准面上恢复参考场。对于调整的大地水准面Stokes边值问题解，则只能在大地水准面上完成移去和恢复参考场运算。对于Molodensky边值问题级数解（包括使用Faye异常的Stokes积分解），也只能在大地水准面上完成移去和恢复参考场运算，因为级数解中的高阶项已经隐含对地面重力异常的向下延拓和大地水准面的向上延拓，如果在此情形下仍采用在地形面上移去和恢复参考场的计算方案，那么就相当于重复考虑了地形效应的影响。对于Bjerhammar边值问题解，则必须在地形面上移去参考场重力异常，在地形面上恢复参考场似大地水准面（高程异常），或在Bjerhammar球面上恢复参考场重力异常。对于Moritz解析延拓解，既可以移去地形面上的参考场重力异常，也可以移去大地水准面上的参考场重力异常，前者意味着分步将地面重力异常向下延拓到大地水准面，第一步使用位模型做近似延拓，第二步使用地面重力异常与地形面位模型重力异常的残差做精细延拓；后者则意味采用一步法即使用地面重力异常与大地水准面位模型重力异常的残差，将地面重力异常向下延拓到大地水准面。相比较而言，在相同的数据条件下，前者的延拓计算精度应当略优于后者。但两者都必须在地形面上恢复参考场似大地水准面（高程异常）。
在大地水准面大规模计算实践中，人们早期应用较多的是传统的Stokes积分解，后来逐步推广使用基于Faye异常的Stokes积分解，也就是Molodensky级数零阶加一阶项解，近期使用更多的是Stokes-Helmert方法和Moritz解析延拓法。美国在20世纪90年代以后，每隔3至5年就更新换代一次国家高程基准（NAVD88）大地水准面模型，具体数值模型包括GEOID90、GEOID93、G9501、GEOID96、GEOID99、USGG2003和USGG2009等[96-100]。其中，USGG2003及其之前的各个大地水准面模型均采用基于Faye异常的Stokes积分解法，不同的只是在计算细节和数据使用上做了一些改进和补充。新的USGG2009模型则采用了“超高阶位模型+剩余地形模型（RTM）”移去恢复技术加解析延拓解，Wang等认为[100]，在各个改正量都得到精确计算的条件下，该解法与严密的Stokes-Helmert方法是等价的，文献[85]也因此认为USGG2009模型采用的就是严密的Stokes-Helmert方法。USGG2009模型分辨率为1¢´1¢,与GPS水准观测量作比较，其差值的标准差为6.3cm，剔除长波长误差后的差值标准差为4.3cm[100]。欧洲地区大地水准面的计算工作始于20世纪80年代初，第一代欧洲大地水准面模型EGG1的精度为dm级，分辨率约为20km。1990年，IAG大地水准面欧洲分委会启动新一轮欧洲大地水准面精化计划，并从1994年开始，先后推出了EGG94、EGG95、EGG96和EGG97系列欧洲大地水准面模型[101-104]。计算这些数值模型使用的基本方法是“位模型+剩余地形模型（RTM）”移去恢复法。澳大利亚也曾在不同时期使用顾及Molodensky一阶项的级数解方法建立不同版本的国家大地水准面模型（如AUSGeoid98、AUSGeoid09等）[105]。我国于20世纪70年代完成第一代对应于1954北京坐标系的似大地水准面模型CQG60的计算，并于80年代初将该模型转换到新建立的1980西安大地坐标系[85]；90年代初我国推出第二代似大地水准面模型WZD94[106]；21世纪初我国推出新一代似大地水准面模型CQG2000[17]；2012年，我国推出最新一代重力似大地水准面模型CNGG2011[85]。构制CQG2000模型使用的计算方案是基于Faye异常的Stokes积分解，CQG2000模型的分辨率为5¢´5¢，在全国范围内与GPS水准点比较的精度为0.44m；构制CNGG2011模型采用的计算方案是严密的Stokes-Helmert方法，CNGG2011模型的分辨率为2¢´2¢，在全国范围内与GPS水准点比较的精度为0.13m。总体上讲，我国精化大地水准面采取的技术途径与美国基本保持一致，但在模型分辨率和精确度方面仍有一定的差距。下一步除了要在困难地区加大重力场信息获取资源投入外，应当继续关注地形和重力等多源数据的融合处理与应用，以及基于地形面的边值问题精细化解算理论和方法研究。数据源是提升我国似大地水准面精细化水平的根本，数学建模和解算方法是提高模型计算精度的关键[107]。
四、结束语
数值模型是海空重力测量数据产品的一种最常用的表现形式，是连接外业测量成果和应用开发系统的重要纽带。地球外部重力场赋值和大地水准面精化是海空重力测量数据最重要的两个应用领域，也是海空重力测量成果应用价值的最终体现。本文分析总结了海空重力测量数据向上和向下延拓技术的研究现状及发展方向，简要论述了海空重力测量多源数据融合处理技术的研究动态及发展前景，分析讨论了海空重力测量数据应用于地球外部重力场赋值和大地水准面精化技术的研究进展及发展思路。上述各个环节都属于海空重力测量技术体系的后端部分，即数值模型构建及综合应用部分，它们与前端的传感器及测量规划设计技术，中端的数据处理及精度评估技术，一起构成完整的海空重力测量技术体系。
参考文献
[1]刘敏,黄谟涛,欧阳永忠,等.海空重力测量及应用技术研究进展与展望（一）:目的意义与技术体系[J].海洋测绘,2017,37(2):1-6.
[2]刘敏,黄谟涛,欧阳永忠,等.海空重力测量及应用技术研究进展与展望（二）:传感器与测量规划设计技术 [J].海洋测绘,2017，37(3):1-6.
[3]刘敏,黄谟涛,欧阳永忠,等.海空重力测量及应用技术研究进展与展望（三）:数据处理与精度评估技术 [J].海洋测绘,2017,37(4):1-6.
[4]Heiskanen W A,Moritz H.Physical Geodesy [M].San Francisco: Freeman W H and Company,1967.
[5]Moritz H. The Computation of the External Gravity Field and the Geodetic Boundary-Value Problem [A].In :Proceedings of the Symposium“Extension of Gravity Anomalies to Unsurveyed Areas”[C], Am. Geophys.Union,Geophysical Monograph Series,1966,No.9:127-136.
[6]王兴涛,夏哲仁,石磐,等.航空重力测量数据向下延拓方法比较[J].地球物理学报，2004,47(6):1017-1022.
[7]Tziavos I N,Andritsanos V D,Forsberg R, et al.Numerical Investigation of Downward Continuation Methods for Airborne Gravity Data[A]. IAG Symposia 129: Gravity, Geoid and Space Missions(GGSM)[C]. Springer,Heidelberg,Berlin,2005:119-124.
[8]吴太旗,邓凯亮,黄谟涛,等.一种改进的不适定问题奇异值分解法[J].武汉大学学报·信息科学版，2011,36（8）:900-903.
[9]Hwang C,Hsiao Y S,Shih H C,et al.Geodetic and Geophysical Results from a Taiwan Airborne Gravity Survey: Data Reduction and Accuracy Assessment[J].Journal Geophysical Research,2007,112(B10), doi:10.1029/2005JB004220.
[10]翟振和,孙中苗,李迎春,等.航空重力测量在近海区域的精度评估与分析[J].测绘学报,2015,44（1）:1-5.
[11]Argeseanu V S.Upward Continuation of Surface Gravity Anomaly Data[C].In: Proceedings of the IAG Symposium on Airborne Gravity Field Determination,1995,95-102.
[12]Cruz J Y,Laskowski P.Upward Continuation of Surface Gravity Anomalies[C]. Report No.360, Dept. of Geodetic Science and Surveying,the Ohio State University,1984.
[13]石磐,王兴涛.利用航空重力测量和DEM确定地面重力场[J].测绘学报,1997,26(2):117-121.
[14]Kern M.An Analysis of the Combination and Downward Continuation of Satellite, Airborne and Terrestrial Gravity Data [D].UCGE Reports，University of Calgary,2003,Number 20172.
[15]翟振和,李超,李红娜.空中重力测量数据代表地面数据的误差分析[J].海洋测绘,2012,32（2):1-3.
[16]Forsberg R,Olesen A V.Airborne Gravity Field Determination [A].in: Sciences of Geodesy[C], Springer-Verlag Berlin Heldelberg,2010,83-103.
[17]李建成,陈俊勇,宁津生,等.地球重力场逼近理论与中国2000似大地水准面的确定[M].武汉:武汉大学出版社,2003.
[18]黄谟涛,翟国君,管铮,等.海洋重力场测定及其应用[M].北京:测绘出版社,2005.
[19]王彦飞.反演问题的计算方法及其应用[M].北京:高等教育出版社,2007.
[20]Novak P,Heck B.Downward Continuation and Geoid Determination Based on Band-limited Airborne Gravity Data [J].Journal of Geodesy,2002,76:269-278.
[21]Martinec Z.Stability Investigations of a Discrete Downward Continuation Problem for Geoid Determination in the Canadian Rocky Mountains[J]. Journal of Geodesy, 1996,70:805-828.
[22]Xu Peiliang.Truncated SVD Methods for Discrete Linear Ill-posed Problems[J]. Geophys. J. Int., 1998, 135:505-514.
[23]Hansen P C, O’Leary D P.The Use of the L-curve in the Regularization of Discrete Ill-posed Problems [J].SIAM J.Sci.Comput.,1993,14(6):1487-1503.
[24]沈云中,许厚泽.不适定方程正则化算法的谱分解式[J].大地测量与地球动力学,2002,23(3):10-14.
[25]王兴涛,石磐,朱非洲.航空重力测量数据向下延拓的正则化算法及其谱分解[J].测绘学报,2004,33（1):33-38.
[26]Mueller F,Mayer-Guerr T.Comparison of Downward Continuation Methods of Airborne Gravimetry Data [A].IAG Symposia 128:A Window on the Future of Geodesy[C].Springer,Heidelberg,Berlin, 2003:254-258.
[27]Alberts B,Klees R.A Comparison of Methods for the Inversion of Airborne Gravity Data[J]. Journal of Geodesy,2004,78:55-65.
[28]王振杰.测量中不适定问题的正则化解法[M].北京:科学出版社,2006.
[29]张辉,陈龙伟,任治新,等.位场向下延拓迭代法收敛性分析及稳健向下延拓方法研究[J].地球物理学报，2009,52(4):1107-1113.
[30]刘东甲,洪天求,贾志海,等.位场向下延拓的波数域迭代法及其收敛性[J].地球物理学报,2009,52(6):1599-1605.
[31]顾勇为,归庆明.航空重力测量数据向下延拓基于信噪比的正则化方法的研究[J].测绘学报,2010,39(5):458-464.
[32]顾勇为,归庆明,韩松辉,等.航空重力向下延拓分组修正的正则化解法[J].武汉大学学报·信息科学版,2013,38(6):720-723.
[33]邓凯亮,黄谟涛,暴景阳,等.向下延拓航空重力数据的Tikhonov双参数正则化法[J].测绘学报,2011,40(6):690-696.
[34]蒋涛,李建成,王正涛,等.航空重力向下延拓病态问题的求解[J].测绘学报,2011,40(6):684-689.
[35]孙文,吴晓平,王庆宾,等.航空重力数据向下延拓的波数域迭代Tikhonov正则化方法[J].测绘学报,2014,43(6):566-574.
[36]刘晓刚,李迎春,肖云,等.重力和磁力测量数据向下延拓中最优正则化参数确定方法[J].测绘学报,2014,43(9):881-887.
[37]Moritz H.Advanced Physical Geodesy [M].Abacus Press,England,1980.
[38]蒋涛,党亚民,章传银,等.基于矩谐分析的航空重力向下延拓[J].测绘学报,2013,42(4):475-480.
[39]黄谟涛,欧阳永忠,翟国君,等.融合多源重力数据的Tikhonov正则化配置法[J].海洋测绘,2013,33(3):6-12.
[40] 黄谟涛，欧阳永忠，刘敏，等. 融合海域多源重力数据的正则化点质量方法[J]. 武汉大学学报·信息科学版，2015，40(2）：170-175.
[41]欧阳永忠.海空重力测量数据处理关键技术研究[D].武汉:武汉大学,2013.
[42]黄谟涛,欧阳永忠,刘敏,等.海域航空重力测量数据向下延拓的实用方法[J].武汉大学学报·信息科学版,2014,39(10):1147-1152.
[43]黄谟涛,宁津生,欧阳永忠,等.联合使用位模型和地形信息的陆区航空重力向下延拓方法[J].测绘学报,2015,44(4):355-362.
[44]Novak P,Kern M,Schwarz K P,et al.On Geoid Determination from Airborne Gravity [J].Journal of Geodesy,2003,76:510-522.
[45]Novak P,Kern M,Schwarz K P.Numerical Studies on the Harmonic Downward Continuation of Band-limited Airborne Gravity[J].Studia Geophysica et Geodaetica,2001,45:327-345.
[46]刘敏,黄谟涛,欧阳永忠,等.顾及地形效应的重力向下延拓模型分析与检验 [J].测绘学报,2016, 45(5):521-530.
[47]Krarup T.A Contribution to the Mathematical Foundation of Physical Geodesy [R].Publ.44,Dan. Geod. Inst.,Copenhagen,1969.
[48]Balmino G.Introdution to Least-Squares Collocation [A].In:Approximation Methods in Geodesy[C],edited by H. Moritz and H. Sunkel,Herbert Wichmann Verlag Karlsruhe,1978.
[49]Schwarz K P.On the Application of Least-Squares Collocation Models to Physical Geodesy [A]. In: Approximation Methods in Geodesy[C],edited by H. Moritz and H. Sunkel, Herbert Wichmann Verlag Karlsruhe,1978.
[50]Rapp R H.Results of the Application of Least-Squares Collocation to Selected Geodetic Problems [A].In:Approximation Methods in Geodesy[C],edited by H. Moritz and H. Sunkel, Herbert Wichmann Verlag Karlsruhe,1978.
[51]Tscherning C C,Rubek F,Forsberg R. Combining Airborne and Ground Gravity Using Collocation[C]. IAG Symposia,Vol.119，Springer-verlag,1997,18-23.
[52]Strykowski G,Forsberg R.Operational Merging of Satellite, Airborne and Surface Gravity Data by Draping Techniques [C].IAG Symposia,Vol.119,Springer-verlag,1997,243-248.
[53]Tziavos I N,Sideris M G,Forsberg R.Combined Satellite Altimetry and Shipborne Gravimetry Data Processing [J].Marine Geodesy,1998,21:299-317.
[54]Kern M,,Schwarz K P,Sneeuw N.A Study on the Combination of Satellite,Airborne and Terrestrial Gravity Data[J].Journal of Geodesy,2003,77:217–225.
[55]邹贤才,李建成.最小二乘配置方法确定局部大地水准面的研究[J].武汉大学学报·信息科学版,2004,29(3):218-222.
[56]成怡.多源海洋重力数据融合技术研究[D].哈尔滨:哈尔滨工程大学,2008.
[57]翟振和.陆海交界区域多源重力数据的融合处理方法研究[D].郑州:解放军信息工程大学,2009.
[58]翟振和,孙中苗.基于配置法的局部重力场延拓模型构建与应用分析[J].地球物理学报,2009,52(7):1700-1706.
[59]刘晓刚.GOCE卫星测量恢复地球重力场模型的理论与方法[D].郑州:解放军信息工程大学,2011.
[60]Wenzel H G.Geoid Computation by Least Squares Spectral Combination Using Integral Kernels[C]. Proceedings of Symposium No.46“Geoid Definition and Determination”,General Meeting of IAG, Tokyo,1982：438-453.
[61]吴晓平,陆仲连.卫星重力梯度向下延拓的最佳积分核组合解[J].测绘学报,1992,21(2):121-131.
[62]Li J, Sideris M G.Marine Gravity and Geoid Determination by Optimal Combination of Satellite Altimetry and Shipborne Gravimetry Data [J].Journal of Geodesy,1997,71:209-216.
[63]郝燕玲,成怡,刘繁明,等.融合多类型海洋重力数据算法仿真研究[J].系统仿真学报,2007,19(21):4897-4900.
[64]黄谟涛,欧阳永忠,翟国君,等.海域多源重力数据融合处理的解析方法.武汉大学学报·信息科学版，2013,38(11):1261-1265.
[65]陆仲连,吴晓平,丁行斌,等.弹道导弹重力学[M].北京:八一出版社,1993.
[66]王建强,李建成,王正涛,等.球谐函数变换快速计算扰动引力[J].武汉大学学报·信息科学版，2013,38(9):1039-1043.
[67]吴晓平.在推求地球外部扰动重力场中数据的采用[J].解放军测绘学院学报,1992,(4):1-10.
[68]Bjerhammar A. A New Theory of Geodetic Gravity[R].Transaction of the Royal Institute of Technology,No.243,Stockholm,1964.
[69]吴晓平.局部重力场的点质量模型[J].测绘学报,1984,13(4):249-258.
[70]赵东明,吴晓平.扰动引力快速确定的替代算法[J].解放军测绘学院学报,2001,18(B09):11-13.
[71]张嗥.快速逼近弹道扰动引力的算法研究[D].郑州：解放军测绘学院，2007.
[72]郑伟,钱山,汤国建.弹道导弹制导计算中扰动引力的快速赋值[J].飞行力学,2007,25(3):42-44.
[73]江东,王庆宾,赵东明.空中扰动引力快速赋值算法的效能分析[J].测绘科学技术学报,2011, 28(6):411-415.
[74]马彪,刘晓刚,张丽萍.几种弹道扰动引力计算模型的分析与比较[J].大地测量与地球动力学,2012, 32(1):105-109.
[75]赵东明,吴晓平.利用有限元方法逼近飞行器轨道主动段扰动引力[J].宇航学报,2003,24(3):309-313.
[76]王建强,李建成,赵国强,等.多项式拟合快速计算扰动引力方法[J].大地测量与地球动力学,2013,33(4):52-55.
[77] 周世昌，王庆宾，张传定. 快速确定扰动引力的广域多项式逼近方法的模拟实验[J]. 测绘科学，2009,34（4）：153-154.
[78]王继平,王明海,张志辉.扰动引力的神经网络逼近算法[J].宇航学报,2008,29(1):385-390.
[79]刘长弘.改进的直接积分方法计算低空扰动引力[D].郑州:解放军信息工程大学,2016.
[80]常岺.空间分层扰动引力场快速构建算法研究[D].郑州:解放军信息工程大学,2016.
[81]翟曜.基于附加参数最小二乘配置法的扰动引力逼近方法研究[D].郑州:解放军信息工程大学,2016.
[82]孙中苗.航空重力测量理论、方法及应用研究[D].郑州:解放军信息工程大学,2004.
[83]孙中苗,翟振和,肖云.渤海湾航空重力及其在海域大地水准面精化中的应用[J].测绘学报,2014,43(11):1101-1108.
[84]李建成.我国现代高程测定关键技术若干问题的研究及进展[J].武汉大学学报·信息科学版,2007,32(11):980-987.
[85]李建成.最新中国陆地数字高程基准模型：重力似大地水准面CNGG2011[J].测绘学报,2012,41(5):651-660,669.
[86]陈俊勇,李建成,宁津生,等.中国新一代高精度、高分辨率大地水准面的研究与实施[J].武汉大学学报·信息科学版,2001,26(4):283-289.
[87]章传银,晁定波,丁剑,等.厘米级高程异常地形影响的算法及特征分析[J].测绘学报,2006,35(11):308-314.
[88]Wang Y M.Comments on Proper Use of the Terrain Correction for the Computation of Height Anomalies[J].Manuscripta Geodaetica,1993,18:53-57.
[89]李姗姗，吴晓平，张传定，等.顾及地形与完全球面布格异常梯度项改正的区域似大地水准面精化[J].测绘学报,2012,41(4):510-516.
[90]Bjerhammar A.Discret Physical Geodesy[R].Dept.of Geodetic Science and Surveying,The Ohio University,Rep.No.380,1987.
[91]Sunkel H.The Generation of a Mass Point Model from Surface Gravity Data[R].Dept.of Geodetic Science and Surveying,The Ohio University,Rep 353,1983.
[92]Antunes C,Pail R,Catalao J.Point Mass Method Applied to the Regional Gravimetric Determination of the geoid[J]Studia geophysica et geodaetica,2003,47:495-509.
[93]许厚泽,朱灼文.地球外部重力场的虚拟单层密度表示[J].中国科学(B辑),1984,6:575-580.
[94]申文斌,宁津生,晁定波.边值问题虚拟压缩恢复原理及其在Bjerhammar理论中的一个应用[J].测绘学报,2005,34(1):14-18.
[95]束蝉方,李斐,李明峰,等.应用Bjerhammar方法确定GPS重力似大地水准面[J].地球物理学报,2011,54(10):2503-2509.
[96]Milbert D G.Computing GPS-derived Orthometric Heights with the GEOID90 Geoid Height Model [C].In:Technical Papers of the 1991 ACSM-ASPRS Fall Convention, American Congress on Surveying and Mapping,1991,46-55.
[97]Smith D A,Milbert D G.The GEOID96 Hight-resolution Geoid Height Model for the United States [J].Journal of Geodesy,1999,73:219-236.
[98]Smith D A,Roman D R.GEOID99 and G99SSS:One Arc-minute Models for the United States [J].Journal of Geodesy,2001,75:469-490.
[99]Roman D R.GEOID99 and G99SSS:One Arc-minute Models for the United States [J].Journal of Geodesy,2001,75:469-490.
[100]Wang Y M,Saleh J,Li X,Roman D R.The US Gravimetric Geoid of 2009 (USGG2009): Model Development and Evaluation[J].Journal of Geodesy,2012,86:165-180.
[101]Denker H,Behrend D,Torge W.European Gravimetric Geoid: Status Report[A].In:Proceedings of IAG Symposium,Vol.113,423-433.
[102]Denker H,Behrend D,Torge W.The European Gravimetric Quasigeoid EGG95[J].IGES Bulletin,No.4, 1996a,3-11.
[103]Denker H,Behrend D,Torge W.The European Gravimetric Quasigeoid EGG96[A].In:Proceedings of IAG Symposium,Vol.117,1996b,532-539.
[104]Denker H,Torge W,Wenzel G.Investigation of Different Methods for the Combination of Gravity and GPS/Leveling Data[A].In:Proceedings of IAG Symposium,Vol.121,1999,137-142.
[105]Featherstone W E,Kirby J E, Hirt C,et al.The AUSGeoid09 Model of Australian Height Datum [J].Journal of Geodesy,2011,85:133-150.
[106]管泽霖,李建成,晁定波,等. WZD94中国大地水准面研究[J].武汉测绘科技大学学报,1994,19(4):292-297.
[107]陈欣.海域重力测量数据的精化处理研究[D].大连:海军大连舰艇学院,2016.
【作者简介】第一作者刘敏，1980出生，男，湖南衡阳人，博士研究生，主要从事海洋重力场测定理论方法及应用研究；本文为基金项目，包括国家重大科学仪器设备开发专项（2011YQ12004503）、国家自然科学基金（41474012,41374018）、国防973计划（613219）、国家重点研发计划（2016YEC0303007,2016YFB0501704)；本文来自《海洋测绘》（2017年第5期），若其他公众平台转载，请备注论文作者，并说明文章来源，版权归《海洋测绘》所有。
相关阅读推荐
论文专区▏海空重力测量及应用技术研究进展与展望（二）：传感器与测量规划设计技术
公众号
溪流之海洋人生
微信号▏xiliu92899
用专业精神创造价值
用人文关怀引发共鸣
您的关注就是我们前行的动力
投稿邮箱▏452218808@qq.com