胡娟 鄭軍 許文龍 邱玉珺 胡方超

(1 南京信息工程大學(xué) 環(huán)境科學(xué)與工程學(xué)院，南京 210044；2 南京信息工程大學(xué) 網(wǎng)絡(luò)信息中心，南京 210044； 3 南京信息工程大學(xué) 大氣物理學(xué)院，南京 210044；4 常州市環(huán)境科學(xué)研究院，江蘇 常州 213022)

引 言

PM2.5是指懸浮在空中的空氣動力學(xué)當(dāng)量直徑小于或等于2.5 μm的細顆粒物，嚴重影響環(huán)境、氣候，甚至危害人體健康[1-2]。盡管環(huán)境監(jiān)測站對近地面PM2.5直接測量的精度很高，但在監(jiān)測樣本點的數(shù)量與分布上有一定的局限性，難以反映PM2.5的時空變化，而利用遙感衛(wèi)星測量則可對其進行大范圍的區(qū)域監(jiān)測。目前諸多研究是在衛(wèi)星遙感反演產(chǎn)品AOD (Aerosol Optical Depth)的基礎(chǔ)上對PM2.5質(zhì)量濃度進行估計[3]?，F(xiàn)今衛(wèi)星探測傳感器可為大氣氣溶膠監(jiān)測提供數(shù)據(jù)產(chǎn)品，尤以MODIS的氣溶膠光學(xué)厚度AOD L2產(chǎn)品Terra MOD04與Aqua MYD04為典型，因其反演算法成熟、產(chǎn)品質(zhì)量可靠、時間連續(xù)性好，所以應(yīng)用最為廣泛[4]，特別是在大氣環(huán)境監(jiān)測領(lǐng)域。

利用衛(wèi)星遙感技術(shù)對PM2.5進行長時間、大范圍、高空間分辨率的監(jiān)測也是當(dāng)前大氣污染治理的一項重要內(nèi)容。已有研究表明，AOD與近地面PM2.5具有強相關(guān)性[5]，利用衛(wèi)星估算近地面PM2.5質(zhì)量濃度幾乎都是基于AOD展開的[6]，主要有線性或多元回歸模型[7-8]、地理加權(quán)回歸及神經(jīng)網(wǎng)絡(luò)的高級統(tǒng)計模型[9-10]。一般而言，衛(wèi)星AOD包括實際大氣整層信息，需結(jié)合氣象數(shù)據(jù)或再分析資料做相應(yīng)的高度與濕度訂正，這樣與PM2.5的相關(guān)性更高。

由于PM2.5的時空變化會受到大氣環(huán)境、氣象、地理等諸多因素的影響，隨著近年來機器學(xué)習(xí)技術(shù)的發(fā)展，已有學(xué)者利用BP(Back Propagation)神經(jīng)網(wǎng)絡(luò)[11]、隨機森林RF(Random Forest)[12-14]、支持向量機SVM(Support Vector Machine)[15]等方法來估算中國東部PM2.5[16]。機器學(xué)習(xí)模型既不需詳細了解大氣污染的物理化學(xué)過程[17]，也不需對衛(wèi)星遙感的機理進行更深入的研究，只要根據(jù)統(tǒng)計學(xué)原理，利用機器學(xué)習(xí)的特定算法即可。理論上，輸入的信息越多，機器學(xué)習(xí)的效果越好。由于地面PM10包含PM2.5，兩者明顯正相關(guān)，且常常是通過AOD一起估算，而能見度(VIS)數(shù)據(jù)包括大氣成分、相對濕度等多因素，其與PM10很少被用來估算PM2.5，但能見度也代表了近地面信息，與AOD聯(lián)合可用來表示氣溶膠的標(biāo)高，但近年來已有研究在機器學(xué)習(xí)算法中應(yīng)用VIS與PM10來估算PM2.5[14，18-19]。

因此，通過機器學(xué)習(xí)進行衛(wèi)星遙感反演PM2.5是一種更具應(yīng)用潛力的監(jiān)測方法?；谶z傳算法的GA-BP(Genetic Algorithm-Back Propagation)神經(jīng)網(wǎng)絡(luò)作為BP神經(jīng)網(wǎng)絡(luò)的發(fā)展，被用于估算PM2.5的研究還不多見[20]。本文嘗試結(jié)合地面站點的氣象、環(huán)境監(jiān)測數(shù)據(jù)，通過GA-BP機器學(xué)習(xí)算法來訓(xùn)練表達衛(wèi)星遙感反演得到的AOD與近地面PM2.5之間復(fù)雜的非線性關(guān)系，從而估計PM2.5的質(zhì)量濃度。

1 數(shù)據(jù)來源及算法原理

1.1 研究區(qū)域

選取南京地區(qū)各站點所代表的區(qū)域作為研究區(qū)域，研究區(qū)域內(nèi)的氣象觀測站點(簡稱氣象站，下同)與環(huán)境監(jiān)測站點(簡稱環(huán)保站，下同)分布如圖1所示。在進行數(shù)據(jù)站點匹配處理時，按照最近距離原則，即每一個環(huán)保站總是要找到一個最近的氣象站與之匹配對應(yīng)。

圖1 南京地面站點分布(單位：km)

1.2 數(shù)據(jù)來源

衛(wèi)星遙感數(shù)據(jù)來自美國地球觀測系統(tǒng)EOS(Earth Observation System)系列衛(wèi)星上的Terra和Aqua衛(wèi)星所攜帶的MODIS提供的Level 2產(chǎn)品數(shù)據(jù)AOD，上午為Terra MOD04_3 km， 下午為Aqua MYD04_3 km，地面分辨率均為3 km，每日上、下午各一次，可提取有關(guān)日期(Date)、經(jīng)度(Lon，(°))、緯度(Lat，(°))及AOD等信息。

氣象數(shù)據(jù)為中國地面氣象站逐小時觀測資料，包括日期(Date)、站號(Station_Id_C)、緯度(Lat，(°))、經(jīng)度(Lon，(°))、氣溫(TEM，℃)、氣壓(PRS，hPa)、相對濕度(RHU，%)、2 min平均風(fēng)向(角度)(WIN_D_Avg_2 min，(°))、2 min平均風(fēng)速(WIN_ S_Avg_2 min，m·s-1)、能見度(VIS，m)、天氣特征(WEP)。其中WEP表示測站所觀測到的各種天氣現(xiàn)象，變化范圍從00～198。

國家環(huán)保部地面監(jiān)測站點數(shù)據(jù)來自中國環(huán)境監(jiān)測總站的全國城市空氣質(zhì)量實時發(fā)布平臺，使用其逐小時空氣質(zhì)量數(shù)據(jù)，包括日期(Date)、PM2.5、PM10質(zhì)量濃度等。

1.3 研究方法

1.3.1 BP神經(jīng)網(wǎng)絡(luò)工作原理

BP神經(jīng)網(wǎng)絡(luò)(Back Propagation neural network)由Rumelhart，et al[21]提出，在神經(jīng)網(wǎng)絡(luò)訓(xùn)練中引入了δ學(xué)習(xí)規(guī)則(用于調(diào)整各層神經(jīng)元的權(quán)值和閾值的方法)，是一種按誤差反向傳播算法訓(xùn)練的多層前向反饋神經(jīng)網(wǎng)絡(luò)，由信息的正向傳播和誤差的反向傳播兩個過程組成。

正向傳播過程指的是輸入層各節(jié)點接收輸入信息(如本文輸入AOD、能見度、相對濕度、氣溫等)，并將信息傳遞給隱含層的各個節(jié)點，最后由輸出層向外界輸出信息處理結(jié)果(圖2)。當(dāng)實際輸出與期望輸出(本文輸出為PM2.5質(zhì)量濃度)沒達到最小均方誤差(Mean Square Error， MSE)條件時，則進入誤差的反向傳播，即誤差通過輸出層，按誤差梯度下降的方式修正各層權(quán)值，向隱含層、輸入層逐層再反向傳遞。就這樣循環(huán)不斷地調(diào)整各層權(quán)值和閾值，進行網(wǎng)絡(luò)學(xué)習(xí)訓(xùn)練一直到MSE減少至可以接受的條件，或?qū)W習(xí)次數(shù)達到預(yù)設(shè)值為止。

圖2 BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)

由于BP神經(jīng)網(wǎng)絡(luò)中的初始權(quán)重和閾值一般是隨機的，且梯度下降法的自變量范圍被限制在單調(diào)區(qū)間內(nèi)，這樣容易導(dǎo)致陷入局部最優(yōu)解而不是全局最優(yōu)解，而且還存在學(xué)習(xí)收斂速度慢、泛化能力弱、存在欠擬合等缺陷[22]。

1.3.2 GA-BP網(wǎng)絡(luò)算法

遺傳算法(Genetic Algorithm，GA)是Holland[23]基于生物進化規(guī)律演變而得到的一種優(yōu)化算法，它對復(fù)雜的、無明確關(guān)系的、多變量的、非線性的逼近效率較高，且能夠按照適應(yīng)度函數(shù)進行智能搜索以求得全局最優(yōu)解。而遺傳算法與BP神經(jīng)網(wǎng)絡(luò)相結(jié)合的GA-BP算法采用優(yōu)化的初始化權(quán)值與閾值，能夠在整個定義域內(nèi)實現(xiàn)變量取值最大可能地接近全局最優(yōu)解。

1.4 GA-BP優(yōu)化算法建模

本文采用GA-BP算法結(jié)合MODIS AOD、地面氣象站數(shù)據(jù)包括氣溫(℃)、氣壓(hPa)、相對濕度(%)、2 min平均風(fēng)向(°)、2 min平均風(fēng)速(m·s-1)、能見度(m)，以及地面環(huán)保站點PM10等數(shù)據(jù)建模估算地面PM2.5質(zhì)量濃度(μg·m-3)。利用GA-BP神經(jīng)網(wǎng)絡(luò)訓(xùn)練，不斷更新BP網(wǎng)絡(luò)的權(quán)值和閾值直至滿足條件判斷并輸出測試結(jié)果。GA-BP神經(jīng)網(wǎng)絡(luò)流程如圖3所示。

圖3 GA-BP神經(jīng)網(wǎng)絡(luò)流程

1.4.1 原始數(shù)據(jù)采集及時空匹配數(shù)據(jù)源(訓(xùn)練集和測試集)

其中，時空匹配站點，獲取數(shù)據(jù)源包括：

(1)地面環(huán)保站點與氣象站點時空匹配

將兩者按時間和空間(站點經(jīng)緯度最近距離)匹配，匹配后用經(jīng)緯度標(biāo)識(圖1)。表1為匹配結(jié)果。選取六合、儀征兩個氣象站點，9個環(huán)保站點，其中5個用于訓(xùn)練，4個用于測試，如表1所示。

表1 國家環(huán)保站點與氣象站點時空匹配數(shù)據(jù)

(2)地面站點與MODIS AOD時空匹配

表1中的中點經(jīng)度、中點緯度是匹配后的環(huán)保站與氣象站兩站點之間中點的經(jīng)緯度。對在此中點30 km范圍內(nèi)， 15 min時間內(nèi)相應(yīng)的MODIS AOD數(shù)據(jù)進行匹配，所用到的地表兩點間(大圓)距離公式為：

L=Re×arccos[sin(lat1)sin(lat2)+cos(lat1)cos(lat2)cos(lon1-lon2)]

(1)

其中：地球半徑Re=6 378.137 km；(lat1，lon1)、(lat2，lon2)分別為A、B兩點的經(jīng)緯度。

(3)PM2.5變量因子的相關(guān)性分析

不同時空范圍內(nèi)，AOD、PM10及氣象條件因子和PM2.5質(zhì)量濃度的相關(guān)性也存在差異，表2為對各站點某些變量因子與PM2.5作相關(guān)分析確定的網(wǎng)絡(luò)輸入層的特征因子。

表2 PM2.5與PM10、AOD和氣象條件的相關(guān)性分析(相關(guān)系數(shù)R)

(4)確定輸入變量組合試驗

為了在模型中考查多源數(shù)據(jù)中地面數(shù)據(jù)的輸入變量對衛(wèi)星AOD輸入變量的影響，構(gòu)建了9個組合試驗(表3)，分別用試驗1，2，…，9表示。

表3 9組試驗的各輸入變量

1.4.2 GA-BP神經(jīng)網(wǎng)絡(luò)訓(xùn)練及仿真測試過程

在網(wǎng)絡(luò)訓(xùn)練前首先要對樣本數(shù)據(jù)統(tǒng)一進行歸一化(Normalization)處理，然后隨機選取樣本進行網(wǎng)絡(luò)訓(xùn)練，經(jīng)過多次訓(xùn)練確定網(wǎng)絡(luò)優(yōu)化參數(shù)。將網(wǎng)絡(luò)估算的PM2.5質(zhì)量濃度和實際值的誤差平方和的倒數(shù)作為遺傳算法的個體適應(yīng)度函數(shù)，通過選擇、交叉、變異最終確定適應(yīng)值最優(yōu)的個體，即確定最優(yōu)的權(quán)值和閾值并將其賦值于BP神經(jīng)網(wǎng)絡(luò)，再經(jīng)過網(wǎng)絡(luò)訓(xùn)練進行仿真測試，并將輸出數(shù)據(jù)反歸一化處理，得到網(wǎng)絡(luò)模型估算的PM2.5。

1.4.3 BP及GA-BP模型性能評價

性能指標(biāo)定義及公式，包括：

(1)決定系數(shù)(Coefficient of Determination)：

(2)

其中:E表示估算值;T表示實際值。

(2)相對誤差(ei)、平均絕對誤差(Mean Absolute Error， MAE)：

ei=|E-T|/T

(3)

(4)

(3)均方根誤差(Root Mean Squared Error， RMSE)，用來衡量預(yù)測值與實際值之間的偏離程度：

(5)

(4)準(zhǔn)確率(Accuracy)：把估算值的相對誤差超過其實際值的20%視為不準(zhǔn)確數(shù)據(jù)。

選取樣本采用隨機選取每個匹配站點數(shù)據(jù)的90%進行神經(jīng)網(wǎng)絡(luò)訓(xùn)練，并用匹配數(shù)據(jù)的10%進行仿真測試，并根據(jù)多次訓(xùn)練仿真測試保存訓(xùn)練最好的網(wǎng)絡(luò)用于南京其他站點的驗證試驗，即在9個環(huán)保站中選用5個站點訓(xùn)練選出最好的網(wǎng)絡(luò)來估算其他4個測試站點的PM2.5。

2 結(jié)果與討論

2.1 影響PM2.5估計的各相關(guān)因子分析

以南京各站點為例進行分析，試驗利用皮爾遜相關(guān)系數(shù)估計變量x、y間的相關(guān)。公式如下：

(6)

就南京地區(qū)2017年和2018年的數(shù)據(jù)分析來看(表2)，各變量因子與PM2.5相關(guān)性系數(shù)R按平均值大小依次為： PM10(0.771)、WEP(0.467)、PRS(0.312)、AOD(0.186)、WIN_D_Avg_2 min(0.061)、RHU(0.008)、WIN_S_Avg_2 min(-0.135)、TEM(-0.370)、VIS(-0.625)。從相關(guān)系數(shù)絕對值的大小劃分相關(guān)強度，PM10、VIS為中度相關(guān)，WEP、PRS、TEM、AOD為低度相關(guān)，其他因子相關(guān)相對較弱。不同的年度、季節(jié)和月份相關(guān)性也存在一定的差別，同一站點PM10與PM2.5相關(guān)性，2017年要明顯高于2018年。表2中WEP多數(shù)為無云、有霧或霾的情況，但由于2017年南京夏季W(wǎng)EP匹配后的數(shù)據(jù)全為0(表示未觀測或觀測不到云的發(fā)展)，故當(dāng)訓(xùn)練模型采取分站點分季節(jié)訓(xùn)練時，不再考慮WEP作為變量因子，最多使用8個輸入變量作為變量因子，如表3試驗6所示。

由于受到匹配樣品數(shù)量的限制，需要把數(shù)據(jù)分為同站點不同年度、不同季節(jié)進行網(wǎng)絡(luò)訓(xùn)練、仿真測試和驗證試驗。

2.2 BP與GA-BP網(wǎng)絡(luò)模型性能比較

2.2.1 BP網(wǎng)絡(luò)與GA-BP網(wǎng)絡(luò)訓(xùn)練性能指標(biāo)評價結(jié)果比較

對南京各環(huán)保站和氣象站按最近空間距離匹配，站點匹配結(jié)果并按“環(huán)保站_氣象站”的格式來表示。匹配后的網(wǎng)絡(luò)訓(xùn)練站點(5組)分別為：邁皋橋_儀征、玄武湖_儀征、草場門_六合、瑞金路_儀征、中華門_六合。匹配后的試驗測試站點(4組)分別為：山西路_六合、浦口_六合、奧體中心_六合、仙林大學(xué)城_儀征。

在輸入數(shù)據(jù)中使用“試驗6”的變量組合，對BP與GA-BP的訓(xùn)練結(jié)果比較見表4，可以看出，GA-BP整體上比BP(除了春季外)性能表現(xiàn)優(yōu)越。盡管在春季時，BP表現(xiàn)較GA-BP好，但在其他大量試驗中，春季的BP常出現(xiàn)比GA-BP表現(xiàn)欠佳的情況。

表4 BP與GA-BP網(wǎng)絡(luò)訓(xùn)練及仿真測試的主要性能評價指標(biāo)值

2.2.2 網(wǎng)絡(luò)訓(xùn)練及仿真測試效果

同上也選擇試驗6組合，對BP網(wǎng)絡(luò)和GA-BP網(wǎng)絡(luò)預(yù)測樣本的估算值與實際值進行比較，分析兩個網(wǎng)絡(luò)的訓(xùn)練及仿真測試的效果，如圖4所示?？梢钥闯?，2018年除了春季外，其他3個季節(jié)GA-BP的R2與RMSE表現(xiàn)均較BP好，尤其是夏秋兩季，特別是在冬季匹配數(shù)據(jù)很少的情況下，BP幾乎不能起到估計的作用，而GA-BP網(wǎng)絡(luò)R2與RMSE的表現(xiàn)，其優(yōu)勢則比較明顯。故下文就采用GA-BP網(wǎng)絡(luò)設(shè)計做進一步的試驗。

圖4 2018年南京匹配站點樣本BP與GA-BP估算PM2.5結(jié)果對比：(a)春季；(b)夏季；(c)秋季；(d)冬季

2.2.3 試驗1—9結(jié)果比較

圖5a、b可以看出，從試驗1到試驗5，每組試驗均增加一個變量因子，每增加一個變量因子時，對輸出PM2.5的R2與RMSE，其影響是不同的。在1—6組試驗中，與其他5組試驗相比，試驗1中AOD的R2較小，RMSE較大。這主要是因為AOD表征的是大氣整層的氣溶膠的影響，而其他變量因子均是在近地面幾個參數(shù)的測量，均一定程度上反映了近地面的信息，所以與PM2.5有著更好的相關(guān)性。在試驗2—5中，由于變量因子PM10的加入，各個季節(jié)試驗2中R2較試驗1有大幅度地增加，而在春季試驗5中由于TEM的加入使R2較試驗4略有下降；在秋季由于PRS(氣壓)加入，試驗3中R2較試驗2也略有下降；在春季試驗4中VIS的加入，較試驗3 2018年的R2增加較大，而2017年時的R2卻略有降低；在試驗6中，除了2017年秋季，加入相對濕度與平均風(fēng)速、平均風(fēng)向后，R2均有所增加。作為試驗1—6的對照試驗，增加了試驗7—9，分別考慮沒有VIS與PM10及兩者都沒有的情況。但從試驗7—9的仿真測試的效果看，對比試驗6，除了2017年秋冬季R2略有增大外，其他季節(jié)大多是R2減小RMSE增大。

圖5 南京站點9組試驗各季節(jié)R2與RMSE的變化：(a)2017年R2；(b)2018年R2；(c)2017年RMSE；(d)2018年RMSE

總體來看，對于RMSE幾乎都隨著R2的增大而減小。四季中秋季相關(guān)性最高，2017年的R2與RMSE總體要好于2018年，而且在9組試驗中試驗6總體也表現(xiàn)較好。

以仙林大學(xué)城_儀征的匹配站點為例，僅選用了試驗1與試驗6進行比較，如圖6所示?？芍?圖6A、C為試驗1， 圖6B、D為試驗6 ，2017年的R2與RMSE均是試驗6較好；四季表現(xiàn)為秋季>冬季>夏季>春季，秋季R2最大為0.91，RMSE=11.79;春季最小R2=0.65，RMSE=8.67。而圖6C、D分別為2018年試驗1與試驗6分季節(jié)情況，2018年的R2與RMSE也是試驗6較好，但四季表現(xiàn)為秋季>春季>冬季>夏季，秋季的最大R2=0.87，RMSE=7.50;夏季最小R2=0.56，RMSE=6.46?？傮w上看2018年R2雖較2017年小，但RMSE也較小。

圖6 仙林大學(xué)城_儀征站點試驗1 (A：2017；C：2018)與試驗6(B：2017；D：2018)網(wǎng)絡(luò)測試效果(a—d分別為春夏秋冬；橫坐標(biāo)為估算樣本的時間序列)

圖5、6均表明，試驗6的效果總體上要較前面試驗表現(xiàn)較好，這是因為，對GA-BP而言，輸入相關(guān)的變量因子越多，用于網(wǎng)絡(luò)學(xué)習(xí)的信息就越多，對諸多變量因子所組成的非線性關(guān)系就表達的越好，對PM2.5的估計效果也較明顯。

3 結(jié)論

(1)在有足夠數(shù)據(jù)樣本的情況下，不同季節(jié)的表現(xiàn)有時差別很大，需要進行分季節(jié)考慮。

(2)影響PM2.5估計的多源數(shù)據(jù)的各輸入變量因子中，PM10與VIS(能見度)變量為相對高相關(guān)因子，其他輸入變量 PRS(氣壓)、TEM(氣溫)、AOD等為低度相關(guān)，并在此基礎(chǔ)上構(gòu)建了9組基于衛(wèi)星AOD作為基本輸入變量的多源數(shù)據(jù)分季節(jié)估算PM2.5的試驗，試驗6的整體表現(xiàn)較好。

(3)利用GA-BP算法估計PM2.5比BP算法總體表現(xiàn)較好，所建立的多源數(shù)據(jù)GA-BP算法模型也為估算PM2.5的提供了一種新的機器學(xué)習(xí)的方法。