热红外航空成像光谱仪(thermal airborne hyperspectral imager, TASI)是加拿大研制的先进的机载热红外高光谱设备。 该设备在8~11.5 μm范围内有32波段, 波段间隔为0.109 5 μm, 半高宽为0.054 8 μm, 总视场角为40°(表1)。 本文利用核工业北京地质研究院遥感信息与图像分析技术国家级重点实验室的TASI航空高光谱成像系统, 在东北黑龙江省海伦地区附近获取了高空间分辨率的高光谱热红外遥感数据, 在经过大气校正等数据预处理后, 采用杨杭等^[15]改进的TES(temperature-emissivity separation)分离方法进行温度与发射率的分离, 获取研究区发射率数据。

表1

Table 1

表1(Table 1)

表1 TASI各通道中心波长 Table 1 TASI channel center wavelength 波段号 波长/μm 波段号 波长/μm 波段号 波长/μm 波段号 波长/μm 1 8.054 8 9 8.930 8 17 9.806 8 25 10.682 8 2 8.164 3 10 9.040 3 18 9.916 3 26 10.792 3 3 8.273 8 11 9.149 8 19 10.025 8 27 10.901 8 4 8.383 3 12 9.259 3 20 10.135 3 28 11.011 3 5 8.492 8 13 9.368 8 21 10.244 8 29 11.120 8 6 8.602 3 14 9.478 3 22 10.354 3 30 11.230 3 7 8.711 8 15 9.587 8 23 10.463 8 31 11.339 8 8 8.821 3 16 9.697 3 24 10.573 3 32 11.449 3

表1 TASI各通道中心波长 Table 1 TASI channel center wavelength

土壤样品采集地点为东北黑龙江省海伦地区, 该地区土壤类型主要为黑土, 土壤腐殖质层较厚, 约30~60 cm, 有机质含量在2.5%~4.5%之间, 粘粒含量在40%~60%之间, 属于粘土, 西南部旱田土壤基本为中性, 中部及东北部旱田土壤为酸性。 在工作区不同位置共采集土壤样本40个(图1)。 测区表层为黑色腐殖质层, 当天同步飞行采集表层0~20 cm的土样, 剔除大的植物残茬、 石砺等杂物, 置于实验室风干研磨, 过0.15 mm筛选用于含量测定。 全钾含量采用X射线荧光光谱法测定(表2)。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 TASI数据采集区域及采样点分布图Fig.1 TASI data area and sampling point distribution map

表2

Table 2

表2(Table 2)

表2 土壤样本钾(K)元素含量信息表 Table 2 The element content information table of soil samples 样本 类别 最小值/ (g·kg-1) 最大值/ (g·kg-1) 平均值/ (g·kg-1) K 建模 2.39 2.57 2.48 预测 2.45 2.62 2.51

表2 土壤样本钾(K)元素含量信息表 Table 2 The element content information table of soil samples

探究多个自变量与因变量之间关系的方法很多, 对于光谱数据, 常用的建模方法有最小二乘法、 多元逐步回归拟合、 BP神经网络、 支持向量机SVM等。 机器学习类方法虽然能较好地训练模型并进行预测, 但是存在较多的人为干预调参问题, 同时对于某一元素在研究区整体的提取具有较大难度。 本研究主要聚焦于偏最小二乘法与多元逐步回归拟合, 尤其利用多元逐步回归方法研究发射率与K元素含量关系时创新性地使用了全二次逐步回归进一步提高模型的精度。

相比常规的n元线性逐步回归仅有的常数项和线性项, 多元全二次逐步回归引入了交叉乘积项和平方项进行回归, 以对回归方程中常数项、 线性项和二次项进行的t检验的p值是否小于等于0.05为判定依据, 依次引进显著项剔除非显著项, 可以有效地解决自变量的多重共线性问题。 但参数不宜引入过多以免数据产生过拟合, 其中, 引入了4个参数进行全二次多元逐步回归方程如下所示, 由于增加了更多的系数, 因此能够更加精确地进行回归模型的建立。 本次研究所有建模方法均利用Matlab编程实现。 TASI影像数据预处理及温度与发射率分离利用ENVI5.3-IDL编程实现。

yi=b0+∑i=14bixi+∑i=13∑j=i+14bijxixj+∑i=14biixi2εi~N(0, δ2), i=1, 2, …, n

式中, y_i为响应变量(预测值); b₀, b_i, b_ij和b_ii分别为回归方程常数项、 线性项、 交叉乘积项、 平方项的系数; x_i为预测变量(输入值), 本文中为输入的四个相关性强波段值; ε _i~N表示数据服从期望值μ为0的正态分布。

物体的发射率除了取决于其材质, 更取决于其所存在的环境。 对于土壤, 其热红外发射率光谱曲线形态主要取决于土壤中所含的矿物种类以及含量、 水分、 有机质含量和温度等因素。 研究所采集的40个土壤样本的发射率光谱曲线在8~9.6 μm变化趋势基本一致(图2), 9.6~11.45 μm后形态有所变化但基本变化趋于平缓, 发射率值在8.38, 8.6以及9.26 μm处出现了三处非常明显的波谷, 其中8.38与9.3 μm的两处波谷出现了类似石英的波谱特征。 夏军等研究发现, 土壤中硅酸盐矿物导致发射率光谱曲线呈现明显的Reststrahlen吸收特征, 即不对称双吸收谷, 两个吸收谷分别位于8.23和9.27 μm波长附近, 且后一个吸收谷较深, 宽度较大^[13]。 9.6 μm后光谱曲线出现不同波动是由于土壤中所含不同矿物的基团内部振动产生的谱带不同所致, 但总体幅度变化较小。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 40个土壤样本的热红外发射率光谱曲线Fig.2 Thermal infrared emissivity spectral curves of 40 soil samples

通过40个土壤样本K元素含量与发射率各波段做皮尔森相关性分析可以看出, K元素与32个波段发射率呈负相关关系, 其中相关系数大于0.6呈强相关的波段有5个, 分别是6, 11, 15, 22和23波段; 相关系数介于0.4~0.6中等相关的波段有17个(图3), 整体具有较强的相关性。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 钾元素含量与发射率相关性Fig.3 Correlation coefficient curve of potassium content and its thermal infrared emissivity

由于多元全二次逐步回归引入的参数数量直接决定了模型建立引入系数的多少和复杂程度, 因此, 为了避免建模时系数过多导致过拟合增大无谓的计算量, 初步选取与K元素含量相关系数最高的4个波段用于模型的建立。 土壤中钾元素有多种赋存形态, 大部分以原生或次生的结晶硅酸盐状态存在于土壤中, 其中云母族矿物参考(白云母、 黑云母)及富钾长石(正长石)中含钾元素最多, 白云母、 黑云母以及正长石特征吸收位置选择6, 11, 15和23波段作为特征波段(表3)。

表3

Table 3

表3(Table 3)

表3 土壤元素含量与发射率显著相关所对应的TASI波段 Table 3 The TASI band corresponding to the significant correlation between soil element content and emissivity 元素名称 波段号 波长/μm 相关系数 对应方程系数 6 8.602 -0.658 x1 K 11 9.150 -0.673 x2 15 9.588 -0.645 x3 23 10.464 -0.627 x4

表3 土壤元素含量与发射率显著相关所对应的TASI波段 Table 3 The TASI band corresponding to the significant correlation between soil element content and emissivity

将40个土壤样本随机分为两组, 其中32个样本用于含量预测模型的建立, 8个样本用来测试模型的精度。 验证K元素含量数据符合正态分布后, 以所选4个特征波段发射率数据作为自变量, K元素含量为因变量, 以对回归方程中常数项、 线性项和二次项进行的t检验的p值是否小于等于显著性水平0.05为判定依据, 依次引进显著项剔除非显著项, 同时对模型总体进行F检验的p值是否小于等于显著性水平0.05来验证样本观测值与总体假设值是否存在显著性差异从而建立模型[图4(a, b)]。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 逐步回归模型预测值与真实值拟合效果图 (a): 常规多元逐步回归; (b): 全二次多元逐步回归Fig.4 Stepwise regression model prediction value and true value fitting effect plot (a): Constant multiple stepwise regression; (b): Quadratic multiple stepwise regression

常规逐步回归建立的回归拟合模型均方根误差RMSE为0.031, 调整后的判定系数R²为0.569, 测试集的均方根误差RMSE为0.031, 调整后的判定系数R²为0.78[见图5(a)]; 全二次多元逐步回归建立的回归拟合模型均方根误差RMSE为0.027, 调整后的判定系数R²为0.667, 测试集的均方根误差RMSE为0.032, 调整后的判定系数R²为0.82[见图5(b)], 所有指标均通过p值小于0.05的显著性验证(表4)。 通过模型拟合结果及评价指标来看, 全二次多元逐步回归比常规多元逐步回归建模精度以及验证精度均有所提高。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 钾元素含量实测值与预测值散点图 (a): 常规多元逐步回归; (b): 全二次多元逐步回归Fig.5 Scatter plot of measured value and predicted value of potassium content (a): Constant multiple stepwise regression; (b): Quadratic multiple stepwise regression

表4

Table 4

表4(Table 4)

表4 逐步回归模型建模结果 Table 4 Stepwise regression model modeling results 建模方法 变量数 建模 测试 Adjusted R2 RMSE/% 模型F检验的p值 Adjusted R2 RMSE/% 模型F检验的p值 常规多元逐步回归 4(参数) 0.569 3.1 6.25×10-6 0.78 3.1 0.002 25 全二次多元逐步回归 7(参数) 0.667 2.7 9.25×10-6 0.82 3.2 0.001 22

表4 逐步回归模型建模结果 Table 4 Stepwise regression model modeling results

为进一步提高建模精度, 利用Matlab回归诊断的学生化残差来进行模型改进。 通过|Se_i|> 2来查找远离数据集中心观测点即异常点, 剔除异常点来进一步提高模型精度。 同时, 以32个波段发射率数据为自变量, K元素含量为因变量进行偏最小二乘法建模, 进一步对比三种模型的优劣。

Sei=ei/MSE(1-hii)(2)

式(2)中, Se_i即学生化残差; e_i=y_i- y^i为第i个观测对应的残差; MSE=SSE/(n-1-p)为均方残差; h_ii为帽子矩阵H=X(X^TX)^-1X^T对角线上的第i个元素。

偏最小二乘法建立的回归拟合模型入选主成分数为2, 均方根误差RMSE为0.033, 判定系数R²为0.45, 测试集的均方根误差RMSE为0.037, 判定系数R²为0.51(图6)。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 钾元素含量实测值与预测值散点图Fig.6 Scatter plot of measured value and predicted value of potassium content

通过评价指标分析改进后的多元逐步回归模型发现, 虽然建模精度有所提高, 但测试集的精度却均有所下降。 与改进前模型相比, 常规多元逐步回归建模样本的均方根误差RMSE降低了0.7%, 判定系数R²提高了0.163; 测试集的均方根误差RMSE提高了0.2%, 判定系数R²降低了0.015; 全二次多元逐步回归均方根误差RMSE降低了0.71%, 判定系数R²提高了0.135; 测试集的均方根误差RMSE提高了0.2%, 判定系数R²降低了0.1, 同时由于剔除了某些不显著的变量, 模型再次引入了新的变量, 参数从7个增加到了10个(表5)。 改进后的模型训练集精度上升而测试集精度下降的原因应该是训练集数据发生了过拟合, 从实验结果分析, 改进前的模型泛化能力更强, 更适用于研究区K元素的反演。 同时对比偏最小二乘法建模, 全二次多元逐步回归各项评价指标均优于其余两种方法。 由于混合像元影响以及样本的选择有所差异, 模型整体拟合精度不是很高, 但本研究提出新的逐步回归方法有效地提高了模型的精度。

表5

Table 5

表5(Table 5)

表5 建模结果对比 Table 5 Comparison of different models 建模方法 变量数 建模 测试 Adjusted R2 RMSE/% 模型F检验的p值 Adjusted R2 RMSE/% 模型F检验的p值 常规多元逐步回归 4参数 0.569 3.1 6.25×10-6 0.78 3.1 0.002 25 全二次多元逐步回归 7参数 0.667 2.7 9.25×10-6 0.82 3.2 0.001 22 学生化残差去除异常值后 建模方法 变量数 建模 测试 Adjusted R2 RMSE/% 模型F检验的p值 Adjusted R2 RMSE/% 模型F检验的p值 常规多元逐步回归 4参数 0.732 2.4 6.37×10-8 0.765 3.3 0.002 78 全二次多元逐步回归 10参数 0.802 1.99 7.6×10-7 0.72 3.4 0.004 97 偏最小二乘 2成分 0.45 3.3 0.51 3.7

表5 建模结果对比 Table 5 Comparison of different models

相比常规多元逐步回归仅考虑常数项和线性项, 全二次多元逐步回归能够引入更多的变量参与到回归模型的建立中, 从而提高模型的反演精度。 说明利用TASI数据的相关波段通过全二次多元逐步回归方法反演元素含量是可行的, 比起传统的化学填图, 遥感反演的方法在损失部分精度的条件下能够高效大面积地反演某个地区元素含量。