劉宇昕,王斌,曹詩頌,張峰,王兆徽
(1.國家衛(wèi)星海洋應用中心,北京 100081;2.自然資源部 空間海洋遙感與應用重點實驗室,北京 100081;3.河海大學 海洋學院,南京 210098;4.中國氣象局武漢暴雨研究所,武漢 430074;5.北京建筑大學 測繪與城市空間信息學院,北京 102616;6.北京集思睿景科技有限公司,北京 102401)
近百年來,伴隨著社會工業(yè)化的發(fā)展以及城市化進程的加速,多種大氣污染源的過度排放,使城市空氣狀況受到極大影響,已經(jīng)遠遠超過城市生態(tài)承載力。大氣污染對城市人群健康以及建筑等各方面產(chǎn)生的影響,已成為我國城市和區(qū)域可持續(xù)發(fā)展需要面對的一個難題。同時,區(qū)域性的大氣污染也已成為眾多學者關注的重點,諸如京津冀、珠三角、長三角等城市群的大氣污染越來越受到人們重視,并已開展了大量相關研究。此外,在中部地區(qū)如武漢及周邊區(qū)域,工業(yè)化和城市化帶來的大氣污染也是環(huán)境管理部門所面臨的棘手問題。在諸多大氣污染物之中,氣溶膠與大氣環(huán)境狀況密切相關,氣溶膠顆粒的增加可引發(fā)霧霾天氣的產(chǎn)生,嚴重影響人們的身體健康和正常生活。
衛(wèi)星遙感可以迅速、準確、動態(tài)地監(jiān)測較大區(qū)域的大氣狀況及變化,用以及時研究相關應對方案和規(guī)劃,降低污染給人們所帶來的影響。氣溶膠光學厚度(aerosol optical depth,AOD) 作為衡量城市空氣質(zhì)量的重要指示因子,對大氣污染情況的判定以及定量遙感研究都具有很大意義。作為氣溶膠觀測的一種有效方法,氣溶膠光學厚度遙感定量反演研究至今已有40多年的歷史,多年來針對不同種類的傳感器,開發(fā)出了多種氣溶膠光學厚度反演算法。如:基于改進的甚高分辨率輻射計(advanced very high resolution radiometer,AVHRR)的單通道和雙通道氣溶膠光學厚度反演算法;基于總臭氧測繪光譜儀(total ozone mapping spectrometer,TOMS)的氣溶膠指數(shù)和近紫外光通道算法;基于寬視場海洋觀測傳感器(sea-viewing wide field-of-view sensor,SeaWiFS)、中分辨率成像光譜輻射計(moderate resolution imaging spectroradiometer,MODIS)、中分辨率成像光譜儀(medium-resolution imaging spectrometer,MERIS)的多通道算法;基于偏振探測器(polarization and directionality of the Earth’s reflectance,POLDER),以及暗目標法、深藍算法、結(jié)構(gòu)函數(shù)法、MAIAC(multi-angle implementation of atmospheric correction)算法[1]、SARA(simplified aerosol retrieval algorithm)算法[2]等氣溶膠光學厚度反演算法。結(jié)構(gòu)函數(shù)法是目前普遍采用的一種氣溶膠反演算法,其在陸地衛(wèi)星專題測圖儀(landsat thematic mapper,Landsat TM)、MODIS、AVHRR以及環(huán)境(HJ)衛(wèi)星數(shù)據(jù)中都得到了廣泛的應用。
使用結(jié)構(gòu)函數(shù)法反演氣溶膠光學厚度至今已有30多年的歷史。在國外研究方面,Tanré等[3]率先提出了結(jié)構(gòu)函數(shù)算法用于干旱地區(qū)的氣溶膠滾光學厚度反演,并使用該方法成功反演出Saharan沙漠地區(qū)的氣溶膠光學厚度;Holben等[4]使用結(jié)構(gòu)函數(shù)法實現(xiàn)了AVHRR等數(shù)據(jù)的氣溶膠光學厚度反演;Levy等[5]使用MODIS數(shù)據(jù)反演氣溶膠光學厚度,并對結(jié)構(gòu)函數(shù)法進行了相應的改進。國內(nèi)研究方面,Liu等[6]使用結(jié)構(gòu)函數(shù)法實現(xiàn)了地球觀測衛(wèi)星系統(tǒng)(systeme probatoire d’observation de la terre,SPOT)數(shù)據(jù)的氣溶膠光學厚度反演并判定大氣污染狀況;李曉靜等[7]利用MODIS數(shù)據(jù)運用結(jié)構(gòu)函數(shù)法,成功反演得到北京和周邊地區(qū)的氣溶膠光學厚度;黃艇[8]以蘭州地區(qū)為例,結(jié)合MODIS數(shù)據(jù),使用結(jié)構(gòu)函數(shù)法,反演得到了該地區(qū)的氣溶膠光學厚度,并與多種方法的結(jié)果進行比較;孫林[9]使用改進的結(jié)構(gòu)函數(shù)法,并結(jié)合雙向反射分布函數(shù)(bidirectional reflectance distribution function,BRDF)模型,針對MODIS數(shù)據(jù)成功反演得到北京地區(qū)的氣溶膠光學厚度,提高了反演精度;周春艷等[10]針對高分辨率HJ衛(wèi)星數(shù)據(jù)使用結(jié)構(gòu)函數(shù)法獲取了北京地區(qū)的氣溶膠光學厚度;朱琳等[11-12]研究了使用已有的地表反射率產(chǎn)品為待反演圖像提供清晰圖像的方法,以推動結(jié)構(gòu)函數(shù)法的業(yè)務化應用,并研究獲取了最佳的像元間隔設置。武漢地處中部,地形包括低山、丘陵、平原等,土地利用類型多樣且利用程度較高,湖泊塘堰眾多,工業(yè)發(fā)達且市內(nèi)較多重工業(yè),大氣污染較嚴重,但武漢這一較具代表性區(qū)域還缺乏氣溶膠光學厚度反演相關研究。相關研究已經(jīng)證明,對于城市等亮地表地區(qū),結(jié)構(gòu)函數(shù)法在反演TM、MODIS等不同分辨率載荷的AOD方面都有較好的結(jié)果。本文選擇武漢市作為研究區(qū)域,使用結(jié)構(gòu)函數(shù)法反演環(huán)境一號(HJ-1)A/B衛(wèi)星的氣溶膠光學厚度,并對結(jié)果進行了驗證。
由于衛(wèi)星傳感器接收到的輻射能量來源于大氣散射以及地表反射,當?shù)乇矸瓷浜苌贂r,衛(wèi)星觀測到的信號主要來自于大氣部分;當?shù)乇矸瓷湓龃髸r,地表反射又成為主要部分。通過衛(wèi)星遙感數(shù)據(jù)反演氣溶膠光學厚度方法的實質(zhì),就是將來自于大氣和地表反射的輻射分離開來,也就是地氣解耦。
在地表反射率較高的區(qū)域,地表反射在衛(wèi)星接收到的能量中占主要部分,大氣氣溶膠的散射較弱,所以大氣氣溶膠信息會因使用地表反射率模型帶來的誤差而丟失。不同地點地表覆蓋類型的差異,通過遙感反射率在空間上的變化率反映出來,就是此區(qū)域的結(jié)構(gòu)函數(shù)。結(jié)構(gòu)函數(shù)法通過大氣透過率得到氣溶膠光學厚度,能夠較好解決地表反射率模型誤差帶來的問題。
使用結(jié)構(gòu)函數(shù)法反演一組遙感數(shù)據(jù),其中包括了一張晴天無云狀態(tài)下氣溶膠光學厚度非常小的影像,通過計算其像元間地表反射率變化率后作為參考圖像,代表該地區(qū)地表反射率的分布規(guī)律。對于該組中其他影像,假設該地表空間分布規(guī)律在此組遙感數(shù)據(jù)對應時間內(nèi)保持不變,則該時間段內(nèi)遙感數(shù)據(jù)表觀反射率在空間分布的變化是由氣溶膠所引起,依據(jù)其分別與參考影像的差別,可以對氣溶膠光學厚度進行計算。結(jié)構(gòu)函數(shù)法由于較少受到地面反射率的限制,因此在干旱和半干旱地區(qū)的應用較為廣泛,原理如下[3]:假設大氣水平均勻,陸地表面為朗伯體,信號用表觀反射率ρ*表示,定義如式(1)所示。
(1)
式中:θs、θv、φv是觀測幾何參數(shù)中的太陽天頂角、觀測天頂角和方位角;ρa(θs,θv,φv)是大氣反射率;T(θs)是大氣下行透過率;td(θv)為環(huán)境輻射大氣上行透過率;S是大氣向下半球反照率;ρc(M)是目標物M的表面反射率;<ρc(M)>是平均地表反射率;τ是大氣光學厚度;μv為cosθv。

(2)

(3)
在上述假設條件下,Δρi,j是一個常數(shù),表觀反射率是光學厚度的一個函數(shù),給出一個地表結(jié)構(gòu)函數(shù)見式(4)。
(4)
式中:N為區(qū)域總像元數(shù),通過式(3)和式(4),可得式(5)。
(5)
式中:M*和M分別表示衛(wèi)星觀測數(shù)據(jù)和真實地表情況的結(jié)構(gòu)函數(shù)。當t1和t2的地表反射率沒有發(fā)生變化,M2(d,t1)=M2(d,t2),根據(jù)公式(5)可得式(6)。
(6)
根據(jù)經(jīng)驗,ln[T(θs)exp(-τ/μv)]與氣溶膠光學厚度成線性關系[3]。因此,當已知一個時刻的氣溶膠光學厚度時,即可以計算其他時刻的氣溶膠光學厚度。
在結(jié)構(gòu)函數(shù)法反演氣溶膠光學厚度的過程中,非常重要的一項內(nèi)容就是確定結(jié)構(gòu)函數(shù)模型,包括結(jié)構(gòu)函數(shù)公式、窗口大小以及距離值的確定。其反演原理是以朗伯體假設為出發(fā)點,使得地表反射特性不受觀測幾何或地形的影響;然而實際中大部分地表是非朗伯體,意味著氣溶膠光學厚度的反演結(jié)果受地形、觀測幾何等影響。公式(4)定義的結(jié)構(gòu)函數(shù)只考慮了一個方向,不能很好地對地表結(jié)構(gòu)進行建模,可以通過增加(ρi,j-ρi+d,j)2和(ρi,j-ρi+d,j+d)2項,將表面結(jié)構(gòu)模型從一維擴展到多維,以減小模型誤差,包含了行、列以及對角線方向,如公式(7)所示[6]。
(7)
式中:m、n確定結(jié)構(gòu)函數(shù)窗口范圍;d確定計算結(jié)構(gòu)函數(shù)值時兩像元之間的距離。很明顯,不同距離d下的氣溶膠光學厚度計算結(jié)果會有差別,計算過程中需考慮不同的結(jié)構(gòu)函數(shù)窗口,并選擇不同的距離值進行匹配。
本文使用的HJ-1A/B衛(wèi)星CCD(charge-coupled device)數(shù)據(jù),來源于環(huán)境減災衛(wèi)星,是我國專門用于環(huán)境與災害監(jiān)測的小衛(wèi)星星座,星座中包含2顆光學衛(wèi)星(HJ-1 A星和HJ-1 B星),這2顆光學衛(wèi)星上分別裝有2臺寬覆蓋多光譜可見光相機(簡稱CCD相機)。另外,HJ-1 A星上裝有一臺超光譜成像儀,HJ-1 B星上裝有一臺紅外掃描儀。環(huán)境一號系列衛(wèi)星的寬覆蓋多光譜相機、紅外相機、超光譜成像儀的重訪觀測周期分別為48、96和96 h,基本具備了在全國范圍內(nèi)開展環(huán)境和災害大范圍、全天候觀測的能力。
結(jié)構(gòu)函數(shù)法反演HJ-1 A/B衛(wèi)星遙感影像的預處理主要有輻射定標、計算表觀反射率、參考影像大氣校正、幾何校正、波段合成以及生成6S查找表等。
1)輻射定標是指建立遙感傳感器的數(shù)字量化輸出值(digital number,DN)與其所對應視場中輻射亮度值之間的定量關系,HJ衛(wèi)星數(shù)據(jù)的定標參數(shù)來自中國資源衛(wèi)星應用中心[13]。首先利用絕對定標系數(shù),將CCD影像DN值轉(zhuǎn)換為輻亮度影像,然后將輻亮度轉(zhuǎn)化為表觀反射率。
2)大氣校正部分借助ENVI(the environment for visualizing images)軟件提供的快速大氣校正工具(quick atmospheric correction,QUAC),通過從遙感圖像上采集不同地物的波譜特征,依靠得到的經(jīng)驗值實現(xiàn)大氣校正。目前QUAC支持的多光譜和高光譜波譜范圍是400~2 500 nm,輸入數(shù)據(jù)可以是輻射亮度值、表觀反射率、無單位的raw數(shù)據(jù),數(shù)據(jù)儲存格式和類型沒有特殊要求,但必須提供多光譜和高光譜傳感器數(shù)據(jù)每個波段的中心波長信息。QUAC可以利用很少的參數(shù)快速地對多光譜和高光譜數(shù)據(jù)進行處理,自動化程度較高,已被很多專家學者采用[14]。
3)幾何校正部分僅針對結(jié)構(gòu)函數(shù)法,原因是結(jié)構(gòu)函數(shù)法是只針對一個時間序列內(nèi)的影像進行氣溶膠光學厚度反演的算法,這就要求該時間序列內(nèi)的所有影像必須在地理位置上可以很好地重合,盡量保證不同影像的幾何偏移在1個像元之內(nèi)。幾何校正部分使用ERDAS 軟件的自動幾何校正模塊實現(xiàn),具體流程是:首先選取一景幾何校正后效果較好的遙感影像作為標準,并對其他影像進行校正。操作過程中,首先要手動選取3個或3個以上地理位置比較明確的控制點,借助自動校正模塊,選取其余的控制點,并將誤差超過1個像元的控制點刪除,以此為準,對該時間序列內(nèi)的所有影像進行幾何校正,誤差控制在0.5個像素以內(nèi)。
4)波段合成的目的是通過遙感影像不同波段進行組合,從而突出與氣溶膠光學厚度反演相關的信息。對于結(jié)構(gòu)函數(shù)法,選擇HJ衛(wèi)星的第1(藍光波段)、第3(紅光波段)、第4(近紅外波段)波段進行合成,然后通過影像裁剪,獲得武漢地區(qū)的多波段合成結(jié)果。
5)查找表法。6S模式是“第二代太陽短波輻射的衛(wèi)星信號模擬(second simulation of satellite signal in the solar spectrum,6S)”的簡稱,由法國大氣光學實驗室和美國馬里蘭大學地理系Veromote E,在太陽光譜波段衛(wèi)星信號模擬程序(simulation of satellite signal in the solar spectrum,5S)模型的基礎上開發(fā)的輻射傳輸模型,目前己成為世界上發(fā)展最完善的大氣輻射校正模型之一[15]。該模式主要用來模擬星載或機載遙測儀器在0.25~4.0 μm光譜上無云情況下衛(wèi)星傳感器理論上應接收到的輻射值,模式的光譜積分步長為2.5 nm[8]。6S大氣輻射傳輸模式主要包括觀測幾何參數(shù)、大氣參數(shù)、氣溶膠參數(shù)、傳感器的光譜特性以及地表反射率。根據(jù)以下6S參數(shù)設置生成查找表:9個太陽天頂角為0°、6°、12°、24°、35.2°、48°、54°、60°和66°;12個觀測天頂角為0°、6°、12°、18°、24°、30°、36°、42°、48°、54°、60°、66°;16個相對方位角為0°、12°、24°、36°、48°、60°、72°、84°、96°、108°、120°、132°、144°、156°、168°、180°;大氣模式為中緯度夏季和中緯度冬季;大氣氣溶膠模式為大陸型、城市型以及自定義類型;氣溶膠光學厚度預設值分別為0、0.25、0.5、1.0、1.5和1.95。
在遙感數(shù)據(jù)預處理結(jié)果的基礎上,反演過程中需要使用到HJ衛(wèi)星CCD影像第1、3、4波段的表觀反射率、遙感影像元數(shù)據(jù)以及生成的查找表,通過編程實現(xiàn)結(jié)構(gòu)函數(shù)法反演氣溶膠光學厚度。通過數(shù)據(jù)查詢,選取2012年9月1日至2012年10月2日中8 d的數(shù)據(jù)作為基礎數(shù)據(jù)集,具體的數(shù)據(jù)狀況如表1所示。

表1 反演數(shù)據(jù)列表
以上數(shù)據(jù)集在時間序列覆蓋整個9月份,數(shù)據(jù)基本無幾何畸變,且云覆蓋量很少。通過在湖北省環(huán)境保護廳網(wǎng)站查詢得到這段時間內(nèi)武漢市的大氣污染狀況,如圖1所示,以此為基礎,獲取相應時間的空氣污染指數(shù)(air pollution index,API)。根據(jù)空氣污染指數(shù),選擇2012年9月15日為“清潔日”,標準是選取一段時間序列內(nèi)空氣質(zhì)量次優(yōu)的一天。即空氣污染指數(shù)第二小的日期作為“清潔日”。其原理是“清潔日”要求大氣潔凈,大氣污染程度較輕,同時為避免觀測數(shù)據(jù)的隨機性和偶然性帶來的影響,因此放棄API最小的日期,選擇API第二小的日期作為“清潔日”。

圖1 2012年9月1日至2012年10月2日武漢大氣污染狀況
以2012年9月20日為例,依次對原始影像進行輻射定標、表觀反射率計算、波段融合、幾何校正和裁剪等操作,對于“清潔日”的影像,裁剪完成后需要進行大氣校正得到參考影像。實驗結(jié)果表明,武漢地區(qū)由于水體面積覆蓋范圍較大,河湖縱橫,水體信號對于氣溶膠光學厚度的反演造成了很大的影響,導致長江、東湖、南湖等水體面積較大區(qū)域的氣溶膠光學厚度明顯偏高。并且由于結(jié)構(gòu)函數(shù)法反演氣溶膠光學厚度是基于窗口進行反演,這樣水陸邊界區(qū)域就存在大量的混合像元,這些地區(qū)的氣溶膠光學厚度值也偏高,因此需要尋求合理的方案解決水體對于氣溶膠光學厚度反演過程的影響。
剔除水體影響的具體方法是使用歸一化植被指數(shù)(normalized difference vegetation index,NDVI)進行水體和陸地的判別,并對影像進行掩膜,從武漢地區(qū)的表觀反射率影像中提取陸地部分。經(jīng)過處理得到9月15日和9月20日的武漢市地區(qū)陸地區(qū)域的表觀反射率影像,如圖2所示。

圖2 2012年9月15日與2012年9月20日武漢地區(qū)陸地區(qū)域的表觀反射率合成影像
以圖2所示的2景影像為基礎,首先要根據(jù)表觀反射率,計算得到2 d的結(jié)構(gòu)函數(shù)值,經(jīng)過比較分析,結(jié)構(gòu)函數(shù)窗口選擇10、距離值選擇5,為最佳的參數(shù)組合,計算得到2天的結(jié)構(gòu)函數(shù)值影像。接下來,對查找表進行插值,結(jié)構(gòu)函數(shù)法查找表插值計算過程中需要使用到查找表中的太陽天頂角、觀測天頂角、大氣下行輻射透過率3項參數(shù),通過插值運算結(jié)合結(jié)構(gòu)函數(shù)值影像進行計算,最終得到氣溶膠光學厚度。
計算得到的2012年9月20日的氣溶膠光學厚度如圖3所示,水體部分已經(jīng)全部掩膜,對應于圖中的白色部分,陸地部分的氣溶膠光學厚度值介于0~2.0之間。

圖3 2012年9月20日氣溶膠光學厚度
對于以上反演結(jié)果,本文通過2種方式進行驗證:一是依據(jù)全自動太陽光度計CE318實測數(shù)據(jù)對反演結(jié)果進行精度驗證;二是依據(jù)湖北省環(huán)保廳提供的武漢地區(qū)大氣污染狀況,從時間上的分布趨勢對反演結(jié)果進行驗證。
1)實測數(shù)據(jù)驗證。2012年9月1日、2012年9月5日、2012年9月18日、2012年9月19日、2012年9月20日、2012年9月28日、2012年10月2日的氣溶膠光學厚度如圖4所示,選取的實測數(shù)據(jù)與影像進行匹配的時間窗口小于1 h,驗證情況如表2所示。
經(jīng)過以上驗證可知,本次反演結(jié)果的絕對誤差均在0.15以內(nèi),相對誤差均在15%以內(nèi),表明反演

圖4 氣溶膠光學厚度反演結(jié)果

表2 氣溶膠光學厚度反演精度驗證
結(jié)果具有較高精度。同時,從空間分布特征可以看出,氣溶膠光學厚度較大的區(qū)域主要包括市區(qū)以及長江沿岸,并且武漢三鎮(zhèn)中漢口和漢陽地區(qū)較高,這一趨勢與武漢市周邊地區(qū)的實際狀況基本一致。
2)時間分布趨勢驗證。針對反演結(jié)果,對武漢市地區(qū)的氣溶膠光學厚度的平均值進行排序,得出氣溶膠光學厚度從大到小的分布情況為:2012年9月19日>2012年9月20日>2012年10月2日>2012年9月18日>2012年9月1日>2012年9月28日>2012年9月5日。結(jié)合湖北省環(huán)保廳查詢到的大氣污染狀況數(shù)據(jù),可得到2012年9月1日、2012年9月5日、2012年9月18日、2012年9月19日、2012年9月20日、2012年9月28日、2012年10月2日的API分別為:88、66、83、119、105、71、85。
以上述數(shù)據(jù)為基礎,得到AOD與API的分布趨勢如圖5所示。其中,橫軸代表時間,縱軸(左)代表API,縱軸(右)代表AOD。可以看出,氣溶膠光學厚度的反演結(jié)果與大氣污染狀況的時間分布趨勢基本一致,從時間序列分布狀況反映了反演結(jié)果的正確性。

圖5 氣溶膠光學厚度時間序列分布趨勢
3)與MODIS氣溶膠光學厚度產(chǎn)品驗證結(jié)果對比。根據(jù)國家衛(wèi)星氣象中心運行的MODIS氣溶膠遙感產(chǎn)品質(zhì)量檢驗分析結(jié)果,1 km分辨率MODIS L1B AOD產(chǎn)品的絕對誤差平均值為0.041 1,均方根誤差為0.27[16]。本文反演結(jié)果的絕對誤差平均值為0.051 1,精度略低于MODIS產(chǎn)品,均方根誤差為0.057 3,準確度高于MODIS產(chǎn)品。綜合考慮,反演結(jié)果與1 km分辨率MODIS氣溶膠產(chǎn)品具有相當?shù)馁|(zhì)量評價,但較高的時空分辨率決定其具有更高的細致性評價。
本文以武漢及周邊為研究區(qū)域,選取適用于干旱、半干旱以及城市等亮地表區(qū)域的結(jié)構(gòu)函數(shù)法,進行氣溶膠光學厚度反演,通過對結(jié)構(gòu)函數(shù)法的理論基礎和模型進行研究,確定了結(jié)構(gòu)函數(shù)公式、反演窗口和距離值。然后編程實現(xiàn)對武漢地區(qū)2012年9月至10月的8景HJ衛(wèi)星CCD影像的氣溶膠光學厚度反演。并結(jié)合CE318實測數(shù)據(jù)、湖北省環(huán)境保護廳網(wǎng)站提供的武漢市大氣污染狀況及分布趨勢數(shù)據(jù)、MODIS氣溶膠產(chǎn)品質(zhì)量檢驗分析結(jié)果,從多方面對反演結(jié)果進行了綜合驗證。HJ-1 A和HJ-1 B雙星聯(lián)合觀測情況下,全球覆蓋的時間分辨率是2 d,空間分辨率30 m,與MODIS和TM等載荷相比,具有較高的時空分辨率,適用于結(jié)構(gòu)復雜、氣溶膠厚度變化較快的城市地區(qū)業(yè)務化監(jiān)測,是一個較有意義的嘗試。且驗證結(jié)果表明,產(chǎn)品具有較高精度,與實際特征相一致,使用結(jié)構(gòu)函數(shù)法進行HJ-1 A/B衛(wèi)星的CCD影像氣溶膠光學厚度反演在武漢地區(qū)具有較高的可靠性,可應用于該地區(qū)的氣溶膠分布和變化趨勢分析,以及大范圍的區(qū)域大氣環(huán)境監(jiān)測等。