Recurrence Relations for Fully Normalized Associated Legendre Functions and Their Derivatives and Integrals
-
摘要: 在地球重力场问题中,常用到完全正常化缔合勒让德函数及其导数、积分的递推关系。当前流行的地球扰动位模型均采用完全正常化的缔合勒让德函数,用此类模型可以高效方便计算各种扰动重力场元。随着本世纪多个新一代卫星重力探测计划成功实施,高阶或超高阶地球重力场模型的研究备受学界的关注。有关完全正常化缔合勒让德函数的递推关系对于高阶重力场模型具有特别意义。本文在前人研究的基础上,用初等微积分导出了若干新的递推关系式。同时还推导了正常化缔合勒让德函数及其导数、积分的检核式,这些检核式涉及地球位的球谐级数的数学性质。
-
关键词:
- 完全正常化缔合勒让德函数 /
- 导数 /
- 积分 /
- 递推关系 /
- 检核式
Abstract: Recurrence relations for fully normalized associated Legendre functions and their derivatives and integrals are often used when studying the Earth gravity field. Fully normalized associated Legendre functions are all adopted in current popular Earth disturbance potential models make the calculations of various disturbance gravity field elements highly efficient and easy when such models and related recurrence relations are used. New generation satellite gravity exploration missions have been successfully implemented, and high or ultrahigh degree Earth gravity field models have been quickly developed. This study area has received much attention in academic community; our investigation has a special meaning as it addresses the recurrence relations of fully normalized associated Legendre functions. Building upon previous research, several new recurrence relation expressions are derived in detail based on elementary calculus as the mathematical tool for the derivation. Formulas for assessing the values of fully normalized associated Legendre functions, their derivatives, and integrals are also deduced. This study involves the mathematical properties of the spherical harmonic expansion series of the geo-potential, thus is related to fundamental theory research. -
北斗卫星导航系统(BeiDou navigation satellite system，BDS)是中国自主研制、独立运行的卫星导航系统，在经历了“北斗一代”到“北斗二代”的发展后，现已基本实现亚太地区的导航、定位、授时等服务，并计划在2020年完成全球系统的建设。现阶段可以正常提供导航定位服务的13颗北斗卫星分别是5颗同步地球轨道(geostationary orbits，GEO)卫星、5颗中倾斜地球同步轨道(inclined geo-synchronous orbits，IGSO)卫星和3颗中圆地球轨道(medium earth orbit，MEO)卫星，经研究发现这13颗卫星的伪距观测值均存在明显的系统偏差，因此针对北斗伪距偏差的分析和改正对提高北斗卫星导航定位精度和应用性能有着重要意义。
自从Hauschild等提出北斗伪距偏差的存在后，国内外学者进行了广泛的研究，Montenbruck等发现了伪距偏差与高度角和频率相关的特性，由于不同的测站和不同的接收机类型计算得到的伪距偏差对于同一颗卫星具有一致性，因此认为伪距偏差的误差源来自于北斗卫星端的多路径[1-5]。张小红等通过对北斗卫星非差、单差和双差宽巷组合(melbourne-wübbena，MW)的系统偏差进行研究进一步确定了伪距偏差的产生归因于卫星端的多路径，其存在对宽巷模糊度的固定有重大的影响[6-9]。为了削弱北斗卫星伪距偏差的影响，Wanninger等根据伪距偏差与高度角相关的特性建立了北斗IGSO和MEO卫星的伪距改正模型，该模型的建立对BDS的单频定位精度有较为明显的提高[10]。但是该模型仅采用了2014年3月10~20日、7月10~20日共20 d的数据进行建模，没有具体分析测站因素和时间因素对建模的影响，建模数据量较少不足以得到适用于全球的高精度模型。
本文通过分析测站数目及分布，以及观测时长对建模的影响后，采用全球18个测站对MEO卫星的伪距偏差进行拟合，利用其中4个可以观测到较为完整的IGSO卫星弧段的测站对IGSO卫星进行拟合，以高度角每隔5°为间隔进行加权分段线性拟合联合抗差估计得到改进后的伪距改正模型，并测试其对不同测站的精密单点定位(precise point positioning，PPP)精度提高情况，证明了模型的有效性。
1 北斗卫星伪距偏差分析
多路径组合(multipath combination，MP)由单频伪距和双频相关观测值组合得到，具体表示为[11]：
$$ \begin{array}{l} {M_{Pi}} = {P_i} - \frac{{f_i^2 + f_j^2}}{{f_i^2 - f_j^2}}{L_i} + \frac{{2f_j^2}}{{f_i^2 - f_j^2}}{L_j}\\ = {M_i} - \frac{{f_i^2 + f_j^2}}{{f_i^2 - f_j^2}}{m_i} + \frac{{2f_j^2}}{{f_i^2 - f_j^2}}{m_j} + {B_i} + {\varepsilon _i} \end{array} $$ (1) 式中，i、j(i≠j)表示不同的信号频率；MPi是伪距多路径组合；P、L是伪距观测值和相位观测值；f对应信号的频率；M是伪距观测值的多路径误差；而m是相位观测值的多路径影响；B包含了相位的模糊度和信号延迟信息，如果在连续观测且没有发生周跳的情况下B保持常数特性；ε是组合观测值的噪声。
由式(1) 可知，MP组合消除了卫星到地球的几何距离、钟差、各项误差改正和大气延迟改正，仅残余模糊度和通道时延、多路径误差和噪声相关。模糊度和通道时延当作常数值平滑后，MP组合仅与多路径误差和噪声有关[12-13]。
选取CUT0站2015年3月1日的观测数据计算得到北斗卫星的MP序列与高度角序列，图 1给出了不同频率下北斗GEO、IGSO、MEO卫星的MP时间序列以及高度角序列，可以看出IGSO和MEO卫星的MP序列中存在随着高度角的增大而减小的系统误差，即伪距偏差。GEO卫星受到伪距偏差的影响在一天的范围内MP序列呈现波动现象，但由于轨道高度角的变化范围较小，伪距偏差的大小与高度角的关系并不显著。图中不同频率对应的MP序列存在着比较明显的区别，说明伪距偏差的大小与频率相关。下文将主要根据IGSO和MEO卫星伪距偏差与高度角相关的特性，针对不同的频率进行模型估计。
![]()
图 1 不同频率上的北斗卫星多路径误差序列与高度角序列Figure 1. BDS MP Series and Elevation Series on B1/B2/B3 Frequency2 北斗IGSO和MEO卫星的伪距偏差模型分析
截止2015年12月，全球定位系统实验中心(MultiGNSS Experiment，MGEX)提供了全球可以接收北斗信号的约70个测站的观测信息，大多数测站使用的Trimble Net R9接收机可以提供三频观测值(B1/B2/B3)，部分测站的接收机如Septentrio POLARX4，仅提供北斗双频观数据(B1/B2)。本文选择了18个观测到较多北斗卫星的测站不同时间长度的数据作为采样数据进行模型分析和建模处理，数据的采样间隔是30 s，在数据处理之前已对观测数据中的周跳进行了探测和修复。
2.1 伪距偏差模型的估计方法
根据北斗IGSO和MEO卫星的伪距偏差与高度角大小相关的特性，文章选择加权分段线性拟合联合抗差估计的方法进行高度角建模，具体方法如下。
首先将MP序列按高度角5°间隔划分弧段，高度角范围是5°~85°，每一个弧段进行最小二乘线性拟合，分段的节点具备连续性，对应的目标函数S为:
$$ \left\{ \begin{array}{l} S = \sum\limits_{j = 1}^{m - 1} {\sum\limits_{i = 1}^{{n_j}} {{P_{j, i}}} {{\left( {{f_{j, i}} - {M_{j, i}}} \right)}^2} = \min } \\ {f_j}\left( {{e_j}} \right) - {f_{j + 1}}\left( {{e_j}} \right) = 0 \end{array} \right. $$ (2) 式中，j=1, 2, …，m-1，其中m是分段节点总数；i=1, 2, …，n，n是对应弧段上MP观测值总数；fj, i、Mj, i和Pj, i分别对应于第j弧段上第i个点的拟合函数值、MP观测值和权重；ej是第j个节点上的高度角；fj、fj+1是第j和j+1弧段上与高度角相关的拟合函数。由于MP序列中包含了测站多路径和伪距偏差，在低高度角情况下测站的多路径尤其严重，考虑到不同高度角下MP观测值精度不同，在分段拟合过程中对不同的高度角对应的MP观测值设置了不同的权重。文章选择的权重设定方案具体为：
$$ {P_{j, i}} = {a^2} + {a^2}/{\sin ^2}\left( {{e_{j, i}}} \right) $$ (3) 式中，a是指北斗观测值的伪距精度，通常取0.3 m；ej, i指第j弧段上第i个点对应的高度角。
另一方面，最小二乘估计后的估值精度R的计算式为：
$$ R = \sqrt {\frac{{\sum\limits_{j = 1}^{m - 1} {\sum\limits_{i = 1}^{{n_j}} {{{\left( {{f_{j, i}} - {M_{j, i}}} \right)}^2}} } }}{{n - 1}}} $$ (4) 式中，n表示样本的总量，其他符号的含义与式(2) 相同。利用式(2)、式(3)、式(4) 可以计算出每一颗IGSO和MEO卫星的伪距偏差模型和模型精度，同类卫星最终的模型精度通过所有卫星取均值得到。选择CUT0测站2015年3月1~10日的观测数据进行拟合建模和精度评估，IGSO和MEO卫星的伪距偏差的模型精度分别是0.35 m和0.41 m，相比于北斗伪距的常用精度信息0.3 m精度偏低, 因此为了控制北斗伪距的精度和抵抗粗差的影响，本文使用了抗差估计进行质量控制，所使用的抗差方法是IGGⅢ方案[14]。
2.2 测站数目和测站分布对模型估计的影响
通常每个测站观测到同一颗卫星的高度角范围并不相同，这与北斗IGSO卫星和MEO卫星的运行轨迹密切相关。IGSO卫星的星下点轨迹主要分布于亚太区域，在此范围的少数测站可以观测到高度角范围为5°~85°的IGSO卫星，大部分测站观测到IGSO卫星的高度角范围较小，在此范围外的测站基本观测不到IGSO卫星。而大多数测站都可以观测到MEO卫星完整的观测弧段，且每个测站的观测到的卫星高度角范围仍有差异。本文统计了18个测站可以观测到的北斗IGSO和MEO卫星最大的高度角范围后发现，对于IGSO卫星，只有CUT0、GMSD、JFNG和ADE2测站能够观测到基本完整的弧段，所有测站均可观测到较大高度角范围的MEO卫星。为了验证测站数目和分布对模型估计的影响，分别用4个测站(CUT0、GMSD、JFNG和ADE2)、12个测站、18个测站在2015年7月1~10日共10 d的数据对北斗卫星建模，然后测试不同改正模型对SIN1测站(SIN测站没有参与建模)MP序列的精度的提高情况。图 2给出了B2频率上的IGSO卫星的建模结果和精度提高情况，从图中可以看出，对于IGSO卫星，用4个测站建立的高度角模型与用12、18个测站计算的结果有较为明显的差异，而12个测站和18个测站的建模结果整体吻合，仅在某些节点处有所差异；从精度提升情况来看，相比12、18个测站的建模结果，4个测站的建模结果更好，主要是因为如果不能在一天的时间范围内连续观测到较为完整的卫星弧段，MP序列中包含的测站多路径和噪声的均值不为0，在进行去均值的平滑处理时会导致MP值的大小产生较大的偏差，从而降低模型估计的精度。图中的4个测站全部可以观测到全弧段的IGSO卫星，12个测站和18个测站中由于加入了大量没有全弧段观测的数据，降低了模型的精确度，因此，测站的分布也会影响到IGSO卫星的模型估计。对于MEO卫星，18个测站建立的改正模型的改正效果优于4个测站和12个测站的结果，主要是因为不同测站对应的高度角范围不同，测站数目太少会导致部分高度角范围缺少相应的观测值，因此测站数目越多, 外符合精度越高。
![]()
图 2 不同测站数目和分布建立的模型对SIN1站MP2序列的精度提高情况Figure 2. The Improved Value of MP2 Series 's Accuracy in SIN1 Station by Different Model Considering the Number and Distribution of Stations2.3 观测时间对模型估计的影响
考虑到测站数目以及测站分布对模型估计的影响，选择的测站基本可以观测到北斗IGSO和MEO的完整弧段，但是受到卫星轨道的影响，单天可以观测到的低高度角和高高度角下的观测值较少，为了能够更加准确的估计两个极端高度角时的伪距偏差，因此需要考虑观测时间对模型估计的影响。图 3给出了不同时段观测数据1个月(2015年3月)、6个月(2015年1~6月)、12个月(2015年1~12月)的B1频率上的建模结果对比，以及对伪距偏差的改正程度。图 3(a)是4个测站IGSO卫星的建模结果，可以看出仅用1个月的数据建模后的估计值波动较大，采用了6个月的数据和12个月的数据后，模型不再展现出强烈的波动性，模型改正值也基本接近。图 3(b)是不同模型对JFNG测站MP序列精度的提高情况，可以看出采用6个月的数据相对于仅采用1个月的数据建模对MP序列的精度提高了1~2 cm，选择12个月的数据建模结果进一步提高了MP序列的精度，但是提高幅度已不明显，接近饱和状态。说明适当增加观测数据的时间可以削弱随机误差的影响，提高模型的估计精度。
![]()
图 3 根据不同时间的数据建立的估计模型及对应模型改正对MP1序列的精度提高情况Figure 3. The Improved Value of MP1 Series's Accuracy by Different Model Considering the Observation Time2.4 建模结果及模型精度
在分析了测站数目与分布以及观测时间对模型建立的影响后，本文选择了18个测站2015全年共12个月的数据对北斗MEO卫星的伪距偏差建模，其中4个测站(CUT0、GMSD、JFNG和ADE2) 的数据对北斗IGSO卫星进行建模，建模的方法选择加权分段线性拟合，得到北斗IGSO卫星和MEO卫星在高度角5°~85°范围内以5°为间隔的伪距改正值, 具体结果如表 1。由于每一段的伪距偏差与高度角具有线性关系，可以内插得到任意高度角下的伪距改正值。图 4给出了CUT0、WHHN、KHAJ、DRAG测站改正前后IGSO卫星和MEO卫星B3频率上MP序列的精度信息对比，其中WHHN、KHAJ、DRAG 3个测站没有参与建模，可以看出，伪距偏差的改正可以修正伪距的系统误差，提高了伪距的精度，将北斗伪距的精度基本控制在0.3 m左右，B1频率和B2频率对应的MP序列结果与之相符，说明模型的改正的是合理的。
表 1 改进模型对应的节点伪距改正值Table 1. Piecewise Linear Correction Values for BeiDou Code Measurements for Improved Model高度角/(°) 改正值/m IGSO卫星 MEO卫星 B1 B2 B3 B1 B2 B3 5 -0.238 -0.246 -0.406 -0.375 -0.400 -0.134 10 -0.540 -0.421 -0.211 -0.405 -0.238 -0.174 15 -0.430 -0.296 -0.169 -0.314 -0.213 -0.112 20 -0.274 -0.272 -0.157 -0.220 -0.228 -0.131 25 -0.255 -0.260 -0.138 -0.180 -0.131 -0.047 30 -0.265 -0.227 -0.253 -0.144 -0.113 -0.052 35 -0.168 -0.154 -0.089 -0.179 -0.114 -0.081 40 -0.137 -0.115 -0.118 -0.091 -0.065 -0.030 45 -0.064 -0.068 -0.073 0.005 0.006 0.044 50 -0.019 -0.041 -0.010 0.081 0.068 0.048 55 0.025 0.037 -0.008 0.222 0.187 0.083 60 0.130 0.079 0.049 0.322 0.224 0.133 65 0.175 0.116 0.084 0.463 0.326 0.201 70 0.238 0.167 0.127 0.631 0.440 0.287 75 0.234 0.202 0.152 0.716 0.485 0.288 80 0.272 0.250 0.207 0.918 0.583 0.367 85 0.302 0.259 0.155 0.955 0.628 0.393 ![]()
图 4 不同测站伪距偏差改正前后的MP3序列精度Figure 4. The Accuracy of MP3 Series of Different Stations After Correction3 北斗卫星伪距偏差改正对精密定位的影响
通过实验分析北斗卫星伪距偏差改正对精密单点定位的影响，验证改进模型的有效性并与Wanninger等建立的伪距偏差改正模型(传统模型)进行比较，采用3种策略处理北斗PPP：(1) 不进行伪距偏差改正，直接进行PPP处理；(2) 利用传统模型改正伪距偏差后进行PPP处理；(3) 利用改进模型改正伪距偏差后进行PPP处理。
选择了5个可以观测到较多北斗卫星的测站(CUT0/GMSD/JFNG/SIN1/WARK)2014年10月15~20日和2015年3月1~5日共10 d的数据进行PPP解算，表 2统计了3种PPP处理策略后的平均定位精度，结果表明改进模型相对于传统模型的定位精度在E、N、U方向有8.4%、10.1%、8.5%的提高，相比于未加伪距改正进行PPP处理的定位精度，E、N、U方向平均提高了17.4%、19.1%、9.5%。由于PPP主要依赖高精度载波相位观测值的作用，伪距观测值的改正对提高PPP的精度能力有限，PPP解算初始阶段模糊度还未收敛，伪距观测值在此时起决定性作用，改正伪距偏差可以有效得提高PPP的收敛速度，当初始阶段伪距偏差的影响剧烈时，甚至会引起结果难以收敛的状态，此时伪距偏差改正模型的改进效果犹为明显[15-17]。图 5分别给出了GMSD测站和JFNG测站在DOY122/2015和DOY293/2014的单天坐标偏差序列。两个测站的结果有明显的差异，主要是由于初始阶段这两个测站的北斗卫星受到伪距偏差的影响程度不同，根据GMSD测站和JFNG测站当天北斗卫星的高度角序列，发现在前4 h(480个历元)，GMSD测站的IGSO和MEO卫星多处于高高度角状态，JFNG测站的IGSO和MEO卫星在初始阶段多处于中低高度角状态。伪距偏差值在低高度角和高高度角下最为明显，尤其是MEO卫星在高高度角情况下偏差值甚至达到1 m，因此造成了GMSD测站在PPP解算时难以收敛的结果，经过传统模型的伪距偏差改正后基本缓解了系统误差，而改进模型在传统模型的基础上大大改善了收敛效果和定位精度。JFNG测站受到伪距偏差的影响较弱，但仍可以看出改进模型相对与传统模型的定位精度已有所提高。
表 2 3种PPP处理策略的定位结果精度Table 2. Positioning Errors for BeiDou PPP with Three Different Strategies定位精度/m 无改正 传统模型 改进模型 E方向 0.189 0.172 0.156 N方向 0.157 0.143 0.127 U方向 0.429 0.414 0.295 ![]()
图 5 GMSD测站和JFNG测站3种PPP解算的结果对比Figure 5. Positioning Result for BeiDou PPP with Three Different Strategies in GMSD and JFNG Station4 结语
与其他卫星系统相比较，北斗卫星的伪距观测值中存在明显的系统偏差，需要通过模型改正的方法消除该系统偏差的影响，文章分析了测站数目和分布以及观测时间的长短对模型估计的影响，选择全球18个测站全年12个月的数据对北斗MEO卫星进行建模，4个可以观测到全弧段IGSO卫星的测站对IGSO卫星建模，在高度角范围5°~85°内，以5°为间隔进行分段，通过加权分段线性拟合联合抗差估计的方法，得到IGSO卫星和MEO卫星的伪距偏差的高度角模型。进一步测试了改正模型对PPP定位结果的影响，结果表明与不加模型改正的结果对比，定位精度在E、N、U方向平均提高了17.4%、19.1%、9.5%，相对传统模型E、N、U方向有8.4%、10.1%、8.5%的提高。
-
[1] Belikov M V. Spherical Harmonic Analysis and Synthesis with the Use of Column-Wise Recurrence Relations[J]. Manuscripta Geodaetica, 1991,16:384-410 [2] Belikov M V, K A Taybatorov. An Efficient Algorithm for Computing the Earth's Gravitational Potential and Its Derivatives at Satellite Altitudes[J]. Manuscripta Geodaetica, 1992,17:104-116 [3] Varshalovich D A, Moskalev A N, Khersonskii V K. Quantun Theory of Angular Momentum[M]. Singapore:World Scientific Publ., 1989 [4] Heiskanen W, Moritz H. Physical Geodesy[M]. San Francisco:W.H. Freeman and Co., 1967 [5] Paul M K. Recurrence Relations for Integrals of Associated Legendre Functions[J]. Bulletin Geodesique, 1978, 52(3):177-190 [6] Wang Jianqiang,Li Jiancheng,Zhao Guoqiang,et al. Geoid Undulation Computed Based on Clenshaw Summation[J]. Geomatics and Information Science of Wuhan University, 2010, 35(3):286-289(王建强,李建成,赵国强,等. 利用Clenshaw求和计算大地水准面差距[J]. 武汉大学学报·信息科学版, 2010, 35 (3):286-289) [7] Tian Jin, Bao Jingyang, Liu Yanchun. Clenshaw Summation in Constructing High Resolution Gravity Field from the Geopotential Coefficient Expression Model[J]. Geomatics and Information Science of Wuhan University, 2005, 30(10):905-908(田晋,暴景阳,刘雁春. 全球位系数模型构建高精度局部重力场的Clenshaw求和[J]. 武汉大学学报·信息科学版, 2005, 30 (10):905-908) -
期刊类型引用(8)
1. 张捍卫，李晓玲，张华，杨永勤. 缔合勒让德函数的列式递推公式. 大地测量与地球动力学. 2023(06): 622-626+640 . 
百度学术
2. 王鑫耀，李荣冰，刘刚，刘建业. 民用航空惯性基准系统重力补偿方法. 导航与控制. 2022(01): 49-56 . 百度学术
3. 邢志斌，李姗姗，田苗，范雕，张驰，马越原. 不同Legendre函数递推公式对计算球谐函数定积分的影响. 海洋技术学报. 2019(05): 37-42 . 百度学术
4. 李玉平. 超高阶重力场模型计算高程异常算法与实现. 海洋测绘. 2018(03): 1-4 . 百度学术
5. 魏子卿，杨正辉. Helmert扰动位及其积分核函数的椭球实用公式. 武汉大学学报(信息科学版). 2018(12): 1768-1774 . 百度学术
6. 万晓云，张润宁，李洋，刘波，眭晓虹. 基于球谐函数的重力异常和垂线偏差误差匹配关系. 测绘学报. 2017(06): 706-713 . 百度学术
7. 朱永超，万晓云，于锦海. 引力和垂线偏差的非奇异公式. 武汉大学学报(信息科学版). 2017(12): 1854-1860 . 百度学术
8. 铁俊波，吴美平，蔡劭琨，张开东，练军想. 基于EGM2008重力场球谐模型的水平重力扰动计算方法. 中国惯性技术学报. 2017(05): 624-629 . 百度学术
其他类型引用(6)
下载:
计量
- 文章访问数: 3901
- HTML全文浏览量: 266
- PDF下载量: 1765
- 被引次数: 14
