于 庆,吴泉源,姚 磊,徐夕博,周 旭
(山东师范大学 地理与环境学院, 山东 济南 250358)
在我国水资源相对匮乏的北方地区,有利用处理后的工业污水灌溉的传统。20世纪60年代以来,北方地区的污水灌溉面积迅速扩大,约占全国污水灌溉总面积的90%。由于污水中含有丰富的有机质悬浮物,采用污水灌溉可以缓解水资源短缺、节省肥料、提高土壤肥力,然而长期的污水灌溉必然会导致Pb、Cr、Cd、Hg等重金属在土壤中累积[1],最终导致污水灌溉的农作物受到重金属的污染[2-3]。冬小麦是北方地区的主要粮食作物,产量约占全国小麦总产量的56%,冬小麦中重金属的富集关系食品安全和人体健康[4-7]。因此,在污水灌溉区进行冬小麦重金属监测具有重要意义。
传统的重金属污染检测是实地采样后进行实验室分析,具有精度高的特点,但工作量大、成本高,检验点不连续,只能做到以点代面,不能推广到大范围,极大地影响了农业决策的全面性、时效性和客观性[8]。高光谱遥感技术使快速进行农作物重金属无损监测成为可能,具有监测面积大、易推广、工作量小的特点,尤其是冠层高光谱[9-12],其不但可以获取作物体内生化组分信息,而且还能反映作物下伏土壤的理化性质。目前,利用高光谱对农田重金属污染进行监测的研究很多[13-23],但多是对作物下伏土壤重金属的高光谱监测,并未直接对作物体内的重金属含量进行反演估算,而重金属在作物中的累积会因各种因素有所差异,作物下伏土壤的重金属含量并不能间接反映作物自身重金属的含量。
本研究选择山东省的典型污水灌溉区龙口市,在龙口市北部平原污水灌溉区设61个采样点,野外实测污灌区冬小麦冠层反射光谱,结合小麦重金属含量数据,运用逐步多元线性回归(Stepwise multiple linear regression,SMLR)以及偏最小二乘回归(Partial least squares regression,PLSR)的方法建立反演模型,选取各种重金属含量的最优反演模型,并通过反演及插值,分析该污水灌溉区冬小麦重金属含量的空间分布特征,以期为该污水灌溉区农作物重金属的实时监测提供模型理论基础,并为该污水灌溉区冬小麦重金属胁迫诊断提供理论依据。
本研究的污灌区位于山东省龙口市(120°13′14″~120°44′46″E、37°27′ 30″~37°47′24″N)北部平原区,属温带季风性气候,冬无严寒,夏无酷暑,气候宜人,总面积409 km2,是胶东半岛的主要粮食产区,农业生产条件发达,主要种植冬小麦和夏玉米等作物。该区污水灌溉历史已有30 a,铝金属冶炼、塑胶、钛业化工等工厂排放的污水以及生活废水经龙口市污水处理厂处理后直接用于农田的灌溉,使得土壤中重金属含量不断积累,重金属主要为Cr、Ni、Pb、Zn、Hg和Cd 6种。
冬小麦冠层光谱的采集采用美国ASD (Analytical spectral device)公司的Field Spec HH便携式手持地物光谱仪,采样间隔为1 nm,光谱分辨率为3 nm,其波长为325~1 075 nm,采样时间为冬小麦的拔节期和灌溉最佳时期的4月中旬(2016年4月13—15日),采样时段是晴朗无风的中午(11:00—14:00)。根据本研究区的土地利用现状图以及实地勘察,在本研究区的冬小麦种植区随机设立61个样点(图1),采集61个样点处的冬小麦冠层光谱数据并随机采集41处冬小麦叶片,将其带回实验室化验重金属含量。样点坐标采用GPS载波相位差分技术RTK方法进行定位。
在进行冠层光谱测量前,先利用白板进行校正,将光谱仪器探头垂直向下,距离冬小麦冠层顶部0.5 m,视角为8°,光谱测量时,每个样点重复测量10次。在ViewSpecPro软件中将每个样点重复测量的10条光谱中与其他光谱曲线有明显区别的曲线去除掉,再将剩下的光谱曲线取平均值得到该采样点实际光谱的反射率。为了从野外高光谱数据中提取有效光谱信息,需对野外获得光谱信息进行平滑处理,并用9点加权移动平均法去除噪声,原始高光谱经处理后可以消除背景噪声,增强差别,突出光谱特征[24],见图2。
图1 污灌区冬小麦冠层光谱数据采样点分布情况
图2 去除噪声以及平滑处理后的光谱
光谱采集后,将采样点处的冬小麦连根拔起,装入样品袋中带回实验室分析。将带回实验室的小麦样品使用去离子水洗净后放入80 ℃烘箱中烘干,将叶子剪下,放入粉碎机进行粉碎,冷却后用1∶1的HNO3溶解灰样,蒸干,用0.1 mol/L的HNO3定容,将定容完成的植物样品溶液用氢化物发生-原子荧光光谱法(HG-AFS)检测Hg含量,用电感耦合等离子体质谱法(ICP-MS)检测Cd含量,用电感耦合等离子体原子发射光谱法(ICP-OES)检测pb、Zn、Ni、Cr含量。
原始光谱反射率(Raw spectra reflectance,R)与各种重金属含量的相关性较差,对重金属含量信息的反映很微弱,而低阶微分处理能平缓背景干扰产生的影响,使光谱数据的轮廓变得更加清晰,提高原始光谱分辨率并且具有放大微小峰值的效果[25]。光谱参数(Spectral parameters,SP)也广泛应用于植被光谱建模反演研究中,它可以综合多个敏感波段的特征,加强对信息的提取能力,避免单一波段的偶然性和不精确性[26]。因此,本研究将原始光谱反射率、反射率一阶微分 (First derivate reflectance,FDR)、反射率二阶微分(Second derivate reflectance,SDR)和光谱参数分别作为反演指标与原始反射光谱进行建模对比,以期得到高精度的反演模型。
1.4.1 低阶微分指标 将原始光谱反射率进行低阶微分处理,反射率一阶微分、反射率二阶微分见图3。反射率一、二阶微分计算公式如公式(1)和公式(2)所示:
(1)
(2)
其中,R(λi+1)为第i+1条波段的反射率,R(λi-1)为第i-1条波段的反射率,Δλ为光谱的波段间隔。
1.4.2 光谱参数 植物在可见光波段具有很强的光吸收能力,根据不同波长的光吸收强弱变化形成波峰和波谷,因此,可见光波段是研究的重点光谱区域[27-29]。本研究在可见光区域内提取光谱参数,提取的光谱参数有各种光谱位置参数、光谱面积参数。各种参数的计算方法见表1。
偏最小二乘回归法是一种多对多的回归建模方法,比较适合处理变量多而样本数少的问题[30],而逐步多元线性回归法是常用的统计建模方法,可以解决数据维数带来的波段冗余和高光谱数据繁多的问题[31],本研究样本较小且光谱数据冗余,因此,选择偏最小二乘回归法及逐步多元线性回归法进行建模,并对比分析,旨在选出各种重金属的最优反演模型。本研究在对4种光谱指标与6种重金属含量进行相关性分析的基础上,分别以各种光谱指标为自变量,重金属含量为因变量进行建模分析,建模过程如下。
图3 原始光谱反射率微分变换
类型参数符号计算方法光谱位置参数红边位置λr红边内(680~750 nm)最大一阶微分对应的波长位置红边峰值Dr红边内(680~750 nm)光谱一阶微分的最大值红谷深度Ro红边内(680~750 nm)包络线去除后光谱的最小值黄边位置λy黄边内(550~650 nm)最小一阶微分值对应的波长位置黄边峰值Dy黄边内(550~650 nm)光谱一阶微分的最小值蓝谷位置λb蓝边内(490~530 nm)最大一阶微分值对应的波长位置蓝谷深度Db蓝边内(490~530 nm)光谱一阶微分的最大值绿峰位置λg510~580 nm波段内最大的波段发射率所在位置绿峰峰值Rg510~580 nm波段内最大的波段发射率光谱面积参数红边面积SDr红边内(680~750 nm)一阶微分值的总和蓝边面积SDb蓝边内(490~530 nm)一阶微分值的总和黄边面积SDy黄边内(550~650 nm)一阶微分值的总和
将获取的61份实测数据分为两部分,一部分是测得冠层光谱并测得对应冬小麦各种重金属含量的41份样本数据,这部分数据进行建模与验证,其中随机选择31个样本数据用于建立反演模型,10个样本为检测样本用于检验模型;另一部分剩余的20份只采集光谱数据而未测重金属含量的样本用来反演。采用模型决定系数(R2)、校正集的均方根误差(RMSEC)、检验集均方根误差(RMSEV)以及相对分析误差(RPD)、相对误差(RE)来检验模型精度,计算公式如式(3)—(5)所示。其中,R2越大,RMSEC越小,建模效果越好;RMSEV、RE越小,RPD越大,模型的预测精度越高。
(3)
(4)
(5)
本研究区内41个采样点的冬小麦重金属含量描述性统计结果见表2,Cr、Ni、Pb、Zn、Hg、Cd含量的平均值分别为8.768、3.663、4.279、35.493、0.011、0.149 mg/kg;变异系数反映数据离散程度,重金属含量变异系数表现为Ni>Cr>Pb>Cd>Hg>Zn,表明本研究区Cr、Ni含量相较Hg、Pb、Zn、Cd含量离散程度大,地区间分布不均。
表2 山东典型污灌区冬小麦叶片重金属含量的描述性统计结果
由图4可知,原始光谱反射率对不同重金属含量的波谱响应不明显,相关性曲线相对平缓,相关系数较低,基本在-0.25~0.25,重金属Cr、Ni、Pb、Zn含量与700~900 nm波段相关系数较大,相关性较强。原始光谱反射率与重金属Hg含量呈负相关,在380~700 nm的可见光范围内的相关性很强,相关系数为-0.45左右,而原始光谱反射率与重金属Cd含量相关系数的高值区在750~1 100 nm,但总体上原始光谱反射率与重金属含量的相关性较低,对重金属含量信息的反映很微弱。
在经过一阶、二阶微分变换后,各波段与重金属含量的相关性明显增强,光谱响应剧烈,相关系数的变化幅度很大,在正负之间来回波动,相关系数的绝对值基本都在0.35以上,重金属Hg含量与变换后各波段的相关性较高,相关系数的绝对值在0.40~0.75。光谱参数与各种重金属含量的相关系数绝对值基本在0.25~0.40(图5),相关性大部分都达到0.05显著性水平,其中重金属Hg含量与黄边面积、红边面积的相关系数分别为0.78、-0.75,相关性较强,达到了极显著水平(P<0.01)。
图4 冬小麦冠层原始光谱反射率、反射率一阶微分以及反射率二阶微分与重金属含量的相关系数
由于光谱波段数目繁多,本研究对原始光谱反射率、反射率一阶微分、反射率二阶微分进行重采样,选出与各种重金属含量相关系数大于0.35且呈显著相关的波段进行建模。建模过程中,以光谱参数及重采样的原始光谱反射率、反射率一阶微分、反射率二阶微分作为模型的自变量,以各种冬小麦叶片重金属含量作为因变量,建模比较结果见表3。从表3可以看出,利用偏最小二乘回归法建立的回归模型的R2都达到0.70以上,建模稳定性较高,建模效果明显优于逐步多元线性回归法,尤其对于光谱参数来说,建立的偏最小二乘回归模型的稳定性远远大于建立的逐步多元线性回归模型,误差也小得多。相比原始光谱反射率建模,利用反射率低阶微分建立的模型的R2总体明显增大,RMSEC较小。对于重金属Cr,SDR-PLSR模型效果最优,R2可达到0.846,RMSEC较小;对于重金属Ni,SP-PLSR模型的R2最高,达0.887,RMSEC最小,为1.313,建模效果很好;对于重金属Pb,FDR-PLSR建模效果最优,R2为0.848,RMSEC较小,为1.964;对于重金属Zn,FDR-PLSR模型的R2达0.790,RMSEC较小,为5.139,建模效果最优;对于重金属Hg,SP-PLSR模型的R2最高,为0.819,RMSEC较小,为0.002,建模效果很好;对于重金属Cd,FDR-PLSR建模效果最优,R2可达到0.868,模型稳定性远远高于其他模型的稳定性,RMSEC为0.085,与其他指标差距较小。
图5 冬小麦冠层光谱参数与重金属含量的相关系数
表3 冬小麦冠层各光谱指标与6种重金属含量的建模结果比较
对回归模型进行检验(表4)发现,对于各种重金属,利用逐步多元线性回归法建立的反演模型的RMSEV较大,RPD较小,在0.509~1.057,RE较大,说明对重金属含量的预测能力一般;而偏最小二乘回归模型的RPD为0.717~2.406,RMSEV、RE较小,说明偏最小二乘回归模型对重金属含量的预测误差较小,预测能力较好。对于重金属Cr,SDR-PLSR模型的RMSEV为1.136,RPD为2.013,RE为26.711%,与其他指标相比,模型预测效果最优;对于重金属Ni,SP-PLSR模型的RMSEV最小,RPD和RE分别为1.872和22.277%,相比其他指标,RE最小,RPD最大,模型预测效果最好;对于重金属Hg,SP-PLSR模型的RPD为1.684,RE为8.182%,模型预测效果最好;而对于重金属Pb、Zn、Cd,FDR-PLSR模型的RMSEV最小,RE也最小,RPD最大,模型预测效果最优。
综合建模与验证结果可以看出,重金属Cr含量的最优反演模型为SDR-PLSR模型,SP-PLSR模型是重金属Ni、Hg含量的最优反演模型,FDR-PLSR模型是重金属Pb、Zn和Cd含量的最优反演模型。
表4 冬小麦冠层各光谱指标与6种重金属含量反演模型的验证结果比较
续表4 冬小麦冠层各光谱指标与6种重金属含量反演模型的验证结果比较
为了更加直观地表现各种重金属的最优模型拟合结果,绘制6种重金属含量实测值和反演模型预测值的散点图,如图6所示。每种重金属含量的实测值与预测值的拟合曲线都在y=x附近,所选取的最优反演模型的拟合效果都较好。
图6 山东污灌区冬小麦叶片重金属含量实测值和最优反演模型预测值的拟合情况
在每种重金属的最优反演模型选出来后,对20个反演样本的重金属含量进行反演并根据反演出的重金属含量以及实测的重金属含量进行普通克里金插值(图7)。由图7可知,冬小麦叶片中Cr含量为2.579~30.237 mg/kg,Ni含量为0.730~15.374 mg/kg,Cr、Ni含量的高值区都位于研究区的东南部,北部及西北部含量较低;Pb、Zn含量分布的高值区主要位于中部及南部地区,北部含量较低;Hg含量为0.002 5~0.021 8 mg/kg,在西北部较低;Cd含量为0.062~0.478 mg/kg,在中部、北部以及西北部较低,其他位置较高。
图7 山东典型污水灌溉区冬小麦叶片重金属的空间分布
经污水灌溉的农作物会受到重金属的污染,重金属在植物叶片内的积累可通过微弱的光谱波动体现出来[32-33]。本研究利用冠层高光谱信息通过建立回归模型对污灌区冬小麦的重金属含量进行估算。结果表明, 对于Pb、Zn、Cd,基于反射率一阶微分的偏最小二乘回归模型(FDR-PLSR)为最优模型(Pb:R2=0.848,RPD=1.598;Zn:R2=0.790,RPD=2.295;Cd:R2=0.868,RPD=2.406);对于Cr,基于反射率二阶微分的偏最小二乘回归模型(SDR-PLSR)为最优模型(R2=0.846,RPD=2.013);对于Ni、Hg,基于光谱参数的偏最小二乘回归模型(SP-PLSR)为最优模型(Ni:R2=0.887,RPD=1.872;Hg:R2=0.819,RPD=1.684)。
从空间插值结果可以看出,冬小麦叶片中Cr、Ni含量在研究区东南部较高,北部及西北部较低;Pb、Zn含量在中部以及南部较高;Hg含量在西北部较低;Cd含量在中部、北部、西北部较低,其他区域较高。研究表明,选取的最优模型能较精确地对污水灌溉区冬小麦叶片重金属进行估算,可实时且大面积地进行冬小麦叶片重金属污染监测。
本研究所用的是野外实测冠层光谱数据,在采样过程中,土壤背景对冠层光谱的影响不可避免,并且由于是区域实地监测,不同区域土壤质地、水分状况等因素差异较大,可能对光谱产生一定影响。在下一步的研究中,将研究土壤质地以及含水量对冠层光谱的影响,并结合高空高光谱影像在更大范围更加精确地进行重金属污染监测。