基于WorldView-2影像估算黃土高原六道溝小流域退耕還林(草)后的切溝發(fā)育速率

黃婷婷, 楊 揚?, 史揚子, 曹 琦, 馮 博, 劉寶元, 劉瑛娜

(1. 北京師范大學(xué)地理科學(xué)學(xué)部 地表過程與資源生態(tài)國家重點實驗室, 100875, 北京;2. 北京師范大學(xué)地理科學(xué)學(xué)部地理學(xué)院, 100875, 北京; 3. 中國科學(xué)院 水利部 水土保持研究所, 712100, 陜西楊凌; 4. 沈陽農(nóng)業(yè)大學(xué)水利學(xué)院 遼寧省水土流失防控與生態(tài)修復(fù)重點實驗室, 110866, 沈陽)

切溝是黃土高原主要的土壤侵蝕類型之一。在黃土丘陵溝壑區(qū)，切溝侵蝕的產(chǎn)沙量通常占小流域總產(chǎn)沙量的50%以上，部分地區(qū)甚至可以達(dá)到70%以上[1-2]。21世紀(jì)以來，GPS技術(shù)發(fā)展迅速，其中實時動態(tài)差分GPS(real-time kinematic GPS, RTK GPS)因其定位精度高、作業(yè)速度快等優(yōu)勢，成為切溝監(jiān)測的重要手段之一[3]；然而，RTK GPS測量成本高、野外負(fù)擔(dān)重，多用于較小時空尺度的監(jiān)測研究。近年來，高分辨率遙感產(chǎn)品不斷涌現(xiàn)，兼具監(jiān)測時間長、覆蓋范圍大等特點，為研究較大時空尺度的切溝分布特征與發(fā)育速率提供便捷可靠的數(shù)據(jù)源[4-5]。利用高精度野外測量儀器如三維激光掃描儀[6]、RTK GPS[7]評估遙感影像的切溝解譯精度，分析誤差來源并總結(jié)解譯經(jīng)驗，對于提高遙感影像解譯精度具有重要意義。

切溝發(fā)育與集水區(qū)面積、降雨強度、土地利用和土壤物理性質(zhì)等密切相關(guān)[8-9]。黃土高原切溝侵蝕嚴(yán)重，退耕還林(草)工程大規(guī)模實施后，土地利用格局改變，必然導(dǎo)致切溝侵蝕的相應(yīng)變化。筆者以六道溝小流域為研究對象，利用RTK GPS進行切溝實地測量，評價WorldView-2影像的切溝解譯精度；在此基礎(chǔ)上，對比2012年和2018年WorldView-2影像的切溝解譯結(jié)果，估算退耕后的切溝發(fā)育速率并分析其影響因素，以期為治溝治坡決策、水土保持措施配置提供科學(xué)依據(jù)。

1 研究區(qū)概況

六道溝小流域(E 110°21′～110°23′，N 38°46′～38°51′)位于陜西省神木市以西14 km(圖1)，屬黃土高原北部水蝕風(fēng)蝕交錯區(qū)，是典型的蓋沙黃土丘陵溝壑地貌。小流域面積689.2 hm2，主溝道長4.2 km，自南向北流入窟野河的一級支流三道河。小流域?qū)僦袦貛О敫珊禋夂颍那锛竟?jié)暴雨集中，冬春季節(jié)干旱多風(fēng)。自2002年退耕還林(草)工程在神木市全面實施以來，小流域內(nèi)大量坡耕地退耕；同時，林草生態(tài)恢復(fù)與人工造林種草等生物措施廣泛開展，植被蓋度增加[10]。

圖1 六道溝小流域與RTK GPS所測切溝的地理位置Fig.1 Geographic locations of Liudaogou catchment and gullies measured by RTK GPS

經(jīng)實地踏查，筆者在六道溝小流域共選取33條切溝開展RTK GPS測量(圖1b)。2018年10月，在道路兩側(cè)及步行易達(dá)的山坡測量10條代表性切溝，切溝長度介于8～71 m之間，面積范圍25～613 m2，主要土地利用類型為草地。2019年5月，在這10條切溝的基礎(chǔ)上，根據(jù)溝緣植被、溝壁特征等另外選取23條切溝，長度范圍擴大至7～200 m，面積20～2 500 m2。

2 材料與方法

2.1 WorldView-2影像的切溝解譯驗證

2.1.1 RTK GPS測量 采用水平測量精度<1 cm、垂直測量精度<2 cm的中海達(dá)iRTK2經(jīng)典版差分GPS對所選切溝開展實地測量。按大致0.5 m的間距測量切溝溝緣線各點的地理坐標(biāo)，同時在溝緣曲折變化劇烈處適當(dāng)增加測量點。將所測各點坐標(biāo)導(dǎo)入ArcGIS 10.2，分別提取各切溝的長度、最大表面寬度、周長和面積。

2.1.2 WorldView-2影像處理與解譯 以2018年10月23日的WorldView-2(WV2)影像為數(shù)據(jù)源，對0.5 m分辨率的全色數(shù)據(jù)和2.4 m分辨率的多光譜數(shù)據(jù)進行正射校正后，利用Gram-schmidt方法融合校正后的全色和多光譜數(shù)據(jù)，得到0.5 m分辨率的真彩色影像。在此基礎(chǔ)上，根據(jù)形狀、尺寸、陰影、色調(diào)和紋理等對影像進行切溝目視解譯。

2.1.3 WorldView-2影像解譯精度評估 為評估WV2解譯各切溝參數(shù)p的精度，計算其相對誤差

(1)

式中：pwv和pgps分別為基于WV2和RTK GPS得到的切溝參數(shù)p。

在此基礎(chǔ)上，采用配對制圖誤差指數(shù)E分析WV2和RTK GPS描繪切溝輪廓的差異[11]：

(2)

式中：Awv∪Agps為基于WV2和RTK GPS得到的切溝并集面積，hm2，Awv∩Agps為相應(yīng)的切溝交集，hm2；E取值范圍為0～100%，0表示二者輪廓完全吻合，100%表示完全不吻合。

2.2 小流域集水區(qū)劃分及切溝解譯

2.2.1 集水區(qū)劃分及其屬性提取 基于1∶1萬地形圖生成5 m分辨率的DEM(digital elevation model)，將六道溝小流域劃分為若干集水區(qū)，并參照Strahler[12]河網(wǎng)分級原則進行集水區(qū)分級，即單條干溝控制的集水區(qū)為Ⅰ級集水區(qū)；2個Ⅰ級集水區(qū)交匯形成Ⅱ級集水區(qū)；不同級別的集水區(qū)交匯，等級為原來集水區(qū)中級別較高者。

利用DEM提取各集水區(qū)的平均坡度、坡長、坡向、剖面曲率和平面曲率作為背景因子。因2012年無可用Landsat數(shù)據(jù)(http://www.gscloud.cn)，以2013—2018年夏季30 m 分辨率的Landsat 8 OLI影像為數(shù)據(jù)源提取六道溝小流域的植被蓋度,計算小流域多年平均植被蓋度后，再根據(jù)各集水區(qū)范圍得到其對應(yīng)植被蓋度。

參照小流域2010年和2017年的遙感影像和土地利用數(shù)據(jù)[10]，建立遙感解譯標(biāo)志，分別對2012年和2018年的WV2影像進行判讀，獲得相應(yīng)的土地利用數(shù)據(jù)。由此計算各集水區(qū)的綜合土地利用動態(tài)度[13]，以量化集水區(qū)d的土地利用在2012—2018年間(Δt)的總體變化速率：

(3)

式中：Kd為集水區(qū)d的綜合土地利用動態(tài)度，%/a；Δt為研究期時間，a；ΔUij為集水區(qū)d在Δt內(nèi)土地利用類型i轉(zhuǎn)變?yōu)橥恋乩妙愋蚸的面積，hm2；n為集水區(qū)d的土地利用類型總數(shù)；Sd為集水區(qū)d的面積，hm2。

2.2.2 集水區(qū)切溝解譯 以2012年11月23日與2018年10月23日的WV2遙感影像為數(shù)據(jù)源，分別解譯各集水區(qū)內(nèi)的切溝，影像處理及解譯方法同2.1.2。在此基礎(chǔ)上，統(tǒng)計各集水區(qū)各切溝的長度、最大表面寬度、周長和面積，分別加和得到對應(yīng)集水區(qū)的切溝總長度、總最大表面寬度、總周長和總面積。

2012—2018年間(Δt)集水區(qū)d的參數(shù)p年均變化率

(4)

式中：Δd,p為集水區(qū)d在Δt內(nèi)的參數(shù)p年均變化率，%/a；pd,wv2018和pd,wv2012分別為2018年和2012年得到的集水區(qū)d參數(shù)p的總和。

3 結(jié)果與分析

3.1 WorldView-2影像的切溝解譯精度評價

3.1.1 WorldView-2與RTK GPS所得切溝參數(shù)對比 相比2018年10月RTK GPS測量的10條切溝數(shù)據(jù)，同時期WV2解譯得到的切溝長度、最大表面寬度、周長和面積相對誤差分別介于0.02%～3.2%、0.1%～9.2%、1.6%～11.7%和0.1%～9.3%之間，平均值分別為1.3%、3.8%、5.4%和3.5%。配對制圖誤差指數(shù)介于7.4%～24.4%之間，平均15.9%，小于Fiorucci等[14]以RTK GPS為基準(zhǔn)得到的18%的WV2制圖誤差指數(shù)，表明同時相WV2影像與RTK GPS得到的切溝輪廓非常接近，可用于切溝解譯。

對比2018年10月和2019年5月利用RTK GPS測得的10條切溝數(shù)據(jù)發(fā)現(xiàn)，2個時期切溝變化不大，切溝長度、最大表面寬度、周長和面積變化率分別介于0.004%～3.5%、0.1%～3.4%、0.7%～3.1%和0.5%～6.7%之間，平均值依次為1.3%、1.4%、1.9%和2.5%。盡管小流域位于黃土高原水蝕風(fēng)蝕交錯區(qū)，4—5月風(fēng)蝕強盛，7—8月水蝕嚴(yán)重，但土壤侵蝕仍以水蝕為主，占全年土壤侵蝕總量的80%以上[15]。此外，切溝是水力侵蝕的結(jié)果，風(fēng)蝕對坡面切溝形態(tài)的影響較小，主要表現(xiàn)為對切溝溝岸的吹蝕與磨蝕。因此，筆者認(rèn)為2018年10月至2019年5月期間切溝發(fā)育不明顯。由于2019年5月所測切溝數(shù)量多、規(guī)模代表性強，以該時期的RTK GPS結(jié)果為基準(zhǔn)，驗證2018年10月WV2影像的切溝解譯精度。

根據(jù)配對t檢驗，2種方法所測切溝長度與最大表面寬度無顯著差異，周長和面積則分別在0.01和0.05的置信水平上差異顯著。統(tǒng)計切溝形態(tài)參數(shù)的相對誤差發(fā)現(xiàn)，切溝長度、最大表面寬度、周長和面積的相對誤差范圍分別為0.1%～10.6%、0.2%～16.7%、0.3%～14.4%和0.2%～17.4%，平均1.9%、2.9%、6.0%和3.8%；且誤差主要集中在5%以內(nèi)，分別占所測切溝總條數(shù)的94%、79%、58%和76%。配對制圖誤差指數(shù)范圍為5.3%～39.6%，平均17.3%，同樣<18%[14]，說明WV2影像解譯可用于六道溝小流域的切溝研究。

盡管WV2與RTK GPS得到的切溝周長與面積差異顯著，但對應(yīng)散點的分布與1∶1線十分接近(圖2)。為準(zhǔn)確描述這2個形態(tài)參數(shù)，以RTK GPS結(jié)果為基準(zhǔn)，分別建立WV2所得切溝周長與面積的轉(zhuǎn)換方程(圖2)。結(jié)果表明，2個方程均顯著，且決定系數(shù)都十分接近1，可用于對應(yīng)參數(shù)的轉(zhuǎn)換。

圖2 WV2影像與RTK GPS提取的切溝周長和面積及對應(yīng)轉(zhuǎn)換方程Fig.2 Linear conversions of gully perimeter and area extracted by WV2 image and RTK GPS

3.1.2 WorldView-2切溝解譯誤差的影響因素 對比分析WV2與RTK GPS提取的各切溝輪廓，WV2的切溝解譯誤差可能主要與溝頭活躍與否、溝緣植被類型及溝壁土體狀況有關(guān)。切溝溝頭可分為不活躍和活躍2種類型[16]。不活躍溝頭具圓角溝緣，無陷穴，溝頭邊壁傾斜，溝壁及溝底有植被覆蓋；活躍溝頭表現(xiàn)為銳利溝緣，有陷穴，溝底有下切、沉積和流水痕跡。相比之下，具不活躍溝頭的切溝溝緣在影像上更不易識別(圖3a)。測量的33條切溝中，此類切溝共5條， WV2提取的各切溝參數(shù)平均誤差均在5%左右，配對制圖誤差指數(shù)平均為22.7%，切溝輪廓存在較大差異。具活躍溝頭的切溝共28條，WV2與RTK GPS得到的輪廓較吻合(圖3b)，配對制圖誤差指數(shù)平均僅16.3%；且除面積平均誤差略高于5%外，其余形態(tài)參數(shù)的誤差皆<3.5%。

植被可有效減緩切溝發(fā)育，但當(dāng)溝緣被樹冠遮蔽時，會給影像解譯帶來極大的困難[17]。六道溝小流域中的切溝溝緣多有草本和灌木植被覆蓋，溝底則為喬灌草混合分布。當(dāng)溝緣被灌木覆蓋時，從影像上難以準(zhǔn)確識別(圖3c)。對9條溝緣以灌木為主的切溝而言，WV2解譯的最大表面寬度和周長都與RTK GPS的對應(yīng)結(jié)果之間存在顯著差異(P<0.05)，各形態(tài)參數(shù)的平均誤差介于2.0%～6.4%之間。余下24條溝緣被草本植物覆蓋的切溝的解譯精度則更高(圖3d)，各參數(shù)的平均誤差與配對制圖誤差指數(shù)均更小。一般來說，溝底的草本、灌木并不會影響切溝溝緣的識別，但溝底喬木在影像上多近似圓形(圖3d)，當(dāng)靠近溝緣時易與植被覆蓋的溝緣混淆，需仔細(xì)分辨。

實地調(diào)查發(fā)現(xiàn)，切溝溝壁因植被稀少或溝岸塌陷出現(xiàn)的土體裸露往往在影像上呈白色或灰白色(圖3e和f)，溝緣較易區(qū)分。出現(xiàn)此類狀況的19條切溝的長度、最大表面寬度和面積相對誤差僅3%；周長誤差略高，為5.5%。其余14條溝壁無裸露的切溝各參數(shù)誤差均較高。

3.2 切溝發(fā)育速率

3.2.1 集水區(qū)基本情況 將六道溝小流域劃分為50個集水區(qū)，包括37個Ⅰ級集水區(qū)、12個Ⅱ級集水區(qū)和1個Ⅲ級集水區(qū)(圖4e)，面積平均23.4 hm2。其中，Ⅰ-17集水區(qū)面積最小，僅為2.7 hm2；Ⅱ-3面積最大，高達(dá)120.5 hm2。集水區(qū)平均坡度介于12.1°～25.1°之間，坡度較大的集水區(qū)主要為分布在小流域北部的Ⅰ級集水區(qū)，如Ⅰ-21和Ⅰ-27，坡度均超過22°。平均坡長范圍12.8～33.0 m，小流域西側(cè)的集水區(qū)坡長總體高于東側(cè)。平均剖面曲率和平面曲率分別介于-0.4～0.3和-0.1～0.1之間。2013—2018年均植被蓋度介于27.9%～78.2%之間，Ⅰ-24土地利用以獨立工礦用地為主，植被蓋度最低；Ⅰ-34草地、灌木林地和有林地面積廣，植被蓋度最高。綜合土地利用動態(tài)度在1.6%～12.8%/a之間變化。土地利用轉(zhuǎn)移趨勢主要表現(xiàn)為草地轉(zhuǎn)出為灌木林地。

圖3 WV2影像與RTK GPS提取的典型切溝輪廓Fig.3 Typical gully boundaries extracted by WV2 image and RTK GPS

圖4 六道溝小流域集水區(qū)劃分及2個典型集水區(qū)2012與2018年解譯切溝分布Fig.4 Division of drainage areas in Liudaogou catchment and interpreted gullies in two typical drainage areas in 2012 and 2018

3.2.2 切溝發(fā)育狀況 根據(jù)2012和2018年的WV2解譯結(jié)果，統(tǒng)計各集水區(qū)的切溝總長度、總最大表面寬度、總周長和總面積，并基于3.1.1建立的轉(zhuǎn)換方程對WV2提取的切溝總周長和總面積進行轉(zhuǎn)換，計算得到各參數(shù)對應(yīng)的年均變化率分別介于-2.1%～1.2%、-1.7%～0.5%、-1.8%～1.4%和-0.2%～1.5%之間，平均-0.02%、0.01%、0.01%和0.2%。這與Li等[8]利用QuickBird影像解譯得到的黃土高原蔡家川流域2003—2010年切溝周長與面積年均變化率十分相近，后者范圍分別為-1.0%～1.2% 、0.1%～1.4%。

溝頭溯源侵蝕與溝岸側(cè)向侵蝕是切溝侵蝕在水平方向上的2種主要形式，其結(jié)果分別體現(xiàn)在切溝總長度和總最大表面寬度。統(tǒng)計各集水區(qū)切溝的溝頭后退與溝岸擴張距離發(fā)現(xiàn)，2012—2018年各集水區(qū)的溝頭后退總距離介于-47.4～29.1 m之間，溝岸擴張總距離介于-15.7～10.4 m。其中，9個集水區(qū)(18%)的切溝總長度和10個集水區(qū)(20%)的總最大表面寬度不增反減，尤以Ⅰ-22和Ⅰ-2集水區(qū)最為劇烈， 溝頭分別后退-47.4和-36.5 m，溝岸分別擴張-5.3和-15.7 m，這分別與集水區(qū)內(nèi)獨立工礦用地的挖方擴張、徑流小區(qū)的填土擴建有關(guān)。另有26個集水區(qū)(52%)的切溝總周長減少，這可能與相鄰切溝的兼并有關(guān)[8]。切溝總長度、總最大表面寬度、總周長和總面積年均增長率主要集中在0～0.5%內(nèi)，分別占集水區(qū)總數(shù)的78%、80%、40%和78%，其中Ⅰ-23切溝擴張最為迅速，溝頭后退9.4 m，溝岸擴張2.5 m，總長度、總周長和總面積年均變化率皆達(dá)到1.2%以上(圖4a和b)。

3.2.3 切溝發(fā)育的影響因素 小流域各集水區(qū)的切溝發(fā)育速率存在明顯差異。由于小流域面積較小，氣候和土壤的空間差異不大，筆者僅探討地形、植被蓋度和土地利用等屬性對切溝發(fā)育的影響，建立各屬性與切溝各參數(shù)變化率的相關(guān)系數(shù)矩陣(表1)。

表1 集水區(qū)切溝總參數(shù)變化率與集水區(qū)屬性相關(guān)系數(shù)

切溝發(fā)育速率與集水區(qū)面積通常呈顯著的正相關(guān)關(guān)系[8]。然而，筆者發(fā)現(xiàn)切溝各參數(shù)變化率與集水區(qū)面積的均無顯著相關(guān)(表1)。原因可能是集水區(qū)面積僅粗略反映了流向溝頭的徑流規(guī)模，實際徑流量還受集水區(qū)形狀、地形、土地利用、土壤性質(zhì)與降雨空間分布的綜合影響[9]。

坡度和坡長對坡面匯水過程及徑流量大小等有重要影響[18]。切溝總長度變化率與集水區(qū)平均坡度呈顯著負(fù)相關(guān)關(guān)系(P<0.01)，切溝總面積變化率與集水區(qū)平均坡長顯著正相關(guān)(P<0.05)，即切溝擴張速率隨坡度減小、坡長增加總體呈增大趨勢。當(dāng)其他條件一致時，坡度可影響坡面的承雨面積，坡度越大，承雨面積越小，降雨量與產(chǎn)流量相應(yīng)減小[19]。陡峭坡地不易受牲畜與人類擾動可能也是切溝擴張速率與坡度呈負(fù)相關(guān)的一個原因[20]。隨著坡長增加，匯流面積增大，徑流量增加[21]，往往會加劇切溝擴張。

在坡向方面，迎風(fēng)坡(E、S、SE、SW)集水區(qū)的切溝擴張速率總體高于背風(fēng)坡(W)(圖5c)。在六道溝小流域，切溝主要發(fā)生在夏季強降雨期間，此時當(dāng)?shù)匾詵|南風(fēng)為主，冬春季則多為西北風(fēng)。

曲率在很大程度上控制著水流和泥沙的運移方向[22]，尤其是剖面曲率，直接影響徑流的加速與減速，進而作用于土壤的侵蝕和沉積。切溝總長度和總面積變化率與集水區(qū)平均剖面曲率呈顯著正相關(guān)(表1)，即剖面曲率越大，切溝擴張越快。平面曲率影響徑流的匯聚和分散，對切溝發(fā)育的影響更為復(fù)雜。對各集水區(qū)的平面曲率進行分級統(tǒng)計發(fā)現(xiàn)，切溝總面積變化率隨其總體無明顯變化(圖5d)。

圖5 切溝總面積變化率隨主要集水區(qū)屬性的變化Fig.5 Change rate of total gully area with main drainage area properties

植被可增加地表粗糙度，促進入滲，減少地表徑流，從而減緩切溝侵蝕[23]。盡管集水區(qū)年均植被蓋度與切溝各形態(tài)參數(shù)變化率之間均無顯著相關(guān)關(guān)系(表1)。但按照SL190—2007《土壤侵蝕分類分級標(biāo)準(zhǔn)》將植被蓋度劃分為4級：<45%、45%～60%、60%～75%、>75%，統(tǒng)計不同植被蓋度的切溝總面積變化率發(fā)現(xiàn)，后者隨植被蓋度的增加而減小(圖5e)。當(dāng)植被蓋度<45%時，切溝總面積變化率平均高達(dá)0.7%/a；當(dāng)植被蓋度>45%時，切溝面積變化率均降至0.2%/a左右，表明蓋度大于45%時，植被可顯著抑制切溝發(fā)育，這與陳一先等[24]利用1976與1997年1∶1萬地形圖估算陜北黃土區(qū)切溝發(fā)育速率得到的結(jié)論一致。

除切溝總周長變化率外，綜合土地利用動態(tài)度與切溝各形態(tài)參數(shù)變化率均無顯著相關(guān)關(guān)系(表1，圖5f)。2012—2018年間土地利用變化多表現(xiàn)為草地轉(zhuǎn)出為灌木林地，但這2種土地利用類型的產(chǎn)流系數(shù)并無顯著差異[25]，在降水特征無明顯趨勢性變化的背景下，徑流特征和切溝侵蝕也不會發(fā)生太大的變化。

4 結(jié)論

1)WV2目視解譯提取切溝形態(tài)參數(shù)具有較高的精度。與RTK GPS實測數(shù)據(jù)相比，WV2解譯切溝長度、最大表面寬度、周長和面積的相對誤差主要集中在5%以內(nèi)，分別占切溝樣本總數(shù)的94%、79%、58%和76%。 當(dāng)切溝溝頭活躍、溝緣以草本植物為主或溝壁存在土體裸露時，WV2對切溝的解譯精度較高，誤差較小。

2)2012—2018年間，小流域各集水區(qū)切溝總長度、總最大表面寬度、總周長和總面積的年均變化率分別介于-2.1%～1.2%、-1.7%～0.5%、-1.8%～1.4%和-0.2%～1.5%之間，平均-0.02%、0.01%、0.01%和0.2%。

3)小流域切溝擴張速率與集水區(qū)平均坡長和平均剖面曲率呈顯著正相關(guān)，與平均坡度呈負(fù)相關(guān)，與坡向、植被蓋度也存在一定的相關(guān)關(guān)系。當(dāng)位于背風(fēng)坡或植被蓋度>45%時，切溝擴張速率較低。