張錦繡，張 諭，王繼河，楊繼坤，逯振坤，宋昱岐

（1.中山大學(xué) 航空航天學(xué)院，深圳 518107；2.中山大學(xué) 物理與天文學(xué)院，珠海 519082）

引 言

為印證廣義相對(duì)論，引力波探測(cè)成為當(dāng)前國(guó)際物理和空間科學(xué)領(lǐng)域的研究熱點(diǎn)。國(guó)內(nèi)外各研究機(jī)構(gòu)圍繞空間引力波探測(cè)技術(shù)展開了深入研究，如：歐洲航天局（European Space Agency，ESA）的“激光干涉測(cè)量空間天線”（Laser Interferometer Space Antenna，LISA）/eLISA（evolved LISA）計(jì)劃[1-2]，中國(guó)的“天琴計(jì)劃”[3]、“太極計(jì)劃”[4]。

上述引力波探測(cè)計(jì)劃大同小異，如圖1所示，它們的空間引力波探測(cè)器均由3個(gè)航天器組成等邊三角形構(gòu)型，C位于構(gòu)型平面的中心，表示地球或者虛擬的點(diǎn)，航天器圍繞太陽周期性運(yùn)動(dòng)的同時(shí)也圍繞C周期性運(yùn)動(dòng)。另外，它們都假設(shè)空間引力波會(huì)造成航天器搭載的檢驗(yàn)質(zhì)量相對(duì)位置發(fā)生改變，并通過測(cè)量該變化量來印證引力波的存在。檢驗(yàn)質(zhì)量之間的相對(duì)位置，是通過安裝在航天器上的激光干涉測(cè)距系統(tǒng)測(cè)量得到的。為了保證高精度激光干涉測(cè)距系統(tǒng)可以長(zhǎng)時(shí)間處于穩(wěn)定對(duì)準(zhǔn)測(cè)量狀態(tài)，航天器編隊(duì)需要運(yùn)行在可以保持長(zhǎng)期穩(wěn)定構(gòu)型的軌道上。不同的是，LISA/eLISA計(jì)劃和“太極計(jì)劃”，構(gòu)型的中心C為虛擬點(diǎn)；“天琴計(jì)劃”是以地球?yàn)橹行奶祗w，構(gòu)型中心C為地球。

圖1 空間引力波探測(cè)器構(gòu)型示意圖Fig.1 Schematic diagram of space gravitational wave detector configuration

空間引力波探測(cè)編隊(duì)構(gòu)型的穩(wěn)定性主要由3個(gè)指標(biāo)表征：臂長(zhǎng)（任意兩個(gè)探測(cè)器之間的距離）、呼吸角（每?jī)蓚€(gè)臂之間的角度）和星間視線相對(duì)速度。臂長(zhǎng)決定了干涉儀的探測(cè)頻帶，不等臂長(zhǎng)將會(huì)在激光干涉測(cè)量中引入激光頻率噪聲，影響核心科學(xué)載荷的設(shè)計(jì)；呼吸角的變化影響探測(cè)器的指向精度；星間視線相對(duì)速度影響超穩(wěn)時(shí)鐘和相位計(jì)等重要載荷的指標(biāo)設(shè)計(jì)。探測(cè)器為了成功探測(cè)到引力波，需要將3個(gè)航天器之間的臂長(zhǎng)、呼吸角以及相對(duì)速度等三角形軌道構(gòu)型參數(shù)長(zhǎng)時(shí)間保持在一定范圍內(nèi)。準(zhǔn)確的航天器動(dòng)力學(xué)模型和軌道優(yōu)化算法是設(shè)計(jì)滿足上述軌道構(gòu)型的兩大關(guān)鍵要素。

對(duì)于如LISA/eLISA計(jì)劃和“太極計(jì)劃”這種軌道構(gòu)型中心為虛擬點(diǎn)的航天器編隊(duì)方案，多數(shù)學(xué)者以笛卡爾坐標(biāo)或軌道要素為狀態(tài)量，在日心坐標(biāo)系下推導(dǎo)航天器間相對(duì)動(dòng)力學(xué)方程，然后將其簡(jiǎn)化為線性的C-W方程進(jìn)行軌道設(shè)計(jì)。但是，當(dāng)探測(cè)器軌道構(gòu)型中心有中心天體時(shí)（例如“天琴計(jì)劃”），上述方式不再適用。這是因?yàn)楹教炱髦g的臂長(zhǎng)和航天器到中心天體之間的距離在同一量級(jí)，不滿足C-W方程的簡(jiǎn)化條件。針對(duì)有中心天體的航天器編隊(duì)情況，本文以對(duì)偶四元數(shù)為工具，在J2000日心坐標(biāo)系下建立航天器之間的相對(duì)運(yùn)動(dòng)動(dòng)力學(xué)模型。值得注意的是，本文所提方法同樣適用于軌道構(gòu)型中心為虛擬點(diǎn)的航天器編隊(duì)方案。

為了得到滿足空間引力波探測(cè)任務(wù)要求的高精度軌道構(gòu)型，通常需要采用優(yōu)化算法對(duì)航天器運(yùn)行軌道進(jìn)行優(yōu)化。文獻(xiàn)[5]使用全局優(yōu)化的禁忌搜索算法對(duì)LISA引力波探測(cè)軌道進(jìn)行了優(yōu)化設(shè)計(jì)，價(jià)值函數(shù)有18個(gè)自變量，找到了滿足引力波探測(cè)要求的軌道，并給出了從地球停泊軌道進(jìn)入引力波探測(cè)實(shí)驗(yàn)軌道的發(fā)射段和分離段的軌道設(shè)計(jì)。文獻(xiàn)[6]分析了影響LISA編隊(duì)構(gòu)型的性能指標(biāo)（呼吸角、滯后角、相對(duì)速度）相對(duì)于初始軌道要素的敏感性，將LISA軌道優(yōu)化問題轉(zhuǎn)化為兩步級(jí)聯(lián)單目標(biāo)優(yōu)化問題，先優(yōu)化滯后角，后優(yōu)化呼吸角，并分析了軌道周期偏差對(duì)構(gòu)型穩(wěn)定性的影響。針對(duì)“天琴計(jì)劃”軌道構(gòu)型優(yōu)化設(shè)計(jì)，文獻(xiàn)[7]對(duì)構(gòu)型發(fā)生長(zhǎng)周期漂移的軌道要素進(jìn)行了分析，采用粒子群算法對(duì)3顆航天器的軌道半長(zhǎng)軸進(jìn)行優(yōu)化，初步滿足軌道構(gòu)型的穩(wěn)定性指標(biāo)。文獻(xiàn)[8]采用組合優(yōu)化方法，討論了由于軌道進(jìn)動(dòng)而導(dǎo)致探測(cè)器指向的緩慢長(zhǎng)期漂移，設(shè)計(jì)了滿足構(gòu)型穩(wěn)定性要求的任務(wù)軌道。

雖然上述研究滿足了引力波探測(cè)構(gòu)型穩(wěn)定性的需求，但是它們均是在已知初始軌道的基礎(chǔ)上，對(duì)航天器的軌道根數(shù)進(jìn)行優(yōu)化，并且軌道優(yōu)化的起始時(shí)刻是確定的。沒有考慮初始軌道以外是否存在更加穩(wěn)定的軌道，使得航天器滿足引力波探測(cè)構(gòu)型穩(wěn)定性的需求。運(yùn)行在甚高軌道的航天器會(huì)受到地球、月球和太陽等主要天體的引力影響，若在軌道優(yōu)化過程中通過星歷得到每一時(shí)刻地球、月球和太陽的位置，這嚴(yán)重影響了優(yōu)化算法的計(jì)算效率。因此，有必要以太陽為中心建立慣性坐標(biāo)系，為地球和月球建立簡(jiǎn)化但足夠精確的動(dòng)力學(xué)模型，從而加快優(yōu)化速度。把地球、月球和航天器看作剛體來描述它們的運(yùn)動(dòng)規(guī)律，這實(shí)際上是描述六自由的運(yùn)動(dòng)問題。而對(duì)偶四元數(shù)[9-13]是目前描述剛體運(yùn)動(dòng)的最合適的數(shù)學(xué)工具，這是因?yàn)樗鼘⑵揭坪托D(zhuǎn)統(tǒng)一在同一框架下，同時(shí)考慮了二者之間的耦合影響，并且可以提供很高的計(jì)算效率。

綜上所述，本文在J2000日心坐標(biāo)系下，研究日–地–月引力場(chǎng)中，空間引力波探測(cè)系統(tǒng)穩(wěn)定構(gòu)型的全局搜索優(yōu)化方案。①基于對(duì)偶四元數(shù)，在J2000日心坐標(biāo)系下建立地球、月球和3個(gè)航天器的姿軌一體化動(dòng)力學(xué)模型；②分析構(gòu)型穩(wěn)定性的表征，確定優(yōu)化的目標(biāo)函數(shù)和待優(yōu)化的變量；③使用遺傳算法對(duì)航天器軌道構(gòu)型進(jìn)行優(yōu)化設(shè)計(jì)，使得地球、月球和3顆航天器在不進(jìn)行軌道控制的前提下盡可能長(zhǎng)時(shí)間地保持穩(wěn)定的構(gòu)型；④通過數(shù)學(xué)仿真驗(yàn)證優(yōu)化方法的有效性和普適性。

1 日心系編隊(duì)動(dòng)力學(xué)表征模型

1.1 自然/人造天體動(dòng)力學(xué)方程

考慮日地月引力場(chǎng)對(duì)引力波探測(cè)器的影響，本文在日心系下，將自然天體（地球、月球）和人造天體（3個(gè)航天器）視為一個(gè)完整的編隊(duì)。這樣可以在建立編隊(duì)動(dòng)力學(xué)模型的時(shí)候，更加方便地對(duì)日地月引力進(jìn)行表述。

地球在J2000日心坐標(biāo)系[14]OsXsYsZs下的運(yùn)動(dòng)學(xué)方程為

地球在J2000日心坐標(biāo)系OsXsYsZs下的動(dòng)力學(xué)方程為

其中：=fe+ετe為作用于地球質(zhì)心的對(duì)偶力旋量；=medI/dε+εJe表示地球的對(duì)偶慣量矩陣；me表示地球的質(zhì)量；Je表示地球的轉(zhuǎn)動(dòng)慣量。

月球在J2000日心坐標(biāo)系OsXsYsZs下的運(yùn)動(dòng)學(xué)方程為

月球在J2000日心坐標(biāo)系OsXsYsZs下的動(dòng)力學(xué)方程為

其中：=fm+ετm為作用于月球質(zhì)心的對(duì)偶力旋量；=mmdI/dε+εJm表示月球的對(duì)偶慣量矩陣；md表示月球的質(zhì)量；Jd表示月球的轉(zhuǎn)動(dòng)慣量。

以航天器的質(zhì)心為坐標(biāo)原點(diǎn)，建立本體坐標(biāo)系Oi?xiyizi，3個(gè)坐標(biāo)軸方向分別與航天器的慣量主軸重合。航天器在日心坐標(biāo)系下運(yùn)動(dòng)學(xué)方程表示為[15]

其中：i=1,2,3分別表示航天器SC1、SC2和SC3；對(duì)偶四元數(shù)qi表示第i個(gè)航天器的本體坐標(biāo)系相對(duì)于日心坐標(biāo)系的四元數(shù)，rii表示航天器的位置矢量在本體坐標(biāo)系下的表示；對(duì)偶速度旋量=wii+εvii，wii表示航天器相對(duì)于太陽的角速度在本體坐標(biāo)系下的表示，表示航天器相對(duì)于太陽的速度矢量在本體坐標(biāo)系下的表示。

動(dòng)力學(xué)方程為

其中：=midI/dε+εJi表示對(duì)偶慣量矩陣，mi表示航天器的質(zhì)量，Ji表示航天器的轉(zhuǎn)動(dòng)慣量；表示航天器受到的對(duì)偶力。當(dāng)圖1中C點(diǎn)為虛擬點(diǎn)時(shí)，航天器以太陽為中心天體，若只考慮太陽引力，那么航天器所受到的對(duì)偶力為fsi表示太陽的引力。當(dāng)圖1中C點(diǎn)以地球?yàn)橹行奶祗w時(shí)，航天器受到地球、月球和太陽的引力和太陽光壓，那么航天器所受到的對(duì)偶力為：分別表示地球引力、月球引力、太陽引力、太陽光壓和重力梯度力矩。

其中：地心引力常數(shù)為μe=398 600.441 90 km3/s2，月心引力常數(shù)為μm=4 902.800 076 km3/s2；日心引力常數(shù)為μs=132 712 440 040.944 00 km3/s2；太陽光壓為P⊙表示太陽光壓強(qiáng)，它與航天器到太陽的距離有關(guān)，A表示航天器表面積，r⊙表示太陽相對(duì)于地心的位置矢量，AU表示地球與太陽的距離，δ 表示衛(wèi)星表面材料的反射系數(shù)。

1.2 編隊(duì)成員間相對(duì)動(dòng)力學(xué)方程

本文主要研究航天器之間的相對(duì)運(yùn)動(dòng)，即兩兩航天器本體坐標(biāo)系Oi?xiyizi和Oj?xjyjzj（i,j=1,2,3,i≠j）之間的相對(duì)運(yùn)動(dòng)。利用乘性誤差對(duì)偶四元數(shù)，該相對(duì)運(yùn)動(dòng)可表示為

坐標(biāo)系Oi?xiyizi相對(duì)于坐標(biāo)系Oj?xjyjzj的運(yùn)動(dòng)學(xué)方程為

航天器之間相對(duì)動(dòng)力學(xué)方程為

本研究關(guān)注的核心指標(biāo)是探測(cè)器之間的相對(duì)運(yùn)動(dòng)參數(shù)，但由于編隊(duì)中包含了自然天體，所以需要建立探測(cè)器和臨近自然天體之間的相對(duì)運(yùn)動(dòng)關(guān)系。

航天器相對(duì)于臨近天體的相對(duì)運(yùn)動(dòng)方程為

航天器相對(duì)于臨近天體的相對(duì)動(dòng)力學(xué)方程為

2 穩(wěn)定構(gòu)型優(yōu)化設(shè)計(jì)

2.1 待優(yōu)化參數(shù)

本文將進(jìn)入預(yù)定軌道的起始時(shí)刻d作為優(yōu)化參數(shù)之一，尋找使得航天器構(gòu)型最穩(wěn)定的時(shí)刻。

在日心坐標(biāo)系下，確定了軌道優(yōu)化的起始時(shí)刻后，此時(shí)地球、月球的初始位置和速度均可通過星歷DE421獲得，通過式（1）～（4）可得到初始時(shí)刻之后地球和月球在日心坐標(biāo)系下的位置和速度。

航天器本體坐標(biāo)系相對(duì)于日心坐標(biāo)系的四元數(shù)表示為

其中：n為單位軸；σ 為繞單位軸轉(zhuǎn)過的角度；β為方位角；φ為極角。單位軸n=[nx,ny,nz]由極坐標(biāo)形式表示為

其中：β,φ ∈[0,2π]。

在航天器軌道坐標(biāo)系中，位置矢量為rii=[r,0,0]T，其中r=p/(1+ecosf)，p=a(1?e2)，a為軌道半長(zhǎng)軸，e為偏心率，f為真近點(diǎn)角。速度矢量為vii=[vr,vt,0]T，為徑向速度,(1+ecosf)為橫向速度。所以當(dāng)航天器本體坐標(biāo)系和質(zhì)心軌道坐標(biāo)系重合，每個(gè)航天器需要6個(gè)參數(shù)p、σ、β、φ、e、f作為待優(yōu)化參數(shù)，優(yōu)化得到航天器的軌道。

2.2 目標(biāo)函數(shù)的確定

根據(jù)表征等邊三角形穩(wěn)定性要素指標(biāo)要求，目標(biāo)函數(shù)CF由衛(wèi)星之間的臂長(zhǎng)變化量CF1和呼吸角變化量CF2兩部分線性組成，并求最小值。

其中：w1、w2為權(quán)重；c1、c2為歸一化常數(shù)；θ0=60?；L為衛(wèi)星之間的臂長(zhǎng)。

干涉臂長(zhǎng)和呼吸角的計(jì)算公式為

其中：i,j,k=1,2,3，i≠j≠k，rij(t)、rik(t)、rjk(t) 為兩航天器之間相對(duì)位置向量。

2.3 遺傳算法優(yōu)化方法

根據(jù)前面的推導(dǎo)結(jié)果可以看出，對(duì)于空間引力波探測(cè)系統(tǒng)，只要確定了構(gòu)型的優(yōu)化參數(shù)，就可以計(jì)算出地球、月球和航天器的位置和速度。取不同的優(yōu)化參數(shù)，就可以得到不同的構(gòu)型，構(gòu)型穩(wěn)定性問題可以抽象為求取目標(biāo)函數(shù)CF(t)全局最小值的問題。

對(duì)于引力波探測(cè)器軌道的優(yōu)化，目標(biāo)函數(shù)有19個(gè)自變量，存在數(shù)量眾多的局部極小值。使用簡(jiǎn)單的局部?jī)?yōu)化算法易陷入局部極小值點(diǎn)而無法找到滿足要求的軌道，所以這里選擇遺傳算法作為全局優(yōu)化算法。遺傳算法只利用目標(biāo)函數(shù)的信息，無需梯度等高價(jià)信息，適用于大規(guī)模、高度非線性的不連續(xù)的目標(biāo)函數(shù)的優(yōu)化，具有很強(qiáng)的通用性，是構(gòu)型優(yōu)化設(shè)計(jì)算法的一個(gè)理想選擇。

2.3.1 編碼和適應(yīng)度函數(shù)

遺傳算法的編碼方式包括二進(jìn)制編碼、浮點(diǎn)數(shù)編碼等，本文采用二進(jìn)制編碼方式，通過控制二進(jìn)制的位數(shù)決定參數(shù)的精度。

選擇合適的適應(yīng)度函數(shù)是設(shè)計(jì)遺傳算法過程中的一個(gè)關(guān)鍵問題。本文以式（17）作為適應(yīng)度函數(shù)

2.3.2 遺傳算子

1）選擇算子

通過選擇算子選出種群中較優(yōu)的個(gè)體形成一個(gè)新的種群，可使得種群中的個(gè)體不斷逼近最優(yōu)解。選擇算子可采用常用的輪盤算法，設(shè)定個(gè)體被選中的概率如下

其中：f′為個(gè)體j的適應(yīng)度值；n為種群數(shù)目。在上述選擇方式下，個(gè)體的適應(yīng)度值越大，被選擇的概率也越大，則其基因構(gòu)造被遺傳到下一代的可能性越大。這種選擇方式使得種群中的個(gè)體都存在被選中的概率，保證了較小個(gè)體的遺傳性。

2）交叉算子

交叉是指按交叉概率從種群中選擇部分個(gè)體，通過交換個(gè)體部分基因形成新的個(gè)體，從而形成新的種群。交叉算子可選擇單點(diǎn)交叉算子。

3）變異算子

變異算子模擬的是生物進(jìn)化過程中的基因突變現(xiàn)象，具體操作為：令變異概率為Pc，則變異的節(jié)點(diǎn)個(gè)數(shù)為[nSPc]，S為單個(gè)個(gè)體的節(jié)點(diǎn)數(shù)；根據(jù)生成的隨機(jī)數(shù)r決定變異方向，選中某些個(gè)體上的某些節(jié)點(diǎn)，然后用新值替代這些節(jié)點(diǎn)上的原值。變異算子維持了群體的多樣性，有利于防止早熟現(xiàn)象的出現(xiàn)。

2.3.3 算法流程

本文算法的計(jì)算步驟如下：

1）生成初始種群；

2）根據(jù)設(shè)計(jì)的適應(yīng)度函數(shù)，計(jì)算該代種群所有個(gè)體的適應(yīng)度；

3）判斷是否滿足終止條件，若不滿足則到步驟4，若滿足則到步驟5；

4）根據(jù)交叉概率和變異概率，經(jīng)過選擇算子、交叉算子和變異算子產(chǎn)生下一代種群，然后返回步驟3；

5）選出最后一代種群中的最優(yōu)個(gè)體，輸出相應(yīng)的最優(yōu)解，算法終止。

3 仿真算例與結(jié)果分析

為了驗(yàn)證所提方法的普適性，以地心軌道和日心軌道空間引力波探測(cè)任務(wù)為應(yīng)用背景，對(duì)構(gòu)型優(yōu)化問題進(jìn)行仿真驗(yàn)證。

3.1 地心軌道空間引力波探測(cè)系統(tǒng)構(gòu)型優(yōu)化仿真

以地心軌道空間引力波探測(cè)系統(tǒng)為應(yīng)用背景，對(duì)所提方法進(jìn)行仿真驗(yàn)證。3個(gè)航天器各自圍繞地球飛行，軌道高度為10萬km，假設(shè)軌道偏心率和真近點(diǎn)角均為0，航天器軌道由4個(gè)參數(shù)p、σ、β、φ唯一確定，參數(shù)范圍如表1所示。軌道起始日期為2034年1月1日，在1年的時(shí)間內(nèi)優(yōu)化得到最優(yōu)的構(gòu)型。遺傳算法的種群規(guī)模為100，交叉概率為0.7，變異概率為0.08，迭代次數(shù)為100代。其余參數(shù)分別為c1=3/104，c2=3/400，w1=0.5，w2=0.5。P⊙=4.56×10?6Nm2，δ=0.21，A=12m2，AU=149 597 870 691km。將優(yōu)化結(jié)果表示為位置速度后如表2所示。

表1 地心軌道編隊(duì)待優(yōu)化參數(shù)及取值范圍Table 1 Geocentric orbit formation parameters to be optimized and their value ranges

表2 地心軌道編隊(duì)優(yōu)化結(jié)果Table 2 Geocentric orbital formations optimization results

圖2是地心軌道編隊(duì)航天器SC1、SC2和SC3之間的干涉臂長(zhǎng)變化率曲線，圖2中綠色曲線表示干涉臂長(zhǎng)L12的變化率，最大值?L12為0.14%；紅色曲線表示干涉臂長(zhǎng)L13的變化率，最大值?L13為0.14%；藍(lán)色曲線表示干涉臂長(zhǎng)L23的變化率，最大值?L23為0.18%。

圖2 地心軌道編隊(duì)航天器之間的干涉臂長(zhǎng)變化率Fig.2 Geocentric orbit formation arm length variation between spacecraft

圖3是地心軌道編隊(duì)航天器SC1、SC2和SC3軌道構(gòu)型呼吸角的變化曲線，圖中綠色曲線表示呼吸角θ1，變化量最大值 ?θ1為0.14°；紅色曲線表示呼吸角θ2，變化量最大值? θ2為0.084°；藍(lán)色曲線表示呼吸角θ3，變化量最大值? θ3為0.13°。

圖3 地心軌道編隊(duì)呼吸角Fig.3 Geocentric orbit formation breathing angle

圖4是地心軌道編隊(duì)航天器SC1、SC2和SC3相對(duì)速度變化曲線，圖中綠色曲線表示SC1和SC2之間的相對(duì)速度，變化量最大值?v1為6.17 m/s；紅色曲線表示SC1和SC3之間的相對(duì)速度，變化量最大值?v2為4.56 m/s；藍(lán)色曲線表示SC2和SC3之間的相對(duì)速度，變化量最大值 ?v3為6.70 m/s。

圖4 地心軌道編隊(duì)航天器之間的相對(duì)速度Fig.4 Geocentric orbit formation relative velocity between spacecraft

上述結(jié)果表明，在考慮了日、地、月引力和太陽光壓攝動(dòng)的情況下，地心空間引力波探測(cè)系統(tǒng)的軌道在一年內(nèi)的變化能滿足引力波探測(cè)的需求。隨著時(shí)間的推移，構(gòu)型逐漸發(fā)散，需要通過軌道機(jī)動(dòng)將航天器編隊(duì)配置恢復(fù)到初始狀態(tài)。

以優(yōu)化得到的結(jié)果作為初始條件，用Monte-Carlo法分析入軌誤差對(duì)地心軌道編隊(duì)構(gòu)型的影響。分別向航天器的徑向、切向和法向位置增加標(biāo)準(zhǔn)差為5 m的位置誤差，對(duì)構(gòu)型的影響如表3所示。分別向航天器的徑向、切向和法向速度增加標(biāo)準(zhǔn)差為0.1 mm/s的速度誤差，對(duì)構(gòu)型的影響如表4所示。位置誤差和速度誤差均符合正態(tài)分布?？梢钥闯觯瑥较蛭恢谜`差和切向速度誤差導(dǎo)致構(gòu)型標(biāo)準(zhǔn)差的變化最大。因此，航天器的入軌誤差中徑向位置誤差和切向速度誤差對(duì)構(gòu)型穩(wěn)定性的影響最大。若位置誤差和速度誤差同時(shí)存在，當(dāng)速度誤差是1 mm/s時(shí)，位置誤差超過5 m 時(shí)構(gòu)型將發(fā)散。

表3 位置誤差對(duì)地心軌道編隊(duì)構(gòu)型的影響Table 3 Effect of position error on configuration of geocentric orbital formation

表4 速度誤差對(duì)地心軌道編隊(duì)構(gòu)型的影響Table 4 Effect of velocity error on configuration of geocentric orbital formation

3.2 日心軌道空間引力波探測(cè)系統(tǒng)構(gòu)型優(yōu)化仿真

以日心軌道空間引力波探測(cè)系統(tǒng)為應(yīng)用背景，對(duì)所提方法進(jìn)行驗(yàn)證。該方案中3顆航天器各自圍繞太陽飛行，所構(gòu)成等邊三角形的臂長(zhǎng)L=2.5×106km，編隊(duì)三角形中心與地球相差約20°。本算例中待優(yōu)化參數(shù)范圍如表5所示，遺傳算法參數(shù)配置與地心軌道編隊(duì)算例相同。將優(yōu)化結(jié)果表示為位置速度后如表6所示。

表5 日心軌道編隊(duì)待優(yōu)化參數(shù)及取值范圍Table 5 Heliocentric orbit formation parameters to be optimized and their value ranges

表6 日心軌道編隊(duì)優(yōu)化結(jié)果Table 6 Heliocentric orbital formations optimization results

圖5是日心軌道編隊(duì)航天器SC1、SC2和SC3之間的干涉臂長(zhǎng)變化率曲線，圖中綠色曲線表示干涉臂長(zhǎng)L12的變化率，最大值?L12為0.48%；紅色曲線表示干涉臂長(zhǎng)L13的變化率，最大值?L13為0.25%；藍(lán)色曲線表示干涉臂長(zhǎng)L23的變化率，最大值?L23為0.32%。

圖5 日心軌道編隊(duì)航天器之間的臂長(zhǎng)變化率Fig.5 Heliocentric orbit formation arm length variation between spacecraft

圖6是日心軌道編隊(duì)航天器SC1、SC2和SC3軌道構(gòu)型呼吸角的變化曲線，圖中綠色曲線表示呼吸角θ1，變化量最大值?θ1為0.36°；紅色曲線表示呼吸角θ2，變化量最大值? θ2為0.38°；藍(lán)色曲線表示呼吸角θ3，變化量最大值? θ3為0.26°。

圖6 日心軌道編隊(duì)航天器呼吸角Fig.6 Heliocentric orbit formation spacecraft breathing angle

圖7是日心軌道編隊(duì)航天器SC1、SC2和SC3相對(duì)速度變化曲線，圖中綠色曲線表示SC1和SC2之間的相對(duì)速度，變化量最大值?v1為1.82 m/s；紅色曲線表示SC1和SC3之間的相對(duì)速度，變化量最大值?v2為1.05 m/s；藍(lán)色曲線表示SC2和SC3之間的相對(duì)速度，變化量最大值?v3為0.92 m/s。

圖7 日心軌道編隊(duì)航天器之間的相對(duì)速度Fig.7 Heliocentric orbit formation relative velocity between spacecraft

上述結(jié)果表明，在考慮了太陽引力作用下，日心軌道空間引力波探測(cè)系統(tǒng)在一年內(nèi)的變化能滿足引力波探測(cè)的需求。日心編隊(duì)構(gòu)型軌道頻率變化緩慢，可以在較長(zhǎng)時(shí)間保持穩(wěn)定的構(gòu)型。

以優(yōu)化得到的結(jié)果作為初始條件，用Monte-Carlo法分析入軌誤差對(duì)日心軌道編隊(duì)構(gòu)型的影響。分別向航天器的徑向、切向和法向位置增加標(biāo)準(zhǔn)差為200 km的位置誤差，對(duì)構(gòu)型的影響如表7所示。分別向航天器的徑向、切向和法向速度增加標(biāo)準(zhǔn)差為1 cm/s的速度誤差，對(duì)構(gòu)型的影響如表8所示。位置誤差和速度誤差均符合正態(tài)分布?？梢钥闯?，徑向位置誤差和切向速度誤差導(dǎo)致構(gòu)型標(biāo)準(zhǔn)差的變化最大。若位置誤差和速度誤差同時(shí)存在，當(dāng)速度誤差為3 cm/s時(shí)，位置誤差最大不超過90 km，否則構(gòu)型將會(huì)發(fā)散。

表7 位置誤差對(duì)日心軌道編隊(duì)構(gòu)型的影響Table 7 Effect of position error on configuration of heliocentric orbital formation

表8 速度誤差對(duì)日心軌道編隊(duì)構(gòu)型的影響Table 8 Effect of velocity error on configuration of heliocentric orbital formation

地心軌道編隊(duì)和日心軌道編隊(duì)的優(yōu)化結(jié)果對(duì)比如表9所示，地心軌道編隊(duì)的臂長(zhǎng)變化率最大值為0.18%，呼吸角變化量最大值為0.14°，相對(duì)速度的變化量最大值為6.70 m/s；日心軌道編隊(duì)的臂長(zhǎng)變化率最大值為0.48%，呼吸角變化量最大值為0.38°，相對(duì)速度的變化量最大值為1.82 m/s。地心軌道編隊(duì)和日心軌道編隊(duì)的優(yōu)化結(jié)果均滿足各自的指標(biāo)要求，入軌誤差中對(duì)構(gòu)型穩(wěn)定性影響較大的均為徑向位置誤差和切向速度誤差。

表9 仿真算例數(shù)據(jù)對(duì)比Table 9 Comparision of simuilation results

4 結(jié) 論

本文以空間引力波探測(cè)器編隊(duì)構(gòu)型設(shè)計(jì)為研究背景，重點(diǎn)研究了探測(cè)器動(dòng)力學(xué)建模和軌道優(yōu)化算法。仿真結(jié)果表明，所提優(yōu)化算法可以使引力波探測(cè)器在持續(xù)一年的時(shí)間里，保持構(gòu)型的穩(wěn)定性，其中，無中心天體的日心軌道探測(cè)器編隊(duì)的臂長(zhǎng)變化率保持在0.48 %以內(nèi)，呼吸角最大變化量保持在0.38%以內(nèi)；地心軌道探測(cè)器編隊(duì)的相對(duì)速度變化量保持在6.70 m/s以內(nèi)。另外，徑向位置誤差和切向速度誤差對(duì)編隊(duì)構(gòu)型穩(wěn)定性的影響較大。上述結(jié)果說明，所提基于對(duì)偶四元數(shù)編隊(duì)動(dòng)力學(xué)模型，可以對(duì)探測(cè)器有無中心天體的情況進(jìn)行統(tǒng)一描述，在保證動(dòng)力學(xué)模型精度的同時(shí)提升模型的適用性。本文的研究成果能夠?yàn)槲磥黹_展空間引力波探測(cè)提供建議與參考。