土壤生源要素是指土壤层中所含的生物生长发育所需的大量元素, 主要包括土壤养分在内的碳(C)、氮(N)、磷(P)、硅(Si)、硫(S)、钾(K)等元素^[1]。土壤层是陆地生态系统生源要素的主要蓄积库, 生源要素的空间分布、富集矿化、分散流失和迁移转化等生物地球化学循环过程, 与山水林田湖草湿陆地生态系统服务功能的维持和发挥密切相关^[2—4]。生态化学计量学主要研究生态系统能量和多重化学元素的平衡及元素平衡对生态系统交互作用的影响^[5—6]。基于生态化学计量的土壤生源要素C、N、P、K和S的空间分布及影响因素研究, 对揭示陆地生态系统空间演变和生源要素的循环和平衡机制, 预测土壤肥力变化, 指导土地利用优化具有重要意义^[7—8]。

承德地处国家生态文明建设重大战略布局、京津冀协同发展和脱贫攻坚三大战略的交汇节点, 2019年获批为国家可持续发展议程创新示范区, 定位于建设京津冀“城市群水源涵养功能可持续发展示范区”, 发挥着京津冀水源涵养功能和防风固沙生态安全绿色屏障作用^[9—11]。承德北部坝上高原属浑善达克沙地防风固沙带, 中部属燕山太行山水土保持生态涵养区, 南部属燕山山地水环境保护区。承德市是重要的农业、蔬菜和水果基地, 同时也是我国第二大钒钛磁铁矿资源基地, 在矿产资源开发趋紧和绿色发展背景下, 生态农业已成为地区经济发展的重要产业^[12—13]。多年来, 承德市土地利用类型和生态服务功能变化显著, 一方面北侧塞罕坝机械林场等地区植被覆盖显著增加, 另一方面全区建设用地增加, 水体面积减少^[14—16], 对区域土壤生源要素空间分布具有重要影响。与此同时, 承德市位于京津唐重要水源地密云水库、潘家口水库和大黑汀水库上游, 水环境氮素、磷素和硫酸盐面源污染一直是区内水源涵养的重要影响因素^[17—18]。土壤生源要素氮磷硫的空间分布和侵蚀风化流失过程研究对承德市“城市群水源涵养功能可持续发展示范区”的建设具有重要意义。生源要素是土壤养分循环的核心和水体富营养化的重要物质来源^[19], 土壤生源要素的生态化学计量, 空间变异及影响因素研究, 对全区宜林宜耕宜草土地利用适宜性评价, 土地质量和林业生产力提升, 指导特色农业和精准农业发展, 可持续发展议程创新示范区建设, 生态屏障和水源涵养功能的维持和改善具有重要意义。因此, 本文选择承德地区为研究区, 通过资料收集、样品采集和分析测试, 基于生态化学计量理论运用描述性统计、地统计学、空间自相关分析和冗余分析, 阐述土壤生源要素碳氮硫磷钾的空间分布格局和垂向迁移特征, 探讨土地利用类型、土壤类型和成土母质类型及人类活动对土壤生源要素空间分布的影响, 为区内国土空间利用优化和生态保护屏障及水源涵养功能维持提供科学参考和基础数据支撑。

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

承德市位于河北省北部, 西南部与京津接壤, 地处内蒙古高原到华北平原过渡地带, 地势总体北高南低, 地貌由北向南从坝上高原、燕山山脉中低山、燕山南麓丘陵逐步过渡为华北平原区, 平均海拔350 m, 最高峰为雾灵山2118 m(图 1)。区内水系发育, 横跨三大流域单元, 东侧属辽河流域, 分布辽河和大凌河水系；中部属滦河流域, 位于天津和唐山水源地潘家口水库和大黑汀水库上游, 占承德全区的72.5%；西侧属潮河、白河和蓟运河流域, 位于北京密云水库上游。本区属典型的中温-暖温带大陆性季风型半干旱半湿润气候, 年平均降水量536 mm, 气象要素具立体分布特征, 南北差异明显。

由于全区横跨三大地貌单元, 南北气候差异明显, 区内土地利用类型、土壤类型和成土母质类型分布亦具有一定分带特征(图 1)。土地利用类型以林地、草地和耕地为主, 其中林地面积占67.97%(乔木林地37.69%, 灌木林地22.80%), 草地面积占12.33%, 耕地面积占10.89%, 另分布有工矿用地占比约0.80%。全区地表植被覆盖度较高处主要有北部塞罕坝机械林场和云雾山, 中南部千层背和雾灵山周边地区。土壤类型以褐土、棕壤、粗骨土和潮土为主, 北部坝上地区分布草甸土、灰色森林土、栗钙土及少量风沙土和沼泽土, 坝上高原至燕山山地过渡带分布棕壤, 中部燕山山地主要分布褐土, 河流沟谷地势低洼处分布潮土^[12]。成土母质类型总体为北部风成沙、玄武岩、安山岩、流纹岩和凝灰岩类, 中部花岗岩、闪长岩、角闪岩和变粒岩类, 北部碳酸盐岩、陆源碎屑岩和片麻岩类, 河谷两侧分布冲洪积物和黄土^[13](图 1)。区内农业耕种方式以旱作为主, 河流干流两侧分布少量水稻田, 种植作物以玉米、高粱、谷子、豆类和水稻为主, 经济树种有苹果、山楂、板栗、核桃、山杏等。

1.2 样品采集与分析

本次研究共采集土壤样品1597件, 其中表层样品1551件, 垂向浅钻样品46件；样品分两个时间段集中采集, 分别为2019和2020年的6—9月。表层土壤样品在考虑全区土地利用类型、成土母质和土壤类型等因素基础上, 总体以网格法进行采集, 保障样品在空间上的相对均匀；采样深度为0—20cm, 选取代表性地块按“S”或“X”形组合4—5处样点混合成一个样品。为保证土壤样品的代表性, 在承德北部坝上农牧交错带, 中南部燕山山地耕地集中区对样点进行适当加密, 并布设5处土壤浅钻, 根据土壤质地结构按10—20 cm间隔分层采集0—160 cm深度样品46件。采样时, 去除碎石、杂物、植物残体后, 将样品用清洁棉布样袋编码保存, 在研磨、过筛加工后进行室内测试分析。土壤全氮(STN)采用半微量凯氏定氮法, 土壤全磷(STP)采用钼锑抗比色法, 土壤全钾(STK)采用火焰光度计法, 土壤全硫(STS)采用燃烧碘量法, 土壤有机碳(SOC)采用氧化热解-电位法, pH采用玻璃电极法测定。室内分析按规范要求加10%空白样与平行样控制, 分析方法准确度和精密度采用国家一级土壤标准物质(GBW07349)控制, 各指标的加标回收率均在国家标准参比物质的允许范围内。

1.3 数据分析方法

在样品分析测试基础上, 通过承德全区土壤生源要素含量统计结果与黄淮海平原、河北省全区、华北玉米种植区生源要素丰度特征及全国第二次土壤普查养分分级标准的对比, 确定区内生源要素的相对丰缺程度。计算各生源要素的比值关系反映各元素的矿化程度和空间耦合变化, 通过典型土壤剖面分析生源要素的垂向迁移转化特征。结合地统计变异函数和全局Moran's I指数描述土壤生源要素的空间自相关特征, 通过土壤生源要素含量与土地利用类型、成土母质和土壤类型的冗余分析(RDA)反映生源要素及其比值与环境因子的相关关系, 识别土壤生源要素空间分布的主要影响因素。

本次研究承德市土壤类型资料来自中国土壤科学数据库1 ∶ 100万土壤类型数据, 土地利用类型采用最新土地三调图斑数据, 归一化差分植被指数(NDVI)采用2019年Landsat 8遥感影像处理和计算。描述性统计和相关分析采用SPSS 21.0完成, 地统计分析采用GS+ 9.0完成, 全局Moran's I指数采用ArcGIS空间自相关工具进行, 生源要素含量与环境因子RDA分析采用Canoco 5.0进行, 各要素含量空间分布图在半方差拟合模型基础上采用ArcGIS普通克里金插值获得。

2 结果与讨论 2.1 土壤生源要素生态化学计量特征

对研究区内1551件表层土壤样品测试数据进行描述性统计和Kolmogorov-Smirnov检验(α = 0.05), 结果如表 1所示。STN、STP、STK、SOC、STS含量和CNPS比值数据经对数转换处理后, 均符合正态分布。表土生源要素平均含量STK＞SOC＞STN＞STP＞STS, 与全国第二次土壤普查养分分级标准^[20—21]相比, STK含量属高(二级)水平, STN和STP属于中上(三级)水平, SOC含量属于中下(四级)水平；土壤pH中位数为7.08。STP、SOC、STS、STN和STK含量变异系数(Cv)均小于1.0, 根据Nielsen标准, 均属于中等程度变异。研究区STK含量变异系数最小, 为0.303, 含量较为稳定, STP和SOC变异系数相对较大, 分别为0.923和0.911, 含量波动范围相对较大。

与土壤生源要素全国平均水平^[20—22]、承德下游黄淮海平原土壤生态地球化学基准值^[23]、华北玉米种植区^[24]及河北省全区表层土壤平均含量^[7]相比, 研究区表土STN平均含量显著高于黄淮海平原基准值0.381 g/kg, 与华北玉米种植区STN含量1.20—1.30 g/kg总体一致, 但低于全国平均水平1.54 g/kg和河北省平均值1.29 g/kg；表土SOC含量显著高于黄淮海平原基准值2.60 g/kg, 略高于全国背景值17.53 g/kg和华北玉米种植区17.00—17.40 g/kg, 显著高于河北省平均值15.43 g/kg；STP、STK和STS则均高于黄淮海平原基准值和全国背景值, STP平均含量显著低于河北省平均值1.23 g/kg^[7]。

相对于土壤生源要素含量描述性统计, 生源要素间的比值关系能更好的阐述各元素的相对富集关系和空间耦合变化, 反映土壤生源要素的矿化和固定生物地球化学过程^[25]。新鲜的外源有机残体多具有较高的碳氮比(C/N Ratio, CNR), SOC累积矿化量随CNR的增加而增加^[4]。研究区表土CNR平均值为15.001, 高于黄淮海平原6.824、河北省平均值13.57^[7]和中国北方半干旱区的其他研究区及全国平均水平^[26—27], 表明区内表土SOC丰度相对高于STN, 土壤碳源和有机质分解速率相对较低。土壤碳磷比(C/P Ratio, CPR)是衡量微生物矿化土壤有机物质释放P或从环境中吸收固持P素潜力的一个指标, CPR较低时表明土壤P矿化程度高, P有效含量相对较高。研究区CPR平均值为33.881, 显著高于黄淮海平原5.029, 高于全国地球化学基准值30.754, 但略低于河北省平均值35.77^[7], 远低于中国陆地平均值52.70^[27], 全球陆地平均水平72.0^[28], 表明研究区土壤P的矿化速率与河北省全区基本一致。区内土壤氮磷比(N/P Ratio, NPR)平均值为2.192, 低于全球陆地平均值5.9及中国平均值3.9^[28], STP含量的改变相对滞后于SOC和STN, 相对稳定。土壤中S的固定是可逆反应, 当SOC与STS比值(C/S Ratio, CSR)小于300—400区间时, 对土壤有机硫矿化有利, CSR高于300—400区间时, 则有利于生物固硫^[29]。研究区表土CSR平均值为92.650, 最大值为378.814, 有利于土壤有机硫的矿化。研究区表土CSR和CNR的变异系数相对较小, 分别为0.517和0.549, 显著低于CPR和NPR的变异系数, 总体较为稳定, 与土壤CNS对人类活动干扰具有一定耐受性, 在特定地理环境下C、N和S元素含量比值往往较为稳定有关^[30]。总体而言, 研究区CNR＜30, CPR＞15, 土壤化学计量比主要受控于C和N, 其中C的影响作用更为明显。

2.2 土壤生源要素的垂向分布

在围场县御道口牧场黄土山(PM01)和大压岱(PM02)两处, 承德中部红旗镇钒钛磁铁矿矿集区(PM03), 滦河河岸滦平县五道河(PM04)和双滦区陈栅子(PM05)处共采集5处土壤剖面, 各剖面生源要素含量随深度分布变化曲线见图 2。

研究区不同土壤类型和土地利用的样品STN、STP、SOC和STS含量具有明显差异, 随着垂向剖面土层深度的增加, STN、SOC和STS含量总体呈逐渐降低特征。不同剖面STP和STK在垂向空间上的分布具有明显差异, PM02和PM04剖面STP含量呈随深度下降趋势, 其它剖面则呈先下降后上升趋势, PM02土壤质地总体为风积粉细砂和亚砂土, STP含量保持在较低水平；PM04和PM05剖面STK含量随深度增加而上升, 其它剖面则呈先升高后下降并趋于稳定的变化特征。PM01和PM05处STN和SOC明显高于其他3处剖面, STN分别为2.815 g/kg和3.014 g/kg, SOC分别达2.721 g/kg和5.086 g/kg, 黄土山表层土壤属湖沼积相成土母质, 天然牧草地枯落物、凋落物和死亡根系等植物组织累积量较高；PM05处属农用地聚集区, 氮肥和有机肥料施用量较高, 使得两处剖面表层STN和SOC含量相对较高。PM01处STN含量在0—60 cm内迅速降低约一个数量级至0.198 g/kg后趋于稳定, SOC则于100 cm处降至0.910 g/kg, 后处于较低水平。PM05处STN和SOC则均于100 cm深度降至最低, 后保持在较低水平。表层土壤N多以硝铵态氮的形式呈聚集状态存在, 在植物根系吸收和重分配作用下, 垂向剖面STN多表现为逐渐降低^[31]。PM02处STN和SOC含量呈波动变化特征, 0—40 cm处于较低水平, 40—60 cm内相对升高, 80—120 cm后又相对升高, 与坝上地区特殊的土壤层沉积韵律有关。御道口地区全新世以来古植被与古气候演变研究表明, 由于气候暖湿和冷干交替变化, 湖沼积亚砂土亚粘土和风积粉细砂交替沉积, 形成特殊的沉积韵律。气候暖湿期, 区域植被覆盖好, 风沙活动弱, 湖沼积层STN和SOC含量则相对较高；反之气候冷干期, 植被覆盖少, 风沙活动强, 风成粉细砂层STN和SOC含量则相对较低^[32—33]。PM03和PM04处深部STP含量逐渐升高, 与其周边成土母质为钒钛磁铁矿赋矿岩体和围岩有关, 磷灰石为钒钛磁铁矿伴生矿物, 使得深部风化程度较低的土壤母质层STP含量相对较高, 其中红旗剖面STP在60—140 cm处含量范围为2.088—2.198 g/kg, 稳定在较高水平。各剖面STK含量在垂向上的空间分异较小, 含量总体较为稳定, PM04处土壤成土母质为黑云斜长片麻岩, 随深度增加土壤风化程度相对降低, 强活动元素K的淋滤流失程度逐渐减小, 故剖面STK含量随深度增加逐渐升高。PM05处STK含量保持在较高且稳定水平, 范围为2.243—2.419 g/kg。

垂向剖面上, CPR、NPR和CSR总体呈逐渐降低趋势, 表土SOC和STN受植被生物量、微生物呼吸作用以及植被根系分布深度影响, 矿化程度相对更高, CPR和NPR亦显著较高。PM04和PM05的CNR值范围分别为13.513—22.266和15.009—19.036, 显著高于坝上地区PM01平均值12.063和PM02平均值15.669, 可能与滦河沿岸耕地氮肥的施用有关。坝上地区PM01和PM02及PM04的CPR值相对PM03和PM05较高, 可能与成土母质P含量较高有关。PM03处CSR明显小于其他剖面, 与硫化物矿床周围土壤S矿化程度较高有关。总体而言, 不同剖面表层土壤受土地利用类型、质地、成土母质等因素影响, 生源要素垂向分布空间分异差别较大。各剖面STN和STS和SOC含量至深层后差异逐渐减小, 其分布受有机质和氮的外源输入影响更为明显；深层STP和STK则主要受成土母质和岩石风化程度影响, 不同成土母岩区剖面STP和STK差异明显。

2.3 土壤生源要素的空间变异特征

地统计学能在有限空间内对区域化变量进行无偏最优估计, 有效解释变量在空间上的分布、变异和相关特征, 通过变异函数理论和结构分析模型定量刻画变量在空间尺度上的随机性、结构性、独立性及自相关性^[34—35]。本次研究在描述性统计基础上, 采用地统计方法得到各生源要素的半方差函数参数见表 2。半方差函数的表征参数主要有拟合理论模型、基台值(Still)、块金值(Co)、块金系数(Co/Still)、变程(range)和决定系数(R²)。基台值表示空间总变异, 块金值表示随机变异。块金系数(Co/Still)描述变量随机变异的占比, 用以度量变量空间相关性的强弱, 揭示空间变异受到自然结构变异和随机人为因素影响的作用程度^[34—35]。块金系数≤0.25时, 变量空间变异以结构性变异为主, 空间自相关性较强；块金系数在0.25—0.75时, 变量的空间自相关程度为中等, 空间变异受到结构性和随机性因素共同影响；块金系数≥0.75时, 变量空间变异以随机性变异为主, 空间自相关程度弱。变程(range)用以表征空间变异的尺度范围, 即空间自相关的作用范围；决定系数(R²)用以衡量半方差函数理论模型的拟合精度。

研究区表土STN和SOC半方差函数最佳拟合模型为指数模型, STK和STS的最佳拟合模型为线性模型, STP的最佳拟合模型为高斯模型。表土STN、STP、STK和STS半方差函数块金系数值均介于0.25—0.75之间, 空间分异受气候、地形、成土母质、土壤类型等结构性因素和土地利用管理、施肥、种植制度和耕作措施随机性因素共同影响。SOC含量的半方差函数块金系数小于0.25, 变程相对最小, 表明SOC具有小范围, 强烈的空间自相关性, 空间分异主导因素为结构性自然因素。STK和STS半方差拟合模型为偏方差模型拟合, 变程相对较大, 决定系数(R²)相对较小, 表明STK和STS空间变异的尺度范围相对较大, 空间相关性相对减弱, 在成土母质和土壤类型等自然因素和矿山开采、耕作措施和农业施肥等人为因素影响下, STK和STS朝均一化方向发展^[36]。STN和SOC半方差函数决定系数(R²)大于0.75, 理论模型的拟合精度相对较高, 能较好的反映出变量的空间自相关性。

由半方差函数获取的STN、STP、STK和STS块金系数相对较为接近, 空间分异影响因素差异不明显, 采用全局Moran's I指数(Global Moran's I index)能进一步阐明生源要素的空间自相关性。全局Moran's I指数常用以揭示变量在有限区域内观测数据分布的相互依赖性, 即空间自相关程度。通过ArcGIS空间自相关分析得到的土壤生源要素Moran's I参数见表 2, 数据标准化后的全局Moran's I指数范围为[-1, 1], 其值为正时表明变量空间正相关, 其值为负时表明变量空间负相关；其绝对值愈接近于1, 空间自相关程度愈高, 愈接近于0则表明变量在空间上愈趋向于随机分布。Moran's I指数Z值用以判断数据标准化后空间自相关性的显著程度, ±1.96、±2.58是空间自相关显著、极显著的分界值。研究区表土SOC存在极显著空间自相关(Moran's I > 0.05, Z > 2.58, P＜0.01), 呈聚集分布特征, 与半方差分析结论相一致。STK全局Moran's I指数为0.146, 相对最大, Z值介于1.96—2.58间(P＜0.05), 空间自相关性显著^[37—38]。STN和STP全局Moran's I指数略小于STK, Z值接近1.96, 具有相对较显著的空间自相关性；STS全局Moran's I指数相对最小, 且Z值小于1.96(P＞0.05), 其空间自相关程度相对最弱。结合半方差分析结果, 表土SOC空间自相关程度极显著, 主导因素为结构性因素；STK空间自相关程度较高, 但自相关作用尺度范围较大, STK在空间分布上较为稳定。STP空间自相关程度相对较显著, 空间异质性较大；STS空间自相关程度相对最弱, 受随机性因素即人为因素影响最为明显。

2.4 土壤生源要素空间分布格局

在半方差分析获得的最佳拟合模型基础上, 采用普通克里金插值得到土壤生源要素和pH值的空间分布图(图 3)。表土STN和SOC的空间分布格局较为一致, 高值区位于北侧塞罕坝、丰宁滦河源山区, 中部隆化东侧茅荆坝自然保护区, 南部雾灵山白草洼、潘家口水库千鹤山和都山自然保护区周边, 均为植被覆盖度相对较高区域。STN和SOC含量总体呈现出北部塞罕坝地区高, 向中部逐渐降低, 至南部兴隆至宽城雾灵山-潘家口水库-都山周边升高且变化较快的分布特征, 其中隆化西北部和平泉西北部STN和SOC含量总体较低、茅荆坝地区相对较高。土壤pH低值区与STN和SOC高值区分布总体较为一致, 植被高覆盖区塞罕坝、茅荆坝、雾灵山白草洼等地土壤枯落物等植物残体累积量大, 微生物群落丰富, 土壤腐殖质层相对较厚, 使得STN和SOC含量相对较高。STP总体呈北低南高分布格局, 主要存在两处高异常区, 一处为围场县牌楼乡至道坝子乡一带, 一处为滦平至双滦的红旗-大庙钒钛磁铁矿矿集区一带, 低值区则主要位于围场御道口牧场等植被覆盖为天然牧草地区域。STK总体呈现出北部和中部高, 至南部后明显降低的分布格局, 高值主要浓集于丰宁西北侧和红旗-大庙一带。全区北部和中部STK含量明显高于南部, 且中北部STK含量较稳定的处于相对高值水平, 与中北部分布火山岩成土母质, 岩土中钾长石等富钾矿物含量较高, 南侧分布贫钾的沉积岩碳酸盐岩成土母质有关。STS高值区多以斑块状分布, 呈现出高值浓集中心较为聚集的分布格局, 中部滦河流域STS由北至南总体呈现出随地势降低逐渐升高特点, 高值区主要为丰宁鱼儿山至黄旗, 滦平至双滦红旗-大庙和南部鹰手营子、兴隆和宽城峪耳崖一带, 多为煤矿、金矿和钒钛磁铁矿等硫化物矿床聚集区(图 1)。

空间分布上, 表土CNR空间自相关程度高, CNR值较为稳定, 与SOC和STN的空间格局耦合程度较高。表土CPR和CSR高值区与SOC高值区分布较为一致, 分布于区内北侧坝上、中部茅荆坝、南侧雾灵山白草洼和潘家口水库周边, SOC高值区外源碳的输入缓解了能源物质对微生物分解速率的限制, 加速了微生物对P和S的矿化, 有利于土壤中有效P和S的增加。表土NPR的空间分异较明显, 其受耕作措施、外源氮肥和磷肥的输入影响较为明显。

2.5 土壤生源要素分布与主要环境因子的关系

地统计半方差分析和全局Moran's I指数仅对土壤生源要素空间分布的结构性和随机性影响因素进行了简要区分, 但对结构性和随机因素的具体内涵及其影响程度未经剖析。为进一步探讨土壤生源要素空间分布的影响因素, 分不同土地利用类型、土壤类型和成土母质类型统计1551件表层样土壤样品5种生源要素含量差异, 得到不同生源要素含量分布箱线和正态分布曲线(图 4)。同时在土壤生源要素含量、比值和pH及各样点处NDVI指数相关分析基础上, 采用冗余分析辨析生源要素含量与环境因子的相关关系, 识别生源要素空间分布格局的主导影响因素。

冗余分析(RDA)可用以描述环境因子对定量属性变量的相关关系, 通过原始变量与典型变量之间的相关性, 分析引起原始变量变异的原因^[39—40]。RDA分析图内指示线表征各变量因子主成分分析(PCA)的特征向量大小, 变量因子指示线间夹角的余弦值表示因子的相关系数^[18]。通过表层土壤生源要素与11种主要土地利用类型、8种主要土壤类型和12种成土母质类型间的RDA分析(图 5)与Pearson相关分析(表 3)表明, RDA分析第一轴和第二轴分别解释了土壤生态化学计量84.90%和7.29%的变量, 生源要素含量及其比值、NDVI和pH值等指标提取的2个主成分特征值方差分别36.68%和12.95%, 具有较好的代表性, 能有效解释变量。STN和SOC呈显著的正相关关系, 相关系数为0.879, 对STN和SOC分布影响相对最大的环境因子为林地(乔木林地、灌木林地及其他林地)和NDVI值, NDVI与STN和SOC的相关系数分别为0.103和0.114, 存在一定的正相关关系。由图 4知, 各土地利用类型STN平均值大小关系为灌木林地＞乔木林地＞其他林地＞草地＞水域滩涂＞交通运输用地＞水浇地和水田＞园地＞住宅用地＞工矿用地＞旱地；其中乔木林地、灌木林地和其他林地STN平均值分别为1.383 g/kg、1.339 g/kg和1.231 g/kg, 其他类型用地STN平均值在1.097—1.198 g/kg之间, 较为稳定；旱地内STN相对最低, 平均为0.941 g/kg。STN正态分布曲线表明, 灌木林地、乔木林地和其他林地STN含量变异系数分别为1.14、0.99和1.02, 波动范围相对最大, 其他用地内则均属于中等程度变异。乔木林地和园地内SOC相对最高, 平均达21.567 g/kg和21.237 g/kg, 在旱地内最低, 平均为13.870 g/kg。NDVI反映了地表植被的覆盖及生长状态, 林地等植被覆盖度相对较高区域凋落物体量明显高于其他用地, 使得土壤腐殖质含量较高, 表土具高STN和SOC。土壤pH受气候、地形、成土母质、植被覆盖和人为活动等众多因素影响, 不同成土母质区, 风成粉细砂、玄武岩、凝灰岩、片麻岩、花岗岩和安山岩区土壤总体呈弱酸性, pH平均为6.59、6.68、6.79、6.86、6.91和6.92, 碳酸盐岩区和第四系风成黄土和冲击洪积物分布区土壤总体呈弱碱性, pH平均为7.31、7.54和7.25。研究区北侧成土母岩以酸性火山岩为主, 与酸性母岩上形成的土壤pH较碳酸盐岩区低有关。不同土地利用类型分区土壤中, 园地、草地、乔木林地、灌木林地和其他林地pH平均为6.64、6.89、6.95和7.07, 相对较低；工矿用地、住宅用地、水域滩涂、水浇地和水田土壤pH平均为7.49、7.42、7.41和7.29, 呈弱碱性。植被覆盖度较高的林地土壤pH值相对较低, pH与SOC和STN呈显著程度较弱的负相关关系, 植物根系对土壤溶液离子的选择性吸收、表层土壤枯落物和腐殖酸含量较高、土壤微生物活性相对较强使得林地土壤pH偏酸性。林地多分布于高坡地形部位, 风化淋溶作用相对较强, pH较低地为低；工矿活动和农业生产生活较强烈区域多分布于宽缓河谷和沟谷地带, 水域滩涂、水田和水浇地相对林地多分布于地势低洼处, 由于承纳高处淋滤成分较多, 盐基离子相对聚集, 土壤pH相对升高^[29]。对STN和SOC分布影响相对较大的土壤类型因子为棕壤和其他土壤, 平均值分别为1.586 g/kg和1.355 g/kg；其他土壤类型主要为沼泽土和新积土, 多位于坝上地区河湖湿地旁, 属自然界生产力最高的生态系统之一, 碳储备丰富^[41]。对STN和SOC分布影响相对最大的成土母质环境因子为碳酸盐岩, 平均为1.543 g/kg, 与其他成土母质相比, 碳酸盐岩含丰富的白云石和方解石等无机碳源, 在微生物和植物根系生物风化作用下, 无机碳源转化为有机碳, 使得区内碳酸盐成土母质区SOC和STN含量相对较高。承德市内塞罕坝林场、茅荆坝和雾灵山自然保护区等林地分布区, 即植被覆盖度NDVI高值区, SOC和STN含量亦相对较高。塞罕坝地区林地面积的增加, 显著增加了土壤碳氮储备, STN和SOC明显高于周边风沙活动区草地土壤。塞罕坝地区沙地和草地向林地的流转^{[15, 42]}, 使得土壤团聚体增加, 对SOC和STN的积累发挥了积极作用, 从风沙治理和水源涵养功能角度来看, 土壤团聚体和植被恢复有利于防风固沙, 增加了滦河源、辽河源及潮白河源头地区水源涵养资源量, 能有效控制土壤氮素随降雨径流的流失, 对坝上高原和燕山山地过渡带农牧交错带退耕还林, 风沙治理具有积极意义。

对区内表土STP分布影响最大的环境因子为成土母质, 不同成土母质STP平均值大小关系为角闪-闪长岩＞片麻岩类＞冲洪积物＞花岗岩类＞流纹岩类＞碳酸盐岩类＞安山岩类＞陆源碎屑岩(砂岩和砂砾岩)＞凝灰岩类＞风成黄土类＞玄武岩类＞风成粉细砂类。角闪岩和闪长岩类、片麻岩类STP平均含量分别为1.034 g/kg和0.837 g/kg, 显著高于其他成土母质区土壤。玄武岩区和风成沙区STP平均含量分别为0.510 g/kg和0.342 g/kg, 坝上草原风沙带STP含量相对较低。承德市大庙-黑山斜长岩、斜长角闪岩岩体为钒钛磁铁矿赋矿岩体, 黑云斜长片麻岩、花岗岩岩体为斜长岩赋矿岩体围岩, 岩石矿物中含有大量的钒钛磁铁矿伴生矿物磷灰石^[43—44], 使得相应成土母岩区STP含量具高地质背景特征。与前文所述红旗-大庙钒钛磁铁矿周围STP具高值区浓集, 土壤垂向分布上具有逐渐升高特征相一致。对STP分布响应程度较高的土壤类型为褐土、粗骨土和潮土, 土地利用类型为水田和水浇地、住宅用地、水域滩涂及工矿用地等, STP高值区多位于河漫滩和沟谷两侧地势较低洼处。一方面, P属于一种沉积性矿物^[45—46], 在矿物风化程度较低的粗骨土中含量相对较高, 粗骨土区STP平均达0.729 g/kg。另一方面, 斜长角闪岩、斜长片麻岩和花岗岩在侵蚀风化过程中淋滤流失, 露天矿山开采使得地表植被覆盖度减小, 一定程度上加速了成土母质中P的流失, 使得地势低洼磷素趋向于汇水处水成土潮土和淋溶褐土STP沉淀聚集, 含量相对升高, 平均分别达0.835 g/kg和0.781 g/kg。土壤pH对STP的淋滤释出存在一定影响, pH与STP含量相关系数为0.119, 总体呈正相关关系, pH值的升高有利于土壤P元素的固定。褐土主要由枯枝落叶层、腐殖质层、黏化层、钙质层和母质组成, 随着土层深度增加, STP含量在垂向分布上呈逐渐增加趋势。STP高值区与工矿用地和农用地分布耦合程度较高, 交通运输用地、住宅用地和旱地、水浇地和水田及工矿用地内STP平均含量相对较高, 平均分别为0.879 g/kg、0.814 g/kg、0.798 g/kg、0.793 g/kg和0.736 g/kg, 交通运输用地和旱地内STP含量变异系数分别为0.920和0.900, 波动相对最大, 表明区内矿业活动和农业活动对STP的分布具有一定影响。

对STK分布响应程度较高的成土母质类型为花岗岩类、变粒岩类和风成黄土类。不同成土母质表土STK含量的变异系数均小于0.1, 属低变异水平, 空间分布相对最为稳定, STK平均值范围为19.757—24.009 g/kg。风成黄土、花岗岩、安山岩和凝灰岩区STK平均含量相对最高, 分别为24.009 g/kg、23.875 g/kg、23.254 g/kg和23.217 g/kg。研究区富钾火山岩广泛分布, 岩石矿物中钾长石含量较高, 如红旗-大庙一带中粗粒钾长花岗岩钾长石成分约占60—65%, 承德中部张家口组凝灰岩和安山岩内钾长成分可达45%—50%^[13], 与此同时风成黄土分布与火山岩分布较为一致, 亦具有较高的K含量, 使得全区北部和中部火山岩分布区STK含量明显高于南部碳酸盐岩地区。对STK分布影响最大的土地利用因子为工矿用地、住宅用地、交通运输用地和草地等植被覆盖相对较低区域, 工矿用地和住宅用地区STK平均含量分别为23.663 g/kg和23.100 g/kg, 其他地类STK平均值则在21.086—22.985 g/kg区间, 较为稳定。全区STK高异常区主要位于大型矿集区周边, 即红旗-大庙黑山一带等植被覆盖度明显降低区域。STK和NDVI呈一定程度的负相关关系, 大型铁矿露天开采裸露基岩和山区交通运输筑路切坡裸露对岩石风化成土过程中K的淋滤释出存在一定影响, 矿区周边和交通运输道路旁等人类活动强烈区域STK含量相对升高。区内栗钙土分布区STK含量显著高于其他类型土壤, 平均达27.622 g/kg。栗钙土主要分布于研究区丰宁西北角, 即浑善达克沙地南缘坝上高原中部, 属温带半干旱大陆气候和草原植被覆盖区, 风沙活动强烈, 昼夜温差大^[42]。与此同时, 栗钙土分布区出露基岩多为钾长花岗岩, 在强烈风蚀、昼夜温差和冻融作用下, 岩石风化程度相对较高, 风化过程中强活动元素K表现为淋滤流失并向地势低洼处迁移聚集。再者, 栗钙土成土过程中腐殖质积累和钙积过程中形成了碳酸钙淀积层或膏盐层, 使得降水难以入渗, 表土由于风化淋溶承纳的K⁺和Na⁺含量较高, 局部土壤具盐渍土特征^[47—48], 故此栗钙土分布区STK含量相对较高。

对STS分布响应程度较高的成土母质类型为碳酸盐岩、片麻岩、角闪岩-闪长岩类和流纹岩。碳酸盐岩区STS含量明显高于其他成土母质区, 平均达0.258 g/kg；片麻岩区STS平均为0.207 g/kg, 角闪-闪长岩类和片麻岩区STS的含量变异系数相对最大, 分别为1.451和1.039, 含量波动相对最大。研究区碳酸盐岩区多含炭质泥页岩和煤质夹层, 含硫量相对较高, 而角闪-闪长岩和片麻岩为钒钛磁铁矿赋矿岩体, 含黄铁矿、黄铜矿、闪锌矿和方铅矿等金属硫化矿物, 故此STS高值区主要位于碳酸盐岩和片麻岩, 角闪-闪长岩和流纹岩成土母岩区。植被覆盖度对STS分布的响应程度亦相对较高, STS和STN、SOC存在一定的正相关关系, 相关系数分别达0.543(P＜0.05)和0.492(P＜0.01), 区内表土CSR值总体均小于400, 有利于有机硫的矿化, SOC和STN含量高值区STS矿化程度亦相对较高。沼泽土-新积土和粗骨土土壤类型区STS含量相对较高, 平均分别为0.288 g/kg和0.232 g/kg, 与沼泽湿地土壤有机硫矿化程度较高, 粗骨土内硫素流失较少有关, 其他土类STS平均在0.163—0.199 g/kg之间, 相对较为稳定。另外, 人类活动强度相对较大的住宅用地、水浇地和水田、工矿用地周围STS含量亦相对较高, 平均分别达0.296 g/kg、0.210 g/kg和0.198 g/kg, 住宅用地周边STS含量变异系数最大, 达1.940, 具有强变异特征。人类活动生活污水, 农业耕种施肥, 矿业活动加速了硫化物氧化淋滤释放过程, 使得住宅用地、农用地和矿山周围STS显著升高。住宅用地和农用地主要分布于河谷地势低洼处, 具有较好的耦合关系, 农业种植施用含硫肥料如过磷酸钙、生活粪肥等有机肥料等, 对硫的矿化有利；矿业活动改变了金属硫化物的原生氧化还原环境, 加速了硫化物的氧化, 且在富含铁、铝氧化物和水化氧化物等矿物的土壤中, 硫酸根吸附量相对最高, 故矿集区周边STS具高异常特征。

山地土壤多为基岩就近风化而成, 成土母质影响着STK、STP和STS的原生背景^[49—50], 丰宁和围场等马铃薯种植区、河谷地带水浇地和水田磷肥的使用, 增加了土壤外源磷素的输入, 钒钛磁铁矿矿山采掘活动加速了裸露基岩侵蚀过程中磷素的淋溶过程, 使得地势低洼处沉积型矿物P相对固定富集。承德北部坝上风沙活动区草甸土STP相对缺乏, 南侧碳酸盐岩分布区STK含量相对缺乏, 在相应地区农业种植活动中可合理施用磷肥和钾肥进行调控。工矿活动造成的基岩裸露和植被破坏一定程度上造成了承德中南部地区水土流失, 加速了岩石侵蚀P和S元素的流失, 是区内水环境和下游水源地总磷和硫酸盐的重要来源, 研究区中南部应注意控制矿山采选, 交通运输筑路等造成的水土流失, 以防区域氮素、磷素和硫素流失对区内水环境及下游水源地水质的影响。

研究区表土CNR值大小主要受SOC含量控制, 其与SOC的相关系数为0.409(表 3), 对CNR分布响应程度较高的环境因子为土地利用类型, 各地类CNR平均大小关系为水浇地和水田＞园地＞乔木林地＞旱地＞住宅用地＞草地＞水域滩涂＞其他林地＞灌木林地＞工矿用地＞交通运输用地, 水浇地和水田CNR平均值为16.342, 变异系数相对最大, 为1.785；园地CNR平均值为17.138, 其他用地内土壤CNR总体均小于15.000。水浇地和水田存在施肥等外源氮的输入^[18], 一定程度上促进了硝化-反硝化微生物的生长, 农作物种植加速了土壤中N的流失；另一方面, 土壤中硝化-反硝化微生物的生长抑制了其他微生物活性, 或者外源氮素与土壤中木质素结合形成了更为稳定的有机物, 降低了土壤中有机碳的利用, 使得CNR值相对升高^[4]。流纹岩区土壤CNR平均为18.840, 明显高于其他成土母岩区, 可能与承德市内流纹岩分布区多形成高陡剥蚀侵蚀山地地貌有关, 流纹岩区地形切割深度相对高于其他成土母岩区, 土壤氮素的流失相对高于其他地区, 使得CNR相对升高。粗骨土区CNP平均值为19.322, 相对高于其他类型土壤, 与粗骨土多分布于流域地势较高处, STN淋滤流失程度相对较高有关。SOC和STN呈显著正相关关系, 故CPR和NRP亦具正相关关系, 对二者空间分布响应程度较高的环境因子与SOC和STP较为一致, 主要为植被覆盖度。与此同时, CSR与SOC相关系数为0.709, 呈显著正相关关系, CSR与SOC高值区分布基本一致, 植被覆盖程度较高区域表层土壤多具有更高的碳氮矿化速率, 故此林地、沼泽土和新积土、玄武岩和碳酸盐岩分布区等STN和SOC高值区CPR、CPR和CSR值亦相对较高。SOC和STN空间分布主要受土地利用类型影响, STP和STS则主要受成土母质影响, SOC和STN含量相对STP和STS含量高出约一个数量级, 植被覆盖和土地利用对STN和SOC的空间分异特征的影响明显高于STP和STS, 区内土壤生源要素化学计量比主要受控于C和N, 其中C的影响作用更明显。

3 结论

承德市表层土壤生源要素平均含量大小关系为STK＞SOC＞STN＞STP＞STS, STK相对丰富, STN和STP属中上水平, SOC则相对缺乏。表土SOC和STN含量显著相关, 空间分布耦合程度高。垂向分布上STN和SOC含量总体呈随土层深度降低特征。表土SOC具有小范围、强烈的空间自相关性, 空间分异主导因素为结构性因素；STN具有较显著的空间自相关性, 空间分异受结构因素和随机因素共同影响。影响SOC空间分异的主要环境因子为土地利用类型, STN分布除受土地利用类型影响较大外, 还受外源氮素的输入等因素影响。STK、STP和STS空间分异受结构和随机因素共同影响, STK空间自相关程度高, 空间分布较为稳定；STP空间自相关程度相对较显著, 空间异质性较大；STS空间自相关程度相对最弱, 受人为影响最为明显。不同地区土壤剖面STK和STP分布具有明显差异, 成土母质高背景区STP和STK含量随深度增加而上升, STS随深度增加逐渐降低。STK和STP空间分异的主要影响因素为成土母质, 受岩石风化影响, 钾长石含量较高的火山岩成土母岩区STK含量相对较高, 钒钛磁铁矿赋矿岩体和围岩成土母岩区角闪-闪长岩和片麻岩区STP含量相对较高。STS空间分异主导因素为成土母质, 其次为工矿活动；硫化物金属矿赋矿岩体角闪岩-闪长岩和片麻岩, 含煤质夹层的碳酸盐岩成土母岩区STS相对较高。