龔秋武 陳鼎 李翔 黃宇

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

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

0 引言

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

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

本文在質(zhì)心運動的基礎上，考慮火箭子級殘骸的姿態(tài)運動，分析地球引力攝動以及氣動力的影響，推導建立質(zhì)心運動與姿態(tài)運動方程，對火箭子級殘骸進行落點預報。

1 運動過程及受力分析

1.1 運動過程

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

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

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

1.2 受力分析

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

1.2.1 地球引力

由于地球是非標準球體，為提高地球引力的描述精度，須將地球引力分為中心引力和非球形攝動力兩部分。其中，中心引力的加速度為

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

地球非球形引力攝動加速度為地球非球形引力位的梯度，即

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

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

1.2.2 大氣作用力

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

（1）氣動力

殘骸受到的氣動力如圖1所示。

圖1 殘骸所受氣動力Fig.1 Aerodynamic force of rocket wreck age

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

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

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

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

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

（2）氣動力矩

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

圖2 殘骸所受氣動力矩Fig.2 Moment of aerodynamic force

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

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

2 動力學建模

2.1 質(zhì)心動力學模型

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

圖3 質(zhì)心運動Fig.3 Motion of barycenter

根據(jù)運動學知識，對位置矢量求絕對導數(shù)即得到速度矢量，對速度矢量求絕對導數(shù)即得到加速度矢量a，整理得到

式中ω為兩坐標系之間的相對轉動角速度矢量。

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

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

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

2.2 姿態(tài)動力學模型

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

圖4 殘骸姿態(tài)運動示意Fig.4 Attitude motion of racket wreckage

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

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

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

選定彈體坐標系為計算坐標系，將矢量形式的姿態(tài)動力學方程轉換為可積分運算的標量形式方程組。

3 數(shù)值仿真

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

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

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

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

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

4 結束語

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