岭北稀土矿区位于江西省赣州市定南县, 面积约213 km²。 该稀土矿持续进行稀土开采30余年, 在长期人工浸溶开采与自然高温多雨的双重作用下, 土壤退化问题突出, 属于国家级水土流失重点预防区。 选取岭北甲子背复垦矿点作为研究区域, 坐标: 115° 2'22″— 115° 3'16″E, 24° 58'21″— 24° 59'15″N, 面积约1.6 km²。 该矿点历经“ 池浸” 、 “ 堆浸” 和“ 原地浸矿” 等多种工艺开采, 遗留大量裸露的尾矿和废弃地。 稀土开采扰动引起的复垦地土壤特殊的物理化学性质变化及稀土开采所导致的环境胁迫, 致使矿区生态恢复进展缓慢。

由于稀土矿区废弃地以土壤沙化酸化及重金属污染为主要表现形式, 复垦植被长势较差, 成活率低, 生态恢复困难。 为此, 其复垦植被应具有易存活、 适应性强等特点。 故选取稀土矿区竹柳、 山茶、 桉树、 松树、 油桐、 红叶石楠6种生存能力强、 耐土壤贫瘠的典型复垦植被作为研究对象, 基于复垦植被在稀土矿区环境胁迫下的光谱特征变化展开研究。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 研究区地理位置及概况Fig.1 Location and overview of the study area

于2020年8月14日— 15日使用美国ASD (Analytical Spectral Devices) Field Spec4地物光谱仪进行叶片光谱测量, 采样间隔为1.4 nm (350~1 000 nm)和2 nm (1 001~2 500 nm)。 具体时间为11:00— 14:00之间, 期间天气晴朗、 无风无云, 植被叶片处于成熟阶段。 选取矿区的竹柳、 山茶、 桉树、 松树、 油桐、 红叶石楠六种复垦植被作为采集对象。 此外, 为获取复垦植被光谱与正常植被光谱的差异, 同时在距离甲子背矿点1km范围正常区域(即未进行稀土开采的区域), 相同环境下根据均匀分布原则进行正常环境生长条件下的叶片光谱采集。

测量时确保样本和参考白板处于同一环境, 仪器探头垂直向下, 与叶片保持垂直, 置于叶片上方。 每隔15 min校准波谱仪, 每个叶片样品测量10次, 取其平均值作为样本光谱反射率。 最终共采集复垦油桐92组、 复垦湿地松87组、 复垦红叶石楠68组、 复垦竹柳92组、 复垦山茶130组、 复垦桉树160组、 正常油桐35组、 正常湿地松36组、 正常红叶石楠33组、 正常竹柳36组、 正常山茶40组、 正常桉树43组。 图2为甲子背矿点复垦植被实地数据采集时拍摄, 该图反映出该矿点复垦植被根细枝稀, 叶片发黄。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 复垦植被实地采集Fig.2 Reclaimed vegetation field collection

野外实测叶片高光谱数据易受土壤环境与人为操作干扰, 特别是在稀土矿区地表大面积裸露的背景下, 光谱数据预处理在光谱分析过程中起着至关重要的作用。 预处理主要包括剔除异常光谱、 去噪和均值。 首先剔除异常光谱, 确定植被光谱曲线, 共剔除复垦植被37组, 正常植被14组。 其次去除噪声间隔, 在野外高光谱数据采集过程中, 植被对水汽和二氧化碳具有很强的吸收带, 对可见光和近红外区植被光谱产生重要影响^[10]。 因此, 最终选择350~1 350 nm的光谱波段进行研究, 将1 350 nm以上的波段进行剔除。 最后将所有合格叶片样本光谱数据的均值作为最终光谱反射率。

导数法是高光谱光谱变换中最常用的方法, 表示为原始光谱的变化率。 原始光谱经过一阶导数变换处理, 能够部分消除太阳高度角、 地形等引起的亮度变化。 有针对性的减弱稀土矿区特殊的裸露红壤背景对植被光谱造成的影响, 提升光谱特征区间的对比细腻度, 一阶导数变换公式如式(1)

D(λi)=R(λi+1)-R(λi)Δλi(1)

式(1)中, D(λ _i)为一阶导数光谱; R(λ _i+1)、 R(λ _i)分别为波长λ _i+1、 λ _i处的反射率; Δ λ _i为两波长距离。

一维曲线分形维数可以表示该波形复杂程度, 矿区环境胁迫会导致复垦植被的光谱曲线波形产生特殊变化。 假设有一条不规则形状的光谱曲线I, 用一个长度为ε 的尺子来衡量该光谱曲线I的长度。 不考虑光谱曲线上小于ε 的不规则部分时, 尺子的测量值为M_ε (I), 当ε → 0时, 如果对常数c和D有

Mε(I)~cε-D(2)

则该光谱曲线I的维数为D, 常数c可以看成光谱曲线I的D维长度, 将公式两边取对数得

logMε(I)~logc-Dlogε(3)

当ε 趋于零时, 有

D=limε→0logMε(I)-logε(4)

当计算某条光谱曲线分形维数时, 可以选取ε 值的一个恰当范围, 绘制logM_ε (I)和logε 双对数图, 其斜率即为光谱曲线的维数D。

离散小波变换(discrete wavelet transformation, DWT)是在一系列离散尺度上对信号进行分析的方法, 凭借多尺度特点和突出细化特征的优点被应用于估算植物叶绿素、 重金属胁迫程度等一系列实际问题^[11]。 “ Db5” 小波函数将植被光谱分解为低频和高频两部分, 低频部分被保留, 高频部分被阈值化, 然后继续在低频的部分进行小波分解, 得到低、 高频部分。 经过多层分解可以将原始信号分解为多个子信号, 分解公式如式(5)

f(λ)=aj(λ)+∑i=1jdi(λ)(5)

式(5)中, f(λ )为光谱信号, j为分解层, a_j为低频部分, d_i为高频部分。

短时傅里叶变换 (short time Fourier transform, STFT)将难以处理的时域信号转换成易于分析的频域信号。 其基本原理是通过窗函数将光谱曲线切分成许多相等的波长间隔, 用傅里叶变化分析每一个波长间隔, 得到频域信号的变化结果, 将所得结果平铺得到二维的表象, 进而通过分析得到波段上的频率特征, 最后实现光谱曲线空频局域化, 具体公式如式(6)

STFTf(ω, τ)=∫-∞+∞[f(t)g(t-τ)e-iwtdt](6)

式(6)中, f(t)为信号函数, g(t-τ )为窗函数, w为信号函数基频率。

此外, 实验通过Specgram函数来实现STFT, Specgram函数的语法如式(7)

[S, F, T, P]=spectrogram(A, window, noverlap, nfft, fs)(7)

式(7)中, A为植被光谱曲线; window为窗函数; noverlap为各段之间重叠采样的点数; nfft为计算离散傅里叶变换的点数; fs为采样频率; S为输入A的短时傅里叶变换; F是在输入变量中使用F频率变量; T是频谱图计算的时刻点; P是能量谱密度。

获取竹柳、 山茶、 桉树、 松树、 油桐、 红叶石楠的复垦植被和正常植被叶片高光谱数据, 对比复垦植被和正常植被的光谱曲线, 其结果图3所示。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 植被原始光谱均值反射率曲线Fig.3 Original mean spectral reflectance curves of vegetation

对照图3中经过预处理的六种复垦植被和正常植被光谱均值反射率曲线, 复垦植被和正常植被的光谱波动趋势基本一致。 400~700 nm(可见光波段), 植被叶片中的叶绿素强烈吸收入射光, 形成” 两谷一峰” [450 nm蓝光附近(“ 蓝谷” )和660 nm红光附近(“ 红谷” )以及550 nm绿光附近(“ 绿峰” )], 呈现出区别于其他地物独特的植被光谱特征。 但是无论在可见光还是近红外波段, 矿区复垦植被光谱呈现出相同波长下的反射率均高于其同种正常植被的现象。 有研究证明, 病害程度的增加或者环境胁迫的加剧, 会减少植被叶片叶绿素含量, 而当叶绿素含量降低时, 会导致该光谱区域的叶片反射率增加^[12], 这也是矿区环境下复垦植被光谱高于同种正常植被的主要原因。 但正常植被与复垦植被光谱曲线都具备植被基本特征, 虽然存在部分差异, 仅基于原始光谱无法获取丰富的有效信息。

导数光谱处理本质上反映了植被中各生理成分的吸收特性, 它可以消除背景噪声的干扰, 并对重叠峰进行分解。 对六种不同胁迫环境下的植被作一阶导数变换, 并根据其光谱特征进行分析。 基于原始光谱一阶导数的植被光谱特征对照结果如图4所示。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 植被的一阶导数光谱均值反射率曲线Fig.4 First derivative mean spectral reflectance curves of vegetation

植被体内的叶绿素、 含水量和其他物质的光谱吸收, 在可见光和近红外波段形成了“ 蓝边位置” (490~530 nm最大导数值对应的波长)、 “ 黄边位置” (560~640 nm最大导数值对应的波长)、 “ 红边位置” (680~760 nm最大导数值对应的波长)等特征光谱变化相关的区域, 这些光谱特征区域是与植被区别于其他地物的独有性质, 其中应用最多的是“ 红边位置” 。 当植被受到环境胁迫或病虫害等因素的影响时而失绿, “ 红边位置” 将会往短波蓝光的方向偏移(即“ 蓝移” ), 当植被生长旺盛, 生物量高“ 红边位置” 将会往长波红外的方向偏移(即“ 红移” )^[13]。 矿区植被光谱“ 红边位置” 变化情况分别为: 正常竹柳-复垦竹柳(724~718 nm)、 正常油桐-复垦油桐(718~701 nm)、 正常松树-复垦松树(718~718 nm)、 正常山茶-复垦山茶(718~700 nm)、 正常红叶石楠-复垦红叶石楠(719~701 nm)、 正常桉树-复垦桉树(706~701 nm)。 除湿地松外, 大部分植被均出现了“ 红边位置” 的“ 蓝移” 现象, 其中复垦红叶石楠、 油桐和山茶出现了尺度较大的红边“ 蓝移” 。 说明复垦红叶石楠、 复垦油桐和复垦山茶在矿区环境下受到的胁迫较为严重。 掌握了能反映出复垦植被受矿区环境胁迫强弱的信息, 有助于实施针对性治理, 为矿区综合治理提供科学依据。

分形维数通常被用来描述复杂图形的整体特征。 上述的光谱特征参数包括植被指数与“ 三边” 参数(“ 三边” 参数是指基于光谱位置特征的相关变量, 即红边、 蓝边和黄边)都只针对某个波段范围内的光谱曲线, 而部分研究表明地物光谱曲线也具有分形特征^[14], 因此引入分形理论综合描述矿区复垦植被和正常植被的光谱差异, 分形维数结果如表1所示。

表1

Table 1

表1(Table 1)

表1 植被分形维数 Table 1 Vegetation fractal dimension 植被/分形 维数D 桉树 红叶 石楠 山茶 松树 油桐 竹柳 复垦植被 0.921 1 0.917 6 0.905 8 0.900 4 0.921 9 0.912 6 正常植被 0.900 8 0.901 3 0.898 9 0.900 0 0.901 3 0.902 2

表1 植被分形维数 Table 1 Vegetation fractal dimension

通过计算矿区植被光谱曲线的分形维数, 表1中用D表示, 可以看出同种植被在不同生长情况条件下, 呈现出复垦植被分形维数高于正常植被的规律, 其中复垦油桐和正常油桐的分形维数差异最大。 不同于其他光谱参数计算, 分形维数针对全段光谱曲线, 可以综合反映植被光谱总体特征。 根据分形维数的大小一定程度可以判别同种复垦植被与正常植被, 为矿区植被识别与生态恢复提供参考依据。

用Daubechies小波函数处理样品的光谱曲线, 探测植被光谱在环境下的弱胁迫响应。 将六种植被不同环境胁迫下原始反射以及一阶导数光谱进行离散小波变换, 如图5所示。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 植被原始光谱和一阶导数光谱不同分解尺度离散小波变换波形图 (a)— (f): 植被原始光谱小波变换; (g)— (l)为植被一阶导数(dv1)光谱小波变换Fig.5 Discrete wavelet transforms at different decomposition scales for the original and first-derivative spectra of vegetation (a)— (f): The wavelet transforms of the original spectra of vegetation; (g)— (l): The wavelet transforms of first-derivative spectra of vegetation

一维离散小波分解的分解尺度必须为2的整数幂。 如图5所示, 1— 8层的分解尺度分别为2、 4、 8、 …、 256 nm, 最终在1— 8层产生8个小波细节系数分别为d1— d8。 从d1— d8, 随着分解层数的增加, 呈现出光谱信号频率逐渐降低, 小波细节系数减少, 且小波细节系数对应的光谱特征减弱的趋势。

图5(a)— (f)为植被叶片原始光谱离散小波变换, 在小波细节系数d1— d4范围内噪声较多, 峰谷特征不明显。 从d5开始, 不同环境生长同种植被光谱曲线开始出现明显的差异, 复垦植被和正常植被的峰谷特征在“ 红边位置” 出现交叉。 但小波细节系数d6— d8范围内, 光谱曲线逐渐平滑, 极值交叉点越来越少, 所包含的信息也就越少。 最终原始光谱离散小波换最佳细节系数为d5。

图5(g)— (l)为植被叶片一阶导数光谱离散小波变换, 不同于原始光谱离散小波变换, 在小波细节系数d1— d5范围内存在较多的噪声, 直到d6复垦植被和正常植被在“ 红边位置” 才出现明显的峰谷特征位置交叉现象。 小波细节系数d7-d8范围内光谱曲线的奇异性开始减弱, 特征点开始减少, 最佳细节系数为d6。

与原始光谱离散小波变换相比, 在光谱特征数量相近的情况下, 一阶导数光谱离散小波变换在更小的尺度下放大了光谱特征细节差异, 取得更好的效果。 叶片光谱离散小波变换有利于在多尺度上放大光谱细节特征, 对植被在矿区环境胁迫下的特征为矿区植被生长准确监测提供理论依据。

短时傅里叶变换(STFT)又称窗口傅里叶变换, 短时傅里叶变换可以全面把握光谱回波的信号, 将空间域和频域结合起来, 通过空频分析同时分析光谱信号在空间和频率的变化情况, 将一维植被光谱曲线转换成二维傅里叶图像, 矿区六种植被叶片原始和一阶导数的短时傅里叶变换光谱如图6所示。

图6

Fig.6

Figure Option ViewDownloadNew Window

图6 植被原始和一阶导数的短时傅里叶变换光谱 (a)— (f): 植被原始光谱傅里叶变换; (g)— (l): 植被一阶导数光谱傅里叶变换, 其中每个子图中从上到下顺序为植被原始光谱, 植被光谱频谱, 正常植被空频图, 复垦植被空频图Fig.6 Short-time Fourier transform spectra of the original and first-derivative spectra of vegetation (a)— (f): Fourier transforms of vegetation raw spectra; (g)— (l): Fourier transforms of vegetation first-derivative spectra. In each subfigure, from top to bottom, original/first-derivative spectra, frequency map, spatial frequency map of normal vegetation and spatial frequency map of reclaimed vegetation

如图6所示, 我们将矿区复垦植被和正常的原始光谱以及一阶导数光谱进行短时傅里叶变换, 短时傅里叶变换在空间域将光谱信号截取为多段, 对每一段光谱信号做傅里叶变换, 并计算出频率特性, 继而组合完成空间-频率域的二维图像(空频图)。 经过STFT处理后的光谱信号具有空间域和频域的局部化特征。

在图6(a)— (f)为植被叶片原始光谱短时傅里叶变换, 各原始植被光谱曲线所对应频谱图具有随波长增加其功率谱下降的趋势, 全波段内波形起伏转折较多, 主要集中在350~750 nm(可见光-近红外波段)。 这些空频特征基本都出现在植被原始光谱的“ 红边” (700~750 nm)与中红外第一个“ 波谷” (900~1 000 nm)处, 这两个位置光谱反射率相近, 但是通过短时傅里叶变换, 在空频图上实现局域化, 突出了“ 红边” 与中红外第一个“ 波谷” 处的空频特征, 放大了复垦植被与正常植被的特征差异。

在图6(g)— (l)为植被叶片一阶导数光谱短时傅里叶变换, 各植被一阶导数光谱曲线所对应频谱图较原始光谱波形变化在全波段内更加剧烈。 其空频图整体波长归一化频率高于原始植被光谱, 并且出现了更多的空频特征波段。 一阶导数光谱经过短时傅里叶变换后, 整体空频特征波段增加, 同时可以观察到波长归一化频率的变化幅度也更大。 一阶导数在更小的尺度上, 更多的波段放大并增加了空频特征, 对于环境胁迫下的矿区植被识别与特征分析提供依据。