在線彈道參數濾波與辨識方法分析

李小元,王中原,常思江

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

炮彈在發射飛行過程中,存在火炮初始擾動、初速或然誤差、炮彈結構參數誤差、氣象條件誤差等,造成射彈散布。

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

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

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

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

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

(1)

控制方程為

(2)

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

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

2.1 卡爾曼濾波方法

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

預測方程:

(3)

(4)

(5)

濾波方程:

(6)

卡爾曼濾波增益:

(7)

協方差矩陣:

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

(8)

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

下面以上述卡爾曼濾波方程為基礎,結合具體的彈道模型,構建一維修正彈的彈道濾波模型。

2.2 彈道濾波模型

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

(9)

則狀態變量取為

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

(10)

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

(11)

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

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

(12)

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

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

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

(13)

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

假設可直接或間接測得的彈道參數為

Z=(vθaψ2xyz)T

(14)

則彈道濾波系統的量測方程可以表達為

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

(15)

式中:n為量測噪聲。

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

(16)

量測噪聲n的方差矩陣為

(17)

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

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

數值迭代求解如下形式的Ricatti方程:

(18)

式中:Pk/k-1為k時刻數據更新前的狀態估計值協方差矩陣;Pk為數據更新后的狀態估計值協方差矩陣,相應地,Pk-1為上一時刻的狀態估計值協方差矩陣;Kk為卡爾曼增益矩陣。

為了使卡爾曼濾波系統正常工作,狀態估計值協方差矩陣的初始值P0根據量測噪聲的標準偏差來選取,有:

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

(19)

卡爾曼濾波的基本方程為

(20)

2.3 彈道預報方法

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

3 在線彈道預報分析

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

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

3.1 在線彈道預報結果

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

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

圖1 彈丸位置實測值與最優估計值

圖2 彈丸速度實測值與最優估計值

圖3 參數kcx的最優估計結果

表1 6發彈丸預報射程與實測射程之差

3.2 基礎阻力系數變化對彈道預報的影響

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

表2 不同基礎阻力系數對應的預報射程與實測射程之差

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

3.3 測量數據對彈道預報的影響

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

表3 利用兩部雷達跟蹤數據預報的射程

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

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

圖4 不同數據長度對應的射程預報誤差

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

3.4 關鍵參數對彈道預報的影響

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

表4 不同Pk,0值對應的預報射程與實測射程之差

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

圖5 不同Pk,0值對應的射程預報誤差極差

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

4 結束語

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