程胡華，王益柏，蔡其發(fā)，趙 亮

(1. 63729部隊(duì)，太原030027；2. 61741部隊(duì)，北京 100094；3. 中國(guó)科學(xué)院大氣物理研究所大氣科學(xué)和地球流體力學(xué)數(shù)值模擬國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100029)

0 引 言

近年來(lái)，中國(guó)在航天領(lǐng)域取得一系列豐碩成果，運(yùn)載火箭作為航天技術(shù)領(lǐng)域里的核心之一，其內(nèi)部結(jié)構(gòu)非常復(fù)雜，并匯聚了無(wú)數(shù)尖端技術(shù)，為持續(xù)提高其完成任務(wù)的成功率，大量科學(xué)家對(duì)運(yùn)載火箭結(jié)構(gòu)特征、飛行控制等方面開(kāi)展了一系列研究[1-4]，并取得了豐富的成果；在發(fā)射之前，除確保運(yùn)載火箭本身滿(mǎn)足發(fā)射、飛行條件外，還需考慮大氣環(huán)境的影響，目前主要關(guān)注高空風(fēng)對(duì)運(yùn)載火箭飛行的影響[5-13]，例如，程鎮(zhèn)煌[5]通過(guò)全面敘述大型火箭風(fēng)載試驗(yàn)的各種問(wèn)題，指出火箭設(shè)計(jì)必須經(jīng)過(guò)靶場(chǎng)風(fēng)載試驗(yàn)和高空風(fēng)載試驗(yàn)，測(cè)出火箭各艙段載荷和儀器艙處風(fēng)激振幅等，準(zhǔn)確確定各種設(shè)備的適應(yīng)能力；余夢(mèng)倫[7]指出發(fā)現(xiàn)高空風(fēng)引起的氣流攻角對(duì)火箭飛行中的氣動(dòng)載荷和飛行環(huán)境有較大的影響，這種影響有時(shí)可造成火箭飛行失敗。

目前，在定量分析大氣環(huán)境對(duì)運(yùn)載火箭飛行影響(即qαmax)的計(jì)算程序中，除了需要大氣高空風(fēng)參數(shù)外，還需要大氣密度參數(shù)，但程序中的大氣密度參數(shù)通常采用參考大氣密度，而非發(fā)射場(chǎng)零日的真實(shí)大氣密度，針對(duì)大氣密度對(duì)運(yùn)載火箭或其它飛行器飛行的影響，廣大科學(xué)工作者已開(kāi)展了一系列工作，并取得大量有意義成果[14-23]，黃世勇和聞悅[14]指出大氣密度對(duì)飛行器的安全及返回有非常重要影響，李健等[20]研究發(fā)現(xiàn)大氣密度是建立飛行動(dòng)力學(xué)模型的重要參數(shù)，不同的大氣密度模型對(duì)彈道終點(diǎn)高度、最大過(guò)載、最大熱流等參數(shù)影響明顯；程旋等[21]指出飛行器的氣動(dòng)力和氣動(dòng)力矩都是大氣密度的函數(shù)，大氣密度的變化直接影響到飛行狀態(tài)下飛行器的受力，對(duì)飛行器的姿態(tài)影響顯著，引起攻角、側(cè)滑角、偏航角等角度的變化，導(dǎo)致再入飛行器偏離預(yù)定軌跡；因此，大量科技工作者針對(duì)大氣密度開(kāi)展一系列研究工作[24-27]。

獲取參考大氣密度的參考大氣模型有CIRA系列、MSIS系列等，由于大氣密度受來(lái)自太陽(yáng)輻射及大氣波動(dòng)的共同影響，變化十分復(fù)雜，人們至今尚未完全掌握其變化的物理機(jī)制，導(dǎo)致已有的參考大氣模型精度不高[19,24,28]，因此，在采用參考大氣密度定量計(jì)算大氣環(huán)境對(duì)運(yùn)載火箭飛行的qαmax時(shí)，由于參考大氣密度與真實(shí)大氣密度之間存在偏差，但針對(duì)其偏差對(duì)運(yùn)載火箭飛行的qαmax精度特征研究尚未見(jiàn)相關(guān)文獻(xiàn)，考慮到目前計(jì)算程序仍采用參考大氣密度，本文利用2017年8月1日-2019年7月31日晉西北地區(qū)探空資料，分別利用MSIS00模式參考大氣密度、真實(shí)大氣密度計(jì)算某型運(yùn)載火箭飛行時(shí)的qαmax，一方面對(duì)MSIS00模式參考大氣密度與真實(shí)大氣密度之間的差異特征進(jìn)行研究，另一方面對(duì)它們之間的qαmax差異進(jìn)行統(tǒng)計(jì)，最后，利用多元線(xiàn)性回歸、BP神經(jīng)網(wǎng)絡(luò)方法對(duì)基于MSIS00模式參考大氣密度得到的qαmax建立訂正模型，并對(duì)訂正效果進(jìn)行檢驗(yàn)。

1 資料及方法

1.1 探空資料

探空資料為晉西北地區(qū)每日08時(shí)(北京時(shí)，下同)，時(shí)間長(zhǎng)度為2017年8月1日-2019年7月31日，所需氣象要素為高度、氣壓、風(fēng)速、風(fēng)向、溫度和相對(duì)濕度，單位分別為m，hPa，m·s-1，(°)，℃和%，考慮到qαmax值常出現(xiàn)在16 km高度以下，故選取探空資料的幾何高度范圍為0～20000 m，其中，0～10000 m高度范圍內(nèi)間隔為200 m，10000～20000 m高度范圍內(nèi)間隔為500 m。

1.2 計(jì)算真實(shí)大氣密度

利用高垂直分辨率探空資料計(jì)算真實(shí)大氣密度，對(duì)于濕空氣，利用溫度、相對(duì)濕度、氣壓可計(jì)算大氣密度，計(jì)算步驟如下[29]：

1)計(jì)算飽和水汽壓

(1)

式中：T為溫度。

2)計(jì)算水汽壓

e=10-2RH·es

(2)

式中：RH為相對(duì)濕度。

3)計(jì)算比濕

(3)

式中：p為氣壓。

4)計(jì)算濕空氣的比氣體常數(shù)

Rv=Rd(1+0.608q)

(4)

式中：Rd為干空氣的比氣體常數(shù)，數(shù)值為287.05，單位為J/(kg·K)。

5)利用濕空氣的狀態(tài)方程，可計(jì)算濕空氣大氣密度

(5)

1.3 MSIS參考大氣模式

MSIS(Mass Spectrometer and Incoherent Scatter empirical model)大氣模式最初版本是由美國(guó)國(guó)家航空航天局空間飛行中心編制的經(jīng)驗(yàn)大氣預(yù)報(bào)模式，通過(guò)利用多顆衛(wèi)星質(zhì)譜計(jì)所測(cè)得的成分資料和非相干散射雷達(dá)探測(cè)資料而建立，目前MSIS00是該系列模式的較新版本[30]，其利用了更新、更多的雷達(dá)、衛(wèi)星和火箭資料，是當(dāng)今國(guó)際上廣泛使用的參考大氣模式[31-35]，MSIS00模式的適用高度范圍為地面至高層大氣，可給出不同高度層的大氣溫度、大氣密度及成分結(jié)構(gòu)。

1.4 統(tǒng)計(jì)方法

以真實(shí)大氣密度為基準(zhǔn)，通過(guò)計(jì)算參考大氣密度與真實(shí)大氣密度之間的偏差、相對(duì)偏差和偏差區(qū)間占有率，反映參考大氣密度誤差特征，其相應(yīng)的計(jì)算公式如下：

(6)

(7)

式中：Bias為偏差，Rerr為相對(duì)偏差，n為樣本數(shù)，x和y分別代表參考大氣密度與真實(shí)大氣密度。

文中的偏差區(qū)間占有率計(jì)算式定義如下：

(8)

式中：若xi-yi<-N結(jié)果成立，則定義xi-yi<-N的值為1，否則值為0;該定義同樣適用于表達(dá)式-N≤xi-yi≤N和xi-yi>N。

2 密度差異特征

2.1 平均差異特征

在2017年8月1日-2019年7月31日，參考大氣密度與真實(shí)大氣密度之間的偏差、相對(duì)偏差隨高度變化特征如圖1所示，偏差隨高度增加呈現(xiàn)正值、負(fù)值的交替變化特征(見(jiàn)圖1(a))，其正偏差值主要位于高度2.00～7.60 km，負(fù)偏差值主要位于高度11.40～19.90 km，整個(gè)高度范圍內(nèi)的值范圍為-0.0083 kg·m-3(14.90 km)～0.0065 kg·m-3(4.20 km)，平均值為5.2843×10-4kg·m-3；類(lèi)似偏差(見(jiàn)圖1(a))，相對(duì)偏差隨高度增加同樣呈現(xiàn)正值、負(fù)值的交替變化特征(見(jiàn)圖1(b))，但高層的相對(duì)偏差值明顯大于低層，其正相對(duì)偏差值主要位于高度2.00～7.60 km，負(fù)相對(duì)偏差主要位于高度11.40～19.90 km，整個(gè)高度范圍內(nèi)的值范圍為-4.2142%(15.90 km)～5.2210%(21.40 km)，平均值為-0.3063%；對(duì)照?qǐng)D1(a)與圖1(b)可知，雖然偏差、相對(duì)偏差隨高度變化趨勢(shì)基本一致，但大值區(qū)所處的高度層存在一定差異，例如，偏差極大值(0.0065 kg·m-3)、極小值(-0.0083 kg·m-3)分別位于4.20 km，14.90 km，而相對(duì)偏差極大值(5.2210%)、極小值(-4.2142%)分別位于21.40 km和15.90 km。

圖1 在2017年8月1日-2019年7月31日，參考大氣密度與真實(shí)大氣密度之間的密度偏差、密度相對(duì)偏差隨高度變化特征Fig.1 The altitude variations of deviations and relative deviations between the reference atmospheric density and the true atmospheric density from August 1, 2017 to July 31, 2019

2.2 年內(nèi)差異特征

在2017年8月1日-2019年7月31日，參考大氣密度與真實(shí)大氣密度之間偏差、相對(duì)偏差的年內(nèi)—高度變化特征如圖2所示；對(duì)于年內(nèi)變化特征，在低層，偏差值呈現(xiàn)負(fù)值、正值、負(fù)值變化特征，其正偏差極大值出現(xiàn)在7月，而在高層則表現(xiàn)為正值、負(fù)值、正值變化特征，其負(fù)偏差極大值出現(xiàn)在8月(見(jiàn)圖2(a))；對(duì)于高度變化特征，1-3月、11-12月偏差隨高度增加呈現(xiàn)出由負(fù)偏差向正偏差轉(zhuǎn)變過(guò)程，而在4-10月則表現(xiàn)為正偏差向負(fù)偏差轉(zhuǎn)變過(guò)程；在圖2(a)中，低層(2.0～4.0 km)的6-8月存在一個(gè)正偏差大值區(qū)，高層(13～17 km)的7-9月存在一個(gè)負(fù)偏差大值區(qū)，最小值為-0.0538 kg·m-3(12月的2.00 km)，最大值為0.0558 kg·m-3(7月的2.0 km)，平均值為5.3219×10-4kg·m-3；類(lèi)似偏差(見(jiàn)圖2(a))，相對(duì)偏差同樣存在明顯的年內(nèi)變化及隨高度變化特征(見(jiàn)圖2(b))，對(duì)于年內(nèi)變化特征，低層(10 km以下)呈現(xiàn)出負(fù)值、正值、負(fù)值的變化特征，而高層(10 km以上)則表現(xiàn)為正值、負(fù)值、正值的變化特征；對(duì)于高度變化特征，1-3月、11-12月隨高度增加呈現(xiàn)出負(fù)相對(duì)偏差向正相對(duì)偏差轉(zhuǎn)變特征，而4-10月隨高度增加呈現(xiàn)出正相對(duì)偏差向負(fù)相對(duì)偏差轉(zhuǎn)變，圖2(b)中的負(fù)相對(duì)偏差最大值為-13.2699%(8月的15.90 km)，正相對(duì)偏差最大值為10.1708%(12月的21.40 km)，平均值為-0.3046%。

圖2 2017年8月1日-2019月7月31日，參考大氣密度與真實(shí)大氣密度之間的密度偏差、密度相對(duì)偏差隨年內(nèi)—高度變化特征Fig.2 The time-altitude cross section of deviations and relative deviations between the reference atmospheric density and the true atmospheric density from August 1, 2017 to July 31, 2019

2.3 季節(jié)差異特征

在2017年8月1日-2019年7月31日，不同季節(jié)參考大氣密度與真實(shí)大氣密度之間的偏差、相對(duì)偏差隨高度變化特征如圖3所示，偏差、相對(duì)偏差隨高度變化特征與季節(jié)有關(guān)，春、夏、秋季的大氣密度偏差在低層以正值為主，偏差值由大到小分別為夏季、春季、秋季，在高層以負(fù)值為主，偏差值由大到小分別為夏季、秋季、春季；而冬季的大氣密度偏差在低層以負(fù)值為主，在高層以正值為主(見(jiàn)圖3(a))，在整個(gè)高度范圍內(nèi)，春、夏、秋、冬季的大氣密度偏差范圍分別為-0.0055 kg·m-3(11.90 km)～0.0120 kg·m-3(2.00 km)，-0.0234 kg·m-3(14.90 km)～0.0515 kg·m-3(2.00 km)，-0.0097 kg·m-3(14.90 km)～0.0047 kg·m-3(4.20 km)，-0.0505 kg·m-3(2.00 km)～0.0134 kg·m-3(11.40 km)，平均值分別為4.8702×10-4kg·m-3，0.0083 kg·m-3，-0.0013 kg·m-3，-0.0054 kg·m-3；類(lèi)似圖3(a)，相對(duì)偏差隨高度變化特征同樣與季節(jié)密切相關(guān)(見(jiàn)圖3(b))，其中，春、夏、秋季的大氣密度相對(duì)偏差在低層以正值為主，相對(duì)偏差值由大到小分別為夏季、春季、秋季，在高層以負(fù)值為主，相對(duì)偏差值由大到小分別為夏季、秋季、春季；而冬季的大氣密度相對(duì)偏差在低層以負(fù)值為主、高層以正值為主，在整個(gè)高度范圍內(nèi)，春、夏、秋、冬季的大氣密度相對(duì)偏差范圍分別為-3.0923%(17.90 km)～5.9970%(21.40 km)，-11.6516%(15.90 km)～5.4388%(2.00 km)，-4.7972%(15.90 km)～5.2188%(21.40 km)，-4.7440%(2.00 km)～10.0290%(21.40 km),平均值分別為-0.2805%，-0.7126%，-0.7011%，0.4637%，對(duì)照?qǐng)D3(a)和圖3(b)及上述分析可知，雖然大氣密度偏差與相對(duì)偏差隨高度變化趨勢(shì)基本一致，但出現(xiàn)極大值、極小值的高度層存在差異，且在整個(gè)高度層，春季和夏季的平均偏差為正值、秋季和冬季的平均偏差為負(fù)值；而春、夏、秋季的相對(duì)偏差均為負(fù)值，冬季為正值。

圖3 2017年8月1日-2019年7月31日，不同季節(jié)參考大氣密度與真實(shí)大氣密度之間的密度偏差、密度相對(duì)偏差隨高度變化特征Fig.3 The altitude variations of deviations and relative deviations between the reference atmospheric density and the true atmospheric density from August 1, 2017 to July 31, 2019 in different seasons

2.4 平均年內(nèi)變化特征

在2017年8月1日-2019年7月31日的1-12月，整個(gè)高度層內(nèi)的平均大氣密度偏差呈現(xiàn)“∧”字符變化特征(見(jiàn)圖4(a))，其中，1-3月、10-12月為負(fù)偏差，4-9月為正偏差，圖4(a)中的值范圍為-0.0060 kg·m-3(11月)～0.0099 kg·m-3(7月)，平均值為5.3219×10-4kg·m-3；整個(gè)高度層內(nèi)的平均大氣密度相對(duì)偏差年內(nèi)變化特征如圖4(b)所示，與圖4(a)基本相反(從圖2(a)與圖2(b)之間的高層、低層差異可以解釋)，圖4(b)中的1-3月和12月為正相對(duì)偏差、4-11月為負(fù)相對(duì)偏差，其值范圍為-0.9569%(8月)～0.7225%(2月)，平均值為-0.3046%。

圖4 2017年8月1日-2019年7月31日，整個(gè)高度層內(nèi)參考大氣密度與真實(shí)大氣密度之間的密度偏差、密度相對(duì)偏差年內(nèi)變化特征Fig.4 The annual variation of deviations and relative deviations between the reference atmospheric density and the true atmospheric density within the entire height from August 1, 2017 to July 31, 2019

3 qαmax差異特征

3.1 qαmax差異統(tǒng)計(jì)

通過(guò)上述研究發(fā)現(xiàn)，參考大氣密度與真實(shí)大氣密度之間的差異特征與高度、季節(jié)均有密切關(guān)系，在2017年8月1日-2019年7月31日，參考大氣密度與真實(shí)大氣密度之間的qαmax偏差特征如圖5所示，qαmax偏差呈現(xiàn)“M”字型變化特征(見(jiàn)圖5(a)，圖中粗虛線(xiàn)分別對(duì)應(yīng)+100 Pa·rad，-100 Pa·rad，圖5(a)中的值范圍為-217.1721～219.0648 Pa·rad，平均值為-9.4684 Pa·rad，其中，正偏差超過(guò)100 Pa·rad的樣本數(shù)有24個(gè)，負(fù)偏差超過(guò)-100 Pa·rad的樣本數(shù)有49個(gè)；qαmax偏差在不同區(qū)間范圍內(nèi)的占有率特征如圖5(b)所示，沿橫坐標(biāo)軸方向呈現(xiàn)增大、減小的變化趨勢(shì)，其中，qαmax偏差在(-20, 0]范圍內(nèi)占有率最大(20.4110%)，在[0，20)范圍內(nèi)占有率為次大值(16.0274%)，在(-∞， -100]、(100， +∞)范圍內(nèi)的占有率分別為6.7123%，3.2877%；類(lèi)似平均大氣密度偏差存在明顯的年內(nèi)變化特征(見(jiàn)圖4(a))，qαmax偏差同樣具有明顯的年內(nèi)變化特征(見(jiàn)圖5(c))，但變化特征與圖4(a)相反，圖5(c)中的qαmax偏差在1-3月、11-12月為正偏差，4-10月為負(fù)偏差，在1-12月的值范圍為-51.3659 Pa·rad(6月)～46.3820 Pa·rad(1月)，平均值為-9.5690 Pa·rad，qαmax偏差同樣與季節(jié)有關(guān)(見(jiàn)圖5(d))，其中，春、夏、秋季均為負(fù)偏差，冬季為正偏差，負(fù)偏差最大值出現(xiàn)在夏季(-42.5615 Pa·rad)，正偏差最大值出現(xiàn)在冬季(37.8531 Pa·rad)。

3.2 qαmax偏差超過(guò)100 Pa·rad

參考大氣密度與真實(shí)大氣密度之間的qαmax偏差超過(guò)100 Pa·rad時(shí)，其分布特征如圖6所示，其中，有24個(gè)樣本的qαmax偏差超過(guò)100 Pa·rad(見(jiàn)圖6(a))，有49個(gè)樣本的qαmax偏差超過(guò)-100 Pa·rad(見(jiàn)圖6(b))，在圖6(a)中，qαmax偏差值范圍為101.1910～219.0648 Pa·rad，其中，基于參考大氣密度得到的qαmax值范圍為1.8522×103～4.7914×103Pa·rad，而基于真實(shí)大氣密度得到的qαmax值范圍為1.7337×103～4.6090×103Pa·rad；在圖6(b)中，qαmax偏差值范圍為-101.4770～-217.1721Pa·rad，其中，基于參考大氣密度得到的qαmax值范圍為873.2646～2.8181×103Pa·rad，而基于真實(shí)大氣密度得到的qαmax值范圍為982.6712～3.0027×103Pa·rad。

在2017年8月1日-2019年7月31日，qαmax偏差值超過(guò)100 Pa·rad時(shí)，對(duì)應(yīng)的風(fēng)場(chǎng)、大氣密度偏差隨高度變化特征如圖7所示。從圖7(a)可以看出，雖然qαmax偏差≤-100、qαmax偏差≥100對(duì)應(yīng)的風(fēng)場(chǎng)隨高度變化趨勢(shì)基本一致，且極大值出現(xiàn)在基本相同的高度層(12 km左右)、極小值均出現(xiàn)在2 km，但qαmax偏差≥100對(duì)應(yīng)的風(fēng)速在各高度層均大于相應(yīng)的qαmax偏差≤-100，在整個(gè)高度范圍內(nèi)，qαmax偏差≤-100，qαmax偏差≥100對(duì)應(yīng)的風(fēng)速范圍分別為5.9592 m·s-1(2.8 km)～40.0294 m·s-1(12.90 km)，7.2083 m·s-1(2.0 km)～58.8748 m·s-1(11.90 km)，平均值分別為18.0332 m·s-1，33.8605 m·s-1，qαmax偏差≤-100、qαmax偏差≥100對(duì)應(yīng)的大氣密度偏差隨高度變化特征存在明顯差異(見(jiàn)圖7(b))，其中，qαmax偏差≤-100對(duì)應(yīng)大氣密度偏差隨高度增加呈現(xiàn)正偏差向負(fù)偏差轉(zhuǎn)變，而相應(yīng)的qαmax偏差≥100為負(fù)偏差向正偏差轉(zhuǎn)變，即兩者表現(xiàn)為相反的變化特征，在整個(gè)高度范圍內(nèi)，qαmax偏差≤-100，qαmax偏差≥100對(duì)應(yīng)的大氣密度偏差范圍分別為-0.0217 kg·m-3(13.90 km)～0.0386 kg·m-3(2.0 km)，-0.0539 kg·m-3(2.0 km)～0.0178 kg·m-3(11.0 km)，平均值分別為0.0058 kg·m-3，-0.0052 kg·m-3。

圖5 2017年8月1日-2019年7月31日，qαmax偏差變化及其在不同區(qū)間范圍內(nèi)的占有率、年內(nèi)變化和季節(jié)變化特征Fig.5 The change characteristics, percentages in different ranges, annual variation and seasonal variation of qαmax deviation from August 1, 2017 to July 31, 2019

圖6 在qαmax偏差≥100 Pa·rad、qαmax偏差≤-100 Pa·rad時(shí)，基于真實(shí)大氣密度與參考大氣密度計(jì)算得到的qαmax值變化特征Fig.6 Under the qαmaxdeviation≥100 Pa·rad, qαmax deviation≤-100 Pa·rad conditions, the variation characteristics of qαmax value calculated based on true atmospheric density and reference atmospheric density

圖7 qαmax偏差值超過(guò)100 Pa·rad時(shí)，風(fēng)速、密度偏差隨高度變化特征Fig.7 The altitude variations of upper wind speed, atmospheric density deviation when the value of qαmaxdeviation exceeds 100 Pa·rad

4 建 模

基于目前計(jì)算大氣環(huán)境對(duì)運(yùn)載火箭飛行的qαmax程序包中，大氣密度為參考大氣密度、而非真實(shí)大氣密度，通過(guò)上述研究可以看出，由于參考大氣密度與真實(shí)大氣密度之間存在差異，導(dǎo)致該程序包計(jì)算得到的qαmax與真實(shí)值之間存在差異，偏差最大值超過(guò)200 Pa·rad，為減小運(yùn)載火箭飛行的風(fēng)險(xiǎn)，一方面，強(qiáng)烈建議采用真實(shí)大氣密度代替程度包中的參考大氣密度，另一方面，可考慮采用多元線(xiàn)性回歸、BP神經(jīng)網(wǎng)絡(luò)等建模方法建立訂正模型，本文采用常用的多元線(xiàn)性回歸方法、BP神經(jīng)網(wǎng)絡(luò)方法建立訂正模型，并對(duì)訂正模型的效果進(jìn)行檢驗(yàn)。

4.1 建模思想

基于2017年8月1日-2019年7月31日，一共730個(gè)樣本，其中:1)前600個(gè)樣本作為訓(xùn)練樣本，后130個(gè)樣本作為檢驗(yàn)樣本；2)基于參考大氣密度計(jì)算得到qαmax對(duì)應(yīng)的運(yùn)載火箭飛行高度h，獲取0～h高度范圍內(nèi)的緯向風(fēng)、經(jīng)向風(fēng)、參考大氣密度平均值，分別記為umean，vmean，ρmean；3)日期以天數(shù)表示，其中，1月1日為1，1月31日為31，…，記為d；4)計(jì)算參考大氣密度與真實(shí)大氣密度之間的qαmax偏差；5)對(duì)2017年8月1日-2019年7月31日(共730個(gè)樣本)的h，umean，vmean，ρmean，d歸一化處理，即將數(shù)值轉(zhuǎn)換為[-1,1]范圍內(nèi)；6)針對(duì)前600個(gè)樣本，利用多元線(xiàn)性回歸、BP神經(jīng)網(wǎng)絡(luò)方法對(duì)h，umean，vmean，ρmean，d與qαmax偏差之間建立模型；7)基于已建立好的訂正模型，利用后130個(gè)樣本的h，umean，vmean，ρmean，d參數(shù)，得到訂正模型的qαmax偏差值(記為Dqαmax)，Dqαmax與參考大氣密度計(jì)算得到的qαmax(記為Rqαmax)之和，即為所需要的結(jié)果。

基于參考大氣密度、真實(shí)大氣密度得到的qαmax分別記為Rqαmax，Tqαmax，BP神經(jīng)網(wǎng)絡(luò)、多元線(xiàn)性回歸方法建立訂正模型得到的qαmax，分別記為Bqαmax，Mqαmax；本文基于多元線(xiàn)性回歸方法建立的模型如下：

49.89umean-40.42vmean-83.18ρmean-9.94d-

103.56h-23.08=Dqαmax

(9)

Mqαmax=Rqαmax-Dqαmax

(10)

在式(9)～(10)中，umean，vmean，ρmean，d和h的變量意義與前面描述一致，為后130個(gè)樣本的已知量，Rqαmax為基于參考大氣密度得到的后130個(gè)樣本qαmax；利用式(9)可計(jì)算得到Dqαmax，Dqαmax反映了Rqαmax與基于多元線(xiàn)性回歸方法得到Mqαmax之間的偏差；利用式(10)得到的Mqαmax，即為基于多元線(xiàn)性回歸方法建立模型得到的qαmax訂正值。

本文使用的BP神經(jīng)網(wǎng)絡(luò)方法，是采用MATLAB軟件自帶的神經(jīng)網(wǎng)絡(luò)工具箱進(jìn)行建模，同樣先得到Dqαmax，其中，Dqαmax反映了Rqαmax與基于BP神經(jīng)網(wǎng)絡(luò)方法得到Bqαmax之間的偏差；利用式(11)即可得到基于BP神經(jīng)網(wǎng)絡(luò)方法的Bqαmax，即

Bqαmax=Rqαmax-Dqαmax

(11)

式中：Rqαmax為基于參考大氣密度得到的后130個(gè)樣本qαmax。

記Rqαmax，Bqαmax，Mqαmax與Tqαmax之間的偏差分別為RqαmaxDev，BqαmaxDev，MqαmaxDev，下面利用數(shù)理統(tǒng)計(jì)方法對(duì)BP神經(jīng)網(wǎng)絡(luò)、多元線(xiàn)性回歸方法建立的訂正模型效果進(jìn)行檢驗(yàn)。

4.2 模型檢驗(yàn)

在后130個(gè)樣本中，RqαmaxDev，BqαmaxDev，MqαmaxDev的變化特征及在不同區(qū)間范圍占有率分布特征如圖8所示，從圖8(a)可以看出，通過(guò)對(duì)參考大氣密度的結(jié)果進(jìn)行BP神經(jīng)網(wǎng)絡(luò)、多元線(xiàn)性回歸建模訂正處理，所得到的qαmax偏差明顯減小，在圖8(a)中，RqαmaxDev，BqαmaxDev，MqαmaxDev的值范圍分別為-201.6939～122.1358 Pa·rad，-75.9311～39.4346 Pa·rad，-52.6029～60.1846 Pa·rad；RqαmaxDev，BqαmaxDev，MqαmaxDev在不同區(qū)間范圍內(nèi)的占有率存在明顯差異(見(jiàn)圖8(b))，其中，RqαmaxDev在區(qū)間(-20，0]占有率最高(19.2308%)，在區(qū)間(0，20]占有率次高(15.3846%)，BqαmaxDev在區(qū)間(-20，0]占有率最高(31.5385%)，在區(qū)間(-40，-20]占有率次高(27.6923%)，MqαmaxDev在區(qū)間(-40，-20]占有率最高(27.6923%)，在區(qū)間(0，20]占有率次高(26.9231%)，從圖8(b)還可以看出，RqαmaxDev，BqαmaxDev，MqαmaxDev占有率為0.00%的區(qū)間范圍也存在明顯差異，RqαmaxDev在區(qū)間(60,80]占有率為0.00%，BqαmaxDev在區(qū)間(-∞，-100]，(-100，-80]，(40，60]，(60，80]，(80，100]，(100，∞]內(nèi)的占有率均為0.00%，MqαmaxDev在區(qū)間(-∞，-100]，(-100，-80]，(-80，-60]，(80，100]，(100，∞]內(nèi)的占有率均為0.00%。Rqαmax，Bqαmax，Mqαmax與Tqαmax之間的平均偏差、平均絕對(duì)差、平均相對(duì)偏差、平均標(biāo)準(zhǔn)差及相關(guān)系數(shù)統(tǒng)計(jì)結(jié)果見(jiàn)表1。

圖8 RqαmaxDev，BqαmaxDev，MqαmaxDev的變化及其在不同區(qū)間范圍占有率分布特征Fig.8 The variation characteristics, and percentages in different ranges of RqαmaxDev，BqαmaxDev，MqαmaxDev

表1 Rqαmax，Bqαmax，Mqαmax與Tqαmax之間的數(shù)理統(tǒng)計(jì)特征Table 1 The mathematical statistics results between Rqαmax，Bqαmax，Mqαmax and Tqαmax

5 結(jié) 論

運(yùn)載火箭飛行前，需確定發(fā)射場(chǎng)大氣環(huán)境對(duì)運(yùn)載火箭飛行的影響，若計(jì)算得到的qαmax超過(guò)閾值，則發(fā)射推遲；在計(jì)算qαmax的程序中，所用的大氣密度通常是參考大氣密度，而非發(fā)射場(chǎng)零日的真實(shí)大氣密度，目前針對(duì)參考大氣密度對(duì)qαmax精度的影響研究較少，本文以2017年8月1日-2019年7月31日高垂直分辨率探空資料為例，分析參考大氣密度與真實(shí)大氣密度以及它們計(jì)算得到某型運(yùn)載火箭飛行qαmax之間的差異特征，最后建立訂正模型，得到主要結(jié)論如下：

1)密度偏差、相對(duì)偏差隨高度增加而呈現(xiàn)正值、負(fù)值的交替變化特征，整個(gè)高度層內(nèi)的平均大氣密度偏差、相對(duì)偏差具有明顯年內(nèi)變化特征，其中，偏差值在1-3月、10-12月為負(fù)值，4-9月為正值，值范圍為-0.0060 kg·m-3(11月)～0.0099 kg·m-3(7月)，相對(duì)偏差在4-11月為負(fù)值，1-3月、12月為正值，值范圍為-0.9569%(8月)～0.7225%(2月)。

2)參考大氣密度與真實(shí)大氣密度之間的qαmax偏差隨時(shí)間變化呈現(xiàn)“M”字型特征，其值范圍為-217.1721～219.0648 Pa·rad，平均值為-9.4684 Pa·rad；qαmax偏差值在不同區(qū)間范圍內(nèi)占有率存在較明顯差異，其最大占有率(20.4110%)位于(-20,0]區(qū)間范圍。

3)參考大氣密度與真實(shí)大氣密度之間的qαmax偏差，具有明顯的年內(nèi)、季節(jié)變化特征；其中，1-3月、11-12月為正偏差，4-10月為負(fù)偏差，值范圍為-51.3659 Pa·rad(6月)～46.3820 Pa·rad(1月)，春、夏、秋季為負(fù)偏差，冬季為正偏差，值范圍為-42.5615 Pa·rad(夏季)～37.8531 Pa·rad(冬季)。

4)通過(guò)對(duì)qαmax偏差≥100 Pa·rad、qαmax偏差≤-100 Pa·rad對(duì)應(yīng)的風(fēng)場(chǎng)、大氣密度偏差進(jìn)行對(duì)比發(fā)現(xiàn)，qαmax偏差≥100 Pa·rad對(duì)應(yīng)的風(fēng)速在各高度層均明顯偏大、大氣密度偏差隨高度增加呈現(xiàn)負(fù)值向正值轉(zhuǎn)變，而qαmax偏差≤-100 Pa·rad對(duì)應(yīng)的風(fēng)速在各高度層均偏小、且大氣密度偏差為正值向負(fù)值轉(zhuǎn)變。

5)基于參考大氣密度得到的qαmax偏差，通過(guò)利用多元線(xiàn)性回歸、BP神經(jīng)網(wǎng)絡(luò)方法建立模型進(jìn)行訂正，結(jié)果表明均能夠較好地提高qαmax精度，且多元線(xiàn)性回歸建立的模型改進(jìn)效果更好；與參考大氣密度得到的qαmax相比，多元線(xiàn)性回歸、BP神經(jīng)網(wǎng)絡(luò)方法建立模型得到的qαmax偏差由-36.1569 Pa·rad減小到-7.7012 Pa·rad，-17.2267 Pa·rad，絕對(duì)差由49.3254 Pa·rad減小到20.3244 Pa·rad，23.9357 Pa·rad，相對(duì)偏差由-2.0201%減小到-0.7741%，-1.2917%，標(biāo)準(zhǔn)差由55.5491 Pa·rad減小到23.2982 Pa·rad，23.8563 Pa·rad，相關(guān)系數(shù)由0.9952提高到0.9960，0.9960。

通過(guò)上述研究結(jié)果可以看出，由于參考大氣密度與真實(shí)大氣密度之間的偏差存在年內(nèi)、季節(jié)變化特征，導(dǎo)致所計(jì)算得到的qαmax偏差同樣存在較明顯年內(nèi)、季節(jié)變化，且正、負(fù)偏差最大值均超過(guò)200 Pa·rad；考慮到目前計(jì)算程序中所用的大氣密度通常是參考大氣密度，而非發(fā)射場(chǎng)零日真實(shí)大氣密度，若程序計(jì)算得到的qαmax接近閾值，則此時(shí)真實(shí)的qαmax可能超過(guò)閾值，導(dǎo)致運(yùn)載火箭飛行存在重大風(fēng)險(xiǎn)。雖然采用多元線(xiàn)性回歸或BP神經(jīng)網(wǎng)絡(luò)建立訂正模型可以提高qαmax精度，但同樣存在一定風(fēng)險(xiǎn)，因此，在分析大氣環(huán)境對(duì)運(yùn)載火箭飛行的影響時(shí)，最好采用發(fā)射場(chǎng)零日真實(shí)大氣密度代替參考大氣密度，從而得到運(yùn)載火箭飛行時(shí)真實(shí)的qαmax，進(jìn)一步提高運(yùn)載火箭飛行的保障能力，避免出現(xiàn)重大事故；考慮到火箭發(fā)射前，需提前幾天確定發(fā)射場(chǎng)區(qū)發(fā)射零日是否滿(mǎn)足發(fā)射條件，然而，此時(shí)沒(méi)有未來(lái)幾天的大氣密度實(shí)況資料、且考慮大氣密度與多種因素有關(guān)(溫度、相對(duì)濕度等)，因此，通過(guò)建立大氣密度預(yù)報(bào)模型，獲取更接近真實(shí)大氣密度的大氣密度預(yù)報(bào)值，代替計(jì)算程序中的參考大氣密度，可能是后期進(jìn)一步提高精度的研究方向。