基于HJ1A-HSI的龙口污水灌溉区土壤重金属含量反演

蔡东全，吴泉源，曹学江，蒙永辉，周历媛，刘 硕，张龙龙

（1.山东师范大学人口·资源与环境学院，山东　济南250014；2.山东省煤田地质规划勘察研究院，山东泰安271000；3.山东省地质环境监测总站，山东济南250014）

基于HJ1A-HSI的龙口污水灌溉区土壤重金属含量反演

蔡东全1，吴泉源1，曹学江2，蒙永辉3，周历媛1，刘 硕1，张龙龙1

（1.山东师范大学人口·资源与环境学院，山东济南250014；2.山东省煤田地质规划勘察研究院，山东泰安271000；3.山东省地质环境监测总站，山东济南250014）

污水灌溉在一定程度上减轻了农业用水问题，但是也造成了重金属元素的积累。以HSI高光谱遥感影像为数据源，利用地理信息系统和偏最小二乘回归（PLSR）分析方法，并结合实地采样检测的土壤重金属含量数据，分析了重金属离子在HSI数据中的光谱特征，建立了土壤重金属含量与高光谱数据的偏最小二乘回归模型，实现了对龙口污水灌溉区土壤重金属含量的反演研究。结果表明:原始光谱数据经过一定的光谱变换后，就建模效果而言，Cu的建模效果最好，判定系数（R2）高达0.965 5，Mn、Ni和Pb的建模效果较好，R2分别为0.927 3、0.910 4 和0.912 7；验证样本中，Mn的R2为0.9128，说明该模型对Mn的反演能力最强，As、Cu和Ni的模型反演能力相当，R2分别为0.727 6、0.780 4和0.733 3。基于HJ1A-HSI数据并结合PSLR方法反演龙口污水灌溉区土壤重金属含量具有方法可行性与结果的可靠性。

土壤重金属；HSI；偏最小二乘回归；污水灌溉区

合理的污水灌溉不仅可以缓解一个地区农业用水的紧张局面，而且可以充分利用污水中所含有的作物营养元素，降低农业生产成本，还可以提高水资源的重复利用率，并可作为污水处理回用的一种有效手段［1］。然而随着灌溉时间的不断延长和灌溉面积的不断扩大，这种方式带来了十分显著的土壤重金属积累问题。受到重金属污染的土壤常常通过农作物和水进入人体，经过富集作用对人类和环境产生危害［2］。

传统的土壤重金属监测方法可以针对不同的重金属元素种类取得较好的测量精度，但是只能对特定的采样点或剖面进行监测，无法做到实时、动态地对重金属污染进行监测研究，且效率低下、耗时费力［3-4］。高光谱遥感数据具有更精细的地物光谱曲线，更强的目标识别和信息提取能力［5］，使得以遥感手段定量反演和研究土壤的化学成分成为可能，也为定量反演土壤重金属含量提供了一条新的思路。近年来，国内外的一些学者尝试利用高光谱数据，对土壤重金属含量进行定性与定量研究［6-8］，但以我国环境减灾卫星HSI高光谱为数据源，定量反演土壤重金属含量的研究则相对较少。因此，本文利用地理信息系统和偏最小二乘回归方法，研究龙口污水灌溉地区土壤几种主要重金属盐分参数与HSI高光谱影像的光谱反射率之间的关系，为快速、动态和大面积定量反演污水灌溉区土壤重金属含量寻求一种新的方法。

1　研究区概况

龙口市地处胶东半岛的西北部，地理坐标为东经120°13′20″～120°44′45″，北纬37°27′45″～37°47′30″。全市地势南高北低，呈阶梯形下降，东部和南部为低山丘陵地带，西部和北部为滨海平原。受成土母质、地形、地貌、气候和水文条件的影响，该地区形成了棕壤、褐土、潮土和姜黑土四大土壤类型。研究区内有效灌溉面积超过8 000 hm2，其中污水灌溉总面积约4 500 hm2，主要灌溉水源为龙口、黄城两座污水处理厂排放的工业废水和生活废水。

2　数据获取与预处理

2.1土壤数据的获取

本研究以Landsat8数据7、4、1波段合成的遥感图像为参考图，并利用GPS定位技术于2013年11月5—15日在研究区内进行野外调查和采点，选择32个100 m×100 m样方区（见图1），并要求样方区在各种土地利用类型中尽可能规则分布，且样方区周围的土壤性质异质性较小、环境因子类似，便于分析统计。在每个样方区内采用网格法采集表层土壤样品10～15个，采样深度为20 cm，并用数码相机对采样区进行拍照，实地记录采样点土壤表面状况及景观特征。将采集的土壤样本在实验室内自然晾干，研磨后过1 mm筛，剔出土壤中的石块等杂质，按水土比5∶1配置土壤浸提液，在样品中加入热硝酸溶解后并采用电感耦合等离子体原子发射光谱法（ICP-AES）测定重金属元素含量，最后得到包括铬（Cr）、锰（Mn）、铁（Ee）、镍（Ni）、铜（Cu）、锌（Zn）、砷（Sb）、镉（Cd）、铅（Pb）9种元素在内的42种土壤盐分离子含量，去掉异常值后，每个样方区采样点重金属元素含量取平均值，即得到所需的土壤重金属含量数据。

图1　采样点示意图Eig.1 Sketch diagram of the sampling location

2.2土壤样本重金属含量分析

本研究对32个土壤样本的重金属含量进行了简单的统计分析，其结果见表1。

由表1可以看出:研究区土壤样本中重金属元素含量的变异系数普遍较小，其中Cu的变异系数最大，为9.79%；从偏度上看，除Zn以外，其他重金属元素都属于中高度偏态分布；另外，所有重金属元素含量的最大值和最小值之间相差较大，说明研究区内重金属元素的分布不均衡。

2.3高光谱遥感数据的获取

遥感数据来源于我国环境小卫星（HJ-1 A）的HSI高光谱数据。我国环境小卫星在2009年3月30日开始正式交付使用，分为A星和B星，简写为HJ-1A/1B。两颗星都配备了CCD相机，其中HJ-1A星搭载的超光谱成像仪（HSI）能够完成对地展宽为50 km、地面像元分辨率为100 m、110～128个光谱谱段的推扫成像，具有±30°侧视功能和星上定标功能，重访周期为4 d，能够迅速、大尺度地观测大气、土壤和植被等状况，可为反演计算提供数据支持［9-10］。

表1　土壤样本重金属含量分析结果Table 1 Contents of the heavy metals in soil-samples

本研究高光谱遥感数据源选择2013年10月覆盖整个研究区的3景HJ-1A高光谱2级产品，要求所选取的影像清晰，无云或少云覆盖，且与采集土壤样本的时间尽量接近，用以提高土壤重金属含量的反演精度。

2.4HSI高光谱数据的预处理

HSI高光谱数据在应用于重金属含量反演之前，需要进行包括几何纠正（误差控制在半个像元内）、大气校正、投影变化、拼接和裁剪等在内的几何及光谱方面的预处理。本数据源属于2级产品，在下载之前已经完成了波谱复原和系统粗校正等处理。产品在生成时选用了100的扩大因子，所以需要将像元DN值还原为标准辐射亮度，即

式中:Lλ为真实辐射亮度值［W/（m2·sr·μm）］；DN为原始图像像元值［W/（m2·sr·μm）］。

通过对不同地物辐射亮度光谱曲线的观察，发现光谱曲线的总体形态与趋势相似，且受大气中水汽的影响，不同地物在760 nm和820 nm附近均存在两个明显的波谷，如果直接利用未经校正的影像，在水汽的影响下将很难体现地物的真实物理特性，从而影响遥感反演的精确度。

ENVI5.1能够充分满足遥感影像在读取、配准和大气校正等方面的需求，因此本研究利用该软件提供的ELAASH模型，并通过读取影像头文件获取的模型输入参数对HJ-1A数据进行大气校正，最终得到115个波段的地表反射率图像，此时的光谱曲线可以较好地体现出真实地物的光谱特征。HSI高光谱数据覆盖450～950 nm波段，而在480 nm波段之前，即前10个波段的噪声非常大，信息的准确性不高，因此信息量主要集中在480～950 nm波段之间［11］。

3　研究过程

3.1模型方法的选择

高光谱遥感拥有较强的连续性波段、丰富的光谱信息和高分辨率等优点，它在特定光谱域内蕴含着近似连续的地物光谱信息，使地物的识别能力极大地提高［12-13］。然而由于高光谱数据的波段较多、信息冗余现象严重，使得各波段间存在着显著的相关性，而利用传统的最小二乘回归方法分析高光谱数据建立的估计模型时，无法避免地会遇到多重共线性和自由度过小的问题，将造成模型方程系数矩阵的严重病态性，使模型参数的最小二乘估计失真。

偏最小二乘回归（PLSR）方法是在传统的最小二乘回归分析方法的基础上，引入主成分分析和方差分析的思想，为一种多变量回归分析方法［14］。基于该方法，利用高光谱数据对土壤重金属含量进行反演计算，其主要思路是从高光谱原始数据及其变换形式中提取能够最大解释土壤盐分含量变化的主成分，提取的主成分是各波段光谱数据的线性组合，各个主成分之间是正交关系，既最大程度地保存了信息量，又消减了这些主成分线性相关性的影响。

PLSR方法并不直接考虑因变量对自变量的回归建模，而是从自变量中逐步提取多个对自变量系统X和因变量系统Y都具有最佳解释能力的新综合变量t1，亦称之为主成分。首先建立yk对主成分的MLR回归方程；然后还原为yk关于原自变量系统的PLS回归方程，PLS关键是提取主成分；最后利用模型的交叉有效性Q2h判定回归方程是否达到满意的精度。为了判定模型的精度，引入判定系数（R2）、均方根误差（RMSE）和相对均方根误差［RMSE（%）］3个判定指标，R2越大、RMSE和RMSE越小，表明该模型稳定性和精度越好。均方根误差（RMSE）亦称标准误差，其计算公式为

相对均方根误差（RMSE）计算公式为

上式中:εi为土壤样本重金属元素含量实测值（mg/kg）；ˆεi为土壤样本重金属元素含量PLSR模型预测值（mg/kg）；¯ε为土壤样本重金属元素含量均值；n为样本个数。

本文从32个采样点中随机抽取20个作为建模样本，余下12个作为检验样本，采用SIMCA-P软件进行统计分析。

3.2光谱指标的提取

为了提高反演效果，将从遥感影像获取的原始光谱反射率进行一系列的光谱变换，即以HSI影像的光谱反射率及多种数学变换形式为特征向量，提取土壤光谱信息。这些数学变换形式包括光谱去除包络线（CR）、光谱反射率倒数的对数变换［Log（1/ R）］、光谱反射率的一阶微分（RD1）和二阶微分（RD2）变换。

光谱去除包络线是一种常用的光谱处理方法，它是将反射率归一化在0～1.0范围之内，并将其归一化到统一的背景值上，突出了光谱的反射和吸收特征［15-16］。将光谱反射率进行倒数的对数变换（log （1/R））是处理非线性问题的方法，它不仅趋向于增强可见光区的光谱差异，而且趋向于减少因光照条件变化引起的乘性因素的影响［17］。对反射率值进行微分变换有助于限制低频噪声对目标光谱的影响［18-19］。对于光谱反射率的微分计算，一般用光谱的差分作为微分的有限近似，光谱数据的一阶微分计算公式近似为

式中:λi为波长（nm）；Δλ为相邻的两个波长λi+1和λi的间隔（nm）；ρ′（λi）为波长λi的一阶微分。

3.3模型分析

为了研究各种重金属离子在HSI影像中的光谱响应，将采样点重金属离子含量与影像的反射率及光谱变换形式建立偏最小二乘回归模型，其建模结果见表2。

由表2可以看出:建模所用的光谱指标全部是光谱反射率的数学变换形式，包括微分变换、去除包络线和反射率倒数的对数变换等；由于在前期对原始数据采用9点加权移动平均法进行平滑处理，所以选择的主成分个数较少；对于建模效果而言，Cu的建模效果最好，判定系数（R2）高达0.965 5，Mn、Ni和Pb的建模效果较好，R2分别为0.927 3、0.910 4和0.912 7，As和Zn的R2分别为0.855 5 和0.832 7，Cr的建模效果较差，R2为0.640 5；验证样本中，Mn的R2最大，为0.912 8，说明该模型对Mn的预测能力最强，As、Cu和Ni的模型预测能力相当，R2分别为0.727 6、0.780 4和0.733 3，Cr和Zn的模型预测能力一般，R2分别为0.598 9 和0.622，Pb的模型预测能力较差，R2为0.515。可见，利用光谱反射率及其变换形式估算土壤重金属离子含量大部分能够获得较大的判定系数，取得了较好的反演效果，说明HSI影像对土壤重金属离子具有较强的探测能力。

表2　PLSR建模结果比较Table 2 Comparison of the PLSR modeling results

图2为利用PLSR方法建立的土壤重金属含量预测值与实测值比较的效果图。实测值与预测值越均匀分布于1∶1直线两侧，说明基于偏最小二乘回归方法建立的土壤重金属含量估测模型反演效果越好。

通过对建模结果的比较，选择相应的敏感波段，获得了土壤重金属元素含量的回归方程（见表3），其中Xn表示重金属元素含量第n个波段经过特定光谱变换形式后的数值。

图2　土壤重金属含量实测值与预测值的比较Eig.2 Comparation of measured and predicted values of heavy metal ions

表3　反射率不同变化形式与土壤重金属含量回归方程Table 3 Regression equations between different variations of the reflectivity and the soil heavy metal content

利用土壤各重金属元素含量的回归模型，并结合HSI遥感影像反演了研究区范围内土壤重金属含量分布，其结果见图3。由图3可以看出:龙口地区土壤的重金属污染主要集中在西北部沿海地区，随着污水灌溉年限的延长，污水中的重金属离子逐渐在龙口地势低洼的地区积累起来。

图3　研究区土壤重金属含量分布图Eig.3 Distribution map of the contents of soil heavy metal ions

4　结 论

本文以HSI高光谱遥感影像为数据源，利用地理信息系统和偏最小二乘回归（PLSR）分析方法，并结合实地采样检测的重金属含量数据，分析了土壤重金属离子在HSI数据中的光谱特征，并建立土壤重金属含量与高光谱数据的偏最小二乘回归模型，实现了对龙口地区土壤重金属含量的反演研究，得到如下结论:

（1）利用HSI原始高光谱数据及其变换形式反演土壤重金属含量，能够获得较好的建模效果，从而证明了环境小卫星HJ-1A影像对土壤重金属离子具有较好的响应能力，其空间覆盖能力可满足将来大面积的土壤重金属含量反演，具有一定的可行性。

（2）偏最小二乘回归（PLSR）方法能够有效地解决样本中多重共线性和自由度过小等问题，在非正态分布、小样本的情况下，仍然很稳定。

（3）从建模效果较好的重金属出发，对龙口污水灌溉地区的土壤重金属含量进行反演计算，初步尝试了由点到面的研究，研究结果表明龙口污水灌溉地区的土壤重金属污染主要集中在西北部沿海地区。

但是由于本次研究样本数据较少，HSI高光谱数据空间分辨率较低，仅为100 m×100 m，在混合像元的影响下必然降低了反演精度，因此在进一步的研究中应该着重考虑光谱重建的问题。

［1］张淑丽.浅析污水灌溉环境污染的研究进展［J］.科技情报开发与经济，2009，19（27）:128-129.

［2］陈建勇，张江山，郑育毅.基于熵权的物元分析法在土壤重金属污染评价中的应用［J］.安全与环境工程，2011（5）:57-60.

［3］龚绍琦，王鑫，沈润平，等.滨海盐土重金属含量高光谱遥感研究［J］.遥感技术与应用，2010，25（2）:169-177.

［4］徐明星，吴绍华，周生路，等.重金属含量的高光谱建模反演:考古土壤中的应用［J］.红外与毫米波学报，2011，30（2）:109-114.

［5］刘焕军，张柏，王宗明，等.基于反射光谱特征的土壤盐碱化评价［J］.红外与毫米波学报，2008，27（2）:138-142.

［6］刘华，张利权.崇明东滩盐沼土壤重金属含量的高光谱估算模型［J］.生态学报，2007，27（8），3427-3434.

［7］王维，沈润平，吉曹翔.基于高光谱的土壤重金属铜的反演研究［J］.遥感技术与应用，2011，26（3），348-354.

［8］Kooistra L，Wehrens R，Leuven R S E W，et al.Possibilities of visible-near-infrared spectroscopy for the assessment of soil contamination in river flood plains［J］.Analytica Chimica Acta，2001，446（1/2）:97-105.

［9］李传荣，贾媛媛，胡坚，等.HJ-1光学卫星遥感应用前景分析［J］.国土资源遥感，2008，9（3）:1-9.

［10］中国资源卫星应用中心.HJ-1A、B卫星介绍［EB/OL］.http://www.cresda.com/n16/n1130/n1582/8384.html.

［11］马驰.基于HJ-1A反演松嫩平原土壤盐分含量［J］.干旱区研究，2014，31（2）:226-230.

［12］浦瑞亮，宫鹏.高光谱遥感及其应用［M］.北京:高等教育出版社，2000.

［13］宫鹏，黎夏，徐冰.高分辨率影像解译理论与应用方法中的一些研究问题［J］.遥感学报，2006，10（1）:1-5.

［14］李世玲.基于PLS成分的变量筛选法［J］.信息与电子工程，2003，1（2）:31-35.

［15］何挺，王静，林宗坚，等.土壤有机质光谱特征研究［J］.武汉大学学报，2006，31（11）:975-979.

［16］张雪红，田庆久.基于连续统去除法的冬小麦叶片氮积累量的高光谱评价［J］.生态学杂志，2010，29（1）:181-186.

［17］Clark R N.Spectral properties of mixtures of montmorillonite and dark carbon grains:Implications for remote sensing minerals containing chemically and physically adsorbed water［J］.Journal of Geophysical Research，1983，88:10635-10644.

［18］胡军杰，陈美兰，吴连喜，等.铀矿区植物的光谱特性分析［J］.东华理工大学学报（自然科学版），2011，34（4）:374-400.

［19］Huang Z，Turner B J，Stephen J D，et al.Estimating foliage nitrogen concentration from HYMAP data using continuum removal analysis［J］.Remote Sensing of Environment，2004，93: 18-29.

Quantitative Retrieval of Soil Heavy Metal Content in Longkou Wastewater Irrigation Area Based on HJ1A-HSI Images

CAI Dongquan1，WU Quanyuan1，CAO Xuejiang2，MENG Yonghui3，ZHOU Liyuan1，LIU Shuo1，ZHANG Longlong1
（1.College of Population，Resources and Environment，Shandong Normal University，Jinan 250014，China；2.Shandong Provincial Coal Geological Planning，Exploration and Research Institute，Taian 271000，China；3.Shandong Monitoring Center of Geological Environment，Jinan 250014，China）

Wastewater irrigation can alleviate water problems in agriculture to some extent，but it causes the accumulation of heavy metals.This paper adopts the HSI hyperpectral remote sensing images as the data source using geographic information system and PLSR analysis method.Combining with concentration data of heavy metals through field sampling and monitoring，this paper analyzes spectral characteristics of heavy metal ions in HSI data，and establishes the equation between concentration data of heavy metals and hyperspectral data to achieve the quantitative retrieval of soil heavy metals contents in Longkou wastewater irrigation area.After spectral transformations of the original spectral data，Cu element is optimal in terms of modeling performance，with the determination coefficient（R2）up to 0.965 5；Mn，Ni and Pb are suboptimal with the R2values of 0.927 3，0.910 4 and 0.912 7，respectively.In validation samples，the R2of Mn is 0.912 8，indicating the best retrieval performance.And then it is followed by As，Cu and Ni，with the R2values of 0.727 6，0.780 4 and 0.733 3，respectively.Therefore，it is feasibleand reliable to retrieve soil heavy metal contents in Longkou wastewater irrigation area by using the data of HJ1AHSI and method of PSLR.

soil heavy metals；HSI；partial least squares regression（PLSR）；wastewater irrigation

X53

A

10.13578/j.cnki.issn.1671-1556.2015.05.006

1671-1556（2015）05-0033-07

2015-03-30

2015-08-11

国家自然科学基金项目（41371395）；国家科技支撑计划项目（2012BAB11B01）；龙口矿区及周边海岸带遥感监测研究项目（鲁勘字［2012］110号）；海（咸）水入侵典型区数值模型与预警预报系统开发建设项目（SDZS-2015-DHJ01）

蔡东全（1989—），男，硕士研究生，主要研究方向为遥感地学应用。E-mail:cdq8906@126.com

吴泉源（1959—），男，博士，教授，主要从事遥感地学应用方面的研究。E-mail:wqy642052582@163.com