基于廣義高斯過程回歸的北京市近五十年能見度分析

李海蓉，朱曉欣，曹春正

（南京信息工程大學 數學與統計學院，南京 210044）

0 引言

大氣能見度的降低逐漸成為氣候和空氣污染研究領域的關鍵問題，引起了政府部門和公眾的高度重視。目前，大氣能見度的研究主要集中在影響因素分析、趨勢研究和預測推斷上。需要指出的是，雖然大氣能見度明確地受到氣象條件和污染物狀況的影響，但這種復雜的影響關系很難通過機理分析或常規的統計分析方法獲取。Lin等[1]研發了一種研究能見度、空氣質量及氣象條件間因果關系的最優經驗回歸模型；Yang等[2]通過PLAM/h指標建立了空間大尺度上低能見度的對數線性模型。此類方法的明顯缺陷是忽視了數據的序列相關性以及變量之間的非線性相關性。各種時間序列模型可以一定程度上刻畫這些特征。例如，Pedrycz等[3]運用模糊時間序列模型對上海的大氣能見度進行了預測；Yang等[4]研究了霧霾變化的影響因子及作用機制，并運用時間序列分析提出一種對北京霧霾的長期預測模型。另外一種處理時空數據的常用方法是通過預報方程進行數值模擬，比如CALPUFF模式[5,6]。但是數值預報方法存在幾個問題：首先，預報方程的建立需要機理支撐，這對于很多情況難以實現；其次，由于預報方程的動態性和復雜性，當數據信息不夠充分時，預報的準確性不足；另外，模式方法往往還會在數據空間尺度范圍上受到限制。支持向量機也是數值預報的常用方法[7]。但是單純使用支持向量機方法進行預報容易出現空報、漏報情況。也有文獻將支持向量機和粒子群算法結合對能見度進行預測[8]，但只適用于低能見度情形，且預測結果不穩健。

鑒于能見度的數據格式以及時空結構關系，本文基于廣義高斯過程回歸（GGPFR）[9]，從函數型數據分析[10]的角度對能見度進行研究。GGPFR模型能夠針對非高斯型函數數據，通過聯系函數和高斯過程，刻畫響應變量與時間以及其他多維函數型協變量之間的非線性相關性，本質上屬于一種非參數回歸方法。另外，相比單一站點的預報模式，該模型能夠對多個站點數據進行整體建模，充分挖掘數據信息之間的關聯性。為了獲得更好的擬合和預測效果，本文首先在模型上做出改進，建立帶有混合效應的GGPFR（M-GGPFR）模型，充分利用參數結構的可釋性和非參數結構的靈活性。本文首先闡述模型結構和估計、預測方法，并通過數值模擬研究展示推斷效果，然后對我國20個城市進行建模，并重點分析北京市近五十年的能見度變化規律和趨勢。

1 M-GGPFR模型

1.1 模型介紹

令{zm（t），t∈T}，m=1，2，…，M為第m批數據的函數型響應變量，服從指數族分布：

其中，αm（t） 和?m（t）都是函數型變量，αm（t）是點則參數，?m（t）為散度參數。

令xm（t）為一組函數型協變量，t表示時間或其他一維協變量。由此，函數型響應變量zm（t）與函數型協變量之間的關系刻畫如下：

其中，βm（t）采用B樣條近似估計，vm（t）、wm（t）和xm（t）為函數型協變量的子集，維數分別為r、k、Q。τm（t）為潛在隨機過程，通過高斯過程回歸（GPR）訓練學習，協方差核函數采用平方指數形式，即

其中，θ=（vc，wc1，…，wcq）為一組超參數，可以通過訓練數據集進行估計。

1.2 參數估計

設訓練集共有M批數據，其中第m批數據觀測的時間節點為為方便起見，記zm（tmi）、τm（tmi） 和xm（tmi） 為zmi、τmi和xmi，記Zm= （zm1，… ，zmNm）T，由此給出離散數據下的對應模型：

其中，EF（，）指代（1）式的指數族分布 ，這里為協方差矩陣。結構上，上述模型比GGPFR模型多了混合效應成分，其參數估計的具體實施仍可采用極大似然估計結合Laplace近似方法，詳見文獻[10]。

1.3 預測

以下考慮兩種不同類型的預測問題。

（1）若預測對象來自訓練樣本，即已獲得預測對象的部分數據，目的是預測未來時刻或缺失時刻數據。記在時刻{tk1，…，tkN}的觀測值為Zk={zki，i=1，…，N} ，對應的輸入向量為xk={xk1，…，xkN} 、vk（t）={vk1（t），…，vkN（t）}和wk（t）={wk1（t），…，wkN（t） } ，希望預測新時刻t*下的z*。令τ*=τk（t*）為t*時的潛在過程值，則給定τ*下z*的條件期望為：

于是有：

而Var（z*|D）可由下式計算

（2）若預測對象為全新樣本，即沒有該樣本的響應變量信息，需要給出t*時刻下預測z*，上述方法不再適用。此時可以考慮一種加權預測的思想：先將該樣本視為訓練集中的已知樣本，由此獲得預測值，再通過z*屬于第m樣本的概率ωm進行如下加權預測：

如果沒有信息推斷ωm時，可以取ωm=1/M，m=1，…，M。

2 模擬研究

通過數值模擬檢測模型估計和預測效果。數據產生自：

其中，對于每個m，xmi（=tmi）為將區間[- 4，4] 等分所得，vmi為取值為0或1的列向量，γ=1.0，η=100分別為vmi和wmi的系數。τmi是均值為0，協方差為（3）式所示的高斯過程，其中vc=0.5，wcq=0.1。響應變量zmi由（12）式生成：

由此生成40條曲線的樣本，每條曲線包含41個樣本點。為方便運算且避免協方差矩陣出現奇異，在模型最后加入一項噪聲擾動。一般來講，對于n類有序數組，若定義模型：

其中，b0=-∞，ba=∞，j=1，…，n-1為需要估計的閾值的密度函數可由下式得出：

進而運用經驗貝葉斯方法估計出B樣條系數、函數型協變量系數、超參數和閾值。

模擬部分n=4且閾值b12和b3均為未知參數。圖1（a）展示了M-GGPFR模型對均值曲線的估計（選取第40條曲線進行展示），實線為真實曲線，虛線為估計曲線，表明估計有效。圖1（b）比較了GGPFR模型（點線）與M-GGPFR模型（虛線）對β（t）（實線）的估計結果，結果表明混合效應模型估計更為準確。

圖1 模型估計效果：（a）均值的擬合,（b）兩種模型對 β（t）的擬合

進一步考慮預測問題，首先生成一條由40個點組成的全新曲線，隨機選取一部分觀測值估計模型，剩下的部分作為測試集進行預測（內含預測），同時比較兩種模型的預測情況，如圖2所示。不難發現，M-GGPFR預測效果明顯優于GGPFR模型。

圖2內含預測效果比較

圖2 中實線（星型）為真實值，虛線（空圈）為預測值。其中：（a）GGPFR模型對ym（t） 的預測；（b）GGPFR模型對zm（t）的預測及95%的置信區間；（c）M-GGPFR模型對ym（t） 的預測；（d）M-GGPFR模型對zm（t）的預測及95%的置信區間。

接下來，考慮外延預測問題，選取[-4,0]時刻對應數據作為觀測值，預測（0,4]時刻部分。圖3展示兩種模型的預測效果。相比而言，混合GGPFR模型優于GGPFR模型。

圖3 外延預測效果比較（記號同前）

3 實例分析

本文分析中國20個城市1960—2012年的日觀測能見度數據，包括能見度、降雨、風速和相對濕度。由于測量標準不同，能見度在1980年之前按0～9等級[11]記錄，而1980年后記錄具體的能見度值，2010年開始采取0～6級標準（見表1）。

表1 能見度等級和劃分標準

表2給出樣本下的能見度、風速、相對濕度及降雨之間相關系數，表明能見度與風速呈正相關，與相對濕度呈負相關。由此，將能見度作為響應變量，風速作為固定效應，降雨為隨機效應，而時間和相對濕度用于訓練（3）式中的協方差核函數。

表2 能見度、風速、相對濕度及降雨之間的相關系數

本文運用1981—2002年20個城市1月份的月平均數據數據作為訓練集，以此預測后十年的能見度。下頁圖4是不同方法下北京市能見度趨勢預測效果，包括普通線性回歸、ARMA模型和M-GGPFR模型預測。

圖4 北京市能見度趨勢預測比較（單位：0.1km）

實線和星型表示真實測量值，虛線為預測值；實心圓表示M-GGPFR模型預測，三角形為時間序列預測，正方形為線性回歸預測。

觀察預測曲線并比較三者的均方根誤差（RMSE）發現，M-GGPFR模型有最好的預測效果，其均方根誤差為0.12，遠小于時間序列模型的61.57和線性回歸的72.02。

接下來，選取北京市1969—1989年冬季的數據展開研究。參照Zhang等[12]提出的能見度評價體系，日能見度20km以上視為“極好”天氣，10～20km為“好”天氣，2～10km為“壞”天氣，低于2km為“特壞”天氣。將北京市1969—1989年冬季的能見度進行頻率劃分，見下頁圖5。結果表明，北京市冬季約30%的天氣為“好”天氣，23%為“壞”天氣。從1969—1989年的20年時間里，“特壞”天氣從無到有，比重保持在5%以下，“壞”天氣有所增加，而“好”天氣在逐步減少，“極好”天氣的頻率基本保持穩定?？梢姡本┦械沫h境狀況從80年代起已經開始受到一定的污染。

圖5 北京市歷年冬季不同可視范圍下能見度頻率柱狀圖

4 結論

本文基于M-GGPFR模型構建能見度與氣象因子的關聯并進行估計和預測分析。模型可以刻畫氣象因子與能見度的非線性相關性，而且可以借助不同站點數據信息，對目標站點進行多類型預測。數值分析說明模型很大程度改進了GGPFR模型的估計和預測效果，并且優于普通的回歸模型和時間序列模型。在應用上，選取了部分城市能見度歷史數據建模，并重點分析了北京市的能見度年際演變規律。由于能見度數據前后記錄格式不同，常規方法很難獲得年際演變規律。本文模型可以將等級數據背后的潛在數值進行恢復，進而在同一尺度下對數據進行比較分析。從結構上看，該模型可以視為廣義線性模型與函數型數據模型的結合。除等級數據外，本文模型也適用于其他非高斯型函數數據，并且可以利用其進行聚類和判別分析。