孫灝，劉偉漢，王艷梅，周偉，蔡創(chuàng)創(chuàng)

(中國礦業(yè)大學(xué)(北京) 地球科學(xué)與測繪工程學(xué)院，北京 100083)

0 引言

京津冀地區(qū)是世界上人口密度最大的地區(qū)之一。改革開放以來，京津冀及周邊地區(qū)經(jīng)濟(jì)發(fā)展迅猛，然而，經(jīng)濟(jì)快速增長伴隨著化學(xué)燃料的快速消耗，也導(dǎo)致了顆粒物濃度的逐年增加，頻繁引起灰霾天氣，影響著氣候環(huán)境與人們的生活質(zhì)量[1]。

灰霾的核心物質(zhì)是空氣中懸浮的灰塵顆粒，氣象學(xué)上稱為氣溶膠顆粒[2]。氣溶膠是懸浮于地球大氣中，沉降速度小，半徑在0.001～100 μm的液態(tài)及固態(tài)粒子，顯著影響著氣候環(huán)境。氣溶膠光學(xué)厚度(aerosol optical depth，AOD)是定量描述氣溶膠對電磁波衰減作用的物理量，是氣溶膠最重要的參數(shù)之一。遙感技術(shù)能夠彌補(bǔ)氣溶膠地基觀測站空間性和連續(xù)性的不足，實(shí)現(xiàn)大范圍實(shí)時性的氣溶膠監(jiān)測。遙感器觀測到的表觀反射率是受大氣和地表共同影響的結(jié)果，且不同的氣溶膠類型對表觀反射率的影響也不同。因此，氣溶膠光學(xué)厚度的衛(wèi)星遙感反演的2個核心問題分別是從衛(wèi)星傳感器接收到的信息中去除地表噪聲以及確定合適的氣溶膠模式。

傳統(tǒng)的暗像元反演算法利用在晴朗無云的濃密植被等低地表反射率上空，氣溶膠使反射率增大這一原理反演AOD，利用對氣溶膠不敏感的短波紅外通道表觀反射率獲得紅藍(lán)波段的地表反射率[3]。研究表明，暗像元算法僅在暗地表能獲得較好的反演結(jié)果[4]，冬季華北平原大部分的植被消失，不能滿足暗像元的假定，監(jiān)測精度大大降低，難以滿足業(yè)務(wù)化需求[5]。針對地表反射率較高的旱季和中高緯度地區(qū)冬季，Hsu等提出的基于地表反射率庫的深藍(lán)算法已得到成功應(yīng)用[6-7]。2013年，改進(jìn)的深藍(lán)算法借助地表反射率數(shù)據(jù)與歸一化植被指數(shù)NDVI，將反演范圍從干旱半干旱亮表面擴(kuò)展到了包括植被等暗地表的廣泛陸地區(qū)域，并進(jìn)一步討論和改進(jìn)了氣溶膠模式與云雪檢測算法[8]。MODIS在V5.2版大氣氣溶膠反演中加入了深藍(lán)算法，并在C6版本提供了改進(jìn)的深藍(lán)算法產(chǎn)品以及暗像元與深藍(lán)算法融合產(chǎn)品[9]，但與MODIS反射率數(shù)據(jù)(1 km)相比，分辨率仍較低，無法滿足區(qū)域氣溶膠的高精度監(jiān)測需求[10]。不同地區(qū)的氣溶膠類型、地表光譜特性各不相同，MODIS氣溶膠產(chǎn)品且作為全球范圍內(nèi)的統(tǒng)一算法產(chǎn)品，對特定研究區(qū)域的針對性較差。本文利用融合互補(bǔ)思路，集成深藍(lán)算法適應(yīng)范圍廣和暗像元算法精度高的優(yōu)勢[11-12]，結(jié)合華北平原的自然地理特征，開展了一種暗像元算法與深藍(lán)算法結(jié)合反演冬季京津冀地區(qū)AOD的方法，并進(jìn)行了連續(xù)的AOD監(jiān)測。依據(jù)監(jiān)測結(jié)果，結(jié)合輔助資料，分析了京津冀地區(qū)冬季的氣溶膠分布特征與其影響因素。

1 研究區(qū)域和數(shù)據(jù)

1.1 研究區(qū)概況

注：該圖基于國家測繪地理信息局標(biāo)準(zhǔn)地圖服務(wù)網(wǎng)站下載的審圖號為GS(2016)1666號的標(biāo)準(zhǔn)地圖制作，底圖無修改。圖1 研究區(qū)域

京津冀是中國“首都圈”，總面積約21.8×104km2，常駐人口約1.1億，其地理位置及氣溶膠地基站點(diǎn)分布情況如圖1所示。作為人口密度集中、工業(yè)發(fā)達(dá)的經(jīng)濟(jì)圈，京津冀地區(qū)地處華北平原，東臨渤海，北有燕山山脈，西有太行山脈，這種復(fù)雜的地理環(huán)境作用形成的海陸風(fēng)環(huán)流、山谷風(fēng)環(huán)流對區(qū)域大氣氣溶膠的擴(kuò)散有著較大影響，灰霾天氣頻繁發(fā)生。對京津冀氣溶膠進(jìn)行全方位、長時間的連續(xù)監(jiān)測是有必要的。

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

本研究采用的數(shù)據(jù)包括AQUA星MODIS L1B數(shù)據(jù)MYD021KM，該數(shù)據(jù)是對經(jīng)過初步處理的一級產(chǎn)品進(jìn)行定位和定標(biāo)后所生成的二級產(chǎn)品，以HDF格式存儲，空間分辨率1 km，時間分辨率1 d，對于太陽反射波段(對應(yīng)MODIS波段1～19和波段26)提供反射率和輻亮度2種類數(shù)據(jù)產(chǎn)品；MYD021KM對應(yīng)的地理定位數(shù)據(jù)MYD03；地表發(fā)射率庫的構(gòu)建選擇8日合成的MODIS地表反射率產(chǎn)品MOD09A1。驗(yàn)證數(shù)據(jù)選擇AERONET(aerosol robotic network)氣溶膠監(jiān)測網(wǎng)提供的北京、香河地基站AOD數(shù)據(jù)。此外，還使用了MODIS 10 km分辨率氣溶膠產(chǎn)品MOD04_L2以及地理空間數(shù)據(jù)云提供下載的SRTMDEM 90 m分辨率原始高程數(shù)據(jù)。MODIS和AERONET數(shù)據(jù)分別通過NASA和AERONET官網(wǎng)下載。1月通常是京津冀地區(qū)冬季的最冷月份，將研究的時間跨度擴(kuò)充到3年，總時間擴(kuò)充到6個月，時間范圍定為2016—2018年1—2月。

2 研究方法

2.1 反演原理

遙感反演AOD的基本原理是氣溶膠粒子對電磁波散射的吸收作用使遙感器接收的輻射能的性質(zhì)和強(qiáng)度發(fā)生變化，通過測量這種變化就可以反演氣溶膠的特性[13]。

假設(shè)大氣水平均一、地表朗伯體，遙感器觀測到的表觀反射率可表示為[5]：

(1)

式中：ρ0是大氣路徑輻射項(xiàng)等效反射率；μS=cosθS，μV=cosθV，θS和θV分別是太陽天頂角和觀測天頂角；φ和ρS分別是相對方位角和地表雙向反射率；S和T分別是為大氣下界半球反射率和大氣透過率，T(μS)T(μV)作為一個大氣參數(shù)考慮。因此，地表反射率和大氣氣溶膠模式的確定是AOD反演的2個重點(diǎn)。

通過輻射傳輸模型計(jì)算S、ρ0、T(μS)T(μV)3個大氣參數(shù)和AOD之間的對應(yīng)關(guān)系，建立包含各參數(shù)的查找表，通過查找表獲取AOD。

2.2 反演過程

1)構(gòu)建查找表。查找表通過6S(second simulation of the satellite signal in the solar spectrum)輻射傳輸模型構(gòu)建，使用IDL調(diào)用6S軟件，設(shè)定7參數(shù)進(jìn)行輻射傳輸計(jì)算，包括9個太陽天頂角，12個觀測天頂角，16個相對方位角。設(shè)定0、0.025、0.5、1.0、1.5和1.95共6個在波長0.55 μm處的AOD以供內(nèi)插。根據(jù)研究區(qū)域與監(jiān)測時間段，大氣模式選擇中緯度冬季大氣。針對氣溶膠模式的選擇問題，國際氣象與大氣物理協(xié)會(IAMAP)由氣溶膠基本組分含量的不同定義了標(biāo)準(zhǔn)輻射大氣(SRA)下的大陸型、城市型及海洋型3種基本氣溶膠類型(表1)[14-15]。其中大陸型和城市型氣溶膠分別表示了陸地自然條件下和嚴(yán)重人為排放條件下的氣溶膠模型。京津冀地區(qū)的實(shí)際氣溶膠類型可以理解為在自然狀態(tài)的大陸型氣溶膠背景下疊加了一定程度的人為排放，介于城市型和大陸型之間。在不考慮海洋性粒子的情況下使用正交試驗(yàn)法將沙塵粒子和水溶粒子占比各分為3個等級，共得到7組候選氣溶膠模型(煤煙粒子占比介于1%到22%之間)，分別使用暗像元算法計(jì)算2017年1月區(qū)間內(nèi)的氣溶膠光學(xué)厚度，通過匹配AERONET北京站實(shí)測數(shù)據(jù)，得出了誤差最小的京津冀地區(qū)冬季氣溶膠類型作為本研究使用的氣溶膠模型。其中，沙塵性、水溶性和煤煙性粒子分別占比44%、37%和19%。

表1 氣溶膠模型及粒子組分

2)構(gòu)建深藍(lán)算法地表反射率庫。針對在高地表反射率地區(qū)，紅藍(lán)波段和近紅外波段地表反射率無法構(gòu)成線性關(guān)系，通過構(gòu)建地表反射率庫實(shí)現(xiàn)地氣解耦。本研究以MODIS地表反射率產(chǎn)品MOD09 8天合成產(chǎn)品為基礎(chǔ)構(gòu)建地表反射率庫，按日期存儲在數(shù)據(jù)文件中。

3)數(shù)據(jù)預(yù)處理。使用MODIS的L1B產(chǎn)品讀取相關(guān)波段的表觀反射率及相應(yīng)的定標(biāo)系數(shù)，使用地理定位數(shù)據(jù)讀取經(jīng)緯度、太陽天頂角、太陽方位角、衛(wèi)星天頂角、衛(wèi)星方位角和高程等信息。由定標(biāo)系數(shù)將相應(yīng)的數(shù)據(jù)轉(zhuǎn)換成真實(shí)物理值。利用海陸掩碼文件進(jìn)行海陸分離，使用多波段閾值法去除云像元，進(jìn)而進(jìn)行AOD反演[16-17]。

4)氣溶膠光學(xué)厚度反演。當(dāng)2.1 μm波段的表觀反射率小于0.4時，該像元可被認(rèn)為是暗像元[5]。通過角度幾何參數(shù)求得散射角，通過受氣溶膠影響較小的中紅外波段計(jì)算出NDVISWIR后[18]，由2.1 μm波段的反射率獲得0.47 μm和0.66 μm波段的地表反射率：

(2)

式中：a、b是散射角與NDVISWIR的函數(shù)。

按照上式計(jì)算出紅藍(lán)波段的地表反射率；在構(gòu)建的查找表中線性內(nèi)插由地理定位數(shù)據(jù)讀取的角度參數(shù)，得到不同波段和不同AOD的3個大氣參數(shù)S、ρ0和T(μS)T(μV)；將反射率和大氣參數(shù)代入式(1)得到表觀反射率，再對真實(shí)的表觀反射率進(jìn)行插值，最終得到AOD。

深藍(lán)算法原理與暗像元原理類似，根據(jù)過境時間提取對應(yīng)的藍(lán)波段地表反射率，利用查找表插值得到AOD。

2.3 合成處理

2種算法通過6S查找表均得到0.55 μm波長處的AOD。使用優(yōu)選的融合方式對2種方法反演的AOD產(chǎn)品進(jìn)行合成，以暗像元算法為主導(dǎo)，對暗像元算法未能反演的區(qū)域使用深藍(lán)AOD進(jìn)行填充，并使用以5×5模板窗口為基礎(chǔ)的掩膜平滑處理。在窗口內(nèi)以中心像素(i，j)為基準(zhǔn)點(diǎn)，制作共3種形狀的9個掩膜窗口，分別計(jì)算各窗口內(nèi)的灰度值方差(圖2)。該方法基于的原則是掩膜窗口方差越小，該窗口的像素值越均勻，中心像素屬于該窗口的可能性越大。取方差值最小的掩膜窗口進(jìn)行平均，用該均值代替中心像素點(diǎn)的灰度值。對于噪聲點(diǎn)，該方法可以對其有效平滑，又能較好地保留面狀地物的邊界[19-20]。

圖2 選擇式掩膜窗口

方差的計(jì)算公式為：

(3)

均值的計(jì)算公式為：

(4)

式中：N為各掩膜對應(yīng)的像素個數(shù)。

最終結(jié)果以Tiff格式輸出，整體流程見圖3。

圖3 氣溶膠反演流程

3 結(jié)果

本研究根據(jù)上文方法，對2016—2018年冬季1—2月的MODIS L1B數(shù)據(jù)對京津冀地區(qū)進(jìn)行了連續(xù)監(jiān)測。

3.1 精度驗(yàn)證

選擇AERONET氣溶膠監(jiān)測網(wǎng)北京站和香河站的氣溶膠觀測數(shù)據(jù)對本文的反演結(jié)果進(jìn)行驗(yàn)證[21]。除AERONET尚未提供下載的2018年香河站相關(guān)數(shù)據(jù)以外，期間共匹配194對有效數(shù)據(jù)。將AERONET得到的AOD轉(zhuǎn)換到與反演一致的550 nm波段，以MODIS氣溶膠產(chǎn)品的誤差Δτ=±0.05±0.2τ[22]作為誤差線，并與只使用暗像元算法的反演結(jié)果作進(jìn)一步對比(圖4、圖5、表2)。

圖4 暗像元-深藍(lán)AOD與地基AOD對比圖

圖5 觀測周期匹配散點(diǎn)變化對比圖

表2 反演結(jié)果驗(yàn)證

由表2可見，暗像元-深藍(lán)反演合成AOD與地基AOD表現(xiàn)出了較高好相關(guān)性，R2約為0.79，整體結(jié)果略高于地基AOD，當(dāng)AOD>1時，監(jiān)測精度走低，絕對誤差和相關(guān)系數(shù)低于整體精度水平；AOD>1.2時，由圖4可見，多數(shù)監(jiān)測結(jié)果明顯偏高而脫離了誤差線范圍。相比之下，暗像元算法的有效匹配數(shù)降低至89對，其中香河站僅有33對，且只集中在AOD較小區(qū)域(<1)，相對誤差與相關(guān)系數(shù)明顯偏低，說明了暗像元算法在冬季的應(yīng)用存在較大局限性，合成AOD結(jié)果匹配點(diǎn)數(shù)量的明顯增加表現(xiàn)了深藍(lán)算法與暗像元算法在冬季的高互補(bǔ)性。

注：該圖基于國家測繪地理信息局標(biāo)準(zhǔn)地圖服務(wù)網(wǎng)站下載的審圖號為GS(2016)1666號的標(biāo)準(zhǔn)地圖制作，底圖無修改。圖6 2018年2月26日AQUA/MODIS真彩影像

以云量較少的一次典型灰霾天氣(2018年2月26日)為例對暗像元算法與深藍(lán)算法監(jiān)測效果進(jìn)行進(jìn)一步比較，從AQUA/MODIS真彩影像(圖6)紅線標(biāo)記區(qū)域附近可明顯看出京津冀中南部的大范圍灰霾現(xiàn)象[23]。圖7分別展示了該日暗像元、深藍(lán)及合成AOD結(jié)果，可見深藍(lán)算法對冬季植被稀疏地表尤其是灰霾區(qū)域AOD監(jiān)測具有極大優(yōu)勢，而暗像元算法難以實(shí)現(xiàn)此類區(qū)域的有效監(jiān)測。

注：該圖基于國家測繪地理信息局標(biāo)準(zhǔn)地圖服務(wù)網(wǎng)站下載的審圖號為GS(2016)1666號的標(biāo)準(zhǔn)地圖制作，底圖無修改。圖7 使用不同算法的AOD反演結(jié)果(2018年2月26日)

3.2 與MODIS AOD產(chǎn)品對比驗(yàn)證

為了進(jìn)一步驗(yàn)證本研究在監(jiān)測京津冀地區(qū)AOD空間分布上的準(zhǔn)確性和適用性，分別選擇了2016—2018年3次較為典型的冬季灰霾天氣(2016年1月10日、2017年2月13日、2018年2月26日)，利用MODIS MYD04_L2暗像元-深藍(lán)AOD融合產(chǎn)品與本研究的暗像元-深藍(lán)AOD反演結(jié)果進(jìn)行對比驗(yàn)證。

由圖8可見，本研究反演合成的AOD與MYD04_L2 AOD產(chǎn)品在空間分布上具有較高的一致性，AOD的高、低值區(qū)均表現(xiàn)出較為一致的分布趨勢。利用ArcGIS生成隨機(jī)點(diǎn)的方式，對2種AOD進(jìn)行了統(tǒng)計(jì)分析，統(tǒng)計(jì)結(jié)果見表3。

注：該圖基于國家測繪地理信息局標(biāo)準(zhǔn)地圖服務(wù)網(wǎng)站下載的審圖號為GS(2016)1666號的標(biāo)準(zhǔn)地圖制作，底圖無修改。圖8 AOD反演結(jié)果與MYD04_L2同期AOD對比圖

針對2016年1月10日、2017年2月13日、2018年2月26日3次AOD監(jiān)測，本研究監(jiān)測結(jié)果在ArcGIS生成的隨機(jī)點(diǎn)中，有效AOD點(diǎn)位占比分別達(dá)到69%、98%、80%，遠(yuǎn)好于MYD04_L2同期產(chǎn)品的22%、74%、62%；在有效監(jiān)測匹配點(diǎn)的對比中，二者呈現(xiàn)出了較為顯著的相關(guān)性，相關(guān)系數(shù)均達(dá)0.9以上。

表3 反演結(jié)果與MODIS產(chǎn)品對比驗(yàn)證

通過本研究與MODIS AOD產(chǎn)品在京津冀地區(qū)冬季的AOD監(jiān)測結(jié)果對比分析中可得，本研究AOD監(jiān)測結(jié)果的空間分布趨勢一致性高，有效監(jiān)測值的相關(guān)性強(qiáng)；監(jiān)測盲區(qū)顯著減少，適用性更高；且空間分辨率提高到了1 km，AOD空間分布的表現(xiàn)更為具體，整體監(jiān)測效果優(yōu)于MODIS AOD產(chǎn)品。

3.2 京津冀AOD的空間特征與影響因素分析

研究表明，輕度灰霾天氣的AOD均值在0.8左右[24]。以0.8為閾值統(tǒng)計(jì)了京津冀地區(qū)2016—2018年1—2月高AOD天數(shù)分布情況(圖9(a))。能夠看出京津冀地區(qū)AOD的空間分布呈現(xiàn)較為明顯的西北低、東南高，且分布圖的均質(zhì)性整體表現(xiàn)較好，能夠?qū)⑽鞅辈康膹埣铱?、承德劃為低AOD區(qū)，將京津唐及以南各市劃為高AOD區(qū)，說明了京津冀地區(qū)的氣溶膠分布與地勢的顯著負(fù)相關(guān)性。從圖9(b)可見，西北部連續(xù)的高大山體形成了氣溶膠擴(kuò)散流通的天然屏障，對北方吹來的冷空氣有著較強(qiáng)的阻擋和削弱作用，南部的山西、河南、山東省部分地區(qū)也向京津冀持續(xù)輸送著大氣污染物，在山谷風(fēng)和城市熱島環(huán)流和海陸風(fēng)環(huán)流的共同作用下，京津冀平原地區(qū)氣溶膠的流通與擴(kuò)散較為緩慢，輸送過程較為復(fù)雜[25]。

注：該圖基于國家測繪地理信息局標(biāo)準(zhǔn)地圖服務(wù)網(wǎng)站下載的審圖號為GS(2016)1666號的標(biāo)準(zhǔn)地圖制作，底圖無修改。圖9 京津冀AOD的空間特征分析圖

燃煤取暖依然是人們過冬采暖的主要手段，采暖燃煤氣溶膠也是冬季氣溶膠的主要類型[26]。除北京、天津兩大直轄市外，保定、石家莊、邯鄲、唐山、滄州等東部南部城市都是河北省的人口大市，燃煤取暖氣溶膠排放量相對于西北部人口稀疏的張家口、承德市明顯更高。北京、天津、唐山等傳統(tǒng)的工業(yè)城市和各類企業(yè)的聚集地也通常是AOD的偏高區(qū)域。表明了能源消耗、人類活動對氣溶膠也分布也有著較大影響[27]。

4 結(jié)束語

本文發(fā)展了冬季氣溶膠光學(xué)厚度遙感反演算法，該算法以暗像元算法和深藍(lán)算法為基礎(chǔ)，對近3年1—2月京津冀地區(qū)大氣氣溶膠進(jìn)行了連續(xù)監(jiān)測，通過與AERONET北京和香河地基站實(shí)測數(shù)據(jù)的對比驗(yàn)證，評價(jià)了算法的性能，并根據(jù)連續(xù)監(jiān)測結(jié)果和輔助資料分析了冬季京津冀地區(qū)AOD的空間特征及其影響因素，主要結(jié)論如下：

本算法相較于傳統(tǒng)的暗像元算法，有效觀測數(shù)明顯增加，且具備較為良好的監(jiān)測精度，說明了深藍(lán)算法與暗像元算的較強(qiáng)的互補(bǔ)性，是暗像元算法的有效補(bǔ)充，對2種算法結(jié)果進(jìn)行融合處理在冬季是一種良好的監(jiān)測手段。在與MYD04_L2氣溶膠產(chǎn)品的對比中，二者對AOD的空間分布趨勢及AOD的監(jiān)測值均表現(xiàn)出了較高一致性，但本算法對冬季京津冀地區(qū)的監(jiān)測結(jié)果擁有更大的有效范圍，且空間分辨率提高到了1 km，監(jiān)測效果更好。

當(dāng)AOD較高時，本算法可能出現(xiàn)AOD高估現(xiàn)象而使監(jiān)測結(jié)果超出Δτ=±0.05±0.2τ的誤差要求，但相較于暗像元算法在冬季的監(jiān)測表現(xiàn)[5]，本算法的監(jiān)測效果在精度與適用性方面均有明顯提高。

通過對連續(xù)監(jiān)測結(jié)果的統(tǒng)計(jì)，分析了冬季京津冀地區(qū)AOD的空間分布情況。京津冀冬季AOD整體分布呈現(xiàn)較為明顯的西北低、東南高，與地勢呈顯著負(fù)相關(guān)，較為復(fù)雜的地理?xiàng)l件使得華北平原氣溶膠積累容易，擴(kuò)散緩慢。此外，冬季的燃煤取暖等人為氣溶膠排放活動也顯著影響著AOD的分布。

致謝：MODIS數(shù)據(jù)產(chǎn)品由NASA提供，地基太陽光度計(jì)站點(diǎn)驗(yàn)證數(shù)據(jù)由AERONET提供，在此表示感謝！