Orekit 航天动力学库调用方法

Orekit 是一个基于 Java 的优秀航天动力学库，世界各大航天机构及一些科研院所皆有使用。本文将介绍其调用方法，并着重阐述其轨道预报方法。

Orekit 是一个基于 Java 语言开发的航天动力学库，采用对商业友好的 Apache 开源许可协议，第一个公开发行版本始于 2003 年。十多年来，Orekit 一直专注于航天动力学底层算法的实现，包含丰富的航天动力学元素：轨道、时间、参考架、姿态和事件等，以及大量处理这些元素的算法：元素转换、航天器状态预报、卫星姿态指向和事件响应等。

Orekit 在世界各大航天机构及一些科研院所中皆有使用，尤其在欧洲十分流行，诸如 CNES，ESA，NASA，SSTL，SSC，ISRO，Eumetsat，Airbus，Thales，Telespazio，Nexeya，U.S. Naval Research Laboratory，Texas University at Austin 等大型航天机构和科研院所都在使用，通常它作为应用软件的底层代码库，为上层实现提供航天动力学算法支持。

Orekit 项目首页

图 1: Orekit 项目首页

1 获取 Orekit 库

Orekit 官方网站 www.orekit.org 提供封装好的 .jar 库文件以及源码，目前最新版为 v9.3（2019 年 2 月 15 日更新），.jar 库文件可直接使用，源码可编译后使用，此处介绍 .jar 库文件的调用。

2 调用 Orekit 库

2.1 添加 Orekit 库

将 orekit-9.2.jar 添加到项目参考库即可调用 Orekit 中的所有类库，由于 Orekit 采用 Hipparchus 库进行数值运算，所以还需下载并添加该库，方法跟添加 Orekit 库一样，下载地址为 www.hipparchus.org。

添加好的 Orekit 与 Hipparchus 库列表

图 2: 添加好的 Orekit 与 Hipparchus 库列表

2.2 添加动力学建模数据

动力学建模数据包含了 Orekit 建立动力学模型、进行时间和参考架转换以及星历计算所需要的重要参数，必须在调用 Orekit 下属类库之前加载，配置代码如下所示。代码中的文件路径根据具体情况修改。

// configure Orekit
File home       = new File("D:/CodeProjects/eclipse_java/Orekit");
File orekitData = new File(home, "orekit-data");
if (!orekitData.exists()) {
	System.err.format(Locale.US, "Failed to find %s folder%n",
			orekitData.getAbsolutePath());
	System.err.format(Locale.US, "You need to download %s from the %s page and unzip it in %s for this tutorial to work%n",
			"orekit-data.zip", "https://www.orekit.org/forge/projects/orekit/files",
			home.getAbsolutePath());
	System.exit(1);
	}
DataProvidersManager manager = DataProvidersManager.getInstance();
manager.addProvider(new DirectoryCrawler(orekitData));

3 轨道预报

作为航天动力学库，Orekit 最基本的功能性算法就是对航天器状态矢量进行预报，该功能由 org.orekit.propagation 包提供。Orekit 有三种预报方式，一是解析法预报，由 org.orekit.propagation.analytical 实现；二是数值法预报，由 org.orekit.propagation.numerical 和 org.orekit.propagation.integration 实现；三是半解析法预报，由 org.orekit.propagation.semianalytical.dsst 实现。

解析法预报速度快，但精度随时间下降快，且需要平根或拟平根作为预报输入，即需要事先定轨；数值法通过求解轨道动力学方程实现，模型精确的情况下，预报精度高，但由于必须递推中间过程量，故需消耗较多的计算资源；半解析法速度和精度都介于解析法和数值法之间，一般用于卫星轨道寿命与离轨分析，本文不作介绍。

3.1 轨道描述

Orekit 中的轨道描述由 org.orekit.orbits 包提供，该包是构建航天动力学相关程序的基础。从 v3.0 开始，Orekit 就支持多种轨道描述，包括：经典开普勒轨道 KeplerianOrbit、圆轨道 CircularOrbit、赤道轨道 EquinoctialOrbit 以及三维状态矢量 CartesianOrbit。

经典开普勒轨道描述：KeplerianOrbit
1). $a$ ：轨道半场轴（ $m$ ）
2). $e$ ：轨道偏心率（ $NAN$ ）
3). $i$ ：轨道倾角（ $rad$ ）
4). $\omega$ ：近地点幅角（ $rad$ ）
5). $\Omega$ ：升交点赤经（ $rad$ ）
6). $\nu$ ， $E$ or $M$ ：真近角、偏近角或平近角（ $rad$ ）
圆轨道描述：CircularOrbit
1). $a$ ：轨道半场轴（ $m$ ）
2). $e_x$ ：偏心率矢量 X 分量（ $NAN$ ）， $e_x = e\times cos(\omega)$
3). $e_y$ ：偏心率矢量 Y 分量（ $NAN$ ）， $e_y = e\times sin(\omega)$
4). $i$ ：轨道倾角（ $rad$ ）
5). $\Omega$ ：升交点赤经（ $rad$ ）
6). $u_\nu$ ， $u_E$ or $u_M$ ：真纬度幅角（ $\omega + \nu$ ）、偏纬度幅角（ $\omega + E$ ）和平纬度幅角（ $\omega + M$ ）（ $rad$ ）
赤道轨道：EquinoctialOrbit
1). $a$ ：轨道半场轴（ $m$ ）
2). $e_x$ ：偏心率矢量 X 分量（ $NAN$ ）， $e_x = e \times cos(\omega)$
3). $e_y$ ：偏心率矢量 Y 分量（ $NAN$ ）， $e_y = e \times sin(\omega)$
4). $i_x$ ：倾角矢量 X 分量（ $NAN$ ）， $h_x = tan(i/2) \times cos(\Omega)$
5). $i_y$ ：倾角矢量 Y 分量（ $NAN$ ）， $h_y = tan(i/2) \times sin(\Omega)$
5). $\Omega$ ：升交点赤经（ $rad$ ）
6). $l_\nu$ ， $l_E$ or $l_M$ ：真纬度幅角（ $\omega +\nu + \Omega$ ）、偏纬度幅角（ $\omega + E + \Omega$ ）和平纬度幅角（ $\omega + M + \Omega$ ）（ $rad$ ）
三维状态矢量：CartesianOrbit
1). $X, Y, Z$ ：位置矢量（ $m$ ）
2). $V_x, V_y, V_z$ ：速度矢量（ $m/s$ ）

3.2 轨道预报流程

3.2.1 输入参数

四种轨道描述方法定义的轨道根数都可以作为轨道预报的输入，采用 $Orbit$ 关键字进行定义。开普勒轨道根数定义方法如下：

1 2	double mu = 3.986004415e+14； Orbit kepOrbit = new KeplerianOrbit(a, e, i, omega, raan, M, PositionAngle.MEAN, Frame, Date, mu);

代码中， $M$ 对应 PositionAngle.MEAN，即表示平近角，真近角和偏近角分别用 PositionAngle.TURE 和 PositionAngle.ECCENTRIC 表示； $mu$ 为中心天体引力常数， $Frame$ 为参考架，J2000 坐标系用 FramesFactory.getEME2000() 表示，遵循 IERS 2010 约定的地球固连坐标系，简称地固系，用 FramesFactory.getTIRF(IERSConventions.IERS_2010) 表示； $Date$ 为时间，采用如下代码初始化：

1	AbsoluteDate Date = new AbsoluteDate(yr, mt, day, hr, min, sec, TimeScalesFactory.getUTC());

类似地，圆轨道根数定义方法如下：

1 2	double mu = 3.986004415e+14； Orbit cirOrbit = new CircularOrbit(a, ex, ey, i, raan, u, PositionAngle.MEAN, Frame, Date, mu);

赤道轨道根数定义方法如下：

1 2	double mu = 3.986004415e+14； Orbit equOrbit = new EquinoctialOrbit(a, ex, ey, hx, hy, l, PositionAngle.MEAN, Frame, Date, mu);

三维状态矢量定义方法如下：

1
2
3

double mu           = 3.986004415e+14；
PVCoordinates catPV = new PVCoordinates(new Vector3D, new Vector3D);
Orbit catOrbit      = new CartesianOrbit(catPV, Frame, Date, mu);

代码中， $catPV$ 为状态矢量。根据轨道数据的来源，选择对应的轨道描述方式。

3.2.2 建立轨道预报器

Orekit 提供多种解析法预报器，由 org.orekit.propagation.analytical 包实现，最简单的为基于 二体模型 的解析法预报器：

1	KeplerianPropagator tbProp = new KeplerianPropagator(Orbit initialOrbit);

较为复杂的解析法预报器基于 EcksteinHechler 模型：

double referenceRadius = 6378136.3;
double mu  =  3.986004415e+14；
double c20 = -0.484165143790815e-03；
double c30 =  0.957161207093473e-06；
double c40 =  0.539965866638991e-06；
double c50 =  0.686702913736681e-07；
double c60 = -0.149953927978527e-06
EcksteinHechlerPropagator ehProp = 
        new EcksteinHechlerPropagator(Orbit initialOrbit, referenceRadius, mu, c20, c30, c40, c50, c60);

Orekit 数值预报器由 org.orekit.propagation.numerical 包实现，定义方法如下：

1	NumericalPropagator propagator = new NumericalPropagator(integrator);

参数 integrator 较具有多种选择，通用型高精度积分器可选择 PD 系列中的 DormandPrince853Integrator，定义方法如下：

final double minStep = 0.001;
final double maxstep = 100.0;
final double positionTolerance = 1.0e-13;
final OrbitType propagationType = OrbitType.KEPLERIAN;
final double[][] tolerances = 
        NumericalPropagator.tolerances(positionTolerance, initialOrbit, propagationType);
AdaptiveStepsizeIntegrator integrator = 
        new DormandPrince853Integrator(minStep, maxstep, tolerances[0], tolerances[1]);

此外，数值轨道预报，还需要设置动力学模型，一般只考虑地球引力场：

final NormalizedSphericalHarmonicsProvider provider =
        GravityFieldFactory.getNormalizedProvider(20, 20);
ForceModel holmesFeatherstone = 
        new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);

更多力学模型的设置参考 API 文档中的 NumericalPropagator 类：orekit-9.2-javadoc/org/
orekit/propagation/numerical/NumericalPropagator.html。

3.2.3 执行轨道预报

执行预报语句，获取指定时间的轨道根数 $currentState$ ：

SpacecraftState currentState = kepler.propagate(extrapDate);
java

也可构建循环语句，以定步长预报指定时间段的轨道根数：

```java
double duration = 600.;
AbsoluteDate finalDate = initialDate.shiftedBy(duration);
double stepT = 60.;
int cpt = 1;
     extrapDate.compareTo(finalDate) <= 0;
     extrapDate = extrapDate.shiftedBy(stepT))  {
    SpacecraftState currentState = kepler.propagate(extrapDate);
    System.out.println("step " + cpt++);
    System.out.println(" time : " + currentState.getDate());
    System.out.println(" " + currentState.getOrbit());
}

数值法预报器于解析法稍有不同，是通过 StepHandler 执行数据操作的：

propagator.addForceModel(holmesFeatherstone);
propagator.setInitialState(initialState);
propagator.setMasterMode(60., new TutorialStepHandler());
SpacecraftState finalState = 
        propagator.propagate(initialDate.shiftedBy(600));

4 总结

本文介绍了 Orekit 航天动力学库的调用方法，并简单阐述了其轨道预报方法，可作为相关应用软件开发的入门级参考。