,,

(1.解放軍航天工程大學 研究生院,北京 101416；2.解放軍航天工程大學 電子與光學工程系,北京 101416)

0 引言

星載高光譜遙感的射傳輸路徑是：太陽光譜輻射經過大氣反射、散射、吸收等影響后，到達地表；在被地物反射后，其反射輻射再次經過大氣影響，并與大氣直接反射的太陽輻射一起構成傳感器入瞳前的輻亮度數(shù)據(jù)，進入傳感器中的輻射能量通過傳感器的輻射維、空間維、光譜維離散采樣后輸出數(shù)字圖像。這一過程可以分為場景系統(tǒng)和傳感器系統(tǒng)。其中，場景系統(tǒng)作為整個高光譜遙感體系中影響因素最復雜多變的部分，尚未建立完備的模型來全面分析。

目前主要有兩種方法對地面場景高光譜特性仿真。第一種是人工描繪場景中各物體的三維形貌，并利用材質光譜結合紋理圖設定各物體的反射特性[1]，第2種利用高分辨率的高光譜數(shù)據(jù)描述場景反射率分布[2]。對于分析單一因素對高光譜成像的影響能力來說，以上兩種方法顯得過于復雜。本文以傳感器入瞳前的大氣輻亮度傳輸模型為依據(jù)，通過簡化場景的幾何形態(tài)，并將其分為地面系統(tǒng)和大氣系統(tǒng)兩部分，利用地面反射率數(shù)據(jù)和MODTRAN大氣傳輸軟件，實現(xiàn)對高光譜場景的仿真模擬。

1 大氣輻亮度傳輸模型

多年來，研究者根據(jù)輻射在地表與大氣之間的傳輸過程出發(fā)，針對不同的地表特性，主要建立了以下5種典型的大氣輻射傳輸模型[3]：

1)不考慮地-氣多次反彈的模型

2)基于地表均一的朗伯面模型

3)基于地表均一的非朗伯面模型

4)基于地表非均一的朗伯面模型

5)基于地表非均一的非朗伯面模型

以上5種模型的理想化程度逐漸降低，計算精度和要考慮的參數(shù)也越來越復雜。

在高光譜遙感成像過程中，太陽光譜輻射經過大氣反射、散射、吸收等影響后，到達地表；在被地物反射后，其反射輻射再次經過大氣傳輸?shù)竭_傳感器入瞳處。

如圖1所示，將入射到傳感器前的總輻射能量分成三部分：

(1)來自太陽輻射經過大氣傳輸?shù)竭_目標表面并被反射至傳感器前能量：圖中路徑a和路徑b(路徑a是太陽輻射經大氣傳輸后直接照射地物目標，路徑b是太陽輻射經大氣多次散射后的輻射能量照射到地物目標后直接進入傳感器)；(2)來自太陽輻射經過大氣傳輸照射到背景景物并被反射至傳感器前能量：圖中路徑c和路徑d(路徑c是太陽輻射經大氣傳輸直接照射到背景景物，路徑d是太陽輻射經過大氣多次散射后照射到背景景物)；(3)路徑e是大氣路徑程輻射，即太陽輻射被大氣多次散射后，不經過地物場景反射，直接進入傳感器視場內的能量。

圖1 大氣輻射傳輸路徑

文獻[4]總結了包含多次散射項的傳感器入瞳處輻亮度。將地面看做是非均勻的，反射率為ρ(λ)的朗伯體時：

(1)

在單一角度的天頂觀測情況下，把場景看作非均一的朗伯面的模型時，上式可以等效為[5-7]：

(2)

式(2)中，L為入瞳前傳感器處像元接收的總輻亮度；ρ為目標像元表面反射率；ρe為鄰近像元表面平均反射率；S為大氣朝下的球面反照率；La為大氣程輻射；A、B為基于大氣條件和幾何光照條件的系數(shù)。

式中所有的變量都為光譜波長的函數(shù)，為方便起見，波長指數(shù)被省略。等式右邊第一部分代表太陽輻射經大氣傳輸?shù)降乇矸瓷渲苯舆M入傳感器前的輻亮度；第二部分為背景景物像元經大氣散射進入傳感器前的輻亮度；第三部分即為大氣程輻射。

等式中ρe需要通過大氣點擴散函數(shù)對ρ處理得到，而A、B、S和La可通過大氣傳輸模型MODTRAN設定觀測視場角、平均海拔高度、大氣模型、氣溶膠類型、能見度范圍等參數(shù)進行仿真反演[8]。

2 鄰近像元反射率計算

被鄰近像元反射的太陽光子通過與大氣粒子的碰撞作用被遙感器接收，不考慮地氣間多次反彈，即認為鄰近像元對遙感器接收到的總輻亮度貢獻率的空間分布函數(shù)等同于光子射向地面并與大氣的隨機碰撞過程，最終在地面形成的光子數(shù)的相對密度空間分布函數(shù)[9]。

蒙特卡洛方法[10]是一種基于“隨機數(shù)”的計算方法，將所求解的問題同一定的概率模型相聯(lián)系，用計算機實現(xiàn)統(tǒng)計模擬或抽樣，以獲得問題的近似解，故又稱隨機抽樣法或統(tǒng)計試驗法。因此，選擇用蒙特卡洛方法來仿真大氣碰撞的隨機過程。

圖2 光子傳輸路徑

光子在大氣中的傳輸路線如圖2所示。其中z為垂直高度，x，y為水平位置坐標。光子初始位置為A點，目標位置為B點，B點在A點水平面方向的投影為B’。d為自由路程，即光子在連續(xù)兩次碰撞之間經歷的路程。θ是線段BA與z軸反方向的夾角，稱為散射極角；Φ為BA與x軸正方向的夾角，稱為散射方位角。散射方位角Φ是在[0，2π]內的隨機數(shù)，散射極角θ是在[0，π]內的隨機數(shù)。

由以上分析可知，求解光子的傳輸路徑可等效為求解d、θ和Φ的值。

引入光學距離l的定義：沿自由路程對d削弱系數(shù)的積分?？捎墒?3)表示：

(3)

上式中，Ks為大氣散射系數(shù)，Ka為大氣吸收系數(shù)。l為無因次量，令非碰撞概率密度函數(shù)P與l的關系為：

P=e-l

(4)

則光子和大氣粒子距離為l時的非碰撞概率可表示為：

(5)

其中：r1為區(qū)間[0,1]的隨機數(shù)。當l趨近無窮時，r1趨近1，即相距無窮遠的光子必然不發(fā)生碰撞，反之同理。則l可由式(6)求得：

l=-ln(1-r1)=|ln(1-r1)|

(6)

大氣光學厚度定義為大氣削弱系數(shù)在垂直路徑自上而下的積分：

(7)

對于水平均勻、垂直分布的大氣，光學厚度τ只和垂直高度Z有關。對應的大氣光學厚度之差：△τ=lcosθ；根據(jù)式(7)中光學厚度與垂直高度的關系，可求得自由路程d=△z/cosθ，光子的水平位移量分別為：△x=dsinθcosΦ和△y=dsinθsinΦ，其中Φ=2πr2，θ=πr3，r2、r3分別為[0，1]內的隨機數(shù)。

3 大氣傳輸計算

3.1 大氣傳輸影響因素

地球表層的大氣主要分布在距離地球表面約80～90 km的低層大氣中，高層大氣對太陽輻射傳輸幾乎沒有影響。底層大氣中存在著多種氣體以及微粒，如塵埃、煙、霧、風、雨、雪等，使光的傳輸特性發(fā)生變化。影響大氣傳輸?shù)淖钪饕蛩厥谴髿夥肿雍蜌馊苣z的影響[11]。

大氣分子主要包括N2、O2、CO2、H2、H2O等，其中占N2對遠紅外和微波會表現(xiàn)出強吸收；O2在可見光近紅外波段具有0.76 μm、0.688 μm、0.628 μm3個強吸收波段。CO2、H2O是對可見光和近紅外吸收的重要吸收分子。在太陽光譜范圍內，2 μm、1.6 μm、1.4 μm是CO2的3個主要吸收通道；水汽對太陽輻射的吸收主要位于近紅外及短波紅外區(qū)，其吸收強度與大氣壓、溫度和光譜波長緊密相關。除上述分子外，大氣中還含有He、Ne、Ar、Xe、O3等少量分子，雖然在可見光和近紅外區(qū)有較強的吸收譜線，但是含量很少，可不考慮吸收作用。

大氣中除大氣分子外，還存在尺寸較大的大量固態(tài)、液態(tài)微粒，包括塵埃、煙粒、微水滴、鹽粒及有機微生物等，這些微粒在大氣中的懸浮成溶膠狀態(tài)，稱為氣溶膠。氣溶膠的光學特性、粒子形態(tài)等特性會因為時空不同發(fā)生劇烈變化，深刻影響大氣輻射傳輸過程。氣溶膠的尺寸一般為0.001～10 μm。氣溶膠一方面使太陽輻射在傳輸過程中發(fā)生散射，使返回大氣層的輻射得到增強，有色氣溶膠的吸收特性還會使到達地表的太陽輻射產生衰減；另一方面，已經到達大氣層內的太陽輻射又難以被反射出大氣層。大氣中存在的氣溶膠主要來自兩方面：一是火山煙、風暴中的揚塵、植物的花粉、海水蒸發(fā)的鹽粒等自然現(xiàn)象；而是諸如化石燃料燃燒、汽車尾氣等人類活動。研究者把常見的氣溶膠粒子按照不同比例混合來定義其類型，例如可用于干旱地區(qū)的沙漠型、用于受海洋季風影響嚴重的海洋型以及具有不同地物屬性的大陸型[12]。

3.2 MODTRAN軟件

MODTRAN是一個中分辨率大氣輻射傳輸模型，由美國空軍研究實驗室(ARFL)與光譜科技公司利用FORTRAN語言開發(fā)。它把0至100 km的整層大氣當成36個均一水平的薄層氣體疊置，每薄層的光學特性通過溫度、氣壓、氣體分子的不同組合來刻畫。內置了多種大氣模型，能在0到50000 cm-1波數(shù)范圍內以1.0 cm-1的分辨率進行輻射傳輸運算。MODTRAN利用二流近似、離散縱標等多次散射方法，計算包含大氣吸收與多次散射作用的透過率、輻照度等。目前已更新至MODTRAN6版本，但由于知識產權原因，國內一般研究人員仍普遍使用MODTRAN4版本。

在查閱相關MODTRAN文獻后，發(fā)現(xiàn)目前絕大部分以性能介紹及其在遙感圖像處理中大氣校正應用為主題。有關于對入瞳處輻亮度計算原理和詳細實現(xiàn)過程的說明則非常少。本文利用MODTRAN在windows環(huán)境下運行的可視化界面軟件PcModwin進行仿真分析。

通過多次實踐以及查看輸出文件MODOUT(包括Modout1、Modout2、Modout3，其中Modout1包含輸出的完整結果)，初步總結了其關于輻亮度的計算過程：(1)選擇內置的大氣廓線數(shù)據(jù)存儲于公共數(shù)據(jù)中；(2)將輸入數(shù)據(jù)中確定的光譜分辨率和光譜范圍等效為參與輻射傳輸計算的波數(shù)范圍；(3)將輸入數(shù)據(jù)的幾何觀測關系等效為由傳感器位置、觀測天頂角、方位角幾何關系；(4)利用輸入的氣溶膠數(shù)據(jù)，利用氣溶膠吸收和非對稱系數(shù)，構建沿“地面到傳感器高度”垂直路徑的單次散射幾何關系，計算單次散射路徑氣溶膠吸收輻射量；(5)在計算沿視線方向的單次散射幾何關系基礎上，得到多次散射參數(shù)；(6)根據(jù)氣體分層理論，利用多次嵌套循環(huán)方法，計算整層大氣的光譜輻亮度等參數(shù)；(7)重采樣1cm-1波數(shù)的光譜輻亮度數(shù)據(jù)，計算結束。

MODTRAN按照計算目的不同可以分為透過率(Transmittance)、熱輻射亮度(Thermal Radiance)、散射輻射亮度(Radiance with Scattering)、太陽直射照度(Direct Solar Irradiance)4種運行模式。

其中，散射輻亮度模式輸出結果包含表1所示內容。

表1 MODTRAN散射輻亮度模式部分輸出結果

利用MODTRAN的散射輻亮度模式可以獲得在特定地表反射率狀態(tài)下的傳感器入瞳處輻亮度，而后通過相關數(shù)學逆推，得到大氣傳輸參數(shù)A、B、S、La。

4 仿真結果

4.1 地面場景仿真

建立理想地面場景如圖3所示，將該場景當作非均勻的朗伯體。設定傳感器空間分辨率20 m，光譜分辨率10 nm，覆蓋0.4～2.4 μm范圍，軌道高度600 km并以推掃方式對地垂直成像。把地表按照傳感器空間分辨率大小分割為若干面元，并假設每一面元的反射率均勻分布。

圖3 地面場景分布情況

地面反射反射率數(shù)據(jù)ρ來自美國USGS光譜庫[13]，如圖4所示，分別為典型植被和混凝土道路的反射率。典型植被的波譜特性是在可見光波段的0.55 μm附近有一個小的反射峰，在0.45 μm和0.65 μm附近有2個吸收谷；在近紅外波段，0.7～0.8 μm處反射率急速升高，并在0.8～1.3 μm處形成一個平緩的高反射率平臺，在1.4 μm、1.9 μm附近各有一個吸收谷。與混凝土的反射率情況相比，在波長0.8～1.3 μm處的反射率大于混凝土道路，在其他波長時情況相反。

圖4 不同地物反射率曲線

圖5、圖6分別描述了地表反射率分別為0.5 μm和1.0 μm波長時的仿真結果，可見該結果與圖4的描述相符，即在0.5 μm時混凝土反射率高于典型植被，在1.0 μm時情況相反。

圖5 波長為0.5 μm時地面反射率

圖6 波長為1.0 μm時地面反射率

4.2 鄰近像元反射率仿真結果

根據(jù)第2節(jié)所述，利用蒙特卡羅方法仿真大氣點擴散函數(shù)的過程如圖7所示。

圖7 蒙特卡洛仿真點擴散函數(shù)流程

經過不斷循環(huán)，直到光子到達地面，得到圖8的光子分布密度圖。

圖8 光子密度圖

蒙特卡洛法逼真描述了光子的傳輸過程，其結果是靠大量的統(tǒng)計數(shù)據(jù)的分析獲得的。得到點擴散函數(shù)后，對4.1節(jié)中的反射率數(shù)據(jù)ρ進行平滑處理可以得到鄰近像元反射率ρe。如圖9所示，ρe(0.5)為波長0.5 μm時的鄰近像元反射率情況。

圖9 波長0.5μm時的鄰近像元反射率

4.3 大氣參數(shù)仿真計算

利用MODTRAN計算大氣傳輸特性的參數(shù)。設置模式如表2所示，在運行MODTRAN時選擇DISORT多次散射算法，通常設置散射算法流數(shù)為8，同時在地面反射率設置中，將地面反射率(Ground Albero)分別設置為0、0.5、1[14]，得到3種情況下總的TOT(傳感器接收的總輻亮度)、PATH(多次散射程輻射)、TGR(總地面反射輻亮度)。

表2 包含多次散射的輻亮度計算步驟

通過MODTRAN仿真得到以上結果后，則關于式(1)的相關參數(shù)S、A、B、La可由式(8)仿真得到：

(8)

式中,PATH0是地面反射率為0的大氣程輻射值；△TGR1是地面反射率分別為1和0時地面總反射輻亮度的差值；△PATH1是地面反射率分別為1和0時的大氣程輻射差值；△TOT1是地面反射率分別為1和0時傳感器入瞳前總輻亮度差值；△TOT0.5地面反射率分別為0.5和0時傳感器入瞳前總輻亮度差值。

由式(8)計算的參數(shù)S、A、B以及程輻射La的結果為圖10，需要注意的是，S中的峰值為算式中分母為0時出現(xiàn)的極化現(xiàn)象，上述求得的幾類參數(shù)均為需求大氣條件時的值，當大氣條件改變時，反演得到的參數(shù)值隨之改變。

圖10 A、B、La、S仿真結果

4.4 地面場景的入瞳處輻亮度

利用4.1節(jié)、4.2節(jié)和4.3節(jié)計算結果，結合大氣傳輸公式，就得到高光譜遙感入瞳處輻亮度數(shù)據(jù)。如圖11、12所示為波長為0.49 μm時，地面場景的入瞳處輻亮度仿真結果。

圖11 波長0.49 μm時入瞳處輻亮度

圖12 波長0.49 μm時輻亮度灰度圖

從圖中可以看出，雖然混凝土道路上像元的入瞳處輻亮度仍然大于兩邊植被對應的像元輻亮度，但與反射率數(shù)據(jù)相比，出現(xiàn)了明顯弱化，這主要是因為在大氣傳輸過程中的衰減造成的。

5 總結

本文采取基于地表非均勻的朗伯面大氣場景模型，利用蒙特卡羅算法仿真大氣點擴散函數(shù)，得到了包含鄰近效應的地表平均反射率；在對MODTRAN軟件原理分析基礎上，考慮多次散射效應，獲得了傳感器入瞳處的輻亮度數(shù)據(jù)。最后對給定的波長為0.4～2.4 μm范圍的高光譜地面場景進行了仿真。根據(jù)結果來看，本文采用的方法實現(xiàn)了對預定場景在光譜尺度的仿真。但是本文的方法忽略了地面幾何和輻射能量入射角度等因素造成的輻亮度變化情況，還具有一定的改進空間。