根据卫星运动矢量计算轨道六根数-Toy模板网

这篇具有很好参考价值的文章主要介绍了根据卫星运动矢量计算轨道六根数。希望对大家有所帮助。如果存在错误或未考虑完全的地方，请大家不吝赐教，您也可以点击"举报违法"按钮提交疑问。

前言

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

计算过程

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

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

}