劉洋， 常思江， 魏偉

(1.南京理工大學 能源與動力工程學院, 江蘇 南京 210094; 2.瞬態(tài)沖擊技術(shù)重點實驗室, 北京 102202;3.北京理工大學 機電學院, 北京 100081)

0 引言

氣動辨識是飛行器系統(tǒng)辨識中的關(guān)鍵問題，是炮彈、火箭彈、導彈、飛機等獲取自身氣動力參數(shù)的重要手段，在各類彈箭的氣動外形設(shè)計、系統(tǒng)參數(shù)確定、設(shè)計定型、射表編制等過程中具有不可替代的作用。氣動辨識是以一定的彈箭運動模型和辨識算法為基礎(chǔ)，從野外靶場或靶道實測彈箭運動數(shù)據(jù)中提取出彈箭氣動力參數(shù)的過程。理論上講，氣動辨識本質(zhì)上就是利用輸入、輸出信息確定系統(tǒng)結(jié)構(gòu)參數(shù)。從工程應(yīng)用角度看，與導彈、飛機等有控飛行器相比，常規(guī)炮彈、火箭彈等無控飛行器的氣動辨識難度更大，其原因主要在于：1)描述無控飛行器的動力學模型往往具有較強的非線性結(jié)構(gòu)[1]，數(shù)學處理難度較大；2)常規(guī)炮彈、火箭彈發(fā)射時的瞬時狀態(tài)往往具有隨機性且難以準確測量，并且由于穩(wěn)定飛行過程中不受控，無法為系統(tǒng)提供有效的輸入；3)某些彈丸(如槍彈)的氣動力參數(shù)非線性特性明顯，對彈道影響較大，而非線性氣動力研究本身就具有很大難度；4)常規(guī)炮彈、火箭彈等開展飛行試驗時，由于體積空間、試驗成本、發(fā)射條件、飛行環(huán)境等各方面的限制，往往造成測量數(shù)據(jù)的種類、精度及采樣頻率等，難以同導彈、飛機等相提并論，從而給氣動辨識帶來困難。

由于常規(guī)炮彈、火箭彈等無控飛行器的射程、密集度、穩(wěn)定性等與氣動力參數(shù)密切相關(guān)，盡可能準確地獲取氣動力參數(shù)，對性能評估、方案改進優(yōu)化等具有重要作用，故有必要深入開展氣動辨識方法與技術(shù)研究。從現(xiàn)有文獻看，相關(guān)氣動辨識研究的重點主要是辨識算法，例如Chapman-Kirk方法(C-K法)[2]、Levenberg-Marquardt方法(L-M法)[3-5]、極大似然法[6-9]以及某些全局優(yōu)化算法[3,10-11]等。實際工程中，為了確定彈箭氣動力參數(shù)，一般都是在不同初速、射角等條件下，成組發(fā)射彈箭并進行測量。對于1組射彈，由于氣象變化、起始擾動、彈箭結(jié)構(gòu)參數(shù)等具有不同程度的隨機性，即便測量設(shè)備具有相同誤差，同一組試驗的每發(fā)彈也不可能測得完全相同的彈道數(shù)據(jù)。已有實踐表明，同一組試驗中每發(fā)彈的氣動力參數(shù)辨識結(jié)果往往有一定差異，對于某些動態(tài)氣動力參數(shù)，如俯仰阻尼力矩系數(shù)、馬格努斯力矩系數(shù)等，發(fā)與發(fā)之間的差異甚至很大。目前，絕大多數(shù)相關(guān)研究關(guān)注的都是單發(fā)彈丸的氣動辨識過程，主要涉及模型和辨識算法，并未考慮如何優(yōu)化處理多發(fā)彈丸辨識結(jié)果的差異。例如，在對氣動力參數(shù)精度要求甚高的射表試驗數(shù)據(jù)處理中[12-13]，就是對多發(fā)彈測量數(shù)據(jù)辨識出的氣動力參數(shù)取平均值，作為該彈的氣動力參數(shù)以供后續(xù)使用。然而實際科研中發(fā)現(xiàn)，采用平均值氣動力參數(shù)去重構(gòu)彈道，經(jīng)常出現(xiàn)計算值與單發(fā)測量值差異較大的現(xiàn)象。這表明，盡管基于單發(fā)測量數(shù)據(jù)的獨立辨識結(jié)果具有最優(yōu)性，但平均值氣動參數(shù)對于多發(fā)測量數(shù)據(jù)來說并不是全局最優(yōu)的。目前常用的C-K法、L-M法、極大似然法等都屬于局部優(yōu)化方法，在用于單發(fā)測量數(shù)據(jù)獨立辨識時具有精度較高、收斂性好等優(yōu)點，但用于多發(fā)彈氣動參數(shù)聯(lián)合辨識卻具有較大局限性，辨識過程極不穩(wěn)定。近年來全局優(yōu)化算法在氣動辨識中的應(yīng)用研究逐漸增多，但大多還是解決單發(fā)數(shù)據(jù)氣動辨識的狀態(tài)初值問題、收斂性問題等[3,10]。在考慮多發(fā)彈氣動參數(shù)聯(lián)合辨識方面，國外有一些初步研究。如文獻[3]針對仿真數(shù)據(jù)，利用L-M法和遺傳算法分別對多發(fā)彈的氣動參數(shù)進行了聯(lián)合辨識，然而，L-M法收斂于局部極值，而遺傳算法的結(jié)果殘差比L-M法更大，這說明兩種算法都未能在多發(fā)彈氣動參數(shù)聯(lián)合辨識中取得較好的全局最優(yōu)解。文獻[4]中給出了L-M法和有限差分法的氣動辨識結(jié)果，比文獻[3]有較大改進，但也未能解決局部極值問題。文獻[14]提出采用多組試驗數(shù)據(jù)聯(lián)合辨識某無控空間探測器的氣動參數(shù)，但其重點在于將氣動參數(shù)表示成攻角和馬赫數(shù)的非線性函數(shù)，未考慮可能存在的局部解問題，并未給出具體的全局策略和算法。

本文針對槍彈、常規(guī)炮彈、火箭彈等無控彈箭的氣動參數(shù)聯(lián)合辨識問題開展研究，提出一個可同時處理多發(fā)彈測量數(shù)據(jù)并給出唯一具體氣動力參數(shù)值的全局優(yōu)化策略。以單發(fā)彈獨立辨識為基礎(chǔ)、以目前國際優(yōu)化領(lǐng)域新興的差分進化算法為工具，建立了利用多發(fā)彈測量數(shù)據(jù)聯(lián)合辨識氣動力參數(shù)的全局優(yōu)化流程，并通過對某大口徑炮彈攻角紙靶試驗數(shù)據(jù)的處理，驗證了本文所提策略的正確性和有效性。

1 多發(fā)彈測量數(shù)據(jù)聯(lián)合辨識策略

1.1 辨識流程設(shè)計

彈丸氣動辨識流程主要由以下5個部分組成[5]：測量數(shù)據(jù)預處理、數(shù)學模型、辨識算法、準則函數(shù)和辨識結(jié)果后處理。測量數(shù)據(jù)預處理和辨識結(jié)果后處理主要為剔除野值與曲線的平滑濾波，數(shù)學模型則根據(jù)辨識問題的不同可采用4自由度、5自由度、6自由度外彈道方程或攻角運動方程等，準則函數(shù)的選取與辨識算法有關(guān)。因此，多發(fā)彈測量數(shù)據(jù)聯(lián)合辨識的流程設(shè)計重點是辨識算法。

辨識算法按尋優(yōu)點個數(shù)可分為局部優(yōu)化算法和全局優(yōu)化算法，按搜索方向可分為梯度法和非梯度法，在氣動辨識中，一般選用基于梯度的局部優(yōu)化算法或非梯度的全局優(yōu)化算法。

文獻[4]中選取L-M法作為多發(fā)彈測量數(shù)據(jù)聯(lián)合辨識的算法，該算法是局部優(yōu)化算法中的一種，具有計算時間較短、計算效率較高的優(yōu)點，在對單發(fā)彈的辨識中得到了廣泛應(yīng)用。但是，利用局部優(yōu)化算法對多發(fā)彈測量數(shù)據(jù)進行辨識，存在以下2個缺點：1)多發(fā)彈聯(lián)合辨識中求得的矩陣可能達到幾十甚至上百階[3]，若矩陣接近奇異矩陣，計算可能發(fā)散；2)局部優(yōu)化算法尋得全局最優(yōu)的概率較低[3-4,15]，當預先給出的待辨識參數(shù)與實際值相差較大時，容易陷入局部最優(yōu)解。此外，由于實測數(shù)據(jù)的誤差種類較多，且有些誤差難以定量甚至難以定性(如紙靶試驗中彈丸穿靶瞬間對其自身的干擾)，可能使局部最優(yōu)解個數(shù)增多，從而導致局部優(yōu)化算法搜索到全局最優(yōu)解的概率更低。

因此，局部優(yōu)化算法不適用于多發(fā)彈測量數(shù)據(jù)聯(lián)合辨識。故考慮使用全局優(yōu)化算法作為多發(fā)彈測量數(shù)據(jù)聯(lián)合辨識的算法，該類算法計算穩(wěn)定性更好，搜索到全局最優(yōu)解的概率更高，但搜索維度較多，計算效率較低。

為改善全局優(yōu)化算法的性能，本文提出一種新的辨識流程：先利用局部優(yōu)化算法分別對各發(fā)彈進行快速辨識(以下簡稱獨立辨識)，利用辨識結(jié)果生成全局優(yōu)化算法中的搜索空間；再使用全局優(yōu)化算法對多發(fā)彈進行聯(lián)合辨識(以下簡稱聯(lián)合辨識)。具體辨識流程如圖1所示。圖1中n表示試驗彈丸的發(fā)數(shù)，1發(fā)彈對應(yīng)1套試驗數(shù)據(jù)、氣象數(shù)據(jù)和結(jié)構(gòu)參數(shù)，n發(fā)彈共有n套數(shù)據(jù)。待辨識參數(shù)為狀態(tài)初值和氣動參數(shù)。狀態(tài)初值為彈道方程中的積分初值或攻角方程的初始幅值、初始相位等，由于每發(fā)彈的發(fā)射條件、起始擾動等都有可能不同，導致每發(fā)彈的狀態(tài)初值差異較大，而某些狀態(tài)初值，如彈丸角運動的初始幅值與初始相位等難以直接測量，故狀態(tài)初值應(yīng)作為待辨識參數(shù)與氣動參數(shù)同時進行辨識。每組狀態(tài)初值中包含N1個待辨識參數(shù)，每組氣動參數(shù)中包含N2種待辨識氣動參數(shù)，如阻力系數(shù)、升力系數(shù)導數(shù)等。

圖1 辨識流程Fig.1 Identification procedure

圖1中的獨立辨識中僅對單發(fā)彈測量數(shù)據(jù)辨識，可根據(jù)實際情況選用C-K法、極大似然法等辨識算法，辨識結(jié)果為聯(lián)合辨識提供搜索空間，該搜索空間用于聯(lián)合辨識中待辨識參數(shù)的初始化，使其位于較小的約束范圍內(nèi)，進而減少迭代次數(shù)，提升計算效率與穩(wěn)定性；聯(lián)合辨識中對多發(fā)彈測量數(shù)據(jù)同時辨識，辨識算法選用全局優(yōu)化算法。

因此，多發(fā)彈測量數(shù)據(jù)聯(lián)合辨識的具體步驟為：

1)對多發(fā)彈測量數(shù)據(jù)進行預處理，剔除野值、平滑測量曲線等，并根據(jù)測量數(shù)據(jù)種類及待辨識氣動參數(shù)構(gòu)建數(shù)學模型；

2)使用局部優(yōu)化算法進行獨立辨識，得到n組狀態(tài)初值與n組氣動參數(shù)；

3)利用獨立辨識結(jié)果生成聯(lián)合辨識的搜索空間，使用全局優(yōu)化算法對多發(fā)彈相關(guān)氣動參數(shù)進行聯(lián)合辨識，得到n組狀態(tài)初值與1組氣動參數(shù)；

4)對辨識結(jié)果進行平滑、插值等后處理。

1.2 準則函數(shù)與最優(yōu)解定義

在整個辨識過程中，獨立辨識與聯(lián)合辨識的準則函數(shù)應(yīng)在形式上保持一致，如：不能同時使用最小二乘準則和極大似然準則。本文使用最小二乘準則作為辨識準則。

最小二乘準則可用(1)式表示：

(1)

在辨識、優(yōu)化等問題中，J值越小，表明全局優(yōu)化性越好[16]。在約束域內(nèi)：若J是最小值，則對應(yīng)的解是全局最優(yōu)解；若J是極小值且非最小值，則對應(yīng)的解是局部最優(yōu)解。將此概念擴展到多發(fā)彈測量數(shù)據(jù)聯(lián)合辨識的問題中，利用辨識出的1套氣動參數(shù)和n套狀態(tài)初值重構(gòu)出n條彈道，若對應(yīng)的J值是約束域內(nèi)的最小值，則該氣動參數(shù)與狀態(tài)初值共同構(gòu)成全局最優(yōu)解，反之則為局部最優(yōu)解。以上概念也作為后續(xù)評價全局性與最優(yōu)解的依據(jù)。

2 差分進化算法的應(yīng)用

2.1 基本算法

差分進化算法是近年來國際優(yōu)化領(lǐng)域新興的一種全局優(yōu)化算法，相比于遺傳算法、粒子群算法等，差分進化算法具有更好的收斂速度、魯棒性和全局尋優(yōu)能力[16]，適合計算量較大的聯(lián)合辨識。

差分進化算法主要是由“變異”、“交叉”和“選擇”3個部分構(gòu)成。但是在聯(lián)合辨識中，種群的搜索空間與初始化也較為重要。

差分進化算法的“變異”遵循(2)式：

(2)

式中：k為迭代次數(shù)；θ是待變異的個體，θr1和θr2是隨機的2個個體；a是[0,1]之間的1個隨機數(shù)。與粒子群算法等具有一定方向性算法不同，差分進化算法的變異方向是完全隨機的，這意味著算法具有更好的全局搜索能力。

“交叉”是指2個個體間隨機交換若干參數(shù)生成新的個體，也可以進行算數(shù)重組或者連續(xù)重組?！斑x擇”可以選用模擬退火的選擇方式，防止算法因為早熟而產(chǎn)生局部最優(yōu)解[17]。

2.2 算法應(yīng)用

在全局尋優(yōu)過程中，較大的搜索空間不僅增加迭代次數(shù)、降低辨識效率，對于某些彈道方程，如果氣動參數(shù)與狀態(tài)初值的迭代初值不在合理區(qū)間內(nèi)，可能導致計算發(fā)散。對于帶約束的優(yōu)化問題，若能將搜索空間預先確定在約束范圍內(nèi)，可大幅減少無效迭代。因此，在聯(lián)合辨識中，獲得合理搜索空間是非常重要的環(huán)節(jié)。

差分進化算法應(yīng)用于圖1中的聯(lián)合辨識部分，具體流程如圖2所示。

圖2 差分進化算法的應(yīng)用流程Fig.2 Flowchart of differential evolution algorithm

如圖2所示，聯(lián)合辨識的搜索空間由獨立辨識結(jié)果確定。由于每發(fā)彈狀態(tài)初值不同，故以狀態(tài)初值的辨識結(jié)果為中心按一定比例擴展生成新的搜索空間。經(jīng)反復調(diào)試發(fā)現(xiàn)擴展比例與狀態(tài)初值類型有關(guān)，當狀態(tài)初值為可測量(如速度、位置等)，擴展比例可取±2%～±5%；當狀態(tài)初值為不可測量(如角運動初始幅值和初始相位)，擴展比例可取±20%～±30%. 需說明的是，在調(diào)試過程中當采用±50%、±30%和±20%這3組不同擴展比例進行辨識，得到了相同的辨識結(jié)果，但迭代次數(shù)分別為 26 541次、9 538次和4 076次。這表明，對于本文所提策略，擴展比例大小對計算速度影響較大，但對計算精度影響較小。擴展比例根據(jù)具體問題調(diào)試選取，調(diào)試的一般原則是在保證多次測試均能成功獲得最優(yōu)解的前提下，盡可能地減小其數(shù)值[16]。氣動參數(shù)的搜索空間可由獨立辨識結(jié)果的最大值和最小值確定，為增大全局最優(yōu)解位于搜索空間的概率，該搜索空間可進行適當放大。因此，狀態(tài)初值的搜索空間共有n×N1維，氣動系數(shù)的搜索空間共有N2維。

由于差分進化算法具有一定的隨機性，為確保搜索到全局最優(yōu)解，當搜索空間維度較多時，應(yīng)增大個體數(shù)量N，使種群在初始化搜索空間內(nèi)的密度更大、分布更均勻，增大搜索到全局最優(yōu)解的概率[16]。此外，由于獨立辨識采用的是局部優(yōu)化算法，其給出的初始搜索空間有可能未包含全局最優(yōu)解，為解決這一問題，本文的策略是利用同一程序進行多次辨識，前幾次辨識不對搜索空間施加邊界約束，可使迭代向初始化搜索空間外進行。若前幾次辨識結(jié)果差異很小且位于搜索空間內(nèi)，該結(jié)果大概率為全局最優(yōu)解，則后續(xù)幾次辨識中可對搜索空間施加“出界反射”的邊界約束，即當“變異”后的個體在某個維度上超出搜索空間時，則在該維度上將其“反射”回搜索空間內(nèi)，以提高計算效率[16]；若前幾次辨識結(jié)果差異較大或不在搜索空間內(nèi)，后續(xù)若干次辨識仍然不對搜索空間施加邊界約束，以增加搜索到全局最優(yōu)解的概率。當多次辨識完成后，若辨識結(jié)果差異很小(如小于0.01%)，可認為辨識結(jié)果是全局最優(yōu)解；若辨識結(jié)果差異較大，應(yīng)選取物理意義正確且準則函數(shù)更小的解作為辨識結(jié)果，但該解是否為真正意義上的全局最優(yōu)解，還需結(jié)合氣動工程計算、計算流體力學數(shù)值計算等進行驗證。

辨識過程中，局部優(yōu)化算法和全局優(yōu)化算法的收斂條件均為

Jk

(3)

式中：Jk為準則函數(shù)的k次迭代值；b1=10-3. 然而，實際中往往達不到這樣的條件，所以可采用以下條件中的任意一種：

|Jk+1-Jk|

(4)

|θk+1-θk|

(5)

式中：Jk+1為準則函數(shù)的k+1次迭代值；b2=10-5，b3=10-7. (4)式代表準則函數(shù)已經(jīng)收斂，而(5)式代表待辨識參數(shù)已經(jīng)收斂。在局部優(yōu)化算法中，達到任意一種收斂條件都說明迭代已經(jīng)收斂，即使該解可能為局部最優(yōu)解。

對于差分進化算法，個體數(shù)量不止1個，所以無法使用(5)式，且由于可能存在局部最優(yōu)解，使得相差較大的θ能計算出很接近的J值，所以(4)式也不能單獨使用。在聯(lián)合辨識中，考慮到整個種群應(yīng)在所有維度上收斂，即有

(6)

式中：D()表示取方差；E()表示均值；i表示維度；b4=10-3.(6)式代表種群中所有個體在第i個維度達到收斂，當n×N1+N2個維度都滿足(6)式時，可認為種群在所有維度上收斂。當可能存在較多局部最優(yōu)解時，可以考慮同時使用多個種群進行辨識。當辨識結(jié)果滿足(5)式而不滿足(6)式時，說明辨識結(jié)果可能是局部最優(yōu)的；只有同時滿足(4)式和(6)式，才能說明辨識結(jié)果是全局最優(yōu)的。

綜上分析，在獨立辨識中，可以同時使用(3)式、(4)式和(5)式，滿足任意1個條件即認為計算收斂；在聯(lián)合辨識的差分進化算法中，需要同時滿足(4)式和(6)式，才能認為收斂到全局最優(yōu)解。b1、b2、b3和b4可視具體情況而定，基本原則是：當彈丸發(fā)數(shù)較少或測量點較少時，為保證充分收斂，b1、b2、b3和b4的取值應(yīng)更小(如比上述所給數(shù)值再小1個數(shù)量級)；當彈丸發(fā)數(shù)較多或測量點較多時，為優(yōu)化計算時間，b1、b2、b3和b4的取值可適當放大，但一般不能大于上述所給數(shù)值。

3 一個典型氣動力矩參數(shù)辨識問題

為了驗證本文建立的多發(fā)彈測量數(shù)據(jù)聯(lián)合辨識策略，本節(jié)考慮一個典型的氣動力矩參數(shù)辨識問題[17]，即利用紙靶試驗(測攻角)辨識彈丸靜力矩系數(shù)導數(shù)、俯仰阻尼力矩系數(shù)導數(shù)及馬格努斯力矩系數(shù)導數(shù)，將上述辨識流程應(yīng)用于實測紙靶數(shù)據(jù)處理。

3.1 氣動辨識用數(shù)學模型

由于紙靶試驗往往采取平射，單發(fā)彈丸飛行時間較短，可認為全彈道馬赫數(shù)基本不變、轉(zhuǎn)速無衰減；同組試驗的多發(fā)彈之間初速相差很小，認為具有相同馬赫數(shù)?？刹捎脧凸ソ沁\動方程作為氣動力矩參數(shù)辨識用數(shù)學模型。

彈丸復攻角運動方程的齊次形式[18]如下：

(7)

式中：Δ為復攻角；s為彈道弧長；i為虛數(shù)單位；H為角運動阻尼作用項；P為轉(zhuǎn)速作用項；M為靜力矩作用項；T為升力和馬格努斯力矩耦合作用項；α為縱向攻角；β為橫向攻角；Clα為彈丸升力系數(shù)導數(shù)；CD為阻力系數(shù)；ρ為大氣密度；d為彈徑；m為彈丸質(zhì)量；l為特征長度；Iy為赤道轉(zhuǎn)動慣量；Ix為極轉(zhuǎn)動慣量；p為彈丸轉(zhuǎn)速；v為彈丸運動速度；CMα為靜力矩系數(shù)導數(shù)；CMqα為俯仰阻尼力矩系數(shù)導數(shù)；CMpα為馬格努斯力矩系數(shù)導數(shù)。

該方程的解析解[18]為

(8)

式中：KS0和KF0分別為慢圓運動和快圓運動的初始幅值;φS0和φF0分別為慢圓運動和快圓運動的初始相位;λS和λF分別為慢圓運動和快圓運動的阻尼指數(shù);φ′S和φ′F分別為慢圓運動和快圓運動的角頻率。

旋轉(zhuǎn)穩(wěn)定彈的陀螺穩(wěn)定因子Sg為

(9)

工程上一般取Sg>1.3為陀螺穩(wěn)定性條件，只有滿足該條件，才能保證彈丸形成周期運動[18]。因此，陀螺穩(wěn)定性條件即可作為氣動辨識的約束條件：

P2-5.2M>0，

(10)

在整個辨識過程中，P和M的取值需始終滿足約束條件(10)式。

因此，計算對應(yīng)s處的α和β時，先利用(7)式計算出H、P、T和M，并檢驗其是否滿足約束條件(10)式：若滿足約束，則將其代入(8)式中進行計算，若不滿足約束，則對CMα、CMqα和CMpα重新取值，直到其滿足約束條件。

3.2 試驗條件及測量數(shù)據(jù)

在某大口徑榴彈的紙靶試驗中，共射擊1組6發(fā)彈丸。試驗前，對每發(fā)彈進行了靜態(tài)參數(shù)測量，且每發(fā)彈的紙靶布置方式相同：第1張紙靶布置在距離炮口25 m處，測量段s為25～125 m，每間隔5 m布置1張紙靶，彈丸平均飛行時間t≈0.11 s，平均馬赫數(shù)Ma=2.679 2，彈丸初速的最大值與最小值之差為3.1 m/s，符合3.1節(jié)中描述的假設(shè)。試驗完成后讀取紙靶數(shù)據(jù)并對其進行曲線平滑和野值剔除，處理后每發(fā)彈有15～17個測量點。每發(fā)彈的m、Ix、Iy、和v均采用實測值；由于沒有直接測量轉(zhuǎn)速p，故將p作為待辨識參數(shù)。彈道測量數(shù)據(jù)為α、β和s，辨識中p和v作為常數(shù)處理。

表1中給出了獨立辨識中的48個待辨識參數(shù)。

表1 獨立辨識的待辨識參數(shù)

而聯(lián)合辨識中的待辨識參數(shù)為

共有33個待辨識參數(shù)。

3.3 辨識的收斂性檢驗

利用數(shù)值仿真MATLAB軟件編程并對多發(fā)彈測量數(shù)據(jù)進行辨識，為充分地利用計算資源，使用軟件中的Parallel Pool模塊進行多線程并行計算。在獨立辨識中，使用不同線程對多發(fā)彈并行計算；在聯(lián)合辨識中，使用不同線程對多個個體并行計算。

在獨立辨識中，使用文獻[19]中方法進行辨識，6發(fā)彈各自的準則函數(shù)收斂曲線如圖3所示。

圖3 獨立辨識收斂性Fig.3 Convergence of independent identification

由圖3可以看出，6發(fā)彈的準則函數(shù)都不滿足(3)式，而是滿足(4)式或(5)式。由于預先給定的初始相位、幅值和氣動參數(shù)與實際值(真值)相差較大，而迭代增量由導數(shù)確定，故圖3中各條曲線的截距不同，而斜率相同。

搜索空間維度D=33，根據(jù)2.2節(jié)中的辨識策略，選取個體數(shù)量N=2 000，并使用同一程序連續(xù)進行10次辨識。前5次辨識沒有進行邊界約束，其辨識結(jié)果均位于初始化搜索空間內(nèi)且相對誤差小于0.01%，可認為該結(jié)果很可能為全局最優(yōu)解，后5次辨識施加“出界反射”的邊界約束以提高計算效率。10次辨識結(jié)果的相對誤差小于0.01%，故認為該辨識結(jié)果即為全局最優(yōu)解。

準則函數(shù)最大值Jmax、最小值Jmin、平均值Jmean和方差D(J)能夠反應(yīng)整個種群在迭代過程中的收斂情況。Jmin所對應(yīng)個體是潛在的最優(yōu)解，而Jmax與Jmean則能代表種群的初始分布情況與收斂速度。圖4給出了最后1次辨識中J和D(J)隨迭代次數(shù)變化的曲線。

圖4 聯(lián)合辨識收斂性Fig.4 Convergence of combined identification

由圖4可以看出：Jmax在迭代之初的值很大，但是下降速度非?？?，k=1 000時幾乎與Jmean和Jmin重合，且D(J)<10-3，準則函數(shù)基本收斂；k=2 000之后3條曲線就完全重合；k=3 000時，D(J)<10-5，整個種群的準則函數(shù)已經(jīng)完全收斂。計算完成后，Jmin對應(yīng)的氣動參數(shù)與狀態(tài)初值即為多發(fā)彈聯(lián)合辨識的全局最優(yōu)解。

除了J需要收斂外，種群中的所有個體θ也需收斂。種群位置的關(guān)系與辨識的收斂性相關(guān)，種群位置越集中，收斂性越好。圖5給出了不同k時種群在氣動參數(shù)搜索空間中的位置，為方便觀察種群位置變化，在圖5(b)～圖5(f)中將其投影到CMα、CMpα和CMqα相互正交形成的3個面上。

圖5 聯(lián)合辨識氣動參數(shù)收斂性Fig.5 Convergence of aerodynamic parameters for combined identification

從圖5(a)～圖5(f)可以看出各氣動參數(shù)的迭代過程：當k=1時，種群的初始分布較為均勻，能充分地填充整個搜索空間；當k=179時，種群已向某一區(qū)域收斂；當k=1 000時，相比于初始搜索空間，整個種群幾乎收斂為一個點；隨著k的增大，種群的收斂性越來越高，但是在k=2 000之后，收斂方向發(fā)生改變。當搜索空間內(nèi)僅有1個最優(yōu)解時，在達到一定迭代次數(shù)后，種群的收斂方向應(yīng)僅向內(nèi)部收縮，而種群位置幾乎不會移動[16]，如果種群位置發(fā)生移動，說明J的梯度方向發(fā)生改變，這種改變可能陷入局部最優(yōu)解。

從圖5(b)和圖5(c)的投影可以看出，從k=179到k=1 000的收斂過程即為整個種群向內(nèi)部收縮的過程；從圖5(d)～圖5(f)的投影可以看出，雖然從k=2 000到k=4 000同樣是整個種群向內(nèi)部收縮，但種群整體的位置也發(fā)生了改變，說明該位置附近可能存在局部最優(yōu)解。在多次辨識中均出現(xiàn)了這種情況，說明當多發(fā)彈同時辨識時，若使用局部優(yōu)化算法或搜索能力不足的全局優(yōu)化算法時，很容易陷入局部最優(yōu)解[3]。

在聯(lián)合辨識中，由于無需雅克比矩陣，且搜索空間的選取較為合理，辨識過程中始終滿足(10)式的約束，未出現(xiàn)計算發(fā)散；而獨立辨識中卻出現(xiàn)了發(fā)散，人為減小步長后才能收斂。這表明全局優(yōu)化算法的穩(wěn)定性要優(yōu)于局部優(yōu)化算法。

3.4 辨識結(jié)果與分析

為便于討論，以單發(fā)彈獨立辨識結(jié)果重構(gòu)的彈道稱為獨立辨識彈道，以聯(lián)合辨識結(jié)果重構(gòu)的彈道稱為聯(lián)合辨識彈道。理論上講，1種彈形對應(yīng)唯一的氣動參數(shù)，實際中由于多發(fā)彈辨識出的氣動參數(shù)在數(shù)值上并不完全一致，在對辨識結(jié)果的處理中，往往取氣動參數(shù)的平均值或者加權(quán)平均值，作為該彈的氣動參數(shù)[12-13]。將獨立辨識結(jié)果中的氣動參數(shù)取平均值，利用該氣動參數(shù)與辨識結(jié)果中的狀態(tài)初值重構(gòu)的彈道稱為平均值計算彈道，其準則函數(shù)計算方法與聯(lián)合辨識彈道相同。

由于單發(fā)彈的獨立辨識不具有全局性，故以平均值計算彈道和聯(lián)合辨識彈道為主要研究對象，獨立辨識彈道作為參考，Jt為單發(fā)彈辨識彈道準則函數(shù)，Jm為平均值計算彈道準則函數(shù)，Jc為聯(lián)合辨識彈道準則函數(shù)，獨立辨識結(jié)果與平均值計算結(jié)果如表2～表4所示。

表2 狀態(tài)初值的獨立辨識結(jié)果

表3 氣動參數(shù)的獨立辨識結(jié)果

表4 氣動參數(shù)平均值計算結(jié)果

從表3和表4中可以看出：CMα的均值與極差分別為3.428與0.46，極差與均值之比為13.4%，CMqα的極差與均值之比為121%，CMpα的極差與均值之比為130%；平均值計算彈道的準則函數(shù)Jm=4.996 5，遠大于獨立辨識彈道的準則函數(shù)之和Jt=1.685 1，這說明獨立辨識彈道的準則函數(shù)雖然較小，但多發(fā)彈辨識結(jié)果相差較大，故平均值氣動參數(shù)計算彈道與測量值相差較大。

利用獨立辨識結(jié)果生成初始搜索空間進行多發(fā)彈聯(lián)合辨識，經(jīng)反復調(diào)試，選取狀態(tài)初值和氣動參數(shù)的擴展比例均為±20%. 結(jié)果如表5和表6所示。

表5 狀態(tài)初值的聯(lián)合辨識結(jié)果

表6 氣動參數(shù)的聯(lián)合辨識結(jié)果

對比表3、表4和表6可以看出，盡管兩種方法辨識出的狀態(tài)初值和氣動參數(shù)相差不大，但聯(lián)合辨識彈道的準則函數(shù)值Jc=2.089 8卻遠小于平均值計算彈道的準則函數(shù)Jm. 需要說明的是，這6發(fā)彈獨立辨識的準則函數(shù)值并不相同，表3中給出的是6發(fā)彈的準則函數(shù)之和。為便于研究，這里選擇獨立辨識結(jié)果中準則函數(shù)最小的第4發(fā)彈作為研究對象，研究攻角變化規(guī)律與辨識誤差。

將第4發(fā)彈的獨立辨識彈道、聯(lián)合辨識彈道、平均值計算彈道與測量值進行對比，結(jié)果如圖6所示，圖6中δ為攻角。測量段s取25～125 m，為便于研究，2條辨識彈道和平均值計算彈道范圍取s為0～225 m，由于彈丸初速和轉(zhuǎn)速較高，該段彈道上的速度變化約為1%，轉(zhuǎn)速變化約為0.5%，工程上可作為常數(shù)處理。

圖6 辨識彈道與平均值計算彈道Fig.6 Identified and calculated trajectories

從圖6中可以看出：獨立辨識彈道和聯(lián)合辨識彈道與測量值較為接近，而平均值計算彈道與測量值相差較大；3條曲線在第1個攻角周期時相差不大，峰值位置和高度也很接近，隨著距離的增加，差異逐漸變大，這說明在待辨識參數(shù)中，KS0、KS0、φS0和φF0的誤差對彈道影響較小，而CMα、CMqα、CMpα和p的誤差對彈道影響很大。

由于測量值間隔較大，且在攻角峰值和零點處缺少測量值，而該發(fā)彈的獨立辨識彈道與測量值最接近，故以圖6中獨立辨識彈道為基準，研究第4發(fā)彈聯(lián)合辨識彈道與平均值計算彈道的運動規(guī)律。

在測量段內(nèi)，聯(lián)合辨識彈道每個周期的運動距離與獨立辨識彈道相差0.9%，前3個周期攻角峰高度值分別相差0.118°、0.016°和0.064°，峰值位置分別相差0.35 m、0.06 m和0.49 m；平均值計算彈道每個周期的運動距離與獨立辨識彈道相差2.3%，并從第2個周期開始，攻角峰值高度、位置開始出現(xiàn)偏離，這種偏離隨著距離增加，前3個周期攻角峰值高度分別相差0.066°、0.174°和0.256°，峰值位置分別相差0.73 m、1.78 m和2.83 m. 在測量段之后，由于誤差逐漸累積，在第5個運動周期時，平均值計算彈道對應(yīng)的峰值位置將提前4.93 m，對于本次試驗，這意味著將提前1張紙靶出現(xiàn)攻角峰值。平均值計算彈道第4、第5個峰值高度與獨立辨識彈道的峰值高度之差達到15%和19%，而聯(lián)合辨識彈道僅有6%和8.8%.

以攻角計算值δ與對應(yīng)彈道弧長處的攻角測量值δe之差作為絕對誤差Δδ，以絕對誤差與測量值之比Δδ/δe作為相對誤差，具體如表7所示。

表7 攻角誤差

從表7中可以看出，獨立辨識彈道的誤差最小，聯(lián)合辨識彈道誤差次之，平均值計算彈道誤差最大。相對誤差最大的點出現(xiàn)在δ=0°附近，因而該點的相對誤差較大，絕對誤差較小。獨立辨識彈道和聯(lián)合辨識彈道的最大絕對誤差分別為0.213°和0.209°，而平均值計算彈道的最大絕對誤差達到了0.631°，由于攻角最大測量值僅為2.9°，故該誤差不容忽視。攻角相對誤差和絕對誤差隨距離的變化如圖7所示。

圖7 攻角誤差Fig.7 Errors of angle of attack

從圖7中可以看出，獨立辨識彈道、聯(lián)合辨識彈道和平均值計算彈道相對誤差的2個峰值分別出現(xiàn)在s=45 m和s=90 m處，這兩處均在零點附近，除此之外，最大相對誤差分別為14.7%、14%和40.6%. 最大絕對誤差分別占攻角最大測量值的7.34%、7.21%和21.76%，結(jié)合表7中的誤差大小可知，平均值計算彈道的誤差約為其他2條彈道的2～3倍。

對于第4發(fā)彈丸，獨立辨識彈道和聯(lián)合辨識彈道與測量值的誤差較小，而平均值氣動參數(shù)計算的彈道誤差則相對較大。其他5發(fā)彈與第4發(fā)彈類似，只是誤差數(shù)值大小不同。對于全部彈丸，Jm比Jc大139.1%，因此，聯(lián)合辨識所得氣動參數(shù)更接近實際，全局性更好。

值得注意的是，在氣動辨識中，只有同時計算出CMα、CMqα、CMpα和轉(zhuǎn)速p，才能得到準確的彈丸角運動周期和阻尼指數(shù)。而KS0、KS0、φS0和φF0僅與試驗條件有關(guān)，每發(fā)彈丸都可能有所不同，故只有同時辨識出以上所有參數(shù)，才能得到全局最優(yōu)解。因此可以認為，只有聯(lián)合辨識的結(jié)果為全局最優(yōu)，平均值氣動參數(shù)不是全局最優(yōu)。

由于加工制造誤差及發(fā)射條件不完全相同，每發(fā)彈丸在實際飛行過程中的真實氣動特性不可能完全相同，且由于飛行環(huán)境差異、測量誤差等影響，每發(fā)彈丸的辨識結(jié)果僅是對該發(fā)彈實際條件下的最優(yōu)估計，不一定是真值，而本文對于多發(fā)彈聯(lián)合辨識結(jié)果，實際上是對這一型或這一批次彈丸氣動特性的最優(yōu)估計，對于未進行飛行試驗的同型彈丸來說，該辨識結(jié)果是對其氣動特性的最可靠預測。

4 結(jié)論

本文提出了一種利用全局優(yōu)化算法對多發(fā)彈測量數(shù)據(jù)進行聯(lián)合辨識的策略，設(shè)計了對應(yīng)的辨識流程。利用某紙靶試驗中多發(fā)彈測量數(shù)據(jù)對該辨識策略進行了驗證，并與現(xiàn)有的氣動辨識方法進行對比，得出以下主要結(jié)論：

1)應(yīng)用本文辨識策略不僅獲得了彈丸氣動參數(shù)與狀態(tài)初值的全局最優(yōu)解，還通過利用局部優(yōu)化算法構(gòu)建初始搜索空間，進一步提高了計算效率和穩(wěn)定性，有效地解決了多發(fā)彈測量數(shù)據(jù)聯(lián)合辨識氣動參數(shù)的局部最優(yōu)問題，所得辨識結(jié)果是對該型彈丸氣動特性的最可靠預測。

2)與現(xiàn)有方法相比，應(yīng)用本文辨識策略針對多發(fā)彈測量數(shù)據(jù)辨識，所得準則函數(shù)更小，重構(gòu)的彈道與測量值更接近，能更為準確反映彈丸運動規(guī)律。

3)差分進化算法應(yīng)用于多發(fā)彈測量數(shù)據(jù)聯(lián)合辨識氣動參數(shù)，具有較好的計算穩(wěn)定性與收斂性，且算法本身無需求導和大矩陣求逆，進一步增強了計算穩(wěn)定性，適于數(shù)據(jù)量較大的氣動辨識問題，具有較好的實用性和通用性，為多發(fā)彈聯(lián)合辨識氣動力參數(shù)提供了新的思路。