高精度轨道传播的国际空间站

这个示例使用:

这个例子展示了如何传播国际空间站(ISS)的轨道使用高精度数字轨道传播与航空航天Blockset™。它使用轨道传播算子块计算24小时国际空间站轨道。例子然后比较位置和速度州公开空间站轨道数据可以从美国国家航空和宇宙航行局轨道操作和规划(威尼斯平底渔船)飞行控制器在约翰逊航天中心。

航天器的轨道传播算子块模型转化动态使用数值积分。计算一个或多个航天器的位置和速度。最准确的结果,使用一个变量一步解决公差设置较低(小于1 e-9)。根据任务要求,您可以使用更大的公差的增加速度。这样做可能会影响解决方案的准确性。传播轨道状态,阻止使用引力模型选择对当前中央的身体。块还包括大气阻力(地球轨道),和外部扰动加速度作为输入提供。

在这个例子中,使用国际空间站的天球参考系(ICRF)位置和速度状态直接从美国宇航局引用公共分配文件[1]来初始化一个空间站轨道的轨道传播算子。使用EGM2008球面谐波重力模型来模拟地球引力,和NRLMSISE-00大气模型来计算大气阻力。计算第三身体重力扰动(月球、太阳和其他行星在太阳系)和太阳辐射压力(SRP)以外的轨道传播算子块并通过结果成块外部应用加速度。24小时运行模型,并针对发布的美国宇航局轨迹数据的比较。

设置任务初始条件

这个示例使用MATLAB®结构组织任务数据。这些结构在稍后的示例中访问数据更直观。他们还帮助整理全球基地工作区。

引用ISS星历数据文件[1]由NASA提供1月3日开始,2022年12:00:00.000 UTC。

的使命。年代tartDate = datetime(2022, 1, 3, 12, 0, 0); mission.Duration = hours(24);

文件包含ICRF位置(公里)和速度(公里/秒)数据采样每4分钟。

国际空间站。X0_icrf = (-1325.896391725290 5492.890955896010 3762.423747679220);%公里国际空间站。V0_icrf = (-4.87470128630892 -4.10251688094599 4.26428812476909);%公里/秒

它还包含质量(公斤),拖动区域(m ^ 2)和阻力系数与时代相应的数据。

国际空间站。质量= 459023.0;%公斤国际空间站。Dr一个gArea = 1951.0;% m ^ 2国际空间站。Dr一个gCoeff = 2.0;

该文件还包括太阳辐射面积(m ^ 2)和太阳辐射系数;然而这些值是零,因为SRP并不突出的低地球轨道(LEO)国际空间站运行的地方。尽管他们微不足道的影响产生的轨迹,我们将包括太阳辐射压力计算在这个例子充分证明高精度轨道传播。SRP更高轨道政权更突出。

国际空间站。年代RP一个re一个=1500;% m ^ 2国际空间站。年代RPCoeff = 1.8; iss.P_sr = 4.5344321e-6;% N / m ^ 2

高精度轨道传播算法

这个例子执行高精度数值轨道传播使用考威尔的方法。惯性(ICRF)框架,地球的重力加速度是总结所有加速度扰动和double-integrated找到速度和位置在每个时间步的仿真。

${\vec{一个}}_{icrf} = {\vec{一个}}_{重力} + {\vec{一个}}_{拖} + {\vec{一个}}_{3 理查德·道金斯身体} + {\vec{一个}}_{srp}$

${\vec{一个}}_{icrf} \underset{集成}{⟹} {\vec{v}}_{icrf}, {\vec{r}}_{icrf}$

地点:

${\vec{一个}}_{重力}$ 是地球重力加速度。

${\vec{一个}}_{拖}$ 是大气阻力加速度。

${\vec{一个}}_{3 理查德·道金斯身体} = {\vec{一个}}_{月亮} + {\vec{一个}}_{太阳} + {\vec{一个}}_{汞} + {\vec{一个}}_{金星} + {\vec{一个}}_{火星} + {\vec{一个}}_{木星} + {\vec{一个}}_{土星} + {\vec{一个}}_{天王星} + {\vec{一个}}_{海王星}$ 月球的重力加速度,太阳和行星的太阳系。

${\vec{一个}}_{srp}$ 加速度是由于太阳辐射压力。

地球引力

地球的引力模型在本例中是EGM2008球面谐波重力模型。这个模型占区域、部门和田形谐波。供参考,二级零阶带谐函数 $J_{2}$ ,占地球的扁率 $- - - - - - C_{2, 0}$ 。球面谐波模型谐波占到最大程度 $l = l_{马克斯}$ 。EGM2008, $l_{马克斯} = 2159年$ 。

球面谐波计算重力在固定框架(ff)协调系统(国际地球参考架(ITRF),在地球的情况下)。然而,数值积分,总是惯性ICRF坐标系统中执行。因此,在每一个步伐,位置和速度状态转换为固定坐标系,计算重力在固定框架,以及由此产生的加速度转换成惯性坐标系。

${\vec{一个}}_{重力} = - - - - - - \frac{μ}{r^{3}} \vec{r_{icrf}} + ff2icrf ({\vec{一个}}_{nonspherical})$

地点:

$\begin{array}{l} {\vec{一个}}_{nonspherical} = {(\frac{1}{r} \frac{\partial}{\partial r} U - - - - - - \frac{{r_{ff}}_{k}}{r^{2} \sqrt{{r_{ff}}_{我}^{2} + {r_{ff}}_{j}^{2}}} \frac{\partial}{\partial ϕ} U] {r_{ff}}_{我} - - - - - - (\frac{1}{{r_{ff}}_{我}^{2} + {r_{ff}}_{j}^{2}} \frac{\partial}{\partial λ} U] {r_{ff}}_{j}} 我 \\ + {(\frac{1}{r} \frac{\partial}{\partial r} U + \frac{{r_{ff}}_{k}}{r^{2} \sqrt{{r_{ff}}_{我}^{2} + {r_{ff}}_{j}^{2}}} \frac{\partial}{\partial ϕ} U] {r_{ff}}_{j} + (\frac{1}{{r_{ff}}_{我}^{2} + {r_{ff}}_{j}^{2}} \frac{\partial}{\partial λ} U] {r_{ff}}_{我}} j \\ + {\frac{1}{r} (\frac{\partial}{\partial r} U) {r_{ff}}_{k} + \frac{\sqrt{{r_{ff}}_{我}^{2} + {r_{ff}}_{j}^{2}}}{r^{2}} \frac{\partial}{\partial ϕ} U} k \end{array}$

鉴于以下偏导数在球坐标:

$\begin{array}{l} \frac{\partial}{\partial r} U = - - - - - - \frac{μ}{r^{2}} \sum_{l = 2}^{l_{马克斯}} \sum_{米 = 0}^{l} {(\frac{R_{cb}}{r})}^{l} (l + 1) P_{l, 米} (罪 (ϕ)] {C_{l, 米} 因为 (米 λ) + {年代}_{l, 米} 罪 (米 λ)} \\ \frac{\partial}{\partial ϕ} U = \frac{μ}{r} \sum_{l = 2}^{l_{马克斯}} \sum_{米 = 0}^{l} {(\frac{R_{cb}}{r})}^{l} {P_{l, 米 + 1} (罪 (ϕ)] - - - - - - (米) 棕褐色 (ϕ) P_{l, 米} (罪 (ϕ)]} {C_{l, 米} 因为 (米 λ) + {年代}_{l, 米} 罪 (米 λ)} \\ \frac{\partial}{\partial λ} U = \frac{μ}{r} \sum_{l = 2}^{l_{马克斯}} \sum_{米 = 0}^{l} {(\frac{R_{cb}}{r})}^{l} (米) P_{l, 米} (罪 (ϕ)] {{年代}_{l, 米} 因为 (米 λ) - - - - - - C_{l, 米} 罪 (米 λ)} \end{array}$

$P_{l, 米}$ 勒让德函数相关联。

$C_{l, 米}$ 和 ${年代}_{l, 米}$ 非规范谐波系数。

大气阻力

轨道传播算子块支持包容大气阻力加速度使用NRLMSIS万博1manbetxE-00大气模型。狮子座和大气阻力是一个占主导地位的扰动导致国际空间站轨道随时间降低没有援助轨道机动。

大气阻力计算为:

${\vec{一个}}_{拖} = - - - - - - \frac{1}{2} ρ (\frac{C_{D} {一个}_{D}}{米}) {v_{rel}}^{2}$

地点:

$ρ$ 是大气密度。

$C_{D}$ 是宇宙飞船阻力系数。

${一个}_{D}$ 宇宙飞船阻力区域,或该地区正常吗 $v_{rel}$ 。

$米$ 航天器的质量。

${\vec{v}}_{rel} = {\vec{v}}_{icrf} - - - - - - Ω \times {\vec{r}}_{icrf}$ 是宇宙飞船速度相对于大气中。

NRLMSISE-00大气模型是用来计算大气密度。它需要空间天气数据(F10.7和F10.7A广播流量值和地磁指数),可以获得一个统一格式从CelesTrak®。关于空间天气数据的更多信息在MATLAB和Simulink仿真®,看到太阳能通量和地磁索引块参考页面和航空航天工具箱万博1manbetxaeroReadSpaceWeatherData函数引用页面。

第三个身体重力

第三个身体重力是包含在轨道传播通过总结加速度由于每个天体的轨道传播算子块使用应用加速度( ${一个}_{icrf}$ )端口。为每个身体计算重力(太阳)所示的例子:

${\vec{一个}}_{⊙} = μ_{⊙} (\frac{{\vec{r}}_{坐, ⊙}}{{r_{坐, ⊙}}^{3}} - - - - - - \frac{{\vec{r}}_{\oplus, ⊙}}{{r_{\oplus, ⊙}}^{3}})$

地点:

$μ_{⊙}$ 太阳的引力参数。

${\vec{r}}_{坐, ⊙}$ 是向量从宇宙飞船到太阳的中心,基于喷气推进实验室的行星DE405星历数据。关于行星的星历表的更多信息,请访问行星历块引用页面或航空航天工具箱planetEphemeris函数引用页面。

${\vec{r}}_{\oplus, ⊙}$ 是向量从地球中心太阳的中心,基于喷气推进实验室DE405星历数据。

太阳辐射压力

太阳辐射压力在狮子座轨道传播影响很小,但是这个例子的计算包括完整性。太阳辐射压力计算为:

${\vec{一个}}_{srp} = - - - - - - ν \frac{C_{r} {一个}_{年代}}{米} \frac{p_{老}}{{r_{坐, ⊙}}^{3}} {\vec{r}}_{坐, ⊙}$

地点:

$ν$ 阴影部分。暗影的值等于0,在半影在0和1之间,否则等于1。Eclipse模型不同程度的忠诚可以用来确定这个值根据任务要求。

$C_{r}$ 是反射系数。

${一个}_{年代}$ 宇宙飞船SRP区域,或太阳所看到的横截面积。

$米$ 航天器的质量。

$p_{老} = \frac{1353年}{3 e8} \frac{W / 米^{2}}{米 / 年代} = 4 。 5344321 e^{- - - - - - 6} \frac{N}{米^{2}}$ 是太阳能通量辐射压力。这个值与太阳能通量对应1 au除以光速。

${\vec{r}}_{坐, ⊙}$ 是向量从航天器太阳中心,基于喷气推进实验室DE405星历数据。

打开轨道传播模型

打开包含仿真软件模型。万博1manbetx这个模型包含一个轨道传播算子块连接到root外港阻塞和3子系统用于扰动加速度模型。该模型还包括一个“任务分析/可视化”部分包含一个仪表板回调按钮。当点击时,此按钮运行模型,创建一个新的satelliteScenario对象包含航天器在轨道中定义的基本工作空间传播算子块,并打开一个卫星场景查看器窗口。查看此操作的源代码,双击按钮回调。“任务分析/可视化”部分是一个独立的工作流来创建一个新的satelliteScenario对象并没有使用这个示例的一部分。

的使命。mdl =“ISSHighPrecisionOrbitPropagationExampleModel”;open_system (mission.mdl);

配置仿真软件模型万博1manbetx

轨道的路径传播算子定义块模型。

国际空间站。黑色=任务。mdl +“/高精度数值轨道传播算子”;

航天器模型的初始条件属性检查器,或通过编程的方式使用万博1manbetxset_param。

set_param (iss.blk…startDate可以=字符串(juliandate (mission.StartDate)),…stateFormatNum =“ICRF状态向量”,…inertialPosition = mat2str (iss.X0_icrf),…inertialVelocity = mat2str (iss.V0_icrf));

航天器轨道传播算子块物理属性。

set_param (iss.blk…质量=字符串(iss.Mass),…dragCoeff =字符串(iss.DragCoeff),…dragArea = string (iss.DragArea));

宇宙飞船在扰动部分物理性质。

的使命。srpBlk = mission.mdl +“/太阳辐射压力”;set_param(任务。srpBlk +“/ SRPCoeff”值=字符串(iss.SRPCoeff));set_param(任务。srpBlk +“/ SRPArea”值=字符串(iss.SRPArea));set_param(任务。srpBlk +“/质量”值=字符串(iss.Mass));

应用模型级解算器设置并保存模型输出端口时间表对象的数据作为一个数据集。

set_param (mission.mdl…SolverType =“变步”,…SolverName =“VariableStepAuto”,…RelTol =“1 e-12”,…AbsTol =“1 e-12”,…StopTime =字符串(秒(mission.Duration)),…SaveOutput =“上”,…OutputSaveName =“你”,…SaveFormat =“数据集”,…DatasetSignalFormat =“时间表”);

运行模型,收集星历表

模拟模型。轨道传播算子的块被配置为输出位置和速度州ICRF坐标系。仿真运行可能需要几分钟时间由于严格的公差设置在前一节中我们应用到模型(1 e-12)。

的使命。年代我米Output = sim(mission.mdl);

从模型中提取位置和速度数据输出的数据结构。

国际空间站。EphPosICRF = mission.SimOutput.yout {1} . values;国际空间站。EphVelICRF = mission.SimOutput.yout {2} . values;

设置启动任务的数据时间表对象将datetime行从持续的时间。

iss.EphPosICRF.Properties。年代tartTime = mission.StartDate; iss.EphPosICRF.Time.TimeZone =“UTC”;iss.EphVelICRF.Properties。年代tartTime = mission.StartDate; iss.EphVelICRF.Time.TimeZone =“UTC”;

星历表导入到一个`satelliteScenario`对象

创建一个卫星场景对象。

场景= satelliteScenario;

添加的航天器satelliteScenario是ICRF位置和速度时间表对象的使用卫星方法。

国际空间站。EphPosICRF_m = iss.EphPosICRF;iss.EphPosICRF_m。D一个ta = iss.EphPosICRF.Data*1e3; iss.EphVelICRF_mps = iss.EphVelICRF; iss.EphVelICRF_mps.Data = iss.EphVelICRF.Data*1e3; iss.obj = satellite(scenario, iss.EphPosICRF_m, iss.EphVelICRF_mps,…CoordinateFrame =“惯性”、名称=“国际空间站”);

的启动和停止时间场景自动调整,以反映的启动和停止时间时间表对象。

disp (“开始:“+字符串(scenario.StartTime));disp (“停止:“+字符串(scenario.StopTime));

开始:03 - 1月- 2022 12:00:00站:04 -简- 2022 12:00:00

查看器打开卫星场景。

satelliteScenarioViewer(场景);

负荷计算的参考轨迹和比较轨迹

负载空间站轨道数据参考[1]。

国际空间站。RefEphemeris = load(“ISS_ephemeris.mat”,“pos_eme2000”,“vel_eme2000”);

重新取样模拟数据匹配4分钟的采样率的参考轨迹。

国际空间站。EphPosICRF =调整时间(iss)。EphPosICRF iss.RefEphemeris.pos_eme2000.Time);国际空间站。EphVelICRF =调整时间(iss)。EphVelICRF iss.RefEphemeris.vel_eme2000.Time);

情节之间的位置误差仿真数据和参考轨迹。

情节(iss.EphPosICRF。时间,iss.EphPosICRF。D一个ta - iss.RefEphemeris.pos_eme2000.Position); title(“空间站高精度轨迹位置错误”);包含(“时间”);ylabel (“位置误差(公里)”);传奇([“\ DeltaX_ {icrf}”,“\ DeltaY_ {icrf}”,“\ DeltaZ_ {icrf}”),位置=“西北”);

仿真和引用数据之间的最大误差超过24小时大约是0.25公里。

情节(iss.EphPosICRF。时间,vecnorm (iss.EphPosICRF。D一个ta - iss.RefEphemeris.pos_eme2000.Position, 2, 2)); title(“空间站高精度轨迹位置误差大小”);包含(“时间”);ylabel (“位置误差(公里)”);

比较高精度轨迹分析传播

相比之下,添加一个新的宇宙飞船卫星场景使用Two-Body-Keplerian传播。

初始位置和速度转换为轨道元素。

[iss.OrbEl.SemiMajorAxis,…iss.OrbEl.Eccentricity,…iss.OrbEl.Inclination,…iss.OrbEl.RAAN,…iss.OrbEl.ArgPeriapsis,…iss.OrbEl。TrueAnomaly] = ijk2keplerian (iss)。X0_icrf * 1 e3, iss.V0_icrf * 1 e3);

添加一个新的空间站使用Two-Body-Keplerian传播对象的场景。

iss.KeplerianOrbit。obj =卫星(场景、iss.OrbEl.SemiMajorAxis…iss.OrbEl.Eccentricity,iss.OrbEl.Inclination,iss.OrbEl.RAAN,…iss.OrbEl.ArgPeriapsis,iss.OrbEl。TrueAnomaly,…OrbitPropagator =“two-body-keplerian”,…Name =“空间站开普勒”);iss.KeplerianOrbit.obj。MarkerColor =“红色”;iss.KeplerianOrbit.obj.Orbit。LineColor =“红色”;

相比之下,使用开普勒二体传播的最大误差。

[iss.KeplerianOrbit。EphPosICRFData ~,…iss.KeplerianOrbit。步伐]= (iss.KeplerianOrbit状态。obj CoordinateFrame =“惯性”);iss.KeplerianOrbit。(iss.KeplerianOrbit EphPosICRF =调整时间(时间表。步伐,iss.KeplerianOrbit。EphPosICRFData’, VariableNames =“数据”),…iss.RefEphemeris.pos_eme2000.Time);情节(iss.EphPosICRF。时间,iss.KeplerianOrbit.EphPosICRF。D一个ta/1e3 - iss.RefEphemeris.pos_eme2000.Position); title(“国际空间站的两个身体开普勒位置误差”);包含(“时间”);ylabel (“位置误差(公里)”);传奇([“\ DeltaX_ {icrf}”,“\ DeltaY_ {icrf}”,“\ DeltaZ_ {icrf}”),位置=“西北”);

Two-Body-Keplerian传播之间的最大误差和参考数据超过450公里,超过24小时方法相比,0.25公里的高精度数值模拟传播。

情节(iss.EphPosICRF。时间,vecnorm (iss.KeplerianOrbit.EphPosICRF。数据/ 1 e3 - iss.RefEphemeris.pos_eme2000。位置2 2));标题(“国际空间站的两个身体开普勒位置误差大小”);包含(“时间”);ylabel (“位置误差(公里)”);