4. 地球非球形引力场建模

在二体问题中，地球被简化为质点或均匀球体，引力加速度沿地心矢径方向，轨道根数全部为常数。但真实地球质量分布并不均匀：赤道隆起、两极略扁，还有山脉、海沟、大陆与海洋的质量差异。这些偏离球对称的部分统称为地球非球形引力摄动，是低地球轨道（LEO）航天器最主要的保守摄动源，量级可达二体引力的千分之一。要定量描述它，必须把地球引力位展开成球谐级数。

4.1 引力位的球谐展开

设地心固连坐标系中航天器位置为 \((r,\lambda,\phi)\)，其中 \(r\) 为地心距，\(\lambda\) 为地心经度，\(\phi\) 为地心纬度。地球引力位 \(U\) 可写成：

U(r,λ,φ) = (μ/r) · Σ_{n=0}^{∞} Σ_{m=0}^{n} (R_e/r)^n · P_n^m(sinφ)
            · [ C_nm·cos(mλ) + S_nm·sin(mλ) ]

式中 \(\mu=GM\) 为地球引力常数，\(R_e\) 为地球赤道平均半径，\(P_n^m\) 为缔合勒让德函数，\(C_{nm}\)、\(S_{nm}\) 为无量纲球谐系数，由地球质量分布决定。\(n\) 称为阶，\(m\) 称为次。\(n=0\) 项即 \(\mu/r\)，对应二体引力；\(n\ge 2\) 各项才是非球形摄动。按 \(m\) 的取值，系数分为三类：

• 带谐项（zonal）：\(m=0\)，系数记作 \(J_n=-C_{n0}\)。与经度无关，只随纬度变化，把地球分成南北对称的“带”，反映旋转对称的扁率与纬向质量分布。
• 田谐项（tesseral）：\(m\ne 0\) 且 \(m\ne n\)，系数 \(C_{nm}\)、\(S_{nm}\) 同时随经度和纬度变化，把球面分割成“田”字形网格，反映非轴对称质量分布。
• 扇谐项（sectorial）：\(m=n\)，是田谐项的特例，沿经度呈扇形分布。

工程上最常用的是 \(J_2\)、\(J_3\)、\(J_4\) 三个带谐项，其典型量级为 \(J_2\approx1.0826\times10^{-3}\)、\(J_3\approx-2.53\times10^{-6}\)、\(J_4\approx-1.61\times10^{-6}\)。可见 \(J_2\) 比 \(J_3\)、\(J_4\) 大约三个数量级，因此 \(J_2\) 是首要摄动，\(J_3\)、\(J_4\) 属于高阶修正。

4.2 J₂ 项的摄动公式与物理意义

\(J_2\) 项对应地球的动力学扁率，是赤道隆起造成的。其摄动函数主项为：

R_J2 = (μ J₂ R_e²)/(2 r³) · (3 sin²φ - 1)

将 \(R_{J2}\) 代入拉格朗日行星方程，对轨道根数取平均后，得到著名的 \(J_2\) 长期摄动公式（以近圆轨道近似）：

dΩ/dt = -(3/2) J₂ (R_e/p)² n cos i
dω/dt =  (3/4) J₂ (R_e/p)² n (5 cos²i - 1)
dM/dt  =  (3/4) J₂ (R_e/p)² n √(1-e²) (3 cos²i - 1)

其中 \(p=a(1-e^2)\) 为半通径，\(n=\sqrt{\mu/a^3}\) 为平均角速度。物理意义十分明确：

• 升交点赤经 Ω 的进动：\(d\Omega/dt\propto-\cos i\)。当 \(i<90°\)（顺行轨道）时 \(\Omega\) 持续西退；当 \(i>90°\)（逆行轨道）时 \(\Omega\) 东进。这正是太阳同步轨道设计的核心——选择合适的高度与倾角，使 \(\Omega\) 每年进动约 \(0.9856°\)，与太阳视运动同步。
• 近地点幅角 ω 的进动：\(d\omega/dt\propto(5\cos^2 i-1)\)。当 \(i=63.43°\) 时 \(\cos^2 i=1/5\)，\(\omega\) 长期项为零，称为临界倾角，Molniya 大椭圆轨道即利用此特性冻结近地点。
• 平近点角 M 的漂移：表现为轨道周期的微小修正，影响星下点重复周期。

以 \(a=7000\ \text{km}\)、\(e=0.001\)、\(i=98°\) 的太阳同步轨道为例，代入 \(J_2=1.0826\times10^{-3}\)、\(R_e=6378.137\ \text{km}\)，可估算出 \(\dot\Omega\approx0.9856°/\text{天}\)，与太阳同步条件吻合。这说明 \(J_2\) 不是“误差”，而是可以被主动利用的轨道设计工具。

4.3 J₃、J₄ 项的作用

\(J_3\) 项反映地球南北半球质量分布的不对称（南极略“梨形”外凸），其摄动函数含 \(\sin\phi\) 的奇次项，因此对轨道根数主要产生长周期项而非长期项。对 \(e\) 和 \(\omega\) 的影响表现为：

Δe_J3 ∝ J₃ (R_e/p)³ sin i · cos ω
Δω_J3 ∝ J₃ (R_e/p)³ sin i · sin ω

可见 \(J_3\) 使偏心率矢量在近地点方向附近产生周期性摆动，对近地点高度有几十米到几百米量级的影响，在精密定轨和再入预报中不可忽略。

\(J_4\) 项是 \(J_2\) 的高阶修正，量级约为 \(J_2\) 的千分之一，主要贡献是对 \(\Omega\)、\(\omega\) 长期漂移率的二阶修正，以及对倾角 \(i\) 的周期项。在太阳同步轨道设计中，若只取 \(J_2\) 会使 \(\dot\Omega\) 偏差约 \(0.1°/\text{天}\)，加入 \(J_4\) 后可将设计精度提升一个量级。

4.4 从系数到加速度：实现要点

在数值积分中，通常不直接对位函数求导，而是采用递推算法计算加速度分量。以 \(J_2\) 为例，在地心惯性系中加速度分量为：

a_J2x = -(3/2) J₂ μ R_e² x / r⁵ · (1 - 5z²/r²)
a_J2y = -(3/2) J₂ μ R_e² y / r⁵ · (1 - 5z²/r²)
a_J2z = -(3/2) J₂ μ R_e² z / r⁵ · (3 - 5z²/r²)

其中 \(r=\sqrt{x^2+y^2+z^2}\)。实现时需注意：坐标系必须与球谐系数定义一致（通常为地心固连系），若在惯性系积分，需先做地球自转旋转；\(J_3\)、\(J_4\) 的加速度表达式更复杂，建议使用标准递推公式（如 Cunningham 递推或 Pines 递推）逐阶计算，避免手工推导出错。

验证建模正确性的一个实用方法是：令所有 \(C_{nm}\)、\(S_{nm}\) 为零，只保留 \(n=0\) 项，积分结果应严格退化为二体轨道；再逐步加入 \(J_2\)、\(J_3\)、\(J_4\)，观察 \(\Omega\)、\(\omega\) 的漂移率是否与解析公式一致。若偏差超过 1%，通常说明坐标系定义或系数符号出现了错误。