随着社会的发展，重金属污染与监测受到广泛关注。由于微生物等很难降解土壤及植物中的重金属离子，因而其会在土壤以及植物体内累积，当这种富集（如Cu²⁺等）达到一定的限度后就会对植物造成危害，危及人类的健康，甚至影响整个生态系统环境，因此如何监测、预防和治理重金属污染是当今社会面临的热点问题^[1-3]。常规的化学、物理化学、生物法等^[4]对重金属污染监测存在过程繁琐、处理成本高、二次污染、监测范围小等缺点，而高光谱遥感技术具有光谱分辨率高、监测范围广、可实时动态获取大量信息以及不接触不损害本体等优点^[5-6]，因此应用高光谱遥感监测植物中的重金属污染具有现实意义。重金属对土壤造成污染并被植物吸收后，会影响和破坏土壤和植物体内的微生物活性及其他结构，并引发植物产生一系列异常，尤其是叶绿素等生化成分含量会发生变化，从而引起光谱曲线变化^[7-8]，所以，对光谱畸变信息的局部特征差异性识别成为一种重金属污染监测的有效手段。王平等^[9]利用光谱吸收谷、红边位置（Red Edge Position，REP）等作为玉米重金属污染监测参数用于反演叶绿素含量，迟光宇等^[10]发现铜污染胁迫下小麦反射光谱会发生变化且红边“蓝移”，李庆亭等^[11]发现植被重金属含量与红边光谱特征等有显著相关性，王维等^[12]发现一阶微分光谱具有估算土壤中重金属铜含量的潜力，刘美玲等^[13]提出可利用高光谱高频组分的盒维数法（Box Dimension Method，BDM）分形特征监测水稻重金属污染等。以上成果都有较好的重金属污染监测效果，但有时对细微的光谱差异敏感性及适应性较差。

近年来，时频分析技术也常被用于高光谱数据分析和信息识别研究，如谐波分析（Harmonic Analysis，HA）技术主要应用于电力行业^[14]，但也有学者将其引用到了高光谱影像小目标探测^[15]等领域；小波变换常用于高光谱影像压缩^[16]、图像去噪^[17]、混合像元分解^[18]等；分形维数（Fractal Dimension，FD）也常被用在土壤颗粒组成及分形维数特征提取^[19]、矿物受力分析及筛选^[20]、高光谱数据特征提取^[21]等方面。而时频分析技术在重金属污染监测领域应用较少，且结合多种时频分析技术优势的综合应用研究就更少。本文拟根据重金属铜（Cu）胁迫下盆栽玉米培育实验及玉米叶片光谱测量等数据，采用谐波分析、经验模态分解（Empirical Mode Decomposition，EMD）、离散小波多层分解（Discrete Wavelet Multi-layer Decomposition，DWMD）等多种时频分析技术，在建立光谱变异信息甄别模型（Screening Model，SM）的基础上，结合量规法分形维数（Divider Method on FD，DMFD）计算方法构建SM-DMFD模型，从而实现不同Cu²⁺浓度胁迫下玉米叶片光谱的差异性变化特征提取与重金属污染信息探测。同时，通过与绿峰高度（Height of Green Peak，GH）、红边位置（REP）、红边最大值（Maximum-value of Red-edge，MR）、红边一阶微分包围面积（First-derivative Area of Red-edge，FAR）和盒维数法（BDM）监测方法应用结果的对比分析，以及监测方法模型应用结果的污染信息可视化表达，验证SM-DMFD模型在重金属污染监测方面具有较好的有效性和优越性。

1 材料与方法 1.1 胁迫实验与数据采集 1.1.1 植株栽培

实验采用“中糯1号”玉米种子进行盆栽培育，采用浓度梯度为0、250、500 μg·g^-1的CuSO₄·5H₂O溶液污染土壤，以模拟重金属Cu²⁺胁迫玉米生长，其胁迫梯度相应记为Cu（ck）、Cu（250）、Cu（500），每个胁迫梯度平行培育3组盆栽，共9组。先对玉米种子做催芽处理，两日后栽种，待玉米出苗后于统一时间向盆土中添加NH₄NO₃、KH₂PO₄和KNO₃营养液。育苗期间保持适宜的温度及湿度，定期养护和浇水。

1.1.2 光谱数据采集

玉米出苗后使用光谱范围为350~2500 nm的SVC HR-1024I地物光谱仪对玉米叶片进行光谱采集，其中测量时采用50 W的卤素灯进行照明，置视场角为4°的探头垂直于玉米叶片表面40 cm处，采集光谱的叶片放置于专用黑色硬板上，对采集后的光谱反射系数经专用白板标准化。测量时分别对不同胁迫梯度下每株玉米老、中、新三种代表性叶片进行光谱采集，最后计算每一胁迫梯度下玉米叶片的平均光谱。

1.1.3 Cu²⁺含量测定

对采集完光谱数据的玉米叶片进行冲洗、烘干、粉碎并装入样品袋封存。保存一段时间后经高纯硝酸、高氯酸消化处理，使用WFX-120原子吸收分光光度计测定Cu²⁺含量。检测方式：使用火焰检测器；仪器条件：波长324.7 nm，狭缝0.7 nm，负电压304 V，灯电流3.00 mA。在相同条件下，对每一胁迫梯度的叶片重复3次测量后取平均值作为该胁迫梯度下的叶片中Cu²⁺含量，如表 1所示。

1.2 谐波分析法光谱预处理

谐波分析（HA）是将时间t序列信息f（t）转换到频率域，并以多个频率不同的正（余）弦波相叠加的方式来表示时域分量的一种方法。高光谱遥感的光谱信息可看作是在一定波段范围内波谱反射率随波长变化而不同的连续信息数据，如果把波长变化看作时间序列，那么也可以用HA处理光谱数据，即将光谱视作总波段数为N的函数，则HA就可以把光谱信息多次分解成由一系列振幅、相位和余项等能量谱特征分量组成的正（余）相叠加的谐波。光谱经HA分解处理后，再通过光谱重构（即HA逆运算），也具有光谱平滑和去噪效果。若记光谱为V（s）=（v₁，v₂，…，v_N），相对应波长的波谱反射率记为v_s，s是波段序号（s=1，2，…，N），则光谱信息的h次HA分解式为：

式中h次谐波分解的能量谱特征分量为：

式中：s表示波长序列；A₀/2表示谐波余项；N表示波段数；h表示分解次数；A_h、B_h、C_h、φ_h分别表示第h次谐波分解的余弦振幅、正弦振幅、谐波分量振幅和谐波分量相位。

1.3 经验模态分解法模型参量提取

具有很高信噪比的经验模态分解（EMD）方法是由Huang等^[22]提出，是用来将信号中不同时间尺度的波动或变化形式逐级分解，进而产生有限个频率由大到小的本征模态函数（Intrinsic Mode Function，IMF）分量。EMD适合于分析非平稳、非线性信号序列，其分解的IMF分量包含了原始信号不同时间尺度的局部特征信息。光谱数据的EMD处理过程主要包括：

（1）任选一条波段序号为s的光谱V（s），确定光谱上所有极值点，依据这些极值点，用三次样条函数拟合光谱曲线的上、下包络线，并计算该上、下包络线均值E₁（s）。通过求取V（s）与E₁（s）的差产生新的数据序列H₁（s），再把H₁（s）视作新的V（s）。重复以上操作，直到H_i（s）满足关于波段序号轴的局部对称，则其变为从原始光谱中选出的第一阶IMF，记为IMF₁，一般情况下IMF₁包含了信号的最高频成分。

（2）将IMF₁从原始光谱中分离出来，进而得到去除高频成分的差值信号Q₁（s），再将Q₁（s）作为新的信号重复步骤（1），直到第n阶光谱模态分解的残余信号数据不能再选出第n阶的IMF_n分量为止。

1.4 SM-DMFD探测模型 1.4.1 光谱变异信息甄别模型（SM）

由于受污染的与健康的植被光谱经过HA去噪与EMD分解处理后，其光谱曲线形态仍具有很高相似度，光谱局部特征信息与差别仍难以甄别。为了从光谱中有效地提取出光谱变异的细微差别信息，可依据EMD分解出的IMF₁分量，构建光谱变异后局部特征信息探测的SM，即用SM对IMF₁分量进行如下处理：

式中：x表示原光谱向量；V（x）表示处理后的光谱向量。该模型最大的优点在于扩大了光谱曲线的局部差异。

1.4.2 离散小波多层分解

小波变换概念由法国工程师Morlet于1974年首次提出。小波变换很适合于非平稳的信号f（t）局部化分析与信息的局部特征提取，在低频部分具有较高的频率分辨率和较低的时间分辨率，在高频部分具有较高的时间分辨率和较低的频率分辨率。光谱信息是一种离散数据，所以对光谱进行小波分解，常采用离散的小波变换技术。离散小波变换（Discrete Wavelet Transform，DWT）指对尺度因子a和平移因子b进行离散化，其离散化方法可取幂级数形式，即：

式中：a₀≠1固定值，通常假定a₀>1；k为系数；b₀为变化因子。若信号f（t）的离散小波表示为：

则，信号f（t）的DWT系数C_j，k为：

式中C_j，k包括细节系数（高频）和近似系数（低频）。选择特定母小波对SM处理结果进行离散小波多层分解后，可提取其高频的细节系数序列，以便于求取量规法分形维数值。

1.4.3 量规法分形维数（DMFD）

分形指整体与局部以某种方式相似，其概念最早由数学家Mandelbrot提出^[23]，它的两个重要原则是自相似原则和迭代生成原则。分形维数（FD）是表征分形的主要参数，反映了复杂形体占空间的有效性。FD有盒维数法（BDM）、方差法、结构函数法、量规法（Divider Method，DM）等多种计算方法^[24-26]，本文首次引用DM进行FD计算，并用于植被重金属污染光谱时频分解后曲线分形维数求取。设光谱长度L（r）由尺子长度r和尺子测量的次数N（r）来决定，表达方式为：

根据Mandelbrot的研究，有下式成立：

式中：L（r）为被测光谱的长度；r为标度；M为待定常数；D为分形维数。对式（8）两边取对数，可得：

式中：C为lgM；K为斜率，K=1-D，D=1-K，即分形维数FD=1-K。对于一维光谱曲线，取1 < FD < 2。

1.4.4 污染信息探测的SM-DMFD模型

基于SM与DMFD所构建的SM-DMFD探测模型处理流程如图 1所示，具体过程为：

（1）先对原始光谱数据进行HA分解和重构预处理；

（2）光谱重构后进行EMD分解，提取其中的IMF₁高频分量；

（3）接着用SM对IMF₁分量进行处理，有利于识别局部差异信息特征；

（4）然后对SM处理后的向量数据进行离散小波多层分解，提取多个高频细节系数序列；

（5）分别对各高频细节系数序列计算DMFD值，根据其值与植物体内所测重金属元素含量间相关系数的大小，择取最优细节系数序列及其对应的DMFD值；

（6）最后，根据DMFD值的大小实现植被重金属污染信息探测。

2 结果与讨论 2.1 基于SM-DMFD模型的污染监测 2.1.1 玉米叶片光谱预处理与IMF₁分量提取

采用式（1）对玉米叶片光谱进行谐波分解，经反复试验得出，当分解次数达到20时去噪效果理想，因此当分解20次后再根据式（1）的反变换重构叶片光谱。以Cu（ck）光谱的谐波分解为例，比较HA前后的光谱曲线变化，如图 2所示，可见玉米叶片光谱得到了较好的平滑与去噪预处理。再对重构后光谱进行EMD处理，根据EMD处理过程提取IMF₁分量，各种不同Cu²⁺胁迫梯度下提取的IMF₁分量如图 3所示。

2.1.2 光谱变异信息甄别

运用式（3）对IMF₁分量进行SM处理，结果如图 4所示。由图 4可见，SM处理后，不同Cu²⁺胁迫梯度下玉米叶片光谱在不同波段的SM值具有较大的局部差异，且SM均变为正值更易发现光谱变异的局部特征细微差异和变化规律；同时，SM值曲线在1263~1373 nm和1522~1649 nm两个波段区间内变化明显，即随着Cu²⁺胁迫梯度增加，波段内峰值呈逐渐降低的趋势，此两个波段区间内SM峰值及其与叶片中Cu²⁺含量的相关系数R如表 2所示。从表 2可发现，两个波段范围内所得峰值与玉米叶片中Cu²⁺含量相关系数均达到0.95以上，认为SM值能明显突出不同Cu²⁺胁迫梯度下叶片光谱曲线的变化，可有效甄别光谱局部差异特征。

2.1.3 SM-DMFD模型的污染探测

由表 2发现，SM模型可对玉米叶片局部光谱波段差异做出有效甄别，为使模型具有充分的有效性和可靠性，需根据式（6），进一步对SM处理后的曲线（图 4）进行离散小波多层分解。小波变换的母小波有多种，如：Haar小波、Symlets小波族、Db小波族等，其中Db小波族常称DbN（N为阶数）。经过多次应用比较，采用Db10母小波对图 4中的数据进行离散小波3层分解，应用效果理想，其分解之后得到高频的细节系数序列为Cd₁、Cd₂、Cd₃，如图 5所示。将各细节系数序列经式（9）计算后可求取DMFD值，如表 3所示，其中R为DMFD值与叶片中Cu²⁺含量的相关系数。

由表 3发现，随着Cu²⁺胁迫梯度的增加，不同细节系数序列中的DMFD值有所变化，但Cd₃序列的DMFD值与玉米叶片中Cu²⁺含量的相关系数达到0.986 0，即相关性最高，且DMFD值随着Cu²⁺胁迫梯度的增大而增大。因此，采用Cd₃序列的DMFD值进行玉米重金属污染信息探测，其监测效果最优，能有效探测玉米重金属污染信息。

2.2 方法应用比较与拟合分析 2.2.1 不同方法应用结果比较与分析

为了验证SM-DMFD模型在甄测玉米叶片光谱变异信息与污染探测方面的有效性和优势，同时采用常规的绿峰高度（GH）、红边位置（REP）、红边最大值（MR）、红边一阶微分包围面积（FAR）等污染光谱响应参量和盒维数法（BDM）进行污染监测应用结果比较分析，GH、REP、MR和FAR的计算方法如表 4所示，BDM计算方法见文献[24]。不同Cu²⁺胁迫梯度下叶片光谱的SM-DMFD模型应用结果DMFD值和其他常规方法的计算值，及其与叶片中Cu²⁺含量的相关系数计算结果如表 5所示。

分析表 5发现，DMFD与玉米叶中所测Cu²⁺含量的相关系数达到0.986 0，优于其他监测参数。同时将SM-DMFD模型应用结果和GH、REP、MR、FAR、BDM计算结果与玉米叶片中Cu²⁺含量进行线性拟合分析，图 6的拟合结果显示，SM-DMFD探测的DMFD值拟合效果最好，其拟合判定系数R²最高，达到了0.972 3，高于其他监测参量与方法的判定系数。同时通过构建的拟合曲线线性方程式，可以推算出玉米叶片中Cu²⁺含量，进而为反演和预测重金属Cu对植被的污染程度奠定基础。分析结果也验证了SM-DMFD探测模型的有效性和优越性。

2.2.2 不同方法应用结果的可视化分析

对玉米叶片的Cu²⁺污染信息探测结果进行可视化表达，可直观地描述叶片所受Cu²⁺污染的程度。

（1）建立一种标准化色条，同时将Cu²⁺胁迫梯度与相应的GH、REP、MR、FAR、BDM、DMFD不同计算值进行向量化表达。

由此可得，所制定的标准色条向量BZ′=[0 50 100 150 200 250 300 350 400 450 500]；而Cu²⁺胁迫梯度向量XP′=[0 250 500]；绿峰高度计算值向量GH′=[18.6744 19.9967 20.4433]；红边位置计算值向量REP′=[715 715 700]；红边最大值计算值向量MR′=[0.7594 0.7620 0.7923]；红边一阶微分包围面积计算值向量FAR′=[34.1344 34.1317 34.4033]；盒维数法计算值向量BDM′=[1.1308 1.1305 1.1332]；以及SM-DMFD模型的应用计算值向量DMFD′=[1.3879 1.3980 1.4045]。

（2）对所表达的各向量进行归一化处理，使所得各向量内的元素值归化到0~1范围，归一化处理公式：

式中：y_i表示归一化后向量的元素值；x_i表示初始向量的元素值；x表示初始向量。

由式（10）对上述标准色条向量以及各玉米叶片污染监测方法的计算值向量进行归一化处理，得到各自的归一化向量为：BZ=[0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1]；XP=[0 0.5 1]；GH=[0 0.7475 1]；REP=[1 1 0]；MR=[0 0.0790 1]；FAR=[0.0099 0 1]；BDM=[0.1111 0 1]；DMFD=[0 0.6084 1]。

（3）采用Matlab R2013a中的imagesc函数，依据所建的标准化色条BZ，对归一化后的各向量中元素进行赋色处理。归一化的标准化色条BZ以及XP、GH、REP、MR、FAR、BDM、DMFD向量赋色处理结果如图 7所示。

图 7中，BZ作为标准色条，分别分割并赋置了clear（无胁迫影响）、0~100、100~200、200~300、300~400、400~500 μg·g^-1的不同Cu²⁺胁迫梯度标准色。XP为本次实验设置的Cu²⁺胁迫梯度0、250、500 μg·g^-1对应的色条，BZ与XP色条为参照色条，通过与BZ标准色条比对，可证实XP色条的合理性。再把DMFD、GH、REP、MR、FAR、BDM色条与BZ、XP色条进行比较，可发现DMFD、GH色条与Cu²⁺胁迫梯度具有较好一致性，而其他色条与Cu²⁺胁迫梯度的一致性较差。而针对DMFD、GH色条，在Cu²⁺胁迫梯度为250 μg·g^-1的左右区间颜色上，与BZ的标准色条比较，可看出DMFD值的赋色在300 μg·g^-1梯度以下，而GH值的赋色在300 μg·g^-1梯度以上，所以，认为SM-DMFD模型更能准确预测Cu²⁺对玉米叶片的污染程度，从而也验证该模型具有有效性和优越性。

3 结论

（1）随着土壤中Cu²⁺胁迫梯度的增加，玉米叶片中Cu²⁺含量逐渐升高，与Cu²⁺胁迫梯度呈正相关。

（2）基于建立的玉米叶片光谱变异SM模型以及污染信息的SM-DMFD探测模型，最终得到不同污染叶片光谱的DMFD值能够有效甄别光谱间差异信息，并能判别叶片的污染程度，且与叶片中Cu²⁺含量的相关系数R和拟合判定系数R²分别达到0.986 0和0.972 3。同时通过GH、REP、MR、BDM、DMFD监测方法应用结果比较和计算值的可视化分析，验证了SM与SM-DMFD模型应用于玉米Cu污染信息监测具有较好的有效性和优越性，可为重金属污染监测研究提供参考方法和技术支持。