高光谱场景辐亮度模型仿真研究

,,

(1.解放军航天工程大学 研究生院,北京 101416；2.解放军航天工程大学 电子与光学工程系,北京 101416)

0 引言

星载高光谱遥感的射传输路径是：太阳光谱辐射经过大气反射、散射、吸收等影响后，到达地表；在被地物反射后，其反射辐射再次经过大气影响，并与大气直接反射的太阳辐射一起构成传感器入瞳前的辐亮度数据，进入传感器中的辐射能量通过传感器的辐射维、空间维、光谱维离散采样后输出数字图像。这一过程可以分为场景系统和传感器系统。其中，场景系统作为整个高光谱遥感体系中影响因素最复杂多变的部分，尚未建立完备的模型来全面分析。

目前主要有两种方法对地面场景高光谱特性仿真。第一种是人工描绘场景中各物体的三维形貌，并利用材质光谱结合纹理图设定各物体的反射特性[1]，第2种利用高分辨率的高光谱数据描述场景反射率分布[2]。对于分析单一因素对高光谱成像的影响能力来说，以上两种方法显得过于复杂。本文以传感器入瞳前的大气辐亮度传输模型为依据，通过简化场景的几何形态，并将其分为地面系统和大气系统两部分，利用地面反射率数据和MODTRAN大气传输软件，实现对高光谱场景的仿真模拟。

1 大气辐亮度传输模型

多年来，研究者根据辐射在地表与大气之间的传输过程出发，针对不同的地表特性，主要建立了以下5种典型的大气辐射传输模型[3]：

1)不考虑地-气多次反弹的模型

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

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

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

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

以上5种模型的理想化程度逐渐降低，计算精度和要考虑的参数也越来越复杂。

在高光谱遥感成像过程中，太阳光谱辐射经过大气反射、散射、吸收等影响后，到达地表；在被地物反射后，其反射辐射再次经过大气传输到达传感器入瞳处。

如图1所示，将入射到传感器前的总辐射能量分成三部分：

(1)来自太阳辐射经过大气传输到达目标表面并被反射至传感器前能量：图中路径a和路径b(路径a是太阳辐射经大气传输后直接照射地物目标，路径b是太阳辐射经大气多次散射后的辐射能量照射到地物目标后直接进入传感器)；(2)来自太阳辐射经过大气传输照射到背景景物并被反射至传感器前能量：图中路径c和路径d(路径c是太阳辐射经大气传输直接照射到背景景物，路径d是太阳辐射经过大气多次散射后照射到背景景物)；(3)路径e是大气路径程辐射，即太阳辐射被大气多次散射后，不经过地物场景反射，直接进入传感器视场内的能量。

图1 大气辐射传输路径

文献[4]总结了包含多次散射项的传感器入瞳处辐亮度。将地面看做是非均匀的，反射率为ρ(λ)的朗伯体时：

(1)

在单一角度的天顶观测情况下，把场景看作非均一的朗伯面的模型时，上式可以等效为[5-7]：

(2)

式(2)中，L为入瞳前传感器处像元接收的总辐亮度；ρ为目标像元表面反射率；ρe为邻近像元表面平均反射率；S为大气朝下的球面反照率；La为大气程辐射；A、B为基于大气条件和几何光照条件的系数。

式中所有的变量都为光谱波长的函数，为方便起见，波长指数被省略。等式右边第一部分代表太阳辐射经大气传输到地表反射直接进入传感器前的辐亮度；第二部分为背景景物像元经大气散射进入传感器前的辐亮度；第三部分即为大气程辐射。

等式中ρe需要通过大气点扩散函数对ρ处理得到，而A、B、S和La可通过大气传输模型MODTRAN设定观测视场角、平均海拔高度、大气模型、气溶胶类型、能见度范围等参数进行仿真反演[8]。

2 邻近像元反射率计算

被邻近像元反射的太阳光子通过与大气粒子的碰撞作用被遥感器接收，不考虑地气间多次反弹，即认为邻近像元对遥感器接收到的总辐亮度贡献率的空间分布函数等同于光子射向地面并与大气的随机碰撞过程，最终在地面形成的光子数的相对密度空间分布函数[9]。

蒙特卡洛方法[10]是一种基于“随机数”的计算方法，将所求解的问题同一定的概率模型相联系，用计算机实现统计模拟或抽样，以获得问题的近似解，故又称随机抽样法或统计试验法。因此，选择用蒙特卡洛方法来仿真大气碰撞的随机过程。

图2 光子传输路径

光子在大气中的传输路线如图2所示。其中z为垂直高度，x，y为水平位置坐标。光子初始位置为A点，目标位置为B点，B点在A点水平面方向的投影为B’。d为自由路程，即光子在连续两次碰撞之间经历的路程。θ是线段BA与z轴反方向的夹角，称为散射极角；Φ为BA与x轴正方向的夹角，称为散射方位角。散射方位角Φ是在[0，2π]内的随机数，散射极角θ是在[0，π]内的随机数。

由以上分析可知，求解光子的传输路径可等效为求解d、θ和Φ的值。

引入光学距离l的定义：沿自由路程对d削弱系数的积分。可由式(3)表示：

(3)

上式中，Ks为大气散射系数，Ka为大气吸收系数。l为无因次量，令非碰撞概率密度函数P与l的关系为：

P=e-l

(4)

则光子和大气粒子距离为l时的非碰撞概率可表示为：

(5)

其中：r1为区间[0,1]的随机数。当l趋近无穷时，r1趋近1，即相距无穷远的光子必然不发生碰撞，反之同理。则l可由式(6)求得：

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

(6)

大气光学厚度定义为大气削弱系数在垂直路径自上而下的积分：

(7)

对于水平均匀、垂直分布的大气，光学厚度τ只和垂直高度Z有关。对应的大气光学厚度之差：△τ=lcosθ；根据式(7)中光学厚度与垂直高度的关系，可求得自由路程d=△z/cosθ，光子的水平位移量分别为：△x=dsinθcosΦ和△y=dsinθsinΦ，其中Φ=2πr2，θ=πr3，r2、r3分别为[0，1]内的随机数。

3 大气传输计算

3.1 大气传输影响因素

地球表层的大气主要分布在距离地球表面约80～90 km的低层大气中，高层大气对太阳辐射传输几乎没有影响。底层大气中存在着多种气体以及微粒，如尘埃、烟、雾、风、雨、雪等，使光的传输特性发生变化。影响大气传输的最主要因素是大气分子和气溶胶的影响[11]。

大气分子主要包括N2、O2、CO2、H2、H2O等，其中占N2对远红外和微波会表现出强吸收；O2在可见光近红外波段具有0.76 μm、0.688 μm、0.628 μm3个强吸收波段。CO2、H2O是对可见光和近红外吸收的重要吸收分子。在太阳光谱范围内，2 μm、1.6 μm、1.4 μm是CO2的3个主要吸收通道；水汽对太阳辐射的吸收主要位于近红外及短波红外区，其吸收强度与大气压、温度和光谱波长紧密相关。除上述分子外，大气中还含有He、Ne、Ar、Xe、O3等少量分子，虽然在可见光和近红外区有较强的吸收谱线，但是含量很少，可不考虑吸收作用。

大气中除大气分子外，还存在尺寸较大的大量固态、液态微粒，包括尘埃、烟粒、微水滴、盐粒及有机微生物等，这些微粒在大气中的悬浮成溶胶状态，称为气溶胶。气溶胶的光学特性、粒子形态等特性会因为时空不同发生剧烈变化，深刻影响大气辐射传输过程。气溶胶的尺寸一般为0.001～10 μm。气溶胶一方面使太阳辐射在传输过程中发生散射，使返回大气层的辐射得到增强，有色气溶胶的吸收特性还会使到达地表的太阳辐射产生衰减；另一方面，已经到达大气层内的太阳辐射又难以被反射出大气层。大气中存在的气溶胶主要来自两方面：一是火山烟、风暴中的扬尘、植物的花粉、海水蒸发的盐粒等自然现象；而是诸如化石燃料燃烧、汽车尾气等人类活动。研究者把常见的气溶胶粒子按照不同比例混合来定义其类型，例如可用于干旱地区的沙漠型、用于受海洋季风影响严重的海洋型以及具有不同地物属性的大陆型[12]。

3.2 MODTRAN软件

MODTRAN是一个中分辨率大气辐射传输模型，由美国空军研究实验室(ARFL)与光谱科技公司利用FORTRAN语言开发。它把0至100 km的整层大气当成36个均一水平的薄层气体叠置，每薄层的光学特性通过温度、气压、气体分子的不同组合来刻画。内置了多种大气模型，能在0到50000 cm-1波数范围内以1.0 cm-1的分辨率进行辐射传输运算。MODTRAN利用二流近似、离散纵标等多次散射方法，计算包含大气吸收与多次散射作用的透过率、辐照度等。目前已更新至MODTRAN6版本，但由于知识产权原因，国内一般研究人员仍普遍使用MODTRAN4版本。

在查阅相关MODTRAN文献后，发现目前绝大部分以性能介绍及其在遥感图像处理中大气校正应用为主题。有关于对入瞳处辐亮度计算原理和详细实现过程的说明则非常少。本文利用MODTRAN在windows环境下运行的可视化界面软件PcModwin进行仿真分析。

通过多次实践以及查看输出文件MODOUT(包括Modout1、Modout2、Modout3，其中Modout1包含输出的完整结果)，初步总结了其关于辐亮度的计算过程：(1)选择内置的大气廓线数据存储于公共数据中；(2)将输入数据中确定的光谱分辨率和光谱范围等效为参与辐射传输计算的波数范围；(3)将输入数据的几何观测关系等效为由传感器位置、观测天顶角、方位角几何关系；(4)利用输入的气溶胶数据，利用气溶胶吸收和非对称系数，构建沿“地面到传感器高度”垂直路径的单次散射几何关系，计算单次散射路径气溶胶吸收辐射量；(5)在计算沿视线方向的单次散射几何关系基础上，得到多次散射参数；(6)根据气体分层理论，利用多次嵌套循环方法，计算整层大气的光谱辐亮度等参数；(7)重采样1cm-1波数的光谱辐亮度数据，计算结束。

MODTRAN按照计算目的不同可以分为透过率(Transmittance)、热辐射亮度(Thermal Radiance)、散射辐射亮度(Radiance with Scattering)、太阳直射照度(Direct Solar Irradiance)4种运行模式。

其中，散射辐亮度模式输出结果包含表1所示内容。

表1 MODTRAN散射辐亮度模式部分输出结果

利用MODTRAN的散射辐亮度模式可以获得在特定地表反射率状态下的传感器入瞳处辐亮度，而后通过相关数学逆推，得到大气传输参数A、B、S、La。

4 仿真结果

4.1 地面场景仿真

建立理想地面场景如图3所示，将该场景当作非均匀的朗伯体。设定传感器空间分辨率20 m，光谱分辨率10 nm，覆盖0.4～2.4 μm范围，轨道高度600 km并以推扫方式对地垂直成像。把地表按照传感器空间分辨率大小分割为若干面元，并假设每一面元的反射率均匀分布。

图3 地面场景分布情况

地面反射反射率数据ρ来自美国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 邻近像元反射率仿真结果

根据第2节所述，利用蒙特卡罗方法仿真大气点扩散函数的过程如图7所示。

图7 蒙特卡洛仿真点扩散函数流程

经过不断循环，直到光子到达地面，得到图8的光子分布密度图。

图8 光子密度图

蒙特卡洛法逼真描述了光子的传输过程，其结果是靠大量的统计数据的分析获得的。得到点扩散函数后，对4.1节中的反射率数据ρ进行平滑处理可以得到邻近像元反射率ρe。如图9所示，ρe(0.5)为波长0.5 μm时的邻近像元反射率情况。

图9 波长0.5μm时的邻近像元反射率

4.3 大气参数仿真计算

利用MODTRAN计算大气传输特性的参数。设置模式如表2所示，在运行MODTRAN时选择DISORT多次散射算法，通常设置散射算法流数为8，同时在地面反射率设置中，将地面反射率(Ground Albero)分别设置为0、0.5、1[14]，得到3种情况下总的TOT(传感器接收的总辐亮度)、PATH(多次散射程辐射)、TGR(总地面反射辐亮度)。

表2 包含多次散射的辐亮度计算步骤

通过MODTRAN仿真得到以上结果后，则关于式(1)的相关参数S、A、B、La可由式(8)仿真得到：

(8)

式中,PATH0是地面反射率为0的大气程辐射值；△TGR1是地面反射率分别为1和0时地面总反射辐亮度的差值；△PATH1是地面反射率分别为1和0时的大气程辐射差值；△TOT1是地面反射率分别为1和0时传感器入瞳前总辐亮度差值；△TOT0.5地面反射率分别为0.5和0时传感器入瞳前总辐亮度差值。

由式(8)计算的参数S、A、B以及程辐射La的结果为图10，需要注意的是，S中的峰值为算式中分母为0时出现的极化现象，上述求得的几类参数均为需求大气条件时的值，当大气条件改变时，反演得到的参数值随之改变。

图10 A、B、La、S仿真结果

4.4 地面场景的入瞳处辐亮度

利用4.1节、4.2节和4.3节计算结果，结合大气传输公式，就得到高光谱遥感入瞳处辐亮度数据。如图11、12所示为波长为0.49 μm时，地面场景的入瞳处辐亮度仿真结果。

图11 波长0.49 μm时入瞳处辐亮度

图12 波长0.49 μm时辐亮度灰度图

从图中可以看出，虽然混凝土道路上像元的入瞳处辐亮度仍然大于两边植被对应的像元辐亮度，但与反射率数据相比，出现了明显弱化，这主要是因为在大气传输过程中的衰减造成的。

5 总结

本文采取基于地表非均匀的朗伯面大气场景模型，利用蒙特卡罗算法仿真大气点扩散函数，得到了包含邻近效应的地表平均反射率；在对MODTRAN软件原理分析基础上，考虑多次散射效应，获得了传感器入瞳处的辐亮度数据。最后对给定的波长为0.4～2.4 μm范围的高光谱地面场景进行了仿真。根据结果来看，本文采用的方法实现了对预定场景在光谱尺度的仿真。但是本文的方法忽略了地面几何和辐射能量入射角度等因素造成的辐亮度变化情况，还具有一定的改进空间。