京津冀冬季氣溶膠光學厚度遙感反演

孫灝，劉偉漢，王艷梅，周偉，蔡創創

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

0 引言

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

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

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

1 研究區域和數據

1.1 研究區概況

注：該圖基于國家測繪地理信息局標準地圖服務網站下載的審圖號為GS(2016)1666號的標準地圖制作，底圖無修改。圖1 研究區域

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

1.2 數據來源

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

2 研究方法

2.1 反演原理

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

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

(1)

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

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

2.2 反演過程

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

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

2)構建深藍算法地表反射率庫。針對在高地表反射率地區，紅藍波段和近紅外波段地表反射率無法構成線性關系，通過構建地表反射率庫實現地氣解耦。本研究以MODIS地表反射率產品MOD09 8天合成產品為基礎構建地表反射率庫，按日期存儲在數據文件中。

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

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

(2)

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

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

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

2.3 合成處理

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

圖2 選擇式掩膜窗口

方差的計算公式為：

(3)

均值的計算公式為：

(4)

式中：N為各掩膜對應的像素個數。

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

圖3 氣溶膠反演流程

3 結果

本研究根據上文方法，對2016—2018年冬季1—2月的MODIS L1B數據對京津冀地區進行了連續監測。

3.1 精度驗證

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

圖4 暗像元-深藍AOD與地基AOD對比圖

圖5 觀測周期匹配散點變化對比圖

表2 反演結果驗證

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

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

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

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

3.2 與MODIS AOD產品對比驗證

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

由圖8可見，本研究反演合成的AOD與MYD04_L2 AOD產品在空間分布上具有較高的一致性，AOD的高、低值區均表現出較為一致的分布趨勢。利用ArcGIS生成隨機點的方式，對2種AOD進行了統計分析，統計結果見表3。

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

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

表3 反演結果與MODIS產品對比驗證

通過本研究與MODIS AOD產品在京津冀地區冬季的AOD監測結果對比分析中可得，本研究AOD監測結果的空間分布趨勢一致性高，有效監測值的相關性強；監測盲區顯著減少，適用性更高；且空間分辨率提高到了1 km，AOD空間分布的表現更為具體，整體監測效果優于MODIS AOD產品。

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

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

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

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

4 結束語

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

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

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

通過對連續監測結果的統計，分析了冬季京津冀地區AOD的空間分布情況。京津冀冬季AOD整體分布呈現較為明顯的西北低、東南高，與地勢呈顯著負相關，較為復雜的地理條件使得華北平原氣溶膠積累容易，擴散緩慢。此外，冬季的燃煤取暖等人為氣溶膠排放活動也顯著影響著AOD的分布。

致謝：MODIS數據產品由NASA提供，地基太陽光度計站點驗證數據由AERONET提供，在此表示感謝！