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

賈興斌，宮響

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

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

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

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

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

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

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

山東省目前有濟(jì)南、莒縣和福山三個(gè)國(guó)家級(jí)輻射觀測(cè)站，其中濟(jì)南站的觀測(cè)序列最長(zhǎng)，且濟(jì)南市屬于內(nèi)陸城市，環(huán)境污染問題較為突出。因此，本文選取濟(jì)南市1961—2016年的地表太陽(yáng)年總輻射數(shù)據(jù)作為研究對(duì)象。該數(shù)據(jù)集下載于國(guó)家氣象科學(xué)數(shù)據(jù)中心[14]，其中1978年和1979年的數(shù)據(jù)缺失，采用月輻射數(shù)據(jù)補(bǔ)全，對(duì)其中缺失月份(1978年7月、1978年8月、1979年3月)進(jìn)一步采用日輻射數(shù)據(jù)補(bǔ)全。而日輻射數(shù)據(jù)中也存在部分?jǐn)?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月缺測(cè)數(shù)據(jù)采用3月份前15天的數(shù)據(jù)與4月份前15天的數(shù)據(jù)取均值計(jì)算，最后再結(jié)合其他月份數(shù)據(jù)，求得年輻射值。

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

1.2 時(shí)間序列分析模型

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

(1)

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

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

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

。

(2)

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

2 結(jié)果與討論

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

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

圖1 1961—2016年濟(jì)南市太陽(yáng)年輻射時(shí)序圖Fig.1 Time series data of solar radiation at Jinan station during 1961 to 2016

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

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

表1 ADF單位根檢驗(yàn)Table 1 ADF unit root test

圖2 1961—2016年濟(jì)南市地表太陽(yáng)輻射差分時(shí)序圖Fig.2 Time series data of differential solar radiation at Jinan during 1961 to 2016

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

,

(3)

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

表2 濟(jì)南市年太陽(yáng)輻射序列的白噪聲檢驗(yàn)Table 2 White noise examination of annual solar radiation series in Jinan

2.2 時(shí)間序列模型的建立

2.2.1 模型的初步識(shí)別

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

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

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

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

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

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

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

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

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

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

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

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

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

進(jìn)一步對(duì)疏系數(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)模型通過檢驗(yàn)，且修正后的模型為：

,

(4)

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

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

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

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

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

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

2.4模型對(duì)比

為對(duì)比分析ARIMA模型預(yù)測(cè)效果，本文使用多元線性回歸方法[27-28]，利用平均氣溫(X1)、平均最高氣溫(X2)、平均最低氣溫(X3)、平均露點(diǎn)溫度(X4)、降水量(X5)、最大單日降水量(X6)、降水天數(shù)(X7)、平均能見度(X8)、平均風(fēng)速(X9)、平均站點(diǎn)氣壓(X10)10個(gè)變量構(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年的太陽(yáng)輻射數(shù)據(jù)進(jìn)行檢驗(yàn)，結(jié)果如圖7所示。整體上，多元線性回歸模型多年平均相對(duì)誤差為4.2%，多年平均均方根誤差約為201 MJ/m2。與ARIMA((1,2,4),1,0)模型效果相比，誤差偏大，但就變化趨勢(shì)的擬合效果而言，多元線性回歸模型優(yōu)于ARIMA模型，這可能是由于多元線性回歸模型考慮了降雨量等因素。

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

3 結(jié)語(yǔ)

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