基於 Tschebyschev 多項式的輕量化業餘衛星多普勒頻移校正演算法
500 位的 Flash 記憶體，與 SGP4 具有相同的精確度。
摘要
在低軌道的衛星通訊中，多普勒頻移校正通常需要運行 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 公里的近圓軌道，速度約 7.5 公/秒，產生的多普勒頻移在 145 MHz（2 米波段）可達 ±3.5 kHz，在 435 MHz（70 厘米波段）可達 ±10.5 kHz。若不進行實時頻率校正，接收訊號數秒內即超出接收機通帶。
最近在開發 Baofeng UV-K61 的韌體時，發現在使用低階業餘無線電微控制器上運行 SGP4 軌道預報器時，需要儲存完整的 TLE 星曆並實作 SGP4 演算法，導致 Flash 使用量超過 30 KB（考慮到軟浮點數程式碼的佔用）。對於只有 64-128 KB 的 Cortex-M0 無線電，這迫使開發者做出困難的選擇：要么放棄其他無線電功能來運行專用的 SGP4 韌體，或使用電腦預先下載過境資料。
1.2 目前的方案及其缺點
目前，主要的開源韌體採用兩種策略，以 Quansheng UV-K5 為例（不考慮高性能 DMR 手機）：
- 完整功能的 SGP4 韌體將所有 Flash 功能用於 SGP4，放棄其他無線電功能。
- 預先計算的查找表（LUT）：每次過境前，透過 PC 透過串接埠上傳校正表。 需要電腦，並且每次過境都需要重新燒寫。
1.3 本文的貢獻
我們證明：在衛星經過時，多普勒曲線可以被…一個階的切比雪夫多項式精確估算，僅需要使用者已知的三個參數：衛星辨識碼、AOS時間、LOS時間。此演算法：
- 採用 Q1.16：定點數，消除浮點數庫依賴，僅需要 大約 500 節 Flash (相較於 SGP4，節省超過 100 倍)
- 達到 397 赫茲（VHF）868 赫茲（UHF） 平均誤差 (與 SGP4 真值的比較)
- 每次頻率更新僅需要 30個 CPU 週期（Cortex-M0 僅支援整數指令）
- 無需每次過境時上傳資料——使用者只需選擇衛星，並輸入其經過時間。
- 歸一化的曲線形狀具有時間不變性——一年後，形狀誤差僅下降 6%。
- 全球性——一套係數適用於所有緯度（UHF 偏差）<13 %）
2. 數據與方法
2.1 資料來源
TLE 星表：2026年8月3日從 星際網路業餘無線電組 取得常用的 10 顆 FM 業餘衛星的最新 TLE 資料。
觀測站：位於山東省泰安市（北緯36.18°、東經117.13°）。
分析衛星：
| 衛星 | NORAD | 向下頻率 | 高度 | 30 天的過境次數 |
|---|
| SO-50（沙特衛星-1C） | 27607 | 436.795 ميجاهرتز | 627 公里 | 163 |
| AO-91（RadFxSat/Fox-1B） | 43017 | 145.960 ميجاهرتز | 533 公里 | 109 |
| AO-85（Fox-1A） | 40967 | 145.980 ميجاهرتز | 601 公里 | 138 |
| AO-95（Fox-1Cliff） | 43770 | 145.940 ميجاهرتز | 490 公里 | 112 |
| LilacSat-2 (CAS-3H) | 40908 | 437.200 ميجاهرتز | 415 公里 | 103 |
| 國際空間站（「雅羅」） | 25544 | 437.800 ميجاهرتز | 426 公里 | 208 |
| PO-101（迪瓦塔-2） | 43678 | 437.500 ميجاهرتز | 572 公里 | 150 |
| AO-73 (FUNcube-1) | 39444 | 145.950 ميجاهرتز | 543 公里 | 118 |
| NO-44 (PCSAT) | 26931 | 145.825 ميجاهرتز | 792 公里 | 157 |
| RS-44 | 44909 | 145.935 ميجاهرتز | 1350 公里 | 162 |
總計：共發生1420次過境（其中VHF 796次，UHF 624次），時間窗口為2026年8月3日至2026年9月2日。
2.2 真值生成
針對每次過境，使用 SGP4 模型進行計算。skyfield Python 函式庫 (封裝標準)sgp4（包）傳播 TLE，在 150 個均勻的時間採樣點計算真實的多普勒頻移：
Δf{\text{真}}(t) = -f載波} · $\frac{v_r(t)}{c}$ (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{載波}}/c$。由於 $fUHF/f<sub>VHF</sub> ≈ 3，UHF 的誤差通常是 VHF 的 3 倍。
2.4 簡略模型
我們評估了六種複雜度遞增的模型：
模型 A：餘弦基準
$$\Delta \hat{f}(t) = \Delta f$$最大值 · cos(π·t)(規範值),
t{\text{規範}} = \frac{t - t}{AOS}{\text{LOS}} - t{\text{AOS}} ∈ \[0, 1] (3)
模型 B：單項式多項式 ( $d$ 階)
$$\Delta \hat{f}(t) = \Delta \hat{f}$$最大值 · Σ{i從0到d} * c\_i * t\_normi \tag{4}
係數 $c_i$ 透過最小二乘法對歸一化的 SGP4 曲線集合進行擬合而得。使用 Horner 方法來計算。
模型 C：切比雪夫多項式 ( $d$ 階)
$$\Delta \hat{f}(t) = \Delta \hat{f}$$最大值 · Σ$\sum_{i=0}^{d} a_i \cdot T_i(x)$, 其中 $x = 2 t_{\text{norm}} - 1 \in [-1, 1]$ (5)
其中，$T_i(x) = \cos(i \cdot \arccos(x))$ 是一個在 $[-1, 1]$ 區間內的 birinci類 Chebyshev 多項式。 系數 $a_i$ 透過最小二乘法進行擬合得到。 計算方式採用克倫沙（Clenshaw）遞迴：
$$\begin{cases} b{d+2} = b{d+1} = 0 \ b\_k = a\_k + 2x · b{k+1} - b{k+2}，其中 k 為 d、d-1、...、1 \ \[
\text{norm} = a_0 + x \cdot b_1 - b_2
\] (6)
在 Cortex-M0 的 float32 精確度下，Clenshaw 遞迴相對於同階的 Horner 方法，其結果更加穩定。
模型 D：3 個隔間的仰角分段 Chebyshev-7
將過境曲線依照最大仰角分為三組（低/中/高，以訓練集的第33和67百分位數作為分界點），並對每組獨立擬合Chebyshev-7多項式：
$$\Delta \hat{f}(t) = \Delta \hat{f}$$最大值 · Σ{i=0}^{7} a\_i^{(\text{bin})} \cdot T\_i(2t){\text{norm}} - 1), \quad \text{bin} = \begin{cases} \text{低} & θ{\max} ≤ θ{33} \ 中點 & θ{33} < θ$\{\max\} \leq \theta${67} \ 高 & θ最大值 > θ{67} \end{cases} (7)$$
最大 Δf 估計模型
所有多項式模型都需要估算峰值多普勒變 ($\Delta \hat{f}_{\max}$). 我們使用高階感知二元模型：
$$\Delta \hat{f}_{\max} = \alpha_1 T2 + α₂T + α₃Th + α₄h + α₅ (式 8)
其中，$T = t${\text{LOS}} - t{\text{AOS}} 代表過境時間（秒），$h$ 代表衛星高度（公里）。此模型捕捉了以下物理事實：過境時間和峰值多普勒頻率都取決於衛星高度和最大仰角，只要已知衛星身份（從而得知其高度），就能消除「不同衛星相同過境時間對應不同的Δf_max」的歧義。
2. 驗證方法
主要實驗：採用使用衛星進行交叉驗證 (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真實值的差異。最大絕對誤差。報告這些超過邊界的最大誤差值的平均值、中位數、第95百分比值以及最差值。
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 曲線與樣本內曲線幾乎重疊，確認即使到 15 次，也沒有發生過度擬合。在相同階數下，Chebyshev 與 Monomial 的精確度完全一致。——兩者都使用了相同的多項式空間，但 Chebyshev 系數的衰減速度較快，且在 float32 下的數值更穩定。
3.2 多普勒曲線歸一化的物理本質
圖 2：物理機制與仰角依賴性
圖 2a 展示了徑向速度 $v_r(t)$ 的幾何：徑向速度向量在觀察者視線方向上的投影。在AOS和LOS時刻（低仰角），速度幾乎與視線垂直，多普勒變化的速率較慢。而在TCA附近（高仰角），速度與視線對齊，多普勒變化的速率非常快。
圖 2b 已驗證三倍縮放定律將 145 MHz 的多普勒曲線乘以 3，使其與在相同過境幾何下的 435 MHz 曲線精確對應。這驗證了 $\text{norm}(t) = v_r(t)/v_r(\text{AOS})$ 在頻率上的不變性。
圖 2c 揭示了一個關鍵發現：針對同一個衛星（SO-50），歸一化曲線的形狀隨最大仰角劇烈變化。9° 的仰角過境高度不對稱（遠離中心點），而 87° 的仰角則幾乎對稱。形狀差異達達到滿量程的 47%。
圖 2d 依照最大仰角對所有歸一化的曲線進行著色，可以直觀地顯示仰角是曲線形狀變異的主要驅動因素。高仰角的過境（以黃色表示）集中在對稱的S型曲線周圍，而低仰角的過境（以紫色表示）則分散開來。
3.3 模型的比較
圖 3：模型比較
圖 3a 呈現了過境時間與多普勒峰值的關係。高頻度二維模型在 70 厘米波段達到 R 值。2 = 0.982 美元 (均方根殘差 397 赫茲)，而不包含高頻資訊的波段級二次模型僅為 $R2 = 0.889 美元。各衛星的 Δf_max 估計殘差，請見表 1。
表格 1： Δf_max 估計殘差 (70 公分波段，高頻度感知模型)
| 衛星 | 平均殘差 | 標準差 | 最大誤差 |
|---|
| ISS | -318 赫茲 | 185 赫茲 | 635 赫茲 |
| LilacSat-2 | +482 赫茲 | 133 赫茲 | 687 赫茲 |
| PO-101 | +340 赫茲 | 260 赫茲 | 782 赫茲 |
| SO-50 | -212 赫茲 | 246 赫茲 | 614 赫茲 |
圖 3b 比較了各個模型尺寸誤差（CV，假設理想的 Δf_max）：
| 模型 | Flash | 2 米平均 | 70 公分 |
|---|
| Cosine | 8 B | 1013 赫茲 | 3305 赫茲 |
| Poly-3 | 16B | 600 赫茲 | 1931 赫茲 |
| 多利-7 (基準) | 32 B | 499 赫茲 | 1440 赫茲 |
| 3 槽 Cheb-7 | 96 B | 313 赫茲 | 714 赫茲 |
| SGP4 (完整) | >10 KB | 0 赫茲 | 0 赫茲 |
從 Cosine 到 Poly-7，形狀誤差減少 50% (2 公尺) / 56% (70 公分)。從 Poly-7 到 3-bin，進一步減少 37% (2 公尺) / 50% (70 公分)。
圖 3c 展示了總誤差的累積分佈函數（形狀 + Δf_max 的估計）。關鍵統計數據請見表 2。
表格 2：總誤差預算 (3-bin Cheb-7 + 深度感知 Δf_max)
| 指標 | 2m（VHF） | 70 公分（UHF） |
|---|
| 平均最大誤差 | 397 赫茲 | 868 赫茲 |
| 中位數的最大誤差 | 373 赫茲 | 905 赫茲 |
| 最壞情況 | 883 赫茲 | 1753 赫茲 |
| P95 | 707 赫茲 | 1464 赫茲 |
| < 500 赫茲 | 79 % | 18 % |
| < 1000 赫茲 | 100 % | 67 % |
| < 2000 赫茲 | 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)
| 數據集 | 2 米平均最大誤差 | 2 米，最差 | 70 公分平均最大誤差 | 70 公分最差 |
|---|
| TRAIN (八月) | 397 赫茲 | 883 赫茲 | 868 赫茲 | 1753 赫茲 |
| 測試一 (九月) | 390 赫茲 (-2%) | 872 赫茲 | 866 赫茲 (–0%) | 1820 赫茲 |
| 測試-2 (10月份) | 381 赫茲 (-4%) | 925 赫茲 | 922 赫茲（+6%） | 1697 赫茲 |
所有衛星在兩個測試月的誤差與訓練月幾乎完全一致。90天內沒有明顯的泛化退化——這證明了歸一化的徑向速度曲線的物理本質是時間不變的：衛星經過時的幾何關係不會隨著 TLE 的老化而改變。
3.4.2 年度概括
使用於2026年8月的訓練參數，在2027年8～9月（TLE年齡為365～395天）的數據上進行測試：
表格 4：一年泛化的結果 (3-bin Cheb-7)
| 數據集 | 總誤差 2 公尺 | 2 米的形狀誤差 | 2 米 Δf 估計誤差 | 總誤差 70 公分 | 70 公分尺寸誤差 | 70 公分 Δf 的估計誤差 |
|---|
| TRAIN（2026.08） | 397 赫茲 | 287 赫茲 | 234 赫茲 | 868 赫茲 | 691 赫茲 | 416 赫茲 |
| +1 年（2027.08） | 413 赫茲（+4%） | 308 赫茲（+7%） | 244 赫茲（+5%） | 1165 赫茲（+34%） | 729 赫茲（+6%） | 730 赫茲（+75%） |
| +1 年 1 月（2027.09） | 440 赫茲（+11%） | 298 赫茲（+4%） | 277 赫茲（+18%） | 1220 赫茲（+41 %） | 725 赫茲（+5%） | 806 赫茲（+94 %） |
表 5：各衛星一年的衰變詳情 (70 公分波段)
| 衛星 | 高度 | 列車平均 | +1 年，平均 | 退化速率 |
|---|
| ISS | 426 公里 | 800 赫茲 | 761 赫茲 | −5 %* |
| LilacSat-2 | 415 公里 | 1049 赫茲 | 2004 赫茲 | +91 % |
| PO-101 | 572 公里 | 767 赫茲 | 977 赫茲 | +27 % |
| SO-50 | 627 公里 | 804 赫茲 | 907 赫茲 | +13 % |
*國際太空站 (ISS) 會定期進行軌道維持 (reboost)，因此高度幾乎不會改變。
結果顯示，
歸一化的曲線形狀幾乎與時間無關。：2公尺和70公分的形狀誤差，在一年後分別僅增加7%和6%。這表示Chebyshev多項式的係數不需要更新因為它們編碼的是過境幾何的純物理關係。
退化的主要原因是 Δf_max 的估計產生偏差。：70 厘米波段的 Δf 估計誤差從 416 Hz 增加到 730 Hz （+75%），佔總衰減的大部分。根本原因是軌道衰減，低軌衛星（特別是<LilacSat-2 的軌道距離地球約 500 公里，高度僅 415 公里（衰減 +91%），因此受到大氣阻力的影響，軌道高度會緩慢下降，進而改變了穿越時長與多普勒效應的關係。ISS (426 公里) 因為定期軌道維持，衰減反而為 -5%。
實際部署的指導：只需每 6-12 個月從最新的 TLE 中重新調整 Δf\_max 的估計係數（5 個浮點數，20 字节），即可維持精確度。多項式係數則可以永久使用。
3.5 最終演算法
圖 4：最終演算法
圖 4a 展示了三組仰角分段的Chebyshev-7擬合疊加在歸一化曲線集合上的效果。各組多項式捕捉了其仰角的範圍內的特徵形狀：高仰角為較平坦的S型曲線，低仰角則為更非對稱的曲線。
圖 4b 呈現了 LOSO CV 下各衛星的總誤差分佈。 衛星間差異不大，驗證了模型的可靠性。
圖 4c 以 SO-50 (CV 中未參與訓練) 為例，展示了 3-bin 模型與 SGP4 真值的多個仰角過境上的優異一致性。
圖 4d 將所有誤差置於 ±5 kHz 的 FM 頻偏範圍內：即使在 70 公分最差情況下 (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 或 t\_norm > 1: 返回 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
每次更新的計算成本：一次除法 (t_norm) + 克林沙遞迴 (7次乘加 + 7次 FMA) ≈ 大約 Cortex-M0 的 30 個時鐘。
3.6 使用定點數實現
Cortex-M0 沒有內建的硬體浮點運算器（FPU）。任何 float 這些操作都與軟體模擬庫連結，顯著增加了 Flash 的佔用空間（見表 5）。
表 5：Arm GCC Cortex-M0 軟體浮點運算函式庫的大小
| 函數 | 用途 | 體積 |
|---|
__addsf3 | 浮點數的加法 | 504 B |
__subsf3 | 浮點數減法 | 428 B |
__mulsf3 | 浮點數的乘法 | 416 B |
__divsf3 | 浮點數除法 | 524 B |
__fixsfsi | 浮點數 → 整數 | 112 B |
__浮動 | 整數 → 浮點數 | 168 B |
__cmpsf2 | 浮點數比較 | 180 B |
| 總計 | | 2,332 B |
這個 2.3 KB 的浮點數函式庫，確實是「Flash」的終結者——其大小遠遠超過原本的演算法本身所包含的資料和程式碼。
3.6.1 Q1.16 固定點數格式
我們採用 Q1.16 採用固定精確度格式，避免使用浮點數。所有值皆以有符號 32 位整數儲存，其中低 16 位為小數位：
$$x{\text{Q16}} = \left\lceil x \cdot 2^{16} \right\rceil, \quad x = \frac{x}{Q16 2的16次方
- 允許的範圍：\[-32768, 32768 - 2^{-16}]（涵蓋歸一化的多普勒 $[-1, 1]$ 以及切比雪夫的中間值）
- 解析度：$2^{-16} \approx 1.53 \times 10^{-5}$ (相當於 10 kHz 的 Δf\_max，即 0.15 Hz)
- 乘法：
(int64_t) a \* b >> 16（Cortex-M0 具有） SMULL/UMULL 硬體指令，無需庫
- 克倫沙爾的中間值：
2\*x\*bk1 三因子的乘積需要 (int64_t) 2 * x\_q * bk1 >> 32，64 位中間結果不會溢出
3.6.2 寬度掃描
圖 5：固定點數量化分析
圖 5a 展示了量化誤差與小數位數之間的關係。在 Q8-Q16 的範圍內，量化雜訊從 46 Hz 單調下降到 0.2 Hz (70 公分)。Q18 以上，最大頻率變化值 (Δf_max)10 千赫茲 × 2^18 > 使用 $2^{31}$ 導致 `int32` 發生溢位，誤差急劇增加。 Q16 在精確度和避免溢位的最佳選擇。
圖 5b 比較了 Flash 的佔用空間。使用軟浮點方案，總共約占用 2.7 KB (其中 2.3 KB 是浮點庫)，而使用定點 Q16 方案則僅需約 500 字节 (不依賴浮點庫)。節省 91%。
圖 5c 展示了 Q16 的量化雜訊分佈（與 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 公分波段)
| bin | c0 | c1 | c2 | c3 | c4 | c5 | c6 | c7 |
|---|
| 低（5–22°） | 15 | −69374 | −61 | 5090 | 22 | −797 | −7 | 178 |
| 中等（22–55°） | −34 | −77378 | 17 | 16074 | −6 | −5089 | 5 | 2407 |
| 高（55–90°） | −7 | −79773 | 16 | 20900 | −10 | −8055 | 11 | 5000 |
將結果乘以 $2^{-16}$ 以恢復浮點數值。完整的 Q16 系數及 C 標頭檔請參閱附件 C。
3.6.4 使用 C 語言實現的定點克倫沙 (Clenshaw) 遞迴
// 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;
}
// 測速更新：在過境期間每秒呼叫，返回整數 Hz
int32_t doppler_update(uint32_t now_unix) {
int32_t elapsed = (int32_t)(now_unix - (uint32_t)g_aos);
if (elapsed <= 0 或 経過時間 >= g_dur_sec) return 0;
int32_t t_q = (int32_t)(((int64_t)経過時間 << 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) // 赫茲
}
每次更新只需要進行 7 次。 SMULL（64 位乘法）+ 約 20 個 ALU 指令，沒有任何函式庫呼叫。
3.7 跨地點概括
為了驗證模型對觀察者位置的依賴性，我們在泰安（36°N）訓練模型，並在全球六個不同緯度的地點進行測試（圖6）。
圖 6：跨地點泛化實驗
表格 7：跨地點泛化的結果 (3-bin Cheb-7 Q16)
| 地點 | 緯度 | 平均 2 公尺 | 最差情況：2公尺 | 平均身高 70 公分 | 70 公分粗紗 | 與泰安（70公分）比較 |
|---|
| 泰安 (基地) | 36 度北緯 | 377 赫茲 | 575 赫茲 | 815 赫茲 | 1428 赫茲 | — |
| 新加坡 | 1 度北緯 | 438 赫茲 | 820 赫茲 | 811 赫茲 | 1507 赫茲 | −0.4 % |
| 赫爾辛基 | 60 度北緯 | 378 赫茲 | 625 赫茲 | 830 赫茲 | 1307 赫茲 | +1.9 % |
| 聖地牙哥 | 33 度南緯 | 469 赫茲 | 872 赫茲 | 881 赫茲 | 1450 赫茲 | +8.1 % |
| 安克雷奇 | 61 度北緯 | 410 赫茲 | 698 赫茲 | 873 赫茲 | 1476 赫茲 | +7.1 % |
| 悉尼 | 34° 南緯 | 465 赫茲 | 836 赫茲 | 916 赫茲 | 1287 赫茲 | +12.4 % |
圖 6a–b 呈現了各地點在 U波段和 V波段中的誤差分佈情況。在 70 公分的波段，所有地點的誤差值皆低於泰安基準的 13%；而在 2 公分的波段，誤差值則低於 25%。圖 6c 顯示偏差百分比。圖 6d 以緯度為橫軸，南半球的誤差略高於北半球（約 +100 Hz），主要原因是南半球經過的地理上的輕微非對稱性，但仍然遠低於 FM 的可接受範圍。
如同 2.3 節所闡述的，將優良的物理原因泛化：歸一化消除了頻率依賴性，而不同緯度主要差異（過境仰角分布）已被 3-bin 仰角分段模型涵蓋。地球自轉在不同緯度的線速度差異（赤道 460 m/s vs 極地 0 m/s）僅產生約 6 % 的視線速度修正，且該差異在歸一化中被基本消除。
4. 結論
4.1 總結
本文證明了3-bin 仰角分段的七次齊柏羅多項式配合高階感知二次 Δ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的遞迴演算法在浮點數32位的情況下，數值更穩定，且係數衰減快速，因此可以安全地截斷。
- 最大仰角是最強的預測因數。：採用仰角分段 (3-bin) 方法，相較於單一多項式，可以減少約 50% 的形狀誤差。
- 加入衛星高度資訊將 Δf\_max 估算的 R2$ 從 0.89 提升至 0.98 (70 公分波段)，已消除不同高度衛星的歧義。
- TCA 資訊增益有限：TCA 的偏移與俯仰角有關，但俯仰角才是形狀的直接驅動因素，而 TCA 只是輔助效果。
- 歸一化的曲線形狀是時間不變的。：一年後，形狀誤差僅下降 6%，Chebyshev 多項式的係數可以永久使用。
- Δf_max 的估計會隨著軌道的衰減而變化。：一年後，Δf的估計誤差增加了75%（70 公分），建議每6-12個月從最新的TLE更新一次（包含5個Q16係數，20位元組）。
- 定點數量化幾乎沒有損害。：在 Q1.16 格式下，量化噪音僅為 0.5 Hz。（UHF）與模型誤差相比，可以忽略不計。移除浮點數庫後，Flash從 2.7 KB 減少到 500 字节。
- 無需在全世界範圍內重新訓練：在泰安（36°N）訓練的係數可以直接應用於從赤道到60°N/S緯度的區域，誤差範圍為70公分。<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 公里，離心率)<0.05）。高橢圓軌道（例如 AO-10、AO-7）需要單獨繪製。
- 3倍的UHF/VHF誤差比是多普勒物理本身的限制，無法透過任何模型來降低，只能透過使用較低的頻率範圍來改善。
- 接近頂點（>（85° 仰角）產生的劇烈多普勒效應對任何多項式模型都構成挑戰。3-bin 方案透過專用的高仰角多項式部分，來緩解這個問題。
- 為了維持 Δf\_max 的估計精確度，需要將 TLE 更新至少每 6-12 個月一次。
- 南半球的誤差略高於北半球（+12%），但絕對值的增加（100 （以赫茲為單位）對 FM 通訊沒有實質影響。
附錄 A：切比雪夫係數表
A.1 70 公分頻寬，3 區 Chebychev-7
// 低仰角 (5–22°), CV平均最大誤差 569 Hz
static const float cheb7_70cm_low[8] = {
+0.0002335556, -1.0585669335, -0.0009278245, +0.0776628764,
+0.0003356745, -0.0121636282, -0.0001125546, +0.0027210660
};
// 中仰角 (22–55°), CV平均最大誤差 752 Hz
static const float cheb7_70cm_mid[8] = {
-0.0005193342, -1.1806953138, +0.0002602060, +0.2452703434,
-0.0000972926, -0.0776532436, +0.0000699958, +0.0367319887
};
// 高仰角 (55–90°), CV平均最大誤差 820 Hz
static const float cheb7_70cm_high[8] = {
-0.0001019328, -1.2172468235, +0.0002460633, +0.3189162096,
-0.0001456953, -0.1229115857, +0.0001671119, +0.0762920805
};
A.2 70 公分頻寬——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 公里
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 仰角閾值
2 ميجاه茲頻段：$\theta$低 ≤ 26.1〇$, $\theta高 > 56.7〇$
70 公分波段：$\theta$低 ≤ 22.1〇$, $\theta高 > 54.7〇$
附錄 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];
bk2 = bk * bk - bk1;
bk1 = bk;
x += bk2;
}
return x;
} > 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 2 ميﮕ赫頻段 — 3個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 70 公分頻寬 — 3 節 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）
// 2公尺：Δf_max = c0·dur² + c1·dur + c2·dur·alt + c3·alt + c4
static const int32_t dur2df_2m[5] = {
-156, 461355, -27, -140350, 81041129
};
// 70公尺：同上公式
static const int32_t dur2df_70cm[5] = {
75, 2339053, -2346, 120613, -64101833
};
C.4 仰角閾值
| 頻寬 | θlow | θ\_高 |
|---|
| 2m | ≤ 26.1° | > 56.7° |
| 70cm | ≤ 22.1° | > 54.7° |
數據：1420次過境，10顆衛星，30天的訓練窗口，使用剩餘一顆衛星進行驗證。
TLE 來源：Celestra 業餘無線電組（2026-08-03）。
一般性實驗：90天（至2026年10月）+ 1年（至2027年9月），TLE傳播時間最長可達395天。
跨地點實驗：6個緯度（1°北—61°北、33°南—34°南），14天的測試期間。
P.S. 近期正在開發寶鋒 UV-K61 的韌體 (是的，就是寶鋒 K61！它應該具備的所有功能都包含在內)，請稍後關注，敬請期待～
P.S. HamCQ 怎麼不支持 LaTeX 呢？
闽ICP备2021006864号-7
闽公网安备35020602002794号