焦博涵，黨朝輝

（1.西北工業(yè)大學(xué) 航天學(xué)院，西安 710072；2.西北工業(yè)大學(xué) 航天飛行動(dòng)力學(xué)技術(shù)重點(diǎn)實(shí)驗(yàn)室，西安 710072）

引 言

2015年9月14日激光干涉引力波天文臺(tái)（Laser Interferometer Gravitational-Wave Observatory，LIGO）合作小組第一次探測(cè)到了引力波信號(hào)，這一發(fā)現(xiàn)填補(bǔ)了廣義相對(duì)論實(shí)驗(yàn)驗(yàn)證的最后一塊缺失的拼圖，同時(shí)也揭開了人們利用引力波進(jìn)行宇宙研究的序幕。而隨著研究的深入，傳統(tǒng)地面引力波探測(cè)裝置由于地面噪聲、臂長(zhǎng)等限制的存在已越發(fā)不能滿足研究所需的探測(cè)需求[1]。為此一系列采用衛(wèi)星編隊(duì)進(jìn)行引力波探測(cè)的空間任務(wù)被陸續(xù)提出。歐洲航天局（European Space Agency，ESA）與美國(guó)國(guó)家航空航天局（National Aeronautics and Space Administration，NASA）于1993年最早提出了空間引力波探測(cè)計(jì)劃LISA[2]項(xiàng)目。隨后，日本、歐洲等也先后提出了類似的計(jì)劃如DECIGO[3]、BBO[4]、ALIA[5]等。中國(guó)從2008年開始探討空間引力波探測(cè)的可行性，先后開展了概念與方案設(shè)計(jì)、關(guān)鍵科學(xué)載荷研制等，并于2014年和2016年提出“天琴計(jì)劃”[6]和“太極計(jì)劃”[7]。

空間引力波探測(cè)任務(wù)是利用3個(gè)航天器構(gòu)成空間正三角形編隊(duì)，根據(jù)邁克爾遜激光干涉原理，通過測(cè)量相鄰兩個(gè)航天器間的臂長(zhǎng)變化來(lái)表征引力波信號(hào)[8-9]。因此如何設(shè)計(jì)一個(gè)可以長(zhǎng)期穩(wěn)定保持的正三角形編隊(duì)便成為了空間引力波探測(cè)任務(wù)的重要問題，對(duì)此相關(guān)學(xué)者進(jìn)行了大量的研究并給出了3類設(shè)計(jì)方法[9]：①共軌星座方式，即將3個(gè)航天器均勻地布置在同一條圍繞中心天體的圓形軌道上；②相對(duì)繞飛方式，即將3個(gè)航天器均勻地布置在圓參考軌道附近的相對(duì)圓形繞飛軌道上，需要額外說明的是：這類方法在具體進(jìn)行設(shè)計(jì)時(shí)又可分為基于絕對(duì)運(yùn)動(dòng)的設(shè)計(jì)方法（通過確定一個(gè)航天器的絕對(duì)運(yùn)動(dòng)軌道并根據(jù)正三角形編隊(duì)的幾何關(guān)系，將其繞垂直軌道平面的方向旋轉(zhuǎn)，得到另外兩個(gè)航天器的軌道[10-12]）和基于相對(duì)運(yùn)動(dòng)的設(shè)計(jì)方法〔利用CW（Clohessy-Wiltshire）方程構(gòu)造參考軌道附近的空間圓形相對(duì)軌道，并將3個(gè)航天器均勻布置在該相對(duì)軌道上從而構(gòu)成三角形編隊(duì)[12-13]〕；③三角平動(dòng)點(diǎn)方式，即將3個(gè)航天器分別布置在日地圓形限制性三體模型中的L3、L4、L5共3個(gè)平動(dòng)點(diǎn)上。針對(duì)上述3類設(shè)計(jì)方法，由于基于三角平動(dòng)點(diǎn)方式的設(shè)計(jì)方法的工程實(shí)現(xiàn)難度較大，因此對(duì)于絕大多數(shù)的空間引力波探測(cè)任務(wù)通常采用前兩類設(shè)計(jì)方法。但不論是共軌星座方式還是相對(duì)繞飛方式，其均為理想化二體引力場(chǎng)下（無(wú)攝動(dòng)或線性引力場(chǎng)）的設(shè)計(jì)結(jié)果。因此會(huì)導(dǎo)致編隊(duì)在高精度引力場(chǎng)運(yùn)行過程中出現(xiàn)較為嚴(yán)重發(fā)散情況（具體表現(xiàn)在正三角形編隊(duì)臂長(zhǎng)、臂長(zhǎng)變化率、呼吸角等指標(biāo)的發(fā)散上）[14-16]，從而無(wú)法滿足引力波探測(cè)任務(wù)所需的編隊(duì)穩(wěn)定性。

為解決上述編隊(duì)構(gòu)形長(zhǎng)時(shí)間運(yùn)行的穩(wěn)定性問題，有不同研究提出了若干解決方法。針對(duì)這些方法，從其所依據(jù)的模型來(lái)看，可將其分為基于Relative Orbital Elements（ROE）方程[17]的編隊(duì)設(shè)計(jì)、基于絕對(duì)軌道動(dòng)力學(xué)方程的編隊(duì)設(shè)計(jì)以及基于二階CW方程[10]的編隊(duì)設(shè)計(jì)共3類。①基于ROE方程的編隊(duì)設(shè)計(jì)：ROE方程是指基于相對(duì)軌道根數(shù)差的相對(duì)運(yùn)動(dòng)方程。清華大學(xué)蔣方華副教授團(tuán)隊(duì)[18]基于ROE方程，通過3個(gè)約束條件，將決策變量由18個(gè)減少到14個(gè)，并提出了一種變量區(qū)域自適應(yīng)調(diào)整算法，實(shí)現(xiàn)了高效優(yōu)化，而Joffre等[19]則進(jìn)一步考慮了太陽(yáng)系內(nèi)主要天體的攝動(dòng)影響。中國(guó)科學(xué)院空間應(yīng)用工程與技術(shù)中心張皓研究員[20]在ROE方程的基礎(chǔ)上提出了基于非奇異軌道根數(shù)的相對(duì)運(yùn)動(dòng)方程，并通過分析優(yōu)化指標(biāo)間的相關(guān)性將多目標(biāo)優(yōu)化問題轉(zhuǎn)化為了多個(gè)單目標(biāo)優(yōu)化問題。②基于絕對(duì)軌道動(dòng)力學(xué)方程：李卓[21]在二體模型的基礎(chǔ)上，考慮多種攝動(dòng)因素，采用辛積分方法實(shí)現(xiàn)了構(gòu)形優(yōu)化。易照華[22]著重關(guān)注LISA任務(wù)編隊(duì)中心與地球同軌的特點(diǎn)，建立了共軌限制性三體問題，并推導(dǎo)了共軌限制性三體模型下的編隊(duì)臂長(zhǎng)的表達(dá)式，但卻并未研究如何利用該表達(dá)式進(jìn)行編隊(duì)優(yōu)化。③基于二階CW方程的編隊(duì)設(shè)計(jì)：Nayak[10]基于二階CW方程推導(dǎo)了編隊(duì)臂長(zhǎng)與編隊(duì)平面傾角的關(guān)系，發(fā)現(xiàn)可通過改變編隊(duì)平面傾角來(lái)降低編隊(duì)臂長(zhǎng)的波動(dòng)峰值。Pucacco[23]和Dhurandhar[24]在二階CW方程的基礎(chǔ)上考慮不同的攝動(dòng)因素進(jìn)行了進(jìn)一步的優(yōu)化。雖然Nayak等基于二階CW方程對(duì)編隊(duì)構(gòu)形進(jìn)行了優(yōu)化設(shè)計(jì)，但其仍存在以下兩點(diǎn)不足：①Nayak所得到的是一種具有二階精度的CW方程周期解而并非二階CW方程的解析解，因此并不能直接地解釋構(gòu)形發(fā)散地深層機(jī)理；②Nayak的方法計(jì)算較為復(fù)雜且不直觀，無(wú)法解釋更優(yōu)穩(wěn)定解的物理機(jī)制，且構(gòu)型穩(wěn)定性的優(yōu)化結(jié)果仍然有限。

針對(duì)Nayak方法的局限性，本文提出了一種基于二階CW方程的編隊(duì)構(gòu)形設(shè)計(jì)方法：基于攝動(dòng)法推導(dǎo)了二階CW方程的近似解析解，得益于這種解析解，成功解釋了編隊(duì)構(gòu)形發(fā)散的主要原因，并依據(jù)分析結(jié)果給出了一種以航天器相位角為優(yōu)化變量的構(gòu)形優(yōu)化方法。該方法相較于Nayak的方法具有更好的優(yōu)化結(jié)果，且在優(yōu)化臂長(zhǎng)指標(biāo)的基礎(chǔ)上能夠同時(shí)實(shí)現(xiàn)對(duì)于呼吸角指標(biāo)的優(yōu)化。

1 空間引力波探測(cè)編隊(duì)標(biāo)稱構(gòu)形設(shè)計(jì)

本小節(jié)首先回顧太極計(jì)劃對(duì)于編隊(duì)構(gòu)形設(shè)計(jì)的主要要求。其次，給出描述相對(duì)運(yùn)動(dòng)的CW方程，并在此基礎(chǔ)上通過構(gòu)建空間圓形繞飛軌道實(shí)現(xiàn)標(biāo)稱構(gòu)形的設(shè)計(jì)。

1.1 太極引力波探測(cè)計(jì)劃

作為中國(guó)的空間引力波探測(cè)計(jì)劃，太極計(jì)劃于2016年由中國(guó)科學(xué)院正式提出，旨在探測(cè)頻段位于0.1～1 mHz范圍內(nèi)的引力波源。太極任務(wù)由3顆完全相同的航天器組成，它們以 20?的角度沿超前或落后地球的日心軌道運(yùn)行。如圖1所示，3顆航天器形成一個(gè)臂長(zhǎng)約300萬(wàn)km的等邊三角形。

圖1 太極引力波探測(cè)計(jì)劃Fig.1 Taiji gravitational wave detection program

與LISA計(jì)劃類似，太極計(jì)劃同樣利用每顆航天器上所攜帶的一對(duì)激光測(cè)距干涉儀和兩個(gè)重力參考傳感器，基于正三角形編隊(duì)構(gòu)形形成三組邁克爾遜干涉儀實(shí)現(xiàn)對(duì)于引力波信號(hào)的探測(cè)。在實(shí)際工程中，若要實(shí)現(xiàn)上述原理的引力波探測(cè)任務(wù)，則需要對(duì)星間所形成的三組邁克爾遜干涉儀提出極高的精度要求，即需要對(duì)正三角形編隊(duì)的整體穩(wěn)定性提出一定的要求。

文獻(xiàn)[21]指出了任務(wù)過程中對(duì)于編隊(duì)的穩(wěn)定性要求（示意圖見圖2），其具體指標(biāo)要求見表1。在本文中，從簡(jiǎn)化問題的角度出發(fā)，將著重關(guān)注臂長(zhǎng)變化這一指標(biāo)。

表1 “太極”任務(wù)編隊(duì)穩(wěn)定性指標(biāo)Table 1 Formation stability index of Taiji mission

圖2 編隊(duì)穩(wěn)定性指標(biāo)示意圖Fig.2 Diagram of formation stability index

1.2 參考坐標(biāo)系

本文采用如下參考坐標(biāo)系用于描述空間引力波探測(cè)編隊(duì)的運(yùn)動(dòng)：

1）日心慣性坐標(biāo)系

本文日心慣性坐標(biāo)系的建立參考國(guó)際天體參考系[25]（International Celestial Reference System），其坐標(biāo)原點(diǎn)位于太陽(yáng)的質(zhì)心，X軸位于黃道面內(nèi)并指向J2000春分點(diǎn)，Z軸與黃道面垂直，與地球公轉(zhuǎn)速度方向一致，Y軸位于黃道面內(nèi)，與X、Z軸構(gòu)成右手坐標(biāo)系。

2）Hill坐標(biāo)系

Hill坐標(biāo)系常用于描述相對(duì)運(yùn)動(dòng)，其坐標(biāo)原點(diǎn)位于正三角形編隊(duì)的虛擬中心，x軸位于虛擬中心的軌道平面內(nèi)并由日心指向虛擬中心，z軸垂直于虛擬中心的軌道平面并指向角動(dòng)量方向，y軸分別與另外兩軸垂直從而構(gòu)成右手坐標(biāo)系。該坐標(biāo)系也常稱為當(dāng)?shù)卮怪碑?dāng)?shù)厮剑↙VLH）坐標(biāo)系。

上述兩坐標(biāo)系如圖3所示。

圖3 坐標(biāo)系示意圖Fig.3 Diagram of coordinate system

1.3 編隊(duì)運(yùn)動(dòng)方程

在不考慮攝動(dòng)力與控制力的情況下，若假設(shè)：①參考航天器的運(yùn)動(dòng)軌道為圓形軌道；②從航天器與參考航天器之間的距離遠(yuǎn)小于其各自的軌道半徑，則二體引力場(chǎng)中的線性相對(duì)運(yùn)動(dòng)動(dòng)力學(xué)方程可在Hill坐標(biāo)系中表示為：

在不考慮攝動(dòng)力與控制力的情況下，由1.1節(jié)可知：①太極計(jì)劃的編隊(duì)虛擬中心位于偏心率接近為0的地球公轉(zhuǎn)軌道；②編隊(duì)虛擬中心距離太陽(yáng)的距離為1 AU遠(yuǎn)大于3×106km的編隊(duì)尺度。因此針對(duì)太極計(jì)劃編隊(duì)的上述兩個(gè)特點(diǎn)，可采用如下建立在Hill坐標(biāo)系中的CW方程來(lái)描述整個(gè)編隊(duì)的運(yùn)動(dòng)[26]

其中：x0、y0、z0、、、分別表示t=0時(shí)刻相對(duì)位置和相對(duì)速度的值。

1.4 基于CW方程的標(biāo)稱構(gòu)形設(shè)計(jì)

在長(zhǎng)時(shí)間的任務(wù)周期中，為了保持編隊(duì)整體的有界性，必須消除式（2）中的長(zhǎng)期項(xiàng)，由此可得到CW方程的周期性運(yùn)動(dòng)條件為

將該條件代入到式（2）中便可得到CW方程的周期解為

正三角形編隊(duì)構(gòu)形設(shè)計(jì)的關(guān)鍵在于空間圓形相對(duì)軌道的構(gòu)建。而所謂空間圓形是指從航天器相對(duì)于參考航天器間的距離為一定值，因此有

進(jìn)而將式（4）代入便可得空間圓形編隊(duì)的初始條件為

由此若要獲得正三角形編隊(duì)，則只需要在上式所確定的空間圓形相對(duì)軌道上布置三個(gè)相位角相差120?的航天器即可，即

而式（6）中z方向與x方向在相對(duì)位置和相對(duì)速度上的正負(fù)比例關(guān)系則代表了順時(shí)針旋轉(zhuǎn)和逆時(shí)針旋轉(zhuǎn)兩類構(gòu)形模式[18-19]。

綜上，式（6）和（7）即為基于CW方程所設(shè)計(jì)的正三角形標(biāo)稱構(gòu)形。

2 二階CW方程及其近似解析解

在上一節(jié)中介紹了描述相對(duì)運(yùn)動(dòng)的CW方程，其是在非線性相對(duì)運(yùn)動(dòng)動(dòng)力學(xué)的基礎(chǔ)上，將非線性的引力梯度項(xiàng)進(jìn)行Taylor級(jí)數(shù)展開并保留一階項(xiàng)所得到的線性方程。若對(duì)非線性的引力梯度項(xiàng)的Taylor級(jí)數(shù)保留至二階項(xiàng)，則可獲得式（8）所示的二階CW方程

由于二階CW方程相較于CW方程保留了一部分二階非線性項(xiàng)，因此其具有更高的精度，能夠體現(xiàn)出更多的非線性規(guī)律。為了充分利用二階CW方程中所蘊(yùn)含的非線性信息，望求取式（8）的解析解。因式（8）為一組非線性微分方程組，無(wú)法直接寫出其通解的形式，因此考慮使用攝動(dòng)法求解式（8）的近似解析解。

攝動(dòng)法（或稱小參數(shù)法）是一種將非線性因素作為對(duì)線性系統(tǒng)的攝動(dòng)，從而在線性解基礎(chǔ)上尋求非線性系統(tǒng)近似解的方法[27]。據(jù)攝動(dòng)法核心思想式（8）的解可以寫為

其中：ε即為攝動(dòng)法中的小參數(shù)，通過選擇保留ε的不同階次項(xiàng)即可獲得不同精度的近似解析解。對(duì)于本文所要研究的問題，為避免使問題過度復(fù)雜，本文僅選擇保留到ε的一次項(xiàng)。

對(duì)于式（8）可選取ε=n2/a作為小參數(shù)，此時(shí)根據(jù)式（9），在保留ε的一次項(xiàng)基礎(chǔ)上將近似解析解的形式代入到式（8）中，則可將式（8）的求解問題轉(zhuǎn)化為如下兩個(gè)微分方程組的求解問題

對(duì)于式（10）的方程組而言，它的解即為CW方程的解析解。進(jìn)一步，將式（10）的解代入到式（11）微分方程組的右端，并根據(jù)線性常微分方程理論，可求得式（11）的解為

其中：α0～α8、β0～β7、γ0～γ6分別為與x0、y0、z0、、、有關(guān)的代數(shù)表達(dá)式，其具體形式為

進(jìn)而，根據(jù)式（2）、式（12）可得到二階CW方程的近似解析解為

為了直觀地說明二階CW方程在精度上相較于CW方程的優(yōu)越性，分別采用CW方程的解析解〔式（2）〕、二階CW方程的近似解析解〔式（36）〕以及二體非線性方程，以標(biāo)稱構(gòu)形設(shè)計(jì)結(jié)果為初始條件進(jìn)行演化計(jì)算，并分別繪制CW方程解析解、二階CW方程近似解析解的演化結(jié)果與二體方程的演化結(jié)果的差的變化情況，其結(jié)果如圖4所示：

圖4 CW方程與二階CW方程的精度比較Fig.4 Accuracy comparison between CW equation and second-order CW equation

可以看出，二階CW方程相較于CW方程具有更好的精度，也因此可以充分利用二階CW方程的這一優(yōu)勢(shì)進(jìn)行編隊(duì)設(shè)計(jì)。

3 編隊(duì)構(gòu)形的分析與優(yōu)化

在第2節(jié)中推導(dǎo)并給出了二階CW方程的近似解析解表達(dá)式，本節(jié)將基于該近似解析解針對(duì)空間引力波探測(cè)編隊(duì)構(gòu)形的設(shè)計(jì)與優(yōu)化問題進(jìn)行分析。

3.1 完美圓形繞飛軌道不存在性證明

值得說明的是，對(duì)于以CW方程為基礎(chǔ)的正三角形編隊(duì)構(gòu)形優(yōu)化設(shè)計(jì)問題，另一種描述為：在不同精度的引力場(chǎng)中是否存在完美的圓形相對(duì)繞飛軌道[9]。對(duì)于這個(gè)問題，盡管是最為簡(jiǎn)單的二體引力場(chǎng)也是難以回答的，而在獲得能夠表征一定二體引力場(chǎng)非線性特征的二階CW方程的基礎(chǔ)上，本節(jié)嘗試在二階CW方程下回答這一問題。

同基于CW方程的標(biāo)稱構(gòu)形思路一致，為了構(gòu)成穩(wěn)定的空間圓形構(gòu)形則需要消除式（36）中的長(zhǎng)期項(xiàng)和漂移項(xiàng)，這要求如下關(guān)系式成立

至此在滿足式（37）的約束下得到了二階CW方程的周期解，進(jìn)一步需要結(jié)合式（5）給出二階CW方程下形成空間圓形編隊(duì)所需要滿足的條件。先將式（37）代入式（36），隨后再將式（36）代入到式（5）中，忽略關(guān)于ε的二次項(xiàng)，且為保證最后的結(jié)果為一常值需令sin2nt和cos2nt外其余項(xiàng)的系數(shù)均為0，由此有下列各式成立

式（37）、（38）的意義在于當(dāng)兩式中的所有表達(dá)式同時(shí)成立時(shí)才可以形成空間圓形編隊(duì)。此外，由式（37）、（38）可以看出，二階CW方程的空間圓形編隊(duì)條件是在CW方程空間圓形編隊(duì)條件的基礎(chǔ)上又增添了若干約束，當(dāng)忽略關(guān)于ε的一次項(xiàng)時(shí)，兩式便退化為CW方程的空間圓形編隊(duì)條件，進(jìn)一步驗(yàn)證了二階CW方程的正確性。

下面利用式（38）中的表達(dá)式進(jìn)行線性組合可得到如下表達(dá)式

進(jìn)一步化簡(jiǎn)可以推出

若要式（40）成立則需要滿足z0==0，此時(shí)考慮到ε為一近似為0的小量，則若式（38）中的η0成立則需要滿足x0≈≈0，而這種情況顯然不會(huì)出現(xiàn)，這說明在考慮 η0、η1、η2、η3、η44個(gè)表達(dá)式均成立的假設(shè)下推出了矛盾的結(jié)果，因此可以說明 η0、η1、η2、η3、η4不可能同時(shí)成立，而這則說明在二階CW方程下不存在嚴(yán)格精確的空間圓形繞飛軌道。

3.2 基于二階CW方程的構(gòu)形發(fā)散性分析

由3.1節(jié)可知，在二階CW方程下不存在嚴(yán)格精確的空間圓形繞飛軌道，因此進(jìn)一步可以說明在二階CW方程的精度下無(wú)法通過完全解析的方法實(shí)現(xiàn)對(duì)于正三角形編隊(duì)的設(shè)計(jì)。為此，較為自然且直接的思路是如何利用已有的有限解析結(jié)果來(lái)輔助編隊(duì)的優(yōu)化設(shè)計(jì)。以此為目的，本小節(jié)旨在通過二階CW方程分析基于CW方程的標(biāo)稱構(gòu)形設(shè)計(jì)結(jié)果在二體非線性引力場(chǎng)下的發(fā)散機(jī)理，從而為后續(xù)的編隊(duì)優(yōu)化設(shè)計(jì)提供依據(jù)。

根據(jù)式（6），首先考慮周期性條件以及無(wú)漂移條件，即將y˙0=?2nx0和y0=/n代入式（36），可得

通過比較式（41）和式（36），并觀察包含t的項(xiàng)的變化，可以發(fā)現(xiàn)：標(biāo)稱構(gòu)形的周期性條件消除了二階CW方程近似解析解中的x方向的發(fā)散項(xiàng)，但并未消除y方向的發(fā)散項(xiàng)；標(biāo)稱構(gòu)形的無(wú)漂移條件消除了二階CW方程近似解析解中x和y方向的主要漂移項(xiàng)。

進(jìn)一步，分析式（6）中x方向與z方向在相對(duì)位置、相對(duì)速度上的關(guān)系，不喪失一般性可將z0=和代入到式（41）中。結(jié)果顯而易見，同樣從包含t的項(xiàng)的變化來(lái)看，上述兩條件并未對(duì)式（41）的形式造成影響。但通過觀察y方向發(fā)散項(xiàng)系數(shù) β1的具體表達(dá)式〔式（42）〕可以發(fā)現(xiàn)：除x0==0的情況外，不論何種相對(duì)初始狀態(tài)，其構(gòu)形在二階CW方程下均會(huì)發(fā)散，然而根據(jù)式（6）中x0和x的表達(dá)式，不會(huì)出現(xiàn)x0==0的情況。綜合上述分析可以發(fā)現(xiàn)基于CW方程的標(biāo)稱構(gòu)形設(shè)計(jì)結(jié)果在二階CW方程下必然會(huì)出現(xiàn)發(fā)散的情況，即二體引力場(chǎng)中的非線性部分會(huì)破壞標(biāo)稱構(gòu)形的周期條件，因此在進(jìn)行編隊(duì)優(yōu)化設(shè)計(jì)的過程中應(yīng)著重考慮周期條件的修正

3.3 基于二階CW方程的編隊(duì)優(yōu)化設(shè)計(jì)

在分析基于CW方程的標(biāo)稱構(gòu)形設(shè)計(jì)結(jié)果在二階CW方程下的發(fā)散機(jī)理的基礎(chǔ)上，本節(jié)主要研究如何通過對(duì)已有的標(biāo)稱設(shè)計(jì)結(jié)果進(jìn)行修正，從而實(shí)現(xiàn)編隊(duì)的優(yōu)化設(shè)計(jì)。

1）周期匹配條件與能量匹配條件

根據(jù)3.2節(jié)的分析結(jié)果，二體引力場(chǎng)下基于CW方程的標(biāo)稱構(gòu)形設(shè)計(jì)結(jié)果的發(fā)散原因是因?yàn)槎w引力場(chǎng)中的非線性部分破壞了標(biāo)稱構(gòu)形的周期條件。因此在研究編隊(duì)的優(yōu)化設(shè)計(jì)問題之前，需要針對(duì)二體引力場(chǎng)下的相對(duì)運(yùn)動(dòng)周期性條件進(jìn)行討論。

對(duì)于二體引力場(chǎng)下的相對(duì)運(yùn)動(dòng)，其周期性相對(duì)運(yùn)動(dòng)條件存在兩種形式上的表達(dá)[28]：一種是由軌道根數(shù)所描述的周期匹配條件，要求編隊(duì)中的從航天器和參考航天器具有相同的半長(zhǎng)軸；另一種是由笛卡爾坐標(biāo)所描述的能量匹配條件，要求從航天器的相對(duì)運(yùn)動(dòng)初始狀態(tài)滿足式（43），其中：R表示參考航天器的軌道半徑。

為了說明上述兩個(gè)周期性相對(duì)運(yùn)動(dòng)條件的有效性，圖5和圖6分別給出了基于周期匹配條件和能量匹配條件修正后太極計(jì)劃編隊(duì)臂長(zhǎng)的演化情況（注：SC1、SC2、SC3的相位角依次為0?、120?、240?）。

圖6 二體引力場(chǎng)下基于能量匹配條件修正后的臂長(zhǎng)演化情況Fig.6 Evolution of arm length based on modified energy matching condition in two-body gravitational field

從整體來(lái)看，經(jīng)周期匹配條件和能量匹配條件修正后，與圖7相比編隊(duì)臂長(zhǎng)在10年任務(wù)周期內(nèi)不再呈現(xiàn)發(fā)散趨勢(shì)，其演化表現(xiàn)出明顯的周期性特征。但具體來(lái)看，經(jīng)周期匹配條件修正后，航天器1與航天器2以及航天器1與航天器3之間的臂長(zhǎng)變化最大值已經(jīng)不滿足太極計(jì)劃的臂長(zhǎng)要求；而經(jīng)能量匹配修正后，3顆航天器兩兩之間的臂長(zhǎng)變化最大值均能夠滿足任務(wù)要求，且臂長(zhǎng)波動(dòng)范圍能夠保持在較細(xì)小的范圍之內(nèi)。

圖7 二體引力場(chǎng)下標(biāo)稱構(gòu)形臂長(zhǎng)演化情況Fig.7 Evolution of arm length of nominal configuration in two-body gravitational field

盡管周期匹配條件和能量匹配條件僅在表達(dá)方式上有所不同，但其在應(yīng)用于編隊(duì)設(shè)計(jì)的過程中卻表現(xiàn)出了不同的效果。通過簡(jiǎn)單分析可以發(fā)現(xiàn)，造成這種現(xiàn)象的原因與標(biāo)稱構(gòu)形的設(shè)計(jì)方法有著緊密的聯(lián)系：由于本文所采用的是基于CW方程的相對(duì)運(yùn)動(dòng)設(shè)計(jì)方法，所得到的設(shè)計(jì)結(jié)果均基于笛卡爾坐標(biāo)所表示。對(duì)于周期匹配修正而言，雖然其修正了由于非線性所破環(huán)的周期條件，但同時(shí)在由笛卡爾坐標(biāo)向軌道根數(shù)進(jìn)行相互轉(zhuǎn)換的過程中同樣破環(huán)了式（6）中的其它條件；但對(duì)于能量匹配修正，由式（43）可以看出，其僅對(duì)式（6）中的周期條件進(jìn)行了修正而并未對(duì)其他條件造成影響。因此，通過上述分析可以得到如下結(jié)論：在基于CW方程所設(shè)計(jì)的標(biāo)稱構(gòu)形的基礎(chǔ)上進(jìn)行優(yōu)化設(shè)計(jì)時(shí)，能量匹配條件是一項(xiàng)行之有效的約束條件。

2）編隊(duì)優(yōu)化問題的建立

結(jié)合上述針對(duì)相對(duì)運(yùn)動(dòng)周期性條件的分析結(jié)果，可以進(jìn)一步構(gòu)建基于二階CW方程的編隊(duì)優(yōu)化模型。具體地，在優(yōu)化指標(biāo)方面，可直接基于式（36）的近似解析解構(gòu)造如下的臂長(zhǎng)變化量的解析表達(dá)式

其中：l0表示編隊(duì)的標(biāo)稱臂長(zhǎng)；下標(biāo)i和j表示航天器的編號(hào)，滿足i=1,2,3，j=1,2,3且i≠j。與現(xiàn)有文獻(xiàn)中常見的直接利用數(shù)值方法求解非線性方程進(jìn)而構(gòu)造臂長(zhǎng)的方法相比，解析的臂長(zhǎng)表達(dá)式雖然在精度上略有不足，但其在計(jì)算效率上卻具有顯著優(yōu)勢(shì)。在約束方面，根據(jù)上文對(duì)于相對(duì)運(yùn)動(dòng)周期性條件的分析結(jié)果，可直接采用式（43）作為約束條件。

在優(yōu)化變量方面，通過分析圖5和圖6的結(jié)果可以發(fā)現(xiàn)，航天器1/航天器2與航天器1/航天器3的臂長(zhǎng)變化始終處于同一量級(jí)，且大于航天器2/航天器3的臂長(zhǎng)。從幾何的角度來(lái)看，說明演化過程中的編隊(duì)構(gòu)形更趨向于等腰三角形的形狀而非正三角形。受這一現(xiàn)象啟發(fā)，在優(yōu)化過程中可針對(duì)式（6）中的航天器相位角α進(jìn)行修正（如圖8所示），以期望使構(gòu)形更趨向于正三角形。

圖8 優(yōu)化變量的幾何意義Fig.8 Geometric meaning of optimization variables

綜上，最終得到基于二階CW方程的編隊(duì)優(yōu)化模型為

其中：[t0,tf]表示任務(wù)周期；δα表示優(yōu)化變量的邊界值；非線性函數(shù)f代表能量匹配條件，即式（43）。通過求解式（45）的優(yōu)化問題，可進(jìn)一步計(jì)算出修正后正三角形編隊(duì)構(gòu)形條件

4 仿真校驗(yàn)與結(jié)果分析

本節(jié)以太極計(jì)劃為任務(wù)背景，對(duì)3.3節(jié)中所構(gòu)建的優(yōu)化問題〔式（45）〕進(jìn)行驗(yàn)證。根據(jù)式（45）可知，其本質(zhì)上是一個(gè)雙層優(yōu)化問題：內(nèi)層是以時(shí)間為優(yōu)化變量，以臂長(zhǎng)變化量最大為性能指標(biāo)的優(yōu)化問題；外層是以航天器相位角為優(yōu)化變量，以任務(wù)周期內(nèi)最大臂長(zhǎng)變化量最小為性能指標(biāo)的優(yōu)化問題。對(duì)于內(nèi)層優(yōu)化問題，根據(jù)圖5～圖7可以看出：對(duì)于不同的初始條件，臂長(zhǎng)的變化均存在一定頻率的波動(dòng)現(xiàn)象，對(duì)于優(yōu)化問題而言其意味著存在若干的局部極值點(diǎn)，因此對(duì)于內(nèi)層優(yōu)化問題應(yīng)采用全局優(yōu)化算法進(jìn)行求解，本文基于MATLAB中的GlobalSearch函數(shù)進(jìn)行求解。對(duì)于外層優(yōu)化問題，由于考慮?α1=0，其實(shí)質(zhì)上退化為了雙變量?jī)?yōu)化問題，同時(shí)考慮到指標(biāo)函數(shù)并非連續(xù)函數(shù)，因此采用模式搜索算法進(jìn)行求解。

圖9給出了經(jīng)優(yōu)化修正后二體引力場(chǎng)下編隊(duì)臂長(zhǎng)的演化情況。通過對(duì)比圖6可以發(fā)現(xiàn)：從總體來(lái)看，經(jīng)優(yōu)化修正后相比于基于能量匹配的修正方法，3個(gè)臂長(zhǎng)的平均最大變化量由0.42%降低至0.32%，而3個(gè)臂長(zhǎng)的最大變化量由0.62%降低至0.44%，且相較于Nayak[10]將臂長(zhǎng)最大變化量?jī)?yōu)化為1%的結(jié)果，可以看出本文所提出方法的優(yōu)勢(shì)。具體來(lái)看，經(jīng)優(yōu)化修正后航天器2和航天器3之間臂長(zhǎng)的變化幅值得到了減小，而航天器1和航天器2以及航天器1和航天器3之間的臂長(zhǎng)變化幅值的減小量并不明顯，而這一現(xiàn)象則正好印證了式（45）中選取航天器相位角作為優(yōu)化變量的有效性。

圖9 二體引力場(chǎng)下通過優(yōu)化修正后的臂長(zhǎng)演化情況Fig.9 Evolution of arm length after optimization correction in two-body gravitational field

除考量編隊(duì)的臂長(zhǎng)變化情況外，圖10和圖11分別給出了基于能量匹配修正和優(yōu)化修正后編隊(duì)呼吸角在10年任務(wù)周期內(nèi)的演化情況?？梢钥闯觯?jīng)優(yōu)化修正后相較于基于能量匹配的修正方法，編隊(duì)呼吸角的平均最大變化量由0.79%降低至0.6%，而呼吸角的最大變化量由0.92%降低至0.65%。雖然在式（45）的優(yōu)化問題中僅將編隊(duì)臂長(zhǎng)最大變化量作為優(yōu)化指標(biāo)，但通過修正航天器的相位角將編隊(duì)構(gòu)形由等腰三角形修正為正三角形，進(jìn)而使得呼吸角的變化幅值同樣得到了抑制，而這也體現(xiàn)出選擇航天器相位角作為優(yōu)化變量的一項(xiàng)優(yōu)勢(shì)。

圖10 二體引力場(chǎng)下基于能量匹配條件修正后的呼吸角演化情況Fig.10 Evolution of breathing angel length based on modified energy matching condition in two-body gravitational field

圖11 二體引力場(chǎng)下通過優(yōu)化修正后的呼吸角演化情況Fig.11 Evolution of breathing angel length after optimization correction in two-body gravitational field

5 結(jié) 論

本文研究了基于二階CW方程的空間引力波探測(cè)正三角形編隊(duì)的構(gòu)形設(shè)計(jì)問題，利用攝動(dòng)法給出了二階CW方程的近似解析解，并基于該解析解分析了基于CW方程設(shè)計(jì)的標(biāo)稱構(gòu)形的發(fā)散機(jī)理，構(gòu)造了編隊(duì)構(gòu)形優(yōu)化問題，最終得到如下結(jié)論：

1）在二階CW方程的模型精度下，不存在完美的圓形相對(duì)繞飛軌道；

2）基于CW方程所設(shè)計(jì)的標(biāo)稱構(gòu)形在二體非線性引力場(chǎng)下發(fā)散的主要原因是：二體引力場(chǎng)中的非線性部分破壞了標(biāo)稱構(gòu)形條件中的周期條件；

3）針對(duì)基于CW方程所設(shè)計(jì)的標(biāo)稱構(gòu)型，能量匹配修正相較于周期匹配修正具有更好的效果；

4）通過構(gòu)造基于二階CW方程的編隊(duì)優(yōu)化模型，在二體引力場(chǎng)下，將編隊(duì)臂長(zhǎng)10年內(nèi)的平均誤差降低至0.32%，最大誤差降低至0.44%，同時(shí)通過考慮三個(gè)臂長(zhǎng)變化的均衡性，以航天器的相位角為優(yōu)化變量，在僅以臂長(zhǎng)為優(yōu)化指標(biāo)的基礎(chǔ)上實(shí)現(xiàn)了對(duì)呼吸角的優(yōu)化。

本文的有關(guān)分析方法與仿真結(jié)果對(duì)太極引力波探測(cè)任務(wù)的編隊(duì)設(shè)計(jì)問題具有一定的借鑒意義。