反射率反演

遥感器接收目标辐射或反射的电磁波所形成的遥感原始图像与目标相比是失真的，这是因为在太阳-大气-目标-大气-遥感器的光线传播路径中，许多因素的影响造成接收的信号不能准确的反映地表物理特征。这些因素归结为以下几个方面：

1）大气内容物的影响。大气主要由大气分子和气溶胶组成，这两者的影响行为是不相同的。大气分子瑞利散射、气溶胶的Mie散射；大气分子与气溶胶的吸收及两者的汪盯做耦合作用。一方面，大气的吸收导致消光，减少了辐射量，降低了图像对比度，使图像变得暗淡；另一方面，大气散射导致的程辐射，增加了辐射量；

图5.16 系统级几何校正效果（660nm波段）

2）表面因素的贡献。在一般的应用中，为了简化计算，假定地表为朗伯体，反射与方向无关。事实上则誉任何物体表面在物理特性与物质结构上都不是理想朗伯体，因此认为地面是朗伯体会带来误差，而当地表方向反射特性突出时，假设地面是朗伯体的大气纠正方法精度受到限制。另一个因素是，由于大气散射的存在，邻近像元的反射光也会进入目标视场从而影响辐射量，即交叉辐射。

3）地形因素的影响。目标高度与坡向会对辐射造成影响。

4）太阳辐射光谱的影响。太阳本身是一个黑体，其光谱辐射按照普朗克定律有一定的形状，这个因素在反射率反演中需要予以考虑。

由以上可知，大气对光学遥感的影响是十分复杂的。为此，学者们尝试着提出不同的大气纠正模型来模拟大气的影响。但是对于任一幅图像，其对应的大气数据几乎是永远变化的，且难以获得，因而应用完整的模型纠正每个像元是不可能的。最早的大气纠正方法是从图像本身来估计大气参数，反复运用大气模拟模型进行纠正。结合地面实况数据进行大气校正是另一种方法，其包括两种类型：一种是通过地面测定大气参数（如可见光近红外的气溶胶的密度及红外区域的水汽浓度），再结合辐射传输方程作近似求解；另一种是测得地面目标物的反射率，再与图像数据进行比较来消除大气的影响。地面同步测量有助于提高精度，但是却需要人力物力，且应用区域也有限。此外还有一些大气纠正的方法。例如在同一平台上，除了安装获取目标图像的遥感器以外，也安装上专门测量大气参数的遥感器，利用这些数据进行大气校正。

综上，大困衡气纠正具体算法大致可归纳为基于图像特征的相对校正法、基于地面的线性回归模型法、基于大气辐射传输模型法和复合模型法这四种。

基于图像特征的相对校正法是在没有条件进行地面同步测量的情况下，借用统计方法进行图像相对反射率转换。从理论上来讲，基于图像特征的大气校正方法都不需要进行实际地面光谱及大气环境参数的测量，而是直接从图像特征本身出发消除大气影响，进行反射率反演，基本属于数据归一化的范畴。精确的大气校正需要精确的测量大气参数和复杂的运算，这些在许多遥感应用中，往往很难满足。并且，在某些应用中不一定需要绝对的辐射校正。此时，这种基于图像的相对校正就能满足其要求。辐射校正的统计模型主要有内部平均法、平场域法、对数残差法等。

基于地面线性回归经验模型法是一个比较简便的定标算法，国内外已多次成功地利用该模型进行遥感定标实验。它首先假设地面目标的反射率与遥感器探测的信号之间具有线性关系，通过获取遥感影像上特定地物的灰度值及其成像时相应的地面目标反射光谱的测量值，建立两者之间的线性回归方程式，在此基础上对整幅遥感影像进行辐射校正。该方法数学和物理意义明确，计算简单，但必须以大量野外光谱测量为前提，因此成本较高，对野外工作依赖性强，且对地面定标点的要求比较严格。这种方法仅适用于地面实况数据特定的地区及时间。

大气辐射传输模型能够较合理地描述大气散射、大气吸收、发射等过程，且能产生连续光谱，避免光谱反演的较大误差，因而得到了最广泛的应用。在遥感实际应用中，大气辐射传输模型需要进一步简化，如：假定大气是水平均匀的、假定地表是朗伯体、排除云的存在、运用各种条件下的标准大气模式及气溶胶模式（由于长期试验数据积累和理论研究归纳而成）等。不同的模型的假定也是有些区别的，比如6S是地表均匀、非朗伯体的模型而5S是地表均匀、朗伯体的模型。

5.5.1 原理与方法

基于图像特征的相对校正法主要有内部平均法、平场域法、对数残差法等。

（1）内部平均法

假定一幅图像内部的地物充分混杂，整幅图像的平均光谱基本代表了大气影响下的太阳光谱信息。因而，把图像DN值与整幅图像的平均辐射光谱值的比值确定为相对反射率光谱，即

ρλ ＝ Rλ /F （5.14）

式中：Rλ为像元在该波段的辐射值；Fλ为整幅图像的平均辐射光谱值；ρλ为该像元的相对反射率。

（2）平场域法

平场域法是选择图像中一块面积大且亮度高而光谱响应曲线变化平缓的区域，利用其平均光谱辐射值来模拟飞行时大气条件下的太阳光谱。将每个像元的DN值与该平均光谱辐射值的比值作为地表反射率，以此来消除大气的影响。

ρλ ＝ Rλ /Fλ （5.15）

式中：Rλ为像元在该波段的辐射值；Fλ为平场域的平均辐射光谱值；ρλ为该像元的相对反射率。

利用平场域消除大气影响并建立反射率光谱图像有两个重要的假设条件：一个是平场域自身的平均光谱没有明显的吸收特征；另一个是平场域辐射光谱主要反映的是当时大气条件下的太阳光谱。

平场域模型已广泛应用于遥感数据处理中，它是在内部平均法模型基础上发展起来的，这种模型克服了内部平均法模型易受像幅内吸收特征影响而出现假反射峰的弱点，而且计算量更小，其不足之处在于选取光谱地理平台单元时，会引入人为的误差，而且需要对研究区内地物光谱有一定的先验了解，当选取具有不同反射率等级的地理平台单元时，会引出不同处理结果。当研究区位于山区或其他地形起伏较大的复杂地区时，选择地理平台单元较为困难。

（3）对数残差法

对数残差法的意义是为了消除光照及地形因子的影响。按照一定的规则调节每个像元值，使其在每一个被选定的波段上的值等于整个图幅的最大值，然后对每一个波段减去其归一化后的平均值。假设有

DNij ＝ TiRijIj （5.16）

式中：DNij为像元i的j波段的灰度值；Ti 为像元i处表征表面变化的地貌因子，对确定的像元所有的波段该值都相同；Rij为像元i波段j的反射率；Ij为波段j的光照因子。

由表5.2我们可以看出，以上三种方法中，只有残差图像法是真正意义上的辐射校正。

表5.2 高光谱基于图像特征的相对校正法对各种影响辐射的物理因素的补偿能力的比较

除上述基于图像特征的相对反射率校正外，还可基于大气辐射传输理论的大气纠正模型开展反射率校正工作。

5.5.2 实例分析

（1）基于6S的反射率反演模型对CHRIS数据进行反射率反演实例

欧空局（European Space Agency，简称ESA）的Proba（Project for On-Board Autonomy）卫星于2001年10月发射成功，是星上自主运行技术的示范，也是新的航天（包括硬、软件）技术的成功范例，可用于测试地球观测和空间环境仪器性能。有效载荷包括一台紧密型高分辨率高光谱成像仪（CHRIS）和一台辐射测量传感器（SREM）及岩屑探测器（DEBIE）和宽视场角地球定位相机和恒星跟踪器及陀螺仪。

经过一年的运行Proba已经完成了它的技术示范任务，它为科学界提供了前所未有的创新性的卫星高光谱多角度CHRIS数据。CHRIS图像提供了可见/近红外高空间、高光谱分辨率的地表反射率数据，利用Proba的定位功能，可以得到试验区五个观测方向上的二向反射率数据（BRDF），五个角度的观测几何见图5.17。CHRIS的主要参数见表5.3：

图5.17 CHRIS/Proba图像获取几何示意

C1～C5为相应的中心时间

表5.3 CHRIS/PROBA 的主要技术参数

CHRIS有五种工作模式，其中模式3和模式5是为陆地应用设计的，模式2则应用于水体，模式4应用于叶绿素反演与监测。模式3可获取18个波段，相应图像空间分辨率为17m，不包含水汽通道；模式5可获取37 个波段，空间分辨率为34m，包含940nm的水汽通道，可用于水汽反演。图5.18为模式5的每个扫描行的组成。

图5.18 CHRIS模式5每扫描行像素组成

我们在本试验中获取的模式5观测天顶角为0°的数据，相应的波段信息见表5.4。

表5.4 PROBA/CHRIS 工作模式5（mode5）对应的波长信息

续表

实验所采用CHRIS数据其他信息描述：

获取方式：MODE5；

波段数：37 个，波长范围 442.49～1025.30nm，包括一个水汽通道（波段 31）：中心波长945.31 nm；

空间分辨率：34 m；

图像行列数：766列× 748行；

图像数据类型：长整型；

图像中心点经纬度：116°52′E，40°10′N（昌平一带）；

图像获取时间：2004年7月8日，3：22（UTC时间）；

卫星平台高度：596 km；

图像物理单位：μW/（m2 ·nm·sr）；

地面平均高程：200 m。

图像处理：

1）去坏行处理，以相邻两行（每侧各两行）取平均，代替坏行。

2）根据反射率反演软件的要求，即图像为整型数据和图像定标后辐射单位为W/（m2 ·μm·sr）将图像单位μW/（m2·μm·sr）转换为W/（m2·μm·sr），从量纲上来看，前者是后者的1/1000。因此根据原图像的数据范围，除以10取整得到整型数据（短整型），然后将所有波段的增益均设为0.01，将得到辐射单位为W/（m2·μm·sr）的整型图像数据。

3）启动反射率反演软件，设置各项参量，运行程序。程序输入参数界面如图5.19所示。

图5.19 CHRIS/Proba反射率反演输入界面

4）由于传感器自身光谱与辐射定标的精度直接制约着反射率转换的可靠性，为了有效去除图像数据和大气辐射传输模型间存在辐射定标不匹配现象需要进行图像反射率光谱去噪平滑。

结果表明：反演得到的反射率在498～760 nm波长区间能够表征植被（玉米）的反射率光谱特征，与图5.20（c）相比，能够去除绝大多数的大气吸收特征，但在760～805nm之间的峰形与标准植被光谱差异较大，这可能与CHRIS仪器本身在760 nm附近的氧气吸收带的光谱定标误差有关。在805 nm以后与被标准植被反射率曲线差异也很大，主要是近红外的高反射率“平台”不明显，反而呈急剧下降趋势，940 nm附近的水汽吸收带也没有反映；对于土壤光谱，众所周知，常见的土壤光谱反射率在＜1140nm波长范围内呈现单调增加的趋势，而图5.20（b）中的土壤光谱反射率在900nm之后递减，这与常识相违背，而事实上即便在土壤的野外光谱测量上940 nm的水汽吸收作用也并不明显。究其原因可能在于两方面：一是CHRIS仪器本身的定标精度，另一方面也与反射率反演模型的校正误差有关。

图5.21是将通常的统计方法IRAA和FF应用于CHRIS图像上得到的同样采样点上的植被和反射率光谱，可以看到，两种方法得到的植被反射率光谱在谱形上非常相似，FF方法得到的反射率更平滑些，但二者在以500 nm为转折，反射率先下降后增大，这与通常的植被在蓝光波长范围的单一上升趋势不一致，与反射率反演结果相比，没有760 nm的凸起变形，这从另一个侧面反映了基于模型的反射率反演对定标要求更为苛刻；从土壤光谱来看，两种统计方法获得的结果差异很大，特别是在＜750 nm波长区域，IRAA起伏变化剧烈，而FF结果总体上要平缓得多，呈平稳微小上升趋势，这与土壤光谱反射率在＜1140 nm呈单调递增趋势的常识相一致，虽然平场域法的有效性与“平场区域”的选择恰当与否直接有关，但无疑在本试验中平场域法得出的结果最切合实际。

图5.20 CHRIS数据反射率反演结果

图5.21 内部平均法（IRAA）和平场域法（FF）得到的CHRIS图像反射率光谱

图5.22 地面实测光谱采样到与CHRIS波长设置相一致

为比较三种方法所得地物反射率光谱与地面实测光谱间谱形上的相似性，我们收集到2004年7月6日本试验区内的地面测量光谱，典型地物为玉米和土壤（裸土）。测量时间为北京时间11：40左右，与图像获取准同步，便于与图像光谱进行比对。光谱测量采用的仪器是ASD Fieldspec FR2500光谱仪，其光谱范围为350～2500nm，采样间隔为1.4nm（350～1000nm 区间）和2nm（1000～2500 nm区间）。每个样本测量10 次取平均作为最终光谱，以避免随机噪声干扰。图5.22 为根据CHRIS中心波长和半高全宽（FWHM）采样后的光谱。

计算得出它们之间的相关性（表5.5），并对可见光和近红外分别进行比较。可以看到，对玉米光谱而言，反射率反演算法订正后的反射率光谱与实测光谱间的一致性最好，特别是在可见光范围的一致性远远高于经验方法；对土壤光谱而言，在可见光范围，基于反射率反演的仍保持相似性最高，但是在近红外波长范围，反射率反演和IRAA都与实测光谱呈负相关，基于模型反射率反演不能很好地表征植被在近红外反射率“平台”（750～900nm）和900～1100nm的水汽吸收特征，特别是平台部分相关系数为-0.43221。相比之下只有FF方法在可见-近红外波段都保持较高的相似性。也印证了上文的分析结论。

表5.5 三种方法得到的植被反射率光谱在可见光区域的相关性（R）

为此，将反射率反演纠正结果与FF相结合，即保留760 nm之前的反射率反演光谱，将760 nm之后的FF光谱做适当平移，然后采用经验平场反射率转换算法（EFFORT）对光谱做进一步平滑处理，可以得到与真实光谱更加一致的光谱：谱形的相似性和特征位置的保持。修正后的玉米光谱见图5.23。

（2）基于MODTRAN的反射率反演模型对Hyperion数据进行反射率反演实例

启动基于MODTRAN的反射率反演模型，其界面如图5.24所示。

美国航天局（National Aeronautics and Space Administration，简称NASA）的EO-1（Earth Observing One）卫星于2000年11月发射成功，其上搭载的Hyperion高光谱成像仪目前已获取了大量高质量的星载高光谱数据（表5.6）。

使用如图5.25所示庆阳地区Hyperion数据进行反射率反演，提取裸土像元反射率光谱，并将其与ASD数据进行比较，如图5.26 所示。由图可见，在可见近红外波段，两光谱在波形及量值上相近，在短波红外波段，反射率反演结果要略低于ASD采集的数据，这主要是因为气溶胶数据不准确造成的。反射率反演的裸土反射率光谱与 ASD 采集的裸土反射率光谱之间的光谱相关系数达到0.9342。

图5.23 与FF 相结合修正后的玉米光谱与真实光谱比较

图5.24 基于MODTRAN的反射率反演模型界面

图5.25 实验所用Hyperion数据

表5.6 Hyperion/EO-1 的主要技术参数

图5.26 反射率反演得到的反射率与ASD反射率比较

图5.27 敦煌实验场Hyperion数据

使用如图5.27所示敦煌实验场地区的Hyperion数据进行反射率反演，提取实验场均一像元的反射率光谱，并将其与准同步的ASD数据进行比较，如图5.28所示。由于该ASD数据在大于1800 nm的谱段噪声非常大，所以只比较450～1800 nm之间的谱段。两光谱在波形及量值上相近，反射率反演的反射率光谱与ASD采集的反射率光谱之间的光谱相关系数达到0.9516。

图5.28 反射率反演得到的反射率与ASD反射率比较

以下是对Visual Modflow和FEFLOW区别的一些浅薄的认识，和大家共同探讨。这里的Visual Modflow指的是目前大家常用的版本，不含Visual Modflow中MODFLOW－SURFACT模块。

Visual Modflow和FEFLOW都是当前世界上十分流行的可用于模拟三维地下水流和溶质运移模拟评价专业软件系统，属商业软件。它们都具有直观的、强有力的图形交互界面，模型剖分、输入参数和模拟结果，都可以用图形显示，并支持三维可视化和，做到了真正的人机对话，在许多行业和部门内得到了广泛的应用。但同时，它们之间却存在着许多差别，各有千秋：

（1）从软件功能上看，FEFLOW要比Visual MODFLOW更为全面一些，FEFLOW除了可以模拟Visual Modflow所能模拟的二维、三维饱和流状态的水流和溶质运移问题之外，还可以模拟多层自由表面含水系（包括滞水模拟）、热转递、可变密度流场（盐水或海水入侵问题）以及非饱和带流场及物质运移问题。

（2）从数值法的计算原理上，Visual Modflow采用的是有限差分法，而FEFLOW采用的是有限单元法。

（3）从离散化方面来说猛姿，由于Visual Modflow采用的是有限差分法，所以对所模拟的地质体采用矩形网格进行剖分，这种网格的优点在于，用户易于准备数据文件，便于输入文件的规范化，但是，当需要在所关注的地点附近（比如井附近）要增加计算单元的密度，就必须同时对经过该点附近区域的所有的行和列都进行加密，这样使得计算量大大增加。而FEFLOW不存在这个问题，由于其采用的是有限单元法，其剖分单元的形状可以是灵活多变（可以是三角形、也可以是矩形），一般采用三角形剖分，加密的时候，可以只对感兴趣的地方加密，相比Visual Modflow来讲，减少了运算量。用三角形剖分的另一个好处是，在刻画模拟区的外部边界时，可以利用三角形的边很好地控制外边界范围，这样刻画出来的边界比用Visual Modflow刻画出的外部边界要精确。以上两点通过两个软件的实例便可看出。另外，有限差分法对于处理复杂地质体中的地下水三维渗流场模拟方面存在着不足，没有有限元三角剖分灵活多变。拿知厅

（4）FEFLOW具备地理信息系统数据接口，可以充分利用已有的 ARC/INFO GIS地理信息系统数据产生有限单元网，设置边界条件和参数。

（5）Visual Modflow采用的是模块化结构，在软件输入 *** 作过程便体现了这一点，如边界条件里面的定水头边界、河流边界、截渗墙边界、排水沟边界、补给边界和蒸发边界等，而FEFLOW的边界条件是按照一类、二类、三类和井流边界划分的。这两种形式各有优缺点，模块化结构对常见的几类边界进行了分类，用户可根据问题直接选择该边界模块进行输入编辑 *** 作，十分方便，但遇到特殊水文地质问题时就显得不足；FEFLOW中边界问题的分类采取了广义的边界条件分类，因此在处理水文地质边界条件时就非常的灵活，但是这种过于集中的输入方式也给输入工作带来了不便，比如源汇项的输入过于集中，用户需要对其数据进行整理或者预处理才能输入；另外，在非稳定流模拟中，对于蒸发输入的 *** 作在Visual Modflow中可以很容易实现，但是在FEFLOW中却要通过FEFLOW的二次开发工具IFM模块编程来实现。

（6）在混合井的模拟方面二者都存在不同程度上的不足。混合井流是在生产过程中十分常见的一种地下水开采方式。但混合井流的模拟一直是MODFLOW的一个缺陷。尽管MODFLOW建议“多层井的流量必须以某种形式人消隐为地分配给每一单层，……把井流量按每一层的导水系数大小分配，即 ”，其中 和 分别为第 层流量和总（井口）流量，Ti和求和公式T分别为第i层导水系数和总导水系数。但实质上，这种方法是不妥的，它不是模拟，而是“处理”，一种与机理不符的“处理”。因此这个问题应引起地下水流数值模拟工作者的重视。FEFLOW以前的版本中在处理该问题时，也有类似的不足，其处理方法是在同一点上布置多口井（每个井开采的层位不同）来实现混合井开采问题，不过在新的Version5.3版本中已有改进,是否完善还有待进一步验证。（7）干涸单元问题。在Visual Modflow的计算过程中，如果计算水头低于该计算单元之地面标高，该计算单元则处于非饱和状态。这时，Visual Modflow并不考虑该计算单元在非饱和态下的渗透系数，而是将这个计算单元列为“干涸计算单元”（Dry cell），并将其导水系数赋值为零。一旦一个计算单元变为干涸单元，他就将被重新定义为不透水或非活动单元计算，从而被摒除出以后的模型计算。计算结果也可能因为干涸计算单元的出现而受到影响。所以，Visual Modflow中对于干枯单元的处理方法还不够完善。为了克服88年版MODFLOW中计算单元“一旦变干就永不能恢复”的问题，美国地调局的McDonald等人于1992年推出了一个新的子程序包，称为BCF2（Block Centered Flow Package)，允许干涸单元重新成为有效的计算单元，即出现了重新变湿（Rewetting）选项，但是其结果并不理想，虽然干涸单元格可以得已重新变湿，但是又造成了运算不易收敛的问题。

而在FEFLOW中，由于其没有采用这种处理方式，也就不存在这一问题。

（8）模型文件的保存形式不同。在Visual Modflow中，一个模型往往生成很多文件，模型构成数据及计算结果分别被保存成不同的文件类型，如各层标高文件保存为VMG文件，边界及参数保存在VMP文件中，井文件保存为VMW文件，输出结果水头数据保存为HDS文件，降深数据保存为DDN文件，水均衡数据保存为ZOT文件等，你可以直接对这些文件 *** 作，或者利用这些文件构建模型等。在FEFLOW中，一个模拟问题全部保存在一个fem文件内，模拟结果也保存为一个文件（DAC文件）。

（9）在利用Visual Modflow模拟非稳定流问题时，模拟计算可以随时停止、暂停，但若问题没有模拟结束就无法显示模拟结果。而在FEFLOW中，非稳定流模拟计算可以随时暂停，以便用户显示和分析中间模拟结果，并且工作窗口可以实时显示地下水非稳流场，温度场及污染物迁移模拟结果。

以上是对两软件的一些比较，大家可以根据自己需要学习或者选用这两款软件，当然，还有其它类似的一些软件，如Processing Modflow、GMS、Visual Groundwater等等。但是，需要注意的是，软件只是我们解决模拟问题的工具而以，最重要的还是对基础理论把握和对实际问题的正确认识，只有正确地认识和分析各类水文地质问题，才能做好水文地质模拟工作。

---