高超聲速飛行器曲線擬合建模方法

崔建峰, 張科, 呂梅柏

(1.西北工業(yè)大學(xué) 航天學(xué)院, 陜西 西安 710072;2.航天飛行動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室, 陜西 西安 710072)

0 引言

機(jī)體/發(fā)動(dòng)機(jī)一體化的設(shè)計(jì)理念及輕質(zhì)材料的應(yīng)用,使得吸氣式高超聲速飛行器的動(dòng)力學(xué)特性十分復(fù)雜,存在著較為顯著的氣動(dòng)-熱-彈性-推進(jìn)耦合現(xiàn)象。構(gòu)建一個(gè)合理綜合上述交叉耦合因素的動(dòng)力學(xué)模型成為研究吸氣式高超聲速飛行器控制系統(tǒng)的基礎(chǔ)。目前,吸氣式高超聲速飛行器控制系統(tǒng)研究廣泛應(yīng)用的動(dòng)力學(xué)模型之一,是由Bolender等人提出的基于斜基波-膨脹波理論建立的二維高超聲速飛行器動(dòng)力學(xué)模型[1-2]。國(guó)內(nèi)一些學(xué)者也借鑒此種解析建模方法構(gòu)建了相應(yīng)的吸氣式高超聲速飛行器二維動(dòng)力學(xué)模型[3-5]。

上述模型所提供的氣動(dòng)數(shù)據(jù)離散,不便于進(jìn)行非線性控制系統(tǒng)綜合設(shè)計(jì)。因此,文獻(xiàn)[6-8]在解析模型的基礎(chǔ)上對(duì)氣動(dòng)數(shù)據(jù)進(jìn)行曲線擬合,得到面向控制的曲線擬合模型(Curve-Fitted Model,CFM)。氣動(dòng)數(shù)據(jù)曲線擬合是一個(gè)較為耗費(fèi)時(shí)間與計(jì)算機(jī)資源的過(guò)程。例如,為得到對(duì)應(yīng)于7個(gè)自變量、每個(gè)自變量10個(gè)水平值的氣動(dòng)數(shù)據(jù),需要消耗將近80 GB的計(jì)算機(jī)存儲(chǔ)容量,如果每次計(jì)算需要5 ms,則獲取上述數(shù)據(jù)需要計(jì)算將近14 h[7]。所以,擬合前需要依據(jù)可供使用的計(jì)算機(jī)資源選取適當(dāng)簡(jiǎn)化形式的模型,這將限制擬合模型的可信度及其與原始物理模型之間的一致程度。因此,建立一種快速、可靠的曲線擬合建模方法是十分必要的。

本文將在現(xiàn)有彈性高超聲速飛行器動(dòng)力學(xué)模型的基礎(chǔ)上,借助均勻設(shè)計(jì)、逐步回歸等統(tǒng)計(jì)學(xué)方法快速高效地獲取飛行器的CFM模型,最后仿真驗(yàn)證模型的可靠性,并與現(xiàn)有CMF模型進(jìn)行對(duì)比分析。

1 曲線擬合建模

1.1 原始物理模型

本文以文獻(xiàn)[1-2]中給出的高超聲速飛行器動(dòng)力學(xué)模型為曲線擬合參考對(duì)象。該飛行器的幾何結(jié)構(gòu)如圖1所示,外形參數(shù)及其他詳細(xì)信息可參見(jiàn)文獻(xiàn)[2,6-8]。

圖1 高超聲速飛行器模型幾何結(jié)構(gòu)Fig.1 Geometry of the hypersonic vehicle model

該模型主要采用理論分析方法構(gòu)建,飛行器表面壓強(qiáng)分布及氣流參數(shù)利用空氣動(dòng)力學(xué)中的斜基波-膨脹波理論,分別在高超聲速飛行器前體斜坡、發(fā)動(dòng)機(jī)燃燒室、后體膨脹面以及控制舵面處進(jìn)行計(jì)算,飛行器超燃沖壓發(fā)動(dòng)機(jī)的推力則通過(guò)瑞利流模型進(jìn)行估算。該飛行器模型還考慮到彈性效應(yīng)。模型的結(jié)構(gòu)彈性模態(tài)基于單根歐拉-伯努利自由梁,通過(guò)假設(shè)模態(tài)法求解,并利用一階活塞理論計(jì)算由剛體模態(tài)及彈性模態(tài)引起的非定常氣動(dòng)力。另外,模型的飛行力學(xué)方程主要在穩(wěn)定軸坐標(biāo)系下,利用拉格朗日方程推導(dǎo)得出,具體如式(1)所示。

(1)

1.2 曲線擬合模型數(shù)據(jù)生成

傳統(tǒng)的曲線擬合模型求解方法需要獲取大量的擬合數(shù)據(jù)集和驗(yàn)證數(shù)據(jù)集,這個(gè)過(guò)程較為耗費(fèi)時(shí)間與計(jì)算機(jī)資源，且因數(shù)據(jù)集龐大,給后期模型擬合帶來(lái)不便。因此,本文將使用均勻試驗(yàn)設(shè)計(jì)來(lái)簡(jiǎn)化數(shù)據(jù)生成的過(guò)程。

均勻試驗(yàn)設(shè)計(jì)是應(yīng)對(duì)探索性試驗(yàn)獲取近似模型要求而產(chǎn)生的近似試驗(yàn)設(shè)計(jì)方法,也是計(jì)算機(jī)試驗(yàn)設(shè)計(jì)的主要方法之一[9]。借助均勻設(shè)計(jì)方法,可以在需要進(jìn)行多水平試驗(yàn)的情況下,極大地減少待試驗(yàn)次數(shù)。利用均勻設(shè)計(jì)方法,首先要確定試驗(yàn)因素及試驗(yàn)設(shè)計(jì)空間。

表1 狀態(tài)變量工作范圍Table 1 Parameter ranges for curve fitting

在確定試驗(yàn)因素及試驗(yàn)設(shè)計(jì)空間后,一般情況下,當(dāng)試驗(yàn)因素s≤7時(shí),可以直接選用文獻(xiàn)[9]及其他相關(guān)文獻(xiàn)中提供的均勻設(shè)計(jì)表開(kāi)始進(jìn)行仿真試驗(yàn)。由于本文試驗(yàn)因素有12個(gè),因此需要根據(jù)均勻設(shè)計(jì)的原理和方法構(gòu)造指定試驗(yàn)因素的均勻設(shè)計(jì)表。

常用的均勻設(shè)計(jì)表構(gòu)造方法有好格子點(diǎn)法、方冪好格子點(diǎn)法及在其基礎(chǔ)上改進(jìn)的修正法與切割法[9]。本文中的試驗(yàn)因素均為連續(xù)變量,對(duì)變量水平?jīng)]有嚴(yán)格限制,故本文使用方冪好格子點(diǎn)法及其修正方法構(gòu)建均勻設(shè)計(jì)表。詳細(xì)的方法與步驟如下:

(1)給定試驗(yàn)次數(shù)n,尋找比n小的整數(shù)h,且使n和h的最大公約數(shù)為1。符合這些條件的正整數(shù)組成一個(gè)數(shù)組h=(h1,h2,…,hm)。其中m由歐拉函數(shù)決定。

(2)設(shè)計(jì)表的第j列利用下式生成:

uij=ihj[modn]

(2)

這里[modn]表示同余運(yùn)算,若ihj超過(guò)n,則用它減去n的一個(gè)適當(dāng)倍數(shù),使差落在[1,n]之中。uij可以遞推來(lái)生成:

(3)

用上述方法生成的設(shè)計(jì)表記為U(nm),如果劃去Un+1表的最后一行,則獲得修正的設(shè)計(jì)表

U*(nm);

(3)均勻設(shè)計(jì)表應(yīng)用需要使用表,方冪好格子點(diǎn)法利用整數(shù)的同余冪生成使用表。

令a為小于n的整數(shù),且a,a2,…,at在[modn]的意義下互不相同,若at+1[modn]=1,則稱(chēng)a對(duì)n的次數(shù)為t。若a對(duì)n的次數(shù)大于或等于s-1,則可用式(4)作為生成向量。

(a0,a,…,as-1)[modn]

(4)

在一切可能的a(最多n-1個(gè))中去比較相應(yīng)生成表的均勻性,均勻性最優(yōu)的一組生成向量即為使用表。于是,給定n和s,便可以獲得均勻性最佳的設(shè)計(jì)表,此時(shí)記為U(n,ns),如是修正的均勻設(shè)計(jì),則記為U*(n,ns)。

均勻設(shè)計(jì)表的構(gòu)造過(guò)程中,需要在不同生成表間進(jìn)行尋優(yōu)。比較兩個(gè)均勻設(shè)計(jì)表的好壞等價(jià)于比較它們所對(duì)應(yīng)兩組點(diǎn)集的均勻性。常用的均勻性測(cè)度有Lp-星偏差、中心化L2-偏差(CD2)、可卷L2-偏差(WD2)等[9]。根據(jù)上述均勻性測(cè)度的準(zhǔn)確性及計(jì)算復(fù)雜程度,本文選用中心化L2-偏差作為均勻設(shè)計(jì)表的均勻性測(cè)度,其計(jì)算公式為:

(5)

本文以CD2(P)≤0.12為評(píng)價(jià)標(biāo)準(zhǔn),最終獲取的均勻設(shè)計(jì)表為U*(240,24012),其中星號(hào)表示由修正方法生成,該表的中心化L2-偏差CD2(P)=0.1172。

1.3 曲線擬合模型求解

利用上節(jié)獲取的均勻設(shè)計(jì)表U*(240,24012)進(jìn)行仿真試驗(yàn),生成相應(yīng)的氣動(dòng)力、推力及廣義力數(shù)據(jù)。在這些數(shù)據(jù)的基礎(chǔ)上,為獲取各變量間的相互關(guān)系,需要進(jìn)行回歸分析。根據(jù)理論及實(shí)際計(jì)算結(jié)果分析可知,各力與試驗(yàn)因素間并非簡(jiǎn)單的線性關(guān)系,需要考慮到各因素之間的交互作用及因素高次項(xiàng)的影響。但將如此多的因素全部引入到擬合模型中,勢(shì)必過(guò)于復(fù)雜,且可能形成“過(guò)擬合”。因此,需要從這些眾多相關(guān)因素中挑選出對(duì)結(jié)果有顯著影響的部分因素,剔除一些不顯著因子,故可以使用逐步回歸分析法。

逐步回歸的基本步驟是:從一個(gè)因素開(kāi)始,視因素對(duì)結(jié)果作用的顯著程度,從大到小逐步引入回歸方程;在引入每個(gè)新因素的同時(shí),要對(duì)已引入因素進(jìn)行顯著性判斷,當(dāng)原引入因素由于后面因素的引入變得不再顯著時(shí),將其剔除;該過(guò)程反復(fù)進(jìn)行,直至既無(wú)不顯著因素從回歸方程中剔除,又無(wú)顯著因素選入回歸方程時(shí)為止。對(duì)因素的顯著性判斷是通過(guò)針對(duì)偏回歸平方和的假設(shè)檢驗(yàn)實(shí)現(xiàn)的[9]。

通過(guò)觀察單變量對(duì)結(jié)果的影響及參考文獻(xiàn)[6-8]的擬合形式,最終獲取的CFM模型為:

(6a)

(6b)

式中,角度單位為rad;空氣密度則根據(jù)美國(guó)1976年標(biāo)準(zhǔn)大氣表計(jì)算得到。

阻力系數(shù)及推力系數(shù)因受剛體參數(shù)影響更為顯著,因此,在逐步回歸過(guò)程中,有關(guān)的彈性氣動(dòng)導(dǎo)數(shù)項(xiàng)作為非顯著因子被剔除。這與文獻(xiàn)[6-8]中的擬合結(jié)果不一致。文獻(xiàn)[6-8]因采用最小二乘法進(jìn)行擬合,受資源限制提前選取了特定形式的模型,從而在擬合結(jié)果中將非顯著因子引入,處理簡(jiǎn)單劃一,夸大了非顯著因子的作用,使擬合得到的CFM模型可靠性下降。而逐步回歸法以偏回歸平方和的假設(shè)檢驗(yàn)為依據(jù),動(dòng)態(tài)引入顯著因子,使擬合結(jié)果更為可靠,這些可通過(guò)后續(xù)的模型驗(yàn)證與對(duì)比分析體現(xiàn)出來(lái)。

2 CFM模型驗(yàn)證與對(duì)比分析

2.1 CFM模型驗(yàn)證

在獲取CFM模型后,需要對(duì)模型進(jìn)行驗(yàn)證。一般可以通過(guò)復(fù)相關(guān)系數(shù)R與殘余標(biāo)準(zhǔn)偏差S來(lái)檢驗(yàn)回歸方程的顯著性與精度。這兩個(gè)統(tǒng)計(jì)量的定義為:

(7)

對(duì)于驗(yàn)證數(shù)據(jù)集的生成,同樣可以使用均勻設(shè)計(jì)表來(lái)獲取。借助前述均勻設(shè)計(jì)表構(gòu)造方法,獲取的驗(yàn)證用均勻設(shè)計(jì)表為U*(1034,103412),該表的中心化L2-偏差CD2(P)=0.055。

在獲取驗(yàn)證數(shù)據(jù)集后,本文CFM模型各統(tǒng)計(jì)量數(shù)值如表2所示。

表2 CFM模型統(tǒng)計(jì)量Table 2 Statistics of curve-fitted model

從表2中可以看出,飛行器各力回歸方程的復(fù)相關(guān)系數(shù)均接近于1,說(shuō)明CFM模型與原物理模型的一致程度較高。

2.2 CFM模型對(duì)比分析

對(duì)于已有的CFM模型,文獻(xiàn)[8]在吸氣式高超聲速飛行器幾何外形與配置上與本文一致。本文以升力、阻力、推力及俯仰力矩為對(duì)象對(duì)比擬合效果。相比于原始物理模型,計(jì)算所得的統(tǒng)計(jì)量數(shù)值如表3所示,各自的殘差如圖2所示(橫坐標(biāo)P為驗(yàn)證數(shù)據(jù)點(diǎn)的序號(hào))。為了更直觀地查看擬合效果,在圖3中還給出了使用真實(shí)物理模型、本文CFM模型、文獻(xiàn)CFM模型計(jì)算得到氣動(dòng)力和力矩隨速度和迎角的變化圖。

從表3、圖2和圖3中可以看出,本文借用均勻設(shè)計(jì)和逐步回歸法所獲取的CFM模型與真實(shí)物理模型更為一致,模型精度也更好,但用于擬合及驗(yàn)證的數(shù)據(jù)點(diǎn)卻不到2000個(gè),建模過(guò)程耗費(fèi)的時(shí)間與計(jì)算機(jī)資源大為減少。

表3 兩種CFM模型統(tǒng)計(jì)量對(duì)比Table 3 Comparison of statistics between two types of CFM

圖2 CFM模型殘差圖Fig.2 Residual chart of curve-fitted model

圖3 氣動(dòng)力和力矩?cái)M合對(duì)比Fig.3 Comparison of curve-fitting results

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

本文提出了一種借助均勻設(shè)計(jì)、逐步回歸等統(tǒng)計(jì)學(xué)理論建立、評(píng)價(jià)及檢驗(yàn)彈性高超聲速飛行器CFM模型的方法。通過(guò)與現(xiàn)有CFM模型的各項(xiàng)數(shù)據(jù)對(duì)比分析表明,基于均勻設(shè)計(jì)和逐步回歸方法進(jìn)行高超聲速飛行器曲線擬合建模是一種高效、可靠的方法,其所建立的CFM模型與原始物理模型的一致程度優(yōu)于現(xiàn)有CFM模型,而所消耗的時(shí)間與計(jì)算機(jī)資源遠(yuǎn)小于現(xiàn)有擬合方法。

參考文獻(xiàn)：

[1] Bolender M A,Doman D B.A non-linear model for the longitudinal dynamics of a hypersonic air-breathing vehicle[R].AIAA-2005-6255,2005.

[2] Bolender M A,Doman D B.Flight path angle dynamics of air-breathing hypersonic vehicles[R].AIAA-2006-6692,2006.

[3] 張棟,唐碩,李世珍.高超聲速飛行器一體化縱向氣動(dòng)特性分析[J].飛行力學(xué),2012,30(4):328-331.

[4] 曲鑫,李菁菁,宋勛,等.考慮推進(jìn)和氣動(dòng)彈性影響的高超飛行器的建模與控制[J].宇航學(xué)報(bào),2011,32(2):303-309.

[5] 張冉,李惠峰,肖進(jìn).高超聲速飛行器剛體/彈性體耦合動(dòng)力學(xué)建模[J].北京航空航天大學(xué)學(xué)報(bào),2012,38(2):160-165.

[6] Parker J T,Serrani A,Yurkovich S,et al.Control-oriented modeling of an air-breathing hypersonic vehicle[J].Journal of Guidance,Control,and Dynamics,2007,30(3):856-869.

[7] Sigthorsson D O.Control-oriented modeling and output feedback control of hypersonic air-breathing vehicles [D].Ohio:The Ohio State University,2008.

[8] Fiorentini L.Nonlinear adaptive controller design for air-breathing hypersonic vehicles [D].Ohio:The Ohio State University,2010.

[9] 方開(kāi)泰,劉民千,周永道.試驗(yàn)設(shè)計(jì)與建模[M].北京:高等教育出版社,2011:148-201.