王丽云,王小燕,王 义,岳本江,屈 创
(1.黄河上中游管理局,陕西 西安 710021; 2.神华神东煤炭集团有限公司,陕西 神木 719300)
神东矿区位于晋陕蒙接壤地区,是我国最重要的煤炭生产和能源基地之一[1]。面对大规模开采与极脆弱生态环境保护之间的矛盾,矿区不断创新生产方式,探索生态治理与环境保护的新理念,提出了以控制外围风沙侵蚀为前提,内外围结合治理,促进矿区整体恢复和改善矿区生态环境的思路,走出了一条资源保护型开采与生态环境治理相协调的绿色矿业之路[2]。为全面体现神东矿区2005—2018年水土流失综合治理成果,客观评价矿区水土保持工作成效,本研究采用中国土壤流失方程(CSLE)和风力侵蚀模型,定量计算矿区2005、2010、2015、2018年4期土壤侵蚀模数,并对不同时期不同侵蚀强度土地面积进行对比分析,希望能为矿区水土流失综合治理提供数据支撑。
神东矿区分布于陕西省榆林市神木市北部、府谷县西部,内蒙古自治区鄂尔多斯市东胜区及伊金霍洛旗的南部、准格尔旗的西南部,地理位置为109°51′~110°46′E、38°52′~39°41′N,包括大柳塔矿、活鸡兔矿、哈拉沟矿等14个矿井,总面积约1 064 km2,黄河一级支流窟野河(乌兰木伦河)纵贯矿区中部。地处黄土高原北部与毛乌素沙地交错区,具有典型的干旱荒漠高原气候特征,多年平均降水量仅362 mm,年水面蒸发量为2 297.4~2 838.7 mm[1]。区内大部分为典型的风成沙丘及沙漠滩地,尤其是北部、东南部为黄土丘陵沟壑区,梁峁起伏,沟壑纵横,地表支离破碎,风水蚀交互作用导致水土流失严重。
本研究主要涉及遥感影像、地形数据、气象资料、土壤资料等基础数据。
(1)遥感影像。采用2005、2010年的SPOT5影像和2015、2018年的ZY-3数据进行2005、2010、2015、2018年矿区土地利用和水土保持措施遥感解译。采用监测年前3年空间分辨率为250 m的MOD13Q1数据计算植被覆盖度。采用监测年逐日AMSR-E Level 2A亮温数据产品(空间分辨率为25 km)计算表土湿度数据。遥感影像使用前进行辐射校正、几何校正、配准、影像融合、镶嵌与裁切等预处理,融合后的影像既提高了多光谱影像的空间分辨率,又保留了其多光谱特性。遥感数据统一采用CGCS2000坐标系、Albers投影。
(2)地形数据。基于矿区1∶10 000比例尺的数字高程模型DEM,用符素华等[3]开发的土壤侵蚀模型地形因子计算工具计算坡长、坡度因子等。
(3)其他数据。收集矿区8个站点2000年以来的逐日降水量资料和逐日整点风速资料,用于计算降雨侵蚀力因子和风力因子。植被覆盖度修订系数采用第一次全国水利普查水土保持情况普查成果数据。
2.2.1 水力侵蚀模型
在水力侵蚀区,采用中国土壤流失方程(CSLE)计算土壤侵蚀模数[4-5]。计算公式为
A=R·K·L·S·B·E·T
(1)
式中:A为土壤侵蚀模数,t/(hm2·a);R为降雨侵蚀力因子,MJ·mm/(hm2· h·a);K为土壤可蚀性因子,t·hm2·h/(hm2·MJ·mm);L为坡长因子;S为坡度因子;B为生物措施因子;E为工程措施因子;T为耕作措施因子。
(1)降雨侵蚀力因子。将8个站点多年平均1~24个半月降雨侵蚀力矢量化,利用普通克里金空间插值方法,生成空间分辨率为2 m的24个半月降雨侵蚀力栅格数据。
(2)土壤可蚀性因子。采用第一次全国水利普查水土保持情况普查土壤可蚀性因子成果作为矿区土壤可蚀性因子值。
(3)坡长、坡度因子。基于区域DEM数据,采用北京师范大学符素华等[3]开发的地形因子计算工具计算坡长、坡度因子值,并对林草地采用缓坡公式修正。
(4)生物措施因子。基于MOD13Q1 NDVI数据采用像元二分法计算植被覆盖度,利用植被覆盖度修订系数,获得每年24个半月30 m空间分辨率的植被覆盖度,再结合24个半月降雨侵蚀力因子比例与土地利用类型计算园地、林地、草地的生物措施因子;参照《区域水土流失动态监测技术规定(试行)》对非园地、林地和草地生物措施因子查表赋值。
(5)工程措施因子。根据野外调查和遥感影像解译,矿区水土保持工程措施主要为土坎水平梯田,工程措施因子赋值为0.084。
(6)耕作措施因子。矿区在全国轮作区中涉及两个一级区、三个二级区。一级区北部中高原半干旱喜凉作物一熟区下的二级区后山坝上晋北高原山地半干旱喜凉作物一熟区和陇中青东宁中南黄土丘陵半干旱喜凉作物一熟区,耕作措施因子取值0.488;一级区北部低高原易旱喜温一熟区下的二级区黄土高原东部易旱喜温作物一熟区,耕作措施因子取值0.417。
2.2.2 风力侵蚀模型
在风力侵蚀地区,根据土地利用类型,选用耕地、草(灌)地、沙地(漠)等风力侵蚀模型,计算风力侵蚀模数[6-7]。耕地、草(灌)地、沙地(漠)风力侵蚀模型计算公式分别为
(2)
(3)
(4)
式中:Qfa为每半个月内耕地风力侵蚀模数,t/(hm2·a);W为每半个月内表土湿度因子,介于0~1之间;Tj为每半个月内各风速等级的累计时间,min;Z0为地表粗糙度,cm;j为风速等级序号;Uj为第j等级的平均风速,m/s;Qfg为每半个月内草(灌)地风力侵蚀模数,t/(hm2·a);V为植被覆盖度,%;Qfs为每半个月内沙地(漠)风力侵蚀模数,t/(hm2·a)。
(1)风力因子。风力因子是指不同等级风速和风向对风力侵蚀的潜在能力,与临界侵蚀风速和各等级风速的累积时间有关。风力因子计算采用收集的矿区逐日24 h整点风速和风向资料,按1 m/s间隔统计全年大于等于临界侵蚀风速的各等级风速的累积时间,利用普通克里金方法,插值生成24个半月各风速等级对应的累计时间。耕地的临界侵蚀风速取值为5 m/s,草(灌)地、沙地(漠)的临界侵蚀风速按照李智广等[7]在我国风力侵蚀抽样调查方法研究中不同植被覆盖度的临界侵蚀风速取值。
(2)表土湿度因子。表土湿度因子表征了土壤表层0~2.5 cm深度范围内含水率抗拒土壤风力侵蚀的潜在能力。基于AMSR-E Level 2A亮温数据计算每天的表土湿度因子。
(3)植被覆盖度。计算方法同水力侵蚀植被覆盖度。
(4)地表粗糙度因子。地表粗糙度因子指因植被、微地形和农田耕作导致的零风速位置的高度。在风力侵蚀区,除耕地因呈斑块状分布,地表粗糙度在边界发生突变外,其他诸如草(灌)地和沙地(漠)等的边界地表粗糙度都呈逐渐过渡特征。根据不同土地利用类型,采用李智广等[7]在我国风力侵蚀抽样调查方法研究中不同下垫面粗糙度的赋值。
神东矿区各时期风力侵蚀和水力侵蚀面积见表1。2005年矿区平均土壤侵蚀模数为3 322.53 t/(km2·a),水土流失面积943.36 km2,占矿区总面积的88.66%,侵蚀强度以轻、中度为主。2010年矿区平均土壤侵蚀模数为2 497.58 t/(km2·a),水土流失面积855.47 km2,占矿区总面积的80.40%,侵蚀强度以轻、中度为主。2015年矿区平均土壤侵蚀模数为1 843.98 t/(km2·a),水土流失面积730.66 km2,占矿区总面积的68.67%,侵蚀强度以轻度为主。2018年矿区平均土壤侵蚀模数为852.17 t/(km2·a),水土流失面积587.61 km2,占矿区总面积的55.23%,侵蚀强度以轻度为主。
表1 神东矿区各时期不同侵蚀强度土地面积统计
对比4期土壤侵蚀数据(表2),2005—2010年神东矿区水土流失面积减少了87.89 km2,中度及以上侵蚀面积均有所减少,其中:中度侵蚀面积减少最明显,减少了89.75 km2;其次为强烈侵蚀,面积减少了64.76 km2。2010—2015年水土流失面积减少了124.81 km2,中度、强烈侵蚀面积减少明显,分别减少了117.17、72.41 km2。2015—2018年水土流失面积减少了143.05 km2,中度侵蚀面积减少最多,减少了139.97 km2。
表2 神东矿区各时期不同侵蚀强度土地面积变化统计
总的来看,2005—2018年神东矿区水土流失面积明显减少,共减少了355.75 km2,中度侵蚀面积减少最多,其次为强烈侵蚀,分别减少了346.89、170.96 km2。对比分析矿区不同时期水土流失面积动态变化可知,2005—2018年矿区水土流失状况为高强度侵蚀向低强度侵蚀转移,水土流失面积明显减少,水土流失强度明显降低,生态环境显著改善。
2005—2018年,神东矿区土壤侵蚀模数由3 322.53 t/(km2·a)降低至852.17 t/(km2·a),水土流失面积减少了355.75 km2,相当于矿区总面积的1/3还多。取得如此成效主要得益于矿区始终坚持“生态矿区、环保矿井、清洁煤炭”的建设理念。矿区不仅实施了传统水土保持工程措施及生态水保林、经济林、灌草等林草植被措施,还创造性地利用柳杆障蔽生态锚固坡等技术防治水土流失。截至2018年底,矿区已累计治理水土流失面积339 km2,水土流失和土地沙化得到有效控制,植被覆盖度大幅度提升,由开发初的3%~11%提高到60%以上[8],生态环境得到显著改善。
本研究采用资料收集、遥感监测、模型计算(风力侵蚀、水力侵蚀模型)等方法,定量推算了神东矿区2005、2010、2015、2018年4期土壤侵蚀模数,并对不同侵蚀强度土地面积进行了对比分析,有效评价了矿区土壤侵蚀状况。但是,神东矿区侵蚀营力存在水力、风力及重力侵蚀叠加情况,重力侵蚀主要是井工采煤引起的滑坡和沉陷,如何率定重力侵蚀模数,更准确地掌握矿区水土流失状况,还需进一步研究。