基于氣動(dòng)參數(shù)辨識(shí)的彈道風(fēng)速估計(jì)算法設(shè)計(jì)

張衍儒 肖練剛 張繼生 田 豐 陳 昌

北京航天自動(dòng)控制研究所，北京 100854

基于氣動(dòng)參數(shù)辨識(shí)的彈道風(fēng)速估計(jì)算法設(shè)計(jì)

張衍儒 肖練剛 張繼生 田 豐 陳 昌

北京航天自動(dòng)控制研究所，北京 100854

首先介紹低空風(fēng)分布，確定風(fēng)速估算范圍。然后，以射程預(yù)測(cè)為主要研究目標(biāo)，通過(guò)彈上傳感器信息分析飛行過(guò)程中的實(shí)時(shí)氣動(dòng)阻力，對(duì)比在無(wú)風(fēng)條件下的氣動(dòng)阻力理論值，采用氣動(dòng)參數(shù)辨識(shí)算法實(shí)時(shí)估算縱向風(fēng)的近似風(fēng)速，從而提高射程預(yù)測(cè)的精度。估算過(guò)程中，采用最小二乘算法計(jì)算不同風(fēng)速下的氣動(dòng)阻力系數(shù)常值項(xiàng)，利用迭代算法求取氣動(dòng)阻力系數(shù)常值項(xiàng)為理論值時(shí)所對(duì)應(yīng)的設(shè)定縱向風(fēng)，即實(shí)際縱向風(fēng)的近似值。最后通過(guò)仿真試驗(yàn)給出了估算算法的分析結(jié)果，該算法可以實(shí)時(shí)估算風(fēng)速用于射程預(yù)測(cè)，并具有運(yùn)算速度快，計(jì)算負(fù)荷小等優(yōu)點(diǎn)，適合于工程實(shí)際應(yīng)用。

制導(dǎo)彈藥;氣動(dòng)參數(shù)辨識(shí);最小二乘算法;風(fēng)速估計(jì)

目前，部分制導(dǎo)彈藥采用預(yù)測(cè)落點(diǎn)的方式進(jìn)行制導(dǎo)控制，如美國(guó)XM1156式制導(dǎo)炮彈［1］、RCGM式制導(dǎo)迫擊炮彈［2］等。預(yù)測(cè)落點(diǎn)過(guò)程中，需要分析彈道風(fēng)速，如美國(guó)XM1156式制導(dǎo)炮彈利用發(fā)射前半小時(shí)的氣象分析數(shù)據(jù)，進(jìn)行落點(diǎn)預(yù)測(cè)。90年代，Mark F.Costello對(duì)有風(fēng)條件下，鴨舵式制導(dǎo)迫擊炮彈道仿真模型進(jìn)行了詳細(xì)的研究［3］，從而在風(fēng)速已知的條件下提高了預(yù)測(cè)落點(diǎn)的精度。2008年，趙慶嵐［4］對(duì)有風(fēng)條件下末制導(dǎo)炮彈彈道模型進(jìn)行了研究，同樣提高了預(yù)測(cè)落點(diǎn)的精度。但是，上述方法需要在射前進(jìn)行風(fēng)速測(cè)定，并且重新生成射表，用于射前數(shù)據(jù)裝訂。2011年Frank Fresconi等人［5］對(duì)彈上實(shí)時(shí)落點(diǎn)預(yù)測(cè)算法進(jìn)行了理論研究，利用彈上GPS提供的位置、速度信息積分求取預(yù)測(cè)落點(diǎn)。該方法由于利用彈上實(shí)時(shí)數(shù)據(jù)進(jìn)行落點(diǎn)預(yù)測(cè)，落點(diǎn)精度有顯著提高，但是需要進(jìn)行大量的積分運(yùn)算，計(jì)算負(fù)荷大，對(duì)于全程實(shí)時(shí)控制的制導(dǎo)彈藥，此種方式并不適用。為保證預(yù)測(cè)的實(shí)時(shí)性，可以采用Jonathan Rogers等人［7］提供的完整的無(wú)風(fēng)仿真模型，利用在標(biāo)準(zhǔn)氣象條件下的氣動(dòng)參考值，建立基于導(dǎo)航系下彈體三維速度、三維位置的落點(diǎn)預(yù)測(cè)插值表。落點(diǎn)預(yù)測(cè)插值表在發(fā)射前裝訂到彈丸中，彈丸在飛行過(guò)程中實(shí)時(shí)查詢預(yù)測(cè)表，采用插值方式進(jìn)行落點(diǎn)預(yù)測(cè)。此方式可以大幅降低處理器計(jì)算負(fù)荷，但由于未考慮風(fēng)速影響，降低了預(yù)測(cè)精度。

本文在未知風(fēng)速的條件下，利用彈上加速度計(jì)、GPS和彈體姿態(tài)角信息實(shí)時(shí)分析飛行過(guò)程中的氣動(dòng)阻力，對(duì)比在無(wú)風(fēng)條件下的氣動(dòng)阻力理論值，采用氣動(dòng)參數(shù)辨識(shí)算法［8-9］實(shí)時(shí)估算縱向風(fēng)的近似值，并利用縱向風(fēng)對(duì)射程的修正模型，完成對(duì)射程預(yù)測(cè)值的修正，從而提高射程預(yù)測(cè)的精度。

1 風(fēng)的分布和縱向風(fēng)對(duì)射程的修正

習(xí)慣上將風(fēng)力大小用風(fēng)力等級(jí)來(lái)描述［10］，如0～10級(jí)風(fēng)分別稱為無(wú)風(fēng)、軟風(fēng)、輕風(fēng)、微風(fēng)、和風(fēng)、清勁風(fēng)、強(qiáng)風(fēng)、疾風(fēng)、大風(fēng)、烈風(fēng)、狂風(fēng)。風(fēng)力等級(jí)K與平均風(fēng)速W可用下式換算:

由式(1)可知，低空風(fēng)速大小范圍在26.4±2m/s以內(nèi)。因此文中考慮的風(fēng)速設(shè)定值范圍為-40 ～40m/s。

可以根據(jù)目前彈體速度v0、射角θ0和高度H，采用插值的方式進(jìn)行彈道的射程預(yù)測(cè)。在有縱風(fēng)wx時(shí)，可由式(2)求取動(dòng)坐標(biāo)系中相對(duì)速度vr與相對(duì)射角θr:

根據(jù)所求相對(duì)速度vr、相對(duì)射角θr和高度H，采用插值的方式求取相對(duì)射程Xr與飛行時(shí)間T，由式(3)求取修正后預(yù)測(cè)射程Xwx:

由式(3)可以看出，修正模型的縱風(fēng)wx是指彈道的平均風(fēng)速，因此估算縱風(fēng)風(fēng)速時(shí)，應(yīng)濾掉短時(shí)隨機(jī)風(fēng)，保證射程修正模型的有效性。

圖1 縱風(fēng)對(duì)射程的修正曲線

圖1介紹了修正過(guò)程包含的3條曲線:“標(biāo)準(zhǔn)”曲線表示未修正時(shí)的預(yù)測(cè)彈道;“相對(duì)”曲線表示在動(dòng)坐標(biāo)系中求取的相對(duì)彈道;“有風(fēng)”曲線表示修正后的預(yù)測(cè)彈道。圖中預(yù)測(cè)彈道參數(shù)包含彈體速度v0，射角θ0和高度H，相對(duì)彈道參數(shù)包含相對(duì)速度vr，相對(duì)射角θr和高度H，修正后的預(yù)測(cè)彈道參數(shù)包含相對(duì)射程Xr，飛行時(shí)間T和縱風(fēng)wx。

2 氣動(dòng)參數(shù)辨識(shí)算法

制導(dǎo)彈藥氣動(dòng)阻力的三維分量可用式(4)表示:

式中:FG_DX，F(xiàn)G_DY，F(xiàn)G_DZ是制導(dǎo)彈藥氣動(dòng)阻力的三維分量，m，c，H，y，G，vr，cs是彈體質(zhì)量、彈形系數(shù)、空氣密度函數(shù)、高度、阻力函數(shù)、相對(duì)速度、聲速;vx，vy，vz是導(dǎo)航系下彈藥速度;wx，wz是導(dǎo)航系下縱向風(fēng)速度、橫向風(fēng)速度;d，ρ是彈體特征長(zhǎng)度、空氣密度;CD是氣動(dòng)阻力系數(shù)。

通過(guò)式(6)，求取攻角α、側(cè)滑角β、彈體速度與彈體軸向之間的夾角αT:

利用加速度計(jì)測(cè)量值aximu，ayimu，azimu和彈體姿態(tài)信息α，β，通過(guò)式(7)可以實(shí)時(shí)計(jì)算與彈丸速度反向的制導(dǎo)彈藥氣動(dòng)阻力FACC_D。

式中，F(xiàn)A_Xb，F(xiàn)A_Yb，F(xiàn)A_Zb表示加速度計(jì)測(cè)得的彈體三維氣動(dòng)力;FACC_D表示通過(guò)加速度計(jì)計(jì)算得到的制導(dǎo)彈藥氣動(dòng)阻力。

通過(guò)式(6)和(7)，可以得到氣動(dòng)阻力FACC_D關(guān)于彈體速度vXb，vYb，vZb和加速度aximu，ayimu，azimu的表達(dá)式:

計(jì)算得到的氣動(dòng)阻力FACC_D與氣動(dòng)阻力三維分量FG_DX，F(xiàn)G_DY，F(xiàn)G_DZ存在關(guān)系式(9):

根據(jù)式(8)和(10)可以得到式(11)，用于計(jì)算氣動(dòng)阻力系數(shù)CD:

氣動(dòng)阻力系數(shù)CD，主要受到彈體馬赫數(shù)、彈體軸向與彈體速度之間夾角αT、氣溫、氣壓、空氣濕度、舵角等項(xiàng)的影響。其中αT的變化對(duì)CD數(shù)值變化影響較大，因此用式(12)近似αT對(duì)CD的影響:

式中，CD0為氣動(dòng)阻力系數(shù)常值項(xiàng)，CDαT，CDαT2表示氣動(dòng)阻力系數(shù)關(guān)于αT的1，2階變化量。

利用計(jì)算得到的氣動(dòng)阻力系數(shù)CD求取氣動(dòng)阻力系數(shù)的常值項(xiàng)CD0，求取過(guò)程如下:

應(yīng)用最小二乘估計(jì)法，得到C的估計(jì)值為

取極小值。

由求取的C的第一行值可以確定氣動(dòng)阻力系數(shù)常值項(xiàng)CD0。

3 風(fēng)速迭代算法

對(duì)于全程控制的制導(dǎo)彈藥，橫向控制可以全程采用比例導(dǎo)引進(jìn)行精確制導(dǎo)，因此橫向風(fēng)的影響可以通過(guò)舵控進(jìn)行實(shí)時(shí)修正。但是制導(dǎo)彈藥在上升段飛行時(shí)縱向控制算法需要實(shí)時(shí)預(yù)測(cè)射程用于精確制導(dǎo)控制，考慮到縱向風(fēng)對(duì)預(yù)測(cè)射程的影響，本文設(shè)計(jì)了風(fēng)速迭代算法實(shí)時(shí)求取縱向風(fēng)的近似值，用于預(yù)測(cè)射程的修正。

通過(guò)第2節(jié)對(duì)氣動(dòng)參數(shù)辨識(shí)算法的介紹，可知?dú)鈩?dòng)阻力相同時(shí)，設(shè)定不同風(fēng)速w'x(k)會(huì)對(duì)應(yīng)不同的氣動(dòng)阻力系數(shù)常值項(xiàng)C'D0(k)，由于氣動(dòng)阻力系數(shù)常值項(xiàng)CD0僅受到彈體馬赫數(shù)、氣溫、氣壓、空氣濕度和舵角等項(xiàng)的影響，因此可以通過(guò)裝訂在彈上由風(fēng)洞試驗(yàn)提供的相關(guān)插值表，快速計(jì)算其實(shí)際的參考值CDref，只有當(dāng)C'D0(k0)≈CDref時(shí)，設(shè)定的w'x(k0)才是實(shí)際縱向風(fēng)的近似值。下面詳細(xì)介紹一下迭代算法的工作原理:

1)根據(jù)風(fēng)的分布，本文設(shè)定縱向風(fēng)速的迭代范圍為-40～40m/s，由于橫向控制算法修正了橫向風(fēng)的影響，因此設(shè)定縱風(fēng)與橫風(fēng)為:

利用導(dǎo)航初始化過(guò)程中的前10組共1s的數(shù)據(jù)，通過(guò)設(shè)定不同風(fēng)速w'x(k)采用氣動(dòng)參數(shù)辨識(shí)算法，求取不同設(shè)定風(fēng)速所對(duì)應(yīng)的氣動(dòng)阻力系數(shù)常值項(xiàng)C'D0(k)，k=1～17。

2)初步確定風(fēng)速w'x(k0)后，縮小縱風(fēng)設(shè)定范圍進(jìn)行精確求取，設(shè)定風(fēng)速為:

式中，t為采樣次數(shù)，單位為s。

每隔1s時(shí)間，求取一次設(shè)定風(fēng)速(t，j)所對(duì)應(yīng)氣動(dòng)阻力系數(shù)常值項(xiàng)(t，j)和實(shí)際參考值CDref(t)之間的差值絕對(duì)值。在介紹縱向風(fēng)對(duì)射程的修正模型時(shí)，已經(jīng)說(shuō)明修正模型中所采用的縱風(fēng)wx是指彈道的平均風(fēng)速，因此估算縱風(fēng)風(fēng)速時(shí)，應(yīng)濾掉短時(shí)隨機(jī)風(fēng)，保證射程修正模型的有效性。為消除零均值噪聲和短時(shí)隨機(jī)風(fēng)的影響，采用求和的方式記錄差值絕對(duì)值的累加值

3)根據(jù)求取的S(t，j)(j=1，2，4，5)變化和初步確定的風(fēng)速w'x(k0)，利用式(20)實(shí)時(shí)求解采樣時(shí)刻為t時(shí)，縱向風(fēng)的近似值wx(t):

根據(jù)式(18)所得縱風(fēng)的近似值wx(t)，可以采用式(2)與(3)進(jìn)行預(yù)測(cè)射程的實(shí)時(shí)修正。

4 仿真試驗(yàn)與誤差分析

通過(guò)軟件流程圖再次介紹利用迭代算法求取縱風(fēng)近似值的過(guò)程，如圖2。

圖2 風(fēng)速估計(jì)算法軟件流程圖

在仿真試驗(yàn)中，設(shè)定實(shí)際縱向風(fēng)速wx=20，設(shè)定不同的風(fēng)速w'x(k)=0，10，20，40 ，通過(guò)最小二乘算法計(jì)算不同設(shè)定風(fēng)速所對(duì)應(yīng)的氣動(dòng)阻力系數(shù)常值項(xiàng)C'D0(k)，如圖3所示。

圖3 最小二乘算法求取氣動(dòng)阻力常值項(xiàng)

根據(jù)圖3可知，不同設(shè)定風(fēng)速w'x(k)代入會(huì)對(duì)應(yīng)不同的氣動(dòng)阻力系數(shù)常值項(xiàng)C'D0(k)。當(dāng)設(shè)定風(fēng)速w'x(k)越接近實(shí)際風(fēng)速wx時(shí)，C'D0(k)與由風(fēng)洞試驗(yàn)提供的實(shí)際參考值CDref之間差值絕對(duì)值的累加項(xiàng)S(t)越小，如圖4所示。

圖4 不同風(fēng)速設(shè)定的S(t)值變化

圖4顯示了不同風(fēng)速設(shè)定對(duì)應(yīng)的S(t)值變化，可以看出設(shè)定風(fēng)速越接近于實(shí)際風(fēng)速，S(t)值越小。

采用風(fēng)速估算算法，進(jìn)行4組仿真試驗(yàn)，設(shè)定實(shí)際縱向風(fēng)速wx=-14，0，14，20 ，求取近似風(fēng)速，求取過(guò)程中的迭代曲線如圖5所示。

圖5 利用迭代算法求取不同實(shí)際風(fēng)速曲線

圖5中，X向表示縱向設(shè)定風(fēng)速，Y向表示差值絕對(duì)值的累加項(xiàng)S(t)的變化，可以看出設(shè)定風(fēng)速?gòu)膶?shí)際風(fēng)速的兩側(cè)逐漸收斂，最后收斂到最低累加值，此位置的設(shè)定風(fēng)速近似為實(shí)際風(fēng)速，從而證明了風(fēng)速估算算法的快速性與準(zhǔn)確性。近似風(fēng)速與實(shí)際設(shè)定風(fēng)速之間存在一定誤差，誤差來(lái)源于舵角與橫向風(fēng)的變化對(duì)風(fēng)速求取過(guò)程的影響。

采用仿真試驗(yàn)，對(duì)比采用風(fēng)速估算算法修正前后的預(yù)測(cè)射程與實(shí)際射程之間偏差變化曲線，如圖6所示。圖中X向表示制導(dǎo)彈藥飛行時(shí)間，Y向表示預(yù)測(cè)射程與實(shí)際射程之間的偏差，通過(guò)對(duì)比可以看出采用修正算法后，飛行控制時(shí)間大幅降低，在20s時(shí)，修正過(guò)程就已經(jīng)完成，而未采用修正算法，由于預(yù)測(cè)射程的誤差需要長(zhǎng)時(shí)間的比例導(dǎo)引進(jìn)行制導(dǎo)修正，這對(duì)彈體結(jié)構(gòu)、強(qiáng)度要求、舵機(jī)功率要求、彈上儀器設(shè)備和減少制導(dǎo)誤差都是不利的。

圖6 修正前后預(yù)測(cè)射程與實(shí)際射程之間偏差對(duì)比圖

5 結(jié)論

在算法實(shí)際使用過(guò)程中氣溫、氣壓、空氣濕度、舵角、橫風(fēng)風(fēng)速等參數(shù)變化都會(huì)對(duì)算法產(chǎn)生影響，因此求取的風(fēng)速僅僅為縱向風(fēng)的近似值。雖然采用近似風(fēng)速應(yīng)用在射程修正時(shí)存在一定的預(yù)測(cè)偏差，但是相對(duì)于不考慮風(fēng)速的預(yù)測(cè)方式卻大幅提高了預(yù)測(cè)精度，同時(shí)相對(duì)于積分方式的預(yù)測(cè)算法，也降低了處理器計(jì)算負(fù)荷，增強(qiáng)了算法的實(shí)時(shí)性與可應(yīng)用性，為制導(dǎo)彈藥提高預(yù)測(cè)射程精度提供了可行性方案。

［1］Peter Burke.XM1156 Precision Guidance Kit(PGK)Overview［C］.NDIA-53rd Annual Fuze Conference，2009.

［2］Yousef Habash.Roll Control Guided Mortar(RCGM)［C］.America:NDIA Joint Armaments Conference，2012.

［3］Mark F.Costello.Range Extension and Accuracy Improvement of an Advanced Projectile Using Canard Control［C］.AIAA Atmospheric Flight Mechanics Conference in Baltimore，Maryland，1995:324-331.

［4］趙慶嵐，宋衛(wèi)東，魯飛，喬梁.有風(fēng)條件下末制導(dǎo)炮彈彈道模型研究［J］.工程設(shè)計(jì)學(xué)報(bào)，2008，5(15):373-376.(Zhao fenglan，Song Weidong，Lu Fei，Qiao Liang.Research on Terminal-guidance Ballistic Model under Windy Condition［J］.Journal of Engineering Design，2008，5(15):373-376.)

［5］Frank Fresconi，Gene Cooper，Mark Costello.Practical Assessment of Real-Time Impact Point Estimators for Smart Weapons［J］.Journal of Aerospace Engineering，2011，1(24):1-11.

［6］王敏忠.炮兵應(yīng)用外彈道學(xué)及仿真［M］.北京:國(guó)防工業(yè)出版社，2009:25-28.(Wang Minzhong.Artillery Applied Exterior Ballistics and Simulation［M］.Beijing:National Defence Industry Press，2009:25-28)

［7］Jonathan Rogers，Mark Costello.Design of a Roll-Stabilized Mortar Projectile with Reciprocating Canards［J］.Journal of Guidance，Control，and Dynamics，2010，4(33):1026-1034.

［8］戴明祥，楊新民，易文俊，何穎.用于衛(wèi)星制導(dǎo)彈藥落點(diǎn)預(yù)測(cè)的卡爾曼濾波算法［J］.彈箭與制導(dǎo)學(xué)報(bào)，2012，5(32):117-120.(Dai Mingxiang，Yang Xinmin，Yi Wenjun，He Ying.Kalman Filtering Algorithm for Impact Point Prediction of Satellite Guided Projectile［J］.Jouranl of Projectiles，Rockets，Missile and Guidance，2012，5(32):117-120.)

［9］Luisa D Fairfax and Frank E Fresconi.Position Estimation for Projectiles Using Low-cost Sensors and Flight Dynamics［R］.America:Army Research Laboratory，2012.

［10］韓子鵬，等.彈箭外彈道學(xué)［M］.北京:北京理工大學(xué)出版社，2008:14，96-98.(Hang Zipeng，et al.Exterior Ballistics of Rockets and Missile［M］.Beijing:Beijing Institute of Technology Press，2008:14，96-98.)

［11］George M Siouris.Missle Guidance and Control Systems［M］.Beijing:National Defence Industry Press，2010:43.

Estimation of Ballistic Wind Speed Base on Aerodynamic Parameter Identification

ZHANG Yanru XIAO Liangang ZHANG Jisheng TIAN Feng CHEN Chang
Beijing Aerospace Automatic Control Institute，Beijing 100854，China

In the predictive process，guided munitions needed analysis of ballistic wind speed to improve the impact point precision．This paper introduced the wind distribution to determine wind speed estimation range．Then，according to the wind tunnel experimental data，the aerodynamic reference value of standard model and ballistic aerodynamic drag model was used for aerodynamic parameter identification by using the least square algorithm．This paper used the iterative algorithm to compute the approximate value of longitudinal wind，and built range correction model．According to the simulation test results，this algorithm could be used in real time correction，it was also suitable for engineering application and had the advantages of fast computing speed and small computation load．

Guided munitions;Aerodynamic parameter identification;Least square algorithm;Wind estimation

TJ412.+1

A

1006-3242(2014)02-0029-06

2013-03-12

張衍儒(1985-)，男，哈爾濱人，博士研究生，主要研究方向?yàn)橹茖?dǎo)彈藥控制系統(tǒng)綜合;肖練剛(1973-)，男，四川資中人，博士，研究員，主要研究方向?yàn)閷?dǎo)航、制導(dǎo)與控制;張繼生(1980-)，男，河北人，工程師，主要研究方向?yàn)榭刂葡到y(tǒng)綜合;田 豐(1985-)，男，江蘇人，工程師，主要研究方向控制系統(tǒng)綜合;陳 昌(1983-)，男，河北人，工程師，主要研究方向?yàn)榭刂葡到y(tǒng)綜合。