卫星轨道问题.docx
《卫星轨道问题.docx》由会员分享，可在线阅读，更多相关《卫星轨道问题.docx（17页珍藏版）》请在冰豆网上搜索。
卫星轨道问题
南京师范大学2015数学建模校内赛问题
卫星轨道问题
自1957年世界第一颗人造卫星发射升空以来，人造卫星在军事、气象预报、车辆导航、大地测量等方面有越来越多的应用。
大部分卫星应用的前提是精确地确定卫星的轨道，根据牛顿定律，人造卫星绕地球沿椭圆轨道周而复始运行，但由于地球的非球形摄动力、日月引力、太阳辐射压、大气阻力等影响，卫星的实际轨道不会重复，确定卫星的精确位置并不容易。
为描述卫星位置，需要建立适当的坐标系。
本题惯性坐标系定义如下，取地球中心O为原点，地球自转轴取为
轴，指向北极为正向，
轴由O指向春分点，再按右手系确定
轴，建立直角坐标系。
地心O在绕太阳椭圆轨道上运动，所以理论上该坐标系是非惯性系。
但地球公转周期远大于卫星的观测弧段时长，故本题在短时间内认定该系为惯性坐标系，该坐标系不随地球旋转。
由于受到岁差章动极移等影响，实际的春分点和地球自转轴都发生变化，国际大地测量学协会（IAG）和国际天文学联合会（IAU）决定，以2000年1月1日12时UTC为标准历元，定义协议天球坐标系，其z轴和x轴分别指向2000年1月1.日12时的瞬时平北天极和瞬时平春分点，这个坐标系称为J2000平赤道平春分点坐标系，将此坐标系记为ICRS，求解卫星运动方程时总是在该坐标系中进行。
卫星的观测和应用是在地球上进行的，因此还需定义和地球固连的坐标系，和天球坐标系类似，国际上有一个协议地球坐标系，最新版本为ITRF2008，而美国国防制图局定义了一个地心地固坐标系WGS84坐标系，由于WGS84坐标系和ITRF相差很小，而GPS卫星使用WGS84坐标系，因此本题使用WGS84坐标系作为地球坐标系，WGS84直角坐标系和ICRS的转换公式见附录。
以卫星观测站P为原点，P点的法线为z轴（指向天顶为正），以子午线方向为x轴（向北为正），y轴与x、z轴垂直（向东为正），可以建立站心地平直角坐标系
以观测测站P为原点，至卫星s的距离r、卫星的方位角A、卫星的高度角h可以建立站心地平球坐标系。
这些坐标系的转换关系见附录中的公式。
请研究如下问题
1.问题一
已知卫星在时刻T0惯性坐标系ICRS中位置向量为（11000.0,39000.0,-4900.0）km,
速度为（-1.5,1.0，-0.1）km/s
计算卫星在该时刻的开普勒轨道根数，以及从T0时刻开始后一小时每次10分钟间隔的惯性系ICRS位置向量。
2.问题二
设一个观测站坐标北纬35.4251度，西经116.8893度，大地高为978.98,（单位km）
测到某卫星的距离和角度
Datautc方位角A（度）俯仰角h（度）距离（km）
2015/03/0208:
00:
00-165.837940.4421969.45
2015/03/0209：
00:
00-158.316170.665120514.34
计算卫星在两次观测时刻的卫星在坐标系ICRS中的位置向量，用这两组位置向量计算卫星轨道根数（第一个观测时刻），分析岁差、章动、极移对计算卫星轨道的影响。
3.问题3
设有卫星轨道观测数据，数据格式如下
日期
时
分
秒
wgs84x
wgs84y
wgs84z
2015/3/2
0
5
13
4273780.417
-3405983.549
-4581931.3
0
10
13
3649554.012
-1759990.152
-5867671.8
第一列为日期，第二、三、四列分别为小时、分、秒（UTC时间），第五、六、七列分别为WGS84坐标系x,y,z坐标，数据文件为data1.xlsx。
请估计卫星的实际轨道，给出卫星在1点到12点处整点处的坐标（ICRS坐标系和WGS84坐标系）
4.问题4
卫星数据同第三问。
设有北纬32.05度，东经118.84度，大地高为200米的卫星观测站，确定卫星可由该观测站最长可观测时段（卫星高于当地地平线），如果仅仅只有时间连续一段可观测数据，
如何外推卫星轨道？
附录
卫星轨道根数
轨道倾角i：
inclinationangle赤道平面与卫星轨道平面间的夹角
升交点卫星由南往北飞行轨迹在赤道平面上交点。
升交点赤经Ω：
therightascensionofascendingnode：
春分点轴向东度量到升交点的角度
轨道半长轴a：
semi-majoraxis
偏心率e：
eccentricity（0≤e<1）焦距与轨道半长轴之比。
近地点幅角ω:
argumentoftheperigee在轨道平面内升交点和近地点与地心连线间的夹角。
平近点角M:
meananomaly
若卫星通过近地点的时刻为tp，卫星的平均角速度为n，则时刻t的平近点角M=n（t-tp）。
坐标系
天球的概念
以地球质心M为球心，以任意长为半径的假想球体称为天球。
天文学中常将天体沿天球半径方向投影到天球面上，再根据天球面上的参考点、线、面来确定天体位置。
（1）天轴与天极
地球自转轴的延伸直线为天轴，天轴与天球面的交点称为天极，交点Pn为北天极，位于北极星附近，Ps为南天极。
（2）天球赤道面与天球赤道
通过地球质心M且垂直于天轴的平面称为天球赤道面，与地球赤道面重合。
天球赤道面与天球面的交线称为天球赤道。
（3）天球子午面与天球子午圈
包含天轴的平面称为天球子午面，与地球子午面重合。
天球子午面与天球面的交线为一大圆，称为天球子午圈。
天球子午圈被天轴截成的两个半圆称为时圈。
（4）黄道
地球绕太阳公转的轨道面称为黄道面。
黄道面与赤道面的夹角ε称为黄赤交角，约为23.5°。
黄道面与天球面相交成的大圆叫黄道，也就是地球上的观测者见到的太阳在天球面上的运行轨道。
由于地球自转，对于地面上的观测者来说，天球赤道面不动而黄道面每日绕天轴旋转一周。
又由于地球绕太阳公转，直观上看，太阳在黄道上每日自西向东运行约1°，每年运行一周。
（5）黄极
通过天球中心且垂直于黄道面的直线与天球面的两个交点称为黄极，靠近北天极Pn的交点Πn称为北黄极，Πs称为南黄极。
（6）春分点
当太阳在黄道上从天球南半球向北半球运行时，黄道与天球赤道的交点称为春分点，也就是春分时刻太阳在天球上的位置，如图中的Γ。
春分之前，春分点位于太阳以东。
春分过后，春分点位于太阳以西。
春分点与太阳之间的距离每日改变约1°。
天球坐标系
常用的天球坐标系有天球空间直角坐标系和天球球面坐标系。
天球空间直角坐标系的坐标原点位于地球质心。
z轴指向北天极Pn，x轴指向春分点Γ，y轴垂直于xMz平面，与x轴和z轴构成右手坐标系。
在天球空间直角坐标系中，任一天体的位置可用天体的三维坐标（x，y，z）表示。
天球球面坐标系的坐标原点也位于地球质心。
天体所在天球子午面与春分点所在天球子午面之间的夹角称为天体的赤经，用α表示；天体到原点M的连线与天球赤道面之间的夹角称为赤纬，用δ表示；天体至原点的距离称为向径，用r表示。
这样，天体的位置也可用三维坐标（α，δ，r）唯一地确定。
天体的天球空间直角坐标系和球面坐标系是同一天体位置的不同表达方式。
两种表达方式可通过下面的式子进行转换。
岁差与章动的影响
地球绕自转轴旋转，在无外力矩作用时，其旋转轴指向应该不变。
但由于日月对地球赤道隆起部分的引力作用，使得地球自转受到外力矩作用而发生旋转轴的进动现象，即从北天极上方观察时，北天极绕北黄极在圆形轨道上沿顺时针方向缓慢运动，致使春分点每年西移50.2″，25800年移动一周。
这种现象叫岁差。
在岁差影响下的北天极称为瞬时平北天极，相应的春分点称为瞬时平春分点。
瞬时平北天极绕北黄极旋转的圆称为岁差圆。
协议天球坐标系
由上可知，北天极和春分点是运动的，这样，在建立天球坐标系时，z轴和x轴的指向也会随之而运动，给天体位置的描述带来不便。
为此，人们通常选择某一时刻作为标准历元，并将标准历元的瞬时北天极和真春分点作章动改正，得z轴和x轴的指向，这样建立的坐标系称为协议天球坐标系。
国际大地测量学协会（IAG）和国际天文学联合会（IAU）决定，从1984年1月1日起，以2000年1月1.5日为标准历元。
也就是说，目前使用的协议天球坐标系，其z轴和x轴分别指向2000年1月1.5日的瞬时平北天极和瞬时平春分点。
将此坐标系记为ICRS.J2000.0
为了便于区别，z轴和x轴分别指向某观测历元的瞬时平北天极和瞬时平春分点的天球坐标系称为平天球坐标系，z轴和x轴分别指向某观测历元的瞬时北天极和真春分点的天球坐标系称为瞬时天球坐标系。
为了将协议天球坐标系的坐标转换为瞬时天球坐标系的坐标，须经过如下两个步骤的坐标转换。
（1）将协议天球坐标系的坐标转换为瞬时平天球坐标系的坐标
以
和
分别表示天体在瞬时平天球坐标系和协议天球坐标系中的坐标，因两坐标系原点同为地球质心，所以只要将协议天球坐标系的坐标轴旋转三次，便可转换为瞬时平天球坐标系的坐标，转换公式如下：
式中ζ、θ、z为坐标系绕z轴和x轴旋转的角度，其值由观测历元与标准历元之间的时间差计算。
为岁差矩阵。
（2）将瞬时平天球坐标系的坐标转换为瞬时天球坐标系的坐标
以
表示瞬时天球坐标系的坐标，则转换公式如下：
式中ε为观测历元的平黄赤交角；Δψ，Δε分别为黄经章动和交角章动。
，称为章动矩阵
协议地球坐标系
地球坐标系
确定卫星位置用天球坐标系比较方便，而确定地面点位则用地球坐标系比较方便。
最常用的地球坐标系有两种，一种是地球空间直角坐标系，另一种是大地坐标系。
地球空间直角坐标系的坐标原点位于地球质心（地心坐标系）或参考椭球中心（参心坐标系），z轴指向地球北极，x轴指向起始子午面与地球赤道的交点，y轴垂直于xoz面并构成右手坐标系。
大地坐标系是用大地经度L、大地纬度B和大地高H表示地面点位的。
过地面点P的子午面与起始子午面间的夹角叫P点的大地经度。
由起始子午面起算，向东为正，叫东经（0°~180°），向西为负，叫西经（0°~-180°）。
过P点的椭球法线与赤道面的夹角叫P点的大地纬度。
由赤道面起算，向北为正，叫北纬（0°~90°），向南为负，叫南纬（0°~-90°）。
从地面点P沿椭球法线到椭球面的距离叫大地高。
同一地面点在地球空间直角坐标系中的坐标和在大地坐标系中的坐标可用如下两组公式转换。
式中e——子午椭圆第一偏心率，可由长短半径按式
算得，
算得。
公式中的B必须用迭代的方法求解。
Matlab中有函数作变换，将经纬度转为直角坐标
[x,y,z]=geodetic2ecef（phi,lambda,h,ellipsoid）
其中phi为纬度，lambda为经度，h为大地高，ellipsoid=[a,e]，a为长半轴，e为偏心率
将直角坐标转为经纬度使用函数
[phi,lambda,h]=ecef2geodetic（x,y,z,ellipsoid）
注意matlab使用的是弧度
地极移动与协议地球坐标系
由于地球不是刚体，在地幔对流以及其它物质迁移的影响下，地球自转轴相对于地球体发生移动，这种现象叫地极移动，简称极移。
在建立地球坐标系时，如果使z轴指向某一观测时刻的地球北极，这样的地球坐标系称为瞬时地球坐标系。
显然，瞬时地球坐标系并未与地球固连，因而，地面点在瞬时地球坐标系中的位置也是变化的。
为了比较简明地描述地极移动规律，国际纬度局根据1900.0至1905.0期间5个国际纬度站的观测结果取平均，定义了协议原点（CIO）。
过CIO作地球切平面，并以CIO为原点建立平面直角坐标系，其中xp轴指向格林尼治方向，yp轴指向西经90°方向。
某一观测时刻的地极位置可用瞬时地极坐标xp和yp表示。
国际地球自转服务组织（IERS）定期公布瞬时地极坐标和各年度的平均地极坐标，
z轴指向协议地球北极的地球坐标系称为协议地球坐标系。
瞬时地球坐标系与协议地球坐标系之间的坐标可通过下式转换
，
称为极移矩阵,该矩阵中的xp和yp为弧度。
在2015年3月2日0h（UTC）,极移数据为x=3.458mas,y=358.634mas
2015年3月3日0h（UTC）,极移数据为x=3.877mas,y=360.608mas
为确定地面点的位置，需要将卫星在协议天球坐标系中的坐标转换为协议地球坐标系中的坐标，转换步骤为：
协议天球坐标系——瞬时平天球坐标系——瞬时天球坐标系——瞬时地球坐标系——协议地球坐标系。
瞬时天球坐标系与瞬时地球坐标系的坐标原点相同，z轴指向相同，只是两坐标系的x轴在赤道上有一夹角，角值为春分点的格林尼治恒星时。
因此只需将瞬时天球坐标系绕z轴旋转春分点的格林尼治恒星时时角GAST即可。
计算公式如下：
，
称为地球自转矩阵
WGS-84坐标系
WGS-84坐标系的原点在地球质心，Z轴指向BIH1984.0定义的协定地球极（CTP）方向，X轴指向BIH1984.0的零度子午面和CTP赤道的交点，Y轴和Z、X轴构成右手坐标系。
它是一个地心地固坐标系。
基本参数
长半径：
a=6378137±2（m）；
地球引力常数：
GM=3986005×D8±0.6×D8；
正常化二阶带谐系数：
C20=-484.16685×D-6±1.3×D-9；
J2=108263×D-8
地球自转角速度：
ω=7292115×D-11rad/s±0.150×D-11rad/s
WGS-84坐标系和国际协议地球系只有厘米级的差异，因此从ICRS到WGS-84的变换可以表示为
注意有的文献上的岁差矩阵的定义和这里的定义互为逆矩阵或者说是相差一个转置。
站心赤道直角坐标系与站心地平直角坐标系
P是测站（经度，纬度，大地高分别为L,B,H），O是球心。
以P为原点建立与球心空间直角坐标系相应坐标轴平行的坐标系叫做站心赤道直角坐标系。
显然，站心赤道直角坐标系与地心直角坐标系坐标系间有简单的平移关系。
站心地平直角坐标系
以P1为原点，P1点的法线为z轴（指向天顶为正），以子午线方向为y轴（向北为正），x轴与y、z轴垂直（向东为正）。
该坐标系又称为ENU（东北上）坐标系。
MATLAB将地心直角坐标系转为地平直角坐标系的方法为
[xl,yl,zl]=ecef2lv（x,y,z,phi0,lambda0,h0,ellipsoid）
其中x,y,z为地心坐标系坐标，phi0,lamda0,h0分别为测站的纬度，经度，大地高（纬度、经度用弧度表示）
X1,y1,z1为地平坐标系坐标。
MATLAB将地平直角坐标系转为地心直角坐标系的方法为
[x,y,z]=lv2ecef（xl,yl,zl,phi0,lambda0,h0,ellipsoid）
卫星方位角，高度角（俯仰角）
方位角是测站—卫星向量在水平面的投影与向北方向的夹角，从北至东为正，高度角是测站—卫星向量与水平面的夹角
由于文献中对地平直角坐标系的x,y,z的解释不同，因此对方位角和高度角的计算公式不同。
将卫星在东向轴的坐标记为SE,北向轴的坐标记为SN,天向轴的坐标记为SU，
则方位角A和高度角E的计算公式为
，其中A的计算要考虑SE和SN的符号，按照这个公式计算总是正确的。
时间系统
1、恒星时（SiderealTime—ST）
以春分点为参考点，由春分点的周日视运动所确定的时间，称为恒星时。
春分点连续两次经过本地子午圈的时间间隔为一个恒星日，含24个恒星小时。
所以恒星时在数值上等于春分点相对于本地子午圈的时角。
因为恒星时是以春分点通过本地子午圈时为原点计算的，同一瞬间对不同测站的恒星时各异，所以恒星时具有地方性，有时也称之为地方恒星时。
由于岁差、章动的影响，严格地讲，地球自转轴在空间的指向是变化的，春分点在天球上的位置并不固定，所以对于同一历元所相应的真北天极和平北天极，也有真春分点和平春分点之分。
因此，相应的恒星时也有真恒星时和平恒星时之分。
格林尼治真恒星时记为GAST,平恒星时记为GMST
2、平太阳时（MeanSolarTime—MT）
由于地球的公转轨道为一椭圆，根据天体运动的开普勒定律已知，太阳的视运动速度是不均均的。
如果以真太阳作为观察地球自转运动的参考点，那将不符合建立时间系统的基本要求。
为此，假设一个参考点的视运动速度，等于真太阳周年运动的平均速度，且其在天球赤道上作周年视运动。
这个假设的参考点，在天文学中称为平太阳。
平太阳连续两次经过本地子午圈的时间间隔，为一个平太阳日，而一个平太阳日包含有24个平太阳时。
与恒星时一样，平太阳时也具有地方性，故常称为地方平太阳时或地方平时。
3、世界时（UniversalTime—UT）
以平子夜为零时起算的格林尼治平太阳时称为世界时。
世界时与平太阳时的尺度基准相同，其差别（12小时）仅在于起算点不同。
世界时系统是以地球的自转为基础的。
但是地球自转轴在地球内部的位置并不是固定的，即有极移现象并且地球的自转速度也不均匀，它不仅包合有长期的减缓趋势，而且还具有一些短周期的变化和季节性的变化，情况甚为复杂。
地球自转的不稳定性，就破坏了建立时间系统的基本条件。
为了弥补这一缺陷，从1956年开始,在世界时UT0中引入了极移改正,得到世界时UT1（含有地球自转速度变化的影响）,在UT1中加入地球自转速度的季节性改正得到世界时UT2。
UT2虽经地球自转季节性变化的改正，但仍含有地球自转速度长期变化和不规则变化的影响，所以世界时UT2仍不是一个严格均匀的时间系统。
4、原子时（AtomicTime——TA）
原子时秒长的定义为：
位于海平面上的铯原子基态两个超精细能级，在零磁场中跃迁辐射振荡9192631770周所持续的时间，为一原子时秒。
该原子时秒作为国际制秒（SI）的时间单位。
国际上大约有100座原子钟，通过相互比对，并经数据处理推算出统一的原子时系统，称为国际原子时TAI
5协调世界时（UTC）
1972年采用了一种以原子时秒长为基础，在时刻上尽量接近于世界时的一种折衷的时间系统，这种时间系统称为协调世界时（UTC）
协调世界时的秒长严格等于原子时的秒长，采用闰秒（或跳秒）的办法使协调时与世界时的时刻相接近。
当协调时与世界时的时刻超过±0.9s时，便在协调时中引入一闰秒（正或负），闰秒一般在12月31日或6月30日末加入。
具体日期由国际时间局安排并通告。
6地球时（TT）
地球时TT,以前称为TDT,和国际原子时的关系为
TT=TDT=TAI+32.184
在2015年3月2日0h（UTC），UT1-UTC=-528.5949ms
2015年3月3日0h（UTC）,UT1-UTC=-529.5678ms
目前TAI-UTC=35s
7GPS时间系统GPST
GPST属于原子时系统，它的秒长即为原子时秒长，GPST与国际原子时TAI相差19s。
有关系式：
TAI-GPST=19（s）
儒略日
儒略日（Julianday，JD）是指由公元前4713年1月1日，世界时中午12时开始所经过的天数,例如，1996年1月1日12:
00:
00的儒略日是2450084,
由于儒略日数字位数太多，国际天文学联合会于1973年采用简化儒略日（MJD），其定义为MJD=JD-2400000.5。
MJD相应的起点是1858年11月17日世界时0时。
例如1979年10月1日零时JD=2,444,147.5,MJD=44147
MATLAB计算JD的函数为
JD=juliandate（y,mo,d,h,mi,s）；
计算10月10日,2004年,12:
21:
00对应的儒略日
JD=juliandate（2004,10,10,12,21,0）
MATLAB计算MJD的函数为
MJD=mjuliandate（y,mo,d,h,mi,s）；
计算公式
GMST=24110.54841+8640184.812866*T_0+1.002737909350795*UT1+（0.093104-6.2e-6*T）*T*T;
其中，T0=（JD（0时ut1）-2451545））/36525T=（JD（ut1）-2451545））/36525
GAST=GMST+cos（ε）ΔΨ
本题中ΔΨ=4.75角秒
其中黄赤交角ε=23.43929111度-46.8150秒*T-0.00059秒*T*T+0.001813秒*T*T*T
岁差矩阵中三个角度计算公式
ζ=2.5976176+2306.0809506*T+0.3019015*T^2+0.0179663*T^3-0.0000327*T^4-0.0000002*T^5（秒）
θ=2004.1917476*T-0.4269553*T^2-0.0418251*T^3-0.0000601*T^4-0.0000001*T^5（秒）
z=2.5976176+2306.0803226*T+1.0947790*T^2+0.0182273*T^3+0.0000470*T^4-0.0000003*T^5（秒）
章动矩阵在本题中为
章动矩阵
升级会员 