超高阶次Legendre函数的计算是地球重力场以及相关领域中的重要研究课题，国内外有许多学 者(Colombo，1981; Fantino and Casotto， 2009; Fukushima，2012; Gleason，1985; Holmes and Featherstone， 2002; Jekeli et al., 2007; Wenzel，1998; Wittwer et al., 2008; 刘缵武等，2011; 彭富清和夏哲仁，2004; 王建国等，2009; 吴星和刘雁雨，2006; 张传定等，2004)都讨论过超高阶次Legendre函数的计算问题.从目前的研究成果来看，计算Legendre函数的主要途径仍然是利用Legendre函数的性质对其进行递推计算.若按使用的Legendre函数的递推公式来分类的话，则递推方法大致分为三类(刘缵武等，2011;王建国等，2009; 吴星和刘雁雨，2006)：第一是按阶(次)数递推方法，该方法的优点是计算方便和用时较少，但缺点是该方法会导致计算数据向下溢出，例如在双精度数范围内大致在计算至1900阶连带Legendre函数时就会产生数据向下溢出现象.第二是Belikov递推方法.第三是跨阶数递推方法.就这三种递推算法而言，吴星和刘雁雨(2006)、王建强等(2009)都曾进行过分析，特别是王建强等(2009)对阶数3060时计算了P_3060，m(0.1)与P_3060，m(0.9)的值，得出了跨阶数递推方法更适合用于超高阶次Legendre函数计算的结论.然而在目前对超高阶次Legendre函数计算的研究中大多算法仍然仅限于按阶(次)数递推方法的研究上.

事实上，目前关于超高阶次Legendre函数计算的研究主要还是集中在按阶(次)数递推的改进方法上，而改进的途径体现在两个方面：第一是在递推过程中分离出奇异因子sin^mθ，从而在双精度数范围内增加递推计算Legendre函数的阶次数，例如Holmes和Featherstone(2002)利用该方法将Legendre函数计算至2700完整阶次；第二是采用扩充双精度数域的方法(Fukushima，2012; Jekeli et al., 2007; Wittwer et al., 2008)，例如可以将双精度数的存储空间扩充四倍从而使得数据下溢向后 推延，但这样做的结果是计算用时也相应地增大(大约增加50倍用时以上).这里要特别指出，Fukushima(2012)引入了X-数域的方法，该方法是将用两个双精度数来表示一个X-数，其中第一个双精度数表示某个双精度数的有效数字，而第二个双精度数则表示该数的指数，这样就构成了一个X-数，显然X-数表示数值的范围就扩大了很多，而且计算用时也比直接将双精度数扩充四倍的方法要节约许多.

此外还有许多关于Legendre函数性质的研究工作.例如Jekeli等(2007)证明了连带Legendre函数P _n，m(cosθ)当m≥nsinθ+μ时可以近似地等于零，这里μ是一个正整数(例如：μ=15)，这样的简化并不会影响球谐级数的计算结果，从而就极大地减少了计算Legendre函数的用时.又如 Cheong等(2012)、Sneeuw和Bun(1996)、Swarztrauber(1993)还研究了连带Legendre函数的Fourier展 开式，特别是Cheong等(2012)给出了连带Legendre函数的 Fourier展开系数的算法，为Legendre函数的计算提供另一种途径.除了递推算法外，Yu等(2012)还研究了Legendre函数的积分表达式，给出了连带Legendre函数一种积分表达式，其中被积函数是有界的，这样可以在双精度数的范围内对任意阶次的Legendre函数进行计算.

以上简要地概括了关于超高阶次Legendre函数计算的一些主要结果，但这些结果存在着一些缺陷：例如，若采用扩充双精度数域的方法，则计算用时增加很多；若仍以双精度数作为基准，则计算过程中难以评价是否会因数据溢出而产生的发散等现象.本文的目的就是要基于跨阶数递推方法来讨论在双精度数范围内超高阶次连带Legendre函数的计算问题.在本文计算Legendre函数时将阶数取为20000，对P _n，m(cosθ)以1°为间隔从θ=1°至θ=89°进行了计算，并通过多种途径对计算结果的准确性进行了验证. 2 Legendre函数的计算公式 2.1 关于Legendre函数的记号

记P_n(x)是n阶Legendre多项式，即：

因为在球谐级数计算中总有x=cosθ，这里θ∈[0，π]是归化余纬，所以使用符号时不再区分P _n，m(x)和P _n，m(cosθ)，并记成P _n，m.又因为连带Legendre函数P _n，m(x)在区间[－1，+1]中当n－m是偶数时是偶函数，当n－m是奇数时是奇函数，因此下面计算时仅需讨论归化余纬θ∈[0，π/2]的情况. 2.2 按阶数递推算法

关于P _n，m的按阶数递推算法的公式为(Fantino and Casotto， 2009; Fukushima，2012; Holmes and Featherstone， 2002; Jekeli et al., 2007)

事实上，按阶数递推公式(5)是计算Legendre函数的基本算法，该算法在θ接近90°(即：赤道附近)具有很好的收敛性. 2.3 跨阶数递推公式

引入记号：当m>n时

由于公式(8)是按跨阶数给出的，所以公式(8)被称为Legendre函数的跨阶数递推公式. 2.4 当m=0和m=1时P _n，m的计算

因为公式(8)仅对m≥2成立，所以还需讨论m=0和m=1时P _n，m的计算方法.事实上，此时P _n，0和P _n，1可由公式(5)给出，即：

第一步，引入初始值，由此利用公式(12)与(13)便可计算P _n，0与P _n，1；第二步，再利用初始值P_0，0，P_1，0，P_0，1，P_1，1以及P _n，0与P _n，1，按公式(8)便可递推计算当n≥2和m≥2时P _n，m的值.称上述递推计算P _n，m的方法为跨阶数递推算法. 3 数值计算与精度分析

利用上述给出的跨阶数递推算法可以对P _n，m(cosθ)在双精度数范围内进行计算，本文对P _n，m(cosθ)的值计算至完整20000阶次，其中余纬以1°为间隔从θ=1°变换至θ=89°.为了与P _n，m的理论值进行区分，下面使用P _n，m^(r)表示利用跨阶数递推算法计算得到的P _n，m的值.为了分析P _n，m^(r)的计算精度，下面将用几种不同的检验途径.

例1 首先引入关于P _n，m的恒等式如下

图1

Fig. 1

图 1 P n，m的平方和的精度 Fig.1 The accuracies for summation of square of P n，m

例2 例1中的η(θ)反映的是同一阶数内Legendre函数值的平方和的精度，该精度并不能用于评价某个特定阶次Legendre函数的计算精度.如何评价每个Legendre函数的计算精度呢？注意到P _n，m^(r)(cosθ)的值是经过递推计算的，因此只需评价P _20000，m^(r)(cosθ)或 ^(r)_19999，m(cosθ)的精度即可，这里0≤m≤20000.为此引入Legendre函数的积分公式如下(Yu et al., 2012)： 当n－m=2k是偶数时，则有

图2

Fig. 2

图 2 δm(θ)的分布 Fig. 2 The distributions of δm(θ)

对于 ^(r)_19999，m(cosθ)而言(m是奇数)，同样能够给出类似图 2的计算结果，这里不再叙述.而对于m是奇数时可利用Legendre函数的性质(党诵诗，1988)

总之，图 2描述了P _20000，m^(r)(cosθ)中每项的精度.由图 2可见，对于中高纬度的项来说P _20000，m^(r)(cosθ)的计算精度均达到了10^－10，但在赤道附近δ_m(θ)的值则要相应地大了一些.虽然如此，但在赤道附近P _20000，m^(r)(cosθ)的计算精度均也能达到10^-9，这样的精度对于超高阶次Legendre函数值的计算来说是能够满足要求的.

例3 是否在赤道附近P _20000，m^(r)(cosθ)的计算精度就会差些呢？或者说，为何δ_m(θ)的值在赤道附近会比在中高纬度地区要大一个量级呢？事实上，产生这样情况的原因是利用积分方法计算的P _20000，m^(r)(cosθ)值在赤道附近的精度较低而引起的.为证实这样的结论，下面在赤道附近使用按阶数递推算法来计算Legendre函数，理由是在赤道附近按阶数递推算法具有较好的收敛性.本文专门针对赤道附近(θ=87°、89.5°)利用按阶数递推公式(5)计算了P _20000，m(cosθ)的值，记为 ^(d)_20000，m(cosθ)，并统计了 P _20000，m^(r)(cosθ)－ ^(d)_20000，m(cosθ)的值(m=0，1，…，20000)，结果由图 3给出.

图3

Fig. 3

图 3 |P 20000，m(r)(cosθ)－P 20000，m(d)(cosθ)|的分布 Fig. 3 The distributions of |P 20000，m(r)(cosθ)－P 20000，m(d)(cosθ)|

从图 3可以断言，在赤道附近P _20000，m^(r)(cosθ)每项的计算精度均可达到10^-11.总之，结合例2与例3可知，跨阶数递推算法计算的P _n，m(cosθ)的绝对精度要优于10^-10.

例4 本文在双精度数范围内研究了超高阶次Legendre函数的计算问题，并针对阶数为20000时对余纬θ以1°为间隔从θ=1°至θ=89°计算了Legendre函数的值.通常情况下，20000阶次的Legendre函数在纬度方向上对应的分辨率为 180/20000 ×1°=32.4″，而在θ接近于0°时P _n，m(cosθ)的值极易产生向下溢出，从而在计算过程中容易导致P _n，m(cosθ)会出现不稳定的现象.是否本文使用的Legendre函数的跨阶次算法在计算至20000阶时也会出现这种不稳定现象呢？为此本文专门针对θ=30″与θ=60″使用跨阶次算法计算了P _n，m(cosθ)的值，并统计了相应δ_m(θ)的分布，结果由图 4给出.

图4

Fig. 4

图 4 θ=30″与θ=60″时δm(θ)的分布 Fig. 4 The distributions of δm(θ) when θ=30″ and θ=60″

从图 4可见，在两极附近使用跨阶次算法计算了P _n，m(cosθ)至20000阶次时在分辨率范围内该算法是稳定的，而且Legendre函数计算值的绝对精度至少优于10^-10量级. 4 球谐级数合成的精度分析

本节讨论球谐级数合成的精度.假设地球引力位的球谐级数展开式为

若采用跨阶数递推算法去计算完整到20000阶次的Legendre函数，则使用公式(21)的实际计算值为

引入关于位系数的Kaula假设(Kaula， 1966，2000)，即：

如果P _n，m^(r)(cosθ)的计算精度ε满足：ε≤10^－9，则利用跨阶数递推算法计算的精度便可达到 ·O(10^－14).由于×10^－9约为1 μGal，因此可得结论：利用跨阶数递推算法从球谐展开式计算引力的精度可达到10^-5 μGal的量级.显然，这样的精度足以满足重力观测的精度要求. 5 计算时间

由于按阶数递推是计算Legendre函数的常用方法，所以本文需要去比较跨阶数递推算法与按阶数递推算法的计算用时.表 1给出了在同一台计算机上用跨阶数递推算法与按阶数递推算法计算不同阶数Legendre函数从θ=1°到89°的时间.由于按阶 数递推算法会遇到数值向下溢出的情况，所以在出现数据向下溢出时总是对溢出变量赋予了一个小数以便能够统计出按阶数递推算法的计算用时.

表1(Table 1)

表 1 计算用时比较 Table 1 The comparison of run-times 阶数 按阶数递推用时(s) 跨阶数递推用时(s) 用时比值 2000 7.2 11.8 1.64 4000 28.3 45.5 1.61 10000 176.2 278.5 1.58 20000 703.2 1080.8 1.54

表 1 计算用时比较 Table 1 The comparison of run-times

从表 1可见，跨阶数递推算法的计算用时大约是按阶数递推算法的1.6倍.其原因是按阶数递推算法公式(5)是从两个变量去推算下一个变量，而跨阶数递推算法公式(8)则是从三个变量去推算下一个变量，因此跨阶数递推算法的用时应该大约是按阶数递推算法的1.6倍左右.

与扩充双精度数域来计算Legendre函数算法的计算用时进行比较的话，明显可以看出，跨阶数递推算法的计算用时是完全可以接受的. 6 结论

本文基于Legendre函数的跨阶数递推公式研究了相应的计算Legendre函数的方法，使用该算法在双精度数范围将Legendre函数计算至完整20000阶次，实现了在双精度数范围内快速计算超高阶次Legendre函数的方法.

使用多种方法检验了使用跨阶数递推算法计算Legendre函数的精度，得出了该算法计算的每个Legendre函数的绝对精度都优于10^-10的结论.

统计了跨阶数递推算法计算Legendre函数的用时，相对于常用的按阶数递推算法而言，跨阶数递推算法的计算用时大约是按阶数递推算法的1.6倍.

在计算Legendre函数的绝对精度是10^-9的前提下，借助于Kaula准则证明了在进行球谐级数合成时其精度大致可以达到计算量的10^-14左右.例如，若使用球谐级数去计算引力值，则引力值的计算精度将会达到10^-5 μGal左右.

最后说明一下，Legendre函数的跨阶数递推公式(8)并非本文首次给出，跟踪的国外最早文献是Swarztrauber(1993)，而国内则由张传定等(2004)最早给出，但将该公式用于超高阶次Legendre函数的计算，并对计算结果的精度进行检验与验证则是本文首次实现的.根据检验与验证的结果，可以得出结论：跨阶数递推算法不仅可以在双精度数范围内用于超高阶次Legendre函数的计算，而且计算结果准确、计算用时还完全可以接受.总之，本文提出的算法完全满足超高阶次Legendre函数的要求.