李魁， 鄧小龍， 楊希祥， 侯中喜， 周新

(國防科技大學(xué)航天科學(xué)與工程學(xué)院, 長沙 410073)

臨近空間是介于傳統(tǒng)航空與航天區(qū)域之間、尚未為人類所開發(fā)利用大氣層，由于其獨(dú)特的資源優(yōu)勢成為人們關(guān)注的焦點(diǎn)[1]。臨近空間底部的平流層具有氣流穩(wěn)定，空氣流動相對緩慢的特點(diǎn)，特別是在一定時間內(nèi)，存在風(fēng)速較小的準(zhǔn)零風(fēng)層(Quasi-Zero Wind Layer，QZWL)[2]。準(zhǔn)零風(fēng)層一般是指平流層下層20 km高度附近的大氣層，上下層緯向風(fēng)風(fēng)向相反，同時經(jīng)向風(fēng)分量亦很小，是由于平流層風(fēng)場中上下層緯向風(fēng)逆轉(zhuǎn)形成的，是存在于平流層風(fēng)場的一種特殊現(xiàn)象。該區(qū)域可用于部署長時間駐空低動態(tài)飛行器進(jìn)行高分辨率對地觀測、通信中繼等任務(wù)。平流層風(fēng)場對低動態(tài)臨近空間飛行器總體設(shè)計與飛行控制有重要影響[3]。

風(fēng)場環(huán)境的建模與預(yù)測，需要對大量歷史數(shù)據(jù)進(jìn)行處理和分析，有效提高建模與預(yù)測效率和精度值得深入研究[4]。本征正交分解(Proper Orthogonal Decomposition，POD)是風(fēng)場建模中常用的一種方法，將風(fēng)場數(shù)據(jù)通過特征值分解為本征模態(tài)和主坐標(biāo)的組合，其中本征模態(tài)主要與空間因素相關(guān)，主坐標(biāo)主要與時間因素相關(guān)，物理意義明確，通過模態(tài)截斷還能提高計算的效率[5]。

目前，POD方法已被廣泛地應(yīng)用于解決實(shí)際問題數(shù)值模擬的降維問題中。Holmes等[6]采用氣動平均的方法對雙坡屋蓋的表面壓力進(jìn)行測量，指出POD方法是描述鈍體風(fēng)壓的有效途徑。Fic等[7]將POD方法的離散形式應(yīng)用于線性和非線性瞬態(tài)熱傳導(dǎo)的降階模型中。Tan等[8]應(yīng)用POD結(jié)合三次樣條插值方法對流場壓力系數(shù)進(jìn)行預(yù)測。對于風(fēng)場建模，胡亮[9]利用POD方法對具有橋塔風(fēng)效應(yīng)的風(fēng)場進(jìn)行了簡化模擬，將POD型譜表示法取代原型譜表示法，物理意義明確且可通過模態(tài)截斷節(jié)省計算量。陶青秋[10]將POD方法用于分析在大跨屋面和高層建筑物上的風(fēng)致響應(yīng)，利用POD與結(jié)構(gòu)模態(tài)相結(jié)合的雙模態(tài)方法來分析結(jié)構(gòu)的動態(tài)響應(yīng)，表明POD方法是一種有效分析多激勵響應(yīng)的方法。鄒垚等[11]通過對POD原理進(jìn)行修正后用于風(fēng)洞實(shí)驗風(fēng)場的重建，主要通過修正第一個本征模態(tài)值，對偏差較大測點(diǎn)的風(fēng)壓場重建進(jìn)行改進(jìn)，使重建風(fēng)場的方差值與實(shí)際風(fēng)場方差值的差距大大縮小。Chen和Kareem[12]利用POD方法對建筑結(jié)構(gòu)的數(shù)值風(fēng)場進(jìn)行了模擬并提出了狀態(tài)空間建模方案，建立了風(fēng)載荷與建筑結(jié)構(gòu)共振響應(yīng)在物理意義上的聯(lián)系，有助于更好地了解結(jié)構(gòu)隨時空變化的動態(tài)負(fù)載響應(yīng)。

本文基于適用于平流層風(fēng)場環(huán)境的POD方法，以長沙地區(qū)風(fēng)場數(shù)據(jù)為例，研究平流層風(fēng)場隨高度和時間的變化特點(diǎn)，通過POD方法對風(fēng)場數(shù)據(jù)進(jìn)行降階處理，實(shí)現(xiàn)在特定時間范圍內(nèi)快速準(zhǔn)確計算出風(fēng)場輪廓線。鑒于所選地區(qū)平流層20 km附近東西方向風(fēng)速呈現(xiàn)明顯的空間相關(guān)性和周期性，在POD方法的基礎(chǔ)上，提出了一種可以對風(fēng)場進(jìn)行預(yù)測的Fourier模型。

1 基于POD方法的風(fēng)場建模

POD方法是指從一組龐大的數(shù)據(jù)中獲得一組低維最優(yōu)基，而這組低維最優(yōu)基是一種能夠表示原復(fù)雜系統(tǒng)的降階模型。也就是說POD方法是一種降維的方法，對給定的數(shù)據(jù)進(jìn)行最優(yōu)的低維逼近，用較小的維數(shù)將原物理模型的主要特征表現(xiàn)出來[13]。本文通過POD方法獲得風(fēng)場信息的最優(yōu)標(biāo)準(zhǔn)正交基，將風(fēng)場數(shù)據(jù)中的每個風(fēng)場信息投影到標(biāo)準(zhǔn)正交基上，獲得每個風(fēng)場數(shù)據(jù)在標(biāo)準(zhǔn)正交基上的投影系數(shù)[14]。并通過最優(yōu)標(biāo)準(zhǔn)正交基與相對應(yīng)的投影系數(shù)計算出一組新的數(shù)據(jù)，將原風(fēng)場數(shù)據(jù)的主要特征表現(xiàn)出來。

1.1 瞬像矩陣的生成

瞬像也稱為快照，由Sirovich提出[15]，是指物理場的數(shù)值解在不同時刻的空間分布值。建立海拔高度h1至hL的風(fēng)場模型,取瞬像矩陣為

(1)

式中：V為風(fēng)速；下標(biāo)L為瞬像時刻數(shù)或快照數(shù)。假設(shè)風(fēng)場模型統(tǒng)計天數(shù)為M，則A∈RL×M，一般由試驗、經(jīng)驗、統(tǒng)計等方法獲得，通常M?L。

1.2 POD模態(tài)的獲取

POD模態(tài)又稱為POD基。POD方法的目的是尋找一組最優(yōu)的POD基φ1,φ2,…,φr(φi為列向量)，使得數(shù)據(jù)模型中任意一天特定時刻的風(fēng)場可以表示為POD基的線性組合[16]，即

V(t)=α1(t)φ1+α2(t)φ2+…+αr(t)φr

(2)

式中：αi(t)為POD基的系數(shù)，i=1,2,…,r。為了尋求滿足式(2)的最優(yōu)POD基，定義相關(guān)矩陣：

R=ATA

(3)

則R∈RL×L，求取相關(guān)矩陣R的非零特征值以及特征向量φi：

Rφi=λiφi

i=1,2,…,n;λ1≥λ2≥…≥λn>0;n≤L

(4)

將特征向量化為最優(yōu)POD基：

(5)

且φi(i=1,2,…,n；n≤L)是一組標(biāo)準(zhǔn)正交基，即

(6)

只選取前r階模態(tài)，POD模態(tài)矩陣[17]記為

(7)

投影系數(shù)矩陣記為

(8)

則降階模型矩陣B=ΦΩ∈RL×M。

1.3 降階模型誤差評估

POD模態(tài)表示捕獲物理場的主要特征，截斷后的前r階POD模態(tài)較全階模態(tài)所捕獲的能量比為

(9)

故定義相對截斷誤差[18]為

(10)

2 基于Fourier級數(shù)的風(fēng)場預(yù)測模型

通過第1節(jié)的POD方法夠獲得風(fēng)場數(shù)據(jù)在最優(yōu)POD基上的投影系數(shù)，若具有明顯的空間相關(guān)性和周期性，可以使用Fourier級數(shù)直接進(jìn)行表示。通過最優(yōu)POD基與Fourier級數(shù)所表示的投影系數(shù)計算出特定高度的風(fēng)速隨時間變化情況[19]。其方程表示如下：

(11)

式中：t為時間，d；m為Fourier級數(shù)展開的項數(shù)；aki和bki為Fourier系數(shù)；T為投影系數(shù)的振蕩周期。

3 實(shí)例分析與結(jié)果

以長沙地區(qū)為例，采用2005年1月1日—2009年12月31日每天00：00所測量的風(fēng)場數(shù)據(jù)，風(fēng)場數(shù)據(jù)由國防科技大學(xué)氣象海洋學(xué)院提供，數(shù)據(jù)記錄了長沙地區(qū)00：00和12：00(世界時)的風(fēng)速、濕度、溫度、氣壓等氣象要素，數(shù)據(jù)覆蓋海拔高度0～60 km之間每隔0.5 km所對應(yīng)的天氣數(shù)據(jù)?？紤]到臨近空間低動態(tài)飛行器駐空期間的工作高度范圍，本文針對海拔高度10～30 km的范圍進(jìn)行研究，并通過3次樣條對風(fēng)場數(shù)據(jù)進(jìn)行了插值處理。

3.1 風(fēng)場建模與預(yù)測

相對模態(tài)能量分布主要用于識別數(shù)據(jù)降階處理所需的模態(tài)數(shù)量。圖1(a)給出了東西方向風(fēng)場各階POD模態(tài)的相對模態(tài)能量，前5種模態(tài)能量所占的比重相對于其他模態(tài)能量具有明顯的區(qū)別，可以捕獲到風(fēng)場總特征性的98.9%。同時，由圖1(b)可知，南北方向風(fēng)場的相對模態(tài)能量分布除前2種模態(tài)之外，其余模態(tài)能量呈逐漸下降的趨勢，若只選取前2種模態(tài)進(jìn)行降階處理，不足以表現(xiàn)出風(fēng)場的主要特征。因此，南北方向風(fēng)場不能高效地進(jìn)行降階建模。

圖1 相對模態(tài)能量分布Fig.1 Relative mode energy distribution

通過選用多種POD模型分別對2007年11月6日0：00時刻東西方向和南北方向的風(fēng)場進(jìn)行建模。由圖2(a)可知，10階和15階POD模型較好地表示出東西方向風(fēng)場的實(shí)際數(shù)據(jù)，5階POD模型能體現(xiàn)出東西方向風(fēng)場的主要特征。圖2(b)中實(shí)際風(fēng)場模型清楚地表明南北方向風(fēng)速隨空間的變化是極其不規(guī)則的，導(dǎo)致南北方向風(fēng)場的低階重建性較差。因此，低階的POD模型(比如5階)在捕獲東西方向風(fēng)場數(shù)據(jù)的總體趨勢上是有效的。

經(jīng)過POD方法降階處理后可以得到相對應(yīng)的投影系數(shù)。以第一階POD模態(tài)的投影系數(shù)Ck1為例，通過Fourier級數(shù)擬合，可以得出其Fourier預(yù)測模型，分別從東西方向(見圖3)和南北方向(見圖4)進(jìn)行分析。從圖3(a)和圖4(a)可看出，兩者的投影系數(shù)Ck1都具有周期性和規(guī)律性，都能夠通過Fourier級數(shù)將其擬合(分別見圖3(b)和圖4(b))。但相對于南北方向風(fēng)場，東西方向風(fēng)場的投影系數(shù)更緊密，周期性更強(qiáng)，F(xiàn)ourier級數(shù)的擬合程度更高，說明東西方向風(fēng)場更適用于建立Fourier預(yù)測模型。

圖2 采用POD方法重建風(fēng)場Fig.2 Reconstruction of wind field using POD method

利用東西方向風(fēng)場的Fourier預(yù)測模型，可以獲得某一天特定時刻下的投影系數(shù)，能對該時刻下平流層20 km附近的東西方向風(fēng)速隨高度的變化情況進(jìn)行預(yù)測。如圖5所示，以2008年11月7日為例，5階POD模型能夠捕獲到實(shí)際風(fēng)場的主要特征，F(xiàn)ourier預(yù)測模型能夠體現(xiàn)出5階POD模型的整體變化趨勢狀況。

圖3 系數(shù)擬合(東西方向)Fig.3 Coefficient fitting (east-west direction)

圖4 系數(shù)擬合(南北方向)Fig.4 Coefficient fitting (north-south direction)

圖5 各種方法下風(fēng)矢量圖的比較Fig.5 Comparison of wind vector plots for various methods

將Fourier預(yù)測模型向外延伸和擴(kuò)展，可以對未來特定時刻下的風(fēng)場進(jìn)行預(yù)測。通過2005—2009年風(fēng)場數(shù)據(jù)的統(tǒng)計和分析，對2010年10月30日0:00平流層20 km附近的東西方向風(fēng)速進(jìn)行預(yù)測，結(jié)果如圖6所示，F(xiàn)ourier預(yù)測模型能預(yù)測到風(fēng)速隨高度的變化情況，并且能夠體現(xiàn)出平流層低速風(fēng)帶所在的位置高度。

圖6 實(shí)際風(fēng)矢量圖與預(yù)測風(fēng)矢量圖的對比Fig.6 Comparison of actual wind vector with predicted wind vector plots

圖8 不同高度的風(fēng)速變化情況Fig.8 Change of wind speed at different altitudes

圖9 不同高度的殘差分析Fig.9 Residual analysis at different altitudes

3.2 Fourier預(yù)測精度分析

由圖5、圖6可知，針對不同的高度，F(xiàn)ourier預(yù)測模型與實(shí)際風(fēng)場模型之間存在不同的偏差，說明Fourier預(yù)測精度與空間分布有關(guān)，下面通過模型之間的比較，對Fourier預(yù)測精度進(jìn)行分析。

將圖5中3種風(fēng)場模型分別進(jìn)行比較，如圖7(a)所示，其中POD降階殘差表示5階POD風(fēng)場模型與實(shí)際風(fēng)場之差，F(xiàn)ourier擬合殘差表示Fourier級數(shù)進(jìn)行擬合時與5階POD模型之差，F(xiàn)ourier預(yù)測殘差表示Fourier預(yù)測模型與實(shí)際風(fēng)場之差。由圖7(a)可知，3種殘差的絕對值都小于5，其中POD降階殘差曲線的差波動性最小，說明5階POD模型與實(shí)際風(fēng)場模型數(shù)據(jù)很接近。由于Fourier預(yù)測模型是基于POD方法降階處理后建立的，所以Fourier預(yù)測殘差等于POD降階殘差與Fourier擬合殘差兩者之和，說明Fourier預(yù)測精度不僅與POD降階模型所捕獲到的能量有關(guān)，而且和Fourier級數(shù)的擬合程度有關(guān)。將圖6中Fourier預(yù)測模型與實(shí)際風(fēng)場模型進(jìn)行比較，如圖7(b)所示，F(xiàn)ourier預(yù)測殘差曲線的波動范圍較小，在海拔高度19～22 km之間的殘差接近于0，說明在此范圍內(nèi)，F(xiàn)ourier預(yù)測精度較高，而在15 km附近波動最大，說明預(yù)測精度較低。

由圖7可知，在20 km處Fourier預(yù)測精度相對較高，在15和24 km處的預(yù)測精度較低。分別對15、20和24 km，3個高度的東西方向風(fēng)速隨時間的變化情況進(jìn)行分析，并分別計算出實(shí)際風(fēng)場模型、5階POD模型與Fourier預(yù)測模型兩兩之間的殘差，結(jié)果如圖8、圖9所示。由圖9可知，在這3個高度位置上，POD降階殘差都接近于零，說明對于5階POD模型，在任意高度上都能夠捕獲到實(shí)際風(fēng)場的主要特征；對于Fourier預(yù)測模型，如圖8所示，在20 km處東西方向的風(fēng)場數(shù)據(jù)相對于其他2個高度的數(shù)據(jù)更加緊湊、平穩(wěn)，F(xiàn)ourier級數(shù)的擬合程度也更高，說明Fourier預(yù)測模型與實(shí)際風(fēng)場隨時間變化的緊湊性與規(guī)律性有著密切聯(lián)系，風(fēng)場數(shù)據(jù)越緊湊，周期性越明顯，F(xiàn)ourier預(yù)測精度就越高。

圖7 殘差分析Fig.7 Residual analysis

4 結(jié) 論

本文通過POD方法將平流層風(fēng)場數(shù)據(jù)進(jìn)行高精度高效率降階處理后，為Fourier預(yù)測模型對風(fēng)場進(jìn)行預(yù)測節(jié)省了大量的計算時間，大大提高了計算效率，得到：

1) 提出了一種對平流層風(fēng)場數(shù)據(jù)進(jìn)行降階處理的POD方法，在POD方法的基礎(chǔ)上，提出了一種可以對平流層風(fēng)場進(jìn)行預(yù)測的Fourier預(yù)測模型。

2) 以長沙地區(qū)為例，選取海拔高度10～30 km的5年風(fēng)場數(shù)據(jù)，采用提出的POD方法與Fourier預(yù)測模型對風(fēng)場進(jìn)行了建模與預(yù)測，并對Fourier預(yù)測精度進(jìn)行了分析。

3) 研究結(jié)果表明：采用POD方法可以對東西方向風(fēng)場進(jìn)行高效高精度的降階建模，由于南北方向風(fēng)場變化極其不規(guī)則，不可以采用低階的POD模型進(jìn)行建模；通過Fourier預(yù)測模型能對東西方向風(fēng)場進(jìn)行準(zhǔn)確預(yù)測，預(yù)測精度與實(shí)際風(fēng)場隨時間變化的規(guī)律性有關(guān)，風(fēng)場數(shù)據(jù)越緊湊，周期性越明顯，F(xiàn)ourier預(yù)測精度越高。