基于高光谱影像分解的土壤含水量反演技术

吴 见，刘民士，李伟涛

（滁州学院 地理信息与旅游学院，安徽 滁州239000）

土壤是农业生产的基础，土壤的水分含量、质地、有机质含量［1－2］等性质会影响作物生长，其中，土壤水分不仅是植物生长和发育的必要条件，也是精准农业中重要的信息。因此，土壤含水量监测一直是学者们研究的重点［3－4］。传统的土壤含水量测定方法，如张力计法、烘干法等均以点测量为基础，准确性较高，但周期长、工作量大。遥感技术主要以可见光—近红外、热红外和微波波段对土壤含水量进行反演，热红外波段反映的是地表很薄的一层信息，但也可利用热惯量和辐射平衡估测表层以下的土壤水分，高光谱遥感具有精细的光谱分辨率，能根据可见光和近红外波段的光谱信息，分析影响土壤含水量变化的光谱特征［5－6］，建立光谱特征与 土壤含水量的 关系模型［7－8］，进而定量反演土壤含水量。

目前，土壤含水量遥感监测研究大多是针对裸土进行的，对于卫星遥感来说，实际情况往往是植被与土壤同时存在于像元中，土壤含水量的监测不可避免地会受到植被光谱信息的影响。实际研究发现土壤含水量变化引起的光谱差异与植被含水量变化引起的光谱差异并不相同［9］，因此剔除植被光谱信息对土壤水分的干扰，是监测土壤含水量的有效途径。本研究对高光谱影像进行光谱分解，剔除植被光谱的干扰，同时对土壤贡献的光谱进一步处理，建立土壤含水量反演模型。

1 研究区概况及影像数据选取

1.1 研究区概况

怀柔区地处北京东北部，总面积达2 128.7km2，是华北平原、内蒙古高原以及燕山山脉的过渡地带。该区北部环山，南部是草原，地形包括深山、浅山、平原3种类型，从南向北绵延128km，且山地面积占整个怀柔区的88.7%，宜林山场森林面积达41%。地势呈北高南低，海拔在34～1 661m，变化较大。该区气候冬天受西伯利亚冷空气影响，寒冷干燥，夏天受海洋气团影响，温热湿润，春秋季节短且干旱多风。年平均气温在6～12℃，常年平均降水量在470～850mm。土壤主要包括草垫土、褐土、棕壤、风砂土等类型。

1.2 数据获取

本研究选取北京怀柔部分地区2001年5月19日的EO-1Hyperion高光谱数据，该数据共有242波段，其中1～70波段是可见光近红外波段（VNIR），71～242波段是短波红外波段（SWIR），光谱分辨率为10nm，空间分辨率为30m。

首先删除2个重复，20个受水汽影响严重及44个未定标波段，剩余176个波段；然后对剩余波段进行处理，包括坏线修复、条纹去除以及smile效应去除；最后，对处理后的图像进行检验，继续删除质量差的波段7个，剩余169个波段。被删除的具体波段为1～7，58～78，121～129，166～180，185～186，224～242。利用FLAASH软件对剩余的169个波段进行大气纠正，得到反射率图像。大气纠正后，采用1∶50 000地形图对影像进行几何纠正，总误差是0.35个像元。

2 基于植被光谱信息剔除的土壤含水量反演

2.1 基于光谱匹配的分解算法

目前，全受限的线性光谱混合分解模型是最为常用的方法，该模型中像元的光谱表达成各基本组分光谱及其面积比率的线性组合［10］，即：

式中：Rb——b波段的光谱信息；N——基本组分数目；fi——基本组分i占所在像元的面积比例；Ri，b——基本组分i在b波段的光谱信息；eb——b波段的误差。应用该模型的关键是尽可能减少模型中每个像元的误差，即使均方根误差RMSE最小化［10］：

式中：M——波段数。解该模型最常用的方法是最小二乘法，即求得式（2）RMSE最小的解。

将影像或光谱库中选取的基本组分光谱（r1，r2，…，rn）根据已设定的比例（f1，f2，…，fn）组合，获取若干测试光谱S1＝f1r1＋f2r2＋…＋fnrn，然后将待分解的像元光谱，即目标光谱S2与测试光谱S1进行匹配［10］，把匹配效果最佳的测试光谱的基本组分比例作为影像分解结果。因此，可将混合像元分解视为带约束的非线性最优问题［10］：

式中：G——目标函数，描述目标与测试光谱差异的匹配模型。由最小二乘法原理可定义G为测试与目标光谱间的欧式距离［10］：

因S1b－S2b＝eb，根据式（2）可知GLS＝ ■M ×RMSE，说明最小二乘法也属于光谱匹配的非线性最优化算法，但GLS采用欧式距离描述光谱匹配程度时，易受局部奇异值等因素影响。基于此，本研究引用杨伟等［10］提出的基于相关系数匹配的分解算法，设目标函数GSCM为：

式中：SCM——测试与目标光谱的相关系数。

采用相关系数SCM描述目标与测试光谱之间的匹配程度，以SCM最大值的测试光谱作为混合像元分解的结果。

2.2 土壤光谱信息处理

高光谱混合像元分解后，得到土壤和植被比例，以及裸土和纯植被基本组分的光谱信息，因此可进一步剔除植被光谱的干扰，计算各像元中土壤贡献的反射率信息。获取土壤贡献的反射率信息后，为突出土壤光谱的吸收特征，本研究采用光谱包络线去除［11］和一阶微分［12］两种方法对光谱曲线做进一步处理。

2.2.1 反射率光谱信息 由于不同像元对应的土壤含水量不同，在不同含水量条件下裸土基本组分光谱是有差异的，因此不能采用土壤比例信息与裸土基本组分光谱直接计算不同比例土壤贡献的反射率信息。由于研究区内植被类型较为单一，因此设定研究区域纯植被基本组分光谱不变，从原始影像混合光谱中将植被光谱剔除，进而获取土壤贡献的反射率光谱信息，由像元二分模型［13］：

可知，剔除植被光谱后的土壤光谱信息为：

式中：Rsb，Rv——土壤和植被在波段b贡献的光谱信息；fs，fv——土壤和植被在像元中的比例。

2.2.2 光谱一阶微分 一阶微分法是常用的光谱处理技术，不仅能够识别重叠光谱，增强曲线坡度上的光谱细微变化，使光谱吸收峰参数便于提取，而且可以有效消除系统误差，降低大气吸收、散射、辐射等噪声的影响［12］。一阶微分公式为［12］：

式中：λb——波段b；FDRλb——波段b和b＋1之间的一阶微分值；Rλb，Rλb＋1——波段b和b＋1的光谱反射率；Δλ——波段b和b＋1的波长之差。

2.2.3 包络线去除法 包络线去除法能够增强光谱曲线的反射、吸收特征，且反射率归一化在0～1.0之间，是一种有效突出目标吸收特征的光谱分析法，同时能够将光谱吸收特征归一化在同一背景上，方便不同光谱曲线之间的比较［11］。“包络线”即以直线逐点连接光谱曲线上的峰值点，且峰值点上的折线外角大于180°。光谱去包络是将原始曲线上的值除以包络线上对应的光谱值，计算公式为［11］：

式中：λb——波段b；Rcb——波段b的包络线去除值；Rb——波段b的反射率值；Rstart，Rend——光谱曲线起始、末端点的反射率值；λstart，λend——光谱曲线起始、末端点的波长；K——光谱曲线起始、末端波段间的斜率。

2.3 土壤含水量与光谱信息相关性分析

2.3.1 土壤含水量与反射率光谱相关性分析 从图1可知，反射率光谱与土壤含水量在450～1 050nm的相关系数为正值，相关系数的峰值分别出现在波段456～465，496～505，786～795，1 004～1 013nm处，相关系数值分别为0.47，0.48，0.37，0.29；在1 050～2 380nm的相关系数为负值，相关系数的峰值分别出现在波段1 134～1 143，1 324～1 333，1 467～1 476，2 012～2 021，2 302～2 311nm处，相关系数值分别为－0.42，－0.38，－0.51，－0.61，－0.56。为避免“过度拟合”现象及共线性问题，并非全部选取相关系数值最大的波段，而是根据光谱曲线特征划分区域，选取各区域中具有代表性或相关系数绝对值最大的波段作为敏感波段，即496～505，786～795，1 134～1 143，1 467～1 476，2 012～2 021nm处。

图1 土壤含水量与反射率光谱相关系数

2.3.2 土壤含水量与一阶微分光谱相关性分析 一阶微分光谱与土壤含水量的相关系数在部分波段比反射率光谱明显增加，部分波段相对减小，且一阶微分光谱与土壤含水量之间相关性时正时负，不呈单一相关（如图2所示）。相关系数的正峰值分别出现在波段1 004～1 013，1 174～1 183，1 194～1 203，1 467～1 476，1 497～1 506，1 787～1 796，1 797～1 806nm处，相关系数值分别为0.37，0.52，0.39，0.34，0.23，0.76，0.69；相关系数的负峰值分别出现在波段1 064～1 073，1 264～1 273，1 324～1 333，1 527～1 536，1 952～1 961，1 972～1 981nm处，相关系数值分别为－0.25，－0.48，－0.40，－0.25，－0.58，－0.47。选取各区域中具有代表性或相关系数绝对值最大的波段作为敏感波段，即1 004～1 013，1 174～1 183，1 324～1 333，1 787～1 796，1 952～1 961nm处。

2.3.3 土壤含水量与包络线去除光谱相关性分析从图3可以看出，包络线去除光谱与土壤含水量在1 560～1 805nm为负相关，其余波长范围内为正相关。相关系数的峰值分别出现在波段656～665，736～745，1 014～1 023，1 274～1 283，1 467～1 476，1 517～1 526，2 062～2 071，2 242～2 251nm处，相关系数值分别为0.66，0.67，0.66，0.62，0.26，0.26，0.54，0.49。选取各区域中具有代表性或相关系数绝对值最大的波段作为敏感波段，即656～665，1 014～1 023，1 517～1 526，2 062～2 071，2 242～2 251nm波段。

图2 土壤含水量与一阶微分光谱相关系数

图3 土壤含水量与包络线去除光谱相关系数

2.4 土壤含水量反演模型建立

选取反射率、一阶微分、包络线去除光谱与土壤含水量敏感的波段进行逐步回归分析。通过比较，建立土壤含水量反演多项式（表1—3），选择R2值最大的反演模型，作为土壤光谱与土壤含水量的关系模型。

表1 反射率光谱与土壤含水量逐步回归分析模型

表2 一阶微分光谱与土壤含水量逐步回归分析模型

表3 包络线去除光谱与土壤含水量逐步回归分析模型

由表1—3可知，预测土壤含水量的最佳模型是以波段X661，X1019和X2067的土壤包络线去除光谱为自变量建立的回归方程：

其预测R2值为0.85。采用实测含水量数据与式（12）得到的预测值进行拟合，检验结果的R2值为0.80（如图4所示）。

图4 土壤含水量预测与真实值比较

若不对原始影像进行分解，直接采用未剔除植被光谱的像元反射率、一阶微分、包络线去除光谱建立土壤含水量估测模型，得到的最佳模型是以波段X541，X979和X1632的一阶微分光谱为自变量建立的回归方程：

其预测R2值为0.36。采用实测含水量数据与式（13）得到的估测值进行拟合，检验结果的R2值仅为0.06。因此，剔除植被光谱干扰能够更有效地反演土壤含水量。

3 结论

（1）将混合像元分解问题归结为一个基于光谱匹配的非线性最优化问题，并针对最小二乘法的不足，引用了一种基于相关系数匹配的混合像元分解技术。高光谱数据本身提供了丰富的光谱信息，可以预见在高光谱数据中利用基于相关系数匹配的混合像元分解算法能够得到更为精确的结果。

（2）剔除植被光谱干扰后，土壤反射率、一阶微分、包络线去除光谱的部分波段能够敏感地反映土壤含水量的变化，最佳模型的估测R2值为0.85；若直接采用原始影像光谱建立估测模型，估测R2值仅为0.36。

（3）通过高光谱影像分解剔除植被光谱干扰估测土壤含水量的方法是可行的，可为今后遥感估测土壤含水量的研究提供参考。

［1］ 周萍，王润生，阎柏琨，等.高光谱遥感土壤有机质信息提取研究［J］.地理科学进展，2008，27（5）：27-34.

［2］ 陈荣毅，张元明，潘伯荣，等.古尔班通古特沙漠土壤养分空间分异与干扰的关系［J］.中国沙漠，2007，27（2）：257-265.

［3］ 宋韬，鲍一丹，何勇.利用光谱数据快速检测土壤含水量的方法研究［J］.光谱与光谱学分析，2009，29（3）：675-677.

［4］ Bedidi A，Cervelle B，Madeira J，et al.Moisture effects on visible spectral characteristics of lateritic soils［J］.Soil Science，1992（153）：129-141.

［5］ 孙建英，李民赞，唐宁，等.东北黑土的光谱特性及其与土壤参数的相关性分析［J］.光谱学与光谱分析，2007，27（8）：1502-1505.

［6］ 刘焕军，张柏，宋开山，等.黑土土壤水分光谱响应特征与模型［J］.中国科学院研究生院学报，2008，25（4）：503-509.

［7］ 王昌佐，王纪华，王锦地，等.裸土表层含水量高光谱遥感的最佳波段选择［J］.遥感信息，2003（4）：33-36.

［8］ 王静，何挺，李玉环.基于高光谱遥感技术的土地质量信息挖掘研究［J］.遥感学报，2005，9（4）：438-445.

［9］ 吴代晖，范闻捷，崔要奎，等.高光谱遥感监测土壤含水量研究进展［J］.光谱学与光谱分析，2010，30（11）：3067-3071.

［10］ 杨伟，陈晋，松下文经，等.基于相关系数匹配的混合像元分解算法［J］.遥感学报，2008，12（3）：454-461.

［11］ 丁丽霞，王志辉，葛宏立.基于包络线法的不同树种叶片高光谱特征分析［J］.浙江林学院学报，2010，27（6）：809-814.

［12］ 尹业彪，李霞，赵钊，等.沙质土壤含水率高光谱预测模型建立及分析［J］.遥感技术与应用，2011，26（3）：355-359.

［13］ 牛宝茹，刘俊蓉，王政伟.干旱半干旱地区植被覆盖度遥感信息提取研究［J］.武汉大学学报：信息科学版，2005，30（1）：27-30.