滩涂围垦区是受海陆交互影响的半人工半自然生态系统, 是我国重要的后备土地资源, 具有很强的开发利用潜力。围垦在缓解人口压力、保持耕地动态平衡等方面起着重要作用^[1-2]。近50 a, 我国滩涂围垦总面积达1.1×10⁴~1.2×10⁴ km², 成为滨海围垦力度最大的国家^[3]。土壤质量是“土壤实现各种功能的能力”^[4-6], 农学上侧重维持自然植物生长的能力^[7], 评价土壤质量是研究农业可持续发展的重点之一。滨海滩涂土壤具有高盐低养、地下水矿化程度高等特点, 不适宜种植农作物, 需要在围垦前期实施脱盐补肥措施人为提高土壤质量^[8]。这种人为干扰极大地改变着垦区环境。垦区土壤质量偏低制约了土地的开发利用^[9], 滩涂围垦之后的土壤质量演变倍受关注^{[3, 10]}。

受土壤空间变异性影响, 评价土壤质量时, 参评因子的取舍至关重要。Andrews等^[11]和Yemefack等^[12]首次提出以主成分分析中最高因子荷载值选取最小数据集(MDS)作为评价最小单元, 广受学者认可。李桂林等^[13-14]在MDS基础上引入矢量常模(Norm值), 既能最大程度代表所选参数又尽可能减少信息损失, 但指标仅包括土壤物理性质和养分, 涵盖信息较单一。贡璐等^[15]、邓绍欢等^[16]引入生物指标丰富了最小数据集, 但未全面涉及土壤基本理化性质。岳西杰等^[17]依托GIS叠加多种地形部位、植被覆盖和施肥等构建了土壤质量评价MDS, 却忽视了时间因素的影响。已有研究建立的最小数据集多集中于土壤元素和地形地貌因素, 而参评因子随围垦时间和土地利用的变化与轨迹等方面的研究甚少。

评价土壤质量多基于隶属函数^{[3, 10]}, 但隶属函数中指标同一化转折点的取值均由专家打分确定, 存在明显的主观性; 且在计算土壤质量综合指数时, 多采用层次分析法或主成分法来确定权重, 其中两两判断矩阵的建立由专家确定, 主成分定权通常只取荷载较高的若干因素, 因此, 定权也具有一定主观性。已有研究对于内部参数的选择仅通过公式换算稍加完善, 鲜有学者引入新的模型进行评价。

本文选用综合荷载Norm值并引入浙江省慈溪市海涂围垦区土地利用和年限变化计算综合得分, 筛选出评价最小单元MDS, 结合滩涂围垦区受自然演替和人为活动影响强烈的特点, 丰富指标体系, 以达到既最大程度代表所选参数又尽可能减少信息损失的目的, 并采用变异系数加权Topsis方法, 从时空二维角度客观评价围垦区土壤质量, 揭示阻碍研究区土壤质量提升的障碍因子。

1 材料与方法 1.1 研究区概况

浙江省慈溪市共修筑11条围垦海堤^[18-19]。本文研究垦区位于慈溪市东南部(32°07′~32°19′N, 121°22′~121°34′E), 面积约114 km², 包含5条主要围垦海堤。该区属北亚热带南缘季风季候, 温和湿润, 光热充足, 年均气温16.3 ℃, 年均日照1 933.5 h, 年均降水量1 325 mm。主要土壤类型为滨海盐土、灰潮土和水稻土, 母质为滨海相沉积物。研究区地理位置及土地利用情况见图 1。

1.2 数据获取

野外采样于2013年7月完成, 按研究区面积通过网格法均匀布点, 共采集276个土壤/沉积物样品。采样时使用手持GPS定位, 在目标点10 m范围内, 五点法采集0~20 cm的表层土壤/沉积物, 混合后形成一个土壤样品。土壤样品经自然风干后过10目筛储存备用。根据遥感影像上不同时期的海堤海岸线确定围垦分区, 按围垦年限自海洋向陆地方向分别标记为Ⅰ区、Ⅱ区、Ⅲ区、Ⅳ区、Ⅴ区、Ⅵ区; 按照土地利用方式, 样品取自农田111个、菜地47个、果园48个、林地17个、养殖27个、荒地21个、潮滩5个(图 1, 表 1)。

土壤样品理化测试项目包括粒径比(粉砂粒/黏粒)、pH、有机质(OM)、N、P、K₂O、CaO、Na₂O、MgO、Al₂O₃、Fe₂O₃、SiO₂、Cl、S、Mn、Br、B、Cu、Zn、Pb、Cd、As、Hg和Cr。采用多目标区域地球化学调查规定的分析方法及质量监控方案^[20], 由江苏省地质调查研究院测试应用研究所、安徽省地质实验研究所承担分析测试。

1.3 土壤质量评价体系建立

根据慈溪市滩涂围垦区土壤理化性质(养分、常量与微量元素、土壤质地等), 通过主成分分析(PCA法)提取因子荷载, 计算矢量常模(Norm值), 通过多变量方差分析(MANOVA)、均值多重比较(S-N-K比较法)等方法结合围垦年限和土地利用方式建立围垦区土壤质量的最小数据集。

综合荷载Norm值可弥补PCA中单一化导致的信息缺失。Norm值的几何意义为变量在多维空间中矢量常模的长度, Norm值越大, 表明变量对所有主成分的综合荷载越大^[11], 涵盖信息就越综合。

$ {N_{ik}} = \sqrt {\sum\limits_1^k {\left( {u_{ik}^2 \times {\lambda _k}} \right)} } $

(1)

式中: N_ik为第i个变量的Norm值, 指在特征值λ > 1的前k个主成分上的综合荷载; u_ik为第i个变量在第k个主成分上的荷载; λ_k为第k个主成分的特征值。

采用多变量方差分析(MANOVA)衡量土地利用变化对各土壤参数的影响, 调整后的决定系数R²表征利用方式对各参数的影响大小。均值多重比较(S-N-K比较法)估算围垦年限对各土壤参数的影响, 用最早产生显著影响的年限T作为时间对研究区土壤理化性质产生影响的界限^[14]。

为使最终因子得分具有可比性, 对Norm值、R²、T进行简单线性变换, 其中Norm值和R²为正指标, 其值越大表明对土壤参数的影响越大, 因此, 最大值赋为1, 其他变量除以相应的最大值; T为逆指标, 进行相反转化。求和得到每个土壤参数的总分值。

特征值λ > 1的各个主成分中, 若某因子荷载 > 0.5则归为该主成分组, 若某参数同时在两个主成分中载荷 > 0.5, 或在每一个主成分中荷载都 < 0.3, 则单独归为新的一组。为进一步减少数据冗余, 抽取分组后在最高得分10%范围内的参数, 观察其与小组间其他参数的相关性(r), 若高度相关(r > 0.5)且在0.01置信水平下显著(^**), 则只选取得分最高的进入MDS, 反之则全部进入MDS。

1.4 土壤质量评价方法

Topsis法是一种逼近理想解的多属性排序法^[21], 它根据比较样本与最佳解和最劣解的相对距离来评价样本的优劣。该方法对样本容量无特殊要求, 且不受参考序列选择的干扰, 已广泛应用于水质综合评价和农业综合开发项目评价^[21-25], 但在土壤质量评价中应用甚少。Topsis法主要步骤为:

① 指标极差标准化

$ {\rm{正指标:}}{a_{ij}} = \frac{{{x_{ij}} - \min \{ {x_{ij}}\} }}{{\max \{ {x_{ij}}\} - \min \{ {x_{ij}}\} }} $

(2)

$ {\rm{逆指标:}}{a_{ij}} = \frac{{\max \{ {x_{ij}}{\rm{\} }} - {x_{ij}}}}{{\max \{ {x_{ij}}\} - \min \{ {x_{ij}}\} }} $

(3)

$ {\rm{适中指标:}}{a_{ij}} = \frac{{\max \{ {x_{ij}}\} - \left| {{x_{ij}} - {x_{0j}}} \right|}}{{\max \{ {x_{ij}}\} }} $

(4)

式中: max{x_ij}、min{x_ij}、x_0j分别为j指标中的最大值、最小值、适中值。

② 确定最优与最劣样本

根据标准化矩阵A所有样本实际数据, 由各指标的最优值与最劣值, 确定最优样本A⁺与最劣样本A^-。

$ {\rm{最优样本:}}{A^ + } = (a_{i1}^ +, a_{i2}^ +, \cdots, a_{in}^ + ) $

(5)

$ {\rm{最劣样本:}}{A^ - } = (a_{i1}^ -, a_{i2}^ -, \cdots, a_{in}^ - ) $

(6)

式中: α⁺、α^-分别为正、负理想解。

③ 变异系数定权

$ {\omega _j} = \frac{{{V_j}}}{{\sum\limits_{j = 1}^n {{V_j}} }}, \;\;{V_i} = \frac{{{S_j}}}{{\overline {{x_j}} }} $

(7)

式中: V_j为第j个指标的变异系数, S_j为第j个指标的标准差, ${\bar x_j}$为第j个指标的平均值。

④ 计算各土壤样本与最优及最劣样本的加权欧氏距离D⁺与D^-

$ D_i^ + = \sqrt {\sum\limits_{j = 1}^n {{\omega _j}{{(a_{ij}^ + - {a_{ij}})}^2}} } $

(8)

$ D_i^ - = \sqrt {\sum\limits_{j = 1}^n {{\omega _j}{{(a_{ij}^ - - {a_{ij}})}^2}} } $

(9)

式中: D⁺与D^-分别表示第i个评价对象与最优方案及最劣方案的距离, a_ij表示第i个评价对象在第j个评价指标上经极差标准化后的取值, ${\omega _j}$表示第j个指标的权重。

⑤ 计算土壤质量综合指数C_i

$ {C_i} = k \times \frac{{D_i^ - }}{{D_i^ + + D_i^ - }} $

(10)

式中: C_i表示第i个样本土壤质量综合指数, 即每个土壤样本与最优土壤样本的接近程度; k取100。C_i取值便在0到100之间, C_i越接近100, 说明该土壤样本越接近最优土壤质量, 土壤质量越高; 反之, 土壤质量越差。

1.5 障碍因子诊断模型

评价土壤质量可直观反映围垦对土壤的改良效果, 有针对性地消除阻碍土质提高的障碍因子对促进围垦区土壤发育有重要意义。本文用Topsis评价公式(8)中“与最优方案的距离D⁺”即可表示第i个土壤样本的障碍度, 几何意义直观, D⁺越大则障碍度越高, 越需对该样本进行改良。

2 结果与分析 2.1 围垦区土壤基本属性

研究区土壤总体呈弱碱性, pH均值8.3, 随围垦时间呈先上升之后陆续下降的趋势, 变异系数(CV)3.8%, 表明在人工围垦后土壤pH变化不大。Al₂O₃、Fe₂O₃、SiO₂的平均含量分别为128.90 g·kg^-1、50.39 g·kg^-1、636.89 g·kg^-1, 符合弱碱性土壤特性。选取粉砂粒与黏粒的比值作为土壤机械组成参数^[13], 其变异系数高达55.11%, 差异较大, 表明研究区土壤机械组成受人为活动影响显著。参考土壤比重换算方法^[26], 土壤容重均值为1.474 g·cm^-3, 略高于全国大多数农用地区, 但仍属于结构较好的土壤, 适宜发展农业。有机质平均含量13.4 g·kg^-1, 参照全国第2次土壤普查土壤肥力分级^[27], 肥力状况属中低等水平。

Cl受海水影响极大, 具有高度变异性(CV > 200%), 且含量与围垦时间高度相关, 沿海堤逐步变化, 说明以农业为主的人类活动加速了滩涂沉积物脱盐过程。受海相动力作用主导, S和Br含量高于全国平均水平, 由海向陆存在较大的差异。

研究区Cu、Zn、Pd、Cd、As、Hg、Cr七大主要重金属含量略高于全国土壤背景值, 虽未达到污染水平, 但Cu、Pd、Cd、As、Hg变异系数超15%。

2.2 最小数据集选取

构建合理的指标体系是Topsis法的关键。在主成分分析(PCA)基础上, 叠加多元方差分析(MANOVA)、均值多重比较(S-N-K比较法)和Spearman相关分析等, 得到最小数据集。

主成分特征值λ > 1的4项主成分累积贡献率高达80.29%, 可解释绝大部分指标的变异性, 按每个参数在各主成分的载荷分成5组。运用公式(1)计算各参数的Norm值, 经线性变化后求和得到综合分值(表 2)。

由表 2可知, 第1组中Pb综合得分最高, 且与组内其他参数相关度高, 因子得分前10%的Al₂O₃(r=0.832^**)、Fe₂O₃(r=0.815^**)与Pb高度相关, 结合Norm值大小需剔除冗余, 故第1组中只选取Pb进入MDS; 第2组类似, 只选取Cd参评; 第3、4组都只含1个参数, 则S和B均进入MDS; OM在第5组中综合得分最高, 与前10%的Cl(r=-0.515^**)高度相关, 但与MgO(r=-0.162^**)和Br(r=0.343^**)的相关度较低, 故在选取OM的同时, 添加MgO和Br进入MDS, 但MgO存在较弱的变异性(CV < 10%), 盐渍土壤脱碱过程中Ca与Mg之间有很高的相似性^[19], 因此用CaO替代MgO更符合围垦区土壤实际状况。综上所述, 最终进入研究区土壤质量评价最小数据集的有效参数为Pb、Cd、S、B、OM、CaO和Br等7项指标。

经围垦年限差异性显著检验, MDS 7项参数的组内变异均小于组间变异, 说明人为围垦对土壤理化性质的影响远大于内部随机误差。结合MDS中不同围垦年限土壤化学性质特征(图 2), 围垦后有机质含量整体呈先下降后上升的变化趋势, 滩涂沉积物有机质平均含量为12.0 g·kg^-1, 属轻度贫瘠^[27], 围垦后, 土壤有机质含量降低至围垦10 a后的8.8 g·kg^-1, 之后逐步增至围垦50 a后的15.2 g·kg^-1且增速趋缓。CaO平均含量整体呈先上升后下降又上升再下降的波动变化趋势, 由沉积物中的40.8 g·kg^-1波动上升至围垦30 a的43.4 g·kg^-1后速降至围垦50 a的29.4 g·kg^-1, 脱钙脱碱高肥现象显著, 适宜作物生长。

常微量元素S和Br的变化相似, 在围垦前后呈先下后升再下的“S型”波动, 组间差异明显, 均值分别从滩涂时期的0.44 g·kg^-1和18.52 mg·kg^-1降至围垦50 a的0.21 g·kg^-1和7.22 mg·kg^-1, 减少了近一半。B的变化恰相反, 呈“反S型”波动, 但总体平稳, 围垦前后从66.81 mg·kg^-1小幅升至68.76 mg·kg^-1。

围垦10 a后, 重金属元素Pb平均含量从滩涂沉积物时32.6 mg·kg^-1下降至22.2 mg·kg^-1, 减少31.9%, 之后的40 a间略升至26.7 mg·kg^-1。Cd虽也在围垦10 a出现最低值0.12 mg·kg^-1, 但在之后40 a间波动上升至0.20 mg·kg^-1, 比潮滩时期增加36.2%, 临近土壤环境质量(GB 15618—2008)1级标准^[28]。

2.3 土壤质量评价

最终进入MDS的7项参数按照土壤理化性质可分成肥力指标(OM)、常量及微量元素(CaO、S、Br、B)和重金属指标(Pb、Cd)(表 3)。其中, OM是正指标, 其值越大, 土壤质量越高; Pb与Cd为逆指标, 其值越高, 污染越严重, 土壤质量越差; 其余指标均为适度指标, 适中值的确定在研究区元素理化性质的基础上综合考虑包括中长江平原区(GBW07429-GSS15)^[29]、土壤成分分析元素含量参考值(GSS1-16)^[29]、浙江省和全国土壤元素背景值^[30]在内的国家标准。

通过Microsoft Excel 2013进行Topsis评价, 先对MDS中的7项指标按式(2)~(4)进行极差标准化, 转化至0~1范围内, 使评价模型低敏感高稳健便于操作。围垦活动受时间影响强烈, 不同海堤内元素含量变化大, 选用变异系数按式(8)定权既能保持客观公正, 又能增强分区内指标的可比性。

采用公式(5)~(9)计算每个样本土壤质量综合评价值C_i, 值越高表示土壤质量越高。研究区276个样本点的土壤质量综合评价值在44.6至74.1之间, 平均值56.1, 最优土壤质量样本是种植水稻(Oryza sativa)的农田, 而最劣土壤质量样本为跳跳鱼养殖塘。

土壤质量整体变化趋势受年限主导, 局部变化特征受利用方式影响(图 3)。研究区土壤质量与年限高度相关(R²=0.98), 50 a间提高20%, 呈“前20 a逐步提高—20~30 a相对平稳—30~50 a急剧提高”的轨迹变化。从利用方式来看, 土壤质量菜地 > 果园 > 农田 > 林地 > 养殖 > 潮滩 > 荒地。农田、菜地和果园的土壤质量综合评价均值达57.6, 由于人为施肥量大, 利用效益明显高于养殖、荒地和潮滩, 说明开发农业对滩涂土壤脱碱增肥有较好的改良作用。

进一步按作物细分为种植瓜果蔬菜[西瓜(Citrullus lanatus)、葡萄(Vitis vinifera)、桃子(Amygdalus persica)、南瓜(Cucurbita moschata)、毛豆(Glycine max)、油菜(Brassica campestris)、韭菜(Allium ordorum)等]的旱地和种植水稻的水田两种。以既有水田又有旱地的Ⅱ区为例, 在时间固定条件下, 水稻田平均土质(C_i=52.1)优于旱地平均土质(C_i=51.9), 恰与研究区整体水田劣于旱地的结果相反(图 3), 且在Ⅳ区出现类似情况, 表明在同一时间内种植水稻的环境效益更优。因此, 当地农民在种植瓜果以追求经济利益的情况下, 应重视保护环境, 合理施用化肥农药, 避免人为输入过量的重金属元素Cd诱发地方性骨痛病等。

林地土壤质量呈现较高均值是因为种植的万年青(Ilex chinensis)、苦楝(Melia azedarach)等绿化树相比于苗圃树花更成熟, 不需要过多的人为照料, 肥力不及农用地, 重金属含量也相对少。土壤质量潮滩略高于荒地, 说明围垦后的土地应充分合理利用, 避免空地荒废, 政府应出台鼓励及补助政策减少土地资源的浪费。

在全国范围内比较, 研究区土壤质量仍呈中度贫瘠状态^[29]。对样本点插值后参考ArcGIS中自然断点法(Jenks), 将土壤质量等级分成5种(图 4):一级土壤最少, 仅8.09 km², 占总面积的7.1%, 大多分布在围垦时间最长的Ⅴ区, 以农田和菜地为主; 二级土壤占19.5%, 广泛分布于最早的3条围垦海堤内, 以菜地、果园、林地为主, 辅之以少量农田和养殖用地; 三级土壤占总面积的20.7%, 主要是农田, 散布于Ⅲ区和Ⅵ区; 占27.4%的四级土壤和25.2%的五级土壤分布在最新修筑的两条海堤内, 为少量农田、林地加上绝大多数养殖、荒地和潮滩。农用地的土壤质量远高于养殖和荒地, 说明近50 a间的围海造田不仅缓解了粮食和人口的压力, 带动了当地的经济, 且对土壤环境以及促进生态系统的发育都有积极影响。

2.4 土壤质量障碍因子诊断

基于Topsis改进后的障碍因子诊断模型, 计算MDS中土壤各项指标的障碍度。研究区整体因子障碍程度排序为OM > CaO > Cd > Br > S > B > Pb, 结合表 4的分区情况可知, 在不同年限和利用方式下, 主要障碍因子不尽相同。总体来看, OM对土壤质量的阻碍占主导; 海相沉积物受OM、CaO和Br影响较大, 围垦前10 a呈现出短暂的低肥高钙变化, OM和CaO阻碍加剧, 土壤逐渐脱离海水环境, Br含量下降不再成为主要障碍; 围垦20~40 a, Cd的负面作用持续加剧, 成为阻碍土质的主要因素; 50 a垦区土壤完成脱盐脱钙达到稳定状态, 肥力持续提高, Cd成为仅次于OM的第2障碍因子, 比上一个10 a提高52.9%, 是未围垦时期的2.9倍, 需引起重视。叠加土地利用方式来看, 旱地虽提高了肥力但Cd含量也随之增加。

3 讨论与结论

采用主成分分析、多元方差分析等方法, 结合Norm值、围垦年限和利用方式, 从24项土壤常规理化性质中筛选出OM、CaO、S、Br、B、Pb、Cd 7项指标作为衡量研究区土壤质量的最小数据集MDS。7项指标全面覆盖肥力、常微量元素和重金属三大方面, 仅用29.2%的指标就包含了80.3%的信息, 既减少了工作量又降低了测试成本, 为后续评价围垦区土壤质量提供了一个简洁全面的模型框架。

滩涂围垦活动下, 土壤性质受人为影响发生明显变化。OM含量先下降后上升的变化趋势, CaO含量在30 a急剧下降后至50 a趋于稳定, 这与已有研究相符。垦区整体仍属轻度贫瘠状态, 说明围垦后开发需经历较长的脱碱增肥阶段, 实现滩涂向农田的转换还需一定的时间积累和人为改良。Cd含量在开发后期快速上升甚至接近污染临界点, 这是已有探究围垦区土壤所忽略的一点, 易富集于水稻的Cd元素需引起广泛的重视, 应从种植区的选择到撒药施肥等阶段全面监测农作物中Cd含量。

相比于传统基于隶属函数计算加权关联度的诊断模型, 变异系数定权的Topsis法减少主观打分且几何意义更加直观。研究区土壤质量与围垦年限高度相关, 呈“前20 a逐步提高—20~30 a相对平稳—30~50 a急剧提高”的趋势。从土地利用方式来看, 土壤质量为菜地 > 果园 > 农田 > 林地 > 养殖 > 潮滩 > 荒地, 说明开发农业对滩涂土壤脱钙增肥有明显改良作用。同一围垦年限下水田土质优于旱地, 表明慈溪滩涂更适宜种植水稻以保护基本农田建设。农业活动导致的增肥脱钙脱盐是围垦区土质提升的根本原因, 但研究区整体土质仍属于中度贫瘠水平, 一级土壤最少, 仅占总面积的7.1%, 质量最差的4、5级土壤超1/2。MDS中, 因子障碍程度排序为OM > CaO > Cd > Br > S > B > Pb, 较低OM含量和较高的CaO、Cd含量是垦区土质提升的主要限制因子。淡水浇灌是初期土壤质量提升的根本原因, 养分提高以及脱钙加剧是后期继续提升的关键, 但50 a间Cd含量增加了近3倍。因此在增加土壤培肥的同时, 合理施用化肥农药是进一步提升滩涂围垦区土壤质量的有效措施。

土壤最小数据集因子的选定目前仍处于摸索阶段, 本文指标选取时未包含生物指标, 也未考虑农业生产过程(如灌溉与施肥)的影响。随年限增加和土地利用变化, 土壤的空间变异性将加强, 因此, 建立一个全面评价围垦区土壤质量演变的指标体系, 以及采样密度对模型稳健性的影响还有待后期深入研究。