基于自回歸滑動平均模型的大口徑天線風速預測方法*

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

(1 新疆大學物理科學與技術學院烏 魯木齊 830046)(2 中國科學院新疆天文臺 烏魯木齊 830011)(3 中國科學院射電天文重點實驗室 烏魯木齊 830011)(4 新疆射電天體物理重點實驗室 烏魯木齊 830011)(5 西安電子科技大學電子裝備結構設計教育部重點實驗室 西安 710071)(6 中國科學院大學 北京 100049)

1 引言

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

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

為降低風擾對天線指向的影響,國內外學者往往采用天線罩、增強天線結構剛度和強度等手段降低天線結構變形,或者通過優化伺服控制器以增強控制系統魯棒性[9]以及利用固定補償等方法降低風擾對天線指向的影響.雖然使用天線罩能降低風擾對天線的影響,但制作天線罩的材料會吸收或反射電磁波而降低天線電性能.且天線罩制作、運輸安裝難度隨口徑增大而急劇增加,僅適用于小口徑天線.伺服控制補償都是針對靜態載荷導致的天線變形或振動.對于時變的風載影響,利用風預測提前得到作用在天線上的風速,為控制系統的計算和執行提供足夠時間,才能更好地降低天線在風擾下執行觀測任務時風載荷導致的指向抖動.

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

本研究以QTT臺址的風速預測為目標,通過采集QTT臺址風場數據分析臺址風場特征,并將不同季度的風場數據劃分為訓練集和測試集.通過對不同季度訓練集數據進行建模,即在訓練集數據平穩性檢驗、模型定階、參數估計、殘差檢驗的基礎上,建立ARMA模型進行風預測,獲取臺址風速預測值,與測試值比較計算預測精度.

2 臺址風特性分析

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

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

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

圖2 臺址2017年4–12月風速概率分布圖.(a)時間序列中提取的分布;(b)時間序列中提取的累積分布.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風速風向玫瑰圖中,每個扇形所指方向表示從外部吹向測風塔的風向,扇形上不同的灰度代表不同風速.例如正南方向的扇形數據顯示,正南方向0–2 m·s-1風速所占時間約為6%,2–4 m·s-1風速所占時間約為2%,4–6 m·s-1風速所占時間約為1.5%,6–8 m·s-1風速所占時間約為1%,8–10 m·s-1、10–12 m·s-1風速所占時間均小于1%.從圖3風速風向玫瑰圖還可以得出QTT臺址風在西南、南、東南和西北4個風向出現較多,且高風速在西南風向出現頻次更高.

圖3 臺址2017年4–12月風速風向玫瑰圖Fig.3 Wind rose diagram from April to December 2017 of QTT site

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

圖4 不同季度的風速風向玫瑰圖.(a)夏季風玫瑰圖;(b)秋季風玫瑰圖;(c)冬季風玫瑰圖.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模型的風速預測方法

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

3.1 時間序列模型

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

(1)AR模型

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

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

(2)MA模型

白噪聲的線性組合構成了MA模型假設序列,通常來說,階數為q的MA模型,記作MA(q),其數學表達式為:

式中,θ1,θ2,···,θq是滑動平均系數.

(3)ARMA模型

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

其中,p、q分別為自回歸階數和滑動平均階數.當q=0時,ARMA(p,q)模型轉變為AR(p)模型;當p=0時,ARMA(p,q)模型轉變為MA(q)模型.

3.2 預測模型的構建

ARMA模型要求風速數據{xt}為平穩、正態、零均值的時間序列,因此在建模之前需要對數據進行預處理,從而得到符合要求的風速時間序列.由圖2可以看出臺址的風速不具有正態分布特性,但ARMA模型預測主要是利用風速數據的相關性,要求風速時間序列的基本特性不變,對其正態分布特性要求不嚴格,本文不再對風速數據進行正態化處理,這里只需要對時間序列進行差分處理和零均值標準化處理.但要注意不能過度差分,否則會導致時間序列的方差增大.

(1)風速時間序列的預處理

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

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

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

圖6 訓練集風速時間序列相關函數.(a)夏季自相關函數;(b)夏季偏自相關函數;(c)秋季自相關函數;(d)秋季偏自相關函數;(e)冬季自相關函數;(f)冬季偏自相關函數.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 差分后訓練集風速時間序列相關函數.(a)夏季自相關函數;(b)夏季偏自相關函數;(c)秋季自相關函數;(d)秋季偏自相關函數;(e)冬季自相關函數;(f)冬季偏自相關函數.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)模型識別

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

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

圖8 不同季度風速預測模型的最優階數.(a)夏季ARMA模型的最優階數是(5,5);(b)秋季ARMA模型的最優階數是(5,4);(c)冬季ARMA模型的最優階數是(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)模型參數估計

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

summer:

autumn:

winter:

(4)模型診斷檢驗

經過階次辨識和參數估計得到的模型,需要進行殘差分析.如果分析結果較差,需要重新進入“辨識-估計-診斷”的流程.對于殘差分析,殘差正態分布的檢驗通過判斷直方圖或分位數-分位數(Quantile-Quantile,Q-Q)圖中模型的殘差分布與標準正態分布的擬合程度來實現;殘差隨機性檢驗通過判斷樣本自相關函數和偏自相關函數中相關系數是否滿足截尾來實現.

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

圖9 夏季ARMA模型殘差檢驗.(a)標準化殘差;(b)殘差分布;(c)Q-Q圖;(d)自相關函數;(e)偏自相關函數.從圖(b)中可以看出模型滿足殘差檢驗的正態分布;從圖(c)中可以看出模型的殘差分布擬合標準正態化分布效果較好;圖(d)自相關函數和圖(e)偏自相關函數中相關系數在95%置信區間下滿足趨近于零的條件,說明殘差序列無明顯相關性且隨機性強,模型滿足殘差檢驗的隨機性.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模型殘差檢驗.(a)標準化殘差;(b)殘差分布;(c)Q-Q圖;(d)自相關函數;(e)偏自相關函數.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模型殘差檢驗.(a)標準化殘差;(b)殘差分布;(c)Q-Q圖;(d)自相關函數;(e)偏自相關函數.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 預測結果及分析

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

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

圖13 不同季度的部分預測值與測試值的對比.(a)夏季部分風速預測值;(b)秋季部分風速預測值;(c)冬季部分風速預測值.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.

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

式中yi和?yi分別表示第i時刻風速觀測值和預測值.

本文先預測30 min風速數據,將其加入訓練集再預測下一個30 min的風速數據,經過多次滾動預測,從而實現108 h風速數據預測.由于誤差隨時間累積,我們發現4 h后預測數據誤差較大,對望遠鏡控制執行時間來說,4 h完全能夠滿足,因此選擇3個季度預測前4 h風速預測值與測試值的偏差衡量模型精度,結果如表1所示.

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

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

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

表2 本文風速預測精度與其他文獻中風速預測精度比較Table 2 The accuracy comparison of wind prediction between this article and other work

在時間尺度上,4 h范圍的風預測數據能夠為射電望遠鏡控制系統運動機構提供足夠的執行時間,并且能指導射電望遠鏡在此段時間觀測任務的安排.故本研究基于不同季度風速數據建立的ARMA模型能夠很好地預測臺址風速,為射電望遠鏡的抗風擾控制提供必要數據支撐.

5 結論

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

在風速變化較大的拐點處,本研究建立的風速預測模型預測精度不夠高,隨著預測時間的增加,預測精度也逐漸變低,因此未來還需要對算法進行優化.另外,由于本研究采用統計預測方法,僅利用風速歷史數據進行建模預測,預測效果有限.因此與數值天氣預報有限結合,進一步提升預測精度是未來研究的一個重點.

臺址風場的風速預測不是最終目的,如何利用預測結果和大口徑射電望遠鏡主副伺服控制系統建立聯系,共同保障風擾下射電望遠鏡的觀測性能,是下一步的研究工作.