朱豪坤,魚小軍,羅艷偉,肖 斌,任延超
(湖南云箭集團有限公司,湖南 長沙 410100)
現代戰機作戰性能不斷增強,作戰空域不斷提高,對制導航空炸彈的投放高度和飛行高度提出了更高的要求。隨著高度的增加,大氣變得更加稀薄,炸彈氣動數據相對其在低空飛行時有所變化。通常可通過CFD數值計算和地面風洞試驗兩種方法提供全面的炸彈氣動數據。然而,CFD數值計算方法在建模準確度、計算精度方面能力有限,風洞試驗難以復現炸彈在高空動態飛行時的實際條件,所得炸彈氣動數據與實際高空飛行時有一定差異,而通過參數辨識的方法從真實飛行數據中可辨識出更準確的炸彈氣動參數[1-3]。
目前研究中,歐陽光等[4]以最小二乘法為辨識準則,通過分布線性搜索相結合的方法,基于飛行數據對常規布局飛機的氣動參數進行了辨識;李寒冰等[5]以最小二乘法為準則,通過逐步回歸的方法對小型飛翼式無人機的相關氣動參數進行了辨識;Pinkelman等[6]利用最小二乘準則顯著減小參數估計中偏差誤差,對線性動態結構系統的時間序列模型進行了參數辨識。Ribeiro等[7]利用時域氣動彈性系統最小二乘法辨識出飛行器的氣動彈性參數,從而建立了氣動彈性模型。Grauer等[8]利用遞歸最小二乘和時域數據確定穩定導數和控制導數的精確不確定度方法,提高了辨識的可信度。此外Dutra[9]利用了輸出誤差法對飛機的氣動參數進行了辨識;Gandhi等[10]分別利用最小二乘回歸、非線性回歸、高斯回歸等多種方法對垂直起降飛行器的氣動參數進行了辨識。
文中參考已有研究成果,以制導航空炸彈為研究對象,針對其在高空稀薄大氣飛行時氣動模型不準的問題,基于方程誤差最小平方準則,利用回歸方法對其氣動參數進行辨識。通過仿真算例驗證了該方法的可行性。
方程誤差法進行參數辨識的原理如圖1所示,主要包括彈體運動、動力學模型、最小二乘辨識參數算法3個部分。涉及到的信號主要有:輸入信號u、測量噪聲v、觀測狀態x、觀測方程z。其中輸入信號u用來激活彈體短周期的運動模態,使飛行數據具備可辨識性。觀測狀態x通常為彈體輸出的線速度、角速度、角加速度、攻角、側滑角等飛行狀態量。觀測方程z為力與力矩系數方程,通常不可直接獲取,需要通過動力學對飛行數據處理得到。

圖1 方程誤差法原理圖Fig.1 Schematic diagram of equation error
所謂方程誤差法就是利用動力學模型將彈體實測飛行數據轉化為飛行過程中的力與力矩系數觀測方程z,同時通過既定回歸變量ξj和待辨識參數構造模型回歸方程y,以觀測方程和回歸方程的誤差作為優化目標的方法。以其誤差平方和最小為目標,通過線性回歸的手段對未知參數θ進行辨識。回歸方程通常事先已知,其一般形式可表示為:
(1)
式中:θj為np=n+1個待辨識參數,j=1,2,…,n;i=1,2,…,N,N為數據觀測的個數。將式(1)用向量的形式表示:
y=Xθ
(2)
式中:θ=[θ0,θ1,…,θn]T;X=[1,ξ1,…,ξn],為N×np的回歸矩陣。
最小二乘法的應用,通常要求系統滿足最小二乘的模型假設:
1)待辨識參數θ是一個未知常數參數的向量;
2)測量噪聲v是一個隨機向量,通常假設為零均值、不相關、方程不變的白噪聲。
最小二乘意義上的最優估計來自于最小化測量值與模型回歸方程輸出值的差的平方和所對應的參數。其指標函數可表示為:
(3)
則最小化指標函數的參數估計必須滿足:

(4)
式(4)表示的np個方程被稱為正則方程,求解正則方程可得到最小二乘估計器:
(5)
矩陣XTX總是對稱的。如果構成X的列的回歸向量是線性獨立的,那么XTX是正定的,XTX的特征值為正實數,相關的特征向量相互正交,對應的逆矩陣(XTX)存在。從式(4)中可以看出,代價函數相對于參數向量的二階梯度是XTX,它是正定的,表示該點的代價函數是一個極小值而不是極大值。
為保證炸彈飛行數據具備高可辨識性,需要對飛行試驗進行設計,試驗設計主要包括激勵輸入信號的設計和飛行數據采集設計。
輸入激勵信號設計是為了保證飛行源數據的可辨識性,因此要求輸入激勵信號能夠充分激勵彈體的短周期運動模態,使氣動數據中的靜態導數、控制導數、穩定導數充分解耦。彈體固有頻率可通過CFD或風洞所得數據采用式(6)計算[11]。文中選用“3211”輸入結構作為控制系統輸入,激勵彈體相關模態,其輸入模型結構如式(7)所示。其中td為機動動作開始的時刻;Ad為輸入的幅值,幅值的設計以可以充分激勵彈體運動動態特性,又使彈體處于可控制范圍之內為準則;ΔT為3211中“1”的時間長度,其選取原則是將“2”的脈沖寬度定為彈體固有模態頻率一半的時間段。
(6)
(7)
飛行數據采集設計是為了獲得足夠充分、精確的飛行狀態信息,作為辨識算法的輸入。在基于最小二乘準則的方程誤差法辨識中,辨識所需采集的數據包括但不限于:線加速度、角加速度、角速度、姿態角、動壓、氣動角。在實際工程中,線加速度和角速度可通過加速度計和陀螺儀測量;角加速度和姿態角分別通過角速度微分和積分獲得;相對風速、動壓和氣動角可通過彈載嵌入式大氣數據傳感系統(FADS)測量。
在通過線性回歸方法進行氣動參數辨識中,采用無量綱氣動力和力矩系數作為線性回歸的因變量。對每個力或力矩系數都當成一個單獨的線性回歸問題來求解,對應于最小化炸彈六個自由度下每個運動方程的方程誤差。由于氣動力和力矩系數不能在飛行中直接測量,必須基于其他測量量通過動力學方程進行轉換計算。必要的方程式如式(8)所示,對應圖1中的觀測方程z。
(8)
辨識用氣動模型回歸方程如式(9)所示:
(9)

基于最小二乘估計器的參數辨識算法由式(5)表示,下面以阻力系數cx的辨識過程舉例進行詳細說明:
(10)
其待辨識參數向量、觀測向量、狀態矩陣如式(10)所示。同樣,根據式(9)和式(8)進行類似處理,可得cy,cz,mx,my,mz的辨識算法。


表1 炸彈物性參數Table 1 Physical parameter of bomb
(11)
炸彈初始投放高度為15 000 m,初始速度為300 m/s,初始姿態角、姿態角速度、氣動角和線加速度均為0。炸彈離機飛行穩定后,高度約為14 500 m,在開環控制下在縱向平面做“3211”機動動作,幅值Ad=1°,單位間隔ΔT為由式(7)確定的自然震蕩頻率關系,選取仿真步長為0.02 s。
以飛行數據加上信噪比為5%信號強度的高斯白噪聲,生成帶噪聲的炸彈縱向通道飛行數據。經過辨識算法所辨識的結果如表2所示。

表2 辨識結果Table 2 Identificational results
利用表2中的辨識結果,可以得到縱向氣動參數的辨識預測輸出,由辨識數據重構飛行彈道與原觀測飛行數據對比如圖2所示,由辨識預測輸出與六自由度仿真觀測輸出的對比如圖3所示。

圖2 觀測飛行數據與重構飛行數據Fig.2 Observed and reconstructed flight data

圖3 觀測數據與辨識預測Fig.3 Observational data and identificatinal prediction data

針對制導航空炸彈氣動模型在高空稀薄環境不準的問題,根據方程誤差法,首先將觀測數據通過動力學方程轉換為氣動參數觀測方程輸出,再通過已知的氣動模型結構構造回歸方程,進而以方程輸出誤差平方最小為準則推導參數尋優算法,最后基于高空縱向通道機動飛行仿真算例對算法進行了驗證。仿真結果表明,該方法能夠較好的辨識高空環境制導炸彈的氣動參數,可為制導炸彈氣動模型的確定提供參考,且方法簡單易于實現,工程可實現性強。