超高阶次勒让德函数递推计算中的压缩因子和Horner求和技术

摘要：超高阶次（如2000阶次）缔合勒让德函数值的递推计算,在接近两极时达到极大的数量级（超过10的数千次方）,这导致现有递推方法在计算缔合勒让德函数值及其导数值时失效。我们通过插入压缩因子技术,提出一个修改的递推算法,并结合使用Horner求和技术计算球谐级数的部分和。试验表明,该算法至少可以计算到3600完全阶次的球谐级数式,优于现有结果。     

多种超高阶次缔合勒让德函数计算方法的比较

扩展高阶和超高阶重力场模型的构制与应用的数值稳定性取决于超高阶次缔合勒让德函数的计算方法。文中详细介绍了现有的多种缔合勒让德函数的递推计算方法:标准前向列推法、标准前向行推法、跨阶次递推法和Belikov列推法。从计算速度、计算精度和计算溢出问题3个角度分析比较了阶次高至2 160阶的各种方法的优劣。通过数值试验证明,Belikov列推法和跨阶次递推法是计算超高阶次缔合勒让德函数较优的方法,而其他几种方法不能用于超高阶次缔合勒让德函数的计算。文中结论为超高阶次球谐综合与球谐分析的数值计算提供了可靠的依据。     

多种超高阶次缔合勒让德函数计算方法的比较

扩展高阶和超高阶重力场模型的构制与应用的数值稳定性取决于超高阶次缔合勒让德函数的计算方法.文中详细介绍了现有的多种缔合勒让德函数的递推计算方法:标准前向列推法、标准前向行推法、跨阶次递推法和Belikov列推法.从计算速度、计算精度和计算溢出问题3个角度分析比较了阶次高至2 160阶的各种方法的优劣.通过数值试验证明,Belikov列推法和跨阶次递推法是计算超高阶次缔合勒让德函数较优的方法,而其他几种方法不能用于超高阶次缔合勒让德函数的计算.文中结论为超高阶次球谐综合与球谐分析的数值计算提供了可靠的依据.     

采用Belikov列推和跨阶次递推方法计算超高阶缔合勒让德函数

超高阶球谐重力场模型的精确构制与快速计算取决于缔合勒让德函数的计算方法。在前人研究的基础上,文中对适合超高阶缔合勒让德函数计算的Belikov列推和跨阶次递推方法进行介绍,为验证精度,通过两种途径对计算结果进行检验,并比较其计算速度。结果表明,采用两种算法得到的每个勒让德函数的绝对精度均优于10-12,在低阶,跨阶次递推方法的计算用时大约是Belikov列推法的2倍,随着阶数的升高,跨阶次递推算法表现出明显的速度优势。     

球谐函数级数式的改化形式

首先将级数式中缔合勒让德函数的计算转化为一般勒让德函数的导数计算，并给出相应的递推公式；其次将球谐函数级数式转化为仅含带谐项的级数式，并给出相应系数的计算公式。经实例计算，用机时间节省约90％。这对于应用部门快速计算重力场的有关数据将十分有利。     

完全正常化缔合勒让德函数及其导数与积分的递推关系

在地球重力场问题中,常用到完全正常化缔合勒让德函数及其导数、积分的递推关系。当前流行的地球扰动位模型均采用完全正常化的缔合勒让德函数,用此类模型可以高效方便计算各种扰动重力场元。随着本世纪多个新一代卫星重力探测计划成功实施,高阶或超高阶地球重力场模型的研究备受学界的关注。有关完全正常化缔合勒让德函数的递推关系对于高阶重力场模型具有特别意义。本文在前人研究的基础上,用初等微积分导出了若干新的递推关系式。同时还推导了正常化缔合勒让德函数及其导数、积分的检核式,这些检核式涉及地球位的球谐级数的数学性质。     

超高阶重力场模型的研究与实现

超高阶重力场模型的提出使大地水准面差距的计算结果越来越精确。本文主要介绍了采用重力场模型计算大地水准面差距的算法,采用标准向前列递推法求解缔合勒让德函数值,增大勒让德函数前三项,避免产生超高阶递推过程的不稳定现象。结果表明:增大勒让德函数前三项能稳定递推至2 700阶,验证了大地水准面差距验证算法、软件的正确性与可行性。     

第二类连带勒让德函数及其一阶、二阶导数的递推计算方法

推导出了地球重力场位模型椭球谐级数表达式中的第二类连带勒让德函数及其一阶、二阶导数的修正Jekeli递推计算方法,并与传统的Jekeli递推计算方法的结果进行比较。结果表明,在递推计算收敛精度相同的条件下,修正Jekeli递推计算比传统Jekeli递推计算需要的收敛项数减少一半,最大约为30项;随着阶次n、m的增大以及收敛项k的增加,修正Jekeli递推计算的第二类连带勒让德微分方程的精度一直保持在1×10-6左右;不同高度下,阶数n与高度h之间的关系与球近似下(R/r)n+1的阶数n与高度h之间的关系相似。     

完全规格化缔合勒让德函数及其导数的递推算法与适用范围比较

完全规格化缔合勒让德函数及其导数常用标准向前列递推算法和标准向前行递推算法进行计算。基于第一、第二相对数值精度标准对两种算法的适用范围进行分析比较,计算结果表明,标准向前列递推算法的适用范围大于标准向前行递推算法,说明前者优于后者;结果同时还表明,完全规格化缔合勒让德函数与其导数同一种算法的适用范围也相同,并指出了二者适用范围相同的原因。     

基于球谐分析的EMM-740的计算及实现

针对国内对高精度增强地磁场模型(enhanced magnetic model,EMM)在基于球谐分析时对阶数存在认知偏差、实现模型的软件研究较少以及计算精度不高等问题,考虑了地磁场7要素的关系,建立了基于球谐分析740阶次的EMM模型,给出了Schmidt半标准化缔合勒让德函数,给出了EMM模型的计算步骤,并用MATLAB实现EMM2015模型(2000～2019年)在地球上任意地磁场7要素的计算和等值线图的绘制,将计算值与美国国家海洋和大气管理局公布模型的运行数据进行误差对比分析,各元素的均方根偏差最大为0.86 nT或0.4′,比720阶的EMM模型的精度提高了3倍。结果表明,提供的EMM2015模型软件实现方法具有较高的计算精度,值得推广和应用。     

A unified approach to the Clenshaw summation and the recursive computation of very high degree and order normalised associated Legendre functions

 Spherical harmonic expansions form partial sums of fully normalised associated Legendre functions (ALFs). However, when evaluated increasingly close to the poles, the ultra-high degree and order (e.g. 2700) ALFs range over thousands of orders of magnitude. This causes existing recursion techniques for computing values of individual ALFs and their derivatives to fail. A common solution in geodesy is to evaluate these expansions using Clenshaw's method, which does not compute individual ALFs or their derivatives. Straightforward numerical principles govern the stability of this technique. Elementary algebra is employed to illustrate how these principles are implemented in Clenshaw's method. It is also demonstrated how existing recursion algorithms for computing ALFs and their first derivatives are easily modified to incorporate these same numerical principles. These modified recursions yield scaled ALFs and first derivatives, which can then be combined using Horner's scheme to compute partial sums, complete to degree and order 2700, for all latitudes (except at the poles for first derivatives). This exceeds any previously published result. Numerical tests suggest that this new approach is at least as precise and efficient as Clenshaw's method. However, the principal strength of the new techniques lies in their simplicity of formulation and implementation, since this quality should simplify the task of extending the approach to other uses, such as spherical harmonic analysis. Received: 30 June 2000 / Accepted: 12 June 2001     

The integral formulas of the associated Legendre functions

A new kind of integral formulas for ${\bar{P}_{n,m} (x)}$ is derived from the addition theorem about the Legendre Functions when n ? m is an even number. Based on the newly introduced integral formulas, the fully normalized associated Legendre functions can be directly computed without using any recursion methods that currently are often used in the computations. In addition, some arithmetic examples are computed with the increasing degree recursion and the integral methods introduced in the paper respectively, in order to compare the precisions and run-times of these two methods in computing the fully normalized associated Legendre functions. The results indicate that the precisions of the integral methods are almost consistent for variant x in computing ${\bar{P}_{n,m} (x)}$ , i.e., the precisions are independent of the choice of x on the interval [0,1]. In contrast, the precisions of the increasing degree recursion change with different values on the interval [0,1], particularly, when x tends to 1, the errors of computing ${\bar{P}_{n,m} (x)}$ by the increasing degree recursion become unacceptable when the degree becomes larger and larger. On the other hand, the integral methods cost more run-time than the increasing degree recursion. Hence, it is suggested that combinations of the integral method and the increasing degree recursion can be adopted, that is, the integral methods can be used as a replacement for the recursive initials when the recursion method become divergent.     

On the computation and approximation of ultra-high-degree spherical harmonic series

Spherical harmonic series, commonly used to represent the Earth’s gravitational field, are now routinely expanded to ultra-high degree (> 2,000), where the computations of the associated Legendre functions exhibit extremely large ranges (thousands of orders) of magnitudes with varying latitude. We show that in the degree-and-order domain, (ℓ,m), of these functions (with full ortho-normalization), their rather stable oscillatory behavior is distinctly separated from a region of very strong attenuation by a simple linear relationship: , where θ is the polar angle. Derivatives and integrals of associated Legendre functions have these same characteristics. This leads to an operational approach to the computation of spherical harmonic series, including derivatives and integrals of such series, that neglects the numerically insignificant functions on the basis of the above empirical relationship and obviates any concern about their broad range of magnitudes in the recursion formulas that are used to compute them. Tests with a simulated gravitational field show that the errors in so doing can be made less than the data noise at all latitudes and up to expansion degree of at least 10,800. Neglecting numerically insignificant terms in the spherical harmonic series also offers a computational savings of at least one third.     

Numerical computation of spherical harmonics of arbitrary degree and order by extending exponent of floating point numbers

By extending the exponent of floating point numbers with an additional integer as the power index of a large radix, we compute fully normalized associated Legendre functions (ALF) by recursion without underflow problem. The new method enables us to evaluate ALFs of extremely high degree as 2³² =  4,294,967,296, which corresponds to around 1 cm resolution on the Earth’s surface. By limiting the application of exponent extension to a few working variables in the recursion, choosing a suitable large power of 2 as the radix, and embedding the contents of the basic arithmetic procedure of floating point numbers with the exponent extension directly in the program computing the recurrence formulas, we achieve the evaluation of ALFs in the double-precision environment at the cost of around 10% increase in computational time per single ALF. This formulation realizes meaningful execution of the spherical harmonic synthesis and/or analysis of arbitrary degree and order.     

Fourier-series representation and projection of spherical harmonic functions

Computations of Fourier coefficients and related integrals of the associated Legendre functions with a new method along with their application to spherical harmonics analysis and synthesis are presented. The method incorporates a stable three-step recursion equation that can be processed separately for each colatitudinal Fourier wavenumber. Recursion equations for the zonal and sectorial modes are derived in explicit single-term formulas to provide accurate initial condition. Stable computations of the Fourier coefficients as well as the integrals needed for the projection of Legendre functions are demonstrated for the ultra-high degree of 10,800 corresponding to the resolution of one arcmin. Fourier coefficients, computed in double precision, are found to be accurate to 15 significant digits, indicating that the normalized error is close to the machine round-off error. The orthonormality, evaluated with Fourier coefficients and related integrals, is shown to be accurate to O(10^?15) for degrees and orders up to 10,800. The Legendre function of degree 10,800 and order 5,000, synthesized from Fourier coefficients, is accurate to the machine round-off error. Further extension of the method to even higher degrees seems to be realizable without significant deterioration of accuracy. The Fourier series is applied to the projection of Legendre functions to the high-resolution global relief data of the National Geophysical Data Center of the National Oceanic and Atmospheric Administration, and the spherical harmonic degree variance (power spectrum) of global relief data is discussed.     

Fourier transform summation of Legendre series and D-functions

The relation between D- and d-functions, spherical harmonic functions and Legendre functions is reviewed. Dmatrices and irreducible representations of the rotation group O(3) and SU(2) group are briefly reviewed. Two new recursive methods for calculations of D-matrices are presented. Legendre functions are evaluated as part of this scheme. Vector spherical harmonics in the form af generalized spherical harmonics are also included as well as derivatives of the spherical harmonics. The special dmatrices evaluated for argument equal to/2 offer a simple method of calculating the Fourier coefficients of Legendre functions, derivatives of Legendre functions and vector spherical harmonics. Summation of a Legendre series or a full synthesis on the unit sphere of a field can then be performed by transforming the spherical harmonic coefficients to Fourier coefficients and making the summation by an inverse FFT (Fast Fourier Transform). The procedure is general and can also be applied to evaluate derivatives of a field and components of vector and tensor fields.     

Methods of harmonic synthesis for global geopotential models and their first-, second- and third-order gradients

Four widely used algorithms for the computation of the Earth’s gravitational potential and its first-, second- and third-order gradients are examined: the traditional increasing degree recursion in associated Legendre functions and its variant based on the Clenshaw summation, plus the methods of Pines and Cunningham–Metris, which are free from the singularities that distinguish the first two methods at the geographic poles. All four methods are reorganized with the lumped coefficients approach, which in the cases of Pines and Cunningham–Metris requires a complete revision of the algorithms. The characteristics of the four methods are studied and described, and numerical tests are performed to assess and compare their precision, accuracy, and efficiency. In general the performance levels of all four codes exhibit large improvements over previously published versions. From the point of view of numerical precision, away from the geographic poles Clenshaw and Legendre offer an overall better quality. Furthermore, Pines and Cunningham–Metris are affected by an intrinsic loss of precision at the equator and suffer from additional deterioration when the gravity gradients components are rotated into the East-North-Up topocentric reference system. Electronic supplementary material  The online version of this article (doi:) contains supplementary material, which is available to authorized users.     

A method for computing the coefficients in the product-sum formula of associated Legendre functions

The product of two associated Legendre functions can be represented by a finite series in associated Legendre functions with unique coefficients. In this study a method is proposed to compute the coefficients in this product-sum formula. The method is of recursive nature and is based on the straightforward polynomial form of the associated Legendre function's factor. The method is verified through the computation of integrals of products of two associated Legendre functions over a given interval and the computation of integrals of products of two Legendre polynomials over [0,1]. These coefficients are basically constant and can be used in any future related applications. A table containing the coefficients up to degree 5 is given for ready reference.     

Non-Singular Expressions for the Gravity Gradients in the Local North-Oriented and Orbital Reference Frames

The conventional expansions of the gravity gradients in the local north-oriented reference frame have a complicated form, depending on the first- and second-order derivatives of the associated Legendre functions of the colatitude and containing factors which tend to infinity when approaching the poles. In the present paper, the general term of each of these series is transformed to a product of a geopotential coefficient and a sum of several adjacent Legendre functions of the colatitude multiplied by a function of the longitude. These transformations are performed on the basis of relations between the Legendre functions and their derivatives published by Ilk (1983). The second-order geopotential derivatives corresponding to the local orbital reference frame are presented as linear functions of the north-oriented gravity gradients. The new expansions for the latter are substituted into these functions. As a result, the orbital derivatives are also presented as series depending on the geopotential coefficients multiplied by sums of the Legendre functions whose coefficients depend on the longitude and the satellite track azimuth at an observation point. The derived expansions of the observables can be applied for constructing a geopotential model from the GOCE mission data by the time-wise and space-wise approaches. The numerical experiments demonstrate the correctness of the analytical formulas.An erratum to this article can be found at     

球冠谐分析中非整阶Legendre函数的性质及其计算

局部重力场的谱方法是当前重力学的研究方向,该方法的核心问题是如何构造合适的谱函数以及如何对谱函数实施快速、有效的计算。当所研究的区域近似一个球冠时,迂冠谐函数是该区域对应的谱函数,它由非整阶勒让德（Ｌｅｇｅｎｄｒｅ）函数和三角函数组成,显然非整阶勒让德函数的构造和计算是研究球冠谐函数的关键。本文研究了非高小阶勒让德函数的性质和实用计算方法,包括如何对非整阶勒让德函数实施规格化处理。