李寒冰, 吳大衛(wèi)

(中國商用飛機有限責任公司 上海飛機設計研究院, 上海 201210)

0 引言

研制高性能無人機需要建立其飛行動力學模型,它是無人機控制系統設計與調校、飛行仿真及飛行品質評價的重要前提和基礎,而獲得準確的各項氣動導數是成功建模的關鍵。

目前無人機的氣動導數主要用以下兩種方法求得:一種是理論計算法,采用傳統工程估算和現代計算流體力學技術,但在計算精度、建模準確程度以及求解時間上往往難以兼顧;另一種是試驗法,包括風洞試驗法和參數辨識試驗法。其中風洞試驗法的缺點在于試驗周期較長,成本較高。參數辨識試驗法是以飛行試驗測得的狀態(tài)數據作為飛行動力學參數化方程的輸入輸出數據,選擇合適的辨識方法在仿真環(huán)境下經過有限次迭代得到參數模型中未知參數的值。辨識結果能真實反映無人機的本身屬性,因此成為求取無人機氣動導數的一種有效手段[1-2]。

本文以某小型飛翼無人機為例,針對飛行試驗測量條件有限,提出了一種根據可測數據估算迎角、側滑角和氣動系數的方法,從而得到完整的飛行數據用于無人機氣動參數辨識中,并驗證了該方法的有效性。

1 飛行試驗限制及數據完善

1.1 小型無人機飛行試驗的特點

相比有人駕駛飛機,無人機采用了不同的控制系統和配套機載設備,在飛行試驗中表現出技術風險高、調控難度大、信息處理復雜、氣動參數測量困難等[3-4],而對于小型無人機而言又具有自身承載能力小的特殊性。

反映無人機系統動態(tài)特性的參數有:線位移、線位移速度和加速度、角位移、角速度和角加速度等。一般可利用大氣數據傳感器、慣性測量單元和無線定位天線等設備測量和采集記錄上述數據。但是,在缺少氣動角數據的情況下,還不足以對無人機的氣動參數進行辨識。氣動角包括迎角和側滑角,主要采用壓差歸零和風標式傳感器進行測量[5]。而對于小型無人機的重量和尺寸限制,這些測量裝置難以在保證精度的條件下實現同步小型化和輕量化,使小型無人機氣動參數辨識工作難以展開。

1.2 迎角和側滑角的估算方法

在現有試驗條件下,可以測得體軸系下三軸過載nx,ny,nz和空速V的準確值。動坐標系下的加速度方程[6]為:

dV/dt=δVB/δt+ωVB

(1)

方程(1)寫成體軸系標量形式,即:

(2)

本文采用四階龍格庫塔算法求解常微分方程組(2),可得到體軸系速度分量u,v,w關于過載nx,ny,nz和角速率p,q,r試驗數據的數值解。解得體軸系速度分量并與空速測量值V聯立可得迎角和側滑角:

(3)

上述方法彌補了試驗條件的不足,合理地估算了迎角和側滑角這兩個重要參數,完善了飛行數據,為無人機氣動參數辨識提供了可行性。

1.3 方法適用范圍

上述估算迎角和側滑角的方法是在無人機試驗條件有限、無法測得準確氣動角的情況下提出的一種完善飛行數據的簡化方法。試驗證明,在無風或微風的飛行條件下,該方法所得估算結果能夠很好地再現氣動角的響應歷程。另外,由于辨識模型使用了實際的飛行數據作為辨識輸入輸出數據,無需復雜的推導,具有通用性。因此,上述完善辨識試驗數據的方法對其他布局形式的固定翼無人機也是普遍適用的。本文還采用三軸過載試驗數據估算氣動系數,詳見下文。

2 研究對象和飛行試驗設計

2.1 研究對象

本文以某小型飛翼布局無人機為研究對象(見圖1),其基本參數如表1所示。

圖1 用于飛行試驗的飛翼無人機Fig.1 The flying wing UAV for flight test

機長L/m翼展b/m起飛質量WTO/kg有效載荷WPL/kg續(xù)航時間t/min1.21.761.545

飛機彈射起飛后采用人工對飛機姿態(tài)進行調整控制,完成配平后輸入激勵指令。飛行數據通過機載設備記錄并傳輸給地面站。采集的數據主要包括:副翼偏角δa、升降舵偏角δe、姿態(tài)角(滾轉角φ、俯仰角θ、偏航角ψ)、角速度(滾轉角速度p、俯仰角速度q、偏航角速度r)、過載(nx,ny,nz)、空速V、氣壓高度H、升降速率、GPS速度和高度、航向及經緯度等。

飛行試驗選擇在初冬進行,微風無持續(xù)風向,低空飛行可近似視為無風條件。

2.2 飛行試驗設計

飛行試驗中需要充分激發(fā)飛機的各種運動模態(tài),獲取足夠的信息量用于氣動參數辨識。

常用的激勵信號有階躍輸入、脈沖輸入、偶極方波輸入、3-2-1-1輸入等。對于本文的研究對象而言,階躍輸入可能會導致飛機向某一方向偏離,對飛行安全不利;脈沖輸入在實際應用中難以實現,而且由于舵機速率限制,過快的脈沖輸入往往不能激勵飛機的飛行模態(tài);3-2-1-1輸入信號的時間歷程遠遠超過了飛機的短周期時間歷程,因而無法激勵飛機的短周期模態(tài)。因此，本文采用偶極方波輸入,其數學模型為:

(4)

該信號頻譜密度較寬,在高頻區(qū)域也含有能量,適當調整Δt可以在感興趣的模態(tài)頻段含有足夠的能量,且克服了其他信號的使用限制,是個簡單而有效的輸入信號[2,7]。

3 縱向參數辨識

3.1 迎角計算結果

在飛機縱向氣動參數辨識試驗中,對升降舵進行偶極方波激勵,為了方便與正負過載對應，本文取拉桿時升降舵偏角為正值。

試驗可以測得體軸系下三軸過載nx,ny,nz和空速V的準確值,由式(1)～式(3)可以得到迎角α的估算值。圖2為升降舵偶極方波激勵及其對應的α的響應曲線估算結果。

圖2 升降舵激勵信號及迎角響應估算曲線Fig.2 Elevator input and response of angle of attack

3.2 逐步回歸法的基本原理

縱向參數辨識采用逐步回歸法,其基本思想是:將變量逐個引入,引入變量的條件是其偏回歸平方和經檢驗是顯著的;同時每引入一個新變量,對已選入的變量要進行逐個檢驗,將不顯著變量剔除。這樣保證最后所得的變量子集中的所有變量都是顯著的,經若干步可得“最優(yōu)”變量子集[8]。具體思想如下:

設有m個自變量xj(j=1,2,…,m),因變量為y,有n個觀測點。根據最小二乘原理,y的估計值為:

(5)

式中,l

3.3 辨識過程及結果

推力T沿體軸系xB軸,此處xB軸向后,zB軸向上。飛機的縱向力平衡方程[9]為:

(6)

(7)

利用飛行試驗數據,可得到相應的氣動系數值用于參數辨識。Cz的回歸方程可寫為[10]:

Cz=CzαΔα+CzqΔq+CzδeΔδe+Czα2Δα2+Cz0

(8)

待辨識參數有Czα,Czq,Czδe,Czα2,Cz0。

在圖2所示激勵下的部分狀態(tài)量動態(tài)響應曲線如圖3所示。

圖3 無人機的縱向動態(tài)響應Fig.3 Longitudinal dynamic response of UAV

采用逐步回歸法,對Cz的各項導數進行參數辨識,得到表2所示的辨識結果。

利用辨識所得參數計算得到Cz的辨識模型預測輸出,與其試驗曲線的對照如圖4所示。

圖4 垂向力系數辨識模型預測值與試驗值對比Fig.4 Identification result and flight test result of vertical force coefficient

可以看出,Cz的辨識模型預測輸出與其實際測量計算曲線能夠很好地吻合。經計算,二者擬合誤差為10.31%,辨識精度較為滿意。由于辨識過程是在無風假設基礎上進行的,大氣擾動可能會對辨識精度帶來一定誤差。

將辨識得到的各項導數值代入式(8)可得到各狀態(tài)量的辨識模型預測輸出。與飛行試驗測量曲線進行對照,可以用來檢驗辨識所得氣動參數的正確性。以俯仰角速度為例,其辨識模型預測輸出和實際飛行測量輸出對比如圖5所示。

圖5 俯仰角速度的辨識模型預測值與試驗值對比Fig.5 Identification result and flight test result of pitch rate

可以看出二者能夠很好地吻合。經檢驗,其他相關狀態(tài)量辨識曲線與試驗測量曲線均能較好地吻合,說明以上辨識結果是可信的。

4 橫側向參數辨識

4.1 側滑角計算結果

在飛機橫側向氣動參數辨識試驗中,對副翼進行偶極方波激勵,取右副翼下偏為正偏角。

同上文所述,可得到側滑角的估算值。圖6為副翼偶極方波激勵及其對應的側滑角響應曲線估算結果。

圖6 副翼激勵信號及側滑角響應估算曲線Fig.6 Aileron input and response of side slip angle

4.2 最小二乘法的基本原理

橫側向參數辨識采用最小二乘法。在系統辨識領域中,這是一種基本的估計方法。隨機環(huán)境下使用最小二乘法時,并不要求觀測數據提供其概率統計方面的信息,而其估計結果卻有相當好的統計特性[11]。最小二乘估計算法的數學形式可以簡要描述如下:

(9)

(10)

(11)

J為極小值的充分條件是:

(12)

即矩陣ΦTΦ為正定矩陣,或者說矩陣ΦTΦ是非奇異的。

4.3 辨識過程及結果

此處取歐美體軸系[6],由力矩平衡方程

(13)

可以得到基于飛行試驗數據的力矩系數計算公式。其中以偏航力矩系數Cn為例:

(14)

同理Cn的回歸方程可寫為[10]:

Cn=CnβΔβ+CnpΔp+CnrΔr+CnδaΔδa+Cn0

(15)

待辨識參數有Cnβ,Cnp,Cnr,Cnδa,Cn0。

在圖6所示激勵下的部分狀態(tài)量動態(tài)響應曲線如圖7所示。

圖7 無人機的橫側向動態(tài)響應Fig.7 Lateral dynamic response of UAV

用上述最小二乘法對Cn的氣動導數和操縱導數進行參數辨識,其結果如表3所示。

表3 最小二乘法辨識結果Table 3 Least square method identification results

利用表3中的辨識結果,可以得到Cn的辨識預測輸出,與實際試驗值的對照如圖8所示。

圖8 偏航力矩系數辨識模型預測值與試驗值對比Fig.8 Identification result and flight test result of yaw moment coefficient

從對比曲線可以看出,Cn的辨識曲線與其試驗曲線吻合較好,計算可知二者擬合誤差為14.36%。與前述Cz的辨識結果相比,Cn的辨識精度略低,其原因在于飛機橫側向運動比縱向復雜,存在橫航向耦合。各狀態(tài)量的測量誤差和大氣擾動都會對辨識精度產生一定的影響。

與縱向參數辨識結果檢驗方法類同,將所得Cn的導數值代入式(15),得到各相關狀態(tài)量的辨識模型預測輸出,與其實際飛行試驗測量輸出進行對照即可檢驗辨識結果的正確性。以偏航角速度r為例,其辨識模型預測輸出和實際飛行測量值對比如圖9所示。

圖9 偏航角速率辨識模型預測值與試驗值對比Fig.9 Identification result and flight test result of yaw rate

從圖9可以看出，r的辨識模型預測輸出與實際飛行測量輸出能夠很好地吻合。經驗證,其他相關狀態(tài)量辨識曲線與試驗測量曲線均能較好地吻合,說明表3中的辨識結果是可信的。

逐步回歸法與最小二乘法均適用于縱向和橫側向參數辨識,得到的結果較為接近,故本文不再贅述。

5 結束語

本文提出了根據試驗可測數據估算迎角和側滑角以及氣動系數進而用于氣動參數辨識的方法。該方法是在無人機試驗條件有限、無法測得準確氣動角的情況下完善飛行數據的一種簡化方法。以某小型飛翼無人機飛行試驗數據為例,驗證了該方法用于無人機氣動參數辨識的有效性,辨識結果滿足一定精度要求。

考慮到辨識建模在數學本質上具有通用性,該方法不僅適用于本文研究的小型飛翼無人機,對其他布局形式的固定翼無人機的氣動參數辨識也是普遍適用的。經檢驗,所得氣動參數辨識結果具有正確性和工程實用性,可以進一步推廣應用于無人機飛行仿真、飛行品質評價和控制系統設計等。

參考文獻：

[1] Younsaeng Lee,Seungjoo Kim,Jinyong Suk,et al.System identification of an unmanned aerial vehicle from automated fight tests[R].AIAA-2002-3493,2002.

[2] 蔡金獅.飛行器系統辨識學[M].北京:國防工業(yè)出版社,2003:247-248.

[3] Brad Seanor,Giampiero Campa,Yu Gu,et al.Formation flight test results for UAV research aircraft models[R].AIAA-2004-6251,2004.

[4] Morelli E A.Identification of low oder equivalent system models from flight test data[R].NASA/TM-2000-210117,2000.

[5] 中國人民解放軍總裝備部軍事訓練班教材編輯工作委員會.模型自由飛試驗[M].北京:國防工業(yè)出版社,2002:181-182.

[6] Etkin B.Dynamics of atmospheric flight[M].Newyork: John Wiley & Sons,Inc.,1972:108-110.

[7] Jategaonkar R V. Flight vehicle system identification-a time domain methodology[M].USA:AIAA,2006:93-94.

[8] 李言俊,張科.系統辨識理論及應用[M].北京:國防工業(yè)出版社,2006:51-63.

[9] Vladislav Klein,Morelli E A.Aircraft system identification-theory and practice[M].USA:AIAA Education Series,2006:47.

[10] Morelli E A.System identification programs for Aircraft[R].AIAA-2002-4704,2002.

[11] Pinkelman J K,Batill S M,Kenoe M W.An investigation of the total least-squares criteria in time domain based,parameter identification for flight flutter testing[R].AIAA-1995-1247-CP,1995.