基于時間序列疏系數(shù)模型的太陽輻射年際變化趨勢預(yù)測

賈興斌，宮響

(青島科技大學(xué) 數(shù)理學(xué)院，山東 青島 266061)

太陽能作為清潔的可再生能源，其有效開發(fā)利用有助于人類生存環(huán)境的改善與經(jīng)濟社會的發(fā)展。但由于地表太陽輻射易受氣候變化、大氣污染、日照時長與云量等因素的影響[1]，可利用的太陽能表現(xiàn)出一定的不穩(wěn)定性和不連續(xù)性。研究發(fā)現(xiàn)，從1957年地面太陽輻射觀測網(wǎng)建立以來，20世紀90年代前后，全球大部分區(qū)域地表太陽輻射經(jīng)歷了減少到增加的變化，即先“變暗”后“變亮”[2-3]。不同區(qū)域由于地理環(huán)境和影響太陽輻射的主要因素不同，太陽輻射還存在“振蕩”現(xiàn)象[4-5]。因此，預(yù)測地表太陽輻射的長期變化趨勢，不僅對研究人類活動在全球氣候變化中的作用有重要意義，也可以為新能源利用如光伏電站的建設(shè)提供參考。

國內(nèi)外學(xué)者對于地表太陽輻射的預(yù)測，主要采用基于經(jīng)驗參數(shù)的統(tǒng)計模型以及基于數(shù)理方程的大氣數(shù)值模式。區(qū)別于上述方法，時間序列分析僅以時間為唯一自變量，根據(jù)已有的歷史數(shù)據(jù)對未來進行預(yù)測。其中，自回歸整合移動平均(auto regression integrated moving average，ARIMA)模型是一種經(jīng)典的時間序列分析方法，具有較高的預(yù)測精度。這一模型在經(jīng)濟學(xué)、醫(yī)療衛(wèi)生、氣象等領(lǐng)域已得到廣泛應(yīng)用[6-9]，近年來在地表太陽輻射的預(yù)測研究中也日益受到重視。如張素寧等[10]發(fā)現(xiàn)在地表太陽輻射的逐時預(yù)測中，ARIMA模型優(yōu)于經(jīng)驗?zāi)Ｐ汀un等[11]利用ARMA-GARCH模型，有效擬合出北京和烏魯木齊兩站位地表太陽輻射的月變化。Shadab等[12]基于ARIMA模型，利用印度馬德里34年的遙感地表太陽輻射數(shù)據(jù)較好地預(yù)測了其未來24個月的變化。但目前仍缺乏ARIMA模型在地表太陽輻射年際變化的應(yīng)用研究。

本文利用濟南站1961—2016年地表太陽輻射的年數(shù)據(jù)，初步識別ARIMA模型，通過對模型參數(shù)及殘差序列進行檢驗確定最優(yōu)ARIMA疏系數(shù)模型，并分析預(yù)測未來10年的太陽輻射年際變化。

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

1.1 數(shù)據(jù)來源及處理

山東省是典型的重工業(yè)經(jīng)濟省份，其較為成熟的產(chǎn)業(yè)集群大都集中于能源、化工等傳統(tǒng)領(lǐng)域，而代表未來經(jīng)濟發(fā)展方向的新能源、新材料、節(jié)能環(huán)保等新興產(chǎn)業(yè)卻沒有形成足夠的規(guī)模。同時山東人口眾多，資源能源消耗強度大，隨著城鎮(zhèn)化進程的加快，城市污染日益凸顯，進而導(dǎo)致地面太陽輻射變化[13]。

山東省目前有濟南、莒縣和福山三個國家級輻射觀測站，其中濟南站的觀測序列最長，且濟南市屬于內(nèi)陸城市，環(huán)境污染問題較為突出。因此，本文選取濟南市1961—2016年的地表太陽年總輻射數(shù)據(jù)作為研究對象。該數(shù)據(jù)集下載于國家氣象科學(xué)數(shù)據(jù)中心[14]，其中1978年和1979年的數(shù)據(jù)缺失，采用月輻射數(shù)據(jù)補全，對其中缺失月份(1978年7月、1978年8月、1979年3月)進一步采用日輻射數(shù)據(jù)補全。而日輻射數(shù)據(jù)中也存在部分數(shù)據(jù)缺失，其中1978年7月有12天、8月份有3天，而1979年3月僅有前15天的數(shù)據(jù)。考慮到1978年8月的日數(shù)據(jù)缺失較少，故采用當(dāng)月數(shù)據(jù)均值作為月數(shù)據(jù)，1978年7月份的缺失數(shù)據(jù)采用線性插值獲得，而1979年3月缺測數(shù)據(jù)采用3月份前15天的數(shù)據(jù)與4月份前15天的數(shù)據(jù)取均值計算，最后再結(jié)合其他月份數(shù)據(jù)，求得年輻射值。

多元線性回歸模型所需氣象數(shù)據(jù)(氣溫、降水量、能見度、風(fēng)速等)來自美國國家海洋和大氣管理局國家環(huán)境信息中心[15]。

1.2 時間序列分析模型

ARIMA是時間序列分析中主要用于非平穩(wěn)時間序列分析和預(yù)測的一種較為成熟的分析方法，又稱為Box-Jenkins方法[16]。一般將滿足如下條件的模型簡記為ARIMA(p,d,q)：

(1)

式中，B為延遲算子；Φ(B)為ARIMA(p,d,q)模型的自回歸系數(shù)多項式,Φ(B)=1-φ1B-φ2B2-…-φpBp；Θ(B)為ARIMA(p,d,q)模型的移動平滑系數(shù)多項式,Θ(B)=1-θ1B-θ2B2-…-θ0Bp；at,as為零均值白噪聲序列，E(at)表示t時刻白噪聲序列值的數(shù)學(xué)期望，E(Ys,at)表示s時刻模型預(yù)測值Ys和t時刻白噪聲序列值at的數(shù)學(xué)期望。特別地，當(dāng)d=0，ARIMA(p,d,q)模型實際上是平穩(wěn)時間序列模型ARMA(p,q)；當(dāng)p=0時，ARIMA(p,d,q)模型退化為差分移動平均模型IMA(d,q)；當(dāng)q=0時，ARIMA(p,d,q)模型退化為差分自回歸模型ARI(p,d)。

ARIMA模型實質(zhì)是將自平穩(wěn)時間序列模型ARMA(p,q)和差分運算相結(jié)合，該模型能夠更好地擬合非平穩(wěn)時間序列。如果自相關(guān)和移動平滑部分有缺省，則ARIMA模型可簡寫為:

ARIMA((p1,…,pm),d,(q1,…,qm))

。

(2)

本文采用的時間序列分析軟件為SAS(statistical analysis system)，SAS系統(tǒng)具有全球一流的數(shù)據(jù)倉庫功能，在進行時間序列分析時具有其他統(tǒng)計軟件無可比擬的優(yōu)勢[17-18]。

2 結(jié)果與討論

2.1 時間序列數(shù)據(jù)預(yù)處理

一般地，需要對時間序列的平穩(wěn)性和純隨機性進行檢驗，根據(jù)檢驗結(jié)果，確定要采用的擬合預(yù)測模型。時間序列的平穩(wěn)性檢驗，一般采取時序圖檢驗和構(gòu)造統(tǒng)計量進行假設(shè)檢驗兩種方法。圖1為1961—2016年濟南市太陽年輻射量的時序圖，由圖1可見原時間序列具有明顯的波動性，自1961—1990年呈顯著的下降趨勢，1990—2016年較為平穩(wěn)，但總體呈上升趨勢，因此需進一步進行統(tǒng)計檢驗。

圖1 1961—2016年濟南市太陽年輻射時序圖Fig.1 Time series data of solar radiation at Jinan station during 1961 to 2016

單位根檢驗是構(gòu)造統(tǒng)計量進行序列平穩(wěn)性檢驗最常用的方法，其統(tǒng)計量有很多，ADF(augmenteddickey-Fuller test)檢驗是其中經(jīng)典、簡單的一種，也稱為增廣Dickey-Fuller檢驗。ADF檢驗有三種類型的單位根檢驗?zāi)Ｐ?，具體結(jié)構(gòu)見表1?？梢姡蛄械臋z驗結(jié)果中雖然零均值回歸結(jié)構(gòu)的P值大于顯著性水平0.05，但單均值和趨勢類型中各種延遲模型的Tau統(tǒng)計量(τ)的P值小于顯著性水平0.05，據(jù)此可判斷，該時間序列平穩(wěn)，且該序列的確定性部分可以用常數(shù)均值或趨勢類的各種延遲模型結(jié)構(gòu)進行擬合。也就是說，對濟南市1961—2015年地表太陽輻射序列的擬合與預(yù)測，可采用平穩(wěn)時間序列模型ARMA或帶有趨勢的非平穩(wěn)時間序列模型ARIMA。

進一步對原始序列做一階差分，發(fā)現(xiàn)序列值在0附近波動，呈現(xiàn)出明顯的平穩(wěn)性特征(圖2)，且差分后ADF單位根檢驗值(表1)顯示，三種類型的檢驗?zāi)Ｐ拖娄咏y(tǒng)計量的P值遠小于顯著性水平0.05，這表明濟南市年太陽輻射序列經(jīng)一階差分，消除線性趨勢后為平穩(wěn)序列。

表1 ADF單位根檢驗Table 1 ADF unit root test

圖2 1961—2016年濟南市地表太陽輻射差分時序圖Fig.2 Time series data of differential solar radiation at Jinan during 1961 to 2016

時間序列的白噪聲檢驗一般采用LB統(tǒng)計量(L)，如式(3)所示。

,

(3)

式中n為序列觀測期數(shù),m為延遲期數(shù)。LB統(tǒng)計量近似服從自由度為m的卡方(χ2)分布，同時計算差分前后的LB統(tǒng)計量，檢驗結(jié)果如表2所示。給定顯著性水平α=0.05，各延遲期數(shù)的LB統(tǒng)計量的P值均小于α，判定該序列在差分前后均是非白噪聲序列。結(jié)合平穩(wěn)性檢驗的結(jié)果，我們可以認為濟南市地表太陽輻射原序列與差分后序列均是平穩(wěn)非白噪聲序列。

表2 濟南市年太陽輻射序列的白噪聲檢驗Table 2 White noise examination of annual solar radiation series in Jinan

2.2 時間序列模型的建立

2.2.1 模型的初步識別

對平穩(wěn)非白噪聲序列建模，通過對該序列的樣本自相關(guān)系數(shù)(ACF)和偏自相關(guān)系數(shù)(PACF)的分析，初步確定模型的階數(shù)，即p，q的取值。

首先對地表太陽輻射原序列進行分析。由圖3(a)(b)可見，ACF基本呈指數(shù)衰減，是一種比較典型的拖尾特征，而PACF值延遲一階以后快速減小至2倍標準差范圍以內(nèi)，但在五階時PACF值突然升高至2倍標準差范圍以外，之后又快速減小至0附近，顯示出截尾特征，因此可以初步判斷該模型為AR(5)。

一階差分后時間序列的相關(guān)分析如圖3(c)(d)所示，自相關(guān)系數(shù)ACF值呈現(xiàn)四階截尾，偏自相關(guān)系數(shù)PACF值呈現(xiàn)一階、二階和四階拖尾，因此可初步判定該差分后的時間序列可用ARIMA(4,1,4)擬合序列。

圖3 1961—2015年濟南市年太陽輻射序列自相關(guān)和偏自相關(guān)圖Fig.3 Autocorrelation and partial autocorrelation of annual solar radiation series in Jinan during 1961 to 2015

進一步選擇貝葉斯信息BIC準則，取p∈[0,5]和q∈[0,5]，選取使BIC達到最小的(p,q)組合來分別確定差分前后最優(yōu)的模型階數(shù)，SAS輸出結(jié)果見圖4。對原樣本時間序列，當(dāng)(p,q)=(5,0)時，BIC值為11.04達到最小值(圖4(a))，故最佳擬合模型為ARMA(5,0)模型，即AR(5)模型。對差分后的時間序列，取p∈[0,5]和q∈[0,5]，各(p,q)組合下BIC值結(jié)果見圖4(b)所示。當(dāng)(p,q)=(4,0)時，BIC值為11.13達到最小值，故最佳擬合模型為ARIMA(4,1,0)。

圖4 不同(p, q)組合下的BIC值Fig.4 BIC value with different values of (p, q) in models

2.2.2 疏系數(shù)模型的建立

對建立的AR(5)模型，使用條件最小二乘法對模型參數(shù)進行檢驗，同時考慮到擬合模型殘差的性質(zhì)，對模型進行殘差檢驗，檢驗結(jié)果如表3所示。取α=0.05，參數(shù)φ2的P值大于0.05，未通過檢驗。

同時，對建立的ARIMA(4,1,0)模型參數(shù)和殘差進行統(tǒng)計檢驗，檢驗結(jié)果如表3所示。取α=0.05，殘差檢驗結(jié)果顯示殘差序列為白噪聲序列，參數(shù)顯著性檢驗結(jié)果顯示φ1、φ2和φ4的P值均小于0.05，通過檢驗，但φ3和μ的P值大于0.05，未通過檢驗。如果ARIMA模型中有部分自相關(guān)系數(shù)φj(1≤j

表3 ARIMA(4,1,0)模型參數(shù)檢驗Table 3 ARIMA(4,1,0) model parameter test

2.2.3 疏系數(shù)模型的檢驗

對疏系數(shù)ARIMA((1,2,4),1,0)模型進行參數(shù)顯著性檢驗，結(jié)果如表4所示。根據(jù)條件最小二乘法估計可知P值小于0.05，故ARIMA((1,2,4),1,0)模型的參數(shù)檢驗通過。

表4 ARIMA((1,2,4),1,0)模型參數(shù)檢驗Table 4 ARIMA((1,2,4),1,0) model parameter test

進一步對疏系數(shù)ARIMA((1,2,4),1,0)模型做殘差正態(tài)診斷，結(jié)果如圖5所示，其中核表示擬合的ARIMA((1,2,4),1,0)模型中殘差的核密度函數(shù)曲線。由圖5中的殘差分布圖及其正態(tài)QQ圖知該殘差序列基本呈零均值正態(tài)分布，滿足殘差假定。

圖5 殘差正態(tài)診斷Fig.5 Residual normal diagnosis

綜上，ARIMA((1,2,4),1,0)模型通過檢驗，且修正后的模型為：

,

(4)

其中B為延遲算子，εt為白噪聲序列。

2.3 地表太陽輻射年際變化的分析及預(yù)測

首先對ARIMA((1,2,4),1,0)模型的擬合效果進行驗證，與濟南市1961—2015年的年輻射觀測值相比，模型擬合結(jié)果與觀測值較為吻合(圖6)，除個別年份，如1986年、1992年、2012年擬合值較觀測值偏低，1964—1965年及1985年擬合值略高于觀測值。需要指出的是，部分年份觀測值和擬合值之間相對大小趨勢存在較大差異，如1970、1971、1975、1985、2000—2005等年份，觀測值為極大(或極小)時，擬合值恰為極小(或極大)。分析原因發(fā)現(xiàn)，當(dāng)天氣狀況不太穩(wěn)定，陰雨天氣比較多，降雨量年間變化比較劇烈的年份，地表太陽年總輻射會產(chǎn)生較大波動，模型的擬合效果也較差，而良好天氣狀況下年總輻射比較穩(wěn)定的年份，擬合結(jié)果比較精確，這與文獻[19]研究結(jié)果相似。如1985年降雨量約為708 mm，而1986年降雨量降為344 mm左右；2000—2004年5年平均降雨天數(shù)為108天，遠高于1980—2016年平均降水天數(shù)(78天)?？傮w上，模型觀測值與擬合值的多年平均相對誤差為3.1%，多年平均的均方根誤差約為192 MJ/m2。這表明該模型可用于濟南市年輻射的預(yù)測。

圖6 濟南市1961—2025年ARIMA((1,2,4),1,0)模型擬合地表年輻射值及預(yù)測結(jié)果Fig.6 ARIMA((1,2,4),1,0) modeling results of annual global solar irradiance at Jinan City during 1961 to 2025

模型預(yù)測結(jié)果顯示：2017—2025年濟南市太陽總輻射年平均值約為4 980 MJ/m2，2017—2020年輻射值均高于2021—2025年輻射值(圖6)。為進一步驗證說明預(yù)測結(jié)果的可靠性，本文采用遙感晴天下行總輻射數(shù)據(jù)對比，結(jié)果顯示，ARIMA預(yù)測值的變化趨勢與遙感數(shù)據(jù)較為吻合。

地表太陽總輻射不僅受人類活動、大氣污染、降水、云量與風(fēng)速等因素影響[4-5, 20]，同時也與觀測站位的遷站、海拔高度、周圍遮擋物環(huán)境等諸多因素有關(guān)[21]，且在不同時段主要影響因素也不同，如華東地區(qū)太陽輻射1961—1989年的主要影響因素是氣溶膠，而在1990—1999年和2000—2008年主要影響是云量[1]。早期研究發(fā)現(xiàn)，山東省地表太陽輻射1961—2012年間整體呈下降趨勢[22-25]。王建源等[24]發(fā)現(xiàn)2001—2007年山東省年總輻射比前30年平均減少72.3 MJ/m2。薛德強[26]認為濟南市1961—1990年大氣污染物的增多對地表太陽輻射量的減少起著決定性的作用。本文的數(shù)據(jù)分析顯示，自1992年以來，濟南市地表太陽總輻射呈持續(xù)增大趨勢(圖1)，這可能與濟南站的遷站有一定原因。此外，我們發(fā)現(xiàn)，2012—2016年濟南市地表總輻射平均值較前20年平均增加287 MJ/m2(圖6)，這可能與近年來山東省各項大氣污染防治措施逐步完成并發(fā)揮作用，空氣質(zhì)量得到極大改善有關(guān)。未來，隨著大氣環(huán)境質(zhì)量的進一步改善，地表總輻射整體將繼續(xù)呈增長趨勢(圖6)。

2.4模型對比

為對比分析ARIMA模型預(yù)測效果，本文使用多元線性回歸方法[27-28]，利用平均氣溫(X1)、平均最高氣溫(X2)、平均最低氣溫(X3)、平均露點溫度(X4)、降水量(X5)、最大單日降水量(X6)、降水天數(shù)(X7)、平均能見度(X8)、平均風(fēng)速(X9)、平均站點氣壓(X10)10個變量構(gòu)建線性回歸模型：

Yt=1 613.30X1t-634.92X2t-1 028.40X3t-73.12X4t-0.039X5t+0.89X6t-2.41X7t+97.31X8t

-55.59X9t+30.35X10t-27 379.10

(6)

選擇1980—2015年的太陽輻射數(shù)據(jù)進行檢驗，結(jié)果如圖7所示。整體上，多元線性回歸模型多年平均相對誤差為4.2%，多年平均均方根誤差約為201 MJ/m2。與ARIMA((1,2,4),1,0)模型效果相比，誤差偏大，但就變化趨勢的擬合效果而言，多元線性回歸模型優(yōu)于ARIMA模型，這可能是由于多元線性回歸模型考慮了降雨量等因素。

圖7 濟南市1980—2015年線性模型地表年輻射值及預(yù)測結(jié)果Fig.7 Linear model annual surface radiation values and prediction results during 1980 to 2015 in Jinan City

3 結(jié)語

本文采用ARIMA模型直接預(yù)測的方法，建立ARIMA((1,2,4),1,0)疏系數(shù)模型，預(yù)測了10年的太陽年總輻射量。對比多元線性回歸模型，誤差分析表明該ARIMA模型的預(yù)測精度更高。預(yù)測結(jié)果顯示，2017—2025年地表太陽總輻射量將保持較為平穩(wěn)的增長趨勢，可加強對太陽能資源的利用。