王德冬，秦 聰

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

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

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

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

1 時空回歸克里格方法

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

在時空隨機場下，定義時空環(huán)境變量：

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

(1)

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

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

(2)

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

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

a5s12+a6s22+a7t+a8t2+a9ts1+a10ts2

(3)

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

1.2 時空克里格方法

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

(4)

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

(5)

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

(6)

式中:p0為待估時空位置;pi(i=1,…,n)為p0周圍的n個時空樣點的位置;λi為需要求取的樣點的權(quán)重系數(shù)。n一般取4～8，因為由于屏蔽效應，過多的鄰近點對插值結(jié)果影響不大，在本文中n取8。解該矩陣，得到λi的值，則p0時空位置處的估計值z*(p0)：

(7)

1.3 精度比較

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

1.4 實驗數(shù)據(jù)

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

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

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

2 結(jié)果與討論

因為使用的是十疊交叉驗證法，因此所有建模和插值被計算了10次，為使文章簡潔，僅展示第一次實驗的計算結(jié)果。

2.1 時空趨勢建模

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

2.2 時空變異模型

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

圖3 2014年蘇南地區(qū)PM2.5時空分布趨勢模型整體結(jié)果及其從西到東的時態(tài)分布趨勢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 殘差與取對數(shù)后的原始樣點的經(jīng)驗時空變異函數(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 時空插值

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

圖5 基于時空回歸克里格和時空普通克里格的PM2.5時空分布圖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 精度評價

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

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

2.5 討論

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

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

3 結(jié)語

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