卫星轨道推演计算相关知识点总结（含欧拉角、旋转矩阵、及各坐标系转化等）

码代码的乔木

已于 2024-07-30 15:02:16 修改

阅读量1.6w

点赞数 27

分类专栏： C++ 文章标签：仿真轨道 C++

于 2021-11-13 15:42:23 首次发布

本文链接：https://blog.csdn.net/nannan7777/article/details/119752362

版权

C++ 专栏收录该内容

1 篇文章

订阅专栏

本文总结了卫星轨道推演计算的相关知识点，包括卫星的运动特性，地心惯性坐标系、地心固定坐标系、地心地固坐标系和地心椭球惯性坐标系的定义，以及如何根据轨道根数计算卫星位置。同时，介绍了欧拉角、旋转向量和旋转矩阵的相互转换在C++中的实现。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

卫星轨道推演计算相关基础知识点总结

近期在补卫星仿真领域的相关基础知识，并总结一些卫星轨道计算的算法与工具类代码函数。
参考内容见百度公开文库文档与PPT：
https://max.book118.com/html/2019/0309/5304232023002020.shtm
https://wenku.baidu.com/view/e231b968a32d7375a417806b.html？rec_flag=default&fr=Recommend_RelativeDoc-60272,60321,40155,40311,40251,60314,40300-kpdrec_doc_pc_view-bf93f0936bec0975f465e2dd&sxts=1629169656038

一、卫星的运动特性

在这里插入图片描述

二、卫星的空间坐标系

（1）地心惯性坐标系(ECI×J2000历元坐标系)(ECSF 地心空间坐标系)

地心惯性坐标系是太阳系内的一个惯性坐标系，不随地球而转动，也不受地球、DAC6574IDGS太阳运行的章动和岁差的影响。其坐标原点位于地心Oe；OeX轴位于赤道平面内，指向特定某一年（历元时刻）的太阳春分点位置，每年春分点均会发生变动；OeZ轴指向某一年（历元时刻）地球北极的平均位置处，即地球平均自转极点（CIO）；OcY轴位于赤道平面内，与0eX轴垂直，且与OeX、OeZ构成满足右手定则的笛卡儿直角坐标系。
由于采用的历元时间不同，可以有各种不同的地心惯性坐标系，目前国际上通用的地心惯性坐标系是J2000历元坐标系，它是以公元2000年的春分点为基准的历元坐标系。

（2）地心固定坐标系(ECF)(WGS84坐标系)

如图所示，地心固定坐标系的坐标原点位于地心Oe，OeZ轴指向地球北极，OeX轴位于赤道平面内指向地理经度的零点，OeY轴根据右手定则确定。地心固定坐标系为笛卡尔直角坐标系，该坐标系在宇宙空间中相对地球静止，伴随地球自转和公转。
在这里插入图片描述

（3）地心地固坐标系(ECEF)

地心地固坐标系的坐标原点位于地心Oe，OeZ轴指向地球平均自转极点（CIO），OeX轴位于赤道平面内指向子午线与赤道交点，OeY轴根据右手定则确定。地心地固坐标系为笛卡尔直角坐标系，该坐标系在宇宙空间中相对地球静止，伴随地球自转和公转。

（4）地心椭球惯性坐标系(ECEI)（地心黄道惯性坐标系）

地心椭球惯性坐标系的坐标原点位于地心Oe，OeZ轴指向椭球极轴（黄道极-地球公转轨迹在地球表面投影轨迹形成），OeX轴位于赤道平面内，指向特定某一年（历元时刻）的太阳春分点位置，每年春分点均会发生变动，OeY轴根据右手定则确定。地心地固坐标系为笛卡尔直角坐标系，该坐标系在宇宙空间中相对地球静止，伴随地球自转和公转。

三、根据轨道根数计算卫星位置

（1）计算卫星在轨道坐标系下的位置

轨道根数（或称轨道要素或轨道参数）是描述在牛顿运动定律和牛顿万有引力定律的作用下的天体或航天器，在其开普勒轨道上运动时，确定其轨道所必要的六个参数。由于运动的方式有许多种的参数表示法，依照选定的测量装置不同，对相同的轨道，有几种不同的方式来定义轨道根数。

轨道半长轴：椭圆轨道长轴的一半，有时可视作平均轨道半径。
轨道离心率：为椭圆扁平程度的一种量度，定义是椭圆两焦点间的距离与长轴长度的比值。就是e=c/a。
轨道倾角：行星轨道面对黄道面的倾角；在升交点处从黄道面逆时针方向量到行星轨道面的角度。
升交点赤经：行星轨道升交点的黄道经度。
近日点幅角：从升交点沿行星运动轨道逆时针量到近日点的角度。
平近点角：行星对应于t0时卫星的平近点角。
在这里插入图片描述

（2）轨道坐标系与大地坐标系间的换算

在这里插入图片描述

（3）轨道坐标系与大地坐标系间的换算

在这里插入图片描述

四、空间几何关系相关知识点及C++代码

（1）欧拉角、旋转向量和旋转矩阵的相互转换

在这里插入图片描述
任何一个旋转可以表示为依次绕着三个旋转轴旋三个角度的组合，这三个角度成为欧拉角。对于在三维空间里的一个参考系，任何坐标系的取向都可以用欧拉角来表示，如图所示，蓝色为起始坐标系，红色为旋转之后的坐标系。
在这里插入图片描述
因此欧拉角转为旋转矩阵的公式如下：

欧拉角和旋转矩阵相互转换的C++代码如下所示：

Mat eulerAnglesToRotationMatrix(Vec3f &theta)
{
    // Calculate rotation about x axis
    Mat R_x = (Mat_<double>(3,3) <<
        1,       0,              0,
        0,       cos(theta[0]),   -sin(theta[0]),
        0,       sin(theta[0]),   cos(theta[0])
    );
    // Calculate rotation about y axis
    Mat R_y = (Mat_<double>(3,3) <<
        cos(theta[1]),    0,      sin(theta[1]),
        0,               1,      0,
        -sin(theta[1]),   0,      cos(theta[1])
    );
    // Calculate rotation about z axis
    Mat R_z = (Mat_<double>(3,3) <<
        cos(theta[2]),    -sin(theta[2]),      0,
        sin(theta[2]),    cos(theta[2]),       0,
        0,               0,                  1
    );
    // Combined rotation matrix
    Mat R = R_z * R_y * R_x;
    return R;
}
Vec3f rotationMatrixToEulerAngles(Mat &R)
{
    float sy = sqrt(R.at<double>(0,0) * R.at<double>(0,0) +  R.at<double>(1,0) * R.at<double>(1,0) );
    bool singular = sy < 1e-6; // If
    float x, y, z;
    if (!singular)
    {
        x = atan2(R.at<double>(2,1) , R.at<double>(2,2));
        y = atan2(-R.at<double>(2,0), sy);
        z = atan2(R.at<double>(1,0), R.at<double>(0,0));
    }
    else
    {
        x = atan2(-R.at<double>(1,2), R.at<double>(1,1));
        y = atan2(-R.at<double>(2,0), sy);
        z = 0;
    }
    #if 1
    x = x*180.0f/3.141592653589793f;
    y = y*180.0f/3.141592653589793f;
    z = z*180.0f/3.141592653589793f;
    #endif
    return Vec3f(x, y, z);
}

这里存一下项目开发过程中常用的一些计算需求，例如根据旋转前后的两个向量求二者之间的旋转矩阵，C++代码如下所示：

AcGeVector3d CrossProduct(AcGeVector3d a, AcGeVector3d b)
{
	AcGeVector3d c;
	c.x = a.y * b.z - a.z * b.y;
	c.y = a.z * b.x - a.x * b.z;
	c.z = a.x * b.y - a.y * b.x;
	return c;
}

double DotProduct(AcGeVector3d a, AcGeVector3d b)
{
	double result;
	//result = a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
	result = a.x * b.x + a.y * b.y + a.z * b.z;
	return result;
}

double Absolute(AcGeVector3d v)
{
	double result;
	result = sqrt(v.x * v.x + v.y * v.y + v.z * v.z);
	return result;
}

AcGeMatrix3d RotationMatrix(double angle, AcGeVector3d u)
{
	double norm = Absolute(u);
	u.x = u.x / norm;
	u.y = u.y / norm;
	u.z = u.z / norm;
	AcGeMatrix3d rotatinMatrix;

	rotatinMatrix.entry[0][0] = cos(angle) + u.x * u.x * (1 - cos(angle));
	rotatinMatrix.entry[0][1] = u.x * u.y * (1 - cos(angle)) - u.z * sin(angle);
	rotatinMatrix.entry[0][2] = u.y * sin(angle) + u.x * u.z * (1 - cos(angle));
	rotatinMatrix.entry[0][3] = 0.0;

	rotatinMatrix.entry[1][0] = u.z * sin(angle) + u.x * u.y * (1 - cos(angle));
	rotatinMatrix.entry[1][1] = cos(angle) + u.y * u.y * (1 - cos(angle));
	rotatinMatrix.entry[1][2] = -u.x * sin(angle) + u.y * u.z * (1 - cos(angle));
	rotatinMatrix.entry[1][3] = 0.0;

	rotatinMatrix.entry[2][0] = -u.y * sin(angle) + u.x * u.z * (1 - cos(angle));
	rotatinMatrix.entry[2][1] = u.x * sin(angle) + u.y * u.z * (1 - cos(angle));
	rotatinMatrix.entry[2][2] = cos(angle) + u.z * u.z * (1 - cos(angle));
	rotatinMatrix.entry[2][3] = 0.0;

	rotatinMatrix.entry[3][0] = 0.0;
	rotatinMatrix.entry[3][1] = 0.0;
	rotatinMatrix.entry[3][2] = 0.0;
	rotatinMatrix.entry[3][3] = 1.0;

	return rotatinMatrix;
}

AcGeMatrix3d CalculationMtrix(AcGeVector3d vectorBefore, AcGeVector3d vectorAfter)
{
	double rotationAngle = acos(DotProduct(vectorBefore, vectorAfter) / Absolute(vectorBefore) / Absolute(vectorAfter));
	AcGeVector3d rotationAxis = CrossProduct(vectorBefore, vectorAfter);
	AcGeMatrix3d rotationMatrix = RotationMatrix(rotationAngle, rotationAxis);

	return rotationMatrix;
}