基于時間序列的山區(qū)風(fēng)向預(yù)測方法研究

武 略 焦瑞莉* 夏江江

（1北京信息科技大學(xué)信息與通信工程學(xué)院，北京 100101；2中國科學(xué)院大氣物理研究所，北京 100029；3中國科學(xué)院大學(xué)，北京 100049）

山嶺地區(qū)地形條件復(fù)雜，受背風(fēng)坡和迎風(fēng)坡的影響，氣溫、風(fēng)和能見度等要素在水平范圍變化劇烈。其中，風(fēng)向作為常見的預(yù)報要素，具有重要的實際意義和很高的研究價值，尤其對山區(qū)的風(fēng)力發(fā)電、基礎(chǔ)設(shè)施建設(shè)等具有重要意義。風(fēng)向是2022年北京冬季奧運組委會十分關(guān)注的氣象條件之一[1]，具有較大的波動性和不可預(yù)測性，對運動員的成績和安全有關(guān)鍵影響。

目前，風(fēng)向主要通過物理建模、統(tǒng)計學(xué)和數(shù)據(jù)驅(qū)動3種方法進(jìn)行預(yù)測。物理方法主要為模式的預(yù)報，例如歐洲中期天氣預(yù)報中心針對全球格點所做的預(yù)測數(shù)據(jù)，原理通常為分析大氣運動規(guī)律，構(gòu)建物理學(xué)方程模擬地理環(huán)境和天氣狀況，從而推演風(fēng)場的變化規(guī)律。此類預(yù)測方法需要掌握大氣的物理相互作用，并建立合理規(guī)范的重現(xiàn)方程，但是由于模式的物理化方案尚存在欠缺，許多計算參數(shù)也不確定，使近地面風(fēng)場預(yù)報存在較大的誤差[2]。在統(tǒng)計學(xué)方法方面，Erdem等[3]利用自回歸滑動平均模型（auto regressive moving average，ARMA）對風(fēng)向序列進(jìn)行了分析和預(yù)測，曾曉青等[4]利用傳統(tǒng)的模式輸出統(tǒng)計法進(jìn)行了風(fēng)向矢量預(yù)測，祝 牧等[5]使用隨機Markov鏈擬合了風(fēng)向的時間序列，并建立模型進(jìn)行預(yù)測。近年來，更高效的數(shù)據(jù)驅(qū)動算法如機器學(xué)習(xí)和深度學(xué)習(xí)開始應(yīng)用于風(fēng)向預(yù)測。Mohandes等[6]提出利用支持向量機對風(fēng)向進(jìn)行預(yù)測，張東東等[7]將BP神經(jīng)網(wǎng)絡(luò)應(yīng)用在風(fēng)力發(fā)電機的風(fēng)向預(yù)測中，唐振浩等[8]基于SWLSTM算法對超短期風(fēng)向進(jìn)行了預(yù)測。

為滿足風(fēng)向預(yù)報的需求，算法應(yīng)有較高的準(zhǔn)確率，還應(yīng)易于操作維護(hù)、復(fù)雜度低、建模快。結(jié)合山區(qū)的地形條件和氣候變化，所選用的預(yù)測方法需對山脈地區(qū)的數(shù)據(jù)變化足夠敏感，既能分析突變的風(fēng)向數(shù)據(jù)，又可解決預(yù)測不準(zhǔn)確的問題，還需及時使用最新數(shù)據(jù)再次搭建模型反饋結(jié)果，實現(xiàn)靈活準(zhǔn)確的預(yù)報目標(biāo)。本文分別使用ARMA統(tǒng)計方法、XGBoost機器學(xué)習(xí)算法、LSTM深度學(xué)習(xí)算法對選定山區(qū)4個氣象站的時間觀測序列進(jìn)行預(yù)測，比較3種算法的優(yōu)劣，并在此基礎(chǔ)上添加風(fēng)速序列，將風(fēng)速和風(fēng)向拆分為U、V風(fēng)向，并對其分別進(jìn)行預(yù)測，將結(jié)果合成后能更好地預(yù)測高山站風(fēng)向，進(jìn)而為2022年冬奧會的風(fēng)向預(yù)測探索一種新方法。

1 數(shù)據(jù)與方法

1.1 數(shù)據(jù)來源

為了更好地模擬山區(qū)風(fēng)向的預(yù)測實況，選擇北京延慶燕山山脈的4個氣象觀測站作為試驗站點，站點信息見表1。觀測站記錄的風(fēng)向數(shù)據(jù)以時間序列的形式每小時記錄1次，將此數(shù)據(jù)作為后續(xù)的分析數(shù)據(jù)。

表1 山區(qū)4個氣象觀測站點信息

1.2 算法模型

1.2.1 自回歸滑動平均模型。自回歸滑動平均模型（ARMA）作為研究時間序列的經(jīng)典方法，由自回歸模型（簡稱AR模型）與滑動平均模型（簡稱MA模型）根據(jù)數(shù)據(jù)自身結(jié)構(gòu)特點組合而成[9]。一般用ARMA模型擬合時間序列，預(yù)測該時間序列的未來值[10]。ARMA模型可表示為ARMA（p，q），其中 p是自回歸階數(shù)，q是移動平均階數(shù)，如式（1）所示：

式中：xt-1、xt-2、 ……、xt-p為不同時間點記錄的指標(biāo)數(shù)值；?1、?2、……、?p為自回歸系數(shù)；μt、μt-1、……、μt-q為不同時間點的白噪聲項；θ1、θ2、……、θq為移動回歸方程系數(shù)；xt表示時間點t指標(biāo)數(shù)值。

使用ARMA時，需要先判斷所使用的序列是否平穩(wěn)：若平穩(wěn)，可直接使用ARMA模型；若非平穩(wěn)，需要對原序列進(jìn)行差分直至平穩(wěn)，再使用ARIMA模型，其中I代表差分次數(shù)。根據(jù)序列自相關(guān)和偏相關(guān)系數(shù)來確定模型的選擇，需要定階來確定p和q。其中：p代表時間序列是否存在周期，比如氣溫和降水隨季節(jié)變化，就存在明顯的年周期性，在氣溫預(yù)測中，往年同期數(shù)據(jù)有著很高的參考價值；q則代表白噪聲項。常用的定階方法為AIC[11]或BIC準(zhǔn)則，AIC準(zhǔn)則（akaike information criterion）是擬合精度和參數(shù)個數(shù)的加權(quán)函數(shù)，BIC（bayesian information criterion）則改善了AIC準(zhǔn)則在大樣本情況下模型不收斂的問題。在確定模型之后還需要檢驗?zāi)Ｐ停唧w包括觀察殘差序列的隨機性、是否符合正態(tài)分布、是否為白噪聲。

在Python的使用中，通常從statsmodels包中調(diào)用ARMA，statsmodels是一個包含多種統(tǒng)計模型、適用于數(shù)據(jù)分析的Python模塊；也可使用其中的plot_acf和plot_pacf繪制序列的自相關(guān)與偏相關(guān)圖。預(yù)測結(jié)束后，對預(yù)測出來的數(shù)據(jù)進(jìn)行逆差分操作，即得到最終預(yù)測結(jié)果。

1.2.2 極端梯度提升算法。極端梯度提升算法（extreme gradient boosting，XGBoost）是 Boosting 算法的其中一種，Boosting作為一種提升方法，通過擬合殘差進(jìn)而優(yōu)化目標(biāo)函數(shù)，從而將一組弱分類器集合成強分類器。為防止過擬合，基于提升算法，XGBoost對損失函數(shù)進(jìn)行二階泰勒展開，并且加入正則項，衡量目標(biāo)函數(shù)的下降和模型的復(fù)雜水平[12]。模型公式如式（2）所示：

式中，y?i為輸出的預(yù)測值，K 為樹的數(shù)量，fk表示第k棵樹模型。

對每一棵樹進(jìn)行訓(xùn)練時，目標(biāo)函數(shù)如式（3）所示：

式中：L為損失函數(shù)；N為樣本數(shù)；t表示訓(xùn)練第t棵樹；ft表示第 t輪所生成的樹模型；Ω（fi）表示正則項。

樹的復(fù)雜度如式（4）所示：

XGBoost的主要參數(shù)：eta，即學(xué)習(xí)率參數(shù)，值越小，模型對數(shù)據(jù)的學(xué)習(xí)越精細(xì)；max_depth，樹的最大深度，用來避免過擬合，值越大，模型對局部樣本的學(xué)習(xí)更具體；min_child_weight，同樣避免過擬合；gamma控制節(jié)點分裂的標(biāo)準(zhǔn)，值越大，算法越保守；subsample，采樣率，用于調(diào)節(jié)模型擬合程度；colsample_bytree，用于選擇所生成樹的特征；n_estimators，迭代次數(shù)，即生成樹的個數(shù)。

調(diào)整參數(shù)時，使Python軟件中的GridSearchCV來尋找模型的最佳參數(shù)：先調(diào)節(jié)n_estimators，范圍從400到800，步長為100，確定最優(yōu)值后縮小步長再次尋找。然后調(diào)節(jié)控制樹結(jié)構(gòu)的min_child_weight和max_depth這2個參數(shù)；之后依次調(diào)節(jié)gamma、subsample和colsample_bytree；最后調(diào)節(jié)學(xué)習(xí)率，因為其數(shù)值較小，所以從0.01開始調(diào)節(jié)。由于分別預(yù)測了風(fēng)速和風(fēng)向，因而需要在每次使用新數(shù)據(jù)前開始調(diào)參，從而獲得最好結(jié)果。

1.2.3 長短期記憶網(wǎng)絡(luò)。長短期記憶網(wǎng)絡(luò)（long short memory network，LSTM）作為一種特殊的循環(huán)神經(jīng)網(wǎng)絡(luò)，用于解決對較長序列訓(xùn)練過程中的出現(xiàn)梯度消失與爆炸問題[13]。LSTM通過輸入門、輸出門和遺忘門的門控狀態(tài)來控制傳輸狀態(tài)，有選擇地記住、遺忘或更新歷史信息。其結(jié)構(gòu)如圖1所示。

LSTM的遺忘門用來決定遺棄細(xì)胞狀態(tài)中的哪些信息，數(shù)學(xué)模型為：

式中：ft為遺忘門的輸出；σ為激活函數(shù)，取值范圍為[0，1]；xt為 t時刻輸入；ht-1為 t-1 時刻輸出；Wf和bf為參數(shù)矩陣，下同。

輸入門的作用是選擇信息放入細(xì)胞狀態(tài)中，數(shù)學(xué)模型如式（6）（7）所示：

式中：it為輸入門輸出；tanh同為激活函數(shù)，取值范圍[0，1]。

在每次輸出之前，需要更新信息，丟棄舊狀態(tài)的某些信息，即：

式中：Ct為t時刻細(xì)胞狀態(tài)。

最后，輸出門控制信息輸出，ot代表輸出門的輸出，數(shù)學(xué)模型為：

基于輸入門、輸出門和遺忘門組成的結(jié)構(gòu)特點，LSTM具有了對歷史信息選擇性遺忘或更新的能力，能夠更好地分析時間序列的變化趨勢，從而實現(xiàn)篩選過去信息并且結(jié)合當(dāng)前信息，預(yù)測未來時刻信息的功能[14]。

1.3 預(yù)測難點

在傳統(tǒng)的時間序列預(yù)測中，多采用單時間序列，即預(yù)測某個要素就采用該要素的時間序列進(jìn)行分析。而風(fēng)向作為記錄風(fēng)吹來方向的氣象要素，如圖2所示，被劃分為16個方位，數(shù)值變化呈圓形循環(huán)，正是由于其360°變化的特征導(dǎo)致了風(fēng)向預(yù)測的困難。例如，0°～22.5°的風(fēng)向和 337.5°～360.0°的風(fēng)向在實際中差別不大，都可稱為北風(fēng)，但從算法的數(shù)值上看相差巨大。

在之前對風(fēng)向預(yù)測的試驗中，嘗試了將標(biāo)識的風(fēng)向作為Label應(yīng)用于機器學(xué)習(xí)，其得出的結(jié)果在實際中不符合邏輯：相差最遠(yuǎn)的標(biāo)簽在實際中卻擁有相距最近的角度。使用時間序列分析，則避免了設(shè)置Label這一尷尬問題，算法分析均是基于歷史數(shù)據(jù)，從變化趨勢得到預(yù)測值。

為了進(jìn)一步解決風(fēng)向預(yù)測的難點并嘗試新的預(yù)測方法，本研究針對風(fēng)向的特殊性加入了風(fēng)速序列，結(jié)合氣象學(xué)中風(fēng)向風(fēng)速的轉(zhuǎn)化公式，將其拆分成U、V風(fēng)并分別進(jìn)行預(yù)測，再將結(jié)果結(jié)合，以避免風(fēng)向數(shù)值差異過大，進(jìn)而與使用單一的風(fēng)向序列預(yù)測結(jié)果進(jìn)行對比。

1.3.1 風(fēng)的U、V分量。傳統(tǒng)的風(fēng)向以東、南、西、北作為標(biāo)識，轉(zhuǎn)化到數(shù)值上為0°表示北風(fēng)、90°表示東風(fēng)。在氣象領(lǐng)域，風(fēng)場是由U、V風(fēng)速分量來組成的二維場[15]。因此，可將傳統(tǒng)的風(fēng)進(jìn)行轉(zhuǎn)換，分別為U風(fēng)和V風(fēng)。U為東西風(fēng)，即用U的正負(fù)來代表風(fēng)為東風(fēng)或者西風(fēng)，絕對值代表此方向上風(fēng)速的大小，U為正，代表傳統(tǒng)意義上的西風(fēng)，為負(fù)則代表東風(fēng)；而V為南北風(fēng)，正為南，負(fù)為北。在風(fēng)速上，則是對傳統(tǒng)風(fēng)進(jìn)行矢量分解。具體公式如下：

式中：S為風(fēng)速；D為風(fēng)向。

1.3.2 風(fēng)向評分。為檢驗風(fēng)向預(yù)測的效果，本文將16個風(fēng)向方位簡化為 8個， 即 0°～22.5°和 337.5°～360.0°記 為 北 風(fēng) 、22.5°～67.5°記 為 東 北 風(fēng) 、67.5°～112.5°記為東風(fēng)，以此類推。這樣既能滿足需求，也簡化了數(shù)據(jù)量，并采用了風(fēng)向的預(yù)報評分FaWD，計算公式為：

式中：SCr為該站點預(yù)報的風(fēng)向預(yù)報得分（表2），例如實況風(fēng)向為90°、預(yù)測風(fēng)向為30°，則得分為0.6，Nf為預(yù)報的總次數(shù)。這樣既可以避免使用傳統(tǒng)預(yù)測中的均方根誤差進(jìn)行評估而導(dǎo)致的結(jié)果混亂，還可以簡化計算過程，在得出結(jié)果后就可以獲得算法的評價。

為了使預(yù)測目標(biāo)更具實際意義，同時更好地模擬冬奧會期間的實況，將風(fēng)向的預(yù)測目標(biāo)定為2019年2月4—20日共17 d每天8：00—17：00共10 h的風(fēng)向，并將每天的預(yù)測得分進(jìn)行平均再比較。

表2 8個風(fēng)向的預(yù)報評分對照

2 結(jié)果與分析

2.1 單時間序列

對4個站點分別進(jìn)行預(yù)測，結(jié)果如圖3所示。將每天的評分進(jìn)行平均（表3），可以看到XGBoost在4個站中都獲得了最好的效果，相對ARMA最高可取得125%的提升效果（A1490），并且XGBoost在連續(xù)17 d中評分曲線較為平滑，相對ARMA在A1489中的評分呈上升趨勢，卻在A1490的2月4—10日中呈現(xiàn)極低的準(zhǔn)確率，表明XGBoost相比傳統(tǒng)的時間序列預(yù)測方法不僅在精度上有了較大提升，也更加穩(wěn)定；LSTM在A1489和A1490明顯優(yōu)于ARMA，但在A1491和A1492的效果并不好，尤其在A1491的2月6日評分更是出現(xiàn)了0。具體分析這一日的預(yù)測結(jié)果，可以看到：此日的真實風(fēng)向集中于230°～290°，屬于設(shè)立的西風(fēng)范圍，而LSTM的預(yù)測結(jié)果包含多個方向，并不集中，從而導(dǎo)致評分為0；反而預(yù)測結(jié)果較為單一的ARMA卻在此站獲得了很好的效果，相對于在A1489和A1490低迷的表現(xiàn)，證明了ARMA方法在單一盛行風(fēng)的站點可以獲得較好的效果，而在復(fù)雜風(fēng)場則效果不佳。

表3 單時間序列下不同算法在各站點的評分

從圖3可以看出，隨著海拔的升高，3種算法的效果越來越貼近，在海拔為2 099.8 m的A1492站點中3種算法的差別很小，折線幾乎重合，證明結(jié)果與海拔存在關(guān)聯(lián)。為進(jìn)一步驗證這種關(guān)系，將2月4—20日每小時的風(fēng)向繪出散點圖（圖4）。可以看出，A1489站點的風(fēng)向分布較為均勻，而A1491和A1492站點則在200°～350°分布較多，即站點盛行西風(fēng)。這與李 炬等[16]在2020年對小海坨冬奧賽場的觀測試驗結(jié)果相符，同時也解釋了3種算法的預(yù)測結(jié)果較為接近的原因：2個站點以西風(fēng)為主導(dǎo)，在此基礎(chǔ)上算法進(jìn)行分析，預(yù)測結(jié)果會更加趨于西風(fēng)，尤其是對ARMA這種非常依賴前期元素的預(yù)測模型。這也為后期的預(yù)測方法提供了思路，即先分析數(shù)據(jù)分布，探明數(shù)據(jù)規(guī)律，在算法預(yù)測基礎(chǔ)上進(jìn)行再次訂正，從而進(jìn)一步提高預(yù)測準(zhǔn)確率。

2.2 雙時間序列

雙時間序列即將風(fēng)向、風(fēng)速序列轉(zhuǎn)化為U、V風(fēng)，使用算法對U、V風(fēng)單獨進(jìn)行預(yù)測，再將其合成為風(fēng)向，預(yù)測結(jié)果如圖5所示。可以看出，采用雙時間序列進(jìn)行預(yù)測后，LSTM在A1491站點的2月6日評分為0的情況得到了解決并且獲得了較高的評分，證明U、V風(fēng)的轉(zhuǎn)換對風(fēng)向的預(yù)測有著積極作用。

將每天的評分進(jìn)行平均得到表4，可以看到：XGBoost在風(fēng)向預(yù)測中依舊擁有很大的優(yōu)勢，在A1489中相對ARMA最大能有115%的提升，4個站點平均后的得分達(dá)到了0.7，具有很好的效果；LSTM的效果優(yōu)于ARMA、遜于XGBoost，這與使用單序列進(jìn)行預(yù)測時相同。

表4 雙時間序列下不同算法在站點的評分

相對于表3的單序列預(yù)測結(jié)果，使用雙序列進(jìn)行預(yù)測，XGBoost和LSTM都有提升，而ARMA基本沒有變化。具體觀察圖3和圖5可以看到，雙序列預(yù)測A1489中ARMA評分在后幾天的表現(xiàn)并不如單序列，但XGBoost在A1489和A1492中的評分相比單序列更為平滑。這表明將風(fēng)速風(fēng)向轉(zhuǎn)化后，傳統(tǒng)的時間序列統(tǒng)計方法過于依賴建模的數(shù)值，并不能分析出數(shù)據(jù)原來的變化規(guī)律，相對于ARMA來說U、V風(fēng)是嶄新的輸入數(shù)據(jù)，而XGBoost和LSTM卻可以采用這一方法得到提升，進(jìn)而驗證了此種方法的可行性。

3 結(jié)論與討論

本文基于山區(qū)4個氣象站點觀測數(shù)據(jù)，結(jié)合時間序列分析方法，使用傳統(tǒng)的ARMA分析、機器學(xué)習(xí)XGBoost方法和神經(jīng)網(wǎng)絡(luò)LSTM對4個站點特定時期的風(fēng)速和風(fēng)向進(jìn)行了預(yù)測。為了獲得更準(zhǔn)確的結(jié)果，嘗試將單一的風(fēng)向風(fēng)速序列轉(zhuǎn)化成U、V風(fēng)，分別預(yù)測，進(jìn)而合成結(jié)果，從而避免風(fēng)向的數(shù)值問題，結(jié)果表明：所使用的3種算法中，XGBoost不論是在預(yù)測精度還是穩(wěn)定度上都優(yōu)于傳統(tǒng)的ARMA分析方法，而神經(jīng)網(wǎng)絡(luò)LSTM在時間序列的預(yù)測中并沒有取得很好的效果，同時訓(xùn)練模型還需要耗費大量時間。風(fēng)向散點圖可以反映山區(qū)海拔較高的站點經(jīng)常盛行一種風(fēng)。由于山谷、山脊地形地勢常形成這種盛行風(fēng)，這也導(dǎo)致了在這幾個站點算法預(yù)測評分的高度重疊。這同時也對風(fēng)向預(yù)測提供了一種思路：在獲得算法的預(yù)測結(jié)果后，可由預(yù)報員根據(jù)當(dāng)?shù)貧夂蛱卣鳎瑢︻A(yù)測結(jié)果進(jìn)行人工訂正，以期獲得更好的預(yù)報效果。將風(fēng)向風(fēng)速轉(zhuǎn)化為U、V風(fēng)并進(jìn)行分析后發(fā)現(xiàn)，XGBoost和LSTM方法都有了提升，證明此種方法對提升風(fēng)向的預(yù)測效果有著積極作用。尤其是XGBoost方法，在風(fēng)的諸多試驗中都表明了其可成為山區(qū)風(fēng)向預(yù)報的可靠方法。