侯現欽,王良明,傅 健
(南京理工大學能源與動力工程學院, 南京 210094)
氣動參數辨識是系統辨識中發展最成熟的領域之一,已經成為彈箭研制、評估和控制過程中重要的組成部分[1]。通過試驗測量數據辨識出的氣動參數更能反映出彈丸真實飛行的狀態,比風洞試驗獲取的結果更加精確,同時還可以檢驗彈丸的飛行穩定性[2]。由于目前國內還沒有大口徑的室內靶道,高旋彈的飛行速度、位置尤其是姿態數據較難獲得,所以高旋彈的氣動參數辨識一直是系統辨識的一個難點。
極大似然估計是目前精度較高的辨識方法,在氣動參數辨識中被廣泛應用。汪清等[3]利用極大似然估計對高速自旋飛行器的氣動參數進行了辨識;郭德龍、周永權[4]使用極大似然估計以威布爾三參數分布為例進行參數估計;王曉鵬等[5]將極大似然估計應用到戰術導彈的非線性氣動力參數辨識中;張天蛟等[6]采用最大似然辨識算法對風洞自由飛試驗數據進行氣動力參數辨識。
隨著智能算法的快速發展,差分進化算法由于其卓越的計算性能在辨識方面得到了廣泛應用。簡兆圣、艾劍良[7],袁瑞俠[8]等使用差分進化算法分別對飛機縱向運動氣動力參數、VTOL飛行器參數進行辨識。
文中針對傳統算法的初值敏感問題,將差分進化算法應用到高旋彈的氣動參數辨識中,以零升阻力系數和極阻尼力矩系數為例驗證該算法的實用性。
高速旋轉彈丸在實際運動中,質心運動和姿態運動實際是相互關聯和影響的,為了更精確的描述高速旋轉彈丸的運動規律,文中采用六自由度剛體彈道方程[9]作為參數辨識的理論模型。其具體表達式如下:
(1)
(2)
式中:x、y、z為彈丸的位置坐標;v為速度數據;θa、ψ2、φa、φ2分別為彈丸的速度高低角、速度方位角、彈軸高低角和彈軸方位角;γ為轉角;ωξ、ωη、ωζ分別為滾轉角速度、俯仰角速度和偏航角速度;Fx2、Fy2、Fz2和Mξ、Mη、Mζ分別為作用在彈丸上的力和力矩在3個坐標軸下的分量;A、C分別為赤道轉動慣量和極轉動慣量;δ1、δ2、β分別為高低攻角、方向攻角和側滑角。
差分進化算法(DE)由Storn和Price于1995年首次提出,該算法采用浮點數編碼,使用“貪婪”選擇策略,對特定優化問題的特征信息不敏感。與其他進化算法相比,更適應于多變量優化,收斂更快。其主要步驟如下所示。
傳統的迭代優化算法需要通過先驗知識來確定待辨識參數的初值,如果初值選取偏差較大,計算過程中可能會出現局部收斂或者發散現象。DE算法無需選擇初值,只需根據以往經驗設置待辨識參數的大致范圍即可,所以設置阻力系數的初始種群的上下限分別為BU=0、BL=1(極阻尼力矩系數類似),在此范圍內隨機生成初始種群,種群規模為10。由于文中只選擇一組數據作為觀測量,所以向量維數D=1。
極大似然準則函數向最小方向發展,文中研究內容中,準則函數J可能出現負數情況,所以設置適應度函數為f=-J+500,故選取適應度函數值較大的待辨識參數繼續變異進化。該算法收斂較快,一般在30~40代之間就會收斂,但是為了避免部分數據點出現局部收斂,將進化代數G設置為100。
差分進化算法也是通過變異、交叉、選擇操作逐步尋求最優解,其優化流程圖如圖1所示[10-11]。

圖1 DE算法優化流程圖
1)差分變異
種群個體中的差分向量經過縮放后,與種群內另外的相異個體相加得到變異向量。變異方程為:
vi,g=xr1,g+F×(xr2,g-xr3,g)
(3)
式中:r1≠r2≠r3≠i;g為進化代數;F為變異因子,主要影響算法的全局尋優能力,F越小算法對局部的搜索能力越好,F越大算法越能跳出局部極小點,但是收斂速度會變慢。設置F=0.7。
2)二項式交叉
通過隨機選擇,使實驗向量至少有一個分量由變異向量貢獻。二項式交叉操作的方程為:
(4)
式中:j=1,2,3,…,D。Cr為交叉因子,主要反映的是在交叉的過程中,子代與父代、中間變異體之間交換信息量的大小程度。Cr的值越大,信息量交換的程度越大。反之,如果Cr的值偏小,將會使種群的多樣性快速減小,不利于全局尋優。文中取Cr=0.6。
3)選擇
差分進化算法采用“貪婪”的選擇策略,根據各向量的適應度值來選擇最優個體,方程如式(5)所示。
(5)
式中:xi,g+1為下一代的目標向量。
(6)
也可以取似然函數為ln(p(Y|θ))。
文中待辨識參數矢量為:
θ=[cx0M′xz]
(7)
實際測量值矢量為:
ym(i)=[xyzvωξ]
(8)
極大似然準則函數具體表示形式如式(9)所示[12]。
(9)
極大似然準則將參數辨識問題轉化為求極值問題,利用優化算法將準則函數達到最小值。文中采用DE優化算法對似然函數進行優化。具體流程如圖2所示。
靈敏度可以衡量觀測數據和待辨識參數之間相關性的大小,靈敏度值越大,則相關性越大,反之越小[13]。下面以零升阻力系數為例,介紹靈敏度求解過程。各狀態量對零升阻力系數的靈敏度如式(10)、式(11)所示。
(10)
因為該偏微分的二階偏導是連續的,速度的靈敏度方程可作如下轉換:
(11)
其他狀態量做相同變換,將所有的靈敏度方程與六自由度彈道方程聯立進行求解,即可計算出每個狀態量對零升阻力系數的靈敏度值。通過對比將速度和滾轉角速度作為辨識零升阻力系數和極阻尼力矩系數的觀測量。

圖2 氣動參數辨識流程圖
以155 mm某高旋彈(初速為930 m/s,質量為45.4 kg,射程為32 000 m左右,)為例進行全彈道氣動參數辨識。分別采用無噪聲數據和加噪聲數據來驗證該算法的有效性、精確性和優越性。
1)利用六自由度彈道方程組計算得到彈道數據作為觀測量進行參數辨識。通過與標準氣動參數對比來驗證該算法的可行性。
采用小區間常數法,將彈道分為多個小區間,每個區間取10個數據點,計算步長為0.005 s,將每個區間的氣動參數看作一個常數。DE算法將每個小區間中10個氣動參數的平均值作為辨識結果。隨著代數的增加,DE算法快速收斂,圖3為辨識阻力系數時第一個小區間的平均適應度值隨代數的變化情況,可以看出收斂非常快,在30~40代之間就已經完成收斂。

圖3 適應度變化曲線
表1、表2分別為兩種方法辨識結果對比情況(隨機選取5個數據點)。圖4、圖5分別為該算法、傳統算法對零升阻力系數和極阻尼力矩系數導數的辨識結果與標準氣動參數對比圖??梢钥闯觯珼E算法對零升阻力系數和極阻尼力矩系數的辨識誤差明顯低于傳統算法。由于傳統梯度優化算法對初值要求高,容易出現局部收斂,導致辨識結果變化幅度較大。而DE算法的辨識結果更加平滑,更加接近標準值。
2)為測試所提方法對含噪聲數據的辨識精度和工程應用性,將彈道數據加上信噪比為100∶5的高斯白噪聲進行辨識。圖6、圖7分別為該算法對零升阻力系數和極阻尼力矩系數導數辨識結果與理論值的對比曲線圖??梢钥闯霰孀R結果在標準值上下震動,隨著馬赫數的增加誤差逐漸減小。

表1 零升阻力系數辨識結果對比

表2 極阻尼力矩系數導數辨識結果對比

圖4 零升阻力系數辨識結果對比圖

圖5 極阻尼力矩系數導數辨識結果對比圖
工程上通常將辨識的氣動參數代入到彈道方程,計算出各觀測量,與測量數據進行比較來驗證辨識結果的準確性[14]。為進一步驗證該方法的性能,文中采用該方法進行反驗算。圖8、圖9分別為速度和滾轉角速度的測量值和用辨識結果反算的計算值的對比??梢钥闯觯此阒蹬c測量值之間的誤差非常小,測量彈道的落點坐標為(31 755.558,1.706,1 383.374),反算得到的彈道的落點坐標為(31 767.304,1.476,1 383.731),兩者之間的距離為11.753 7 m。相對于高旋彈(如榴彈),辨識結果已具有相當高的精度。

圖6 零升阻力系數辨識結果對比圖

圖7 極阻尼力矩系數導數辨識結果對比圖

圖8 速度對比圖

圖9 滾轉角速度對比圖
對高旋彈的零升阻力系數和極阻尼力矩系數進行了辨識,通過將差分進化-極大似然方法應用到該類型彈的辨識中取得了較好的辨識結果,得到如下結論:
1)提出的差分進化-極大似然方法的辨識精度相對于傳統方法精度更高。
2)該方法成功解決了傳統算法對初值選取敏感問題。
3)也可以嘗試將DE算法與其他傳統辨識方法相結合(解決初值問題),進一步提高辨識精度。