郭飛

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

關鍵詞： 大氣校正； MODTRAN； APDA； 暗像元法

中圖分類號： TN911?34 文獻標識碼： A 文章編號： 1004?373X（2015）10?0023?05

0 引 言

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

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

1 概 述

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

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

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

圖1 大氣校正流程圖

2 MODTRAN快速查找表并行處理

2.1 調用MODTRAN計算大氣參數(shù)

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

利用這些參數(shù)可以計算如下的大氣校正所需的參數(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 并行處理

并行處理是一個對數(shù)據(jù)進行分片很有效的方法，在調用MODTRAN時可以將tape5文件進行分片，獨立生成tape7文件。

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

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

3 水汽含量反演

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

3.1 大氣預處理微分吸收法

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

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

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

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

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

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

3.2 水汽初始含量選擇

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

3.3 水汽含量反演流程圖

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

圖3 水汽含量反演流程圖

4 能見度反演

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

4.1 暗像元的選取

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

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

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

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

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

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

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

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

[ρ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ù)一定的情況下，反射率與輻亮度有如下關系：

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

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

5 大氣輻射校正

當考慮鄰近效應且地面為朗伯面時，反射率與輻亮度有如下關系：

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

式中：[Lp]表示對應波段的程輻射亮度；[S]為大氣球面反照率；[Q]為目標反射率；[Qe]為背景反射率。

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

大氣輻射校正過程如下：

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

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

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

6 實驗分析

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

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

6.1 偽彩色合成對比

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

圖4 校正得到的偽彩色圖像對比結果

6.2 對應坐標光譜特征比較

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

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

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

圖5 對應圖像坐標處ENVI校正結果與本系統(tǒng)校正結果比較

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

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

7 結 語

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

參考文獻

[1] JENSEN J R.遙感數(shù)字影像處理導論[M].陳曉玲，譯.北京：機械工業(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] 亓雪勇，田慶久.光學遙感大氣校正研究進展[J].國土資源遙感，2005（4）：1?6.

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

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

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

[7] 彭妮娜，易維寧，方勇華.400～1 000 nm波段反演氣溶膠光學厚度的暗像元法[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.