面向多光谱生态遥感观测目标的波段配组方法

(1.同济大学 测绘与地理信息学院，上海 200092；2.上海同繁勘测工程科技有限公司，上海 201900; 3.福建海洋研究所 福建省海岛与海岸带管理技术研究重点实验室，福建 厦门 361013)

自然生态资源的高精度定量监测与演化分析是国家自然资源统一管理的重要科学基础，可为区域生态红线划定和资源环境承载能力研究提供有效支撑[1]，而发展衔接天地观测的机载多光谱成像科学系统是有效解决途径[2]。目前，自然资源研究逐渐由宏观转向中微观，自然资源调查要素的数量也在不断攀升，针对特定目标生态参量的精确观测，发展可定制波段、可灵活组配波段的多光谱成像设备，是一种可行且有效的方案。多光谱生态遥感的目的是构建特征波段和生态参量之间的精确映射关系[3]，并基于此揭示遥感观测目标的生态指标与光谱之间的交互机理。选择针对观测目标的特定波段配组方案，是在进行遥感监测与反演前需要重点解决的问题。波段配组是指面向特定的多光谱遥感生态观测任务预先组配好适用的多光谱波段参量，建立不同观测目标和观测设备光学参数之间的特定关联。

生态遥感观测目标的光谱信息是由各个波长的连续数据组成，波段选择的目的是从高维的光谱数据中选择一小部分波段，在保留地物重要光谱信息的同时，消除光谱的冗余[4]。波段选择的常用方法有基于排序的方法[5]、基于降维的方法[6]、基于聚类的方法[7]以及基于机器学习的方法[8]等，这些方法主要是从高维的光谱数据中提取高信息量的波段，但只对光谱数据进行处理，并未考虑到与生态参量之间的关系。针对生态遥感，分析生态参量和光谱之间的相关性，提取高相关性的波段也是非常重要的波段选择方法[9]。因此，面向特定的生态观测目标，需要综合考虑高信息量波段及其与生态参量之间的关系提出一套混合的波段选择方法。

针对多光谱传感器的波段配组，除了考虑选择的波段，带宽的选择也是需要重点考虑的因素[10]。在高维的光谱数据中，大部分波段的数据是冗余的[11]，冗余波段的信息如果与被选择的波段信息混合会造成生态反演精度降低。目前针对水质遥感和植被遥感，带宽已被证明会影响反演精度，并且这些研究也提出了确定最大允许带宽的方法[12]。但对于多光谱传感器，最小带宽也需要计算以确定最佳的带宽范围，因此需要提出一套综合考虑最大最小带宽的方法。

图1 本方法总技术路线图Fig. 1 Flowchart of the method

为此，本研究提出一种面向多光谱生态遥感观测目标的波段配组方法，包括混合的波段选择方法以及最佳带宽范围的分析方法，为可定制波段的多光谱传感器提供波段配组依据。

1 多光谱波段配组方法

本方法的总技术路线如图1所示，共包含两大模块，分别是配组波段中心波长提取和带宽分析。配组波段中心波长提取主要包含光谱自身高信息量波段提取以及光谱与生态指标高相关性波段提取，最后通过交集分析的方式选择优先提取的波段。带宽分析主要包含最大带宽分析和最小带宽分析，最佳带宽为两者之间的阈值。

1.1 配组波段中心波长提取

提取光谱自身高信息量的波段可以从信息熵和降维两方面进行分析[13]。

信息熵可以分析出各个波段信息的含量，信息熵越大的波段所含信息量越大，因此选择信息熵相对较大的波段。信息熵的计算公式如下所示：

(1)

其中：H(band)代表某一波段的信息熵，p(Rrs'i)代表某一波段出现某一反射率Rrsi的概率。通过计算每个波段的信息熵可获得信息熵值的曲线，最终选取信息熵相对最大的前20%波段作为提取波段[14]。

降维的目的是提取高维数据中对信息贡献量最大的部分，因子分析(factor analysis, FA)是降维方法的一种，可以分析数据中所包含的公共因子，并且给出各个因子中对信息贡献最大的部分。公共因子数量的确定可以通过碎石分析的方式确定，碎石分析通过计算各波段间的相关系数获得相关系数矩阵，再计算相关系数矩阵的特征值，一般选取大于1的特征值数量作为公共因子的数量[15]。相关系数的计算公式如下：

(2)

其中：R是相关系数值，Xa和Xb分别是两个波段的反射率数据一维矩阵，Cov(Xa,Xb)是Xa和Xb的协方差，Var[Xa]是Xa的方差，Var[Xb]是Xb的方差。在确定了公共因子数后可以计算得到公共因子向量、特殊因子向量以及因子载荷矩阵，对载荷矩阵进行旋转，以对公共因子做更好的解释。从旋转后的因子载荷矩阵中可以获取各个波段对于各个公共因子的贡献率，对每个公共因子的贡献率进行[0，1]归一化，选取归一化后贡献率大于0.95的波段作为因子分析所提取的波段。

由于波段反射率比值是遥感生态参量反演回归分析中常用的变量[16]，因此本研究的光谱与生态指标高相关性波段提取主要采用波段反射率比值与目标生态参量相关性分析的方法，分析公式如下：

(3)

其中：Ratab是两个波段反射率比值所构成的一维矩阵，I是观测目标的某一种生态参量值所构成的一维矩阵。通过公式计算可以获得相关系数所组成的二维矩阵，同样将矩阵中的数值进行[0，1]归一化，取相关系数大于0.95组成的比值的两个波段作为相关性分析所提取的波段[17]。

光谱自身高信息量波段提取以及光谱与生态指标高相关性波段提取所获取的波段均为最终确定配组的备选波段，最后选择至少两种方法结果的交集部分作为优先选择的配组波段波长范围。若出现两种方法以上均存在交集的情况时，则交集部分的波段优先级比两种方法交集部分更高。

1.2 配组波段带宽分析

最大带宽可通过计算不同带宽的反射率与中心波长反射率的平均百分比误差(absolute percent difference,APD)来决定。从提取的配组波段波长范围中选择一个波长，将该波长的反射率设置为准确的反射率值，然后将带宽从0向左右两侧扩大，每次带宽扩展2 nm。通过遥感反射率计算公式可得到扩大带宽后的等效反射率，遥感反射率的一般公式如下：

(4)

其中：Rrs为遥感反射率，L为目标物反射的辐亮度，Ed为入射的辐照度。扩大带宽后的反射辐亮度和入射辐照度为扩大的带宽范围中所有波段反射辐亮度和入射辐照度的总和。

计算反射率准确值和等效反射率的绝对百分比误差，公式如下：

(5)

其中：APD为绝对百分比误差，n为测量样本反射率曲线的数量，Rrsacc为反射率准确值，Rrsnew为等效反射率值。最终可以获取不同带宽的绝对百分比误差曲线，百分比误差应在0.25%以内[12]，因此选取绝对百分比误差为0.25%的带宽值作为最大允许的带宽。

最小带宽则从多光谱传感器的辐射分辨率角度分析，需要保证传感器能够分辨出观测目标在被观测时可能反射的最小辐亮度，需要满足的条件公式如下：

(6)

其中：bc为某波段带宽的中心波长，bw为带宽，fLen(b)为传感器在各个波段的透射率，L(b)为目标物在各个波段反射的辐亮度，Radmin为多光谱传感器在以bc为中心波长的波段所能识别的最小辐亮度。

最终，配组波段的带宽范围应在允许的最大及最小带宽之间。若带宽允许范围>10 nm，多选择10的倍数作为带宽；若范围小于10 nm，则宜选择整数带宽。

2 面向遥感水质监测的波段配组案例

2.1 数据采集

本研究使用美国ASD公司产HandHeld2 Pro手持地物光谱仪和标准反射板白板进行水体的光谱数据采集，其中光谱仪波长范围为325～1 075 nm，采样带宽为1.4 nm，采集步骤严格遵循水面以上测量法[18]。在采集水体光谱样本的同时采集水质样本，于实验室进行分析，获取水质指标的浓度，本研究的水质指标浓度数据由污染控制与资源化研究国家重点实验室分析获得。

本研究的光谱数据与水质样本的采集地点位于上海市宝山区杨行镇，该区域为城市典型的建成区，主要为工业区和居住区，共计16条河道，其中有2条为运河，其余为景观河流。采集点位共114个，采集时间2018年8月4日—9月6日，每个点位采集5组光谱信息，取平均值作为最终结果，水质指标选择化学需氧量和氨氮，这两个指标是日常城市水质监测中的重要指标，也是国家标准《地表水环境质量标准(GB 3838—2002)》[19]中重要的水环境指标。采集点位的分布如图2所示。

图2 光谱数据与水质样本数据采集点位Fig. 2 Sampling points of spectral data and water quality samples

图3 光谱数据信息熵分析结果Fig. 3 Information entropy analysis result of spectral data

2.2 波段配组结果

2.2.1 配组波段中心波长提取结果

信息熵分析结果如图3所示，信息熵的最大值为3.51，波长为591 nm，信息熵值相对最大的前20%的波长范围为558～658 nm，集中在绿色、黄色、红色波段，为该范围内的波段及为信息熵分析所提取的波段范围。

在使用因子分析波段提取方法前，需要采用碎石分析得到公共因子的数量，碎石分析的结果如图4所示，其中前4个特征值>1，分别为252.62、12.16、4.19和1.31，因此公共因子的数量为4。因子分析的结果如图5所示，得到4个含有最大信息量的成分，归一化后成分1贡献率>0.95的范围为560～680 nm，主要为黄绿波段；成分2为825～892 nm，主要为近红外波段；成分3为400～410 nm，主要为紫色波段；成分4为703～717 nm，主要为红色波段。

图4 碎石图Fig. 4 Scree plot

反射率比值与水质指标的相关性分析结果如图6所示，图中的相关系数值均已归一化。图6(a)为化学需氧量浓度与反射率比值的相关系数图，其中相关系数>0.95的范围是600～690 nm与703～715 nm反射率的比值部分，为黄色、红色波段与红色波段的比值。图6(b)为氨氮浓度与反射率比值的相关系数图，其中相关系数>0.95的范围是480～520 nm与662～708 nm反射率比值部分，主要为蓝绿波段与红色波段的比值。

通过以上3种方法的分析，分别获得4组以化学需氧量和氨氮为观测目标的配组波段中心波长的提取结果，这些结果均为备选波段。图7将3种方法得到的备选波段范围标出，最终遵循优先提取波段的方法，选择至少3组结果交集的部分，分别为600～658 nm、662～680 nm、703～708 nm，这3个波长的范围即为配组波段中心波长提取的最终结果。

图5 因子分析波段提取结果Fig. 5 Band extraction result of factor analysis

图6 反射率比值与水质指标浓度相关性分析结果Fig. 6 Correlation analysis result between reflectance ratio and water quality index concentration

2.2.2 配组波段带宽分析结果

依据波段中心波长提取的波长范围，分别选择640、675、705 nm作为波段的中心波长进行最大及最小带宽分析。

图7 配组波段中心波长范围提取结果Fig. 7 Central wavelength ranges extraction result of band group matching

从图8可以看出，随着带宽的不断扩大，等效反射率与反射率准确值的绝对百分比误差呈不断上升的趋势，说明带宽会影响遥感对地物反射率的计算精度，带宽越大，与中心波长准确的反射率值偏差越大。选取绝对百分比误差为0.25%的带宽值作为最大允许带宽，因此以640 nm为中心波长，最大允许带宽为68 nm；以675 nm为中心波长，最大允许带宽为42 nm；以705 nm为中心波长，最大允许带宽为34 nm。

图8 反射率准确值与等效反射率的绝对百分比误差曲线Fig. 8 APD curve between accurate reflectance and equivalent reflectance

图9 数据集中相对最低的水体辐亮度曲线Fig. 9 Relative lowest water radiance curve in dataset

在114组光谱数据中心选取相对最低的水体辐亮度曲线，如图9所示，其中640 nm波段的水体辐亮度为0.000 96 W/(m2×nm×sr)，675 nm波段为0.000 83 W/(m2×nm×sr)，705 nm波段为0.000 97 W/(m2×nm×sr)。以同济大学测绘与地理信息学院自主研发的多光谱相机“极视1号”为例，该相机最低可以分辨的辐亮度值为0.1 μW/(cm2×nm×sr)。依据公式(6)，以传感器透射率的理想状态即100%来计算，由于带宽一般为整数，因此可以得到3个配组波段的带宽均需要大于2 nm。

表1为配组带宽分析的最终结果，其中640 nm波段的带宽范围应在2～68 nm之间，675 nm应在2～42 nm之间，705 nm应在2～34 nm之间。

2.3 波段配组结果验证

为验证本研究配组波段结果在实际水质遥感中的可行性，基于114组水质及光谱样本数据，进行水质反演的验证实验。选择640、675、705 nm为中心波长，10 nm为带宽，以3个波段的等效反射率为输入数据。此外，选取450、750、800 nm作为对比波段，这3个波段不在任何一种方法提取的特征波段范围内，同样以10为带宽，以3个波段的等效反射率为输入数据。 依据国家标准《地表水环境质量标准(GB 3838—2002)》[19]，将水质指标浓度分为6个级别，分别为Ⅰ类、Ⅱ类、Ⅲ类、Ⅳ类、Ⅴ类以及劣Ⅴ类，水质分级的标准如表2所示，本研究样本数据集中各个类别的数量如表3所示，反演模型的输出数据即为水质指标的水质等级。水质反演模型选择支持向量机(support vector machine, SVM)，该方法已实际应用于遥感水质反演中[20]，且对于小样本数量有较好的分类精度。

表1 配组波段带宽分析结果Tab. 1 Bandwidth analysis result of band group matching nm

将样本数据集随机分为训练集和测试集，其中训练集数据量为91组，测试集为23组，支持向量机水质反演模型采用多分类方法，核函数选择径向基核函数。模型训练的结果精度检验采用多分类常用的指标宏平均，主要包括宏平均的准确率、召回率及F1分数，计算公式如式(7)～(12)，此外，通过混淆矩阵验证各个类别自身的准确分类精度。

表2 水质分级标准Tab. 2 Water quality grading standards mg/L

表3 数据集各等级样本数量Tab. 3 Sample number of each grade of the dataset

表4 反演模型宏平均准确率、召回率和F1分数Tab. 4 Macroprecision, macrorecall and macro F1-score of the inversion models

(7)

(8)

(9)

(10)

(11)

(12)

其中：TP为将正类预测为正类的数量，FP为将负类预测为正类的数量，FN为将正类预测为负类的数量。公式(7)～(9)分别为准确率、召回率和F1分数的计算公式，公式(10)～(12)分别为宏平均准确率、召回率和F1分数的计算公式。

表4为化学需氧量和氨氮反演模型的精度验证结果，采用本研究配组波段构建的化学需氧量反演模型的宏平均准确率、召回率和F1分数分别为0.77、0.74、0.73，氨氮反演模型的宏平均准确率、召回率和F1分数分别为0.81、0.88、0.83。采用对比波段构建的化学需氧量反演模型的宏平均准确率、召回率和F1分数分别为0.51、0.55、0.52，氨氮反演模型的宏平均准确率、召回率和F1分数分别为0.73、0.58、0.57。采用本研究配组波段模型的混淆矩阵如图10所示，化学需氧量反演模型的Ⅱ类、Ⅲ类、Ⅳ类分级精度较好，均达到了0.8以上，Ⅴ类以及劣Ⅴ类由于样本数量较少分级精度为0.5；氨氮反演模型各个类别分级精度均较好，Ⅱ类、Ⅳ类的分级精度在0.6以上，其余分级精度均为1。采用对比波段模型的混淆矩阵如图11所示，可以看到化学需氧量和氨氮的分级精度都相对较差。因此本研究针对化学需氧量和氨氮的波段配组方案实际可行。

图10 本研究配组波段构建的反演模型混淆矩阵Fig. 10 Confusion matrices of the inversion models built with the selected bands

图11 对比波段构建的反演模型混淆矩阵Fig. 11 Confusion matrices of the inversion models built with the bands for comparison

3 结论

为揭示遥感观测目标的生态指标与光谱之间的交互机理，建立不同观测目标和观测设备光学参数之间的特定关联，需要针对特定观测目标选择特定的波段配组方案。本研究提出一种面向多光谱生态遥感观测目标的波段配组方法，该方法包括波段提取和带宽分析。通过综合考虑高信息量波段和与生态参量之间的关系，提出一种混合的波段提取方法。为确定各波段的最佳带宽范围，提出一种综合考虑最大最小带宽的带宽分析方法。

以水质遥感为例，通过信息熵分析、因子分析、相关性分析，提取针对化学需氧量和氨氮的配组波段，基于平均百分比误差和仪器辐射分辨率的约束，确定了配组波段的带宽允许范围。以该配组方案的波段为中心波长，选取带宽允许范围中的某一值为带宽，计算等效反射率并结合水质指标等级作为水质反演模型的训练数据，通过模型的精度验证得到化学需氧量的总体分级精度在0.7以上，氨氮的总体分级精度在0.8以上，证明该波段配组方案实际可行。

本研究的波段配组方法可以为特定目标生态参量的精确观测提供理论依据，支撑自然资源生态要素的定量监测与演化分析，并为可定制波段的多光谱传感器提供波段配组依据。