王麗云，王小燕，王 義，岳本江，屈 創(chuàng)

(1.黃河上中游管理局，陜西 西安 710021； 2.神華神東煤炭集團(tuán)有限公司，陜西 神木 719300)

神東礦區(qū)位于晉陜蒙接壤地區(qū)，是我國最重要的煤炭生產(chǎn)和能源基地之一[1]。面對(duì)大規(guī)模開采與極脆弱生態(tài)環(huán)境保護(hù)之間的矛盾，礦區(qū)不斷創(chuàng)新生產(chǎn)方式，探索生態(tài)治理與環(huán)境保護(hù)的新理念，提出了以控制外圍風(fēng)沙侵蝕為前提，內(nèi)外圍結(jié)合治理，促進(jìn)礦區(qū)整體恢復(fù)和改善礦區(qū)生態(tài)環(huán)境的思路，走出了一條資源保護(hù)型開采與生態(tài)環(huán)境治理相協(xié)調(diào)的綠色礦業(yè)之路[2]。為全面體現(xiàn)神東礦區(qū)2005—2018年水土流失綜合治理成果，客觀評(píng)價(jià)礦區(qū)水土保持工作成效，本研究采用中國土壤流失方程(CSLE)和風(fēng)力侵蝕模型，定量計(jì)算礦區(qū)2005、2010、2015、2018年4期土壤侵蝕模數(shù)，并對(duì)不同時(shí)期不同侵蝕強(qiáng)度土地面積進(jìn)行對(duì)比分析，希望能為礦區(qū)水土流失綜合治理提供數(shù)據(jù)支撐。

1 研究區(qū)概況

神東礦區(qū)分布于陜西省榆林市神木市北部、府谷縣西部，內(nèi)蒙古自治區(qū)鄂爾多斯市東勝區(qū)及伊金霍洛旗的南部、準(zhǔn)格爾旗的西南部，地理位置為109°51′～110°46′E、38°52′～39°41′N，包括大柳塔礦、活雞兔礦、哈拉溝礦等14個(gè)礦井，總面積約1 064 km2，黃河一級(jí)支流窟野河(烏蘭木倫河)縱貫礦區(qū)中部。地處黃土高原北部與毛烏素沙地交錯(cuò)區(qū)，具有典型的干旱荒漠高原氣候特征，多年平均降水量?jī)H362 mm，年水面蒸發(fā)量為2 297.4～2 838.7 mm[1]。區(qū)內(nèi)大部分為典型的風(fēng)成沙丘及沙漠灘地，尤其是北部、東南部為黃土丘陵溝壑區(qū)，梁峁起伏，溝壑縱橫，地表支離破碎，風(fēng)水蝕交互作用導(dǎo)致水土流失嚴(yán)重。

2 監(jiān)測(cè)方法

2.1 數(shù)據(jù)來源與預(yù)處理

本研究主要涉及遙感影像、地形數(shù)據(jù)、氣象資料、土壤資料等基礎(chǔ)數(shù)據(jù)。

(1)遙感影像。采用2005、2010年的SPOT5影像和2015、2018年的ZY-3數(shù)據(jù)進(jìn)行2005、2010、2015、2018年礦區(qū)土地利用和水土保持措施遙感解譯。采用監(jiān)測(cè)年前3年空間分辨率為250 m的MOD13Q1數(shù)據(jù)計(jì)算植被覆蓋度。采用監(jiān)測(cè)年逐日AMSR-E Level 2A亮溫?cái)?shù)據(jù)產(chǎn)品(空間分辨率為25 km)計(jì)算表土濕度數(shù)據(jù)。遙感影像使用前進(jìn)行輻射校正、幾何校正、配準(zhǔn)、影像融合、鑲嵌與裁切等預(yù)處理，融合后的影像既提高了多光譜影像的空間分辨率，又保留了其多光譜特性。遙感數(shù)據(jù)統(tǒng)一采用CGCS2000坐標(biāo)系、Albers投影。

(2)地形數(shù)據(jù)?；诘V區(qū)1∶10 000比例尺的數(shù)字高程模型DEM，用符素華等[3]開發(fā)的土壤侵蝕模型地形因子計(jì)算工具計(jì)算坡長(zhǎng)、坡度因子等。

(3)其他數(shù)據(jù)。收集礦區(qū)8個(gè)站點(diǎn)2000年以來的逐日降水量資料和逐日整點(diǎn)風(fēng)速資料，用于計(jì)算降雨侵蝕力因子和風(fēng)力因子。植被覆蓋度修訂系數(shù)采用第一次全國水利普查水土保持情況普查成果數(shù)據(jù)。

2.2 土壤侵蝕模型

2.2.1 水力侵蝕模型

在水力侵蝕區(qū)，采用中國土壤流失方程(CSLE)計(jì)算土壤侵蝕模數(shù)[4-5]。計(jì)算公式為

A=R·K·L·S·B·E·T

(1)

式中：A為土壤侵蝕模數(shù)，t/(hm2·a)；R為降雨侵蝕力因子，MJ·mm/(hm2· h·a)；K為土壤可蝕性因子，t·hm2·h/(hm2·MJ·mm)；L為坡長(zhǎng)因子；S為坡度因子；B為生物措施因子；E為工程措施因子；T為耕作措施因子。

(1)降雨侵蝕力因子。將8個(gè)站點(diǎn)多年平均1～24個(gè)半月降雨侵蝕力矢量化，利用普通克里金空間插值方法，生成空間分辨率為2 m的24個(gè)半月降雨侵蝕力柵格數(shù)據(jù)。

(2)土壤可蝕性因子。采用第一次全國水利普查水土保持情況普查土壤可蝕性因子成果作為礦區(qū)土壤可蝕性因子值。

(3)坡長(zhǎng)、坡度因子?；趨^(qū)域DEM數(shù)據(jù)，采用北京師范大學(xué)符素華等[3]開發(fā)的地形因子計(jì)算工具計(jì)算坡長(zhǎng)、坡度因子值，并對(duì)林草地采用緩坡公式修正。

(4)生物措施因子?；贛OD13Q1 NDVI數(shù)據(jù)采用像元二分法計(jì)算植被覆蓋度，利用植被覆蓋度修訂系數(shù)，獲得每年24個(gè)半月30 m空間分辨率的植被覆蓋度，再結(jié)合24個(gè)半月降雨侵蝕力因子比例與土地利用類型計(jì)算園地、林地、草地的生物措施因子；參照《區(qū)域水土流失動(dòng)態(tài)監(jiān)測(cè)技術(shù)規(guī)定(試行)》對(duì)非園地、林地和草地生物措施因子查表賦值。

(5)工程措施因子。根據(jù)野外調(diào)查和遙感影像解譯，礦區(qū)水土保持工程措施主要為土坎水平梯田，工程措施因子賦值為0.084。

(6)耕作措施因子。礦區(qū)在全國輪作區(qū)中涉及兩個(gè)一級(jí)區(qū)、三個(gè)二級(jí)區(qū)。一級(jí)區(qū)北部中高原半干旱喜涼作物一熟區(qū)下的二級(jí)區(qū)后山壩上晉北高原山地半干旱喜涼作物一熟區(qū)和隴中青東寧中南黃土丘陵半干旱喜涼作物一熟區(qū)，耕作措施因子取值0.488；一級(jí)區(qū)北部低高原易旱喜溫一熟區(qū)下的二級(jí)區(qū)黃土高原東部易旱喜溫作物一熟區(qū)，耕作措施因子取值0.417。

2.2.2 風(fēng)力侵蝕模型

在風(fēng)力侵蝕地區(qū)，根據(jù)土地利用類型，選用耕地、草(灌)地、沙地(漠)等風(fēng)力侵蝕模型，計(jì)算風(fēng)力侵蝕模數(shù)[6-7]。耕地、草(灌)地、沙地(漠)風(fēng)力侵蝕模型計(jì)算公式分別為

(2)

(3)

(4)

式中：Qfa為每半個(gè)月內(nèi)耕地風(fēng)力侵蝕模數(shù)，t/(hm2·a)；W為每半個(gè)月內(nèi)表土濕度因子，介于0～1之間；Tj為每半個(gè)月內(nèi)各風(fēng)速等級(jí)的累計(jì)時(shí)間，min；Z0為地表粗糙度，cm；j為風(fēng)速等級(jí)序號(hào)；Uj為第j等級(jí)的平均風(fēng)速，m/s；Qfg為每半個(gè)月內(nèi)草(灌)地風(fēng)力侵蝕模數(shù)，t/(hm2·a)；V為植被覆蓋度，%；Qfs為每半個(gè)月內(nèi)沙地(漠)風(fēng)力侵蝕模數(shù)，t/(hm2·a)。

(1)風(fēng)力因子。風(fēng)力因子是指不同等級(jí)風(fēng)速和風(fēng)向?qū)︼L(fēng)力侵蝕的潛在能力，與臨界侵蝕風(fēng)速和各等級(jí)風(fēng)速的累積時(shí)間有關(guān)。風(fēng)力因子計(jì)算采用收集的礦區(qū)逐日24 h整點(diǎn)風(fēng)速和風(fēng)向資料，按1 m/s間隔統(tǒng)計(jì)全年大于等于臨界侵蝕風(fēng)速的各等級(jí)風(fēng)速的累積時(shí)間，利用普通克里金方法，插值生成24個(gè)半月各風(fēng)速等級(jí)對(duì)應(yīng)的累計(jì)時(shí)間。耕地的臨界侵蝕風(fēng)速取值為5 m/s，草(灌)地、沙地(漠)的臨界侵蝕風(fēng)速按照李智廣等[7]在我國風(fēng)力侵蝕抽樣調(diào)查方法研究中不同植被覆蓋度的臨界侵蝕風(fēng)速取值。

(2)表土濕度因子。表土濕度因子表征了土壤表層0～2.5 cm深度范圍內(nèi)含水率抗拒土壤風(fēng)力侵蝕的潛在能力。基于AMSR-E Level 2A亮溫?cái)?shù)據(jù)計(jì)算每天的表土濕度因子。

(3)植被覆蓋度。計(jì)算方法同水力侵蝕植被覆蓋度。

(4)地表粗糙度因子。地表粗糙度因子指因植被、微地形和農(nóng)田耕作導(dǎo)致的零風(fēng)速位置的高度。在風(fēng)力侵蝕區(qū)，除耕地因呈斑塊狀分布，地表粗糙度在邊界發(fā)生突變外，其他諸如草(灌)地和沙地(漠)等的邊界地表粗糙度都呈逐漸過渡特征。根據(jù)不同土地利用類型，采用李智廣等[7]在我國風(fēng)力侵蝕抽樣調(diào)查方法研究中不同下墊面粗糙度的賦值。

3 結(jié)果與分析

3.1 水土流失狀況

神東礦區(qū)各時(shí)期風(fēng)力侵蝕和水力侵蝕面積見表1。2005年礦區(qū)平均土壤侵蝕模數(shù)為3 322.53 t/(km2·a)，水土流失面積943.36 km2，占礦區(qū)總面積的88.66%，侵蝕強(qiáng)度以輕、中度為主。2010年礦區(qū)平均土壤侵蝕模數(shù)為2 497.58 t/(km2·a)，水土流失面積855.47 km2，占礦區(qū)總面積的80.40%，侵蝕強(qiáng)度以輕、中度為主。2015年礦區(qū)平均土壤侵蝕模數(shù)為1 843.98 t/(km2·a)，水土流失面積730.66 km2，占礦區(qū)總面積的68.67%，侵蝕強(qiáng)度以輕度為主。2018年礦區(qū)平均土壤侵蝕模數(shù)為852.17 t/(km2·a)，水土流失面積587.61 km2，占礦區(qū)總面積的55.23%，侵蝕強(qiáng)度以輕度為主。

表1 神東礦區(qū)各時(shí)期不同侵蝕強(qiáng)度土地面積統(tǒng)計(jì)

3.2 水土流失動(dòng)態(tài)變化

對(duì)比4期土壤侵蝕數(shù)據(jù)(表2)，2005—2010年神東礦區(qū)水土流失面積減少了87.89 km2，中度及以上侵蝕面積均有所減少，其中：中度侵蝕面積減少最明顯，減少了89.75 km2；其次為強(qiáng)烈侵蝕，面積減少了64.76 km2。2010—2015年水土流失面積減少了124.81 km2，中度、強(qiáng)烈侵蝕面積減少明顯，分別減少了117.17、72.41 km2。2015—2018年水土流失面積減少了143.05 km2，中度侵蝕面積減少最多，減少了139.97 km2。

表2 神東礦區(qū)各時(shí)期不同侵蝕強(qiáng)度土地面積變化統(tǒng)計(jì)

總的來看，2005—2018年神東礦區(qū)水土流失面積明顯減少，共減少了355.75 km2，中度侵蝕面積減少最多，其次為強(qiáng)烈侵蝕，分別減少了346.89、170.96 km2。對(duì)比分析礦區(qū)不同時(shí)期水土流失面積動(dòng)態(tài)變化可知，2005—2018年礦區(qū)水土流失狀況為高強(qiáng)度侵蝕向低強(qiáng)度侵蝕轉(zhuǎn)移，水土流失面積明顯減少，水土流失強(qiáng)度明顯降低，生態(tài)環(huán)境顯著改善。

4 結(jié) 語

2005—2018年，神東礦區(qū)土壤侵蝕模數(shù)由3 322.53 t/(km2·a)降低至852.17 t/(km2·a)，水土流失面積減少了355.75 km2，相當(dāng)于礦區(qū)總面積的1/3還多。取得如此成效主要得益于礦區(qū)始終堅(jiān)持“生態(tài)礦區(qū)、環(huán)保礦井、清潔煤炭”的建設(shè)理念。礦區(qū)不僅實(shí)施了傳統(tǒng)水土保持工程措施及生態(tài)水保林、經(jīng)濟(jì)林、灌草等林草植被措施，還創(chuàng)造性地利用柳桿障蔽生態(tài)錨固坡等技術(shù)防治水土流失。截至2018年底，礦區(qū)已累計(jì)治理水土流失面積339 km2，水土流失和土地沙化得到有效控制，植被覆蓋度大幅度提升，由開發(fā)初的3%～11%提高到60%以上[8]，生態(tài)環(huán)境得到顯著改善。

本研究采用資料收集、遙感監(jiān)測(cè)、模型計(jì)算(風(fēng)力侵蝕、水力侵蝕模型)等方法，定量推算了神東礦區(qū)2005、2010、2015、2018年4期土壤侵蝕模數(shù)，并對(duì)不同侵蝕強(qiáng)度土地面積進(jìn)行了對(duì)比分析，有效評(píng)價(jià)了礦區(qū)土壤侵蝕狀況。但是，神東礦區(qū)侵蝕營(yíng)力存在水力、風(fēng)力及重力侵蝕疊加情況，重力侵蝕主要是井工采煤引起的滑坡和沉陷，如何率定重力侵蝕模數(shù)，更準(zhǔn)確地掌握礦區(qū)水土流失狀況，還需進(jìn)一步研究。