随着空间测量技术的发展，对卫星轨道精度的要求越来越高，现有的精密定轨软件程序在解算参数时，都是假设计算误差服从正态分布，采用最小二乘估计进行解算。文献^[1]指出，在数据处理过程中，由于观测数据中存在粗差、模型不准确以及计算方法不同等，使得观测误差并不服从正态分布。文献^[2]从观测误差的实际分布出发，利用LP估计来确定卫星的精密轨道，通过计算，认为LP估计的p值应该在[1, 2]之间；文献^[3]指出，对GPS观测值误差处理时，当p=1.4左右时精度最高，效率最好。

在实际工作中，对于一组给定的观测误差，往往并不知道这个p取什么值最合适，这就需要利用观测数据提供的信息，正确确定p值，进行相应的估计。此前，许多学者采用极大似然估计法对P范分布的多种形式，如从对称到非对称、一元到多元等情况进行了研究，得到了一些有用的结果^{[4, 5, 6, 7, 8, 9]}。当观测样本足够大且分布的形式已知时，极大似然估计是最优估计，其渐近分布为正态分布。文献^[10]将半参数回归模型成功地应用到极大似然平差法的参数估计理论中，但形状参数p值的获取是半参数P范极大似然回归的应用前提。然而，P范分布的形状参数值未知，极大似然估计方程是非线性方程，需要迭代计算，迭代计算需要较大的数据量及合适的初值才能保证估计精度，此时自适应LP估计等价于相应的最或然估计^[11]，而形状参数的不正确估计又会使得LP估计的效率变低^[12]。理论证明，当p=1时，LS估计将比L₁估计损失效率50%；当p=2时，L₁估计又将比LS估计损失效率27%，因此，选择合适的p非常重要^[12]。文献^{[2, 13, 14, 15]}以SLR数据处理为例说明了自适应LP估计在应用中的稳健性优势，若能获得更可靠的p值，其优势将更为显著。

合理的p值可以体现观测误差最符合的分布，由于p为小数，常用的矩估计法为整数阶估计，本文引入实数阶和对数矩估计方法，导出了相应估计方法的P范分布的参数解算的基础方程，另外，简化了基础方程，消除了普西函数和伽玛函数的导数的计算误差对参数计算的影响，使得参数计算的复杂度降低。

设L_i(i=1,2,…，n)是均值μ的一组独立的观测值，服从P范分布_i^[1]，L_i的先验方差为σ²=σ₀²q_i，其中σ₀²为单位权方差，σ_i²为L_i的先验方差，令ω_i=σ₀/σ_i=1/$\sqrt {{q_i}} $，则L_i的密度函数为

点击浏览公式

点击浏览公式

令t=λ=，则dX=λ=，将dX和t代入式(2)，计算得

点击浏览公式

式(3)包含形状参数p、中误差σ₀和阶数r，为了减少参数之间的影响，可以构造如式(4)所示的函数，消除中误差σ₀对p的影响

点击浏览公式

点击浏览公式

点击浏览公式

将式(5)、式(6)代入式(4)，化简得到形状参数的估计式为

点击浏览公式

由式(7)得到形状参数p的估计值后，在式(3)中令r=p，得

点击浏览公式

确定一个r值后，由式(5)、式(7)、式(8)就可以逐步得到参数的估计值${\hat u}$、 ${\hat p}$、${\hat \sigma }$ ₀。从后面的模拟算例可以看出，当阶数r的取值越接近参数p，由实数阶矩估计法计算得到的参数估计量相对偏差越小，其估计精度越高。而r的选取方式比较复杂，一般来说，可以选择多个r值，分别进行平差计算，然后对计算结果进行比较，从中选取一个较好的结果。

从上面的推导过程可知，利用实数阶矩估计法得到的参数估值精度与阶数的选取有较大影响，在实际计算过程中，无法选取合适的阶数r，因此，不能确定参数p的精确取值。在计算过程中，伽玛函数、普西函数的截断误差对参数估计值的计算影响也比较大，为了消除阶数及计算误差对参数估计值精度的影响，采用对数矩估计方法来估计相应参数的估计值${\hat u}$、 ${\hat p}$、${\hat \sigma }$₀。

对式(3)求导数得

点击浏览公式

对式(9)化简得

点击浏览公式

在式(9)中，令r=0 可得观测值的一阶对数绝对矩为

点击浏览公式

用样本矩代替总体矩，综合式(10)、式(11)，化简得参数p的估计表达式为

点击浏览公式

为了消除ψ(·)计算误差及阶数对参数计算的影响，在式(12)中，取r=p，由ψ(·)的性质得

点击浏览公式

由式(8)得到参数u的估计表达式为

点击浏览公式

求解式(8)、式(13)、式(14)，就可得到未知参数的相应估计值${\hat u}$、 ${\hat p}$、${\hat \sigma }$。

由于式(8)、式(13)、式(14)是非线性方程组，可以用迭代法解算。首先联立求解式(13)、式(14)，得${\hat u}$、 ${\hat p}$，然后将得到的估计值代入式(8)，直接计算即得到单位权方差的估计值。采用迭代法求解式(13)、式(14)时，一般取ω_i=1，实数阶矩估计方法计算得到的参数估计值作为μ、p的初始值，令

点击浏览公式

点击浏览公式

将f₁(u+Δu,p+Δp)=0，f₂(u+Δu,p+Δp)=0在μ、p处展开成泰勒级数，略去二次以上各项，则有

点击浏览公式

点击浏览公式

联立求解式(15)—式(18)，将μ+Δμ、p+Δp作为新的近似值，即可进行下一次迭代计算。其中，式(15)、式(16)分别对μ、p的偏导数为

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

令B₁=则式(17)、式(18)变为

点击浏览公式

求解方程组(23)，如果Δμ、Δp小于给定的迭代阀值ε_μ及ε_p，则停止迭代，否则用μ+Δμ、p+Δp作为未知参数的近似值，重新计算式(23)系数阵式及常数项，进而求解该方程。最后可由式(8)求出观测值的单位权方差估计值。

为了检验理论的正确性，便于进行对比分析，利用文献^[1]的两组来自不同母体的观测数据验证本文方法的估计效果。为了分析阶数对精度的影响，不妨取阶数为0.05、0.1、0.5、1、1.5、2、5，利用估计量的相对偏差(${\hat p}$|－p|/p)表示参数估计的精度，拉普拉斯分布和正态分布的参数估计量的相对偏差曲线如图 1、图 2所示。

图 1 p=1时，不同r对应估计量p的相对偏差曲线 Fig. 1 The relative deviation curve of estimation p with different r when p=1

图选项

图 2 p=2时，不同r对应估计量p的相对偏差曲线 Fig. 2 The relative deviation curve of estimation p with different r when p=2

图选项

从图中可以看出，阶数r的取值越接近形状参数实际值，参数估计值越接近其真实值，其相对偏差也越小。因此，可以认为当阶数r近似取p值时，实数阶矩估计法的估计精度最高。不妨以拉普拉斯分布为例，其样本个数、阶数及形状参数的关系如表 1所示。

由表 1可以看出，当子样较少时，无论阶数r取何值，计算得到的参数估值均偏离真实值，阶数与形状参数的差值越大，精度越低。随着参与平差的子样数逐步增多，对应各阶数计算得到的参数估值逐渐趋近于真实值。

为了分析对数矩和极大似然估计方法^[1]的稳定性、收敛速度和精度，利用上述模拟观测数据进行对比分析。由式(23)和式(8)，得到不同样本情况下的参数估值${\hat u}$、${\hat p}$和${\hat \sigma }$，并与文献^[1]的估计结果进行对比，如表 2、表 3所示。

从表 2、表 3可以看出，当子样较少时，对数矩估计法与极大似然平差法估计结果均与实际值相差较大，随着观测子样的增加，两种方法所得到的参数估计值呈现出越来越接近真实值的趋势。当子样数较小，在估计3个参数时，对数矩估计法在收敛速度、稳定性和精度上明显优于极大似然估计法。

根据表 2、表 3中的计算结果，得出运用对数矩估计法和极大似然估计法获得${\hat p}$的${\hat \sigma }$和参数估计量的相对偏差曲线，如图 3、图 4所示。

图 3 p=2时，估计量p的相对偏差曲线 Fig. 3 The relative deviation curve of estimation p when p=2

图选项

图 4 p=2时，估计量σ的相对偏差曲线 Fig. 4 The relative deviation curve of estimation when p=2

图选项

4.3 实际算例1

由文献^[16]可知，地图数字化后的数据，经过标准化处理之后得到的误差，经检验得知，该母体服从p=1.6的P范分布。本文数据来自文献^[17]，利用对数矩估计法计算地图数字化后坐标误差分布的参数估值${\hat u}$、${\hat p}$和${\hat \sigma }$，同时与极大似然估计法得到的估计结果进行对比，计算结果如表 4所示。

通过表 4可以得出结论，当样本服从p=1.6的P范分布时，采用对数矩估计法估计形状参数p值的估计精度要高于极大似然估计法，而对均值u和方差σ的估计效果相当，同时对数矩估计法的迭代次数要远小于极大似然估计法，对数矩估计法的收敛速度较快，表明对数矩估计法能快速定位形状参数p。

由文献^[3]可知，经消除GPS观测值中的粗差，完成电离层折射误差、对流层折射误差等多项改正后，GPS输出的最终观测值误差，经检验得知，该母体服从p=1.4的P范分布。

本文算例所用数据来自加拿大Algonquin Park的ALGO测站点，利用TPS NET-G3A接收机，采集获得2013年4月28日全天的观测数据。在获得的32颗卫星数据中，选取某颗卫星伪距的PPP双频无电离层组合观测值残差进行分析。取200个误差值作为原始数据，分别利用对数矩估计法与极大似然估计法对参数u、p和σ进行计算，得到的结果如表 5所示。

从表 5可以看出，针对本组试验数据，当形状参数p为1.4时，两种估计方法的效果均较好，其中对数矩估计法的精度稍高于极大似然估计法，且其迭代次数明显少于极大似然估计，因此采用对数矩估计法能够更快速精确地估计出误差的实际分布参数，再次表明对数矩估计法优于极大似然估计法。

本文研究了P范分布的实数阶矩估计方法，利用模拟算例讨论了阶数r与参数估计精度的关系，推导出形状参数、期望和中误差的估计公式，引入对数矩估计方法，消除了伽玛函数、普西函数的截断误差对参数估计值计算造成的影响，给出了相应参数的求解方法和计算过程，用模拟算例和实测数据验证了本文理论的正确性，得到了以下结论：

(1) P范分布能够精确地确定观测数据的最符合分布，本文给出了任意实数阶矩估计公式。利用实数阶矩估计法得到的参数估计值精度与阶数的关系较大。当阶数r的取值越接近p，由实数阶矩估计法计算得到的参数估计量相对偏差越小，其估计精度越高。然而在参数估计过程中，r的选取方式比较复杂，合理的阶数的选取方法值得进一步研究。

(2) 对数阶矩估计方法降低了算法的复杂度，计算简单。从模拟和实测数据可以看出，对数估计方法在参数估计的精度、稳定性和收敛速度以及算法的复杂度等方面优于极大似然估计方法。对于模拟数据，实数阶和对数矩估计法略优于极大似然估计法，对于实测数据，对数矩估计法优势较明显，是一个较好的参数估计方案。

(3) 观测数据的分布受外界因素和样本数量的影响，采用P范分布能够较好地估计参数，接近误差的真实分布。观测样本数越大，计算得到的分布越接近它的真实分布，估计精度也越高。