基于影錐特征角的近圓軌道衛星進出地影時間的快速判定法

謝文杰，邵圣祥，王 霞，宋易敏

(1. 中國人民解放軍63611 部隊，庫爾勒841001； 2. 中國人民解放軍63610 部隊，庫爾勒841001)

1 引言

當太陽光線受到地球遮擋后會在地球背后形成陰影區域，稱為地影。 許多繞地運行的衛星都要周期性地進出地影。 衛星位于地影區時，由于星體受到的太陽輻射大幅減少或消失，不僅影響衛星服務功能的提供[1-3]，而且限制航天發射和測控任務的實施[4-7]。 因此，準確、快速計算衛星進出地影的時刻及陰影時間對制定衛星溫控策略、確定太陽輻射壓大小、規劃航天任務等至關重要[8-11]。

在工程應用中主要有2 種地影模型：柱形地影模型和錐形地影模型[12]。 柱形地影模型比較簡單，只能描述衛星被完全照亮或遮擋2 種狀態，適用于低軌衛星的陰影預報[1]，高軌衛星建議謹慎采用[2]。 錐形地影模型考慮到太陽被遮掩的程度，構造了相對復雜和真實的陰影結構，導致陰影判斷條件增多、計算量加大[13-14]，但是卻成為各種高精度應用的必然選擇[3]。 目前關于錐形地影模型陰影判定依據，一種是利用衛星的角度得出太陽、地球的視半徑、角距和視面積來判斷[2]；另一種是根據衛星到陰影軸的距離與影錐在基本面上的投影半徑進行判斷[4]。 由于衛星時刻處于運動之中，根據衛星位置確定的天體視半徑、投影基本面也隨之不斷改變，所以在迭代判定衛星陰影狀態的過程中需要更新每一時刻的陰影條件相關參數。 現有文獻沒有對錐形地影模型不同陰影計算方法的預報精度、計算效率進行定量分析。

本文借鑒賈向華[1]論述的陰影特征角和判定角的概念，把變化率很小的影錐半錐角作為比較的參考基準，提出了一種新的錐形陰影判定條件。 只對背向太陽側的衛星軌道進行地影計算，并給出了基于時間窗口數據篩選用于近圓軌道衛星的地影快速計算方法；使用STK 和MATLAB 工具對3 種軌道高度的衛星進行地影仿真計算，將不同算法計算結果與STK 結果進行比較和分析，以驗證本文算法的精度和效率。

2 錐形地影模型及其判定條件

將太陽視為有一定大小的發光球體，忽略地球扁率和大氣的影響，衛星為一個不計尺寸的點目標，太陽光以直線形式傳播，這樣因地球遮擋在空間形成的圓錐形陰影如圖1 所示。 在地球的背后形成了太陽光完全照射不到的區域和受到部分太陽光照射的區域，前者稱為本影，后者稱為半影。 以日心Os和地心Oe連線為軸，整個地影由本影圓錐和半影圓錐同軸互相嵌套構成。V1為半影圓錐的頂點，α為半影錐的半錐角；V2為本影圓錐的頂點，β為本影錐的半錐角。 本文結合圓錐形陰影幾何特征定義陰影特征角和判定角，根據空間星體幾何關系推導出它們的求解方法和衛星進出地影的判定條件。

圖1 錐形地影模型Fig.1 Conical Earth shadow model

定義rs、re、rt分別為太陽、地球和衛星在慣性坐標系(如J2000.0 坐標系)中的位置矢量，則太陽、衛星相對于地球的坐標s、st分別為式(1)、式(2)：

Rs、Re分別為太陽和地球的半徑，由相似三角形和比例的性質可以計算得到式(3)、式(4)：

半影判定角是影錐頂點V1和衛星位置Ot的連線與陰影軸的夾角，角的開口方向與α一致，記作θ；本影判定角是影錐頂點V2和衛星位置Ot的連線與陰影軸所夾的非鈍角，記作φ。s1、s2分別為地心Oe指向影錐頂點V1、V2的矢量，φ為輔助變量。 以上夾角和矢量由式(5)～式(9)確定：

在地球一個公轉周期中，地球與太陽之間的距離在1.47×108km 到1.52×108km 之間變化[8]，則α的變化區間為0.264 8°～0.273 8°，β的變化區間為0.260 0°～0.268 8°，α的變化率為5.707 8×10-10°/s、β的變化率為5.580 9×10-10°/s，可見α、β隨時間變化很小，在短期內可以視為常值。 由圖1可以看出，地球的本影區實際上內嵌在半影圓錐空間之中，所以確定衛星是否位于地球的陰影內，首先要判斷衛星是否在外層圓錐范圍內，然后再判斷衛星是否在內層圓錐區域內。下面給出地影的判定條件，如式(10)所示：

式中，θ＝α或者φ＝β時，表示衛星正好處于進/出半影或者本影位置。

3 衛星地影快速計算方法

太陽、地球和衛星在空間沿各自的軌道運動，只有在滿足特定的位置關系時，衛星才可能進入到地影中。 因此，實際上僅有一部分衛星軌道穿越地球陰影區，顯然對所有軌道數據進行計算不是合理的做法。 本文提出一種數據篩選方法，先剔除衛星軌道中不可能位于陰影區的數據，然后根據第2 節提出的錐形陰影判定條件對剩余的軌道數據進行計算，從而大幅度減少數據量，以獲得較高的計算效率。

太陽、地球和衛星的空間幾何關系可以用太陽地心矢量s和衛星地心矢量st之間的夾角σ來描述，如圖2 所示。 衛星從地球向陽的一側朝著地球背陽一側運動，當衛星位于位置O′t時(圖1)，σ為90°，衛星與半影圓錐母線交于點M，M到地心的距離為d，其計算公式如式(11)所示：

圖2 太陽、地球和衛星的空間位置關系Fig.2 Relative position of the Sun，the Earth and the satellite

圖3 夾角σ 與衛星地影的關系Fig.3 Relationship between angle σ and Earth shadow of the satellite

把上述星體幾何約束關系轉換為對應的衛星軌道數據時間窗口，可以剔除不可能發生地影的時段。 如圖2 所示，以地球為中心，太陽在黃道面上運動，衛星在其軌道面上運動，太陽地心矢量s在衛星軌道面上的投影與衛星地心矢量st所成的夾角η和σ的變化趨勢一致。 根據幾何中的三垂線定理，當σ＝90°時，η＝90°。 這樣可以用夾角η在同一個軌道平面中等價表征太陽、地球和衛星的位置關系，并且方便運用天體運動規律計算衛星從某一時刻初始位置到達特定位置所需經歷的時間。

為方便計算，太陽、衛星軌道數據時間采用約簡儒略日表示。 地心衛星軌道坐標系定義為：坐標原點位于地心，x軸在軌道面內且指向近地點，y軸在軌道面內且與x軸垂直，指向由x軸逆時針旋轉90°確定，z軸與x軸、y軸構成右手螺旋系。 地心赤道坐標系到地心衛星軌道坐標系的轉換如式(12)所示：

式中Xe為地心赤道坐標系矢量，Xo為地心衛星軌道坐標系矢量，RX、RZ分別為繞x、z軸旋轉的算子，T 表示轉置矩陣[8]，Ω是升交點赤經，i是傾角，ω是近地點輻角。

目前大多數衛星都采用近圓軌道[8，15]。 這樣可以近似認為衛星作勻速圓周運動，同時把太陽也作勻速運動處理。 太陽和衛星在圓形軌道上的角速度分別如式(13)、(14)所示：

式中365.25 代表地球公轉周期，地球引力系數GM＝398 600.441 5 km3/s2，a為衛星軌道半長軸。 在衛星軌道坐標系中，使用式(15)計算太陽投影和衛星的真近點角ν，以確定某一時刻太陽、衛星在軌道面的相對位置及其夾角。

式中x、y分別為太陽或衛星在地心衛星軌道平面中的位置坐標。 假設衛星從某一位置運行到另一位置繞地心旋轉了Δν弧度，則所需時間t為式(16)：

當衛星使用逆行軌道，即與太陽運動方向相反時，式分母取＋號。

綜上，依據初始時刻太陽、衛星的相對位置關系可以計算出衛星在之后周期運動中可能會穿越地影的時間窗口，并對這些時間段內的軌道數據進行陰影判定，此時陰影判定條件中的第2 個不等式可以忽略。

4 仿真計算與分析

在適當的近似假設下，基于錐形地影模型，利用衛星處的太陽、地球的視半徑和視面積判定衛星陰影狀態的方法為天體視半徑判定法[2，4，8](算法1)；利用半影錐、本影錐在基本面上的投影半徑及衛星到陰影軸的距離進行衛星地影解算的方法為影錐基本面投影半徑判定法[4，8](算法2)。 這2 種算法簡單實用，輸入高精度太陽和衛星星歷，能夠計算得到滿足工程應用要求的解。 定義本文提出的基于影錐特征角的陰影判定條件和快速計算方法為影錐特征角快速判定法(算法3)。 使用MATLAB R2015b軟件編寫計算機程序對衛星進出地影的時間進行預報和分析。

選擇低軌、中軌和高軌3 顆近圓軌道衛星進行地影時間計算和分析，衛星數據如表1 所示。首先使用STK11.0.1 計算2019 年7 月30 日0 時(UTC)至2019 年8 月2 日0 時(UTC)期間3 顆衛星的地影狀態，然后利用其高精度軌道外推模型獲取該時間范圍內間隔為1 s 的太陽和衛星J2000.0 坐標系位置矢量數據。

STK(Satellite Tool Kit)是一個經過航天任務驗證的軟件，在仿真分析方面具有極高的準確性[16]。 表2 為3 種算法計算得出的衛星繞地球飛行一圈，經歷一次半影、本影的平均持續時間，并以STK 地影計算結果為基準，計算每種算法的相對誤差。 數據表明3 種算法地影時間計算精度基本相當，半影時間最大誤差5.447%，本影時間最大誤差0.972%，陰影時間越短、計算的相對誤差也越大。

表1 衛星數據Table 1 Satellite data

表2 衛星穿越地影的平均持續時間及相對誤差Table 2 Average duration of satellite through penumbra/umbra and relative errors s

由圖1 可知，衛星在周期運動中按照半影-本影-半影的順序進出地影。 以衛星第一次進入半影的時刻和最后一次離開半影的時刻分別作為衛星進地影的時刻和出地影的時刻，將3 種算法的計算結果和STK 預報的地影時刻進行比較得到衛星進出地影時間差，圖4～圖6 分別是3 顆衛星在仿真計算期間進出地影時間差的變化。 從表2和圖4～圖6 可以看出，算法1 和算法2 對3 種衛星的計算結果完全相同；對于低軌、中軌、高軌衛星，算法3與算法1、算法2 預報的衛星進/出地影時刻的最大時間差分別為1 s、2 s、4 s，且算法3 預報的地影時刻準確度稍高于算法1、算法2。

在處理相同的數據量下，比較3 種算法的計算效率。 計算機的軟硬件配置為Intel Core i7-9700 處理器、32 GB 內存、Windows 7 Service Pack 1 64 位操作系統。 表3 記錄了輸入3 組時間間隔為1 s、共計259 201 點太陽和衛星位置矢量數據時，MATLAB R2015b 軟件測量得到的每一種算法的平均完成時間。 由表3 可知，3 種算法的計算效率高低依次為：算法3、算法1 和算法2，算法3 比算法1 效率提高25.305%，比算法2 提高52.172%，算法3 具有非常明顯的時間性能優勢。

圖4 32060 衛星進出地影時間差值變化Fig.4 Changes of time difference of satellite 32060 entry/exit the Earth shadow

圖5 29061 衛星進出地影時間差值變化Fig.5 Changes of time difference of satellite 29061 entry/exit the Earth shadow

圖6 26880 衛星進出地影時間差值變化Fig.6 Changes of time difference of satellite 26880 entry/exit the Earth shadow

表3 3 種算法計算時間的比較Table 3 Comparison of time calculated by three methods s

由仿真結果可以看出，算法3 地影判定依據正確，預報時間窗口設計合理，在略優于現有其他算法預報精度的條件下，具有較高的計算效率，適用于航天工程中大數據量、高強度地影時間計算。

5 結論

1)新的錐形地影條件以影錐半錐角為特征角和判定依據，在迭代過程中按時間段計算特征角，可大幅減少條件中參數更新的次數；根據發生地影時太陽、地球和衛星存在的幾何邊界條件，確定與這些特定空間位置對應的軌道數據的時間窗口，可大大減少參與地影計算的數據量。

2)仿真計算結果表明，本文建立的地影快速計算方法預報精度略優于天體視半徑判定法和影錐基本面投影半徑判定法，并且計算效率分別提高25.305%和52.172%。

3)盡管影錐特征角快速判定法以地球作為遮擋物進行分析和計算，但該算法同樣適用于其他天體(如月球)遮擋的情況。 下一步應研究橢圓軌道幾何特征和運動規律，進一步拓寬地影快速計算方法應用范圍。