丁傳炳

（中國艦船研究設(shè)計中心，上海 201108）

基于雷達(dá)測量數(shù)據(jù)的彈箭擾動源最優(yōu)估計

丁傳炳

（中國艦船研究設(shè)計中心，上海 201108）

為精確計算艦載武器的飛行狀態(tài)參數(shù)，以彈體縱向運(yùn)動過程為研究對象，推導(dǎo)了包含誤差干擾源在內(nèi)的縱向擾動運(yùn)動學(xué)方程，利用“系數(shù)凍結(jié)法”及拉普拉斯變換得到解析解，擬合出彈體被動段縱向運(yùn)動的氣動參數(shù)公式，采用三坐標(biāo)雷達(dá)測量量作為系統(tǒng)量測方程，從而對氣動參數(shù)誤差干擾源進(jìn)行最優(yōu)估計。計算結(jié)果表明：該算法可以使俯仰操縱力矩系數(shù)導(dǎo)數(shù)誤差的精度穩(wěn)定在±0.007 (°)-1范圍內(nèi)，阻力系數(shù)誤差的精度趨于±0.025之間；升力系數(shù)誤差基本穩(wěn)定在±0.12之間；俯仰力矩系數(shù)對攻角的導(dǎo)數(shù)偏差精度穩(wěn)定在±0.011 (°)-1之間；俯仰阻尼力矩系數(shù)誤差的精度收斂于±0.009 (°)-1之間，且算法收斂速度快，可為重構(gòu)高精度的彈道參數(shù)提供支撐。

最優(yōu)估計；縱向運(yùn)動；擾動源；雷達(dá)

針對彈箭飛行狀態(tài)參數(shù)最優(yōu)估計問題，目前大多采取的方法是通過各種外彈道試驗測量數(shù)據(jù)進(jìn)行彈道濾波[1-3]。為了從源頭上提高外彈道飛行狀態(tài)參數(shù)的估計精度，本文提出采用雷達(dá)觀測量來對氣動系數(shù)誤差干擾源進(jìn)行最優(yōu)估計，進(jìn)而為重構(gòu)高精度的彈箭飛行狀態(tài)參數(shù)創(chuàng)造條件?？紤]到有控飛行器的空間運(yùn)動可以看成是鉛直面內(nèi)的運(yùn)動與側(cè)向平面內(nèi)運(yùn)動的合成，并且在許多情況下，主要是在一個鉛直平面內(nèi)的飛行，因此本文以彈箭縱向運(yùn)動面為研究對象，運(yùn)用三坐標(biāo)雷達(dá)觀測量對縱向氣動參數(shù)誤差干擾源進(jìn)行最優(yōu)估計，估計出來的高精度氣動參數(shù)可為重構(gòu)飛行外彈道最優(yōu)狀態(tài)參數(shù)提供數(shù)據(jù)支持。

1 縱向擾動運(yùn)動方程

如果飛行器的氣動對稱面與鉛直面重合，并且質(zhì)心也在鉛直面內(nèi)運(yùn)動，則稱為縱向運(yùn)動。對于制導(dǎo)控制系統(tǒng)工作正常的彈箭，其實際飛行運(yùn)動參數(shù)也總在理想彈道運(yùn)動參數(shù)附近變化，即彈箭受到控制或干擾產(chǎn)生的擾動可認(rèn)為是加在理想運(yùn)動上的小擾動，因此就可將實際運(yùn)動參數(shù)看作是理想彈道運(yùn)動參數(shù)與對應(yīng)偏差量之和，因偏差量很小，故可將縱向運(yùn)動方程在基準(zhǔn)運(yùn)動附近進(jìn)行線性化，組成關(guān)于偏差量的線性運(yùn)動方程[4-6]，即：

式中：Ρ表示發(fā)動機(jī)推力，V表示彈丸飛行速度，X表示阻力，Y表示升力，1δ表示攻角，δz表示升降舵偏角，m表示彈丸質(zhì)量[6-8]。由縱向擾動運(yùn)動方程可以看出縱向運(yùn)動方程的擾動變量主要有升力系數(shù)誤差阻力系數(shù)誤差ΔCx、俯仰力矩系數(shù)導(dǎo)數(shù)偏差俯仰操縱力矩系數(shù)導(dǎo)數(shù)誤差俯仰阻尼力矩系數(shù)誤差五項。

為了方便采用卡爾曼濾波進(jìn)行彈道重構(gòu)，可將動力學(xué)形式的小擾動方程轉(zhuǎn)化為矩陣的形式，可以得到[8-9]：

式(2)中縱向動力系數(shù)系數(shù)采用aij的符號表示，的第一個腳注i表示運(yùn)動方程的序號，第二個腳注j表示運(yùn)動參數(shù)偏量的順序號[6]。由于本文研究的是擾動因素對彈道的影響，即彈箭對各固有飛行參數(shù)誤差的反應(yīng)，而不糾纏于飛行器對舵面偏轉(zhuǎn)的影響，因此可以認(rèn)為

式(3)是變系數(shù)線性微分方程，采用“系數(shù)凍結(jié)法”和拉氏變換可以解得[10]：因此式(2)可表示為[5-6]：

將式(4)代入式(3)得：

2 擾動源的狀態(tài)建模

2.1 擾動源系統(tǒng)狀態(tài)模型

假設(shè)擾動源服從下列方程[4-5]：

式中，F(xiàn)M為5× 5的矩陣，WM假設(shè)為零均值的隨機(jī)白噪聲。

由于氣動系數(shù)RM是飛行狀態(tài)XM的函數(shù)（如圖1～5），它們之間的具體表達(dá)式可由實驗數(shù)據(jù)擬合而來或由經(jīng)驗公式簡化而來[6]。

圖1 阻力系數(shù)隨攻角和馬赫數(shù)的變化曲線Fig.1 Variation of resistance coefficient with angle of attack and Mach number

圖2 隨馬赫數(shù)變化曲線Fig.2 Variation of with Mach number

圖3 隨馬赫數(shù)變化曲線Fig.3 Variation ofwith Mach number

圖4 隨馬赫數(shù)變化曲線Fig.4 Variation of with Mach number

圖5 隨馬赫數(shù)變化曲線Fig.5 Variation of with Mach number

本文根據(jù)一組給定的某彈箭氣動參數(shù)數(shù)據(jù)，用數(shù)值擬合的辦法得到一組被動段縱向氣動參數(shù)的計算公式，然后將擬合出來的公式進(jìn)行泰勒級數(shù)展開，截取一階項得到小擾動方程，寫成矩陣的形式，如下[6-7]：

將式(4)(5)代入式(7)，可得：

根據(jù)式(8)可以得到矩陣FM中的各元素，式(8)即為系統(tǒng)的狀態(tài)模型。

由于采用卡爾曼濾波算法時測量值必須是系統(tǒng)狀態(tài)參數(shù)的線性組合，因此考慮利用ΔXM與ΔR之間的關(guān)系，尋找δρ與ΔR之間的線性關(guān)系，將Δx、Δy、Δz擴(kuò)充到縱向擾動方程(3)中的狀態(tài)量中[9-11]，則：

求解上述方程(15)，得：

2.2 系統(tǒng)量測模型

三坐標(biāo)雷達(dá)可以實時測定彈丸飛行時的空間三維坐標(biāo)（見圖6示意圖）。

圖6 雷達(dá)測量飛行彈丸示意圖Fig.6 Sketch map of radar measuring flight projectile

設(shè)彈丸的真實位置坐標(biāo)為(x,y,z)，運(yùn)動方程計算出的彈丸位置為(xE,yE,zE)，三坐標(biāo)雷達(dá)測量的坐標(biāo)為(xL,yL,zL)，彈丸至發(fā)射原點的真實距離為r，運(yùn)動方程計算的距離為Eρ，則有[12]：

相對真實位置(x,y,z)，將Eρ進(jìn)行一階泰勒展開，得：

雷達(dá)測量誤差為雷達(dá)測量值與彈丸真實位置的差值，即

雷達(dá)測量的彈丸至發(fā)射點的距離為

則雷達(dá)與運(yùn)動方程計算得到的距離之差為：

式(22)即為系統(tǒng)的量測模型。

3 擾動源的最優(yōu)估計算法

由第2節(jié)公式(9)和(22)可知，擾動源系統(tǒng)的狀態(tài)方程和量測方程具有如下形式：

式中，W(t)、V(t)分別稱為連續(xù)系統(tǒng)的系統(tǒng)噪聲矩陣和量測噪聲矩陣。將上述方程(23)進(jìn)行離散化，可得：

根據(jù)離散后的方程，可以設(shè)計卡爾曼濾波器，對相關(guān)狀態(tài)變量進(jìn)行最優(yōu)估計[13]。

4 仿真計算

為進(jìn)行動態(tài)仿真，需要給定火箭彈的飛行軌跡（見圖7）。仿真初始條件如下：

圖7 被動段縱向彈道平面Fig.7 Longitudinal trajectory plane of passive segment

采樣周期為0.05 s。仿真結(jié)果如圖7～12所示。

由仿真圖 7～12可以看出，系統(tǒng)各狀態(tài)的估計值精度較高，基本上在濾波開始后的20s以后，其狀態(tài)估值均趨于穩(wěn)定，俯仰操縱力矩系數(shù)導(dǎo)數(shù)誤差的精度穩(wěn)定在±0.007 (°)-1范圍內(nèi)，阻力系數(shù)誤差ΔCx的精度趨于±0.025之間，升力系數(shù)誤差ΔCy基本穩(wěn)定在±0.12之間，俯仰力矩系數(shù)對1δ的導(dǎo)數(shù)偏差精度穩(wěn)定在±0.011 (°)-1之間，俯仰阻尼力矩系數(shù)誤差的精度基本上收斂于±0.009 (°)-1之間。由此可以看出：采用雷達(dá)測量數(shù)據(jù)對彈箭氣動系數(shù)誤差干擾源進(jìn)行最優(yōu)估計，可以得到精度較高的氣動參數(shù)，進(jìn)而為重構(gòu)高精度的飛行外彈道狀態(tài)參數(shù)創(chuàng)造條件。

圖8 的狀態(tài)估值Fig.8 The state valuation of

圖9 ΔCx的狀態(tài)估值Fig.9 The state valuation ofΔCx

圖10 ΔCy的狀態(tài)估值Fig.10 The state valuation ofΔCy

圖11 的狀態(tài)估值Fig.11 The state valuation of

圖12 的狀態(tài)估值Fig.12 The state valuation of

5 結(jié) 論

本文以彈體的縱向運(yùn)動面為研究對象，推導(dǎo)了包含干擾源在內(nèi)的縱向擾動方程，建立了以氣動系數(shù)誤差干擾源為狀態(tài)量的系統(tǒng)狀態(tài)模型和以雷達(dá)觀測量為量測量的系統(tǒng)量測模型，并對模型進(jìn)行仿真計算。計算結(jié)果表明：利用雷達(dá)觀測量對被動段縱向平面的氣動系數(shù)誤差干擾源進(jìn)行最優(yōu)估計是有效的。該算法對提高彈箭的導(dǎo)航控制精度具有較大的現(xiàn)實意義，可為工程應(yīng)用提供一定的技術(shù)參考。

（References）：

[1]Li Yong-chen, Li Jian-xun. Robust adaptive kalman filtering for target tracking with unknown observation noise[C]//Proc of 24th Chinese Control and Decision Conference. Taiyuan, 2012: 2075-2080.

[2]楊榮軍, 王良明, 修觀, 等. 利用雷達(dá)測量數(shù)據(jù)的實際彈道重建[J]. 彈道學(xué)報, 2011, 23(3): 43-46.Yang Rong-jun, Wang Liang-ming, Xiu Guan, et al.Trajectory reconstruction using radar measured data[J].Journal of Ballistics, 2011, 23(3): 43-46.

[3]Park D B, Shin D H, Oh S H. Development of a GPS/INS system for precision GPS guided bombs[J]. IEEE Aerospace and Electronic Systems Magazine, 2012, 27(3):31-39.

[4]Kim J, Vaddi S S, Menon P K. Comparison between nonlinear filtering techniques for spiraling ballistic missile state estimation[J]. IEEE Transactions on Aerospace and Electronic Systems, 2012, 48(1): 313-328.

[5]Gupta K K, Voelker L S. Aeroelastic simulation of hypersonic flight vehicles[J]. AIAA Journal, 2012, 50(3):717-723.

[6]韓子鵬. 彈箭外彈道學(xué)[M]. 北京: 北京理工大學(xué)出版社, 2008.

[7]Danowsky B P, Chrstos J R, Klyde D H, et al. Evaluation of aeroelastic uncertainty analysis methods[J].Journal of Aircraft, 2010, 47(4): 1266-1273.

[8]Babbar Y, Suryakumar V S, Strganac T W. Experiments in free and forced aeroelastic response[C]//51st AIAA Aerospace Sciences Meeting. Texas, 2013.

[9]Brown R L, Das K, Whitcomb J D, et al. Aeroelastic simulation of structures in hypersonic flow[C]//Proceedings of 53rd AIAA/ASME/ASCE/AHS/ASC Structures,Structural Dynamics, and Materials Conference. Honolulu,Hawaii, USA, 2012.

[10]夏紅偉, 李秋實, 李莉, 等. 基于 hp 自適應(yīng)偽譜法的飛行器再入軌跡優(yōu)化與制導(dǎo)[J]. 中國慣性技術(shù)學(xué)報,2015, 23(6): 818-823.Xia Hong-wei, Li Qiu-shi, Li Li, et al. Trajectory optimazation and guidance for reentry craft based on hp-adaptive pseudospectral method[J]. Journal of Chinese Inertial Technology, 2015, 23(6): 818-823.

[11]梁明. 靶彈動力學(xué)模型參數(shù)辨識與彈道仿真[J]. 戰(zhàn)術(shù)導(dǎo)彈技術(shù), 2012(3): 48-53.Liang Ming. Dynamic model parameter identification and trajectory simulation of target missile[J]. Tactical Missile Technology, 2012(3): 48-53.

[12]史金光, 劉猛, 曹成壯, 等. 艦船運(yùn)動對彈道修正彈預(yù)報彈道落點的影響分析[J]. 南京理工大學(xué)學(xué)報, 2014,38(3): 366-374.Shi Jin-guang, Liu Meng, Cao Cheng-zhuang, et al.Effects of warship motion on predicted trajectory falling points for trajectory correction projectiles[J]. Journal of Nanjing University of Science and Technology, 2014,38(3): 366-374.

[13]周宏宇, 王小剛, 崔乃剛, 等. 基于 hp 自適應(yīng)偽譜法的組合動力可重復(fù)使用運(yùn)載器軌跡優(yōu)化[J]. 中國慣性技術(shù)學(xué)報, 2016, 24(6): 832-837.Zhou Hong-yu, Wang Xiao-gang, Cui Nai-gang, et al.Trajectory optimization of reusable vehicle with combined power based on hp adaptive pseudospectral algorithm[J]. Journal of Chinese Inertial Technology,2016, 24(6): 832-837.

Optimal estimation of missile disturbance source based on radar measurement data

DING Chuan-bing
(China Ship Development and Design Center, Shanghai 201108, China)

In order to calculate the flight state parameters of shipborne weapons accurately, the equations of longitudinal perturbation kinematics with error disturbance source are deduced based on the longitudinal motion process of the missile. The “coefficient freezing method” and the Laplace transform are used to analyze the motion equation solution. The aerodynamic parameter formula of the longitudinal motion during passive phase is fitted out, and the measurements of the coordinate radar are used as the system measurement equation to optimally estimate the aerodynamic parameter error source. Experiment results show that: the algorithm can stabilize the error of the steering torque coefficient to within ±0.007 (°)-1, and the accuracy of the resistance coefficient tends to be within ±0.025; the lift coefficient error is basically stable to within ±0.12; the accuracy of the torque coefficient is within ±0.011 (°)-1; the accuracy of the torque coefficient of the pitching damping converges to within ±0.009 (°)-1, and the convergence speed of the algorithm is fast. The proposed algorithm can provide technical references for the reconstruction of high-precision trajectory parameters.

optimal estimation; longitudinal motion; perturbation source; radar

1005-6734(2017)04-0544-06

10.13695/j.cnki.12-1222/o3.2017.04.021

TJ413.6

A

2017-05-20；

2017-07-24

兵科院重點預(yù)研項目（20402020101）

丁傳炳（1984—），男，工程師，博士，從事水面艦船作戰(zhàn)系統(tǒng)總體技術(shù)研究。E-mail: sfdcb802@163.com