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

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

1 引言

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

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

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

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

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

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

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

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

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

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

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

3 衛(wèi)星地影快速計算方法

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

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

圖2 太陽、地球和衛(wèi)星的空間位置關(guān)系Fig.2 Relative position of the Sun，the Earth and the satellite

圖3 夾角σ 與衛(wèi)星地影的關(guān)系Fig.3 Relationship between angle σ and Earth shadow of the satellite

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

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

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

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

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

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

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

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

4 仿真計算與分析

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

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

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

表1 衛(wèi)星數(shù)據(jù)Table 1 Satellite data

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

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

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

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

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

圖6 26880 衛(wèi)星進出地影時間差值變化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

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

5 結(jié)論

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

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

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