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

第二类勒让德函数的重规格化及其优化

张捍卫,1,2, 杨永勤,2, 李晓玲2, 张华2

1.山东理工大学建筑工程与信息工程学院，山东 淄博 255000

2.河南理工大学测绘与国土信息工程学院，河南 焦作 454003

基金项目

国家自然科学基金(42074002；41931075)

作者简介

第一作者：张捍卫（1967—），男，博士，教授，博士生导师，主要从事大地测量学的教学与研究。E-mail：zhanwei800@163.com通讯作者: 杨永勤 E-mail：212104010025@home.hpu.edu.cn

摘要

椭球谐函数级数展开是地球重力场椭球谐建模的基础。然而，处理椭球谐函数级数的主要困难在于计算第二类勒让德函数。而Jekeli的重规格化方法简化了这一计算过程。本文在Jekeli重规格化的基础上，详细推导了两种基于高斯超几何函数变换的优化递归方法，同时利用这两种优化递归的方法计算了第二类勒让德函数，并将其展开到二阶导数。通过数值计算，证明了优化递归方法可以有效加快收敛，缩短计算时间，并且适用阶数更高，这使得椭球谐函数级数在实际应用中更加方便、可行。

关键词

第二类勒让德函数; 缔合勒让德微分方程; 重规格化; 高斯超几何函数

本文引用格式

张捍卫, 杨永勤, 李晓玲, 张华. 第二类勒让德函数的重规格化及其优化[J]. 测绘学报, 2024, 53(7): 1298-1307 doi:10.11947/j.AGCS.2024.20230283

ZHANG Hanwei, YANG Yongqin, LI Xiaoling, ZHANG Hua. Renormalization and its optimization of the Legendre function of the second kind[J]. Acta Geodaetica et Cartographica Sinica, 2024, 53(7): 1298-1307 doi:10.11947/j.AGCS.2024.20230283

阅读全文

球谐函数级数和椭球谐函数级数展开是科学技术领域强大的数学工具[1]。在地球重力场建模中，作为拉普拉斯方程的解，通常使用球谐函数级数来近似重力场。然而，实际测量数据近似分布于椭球表面，因此椭球谐函数级数更适合展开地球重力场[2]。但是，由于涉及复数参数[3]，椭球谐函数级数的计算较为困难。

目前，计算（截断）椭球谐函数级数主要有两种方法。一种方法是在e2（或E2）方面展开Qn，m，但随着最大阶数的增加，需要更多的微扰项[4]；另一种方法是根据其倒数（使用超几何级数）展开Qn，m，这实际上已经解决了Qn，m的计算问题。关于超几何函数的递归，已有大量的研究可供参考。文献[5]将第二类勒让德函数的复数表达式Qn，m重规格化为实数表达式Sn，m。文献[6]提出并验证了超几何函数Sn，m的递归，并研究了Jekeli理论的实际应用。文献[7]研究了椭球谐函数级数高阶计算时数值的不稳定性。文献[8]基于分布在旋转椭球体上的重力数据，使用椭球谐函数级数方法构造了边值问题的解。文献[9]进一步发展了这项研究，并更明确地表示椭圆体校正。文献[10]在解决边值问题的椭圆调和方法中引入格林函数，将复杂的非球形边界问题转化为简单边界上的问题。文献[11]为地球磁场的位及其一阶导数构造了椭球调和级数，其中计算超几何函数Sn，m的递推关系取自文献[6]。文献[12]推导了用于计算超几何函数及其导数的替代递推关系。文献[13]计算了地球引力位的椭球谐函数级数及其一阶和二阶导数，并应用不同形式的高斯超几何函数实现传统的Jekeli函数的优化计算。国内，文献[14]采用Laurent级数形式计算第二类勒让德函数，简化了第二类勒让德函数的计算。文献[15]推导并计算了第二类勒让德函数的一阶和二阶导数，并表明修正递归后计算效率更高。

当前的研究进一步修正并优化了第二类勒让德函数的递归方法，其中对高斯超几何函数加减项的处理成为优化递归方式的重要途径。本文在Jekeli重规格化的基础上，推导了两种基于高斯超几何函数变换的优化递归方法，通过数值计算，确定了不同优化方式的适用阶数范围，实现了对第二类勒让德函数的高效计算。

1 第二类勒让德函数的重规格化1.1 椭球谐函数级数

首先引入共焦椭球，图1给出了用球坐标（r，θ，λ）和椭球坐标（u，ϑ，λ）表示的点P的位置。坐标系的原点与地球质心重合；其中半径v被定义为。图1中θ是极距，在椭球坐标中，ϑ表示归化余纬度，λ是地心经度，r是地心半径。图1中蓝色实线b-椭球是参考椭球，表示为∑R，b，b是其短半轴，R是其对应的边界球半径；红色实线u-椭球是共焦椭球，表示为∑R，u，u是通过点P的共焦椭球面的短半轴，Ru是其对应的边界球半径。线性偏心率为E，两个椭球具有相同的线性偏心率。

图1

图1 共焦椭球的定义

Fig.1 Definition of confocal ellipsoid

对于扁椭球域，由椭球坐标（u，ϑ，λ）可以得到引力位V的扁椭球球向近似[16]

(1)

式中，i2=-1，引力位V随高度的增加而减小，由第二类缔合勒让德函数Qn，m（i u/E）/Qn，m（i b/E）的比值控制。其中，GM表示地心引力常数，R为参考椭球体的半长轴（或边界球体半径）。是第一类n阶m次的完全归一化缔合勒让德函数。、是完全规格化椭球谐系数。

在实际应用中，通过使用椭球谐系数的有限值nmax来代替对n的无穷大求和，最大值nmax决定了引力位V所表示的重力场频谱和空间分辨率。

1.2 第二类勒让德函数的定义

直接计算第二类勒让德函数Qn,m较为困难，因为它涉及虚数参数，并且可能会出现数值不足或溢出现象。而球谐系数与椭球谐系数之间的转换及重力位的椭球谐函数级数展开都需要用到第二类勒让德函数，目前，已经提出了许多用于高效计算第二类勒让德函数的策略，如递归[16-17]、递归的进一步简化[7,18]、超几何函数形式[7]，或以重规格化的方式重新定义函数Qn，m[5-6,19-22]。

第二类勒让德函数有几种定义存在，在这里使用以下定义[3]

(2)

式中，n为阶数；m为次数；z=i（u/E），当|z|>1时定义Qn，m（z）；2F1是高斯超几何函数，定义为[23]

(3)

式中，α，β，γ∈R且γ∉Z≤0，δ为变量，当|δ|<1时超几何函数2F1收敛[24]，k表示超几何函数的整数指数，（x）k为Pochhammer符号，定义如下

(4)

式中，Γ表示Gamma函数。本文中的超几何函数须求和为最大kmax，以保证所需的精度。式（3）收敛的速度取决于项α、β、γ和δ的相对大小。增大γ或者减小α、β、δ可以加速2F1的收敛。

1.3 第二类勒让德函数的Jekeli重规格化方法

处理椭球谐函数级数的一种不同的方法是将第二类勒让德函数式（2）重规格化，它最初是由Jekeli引入[5]，用于椭球谐系数与球谐系数之间的变换。其方法是对函数Qn，m（z）添加乘数即式（5），）椭球谐展开式（1）中的分数保持不变，通过添加乘数得到实函数为[5]

(5)

这是Jekeli重规格化函数或Jekeli函数，其中m=0时，δm=1；m≠0时，δm=1/2。重规格化函数与非重规格化函数之间的关系由式（6）给出

(6)

式中，Qn，m（z）是式（2）定义的函数，参见文献[5]。

用式（5）和式（6）化简式（1），可以将引力位表达式变换为

(7)

此外，结合式（2）和式（5）可以得到超几何函数的函数表达式为[5]

(8)

式中，递归为[5-6]

(9)

式中

Jekeli函数对u的一阶偏导数为

(10)

式中

对u的二阶偏导数，以同样的方式展开递归式，得到

(11)

式中

利用上述、、递归计算方法，可以使用椭球谐函数级数计算重力梯度。

2 Jekeli函数递归计算的精度评定方法2.1 递归的相对精度

在递归计算的过程中，需要设定相对精度以控制求和项数。在、、递归计算中，即是要确定高斯超几何函数的求和精度。关于项数k的无限超几何函数进行求和，在实际计算中只需要对收敛项k进行有限项求和就可以得到所需精度。本文使用的精度是16位有效数字即ε=10-16，其方法是对函数Sn，m的0、1、2阶中的项αn，m，k、（∂αn，m，k/∂u）、（∂2αn，m，k/∂u2）求和，直到达到精度ε=10-16。

计算使用的是GRS80椭球参数，其中u=b=6 356 752.314 m。验证了n=m=2500的递归。图2中填充圆表示独立项αn，m，k及其对u的导数，实线表示它们对收敛项k的部分和。当k=60时，∑kαn，m，k、∑k（∂αn，m，k/∂u）、∑k（∂2αn，m，k/∂u2）开始变为常数，但如果想要达到设定精度仍需要继续求和，当独立项变得小于其常数值的大约16个数量级可以停止求和，因为这时达到了16位有效数字的相对精度，由图2可以看出在当前情况下大约需要80项。

图2

图2 当n=m=2500时，超几何级数的收敛性

Fig.2 Convergence of super geometric series when n=m=2500

2.2 递归的绝对精度

除了要满足在计算前设定的相对精度外，由于Jekeli重规格化的实质是将第二类缔合勒让德函数和依赖于n阶和m次的因子相乘，所以计算得到的值、、之间还要满足由缔合勒让德微分方程所决定的绝对精度，即函数、、之间必须满足以下缔合勒让德微分方程（the associated Legendre differential equation，LDE）

(12)

上述计算得到的|LDE|的值必须尽可能为零，有效数字的相对精度与绝对精度之间的差异是由舍入误差引起的。

3 基于超几何函数变换的优化递归3.1 使用同时减小α、β项的优化递归方式计算Jekeli函数

为了减少舍入误差对计算结果的影响，可以通过加快递归函数中超几何函数的收敛速度来实现，即减少求和项的数量，这可以通过对超几何函数进行适当的变换来实现。为此进一步优化式（8），优化递归后收敛速度更快[13]。通过以下等式进行优化[23]

(13)

结合式（8）和式（13）得到

(14)

式中，上标“*”表示使用同时减小α、β项的优化递归方式计算得到的Jekeli函数。

当n接近m（扇形）时，α、β分别减小到极限1和1/2，而γ、δ不是m的函数。因此，对于特定的输入组合和所用的数字精度，式（14）收敛速度比式（8）更快。通过对式（14）的计算，得到一组优化递归式[13]

(15)

式中

Jekeli函数对u的一阶偏导数为

(16)

式中

对u的二阶偏导数，以同样的方式展开递归式，得到

(17)

式中

3.2 使用只减小β项的优化递归方式计算Jekeli函数

通过应用高斯超几何函数的线性变换性质，对式（3）进行以下变换

(18)

以式（18）为基础进一步优化式（8），用于高阶次计算[25]。进一步优化后，计算阶次可达到10 000阶，优化后的Jekeli函数为

(19)

式中，上标“**”表示使用只减小β项的优化递归方式计算得到的Jekeli函数。通过对式（19）的计算得到一组优化递归式为[25]

(20)

式中

Jekeli函数对u的一阶偏导数为

(21)

式中

对u的二阶导数，以同样的方式展开递归式，得到

(22)

式中

4 数值计算4.1 优化前后的项数对比

图3中的统计值表示在满足相对精度ε=10-16的条件下，优化递归前后对Jekeli函数表达式中超几何函数求和时所需要的项数。当n=m=2500时，与原始递归相比（第一行）使用同时减小α、β项和使用只减小β项的优化递归方式进行计算所需的项数如图3中的第2行、第3行所示。第一行（原始递归）图3（a）表示计算时所需的收敛项K0；图3（b）表示计算时所需的收敛项K1；图3（c）表示计算时所需的收敛项K2。可以看出，在原始递归方式下，对于每一个递归式，接近最大阶数nmax=2500时最多需要80项；而使用同时减小α、β项的优化递归方式（第二行）接近最大阶数nmax=2500时所需的最大项数为30；同样，使用只减小β项的优化递归方式（第三行）接近最大阶数nmax=2500时所需的最大项数也为30。此外，通过图3可以明显看出，使用这两种优化递归方式计算相同阶数的Jekeli函数时，所需的项数约为原始递归方式的一半。

图3

图3 当n=m=2500，u=b=6 356 752.314 m时，优化前后项数对比

Fig.3 When n=m=2500, u=b=6 356 752.314 m, the comparison of the number of items before and after optimization

在图4中，当u=b=6 356 752.314 m时，图4（a）是使用原始递归方式计算、、时得到的|LDE|，且图4（a）计算LDE时所需的收敛项k和图3第一行中的项数相对应，其项数k的最大值为80。图4中的色带的值表示每个阶数和次数对应的有效数字数。显然，图4（a）中绝对精度的下降与图3最顶部一行的主要特征相对应。需要的条件越多，结果越差。且由图4（a）可知，原始递归方式满足第二类缔合勒让德微分方程，但是随着计算阶数n和次数m的不断增大及收敛项k增加，收敛精度降低特别显著。

图4

图4 当n=m=2500，u=b=6 356 752.314 m时，优化前后LDE对比

Fig.4 When n=m=2500, u=b=6 356 752.314 m, the comparison of LDE before and after optimization

使用同时减小α、β项的优化递归方式计算后，其所需要的项数k约为原始递归的一半，所有n、m组合计算得到的|LDE|的平均值在10-9数量级，如图4（b）所示，其中计算LDE所需的收敛项k和图3第二行中的项数相对应，其项数k的最大值为30。

使用只减小β项的优化递归方式的计算结果与使用同时减小α、β项的优化递归方式的计算结果相似，其所需项数也约为原始递归的一半，与原始递归相比，其计算得到的|LDE|的平均值集中在10-9数量级，也没有精度降低的现象，如图4（c）所示。图4（c）中计算LDE所需的收敛项k和图3中第三行的项数相对应，其项数k的最大值也为30。

4.2 使用同时减小α、β项的优化递归方式计算Jekeli函数在不同高度下的幅值

为了比较在不同高度下Jekeli函数的适用情况，选择在椭球表面及椭球上一定高度处对Jekeli函数进行数值计算，这里n=m=3000。

图5表示使用同时减小α、β项的优化递归方式计算得到的Jekeli函数在两个高度下的幅值。顶部为零高度，即在椭球面上；底部一行为GOCE高度，h=250 km。由单独的每一幅图中可以看出，随着阶数和次数的增加，函数开始表现为一个与阶数相关的函数。此外，由图5中不同高度下的计算结果可以看出，距离椭球体越远，Jekeli函数的非各向同性就越小；同时，随着阶次和高度的增加，图中的色带值表明了递归计算过程中没有出现溢出现象。

图5

图5 同时减小α、β项的优化递归方式下Jekeli函数在两个高度的值（顶部（u=b=6 356 752.314 m），底部（u=b+250 000 m））

Fig.5 The value of the Jekeli function at two heights in the optimized recursive mode of increasing the γ item (the top (u=b=6 356 752.314 m), the bottom (u=b+250 000 m))

4.3 使用只减小β项的优化递归方式计算Jekeli函数在不同高度下的幅值

与上述情况相同使用只减小β项的优化方式计算Jekeli函数在不同高度的幅值，结果如图6所示。

图6

图6 只减小β项的优化递归方式下，Jekeli函数在两个高度的值（顶部（u=b=6 356 752.314 m），底部（u=b+250 000 m））

由图6可知，随着阶数和次数的增加，Jekeli函数也表现为一个与阶数相关的函数。图6中的色带同样表明，随着阶次和高度的增加递归计算过程中没有出现溢出现象。

4.4 不同递归方式的计算时间及适用阶数

当计算高度为参考椭球表面高度，阶数n取不同值时（m≤n），采用原始递归方式计算的结果见表1，可以看出，原始递归方式的计算阶数不能超过2500阶，因为当n=3000阶时计算得到的|LDE|的平均值为2.06×103，已经远远大于0，不能满足绝对精度无限接近0的要求，所以原始递归不适用于高阶计算。

表1 原始递归|LDE|的计算结果The results of computing the primitive recursion |LDE|原始递归的|LDE|n=1500n=2000n=2500n=3000最大值1.90×10-23.453.75×1034.20×106平均值3.79×10-63.10×10-22.522.06×103

新窗口打开| 下载CSV

当计算高度为参考椭球表面高度，阶数n取不同值时（m≤n），由表2可知，当2500<n≤5000时，使用同时减小α、β项或使用只减小β项的优化递归方式都可以提供较高的计算精度，因为当n≤5000时，使用这两种优化递归方式计算不同高度和阶次的Jekeli函数时都没有出现溢出现象，并且使用这两种优化递归方式计算得到的|LDE|的值也满足精度要求；当计算阶数n达到10 000时，同时减小α、β项的优化递归方式不适用，因为其计算得到的|LDE|平均值为4.95×107远远大于0，不满足绝对精度的要求，而使用只减小β项的优化递归方式计算得到的|LDE|平均值为5.00×10-2是以10-2数量级接近于0，所以当计算较高阶时，即当阶数n接近10 000时应使用只减小β项的优化递归方式。

表2 优化递归后|LDE|的计算结果

Calculation results of |LDE| after optimized recursion

阶数n同时减小α、β项的优化递归的|LDE|只减小β项的优化递归的|LDE|最大值平均值最大值平均值15001.86×10-83.56×10-101.12×10-81.96×10-1030001.81×10-54.16×10-87.15×10-71.56×10-850008.95×10-16.12×10-46.10×10-57.60×10-710 0003.38×10114.95×1071.255.00×10-2

新窗口打开| 下载CSV

在比较不同阶数优化前后的计算时间时，为了保持其他变量恒定，本文所采用的计算机硬件配置（Windows11 64位操作系统，R7-7735H CPU@3.20 GHz，16 GB内存）保持不变，当计算高度为u=b=6 356 752.314时，由表3可知，当计算阶数依次为n=1500、2000、2500时，其中m≤n，使用原始递归方式计算、、的时间分别为15.247、32.961、57.399 s，而使用同时减小α、β项或使用只减小β项的优化递归方式，其计算时间比使用原始递归方式减小了约一半。

表3 优化前后针对不同阶数，计算函数Sn，m的0、1、2阶所需时间对比

Tab.3 Comparison of the time required to calculate the 0th, 1st, and 2nd order derivatives of the function Sn,m for different degrees before and after optimization

阶数n原始递归同时减小α、β项的优化递归只减小β项的优化递归150015.2477.3249.100200032.96114.12617.589250057.39923.73329.7915000—124.493151.84410 000—724.021820.401

新窗口打开| 下载CSV

综合考虑计算时间及计算结果的绝对精度，对于同时减小α、β项的优化递归方式应考虑阶数n≤5000，而当阶数n在10 000附近时应使用只减小β项的优化递归方式。

5 结论

首先关于不同递归方式的适用阶数，考虑到Jekeli函数递归计算必须满足的绝对精度的要求，可以得出结论：原始递归方式适用于n≤2500，当2500<n≤5000时，只有减小β项的优化递归方式和同时减小α、β项的优化递归方式可以提供较高的计算精度，而当n接近10 000时，须使用只减小β项的优化递归方式计算才能满足精度要求。

在Jekeli函数中，高斯超几何函数求和的项数基本控制了递归的速度，本文充分利用这一特性，对高斯超几何函数进行线性变换，并详细推导了两种优化递归方式。这两种方式都能有效减少递归计算所需的项数，从而大幅缩短了计算时间。同时，通过数值计算证明了这两种优化方式在计算精度方面均优于原始递归方式，为第二类勒让德函数的计算提供了更为高效的解决方案。

初审：张 琳复审：宋启凡
终审：金 君

资讯

○