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

前言

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

计算过程

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

利用卫星动量矩计算轨道倾角和升交点赤径
计算卫星相对于地心的动量矩，该动量矩等于卫星地心矩矢量和速度矢量的矢积：h=r×v\textbf{h}=\textbf{r}×\textbf{v}h=r×v，动量矩的方向和卫星轨道面的法线是平行的，动量矩和Z轴夹角为轨道倾角iii，轨道平面和地球赤道平面的交线为节线ON；节线ON与X轴夹角为升交点赤径Ω\OmegaΩ，(i,Ω)(i,\Omega)(i,Ω)确定了轨道平面在空间坐标系中的方位。
i=arccos(hx/h),Ω=arctan(−hx/hy)i=arccos(h_x/h), \Omega=arctan(-h_x/h_y)i=arccos(hx/h),Ω=arctan(−hx/hy)
利用卫星机械能计算轨道半长轴
E=v2/2−μ/r,E=−μ/2aE=v^2/2-\mu/r, E=-\mu/2aE=v2/2−μ/r,E=−μ/2a
其中hhh为动量矩模值，μ\muμ为引力常量：398600.44
km3/s2{km^3}/s^2km3/s2，vvv为速度矢量模值，rrr为位置矢量模值，aaa为椭圆轨道半长轴。
利用轨道半通经和轨道半长轴计算椭圆轨道偏心率
p=h2/μ,e=(1−(p/a))p=h^2/\mu, e=\sqrt{(1-(p/a))}p=h2/μ,e=(1−(p/a))
其中，ppp为半通径，eee为偏心率。
利用偏心率、半通经和位置矢量模值计算真近点角
f=arccos(p−r)/ref=arccos{(p-r)/re}f=arccos(p−r)/re
利用真近点角和升交点幅角计算近地点辐角
ω=u−f,u=arccos(ON⋅r/(r∗ON))\omega=u-f, u=arccos(\textbf{ON} \cdot \textbf{r} /({r*ON}))ω=u−f,u=arccos(ON⋅r/(r∗ON))
其中，升交点幅角为节线ON矢量与卫星位置矢量的夹角。
ON=(cosΩ,sinΩ,0)\textbf{ON} =(cos\Omega,sin\Omega,0)ON=(cosΩ,sinΩ,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{ };       //轨道倾角:degdouble RAAN{ };             //升交点赤经:deg    :计算结果差180double semimajorAxis{};     //半长轴:km  6917.21double Eccentricity{};    //偏心率:double argumentOfPerigee{};  //近地点辐角:deg  :double trueAnomaly{  };  //真近点角 ：主要考虑什么时候要对称变换，因为acos只能输出0—pi，而目标区间范围0—2pi
};
struct Motion {//J2000double location[3]{};   //位置:x、y、z  kmdouble 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—2piopa.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; //卫星的机械能Eopa.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—2piopa.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（矢量的矢积）//J2000double location[3]={ sat_x_g,sat_y_g,sat_z_g };   //位置:x、y、z  kmdouble 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—2piopa.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; //卫星的机械能Eopa.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—2piopa.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;}

总结

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

根据卫星运动矢量计算轨道六根数相关推荐

轨道六根数的含义汇总
文章目录轨道六根数概述其他表示 (1)半长轴椭圆抛物线双曲线与速度位置的转化椭圆双曲线 (2)离心率离心率标量椭圆抛物线双曲线离心率矢量 (3)轨道倾角 (4)近心点辐角 ...
轨道六根数 matlab,轨道六根数
在二体问题中,轨道根数(orbital factors)是描述物体运动轨迹的简便形式.三维空间中,唯一确定物体轨迹需要六个参数,如位置矢量和速度矢量(均为三维)可共同确定物体轨迹.此外,用六个轨道根数 ...
用c++根据轨道六根数计算卫星位置
轨道六根数是描述卫星轨道的一组参数,包括: 轨道长半径(a):卫星轨道的半径,表示卫星到地球中心的平均距离. 轨道偏心率(e):卫星轨道的偏心率,表示轨道的椭圆程度. 轨道倾角(i):卫星轨道与地球赤 ...
TLE两行数与轨道六根数转换
TLE与轨道六根数转换方法一.TLE格式讲解二.轨道六根数三.TLE与六根数转换 1."每天环绕地球的圈数"与"轨道半长轴"转换 2.平近点角与真近点角的 ...
GNSS之轨道六根数及常见轨道类型
文章目录 GNSS之轨道六根数及常见轨道类型一.轨道六根数二.轨道类型 1.与赤道面成64°角的椭圆轨道 2.圆形LEO 3.圆形MEO 4.地球同步轨道 5.临界倾斜轨道 Critically ...
cesium中轨道六根数的参数命名
cesium中轨道六根数的参数命名: Semimajor Axis(SMA) 半长轴:是椭圆长轴的一半.对于圆,也就是半径,另外根据开普勒第三定律,半长轴与运行周期之间有确定的换算关系. Eccent ...
MATLAB与STK互联28：仿真案例3—读取轨道六根数（DataProviders使用示例）
这个案例比较简单,生成20个卫星,然后获取第一个轨道历元的卫星轨道六要素数据. 约束:卫星轨道高度1000~6000km,偏心率0-0.2,轨道倾角0-50°,升交点赤经0-360,近地点辐角.真近点 ...
轨道六根数（开普勒六参数）
名称描述半长轴椭圆轨道长轴的一半,有时可视作平均轨道半径. 离心率椭圆轨道两焦点距离与长轴长度的比值,是椭圆轨道扁平程度的一种量度. 轨道倾角行星轨道面与黄道面的倾角. 升交点赤经 / 黄经 ...
航天器轨道六要素和TLE两行轨道数据格式
航天器轨道要素椭圆轨道六根数指的是:半长轴 a a a,离心率e,轨道倾角 i i i.升交点赤经 Ω \Omega Ω.近地点辐角 ω \omega ω.和过近地点时刻 t 0 t_0 t0(或 ...
java版通过轨道6根数实现计算出经纬度坐标
近期公司有个项目,实现卫星六根数实现计算出经纬度坐标,因为在网上找不到java资源.翻阅了大量的文章.然后自己着手去根据matlab代码实现了java语言的转换. 卫星轨道6根数主要有半长轴a.离心率 ...

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

前言

计算过程

代码实现

运行结果

更新

总结

根据卫星运动矢量计算轨道六根数相关推荐

最新文章

热门文章