土壤盐渍化是一个世界性的生态问题。中国西北、华北、东北以及滨海地区是盐渍土的主要集中分布区^{[1, 2, 3]}。天津位于中国华北平原东北部,其滨海地区属于海退之地,地势低平,土壤长期遭受海水及高矿化度潜水共同作用,土壤盐渍化现象分布广泛^{[4, 5]}。中国目前正在开展土地质量生态地球化学调查与评价工作,为滨海地区盐渍化的防治等提供了丰富的数据资料^{[6, 7, 8]}。但目前该类工作的调查精度无法实现每个农用地块均有实测数据,这就需要利用数据插值来覆盖研究区的全部地块,可是采样密度不同会影响到土壤盐分特征的插值结果^[9],因此需要对采样密度的合理性进行研究。

本文选取天津滨海新区太平镇作为典型地区,通过在工作区进行土壤样品采集、测试与统计分析等工作,研究了天津滨海地区土壤主要盐分离子的相关性; 另外还通过对工作区采样点进行随机均匀的抽稀,探讨了不同采样密度对土壤盐渍化插值结果和评价结果的影响,为合理选择采样密度提供理论与实践基础。

天津滨海地区属暖温带半湿润大陆性季风气候,85%的降水量集中在6—9 月份^[10]。工作区太平镇位于天津滨海新区南端(图1),构造单元处于华北断坳黄骅坳陷中^[11],第四纪以来的构造运动形成一套以河流相和洪泛平原相并夹海相堆积为主的砂、泥质松散沉积物。工作区地势低平,多潟湖、洼淀、湿地等,多次海侵使浅层地下水成为咸水,咸水体底界深180~200 m,潜水水位埋深2.5 m左右,地下水化学类型为Cl-Na型,矿化度大于5 g/L^{[12, 13]}。

图1(Fig.1)

图1 研究区位置图 Fig.1 Location of the study area

根据天津市第二次土地利用类型调查成果,太平镇农用地面积约为95.28 km²,共有农用地块8493 块,农用地主要包括果园、其他草地及设施农用地、水浇地和林地,其中水浇地所占面积最大,为农用地面积的71.10%; 其次为其他草地及设施农用地,占农用地面积的23.26%; 果园占农用地面积的5.34%; 林地仅占农用地面积的0.03%。

土壤样品主要根据天津市第二次土地利用类型调查成果图斑,结合采样网格,尽量在工作区内面积较大的农用地地块中央采集,调查采样精度为1∶5 万。由于降水会淋滤土壤中的易溶盐成分,故本文野外样品采集工作选择在降水量很小的春季完成,时间为2013 年5 月22—31 日。采样时,在预定采样点周围20 m半径范围内,采用十字形确定5处分样点采集,并将分样等份组合为一个混合样,以提高样品的代表性; 若遇新近堆积土、田埂、沟渠、路边、旧房基、粪堆等,或采集的土壤中见有肥料团块等,或土地利用类型发生变化时,要适当移点; 采样时去除地表的杂草、树叶等,但不刮除表层土壤; 垂直采集土壤,保证上下均匀采集,采样深度一般为0~20 cm,果园0~40 cm; 每件土壤样品在室内无污染处自然阴干,用木锤压碎,过20 目尼龙筛后送实验室测试^{[14, 15, 16]}。本次工作共采集土壤样品568 件(包括重复样18 件),采样地块550 块,农用地平均采样密度为6 点/km²。

本文所及样品测试项目为土壤全Cl、全盐量及8 项易溶盐离子(Na⁺、K⁺、Ca²⁺、Mg²⁺、HCO₃^-、CO₃^2-、SO₄^2-和Cl^-),样品测试单位为国土资源部天津矿产资源监督检测中心。土壤全Cl 利用X射线荧光光谱法(XRF)测定; 全盐量采用水浴蒸干法测定; Na⁺、K⁺、Ca²⁺、Mg²⁺和SO₄^2-采用等离子体发射光谱法(ICP-AES) 测定; HCO₃^- 、CO₃^2- 和Cl^- 采用容量法(VOL)测定。测试过程中均采用了标准样、密码样、监控样等多种监控手段保障测试质量。各指标测试结果统计特征见表1。

表1(Table 1)

表1 研究区土壤全Cl和易溶盐特征Table 1 Soil total Cl and soluble salt content characteristics of the study area

表1 研究区土壤全Cl和易溶盐特征 Table 1 Soil total Cl and soluble salt content characteristics of the study area

频数分布反映了样本在总体中的分布情况^[17],利用GeoExpl 软件对该地区土壤盐分离子的含量进行统计分析和分布检验(图2),样本数为568 件。由于区域土壤中CO₃^2-含量几乎为零,因此不进行统计分析和分布检验。从盐分离子频数分布图中可以看出不同盐分离子含量的变异规律,各盐分离子之间有一定的差异,频数分布大致分为近似正态型和正偏态型2 种(不剔除异常点),其中HCO₃^-为近似正态型,K⁺、Cl^-、Na⁺、Mg²⁺、SO₄^2-、Ca²⁺均为正偏态型; 阴离子主要是Cl^-含量为3.687 mmol/kg ± 1.022 mmol/kg,平均占阴离子总量的43.90%,其次是SO₄^2-,含量为9.072 mmol/kg ± 11.356 mmol/kg,占阴离子总量的27.52%; 阳离子主要是Na⁺,含量为25.171 mmol/kg ± 45.627 mmol/kg,占阳离子总量的43.30%。土壤中这种盐分离子组成的差异,主要与该地区土壤母质和地下水中盐分离子的组成有关。

图2(Fig.2)

图2 盐分离子频数分布图 Fig.2 Frequency distribution of salt ions

化学元素的离子半径、化合价、存在形态等具有相似性,它们在植物、土壤、沉积物等生命和非生命体中存在一定的相关性^[18]。本文采用R型聚类分析对土壤盐分离子的8 个量化指标之间的关系进行了研究(图3),可以看出,按照相关系数值为0.45 作为分类标准,8 个盐分离子可以分为3 类,Ⅰ 类:CO₃^2-、HCO₃^-; Ⅱ 类: K⁺; Ⅲ 类: Cl^-、Na⁺、Mg²⁺、SO₄^2-、Ca²⁺。其中Cl^-和Na⁺相关性达0.99,说明天津滨海地区土壤中阳离子以Na⁺为主,阴离子以Cl^-为主,尤其在浅层地下水矿化度显著偏高(达盐水程度)的地区,土壤易溶盐中Na⁺、Cl^-所占比例具有绝对优势,具有滨海盐土的基本发生特征。

图3(Fig.3)

图3 聚类分析龙骨图 Fig.3 Dendogram of clustering analysis

土壤全Cl 含量与降水量、地势以及土壤盐渍化程度有密切关系,一般来说,降水量低、地势地平、土壤淋溶性差的盐渍化土壤,全Cl 的含量高^[19]。绘制工作区土壤全Cl 与易溶盐离子Cl^-含量的散点图,可以看出土壤全Cl 与Cl^-相关性非常好,相关系数R可达0.986(图4)。由于滨海地区土壤易溶盐中Cl^-所占比例具有绝对优势,将本次1∶5 万土壤调查所取得的568 个采样点随机均匀抽稀成30 个样点(约全样点的5%,密度为0.3 点/km²),绘制这30 个样点土壤全Cl 与全盐量含量的散点图,可见即使抽稀,土壤全Cl 与全盐量的相关性也很好,相关系数R为0.974(图5)。根据图5 可建立全Cl 与全盐量的相关方程,将568 个样点的全Cl 实测值代入此方程,推算出对应的全盐量,计算全盐量推算值与实测值的相对误差(RE),若认为RE≤40%为合格,并要求相对误差合格率≥ 80%,则568 个样本相对误差合格率可达92.6%,完全满足相对误差合格率的要求。因此对于滨海地区大规模采集的区域土壤地球化学样品,在进行全量分析的基础上,可随机均匀抽取少量样本加测全盐量,并建立全Cl 与全盐量的相关方程,利用该方程可以比较准确地推算出所有样本的全盐量,这样既可以达到较高的精度,又不会因为过多的测试全盐量而浪费时间、人力和财力。

图4(Fig.4)

图4 土壤全Cl与易溶盐Cl-含量散点图 Fig.4 The content scatter diagram of soil total Cl and soluble salt Cl-

图5(Fig.5)

图5 抽稀试验中土壤全Cl与全盐量含量的散点图 Fig.5 The content scatter diagram of soil total Cl and total salt content in thinning test

本次工作共研究7 种农用地采样密度(表2),其中采样密度为6 点/km²的数据是本次1∶5 万农用地地球化学调查获得的太平镇568 个样点数据,其余采样密度的数据是利用网格法,分别对1∶5 万农用地地球化学调查采样点进行随机的均匀抽稀而获得。

表2(Table 2)

表2 不同采样密度全盐量变异系数统计Table 2 Statistics of variation coefficients of total salinities in different sampling densities

表2 不同采样密度全盐量变异系数统计 Table 2 Statistics of variation coefficients of total salinities in different sampling densities

全盐量的变异系数(CV)可以描述工作区土壤盐分在空间上变化的强度,空间变异性差异与采样密度有一定的关系,农田土壤样本的CV>100%,则认为是强变异^{[20, 21]}。对工作区不同采样密度下全盐量的变异系数进行统计(表2),结果表明工作区土壤盐分数据的离散程度很高,空间变异为强变异。但工作区不同采样密度下的全盐量变异系数之间的差距并不大,仅当采样密度降至1 点/km²及以下时,全盐量变异系数才有所变化,说明本次样点的随机均匀抽稀,对全盐量数据的空间离散性影响不大。

土壤属性空间插值的常用方法包括克里格法(Kriging)、反距离权重法(IDW)和最近邻点插值法(NN) 等,但对于常见插值方法的结果比较,Panagopoulos 等认为Kriging 法对于每种土壤属性来说均是最好的方法,但也有学者研究表明IDW法优于包括Kriging 法在内的其他方法^{[22, 23, 24, 25, 26, 27, 28, 29]}。这说明对于特定工作区的特定数据,哪种插值方法更适合,还需要具体问题具体分析。

本文采用金维软件GeoIPAS进行插值计算,选择采样密度为4.1 点/km²(394 点)和6.0 点/km²(568点)的全盐量数据,在同样的像元分辨率下分别选择Kriging 法、IDW 法(幂值分别为2、4)和NN法插值,计算每块地块在2 种采样密度下全盐量插值平均值之间的绝对误差平均值(AE_MN)、相对误差(RE)和相对误差平均值(RE_MN),并统计通过各种插值方法所获得的RE≥40%地块数量的比例(相对误差合格率)。结果表明,Kriging 法的相对误差合格率明显低于其他方法; IDW 法中,较低的幂值产生了较小的绝对误差; NN 法的相对误差合格率与低幂值IDW法基本相同,但NN法产生的误差平均值却明显比低幂值IDW法高,说明在同样的样本条件和像元分辨率下,这两种插值方法均可以很好地反映全盐量的空间变化特征,但低幂值IDW法的插值效果要好于NN法(表3)。综合上述分析,本文选择反距离权重插值法(IDW)作为研究采样密度适宜性的插值法,幂值为2。

表3(Table 3)

表3 不同插值方法插值精度比较Table 3 Comparison of interpolation precisions of different interpolation methods

表3 不同插值方法插值精度比较 Table 3 Comparison of interpolation precisions of different interpolation methods

对表2 中7 种采样密度的全盐量数据进行插值,确保工作区每个农用地地块内均有插值数据,提取各地块内插值数据的平均值,计算各采样密度下每块地块内不同采样密度全盐量插值数据平均值与1∶5 万全样点(568 点)数据及插值数据平均值的相对误差(RE),统计所获得相对误差的平均值和变异系数,并计算RE≤40%地块数量的比例(相对误差合格率)。结果表明,随着采样密度的下降,相对误差在上升,当采样密度降至1.7 点/km²以下时,相对误差的平均值超过100%; 相对误差变异系数随着采样密度的下降而下降,表明较低的采样密度使相对误差的分散性变小(图6); 较高的采样密度可以获得较大的相对误差合格率,当采样密度大于4 点/km²时,相对误差合格率在80%以上(图7),因此若要满足相对误差合格率≥ 80%的要求,采样密度不应低于4 点/km²。

图6(Fig.6)

图6 不同采样密度RE分散度统计 Fig.6 Statistics of RE dispersions in different sampling densities

图7(Fig.7)

图7 不同采样密度RE合格率变化图 Fig.7 The chart of RE acceptance rate in different sampling densities

由于各地情况差异很大,故存在许多根据各地区自身特点所制定的土壤盐渍化分级标准^[30]。本文依据国家环境保护总局2003 年实施的《生态功能区划暂行规程》中的华北地区土壤盐渍化分级标准,将研究区土壤盐渍化程度划分为5 个等级。根据该标准,以工作区农用地地块为成图单元,可利用本次1∶5 万土壤调查所取得的568 个全盐量数据对工作区土壤盐渍化程度进行评价(表4,图8)。

表4(Table 4)

表4 土壤盐渍化分级标准及评价Table 4 The classification standard and evaluation results of soil salinization

表4 土壤盐渍化分级标准及评价 Table 4 The classification standard and evaluation results of soil salinization

图8(Fig.8)

图8 各采样密度下样点分布及盐渍化评价图 Fig.8 The chart of sample distribution and salinization evaluation in different sampling densities

分别采用表2 中各采样密度数据依据表4 中的分级标准,插值绘制工作区土壤盐渍化评价图,计算不同采样密度各等级面积与1∶5 万全样点数据所获得的图件中各等级面积的相对误差(RE),以对比分析不同采样密度对土壤盐渍化评价的影响。结果表明: 在各种采样密度下,工作区的盐渍土空间分布趋势基本是一致的,但采样密度较低时,样点数据的偶然性对图面效果影响较大,不能体现局部地区盐渍化等级的差异性,随着采样密度的提高,评价图面层次越来越明显,插值的精细度也越来越高(图8); 工作区盐渍化等级为一级的面积约占农用地总面积的50%,在各级别中占有绝对优势,采样密度的变化对其面积变化影响相对较小,其他等级的面积分别占农用地总面积的10%左右,采样密度的下降使得这些等级面积的相对误差变化较大(表5); 当采样密度为1 点/km²时,已经出现个别等级面积RE >40%的情况,因此用于评价盐渍化等级的采样密度不应低于1 点/km²。

表5(Table 5)

表5 不同采样密度各盐渍化等级面积的相对误差Table 5 The relative errors of the salinization level area in different sampling densities

表5 不同采样密度各盐渍化等级面积的相对误差 Table 5 The relative errors of the salinization level area in different sampling densities

通过在滨海地区太平镇进行土壤盐分离子的相关性及采样密度研究表明:

(1)天津滨海地区土壤易溶盐中Na⁺、Cl^-所占比例具有绝对优势,土壤全Cl 与Cl^-、全盐量具有很好的相关性,故可以仅均匀抽取少量样本测试全盐量,而其余样本全盐量利用全Cl 推出,这对于有效解决土壤盐渍化样品采集和测试中研究精度与经济性、效率性的平衡问题具有重要的意义。

(2)工作区土壤全盐量空间变异为强变异,利用网格法随机均匀抽稀样点对全盐量数据的空间离散性影响不大; 最近邻点插值法和幂值为2的反距离权重插值法均可以很好地反映本文工作区全盐量的空间变化特征,但后者的插值效果要好于前者。

(3)各地块全盐量插值结果的误差随采样密度的减少而增加,若要满足相对误差合格率≥80%的要求,采样密度不应低于4点/km²。

(4)大密度采样有利于提高土壤样品的代表性,其评价结果能够更好地反映土地盐渍化差异性的空间分布规律,随着采样密度的降低,样点数据的偶然性对图面效果影响较大,用于评价盐渍化等级的采样密度不宜低于1点/km²。

致谢: 本文是在“天津市重点地区农用地质量地球化学调查评价(一期)”项目所获得的部分资料基础上加工整理而成,直接参加该项目具体工作的还有刘金成、侯佳渝、董路阳、李硕、侯建涛、胡志明、吴海燕、张岩等,市地调院、矿产资源监督检测中心、地质事务中心、滨海新区规划和国土资源局、天津国土资源局有关领导和专家对该项目工作给予了悉心指导与大力支持; 审稿专家、编辑老师杨艳对论文提出了宝贵的修改意见,在此一并致以诚挚的谢意!