王德冬，秦 聰

1. 山東省遙感技術(shù)應(yīng)用中心，山東 濟(jì)南 250013 2. 山東光庭信息技術(shù)有限公司，山東 煙臺(tái) 264000

近年來(lái)，大氣環(huán)境質(zhì)量已成為社會(huì)關(guān)注的熱點(diǎn)問(wèn)題之一。從2013年開(kāi)始，我國(guó)政府陸續(xù)向大眾發(fā)布從各城市監(jiān)測(cè)站點(diǎn)測(cè)得的實(shí)時(shí)大氣顆粒物濃度數(shù)據(jù)，到目前為止，全國(guó)已布設(shè) 1 600個(gè)監(jiān)測(cè)站點(diǎn)，越來(lái)越多的研究開(kāi)始關(guān)注區(qū)域PM2.5的污染狀況及其時(shí)空分異特征，如王振波等[1]基于全國(guó)190個(gè)城市中945個(gè)監(jiān)測(cè)站的觀測(cè)數(shù)據(jù)，采用空間數(shù)據(jù)統(tǒng)計(jì)模型揭示了2014年中國(guó)城市PM2.5濃度時(shí)空變化規(guī)律，發(fā)現(xiàn)PM2.5濃度在時(shí)間上呈現(xiàn)顯著的冬秋高、春夏低“U”型變化規(guī)律與顯著的空間分異與集聚規(guī)律；楊冕等[2]運(yùn)用地理學(xué)時(shí)空分析與GIS可視化方法探索并呈現(xiàn)了2015年長(zhǎng)江經(jīng)濟(jì)帶PM2.5的時(shí)空分布特征及其演變規(guī)律；徐偉嘉等[3]利用空間地統(tǒng)計(jì)方法，分析了珠三角區(qū)域PM2.5時(shí)空分異特征；徐建輝等[4]利用MODIS遙感數(shù)據(jù)對(duì)長(zhǎng)三角PM2.5濃度進(jìn)行了建模和估算，并分析了其時(shí)空分布特征；劉永林等[5]利用17個(gè)空氣質(zhì)量監(jiān)測(cè)站數(shù)據(jù)分析了重慶市主城區(qū)PM2.5時(shí)空分布特征。

大部分研究均是根據(jù)原始監(jiān)測(cè)點(diǎn)數(shù)據(jù)對(duì)區(qū)域PM2.5進(jìn)行時(shí)空分異特征研究，當(dāng)需要獲得某段時(shí)期內(nèi)(如月份、季節(jié)和年)空間連續(xù)型分布數(shù)據(jù)時(shí)，往往采用空間插值方法對(duì)樣點(diǎn)均值進(jìn)行空間插值得到。這種處理方法忽略了樣點(diǎn)本身存在的時(shí)空相關(guān)性和變異性，丟失了大量時(shí)空變化信息，而且得到的只是某時(shí)期內(nèi)均值的空間分布信息，難以進(jìn)行更為深入的時(shí)空分析。針對(duì)這一點(diǎn)，梅楊等[6]利用監(jiān)測(cè)數(shù)據(jù)對(duì)山東省2014年P(guān)M2.5濃度進(jìn)行了時(shí)空建模及插值，得到了區(qū)域PM2.5時(shí)空立方體數(shù)據(jù)，進(jìn)而進(jìn)行了時(shí)空分析，總結(jié)了山東省PM2.5時(shí)空分布特征。此研究為獲取區(qū)域PM2.5時(shí)空立方體數(shù)據(jù)提供了一條途徑。同時(shí)，注意到雖然時(shí)空地統(tǒng)計(jì)方法插值結(jié)果會(huì)高于只使用同一時(shí)期樣點(diǎn)的空間地統(tǒng)計(jì)方法，但由于一般區(qū)域的監(jiān)測(cè)點(diǎn)數(shù)量相對(duì)于區(qū)域面積來(lái)講，較為稀疏，因此時(shí)空估值精度仍然較低，如上述研究中，時(shí)空克里格估值的均方根誤差為22.08 μg/m3，相對(duì)于區(qū)域時(shí)空樣點(diǎn)均值74.84 μg/m3仍然較高。因此，如何提高區(qū)域PM2.5時(shí)空估值精度，進(jìn)而提高時(shí)空分析的可靠性，十分值得研究。

同時(shí)，注意到上述關(guān)于區(qū)域PM2.5的研究，均揭示出由于地理位置特征和季節(jié)的氣候差異，PM2.5存在顯著的時(shí)空分布趨勢(shì)，而對(duì)其他環(huán)境屬性的地統(tǒng)計(jì)建模和分析研究中，往往考慮了空間或時(shí)空趨勢(shì)的估值結(jié)果要高于忽略了這些趨勢(shì)的結(jié)果[7-9]。基于上述考慮，本文擬利用時(shí)空回歸克里格對(duì)區(qū)域PM2.5進(jìn)行時(shí)空建模及插值，目的基于監(jiān)測(cè)數(shù)據(jù)提出PM2.5時(shí)空分布趨勢(shì)模型；對(duì)剝離趨勢(shì)后的時(shí)空樣點(diǎn)殘差建立時(shí)空變異模型，并利用時(shí)空普通克里格方法進(jìn)行插值；將時(shí)空分布趨勢(shì)與殘差的克里格插值結(jié)果，作為最后的預(yù)測(cè)結(jié)果；與不考慮時(shí)空趨勢(shì)模型的時(shí)空普通克里格插值進(jìn)行精度對(duì)比，分析所提方法優(yōu)勢(shì)與不足。

1 時(shí)空回歸克里格方法

1.1 時(shí)空環(huán)境變量定義

在時(shí)空隨機(jī)場(chǎng)下，定義時(shí)空環(huán)境變量：

z(p):p=(s,t),s=(s1,s2)∈S,t∈T

(1)

式中：p表示時(shí)空位置；z(p)表示p位置處時(shí)空變量；s和t分別表示空間坐標(biāo)和時(shí)間；s1和s2分別表示空間坐標(biāo)中2個(gè)橫、縱坐標(biāo)分量；s和T分別表示空間域和時(shí)間域。當(dāng)s和t確定時(shí)，就認(rèn)為取得了一次時(shí)空變量的實(shí)現(xiàn)(即觀測(cè)值)。同時(shí)，為了研究環(huán)境變量的時(shí)空趨勢(shì)及其對(duì)時(shí)空預(yù)測(cè)及制圖的影響，環(huán)境變量被定義為趨勢(shì)項(xiàng)和殘差項(xiàng)之和：

z(p)=m(p)+R(p)

(2)

式中m(p)和R(p)分別表示趨勢(shì)項(xiàng)和殘差項(xiàng)。在沒(méi)有任何輔助數(shù)據(jù)的支持下，根據(jù)KYRIAKIDIS等[10]的研究，關(guān)于時(shí)空位置的2次多項(xiàng)式已經(jīng)足夠表達(dá)地理屬性的時(shí)空分布趨勢(shì)。因此本研究中，將PM2.5的時(shí)空分布趨勢(shì)表示為：

m(p)=m(s1,s2,t)=a1+a2s1+a3s2+a4s1s2+

a5s12+a6s22+a7t+a8t2+a9ts1+a10ts2

(3)

式中的系數(shù)ai將根據(jù)建模樣點(diǎn)、利用matlab的趨勢(shì)分析技術(shù)擬合。

1.2 時(shí)空克里格方法

對(duì)所有參與建模的樣點(diǎn)，減去樣點(diǎn)位置處的時(shí)空趨勢(shì)后，得到樣點(diǎn)殘差R(p)，基于此，計(jì)算各時(shí)空滯后距的經(jīng)驗(yàn)時(shí)空半方差值：

(4)

式中：hs和ht表示空間和時(shí)間上的滯后距；N(hs,ht)表示時(shí)空樣點(diǎn)中符合上述滯后距的點(diǎn)對(duì)數(shù)量；γ(hs,ht)表示在時(shí)空滯后距hs和ht上的經(jīng)驗(yàn)時(shí)空半方差值。計(jì)算若干組不同時(shí)空滯后距的時(shí)空實(shí)驗(yàn)變異函數(shù)后，形成時(shí)空實(shí)驗(yàn)變異函數(shù)散點(diǎn)圖。為了獲取任意時(shí)空滯后距上的時(shí)空變異函數(shù)值，需利用該散點(diǎn)圖，擬合理論時(shí)空變異函數(shù)模型，本文選取的模型形式：

(5)

上述模型中各參數(shù)(c0,c,v,w,ξ)擬合后，聯(lián)合時(shí)空樣點(diǎn)，基于時(shí)空克里格法對(duì)待估時(shí)空位置的PM2.5殘差進(jìn)行估值，克里格矩陣：

(6)

式中:p0為待估時(shí)空位置;pi(i=1,…,n)為p0周?chē)膎個(gè)時(shí)空樣點(diǎn)的位置;λi為需要求取的樣點(diǎn)的權(quán)重系數(shù)。n一般取4～8，因?yàn)橛捎谄帘涡?yīng)，過(guò)多的鄰近點(diǎn)對(duì)插值結(jié)果影響不大，在本文中n取8。解該矩陣，得到λi的值，則p0時(shí)空位置處的估計(jì)值z(mì)*(p0)：

(7)

1.3 精度比較

為比較趨勢(shì)分析和回歸克里格的插值精度，同時(shí)使用原始樣點(diǎn)數(shù)據(jù)按照式(4)的方式計(jì)算實(shí)驗(yàn)變異函數(shù)，并擬合式(5)的理論模型參數(shù)，最后用式(6)對(duì)各待估時(shí)空位置進(jìn)行時(shí)空估值。對(duì)實(shí)驗(yàn)數(shù)據(jù)，使用十疊交叉驗(yàn)證法對(duì)比2種方法所得精度，即每次預(yù)留10%的監(jiān)測(cè)站(即5個(gè))，將其所有監(jiān)測(cè)數(shù)據(jù)作為驗(yàn)證數(shù)據(jù)集，其余作為建模數(shù)據(jù)集，對(duì)比2種方法所得的均方根誤差，以此來(lái)評(píng)價(jià)2種方法所得精度。

1.4 實(shí)驗(yàn)數(shù)據(jù)

以蘇南地區(qū)2014年P(guān)M2.5日監(jiān)測(cè)數(shù)據(jù)為實(shí)驗(yàn)數(shù)據(jù)，該地區(qū)共有大氣監(jiān)測(cè)站53個(gè)，監(jiān)測(cè)站空間分布如圖1所示，其統(tǒng)計(jì)特征如圖2所示。2014年，蘇南地區(qū)PM2.5的最大值和最小值分別為5 μg/m3和531 μg/m3，均值為65.6 μg/m3。

圖1 蘇南地區(qū)大氣監(jiān)測(cè)站空間分布圖Fig.1 Spatial locations of the monitoring sites in southern Jiangsu during 2014

圖2 蘇南地區(qū)大氣監(jiān)測(cè)數(shù)據(jù)統(tǒng)計(jì)特征Fig.2 The summary statistics and histogram of the PM2.5 concentrations of all the collected data in southern Jiangsu during 2014

2 結(jié)果與討論

因?yàn)槭褂玫氖鞘B交叉驗(yàn)證法，因此所有建模和插值被計(jì)算了10次，為使文章簡(jiǎn)潔，僅展示第一次實(shí)驗(yàn)的計(jì)算結(jié)果。

2.1 時(shí)空趨勢(shì)建模

依據(jù)式(3)，建立研究區(qū)PM2.5與時(shí)空位置間的趨勢(shì)模型m(s1,s2,t)=49.56-1.879×10-4s2+5.448×10-5s1+4.298×10-11s1s2-3.788×10-11s22-3.879×10-12s12-1.487×10-5t+9×10-4t2+3.468×10-7s2-2.499×10-7ts1。

圖3展示了該趨勢(shì)模型的整體結(jié)果及其從西到東的三維效果。從圖3可以看出，時(shí)間上該地區(qū)PM2.5呈現(xiàn)明顯冬、春季高，夏、秋季低的“U”型變化規(guī)律，空間上從東到西呈遞增趨勢(shì)。這是因?yàn)?，蘇南地區(qū)東部靠海，受到風(fēng)力和溫度的影響，東部地區(qū)PM2.5濃度相對(duì)西部較低。同時(shí)，用建模點(diǎn)趨勢(shì)值與實(shí)際值的平均絕對(duì)誤差(MAE)、平均誤差(ME)和相關(guān)系數(shù)(r)來(lái)量化建模效果。其中ME為0.003，說(shuō)明建模點(diǎn)實(shí)際值均勻地分布在所得趨勢(shì)(體)的上下方；MAE為18.76，相對(duì)于該地區(qū)PM2.5樣點(diǎn)均值65.63還是較高，說(shuō)明有的監(jiān)測(cè)點(diǎn)實(shí)測(cè)值距離所得趨勢(shì)(體)較遠(yuǎn)，誤差較大。

2.2 時(shí)空變異模型

對(duì)建模組的樣點(diǎn)殘差，計(jì)算K-S的值為2.79，表明樣點(diǎn)殘差接近正態(tài)分布，因此可用殘差直接進(jìn)行變異函數(shù)的計(jì)算和后續(xù)的時(shí)空插值。而原始樣點(diǎn)經(jīng)檢驗(yàn)不符合正態(tài)分布，但對(duì)數(shù)轉(zhuǎn)換后的數(shù)據(jù)接近正態(tài)分布。按照式(4)計(jì)算各時(shí)空滯后距上殘差和lgPM2.5的變異函數(shù)值，并擬合式(5)的理論模型，其結(jié)果如圖4所示。從圖4可看出，代表經(jīng)驗(yàn)變異值的黑色散點(diǎn)緊緊圍繞在代表理論時(shí)空變異模型的曲面周?chē)?，說(shuō)明理論模型擬合較好，能夠表達(dá)區(qū)域PM2.5殘差和原始數(shù)據(jù)的時(shí)空變異特征。另外，從圖4可看出，時(shí)間上，殘差和 lgPM2.5的變異值在0～4 d呈遞增趨勢(shì)，但4 d以后變異值區(qū)域平穩(wěn)，說(shuō)明時(shí)間上有效相關(guān)范圍為4 d左右。而空間上，這種遞增趨勢(shì)持續(xù)到150 km，說(shuō)明空間上的有效相關(guān)范圍至少有150 km。

圖3 2014年蘇南地區(qū)PM2.5時(shí)空分布趨勢(shì)模型整體結(jié)果及其從西到東的時(shí)態(tài)分布趨勢(shì)Fig.3 (a) ST trend model of PM2.5 concentrations in southern Jiangsu province, China, during 2014, and (b) the temporal trend for PM2.5 concentrations during 2014 in southern Jiangsu province, China.

圖4 殘差與取對(duì)數(shù)后的原始樣點(diǎn)的經(jīng)驗(yàn)時(shí)空變異函數(shù)值(散點(diǎn))與擬合的理論時(shí)空變異函數(shù)模型(曲面)Fig.4 Empirical variogram values (black dots) of the PM2.5 residuals (a) and logPM2.5 (b), and the theoretical variogram model (surface) fitted to these empirical values.

2.3 時(shí)空插值

本研究中設(shè)置得待插值時(shí)空網(wǎng)格為2 km × 2 km × 1 d，基于上述的理論時(shí)空變異模型，按照式(6)和式(7)對(duì)每個(gè)時(shí)空網(wǎng)格的殘差和lgPM2.5進(jìn)行估值。其中l(wèi)ogPM2.5的估值結(jié)果進(jìn)行冪函數(shù)轉(zhuǎn)換后，得到不考慮時(shí)空趨勢(shì)的PM2.5時(shí)空分布立方體數(shù)據(jù)，如圖5(b)所示。而對(duì)于殘差的估值結(jié)果，加上每個(gè)時(shí)空網(wǎng)格出由式(8)計(jì)算而來(lái)趨勢(shì)值，得到考慮時(shí)空趨勢(shì)的時(shí)空回歸克里格估值結(jié)果，如圖5(a)所示。

圖5 基于時(shí)空回歸克里格和時(shí)空普通克里格的PM2.5時(shí)空分布圖Fig.5 Plot of the ST map of predicted PM2.5 concentration generated by the spatiotemporal kriging with trend and spatiotemporal ordinary kriging

2.4 精度評(píng)價(jià)

使用十疊交叉驗(yàn)證進(jìn)行了精度驗(yàn)證，即每次預(yù)留5個(gè)監(jiān)測(cè)站的觀測(cè)數(shù)據(jù)作為驗(yàn)證數(shù)據(jù)集，其均方根誤差的結(jié)果如圖6所示。其中時(shí)空普通克里格方法(STOK)的平均均方根誤差為14.27，而時(shí)空回歸克里格(STRK)的平均均方根誤差為12.43。因此，STRK的插值精度優(yōu)于STOK的插值精度，精度平均提高12.9%。

圖6 STKT和STOK的十疊交叉驗(yàn)證結(jié)果Fig.6 The results of 10-fold cross validation of STKT and STOK

2.5 討論

本研究的創(chuàng)新之處在于構(gòu)建了區(qū)域PM2.5的時(shí)空趨勢(shì)模型，它很好地表達(dá)了區(qū)域PM2.5的時(shí)空分布趨勢(shì)特征，并且有效地消除了PM2.5的時(shí)空非平穩(wěn)性。而平穩(wěn)性假設(shè)或二階平穩(wěn)性是地統(tǒng)計(jì)學(xué)的基本假設(shè)，預(yù)測(cè)結(jié)果的誤差正是來(lái)源于地理屬性或環(huán)境變量對(duì)這些假設(shè)的背離，背離程度越高，精度越低。而關(guān)于PM2.5的區(qū)域性研究或插值，往往涉及的空間和時(shí)間范圍很大，這種地理上和時(shí)態(tài)上的差異使得原始數(shù)據(jù)很難符合二階平穩(wěn)性假設(shè)的2個(gè)條件，即區(qū)域化變量的均值存在且平穩(wěn)，增量的方差存在且平穩(wěn)。而趨勢(shì)模型在一定程度上表達(dá)了區(qū)域PM2.5的均值分布，將其剝離后，能夠使殘差的均值較原始數(shù)據(jù)更為平穩(wěn)(接近于0)，這可能是造成時(shí)空回歸克里格精度更高的原因之一。

當(dāng)然，本研究也存在不足之處：一是所得的趨勢(shì)模型精度仍然不高，這是因?yàn)樵斐蒔M2.5時(shí)空差異的因素眾多，其分布不但與時(shí)空位置相關(guān)，還與與社會(huì)經(jīng)濟(jì)及氣候等多種因素相關(guān)[11-13]。因此，如何綜合更多數(shù)據(jù)(如社會(huì)經(jīng)濟(jì)統(tǒng)計(jì)數(shù)據(jù)、氣候數(shù)據(jù)、地形地貌數(shù)據(jù)等)來(lái)建立更為精確的趨勢(shì)模型，將對(duì)提升時(shí)空估值精度有較大幫助；二是大氣監(jiān)測(cè)站的位置一般在城市，因此可能會(huì)造成農(nóng)村或山區(qū)估值精度降低。而且由于城市人口、交通、工業(yè)更為密集，可能會(huì)造成農(nóng)村或山區(qū)估計(jì)值大大高于其真實(shí)值。這就需要尋求與PM2.5密切相關(guān)且能全區(qū)域大范圍覆蓋的輔助要素參與建模與估值。近年來(lái)的研究發(fā)現(xiàn)，基于遙感的氣溶膠光學(xué)厚度數(shù)據(jù)與PM2.5具有顯著相關(guān)性[14-17]。因此，可以考慮將其納入PM2.5時(shí)空建模和估計(jì)中去，提升未設(shè)監(jiān)測(cè)站區(qū)域的估值精度。但目前實(shí)施起來(lái)尚有阻礙，原因是這些輔助數(shù)據(jù)在區(qū)域時(shí)空范圍內(nèi)有缺失，如氣溶膠數(shù)據(jù)為10 km×10 km的形式呈現(xiàn)，且并不是所有位置和所有時(shí)間都有數(shù)據(jù)存在，即這些數(shù)據(jù)并不能完全遍布整個(gè)區(qū)域的時(shí)空范圍。那么，如果氣溶膠參與了趨勢(shì)模型的構(gòu)建，那么在無(wú)氣溶膠觀測(cè)數(shù)據(jù)的時(shí)空位置，則無(wú)法依據(jù)該模型得到PM2.5的趨勢(shì)值。因此，如何融合這些數(shù)據(jù)進(jìn)行PM2.5時(shí)空建模與估值，需要進(jìn)一步研究。

3 結(jié)語(yǔ)

本文利用時(shí)空回歸克里格對(duì)區(qū)域PM2.5進(jìn)行了時(shí)空建模及估值，建立了PM2.5與時(shí)空位置的多元非線(xiàn)性關(guān)系，得到了區(qū)域PM2.5時(shí)空分布趨勢(shì)；計(jì)算樣點(diǎn)殘差的實(shí)驗(yàn)變異函數(shù)并擬合了理論變異函數(shù)模型；對(duì)殘差進(jìn)行時(shí)空克里格插值；將殘差的插值結(jié)果與趨勢(shì)模型結(jié)果相加，得到區(qū)域PM2.5時(shí)空回歸克里格估值結(jié)果。與不考慮趨勢(shì)的時(shí)空普通克里格結(jié)果相比，其預(yù)測(cè)精度提升了18.58%，為更精確地進(jìn)行區(qū)域PM2.5時(shí)空分析提供了更為可靠的數(shù)據(jù)基礎(chǔ)。