基于137Cs示踪法的青藏高原东南部土壤侵蚀状况

张哲茵，衣鹏，c，d，陈鹏

（河海大学，a.水文水资源学院；b.水文水资源与水利工程科学国家重点实验室；c.长江保护与绿色发展研究院；d.全球变化与水循环国际合作联合实验室；e.地球科学与工程学院，南京 210098）

青藏高原作为“世界第三极”，平均海拔超过4 000 m，被称作亚洲水塔，是亚洲众多河流的发源地，为超过14 亿人口提供水源［1］。由于青藏高原处于海拔高、温度低的极端环境，所以其生态环境非常脆弱且自我恢复能力很差［2］。全球气候变暖以及愈加频繁的人类活动对青藏高原本就脆弱的生态环境产生了巨大威胁。此外，青藏高原变暖速率是全球平均速率的2倍，加上降水集中在春夏，青藏高原极易发生土壤侵蚀［1，3-5］。青藏高原东南部是长江、怒江上游源区所在地［6］，相比青藏高原的其他地区，青藏高原东南部水系更多、降水更多、地形复杂，土壤极易发生侵蚀。有学者预测至2050年，青藏高原东南部将是中国潜在土壤侵蚀量最高的地区之一［7］。

土壤碳循环是全球碳循环的重要组成部分，土壤中的碳可能会对全球气候产生影响［8］，有学者认为土壤侵蚀会加速土壤中碳的排放［9］。土壤侵蚀会带来许多危害，不仅会引起生态功能的减弱以及土地生产力的消退［10］，还会对当地的社会、经济、文化方面产生影响［11］。除此之外，还可能产生更严重的后果，例如由于土壤侵蚀河流中的含沙量增加，大量泥沙在下游河道淤积，河床抬高，河道的泄洪能力降低，洪涝灾害发生的风险增加［12］。由此可知，对土壤侵蚀的研究及对土壤动态的科学管理极为重要。

目前主要的土壤侵蚀测量方法有传统测量学方法、地球化学法和GIS 技术方法三大类。传统测量法存在着测量周期长、成本高、监测范围小的问题，而GIS 技术方法也存在着范围不精确、误差大的问题［13］。地球化学法即通过同位素来测量土壤的侵蚀量，目前常用的同位素有10Be、137Cs、210Pbex等。与传统测量法相比，地球化学法测量有着周期短、投入少、不改变原始地貌的优点［14，15］。137Cs 产生于20 世纪的核试验，半衰期为30.08年，对土壤细颗粒有很强的亲和力，不易被植被吸收，只随土壤颗粒运动，因此被作为估算土壤侵蚀量的示踪剂［16，17］。

目前许多学者对青藏高原开展了土壤侵蚀的研究，各类侵蚀测量方法均被使用，有计算总侵蚀量的，也有单独计算水力侵蚀、风力侵蚀和冻融侵蚀的［5，18］。但是大部分研究都是针对流域尺度或是一个小范围的区域（一个盆地或几个县或一段公路途经的区域），很少有学者能针对大范围的区域进行侵蚀量的研究［1，19，20］。有部分学者使用GIS 技术和相关模型对整个青藏高原的水力侵蚀和风力侵蚀状况进行了计算［1，21］，但模型模拟的结果误差较大，不同学者使用同一模型计算同一时间段内同一区域的侵蚀模数的结果也存在较大差异［2］。也有许多学者使用同位素示踪法对青藏高原进行了侵蚀计算，但多位于三江源地区和青藏高原的中部及北部［5，18，22，23］，目前并没有学者对青藏高原东南部的土壤侵蚀状况使用同位素手段进行分析。本研究使用137Cs 示踪法分析青藏高原东南部土壤的侵蚀强度和分布。

1 研究区域概况

研究区域位于青藏高原东南部（91°27′36″E—103°31′12″E，27°8′24″N—34°6′0″N），涉及西藏的东部、青海省的西南部以及四川省的西北部，见图1。研究区域的平均海拔高于4 000 m，多年平均降水量在250～620 mm，多年平均地表温度在5.5～12.6 ℃。区域内水系发达，主要河流有雅鲁藏布江、澜沧江、金沙江以及众多支流河道。区域内的主要土壤类型是高山土和淋溶土，按冻土类型分类则有季节性冻土和少量的多年冻土。研究区域内植被种类以草甸、灌木、针叶林为主。研究区域覆盖了多个生态保护区的部分区域，生态环境比较脆弱，尤其是在海拔超过5 000 m 的地区［6］。

图1 研究区域样点分布

根据中国的三级流域空间分布（https：//www.resdc.cn/data.aspx？DATAID=278），把研究区域分为5 个区域进行侵蚀强度和侵蚀分布的分析。从东到西分别是区域一至区域五，每个区域中分别包含1条一级河流，从左到右分别为怒江、扎曲、金沙江、雅砻江、大渡河。区域一的平均高程为4 643 m，多年平均降水量为400 mm，多年平均地温为8.9 ℃；区域二的平均高程为4 475 m，多年平均降水量为344 mm，多年平均地温为8.2 ℃；区域三的平均高程为4 088 m，多年平均降水量为398 mm，多年平均地温为10.2 ℃；区域四的平均高程为4 147 m，多年平均降水量为448 mm，多年平均地温为8.5 ℃；区域五的平均高程为3 841 m，多年平均降水量为486 mm，多年平均地温为9.7 ℃。

2 数据来源与方法

2.1 样品的采集和测量方法

2.1.1 样点设置及样品采集 采样点多位于地面较平整受人类活动影响较小的区域，以此来反映土壤在自然状态下的侵蚀状况。此外，所设置采样点在研究区域内尽量均匀分布，使其能够较为准确地反映区域侵蚀状况。采样点设置情况具体见图1，共43 个点。2019 年夏季对研究区域内43 个采样点进行了表层0～5 cm土壤的全样采集，在每个点采集4～5个土样进行混合装样，土样的间距在1 m 左右。采样时去除采样点表层植物根系后，使用环刀采样，装入塑封袋并带回实验室分析。对样品进行137Cs 的测量（包含图1 中1～26 号样点），同时测定有机碳（TOC）（包含图1 中1～29 号样点）。此外，对所有点进行土壤粒径的测量。

2.1.2 测定指标与方法

1）137Cs 比活度的测量。在中国科学院南京土壤研究所进行，采用美国ORTEC 公司的高纯锗低本底伽玛能谱仪（探头型号：GWL-120-15；多道型号：DSPEC jr2.0）进行测定。具体步骤为将土样烘干，去除有机质、石块，研磨过200 目筛子，装入测试样品盒，密封10 d，直接放入探测器井，测量时间为80 000 s，137Cs 标准由中国原子能研究院提供，最低检测限为0.12 Bq/kg。

2）TOC 的测量。在河海大学水文水资源与水利工程国家重点实验室进行，采用德国耶拿分析仪器股份公司生产的Multi N/C2100 型总有机碳/总氮分析仪测量。具体步骤为将样品冻干，去除石块草根，研磨过200 目筛子后，称取约100 mg 的土样，加入足量的稀盐酸，经过仪器超过1 000 ℃的燃烧测得样品中TOC 的含量，每隔15 个样质检一次，误差控制在10%以内。

3）土壤粒径的测量。使用LS 13320 型粒径仪，测量范围在20 nm 至2 000 μm。将土样烘干，去除植物根系、石子，经简单研磨后过10 目筛子，加入足量的过氧化氢溶液（30%）进行反应，静置10 h 以上，测量前搅拌均匀，加入粒径仪中测量，测量3 次取平均值。

2.1.3 其他数据来源 多年平均降水量、多年平均地表温度数据来自中国气象数据网（http：//data.cma.cn/data/detail/dataCode/A.0012.0001.html）中研究 区域附近39 个气象站2000—2019 年的均值。研究区域内的土地覆盖数据、归一化植被指数（NDVI）来自中国科学院资源环境科学与数据中心（https：//www.resdc.cn/）2000—2019 年的均值，分辨率均为1 km。数字高程数据（DEM）来自地理空间数据云（https：//www.gscloud.cn/），分辨率为30 m。

2.2 基于137Cs 示踪法的土壤侵蚀计算方法

通过与未发生侵蚀或堆积的点（背景点）的同位素活度比较来计算侵蚀量，若采样点137Cs 的比活度大于背景值，说明该点发生堆积；若采样点137Cs 的比活度小于背景值，说明发生侵蚀［5］。137Cs 随大气沉降后进入土壤，被土壤表层的黏土矿物吸附，随着土壤细颗粒在重力和雨水淋溶作用向下迁移，最终分布深度一般在10 cm 以内［24］。部分学者在青藏高原开展137Cs 的研究，该区域内背景点的137Cs 比活度在1 969～2 376 Bq/m2［18，22，25］，与齐永青等［26］得到青藏高原的背景值812～2 571 Bq/m2相吻合。本研究采用的背景值为实际测定值。对背景点处0～10 cm的土样进行测量，背景点处137Cs 的平均比活度为25.19 Bq/kg，根据土壤干容重一般在1.1～1.6 g/cm3，折算出背景点处的137Cs 的面积活度为2 771～4 030 Bq/m2，考虑到藏东南降水充沛，137Cs 的沉降量较大，因此该背景点的137Cs 比活度较为合理，可以用于后续分析。

研究区域内137Cs 的分布很浅，在背景点处的分布深度只有8 cm，0～5 cm 的土样中的137Cs 占总量的90%以上，故样点处0～5 cm 土壤中的137Cs 能够反映土壤的侵蚀程度。参考李元寿等［5］的模型，假设137Cs在非耕作土中呈指数分布［27］且侵蚀点的剩余部分137Cs分布与未侵蚀点的对应部分的分布相同。

基于上述假设，标准剖面中137Cs 比活度的垂直分布函数如式（1）所示。

式中，f(z)为z（cm）样地深度处137Cs 的比活度（Bq/kg）；a为地表处（z=0）的137Cs 的比活度（Bq/kg），b为137Cs的深度衰减系数（cm-1）。

假设土壤容重沿深度方向不变且各样点的土壤容重相同，则存在一个面积S1，使得面积为S1、高度为1 cm 的土柱的质量为1 kg，则该土柱0～zcm的137Cs比活度的表达式如式（2）所示。

采样点处面积S1的土柱0～5 cm的137Cs 比活度如式（3）所示。

可得：

式中，C0为采样点0～5 cm的137Cs 平均比活度（Bq/kg），h为侵蚀点的侵蚀深度（cm），参数a、b由背景点处137Cs 比活度随深度的分布确定，背景点处137Cs 随深度的分布具体为：0～2 cm 的平均比活度为34.2 Bq/kg；2～4 cm 的平均比活度为65.85 Bq/kg；4～6 cm 的平均比活度为17.3 Bq/kg；6～8 cm 的平均比活度为8.6 Bq/kg；8～10 cm 的平均比活度为0。从而得到C（0）=0，C（2）=68.4 Bq，C（4）=200.1 Bq，C（6）=234.7 Bq，C（8）=251.9 Bq，C（9）=251.9 Bq，C（10）=251.9 Bq。使用最小二乘法进行曲线拟合，可得a=96、b=0.37，见图2。在背景点处，令h=0，可得C0=43.7 Bq/kg。即背景点处0～5 cm 土壤的137Cs 比活度为43.7 Bq/kg。

图2 C（z）-z 的最小二乘拟合

3 结果与分析

3.1 青藏高原东南部137Cs 的分布

研究区域内137Cs 的比活度的范围为0～23.75 Bq/kg，由于低于检测下限而未测出137Cs 的样点共6个，分别是样点4、7、13、18、20、26。测出137Cs 的样点的137Cs 比活度 范围为2.6～23.75 Bq/kg，平均为10.48 Bq/kg，小于背景值43.7 Bq/kg，说明研究区域的土壤整体上呈侵蚀状态。

区域一的137Cs 比活度范围为0～17.96 Bq/kg，平均为7.31 Bq/kg；区域二的137Cs 比活度范围为0～13.67 Bq/kg，平均为8.48 Bq/kg；区域三的137Cs 比活度范围为0～5.63 Bq/kg，平均为2.82 Bq/kg；区域四的137Cs比活度范围为0～18.94 Bq/kg，平均为8.54 Bq/kg；区域五的137Cs 比活度的范围为7.38～23.75 Bq/kg，平均为15.57 Bq/kg。

3.2 基于137Cs 示踪法计算的侵蚀结果

研究区域内，可检测出137Cs 比活度的样点处的137Cs 比活度均低于背景点值，表明所有采样点均存在土壤侵蚀现象。未测出137Cs 的样点是由于该点处侵蚀强度较大，137Cs 比活度低于监测下限，除未能测出137Cs 的样点，区域一的土壤侵蚀厚度范围为2.41～7.63 cm，平均为4.52 cm；区域二的土壤侵蚀厚度范围为3.14～4.82 cm，平均为3.90 cm；区域三的土壤侵蚀厚度为5.54 cm；区域四的侵蚀厚度范围为2.26～7.05 cm，平均为4.29 cm；区域五的侵蚀厚度范围为1.65～4.91 cm，平均为3.23 cm。区域一、区域二、区域三、区域四均存在低于137Cs 测定下限的样点，故实际侵蚀厚度会比计算的结果大。

对于未测出137Cs 的样点，由背景点处137Cs 的最大分布深度为8 cm，未测出的样点由于侵蚀强度过大导致样点处137Cs 的比活度低于监测下限，可得样点处的侵蚀厚度大于8 cm，这里假设其侵蚀厚度为8 cm。土壤干容重一般在1.1～1.6 g/cm3，本试验假设研究区域所有采样点的土壤干容重相同，即都为1.35 g/cm3。通过换算可将137Cs 示踪法计算的侵蚀厚度转换为1963 年以来的平均侵蚀模数，区域一的平均侵蚀模数为13.0 t/（hm2·年），区域二的平均侵蚀模数为11.4 t/（hm2·年），区域三的平均侵蚀模数为16.3 t/（hm2·年），区域四的平均侵蚀模数为12.1 t/（hm2·年），区域五的平均侵蚀模数为7.8 t/（hm2·年）。根据137Cs 示踪法计算的侵蚀模数，5 个区域按平均侵蚀模数从大到小依次为区域三、区域一、区域四、区域二、区域五。

4 土壤侵蚀的影响因素分析

研究区域的年平均温度大于张建国等［28］提出的西藏地区冻融侵蚀区的温度下界-2.5 ℃，故研究区域内冻融侵蚀力不是主要侵蚀力。青藏高原的风蚀力空间上从东南向西北逐渐增大［29］，研究区域位于青藏高原东南部，风蚀力较弱，加上区域内降水丰富，多年平均降水量为250～620 mm，故本区域的主要侵蚀类型为水力侵蚀。

水力侵蚀的影响因素有地形因素、气象因素、植被覆盖和土壤性质［4，30］。将26 个计算侵蚀厚度的点根据植被类型分为草甸和林地（针叶林、灌木）两类。去除未测出137Cs 的样点，余下的样点进行侵蚀厚度与其影响因子的相关性分析，对满足正态分布的因子使用Pearson 相关系数，对不满足正态分布的因子使用Spearman 相关系数，曲线拟合计算相关系数（R），具体结果见图3。

4.1 地形因素

对坡度、高程与土壤侵蚀厚度的关系进行分析，见图3a 和图3b。有研究表明，拉萨河流域的水蚀强度在坡度为15°～25°的区域最大，其次为25°～35°［31］。在本研究区域内也发现了类似的规律，如图3a，虚线椭圆外的3 个点离点群较远，去除这3 个点后使用二次多项式进行拟合（R=0.45），在土壤侵蚀厚度最大的点所对应的坡度约为18°。也就是说，在研究区域内，土壤侵蚀厚度随着坡度的增加先增大后减小，在坡度约为18°的区域土壤侵蚀厚度最大。土壤侵蚀厚度与高程之间没有明显的相关性（图3b）。

4.2 气象因素

对多年平均降水量、地温、风速与土壤侵蚀厚度进行相关性分析。多年平均降水量与土壤侵蚀厚度之间的正相关关系并不显著（图3c），与研究区域以水力侵蚀为主的情况不符。这是因为并非所有的降水都会产生土壤侵蚀［32］，只有产生径流的降雨才会产生侵蚀［33］，所以多年平均降水量与侵蚀量之间没有显著的相关性。在植被类型为草甸的样点处，土壤侵蚀厚度与多年平均地温存在显著的正相关性（R=0.76，P＜0.05）；而在植被类型为林地的样点处，土壤侵蚀厚度与地温之间的正相关性并不显著（图3d）。这可能与融水有关，也有可能与风化作用相关，但根据本次研究的数据很难判断。土壤侵蚀厚度与多年平均风速没有发现明显的关系（图3e），原因是区域内风力侵蚀作用较弱。

4.3 植被覆盖

植被覆盖度的提升能显著加强土壤的水土保持能力［34］，有学者指出在青藏高原的黄河源区，草甸土的土壤侵蚀模数与植被覆盖度之间有显著的负相关关系［5］。这里采用2000—2019 年归一化植被指数（NDVI）的平均值计算植被覆盖度（fc），覆盖度计算见公式（5）。

式中，NDVImin为NDVI的最小值，即纯裸土像元的NDVI；NDVImax指NDVI的最大值，即纯植被像元的NDVI。

在本研究区域内草甸和林地的侵蚀厚度与植被覆盖度间的负相关性都不显著（图3f）。其原因可能有：基于分辨率为1 km的NDVI数据计算的植被覆盖度存在误差；植被分类不够具体，不同植被的根系分布不同，保持水土的能力也不同。

4.4 土壤性质

探讨土壤侵蚀厚度与土壤粒径、有机质含量之间的关系。研究区域内，林地处土壤侵蚀厚度与中值粒径有显著的负相关性，而在草甸处，二者的负相关性不显著，见图3g。本研究采用的是0～5 cm 的土壤测量粒径，这部分土壤会受到侵蚀的影响，侵蚀作用会使得土体颗粒破碎、分散，破坏表层土壤的结构［35］。这部分土不能准确地反映土壤在粒径方面的差异，因此这里不能判断土壤粒径与侵蚀厚度之间的关系。

土壤中的有机质主要存在于黏粒和粉粒中，有机质会降低黏粒和粉粒的胶结能力，从而使土壤团聚体的破坏率增高［4］。因此，土壤有机质含量越高，土壤越容易发生侵蚀。由图3h 可知，研究区域内，在林地，TOC 与侵蚀厚度之间存在显著的负相关性（R=-0.62，P＜0.05），而在草甸处二者的负相关性不显著。在草甸处，TOC 与侵蚀厚度间负相关性不显著的原因可能是样点较少，TOC 的变幅比较小，不容易发现TOC 与侵蚀厚度间的相关性。

5 小结

研究区域内137Cs比活度的范围为0～23.75 Bq/kg，所有采样点都发生土壤侵蚀。研究区域内5 个区域均发生轻度侵蚀，按侵蚀模数从大到小分别为区域三、区域一、区域四、区域二、区域五。

研究区域内，侵蚀厚度随坡度的增加先增大后减小，在坡度为18°附近侵蚀厚度最大。在侵蚀厚度与多年平均降水量、风速、植被覆盖度间没有观察到显著的相关性。林地土壤的侵蚀厚度与TOC 间存在显著的负相关性（P＜0.05），草甸土壤处二者的负相关性不显著。