根据卫星运动矢量计算轨道六根数

前言

STK软件在给定六根数时，可求得卫星位置和速度矢量，但有时我们通过星历参数得到卫星的位置和速度矢量，希望能够反演得出卫星轨道的六根数。从而方便对该卫星轨道进行仿真模拟。

计算过程

给定卫星在J2000坐标系下的的位置矢量r和速度矢量v

1. 利用卫星动量矩计算轨道倾角和升交点赤径
   计算卫星相对于地心的动量矩，该动量矩等于卫星地心矩矢量和速度矢量的矢积：$\text{h}=\text{r}\times \text{v}$，动量矩的方向和卫星轨道面的法线是平行的，动量矩和Z轴夹角为轨道倾角$i$，轨道平面和地球赤道平面的交线为节线ON；节线ON与X轴夹角为升交点赤径$\Omega$，$(i,\Omega )$确定了轨道平面在空间坐标系中的方位。
   $$i=arccos(h_x/h),\Omega =arctan(-h_x/h_y)$$
2. 利用卫星机械能计算轨道半长轴
   $$E=v^2/2-\mu /r,E=-\mu /2a$$
   其中$h$为动量矩模值，$\mu$为引力常量：398600.44
   $k{m}^3/s^2$，$v$为速度矢量模值，$r$为位置矢量模值，$a$为椭圆轨道半长轴。
3. 利用轨道半通经和轨道半长轴计算椭圆轨道偏心率
   $$p=h^2/\mu ,e=\sqrt{(1-(p/a))}$$
   其中，$p$为半通径，$e$为偏心率。
4. 利用偏心率、半通经和位置矢量模值计算真近点角
   $$f=arccos(p-r)/re$$
5. 利用真近点角和升交点幅角计算近地点辐角
   $$\omega =u-f,u=arccos(\text{ON}\cdot \text{r}/(r\ast ON))$$
   其中，升交点幅角为节线ON矢量与卫星位置矢量的夹角。
   $\text{ON}=(cos\Omega ,sin\Omega ,0)$。

代码实现

具体计算时，需要考虑反三角函数的值域与实际情况对应。

#include <string>
#include <iostream>
#include <iomanip>
#include <cmath>
using namespace std;
const double mu{ 398600.44 };//引力常数:(km)3/s2
const double PI{ acos(-1) };//PI
const double rad2deg{ 180.0 / PI };//PI
#define ABS(x) (sqrt((x)[0]*(x)[0]+(x)[1]*(x)[1]+(x)[2]*(x)[2]))
struct OrbitParm {
    double inclination{ };       //轨道倾角:deg
    double RAAN{ };             //升交点赤经:deg    :计算结果差180
    double semimajorAxis{};     //半长轴:km  6917.21
    double Eccentricity{};    //偏心率:

    double argumentOfPerigee{};  //近地点辐角:deg  :
    double trueAnomaly{  };  //真近点角 ：主要考虑什么时候要对称变换，因为acos只能输出0—pi，而目标区间范围0—2pi
};
struct Motion {//J2000
    double location[3]{};   //位置:x、y、z  km
    double speed[3]{};                //速度:x、y、z  km/sec  
};
OrbitParm motionOrbitParmConvert(Motion mot) {//to be done
    //卫星相对于地心的动量矩:h=r*v（矢量的矢积）
    double h[3]{ mot.location[1] * mot.speed[2] - mot.location[2] * mot.speed[1], \
    - mot.location[0] * mot.speed[2] + mot.location[2] * mot.speed[0], \
    mot.location[0] * mot.speed[1] - mot.location[1] * mot.speed[0] };
    double absH{ ABS(h) };
    OrbitParm opa{};
    opa.inclination = acos(h[2] / absH) * rad2deg;
    opa.RAAN = atan2(h[0], -h[1]) * rad2deg;
    if ((opa.RAAN) < 0)//目标区间为0—2pi
        opa.RAAN = opa.RAAN + 360;
    double p{ absH * absH / mu };     //椭圆轨道的半通径
    double absR{ ABS(mot.location) };
    double absV{ ABS(mot.speed) };
    double E = absV * absV / 2.0 - mu / absR; //卫星的机械能E
    opa.semimajorAxis = -mu / E / 2.0;//半长轴由机械能决定
    opa.Eccentricity = sqrt(1 - p / opa.semimajorAxis);//偏心率可通过半长轴和半通径联合求得
    opa.trueAnomaly = acos((p - absR) / absR / opa.Eccentricity) * rad2deg;
    if (1)//主要考虑什么时候要对称变换，因为acos只能输出0—pi，而目标区间范围0—2pi
        opa.trueAnomaly = 360 - opa.trueAnomaly;
    double u[3] = { cos(opa.RAAN / rad2deg),sin(opa.RAAN / rad2deg),0 };
    opa.argumentOfPerigee = acos((u[0] * mot.location[0] + u[1] * mot.location[1]) / absR) * rad2deg;
    opa.argumentOfPerigee -= opa.trueAnomaly;
    if (opa.argumentOfPerigee < 0)
        opa.argumentOfPerigee += 360;
    return opa;
}
Motion motionOrbitParmConvert(OrbitParm opa) {//暂时不编
    return{};
}
int main()
{
    //输入示例，第一个大括号依次填入J2000坐标系下的xyz位置，第二个括号依次填入J2000坐标系下的xyz速度
    auto opa = motionOrbitParmConvert({ {-3904.3,-4663.0,3290.863664} , {1.4,3.4,6.6} });
    cout <<right <<fixed << setprecision(6)<<setfill('0');
    cout << setw(11) << opa.inclination << endl;
    cout << setw(11) << opa.Eccentricity << endl;
    cout << setw(11) << opa.semimajorAxis << endl;
    cout << setw(11) << opa.RAAN << endl;
    cout << setw(11) << opa.trueAnomaly << endl;
    cout << setw(11) << opa.argumentOfPerigee << endl;
    return 0;
}

运行结果

更新

除了计算六根数，还计算了平近点角、偏近点角

class Orbit_Para_Object
{
public:
	//卫星半长轴
	double dOrbit_a;

	//计算轨道偏心率
	double dOrbit_e;

	//计算轨道偏心角
	double dOrbit_E1;

	//计算真近心角
	double dOrbit_Theta;

	//计算平均近心角
	double dOrbit_M;
	//计算轨道倾角
	double dOrbit_Angle_Inclination;
	//升交点赤经
	double dOrbit_Angle_Omig;

	//近地点幅角
	double dOrbit_Angle_W;

protected:

private:

};


#define ABS(x) (sqrt((x)[0]*(x)[0]+(x)[1]*(x)[1]+(x)[2]*(x)[2]))
const double rad2deg{ 180.0 / pi };//PI
Orbit_Para_Object Cal_orbit_info(double sat_x_g, double sat_y_g, double sat_z_g, double sat_vx_g, double sat_vy_g, double sat_vz_g, double Gravitation_P) {//to be done
	
	//  //卫星相对于地心的动量矩:h=r*v（矢量的矢积）
	//J2000
	double location[3]={ sat_x_g,sat_y_g,sat_z_g };   //位置:x、y、z  km
	double speed[3]={ sat_vx_g,sat_vy_g,sat_vz_g };                //速度:x、y、z  km/sec  

	double h[3]={location[1] * speed[2] - location[2] * speed[1], \
	- location[0] * speed[2] + location[2] * speed[0], \
	location[0] * speed[1] - location[1] * speed[0] };
	double absH{ ABS(h) };
	Orbit_Para_Object opa{};
	opa.dOrbit_Angle_Inclination = acos(h[2] / absH) * rad2deg;
	opa.dOrbit_Angle_Omig = atan2(h[0], -h[1]) * rad2deg;
	if ((opa.dOrbit_Angle_Omig) < 0)//目标区间为0—2pi
		opa.dOrbit_Angle_Omig = opa.dOrbit_Angle_Omig + 360;
	double p = { absH * absH / Gravitation_P };     //椭圆轨道的半通径
	double absR = { ABS(location) };
	double absV = { ABS(speed) };
	double E = absV * absV / 2.0 - Gravitation_P / absR; //卫星的机械能E
	opa.dOrbit_a = -Gravitation_P / E / 2.0;//半长轴由机械能决定
	opa.dOrbit_e = sqrt(1 - p / opa.dOrbit_a);//偏心率可通过半长轴和半通径联合求得
	opa.dOrbit_Theta = acos((p - absR) / absR / opa.dOrbit_e) * rad2deg;
	if (1)//主要考虑什么时候要对称变换，因为acos只能输出0—pi，而目标区间范围0—2pi
		opa.dOrbit_Theta = 360 - opa.dOrbit_Theta;
	double u[3] = { cos(opa.dOrbit_Angle_Omig / rad2deg),sin(opa.dOrbit_Angle_Omig / rad2deg),0 };
	opa.dOrbit_Angle_W = acos((u[0] * location[0] + u[1] * location[1]) / absR) * rad2deg;
	opa.dOrbit_Angle_W -= opa.dOrbit_Theta;
	if (opa.dOrbit_Angle_W < 0)
		opa.dOrbit_Angle_W += 360;

	double n=sqrt( Gravitation_P/( opa.dOrbit_a * opa.dOrbit_a * opa.dOrbit_a));//卫星沿椭圆轨道运行的平均速率

	//计算偏近点角
	//opa.dOrbit_E1 = acos(absR * cos(opa.dOrbit_Theta) / opa.dOrbit_a + opa.dOrbit_e);
	opa.dOrbit_E1 = atan2(sqrt(1-opa.dOrbit_e* opa.dOrbit_e* opa.dOrbit_e) *sin(opa.dOrbit_Theta)/(1+opa.dOrbit_e*cos(opa.dOrbit_Theta)), (opa.dOrbit_e + cos(opa.dOrbit_Theta)) / (1 + opa.dOrbit_e * cos(opa.dOrbit_Theta)));
	//计算平近点角
	opa.dOrbit_M =fmod((opa.dOrbit_E1- opa.dOrbit_e*sin(opa.dOrbit_E1)) * rad2deg,360 );
	opa.dOrbit_E1 *= rad2deg;
	if (opa.dOrbit_E1 < 0)
		opa.dOrbit_E1 += 360;

	if (opa.dOrbit_M < 0)
		opa.dOrbit_M += 360;
	return opa;

}

总结

该文实现了通过卫星星历参数反演得出卫星轨道的六根数。