反演理论实验

采用三层H型模型进行直流电测深正演计算。模型参数为:第一层电阻率为100Ω·m，第二层电阻率为10Ω·m，第三层电阻率为300Ω·m，各层厚度均为5 m。正演计算所得理论电测深曲线如图3.2所示，将理论电测深曲线作为实测数据。

图3.2 理论视电阻率及初始模型的确定

用最小二乘法反演程序对图3.2电测深曲线进行反演。根据图3.2将初始模型选择为三层大地，利用勘探深度为AB/5选择层厚度，利用曲线的特征分层及设置层电阻率。

从图3.2可知为三层大地，作两条垂线把大地分为三层。第一条垂线对应AB/2为7m。根据勘探深度为AB/5的原则把第一层厚度定为h1=d1=14/5=2.8m，第二层深度d2=100m/5=20m，所以第二层厚度为h2=d2-d1=20m-2.8m=17.2m。从视电阻率曲线的首支可知第一层电阻率为100Ω·m从尾支可知第三层电阻率为300Ω·m中间层电阻率定为极小值40Ω·m。建立如下初始模型文本文件csmx.dat:

模型层数

3

单个测深点数据个数

20

各层电阻率

100

40

300

各层厚度

2.8

17.2

是否约束反演(0=否，1=是)

0

约束层厚度数 0

层号 厚度

约束层电阻率数 0

层号 电阻率

如果已知第二层电阻率，可进行约束反演。反演结果如图3.3、图3.4所示。

图3.3 无约束反演结果

图3.4 有约束反演结果

图3.3为无约束反演结果，给定反演中止精度为0.01，共迭代4次，电测深曲线拟合误差为0.6%。反演模型参数见表3.2。

表3.2 无约束反演结果

图3.4为有约束反演结果，给定反演中止精度为0.01，共迭代6次，电测深曲线拟合误差为0.86%。反演模型参数见表3.3。

表3.3 约束反演结果

从反演结果图3.3可以看出，当没有先验信息时，由于S2等值现象，反演的模型与理论模型误差较大。而从图3.4可以看出，加入了先验信息进行约束反演的模型与理论模型更接近，但图中所示电测深曲线的拟合误差大于无约束反演。可见不能把电测深曲线拟合精度作为反演好坏的唯一标准，一定要注重标本的采集、先验信息的收集，这样才能获得更接近实际的反演结果。

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

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反射率比较

数值模拟反问题常常转化为优化问题，函数优化就是求一个函数的最优值以及达到该最优值的最优点，而最优化算法本质上是一个最优值的搜索过程。经典的优化算法如牛顿法、单纯形法、共轭方向法、最速下降法和罚函数法等，一般对目标函数要求连续、可微甚至于高阶可微、单峰等；需要对函数求一阶、二阶导数；受初值影响较大，算法容易陷入局部最小值，对于多峰函数优化问题具有较大局限性。

20世纪80年代初期以来，地下水水流与溶质迁移模型和数值优化方法相结合越来越普遍，目前常用的主要有以下两种方法。

3.4.7.1 数学规划方法

主要包括线性规划（LP），该方法广泛应用于线性目标函数及流量约束的地下水管理问题，解线性规划的软件主要有AQMAN，MODMAN，MODOFC，MODFLIP；非线性规划（NLP）；混合整数线性规划（MILP）；混合整数非线性规划（MINLP）。其中线性规划法计算效率较高，但仅适用于承压含水层，通常不能有效地处理溶质运移问题。非线性规划与动态规划的应用较广泛，计算效率上有优势，但需要计算目标函数对决策变量的导数即梯度，因此，该方法又被称为梯度法，在目标函数很复杂，而且为非线性时，结果往往会陷于一个局部最优解而不能识别全局最优解。

3.4.7.2 全局优化方法

主要以启发式搜索技术为根据的一类优化方法，包括模拟退火法、遗传算法、禁忌搜索法、人工神经网络法、外围近似法等，这些方法有识别全局或接近全局范围内最优解的能力。全局优化法能够模仿一定的自然系统，通常计算量很大。本书主要介绍4种现阶段应用广泛发展较为迅速的优化算法。

遗传算法（Genetic Algorithms，GA）是一类借鉴生物界自然选择（Natural Selection）和自然遗传机制的随机搜索算法（Random Searching Algorithms），求解问题一般包括编码、计算适应度、选择、交叉、变异、循环回到计算适应度，反复进行直到满足终止条件。该算法是处理一般非线性数学模型优化的一类新的优化方法，对模型是否线性、连续、可微等不作限制，也较少受优化变量数目和约束条件的束缚，其本质是一种高效、并行、全局搜索的方法，能在搜索过程中自动获取和积累相关搜索空间的知识，并自适应地控制搜索过程以求得最优解。目前已广泛用于函数优化、参数辨识、机器学习、神经网络训练、结构设计和模糊逻辑系统等方面。常用的GA计算程序有MGO（Modular Groundwater Optimizer），模块化地下水优化程序，该程序是地下水水质管理的通用优化模型。将水流和迁移模拟程序与遗传算法相结合，能适应非线性复杂目标函数，能够处理水头、梯度、水流以及浓度等约束条件。SOMOS程序，实现了包括遗传算法和人工神经网络的优化算法，能处理经济、环境以及地下水管理体积等问题，同时SOMOS可以将MODFLOW和MT3DMS作为模型的组成部分进行运算。但是目前遗传算法的应用还存在明显的不足，主要表现为以下几点：

1）GA的算法设计和关键控制参数选择对优化性能的影响明显，直接影响算法的搜索效率和优化性能，甚至导致“早熟”收敛；

2）参数识别研究中的编码方案以二进制编码为主，计算量和存储量大。

人工神经网络（Artificial Neural Network，ANN）是由大量神经元通过极其丰富和完善的联结而构成的自适应非线性动态系统，它使用大量简单的相连的人工神经元来模仿生物神经网络的能力，从外界环境或其他神经元获得资讯，同时加以简单的运算，将结果输出到外界或其他人工神经元。神经网络在输入资讯的影响下进入一定状态，由于神经元之间相互联系以及神经元本身的动力学特性，这种外界刺激的兴奋模式会自动地迅速演变成新的平衡状态。人工神经网络是一种计算系统，包括软件与硬件，它使用大量简单相连的人工神经元来模仿生物神经网络的能力。人工神经网络是生物神经元的简单模拟，它从外界环境或者其他神经元取得资讯，同时加以非常简单的运算，输出其结果到外界环境或者其他人工神经元。人工神经网络系统反映了人脑功能的许多基本特性，但它并不是人脑神经系统的真实写照，而只是对其作某种简化、抽象和模拟，这也是当前的现实情况。是目前对人脑神经及其智能机理的研究水平所能做到的，对人脑智能机理的简化、抽象和模拟是人工神经网络研究的基本出发点。

支持向量机是基于统计学理论的VC维理论和结构风险最小化原理而提出的一种新的机器学习方法。与传统的神经网络学习方法相比，支持向量机从结构风险最小化原则出发，求解的是一个二次规划问题而得到全局最优解，有效地解决了模型选择与过学习问题、非线性和维数灾难以及局部极小等问题，在解决小样本、非线性、高维模式识别问题中表现出许多特有的优势。

模拟退火算法是对固体退火过程的模拟。在金属热加工工艺中，将金属材料加热到某一高温状态后，让其慢慢冷却，随着温度的降低，物质的能量将逐渐趋近于一个较低的状态，并最终达到某种平衡。模拟退火算法是基于金属退火的机理而建立的一种全局最优化方法，它能够以随机搜索技术从概率的意义上找出目标函数的全局最小点。模拟退火算法的主要缺点是解的质量与求解时间之间存在矛盾，该算法对于多应力期模型和大量水文地质参数的反演，收敛缓慢，得不到满意的结果。

---