轨道计算
轨道计算
一种粗略测定天体轨道的方法。在轨道计算中﹐人们事先不必对天体轨道作任何初始估计﹐而是从若干观测资料出发﹐根据力学和几何条件定出天体的初始轨道﹐以便及时跟踪天体﹐或作为轨道改进的初值。为了计算六个轨道要素(见二体问题)﹐至少必须有三次光学观测﹐因为每次观测只能得到天体坐标的两个分量。
轨道计算是从研究彗星的运动开始的。在牛顿以前﹐对天体运动的研究基本上带有几何描述的性质。第谷首先试图计算彗星轨道﹐但未获成功。困难在于只能观测彗星的方向﹐而不知道它同地球的距离﹐由于缺少力学规律的指引﹐无法根据这些定向资料求得天体的空间轨道。在牛顿运动定律和万有引力定律发现螬o开普勒定律有了力学解释﹐得到了椭圆运动的严格数学表达式﹐终于能利用少数几次时间相隔不长的观测来测定彗星的轨道。
拉普拉斯方法 第一个正式的轨道计算方法是牛顿提出的。他根据三次观测的资料﹐用图解法求出天体的轨道。哈雷用这个方法分析了1337～1698年间出现的24颗彗星﹐发现1531年﹑1607年和1682年出现的彗星是同一颗彗星﹐它就是有名的哈雷彗星。在这以后﹐欧拉﹑朗伯和拉格朗日等人也在轨道计算方面做了不少研究。拉普拉斯于1780年发表第一个完整的轨道计算的分析方法。这个方法不限制观测的次数﹐首先根据几次观测﹐定出某一时刻天体在天球上的视位置(例如赤经﹑赤纬)及其一次﹑二次导数﹐然后从这六个量严格而又简单地求出此时天体的空间坐标和速度﹐从而定出圆锥曲线轨道的六个要素。这样﹐拉普拉斯就将轨道计算转化为一个微分方程的初值测定问题来处理。从分析观点来看这是一个好方法﹐然而轨道计算是一个实际问题﹐要考虑结果的精确和计算的方便。拉普拉斯方法在实用上不甚方便。由于数值微分会放大误差﹐这就需要用十分精确的观测资料才能求出合理的导数。尽管许多人曾取得一定进展﹐但终究由于计算繁复﹐在解决实际问题时还是很少使用。
奥伯斯方法和高斯方法[1]与拉普拉斯不同﹐奥伯斯和高斯则认为﹐如果能根据观测资料确定天体在两个不同时刻的空间位置﹐那么对应的轨道也就可以确定了。也就是说﹐奥伯斯和高斯把轨道计算转化为一个边值测定问题来处理。因此﹐问题的关键是如何根据三次定向观测来定出天体在空间的位置。这既要考虑轨道的几何特性﹐又要应用天体运动的力学定律。这些条件中最基本的一条是天体必须在通过太阳的平面上运动。由于从观测掌握了天体在三个时刻的视方向﹐一旦确定了轨道平面的取向﹐除个别特殊情况外﹐天体在三个时刻的空间位置也就确定了。轨道平面的正确取向的条件是所确定的三个空间位置能满足天体运动的力学定律﹐例如面积定律。
彗星轨道大都接近抛物线﹐所以在计算轨道时﹐常将它们作为抛物线处理。完整的抛物线轨道计算方法是奥伯斯于1797年提出的。他采用牛顿的假设﹐得到了彗星地心距的关系式﹔再结合表示天体在抛物线轨道上两个时刻的向径和弦关系的欧拉方程﹐求出彗星的地心距﹔从而求出彗星的抛物线轨道。到现在为止﹐奥伯斯方法虽有不少改进﹐但基本原理并没有变﹐仍然是一个常用的计算抛物线轨道的方法。
1801年1月1日﹐皮亚齐发现了第一号小行星(谷神星)﹐不久高斯就算出了它的椭圆轨道﹐他的方法发表于1809年。高斯使用逐次近似法﹐先求出天体向径所围成的扇形面积与三角形面积之比﹐然后利用力学条件求得天体应有的空间位置﹐再从空间位置求得轨道。高斯不仅从理论上﹑而且从实际上解决了轨道计算问题。可以说﹐用三次观测决定轨道的实际问题是高斯首先解决的。高斯以后﹐虽然有人提出一些新方法﹐但基本原理仍没有变。
人造卫星轨道计算 计算小行星轨道的经典方法﹐原则上都能用来计算人造卫星的轨道。在考虑到人造卫星的运动特点之后﹐又提出了一些新的方法。人造卫星运动快﹐周期短﹐记时误差对轨道计算结果影响显著。巴特拉科夫在高斯方法的基础上﹐用增加观测资料的办法﹐对记时有误差的轨道计算法作了改进。近地卫星一天绕地球飞行十多圈﹐容易从观测定准它的周期﹐因而也就知道了轨道半长径﹐相应地提出了已知半长径的轨道计算法。人造卫星离地球近﹐视差现象明显﹐利用两站或多站同步观测容易求得卫星地心距﹐可以简化经典计算方法。针对卫星摄动影响大的情况﹐又出现了考虑摄动的轨道计算法。尽管这些方法多种多样﹐仍不外乎从观测资料求得两个点的向径﹐或一个点的向径和速度﹐从而得到轨道要素。
通过对人造卫星激光测距和多普勒测速﹐利用多站同步观测﹐或结合光学观测等方法﹐可以直接得到卫星的向径和速度﹐从而求得卫星的轨道。应用高速电子计算机﹐可以进行复杂的迭代运算。因此﹐目前更多的是综合各种类型的观测资料作轨道改进﹐而不把精力放在初始轨道的计算上。现代技术条件已能使入轨后的卫星轨道同预定轨道相差不大。这样﹐预定轨道就能作为初始轨道使用。
参考书目
P.R. Escobal﹐Methods of Orbit Determination﹐J.Wiley and Sons﹐New York﹐1965.
A. D. Dubyago﹐ The Determination of Orbits﹐Macmillan Co.﹐New York﹐1961.
月球轨道的计算： 计算机可精确地求解星球运动微分方程。本算法适用于用计算机求解任意多体问题的星球运动微分方程，也能精确地确认万有引力公式本身和实际情况的相差程度。
一．计算机计算星球轨道原理。 以地球绕日运行为例，如图1-1。地球受力如图1-2。
为叙述方便，以下用矢量形式表达。
因为其中M太为太阳质量，M地为地球质量，G为万有引力恒量，a为地球加速度,为单位矢量。所以我们可按等时间间隔dt（即等步长），以微分形式从地球的初值点逐点向下推算。设t=0时，地球的初值点为r0,v0和于是，地球经dt时间从初值点到达第一点,递推式为,
由于dt是人为设定的，是已知的，因此地球到达1点的近似值v,r和a可由上式算出，算出1点值后，可把1点值作为初值，按步长dt继续推算出下一点的值，如此，可推算到第n点。由于dt值取得越小，递推的精度越高，我们可据此来控制计算误差。
设要计算地球在t=T时的r值，要求计算误差为e，t=0时的初值r0,a0,v0为已知。我们可将0到T的时间间隔划分为n个dt，即令计算步长dt=T/n,然后根据上述，按步长dt从t=0时的初值点推算到t=n·dt=T时的r值。然后将dt二分，即令计算步长dt1=dt/2,再按此新步长值dt1从t=0时的初值点算到t=2·n·dt1=T时的r值为r2，比较一下二分前后的r值，即看一看是否满足条件r2 - r<e,若满足条件，则r2的值即为所要求的r值，若不满足条件，则继续二分，按新的步长值再从t=0时的初值点算到t=T时刻，直到新算出的r值满足条件r2 - r<e，其中r和r2分别为二分前后的r值。注意，上面所说的矢量比较包含其大小比较和其方向比较。
以上为矢量表达，实际计算中，可将矢量v,r，a分别在x,y轴上投影，可得vx,vy,rx,ry,ax,ay。于是，它们的初始值分别为v0x,v0y,r0x,r0y,a0x,a0y。下面给出地球从初值点经dt时间运行到下一点的递推式，
控制误差范围的条件为分别为二分前后在x,y方向的r值。
以上仅为计算机计算地球轨道的原理。实际上，每二分一次，从0到T时间范围内的dt数量将增加一倍，计算机计算的工作量也将增加一倍。由于计算机的计算速度有限，因此二分次数也是有限的。为提高计算精度，减少计算机的计算工作量，有一些标准化的方法（注1），在此不再熬述。
由上述可知，计算机计算星球轨道主要有两个要点。一是列出递推式，二是确定误差范围的条件。
月球轨道计算见下页。
注释1：参见“计算机数值计算方法及程序设计”一书。该r书由周煦编著。于2004年10月由机械工业出版社出版。
二．月球轨道计算
由于地球的运动直接影响月球的运动，因此，先来分析一下地球的受力，如图1-3所示。
在图1-3中，o2x2y2z2坐标系是动坐标系，原点在地球中心。该坐标系跟随地球作平动，且三个坐标轴x2,y2,z2始终分别平行于x,y,z三个坐标轴。r1 是地球的位置矢量，r是月球的位置矢量, r2 是月球相对地球的位置矢量。
F月地是月球对地球的引力，F太地是太阳对地球的引力。设r1 与x,y,z轴的夹角分别为α1，β1，γ1，r与x,y,z轴的夹角分别为α，β，γ，r2 与x2,y2,z2轴的夹角分别为α2,β2,γ2,则，地球在x,y,z方向所受合力为：
因此，地球在x,y,z方向的加速度：
月球的受力如图1-4所示。月球在x,y,z方向所受合力为:
其中，F太月为太阳对月球的引力，F地月为地球对月球的引力。因此，月球的加速度为：
设a的初值为的初值为这样，地球和月球从各自的初值点同时出发，经dt时间后，地球就到达了它的下一点于是可得如下递推式：
(见下页)
控制计算误差的6个条件为：
其中分别为二分前后算出的地球坐标。再次说明一下，以上月球轨道的计算仅是计算机计算原理，实际编程应采取一些标准化方法，以提高计算精度，减少计算机的计算工作量。
目前，在月球轨道计算上，我已做到了，一天的计算误差e<0.001米（即在x,y,z轴方向的计算误差e）,也就是说一年的计算误差e<365×0.001=0.365米。要核实万有引力公式本身和实际情况的相差程度，可取两组实际观测值，一组观测值作为计算的初值，另一组观测值作核实之用，即核实用万有引力公式来计算的星球轨道的准确程度。下面采用一组实际观测值（注2）作为计算初值，让计算机来计算一下月球的轨道。初值为：
下面给出计算结果。
以上计算的时间范围是2006年。月球轨道半径最大和最小值都是指平均值。观测值由紫金山天文台的工作人员提供。以上计算取计算时间t=366天，坐标计算误差e< 0.366米。计算周期时，时间计算误差小于0.05秒。由于天文台提供的月球数据是相对地球坐标系的，地球坐标系和本文所述的动坐标系的关系是，将地球坐标系的yz平面以x轴为转轴旋转一个黄赤交角，就是本文所述的动坐标系。本文取黄赤交角为2326′。
注释2：作为计算初值的观测值的对应时刻为，2006年北京时间3月15日7时47.5分。该时刻恰好为半影月食的食甚时刻。为方便计算， 本文把该时刻定为零时刻。
下面给出从2006年3月月食食甚时刻计算到9月月食食甚时刻的地球和月球坐标。数据如下。
地球坐标（单位：米）