李 琳 許 謙 王文娟 李 帥 薛 松 連培園 王從思 何飛龍

(1 新疆大學(xué)物理科學(xué)與技術(shù)學(xué)院烏 魯木齊 830046)(2 中國科學(xué)院新疆天文臺(tái) 烏魯木齊 830011)(3 中國科學(xué)院射電天文重點(diǎn)實(shí)驗(yàn)室 烏魯木齊 830011)(4 新疆射電天體物理重點(diǎn)實(shí)驗(yàn)室 烏魯木齊 830011)(5 西安電子科技大學(xué)電子裝備結(jié)構(gòu)設(shè)計(jì)教育部重點(diǎn)實(shí)驗(yàn)室 西安 710071)(6 中國科學(xué)院大學(xué) 北京 100049)

1 引言

大型反射面天線具有高增益、窄波束的特點(diǎn),被廣泛應(yīng)用于射電天文、深空探測等領(lǐng)域[1].隨著引力波探測、恒星形成、星系起源等科學(xué)研究的深入,需要進(jìn)一步提高天線的增益和分辨率,為此需要增大口徑或提升工作頻段[2].但是天線口徑增大會(huì)導(dǎo)致天線波束變窄,從而對(duì)天線指向精度提出更苛刻的要求.例如我國即將建造的新疆奇臺(tái)110 m口徑全向可動(dòng)射電望遠(yuǎn)鏡(QiTai Telescope,QTT),為滿足科學(xué)目標(biāo)需求,當(dāng)其工作頻率為115 GHz時(shí),指向精度要達(dá)到1.5′′(0.000416°)[3].對(duì)于大口徑天線而言,其工作往往處于露天環(huán)境中,面臨重力、風(fēng)擾、溫差、積雪等環(huán)境載荷導(dǎo)致指向精度下降的問題,尤其是具有隨機(jī)性和時(shí)變性的風(fēng)擾,是環(huán)境載荷中最復(fù)雜的問題之一[4–5].

國內(nèi)外許多學(xué)者已經(jīng)發(fā)現(xiàn)風(fēng)擾對(duì)天線指向精度的影響.西安電子科技大學(xué)Duan等人利用數(shù)值分析對(duì)7.3 m天線變形引起的指向誤差進(jìn)行研究,得到天線仰角為0°時(shí),方位軸受到靜載引起指向誤差為0.0001°,當(dāng)天線受到20 m·s-1的靜態(tài)風(fēng)載時(shí),其指向誤差增大為0.0135°[6].美國國家航天局通過對(duì)天線伺服系統(tǒng)建模,將風(fēng)擾動(dòng)作為干擾力矩進(jìn)行分析,發(fā)現(xiàn)當(dāng)平均風(fēng)速為3.47 m·s-1時(shí),編碼器檢測到34 m天線的指向誤差為0.00048°,70 m天線的指向誤差為0.0013°[7–8].日本國立天文臺(tái)研究風(fēng)致10 m天線指向誤差,得到在10 m·s-1的風(fēng)速下,10 m天線的指向誤差達(dá)到0.00047°[9].

為降低風(fēng)擾對(duì)天線指向的影響,國內(nèi)外學(xué)者往往采用天線罩、增強(qiáng)天線結(jié)構(gòu)剛度和強(qiáng)度等手段降低天線結(jié)構(gòu)變形,或者通過優(yōu)化伺服控制器以增強(qiáng)控制系統(tǒng)魯棒性[9]以及利用固定補(bǔ)償?shù)确椒ń档惋L(fēng)擾對(duì)天線指向的影響.雖然使用天線罩能降低風(fēng)擾對(duì)天線的影響,但制作天線罩的材料會(huì)吸收或反射電磁波而降低天線電性能.且天線罩制作、運(yùn)輸安裝難度隨口徑增大而急劇增加,僅適用于小口徑天線.伺服控制補(bǔ)償都是針對(duì)靜態(tài)載荷導(dǎo)致的天線變形或振動(dòng).對(duì)于時(shí)變的風(fēng)載影響,利用風(fēng)預(yù)測提前得到作用在天線上的風(fēng)速,為控制系統(tǒng)的計(jì)算和執(zhí)行提供足夠時(shí)間,才能更好地降低天線在風(fēng)擾下執(zhí)行觀測任務(wù)時(shí)風(fēng)載荷導(dǎo)致的指向抖動(dòng).

在風(fēng)預(yù)測領(lǐng)域,Brown等[10]利用自回歸模型(Autoregressive,AR)考慮了太平洋西北地區(qū)某風(fēng)場數(shù)據(jù)的非高斯平穩(wěn)性和晝夜非平穩(wěn)性,采用冪變換和標(biāo)準(zhǔn)化進(jìn)行了建模預(yù)測.采用AR(2)模型,Poggi等[11]對(duì)法國科西地區(qū)的風(fēng)速數(shù)據(jù)進(jìn)行了仿真與預(yù)測,結(jié)果表明,AR(2)模型能夠很好地復(fù)現(xiàn)原始序列,反映風(fēng)場數(shù)據(jù)的統(tǒng)計(jì)特性.丁明等[12]采用自回歸滑動(dòng)平均模型(Autoregressive Moving Average Model,ARMA)對(duì)風(fēng)速進(jìn)行預(yù)測,比較了預(yù)測風(fēng)速和實(shí)際風(fēng)速的分布特性.研究重點(diǎn)關(guān)注如何使用時(shí)間序列模型模擬風(fēng)速序列,使其具有風(fēng)的一般特征,盡管這些模型都涉及到風(fēng)預(yù)測,但在實(shí)際預(yù)測中,并沒有使用單獨(dú)預(yù)測集來驗(yàn)證模型的預(yù)測性能.即以上使用滑動(dòng)自回歸模型研究的重點(diǎn)是模型在訓(xùn)練數(shù)據(jù)上的逼近能力,而非模型預(yù)測能力.

本研究以QTT臺(tái)址的風(fēng)速預(yù)測為目標(biāo),通過采集QTT臺(tái)址風(fēng)場數(shù)據(jù)分析臺(tái)址風(fēng)場特征,并將不同季度的風(fēng)場數(shù)據(jù)劃分為訓(xùn)練集和測試集.通過對(duì)不同季度訓(xùn)練集數(shù)據(jù)進(jìn)行建模,即在訓(xùn)練集數(shù)據(jù)平穩(wěn)性檢驗(yàn)、模型定階、參數(shù)估計(jì)、殘差檢驗(yàn)的基礎(chǔ)上,建立ARMA模型進(jìn)行風(fēng)預(yù)測,獲取臺(tái)址風(fēng)速預(yù)測值,與測試值比較計(jì)算預(yù)測精度.

2 臺(tái)址風(fēng)特性分析

QTT臺(tái)址位于新疆昌吉自治州奇臺(tái)縣半截溝鎮(zhèn)石河子村,臺(tái)址海拔約1730–1830 m,位于東天山北麓一處南北約2 km、東西約1.5 km且四面環(huán)山的矩形盆地[3].臺(tái)址特征是盆地東南高西北低,西側(cè)是落差在110 m左右的中葛根河谷,四周分布的山體海拔高度約為1860–2250 m.臺(tái)址內(nèi)已建一個(gè)60 m高的梯度風(fēng)塔,臺(tái)址地形和測風(fēng)塔的位置如圖1所示,圖中線框所圍區(qū)域?yàn)镼TT臺(tái)址園區(qū),F為測風(fēng)塔位置,T為QTT位置[13].

圖1 QTT射電望遠(yuǎn)鏡與測風(fēng)塔位置.(a)測風(fēng)塔;(b)QTT天線.Fig.1 QTT site and position of anemometer tower.(a)anemometer tower;(b)QTT antenna.

風(fēng)塔數(shù)據(jù)采集設(shè)備包括2D超聲波風(fēng)速風(fēng)向傳感器和CR3000微型采集器,風(fēng)塔采集頻率為每分鐘1次,采集時(shí)長為1 s[13],這樣按照一定時(shí)間間隔采集和記錄的風(fēng)數(shù)據(jù)構(gòu)成臺(tái)址風(fēng)場的時(shí)間序列,該時(shí)間序列包含產(chǎn)生該序列系統(tǒng)的歷史行為全部信息.研究所用數(shù)據(jù)為在QTT臺(tái)址利用測風(fēng)塔采集的2017年4–12月10 m高度的風(fēng)數(shù)據(jù).通過對(duì)4–12月風(fēng)速時(shí)間序列中不同風(fēng)速出現(xiàn)的頻率分析,得出如圖2所示臺(tái)址風(fēng)速概率分布圖以及圖3所示風(fēng)速風(fēng)向玫瑰圖.在圖2風(fēng)速概率分布圖中可以看出在2017年4–12月臺(tái)址風(fēng)速集中在1–8 m·s-1的范圍.

圖2 臺(tái)址2017年4–12月風(fēng)速概率分布圖.(a)時(shí)間序列中提取的分布;(b)時(shí)間序列中提取的累積分布.Fig.2 Probability distribution of wind speed from April to December 2017 of QTT site.(a)Distribution extracted from the time series;(b)Cumulative distribution extracted from the time series.

在圖3風(fēng)速風(fēng)向玫瑰圖中,每個(gè)扇形所指方向表示從外部吹向測風(fēng)塔的風(fēng)向,扇形上不同的灰度代表不同風(fēng)速.例如正南方向的扇形數(shù)據(jù)顯示,正南方向0–2 m·s-1風(fēng)速所占時(shí)間約為6%,2–4 m·s-1風(fēng)速所占時(shí)間約為2%,4–6 m·s-1風(fēng)速所占時(shí)間約為1.5%,6–8 m·s-1風(fēng)速所占時(shí)間約為1%,8–10 m·s-1、10–12 m·s-1風(fēng)速所占時(shí)間均小于1%.從圖3風(fēng)速風(fēng)向玫瑰圖還可以得出QTT臺(tái)址風(fēng)在西南、南、東南和西北4個(gè)風(fēng)向出現(xiàn)較多,且高風(fēng)速在西南風(fēng)向出現(xiàn)頻次更高.

圖3 臺(tái)址2017年4–12月風(fēng)速風(fēng)向玫瑰圖Fig.3 Wind rose diagram from April to December 2017 of QTT site

為觀察臺(tái)址風(fēng)數(shù)據(jù)在不同季度的差異,本研究按照傳統(tǒng)四季分布即4、5、6月為夏季,7、8、9月為秋季,10、11、12月為冬季,得到3個(gè)季度的風(fēng)玫瑰圖,如圖4所示.從3個(gè)季度的風(fēng)玫瑰圖可得出風(fēng)向多集中于南、北兩個(gè)風(fēng)向,相較于南向的風(fēng),在北向風(fēng)中4–6 m·s-1風(fēng)速所占時(shí)間較多,但大于6 m·s-1的風(fēng)速多集中于南向風(fēng).夏秋兩個(gè)季度風(fēng)速風(fēng)向特征變化并不明顯,冬季臺(tái)址的風(fēng)向多集中在東南偏南方向,大于6 m·s-1的風(fēng)速基本不出現(xiàn).3個(gè)季度的風(fēng)玫瑰圖對(duì)比來看,QTT臺(tái)址的風(fēng)場數(shù)據(jù)具有明顯季節(jié)差異,應(yīng)按照不同季度數(shù)據(jù)建立不同風(fēng)預(yù)測模型.

圖4 不同季度的風(fēng)速風(fēng)向玫瑰圖.(a)夏季風(fēng)玫瑰圖;(b)秋季風(fēng)玫瑰圖;(c)冬季風(fēng)玫瑰圖.Fig.4 Wind rose diagram of different quarters.(a)Wind rose diagram in summer;(b)Wind rose diagram in autumn;(c)Wind rose diagram in winter.

3 基于ARMA模型的風(fēng)速預(yù)測方法

根據(jù)預(yù)測的時(shí)間尺度不同,風(fēng)預(yù)測可分為超短期預(yù)測、短期預(yù)測、中期預(yù)測和長期預(yù)測[14–15].把提前幾分鐘至30 min左右的預(yù)測定義為超短期預(yù)測,而短期預(yù)測則是指提前30 min至72 h左右的預(yù)測,而把提前幾天或數(shù)周、數(shù)月的預(yù)測稱之為中期預(yù)測,長期預(yù)測是以年為預(yù)測單位.對(duì)于望遠(yuǎn)鏡來說,短期預(yù)測可以給望遠(yuǎn)鏡觀測任務(wù)的規(guī)劃提供參考,超短期預(yù)測可以給望遠(yuǎn)鏡伺服控制提供輸入數(shù)據(jù),為控制系統(tǒng)的運(yùn)動(dòng)機(jī)構(gòu)爭取執(zhí)行時(shí)間,在風(fēng)作用到望遠(yuǎn)鏡時(shí)完成姿態(tài)調(diào)整.以輸入數(shù)據(jù)分類風(fēng)預(yù)測主要包括物理方法和統(tǒng)計(jì)方法.物理方法的輸入數(shù)據(jù)通常為各種物理描述量,包括地形地貌特征、氣象信息、地表粗糙度、障礙物等地理信息.其優(yōu)點(diǎn)是不會(huì)過度依賴風(fēng)歷史數(shù)據(jù),且考慮了氣象、地理等復(fù)雜物理量,可獲得較精確的中長期預(yù)測值,缺點(diǎn)是需要建立精確的地理模型,且計(jì)算也十分復(fù)雜,需借助超算才可能實(shí)現(xiàn).統(tǒng)計(jì)方法則是將風(fēng)場歷史數(shù)據(jù)作為輸入數(shù)據(jù),通過提取歷史數(shù)據(jù)中的輸入輸出映射關(guān)系構(gòu)建統(tǒng)計(jì)模型,從而對(duì)將來的風(fēng)速進(jìn)行預(yù)測.其優(yōu)點(diǎn)是不考慮風(fēng)速產(chǎn)生的復(fù)雜物理過程,計(jì)算簡單[16].故本研究采用統(tǒng)計(jì)方法預(yù)測超短期風(fēng)速.

3.1 時(shí)間序列模型

經(jīng)典時(shí)間序列模型包括:自回歸模型、滑動(dòng)平均(MA)模型以及自回歸滑動(dòng)平均模型3種.時(shí)間序列模型的建模步驟主要包括模型階次辨識(shí)、參數(shù)估計(jì)與診斷.

(1)AR模型

AR模型的基本思想是當(dāng)前風(fēng)速值xt可以通過p個(gè)自身的過去值xt-1,xt-2,···,xt-p來解釋.其中p表示預(yù)測當(dāng)前值所需要的歷史步長,該模型假定序列當(dāng)前值是過去值的線性函數(shù).階數(shù)為p的AR模型,一般記作AR(p),數(shù)學(xué)表達(dá)式為:

式中φ1,φ2,···,φp是自回歸系數(shù),εt是白噪聲序列.

(2)MA模型

白噪聲的線性組合構(gòu)成了MA模型假設(shè)序列,通常來說,階數(shù)為q的MA模型,記作MA(q),其數(shù)學(xué)表達(dá)式為:

式中,θ1,θ2,···,θq是滑動(dòng)平均系數(shù).

(3)ARMA模型

ARMA模型假設(shè)序列的當(dāng)前值由其過去值和系統(tǒng)噪聲的過去值線性組合而成,即時(shí)間序列的當(dāng)前值,其不僅取決于過去某一特定時(shí)間內(nèi)的歷史數(shù)據(jù),還取決于當(dāng)前和過去時(shí)刻引入系統(tǒng)中的噪聲.ARMA(p,q)的數(shù)學(xué)表達(dá)式為:

其中,p、q分別為自回歸階數(shù)和滑動(dòng)平均階數(shù).當(dāng)q=0時(shí),ARMA(p,q)模型轉(zhuǎn)變?yōu)锳R(p)模型;當(dāng)p=0時(shí),ARMA(p,q)模型轉(zhuǎn)變?yōu)镸A(q)模型.

3.2 預(yù)測模型的構(gòu)建

ARMA模型要求風(fēng)速數(shù)據(jù){xt}為平穩(wěn)、正態(tài)、零均值的時(shí)間序列,因此在建模之前需要對(duì)數(shù)據(jù)進(jìn)行預(yù)處理,從而得到符合要求的風(fēng)速時(shí)間序列.由圖2可以看出臺(tái)址的風(fēng)速不具有正態(tài)分布特性,但ARMA模型預(yù)測主要是利用風(fēng)速數(shù)據(jù)的相關(guān)性,要求風(fēng)速時(shí)間序列的基本特性不變,對(duì)其正態(tài)分布特性要求不嚴(yán)格,本文不再對(duì)風(fēng)速數(shù)據(jù)進(jìn)行正態(tài)化處理,這里只需要對(duì)時(shí)間序列進(jìn)行差分處理和零均值標(biāo)準(zhǔn)化處理.但要注意不能過度差分,否則會(huì)導(dǎo)致時(shí)間序列的方差增大.

(1)風(fēng)速時(shí)間序列的預(yù)處理

為驗(yàn)證預(yù)測效果,本研究將每個(gè)季度風(fēng)速時(shí)間序列分為95%的訓(xùn)練集和5%的測試集.即利用一個(gè)季度2052 h的風(fēng)速數(shù)據(jù)預(yù)測108 h的風(fēng)速數(shù)據(jù),然后將預(yù)測數(shù)據(jù)與測試集數(shù)據(jù)進(jìn)行對(duì)比驗(yàn)證.3個(gè)季度訓(xùn)練集風(fēng)速數(shù)據(jù)如圖5所示,訓(xùn)練集風(fēng)速時(shí)間序列的自相關(guān)函數(shù)和偏自相關(guān)函數(shù)如圖6所示,其中自相關(guān)函數(shù)(Autocorrelation Function,ACF)是表示當(dāng)前時(shí)刻的數(shù)值與過去時(shí)刻的數(shù)值之間的相關(guān)程度,包括直接和間接的相關(guān)性信息;偏自相關(guān)函數(shù)(Partial Autocorrelation Function,PACF)是表示當(dāng)前時(shí)刻的數(shù)值與過去時(shí)刻的數(shù)值的直接相關(guān)程度.相關(guān)函數(shù)中自相關(guān)系數(shù)隨滯后值的增加迅速衰減為零的時(shí)間序列是平穩(wěn)序列,反之是非平穩(wěn)序列.由于ARMA模型要求時(shí)間序列數(shù)據(jù)必須平穩(wěn),因此本研究利用ADF(Augmented Dickey-Fuller Test)檢驗(yàn)和KPSS(Kwiatkowski-Phillips-Schmidt-Shin Test)檢驗(yàn)等單位根檢驗(yàn)方法確定訓(xùn)練集風(fēng)速時(shí)間序列的平穩(wěn)性,自相關(guān)函數(shù)中相關(guān)系數(shù)是否衰減為零也可以作為判斷時(shí)間序列是否平穩(wěn)的依據(jù),可以與單位根檢驗(yàn)平穩(wěn)性的方法相互驗(yàn)證,增加平穩(wěn)性檢驗(yàn)的準(zhǔn)確性.

圖5 不同季度的訓(xùn)練集風(fēng)速時(shí)間序列圖.(a)夏季;(b)秋季;(c)冬季.Fig.5 Wind speed time series of training set in different quarters.(a)Summer;(b)Autumn;(c)Winter.

圖6中自相關(guān)系數(shù)沒有隨著滯后值迅速衰減至0,因此可以判定該時(shí)序數(shù)據(jù)屬于非平穩(wěn)序列,需要對(duì)訓(xùn)練集數(shù)據(jù)進(jìn)行一階差分,差分公式為Δxt=xt+1-xt,差分后每個(gè)季度風(fēng)數(shù)據(jù)時(shí)間序列的自相關(guān)函數(shù)和偏自相關(guān)函數(shù)如圖7所示.可以得到每個(gè)季度自相關(guān)系數(shù)隨著滯后值增加很快衰減至0,判斷3個(gè)季度風(fēng)速時(shí)間序列屬于平穩(wěn)序列,故而不再進(jìn)行差分.

圖6 訓(xùn)練集風(fēng)速時(shí)間序列相關(guān)函數(shù).(a)夏季自相關(guān)函數(shù);(b)夏季偏自相關(guān)函數(shù);(c)秋季自相關(guān)函數(shù);(d)秋季偏自相關(guān)函數(shù);(e)冬季自相關(guān)函數(shù);(f)冬季偏自相關(guān)函數(shù).Fig.6 Correlation function of wind speed time series of training set.(a)Autocorrelation function in summer;(b)Partial autocorrelation function in summer;(c)Autocorrelation function in autumn;(d)Partial autocorrelation function in autumn;(e)Autocorrelation function in winter;(f)Partial autocorrelation function in winter.

圖7 差分后訓(xùn)練集風(fēng)速時(shí)間序列相關(guān)函數(shù).(a)夏季自相關(guān)函數(shù);(b)夏季偏自相關(guān)函數(shù);(c)秋季自相關(guān)函數(shù);(d)秋季偏自相關(guān)函數(shù);(e)冬季自相關(guān)函數(shù);(f)冬季偏自相關(guān)函數(shù).Fig.7 Correlation function of wind speed time series of training set after difference.(a)Autocorrelation function in summer;(b)Partial autocorrelation function in summer;(c)Autocorrelation function in autumn;(d)Partial autocorrelation function in autumn;(e)Autocorrelation function in winter;(f)Partial autocorrelation function in winter.

(2)模型識(shí)別

模型識(shí)別包括類別判斷與階數(shù)選取,通常通過風(fēng)速數(shù)據(jù)中相關(guān)函數(shù)的自相關(guān)性進(jìn)行模型類別判斷,但實(shí)際時(shí)間序列數(shù)據(jù)具有很強(qiáng)的隨機(jī)性,序列中會(huì)隨機(jī)出現(xiàn)異常值波動(dòng)從而導(dǎo)致相關(guān)系數(shù)的截尾性不明顯,所以常利用準(zhǔn)則函數(shù)進(jìn)行模型階數(shù)選取.常用的準(zhǔn)則函數(shù)是基于信息論的赤池信息準(zhǔn)則(Akaike Information Criterion,AIC)[17]和貝葉斯信息準(zhǔn)則(Bayesian Information Criterion,BIC)[18],通過如下公式計(jì)算每個(gè)模型的AIC值或BIC值:

式中,L是模型的似然函數(shù),N是序列的長度,w是噪聲誤差的參數(shù),?σ2w是噪聲誤差的最大似然估計(jì),r=p+q+1表示模型中待估計(jì)參數(shù)的數(shù)目.AIC和BIC準(zhǔn)則的表達(dá)式均由兩項(xiàng)組成,第1項(xiàng)表示模型的擬合程度,一般模型的階次越大其值越小;第2項(xiàng)表示對(duì)過多參數(shù)的懲罰,模型階次越大其值越大.對(duì)各種可能的模型,具有最小值的AIC或BIC即為最終選定的模型,每個(gè)季度的不同階數(shù)預(yù)測模型的AIC和BIC值如圖8所示,每個(gè)灰度格子的數(shù)值對(duì)應(yīng)不同階數(shù)的AIC和BIC值之和,該數(shù)值結(jié)果越小說明模型越優(yōu)異,從而確定預(yù)測模型的階數(shù).其中夏季最佳預(yù)測模型為ARMA(5,5),秋季的最佳預(yù)測模型為ARMA(5,4),冬季的最佳預(yù)測模型為ARMA(3,5).

圖8 不同季度風(fēng)速預(yù)測模型的最優(yōu)階數(shù).(a)夏季ARMA模型的最優(yōu)階數(shù)是(5,5);(b)秋季ARMA模型的最優(yōu)階數(shù)是(5,4);(c)冬季ARMA模型的最優(yōu)階數(shù)是(3,5).Fig.8 Optimal order of wind speed prediction model in different quarters.(a)The optimal order of summer ARMA model is(5,5);(b)The optimal order of autumn ARMA model is(5,4);(c)The optimal order of winter ARMA model is(3,5).

(3)模型參數(shù)估計(jì)

模型階次確定后,可以根據(jù)風(fēng)速時(shí)間序列x1,x2,···,xt進(jìn)行參數(shù)估計(jì),得到φ1,φ2,···,φp和ε1,ε2,···,εq的估計(jì)值.常用的估計(jì)方法如:矩估計(jì)、最小二乘估計(jì)和最大似然估計(jì).矩估計(jì)是讓樣本矩和相應(yīng)的理論矩相等,通過求解方程組得到未知參數(shù)的估計(jì),對(duì)于ARMA模型而言,矩估計(jì)計(jì)算復(fù)雜且無法得到最優(yōu)估計(jì).最小二乘估計(jì)基于誤差平方和最小原則得到參數(shù)估計(jì)值,當(dāng)時(shí)間序列數(shù)據(jù)數(shù)目較大時(shí),使用一些初始值對(duì)最終的參數(shù)估計(jì)影響不大,數(shù)目較少時(shí),采用無條件最小二乘或最大似然進(jìn)行參數(shù)估計(jì).最大似然估計(jì)通過選取模型參數(shù)使得似然函數(shù)最大化,似然函數(shù)是模型參數(shù)的函數(shù),代表了時(shí)間序列值出現(xiàn)的可能性.本研究利用estimate函數(shù)對(duì)3個(gè)季度定階的模型進(jìn)行最大似然估計(jì),分別得到以下模型:

summer:

autumn:

winter:

(4)模型診斷檢驗(yàn)

經(jīng)過階次辨識(shí)和參數(shù)估計(jì)得到的模型,需要進(jìn)行殘差分析.如果分析結(jié)果較差,需要重新進(jìn)入“辨識(shí)-估計(jì)-診斷”的流程.對(duì)于殘差分析,殘差正態(tài)分布的檢驗(yàn)通過判斷直方圖或分位數(shù)-分位數(shù)(Quantile-Quantile,Q-Q)圖中模型的殘差分布與標(biāo)準(zhǔn)正態(tài)分布的擬合程度來實(shí)現(xiàn);殘差隨機(jī)性檢驗(yàn)通過判斷樣本自相關(guān)函數(shù)和偏自相關(guān)函數(shù)中相關(guān)系數(shù)是否滿足截尾來實(shí)現(xiàn).

圖9展現(xiàn)了夏季ARMA模型的殘差檢驗(yàn),圖9(a)是對(duì)夏季模型殘差進(jìn)行標(biāo)準(zhǔn)化處理,得到了零均值且方差為一的標(biāo)準(zhǔn)化殘差序列;圖9(b)是對(duì)標(biāo)準(zhǔn)化殘差序列求得殘差的數(shù)值分布,分布呈現(xiàn)正態(tài)分布,說明滿足殘差檢驗(yàn)的正態(tài)分布;圖9(c)中可以看出代表模型殘差分布的十字標(biāo)與代表標(biāo)準(zhǔn)正態(tài)分布的虛線擬合效果較好,模型滿足殘差檢驗(yàn)的正態(tài)分布;圖9(d)自相關(guān)函數(shù)和圖9(e)偏自相關(guān)函數(shù)中,點(diǎn)表示對(duì)應(yīng)滯后值的自相關(guān)系數(shù)或偏自相關(guān)系數(shù),線表示自相關(guān)系數(shù)或偏自相關(guān)系數(shù)滿足相關(guān)性截尾的95%置信區(qū)間上下限,可以看出相關(guān)系數(shù)在95%置信區(qū)間下滿足趨近于零的條件,由此說明殘差序列無明顯相關(guān)性且隨機(jī)性強(qiáng),模型滿足殘差檢驗(yàn)的隨機(jī)性.夏季ARMA模型滿足模型殘差檢驗(yàn)的要求,模型的殘差序列是白噪聲信號(hào),原始風(fēng)速序列中有用信息已被提取到預(yù)測模型中.圖10展現(xiàn)了秋季ARMA模型的殘差檢驗(yàn),圖11展現(xiàn)了冬季ARMA模型的殘差檢驗(yàn),從圖10和圖11中可以看出秋季模型和冬季模型殘差也是隨機(jī)正態(tài)分布且不自相關(guān).所以,可以判定3個(gè)季度模型類別和階數(shù)選用合適.

圖9 夏季ARMA模型殘差檢驗(yàn).(a)標(biāo)準(zhǔn)化殘差;(b)殘差分布;(c)Q-Q圖;(d)自相關(guān)函數(shù);(e)偏自相關(guān)函數(shù).從圖(b)中可以看出模型滿足殘差檢驗(yàn)的正態(tài)分布;從圖(c)中可以看出模型的殘差分布擬合標(biāo)準(zhǔn)正態(tài)化分布效果較好;圖(d)自相關(guān)函數(shù)和圖(e)偏自相關(guān)函數(shù)中相關(guān)系數(shù)在95%置信區(qū)間下滿足趨近于零的條件,說明殘差序列無明顯相關(guān)性且隨機(jī)性強(qiáng),模型滿足殘差檢驗(yàn)的隨機(jī)性.Fig.9 Residual test of ARMA model in summer.(a)Standardized residual;(b)Residual distribution;(c)Q-Q plot;(d)Autocorrelation function;(e)Partial autocorrelation function.It can be seen from panel(b)that the model satisfies the normal distribution of residual test;It can be seen from panel(c)that the residual distribution of the model fits the standard normal distribution well.The correlation coefficient in the ACF of panel(d)and PACF of panel(e)satisfies the condition of approaching zero under the 95% confidence interval,indicating that the residual sequence has no obvious correlation and strong randomness,and the model satisfies the randomness of the residual test.

圖10 秋季ARMA模型殘差檢驗(yàn).(a)標(biāo)準(zhǔn)化殘差;(b)殘差分布;(c)Q-Q圖;(d)自相關(guān)函數(shù);(e)偏自相關(guān)函數(shù).Fig.10 Residual test of ARMA model in autumn:(a)Standardized residual;(b)Residual distribution;(c)Q-Q plot;(d)Autocorrelation function;(e)Partial autocorrelation function.

圖11 冬季ARMA模型殘差檢驗(yàn).(a)標(biāo)準(zhǔn)化殘差;(b)殘差分布;(c)Q-Q圖;(d)自相關(guān)函數(shù);(e)偏自相關(guān)函數(shù).Fig.11 Residual test of ARMA model in winter.(a)Standardized residual;(b)Residual distribution;(c)Q-Q plot;(d)Autocorrelation function;(e)Partial autocorrelation function.

4 預(yù)測結(jié)果及分析

利用提出的不同季度ARMA模型進(jìn)行預(yù)測,給出95%置信區(qū)間得到預(yù)測數(shù)據(jù).圖12顯示了3個(gè)季度的預(yù)測值與測試值的對(duì)比圖,圖中灰色線代表訓(xùn)練集數(shù)據(jù),點(diǎn)線代表測試集數(shù)據(jù),黑線代表模型預(yù)測值,可得ARMA模型能夠很好預(yù)測風(fēng)速.圖13選取3個(gè)季度部分預(yù)測值與測試值,點(diǎn)線為測試值,黑線為預(yù)測值,能清晰看出預(yù)測值與測試值的吻合程度,證明本研究建立的不同季度ARMA模型能夠很好擬合真實(shí)值,預(yù)測精度較高.但模型對(duì)拐點(diǎn)處的值預(yù)測還不夠準(zhǔn)確,還存在一定誤差,預(yù)測也存在一定滯后現(xiàn)象.

圖12 不同季度的預(yù)測值與測試值的對(duì)比.(a)夏季;(b)秋季;(c)冬季.Fig.12 Comparison of prediction data and test data in different quarters.(a)Summer;(b)Autumn;(c)Winter.

圖13 不同季度的部分預(yù)測值與測試值的對(duì)比.(a)夏季部分風(fēng)速預(yù)測值;(b)秋季部分風(fēng)速預(yù)測值;(c)冬季部分風(fēng)速預(yù)測值.Fig.13 Comparison of prediction data and test data in different quarters:(a)Prediction data of partial wind speed in summer;(b)Prediction data of partial wind speed in autumn;(c)Prediction data of partial wind speed in winter.

為定量比較預(yù)測精度,使用均方根誤差(Root Mean Square Error,RMSE)、平均絕對(duì)誤差(Mean Absolute Error,MAE)與平均絕對(duì)百分誤差(Mean Absolute Percentage Error,MAPE)衡量預(yù)測值與測試值之間偏差,這3個(gè)指標(biāo)取值越小代表預(yù)測精度越高,他們被廣泛用于風(fēng)速預(yù)測領(lǐng)域[19–20].其定義如下:

式中yi和?yi分別表示第i時(shí)刻風(fēng)速觀測值和預(yù)測值.

本文先預(yù)測30 min風(fēng)速數(shù)據(jù),將其加入訓(xùn)練集再預(yù)測下一個(gè)30 min的風(fēng)速數(shù)據(jù),經(jīng)過多次滾動(dòng)預(yù)測,從而實(shí)現(xiàn)108 h風(fēng)速數(shù)據(jù)預(yù)測.由于誤差隨時(shí)間累積,我們發(fā)現(xiàn)4 h后預(yù)測數(shù)據(jù)誤差較大,對(duì)望遠(yuǎn)鏡控制執(zhí)行時(shí)間來說,4 h完全能夠滿足,因此選擇3個(gè)季度預(yù)測前4 h風(fēng)速預(yù)測值與測試值的偏差衡量模型精度,結(jié)果如表1所示.

表1 不同季度ARMA模型預(yù)測效果Table 1 Prediction effect of ARMA model in different quarters

從不同季度ARMA模型預(yù)測效果來看,夏季平均絕對(duì)誤差最小,平均絕對(duì)百分誤差最大,這是由于夏季的測試集前4 h的平均風(fēng)速較冬季和秋季小.整體來看,3個(gè)季度預(yù)測效果較為接近,秋季風(fēng)速預(yù)測精度最高.

為比較本文風(fēng)速預(yù)測精度的優(yōu)劣,將本文3個(gè)季度模型誤差值取平均,并選取3篇文獻(xiàn)的風(fēng)速預(yù)測精度作為參考.其中Jiang等[21]對(duì)山東蓬萊某風(fēng)電場的風(fēng)速進(jìn)行了預(yù)測,Song等[22]對(duì)中國某風(fēng)電場的風(fēng)速短期預(yù)測進(jìn)行了研究,Zhang等[23]對(duì)西班牙某風(fēng)電場的風(fēng)速進(jìn)行了預(yù)測,以上文獻(xiàn)中風(fēng)速預(yù)測精度與本文風(fēng)速預(yù)測精度如表2所示.可以看出本文風(fēng)速預(yù)測的平均絕對(duì)誤差和均方根誤差均優(yōu)于以上文獻(xiàn).本文的風(fēng)速預(yù)測模型具有一定的優(yōu)越性.

表2 本文風(fēng)速預(yù)測精度與其他文獻(xiàn)中風(fēng)速預(yù)測精度比較Table 2 The accuracy comparison of wind prediction between this article and other work

在時(shí)間尺度上,4 h范圍的風(fēng)預(yù)測數(shù)據(jù)能夠?yàn)樯潆娡h(yuǎn)鏡控制系統(tǒng)運(yùn)動(dòng)機(jī)構(gòu)提供足夠的執(zhí)行時(shí)間,并且能指導(dǎo)射電望遠(yuǎn)鏡在此段時(shí)間觀測任務(wù)的安排.故本研究基于不同季度風(fēng)速數(shù)據(jù)建立的ARMA模型能夠很好地預(yù)測臺(tái)址風(fēng)速,為射電望遠(yuǎn)鏡的抗風(fēng)擾控制提供必要數(shù)據(jù)支撐.

5 結(jié)論

風(fēng)因其間歇性、波動(dòng)性以及隨機(jī)性等特點(diǎn),給大口徑射電望遠(yuǎn)鏡的抗風(fēng)擾控制補(bǔ)償帶來嚴(yán)峻挑戰(zhàn),克服這一難題的有效途徑是對(duì)臺(tái)址風(fēng)場進(jìn)行預(yù)測.本研究根據(jù)臺(tái)址風(fēng)塔采集到2017年4-12月風(fēng)速風(fēng)向數(shù)據(jù),分析了臺(tái)址風(fēng)場特征,發(fā)現(xiàn)不同季度風(fēng)速風(fēng)向有較大差異.在此基礎(chǔ)上通過對(duì)不同季度風(fēng)速數(shù)據(jù)的平穩(wěn)性檢驗(yàn)和數(shù)據(jù)處理,利用AIC準(zhǔn)則確定模型階數(shù),最大似然法估算模型參數(shù),最終建立了不同季度的ARMA模型,并進(jìn)行了風(fēng)速預(yù)測,從均方根誤差、平均絕對(duì)誤差與平均絕對(duì)百分誤差的值來看,風(fēng)速預(yù)測取得很好的效果,能夠滿足大口徑射電望遠(yuǎn)鏡抗風(fēng)擾控制在大多數(shù)情況下的使用.

在風(fēng)速變化較大的拐點(diǎn)處,本研究建立的風(fēng)速預(yù)測模型預(yù)測精度不夠高,隨著預(yù)測時(shí)間的增加,預(yù)測精度也逐漸變低,因此未來還需要對(duì)算法進(jìn)行優(yōu)化.另外,由于本研究采用統(tǒng)計(jì)預(yù)測方法,僅利用風(fēng)速歷史數(shù)據(jù)進(jìn)行建模預(yù)測,預(yù)測效果有限.因此與數(shù)值天氣預(yù)報(bào)有限結(jié)合,進(jìn)一步提升預(yù)測精度是未來研究的一個(gè)重點(diǎn).

臺(tái)址風(fēng)場的風(fēng)速預(yù)測不是最終目的,如何利用預(yù)測結(jié)果和大口徑射電望遠(yuǎn)鏡主副伺服控制系統(tǒng)建立聯(lián)系,共同保障風(fēng)擾下射電望遠(yuǎn)鏡的觀測性能,是下一步的研究工作.