重力异常垂向梯度严密改化模型及应用

黄谟涛，邓凯亮，吴太旗，欧阳永忠，陈欣，翟国君，刘敏，王许

海军研究院，天津 300061

0 引言

重力异常场是地球内部质量分布不均匀性及外部地形变化形态的综合反映（Dehlinger,1978；Torge,1989）.利用地表重力观测研究确定地球形状大小及变化特性，是重力与大地测量学的核心任务,其研究内容包含大地测量边值问题解算、局部和外部重力场逼近（Heiskanen and Moritz,1967；李建成等,2003；Sansò and Sideris,2013）；依据地表重力观测研究确定地球内部物质结构、形态及物理特性，则是地球物理学的重要研究主题之一，其研究内容涉及利用数值计算方法解决地球重力场的正演和反演问题（Moritz,1990；马在田等,1997）.由此可见，地球重力场研究不仅跨越了地球物理学和大地测量学两个学科的研究领域，同时跨越了地球内部、地表及外部空间三个不同特征的研究空域，其研究内容除了涉及不同空域的多维度重力场建模技术外，还包含地表重力观测数据的向上和向下延拓归算技术（Moritz,1980；Cruz and Laskowski,1984；梁锦文,1989；王兴涛等,2004；刘敏等,2018）.实际上，大地测量边值问题解算和地球外部重力场逼近计算都可以归结为广义的重力位场向上和向下延拓问题（Moritz,1966；Heiskanen and Moritz,1967；Moritz,1980；于锦海等,2001；Sansò and Sideris,2017）.可见，重力向上和向下延拓技术在地球重力场研究中具有非常重要的应用价值.实施向上和向下延拓解算除了可以依托传统的Poisson积分方程外（Heiskanen and Moritz,1967），还可以采用计算稳定性更好的Taylor级数展开模型（Peters,1949）.精确求取位场各阶垂向导数（也称垂向梯度）是运用Taylor级数展开延拓法的关键，国内外诸多学者为此提出了许多解决该问题的方法（Moritz,1980；刘保华等,1995；姚长利等,1997；Fedi and Florio,2002；Wei,2014；张冲等,2017；黄谟涛等,2018；Tran and Nguyen，2020）.实际上，除了向上和向下延拓解算，位场垂向梯度还可以直接应用于地质结构反演和地球资源勘探，因为重力梯度数据更能精细反映地球内部浅层地质结构密度的变化形态（王虎彪等,2009；刘金钊等,2013,2020）.

在前述各种应用中，重力异常一阶垂向梯度是应用面最广泛的重力场特征参数之一，在各类重力归算和大地测量边值解算中发挥着关键性作用（Heiskanen and Moritz,1967；Moritz,1980；Wang et al.,2012）.求解观测重力异常的全球积分是获取重力异常垂向梯度的主要手段，但实施此类积分模型数值计算时均涉及全球积分模型改化问题，一方面需要对全球积分计算式进行积分域分割处理，通常做法是将全球积分划分为近区和远区，近区采用实测数据进行数值积分计算，远区则采用地球位模型进行补偿（黄谟涛等,2020）；另一方面需要对积分核函数进行去奇异性处理，传统做法是通过引入合适的积分恒等式变换，将原积分计算模型转换为具有稳定数值解的连续函数模型（Heiskanen and Moritz,1967）.需要指出的是，在实际应用中，人们在实施全球积分域分割处理时，往往会忽视积分恒等式成立的全球积分条件，不再关注采用局域积分对积分恒等式带来的数值影响（Heiskanen and Moritz,1967；Moritz,1980；欧阳明达等,2014；翟振和等,2015；刘长弘,2016），从而引起不可忽略的计算误差（黄谟涛等，2019）.考虑到重力异常一阶垂向导数是计算二阶及其更高阶导数的基础，本文主要针对球面及外部重力异常一阶垂向导数全球积分模型改化问题进行分析研究和试验验证，依据实测数据保障条件和积分域分割处理方式，分别推出两类参数全球积分模型的严密改化公式，同时通过数值计算和向下延拓应用对比分析，进一步验证采用严密改化模型的可行性和有效性.

1 计算模型分析与改化

1.1 地球外部重力异常垂向梯度严密改化模型

由Heiskanen和Moritz（1967）得知，地球外部重力异常全球积分计算式可表示为

（1）

（2）

由式（1）不难看出，当计算点趋近于数据点时，即当r→R和ψ→0时，会出现分子分母项同时为零的不定式（0/0）情况，积分核函数K（r,ψ）发生奇异.为了消除积分式（1）的奇异性，通常做法是从积分域中直接扣除计算点所在的网格数据块，从而避免出现积分核函数分母项l→0现象.这种处理方法看似比较简单明了，但往往会给计算结果带来不可忽略的误差.为此，Heiskanen和Moritz（1967）建议采用式（3）表示的积分恒等式对式（1）进行改化：

（3）

略去推导过程，直接写出改化公式如下：

（4）

式中，ΔgRp代表空间计算点P（r,φ,λ）在地面投影点P0（R,φ,λ）处的已知观测值.式（4）即为人们常用的经去奇异性改化后的地球外部重力异常全球积分计算公式.从形式上不难看出，引入积分恒等式（3）变换后，计算式（4）不再存在积分奇异性问题，同时确保当r→R时，积分计算值Δgp收敛于球面观测量ΔgRp.Heiskanen和Moritz（1967）认为，经式（3）变换后，至少可以说式（4）核函数的奇异性被中和了.于锦海等（2001）从理论上证明了式（4）右端奇异积分项在Chauchy主值意义下的存在性，从而为实施式（4）及其微分算子的数值计算提供了理论依据.

地球外部重力异常径向偏导数计算模型可直接由式（4）求微分得到

（5）

依据式（5）可推出重力异常一阶径向偏导数Δg′rp计算模型为

（6）

=K1（r,ψ）.

（7）

如前所述，受观测数据覆盖范围的限制，在实际使用式（6）计算地球外部重力异常一阶径向偏导数时，通常需要将全球积分域划分为近区和远区处理，近区定义为以计算点为中心、ψ0为半径的球冠区域σ0，剩下的部分称为远区（σ-σ0）.一般以一定阶次（比如N阶）的重力位模型作为参考场，联合采用实测重力数据和移去恢复技术，对近区数据影响进行数值积分计算；远区效应则采用更高阶次（比如L阶）的重力位模型进行补偿（黄谟涛等，2020）.引入基于参考场的“移去-恢复”处理模式后，还需要对积分核函数作相应的改化处理，以满足积分核函数与观测重力异常信息之间的频谱匹配要求（Novák and Heck,2002；刘敏等，2016）.这里统一使用简单实用的Wong和Gore（1969）方法对积分核函数进行改化，即从原核函数中截断掉与位模型参考场相同阶次的球谐展开式.经分区处理和核函数改化后，计算式（6）从全球积分模型转换为局域积分模型，也就是第一步改化模型：

（8）

×Pn（cosψ），

（9）

（10）

（11）

（12）

（13）

需要指出的是，式（8）并不是严密的改化模型，还不能作为最终的实用化公式使用.这是因为，式（8）是由全球积分式（6）改化为局域积分得来的，式（6）又是从积分恒等式（3）变换过来的，而恒等式（3）成立的前提条件是积分域覆盖全球，由分区改化得来的计算式（8）对应的是局域积分，显然不满足积分恒等式（3）的假设条件要求.尽管在式（8）的右端已经通过超高阶位模型计算值Δg′rq（σ -σ0）顾及了远区效应的影响，但Δg′rq（σ -σ0）只代表式（8）右端积分项Δgq在远区（σ-σ0）的补偿，并未顾及另一积分项ΔgRp在远区（σ-σ0）的影响.这就意味着，当我们采用局部区域观测数据完成式（8）计算时，其计算结果必然存在一定大小的系统性模型偏差.数值试验结果表明，要想获得高精度的外部重力异常垂向梯度计算结果，必须消除该项误差的影响.

从前面的分析得知，由全球积分过渡到局域积分引起计算式（8）的模型误差可用公式表示为

（14）

（15）

（16）

式中，Δg′rp（σ-σ0）代表ΔgRp在积分远区（σ-σ0）对计算参量Δg′rp的影响.将式（7）代入式（15），并进行积分运算，可推得

+rRcosψ0（r2-5R2）]，

（17）

（18）

在式（8）右端加入模型误差修正项Δg′rp（σ -σ0），可得到计算外部重力异常一阶径向偏导数的第二步改化模型：

+Δg′rq（σ -σ0）+Δg′rp（σ -σ0）.

（19）

按照同样的思路可推得外部重力异常二阶及以上高阶偏导数相对应的改化公式，但考虑到由式（5）定义的高阶偏导数涉及复杂的观测函数连续性和核函数强奇异性问题，我们拟另文作专题研讨，这里不再展开讨论.需要指出的是，在重力异常场变化比较剧烈的区域，使用式（19）计算重力异常垂向梯度还会带来一定的模型误差.这是因为，在式（19）右端的积分项中，我们是把计算点所在的网格数据块重力异常当作常值ΔgRp处理的，该数据块对计算参量Δg′rp的影响已经在式（19）右端的第一项和最后一项中得到补偿，在积分项中不再体现该数据块的影响.但当计算点附近的重力异常场变化比较剧烈时，再将计算点所在数据块当成常值处理可能带来不可忽略的误差，必须对其作相应的补偿.假设与计算点重合的网格数据块半径为ψ00，考虑到该数据块是一个很小的区域，故可采用极坐标系（s,α）对积分核函数作平面近似处理.取

略去（h/R）2及以上高阶项影响，可将由式（7）表示的积分核函数近似表示为

（20）

此时，计算点所在数据块的积分式可写为

（21）

式中，Δg′rp0代表计算点所在数据块重力变化特征对计算参量Δg′rp的影响；s0代表数据网格大小的一半，当数据网格为1′×1′时，s0=0.5′.由式（21）得知，如果把中心数据块的重力异常当成常值ΔgRp看待，即认为在计算点所在的数据网格内处处满足Δgq=ΔgRp，则有Δg′rp0=0.当计算点附近的重力异常场变化比较剧烈时，可参照Heiskanen和Moritz（1967）的思路，将重力异常Δgq在空间计算点P的球面投影点Rp处展开为泰勒（Taylor）级数：

+…，

（22）

式中，x轴指向正北，y轴向东，x=scosα，y=ssinα.并且有

将式（22）代入式（21），同时考虑到修正项Δg′rp0本身的量值一般比较小，Taylor级数展开式（22）中的三阶及以上各项的综合影响可以忽略不计（Heiskanen and Moritz，1967），由此不难推得

（23）

假设与计算点重合的数据格网为（i,j），可按式（24）、（25）计算gxx和gyy：

（24）

（25）

将式（23）加入式（19）的右端，就得到计算外部重力异常垂向梯度的严密改化模型：

+Δg′rq（σ-σ0）+Δg′rp（σ-σ0）+Δg′rp0.

（26）

1.2 地面重力异常垂向梯度严密改化模型

相比地球外部重力异常垂向梯度，在地球重力场逼近计算实际应用中，人们更关注的是地球表面或重力观测面（近似为球面）上的垂向梯度.为此，这里专门给出球面重力异常一阶径向偏导数严密改化公式.

实际上，球面重力异常垂向梯度只是外部重力异常垂向梯度的一个特例.在式（6）和式（7）中，令r=R，可推得球面上的重力异常一阶径向偏导数Δg′Rp计算模型：

（27）

（28）

l0=2Rsin（ψ/2），

（29）

将式（28）代入式（27）得

（30）

式（30）就是重力归算应用中最常见的球面重力异常一阶径向偏导数计算模型（Heiskanen and Moritz,1967）.

显然，在实际应用中，同样需要对式（30）右端的全球积分做分区改化处理，同时需要引入位模型参考场对已知参量和待求参量进行移去恢复计算和积分核函数截断处理.基于与式（6）同样的改化流程，首先将式（30）右端的全球积分改化为近区积分和远区球谐函数展开两部分：

（31）

（32）

（33）

式中，Δg′Rq（σ-σ0）代表球面重力异常一阶径向偏导数远区效应位模型计算值；Qn（Δg′R）为相对应的Poisson积分核截断系数；其他符号意义同前.很显然，式（31）—（33）可以直接由式（8）—（11）令r=R得到.类似于式（8），经第一步改化后的式（31）同样存在由全球积分过渡到局域积分引起的模型误差，该误差大小可由式（34）确定：

（34）

l0（ψ0）=2Rsin（ψ0/2），

（35）

式中，Δg′Rp（σ-σ0）代表ΔgRp在积分远区（σ-σ0）对计算参量Δg′Rp的影响.同样可以证明，式（34）可直接由式（14）令r=R得到.

在式（31）的右端加入模型误差修正项Δg′Rp（σ -σ0），可得到计算球面重力异常一阶偏导数的第二步改化模型：

显然，计算球面重力异常一阶偏导数的严密改化公式还应当包含与式（23）相对应的计算点所在数据块的影响（Heiskanen and Moritz,1967）.在式（23）中令h=0（即r=R），可直接求得

（37）

式中，Δg′Rp0代表计算点所在数据块重力变化特征对计算参量Δg′Rp的影响.式（37）与Heiskanen和Moritz（1967）的推导结果完全一致，将其加入式（36）的右端，就得到计算球面重力异常垂向梯度的严密改化模型：

+Δg′Rq（σ-σ0）+Δg′Rp（σ-σ0）+Δg′Rp0，

（38）

式中，Δg′Rq（σ-σ0）和Δg′Rp（σ-σ0）分别代表Δgq和ΔgRp在积分远区（σ-σ0）对计算参量Δg′Rp的影响，分别由式（32）和式（34）计算；Δg′Rp0代表计算点所在数据块对计算参量Δg′Rp的影响，由式（37）计算.后面的数值计算检验，将进一步验证第一步改化公式（31）增加模型误差两个修正项Δg′Rp（σ-σ0）和Δg′Rp0的必要性及有效性.

1.3 重力异常垂向梯度向下延拓应用

如前所述，向下延拓是重力异常垂向梯度最重要的应用方向之一.除了我们常见的航空重力向地面延拓计算外（王兴涛等,2004；黄谟涛等,2018），重力异常垂向梯度最具标志性的应用场景是，将地面重力异常延拓归算到海平面或过计算点的水准面，进而用于大地测量边值问题解算（Moritz,1980；于锦海等,2001）.实际上，基于现代边值问题理论的Molodensky零阶加一阶项级数解可解释为，首先将地面重力异常向下解析延拓到海平面，用Stokes积分求得海平面上的高程异常，再将该结果向上延拓到地面（Heiskanen and Moritz,1967；李建成等，2003）.其中，地面空间重力异常Δg向海平面延拓的计算公式可表达为

（39）

式中，Δg*代表海平面上的重力异常；h代表地面点的正常高；（∂Δg/∂h）为地面重力异常的垂向梯度，通常使用前面介绍的一阶径向导数（∂Δg/∂r）来替代.美国从20世纪90年代开始，每隔3至5年就会对作为国家高程基准的大地水准面模型进行更新换代，新发布的USGG2009模型在构建过程中，就使用了式（39）作为地面重力异常的归算模型（Wang et al.,2012）.显然，与完整的向下解析延拓模型相比较（Moritz,1980），式（39）已经事先省略了二次及更高次项的影响，关于这些高次项影响的讨论已经超出本文的研究范围，故不再做更多的评述，这里主要就一阶径向导数（∂Δg/∂r）计算模型的完备性对重力异常向下延拓精度的影响进行分析和验证，具体见后面的数值计算检验环节.

2 数值计算检验与分析

2.1 数值检验使用的数据及区域

为了分析比较前述不同阶段改化模型的计算效果，本文采用超高阶位模型EGM2008作为数值计算检验的标准场（Pavlis et al.,2012），用于模拟产生球面1′×1′网格重力异常观测量“真值”（这里使用1′×1′而非5′×5′网格数据是为了减弱积分离散化误差的影响），同时产生球面及外部设定高度的重力异常垂向梯度理论“真值”.由重力位模型计算地球外部重力异常及一阶径向偏导数的公式为（Heiskanen and Moritz,1967；黄谟涛等，2005）

（40）

（41）

式中各个符号意义同前.在式（40）和（41）中令r=R，可得到计算地球表面（球面）重力异常及其一阶径向偏导数的公式.

为了体现检验结果的代表性，这里特意选取重力场变化比较剧烈的马里亚纳海沟作为试验区，具体覆盖范围为：6°×6°（φ：10°N—16°N；λ:142°E—148°E）.首先选取截断到360阶次的位模型EGM2008作为参考场，即取N=360，然后选取361～2160阶次的位模型EGM2008作为数值计算检验的标准场，即取L=2160，进而选取ri=R+hi，R=6371 km，使用EGM2008模型（361～2160阶次）分别计算标准场7个高度面上的1′×1′网格重力异常观测量“真值”Δgti及相对应的一阶径向偏导数理论“真值”Δg′ti（i=1,2，…，7）,每一个高度面对应360×360=129600个网格点数据，7个高度分别为：hi=0 km,0.1 km,0.3 km,1 km,3 km,5 km,10 km.表1列出了其中5个高度面的理论“真值”统计结果，图1和图2分别给出了对应于零高度面的重力异常及其一阶径向偏导数理论“真值”的分布态势.

表1统计结果和图1、图2显示的重力异常及一阶径向偏导数变化形态说明，尽管已经扣除掉2～360阶次频段的参考场，本试验区域标准场重力变化的激烈程度仍然非常显著，可在一定程度上代表真实地球大部分局部重力场的变化特征.

表1 由EGM2008模型（361～2160阶次）计算得到的重力异常及其一阶导数统计结果Table 1 Gravity anomalies and their first order derivatives obtained by the EGM2008 model of degree 361～2160

2.2 垂向梯度改化模型检验结果分析

为了对比分析不同改化模型的计算效果，这里首先采用零高度面上的1′×1′网格重力异常“真值”Δgt0作为观测量，同时使用前述4种地球外部重力异常垂向梯度（Δg′rp）改化模型，对前面选定的试验区对应于7个高度面上的1′×1′网格重力异常一阶偏导数进行数值计算检验和分析，其中，第1模型是指直接对式（1）求径向偏导数作为基础计算模型，并对全球积分域作了分区处理，但在实施近区计算时，扣除掉与计算点重合的那个1′×1′数据块，以避免出现奇异积分；第2模型对应于公式（8）；第3模型对应于公式（19）；第4模型对应于公式（26）.将4种改化模型的计算值分别与相对应高度面的理论“真值”Δg′ti作比较，可获得不同改化模型的精度评估信息，具体比对结果列于表2.这里积分半径统一取为ψ0=2°，为了减小积分边缘效应对评估结果的影响，表2只列出中心区2°×2°方块内的比对结果（下同）.为了定量评估由全球积分过渡到局域积分引起的模型误差影响，这里同时给出了采用式（14）计算得到的两组分别对应于ψ0=2°和ψ0=5°的误差补偿量Δg′rp（σ-σ0）统计结果，具体见表3，该补偿量是本文推出的严密改化公式中最重要的修正项.

图1 球面重力异常分布Fig.1 The gravity anomalies on the sphere

图2 球面重力异常一阶径向偏导数分布Fig.2 The first order radial partial derivatives of the gravity anomaly on the sphere

在前述试验基础上，我们进一步采用同一高度面上的1′×1′网格重力异常“真值”Δgti作为观测量，同时使用前述4种地面重力异常垂向梯度（Δg′Rp）改化模型，对相对应7个高度面上的1′×1′网格重力异常一阶偏导数进行数值计算检验和分析，其中，第1模型是指直接对式（1）求径向偏导数并令r=R作为基础计算模型，且对全球积分域作了分区处理，但在实施近区计算时，扣除掉与计算点重合的1′×1′数据块，以避免出现奇异积分；第2模型对应于公式（31）；第3模型对应于公式（36）；第4模型对应于公式（38）.将4种改化模型的计算值分别与相对应高度面的理论“真值”Δg′ti作比较，可获得不同改化模型的精度评估信息，具体比对结果列于表4.

首先，从表2互比结果可以看出，我们对地球外部重力异常垂向梯度积分模型所作的分阶段改化处理，取得了符合预期的解算效果.第1模型的误差看似主要源于直接扣除了计算点所在数据块的影响，实质上是由于该积分模型在边界面存在不连续性所致.对比表2和表1结果可以看出，第1模型在1 km以下超低空高度段的误差量值远远超过了垂向梯度自身大小，显然，这不是忽略计算点所在数据块影响所能引起的量值，而是第1模型原始计算式（直接对式（1）求径向偏导数得到）在边界面存在比较显著的类似于质面和质体位那样的数值跳跃所致，这是由地球重力位在边界面存在不连续特性所决定的（Heiskanen and Moritz,1967）.这个结果说明，重力异常垂向梯度原始计算模型在超低空高度段是失效的，只有在5 km以上计算高度才是可用的.第2模型从理论上消除了第1模型的积分奇异性和数值不连续性影响，计算精度得到显著提升，其相对检核精度（指互差均方根/垂向梯度自身）都控制在20%以内，但由于该模型的改化过程存在不可忽略的理论缺陷，在10 km以上高度，该模型的计算精度反而不及第1模型.第3模型从理论上弥补了第2模型的缺陷，使得该模型的计算精度得到进一步改善，在超低空高度段，该模型的相对计算精度不低于7%，在1 km以上高度段，相对精度优于3%.这个结果说明，我们对第2模型所作的修正和补偿处理是正确且有效的.第4模型是在第3模型基础上，增加了计算点所在数据块重力变化特征对计算参量的影响，表2结果显示，相比第3模型，第4模型计算精度在300 m以下超低空高度段又得到一定程度的提升，在零高度面，其相对计算精度从6.7%提升到4.0%，提升幅度超过40%，充分体现了该模型的改化效果.可以预见，当采用的数据网格间距加大（比如从1′×1′增大到2′×2′）且计算点周围的重力异常场变化更为剧烈时，第4模型的改化效果会更加显现.

表2 由不同外部改化模型计算得到的7个高度面重力异常一阶导数与“真值”比较（单位：mGal·km-1）Table 2 Comparisons between the first order radial partial derivatives of gravity anomalies,obtained by different modified models outside the earth,and the “true values”on 7 altitude surfaces （unit:mGal·km-1）

表3 模型误差补偿量Δg′rp（σ-σ0）计算结果统计（单位：mGal·km-1）Table 3 Statistics of the computational results of model error compensation Δg′rp（σ-σ0） （unit:mGal·km-1）

表4 由不同地面改化模型计算得到的7个高度面重力异常一阶导数与“真值”比较（单位：mGal·km-1）Table 4 Comparisons between the first order radial partial derivatives of gravity anomalies,obtained by different modified models on the surface,and the “true values”on 7 altitude surfaces （unit:mGal·km-1）

由表3计算结果可进一步看出，尽管第3模型对第2模型的补偿量均随参考场阶数N、积分半径ψ0和计算高度h的增大而减小，但当参考场阶数取N=360时，即使积分半径增大到ψ0=5°，7个高度面的误差补偿量均方根值仍然接近甚至超过垂向梯度自身大小的10%.这样的结果再次说明，对于高精度要求的地球重力场逼近计算，对重力异常垂向梯度传统积分模型进行精细改化和校正是非常必要的.

对比表4和表2计算结果可以看出，在相同的数据分辨率和精度保障条件下，利用某一高度面的重力异常观测数据计算该高度面外部的重力异常垂向梯度，其精度都要比利用本高度面观测数据计算本高度面的垂向梯度精度高.这个结果显然跟重力异常垂向梯度积分计算模型误差随计算高度升高而衰减有关.相比较而言，因第1模型在边界面存在较大的数值跳跃问题，故利用该模型和本高度面观测数据计算本高度面垂向梯度的结果偏离理论“真值”的幅度最大，完全失去了其使用价值.第2至第4模型两种方式计算得到的垂向梯度精度差异相对较小，但在条件允许情况下，仍应优先采用前一种脱离边界面的方式进行垂向梯度计算.

2.3 改化模型向下延拓应用效果分析

为了考察垂向梯度计算模型改化误差对重力异常向下延拓解算结果的影响，这里特别设计如下试验流程：

步骤一：使用地球外部6个高度面上的重力异常一阶径向偏导数理论“真值”Δg′ti，依据公式（39）将对应于6个高度面上的重力异常Δgti向下延拓到零高度面，分别将各个高度面的延拓计算值与零高度面的理论“真值”Δgt0作比较，计算互差统计结果.

步骤二：将式（39）中的垂向梯度替换为与表2统计结果相对应的外部重力异常一阶径向偏导数（Δg′rp）4种改化模型的计算结果，重复步骤一的试验.

步骤三：将式（39）中的垂向梯度替换为与表4统计结果相对应的地面重力异常一阶径向偏导数（Δg′Rp）4种改化模型的计算结果，重复步骤一的试验.

前述三步骤计算统计结果列于表5，为节省篇幅，表中只列出其中的对比互差均方根值（RMS）.

从表5统计结果可以看出，垂向梯度计算模型误差直接影响重力异常向下延拓的解算精度，对积

表5 不同改化模型计算一阶导数用于重力异常向下延拓与“真值”比较（互差均方根值）（单位：mGal）

分计算模型的修正和改化效果已经在向下延拓的解算结果中得到充分体现.不难看出，使用球面外部第3和第4模型计算得到的垂向梯度（Δg′rp）实施重力异常向下延拓解算的效果，已经完全等同于使用垂向梯度理论“真值”（Δg′ti）获得的解算效果，两个模型计算值与使用“真值”计算结果的差异不超过0.1 mGal；而使用第1模型时，两者的差异超过10 mGal；使用第2模型时，两者的差异也超过1 mGal,这些结果再次验证了严密改化模型的有效性.需要指出的是，重力异常向下延拓的解算精度除了与垂向梯度计算精度水平密切相关外，还取决于向下延拓模型自身的完备性、延拓高度大小及计算区域重力场变化的剧烈程度.就本试验而言，由表5可以看出，将第1模型排除在外，使用只顾及到一阶项的向下延拓模型即公式（39）进行重力异常归算，要想得到优于1 mGal的解算精度，必须将向下延拓高度控制在1 km以内，否则，需要将向下延拓模型拓展到更高阶次（黄谟涛等,2018）.这方面的内容已经超出本文的研究范围，这里不再做深入讨论.

3 结论

将全球积分模型改化为局域模型是实现重力异常垂向梯度精密计算的前提条件.本文分析研究了重力异常垂向梯度全球积分计算模型的技术特点和适用条件，指出了开展积分模型精密改化的必要性和可行性.针对全球积分模型向局域积分转换中遇到的积分奇异性和不连续性问题，综合采用积分恒等式变换和移去恢复换算技术，同时依据实测数据保障条件，分别推出了计算地球外部及地面重力异常垂向梯度全球积分模型的分步改化公式，提出了补偿传统改化模型理论缺陷的修正公式.采用超高阶地球位模型EGM2008作为模型比对标准重力异常场，同时选择在重力异常场变化比较剧烈的马里亚纳海沟区块开展数值计算试验，分别对本文推出的重力异常垂向梯度两类8种分步改化模型的计算精度及向下延拓应用效果进行了检核分析和评估.试验结果表明，采用最终的严密改化模型不仅可以有效消除原计算模型固有的积分奇异性和数值跳跃问题，又可显著提高超低空重力异常垂向梯度的计算精度和稳定性，有效提升重力异常向下延拓的解算精度水平.因此，新的严密改化模型具有较高的推广应用价值，可用于地球表面及外部重力场的高精度逼近计算.