在線彈道參數(shù)濾波與辨識方法分析

李小元,王中原,常思江

(1.海軍研究院,北京 100161;2.南京理工大學(xué) 能源與動力工程學(xué)院,江蘇 南京 210094)

炮彈在發(fā)射飛行過程中,存在火炮初始擾動、初速或然誤差、炮彈結(jié)構(gòu)參數(shù)誤差、氣象條件誤差等,造成射彈散布。

彈道修正技術(shù)是指對發(fā)射出去的炮彈的一段飛行彈道參數(shù)進行實時測量,同預(yù)定彈道參數(shù)進行比較,進行彈道偏差計算,形成控制指令,通過彈上控制機構(gòu)作用調(diào)節(jié)彈道向預(yù)定彈道逼近,從而大幅度減小炮彈散布的一種技術(shù)。如果僅對縱向彈道進行調(diào)節(jié),稱為一維彈道修正彈。當彈上控制機構(gòu)采用阻力環(huán)裝置來調(diào)節(jié)阻力,則為打遠修近、阻力環(huán)體制的一維彈道修正彈,它是目前廣泛采用、結(jié)構(gòu)簡潔的一種一維彈道修正技術(shù)[1-3]。

一維彈道修正技術(shù)在實際中要取得好的應(yīng)用效果,一個重要問題是如何根據(jù)實測的一段飛行彈道參數(shù)快速準確地辨識出該彈的彈道特征參數(shù),進而準確地預(yù)報后續(xù)彈道、落點,并準確地知道這發(fā)彈同預(yù)定彈道的偏差,從而實施彈道修正。目前,采用卡爾曼濾波、特征參數(shù)辨識方法是一種較為有效的方法。文獻[4-5]分別采用卡爾曼濾波、速度-時間序列匹配等方法對該問題進行了研究。文獻[6]通過建立彈道濾波方程,采用卡爾曼濾波法外推彈道,對彈道修正彈的落點進行了預(yù)測。為了進一步提高彈箭密集度,文獻[7]比較了線性卡爾曼模型和擴展卡爾曼模型在GPS量測數(shù)據(jù)下的差異,并分析了2種模型在不同條件下的優(yōu)缺點。但是上述文獻在數(shù)據(jù)測量精度、基礎(chǔ)阻力系數(shù)變化等因素對預(yù)報彈道落點的影響方面研究較少。

在卡爾曼濾波、特征參數(shù)辨識方法中,在對應(yīng)整個可提供的實時測量彈道段上,分為預(yù)報和校正2個環(huán)節(jié)交替實施,核心目的是快速、準確辨識出廣義氣動力系數(shù)。在預(yù)報環(huán)節(jié),彈的阻力系數(shù)狀況對預(yù)報結(jié)果極為敏感,影響很大;在校正環(huán)節(jié),雷達實時測量的彈道參數(shù)的精度對校正結(jié)果影響極大。它們對最終辨識出的廣義氣動力系數(shù)影響很大,直接影響在線彈道預(yù)報精度。本文對阻力環(huán)修正體制下的這些問題進行分析、討論,給出相關(guān)建議。

1 一維彈道修正彈的飛行彈道控制方程

針對本文研究的一維彈道修正彈,根據(jù)其工作原理,可建立一維修正彈道方程:

(1)

控制方程為

(2)

式中:Fx2,Fy2,Fz2為作用在彈丸上的合力F在彈道坐標系x2y2z2中的投影,在計算Fx2,Fy2時需要使用阻力系數(shù)符合系數(shù)kcx和升力系數(shù)符合系數(shù)kcy,上述2個系數(shù)為待辨識的參數(shù);Mξ,Mη,Mζ為作用在彈丸質(zhì)心處的合力矩M在彈軸坐標系ξηζ中的投影;MD為阻力環(huán)附加的極阻尼力矩;FD為阻力環(huán)張開后增加的阻力;ΔCx為阻力系數(shù)增量;v為速度;θa為彈道傾角;ψ2為彈道偏角;ωξ,ωη,ωζ為彈丸角速度在彈軸坐標系ξηζ中的分量;φa為彈軸高低擺動角;φ2為彈軸側(cè)向擺動角;γ為滾轉(zhuǎn)角;x,y,z為彈丸在地面坐標內(nèi)的位置分量;JA,JC分別為極轉(zhuǎn)動慣量和赤道轉(zhuǎn)動慣量;m為彈丸質(zhì)量;ρ為大氣密度;S為彈丸特征面積。各坐標系的定義參見文獻[8]。

2 卡爾曼濾波與特征參數(shù)辨識的彈道預(yù)報模型和方法

2.1 卡爾曼濾波方法

如前所述,如何根據(jù)實測的一段飛行彈道參數(shù)快速、準確地辨識出該彈的彈道特征參數(shù)是一個核心問題。本文選取阻力系數(shù)符合系數(shù)kcx和升力系數(shù)符合系數(shù)kcy作為彈道特征參數(shù),分別用于縱向射程與橫向距離的彈道、落點預(yù)報。為了能夠利用一段實測彈道數(shù)據(jù)盡可能準確地估計出這2個彈道特征參數(shù),本文采用卡爾曼濾波方法,將kcx和kcy作為狀態(tài)變量進行估計?？柭鼮V波方法的優(yōu)點是考慮了物理系統(tǒng)本身的動態(tài)模型和測量系統(tǒng)的量測誤差,利用從某時刻k-1以前直到現(xiàn)在的測量數(shù)據(jù)Y(j個)逐次遞推,可估計出k時刻最優(yōu)的系統(tǒng)狀態(tài)X?？柭鼮V波基本方程簡介如下。

預(yù)測方程:

(3)

(4)

(5)

濾波方程:

(6)

卡爾曼濾波增益:

(7)

協(xié)方差矩陣:

Pk=(I-KkHk)Pk/k-1

(8)

式中:Δt為卡爾曼濾波的采樣間隔;Φk為基本矩陣,H為量測矩陣,Qk為過程噪聲方差矩陣,Rk為量測噪聲方差矩陣,其他各符號定義參見文獻[9]。

下面以上述卡爾曼濾波方程為基礎(chǔ),結(jié)合具體的彈道模型,構(gòu)建一維修正彈的彈道濾波模型。

2.2 彈道濾波模型

由于在利用卡爾曼濾波辨識彈道特征參數(shù)的過程中,不涉及阻力環(huán)機構(gòu)作用。因此,在構(gòu)建一維修正彈彈道濾波模型時,令方程(1)中FD=0,并將參數(shù)kcx和kcy分別置于阻力系數(shù)和升力系數(shù)之前作為乘子,并補充如下微分方程:

(9)

則狀態(tài)變量取為

X=(vθaψ2xyzωξωηωζφaφ2γkcxkcy)T

(10)

將方程(1)改寫成狀態(tài)空間形式為

(11)

式中:fi(i=1,2,…,14)為方程(1)和方程(9)對應(yīng)的右端函數(shù)。

對方程(11)進行線性化,可得:

(12)

式中:F為動力矩陣,其元素為fi(i=1,2,…,14)對每個狀態(tài)變量的偏導(dǎo)數(shù),由于具體形式冗長,為節(jié)約篇幅,這里略去。

在求出了動力矩陣F的具體表達式后,可求得彈道濾波系統(tǒng)的基本矩陣Φk為

Φk≈I14×14+F14×14·Δt

(13)

式中:I14×14為單位矩陣。

假設(shè)可直接或間接測得的彈道參數(shù)為

Z=(vθaψ2xyz)T

(14)

則彈道濾波系統(tǒng)的量測方程可以表達為

Z=h(X)+n=H·X+n

(15)

式中:n為量測噪聲。

由于狀態(tài)變量X為14×1維矩陣,Z為6×1維矩陣,故量測矩陣H為6×14維矩陣,n為6×1維矩陣。

(16)

量測噪聲n的方差矩陣為

(17)

式中:σ為測量值的標準差。

另外,與狀態(tài)方程相對應(yīng)的還有過程噪聲方差矩陣Qk,由于很難描述彈道模型與實際彈丸飛行狀態(tài)的差異,故一般取Qk=O14×14。

數(shù)值迭代求解如下形式的Ricatti方程:

(18)

式中:Pk/k-1為k時刻數(shù)據(jù)更新前的狀態(tài)估計值協(xié)方差矩陣;Pk為數(shù)據(jù)更新后的狀態(tài)估計值協(xié)方差矩陣,相應(yīng)地,Pk-1為上一時刻的狀態(tài)估計值協(xié)方差矩陣;Kk為卡爾曼增益矩陣。

為了使卡爾曼濾波系統(tǒng)正常工作,狀態(tài)估計值協(xié)方差矩陣的初始值P0根據(jù)量測噪聲的標準偏差來選取,有:

P0=diag(P0,11,P0,22,…,P0,1414)

(19)

卡爾曼濾波的基本方程為

(20)

2.3 彈道預(yù)報方法

根據(jù)卡爾曼濾波方程(20)可以實時估計出彈道方程組對應(yīng)的狀態(tài)變量,將測量數(shù)據(jù)段最后一點的狀態(tài)變量最優(yōu)估計值作為彈道預(yù)報的起點,對彈道方程組進行數(shù)值積分直至落點,可預(yù)報出阻力環(huán)機構(gòu)全彈道不作用條件下的該發(fā)彈的無控彈道及落點。將事先獲取的阻力環(huán)阻力系數(shù)ΔCx代入彈道方程組中(即FD≠0),并預(yù)設(shè)一個阻力環(huán)啟動作用時刻,可預(yù)報出該發(fā)彈在某一時刻張開阻力環(huán)后的修正彈道及落點。如果將無控彈道預(yù)測落點與目標彈的射程相減,取得該發(fā)彈的射程修正量,利用該修正量進行迭代計算,可很容易求出該射程修正量對應(yīng)的阻力環(huán)機構(gòu)啟動時刻。

3 在線彈道預(yù)報分析

為了驗證上述模型與方法的正確性和有效性,本節(jié)以某大口徑一維彈道修正彈的實際數(shù)據(jù)為例,開展在線彈道預(yù)報分析。在該一維彈道修正彈的某次試驗中,連續(xù)射擊6發(fā)無控彈,采用兩部彈道跟蹤雷達同時對每發(fā)彈丸實時跟蹤,分別測得對應(yīng)的飛行彈道數(shù)據(jù)。這兩部雷達的跟蹤精度有差異,其中一部為高精度彈道跟蹤雷達(記為HA雷達),而另一部的跟蹤精度相對較低(記為LA雷達)。采用不同精度雷達同時跟蹤的目的是為了研究測量數(shù)據(jù)精度對彈道預(yù)報的影響。需要說明的是,由于一維彈道修正彈的彈道預(yù)報主要關(guān)注縱向距離x,故下面分析均對應(yīng)預(yù)報縱向距離x。

具體分析包括幾個方面:研究不同阻力系數(shù)對彈道預(yù)報精度的影響;分析測量飛行彈道數(shù)據(jù)對彈道預(yù)報的影響。

3.1 在線彈道預(yù)報結(jié)果

利用HA雷達實測5 s數(shù)據(jù),包括彈丸速度和位置,進行在線彈道預(yù)報,所得結(jié)果如圖1～圖3和表1所示。

其中,圖1～圖3分別為其中1發(fā)彈丸的位置、速度及阻力系數(shù)符合系數(shù)kcx濾波曲線,表1為6發(fā)彈丸預(yù)報射程與實測射程之差。從圖1～圖3可以看出,采用卡爾曼濾波方法得到了彈丸位置、速度及阻力系數(shù)符合系數(shù)的最優(yōu)估計值。由表1所示的在線射程預(yù)報結(jié)果看,6發(fā)預(yù)報射程與實測射程之差出現(xiàn)了正負號,即對于這一組射彈,有的射程預(yù)報值大于射程實測值,有的則小于實測值。這表明由于各種隨機誤差因素的存在,在線彈道預(yù)報誤差本身也存在隨機性,而一維彈道修正彈的主要功能是提高射程密集度,因此應(yīng)盡可能減小彈道預(yù)報誤差的跳差范圍,即盡量控制極差。從數(shù)值上看,試驗彈的平均射程為45 700 m,而6發(fā)彈的射程預(yù)報極差為73 m(約為平均射程的0.16%),射程預(yù)報誤差平均值為21 m(不到平均射程的0.05%),能夠滿足實際工程應(yīng)用的要求。

圖1 彈丸位置實測值與最優(yōu)估計值

圖2 彈丸速度實測值與最優(yōu)估計值

圖3 參數(shù)kcx的最優(yōu)估計結(jié)果

表1 6發(fā)彈丸預(yù)報射程與實測射程之差

3.2 基礎(chǔ)阻力系數(shù)變化對彈道預(yù)報的影響

在彈丸阻力系數(shù)標準值基礎(chǔ)上變化±2.5%和±5%,將變化后的阻力系數(shù)代入前面建立的在線彈道預(yù)報模型中,所得預(yù)報射程與實測射程之差ΔXc列于表2中。

表2 不同基礎(chǔ)阻力系數(shù)對應(yīng)的預(yù)報射程與實測射程之差

表2中結(jié)果說明,彈丸阻力系數(shù)標準值變化對預(yù)報結(jié)果有一定影響。按照相對于標準值的變化百分比,預(yù)報射程誤差呈線性變化。從一組數(shù)據(jù)來看,預(yù)報射程誤差相對于標準值條件下的預(yù)報誤差產(chǎn)生了一個系統(tǒng)偏差,6發(fā)彈丸預(yù)報射程誤差的極差變化很小。從這個意義上講,彈丸阻力系數(shù)標準值發(fā)生變化時,主要影響在線彈道預(yù)報結(jié)果的準確度,而不影響反映彈丸發(fā)與發(fā)之間差異的散布,即對一維修正落點的密集度無影響。從武器使用角度,應(yīng)兼顧準確度和密集度,故有必要采用盡可能準確的阻力系數(shù)。由于本文通過大量試驗取得了較為精確的彈丸阻力系數(shù),因此預(yù)報準確度也相對較高。

3.3 測量數(shù)據(jù)對彈道預(yù)報的影響

試驗中LA雷達斜距離標準差比HA雷達斜距離標準差大2 m,LA雷達高低角標準差比HA雷達高低角標準差大0.65 mrad,LA雷達方位角標準差比HA雷達方位角標準差大0.5 mrad。表3所示為采用HA和LA兩部雷達跟蹤數(shù)據(jù)預(yù)報出的射程值,其他條件參數(shù)均取為相同。Xm為實測射程,Xc,HA為HA雷達預(yù)報射程,Xc,LA為LA雷達預(yù)報射程,ΔXc,HA=Xc,HA-Xm,ΔXc,LA=Xc,LA-Xm。

表3 利用兩部雷達跟蹤數(shù)據(jù)預(yù)報的射程

從平均值上看,基于LA雷達跟蹤數(shù)據(jù)的預(yù)報誤差要小于基于HA雷達跟蹤數(shù)據(jù)的預(yù)報誤差,特別是對于2號彈、3號彈以及6號彈。但對于一維彈道修正彈,其原理主要是減小彈丸的縱向散布,故對于在線彈道預(yù)報而言,更重要的是盡可能減小一組射彈射程預(yù)報誤差的極差。這一極差主要受跟蹤數(shù)據(jù)精度的影響,由表3可知,HA雷達對應(yīng)的射程預(yù)報誤差極差為73 m,而LA雷達對應(yīng)的射程預(yù)報誤差極差為116 m。因此,應(yīng)盡可能采用跟蹤精度高的雷達。

測量數(shù)據(jù)的長度對預(yù)報結(jié)果也有影響。前面計算中,都是采用了5 s跟蹤數(shù)據(jù)。下面針對HA雷達數(shù)據(jù),以1號彈為例,研究不同數(shù)據(jù)長度對預(yù)報結(jié)果的影響,如圖4所示。

圖4 不同數(shù)據(jù)長度對應(yīng)的射程預(yù)報誤差

由圖4可以看出跟蹤數(shù)據(jù)長度對預(yù)報結(jié)果的影響。當數(shù)據(jù)時長大于5 s,預(yù)報誤差略有波動,但基本上維持在25 m左右,數(shù)據(jù)時長越長對預(yù)報誤差的影響較小,而數(shù)據(jù)越多,在線計算量也越大;當數(shù)據(jù)時長小于5 s,預(yù)報誤差隨著數(shù)據(jù)時長的減小而迅速增大,3 s時長對應(yīng)的預(yù)報誤差已接近70 m。由于圖4所示是1號彈對應(yīng)的結(jié)果,如果是對于6號彈,3 s時長對應(yīng)的預(yù)報誤差接近150 m,對修正彈最終密集度指標的實現(xiàn)極為不利。上述結(jié)果表明,當綜合考慮在線彈道預(yù)報的精度和計算量(實時性),選取5 s數(shù)據(jù)是比較合適的。

3.4 關(guān)鍵參數(shù)對彈道預(yù)報的影響

前面所建基于卡爾曼濾波的在線彈道預(yù)報模型中,參數(shù)kcx也是作為一個狀態(tài)變量。實際研究中發(fā)現(xiàn),該狀態(tài)變量對應(yīng)的協(xié)方差矩陣元素初始值Pk,0對參數(shù)kcx的最優(yōu)估計結(jié)果有影響,并且較為敏感,進而影響射程預(yù)報的結(jié)果。本節(jié)仍采用HA雷達跟蹤數(shù)據(jù),研究不同Pk,0值對預(yù)報結(jié)果的影響,結(jié)果如表4所示。表中,ΔXc為預(yù)報射程與實測射程之差。

表4 不同Pk,0值對應(yīng)的預(yù)報射程與實測射程之差

圖5為不同Pk,0值對應(yīng)的射程預(yù)報誤差極差(ΔXc,m)曲線圖。

圖5 不同Pk,0值對應(yīng)的射程預(yù)報誤差極差

從表4和圖5中不難看出,盡管Pk,0的量值不大,但對預(yù)報結(jié)果卻具有較大影響。前面研究表明,對于一維彈道修正彈,一組射程預(yù)報誤差的極差大小可在一定程度上反映在線彈道預(yù)報的性能。當Pk,0>0.01,隨著Pk,0值的增大,射程預(yù)報誤差的平均值和極差都增大;當Pk,0<0.01,隨著Pk,0值的減小,射程預(yù)報誤差的平均值和極差也增大。圖5所示極差曲線中,顯然Pk,0=0.01是一個最優(yōu)值。需要指出的是,對于不同彈丸、不同測量設(shè)備或不同模型,該最優(yōu)值也不盡相同,可采用試驗手段繪制類似圖5的極差曲線,從而確定Pk,0的最優(yōu)值。工程實踐表明,這一方法是比較有效的。

4 結(jié)束語

本文針對阻力環(huán)體制的一維彈道修正彈,建立了飛行控制彈道模型,并基于卡爾曼濾波理論建立了彈道濾波模型,給出了彈道預(yù)報的方法,并通過算例對本文研究結(jié)果的正確性和有效性進行了分析和驗證。通過對實際工程應(yīng)用結(jié)果的分析, 提出以下建議:為了提高無控彈道和修正彈道的預(yù)報精度,應(yīng)盡可能準確地獲取全彈的阻力系數(shù),系數(shù)值應(yīng)覆蓋全彈道馬赫數(shù)。在確定某一實時測量彈道參數(shù)的基礎(chǔ)上,進行若干發(fā)炮射試驗,通過反復(fù)調(diào)整預(yù)報-校正分段、測量參數(shù)方差同實測彈道的對比,確定出較佳匹配(協(xié)方差、分段等)關(guān)系,來提高“廣義阻力系數(shù)”的辨識效果,提高后續(xù)彈道預(yù)報精度。