基于粒子群-牛頓算法的彈丸阻力系數(shù)辨識(shí)

史繼剛,易文俊,管 軍,劉世平

(南京理工大學(xué) a.瞬態(tài)物理國(guó)家重點(diǎn)實(shí)驗(yàn)室; b.能源與動(dòng)力工程學(xué)院，南京 210094)

【裝備理論與裝備技術(shù)】

基于粒子群-牛頓算法的彈丸阻力系數(shù)辨識(shí)

史繼剛a,易文俊a,管 軍a,劉世平b

(南京理工大學(xué) a.瞬態(tài)物理國(guó)家重點(diǎn)實(shí)驗(yàn)室; b.能源與動(dòng)力工程學(xué)院，南京 210094)

針對(duì)傳統(tǒng)牛頓迭代法在辨識(shí)彈丸氣動(dòng)參數(shù)時(shí)需要精確估計(jì)參數(shù)初值的問(wèn)題，提出了基于粒子群初值選取的牛頓迭代優(yōu)化算法辨識(shí)彈丸的零升阻力系數(shù)。采用彈丸的六自由度模型作為系統(tǒng)模型，以最大似然準(zhǔn)則作為辨識(shí)判據(jù)，結(jié)合粒子群算法的群體搜索性以及牛頓迭代法的局部細(xì)致搜索性，對(duì)辨識(shí)判據(jù)進(jìn)行了優(yōu)化，并且根據(jù)靈敏度計(jì)算分析了參數(shù)的可辨識(shí)性。通過(guò)仿真和實(shí)際數(shù)據(jù)辨識(shí)對(duì)算法的精確性和可靠性進(jìn)行了驗(yàn)證。仿真和實(shí)際辨識(shí)結(jié)果表明，該方法能有效地辨識(shí)旋轉(zhuǎn)彈丸零升阻力系數(shù)，可為進(jìn)一步提高射表精度、節(jié)省用彈量提供參考價(jià)值。

零升阻力系數(shù)辨識(shí)；粒子群；牛頓迭代；靈敏度

氣動(dòng)參數(shù)決定著彈箭飛行的彈道軌跡與穩(wěn)定性，是彈箭飛行控制領(lǐng)域的重要研究課題。目前，確定彈箭氣動(dòng)系數(shù)的方法主要有：理論計(jì)算法、風(fēng)洞實(shí)驗(yàn)法和射擊實(shí)驗(yàn)法。理論計(jì)算法過(guò)程復(fù)雜，不能完全考慮彈箭飛行的實(shí)際情況，計(jì)算結(jié)果和實(shí)際結(jié)果有一定誤差。風(fēng)洞實(shí)驗(yàn)法成本高，試驗(yàn)周期長(zhǎng)，難以模擬彈丸的各種運(yùn)動(dòng)，且模型支撐桿件和洞壁會(huì)影響測(cè)算結(jié)果。射擊實(shí)驗(yàn)法是針對(duì)彈箭的真實(shí)自由飛行數(shù)據(jù)進(jìn)行相關(guān)的參數(shù)辨識(shí)，是工程實(shí)踐中非常準(zhǔn)確、非常符合實(shí)際彈箭氣動(dòng)特性的方法。

目前射擊實(shí)驗(yàn)法中主要的參數(shù)辨識(shí)方法是基于最大似然準(zhǔn)則的辨識(shí)，優(yōu)化算法常為牛頓迭代法。例如汪清等[1]采用傳統(tǒng)牛頓迭代法對(duì)最大似然準(zhǔn)則優(yōu)化，辨識(shí)了高速自旋飛行器相關(guān)參數(shù)；崔平遠(yuǎn)等[2]用最大似然法對(duì)有控飛行器進(jìn)行了參數(shù)辨識(shí)；夏智勛等[3]采用最大似然法完成了軸對(duì)稱(chēng)飛行器的非線(xiàn)性氣動(dòng)參數(shù)辨識(shí)。另外，智能優(yōu)化算法也越來(lái)越多地被用于氣動(dòng)參數(shù)辨識(shí)研究中。王曉鵬[4]研究了自適應(yīng)遺傳算法辨識(shí)飛行器氣動(dòng)參數(shù)；錢(qián)煒祺等[5]用混合遺傳算法對(duì)某飛機(jī)橫向氣動(dòng)力參數(shù)進(jìn)行了辨識(shí)。杜昌平等[6]研究了蟻群算法辨識(shí)彈道參數(shù)；黃炯等[7]討論了模擬退火粒子群算法在飛機(jī)氣動(dòng)參數(shù)辨識(shí)中的應(yīng)用。本文將粒子群算法與傳統(tǒng)牛頓迭代法相結(jié)合，優(yōu)化最大似然準(zhǔn)則函數(shù)，利用彈丸速度數(shù)據(jù)辨識(shí)出高速旋轉(zhuǎn)彈丸的零升阻力系數(shù)，其中速度數(shù)據(jù)由全彈道坐標(biāo)跟蹤雷達(dá)測(cè)得。

1 粒子群-牛頓算法基本原理

最大似然法是目前在參數(shù)辨識(shí)領(lǐng)域應(yīng)用非常廣泛的方法之一，本文將粒子群-牛頓法與最大似然準(zhǔn)則相結(jié)合來(lái)獲取待辨識(shí)參數(shù)的最優(yōu)解[8]。圖1為粒子群-牛頓算法結(jié)合最大似然準(zhǔn)則辨識(shí)零升阻力系數(shù)的基本原理框圖。

圖1 辨識(shí)基本原理

下面將針對(duì)彈丸的自由飛行試驗(yàn)，分別給出氣動(dòng)參數(shù)辨識(shí)所需要的數(shù)學(xué)模型、辨識(shí)準(zhǔn)則以及尋優(yōu)方法。

1.1 高速旋轉(zhuǎn)彈丸運(yùn)動(dòng)模型

為了更加準(zhǔn)確地描述彈丸在空中的運(yùn)動(dòng)狀態(tài)并考慮高速旋轉(zhuǎn)彈丸的縱橫向運(yùn)動(dòng)相互耦合，采用6自由度彈道方程[9]作為參數(shù)辨識(shí)的基本模型：

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

式中各物理量的含義，參考文獻(xiàn)[9]。

1.2 參數(shù)辨識(shí)準(zhǔn)則

本文采用最大似然準(zhǔn)則作為目標(biāo)函數(shù)，其形式可表示為

(13)

1.3 基于粒子群-牛頓算法的尋優(yōu)

由旋轉(zhuǎn)彈丸的運(yùn)動(dòng)方程可知，當(dāng)彈丸參數(shù)、發(fā)射諸元、氣象諸元以及氣動(dòng)參數(shù)等確定以后，彈丸的運(yùn)動(dòng)狀態(tài)是確定的。反之，若已知彈丸的運(yùn)動(dòng)狀態(tài)，通過(guò)優(yōu)化的方法，總可以找到一個(gè)合適的阻力系數(shù)，使得計(jì)算值逼近實(shí)測(cè)數(shù)據(jù)。

1.3.1 粒子群算法

粒子群算法的數(shù)學(xué)描述是：假設(shè)在一個(gè)D維的目標(biāo)搜索空間中，有m個(gè)粒子組成一個(gè)種群，第i個(gè)粒子的位置向量用xi=[xi1,xi2,…,xiD]表示，飛行速度用vi=[vi1,vi2,…,viD]表示，第i個(gè)粒子搜索到的最優(yōu)位置為pi=[pi1,pi2,…,piD]，整個(gè)粒子群迄今搜索到的最優(yōu)位置為pg=[pg1,pg2,…,pgD]。粒子根據(jù)以下公式更新其速度和位置[10]：

vi(n+1)=ωvi(n)+c1r1(pi-xi(n))+c2r2(pg-xi(n))

(14)

xi(n+1)=xi(n)+vi(n)

(15)

ω=ωmax-t×[(ωmax-ωmin)/T]

(16)

式(14)和式(15)中：i=1,2,…,m；c1,c2是學(xué)習(xí)因子，為非負(fù)常數(shù)；r1,r2是均勻分布在(0,1)上的隨機(jī)數(shù)；ω為慣性權(quán)重；ωmax為最大權(quán)重；ωmin為最小權(quán)重；t為當(dāng)前迭代次數(shù)；T為總的迭代次數(shù)。粒子在解空間內(nèi)不斷跟蹤個(gè)體極值與全局極值進(jìn)行搜索，直至達(dá)到規(guī)定的迭代次數(shù)或滿(mǎn)足規(guī)定的誤差允許值為止。

粒子群算法主要的優(yōu)點(diǎn)是收斂速度快，便于工程實(shí)現(xiàn)，特別是在算法的早期。但同時(shí)算法后期也存在搜索精度較低、易發(fā)散等缺點(diǎn)。

1.3.2 牛頓迭代法

參數(shù)辨識(shí)的實(shí)質(zhì)，就是求辨識(shí)準(zhǔn)則的最值。根據(jù)極值原理，求出最大似然函數(shù)導(dǎo)數(shù)為零的值，也就辨識(shí)出了氣動(dòng)參數(shù)的值。在本文研究的問(wèn)題中，就是求解方程：

(17)

這就轉(zhuǎn)化為求非線(xiàn)性方程解的問(wèn)題。

牛頓法是一種常用的應(yīng)用導(dǎo)數(shù)求解非線(xiàn)性方程的算法，每一步的迭代方向都是沿著當(dāng)前點(diǎn)函數(shù)值下降的方向。因此牛頓迭代法十分依賴(lài)初值的選擇，在計(jì)算初值選定合理的情況下，能較快較準(zhǔn)確地計(jì)算出方程的解。本文中牛頓迭代公式為

(18)

其中涉及到求解?v/?cx0(靈敏度)，在下文可辨識(shí)性分析中將作具體介紹。

1.3.3 粒子群-牛頓算法步驟

針對(duì)粒子群優(yōu)化算法和牛頓迭代法的不足，綜合各自的優(yōu)勢(shì)，提出了基于粒子群初始位置選取和牛頓法局部精確迭代的混合算法，其流程如下：

1) 在初始化范圍內(nèi)，對(duì)粒子群進(jìn)行隨機(jī)初始化，包括粒子的位置和速度。

2) 評(píng)價(jià)每個(gè)微粒的適應(yīng)度，將當(dāng)前各微粒的位置和適應(yīng)值存儲(chǔ)在各微粒的pbest中，將pbest中適應(yīng)值最優(yōu)個(gè)體的位置和適應(yīng)值存儲(chǔ)于gbest中。

3) 用式(14)、式(15)進(jìn)行粒子速度和位移更新。

4) 對(duì)每個(gè)微粒，將其適應(yīng)值與其經(jīng)歷過(guò)的最好位置作比較，如果較好，則將其作為當(dāng)前的最好位置。

5) 比較當(dāng)前所有pbest和gbest的值，更新gbest。

6) 若達(dá)到終止條件(達(dá)到最大迭代次數(shù))，返回當(dāng)前全局最優(yōu)個(gè)體Pg，轉(zhuǎn)向步驟7)；否則，k=k+1，轉(zhuǎn)步驟2)。

7) 牛頓迭代法。以步驟5)返回的當(dāng)前全局最優(yōu)個(gè)體Pg為牛頓法的初始點(diǎn)進(jìn)行迭代。

8) 若達(dá)到終止條件則結(jié)束，輸出的當(dāng)前結(jié)果作為所求問(wèn)題的最優(yōu)解；否則將解算值作為初始值繼續(xù)執(zhí)行牛頓算法。

2 可辨識(shí)性分析

參數(shù)辨識(shí)研究領(lǐng)域中，參數(shù)可辨識(shí)性分析的重要依據(jù)就是靈敏度分析(即影響度分析)[11-12]。觀(guān)測(cè)量關(guān)于待辨識(shí)參數(shù)的靈敏度越高，二者的相關(guān)度越高，即用觀(guān)測(cè)量辨識(shí)出該參數(shù)的可能性越大、可靠性越高；反之，二者的相關(guān)度低，即很難用觀(guān)測(cè)量有效辨識(shí)出該參數(shù)。本節(jié)分別以速度數(shù)據(jù)和側(cè)偏數(shù)據(jù)作為觀(guān)測(cè)量，對(duì)相應(yīng)的狀態(tài)方程關(guān)于cx0求偏導(dǎo)得到靈敏度方程，聯(lián)合狀態(tài)方程組，計(jì)算并分析不同觀(guān)測(cè)量關(guān)于零升阻力系數(shù)的靈敏度。靈敏度方程的推導(dǎo)過(guò)程如下：

定義：式(1)、式(12)中各微分方程中左側(cè)的物理量對(duì)cx0的敏感系數(shù)分別為：

以速度v對(duì)零升阻力系數(shù)cx0的靈敏度P11的靈敏度函數(shù)推導(dǎo)如下

(19)

以側(cè)偏z對(duì)零升阻力系數(shù)cx0的靈敏度P121的靈敏度函數(shù)推導(dǎo)如下

(20)

式(19)中：

A17=-mg(cosθacosψ2P21-sinθasinψ2P31)

A11～A17中的一些中間變量，由于篇幅有限，文中不再給出，讀者如需要，可自行推導(dǎo)。只要推出了所有的靈敏度方程，然后和六自由度彈道方程聯(lián)合求解，即可得到各靈敏度隨時(shí)間變化的數(shù)值解。圖2分別為靈敏度P11、P121隨時(shí)間變化的數(shù)值解。

圖2 靈敏度P11和 P121

從圖2可以看出，不同觀(guān)測(cè)量辨識(shí)參數(shù)的靈敏度差別很大。側(cè)偏z對(duì)cx0的靈敏度P121的數(shù)值過(guò)小，所以基本不能用側(cè)偏z來(lái)對(duì)cx0進(jìn)行辨識(shí)。速度v對(duì)cx0的敏感系數(shù)P11的數(shù)值很大，可以通過(guò)v對(duì)cx0進(jìn)行辨識(shí)，后續(xù)的實(shí)驗(yàn)也證明了該結(jié)論的正確性。

3 實(shí)驗(yàn)驗(yàn)證

3.1 仿真驗(yàn)證

已知某彈丸隨馬赫數(shù)變化的氣動(dòng)參數(shù)，代入彈丸的運(yùn)動(dòng)模型中計(jì)算得到彈丸飛行數(shù)據(jù)，在計(jì)算數(shù)據(jù)基礎(chǔ)上加上信噪比為100∶5隨機(jī)白噪聲(模擬測(cè)量誤差)作為模擬彈道，根據(jù)模擬彈道數(shù)據(jù)，利用本文介紹的方法辨識(shí)零升阻力系數(shù)。最后比較給定的零升阻力系數(shù)和計(jì)算的零升阻力系數(shù)。

算例中選取了4段速度數(shù)據(jù)，例如在1馬赫時(shí)，選取1馬赫對(duì)應(yīng)時(shí)間點(diǎn)前后共50個(gè)點(diǎn)的速度數(shù)據(jù)辨識(shí)零升阻力系數(shù)。其他速度的數(shù)據(jù)選取同理。表1為不同馬赫數(shù)下的測(cè)試結(jié)果。

表1 零升阻力系數(shù)辨識(shí)結(jié)果

上述測(cè)試結(jié)果表明，利用速度數(shù)據(jù)，基于最大似然準(zhǔn)則的粒子群-牛頓算法可以有效地對(duì)零升阻力系數(shù)進(jìn)行辨識(shí)。

3.2 實(shí)際數(shù)據(jù)處理

利用小區(qū)間常數(shù)法[13]，認(rèn)為在小區(qū)間內(nèi)零升阻力系數(shù)是不變的，所以將全彈道的測(cè)量數(shù)據(jù)分成若干個(gè)數(shù)據(jù)段，分別對(duì)每段進(jìn)行辨識(shí)，然后再插值得到全彈道的零升阻力系數(shù)。將辨識(shí)結(jié)果代入彈丸運(yùn)動(dòng)數(shù)學(xué)模型計(jì)算得到彈丸速度數(shù)據(jù)，再與雷達(dá)測(cè)得的實(shí)際速度數(shù)據(jù)進(jìn)行比較。

圖3中，橫坐標(biāo)是馬赫數(shù)，縱坐標(biāo)是零升阻力系數(shù)。

圖3 實(shí)際辨識(shí)的零阻隨馬赫數(shù)變化曲線(xiàn)

用該零阻系數(shù)代入彈道方程計(jì)算得到的速度數(shù)據(jù)和實(shí)測(cè)雷達(dá)數(shù)據(jù)的對(duì)比如圖4所示。

圖4 計(jì)算值和測(cè)量值的速度曲線(xiàn)

通過(guò)數(shù)據(jù)的分析，全彈道的擬合概率誤差在6.34 m/s左右，有效地驗(yàn)證了算法和氣動(dòng)辨識(shí)結(jié)果的正確性。

4 結(jié)論

本文提出粒子群-牛頓算法將粒子群算法與傳統(tǒng)牛頓迭代法相結(jié)合，并根據(jù)最大似然準(zhǔn)則，利用彈丸速度數(shù)據(jù)辨識(shí)出高速旋轉(zhuǎn)彈丸的零升阻力系數(shù)。通過(guò)仿真驗(yàn)證和實(shí)際數(shù)據(jù)辨識(shí)的結(jié)果可以看出，基于最大似然準(zhǔn)則的粒子群-牛頓算法能夠有效地辨識(shí)出旋轉(zhuǎn)彈丸的零升阻力系數(shù)。本課題的研究，對(duì)進(jìn)一步提高射表精度、節(jié)省用彈量有著重要的參考價(jià)值。

[1] 汪清,錢(qián)煒祺,萬(wàn)宗國(guó),等.高速自旋飛行器氣動(dòng)參數(shù)辨識(shí)[J].空氣動(dòng)力學(xué)學(xué)報(bào),2004,22(1):1-6.

[2] 崔平遠(yuǎn),楊滌,吳瑤華,等.極大似然法及在有控飛行器氣動(dòng)參數(shù)辨識(shí)中的應(yīng)用[J].航空學(xué)報(bào),1991,12(11):644-649.

[3] 夏智勛,盛湘饒,唐國(guó)金.用極大似然法辨識(shí)軸對(duì)稱(chēng)飛行器的氣動(dòng)參數(shù)[J].航空學(xué)報(bào),1998,19(5):564-566.

[4] WANG Xiaopeng.Aerodynamic Parameter Identification of Flight Vehicles Based on Adaptive Genetic Algorithm[J].Journal of Astronautics,2003,24(3):303-307.

[5] 錢(qián)煒祺,汪清,何開(kāi)鋒,等.混合遺傳算法在氣動(dòng)力參數(shù)辨識(shí)中的應(yīng)用[J].飛行力學(xué),2004,22(1):33-36.

[6] 杜昌平,孫化東,周德云,等.基于蟻群算法的彈道參數(shù)辨識(shí)方法研究[J].電光與控制,2008,15(7):4-6.

[7] 黃炯,艾劍良,萬(wàn)婧.基于模擬退火粒子群算法的飛機(jī)氣動(dòng)參數(shù)辨識(shí)[J].復(fù)旦學(xué)報(bào),2016,55(3):336-341.

[8] 張?zhí)戽?錢(qián)煒祺,何開(kāi)鋒,等.基于最大似然法的風(fēng)洞自由飛試驗(yàn)氣動(dòng)力參數(shù)辨識(shí)技術(shù)研究[J].實(shí)驗(yàn)流體力學(xué),2015,29(5):8-13.

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

[10]姚金杰,韓焱,王黎明,等.基于分組粒子群的飛行目標(biāo)定位方法[J].探測(cè)與控制學(xué)報(bào),2012,34(5):46-49.

[11]王欣,張可菊,張鳳閣,等.彈箭彈道參數(shù)相對(duì)于線(xiàn)性阻力系數(shù)的靈敏度計(jì)算[J].四川兵工學(xué)報(bào),2008,29(4):1-4.

[12]張健,姚俊,王欣,等.基于彈道參數(shù)靈敏度函數(shù)一種阻力系數(shù)辨識(shí)方法[J].彈箭與制導(dǎo)學(xué)報(bào),2004,25(1):47-49.

[13]閆章更，祁載康.射表技術(shù)[M].北京：國(guó)防工業(yè)出版社， 2000.

(責(zé)任編輯 周江川)

Drag Coefficient Identification of Spinning Projectile Using Particle Swarm Newton Iteration Method

SHI Ji-ganga，YI Wen-juna，GUAN Juna， LIU Shi-pingb

(a.National Key Laboratory of Transient Physics; b.School of Energy and Power Engineering, Nanjing University of Science and Technology, Nanjing 210094, China)

Traditional aerodynamic parameter identification using Newton algorithm has strict requirement for initial parameter. Combining the advantages of particle swarm algorithm in initial value selection and the advantages of Newton iteration method in precise iteration, this paper presented a new Newton iteration algorithm using the initial value from particle swarm algorithm to identify the spinning projectile’s drag coefficient. Then, using the 6 degree model and maximum likelihood rule, we analyzed the sensitivity. Finally processing the actual speed data of spinning projectile, we acquired a satisfied result. According to the analysis of the result, it’s concluded that the algorithm can effectively get a satisfied result in the process of parameter identification, which is of great value in compiling fire table and saving the bomb.

drag coefficient identification; particle swarm algorithm; Newton iteration; sensitivity

2016-09-25；

2016-10-29

國(guó)家自然科學(xué)基金項(xiàng)目(11472136,11402117)

史繼剛(1992—)，男，碩士，主要從事彈箭飛行控制研究。

易文俊(1970—)，男，博士，教授，主要從事彈箭飛行控制研究。

10.11809/scbgxb2017.02.006

史繼剛,易文俊,管軍,等.基于粒子群-牛頓算法的彈丸阻力系數(shù)辨識(shí)[J].兵器裝備工程學(xué)報(bào)，2017(2):23-26.

format:SHI Ji-gang，YI Wen-jun，GUAN Jun,et al.Drag Coefficient Identification of Spinning Projectile Using Particle Swarm Newton Iteration Method[J].Journal of Ordnance Equipment Engineering,2017(2):23-26.

TJ3

A

2096-2304(2017)02-0023-04