基于MatlabR2011b的卫星轨道计算

X姨
姨
k
姨 姨
姨姨xkcosΩk-ykcosiksinΩk
姨 姨
姨 姨姨
姨
Y姨
姨 姨
k
=姨
姨 姨
姨姨姨xksinΩk+ykcosikcosΩk
。 姨
姨 姨
Z姨姨 姨姨 姨姨
姨 k姨 姨
yksinik
姨姨 姨
輥輰訛卫星在协议地球坐标系中的位置：
X姨 姨
姨姨
X 1 姨
姨
k
姨 姨
姨 姨
姨姨
姨 姨姨
Y来自百度文库 姨
测绘工程。
设以O为原点的直角坐标系为O-XYZ，卫星的坐标为
（X，Y，Z），则卫星的地心向量r（X，Y，Z），加速度r（X，Y，
Z），代入上式得：
r=g（a，e，I，ω，τ，t）
→ →
→→
dr dt
→
=g（′ a，e，I，ω，τ，t）→→ → →
从上式可以看出，在二体问题情况下，给定6个轨道
参数，即可确定任意时刻t的卫星位置及其运动速度。
① rinexe（ephemerisfile，outputfile）。ephemerisfile—n 文件名 （n文件放在rinexe函数所在目录下）；outputfile— 文本文件名，将读到n文件的数据存放在该文本文件中。
② eph =get_eph（ephemeridesfile）。ephemeridesfile - 文 件名，将存放在该文本中的导航电文数据提取出来，为计 算卫星位置做准备；返回值—eph星历矩阵，将读到的数 据存入内存。 3.2.3 计算卫星位置模块
牛晓楠，蒋丹丹，霍红娟
（中国矿业大学 环境与测绘学院，江苏 徐州 221116）
摘 要：文章阐述了计算卫星在地心坐标系中的位置，利用开普勒方程，根据广播星历的 21 个参数计算出开普勒轨道参
数，再结合卫星运动学方程计算出卫星的坐标。在此基础上，利用 MALAB R2011b 软件读取卫星导航（O，N）文件，在 MALAB
姨姨
=RY(-xP)R（Y -yP）=
Y姨
姨k
=姨
姨
姨 姨
0
姨姨
姨 姨姨
Z姨姨
姨
姨姨姨CTS
Z -X 姨姨 姨姨 姨姨
姨 k姨 姨 P
3.2 重要MATLAB函数说明
0
X X P
姨姨 姨姨
k
姨 姨
姨姨 姨
1
-Y Y 姨 姨 姨 P 姨姨 k 姨
姨姨 姨
XP
1
Z 姨姨 姨姨 姨姨
姨姨 k 姨
3.2.1 时间转换模块
k2M a3
=
姨μ ， 姨a 3
n=n0+△n。
②计算归化时间：t=t′-△t，tk=t-toe。
③观测时刻卫星平近点角MK的计算：MK=MO+ntk。
④计算偏近点角EK，EK=MK+esinE（K EK，MK以弧度计算）。
⑤真近点角fk的计算：fk=arctan
η ξ
=arctan
姨1-e2 sinEk cosEk-e
主要用到两个函数，分别如下： ① jd =julday（year，month，day，hour）year—年 ，四 位 ，
如2011；month—月，1-12，day—日；hour—时间，小时+不
足一小时部分（用小数表示）返回值-儒略日。
② [week，sow]=gps_time （jd）jd—儒略日；返回值— week，GPS周；sow，GPS秒。 3.2.2 读星历文件模块
则a=as-ae=-G（M+m）/r2
上式即为卫星S相对于地球质心O的运动方程。
即卫星相对于地球质心的加速度如下：
r=re-rs=-G（M
+m）
→r r3
。
地球质量远远大于卫星质量，通常略去卫星质量m项
r=re-rs=-GM
→r r3
。
收稿日期：2012 - 05 - 13 作者简介：牛晓楠（1990—），女，河北石家庄人，大学本科，研究方向：
② col_Eph（t）=find_eph（Eph，sats（t），time）。Eph—星 历矩阵，存放计算卫星位置所需的n文件数据；sats（t）— sats中存放o文件中某个历元观测到的卫星PRN号，t为循 环 控 制 ，sats（t）为 sats 中 某 颗 卫 星 ；time—o 文 件 中 某 个 历 元的历元时刻；返回值-col_Eph矩阵存放从Eph星历矩阵 中选出的用于计算的列数，如col_Eph=[1，2]，代表Eph中 第1，第2列用于计算卫星位置。
R2011b 环境中编写计算卫星坐标的程序过程及其重要函数说明和重要代码标注。
关键词：开普勒轨道参数；卫星运动力学方程；卫星轨道计算；程序语言
中图分类号：V412.41
文献标识码：A
文章编号：1006- 8937（2012）22- 0020- 02
Based on the satellite orbit calculation Matlab R2011b NIU Xiao- nan，JIANG Dan- dan，HUO Hong- juan
③X=satpos（tx_GPS，Eph（：，k））。tx_GPS—上节所述的 归化时间，用儒略日表示的；Eph（：，k）—Eph星历矩阵中 的某一列数据；返回值—卫星在地心地固坐标系中坐标。 4 程序的实现过程
程 序 的 开 发 基 于 MATLAB R2011b 平 台 实 现 的 。 MATLAB是矩阵实验室（Matrix Laboratory）的简称，是美 国Math Works公司出品的商业数学软件，用于算法开发、 数据可视化、数据分析以及数值计算的高级技术计算语 言和交互式环境。MATLAB R2011b界面，如图1所示。
GPS卫星导航与定位技术在测量工作中得到越来越 广泛的应用。在利用GPS信号进行导航定位时，为计算用 户在地心坐标系中的位置，GPS接收机需要测定测站到 卫星的距离并且要知道同一卫星在同一时刻的地心坐 标。卫星的地心坐标是从卫星的导航电文中提供的开普 勒轨道参数和轨道摄动修正量按一定的公式计算的。本 文介绍如何利用matlab软件实现卫星坐标的计算。
ξ
cosE-e
2.3 卫星坐标计算
GPS卫星在椭圆轨道上，以3.87 km/s的速度每天旋转
两圈，k卫星在tj时刻的地心坐标可表示为：
姨姨X（K t）j 姨姨
姨姨rkj cosfkj
姨 姨
姨
姨
姨
姨
姨姨姨Y（K t）j 姨姨姨=R（3 -Ωk）j R（3 -ik）j R（3 -ωk）j 姨姨姨rkj sinfkj
。
⑥升交距角Φk的计算：Φk=Vk+ω。
⑦摄动改正项δu，δr，δi的计算。
⑧计算经过摄动改正的升交距角uk、卫星矢径rk和轨
道倾角ik。
姨 ⑨计算卫星在轨道平面坐标系的坐标 xk=rkcosuk。 yk=rksinuk
⑩计算观测时刻升交点经度 Ωk：Ωk=ΩO+（Ω-ωe）tk-
ωetoe。
輥輯訛计算卫星在地心固定坐标系中的直角坐标：
①读观测值文件（*.o文件）函数。在求解卫星位置时， 第一需要利用o文件中每个历元的历元时刻t。在计算某时 刻卫星位置时，这里的某时刻便是o文件历元时刻t。第二 需要利用读取的每个历元不同的卫星PRN号。根据PRN号 和历元时刻t在广播星历n文件中搜索相同的卫星PRN号、 合适的历元时刻，利用其对应的数据，计算卫星位置。
1 卫星运动力学方程
研究卫星S绕地球O的运动，主要是研究卫星运动状
态 随 时 间 的 变 化 规 律 。根 据 牛 顿 定 律 ，可 以 方 便 地 得 到
二体问题的卫星运动方程。
用FS、Fe表示 卫星 与地 球所 受到的 引力 作用 力，则 FS=-GmM/r2，式中，G为万有引力常数。 Fe=GmM/r2
（School of Environment Science and Spatial Informatics，China University of Mining&Technology，Xuzhou，Jiangsu 221116，China）
Abstract:This paper expounds how to calculate the position of the satellite in the geocentric coordinate system：using Kepler equation， calculating Kepler orbit parameters according to 21 parameters of broadcast ephemeris，then calculate coordinate of satellite combining satellites kinematics equations.On these basis，reading satellite navigation（O，N）files using MATLAB R2011b software，and then design the program of calculating the position of satellite and some instructions of important function and important code label in the environment of matlab R2011b. Keywords：kepler orbit parameters；satellites kinematics equations；satellite orbit calculation；programming language
出编码器需采集的数据M=Yθ/360。
以取料量80 t/h，物料密度ρ=0.8～1.0 t/m3，行走速度
υ=2 m/min，电动葫芦直径D=510.5 mm，编码器线数Y=900，
得出表1数据。
表1 取料机每层下降数值表
共50层
0～10
姨 姨 姨
Z（t） 姨姨
姨
K
j
姨姨 姨
0姨姨
姨
姨姨 姨
其中，Ωkj ，ikj ，ωkj 分别为k卫星tj时刻的升交点赤径，卫
星轨道面倾角，近地点角距。R（3 -Ωk）j ，R（3 -ik）j ，R（3 -ωk）j 为
旋转矩阵。
3 卫星坐标计算算法步骤
3.1 算 法
姨 ①计算卫星运行的平均角速度：nn0=
DOI:10.14165/j.cnki.hunansci.2012.22.001
第 31 卷第 22 期 Vol.31 No.22
企业技术开发 TECHNOLOGICA企L D业EVE技LO术PME开NT发OF ENTERPRISE
2012 年 8 月 20A1u2g年.20812月
基于 Matlab R2011b 的卫星轨道计算
利用绘图法确定l1、l2、l3在图面2的位置，确定交接点，利用
CAD绘图工具标注出交接点到定滑轮的距离Hn，计算出
每一层钢丝绳的长度差δl，根据交点到卷扬滚筒钢丝绳的
倍率N，算出钢丝绳的释放长度L，用公式n=L/πD、θ=n/360
得 出 每 层 悬 臂 下 降 角 度 θ，已 知 编 码 器 线 数 Y，就 可 以 算
2 卫星轨道计算
2.1 卫星轨道方程
以卫星轨道平面为基准面的三维坐标系统分别为：
o-ξ、η、ζ。
以下是以开普勒第一定律为基础推导卫星相对地
球的位置表达式结果：
①由开普勒第一定律，卫星的运行轨迹为椭圆，在p-
o-q坐标内有：pa22
+
q2 b2
=1；
②卫星S的位置在ξ-c-η坐标系中（ξ，η）为ξ＝rcosf=a
姨 姨 姨
第 31 卷第 22 期
牛晓楠，等：基于 Matlab R2011b 的卫星轨道计算
21
⑤卫星与地球质心的几何距离为：‖r‖=a（1-ecosE）。 2.2 开普勒方程
姨 开普勒第三定律可表示为：n=
GM a3
。
真近点角f可表示为：
f=arctan η =arctan 姨1-e2 sinE 。
图1 MATLAB R2011b界面图
5 结语 文章利用开普勒定律，卫星运动学 （下转第 36 页）
36
企业技术开发
2012 年 8 月
M=Yθ/360（编码器记数值M，编码器工作线数Y）
取料量Q是恒定的，物料密度ρ和行走速度υ是定值，
取料截面S保持恒定就能满足上述要求。如图1取料截面
所 示 ，S1=h1l1/2，S2=h2l2/2，S3=h3l3/2 … Sn=hnln/2。S1=S2=S3… Sn
cos E-ae=a（cosE-e）；
③可以得到：η=q=b sin E=a 姨1-e2 sinE；
④则在ξ-η-ζ坐标系统下，卫星相对地球质心的距离
ξ姨 姨
姨姨
姨姨a(cosE-e)
姨
姨 姨 姨
姨姨姨
姨
向量r=
η姨
姨 姨
=姨
姨 姨
a姨
姨
姨1-e2
姨
sinE
； 姨
姨 姨
ζ 0 姨姨 姨姨 姨
姨姨姨 姨