PMOE 2003行星历表框架(I)数学模型.pdf

紫金山天文台台刊 PUBLICATIONS OF PURPLE MOUNTAIN OBSERVATORY 2007 1 18 1 PMOE 2003 行星历表框架行星历表框架 I 数学模型 数学模型 李广宇 1 2 倪维斗1 2 田兰兰1 2 1 中国科学院紫金山天文台 210008 2 中国科学院国家天文台 100012 摘要 输入历表研究是单飞行器激光天文动力学空间计划预研究工 作的一项重要内容 我们在张家祥太阳系动力学模型和汤健仁 倪维斗等 CGC 2 星历框架的基础上建立了 PMOE 2003 精密行星历表框架 本框架对 天体后牛顿效应 日月地有限体形状摄动 地球潮汐效应等问题的处理 均较前有大的改进 预报天体位置的精度与 DE 405 历表接近 本文主要 介绍框架涉及的数学模型 数值积分方法和程序结构及文本等 将另文给 出 关键词 行星历表 太阳系动力学 1 引言 建立独立的精密太阳系行星历表是单飞行器激光天文动力学 ASTROD I or Mini ASTROD One Spacecraft Astrodynamical Space Test of Relativity using Optical Devices 和 双 飞 行 器 激 光 天 文 动 力 学 ASTROD II Two Spacecraft Astrodynamical Space Test of Relativity using Optical Devices 空间任务预研究的一项重要工作 独立的具有可调变参数的行星历表不仅是处理本 空间任务行星际激光测距数据 获取和改进各天文物理常数必不可少的工作平台 而且对于本空间任务飞行器轨道的设计也是至关重要的 1 1960 年代以前 太阳系行星历表是以 Newcomb 的分析理论为基础的 60 年代的 两件事情改变了这种状况 一是电子计算机和计算技术的迅速发展 一是太阳系天 体雷达测距观测的实现 这两方面的进展使得发展新的以运动方程数值积分为基本 方法的精密行星历表成为可能 JPL 的 E M Standish 等在此基础上推出了 DE 系 列行星历表 其中影响较大的有 1975 年的 DE 96 1977 年的 DE 102 2 1982 年的 DE 200 3 1995 年的 DE 403 和 1998 年的 DE 405 4 DE 200 历表是拟合处理 50 000 个以上位置观测数据的结果 在拟合中考虑了广义相对论时延 太阳日冕和对流层 电子密度等因素引起的改正 根据 IAU 1976 1979 和 1982 年大会的决议 从 1984 年起 DE 200 数值历表取代 Newcomb 的分析理论 成为世界各国编算天文年历的基 础 DE 403 和 DE 405 历表进一步改进了数学模型 增加了更多的观测数据 DE 405 还考虑了影响较大的 300 颗小行星的摄动作用 成为当今通用的行星历表 5 DE 405 历表以 32 天一组 Chebychev 多项式的形式给出了 1600 AD 2200 AD 期 间太阳 大行星和月球的位置速度 以及地球章动和月球天平动的数据 还以论文 形式部分公布了数学模型和数值积分的初始数据 由于该历表所用的参数对于我们 紫金山天文台台刊 PUBLICATIONS OF PURPLE MOUNTAIN OBSERVATORY 2007 1 18 2 是不可调变的 不能用作激光天文动力学空间任务飞行器轨道的设计的平台 同时 尽管 DE 405 历表就现今行星际雷达测距数据的拟合而言是最精密的 但激光天文动 力学空间任务所要实现的行星际激光测距的精度将要比雷达测距提高几个数量级 要对这些数据进行处理 还必须考虑 DE 405 历表的数学模型中没有考虑到的数量级 更小的摄动因素 例如二阶后牛顿效应 因此在DE 405历表的基础上发展新的历表 是激光天文动力学空间任务预研究的重要一环 1 在 DE 历表之外 国外还有一些重要的历表工作 法国经度局长期进行分析历表 研究 发展了半分析半数值历表 VSOP 6 俄国科学院应用天体物理研究所的 E V Pitjeva 等发展了与 DE 403 相当的 EPM 1998 历表和与 DE 405 相当的 EPM 2000 历 表 并在此基础上进行了确定天文常数的研究 最近 他们更精密地考虑小行星环 带的摄动 拟合了新的测距数据 发表了 EPM 2002 历表 其精度已较 DE 405 略好 7 9 新竹清华大学重力实验室的倪维斗和他的研究生邱大维 汤健仁 张中豪等发 展了 CGC 1 和 CGC 2 历表框架 两者的精度分别与 DE 403 和 DE 405 历表相当 他 们在此基础上进行了激光天文动力学空间任务飞行器轨道的初步设计和天文常数的 模拟计算 10 13 紫金山天文台张家祥为了研究太阳系小天体的运动 在 1970 年代开始太阳系动 力学模型的研究 他和他的学生李广宇用改进的 Cowell 积分方法研究了紫金山彗星 14 哈雷彗星 15 16 小行星伊卡鲁斯 Icarus 等 17 18 的长期轨道演化 该模型给出 的太阳系大行星的位置 在100年内 误差不大于1角秒 内符合计算精度高于10 12 19 1984 年 张家祥用自己建立的历表预报彗木碰撞事件中各个碎块撞击木星的时刻 与伽利略飞船测得的实际碰撞时刻比较 均方差为 8 46 分钟 精度与 JPL 的预报相 当 20 李广宇和赵海斌还用同一历表研究了 1998 年狮子座流星雨第二峰粒子的来源 问题 21 2001 年我们为适应太阳系空间探测 特别是单飞行器激光天文动力学 ASTORD I 空间计划预研究的需要 在张家祥的太阳系动力学模型 19 和 CGC 2 历表框架 13 的基础上 建立了 PMOE 2002 历表框架 22 该框架预报大行星的位置与 DE 405 历表比较 从历元时刻 JD 2440400 5 1969 年 6 月 28 0 日 开始 向后积 分 1 200 天和 12 000 天 约 33 年 以后的差值分别如表 1 表 2 所示 框架所用 天体初始数据和物理 天文常数都取自 DE405 历表 积分 1 200 天后 内行星日心距之差在 20 米内 日心经度和日心纬度之差在 0 5 毫角秒内 PMOE 2002 历表框架预报月球位置的精度较差 积分 1 200 天后 地心距之差超过 150 米 地心经度之差超过 1 700 毫角秒 图 1 地球的日心位置 也相应变差 取不同的积分起始时刻 这些差值会有同数量级的变化 为了确定起 见 我们约定本文中所讨论的差值的起始时刻 全部为 JD 2440400 5 PMOE 2002 历 表框架预报月球位置精度较差的原因如下 i 对于天体之间的后牛顿效应 只考虑了太阳质量对其他天体的影响 没有 考虑其他天体质量的影响 特别是太阳 月亮和地球三者之间的影响 ii 只考虑了地球有限体的二次带谐 J2 项对月球的影响 没有考虑其他高次 带谐项和田谐项的影响 iii 未考虑月球有限体及月球天平动对地球和太阳的影响 iv 未考虑地球潮汐对月球运动的影响 2002 年我们针对以上各项对框架进行了改进 发展为 PMOE 2003 历表框架 该 框架预报大行星的位置与 DE 405 比较 从历元时刻 JD 2440400 5 开始 向后积分 紫金山天文台台刊 PUBLICATIONS OF PURPLE MOUNTAIN OBSERVATORY 2007 1 18 3 1 200 天和 36 000 天 约 98 年 以后的差值分别如表 3 表 4 所示 积分 1 200 天 后 月球的地心距之差减小至 2 2 米 地心经度之差减小至 5 4 毫角秒 图 2 各 天体积分 36 000 天时差值的变化见本文附录所列各图 本文以PMOE 2003框架为主介绍建立这两个框架所用的数学模型 关于数值积分 方法 程序设计要点 以及做为模型的初步检验 应用 Kalman 滤波方法拟合火星中 天观测数据和 Viking 火星雷达测距数据的方法和初步结果 将另文介绍 表 1 积分 1 200 天 PMOE 2002 与 DE 405 的差值 天体 日心距 r 米 日心经度 毫角秒 日心纬度 毫角秒 水星 20 0 0 5 0 2 金星 18 15 0 0 2 0 1 地月质心 15 5 0 0 2 0 1 火星 79 26 0 0 3 0 1 木星 552 0 0 0 2 0 02 土星 124 0 0 006 0 001 0 001 0 天王星 169 0 0 0 002 0 001 0 海王星 94 3 0 0 004 0 001 0 冥王星 131 0 0 001 0 0 001 表 2 积分 12 000 天 PMOE 2002 与 DE 405 的差值 天体 日心距 r 米 日心经度 毫角秒 日心纬度 毫角秒 水星 190 0 5 7 1 6 2 1 金星 25 0 1 7 0 6 地月质心 40 32 0 2 6 0 9 火星 512 421 0 5 3 1 5 2 1 木星 2388 1686 0 5 1 1 0 1 7 土星 1612 0 0 1 9 0 2 0 6 天王星 3687 0 0 0 5 0 0 2 海王星 1682 3 0 0 1 0 0 02 冥王星 2486 0 0 0 1 0 04 0 表 3 积分 1 200 天 PMOE 2003 与 DE 405 的差值 天体 日心距 r 米 日心经度 毫角秒 日心纬度 毫角秒 水星 4 9 4 8 0 12 0 0 10 0 12 金星 12 4 11 6 0 04 0 05 0 01 0 03 地球 2 4 2 2 0 0 01 0 01 0 01 月球 2 2 2 2 1 35 5 39 1 66 1 78 地月质心 2 3 2 2 0 0 01 0 01 0 01 火星 27 8 18 9 0 03 0 07 0 0 02 木星 482 4 0 0 0 18 0 01 0 02 土星 132 3 0 0 0 01 7 7 10 8 0 天王星 15 1 0 0 0003 0 0 00002 0 00009 海王星 0 18 0 0 0 0003 0 00005 0 紫金山天文台台刊 PUBLICATIONS OF PURPLE MOUNTAIN OBSERVATORY 2007 1 18 4 冥王星 0 20 3 0 0004 0 0 0 0003 月球为地心距 地心经度和纬度 表 4 积分 36 000 天 PMOE 2003 与 DE 405 的差值 天体 日心距 r 米 日心经度 毫角秒 日心纬度 毫角秒 水星 155 1 157 4 6 5 0 2 6 1 9 金星 22 6 22 4 0 4 0 1 0 4 0 4 地球 24 0 24 9 0 0 2 0 2 0 2 月球 64 3 65 7 663 3 4 8 272 4 265 3 地月质心 18 4 18 4 0 0 2 0 2 0 2 火星 1223 1270 0 9 4 2 3 3 6 木星 3570 2128 0 16 8 5 4 6 0 土星 950 5 81 20 3 4 1 2 1 1 天王星 657 5 447 9 0 1 0 0 0 1 海王星 0 7539 1 7 0 0 3 0 冥王星 0 12820 0 9 0 0 3 0 1 月球为地心距 地心经度和纬度 图 1 PMOE 2002 框架积分 1 200 日后 月球的地心距离 经度与 DE 405 的比较 距离的单位 为米 经度的单位为毫角秒 横坐标为太阳系质心力学时 TDB 单位为日 图 2 PMOE 2003 框架积分 1 200 日后 月球的地心距离 经度与 DE 405 的比较 距离的单位 为米 经度的单位为毫角秒 横坐标为太阳系质心力学时 TDB 单位为日 2 PMOE 历表框架 在任一惯性系中 太阳系天体的牛顿运动方程可以表达为 23 紫金山天文台台刊 PUBLICATIONS OF PURPLE MOUNTAIN OBSERVATORY 2007 1 18 5 n 0 1 i 3 ij ij ij j i r r r 用于存放右函数 f 在 t t0 kh k 4 3 4 共 9 个步点上的值 Pascal 的 extended 类型与 C 的 long double 类型相同 为 10 字节浮点数 表示为十进制实 数时的有效数字为 19 20 位 另一个是数组 sf sf array 1 2 1 2 of extended 用于存放连续两个二次和分 每个和分又由两个 extended 类型的实数分量构成 和 分有效数字的长度对于减少积分过程中误差的累积起着重要作用 如果采用以上数 据结构 对于和分乘以一个因数后分离整数部分和小数部分 分别存入两个分量 即 可以增加有效数字的长度 在有和分参加的加减运算中 注意避免小量和大量的直 接加减 先分离大量的整数部分和小数部分 然后再以其小数部分和小量加减 由于 数值积分每一步外推生成和分的过程中 总要将做为小量的函数值的和数加到原来 的和分上去 直接相加会损失小量的有效位数 产生不容忽视的误差累积 本积分器取步长为 0 125 日 外推 288 000 步后 水星坐标的最大累积误差不 超过 5 10 15 AU 0 08 cm 月球坐标的最大累积误差不超过 2 10 14 AU 0 3 cm 38 累积误差主要有三个来源 第一个来源是上述外推生成和分的过程中 由于小量 右 函数 和大量 和分 相加而产生的舍入误差 根据数值积分的理论 外推步数为 n 时 随机累积误差的大小约为 0 1124 n 3 2 36 单位同末位有效数字 取步长为 0 25 日 外推 144 000 步 相当于 98 年 后 随机误差约为 6142003 即误差向前传递 7 位数字 如步长为 0 125 日 积分同样的时间 误差向前传递 8 位数字 10 字节 浮点数的有效位数为十进制的 19 位 尚不足避免此种误差累积 引入上述数据结构 并取因数为 2 20时 和分的有效数字增加到十进制的 26 位 可以完全避免因字长不 够而引起的舍入误差的积累 第二个来源是所用近似积分公式的截断误差 我们的 积分器中 每一步积分坐标时的误差正比于右函数 f 的8阶差分与步长8次幂之积 紫金山天文台台刊 PUBLICATIONS OF PURPLE MOUNTAIN OBSERVATORY 2007 1 18 16 选取适当的步长 即可将这一部分误差的累积减至最小 最后一个来源是计算右函 数过程中的舍入误差 减小这部分误差的途径是增加积分初始数据的有效位数和程 序变量的字长 我们现在使用 10 字节浮点数变量和 16 位有效数字初值 即将总的 累积误差控制在了前述范围之内 积分器的基本结构用 Warnier 图 37 表示如下 f f 1 1 计算和分 数组函数值送 计算加速度 计算状态 计算时间 计算函数表 外推 外推步数置初值 迭代逼近 数组函数值送 计算加速度 按无摄运动计算状态 置时间 处理一个步点生成函数表 开头 N N 数值积分由开头和外推两步操作完成 开头操作的目的是生成开始 9 个步点上 的函数值表 f 对于天体 先以无摄运动求得 f 的近似值 对于月球的 Euler 角 先以 Euler 法求得 f 的近似值 然后迭代逼近得到精确值 每一次迭代的处理 与 下面外推的处理基本相同 外推是一个循环过程 每一次循环均先预估计算坐标速 度 更新函数值表 再用校正计算坐标速度 最后更新和分表 各天体的初始状态和有关天文物理常数 均尽可能采用 DE 405 历表的相同数据 4 DE 405 历表文献没有给出的 则取自其他最近有关文献 积分步长一般为 0 25 日 在此步长下 对于运动最快的天体月球 每个轨道周期外推约 120 步 仍能保证 积分精度 本框架的程序由基本运算 轨道计算和主控程序等三个主要单元组成 有 Object Pascal 39 和 C 程序设计语言 40 两种版本 分别在 Delphi 程序设计平台 和 C Builder 程序设计平台上调试完成 致谢 致谢 作者衷心感谢国家天文台优秀空间天文项目基金 紫金山天文台小行星基金 会和紫金山天文台台长基金等方面经费的支持 作者还衷心感谢紫金山天文台张家 祥研究员和南京大学黄天衣教授的关心与支持 参考文献 1 W T Ni et al Mini ASTROD Mission Concept Int J Mod Phys D11 7 1035 2002 W T Ni ASTROD An Overview Int J Mod Phys D11 7 947 2002 2 X X Newhall E M Standish Jr and J G Williams DE102 a numerically integrated ephemeris of the Moon and planets spanning forty four centuries Astron Astrophys 125 150 1983 3 E M Standish Jr The observational basis for JPL s DE200 the planetary ephemerides of the Astronomical Almanac Astron Astrophys 233 252 1990 紫金山天文台台刊 PUBLICATIONS OF PURPLE MOUNTAIN OBSERVATORY 2007 1 18 17 4 E M Standish JPL Planetary and Lunar Ephemerides DE 405 LE405 JPL IOM 312 F 98 048 1998 DE 403 LE403 JPL IOM 314 10 127 1995 5 E M Standish and A Fienga Accuracy limit of modern ephemerides imposed by the uncertainties in asteroid masses A A 384 322 328 2002 6 P Bretagnon Analytical Planetary Solution VSOP 2000 Celest Mech Dyn Ast 80 3 205 2001 7 E V Pitjeva Progress in the determination of some astronomical constants from radiometric observations of planets and spacecraft A A 371 760 2001 8 E V Pijeva EPM2002 and EPM2002C two versions of high accuracy numerical planetary ephemerides constructed for TDB and TCB time scales Tpy N 10 2003 9 E V Pijeva Modern numerical ephemerides of planets and importance of ranging observations for their creation Celest Mech Dyn Ast 80 N3 4 249 271 2001 10 D W Chiou and W T Ni ASTROD orbit simulation and accuracy of relativistic parameter determination Adv Space Res 25 1259 1262 2000 11 汤健仁 小行星微扰 星历与 ASTROD 太空任务概念测定小行星质量之模拟 新竹清华大学物理研究所硕士论文 2000 12 张中豪 小型雷射太空力学任务概念之轨道设计与模拟 新竹清华大学物理研 究所硕士论文 2001 13 汤健仁 倪维斗 云南天文台台刊 3 21 2002 14 张钰哲 张家祥等 天文学报 15 1 1974 15 张钰哲 天文学报 19 1 109 1978 16 李广宇 紫金山天文台台刊 1 2 1982 17 李广宇 紫金山天文台台刊 2 2 1983 18 张家祥等 天文学报 29 3 1988 19 张家祥 紫金山天文台台刊 8 3 1989 20 J X Zhang Q Wang J X Yang S C Wang X Z Chen Science in China Series A 39 2 207 1996 21 G Y Li H B Zhao Int J Mod Phys D11 7 1021 2002 22 李广宇 倪维斗 田兰兰 PMOE 精密行星历表研究进展 天文学会第二届学术 年会 南京 2002年11月 紫金山天文台台刊 2003 23 易照华 天体力学基础 南京大学出版社 1993 24 Dennis D McCarthy IERS Technical Note 21 IERS Conventions 1996 25 V A Brumberg Essential Relativistic Celestial Mechenics Adam Hilger Bristol England 1991 26 D H Eckhardt Theory of the Libration of the Moon The Moon and the Planets 25 3 49 1981 27 刘林 天体力学方法 南京大学出版社 1998 28 刘林 航天器轨道理论 国防工业出版社 2000 29 Allen s Astrophysical Quantities 4th edition p 362 ed A N Cox AIP 2000 30 夏一飞 黄天衣 球面天文学 南京大学出版社 南京 1995 31 紫金山天文台 中国天文年历 科学出版社 北京 2002 32 X X Newhall and J G Williams Estimation of the Lunar Physical Librations Celes Mech 66 21 1997 紫金山天文台台刊 PUBLICATIONS OF PURPLE MOUNTAIN OBSERVATORY 2007 1 18 18 33 G A Krasinsky Tidal Effects in the Earth Moon System and the Earth s Rotation Celes Mech 75 39 1999 34 G A Krasinsky Dynamical History of the Earth Moon System Celes Mech 84 27 2002 35 E Bowell The Asteroid Orbital Elements Database DB ftp ftp lowell edu pub elgb astorb html 36 D Brouwer and G M Clemence Methods of Celestial Mechanics Academic Press New York and London 1961 37 王选 软件设计方法 清华大学出版社 北京 1992 38 李广宇 田兰兰 PMOE精密行星历表框架 II 积分器和程序设计 紫金山 天文台台刊 22 3 4 2003 39 Steve Teixeira Xavier Pacheco Delphi 5 Developer s Guide Sams Publishing U S A 2000 Delphi 5开发人员指南 任旭钧等译 机械工业出版社 北京 2000 40 Stanley B Lippman Josee Lajoie C Primer Third edition Addison Wesley Longman Inc U S A 2000 Sams Publishing U S A 1998 潘爱民等译 中国电力出版社 北京 2002 附录1 PMOE 2003预报行星和月球位置与DE 405的比较 积分起始时刻为 JD 2440400 5 1969 年 6 月 28 0 日 紫金山天文台台刊 PUBLICATIONS OF PURPLE MOUNTAIN OBSERVATORY 2007 1 18 19 图 12 PMOE 2003 框架积分 1 200 日后 月球的地心轨道根数与 DE 405 的比较 半长径的 单位为米 角度根数的单位为毫角秒 横坐标为太阳系质心力学时 TDB 单位为日 图 13 PMOE2003 框架积分 36 000 日后 月球的地心距离 经度与 DE 405 的比较 距离的单 位为米 经度的单位为毫角秒 横坐标为太阳系质心力学时 TDB 单位为日 图 14 PMO E2003 框架积分 36 000 日后 地球的日心距离 经度与 DE 405 的比较 距离的 单位为米 经度的单位为毫角秒 横坐标为太阳系质心力学时 TDB 单位为日 图 15 PMOE 2003 框架积分 36 000 日后 水星的日心距离 经度与 DE 405 的比较 距离的 单位为米 经度的单位为毫角秒 横坐标为太阳系质心力学时 TDB 单位为日 紫金山天文台台刊 PUBLICATIONS OF PURPLE MOUNTAIN OBSERVATORY 2007 1 18 20 图 16 PMOE 2003 框架积分 36 000 日后 火星的日心距离 经度与 DE 405 的比较 距离的 单位为米 经度的单位为毫角秒 横坐标为太阳系质心力学时 TDB 单位为日 图 17 PMOE 2003 框架积分 36 000 日后 木星的日心距离 经度与 DE 405 的比较 距离的 单位为米 经度的单位为毫角秒 横坐标为太阳系质心力学时 TDB 单位为日 图 18 PMOE 2003 框架积分 36 000 日后 海王星的日心距离 经度与 DE 405 的比较 距离 的单位为米 经度的单位为毫角秒 横坐标为太阳系质心力学时 TDB 单位为日 附录 2 潮汐摄动加速度的推导 表达式 29 给出了月球由于地球潮汐摄动获得的加速度 该式推导稍长 但 技巧可资借鉴 因此以附录的形式给出 由 2 5 节 设月球向径为 mmmm zyx r 则试验点 zyx r处的附 加势为 5 05 1 2 22 5 5 5 2 rr rr a V mm m m m rrrr 27 紫金山天文台台刊 PUBLICATIONS OF PURPLE MOUNTAIN OBSERVATORY 2007 1 18 21 对x微分之 mm m m m m xx r r r r a x V rr rr 35 1 5 71 2 2 2 5 5 2 31 对y和z的微分有同样的形式 因此 mmm m m mr r rr a V zyx V rrrr rr 35 1 5 71 2 2 2 5 5 2 32 我们要求取的是月球因潮汐作用获得的加速度 因此试验点向径即是月球在时 刻t的向径 r r 考虑到潮汐隆起的时间延迟 m r则应是月球在时刻t 的向径 但要修正到时刻t隆起实际发生时所在的方向 r rrrr rrr tt m 28 准确到 r r 的一阶小量时有 rrrrr r r r rrr r rrrr rr rr rr r mm m m m r r r r rO r rr 22 24 2 2 5 5 2 55 2 2 5 11 33 代入 32 式 并略去 r r 的高阶小量得 rr rr r rr rr r rr 23 5 2 223 5 2 21 3 35 1 1 3 rrr a rrrr a V m m 34 由 28 式有 rrrr rrr 35 代入 34 式即得所证 紫金山天文台台刊 PUBLICATIONS OF PURPLE MOUNTAIN OBSERVATORY 2007 1 18 22 r rr rrr rrr rr r r