胡錫健, 李妍琳, 石小平
(新疆大學(xué)數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院, 烏魯木齊 830046)
在中國(guó)經(jīng)濟(jì)由高速增長(zhǎng)階段轉(zhuǎn)向高質(zhì)量發(fā)展階段過(guò)程中,環(huán)境治理和污染防治是一道重要關(guān)口。近年來(lái),推進(jìn)大氣污染防治工作卓有成效,京津冀地區(qū)空氣質(zhì)量顯著提升,長(zhǎng)三角、珠三角地區(qū)空氣質(zhì)量持續(xù)達(dá)標(biāo),汾渭平原卻成為中國(guó)空氣污染最嚴(yán)重的區(qū)域之一,引起國(guó)家和社會(huì)的高度重視[1]。2018年7月,在《打贏藍(lán)天保衛(wèi)戰(zhàn)三年行動(dòng)計(jì)劃》中,汾渭平原被納入環(huán)境污染三大重點(diǎn)防控區(qū)域之一,成為藍(lán)天保衛(wèi)戰(zhàn)的“主戰(zhàn)場(chǎng)”。
汾渭平原能源結(jié)構(gòu)以煤為主,煤炭在能源消費(fèi)中占約90%,遠(yuǎn)高于全國(guó)平均水平(60%),煤煙型污染特征明顯,大氣中的主要污染物有二氧化硫、二氧化氮、可吸入顆粒物等,這不僅對(duì)健康帶來(lái)了很大的危害,還影響著人們的正常生活。
大氣污染狀況與當(dāng)?shù)氐臍庀髼l件有著密切的關(guān)系??娒鏖诺萚2]分析了南通市2018年大氣污染物與氣象因素的關(guān)系。張?jiān)儡姷萚3]、張連華等[4]分別對(duì)汾渭平原SO2和NO2進(jìn)行時(shí)空分布特征分析,研究該區(qū)域的SO2、NO2的污染狀況,全面了解汾渭平原SO2、NO2的時(shí)空變化特征,對(duì)有效治理大氣污染具有重要意義。
氣象數(shù)據(jù)在時(shí)間尺度上有明顯的函數(shù)特征,汾渭平原11個(gè)城市的逐時(shí)氣溫?cái)?shù)據(jù)一年有9萬(wàn)多條,面對(duì)如此龐大的數(shù)據(jù)集,傳統(tǒng)的數(shù)據(jù)分析方法不能對(duì)這些數(shù)據(jù)直接處理。函數(shù)型數(shù)據(jù)分析(functional data analysis,F(xiàn)DA)方法專門(mén)針對(duì)此類數(shù)據(jù),采用平滑方法將其擬合成曲線再進(jìn)行分析。Ramsay[5]于1982年率先提出一種全新的數(shù)據(jù)分析思路——函數(shù)型數(shù)據(jù)分析;Ramsay等[6]對(duì)函數(shù)型數(shù)據(jù)做了進(jìn)一步詳細(xì)的描述并講述了諸多關(guān)于FDA的應(yīng)用。與傳統(tǒng)方法相比,函數(shù)型數(shù)據(jù)分析方法不僅在處理高維觀測(cè)數(shù)據(jù)上能給出更加合理的直觀解釋,而且在分析數(shù)據(jù)時(shí)能保留更多的數(shù)據(jù)信息,從而得到更精確的分析結(jié)果。因此,F(xiàn)DA應(yīng)用于許多學(xué)科領(lǐng)域,如生物學(xué)、醫(yī)學(xué)、氣象學(xué)、計(jì)量經(jīng)濟(jì)學(xué)、金融、化學(xué)計(jì)量學(xué)和地球物理學(xué)等[7-8]。函數(shù)型回歸分析是FDA最重要的應(yīng)用之一,經(jīng)典函數(shù)型線性回歸模型(functional linear model,F(xiàn)LM)是函數(shù)型線性回歸的最簡(jiǎn)單形式。Paganoni等[9]給出了概述,在FDA研究領(lǐng)域受到廣泛關(guān)注。在對(duì)FLM進(jìn)行參數(shù)估計(jì)時(shí);Cardot等[10]首次研究函數(shù)型線性回歸模型,利用函數(shù)型主成分分析方法給出斜率函數(shù)估計(jì)的收斂速度;Hall等[11]基于函數(shù)型主成分分析方法和最小二乘估計(jì)方法給出函數(shù)型線性回歸模型的估計(jì)。在做FLM分析時(shí)假設(shè)所有個(gè)體都是相互獨(dú)立的,但是在信息技術(shù)的飛速發(fā)展今天,區(qū)域之間的相關(guān)性信息逐漸顯現(xiàn),具有空間結(jié)構(gòu)的數(shù)據(jù)越來(lái)越普遍。在這些數(shù)據(jù)中,響應(yīng)變量由于空間結(jié)構(gòu)而存在較強(qiáng)的相關(guān)性,如果使用FLM來(lái)分析帶有空間相關(guān)性的數(shù)據(jù),可能無(wú)法充分利用數(shù)據(jù)中包含的信息。當(dāng)協(xié)變量為標(biāo)量時(shí),通常使用空間自回歸(spatial autoregression,SAR)模型來(lái)擬合這類具有空間相關(guān)性的數(shù)據(jù)。SAR模型是處理空間相關(guān)問(wèn)題的一種有效方法。
Cliff等[12]在空間自回歸的一般模型、參數(shù)估計(jì)和假設(shè)檢驗(yàn)方面有了開(kāi)拓性的進(jìn)展;Anselin[13]提出了空間自回歸的一般模型。近年來(lái),SAR模型的應(yīng)用越來(lái)越廣泛;Kanaroglou等[14]將SAR模型用于估計(jì)SO2的空氣污染濃度.在環(huán)境學(xué)領(lǐng)域與氣象學(xué)領(lǐng)域中,觀測(cè)數(shù)據(jù)會(huì)有曲線形式。因此,將SAR模型的數(shù)據(jù)類型擴(kuò)展到函數(shù)型數(shù)據(jù)越來(lái)越迫切。Ahmed[15]首次提出了函數(shù)型空間自回歸模型(functional spatial autoregressive model,F(xiàn)SAR),并利用極大似然估計(jì)方法對(duì)FSAR模型進(jìn)行估計(jì),并證明了該模型中參數(shù)的漸進(jìn)性質(zhì);Huang等[16]通過(guò)借鑒SAR模型的概念,利用空間相關(guān)系數(shù)和權(quán)重矩陣將具有空間結(jié)構(gòu)的數(shù)據(jù)融入FLM模型中,稱該模型為空間函數(shù)型線性模型(spatial functional linear model,SFLM),并利用該模型分析了2005—2007年中國(guó)34個(gè)主要城市年均降水量與月均溫度的關(guān)系;應(yīng)彩云[17]將函數(shù)型變量引入一般的SAR模型中,構(gòu)建半泛函SAR模型,并研究相應(yīng)的參數(shù)估計(jì)方法。汾渭平原作為中國(guó)的四大平原之一,有著嚴(yán)重的大氣污染問(wèn)題,但目前還沒(méi)有學(xué)者函數(shù)型回歸模型分析該區(qū)域的空氣質(zhì)量數(shù)據(jù)與氣象數(shù)據(jù)之間的關(guān)系。從汾渭平原11個(gè)城市的空間分布來(lái)看,空氣質(zhì)量數(shù)據(jù)中必然有部分?jǐn)?shù)據(jù)存在空間相關(guān)性,結(jié)合SAR和FLM分析方法,采用FSAR模型探討該地區(qū)空氣質(zhì)量與氣象的關(guān)系不僅能充分利用數(shù)據(jù)的空間自相關(guān)和高維連續(xù)性特征,而且提供分析汾渭平原大氣污染情況的新途徑。
采用2019年汾渭平原11個(gè)城市的空氣質(zhì)量數(shù)據(jù)[PM2.5、PM10.0、SO2、NO2、O3和CO的監(jiān)測(cè)數(shù)據(jù)及空氣質(zhì)量指數(shù)(air quality index,AQI)],運(yùn)用空間自相關(guān)Moran’sI檢驗(yàn)對(duì)這7項(xiàng)數(shù)據(jù)做空間自相關(guān)分析,選擇檢驗(yàn)結(jié)果最顯著指標(biāo)作為響應(yīng)變量,并采用2019年該區(qū)域11個(gè)城市的逐小時(shí)氣溫,應(yīng)用函數(shù)型空間自回歸方法構(gòu)建空氣質(zhì)量數(shù)據(jù)與氣溫的FSAR模型。定量分析汾渭平原11個(gè)城市空氣質(zhì)量數(shù)據(jù)與氣溫的相關(guān)性,并與普通的函數(shù)型回歸模型進(jìn)行比較分析。
選取汾渭平原(包括河南省洛陽(yáng)市、三門(mén)峽市,陜西省西安市、咸陽(yáng)市、寶雞市、銅川市、渭南市,山西省呂梁市、晉中市、臨汾市、運(yùn)城市) 11 個(gè)市作為研究對(duì)象, 整理了2019年汾渭平原地區(qū)11個(gè)城市的7項(xiàng)空氣質(zhì)量數(shù)據(jù)[PM2.5、PM10.0、SO2、NO2、O3和CO的監(jiān)測(cè)數(shù)據(jù)及空氣質(zhì)量指數(shù)(AQI)]和氣象數(shù)據(jù)(地表氣溫),數(shù)據(jù)分別來(lái)自中國(guó)空氣質(zhì)量在線監(jiān)測(cè)分析平臺(tái)和美國(guó)國(guó)家航空航天局MERRA2中的inst1_2d_asm_Nx文件集。汾渭平原11個(gè)城市在空間上的分布如圖1所示。
圖1 汾渭平原地形圖Fig.1 Topographic map of Fenwei plain
FSAR的一般形式為
(1)
式(1)中:x(t)=[x1(t),x2(t),…,xn(t)]T,t∈T為函數(shù)型協(xié)變量,t表示某一時(shí)間段,T為整個(gè)時(shí)間段;y=(y1,y2,…,yn)T為標(biāo)量響應(yīng)變量;β(t)表示斜率函數(shù);ρ為空間相關(guān)系數(shù);W為空間權(quán)重矩陣;ε為殘差項(xiàng),并假設(shè)ε服從多元正態(tài)分布N(0,σ2I),其中σ2為殘差的方差,I為單位矩陣。
1.2.1 函數(shù)型主成分基展開(kāi)
假設(shè)總體為x(t),記K(s,t)=Cov[x(s),x(t)]為x(t)的協(xié)方差函數(shù),根據(jù)Mercer’s 引理有協(xié)方差函數(shù)的譜分解,可表示為[16]
(2)
第i個(gè)觀測(cè)值xi(t)的Karhumen-Loéve展開(kāi)式可表示為
(3)
在實(shí)際應(yīng)用中,φj(t)是未知的,因此通過(guò)樣本協(xié)方差函數(shù)的譜分解來(lái)估計(jì)特征函數(shù)φj(t),即考慮樣本協(xié)方差陣,可表示為
(4)
(5)
式(5)中:aj=(a1j,a2j,…,anj)T。
(6)
1.2.2 FSAR模型的極大似然估計(jì)
根據(jù)近似FSAR模型[式(6)],用極大似然方法估計(jì)未知參數(shù),令A(yù)=(aij)n×m,b=(b1,b2,…,bm)T,從而式(6)可寫(xiě)為
y≈ρWy+Ab+ε
(7)
基于前面的假設(shè)ε服從多元正態(tài)分布N(0,σ2I),y的似然函數(shù)為
(8)
式(8)中:e=y-ρWy-Ab,先給定ρ,可得到b和σ2的極大似然估計(jì)分別為
(9)
(10)
(11)
通過(guò)最大化式(11),可求得ρ的估計(jì)量為
(12)
(13)
分別對(duì)2019年汾渭平原地區(qū)11個(gè)城市的7項(xiàng)空氣質(zhì)量數(shù)據(jù)[PM2.5、PM10.0、SO2、NO2、O3和CO的監(jiān)測(cè)數(shù)據(jù)及空氣質(zhì)量指數(shù)(AQI)],運(yùn)用R、Geoda等軟件進(jìn)行空間自相關(guān)分析,采用rook矩陣為空間權(quán)重矩陣。首先對(duì)年均空氣質(zhì)量進(jìn)行空間自相關(guān)分析,年均空氣質(zhì)量的空間相關(guān)性是否顯著,由Moran’sI指數(shù)驗(yàn)證,結(jié)果如表1所示。
表1結(jié)果顯示:這7項(xiàng)空氣質(zhì)量數(shù)據(jù)中SO2濃度和O3具有較強(qiáng)的空間自相關(guān)性,且SO2濃度的空間自相關(guān)性非常顯著,P遠(yuǎn)小于0.05。選取7項(xiàng)空氣質(zhì)量數(shù)據(jù)中的SO2濃度數(shù)據(jù)與氣溫?cái)?shù)據(jù)進(jìn)行函數(shù)型空間自回歸分析。由于所研究的是2019年每月SO2濃度與該月氣溫之間的關(guān)系,故對(duì)每月的SO2濃度進(jìn)行Moran’sI檢驗(yàn),如表2所示,結(jié)果顯示每個(gè)月這11個(gè)城市的SO2濃度都具有很強(qiáng)的空間自相關(guān)性,這說(shuō)明運(yùn)用FASR模型分析方法是合理有效的。
表1 2019年汾渭平原11個(gè)城市的年均空氣 質(zhì)量Moran’s I值Table 1 Moran’s I value of annual average air quality of 11 cities in Fenwei Plain in 2019
表2 2019年汾渭平原11個(gè)城市SO2月均Moran’s ITable 2 Monthly average of Moran’s I in SO2 in 11 cities in the Fenwei Plain in 2019
選取汾渭平原11個(gè)城市2019年的逐時(shí)氣溫?cái)?shù)據(jù)和逐月SO2濃度數(shù)據(jù),觀察發(fā)現(xiàn)收集到的逐時(shí)氣溫?cái)?shù)據(jù)具有周期性,故采用傅里葉基樣條函數(shù)擬合2019年每月的逐時(shí)氣溫曲線xi(t),i=1,2,…,11,并將其作為函數(shù)型空間自回歸模型的協(xié)變量,平滑曲線如圖2所示,為驗(yàn)證原始數(shù)據(jù)與平滑曲線之間的擬合效果,計(jì)算各月擬合曲線的均方根誤差(root mean squared error,RMSE),結(jié)果如表3所示。
表3 各月份擬合曲線的均方根誤差Table 3 RMSE of the fitted curve for each month
從圖2可以看出,2019年汾渭平原11個(gè)城市的溫度變化趨勢(shì)大體相似,每日氣溫呈周期性變化。汾渭平原一年四季溫差較大,冬冷夏熱,最高溫出現(xiàn)在8月份,最低溫出現(xiàn)在1月份。1—7月溫度持續(xù)上升,8月達(dá)到最大,之后溫度又持續(xù)下降。通過(guò)導(dǎo)函數(shù)刻畫(huà)了氣溫的動(dòng)態(tài)變化特征,如一階導(dǎo)函數(shù)(圖3)及二階導(dǎo)函數(shù)(圖4)所示。
圖2 2019年每月逐時(shí)氣溫變化曲線Fig.2 Hourly temperature change curve of each month in 2019
圖4 11個(gè)城市每月溫度曲線的二階導(dǎo)函數(shù)Fig.4 The second derivative function of monthly temperature curve in 11 cities
圖3為溫度曲線的一階導(dǎo)函數(shù)圖,表示的是溫度的上升或下降狀態(tài),一階導(dǎo)函數(shù)大于零時(shí)溫度上升,小于零時(shí)溫度下降。圖4為溫度曲線的二階導(dǎo)函數(shù)圖,表示溫度上升或下降速度的快慢。從圖中可以發(fā)現(xiàn),11個(gè)城市每日的溫度上升和下降的變化相似,3、4、5月的一階導(dǎo)函數(shù)圖波動(dòng)較大,這說(shuō)明這3個(gè)月的每日溫度變化較大,且當(dāng)月的最低溫與最高溫相差30 ℃左右。
利用前面所述的函數(shù)型空間自回歸分析方法,基于R語(yǔ)言編寫(xiě)的程序?qū)崿F(xiàn)汾渭平原11個(gè)城市的SO2濃度與氣溫曲線的函數(shù)型空間自回歸實(shí)證分析。以2019年每月的逐時(shí)氣溫曲線xi(t),i=1,2,…,11為函數(shù)型空間自回歸模型的協(xié)變量,月均SO2濃度yi為響應(yīng)變量,建立函數(shù)型空間自回歸(FSAR)模型為
(14)
式(14)中:yi′為第i′個(gè)城市的月均SO2濃度;wii′為城市i和i′之間的權(quán)重;εi為第i個(gè)等式對(duì)應(yīng)的殘差。
為了分析汾渭平原11個(gè)城市之間的SO2濃度存在空間自相關(guān)性,同時(shí)考慮普通的函數(shù)型線性模型(functional linear model,F(xiàn)LM)為
(15)
式(15)中:響應(yīng)變量yi和函數(shù)型協(xié)變量xi(t)分別表示2019年第i個(gè)城市的月均SO2濃度和該月對(duì)應(yīng)的逐時(shí)氣溫曲線;β(t)、γ(t)為斜率函數(shù)。
為了估計(jì)式(14)、式(15)中的斜率函數(shù)β(t)、γ(t),使用函數(shù)型主成分分析降維, 基于主成分基函數(shù)的FSAR模型和FLM分別表示為
(16)
(17)
首先利用函數(shù)型主成分分析方法對(duì)2019年汾渭平原11個(gè)城市12個(gè)月的氣溫變化曲線分別提取了前4個(gè)主成分,特征值及其累積方差貢獻(xiàn)率如表4所示。
根據(jù)表4發(fā)現(xiàn)當(dāng)取前4個(gè)主成分時(shí)累積方差貢獻(xiàn)率均在94%以上,包含了數(shù)據(jù)足夠多的信息?;诤瘮?shù)型主成分分析對(duì)氣溫函數(shù)完成降維,將無(wú)限維的FSAR模型轉(zhuǎn)換為了空間自回歸模型和函數(shù)型線性模型轉(zhuǎn)換為了多元線性模型,可表示為
表4 函數(shù)型主成分分析特征值及其累積貢獻(xiàn)率Table 4 Functional principal component analysis characteristic values and their cumulative contribution rate
(18)
(19)
采用1.2節(jié)介紹的估計(jì)方法進(jìn)行函數(shù)型空間自回歸分析,并對(duì)β(t)進(jìn)行估計(jì),斜率函數(shù)的估計(jì)結(jié)果如圖5所示。
圖5 每月FSRA模型的斜率函數(shù)β(t)Fig.5 Slope function β(t) of monthly FSRA model
觀察圖5可知,每月β(t)為負(fù)的區(qū)域均大于β(t)為正的區(qū)域,這說(shuō)明從整體上看氣溫與SO2濃度呈負(fù)相關(guān),并結(jié)合圖3每月氣溫一階導(dǎo)函數(shù)圖像可以發(fā)現(xiàn),氣溫上升時(shí),對(duì)SO2濃度的負(fù)影響較小,氣溫下降時(shí),對(duì)SO2濃度的負(fù)影響較大,2月的β(t)為負(fù)的面積在全年中是最小的,其次是1月、8月的β(t)為負(fù)的區(qū)域在全年中是最大的,其次是7月。從每天的溫度變化趨勢(shì)來(lái)看,氣溫每日呈周期性變化,其對(duì)SO2濃度的負(fù)影響也呈周期性變化,在每日氣溫上升的時(shí)間段內(nèi),對(duì)SO2濃度的負(fù)影響比較小,在溫度下降的時(shí)間段內(nèi),對(duì)SO2濃度的負(fù)影響較大。
根據(jù)表5均方誤差(mean absolute error,MSE)可知,F(xiàn)SAR模型擬合的響應(yīng)變量MSE明顯低于FLM的MSE,這說(shuō)明FSAR的擬合效果更好,這兩種方法的MSE平均相差7倍以上,12月份的差別最大。從圖6可以看出,F(xiàn)SAR與FLM之間的差異,同時(shí)觀察到FLM模型的MSE在全年呈現(xiàn)明顯的季節(jié)性差異,在冬季MSE很大,夏季很小,很大一部分原因是由于汾渭平原冬季供暖[19]這一人為因素造成,而FSAR并沒(méi)有很明顯的這類變化規(guī)律,這說(shuō)明FSAR模型降低了人為因素的誤差影響。冬季FSAR模型效果明顯,而FLM擬合的MSE 全年遠(yuǎn)大于20,這說(shuō)明FSAR 模型能夠很好地反映汾渭平原11 個(gè)城市的2019每月SO2濃度之間的空間自相關(guān)性,而FLM不能有效地解決這種相關(guān)性。可見(jiàn),面對(duì)具有空間屬性的數(shù)據(jù)時(shí),函數(shù)型空間自回歸模型擬合效果更好。
表5 FSAR模型與FLM模型的響應(yīng)變量均方誤差比較Table 5 Comparison of mean square error of response variables between FSAR model and FLM model
圖6 每月FSAR模型與FLM模型SO2濃度均方 誤差直方圖Fig.6 Monthly mean square error histogram of SO2 concentration in FSAR model and FLM model
借助函數(shù)型數(shù)據(jù)分析的思想并與空間自回歸模型結(jié)合,針對(duì)汾渭平原特殊的地理位置和空氣質(zhì)量數(shù)據(jù)與氣象數(shù)據(jù)的關(guān)系,運(yùn)用函數(shù)型空間自回歸模型分析該區(qū)域2019每月SO2濃度與當(dāng)月氣溫曲線的關(guān)系。由于氣溫?cái)?shù)據(jù)是每小時(shí)觀測(cè)一次,數(shù)據(jù)觀測(cè)頻率高且密集,而SO2濃度數(shù)據(jù)為月均數(shù)據(jù),具有顯著的空間相關(guān)性,為了不丟失數(shù)據(jù)中包含的大量信息,將溫度數(shù)據(jù)擬合成函數(shù)曲線,把它看成隨機(jī)向量,此時(shí)選擇函數(shù)型線性模型和函數(shù)型空間自回歸模型來(lái)處理這組數(shù)據(jù)。通過(guò)函數(shù)型數(shù)據(jù)分析與空間自回歸模型結(jié)合的方式對(duì)汾渭平原的SO2與氣溫進(jìn)行函數(shù)型空間自回歸分析,得到如下結(jié)論。
(1)從空間自相關(guān)分析中可以看出,汾渭平原地區(qū)各市的SO2濃度和氣溫除了受本市影響之外,還均受相鄰城市SO2濃度的影響,甚至比本市的影響更大。
(2)在運(yùn)用函數(shù)型數(shù)據(jù)分析方法對(duì)逐時(shí)氣溫?cái)?shù)據(jù)進(jìn)行擬合成氣溫曲線后,對(duì)氣溫曲線求一階導(dǎo)數(shù)和二階導(dǎo)數(shù),一階導(dǎo)數(shù)顯示汾渭平原氣溫呈周期性上升和下降,個(gè)別城市出現(xiàn)氣溫急劇變化情況。
(3)通過(guò)FSAR模型中系數(shù)函數(shù)β(t)的圖像可以看出,氣溫曲線與SO2濃度之間有顯著的負(fù)相關(guān)性。對(duì)比函數(shù)型空間自回歸模型和函數(shù)型線性模型,可以發(fā)現(xiàn)FSAR模型能夠很好地反映汾渭平原11個(gè)城市的2019每月SO2濃度之間的空間自相關(guān)性,擬合效果更好。
(4)函數(shù)型空間自回歸模型不僅能在宏觀上分析全年氣溫的變化對(duì)SO2濃度,而且在微觀上針對(duì)每小時(shí)的氣溫變化對(duì)SO2濃度的影響也清晰可見(jiàn)。并且通過(guò)對(duì)比FLM和FSAR的MSE,發(fā)現(xiàn)FSAR模型大大降低了人為因素所帶來(lái)的誤差。
氣溫曲線和SO2濃度是研究大氣污染的兩個(gè)重要指標(biāo)。隨著收集數(shù)據(jù)的方法越來(lái)越精細(xì)化,學(xué)者們不僅在意宏觀上數(shù)據(jù)信息,更重視微觀上的分析所帶來(lái)更多更精確的信息,因此高頻數(shù)據(jù)的處理方法變得至關(guān)重要,函數(shù)型數(shù)據(jù)分析方法也受到人們的普遍關(guān)注。將函數(shù)型數(shù)據(jù)與帶有地理信息的數(shù)據(jù)結(jié)合,分析它們的空間分布特征和規(guī)律,函數(shù)型空間自回歸將成為近代氣候工作者關(guān)心的一個(gè)重要方向。