王志信，黃鵬，林友明，賈秀鵬

(1.中國科學院大學，北京 100049；2.中國科學院遙感與數(shù)字地球研究所，北京 100094)

1 引 言

為滿足全球礦產(chǎn)資源、作物估產(chǎn)、水資源、森林生態(tài)和城市變化檢測等應(yīng)用領(lǐng)域的特定數(shù)據(jù)獲取需求，遙感技術(shù)特別是光學遙感技術(shù)在其中發(fā)揮了重要作用[1-4]。在光學遙感數(shù)據(jù)獲取應(yīng)用中，對于特定目標區(qū)域，如果在云量覆蓋較多的時候進行衛(wèi)星成像規(guī)劃，則容易得到利用價值較低的多云數(shù)據(jù)，同時也會造成衛(wèi)星資源的浪費，而如果能提前預(yù)知云量覆蓋情況，選擇云量覆蓋較少的時段進行衛(wèi)星成像規(guī)劃，則更有利于獲取到滿足遙感應(yīng)用需求的少云或無云數(shù)據(jù)。因此，如何進行云量預(yù)測，以獲取優(yōu)質(zhì)的能滿足遙感應(yīng)用需求的無云或少云數(shù)據(jù)，便成為了一個重要問題。

然而，目前還沒有直接面向衛(wèi)星成像規(guī)劃的云量預(yù)測方法。在云量預(yù)測相關(guān)模型方法方面，應(yīng)用較多的是與氣象相關(guān)的短期預(yù)測模型，此類模型主要是利用一些模型方法，如模式識別、交叉相關(guān)法等來對云變化進行短期預(yù)測[5-6]，預(yù)測期通常為幾小時。例如，利用衛(wèi)星云圖對一個熱帶颶風進行為期幾小時的短期預(yù)測。而一般衛(wèi)星成像規(guī)劃所需要的預(yù)測時間長度要遠大于幾小時，此類短期預(yù)測模型對于衛(wèi)星成像規(guī)劃來說實際應(yīng)用價值不大。還有一類有代表性的模型即人工神經(jīng)網(wǎng)絡(luò)模型(Artificial Neural Network，ANN)，這是一種應(yīng)用于對大腦神經(jīng)突觸聯(lián)接的結(jié)構(gòu)進行信息處理的數(shù)學模型。它具有很多優(yōu)點，如自學習功能，這對于預(yù)測有特別重要的意義，還具有聯(lián)想存儲功能、具有高速尋找優(yōu)化解的能力等[7-8]。雖然神經(jīng)網(wǎng)絡(luò)模型的這些優(yōu)點使它可以應(yīng)用于云量預(yù)測，但是由于影響云分布的因素眾多，且衛(wèi)星規(guī)劃的區(qū)域范圍不定，大到國家小到城市都有可能，這中間牽涉到的眾多因素將導致神經(jīng)網(wǎng)絡(luò)中節(jié)點的確定和學習過程非常復雜，因而，很難建立一個能夠全面描述并準確預(yù)測大范圍云量的模型。這使得神經(jīng)網(wǎng)絡(luò)模型也不太適用于衛(wèi)星數(shù)據(jù)獲取所要求的云量預(yù)測。

為了解決傳統(tǒng)云量預(yù)測模型不能滿足遙感實際應(yīng)用需求的問題，本文采用時間序列分析預(yù)測的方法，對云量進行預(yù)測。首先，針對云量數(shù)據(jù)本身具有時間序列的特點，對云量進行分類。然后，在分類的基礎(chǔ)上對不同的類型采取不同的模型進行預(yù)測，最終得到云量預(yù)測結(jié)果，從而避免考慮眾多的云量影響因素導致建模困難的問題，利用云量預(yù)測結(jié)果信息可為衛(wèi)星成像規(guī)劃提供參考，有利于更好地獲取無云或少云數(shù)據(jù)，減少衛(wèi)星資源的浪費，滿足遙感數(shù)據(jù)應(yīng)用需求具有重要的現(xiàn)實意義。

2 時間序列分析建模方法

1927年，英國統(tǒng)計學家G.U.Yule(1871～1951)提出了自回歸(Auto Regressive，AR)模型。不久之后，英國數(shù)學家、天文學家G.T.Walker爵士在分析印度大氣規(guī)律時使用了滑動平均(Moving Average，MA)模型和自回歸滑動平均(Auto Regressive Moving Average，ARMA)模型(1931年)。這些模型奠定了時間序列時域分析方法的基礎(chǔ)。1970年，美國統(tǒng)計學家G.E.P.Box和英國統(tǒng)計學家G.M.Jenkins聯(lián)合出版了《時間序列分析—預(yù)測與控制》一書(文獻[9])。在書中，Box和Jenkins在總結(jié)前人研究的基礎(chǔ)上，系統(tǒng)地闡述了對求和自回歸滑動平均(Autoreg Ressive Integrated Moving Average，ARIMA)模型的識別、估計、檢驗及預(yù)測的原理及方法。這些時間序列分析的方法，已經(jīng)成為時域分析方法的核心內(nèi)容。

2.1 ARMA模型

自回歸滑動平均模型(Auto Regression Moving Average，ARMA)是目前最常用的擬合平穩(wěn)序列的模型。其模型結(jié)構(gòu)如下：

(1)

ARMA(p，q)模型可以簡寫為：

xt=φ1xt-1+φ2xt-2+…+φpxt-p+εt-θ1εt-1-θ2εt-2-…-θqεt-q

(2)

引進延遲算子，ARMA(p，q)模型簡記為Φ(B)xt=Θ(B)εt。其中，Φ(B)=1-φ1B-φ2B2-…-φpBp，為p階自回歸系數(shù)多項式。Θ(B)=1-θ1B-θ2B2-…-θqBq，為q階滑動平均系數(shù)多項式。當q=0時，ARMA(p，q)模型就退化成了AR(p)模型；當p=0時，ARMA(p，q)模型就退化成了MA(q)模型。所以，AR(p)模型和MA(q)模型實際上是ARMA(p，q)模型的特例，它們都統(tǒng)稱為ARMA模型。

2.2 ARIMA模型

ARMA(p，q)模型描述的是平穩(wěn)時間序列，然而在實際工作中遇到的很多都是非平穩(wěn)時間序列。為了將其平穩(wěn)化，通常利用一次差分或多次差分運算把原時間序列轉(zhuǎn)化成平穩(wěn)時間序列。一個時間序列如果能通過差分運算轉(zhuǎn)化為平穩(wěn)時間序列，則稱這個非平穩(wěn)時間序列為差分平穩(wěn)時間序列。對差分平穩(wěn)時間序列可以使用ARIMA模型進行擬合。

2.2.1 差分運算

2.2.2 ARIMA模型

具有如下結(jié)構(gòu)的模型稱為求和自回歸移動平均(Autoreg Ressive Integrated Moving Average)模型，簡記為ARIMA(p，d，q)模型：

(3)

(4)

2.3 SARIMA模型

在某些時間序列中，存在著明顯的周期性變化。這種周期性是由于季節(jié)性變化(包括季度、月度、周度等變化)或其他一些固有因素引起的，這類序列稱為季節(jié)性序列。描述這類序列，本文用季節(jié)時間序列模型(Seasonal ARIMA model，SARIMA)表示。當序列具有短期相關(guān)性時，通?？梢允褂玫碗AARMA(p，q)模型提取。當序列具有季節(jié)效應(yīng)，季節(jié)效應(yīng)本身還具有相關(guān)性時，季節(jié)相關(guān)性可以使用以周期步長為單位的ARMA(P，Q)模型提取。由于短期相關(guān)性和季節(jié)效應(yīng)之間具有乘積關(guān)系，所以擬合模型實質(zhì)為ARMA(p，q)和ARMA(P，Q)的乘積。綜合前面d階趨勢差分和D階以周期S為步長的季節(jié)差分運算，對原觀察值序列擬合的乘積模型完整的結(jié)構(gòu)如下：

(5)

式中：Θ(B)=1-θ1B-…-θqBq，ΘS(B)=1-θ1BS-…-θQBQS，Φ(B)=1-φ1B-…-φpBp，ΦS(B)=1-φ1B-…-φPBPS。該乘積模型簡記為ARIMA(p,d,q)×(P,D,Q)S。

3 實驗方法設(shè)計

3.1 數(shù)據(jù)源

數(shù)據(jù)源采用ISCCP的D2云數(shù)據(jù)集。國際衛(wèi)星云氣候研究計劃(ISCCP)成立于1982年，作為世界氣候研究計劃(WCRP)之一，它主要通過收集和分析氣象衛(wèi)星輻射測量值來研究全球云分布、屬性以及每日、每季和每年的變化。其生產(chǎn)的數(shù)據(jù)集不僅可以用來研究云在氣候變化中所扮演的角色，還可以用來研究云在輻射能量交換中所起的作用以及它在全球水循環(huán)中所起的作用[10-11]。ISCCP提供的云量數(shù)據(jù)分為不同的類別，主要包括：B3和BT級別，大氣數(shù)據(jù)，海冰和雪數(shù)據(jù)，DX、D1、D2級別。其中，D2數(shù)據(jù)是所有級別數(shù)據(jù)當中具有最優(yōu)質(zhì)云量信息的數(shù)據(jù)，在云量研究中應(yīng)用廣泛，因而數(shù)據(jù)源選擇D2數(shù)據(jù)集。

3.2 實驗方法

3.2.1 數(shù)據(jù)讀取及分類

讀取ISCCP的D2數(shù)據(jù)集中與云量有關(guān)的數(shù)據(jù)并提取云量網(wǎng)格的歷史云量特征，對全球云量網(wǎng)格進行分類按照ISCCP數(shù)據(jù)集的說明讀取云氣候?qū)W數(shù)據(jù)中的月平均云量數(shù)據(jù)后，針對云量數(shù)據(jù)的時空分布各異的特征首先對全球云量網(wǎng)格進行分類，在此基礎(chǔ)上針對不同特征進行相應(yīng)的建模預(yù)測。

首先，判斷每個網(wǎng)格的月平均云量時間序列的平穩(wěn)性。判斷某個序列Xi的平穩(wěn)性主要是通過單位根檢驗，判斷標準為ADF檢驗(Augmented Dickey-Fuller)。構(gòu)造ADF檢驗統(tǒng)計量：是參數(shù)ρ的樣本標準差。當τ 大于DW臨界值，則序列存在單位根，是非平穩(wěn)序列。反之則為平穩(wěn)序列。在此基礎(chǔ)上，對于非平穩(wěn)序列，進一步判斷序列是否具有季節(jié)性；主要判斷標準為自相關(guān)函數(shù)ACF(AutoCorrelation function)。自相關(guān)函數(shù)ρt,sρt,s=Corr(xt,xs),t,s=0,±1,±2,…，其中，{xt,t=0,±1,±2,…}是隨機過程。具有季節(jié)性的云量時間序列的ACF具有季節(jié)周期為12的周期性特征，即若某序列的ACF在滯后期12的整數(shù)倍出現(xiàn)峰值，則該序列存在季節(jié)性特征，反之，則序列不存在季節(jié)性特征。

其次，根據(jù)分類結(jié)果選擇合適的模型。采用分類樹分類方法對全球云量網(wǎng)格進行分類：對時間序列是否平穩(wěn)進行判斷；進一步對時間序列是否存在季節(jié)性進行判斷，將全球云量網(wǎng)格分為以下3類：

分類①：月平均云量時間序列平穩(wěn)類型；

分類②：月平均云量時間序列不平穩(wěn)，同時具備季節(jié)周期性變化，即月平均云量時間序列季節(jié)非平穩(wěn)類型；

分類③：月平均云量時間序列不平穩(wěn)，不具備季節(jié)周期性變化，即月平均云量時間序列普通非平穩(wěn)類型。

3.2.2 選取模型

針對不同的云量類型，選取相應(yīng)的模型進行預(yù)測：

對于分類①，選擇建立平穩(wěn)時間序列模型——ARMA模型；

對于分類②，選擇建立非平穩(wěn)季節(jié)性時間序列模型——SARIMA模型；

對于分類③，選擇建立非平穩(wěn)時間序列模型——ARIMA模型。

3.2.3 數(shù)據(jù)預(yù)處理

選取模型之后，根據(jù)不同的模型要求，對歷史云量數(shù)據(jù)進行相應(yīng)的處理。ARMA模型可以直接應(yīng)用原始數(shù)據(jù)；ARIMA模型則需要在結(jié)合檢驗平穩(wěn)性的基礎(chǔ)上，對數(shù)據(jù)進行差分處理；季節(jié)性SARIMA模型需要結(jié)合檢驗平穩(wěn)性的基礎(chǔ)上，根據(jù)需要對數(shù)據(jù)進行季節(jié)差分和普通差分處理。

3.2.4 參數(shù)估計

根據(jù)歷史云量數(shù)據(jù)和相應(yīng)的模型，需要對模型進行參數(shù)估計。一般情形下，參數(shù)估計分兩步工作，第1步，找出參數(shù)的初步估計或稱參數(shù)的初始估計，利用樣本自相關(guān)函數(shù)對模型參數(shù)進行初步估計，例如對ARMA模型的特例AR(P)模型具體可利用Yule-Walker方程進行參數(shù)估計。第2步，在初步估計的基礎(chǔ)上，按照一定的估計準則，求得模型參數(shù)的精細估計。模型參數(shù)的估計計算復雜，需要由計算機完成。

3.2.5 模型診斷

為確定模型擬合是否充分，對殘差進行分析，檢驗殘差的同方差性、正態(tài)性和獨立性，檢驗殘差是否接近白噪聲，同時考慮過度擬合和參數(shù)冗余，對模型做出相應(yīng)修正。

3.2.6 預(yù)測

基于最小化的均方預(yù)測誤差方法，利用第i個網(wǎng)格云量數(shù)據(jù)和相應(yīng)模型對云量進行預(yù)測，得出該網(wǎng)格未來數(shù)月的云量預(yù)測結(jié)果數(shù)據(jù)。對于某序列Xi來說可獲得直到時間t的歷史數(shù)據(jù)，即x1，x2，…xT-1，xT，預(yù)測未來e期的值xe+t，稱時間t為預(yù)測起點，e為預(yù)測前置時間，用xt(e)代表預(yù)測值，最小均方誤差預(yù)測如下：xt(e)=E(xe+t|x1，x2，…xT-1，xT)。它的原理是用x1，x2，…xT-1，xT(為方便起見記為y)的函數(shù)來預(yù)測x，標準是最小化均方誤差，這里需要選擇函數(shù)h(y)，使得下式達到最?。篍[x-h(y)]2根據(jù)概率論條件期望性質(zhì)，上述式子可寫為：E{[x-(y)]2|Y=y}=E{[x-h(y)]2Y=y}，因此對于每個y，h(y)的最優(yōu)選擇都是h(y)=E(x|Y=y)，因而h(y)的這一選擇可使得E[x-h(y)]2整體達到最小，因此，h(y)=E(x|Y)是基于所有y的函數(shù)所得的x的最優(yōu)預(yù)測[12]。

4 實驗及結(jié)果

4.1 2007年平均云量預(yù)測實驗

選取北京地區(qū)為感興趣區(qū)，在時間上選取1984年至2006年間22年的月平均云量數(shù)據(jù)。利用此歷史云量數(shù)據(jù)對2007年12個月的月平均云量進行預(yù)測，得到12個月份的月平均云量預(yù)測值并與2007年份真實云量值進行對比驗證。

提取北京地區(qū)所在網(wǎng)格的月平均云量時間序列；判斷月平均云量時間序列的平穩(wěn)性。

圖1 1984年到2006年北京地區(qū)云量序列圖

圖2 北京地區(qū)云量序列自相關(guān)函數(shù)圖

北京地區(qū)云量序列的ACF如圖2所示。從圖中可知該序列是非平穩(wěn)的，并且具有明顯的季節(jié)性特征。采用ADF檢驗，h=0檢驗的P值：pValue=0.3543，樣本統(tǒng)計量：stat=-0.8112，拒絕臨界值：cValue=-1.9417，由此可見屬于第2類，模型識別為SARIMA模型。季節(jié)性SARIMA模型需要在檢驗平穩(wěn)性的基礎(chǔ)上，根據(jù)云量數(shù)據(jù)本身特點對數(shù)據(jù)進行一次季節(jié)差分和一次普通差分處理，采用SARIMA模型進行分析預(yù)測后得出預(yù)測值，并跟真實值進行驗證。

如圖3所示，用1984年到2006年22年數(shù)據(jù)預(yù)測2007年12個月云量。得出2007年12個月份的預(yù)測值，預(yù)測值的80%、95%置信區(qū)間。與2007年實際值比較，實際真實值均落在80%置信區(qū)間內(nèi)，平均相對誤差3.60%，最大相對誤差5.29%，均方誤差3.2199。

表1 北京地區(qū)2007年預(yù)測值與實際值

圖3 利用1984年到2006年22年數(shù)據(jù)預(yù)測2007年12個月份的預(yù)測圖

4.2 1997年平均云量預(yù)測

仍選取北京為感興趣區(qū)，但在時間上選取1984年到1996年的12年月平均云量數(shù)據(jù)，對1997年進行預(yù)測值并與1997年實際值進行比較，結(jié)果顯示，9月份實際值落在95%置信區(qū)間，其余月份實際值均落在80%置信區(qū)間內(nèi)，平均相對誤差4.18%，最大相對誤差15.29%，均方誤差4.2902。

表2 北京地區(qū)1997年預(yù)測值和實際值

圖4 利用1984年到1996年12年數(shù)據(jù)預(yù)測1997年12個月份的預(yù)測圖

4.3 2009年平均云量預(yù)測實驗

選取烏魯木齊為感興趣區(qū)，在時間上選取1984年到2008年間24年月平均云量數(shù)據(jù)對2009年份12個月份月平均云量進行預(yù)測和驗證，結(jié)果顯示，2009年12個月份的預(yù)測值，預(yù)測值的80%、95%置信區(qū)間。與2009年實際值比較，實際值除一月份落在95%置信區(qū)間外，其余月份均落在80%置信區(qū)間內(nèi)，平均相對誤差2.80%，最大相對誤差6.55%，均方誤差2.8871。

表3 烏魯木齊地區(qū)2009年預(yù)測值與實際值

圖5 利用1984年到2008年24年數(shù)據(jù)預(yù)測2009年12個月份的預(yù)測圖

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

由以上實驗結(jié)果可知，絕大部分云量實際值落在云量預(yù)測80%置信區(qū)間之內(nèi)，而且，不論是利用10年還是20年的歷史數(shù)據(jù)，均能得到精度較好的結(jié)果，說明此方法的時間適應(yīng)性較好。更為重要的是，云量預(yù)測結(jié)果較準確地反映了云量的變化趨勢，這對于衛(wèi)星成像規(guī)劃具有很重要的意義，利用云量變化趨勢信息，可在長周期衛(wèi)星成像規(guī)劃中發(fā)揮重要參考作用。以北京為例，根據(jù)云量預(yù)測參考信息，要想獲取此地區(qū)的少云數(shù)據(jù)，應(yīng)該盡量選取云量低的7、8和9月份來進行規(guī)劃，相比在12、1和2月則更容易獲取到滿足遙感應(yīng)用需求的數(shù)據(jù)。

5 結(jié)束語

針對目前沒有能夠直接滿足遙感數(shù)據(jù)獲取需求相關(guān)的云量預(yù)測方法的問題，本文提出了一種利用時間序列分析預(yù)測方法對云量進行預(yù)測的方法。在云量特征分類的基礎(chǔ)上，整合ARMA、ARIMA和SARIMA 3種模型對云量進行預(yù)測，并得到了滿意的預(yù)測結(jié)果。根據(jù)云量預(yù)測結(jié)果參考信息，選擇云量覆蓋較少的時間段進行衛(wèi)星成像規(guī)劃，有利于更合理的進行衛(wèi)星成像規(guī)劃，對滿足遙感數(shù)據(jù)獲取需求有重要意義。由于云量本身的復雜性以及當前云量預(yù)測研究還處在初步研究階段，本方法仍存在一定的不足，如為了更好地滿足衛(wèi)星成像規(guī)劃的實際需求，以月為單位的預(yù)測結(jié)果在時間分辨率上還有待提高，更精確的預(yù)測方法還有待進一步研究。

參考文獻：

[1] SANKEY J B，WALLACEB C S，RAVIC S.Phenology-based remote sensing of post-burn disturbance windows in rangelands[J].Ecological Indicators，2013(30)：35-44.

[2] TULBURE M G，BROICH M.Spatiotemporal dynamic of surface water bodies using Landsat time-series data from 1999 to 2011[J].Journal of Photogrammetry and Remote Sensing，2013，79(2013)：44-52.

[3] JAWAK S D，LUIS A J.Improved land cover mapping using high resolution multiangle 8-band WorldView-2 satellite remote sensing data[J].Journal of Applied Remote Sensing，2013，7(1)：73-86.

[4] HANSSAN Q K，RAHMAN K M.Remote sensing-based determination of understory grass greening stage over boreal forest[J].Journal of Applied Remote Sensing，2013，7(1)：32-45.

[5] ENDLIEH R M，WOLF D E，HALL D L，et a1.Use of a pattern recognition technique for determining cloud motions from Sequences of Satellite photographs[J].Journal of Applied Meteorology，1971，10(2)：105-117.

[6] THOMASM H，THOMASM N.A short-term cloud forecasts scheme using cross correlation[J].Weather and Forecasting，1993，8(4)：401-411.

[7] LIU H，TIANA H Q，PANA D，et al.Forecasting models for wind speed using wavelet，wavelet packet，time series and Artificial Neural Networks[J].Applied Energy，2013，107(C)：191-208.

[8] MALIK M H，ARIF A F.ANN prediction model for composite plates against low velocity impact loads using finite element analysis[J].Composite Structures，2013，101(2013)：290-300.

[9] BOX G，JENKINS G M.Time series analysis：forecasting and control[M].Beijing：China Statistics Press，1997：101-135.

[10] HAN Q Y，ROSSOW W B，LACIS A A.Near-global survey of effective cloud droplet radii in liquid water clouds using ISCCP data[J].Journal of Climate，1994，7(1)：465-497.

[11] HAN Q，ROSSOW W B，CHOU J，et al.Global variation of droplet column concentration in low-level clouds[J].Geophysical Research Letters，1998，25(9)：1419-1422.

[12] JONATHAN D C，CHAN K S.Time series analysis with applications in R[M].2nd ed.Beijing：China Machine Press，2011：137-152.