1 引言(Introduction)

我国频繁发生的重金属污染事件已成为现今社会的重大环境问题.近年来, 云南阳宗、山东沂南、湖南浏阳、陕西凤翔、福建上杭、河南济源、江苏盐城、广东清远、广西龙江等地爆发的几十起重特大重金属污染事件令全国震惊, 给当地生态与人体健康等带来了严重危害.铜(Cu)虽然是人类生活中广泛应用的一种重金属, 但矿山开采、工业生产甚至农业施肥用药等人类活动所产生的排弃物中常含有铜离子(Cu²⁺), 并能在土壤中富集, 如果植物从土壤中吸收了过量的Cu²⁺就会对其生长发育产生直接影响, 并通过食物链危害人体健康.因此, 重金属污染监测与治理方法研究已成为当今学术科研的一个热点(Vincent et al., 2012；Masoudi et al., 2012；Nayan et al., 2012；迟光宇等, 2006), 其中, 针对植物污染状况监测及植物体内重金属含量估算等是其研究的重要内容(朱叶青等, 2014；任红艳等, 2010).高光谱遥感因具有较高的光谱分辨率、可实时动态大面积监测、能进行理化性质反演等优点而成为重金属污染监测的重要技术手段, 也克服了传统定点采样监测耗时费力、繁杂及具有破坏性等缺点.

植物受重金属的污染程度可通过反射光谱特性表现出来(李庆亭等, 2008；刘素纯等, 2006；刘厚田等, 1984).近年来对植被重金属污染的监测方法研究越来越多, 如刘美玲等(2010)运用小波变换分析水稻受重金属污染的光谱奇异性；张龙等(2013)将近红外光谱与模式识别技术相结合构建监测重金属铅、汞和镉对水稻的胁迫模型；此外, 还有学者研究了重金属胁迫程度与植被指数、绿峰高度、红谷深度、红边最大值、红边面积、蓝边等光谱特征参量的关系(杨可明等, 2016), 如植被受锌胁迫程度与NDVI、叶片叶绿素含量呈负相关(Schuerger et al., 2003)；植物受重金属污染程度越小, “红边位置”会出现“红移”现象, 反之出现“蓝移”, 红边最大值也随之减小(梁守真等, 2010；谢晓金等, 2010；Bianchi et al., 2015)等.以上研究中, 小波变换因具有光谱时间频率的局部化分析能力, 而常被用于光谱变异的局部特征信息分析与提取, 因而在植物重金属污染信息的高光谱遥感探测方面具有优势.目前在光谱分析中, 小波变换主要用于光谱去噪(Dubot et al., 2015)、光谱压缩(田高友等, 2006)、化学信息提取(Zhang et al., 2015)等, 而在植被重金属污染监测方面应用较少.

本文拟基于不同浓度Cu²⁺胁迫梯度下的玉米盆栽实验, 依据盆栽玉米叶片的光谱测量数据与叶片中Cu²⁺含量的测试数据, 采用小波分析与奇异熵理论相结合以扩大光谱变异信息差异的方式, 构建小波奇异熵(Wavelet Singular Entropy, WSE)的玉米Cu²⁺污染监测模型, 并通过与绿峰高度(Green-peak Height, GH)、红谷吸收深度(Red-valley Depth, RD)等常规监测方法应用结果进行比较, 以验证WSE方法的监测效果.

2 理论与方法(Theories and methods) 2.1 光谱一阶微分

光谱微分技术能够消除部分大气效应、植被环境背景(阴影、土壤等)的影响, 反映植物的本质特征.已有研究表明, 对植物反射光谱进行一阶微分处理后, 其微分光谱值与植被体内生化成分的相关性会更高, 因此, 采用光谱一阶微分技术对植被光谱曲线进行微分处理, 可达到增大光谱间差异的效果, 光谱一阶微分公式见式(1).

(1)

式中, λ_i和λ_i+1为波长；R(λ_i)和R(λ_i+1)是波长λ_i、λ_i+1对应的光谱反射率值；R′(λ_i)是波长λ_i处对应的微分值；i=1, 2, 3, …, n(n为总的波段数).

2.2 小波理论

小波分析是在傅立叶变换基础上发展起来的.傅立叶变换虽在频域上具有任意局域性, 但无法提供时域上的任何局部信息.而在实际光谱分析中希望能同时了解在时域和频域上的光谱局部范围内信息, 小波的出现使此问题迎刃而解.小波变换是一种时频变换, 能通过伸缩平移运算对光谱逐步进行多尺度细化, 最终达到高频处时间细分, 低频处频率细分, 可聚焦到光谱的任意细节, 从而能自动适应反射光谱分析要求(Cheng et al., 2010).

小波变换可时-频分析光谱的局部特性, 能将高光谱信息分解成不同尺度的小波系数.对于平方可积函数f(t)∈L²(R), 如果用波长λ代替时间变量t, 则光谱曲线的连续小波母函数可表达为式(2).

(2)

式中, a表示伸缩因子, b表示平移因子.光谱经小波分解后, 得到的结果是小波系数.小波系数是伸缩因子a和平移因子b的函数：

(3)

式中, f(λ)为光谱反射率数据；n为波段数；小波系数W_f(a_i, b_i) (i=1, 2, …, m; j=1, 2, …, n)是由分解尺度(i=1, 2, …, m)和波段数(j=1, 2, …, n)组成的m×n二维矩阵.当检测高频信息时, 时窗会自动变窄, 提供分辨率高的波长信息；对于低频信息, 时窗会自动变宽, 从而符合理想时-频分析的要求.

2.3 奇异值分解理论

线性代数中, 奇异值分解(Singular Value Decompositi-on, SVD)是一种重要的矩阵分解, 是矩阵分析中正规矩阵酉对角化的推广, 在统计学、信号处理等领域应用广泛.设A为p×q阶复矩阵, 则存在p阶酉阵U和q阶酉阵V, 使得任何p×q矩阵A的奇异值矩阵λ满足式(4).

(4)

式中, U和V分别是p×p阶和q×q阶正交矩阵；λ=diag(α₁, α₂, α₃, …, α_s) (s=min(p, q))是对角矩阵, 其对角元素是按降序排列的奇异值.SVD的物理意义在于, 若矩阵A表示信号的时频信息, 那么奇异值矩阵λ则表示了A的基本特征.因此, 对光谱数据的小波变换系数矩阵进行奇异值分解后, 得到的λ能简捷反映光谱局部信息的时频分析特征.

2.4 小波奇异熵

信息熵通过对数据的统计分析来量化反映信息的不确定性和复杂度.奇异熵是奇异值理论和信息熵的有机结合, 能反映光谱变化的奇异性.对奇异值矩阵进行统计分析, 若得到的熵值越大, 则说明光谱的突变奇异性越大, 反之越小.小波奇异熵(WSE)是小波变换、SVD理论和信息熵原理的有机结合, 最终能对被分析光谱的突变性程度给出确定的量度, 因此, WSE能区分受不同程度的植物重金属污染.WSE定义如式(5)所示.

(5)

式中, WSE(k)表示小波分解1~k层的小波奇异熵总和, △p_i为增量小波奇异熵, 由式(6)计算得到.

(6)

式中, α_i表示对角矩阵上的奇异值.

3 实验部分(Experimentals)

选取玉米种子“中糯1号”为实验对象, CuSO₄·5H₂O作为分析纯溶液；测定设备包括波普范围350~2500 nm的美国SVC HR-1024I高性能地物光谱仪、WFX-120原子吸收分光光度计.

3.1 植株培养

采用有底漏花盆进行玉米种子发育.胁迫实验设置0、100、200、300、400、500、600、700、800和900 μg·g^-1 10个剂量的Cu²⁺胁迫梯度, 其相应梯度分别标注为CK(0)、Cu(100)、Cu(200)和Cu(300)等, 每个胁迫梯度平行培育3组盆栽, 共30盆.2016年5月6日进行玉米种子催芽, 5月8日在含有不同CuSO₄·5H₂O浓度的盆栽土壤中种下玉米种子, 5月9日在出苗后的各盆中加入适量的等量NH₄NO₃、KH₂PO₄和KNO₃营养液.培育期内保持各盆栽生长条件一致, 大棚通风良好并定期给玉米浇水.因400~900 μg·g^-1的Cu²⁺胁迫梯度下玉米植株未能正常生长发育, 故本次实验仅采用CK(0)、Cu(100)、Cu(200)和Cu(300) 4个胁迫梯度作为研究对象(表 1).

3.2 反射光谱测量

2016年7月17日使用SVC HR-1024I地物光谱仪采集玉米叶片反射光谱, 为了准确测定叶片光谱数据, 光谱测量时使用光源为50 W卤素灯, 探头视场角为4°并垂直于叶片表面40 cm, 叶片放于黑色硬板片上, 用专用白板标准化光谱反射系数.对不同浓度Cu²⁺胁迫梯度的玉米叶片分别进行光谱测量时, 每一胁迫梯度的玉米叶片光谱是由3条原始扫描输出光谱自动平均获得.

3.3 Cu²⁺含量测定

光谱采集完成后, 对不同胁迫梯度的玉米叶片进行冲洗、烘干、筛选、粉碎等处理并分装保存.2016年9月16日, 经高氯酸、高纯硝酸消化处理后, 采用原子吸收分光光度计测定玉米叶片中的Cu²⁺含量.仪器条件：波长324.7 nm, 狭缝0.7 nm, 负电压304 V, 灯电流3.00 mA；分析方式：火焰原子吸收.在相同条件下, 每一胁迫梯度样品进行3组平行实验, 最后取平均得到叶片中Cu²⁺含量(表 1).

4 结果与分析(Results and analysis) 4.1 胁迫梯度与光谱特征参数的污染关系分析

表 1表明, 健康的与受Cu²⁺污染的玉米叶片中Cu²⁺含量存在显著性差异, 且随着土壤中Cu²⁺含量的增加, 叶片中的Cu²⁺含量也随之增加.相关研究表明, 植被光谱在1300~2500 nm波长范围内受水汽强吸收、噪声等影响较大, 不适合用于监测植被重金属污染(修丽娜等, 2011).因而选取紫谷、蓝边、绿峰(GH)、黄边、红谷(RD)、红边(Red Edge, RE)及红肩等光谱特征参数(表 2), 用于分析玉米叶片污染光谱在380~1100 nm波长范围内的波形及变化特征.由表 2各参数的计算结果可以看出, GH、RD和RE参数值与叶片中测定的Cu²⁺含量呈正相关, 而紫谷、蓝边、黄边及红肩参数值与叶片中测定的Cu²⁺含量无明显相关性.

根据表 2, 用于计算GH、RD、RE参数的波长范围内光谱与一阶微分曲线如图 1所示.GH参数反映了植物在可见光范围内对绿色光的反射作用, 由图 1a可知, 随着铜胁迫梯度增加, GH位置从552.6 nm偏移到558.1 nm, 发生5.5 nm“红移”, 说明受Cu²⁺污染的玉米叶片对绿光反射作用降低, 导致了GH位置发生“红移”现象.由图 1b可知, RD位置基本没有偏移.RE是玉米受Cu²⁺污染在可见光范围内的主要反射特征, 由图 1c可知, RE位置从698.5 nm偏移到694.5 nm, 发生4 nm“蓝移”现象.

4.2 胁迫梯度与WSE的相关性分析

对玉米叶片的光谱进行小波分析后, 再分解其高频信息的奇异值, 通过计算获得奇异值的信息熵, 得到小波奇异熵(WSE), 从而可寻找玉米叶片光谱受污染影响的畸变性, 以反映Cu²⁺对玉米生长的胁迫程度.在光谱曲线500~600 nm(对应GH)、600~720 nm(对应RD)和670~780 nm(对应RE)波长范围内, 根据式(5)分别计算不同浓度Cu²⁺胁迫梯度下玉米叶片的WSE值, 结果如表 3所示.

从表 3可以看出, 在500~600 nm和600~720 nm波长范围内WSE值与胁迫梯度显著相关, 即随着Cu²⁺含量的增加, 两波长范围内的WSE值也都随之增加, 这表明玉米生长受胁迫梯度越大, 其叶片光谱变异程度也越大.并且在500~600 nm范围内WSE值与胁迫梯度相关系数r为0.8948, 在600~720 nm范围WSE值与胁迫梯度相关系数r达到0.9643;而在670~780 nm范围WSE值与胁迫梯度呈负相关, 相关系数r为-0.4624, 也即该波长范围内的WSE值不适合用于玉米叶片的污染监测.从WSE值与胁迫梯度的线性关系图(图 2)上来看, 600~720 nm波长范围内所计算的WSE值与胁迫梯度线性关系的拟合程度更高, 叶片光谱变化信息更为明显.

4.3 叶片中Cu²⁺含量与WSE、GH、RD的相关性分析

对不同胁迫梯度下玉米叶片中的Cu²⁺含量与WSE、GH、RD和RE值进行统计分析, 所得到的相关系数(r)如表 3所示, WSE值与叶片中Cu²⁺含量的线性关系如图 3所示.由表 3和图 2可见, 在RE所对应的670~780 nm波长范围内WSE值与玉米胁迫梯度无明显相关性；由图 3可以看出, 在该波长范围内WSE值与玉米叶片中的Cu²⁺含量线性关系的拟合程度较低, 因此, 不再分析RE值与叶片中Cu²⁺含量的相关性.

由表 3可以看出, 随着叶片中Cu²⁺含量增加, WSE、GH、RD值都随之增加, 其中, WSE值与叶片中Cu²⁺含量在500~600 nm范围内相关系数r为0.9817, 在600~720 nm范围内相关系数r为0.9824, 相关性显著.从图 3来看, 随着叶片中的Cu²⁺含量增加, 在500~600 nm和600~720 nm波长范围内所计算的WSE值与叶片中的Cu²⁺含量线性关系显著, 叶片光谱发生变化幅度较为明显, 光谱信息的复杂度增加, 所以WSE值变大.可见, WSE与玉米叶片中的Cu²⁺含量密切相关.

5 结论(Conclusions)

1) 随着铜胁迫梯度的增加, 玉米叶片中的Cu²⁺含量也相应增加.

2) WSE值与铜胁迫梯度呈正相关, 当胁迫梯度增加时, 在500~600 nm和600~720 nm波长范围内的WSE也随之增加, 并与胁迫梯度呈正相关, 其最高相关系数r达到0.9643, 可以证明WSE值大小与胁迫梯度有关.

3) WSE值与玉米叶片中Cu²⁺含量的相关系数值可达0.98以上, 叶片中的Cu²⁺含量越高则WSE值越大, 光谱的复杂度越高.

4) 通过3种波长范围的WSE值与所对应波长范围内绿峰(GH)、红谷(RD)和红边(RE)值的对比分析, 得出GH和RD对应波长范围内的WSE值效果较好, 与GH、RD具有几乎同等的污染监测效果, 表明WSE可以作为玉米Cu²⁺污染监测的一种新方法.