龔秋武 陳鼎 李翔 黃宇

（1 中國酒泉衛(wèi)星發(fā)射中心，酒泉 732750）

（2 中國人民解放軍95983部隊(duì)，酒泉 732750）

0 引言

航天發(fā)射任務(wù)中，火箭子級(jí)、整流罩等殘骸的搜索回收關(guān)系到地面人員、建筑和設(shè)備等的安全，必須嚴(yán)加防控，嚴(yán)密組織，嚴(yán)格執(zhí)行，確保不傷人員，不損建筑，不毀設(shè)備。

而殘骸搜索任務(wù)的前提就是要準(zhǔn)確計(jì)算殘骸的落點(diǎn)，以便在發(fā)射任務(wù)實(shí)施前做好危險(xiǎn)區(qū)防控和人員疏散，任務(wù)實(shí)施時(shí)引導(dǎo)搜索人員和車輛等前往安全目標(biāo)點(diǎn)。關(guān)于落點(diǎn)預(yù)報(bào)，國內(nèi)外學(xué)者做了許多研究，如：文獻(xiàn)[1-3]在載人飛船返回艙落點(diǎn)預(yù)報(bào)理論的研究中，考慮了地球自轉(zhuǎn)、高空風(fēng)以及其他攝動(dòng)因素對(duì)落點(diǎn)精度的影響；文獻(xiàn)[4-9]對(duì)彈道導(dǎo)彈落點(diǎn)預(yù)報(bào)進(jìn)行了研究，包括軌道外推、氣動(dòng)力系數(shù)辨識(shí)、地球引力攝動(dòng)等計(jì)算方法；文獻(xiàn)[10-17]研究了火箭殘骸的落點(diǎn)預(yù)報(bào)與控制算法；文獻(xiàn)[18-20]研究了濾波算法在殘骸落點(diǎn)預(yù)報(bào)上的應(yīng)用。上述這些研究，都是基于質(zhì)點(diǎn)的動(dòng)力學(xué)模型，并未考慮返回器的姿態(tài)運(yùn)動(dòng)，而對(duì)于高空墜落地面的物體，姿態(tài)變換會(huì)影響其迎風(fēng)面，改變其所受的氣動(dòng)力，進(jìn)而干擾質(zhì)心運(yùn)動(dòng)，影響最后的落點(diǎn)精度。

本文在質(zhì)心運(yùn)動(dòng)的基礎(chǔ)上，考慮火箭子級(jí)殘骸的姿態(tài)運(yùn)動(dòng)，分析地球引力攝動(dòng)以及氣動(dòng)力的影響，推導(dǎo)建立質(zhì)心運(yùn)動(dòng)與姿態(tài)運(yùn)動(dòng)方程，對(duì)火箭子級(jí)殘骸進(jìn)行落點(diǎn)預(yù)報(bào)。

1 運(yùn)動(dòng)過程及受力分析

1.1 運(yùn)動(dòng)過程

火箭子級(jí)在完成推力任務(wù)時(shí)，在指令控制下與火箭主體或者載荷分離，由于飛行程序的不同，各子級(jí)分離時(shí)的位置與姿態(tài)也各不相同，但是基本的物理過程卻是一致的。

一般火箭一、二子級(jí)和整流罩分離時(shí)的飛行高度較低，飛行速度還未達(dá)到繞地飛行最小速度，所以會(huì)在分離后很短的時(shí)間內(nèi)墜落地球?；鸺蛹?jí)和更上面的部分在分離時(shí)已經(jīng)進(jìn)入空間環(huán)境，速度也達(dá)到了繞地飛行的基本條件，因而這些部分會(huì)在太空繞地飛行一段時(shí)間，由于它們處于無控狀態(tài)，不能進(jìn)行軌道保持和控制，一般繞地飛行數(shù)圈后便墜入大氣層被燒蝕掉。

墜落地球表面的火箭一、二子級(jí)和整流罩的著陸過程大致分為三個(gè)階段：1）減速上升段，此階段子級(jí)借助分離時(shí)向上的速度分量繼續(xù)向上飛行，由于受到地球引力和大氣作用力的影響，向上的速度分量逐漸減小至零，因而做減速上升運(yùn)動(dòng)；2）加速下降段，此階段在地球引力的作用下，子級(jí)加速下墜，大氣作用力隨著子級(jí)速度的不斷增大而增大；3）平衡段，此階段子級(jí)受到的地球引力與大氣作用力基本平衡，子級(jí)以一個(gè)相對(duì)穩(wěn)定的速度下墜直至墜落地面。

1.2 受力分析

殘骸在下墜過程中，一般處于無控狀態(tài)，主要受到地球引力和大氣作用力的影響。本節(jié)即對(duì)子級(jí)殘骸受到的兩種作用力進(jìn)行分析。

1.2.1 地球引力

由于地球是非標(biāo)準(zhǔn)球體，為提高地球引力的描述精度，須將地球引力分為中心引力和非球形攝動(dòng)力兩部分。其中，中心引力的加速度為

式中 GME為地球引力常量；r為殘骸質(zhì)心的地心系位置矢量，r為地心距，其中X、Y、Z為地心系位置坐標(biāo)。

地球非球形引力攝動(dòng)加速度為地球非球形引力位的梯度，即

式中U為地球引力位，為締合勒讓德多項(xiàng)式，φ為目標(biāo)點(diǎn)位的緯度，λ為目標(biāo)點(diǎn)位的經(jīng)度，為引力位田諧系數(shù)，RE為地球平均半徑。

從而地球引力加速度g可表示為中心引力加速度g0與地球非球形引力攝動(dòng)加速度g0′兩部分之和，即

1.2.2 大氣作用力

火箭殘骸在下墜過程中，大氣作用力影響很大，一方面氣動(dòng)力影響殘骸質(zhì)心運(yùn)動(dòng)，另一方面氣動(dòng)力產(chǎn)生的氣動(dòng)力矩影響殘骸姿態(tài)運(yùn)動(dòng)。

（1）氣動(dòng)力

殘骸受到的氣動(dòng)力如圖1所示。

圖1 殘骸所受氣動(dòng)力Fig.1 Aerodynamic force of rocket wreck age

圖1中O1-x1y1z1為彈體坐標(biāo)系，O1-xvyvzv為速度坐標(biāo)系，f為氣動(dòng)力合力，O1為殘骸質(zhì)心，Op為壓心。為便于計(jì)算，可將氣動(dòng)力分解到彈體坐標(biāo)系或者速度坐標(biāo)系。當(dāng)分解到彈體坐標(biāo)系時(shí)，可得到軸向力fx1、法向力fy1和橫向力fz1；當(dāng)分解到速度坐標(biāo)系時(shí)，可得到阻力fD、升力fL和側(cè)力fC。根據(jù)氣動(dòng)力與飛行速度、大氣密度等參數(shù)的關(guān)系，可計(jì)算得到氣動(dòng)力的大小為

式中0ρ為大氣密度；Sr為殘骸參考面積，一般取截面積；C為氣動(dòng)力系數(shù)；v為殘骸相對(duì)于大氣的速度，而非自身飛行速度，這里主要考慮高空風(fēng)的影響。一般地，高空風(fēng)數(shù)據(jù)通過氣象測站釋放的探空氣球測得，數(shù)據(jù)格式為對(duì)應(yīng)高度處的風(fēng)向和風(fēng)速，風(fēng)向以當(dāng)?shù)卣睘?°，順時(shí)針方向至360°。假設(shè)某高度處的高空風(fēng)風(fēng)向?yàn)棣罻，風(fēng)速為vW，其在測站坐標(biāo)系中的速度矢量為

經(jīng)過坐標(biāo)轉(zhuǎn)換，可將式（5）所示的高空風(fēng)在測站坐標(biāo)系的速度轉(zhuǎn)換到發(fā)射坐標(biāo)系，進(jìn)而得到發(fā)射系中殘骸相對(duì)于大氣的速度，即

式中vG為殘骸在發(fā)射系中的速度矢量；GE為地心坐標(biāo)系到發(fā)射坐標(biāo)系的轉(zhuǎn)換矩陣；EC為測站坐標(biāo)系到地心坐標(biāo)系的轉(zhuǎn)換矩陣。以上坐標(biāo)系的定義可參考文獻(xiàn)[21-23]，轉(zhuǎn)換矩陣GE和EC的計(jì)算公式分別為：

其中，0α為發(fā)射方位角；0λ為發(fā)射點(diǎn)經(jīng)度；0φ為發(fā)射點(diǎn)緯度；M1、M2和M3為初等轉(zhuǎn)換矩陣，具體定義和表達(dá)式可參考文獻(xiàn)[21]；cλ為測站點(diǎn)位的經(jīng)度；cφ為測站點(diǎn)位的緯度。

（2）氣動(dòng)力矩

由于氣動(dòng)力作用點(diǎn)在壓心Op，而壓心一般不與質(zhì)心重合，所以火箭殘骸受到的氣動(dòng)力會(huì)對(duì)殘骸質(zhì)心產(chǎn)生氣動(dòng)力矩，如圖2所示。圖2中rp為壓心Op的位置矢量。

圖2 殘骸所受氣動(dòng)力矩Fig.2 Moment of aerodynamic force

假設(shè)彈體為靜穩(wěn)定的，用氣動(dòng)力矩系數(shù)來表征氣動(dòng)力矩，進(jìn)而可將氣動(dòng)力矩表示為

式中My1st為偏航力矩；Mz1st為俯仰力矩；my1為偏航力矩系數(shù)；mz1為俯仰力矩系數(shù)；q為動(dòng)壓，l為殘骸特征長度；α為攻角；β為側(cè)滑角。由于壓心與質(zhì)心一般均處于彈體軸線上，所以氣動(dòng)力不產(chǎn)生滾轉(zhuǎn)力矩，即Mx1st=0。

2 動(dòng)力學(xué)建模

2.1 質(zhì)心動(dòng)力學(xué)模型

殘骸在墜落過程中任意時(shí)刻的質(zhì)心運(yùn)動(dòng)狀態(tài)如圖3所示，其中OE-XEYEZE為地心赤道坐標(biāo)系（簡稱地心系），OE為地心，OG-xyz為發(fā)射坐標(biāo)系，OG為發(fā)射點(diǎn)，Rwc為殘骸質(zhì)心在發(fā)射系中的位置矢量、Rfp為發(fā)射點(diǎn)的地心系位置矢量，殘骸質(zhì)心的地心系位置矢量r=Rfp+Rwc。

圖3 質(zhì)心運(yùn)動(dòng)Fig.3 Motion of barycenter

根據(jù)運(yùn)動(dòng)學(xué)知識(shí)，對(duì)位置矢量求絕對(duì)導(dǎo)數(shù)即得到速度矢量，對(duì)速度矢量求絕對(duì)導(dǎo)數(shù)即得到加速度矢量a，整理得到

式中ω為兩坐標(biāo)系之間的相對(duì)轉(zhuǎn)動(dòng)角速度矢量。

假設(shè)殘骸質(zhì)量為m，結(jié)合前面的受力分析，則根據(jù)牛頓第二定律，可得到矢量形式的質(zhì)心動(dòng)力學(xué)方程

式中f為氣動(dòng)力；Fe為牽連力，為科氏力，F(xiàn)k=-mω×（ω×r）。

矢量形式的動(dòng)力學(xué)方程無法進(jìn)行積分運(yùn)算，還需選定計(jì)算坐標(biāo)系，將其轉(zhuǎn)換為標(biāo)量形式的方程，此處以發(fā)射坐標(biāo)系作為計(jì)算坐標(biāo)系，將加速度和各力向各軸進(jìn)行分解后便可得到標(biāo)量形式的方程組。

2.2 姿態(tài)動(dòng)力學(xué)模型

如前所述，殘骸墜落過程的質(zhì)心動(dòng)力學(xué)方程中，含有歐拉角參數(shù)，并且將矢量形式的動(dòng)力學(xué)方程轉(zhuǎn)換為標(biāo)量形式的動(dòng)力學(xué)方程組時(shí)，涉及到力、加速度的分解，同樣也含有歐拉角參數(shù)，而質(zhì)心動(dòng)力學(xué)方程是無法解算角度信息的，所以除了分析質(zhì)心動(dòng)力學(xué)外，還要考慮姿態(tài)動(dòng)力學(xué)問題。如圖4所示，O1-x0y0z0為慣性系，彈體系O1-x1y1z1相對(duì)于慣性系的姿態(tài)角分別為俯仰角φ、偏航角ψ和滾轉(zhuǎn)角γ，對(duì)應(yīng)的姿態(tài)角速度為φ˙、ψ˙和γ˙。

圖4 殘骸姿態(tài)運(yùn)動(dòng)示意Fig.4 Attitude motion of racket wreckage

根據(jù)姿態(tài)動(dòng)力學(xué)理論，可得到姿態(tài)角速度與姿態(tài)角之間的關(guān)系，以及矢量形式的姿態(tài)動(dòng)力學(xué)方程，即

其中，ωT為體坐標(biāo)系下的姿態(tài)角速度矢量，為殘骸的慣量張量（假定殘骸為軸對(duì)稱體），其中Ixx、Iyy和Izz分別為殘骸繞體坐標(biāo)系三個(gè)坐標(biāo)軸的轉(zhuǎn)動(dòng)慣量；M為殘骸受到的合外力矩矢量。

由于殘骸處于無控狀態(tài)，控制力矩為零，另外阻尼力矩、附加科氏力矩和附加相對(duì)力矩是小量，忽略不計(jì)，最終殘骸受到的合外力矩只考慮氣動(dòng)力矩，對(duì)于靜穩(wěn)定彈體則為穩(wěn)定力矩。

選定彈體坐標(biāo)系為計(jì)算坐標(biāo)系，將矢量形式的姿態(tài)動(dòng)力學(xué)方程轉(zhuǎn)換為可積分運(yùn)算的標(biāo)量形式方程組。

3 數(shù)值仿真

以某次航天發(fā)射試驗(yàn)任務(wù)為例，對(duì)本文提出的落點(diǎn)預(yù)報(bào)方法進(jìn)行仿真驗(yàn)證，地球參數(shù)參考文獻(xiàn)[21]，發(fā)射點(diǎn)位參數(shù)、一子級(jí)殘骸參數(shù)以及殘骸分離點(diǎn)參數(shù)分別根據(jù)實(shí)際數(shù)據(jù)設(shè)置，高空風(fēng)采取插值得到。將上述相關(guān)參數(shù)帶入到建立的動(dòng)力學(xué)模型中，利用四階 Runge-Kutta積分算法進(jìn)行數(shù)值積分，得到殘骸從分離點(diǎn)到落地的三維彈道軌跡、飛行速度、飛行姿態(tài)角和姿態(tài)角速度等仿真結(jié)果（如圖5～9所示），殘骸預(yù)報(bào)落點(diǎn)與實(shí)際落點(diǎn)誤差見表1。

表1 殘骸預(yù)報(bào)落點(diǎn)與實(shí)際落點(diǎn)誤差Tab.1 The error between the prediction falling point and real falling point

圖5 殘骸三維軌跡Fig.5 Three dimension trajectory of rocket wreckage

圖8 殘骸姿態(tài)角速度變化曲線Fig.8 Attitude rate curves of rocket wreckage

分析圖5所示的仿真軌跡及實(shí)際落點(diǎn)并結(jié)合表1和圖9，可以看出預(yù)報(bào)落點(diǎn)與實(shí)際落點(diǎn)相差很小，在高度 18km位置處進(jìn)行落點(diǎn)預(yù)報(bào)得到的預(yù)報(bào)落點(diǎn)與實(shí)際落點(diǎn)之間的射向偏差約為 1.2km，橫向偏差約為0.1km，射向和橫向偏差均小于2km，以預(yù)報(bào)落點(diǎn)為圓心，2km為半徑，得到的搜索范圍約占工業(yè)部門提供的理論落區(qū)范圍的10%；由圖6可知，殘骸下墜過程中，在地球引力和氣動(dòng)力的共同作用下，殘骸先減速后加速，再減速至相對(duì)穩(wěn)定；觀察圖7～8殘骸姿態(tài)變化的仿真結(jié)果可發(fā)現(xiàn)，殘骸的俯仰角φ逐漸減小至負(fù)方向至穩(wěn)定，偏航角ψ在下墜過程中呈規(guī)律性波動(dòng)，由于未考慮滾轉(zhuǎn)力矩，滾轉(zhuǎn)角γ以初始滾轉(zhuǎn)角速度持續(xù)滾轉(zhuǎn)。還可以看出，因?yàn)闅埡∈庆o穩(wěn)定的，所以姿態(tài)角速度ωTy和ωTz維持在相對(duì)穩(wěn)定的狀態(tài)。

圖6 殘骸速度變化曲線Fig.6 Velocity curves of rocket wreckage

圖7 殘骸姿態(tài)變化曲線Fig.7 Attitude curves of wreckage

圖9 殘骸落點(diǎn)與落區(qū)Fig.9 The falling point and falling area of rocket wreckage

4 結(jié)束語

火箭殘骸的姿態(tài)運(yùn)動(dòng)對(duì)質(zhì)心運(yùn)動(dòng)具有一定的影響，進(jìn)而影響殘骸飛行軌跡和最終落點(diǎn)，因而本文在質(zhì)心運(yùn)動(dòng)的基礎(chǔ)上，對(duì)姿態(tài)運(yùn)動(dòng)進(jìn)行了建模，進(jìn)而提出一種火箭子級(jí)殘骸六自由度落點(diǎn)預(yù)報(bào)方法，通過數(shù)值仿真，計(jì)算得到了某次航天發(fā)射任務(wù)中一子級(jí)殘骸的墜落軌跡和各時(shí)刻的位置、速度和姿態(tài)等狀態(tài)參數(shù)，還計(jì)算得到了殘骸的預(yù)報(bào)落點(diǎn)。對(duì)比實(shí)際落點(diǎn)，預(yù)報(bào)落點(diǎn)的射向和橫向偏差，均在2km以內(nèi)，對(duì)比工業(yè)部門提供的理論落區(qū)，搜索范圍縮小了約90%，增強(qiáng)了搜索的針對(duì)性，提升了搜索效率。