殷利平 劉宵瑜 盛紹學(xué) 溫華洋 邱康俊

0 引言

為了將氣象研究對社會的積極作用融入到公共服務(wù)中,中國氣象局于2002年投資建設(shè)“三站四網(wǎng)”的大氣監(jiān)測工程,在全國各地陸續(xù)建立自動氣象站.這些自動氣象站引進許多高精度的氣象觀測設(shè)備,大大提高了對氣象要素進行實時探測的能力[1-2].安徽省大多數(shù)氣象站采用散射式能見度儀采樣,但是在日常工作中能見度儀會出現(xiàn)采樣數(shù)據(jù)缺失的情況,一般由以下幾種情況造成:1)能見度儀的鏡頭前或兩個鏡頭之間有異物堵塞,如蜘蛛結(jié)網(wǎng)、小鳥做窩等;2)在一些施工區(qū),或省道縣道等公路旁,灰沙和揚塵可能導(dǎo)致能見度采樣區(qū)內(nèi)顆粒物變化不定;3)惡劣天氣下,局部地區(qū)的風(fēng)速、風(fēng)向變化大且快,導(dǎo)致樹葉、雜物被吹起恰好位于能見度儀的采樣區(qū)內(nèi),雨雪天氣和天氣寒冷凝結(jié)的冰霜也可能使能見度儀鏡頭表面受污染嚴(yán)重,導(dǎo)致能見度數(shù)據(jù)不準(zhǔn)確;4)傳感器各接線端出現(xiàn)接觸不良、松動,以及傳感器的某一單元模塊發(fā)生故障、儀器年久失修得不到有效的維護等情況[3-4].有些自動氣象站建立在高山丘陵地帶,人工維護難度大、成本高,迫切需要一種既可以及時得到氣象站所測的完備氣象信息,又可以減輕工作人員對問題氣象站進行維護的工作量的方法.

目前,處理能見度儀數(shù)據(jù)缺失的方法主要可以分為基于統(tǒng)計的修補算法、基于鄰近性的修補算法、基于機器學(xué)習(xí)的修補算法三大類.基于統(tǒng)計的修補算法包括均值插補[5]、回歸插補[6]、多重插補[7]等,其中均值插補以數(shù)據(jù)序列的平均值作為填充缺失值;回歸插補是把缺失屬性作為因變量,其他相關(guān)屬性作為自變量,利用它們之間的關(guān)系建立回歸模型來預(yù)測缺失值的;多重插補是用一組近似值替換每個缺失值,再用標(biāo)準(zhǔn)的統(tǒng)計分析過程對多次替換后產(chǎn)生的若干數(shù)據(jù)進行分析、比較,從而得到缺失值的估計值.基于統(tǒng)計的插補方法雖然簡單易操作,但容易扭曲數(shù)據(jù)分布,且該類算法需要預(yù)先知道數(shù)據(jù)分布特征,但很多實際應(yīng)用場景中卻無法得到.基于鄰近性的修補算法中最具有代表性的是K近鄰算法(K-Nearest Neighbor,KNN)[8-9].K近鄰算法首先要找出數(shù)據(jù)集中與缺失數(shù)據(jù)的歐式距離最小的K個點,然后用這K個點的平均值替換缺失值,其修補效果易受到鄰近閾值的影響,且容易受到噪聲數(shù)據(jù)的干擾,若對數(shù)據(jù)集未做初步預(yù)處理,修補精度容易產(chǎn)生較大的偏差.基于機器學(xué)習(xí)的修補算法能夠直接處理缺失數(shù)據(jù),并對缺失數(shù)據(jù)集進行訓(xùn)練,該類方法的優(yōu)點是可以直接處理完全隨機缺失模式下的數(shù)據(jù)集.該類算法主要包括:集成方法(以神經(jīng)網(wǎng)絡(luò)集成方法為主)[10]、多層感知機插補[11]、決策樹、貝葉斯[12-13]、支持向量機(Support Vector Machine,SVM)[14-15]等,其中集成方法修補缺失數(shù)據(jù)以BP神經(jīng)網(wǎng)絡(luò)應(yīng)用最為廣泛.BP神經(jīng)網(wǎng)絡(luò)是指利用誤差逆?zhèn)鞑?Error Back Propagation)算法訓(xùn)練的多層網(wǎng)絡(luò),BP算法是將誤差反向傳播使神經(jīng)元各層權(quán)值不斷調(diào)整,直到網(wǎng)絡(luò)輸出的誤差減少到可接受的程度,其優(yōu)點是在處理不完整大規(guī)模數(shù)據(jù)時速度快、泛化能力強.SVM也是一種通用的機器學(xué)習(xí)算法,它以統(tǒng)計學(xué)習(xí)理論為基礎(chǔ),廣泛應(yīng)用在函數(shù)回歸、時間序列預(yù)測等領(lǐng)域[16-17].SVM算法首先是通過非線性映射函數(shù)把樣本向量映射到高維特征空間,使得在特征空間中,原空間數(shù)據(jù)的像具有線性關(guān)系,然后在特征空間中構(gòu)造線性最優(yōu)決策函數(shù),從而解決分類與回歸問題.在處理缺失值回歸分析時,SVM算法可以修補任意缺失模式的數(shù)據(jù),減少計算復(fù)雜度[18].

本文利用機器學(xué)習(xí)的相關(guān)算法在缺失值插補方面的優(yōu)異性,綜合運用SVM和BP神經(jīng)網(wǎng)絡(luò)預(yù)估能見度缺失值.首先選用安徽氣象局歷年來不同地區(qū)氣象站的歷史數(shù)據(jù)進行分析,然后建立數(shù)據(jù)填充模型,再運用權(quán)重優(yōu)化不同模型對缺失值的預(yù)估值.實驗結(jié)果表明,運用組合模型對不同地形的能見度缺失數(shù)據(jù)進行預(yù)估,預(yù)估結(jié)果可以有效地代替真實值,實現(xiàn)了對自動氣象站的缺失數(shù)據(jù)的高精度填補.

1 數(shù)據(jù)的來源和模型構(gòu)建

1.1 數(shù)據(jù)處理

本文中氣象數(shù)據(jù)全部來自安徽氣象局歷年來匯總的氣象站觀測資料.如圖1所示,安徽地形復(fù)雜多樣,不同地形氣候不一,因此所得到的觀測數(shù)據(jù)差值較大.考慮到地形因素對模型處理缺測數(shù)據(jù)的影響,本文以高山、山谷、平原、水源地地形作為特征,分別選取黃山站(高山地形)、山南溪谷站(山谷地形)、靈璧站(平原地形)和白澤湖站(水源地地形)的歷史數(shù)據(jù)進行試驗(圖2).早期的自動氣象站由于設(shè)備質(zhì)量參差不齊,傳感器檢測精度低,得到的數(shù)據(jù)不完整.為了保證數(shù)據(jù)的有效性,本文選取從2015—2019年安徽省氣象局記錄較為完整的小時時序數(shù)據(jù)資料作為總樣本.對于每種地形,將相應(yīng)的樣本分為10個樣本集,其中7個樣本作為訓(xùn)練集,3個樣本集作為驗證集.測試集選取各站點2019年春季3—4月小時觀測資料,一些時間區(qū)間內(nèi)如果能見度數(shù)值變化較小,選取的數(shù)據(jù)量也相應(yīng)減少.

圖1 安徽省地形圖Fig.1 Topographic map of Anhui Province

圖2 安徽省各區(qū)域小型氣象站分布Fig.2 Distribution map of meteorological stations in Anhui Province

1.2 輸入要素的選擇

氣象觀測要素很多,有些氣象觀測要素對能見度的影響很小,如果將一個時序的全部觀測要素作為輸入不僅計算量大,而且會影響預(yù)估結(jié)果的準(zhǔn)確性.人為篩選輸入要素的方法具有很大的主觀性,缺乏理論依據(jù),因此本文選用灰色關(guān)聯(lián)分析法進行輸入變量的選擇[19-20].該方法是根據(jù)各因素之間數(shù)值變化趨勢的程度來確定關(guān)聯(lián)大小,這種方法對數(shù)據(jù)要求較低,步驟清晰且計算量小.灰色關(guān)聯(lián)法中一個重要指標(biāo)是灰色關(guān)聯(lián)度,灰色關(guān)聯(lián)度以數(shù)值的形式表征各變量間關(guān)系的強弱.本文對氣象各要素之間的灰色關(guān)聯(lián)分析步驟如下:

1)為研究能見度要素與其他氣象要素之間的關(guān)系,先對氣象要素進行編號,記為Ai(i=1,2,3,4,…,11),并將各氣象要素數(shù)據(jù)換算成標(biāo)準(zhǔn)單位制數(shù)值,Ai與各氣象要素之間的對應(yīng)關(guān)系如表1所示.

表1 各氣象觀測要素的編號序列

2)求表征關(guān)聯(lián)度的關(guān)聯(lián)系數(shù).以能見度數(shù)據(jù)序列為參考數(shù)據(jù),其他觀測要素的數(shù)據(jù)序列作為比較數(shù)據(jù).參考數(shù)列為A0={A0(1),…,A0(d),…,A0(N)},比較數(shù)列為Ai={Ai(1),…,Ai(d),…,Ai(N)},其中d代表各要素序列中的元素個數(shù),N是選取數(shù)據(jù)序列的總數(shù),1≤d≤N.

數(shù)據(jù)序列Ai與A0在第d點的關(guān)聯(lián)系數(shù)εi(d)為

εi(d)=

(1)

式中:ρ∈(0,+∞) 稱為分辨系數(shù),通常在0到1之間選取,一般取ρ=0.5;i代表氣象要素序號,1≤i≤11.

3)求各氣象要素之間的關(guān)聯(lián)度γ(A0,Ai):

(2)

其中關(guān)聯(lián)度γ∈(0,1),數(shù)值越大表明該氣象要素與能見度的關(guān)聯(lián)度越高.本文按照關(guān)聯(lián)度數(shù)值從大到小的順序選擇輸入要素,不妨設(shè)選擇的輸入要素為M個.

1.3 能見度數(shù)據(jù)預(yù)估模型

1.3.1 SVM能見度缺失值預(yù)估模型

SVM是把線性不可分的樣本通過核函數(shù)映射到特征空間,進而在特征空間中構(gòu)造最優(yōu)分類平面,使樣本到平面的總距離最小,由此實現(xiàn)擬合的[21].對于模型給定的訓(xùn)練數(shù)據(jù)總樣本D={(xi(j),y(j)),i=1,2,…,M,j=1,2,…,N1},其中xi(j)為第i個氣象輸入要素的第j個樣本,y(j)為對應(yīng)的能見度實測值,N1為總樣本容量.記x=[x1,x2,…,xM]T.首先利用一個非線性映射函數(shù)φ(x)將樣本x從原空間RM映射到特征空間,然后在高維特征空間中構(gòu)造最優(yōu)決策函數(shù):

y(x)=wT·φ(x)+b,

(3)

式中:φ(x)為映射函數(shù);w為權(quán)向量;b為偏置量.權(quán)向量w與b通過優(yōu)化下式得到:

(4)

其約束條件為

y(j)-wTφ(x(j))+b+ξ(j)=0,

(5)

式(4)中:C為懲罰因子,為給定值,其數(shù)值越大表示對訓(xùn)練誤差大于設(shè)定誤差的樣本懲罰越大;ξ(j)為松弛變量,定義為ξ(j)=1-y(j),ξ(j)數(shù)值越大表示對樣本訓(xùn)練誤差的容忍程度越大.

在求解最小化問題(4)和(5)之前,首先要找到合適的非線性函數(shù)φ(x),為此引入徑向基核函數(shù):

j=1,2,…,N1,

(6)

并令K(x(β),x(j))=φT(x(β))·φ(x(j)),β=1,2,…,N1.進一步引入Lagrange方程,從而可以求解(4)和(5),得出SVM最優(yōu)決策函數(shù)的估計函數(shù)為

(7)

式中K(x,x(j))為核函數(shù),拉格朗日乘子α(j)∈R,R為實數(shù)集.

本文根據(jù)以上原理,通過能見度與其他氣象觀測要素之間復(fù)雜的非線性關(guān)系進行能見度缺失值預(yù)估.具體步驟如下:

2)運用網(wǎng)絡(luò)搜索法來分別對式(4)和式(6)中的C,δ兩個參數(shù)尋優(yōu),其中懲罰因子C的搜尋范圍在0.1～100,核參數(shù)δ的搜索范圍在0.001～1,利用交叉驗證法可獲得最優(yōu)參數(shù)[22].

3)利用建立的SVM能見度數(shù)據(jù)預(yù)估模型,對預(yù)處理后的樣本數(shù)據(jù)進行訓(xùn)練,并對模型的預(yù)估結(jié)果進行評價.選用平均相對誤差(MAPE,其量值記為ηMAPE)和均方根誤差(RMSE,其量值記為ηRMSE)來評價:

(8)

(9)

1.3.2 BP神經(jīng)網(wǎng)絡(luò)的能見度缺失值預(yù)估模型

BP神經(jīng)網(wǎng)絡(luò)能見度預(yù)估的基本結(jié)構(gòu)如圖3所示,其中輸入層有m個節(jié)點,隱含層有p個節(jié)點,輸出層有1個節(jié)點,Wig(i=1,2,…,m;g=1,2,…,p)為輸入層到隱含層的權(quán)值,Wgk(g=1,2,…,p;k=1)為隱含層到輸出層的權(quán)值,θg(g=1,2,…,p)為隱含層的閾值,σ1為輸出層閾值,(X1,X2,…,Xm)為神經(jīng)網(wǎng)絡(luò)輸入向量,Y1為神經(jīng)網(wǎng)絡(luò)輸出量,Yh為期望輸出,e為神經(jīng)網(wǎng)絡(luò)期望輸出與實際輸出的誤差.

圖3 BP神經(jīng)網(wǎng)絡(luò)能見度預(yù)估的結(jié)構(gòu)Fig.3 Structure of BP neural network for visibility estimation

三層結(jié)構(gòu)BP神經(jīng)網(wǎng)絡(luò)可用于預(yù)估氣象站能見度缺失值,其中輸入層對應(yīng)與能見度相關(guān)性大的氣象要素序列,輸出層是能見度預(yù)估值.隱含層的神經(jīng)元數(shù)量對模型預(yù)估結(jié)果的好壞產(chǎn)生直接的影響,但是目前沒有能直接確認(rèn)最優(yōu)隱含層個數(shù)的方法,只有根據(jù)以下經(jīng)驗公式來計算：

(10)

設(shè)定不同隱含層神經(jīng)元個數(shù),然后采用“試湊法”逐步增大和減少隱含層神經(jīng)元數(shù)目使網(wǎng)絡(luò)誤差最小.式(10)中:m為輸入層節(jié)點個數(shù);n為輸出層節(jié)點個數(shù);a為常數(shù),取值范圍一般為3至10.K為隱含層神經(jīng)元估算個數(shù).

用于能見度缺失值預(yù)估的BP算法各步驟如下:

1)權(quán)值初始化:(wig∪wgk)=random(·),其中random(·)表示權(quán)值在[0,1]之間的均勻分布.

2) 依次輸入訓(xùn)練集中的樣本,設(shè)當(dāng)前輸入第q個樣本.

3) 依次計算各層的輸出:X′g,X″k及Y1.

4) 求各層的反傳誤差,并記下各個X″k(q),X′g(q),Xi(q)的值.

5)記錄已學(xué)習(xí)過的樣本個數(shù)q.如果q

6)按新的權(quán)值再計算X′g,X″k及Y1和學(xué)習(xí)樣本數(shù)的總誤差E,若E<ε(ε為預(yù)估給定誤差),或達到最大學(xué)習(xí)次數(shù),則終止學(xué)習(xí).否則,轉(zhuǎn)步驟2)繼續(xù)新一輪學(xué)習(xí).

1.3.3 組合模型

用不同的機器學(xué)習(xí)算法得出的能見度預(yù)估值與實測值都有誤差.為了減小預(yù)估值與實測值之間的誤差,可以整合不同模型的優(yōu)點,對不同方法的預(yù)估結(jié)果進行加權(quán)組合,以提高預(yù)估精度.在組合模型預(yù)估中最關(guān)鍵的步驟是確定不同預(yù)估方法的權(quán)重.目前,針對多模型組合權(quán)重確定,常用的方法主要有以下幾種:算術(shù)平均法、方差倒數(shù)法、均方倒數(shù)法以及最小二乘法.本文采用方差倒數(shù)法判斷單項模型系數(shù),即對誤差平方和小的模型賦予較高的權(quán)重,誤差平方和大的賦予較小的權(quán)重,使組合模型的誤差和盡可能小.具體方法如下:

設(shè)F為觀測對象,其實際觀測值向量為(F1,F2,…,Fn),U1,U2,…,Ur為r種不同預(yù)估方法得出的預(yù)估值,向量S=(S1,S2,…,Sr)T中元素分別是它們在組合模型中的權(quán)重,第l個預(yù)估方法Ul的預(yù)估值為(U1l,U2l,…,Unl).則組合模型的估計值為

(11)

其中

(12)

2 實驗與分析

2.1 能見度影響因子的選擇

一般關(guān)聯(lián)度大于等于0.8時,子序列與母序列關(guān)聯(lián)度很好.根據(jù)1.2節(jié)的理論,可以計算得到其他觀測要素與能見度之間的關(guān)聯(lián)度,本次實驗選擇與能見度的關(guān)聯(lián)度在0.8以上的氣象觀測要素,如表2所示.

表2 部分氣象要素與能見度之間的關(guān)聯(lián)度

2.2 SVM與BP方法組合模型預(yù)估

本文采用的SVM模型預(yù)估能見度實驗,借助的是Pycharm軟件的libsvm工具箱,其實驗精度主要取決于參數(shù)選取是否合適,本文各參數(shù)設(shè)定值如表3所示.

表3 SVM最優(yōu)參數(shù)值設(shè)置

在BP預(yù)估實驗中,BP神經(jīng)網(wǎng)絡(luò)模型以與能見度要素關(guān)聯(lián)度高的8個氣象要素作為輸入,隱含層選用單層結(jié)構(gòu),依據(jù)估算最優(yōu)隱含層神經(jīng)元個數(shù)的經(jīng)驗公式,推算出K值在[6,13]之間.為了保證隱含層神經(jīng)元個數(shù)對模型預(yù)估結(jié)果的準(zhǔn)確性,設(shè)定隱含層神經(jīng)元個數(shù)K在[5,20]區(qū)間.將平原組訓(xùn)練集歸一化處理后輸入BP模型中,取不同隱含層神經(jīng)元個數(shù),用一組驗證集記錄相對誤差均值.由于初始權(quán)值隨機分配,相同個數(shù)的隱含層神經(jīng)元運行的結(jié)果也有不同,所以BP網(wǎng)絡(luò)中每個K值的設(shè)定都運行10次,誤差結(jié)果算平均值,尋優(yōu)過程如圖4所示.

圖4 隱含層神經(jīng)元個數(shù)尋優(yōu)Fig.4 Optimization of the number of neurons in hidden layer

由隱含層神經(jīng)元尋優(yōu)結(jié)果可知,BP神經(jīng)網(wǎng)絡(luò)預(yù)估模型選用10個隱含層神經(jīng)元最佳.其中BP神經(jīng)網(wǎng)絡(luò)的訓(xùn)練最大迭代次數(shù)設(shè)定為50 000次,學(xué)習(xí)率取0.1,迭代循環(huán)次數(shù)上限值取20,訓(xùn)練最終誤差設(shè)定為0.001,激活函數(shù)選擇雙曲正切函數(shù),訓(xùn)練函數(shù)及學(xué)習(xí)函數(shù)均采用Levenberg-Marquardt算法.

對不同地形代表站的測試集數(shù)據(jù)分別采用已訓(xùn)練好的SVM和BP模型進行預(yù)估,測試集預(yù)估結(jié)果和實測結(jié)果進行對比分析,求出各地形預(yù)估值與實測值的誤差平方和,運用方差倒數(shù)法,得各自的權(quán)重系數(shù)如表4所示.

表4 兩種方法單一預(yù)估結(jié)果

2.3 組合模型預(yù)估實驗結(jié)果

將兩種方法的訓(xùn)練集預(yù)估結(jié)果進行對比分析,由表4可知,在這四種地形中,SVM缺失值預(yù)估模型要比BP神經(jīng)網(wǎng)絡(luò)的更加穩(wěn)定,誤差也更小,但是在水源地和平原地形中BP神經(jīng)網(wǎng)絡(luò)的預(yù)估結(jié)果準(zhǔn)確度相對更高,結(jié)合兩種模型預(yù)估的結(jié)果,可以提高能見度預(yù)估的精度.實驗輸入測試集數(shù)據(jù)得到兩種模型的預(yù)估結(jié)果,運用上文所述的方差倒數(shù)法,加權(quán)組合求出組合模型的預(yù)估值,并計算組合模型預(yù)估值和實測值的平均相對誤差、誤差均值和最大誤差.測試集能見度組合模型預(yù)估性能指標(biāo)結(jié)果如表5所示,預(yù)估效果如圖5所示.

表5 組合方法的性能參數(shù)對比

圖5 各地形能見度數(shù)據(jù)組預(yù)估效果Fig.5 Comparison of observed visibility and estimation by the combined SVM-BP neural network method for mountainous (a),plain (b),valley (c),and water source (d) areas

從表5中的實驗結(jié)果數(shù)據(jù)可以看出,無論是哪種地形,組合模型預(yù)估的平均相對誤差更低,整體誤差均值小,效果要明顯好于單一模型.由此可知,組合方法可以保留單一模型的預(yù)估優(yōu)勢,增加對缺失數(shù)據(jù)預(yù)估的可靠性.

3 結(jié)論

為解決自動氣象站能見度要素缺測的問題,本文利用組合模型對缺測數(shù)據(jù)進行精確的預(yù)估,并以預(yù)估值代替實測值來保證數(shù)據(jù)的完備性.首先通過灰色關(guān)聯(lián)分析方法精簡預(yù)估模型輸入,篩選出與能見度相關(guān)度較大的氣象要素,再從氣象信息的多種要素中建立能見度數(shù)據(jù)預(yù)估模型.在實驗部分,兩種模型對能見度都具有良好的預(yù)估能力,SVM模型在四種地形中對能見度數(shù)據(jù)的預(yù)估結(jié)果比較好且穩(wěn)定,而BP神經(jīng)網(wǎng)絡(luò)則對平原和水源地的數(shù)據(jù)預(yù)估能力突出.對不同模型預(yù)估的結(jié)果加權(quán)組合,結(jié)合測試集的RMSE和MAPE的數(shù)據(jù),將這些數(shù)據(jù)與單一預(yù)估的方法對比,結(jié)果表明組合模型預(yù)估的方法更加接近實測值,更能充分地利用數(shù)據(jù)信息,從而提高對缺失數(shù)據(jù)預(yù)估的準(zhǔn)確性.