Orekit，一个用 Java 编写的低级空间动力学库，已经获得了 自从它在开源下发布以来，它得到了广泛的认可 2008年获得许可。
maven依赖
<dependency>
<groupId>org.orekit</groupId>
<artifactId>orekit</artifactId>
<version>11.3.3</version>
</dependency>
作为航天动力学库，Orekit 最基本的功能性算法就是对航天器状态矢量进行预报，该功能由 org.orekit.propagation 包提供。Orekit 有三种预报方式，一是解析法预报，由 org.orekit.propagation.analytical 实现；二是数值法预报，由 org.orekit.propagation.numerical 和 org.orekit.propagation.integration 实现；三是半解析法预报，由 org.orekit.propagation.semianalytical.dsst 实现。
解析法预报速度快，但精度随时间下降快，且需要平根或拟平根作为预报输入，即需要事先定轨；数值法通过求解轨道动力学方程实现，模型精确的情况下，预报精度高，但由于必须递推中间过程量，故需消耗较多的计算资源；半解析法速度和精度都介于解析法和数值法之间，一般用于卫星轨道寿命与离轨分析
Orekit 中的轨道描述由 org.orekit.orbits 包提供，该包是构建航天动力学相关程序的基础。从 v3.0 开始，Orekit 就支持多种轨道描述，包括：经典开普勒轨道 KeplerianOrbit、圆轨道 CircularOrbit、赤道轨道 EquinoctialOrbit 以及三维状态矢量 CartesianOrbit。
- 经典开普勒轨道描述：
KeplerianOrbit1). a：轨道半场轴（m） 2). e：轨道偏心率（NAN） 3). i：轨道倾角（rad） 4). ω：近地点幅角（rad） 5). Ω：升交点赤经（rad） 6). ν，E or M：真近角、偏近角或平近角（rad） - 圆轨道描述：
CircularOrbit1). a：轨道半场轴（m） 2). ex：偏心率矢量 X 分量（NAN），ex=e×cos(ω) 3). ey：偏心率矢量 Y 分量（NAN），ey=e×sin(ω) 4). i：轨道倾角（rad） 5). Ω：升交点赤经（rad） 6). uν，uE or uM：真纬度幅角（ω+ν）、偏纬度幅角（ω+E）和平纬度幅角（ω+M）（rad） - 赤道轨道：
EquinoctialOrbit1). a：轨道半场轴（m） 2). ex：偏心率矢量 X 分量（NAN），ex=e×cos(ω) 3). ey：偏心率矢量 Y 分量（NAN），ey=e×sin(ω) 4). ix：倾角矢量 X 分量（NAN），hx=tan(i/2)×cos(Ω) 5). iy：倾角矢量 Y 分量（NAN），hy=tan(i/2)×sin(Ω) 5). Ω：升交点赤经（rad） 6). lν，lE or lM：真纬度幅角（ω+ν+Ω）、偏纬度幅角（ω+E+Ω）和平纬度幅角（ω+M+Ω）（rad） - 三维状态矢量：
CartesianOrbit1). X,Y,Z：位置矢量（m） 2). Vx,Vy,Vz：速度矢量（m/s）
轨道预报流程
1.添加动力学建模数据
动力学建模数据包含了 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));
2.创建轨道
四种轨道描述方法定义的轨道根数都可以作为轨道预报的输入，采用 Orbit 关键字进行定义。开普勒轨道根数定义方法如下：
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 为时间，采用如下代码初始化：
AbsoluteDate Date = new AbsoluteDate(yr, mt, day, hr, min, sec, TimeScalesFactory.getUTC());
类似地，圆轨道根数定义方法如下：
double mu = 3.986004415e+14；
Orbit cirOrbit = new CircularOrbit(a, ex, ey, i, raan, u, PositionAngle.MEAN, Frame, Date, mu);
赤道轨道根数定义方法如下：
double mu = 3.986004415e+14；
Orbit equOrbit = new EquinoctialOrbit(a, ex, ey, hx, hy, l, PositionAngle.MEAN, Frame, Date, mu);
三维状态矢量定义方法如下：
double mu = 3.986004415e+14；
PVCoordinates catPV = new PVCoordinates(new Vector3D, new Vector3D);
Orbit catOrbit = new CartesianOrbit(catPV, Frame, Date, mu);
代码中，catPV 为状态矢量。根据轨道数据的来源，选择对应的轨道描述方式。
3.建立轨道预报器
Orekit 提供多种解析法预报器，由 org.orekit.propagation.analytical 包实现，最简单的为基于 二体模型 的解析法预报器：
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 包实现，定义方法如下：
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。
轨道预报流程那部分内容出自spacefan.github.io/2019/01/25/…
轨道六根数转换成笛卡尔积
代码：
File file = new File("C:/SoftwareInstallationPackage/orekit-data-master/");
DataProvidersManager manager = DataContext.getDefault().getDataProvidersManager();
manager.addProvider(new DirectoryCrawler(file));
double a=7275435.324; //轨道半长轴
double e=0.0019976693; //轨道偏心率
double i=89.0603273252; //轨道倾角
double omega=72.2698073713; //近地点幅角
double raan=26.9856045949; //升交点赤经
double m=174.4969934001; //平近点角
//double mu = 3.986004415e+14;
double mu = 398600441500000.000;
//AbsoluteDate absoluteDate = new AbsoluteDate(new Date(), TimeScalesFactory.getUTC());
LocalDateTime time = LocalDateTime.now(Clock.systemUTC());
String format = time.format(DateTimeFormatter.ofPattern("yyyy-MM-dd'T'HH:mm:ss'Z'"));
AbsoluteDate absoluteDate = new AbsoluteDate(format, TimeScalesFactory.getUTC());
AbsoluteDate end = new AbsoluteDate("2023-08-18T18:04:00", TimeScalesFactory.getUTC());
FactoryManagedFrame frame = FramesFactory.getEME2000(); //J2000坐标系
KeplerianOrbit keplerianOrbit = new KeplerianOrbit(a, e, i, omega, raan, m, PositionAngle.MEAN, frame, absoluteDate, mu);
double[][] tolerances = NumericalPropagator.tolerances(1, keplerianOrbit, OrbitType.KEPLERIAN);
DormandPrince853Integrator integrator = new DormandPrince853Integrator(0.001, 100, tolerances[0], tolerances[1]);
NumericalPropagator propagator = new NumericalPropagator(integrator); //创建数值预报器
final NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(20, 20);
HolmesFeatherstoneAttractionModel holmesFeatherstoneAttractionModel = new HolmesFeatherstoneAttractionModel(FramesFactory.getTIRF(IERSConventions.IERS_2010, true), provider);
propagator.setOrbitType(OrbitType.KEPLERIAN); //设置轨道类型
propagator.addForceModel(holmesFeatherstoneAttractionModel);
propagator.addForceModel(new NewtonianAttraction(keplerianOrbit.getMu()));
propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getEarth()));
propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));
propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getJupiter()));
propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMars()));
propagator.setInitialState(new SpacecraftState(keplerianOrbit));
//EphemerisGenerator generator = propagator.getEphemerisGenerator();
//propagator.propagate(absoluteDate.shiftedBy(Time.daysToSeconds(Time.DAYS_PER_YEAR) * 1))
//BoundedPropagator generatedEphemeris = generator.getGeneratedEphemeris();
/*while (absoluteDate.isBefore(end)) {
TimeStampedPVCoordinates pvc = propagator.propagate(absoluteDate).getPVCoordinates();
Vector3D position = pvc.getPosition();
Vector3D velocity = pvc.getVelocity();
System.out.println("x: "+position.getX()+" y: "+position.getY()+" z:"+position.getZ());
System.out.println("vx: "+velocity.getX()+" vy: "+velocity.getY()+" vz:"+velocity.getZ());
System.out.println("--------------------------------------------------");
}*/
TimeStampedPVCoordinates pvc = propagator.propagate(absoluteDate).getPVCoordinates();
Vector3D position = pvc.getPosition();
Vector3D velocity = pvc.getVelocity();
System.out.println("x: "+position.getX()+" y: "+position.getY()+" z:"+position.getZ());
System.out.println("vx: "+velocity.getX()+" vy: "+velocity.getY()+" vz:"+velocity.getZ());
System.out.println("--------------------------------------------------");
预报某一地面站能观测到某一卫星的时间
代码如下：
File file = new File("C:/SoftwareInstallationPackage/orekit-data-master/");
DataProvidersManager manager = DataContext.getDefault().getDataProvidersManager();
manager.addProvider(new DirectoryCrawler(file));
// Initial state definition : date, orbit
final AbsoluteDate initialDate = new AbsoluteDate(2023, 9, 15, 4, 0, 00.000, TimeScalesFactory.getUTC());
final double mu = 3.986004415e+14; // gravitation coefficient
final Frame inertialFrame = FramesFactory.getEME2000(); // inertial frame for orbit definition
final Vector3D position = new Vector3D(2797914.567, -2288195.171, 6012468.374);
final Vector3D velocity = new Vector3D(-6089.132, 2403.774, 3732.121);
final PVCoordinates pvCoordinates = new PVCoordinates(position, velocity);
final Orbit initialOrbit = new KeplerianOrbit(pvCoordinates, inertialFrame, initialDate, mu);
// Earth and frame
final Frame earthFrame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
final BodyShape earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
Constants.WGS84_EARTH_FLATTENING, earthFrame);
// Station
final double longitude = FastMath.toRadians(75.9797);
final double latitude = FastMath.toRadians(39.4547);
final double altitude = 0.;
final GeodeticPoint station1 = new GeodeticPoint(latitude, longitude, altitude);
final TopocentricFrame sta1Frame = new TopocentricFrame(earth, station1, "喀什");
// Defining rectangular field of view
double halfApertureAlongTrack = FastMath.toRadians(50);
double halfApertureAcrossTrack = FastMath.toRadians(50);
FieldOfView fov = new DoubleDihedraFieldOfView(Vector3D.MINUS_I, // From satellite to body center
Vector3D.PLUS_K, halfApertureAcrossTrack, // Across track direction
Vector3D.PLUS_J, halfApertureAlongTrack, // Along track direction
0); // Angular margin
// Defining attitude provider
AttitudeProvider attitudeProvider = new LofOffset(inertialFrame, LOFType.QSW);
// Defining your propagator and setting up the attitude provider
Propagator propagator = new KeplerianPropagator(initialOrbit);
propagator.setAttitudeProvider(attitudeProvider);
// Event definition
final double maxcheck = 60.0;
final double threshold = 0.001;
final double elevation = FastMath.toRadians(5.0);
final EventDetector sta1Visi = new ElevationDetector(maxcheck, threshold, sta1Frame)
.withConstantElevation(elevation).withHandler((s, detector, increasing) -> {
if (increasing) {
System.out.println("Visibility begins on " + detector.getTopocentricFrame().getName() + " at "
+ s.getDate());
} else {
System.out.println("Visibility ends on " + detector.getTopocentricFrame().getName() + " at "
+ s.getDate());
}
return Action.CONTINUE; // Continue processing the event
});
// Add event to be detected
propagator.addEventDetector(sta1Visi);
// Propagate from the initial date to the first raising or for the fixed
// duration
final SpacecraftState finalState = propagator.propagate(initialDate.shiftedBy(14400.));
System.out.println(" Final state : " + finalState.getDate().durationFrom(initialDate));
说明
final Frame earthFrame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
final BodyShape earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
Constants.WGS84_EARTH_FLATTENING, earthFrame);
上面代码是创建一个地球模型，参考系是地固系，OneAxisEllipsoid的参数一代表赤道半径，参数二代表扁平化，参数三代表参考系。
final double longitude = FastMath.toRadians(75.9797); //经度
final double latitude = FastMath.toRadians(39.4547); //纬度
final double altitude = 0.; //高度
final GeodeticPoint station1 = new GeodeticPoint(latitude, longitude, altitude);
final TopocentricFrame sta1Frame = new TopocentricFrame(earth, station1, "喀什");
上面代码是创建地面站
double halfApertureAlongTrack = FastMath.toRadians(50);
double halfApertureAcrossTrack = FastMath.toRadians(50);
FieldOfView fov = new DoubleDihedraFieldOfView(Vector3D.MINUS_I, // From satellite to body center
Vector3D.PLUS_K, halfApertureAcrossTrack, // Across track direction
Vector3D.PLUS_J, halfApertureAlongTrack, // Along track direction
0); // Angular margin
上面的代码是创建矩形的传感器
// Defining attitude provider
AttitudeProvider attitudeProvider = new LofOffset(inertialFrame, LOFType.QSW);
// Defining your propagator and setting up the attitude provider
Propagator propagator = new KeplerianPropagator(initialOrbit);
propagator.setAttitudeProvider(attitudeProvider);
上面代码是预报器设置姿态信息
// Event definition
final double maxcheck = 60.0;
final double threshold = 0.001;
final double elevation = FastMath.toRadians(5.0);
final EventDetector sta1Visi = new ElevationDetector(maxcheck, threshold, sta1Frame)
.withConstantElevation(elevation).withHandler((s, detector, increasing) -> {
if (increasing) {
System.out.println("Visibility begins on " + detector.getTopocentricFrame().getName() + " at "
+ s.getDate());
} else {
System.out.println("Visibility ends on " + detector.getTopocentricFrame().getName() + " at "
+ s.getDate());
}
return Action.CONTINUE; // Continue processing the event
});
// Add event to be detected
propagator.addEventDetector(sta1Visi);
轨道预报器设置卫星事件。ElevationDetector用来设置事件，其参数一代表最大间隔时间（秒），参数二代表收敛阈值（秒）,参数三代表事件时间搜索的最大迭代次数.
使用HolmesFeatherstoneAttractionModel模型进行轨道预报
代码如下：
File file = new File("C:/SoftwareInstallationPackage/orekit-data-master/");
DataProvidersManager manager = DataContext.getDefault().getDataProvidersManager();
manager.addProvider(new DirectoryCrawler(file));
UTCScale utc = TimeScalesFactory.getUTC();
AbsoluteDate initDate = new AbsoluteDate(2023, 9, 12, 04, 00, 00.000, utc);
double a = 6675900;
double e = 0.01256;
double i = FastMath.toRadians(28.482);
double omega = FastMath.toRadians(0);
double rann = FastMath.toRadians(0);
double mean = FastMath.toRadians(0);
Orbit initOrbit = new KeplerianOrbit(a,e,i,omega,rann, mean, PositionAngle.MEAN,
FramesFactory.getEME2000(), initDate,Constants.IERS2010_EARTH_MU);
double dp = 1.0E-6;
double minStep = 0.00001;
double maxStep = 10000.0;
ODEIntegrator integrator = new DormandPrince853IntegratorBuilder(minStep, maxStep, dp)
.buildIntegrator(initOrbit, OrbitType.KEPLERIAN);
// Initialize a NumericalPropagator
NumericalPropagator propagator = new NumericalPropagator(integrator);
propagator.setMu(Constants.EGM96_EARTH_MU);
SpacecraftState initState = new SpacecraftState(initOrbit, 1000);
propagator.setInitialState(initState);
propagator.setPositionAngleType(PositionAngle.MEAN);
// Create a HolmesFeatherstoneAttractionModel
NormalizedSphericalHarmonicsProvider nshp = GravityFieldFactory.getNormalizedProvider(2, 0);
HolmesFeatherstoneAttractionModel gravityField = new HolmesFeatherstoneAttractionModel(
CelestialBodyFactory.getEarth().getBodyOrientedFrame(), nshp);
// Add the force model
propagator.addForceModel(gravityField);
AbsoluteDate start = new AbsoluteDate(2023, 9, 12, 04, 00, 00.000, utc);
AbsoluteDate end = new AbsoluteDate(2023, 9, 12, 05, 00, 00.000, utc);
// Start forecasting 一分钟执行一次预报
while (start.isBefore(end)){
TimeStampedPVCoordinates pvCoordinates = propagator.getPVCoordinates(start, FramesFactory.getEME2000());
System.out.println(pvCoordinates);
start = start.shiftedBy(60);
}
以上仅是自己的学习笔记
