-
摘要: 基于$\frac{{{{\bar{P}}}_{nm}}\left( \cos \theta \right)}{\sin \theta }\left( m>0 \right)$的非奇异递推公式，给出了基于球坐标的引力矢量和垂线偏差非奇异计算公式；针对极点λ可任意取值引起的地方指北坐标系方向的不确定性问题，证明了引力矢量在转换到地心直角坐标系后不随λ的变化而变化，即与λ的取值无关。最终的数值计算结果表明，直角坐标系下的非奇异计算公式与本文提出的球坐标下的非奇异计算公式所得计算结果绝对值差异小于10-16m/s2，证明了该非奇异公式的正确性。最后总结了所有引力位球函数一阶导、二阶导非奇异性计算的一般思路。Abstract: When computing gravity vector and vertical deviation using spherical harmonic function, singular problem exists in the formulas expressed by spherical coordinates. This will cause some errors in gravity vector and vertical deflection data and influence their application. This paper aims at proposing an alternative method to solve this problem. Based on the non-singular expression of $ \frac{{{{\bar{P}}}_{nm}}\left( \cos \theta \right)}{\sin \theta } $(m > 0), the paper gives the non-singular formulas expressed by spherical coordinates for computing gravity vector and vertical deviation. At North and South Poles, the paper proves that even values of λ are arbitrary, the values of gravity vector are sole when the values are transferred to Earth fixed rectangular coordinate system. In order to show the validity of our method, the paper computes gravity vectors at points θ=0, λ=$ \frac{i}{360} $2π(i=0, 1...359) using the former 100 degrees and orders of EGM2008. The absolute differences between the computing results by our method and the non-singular formula expressed in Cartesian coordinates are smaller than 10-16 m/s2, which show the validity of our method. The non-singular expression based on spherical function derived by the paper can make full use of the high accurate algorithms of Legendre function, so the proposed method has better generality ability compared with the non-singular formula expressed by Cartesian coordinates. Finally the methods for non-singular computing of all the first or second derivations of gravity field potential are summarized. The method of this paper can also be directly applied to the non-singular calculation of the spherical harmonic model of magnetic field, and the basic idea is similar to that of this paper.
-
Keywords:
- gravity /
- vertical deviation /
- North and South Poles /
- non-singular /
- spherical coordinates
-
国际上一般用球函数来表示地球引力位，为了方便使用引力位模型，若已知地球引力位系数，可计算出与之相关的任意物理量，如扰动位、大地水准面、引力矢量、垂线偏差、引力梯度等。CHAMP(CHAllenging Minisatellite Payload)、GRACE(Gravity Recovery and Climate Experiment)和GOCE(Gravity field and steady-state Ocean Circulation Explorer)计划等也习惯如此表示。然而，利用引力位球谐模型来表示部分引力场元时，会出现因子$ \frac{1}{\sin \theta }\left( m>0 \right) $，使得许多计算在极区无法进行。极区的数据异常重要，但GOCE卫星由于被设计成太阳同步轨道，倾角为96.7°，因此在极区无观测值。解决方法之一就是用其他模型模拟出该区域的观测值，然后作迭代计算；而GRACE轨道为极地轨道, 在靠近极区时，传统的引力分量球谐表达公式将出现奇异而无法直接计算。因此研究非奇异的计算公式显得很有必要。
文献[1]给出了一组基于笛卡尔坐标的递推公式，避免了球坐标下计算的奇异性，在当前关于引力及引力梯度的各类计算中，使用比较广泛[2-5]。但由于地球引力位一般用球函数表示，且球函数计算的核心是勒让德函数的计算[6-9]，可利用的方法较多，若用直角坐标来表示，则可利用的方法较少，对递推计算的稳定性问题也少有研究成果。文献[10-13]均结合勒让德函数的表示方法讨论引力位及其导数的计算；文献[14-15]研究了引力矢量及其梯度的非奇异表示方法，但许多公式仍然比较复杂；文献[16]讨论了勒让德导数的非奇异计算；文献[17-19]讨论了引力梯度的非奇异算法，但未讨论其他引力场元的非奇异算法，如引力矢量、垂线偏差等。引力矢量的非奇异算法对极轨卫星在极地区域的定轨十分重要。测高卫星可得到全球海洋的垂线偏差数据，陆地区域的垂线偏差可由已有重力场模型计算得到，传统的极区垂线偏差计算公式由于奇异因子的存在使得计算精度不可靠。另外，极点λ可任意取值，地方指北坐标系的x、y轴方向不唯一，这对计算有何影响也鲜有文献进行讨论。
本文的目的是给出引力及垂线偏差的非奇异计算公式，并对极点λ可任意取值带来的坐标系方向不确定性问题进行分析。首先简要分析了奇异性产生的原因和特点，然后给出了$ \frac{{{{\bar{P}}}_{nm}}\left( \cos \theta \right)}{\sin \theta }\left( m>0 \right) $、$ \frac{{{{\bar{P}}}_{nm}}\left( \cos \theta \right)}{{{\sin }^{2}}\theta }\left( m>0 \right) $的非奇异计算公式；接着给出了引力和垂线偏差的非奇异计算公式，并分析了球坐标下极点计算结果的唯一性；最后通过数值计算验证了公式的正确性，并给出了相应结论。
1 $\frac{{{{\bar{P}}}_{nm}}\left( \cos \theta \right)}{\sin \theta }\left( m>0 \right) $、$ \frac{{{{\bar{P}}}_{nm}}\left( \cos \theta \right)}{{{\sin }^{2}}\theta }\left( m>1 \right) $的非奇异计算公式
地球的引力位可表示为[19]：
$$ \begin{array}{*{20}{c}} {v = \frac{{GM}}{r} + \frac{{GM}}{r}\sum\limits_{n = 2}^\infty {{{\left( {\frac{a}{r}} \right)}^n}\sum\limits_{m = 0}^n {\left( {{{\bar C}_{nm}}\cos m\lambda + } \right.} } }\\ {\left. {{{\bar S}_{nm}}\sin m\lambda } \right){{\bar P}_{nm}}\left( {\cos \theta } \right)} \end{array} $$ (1) 式中，(r, θ, λ)为计算点的球坐标；a为地球长半轴；$ \bar{C} $nm、$ \bar{S} $nm为引力位系数；$ \bar{P} $nm(cosθ)为n阶m次的归一化连带勒让德函数；GM为地心引力常数。
为表示方便，首先定义如下坐标系，即沿经线圈朝北为x轴，沿纬线圈朝西为y轴，沿径向朝向天顶为z轴，本文称之为地方指北坐标系。利用式(2)，不考虑0阶项(该项计算不存在奇异)，在地方指北坐标系下引力及垂线偏差为:
$$ \left\{ \begin{array}{l} {g_x} = - \frac{{GM}}{{{r^2}}}\sum\limits_{n = 2}^\infty {{{\left( {\frac{a}{r}} \right)}^n}\sum\limits_{m = 0}^n {\left[ {{C_{nm}}\cos m\lambda + {S_{nm}}\sin m\lambda } \right]\frac{{\partial {{\bar P}_{nm}}}}{{\partial \theta }}} } \\ {g_y} = - \frac{{GM}}{{{r^2}\sin }}\sum\limits_{n = 2}^\infty {{{\left( {\frac{a}{r}} \right)}^n}\sum\limits_{m = 0}^n {m\left[ { - {C_{nm}}\sin m\lambda + {S_{nm}}\cos m\lambda } \right]{{\bar P}_{nm}}\left( {\cos \theta } \right)} } \\ {g_z} = - \frac{{GM}}{{{r^2}}}\sum\limits_{n = 2}^\infty {\left( {n + 1} \right){{\left( {\frac{a}{r}} \right)}^n}\sum\limits_{m = 0}^n {\left[ {{C_{nm}}\cos m\lambda + {S_{nm}}\sin m\lambda } \right]{{\bar P}_{nm}}\left( {\cos \theta } \right)} } \end{array} \right. $$ (2) $$ \left\{ \begin{array}{l} \xi = \frac{a}{r}\sum\limits_{n = 2}^\infty {{{\left( {\frac{a}{r}} \right)}^n}\sum\limits_{m = 0}^n {\left[ {\Delta {{\bar C}_{nm}}\cos m\lambda + \Delta {{\bar S}_{nm}}\sin m\lambda } \right]\frac{{{\rm{d}}{{\bar P}_{nm}}\left( {\cos \theta } \right)}}{{{\rm{d}}\theta }}} } \\ \eta = \frac{a}{{r\sin \theta }}\sum\limits_{n = 2}^\infty {{{\left( {\frac{a}{r}} \right)}^n}\sum\limits_{m = 0}^n {m\left( {\Delta {{\bar C}_{nm}}\sin m\lambda - \Delta {{\bar S}_{nm}}\cos m\lambda } \right){{\bar P}_{nm}}\left( {\cos \theta } \right)} } \end{array} \right. $$ (3) 式(3)中, Δ$ \bar{C} $nm、Δ$ \bar{S} $nm为扰动位系数。根据文献[20]，$ \bar{P} $nm(cosθ)的计算公式不存在奇异，文献[16-17]给出了$ \frac{\text{d}{{{\bar{P}}}_{nm}}\left( \cos \theta \right)}{\text{d}\theta } $的非奇异计算公式，因此gx、gz和ξ的计算均不存在奇异性问题。然而，引力gy及垂线偏差的东西分量η因含有因子$ \frac{1}{\sin \theta } $，在极区难以直接计算。经简单分析可知，上述因子最终的形式可归结为$ \frac{{{{\bar{P}}}_{nm}}\left( \cos \theta \right)}{{{\sin }^{k}}\theta }\left( k=1, 2 \right) $。根据文献[19]，$ \frac{{{{\bar{P}}}_{nm}}\left( \cos \theta \right)}{\sin \theta } $、$ \frac{{{{\bar{P}}}_{nm}}\left( \cos \theta \right)}{{{\sin }^{2}}\theta } $分别在m=0，m=0和1时奇异，而$ \frac{{{{\bar{P}}}_{nm}}\left( \cos \theta \right)}{\sin \theta }\left( m>0 \right) $、$ \frac{{{{\bar{P}}}_{nm}}\left( \cos \theta \right)}{{{\sin }^{2}}\theta }\left( m>1 \right) $非奇异，有关计算公式为[19]：
$$ \frac{{{P_{nm}}}}{{\sin \theta }} = \sin \theta {P_{nm}} + \frac{{\cos \theta }}{{2m}}\left[ {{P_{n,m + 1}} + \left( {n + m} \right)\left( {n - m + 1} \right){P_{n,m - 1}}} \right],m > 0 $$ (4) $$ \frac{{{P_{nm}}\left( {\cos \theta } \right)}}{{{{\sin }^2}\theta }} = {P_{nm}} + \frac{{\cos \theta }}{2}\left[ {\frac{1}{{m\sin \theta }}{P_{n,m + 1}} + \left( {n + m} \right)\left( {n - m + 1} \right)\frac{1}{{m\sin \theta }}{P_{n,m - 1}}} \right],m > 1 $$ (5) 正规化后可得：
$$ \frac{{{{\bar P}_{nm}}}}{{\sin \theta }} =\\ \left\{ \begin{array}{l} \sin \theta {{\bar P}_{nm}} + \frac{{\cos \theta }}{{2m}}\left[ {\sqrt {\left( {n + m + 1} \right)\left( {n - m} \right)} {{\bar P}_{n,m + 1}} + \sqrt {2\left( {n + m} \right)\left( {n - m + 1} \right)} {{\bar P}_{n,m - 1}}} \right],m = 1 < n\\ \sin \theta {{\bar P}_{nm}} + \frac{{\cos \theta }}{{\sqrt n }}{{\bar P}_{n,m - 1}},m = 1 = n\\ \sin \theta {{\bar P}_{nm}} + \frac{{\cos \theta }}{{2m}}\left[ {\sqrt {\left( {n + m + 1} \right)\left( {n - m} \right)} {{\bar P}_{n,m + 1}} + \sqrt {\left( {n + m} \right)\left( {n - m + 1} \right)} {{\bar P}_{n,m - 1}}} \right],1 < m < n\\ \sin \theta {{\bar P}_{nm}} + \frac{{\cos \theta }}{{\sqrt {2n} }}{{\bar P}_{n,m - 1}},m = n \end{array} \right. $$ (6) $$ \frac{{{{\bar P}_{nm}}}}{{{{\sin }^2}\theta }} =\\ \left\{ \begin{array}{l} {{\bar P}_{nm}} + \frac{{\cos \theta }}{{2m}}\left[ {\sqrt {\left( {n + m + 1} \right)\left( {n - m} \right)} \frac{{{{\bar P}_{n,m + 1}}}}{{\sin \theta }} + \sqrt {\left( {n + m} \right)\left( {n - m + 1} \right)} \frac{{{{\bar P}_{n,m - 1}}}}{{\sin \theta }}} \right],1 < m < n\\ {{\bar P}_{nm}} + \frac{{\cos \theta }}{{\sqrt {2n} }}\frac{{{{\bar P}_{n,m - 1}}}}{{\sin \theta }},m = n \end{array} \right. $$ (7) 为便于后面的推导和计算，可令$ {{\overline{PS}}_{nm}}=\frac{{{{\bar{P}}}_{nm}}}{\sin \theta } $，$ \overline{PS}_{nm}^{2}=\frac{{{{\bar{P}}}_{nm}}}{{{\sin }^{2}}\theta } $。利用式(6)和式(7)即可实现$ {{\overline{PS}}_{nm}}\left( m>0 \right) $、$ \overline{PS}{{2}_{nm}}\left( m>1 \right) $的非奇异计算。
2 引力及垂线偏差的非奇异公式
2.1 gy和η的非奇异表达
式(2)可写为：
$$ \begin{array}{*{20}{c}} {{g_y} = - \frac{{GM}}{{{r^2}}}\sum\limits_{n = 2}^\infty {{{\left( {\frac{a}{r}} \right)}^n}\sum\limits_{m = 0}^n {\left( { - {C_{nm}}\sin m\lambda + } \right.} } }\\ {\left. {{S_{nm}}\cos m\lambda } \right)m{{\overline {PS} }_{nm}}} \end{array} $$ (8) 当m=0时，不需要计算$ {{\overline{PS}}_{nm}} $；而当m > 0时，可利用式(6)计算$ {{\overline{PS}}_{nm}} $，所以：
$$ \begin{array}{*{20}{c}} {{g_y} = - \frac{{GM}}{{{r^2}}}\sum\limits_{n = 2}^\infty {{{\left( {\frac{a}{r}} \right)}^n}\sum\limits_{m = 1}^n {\left( { - {C_{nm}}\sin m\lambda + } \right.} } }\\ {\left. {{S_{nm}}\cos m\lambda } \right)m{{\overline {PS} }_{nm}}} \end{array} $$ (9) 因此整个计算过程非奇异，计算稳定。同理，垂线偏差的东西分量也可类似计算：
$$ \begin{array}{*{20}{c}} {\eta = \frac{{GM}}{r}\sum\limits_{n = 2}^\infty {{{\left( {\frac{a}{r}} \right)}^n}\sum\limits_{m = 1}^n {\left( {\Delta {{\bar C}_{nm}}\sin m\lambda - } \right.} } }\\ {\left. {\Delta {{\bar S}_{nm}}\cos m\lambda } \right)m{{\overline {PS} }_{nm}}} \end{array} $$ (10) 虽然垂线偏差的分量习惯称为南北和东西分量，但事实上其划分的依据是经线圈和纬线圈的切线方向。尽管在正北极或南极，不存在东西方向，但经线圈的切线方向和径向方向依然存在，利用右手准则也可确定y轴方向。因此式(10)的计算直接代入坐标即可。若要进一步地将垂线偏差分解到地心地固坐标系下的x、y方向，可将垂线偏差当做二维矢量进行坐标变换。
2.2 垂线偏差在地心地固直角坐标系下的表示
以引力矢量为例，首先按式(2)和式(9)计算出地方指北坐标系下引力矢量，用脚标x、y、z表示，地心地固坐标系下用X、Y、Z表示，然后利用式(11)进行坐标转换：
$$ \left[ {\begin{array}{*{20}{c}} {{g_X}}\\ {{g_Y}}\\ {{g_Z}} \end{array}} \right] = \mathit{\boldsymbol{N}}\left[ {\begin{array}{*{20}{c}} {{g_x}}\\ {{g_y}}\\ {{g_z}} \end{array}} \right] $$ (11) 式中，N = Rz(π－λ)Ry(θ)，表示先绕y轴逆时针旋转θ角，后绕z轴逆时针旋转(π－λ)。具体表达为：
$$ {\mathit{\boldsymbol{R}}_z}\left( {\pi - \lambda } \right) = \left[ {\begin{array}{*{20}{c}} { - \cos \lambda }&{\sin \lambda }&0\\ { - \sin \lambda }&{ - \cos \lambda }&0\\ 0&0&1 \end{array}} \right] $$ (12) $$ {\mathit{\boldsymbol{R}}_y}\left( \theta \right) = \left[ {\begin{array}{*{20}{c}} {\cos \theta }&0&{ - \sin \theta }\\ 0&1&0\\ {\sin \theta }&1&{\cos \theta } \end{array}} \right] $$ (13) 2.3 极点引力计算
利用式(9)、式(10)及文献[19]可实现极区附近引力、垂线偏差及引力梯度的非奇异计算，但对于极点，由于λ可以取任意值，最终的计算结果是否会有不同是一个问题。很显然，极点作为独立的一点，其值必是唯一的。为此本节将通过理论推导以验证式(9)、式(10)等在极点计算的可行性。
对于正极点，以北极为例，θ=0，有：
$$ \left[ {\begin{array}{*{20}{c}} {{g_X}}\\ {{g_Y}}\\ {{g_Z}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} { - \cos \lambda }&{\sin \lambda }&0\\ { - \sin \lambda }&{ - \cos \lambda }&0\\ 0&0&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{g_x}}\\ {{g_y}}\\ {{g_z}} \end{array}} \right] $$ (14) 此时，x方向依然定义为经线圈切线方向，z为地球径向方向，y轴与x、z轴符合右手法则。
根据勒让德函数的性质，易知θ=0时，当且仅当m=0，$ \bar{P} $nm≠0时，有：
$$ \frac{{\partial {{\bar P}_{nm}}}}{{\partial \theta }} = \left\{ \begin{array}{l} \sqrt {\frac{{n\left( {n + 1} \right)}}{2}} {{\bar P}_{n,0}},m = 1\\ 0,其他 \end{array} \right., $$ $$ {\overline {PS} _{nm}} = \left\{ \begin{array}{l} \sqrt {\frac{{n\left( {n + 1} \right)}}{2}} {{\bar P}_{n,0}},m = 1\\ 0,m > 1 \end{array} \right. $$ 所以有：
$$ \left\{ \begin{array}{l} {g_x} = \sum\limits_{n = 2}^\infty {\left( {{C_{n1}}\cos \lambda + {S_{n1}}\sin \lambda } \right){b_n}} \\ {g_y} = \sum\limits_{n = 2}^\infty {\left( { - {C_{n1}}\sin \lambda + {S_{n1}}\cos \lambda } \right){b_n}} \\ {g_z} = \sum\limits_{n = 2}^\infty {{C_{n0}}{c_n}} \end{array} \right. $$ (15) 式中，${{b}_{n}}=-\frac{GM}{{{r}^{2}}}{{\left( \frac{a}{r} \right)}^{n}}\sqrt{\frac{n\left( n+1 \right)}{2}}{{\bar{P}}_{n, 0}} $；$ {{c}_{n}}=-\frac{GM}{{{r}^{2}}}{{\left( \frac{a}{r} \right)}^{n}}\left( n+1 \right){{\bar{P}}_{n, 0}} $。进一步地，有：
$$ {g_x} = \left[ {\begin{array}{*{20}{c}} {\cos \lambda }&{\sin \lambda } \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\sum\limits_{n = 2}^\infty {{C_{n1}}{b_n}} }\\ {\sum\limits_{n = 2}^\infty {{S_{n1}}{b_n}} } \end{array}} \right] $$ (16) $$ {g_y} = \left[ {\begin{array}{*{20}{c}} { - \sin \lambda }&{\cos \lambda } \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\sum\limits_{n = 2}^\infty {{C_{n1}}{b_n}} }\\ {\sum\limits_{n = 2}^\infty {{S_{n1}}{b_n}} } \end{array}} \right] $$ (17) 结合式(14)可知：
$$ \begin{array}{*{20}{c}} {\left[ {\begin{array}{*{20}{c}} {{g_X}}\\ {{g_Y}}\\ {{g_Z}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} { - \cos \lambda }&{\sin \lambda }&0\\ { - \sin \lambda }&{ - \cos \lambda }&0\\ 0&0&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\cos \lambda }&{\sin \lambda }&0\\ { - \sin \lambda }&{\cos \lambda }&0\\ 0&0&1 \end{array}} \right]}\\ {\left[ {\begin{array}{*{20}{c}} {\sum\limits_{n = 2}^\infty {{C_{n1}}{b_n}} }\\ {\sum\limits_{n = 2}^\infty {{S_{n1}}{b_n}} }\\ {\sum\limits_{n = 2}^\infty {{C_{n0}}{c_n}} } \end{array}} \right]} \end{array} $$ (18) 最终得：
$$ \left[ {\begin{array}{*{20}{c}} {{g_X}}\\ {{g_Y}}\\ {{g_Z}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} { - 1}&0&0\\ 0&{ - 1}&0\\ 0&0&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\sum\limits_{n = 2}^\infty {{C_{n1}}{b_n}} }\\ {\sum\limits_{n = 2}^\infty {{S_{n1}}{b_n}} }\\ {\sum\limits_{n = 2}^\infty {{C_{n0}}{c_n}} } \end{array}} \right] $$ (19) 由此可知，对于正极点，虽然λ可取任意值，但经坐标转换到地心直角坐标系时，其值唯一，不随经度改变而改变。此处λ取值仅影响地方指北坐标系方向的定义，对地心直角坐标系下的分量大小并无影响。事实上垂线偏差的计算也如此。虽然北极点无东西方向的划分，只要λ值一定，均可定义出相应的方向，并计算出相应的分量，同样可类似式(18)经坐标变换得到直角坐标系下X、Y方向的分量。由于公式形式相似，此处也不再作单独推导。
不同于非两极点，尽管极点地方指北坐标系的方向随着λ的不同而有所差别，但通过上面的推导发现，极点的计算值转换到直角坐标系后，值是唯一的，与λ的取值无关。这不仅意味着§2.1、§2.2所给公式可用于极点的计算，也再次表明极点计算的奇异性是假奇异。
3 算例
为了验证论文非奇异公式的正确性，特别是极点附近计算的稳定性，现对引力矢量采用如下方法进行计算。
方法1 利用式(9)计算gy，其余计算同式(2)，然后按式(11)进行坐标变换；
方法2 采用直角坐标递推公式进行计算[21]。
引力位模型采用EGM08(Earth Gravitational Model 2008)前2~100阶完整阶次。为了验证计算效果，现利用方法1和方法2分别计算θ=0，λ=$ \frac{i}{360} $2π(i=0, 1，…，359)所有的引力值，统计结果如图 1所示。
在图 1中横轴代表角度，变化范围为0°~359°，纵轴代表两方法计算结果的差异。由图 1可知，利用本文所给出的基于球坐标的非奇异计算公式(方法1)和文献[21]所给出的直角坐标系下的非奇异计算公式(方法2)所得计算结果的绝对值差异小于10-16m/s2，该差异主要由计算机浮点数的精度所导致。这表明本文给出的公式能够解决极区计算的奇异性问题。值得说明的是，此试验的计算点在直角坐标系下为同一点，因此这组算例也验证了§2.3中结论，即采用球坐标表示时，尽管极点的λ取值不唯一，但所对应的地心直角坐标系下的引力矢量却是唯一的。
由于垂线偏差的计算形式和引力矢量类似，因此此节不再安排专门算例。
4 结语
本文可看作文献[19]的后续工作，基本思路均是首先计算出$ \frac{{{{\bar{P}}}_{nm}}\left( \cos \theta \right)}{\sin \theta }\left( m>0 \right) $、$ \frac{{{{\bar{P}}}_{nm}}\left( \cos \theta \right)}{{{\sin }^{2}}\theta }\left( m>1 \right) $，然后将各引力场元中的奇异因子消掉，现对相应成果作如下总结。
1) $ \frac{{{{\bar{P}}}_{nm}}\left( \cos \theta \right)}{\sin \theta } $、$ \frac{{{{\bar{P}}}_{nm}}\left( \cos \theta \right)}{{{\sin }^{2}}\theta } $分别在m=0、m=0和1时奇异，$ \frac{{{{\bar{P}}}_{nm}}\left( \cos \theta \right)}{\sin \theta }\left( m>0 \right) $、$ \frac{{{{\bar{P}}}_{nm}}\left( \cos \theta \right)}{{{\sin }^{2}}\theta }\left( m>1 \right) $的非奇异计算公式可用式(6)和式(7)完成。
2) 在计算与引力位的一阶导数有关的量时，如引力矢量、垂线偏差等，尽管$ \frac{{{{\bar{P}}}_{nm}}\left( \cos \theta \right)}{\sin \theta } $在m=0时奇异，但实际需要计算的因子是$ \frac{m{{{\bar{P}}}_{nm}}\left( \cos \theta \right)}{\sin \theta } $，因此$ \frac{{{{\bar{P}}}_{nm}}\left( \cos \theta \right)}{\sin \theta }\left( m=0 \right) $无需计算；在计算与引力位的二阶导数有关的量时(如引力梯度)，会出现$ \frac{{{{\bar{P}}}_{nm}}\left( \cos \theta \right)}{{{\sin }^{2}}\theta } $，尽管其在m=0和1时奇异，但需要计算的因子是$ \frac{m\left( m-1 \right){{{\bar{P}}}_{nm}}\left( \cos \theta \right)}{{{\sin }^{2}}\theta } $(文献[19])，因此$ \frac{{{{\bar{P}}}_{nm}}\left( \cos \theta \right)}{{{\sin }^{2}}\theta }\left( m=0、1 \right) $也无须计算，这是球坐标表达式非奇异性可以消除的根本原因。
3) 本文推导的基于球函数的非奇异表达可充分利用勒让德函数的有关研究成果，相对于直角坐标表达下的非奇异表达通用性更好。
4) 尽管在极点λ可取任意值，但其主要影响地方指北坐标系方向的定义，在利用球函数计算出引力矢量并转换到地心直角坐标系后，值与λ无关。
本文的研究成果也可以直接用于磁场球谐模型的非奇异计算，基本思路和本文类似，即首先计算出$ \frac{{{{\bar{P}}}_{nm}}\left( \cos \theta \right)}{\sin \theta }\left( m>0 \right) $、$ \frac{{{{\bar{P}}}_{nm}}\left( \cos \theta \right)}{{{\sin }^{2}}\theta }\left( m>1 \right) $，然后代入磁场分量的计算公式中，进行非奇异计算。
-
[1] Cunningham L. On the Computation of the Spherical Harmonic Terms Needed During the Numerical Integration of the Orbital Motion of an Artificial Satellite[J]. Celestial Mechanics, 1970, 2:207-216 doi: 10.1007/BF01229495
[2] 王正涛. 卫星跟踪卫星测量确定地球重力场的理论与方法[D]. 武汉: 武汉大学, 2005 http://cdmd.cnki.com.cn/Article/CDMD-10486-2006031327.htm Wang Zhengtao. Theory and Methodology of Earth Gravity Field Recovery by Satellite-to-Satellite Tracking Data[D]. Wuhan:Wuhan University, 2005 http://cdmd.cnki.com.cn/Article/CDMD-10486-2006031327.htm
[3] 邹贤才. 卫星轨道理论与地球重力场模型的确定[D]. 武汉: 武汉大学, 2007 Zou Xiancai. Theory of Satellite Orbit and Earth Gravity Field Determination[D]. Wuhan:Wuhan University, 2005
[4] 钟波. 基于GOCE卫星重力测量技术确定地球重力场的研究[D]. 武汉: 武汉大学, 2010 http://cdmd.cnki.com.cn/Article/CDMD-10486-2010167102.htm Zhong Bo. Study on the Determination of the Earth's Gravity Field from Satellite Gravimetry Mission GOCE[D]. Wuhan:Wuhan University, 2010 http://cdmd.cnki.com.cn/Article/CDMD-10486-2010167102.htm
[5] Yi Weiyong. An Alternative Computation of a Gravity Field Model from GOCE[J]. Advances in Space Research, 2012, 50(3):371-384 doi: 10.1016/j.asr.2012.04.018
[6] Holmes S, Featherstone W. A Unified Approach to the Clenshaw Summation and the Recursive Computation of very High Degree and Order Normalised Associated Legendre Functions[J]. Journal of Geodesy, 2002, 76(5):279-299 doi: 10.1007/s00190-002-0216-2
[7] Yu Jinhai, Wan Xiaoyun, Zeng Yanyan. The Integral Formulas for Legendre Functions and Associated Legendre Functions[J]. Journal of Geodesy, 2011, DOI: 10.1007/s00190-011-0529-0
[8] 吴星, 刘雁雨.多种超高阶次缔合勒让德函数计算方法的比较[J].测绘科学技术学报, 2006, 23(3):188-191 http://www.docin.com/p-885217372.html Wu Xing, Liu Yanyu. Comparison of Computing Methods of the Ultra-High Degree and Order[J]. Journal of Geomatics Science and Technology, 2006, 23(3):188-191 http://www.docin.com/p-885217372.html
[9] 魏子卿.完全正常化缔合勒让德函数及其导数与积分的递推关系[J].武汉大学学报·信息科学版, 2016, 41(1):27-36 http://ch.whu.edu.cn/CN/abstract/abstract3430.shtml Wei Ziqing. Recurrence Relations for Fully Normalized Associated Legendre Functions and Their Derivatives and Integrals[J]. Geomatics and Information Science of Wuhan University, 2016, 41(1):27-36 http://ch.whu.edu.cn/CN/abstract/abstract3430.shtml
[10] Casotto S, Fantino E. Evaluation of Methods for Spherical Harmonic Synthesis of the Gravitational Potential and its Gradients[J]. Advances in Space Research, 2007, 40(1):69-75 doi: 10.1016/j.asr.2007.01.021
[11] Fantino E, Casotto S. Methods of Harmonic Synthesis for Global Geopotential Models and Their First-, Second-and Third-Order Gradients[J]. Journal of Geodesy. 2009, 83(7):595-619 doi: 10.1007/s00190-008-0275-0
[12] Casotto S, Fantino E. Gravitational Gradients by Tensor Analysis with Application to Spherical Coordinates[J]. Journal of Geodesy, 2009, 83(7):621-634 doi: 10.1007/s00190-008-0276-z
[13] Tscherning C. Computation of Second-Order Deriv-atives of the Normal Potential Based on the Repre-sentation by a Legendre Series[J]. Manuscr Geod, 1976, 1:71-92 http://www.researchgate.net/publication/259055140_Computation_of_second-order_derivatives_of_the_normal_potential_based_on_a_representation_by_a_Legendre_series
[14] Balmino G, Barrriot J, Vales N. Non-singular Formulation of the Gravity Vector and Gravity Gradient Tensor in Spherical Harmonics[J]. Manuscr Geod, 1989, 15(1):11-16 http://www.researchgate.net/publication/241346652_Nonsingular_formulation_of_the_gravity_vector_and_gravity_gradient_tensor_in_spherical_harmonics
[15] Petrovskaya M, Vershkov A. Non-singular Expressions for the Gravity Gradients in the Local North-Oriented and Orbital Reference Frames[J]. Journal of Geodesy, 2006, 80(3):117-127 doi: 10.1007/s00190-006-0031-2
[16] 于锦海, 万晓云.计算Legendre函数导数的非奇异方法[J].测绘科学技术学报, 2010, 27(1):1-3 http://kns.cnki.net/KCMS/detail/detail.aspx?filename=jfjc201001004&dbname=CJFD&dbcode=CJFQ Yu Jinhai, Wan Xiaoyun. Non-singular Formulae for Computing Derivatives of Legendre Functions[J]. Journal of Geomatics Science and Technology, 2010, 27(1):1-3 http://kns.cnki.net/KCMS/detail/detail.aspx?filename=jfjc201001004&dbname=CJFD&dbcode=CJFQ
[17] 刘晓刚, 吴晓平, 赵东明, 等.扰动重力梯度的非奇异表示[J].测绘学报, 2010, 39(5):450-457 http://www.cnki.com.cn/Article/CJFDTotal-CHXB201005005.htm Liu Xiaogang, Wu Xiaoping, Zhao Dongming, et al. Non-singular Expression of the Disturbing Gravity Gradients[J]. Acta Geodaetica et Cartographica Sinica, 2010, 39(5):450-457 http://www.cnki.com.cn/Article/CJFDTotal-CHXB201005005.htm
[18] 刘晓刚. GOCE卫星测量恢复地球重力场模型的理论与方法[D]. 郑州: 信息工程大学, 2011 http://cdmd.cnki.com.cn/article/cdmd-90008-1012325190.htm Liu Xiaogang. Theory and Methods of the Earth's Gravity Field Model Recovery from GOCE Data[D]. Zhengzhou:Information Engineering University, 2011 http://cdmd.cnki.com.cn/article/cdmd-90008-1012325190.htm
[19] 万晓云.引力场梯度张量的非奇异公式推导[J].武汉大学学报·信息科学版, 2011, 36(12):1486-1489 http://ch.whu.edu.cn/CN/abstract/abstract732.shtml Wan Xiaoyun. New Derivation of Non-singular Expression for Gravitational Gradients Calculation[J]. Geomatics and Information Science of Wuhan University, 2011, 36(12):1486-1489 http://ch.whu.edu.cn/CN/abstract/abstract732.shtml
[20] 海斯卡涅, 莫里兹. 卢福康, 胡国理, 译. 物理大地测量学[M]北京; 测绘出版社, 1979 Heiskanen W A, Moritz H. Lu Fukang, Hu Guoli, Trans.Physical Geodesy[M]. Beijing:Surveying Press, 1979
[21] Oliver M, Eberhard G. Satellite Orbits:Models, Methods, and Applications[M]. Berlin:Springer, 2001
-
期刊类型引用(0)
其他类型引用(2)
下载:
计量
- 文章访问数: 1772
- HTML全文浏览量: 190
- PDF下载量: 301
- 被引次数: 2
