土壤中重金属元素的迁移、转化及其对植物的毒害和环境的影响程度, 除了与土壤中重金属的含量有关外, 还与重金属元素在土壤中的存在形态有很大关系^[1].铜是植物生长必需的微量元素.土壤中铜主要以固定态存在, 仅有少量的铜以相对活动态存在^{[1, 2]}.土壤有效铜可直接被植物吸收, 是限制作物产量和品质的重要因子^[3].摄入过量的铜元素会对人体健康产生一定危害.过量的铜元素会对植物产生毒害作用^[4], 也可能会通过食物链间接危害人类健康^[5~7].土壤铜有效性比率(即有效铜/全铜)能在一定程度上反映土壤铜从固定态转换成有效态的比率^[8], 结合土壤有效铜的空间分布信息, 可以更有针对性地划定土壤有效铜的风险调控区域.

在土壤和农业生态系统中, 土壤铜从固定态向有效态转化的比例不是恒定的^[9], 其通常受到多种土壤因子的影响, 如土壤全铜、土壤pH和土壤有机质(SOM)等^[10~12].土壤pH通常会影响土壤中重金属元素的形态变化^[13~15], pH越低, 有效态铜的含量越高.而SOM作为与土壤生产力密切相关的指标^[16], 它主要通过与土壤铜形成络合物, 进而影响土壤铜的有效性^[17].因此, 本研究讨论了这3个主要影响因子(即土壤全铜、pH和SOM)对土壤有效铜的影响.

以往针对土壤重金属有效态与影响因素之间关系的研究通常基于传统的全局回归模型, 如普通最小二乘回归(OLSR)^[18~21].但土壤环境因子都是空间变量, 且具有一定的空间变异性.这种总体回归模型, 没有考虑地理坐标^[22], 所估计的回归系数难以反映局部空间区域内土壤重金属有效态与其影响因子之间的空间局部回归关系.地理加权回归(GWR)是一种空间局部回归技术, 该模型通过将数据的地理位置嵌入到线性回归模型中, 并采用核函数为每个空间位置处的解释变量分配不同的权重, 以此来估计因变量和解释变量之间的空间局部回归关系.

本研究以张家港市为研究区域, 主要目的在于: ①探讨该区域土壤类型和土地利用类型对土壤有效铜的影响; ②采用普通克里格预测土壤有效铜及有其有效性比率的空间分布状况, 并据此划定土壤有效铜的风险调控区域; ③采用GWR模型探索研究区内相关影响因子(即, 土壤全铜、pH和SOM)对土壤有效铜的空间非平稳影响, 以期为区域土壤有效铜的精确调控提供空间决策支持.

研究区域位于江苏省南部的张家港市(120°22′~120°49′E, 31°43′~32°01′N)(图 1).该市地处长江三角洲冲积平原, 地势平坦.张家港市隶属于北亚热带季风气候, 四季分明, 雨水充沛, 年均降雨量和年均气温分别为1 039 mm和15.5℃.全市的土壤类型南北分异明显, 北部沿江的圩田地区主要为潮土, 由长江冲积物母质发育而成; 南部的平田地区主要为水稻土, 由泻湖相沉积母质发育而成.另外其他类型, 如黄棕壤、堆叠土和园田土等在该市亦有零星分布.经调查发现张家港市的土地利用类型主要有水田、旱地和城乡聚落(包括城镇建设用地和农村聚落).全市总占地面积为999 km², 其中耕地面积约410 km², 主要以稻-麦或稻-油轮作为主要的种植模式.

图 1

Fig. 1

图 1 土壤和土壤利用的空间分布示意 Fig. 1 Soil sample locations, soil type, and land-use type

采样点是根据张家港市的土壤类型、土地利用现状及采样点在研究区分布均匀性原则布设的.全市共采集了357个土壤样点(图 1), 用GPS详细记录每个采样点的经纬度, 以及采样点周围的环境信息.在每个采样点400 m²范围内, 采集5~6处耕作层(0~15 cm)土壤, 混合均匀后再缩分至1~2 kg.所有的土样在室温下风干, 用木棒捣碎后过10目筛, 另再取部分土样(约50 g), 用玛瑙研钵研磨使其通过100目筛.然后将处理好的样品存于聚乙烯瓶中, 以备后续使用.主要测定土壤有机质、pH、全铜和有效铜.土壤有机质采用重铬酸钾外加热法测定, 土壤pH值采用玻璃电级法在土:水=1:2.5时测定, 土壤全铜含量用反王水消化后用原子吸收分光光度法测定, 土壤有效铜的含量在用DPTA浸提后用原子吸收光谱法测定(3510, 安捷伦-上海分析仪器有限公司), 具体的测定方法参见文献[23].

普通克里格(OK)是一种对区域化变量进行空间预测的方法.该方法利用已知土壤样本值并考虑样本点的空间相关性来预测未采样位置处的土壤属性^[24], 具有无偏、最优的特点^[25].用OK来预测土壤属性的空间分布, 可以有效减小估计误差和调查成本^[26]. OK模型的公式如下：

(1)

且满足：

(2)

式中, 其中Z^*(x₀)是位置x₀处的预测值; Z(x_i)是位置x_i处的观测值; λ_i是从普通克里格系统中位置x₀邻域范围内的第i个样本点对应的权重; n是x₀邻域范围内用于插值的样点数目.

传统的OLSR模型假定观测变量具有随机性并有相同的频率分布, 其因变量是自变量的多元线性函数, 可表示为：

(3)

式中, y为因变量, x_i为第i个自变量, β₀为模型截距, β_i为区域内第i个自变量的回归系数, p为回归项的个数, ε为高斯误差项.这类模型具有完备的理论体系和统计推断方法, 可以用来分析目标变量和解释变量之间的关系^[27].然而, 这种类型的回归模型是非空间的, 也就是说, 这类回归模型在对回归参数估计时不考虑地理坐标^{[22, 28]}.因此, 所估计的回归参数无法直观地反映因变量和自变量之间的局部空间回归关系.

地理加权回归(GWR)模型是对OLSR模型的一种扩展, 它将数据的地理位置嵌入到回归参数中, 通过每个空间位置处的变量分配不同的权重, 在探索空间数据之间非平稳关系时具有明显的优越性^{[29, 30]}.其公式为：

(4)

式中, μ为区域内不同的空间位置, y(μ)为位置μ处因变量的值, x_i(μ)为位置μ处第i个自变量的值, β₀(μ)为位置μ处的模型截距, β_i(μ)是位置μ处第i个自变量的回归系数, p是回归项的个数, ε(μ)是位置μ处的随机误差项.其位置μ处的回归系数的估值公式为：

(5)

式中, X是自变量的一个m×(n+1)矩阵, 其中m是观测数据的数量, n是自变量的数量; Y是因变量的一个m×1的矩阵; W(μ)是一个m×m对角矩阵, 它是通过对空间上更接近校准位置的数据位置分配更多权重来计算的核函数.本研究采用自适应核函数进行权重估计.基于更正后的Akaike信息准则(AICc)选择最优自适应带宽, 并且最终使用71个周围数据点来校准每个位置处的GWR模型.关于GWR模型更多的介绍参见文献[29].

本研究采用ArcGIS 9.2对空间数据进行数字化、空间插值和制图等.土壤基本理化性质的描述性统计分析在SPSS 19.0中完成.土壤有效铜与其相关控制因子的GWR分析在GWR 4.0中完成.

张家港市2种土壤类型(潮土和水稻土)和3种土地利用类型(水田、旱地和城乡聚落)下土壤有效铜、全铜、pH和SOM的描述性统计量见表 1、2.可以看出, 张家港市有80%土壤样品的全铜含量超过了江苏省土壤背景值23.4 mg·kg^{-1 [30]}, 说明土壤铜积累广泛存在于研究区域的农田土壤中.根据中国土壤环境质量标准的二级标准^[31], 仅有一个土样的全铜含量超过了其阈值, 这表明研究区大多数农田土壤铜低于污染水平, 但仍存在潜在污染的风险, 应重点加强监测.

表 1 (Table 1)

表 1 不同土壤类型下土壤有效铜与其影响因子的描述性统计量 Table 1 Descriptive statistics for available copper and its impact factors for different soil types 项目 潮土(n=247) 水稻土(n=112) pH SOM /g·kg-1 有效铜 /mg·kg-1 全铜 /mg·kg-1 pH SOM /g·kg-1 有效铜 /mg·kg-1 全铜 /mg·kg-1 极小值 4.40 11.00 2.00 13.90 4.03 10.60 2.29 7.55 极大值 8.35 31.00 19.00 68.90 8.25 34.40 17.51 226 均值 8.05 19.11 5.91 31.76 6.37 25.10 6.65 31.61 标准差 0.29 3.61 1.64 7.73 0.85 4.33 1.95 20.36 变异系数 0.04 0.19 0.28 0.24 0.13 0.17 0.29 0.64

表 1 不同土壤类型下土壤有效铜与其影响因子的描述性统计量 Table 1 Descriptive statistics for available copper and its impact factors for different soil types

表 2 (Table 2)

表 2 不同土地利用类型土壤有效铜与其影响因子的描述性统计量 Table 2 Descriptive statistics of available copper and its impact factors for different land-use types 项目 极小值 极大值 均值 标准差 变异系数 pH   4.03   8.35   6.62   0.88 0.11 水田(n=272) SOM/g·kg-1  11.00  34.00  20.58   4.59 0.22 有效铜/mg·kg-1   2.00  19.00   6.09   1.77 0.29 全铜/mg·kg-1  13.90 226.00  32.59  14.35 0.44 pH   5.36   8.24   7.35   1.02 0.14 旱地(n=59) SOM/g·kg-1  11.00  31.00  21.68   4.76 0.22 有效铜/mg·kg-1   3.00  18.00   6.52   1.97 0.30 全铜/mg·kg-1   7.55  47.50  28.63   7.13 0.25 pH   4.40   8.25   7.85   1.14 0.17 城乡聚落(n=26) SOM/g·kg-1  11.00  31.00  23.87   5.27 0.22 有效铜/mg·kg-1   4.00   8.00   5.95   1.18 0.20 全铜/mg·kg-1  20.54  42.90  30.11   6.41 0.21

表 2 不同土地利用类型土壤有效铜与其影响因子的描述性统计量 Table 2 Descriptive statistics of available copper and its impact factors for different land-use types

研究区域内, 潮土中土壤全铜和有效铜的平均值分别为31.76 mg·kg^-1和5.91 mg·kg^-1, 水稻土中土壤全铜和有效铜的平均值分别为31.61 mg·kg^-1和6.65 mg·kg^-1(表 1).水稻土中全铜的含量与潮土中的含量相当, 但有效铜的含量要高于潮土.其原因可能在于, 相对潮土, 研究区域内的水稻土具有更低的土壤pH和更高的SOM含量(表 1).

表 2分析了不同土地利用类型(水田、旱地和城乡聚落)下土壤pH、SOM和全铜对土壤有效铜的影响.研究区域内, 水田和旱地中土壤有效铜的平均含量分别为6.09 mg·kg^-1和6.52 mg·kg^-1, 均高于城乡聚落(5.95 mg·kg^-1), 这可能与研究区域内水田和旱地中较低的土壤pH有关(表 2).另外, 农田土壤(水田和旱地)的有效铜含量要高于城乡聚落(表 2), 这可能与农田土壤较大的人为活动有关^[32], 也可能与具体的土壤性质有关.因此, 本研究区域内, 土壤类型和土地利用可能通过对其它土壤属性, 如土壤pH、SOM等, 也一定程度上影响了土壤铜从固定态到有效态的转化.

普通克里格(OK)预测的土壤全铜和有效铜的空间分布如图 2(a)和2(b)所示.土壤全铜的高值区主要集中在研究区域的中部和西北部, 低值区主要集中在研究区的东北部和西南部.土壤有效铜与全铜有相似的空间分布.这表明土壤全铜是影响土壤有效铜含量的一个重要因子.然而, 在局部区域两指标的空间分布模式也存在一定的差别, 如在土壤全铜含量较高的中部, 有效铜的含量却较低.这表明土壤有效铜还受到其它一些因素的影响, 比如中部地区多样的土地利用类型.另外, 土壤有效铜在研究区内呈现出一种南高北低的空间分布规律, 这与土壤类型的空间分布相一致(图 1), 说明土壤类型也会对该地土壤有效铜的空间分布存在一定的影响.本研究区域内, 土壤铜的有效性比率在空间上并不是恒定的[图 2(c)], 存在明显的空间非平稳性.其高值区主要集中在研究区的东北、东南和西北部, 低值主要集中在研究区的中部和西南部.土壤铜有效性比率的变化范围为13.56%~29.15%.

图 2

Fig. 2

图 2 土壤全铜、有效铜和有效性比率的空间分布 Fig. 2 Spatial distribution of total copper, available copper, and copper available ratio

较高的土壤有效铜和有效性比率, 说明研究区域存在较大程度的土壤铜污染风险.本研究以土壤有效铜和有效性比率的平均值(分别为5.91 mg·kg^-1和19%)为阈值, 结合张家港市的土壤有效铜和有效性比率的空间分布图[图 2(b)、2(c)], 划定了该市土壤有效铜的风险调控区域如图 3所示, 其分别代表有效铜含量高于5.91 mg·kg^-1且有效性比率大于19%的子区域1和有效铜含量高于5.91 mg·kg^-1且有效性比率小于或等于19%的子区域2.对于子区域1, 可以通过改变有效性转化比率来降低该子区域内土壤有效铜的含量; 而对于子区域2[图 3(b)], 由于有效性比率较低(≤19%), 则不容易通过降低土壤铜转化比率来降低土壤有效铜含量.为了对本研究区内的土壤有效铜进行精确调控, 有必要调查相关土壤因子(即土壤全铜、pH和SOM)对土壤有效铜含量的影响程度.

图 3

Fig. 3

图 3 张家港市土壤有效铜的风险调控区域 Fig. 3 Potential risk areas for available soil copper in Zhangjiagang City

为了调查了3个主要影响因子(即土壤全铜、pH和SOM)对土壤有效铜含量的影响, 首先计算了土壤有效铜与相关土壤因子之间的皮尔逊相关系数(表 3).结果表明, 土壤有效铜与SOM呈极显著相关, 相关系数为r=0.43.相比以前的研究^[33], 本研究中SOM和有效铜呈正相关, 这可能是因为随着土壤全铜含量的增加, 仅有少量铜与SOM发生了螯合.本研究中土壤pH与土壤有效铜呈显著负相关(r=-0.13), 表明pH也是影响土壤有效铜的主要土壤因子之一.与前人^[8]研究所不同的是, 土壤全铜和有效铜没有表现出相关性, 这可能与研究区内土壤pH和SOM较大的空间变异性有关.

表 3 (Table 3)

表 3 土壤全铜和有效铜与相关土壤属性的皮尔逊相关系数1) Table 3 Pearson correlation between total copper, available copper, and related soil properties pH SOM 全铜 有效铜 pH 1 SOM -0.35** 1 全铜 0.02 -0.02 1 有效铜 -0.13* 0.43** 0.04 1 1)*表示P < 0.05, **表示P < 0.01

表 3 土壤全铜和有效铜与相关土壤属性的皮尔逊相关系数1) Table 3 Pearson correlation between total copper, available copper, and related soil properties

土壤有效铜(A_Cu)和其3个主要的影响因子(即, 土壤全铜(T_Cu)、pH和SOM)之间的全局OLSR关系如式(6)所示：

(6)

其中R²=0.20, 这意味着研究区域中约20%的土壤有效铜的变化可以通过这3个独立变量(即, 土壤全铜、pH和SOM)来解释. OLSR分析是非空间的, 忽略了土壤有效铜和其3个影响因子之间关系的局部空间特征.而土壤有效铜及其影响因子都是空间变量, 往往具有一定的空间相关性和空间异质性.传统OLSR方法的前提假设是样本观测之间彼此相互独立, 因此OLSR模型无法对土壤有效铜与相关土壤因子之间的关系进行精确建模.本研究对OLSR模型和GWR模型的拟合精度进行了对比, 结果表明：GWR有较大的R², 较小的AICc和残差平方和(表 4).因此, 相比传统的OLSR总体回归模型, GWR模型能更好的拟合本研究区域中土壤有效铜和相关影响因子之间的回归关系.

表 4 (Table 4)

表 4 普通最小二乘回归(OLSR)模型和地理加权回归(GWR)模型的比较1) Table 4 Comparison between ordinary least squares regression (OLSR) and geographically weighted regression (GWR) models 模型 AICc R2 残差的平方和 OLSR 1 394.37 0.20 1 009.45 GWR 1 359.87 0.47 706.25 1)AICc表示更正后的Akaike信息准则; R2表示决定系数

表 4 普通最小二乘回归(OLSR)模型和地理加权回归(GWR)模型的比较1) Table 4 Comparison between ordinary least squares regression (OLSR) and geographically weighted regression (GWR) models

GWR分析结果显示, 土壤有效铜与其3个主要影响因子(即土壤全铜、pH和SOM)相关的土壤属性之间的关系在空间上是非平稳的, 随着空间位置的变化而变化(图 4).土壤全铜、pH和SOM对土壤有效铜的影响可以通过它们相应的回归系数来解释.通常, 正的回归系数代表正的影响或相关, 负的回归系数代表负的影响或相关^[34], 而回归系数的绝对值的大小反映各属性因子对土壤有效铜影响的强烈程度.本研究中土壤pH对土壤有效铜有负的影响, 且其在研究区域南部对土壤有效铜的影响要高于北部.土壤全铜和SOM对土壤有效铜的影响相似, 高值区主要集中在研究区的西北和东南, 而低值区主要集中在研究区域的东北和西南[图 4(b)、4(d)].这说明相对其它子区域, 在研究区的西北和东南部地区减少土壤全铜和SOM可能会更有效地减小土壤有效铜的含量.值得注意的是, 本研究探索了3个主要的土壤因子(即, 土壤全铜、pH和SOM)对土壤有效铜的空间非平稳影响.可能其它因子的影响并没考虑进回归模型, 如土地利用类型、土壤类型、土壤阳离子交换量、土壤含水量等等.然而, 由于这些因子与本研究所涉及的相关因子往往具有一定的信息冗余, 回归模型中过多的冗余信息容易产生多重共线性.同时, 考虑过多次要的其它土壤因子往往也会增加区域土壤及环境调查的成本.相比传统的OLSR总体回归模型, GWR模型具有更高的拟合精度, 且揭示了区域尺度内主要的土壤因子(即土壤全铜、pH和SOM)对土壤有效铜的空间非平稳影响.

图 4

Fig. 4

图 4 土壤有效铜与3个主要的土壤因子之间的GWR回归系数分布 Fig. 4 Spatial distribution maps of the regression parameters from the geographically weighted regression analysis of available soil copper versus the three main impact factors

决定系数(R²)常用来反映各自变量对因变量的解释情况, 由GWR产生的局部R²的空间分布如图 4(e)所示.从中可以看出, 相对传统的OLSR回归模型, 这3个主要土壤因子能更好地解释土壤有效铜累积的原因.研究区域东北部局部决定系数相对偏低(R² < 0.30), 其可能与该区域土壤有效铜及相关环境因子具有较大的人为扰动有关, 应重点加以监测.

GWR回归分析结果显示, 各主要因子对土壤有效铜的影响存在一定的空间非平稳性.根据各土壤属性因子对土壤有效铜的空间非平稳影响(图 4), 结合划定的土壤有效铜的风险调控区域(图 3), 可为研究区内土壤有效铜调控措施的制定提供相应的空间决策支持.如在研究区的中部和南部提高土壤pH较研究区的北部会更容易降低土壤有效铜的含量; SOM在研究区的东南部对土壤有效铜的影响要明显高于其他子区域[图 4(d)], 因此在东南部减小SOM会较其它子区域更容易降低土壤有效铜的含量; 土壤全铜在研究区域的东北对有效铜的影响较其它子区域大, 抑制该子区域全铜含量的增加, 能够更有效地控制土壤有效铜含量.

(1) 张家港市有80%土样土壤铜含量超过当地背景值, 土壤铜积累广泛存在于研究区域的农田土壤中; 土壤类型和土地利用类型可能通过土壤pH和SOM间接影响土壤有效铜的含量. OK预测的土壤全铜的高值位于研究区的中部和西北部, 低值区主要集中在研究区的东北和西南; 通常, 土壤有效铜含量较高的区域有效性比率也较大, 如研究区的东南、东北.但局部地区又存在一定差异, 如研究区的西南.土壤铜有效性比率的变化范围在13.56%~29.15%.根据土壤有效铜和有效性比率划定相应的风险区域, 可以更有针对性地调控土壤有效铜的含量.

(2) 相比传统的OLSR模型, GWR模型具有更好的拟合效果(即较大的决定系数R², 较小的AICc和残差平方和), 且能更有效揭示相关土壤属性(即土壤全铜、pH和SOM)对土壤有效铜的空间非平稳影响.根据划定的土壤有效铜的风险调控区域, 结合GWR回归系数空间分布图, 可以在土壤有效铜含量较高的区域, 如受pH影响较明显的东北部, 受SOM影响较显著的东南部以及受全铜影响较大的西南部, 通过向土壤中施加碱性肥料、减少含铜肥料和含铜杀菌剂等农艺措施来调控各子区域土壤有效铜的含量.

(3) 本研究不仅划定了区域农田土壤重金属有效性风险区域, 而且定量解释了相关因子对其有效性的空间非稳态影响, 所得结果可以为区域农田土壤重金属有效态含量的调控提供空间决策依据.