空间目标二体运动模型及其CPU实现

空间目标二体运动模型及其CPU实现

基于CUDA的空间目标轨道并行计算技术，共四节，其中第一、二、三节目录如下

1 基于CUDA的空间目标轨道计算需求与任务分析

2 基于CUDA的空间坐标系变换矩阵计算
2.1 基于SOFA的空间坐标系变换
2.2 采用CUDA Thrust实现的空间坐标系变换
2.3 采用CUDA Runtime实现的空间坐标系变换**

3 基于CUDA的标准二体模型轨道预推算法
3.1 空间目标二体运动模型及其CPU实现
3.2 采用CUDA Thrust实现的标准二体模型
3.1 采用CUDA Runtime实现的标准二体模型

---

一、基本原理

根据开普勒三定律可知，航天器的运动轨迹是以地心为焦点的椭圆。需指出这是指在惯性坐标系的情况。要想确定航天器在任一给定时刻的位置、速度，需要一些给定参数条件，这就是所谓轨道根数，根据轨道根数计算得到惯性坐标系下的位置、速度，再将其变换到地固坐标系。

（1）确定形状
航天器运行轨道为椭圆，椭圆形状采用半长轴a和偏心率e来表示，这2个参数可以确定轨道的形状。

（2）确定平面
如从数学角度确定三维空间中的平面，需4个参数的平面方程。对于处于惯性坐标系的航天器轨道，假定旋转中心和轴固定，只需描述其与中心和轴的关系，这通过轨道倾角i和升交点赤经Ω确定。

倾角是轨道面与赤道面之间的夹角，如图所示。倾角的范围在0°-180°之间：倾角为0°为赤道轨道；倾角在0°-90°之间，表明卫星运动方向与地球自转方向一致，称为顺行轨道；倾角为90°为极地轨道；倾角为90°-180°之间则表示卫星运动方向与地球自转方向相反，称为逆行轨道。

可以看出，过地心且给定倾角的平面有无数个，要把平面固定，还需一个参数，即升交点赤经。轨道面与赤道面相交于一直线，该直线穿过地心并向空间延伸，如图所示。该直线与卫星轨道有2个交点：一是在卫星由南半球进入北半球的上升段，称为升交点；二是在卫星由北半球进入南半球的下降段，称为降交点。地心到升交点方向与地心到平春分点方向的夹角即为升交点赤经，如图中Ω所示。

（3）确定平面上的轨道
通过轨道倾角和升交点赤经确定了空间平面，利用地心、半长轴和偏心率确定了平面上椭圆的形状，但符合以上条件的椭圆仍存在无数个，可以想象一下：以地心为中心，将一个符合条件的椭圆绕过地心的轨道平面法向旋转，可得到无数个椭圆。为此通过近地点幅角ω确定平面上的椭圆。

椭圆半长轴方向与椭圆有两个交点，分别是距离地心最近的点和最远的点，称为近地点和远地点。地心到近地点方向与地心到升交点方向的夹角即为近地点幅角ω，。对于每个确定的ω值，椭圆在平面上的位置固定。

（4）确定卫星位置
通过上面5个参数，已经唯一地确定空间椭圆，但并没有确定卫星的位置。可以想象：给定初始时刻，卫星可以从椭圆上任一位置出发，然后按动力学规律运行。

为进一步确定卫星位置，一个方法是指定卫星过近地点的时刻，称为过近地点时刻τ。但更为常用的方法是确定给定时刻卫星与近地点的角度差，即“近点角”。有2个可互相转化的量，一是真近点角，给出了轨道历元时刻地心到卫星位置方向与与地心到近地点方向的角度差）所示；还有一个量称为平近点角，后面再讨论。

综上所述，航天器轨道根数通常可采用（a,e,i,Ω,ω,τ）、（a,e,i,Ω,ω,f）、（a,e,i,Ω,ω,M）等形式表达，半长轴、偏心率、轨道倾角、升交点赤经、近地点幅角定义了空间的椭圆轨道，最后一参数分别以过近地点时刻或给定历元的真近点角、平近点角来表示时间的影响。

二、二体模型计算过程

下面介绍如何根据给定轨道根数计算卫星位置

所谓二体模型，就是只考虑地球和航天器之间引力作用的情况下，航天器位置、速度的计算方法。

如果计算出任一给定时刻卫星的真近点角，则根据椭圆方程可以得到卫星位置。由于卫星在轨道上并非匀速匀度，因此引入偏近点角和平近点角概念。

如图所示，以椭圆长半轴为半径做辅助圆。地心到卫星当前位置方向与地心到近地点方向所形成夹角f为真近点角；过卫星当前位置的垂线与辅助圆交于Q点，则OQ方向与地心到近地点方向的夹角即为偏近点角E。定义平近点角为

其中μ＝3.9860044181e^5，为地球引力常数；τ为卫星过近地点时刻。

平近点角的含义是假设航天器从近地点开始以平均角速度运行，经过(t-τ)时间后运行的角度。

如以（a,e,i,Ω,ω,M）形式给定的轨道根数，则通过下式计算任一时刻的偏近点角

式中，E为偏近点角，t为待计算时刻，e为偏心率，a为轨道半长轴，M为平近点角，t_0为轨道根数对应的时刻。该值需迭代计算

利用上式反复迭代计算，直到ε小于给定阈值，也可以迭代固定次数，如100次，但其效率往往低于前者。得到偏近点角后，真近点角采用下式计算

对于轨道根数以（a,e,i,Ω,ω,f）形式给定的情况，可先计算对应的平近点角，再按上述过程计算，平近点角计算为

建立卫星轨道坐标系：以地心为原点，地心到卫星当前位置为x轴，与x轴垂直的卫星运动方向为y轴，过原点的轨道面法线为z轴，则卫星轨道坐标系相当于位于惯性坐标系中的一个局部坐标系。轨道坐标系中的xy平面如图所示。

则在轨道坐标系中，卫星的位置矢量为

即矢量的x分量相当于地心到卫星的距离r，其他两个分量为0。其中，a为半长轴，e为偏心率，f为真近点角，x、y为卫星在轨道面平面坐标系中的平面坐标。

然后将轨道坐标系中的矢量变换到惯性坐标系中
①将轨道坐标系绕z轴顺时针旋转真近点角f，使得其x轴与地心到近地点方向重合

②将坐标系绕z轴按顺时针旋转近地点幅角ω，使得x轴与地心升交点方向重合

③绕x轴旋转轨道倾角i，使得轨道坐标系的xy平面与惯性坐标系的xy平面重合

④绕z轴旋转升交点赤经Ω，使得x轴指向平春分点方向，

最后再将惯性坐标系中矢量转换到地固系，这利用前面第2节讨论的空间坐标系变换矩阵实现。

三、CPU实现

二体模型的CPU串行代码如下

int		propagateOrbitWithTwoBody(timeOfSpace epoTime, OrbitElement orbitElement, timeOfSpace t0, timeOfSpace t1, int stepMilliSecond, std::vector<mat3x3>& mats,std::vector<vec3d>& ECIps, std::vector<vec3d>& ECFps){
	double epoJD0, epoJD1;
	iauCal2jd(epoTime._year, epoTime._month, epoTime._day, &epoJD0, &epoJD1);
	double tmp = double(60 * (60 * epoTime._hour + epoTime._minute) + epoTime._second+(double)epoTime._millisecond/1000) / 86400.0;
	epoJD1 += tmp;
	double n = sqrt(MIU_EARTH / orbitElement._semimajorAxis) / orbitElement._semimajorAxis;
	double timeOfCycle = PI * 2 / n;
	double tmp1 = sqrt(MIU_EARTH / orbitElement._semimajorAxis) / orbitElement._semimajorAxis;
	double tmp2 = sqrt((1 + orbitElement._eccentricity) / (1 - orbitElement._eccentricity));
	double rmat[3][3];
	iauIr(rmat);
	iauRz(-orbitElement._RAAN*PI / 180.0, rmat);
	iauRx(-orbitElement._inclination*PI / 180.0, rmat);
	iauRz(-orbitElement._argumentOfPerigee*PI / 180.0, rmat);
	double jd00, jd01 , jd10 , jd11 ;
	iauCal2jd(t0._year, t0._month, t0._day, &jd00, &jd01);
	tmp = (60 * (60 * t0._hour + t0._minute) + t0._second + (double)t0._millisecond / 1000) / 86400.0;
	jd01 += tmp;
	iauCal2jd(t1._year, t1._month, t1._day, &jd10, &jd11);
	tmp = (60 * (60 * t1._hour + t1._minute) + t1._second + (double)t1._millisecond / 1000) / 86400.0;
	jd11 += tmp;
	double dt = (jd00 + jd01) - (epoJD0 + epoJD1);
	int num = mats.size() ;
	for (int i = 0; i < num; i++)	{
		double t = dt + double(stepMilliSecond / 1000.0) / 86400.0*i;
		t = fmod(t, timeOfCycle);
		t *= 86400;
		double M = t * tmp1;
		M = fmod(M, PI * 2);
		double	radM = orbitElement._meanAnomaly*PI / 180 + M;	//加上初始值
		double E0 = radM;
		double E1;
		while (1){
			E1 = radM + orbitElement._eccentricity * sin(E0);
			if (fabs(E1 - E0) < 1e-10)break;
			E0 = E1;		}
		E1 = fmod(E1, PI * 2);
		double f = 2 * atan(tmp2 *tan(E1 / 2));
		if ((E1 - f) > PI)
			f = f + PI * 2;
		else if ((E1 - f) < -PI)
			f = f - PI * 2;
		double  tmpcosf = cos(f);
		double  tmpsinf = sin(f);
        r=orbitElement._semimajorAxis*(1-orbitElement._eccentricity*orbitElement._eccentricity)/(1 + orbitElement._eccentricity * tmpcosf);
		double h = sqrt(MIU_EARTH*(r + r * orbitElement._eccentricity*tmpcosf));
		double pt[3];
		pt[0] = r;pt[1] = 0.0;pt[2] = 0.0;
		double mat[3][3];
		iauIr(mat);
		iauRz(-f, mat);
		double ECIp[3];
		iauRxp(mat, pt, ECIp);
		iauRxp(rmat, ECIp, pt);
		iauCp(pt, ECIp);
		ECIps.push_back(ECIp);
		class vec3d ECFp;
		iauRxp(c2tmat._mat, ECIp, ECFp._vec);
		iauPmp(rpef, omgxr, ECFv._vec);
		ECFps.push_back(ECFp);
		ECFvs.push_back(ECFv);	}}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62

算法的输入包括轨道根数、星历时间、计算开始和结束时间、与采样点对应的空间坐标系变换矩阵数组（前述基于CUDA的空间坐标系变换矩阵计算得到，由于空间坐标系变换矩阵对于不同的航天器是无需重复计算的，因此合理的方法是将其单独出来，而不是将算法耦合在一起）。输出是计算得到的惯性系和地固系坐标数组（速度矢量的计算略）。

算法首先计算一些共用的数据，如开始和结束时间的儒略时、周期，计算前述轨道坐标系中的矢量变换到惯性坐标系所需的后3个变换矩阵并级联。

然后循环处理每个采样点，此时只需计算平近点角，然后转化为真近点角，计算得到轨道坐标系中的位置矢量，然后再变换到惯性坐标系，再进一步变换到地固坐标系。

可以看出，上述过程中可以并行的部分主要是针对各采样点计算部分，公共部分应该提前在CPU端计算。下面也将分别讨论基于Thrust和Runtime的并行计算方法。