卫星轨道常识

1、5.14常微分方程5.14.2 ODE解算指令的使用演示5.14.2.1解算指令简洁格式使用示例 【*例5.14.2.1-1】采用最简洁格式的ODE文件和解算指令，研究围绕地球 旋转的卫星轨道。图5.14.2.1-1地球轨道卫星运动加速度问题的形成轨道上运动的卫星，在Newton第二定律，和万有引力定律作用下，有 。即 ， ，而 。这里是引力常数，是地球的质量。又假定卫星以初速度在处进入轨道。构成一阶微分方程组 令 ，则(5.14.2 .1-5)初始条件为(5.14.2 .1-6)根据式(5.14.2.1-5)编写最简洁格式的ODE文件dYdt.mfunction Yd=DYdt(t,Y)%

2、t 一定是标量形式的自变量% Y必须是列向量global G ME %在函数中定义全局变量传递参数xy=Y(1:2);Vxy=Y(3:4); %前两个元素是位移量，后两个是速度量 Or=sqrt(sum(xy.”2);Yd=Vxy;-G*ME*xy/r3; %Yd 必须按式(5.14.2.1-5)编写，是与 Y 同维的列 向量。对微分方程进行解算global G ME %在主程序中定义全局变量传递参数G=6.672e-11;ME=5.97e24;vy0=4000;x0=-4.2e7;t0=0;tf=60*60*24*9;tspan=t0,tf; %指定解算微分方程的时间区间Y0=x0;0;0;

3、vy0; % 按式(5.14.2.1-6)给定初值向量 t,YY=ode45(DYDt,tspan,Y0); %X=YY(:,1); %输出Y的第一列是位移数据x(t)Y=YY(:,2); %输出Y的第二列是位移数据y(t) plot(X,Y,b,Linewidth,2); hold on axis(image) %保证x、y轴等长刻度，且坐标框恰包容图形 XE,YE,ZE = sphere(10); %产生单位球面数据RE=0.64e7; %地球半径XE=RE*XE;YE=RE*YE;ZE=0*ZE; %坐标纸上的地球平面数据 mesh(XE,YE,ZE),hold off % 绘地球示意图

4、图5.14.2.1-2卫星轨道【*例5.14.2.1-2】上例中，程序间的参数(如G和ME)传送，是依靠全 局变量形式实现的。一般说来，编写程序时，应尽量少用全局变量，以免引起混 乱。本例演示参数如何在指令间直接传送。要实现参数直接传送，必须对上例 中的程序进行修改，具体如下：重写ODE文件DYDt2.mfunction Yd二DYDt2(t,Y,flag,G,ME)% flag按ODE文件格式规定，必须是第三输入宗量。对它的赋值由ode45指 令自动产生。%第4、5宗量是被传递的参数switch flagcase %按规定：这里必须是空串。在此为真时，完成以下导数计算。 X=Y(1:2);V

5、=Y(3:4);r=sqrt(sum(X.2);Yd=V;-G*ME*X/r3;otherwiseerror( Unknown flag flag .);end按以下方法修改上例第(4)步中的程序，并运行之，可得相同结果。删去原主程序第1条指令；把原主程序的第8条指令改写为t,YY=ode45(DYDt2,tspan,Y0,G,ME); % 第 4 宗量取缺省设置%第5、6宗量是被传递参数5.14.2.2解算指令较复杂格式的使用示例【*例5.14.2.2-1】带事件设置的ODE文件及主程序编写演示。本例将以较 高精度计算卫星经过近地点和远地点的时间，并在图上标志。ODE文件的编写DYDt3.m

6、function varargout二DYDt3(t,Y,flag,G,ME,tspan,Y0)% DYDt3.m供主程序调用的ODE函数文件%本文件自带三个子函数：f,fi,fev。% t, Y分别是自变量和一阶函数向量，是最基本的输入宗量。% flag第三输入宗量，它专供解算指令(如ode45)作调用通知。%在运行中，解算指令会根据需要向flag发不同的字符串。% varargout是变输出宗量。它由变维的元胞数组构成。每个元胞中可以 存放指令%所产生的任意形式的数据。switch flagcase %必须用空串符。情况为真时，计算导数dY/dt = f(t,Y)。varargout1 =

7、 f(t,Y,G,ME); %输出为一个元胞，容纳f子函数的一个输出 Ydcaseinit %必须用init，情况为真时，传送计算区间、初值、设 直参数。varargout(1:3 = fi(tspan,Y0); %输出为三个元胞，容纳fi子函数的三 个输出量case events %必须用events，情况为真时，设置事件性质。varargout(1:3 = fev(t,Y,Y0); %输出三个元胞，容纳fev子函数的三个 输出量otherwiseerror( Unknown flag flag .);end% function Yd = f(t,Y,G,ME)%计算导数子函数，被父函数DY

8、Dt3调用。X=Y(1:2);V=Y(3:4);r=sqrt(sum(X.2);Yd=V; -G*ME*X/r3;% function ts,y0,options = fi(tspan,Y0)%设置时间区间、初值、算法参数子函数，被父函数DYDt3调用。ts=tspan;y0 = Y0;options = odeset( Events , on , Reltol ,1e-5, Abstol ,1e-4);%开动ode45的事件判断功能，设置相对误差和绝对误差。% function value,isterminal,direction = fev(t,Y,Y0)%事件判断子函数，被父函数DYDt

9、3调用。dDSQdt = 2 * (Y(1:2)-Y0(1:2) * Y(3:4);%dDSQdt是离初始点的二乘距离u的时间导数du/dt，而 u=(x-x0)”2+(y-y0)”2。value = dDSQdt; dDSQdt; %定义两个穿越0的事件direction = 1; -1; %第一事件：以渐增方式穿越0。第二事件：以渐减方 式穿越0isterminal = 1; 0; %第一事件发生后，终止计算；而第二事件发生后，继 续计算。运行以下主程序G=6.672e-11;ME=5.97e24;vy0=4000; x0=-4.2e7;t0=0;tf=60*60*24*9;tspan=t

10、0,tf;Y0=x0;0;0;vy0;t,YY,Te,Ye,Ie=ode45(DYDt3, ，口，口，G,ME,tspan,Y0); % X=YY(:,1);Y=YY(:,2);plot(X,Y,b,Linewidth,2);hold ontext(0,6e7,轨道，Color,b) %产生蓝色文字注释axis(image); %保证x、y轴等长刻度，且坐标框恰包容图形%在三个事件发生点上画标记plot(Ye(1,1),0.4e7+Ye(1,2),r,MarkerSize,10)plot(Ye(2,1),0.4e7+Ye(2,2),bv,MarkerSize,10)plot(Ye(3,1),-

11、0.4e7+Ye(3,2),b,MarkerSize,10)%把轨道的半周期和全周期标在图上text(0.8*Ye(3,1),-2e7+Ye(3,2),t3=sprintf(%6.0f,Te(3)text(0.8*Ye(2,1),1.5e7+Ye(2,2),t2=sprintf(%6.0f,Te(2)%在x-y坐标上画地球XE,YE,ZE = sphere(10);RE=0.64e7;XE=RE*XE;YE=RE*YE;ZE=0*ZE;mesh(XE,YE,ZE)text(1e7,1e7,地球，Color,r), hold off % 产生红色文字注释图5.14.2.2-1带事件标注的卫星轨道

12、图5.14.3关于ODE文件的说明 （1）ODE文件的模板下面就是MATLAB的模板文件odefile.m，为便于读者阅读，本书作者适当地 给以注解说明。odefile.mfunction varargout = odefile(t,y,flag,p1,p2)% odefile.m是MATLAB提供的模板文件。用户可根据需要，使用任何其他文件 名。% varargout是MATLAB专门设计的变长度输出宗量。它由元胞数组构 成，因此%可适应任意多个、任意形式的输出宗量。请参看第8.5.2节。% t自变量。不管在下面指令中是否出现，自变量t必须作为第一输入宗量。% y 一阶微分方程组的列向量形式

13、函数。它必须作为第二输入宗量。% flag切换变量。它必须处在第三输入宗量位置。用户无须也不要对它直接赋%值；它的赋值由微分方程解算指令(solver)，自动产生。% pl, p2是被传递的参数。这里，作为示意，仅列出两个。%根据需要，传递参数的数目不受限制。%以下是switch-case构成的多向选择控制。应注意：% ( A)各情况的情况表达式是MATLAB规定的，不得任意改变。%( B)各情况所执行的任务也是规定的，不得任意改变。% ( C)各情况内，变长输出宗量元胞数是规定的，不得任意改变。% ( D)各情况内，(等式右边的)函数名称可以改变，但相应子函数名称要 一致。%( E)各情况内

14、，(等式右边的)函数中所包含的t , y是必须的，%不管它们是否以显式出现相应子函数体中。%( F)各情况内，(等式右边的)函数中所包含的pl, p2是传递参数的示 意性表示。%具体参数视需要而定，只要该参数已经传入odefile.m函数内存。switch flagcase %规定空字符串情况：专管一阶导数dy/dt = f(t,y)的计算 varargout1 = f(t,y,p1,p2);case init % 规定init情况：专管三个宗量tspan , y0 , options的 设置。varargout(1:3 = init(p1,p2);case jacobian %规定jacob

15、ian情况：专管计算解析的Jacobian矩阵 df/dy 。varargout1 = jacobian(t,y,p1,p2);case jpattern %规定jpattern情况：专管计算稀疏的数值Jacobian矩 阵 df/dy 。varargout1 = jpattern(t,y,p1,p2);case mass % 规定 mass 情况：专管计算质量矩阵(mass matrix)。varargout1 = mass(t,y,p1,p2);case events %规定events情况：专管事件定义和判断varargout(1:3 = events(t,y,p1,p2);otherw

16、iseerror( Unknown flag flag .);end %以下是odefile.m的子函数，它们分别与父函数中各情况对应。% function dydt = f(t,y,p1,p2)%空字符串情况内调用的子函数：专管一阶导数dy/dt的计算% dydt 一阶导数。该变量名可由用户按需要取名。%下面函数内容，由用户自己编写，最后产生输出dydt。dydt =由用户编写% function tspan,y0,options = init(p1,p2)% init情况内调用的子函数：专管三个宗量tspan , y0 , options的设 置。%三个输出宗量的名称可由用户自起，但各位置

17、上宗量的性质不能改变。%下面函数内容，由用户自己编写，最后产生输出三个输出。tspan =积分的时间区段或时间点集y0 =初值options =或用odeset设置算法参数，或用空阵符表示缺省设置% function dfdy = jacobian(t,y,p1,p2)% jacobian情况调用的子函数：专管计算解析的Jacobian矩阵df/dy。% 输出宗量的名称可由用户自起。%下面函数内容，由用户自己编写，最后产生输出输出dfdy。dfdy = Jacobian 矩阵% function S = jpattern(t,y,p1,p2)% jpattern情况调用的子函数：专管计算解析的

18、Jacobian矩阵df/dy。% 输出宗量的名称可由用户自起。% 下面函数内容，由用户自己编写，最后产生输出 S。S =稀疏形式的数值Jacobian矩阵% function M = mass(t,y,p1,p2)% mass情况调用的子函数：专管计算质量矩阵(mass matrix)。解刚性方 程时用。%输出宗量的名称可由用户自起。%下面函数内容，由用户自己编写，最后产生输出M。M =质量矩阵% function value,isterminal,direction = events(t,y,p1,p2)% ,events,情况内调用的子函数：专管三个宗量value,isterminal,

19、direction 的设置。%三个输出宗量的名称可由用户自起，但各位置上宗量的性质不能改变。%下面函数内容，由用户自己编写，最后产生输出三个输出。value =各元素是标量事件函数，即函数穿越0点为事件direction =各元素定义事件函数穿越方式：1为向正穿越0 ； -1为向负穿越0 ； 0不管方向isterminal =各元素定义事件发生后计算是否继续：1为终止计算；0为 继续计算%(2)ODE模板的使用方法5.14.4关于解算指令选项options的属性设置5.14.4.2 options属性处理和输出函数使用演示【*例5.14.4 .2-1】仍以卫星轨道问题为例(原题见例5.14.2

20、.1-1, 5.14.2.1-2 , 5.14.2.2-1)。本例演示如何通过对options域的直接设置， 借助微分方程解算输出指令，表现解算的中间结果。具体目标是：画出解向量中由构成的相平面。相平面的绘制是在微分方程解算中间逐步完成的。(1)编写ODE文件DYDt4.m (在DYDt3.m基础上删改而成) DYDt4.mfunction varargout二DYDt4(t,Y,flag,G,ME,tspan,Y0) switch flag case varargout1 = f(t,Y,G,ME);case initvarargout(1:3 = fi(tspan,Y0); otherwi

21、seerror( Unknown flag flag .);end%function Yd = f(t,Y,G,ME)X=Y(1:2);V=Y(3:4);r=sqrt(sum(X.2);Yd=V; -G*ME*X/r3;%function ts,y0,options = fi(tspan,Y0)ts=tspan;y0 = Y0;%采用向域直接赋值法，设置options属性。以供与odeset使用方法对照。 options.RelTol=1e-5;options.AbsTol=1e-4;options.OutputFcn=odephas2 ; %在积分进程中，绘制相平面图。options.Out

22、putSel=1 3; %解向量的第1、3分量分别为相平面图的横、 纵坐标量。编写脚本文件odeexp4.modeexp4.m%odeexp4.mG=6.672e-11;ME=5.97e24;vy0=4000;x0=-4.2e7;t0=0;tf=60*60*24*9;tspan=t0,tf;Y0=x0;0;0;vy0;%以下四条指令，为相平面图预置一个坐标轴范围clf,set(gca, xlim ,-5 25*1e7, ylim ,-3 3*1e3); % 设置适当坐标范围box on %框形坐标hold on; %使中间结果绘在同一幅图上ode45( DYDt4 ,G,ME,tspan,Y0

23、);hold off在MATLAB指令窗中运行以下指令，即可在图形窗中见到相平面图的逐步 绘制。shg,odeexp4图5.14.4.2-1微分方程输出函数相平面图的实时绘制卫星轨道常识一一轨道参数和简单计算【转贴】卫星轨道和TLE数据转自虚幻天空最近由于Sino-2和北斗的关系，很多网友贴了表示卫星运行轨道的TLE数据。 这里想对卫星轨道参数和TLE的格式做一个简单介绍。虽然实际上没有人直接读 TLE数据，而都是借助软件来获得卫星轨道和位置信息，但是希望这些介绍可以 对于理解卫星轨道的概念有所帮助。由于匆匆写成，可能有一些错误，如果看到 还请指出。前面关于轨道一部分写得较早，后来发现和杂志上

24、关于我国反卫的一篇文章里的 相应部分类似。估计都参考类似的资料，这个东西本身也是成熟的理论了。首先来看一下卫星轨道。太空中的卫星在地球引力等各种力的作用下做周期运 动，一阶近似就是一个开普勒椭圆轨道。由于其他力的存在（比如地球的形状， 大气阻力，其他星球的引力等等），实际的轨道和理想的开普勒轨道有偏离，这 个在航天里称为“轨道摄动”。这里我们暂时不看摄动，就先说说理想开普勒轨 道时的情况。为了唯一的确定一个卫星的运行轨道，我们需要6个参数，参见下面的示意图：轨道半长轴，是椭圆长轴的一半。对于圆，也就是半径轨道偏心率，也就是椭圆两焦点的距离和长轴比值。对于圆，它就是0.这两个要素决定了轨道的形状

25、轨道倾角，这个是轨道平面和地球赤道平面的夹角。对于位于赤道上空的同 步静止卫星来说，倾角就是0。升交点赤经：卫星从南半球运行到北半球时穿过赤道的那一点叫升交点。这 个点和春分点对于地心的张角称为升交点赤经。这两个量决定了卫星轨道平面在空间的位置。近地点幅角：这是近地点和升交点对地心的张角。前面虽然决定了轨道平面在空间的位置，但是轨道本身在轨道平面里还可以转 动。而这个值则确定了轨道在轨道平面里的位置。过近地点时刻，这个的意义很显然了。卫星位置随时间的变化需要一个初值。 有一点要指出的是，上面的6个参数并不是唯一的一组可以描述卫星轨道情况的 参数，完全也可以选取其他参数，比如轨道周期。但是由于完

26、备的描述也只需要 6个参数，所以他们之间存在着固定的换算关系。比如轨道周期就可以由半长 轴唯一来确定（这在下面讲TLE的时候也会涉及到），反之亦然。上面选取的这组 是比较自然的一组。下面讲讲TLE（Two-Line Element）两行数据。以北斗最近的数据为例 BEIDOU 2A 1 30323U 07003A 07067.68277059 .00069181 13771-5 44016-2 0 587 2 30323 025.0330 358.9828 7594216 197.8808 102.7839 01.92847527 650 真正的数据实际上是下面2行，但是上面有一行关于空间物体

27、其他情况的一些信 息（空间物体可以是卫星，可以是末级火箭，可以是碎片。这里简单起见，就叫 卫 星）。头一个是卫星名称。注意这个是会变的，而且不一定准确。卫星发射后 的头几个TLE数据里，往往只叫Object A, B, C.慢慢的会搞清楚哪个是卫 星，哪个是末级火箭，哪个是分离时的碎片，并且给予相应的名称。但是如果这猜测了，比如我们关于卫星的尺度，44016-2 0 587个是其他国家的保密卫星，则这个卫星名字就纯粹是美国的 的这个北斗。有些情况下，名称这一行里还包含了一些数字 亮度等等。TLE第一行数据1 30323U 07003A 07067.68277059 .00069181 1377

28、1-530323U30323是北美防空司令部（NORAD）给出的卫星编号。U代表不保密。我们看到的都是U，否则我们就不会看到这组TLE 了07003A国际编号，07表示2007年（2位数字表示年份在50年以后会出问题， 因为1957年人类发射了第一个轨道物体），003表示是这一年的第3次发射。A 则表示是这次发射里编号为A的物体，其他还有B，C，D等等。国际编号就是 2007-003A.07067.68277059这个表示这组轨道数据的时间点。07还是2007年，067表示 第67天，也就是3月8日。68277059表示这一天里的时刻，大约是16时22分左右。.000069181平均运动的对时

29、间一阶导数除2。注意这个并不是瞬时角速度13771-5平均运动对时间的二阶导数除6。44016-2 BSTAR阻力系数。这3个量都是用于轨道摄动模型里面的。其实上 前2个并没有真正被采用。0轨道模型。他们内部有不同数字代表不同模型，但是公布的都是0，也就是采用了 SGP4/SDP4轨道模型58表示这是关于这个空间物体的第58组TLE7最后一位是校验位TLE第二行数据2 30323 025.0330 358.9828 7594216 197.8808 102.7839 01.92847527 650 30323 NORAD卫星编号。025.0191轨道倾角。这个和前面讲的轨道倾角完全对应358.

30、9828升交点赤经，这个和前面讲的升交点赤经也完全对应7594216轨道偏心率，0.7597678，表示这是一个椭圆197.8808近地点幅角，这个和前面讲的也一样102.7839 平近点角。这个表示这组TLE对应的时刻时，卫星在轨道的什么位 置，具体细节有点复杂，就不赘述了。这个和前面讲的“过近地点时刻”可以互 相推导。01.92847527每天环绕地球的圈数。这个的倒数就是周期。可以看出北斗目前 的周期大约是12小时。而周期和轨道的半长轴有简单的换算关系。因此TLE 的关于轨道的6要素和我们前面说的6要素是完全可以互相推导的。65发射以来飞行的圈数0 校验位以上为 shh 原创。-darklighter轨道周期和半长轴的换算人造地球卫星运转周期T （秒）与圆轨道半径或椭圆轨道半长轴R （米）之间的 关系可用下列公式计算：R= （GMT2/4n2） （1/3）其中，GM=398.60047X1012，代入各常数后计算得知，R=21613.546XT2/3已知地球自转周期为86164.09053秒，卫星每天绕地球运转16圈，周期为地球自转周期的十六分之一（约1.5小时），