面向亞軌道解體事故的碎片擴(kuò)散分布建模研究

陳萬通, 田書雨, 張巨聯(lián), 劉 慶, 任詩雨

(1. 中國民航大學(xué)民航航班廣域監(jiān)視與安全管控技術(shù)重點(diǎn)實(shí)驗(yàn)室, 天津 300300;2. 中國民航大學(xué)電子信息與自動化學(xué)院, 天津 300300;3. 上海飛機(jī)設(shè)計(jì)研究院, 上海 201109; 4. 上海航天電子技術(shù)研究所, 上海 201109)

0 引 言

近年來,伴隨著航天活動驅(qū)動者由政府為主轉(zhuǎn)向市場為主,商業(yè)航天成為全球航天活動的重要驅(qū)動力,亞軌道飛行的商業(yè)化進(jìn)程發(fā)展迅速。與民用航空飛行方式不同,亞軌道飛行是介于萬有引力和空氣動力雙重作用下的飛行,由于未達(dá)到繞地飛行速度,飛行器相當(dāng)于在太空中做個(gè)拋物線運(yùn)動回到地球(最高點(diǎn)一般高于卡門線)[1-2]。相比于普通的航班,亞軌道飛行時(shí)間會縮短數(shù)倍;相比于飛船,亞軌道飛行器造價(jià)低廉且可重復(fù)使用,因此呈現(xiàn)出巨大的經(jīng)濟(jì)價(jià)值。其中,以美國縮尺復(fù)合體公司與英國維珍銀河公司聯(lián)合研制的“太空船”、美國藍(lán)色起源公司的“山貓”、瑞士空間系統(tǒng)公司的“亞軌道飛機(jī)可復(fù)用”飛行器[3]為代表性的亞軌道飛行器已開展多次亞軌道飛行試驗(yàn)。據(jù)權(quán)威機(jī)構(gòu)預(yù)測:亞軌道交通運(yùn)輸極有可能成為下一代空天旅行的主要方式之一[4]。

亞軌道飛行器以高超音速飛行,在強(qiáng)烈的氣動載荷作用下易發(fā)生解體并產(chǎn)生大量碎片,碎片在大氣環(huán)境中所受的氣動力呈現(xiàn)一定的隨機(jī)性,其分布和落點(diǎn)較難預(yù)測,一旦與民航客機(jī)發(fā)生碰撞,將瞬時(shí)產(chǎn)生災(zāi)難性后果。全世界對亞軌道事故碎片碰撞風(fēng)險(xiǎn)的廣泛關(guān)注始于2003年美國哥倫比亞號航天飛機(jī)失事事件,碎片橫貫美國德克薩斯州1 000 km×40 km的陸地范圍,覆蓋人口21.6萬[5]。2018年6月30日,日本星際科技公司于北海道航天試驗(yàn)基地試射自行研制的Momo2小型亞軌道火箭,升空后不久發(fā)生爆炸;同年11月28日韓國在亞軌道發(fā)射試驗(yàn)中,火箭飛行10 min后部分殘骸意外落入日本沖繩島海域,引發(fā)日本民眾強(qiáng)烈抗議。據(jù)美國聯(lián)邦航空管理局(Federal Aviation Administration, FAA)計(jì)算,任何超過300 g碎片的撞擊都將導(dǎo)致商用飛機(jī)100%的損毀,不足300 g的碎片如果撞擊到較脆弱的機(jī)身區(qū)域,例如飛機(jī)操縱面或燃油箱,也將導(dǎo)致被撞飛機(jī)損毀的嚴(yán)重后果[6]。

為了保障航天器的在軌安全,國內(nèi)外對空間碎片的研究已經(jīng)長達(dá)數(shù)十年,研究內(nèi)容涵蓋了空間碎片的模型和風(fēng)險(xiǎn)分析[7-8]、碰撞預(yù)警與規(guī)避策略[9]、軌道碎片環(huán)境建模[10]等多個(gè)方面。然而,空間碎片理論主要針對外層空間,重點(diǎn)關(guān)注空間碎片對在軌航天器的碰撞損傷,不涉及地面安全風(fēng)險(xiǎn)。繼哥倫比亞號事件之后,飛行器因解體事故產(chǎn)生的亞軌道碎片引起了全世界廣泛關(guān)注。相比于空間碎片的研究,亞軌道碎片與地面人群的安全利益、財(cái)產(chǎn)利益直接相關(guān)。國外用來預(yù)測亞軌道碎片位置的軟件主要包括法國的發(fā)射和再入安全評估工具ELECTRA?[11]、美國國家航空航天局(National Aeronautics and Space Administratio, NASA)的通用實(shí)時(shí)碎片足跡(common real time debris footprint, CRTF)[12]、美國APT公司的碎片風(fēng)險(xiǎn)評估(debris risk assessment, DeBRA)[13]、斯坦福大學(xué)的靶場安全評估工具(range safety assessment tool, RSAT)[14]、FAA開發(fā)的航天飛機(jī)危險(xiǎn)區(qū)域計(jì)算(shuttle hazard area to aircraft calculator, SHAAC)工具[15]。上述軟件涉及專有的源代碼,僅限特殊用戶群體使用。

在亞軌道碎片分布建模方面,文獻(xiàn)[16-18]根據(jù)哥倫比亞號事件的碎片分布統(tǒng)計(jì)數(shù)據(jù),將亞軌道碎片危險(xiǎn)區(qū)近似等效為矩形,矩形長度為解體高度除以1 000,矩形寬度為長度的1/8,但這種近似較為粗糙,缺乏嚴(yán)謹(jǐn)?shù)目茖W(xué)論證,不具有普遍意義;文獻(xiàn)[19]在亞軌道事故DeBRA的可視化仿真系統(tǒng)中,采用了軌跡分析模型(trajectory analysis program, TAP)仿真碎片分布,但該模型是針對航空特技表演中飛機(jī)發(fā)生意外解體事故時(shí),對碎片的散布范圍進(jìn)行預(yù)測,本質(zhì)上不適用于亞軌道解體事故碎片的分布研究;文獻(xiàn)[20-21]利用不同協(xié)方差傳播算法對空間碎片再入解體過程中碎片的傳播軌跡進(jìn)行預(yù)測,但目前只適用于單個(gè)碎片求解,對影響碎片傳播過程的關(guān)鍵因素,如彈道系數(shù)、阻力系數(shù)等缺乏詳細(xì)闡述,碎片預(yù)測模型不能滿足民航對安全性的高要求。

綜上所述,亞軌道碎片分布建模是實(shí)現(xiàn)危險(xiǎn)區(qū)預(yù)測的關(guān)鍵技術(shù),而有效的危險(xiǎn)區(qū)預(yù)測是實(shí)施安全的空中交通管制的前提。本文在協(xié)方差傳播算法的基礎(chǔ)上,對亞軌道飛行器的解體模型進(jìn)行建模,對解體碎片數(shù)量、尺寸分布以及阻力系數(shù)等關(guān)鍵因素進(jìn)行統(tǒng)計(jì)分析,將有效提高碎片危險(xiǎn)區(qū)的預(yù)測準(zhǔn)確度,使空中交通管制人員能夠在碎片下落到達(dá)民航飛行高度之前將受影響的飛機(jī)引導(dǎo)出危險(xiǎn)區(qū)域,為空中交通管制人員提供安全、高效、魯棒的決策支持,解決商業(yè)亞軌道飛行活動和傳統(tǒng)民航之間的潛在沖突,具有良好的應(yīng)用前景。

1 亞軌道碎片運(yùn)動方程

1.1 坐標(biāo)系建立

利用xyz坐標(biāo)系表示飛行器解體瞬間碎片在站心坐標(biāo)系中的位置,坐標(biāo)系的原點(diǎn)在地球表面的(θ0,φ0)處,θ0、φ0分別表示航空器解體瞬間的初始經(jīng)、緯度。x、y軸所構(gòu)成的平面表示在原點(diǎn)處與地球相切的平面,x軸指向東方向,y軸指向北方向,z軸垂直于xy平面指向天頂方向,xyz坐標(biāo)系也可稱為東北天(east-north-up, ENU)坐標(biāo)系。

1.2 單個(gè)碎片擴(kuò)散機(jī)理和運(yùn)動模型

(1)

式中:r=[x,y,z]T和v=[vx,vy,vz]T分別表示碎片在ENU坐標(biāo)系中的位置矢量和速度矢量;ω=[0,ωecosφ0,ωesinφ0]T表示ENU坐標(biāo)系的角速度矢量,其中ωe=7.292×10-5rad/s表示地球自轉(zhuǎn)的平均角速度;e3=[0,0,1]T表示第三標(biāo)準(zhǔn)單位向量;Re=6.372×106m表示地球半徑;重力加速度表示為g=ge[Re/(Re+z)]2,其中g(shù)e=9.81 m/s2;aD表示與大氣密度相關(guān)的瞬時(shí)加速度,可表示為aD(t)=-1/2(ρ|vrel|vrel)/β,其中β表示彈道系數(shù),vrel表示碎片的速度矢量v與風(fēng)矢量w之差,即vrel=v-w;ξ表示由建模不確定性和干擾所引起的隨機(jī)加速度矢量。

2 碎片傳播預(yù)測算法

空中交通管制對算法的計(jì)算效率要求很高,采用傳統(tǒng)蒙特卡羅方法模擬亞軌道解體事故的大量碎片雖然近似程度高但運(yùn)算量巨大,不能滿足空中管制實(shí)時(shí)響應(yīng)的高要求。因此,可通過解析形式的概率分布模型描述亞軌道解體事故碎片的分布特征。

2.1 協(xié)方差傳播算法

定義擴(kuò)充狀態(tài)矢量s=[rT,vT]T,則碎片運(yùn)動方程式(1)可重新表述[23]為

(2)

(3)

2.2 碎片概率橢球

通過引入關(guān)于時(shí)間t的由中心和形狀矩陣參數(shù)化后構(gòu)造的橢球集,將碎片擴(kuò)散特征的描述推廣到四維(four dimensions, 4D),由于時(shí)間變量的無限性,所得到的4D概率約束問題比三維(three dimensions, 3D)問題復(fù)雜得多。因此,解決該問題需要將時(shí)間離散化,并將每個(gè)樣本時(shí)間與一個(gè)橢球相關(guān)聯(lián)[24]。此時(shí),單一時(shí)刻碎片的4D足跡是概率橢球的一部分,而一系列瞬時(shí)的概率橢球構(gòu)成了一個(gè)“碎片通道”,碎片位置概率密度函數(shù)[21]為

(4)

(5)

式中:

(6)

其中,λ∈[-π/2,π/2],α∈[0,2π]。

3 碎片彈道系數(shù)

當(dāng)亞軌道飛行器發(fā)生解體事故時(shí),碎片的運(yùn)動軌跡主要取決于彈道系數(shù)β,計(jì)算公式如下:

(7)

式中:m表示單個(gè)碎片的質(zhì)量;CD表示阻力系數(shù);A表示碎片橫截面積。接下來分別從面質(zhì)比分布和阻力系數(shù)取值兩個(gè)方面來計(jì)算彈道系數(shù)。

3.1 面質(zhì)比分布

亞軌道飛行器再入過程中發(fā)生解體事故時(shí)會產(chǎn)生大量碎片,本文采用廣泛接受的NASA標(biāo)準(zhǔn)解體模型EVOLVE4.0[25],利用該模型可以求解亞軌道碎片的質(zhì)量、面質(zhì)比、速度增量等關(guān)鍵參數(shù),算法流程圖如圖1所示。

碎片的特征尺寸可表示為

(8)

式中:dx、dy、dz分別表示在三軸上的投影尺寸。對于亞軌道飛行器爆炸解體事故,解體碎片中尺寸大于dc的碎片分布數(shù)量關(guān)系可表示為

(9)

式中:cs的取值與爆炸類型密切相關(guān),對于亞軌道飛行器解體事故來說,cs通常取值為1[26]。根據(jù)式(9)可將亞軌道碎片的尺寸分布函數(shù)表示為

(10)

在標(biāo)準(zhǔn)解體模型中,根據(jù)亞軌道碎片的特征尺寸劃分不同的區(qū)間范圍,亞軌道碎片在各個(gè)區(qū)間范圍內(nèi)的面質(zhì)比分布函數(shù)如下所示。

(1) 對于特征尺寸大于11 cm的亞軌道碎片,面質(zhì)比分布函數(shù)可由雙正態(tài)分布方程確定:

q(λc,η)=δ(η)q1(λc)+(1-δ(η))q2(λc)

(11)

式中:λc=lg(A/m);η=lgd;δ(η)表示權(quán)系數(shù),根據(jù)不同的解體類型確定;q1(λc)和q2(λc)表示正態(tài)分布概率密度函數(shù):

(12)

式中:μi和σi分別表示分布函數(shù)的均值和標(biāo)準(zhǔn)差。

亞軌道碎片面質(zhì)比分布函數(shù)的參數(shù)可通過下式計(jì)算:

(13)

σ1(η)=0.55

(14)

(15)

(16)

μ2(η)=-0.9

(17)

(2) 對于特征尺寸位于3～8 cm的亞軌道碎片,其面質(zhì)比分布函數(shù)可表示為

q(λc,η)=q1(λc)

(18)

即權(quán)系數(shù)δ(η)=1。分布函數(shù)的參數(shù)可通過下式計(jì)算:

(19)

(20)

(3) 對于特征尺寸位于8～11 cm的亞軌道碎片,其面質(zhì)比分布可通過插值的方法獲得。對于特征尺寸小于3 cm的碎片,考慮在下落過程中充分燃燒,本文中不做進(jìn)一步考慮。

在亞軌道飛行器解體時(shí)刻,解體碎片獲得相對于飛行器的速度增量Δv滿足正態(tài)分布,即

(21)

式中:v=lg Δv;均值μv=0.2λc+1.85;標(biāo)準(zhǔn)差σv通常取0.4。對于亞軌道飛行器來說,考慮在解體瞬間各個(gè)碎片獲得的速度增量的方向是隨機(jī)的,其方向余弦可表示為

(22)

式中:θ1和θ2為方向角,分別在[-π/2,π/2]和[0,2π]取值區(qū)間內(nèi)均勻分布。

3.2 阻力系數(shù)

亞軌道解體碎片再入過程中,努森數(shù)(Knudsen number, Kn)的大小對于計(jì)算亞軌道碎片在下落過程中所受到的空氣阻力十分重要[27],可表示為

(23)

式中:d表示亞軌道碎片特征尺寸;χ表示氣體分子平均自由程,數(shù)值大小與空氣分子的數(shù)量密度和分子大小相關(guān):

(24)

式中:l表示空氣的分子直徑,取值為l=3.65×10-10m。N表示1 m3體積內(nèi)空氣分子的總數(shù),取值大小跟海拔高度h相關(guān),計(jì)算公式為

N=(2.0×1025)·exp(-0.14h)+
(4.0×1017)·exp(-0.02h)

(25)

在空氣動力學(xué)中,根據(jù)Kn的大小,可劃分3個(gè)區(qū)間來計(jì)算亞軌道碎片的阻力系數(shù),依次可劃分為連續(xù)流區(qū)、過度流區(qū)和自由分子流區(qū)。對于解體碎片,其在大氣不同流區(qū)的阻力系數(shù)CD如下所示。

(1) 在連續(xù)流區(qū)(Kn≤0.01),亞軌道碎片的阻力系數(shù)為

CD-cont=0.91

(26)

(2) 在過渡流區(qū)(0.01

(27)

(3) 在自由分子流區(qū)(Kn≥1.0),亞軌道碎片的阻力系數(shù)為

CD-fm=2.1

(28)

4 仿真結(jié)果分析

(29)

通過上述航天器標(biāo)準(zhǔn)解體模型,首先根據(jù)不同碎片的特征尺寸劃分區(qū)間范圍,確定碎片的最小特征尺寸dc,利用式(9)計(jì)算碎片總數(shù)目。然后根據(jù)分布律式(10)、式(11)、式(18)、式(21)確定解體碎片的尺寸大小、面質(zhì)比、速度增量,通過采用反函數(shù)的方法對隨機(jī)變量進(jìn)行抽樣取均值計(jì)算,數(shù)據(jù)統(tǒng)計(jì)如表1所示。

表1 標(biāo)準(zhǔn)解體模型數(shù)據(jù)統(tǒng)計(jì)表Table 1 Statistical table of standard disintegration model data

在此次仿真模擬中使用廣泛接受的經(jīng)驗(yàn)密度模型質(zhì)譜儀非相干散射雷達(dá)擴(kuò)展2000[29](mass spectrometer incoherent scatter radar extended 2000, MSISE-00)和水平風(fēng)模型2014[30](horizontal wind model 2014, HWM14)。根據(jù)不同碎片特征尺寸和下落高度變化,計(jì)算碎片在不同高度區(qū)間內(nèi)的Kn,如表2所示,根據(jù)式(26)～式(28)確定阻力系數(shù)值。在此基礎(chǔ)上結(jié)合表1彈道系數(shù)計(jì)算值,根據(jù)碎片特征尺寸區(qū)間將亞軌道碎片數(shù)據(jù)信息分為5組,如表3所示。利用協(xié)方差傳播算法計(jì)算碎片云的傳播軌跡,通過設(shè)置不同置信度參數(shù)分別顯示3個(gè)不同時(shí)刻的碎片傳播軌跡,如圖2和圖3所示。

表2 Kn計(jì)算表Table 2 Kn calculation table

表3 碎片數(shù)據(jù)信息分析表Table 3 Debris data information analysis table

碎片危險(xiǎn)區(qū)邊界的劃定需要依據(jù)民航可接受的風(fēng)險(xiǎn)概率來確定。圖2和圖3分別顯示的是置信度為99.99%和95%時(shí)的橢球邊界范圍。置信度越高時(shí),橢球邊界范圍越大,越接近真實(shí)碎片下落傳播范圍。此外,不同彈道系數(shù)值對應(yīng)碎片群的下落時(shí)間也不相同,如圖4所示。彈道系數(shù)β越小,碎片下落時(shí)間越長,而β越大碎片下落的越快。

當(dāng)亞軌道解體事故發(fā)生時(shí),碎片危險(xiǎn)區(qū)會對飛機(jī)的安全飛行造成嚴(yán)重影響,在改航路線規(guī)劃之前首先要確定碎片危險(xiǎn)區(qū)域,碎片危險(xiǎn)區(qū)的預(yù)測精準(zhǔn)度將對后續(xù)改航路線的規(guī)劃產(chǎn)生重要影響。

5 結(jié) 論

本文針對亞軌道發(fā)射活動中飛行器潛在解體風(fēng)險(xiǎn),利用協(xié)方差傳播算法來預(yù)測亞軌道飛行器在高空解體后碎片的傳播軌跡,將位置概率橢球的概念用于仿真結(jié)果的可視化。在此基礎(chǔ)上,參照標(biāo)準(zhǔn)解體模型,對影響碎片傳播軌跡的關(guān)鍵因素彈道系數(shù)進(jìn)行推導(dǎo)計(jì)算,可有效提高碎片危險(xiǎn)區(qū)的預(yù)測準(zhǔn)確度。文章中充分考慮民航安全置信度的高要求,可通過設(shè)置置信度參數(shù)來適應(yīng)民航對生命安全的風(fēng)險(xiǎn)要求,為航班規(guī)避構(gòu)建安全、高效的決策支持系統(tǒng)。上述研究成果可以為推進(jìn)民航風(fēng)險(xiǎn)防控體系建設(shè)、促進(jìn)商業(yè)亞軌道飛行與民航協(xié)調(diào)一致發(fā)展提供技術(shù)參考。