,,

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

0 引言

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

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

1 大氣輻亮度傳輸模型

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

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

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

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

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

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

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

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

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

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

圖1 大氣輻射傳輸路徑

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

(1)

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

(2)

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

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

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

2 鄰近像元反射率計(jì)算

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

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

圖2 光子傳輸路徑

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

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

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

(3)

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

P=e-l

(4)

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

(5)

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

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

(6)

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

(7)

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

3 大氣傳輸計(jì)算

3.1 大氣傳輸影響因素

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

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

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

3.2 MODTRAN軟件

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

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

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

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

其中，散射輻亮度模式輸出結(jié)果包含表1所示內(nèi)容。

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

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

4 仿真結(jié)果

4.1 地面場景仿真

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

圖3 地面場景分布情況

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

圖4 不同地物反射率曲線

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

圖5 波長為0.5 μm時(shí)地面反射率

圖6 波長為1.0 μm時(shí)地面反射率

4.2 鄰近像元反射率仿真結(jié)果

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

圖7 蒙特卡洛仿真點(diǎn)擴(kuò)散函數(shù)流程

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

圖8 光子密度圖

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

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

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

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

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

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

(8)

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

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

圖10 A、B、La、S仿真結(jié)果

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

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

圖11 波長0.49 μm時(shí)入瞳處輻亮度

圖12 波長0.49 μm時(shí)輻亮度灰度圖

從圖中可以看出，雖然混凝土道路上像元的入瞳處輻亮度仍然大于兩邊植被對(duì)應(yīng)的像元輻亮度，但與反射率數(shù)據(jù)相比，出現(xiàn)了明顯弱化，這主要是因?yàn)樵诖髿鈧鬏斶^程中的衰減造成的。

5 總結(jié)

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