基于MatlabR2011b的卫星轨道计算

DOI:10.14165/j.cnki.hunansci.2012.22.001

第31卷第22期Vol.31No.22

企业技术开发DEVELOPMENT TECHNOLOGICAL OF ENTERPRISE

2012年8月

2012年8月Aug.2012

基于Matlab R2011b 的卫星轨道计算

牛晓楠，蒋丹丹，霍红娟

（中国矿业大学环境与测绘学院，江苏徐州221116）

摘

要：文章阐述了计算卫星在地心坐标系中的位置，利用开普勒方程，根据广播星历的21个参数计算出开普勒轨道参

数，再结合卫星运动学方程计算出卫星的坐标。在此基础上，利用MALAB R2011b 软件读取卫星导航（O ，N ）文件，在MALAB

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

（School of Environment Science and Spatial Informatics ，China University of M ining&Technology，Xuzhou ，Jiangsu 221116，China ）

Abs tr act: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 21parameters of broadcast ephemeris ，then calculate coordinate of satellite combining satellites kinematics equations.On these basis ，reading satellite navigation （O ，N ）files using M ATLAB R2011b software ，and then de-sign the program of calculating the position of satellite and some instructions of important function and important code label in the envi-ronment of matlab R2011b. Keywords ：kepler orbit parameters ；satellites kinematics equations ；satellite orbit calculation ；programming language

GPS 卫星导航与定位技术在测量工作中得到越来越广泛的应用。在利用GPS 信号进行导航定位时，为计算用户在地心坐标系中的位置，GPS 接收机需要测定测站到卫星的距离并且要知道同一卫星在同一时刻的地心坐标。卫星的地心坐标是从卫星的导航电文中提供的开普勒轨道参数和轨道摄动修正量按一定的公式计算的。本文介绍如何利用matlab 软件实现卫星坐标的计算。

设以O 为原点的直角坐标系为O-XYZ ，卫星的坐标为（X ，Y ，Z ），则卫星的地心向量X ，Y ，Z ），加速度X ，Y ，Z ），代入上式得：

r=g（a ，e ，I ，ω，τ，t ）dr =g（′a ，e ，I ，ω，τ，t ）

→→→→→→→→→

从上式可以看出，在二体问题情况下，给定6个轨道参数，即可确定任意时刻t 的卫星位置及其运动速度。

1卫星运动力学方程

研究卫星S 绕地球O 的运动，主要是研究卫星运动状

2

2.1

卫星轨道计算

卫星轨道方程

以卫星轨道平面为基准面的三维坐标系统分别为：

态随时间的变化规律。根据牛顿定律，可以方便地得到二体问题的卫星运动方程。

用F S 、F e 表示卫星与地球所受到的引力作用力，则

o-ξ、η、ζ。

以下是以开普勒第一定律为基础推导卫星相对地球的位置表达式结果：

①由开普勒第一定律，卫星的运行轨迹为椭圆，在p-p 2+q 2=1；o-q ②卫星S 的位置在ξ-c-η坐标系中（ξ，η）为ξ＝rcosf=acos E-ae=a（cosE-e ）；

③可以得到：η=q=bsin E=a姨sinE ；

④则在ξ-η-ζ坐标系统下，卫星相对地球质心的距离向量r=η=a 姨sinE ；ζ0

姨姨姨姨姨姨姨姨姨

F S =-GmM/r2

，式中，G 为万有引力常数。

F e =GmM/r2

则a=as -a e =-G（M+m）/r2

上式即为卫星S 相对于地球质心O 的运动方程。即卫星相对于地球质心的加速度如下：

→

e -r s =-G（M+mr 。

地球质量远远大于卫星质量，通常略去卫星质量m 项

→

r=re -r s =-GMr 。

收稿日期：2012-05-13作者简介：牛晓楠（1990—），女，河北石家庄人，大学本科，研究方向：

测绘工程。

ξ

姨姨姨姨姨姨姨姨姨

姨姨姨姨姨姨姨姨姨姨

a(cosE-e)

姨姨姨姨姨姨姨姨姨姨

第31卷第22期牛晓楠，等：基于M atlab R2011b 的卫星轨道计算

21

⑤卫星与地球质心的几何距离为：‖r ‖=a（1-ecosE ）。2.2

开普勒方程

开普勒第三定律可表示为：n=真近点角f 可表示为：

f=arctanη=arctansinE 。

②[week，sow]=gps_time（jd ）jd —儒略日；返回值—week ，GPS 周；sow ，GPS 秒。3.2.2

读星历文件模块

①rinexe （ephemerisfile ，outputfile ）。ephemerisfile —n 文件名（n 文件放在rinexe 函数所在目录下）；outputfile —文本文件名，将读到n 文件的数据存放在该文本文件中。

②eph =get_eph（ephemeridesfile ）。ephemeridesfile -文件名，将存放在该文本中的导航电文数据提取出来，为计算卫星位置做准备；返回值—eph 星历矩阵，将读到的数据存入内存。3.2.3

计算卫星位置模块

①读观测值文件（*.o文件）函数。在求解卫星位置时，第一需要利用o 文件中每个历元的历元时刻t 。在计算某时刻卫星位置时，这里的某时刻便是o 文件历元时刻t 。第二根据PRN 号需要利用读取的每个历元不同的卫星PRN 号。

和历元时刻t 在广播星历n 文件中搜索相同的卫星PRN 号、合适的历元时刻，利用其对应的数据，计算卫星位置。

②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 中第2列用于计算卫星位置。第1，

③X=satpos（tx_GPS，Eph （：，k ））。tx_GPS—上节所述的Eph （：，k ）—Eph 星历矩阵中归化时间，用儒略日表示的；的某一列数据；返回值—卫星在地心地固坐标系中坐标。

姨

。

2.3卫星坐标计算

GPS 卫星在椭圆轨道上，以3.87km/s的速度每天旋转

k 卫星在t j 时刻的地心坐标可表示为：两圈，

姨姨姨姨姨姨姨姨姨

K 姨X （t ）j 姨

K

姨

姨姨姨姨姨姨

k k k

=R（（（Y （t ）r k j sinf k j j R j R j 3-Ω）3-i ）3-ω）j

K

Z （t ）j

姨

姨姨姨姨姨姨姨姨

r k j cosf k j 姨姨0

姨

姨姨姨姨姨姨

其中，Ωk j ，i k j ，ωk j 分别为k 卫星t j 时刻的升交点赤径，卫

k k k

星轨道面倾角，近地点角距。R （R （R （j ，j ，j 为3-Ω）3-i ）3-ω）

旋转矩阵。

3

3.1

卫星坐标计算算法步骤

算法

①计算卫星运行的平均角速度：nn 0=

姨

=，

a 3

姨3

n=n0+△n 。

②计算归化时间：t=t′-△t ，t k =t-toe 。

③观测时刻卫星平近点角M K 的计算：M K =MO +ntk 。④计算偏近点角E K ，E K =MK +esinE（M K 以弧度计算）。K E K ，

⑤真近点角f k 的计算：f k =arctanη=arctansinE k 。

k ⑥升交距角Φk 的计算：Φk =Vk +ω。⑦摄动改正项δu ，δr ，δi 的计算。

⑧计算经过摄动改正的升交距角u k 、卫星矢径r k 和轨道倾角i k 。

⑨计算卫星在轨道平面坐标系的坐标

4程序的实现过程

程序的开发基于MATLAB R2011b 平台实现的。

姨

x k =rk cosu k

。

y k =rk sinu k

MATLAB 是矩阵实验室（Matrix Laboratory ）的简称，是美国Math Works 公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语MATLAB R2011b 界面，如图1所示。言和交互式环境。

⑩计算观测时刻升交点经度Ωk ：Ωk =ΩO +（Ω-ωe ）t k -ωe t oe 。

訛计算卫星在地心固定坐标系中的直角坐标：輯輥

姨

姨姨姨姨姨姨姨姨

姨

X k 姨x k cos Ωk -y k cosi k sin Ωk 姨姨姨姨姨姨姨=姨。Y k 姨x k sin Ωk +yk cosi k cos Ωk 姨

姨姨姨

姨

姨

姨

Z k

姨姨姨姨姨姨

y k sini k

姨姨姨姨姨姨姨姨姨

姨姨姨

訛卫星在协议地球坐标系中的位置：輰輥

姨姨姨姨姨姨姨姨姨

X 姨姨Y Z

姨

姨姨姨姨姨CTS 姨

X k 姨姨Z k

姨

姨姨姨姨姨姨

=RY (-xP )R （）=Y k =Y -y P

姨

姨姨姨姨姨姨姨姨

10

01X P

X P 1

姨-Y P Y k 姨姨

-X P

姨姨姨姨姨姨姨姨姨姨姨姨姨姨姨姨姨姨

X k 姨姨

姨

Z k

姨姨

姨

3.23.2.1

重要MATLAB 函数说明时间转换模块

图1

MATLAB R2011b 界面图

主要用到两个函数，分别如下：

①jd =julday（year ，month ，day ，hour ）year —年，四位，如2011；month —月，1-12，day —日；hour —时间，小时+不足一小时部分（用小数表示）返回值-儒略日。

5结语

文章利用开普勒定律，卫星运动学（下转第36页）

36

术水平也更进一步。

编码器数值A A

31～400.697769.2172.7432

41～500.608642144.1360

编码器数值A A

取料结束

图2流程图A

2012年8月

M=Yθ/360（编码器记数值M ，编码器工作线数Y ）取料量Q 是恒定的，物料密度ρ和行走速度υ是定值，取料截面S 保持恒定就能满足上述要求。如图1取料截面S 1=h1l 1/2，S 2=h2l 2/2，S 3=h3l 3/2…S n =hn l n /2。S 1=S2=S3…S n 所示，

利用绘图法确定l 1、l 2、l 3在图面2的位置，确定交接点，利用CAD 绘图工具标注出交接点到定滑轮的距离Hn ，计算出每一层钢丝绳的长度差δl ，根据交点到卷扬滚筒钢丝绳的倍率N ，算出钢丝绳的释放长度L ，用公式n=L/πD 、θ=n/360得出每层悬臂下降角度θ，已知编码器线数Y ，就可以算出编码器需采集的数据M=Yθ/360。

以取料量80t/h，物料密度ρ=0.8～1.0t/m3，行走速度υ=2m/min，电动葫芦直径D=510.5mm ，编码器线数Y=900，得出表1数据。

表1

共50层

每层悬臂下降角度θ钢丝绳释放长度Lmm 电动葫芦旋转角度n 编码器记数值M

取料机每层下降数值表0～100.9571179264.6662

11～200.8781051.2236590

21～300.[1**********]0

下降M

下降M

下降M

下降M

下降M

取料机分50层取料，每10层为一组数据。初始状态下编码器复位为零，第1组下降记数到∑1=662×10=6620，第2组下降记数到∑2=∑1+590×10=12520，第3组下降记数到∑3=∑2+510×10=17620，第4组下降记数到∑4=∑3+432×10=21940，第5组下降记数到∑5=∑4+360×10=25540，计算下降总数是否在所使用的高速记数模块记数范围之内，之后在PLC 程序中，用当前记数值A 比较每组数据∑N ，根据比较结果下降相应的编码器记数值M 。其流程图，如图2所示。

侧式刮板取料机就是通过以上计算编程，由编码器准确地控制着物料取料量的恒定。由于增量型编码器的应用，使得取料工艺更加完善，侧式刮板取料机整体技（上接第21页）方程和卫星轨道进行计算。在MATLAB R2011b 良好的工作平台和编程环境（提供了比较完备的调试系统，程序不必经过编译就可以直接运行，而且能够及时地报告出现的错误及进行出错原因分析）上进行编程开发。其程序语言简单易用，语法特征与C++语言极为相似，而且更加简单，更加符合科技人员对数学表达式的书写格式。

用户可以在命令窗口中将输入语句与执行命令同步，也可以先编写好一个较大的复杂的应用程序（M 文件）可拓展性强，具有强大后再一起运行，语言可移植性好、的科学计算机数据处理能力，并且运算速度快、操作简单。

参考文献：

[1]晁代坤，孟红莉，等. 旋转编码器在过程控制系统中的应

用[J].钢管，1999，（6）.

[2]姜守成，公立滨. 旋转编码器在PLC 控制电梯中的应用[J].

应用能源技术，2005，（6）.

[3]于晓东. 旋转编码器在回转体回转角度测量中的应用[J].

林业机械与木工设备，2005，（2）.

[4]李震，刑忠南. 旋转编码器在装箱机中的应用[J].包装与食

品机械，2005，（4）.

参考文献：

[1]李征航，黄劲松.GPS 测量与数据处理[M].武汉：武汉大学

出版社，2005.

[2]刘基余，李征航，王跃虎，等. 全球定位系统原理及其应用

[M].北京：测绘出版社，1993.

[3]徐绍铨，张华海，杨志强，等.GPS 测量原理及应用[M].武

汉：武汉大学出版社，2011.

（第2版）[M].北京：高等[4]刘卫国.MATLAB 程序设计与应用

教育出版社，2006.

（第3版）[M].北京：电子工业出版社，[5]陈杰.MATLAB 宝典

2011.