最小二乘估计

最小二乘估计
随着空间技术的发展,人类的活动开始进入了太空,对航天器(包括人造地球卫星、宇宙飞船、空间站和空间探测器等)的观测手段和轨道确定提出了很高的精度要求。

在计算技术高速发展的推动下,各种估计理论也因此引入到轨道估计方法中。

大约在1795年高斯在他那著名的星体运动轨道预报研究工作中提出了最小二乘法。

最小二乘法就成了估计理论的奠基石。

最小二乘估计不涉及观测数据的分布特性,它的原理不复杂,数学模型和计算方法也比较简单,编制程序不难,所以它颇受人们的重视,应用相当广泛。

对于严格的正态分布数据,最小二乘估值具有最优一致无偏且方差最小的特性。

实践证明,在没有粗差的情况下,大部分测量数据基本上符合正态分布。

这是最小二乘估计至今仍作为估计理论核心的基础。

最早的轨道确定就是利用最小二乘法,用全部观测数据确定某一历元时刻的轨道状态的“最佳”估值,即所谓的批处理算法定轨。

长期以来,在整个天体力学领域之中,各种天体的定轨问题,几乎都是采用这一方法。

卫星精密定轨的基本原理为：利用含有误差的观测资料和不精确的数学模型,通过建立观测量与卫星状态之间的数学关系，参数估计得到卫星状态及有关参数的最佳估值。

参数估计的基本问题就是对一个微分方程并不精确知道的动力学过程，用不精确的初始
状态*0X 和带有误差的观测资料，求解其在某种意义下得卫星运动状态的“最佳”估值0
ˆX 。

常用的参数估计方法有两种，最小二乘法和卡尔曼滤波方法。

最小二乘法是在得到所有的观测数据之后，利用这些数据来估计初始时刻状态量的值，由于用到的观测数据多、计算方法具有统计特性，因此该方法精度高。

卡尔曼滤波在观测数据更新后，利用新的观测数据对状态量进行改进得到这一观测时刻的状态量，卡尔曼滤波适用于实时处理。

卫星精密定轨输运高精度的事后数据处理，通常采用最小二乘法进行参数估计。

记观测量的权阵为 P 。

利用加权最小二乘法计算总的观测方程方程0y Hx ε=+，得
1()T T x H PH H Py -=
卫星的参考状态为**000
ˆX X x =+ 在精密定轨的过程中，由于状态方程和观测方程在线性化过程中会产生误差，上式的解算需要通过不断的迭代。

迭代的收敛准则为：
（1） 卫星位置的改正量小于某个事先指定的值，比如 0.01m ；
（2） 两次迭代的残差满足
1
j j j
RMS RMS RMS
μ--<，其中μ为小量，可取0.01。

由于观测资料较多，在参考状态处()*i X t 线性化得线性方程0i i i y H x ε=+经常严重“病态”，在实际计算中需要特别注意。

卫星轨道精密确定需要高精度的观测数据。

尽管目前的观测精度已大大提高,但在实际观测中,由于各种因素的影响,在观测数据中,往往存在着一定数量粗差。

所以,在轨道确定之前需要对观测数据进行预处理。

但在一般预处理中,并未考虑观测数据与轨道动力学模型的吻合情况,只能剔除大野值,仍有相当一部分质量不好的数据保留了下来,然后进行初步定轨,并由此进行残差分析。

由于观测数据存在着距离偏差和时间偏差,所以残差分析的过程为:首先对初步定轨所得的残差进行距离偏差和时间偏差改正;然后对改正后的残差进行多项式拟合;最后对拟合后剩余残差较大者予以剔除,以消除其中的粗差。

因观测误差中存在着偶然误差和系统误差,而且由于模型不准,上述方法还难以凑效,对预处理过的卫星激光测距数据进行 统计分析表明还有大量的粗差存在。

基于卫星激光测距的数据量十分庞大,经过预处理后还存在着大量的粗差,美国航空航天管理局(NASA)激光测卫工作组建立了一个激光测卫数据质量自动控制系统(AQC),负责对全球的激光测卫数据进行检验,发现观测数据的问题,查出问题原因并进行改正,对粗差予以剔除。

数据质量自动控制系统对Lageos 卫星的激光测距实测数据的检验表明,15%的激光测卫数据存在问题,仅5%能被该系统探测到,而10%有问题的数据要靠其它手段探测。

数据预处理和激光测卫数据质量自动控制系统都难以完全剔除粗差,美国德克萨斯大学空间研究中心卫星精密定轨与动力学测地软件“UTOPIA ”对各卫星观测台站软硬件进行综合评价,对精度较低的台站予以降权;对误差较大的观测数据或观测弧段予以剔除。

其它的卫星精密定轨软件数据预处理与“UTOPIA ”类似。

但对观测台站予以降权或者剔除一些观测数据存在人为性,理论不够完备,容易剔除一些观测精度很高的观测数据,付出的代价大,数据利用率低。

为此我们可以考虑在卫星精密定轨中引入抗差估计来减弱观测资料中粗差的影响。

应用抗差估计理论,开展卫星精密轨道确定的研究,可望更进一步提高定轨的精度。

在卫星的精密定轨中,法方程的解一般采用经典的最小二乘估计。

经典最小二乘估计在观测数据服从正态分布时,具有无偏、一致、优效性,但是当观测数据遭受异常污染时,最小二乘法又具有负面效应,即估值不具有抗干扰性。

抗差估计是在粗差不可避免的情况下,选择适当的估计方法使未知参数估值尽可能避免粗差的影响,得出正常模式下的最佳估值。

抗差估
计的原则是要充分利用有效信息,限制利用可疑信息,排除有害信息。

通过等价权,抗差估计原理与最小二乘形式有机结合起来形成了抗差最小二乘估计。

这样,即给经典最小二乘估计赋予了抗差能力,又保留了经典最小二乘估计算式简洁、计算方便的特点。

为抗差估计用于确定卫星定轨铺平了道路,因为己发展成熟的卫星定轨模型和算法,在形式上仍可保持,只是要加入抗差估计的算法即可。

下面简要介绍抗差最小二乘估计。

抗差最小二乘估计简介
设有一组相互独立的观测值
{}
i L ，1,2,,i m = ，观测值的权为
[]12,,,m P diag P P P = 相应的误差方程式可简单表示为：
V AX L =-
式中L 为m 维观测向量，A 为m ×n 系数矩阵,X 为n 维未知参数向量，
[]12,,,T
m V v v v = 为观测值残差向量，最小二乘估计的准则函数为：
2
min i i
Pv
=∑
可以构成法方程式：ˆ0T
T
A PAX
A PL -= 其解为：
()1ˆT T X
A PA A PL -= ()1
20ˆT x
A PA σ-=∑
()2
0ˆ/T V PV m n σ=-
其中m 为观测个数，n 为未知参数个数。

如果观测值{}i L 含有粗差,则参数估值及方差估值会被扭曲。

为了抵制异常值对ˆX
和2
0ˆσ
的影响,我们可建立如下准则函数： ()1
min m
i
i
i P v ρ==∑
或()()11
0m
m
i
i i i i ij i i i j v P v P v a v x ρ==⎡⎤∂∂=ψ=⎢
⎥∂∂⎣⎦∑∑，()1,2,,j n = 其中()ρ∙是一个凸的、连续函数， ij a 是系数矩阵A 的第i 行第j 列元素。

()()i i i
v v v ρ∂
ψ=
∂ 如果令()()1i i i i i i i i
i i v p p v p p w v v v ρψ⎛⎫
∂===
⎪∂⎝⎭
称i p 为等价权，i w 为权因子。

于是可以得到抗差最小二乘的法方程式：
ˆ0T T A PAX
A PL -= 其他解为：
()1ˆT T X
A PA A PL -= ()1
20ˆT x A PA σ
-=∑ ()2
0ˆ/T V PV m n σ=-
其中P 为等价权矩阵,m 为不计0ii P =的测值个数。

可以看出,抗差最小二乘估计的抗差实质体现在等价权上,其关键是建立恰当的权函数。

为了得到既有较强抗差性,又有较高效率的估值,权函数应包含两方面的内容: ①观测值的信息区间划分为正常观测值(有效信息):可利用观测值(可利用信息):粗差观测值(有害信息)。

②根据这三部分观测值,可将权划分为:保权区(保持原观测值不变);降权区(对观测值作抗差限制):拒绝区(权为零)。

抗差最小二乘估计的基本效率由保权区的观测值来保证,它们应是观测数据的主体。

抗差最小二乘估计的效率和可靠性通过降权区的观测值的权函数得到加强。

它的抗差能力更体现在拒绝区。

故抗差最小二乘估计基本上可称为“三段法”。

权函数的选择应由观测数据的结构决定,下面简要介绍适合卫星激光测距的IGGIII 观测抗差方案。

IGGIII 选取如下等价权函数：
100110
1
ij
j j
ij ij
j j j p u k k u k p p k u k k k u u k ⎧≤⎪⎪⎛⎫
-⎪ ⎪=<≤⎨ ⎪-⎪⎝⎭
⎪>⎪⎩
其中0/j j u v σ=，0σ为单位权中误差，0k 的取值通常在1.0到1.5之间，1k 的取值通常在2.5到3.0之间。

令()
()110/j j d k u k k =--，称为平滑因子，01j d ≤≤，于是
20011
ij j ij ij j
j j j p u k k p p d k u k u u k ⎧≤⎪⎪⎪
=<≤⎨⎪⎪>⎪⎩
IGGIII 方案对应的权函数采用三段法,即正常段采用最小二乘估计，可利用段采用平滑因子降低观测值的权，拒绝段通过给观测值的权赋零值来淘汰粗差。

抗差批处理
设有一组观测值L ，它的权阵P ，其中，L 是n ×l 维列向量，P 是的方阵(各观测值相互独立时为对角阵)，则误差方程式可简单表示为
V AX L =-
式中0c L L L =-，*
X X X =-，下标O 、C 分别代表观测值和计算值，X 、*
X 为状态矢量(包括卫星的位置、速度，以及待估的动力学和运动学参数)的真值和参考解，A 为观测值对状态矢量的偏导数矩阵。

采用抗差最小二乘法可得：
()1ˆT T X
A PA A PL -= ()1
20ˆT X A PA σ-=∑
()2
0ˆ/T V PV m n σ=-
其中P 为等价权矩阵，m 为不计0ii p =的测值个数,n 为待估参数的个数。

批处理的抗差估计与最小二乘估计收敛准则相同，只是观测残差均方差(RMS)和观测残差均方差的线性预报值(RMSP)稍有不同,它们分别为:
()
12
1T k k
i m
RMS v p v p ⎡
⎤
⎢⎥⎡⎤=⎣⎦⎢⎥⎢⎥⎣⎦
∑ ()()1
2
111
ˆˆT k k k k k k m RMSP v H x p v H x p ++⎡⎤
⎢⎥=--⎢⎥
⎢⎥⎣⎦
∑ 式中k
v 第k 次迭代中观测数据的残差，
m
p ∑表示对参加估值计算的观测数据的等价权之和。

以上可以看出采用等价权很容易将批处理转化为抗差批处理算法。

但因等价权是残差的函数,计算开始时一般得不到精确的残差,故需迭代使等价权的准确度得到改善。

在采用抗差批处理进行卫星精密轨道确定时,需采用尽量准确的力学模型及测量模型,同时初始轨道根数也需尽可能的准确。

SLR 的观测资料包括测站与卫星之间的距离、卫星的方位角和高度角。

虽然SLR 的观测站是全球分布的，但是观测站的数目比较少，进行 SLR 观测的卫星数目比较多，北斗卫星的观测资料比较少，没有多站同步观测数据。

单站混合资料定初轨的基本原理是：在测站坐标已知的前提下，利用 SLR 的距离观测值和角度观测值可以计算出卫星在观测时刻的位置矢量；然后利用多项式拟合法求出卫星初
始时刻0t 的位置0r 和速度矢量的初始值0r 、0
r
；假设卫星轨道为椭圆，则任意时刻t 的卫星位置r 和速度r 可以由0t 时刻卫星位置0r 和速度0
r
表示，从而得到t 时刻卫星位置与初始时刻卫星状态之间的函数关系，最后由参数估计方法计算卫星0t 时刻的状态。

记观测站在地固坐标系下的坐标为(,,)T i i i i R X Y Z
，观测时刻i t 时测站与卫星之间的距
离i ρ、站心坐标系中卫星的方位角i A 、高度角i h ，则 SLR 单站混合资料定初轨的具体算法如下：
（1） 将 SLR 原始完整观测数据进行预处理，具体的预处理方法见第一节； （2） 计算i t 时刻导航卫星在 J2000.0 惯性系中的坐标(,,)T i i i i r x y z
： 计算卫星在站心坐标系中坐标，具体公式为：
cosh cos cosh sin sinh si i i si si i i i si i x A r y A z ρ⎡⎤⎡⎤⎢⎥⎢⎥==⎢⎥
⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦
将站心坐标系下卫星的位置矢量si r 转换成地固坐标系下卫星的位置矢量ei r
：
()()18090ei z y si i r R R r R λϕ=--+
最后将地固坐标系下的位置矢量ei r 转换为 J2000.0 坐标系下的位置矢量i r
：
()i z ei r R GAST r =-
其中，λ为测站经度，ϕ为测站纬度；z R 、y R 由公式()1000cos sin 0sin cos x R θθθθθ⎡⎤
⎢⎥=⎢⎥⎢⎥-⎣⎦
、()cos 0sin 0
10sin 0cos y R θ
θθθ
θ-⎡⎤
⎢⎥=⎢⎥⎢⎥⎣⎦
计算；GAST 表示格林尼治真恒星时，具体参考GMST 与儒略日 JD 之间转换公式。

（3） 估算0t 时刻卫星的位置0r 和速度0
r
： 假设i t 时测卫星的位置矢量(,,)T i i i i r x y z
是时间i t 的 m 次多项式：
()()0
1,2,,j
m
i j i j r d i n τ===∑
其中，
0i i t t τ=-，n 为观测值个数，m n ≤。

由最小二乘法求解
()()01,2,,j
m i j i j r d i n τ===∑ ，得01ˆˆˆ,,,m
d d d ，则r 的拟合值为()0
ˆˆm
j i j i j r d τ==∑，0r 的初值为0ˆd ，0r
的初值为1
ˆd 。

（4） 精化0t 时刻卫星的位置0r 和速度0
r
： 初轨计算时，只考虑二体问题对卫星的影响，卫星的运动轨道为椭圆，由开普勒定律可
以推导出i t 时刻卫星的位置i r 可以由0t 时刻卫星的位置0r 和速度0
r
来表示： 00
i i i r f r g r =+ 其中：
()()03211cos sin i i i
i i i a f E r g a E E τ⎧
=--∆⎪⎨
⎪=-∆-∆⎩
上式中0r =
，0v =0
2
00
2r a r v =
-，i E ∆而需要由下式迭代 计算：
)3201sin 1cos i i i i r E a E E a τ-⎛⎫
∆=+-∆--∆ ⎪⎝⎭
式中的00000000r r
x x y y z z =++ ，迭代的初始值可设为3i i E a τ-∆=，收敛条件为连续 两次迭代结果之差的绝对值小于某个设定的极小值，即,,1i k i k e E E ε-∆-∆<。

在计算得到i f 和i g 后，通过最小二乘法解公式00
i i i r f r g r =+ ，求解0r 和0r
，计算公式为：[][][][][][][][][][][][][][]02
2gg fr fg gr r ff gg fg ff gr fg fr r ff gg fg ⎧-=⎪-⎪
⎨-⎪=⎪-⎩
上式中[]1
n i i
i ff f f ==
∑，[]1
n i i
i gg g g ==∑，[]1
n
i i
i fg f g ==∑。

将得到的0r 和0r 代入公式[][][][][][][][][][][][][][]02
2gg fr fg gr r ff gg fg ff gr fg fr r ff gg fg ⎧-=⎪-⎪⎨-⎪=⎪-⎩
，重新计算此时的i f 和i g ，然后
将i g 和i f 代入公式[][][][][][][][][][][][][][]02
2gg fr fg gr r ff gg fg ff gr fg fr r ff gg fg ⎧-=⎪-⎪
⎨-⎪=⎪-⎩
，重复上述过程进行迭代，直到连续两次的
迭代结果满足()()0011k k r r ε--<，()()0021k k r
r ε--< 且()()0031k k a a ε--<时，停止迭代，得到满足精度要求的0r 和0
r。

对于导航卫星，可以利用卫星的伪距/载波相位观测数据进行精密轨道确定。

与 SLR 观测数据不同的是，伪距/载波相位观测值只有距离观测量，可以进行多站同步观测。

多站同步观测资料定初轨的基本原理是：在测站坐标已知的前提下，利用多站伪距/载波相位观测值计算卫星在观测时刻的位置矢量；之后的算法与单站混合资料定初轨的算法相同。

具体算法如下：
已知 n 个观测时刻 m 站同步距离观测资料i t 、
ij ρ，地固系下测站坐标为
(,,)j j j j R X Y Z ，其中i=1，2，…，n ， j=1,2，…，m 。

则利用 m 站伪距/载波资料定初轨的具体算法如下：
（1） 将伪距/载波观测数据进行预处理，重点剔除钟差、电离层等误差项的影响； （2） 计算i t 时刻卫星 J2000.0 惯性系中的坐标(,,)T i i i i r x y z
： 将i t 时刻测站的坐标从地固坐标系转换到 J2000.0 惯性系下：
()ij z j R R GAST R =-⋅
其中， GAST 表示格林尼治真恒星时，由GAST 与儒略日 JD 之间转换公式计算。

根据几何关系，有:
ij ρ将上式两边平方得
()2222ij i ij i ij i ij i ij r R x X y Y z Z ρ=+-++
其中，2222i i i i r x y z =++，2222
ij ij ij ij R X Y Z =++。

()222
2ij i ij i ij i ij i ij r R x X y Y z Z ρ=+-++式中有三个未知参数i x 、i y 、i z ，理论上说，
只要三个方程便可求解上述方程。

但由于该方程组非线性，求解方程组过程很复杂。

考虑用
迭代的算法求解i x 、i y 、i z ，将()
222
2ij i ij i ij i ij i ij r R x X y Y z Z ρ=+-++变换成如下形式：
()222
12i ij i ij i ij ij i ij x X y Y z Z r R ρ++=
-++ 假设i r 的估计值为0i r ，带入公式()222
12
i ij i ij i ij ij i ij x X y Y z Z r R ρ++=-++得到m 个线性方
程，用最小二乘法解出i x 、i y 、i z ，然后再利用公式2222i i i i r x y z =++计算1i r ；将1i r 带入公式()222
12
i ij i ij i ij ij i ij x X y Y z Z r R ρ++=
-++重新求解此时的i x 、i y 、i z ，再将其带入公式2222i i i i r x y z =++计算2i r ，这样就完成两个迭代过程。

得到连续三次的0i r 、1i r 、2i r 的值，
可以利用加速手链公式
()
10*00210
2i i i i i i i r r r r r r r -=-
-+
用上式可以得到较好的初值*0i r ，用*
0i r 替代0i r
，重新进行上述两个计算过程，然后用()
10*0
0210
2i i i i i i i r r r r r r r -=-
-+迭代，得到满足精度的i r 和i x 、i y 、i z 。

（3）估算0t 时刻卫星的位置0r 和速度0
r
：
假设i t 时测卫星的位置矢量(,,)T i i i i r x y z
是时间i t 的 m 次多项式：
()()0
1,2,,j
m
i j i j r d i n τ===∑
其中，
0i i t t τ=-，n 为观测值个数，m n ≤。

由最小二乘法求解
()()01,2,,j
m i j i j r d i n τ===∑ ，得01ˆˆˆ,,,m d d d ，则r 的拟合值为()0
ˆˆm
j i j i j r d τ==∑，0r 的初
值为0ˆd ，0r
的初值为1
ˆd 。

（4）精化0t 时刻卫星的位置0r 和速度0
r
： 初轨计算时，只考虑二体问题对卫星的影响，卫星的运动轨道为椭圆，由开普勒定律可
以推导出i t 时刻卫星的位置i r 可以由0t 时刻卫星的位置0r 和速度0
r
来表示： 00
i i i r f r g r =+ 其中：
()()03211cos sin i i i
i i i a f E r g a E E τ⎧
=--∆⎪⎨
⎪=-∆-∆⎩
上式中0r =
，0v =0
2
00
2r a r v =
-，i E ∆而需要由下式迭代 计算：
)3201sin 1cos i i i i r E a E E a τ-⎛⎫
∆=+-∆--∆ ⎪⎝⎭
式中的00000000r r
x x y y z z =++ ，迭代的初始值可设为3i i E a τ-∆=，收敛条件为连续 两次迭代结果之差的绝对值小于某个设定的极小值，即,,1i k i k e E E ε-∆-∆<。

在计算得到i f 和i g 后，通过最小二乘法解公式00
i i i r f r g r =+ ，求解0r 和0r
，计算公式为：[][][][][][][][][][][][][][]02
2gg fr fg gr r ff gg fg ff gr fg fr r ff gg fg ⎧-=⎪-⎪
⎨-⎪=⎪-⎩
上式中[]1n i i i ff f f ==∑，[]1n i i i gg g g ==∑，[]1n
i i
i fg f g ==∑。

将得到的0r 和0r 代入公式[][][][][][][][][][][]
[][][]0202gg fr fg gr r ff gg fg ff gr fg fr r ff gg fg ⎧-=⎪-⎪⎨-⎪=⎪-⎩
，重新计算此时的i f 和i g ，然后
将i g 和i f 代入公式[][][][][][][][][][][]
[][][]0202gg fr fg gr r ff gg fg ff gr fg fr r ff gg fg ⎧-=⎪-⎪⎨-⎪=⎪-⎩
，重复上述过程进行迭代，直到连续两次的迭代结果满足()()0011k k r r ε--<，()()0021k k r
r ε--< 且()()0031k k a a ε--<时，停止迭代，得到满足精度要求的0r 和0
r 。