大氣氣溶膠光學厚度反演軟件系統設計和實現

劉翔龍，楊再昕

（河海大學 計算機與信息學院，江蘇 南京 211100）

大氣氣溶膠是由大氣介質和混合于其中的固體和液體顆粒物等組成的體系，大氣氣溶膠主要通過直接或間接改變地氣系統的輻射收支來影響氣候和環境。由于大氣氣溶膠光學厚度AOD（Aerosol Optical Depth）是反映大氣渾濁度和確定氣溶膠氣候效應的重要依據，所以對氣溶膠光學厚度的研究具有很現實的意義。在氣溶膠光學厚度反演這一研究課題上國內外專家和學者已經做了很多工作。當前的氣溶膠光學厚度主要是以地基監測為主，但是地基監測方法會受到很多方面因素的影響，例如氣溶膠的生命周期較短，濃度的空間變化很大等。這樣只能監測整層大氣柱的氣溶膠的總量和單點的信息，同時由于各地區地理環境條件對儀器監測有所限制，所以只能在有限的區域進行，所設的網點分布不能全部覆蓋和反映各地信息。衛星遙感反演方法與傳統的地基探測方法不同，它具有覆蓋面積廣、信息獲取方便等特點，能更有效地獲取氣溶膠的信息，擺脫了地基探測局限性的缺點，使人們能夠實時掌握大面積范圍內的氣溶膠的變化提供了可能。目前，在使用衛星遙感資料進行氣溶膠光學厚度的反演是借助不同軟件的部分功能模塊加以整合得以實現的，而專門用來實現氣溶膠光學厚度反演的一體化的軟件幾乎沒有，這樣就使研究者需要對反演方法的全面掌握和對不同軟件的各部分功能模塊有所了解，增加了不必要的工作，使反演過程復雜不便[1-2]。

為了實現對反演軟件的開發，根據目前已有反演的理論基礎，主要提出和設計利用GDAL開源庫對MODIS L1B遙感數據提取和進行相關處理，選取暗像元后利用其參數信息結合6S模型建立查找表。使用QT設計用戶操作界面，實現反演可視化。本軟件能夠實現氣溶膠光學厚度的反演的功能。

1 反演原理

一般來說，當陸地表面是均勻朗伯表面，大氣垂直均勻變化時，衛星測量值可用表觀反射率ρ*表示：

式中L是衛星測量輻亮度，Es是大氣頂的太陽輻射通量，τ0是整層大氣光學厚度，（μs，Φs）表示太陽入射光的方向，（μv，Φv）表示衛星觀測方向，μs，Φs，μv，Φv分別為太陽入射光和衛星觀測天頂角的余弦和方位角。如果不考慮氣體吸收，衛星觀測的表觀反射率為：

式中θs為太陽天頂角，θv為觀測天頂角，Φ為相對方位角，Φs分別為太陽到地面和地面到衛星的整層大氣 T（θv）透過率，S 為大氣的球面反照率，ρ為地表反射率，Fd（θs）為表觀反射率為零時歸一化總的向下輻射通量。當地表反射率很小時，衛星觀測反射率主要取決于大氣貢獻項。當地表反射率很大時，地表的貢獻成為主要貢獻項。這樣可以得出反演的基礎：如果已知地表反射率，確定大氣氣溶膠模型就可以得到氣溶膠光學厚度。反之，如果已知氣溶膠光學厚度和相應的大氣參數，也可以得到地表反射率。

該軟件的設計主要依據暗像元方法DDV（Dense Dark Vegetation），它是由Kaufman和Sendra反演稠密植被上空氣溶膠光學厚度建立的。這種算法主要依據濃密植被和湖泊的反射率很低（約0.01～0.02），因此可以將有著大面積森林或水域的地方可以作為暗像元，暗像元算法是基于表觀反射率的大氣貢獻項，即利用衛星觀測的路徑輻射反演氣溶膠光學厚度，是氣溶膠遙感應用比較常用的算法。利用MODIS圖像反演氣溶膠光學厚度主要是利用了空間分辨率為250 m的1波段（620～670 nm）、500 m 的 3 波段（459～479 nm）、500 m 的7波段（2 105～2 155 nm）3個波段。通過大量實驗證明：中紅外通道（2 100 nm）通道表觀反射率除了受塵粒氣溶膠的影響外，幾乎不受其他氣溶膠的影響而且其值接近地表反射率。因此暗像元的確定主要是在幾何校正后通過第7波段在MODIS圖像中尋找表觀反射率小于0.05的點作為暗像元。

Kaufman等人通過研究大量的資料，得到紅（670 nm）、藍（470 nm）和中紅外通道（2 100 nm）的地表反射率存在如下關系：

利用中紅外通道的表現反射率和地表反射率近似相等的前提，可以采用7波段的表觀反射率代替地表反射率，代入式子（3）、（4）估算求出1、3波段的地表反射率。并且選取合理的氣溶膠模型，通過衛星觀測的表現反射率根據公式（2）計算可得出氣溶膠光學厚度[3-4]。

通過暗像元反演流程圖如圖1所示。

圖1 反演流程圖Fig.1 Flow chart of the inversion experiment

2 軟件設計

2.1 軟件設計的問題

該氣溶膠光學厚度反演軟件是基于windows 7系統操作環境下，編程采用了C++語言，開發工具選用QT C++，調用GDAL開源庫、HDF驅動進行讀取、處理HDF格式文件。系統結構功能圖如圖2所示。

圖2 系統結構功能圖Fig.2 Structure diagram of the system’function

對MODIS遙感圖像數據需要解決如下問題：

1）MODIS遙感影像數據采用HDF文件格式儲存。利用開源GDAL庫的API函數提取HDF文件的柵格數據及屬性實現對MODIS數據的讀取和寫入。

2）對MODIS數據進行幾何糾正，其中包括一般的遙感圖像校正和MODIS數據特有的“Bowtie”（蝴蝶結效應糾正）[5]。同時能夠實現對遙感圖像的一般操作，如圖像的放大，縮小以及漫游等功能。

3）軟件系統集成嵌入6S模塊，實現利用6S大氣輻射傳輸模式，建立起表觀反射率—地表反射率—氣溶膠光學厚度的查找表。

2.2 軟件設計的方法

GDAL（Geospatial Data Abstraction Library）是一個在 X/MIT許可協議下的開源柵格空間數據轉換庫。它主要可以對柵格數據格式進行相關數據處理，它使用抽象數據模型來分析數據格式，抽象的數據模型包括元數據、柵格波段、子數據集域圖、圖像結構域、仿射地理坐標轉換等。MODIS遙感數據采用HDF格式進行存儲，這種文件格式可以存儲不同類型的圖像和數碼數據。HDF包含有6種主要數據類型：柵格圖象，調色板，科學數據庫，注釋，Vdata和Vgroup。在GDAL中每個文件有一個數據集，HDF文件中包含多個數據集，可以將這些數據集提取到抽象模型中的子數據域，將屬性數據提取到子數據域里的元數據中[6]。利用GDAL對HDF的數據處理可以概括為：首先打開文件，獲得子數據名列表，根據列表打開指定的數據集。然后獲取子數據集的屬性信息。打開需要的波段，提取波段數據。將數據信息轉換后關閉數據集并顯示。

1）運行環境配置

MODIS衛星數據一個文件包含多個SUBDATASETSGDAL，每個SUBDATASETSGDAL又包含很多波段的數據。GDAL庫的默認支持文件不包含HDF格式，因此需要向GDAL源文件中添加HDF的驅動重新編譯。在GDAL文件中的nmake.opt中修改GDAL_HOME和#Uncomment the following and update to enable NCSA HDF Release 4 support，添加HDF的函數庫路徑。然后用命令工具進行編譯，將編譯好的GDAL中的bin文件夾添加到系統環境變量中。

2）圖像讀寫和處理模塊

在 GDAL讀取柵數據格式文件時首先用GDALAllRegister（）注冊HDF文件驅動。打開指定的HDF文件的子數據集并獲取該子數據集里的波段數據指針。GDALRasterBand::RasterIO（）函數可以用來提取多波段信息，可以根據情況自動進行數據轉換，以及對數據窗口放大和縮小等。需要指出的是，波段數據是以行為單位進行讀取的。通過以上步驟可以對HDF文件指定的波段進行處理，完成處理工作后，可以通過GDALClose（）函數關閉打開的子數據集和HDF文件。GDAL讀取HDF關鍵代碼實現如下：

MODIS數據雖然已經過輻射糾正，但是MODIS L1B數據仍存在自身缺陷，即具有蝴蝶結效應。所以在使用MODIS數據之前，必需要先進行幾何校正。在幾何糾正功能設計上采用投影—插值的方法進行相關編程，這種方法實際上是在投影變換的同時進行了去除蝴蝶結效應。軟件的讀寫模塊主要實現打開、保存和退出，相應的功能按鈕為“打開影像”，“保存”和“退出”。圖像處理模塊主要實現對圖像進行幾何糾正（包括蝴蝶結效應糾正和一般的幾何糾正），相應的功能按鈕為“幾何糾正”。

3）圖像操作模塊

操作模塊包括在對當前地圖窗口的圖形進行縮放，將圖像置為全圖模式，漫游實現手動調整窗口位置。前后視圖實現打開上下視圖。在鷹眼圖中能夠實現對圖像的選擇和操作。 相應的功能按鈕為“放大”，“縮小”，“漫游”，“全圖”，“前視圖”和“后視圖”。

4）6S 模塊

6S大氣輻射模式能準確模擬太陽到目標物，目標物到傳感器路徑上的大氣影響，是發展比較成熟的大氣修正模式。由于6S模式本身相對精確和復雜，并且由FORTRAN語言編寫。所以在添加6S模式時并沒有進行C++語言的重新編寫，而是將FORTRAN語言的源程序編譯運行后生成的debug文件夾下的main.exe用ShellExecute函數調用，將輸入界面的值保存到輸入文件inputfile.txt中。在運行6S模型后，調用debug文件夾下的main.exe后，生成輸出文件output.txt。6S功能模塊相應的功能按鈕為“6S”。

3 功能實現

3.1 讀取和處理遙感圖像

打開影像文件，利用軟件對數據中1、3、7波段進行幾何校正，目的是去除上文所述的蝴蝶結效應和進行一般的幾何糾正。第7波段幾何校正前后圖如圖3，4所示。

圖3 第7波段幾何校正前Fig.3 Band7 before the geometric correction

圖4 第7波段幾何校正后Fig.4 7 band after the geometric correction

經過幾何校正后，可在第7波段圖中尋找可以用于暗像元法反演的點，該點要滿足表觀反射率小于0.05，并且該點附近的暗像元點的邊關反射率變化不大。本論文所采用暗像元點為北緯19.056 461，東經126.009 181的點。在該點處獲得相關參數。

3.2 調用6S大氣輻射傳輸模型

在6S選擇項中跟著提示操作在幾何參數中依次輸入太陽天頂角、衛星天頂角、太陽方位角、衛星方位角4個變量。大氣模式包括7種，依據暗像元所在地區選擇中緯度夏季大氣模式。氣溶膠模式包括8種，將氣溶膠類型設置為海洋性氣溶膠。目標物和傳感器高度選擇100 m。地表狀況處理設為地物面為均勻朗伯表面，無方向效應[7]。

3.3 建立查算表

在6S模型利用上述的幾何參數，結合了不同大氣模式、氣溶膠模式，根據550 nm氣溶膠光學厚度值和暗像元1、3波段的地表反射率反演出該暗像元點的地表反射率。1、3波段查算表建立完成后，可根據該暗像元點的實際地表反射率、表觀反射率分別對照1、3波段查算表找出最匹配的AOD值，兩波段對應AOD的算數平均值可近似等于該暗像元點550 nm氣溶膠光學厚度。第1、3波段表觀反射率查算表的部分如表格1、2所示。

表1 第1波段部分查算表Tab.1 Part of the band1’s lookup-table

4 結 論

文中設計和開發了一種基于暗像元法的大氣氣溶膠光學厚度反演的一體化軟件，軟件設計主要采用GDAL庫等對MODIS遙感數據進行提取和相關處理。通過實際應用表明軟件運行穩定可靠，GDAL和HDF能夠實現完美匹配和兼容，最終能夠實現氣溶膠光學厚度的反演目的。這種軟件的開發能夠有效的簡化了使用者的操作過程，降低了對使用者的要求，提高了效率。

[1]毛節泰，李成才，張軍華，等.MODIS衛星遙感北京地區氣溶膠光學厚度及與地面光度計遙感的對比 [J].應用氣象學報，2002，13（特刊）:127-135.MAO Jie-tai，LICheng-cai，ZHANG Jun-hua，etal.The comparison of remote sensing aerosol optical depth from MODIS data and ground sun-photometer observations[J].Journal of Applied Meterological Science，2002，13 （specialissue）:127-135.

表2 第3波段部分查算表Tab.2 Part of the band3’s lookup-table

[2]張巖岫，王志清，劉立武，等.大氣散射對激光制導武器對抗態勢構建影響研究[J].現代電子技術，2012（21）：38-40，44.ZHANG Yan-xiu，WANG Zhi-qing，LIU Li-wu，et al.Influence of atmospheric scattering on confrontation building in laser guided weapon[J].Modern Electronics Technique，2012（21）：38-40，44.

[3]李曉靜，劉玉潔，邱紅，等.利用MODIS資料反演北京及其周邊地區氣溶膠光學厚度方法研究[J].氣象學報，2003，61（5）:580-591.LI Xiao-jing，LIU Yu-jie，QIU Hong，et al.Retrieval method for optical thickness of aerosds over Beijing and it’s vicinity by using the MODIS data[J].Acta Meteorological Sinica，2003，61（5）:580-591.

[4]Kaufman Y J，Holben B N.Will aerosol measurements from Terra and Aqua polar orbiting satellites represent the daily aerosol abundance and properties[J].Geophysics Research Letters，2000，27（23）:3861-3864.

[5]郭廣猛.關于MODIS衛星數據的幾何糾正方法 [J].遙感信息，2002（3）:26-28.GUO Guang-meng.Geometric calibration of MODIS data[J].Remote Sensing Information，2002（3）:26-28.

[6]張莉，曾致遠.基于HDF4文件格式的MODIS 1B影像數據提取的研究與實現[J].國土資源感，2004（4）:27-32.ZHANG Li，ZENG Zhi-yuan.The study and implementation of extraction MODIS level 1B image data based on a HDF4 file[J].Remote Sensing for Land&Resources，2004（4）:27-32.

[7]Verraote E F，Tanre D，Denze J L，et al.Second simulation of the satellite signal in the solar spectrum,6S:anoverview[J].IEEE Trans Geosci Remote Sens，1997（35）:675-686.