导弹与航天运载技术
MISSILES AND SPACE VEHICLES
1999年 第5期 No.5 1999



基于延拓法的卫星初轨算法
——应用于单个观测物体确定低轨卫星初始轨道
洪炳熔　刘艳芳　胡芳
　　摘要　研究了用单个观测物体进行低轨卫星初轨计算的问题，给出了一种大范围收敛的解算方法：Newton割线法和中点求积法（延拓法的一种)结合起来的综合法。计算机仿真结果说明：整个算法能根据较短弧段的测量数据计算用户星的初轨轨道，并保证一定的精度。
　　主题词　轨道计算，全球定位系统，非线性方程。
Preliminary Orbit Determination of Satellite Based on Secant-Continuation Algorithm
——Applied in Preliminary Orbit Determination Using Single Observer
Hong Bingrong　Liu Yanfang
（Harbin lnstitute of Technology ，Harbin,150001）
　　Hu Fang
(Beijing University of Aeronautics and Astronautics,Beijing,100083)
　　Abstract　Preliminary orbit determination of low orbit satellite using single observer is studied.An algorithm combining secant and continuation is presented, which has property of global convergence.The result of computer simulation shows that the algorithm can calculate the preliminary orbit of user satellite correctly with finite measurement data in short orbit arc.
　　 Key Words　Orbit calualation ,Global positioning system，Nonlinear equation.
1　前　言
　　轨道计算的任务在于根据观测数据和卫星的自身运动规律给出真实卫星轨道的近似［1］。早期的轨道计算多用角度作为测量量，而今天的激光测距数据的精度已达到20 mm，多普勒测速精度可达到3 mm／s［2］，所以将距离和距离变化率作为观测量的卫星定轨系统越来越多。观测物体可以是卫星地面站，也可以是高轨道的卫星，如：全球定位系统（GPS）卫星、数据中继卫星等。
　　1996年11月，北京卫星信息工程研究所信息室采用的高动态GPS接收机在某卫星上搭载成功，这是第1次利用星载GPS接收机对低轨飞行器进行定轨。目前，这也是中国唯一的星载GPS接收机。利用GPS为低轨卫星定轨，对于某一确定时刻t，通常需要至少4颗GPS卫星的信息才能计算出用户星的位置和速度［3］，但在定轨初期，星载GPS接收机是一颗一颗搜索GPS卫星的，即先跟踪到一颗，然后是两颗GPS卫星。这段时间内，接收机导航处理器无工作可做。如果在跟踪到一颗导航星时，能利用空闲的导航处理器计算出用户星的粗略状态，则可以缩小多普勒频移的搜索范围，利用这些初始信息，就可以缩短搜索其它GPS卫星的时间，尽快进入精密定轨的状态。
　　利用高轨道的中继星或地面站为低轨道的用户星定轨（观测量为距离和距离变化率）的系统中，由于中继星和地面站数量较少，在某一时刻，用户星一般只能和1个～2个观测物体建立观测关系，所以，轨道计算大多采用动力学方程和测量方程相结合的卡尔曼滤波方式。卡尔曼滤波需要卫星的初始状态启动递推算法，好的初值可以提高滤波器的收敛性，加快收敛速度，从而使系统尽快进入精密定轨状态［1，4］。
　　上面所述的就是卫星的初轨计算问题：在卫星轨道确定初期，用尽量短的时间和尽量少的观测数据（较短的卫星弧段）计算出卫星的粗略轨道状态［5］。文献［1］对这个问题作了较深入的研究，用改进的Newton-Raphson方法求解非线性的测量方程，但它所用的数据都是基于二体问题下的理论卫星轨道。对于包含各种摄动项的实际卫星轨道，用这种方法难以得到收敛的解。
　　这类问题以GPS作为导航星的情况最为复杂（除了求三维位置和速度外，还要求钟差和频差），其它情况可以作为它的特例，所以，本文以跟踪到一颗GPS卫星为例，使用基于实际卫星轨道的测量数据。为了使算法收敛，这里采用了求解非线性方程组的割线法和延拓法结合的综合算法。作者对整个算法进行了计算机仿真，实验结果证明，割线延拓法（具体采用了中点求积法）具有大范围的收敛性，对于解决上述问题比较有效。
2　测量方程的建立
　　当用GPS为用户星进行定轨时，其任务是通过观测量，即伪距离和伪距率计算出卫星的三维位置、速度、用户星钟差和频差共8个未知参数。在某一指定时刻点，只要知道4颗GPS卫星的轨道状态信息和观测量，就可以建立8个方程，仅用几何方法计算即可。但在初轨计算中，通常跟踪到的GPS卫星少于4颗，在一个时刻点上的观测量是不足以解算8个未知参数的［3］。以用户星只跟踪到一颗GPS导航星的情况为例，从理论上讲，需要4个时刻点的观测量，才能建立8个方程，所建立的具体的测量方程为


(1)
在所建立的测量方程组中，为所要求的未知向量，即卫星状态，向量中的元素分别为参考时刻（第0时刻）用户星的三维位置、钟差、三维速度和频差；为第τ时刻用户星的状态；X2τ，X3τ与之类似。为导航星的三维位置和速度；为第τ时刻GPS导航星的状态。为测量到的导航星到用户星之间的伪距和伪距率；为第τ时刻的伪距和伪距率。对式(1)进行解算，就可求出用户星的状态。虽然8个方程中有32个未知量，但X,Xτ，X2τ,X3τ之间不是独立的，根据它们的内在关系，可以将所有的未知量归结到参考时刻。
3　卫星状态预测
　　在第2节讨论的非线性方程组中，8个方程包含32个未知量，显然无法直接求解。但是卫星运动有其自身的规律，不同时刻的卫星状态是由卫星的动力学方程来描述的。
　　一般情况，卫星的运动方程是以惯性系为参考坐标系给出的，GPS卫星采用了WGS-84地固坐标系，所以在计算过程中必须时刻考虑坐标系的转换，以及地固系作为旋转坐标系所产生的柯氏速度和柯氏加速度。
　　由于卫星的动力学方程是微分方程的形式，所以不可能直接用X来计算Xτ。我们选择了F，G的封闭表达式［5］来进行轨道预测。F，G的封闭表达式是由轨道根数表示的，轨道根数表示的是惯性系中的卫星轨道。因此，必须将卫星的位置和速度从WGS-84坐标系转换到惯性系中，在惯性系中求出预测的卫星状态。显然，这使得编程较为复杂。但实验的结果证明：F，G方法精度高，计算量小。
　　F，G展开式为

(2)
其中　为平均角速度；E为偏近点角；E0为第0时刻的偏近点角。
　　对F，G求导可得：

(3)
其中　

(4)
　　对于钟差和频差，下面简单的递推公式显然是合理的：

(5)
4　求解非线性方程组的割线延拓算法
　　合适的非线性方程组的解法将是解决本问题的关键。
　　在第2节中建立的测量方程组是非线性的。非线性方程组通常采用数值解法，这是一个迭代过程。Newton法是求解非线性方程组的一个最基本而且十分重要的方法，目前使用的很多有效的迭代法都是以Newton法为基础，并由它发展而得到的。但Newton法具有局部收敛性，只有初始值与真解充分靠近时，才能使迭代序列收敛。
　　从第3节可知，从X推导Xτ,X2τ,X3τ的过程中，为了使算法实用，采用了理想的二体物体轨道模型，它和包含各种摄动项的实际轨道之间有一定的误差，这种轨道递推误差随递推时间的增加而以指数级增大。同时，测量采样时间间隔τ也不能无限制地缩小，因为过小的τ将使方程之间产生相关性，从而使方程组无法解算。仿真实验说明，在求解非线性方程组的过程中，3τ的时间间隔所产生的轨道递推误差非常大，严重影响迭代过程的收敛方向，使基本的Newton法，甚至使具有大范围收敛性的Newton下降法发散。
　　为了解决上面的问题，引入了延拓法［6］。延拓法的优势在于，它对初始值X0没有严格限制，是一种超大范围收敛的算法。对于非线性方程组

(6)
X*是方程组的解。延拓法的思想就是构造一族映象H∶D×［0,1］Rn＋1Rn,代替单个映象F，使H满足条件：

（7）
其中　F0(X)=0的解X0为已知，H（X，1)=0就是非线性方程组（6）的解，这样问题就转化为求同伦方程组。

（8）
的解X=X(τ)的解。
　　中点求积Newton法是根据上面的延拓思想构造的一种具体实用的算法，它的迭代公式是

(9)
　　对于测量方程组中含Xτ,X2τ,X3τ的方程，无法计算出其导数值F′(X)，也就无法给出整个方程组的Jacobi矩阵，所以具体的计算需要使用离散牛顿法，即割线法，其公式为［6］

(10)
其中　ej为第j个分量为1的单位坐标向量。另外，为保持Newton法的平方敛速，可取hkj=这里，为不等于0的常数。
　　中点求积法中，步长N的选取很关键，它不但关系迭代过程是否收敛，也关系到整个算法的效率。虽然，文献［7］在理论上讨论了最优步长的选取，但计算量太大，很不实用。为了减少计算量，采用了循环中点求积法处理：先固定N，根据式（9）由X0计算到XN,然后以XN作为新初值X0重新用式（9）作参数微分法；当方程组的向量2范数（Euclid范数）小于某个门限值，或第i次参数微分法所得的Euclid范数大于i＋1次的Euclid范数时，停止循环；最后用第i次循环所得的XN作为初值启动Newton法，以求得方程组的解。
5　循环中点求积算法
第1步　赋初值：给定初值X0（用户星的位置、速度、钟差和频差八维向量），中点求积步数N，最大循环次数M，Newton法迭代的最大循环次数K，精度门限值esp1和esp2，循环指针i=0。
第2步　计算F(X0)，并用割线法求Jacobi矩阵F′（X0），解方程组F(X0)b=F′(X0),得b0=［F′（X0）］-1F（X0）。
第3步　由式（9）计算
第4步　逐步循环求得Xk-1,Xk,由式(9）计算其中令
共循环N次就可以完成一个中点求积过程。
第5步　如果或i=M-1，执行第8步，否则作第6步。
第6步　当i≥1时，如果第i+1次的‖F(XN)‖小于i次的‖F（XN)‖，执行下一步，否则执行第8步。
第7步　，i＋＋，循环执行第3步到第5步。
第8步　，根据式（10)，采用Newton割线法计算Xk+1=Xk-bk，其中令bk=［F′(Xk)］-1F(Xk),循环执行，同时令k＋＋，当时，执行下一步；当k=K-1时，重新设置初值，从第1步开始执行。
第9步　打印Xk,F（Xk)，和理论轨道的误差等，结束。
　　需要说明的是，第1步所说的初值可根据猜测法和线性变换法求得［8］。
6　计算机仿真结果
　　试验说明：已知导航星在地固坐标系中的位置和速度、伪距和伪距率，求解用户星在地固系中的粗略位置和速度（初值）。时间跨度表示从跟踪到一颗GPS卫星(t=0）到开始解算（t=3τ)的时间间隔，所用采样点测量数据在时间上是等间隔的，为了不使方程组之间线性相关，两个采样点的时间间隔τ取45 s，测量数据的采样点个数为4，显然时间跨度为45 s×3=135 s。所给出的每组数据的顺序为：,位置单位为m，速度单位为m/s;误差数据亦与此相对应。
　　本文所用数据由GPS仿真器提供。由于本文所研究的一种解算方法，适用于所有将距离和距离变化率作为观测量的初轨计算过程，所以没有考虑SA、白噪声等的影响。这里的误差只是数学方法和轨道模型产生的误差。
　　t=0时刻用户星理论状态（地固坐标系）：
　　-2 316 583.253 000, 6 499 495.521 000, -70.306 491, 0.000 000, 
　　-6 685.420 898, -2 382.851 807，-0.350 270，0.000 000。
　　第1组初值设置：
　　-3 364 156.763 303，6 593 386.279 729，58 139.921 443，-51 733.977 767，
　　-6 683.351 845，-2 388.239 170,2.227 688，-6.959 785。
　　轨道计算结果：
　　-2 283 038.849 317，6 508 687.023 723，-63 209.623 812，12 968.581 227，
　　-6 695.727 940，-2 356.172 695，63.580 606，0.000 006。
　　相对于理论轨道的误差：
　　33 544.403 683，9 191.502 723，-63 139.317 320，12 968.581 227，-10.307 042，
　　26.679 112，63.930 876，0.000 006。
　　第2组初值设置：
　　-4 090 007.993 856，4 180 229.025 468， -2 276 755.318 757，147 809.209 786，
　　-7 750.189 630,-2 187.382 733，1 700.621 400，0.0。
　　最终的轨道计算结果：
　　-2 288 488.440 407，6 508 634.111 854，-59 223.153 740，14 533.638 719，
　　-6 693.054 309，-2 360.894 388，50.692 479，0.000 124。
　　相对于理论轨道的误差：
　　28 094.812 593，9 138.590 854,-59 152.847 249，14 533.638 719,-7.633 411,
　　21.957 419,51.042 749,0.000 124。
　　由实验结果可知，这种情况需要积累2 min 15 s时间的测量数据，初轨计算的位置误差约为60 km，速度误差约为50 m/s。
7　结　论
　　一般的初轨方法多是基于卫星轨道的先验知识和多弧段测量的。本文将Newton割线法和中点求积法（延拓法的一种)结合起来应用于低轨卫星的初轨计算。通过计算机仿真实验可以看出，这种方法的收敛范围很大，和基于线性变换的初值方法相结合［8］，可以在极恶劣的前提条件下（单个观测物体，较短的弧段，几乎不需要轨道的先验知识）计算出一定精度的初轨轨道。为进一步的轨道测控提供有用的信息。虽然本文以GPS系统为背景，但其算法适用于TDRSS系统、地面观测系统等只有一个观测物体的情况。由于这些系统没有钟差和频差，可以推测这种算法能提供更好的初轨计算结果；同时，这些系统的观测站数因远远少于GPS系统，所以，更需要这种适用于恶劣条件的算法。本文主要从理论上对这类问题进行了探讨，工程的应用需做更多的工作。
作者单位：洪炳熔　刘艳芳　(哈尔滨工业大学，哈尔滨 150001）
　　　　　胡芳　(北京航空航天大学，北京，100022)
参考文献
　1　Jin Linsheng.Preliminary orbit determination using satellite-to-satellite limited range and range-rate data.University　of Colorado at Boulder，Ph.D.Dissertation，1986.
　2　李济生.人造卫星精密轨道确定.北京：解放军出版社，1995.
　3　王广运.差分GPS定位技术与应用.北京：电子工业出版社，1996.
　4　J O Cappellari，et al.Mathematical theory of goddard trajectory determination system.1976.
　5　刘林.人造地球卫星轨道力学.北京：高等教育出版社，1992.
　6　李庆扬，莫孜中，祁力群.非线性方程组的数值解法.北京：科学出版社，1997:83～87
　7　H Wacker.Continuation methods.New York:Acedemic Press.1978:249～276
　8　刘艳芳，洪炳熔.用一颗GPS卫星确定低轨卫星轨道的初值方法.中国空间科学技术，1999(4).
