基于Chebyshev多项式的轻量化业余卫星多普勒频移校正算法
摘要
低轨卫星通信中的多普勒频移校正通常需要运行SGP4轨道预报器，其代码和星历数据占用超过30 KB Flash，对于资源受限的Cortex-M0微控制器（常用业余手台固件）难以承受。本文提出一种轻量级替代方案：利用Chebyshev多项式族拟合归一化径向速度曲线，以过境时长、卫星高度和最大仰角为参数。采用Q1.16定点数格式消除浮点库依赖，算法+系数总计仅需约500字节Flash存储。在10颗FM业余卫星、1420次过境、30天窗口上的留一卫星交叉验证表明，本方法VHF波段平均最大误差397 Hz，UHF波段868 Hz，完全满足窄带FM ±5 kHz频偏容限。进一步的泛化性实验表明：（1）90天内误差无明显退化，一年后归一化曲线形状误差仅增加6 %；（2）在泰安训练的模型直接应用于全球6个纬度地点，UHF误差偏差在13 %以内；（3）Δf_max估计误差一年后增加75 %（主要由轨道衰减导致，需定期更新星历），总误差仍在FM可用范围内。算法仅需用户输入三参数：卫星选择、AOS时间和LOS时间。
1. 引言
1.1 问题背景
业余无线电爱好者使用FM卫星（如SO-50、AO-91、ISS跨段中继）进行语音通信。这些卫星运行在400–1400 km高度的近圆轨道，速度约7.5 km/s，产生的多普勒频移在145 MHz（2m波段）可达±3.5 kHz，在435 MHz（70 cm波段）可达±10.5 kHz。若不进行实时频率校正，接收信号数秒内即漂出接收机通带。
最近开发Baofeng UV-K61固件时发现，在低端业余手台微控制器上运行SGP4轨道预报器，需要存储完整的TLE星历并实现SGP4算法，Flash占用超过30 KB（考虑到soft-float代码占用）。对于总Flash仅64-128 KB的Cortex-M0电台，这迫使开发者做出困难的取舍：要么放弃其他电台功能运行专用SGP4固件，要么使用电脑预先下载过境数据。
1.2 现有方案及不足
当前开源固件主要采用两种策略，以Quansheng UV-K5为例（不考虑高性能DMR手台）：
- 全功能SGP4固件：将所有Flash用于SGP4，放弃其他电台功能。
- 预计算查找表（LUT）：每次过境前从PC通过串口上传校正表。需要电脑且每次过境需重新烧写。
1.3 本文贡献
我们证明：卫星过境期间的多普勒曲线可由一族低阶Chebyshev多项式精确逼近，仅需用户已知的三个参数——卫星身份、AOS时间、LOS时间。本算法：
- 采用 Q1.16 定点数，消除浮点库依赖，仅需 约500字节 Flash（比SGP4节省100倍以上）
- 达到 397 Hz（VHF）/ 868 Hz（UHF） 平均误差（vs. SGP4真值）
- 每次频率更新仅需 30个CPU周期（Cortex-M0纯整数指令）
- 无需每次过境上传数据——用户仅需选择卫星并输入过境时间
- 归一化曲线形状具有时间不变性——一年后形状误差仅退化6 %
- 全球泛化——一组系数适用于所有纬度（UHF偏差<13 %）
2. 数据与方法
2.1 数据来源
TLE星历：2026年8月3日从 Celestrak业余无线电组 获取10颗常用FM业余卫星的最新TLE。
观测站：山东省泰安市（36.18 °N，117.13 °E），最低仰角5 °。
分析卫星：
| 卫星 | NORAD | 下行频率 | 高度 | 30天过境次数 |
|---|
| SO-50 (SaudiSat-1C) | 27607 | 436.795 MHz | 627 km | 163 |
| AO-91 (RadFxSat/Fox-1B) | 43017 | 145.960 MHz | 533 km | 109 |
| AO-85 (Fox-1A) | 40967 | 145.980 MHz | 601 km | 138 |
| AO-95 (Fox-1Cliff) | 43770 | 145.940 MHz | 490 km | 112 |
| LilacSat-2 (CAS-3H) | 40908 | 437.200 MHz | 415 km | 103 |
| ISS (Zarya) | 25544 | 437.800 MHz | 426 km | 208 |
| PO-101 (Diwata-2) | 43678 | 437.500 MHz | 572 km | 150 |
| AO-73 (FUNcube-1) | 39444 | 145.950 MHz | 543 km | 118 |
| NO-44 (PCSAT) | 26931 | 145.825 MHz | 792 km | 157 |
| RS-44 | 44909 | 145.935 MHz | 1350 km | 162 |
总计：1420次过境（VHF 796次 + UHF 624次），时间窗30天（2026-08-03至2026-09-02）。
2.2 真值生成
对每次过境，以SGP4模型通过skyfield Python库（封装标准sgp4包）传播TLE，在150个均匀时间采样点计算真实多普勒频移：
$$\Delta f{\text{true}}(t) = -f{\text{carrier}} \cdot \frac{v_r(t)}{c} \tag{1}$$
其中 $v_r(t)$ 为径向速度（远离观测者时为正），由卫星与观测者在地心坐标系中的相对位置差分得到，$c = 299\;792\;458$ m/s为光速。
2.3 核心洞察：归一化与频率无关性
本方法的关键物理洞察为：多普勒曲线的形状与载波频率无关。以AOS时刻的多普勒频移归一化：
$$\text{norm}(t_{\text{norm}}) = \frac{\Delta f(t)}{\Delta f(\text{AOS})} = \frac{v_r(t)}{v_r(\text{AOS})} \tag{2}$$
载波频率 $f_{\text{carrier}}$ 在归一化中精确抵消。归一化曲线仅取决于卫星轨道几何——具体而言：
- 卫星高度 $h$：决定轨道角速度 $\omega = \sqrt{GM/(R+h)3}$
- 最大仰角 $\theta_{\max}$：决定过境的对称性和TCA附近S曲线的陡峭程度
推论：对归一化径向速度建模的任何误差 $\epsilon$，转换到频率误差为 $\epsilon \cdot f{\text{carrier}}/c$。由于 $f{\text{UHF}}/f_{\text{VHF}} \approx 3$，UHF误差固有地是VHF误差的3倍。
2.4 近似模型
我们评估了六种复杂度递增的模型：
模型A：余弦基线
$$\Delta \hat{f}(t) = \Delta f{\max} \cdot \cos(\pi \cdot t{\text{norm}}), \quad t{\text{norm}} = \frac{t - t{\text{AOS}}}{t{\text{LOS}} - t{\text{AOS}}} \in [0, 1] \tag{3}$$
模型B：单项式多项式（$d$阶）
$$\Delta \hat{f}(t) = \Delta \hat{f}{\max} \cdot \sum{i=0}^{d} c_i \cdot t_{\text{norm}}i \tag{4}$$
系数 $c_i$ 通过最小二乘拟合归一化SGP4曲线集合得到。用Horner方法求值。
模型C：Chebyshev多项式（$d$阶）
$$\Delta \hat{f}(t) = \Delta \hat{f}{\max} \cdot \sum{i=0}^{d} a_i \cdot T_i(x), \quad x = 2 t_{\text{norm}} - 1 \in [-1, 1] \tag{5}$$
其中 $T_i(x) = \cos(i \cdot \arccos(x))$ 是 $[-1, 1]$ 上的第一类Chebyshev多项式。系数 $a_i$ 通过最小二乘拟合得到。求值采用Clenshaw递推：
$$\begin{cases} b{d+2} = b{d+1} = 0 \ b_k = a_k + 2x \cdot b{k+1} - b{k+2}, \quad k = d, d-1, \ldots, 1 \ \text{norm} = a_0 + x \cdot b_1 - b_2 \end{cases} \tag{6}$$
Clenshaw递推在Cortex-M0的float32精度下比同阶Horner方法数值更稳定。
模型D：3-bin仰角分段Chebyshev-7
将过境曲线按最大仰角分为三组（低/中/高，以训练集第33和67百分位数作为分割点），每组独立拟合Chebyshev-7多项式：
$$\Delta \hat{f}(t) = \Delta \hat{f}{\max} \cdot \sum{i=0}^{7} a_i^{(\text{bin})} \cdot T_i(2t{\text{norm}} - 1), \quad \text{bin} = \begin{cases} \text{low} & \theta{\max} \leq \theta{33} \ \text{mid} & \theta{33} < \theta{\max} \leq \theta{67} \ \text{high} & \theta{\max} > \theta{67} \end{cases} \tag{7}$$
Δf_max估计模型
所有多项式模型都需要峰值多普勒 $\Delta \hat{f}_{\max}$ 的估计。我们使用高度感知二次模型：
$$\Delta \hat{f}_{\max} = \alpha_1 T2 + \alpha_2 T + \alpha_3 T h + \alpha_4 h + \alpha_5 \tag{8}$$
其中 $T = t{\text{LOS}} - t{\text{AOS}}$ 为过境时长（秒），$h$ 为卫星高度（km）。该模型捕捉了如下物理事实：过境时长和峰值多普勒均依赖于卫星高度和最大仰角，已知卫星身份（从而知其高度）即可消除"不同卫星相同过境时长对应不同Δf_max"的歧义。
2.5 验证方法
主实验：采用留一卫星交叉验证（LOSO CV）：对数据集中每颗卫星 $s$：
- 将卫星 $s$ 的所有过境作为测试集
- 用其余9颗卫星的数据训练所有模型
- 在卫星 $s$ 的每条过境曲线上评估
这确保了测试卫星的数据从未出现在训练中，提供了对未知卫星泛化能力的真实评估。
泛化性实验：
- 90天泛化：在2026年8月（第1–30天）数据上训练，在9–10月（第31–90天）数据上测试
- 一年泛化：在2026年8月数据上训练，在2027年8–9月（第365–395天）数据上测试，检验多项式系数对TLE老化和轨道衰减的鲁棒性
评估指标：对每条测试过境，计算150个采样点上预测与SGP4真值的最大绝对误差。报告这些逐过境最大误差的均值、中位数、P95和最差值。
3. 实验与结果
3.1 余弦模型为何失败
图1：余弦模型与SGP4对比
图1a–b 展示了SO-50（436.8 MHz）和AO-91（146.0 MHz）高仰角过境的余弦模型与SGP4真值对比。余弦曲线与SGP4形状存在本质错配：它在TCA附近不够陡峭（卫星最近、天空角速度最快时），而在地平线附近过于陡峭（卫星视运动缓慢时）。误差曲线（绿色）呈现典型的"双峰"残差，SO-50最差达到3.3 kHz——占FM频偏的66 %。
图1c 叠加了全部141条归一化多普勒曲线。灰色曲线间的散布代表了任何单一多项式无法捕捉的不可约形状变异。黑色粗线为Chebyshev-7均值。
图1d 展示了多项式阶数平台效应：对VHF和UHF波段，CV误差在7阶后停止下降。误差地板（VHF约500 Hz，UHF约1440 Hz）代表曲线间形状散布——而非多项式灵活性不足。CV曲线与in-sample曲线几乎重合，确认即使到15阶也未发生过度拟合。Chebyshev与Monomial在同等阶数下精度完全一致——两者张成相同的多项式空间，但Chebyshev系数衰减更快，在float32下数值更稳定。
3.2 归一化多普勒曲线的物理本质
图2：物理机制与仰角依赖性
图2a 展示了过境几何：径向速度 $v_r(t)$ 是卫星轨道速度矢量在观测者视线方向上的投影。在AOS和LOS时刻（低仰角），速度几乎垂直于视线，多普勒变化缓慢。在TCA附近（高仰角），速度与视线对齐，多普勒急剧变化。
图2b 验证了3倍缩放定律：将145 MHz的多普勒曲线乘以3，精确匹配同一过境几何下435 MHz的曲线。这确认了 $\text{norm}(t) = v_r(t)/v_r(\text{AOS})$ 的频率无关性。
图2c 揭示了一项关键发现：对同一颗卫星（SO-50），归一化曲线形状随最大仰角剧烈变化。9 °仰角过境高度非对称（TCA远偏离中点），而87 °过境几乎对称。形状差异达满量程的47 %。
图2d 按最大仰角着色所有归一化曲线，直观展示了仰角是曲线形状变异的主要驱动因素。高仰角过境（黄色）紧密聚集在对称S曲线周围，低仰角过境（紫色）则散布广泛。
3.3 模型对比
图3：模型对比
图3a 展示了过境时长到峰值多普勒的映射关系。高度感知二次模型在70 cm波段达到 $R2 = 0.982$（RMS残差397 Hz），而不含高度信息的波段级二次模型仅 $R2 = 0.889$。各个卫星的Δf_max估计残差见表1。
表1：Δf_max估计残差（70 cm波段，高度感知模型）
| 卫星 | 均值残差 | 标准差 | 最大残差 |
|---|
| ISS | −318 Hz | 185 Hz | 635 Hz |
| LilacSat-2 | +482 Hz | 133 Hz | 687 Hz |
| PO-101 | +340 Hz | 260 Hz | 782 Hz |
| SO-50 | −212 Hz | 246 Hz | 614 Hz |
图3b 对比了各模型的形状误差（CV，假设完美Δf_max）：
| 模型 | Flash | 2m平均 | 70 cm平均 |
|---|
| Cosine | 8 B | 1013 Hz | 3305 Hz |
| Poly-3 | 16 B | 600 Hz | 1931 Hz |
| Poly-7（基线） | 32 B | 499 Hz | 1440 Hz |
| 3-bin Cheb-7 | 96 B | 313 Hz | 714 Hz |
| SGP4（完整） | >10 KB | 0 Hz | 0 Hz |
从Cosine到Poly-7，形状误差减少50 %（2m）/ 56 %（70 cm）。从Poly-7到3-bin，进一步减少37 %（2m）/ 50 %（70 cm）。
图3c 展示了总误差CDF（形状 + Δf_max估计）。关键统计量见表2。
表2：总误差预算（3-bin Cheb-7 + 高度感知Δf_max）
| 指标 | 2m（VHF） | 70 cm（UHF） |
|---|
| 平均最大误差 | 397 Hz | 868 Hz |
| 中位数最大误差 | 373 Hz | 905 Hz |
| 最差情况 | 883 Hz | 1753 Hz |
| P95 | 707 Hz | 1464 Hz |
| < 500 Hz | 79 % | 18 % |
| < 1000 Hz | 100 % | 67 % |
| < 2000 Hz | 100 % | 100 % |
图3d 将各方案置于Flash-精度Pareto前沿。3-bin Chebyshev-7（96字节）位于前沿的"膝点"
3.4 泛化性实验
3.4.1 90天泛化
使用2026年8月（第1–30天）数据训练的系数，在9–10月（第31–90天）数据上测试：
表3：90天泛化结果（3-bin Cheb-7）
| 数据集 | 2m 平均最大误差 | 2m 最差 | 70 cm 平均最大误差 | 70 cm 最差 |
|---|
| TRAIN（8月） | 397 Hz | 883 Hz | 868 Hz | 1753 Hz |
| TEST-1（9月） | 390 Hz (−2 %) | 872 Hz | 866 Hz (−0 %) | 1820 Hz |
| TEST-2（10月） | 381 Hz (−4 %) | 925 Hz | 922 Hz (+6 %) | 1697 Hz |
各卫星在两个测试月的误差与训练月几乎完全一致。90天内无明显泛化退化——证明了归一化径向速度曲线的物理本质是时间不变的：卫星过境的几何关系不随TLE老化而改变。
3.4.2 一年泛化
使用2026年8月训练的系数，在2027年8–9月（TLE年龄365–395天）数据上测试：
表4：一年泛化结果（3-bin Cheb-7）
| 数据集 | 2m 总误差 | 2m 形状误差 | 2m Δf估计误差 | 70 cm 总误差 | 70 cm 形状误差 | 70 cm Δf估计误差 |
|---|
| TRAIN（2026.08） | 397 Hz | 287 Hz | 234 Hz | 868 Hz | 691 Hz | 416 Hz |
| +1年（2027.08） | 413 Hz (+4 %) | 308 Hz (+7 %) | 244 Hz (+5 %) | 1165 Hz (+34 %) | 729 Hz (+6 %) | 730 Hz (+75 %) |
| +1年1月（2027.09） | 440 Hz (+11 %) | 298 Hz (+4 %) | 277 Hz (+18 %) | 1220 Hz (+41 %) | 725 Hz (+5 %) | 806 Hz (+94 %) |
表5：各卫星一年退化详情（70 cm波段）
| 卫星 | 高度 | TRAIN 平均 | +1年 平均 | 退化率 |
|---|
| ISS | 426 km | 800 Hz | 761 Hz | −5 %* |
| LilacSat-2 | 415 km | 1049 Hz | 2004 Hz | +91 % |
| PO-101 | 572 km | 767 Hz | 977 Hz | +27 % |
| SO-50 | 627 km | 804 Hz | 907 Hz | +13 % |
*ISS有定期轨道维持（reboost），因此高度几乎不变。
结果表明，
归一化曲线形状近乎时间不变：2m和70 cm的形状误差一年后分别仅增加7 %和6 %。这意味着Chebyshev多项式系数不需要更新，因为它们编码的是过境几何的纯物理关系。
退化主要来源于Δf_max估计的漂移：70 cm波段Δf估计误差从416 Hz增至730 Hz（+75 %），占总退化的绝大部分。根因是轨道衰减，低轨卫星（特别是<500 km的LilacSat-2，高度仅415 km，退化+91 %）的大气阻力使其轨道高度缓慢下降，从而改变了过境时长与峰值多普勒的映射关系。ISS（426 km）因定期轨道维持，退化反而为−5 %。
对实际部署的指导：每6–12个月从最新TLE重新拟合Δf_max估计系数（5个float，20字节）即可保持精度。形状多项式系数可永久使用。
3.5 最终算法
图4：最终算法
图4a 展示了三组仰角分段Chebyshev-7拟合叠加在归一化曲线集合上的效果。各组多项式捕获了其仰角范围内的特征形状：高仰角为较平坦的S曲线，低仰角为更非对称的曲线。
图4b 展示了LOSO CV下各卫星的总误差分布。卫星间差异不大，验证了模型的鲁棒性。
图4c 以SO-50（CV中未参与训练）为例，展示了3-bin模型与SGP4真值在多个仰角过境上的优秀一致性。
图4d 将所有误差置于±5 kHz FM频偏的上下文中：即使70 cm最差情况（1753 Hz）也仅占FM频偏的35 %，完全在典型FM接收机AFC捕获范围（±3 kHz）之内。
算法伪代码
// 一次性初始化（烧录固件时写入Flash）
// Chebyshev-7系数（每个仰角bin × 两个波段 = 6组 × 8float = 192 bytes）
// dur→Δf_max估计系数（每个波段5float = 40 bytes）
// 卫星高度表（10颗卫星 × 1float = 40 bytes）
// 每次过境初始化
function doppler_init(aos_unix, los_unix, sat_id):
g_aos = aos_unix
g_dur = los_unix - aos_unix // 过境时长
h = sat_altitude_table[sat_id] // 查卫星高度
d = g_dur
g_df_max = DUR_COEFF[0]*d² + DUR_COEFF[1]*d // 公式(8)
+ DUR_COEFF[2]*d*h + DUR_COEFF[3]*h
+ DUR_COEFF[4]
θ_est = estimate_elevation(g_dur, h) // 估算最大仰角
if θ_est ≤ θ_LOW: g_cheb = cheb_coeffs_low
elif θ_est ≤ θ_HIGH: g_cheb = cheb_coeffs_mid
else: g_cheb = cheb_coeffs_high
// 每秒调用（过境期间）
function doppler_update(now_unix):
t_norm = (now_unix - g_aos) / g_dur
if t_norm < 0 or t_norm > 1: return 0 // 过境尚未开始或已结束
x = 2 × t_norm - 1 // [0,1] → [-1,1]
// Clenshaw递推，公式(6)
bk2 = 0; bk1 = 0
for k = 7 down to 1:
bk = g_cheb[k] + 2 × x × bk1 - bk2
bk2 = bk1; bk1 = bk
norm = g_cheb[0] + x × bk1 - bk2 // ∈ [-1, 1]
return norm × g_df_max // Hz
每次更新的计算成本：1次除法（t_norm）+ Clenshaw递推（7次乘加+7次FMA）≈ Cortex-M0约30个周期。
3.6 定点数实现
Cortex-M0 无硬件浮点单元（FPU）。任何 float 操作都链接软件模拟库，大幅增加 Flash 占用（表5）。
表5：Arm GCC Cortex-M0 软件浮点库体积
| 函数 | 用途 | 体积 |
|---|
__addsf3 | float加法 | 504 B |
__subsf3 | float减法 | 428 B |
__mulsf3 | float乘法 | 416 B |
__divsf3 | float除法 | 524 B |
__fixsfsi | float→int | 112 B |
__floatsisf | int→float | 168 B |
__cmpsf2 | float比较 | 180 B |
| 合计 | | 2,332 B |
这 2.3 KB 的浮点库是真正的 Flash 杀手——比算法本身的数据和代码大一个数量级。
3.6.1 Q1.16 定点数格式
我们采用 Q1.16 定点数格式消除浮点依赖。所有值以有符号32位整数存储，低16位为小数部分：
$$x{\text{Q16}} = \lfloor x \cdot 2^{16} \rceil, \quad x = \frac{x{\text{Q16}}}{2^{16}}$$
- 取值范围：$[-32768, 32768 - 2^{-16}]$（覆盖归一化 Doppler $[-1, 1]$ 及 Chebyshev 中间值）
- 分辨率：$2^{-16} \approx 1.53 \times 10^{-5}$（相当于 10 kHz Δf_max 下 0.15 Hz）
- 乘法：
(int64_t)a * b >> 16（Cortex-M0 有 SMULL/UMULL 硬件指令，无需库）
- Clenshaw 中间值：
2*x*bk1 的三因子乘积需要 (int64_t)2 * x_q * bk1 >> 32，64位中间结果不溢出
3.6.2 位宽扫描
图5：定点数量化分析
图5a 展示了量化误差与小数位数的关系。Q8–Q16 范围内，量化噪声从 46 Hz 单调降至 0.2 Hz（70 cm）。Q18 及以上，Δf_max 值（10 kHz $\times 2^{18} > 2^{31}$）导致 int32 溢出，误差急剧上升。Q16 是平衡精度与不溢出的最优选择。
图5b 对比了 Flash 占用。软浮点方案总计约 2.7 KB（其中 2.3 KB 是浮点库），而定点 Q16 方案仅需约 500 字节（无浮点库依赖）。节省 91 %。
图5c 展示了 Q16 的量化噪声分布（vs float32 参考）。VHF 标准差 0.15 Hz，UHF 标准差 0.50 Hz——与模型本身的 400–900 Hz 误差相比完全可忽略。
图5d 将定点方案置于 Flash-精度 Pareto 前沿，Q16 3-bin Cheb-7 位于最优膝点。
3.6.3 Q16 定点系数表
表6：Q1.16 格式 Chebyshev-7 系数（70 cm 波段）
| bin | c0 | c1 | c2 | c3 | c4 | c5 | c6 | c7 |
|---|
| low (5–22°) | 15 | −69374 | −61 | 5090 | 22 | −797 | −7 | 178 |
| mid (22–55°) | −34 | −77378 | 17 | 16074 | −6 | −5089 | 5 | 2407 |
| high (55–90°) | −7 | −79773 | 16 | 20900 | −10 | −8055 | 11 | 5000 |
乘以 $2^{-16}$ 恢复浮点值。完整 Q16 系数及 C 头文件见附录 C。
3.6.4 定点 Clenshaw 递推（C 语言实现）
// Q1.16 Clenshaw 无浮点依赖，Cortex-M0 原生硬件指令
static inline int32_t clenshaw_q16(const int32_t c[8], int32_t x_q) {
int32_t bk2 = 0, bk1 = 0, bk;
int32_t two_q = 2 << 16;
for (int k = 7; k > 0; k--) {
// (int64)2 * x_q * bk1 >> 32 SMULL 硬件指令，无库调用
int32_t term = (int32_t)(((int64_t)two_q * x_q * bk1) >> 32);
bk = c[k] + term - bk2;
bk2 = bk1; bk1 = bk;
}
int32_t term_f = (int32_t)(((int64_t)x_q * bk1) >> 16);
return c[0] + term_f - bk2;
}
// Doppler update 过境期间每秒调用，返回整数 Hz
int32_t doppler_update(uint32_t now_unix) {
int32_t elapsed = (int32_t)(now_unix - (uint32_t)g_aos);
if (elapsed <= 0 || elapsed >= g_dur_sec) return 0;
int32_t t_q = (int32_t)(((int64_t)elapsed << 16) / g_dur_sec);
int32_t x_q = (t_q << 1) - (1 << 16); // x = 2t - 1
int32_t norm_q = clenshaw_q16(g_cheb, x_q);
return (int32_t)(((int64_t)norm_q * g_df_max_q) >> 16); // Hz
}
每次更新仅需 7 次 SMULL（64位乘）+ 约 20 条 ALU 指令，无任何库调用。
3.7 跨地点泛化
为检验模型对观测者位置的依赖性，在泰安（36°N）训练模型，于全球6个不同纬度地点测试（图6）。
图6：跨地点泛化实验
表7：跨地点泛化结果（3-bin Cheb-7 Q16）
| 地点 | 纬度 | 2m avg | 2m worst | 70cm avg | 70cm worst | vs Tai'an (70cm) |
|---|
| 泰安（基线） | 36°N | 377 Hz | 575 Hz | 815 Hz | 1428 Hz | — |
| 新加坡 | 1°N | 438 Hz | 820 Hz | 811 Hz | 1507 Hz | −0.4 % |
| 赫尔辛基 | 60°N | 378 Hz | 625 Hz | 830 Hz | 1307 Hz | +1.9 % |
| 圣地亚哥 | 33°S | 469 Hz | 872 Hz | 881 Hz | 1450 Hz | +8.1 % |
| 安克雷奇 | 61°N | 410 Hz | 698 Hz | 873 Hz | 1476 Hz | +7.1 % |
| 悉尼 | 34°S | 465 Hz | 836 Hz | 916 Hz | 1287 Hz | +12.4 % |
图6a–b 展示了各地点 UHF 和 VHF 波段的误差分布。70 cm 波段所有地点均在泰安基线的 13 % 以内，2m 波段在 25 % 以内。图6c 显示偏差百分比。图6d 以纬度为横轴，南半球误差略高于北半球（约 +100 Hz），主要由南半球过境几何的轻微非对称性所致，但仍远在 FM 可接受范围内。
泛化优良的物理原因已在 2.3 节阐明：归一化消除了频率依赖性，而不同纬度的主要差异（过境仰角分布）已被 3-bin 仰角分段模型覆盖。地球自转在不同纬度的线速度差异（赤道 460 m/s vs 极地 0 m/s）仅产生约 6 % 的视线速度修正，且该差异在归一化中被基本消除。
4. 结论
4.1 总结
本文证明了3-bin仰角分段Chebyshev-7多项式配合高度感知二次Δf_max估计器，可在资源受限的微控制器上实现亚kHz精度的FM卫星多普勒校正。采用Q1.16定点数格式消除浮点库依赖后，算法+系数总计仅需约500字节Flash（比SGP4节省100倍以上），无需每次过境上传数据，每次频率更新仅需约30个CPU周期。模型在泰安训练后可全球泛化，6个纬度地点的UHF误差偏差均<13 %。一年泛化实验表明归一化曲线形状近乎时间不变，Δf_max估计需每6–12个月用最新TLE更新一次。
4.2 核心发现
- 多普勒曲线的形状就是归一化径向速度：$f_{\text{carrier}}$ 在归一化中抵消，一组多项式系数适用于所有频率。
- UHF误差 = 3 × VHF误差：这是多普勒物理定律的必然结果，任何模型（包括完整SGP4）都无法规避。
- Chebyshev与Monomial精度完全一致：推荐Chebyshev是因为Clenshaw递推在float32下数值更稳定，且系数衰减快可安全截断。
- 最大仰角是曲线形状的最强预测因子：按仰角分段（3-bin）比单一多项式减少形状误差50 %。
- 加入卫星高度信息将Δf_max估计的 $R2$ 从0.89提升至0.98（70 cm波段），消除了不同高度卫星的歧义。
- TCA信息增益有限：TCA偏移与仰角相关，但仰角是形状的直接驱动力，TCA仅是副产品。
- 归一化曲线形状是时间不变的：一年后形状误差仅退化6 %，Chebyshev多项式系数可永久使用。
- Δf_max估计随轨道衰减退化：一年后Δf估计误差增加75 %（70 cm），建议每6–12个月从最新TLE更新一次（5个Q16系数，20字节）。
- 定点数量化几乎无损：Q1.16格式下量化噪声仅0.5 Hz（UHF），与模型误差相比可忽略。消除浮点库后Flash从2.7 KB降至500字节。
- 全球泛化无需重训练：在泰安（36°N）训练的系数直接用于赤道至60°N/S纬度，70 cm误差偏差<13 %.
4.3 实际部署建议
算法适合集成到业余电台固件中：
- 一次性烧录：Chebyshev Q16系数（6组×8×4=192字节）、dur2Δf Q16系数（2组×5×4=40字节）、卫星高度表（10颗×16=160字节）、代码（180字节）。总计约572字节（含卫星表），无浮点库依赖。
- 每次过境设置：用户从菜单选择卫星，输入AOS和LOS时间。MCU从时长和高度估算Δf_max和最大仰角，选择对应Chebyshev系数组。
- 过境期间：每秒计算 $t_{\text{norm}}$（Q16整数除法），通过Clenshaw递推求值Chebyshev多项式（7次64位乘加），乘以Δf_max（1次64位乘+移位），施加频率校正。全部使用Cortex-M0硬件整数指令。
4.4 局限性
- 本模型针对近圆LEO轨道验证（高度400–1400 km，偏心率<0.05）。高椭圆轨道（如AO-10、AO-7）需单独刻画。
- 3倍UHF/VHF误差比是多普勒物理的固有约束，无法通过任何模型降低，仅可通过使用更低频段改善。
- 近顶过境（>85 °仰角）产生的急剧多普勒变化对任何多项式模型构成挑战。3-bin方案通过专用高仰角多项式部分缓解了该问题。
- TLE每6–12个月需更新一次以维持Δf_max估计精度。
- 南半球误差略高于北半球（+12 %），但绝对值增量（100 Hz）对FM通信无实质影响。
附录A：Chebyshev系数表
A.1 70 cm波段，3-bin Chebyshev-7
// 低仰角 (5–22°), CV平均最大误差 569 Hz
static const float cheb7_70cm_low[8] = {
+0.0002335556f, -1.0585669335f, -0.0009278245f, +0.0776628764f,
+0.0003356745f, -0.0121636282f, -0.0001125546f, +0.0027210660f
};
// 中仰角 (22–55°), CV平均最大误差 752 Hz
static const float cheb7_70cm_mid[8] = {
-0.0005193342f, -1.1806953138f, +0.0002602060f, +0.2452703434f,
-0.0000972926f, -0.0776532436f, +0.0000699958f, +0.0367319887f
};
// 高仰角 (55–90°), CV平均最大误差 820 Hz
static const float cheb7_70cm_high[8] = {
-0.0001019328f, -1.2172468235f, +0.0002460633f, +0.3189162096f,
-0.0001456953f, -0.1229115857f, +0.0001671119f, +0.0762920805f
};
A.2 70 cm波段——dur→Δf_max估计系数
// Δf_max = c[0]·dur² + c[1]·dur + c[2]·dur·alt + c[3]·alt + c[4]
static const float dur2df_70cm[5] = {
0.00114662f, 35.69112240f, -0.03579820f, 1.84041431f, -978.11635070f
};
A.3 卫星高度参考表
// 从TLE平均运动反算: a = (GM/(n/60)²)^(1/3), h = a - 6371 km
static const float sat_altitudes[10] = {
// SO-50 AO-91 AO-85 AO-95 LilacSat-2 ISS PO-101 AO-73 NO-44 RS-44
627, 533, 601, 490, 415, 426, 572, 543, 792, 1350
};
A.4 仰角阈值
2m波段：$\theta{\text{low}} \leq 26.1\circ$，$\theta{\text{high}} > 56.7\circ$
70 cm波段：$\theta{\text{low}} \leq 22.1\circ$，$\theta{\text{high}} > 54.7\circ$
附录B：Clenshaw递推C语言实现
// 在 t_norm ∈ [0, 1] 上求值Chebyshev级数
// coeffs[0..7] 对应于 x ∈ [-1, 1] 上的 T_0 至 T_7
float chebyshev_eval(const float coeffs[8], float t_norm) {
float x = 2.0f * t_norm - 1.0f; // [0,1] → [-1,1]
float bk2 = 0.0f, bk1 = 0.0f, bk;
for (int k = 7; k > 0; k--) {
bk = coeffs[k] + 2.0f * x * bk1 - bk2;
bk2 = bk1;
bk1 = bk;
}
return coeffs[0] + x * bk1 - bk2; // ∈ [-1, 1]
}
// 完整多普勒更新 —— 过境期间每秒调用
float doppler_update(uint32_t now_unix) {
float t_norm = (float)(now_unix - g_aos_time) / g_duration_sec;
if (t_norm < 0.0f || t_norm > 1.0f) return 0.0f; // 非活跃
// 根据仰角bin选择Chebyshev系数组
const float *cheb;
if (g_max_el <= EL_THRESH_LOW)
cheb = cheb7_low;
else if (g_max_el <= EL_THRESH_HIGH)
cheb = cheb7_mid;
else
cheb = cheb7_high;
float norm = chebyshev_eval(cheb, t_norm);
return norm * g_df_max; // Hz
}
附录C：Q1.16 定点系数表
C.1 2m 波段 — 3-bin Chebyshev-7（Q1.16）
// 低仰角 (5–26°), CV平均最大误差 215 Hz
static const int32_t cheb7_2m_low[8] = {
-869, -70013, 376, 5967, -58, -944, 10, 216
};
// 中仰角 (26–57°), CV平均最大误差 328 Hz
static const int32_t cheb7_2m_mid[8] = {
-204, -77140, 280, 15652, -88, -4802, 24, 2187
};
// 高仰角 (57–90°), CV平均最大误差 406 Hz
static const int32_t cheb7_2m_high[8] = {
-224, -79212, 467, 19698, -226, -7196, 159, 4131
};
C.2 70cm 波段 — 3-bin Chebyshev-7（Q1.16）
// 低仰角 (5–22°), CV平均最大误差 569 Hz
static const int32_t cheb7_70cm_low[8] = {
15, -69374, -61, 5090, 22, -797, -7, 178
};
// 中仰角 (22–55°), CV平均最大误差 752 Hz
static const int32_t cheb7_70cm_mid[8] = {
-34, -77378, 17, 16074, -6, -5089, 5, 2407
};
// 高仰角 (55–90°), CV平均最大误差 820 Hz
static const int32_t cheb7_70cm_high[8] = {
-7, -79773, 16, 20900, -10, -8055, 11, 5000
};
C.3 dur→Δf_max 估计系数（Q1.16）
// 2m: Δf_max = c0·dur² + c1·dur + c2·dur·alt + c3·alt + c4
static const int32_t dur2df_2m[5] = {
-156, 461355, -27, -140350, 81041129
};
// 70cm: 同上公式
static const int32_t dur2df_70cm[5] = {
75, 2339053, -2346, 120613, -64101833
};
C.4 仰角阈值
| 波段 | θlow | θ_high |
|---|
| 2m | ≤ 26.1° | > 56.7° |
| 70cm | ≤ 22.1° | > 54.7° |
数据：1420次过境，10颗卫星，30天训练窗口，留一卫星交叉验证。
TLE来源：Celestrak业余无线电组（2026-08-03）。
泛化实验：90天（至2026年10月）+ 1年（至2027年9月），TLE传播最长395天。
跨地点实验：6个纬度（1°N–61°N, 33°S–34°S），14天测试窗口。
ps. 最近在开发Baofeng UV-K61的固件（对，是宝锋 K61！该有的功能都有），预告一下，敬请期待喵~
pps. HamCQ怎么不支持LaTeX哇
闽ICP备2021006864号-7
闽公网安备35020602002794号