張鈺萌，陳輝，馬鵬飛，厲青，王中挺

(1.河南理工大學(xué) 測(cè)繪與國(guó)土信息工程學(xué)院，河南 焦作 454000；2.生態(tài)環(huán)境部衛(wèi)星環(huán)境應(yīng)用中心，北京 100094)

0 引言

華北地區(qū)是指位于中國(guó)北部的區(qū)域，包括北京、天津、河北、山西和內(nèi)蒙古中部。該地區(qū)面積約84×104km2，人口約1.68 億，2015年GDP為10 萬(wàn)億元。該區(qū)域聚集有大量鋼鐵、化工、水泥等重污染產(chǎn)業(yè)，大量排放著氮氧化物、粉塵等大氣污染物，而普遍缺少治污措施的“散亂污”企業(yè)也在排放大氣污染物，同時(shí)，華北地區(qū)機(jī)動(dòng)車(chē)總數(shù)超過(guò)了2 600萬(wàn)，成為大氣污染物的主要來(lái)源之一。研究表明[1]，較為嚴(yán)重的大氣污染，使得我國(guó)北方的人均預(yù)期壽命將減少5.5 年，同時(shí)也大大提高了肺癌、心臟病等重大疾病的發(fā)病率。

利用遙感數(shù)據(jù)獲取華北地區(qū)氣溶膠的空間分布，不僅能獲得氣溶膠污染的空間分布狀況，而且能為顆粒物指標(biāo)的反演提供數(shù)據(jù)源[2]。針對(duì)陸地地表，由于衛(wèi)星接收到的信號(hào)中地表反射較強(qiáng)，需要去除地表反射率貢獻(xiàn)得到大氣信息，從而反演得到氣溶膠，目前主要的算法有：假定濃密植被在紅藍(lán)波段存在線(xiàn)性關(guān)系的暗目標(biāo)法[3-5]、假定地表反射隨時(shí)間變化較慢的深藍(lán)算法[6-8]、利用氣溶膠散射的角度變化的多角度算法[9-11]、采用偏振信號(hào)的多角度偏振法[12-15]。在暗目標(biāo)法基礎(chǔ)上，可以將紅藍(lán)等波段存的線(xiàn)性關(guān)系擴(kuò)展到低植被覆蓋的地表，從而實(shí)現(xiàn)多種地表覆蓋類(lèi)型的氣溶膠反演[16-17]。

搭載于美國(guó)的Terra以及Aqua星的MODIS傳感器已持續(xù)獲取數(shù)據(jù)近20 年，擁有覆蓋可見(jiàn)、近紅外、短波紅外、熱紅外等36個(gè)波段，目前該數(shù)據(jù)已廣泛應(yīng)用于氣象、環(huán)境、農(nóng)業(yè)等方面。據(jù)不完全統(tǒng)計(jì)，我國(guó)已在氣象、環(huán)境、農(nóng)業(yè)等多個(gè)領(lǐng)域建立有數(shù)十個(gè)能直接接收MODIS數(shù)據(jù)的直收站。

然而，在利用MODIS數(shù)據(jù)進(jìn)行華北地區(qū)氣溶膠反演時(shí)，還存在著以下兩個(gè)問(wèn)題：一是冬季時(shí)，我國(guó)北方地區(qū)植被稀疏，暗目標(biāo)較少，采用暗目標(biāo)法較為困難，且精度偏低[18]；二是MODIS傳感器已超期服役，器件開(kāi)始老化，直收站接收到的數(shù)據(jù)在藍(lán)光波段存在著明顯的條帶。因此，本文在Wang等[17]提出的算法基礎(chǔ)上，結(jié)合短波紅外波段從MODIS數(shù)據(jù)監(jiān)測(cè)華北地區(qū)秋陸地氣溶膠。

1 基本原理與算法

假定大氣層中的氣溶膠、大氣分子水平分層分布，傳感器接收到的反射率信號(hào)ρTOA為[19]：

(1)

式中：S、ρ0和T是3個(gè)代表大氣狀況參數(shù)，去除大氣分子影響即可反演氣溶膠光學(xué)厚度(AOD)。

除水體與雪外，大部分地表在短波紅外波段(2.12 μm)的反射率與紅波段(0.66 μm)反射率存在著線(xiàn)性關(guān)系[20]：

rSWIR=a×rred+b

(2)

式中：a、b為經(jīng)驗(yàn)參數(shù)，根據(jù)隨著地物類(lèi)型而變化。Kaufman等[21]認(rèn)為，b是大氣校正中對(duì)氣溶膠影響校正不足或校正過(guò)多帶來(lái)的，因此，本研究中將b統(tǒng)一設(shè)置為0。

結(jié)合式(1)，可以得到方程組：

(3)

解上述方程組反演得到AOD。具體為：使用輻射傳輸軟件計(jì)算各個(gè)AOD條件下短波紅外波段和紅波段的大氣參數(shù)S、ρ0和T；利用式(3)和式(4)，將MODIS觀測(cè)到的短波紅外波段表觀反射率ρTOA，SWIR和紅波段ρTOA，red代入，得到rSWIR和rred；將rred代入式(5)，計(jì)算得到rSWIR；計(jì)算利用式(3)和式(5)得到的rSWIR之間的差，以差為自變量，AOD為因變量，進(jìn)行線(xiàn)性擬合，求取差為0時(shí)的AOD值，即為所求。

在實(shí)際反演中，為加快運(yùn)算速度，本文預(yù)先使用6S vector1.0[22]設(shè)定AOD從0～4變化，進(jìn)行輻射傳輸運(yùn)算存儲(chǔ)在查找表中。根據(jù)MODIS的過(guò)境時(shí)間和觀測(cè)特征，設(shè)置太陽(yáng)天頂角為0～60°，觀測(cè)天頂角0～60°，相對(duì)方位角為0～180°。氣溶膠類(lèi)型的設(shè)定參照了地面觀測(cè)的結(jié)果[23-25]。華北地區(qū)的氣溶膠由3種基本氣溶膠混合而成，以水溶型為主，沙塵型和煤煙型為輔，設(shè)定各氣溶膠組分比例為：沙塵型15%、水溶型83%、黑炭型2%。

2 地表反射率分析

為分析華北地區(qū)各地表類(lèi)型的短波紅外波段的反射率與紅波段反射率的情況，本文收集了2015年MODIS的地表反射率產(chǎn)品，該產(chǎn)品校正了氣溶膠、大氣分子等對(duì)衛(wèi)星觀測(cè)反射率的影響[26]。提取2015年京津冀區(qū)域的紅光、近紅外和短波紅外波段數(shù)據(jù)共計(jì)46景。然后，計(jì)算每景數(shù)據(jù)的紅光和短波紅外波段的比率和歸一化植被指數(shù)(NDVI)。為消除云、雪等因素的影響，去除了地表反射率大于0.2以及NDVI小于0的數(shù)據(jù)。

隨機(jī)選取了2015年度所有時(shí)間的城市、農(nóng)田、森林、草地共4種典型土地類(lèi)型的反射率比率、NDVI、紅光波段反射率，分析它們隨時(shí)間變化的情況。

圖1表明紅光和短波紅外地表反射率比率整體較為穩(wěn)定，隨著時(shí)間變化較小，城市最高，在0.8左右；農(nóng)田次之，在0.7左右，7—9月有明顯的低值出現(xiàn)；草地更低，為0.6左右，隨時(shí)間變化最小；森林最小，在0.4左右。

圖1 城市、農(nóng)田、森林、草地反射率比率隨時(shí)間變化圖

從圖2可以看出，各典型地物的NDVI有較強(qiáng)季節(jié)變化規(guī)律，夏季植物茂盛達(dá)到最高，冬季植被凋零降到最低。森林整體最高，5—10月最高，在0.6以上，其他時(shí)間在0.4以下；農(nóng)田次之，7—10月最高，在0.6以上，4—6月在0.4左右，其他在0.2左右；草地和城市較為接近，6—10月較高，在0.2以上。

圖2 2015年度的城市、農(nóng)田、森林、草地NDVI

從圖3可以看出，各典型地物在紅光波段的反射率也體現(xiàn)了季節(jié)變化規(guī)律。草地整體最高，在0.148左右；農(nóng)田次之，7—9月最低，3—5月也較低；城市更低，7—9月達(dá)到最低；森林最低，2—4月最高，8—10月最低。

圖3 2015年度的城市、農(nóng)田、森林、草地紅波段反射率

表1為各典型地物反射率比值、NDVI、紅波段反射率統(tǒng)計(jì)結(jié)果。本文用標(biāo)準(zhǔn)差占平均值的百分比來(lái)衡量不同指標(biāo)隨季節(jié)變化的離散程度。可以看出，反射率比率的離散程度最小，百分比在15%以下；NDVI離散程度最大，百分比都在30%以上。因此，可以將紅光和短波紅外波段反射率比率的年平均值作為經(jīng)驗(yàn)系數(shù)實(shí)現(xiàn)陸地氣溶膠的反演。

表1 各典型地物反射率比值、NDVI、紅波段反射率統(tǒng)計(jì)表

因此，本文將2015年的紅波段與短波紅外波段地表反射率比率作年平均，得到華北地區(qū)的反射率比率分布圖，如圖4所示?？梢钥闯?，華北地區(qū)的反射率比率隨地物類(lèi)型有明顯的不同，河北北部的燕山、山西與河北交界的太行山等以森林覆蓋為主等山區(qū)較低，在0.6以下，占據(jù)華北平原大部分地區(qū)的農(nóng)田區(qū)域則在0.6～0.8之間，部分城市密集區(qū)以及河流高于0.8，內(nèi)蒙古的草原區(qū)域在0.6左右。

注：該圖基于國(guó)家測(cè)繪地理信息局標(biāo)準(zhǔn)地圖服務(wù)網(wǎng)站下載的審圖號(hào)為GS(2017)1268號(hào)的標(biāo)準(zhǔn)地圖制作，底圖無(wú)修改。圖4 華北地區(qū)年均反射率比率分布圖

根據(jù)圖4，本文統(tǒng)計(jì)了紅波段與短波紅外波段地表反射率比值出現(xiàn)的頻率，如圖5所示?？梢钥闯觯饕植荚?.4～0.7之間，其中，值在0.5和0.6出現(xiàn)的頻率最多，超過(guò)了40%。

圖5 華北地區(qū)年均反射率比值分布圖

3 結(jié)果與驗(yàn)證

本文將上述分析得到的紅波段與短波紅外波段地表反射率比率年均值作為經(jīng)驗(yàn)系數(shù)代入公式(3)至公式(5)，從而完成氣溶膠的反演。具體的反演過(guò)程采用IDL語(yǔ)言實(shí)現(xiàn)。

3.1 反演實(shí)驗(yàn)

本文選取了2016年9月—2017年8月期間過(guò)境京津冀地區(qū)的MODIS數(shù)據(jù)進(jìn)行了算法反演實(shí)驗(yàn)。圖6為2016年10月13日、2017年1月17日、2017年3月19日、2017年5月4日4天比較典型的氣溶膠狀況反演獲得的反演結(jié)果。由于可見(jiàn)紅外波段無(wú)法穿透云層，本算法在云覆蓋區(qū)域失效；同時(shí)，在冰雪覆蓋區(qū)域，由于地表反射率過(guò)高，本算法也無(wú)法獲得有效結(jié)果。

注：該圖基于國(guó)家測(cè)繪地理信息局標(biāo)準(zhǔn)地圖服務(wù)網(wǎng)站下載的審圖號(hào)為GS(2017)1268號(hào)的標(biāo)準(zhǔn)地圖制作，底圖無(wú)修改。圖6 本算法的AOD反演結(jié)果圖

如圖6所示，2016年10月13日，華北地區(qū)的空氣狀況整體較好，僅在北京東南部、天津和河北中部出現(xiàn)了較高的氣溶膠，部分地區(qū)高于1。2017年1月17日，華北地區(qū)的空氣狀況東南部出現(xiàn)了污染，較高的氣溶膠主要分布在河北南部、河南北部、山東西北部，在河北河南交界處出現(xiàn)大片AOD大于2的區(qū)域，在河北北部有部分地區(qū)沒(méi)有反演結(jié)果，可能是該區(qū)域?qū)儆谏絽^(qū)，被雪覆蓋后沒(méi)有融化，雪的反射率過(guò)高導(dǎo)致算法失效。2017年3月19日，華北地區(qū)大部地區(qū)都出現(xiàn)了較重的空氣污染，整體AOD在1以上，河北中南部AOD較高，在2以上，河北南部和河南由于被云覆蓋沒(méi)有獲得反演結(jié)果。2017年5月4日，華北地區(qū)大部分地區(qū)受到沙塵天氣，整體上環(huán)境大氣質(zhì)量較差，在河北中部出現(xiàn)了污染聚集團(tuán)，AOD在2以上，但在中間最重區(qū)域算法失效，存在部分空白區(qū)域，原因是高濃度氣溶膠與云的反射特征較為相近，本文將該處誤識(shí)別為云。

3.2 結(jié)果驗(yàn)證

本文收集了衛(wèi)星過(guò)境時(shí)AERONET北京Radi站和香河站氣溶膠產(chǎn)品以驗(yàn)證本算法精度。北京Radi站位于北京市區(qū)北部，周?chē)写罅烤用顸c(diǎn)和商業(yè)區(qū)，屬于典型的城市區(qū)域，香河站位于北京東南的香河縣，周?chē)赞r(nóng)田為主。AERONET是對(duì)氣溶膠特性進(jìn)行地面觀測(cè)的網(wǎng)絡(luò)，提供全球的氣溶膠光學(xué)厚度、譜分布和散射相函數(shù)[27]觀測(cè)結(jié)果。AERONET一般提供440 nm、675 nm等波段的氣溶膠光學(xué)厚度，本文采用Angstrom公式[28]將AOD轉(zhuǎn)化為550 nm結(jié)果，與本文結(jié)果相比較。

去除云、霧、雨、雪等天氣影響的數(shù)據(jù)，在2016年9月1日至2017年8月31日TERRA和AQUA衛(wèi)星過(guò)境的365 d中，本算法在北京Radi站共178 d獲得有效驗(yàn)證數(shù)據(jù)，本算法在香河站共191 d獲得有效驗(yàn)證數(shù)據(jù)(圖7、圖8)。

圖7 本文結(jié)果與AERONET氣溶膠產(chǎn)品時(shí)間變化(香河站)

圖8 本文結(jié)果與AERONET氣溶膠產(chǎn)品時(shí)間變化(北京Radi站)

從圖7和圖8可以看出，本算法獲得的結(jié)果與地面結(jié)果較好地貼合在一起，但在AOD較高時(shí)部分反演結(jié)果與地面觀測(cè)的結(jié)果有較大差異。而時(shí)間序列的觀測(cè)結(jié)果也部分體現(xiàn)了季節(jié)變化的特征，香河站5—8月AOD明顯低于其他月份，北京站6—8月的AOD也低于其他季節(jié)，但中間存在若干高AOD日。

圖9 本文反演AOD結(jié)果與AERONET散點(diǎn)圖(香河站)

圖10 本文反演AOD結(jié)果與AERONET散點(diǎn)圖(北京Radi站)

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

注：r1為誤差線(xiàn)范圍內(nèi)反演結(jié)果的比例；r2為比地面結(jié)果大同時(shí)不在誤差線(xiàn)范圍內(nèi)反演結(jié)果的比例；r3為比地面結(jié)果小同時(shí)不在誤差線(xiàn)范圍內(nèi)反演結(jié)果的比例。

3.3 分季節(jié)驗(yàn)證

本文分春、夏、秋、冬4個(gè)季節(jié)開(kāi)展了算法驗(yàn)證，結(jié)果見(jiàn)圖11、圖12和表3。從圖中可以看出，本算法在各個(gè)季節(jié)獲得的結(jié)果都較為理想，相關(guān)系數(shù)都在0.9以上，RMSE也基本在0.1左右，在算法整體低估的情況下，秋季低估的狀況最為明顯；其中，香河站春季效果最好，夏季效果最差；北京站秋季和冬季效果較好，春季效果最差?？赡苁潜疚乃玫臍馊苣z模型，在春季時(shí)，跟香河站的情況較為一致，但低估了北京地區(qū)的沙塵含量。

圖11 本文反演AOD與AERONET產(chǎn)品之間分季節(jié)比較散點(diǎn)圖(香河站)

圖12 本文反演AOD與AERONET產(chǎn)品之間分季節(jié)比較散點(diǎn)圖(北京Radi站)

表3 反演結(jié)果分季節(jié)驗(yàn)證表

4 結(jié)束語(yǔ)

本文探討了通過(guò)假定短波紅外波段與紅光波段的地表反射率之間比率不變來(lái)獲取華北地區(qū)陸地上空的氣溶膠，主要取得如下結(jié)論：①與NDVI、紅波段地表反射率相比，短波紅外波段與紅波段的地表反射率的比值較為穩(wěn)定，利用該比值能夠去除地表影響反演得到陸地氣溶膠。②利用本文算法能夠較好地獲得華北地區(qū)城市和鄉(xiāng)村區(qū)域的氣溶膠分布，地面驗(yàn)證表明70%的結(jié)果誤差范圍在±0.05±0.15τ之內(nèi)，且與地面結(jié)果的相關(guān)系數(shù)在0.9以上。

但本算法還需要在2個(gè)方面做進(jìn)一步改進(jìn)：①氣溶膠各組分的比例固定，雖然整體也取得了不錯(cuò)的反演效果，但隨著地表類(lèi)型、季節(jié)的變化，會(huì)對(duì)反演結(jié)果帶來(lái)±10%左右的偏差(圖10、圖11)，在今后的反演中，還需借助紫外、偏振等多源衛(wèi)星數(shù)據(jù)提取氣溶膠組分信息以提高反演精度。②本文算法將每個(gè)像元的紅波段與短波紅外波段地表反射率比率固定為其年均值，但還存在著一定的變化，統(tǒng)計(jì)表明(表1)，其與年均值偏移量在15%左右，考慮到紅波段的反射在0.1左右，反射率比值在0.5～0.6，則其對(duì)地表反射率的誤差在0.01左右，會(huì)造成約0.1氣溶膠光學(xué)厚度誤差[3]。

致謝：本文的MOD09數(shù)據(jù)由美國(guó)國(guó)家航空航天局提供，AERONET北京Radi站數(shù)據(jù)由中科院遙地所提供，AERONET香河站數(shù)據(jù)由中科院大氣所提供，6S Vector1.1程序由Vermote E等提供，在此表示感謝。