申 原 陳朝亮 錢 靜 劉 軍

(中國科學(xué)院深圳先進(jìn)技術(shù)研究院 深圳 518055)

1 引 言

隨著我國工業(yè)化、城市化進(jìn)程不斷加快，環(huán)境污染越來越嚴(yán)重，與居民生活密切相關(guān)的空氣質(zhì)量指數(shù)成為社會關(guān)注的熱點話題。細(xì)顆粒物(PM2.5)是衡量空氣質(zhì)量的重要參數(shù)指標(biāo)，采用科學(xué)的方法監(jiān)測 PM2.5的分布和濃度，對研究其本身的理化特性、揭示霧霾成因及采取正確的防治措施具有重要意義。

目前，PM2.5的監(jiān)測方法主要包括地面監(jiān)測和衛(wèi)星遙感監(jiān)測兩種[1]。雖然地面監(jiān)測結(jié)果較精確，但由于成本較高、地面監(jiān)測站數(shù)量少，導(dǎo)致監(jiān)測結(jié)果時空不連續(xù)，無法獲得足夠多的數(shù)據(jù)來研究整個區(qū)域 PM2.5的擴(kuò)散方式和傳輸特性[2,3]。衛(wèi)星遙感監(jiān)測具有數(shù)據(jù)獲取方便、監(jiān)測范圍廣等優(yōu)勢，能很好地彌補地面監(jiān)測的不足?，F(xiàn)有的 PM2.5監(jiān)測反演方法都是先反演大氣氣溶膠光學(xué)厚度(Aerosol Optical Depth，AOD)，再對氣溶膠光學(xué)厚度與地面實測 PM2.5的關(guān)系進(jìn)行統(tǒng)計分析，用統(tǒng)計得到的關(guān)系來推出無地面監(jiān)測點區(qū)域的 PM2.5值。國內(nèi)外很多學(xué)者用此方法進(jìn)行了大量的研究，其基本假設(shè)是，AOD 與 PM2.5具有良好的、穩(wěn)定的統(tǒng)計關(guān)系。因此，現(xiàn)有大量研究都集中在提高 AOD 反演精度上，如引入各種訂正、加入更多輔助數(shù)據(jù)、結(jié)合數(shù)值預(yù)報模式等，也在特定研究區(qū)域內(nèi)取得了比較好的結(jié)果。如陳輝等[4]利用地理加權(quán)回歸模型建立了我國冬季的 AOD-PM2.5模型；王子峰[5]系統(tǒng)地研究了衛(wèi)星遙感估算近地面顆粒物濃度的算法；王中挺等[6]利用暗目標(biāo)和高分一號衛(wèi)星 16m 相機數(shù)據(jù)反演了京津地區(qū)的氣溶膠光學(xué)厚度；王靜等[7]研究了北京市中分辨率成像光譜儀(Moderate Resolution Imaging Spectroradiometer，MODIS)氣溶膠光學(xué)厚度和 PM2.5質(zhì)量濃度的特征及其相關(guān)性；Song等[8]利用 MODIS C5AOD 產(chǎn)品反演出了珠江三角洲地區(qū)的 PM2.5濃度；張?zhí)炱錥9]則驗證了基于MODIS 影像的中國地區(qū)氣溶膠產(chǎn)品與 PM2.5反演的關(guān)系；Chu 等[10]基于 MODIS AOD 產(chǎn)品，反演出中國臺灣北部的 PM2.5濃度等。但由于在反演氣溶膠光學(xué)厚度 AOD 過程中會產(chǎn)生誤差，用帶有誤差的 AOD 與地面實測 PM2.5建立統(tǒng)計關(guān)系時會導(dǎo)致誤差的傳遞，從而影響最終 PM2.5的反演精度。

在反演 AOD 方面，一般用 MODIS 數(shù)據(jù)?；舅悸肥峭ㄟ^假設(shè)不同的氣溶膠模式和觀測幾何狀況，再計算氣溶膠光學(xué)厚度與大氣下界的半球反射率、大氣反射率和衛(wèi)星天頂角與太陽天頂角的 cos 余弦值之間的對應(yīng)關(guān)系。據(jù)此建立查找表，通過動態(tài)氣溶膠模式的輸入來查算氣溶膠光學(xué)厚度。Griggs[11]在大氣層頂平行且無云等假設(shè)前提下，根據(jù)模擬大氣輻射傳輸模型，發(fā)現(xiàn)了 AOD 與紅外波段和可見光波段的相關(guān)關(guān)系。Levy 等[12]將中紅外波段的氣溶膠信息加入到反演過程中，將紅藍(lán)兩波段獨立反演改進(jìn)為紅、藍(lán)、中紅外三波段同時反演，并更新了原有的大氣氣溶膠模式和 AOD 反演查找表。除了 MODIS數(shù)據(jù)，也可借助于其他數(shù)據(jù)反演氣溶膠光學(xué)厚度，如 Holben 等[13]用先進(jìn)超高分辨率輻射計數(shù)據(jù)和對比法反演了馬里薩赫勒地區(qū)的 AOD，反演誤差在 0.1左右；Isakov 等[14]利用機載可見紅外成像光譜儀數(shù)據(jù)，反演了美國俄克拉荷馬州和拉皮德城兩個地點的 AOD，結(jié)果發(fā)現(xiàn)，當(dāng)?shù)乇矸瓷渎什町惔笥?0.5時，AOD 反演精度可達(dá)0.9。上述反演方法在實際運用中均取得了較好的成果，但由于用到了很多輔助數(shù)據(jù)，使得計算精度難以控制。另外，AOD 數(shù)據(jù)在反演過程中本身就存在誤差，用上述方法無法避免誤差的傳遞[15]。因此，如何減小誤差、獲得更高的反演精度，一直是近年的研究熱點[16,17]。

但是，目前 PM2.5遙感反演方法存在以下三個方面的問題：

(1)AOD 與 PM2.5關(guān)系的穩(wěn)定性。這是通過AOD 反演 PM2.5的基本假設(shè)，大量研究結(jié)果表明，AOD 與 PM2.5存在一定的統(tǒng)計相關(guān)性，但是這個相關(guān)性在不同的區(qū)域有所不同，在同一區(qū)域的不同時間也可能不同。因此，針對特定區(qū)域、特定時間的 PM2.5反演，該關(guān)系的穩(wěn)定性起著至關(guān)重要的作用。

(2)誤差傳遞過程。通過建立各種精細(xì)的物理模型，提高 AOD 反演精度，從而能夠更精確地反演 PM2.5，但仍然存在一個誤差傳遞的過程。反演 AOD 是有誤差的，從帶誤差的 AOD 反演 PM2.5仍然會存在誤差，因此，誤差傳遞過程可能會導(dǎo)致某些區(qū)域的 PM2.5反演精度偏低?，F(xiàn)在也有大量學(xué)者通過日平均、月平均、季平均、年平均等尺度研究 AOD 與 PM2.5的關(guān)系，這在一定程度上能夠抵消誤差傳遞過程帶來的偏差。但在某一時刻尺度上，這個誤差傳遞過程造成的影響可能更大。

(3)模型的適用性。通過引入各種訂正、加入更多的輔助數(shù)據(jù)、結(jié)合數(shù)值模式等能夠提高AOD 反演的精度。但是，加入更多的因子，就意味著引入了更多的不確定性，對模型的適用性提出了更嚴(yán)格的要求。因此，同樣的方法，換一個研究區(qū)，效果可能會變得很差。

針對上述問題，本文提出了一種基于隨機森林機器學(xué)習(xí)法與 MODIS 影像相結(jié)合的 PM2.5遙感反演方法。從 MODIS 遙感數(shù)據(jù)出發(fā)，通過機器學(xué)習(xí)的手段，直接建立遙感影像本身與實測PM2.5的關(guān)系，以避免誤差的傳遞。初步實驗結(jié)果表明，反演的結(jié)果與地面實測 PM2.5具有較好的相關(guān)性。

2 MODIS 的氣溶膠光學(xué)厚度反演算法

TERRA 和 AQUA 是美國地球觀測系統(tǒng)計劃中的兩顆重要衛(wèi)星，它們搭載的 MODIS 掃描寬度為 2330km，具有 36個光譜波段，波長范圍0.14～14μm，空間分辨率 0.25～1km。MODIS以其高時空分辨率、多通道、覆蓋范圍廣等優(yōu)點被廣泛應(yīng)用于氣溶膠光學(xué)厚度的反演。

MODIS AOD 算法自問世以來經(jīng)過多次改進(jìn)，現(xiàn)已更新到 C6版本。在 2008年發(fā)布的C5版本中，暗目標(biāo)法(Dark Target Algorithm，DT)，又稱暗像元法，只用于暗目標(biāo)地區(qū)的反演；深藍(lán)算法(Deep Blue Algorithm，DB)只用于反演亮目標(biāo)區(qū)域。DT 與 DB 反演結(jié)果不做融合，只提供分辨率 10km 的氣溶膠產(chǎn)品。在2012年發(fā)布的 C6版本中，DT 與 DB 反演結(jié)果進(jìn)行融合，并且基于 DT 算法反演的 AOD 產(chǎn)品分辨率可達(dá) 3km。

反演氣溶膠光學(xué)厚度的基本原理[1]是，假定觀測表面是均勻的朗伯面，建立大氣頂層輻射亮度值與表面反射率關(guān)系，在不考慮大氣吸收情況下，衛(wèi)星接收的輻射值為：

其中，為衛(wèi)星觀測表面的反射率；為大氣反射率；分別為衛(wèi)星天頂角和太陽天頂角；為衛(wèi)星與太陽之間的相對方位角；S 為大氣下界的半球反射率；T 為大氣透過率。

當(dāng)?shù)乇矸瓷渎屎苄r，衛(wèi)星觀測反射率主要取決于大氣分子和氣溶膠散射發(fā)生的反射率。在反演過程中，先假設(shè)不同的氣溶膠模型和不同的觀測幾何狀況，再計算氣溶膠光學(xué)厚度與大氣下界的半球反射率、大氣反射率及之間的對應(yīng)關(guān)系。據(jù)此建立查找表，通過動態(tài)氣溶膠模式的輸入來查算氣溶膠光學(xué)厚度。

陸地上的植被、濕土壤和水體在可見光波段發(fā)射率都很低，在衛(wèi)星圖像上被稱為暗像元。在無云的暗像元上空區(qū)域，衛(wèi)星觀測反射率與氣溶膠光學(xué)厚度之間呈正比例關(guān)系，利用這種關(guān)系反演 AOD 的算法被稱為暗像元法。暗像元法根據(jù)遙感圖像的歸一化植被指數(shù)(Normalized Difference Vegetation Index，NDVI)值和短波紅外通道(2.13μm 和 3.8μm)觀測值進(jìn)行暗像元識別，再依據(jù)一定的關(guān)系假定這些暗像元在可見光紅藍(lán)通道的地表反射率，之后基于表觀反射率的大氣貢獻(xiàn)項和大氣輻射傳輸模型建立氣溶膠查找表，以此來反演氣溶膠光學(xué)厚度。該方法基于表觀反射率的大氣貢獻(xiàn)項及利用衛(wèi)星觀測的路徑輻射反演氣溶膠光學(xué)厚度。利用暗像元法反演AOD 時，對茂密的綠色植被、濕土壤和水體等低地表反射率區(qū)域反演效果明顯，但對中高緯度的冬季和干旱地區(qū)等高反射率區(qū)域，不能采用暗像元法反演 AOD。因為當(dāng)?shù)乇矸瓷渎试龃髸r，傳感器接收到的輻射值與氣溶膠厚度的正比例值將減小，甚至隨著反射率的增加，輻射值與氣溶膠厚度之間不再存在比例關(guān)系。

在中高緯度的冬季、裸地和沙漠等區(qū)域，晴天無云時地表反射率很高，但藍(lán)光波段對高亮地表具有低反射率的特征。深藍(lán)算法原理就是利用藍(lán)光波段這一特性，構(gòu)建高亮地表反射率數(shù)據(jù)庫，再通過查找表來構(gòu)建與最優(yōu)衛(wèi)星觀測輻射值的對應(yīng)關(guān)系，以此確定氣溶膠光學(xué)厚度。由于該算法只針對反射率高的地表物體，對海洋等低反射率地區(qū)不能達(dá)到很好的反演效果，因此開發(fā)了融合 DT/DB 算法的融合 AOD 產(chǎn)品。

由于暗像元法不能反演高反射率區(qū)域的氣溶膠厚度，而深藍(lán)算法無法反演海洋區(qū)域的AOD，鑒于此種弊端，C6版本將二者融合得到了 DT/DB 融合算法。該算法的核心思想是：對于海洋區(qū)域，選用暗目標(biāo)法反演 AOD；對于陸地上的高亮區(qū)域用深藍(lán)算法反演，陸地上的暗地表則用暗目標(biāo)法反演。陸地上用 NDVI 值來劃分暗地表和高亮區(qū)域。

3 本文的方法

3.1 隨機森林

隨機森林(Random Forest，RF)是并行式集成學(xué)習(xí)法的一個擴(kuò)展變體。RF 在以決策樹為基學(xué)習(xí)器構(gòu)建并行式集成學(xué)習(xí)法的基礎(chǔ)上，進(jìn)一步在決策樹的訓(xùn)練過程中引入隨機屬性選擇。它通過自助法重采樣技術(shù)，從原始訓(xùn)練樣本集 N 中有放回地重復(fù)隨機抽取 k 個樣本生成新的訓(xùn)練樣本集合，然后根據(jù)自助樣本集生成 k 個分類樹組成隨機森林，新數(shù)據(jù)的分類結(jié)果按分類樹投票多少形成的分?jǐn)?shù)而定。其實質(zhì)是對決策樹算法的一種改進(jìn)，將多個決策樹合并在一起，每棵樹的建立依賴于一個獨立抽取的樣品，森林中的每棵樹具有相同的分布，分類誤差取決于每一棵樹的分類能力和它們之間的相關(guān)性。特征選擇采用隨機的方法去分裂每一個節(jié)點，然后比較不同情況下產(chǎn)生的誤差。能夠檢測到的內(nèi)在估計誤差、分類能力和相關(guān)性決定選擇特征的數(shù)目。單棵樹的分類能力可能很小，但在隨機產(chǎn)生大量的決策樹后，一個測試樣品可以通過每一棵樹的分類結(jié)果經(jīng)統(tǒng)計后選擇最可能的分類。具體來說，傳統(tǒng)決策樹在選擇劃分屬性時是在當(dāng)前結(jié)點的屬性集合(假定有 d 個屬性)中選擇一個最優(yōu)屬性。而在 RF中，對基決策樹的每個結(jié)點，先從該結(jié)點的屬性集合中隨機選擇一個包含 k 個屬性的子集，然后再從這個子集中選擇一個最優(yōu)屬性用于劃分。這里的參數(shù) k 控制了隨機性的引入程度，一般取值為。

隨機森林是以 K 個決策樹為基本分類器，進(jìn)行集成學(xué)習(xí)后得到一個組合分類器。當(dāng)輸入待分類樣本時，隨機森林輸出的分類結(jié)果由每個決策樹的分類結(jié)果簡單投票決定。這里的 θk(k＝1,2,…,K)是一個隨機變量序列，它是由隨機森林的兩大隨機化思想決定的：

(1)Bagging 思想。從原樣本集 X 中有放回地隨機抽取 K 個與原樣本集同樣大小的訓(xùn)練樣本集(每次約 37% 的樣本未被抽中)，每個訓(xùn)練樣本集構(gòu)造一個對應(yīng)的決策樹。

(2)特征子空間思想。在對決策樹每個節(jié)點進(jìn)行分裂時，從全部屬性中等概率隨機抽取一個屬性子集(通常取 log2M＋1個屬性，M 為特征總數(shù))，再從該子集中選擇一個最優(yōu)屬性來分裂節(jié)點。

由于構(gòu)建每個決策樹時，隨機抽取訓(xùn)練樣本集和屬性子集的過程都是獨立的，且總體都是一樣的，因此是一個獨立同分布的隨機變量序列。訓(xùn)練隨機森林的過程就是訓(xùn)練各個決策樹的過程，由于各個決策樹的訓(xùn)練是相互獨立的，因此隨機森林的訓(xùn)練可以通過并行處理來實現(xiàn)，這將大大提高生成模型的效率。將以同樣的方式訓(xùn)練得到 K 個決策樹組合起來，就可以得到一個隨機森林。當(dāng)輸入待分類的樣本時，隨機森林輸出的分類結(jié)果由每個決策樹的輸出結(jié)果進(jìn)行投票決定。

3.2 本文方法

如前所述，本文試圖越過反演 AOD 的過程，基于隨機森林的機器學(xué)習(xí)方法，直接建立PM2.5與 MODIS 影像本身的關(guān)系。具體而言，分為以下幾個步驟。

3.2.1 時空匹配預(yù)處理

本文使用的數(shù)據(jù)是 MOD021KM，空間分辨率為 1km，包含 16個發(fā)射率波段、22個輻射率波段和 22個反射率波段。AOD 是 MODIS 提供的 3km 產(chǎn)品，地面實測 PM2.5數(shù)據(jù)使用了 102個站點的每小時觀測數(shù)據(jù)。

在時間匹配方面，由于 MODIS Terra 衛(wèi)星過境時間是上午十點半，因此選取過境當(dāng)天上午十點與十一點的 PM2.5監(jiān)測數(shù)據(jù)，并計算其平均值，作為衛(wèi)星過境時的地面觀測值。

在空間匹配方面，由于 AOD 數(shù)據(jù)空間分辨率為 3km，MODIS 數(shù)據(jù)的空間分辨率為 1km。因此，通過地面監(jiān)測站點的經(jīng)緯度實現(xiàn)監(jiān)測站點與影像數(shù)據(jù)的空間匹配。同時，為了直觀顯示PM2.5的真實空間分布，將所有站點的 PM2.5監(jiān)測值通過克里金插值法，插值成空間分辨率為 1km的數(shù)據(jù)。另外，在研究區(qū)內(nèi)，受云層及其他因素的影響，AOD 數(shù)據(jù)經(jīng)常出現(xiàn)大量數(shù)據(jù)缺失，因此，采用克里金插值法將缺失的數(shù)據(jù)進(jìn)行插值。

3.2.2 樣本選擇

本文研究區(qū)內(nèi)有 102個 PM2.5地面監(jiān)測站點，大致按照 7∶3的比例，隨機將 70個站點用于訓(xùn)練，32個站點用于測試。在生成訓(xùn)練樣本時需要考慮云層的影響。將云產(chǎn)品疊加到 MODIS數(shù)據(jù)上，如果站點位置有云，則該站點對應(yīng)位置的像素不作為訓(xùn)練樣本。對于第 i 個站點，如果其對應(yīng)的 MODIS 影像像素不為云，則該像素為有效像素，其對應(yīng)的訓(xùn)練樣本格式如公式(2)所示。

其中，x1～x16為 16個波段的發(fā)射率值；x17～x38為 22個波段的輻射率值；x39～x60為 22個波段的反射率值；yi為當(dāng)前站點對應(yīng)的 PM2.5實測值。

為了提高模型的預(yù)測能力，本文對訓(xùn)練樣本做了一定的增強處理，即除了選取當(dāng)前站點對應(yīng)的 MODIS 影像像素外，同時也選取了該像素5×5鄰域內(nèi)的所有像素，連同這些像素對應(yīng)于插值后的 PM2.5實測數(shù)據(jù)，一起構(gòu)成新的訓(xùn)練樣本。這樣做的理由是，根據(jù)地理學(xué)第一定律，對于插值后的 PM2.5實測數(shù)據(jù)，本文認(rèn)為站點附近5×5鄰域內(nèi)的插值數(shù)據(jù)具有很大的可信度，可以認(rèn)為是真實值。通過這種方式，在沒有云的情況下，一個站點最多可以生成 25個訓(xùn)練樣本。

對于測試樣本，則只選取當(dāng)前站點對應(yīng)的MODIS 影像像素的值。因此，在沒有云的情況下，最多可以有 32個測試樣本。

3.2.3 模型訓(xùn)練

一般而言，訓(xùn)練樣本的分布越均勻，訓(xùn)練的模型越具有代表性。由于 102個站點是按照大約7∶3的比例隨機分配，因此，為了達(dá)到訓(xùn)練樣本分布的均勻性，本文將此隨機分配過程重復(fù) 150次，得到 150組訓(xùn)練樣本及其配套的測試樣本。將每一組樣本輸入到隨機森林算法中，得到 150個訓(xùn)練模型。然后將每組樣本中的測試樣本輸入到對應(yīng)的模型中，得到對應(yīng)的預(yù)測值。選擇預(yù)測值表現(xiàn)指標(biāo)最優(yōu)的模型作為最終的模型，其對應(yīng)的訓(xùn)練樣本和測試樣本作為最終選出的樣本。

表現(xiàn)指標(biāo)以均方根誤差(Root Mean Square Error，RMSE)最小來判定，即：

其中，j 為訓(xùn)練的組數(shù)；Nj為該組的有效測試樣本數(shù)；分別為預(yù)測值和實際觀測值。選擇 RMSE 最小的一個模型為最優(yōu)模型。

3.2.4 模型測試

將上一步中選出的最優(yōu)模型用于整幅MODIS 影像的無云區(qū)域，對于影像中有云區(qū)域PM2.5的值以 0代替，從而得到整幅 MODIS 影像的 PM2.5反演結(jié)果。

4 結(jié)果與分析

4.1 研究區(qū)

廣東省地處中國大陸最南部，東鄰福建，北接江西、湖南，西連廣西，南臨南海，珠江口東西兩側(cè)分別與香港、澳門特別行政區(qū)接壤，西南部雷州半島隔瓊州海峽與海南省相望。全境位于北緯 20°13′～25°31′ 和東經(jīng) 109°39′～117°19′，東西跨度約 800km，南北跨度約 600km。全省陸地面積為 179800km2。廣東省屬于東亞季風(fēng)區(qū)，從北向南分別為中亞熱帶、南亞熱帶和熱帶氣候，是中國光、熱和水資源最豐富的地區(qū)之一。以廣州為核心的珠三角地區(qū)是中國城市化進(jìn)程最快的區(qū)域之一，伴隨而來的大氣污染問題也比較突出。本文的研究區(qū)域如圖1所示，圖中三角形標(biāo)示點為 102個環(huán)境監(jiān)測站，逐小時發(fā)布PM2.5監(jiān)測數(shù)據(jù)。

為了驗證本文方法的有效性，采用了MODIS 的 L2級 1km 數(shù)據(jù)(MOD021KM)對PM2.5進(jìn)行反演，數(shù)據(jù)時間分別為 2015.04.15、2015.04.17、2015.08.08、2015.08.26、2015.10.15、2015.12.20、2016.02.06、2016.02.09、2016.03.20，時間跨越 2個年份，包含了 4個季節(jié)。數(shù)據(jù)來源于美國國家航天宇航局(https://ladsweb.modaps.eosdis.nasa.gov/)。該數(shù)據(jù)包含 16個波段的發(fā)射率數(shù)據(jù)、22個波段的反射率數(shù)據(jù)和 22個波段的輻射率數(shù)據(jù)。作為對比，同時采用 MODIS 產(chǎn)品中分辨率最高的 3km 氣溶膠產(chǎn)品(AOD)進(jìn)行試驗分析(http://modis-atmos.gsfc.nasa.gov/products.html)，該產(chǎn)品采用最新C6版本中的 DT 與 DB 融合算法。本文使用的隨機森林算法通過 Weka[18]來實現(xiàn)，網(wǎng)址為 https://www.cs.waikato.ac.nz/ml/weka/index.html。

圖1 研究區(qū)域Fig.1Study area

受廣東省氣候環(huán)境的影響，MODIS 數(shù)據(jù)經(jīng)常被大量云層覆蓋，導(dǎo)致 AOD 產(chǎn)品上經(jīng)常出現(xiàn)大面積數(shù)據(jù)缺失。因此，在使用時，需要使用插值算法彌補這些數(shù)據(jù)缺失。本文采用克里金插值方法。也有很多研究者自行反演 AOD，但是其精度往往取決于引入的更多輔助數(shù)據(jù)和特別操作。因此，為了消除其他因素的影響，本文僅僅使用 MODIS 發(fā)布的最高分辨率的 AOD 產(chǎn)品，通過最經(jīng)典的線性回歸方法反演 PM2.5。線性回歸模型基本形式如公式(4)所示。

決定系數(shù)是指在表征因變數(shù)的總平方和中，由自變數(shù)引起的平方和所占的比例，稱為R平方，記為R2。由于R2＜R 可以防止對相關(guān)系數(shù)所表示的相關(guān)做夸張的解釋，因此決定系數(shù)的大小決定了相關(guān)的密切程度。當(dāng)R2越接近1時，表示相關(guān)方程式參考價值越高；相反，越接近 0時，表示參考價值越低。表達(dá)公式為：

其中，y 為待擬合數(shù)值；為其均值；為其擬合值。

4.2 2015年8月8日反演結(jié)果對比

圖2給出了 2015年8月8日的 MOD021KM數(shù)據(jù)。為顯示方便，發(fā)射率和輻射率采用 123波段合成，反射率采用 456波段合成。當(dāng)日的數(shù)據(jù)中有部分云層，采用 MODIS 的云檢測產(chǎn)品構(gòu)建掩膜，如圖2(f)所示的黑色區(qū)域即為云區(qū)。當(dāng)日的 AOD 數(shù)據(jù)存在大量缺失，如圖2(d)中的黑色區(qū)域。將此數(shù)據(jù)進(jìn)行克里金插值，并用假彩色顯示，具體如圖2(e)所示，表現(xiàn)出明顯的區(qū)塊效應(yīng)。PM2.5的地面觀測值是點狀數(shù)據(jù)，本文利用克里金插值將點狀數(shù)據(jù)插值為面狀數(shù)據(jù)，如圖2(f)、(g)所示，其中圖2(f)加了云掩膜。圖2(h)是基于克里金插值后的 AOD 數(shù)據(jù)經(jīng)過線性回歸反演得到的 PM2.5。圖2(i)是本文方法得到的 PM2.5結(jié)果，其中顏色越紅，表示 PM2.5濃度越大；顏色越藍(lán)，表示 PM2.5濃度越低。

圖2 2015年8月8日實驗結(jié)果Fig.2Experimental results on 2015.8.8

從圖2(f)、(g)所示的地面觀測值可以看出，中間區(qū)域的 PM2.5濃度很高，東北和西南兩個區(qū)域的濃度較低。AOD 反演的結(jié)果與地面觀測結(jié)果差異很大，這是由 AOD 數(shù)據(jù)缺失導(dǎo)致的。而本文方法在整體趨勢上與地面觀測結(jié)果非常一致，表現(xiàn)出中間高、東北和西南低的趨勢。受云層影響，32個驗證站中只有 26個站有數(shù)據(jù)，因此給出了 26個地面觀測站的統(tǒng)計結(jié)果。

由圖3(a)可以看出，本文方法在各個觀測站上的預(yù)測值都能與實際觀測值有較好的匹配，而AOD 方法匹配度較差。由散點圖和線性擬合結(jié)果(圖3(b)、(c))可以看出，本文方法的 R2達(dá)到0.97，RMSE 小于 2，表現(xiàn)出了極強的相關(guān)性；而 AOD 方法表現(xiàn)非常差，這也說明 AOD 的數(shù)據(jù)缺失對 PM2.5的反演有極大的負(fù)面影響。

4.3 2015年4月17日反演結(jié)果對比

另外選擇云量更少的一天的數(shù)據(jù)進(jìn)行實驗，結(jié)果如圖4所示。圖中少量的黑色區(qū)域即為MODIS 云檢測產(chǎn)品提供的云掩膜。AOD 產(chǎn)品的數(shù)據(jù)缺失程度較上一個實驗有了明顯改善，采用克里金插值后沒有表現(xiàn)出明顯的區(qū)塊效應(yīng)。AOD和本文方法反演的結(jié)果如圖4(h)、(i)所示。

圖3 2015年8月8日實驗結(jié)果Fig.3Experimental results on 2015.8.8

圖4 2015年4月17日實驗結(jié)果Fig.4Experimental results on 2015.4.17

從圖4(f)、(g)所示的地面觀測值插值結(jié)果可以看出，中部 PM2.5濃度較大，東北角區(qū)域次之，西南角區(qū)域最小。AOD 反演結(jié)果在中部區(qū)域跟地面觀測值較一致，但是西南角區(qū)域明顯偏大。而本文方法依然與地面觀測結(jié)果表現(xiàn)出明顯的一致性。由于云量較少，32個驗證站點中有31個屬于有效站點。AOD 和本文方法反演的結(jié)果與地面觀測結(jié)果如圖5(a)所示，可以看出，大部分站點上本文方法與地面觀測結(jié)果吻合度非常好，而 AOD 則有明顯的偏差。

由散點圖和線性擬合統(tǒng)計結(jié)果(圖5(b)、(c))可以看出，本文方法的 R2遠(yuǎn)比 AOD 反演的高，RMSE 更低，表現(xiàn)出更強的線性關(guān)系；而 AOD 反演結(jié)果則表現(xiàn)得比較離散。由于云量比上一個實驗少，所以 AOD 反演的 R2有了明顯提高。

用同樣方法對其他幾個日期的數(shù)據(jù)進(jìn)行了實驗，計算 R2和 RMSE 的平均值，結(jié)果如表 1所示。

由表 1可以看出，本文方法(基于 RF 反演)的 PM2.5均值比基于 AOD 法要高，而 RMSE 更低。從 RMSE 來看，本文方法也具有比較明顯的優(yōu)勢。R2和 RMSE 波動的主要因素是云量的影響，以及 AOD 數(shù)據(jù)本身的缺失問題。因為在本研究區(qū)，AOD 數(shù)據(jù)缺失現(xiàn)象有時比較嚴(yán)重，通過克里金插值補齊的數(shù)值并不能完全反映真實的AOD 空間分布。

圖5 2015年4月17日實驗結(jié)果Fig.5Experimental results on 2015.4.17

表1 評價指標(biāo)的平均值Table 1The average values of assessment indices

5 結(jié) 論

MODIS AOD 算法的不斷改進(jìn)，目的是得到結(jié)果更為精確的 AOD 產(chǎn)品。但在反演 AOD 過程中不可避免地會有誤差的存在，因此現(xiàn)在常用的 AOD 反演 PM2.5方法無法避免誤差。本文結(jié)合隨機森林的機器學(xué)習(xí)算法，從遙感影像本身數(shù)據(jù)出發(fā)，直接建立遙感影像與實測 PM2.5的關(guān)系，從而避免了誤差傳遞。選取了 MODIS 數(shù)據(jù)分辨率 3km 的 AOD 產(chǎn)品和廣東省 102個環(huán)境監(jiān)測站點的 PM2.5數(shù)據(jù)進(jìn)行試驗。試驗結(jié)果表明，本方法能夠取得更好的 PM2.5反演效果，同時將 PM2.5反演的空間分辨率提高到 1km。下一步研究將擴(kuò)大研究區(qū)域，采用更多的數(shù)據(jù)，進(jìn)一步提高算法的可用性。另外，采用其他更好的機器學(xué)習(xí)方法來確定反演模型也是今后研究的重點。

參 考 文 獻(xiàn)

[1]Lin CQ,Li Y,Yuan ZB,et al.Using satellite remote sensing data to estimate the high-resolution distribution of ground-level PM2.5[J].Remote Sensing of Environment,2015,156: 117-128.

[2]Gong W,Huang S,Zhang TH,et al.Impact and suggestion of column-to-surface vertical correction scheme on the relationship between satellite AOD and ground-level PM2.5in China [J].Remote Sensing,2017,9(10): 1038.

[3]Chen ZY,Chen DL,Zhuang Y,et al.Examining the infl uence of crop residue burning on local PM2.5concentrations in Heilongjiang province using ground observation and remote sensing data [J].Remote Sensing,2017,9(10): 971.

[4]陳輝,厲青,張玉環(huán),等.基于地理加權(quán)模型的我國冬季 PM2.5遙感估算方法研究 [J].環(huán)境科學(xué)學(xué)報,2016,36(6): 2142-2151.

[5]王子峰.衛(wèi)星遙感估算近地面顆粒物濃度的算法研究 [D].北京: 中國科學(xué)院遙感應(yīng)用研究所,2010.

[6]王中挺,辛金元,賈松林,等.利用暗目標(biāo)法從高分一號衛(wèi)星 16m 相機數(shù)據(jù)反演氣溶膠光學(xué)厚度[J].遙感學(xué)報,2015,19(3): 530-538.

[7]王靜,楊復(fù)沫,王鼎益,等.北京市 MODIS 氣溶膠光學(xué)厚度和 PM2.5質(zhì)量濃度的特征及其相關(guān)性[J].中國科學(xué)院研究生院學(xué)報,2010,27(1): 10-16.

[8]Song WZ,Jia HF,Huang JF,et al.A satellitebased geographically weighted regression model for regional PM2.5,estimation over the pearl river delta region in China [J].Remote Sensing of Environment,2014,154: 1-7.

[9]張?zhí)炱?基于 MODIS 影像的中國地區(qū)氣溶膠產(chǎn)品驗證與 PM2.5反演 [D].成都: 西南交通大學(xué),2017.

[10]Chu DA,Tsai TC,Chen JP,et al.Interpreting aerosol lider profiles to better estimate surface PM2.5,for columnar AOD measurements [J].Atmospheric Environment,2013,79(11): 172-187.

[11]Griggs M.Measurements of atmospheric aerosol optical thickness over water using ERTS-1data [J].Journal of the Air Pollution Control Association,1975,25(6): 622-626.

[12]Levy RC,Kozak GM,Wadsworth CB,et al.Towards a long-term global aerosol optical depth record: applying a consistent adrosol retrieval algorithm to MODIS and VIIRS-observed reflectance [J].Atmospheric Measurement Techniques,2015,10(8): 4083-4110.

[13]Holben B,Vermote E,Kaufman YJ,et al.Aerosol retrieval over land from AVHRR data application for atmospheric correction [J].IEEE Transactions on Geoscience and Remote Sensing,1992,30(2):212-222.

[14]Isakov VY,Feind RE,Vasilyev OB,et al.Retrieval of aerosol spectral optical thic-kness from AVIRIS data [J].International Journal of Remote Sensing,1996,17(11): 2165-2184.

[15]Shi YS,Matsunaga T.Long-term trends and sptial patterns of satellite-retrieved PM2.5concentrations in South and Southeast Asia from 1999to 2014[J].Science of the Total Environment,2018,615: 177-186.

[16]Jung CR,Hwang BF,Chen WT.Incorporating longterm satellite-based aerosol optical depth,localized land use data,and meteorological variables to estimate ground-level PM2.5concentrations in Taiwan from 2005to 2015[J].Environmental Pollution,2017,doi: 10.1016/j.envpol.2017.11.016.

[17]Mao X,Shen T,Feng X.Prediction of hourly ground-level PM2.5concentrations 3days in advance using neural networks with satellite data in eastern China [J].Atmospheric Pollution Research,2017,8(6): 1005-1015.

[18]The WEKA Workbench.Data mining: practical machine learning tools and techniques [EB/OL].[2018-3-12].https://www.cs.waikato.ac.nz/ml/weka/citing.html.