基于地面光谱水稻重金属胁迫监测光谱特征尺度识别

黄芝,刘湘南,赵爽,张仙

(1.衡阳师范学院城市与旅游学院， 湖南 衡阳 421002；2.中国地质大学(北京)信息工程学院， 北京 100083；3.天津城建大学地质与测绘学院， 天津 300384；4.中国自然资源航空物探遥感中心， 北京 100083)

工业化和城镇化进程加速，但环保措施相对滞后，我国农田重金属污染呈现越演越烈的趋势，研究表明，我国约有16.67%耕地受到重金属污染[1]。重金属具有隐蔽性强、不能被生物降解等特点，可由食物链进入人体，并在人体内富集，进而导致食物中毒[2-3]。解决重金属污染，尤其是大面积农田土壤重金属污染，从而保证农产品质量安全和农田生态系统健康已成为目前面临的严峻且棘手的问题[4]。治理和防治的前提在于有效地监测农田重金属污染，掌握受污染作物的响应特征以及有效地捕获其特征。高光谱遥感所提供的精细光谱可有效地分析作物叶片的生化成分，为食品安全与污染治理提供数据支撑且监测效果显著[5-6]。

已有研究集中探索了重金属胁迫下作物的生理生态参数变化及其相应的光谱响应特征[7-8]，分析了利用高光谱遥感手段进行农作物重金属胁迫监测的可行性和光谱弱信息提取方法的适宜性[9-10]，但都基于特定的尺度，对于遥感监测中所采用尺度的合理性和适宜性缺乏深入探讨。作物光谱信息的采集通常有两种尺度：一是田间尺度，一般基于地面实测的高光谱数据获取；二是较大面积的区域尺度，基于卫星传感器的多/高光谱数据获取。不同尺度获取的同一重金属污染胁迫水平的高光谱信息有何差异？利用高光谱遥感数据能明显监测到重金属污染胁迫，体现不同胁迫水平差异需要多大的尺度？这些问题都有待深入研究，因此，为保证作物重金属胁迫监测的有效性，不仅需充分探究作物受污染后的光谱特征响应，还需研究所选光谱数据的适宜性。

特征尺度是较复杂的概念，没有明确的定义，有学者指出：“特征尺度定义了地表过程的空间间隔以及明确了在给定应用中利用遥感数据进行监测的空间范围”[11]，即特征尺度是相对于特定应用目标而言的，取决于目标的自身特性。农田重金属胁迫中涉及的光谱特征尺度问题包括两个方面：首要是正确定义观测对象特征，即寻找重金属胁迫导致的作物光谱响应特征，已有的研究多是基于此进行[12-13]。其次是范围的框定，具体在于利用获取到的光谱数据找到基于作物光谱响应特征下的明显监测到重金属胁迫的适宜尺度。与“空间尺度”不同，“光谱尺度”在遥感领域的提法相对较少，但通过不同光谱分辨率的遥感数据所得到的各类信息存在一定的差异[14]，因而尺度效应的概念也可广义上推及到光谱方面[15-16]。在光谱尺度选择问题上，已有研究大多集中探讨各种光谱特征提取技术[17-18]，对光谱尺度适宜性研究相对匮乏。

鉴于此，本文选取受重金属胁迫较严重的湘江流域附近的农田为研究区，基于所测两种不同污染胁迫水平下的地面高光谱数据，通过分析它们各自的光谱响应特点，利用小波分析方法对高光谱信息进行分解，提取小波多尺度特征参数和分形维数，基于拐点思想探讨水稻重金属胁迫的有效光谱特征尺度，为进一步采用合适光谱分辨率的高光谱卫星影像进行水稻重金属胁迫监测提供依据，以期提高大范围农田重金属污染监测的准确性，为综合治理我国土壤污染提供科学的支撑。

1 材料与方法

1.1 研究区概况

研究区位于湖南省株洲市辖区内，地处湖南省东部，介于26°03′05″N—28°01′07″N和112°57′30″E—114°07′15″E之间。该区土壤类型以红壤为主，有机质含量2%～3%，土壤PH介于5.0～8.5之间[19]，适宜多种农作物生长，是湖南省著名的粮食高产区，其粮食品种以水稻为主，田间管理以农民为主，当地政府提供技术支持为辅。该区水稻通常5月下旬移栽，9月下旬收割，整个生长期约为120 d左右。由于多年来环保手段相对不足，流经株洲市的湘江及其支流长期以来受到工业排放物的污染，由于直接引江水灌溉导致该地区农田重金属含量严重超标。

本文选取不同重金属胁迫等级的两个水稻区作为研究区，研究区A1靠近湘江，直接引流灌溉，重金属污染较重；而研究区A2污染相对较轻。研究区内的土壤和水稻重金属含量通过实地采样获取，根据国家标准和当地土壤背景值将研究区A1和A2的重金属胁迫等级分别定为“重度胁迫”和“轻度胁迫”(表1)。由于研究区内地势起伏不大，气候条件趋于一致，且田间管理手段都大致相同，因此假设整个生长期内研究区A1和A2的水稻都有充足的水分和养料供其生长，从而认为研究区内的水稻不受除重金属胁迫外其他环境因素的影响。

表1 研究区重金属含量表Table 1 Rice heavy metal concentrations of the two study areas

1.2 数据获取及处理

在研究区A1和A2分别选择40个采样点测量光谱数据(样本均匀分布)，采集日期为2014年6月15日、7月18日、8月28日，测量时间在北京时间10：00—12：00之间。水稻冠层及土壤光谱测量选择在风力较小，天气晴朗条件下进行。测量仪器为FieldSpec Pro光谱仪(美国ASD公司生产)，其波段范围：350～2 500 nm(其中350～1 000 nm的光谱采样间隔为1.4 nm，1 000～2 500 nm的光谱采样间隔为2 nm)。在进行光谱测量前首先利用标准白板进行光谱纠正并设定视场角为10°，此后每隔2°三个采样点进行一次白板纠正，光谱测量时传感器探头距水稻冠层顶部大约80～100 cm，并尽可能保持与地面垂直，最后取所有采样点的平均值表示这个时期光谱测量值。

1.3 研究方法

1.3.1小波离散变换 利用MATLAB小波工具箱的DB5小波函数对不同胁迫水平下的水稻原始反射光谱(450～900 nm)进行一维离散小波分解，设置小波分解尺度为8层。离散小波变换(discrete wavelet transformation，DWT)是在一系列离散尺度上对信号进行分析的方法。离散尺度一般为二阶(2,4,6,8……)，利用多种快速计算方法可以实现小波的变换。小波分解每一层的结果都是针对前一次尺度分解的信号再进行分解，将上一层所得的低频信号再次分解成一个低频信号和一个高频信号。在小波分解后可以得到一个小波近似系数和与相应尺度数目相等的小波细节系数，小波近似系数表征了信号分解重构的低频部分信息，是信号与尺度的内积，高频部分信息由小波细节系数表示，是信号与小波函数的内积。将A1、A2区域每个波长的40个采样点取平均值作为变量导入matlab R2011b软件中，使用wavedec函数对小波进行分解，将一个信号分解成指定数量n层然后返回每个层的小波系数。

[C,L]=wavedec(X,N,‘wname’)

(1)

式中，wavedec为使用的DB5小波函数，X为输入的光谱反射信号，N为小波分解的层数。小波“wname”将信号X分解成N层，C是一个列向量，存储经不同尺度分解所得的每层相应高频小波细节系数和最后一次分解所得到低频的小波近似系数，L也是列向量，存储了C中不同尺度分解的小波近似长度和相应各阶细节系数的长度。

1.3.2小波信息熵 小波信息熵可以度量多分解尺度下各子空间包含的信息量，水稻重金属胁迫遥感监测特征尺度是指在一定谱段范围内能区分胁迫水平的最佳光谱率大小，不同分解尺度下对应的信息量差异可作为胁迫水平的区分指标，因此，本研究将小波信息熵作为光谱特征尺度的识别参数，其可根据各尺度小波变换系数的概率分布计算。

(2)

对于任意尺度a，可先计算小波变换系数的统计直方图，再根据直方图确定该尺度的小波系数近似概率Pi(i=1，2，…M)，M为小波系数的区间数目。

1.3.3小波分形维数 为准确且全面识别不同胁迫水平的光谱特征尺度，本文采用盒维数法计算8层小波分解尺度下的小波分形维数，通过对比不同尺度下的分形维数变化点来判别区分不同胁迫水平的最佳尺度。盒维数法的计算过程为：用边长为ε的正方形格网去覆盖小波分解尺度j(j=1，2，…，J，J为最大分解尺度)的光谱小波近似系数重构曲线，覆盖有限体的网格数据记为N。边长ε与网格数存在以下形式的幂指数关系。

N(ε)=ε-D

(3)

当正方形网格边长为ε1,ε2,ε3,…,εk，时，覆盖有被测对象的正方形网格数据相应为N(ε1),N(ε2),N(ε3),…,N(εk)，两边同时取对数可得式(4)。

Log2N(ε)=-Dlog2ε+A

(4)

式中，A为待定常数，D为被测曲线的分维，其值等于公式(3)斜率的绝对值。

1.3.4特征光谱尺度验证分析 为验证小波信息熵和小波分形维数所识别出的光谱特征尺度的可行性，本文利用实测的ASD光谱特征波段构建对叶绿素含量存在敏感性的综合型叶绿素光谱指数MCARI/OSAVI(modified chlorophyll absorption ratio index)[20]、NDSI_R[21]、Depth[22]来探讨不同尺度下的光谱指数的敏感程度。

MCARI/OSAVI=

(5)

NDSI_R=(R598-R508)/(R598+R508)

(6)

Depth=

[(R700-R670)-0.2(R700-R550)]R700/R670)

(7)

利用高斯响应函数对水稻反射率光谱进行重采样，获取不同光谱尺度下的水稻重金属胁迫光谱敏感指数数据集，探究重金属胁迫高光谱指数随光谱分辨率的变化情况。具体而言，对所采集的A1、A2区域40个水稻样本的原始光谱数据进行高斯重采样，采样间隔一次为2、4、8、16、32、64、128和256 nm。高斯函数表示如下。

(8)

式中，μ为期望值，代表每个函数的中心波长；σ为标准差，用来计算每一个高斯函数的两端到中心的距离t。

高斯函数的个数为采样后的波段数，采样后的光谱值为每个高斯函数在区间[μ-t，μ+t]内的积分。

2 结果与分析

2.1 重金属胁迫敏感光谱波段分析

重金属污染胁迫遥感监测特征光谱尺度识别的首要前提在于正确定义观测对象的特征，即寻找适宜指示胁迫的敏感波段，有利于提高数据的利用效率和监测的准确度。水稻受重金属胁迫后，体内生理生化参数会发生改变，如叶绿素含量、叶面积指数等，且光谱反射特性也将发生相应变化。有研究表明，植株受叶片色素含量影响的高光谱波段为400～900 nm，此波段范围内水稻叶片的反射光谱变化最为敏感[13]。不同污染水平的ASD光谱波段中，波长1 250～1 500、1 750～2 000 nm处存在异常噪声点，数据质量较差，受水分和大气影响严重[23-24]，故剔除1 300～2 500 nm区间的光谱曲线，同时，1 000～1 200 nm波段存在断开现象(图1)，因此，本研究将450～900 nm波段作为指示水稻重金属胁迫的敏感特征。

图1 A1、A2区域原始水稻高光谱反射率Fig.1 Primitive hyperspectral reflectance of rice in A1 and A2 regions

2.2 小波分解细节系数分析

图2为不同污染程度下所选择的两个样本点进行DB5小波分解示例，研究区A1、A2的光谱信号经过8层小波分解，1～8层对应的分解尺度分别为2、4、8、16、32、64、128和256 nm。每层光谱信号经离散小波分解后，在第1～8层能够产生8个细节系数向量，最后八个小波细节系数被记录为d1、d2、d3、d4、d5、d6、d7和d8。由图3可知，随着小波分解层数量逐渐增加，小波细节系数减少,对应的小波细节系数光谱曲线的峰谷特征减弱。小波细节系数分解曲线在相同胁迫水平下呈现相似特征,不同胁迫水平差异明显。总结来说，小波分解可以反映出分解的多分辨率特性。小波系数能够反映波形的整体信息，在d1～d4小尺度分解下的高频信号可以明显看到光谱波段噪声的存在(图2 B、C、D、E)，在d5～d8大尺度分解下的低频信号噪声减少(图2 F、G、H、I)，即随着分解层数的增加，光谱特征和低频信号能够更多地显现出来。

注：A为A1和A2区域两个样本点对应的光谱敏感曲线；B～I分别为DB5小波分解的d1～d8细节系数。Note: A shows the original sensitive rice spectral curve in A1 and A2，and B～I show the detail coefficient of d1～d8 of DB5 wavelet decomposition, respectively.图2 研究区域分解后的小波细节系数和原始水稻光谱敏感曲线Fig.2 Wavelet detail coefficients of decomposition in research area and the original sensitive rice spectral curve

2.3 不同光谱尺度小波信息熵分析

小波信息熵能描述小波各分解尺度下不同胁迫水平光谱曲线携带的信息量差异，熵大所携带的信息量大；熵小所携带的信息量小，特征尺度对应于信息量较大的尺度[25]。图3为A1、A2区域40个实测点的小波信息熵的平均值柱状图，可以看出，不同胁迫水平下小波信息熵随尺度的变化走势大体相同，A1区(重度胁迫)的小波信息熵随尺度增加总体呈上升趋势，且在1～6尺度下上升速度较快，A1区小波信息熵的极大值集中出现在6、7尺度，A2区的小波信息熵值虽也呈上升趋势，但速度相对较缓。A1、A2 区小波信息熵的值有较大的差异，分解的第1尺度为不同胁迫水平差异的最大值，随着小波分解尺度增大，小波信息熵的差异值在1～5尺度内都逐渐减小，第6尺度A1区的小波信息熵值大于A2区，在第7尺度附近A1、A2区小波信息熵值持平，第6尺度为小波信息熵值差异变化的拐点，也可看作不同胁迫状态区分的拐点，即在该点之前不同重金属胁迫水平的小波信息熵差异较大，而该点之后变小。

图3 各尺度下不同污染水平小波信息熵Fig.3 Histograms of wavelet information entropy of different pollution levels

2.4 不同水稻胁迫水平多尺度分维数分析

为挖掘小波分形维数对水稻重金属高光谱弱信息提取潜力，分析小波分解不同尺度的分维数差异从而识别区分水稻重金属胁迫不同水平的特征光谱尺度，本文分别计算野外采集的高光谱测试数据小波变化和分维数，结果如表2所示。由表2可知，同一胁迫水平的小波分形维数值在不同尺度的波动范围有明显差异，重度胁迫和轻度胁迫在小波分解前5个尺度小波分维数波动值都在0.2以上，第5尺度后，同一胁迫水平的波动差异较小，其波动值位于0.1以内。不同胁迫水平的小波分形维数范围值在同一尺度有些许差异，但频率集中值差异不大。

本文引入拐点思想利用小波分形维数的方法，通过统计重构后小波分形值，在不同小波分解尺度上绘制离散小波变换的光谱效应图，图中发生显著变换的点所对应的的尺度判断为光谱特征尺度。丰富的光谱特征可以用小波分形维数表示，原始光谱曲线分形特征用小波变换的低频分形维数表示，高光谱细节特征用小波系数的高频分形维数表示。图4反映不同胁迫水平的小波分形维数在各分解尺度下的波动规律，不同胁迫水平的分维数随尺度增加均呈下降趋势，在小尺度下(1～4)重度胁迫和轻度胁迫的小波分维数值均波动较大，大尺度下(6～8)其值出现重叠，没有明显差异，但小波分形维数下降剧烈，在第8尺度时其值接近于1。结果表明，随着尺度增大，不同污染水平的小波分维数差异值逐渐减小，在第6尺度时，重度胁迫和轻度胁迫曲线的分维数值出现相交情况，之后6～8尺度下的小波分维数值相同且呈现同步下降趋势，结合表2结果分析，可认为小波分解第5尺度为分形维数的拐点，即明显区分胁迫状态的拐点。

表2 不同污染水平的水稻光谱各尺度小波分维数Table 2 Wavelet fractal dimension of rice spectrum at different pollution levels

图4 各尺度不同污染水平小波分形维数Fig.4 Wavelet fractal dimension of different pollution levels at different scales

2.5 水稻重金属胁迫特征尺度验证分析

为了检验定量计算的小波信息熵和小波分形维数在识别不同尺度下区分不同污染胁迫水平的可靠性和普适性，全面描述不同尺度对水稻重金属胁迫监测的影响，计算实测光谱曲线在不同尺度下的叶绿素光谱指数，图5为40个实测样本点的高光谱指数平均值随不同尺度的变化情况。由图5可知，MCARI/OSAVI、NDSI_R、Depth三类叶绿素植被指数在相同胁迫水平下随不同尺度的变化情况基本趋于一致，重度胁迫下的高光谱植被指数在第5尺度内总体呈上升趋势，第6尺度起三种植被指数值快速下降，轻度胁迫变化趋势存在差异，1～5尺度的指数值为波动性下降，之后高光谱植被指数值上升，两种不同污染胁迫水平下的 MCARI/OSAVI、NDSI_R、Depth的差值在第5尺度时最大，而在第4和第5尺度下重度胁迫的高光谱指数值均高于轻度胁迫，其他分辨率下两种胁迫水平的高光谱指数值表现为相反的趋势，即轻度胁迫高于重度胁迫。因此，可将这两个光谱分辨率作为水稻重金属胁迫水平区分的拐点，即可将第5尺度作为利用地面高光谱数据监测水稻重金属胁迫水平的特征尺度，同时，也很好地验证了基于小波信息熵和小波分形维数得出的特征尺度为第5尺度的结论。

图5 各光谱分辨率不同胁迫水平水稻叶绿素高光谱指数 Fig.5 Chlorophyll hyperspectral index of rice at different spectral resolutions

3 讨论

本研究基于地面成像光谱技术获得水稻不同污染胁迫状态下的详细光谱信息和每一波长下的反射率信息，以特征光谱尺度为识别目标，提出一种综合光谱曲线信息量和分形特征进行水稻重金属胁迫地面遥感技术特征光谱尺度识别的方法，为进一步采用合适光谱分辨率的高光谱卫星影像进行水稻重金属胁迫监测提供依据。

水稻重金属胁迫遥感监测特征光谱尺度是指在一定的谱段范围内能区分胁迫水平差异的最适尺度，高光谱数据中蕴含大量地物光谱信息，如何筛选重金属胁迫导致的作物光谱响应敏感波段对于光谱特征尺度的研究至关重要，本文通过研究不同胁迫水平下实测ASD水稻高光谱曲线发现波长在1 000～1 200、1 250～1 500、1 750～2 000和2 400～2 500 nm范围内的光谱反射率出现断层现象，这是由于在近红外和短波红外区域水汽吸收波段的强烈噪声的影响，任红艳等[13]以受铅胁迫的水稻为研究对象，发现水稻受到铅重金属污染的冠层反射高光谱敏感的特征波段在可见光520～560、690 ～710和760～810 nm处；田国良等[26]以受Cd胁迫的水稻为研究对象，发现监测Cd污染的最佳波段范围为520～580、610～680和730～800 nm。本研究结合已有结果和水稻生理生化特征考虑，最终选取A1、A2研究区域波长范围为450～900 nm为研究对象。

植被受重金属污染引起的光谱变化为弱信息，已有学者基于数学手段对重金属胁迫下农作物的原始光谱信号进行变换和增强处理，以提取重金属胁迫下隐含的光谱弱信息，如刘美玲等[27]利用db5小波系数可以探测到光谱奇异点从而能识别不同重金属污染胁迫水平;Liu等[28]利用小波分形方法对吉林省长春市城区受重金属污染的水稻光谱进行分析，发现对可见光/近红外波段范围的光谱曲线进行小波变换可以降低噪音并放大胁迫信息，小波分析是用一簇函数去逼近待分析的信号，其时域-频域局部性质不仅能有效地从原始信号中提取瞬态突变信息，而且能够有效地处理和分析多分辨率、多层次、多尺度问题，本研究基于弱信息方法从寻找水稻胁迫监测的特征尺度角度出发，利用小波变换的多尺度模拟的特点探究不同光谱分辨率的尺度效应，从而识别区分不同污染水平重金属胁迫特征尺度。

在基于小波变换的光谱特征提取中确定分解尺度是关键因素，尺度选择的过大或过小都不能获取较高的分类精度。本研究中的特征尺度识别研究受限于小波分解的层数，因不同分解层数下的波形曲线的波峰波谷特征差异明显。1维离散小波分解时要求分解尺度是2的整次幂，根据不同光谱尺度的长度设置相应的分解层数[25]。Bruce等[29]提出,确定最佳的分解尺度依据母小波函数和光谱信号的长度两者。依据母小波的滤波器长度确定了最优分解尺度，在本研究中选择进行多尺度分解的母小波是db5小波函数，小波分解的最佳分解尺度应为5，但是另一方面根据本实验所选的敏感特征光谱区间，最佳的分解尺度应为8。虽然母小波和光谱信号长度确定的最佳分解尺度不同，但是最终的分解尺度取值通常以两者之间的最大值为准，故本研究的最终小波分解尺度为8层。不同胁迫水平下的ASD光谱曲线经通过8层小波分解后生成的不同小波特征系数可近似模拟不同光谱分辨率工作模式，分解后的细节系数随分解层数增加递减，频率降低，可反映出光谱的多分辨率特性，与刘美玲等[30]研究结果一致。

本研究采用的特征尺度度量参数主要为小波信息熵和小波分形维数，小波信息熵可表示不同信号分量的信息大小，熵值大则所携带的信息量大，特征尺度对应信息量大的尺度，量化数据的随机性和不确定性。苏俊英[31]提出，丰富的光谱特征可以用小波分形维数表示，原始光谱曲线分形特征用小波变换的低频分形维数表示，高光谱细节特征用小波系数的高频分形维数表示。不同参数研究结果均表明,不同光谱曲线特征转折点出现在第5尺度，鉴于从第5尺度开始光谱曲线中具有特征意义的峰谷细节信息出现明显弱化，且在1～4尺度不同胁迫下的小波信息熵和小波分形维数区分较为明显，因此，本研究将尺度5即对应的光谱分辨率为32 nm时为水稻重金属胁迫特征尺度拐点。但本文的特征尺度是针对所选的特定特征而言的，不同尺度间污染胁迫表现与所选的特征密切相关，农田重金属污染最佳观测尺度的选择依然要综合考虑所采用的方法与具体的观测对象，特别需要结合重金属胁迫下的作物生长机理分析。