鄭宇程,易文俊,余春華,管 軍
(1.南京理工大學(xué) 瞬態(tài)物理國家重點實驗室, 南京 210094; 2.南京理工大學(xué) 理學(xué)院, 南京 210094)
火炮在現(xiàn)代戰(zhàn)爭中具有重要的使命,在對敵火力壓制過程中能夠起到關(guān)鍵性作用。獲取傳統(tǒng)高速旋轉(zhuǎn)穩(wěn)定彈丸準(zhǔn)確的氣動參數(shù),對于提高火炮射表精度、減小落點散布、增強打擊精度具有重要的意義。目前,國內(nèi)外獲取彈丸氣動參數(shù)的方法主要有三種[1]:第一種是理論計算,第二種是風(fēng)洞吹風(fēng),第三種是利用彈丸的自由飛行數(shù)據(jù)對彈丸的氣動參數(shù)進行離線辨識。理論計算法雖然簡單,但由于模型中的未建模因素和不確定因素導(dǎo)致計算結(jié)果存在誤差。風(fēng)洞吹風(fēng)法作用于彈丸模型,結(jié)果較為準(zhǔn)確但成本較高,不能精準(zhǔn)模擬高速旋轉(zhuǎn)狀態(tài)。利用彈丸自由飛行數(shù)據(jù)辨識彈丸的氣動參數(shù),不僅符合實際情況,還能根據(jù)辨識結(jié)果,及時對彈丸運動狀態(tài)進行調(diào)整,從而提高炮彈的打擊精度。
用于參數(shù)辨識的方法通常有遞推最小二乘法、遞推極大似然法、卡爾曼(kalman)濾波法等[2-5]。近年來,許多學(xué)者對此進行了研究。史繼剛[2]基于粒子群初值選取的牛頓迭代優(yōu)化算法辨識彈丸的零升阻力系數(shù),該方法基于最大似然準(zhǔn)則,相比于卡爾曼濾波法實現(xiàn)更加困難;管軍等[3]提出一種新的自適應(yīng)混沌變異粒子群算法求解該準(zhǔn)則下的氣動參數(shù)最優(yōu)解,得到彈丸的氣動參數(shù),但是在工程上難以實現(xiàn);Rogers等[4]提出了一種基于證據(jù)理論的參數(shù)估計方法,偏于理論計算;史金光等[6]利用擴展卡爾曼濾波法對彈道修正彈的阻力和升力符合系數(shù)進行了辨識,并且由此對后續(xù)彈道進行修正,然而該方法技術(shù)要求較高,難以在實際應(yīng)用中實現(xiàn);楊靖等[7]利用改進的混雜擴展卡爾曼濾波法對彈丸的阻力系數(shù)進行辨識,其對EKF的改進不大。
相比于最小二乘法和極大似然法,卡爾曼濾波法應(yīng)用最為廣泛。在處理非線性問題上,擴展性卡爾曼濾波(EKF)和無跡卡爾曼濾波(UKF)已經(jīng)成功地應(yīng)用于各種參數(shù)辨識問題[8-11]。本文主要利用擴展卡爾曼濾波和無跡卡爾曼濾波對彈丸的阻力系數(shù)進行辨識,并進行了對比研究。為了進行氣動參數(shù)的辨識,本文將未知的氣動參數(shù)看成系統(tǒng)的狀態(tài)量,分別用EKF方法和UKF方法對氣動參數(shù)進行辨識。為了檢驗實驗的正確性,又分為氣動參數(shù)不變和氣動參數(shù)變化兩種情況進行研究。
參數(shù)辨識的本質(zhì)問題就是如何通過輸出數(shù)據(jù)準(zhǔn)確辨識相關(guān)輸入?yún)?shù)或系統(tǒng)參數(shù)。如何選取系統(tǒng)的模型是重要的步驟之一。
阻力系數(shù)是傳統(tǒng)無控旋轉(zhuǎn)彈丸所有氣動參數(shù)中最重要的氣動參數(shù)之一,其直接影響了彈丸的射程和落點散布。本文的主要研究內(nèi)容是彈丸阻力系數(shù)辨識,雖然外彈道學(xué)中的彈道模型多種多樣,但2D彈道模型已經(jīng)能夠描述阻力系數(shù)對彈道性能的絕大部分影響。且2D彈道模型簡單、便捷,能夠直接揭示阻力系數(shù)和彈道參數(shù)之間的關(guān)系。因此,本文選用2D彈道模型作為系統(tǒng)的模型。
2D彈道模型[12]的方程為
(1)
其中,ρ為當(dāng)前高度的大氣密度,S為參考面積,m為彈體質(zhì)量。
卡爾曼濾波法最初只用于離散線性系統(tǒng),后來經(jīng)過改進,產(chǎn)生了許多種不同的方法。擴展性卡爾曼濾波(EKF)和無跡卡爾曼濾波(UKF)是應(yīng)用得比較廣泛的兩種方法。EKF方法在線性化處理的過程中對原系統(tǒng)進行泰勒展開,略去了二階以及二階以上的高階項,只保留了線性項。
擴展卡爾曼濾波的基本原理為:
1) 狀態(tài)一步預(yù)測方程
(2)
2) 求一步預(yù)測均方誤差
(3)
式中:Qk-1是系統(tǒng)噪聲的協(xié)方差矩陣,Γk-1是系統(tǒng)噪聲驅(qū)動陣,Φk,k-1為一步轉(zhuǎn)移矩陣,計算式如下:
Φk,k-1=I+A·Δt
(4)
3) 求濾波增益矩陣
(5)
4) 估計均方誤差
Pk=[I-KkHk]Pk,k-1
(6)
5) 得到最新一步的狀態(tài)估計
(7)
式中,zk為k時刻的觀測向量。
由于EKF是一種非線性估計問題的線性化方法,其線性化過程會產(chǎn)生一定的誤差,且其需要進行復(fù)雜的雅克比矩陣求解;而UKF通過UT變換,直接對非線性問題進行處理,不需進行復(fù)雜的雅克比矩陣求解,且其精度可以達到二階泰勒展開的精度。因此,本文同樣使用UKF對阻力系數(shù)進行辨識,并同EKF進行比較研究,其基本計算步驟為:
1) 初始化濾波參數(shù)
(8)
式中:n為增廣狀態(tài)向量的維數(shù);α是很小的一個正數(shù),其取值范圍為10-4~1。本文取為0.01;κ的取值為3-n;β的值與x的分布形式有關(guān),對于正態(tài)分布的情況,β取2最優(yōu)。
2) 計算k-1時刻的sigma樣本點
(9)
3) 計算k時刻的一步預(yù)測值
(10)
4) 計算k時刻的一步預(yù)測樣本點
(11)
5) 計算P(XZ)k,k-1和P(ZZ)k,k-1
(12)
6) 更新濾波值
(13)
由于待辨識參數(shù)是阻力系數(shù)CD,本文通過增廣的方法將其添加到式(1),以建立包含待辨識參數(shù)的狀態(tài)方程。
利用系數(shù)凍結(jié)法,設(shè)在小區(qū)間內(nèi)彈丸的阻力系數(shù)是常數(shù),因此有下式成立:
(14)
則增廣狀態(tài)變量x為
x=[u,θ,x,y,CD]T
(15)
增廣之后的狀態(tài)方程可寫為
(16)
ω為系統(tǒng)噪聲向量。
通常情況下,彈道跟蹤雷達只能得到彈丸的速度和位置信息,因此量測向量取為
z=[u,x,y]T
(17)
量測方程為
z=h(x)+v
(18)
v為觀測噪聲向量。
彈丸的物理參數(shù)為:彈重45.5 kg;彈徑0.155 m;彈長0.909 m。氣象條件為標(biāo)準(zhǔn)氣象條件[13]。射擊初始值為:初速930 m/s,射角45°。
狀態(tài)初值為
x0=[930,45°,0,0,CD]Τ
(19)
初始狀態(tài)誤差協(xié)方差矩陣為
(20)
系統(tǒng)噪聲的協(xié)方差矩陣為
(21)
量測噪聲的協(xié)方差矩陣為
(22)
算例1
假設(shè)CD為常數(shù),其真實值為0.35。用CD為0.35計算出來的理論值(速度,位置坐標(biāo))加上高斯白噪聲v來模擬觀測值。對于CD為常數(shù)的情況,用EKF和UKF辨識結(jié)果如圖1~圖2所示。從圖1和圖2可以看出,兩種方法辨識的CD最終都收斂在0.35,但UKF具有更快的收斂速度。
算例2
實際情況中,阻力系數(shù)CD是馬赫數(shù)的函數(shù),本算例所采用的數(shù)值如表1所示。
根據(jù)變化的CD計算出的理論值,加上高斯白噪聲v的值,將其用來模擬觀測值。將EKF與UKF辨識的馬赫數(shù)與真實的馬赫數(shù)之差畫成曲線在一起對比,如圖3所示。
圖1 EKF辨識CD不變
圖2 UKF辨識CD不變
Ma0.50.751.01.251.51.75CD0.159 20.1570.299 50.321 40.296 20.275 7Ma2.02.252.52.753.0CD0.2580.239 90.225 70.207 90.191 5
圖3 馬赫數(shù)辨識誤差曲線
圖3不但表示了CD隨馬赫數(shù)變化的情況,而且可看出UKF辨識的馬赫數(shù)誤差較小,EKF辨識的馬赫數(shù)誤差較大。但是最大的誤差不超過10-3,為原馬赫數(shù)的0.05%,因此認為兩種方法辨識的馬赫數(shù)是準(zhǔn)確的。CD雖然是馬赫數(shù)的函數(shù),但是為了便于比較研究,將真實CD曲線、EKF辨識CD曲線和UKF辨識CD曲線對比,如圖4所示。
從圖4可以看出,前后兩端三者差別不大,現(xiàn)在選取中間20~50 s的時間段,參見圖5。
圖4 EKF與UKF辨識結(jié)果對比(1)
圖5 EKF與UKF辨識結(jié)果對比(2)
圖5中,從上至下分別是真實CD曲線、UKF辨識曲線和EKF辨識曲線。經(jīng)過比較,UKF辨識曲線的精度較EKF高,兩者的差距在0.002 5左右。UKF辨識的誤差為0.8%,EKF辨識的誤差為1.6%。由此可見,參數(shù)辨識的結(jié)果具有可行性。
1)CD不變的情況下,EKF方法與UKF方法都能快速收斂到真值,且UKF具有較快的收斂速度。
2)CD隨馬赫數(shù)變化的情況下,EKF方法與UKF方法都能有效辨識出阻力系數(shù),且UKF較EKF具有更高的辨識精度。