从20世纪80年代开始，学者们认识到地球扁率在求解重力场边值问题中的重要性，于是在建立地球重力场球谐级数模型时，增加了椭球改正项，并开展了对椭球谐级数的研究。当前国际上关于各种地球重力场模型的计算都是以球近似问题的解加上椭球改正来得到^[1-3]。为了确定1 cm级精度大地水准面^[4]，越来越多的学者致力于椭球域大地边值问题的研究^[5-8]。由于地球的几何形状更接近于椭球，如果利用椭球谐直接逼近计算，逼近误差应该小于球谐逼近所带来的误差^[9]，这样会提高椭球面边值问题的解算精度，使全球1 cm级精度的大地水准面的实现成为可能。

在地球重力场位模型椭球谐级数展开式中，最主要的问题是关于第二类连带勒让德函数$ {Q_{nm}}$的计算方法。目前，大多数文献对第二类连带勒让德函数$ {Q_{nm}}$的计算主要有两种解决办法：一是将$ {Q_{nm}}$以${{\rm{e}}^2} $的形式直接进行展开^[10-11]；二是将$ {Q_{nm}}$用Gauss超几何函数的展开形式进行计算^[12-13]。

本文基于Gauss超几何函数的定义，详细推导了第二类连带勒让德函数及其一阶、二阶导数的修正Jekeli递推计算公式，并与传统的Jekeli递推计算方法^[14]进行比较，证明了修正Jekeli递推计算方法的收敛速度更快，收敛精度更稳定，为解决构建地球重力场的椭球谐位模型及其泛函计算有重要的理论指导作用。

1 位椭球谐理论

对于地球外空间的任意一点P的引力位$ V(u, \vartheta , \lambda )$，可以用椭球谐级数的形式表示为^[15]：

$ \begin{array}{c}V(u, \vartheta , \lambda )=\frac{\mathrm{G}\mathrm{M}}{R}\stackrel{\mathrm{\infty }}{\sum \limits_{n=0}}\stackrel{n}{\sum \limits_{m=0}}\frac{{Q}_{nm}\left(\mathrm{i}\frac{u}{E}\right)}{{Q}_{nm}\left(\mathrm{i}\frac{b}{E}\right)}\times \\ \left({\overline{C}}_{nm}^{e}\mathrm{c}\mathrm{o}\mathrm{s}\right(m\lambda )+{\overline{S}}_{nm}^{e}\mathrm{s}\mathrm{i}\mathrm{n}(m\lambda \left)\right){\overline{P}}_{nm}\left(\mathrm{c}\mathrm{o}\mathrm{s}\vartheta \right)\end{array} $

(1)

式中，$ \mathrm{G}\mathrm{M} $为地心引力常数；$ {\mathrm{i}}^{2}=-1$；$ E=\sqrt{{a}^{2}-{b}^{2}} $为线性偏心率，$ a$为参考椭球的长半轴长；$ b $为参考椭球的短半轴长；$ \lambda $为经度；$ \vartheta$为归化纬度的余角；$ R$为地球的平均半径；$ u$为与参考椭球共焦的椭球短半轴长；$ \left\{{\overline{C}}_{nm}^{e}, {\overline{S}}_{nm}^{e}\right\}$为椭球谐系数，$ n\mathrm{、}m$分别为阶和次；$ {\overline{P}}_{nm}\left(\mathrm{c}\mathrm{o}\mathrm{s}\vartheta \right) $为第一类缔合勒让德函数；$ {Q}_{nm}\left(\mathrm{i}\frac{u}{E}\right) $为第二类连带勒让德函数。

2 第二类连带勒让德函数及其一阶、二阶导数 2.1 第二类连带勒让德函数的定义

第二类连带勒让德函数的定义有许多种，本文使用如下定义^[16]：

$ \begin{array}{c}{Q}_{nm}(n, m)={(-1)}^{m}\frac{{2}^{n}n!(n+m)!}{(2n+1)!}\frac{({z}^{2}{-1)}^{\frac{m}{2}}}{{z}^{n+m+1}}\times \\ {}_{2}{F}_{1}(\frac{n+m+2}{2}, \frac{n+m+1}{2}, n+\frac{3}{2}, \frac{1}{{z}^{2}})\end{array} $

(2)

式中，$ z=\mathrm{i}\frac{u}{E}$，$ \left|z\right|>1$；$ {}_{2}{F}_{1}\left(\mathrm{ }\right) $表示Gauss超几何函数，形式如下：

$ {}_{2}{F}_{1}(\alpha , \beta , \gamma , \delta )=\stackrel{\mathrm{\infty }}{\sum \limits_{k=0}}\frac{{\left(\alpha \right)}_{k}{\left(\beta \right)}_{k}}{{\left(\gamma \right)}_{k}}\frac{{\delta }^{k}}{k!} $

(3)

其中，$ {\left(x\right)}_{k}=\frac{\mathit{\Gamma} (x+k)}{\mathit{\Gamma} \left(x\right)}=\frac{(x+k-1)!}{(x-1)!}$称为Pochhammer符号，$ \mathit{\Gamma} \left(\mathrm{ }\right)$为Gamma函数，收敛项$ k $为正整数。显然当$ \left|\delta \right| <1 $时，式（2）中第二类连带勒让德函数$ {Q}_{nm}\left(z\right)$的级数展开式的收敛性是依式（3）Gauss超几何函数中的参数$ \alpha $、$ \beta $、$ \gamma $和$ \delta $之间的关系而定。

因此，为了得到第二类连带勒让德函数及其一阶、二阶导数更加快速、稳定的递推计算方法，可以对式（3）Gauss超几何函数中增大$ \gamma $或减小$ \alpha $、$ \beta $、$ \delta $的值来实现。

2.2 传统Jekeli递推计算

传统Jekeli递推计算的思路是利用式（2）定义的第二类连带勒让德函数，对其乘以一个常数，使其转化为第二类连带勒让德实函数$ {\overline{S}}_{n, m}\left(\frac{u}{E}\right) $^[14]：

$ \begin{array}{c}{\overline{S}}_{n, m}\left(\frac{u}{E}\right)=\frac{{\left(\frac{R}{E}\right)}^{n+1}{\mathrm{i}}^{n+1}(2n+1)!}{{2}^{n}n!}\times \\ \sqrt{\frac{{\varepsilon }_{m}}{(2n+1)(n-m)!(n+m)!}}{\overline{Q}}_{n, m}\left(z\right)\end{array} $

(4)

其中，

$ {\overline{Q}}_{n, m}\left(z\right)=\sqrt{\frac{(2n+1)(n-m)!}{{\varepsilon }_{m}(n+m)!}}{Q}_{n, m}\left(z\right) $

(5)

式中，当$ m=0 $时，$ {\varepsilon }_{m}=1 $；当$ m\ne 0$时，$ {\varepsilon }_{m}=\frac{1}{2} $。这样式（1）的椭球谐引力位可以用实函数的形式表示为：

$ \begin{array}{c}V(u, \vartheta , \lambda )=\frac{\mathrm{G}\mathrm{M}}{R}\stackrel{\mathrm{\infty }}{\sum \limits_{n=0}}\stackrel{n}{\sum \limits_{m=0}}\frac{{\overline{S}}_{n, m}\left(\frac{u}{E}\right)}{{\overline{S}}_{n, m}\left(\frac{b}{E}\right)}\times \\ \left({\overline{C}}_{nm}^{e}\mathrm{c}\mathrm{o}\mathrm{s}\right(m\lambda )+{\overline{S}}_{nm}^{e}\mathrm{s}\mathrm{i}\mathrm{n}(m\lambda \left)\right){\overline{P}}_{nm}\left(\mathrm{c}\mathrm{o}\mathrm{s}\vartheta \right)\end{array} $

(6)

式（2）结合式（4），得到$ {\overline{S}}_{n, m}\left(\frac{u}{E}\right) $的Gauss超几何函数表示形式为：

$ \begin{array}{c}{\overline{S}}_{n, m}\left(\frac{u}{E}\right)={(1+\frac{{E}^{2}}{{u}^{2}})}^{\frac{m}{2}}{\left(\frac{R}{u}\right)}^{n+1}\times \\ {}_{2}{}^{}{F}_{1}(\frac{n+m+2}{2}, \frac{n+m+1}{2}, n+\frac{3}{2}, -\frac{{E}^{2}}{{u}^{2}})\end{array} $

(7)

接下来给出$ {\overline{S}}_{n, m}\left(\frac{u}{E}\right)$及其一阶偏导数$ {\overline{S}}_{n, m}^{1}$、二阶偏导数$ {\overline{S}}_{n, m}^{2} $的Jekeli递推计算公式。

1）$ {\overline{S}}_{n, m}^{0}$的Jekeli递推计算公式为^{[14, 17]}：

$ {\overline{S}}_{n, m}^{0}={\overline{S}}_{n, m}\left(\frac{u}{E}\right)={\beta }_{n, m}\stackrel{\mathrm{\infty }}{\sum \limits_{k=0}}{\alpha }_{n, m, k} $

(8)

其中，

2）$ {\overline{S}}_{n, m}^{1}$的Jekeli递推计算公式为^[17]：

$ {\overline{S}}_{n, m}^{1}=\frac{\partial {\overline{S}}_{n, m}^{0}}{\partial u}=\frac{\partial {\beta }_{n, m}}{\partial u}\stackrel{\mathrm{\infty }}{\sum \limits_{k=0}}{\alpha }_{n, m, k}+{\beta }_{n, m}\stackrel{\mathrm{\infty }}{\sum \limits_{k=0}}\frac{\partial {\alpha }_{n, m, k}}{\partial u} $

(9)

其中，

3）$ {\overline{S}}_{n, m}^{2}$的Jekeli递推计算公式为：

$ \begin{array}{c}{\overline{S}}_{n, m}^{2}=\frac{{\partial }^{2}{\overline{S}}_{n, m}^{0}}{\partial u}=\frac{{\partial }^{2}{\beta }_{n, m}}{\partial {u}^{2}}\stackrel{\mathrm{\infty }}{\sum \limits_{k=0}}{\alpha }_{n, m, k}+\\ 2\frac{\partial {\beta }_{n, m}}{\partial u}\stackrel{\mathrm{\infty }}{\sum \limits_{k=0}}\frac{\partial {\alpha }_{n, m, k}}{\partial u}+{\beta }_{n, m}\stackrel{\mathrm{\infty }}{\sum \limits_{k=0}}\frac{{\partial }^{2}{\alpha }_{n, m, k}}{\partial {u}^{2}}\end{array} $

(10)

其中，

2.3 修正Jekeli递推计算

为了使第二类连带勒让德函数及其一阶、二阶导数的级数展开式收敛速度更快、收敛更加稳定，本文对式（3）中定义的Gauss超几何函数$ {}_{2}{}^{}{F}_{1}(\alpha , \beta , \gamma , \delta ) $进行改进。当$ \left|\delta \right| <1$时，对各参数进行调整，通过减小$ \alpha $、$ \beta $的值，同时$ \gamma $保持不变，由此给出新的Gauss超几何函数$ {}_{2}{}^{}{F}_{1}(\alpha , \beta , \gamma , \delta ) $定义，并推导出相应的第二类连带勒让德实函数及其一阶、二阶$\mathrm{*}{\overline{S}}_{n, m}^{0} $、$ \mathrm{*}{\overline{S}}_{n, m}^{1} $和$ \mathrm{*}{\overline{S}}_{n, m}^{2} $的修正Jekeli递推计算公式。新的Gauss超几何函数定义如下：

$ \begin{array}{c}{}_{2}{}^{}{F}_{1}(\alpha , \beta , \gamma , \delta )={(1-\delta )}^{\gamma -\alpha -\beta }\times \\ {}_{2}{}^{}{F}_{1}(\gamma -\alpha , \gamma -\beta , \gamma , \delta )\end{array} $

(11)

第二类连带勒让德实函数的修正Jekeli递推计算公式为：

$ \begin{array}{c}\mathrm{*}{\overline{S}}_{n, m}\left(\frac{u}{E}\right)={(1+\frac{{E}^{2}}{{u}^{2}})}^{-\frac{m}{2}}{\left(\frac{R}{u}\right)}^{n+1}\times \\ {}_{2}{}^{}{F}_{1}(\frac{n-m+2}{2}, \frac{n-m+1}{2}, n+\frac{3}{2}, -\frac{{E}^{2}}{{u}^{2}})\end{array} $

(12)

1）$ \mathrm{*}{\overline{S}}_{n, m}^{0} $的修正Jekeli递推计算公式为：

$ \mathrm{*}{\overline{S}}_{n, m}^{0}={\overline{S}}_{n, m}\left(\frac{u}{E}\right)={\beta }_{n, m}\stackrel{\mathrm{\infty }}{\sum \limits_{k=0}}{\alpha }_{n, m, k} $

(13)

其中，

$ \frac{\partial {\beta }_{n, m}}{\partial u}=-\frac{{E}^{2}{\beta }_{n, m-1}}{({u}^{2}+{E}^{2}{)}^{\frac{3}{2}}}+\frac{u}{\sqrt{{u}^{2}+{E}^{2}}}\frac{\partial {\beta }_{n, m-1}}{\partial u} $

(14)

2）$ \mathrm{*}{\overline{S}}_{n, m}^{1}$的修正Jekeli递推计算公式为：

$ \mathrm{*}{\overline{S}}_{n, m}^{1}=\frac{\partial \mathrm{*}{\overline{S}}_{n, m}^{0}}{\partial u}=\frac{\partial {\beta }_{n, m}}{\partial u}\stackrel{\mathrm{\infty }}{\sum \limits_{k=0}}{\alpha }_{n, m, k}+{\beta }_{n, m}\stackrel{\mathrm{\infty }}{\sum \limits_{k=0}}\frac{\partial {\alpha }_{n, m, k}}{\partial u} $

(15)

其中，

$ \frac{\partial {\beta }_{n, m}}{\partial u}$的计算公式同式（14）。

3）$ \mathrm{*}{\overline{S}}_{n, m}^{2} $的修正Jekeli递推计算公式为：

$ \begin{array}{c}\mathrm{*}{\overline{S}}_{n, m}^{2}=\frac{{\partial }^{2}\mathrm{*}{\overline{S}}_{n, m}^{0}}{\partial u}=\frac{{\partial }^{2}{\beta }_{n, m}}{\partial {u}^{2}}\stackrel{\mathrm{\infty }}{\sum \limits_{k=0}}{\alpha }_{n, m, k}+\\ 2\frac{\partial {\beta }_{n, m}}{\partial u}\stackrel{\mathrm{\infty }}{\sum \limits_{k=0}}\frac{\partial {\alpha }_{n, m, k}}{\partial u}+{\beta }_{n, m}\stackrel{\mathrm{\infty }}{\sum \limits_{k=0}}\frac{{\partial }^{2}{\alpha }_{n, m, k}}{\partial {u}^{2}}\end{array} $

(16)

其中，

为了检验推导的修正Jekeli递推方法计算第二类连带勒让德函数及其一阶、二阶导数的正确性，可以构建类似于第一类连带勒让德函数满足的微分方程，在$ \mathrm{*}{\overline{S}}_{n, m}^{0}$、$ \mathrm{*}{\overline{S}}_{n, m}^{1}$、$ \mathrm{*}{\overline{S}}_{n, m}^{2}$之间也应当满足第二类连带勒让德微分方程(the second kind of connected Legendre differential equation，LDE)，即存在微分方程关系式：

$ \begin{array}{c}({u}^{2}+{E}^{2})\mathrm{*}{\overline{S}}_{n, m}^{2}+2u\mathrm{*}{\overline{S}}_{n, m}^{1}-\\ \left[n\right(n+1)-\frac{{m}^{2}{E}^{2}}{{u}^{2}+{E}^{2}}]\mathrm{*}{\overline{S}}_{n, m}^{0}=0\end{array} $

(17)

3 数值计算与分析

图 1为传统Jekeli和修正Jekeli递推方法的收敛性比较。由图 1可知，当计算特殊情形，即在参考椭球边界面上$ u=b=6\mathrm{ }356\mathrm{ }752.314\mathrm{ }1 $ m，递推计算收敛精度要求为$ \varepsilon =1\times {10}^{-12}$时，随着计算阶数$ n$和次数$ m$的不断增大，传统Jekeli递推方法计算的$ {\overline{S}}_{n, m}^{0}$、$ {\overline{S}}_{n, m}^{1}$、$ {\overline{S}}_{n, m}^{2}$直到最高阶次$ n=m=2\mathrm{ }190$时，收敛项$ k=70$；同样，修正Jekeli递推方法计算的$ \mathrm{*}{\overline{S}}_{n, m}^{0}$、$ \mathrm{*}{\overline{S}}_{n, m}^{1}$、$ \mathrm{*}{\overline{S}}_{n, m}^{2}$直到最高阶次为$ n=m=2\mathrm{ }190$时，收敛项$ k=30 $。可以得到，在同样的递推计算收敛精度要求下，修正Jekeli递推方法比传统Jekeli递推方法计算需要的项数$k $减少一半，即需要的收敛项更少，收敛速度更快，可使得第二类连带勒让德函数的各阶偏导数计算效率提高。

图 2为传统Jekeli方法和修正Jekeli方法计算LDE方程的比较。由图 2可知，传统Jekeli递推方法计算的$ {\overline{S}}_{n, m}^{0}$、$ {\overline{S}}_{n, m}^{1} $、$ {\overline{S}}_{n, m}^{2} $满足第二类连带勒让德微分方程，随着计算阶数m和次数n的不断增大以及收敛项$k $增加时，收敛精度退化特别显著；而修正Jekeli递推方法计算的$ \mathrm{*}{\overline{S}}_{n, m}^{0}$、$ \mathrm{*}{\overline{S}}_{n, m}^{1} $、$ \mathrm{*}{\overline{S}}_{n, m}^{2} $也满足第二类连带勒让德微分方程（LDE），收敛精度始终保持在$ 1\times {10}^{-6}$左右。

图 3为修正Jekeli方法计算的各阶导数。由图 3可知，对于不同的高度h，即在参考椭球边界面$ u=b=6\mathrm{ }356\mathrm{ }752.314\mathrm{ }1$ m和$ u=b+255\mathrm{ }000 $ m，分别计算修正Jekeli递推方法的$ \mathrm{*}{\overline{S}}_{n, m}^{0}$、$ \mathrm{*}{\overline{S}}_{n, m}^{1} $、$ \mathrm{*}{\overline{S}}_{n, m}^{2}$值，可以发现，其阶数$n $与高度$ h$之间的关系与球近似下$ {\left(\frac{R}{r}\right)}^{n+1}$的阶数$ n$与高度$ h$之间的情况相似，且随着阶数$n $和高度$ h $增大，计算结果并没有出现溢出现象。

4 结语

本文推导出了第二类连带勒让德函数及其一阶、二阶导数的修正Jekeli递推计算方法，并与传统的Jekeli递推计算方法进行比较，结果表明，该递推方法计算收敛速度更快，收敛精度更加稳定。在具体计算过程中，推导了关于$ \mathrm{*}{\overline{S}}_{n, m}^{0}$、$ \mathrm{*}{\overline{S}}_{n, m}^{1}$、$ \mathrm{*}{\overline{S}}_{n, m}^{2} $的具体计算公式，对后面开展椭球域第二大地边值问题的研究具有很好的理论指导作用。