李金晟,常思江,陳升富

(1.北京航天長(zhǎng)征飛行器研究所,北京 100076;2.南京理工大學(xué) 能源與動(dòng)力工程學(xué)院,江蘇 南京 210094)

氣動(dòng)參數(shù)對(duì)于常規(guī)彈箭的飛行穩(wěn)定性及彈道軌跡具有重要影響。目前,彈箭的氣動(dòng)參數(shù)的估計(jì)方法已有許多,其中經(jīng)典常用的有基于最小二乘準(zhǔn)則的C-K法[1]、基于極大似然準(zhǔn)則的輸出誤差法[2]以及結(jié)合了濾波思想的濾波誤差法[3]。但是,這些方法在應(yīng)用于實(shí)測(cè)飛行數(shù)據(jù)時(shí),無(wú)一例外都會(huì)出現(xiàn)一些復(fù)雜的問(wèn)題。其中較為普遍的一個(gè)難題便是對(duì)于真實(shí)系統(tǒng)缺少精確的數(shù)學(xué)模型。大多數(shù)傳統(tǒng)的方法利用動(dòng)力學(xué)模型來(lái)近似描述真實(shí)系統(tǒng),這種近似會(huì)帶來(lái)建模誤差的問(wèn)題。但對(duì)于建模誤差的處理目前還沒(méi)有成熟的理論可用,通常是將其簡(jiǎn)單的歸為狀態(tài)或測(cè)量噪聲來(lái)處理,這種處理方式?jīng)]有考慮到建模誤差可能是確定性的而非隨機(jī)性的[4]。為了解決上述傳統(tǒng)方法存在的建模誤差問(wèn)題,相關(guān)人員開(kāi)始采用人工神經(jīng)網(wǎng)絡(luò)(artificial neural networks,ANN)進(jìn)行氣動(dòng)參數(shù)辨識(shí),其完全基于對(duì)輸入輸出的測(cè)量而不考慮系統(tǒng)內(nèi)部的物理機(jī)制。

基于神經(jīng)網(wǎng)絡(luò)的氣動(dòng)參數(shù)辨識(shí)方法主要有三類:一是基于動(dòng)力學(xué)模型的辨識(shí)方法[5],其優(yōu)點(diǎn)是無(wú)需預(yù)先計(jì)算基準(zhǔn)氣動(dòng)數(shù)據(jù),缺點(diǎn)是所辨識(shí)的氣動(dòng)模型仍是一個(gè)“黑箱”,無(wú)法得到具體表達(dá)式;二是基于穩(wěn)定或控制導(dǎo)數(shù)的Delta法和Zeros法[6],其優(yōu)點(diǎn)是可以辨識(shí)出氣動(dòng)模型,缺點(diǎn)是需預(yù)先得知較為準(zhǔn)確的氣動(dòng)參數(shù)數(shù)據(jù),從而利用其訓(xùn)練神經(jīng)網(wǎng)絡(luò);三是將前饋神經(jīng)網(wǎng)絡(luò)(feed forward neural networks,FFNNs)與優(yōu)化算法高斯-牛頓法(Gauss-Newton)相結(jié)合的神經(jīng)網(wǎng)絡(luò)-牛頓(neural Gauss-Newton,NGN)辨識(shí)方法[7],其獨(dú)特之處在于無(wú)需預(yù)先對(duì)系統(tǒng)進(jìn)行建模,不需求解運(yùn)動(dòng)方程。這三類方法已經(jīng)成功應(yīng)用于飛機(jī)、導(dǎo)彈、再入體等飛行器的氣動(dòng)參數(shù)辨識(shí),但尚沒(méi)有應(yīng)用于彈丸。

本文針對(duì)某無(wú)控旋轉(zhuǎn)彈丸的阻力系數(shù)辨識(shí)問(wèn)題,采用基于誤差反向傳播(back propagation,BP)神經(jīng)網(wǎng)絡(luò)的NGN法進(jìn)行研究,并對(duì)結(jié)果進(jìn)行了分析。

1 BP網(wǎng)絡(luò)辨識(shí)動(dòng)力學(xué)模型

1.1 BP網(wǎng)絡(luò)原理簡(jiǎn)介

BP神經(jīng)網(wǎng)絡(luò)是一種單向傳播的多層前向網(wǎng)絡(luò),其分為輸人層、隱含層和輸出層。隱含層可以有若干層,每一層的神經(jīng)元只接收前一層神經(jīng)元的輸出,同一層神經(jīng)元之間不連接。該網(wǎng)絡(luò)的拓?fù)浣Y(jié)構(gòu)如圖1所示[8]。

1.2 辨識(shí)動(dòng)力學(xué)模型

NGN算法區(qū)別于傳統(tǒng)辨識(shí)算法的獨(dú)特之處在于,利用神經(jīng)網(wǎng)絡(luò)強(qiáng)大的非線性逼近能力來(lái)減小建模誤差。由于本文主要是辨識(shí)彈丸的零升阻力系數(shù)Cx0,且文獻(xiàn)[2]已表明通過(guò)量測(cè)速度數(shù)據(jù)v便可實(shí)現(xiàn)對(duì)Cx0的辨識(shí)。故利用網(wǎng)絡(luò)辨識(shí)出速度狀態(tài)滿足的關(guān)系式即可。

BP神經(jīng)網(wǎng)絡(luò)在對(duì)狀態(tài)進(jìn)行預(yù)測(cè)之前首先要對(duì)網(wǎng)絡(luò)進(jìn)行多次訓(xùn)練,使其具備聯(lián)想記憶和逼近能力,具體的訓(xùn)練過(guò)程包括以下幾個(gè)步驟[9]。

第1步:網(wǎng)絡(luò)初始化。由質(zhì)點(diǎn)彈道模型可得,加速度主要與速度、密度及零升阻力系數(shù)有關(guān),而在標(biāo)準(zhǔn)條件下密度又是高度的函數(shù)。因此本文的輸入向量為

X(t)=(vx(t)vy(t)vz(t)y(t)Cx0(t))

(1)

式中:t為離散的采樣時(shí)間,Cx0(t)是通過(guò)對(duì)風(fēng)洞吹風(fēng)或理論計(jì)算所得的基準(zhǔn)氣動(dòng)參數(shù)進(jìn)行線性插值所得。基于X(t),輸出向量為

Y(t+1)=(vx(t+1)vy(t+1)vz(t+1)y(t+1))

(2)

由輸入、輸出向量可得輸入層節(jié)點(diǎn)n=5,輸出層節(jié)點(diǎn)數(shù)m=4。

第2步:求隱含層輸出H。由X(t),輸入層和隱含層間連接權(quán)值wij及隱含層閾值a,求H。

(3)

式中:l為隱含層節(jié)點(diǎn)數(shù);f(*)為隱含層激勵(lì)函數(shù),本文設(shè)置其為Sigmoid函數(shù),表達(dá)式為

(4)

式中:β>0,本文取β=1。

第3步:求輸出層輸出Y(t+1)。由H,隱含層與輸出層連接權(quán)值wjk及輸出層閾值b,求預(yù)測(cè)輸出Y(t+1)。

(5)

第4步:誤差計(jì)算。由Y(t+1)和期望輸出O(t+1),求出預(yù)測(cè)誤差e及誤差性能指標(biāo)函數(shù)E。

(6)

第5步:權(quán)值與閾值更新。基于最速下降法,由e更新連接權(quán)值wij,wjk以及閾值a和b。

(7)

(8)

η(g)=ηmax-[g(ηmax-ηmin)/gmax]

(9)

式中:ηmax,ηmin分別為最大、最小學(xué)習(xí)速率;gmax為最大迭代次數(shù);g為當(dāng)前迭代次數(shù)。

第6步:判斷算法迭代是否結(jié)束,若沒(méi)有結(jié)束,返回第2步。為了提高迭代效率,本文設(shè)置了2個(gè)迭代終止準(zhǔn)則,準(zhǔn)則一為達(dá)到最大迭代次數(shù)gmax,準(zhǔn)則二為誤差性能函數(shù)E小于給定精度。

神經(jīng)網(wǎng)絡(luò)訓(xùn)練結(jié)束后,所得到的BP網(wǎng)絡(luò)就可以看成一個(gè)非線性動(dòng)力學(xué)模型,其表達(dá)了從n個(gè)輸入自變量到m個(gè)輸出因變量的函數(shù)映射關(guān)系。

2 NGN算法應(yīng)用過(guò)程

當(dāng)神經(jīng)網(wǎng)絡(luò)的相關(guān)參數(shù)調(diào)試恰當(dāng)后,就可以對(duì)t時(shí)刻的任意測(cè)量輸入變量X(t)進(jìn)行預(yù)測(cè),得到t+1時(shí)刻的輸出變量Y(t+1)。接下來(lái)便可根據(jù)神經(jīng)網(wǎng)絡(luò)模型辨識(shí)所需的氣動(dòng)參數(shù)。

2.1 分段辨識(shí)思想

本文基于分段辨識(shí)思想提取阻力系數(shù),即將雷達(dá)所測(cè)的速度數(shù)據(jù)劃分為若干個(gè)小區(qū)間,在每一個(gè)小區(qū)間內(nèi)由于速度變化很小,故可把該區(qū)間內(nèi)的Cx0視為常數(shù)。

2.2 辨識(shí)準(zhǔn)則

本文采用極大似然準(zhǔn)則做為目標(biāo)準(zhǔn)則函數(shù)。在實(shí)際應(yīng)用中由于彈丸自由飛行試驗(yàn)要求較為嚴(yán)格,這時(shí)系統(tǒng)的過(guò)程噪聲很小,故目標(biāo)函數(shù)可表示為

(10)

V(q)=ym(q)-y(q)

(11)

式中:ym(q)=(vm,x(q)vm,y(q)vm,z(q))T為已經(jīng)過(guò)坐標(biāo)轉(zhuǎn)換的雷達(dá)量測(cè)速度數(shù)據(jù);y(q)=(vx(q)vy(q)vz(q))T是基于當(dāng)前輸入,通過(guò)BP神經(jīng)網(wǎng)絡(luò)所計(jì)算的理論速度值。

2.3 靈敏度參數(shù)計(jì)算

θ(r+1)=θ(r)+Δθ(r)

(12)

傳統(tǒng)極大似然法是采用數(shù)值積分法聯(lián)立求解彈道微分方程組與靈敏度方程組,從而獲得靈敏度系數(shù)值。由于需要建立靈敏度方程組,故其較為繁瑣,尤其是當(dāng)模型比較復(fù)雜時(shí),計(jì)算量會(huì)急劇增大。而NGN憑借其獨(dú)特優(yōu)勢(shì),極大地簡(jiǎn)化了計(jì)算靈敏度系數(shù)的過(guò)程,其方式如下:

(13)

式中:d為辨識(shí)系數(shù)個(gè)數(shù);對(duì)于本文而言d=1;Ypk(t+1)為擾動(dòng)輸出,即給當(dāng)前輸入θd加一任意小擾動(dòng)Δθd,利用已訓(xùn)練好的神經(jīng)網(wǎng)絡(luò)計(jì)算得到的輸出值。

至此,利用式(13)計(jì)算的靈敏度值可求出修正量ΔCx0。重復(fù)以上的計(jì)算過(guò)程,每迭代一次都需要由似然準(zhǔn)則函數(shù)計(jì)算Qr,直到當(dāng)|1-Qr/Qr-1|≤0.01時(shí),迭代收斂,計(jì)算結(jié)束。此時(shí)輸出的Cx0就是該彈道段內(nèi)辨識(shí)計(jì)算所得到的阻力參數(shù)。

3 數(shù)值結(jié)果分析

3.1 仿真驗(yàn)證

為了檢驗(yàn)NGN算法的可行性,本文針對(duì)無(wú)控旋轉(zhuǎn)穩(wěn)定彈對(duì)其進(jìn)行仿真分析,相關(guān)參數(shù)設(shè)計(jì)如下。①訓(xùn)練神經(jīng)網(wǎng)絡(luò)所需數(shù)據(jù)是由質(zhì)點(diǎn)彈道模型計(jì)算所得,初始彈道諸元設(shè)置為:初速v0=1 100 m/s,射角θa0=45°,射向角ψ20=0,積分步長(zhǎng)Δt=0.01 s;②BP神經(jīng)網(wǎng)絡(luò)調(diào)整參數(shù):隱層層數(shù)為1,隱層節(jié)點(diǎn)數(shù)為l=10,學(xué)習(xí)速率ηmax=0.7,ηmin=0.4,動(dòng)量因子α=0.1,權(quán)值wij、wjk和閾值a、b在[-1,1]間隨機(jī)初始化,迭代終止準(zhǔn)則一為最大迭代次數(shù)gmax=1 000,準(zhǔn)則二為誤差性能函數(shù)E≤0.001。

基于上述參數(shù)的設(shè)計(jì),首先分析BP神經(jīng)網(wǎng)絡(luò)對(duì)動(dòng)力學(xué)模型的逼近精度。本文基于質(zhì)點(diǎn)彈道模型共生成7 713組彈道數(shù)據(jù),隨機(jī)選擇其中5 000組數(shù)據(jù)來(lái)訓(xùn)練神經(jīng)網(wǎng)絡(luò),利用其余的數(shù)據(jù)來(lái)測(cè)試神經(jīng)網(wǎng)絡(luò)的逼近能力。質(zhì)點(diǎn)彈道的速度規(guī)律如圖2所示。

本文采用均方誤差(mean square error,MSE)及均方根誤差(root mean square error,RMSE)來(lái)反映神經(jīng)網(wǎng)絡(luò)的預(yù)測(cè)精度,表達(dá)式為

(14)

BP神經(jīng)網(wǎng)絡(luò)的MSE及RMSE值如表1所示。

表1 神經(jīng)網(wǎng)絡(luò)模型的MSE及RMSE值

由表1可得,本文所建立的神經(jīng)網(wǎng)絡(luò)模型預(yù)測(cè)的3個(gè)速度分量值、高度值與真實(shí)值基本一致,且RMSE值相比于狀態(tài)值而言非常小,這表明辨識(shí)效果較好。此外,神經(jīng)網(wǎng)絡(luò)的性能,即預(yù)測(cè)能力及收斂速度,則取決于網(wǎng)絡(luò)的調(diào)整參數(shù)。經(jīng)過(guò)調(diào)試分析可知,隱層層數(shù)、隱層節(jié)點(diǎn)數(shù)、最大迭代次數(shù)主要決定了BP神經(jīng)網(wǎng)絡(luò)的預(yù)測(cè)精度,而初始權(quán)值及閥值、學(xué)習(xí)速率和動(dòng)量因子主要影響B(tài)P神經(jīng)網(wǎng)絡(luò)的收斂速度。

網(wǎng)絡(luò)訓(xùn)練成功后,便可采用NGN法提取零升阻力系數(shù)。本文選擇7 500組仿真雷達(dá)測(cè)速數(shù)據(jù)(其是由質(zhì)點(diǎn)彈道模型與高斯白噪聲疊加而成)。為了滿足每一彈道段內(nèi)阻力系數(shù)可以近似為常數(shù)這一條件,本文共設(shè)置500個(gè)彈道段,每個(gè)彈道段內(nèi)最大速度降約為20 m/s。選定雷達(dá)量測(cè)噪聲為σvx=σvy=σvz=0,辨識(shí)結(jié)果如圖3所示。圖中,Cd為阻力系數(shù),Ma為馬赫數(shù)。

由圖2可得,基于NGN法辨識(shí)的阻力系數(shù)與基準(zhǔn)阻力系數(shù)基本一致,辨識(shí)精度較高,可滿足要求。相對(duì)而言,在跨聲速段的辨識(shí)誤差較大,這是分段小區(qū)間常數(shù)法的固有缺陷。經(jīng)過(guò)調(diào)試可知,神經(jīng)網(wǎng)絡(luò)的預(yù)測(cè)精度決定了NGN算法的辨識(shí)精度,而靈敏度參數(shù)的設(shè)置決定了NGN算法的收斂速度。

3.2 實(shí)測(cè)數(shù)據(jù)處理

接下來(lái)分析NGN法對(duì)雷達(dá)實(shí)測(cè)飛行數(shù)據(jù)的辨識(shí)效果。本文所用雷達(dá)的測(cè)量誤差:σγ=5 m,σβ=0.001 rad,σε=0.001 rad,σv=2 m/s。首先直接采用雷達(dá)量測(cè)數(shù)據(jù)訓(xùn)練BP神經(jīng)網(wǎng)絡(luò),訓(xùn)練并調(diào)試成功后,其MSE及RMSE值如表2所示。

表2 神經(jīng)網(wǎng)絡(luò)模型的MSE及RMSE值

將表1與表2進(jìn)行對(duì)比可得,BP網(wǎng)絡(luò)基于實(shí)測(cè)數(shù)據(jù)的預(yù)測(cè)精度較低,其原因可能如下:一是神經(jīng)網(wǎng)絡(luò)參數(shù)設(shè)置仍不是最優(yōu);二是由于訓(xùn)練數(shù)據(jù)是雷達(dá)實(shí)測(cè)數(shù)據(jù),結(jié)合六自由度剛體彈道模型可得,此時(shí)的狀態(tài)量還受到其他參量的影響,如速度狀態(tài)還受到攻角、彈道傾角、升力系數(shù)導(dǎo)數(shù)等的影響,而這些參量由于傳感器的原因,并沒(méi)有被作為輸入輸出量進(jìn)行考慮;三是訓(xùn)練網(wǎng)絡(luò)使用的基準(zhǔn)阻力系數(shù)并不是該發(fā)試驗(yàn)彈丸最精確的系數(shù)。

NGN法辨識(shí)的零升阻力系數(shù)Cx0與馬赫數(shù)Ma的關(guān)系如圖4所示。

由圖4可得,NGN法辨識(shí)的阻力系數(shù)與基準(zhǔn)阻力系數(shù)在馬赫數(shù)為[2,2.8]之間時(shí),曲線基本一致,偏差較小;在馬赫數(shù)為[1.2,2]之間時(shí),辨識(shí)曲線已難以跟隨基準(zhǔn)曲線,誤差較大;在馬赫數(shù)為[1,1.2]之間時(shí),辨識(shí)曲線圍繞基準(zhǔn)曲線振蕩不斷。出現(xiàn)上述現(xiàn)象的原因可能是:一是神經(jīng)網(wǎng)絡(luò)辨識(shí)的模型精度不高,這一點(diǎn)已在前文進(jìn)行了分析;二是分段小區(qū)間常數(shù)法的固有缺陷使然;三是雷達(dá)量測(cè)數(shù)據(jù)存在野值,可考慮針對(duì)濾波數(shù)據(jù)進(jìn)行訓(xùn)練和辨識(shí)。

為了進(jìn)一步分析NGN法的應(yīng)用特點(diǎn),本文又采用目前工程應(yīng)用已成熟的擴(kuò)展卡爾曼濾波(extended Kalman filter,EKF)方法處理該發(fā)彈丸的實(shí)測(cè)數(shù)據(jù),將分段辨識(shí)所得的阻力符合系數(shù)與對(duì)應(yīng)彈道段內(nèi)的基準(zhǔn)系數(shù)相乘,便得到全彈道的阻力系數(shù)圖,如圖5所示。

由圖5可得,EKF前期辨識(shí)精度較差,其原因是EKF法需要一定的數(shù)據(jù)量才能收斂于穩(wěn)定值。

綜上所述,在實(shí)際應(yīng)用中可將NGN和EKF這2種方法進(jìn)行結(jié)合,取長(zhǎng)補(bǔ)短,從而達(dá)到工程要求?；谏鲜鏊枷?本文在馬赫數(shù)為[2.2,2.8]時(shí),采用NGN法辨識(shí)的阻力系數(shù);在馬赫數(shù)為[1,2.2]時(shí),采用EKF辨識(shí)的結(jié)果。然后對(duì)雷達(dá)測(cè)量彈道、基于理論基準(zhǔn)阻力系數(shù)計(jì)算的彈道、基于辨識(shí)阻力系數(shù)計(jì)算的彈道三者的速度進(jìn)行比較,如圖6所示。

理論參數(shù)計(jì)算彈道、辨識(shí)參數(shù)計(jì)算彈道與雷達(dá)測(cè)量彈道的速度及位置的偏差均值如表3所示。

表3 狀態(tài)偏差均值表

由圖6及表3可得,分段辨識(shí)參數(shù)計(jì)算的彈道相比理論基準(zhǔn)參數(shù)計(jì)算的彈道而言,其與雷達(dá)測(cè)量彈道的接近程度更高。因此,可認(rèn)為分段辨識(shí)所得的阻力系數(shù)比基準(zhǔn)阻力系數(shù)更為準(zhǔn)確,從而說(shuō)明結(jié)合使用NGN法與EKF法具有一定的應(yīng)用意義。

4 結(jié)束語(yǔ)

本文首次將基于神經(jīng)網(wǎng)絡(luò)的參數(shù)辨識(shí)方法應(yīng)用于彈丸的氣動(dòng)參數(shù)辨識(shí)。通過(guò)對(duì)仿真結(jié)果和實(shí)際辨識(shí)結(jié)果進(jìn)行分析可得,當(dāng)輸入輸出量測(cè)數(shù)據(jù)較為完整且網(wǎng)絡(luò)的調(diào)整參數(shù)設(shè)置恰當(dāng)時(shí),NGN法能夠以較高的精度及可靠性辨識(shí)出阻力系數(shù);但當(dāng)所獲量測(cè)狀態(tài)數(shù)據(jù)并不完整時(shí),NGN法可用于補(bǔ)償EKF法,提高其前期的辨識(shí)精度,從而滿足工程應(yīng)用要求。本文研究結(jié)果可為神經(jīng)網(wǎng)絡(luò)在彈丸氣動(dòng)參數(shù)辨識(shí)中的實(shí)際應(yīng)用提供一定的理論參考。