郭 強,唐 家 奎*,何 文 通,田 媛,于 新 菊,2
(1.中國科學院大學,北京 100049;2.中國科學院遙感與數字地球研究所,北京 100101)
?
利用MODIS可見光波段反演陸地氣溶膠光學厚度
郭 強1,唐 家 奎1*,何 文 通1,田 媛1,于 新 菊1,2
(1.中國科學院大學,北京 100049;2.中國科學院遙感與數字地球研究所,北京 100101)
利用Herold等建立的地表反射率庫及MODIS遙感影像研究城市區和非城市區典型地物在可見光紅藍波段地表反射率的比值特性。在此基礎上,利用MODIS 1 km分辨率遙感影像紅藍可見光波段實現了氣溶膠光學厚度的反演,采用衛星過境時間前后半小時北京和香河AERONET站的氣溶膠光學厚度觀測平均值作為驗證參考。結果顯示,66.67%的反演結果處于±0.05±0.15 τ的誤差界限內,反演算法不受地表反射率的限制,而且只利用了可見光紅藍波段,避免缺少近紅外波段數據的限制。
氣溶膠光學厚度;遙感;MODIS;地表反射率
大氣氣溶膠作為大氣中最活躍的3種因素之一,它的輻射強迫效應深刻影響著全球的輻射能量平衡及氣候變化[1]。氣溶膠顆粒在大氣中的停留時間,一般為幾天到幾周,而且其特征隨空間和時間都有顯著的變化。目前研究確定大氣氣溶膠對氣候系統的精確輻射強迫作用仍十分困難,而準確的氣溶膠監測對于研究全球輻射能量平衡及環境監測保護具有重要意義。
衛星遙感具有大尺度同步觀測、高時效性等特點,使之成為監測大氣氣溶膠的有效手段。衛星傳感器接收的信號來源于地表和大氣兩部分。利用衛星影像反演氣溶膠特性首先需要從衛星接收信號中區分地表貢獻和大氣貢獻,即實現地氣解耦。不同衛星傳感器利用衛星觀測值的光譜、角度及極化等特性實現地氣解耦。Kaufman等[2,3]基于地表反射率特征提出暗目標法,其基本原理是在濃密植被區及暗色土壤區,可見光的紅、藍波段地表反射率較低且與短波紅外波段存在線性相關關系。在短波紅外波段氣溶膠的散射吸收作用很小,假設氣溶膠光學厚度為0,衛星觀測到的大氣層頂表觀反射率就是地表真實反射率。通過短波紅外波段的反射率推斷可見光波段的地表反射率,實現地氣解耦。該方法局限于暗目標區域,對于城市、干旱和半干旱等高地表反射率區域無法獲得氣溶膠光學厚度(Aerosol Optical Depth,AOD)。此外,算法需要短波紅外通道,而目前在軌的遙感衛星大多沒有短波紅外通道,尤其是高分衛星,限制了該算法的應用。Tanré等[4,5]假定一定范圍內復雜地表上空大氣透過率相同,提出了結構函數法以解決高地表反射率對于氣溶膠監測的限制。結構函數法較少受到地表反射率大小的限制,但由于對幾何校正的要求較高,業務化運行比較困難。Deuzé等[6]從POLDER衛星觀測獲得的大氣后向散射偏振信息中分離出氣溶膠貢獻,進而反演得到氣溶膠光學厚度,但偏振的方法只能反演細粒子氣溶膠。Tang等[7]基于AQUA和TERRA數據提出了雙星協同反演算法,實現了氣溶膠光學厚度和地表反射率的同時反演,部分解決了亮地表區域反演無效的難題,但使用兩顆衛星的數據反演的光學厚度會受到傳感器的性能退化、不同傳感器影像之間圖像匹配精度等影響。Knapp等[8]假設地表反射率特征在一段時間內保持不變,而且至少有一天是不受氣溶膠影響的清潔日。考慮到殘云的影響,選取每個像元位置第二最低的反射率值合成一幅參考影像。假設一個背景AOD并用其校正參考影像,得到參考的地表反射率用于反演AOD,算法的誤差主要來源于地表反射率和背景AOD的假設。Hsu等[9,10]根據在紅、藍波段AOD對天頂輻亮度有顯著的貢獻提出了基于地表反射率庫的深藍(Deep Blue)算法,并成功應用于撒哈拉沙漠、阿拉伯半島等干旱半干旱地區,但該方法需要事先建立地表反射率庫作為先驗知識。Diner等[11]假設地表輻射亮度的角度分布與波長無關,利用多角度觀測數據消除地表反射率的影響,反演氣溶膠光學厚度。Wong等[12]利用模型模擬將衛星觀測到的表觀反射率解析為氣溶膠反射率和大氣瑞利散射引起的路徑反射率,在此基礎上利用最小反射率技術(MRT)確定季節性的地表反射率,然后利用SBDART輻射傳輸模型構建查找表LUT,反演得到香港地區500 m分辨率的AOD。該反演算法有希望應用于其他區域的城市/局部尺度高分辨率AOD反演。Remer等[13]提出了3 km分辨率的AOD反演算法,并將加入MODIS C006版本的AOD反演算法中。MODIS 3 km反演算法相比10 km反演算法縮小了反演窗口,精度也略差于10 km分辨率算法。但提高分辨率的AOD產品可以為相關領域的研究提供重要的工具。
本文基于MODIS數據探討北京及周邊不同地表覆蓋類型區域的可見光紅藍波段地表反射率比值特征,并嘗試利用可見光紅藍波段衛星觀測值反演氣溶膠光學厚度,從而擺脫反演算法對近紅外波段的依賴。
1.1 MODIS數據
采用的數據包括衛星數據和地面站點數據。衛星數據包括MODIS Level 1B地表反射率數據、Level 2氣溶膠光學厚度數據、Level 3合成的地表反射率數據和Level 3土地覆蓋類型產品。MODIS是搭載于EOS/TERRA和EOS/AQUA太陽同步軌道系列衛星的主要傳感器。TERRA的過境時間為上午10∶30,AQUA的過境時間為下午13∶30。MODIS具有從可見光到遠紅外的36個光譜通道,地面分辨率分別為250 m、500 m和1 000 m,掃描幅寬2 330 km。反演采用1 km分辨率的地表反射率數據,時間范圍為2009年4-9月。Kaufman等[2]指出,遙感影像的空間分辨率為250 m及更高時,大氣點擴散函數引起的臨近像元影像顯著。當空間分辨率為1 km甚至更低時,臨近像元效應很小,在大氣校正過程中臨近像元效應帶來的誤差遠小于氣溶膠特征的不確定性帶來的誤差。本文選取1 km分辨率的數據,因此可以不考慮臨近像元效應的影響。
1.2 AERONET數據
地面站點數據來源于AERONET站點AOD數據。AERONET是NASA建立的覆蓋全球的地基氣溶膠監測網絡。站點觀測設備主要采用法國CIMEL公司的偏振太陽光度計CE318-II和標準太陽光度計CE318-I,觀測波段主要有8個:0.340 μm、0.380 μm、0.440 μm、0.500 μm、0.675 μm、0.870 μm、0.936 μm、1.020 μm。AERONET氣溶膠光學厚度數據分為3個等級:1級數據未做云去除;1.5級數據僅做了云去除;2級數據既經過云去除又經過人工檢查。本文采用其官網(http://aeronet.gsfc.nasa.gov)發布的2級AOD數據對衛星反演的氣溶膠光學厚度進行驗證。AERONET 2級AOD數據的精度為±0.02[14],因此可以認為是AOD的真值。AERONET站點在0.55 μm波長沒有AOD觀測值,本文采用Angstrom公式擬合得到。
τλ=β·λ-α
(1)
式中:β為大氣渾濁度指數,α為Angstrom指數,λ為波長,τλ為λ對應的氣溶膠光學厚度。
采用北京站(位于城市區域)、香河站(位于非城市區域)2個AERONET站點驗證反演的結果。
1.3 反演原理與算法
Flowerdew等[15]研究指出,地物對于不同波長的反射特征取決于地物的微觀結構,即地物的分子、原子與電磁波的作用機理。因此對于不同的波長,地物BRDF的大小不同,但地物反射率的角度分布特征即BRDF的形狀特征主要取決于地物的宏觀結構,與波長的相關性很小,幾乎可以忽略。因此有:
(2)
式中:rm、rn分別表示兩個不同觀測幾何條件的地表反射率,m、n分別表示兩個不同的觀測幾何條件;b1、b2表示兩個不同的波段。
由上式可以導出:
(3)
式中:比例系數k為常數。
通過以上分析可知,對于同一種地物,在兩個觀測幾何條件下兩個波段反射率的比值相同。Levy等[16]通過對全球AERONET站點地表反射率的研究發現,在暗像元區可見光紅藍波段反射率的相關性比可見光藍波段與短波紅外波段的相關性高,而且更加穩定。在MODIS 官方C005及以后版本的反演算法中用以下兩個等式實現地氣解耦:
(4)

Levy等[16]對AERONET站點的研究發現,全球范圍內紅藍波段地表反射率比值恒定在1.923附近,夏季城市區比較低(1.305左右)。因此考慮利用這個特性實現AOD的反演。
根據Vermote等[17]的理論,在地表朗伯體、大氣水平均一的假設條件下,衛星平臺觀測到的大氣頂部反射率可以表示為:

(5)
式中:ρ0是大氣路徑輻射項等效反射率,ρs是地表反射率,T(μs)、T(μv)分別是入射方向和觀測方向的大氣透過率,S為大氣下界面的半球反射率。T(μs)、T(μv)總是同時出現,因此將T(μs)T(μv)作為一個參數考慮;未知數ρ0、S、T(μs)T(μv)與AOD、氣溶膠粒子單次散射率、相函數、不對稱指數等有關。
在實際反演過程中,首先針對試驗區域確定一個氣溶膠模型,然后用6s、MODTRAN等輻射傳輸模型模擬在不同觀測幾何情況下AOD與上述3個參數之間的對應關系,據此建立查找表得到AOD。MODIS紅、藍波段地表反射率之間存在如下關系:
(6)

(7)
不同區域地表覆蓋類型不同,假設可見光紅波段與藍波段地表反射率比值在城市區、植被區固定。采取3種方法驗證:1)根據Herold等[18]提出的地表反射率庫,選取包括植被區、高亮地表、城市區等典型地物,用MODIS波段響應函數擬合地表反射率,考察紅、藍波段的反射率比值;2)選取一系列無云的MODIS影像,提取香河和北京站點位置的表觀反射率,根據AERONET的觀測值插值得到0.55μmAOD,然后用6s校正,得到香河和北京站點的地表反射率,最后考察反射率比值;3)選取春、夏、秋、冬4幅無云影像的MODIS地表反射率產品MOD09,考察地表反射率比值。具體反演流程見圖1。
2.1 紅、藍波段地表反射率比值估計

圖1 算法反演流程
Fig.1Flowchartofalgorithmprocessing
Herold等[18]利用AVIRIS數據研究了城市地區典型地物的光譜特征,并建立了城市地區包括植被、土壤、瀝青道路、水泥道路、屋頂等多種典型地物的光譜反射率數據庫。AVIRIS數據具有0.4~2.5 μm范圍內224個光譜波段,空間分辨率為4 m,光譜分辨率為10 nm。從上述光譜數據庫中選取了9種典型地物類型,用MODIS紅藍波段響應函數擬合得到反射率,計算紅、藍波段地表反射率比值與藍波段地表反射率的關系,結果如圖2所示(見封3)。
從圖2中藍色橢圓內的地物可以看出,裸土、綠色植被、不進行光合作用的枯萎植被大多集中在1.4~2.1的范圍內。城市周邊主要地物是植被和高亮的裸土,因此可以推斷在城市周邊區域像元的紅、藍波段反射率比值處于1.4~2.1。城市區域的典型地物選取了3種典型的屋頂、3種道路。城市地區地物在藍波段的地表反射率分散很大,例如水泥道路的反射率分散在0.13~0.46的范圍內。紅、藍波段地表反射率比值分散也較大,分布在1.0~3.2之間。城市地物反射率比值的差異性反映了城市區地物的復雜性。由于選取的MODIS影像分辨率為1 km,在城市區多為混合像元,且受到臨近像元的影響,因此需針對MODIS數據特點考察不同地物地表反射率特征。
為了考察MODIS影像上地物地表反射率特征,選取2009年北京周邊61景無云的MODIS 1 km分辨率影像。為了降低地表BRDF的影響,選取的影像在AERONET北京站和香河站的太陽天頂角需小于30°。對于每一景影像,提取AERONET香河站和北京站的表觀反射率,用衛星過境時間AERONET站點觀測的氣溶膠光學厚度進行大氣校正,得到紅、藍波段地表反射率,進而得到地表反射率的比值(圖3),可以看出,AERONET站點多數比值處在1.65~2.25之間。

圖3 AERONET站點紅藍波段地表反射率的比值
Fig.3 The ratio of red/blue band at the AERONET sites
為了進一步考察研究區域地表反射率比值特性,選取了四季的MOD09地表反射率產品,共4景影像。影像的成像時間分別為春分、夏至、秋分和冬至日,空間范圍為北京周邊400×400 km。用MODIS Level 3土地覆蓋類型數據MOD12將研究區分為城市區和非城市區,分別考察不同區域地表反射率特征。對于每種區域,在影像上隨機選取2 000個樣點,考察地表反射率比值,結果如圖4所示。
圖4a-圖4d分別為城市區四季紅藍波段地表反射率,圖4e-圖4h分別為非城市區四季紅藍波段地表反射率??梢钥吹剑鞘袇^和非城市區可見光紅藍波段地表反射率都具有很強的相關性,相關系數都在0.86以上。圖4中的直線是回歸線,顏色條顯示的是散點的密度??梢钥吹浇^大多數點分布在回歸直線兩側很小的范圍,顯示了回歸關系的穩定性。此外,除了春季非城市區相關性都大于城市區,這與城市區地表覆蓋復雜有關。四季城市區和非城市區的回歸直線的斜率如表1。

圖4 MOD04紅、藍波段地表反射率關系
Fig.4 The relationship of red/blue band from MOD04
表1 紅藍波段地表反射率回歸直線的斜率
Table 1 Slope of the regression line between red and blue band reflectance

春夏秋冬城市區1.93741.90921.94581.9105非城市區1.98871.92011.95571.9520
2.2 氣溶膠光學厚度反演實例
為了驗證算法的可行性,選取2009年北京周邊4-9月AERONET香河站和北京站上空無云的53景MODIS 1 km影像進行反演實驗。2009年5月3日和6月20日兩景遙感影像中研究區域上空云量很少。圖5(見封3)是兩個日期研究區域影像真彩色合成。
反演前需對遙感影像進行輻射定標、重投影、云去除、裁剪等操作。云像元的去除采用紅波段閾值法。Ackerman等[19]通過研究發現,在MODIS影像紅波段云像元的表觀反射率一般大于0.18,并將該值作為MODIS云檢測產品MOD06的陸地區域默認閾值。因此采用0.18作為云像元去除的閾值。查找表是基于6s輻射傳輸模型設定不同的觀測幾何條件計算得到的。6s輻射傳輸模型可以模擬0.25~4.0 μm波長范圍內衛星接收的信號。它采用連續階散射(Successive Orders of Scattering)方法求解輻射傳輸方程,考慮了吸收氣體(水汽、二氧化碳、甲烷、臭氧、二氧化氮等)的吸收、散射及氣溶膠粒子的多次散射。設定條件包括幾何條件、光譜波段、大氣模式、氣溶膠模式。選擇數據的時間都是4-9月之間,因此大氣模式選取中緯度夏季,氣溶膠模式選擇大陸性氣溶膠。
圖6(見封3)是兩個日期反演得到的1 km分辨率氣溶膠光學厚度??梢钥闯觯瑑蓚€日期北京城區上空的氣溶膠光學厚度都明顯大于周邊地區,另外,天津也是氣溶膠光學厚度的高值帶。從氣溶膠光學厚度的區域分布看,兩個日期都呈現西北低、東南高的分布趨勢。北京西北部山區植被茂密、人口稀少,氣溶膠光學厚度小,大氣污染相對較輕;北京城區及其東南方向平原地區人口密集,氣溶膠光學厚度相對較大,反映了這些地區大氣污染嚴重。6月20日北京城區東北部的山區上空有一自東向西的氣溶膠光學厚度高值分布帶,可能是由周邊高污染地區的污染物通過風力的作用漂移、聚集形成。研究區域地表覆蓋多樣,包含了城市、濃密植被、高反射裸土等地物,反演結果較好地反映了不同類型區域上空氣溶膠粒子的空間分布特征。
為了驗證反演結果,獲取兩景影像衛星過境時間北京和香河AERONET站點觀測的氣溶膠光學厚度,與反演結果對比,結果如表2所示。5月3日的反演結果與AERONET觀測值很接近,兩個站點最大的誤差為0.088。6月20日的反演結果與AERONET結果有一定的偏差,兩個站點最大的誤差達0.173。
表2 反演結果與AERONET觀測結果對比
Table 2 Retrieved AOD and observed value from AERONET sites

北京站香河站反演AERONET誤差反演AERONET誤差2009-05-030.2370.2630.0260.2020.2900.0882009-06-200.3600.4980.1380.2560.4290.173
為了進一步驗證算法的精度,選取2009年4-9月香河和北京AERONET站點上空無云的75景影像進行了反演試驗。同時選取衛星過境時間前后半小時之內AERONET站點的觀測值取平均,得到光學厚度值。兩者的對比如圖7(見封3)所示。
圖7中橙色實線為回歸線,黃色虛線為誤差界限±0.05±0.15τ,灰色虛線為1∶1線??偟膩碚f,算法反演結果與AERONET觀測結果有較好的相關性,75個有效反演值中,有66.67%的反演結果處于±0.05±0.15 τ的誤差界限內。圖8為反演值和AERONET站點觀測值的時間序列散點圖,可以看到,AERONET北京站和香河站上空反演結果與觀測值的時間序列比較吻合。

圖8 反演結果與AERONET觀測值對比散點圖
Fig.8 Scatter plot of AOD retrieved and observed value from AERONET
本文提出了一種利用MODIS影像可見光波段數據反演高分辨率氣溶膠光學厚度的算法,利用地表反射率庫與MODIS遙感影像分析了北京周邊不同典型地物類型可見光紅、藍波段地表反射率的比值特征,實現了氣溶膠光學厚度的反演。選取了2009年4-9月北京和香河AERONET站點上空無云的多景遙感影像進行了反演實驗,并將反演結果與AERONET站點觀測值進行了對比驗證。結果顯示,有66.67%的反演結果處于±0.05±0.15 τ的誤差界限內。
本文反演算法不受地表反射率的限制,可應用于包含城市區、高亮裸土等典型地物的廣大區域。另外,反演算法只利用了可見光紅藍波段,有潛力應用于包含可見光紅藍波段的高分辨率遙感衛星數據。但需要指出的是,算法的反演需要提前估計大氣模型和氣溶膠模型,這是算法的誤差來源之一。
[1] SOLOMON S,QIN D,MANNING M,et al.Climate Change 2007——The Physical Science Basis.Contribution of Working Group I to the Fourth Assessment Report of the IPCC[M].London:Cambridge University Press,2007.
[2] KAUFMAN Y J,WALD A E,REMER L A,et al.The MODIS 2.1 μm channel-correlation with visible reflectance for use in remote sensing of aerosol[J].Geoscience and Remote Sensing,1997,35(5):1286-1298.
[3] KAUFMAN Y J,TANRé D,GORDON H R,et al.Passive remote sensing of tropospheric aerosol and atmospheric correction for the aerosol effect[J].J.Geophys.Res.,1997,102(D14):16815-16830.
[4] TANRé D,DESCHAMPS P Y,DEVAUX C,et al.Estimation of Saharan aerosol optical thickness from blurring effects in thematic mapper data[J].J.Geophys.Res,1988,93(D12):15955-15964.
[5] HOLBEN B N,VERMOTE E,KAUFMAN Y J,et al.Aerosol retrieval over land from AVHRR data——Application for atmosphere correction[J].IEEE Transactions on Geoscience and Remote Sensing,1992,30(2):212-222.
[6] DEUZé J L,BRéON F M,DEVAUX C,et al.Remote sensing of aerosols over land surfaces from POLDER/ADEOS-1 polarized measurments[J].Journal of Geophysical Research,2001,106(D5):4913-4926.
[7] TANG J,XUE Y,YU T,et al.Aerosol optical thickness determination by exploiting the synergy of TERRA and AQUA MODIS[J].Remote Sens.Environ.,2005,94(3):327-334.
[8] KNAPP K R,FROUIN R,KONDRAGUNTA S,et al.Toward aerosol optical depth retrievals over land from GOES visible radiances:Determining surface reflectance[J].International Journal of Remote Sensing,2005,26(18):4097-4116.
[9] HSU N C,TSAY S C,KING M D,et al.Aerosol properties over bright-reflecting source regions[J].IEEE Trans.Geosci.Remote Sens.,2004,42(3):557-569.
[10] HSU N C,TSAY S C,KING M D,et al.Deep blue retrievals of Asian aerosol properties during ACE-Asia[J].IEEE Trans.Geosci.Remote Sens.,2006,44(11):3180-3195.
[11] DINER D J,MARTONCHIK J V,KAHN R A.Using angular and spectral shape similarity constrains to improve MISR aerosol and surface retrievals over land[J].Remote Sensing of Environment,2005,94(2):155-171.
[12] WONG M S,NICHOL J E,LEE K H.An operational MODIS aerosol retrieval algorithm at high spatial resolution,and its application over a complex urban region[J].Atmospheric Research,2011,99(3):579-589.
[13] REMER L A,MATTOO S,LEVY R C,et al.MODIS 3 km aerosol product:Algorithm and global perspective[J].Atmos.Meas.Tech.Discuss.,2013,6(1):69-112.
[14] HOLBEN B N,ECK T F,SLUTSKER I,et al.AERONET——A federated instrument network and data archive for aerosol characterization[J].Remote Sens.Environ.,1998,66(1):1-16.
[15] FLOWERDEW R J,HAIGH J D.An approximation to improve accuracy in the derivation of surface reectances from multi-look satellite radiometers[J].Geophys.Research Letters,1995,22(13):1693-1696.
[16] LEVY R,REMER L,TANRE D,et al.Algorithm for remote sensing of tropospheric aerosol over dark targets from MODIS:Collections 005 and 051:Revision 2[DB/OL].MODIS Algorithm Theoretical Basis Document,2009.
[17] VERMOTE E,TANRE D,DEUZJ L,et al.Second simulation of the satellite signal in the solar spectrum:An overview[J].IEEE Trans.Geosci.Remote Sensing,1997,35(3):675-686.
[18] HEROLD M,GARDNER M,ROBERTS D A.Spectral resolution requirements for mapping urban areas[J].IEEE Transactions on Geoscience and Remote Sensing,2003,41(9):1907-1919.
[19] ACKERMAN S A,STRABALA K I,MENZEL W P,et al.Discriminating clear-sky from clouds with MODIS[J].J.Geophys.Res,1998,103(D24):32141-32157.
Aerosol Optical Depth Retrieval over Land Using MODIS Visible Bands Imagery
GUO Qiang1,TANG Jia-kui1,HE Wen-tong1,TIAN Yuan1,YU Xin-ju1,2
(1.UniversityofChineseAcademyofSciences,Beijing100049; 2.InstituteofRemoteSensingandDigitalEarthofChineseAcademyofSciences,Beijing100101,China)
A surface reflectance library proposed by Herold et.al.(2003)and MODIS imagery are used to estimate the ratio of surface reflectance between visible red and blue bands in both urban and non-urban areas.On this basis,blue and red bands of MODIS 1 km resolution imagery are used to retrieve aerosol optical depth.The observations of Beijing and Xianghe AERONET sites are used to validate the results.In order to minimize the error,the averaged measurements within half an hour before and after the MODIS overpass are used in the validation.The result shows that 66.67% of the retrievals are within the expected error envelope of ±0.05±0.15 τ compared to the AERONET observed AOD.The application of this algorithm is not subject to any limitations of surface reflectance.Moreover,the algorithm only utilizes the visible blue and red bands,without using NIR band.
aerosol optical depth;remote sensing;MODIS;surface reflectance
2014-04-17;
2014-09-08
國家973支撐項目(2013CB733402);中國科學院數字地球重點實驗室開放基金項目(2011LDE015);環保公益行業科研專項(201309011);中國科學院大學校長基金項目;浙江省重大科技專項(2012C13011-2)
郭強(1989-),男,碩士研究生,從事大氣遙感研究。 *通訊作者E-mail:jktang@ucas.ac.cn
10.3969/j.issn.1672-0504.2015.02.009
P407
A
1672-0504(2015)02-0038-06