基于改進混雜擴展卡爾曼濾波的炮彈阻力系數(shù)辨識

楊 靖,常思江,王中原

(南京理工大學(xué) 能源與動力工程學(xué)院,江蘇 南京 210094)

基于改進混雜擴展卡爾曼濾波的炮彈阻力系數(shù)辨識

楊 靖,常思江,王中原

(南京理工大學(xué) 能源與動力工程學(xué)院,江蘇 南京 210094)

針對利用炮彈飛行數(shù)據(jù)辨識阻力系數(shù)的問題,提出了一種改進形式的混雜擴展卡爾曼濾波(IHEKF)方法。引入虛擬過程噪聲并且在每個時間步長對稱化協(xié)方差矩陣,避免了混雜擴展卡爾曼濾波在實際應(yīng)用中由于系統(tǒng)建模誤差及彈載計算機的精度有限引起的發(fā)散問題,提高了阻力系數(shù)辨識的魯棒性。數(shù)值仿真結(jié)果表明:當建模誤差可以忽略不計時,采用IHEKF能以高精度辨識出阻力系數(shù);當建模誤差較大時,采用IHEKF具有較好的魯棒性。

炮彈;無控飛行;參數(shù)辨識;擴展卡爾曼濾波

火炮武器系統(tǒng)在戰(zhàn)場中反應(yīng)迅速,能夠提供強大且持續(xù)的火力支援,是未來部隊聯(lián)合作戰(zhàn)的重要組成單元。隨著現(xiàn)代戰(zhàn)爭理念的發(fā)展,遠程壓制、精確打擊及高效毀傷等成為適配彈箭的發(fā)展方向。受炮射平臺的限制,新型適配的彈箭往往修正能力或控制能力有限。因此,提高參數(shù)的準確度,降低模型誤差,對充分發(fā)揮新型彈箭的潛在能力具有重要意義。

氣動參數(shù)的獲得方法分為三大類[1]:計算流體力學(xué)方法,風洞吹風法,利用飛行數(shù)據(jù)進行參數(shù)辨識的方法。數(shù)值計算得到的結(jié)果可作為參數(shù)的先驗知識,風洞試驗與飛行試驗可用來校正參數(shù)。3種方法相輔相成,可以獲得更為準確的氣動參數(shù)。本文研究在已獲得先驗氣動參數(shù)的基礎(chǔ)上,利用飛行試驗數(shù)據(jù)辨識彈箭阻力系數(shù)的問題。

關(guān)于利用飛行試驗數(shù)據(jù)提取彈箭阻力系數(shù)的問題,許多學(xué)者開展了研究。祁載康[2]采用多段樣條曲線擬合的方式,設(shè)計了一種通過雷達遙測數(shù)據(jù)提取彈丸阻力系數(shù)的方法。杜昌平等[3]提出了一種基于遺傳算法的分段多項式擬合阻力系數(shù)辨識方法。Dutta等[4]采用最大似然估計理論,設(shè)計了一種阻力系數(shù)的辨識方法。Rogers等[5]針對智能彈藥,提出了一種基于證據(jù)理論的參數(shù)估計方法。史金光等[6]將阻力系數(shù)辨識問題轉(zhuǎn)化為狀態(tài)估計問題,應(yīng)用擴展卡爾曼濾波方法,通過地面雷達測量的數(shù)據(jù)對彈道修正彈進行了阻力系數(shù)辨識。相比于多項式擬合及最大似然估計等離線辨識方法,將參數(shù)辨識問題轉(zhuǎn)化為狀態(tài)估計問題,可以實現(xiàn)在線遞歸估計。對于某些應(yīng)用問題,實時的參數(shù)辨識具有重要意義。因此,本文將阻力系數(shù)的辨識問題轉(zhuǎn)化為狀態(tài)估計問題。

阻力系數(shù)辨識需要的彈道信息可采用北斗二代接收機獲得。北斗導(dǎo)航系統(tǒng)是我國具有完全自主知識產(chǎn)權(quán)的全球?qū)Ш蕉ㄎ幌到y(tǒng),可以提供三維位置和速度信息,其采用載波相位測量原理的二代產(chǎn)品可實現(xiàn)厘米級的定位精度[7]。

實際中,描述彈箭運動軌跡的彈道模型是連續(xù)時間的,而采用的測量數(shù)據(jù)往往是離散的。采用擴展卡爾曼濾波需要對彈道模型進行離散化處理,而混雜擴展卡爾曼濾波方法(HEKF)[8]可直接應(yīng)用于這種情形。但是,彈道建模與實際的彈箭飛行過程總是存在一定的偏差,且彈載計算機的精度有限,存在截斷誤差,直接應(yīng)用標準形式的HEKF可能出現(xiàn)狀態(tài)估計發(fā)散的問題。

基于上述考慮,本文將利用北斗測量數(shù)據(jù)辨識炮彈阻力系數(shù)的問題轉(zhuǎn)化為狀態(tài)估計問題。引入虛擬過程噪聲并且在每個時間步長對稱化協(xié)方差矩陣,增強混雜擴展卡爾曼濾波方法的魯棒性和狀態(tài)估計精度。應(yīng)用改進的方法進行阻力系數(shù)辨識。

1 問題描述

采用北-天-東地面坐標系作為參考慣性坐標系,記為OxEyEzE。由外彈道理論[9]可知,六自由度剛體彈道模型可以較為全面地描述彈箭飛行過程的運動特征。但是,該模型較為復(fù)雜,系統(tǒng)階數(shù)較高,包含較多參數(shù)以及不可測量的狀態(tài)信息,不便于直接用于參數(shù)辨識。因此,將六自由度彈道模型作為基準彈道,而辨識過程采用下述簡化模型。

為得到合理的簡化模型,采用分層思想,在高度方向?qū)椀肋M行分層,如圖1所示,在每一層內(nèi),炮彈的高度與速度變化不大,從而隨高度變化的空氣密度以及隨彈體質(zhì)心與氣流相對速度變化的阻力系數(shù)可近似為常數(shù)。若炮彈具有較好的動態(tài)穩(wěn)定性,則在每一彈道層內(nèi)其繞心運動對質(zhì)心運動造成的影響相對較小。因此,可以用質(zhì)點彈道模型近似描述在每一彈道層內(nèi)的彈體運動過程。

圖1 彈道分層示意圖

取系統(tǒng)狀態(tài)為

x=(xyzvxvyvz)T

(1)式中:x,y,z分別為彈體質(zhì)心位置分量;vx,vy,vz分別為彈體質(zhì)心速度在地面坐標系OxEyEzE中三軸的分量。

假設(shè)當?shù)仫L速為常值,將隨機風及其他隨機因素對炮彈運動的作用描述為高斯白噪聲,即ω～N(0,Q),可得隨機質(zhì)點彈道模型:

(2)

綜上所述,阻力系數(shù)辨識問題描述為采用式(2)的隨機質(zhì)點彈道模型,利用每一彈道層的飛行數(shù)據(jù)(位置與速度測量值),獲取對應(yīng)的阻力符合系數(shù)kD。

2 改進的混雜擴展卡爾曼濾波器

隨機質(zhì)點彈道模型(2)是連續(xù)時間的,而采用的測量數(shù)據(jù)往往是離散的。混雜卡爾曼濾波方法(HEKF)適用于這種情形的狀態(tài)估計問題。然而,在實際應(yīng)用中,由于彈載計算機的精度有限以及建模誤差等因素,標準的HEKF可能會出現(xiàn)發(fā)散問題。因此,本文在標準HEKF的基礎(chǔ)上引入2種改進策略,得到了一種實用的改進形式的HEKF。

2.1 標準形式的HEKF

考慮狀態(tài)方程由微分方程描述而量測方程是離散的一般形式的動力學(xué)系統(tǒng):

(3)

式中:狀態(tài)向量x∈Rn,輸出向量y∈Rm,向量yk為輸出向量y的第k次測量值,向量xk為對應(yīng)采樣時刻的系統(tǒng)狀態(tài),ω∈Rn為系統(tǒng)過程噪聲,vk∈Rm為量測噪聲,將ω及vk假設(shè)為高斯白噪聲,則

(4)

時間更新方程表示為

(5)

(6)

量測更新方程表示為

(7)

(8)

2.2 改進策略

時間更新方程(5)中的第二式等價于量測協(xié)方差矩陣趨于無窮的微分Riccati方程,要求協(xié)方差矩陣P始終為對稱矩陣。然而,數(shù)值求解該微分方程時,由于計算機精度有限,存在舍入誤差。因此,在時間更新過程中的每一次積分后,對稱化協(xié)方差矩陣P,即令

P=(P+PT)/2

(9)

為后續(xù)表述方便,將采用上述改進策略的標準HEKF記為IHEKF。

3 基于IHEKF的阻力系數(shù)辨識

采用濾波方法進行參數(shù)辨識的思想是將未知常值參數(shù)增廣為系統(tǒng)狀態(tài),利用濾波方法對增廣后的系統(tǒng)進行狀態(tài)估計,從而獲得待辨識的未知常值參數(shù)。本節(jié)以式(2)描述的隨機質(zhì)點模型為基礎(chǔ),將阻力符合系數(shù)kD擴展為系統(tǒng)狀態(tài),得到用于阻力系數(shù)辨識的增廣隨機質(zhì)點模型。然后,利用IHEKF對所得到的增廣模型進行狀態(tài)估計,即實現(xiàn)了阻力系數(shù)的辨識。

3.1 用于阻力系數(shù)辨識的增廣質(zhì)點模型

將阻力系數(shù)增廣為系統(tǒng)狀態(tài),可得:

xa=(xyzvxvyvzkD)T

(10)

結(jié)合式(2),可得新的隨機質(zhì)點彈道模型:

(11)

(12)

對照式(6),可得噪聲增益矩陣為La=I7×7,I為單位矩陣,增廣狀態(tài)偏導(dǎo)數(shù)矩陣為

其中:

(13)

(14)

(15)

3.2 量測模型

北斗導(dǎo)航系統(tǒng)是我國具有完全自主知識產(chǎn)權(quán)的全球?qū)Ш蕉ㄎ幌到y(tǒng),可以提供三維位置和速度信息。因此,本文以北斗定位裝置的測量信息進行阻力系數(shù)的辨識。

北斗導(dǎo)航系統(tǒng)采用東-北-天導(dǎo)航坐標系OxNyNzN。該坐標系與本文采用的地面坐標系OxEyEzE是相對靜止的,其轉(zhuǎn)換易于實現(xiàn)。因此,直接在地面坐標系下描述北斗二代系統(tǒng)的測量模型。

量測向量取為

yk=(xyzvxvyvz)T

(16)

量測方程可表示為

yk=Hkxa,k+Mkvk

(17)

式中:Hk=(I6×6O6×1),Mk=I6×6。

4 數(shù)值仿真分析

為驗證本文提出的改進型混雜卡爾曼濾波器(IHEKF)在實際應(yīng)用中的有效性與優(yōu)越性,本節(jié)分2種情形進行了仿真分析:其一為IHEKF在無系統(tǒng)建模誤差情況下的阻力系數(shù)辨識;其二為存在系統(tǒng)建模誤差情況下的IHEKF與標準形式的HEKF的阻力系數(shù)辨識性能對比。

1)無系統(tǒng)建模誤差(情形1)。

圖2～圖4展示了采用IHEKF在無系統(tǒng)建模誤差情形下的濾波過程。可以看出,位置濾波誤差Δx,Δy,Δz,速度濾波誤差Δvx,Δvy,Δvz與阻力符合系數(shù)辨識誤差ΔkD都很快趨于穩(wěn)定。狀態(tài)估計協(xié)方差矩陣很好地預(yù)測了濾波結(jié)果。表1給出了整個過程的各狀態(tài)估計偏差的標準差S(·)。

圖2 IHEKF方法的位置濾波誤差(情形1)

圖3 IHEKF方法的速度濾波誤差(情形1)

圖4 IHEKF方法的阻力符合系數(shù)辨識誤差(情形1)

表1 估計誤差標準差

2)有系統(tǒng)建模誤差(情形2)。

如前所述,六自由度剛體彈道模型可以較為全面地描述炮彈在空中的實際飛行過程。因此,本節(jié)將六自由度剛體彈道仿真數(shù)據(jù)作為實際的測量數(shù)據(jù)。顯然,式(2)的隨機質(zhì)點彈道模型與六自由度剛體彈道模型相比存在一定的系統(tǒng)建模誤差。

圖5～圖7為HEKF未引入過程噪聲時的濾波過程。

圖5 HEKF方法的位置濾波誤差(情形2)

圖6 HEKF方法的速度濾波誤差(情形2)

圖7 HEKF方法的阻力符合系數(shù)辨識誤差(情形2)

由圖可見,位置及速度濾波誤差出現(xiàn)發(fā)散現(xiàn)象,盡管阻力符合系數(shù)的辨識結(jié)果收斂,但其收斂過程非常緩慢。狀態(tài)估計協(xié)方差矩陣不再具有預(yù)測濾波結(jié)果準確度的性能。

圖8～圖10為引入過程噪聲的IHEKF的濾波過程。

圖9 IHEKF方法的速度濾波誤差(情形2)

圖10 IHEKF方法的阻力符合系數(shù)辨識誤差(情形2)

可以看出,IHEKF方法在存在建模誤差的情況下未出現(xiàn)濾波結(jié)果發(fā)散的現(xiàn)象。狀態(tài)估計協(xié)方差矩陣對濾波結(jié)果具有較好的預(yù)測性。盡管阻力系數(shù)辨識的精度有所下降,但是提高了辨識過程的魯棒性。

5 結(jié)束語

本文針對利用炮彈飛行數(shù)據(jù)辨識阻力系數(shù)的問題,提出了一種改進形式的混雜擴展卡爾曼濾波方法。在濾波過程中,引入虛擬過程噪聲并且在每個時間步長對稱化協(xié)方差矩陣。仿真結(jié)果表明:

①當建模誤差可以忽略不計時,采用改進形式的混雜擴展卡爾曼濾波器能夠以高精度辨識出阻力

系數(shù);

②當建模誤差較大時,采用改進形式的混雜擴展卡爾曼濾波器具有較好的魯棒性。虛擬噪聲過程的調(diào)整需要綜合考慮實際應(yīng)用中的炮彈飛行條件。

[1] KLEIN V,MORELLI E A.Aircraft system identification theory and practice[M].Reston:American Institute of Aeronautics and Astronautics,2006.

[2] QI Z.Multi-spline technique for the extraction of drag coefficients from radar data[J].Journal of Beijing Institute of Technology,1994(1):33-42.

[3] 杜昌平,周德云,宋筆鋒,等.基于遺傳算法的彈道參數(shù)辨識方法研究[J].西北工業(yè)大學(xué)學(xué)報,2008,26(3):373-376. DU Chang-ping,ZHOU De-yun,SONG Bi-feng,et al.Weapon parameter identification algorithm based on genetic algorithm[J].Journal of Northwestern Polytechnical University,2008,26(3):373-376.(in Chinese)

[4] DUTTA G G,SINGHAL A,GHOSH A K.Estimation of drag coefficient from flight data of a cargo shell:AIAA 2006-6149[R].Keystone:AIAA,2006.

[5] ROGERS J,COSTELLO M.Smart projectile state estimation using evidence theory:AIAA 2011-6337[R].Portland:AIAA,2011.

[6] 史金光,劉猛,曹成壯,等.彈道修正彈落點預(yù)報方法研究[J].彈道學(xué)報,2014,26(2):29-33. SHI Jin-guang,LIU Meng,CAO Cheng-zhuang,et al.Research on method of falling point prediction for trajectory correction projectile[J].Journal of Ballistics,2014,26(2):29-33.(in Chinese)

[7] 白征東,元榮,過靜珺,等.北斗二代衛(wèi)星厘米級相對定位[J].測繪通報,2014(7):11-13. BAI Zheng-dong,YUAN Rong,GUO Jing-jun,et al.Performing centimeter-level relative positioning using BeiDou navigation satellite system[J].Bulletin of Surveying and Mapping,2014(7):11-13.(in Chinese)

[8] SIMON D.Optimal state estimation:Kalman,H∞,and nonlinear approaches[M].New Jersey:John Wiley & Sons,Inc.,2006.

[9] 韓子鵬.彈箭外彈道學(xué)[M].北京:北京理工大學(xué)出版社,2008. HAN Zi-peng.Exterior ballistics of projectiles and rockets[M].Beijing:Beijing Institute of Technology Press,2008.(in Chinese)

Drag Coefficient Estimation for Projectiles Based on Improved Hybrid Extended Kalman Filter

YANG Jing,CHANG Si-jiang,WANG Zhong-yuan

(School of Energy and Power Engineering,Nanjing University of Science and Technology,Nanjing 210094,China)

An improved hybrid extended Kalman filter(IHEKF)was proposed aiming at drag coefficient estimation from flight data of projectiles.The fictitious process noise was added to the system model,and the covariance matrix was symmetrized at each time-step.The divergence caused by modeling errors and finite precision of digital microprocessors onboard the projectiles in the implementation of HEKF is avoided,and the robustness is improved.While the modeling errors can be ignored,the IHEKF can estimate the drag-coefficient with high-precision;while the modeling errors are great,IHEKF has great robustness.

projectiles;free flight;parameter estimation;extended Kalman filter

2016-10-11

國家自然科學(xué)基金項目(11402117)

楊靖(1988- ),男,博士研究生,研究方向為彈箭飛行動力學(xué),制導(dǎo)與控制。E-mail:JingYangNUST@163.com。

TJ303.4

A

1004-499X(2017)01-0028-06