基于MODTRAN的高光譜快速大氣校正方法

郭飛

摘 要： 高光譜遙感成像過程中一般會(huì)受大氣成分的影響，為了利用高光譜數(shù)據(jù)進(jìn)行有效的地物識(shí)別，這些影響必須消除。詳細(xì)介紹大氣校正的一個(gè)過程：首先采用大氣預(yù)處理微分吸收法（APDA）和暗像元法，通過并行調(diào)用MODTRAN先后建立水汽和能見度的查找表，從而計(jì)算出圖像中水汽含量和能見度，并利用得出的這些參數(shù)建立反射率查找表，最終對(duì)高光譜圖像進(jìn)行逐像元校正。該過程結(jié)合輻射傳輸模型法和基于圖像信息的校正方法的優(yōu)點(diǎn)，不需要實(shí)時(shí)的大氣參數(shù)也能進(jìn)行絕對(duì)大氣校正。同時(shí)，加入并行處理使運(yùn)行速度顯著提高。通過對(duì)校正后的數(shù)據(jù)進(jìn)行對(duì)比分析，結(jié)果表明校正后的光譜曲線與地面實(shí)物和其他軟件校正后的光譜特征相吻合，該方法能夠滿足大氣輻射校正的基本要求。

關(guān)鍵詞： 大氣校正； MODTRAN； APDA； 暗像元法

中圖分類號(hào)： TN911?34 文獻(xiàn)標(biāo)識(shí)碼： A 文章編號(hào)： 1004?373X（2015）10?0023?05

0 引 言

20世紀(jì)80年代初期成像光譜（高光譜）技術(shù)的出現(xiàn)，使光學(xué)遙感進(jìn)入了嶄新的階段。高光譜遙感能在電磁波譜的可見光和紅外區(qū)域， 獲取許多非常窄的光譜圖像數(shù)據(jù)，為每個(gè)像元提供數(shù)十至數(shù)百個(gè)窄波段（一般小于10 nm）的光譜信息，產(chǎn)生一條完整而連續(xù)的光譜曲線，以便于識(shí)別出地物特征，具有波段多、光譜分辨率高、圖譜合一等特點(diǎn)。大氣校正是遙感信息定量化過程中不可缺少的環(huán)節(jié)。在高光譜圖像的拍攝過程中，由于大氣的存在，進(jìn)入大氣的太陽輻射會(huì)發(fā)生反射、折射、吸收、散射[1]，其中對(duì)傳感器接受影響較大的是吸收和散射。由于大氣的吸收和散射，使得進(jìn)入到傳感器中的信號(hào)減弱[2]。同時(shí)，大氣的散射光也有一部分直接或經(jīng)過背景地物反射進(jìn)入到傳感器時(shí)卻又使得原型號(hào)增強(qiáng)了。為了消除這些因素的影響，必須對(duì)高光譜數(shù)據(jù)進(jìn)行大氣輻射校正。

目前大氣校正的方法主要有基于圖像特征的經(jīng)驗(yàn)線性法、相對(duì)校正方法和基于物理的輻射傳輸模型法[3]。經(jīng)驗(yàn)線性法需要做大量的地物反射率測(cè)量工作，通過統(tǒng)計(jì)的方法建立傳感器接收到的輻亮度與地物反射率的線性回歸方程。基于圖像特征的相對(duì)校正方法得到的不是地物的真實(shí)反射率，而是相對(duì)反射率。它不需要實(shí)際地面光譜和大氣參數(shù)的測(cè)量，可以應(yīng)用于歷史數(shù)據(jù)。它校正后的結(jié)果不是具體的反射率或輻亮度，相同的結(jié)果只表示相同的反射率，至于多少不能給定[4]。輻射傳輸模型法建立在大氣輻射傳輸方程的基礎(chǔ)上，在進(jìn)行大氣輻射校正時(shí)需要輸入各種大氣參數(shù)（例如：衛(wèi)星過境時(shí)的地物反射率、大氣的能見度、太陽天頂角和衛(wèi)星傳感器的標(biāo)定參數(shù)等），獲得這些參數(shù)的代價(jià)昂貴，又不易實(shí)現(xiàn)，尤其對(duì)于歷史存檔數(shù)據(jù)，獲取這些參數(shù)更是不可能的[5]。陶東興，趙慧潔通過調(diào)用MODTRAN分別建立水汽吸收通道和2個(gè)參考通道的2個(gè)查找表[6]，搜索這2個(gè)查找表使用三次樣條插值方法計(jì)算水汽含量，實(shí)驗(yàn)結(jié)果表明生成的水汽含量分布圖精度與MODTRAN計(jì)算結(jié)果只有0.1%的誤差。另外還有暗像元法，暗目標(biāo)方法通過確定暗目標(biāo)在可見光波段的地表反射率來估計(jì)氣溶膠參數(shù)。趙祥等選取小麥作為新的暗目標(biāo)[4]，使用波段間的線性回歸與插值的方法對(duì)常規(guī)暗目標(biāo)方法進(jìn)行改進(jìn)，從而將暗目標(biāo)方法擴(kuò)展應(yīng)用于高光譜遙感數(shù)據(jù)的大氣校正。彭妮娜等提出了濃密植物像元660 nm與830 nm波段處的線性關(guān)系進(jìn)行氣溶膠含量反演[7]，為獲取氣溶膠含量提供了一種新途徑。

1 概 述

為了保證大氣校正的精確性和程序運(yùn)行速度，結(jié)合前人的工作經(jīng)驗(yàn)，總結(jié)出一套大氣輻射校正方法，形成一整套大氣輻射校正的具體流程，在流程中加入并行計(jì)算提高運(yùn)行效率。

該方法具體思路是通過大氣預(yù)處理微分吸收法（APDA法）并調(diào)用MODTRAN建立的水汽查找表來獲得水汽含量，然后采用暗像元法并調(diào)用MODTRAN建立的能見度查找表獲得能見度，最后利用水汽含量和能見度參數(shù)調(diào)用MODTRAN建立反射率查找表，計(jì)算FLAASH模型中的大氣輻射校正的參數(shù)從而實(shí)現(xiàn)大氣輻射校正的目標(biāo)。

運(yùn)用該方法開發(fā)了用于大氣輻射校正的軟件。并對(duì)調(diào)用MODTRAN建立查找表進(jìn)行并行處理，提高運(yùn)行效率。該軟件包括水汽反演、能見度反演、大氣輻射校正這三大模塊。為了驗(yàn)證這三大模塊的精確性，對(duì)東天山地區(qū)高光譜數(shù)據(jù)（已完成輻射定標(biāo)和幾何校正）進(jìn)行大氣輻射校正，將校正后的光譜和地面采集的光譜進(jìn)行對(duì)應(yīng)像元對(duì)比，分析校正結(jié)果的準(zhǔn)確性。本文方法的大體思路流程見圖1。

圖1 大氣校正流程圖

2 MODTRAN快速查找表并行處理

2.1 調(diào)用MODTRAN計(jì)算大氣參數(shù)

在調(diào)用MODTRAN后會(huì)生成tape7文件，該文件得到一系列參數(shù)（FREQ，TOTAL RAD，GRND RFLT，DRCT RFLT等）。

利用這些參數(shù)可以計(jì)算如下的大氣校正所需的參數(shù)：

[λ=10×107FREQ] （1）

[Lp=（TOTAL RAD）0] （2）

[S=（（GRND RFLT）1.0-（GRND RFLT）0）（（GRND RFLT）1.0-（GRND RFLT）0）-→ ←-2?（（GRND RFLT）1.0-（GRND RFLT）0）（（GRND RFLT）1.0-（GRND RFLT）0）] （3）

[Gb=（（TOTAL RAD）1.0-（TOTAL RAD）0）·（1-S）] （4）

[Gt=（（DRCT RFLT）1.0-（DRCT RFLT）0）·（1-S）] （5）

式中：[λ]表示波長；[Lp]表示程輻射；[S]為大氣球面反照率；[Gb]，[Gt]為輻射校正參數(shù)。

2.2 并行處理

并行處理是一個(gè)對(duì)數(shù)據(jù)進(jìn)行分片很有效的方法，在調(diào)用MODTRAN時(shí)可以將tape5文件進(jìn)行分片，獨(dú)立生成tape7文件。

以建立水汽查找表為例，為了提高運(yùn)行速度，在寫tape5文件時(shí)，將12個(gè)不同的水汽含量和3種不同的反射率共36種情形分別寫入到12個(gè)文件夾中，在這12個(gè)文件夾中獨(dú)立調(diào)用以達(dá)到減少程序運(yùn)行時(shí)間的目的，在建立查找表需要整個(gè)數(shù)據(jù)時(shí)再將這些數(shù)據(jù)合并。而且高光譜數(shù)據(jù)處理的典型特點(diǎn)是像素間、波段間依賴關(guān)系比較小，相互獨(dú)立性強(qiáng)，在不同波段數(shù)據(jù)間也可以進(jìn)行分片。

因此在3個(gè)模塊中，涉及了逐個(gè)波段計(jì)算時(shí)，一般都對(duì)每個(gè)波段進(jìn)行并行計(jì)算處理。以建立水汽查找表為例，并行操作的流程圖如圖2所示。

3 水汽含量反演

水汽含量的影響包括對(duì)大氣透過率的影響和程輻射的影響，但是它卻不是惟一的因素，氣溶膠和大氣分子的散射也對(duì)大氣透過率有不可忽視的影響。本系統(tǒng)為了區(qū)分這兩種因素的影響，根據(jù)水汽吸收具有波段選擇性，而氣溶膠和大氣分子的散射是關(guān)于波長的緩變函數(shù)這一特點(diǎn)，選用大氣預(yù)處理微分吸收法（Atmospheric Precorrected Differential Absorption，APDA）[7]。

3.1 大氣預(yù)處理微分吸收法

該方法利用1 135 nm處的水汽吸收波段及其鄰近的非水汽吸收波段的比值來獲取大氣水汽柱含量。在特定的大氣條件下，水汽含量值與參數(shù)[RAPDA]有如下經(jīng)驗(yàn)關(guān)系：

[RAPDA=e-（γ+α（PW）β）] （6）

式中：PW是水汽含量；[α，β，γ]是未知的參數(shù)，在特定不變的大氣條件下[α，β，γ]的值不變；[RAPDA]的計(jì)算方法如下：

[RAPDA=[Lm-Lp，m]LIR（[λr]j，[Lr-Lp，r]j）λm] （7）

式中：[Lm]是水汽吸收通道的輻射值；[Lp，m]是水汽吸收通道的程輻射，在MODTRAN中通過設(shè)置地面反射率為0來獲取各波段的程輻射；[Lr]是水汽參考通道的輻射值；[Lp，r]是水汽參考通道的程輻射；[LIR（[λr]j，[Lr-Latm，r]j）λm]是指通過線性擬合參考通道的輻射值，估計(jì)在水汽吸收通道應(yīng)有的輻射值。

圖2 建立水汽查找表流程圖

3.2 水汽初始含量選擇

水汽含量的初始值的選擇直接關(guān)系著最終結(jié)果的正確性。關(guān)于如何設(shè)置水汽含量初值，ENVI的幫助文檔中給出了設(shè)置方法：如果沒可用的水汽信息，則根據(jù)已知的或預(yù)期的大氣溫度來設(shè)置水汽含量。如果溫度是未知的，則根據(jù)數(shù)據(jù)獲取時(shí)間的季節(jié)性來設(shè)置?？偟膩碚f，它是基于一個(gè)季節(jié)性緯度表面溫度模型的設(shè)置方法[8]。

3.3 水汽含量反演流程圖

利用大氣預(yù)處理微分吸收法求解水汽含量具體流程圖如圖3所示。

圖3 水汽含量反演流程圖

4 能見度反演

大氣中的氣溶膠粒子將會(huì)吸收和散射太陽輻射，從而使衛(wèi)星傳感器接收到的輻射信號(hào)產(chǎn)生衰減，是影響能見度的主要因素。本項(xiàng)目基于精確度需要以及避免地面輻射定標(biāo)的巨大工作量，選用暗像元法反演能見度。所謂暗像元是指某些波段反射率接近為0的像元。暗像元可以選擇水體、密度較高的植被森林或陰影地區(qū)。

4.1 暗像元的選取

該方法的關(guān)鍵點(diǎn)在于暗像元的選取，暗像元選擇的科學(xué)、正確與否關(guān)系到算法的精度。Kaufman等首先采用歸一化植被指數(shù)[9]（Normalized Difference Vegetation Index，NDVI）進(jìn)行暗目標(biāo)確認(rèn)：

[NDVI=ρ860-ρ660ρ860+ρ660>0.7ρ2100<0.05] （8）

式中：[ρ860，ρ660，ρ2100]分別表示在860 nm波段、660 nm波段、2 100 nm波段的地物反射率，若滿足式（8）則是暗像元，否則不是暗像元。

但是由于NDVI對(duì)氣溶膠的變化較為敏感[10]，采用對(duì)大氣不敏感的大氣耐抗植被指數(shù)[11]（Atmospherically Resistant Vegetation Index，ARVI），在已得到的暗像元位置上繼續(xù)進(jìn)行再次選擇得出最終的暗像元位置。該指數(shù)是通過地面調(diào)查得到地物類別在影像中的精確位置，經(jīng)過分析后確定暗目標(biāo)選取方式的經(jīng)驗(yàn)公式如下：

[ARVI=ρ*884+ρ*467-2ρ*660ρ*884-ρ*467+2ρ*660>0.5ρ*2042<0.04ρ*854>0.4] （9）

式中：[ρ*i]表示中心波長為i（單位：nm）的波段的表觀反射率。

4.2 基于輻射傳輸方程的大氣校正參數(shù)計(jì)算

Kaufman等人經(jīng)過大量研究發(fā)現(xiàn)，由于2.1 μm波長比大部分氣溶膠微粒的直徑要大，故該波段受氣溶膠影響可以忽略，當(dāng)?shù)乇矸瓷渎瘦^低時(shí)，即為暗像元時(shí)，2.1 μm表觀反射率與0.47 μm及0.66 μm的地表反射率存在很好的線性關(guān)系：

[ρ470=ρ*2 1004≈ρ2 1004ρ660=ρ*2 1002≈ρ2 1002] （10）

式中：[ρ*2 100]表示中心波長為2 100 nm的波段的表觀反射率；[ρ470，ρ660]表示中心波長分別為470 nm，660 nm處的地表反射率。因此根據(jù)式（10）可以得出暗像元在470 nm和660 nm處的地物反射率。在MODTRAN模型中，在校正參數(shù)一定的情況下，反射率與輻亮度有如下關(guān)系：

[L=Lp+（Gb+Gt）ρ1-Sρ] （11）

式中：[L]表示某一波段傳感器接收到的輻亮度；[Lp]表示對(duì)應(yīng)波段的程輻射亮度；[S]為大氣球面反照率；[ρ]為地表反射率，[Gb]，[Gt]為輻射校正參數(shù)，[Gb+Gt]可以通過能見度查找表查找出[Lp]，[S]，根據(jù)式（11）以及在470 nm和660 nm處的地物反射率和輻亮度來計(jì)算。

5 大氣輻射校正

當(dāng)考慮鄰近效應(yīng)且地面為朗伯面時(shí)，反射率與輻亮度有如下關(guān)系：

[L=Lp+Q1-SQeGb+Qe1-SQeGt] （12）

式中：[Lp]表示對(duì)應(yīng)波段的程輻射亮度；[S]為大氣球面反照率；[Q]為目標(biāo)反射率；[Qe]為背景反射率。

利用水汽含量反演模塊和能見度反演模塊中的方法的到的水汽含量和能見度建立反射率查找表，根據(jù)式（12）的關(guān)系完成反射率與輻亮度的轉(zhuǎn)換。

大氣輻射校正過程如下：

（1） 將已得到的平均水汽含量、能見度作為已知建立反射率查找表；

（2） 通過查找表求解出式（8）中的參數(shù)：[Lp]，[S]，[Gt]，[Gb]等。

（3） 用計(jì)算到的參數(shù)和式（8）校正高光譜圖像得到大氣輻射校正后結(jié)果圖像。

6 實(shí)驗(yàn)分析

為了驗(yàn)證所設(shè)計(jì)的大氣校正模塊的準(zhǔn)確性，將在東天山土屋—延?xùn)|地區(qū)拍攝的高光譜數(shù)據(jù)作為實(shí)驗(yàn)對(duì)象，做大氣輻射校正。把處理后結(jié)果與地面拍攝的實(shí)際的光譜曲線和其他軟件的校正結(jié)果相比較，計(jì)算同一坐標(biāo)下各校正結(jié)果的光譜特征，驗(yàn)證本模塊的有效性。

本次遙感數(shù)據(jù)獲取使用了澳大利亞HyVista公司研制的HyMap成像光譜儀，于2002年10月10日格林威治時(shí)間8時(shí)28分30秒獲取東天山土屋?延?xùn)|地區(qū)（北緯42.156 174°，東經(jīng)94.245 092°，平均地面海拔910.433 960 m，飛行高度為3 377.978 153 m）的航空成像光譜信息，波長范圍為可見光、近紅外（VNIR）、短波紅外（SWIR）和中熱紅外（MIR和TIR），波段共124個(gè)，光譜分辨率約15 nm，圖像大小651×11 866，像元大小4.9 m×4.9 m，原始數(shù)據(jù)大小約1.3 GB。該遙感數(shù)據(jù)還提供了地物真實(shí)反射率數(shù)據(jù)，為評(píng)價(jià)大氣輻射校正結(jié)果提供了有效的參考。

6.1 偽彩色合成對(duì)比

利用本方法校正后的圖像與反射率圖像以及ENVI中的FLAASH模塊校正后的圖像對(duì)比如下：

圖4 校正得到的偽彩色圖像對(duì)比結(jié)果

6.2 對(duì)應(yīng)坐標(biāo)光譜特征比較

校正結(jié)果與ENVI中FLAASH模塊校正結(jié)果比較如圖5所示。校正結(jié)果與地面實(shí)際光譜比較如圖6所示。

在校正前后的結(jié)果圖像中，同一坐標(biāo)下兩個(gè)光譜曲線的光譜角和光譜相關(guān)系數(shù)如表1所示。在表1中“A”表示ENVI的處理結(jié)果與地面測(cè)量出的輻亮度光譜曲線做比較，“B”表示本文所用方法的處理結(jié)果與地面測(cè)量出的輻亮度光譜曲線做比較。從表1可以看出，本項(xiàng)目校正后的光譜與地面實(shí)際光譜相比，光譜角在0.1～0.24之間，相關(guān)系數(shù)在0.90～0.99之間。說明校正結(jié)果穩(wěn)定，誤差在可接受范圍內(nèi)，與地面實(shí)際光譜較相識(shí)，達(dá)到了預(yù)期目標(biāo)。

表1 校正結(jié)果與地面反射率光譜比較

圖5 對(duì)應(yīng)圖像坐標(biāo)處ENVI校正結(jié)果與本系統(tǒng)校正結(jié)果比較

同時(shí)，我們注意到校正過程中存在難免的誤差，為此做如下可能性分析：系統(tǒng)如果在調(diào)用MODTRAN時(shí)建立的查找表數(shù)目偏少，會(huì)造成一定的誤差，如水汽反演過程中得到了經(jīng)過第一次迭代計(jì)算的水汽含量，需要在水汽查找表中查找對(duì)應(yīng)的程輻射，通過線性插值得到，而如果建立更加細(xì)致的水汽含量查找表則減小誤差，但必須付出更多的時(shí)間代價(jià)；另外在水汽含量迭代計(jì)算過程開始時(shí)有一個(gè)水汽含量初值的選擇，若初值選擇不當(dāng)則會(huì)使誤差擴(kuò)大；在能見度反演模塊中，將水汽含量的平均值作為每一個(gè)像素點(diǎn)下的水汽含量，而在一幅高光譜圖像中，如果水汽含量分布極不均勻，那么必定會(huì)造成一定的誤差，因此在必要時(shí)，需要將圖像拆分成更小的數(shù)據(jù)塊進(jìn)行單獨(dú)的大氣校正。

圖6 對(duì)應(yīng)圖像坐標(biāo)處本系統(tǒng)校正后光譜和地面拍攝反射率光譜

7 結(jié) 語

從以上實(shí)驗(yàn)可以得出本系統(tǒng)的結(jié)果與ENVI校正結(jié)果較為吻合，校正后輻亮度值與ENVI的結(jié)果或者地面拍攝數(shù)據(jù)相比都比較吻合，只有在受水汽影響的波段誤差偏大，在實(shí)際應(yīng)用時(shí)只需要去除受水汽影響嚴(yán)重的波段即可。而且利用并行加速，可以大大減少輻射校正的運(yùn)行時(shí)間，提高運(yùn)行效率，使大氣輻射校正達(dá)到了既快速又精確的目標(biāo)。最后通過光譜特征分析，驗(yàn)證了大氣輻射校正的有效性。

參考文獻(xiàn)

[1] JENSEN J R.遙感數(shù)字影像處理導(dǎo)論[M].陳曉玲，譯.北京：機(jī)械工業(yè)出版社，2007.

[2] FRASER R S， FERRARE R A， KAUFMAN Y J， et al. Algorithm for atmospheric corrections of aircraft and satellite imagery [J]. Int J Remote Sens， 1992， 13： 541?557.

[3] 亓雪勇，田慶久.光學(xué)遙感大氣校正研究進(jìn)展[J].國土資源遙感，2005（4）：1?6.

[4] 丁麗霞，周斌，王人潮.遙感監(jiān)測(cè)中5種相對(duì)輻射校正方法研究[J].浙江大學(xué)學(xué)報(bào)：農(nóng)業(yè)與生命科學(xué)版，2005，31（3）：269?276.

[5] 陶東興，趙慧潔，賈國瑞.基于MODTRAN的雙查找表法反演高光譜MODTRAN數(shù)據(jù)的水汽含量[J].光電學(xué)報(bào)，2012，32（10）：1?6.

[6] 趙祥，梁順林，劉素紅，等.高光譜遙感數(shù)據(jù)的改正暗目標(biāo)大氣校正方法研究[J].地球科學(xué)，2007，37（12）：1653?1659.

[7] 彭妮娜，易維寧，方勇華.400～1 000 nm波段反演氣溶膠光學(xué)厚度的暗像元法[J].紅外與激光工程，2008，37（5）：878?883.

[8] SCHLAPFER D， BOREL C C， KELLER J， et al. Atmospheric precorrected differential absorption technique to retrieve colummar water vapor [J]. Remote Sens Environ，1998，65： 353?366.

[9] KAUFMAN Y J， SENDRA C. Algorithm for automatic atmospheric corrections to visible and near?IR satellite imagery [J]. Int J Remote Sens，1988， 9（8）： 1357?1381.

[10] KAUFMAN Y J. Atmospherically resistant vegetation index （ARVI） for EOS?MODIS [J]. IEEE Transactions on Geoscience and Remote Sensing ，1992， 30（2）： 261?270.