基于影錐特征角的近圓軌道衛(wèi)星進(jìn)出地影時(shí)間的快速判定法

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

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

1 引言

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

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

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

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

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

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

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

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

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

在地球一個(gè)公轉(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，可見α、β隨時(shí)間變化很小，在短期內(nèi)可以視為常值。 由圖1可以看出，地球的本影區(qū)實(shí)際上內(nèi)嵌在半影圓錐空間之中，所以確定衛(wèi)星是否位于地球的陰影內(nèi)，首先要判斷衛(wèi)星是否在外層圓錐范圍內(nèi)，然后再判斷衛(wèi)星是否在內(nèi)層圓錐區(qū)域內(nèi)。下面給出地影的判定條件，如式(10)所示：

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

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

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

太陽、地球和衛(wèi)星的空間幾何關(guān)系可以用太陽地心矢量s和衛(wèi)星地心矢量st之間的夾角σ來描述，如圖2 所示。 衛(wèi)星從地球向陽的一側(cè)朝著地球背陽一側(cè)運(yùn)動(dòng)，當(dāng)衛(wèi)星位于位置O′t時(shí)(圖1)，σ為90°，衛(wèi)星與半影圓錐母線交于點(diǎn)M，M到地心的距離為d，其計(jì)算公式如式(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)換為對(duì)應(yīng)的衛(wèi)星軌道數(shù)據(jù)時(shí)間窗口，可以剔除不可能發(fā)生地影的時(shí)段。 如圖2 所示，以地球?yàn)橹行?，太陽在黃道面上運(yùn)動(dòng)，衛(wèi)星在其軌道面上運(yùn)動(dòng)，太陽地心矢量s在衛(wèi)星軌道面上的投影與衛(wèi)星地心矢量st所成的夾角η和σ的變化趨勢(shì)一致。 根據(jù)幾何中的三垂線定理，當(dāng)σ＝90°時(shí)，η＝90°。 這樣可以用夾角η在同一個(gè)軌道平面中等價(jià)表征太陽、地球和衛(wèi)星的位置關(guān)系，并且方便運(yùn)用天體運(yùn)動(dòng)規(guī)律計(jì)算衛(wèi)星從某一時(shí)刻初始位置到達(dá)特定位置所需經(jīng)歷的時(shí)間。

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

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

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

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

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

當(dāng)衛(wèi)星使用逆行軌道，即與太陽運(yùn)動(dòng)方向相反時(shí)，式分母?。?hào)。

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

4 仿真計(jì)算與分析

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

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

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

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

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

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

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

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

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

圖6 26880 衛(wèi)星進(jìn)出地影時(shí)間差值變化Fig.6 Changes of time difference of satellite 26880 entry/exit the Earth shadow

表3 3 種算法計(jì)算時(shí)間的比較Table 3 Comparison of time calculated by three methods s

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

5 結(jié)論

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

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

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