丹江口水库入库河口土壤肥力综合评价

薛盈杉 ，雷俊山

（1.中国地质大学（北京）信息工程学院，北京 100083；2.长江水资源保护科学研究所，武汉 430051；3.长江水利委员会湖库水源地面源污染生态调控重点实验室，武汉 430051）

我国南方地区水资源量占全国的81%，北方地区仅占19%，水资源分布极不平衡。南水北调中线工程每年向北方调水近百亿方，可有效削峰补枯，改善京津冀等受水区域的生态环境。丹江口水库是南水北调中线工程的水源地，是全国最大的饮用水源保护区。

入库河口良好的植被生态能有效减少污染物进入水库[1]，对保持水土、涵养水源、稳定库岸和拦截地表径流污染物效果明显[2-5]。然而，随着丹江口水库大坝加高蓄水，入库河口反复被淹没，河口植被生态受到严重影响。而入库河口土壤肥力的高低及空间分异对河口植被恢复、湿地空间结构及生态修复具有直接影响[6-10]。

国内外关于土壤肥力的学说与观点很多，尚不统一。多数观点认为，土壤肥力是土壤为植物生长提供营养条件和环境条件的能力，是土壤物理、化学、生物以及环境条件综合作用下形成的总体表现[11-15]。土壤养分是土壤肥力的核心部分[13-15]。国内外许多学者对区域土壤养分进行了大量研究，R.Berndtsson等[16]研究了土壤母质类型、质地、有机质含量、土壤pH和土地利用方式等因素对土壤养分空间分布的影响。J.O.Adejuwon等[17]定量化描述了土壤质量的变化。很多研究采用传统统计方法对土壤肥力进行评价[18-20]。R.Webster，张阳等[21-25]根据土壤养分空间变异性理论和点状克里格空间内插方法获得了田间土壤养分状况连续分布图。

近年来关于流域土壤养分的研究[26-29]，大多侧重于不同土地利用方式对土壤养分含量的影响，以期能够通过合理规划土地利用来实现流域植被的重建以及生态环境的改善。然而，本研究涉及的入库河口区域，在丹江口水库大坝加高蓄水后，水位会出现与自然洪枯规律逆反的反复周期性涨落[30]，新的库岸带还会面临巨大的河水冲刷，严重威胁岸坡的安全[31]，导致该区域不适合通过人工改变土地利用来改善土壤肥力状况。因此，本研究通过分析入库河口土壤养分指标的方式，对土壤综合肥力进行空间分布特征的研究，从而指导区域植被的恢复。

本研究的目的是通过改善入库河口植被，达到涵养水源、稳定库岸和拦截污染物，从而最终实现丹江口水源地生态恢复、持续确保优良水质。对丹江口库区入库河口土壤养分指标进行描述性统计分析，并构建评价土壤综合肥力的理论模型，对入库河口土壤肥力进行评价。研究一方面运用克里金插值将离散的样点转为空间上的连续分布，弥补传统统计学的不足；另一方面对每个因子都进行了空间分布研究，从单因子维度细化了土壤肥力的空间分布。因此，本研究对河口植被恢复、湿地生态构建及生态修复具有参考意义。

1 研究区概况

本文选取青塘河小流域入库河口为研究区，面积约为0.378 km2。高程为147～175 m，土地利用有河漫滩荒地、河漫滩沼泽、人工林地、农田、撂荒地、河岸草地和经果林。研究区植被条件较好，无明显水土流失。青塘河小流域位于湖北省丹江口市习家店镇，系汉江河的一级小支流，从北向南汇入丹江口水库。流域属汉江北岸丘陵岗地，地势北高南低。流域内地质构造比较发育，地形复杂，地貌属于构造侵蚀低山区丘陵，山高谷深，呈树枝状水系。流域内土壤类型主要有石灰土、紫色土和黄棕壤。流域内基本农田少，坡耕地面积大，农业生产单一，肥力和土质差。林地以柏树为主，灌木林主要有杨树、刺槐和紫穗槐等，经济林以柑橘和山楂为主，无成片草场，林草覆盖率35.40%。青塘河流域属北亚热带半湿润季风气候，多年平均气温16.1℃，年均降水量797.6 mm。

2 研究方法

2.1 数据来源

研究区布设31个采样点位，经度范围覆盖111°11′～111°12′E，纬度范围覆盖 32°43′～32°45′N，采样点分布（图1）。采样分析工作有长江水资源保护科学研究所完成，布设采样点时，垂直于河流两岸均匀布设6条采用条带，每条采样条带根据地形均匀布设采样点，兼顾包含各种土地利用（表1）。采样时，根据采样点距离水库库岸的距离，在青塘河主河道两侧定点采样，采集不同植被覆盖下土壤表层0～20 cm的土壤，记录每个采样点的坐标，在样点附近再取3～5个点，以几个点的混合土样作为该点的土壤样品，记录每个采样点的坐标与高程。将所采集土壤分拣出杂物，磨碎后过2 mm和0.149 mm孔筛，用于测定pH值、有机质、碱解氮、速效磷和速效钾等指标。pH值采用pH值计法（1∶2.5），碱解氮采用碱解扩散法，速效磷采用碳酸氢钠法，速效钾采用醋酸铵浸提火焰光度计法，有机质采用外加热重铬酸钾溶液法。为减小分析过程中的误差，增强数据的精确度，每个样品各指标分别进行3次重复测试。

图1 采样点分布图Figure 1 Distribution of sampling points

2.2 研究方法

兼顾指标对植物生长发育的重要性和可度量等因素，选取全氮、全磷和全钾等8个因子作为评价指标（表2）。常规数据统计分析采用SPSS 25.0软件完成，采用因子分析提取主成分分析法（PCA），综合各土壤养分指标对样点区土壤养分进行综合评价，利用综合评价模型计算出各样点的综合土壤肥力IFI值。再选择系统聚类分析法，将研究区样本点的IFI值进行等级划分，确定等级分界点。

表2 研究区土壤养分指标描述性统计分析结果Table 2 Results of descriptive statistical analysis of soil nutrient indexes in the study area

地统计分析，采用空间自相关性分析和半方差函数在GS+9.0软件中计算样区土壤各养分指标空间变异特性，并对理论模型进行拟合，按照决定系数（R）最大、残差平方和（RSS）最小的原则求出各指标的最优半方差函数及其相关参数。在此基础上，利用ArcGIS 10.2软件中的半方差函数对经过模型拟合后的各指标进行普通Kriging插值分析，绘出插值结果栅格图。然后利用ArcGIS 10.2软件中的栅格计算工具，根据主成分分析的建模结果，求出样区IFI值的空间分布图。最后，依据聚类分析结果，对其进行等级划分，并计算各肥力等级所占面积的比例。

3 结果与分析

3.1 土壤养分指标描述性统计

采用SPSS 25.0软件对研究区采样点土壤养分数据进行描述性统计（表2）。根据全国第二次土壤普查养分分级标准[32]，研究区土壤pH值处于中性水平（6.5～7.5）。有机质含量符合三级标准（20～30 g/kg）。全氮含量符合三级标准（1～1.5 g/kg）。全磷含量符合三级标准（0.6～0.8 g/kg）。全钾含量符合四级标准（10～15g/kg）。碱解氮含量符合四级标准（60～90mg/kg）。速效磷含量符合六级标准（＜3 mg/kg）。速效钾量符合一级标准（＜3 mg/kg）。

从变异情况可以看出，pH值变异系数仅为4.99%，呈弱变异，说明研究区内土壤酸碱性变化程度不大，属于中性土壤。有机质含量、全氮含量、全磷含量、全钾含量、碱解氮含量、速效钾和速效磷含量的变异系数为23.43%～86.51%，属于中等变异。可见，研究区土壤养分为中等及以下变异。变异程度从大到小的顺序为：速效磷＞速效钾＞有机质＞碱解氮＞全钾＞全氮＞全磷，其中速效磷数值变化范围较大，变异系数也最大，说明速效磷属性值存在较大的波动，受个别极大值干扰作用强烈。

3.2 土壤养分指标的相关性分析

采用Pearson相关分析对表示土壤养分的各属性指标进行相关性检验表明（表3），pH值除了与速效钾呈现较弱的负相关性外，与其余各养分指标均呈现极显著的负相关关系。其中，pH值与有机质含量的相关性最显著，相关系数达到-0.637。说明在土壤呈酸性的条件下，有机质含量较高。有机质含量除了与全钾含量呈较弱正相关性外，与其他养分含量均呈现极显著的正相关性。全氮含量与全磷含量、碱解氮含量、速效磷含量以及速效钾含量呈极显著的正相关性，与全钾含量的正相关性较弱。全磷含量与全钾含量、碱解氮含量、速效磷含量以及速效钾含量呈极显著的正相关性。全钾含量与碱解氮含量、速效磷含量与速效钾含量分别具有显著、极显著、极显著的正相关性。碱解氮含量与速效磷、速效钾含量分别呈现极显著、显著的正相关关系。速效磷含量与速效钾含量具有极显著的正相关性。

表3 土壤各养分指标Pearson相关性分析结果Table 3 Pearson correlation analysis results of soil nutrient indexes

3.3 土壤养分的主成分分析

对各变量进行KMO和Bartlett球状检验，结果表明，KMO统计量为0.778，同时Bartlett球状检验表明相伴概率P值小于显著性水平0.01，因此拒绝Bartlett球状检验的零假设。检验表明各变量间存在较强的相关性，可能存在冗余信息，需要通过因子分析进行主成分提取[33]。

采用SPSS 25.0软件，将研究区土壤的pH值与其余7个养分含量值进行主成分分析，标准化计算后原始数据矩阵的相关性矩阵的特征值与特征向量，基于特征值大于1的原则，提取出前两个主要因子（表4）。可以看出，成分1、成分2的特征值分别为5.355及1.167，符合提取原则，方差贡献百分比分别为66.942%及14.585%，累计方差贡献率达到81.527%，说明前两个主成分可以很好地表示原始数据的绝大部分信息，反映研究区各采样点土壤的肥力特征。

表4 养分指标提取主要因子结果Table 4 Results of main factors of nutrient index extraction

通过计算获得主成分因子载荷矩阵和成分得分系数矩阵（表5）。可以看出，有机质、全氮、全磷、全钾、碱解氮、速效磷以及速效钾含量在成分1中具有较大的正向载荷，说明这7个土壤养分指标在第一主成分上占有较高的比重。

表5 主成分因子载荷与得分矩阵Table 5 Principal component factor load and score matrix

3.4 土壤养分的IFI空间分布特征

用 x1、x2、x3、x4、x5、x6、x7和 x8分别表示研究区土壤的pH值、有机质含量、全氮含量、全磷含量、全钾含量、碱解氮含量、速效磷含量和速效钾含量这8个因子，根据主成分因子载荷矩阵和成分得分系数矩阵可以得出前两个主要因子的拟合模型[34]：

再根据指数和法的综合肥力计算表达式

Fn（权重系数）计算土壤养分综合得分，通过各主成分的方差贡献率数值，经过归一化，计算其贡献率占总方差贡献率的比例，得到线性表达式的系数矩阵，从而得到本研究最终使用的土壤肥力综合评价模型：

利用上述模型，求出31个样点的综合土壤肥力指标值IFI（表6）。对IFI进行描述性统计分析。可以看出，研究区土壤综合肥力指数IFI的范围为-25.30～39.47，平均值为 7.18，变异系数为173.50%，属于极强变异性，说明研究区土壤综合肥力指数差异较大，土壤肥力受人类活动以及土壤酸碱性影响显著。综合所求出的各样点土壤的肥力指数IFI，可以推断出，IFI综合肥力指标值越高，表明土壤肥力越高，反之则表明土壤肥力低。

表6 土壤肥力IFI的描述性统计Table 6 Descriptive statistics of soil fertility IFI

以欧氏距离为衡量样本间差异大小的依据，采用类平均法进行系统聚类分析，将31个样点按照其综合肥力指数大小分为3类，得到3种肥力级别的样本数及对应的IFI指标值取值范围（表7），3种级别分别代表土壤肥力高、中和低3个等级，作为划分研究区土壤肥力级别的模型依据。

表7 3种肥力级别划分依据及所占比例Table 7 Basis and proportion of three fertility grades

研究表明，研究区31个土壤样本点中，土壤养分属于三级肥力（IFI≤8.38）的样本点有18个，占总样点数的58.06%，表明青塘河入库河口新增淹没地有较大面积区域的土壤肥力属于低水平级别。有22.58%的采样点土壤养分处于二级（8.38＜IFI≤15.82）肥力，19.36%处于一级（IFI＞15.82）肥力，属于中、高水平，表明部分土壤具有一定肥力，可以通过恢复植被实现生态修复。

3.5 研究区土壤肥力指标的地统计分析

用GS+9.0软件，对影响土壤综合肥力指数的8个决定性因子进行半方差分析，计算出各拟合模型所对应的块金值、基台值、变程、块金效应、残差平方以及决定系数，根据残差平方和最小以及决定系数最大的原则选定各因子的最佳理论插值拟合模型。

检验结果表明，pH值、经过对数变换的碱解氮含量和速效钾含量，3个指标的最佳拟合模型均为高斯模型；有机质含量、全氮含量、全磷含量的最佳拟合模型均为线性模型；全钾含量的最佳拟合模型为球状模型，经过对数变换后的速效磷含量最佳拟合模型为指数模型。

依据上述结论，采用ArcGIS 10.2软件，对模型（1）模型（3）计算的土壤肥力综合指标IFI进行普通Kriging插值，获得研究区域土壤综合肥力IFI值的空间分布图（图2～图9）。总体来看，研究区养分含量的总体分布特点为西北方向较低，东南方向较高。这是因为研究区耕作、施肥、水分等原因造成，东南部分地势平坦，水分充足，耕作活动频发，土壤养分偏高。土壤养分与有机质等有关，土壤有机质来源主要依靠作物的根茬、秸秆还田、耕地利用及施肥管理等，同时，其他肥料的施用也会对有机质含量造成一定的影响。

图2 pH值分布图Figure 2 pH distribution图3 有机质含量分布图Figure 3 Distribution of organic matter content

图4 全氮含量分布图Figure 4 Distribution of total nitrogen content图5 全磷含量分布图Figure 5 Distribution of total phosphorus content

图6 全钾含量分布图Figure 6 Distribution of total potassium content图7 碱解氮含量分布图Figure 7 Distribution of alkali hydrolyzable nitrogen content

图8 速效磷含量分布图Figure 8 Distribution of available phosphorus content图9 速效钾含量分布图Figure 9 Distribution of available potassium content

依据表6对IFI空间分布图进行重分类，得到研究区土壤综合肥力等级划分图（图10）。可以看出，研究区综合土壤肥力指数IFI的总体变化趋势为由西北到东南方向增高。其中，达到一级肥力标准的区域集中在东南部。不同等级肥力面积占比分别为：一级肥力占14%、二级肥力占46%、三级肥力占40%。研究区土壤综合肥力主要集中在中等等级，总体肥力属于中下水平。

图10 研究区土壤肥力等级空间分布图Figure 10 Spatial distribution of soil fertility grades in the study area

4 结论

研究区土壤大多属于中等偏低的肥力水平。其中，速效钾含量达到了一级标准，有机质含量、全氮含量和全磷含量均达到三级标准，全钾含量、碱解氮含量仅达到四级标准。

研究区土壤基本呈中性，酸性区域土壤养分含量较高。速效磷变异系数高达86.51%，其余养分变异处于中等强度。

研究区土壤肥力分布由西北到东南方向呈上升趋势。一、二级土壤肥力主要分布在农田、人工林地、河岸草地和经果林等植被覆盖较好的区域，表明研究结果合理可信。在入库河口湿地生态构建中，二、三级肥力土壤应选择以恢复植被生态为主的耐贫瘠、耐淹耐旱的植物种类。并结合不同肥力因子的分布差异，依据不同植物生长对不同肥力因子相应机制的差异，有针对性地选取植物种类。一级区肥力土壤条件较好，选择以具有较高的污染阻控和拦截泥沙能力的植物种类。