基于蒙特卡洛法的衛(wèi)星分離姿態(tài)計(jì)算方法

康士朋，江 濤，耿盛韋，王金童

?

基于蒙特卡洛法的衛(wèi)星分離姿態(tài)計(jì)算方法

康士朋1，江 濤1，耿盛韋2，王金童1

（1. 上海宇航系統(tǒng)工程研究所，上海，201109；2. 南京航空航天大學(xué)航空宇航學(xué)院，南京，201106）

為了計(jì)算以彈簧作為分離動(dòng)力源的衛(wèi)星在星箭分離時(shí)刻的分離姿態(tài)，建立了六自由度動(dòng)力學(xué)方程，并采用龍格庫(kù)塔法求解了該方程，算例結(jié)果表明：通過(guò)該方法計(jì)算得到的分離角速度與ADAMS仿真計(jì)算結(jié)果相比，最大偏差小于0.77%。基于以上研究，針對(duì)衛(wèi)星、彈簧裝置多個(gè)設(shè)計(jì)偏差對(duì)衛(wèi)星分離姿態(tài)的影響，提出基于蒙特卡洛法的衛(wèi)星分離姿態(tài)分析方法。飛行試驗(yàn)結(jié)果表明，由該方法得到的分離角速度與常規(guī)保守算法相比，計(jì)算精度提高35%左右，解決了工程中由于采用將各個(gè)偏差取最大值進(jìn)行保守計(jì)算而導(dǎo)致計(jì)算結(jié)果偏差較大的問(wèn)題，可為衛(wèi)星設(shè)計(jì)、彈簧裝置設(shè)計(jì)提供參考。

衛(wèi)星；分離姿態(tài)；彈簧裝置；蒙特卡洛法；龍格庫(kù)塔法

0 引 言

隨著計(jì)算機(jī)技術(shù)、微納技術(shù)及微機(jī)電系統(tǒng)技術(shù)的迅速發(fā)展，質(zhì)量在10~100 kg范圍的微納衛(wèi)星以其研制周期短、低成本、高性能等優(yōu)勢(shì)，在通信、軍事、遙感與對(duì)地測(cè)繪、空間科學(xué)試驗(yàn)等領(lǐng)域表現(xiàn)出較好的發(fā)展前景[1]。

微納衛(wèi)星通常以搭載或一箭多星的方式進(jìn)行發(fā)射，以彈簧裝置作為分離動(dòng)力源，以一定的分離速度、分離角速度與運(yùn)載火箭進(jìn)行分離[2]。衛(wèi)星分離過(guò)程中，除了不能與火箭發(fā)生碰撞，還必須將分離角速度控制在允許范圍內(nèi)，過(guò)大的分離角速度甚至?xí)?dǎo)致衛(wèi)星無(wú)法控制姿態(tài)，直接影響飛行任務(wù)成敗。

衛(wèi)星分離姿態(tài)計(jì)算方法主要包括兩種：數(shù)值計(jì)算方法和基于ADAMS軟件的仿真計(jì)算方法。

數(shù)值計(jì)算方法包括兩種：a）建立常質(zhì)量航天器動(dòng)力學(xué)方程，并采用龍格庫(kù)塔方法進(jìn)行數(shù)值求解，如：付碧紅[3]、盧麗穎[4]分別以搭載星、子星為研究對(duì)象，建立并求解了三自由度動(dòng)力學(xué)方程，Jeyakumar，BN Rao[5]也通過(guò)龍格庫(kù)塔法得到了分離體的非線性常微分方程數(shù)值解；b）建立拉格朗日動(dòng)力學(xué)方程，如：王鑫等[6]采用拉格朗日動(dòng)力學(xué)原理建立了分離體姿態(tài)動(dòng)力學(xué)模型，探索了分離體姿態(tài)動(dòng)力學(xué)方程組數(shù)值求解方法。而隨著計(jì)算機(jī)仿真技術(shù)的發(fā)展，眾多學(xué)者[7~11]提出了采用ADAMS軟件進(jìn)行了衛(wèi)星分離仿真分析，ADAMS軟件基于多剛體動(dòng)力學(xué)方程，采用自動(dòng)變階、變步長(zhǎng)方法進(jìn)行計(jì)算[12]。

但相比上述計(jì)算模型，實(shí)際情況中存在兩類偏差：一類是衛(wèi)星重量偏差、質(zhì)心偏差及轉(zhuǎn)動(dòng)慣量偏差；另一類是彈簧裝置推力偏斜、推力大小偏差及安裝誤差。這些偏差因素都將影響衛(wèi)星的分離姿態(tài)。

在工程應(yīng)用中，為考慮以上偏差因素對(duì)分離姿態(tài)的影響，通常采用的方法是將每一個(gè)偏差量都取到最大，計(jì)算極限的分離姿態(tài)。然而，用所有偏差因素的極限值計(jì)算得到的極限分離姿態(tài)產(chǎn)生概率非常小，遠(yuǎn)遠(yuǎn)超出了航天可靠性要求。因此，使用以上方法計(jì)算結(jié)果偏于保守。

使用以上保守的計(jì)算方法時(shí)，為了滿足衛(wèi)星分離指標(biāo)要求，通常采取兩種措施：一種是減小以上因素的設(shè)計(jì)偏差值；另一種是提高彈簧裝置生產(chǎn)精度，并大量投產(chǎn)，選配出彈簧推力偏差較小的彈簧裝置。以上方法在限制了衛(wèi)星與彈簧裝置結(jié)構(gòu)設(shè)計(jì)的同時(shí)也提高了產(chǎn)品生產(chǎn)成本。

因此，本文以采用彈簧裝置進(jìn)行分離的微納衛(wèi)星為例，從概率分布的角度，提出基于蒙特卡洛法的考慮設(shè)計(jì)偏差影響的衛(wèi)星分離姿態(tài)計(jì)算方法，并通過(guò)飛行試驗(yàn)驗(yàn)證了該方法的有效性。

1 動(dòng)力學(xué)方程建立與求解

1.1 建立動(dòng)力學(xué)方程

以CZ-4B運(yùn)載火箭搭載發(fā)射微納衛(wèi)星（X衛(wèi)星）為例，建立六自由度常質(zhì)量航天器動(dòng)力學(xué)方程。該衛(wèi)星在4個(gè)彈簧裝置作用下與運(yùn)載火箭進(jìn)行分離，分離方案如圖1所示。

圖1 衛(wèi)星分離方案

由于衛(wèi)星質(zhì)量、轉(zhuǎn)動(dòng)慣量遠(yuǎn)遠(yuǎn)小于運(yùn)載火箭末子級(jí)質(zhì)量、轉(zhuǎn)動(dòng)慣量，因此分離過(guò)程中，不考慮運(yùn)載火箭對(duì)衛(wèi)星分離姿態(tài)的影響。

將坐標(biāo)系設(shè)定于星箭分離面中心，固定在運(yùn)載火箭上，以圖1模型分離系統(tǒng)為研究對(duì)象，建立4個(gè)彈簧裝置作用下的分離姿態(tài)計(jì)算分析模型，如圖2所示。

圖2 分離姿態(tài)計(jì)算模型

Δ，Δ，Δ—衛(wèi)星質(zhì)心在，，軸的偏移量；1，2，3，4—考慮安裝偏差后的4個(gè)彈簧的安裝距離，其中3，4為負(fù)值；F第個(gè)彈簧裝置推力在向分力（=1,2,3,4）；F第個(gè)彈簧裝置推力在向分力；F第個(gè)彈簧裝置推力在向分力

考慮彈簧推力偏斜影響，以單個(gè)彈簧為研究對(duì)象，根據(jù)總體坐標(biāo)系，將每一個(gè)彈簧裝置的推力進(jìn)行分解，如圖3所示。

圖3 彈簧推力分解示意

F—第個(gè)彈簧裝置推力；—第個(gè)彈簧推力與軸夾角；—第個(gè)彈簧力在平面投影力與軸夾角

根據(jù)圖3可以得到每個(gè)彈簧推力在坐標(biāo)方向的分解力為

根據(jù)圖2和式（1）可以得到使衛(wèi)星繞三軸旋轉(zhuǎn)力矩M，M，M的計(jì)算公式為

根據(jù)牛頓第二定律和動(dòng)量矩定理，建立衛(wèi)星與運(yùn)載火箭末子級(jí)兩個(gè)剛體之間的六自由度動(dòng)力學(xué)方程：

1.2 動(dòng)力學(xué)方程求解

上述方程維數(shù)高、變量多、且互相耦合，無(wú)法求解其解析解，可采用龍格庫(kù)塔方法進(jìn)行數(shù)值計(jì)算求解。

在計(jì)算過(guò)程中，彈簧推力、力矩為變化值，為了解決這一問(wèn)題，本文提出將彈簧工作時(shí)間劃分為多個(gè)連續(xù)子區(qū)間，在每一個(gè)子區(qū)間內(nèi)，將該區(qū)間中間時(shí)刻的彈簧力作為該子區(qū)間內(nèi)的彈簧推力。

首先計(jì)算彈簧裝置工作時(shí)間，以4個(gè)彈簧組成的分離系統(tǒng)為研究對(duì)象，由于彈簧質(zhì)量遠(yuǎn)小于衛(wèi)星質(zhì)量，衛(wèi)星運(yùn)動(dòng)為彈簧振子簡(jiǎn)諧運(yùn)動(dòng)。

在初始時(shí)刻0=0，初始速度為0，彈簧變形量達(dá)到最大值max，max即為彈簧振子振幅；在彈簧工作結(jié)束時(shí)刻s，彈簧變形量降低至最小值min：

式中為分離系統(tǒng)簡(jiǎn)諧運(yùn)動(dòng)角頻率。

根據(jù)以上公式計(jì)算彈簧工作時(shí)間為

式中k為第個(gè)彈簧的剛度。

式中為第個(gè)彈簧力角頻率；為彈簧效率，根據(jù)工程經(jīng)驗(yàn)取值。

本文基于MATLAB軟件，采用龍格庫(kù)塔數(shù)值計(jì)算方法來(lái)求解公式，求解過(guò)程如圖4所示。

圖4 動(dòng)力學(xué)方程求解流程

1.3 計(jì)算方法驗(yàn)證

為了驗(yàn)證本文計(jì)算方法的正確性，針對(duì)同一組計(jì)算參數(shù)，采用該方法和ADAMS仿真計(jì)算兩種方法進(jìn)行計(jì)算，并將計(jì)算結(jié)果進(jìn)行比對(duì)。

參照X衛(wèi)星給出一組設(shè)計(jì)參數(shù)，如表1所示，進(jìn)行衛(wèi)星分離姿態(tài)計(jì)算。

表1 計(jì)算參數(shù)匯總

續(xù)表1

4#彈簧力/NFmax=250；Fmin=80 1#彈簧位置/mm107.05 2#彈簧位置/mm107.08 3#彈簧位置/mm106.99 4#彈簧位置/mm107.0 1#彈簧推力角度/(°)θ1=0.762；φ1=30 2#彈簧推力角度/(°)θ2=0.742；φ2=45 3#彈簧推力角度/(°)θ3=0.758；φ3=60 4#彈簧推力角度/(°)θ4=0.750；φ4=90 彈簧效率η0.85

注：max—彈簧裝置最大推力；min—彈簧裝置最小推力

根據(jù)式（1）計(jì)算得到彈簧工作時(shí)間為0.074 2 s，再將該時(shí)間劃分為50份，根據(jù)1.2節(jié)公式，采用MATLAB軟件對(duì)衛(wèi)星分離姿態(tài)進(jìn)行計(jì)算，計(jì)算結(jié)果以分離角速度為例，如圖5和表2所示。

圖5 數(shù)值計(jì)算結(jié)果

表2 數(shù)值計(jì)算結(jié)果匯總

在ADAMS軟件中，建立仿真分析模型，衛(wèi)星的質(zhì)量、慣量和質(zhì)心位置同表1。在地面建立總坐標(biāo)系，在彈簧推力作用點(diǎn)施加彈簧推力，該采用行程與剛度函數(shù)表示。

設(shè)置計(jì)算時(shí)間為0.1 s、載荷步數(shù)量為50，進(jìn)行計(jì)算，結(jié)果如圖6、表3所示。

圖6 ADAMS姿態(tài)角速度結(jié)果

表3 仿真計(jì)算結(jié)果匯總

將由兩種方法得到的計(jì)算結(jié)果進(jìn)行對(duì)比，得到兩種方法相對(duì)差如表4所示。

表4 計(jì)算偏差匯總

經(jīng)過(guò)表4對(duì)比可得，兩種計(jì)算方法得到的計(jì)算結(jié)果最大相對(duì)差為0.77%，這是由于將子區(qū)間中間時(shí)刻彈簧推力簡(jiǎn)化為該子區(qū)間彈簧推力而造成的。本算例證明了1.2節(jié)計(jì)算方法的正確性。

2 考慮設(shè)計(jì)偏差的分離姿態(tài)計(jì)算方法

2.1 計(jì)算方法

蒙特卡洛方法不同于確定性的數(shù)值方法，它是通過(guò)模擬隨機(jī)變量來(lái)解決數(shù)學(xué)問(wèn)題的非確定性的（概率統(tǒng)計(jì)的或隨機(jī)的）數(shù)值方法[13]。

因此，本文提出基于蒙特卡洛方法的衛(wèi)星分離姿態(tài)評(píng)估方法：首先根據(jù)各個(gè)設(shè)計(jì)變量的偏差范圍隨機(jī)產(chǎn)生大量的設(shè)計(jì)參數(shù)；再根據(jù)第1節(jié)方法進(jìn)行衛(wèi)星分離姿態(tài)計(jì)算，將得到大量衛(wèi)星分離姿態(tài)計(jì)算結(jié)果；再對(duì)這些計(jì)算結(jié)果進(jìn)行統(tǒng)計(jì)，得到其分布規(guī)律；最后根據(jù)概率分布范圍，找出合理的極限分離姿態(tài)。分析方法流程如圖7所示。

圖7 衛(wèi)星分離姿態(tài)分析方法

2.2 計(jì)算方法驗(yàn)證

以CZ-4B運(yùn)載火箭搭載發(fā)射X衛(wèi)星為例，驗(yàn)證該方法有效性。

衛(wèi)星、彈簧裝置完成生產(chǎn)后，對(duì)衛(wèi)星的質(zhì)量、質(zhì)心位置、轉(zhuǎn)動(dòng)慣量進(jìn)行實(shí)測(cè)；對(duì)彈簧裝置的推力大小、安裝位置進(jìn)行實(shí)測(cè)，如表5所示。

表5 產(chǎn)品實(shí)測(cè)參數(shù)

由于彈簧裝置在設(shè)計(jì)過(guò)程中為防止彈簧推桿與彈簧殼體發(fā)生卡滯，兩者之間為間隙配合，因此彈簧推力是在一定角度范圍偏斜。根據(jù)彈簧推桿、彈簧殼體的加工參數(shù)，計(jì)算得到彈簧最大推力偏斜為0.622°。

針對(duì)這一項(xiàng)設(shè)計(jì)偏差，將采用傳統(tǒng)保守算法以及第2.1節(jié)提出的算法對(duì)X衛(wèi)星的分離姿態(tài)進(jìn)行計(jì)算，并將兩種方法的計(jì)算結(jié)果與飛行試驗(yàn)結(jié)果進(jìn)行對(duì)比。

按照保守算法，將推力偏斜設(shè)置為最大值，確定衛(wèi)星分離最嚴(yán)酷工況，如表6所示。

表6 保守算法計(jì)算參數(shù)

根據(jù)表5、表6的設(shè)計(jì)參數(shù)，在ADAMS軟件中進(jìn)行衛(wèi)星分離姿態(tài)仿真分析，建模方法同1.3節(jié)，以分離角速度為例，計(jì)算結(jié)果為：2.892 (°)/s，-3.0 (°)/s，0.970 (°)/s。

根據(jù)2.1節(jié)提出的評(píng)估方法，設(shè)置彈簧推力偏斜角度為：∈[0,0.622°]，∈[0,360°]，在MATLAB計(jì)算程序中產(chǎn)生100 000個(gè)隨機(jī)變量，通過(guò)計(jì)算得到100 000組計(jì)算結(jié)果，以分離角速度為例，對(duì)這些計(jì)算結(jié)果進(jìn)行處理，通過(guò)程序判斷計(jì)算結(jié)果服從正態(tài)分布規(guī)律，如圖8至圖10所示，并得到繞三軸分離角速度均值與方差，如表7所示。

圖8 繞X軸分離角速度分析結(jié)果

圖9 繞Y軸分離角速度分析結(jié)果

圖10 繞Z軸分離角速度分析結(jié)果

表7 分離角速度均值與方差

據(jù)表7結(jié)果，按照3原則，得到置信度為99.73%的衛(wèi)星繞三軸分離角速度：1.822 (°)/s，-1.940 (°)/s，0.594 (°)/s。

2014年，該衛(wèi)星通過(guò)CZ-4B運(yùn)載火箭搭載發(fā)射成功，根據(jù)衛(wèi)星回傳數(shù)據(jù)，星箭分離時(shí)刻，衛(wèi)星繞三軸分離角速度分別為：0.8 (°)/s，-0.6 (°)/s，0.17 (°)/s。

對(duì)以上數(shù)據(jù)進(jìn)行匯總，如表8所示。

表8 計(jì)算、試驗(yàn)結(jié)果數(shù)據(jù)匯總

對(duì)比表8數(shù)據(jù)可得，由第2.1節(jié)方法得到的計(jì)算結(jié)果與常規(guī)保守法計(jì)算結(jié)果相比，計(jì)算精度分別提高37%，35%，39%，該方法可以代替?zhèn)鹘y(tǒng)的保守算法。

綜合表7和飛行結(jié)果，如果按照1σ原則，得到置信度為68.27%的分析結(jié)果，0.7 (°)/s、-0.806 (°)/s、0.198 (°)/s，不能包絡(luò)飛行試驗(yàn)結(jié)果；如果按照2σ原則，得到置信度為95.45%的分析結(jié)果，1.261 (°)/s、-1.373 (°)/s、0.396 (°)/s，可以包絡(luò)飛行試驗(yàn)結(jié)果；說(shuō)明至少按照2σ原則給出的計(jì)算結(jié)果才可以滿足使用要求，因此本文按照3σ原則給出的計(jì)算結(jié)果。

3 結(jié) 論

本文建立了考慮推力偏斜的六自由度動(dòng)力學(xué)方程，采用龍格庫(kù)塔方法求解該方程，算例結(jié)果表明，由該方法得到的分離角速度結(jié)果與仿真計(jì)算結(jié)果比最大相差為0.77%，驗(yàn)證了該方法正確性。

基于以上研究，提出基于蒙特卡洛法的考慮設(shè)計(jì)偏差影響的衛(wèi)星分離姿態(tài)計(jì)算方法。算例結(jié)果表明，衛(wèi)星分離角速度計(jì)算結(jié)果服從正態(tài)分布規(guī)律。飛行試驗(yàn)結(jié)果表明，要至少按照2σ原則給出的計(jì)算結(jié)果才可以滿足使用要求，按照本文提出的3σ原則給出的計(jì)算結(jié)果與常規(guī)保守算法相比其計(jì)算精度提高35%左右。

該方法解決了工程中由于采用將各個(gè)偏差取最大值進(jìn)行計(jì)算而導(dǎo)致計(jì)算結(jié)果偏差較大的問(wèn)題,該方法可為衛(wèi)星設(shè)計(jì)、彈簧裝置設(shè)計(jì)提供參考。

[1] 石榮, 李瀟, 鄧科. 微納衛(wèi)星發(fā)展現(xiàn)狀及在光學(xué)成像偵察中的應(yīng)用[J]. 航天電子對(duì)抗, 2016: 32(1): 8-13.

[2] Liu L K, Zheng G T. Parameter analysis of PAF for whole-spacecraft vibration isolation[J]. Aerospace Science and Technology, 2007(11): 464-472.

[3] 付碧紅, 杜光華. 搭載星與運(yùn)載火箭分離的動(dòng)力學(xué)研究[J]. 飛行力學(xué), 2006, 24(1): 55-58.

[4] 盧麗穎, 孟憲紅, 邢依琳. 衛(wèi)星空間分離動(dòng)力學(xué)研究[J]. 動(dòng)力學(xué)與控制學(xué)報(bào)，2014, 12(2): 165-169.

[5] Jeyakumar D, Rao B N. Dynamics of satellite separation system[J]. Journal of Sound and Vibration, 2006, 297(1-2): 444-455.

[6] 王鑫, 袁曉光, 楊星. 基于拉格朗日方法的飛行器多體分離姿態(tài)動(dòng)力學(xué)分析研究[J]. 西北工業(yè)大學(xué)學(xué)報(bào): 2014, 32(1): 18-22.

[7] 王秋梅, 孟憲紅, 楊慶成. 衛(wèi)星二次分離方案仿真研究[J]. 系統(tǒng)仿真學(xué)報(bào), 2010, 22(9): 2217-2222.

[8] 沈曉鳳, 肖余之, 杜三虎, 等. 基于蒙特卡羅方法的小衛(wèi)星偏心分離動(dòng)力學(xué)分析[J]. 上海航天, 2014, 31(1): 12-17.

[9] 沈曉鳳, 肖余之, 康志宇. 小衛(wèi)星偏心分離動(dòng)力學(xué)仿真模型的建立與驗(yàn)證[J]. 飛行力學(xué), 2012, 30(3): 258-262.

[10] 曾文花, 吳琳娜. 多星與上面級(jí)非對(duì)稱分離技術(shù)研究[J]. 上海航天, 2014, 31(2): 30-36.

[11] 張兵, 岑拯. 多星分離的ADAMS仿真[J]. 導(dǎo)彈與航天運(yùn)載技術(shù), 2004(2): 1-6.

[12] 胡星志. 小衛(wèi)星星箭分離系統(tǒng)設(shè)計(jì)、分析與優(yōu)化研究[D]. 長(zhǎng)沙: 國(guó)防科技大學(xué), 2012.

[13] 金暢. 蒙特卡洛方法中隨機(jī)數(shù)發(fā)生器和隨機(jī)抽樣方法的研究[D]. 大連: 大連理工大學(xué), 2005.

Calculation Method of Satellite Separation Attitude Based on Monte Carlo Method

Kang Shi-peng1, Jiang Tao1, Geng Sheng-wei2, Wang Jin-tong1

(1. Aerospace System Engineering Shanghai, Shanghai, 201109; 2. College of Aerospace Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing, 201106）

In order to calculate the separation attitude of the satellite which is separated by spring device, the 6 degrees of freedom dynamic equation is established and solved by the Runge-Kutta Method. The numerical examples indicate that the maximum deviation of separation angular velocity obtained by this method and the ADAMS simulation is less than 0.77%. Based on the above studies, a analysis method based on Monte Carlo Method is proposed to calculate satellite separation attitude considering the satellite and spring device design deviations. The results of flight test indicate that the proposed method can improve the accuracy by 35% compared with the conventional conservative algorithm. The method solves the problem that the deviation of the calculation result is large due to the conservative calculation of the maximum value of each deviation. This method can provide a reference for satellite design and spring device design.

Satellite; Separation attitude; Spring device; Monte Carlo method; Runge-Kutta method

1004-7182(2017)06-0036-06

10.7654/j.issn.1004-7182.20170609

V412.4+2

A

2017-06-04；

2017-07-11

康士朋（1983-），男，博士研究生，工程師，主要研究方向?yàn)榧w結(jié)構(gòu)和星箭連接分離裝置設(shè)計(jì)