李 卓，王有亮，鄭建華，李明濤*

(1.中國科學(xué)院 國家空間科學(xué)中心，北京 100190；2.中國科學(xué)院大學(xué)，北京 100039)

1 引 言

2016年，美國地基激光干涉引力波天文臺(tái)(LIGO)第一次完成地面引力波信號探測，檢驗(yàn)了廣義相對論的正確性，也預(yù)示人類可以通過引力波探測來認(rèn)識(shí)宇宙結(jié)構(gòu)和天體的演化。為了擴(kuò)展引力波探測頻段，對0.01 mHz～1 Hz左右中低頻引力波進(jìn)行探測，需要開展空間引力波探測活動(dòng)?？臻g引力波探測任務(wù)具有更寬廣的視野和大量的波源，涉及的關(guān)鍵技術(shù)包括無拖曳控制，高精度激光測距，微推進(jìn)技術(shù)，大尺度高穩(wěn)定編隊(duì)構(gòu)型設(shè)計(jì)等[1]。

20世紀(jì)90年代開始美國國家航空航天局(NASA)和歐洲航天局(ESA)共同提出了空間激光干涉引力波探測項(xiàng)目LISA(Laser Interferometer Space Antenna)，旨在觀測由超大質(zhì)量黑洞合并和由中子星或黑洞組成的雙星系統(tǒng)所產(chǎn)生的引力波[2]。LISA任務(wù)包括3個(gè)航天器，在地球軌道后方20°左右位置共同構(gòu)成臂長為2.5×106的等邊三角形。利用空間自由懸浮測試質(zhì)量塊作為傳感器，將引力波信號轉(zhuǎn)化為測試質(zhì)量塊間距變化的信號，也是干涉儀臂長的變化[3]。

中國科學(xué)家也提出了類似的引力波探測項(xiàng)目，其中包括太極計(jì)劃和天琴計(jì)劃。太極計(jì)劃由位于等邊三角形頂端的三顆衛(wèi)星組成，旨在探測中低頻波段的引力波。太極計(jì)劃的主要科學(xué)目標(biāo)是觀測雙黑洞并合和極大質(zhì)量比天體并合時(shí)產(chǎn)生的引力波輻射，以及其它的宇宙引力波輻射過程[1]。中山大學(xué)也提出了天琴計(jì)劃，天琴計(jì)劃在距地球105km的近圓軌道內(nèi)放置三顆航天器形成等邊三角形，來驗(yàn)證引力波在中低頻的存在和地球重力測量等其他觀測內(nèi)容[4-5]。

空間引力波探測通常會(huì)維持四年以上任務(wù)周期，由于初始入軌誤差等不確定擾動(dòng)的存在，使構(gòu)型發(fā)生變化，可能導(dǎo)致其不再滿足任務(wù)要求。因此需要對誤差傳播進(jìn)行分析。已知常用的方法有Monte-Carlo法，通過抽樣實(shí)驗(yàn)統(tǒng)計(jì)狀態(tài)量的均值和標(biāo)準(zhǔn)差。Monte-Carlo法原理簡單，對各種誤差普遍適用，然而為得到準(zhǔn)確值需要進(jìn)行多次抽樣，當(dāng)任務(wù)周期長時(shí)，計(jì)算量大，耗時(shí)長。

對于非線性誤差傳播分析問題，可以運(yùn)用協(xié)方差分析方法，混沌多項(xiàng)式展開法，狀態(tài)轉(zhuǎn)移張量法，微分代數(shù)方程法，高斯混合模型法，解Fokker-Plank 方程等方法解決[6]，工程實(shí)際中非線性誤差通常符合正態(tài)分布形式，需要根據(jù)實(shí)際情況選擇合適的方法。線性協(xié)方差方法多用于分析動(dòng)力學(xué)模型為線性函數(shù)的軌道預(yù)報(bào)、交會(huì)、轉(zhuǎn)移等問題，可快速獲得誤差傳播結(jié)果[7]。已有的研究結(jié)果可以解決線性模型下的空間引力波探測誤差傳播問題[8]，對于非線性模型，可用協(xié)方差分析描述函數(shù)法(CADET)解決。CADET法對非線性函數(shù)進(jìn)行統(tǒng)計(jì)線性化，然后通過積分計(jì)算狀態(tài)變量的均值和協(xié)方差矩陣。不同于Monte-Carlo法的耗時(shí)長和線性協(xié)方差法的精度低，CADET方法計(jì)算速度快且精度高，適用于空間引力波探測軌道誤差分析問題。

本文運(yùn)用CADET方法分析了空間引力波探測任務(wù)構(gòu)型在入軌誤差下的穩(wěn)定性，并研究了位置誤差和速度誤差的大小和方向?qū)?gòu)型穩(wěn)定性的影響。

2 空間引力波探測任務(wù)動(dòng)力學(xué)模型

選取J2000日心慣性坐標(biāo)系作為參考坐標(biāo)系，航天器受到太陽引力及其他天體的攝動(dòng)力。衛(wèi)星在空間運(yùn)動(dòng)中主要受太陽引力及水星、金星、地球、火星、木星、土星、天王星、海王星和冥王星的攝動(dòng)力作用。以LISA為例，對于相同初始條件，考慮以上所有力的作用時(shí)，其臂長變化約為2.9×105km，當(dāng)去掉金星、地球、木星的作用時(shí)，臂長變化分別變?yōu)?.4×105km，2.0×105km，2.7×105km；而去掉其他天體作用時(shí)臂長變化仍然保持在2.9×105km左右。因此為簡化模型，只考慮金星、地球、木星的攝動(dòng)力[8]?？臻g引力波探測任務(wù)動(dòng)力學(xué)模型表示為[9]:

(1)

其中，Rk表示三顆衛(wèi)星的位置矢量，k=1,2,3。Ri表示金星、地球和木星的位置矢量，i=1,2,3。μ0表示太陽引力常數(shù)，μi分別表示金星、地球和木星的引力常數(shù)，i=1,2,3。

定義Rij(t)=Rj(t)-Ri(t),(ij=12,23,31)為衛(wèi)星之間的相對位置矢量。衛(wèi)星編隊(duì)動(dòng)力學(xué)指標(biāo)分別為臂長L，呼吸角θ，臂長變化率V，星地距離D，可用公式(2)～(5)表達(dá)。

Lij(t)=‖Rij(t)‖ ,

(2)

(3)

(4)

D(t)=‖Re(t)-Rc(t)‖ ,

(5)

其中：Re是地球位置矢量，Rc是三衛(wèi)星中心位置矢量。

3 CADET在空間引力波探測任務(wù)軌道誤差分析中的應(yīng)用

若x(t)是系統(tǒng)狀態(tài)向量；F(t)是系統(tǒng)狀態(tài)矩陣，無過程擾動(dòng)的線性連續(xù)隨機(jī)系統(tǒng)微分方程是:

(6)

(7)

p(t)=E[δx(t)δxT(t)] ,

(8)

(9)

根據(jù)第二節(jié)，本文的動(dòng)力學(xué)模型是非線性系統(tǒng)，可以表示為：

(10)

其中，x(t)是方程的狀態(tài)變量，x1～x18為衛(wèi)星1,2,3的位置速度；x19～x21為臂長，x22～x24為呼吸角，x25～x27為臂長變化率；x28為星地距離。運(yùn)用CADET方法分析入軌誤差對空間引力波探測任務(wù)的影響，需對系統(tǒng)進(jìn)行統(tǒng)計(jì)線性化。

(11)

(12)

式(12)的計(jì)算比較復(fù)雜，工程中常常用如下方法簡化：

已知f(x)關(guān)于x一階可微，在m進(jìn)行f(x)的一階泰勒展開，得

(13)

忽略一階泰勒余項(xiàng)R1(x)，求式(13)期望值，f(x)的線性化形式是

(14)

動(dòng)力學(xué)方程的描述函數(shù)N為[12]

(15)

其中，I為28×28的單位矩陣，n為

(16)

式中，f1至f28包括下列方程：f1至f18為衛(wèi)星1，2，3的位置速度；f19至f21為臂長L12，L13，L23；f22至f24為呼吸角θ1，θ2，θ3；f25至f27為臂長變化率V12，V13，V23；f28為星地距離D。

n的具體形式是

(17)

其中，

(18)

(19)

i(i=1,2,3)的速度對位置的偏導(dǎo)。

(20)

(21)

其中:K包括臂長L、呼吸角θ、臂長變化率V和星地距離D，其中ij=12,13,23，k=1,2,3。

根據(jù)CADET方法知，具有初始誤差，空間引力波探測任務(wù)航天器在J2000日心慣性坐標(biāo)系中的狀態(tài)均值與協(xié)方差的微分方程如下[11]:

(22)

其中，初始條件m0包括三衛(wèi)星初始位置速度和臂長、呼吸角、臂長變化率和星地距離，p0為對角線矩陣，數(shù)值由三衛(wèi)星初始位置和速度誤差決定，N(t)由式(15)～(21)求出。根據(jù)以上條件，可通過積分求出任意時(shí)刻的均值和協(xié)方差矩陣。

4 CADET仿真驗(yàn)證

基于動(dòng)力學(xué)模型，以任務(wù)周期內(nèi)一直滿足指標(biāo)要求為標(biāo)準(zhǔn)選擇合適的位置速度初始條件，并加入正態(tài)分布誤差。分別用CADET法和Monte-Carlo法求任務(wù)周期10年內(nèi)構(gòu)型的臂長，呼吸角，臂長變化率和星地距離的均值和標(biāo)準(zhǔn)差。

圖1給出了CADET法和Monte-Carlo法求出的均值和標(biāo)準(zhǔn)差隨時(shí)間變化的結(jié)果。其中圓圈為Monte-Carlo法抽樣仿真1 000次結(jié)果，實(shí)線為CADET法結(jié)果。由圖1可知，兩種方法結(jié)果一致，以Monte-Carlo法的仿真結(jié)果作為標(biāo)準(zhǔn)，可以認(rèn)為CADET法是有效的。

圖1 CADET法和Monte-Carlo法的對照結(jié)果Fig.1 Comparison diagrams of CADET and Monte-Carlo

任務(wù)過程中將產(chǎn)生大量數(shù)據(jù)，所以選取一個(gè)特殊點(diǎn)(10年周期內(nèi)的最后時(shí)刻的數(shù)據(jù))進(jìn)行分析?？紤]初始誤差時(shí)CADET法和Monte-Carlo法均值、標(biāo)準(zhǔn)差誤差和相對誤差如表1所示。

由表1可知，臂長均值相對誤差不超過0.017 8%，標(biāo)準(zhǔn)差相對誤差不超過5.600 4%；呼吸角均值相對誤差不超過0.014 5%，標(biāo)準(zhǔn)差相對誤差不超過4.780 2%；臂長變化率均值相對誤差不超過1.633 6%，標(biāo)準(zhǔn)差相對誤差不超過5.339 7%；星地距離均值相對誤差不超過0.028 8%，標(biāo)準(zhǔn)差相對誤差不超過2.833 7%。說明CADET方法對于空間引力波探測誤差分析問題確實(shí)準(zhǔn)確有效。

表1 CADET法和Monte-Carlo法結(jié)果比較Tab.1 Results comparison of CADET and Monte-Carlo methods

5 基于CADET的空間引力波探測軌道誤差傳播分析

在空間引力波探測任務(wù)過程中，航天器面臨的入軌誤差是不確定的，不同的入軌誤差會(huì)對構(gòu)型產(chǎn)生不同程度的影響。當(dāng)使用Monte-Carlo法分析入軌誤差對構(gòu)型影響的不確定性時(shí)，需要多次抽樣保證結(jié)果的可靠性，因此解決問題的過程將花費(fèi)大量時(shí)間?，F(xiàn)在已知CADET方法對于空間引力波探測誤差分析問題準(zhǔn)確有效，可以運(yùn)用CADET方法進(jìn)行誤差傳播分析。

航天器的入軌誤差主要分為位置誤差和速度誤差，下面將對這兩類誤差對構(gòu)型指標(biāo)產(chǎn)生的影響進(jìn)行分析。以空間引力波探測任務(wù)為例，臂長要求為3×106km，呼吸角變化小于1°，臂長變化率小于10 m/s，星地距離小于6.5×107km時(shí)可認(rèn)為構(gòu)型滿足空間引力波探測任務(wù)要求。初始 構(gòu)型如圖2所示，1、2、3曲線分別代表衛(wèi)星1、2，衛(wèi)星1、3，衛(wèi)星2、3之間的指標(biāo)變化情況，其中臂長最大值為3.014×106km，呼吸角最大值為60.37°，臂長變化率最大值為3.325 m/s，標(biāo)準(zhǔn)差均為0。

圖2 初始構(gòu)型Fig.2 Initial constellations

5.1 誤差方向?qū)?gòu)型的影響

(1)不同方向的位置誤差

當(dāng)同時(shí)向3顆衛(wèi)星增加大小相等的位置誤差時(shí)，誤差的方向?qū)?gòu)型產(chǎn)生顯著影響。徑向誤差使構(gòu)型標(biāo)準(zhǔn)差產(chǎn)生的變化遠(yuǎn)大于其他兩個(gè)方向的誤差。徑向位置誤差增加了衛(wèi)星的勢能，所以構(gòu)型產(chǎn)生的擾動(dòng)更大。因此，航天器發(fā)射時(shí)需優(yōu)先考慮減小徑向位置誤差。

(2)不同方向的速度誤差

同時(shí)向衛(wèi)星1,2,3增加均值為1 cm/s，方向分別沿徑向，切向和法向的正態(tài)分布的速度誤差，構(gòu)型變化如表3所示。

當(dāng)同時(shí)向3顆衛(wèi)星增加大小相等的速度誤差

表2 位置誤差方向?qū)?gòu)型的影響Tab.2 Effect of position error direction on constellation

表3 速度誤差方向?qū)?gòu)型的影響Tab.3 Effect of velocity error direction on constellation

時(shí)，誤差的方向?qū)?gòu)型產(chǎn)生顯著影響。切向誤差使構(gòu)型標(biāo)準(zhǔn)差產(chǎn)生的變化遠(yuǎn)大于其他兩個(gè)方向的誤差。切向速度誤差增加了衛(wèi)星的動(dòng)能，所以構(gòu)型產(chǎn)生的擾動(dòng)更大。因此，航天器發(fā)射時(shí)需優(yōu)先考慮減小切向速度誤差。

5.2 誤差大小對構(gòu)型的影響

(1)不同大小的位置誤差

因?yàn)閺较蛭恢谜`差對構(gòu)型影響最大，因此在本節(jié)中僅以徑向位置誤差為例分析誤差大小對構(gòu)型的影響。

首先以任務(wù)周期4年為例，同時(shí)向衛(wèi)星1,2,3增加正態(tài)分布的徑向位置誤差，誤差均值分別為100 km，1 000 km，構(gòu)型變化表4所示。

同時(shí)向三顆衛(wèi)星增加徑向位置誤差時(shí)，位置誤差對構(gòu)型的影響不明顯。當(dāng)誤差均值為100 km和1 000 km時(shí)，只考慮均值時(shí)構(gòu)型指標(biāo)滿足要求，但觀察標(biāo)準(zhǔn)差發(fā)現(xiàn)，1 000 km時(shí)標(biāo)準(zhǔn)差增大，構(gòu)型發(fā)散不滿足指標(biāo)要求。

表4 誤差大小對構(gòu)型的影響(周期4年)Tab.4 Effect of error magnitude on constellation(T=4 years)

(2)不同大小的速度誤差

因?yàn)榍邢蛩俣日`差對構(gòu)型影響最大，因此在本節(jié)中僅以切向速度誤差為例分析誤差大小對構(gòu)型的影響。

同時(shí)向衛(wèi)星1,2,3增加正態(tài)分布的切向速度誤差，誤差均值分別為1 cm/s，10 cm/s，構(gòu)型變化如表4所示。

向三顆衛(wèi)星同時(shí)增加相同方向速度誤差時(shí)，速度誤差對構(gòu)型的影響不明顯。當(dāng)誤差均值為1 cm/s和10 cm/s時(shí)，只考慮均值時(shí)構(gòu)型指標(biāo)滿足要求，當(dāng)誤差均值為10 cm/s時(shí)雖然標(biāo)準(zhǔn)差增大，但由于任務(wù)周期短，構(gòu)型發(fā)散程度小，仍然滿足指標(biāo)要求。

為比較周期4年和10年時(shí)誤差大小對構(gòu)型的影響，對周期為10年的情況重復(fù)上述仿真，結(jié)果如表5所示。當(dāng)位置誤差為100 km和速度誤差為1 cm/s時(shí)構(gòu)型滿足指標(biāo)要求；當(dāng)位置誤差為1 000 km和速度誤差為10 cm/s時(shí)構(gòu)型不滿足指標(biāo)要求。當(dāng)比較任務(wù)周期為4年和10年兩種情況，發(fā)現(xiàn)標(biāo)準(zhǔn)差大小隨時(shí)間增加，任務(wù)周期延長時(shí)，大小相等的入軌誤差會(huì)使構(gòu)型產(chǎn)生更大擾動(dòng)。

表5 誤差大小對構(gòu)型的影響(周期10年)Tab.5 Effect of error magnitude on constellation(T=10 years)

5.3 三星間相對誤差方向?qū)?gòu)型的影響

當(dāng)航天器發(fā)射存在入軌誤差時(shí)，三衛(wèi)星誤差的相對方向可能會(huì)對構(gòu)型產(chǎn)生不同的影響。下文將探究相對方向相反的入軌誤差對構(gòu)型產(chǎn)生的影響。

(1)相對方向相反的位置誤差

以任務(wù)周期4年為例，分別向三衛(wèi)星徑向、切向、法向施加均值為100 km的正態(tài)分布的，相對方向相同和相反的位置誤差，構(gòu)型變化如表6所示。

表6 相對位置誤差方向?qū)?gòu)型的影響Tab.6 Effect of direction of relative position error on constellation

由于增加的位置誤差大小相等，方向相反，構(gòu)型標(biāo)準(zhǔn)差相同，反向誤差的擾動(dòng)大于同向誤差。當(dāng)航天器受到的位置誤差相對方向不同時(shí)，構(gòu)型產(chǎn)生的擾動(dòng)更大。

(2)相對方向相反的速度誤差

分別向三衛(wèi)星徑向、切向、法向施加均值為10 cm/s的正態(tài)分布的，相對方向相同和相反的速度誤差，構(gòu)型變化如表7所示。

表7 相對速度誤差方向?qū)?gòu)型的影響Tab.7 Effect of direction of relative velocity error on constellation

由于增加的速度誤差大小相等，方向相反，構(gòu)型標(biāo)準(zhǔn)差相同，反向誤差的擾動(dòng)大于同向誤差。當(dāng)航天器受到的速度誤差相對方向不同時(shí)，構(gòu)型產(chǎn)生的擾動(dòng)更大。

無論是位置誤差還是速度誤差，當(dāng)航天器受到的誤差相對方向相同時(shí)，誤差對構(gòu)型的影響明顯小于誤差相對方向相反時(shí)的情況。這也為任務(wù)入軌工作帶來一些啟示，如果無法減小入軌誤差的絕對值，可以通過使三衛(wèi)星受到的誤差保持相同方向，減小入軌誤差對構(gòu)型的影響，保證后續(xù)任務(wù)順利完成。

5.4 位置速度誤差對構(gòu)型的影響

入軌位置誤差和速度誤差通常同時(shí)存在，方向也可能不同，因此下文將分析兩種不同方向的誤差同時(shí)存在對構(gòu)型產(chǎn)生的影響。以周期4年為例，使切向速度誤差由0.5 cm/s逐漸增加至3 cm/s，并調(diào)整徑向位置誤差的大小，以臂長在3×106±3.5×104km內(nèi)為標(biāo)準(zhǔn)，記錄使構(gòu)型發(fā)散的位置誤差臨界點(diǎn)，如圖3所示。

圖3 位置誤差與速度誤差的關(guān)系Fig.3 Relationship between position error and velocity error

由圖3可知，速度誤差的影響大于位置誤差，對于同一構(gòu)型，當(dāng)速度誤差增大時(shí)，為滿足指標(biāo)要求需減小位置誤差。當(dāng)速度誤差是0.5 cm/s，位置誤差將超過595 km時(shí)構(gòu)型發(fā)散；當(dāng)速度誤差是3 cm/s，位置誤差超過160 km時(shí)構(gòu)型發(fā)散。

6 結(jié) 論

本文首先介紹了空間引力波探測任務(wù)的動(dòng)力學(xué)模型；運(yùn)用CADET法提出了空間引力波探測任務(wù)入軌誤差傳播方程；進(jìn)行CADET方法仿真和Monte-Carlo法仿真，并比較二者結(jié)果；以本文構(gòu)型為例，運(yùn)用CADET方法對入軌誤差的影響進(jìn)行分析。仿真結(jié)果表明：

(1)以Monte-Carlo法結(jié)果為標(biāo)準(zhǔn)，CADET法在10年任務(wù)周期內(nèi)與其相對誤差不超過6%， CADET法可以進(jìn)行有初始誤差情況下的任務(wù)過程的誤差傳播分析；

(2)在相同環(huán)境下完成Monte-Carlo法1 000次抽樣運(yùn)算平均用時(shí)約為15 min，而CADET法平均用時(shí)不超過1 min， CADET法可以實(shí)現(xiàn)任務(wù)過程中對誤差傳播的快速預(yù)測，用時(shí)遠(yuǎn)小于Monte-Carlo法。

(3)徑向位置誤差和切向速度誤差分別改變了衛(wèi)星的勢能和動(dòng)能，對空間引力波探測任務(wù)構(gòu)型影響較大。當(dāng)誤差的相對方向相同時(shí)，對構(gòu)型的影響小于相對方向不同時(shí)，可以通過保持誤差方向相同減小入軌誤差對構(gòu)型的影響。當(dāng)周期為4年且同時(shí)存在兩種誤差時(shí)，位置誤差不超過160 km，速度誤差不超過3 cm/s時(shí)，構(gòu)型仍滿足指標(biāo)要求。