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

(南京理工大學 a.瞬態(tài)物理國家重點實驗室; b.能源與動力工程學院，南京 210094)

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

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

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

(南京理工大學 a.瞬態(tài)物理國家重點實驗室; b.能源與動力工程學院，南京 210094)

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

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

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

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

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

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

圖1 辨識基本原理

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

1.1 高速旋轉(zhuǎn)彈丸運動模型

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

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

式中各物理量的含義，參考文獻[9]。

1.2 參數(shù)辨識準則

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

(13)

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

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

1.3.1 粒子群算法

粒子群算法的數(shù)學描述是：假設(shè)在一個D維的目標搜索空間中，有m個粒子組成一個種群，第i個粒子的位置向量用xi=[xi1,xi2,…,xiD]表示，飛行速度用vi=[vi1,vi2,…,viD]表示，第i個粒子搜索到的最優(yōu)位置為pi=[pi1,pi2,…,piD]，整個粒子群迄今搜索到的最優(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是學習因子，為非負常數(shù)；r1,r2是均勻分布在(0,1)上的隨機數(shù)；ω為慣性權(quán)重；ωmax為最大權(quán)重；ωmin為最小權(quán)重；t為當前迭代次數(shù)；T為總的迭代次數(shù)。粒子在解空間內(nèi)不斷跟蹤個體極值與全局極值進行搜索，直至達到規(guī)定的迭代次數(shù)或滿足規(guī)定的誤差允許值為止。

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

1.3.2 牛頓迭代法

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

(17)

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

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

(18)

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

1.3.3 粒子群-牛頓算法步驟

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

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

2) 評價每個微粒的適應(yīng)度，將當前各微粒的位置和適應(yīng)值存儲在各微粒的pbest中，將pbest中適應(yīng)值最優(yōu)個體的位置和適應(yīng)值存儲于gbest中。

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

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

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

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

7) 牛頓迭代法。以步驟5)返回的當前全局最優(yōu)個體Pg為牛頓法的初始點進行迭代。

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

2 可辨識性分析

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

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

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

(19)

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

(20)

式(19)中：

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

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

圖2 靈敏度P11和 P121

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

3 實驗驗證

3.1 仿真驗證

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

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

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

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

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

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

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

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

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

圖4 計算值和測量值的速度曲線

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

4 結(jié)論

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

[1] 汪清,錢煒祺,萬宗國,等.高速自旋飛行器氣動參數(shù)辨識[J].空氣動力學學報,2004,22(1):1-6.

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

[3] 夏智勛,盛湘饒,唐國金.用極大似然法辨識軸對稱飛行器的氣動參數(shù)[J].航空學報,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] 錢煒祺,汪清,何開鋒,等.混合遺傳算法在氣動力參數(shù)辨識中的應(yīng)用[J].飛行力學,2004,22(1):33-36.

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

[7] 黃炯,艾劍良,萬婧.基于模擬退火粒子群算法的飛機氣動參數(shù)辨識[J].復(fù)旦學報,2016,55(3):336-341.

[8] 張?zhí)戽?錢煒祺,何開鋒,等.基于最大似然法的風洞自由飛試驗氣動力參數(shù)辨識技術(shù)研究[J].實驗流體力學,2015,29(5):8-13.

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

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

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

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

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

(責任編輯 周江川)

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

國家自然科學基金項目(11472136,11402117)

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

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

10.11809/scbgxb2017.02.006

史繼剛,易文俊,管軍,等.基于粒子群-牛頓算法的彈丸阻力系數(shù)辨識[J].兵器裝備工程學報，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