大氣氣溶膠光學(xué)厚度反演軟件系統(tǒng)設(shè)計和實現(xiàn)

劉翔龍，楊再昕

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

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

為了實現(xiàn)對反演軟件的開發(fā)，根據(jù)目前已有反演的理論基礎(chǔ)，主要提出和設(shè)計利用GDAL開源庫對MODIS L1B遙感數(shù)據(jù)提取和進(jìn)行相關(guān)處理，選取暗像元后利用其參數(shù)信息結(jié)合6S模型建立查找表。使用QT設(shè)計用戶操作界面，實現(xiàn)反演可視化。本軟件能夠?qū)崿F(xiàn)氣溶膠光學(xué)厚度的反演的功能。

1 反演原理

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

式中L是衛(wèi)星測量輻亮度，Es是大氣頂?shù)奶栞椛渫?，?是整層大氣光學(xué)厚度，（μs，Φs）表示太陽入射光的方向，（μv，Φv）表示衛(wèi)星觀測方向，μs，Φs，μv，Φv分別為太陽入射光和衛(wèi)星觀測天頂角的余弦和方位角。如果不考慮氣體吸收，衛(wèi)星觀測的表觀反射率為：

式中θs為太陽天頂角，θv為觀測天頂角，Φ為相對方位角，Φs分別為太陽到地面和地面到衛(wèi)星的整層大氣 T（θv）透過率，S 為大氣的球面反照率，ρ為地表反射率，F(xiàn)d（θs）為表觀反射率為零時歸一化總的向下輻射通量。當(dāng)?shù)乇矸瓷渎屎苄r，衛(wèi)星觀測反射率主要取決于大氣貢獻(xiàn)項。當(dāng)?shù)乇矸瓷渎屎艽髸r，地表的貢獻(xiàn)成為主要貢獻(xiàn)項。這樣可以得出反演的基礎(chǔ)：如果已知地表反射率，確定大氣氣溶膠模型就可以得到氣溶膠光學(xué)厚度。反之，如果已知氣溶膠光學(xué)厚度和相應(yīng)的大氣參數(shù)，也可以得到地表反射率。

該軟件的設(shè)計主要依據(jù)暗像元方法DDV（Dense Dark Vegetation），它是由Kaufman和Sendra反演稠密植被上空氣溶膠光學(xué)厚度建立的。這種算法主要依據(jù)濃密植被和湖泊的反射率很低（約0.01～0.02），因此可以將有著大面積森林或水域的地方可以作為暗像元，暗像元算法是基于表觀反射率的大氣貢獻(xiàn)項，即利用衛(wèi)星觀測的路徑輻射反演氣溶膠光學(xué)厚度，是氣溶膠遙感應(yīng)用比較常用的算法。利用MODIS圖像反演氣溶膠光學(xué)厚度主要是利用了空間分辨率為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）、藍(lán)（470 nm）和中紅外通道（2 100 nm）的地表反射率存在如下關(guān)系：

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

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

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

2 軟件設(shè)計

2.1 軟件設(shè)計的問題

該氣溶膠光學(xué)厚度反演軟件是基于windows 7系統(tǒng)操作環(huán)境下，編程采用了C++語言，開發(fā)工具選用QT C++，調(diào)用GDAL開源庫、HDF驅(qū)動進(jìn)行讀取、處理HDF格式文件。系統(tǒng)結(jié)構(gòu)功能圖如圖2所示。

圖2 系統(tǒng)結(jié)構(gòu)功能圖Fig.2 Structure diagram of the system’function

對MODIS遙感圖像數(shù)據(jù)需要解決如下問題：

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

2）對MODIS數(shù)據(jù)進(jìn)行幾何糾正，其中包括一般的遙感圖像校正和MODIS數(shù)據(jù)特有的“Bowtie”（蝴蝶結(jié)效應(yīng)糾正）[5]。同時能夠?qū)崿F(xiàn)對遙感圖像的一般操作，如圖像的放大，縮小以及漫游等功能。

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

2.2 軟件設(shè)計的方法

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

1）運行環(huán)境配置

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

2）圖像讀寫和處理模塊

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

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

3）圖像操作模塊

操作模塊包括在對當(dāng)前地圖窗口的圖形進(jìn)行縮放，將圖像置為全圖模式，漫游實現(xiàn)手動調(diào)整窗口位置。前后視圖實現(xiàn)打開上下視圖。在鷹眼圖中能夠?qū)崿F(xiàn)對圖像的選擇和操作。 相應(yīng)的功能按鈕為“放大”，“縮小”，“漫游”，“全圖”，“前視圖”和“后視圖”。

4）6S 模塊

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

3 功能實現(xiàn)

3.1 讀取和處理遙感圖像

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

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

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

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

3.2 調(diào)用6S大氣輻射傳輸模型

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

3.3 建立查算表

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

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

4 結(jié) 論

文中設(shè)計和開發(fā)了一種基于暗像元法的大氣氣溶膠光學(xué)厚度反演的一體化軟件，軟件設(shè)計主要采用GDAL庫等對MODIS遙感數(shù)據(jù)進(jìn)行提取和相關(guān)處理。通過實際應(yīng)用表明軟件運行穩(wěn)定可靠，GDAL和HDF能夠?qū)崿F(xiàn)完美匹配和兼容，最終能夠?qū)崿F(xiàn)氣溶膠光學(xué)厚度的反演目的。這種軟件的開發(fā)能夠有效的簡化了使用者的操作過程，降低了對使用者的要求，提高了效率。

[1]毛節(jié)泰，李成才，張軍華，等.MODIS衛(wèi)星遙感北京地區(qū)氣溶膠光學(xué)厚度及與地面光度計遙感的對比 [J].應(yīng)用氣象學(xué)報，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]張巖岫，王志清，劉立武，等.大氣散射對激光制導(dǎo)武器對抗態(tài)勢構(gòu)建影響研究[J].現(xiàn)代電子技術(shù)，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資料反演北京及其周邊地區(qū)氣溶膠光學(xué)厚度方法研究[J].氣象學(xué)報，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]郭廣猛.關(guān)于MODIS衛(wèi)星數(shù)據(jù)的幾何糾正方法 [J].遙感信息，2002（3）:26-28.GUO Guang-meng.Geometric calibration of MODIS data[J].Remote Sensing Information，2002（3）:26-28.

[6]張莉，曾致遠(yuǎn).基于HDF4文件格式的MODIS 1B影像數(shù)據(jù)提取的研究與實現(xiàn)[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.