基于GWRK的富硒区表层土壤镉含量预测及农业产业规划研究

项剑桥　戴云哲

摘要：構建回归克里格模型（RK）和地理加权回归克里格（GWRK）模型，以土壤硒含量为基础，构建对土壤镉含量进行预测的模型，探索在只具备土壤硒含量数据的条件下进行土壤镉含量预测的方法，以及镉污染风险评估防治和富硒农业规划。结果表明，①柴湖镇富硒耕地面积达到5 316.57 hm2，具备发展富硒农业的条件；②耕地土壤镉背景值处于《土壤环境质量标准》Ⅱ类中含量偏低的水平，基本不会对植物和环境造成危害和污染；③土壤硒含量与镉含量之间存在显著的正相关关系；④利用GWRK预测的土壤镉含量与检测值较为接近，精度明显高于RK的预测结果，可作为区域土壤镉含量评估的依据；⑤RK和GWRK模型的预测结果都存在值域向中值收缩的现象，是造成模拟精度下降的主要原因，而在中值附近预测的精度明显更高。GWRK模型可以用于在已知土壤硒含量的情况下镉含量的预测，根据不同区域耕地土壤硒含量和镉含量的差异，进行差别化的农业规划。

关键词：地理加权回归克里格（GWRK）模型；富硒；镉；农业规划

中图分类号：S153.6+1 文献标识码：A 文章编号：0439-8114（2017）10-1832-07

DOI：10.14088/j.cnki.issn0439-8114.2017.10.009

Prediction for Cadmium in Surface Soil and Agricultural Industrial Planning in Selenium Enriched Area Based on GWRK

XIANG Jian-qiao1， DAI Yun-zhe2

（1.Hubei Selenium Industrial Research Institute，Wuhan 430034，China；

2.College of Public Administration，China University of Geosciences （Wuhan）， Wuhan 430074，China）

Abstract： Geographically weighted regression Kriging（GWRK） model and regression Kriging（RK） model were employed. On the basis of soil selenium content， the model to predict the soil cadmium content was constructed， which was explored the soil cadmium content prediction method under the condition of only having soil selenium content data. And cadmium pollution prevention and control of risk assessment and selenium-rich agricultural planning were explored. The results indicated that ①Chaihu is qualified for establishing selenium industry as the area of selenium-rich cultivated land reached 5 316.57 hm2. ②The cadmium background value in soil fits the class Ⅱ standard of Soil Environmental Quality Standard of China，which is safe for plants and the environment. ③There exists significant positive correlation between selenium content and cadmium content in soil. ④GWRK model is qualified for predicting cadmium content in soil based on the data of selenium content in soil， and which is much better than RK does. ⑤Both RK and GWRK model have a similar issue of value range shrinkage， which leads the decline in simulation accuracy， and the accuracy of prediction is significantly higher in the vicinity of the median. Overall， GWRK model could be used for predicting cadmium content based on selenium content， and different plating plans should be carried out for cultivated land according to the differentiation of selenium content and cadmium content in soil.

Key words： geographically weighted regression Kriging（GWRK） model； selenium enrichment； cadmium； agricultural industrial planning

硒是人体和动物必需的微量元素之一，而中国72%的地区属于缺硒或低硒地区，2/3的人口存在不同程度的硒摄入不足，普遍存在潜在的健康危机[1]。富硒食品具有提高人体免疫力、预防各种疾病的作用，而富硒土壤是开发富硒农产品的天然载体[2]，因此，科学合理利用富硒土壤资源，是提高富硒农产品品质的主要途径。镉是一种非必需且生物毒性极强的重金属元素[3]，土壤中自然存在的镉一般不会对人类造成危害，但长期大量摄入镉会影响钙和磷的代谢，诱发骨质疏松等疾病[4]。

富硒土壤往往伴生重金属[5]，而硒对金属元素具有较强的亲和能力，可能会影响到重金属在土壤-农作物系统中的分布和富集转移性[6]，富硒区农业规划中应当高度关注重金属污染的潜在风险。土壤样品多因子的检测工作会大量增加财力、人力和时间上的消耗。在完成土样硒含量的测定之后，如果能够利用硒含量的检测结果直接对土壤镉含量进行达到一定精度的预测，将能在不增加财力、人力和时间成本的情况下完成富硒地区镉含量评估以及响应农业规划策略的研究工作。因此，本研究基于回归克里格（RK）和地理加权回归克里格（GWRK）两种方法，分别以土壤硒含量为基础，对土壤镉含量进行预测模型的构建，探索在只具备土壤硒含量数据的条件下进行土壤镉含量预测的方法，以及镉污染风险评估防治和富硒农业规划等。

1 研究区概况与数据来源

1.1 研究区概况

柴湖镇位于湖北省钟祥市南20 km，地处江汉平原，汉江由北绕西南环抱柴湖。研究区位于汉水以西主要的生产生活区域，面积为18 575.94 hm2。柴湖镇地质分布属江汉平原一级沉降区晚近期构造带。由于受汉江上游的砂岩、页岩、片麻岩等风化物的影响，地层上部有灰青、灰黑色的亚黏土；下部为砂和砂砾石层，并夹有淤泥质粉细沙和淤泥层。地貌为北高南低，海拔高程多在38～46 m。柴湖鎮气候为北亚热带季风气候，四季分明，寒暖干湿较为明显，雨量在时空分布上变差较大，并从南向北递增。年平均气温16 ℃，年平均降雨量811.2 mm。根据柴湖镇第二次全国土地调查变更成果（2013）统计得出，区内耕地面积为12 284.13 hm2，占全区总面积的66.13%，其中旱地11 505.19 hm2，水浇地6.83 hm2，水田772.11 hm2。

1.2 数据来源与处理

依托湖北省“金土地”工程——高标准基本农田地球化学调查项目（JTD20140101），于2014年5月在柴湖镇采集表层土壤土样共计1 244份。在釆样的同时，以1980西安坐标系（西安大地原点）作为记录样点的实际大地坐标。土壤样品采集完成后，送至土壤重构实验室，由专业检测人员对土壤样品进行测试分析，得出土壤中各元素的化验数据，作为土壤硒含量和镉含量评价及后续农业规划研究的基础数据。

对土壤样品数据进行测试统计，结果见表1。将有效土样的土壤硒含量和土壤镉含量与1990年全国土壤背景值[7]进行对比，硒含量大于0.4 mg/kg的样品数量为616个，最大值为0.600 mg/kg，最小值为0.170 mg/kg，平均值为0.389 mg/kg，明显高于全国背景值，变异系数为19.36%，属于低变异。镉含量大于0.4 mg/kg的样品平均值为0.405 mg/kg，变异系数为19.48%，属于低变异。此外，硒和镉样品的峰度分别为2.742和3.172，说明观察量较为集中；两者的偏度都较低，因此数据在后续的计算和应用中无需进行指数变换或Box-Cox变换。

依据GB15618-1995[8]，对柴湖镇1 244个土样的镉含量进行测试，从表1可以看出，土壤镉含量总体比1990年全国背景值略高，处在《土壤环境质量标准》二级标准的范围内，可作为一般农田、蔬菜地、茶园、果园、牧场等土壤，基本上对植物和环境不造成危害和污染。

2 方法

2.1 回归克里格（RK）

当目标变量与辅助变量有相关关系时，先建立目标变量与环境变量的多元（或一元）回归关系，并根据无数据位置0处环境变量的值和回归关系推算x0处目标变量的确定性趋势项m（x0）和残差e（x0）；最后由式（1）得到x0处的估计值z（x0），计算公式如式（1）所示。

z（x0）=m（x0）+e（x0） （1）

2.2 地理加权回归克里格（GWRK）

回归克里格（RK）使用的线性回归方法是基于普通最小二乘法（Ordinary least squares，OLS）的全局回归技术，没有考虑到回归关系的空间非平稳性。而GWR是一种空间局部回归方法，可以用来探测和建模空间关系的非平稳性[9]。地理加权回归克里格（GWRK）则是用GWR的局部回归部分代替回归克里格中的全局回归部分，而残差仍然采用克里格插值的一种混合插值方法[10]，它将确定性和不确定性的两个部分分开进行模拟，并且相较于回归克里格方法更加准确高效[11]。它在普通线性回归模型中增加了一个权函数。该函数仅考虑在特定空间邻域范围内的样点，待测点离样点越近将被赋予更高的权重，反之权重较小[12]。GWR模型为：

式中，yi是因变量，xij是因子变量，？着i为误差项；（μi，vi）为第i个样点的坐标；？茁k（μi，vi）是第i个样点上的第k个回归系数，是地理空间位置函数。

2.3 半方差函数分析

半方差函数是地统计学解释土壤理化性质空间变异结构的重要理论基础，包含3个主要参数：块金值（Nugget）、变程（Major range）和基台值（Sill），其中变程反映土壤性质的空间变异特性，在变程范围外土壤性质表现为空间不相关，在变程范围内表现为空间相关。块金值是由测量误差和最小取样间距内土壤性质的随机性因子（耕作条件、施肥剂量与频率、作物布局等）引起。基台值反映因子区域化变量受结构性因子（地形、气候、成土母质、植被等）影响的程度。土壤性质的空间相关性可根据块金值与基台值之比C0/C进行划分[13]，C0/C值越低，则说明因子的空间相关性越强烈[14]：①当C0/C≤25%时，因子在其最大变程范围内具有十分强烈的空间相关性，结构性因素是主导空间变异大小程度的主要因素；②当25%

克里格插值采用的半方差函数相关参数如表2所示，插值的模型全部采用球面模型，硒检测值的C0/C为47.46%，具有较强的空间相关性；镉的检测值的C0/C为60.66%，空间相关性不强，但仍能反映出一定的空间相关性。RK和GWRK方法下镉的C0/C分别为47.37%和36.84%，具有较强的空间相关性。

2.4 模型拟合精度评价

对模型预测精度的评价一般有R2、MEE、MAEE、RMSE、RPD等指标[10，16，17]，计算公式见式（3）至式（7）。

从表3中MEE、MAEE和RMSE这几个指标来看，RK和GWRK的预测结果都有一定的偏差。就MEE而言，RK小于GWRK，说明检测值在RK拟合模型的两侧分布更为均匀，而在GWRK模型的一侧分布更多。GWRK的MAEE和RMSE都比RK小，表明检测值更加接近GWRK模型。RK模型的R2为0.56，在模型质量中等的范畴中精度相对较差，而GWRK模型的R2为0.76，虽然质量也是中等水平，但明显高于RK模型，已经接近高质量模型R2≥0.80的标准。另外，RK模型的RPD低于1.4，拟合精度较差，GWRK模型的RPD=1.717，处于中高水平。综上所述，GWRK方法明显优于RK方法，且预测的结果是有效而可靠的。

3 结果与分析

3.1 土壤硒分布情况及评价

李家熙等[18]将土壤硒含量分为缺硒土壤（<0.1 mg/kg）、低硒土壤（0.1～0.2 mg/kg）、中硒土壤（0.2～0.4 mg/kg）、富硒土壤（0.4～3.0 mg/kg）和硒过剩土壤 （>3.0 mg/kg），并以土壤硒含量低于0.05 mg/kg为临界值，表示當植物性饲料硒含量低于0.05 mg/kg时，可引发动物缺硒的相关症状。依据此标准对柴湖镇土壤富硒研究的指标同样分为5级。土壤硒含量分布如图1所示，统计结果见表4。由表4可知，硒含量为0.4～3.0 mg/kg的比例达到样品总量的49.60%；含量在0.2～0.4 mg/kg之间的样品所占比例最大，占样品总量的50.40%。表明柴湖镇土壤硒含量以中硒为主，同时具有规模可观的富硒土壤。进一步统计不同硒含量的土地面积可得出，富硒区总面积达到6 869.90 hm2，占研究区土地面积的36.98%，其中耕地面积为5 316.57 hm2。主要分布在柴湖镇的中北部和南部。柴湖镇富硒区土壤含硒量平均为0.45 mg/kg，与国内其他富硒地区相比，属于足硒或中等富硒的水平[1]，具备富硒农业产业园区建设的条件。

3.2 土壤镉分布预测及评价

参照国家环境保护局1995年颁布的国家标准《土壤环境质量标准》，将土壤镉含量分为Ⅰ类（0～0.2 mg/kg）、Ⅱ类（0.2～1.0 mg/kg）、Ⅲ类（1.0～1.5 mg/kg）三个类别[8]。Ⅰ类主要适用于国家规定的自然保护区（原有背景重金属含量高的除外）、集中式生活饮用水源地、茶园、牧场和其他保护地区的土壤，土壤质量基本保持自然背景水平。Ⅱ类主要适用于一般农田、蔬菜地、茶园、果园、牧场等土壤，土壤质量基本对植物和环境不造成危害和污染。Ⅲ类主要适用于林地土壤及污染物容量较大的高背景值土壤和矿产附近等地的农田土壤（蔬菜地除外），土壤质量基本对植物和环境不造成危害和污染。土壤镉含量统计结果见表5。由表5可知，土壤镉含量为Ⅰ类的样品仅有11个，比例为样品总量的0.884%；绝大部分样品（1 233个）的镉含量在0.2～1.0 mg/kg之间，等级为Ⅱ类，平均值为0.405 mg/kg，占样品总量的99.116%。研究区全部表层土壤的镉元素含量均为0.2～1.0 mg/kg，柴湖镇土壤镉含量基本为Ⅱ类，符合农作物种植和供给畜牧养殖的牧草生长的土壤条件，对植物和环境不会造成污染和危害。

根据RK和GWRK的预测值进行克里格插值得到的土壤镉含量分布见图2。图2显示的值域均存在缩短的现象，RK的预测结果相较于检测值而言值域明显缩短，从[0.248，0.588]缩短到[0.280，0.502]，相比而言，GWRK预测结果的值域则没有明显变化，仅缩短为[0.244，0.540]。按每隔0.1 mg/kg为一个区间统计，RK对土壤镉含量面积分布的预测正确率为75%。其中，利用检测值对土壤镉含量<0.3 mg/kg的预测面积为691.44 hm2，而RK的预测面积为274.88 hm2，其中重合区域面积为273.80 hm2，正确率为40%；检测值预测的镉含量在0.3～0.4 mg/kg之间的区域面积为7 619.00 hm2，RK预测面积为8 574.56 hm2，重合区域面积为5 997.54 hm2，预测正确率达79%；根据检测值预测的镉含量在0.4～0.5 mg/kg之间的区域面积为9 731.35 hm2，RK预测面积为9 724.53 hm2，重合区域面积为7 570.72 hm2，预测正确率达78%；检测值预测的镉含量>0.5 mg/kg的区域面积为534.24 hm2，RK的预测结果中则没有面积分布，预测正确率为0（表6）。

相比RK，GWRK的预测精度要高出很多，合计预测正确率达到89%。其中，土壤镉含量<0.3 mg/kg的区域的GWRK预测面积为705.52 hm2，重合区域面积为586.36 hm2，正确率为85%；含量在0.3～0.4 mg/kg之间的区域预测面积为7 769.39 hm2，重合区域面积为6 836.39 hm2，预测正确率达90%；含量在0.4～0.5 mg/kg之间的区域预测面积为9 771.07 hm2，重合区域面积为8 874.78 hm2，预测正确率达91%；含量>0.5 mg/kg的区域预测面积为329.99 hm2，重合区域面积为301.36 hm2，预测正确率为56%（表7）。

总体来说，无论是采用RK方法还是GWRK方法，都表现出在检测值域中部预测精度较高，而在检测值域两侧预测精度较低的趋势，这是由于处于正态分布两侧的样点距离拟合的函数较远造成的。各含量梯度的面积分布预测结果也显示出GWRK优于RK，这与之前两种方法的评价结果是一致的。

3.3 耕地利用分区

由于GWRK方法的精度高于RK方法，且本研究的重点是在有土壤硒含量数据而缺少镉含量数据的情况下的农业规划，因此，耕地等别的区划基于GWRK方法获得的土壤镉含量分布情况见图3。

耕地资源利用的合理性是由多方面因素共同决定的，针对不同区域土壤硒和镉的含量，需要采取不同的耕地利用策略。低于0.4 mg/kg的土壤镉含量的耕地在《土壤环境质量标准》中属于Ⅱ类标准下含量偏低的区域，通过合理的农业工程措施可以在一定程度上降低土壤镉含量，甚至有可能使部分耕地的镉含量降低到0.2 mg/kg以下，从而成为符合Ⅰ标准的高质量耕地。结合柴湖镇种植小麦、油菜、麦冬、花生、棉花、大蒜等作物的農业生产结构和土壤硒、镉的分布情况，将耕地划分为三类：镉含量较低且非富硒的耕地、富硒耕地、其他耕地，结果见表8，分别采取不同的规划管理方式。

1）镉含量较低（<0.4 mg/kg）且非富硒（<0.4 mg/kg）的耕地，这类耕地是需要重点改良的部分。从统计结果来看，镉含量低于0.4 mg/kg的耕地共5 044.83 hm2，能够满足普通中端农产品生产的需要，而硒含量低于0.4 mg/kg的耕地有6 966.74 hm2，两者重叠的区域达到4 611.70 hm2。施用硒肥并增施石灰能够提高土壤中有效硒的含量，并能明显降低土壤可溶性镉含量[19]。因此，在此类耕地使用增硒降镉技术、提高土壤硒含量的同时降低土壤镉含量是具有实施空间的。通过一定时间的土壤改良，有可能在部分区域达到富硒和镉含量低于0.2 mg/kg的效果。

2）富硒（>0.4 mg/kg）耕地面积达到5 316.57 hm2，是发展富硒农业的根基所在。有研究表明，鄂西南高镉环境同时又高硒，使得镉的毒性大为减小，至今未见镉中毒的报道[20]。这是由于硒具有一定的解毒作用，硒拮抗镉毒性的作用可能与形成Se-Cd复合物有关[21]。研究区土壤硒和镉的含量都远低于鄂西南地区，Ⅱ类的土壤镉含量使生产出的农作物导致食用者发生镉中毒的概率几乎可忽略不计，因此，可以安全地发展富硒农业。一般情况下，富硒农产品较常规农产品价格至少高出1倍[22]，结合柴湖镇目前的农作物种植情况，着力发展富硒水稻、小麦、大蒜、花生和麦冬生产基地。其中，大蒜具有比其他绝大部分作物更强的聚硒能力[23]，而且柴湖镇本来就是远近闻名的大蒜之乡[24]，应考虑成为重点发展对象。

3）其他耕地面积为2 355.85 hm2，虽然不具备富硒禀赋，但作为普通粮食产区以及种植油菜、麦冬、花生、棉花、大蒜等作物是可行的。尤其在这部分耕地集中种植棉花等经济作物可以将原本占用的优质富硒耕地改为进行水稻、花生、大蒜等聚硒能力较强作物的生产，从而能够创造更大的经济效益。

4 小结与讨论

通过构建回归RK和GWRK模型，分析了柴湖镇土壤硒含量和镉含量之间的关系，得到如下结论。

1）柴湖镇富硒耕地面积达5 316.57 hm2，占耕地总面积的43.28%，具备发展富硒农业的土壤地球化学背景的条件。耕地土壤镉背景值处于《土壤环境质量标准》Ⅱ类中含量偏低的水平，基本不会对植物和环境造成危害和污染。

2）通过RK和GWRK模型的建立发现，土壤硒含量与镉含量之间存在显著的正相关关系，利用GWRK预测的土壤镉含量与检测值较为接近，精度明显高于RK的预测结果，可作为区域土壤镉含量评估的依据。

3）RK和GWRK模型的预测结果都存在值域向中值收缩的现象，是造成模拟精度下降的主要原因，而在中值附近预测的精度明显更高。

基于上述结论，对各区域的发展分别提出建议：①利用土样检测数据和分析得出的结论，在低硒低镉的耕地上合理施用硒肥，在提高土壤有效硒含量的同时降低镉的有效性，提升耕地的土壤环境质量，有助于从数量上扩大优质富硒农产品的产区来源。②在耕地利用规划工作中，将富硒耕地用于种植聚硒能力强的食用类作物，非食用类经济作物逐步转移到非富硒且土壤镉背景值相对较高的耕地上，科学高效地利用土地资源。③鉴于在取值较大和较小时对镉含量的预测精度相对偏低的情况，可以在预测结果的高值区和低值区的核心区域对土样的镉含量进行适当数量的测定，用作修正预测结果，从而使其达到更高的拟合精度。

参考文献：

[1] 郦逸根，董岩翔，郑 洁，等.浙江富硒土壤资源调查与评价[J].第四纪研究，2005，25（3）：323-330.

[2] 李家煕，黄怀曾，刘晓端.环境地球化学在农业和生命科学上的应用研究[J].第四纪研究，1999，19（3）：224-230.

[3] MORENO C J，MORAL R，PERREZ E A，et al. Cadmium accumulation and distribution in cucumber plant[J].Plant Nutrition，2000，23（2）：243-250.

[4] 崔力拓，耿世刚，李志伟.我国农田土壤镉污染现状及防治对策[J].现代农业科技，2006（11）：184-185.

[5] 张运强，毛德玲.金华富硒区菜地土壤重金属污染程度模糊综合评判[J].湖南农业科学，2012，12（1）：47-49.

[6] 于 洋，罗盛旭，肖钰杰，等.富硒土壤-蔬菜中硒、镉含量和镉形态的分布及其相关性[J].环境化学，2015，34（4）：798-800.

[7] 魏复盛，陈静生，吴燕玉，等.中国土壤环境背景值研究[J].环境科学，1991，12（4）：12-19.

[8] GB 15618-1995，土壤环境质量标准[S].

[9] BRUNSDON C，FOTHERINGHAM A S，CHARLTON M E. Geographically weighted regression：A method for exploring spatial non-stationarity[J].Geographical Analysis，1996，28（4）：281-298.

[10] KUMAR S，LAL R，LIU D S. A geographically weighted regression Kriging approach for mapping soil organic carbon stock[J].Geoderma，2012，（189-190）：627-634.

[11] HARRIS P，FOTHERINGHAM A S，CRESPO R，et al. The use of geographically weighted regression for spatial prediction： An evaluation of models using simulated data sets[J].Mathematical Geosciences，2010，42（6）：657-680.

[12] JAIMES N B P，SENDRA J B，DELGADO M G，et al. Exploring the driving forces behind deforestation in the state of Mexico （Mexico） using geographically weighted regression[J].Applied Geography，2010，30：576-591.

[13] 郑袁明，陈 煌，陈同斌，等.北京市土壤中Cr、Ni含量的空间结构与分布特征[J].第四纪研究，2003，23（4）：436-445.

[14] LADO L R，POLYA D，WINKEL L，et al. Modelling arsenic hazard in Cambodia： A geostatistical approach using ancillary data[J].Applied Geochemistry，2008，23：3010-3018.

[15] 廖桂堂，李廷轩，王永东，等.基于GIS和地统计学的低山茶园土壤肥力质量评价[J].生態学报，2007，27（5）：1978-1986.

[16] CHANG C C，LAIRD D，MAUSBACH J M，et al. Near-infrared reflectance spectroscopy-principal components regression analyses of soil properties[J].Soil Sci Soc Am，2001，65（3）：480-490.

[17] YANG Y，WU J P，CHRISTAKOS G. Prediction of soil heavy metal distribution using Spatiotemporal Kriging with trend model[J].Ecological Indicators，2015，56：125-133.

[18] 李家熙，张光弟，葛晓立，等.人体硒缺乏与过剩的地球化学环境特征及其预测[M].北京：地质出版社，2000.

[19] 颜送贵，邓正春，郝界州，等.栽培措施对稻米增硒降镉效果研究[J].中国稻米，2014，20（6）：55-58.

[20] 毛大均，郑宝山，严文祥.鄂西南石煤和石煤出露区土壤中硒与镉的含量[J].湖北预防医学杂志，1999，10（2）：1-2.

[21] 杨克敌.微量元素与健康[M].北京：科学出版社，2003.

[22] 张 洋，张 荣，孙小凤，等.富硒大蒜的功效及开发利用研究[J].园艺与种苗，2011（6）：97-98，102.

[23] IZGI B，GUCER S，JACIMOVIC R. Determination of selenium in garlic （Allium sativum） and onion（Allium cepa） by electro thermal atomic absorption spectrometry[J].Food Chemistry，2006， 99：630-637.

[24] 戴光忠.湖北省钟祥市柴湖硒产业示范区建设对策[J].安徽农业科学，2013，41（29）：11861-11862，11865.