劉翼遙，吳太夏

（河海大學(xué)，南京 211100）

黃土高原是黃河流域的重要組成部分，同時也是我國水土流失問題最嚴(yán)重的地區(qū)。根據(jù)資料顯示，黃土高原總面積為64萬km2，但是截至2017年，其水土流失面積高達(dá)47.2 km2，侵蝕模數(shù)大于8 000 t/km2·a的極強(qiáng)度以上的水蝕面積占全國同類面積的64.95%[1]。黃土高原的土壤流失問題不僅與土質(zhì)松軟、植被覆蓋度過低有關(guān)[2-4]，區(qū)域內(nèi)存在大量坡耕地也是導(dǎo)致水土流失的重要原因[5]。

由于坡耕地的存在，使得黃土高原跑水、跑土、跑肥現(xiàn)象時有發(fā)生[5]。坡耕地與其他耕地類型本質(zhì)的區(qū)別是其具有一定坡度，坡度屬性使得坡耕地水土流失和面源污染問題十分嚴(yán)重。截至2012年，黃土高原坡耕地總面積占耕地總面積的2/3，平均土壤侵蝕模數(shù)2.5萬t/km2·a[6]。治理坡耕地是治理水土流失最重要的步驟之一。1999年我國實(shí)行退耕還林，將易水土流失的坡耕地停耕恢復(fù)植被，黃土高原水土流失得到有效改善，效果顯著[7]。但是2017年黃土高原無定河流域發(fā)生了“7·26”特大暴雨，造成陜西省子洲縣等多個地區(qū)嚴(yán)重內(nèi)澇和土壤侵蝕。子洲縣坡耕地水土流失，形成了許多細(xì)溝和小切溝，嚴(yán)重的土壤侵蝕強(qiáng)度可達(dá)110 000 t/km2，坡耕地多的支溝的侵蝕強(qiáng)度至少為46 000 t/km2，坡耕地少的支溝為12 000 t/km2，前者是后者的3.8倍[8]。本次特大暴雨事件也反映了子洲縣坡耕地面積占比依舊較大，并且分布在陡峭的溝坡，溝蝕嚴(yán)重，造成大量水土流失。因此，準(zhǔn)確掌握坡耕地以及在每個坡度級的面積大小、分布位置，并在相應(yīng)區(qū)域做出防護(hù)措施，就能防止悲劇再次發(fā)生，對減少土壤侵蝕和水土流失也具有重大意義。

然而坡耕地具有較為分散、不夠集中、地形破碎的特點(diǎn)[9]，并且坡耕地普遍分布在較為陡峭的區(qū)域，這無疑給人工調(diào)查增加了難度。而遙感技術(shù)具有大范圍觀測地表的能力，具有全天候、多時相、涵蓋信息多且連續(xù)觀測的特點(diǎn)，并且在耕地等地物信息分類提取等行業(yè)內(nèi)有大量應(yīng)用，已經(jīng)取得良好成果[10-12]。王禎等[6]使用Landsat 4/5結(jié)合Landsat 8數(shù)據(jù)對延安市坡耕地進(jìn)行提取。陳正發(fā)等[9]使用資源環(huán)境數(shù)據(jù)云平臺的土地利用數(shù)據(jù)對云南坡耕地進(jìn)行識別，并對其質(zhì)量進(jìn)行評價。時亞坤等[13]使用GFSAD數(shù)據(jù)提取坡耕地，研究退耕還林對糧食安全的影響。這些坡耕地提取結(jié)果分辨率大多都為30 m，對于那些破碎、分散的坡耕地的提取結(jié)果不高。楊萌等[14]使用無人機(jī)航拍獲取高分辨率影像提取岔巴溝和清水溝小流域的坡耕地。使用無人機(jī)航拍雖然分辨率和精度較高，但是只針對于面積較小的區(qū)域，無法對大區(qū)域做到快速提取。從2015年起，歐空局發(fā)射哨兵（Sentinel）系列衛(wèi)星，其最高分辨率可達(dá)10 m，重訪周期為5 d。哨兵系列衛(wèi)星普遍應(yīng)用于地物的提取，取得了不錯的成果。

本研究融合Sentinel-2遙感影像和ASTER GDEMV2數(shù)據(jù)在子洲縣建立一種大范圍的坡耕地快速提取模型，對坡耕地進(jìn)行大范圍的識別。得到坡耕地的空間分布情況，完成研究區(qū)坡耕地時空特征的分析，為黃土高原依然存在的坡耕地的退耕還林和水土流失防治工作提供參考依據(jù)。

1 研究區(qū)概況與數(shù)據(jù)源

1.1 研究區(qū)概況

子洲縣位于陜西省榆林市南部地區(qū)，全縣總面積達(dá)2 024 km2，位于北緯37°15′30″~37°50′00″，東經(jīng)109°29′08″~110°07′30″之間，其地理位置如圖1所示。子洲縣屬于典型的陜北黃土高原丘陵溝壑區(qū)，地貌主要以梁、峁、溝和川為主，研究區(qū)內(nèi)山區(qū)占總面積的95%，剩下的5%為川區(qū)[15]。子洲縣整體地勢較高，海拔最高為1 339 m，最低為886 m?？傮w來說，子洲縣地貌特征較為統(tǒng)一，便于進(jìn)行大范圍的整體研究。

圖1 研究區(qū)地理位置

子洲縣是農(nóng)業(yè)縣，農(nóng)業(yè)人口占總?cè)丝诘?3%，縣內(nèi)黃土性土壤廣泛分布，占總土地面積的89.97%，宜于發(fā)展農(nóng)業(yè)。子洲縣主要作物為秋糧收獲的玉米、大豆、谷子等，占比可達(dá)全年糧食產(chǎn)量的90.82%。根據(jù)《子洲縣土地利用總體規(guī)劃（2006—2020年）》文件顯示[16]，全縣基本農(nóng)田約為551 km2。并且大部分耕地地理位置所處坡度較高，受到地形地勢影響，耕地較為分散，不易組成統(tǒng)一、集中的區(qū)域，所以需要從縣域角度進(jìn)行坡耕地研究。

研究區(qū)位于溫帶半干旱區(qū)，屬于大陸性季風(fēng)氣候，氣候干燥，日照充足，冷熱溫差較大。子洲縣年均氣溫9.1℃，年平均降雨量為430.8 mm，降水分配極不均勻，5—9月的降雨量占全年總降雨量的90%以上，多集中于幾場高強(qiáng)度、短歷時的暴雨中。子洲縣幾乎每年都會有干旱、霜凍、暴雨、大風(fēng)和沙塵暴災(zāi)害性天氣發(fā)生，土壤侵蝕嚴(yán)重，加之研究區(qū)的坡耕地較多，極易造成水土流失，使得子洲縣成為榆林市水土流失的重點(diǎn)區(qū)域。

1.2 遙感數(shù)據(jù)來源

為滿足子洲縣坡耕地大范圍信息提取的空間分辨率要求（10 m）和研究坡耕地變化的時間分辨率要求，本文選擇Sentinel-2數(shù)據(jù)進(jìn)行研究。Sentinel-2衛(wèi)星是歐空局發(fā)射的雙星衛(wèi)星，Sentinel-2A和Sentinel-2B兩顆衛(wèi)星互為補(bǔ)充，重訪周期可達(dá)到5 d。Sentinel-2衛(wèi)星可覆蓋從可見光波段到短波紅外的13個光譜波段，不同波段的空間分辨率分別為10 m、20 m和60 m。哨兵2號衛(wèi)星為總作物的提取和制圖提供了較為充足且質(zhì)量較好的數(shù)據(jù)源，其具體參數(shù)見表1。

表1 Sentinel-2衛(wèi)星具體參數(shù)

由于研究區(qū)諸如玉米、大豆和谷子等主要作物于每年的9月至10月開始收割，故本文選擇2018年至2021年每年8月前后的四幅影像進(jìn)行耕地的提取。研究區(qū)范圍包含了2張哨兵2號影像，在保證影像質(zhì)量的前提下選擇了2幅S2A影像和2幅S2B影像進(jìn)行實(shí)驗(yàn)，數(shù)據(jù)獲取的信息見表2。

表2 Sentinel-2數(shù)據(jù)獲取信息

本文所采用的Sentinel-2數(shù)據(jù)下載于歐洲航天局（https：//scihub.copernicus.eu/）。接著使用歐空局官方提供的Sen2cor插件對數(shù)據(jù)進(jìn)行大氣校正，最后使用SNAP對校正好的數(shù)據(jù)進(jìn)行拼接、按掩膜裁剪以及格式轉(zhuǎn)換等操作，以供后續(xù)使用。

1.3 DEM數(shù)據(jù)

本文所使用的DEM數(shù)據(jù)是先進(jìn)星載熱發(fā)射和反射輻射儀全球數(shù)字高程模型（Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation）。ASTER GDEM數(shù)據(jù)由日本METI和美國NASA聯(lián)合研制并免費(fèi)面向公眾分發(fā)，是目前唯一覆蓋全球陸地表面的高分辨率高程影像數(shù)據(jù)，并在全球?qū)Φ赜^測研究中取得了廣泛的應(yīng)用。V2版則采用了一種先進(jìn)的算法對V1版GDEM影像進(jìn)行了改進(jìn)，對V1版中存在的錯誤做了很好的矯正，提高了數(shù)據(jù)的空間分辨率精度和高程精度。本次研究所使用的30 m分辨率的ASTER GDEMV2數(shù)據(jù)下載于地理空間數(shù)據(jù)云網(wǎng)站（http：//www.gscloud.cn/home）。

1.4 樣本與驗(yàn)證數(shù)據(jù)

樣本數(shù)據(jù)的好壞直接影響到提取結(jié)果的精度，應(yīng)選擇那些具有代表性的純凈像元進(jìn)行實(shí)驗(yàn)[17]。為驗(yàn)證子洲縣耕地遙感提取的精度，我們實(shí)地采集了耕地及非耕地（包括水體、不透水體、林地和其他）的樣本，之后，再結(jié)合Google Earth高清衛(wèi)星地圖分辨率產(chǎn)品共選擇了206個耕地樣本和428個非耕地樣本，各類樣本均勻分布于整個研究區(qū)，它們的位置如圖1所示。以上兩部分?jǐn)?shù)據(jù)用于耕地提取模型的建立和精度驗(yàn)證。其中，Google Earth高清衛(wèi)星地圖分辨率產(chǎn)品為17級，其空間分辨率為2.39 m。

2 研究方法

2.1 地物分類體系

本次研究參照GB/T 21010—2017《土地利用現(xiàn)狀分類》的12個一級分類、73個二級分類規(guī)定對研究區(qū)地物進(jìn)行分類。由于使用的衛(wèi)星產(chǎn)品分辨率有所不同，在真彩色圖像下，同一處耕地的情況存在著區(qū)別，展示了不同的耕地特征。因此在研究中，將耕地具體分為未種植作物、種植作物兩種。研究區(qū)的地物類型劃分出種植作物耕地、未種植作物耕地、不透水體、水體、林地及其他幾類，各種地物類型的衛(wèi)星影像圖如圖2所示。

2.2 耕地提取方法

本文采用基于規(guī)則的面向?qū)ο筇崛》椒?，建立耕地提取決策樹。面向?qū)ο蟮臎Q策樹遙感分類方法首先對影像進(jìn)行分割，合并，將很多對象合并為一個整體；然后針對每一個對象構(gòu)建分類規(guī)則，建立分類決策樹，從而實(shí)現(xiàn)對地物的分類提取[18]。

本次研究以分割閾值45.0和合并閾值85.0進(jìn)行實(shí)驗(yàn)，再結(jié)合由光譜特征、紋理特征、幾何特征、歸一化差異植被指數(shù)及水體指數(shù)分析得到的劃分標(biāo)準(zhǔn)，進(jìn)行適宜的耕地提取規(guī)則設(shè)定，構(gòu)建子洲縣耕地提取決策樹，接著再根據(jù)DEM影像進(jìn)行坡度分級，達(dá)到提取研究區(qū)坡耕地的結(jié)果。歸一化差異水體指數(shù)（NDWI）可以實(shí)現(xiàn)水體與陸地的分離[19]。與其他水體指數(shù)相比，NDWI能有效抑制植被信號，增強(qiáng)水分信號，其公式如公式（1）所示。歸一化差異植被指數(shù)（NDVI）常常被用于對農(nóng)作物進(jìn)行分類，以及對半干旱地區(qū)降水量的預(yù)測研究，占有相當(dāng)重要的位置[20-21]，其公式如公式（2）所示。

式中：NIR為影像的近紅外波段；R為影像紅外波段；G為綠波段。

利用上述指數(shù)區(qū)別出水體和有植被地，再通過紋理特征和光譜特征，可以進(jìn)一步區(qū)分出種植作物耕地和未種植作物耕地，決策樹如圖3所示。得到耕地所在位置后，根據(jù)DEM影像進(jìn)行坡度分級，與提取出來的耕地進(jìn)行相交處理，進(jìn)而得到坡耕地的分布位置。

圖3 耕地提取決策樹

圖3所示閾值中，NDWI為歸一化水體指數(shù)，NDVI為歸一化植被指數(shù)，Spectral Mean（Bi）為第i波段的光譜均值，Bi_ME為第i波段的紋理均值，Retangular Fit為矩形化擬合程度。

2.3 精度評價

在地物提取工作結(jié)束后，必須對得到的提取結(jié)果進(jìn)行精度評價，驗(yàn)證其是否有令人滿意的準(zhǔn)確性。學(xué)者常用混淆矩陣和Kappa系數(shù)進(jìn)行精度評價工作。矩陣行與列分別表示分類得到結(jié)果與實(shí)際地物情況，對角線上數(shù)據(jù)表示得到正確分類的數(shù)量，利用對角線數(shù)據(jù)、行總值、列總值等數(shù)據(jù)，可以得到總體精度、Kappa系數(shù)等評價指標(biāo)，混淆矩陣見表3，總體精度、Kappa系數(shù)計算公式如公式（3）、公式（4）所示。

表3 混淆矩陣

式中：N為樣本的總數(shù)，n為一共有多少種分類，Xii為對角線數(shù)據(jù)，Xi+與X+i分別為某種分類的列/行總值。Kappa系數(shù)的范圍為0~1，得到系數(shù)越大，表明分類一致性越高，精度相較而言更好。當(dāng)Kappa系數(shù)小于0.2，表示該分類結(jié)果較差；當(dāng)Kappa系數(shù)大于0.2小于0.6時，分類一致性一般；當(dāng)Kappa系數(shù)大于0.6小于0.8時，表示分類精度較高；當(dāng)Kappa系數(shù)大于0.8時，該情況下可以近似于分類結(jié)果與實(shí)際趨于一致。

3 結(jié)果與分析

3.1 精度評價結(jié)果

按上述方法對子洲縣2018年至2021年耕地進(jìn)行提取，并計算提取結(jié)果的精度。精度評價結(jié)果如圖4和表4所示。

表4 耕地提取結(jié)果精度評價

圖4 子洲縣2018年至2021年耕地提取結(jié)果精度評價

從圖表中的結(jié)果可知，2018年至2021年耕地提取精度都較好，提取結(jié)果的總體精度都大于80%，并且每年的Kappa系數(shù)都大于0.6，表明提取的耕地與實(shí)際較為一致，結(jié)果可結(jié)合DEM數(shù)據(jù)進(jìn)行坡度分級，提取子洲縣2018年至2021年坡耕地。

3.2 坡耕地提取結(jié)果

以1984年《土地利用現(xiàn)狀調(diào)查技術(shù)規(guī)程》為技術(shù)依據(jù)，可將耕地以坡度為標(biāo)準(zhǔn)，劃分為5種等級，分別為小于等于2°；2~6°；6~15°；15~25°；大于25°進(jìn)行坡度分級。目前多少坡度范圍的耕地屬于坡耕地，不同研究人員給出的標(biāo)準(zhǔn)也不盡相同[22-24]。本次研究根據(jù)羅光杰等[25]的研究，將坡度大于6°的耕地認(rèn)為坡耕地，進(jìn)行坡耕地的提取，其結(jié)果如圖5所示。

圖5 2018年至2021年子洲縣坡耕地提取結(jié)果

使用ArcGIS的相交工具計算4年中每個坡度級上坡耕地的面積和占比，其結(jié)果如表5、圖6所示。

從表5和圖6結(jié)果可知，子洲縣在2018年至2020年間，耕地分布情況變化并不明顯，基本都維持在520~530 km2左右，仍然維持原來的高程、坡度等級進(jìn)行耕種活動。2020年至2021年，坡耕地面積從519.02 km2增加到569.48 km2，變化較為明顯，總體呈增加趨勢。但從每個坡度級的坡耕地面積占當(dāng)年坡耕地總面積來說，變化較不明顯，4年都基本維持在同一數(shù)量級上。從每個坡度級占比來看，4年中6~15°坡度級的占比最大，都大于50%，最高可達(dá)57.28%，其次是15~25°和大于25°坡度級。該種現(xiàn)象與研究區(qū)溝壑起伏變化大等情況有關(guān)，該因素也直接影響子洲縣耕地分布總體較為零散，以未來發(fā)展角度，不宜進(jìn)行更高效率的大機(jī)器耕作方法，對于農(nóng)業(yè)縣來說，對其發(fā)展并不友善。根據(jù)陳正發(fā)[24]的研究，大于25°的坡耕地極易發(fā)生水土流失，應(yīng)當(dāng)嚴(yán)格實(shí)行退耕還林。但是這4年中子洲縣大于25°的坡耕地依然存在，并且面積都大于30 km2，在2021年甚至達(dá)到了44 km2，需要當(dāng)?shù)卣谄胶飧丶t線的情況下，更好更快落實(shí)退耕還林、還草的相關(guān)政策。

表5 2018年至2021年每個坡度級中坡耕地面積占比

圖6 2018年至2021年每個坡度級中耕地面積

4 結(jié)論

本研究以陜西省榆林市子洲縣作為研究區(qū)域，利用Sentinel-2和ASTER GDEMV2作為數(shù)據(jù)來源，基于規(guī)則的面向?qū)ο蠓诸惙椒?，確定影像分割、影像合并閾值大小，以及各種地物類型的提取規(guī)則，并最終建立關(guān)于種植耕地、未種植耕地及非耕地的提取決策樹，基于實(shí)地勘察數(shù)據(jù)和googleearth影像，對提取結(jié)果進(jìn)行驗(yàn)證，得到Kappa系數(shù)高于0.65的分類結(jié)果，可認(rèn)為一致性較高，精度較好。接著對2018至2021年，子洲縣同期耕地時空變化情況進(jìn)行分析，發(fā)現(xiàn)子洲縣在2018年至2020年間，耕地分布情況變化并不明顯。2021年，坡耕地面積增加到569.48 km2，變化較為明顯。4年中6~15°坡度級的占比最大，都大于50%，最高可達(dá)57.28%，其次是15~25°和大于25°坡度級。并且這4年中子洲縣大于25°的坡耕地依然存在，并且面積都大于30 km2，在2021年甚至達(dá)到了44 km2。