J4解析预报器
2026/1/21大约 3 分钟
J4解析预报器
概述
J4解析预报器考虑地球扁率J2和J4项的解析预报器，提供比J2预报器更高精度的轨道预报。该方法在保持较高计算效率的同时，能够更准确地模拟地球非球形引力场的摄动效应。
核心概念
J4摄动项
J4项是地球引力场球谐展开中的高阶项，J4摄动与J2摄动共同作用，影响：
- 轨道面进动：更精确的轨道平面进动计算
- 近地点进动：更准确的近地点旋转速率
预报原理
J4解析预报器扩展了J2预报器的模型：
- 计算J2和J4摄动引起的轨道根数长期变化率
- 对轨道根数进行线性预报
- 计算对应的位置和速度向量
主要接口
类接口
class J4Analytical: public J2J4Analytical构造函数：
J4Analytical(const ModOrbElem& modOrbElem, const TimePoint& epoch, double gm, double j2, double j4, double re);
J4Analytical(const Vector3d& pos, const Vector3d& vel, const TimePoint& epoch, double gm, double j2, double j4, double re);函数接口
AST_CORE_API errc_t aJ4AnalyticalProp(double duration, double gm, double j2, double j4, double re, Vector3d& r, Vector3d& v);
AST_CORE_API errc_t aJ4AnalyticalProp(double duration, double gm, double j2, double j4, double re, ModOrbElem& modOrbElem);参数说明：
duration：传播时间间隔（秒）gm：引力常数（m³/s²）j2：J2项系数（无量纲）j4：J4项系数（无量纲）re：天体半径（米）r：输入/输出位置向量（米）v：输入/输出速度向量（米/秒）modOrbElem：输入/输出改进轨道根数
用法示例
1. 位置速度传播
#include "AstCore/J4Analytical.hpp"
#include "AstCore/TimePoint.hpp"
#include "AstMath/Vector.hpp"
#include "AstUtil/Literals.hpp"
#include "AstUtil/Constants.h"
#include <iostream>
#include <iomanip>
AST_USING_NAMESPACE
using namespace _AST literals;
int main()
{
// 地球参数
const double gm_earth = kEarthGrav;
const double j2_earth = 1.08262668e-3; // 地球J2系数
const double j4_earth = -1.620e-6; // 地球J4系数
const double re_earth = kEarthRadius; // 地球半径
// 初始轨道状态 (近地轨道)
Vector3d position{7000_km, 0.0, 0.0};
Vector3d velocity{0.0, 7.5_km_s, 1.0_km_s};
std::cout << std::fixed << std::setprecision(3);
std::cout << "J4解析传播器 - 位置速度传播示例" << std::endl;
std::cout << "==================================" << std::endl;
std::cout << "初始轨道状态:" << std::endl;
std::cout << "位置: (" << position[0]/1000.0 << ", "
<< position[1]/1000.0 << ", "
<< position[2]/1000.0 << ") km" << std::endl;
std::cout << "速度: (" << velocity[0]/1000.0 << ", "
<< velocity[1]/1000.0 << ", "
<< velocity[2]/1000.0 << ") km/s" << std::endl;
std::cout << std::endl;
// 传播时间间隔 (3天)
double duration = 3 * 24 * 3600.0; // 秒
// 执行J4解析传播
errc_t result = aJ4AnalyticalProp(duration, gm_earth, j2_earth, j4_earth, re_earth, position, velocity);
if (result == eNoError) {
std::cout << "传播后轨道状态 (" << duration/(24*3600.0) << " 天后):" << std::endl;
std::cout << "位置: (" << position[0]/1000.0 << ", "
<< position[1]/1000.0 << ", "
<< position[2]/1000.0 << ") km" << std::endl;
std::cout << "速度: (" << velocity[0]/1000.0 << ", "
<< velocity[1]/1000.0 << ", "
<< velocity[2]/1000.0 << ") km/s" << std::endl;
// 计算轨道倾角变化
double initial_inclination = std::atan2(1.0, 7.5) * kRadToDeg; // 近似计算
double final_inclination = std::atan2(velocity[2], std::sqrt(velocity[0]*velocity[0] + velocity[1]*velocity[1])) * kRadToDeg;
std::cout << std::endl;
std::cout << "轨道参数:" << std::endl;
std::cout << "初始倾角: " << initial_inclination << " 度" << std::endl;
std::cout << "最终倾角: " << final_inclination << " 度" << std::endl;
std::cout << "倾角变化: " << (final_inclination - initial_inclination) << " 度" << std::endl;
} else {
std::cout << "传播失败，错误码: " << result << std::endl;
}
return 0;
}2. 轨道根数传播
#include "AstCore/J4Analytical.hpp"
#include "AstCore/TimePoint.hpp"
#include "AstCore/OrbitElement.hpp"
#include "AstMath/Vector.hpp"
#include "AstUtil/Literals.hpp"
#include "AstUtil/Constants.h"
#include <iostream>
#include <iomanip>
#include <cmath>
AST_USING_NAMESPACE
using namespace _AST literals;
int main()
{
// 地球参数
const double gm_earth = kEarthGrav;
const double j2_earth = 1.08262668e-3;
const double j4_earth = -1.620e-6;
const double re_earth = kEarthRadius;
// 创建改进轨道根数 (大椭圆轨道示例)
ModOrbElem orb_elem;
orb_elem.rp_ = 26500_km; // 近拱点半径
orb_elem.e_ = 0.73; // 偏心率
orb_elem.i_ = 28.5 * kDegToRad; // 倾角
orb_elem.raan_ = 0.0;
orb_elem.argper_ = 0.0;
orb_elem.trueA_ = 0.0;
std::cout << std::fixed << std::setprecision(6);
std::cout << "J4解析传播器 - 轨道根数传播示例" << std::endl;
std::cout << "==================================" << std::endl;
std::cout << "初始轨道根数 (大椭圆轨道):" << std::endl;
std::cout << "半长轴: " << orb_elem.getA() / 1000.0 << " km" << std::endl;
std::cout << "偏心率: " << orb_elem.e_ << std::endl;
std::cout << "倾角: " << orb_elem.i_ * kRadToDeg << " 度" << std::endl;
std::cout << "升交点赤经: " << orb_elem.raan_ * kRadToDeg << " 度" << std::endl;
std::cout << "近地点幅角: " << orb_elem.argper_ * kRadToDeg << " 度" << std::endl;
std::cout << std::endl;
// 传播多个时间点，观察J4摄动效应
double duration_day = 7.0; // 7天
double duration_sec = duration_day * 24 * 3600.0;
std::cout << "轨道根数随时间变化 (" << duration_day << " 天):" << std::endl;
std::cout << "时间(天)\t升交点赤经(度)\t近地点幅角(度)\t偏心率" << std::endl;
std::cout << "----------------------------------------------------------------" << std::endl;
for (int i = 0; i <= 7; i++) {
double time = i * duration_sec / 7.0;
// 复制初始轨道根数
ModOrbElem current_orb = orb_elem;
// 执行J4解析传播
errc_t result = aJ4AnalyticalProp(time, gm_earth, j2_earth, j4_earth, re_earth, current_orb);
if (result == eNoError) {
std::cout << std::setw(6) << time/(24*3600.0) << "\t"
<< std::setw(12) << current_orb.raan_ * kRadToDeg << "\t"
<< std::setw(12) << current_orb.argper_ * kRadToDeg << "\t"
<< std::setw(10) << current_orb.e_ << std::endl;
} else {
std::cout << "传播失败，错误码: " << result << std::endl;
break;
}
}
// 计算J4对近地点进动的修正
double semi_major_axis = orb_elem.getA();
double eccentricity = orb_elem.e_;
double inclination = orb_elem.i_;
// J2引起的近地点进动率
double n = std::sqrt(gm_earth / (semi_major_axis * semi_major_axis * semi_major_axis));
double p = semi_major_axis * (1 - eccentricity * eccentricity);
double aop_rate_j2 = 1.5 * j2_earth * n * std::pow(re_earth / p, 2) * (2 - 2.5 * std::sin(inclination) * std::sin(inclination));
// J4引起的近地点进动率修正
double aop_rate_j4 = 0.0; // 简化计算，实际需要复杂公式
std::cout << std::endl;
std::cout << "摄动理论值:" << std::endl;
std::cout << "J2近地点进动率: " << aop_rate_j2 * kRadToDeg * 86400.0 << " 度/天" << std::endl;
std::cout << "J4提供更高精度修正" << std::endl;
return 0;
}3. 类接口使用
#include "AstCore/J4Analytical.hpp"
#include "AstCore/TimePoint.hpp"
#include "AstCore/OrbitElement.hpp"
#include "AstMath/Vector.hpp"
#include "AstUtil/Literals.hpp"
#include "AstUtil/Constants.h"
#include <iostream>
#include <iomanip>
AST_USING_NAMESPACE
using namespace _AST literals;
int main()
{
std::cout << std::fixed << std::setprecision(6);
std::cout << "J4解析传播器 - 类接口使用示例" << std::endl;
std::cout << "================================" << std::endl;
// 地球参数
const double gm_earth = 3.986004418e14; // 地球引力常数 (m³/s²)
const double j2_earth = 1.08262668e-3; // 地球J2系数
const double j4_earth = -1.620e-6; // 地球J4系数
const double re_earth = 6378137.0; // 地球半径 (m)
// 示例1: 使用轨道根数构造J4传播器
std::cout << "示例1: 使用轨道根数构造J4传播器" << std::endl;
std::cout << "--------------------------------" << std::endl;
TimePoint epoch1 = TimePoint::FromUTC(2026, 1, 1, 0, 0, 0);
// 创建改进轨道根数 (太阳同步轨道)
ModOrbElem modOrbElem1;
modOrbElem1.rp_ = 7178.137_km; // 近拱点半径 (km)
modOrbElem1.e_ = 0.001; // 偏心率
modOrbElem1.i_ = 98.0_deg; // 轨道倾角
modOrbElem1.argper_ = 0.0_deg; // 近拱点角
modOrbElem1.raan_ = 0.0_deg; // 升交点赤经
modOrbElem1.trueA_ = 0.0_deg; // 真近点角
// 使用轨道根数构造J4传播器
J4Analytical j4_prop1(modOrbElem1, epoch1, gm_earth, j2_earth, j4_earth, re_earth);
std::cout << "轨道根数构造成功" << std::endl;
std::cout << "半长轴: " << modOrbElem1.getA()/1000.0 << " km" << std::endl;
std::cout << "偏心率: " << modOrbElem1.e_ << std::endl;
std::cout << "轨道倾角: " << modOrbElem1.i_/kDegToRad << " 度" << std::endl;
std::cout << std::endl;
// 示例2: 使用位置速度构造J4传播器
std::cout << "示例2: 使用位置速度构造J4传播器" << std::endl;
std::cout << "--------------------------------" << std::endl;
TimePoint epoch2 = TimePoint::FromUTC(2026, 1, 1, 0, 0, 0);
// 初始位置速度 (近圆轨道)
Vector3d position{7000_km, 0.0, 0.0};
Vector3d velocity{0.0, 7.5_km_s, 0.0};
// 使用位置速度构造J4传播器
J4Analytical j4_prop2(position, velocity, epoch2, gm_earth, j2_earth, j4_earth, re_earth);
std::cout << "位置速度构造成功" << std::endl;
std::cout << "位置: (" << position[0]/1000.0 << ", "
<< position[1]/1000.0 << ", "
<< position[2]/1000.0 << ") km" << std::endl;
std::cout << "速度: (" << velocity[0]/1000.0 << ", "
<< velocity[1]/1000.0 << ", "
<< velocity[2]/1000.0 << ") km/s" << std::endl;
std::cout << std::endl;
// 示例3: 完整传播流程演示
std::cout << "示例3: 完整传播流程演示" << std::endl;
std::cout << "------------------------" << std::endl;
// 创建时间点
TimePoint start_time = TimePoint::FromUTC(2026, 1, 1, 0, 0, 0);
TimePoint end_time = start_time + 24 * 3600.0; // 1天后
// 初始轨道状态 (大椭圆轨道)
Vector3d pos{26560_km, 0.0, 0.0}; // 约20000km高度
Vector3d vel{0.0, 3.87_km_s, 0.0}; // 约12小时周期
std::cout << "传播前状态:" << std::endl;
std::cout << "起始时间: 2026-01-01 00:00:00 UTC" << std::endl;
std::cout << "位置: (" << pos[0]/1000.0 << ", "
<< pos[1]/1000.0 << ", "
<< pos[2]/1000.0 << ") km" << std::endl;
std::cout << "速度: (" << vel[0]/1000.0 << ", "
<< vel[1]/1000.0 << ", "
<< vel[2]/1000.0 << ") km/s" << std::endl;
// 使用函数接口进行传播
errc_t result = aJ4AnalyticalProp(24 * 3600.0, gm_earth, j2_earth, j4_earth, re_earth, pos, vel);
if (result == eNoError) {
std::cout << std::endl;
std::cout << "传播后状态:" << std::endl;
std::cout << "结束时间: 2026-01-02 00:00:00 UTC" << std::endl;
std::cout << "位置: (" << pos[0]/1000.0 << ", "
<< pos[1]/1000.0 << ", "
<< pos[2]/1000.0 << ") km" << std::endl;
std::cout << "速度: (" << vel[0]/1000.0 << ", "
<< vel[1]/1000.0 << ", "
<< vel[2]/1000.0 << ") km/s" << std::endl;
// 计算轨道参数变化
double initial_radius = 26560.0; // km
double final_radius = pos.norm() / 1000.0;
double radius_change = final_radius - initial_radius;
std::cout << std::endl;
std::cout << "轨道参数变化:" << std::endl;
std::cout << "初始轨道半径: " << initial_radius << " km" << std::endl;
std::cout << "最终轨道半径: " << final_radius << " km" << std::endl;
std::cout << "轨道半径变化: " << radius_change << " km" << std::endl;
std::cout << "注: 变化主要由J2和J4摄动引起" << std::endl;
} else {
std::cout << "传播失败，错误码: " << result << std::endl;
}
std::cout << std::endl;
std::cout << "J4解析传播器类接口使用总结:" << std::endl;
std::cout << "1. 支持使用轨道根数或位置速度构造传播器" << std::endl;
std::cout << "2. 提供函数接口和类接口两种使用方式" << std::endl;
std::cout << "3. 考虑J2和J4摄动项，精度高于J2传播器" << std::endl;
std::cout << "4. 适用于需要较高精度的近地轨道传播" << std::endl;
return 0;
}精度与限制
精度特点
- 计算精度：考虑J2和J4项长期摄动，精度优于J2传播器
- 适用范围：适用于需要较高精度地球扁率效应的近地轨道
- 计算效率：较高，适合需要较高效率的应用
使用限制
- 仅考虑J2和J4项摄动，忽略更高阶项和其他摄动力
- 对于高精度轨道预报需求，建议使用HPOP等数值方法
- 长时间传播时，由于模型简化，误差会累积
地球J2/J4参数
对于地球轨道计算，推荐使用以下参数：
- J2系数：1.08262668e-3
- J4系数：-1.620e-6
- 地球半径：6378137.0米
- 引力常数：3.986004418e14 m³/s²
精度对比
与J2传播器相比，J4传播器在以下方面具有优势：
- 近地点进动：更准确的近地点旋转速率计算
- 长期预报：长时间轨道预报精度更高
- 高倾角轨道：对极轨道和临界倾角轨道的模拟更精确
注意事项
- 输入的位置和速度向量应在同一惯性坐标系下
- J2/J4系数和天体半径应与中心天体匹配
- 传播时间间隔不宜过长，避免数值误差累积
- 对于科学任务和精密轨道确定，建议验证模型精度
API参考
///
/// @file J4Analytical.hpp
/// @brief ~
/// @details ~
/// @author axel
/// @date 2026-01-05
/// @copyright 版权所有 (C) 2026-present, ast项目.
///
/// ast项目（https://github.com/space-ast/ast）
/// 本项目基于 Apache 2.0 开源许可证分发。
/// 您可在遵守许可证条款的前提下使用、修改和分发本软件。
/// 许可证全文请见：
///
/// http://www.apache.org/licenses/LICENSE-2.0
///
/// 重要须知：
/// 软件按"现有状态"提供，无任何明示或暗示的担保条件。
/// 除非法律要求或书面同意，作者与贡献者不承担任何责任。
/// 使用本软件所产生的风险，需由您自行承担。
#pragma once
#include "AstGlobal.h"
#include "J2J4Analytical.hpp"
AST_NAMESPACE_BEGIN
/*!
@addtogroup Propagator
@{
*/
class ModOrbElem;
class AST_CORE_API J4Analytical: public J2J4Analytical
{
public:
/// @brief 构造函数
/// @details ~
/// @param modOrbElem 改进轨道根数
/// @param epoch 历元
/// @param gm 引力常数
/// @param j2 J2项
/// @param j4 J4项
/// @param re 天体半径
J4Analytical(const ModOrbElem& modOrbElem, const TimePoint& epoch, double gm, double j2, double j4, double re);
/// @brief 构造函数
/// @details ~
/// @param pos 位置向量
/// @param vel 速度向量
/// @param epoch 历元
/// @param gm 引力常数
/// @param j2 J2项
/// @param j4 J4项
/// @param re 天体半径
J4Analytical(const Vector3d& pos, const Vector3d& vel, const TimePoint& epoch, double gm, double j2, double j4, double re);
private:
// double j2_;
// double j4_;
// double re_;
};
/// @brief J4轨道预报
/// @details
/// @param duration 时间间隔
/// @param gm 引力常数
/// @param j2 J2项
/// @param j4 J4项
/// @param re 天体半径
/// @param r 位置向量
/// @param v 速度向量
/// @return 错误码
AST_CORE_API errc_t aJ4AnalyticalProp(double duration, double gm, double j2, double j4, double re, Vector3d& r, Vector3d& v);
/// @brief J4轨道预报
/// @details
/// @param duration 时间间隔
/// @param gm 引力常数
/// @param j2 J2项
/// @param j4 J4项
/// @param re 天体半径
/// @param modOrbElem 改进轨道根数
/// @return 错误码
AST_CORE_API errc_t aJ4AnalyticalProp(double duration, double gm, double j2, double j4, double re, ModOrbElem& modOrbElem);
/*! @} */
AST_NAMESPACE_END