基于軌道可達域的機動航天器接近威脅規避方法

2024-04-07 02:23:20張賽楊震杜向南羅亞中

航空學報 2024年4期

張賽，楊震，*，杜向南，羅亞中

1.國防科技大學空天科學學院，長沙 410073

2.空天任務智能規劃與仿真湖南省重點實驗室，長沙 410073

3.上海機電工程研究所，上海 201109

隨著技術的進步，在軌航天器數量急劇增加，并且具有了更加強大的機動能力。一方面，軌道機動使得軌道攔截、交會、飛越等空間任務成為可能；另一方面，在軌航天器面臨被非合作目標機動接近甚至打擊的風險。已有的不確定性事件［1］和X-37B、MEV 等空間實驗［2-3］都表明太空安全形式已經變得非常復雜，因此開展來襲航天器威脅評估與規避方法的研究十分必要。

考慮實際工程約束，為了評估空間任務的可行性和安全性，有必要研究機動能力有限的非合作目標軌道機動可達域計算方法。航天器軌道機動可達域（Reachable Domain， RD）是在給定軌道機動約束條件下航天器所有可能到達位置的集合［4］，是機動航天器進行態勢分析的有效工具。Caruso 等［5］首先將可達域的概念應用于彈道武器的可達性分析，Vinh 等［6］和Battin［7］隨后將此概念擴展到彈道武器的攔截問題和具有固定速度的航天器軌道轉移問題。在此基礎上，Xue 等［8］研究了具有固定脈沖幅值的航天器單脈沖機動（三維）可達域，但求解精度有待提高。Wen 等［9］提出的將可達域包絡求解問題轉化為可達方向上矢徑極值求解問題的方法，能夠精確求解機動位置固定的單脈沖機動可達域。該方法被杜向南和楊震［10］所采用并作了進一步優化，求解流程更加簡單。可達域方法已經在航天器地面軌跡覆蓋區域分析［11］、目標飛越［12-13］、航天器軌道安全性分析與碰撞預警［14-15］以及軌道博弈［16-17］等空間任務的分析中得到應用。在軌道安全性分析中，可達域能夠直觀描述航天器軌道在非合作目標可達域范圍內的危險區段，并得到航天器穿過危險區段時間［15］，即得到目標的威脅告警信息。

基于可達域獲得來襲目標的潛在威脅告警信息后，航天器可以進行變軌以規避目標威脅。已有的航天器威脅規避方法主要分為針對無機動目標的碰撞規避策略和針對來襲機動目標的主動規避策略2 類［18］。對于無機動空間目標的碰撞規避，現有研究以碰撞概率計算為基礎［19-20］，研究最優機動的計算方法和實施策略等問題［21］。相比于無機動空間目標，可機動空間目標產生的打擊威脅對航天器來說更加危險且難以預測。目前對于航天器機動規避策略的研究較為充分［22-25］。在軌航天器對來襲機動航天器進行機動規避時，兩航天器均主動施加控制但追求相反的相對位置狀態，因此該過程屬于軌道博弈的范疇［26］。部分學者研究了追逃博弈雙邊規劃問題［27-30］，分別給出了博弈雙方的追逃策略，上述文獻研究對于航天器威脅規避策略的形成具有一定的參考價值。然而在空間攻防中，攻擊方大多具備較高機動能力，地面站難以及時更新其異常機動信息，使得在軌航天器不能及時探測到被打擊的威脅。因此如何通過有限信息分析來襲目標潛在威脅能力并設計主動規避機動策略對潛在威脅域進行避讓具有一定的研究意義。

本文在現有研究基礎上，提出了基于軌道可達域的機動航天器接近威脅計算、評估和規避方法。首先給出了航天器單脈沖軌道機動可達域的一般求解方法；其次，通過判斷在軌航天器軌道在來襲航天器可達域內的區段得到危險區域，并基于兩航天器進出危險區的時間定義危險區比例和時間窗口匹配性威脅評價指標；最后，給出了基于最小化危險區的航天器多脈沖機動主動規避策略，并通過仿真驗證了所提基于軌道可達域的機動航天器接近威脅規避方法的有效性。

1 單脈沖軌道機動可達域求解

1.1 坐標系定義及軌道狀態描述

1.1.1 坐標系定義

如圖1所示，定義J2000 地心慣性坐標系SJ(OJxJyJzJ)：坐標原點OJ在地球中心；xJ軸沿地球赤道面和黃道面的交線，指向春分點?；zJ軸和地球自轉軸重合并指向北極；yJ軸根據右手法則得到。定義近心點坐標系Sp(Opxpyp)zp：原點Op位于地心；xp軸指向初始軌道近地點；zp軸指向初始軌道角動量方向；yp軸位于初始軌道平面內，并與xp、zp軸構成右手直角坐標系。

圖1 坐標系位置關系示意圖Fig.1 Illustration of defined coordinates

如圖2所示，定義航天器LVLH 坐標系：原點O位于航天器質心；x軸與r重合；z軸指向軌道角動量方向；y軸位于軌道平面內，并與x、z軸構成右手直角坐標系。初始軌道上的LVLH 坐標系記為S0(O0x0y0z0)，轉移軌道上的LVLH 坐標系記為S1(O1x1y1z1)。

圖2 坐標系關系與速度機動球示意圖Fig.2 Illustration of defined coordinates and impulsive maneuvering ball

1.1.2 軌道狀態描述

如圖2所示，將航天器初始位置矢量和終端位置矢量分別記為r0和rf；定義航天器初始軌道平面為M0；轉移軌道平面為M1；記M0和M1的夾角為β，r0和rf的夾角為θ。

設航天器在初始軌道上機動位置處的軌道根數為[a0e0i0Ω0ω0f0]，則r0在坐標系Sp下可表示為

初始速度v0在坐標系S0中可表示為

式中：μ=3.986×1014m3s2為地球引力常數；為軌道半通徑。

初始速度v0在坐標系S1中可表示為［4］

將航天器的最大機動能力記為Δvmax，機動方向任意，因此其機動速度增量可以用速度機動球建模表示。若航天器的機動速度增量為，其在轉移軌道面內的投影為ΔvM，其大小為［4］

1.2 單脈沖軌道機動可達域求解方法

1.2.1 可達域問題描述

航天器單脈沖軌道機動可達域是圍繞初始軌道的有限三維區域，該區域被一個邊界面所覆蓋，稱為可達域包絡［15］。通常情況下，二維平面上的包絡由內包絡和外包絡構成，可達域及其包絡的幾何關系如圖3所示，求解可達域的核心工作是確定空間中全部可達區域的包絡。

圖3 初始軌道、可達域及其包絡的幾何關系Fig.3 Geometry of initial orbit， RD， and envelope of RD

為簡化計算，可以將連續包絡的求解離散為若干包絡點的求解。在任意一個以地心為坐標原點的笛卡爾坐標系中，方位角γ和φ能夠確定空間中的一個方向d(γ，φ)。若d(γ，φ)可達，可求得該可達方向上最近距離dmin和最遠距離dmax，因此A(γ，φ，dmin)、B(γ，φ，dmax)兩點分別為P上內包絡和外包絡的點。通過遍歷參數γ、φ的取值能夠得到所有包絡點，最終得到可達域。

1.2.2 軌道機動可達性判據

將圖3 中的坐標系選作近心點坐標系Sp，為使得Sp中的方向d(γp，φp)可達，需保證機動后的速度矢端能夠落在邊界速度雙曲線（Terminal Velocity Hyperbola， TVH）上。此時，可達的判定條件為［4］

1.2.3 可達域求解算法

已知坐標系Sp中的一個可達方位d(γp，φp)，rf為該方位上任意終端矢徑。航天器從r0以速度機動至rf，則基于二體軌道動力學模型可算得rf的大小為［10］

定義表達式［10］：

由文獻［10］可知，對于機動位置固定的可達域求解問題，rf的方位角[γp，φp]確定后，rf矢徑表達式（7）是關于變量α的一元函數，的零點可以通過求解P(α)=0 得到。

易知P(α) 至多存在2 個零點，記為和，同時將區間兩端邊界值設為和2π。將和分別代入式（7），得到、和。記4 個值中較大的為dmax，較小的為dmin，可得Sp坐標系下rf方向上的極值點位置坐標為

綜上所述，可達域計算步驟為：

步驟1首先給出γp和φp用于確定空間中的一個方位。

步驟2根據式（6）判斷該方位是否可達，若可達進入步驟3；若不可達，返回步驟1。

步驟3通過求解式（8）的零點求得該可達方位矢徑極值，確定極值點坐標，即式（9），即可得到該方位上可達域包絡。

步驟4重復步驟1～步驟3，通過遍歷γp和φp的取值得到所有可達方位上的內外包絡點，進而獲得Sp中的近似包絡，求得可達域。

2 在軌航天器軌道危險區求解及威脅評估

來襲航天器軌道可達域將會是圍繞其初始軌道平面形成的扁平形區域。如圖4所示，將來襲航天器軌道可達域與在軌航天器軌道的交集部分定義為在軌航天器軌道危險區。危險區會出現在來襲航天器初始軌道平面與在軌航天器軌道平面的交點附近，在危險區內，在軌航天器存在被接近的風險。

圖4 在軌航天器軌道危險區示意圖Fig.4 Diagram of danger area for on-orbit spacecraft

2.1 軌道危險區域求解方法

由于來襲航天器機動能力有限，因此要進行威脅規避機動策略設計，首先要判斷危險區的存在性。

設來襲航天器機動前的初始軌道根數為[acecicΩcωcfc]；在軌航天器軌道的初始軌道根數為[atetitΩtωtft]。如圖5所示，綠色軌道面為來襲航天器初始軌道，天球面交點緯度幅角為uc；藍色軌道面為目標軌道，天球面緯度幅角為ut。

圖5 軌道面天球交點Fig.5 Orbital celestial sphere node

在慣性坐標系SJ下，兩軌道天球面的交點在慣性系中存在等式：

式中：rs表示天球的半徑；M1(·)表示繞OJxJ軸旋轉的轉移矩陣；M3(·)表示繞OJzJ軸旋轉的轉移矩陣。

進一步可以求得

從式（11）中可以求解得到uc和ut的值，且由圖5 可知，uc和ut有2 組解關于引力中心對稱，且是一一對應的，即

uc1和ut1對應的交點在來襲航天器初始軌道及在軌航天器軌道上的真近點角fc1、ft1表示為

若來襲航天器在fc1處采取軌道機動，則兩航天器軌道在天球面的交點的方向在來襲航天器近心點坐標系中可表示為

若來襲航天器軌道可達域與在軌航天器軌道存在交點，如圖6所示，虛線表示軌道面交線，綠色軌跡表示目標軌道，在右側在軌航天器軌道存在危險點如圖中紅點所示。

圖6 在軌航天器軌道危險點Fig.6 Danger point of on-orbit spacecraft

由于來襲航天器機動可達域具備一定厚度，因此在軌航天器軌道穿過該威脅域的危險區段應是圖6所示危險點附近一段延伸的軌道區段。

二體假設下航天器任意時刻的位置速度可通過初始時刻t0=0 時的位置速度r0、v0外推得到。設外推時間步長為Δt，ti時刻表示為

通過二體軌道外推確定的ti時刻在軌航天器位置ri［7］在基于來襲航天器初始軌道建立的近心點坐標系Sp中可以表示為

進一步，在來襲航天器近心點坐標系中確定(ri)p的空間指向[γi，φi]為

將式（17）所得ti時刻在軌航天器的空間指向[γi，φi]代入1.2.3 節中可獲得來襲航天器在該指向上的可達范圍[ri，min，ri，max]，此時在軌航天器的地心距為，分析目標地心距與來襲航天器可達域之間的關系可能有如下2 種情況：

按照式（15）給定的時間步長對目標軌道外推并記錄所有處于危險區的時間范圍[tmin，tmax]，如圖7所示。

圖7 在軌航天器軌道危險區的時間范圍Fig.7 Time window of on-orbit spacecraft’s danger area

圖8 來襲航天器與危險區關系Fig.8 Relationship between incoming maneuvering spacecraft and danger area

由于在軌航天器軌道是否處于危險區是按照步長Δt外推進行遍歷搜索的，因此圖7所示危險區的2 個端點應分別在[tmin-Δt，tmin] 和[tmax，tmax+Δt]時間內獲得。危險區端點時刻的精確區間可以通過二分法求得。以[tmin-Δt，tmin]區間內的端點為例，具體步驟為：

步驟1根據式（16）確定t=(tmin+tmin-Δt)2 處目標軌道地心距(ri)p∈[ri，min，ri，max]。

步驟2若區間長度ε＞ε0（一般可取為0.001s），則進入步驟3；否則結束計算，所得區間窗口記為危險區時間窗口。

步驟3判斷該地心距是否位于危險區內，若是，則用t替換tmin，組成新的時間范圍；若不是，則用t替換tmin-Δt，組成新的時間范圍；重復步驟1 和步驟2。

對危險區的2 個端點時刻分別執行上述步驟從而得到滿足高精度要求的危險區時間窗口[t'min，t'max]。

2.2 威脅評價指標

當在軌航天器的軌道存在危險區時，其受威脅程度取決于兩航天器軌道關系及來襲航天器的機動能力。下面給出指標描述在軌航天器受來襲航天器的威脅程度。

2.2.1 基于危險區比例的評價指標

由圖 7 可知，可采用危險區時間占軌道周期的比例作為評價指標。設2.1 節所述方法判定危險區存在，且經過二分法求得在軌航天器處于危險區的時間范圍為[t'min，t'max]，則危險系數ξD表示為

式中：危險系數ξD越大，在軌航天器軌道位于威脅可達域內的比例越高，被來襲航天器打擊的可能性就越大。

2.2.2 考慮時間窗口匹配性的評價指標

除了考慮上述危險系數ξD外，還需要考慮來襲航天器到達危險區時間和在軌航天器本身到達危險區時間窗口的匹配性。下面給出來襲航天器到達危險區的時間窗口估計方法。

圖 8所示藍色平面為危險區的2 個端點垂直于來襲航天器初始軌道平面在威脅域內形成的2 個截面。2 個截面在基于來襲航天器初始軌道建立的近心點坐標系Sp中，所對應的方位角γ分別為γ1、γ2。危險區2 個端點的位置分別為rtmin、rtmax。當來襲航天器進行機動后，轉移軌道的平均角速度可以表示為

式中：atran為轉移軌道的半長軸。

由于來襲航天器采取的脈沖機動限制在速度機動球內，其到達危險區的時間窗口可以根據軌道轉移時間取值范圍進行估計。由式（19）可知，在轉移相位相同的情況下，轉移時間僅與轉移軌道的半長軸atran相關。又由于當來襲航天器的最大機動能力Δvmax全部用于進行面內機動時，可以使atran的變化幅度最大，因此本文采用來襲航天器進行軌道面內機動到達γ1、γ2截面的時間極值作為到達危險區時間窗口的保守估計。由于γ1截面距離出發點更近，對應時間窗口估計的下界；γ2截面距離出發點更遠，對應時間窗口估計的上界。

設來襲航天器在fw0處脈沖機動大小為Δvmax，方向與原速度方向相反，則由式（2）可知，機動后速度為

由rw1和vw1可算得機動后轉移軌道的半長軸aw1，偏心率ew1以及機動位置真近角fw1，則截面γ1處對應的真近角可以表示為

由真近角與偏近點角轉換關系，能夠得到來襲機動平臺在機動位置和到達γ1截面處對應的偏近點角Ew1、Eγ1。進而可得按照最大軌道平均角速度到達危險區的最小時間為

同理為了得到來襲航天器到達危險區的時間窗口上界，在fw0處施加脈沖大小Δvmax，方向與原速度方向相同。

與求解窗口下界方法類似，可以求得來襲機動平臺在機動位置和到達截面γ2處的偏近角Ew2、Eγ2，則按照最小軌道平均角速度到達危險區的最大時間為

式（22）和式（24）確定的來襲航天器到達危險區的時間窗口比來襲航天器實際到達危險區具體危險點的時間范圍更廣，可以用于判斷來襲航天器與在軌航天器時間窗口的匹配性，從而評價危險區有效性。

至此能夠得到在軌航天器處于危險區的時間范圍[t'min，t'max]和來襲航天器到達危險區的時間窗口[tw，min，tw，max]。記tp為2 個時間窗口的交集大小，定義危險系數ξT描述來襲機動平臺與在軌航天器的時間窗口匹配性

3 接近威脅規避策略

接近威脅規避策略即在軌航天器的脈沖機動方案，其目的是實現對來襲航天器的威脅規避。本文通過優化方法求解滿足約束條件的燃料最省威脅規避策略，最小化危險區以實現威脅規避。進一步的，考慮到在軌航天器執行任務需求，要求在軌航天器在完成威脅規避的同時能夠回到原軌道，求解最優脈沖機動方案。

設來襲航天器最大機動能力（最大速度脈沖）為Δvmax，c，在軌航天器不采取任何機動時存在至少一段危險軌道區段，且危險軌道區段的危險系數為ξD，時間窗口匹配性為ξT。

在軌航天器采用多脈沖機動變軌對危險區進行避讓，相鄰脈沖之間最小時間間隔為δt。進行機動避讓后的在軌航天器軌道危險系數為ξ'D，時間窗口匹配性為ξ'T。

3.1 優化變量

多脈沖機動規避策略優化變量取為脈沖施加時刻ti(i=1，2，…，n)和施加脈沖矢量Δvi(i=1，2，…，n)，在在軌航天器LVLH 坐標系中，Δvi可以表述為

式中：n為脈沖次數，優化變量共4n個。

假設在軌航天器在機動軌道上通過兩脈沖機動與原軌道交會，則2 次變軌的速度增量可通過Lambert 算法求得。

3.2 約束條件

本文只考慮在軌航天器一個軌道周期T=內的威脅規避機動策略求解。對于t＞T時的規避機動，可基于T時刻的軌道狀態重新求解軌道危險系數和時間窗口匹配性，再進行優化求解。

由于在軌航天器機動能力受限，因此脈沖施加時刻ti滿足

式中：ti表示相對初始時刻脈沖時間。

由于危險區的威脅程度不僅與危險區比例ξD有關，其實際攔截可行性還與時間窗口匹配性ξT相關。為實現對來襲航天器的威脅規避，將危險系數和時間窗口匹配性設為約束條件：

若在軌航天器在機動軌道上通過兩脈沖變軌回到原軌道，則在式（27）和式（28）的基礎上需增加軌道瞄準約束，即最后一次脈沖后在軌航天器軌道根數滿足

式中：Ei(tn)為在軌航天器進行i次機動后在tn時刻的軌道根數。

3.3 目標函數

優化目標為最小化燃料消耗，即

3.4 優化模型和算法

接近威脅規避策略優化模型為

選擇差分進化（Differential Evolution，DE）算法［31］作為多脈沖機動主動規避策略的優化方法。

4 仿真分析

根據軌道危險區和威脅評價指標的定義可知，在軌航天器軌道危險區的存在性及受威脅程度由兩航天器軌道關系及來襲航天器機動能力決定，因此本文提出的方法適用于任意兩航天器軌道為圓軌道或橢圓軌道的情況。下面給出仿真結果以驗證該方法的有效性。

4.1 在軌航天器軌道危險區仿真算例

本小節給出在軌航天器軌道危險區仿真結果，分別考慮在軌航天器與來襲航天器不在同軌道面和處于同軌道面的情況。

設來襲航天器初始軌道根數如表1所示，航天器在當前位置采取脈沖機動，機動量最大值=300 m s，機動方向任意，所得單脈沖機動可達域仿真結果如圖9所示。圖9 中，深綠色網格區域為來襲航天器機動可達域。

表1 來襲航天器初始軌道根數Table 1 Initial orbital elements of chaser

圖9 在軌航天器軌道與來襲航天器可達域位置關系Fig.9 Geometry of chaser’s RD and target’s orbit

設在軌航天器初始軌道根數如表2所示，其中在軌航天器2 與來襲航天器處于相同軌道面上。

表2 在軌航天器初始軌道根數Table 2 Initial orbital elements of target

通過二體軌道外推，可得在軌航天器1 軌道如圖9 和圖10 中淺綠色軌道所示。圖10 中，在軌航天器1 的2 段危險區的位置分別如圖中紅色區段和藍色區段所示。

圖10 在軌航天器1 軌道危險區Fig.10 Orbital danger-zone of Target 1

可知在軌航天器1 經過來襲航天器機動可達域且存在2 段危險區，其進出危險區的時間窗口和來襲航天器進出危險區的時間窗口估計見表3。基于危險區比例的危險系數ξD和時間窗口匹配性系數ξT計算結果見表4。

表3 到達危險區時間窗口Table 3 Time-window of danger-zone

表4 危險系數ξD 和時間窗口匹配性系數ξTTable 4 Danger index ξD and time-window matching index ξT

分析圖10 可知，2 段危險區在來襲機動平臺可達域內的分布具備一定對稱性，間隔約半個軌道周期。經計算來襲航天器初始軌道的運行周期為Tv=9 952.01 s，結合表3 數據可知，來襲航天器2 次到達危險區的時間間隔為7 008.32 s-2 423.35 s=4 584.97 s，驗證了所提危險區求解算法的正確性。

圖11所示為在軌航天器2 的軌道危險區。此時，兩航天器進出危險區的時間窗口及基于危險區比例的危險系數和時間窗口匹配性系數計算結果分別如表5 和表6所示。

表5 到達危險區時間窗口（共面情況）Table 5 Time-window of danger-zone（coplane case）

表6 危險系數ξD 和時間窗口匹配性系數ξT（共面情況）Table 6 Danger index ξD and time-window matching index ξT（coplane case）

圖11 在軌航天器2 軌道危險區Fig.11 Orbital danger-zone of Target 2

由表5 和表6 可知，由于在軌航天器2 與來襲航天器處于相同的軌道平面，其軌道危險區范圍遠大于在軌航天器1 的軌道危險區范圍。

4.2 多脈沖機動主動規避仿真算例

設來襲航天器初始軌道根數如表1所示，在軌航天器初始軌道根數如表2所示。來襲航天器在t=0 s 時刻采取脈沖機動，脈沖機動最大值=300 m s，機動方向任意。經計算，優化前在軌航天器存在2 段危險區，基于危險區比例的危險系數ξD=0.023 613，時間窗口匹配性系數ξT=1。

設在軌航天器相鄰脈沖最小時間間隔Δt=10 s。按照上述設置，可求解在軌航天器多脈沖機動威脅規避策略。

4.2.1 不考慮與原軌道的交會約束

不考慮與原軌道交會的約束時，基于式（31）的優化模型可以求得多脈沖機動威脅規避優化策略，結果如表7所示。以兩脈沖機動規避策略為例，優化后的威脅規避機動軌道如圖12 和圖13所示。

表7 多脈沖機動優化結果（不考慮交會約束）Table 7 Optimal results of multi-impulses （without rendezvous constraint）

圖12 在軌航天器機動軌道與來襲航天器可達域位置關系（不考慮交會約束）Fig.12 Geometry of chaser’s RD and target’s maneuvering orbit （without rendezvous constraint）

圖13 Y-Z 平面在軌航天器規避軌道與來襲航天器可達域位置關系（不考慮交會約束）Fig.13 Geometry of chaser’s RD and target’s maneuvering orbit in Y-Z plane（without rendezvous constraint）

從圖12 和圖13 中可以看出，不考慮交會約束時，在軌航天器威脅規避機動軌道與來襲航天器可達域不存在重合區段，即危險系數ξ'D=0。對比表7 中多脈沖機動規避策略求解結果，有，即當脈沖次數＞2 時，在軌航天器進行威脅規避所需速度增量相差不大，因此可采取兩脈沖機動作為威脅規避方案。

此外，通過對比多個機動方案中所求速度增量發現，其垂直軌道面方向的速度分量遠小于軌道面內的速度分量，這是由于軌道危險區在在軌航天器軌道面內導致的。

4.2.2 考慮與原軌道的交會約束

在4.2.1 節的基礎上，考慮與原軌道交會的約束，多脈沖機動規避策略求解結果如表8所示。以四脈沖機動為例，優化后的轉移軌道如圖14 和圖15所示。從圖14 和圖15 中可以看出，考慮交會約束時，在軌航天器威脅規避機動軌道與來襲航天器可達域仍存在重合區段，但雙方不滿足時間匹配關系，即危險系數ξ'D≠0、時間窗口匹配性系數ξ'T=0。