張登科 呂偉才 劉 歡 仲 臣

(1. 安徽理工大學 空間信息與測繪工程學院, 安徽 淮南 232001; 2. 河南理工大學 測繪與國土信息工程學院, 河南 焦作 454002)

0 引言

霧霾污染較為嚴重,所以關于霧霾的大量研究層出不窮,在這些研究中涉及霧霾識別的研究也比較多。一些學者認為氣溶膠是表征霧霾污染的一個重要特征,因此,他們從氣溶膠與霧霾天氣的內部聯(lián)系出發(fā),以中分辨率成像光譜儀(Moderate Resolution Imaging Spectroradio meter,MODIS)遙感圖像為數(shù)據(jù)源,利用當前在氣溶膠厚度反演方面最成熟的美國國家航空航天局(National Aeronautics and Space Administration,NASA)的V5.2算法對嚴重霧霾污染地區(qū)進行識別[1-2];一些學者利用地面氣象觀測霧和霾,并在已有綜合識別算法基礎上,提出利用測云激光雷達的后向散射信號信息,并設置合理的閾值,結合綜合識別算法來識別霧霾[3];一些學者提出了一種基于K均值聚類算法的霧霾天氣識別方法,該方法在傳統(tǒng)K均值聚類算法的基礎上,提取圖像飽和度分量作為識別霧霾天氣的特征,采用馬氏距離作為圖像間的相似性度量,將相鄰迭代次數(shù)內聚類中心的變化距離作為聚類準則函數(shù),對圖像進行聚類,實現(xiàn)霧霾天氣的識別[4-5]。而本文進行重霧霾區(qū)域識別算法的研究是著重分析京津冀等霧霾比較嚴重的區(qū)域的霧霾天氣特性,并通過分析地物光譜反射率借助于算法計算出一個能有效識別霧霾區(qū)域的新的霧霾指標。

本文中霧霾識別的原理主要是對2007年法國小衛(wèi)星PARASOL (Polarization and Anisotropy of Reflectances for Atmospheric Sciences Coupled with Observations from a LiDAR)獲取的數(shù)據(jù)進行處理和提取,為了研究重霧霾區(qū)域下的天氣特征,從數(shù)據(jù)中提取了與研究相關的光譜反射率,通過研究區(qū)域內的地物光譜反射率的值[6],并對該值作相應的處理計算,最后得出一個能夠識別非霧霾天氣與霧霾天氣的指標,即修正歸一化差異霾指數(shù)(Modify the Normalized Differential Haze Index，M-NDHI)。

1 數(shù)據(jù)與方法

1.1 研究區(qū)域

本次研究選取的區(qū)域位于華北中部平原,是中國霧霾影響最嚴重的四個地區(qū)之一——北京。由于在數(shù)據(jù)讀取以及提取過程中只需要研究北京區(qū)域,因此,需要在程序中設置北京區(qū)域的經(jīng)緯度,北京市界的地理坐標為:39°26′～41°03′N,115°25′～117°30′E。在中國的這些城市中,北京市是最容易形成霧霾的城市。造成這一狀況的原因主要有兩個:其特殊的地理位置;近年來其工業(yè)化和城市化的快速發(fā)展。從地理位置上看,北京三面環(huán)山,污染物難以及時擴散,又毗鄰東海,容易形成霧,造成污染物積累[7]。這就為霧霾的長期繁榮創(chuàng)造了一個最佳地理條件。其次,汽車尾氣和工業(yè)廢氣的迅速增加使得這個城市的霧霾狀況更加嚴重,大量污染物集聚到空中,促進了霧霾的形成,破壞了城市的環(huán)境質量[8]。因此,監(jiān)測空氣污染狀況的新方法的性能極為重要。

1.2 數(shù)據(jù)來源與處理方法

1.2.1數(shù)據(jù)選取

本次研究使用的數(shù)據(jù)來源于iCare數(shù)據(jù)和服務中心,主要分析了大量2007年以后的云、氣溶膠等的輻射交互作用數(shù)據(jù),通過提取數(shù)據(jù)得到各種植被對霧霾的光譜反射特性和光譜反射率數(shù)據(jù),將所得到的數(shù)據(jù)代入霧霾指數(shù)求算公式,從而將所得到數(shù)值與非霧霾天氣下所得到的數(shù)值閾值做比較,能夠作為評判一個地區(qū)是否有霧霾的依據(jù)。本文選取了2007年5、6月以及少量其他月份的數(shù)據(jù)進行分析,在地球觀測系統(tǒng)數(shù)據(jù)和信息系統(tǒng)世界觀(Earth Observing System Data and Information System,EOSDIS Worldview)網(wǎng)站上查找研究了這幾個月的霧霾天氣信息,并將霧霾天、非霧霾天分別記錄下來,然后對其進行數(shù)據(jù)提取和數(shù)據(jù)處理,從而得到計算公式所需要的反射率的值。

在數(shù)據(jù)提取和處理的過程中,主要提取計算了490、670、865 nm三個波段在不同天氣下的光譜反射率,并計算出了日光譜反射率的平均值。光譜反射率=I(非偏振光強)/cos(太陽天頂角),在程序中的計算公式為[9]

R1=I1/cosθ
R2=I2/cosθ
R3=I3/cosθ

(1)

式中，R1、R2、R3分別為490、670、875 nm波段光譜反射率;I1、I2、I3分別為490、670、875 nm波段非偏振光強;θ為太陽天頂角。

1.2.2數(shù)據(jù)處理

在建立表征地表類型特征的指標時,從光譜波段提取的元素值要高于得到的圖像中的其他特征值。比如,在歸一化植被指數(shù)(Normalized Difference Vegetation Index,NDVI)中,植被的值高于裸地、建筑物和水。同樣地,在M-NDHI指數(shù)中,霧霾的值比其他元素要高[10]。但值得注意的是,M-NDHI指數(shù)不同于歸一化水指數(shù)(Norma-lized Difference Water Index,NDWI)、歸一化建筑物指數(shù)(Normalized Difference Build-up Index,NDBI)、NDVI等其他指數(shù),因為霧霾是在地表形成的物質,所以衛(wèi)星圖像上的光譜值是霧霾加上地面元素的和。因此,為了有效地識別霧霾,必須抑制研究區(qū)域的下墊面因素。本文主要技術流程如圖1所示。

通過比較在三種天氣下三個波段的光譜反射率,可以得到以下三個結論:(1)在同一種天氣條件下,865波段(865 nm;近紅外波段)的光譜反射率值最大,670波段(670 nm;可見光紅光波段)的光譜反射率值最小,而490波段(490 nm;可見光青光區(qū)域)的光譜反射率值介于兩者之間;(2)同一波段在霧霾和非霧霾天氣下的光譜反射率大小也有較大的區(qū)別,在非霧霾天氣下光譜反射率的值較小,而有霧霾天氣下的光譜反射率的值較大;(3)在表1中除了統(tǒng)計不同天氣下、不同波段的光譜反射率大小,還計算出一組差異分析數(shù)據(jù),該數(shù)據(jù)是由霧霾天氣下的光譜反射率值與非霧霾天氣下的光譜反射率值求差值所得到的,在所求的三個數(shù)值中,490 nm波段的數(shù)值差異較大,865 nm波段的數(shù)值差異較小,670 nm波段的數(shù)值差異則介于兩者之間。

通過比較非霧霾天和霧霾天的數(shù)據(jù)處理結果,發(fā)現(xiàn)490 nm波段(可見光青光區(qū)域)對霧霾最為敏感。在三個波段間490 nm波段和865 nm波段(近紅外波段)之間的光譜變化最大(表1),但是由于865 nm波段無論是在霧霾天還是非霧霾天,光譜反射率的值相差不大,而且經(jīng)過計算670 nm波段對于有云的天氣下,光譜反射率的值較大,具有很好的效果。因此,選擇490 nm波段和670 nm波段來構建M-NDHI,即(490 nm波段-670 nm波段)/(490 nm波段+670 nm波段),則分子(490 nm波段-670 nm波段)在云層和云層附近較其他位置大得多;分母同樣是(490 nm波段+670 nm波段),如式(2)。因此,這種推導的M-NDHI對于較厚的霧霾具有更大的值。

(2)

式中，M為修正歸一化差異霾指數(shù)(M-NDHI)值;R1、R2分別為490、670 nm波段光譜反射率。

從圖2中很容易看出兩種不同天氣下的北京市,在M-NDHI計算公式中涉及的數(shù)據(jù)主要是霧霾天和非霧霾天的數(shù)據(jù),不考慮有云天氣下的北京市光譜反射率,因為在有云的天氣下,光譜反射率的值呈現(xiàn)不規(guī)律的變化,容易對研究造成誤導,因此,在數(shù)據(jù)的提取、分析以及計算的過程中都是無云天氣下的數(shù)據(jù),而無云天氣又分為霧霾天和非霧霾天。

圖2 分別為霧霾天和非霾天的北京云圖[審圖號為GS(2019)3266號]

1.2.3閾值的取定

為了確定閾值的取值范圍,本文對提取的幾個月內北京市霧霾天與非霧霾天的數(shù)據(jù)進行分析和計算,借助修正歸一化霾指數(shù)分別計算出來各天所對應的指數(shù)值,為了能夠更好地看出霧霾天與非霧霾天M-NDHI值的變化情況,以便于更好地確定閾值的范圍,將霧霾天與非霧霾天求取的M-NDHI值按照時間的順序做成折線圖,如圖3所示。

圖3 北京區(qū)域霧霾天氣下和非霧霾天氣下M-NDHI值的變化折線圖

圖3分別為非霧霾和霧霾天氣下北京市光譜反射率的值,通過分析、統(tǒng)計和計算2007年北京市的霧霾天和非霧霾天,選取了天氣較為明顯的一些時間作為參考,根據(jù)推導得出的修正歸一化霾指數(shù)(M-NDHI)計算得出了某些天的指數(shù)值,并將這些指數(shù)進行排序,最終得出了修正歸一化霾指數(shù)在霧霾天和非霧霾天的指數(shù)范圍。根據(jù)計算出的指數(shù)范圍,最初確定了霧霾天氣下修正歸一化霾指數(shù)的取值范圍為0.082～0.220,非霧霾天氣下的修正歸一化霾指數(shù)的取值范圍為0.020～0.082(不包含0.082)。因為研究數(shù)據(jù)的局限性沒有準確地得到指數(shù)的閾值,只是通過計算得出的數(shù)值推測的一個指數(shù)范圍,因此,具有一定的模糊性。

1.2.4M-NDHI閾值的驗證

為了使推導得到的指數(shù)公式更有說服性,也為了驗證指數(shù)閾值是否能夠作為評判霧霾天與非霧霾天的標準,于是又對天津市做了相似的研究,北京市與天津市都是位于京津冀霧霾區(qū)的兩個城市,因此對于天津市的研究可以作為驗證北京區(qū)域得到的霾指數(shù)和閾值的依據(jù)。

表1 三個波段在霧霾天和非霧霾天的光譜反射率變化

同樣的將提取并計算出的霧霾天與非霧霾天的光譜反射率數(shù)據(jù)分別代入M-NDHI計算公式中,為了能夠更好地看出霧霾天與非霧霾天M-NDHI值的變化情況,以便于確定取定的閾值是否合理,遂將天津區(qū)域霧霾天與非霧霾天求取的M-NDHI值按照時間的順序做成了折線圖如圖4所示。

通過圖4可以發(fā)現(xiàn),天津市霧霾天與非霧霾天所得到的修正歸一化霾指數(shù)也在所得到的閾值之內。由此可以得到,在與北京市地貌和地類相差不大的區(qū)域可以利用推導得到的M-NDHI來判斷該區(qū)域在相應的時間是霧霾天還是非霧霾天。

2 實驗結果與分析

本文提出了一個識別北京和天津霧霾區(qū)域的修正歸一化差異霾指數(shù),該指數(shù)能夠有效識別京津冀霧霾區(qū)的霧霾狀況。為了使研究成果更加真實可靠,選取了2007年6月4日的北京,并截取了當天的霧霾圖。從圖5中可以看出,當天的北京市有一半處于霧霾區(qū),而另一半?yún)^(qū)域并沒有發(fā)生霧霾。

圖5 6月4日北京市霧霾云圖[審圖號為GS(2019)3266號]

本次研究提取了北京當天的490波段數(shù)據(jù)、670波段數(shù)據(jù)以及該區(qū)域經(jīng)緯度的數(shù)據(jù),并將提取的數(shù)據(jù)代入式(2)中,求出了霧霾區(qū)域內相應位置的M-NDHI離散值,如表2所示, 其中M-NDHI為0時,表示該點數(shù)據(jù)為空。

通過對照圖5和表2中的數(shù)據(jù),能夠從兩方面區(qū)分霧霾區(qū)與非霧霾區(qū)。首先在圖上能夠直觀地看出霧霾區(qū)與非霧霾區(qū)的分界線,從而判斷出一個區(qū)域是否發(fā)生霧霾,其次,可以通過觀測表中的數(shù)據(jù)發(fā)現(xiàn),在M-NDHI閾值的附近有明顯的數(shù)據(jù)差,通過與閾值的比較也能容易地區(qū)分霧霾區(qū)域。

在表2中可以看到左下角的數(shù)據(jù)略微偏大,因為這是有云的情況下,當既有云又有霾的情況下,光譜反射率的變化比較大,因此忽略了有云的天氣。

表2 北京區(qū)域內離散M-NDHI值

3 結束語

本文主要統(tǒng)計分析了北京與天津區(qū)域的霧霾狀況,因為這兩個城市是以京津冀為主的重大霾區(qū)的兩個重要城市,通過分析2007年北京市與天津市霧霾天與非霧霾天的地物光譜反射率的變化情況,從提取的三個波段中選取兩個波段來構建霧霾指數(shù)公式,并分別計算了數(shù)月的霧霾天與非霧霾天的M-NDHI值,從而選取了修正歸一化霧霾指數(shù)的閾值。提出的M-NDHI可為京津冀主要霧霾區(qū)域的識別提供一種新的算法和理論依據(jù)。

因為本次霧霾研究是在京津冀霧霾區(qū)域進行的研究,而該區(qū)域是一個地勢較為平坦并且地物類型相對較少的區(qū)域,而在一些高原、裸土、荒漠以及雪地等區(qū)域該指標與閾值可能不具有適應性。在后續(xù)的霧霾研究工作中,需要尋找一個適應范圍廣,能夠完美表征水域、平原、高山、裸土、雪地以及荒漠等各種地類的霧霾指數(shù),尋找更加精簡的霧霾識別算法。