本文内容来源于《测绘学报》2023年第10期（审图号GS京(2023)1931号）

基于改进径向导数的解析向下延拓优化算法

马健1,2,3, 翟振和1,2, 冯长强1,2, 管斌1,2, 王云鹏1,2, 李端1,2

1. 西安测绘研究所, 陕西 西安 710054;

2. 地理信息工程国家重点实验室, 陕西 西安 710054;

3. 深圳大学建筑与城市规划学院, 广东 深圳 518060

基金项目：国家自然科学基金(42204002);中国博士后科学基金(2022M712162)

摘要：解析延拓算法是一种重要的位场延拓方法。与Poisson延拓算法相比, 解析延拓算法具有算法简单、运算速度快等特点, 在基于第二或第三大地边值理论构建(似)大地水准面或垂线偏差模型的实践中具有重要的应用价值。径向导数的计算是解析延拓算法的关键, 而传统解析延拓算法中径向导数的计算存在一定程度的近似误差。本文提出了一种径向导数的改进方法, 在此基础上实现了解析延拓算法的优化。本文利用球谐方法从理论上证明了改进径向导数比传统径向导数更接近理论严密的径向导数, 并通过绘制经度剖面和纬度剖面图直观地印证了这一结论。在我国中部地形起伏较大的区域进行的试验表明, 改进径向导数精度较传统径向导数精度提高了32.45%。在1、2、3、4 mGal的误差条件下, 改进径向导数的解析延拓算法较传统解析延拓算法向下延拓精度分别提高了29.04%、19.48%、10.12%、2.65%。本文通过理论分析和算法试验证明了基于改进径向导数的解析向下延拓优化算法的有效性, 该方法在大地边值解算中具有一定的应用价值。

关键词：解析延拓算法 传统径向导数 改进径向导数 向下延拓 球谐方法

引文格式：马健, 翟振和, 冯长强, 等. 基于改进径向导数的解析向下延拓优化算法[J]. 测绘学报，2023，52(10)：1650-1660. DOI: 10.11947/j.AGCS.2023.20220625

MA Jian, ZHAI Zhenhe, FENG Changqiang, et al. Optimization algorithm for analytical downward continuation based on improved radial derivative[J]. Acta Geodaetica et Cartographica Sinica, 2023, 52(10): 1650-1660. DOI: 10.11947/j.AGCS.2023.20220625

阅读全文：引 言

Poisson延拓和解析延拓是两种重要的位场延拓算法[1]，延拓算法在物理大地测量领域具有重要的应用价值。一方面，航空重力数据[2-3]需通过向下延拓算法转化为地面重力数据；另一方面，基于大地边值理论构建(似)大地水准面或垂线偏差模型时需将地面重力数据延拓到边界面上。

过去，众多国内外学者对Poisson延拓算法的优化方法进行了深入的研究，取得了丰硕的研究成果。进行Poisson延拓时需对Poisson积分离散化，文献[4]通过对几种Poisson离散化模式进行试验得出单平均模式的离散化误差最小的结论。针对传统Poisson算法存在的缺陷，文献[5]和文献[6]分别采用顾及远区影响的Poisson延拓算法、重力异常向上严密改化模型实现了Poisson算法的优化，在一定程度上提高了Poisson延拓算法的精度。病态问题的存在严重影响了Poisson向下延拓精度，文献[7—11]通过优化Poisson向下延拓的正则化算法以减弱病态性的影响。为规避传统Poisson积分向下延拓中的不适定问题，文献[12—14]提出了顾及地形效应的延拓方案，通过试验表明位模型加地形改正的延拓模型具有良好的计算稳定性。在不同延拓算法延拓效果方面，文献[15]通过试验表明直接代表法和正则化法的延拓算法具有较高的延拓精度和可靠性，文献[16]表明改进的Poisson迭代法精度略优于改进最小二乘法，略低于Tikhonov正则化延拓法。

解析延拓算法可分为一阶和高阶解析延拓算法，目前普遍应用的算法为一阶解析延拓算法。在高阶解析延拓径向导数算法推导方面，文献[17]对高阶解析延拓径向导数推导方法进行研究，给出了直接求导法与递推法、迭代求导法推得的解析延拓算式存在一定的差异的结论。文献[18]基于直接推导法首次推导了任意调和函数和重力异常一至九阶的径向导数具体算式，其给出的高阶径向导数算式可直接应用于高阶解析延拓算法。在高阶解析延拓径向导数算法优化方面，文献[19]推导了一组与带限核函数相对应的重力异常高阶径向导数带限计算公式，通过试验表明使用带限径向导数的高阶解析延拓算法具有良好的可靠性和有效性。文献[20]通过对解析延拓高阶径向导数进行频域Tikhonov正则化，提高了高阶解析延拓算法的稳定性。文献[21]针对传统重力异常垂向梯度的理论缺陷，提出了重力异常垂向梯度的修正公式，垂向梯度的实质为一阶径向导数，将垂向梯度修正算法应用于解析延拓算法即可实现一阶解析延拓算法的优化。

近年来，有学者开展了向下延拓的谱方法研究，如文献[22]提出了基于谱融合的延拓算法，该算法的优势可实现稳定的向下延拓。文献[23]针对卫星和航空重力数据向下延拓到不规则的地面的应用场景，提出了不同类型重力场元的谱域向下延拓方法。

根据大地边值理论[24-28]构建区域(似)大地水准面或垂向偏差模型时，需要将地面格网重力数据延拓到边界面上。与航空重力数据相比，地面格网重力数据通常区域范围更大、分辨率更高，且其误差受地理条件影响较大。由于Poisson向下延拓算法计算复杂，耗时长，对重力数据噪声敏感导致其可能存在病态问题等，因而在边值解算中的应用效果并不理想。与Poisson延拓算法相比，解析延拓算法具有算法简单、计算速度快、无须处理病态问题等优点，因而在基于边值理论构建区域(似)大地水准面或垂线偏差模型时解析延拓算法具有Poisson算法无法比拟的优势。径向导数的计算是解析延拓算法的关键，而传统解析延拓算法中径向导数的计算存在一定程度的近似误差。为了提高解析延拓算法的精度，本文对解析延拓算法展开深入研究，提出一种径向导数的改进方法，并通过理论分析和算法试验对基于改进径向导数的解析向下延拓算法的有效性进行了说明。

1 Poisson算法与解析延拓算法Poisson算法和解析延拓算法是位场延拓的两大理论工具。位场延拓的数据对象通常为重力异常或重力扰动，实际上两种重力数据的延拓算法是相同的。以重力扰动作为延拓数据为例，Poisson算法的实用公式为[1] (1)式中，Ω0代表全球域；dΩ′为球面微分单元；下标P表示重力数据测量面计算点；下标Q表示测量面流动积分点；上标*表示边界面上的值；rP为计算点的地心向径；为计算点与边界面积分点之间的距离(ψ为计算点与积分点之间的地心角距)；δgP为计算点处的重力扰动；δgQ*为边界面上流动积分点的重力扰动；δgP*为计算点在边界面上投影点的重力扰动。解析延拓算法的数学本质是泰勒级数展开式，其计算公式为[28] (2)式中，∂kδgP*/∂rPk为δgP*沿计算点地心向径rP的k阶径向导数值；hP为计算点高程。实际应用解析延拓算法进行向下延拓时通常采用一阶解析延拓算法，即 (3)式中，∂δgP*/∂rP为δgP*的一阶径向导数，又称垂向梯度。解析延拓算法中径向导数算式是在对Poisson积分公式求导的基础上得到的，为 (4)式中, 为边界面上计算点和流动积分点的距离，当地球平均半径R为半径的球面为边界面时，。当解析延拓算法应用于向下延拓时，式(4)右端积分中边界面重力扰动值为待求量，因而无法直接进行此积分计算。为了能够计算重力扰动径向导数，实际计算时通常将测量面重力扰动值作为边界面重力扰动的近似值，即[1] (5)式中，δgQ为测量面上流动积分点处的重力扰动。理论上，式(5)右端积分中流动点和计算点应在同一延拓球面上。由于通常相邻网格的地形相差不大，且式(5)中积分核随着计算点与积分点距离增大迅速减小，因此虽然对地面重力数据进行向下延拓时计算点和流动积分点不在同一球面上，但造成的误差非常小。将式(5)代入式(3)即可求得边界面上的重力扰动值。2 解析延拓算法的径向导数的改进方法径向导数的实质是重力梯度，重力梯度的计算是解析延拓算法的关键。理论严密的径向导数值应为测量面与边界面径向导数值(即重力梯度)的平均值。在传统的解析延拓算法中，径向导数值是通过测量面重力数据代入梯度算式中得到的，因而其值实际为测量面重力梯度值。由于重力扰动梯度并非固定值，其绝对量值通常随计算点高度的增大而衰减，因而传统解析延拓算法的径向导数的计算存在一定的误差。下文利用球谐方法进行理论的说明，并给出改进方法。重力扰动用球谐方法可表示为[1] (6)式中，GM为地心引力常数；a代表参考椭球的长半轴；球谐项Yn的表达式为 (7)式中，Cnm*、Snm为n阶m次的模型扰动位球谐系数，Cnm*由重力场模型球谐系数中的C20、C40、C60、C80、C10, 0项分别除去正常引力位中相应的2、4、6、8、10阶带谐系数得到的；φ、λ分别为计算点的球心纬度和经度；Pnm是完全正常化的伴随Legendre函数。传统解析延拓算法中，重力扰动梯度用球谐方法可表示为 (8)假设延拓面是半径为R的球面，则延拓后的边界面重力扰动真值为 (9)假设理论严密的径向导数值为，则有 (10)将式(6)、式(9)代入式(10)，可得理论严密的径向导数的公式为 (11)令理论严密的径向导数对应的地心径向为，根据式(8)，有 (12)根据式(11)、式(12)，可得 (13)传统解析延拓算法中径向导数是地心向径r处，而非理论严密的处重力值计算的，因而传统算法存在一定的近似误差。为此，本文提出了一种径向导数的改进算法，其公式为 (14)式中，δgQ、δgP分别为测量面流动积分点和计算点沿径向方向在高度为h/2处的近似重力扰动值，其计算公式为 (15)将改进径向导数替换传统解析算法中的径向导数，即可得到基于改进径向导数的解析向下延拓优化算法公式为 (16)下面用球谐方法证明改进径向导数的合理性和有效性。δg用球谐方法可表示为 (17)令改进径向导数对应的地心径向为r，根据式(6)，有 (18)根据式(17)、式(18)，可得 (19)对比式(13)、式(19)可以看出，解析延拓算法中改进径向导数对应的地心向径r与传统径向导数对应的地心向径r相比，更接近理论严密的径向导数对应的地心向径，即，改进径向导数值比传统径向导数值更接近理论严密的径向导数值。因此，基于改进径向导数的解析延拓算法理论上比传统解析延拓算法延拓精度更高。重力异常与重力扰动具有相同的物理特性，采用相同的推导方法可得重力异常数据的改进径向导数算式为 (20)式中，ΔgQ、ΔgP分别为流动积分点和计算点在改进地心向径处的重力异常，其计算公式分别为 (21)式中，ΔgQ、ΔgP分别为测量面流动积分点和计算点重力异常；∂ΔgQ/∂rQ、∂ΔgP/∂rP由传统梯度算式计算。面向重力异常数据的基于改进径向导数的解析向下延拓优化算法为 (22)式中，ΔgP*为边界面上计算点位置的重力异常延拓值。以上是对解析延拓向下算法的优化方法的推导和说明。对于向上延拓的应用，由于Poisson积分本身算法稳定、计算快速方便、延拓精度高，因而本文未对解析向上延拓算法进行优化方法研究。3 试验与分析3.1 径向导数算法分析为了直观反映重力扰动梯度随高度的变化规律，选取A(102°E、30°N)和B(102.3°E、30.5°N)两点，利用2190阶EIGEN-6C4重力场模型分别计算不同高度下A、B两点的重力扰动梯度，如图 1所示。图 1 重力扰动梯度随高度变化Fig. 1 Variation of gradients of gravity disturbance with heights图选项

图 1中，A点重力扰动梯度值均为负，B点重力扰动梯度值均为正，说明重力扰动梯度值可正可负，这与重力扰动取值有正有负的现象是一致的。通常情况下，随着高度的增加，重力扰动梯度的绝对值呈现衰减的趋势。为了反映解析延拓算法的传统径向导数、改进径向导数与理论严密的径向导数的差异，选取经度为102.3°E、纬度范围为29.2°N—30.8°N的经度剖面，以及纬度为30.0°N，经度范围为101.5°E—103.2°E的纬度剖面，利用2190阶EIGEN-6C4重力场模型分别通过不同算法计算径向导数值。图 2为两剖面的地形高度。图 2 地形高度Fig. 2 Topographical heights图选项

延拓高度越高，传统径向导数与改进径向导数的差异越明显。由图 2可以看出，本文选择的经度剖面和纬度剖面的地形起伏较大，因而有利于明显地体现两种径向导数算法的差异。图 3为两剖面传统径向导数、改进径向导数及理论径向导数。图 3 径向导数Fig. 3 Radical derivatives图选项

图 3中，传统径向导数、改进径向导数与理论严密的径向导数的变化趋势是一致的。对比图 2、图 3可知，径向导数与地形高度存在一定的相关性。由于陆地上地形高通常为正值，而重力梯度取值有正有负，因而二者之间的相关关系并非简单的线性相关关系。由图 3可知，改进径向导数比传统径向导数更接近理论上严密的径向导数，因而理论上精度更高。以理论严密的径向导数为真值，图 4反映了传统和改进径向导数的误差。图 4 径向导数误差Fig. 4 Errors of radical derivatives图选项

由图 4可知，改进径向导数普遍比传统径向导数误差更小。正因如此，基于改进径向导数的解析延拓算法的向下延拓解精度理论上应优于传统解析延拓算法。

3.2 地面模型重力延拓试验

本文选择我国中部101°E—118°E，24°N—39°N(17°×15°)范围的试验区进行向下延拓试验。采用的地形数据为SRTM3数据，SRTM3数据是以EGM96大地水准面为基准面的正高数据，因此通过将SRTM3数据与EGM96模型大地水准面高相加得到大地高数据。试验区分辨率为2′×2′，大地高最低为-7.28 m，最高为6 112.42 m，平均高程为1 019.34 m，如图 5所示。本文试验区范围较大，地形复杂，具有不同的地理特征，能够较好地反映不同地形高度下解析延拓算法的延拓效果。图 5 试验区地形Fig. 5 Topographical map of the test area图选项

试验采用的重力数据是利用2190阶的EIGEN-6C4重力场模型计算的地面和参考椭球面上2′×2′格网的重力扰动，其分布如图 6所示。表 1为地面重力扰动、参考椭球面重力扰动及二者差异的统计信息。图 6 2190阶模型重力扰动Fig. 6 2190-degree model gravity disturbances图选项

表 1 重力扰动统计Tab. 1 Statistics of gravity disturbances mGal

表选项 对比图 5、图 6可知，两种边界面的重力扰动的起伏与地形具有一定的相关性。由于本文试验区覆盖最高高度可达6000多米的高山区，参考椭球面重力扰动与地面重力扰动的差异明显，最大达280.19 mGal。由于参考椭球面与地面重力扰动的差异可作为向下延拓理论上的真值，因此本文试验区能够较好地反映延拓算法的延拓精度。解析延拓算法中径向导数的计算理论上需全球积分，由于远区积分值极小(正因如此，本文试验采用参考椭球面而非标准球面作为延拓面引起的误差非常小)，实际计算时通常只需在近区范围进行积分即可，本文选取的积分半径为1°。为了提高数据有效利用率，计算径向导数时，边缘网格的积分域未包含在试验区内的部分可补充低阶模型值或简单地设置为0。统计时舍弃边缘0.5°范围的延拓值以较大程度削弱边缘效应的影响。为了说明地面和延拓面重力扰动梯度量值的具体差异，利用2190阶的EIGEN-6C4重力场模型分别计算两种界面的重力扰动梯度值，如图 7所示。表 2为地面重力扰动梯度、参考椭球面重力扰动梯度及二者差异的统计信息。图 7 2190阶模型重力扰动梯度Fig. 7 2190-degree model gradients of gravity disturbances图选项

表 2 重力扰动梯度统计Tab. 2 Statistics of gradients of gravity disturbances mGal/m

表选项

根据图 7可看出，重力扰动梯度变化与地形起伏的趋势一致，根据前文分析，二者之间并非线性相关关系。结合图 7、表 2可看出，地面重力扰动梯度变化幅值明显小于参考椭球面重力扰动梯度变化幅值，这与图 1的结论是一致的。由于重力扰动梯度值在地面和参考椭球面有较大的差异，因此传统解析延拓算法将地面重力扰动梯度值作为径向导数值，理论上存在不可避免的误差。表 3统计了试验区传统与改进径向导数误差。

表 3 径向导数误差Tab. 3 Errors of radical derivatives mGal/m

表选项

由表 3可看出，改进径向导数精度(均方根)为1.27×10-3 mGal/m，较传统径向导数精度(1.88×10-3 mGal/m)提高了32.45%，该试验定量地体现了改进径向导数的改善效果。为了验证传统解析延拓算法、改进解析延拓算法及Poisson延拓算法对不同误差的重力数据的向下延拓效果，生成均值为0、均方差为1且服从标准均匀分布的随机数，将其作为1 mGal的随机噪声，在其基础上分别乘以2、3、4可获得2、3、4 mGal的随机噪声。将含有不同噪声的地面重力扰动分别通过各延拓算法向下延拓到参考椭球面，表 4为不同延拓算法的计算时长。

表 4 不同向下延拓算法的向下延拓计算时长Tab. 4 Calculation time of different algorithms of downward continuation min

表选项

表 4明显地反映出解析延拓算法的优势，即解析延拓算法的计算效率高，而Poisson向下延拓计算效率低。本文试验中Poisson向下延拓采用迭代法计算，重力数据误差越大，迭代达到阈值需要的时间越长。而传统和改进解析延拓算法的计算效率与重力数据误差的无关。根据表 4，改进径向导数的解析延拓法的计算时间约为传统延拓算法的2倍，但比Poisson向下延拓算法效率高15倍以上。表 5统计了各延拓算法的向下延拓误差。

表 5 含有标准均匀分布误差的重力数据向下延拓误差Tab. 5 Downward continuation errors of gravity data including standard uniform distribution errors mGal

表选项

由表 5可以看出，在1~4 mGal的误差条件下，改进径向导数的解析延拓算法的向下延拓精度均优于传统解析延拓算法，精度分别提高了29.04%、19.48%、10.12%、2.65%。向下延拓过程在放大重力数据信号的同时也放大了重力数据包含的误差，传统解析延拓算法使用量值较小的地面径向导数因而放大作用小，优化的解析延拓算法使用改进径向导数因而放大作用较大。因此，随着重力数据误差的增大，改善的效果呈现降低的趋势。表 5中，当重力数据包含1 mGal误差时，Poisson算法的延拓精度最高。但随着重力数据误差增大，Poisson向下延拓精度迅速下降，且低于改进的解析延拓算法，由此可看出，Poisson算法的放大数据误差的特性与解析延拓算法相比更为显著。需要说明的是，由于在2′×2′分辨率的本试验区，Poisson向下延拓的条件数均不超过1000，当重力数据含有1~4 mGal的误差时，均能得到稳定的向下延拓解，因此本文进行Poisson向下延拓时并未采用正则化等病态问题处理方法。为了说明不同的加噪方式对向下延拓精度的影响，生成均值为0、均方差为1且服从正态分布的随机数，将其作为1 mGal的随机噪声，在其基础上分别乘以2、3、4可获得2、3、4 mGal的随机噪声。将含有不同噪声的地面重力扰动分别通过各延拓算法向下延拓到参考椭球面，表 6统计了各延拓算法的向下延拓误差。

表 6 含有正态分布误差的重力数据向下延拓误差Tab. 6 Downward continuation errors of gravity data including normal distribution errors mGal

表选项

对比表 5、表 6可看出，对于含有两种误差分布的相同误差量级的重力数据，采用同一延拓算法的向下延拓精度基本一致，二者误差差异最大仅为0.03 mGal。表 5、表 6的试验结果说明不同的加噪处理方式对向下延拓精度的影响非常小，可忽略不计。根据表 5、表 6的试验结果可知，当对山区地面重力数据进行向下延拓时，由于此时格网重力数据精度通常不高，适宜采用改进径向导数的解析延拓优化算法。上文延拓试验采用的重力扰动为2190阶模型重力扰动，为了验证移去-恢复在向下延拓算法中的应用效果，本文最后进行了基于移去-恢复的向下延拓试验。试验过程为：将地面重力扰动除去360阶的低阶模型重力扰动，然后对剩余高频重力扰动进行向下延拓，最后将得到的延拓值与360阶参考椭球面模型重力扰动相加即为最终的延拓值。表 7统计了基于移去恢复技术的向下延拓误差，使用的噪声为表 5使用的标准均匀分布随机噪声。

表 7 基于移去-恢复的向下延拓误差Tab. 7 Downward continuation errors based on remove-compute-recovery technology mGal

表选项

对比表 5、表 7可以看出，采用移去-恢复可使不同延拓算法的延拓精度均得到改善，但改善的程度非常小(不足0.2 mGal)。这是因为向下延拓算法的物理特性是放大高频信号和误差，而移去-恢复的实质是避免低频信息的传播误差，因此移去-恢复对延拓算法的改善程度较小。4 结论解析向下延拓算法是一种重要的延拓方法，与Poisson算法相比具有算法简单、计算速度快等优点，在基于Stokes-Helmert/Hotine-Helmert理论构建区域(似)大地水准面等实践中具有重要的应用价值。径向导数的计算是解析延拓算法的关键，传统的径向导数算法存在一定的近似误差，本文在对其进行深入分析的基础上给出了改进方法，具有较高的应用价值。(1) 本文利用球谐方法证明了传统解析延拓算法径向导数算法存在一定的近似误差。提出了解析延拓算法径向导数的改进方法，并用球谐方法从理论上证明了改进径向导数比传统径向导数更接近理论严密的径向导数，通过经度剖面和纬度剖面图直观地反映了这一观点。(2) 本文选择我国中部较大的区域进行试验，改进径向导数精度较传统径向导数精度提高了32.45%。在1、2、3、4 mGal的误差条件下，改进径向导数的解析延拓算法较传统解析延拓算法向下延拓精度分别改善29.04%、19.48%、10.12%、2.65%，从而说明了本文提出的基于改进径向导数的解析延拓优化算法的有效性。作者简介

第一作者简介：马健(1988-), 女, 博士, 助理研究员, 研究方向为物理大地测量。E-mail: majian_geodesy@163.com

初审：张艳玲复审：宋启凡
终审：金 君

资讯