基于TM圖像的南京市氣溶膠光學厚度反演

程 晨，陳 健，李鑫慧

(南京信息工程大學遙感學院，南京 210044)

0 引言

大氣氣溶膠是大氣與懸浮在其中的固體和液體微粒共同組成的多相體系，是指懸浮在大氣中的具有一定穩定性、沉降速度小的，尺度范圍在10-3～102μm之間的分子團、液態或固態粒子所組成的混合物[1]。描述氣溶膠性質的基本參數有光學厚度、尺度分布、折射率、散射相函數。在這些參數中，氣溶膠光學厚度(aerosol optical depth，AOD)是各個特征參數的總體效應，集中體現了氣溶膠的基本狀況。傳統的地面實時觀測仍是研究AOD的重要手段[2-4]，但由于觀測站點數量有限，且分布比較離散，無法反映整個區域的氣溶膠空間分布情況。遙感技術連續、宏觀、動態、快速的特點，為大范圍的氣溶膠觀測和研究提供了可能。

20世紀90年代，Kaufman等利用綠色植被在紅光波段和藍光波段的反射率與2.1 μm處的反射率存在線性關系的這一性質，通過提取暗像元對陸地上的氣溶膠光學厚度進行了反演[5];隨后Robert等人針對“暗像元”的缺陷進行了改進，研究出了更加成熟的MODIS v5.2氣溶膠產品業務化算法[6]。在國內，李成才等通過分析MODIS光學厚度遙感產品討論了北京及周邊地區的大氣污染情況[7]。隨著城市大氣污染的監測需求日益增大，空間分辨率為30 m的Landsat TM數據在氣溶膠監測中得到廣泛應用。閔祥軍[8]、宋巍巍[9]等利用暗像元法與 DEM模型結合的方法分別反演了敦煌地區和廣州市的AOD，得到了較好的效果，但該方法對DEM的空間分辨率要求較高;劉小平[10]等通過改進暗像元獲得AOD信息，從而實現了TM圖像的快速大氣校正，但是對于暗像元的選取要求較為苛刻;王耀庭[11]等利用兩幅TM圖像對比的結果對北京市的氣溶膠空間分布進行了研究，但其前提是需要找到一幅“無污染”圖像作為基準，使該方法的使用受到了限制。

本文以南京市為研究區，以Landsat5 TM圖像為數據源，通過6S大氣輻射傳輸模型建立查找表，在此基礎上建立暗像元AOD的多元回歸模型，最后利用克里金插值得到南京市的氣溶膠光學厚度的空間分布數據，并對南京市AOD的空間分布特征進行了分析。

1 研究區概況和數據源

1.1 研究區概況

以整個南京市地區為研究區。南京市位于N31°14′～32°37′，E118°22′～119°14′之間。圖 1 為南京市的TM圖像。南京地區屬于亞熱帶季風氣候，雨量充沛，四季分明。全市擁有常住人口810萬，南京市是一座人類活動密集、偏重于工業生產的大城市，每天會排放出大量的廢氣，特別是機動車排放出的尾氣，礦產采集和化工生產的過程中揚起的粉塵都會進入大氣中，成為氣溶膠的一部分，嚴重影響當地的空氣質量，以致影響居民的正常生活。

圖1 研究區TM圖像Fig.1 TM image of the study area

1.2 遙感數據與預處理

本研究選用的是一幅2006年5月20日的Landsat5 TM圖像，軌道號120-38，從國際科學數據服務平臺上下載獲得。當天天氣晴朗，少云。首先，通過已知的控制點對圖像進行幾何糾正，重采樣方式為最鄰近插值法;然后按照南京市的行政區劃裁切出研究區范圍;最后進行輻射定標，將圖像中每個像元對應的DN值轉換為輻亮度，進行AOD反演，詳細計算過程和主要參數參見文獻[12]。

2 氣溶膠光學厚度反演原理與方法

在地表為朗伯體、大氣水平均一的假設條件下，大氣頂部的表觀反射率[5]可以表達為

式中:ρ*為表觀反射率(衛星觀測到的反射率);ρa為大氣的路徑輻射對表觀反射率的貢獻;μ，φ，μ0，φ0分別為太陽天頂角的余弦、方位角和觀測天頂角的余弦、方位角;T(μ0)和T(μ)分別為上行輻射和下行輻射的透過率;s為大氣的球面反照率;ρ為地面的真實反射率。在實際反演AOD時，往往利用輻射傳輸模型計算出 AOD 與 s，ρa，T(μ0)和 T(μ)等參數的對應關系，據此建立查找表，通過查找表獲取氣溶膠光學厚度。本文的技術路線如圖2所示。

圖2 反演氣溶膠光學厚度流程Fig.2 Flow chart of AOD retrieveing

2.1 計算表觀反射率和地表反射率

圖像在經過預處理后，通過計算將暗像元的表觀輻亮度轉化為表觀反射率，即

式中:ρ為表觀反射率;L為表觀輻亮度，W·m-2·μm-1sr-1;D為日地之間距離(天文單位);ESUN為大氣層頂的平均太陽輻照度，W·m-2;θ為太陽天頂角。

根據Kaufman的研究，在潔凈大氣的條件下，綠色植被在可見光的紅、藍波段和2.1 μm紅外波段的反射率在圖像上有相似的分布情況，進而推測紅、藍波段和2.1 μm的中紅外波段可能存在著一定的線性關系，而這種線性關系對于紅外波段波長為2.08 ～2.35 μm 的 TM 圖像同樣適用[13]。所以，如果TM圖像上植被覆蓋區在可見光的紅、藍波段的地表反射率可以通過紅外通道的反射率估算出來，那么圖像上得到的紅、藍波段的表觀反射率和真實地表反射率之間的差異就是氣溶膠造成的影響。TM圖像的紅光和藍光波段地表反射率由以下經驗公式[13]獲得，即

式中，Rred，Rblue分別為TM圖像紅、藍波段的地表反射率;MIR為TM圖像中紅外波段的表觀反射率。

2.2 確定暗像元

在紅、藍波段遙感圖像上，密集的植被相對于其他地物背景較暗，所以被稱為“暗像元”。本文采用閾值法，通過設置NDVI和MIR的閾值來提取暗像元[14]，當像元滿足 NDVI＞0.2，MIR ＜0.1 時，該像元被認定為暗像元。暗像元提取結果見圖3。由圖上可以看出，暗像元基本呈現均勻分布。

圖3 暗像元提取結果Fig.3 Results of dark pixel extraction

2.3 確定6S輻射傳輸模型及其參數

6S輻射傳輸模型是一個被廣泛接受和使用的大氣輻射傳輸模型，用來模擬地氣系統中太陽輻射的傳輸過程。其輸入參數主要包括太陽和衛星幾何路徑參數、大氣模式、氣溶膠模式、光譜條件和地面反射率類型。

TM圖像的太陽和衛星幾何路徑參數及光譜條件由6S模型提供;大氣模式選擇為中緯度夏季大氣(確定大氣溫、壓、濕廓線、臭氧、水汽含量);下墊面假定為均一表面，無方向性反射。而根據MODIS AOD業務化算法[15]，本文研究區域被劃分為中等吸收氣溶膠類型(單次散射反照)。陳驍強計算AERONET太湖站點2005年9月—2009年10月月平均單次散射反照率的平均值后發現，太湖站的單次散射反照率平均值為0.907 5，與大陸型氣溶膠中等吸收型氣溶膠對應的單次散射反照率值最為接近[16]，所以本文在6S輻射傳輸模型中氣溶膠類型選擇為大陸型氣溶膠。

2.4 構建查找表

本文借助6S輻射傳輸模型來模擬TM圖像紅、藍波段在一定的氣溶膠光學厚度和氣溶膠模式的條件下，表觀反射率隨AOD和地表反射率變化而發生的變化，從而建立查找表。建立查找表時，由于其他參數已經確定，主要考慮AOD和表觀反射率的變化。通過對TM1和TM3這2個波段的像元統計發現，超過99.9%的像元反射率分別集中于0.09～0.15和0.06 ～0.13 區間內，所以在查找表中 TM1和TM3波段的表觀反射率分別設定為0.09～0.15和0.06～0.13(都是以0.01為步長)。由于晴空下一般AOD＜2.0，且6S模型在氣溶膠光學厚度過大(能見度小于5.0 km)時會出現程序錯誤，所以TM1和TM3波段的氣溶膠光學厚度設定為0.21～1.21(以0.1為步長)，對于超出范圍的采用外推計算，建立起了利用TM數據反演AOD的查找表。

2.5 多元回歸建模

在得到TM1和TM3波段關于氣溶膠光學厚度、地表反射率、表觀反射率3者關系的查找表后，利用SPSS軟件分別對兩組數據進行多元線性回歸分析。TM1和TM3波段的回歸方程分別為

式中:z為氣溶膠光學厚度;x為表觀反射率;y為地表反射率。

2.6 反演氣溶膠光學厚度

通過查找表分別得到TM1和TM3波段在550 nm上的暗像元AOD值以后，求平均獲得南京市AOD空間分布數據。在ArcGIS中進行插值，為了減少計算量，在插值前將研究區TM圖像重采樣成60 m分辨率，重采樣方式為最鄰近像元法。由于在不均勻離散點分布的插值中，克里金插值比其他插值方法具有明顯的優越性[17]，所以本研究選用克里金插值法。得到的結果如圖4所示。

圖4 2006年5月20日南京市氣溶膠光學厚度空間分布Fig.4 AOD spatial distribution in Nanjing City，May 20，2006

3 結果與分析

3.1 AOD反演結果合理性檢驗

趙英時等認為，氣溶膠散射會對TM圖像的可見光波段產生較大的影響[18]，通過對比大氣校正前后南京市真彩色合成圖像和AOD反演結果，可以定性地檢驗AOD反演結果的合理性。因此，對TM數據進行了FLAASH大氣校正。在大氣校正前后的圖像上各選取茂密植被區和裸地區，獲取其隨波長變化而發生變化的光譜曲線，分析大氣校正前后地物光譜的變化。如圖5所示。

圖5 大氣校正前(左)后(右)植被和裸地光譜反射率曲線Fig.5 Spectral signature of bare soil，vegetation before(left)and after(right)atmospheric correction

從圖5可以看到，在可見光波段植被和裸地的表觀反射率在校正前后發生了明顯變化，其變化符合實際地物反射率變化規律，這說明大氣校正取得了較好的效果。

圖6分別為大氣校正前后的真彩色合成圖像以及把AOD反演結果取代校正后圖像紅波段進行彩色合成后的圖像。

圖6 南京市大氣校正前(左)、后(中)圖像與AOD反演結果(右)的比較Fig.6 Comparison of atmospheric correction before(left)and after(middle)and AOD inversion results(right)of Nanjing City

對比圖6(左)和圖6(中)可以看出，大氣校正后南京市城區大部分氣溶膠被去除，圖像對比明顯增強，地物的細節更加清晰(圖6(中));對比圖6(右)和圖6(左)則可以發現，AOD空間分布與圖6(左)中的氣溶膠(相對模糊的區域)的強弱分布基本一致。綜上所述，本研究中南京市當日的AOD反演結果較合理。

3.2 AOD空間分布特征分析

從圖4中可以發現，南京市AOD分布呈現北部相對較高，而中部和南部較低的基本空間分布特征，尤其在長江兩岸地勢平坦、人類生活頻繁、制造業發達的地區甚至出現了超過1.0的高值。在市區范圍，除了燕子磯、紫金山周圍茂密的植被覆蓋地區AOD較低以外，其他地區的AOD都在0.7～0.9之間，處于相對較高的水平，這很可能跟汽車尾氣排放等較密集的人類活動有關。從市區往南是江寧區、祿口鎮、溧水縣等區域，可以發現這一帶是整個南京市AOD整體最低的一個區域，主要原因是這些區縣相對市區人口較少，植被茂密。但是有小部分區域也出現了超過0.6的相對高值，如溧水縣愛景山地區。這是因為該地區有全國最大的鍶礦在開采，采礦過程中必然會產生很多粉塵，當粉塵散布到大氣中時就會對氣溶膠光學厚度產生顯著影響。根據上述描述可以作出初步判斷，植被、城市建成區、地形是影響AOD分布的主要因素。這個認識與李成才[19]、段婧[20]等的研究結果一致。

為了深入研究植被、城市建成區、地形對AOD空間分布的影響，本文引入了3個輔助數據，即歸一化植被指數(normalized difference vegetation index，NDVI)、歸一化城鎮指數(normalized difference building index，NDBI)[21]和數字高程模型 (digital elevation model，DEM)，將 AOD 以每 0.01 步長進行分割，然后分別計算每0.01步長中所有像元的NDVI，NDBI，DEM的平均值。計算結果如圖7所示。這里需要說明的是，在提取暗像元時，由于提取條件的限制，水體區并沒有暗像元提出，水體上空的AOD由周圍區域暗像元AOD插值得來，所以為了不影響準確性，在統計的過程中剔除了水體。

圖7 AOD與NDVI(左)、NDBI(中)、DEM(右)間的變化規則Fig.7 Scatter plot of AOD corresponding with the average of NDVI(left)，NDBI(middle)and DEM(right)

從圖上可以發現，AOD 與 NDVI，NDBI，DEM 都存在著良好的相關性，隨著AOD的逐漸升高，3者都呈現了非常明顯的變化。

從圖7(左)中可以看出，AOD與NDVI呈現明顯的負相關關系。隨著AOD的增加，NDVI值相應迅速降低，然后趨于平緩，最后出現小幅上升。AOD低值的區域主要對應了NDVI較高的茂密森林區域，這說明植被茂密的程度與氣溶膠的空間分布存在緊密的聯系。當AOD處于高值時，對應的NDVI值變化不大，甚至有小幅的上升，這有可能是因為在一些AOD高值區域，氣溶膠顆粒的來源不再跟下墊面類型有緊密的關系，而是主要來自于周圍區域的擴散，這也可能就是南京市城區一些植被覆蓋很好的公園、綠地AOD偏高的原因。

從圖7(中)中可以看出，AOD與NDBI呈現明顯的正相關關系。隨著AOD的增加，NDBI值迅速升高，在0.7～0.8處出現了震蕩，在趨于平緩后又出現小幅下降。這說明城鎮建筑的集中程度與氣溶膠的空間分布也存在緊密的聯系。在AOD為0.7～0.8時出現的震蕩，經過觀察發現六合區東南部AOD相對較高，而這部分區域下墊面主要為農田，NDBI處于相對低值。造成這樣情況的原因可能是這里地勢平坦，向西是處于AOD高值區的大廠區域，受到擴散的作用;而南部是長江，水體周圍的區域會受到水體蒸發產生液態懸浮性氣溶膠顆粒的影響。當AOD處于高值時，NDBI值變化不大，甚至有小幅的下降，這可能由于在南京城市規劃中，化工廠、發電廠等能夠產生大量大氣污染物的區域并不是位于市區中心區域，而是位于大廠鎮等城市邊緣鄉鎮。可見，造成一座城市AOD最高的區域常常不是人口和商業活動最集中的區域。

從圖7(右)中可以看出，隨著地形的升高，AOD呈現負指數衰減，這與徐希孺[22]等的觀點一致。圖7中AOD的低值對應的主要是南京市的一些地勢較高的區域，例如紫金山、老山等海拔超過200 m的區域，而AOD的高值區主要位于海拔50 m以下的區域，這主要是由于地球重力作用，氣溶膠顆粒密度隨高度呈指數衰減[22]。對于地形較平坦的南京來說，地形的微小變化無法對氣溶膠的擴散產生影響，但邊界層的湍流運動導致城區氣溶膠顆粒向四周大范圍的擴散，這可能就是江寧北部和六合南部地區雖然工業較少，依然很難看到藍天的緣故。

綜上所述，利用TM圖像獲取60 m空間分辨率的AOD空間分布信息，對于在城市尺度上AOD的定量化研究會有所幫助，并且對于類似南京市這樣的空間污染相對嚴重且分布具備明顯規律的城市非常有應用價值。

4 結論與展望

本文利用TM數據，對南京市2006年5月20日的氣溶膠光學厚度進行了衛星遙感反演，得到了當天南京市氣溶膠光學厚度的空間分布圖。主要獲得以下結論:

1)基于TM圖像，通過暗像元法和輻射傳輸模型建立查找表反演氣溶膠光學厚度的方法是可行的，達到了氣溶膠光學厚度精細化反演的目的。在提取足夠的濃密植被作為暗像元過程中，針對南京市的植被覆蓋情況和TM圖像的特點，NDVI＞0.2，MIR＜0.1的閾值設定是有效的。

2)南京市的AOD與植被覆蓋呈現明顯的負相關，與城市建筑密集程度呈現明顯的正相關，與地形高低呈現明顯的負指數相關。氣溶膠主要聚集于植被覆蓋較低、城市建筑密集、人類活動較頻繁、地形較低的區域。高空間分辨率圖像比一般分辨率圖像能夠提供的更多的氣溶膠空間分布信息。

結合輻射傳輸模型進行AOD反演的暗像元法針對植被覆蓋條件較好的中國東南地區是可行的，但是也正是需要濃密植被覆蓋這一條件局限了這一方法在其他植被覆蓋稀少地區的應用。如果能夠在以后的研究中減少對于濃密植被像元的依賴，并且更好地解決水體上空暗像元不足的問題，那么城市AOD反演的精細化研究一定可以得到進一步的發展。

[1]章澄昌，周文賢.大氣氣溶膠教程[M].北京:氣象出版社，1995:2.Zhang C C，Zhou W X.Atmospheric aerosol tutorial[M].Beijing:China Meteorological Press，1995:2.

[2]夏祥鰲，王普才，陳洪濱，等.中國北方地區春季氣溶膠光學特性地基遙感研究[J].遙感學報，2005，9(4):429-437.Xia X A，Wang P C，Chen H B，et al.Ground-based remote sensing of aerosol optical properties over North China in spring[J].Journal of Remote Sensing，2005，9(4):429-437.

[3]趙柏林，王 強，毛節泰，等.光學遙感大氣氣溶膠和水汽的研究[J].中國科學:B 輯，1983，13(10):951-962.Zhao B L，Wang Q，Mao J T，et al.Study on remote sensing of aerosol optical depth and vapor[J].Scientia Sinica Chimica:Serier B，1983，13(10):951-962.

[4]毛節泰，王 強，趙柏林.大氣透明度光譜和渾濁度的觀測[J].氣象學報，1983，41(3):322-331.Mao J T，Wang Q，Zhao B L.The observation of the atmospheric transparency spectrum and the turbidity[J].Acta Meteorologica Sinica，1983，41(4):322-331.

[5]Kaufman Y J，Tanré D，Remer L A，et al.Operational remote sensing of tropospheric aerosol over land from EOS moderate resolution imaging spectroradiometer[J].Journal of Geophysical Research，1997，102(D14):17051-17067.

[6]Levy R C，Remer L A，Matto S，et al.A new algorithm for retrieving aerosol properties over land from MODIS spectral reflectance[J].Journal of Geophysical Research，2006.

[7]李成才，毛節泰，劉啟漢，等.利用MODIS光學厚度遙感產品研究北京及周邊地區的大氣污染[J].大氣科學，2003，27(5):869-880.Li C C，Mao J T，Liu Q H，et al.Research on the air pollution in Beijing and its surroundings with MODIS AOD products[J].Chinese Journal of Atmospheric Sciences，2003，5(27):869-880.

[8]宋巍巍，管東生.利用TM影像反演廣州市氣溶膠光學厚度空間分布[J].環境科學學報，2008，28(8):1638-1645.Song W W，Guan D S.The distribution of aerosol optical depth retrieved by TM imagery over Guangzhou，China[J].Acta Scientiae Circumstantiae，28(8):1638-1645.

[9]閔祥軍，朱永豪.基于Landsat-TM圖像自身的反射率反演方法[J].遙感技術與應用，1997，12(3):2-9.Min X J，Zhu Y H.A method of reflectance retrieval based on image of Landsat-5 TM[J].Remote Sensing Technology and Application，1997，12(3):2-9.

[10]劉小平，鄧孺孺，彭曉鵑.基于TM影像的快速大氣校正方法[J].地理科學，2005，25(1):88-92.Liu X P，Deng R R，Peng X J.A fast atmospheric correction method based on TM imagery[J].Scientia Geographica Sinica，2005，25(1):88-92.

[11]王耀庭，王 橋，楊一鵬，等.利用Landsat/TM影像監測北京地區氣溶膠的空間分布[J].地理與地理信息科學，2005，21(3):20-22.Wang Y T，Wang Q，Yang Y P，et al.Monitoring spatial distribution of aerosols in Beijing based on Landsat TM Data[J].Geography and Geo-Information Science，2005，21(3):20-22.

[12]Chander G，Markham B.Revised Landsat 5 TM radiometric calibration procedures and post calibration dynamic ranges[J].IEEE Transactions on Geoscience and Remote Sensing，2003，41(11):2674-2677.

[13]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].IEEE Transactions on Geoscience and Remote Sensing，1997，35(5):1286-1298.

[14]Guo Y W，Tsay S C，Cahalan R F，et al.Path radiance technique for retrieving aerosol optical thickness over land[J].Journal of Geophysical Research，1999，104(D24):31321-31332.

[15]Remer L A，Tanré D，Kaufman Y J，et al.Algorithm for remote sensing of tropospheric aerosol from MODIS:Collection 5[R].National Aeronautics and Space Administration，2006.

[16]陳驍強.基于環境一號衛星陸地氣溶膠光學厚度反演方法研究[D].南京:南京師范大學，2011.Chen X Q.Study on the retrieval of aerosol optical thickness over land based on HJ-1 data[D].Nanjing:Nanjing Normal University，2011.

[17]劉光孟，汪云甲，張海榮，等.空間分析中幾種插值方法的比較研究[J].地理信息世界，2011，9(3):41-45.Liu G M，Wang Y J，Zhang H R，et al.Comparative study of several interpolation methods on spatial analysis[J].Geomatics World，2011，9(3):41-45.

[18]趙英時.遙感應用分析原理與方法[M].北京:科學出版社，2003:21-30.Zhao Y S.The principle and method of analysis of remote sensing application[M].Beijing:Science Press，2003:21-30.

[19]李成才，毛節泰，劉啟漢.利用MODIS資料遙感香港地區高分辨率氣溶膠光學厚度[J].大氣科學，2005，29(3):335-342.Li C C，Mao J T，Liu Q H.Remote sensing of high spatial resolution aerosol optical depth with MODIS data over Hong Kong[J].Chinese Journal of Atmospheric Sciences，2005，29(3):335-342.

[20]段 婧，毛節泰.長江三角洲大氣氣溶膠光學厚度分布和變化趨勢研究[J].環境科學學報，2007，27(4):537-543.Duan J，Mao J T.Study on the distribution and variation trends of atmospheric aerosol optical depth over the Yangtze River Delta[J].Acta Scientiae Circumstantiae，2007，27(4):537-543.

[21]查 勇，倪紹祥，楊 山.一種利用TM圖像自動提取城鎮用地信息的有效方法[J].遙感學報，2003，7(1):38-39.Cha Y，Ni S X，Yang S.An effective approach to automatically extract urban Land-use from TM imagery[J].Journal of Remote Sensing，2003，7(1):38-39.

[22]徐希孺.遙感物理[M].北京:北京大學出版社，2003:198-312.Xu X R.Physics of remote sensing[M].Beijing:Peking University Press，2003:198-312.