曹樹新,李 巖,李帥孝,徐樹威

(1.南京理工大學(xué) 能源與動(dòng)力工程學(xué)院,江蘇 南京 210094;2.陸軍裝備部 沈陽地區(qū)軍事代表局駐沈陽地區(qū)第二軍事代表室,遼寧 沈陽 110004;3.齊齊哈爾建華機(jī)械有限公司,黑龍江 齊齊哈爾 161006)

隨著快速機(jī)動(dòng)能力的不斷增強(qiáng)和偵察校射技術(shù)的不斷發(fā)展,現(xiàn)代戰(zhàn)爭(zhēng)對(duì)火炮的快速反應(yīng)能力提出了更高要求,氣象保障作為火炮射擊保障的核心環(huán)節(jié)之一,有必要在提升氣象保障的速度和精度方面加以研究[1]。近年來,數(shù)值天氣預(yù)報(bào)在火炮氣象保障中的應(yīng)用極大縮短了氣象保障的準(zhǔn)備時(shí)間[2-5],但是,初始場(chǎng)測(cè)量處理誤差和數(shù)值模式建模誤差,使得數(shù)值預(yù)報(bào)結(jié)果中含有較大的預(yù)報(bào)誤差,且尤以風(fēng)的預(yù)報(bào)誤差最為顯著[6]?，F(xiàn)行有效的數(shù)值修正方法以數(shù)值天氣預(yù)報(bào)的預(yù)報(bào)誤差及其變化規(guī)律為依據(jù),對(duì)臨近天氣預(yù)報(bào)的預(yù)報(bào)誤差進(jìn)行數(shù)值修正[7],提升氣象要素的預(yù)報(bào)精度。因此,如何快速準(zhǔn)確獲取預(yù)報(bào)誤差已成為提升彈道解算用數(shù)值天氣預(yù)報(bào)的預(yù)報(bào)精度需首要解決的問題之一。

目前,利用成熟的火炮外彈道跟蹤觀測(cè)技術(shù)獲取的實(shí)測(cè)彈道數(shù)據(jù),為快速準(zhǔn)確獲取彈道解算用數(shù)值天氣預(yù)報(bào)的預(yù)報(bào)誤差提供了新途徑。本文為獲取風(fēng)的預(yù)報(bào)誤差,使用風(fēng)修正系數(shù)描述風(fēng)的預(yù)報(bào)誤差,提出一種采用無跡卡爾曼濾波算法,直接從實(shí)測(cè)彈道數(shù)據(jù)和原始預(yù)報(bào)氣象數(shù)據(jù)中辨識(shí)風(fēng)修正系數(shù)的方法,并通過數(shù)值比對(duì)和彈道重構(gòu),對(duì)風(fēng)修正系數(shù)的辨識(shí)精度進(jìn)行分析。

1 系統(tǒng)模型

1.1 WRF模式簡(jiǎn)介及參數(shù)設(shè)置

WRF(weather research and forecast,WRF)模式是由美國多家科研機(jī)構(gòu)于2000年開發(fā)的新一代中尺度數(shù)值模式和數(shù)據(jù)同化系統(tǒng),其廣泛應(yīng)用于中小尺度到全球尺度的數(shù)值預(yù)報(bào)和大氣模擬,并可提供指定區(qū)域在未來時(shí)段內(nèi)的時(shí)空網(wǎng)格化氣象數(shù)據(jù)。

本文中WRF模式預(yù)報(bào)區(qū)域位于我國東北地區(qū)的松嫩平原,該區(qū)域?yàn)闇貛Т箨懶约撅L(fēng)氣候,太陽輻射充足。模式水平方向上采用單重網(wǎng)格結(jié)構(gòu),格點(diǎn)數(shù)為121×121,格點(diǎn)間距為4 km,垂直方向上分為81層,模式頂層氣壓為1 000 Pa,初始場(chǎng)使用NCEP發(fā)布的分辨率為0.25°×0.25°的再分析資料,積分步長(zhǎng)為15 s,預(yù)報(bào)時(shí)長(zhǎng)為4 h。

1.2 假設(shè)條件

影響彈丸飛行狀態(tài)的氣象要素包括氣壓、虛溫、縱風(fēng)和橫風(fēng),為了具體分析彈道解算用數(shù)值天氣預(yù)報(bào)的氣象保障效果,本文以WRF模式的短時(shí)預(yù)報(bào)為例,一方面使用均方根誤差、平均絕對(duì)誤差和誤差標(biāo)準(zhǔn)差[8]對(duì)WRF模式預(yù)報(bào)的氣象要素進(jìn)行檢驗(yàn),另一方面,基于外彈道仿真,分析WRF模式中單一氣象要素的預(yù)報(bào)誤差對(duì)彈丸飛行狀態(tài)的影響。其中,預(yù)報(bào)氣象數(shù)據(jù)和實(shí)測(cè)氣象數(shù)據(jù)描述的均為該地2019年12月14日12時(shí)的高空大氣狀態(tài),外彈道仿真則使用某型榴彈的六自由度剛體彈道模型,射角為51°,射向?yàn)?00.487°,射程約為30 km。WRF模式預(yù)報(bào)氣象要素的檢驗(yàn)結(jié)果如表1所示,WRF模式單一氣象要素的預(yù)報(bào)誤差對(duì)彈丸飛行狀態(tài)影響的分析結(jié)果如表2所示。

表1 WRF模式預(yù)報(bào)氣象要素檢驗(yàn)結(jié)果

表2 彈丸飛行狀態(tài)偏差均值表

由表1和表2可得,相比于風(fēng)的預(yù)報(bào)誤差,氣壓和虛溫的預(yù)報(bào)誤差及其離散程度以及預(yù)報(bào)誤差引起的彈道偏差均相對(duì)較小,因此,本文只將風(fēng)修正系數(shù)作為待辨識(shí)參數(shù),并忽略氣壓和虛溫的預(yù)報(bào)誤差對(duì)辨識(shí)風(fēng)修正系數(shù)的影響,故本文做出如下假設(shè):

①忽略氣壓和虛溫的預(yù)報(bào)誤差對(duì)辨識(shí)風(fēng)修正系數(shù)的影響;

②風(fēng)修正系數(shù)在辨識(shí)區(qū)間內(nèi)為常值。

1.3 彈道模型

為了從實(shí)測(cè)彈道數(shù)據(jù)和原始預(yù)報(bào)氣象數(shù)據(jù)中辨識(shí)風(fēng)修正系數(shù),需構(gòu)建彈丸運(yùn)動(dòng)模型,建立實(shí)測(cè)彈道數(shù)據(jù)與氣象數(shù)據(jù)之間的聯(lián)系?？紤]到目前彈道跟蹤觀測(cè)數(shù)據(jù)主要以位置和速度為主,本文選用質(zhì)點(diǎn)彈道模型,此彈道模型在有風(fēng)條件下的彈道方程[9]為

(1)

式中:

x、y、z,vx、vy和vz分別為彈丸在發(fā)射坐標(biāo)系下的坐標(biāo)分量和速度分量;m,d,l和JC分別為彈丸質(zhì)量、彈丸直徑、彈丸長(zhǎng)度和彈丸極轉(zhuǎn)動(dòng)慣量;Cx,C′y,m′z和m′xz分別為空氣阻力系數(shù)、升力系數(shù)導(dǎo)數(shù)、靜力矩系數(shù)導(dǎo)數(shù)和極阻尼力矩系數(shù)導(dǎo)數(shù);p,Tv,wx,wz,vs,Rd和k分別為當(dāng)?shù)貧鈮骸?dāng)?shù)靥摐?、?dāng)?shù)乜v風(fēng)、當(dāng)?shù)貦M風(fēng)、當(dāng)?shù)芈曀佟⒏煽諝鈿怏w常數(shù)和比熱比;kL和kC分別為縱風(fēng)修正系數(shù)和橫風(fēng)修正系數(shù);g,v0,s和η分別為當(dāng)?shù)刂亓铀俣?、?biāo)定初速、彈道弧長(zhǎng)和膛線纏度。

2 無跡卡爾曼濾波算法辨識(shí)風(fēng)修正系數(shù)

為了實(shí)現(xiàn)從非線性彈道模型中辨識(shí)風(fēng)修正系數(shù),本文將無跡卡爾曼濾波算法作為辨識(shí)算法。由卡爾曼濾波穩(wěn)定性原理[10]可知,如果濾波穩(wěn)定,估計(jì)值和估計(jì)均方誤差將隨濾波時(shí)間的增長(zhǎng)逐漸不受濾波初值的影響,并最終收斂,實(shí)現(xiàn)無偏估計(jì)。但是,每個(gè)估計(jì)值和每個(gè)估計(jì)均方誤差的收斂速度不同,若濾波時(shí)間較短,可能存在估計(jì)值和估計(jì)均方誤差沒有完成收斂的情況。因此,本文將原始預(yù)報(bào)氣象數(shù)據(jù)作為彈道模型的基礎(chǔ)氣象條件,采用無跡卡爾曼濾波算法對(duì)實(shí)測(cè)彈道數(shù)據(jù)進(jìn)行迭代處理,并將每次迭代所得風(fēng)修正系數(shù)和狀態(tài)變量的估計(jì)均方誤差作為下次迭代的初值,直至滿足迭代終止條件,其流程如圖1所示。

圖1 無跡卡爾曼濾波算法辨識(shí)風(fēng)修正系數(shù)流程圖

由式(1)和實(shí)測(cè)彈道數(shù)據(jù)的種類可知,無跡卡爾曼濾波算法中的觀測(cè)變量為x、y、z,vx、vy、vz,狀態(tài)變量為x、y、z,vx、vy、vz、kL和kC,即:

Z=(xyzvxvyvz)T

(2)

X=(x1x2x3x4x5x6x7x8)T=

(xyzvxvyvzkLkC)T

(3)

將式(3)代入式(1),得:

(4)

結(jié)合圖1,無跡卡爾曼濾波算法辨識(shí)風(fēng)修正系數(shù)的主要步驟如下。

①第1步。選定初值。

②第2步。無跡卡爾曼濾波計(jì)算,詳細(xì)計(jì)算步驟可參考文獻(xiàn)[10]。

③第3步。判斷是否滿足迭代終止條件,若沒有滿足,則返回第2步。為了提高迭代效率,本文設(shè)置2個(gè)迭代終止條件,條件一為達(dá)到最大迭代次數(shù)nmax,條件二為連續(xù)2次迭代所得縱風(fēng)修正系數(shù)之差εL和橫風(fēng)修正系數(shù)之差εC均小于給定閾值。

為了提升辨識(shí)結(jié)果精度,減小包括氣動(dòng)參數(shù)誤差和起始擾動(dòng)誤差在內(nèi)的其他誤差對(duì)辨識(shí)結(jié)果的影響,本文將彈丸風(fēng)洞吹風(fēng)數(shù)據(jù)作為彈道模型中氣動(dòng)參數(shù)的真值,并將多組辨識(shí)結(jié)果的平均值作為最終的辨識(shí)結(jié)果。

3 數(shù)值計(jì)算與分析

為了檢驗(yàn)無跡卡爾曼濾波算法辨識(shí)風(fēng)修正系數(shù)的效果,本文基于某實(shí)測(cè)彈道數(shù)據(jù)和WRF模式原始預(yù)報(bào)氣象數(shù)據(jù)對(duì)其進(jìn)行數(shù)值分析,相關(guān)參數(shù)設(shè)計(jì)如下。

①實(shí)測(cè)彈道數(shù)據(jù)的位置誤差(DRMS,68%):水平10 m,垂直10 m;速度精度(DRMS,68%):水平0.2 m/s,垂直0.5 m/s;

②彈道數(shù)據(jù)選擇距地面5 200～6 800 m的降弧段彈道數(shù)據(jù);

③WRF模式參數(shù)設(shè)置與第1.1節(jié)中的參數(shù)設(shè)置相同;

④無跡卡爾曼濾波算法迭代終止條件一為最大迭代次數(shù)nmax=200,迭代終止條件二為連續(xù)2次迭代差值εL<0.000 1且εC<0.000 1。

基于上述參數(shù),本文以其中一組實(shí)測(cè)彈道數(shù)據(jù)為例,在辨識(shí)過程中x、y、z的估計(jì)均方誤差隨迭代次數(shù)n(取前10次,下同)的變化如圖2所示,vx、vy、vz的估計(jì)均方誤差隨迭代次數(shù)的變化如圖3所示,kL、kC的估計(jì)均方誤差隨迭代次數(shù)的變化如圖4所示,kL、kC隨迭代次數(shù)的變化如圖5所示。

圖2 迭代過程中的位置估計(jì)均方誤差

圖3 迭代過程中的速度估計(jì)均方誤差

圖4 迭代過程中的風(fēng)修正系數(shù)估計(jì)均方誤差

將3組辨識(shí)結(jié)果取平均值作為該辨識(shí)區(qū)段內(nèi)的辨識(shí)結(jié)果,即辨識(shí)的縱風(fēng)修正系數(shù)為2.084 056,橫風(fēng)修正系數(shù)為1.254 748。實(shí)測(cè)的風(fēng)修正系數(shù)如表3所示,使用加權(quán)平均值(權(quán)重由彈丸在每個(gè)高度段內(nèi)的飛行時(shí)間確定,下同)代表實(shí)測(cè)的風(fēng)修正系數(shù),即實(shí)測(cè)的縱風(fēng)修正系數(shù)為1.944 045,橫風(fēng)修正系數(shù)為1.206 837。

表3 實(shí)測(cè)的風(fēng)修正系數(shù)

將實(shí)測(cè)的風(fēng)修正系數(shù)和辨識(shí)的風(fēng)修正系數(shù)進(jìn)行對(duì)比可得,辨識(shí)的縱風(fēng)修正系數(shù)和橫風(fēng)修正系數(shù)的相對(duì)誤差分別為7.20%和3.97%,產(chǎn)生上述誤差的原因可能如下:一是辨識(shí)過程中未考慮氣壓和虛溫的預(yù)報(bào)誤差對(duì)辨識(shí)風(fēng)修正系數(shù)的影響,將氣壓和虛溫的預(yù)報(bào)誤差折算入風(fēng)修正系數(shù)當(dāng)中;二是彈道模型誤差,結(jié)合六自由度剛體彈道模型可知,狀態(tài)變量還受彈丸姿態(tài)角等因素的影響,而質(zhì)點(diǎn)彈道模型未考慮這些因素;三是實(shí)測(cè)氣象數(shù)據(jù)與實(shí)測(cè)彈道數(shù)據(jù)時(shí)空不匹配。

原始預(yù)報(bào)氣象數(shù)據(jù)和實(shí)測(cè)氣象數(shù)據(jù)描述的均為同一時(shí)刻且同一地點(diǎn)的高空大氣狀態(tài),利用前文辨識(shí)所得的縱風(fēng)和橫風(fēng)的修正系數(shù)對(duì)原始預(yù)報(bào)氣象數(shù)據(jù)中的風(fēng)進(jìn)行修正,得到修正預(yù)報(bào)氣象數(shù)據(jù),3種氣象數(shù)據(jù)中風(fēng)的加權(quán)平均值如表4所示。

表4 彈道區(qū)段內(nèi)風(fēng)的加權(quán)平均值

由表4可得,原始預(yù)報(bào)氣象數(shù)據(jù)中的縱風(fēng)和橫風(fēng)經(jīng)修正后,其誤差分別減小85.19%和79.27%。

分別使用原始預(yù)報(bào)氣象數(shù)據(jù)和修正預(yù)報(bào)氣象數(shù)據(jù)對(duì)位于距地面5 200～6 800 m的降弧段彈道進(jìn)行重構(gòu),2種重構(gòu)的彈道段和濾波處理的實(shí)測(cè)彈道段在x-y平面的投影如圖6所示,在x-z平面的投影如圖7所示,在z-y平面的投影如圖8所示,2種重構(gòu)的彈道段與濾波處理的實(shí)測(cè)彈道段間的狀態(tài)均方根誤差如表5所示。

圖6 重構(gòu)彈道段在x-y平面的投影

圖7 重構(gòu)彈道段在x-z平面的投影

圖8 重構(gòu)彈道段在z-y平面的投影

表5 重構(gòu)彈道段的狀態(tài)均方根誤差

由圖6～圖8及表5可得,修正預(yù)報(bào)氣象重構(gòu)的彈道段相比原始預(yù)報(bào)氣象重構(gòu)的彈道段而言,更加接近于濾波處理的實(shí)測(cè)彈道段,因此,修正預(yù)報(bào)氣象比原始預(yù)報(bào)氣象對(duì)真實(shí)氣象的描述更為準(zhǔn)確,從而說明辨識(shí)所得風(fēng)修正系數(shù)的精度較高以及修正彈道解算用數(shù)值天氣預(yù)報(bào)的預(yù)報(bào)誤差是必要的。

4 結(jié)束語

本文將原始預(yù)報(bào)氣象數(shù)據(jù)作為基礎(chǔ)氣象條件,采用無跡卡爾曼濾波算法從實(shí)測(cè)彈道數(shù)據(jù)中辨識(shí)風(fēng)修正系數(shù),通過對(duì)實(shí)際辨識(shí)結(jié)果進(jìn)行分析,得出以下結(jié)論:該方法辨識(shí)所得風(fēng)修正系數(shù)具有較高的精度,可作為彈道解算用數(shù)值天氣預(yù)報(bào)的預(yù)報(bào)誤差,實(shí)現(xiàn)對(duì)臨近天氣預(yù)報(bào)的數(shù)值修正,滿足工程應(yīng)用需求;無跡卡爾曼濾波算法在辨識(shí)風(fēng)修正系數(shù)的過程中兼有濾波功能,能夠在一定程度上減小噪聲對(duì)辨識(shí)結(jié)果的影響。本文研究結(jié)果可以為提升彈道解算用數(shù)值天氣預(yù)報(bào)的預(yù)報(bào)質(zhì)量提供一定的理論參考。