基于地形校正后Landsat 8的土壤重金属反演研究

王海潇，王勇辉,2

(1.新疆师范大学地理科学与旅游学院，新疆 乌鲁木齐 830054；2.新疆干旱区湖泊环境与资源实验室，新疆 乌鲁木齐 830054)

土壤中的重金属滞留时间较长而且很难降解，因而易在植物体内累积，进而通过整个食物链在人体中富集，造成潜在的健康风险[1-2]。随着社会经济的快速发展，人为活动也成为土壤污染最主要的原因，工业生产、交通运输、矿业开采、农药滥用等均为增加土壤重金属含量的主要途径[3-4]。近几年来，国内外学者对土壤重金属污染开展了广泛研究[5-6]。传统的重金属原位监测技术精度较高，但需耗费大量的人力物力，无法适用于大范围重金属污染监测需求，而遥感技术具有快速、宏观的特性，表现出其他方法无法比拟的优点。如Malley和Williams利用湖泊沉积物的反射光谱快速预测重金属含量[7]；Ren等[8]对矿区农用地中的反射光谱与重金属进行研究，发现反射光谱对As和Cu元素有一定的指示性；吴昀昭等[9]应用农田土壤光谱反射率估算土壤表层Hg元素的浓度；李巨宝等[10]采用偏最小二乘方法，构建基于土壤反射光谱的Se，N，Fe含量预测模型。地表反射率表征地球表面对太阳辐射的反射能力，是揭示地物目标本质的最有用信息[11]，而不同地区的像元由于受太阳、大气和地形等多种因素的影响，接收到的太阳辐射有很大的差别，造成有些影像区处于阴影覆盖下，而另一些却处于过度感光状态[12]，导致了“同物异谱”和“异物同谱”现象的出现，这对地物信息提取造成了不利影响。地形校正可以消除由地形起伏而引起的像元反射率的变化，使像元反射率更好地反映地物光谱特性[13]。目前利用地形校正模型消除地形影响，已成为目前高光谱遥感图像步入应用前的有效处理手段，郭云开等[14]和穆悦等[15]通过地形校正对地表反射率进行校正，并取得较好结果。因此本研究在Landsat 8 OLI数据的基础上引入地形校正，旨在提高研究区重金属反演精度。

玛纳斯湖湿地是干旱区典型的湿地，而土壤重金属严重地危害了湿地的生态安全。为了更好地提高重金属反演精度，本研究以Landsat8 OLI数据为基础，利用DEM数据进行地形校正，并且对校正后的反射率进行倒数、导数、对数等数学变换，从而分别建立8种重金属最优的PLSR预测模型，并对土壤重金属的空间分布及生态风险进行评价，进而更好地利用遥感数据对土壤重金属进行预测。

1 材料与方法

1.1 研究区概况

玛纳斯湖湿地位于欧亚腹地新疆维吾尔自治区昌吉回族自治州境内，其所处位置属于温带大陆性干旱气候区，年平均气温6℃，年降水量100～150 mm，年蒸发量3 500 mm左右，是准噶尔盆底西部的一个构造沉降中心，是玛纳斯河地表径流的尾闾，湖体受玛纳斯河季节性水量补给变化而游移变动。玛纳斯湖是一个咸水湖，湖区位于坐标86.0695°～86.3088°E，44.3741°～44.5196°N，湖的南北跨度有28.62 km，东西跨度21.40 km，湖区面积约550 km2，平均深度约6 m。在气候、地形和水文主要因素作用下，玛纳斯湿地地区形成以湖泊沼泽、冲洪积平原、湖积平原为主的地质地貌类型。

1.2 数据来源及数据处理

于2015年11月10日至15日在玛纳斯湖湿地区域内，选择了3个有代表性的区域，分别是湖泊入湖口(A区)、退化区(B区)和恢复区(C区)，每个区布设9个样点，共布设27个采样点，每个采样点采集0～20 cm土层土壤样本2个，共采集到54个土壤样本。将采集的样品在室内风干，剔除枯枝落叶、根茎、石子、动物残体等杂质，碾碎，用100目尼龙筛过筛后供测试用。土壤样品中As 和 Hg采用原子荧光法(AFS)测定；重金属Cu、Mn、Cr、Pb、Co、Ni、Zn采用离子体光谱(ICP)测定，将每个采样点的测定数据取平均值。为了与野外实测数据相结合，选择的影像数据为野外采样时间段内的Landsat 8 OLI遥感数据，获取时间为2015年11月15日，DEM数据的来源为中国科学院数据中心的ASTER GDEM 数据，分辨率为30 m。

地形校正法是目前应用最为广泛的C校正模型，前期数据处理包括几何精校正和COST模型大气校正，几何精校正误差(RMS)小于0.5，然后用矢量边界进行裁剪。利用DEM数据在ENVI软件中计算得到坡度、坡向数据，同时为了消除大气散射和地表相邻点反射光折射造成的像元亮度值(DN值)和太阳入射角α之间的余弦关系，本研究在二阶校正模型的基础上利用IDL(Interactive data language)语言进行改进，并经过二次校正才能得到校正结果。

1.3 PLSR模型的建立与验证

PLSR集成了主成分分析、典型相关分析和普通多元线性回归3种方法的优点，克服了自变量之间的多重线性相关和样本数量小于波长变量的问题，使构建的模型更稳定，有助于多元数据统计分析[16-17]。

1.4 生态风险预警评估方法

生态风险预警评估源于生态风险评价，更强调对生态系统可能存在风险的警示研究，狭义预警仅指对自然资源或生态风险可能出现的衰竭或危机而建立的报警，而广义预警则涵盖了生态风险的维护及防止危机发展的过程[19-20]。对研究区土壤重金属进行生态风险预警评估采用Rapant等[21]提出的生态风险指数法，属狭义预警，其表达式为:

式中，IERi表示超过临界限量的第i种重金属生态风险指数,CAi表示第i种重金属的实测含量(mg·kg-1),CRi表示第i种重金属的临界限量(mg·kg-1),IER表示待测样品生态风险指数。其中CRi采用新疆土壤背景值[22]。

Rapant等[21]给出了相应生态风险划分标准，参考相关研究提出的生态风险程度描述[23-24]，将生态风险预警指数值与警度进行了概念关联[25]，如表1所示。

表1土壤生态风险预警综合判别标准

Table 1 Standard of soil ecological risk warning

风险等级Risk level风险指数Risk index预警类型Warning type风险程度描述Description of risk level1IER≤0无警No warning生态系统服务功能基本完整,生态环境基本未受干扰,生态系统结构完整,功能性强,系统恢复再生能力强,生态问题不显著,生态灾害少。The ecosystem service function is basically complete; the ecological environment is basically un-disturbed; the ecosystem structure is complete; the function is strong; the system has strong re-covery and regeneration ability; the ecological problem is not significant; and the ecological dis-aster is less.20

2 结果与分析

2.1 反射率的地形校正

玛纳斯湖流域地势南高北低、沟壑相间的特点决定了局部微地形是影响研究区表观反射率、植被、土壤、气候等各方面地域分异的主要控制因素。由图1可知，以B1波段反射率和重金属Cu为例，经过地形校正后的B1波段的反射率值与27个实测土壤表层Cu含量的决定系数高于地形校正前，R2从0.46提高至0.52，拟合结果表明，利用DEM校正后的反射率反映土壤重金属Cu的效果较DEM校正前的好。

2.2 土壤反射率(反射率变换形式)与土壤各重金属含量相关性分析

将进行地形校正后的各波段反射率及其各种数学变换(倒数1/R、对数LogR、一阶微分R′)与各土壤重金属含量进行相关性分析，并且从各种数学变换中筛选出与各土壤重金属含量决定系数最大的波段，其相关性均达到0.01显著水平，结果如表2所示。从表2可知，相对于地形校正后原始反射率(R)而言，其倒数(1/R)、对数LogR与各重金属(除Pb以外)的相关性都有所提高，而一阶微分R′与Cu和Pb含量的相关性有所降低，与其余6种重金属的相关性也有所提高。

2.3 土壤重金属反演模型建立与检验

图1 土壤重金属Cu实测值和地形校正前(a)后(b)B1反射率拟合效果对比Fig.1 Relationships between Cu content and reflectance of B1 before (a) and after (b) topographic correction

变换形式Transformation form波段及决定系数Band and coefficientof determinationCrCuMnNiPbZnHgAsR波段 Band21434124决定系数 R20.410.520.410.240.510.420.430.441/R波段 Band11312243决定系数 R20.450.590.510.550.480.450.470.56LogR波段 Band21431356决定系数 R20.420.570.430.460.350.520.490.51R'波段 Band24132512决定系数 R20.460.470.530.480.460.480.460.47

表3 基于单波段的土壤各重金属预测模型

表4 基于多波段的土壤各重金属PLSR预测模型

2.4 土壤重金属空间分布格局及生态风险预警评估

2.4.1 土壤重金属风险空间分布 利用表4中的各重金属的预测模型分别反演各重金属空间分布，同时根据表1将IER值按其大小分为5个等级以表示不同程度的风险程度，并利用ARCGIS制作了研究区2015年11月的土壤重金属风险等级分布图，如图2所示。

通过图2可知，研究区土壤重金属风险等级总体上呈现从西南方向至东北方向逐渐减弱的趋势，并且研究区土壤重金属污染总体上处于轻警以上级别，生态服务功能已开始退化，生态环境受到一定破坏，生态系统结构有所变化，但尚可维持基本功能，受干扰后易恶化，生态问题显现，应该加强对该地区的重金属污染的治理。

图2 土壤重金属风险等级分布Fig.2 Distribution map of risk grade of soil heavy metals

2.4.2 基于实测数据生态风险预警评估 为了验证土壤重金属遥感反演的精度，通过研究区采集的实测样点实测数据进行验证，并且通过对实测数据进行不同分区(A、B、C三个区)来分析研究区重金属生态风险，基于实测数据的土壤重金属风险评价如表5、表6所示。由表5可知，从单项指数看Cr、Cu、Mn、Ni、Pb、Zn、Hg、As处于无警级的概率分别为92.59%、77.78%、96.30%、3.70%、11.11%、48.15%、3.70%、25.93%，处于预警级的概率分别为7.41%、7.41%、3.70%、96.30%、33.33%、40.74%、3.70%、74.07%，Cu、Pb、Zn和Hg处于轻警级的概率分别为11.11%、55.56%、3.70%、3.70%，其他元素均未出现轻警级的样点，处于中警级的重金属元素主要有Cu和Hg，其概率分别为3.70%、37.04%，处于重警级的重金属元素主要有Zn和Hg，其概率分别为4.41%、51.85%。其单项生态风险指数处于无警以上级别的风险概率大小顺序为Hg=Ni>Pb>As>Zn>Cu>Cr>Mn。从综合指数看研究区11.11%的样点处在无警级别，属于最低生态风险，有33.33%的样点处于预警级风险，而样点落在轻警、中警、重警级风险的概率分别为25.93%、18.52%、11.11%，由表5、表6可知，属于重警级污染的样点主要分布在恢复区(C区)，并且研究区土壤重金属生态风险排序为：恢复区(C区)>退化区(B区)>湖泊入湖口(A区)。

通过比较基于遥感的土壤重金属生态风险预警与基于实测值的土壤重金属生态风险预警可知，两者的分析结果较为一致，结果表明可以用遥感的手段来反演该研究区的重金属分布情况。

表5 研究区重金属生态风险分级评价

表6 研究区不同分区重金属生态风险分级评价

3 结论与讨论

目前已有不少学者利用遥感技术对土壤重金属污染监测进行研究探索，并取得一定的成果[26]，但是很少有研究会考虑去除地形对反射率的影响。本研究在Landsat 8的基础上，引入DEM数据对反射率进行地形校正，经过地形校正后的B1波段的反射率值与27个实测土壤表层Cu含量的决定系数高于地形校正前，R2从0.46提高至0.52，拟合结果表明，利用DEM校正后的反射率反映土壤重金属Cu的效果较DEM校正前的好。诸多学者的研究表明，将反射率进行适当的数学变换可以有效抑制遥感影像中噪声对地物光谱的影响，改善反射率与土壤某些元素的相关性，提高模型反演精度[27]，因此本研究在地形校正后原始反射率的基础上进行倒数、对数和导数等变换，研究发现相对于地形校正后原始反射率(R)而言，其倒数(1/R)、对数LogR与各重金属(除Pb以外)的相关性都有所提高，这与郭云开等[28]的研究结果较为一致。

在重金属的反演模型中，基于多波段的PLSR预测模型整体上相对于单波段要好，同时利用土壤各重金属的最佳预测模型分别反演相应的土壤重金属含量，并引入土壤重金属生态风险指数用于评价研究区的土壤重金属生态风险并进行预警，研究发现研究区土壤重金属风险等级总体上呈现从西南方向至东北方向逐渐减弱的趋势，土壤重金属生态风险排序为恢复区(C区)>退化区(B区)>湖泊入湖口(A区)，这与实地验证情况较为一致，原因可能为C区为恢复区，靠近克拉玛依市，其受到克拉玛依市生活废水、工业生产废水、农业污水灌溉的影响，使得该区的土壤重金属风险高于A区和B区，而A区由于为湖泊如湖口，受到湖水对土壤的冲刷作用，使得土壤重金属风险最低。研究区土壤重金属主要来源可能是来自工业生产、生活用水中的重金属成分，这与刘娜等[29]、马辉英等[30]的研究结论基本一致。

为了验证基于遥感的土壤重金属生态风险预警的预测精度，将研究区土壤重金属含量实测数据也通过重金属生态风险指数进行计算，通过比较发现两者的分析结果较为一致，结果表明可以用遥感的手段来反演该研究区的重金属分布情况，同时研究区土壤重金属污染总体上处于轻警以上级别，生态服务功能已开始退化，应该加强该地区的重金属污染治理工作。