陳余軍，魏 釗，龍 騰，史人赫

（1. 中國(guó)空間技術(shù)研究院通信與導(dǎo)航總體部，北京100094；2. 北京理工大學(xué)宇航學(xué)院，北京100081；3. 飛行器動(dòng)力學(xué)與控制教育部重點(diǎn)實(shí)驗(yàn)室，北京100081；4. 清華大學(xué)航天航空學(xué)院，北京100084）

1 引 言

地球靜止軌道（Geostationary Earth Orbit，GEO）衛(wèi)星在軌運(yùn)行時(shí)會(huì)受到地球非球形、日月引力和太陽(yáng)光壓等攝動(dòng)的影響，使其逐漸偏離標(biāo)稱軌道，因此，GEO 衛(wèi)星在軌運(yùn)行期間需要進(jìn)行位置保持操作［1-3］。傳統(tǒng)的GEO 衛(wèi)星通常采用化學(xué)推力器來(lái)實(shí)現(xiàn)位置保持，然而，其存在燃料消耗較大、位置保持精度較低等問(wèn)題，嚴(yán)重影響了GEO 衛(wèi)星的壽命［4］。相比之下，電推力器具有高比沖、小推力的特點(diǎn)，基于電推力器的位置保持可有效減少燃料消耗并提高位置保持精度，近年來(lái)得到了廣泛的應(yīng)用［5-7］，例如美國(guó)波音公司的ABS 3A 全電推進(jìn)衛(wèi)星、德國(guó)OHB 公司的Electra 全電推進(jìn)衛(wèi)星以及歐洲通信衛(wèi)星公司的172B 全電推進(jìn)衛(wèi)星［8］等。由于電推力器產(chǎn)生的推力大小僅是化學(xué)推力器的百分之一，且位置保持周期較長(zhǎng)，因此，實(shí)現(xiàn)有效的GEO 位置保持并提高位置保持效率仍面臨巨大挑戰(zhàn)。

對(duì)于基于電推進(jìn)系統(tǒng)的GEO衛(wèi)星，其背地板通常安裝有4 臺(tái)電推力器，并形成一個(gè)矩形構(gòu)型［9-10］。然而，當(dāng)其中一臺(tái)電推力器發(fā)生故障時(shí)，其對(duì)角線上另一臺(tái)電推力器單獨(dú)不能產(chǎn)生期望的速度增量來(lái)抵消攝動(dòng)的影響，只有剩下的兩臺(tái)電推力器可以正常用于衛(wèi)星的位置保持，即故障模式［11］。在完備模式下，4 臺(tái)電推力器可以依次在升降交點(diǎn)開(kāi)機(jī)來(lái)修正軌道傾角、偏心率等，能有效實(shí)現(xiàn)GEO 定點(diǎn)位置保持并最大化位置保持效率［12］。然而在故障模式下，2 臺(tái)電推力器依次在升降交點(diǎn)開(kāi)機(jī)不能抵消攝動(dòng)的影響，因此必須針對(duì)電推力器的小推力特性，對(duì)位置保持策略進(jìn)行定制。

目前，有不少學(xué)者開(kāi)展了故障模式下的電推力器位置保持策略研究工作，例如，Zhang等［11］考慮地影等約束，提出了一種故障模式下的最優(yōu)電推力器位置保持控制策略，通過(guò)優(yōu)化點(diǎn)火位置和點(diǎn)火時(shí)長(zhǎng)，可有效降低燃料消耗；Li 等［13］基于故障診斷算法設(shè)計(jì)了一套實(shí)時(shí)閉環(huán)位置保持策略，用于實(shí)現(xiàn)故障模式下GEO衛(wèi)星的位置保持，并通過(guò)仿真驗(yàn)證了該方法的有效性；劉宇鑫［14］提出了三次點(diǎn)火故障模式下的電推力器位置保持策略，通過(guò)仿真驗(yàn)證了該策略可有效實(shí)現(xiàn)故障模式下的GEO 位置保持。盡管已有大量文獻(xiàn)對(duì)故障模式下的GEO 衛(wèi)星位置保持策略進(jìn)行了研究，但大部分研究工作沒(méi)有考慮電推力器布局的影響。由于GEO 衛(wèi)星在位置保持期間，南北位保所消耗的燃料占整個(gè)位保過(guò)程所消耗燃料的80%以上［15］，而電推力器的布局直接決定了產(chǎn)生的推力的方向，因此，對(duì)電推力器布局進(jìn)行合理的設(shè)計(jì)與優(yōu)化，對(duì)于提高電推力沿南北方向的分量系數(shù)，進(jìn)而提高位置保持效率、降低燃料消耗具有非常重要的意義。

本文為進(jìn)一步提高故障模式下基于電推力器的GEO衛(wèi)星的位置保持效率，提出了一種故障模式下電推力器布局設(shè)計(jì)優(yōu)化方案。首先建立了故障模式下電推力器位置保持模型，然后以電推力器安裝位置、開(kāi)機(jī)赤經(jīng)和開(kāi)機(jī)時(shí)長(zhǎng)為設(shè)計(jì)變量，考慮位置保持精度等約束，建立了故障模式下電推力器布局設(shè)計(jì)優(yōu)化模型，通過(guò)優(yōu)化顯著提高了位置保持效率、降低了燃料消耗。

2 故障模式下電推力器位置保持模型

2.1 軌道動(dòng)力學(xué)模型

由于GEO 的軌道傾角和偏心率趨近于0，為防止動(dòng)力學(xué)方程發(fā)生奇異，本文采用非奇異的軌道根數(shù)來(lái)描述衛(wèi)星的運(yùn)動(dòng)，如式（1）所示：

其中，a為半長(zhǎng)軸，Ω為升交點(diǎn)赤經(jīng)，ω為近心點(diǎn)角距，i為軌道傾角，M為平近點(diǎn)角，λ為平經(jīng)度，S為格林尼治恒星時(shí)。

衛(wèi)星在軌期間通過(guò)電推力器開(kāi)機(jī)產(chǎn)生徑向、切向和法向速度增量來(lái)抵消空間攝動(dòng)的影響，電推力器造成的經(jīng)度漂移率改變量、偏心率矢量改變量和軌道傾角矢量改變量可表示為：

其中，ΔVR、ΔVT和ΔVN分別為電推力器開(kāi)機(jī)產(chǎn)生的徑向、切向和法向速度增量，ΔD為經(jīng)度漂移率改變量，Δex和Δey為偏心率矢量改變量，Δix和Δiy為軌道傾角矢量改變量，as為GEO 半長(zhǎng)軸，Vs為GEO 衛(wèi)星運(yùn)行速度，l為電推力器開(kāi)機(jī)赤經(jīng)。

GEO 衛(wèi)星在軌運(yùn)動(dòng)時(shí)，主要受到日月引力、地球非球形引力和太陽(yáng)光壓攝動(dòng)的影響。

由日月引力造成的軌道傾角矢量改變量如式（3）所示［16］：

其中，Ωsm為白道升交點(diǎn)黃經(jīng)。

由地球非球形引力造成的經(jīng)度改變量如式（4）所示［16］：

其中，J22、λ22、J31、λ31、J33和λ33為地球田諧項(xiàng)系數(shù)；Re是地球赤道半徑；r是衛(wèi)星與地心之間的距離；ωe是地球自轉(zhuǎn)角速度；TD是一天的秒數(shù)。

由太陽(yáng)光壓造成的偏心率矢量改變量如式（5）所示［16］：

其中，K是太陽(yáng)光壓系數(shù)；A是衛(wèi)星正對(duì)太陽(yáng)光的面積；m是衛(wèi)星質(zhì)量；ε是黃赤交角；T0是初始儒略世紀(jì)數(shù)；Te是最終的儒略世紀(jì)數(shù)；Ls為歷元時(shí)刻太陽(yáng)平黃經(jīng)，可通過(guò)式（6）計(jì)算得到：

其中，JD為儒略日。

基于脈沖假設(shè)，電推力器產(chǎn)生的速度增量和燃料消耗可通過(guò)式（7）計(jì)算得到：

其中，Δv為電推力器產(chǎn)生的速度增量，Δmfuel為燃料消耗，Isp為電推力器比沖，g為地球引力系數(shù)，m0為衛(wèi)星初始質(zhì)量，F(xiàn)p為電推力器推力，Δt為開(kāi)機(jī)時(shí)長(zhǎng)。

2.2 小推力位置保持控制策略

衛(wèi)星背地板共裝有4 臺(tái)電推力器，其中電推力器1和2安裝在背地板北側(cè)，電推力器3和4安裝在背地板南側(cè)，對(duì)稱安裝。在故障模式下，當(dāng)其中一臺(tái)電推力器發(fā)生故障時(shí)，則該電推力器與對(duì)角線上的電推力器不能正常用于衛(wèi)星位置保持，以電推力器2 發(fā)生故障為例，則電推力器1 和4 可以正常用于衛(wèi)星位置保持，而電推力器2 和3 不能正常用于衛(wèi)星位置保持，如圖1所示。

圖1 故障模式下電推力器布局Fig.1 Configuration of electric thrusters in failure mode

衛(wèi)星正常在軌運(yùn)行時(shí)，電推力器產(chǎn)生的推力應(yīng)指向衛(wèi)星質(zhì)心，則推力沿徑向R、切向T和法向N的分量系數(shù)為：

其中，kR，kT和kN分別為推力沿徑向、切向和法向的分量系數(shù)，dR，dT和dN分別為電推力器安裝位置沿軌道坐標(biāo)系R，T，N三軸坐標(biāo)的絕對(duì)值。

圖2 電推力器開(kāi)機(jī)赤經(jīng)示意圖Fig.2 Electric thruster firing right ascension diagram

一個(gè)完整的位置保持周期是14天，由1個(gè)為期2 天的軌道確定周期和6 個(gè)為期2 天的小控制周期組成。其中，在軌道確定周期內(nèi)測(cè)定衛(wèi)星當(dāng)前的軌道根數(shù)，在小控制周期內(nèi)實(shí)現(xiàn)衛(wèi)星的GEO定點(diǎn)位置保持。在完備模式下，4 臺(tái)電推力器可以依次在升降交點(diǎn)開(kāi)機(jī)以最大化位置保持效率，并有效實(shí)現(xiàn)GEO 定點(diǎn)位置保持。但在故障模式下，有2 臺(tái)電推力器不能正常用于衛(wèi)星位置保持，剩余2 臺(tái)電推力器依次在升降交點(diǎn)開(kāi)機(jī)不能有效抵消攝動(dòng)的影響，因此必須找到電推力器合適的開(kāi)機(jī)位置。本研究采用文獻(xiàn)［11］中提出的位置保持控制策略，即在每個(gè)小控制周期內(nèi)，電推力器1和4各開(kāi)機(jī)兩次，電推力器1 開(kāi)機(jī)赤經(jīng)分別在0°～90°和90°～180°區(qū)間內(nèi)，電推力器4開(kāi)機(jī)赤經(jīng)分別在180°～270°和270°～360°區(qū)間內(nèi)，電推力器開(kāi)機(jī)赤經(jīng)如圖2 所示，通過(guò)優(yōu)化電推力器的開(kāi)機(jī)赤經(jīng)和開(kāi)機(jī)時(shí)長(zhǎng)，來(lái)實(shí)現(xiàn)GEO定點(diǎn)位置保持。電推力器在上述的赤經(jīng)范圍內(nèi)開(kāi)機(jī)，把產(chǎn)生的速度增量代入式（2），可得電推力器開(kāi)機(jī)造成的經(jīng)度漂移率改變量如式（9）所示：

其中，ΔDt為電推力器開(kāi)機(jī)造成的經(jīng)度漂移率改變量，kT1和kT4分別為電推力器1 和電推力器4 的推力切向分量系數(shù)，ΔvN1、ΔvN2、ΔvS1和ΔvS2分別為電推力器1 第一次開(kāi)機(jī)、電推力器1 第二次開(kāi)機(jī)、電推力器4第一次開(kāi)機(jī)和電推力器4 第二次開(kāi)機(jī)產(chǎn)生的速度增量。

2015年12月在宣布了中國(guó)承辦2016年G20峰會(huì)之后，習(xí)總書(shū)記在土耳其是這樣向世界介紹杭州的：“杭州是歷史文化名城，也是創(chuàng)新活力之城，相信2016年峰會(huì)將會(huì)給大家呈現(xiàn)一種歷史和現(xiàn)實(shí)交匯的獨(dú)特韻味?！薄昂贾菁瘸錆M濃郁的中華文化韻味，也擁有面向世界的寬廣視野?！?/p>

電推力器開(kāi)機(jī)造成的軌道傾角矢量改變量如式（10）所示：

其中，Δixt和Δiyt為電推力器開(kāi)機(jī)造成的軌道傾角矢量改變量；kN1和kN4分別為電推力器1 和電推力器4的推力法向分量系數(shù)；lN1、lN2、lS1和lS2分別為電推力器1第一次開(kāi)機(jī)、電推力器1第二次開(kāi)機(jī)、電推力器4第一次開(kāi)機(jī)和電推力器4第二次開(kāi)機(jī)的赤經(jīng)。

電推力器開(kāi)機(jī)造成的偏心率矢量改變量如式（11）所示：

其中，Δext和Δeyt為電推力器開(kāi)機(jī)造成的偏心率矢量改變量；kR1和kR4分別為電推力器1 和電推力器4 的推力徑向分量系數(shù)。

為有效實(shí)現(xiàn)位置保持，在每個(gè)小控制周期內(nèi)，電推力器產(chǎn)生的軌道根數(shù)改變量應(yīng)等于攝動(dòng)造成的軌道根數(shù)改變量，如式（12）所示：

其中，ΔDd為攝動(dòng)造成的經(jīng)度漂移率改變量；Δixd和Δiyd為攝動(dòng)造成的軌道傾角矢量改變量；Δexd和Δeyd為攝動(dòng)造成的偏心率矢量改變量。

按照式（12），使電推力器逐次開(kāi)機(jī)來(lái)抵消攝動(dòng)的影響，即可實(shí)現(xiàn)GEO定點(diǎn)位置保持。

3 故障模式下電推力器布局設(shè)計(jì)優(yōu)化問(wèn)題描述

3.1 電推力器布局設(shè)計(jì)優(yōu)化模型

故障模式下電推力器布局設(shè)計(jì)優(yōu)化問(wèn)題主要包括兩部分，一是每個(gè)小控制周期的位保參數(shù)優(yōu)化，二是電推力器安裝位置優(yōu)化。其中，電推力器安裝位置優(yōu)化在外層，每個(gè)小控制周期的位保參數(shù)優(yōu)化在內(nèi)層。

其中，tN1、tN2、tS1和tS2分別為電推力器1 第一次開(kāi)機(jī)、電推力器1 第二次開(kāi)機(jī)、電推力器4 第一次開(kāi)機(jī)和電推力器4 第二次開(kāi)機(jī)的時(shí)長(zhǎng)，mfuel為每個(gè)小控制周期的燃料消耗；λf和φf(shuō)分別為每個(gè)小控制周期末端的經(jīng)度和緯度；λd為定點(diǎn)經(jīng)度。

電推力器安裝位置優(yōu)化問(wèn)題本質(zhì)上是在保證每個(gè)小控制周期的優(yōu)化過(guò)程均能找到可行解的前提下，探索最優(yōu)的電推力器1 和電推力器4 的安裝位置，使得整個(gè)位置保持周期的燃料消耗最少。電推力器安裝位置優(yōu)化問(wèn)題的數(shù)學(xué)模型如式（14）所示：

其中，T1和N1分別為電推力器1 的切向和法向安裝坐標(biāo)；mfuel_all為整個(gè)位置保持周期的總?cè)剂舷?；g1用于判斷每個(gè)小控制周期的位保參數(shù)優(yōu)化是否找到可行解，若沒(méi)有找到可行解，則g1＞0，反之，g1≤0。

3.2 電推力器布局設(shè)計(jì)優(yōu)化流程

故障模式下電推力器布局設(shè)計(jì)優(yōu)化流程圖如圖3 所示，其中，為提高全局搜索性，外層的電推力器安裝位置通過(guò)遺傳算法（Genetic Algorithm，GA）進(jìn)行優(yōu)化；為提高找到可行解的速度，每個(gè)小控制周期的位保參數(shù)通過(guò)多初始點(diǎn)序列二次規(guī)劃算法（Sequential Quadratic Programming，SQP）進(jìn)行優(yōu)化，具體步驟如下所示。

步驟一：通過(guò)GA優(yōu)化電推力器安裝位置，計(jì)算推力沿徑向、切向和法向分量系數(shù)。

步驟二：在軌道確定周期內(nèi)，確定衛(wèi)星初始軌道根數(shù)以及攝動(dòng)對(duì)軌道根數(shù)的改變量。

步驟三：確定每一個(gè)小控制周期內(nèi)所需的電推力器開(kāi)機(jī)產(chǎn)生的軌道根數(shù)改變量，通過(guò)多初始點(diǎn)SQP 優(yōu)化電推力器的開(kāi)機(jī)赤經(jīng)和開(kāi)機(jī)時(shí)長(zhǎng)，由式（7）計(jì)算速度增量以及燃料消耗量。若優(yōu)化沒(méi)有找到可行解，則返回步驟一；否則，循環(huán)步驟二和步驟三，直到達(dá)到最大位置保持天數(shù)。

步驟四：輸出最優(yōu)電推力器安裝位置，以及每個(gè)小控制周期內(nèi)的電推力器開(kāi)機(jī)赤經(jīng)和開(kāi)機(jī)時(shí)長(zhǎng)。

圖3 故障模式下電推力器布局設(shè)計(jì)優(yōu)化流程圖Fig.3 Electric thruster configuration design optimization flowchart in failure mode

4 仿真案例

本節(jié)通過(guò)某故障模式下的GEO 衛(wèi)星電推力器布局設(shè)計(jì)優(yōu)化實(shí)例來(lái)驗(yàn)證提出方案的有效性。本案例中，GEO 衛(wèi)星入軌質(zhì)量為2000kg，衛(wèi)星尺寸為2m×2m×2m，衛(wèi)星面質(zhì)比為0.05，光壓系數(shù)為1.5，電推力器額定推力為100mN，比沖4000s，仿真開(kāi)始時(shí)間為2020年12月1日，仿真時(shí)長(zhǎng)為1年，GEO 衛(wèi)星定點(diǎn)經(jīng)度為120°E。采用MATLAB R2017b ga 和fmincon 工具包進(jìn)行優(yōu)化，其中，GA 種群規(guī)模設(shè)為200，最大迭代次數(shù)設(shè)為30，SQP 初始點(diǎn)個(gè)數(shù)設(shè)為50，其他參數(shù)設(shè)為默認(rèn)。衛(wèi)星初始軌道根數(shù)如表1所示，設(shè)計(jì)變量取值范圍如表2所示。

表1 衛(wèi)星初始軌道根數(shù)Table 1 Initial orbital elements of the satellite

表2 設(shè)計(jì)變量取值范圍Table 2 Range of design variables

通過(guò)優(yōu)化，可以得到電推力器開(kāi)機(jī)赤經(jīng)和速度增量示意圖如圖4 所示，經(jīng)度曲線和緯度曲線如圖5 和圖6 所示，優(yōu)化前后電推力器的安裝位置對(duì)比如表3 所示，優(yōu)化前后電推力器布局示意圖如圖7所示。

圖4 電推力器開(kāi)機(jī)赤經(jīng)和速度增量Fig.4 Firing right ascensions and velocity increments of electric thrusters

圖5 經(jīng)度曲線Fig.5 Longitude curve

圖6 緯度曲線Fig.6 Latitude curve

表3 電推力器的安裝位置對(duì)比Table 3 Electric thruster configuration comparison

表4 初始方案與優(yōu)化方案燃料消耗對(duì)比Table 4 Fuel consumption comparison between the initial scheme and optimized schemes

從圖4 可以看出，電推力器1 第二次開(kāi)機(jī)和電推力器4 第一次開(kāi)機(jī)產(chǎn)生的速度增量相對(duì)較大，因?yàn)閺膱D2 可以看出，只有電推力器1 第二次開(kāi)機(jī)和電推力器4 第一次開(kāi)機(jī)能產(chǎn)生向右的推力，而電推力器1和4的四次開(kāi)機(jī)均能產(chǎn)生向左的推力，因此，電推力器1第二次開(kāi)機(jī)和電推力器4第一次開(kāi)機(jī)時(shí)間相對(duì)較長(zhǎng)，來(lái)抵消四次開(kāi)機(jī)產(chǎn)生的反向推力的影響。從圖5 和圖6 可以看出，GEO 衛(wèi)星的經(jīng)緯度均控制在［-0.05°，0.05°］，滿足控制精度要求。從表3可以計(jì)算得到優(yōu)化前后的推力法向分量系數(shù)分別為0.43 和0.69，且從計(jì)算結(jié)果和圖7 可以看出，優(yōu)化后電推力器的法向坐標(biāo)到達(dá)邊界，切向坐標(biāo)的絕對(duì)值減小，使得優(yōu)化布局的推力法向分量系數(shù)相比于初始布局增加了60.47%。因?yàn)槟媳蔽槐Ｋ牡娜剂险伎側(cè)剂系?0%以上，因此，優(yōu)化后的電推力器布局更有利于提高位置保持效率。

優(yōu)化前后燃料消耗對(duì)比如表4 所示，從表4 可以看出基于優(yōu)化布局的一年的位置保持燃料消耗相比于初始布局減少了3.93kg（約為41.9%），驗(yàn)證了本文提出的故障模式下電推力器布局設(shè)計(jì)優(yōu)化方案的有效性與工程實(shí)用性。

為進(jìn)一步說(shuō)明提出的故障模式下電推力器布局設(shè)計(jì)優(yōu)化方案的優(yōu)勢(shì)，基于上述優(yōu)化框架，對(duì)完備模式下的電推力進(jìn)行布局設(shè)計(jì)優(yōu)化，一年的位置保持燃料消耗對(duì)比如表5所示。

表5 基于故障和完備模式下電推力器布局的燃料消耗對(duì)比Table 5 Fuel consumption comparison between the failure mode and nominal mode based electric thruster configurations

從表5 可以看出，盡管只有兩臺(tái)電推力器可以用于位置保持，但基于故障模式下電推力器布局的一年位置保持燃料消耗相比于完備模式僅增加了0.46kg，因此說(shuō)明了提出的故障模式下電推力器布局設(shè)計(jì)優(yōu)化方案具有較強(qiáng)的在軌應(yīng)急能力。

5 結(jié) 論

本文提出了一種故障模式下電推力器布局設(shè)計(jì)優(yōu)化方案，可以有效實(shí)現(xiàn)故障模式下GEO衛(wèi)星的定點(diǎn)位置保持，并顯著提高位置保持效率。首先針對(duì)電推力器的小推力特性，建立了故障模式下電推力器位置保持模型。然后考慮電推力器安裝位置、開(kāi)機(jī)赤經(jīng)、開(kāi)機(jī)時(shí)長(zhǎng)以及位置保持精度等，建立了故障模式下電推力器布局設(shè)計(jì)優(yōu)化模型，并通過(guò)GA 與SQP 算法進(jìn)行求解。優(yōu)化結(jié)果表明，基于優(yōu)化布局的位置保持燃料消耗相比于初始布局減少了41.9%，因此說(shuō)明了提出的故障模式下電推力器布局設(shè)計(jì)優(yōu)化方案的有效性與工程實(shí)用性。