基于SNIC-CNN-SVM模型的京津風沙源二期工程區土地利用/土地覆蓋遙感識別研究

李長龍,李增元,高志海,孫 斌,*,王絲絲

1 廣州商學院信息技術與工程學院, 廣州 511363 2 中國林業科學研究院資源信息研究所, 北京 100091 3 國家林業和草原局林業遙感與信息技術重點實驗室, 北京 100091 4 國家遙感中心, 北京 100084

土地覆被信息不僅是單一地表類型的表達,而且體現一定區域內生態環境和人類活動等自然和非自然屬性的綜合特征[1]。因此,不同的土地覆蓋類型代表著地球表面豐富多樣的生態系統和自然人文景觀,它是生態系統內部長期動態演變的積累,也是量變到質變的外部反映[2—3]。通過對土地覆被長時間序列變化研究能夠分析不同生態系統內部生物和非生物的綜合生態效應,從而掌握生態系統的長期演變規律,為全球氣候變化及其人類活動影響等提供重要信息支撐[4—6]。

隨著遙感技術的飛速發展,特別是深度學習算法的不斷優化,依靠遙感技術自動半自動地提取土地覆被信息已然成為當前的重要手段[7—9]。近年來,在國家和區域尺度,部分學者研發了多套地表覆蓋遙感產品,例如劉良云團隊GLC_FCS30[10]、宮鵬團隊FROM_GLC30[11]和FROM_GLC10[12]、國家自然資源部GlobeLand30[13]、歐空局全球陸地覆蓋數據(ESA GlobCover)[14]等,這些都極大地提高了氣候變化和生態環境評估等綜合性、大尺度深度研究的效率,為地理國情監測等提供了重要的基礎信息。

稀疏植被區主要是指植被覆蓋度在10%—50%的地表區域,包括草原、苔原、荒漠草原、荒漠、巖溶地區和部分高海拔植被區等。其中,零星的植被由草本和(或)木質和半木質植物組成,其余區域為自然裸地[15—16]。京津風沙源治理二期工程區從東南向西北降水逐步減少,稀疏植被分布廣泛。由于稀疏植被區內的自然裸地多為反射率較高的沙礫或巖石,而植被又零散分布在自然裸地內部,這就造成稀疏植被的光譜信息易被淹沒,提取難度較大[17—18],而針對不同地類的有效性特征是提高地表覆蓋制圖精度的重要因素[19]。目前,全球或全國的地表覆蓋產品尺度較大,針對稀疏植被區的特異性較差,類型識別精度較低,而且缺乏對不同沙化土地類型的精細劃分[20—21]。這些都極大地限制了稀疏植被區生態效應評價、系統管理和科學決策等方面的應用需求,因此迫切需要建立基于大尺度的稀疏植被區土地覆被信息的遙感提取方法。

根據遙感分類的影像基礎可以將土地覆蓋分類方法細化為基于像元、基于超像素和面向對象三類[22]。基于像元的方法較為傳統,主要依靠影像的光譜信息對不同地物進行監督或非監督的聚類,但是缺乏對幾何和周圍影像信息考慮[23—25]。面向對象的方法主要有分水嶺變換算法[26]、均值漂移算法[27]、分形網絡進化算法[28]和馬爾可夫隨機場分割算法[29]等,分類的基礎是對遙感圖像信息提取后分割成的對象,圖像分割優劣直接影響分類精度,不同地物類型對應不同分割尺度,計算量較大,效率較低,最優分割尺度的選擇多以人工或人機交互方式進行,不利于建立自動化提取模型[30—31]。超像素分割則彌補兩者缺點,通過對像素進行相似性分組,從而形成相對規則和均勻大小對象,以類似過分割的方式對圖像信息進行解析[32],去掉了最優尺度選擇和對象形狀特征提取等工作量,分析作用更強,更有利于沙化土地類型中植被和沙地邊界模糊或漸變性邊界的提取和劃分,目前常用的算法包括基于圖論的算法和基于梯度下降的算法等[33]。以卷積神經網絡(Convolutional Neural Networks,CNN)為代表的深度學習方法解決了最鄰近、極大似然等傳統遙感分類方法在遙感圖像中隱藏的深層次語義信息缺乏的問題,能夠從本身的深層結構中學習高度抽象的圖像語義區別特征[34—35],對于稀疏植被區土地覆蓋分類具有較強的適用性,能夠一定程度解決該區域類型復雜,變化多樣,植被分布很不均一,漸變型或混雜型地類較多等難題。

綜上,本研究以京津風沙源治理二期工程區為研究區,基于多源遙感數據,通過建立懲罰性機制優化超像素分割算法,研究構建用于稀疏植被區典型土地利用/土地覆蓋分類識別的SCS模型,以期提高稀疏植被的遙感信息提取能力,獲取科學準確的土地覆被類型信息,為研究區生態工程效應評價奠定基礎。

1 研究區概況

本研究以京津風沙源治理二期工程區為研究區,該區域東西距離約1300km,南北距離約1000km,總面積約為70.98萬km2,包括陜西省、內蒙古自治區、北京市、天津市、山西省和河北省的共計138個縣(旗、市、區),地理坐標為105°8′—120°59′E,36°48′—46°47′N,具體位置如圖1所示。

研究區地貌由平原、山地、高原三大部分組成,地形起伏較為復雜,最低處僅有幾十米,最高處超過了3000m。區域內的沙化土地聚集分布,形成了“三沙地一沙漠”的格局(科爾沁沙地、渾善達克沙地、毛烏素沙地和庫布齊沙漠)。研究區受到降水、氣溫和太陽直射時間等自然條件的影響,從東南向西北的氣候分異明顯,依次包括有亞濕潤區、干旱亞濕潤區、半干旱區和干旱區,植被則依次分布有森林草原、典型草原、荒漠草原、草原化荒漠和荒漠。由于降水和地形等自然因素的影響,研究區中部和西北部呈現明顯植被和荒漠交錯的破碎化自然景觀,屬于典型的稀疏植被區,其中植被沿低洼處或河流周圍聚集性分布,種類多樣且林灌草多混雜生長；裸露地表的差異性也較大,依次分布有栗鈣土、砂質土、粗沙、礫石等。另外,再疊加人類活動和林農牧交錯區影響,多種地表類型在一定區域內共存或多次更迭,光譜信息均一性較差,遙感信息提取難度較大[36—37]。

2 研究方法

2.1 數據獲取與處理

本研究主要獲取的遙感數據為Landsat影像,其中,2000年和2010年采用的是Landsat 5影像,總景數分別為274景和255景；2020年采用的是Landsat 8影像,總景數為322景,數量分布如圖2所示。每一年影像均選取研究區6—9月生長季影像,然后基于像元的最小云量進行影像合成,從而最大程度減少云對影像分類的干擾。影像的波段選擇為藍、綠、紅、近紅外、短波紅外和熱紅外,均采用經過大氣校正和幾何校正等預處理后的數據。

為了能夠更好地識別沙化土地,需要去除表層植被的干擾和影響,因此本研究采用SAR數據穿透表層植被,強化沙化土地的識別能力。根據實驗分析,極化方式為VV+VH組合的雙極化影像對沙化土地的敏感性更強,因此數據選擇Sentinel- 1A中的GRD一級產品,極化方式為VV+VH組合的雙極化,成像時間為2020年7月15日。數字高程模型(DEM)數據來自NASA的地球觀測系統數據信息系統,空間分辨率30m,成像時間為2010年,包含高程、坡度、坡向三類數據。為了與Landsat數據進行匹配,SAR數據經正射校正和重采樣,DEM數據經過尺度轉換,均統一到空間分辨率為30m。真值數據的來源主要分為兩部分：野外調查和高分辨率遙感影像判讀,各數據的數量和年份分布如表1所示,空間分布如圖1所示。每年份的野外調查數據采樣時間均為植被的生長季(7月或8月),影像判讀數據則是基于Google Earth上更高分辨率遙感影像進行。

表1 真值數據來源和年份分布

2.2 研究方法

本研究在原有超像素分割算法——簡單非迭代聚類(Simple Non-Iterative Clustering,SNIC)的基礎上,通過建立懲罰性機制,將更多的植被光譜信息引入影像分割中,增強稀疏植被遙感信息提取能力,然后利用SAR數據本身的穿透能力建立不同紋理特征對沙化土地進行識別,最后基于CNN機器學習方法和SVM分類器對整個研究區進行不同類別的確定,這樣就形成了SCS稀疏植被區遙感分類模型,系統地將(稀疏)植被、沙化土地和其他類別進行精準劃分。具體主要分為五步：分類體系的建立,特征選擇,超像素分割,影像分類和精度評價。

第一,分類體系的建立。參考自然資源部《土地利用現狀分類》(GB/T 21010—2017)和《第三次全國國土調查技術規程》(TD/T 1055—2019)分類標準,綜合Landsat遙感影像土地覆蓋信息提取能力,研究區分類體系定為8大類,分別是：喬木林、灌木林、草地、沙地、耕地、水體、戈壁和其他,含義如表2所示：

表2 土地覆蓋分類體系及含義

第二,特征選擇。主要分為光譜特征、指數特征和紋理特征,其中光譜特征為影像數據的波段信息；指數特征除了歸一化水體指數(NDWI)和歸一化植被指數(NDVI)外,增加了土壤調節植被指數(SAVI),其目的是去除稀疏植被區(特別是戈壁)背景土壤亮度較高的影響,三者的計算公式如下：

(1)

(2)

(3)

式中,NIR、Red和Green分別為近紅外、紅光和綠光波段的光譜值；L為植被密度參數。紋理特征是基于SAR數據(Sentinel- 1A)形成的七種,分別為差異性、對比度、方差、角二階矩、逆差矩、熵和相關性,示意圖如圖3所示。

圖3 SAR影像的7種紋理特征示意圖Fig.3 The seven texture features of SAR imageSAR：合成孔徑雷達 Synthetic Aperture Radar

第三,超像素分割。SNIC算法是將影像轉化為CIELAB顏色空間的五維特征向量,表征空間距離,即位置信息X、Y的兩維數據；表征顏色距離,即色彩信息的三維數據,最后基于這五維特征向量建立距離度量標準,進而對圖像進行局部聚類的過程。算法的實現過程主要有四步：生成種子點；種子點優化；最優距離的求解；優先隊列算法像元聚類。其中最優距離的求解如公式4—6所示,優先隊列算法則是基于K個種子點建立K個元素ei={(xi,yi),(li,ai,bi),k,Di,k},這K個元素建立最原始的優先隊列Q,然后Q每次返回隊列中到第k個種子點的距離Di,k最小的元素ei,這樣可以將多次迭代優化成一次迭代,極大地提高運行效率。

(4)

(5)

(6)

式中,Ds、Dc和D分別代表顏色距離、空間距離和最優距離,s和m為顏色距離和空間距離的歸一化參數。

原有SNIC算法將空間距離和顏色距離作為同等重要參數進行計算,更易形成較為規則的分割圖像,這樣在稀疏植被區土壤和植被的交界處,高亮的背景土壤會淹沒稀疏植被本身的光譜信息,造成影像分割誤差。因此,本研究在原有算法后增加第五步,懲罰性機制和連通性檢查。即對于某個屬于種子點a的像元jk,假設它的空間距離參數屬于種子點a,而顏色距離參數屬于種子點b,且滿足下面公式：

Da≥ε1,Db≤ε2

(7)

式中,Da和Db分別為jk屬于種子點a和種子點b的類內距離,ε1和ε2為最大類內距離。經連通性檢查,如果像元jk與種子點b的類內像元相連接,那么該像元jk將被標簽為種子點b的類內元素。

第四,影像分類。運用CNN算法挖掘深層特征,支持向量機(Support Vector Machine,SVM)進行精確分類,分類過程如圖4所示。CNN算法通過兩次卷積,兩次采樣和一次全連接算法,獲得卷積層權重參數,將一維特征影像轉換成256維。假設在CNN的卷積層中,輸入的一維特征向量表示為xk-1,則其經過卷積計算加上偏置參數的非線性變換,輸出特征向量為xk,計算過程可以表示為：

xk=f(wk×xk-1+bk)

(8)

式中,wk和bk表示卷積層k的卷積核參數(5×5)和偏置參數；*表示卷積計算符號,f()表示激活函數,本研究采用ReLU函數f(x)=max(0,x)進行非線性計算。

圖4 CNN-SVM算法分類過程示意圖Fig.4 The classification process of CNN+SVM algorithmCNN：卷積神經網絡 convolutional neural networks；SVM：支持向量機 support vector machine

第五,精度評價。采用混淆矩陣(Confusion matrix),包括生產者精度(Producer accuracy)、用戶者精度(User accuracy)、總體精度(Overall accuracy)和Kappa系數(Kappa coefficient)。

3 結果與分析

3.1 超像素分割結果

為了呈現超像素分割過程優化效果,選取渾善達克沙地腹地部分區域,對比了原始影像、SNIC算法優化前后的分割結果(圖5)。結果顯示,優化的SNIC算法在地類內部分割效果與未優化前相比,變化不大。但在不同地類間,特別是在稀疏植被和沙化土地類別的邊界(圖5的紅圈部分)區域,原有算法區分度不夠,在多個位置的地類界限具有明顯的混雜現象,而優化算法有效解決了這一問題,使分割效果更加清晰,能夠準確細致地對稀疏植被和沙化土地的邊界進行刻畫,有利于稀疏植被和沙化土地信息的區分。

整體而言,對于整個工程區中諸如沙地和草地、沙地和戈壁、草地和灌木林地、草地和森林等易混分的土地覆蓋類型,優化算法的遙感信息提取能力更強,分割效果更好,特別是各易混類別的邊界區分度更高,這樣更有利于下一步影像分類,進而提升最終的信息提取精度。

圖5 SNIC算法優化前后的影像分割對比Fig.5 Comparison of image segmentation before and after SNIC algorithm optimizationSNIC：簡單非迭代聚類 Simple non-iterative clustering

3.2 分類精度評價

優化前后SCS分類結果的混淆矩陣分別如表3和表4所示。結果表明,基于原有SCS模型的工程區典型土地利用/土地覆蓋的分類識別總體精度為78.24%,Kappa系數為0.7452,而經過優化,模型的總體分類精度顯著提高,達到89.41%,Kappa系數為0.8760,總體提高了11.17%,Kappa系數提高了0.1308。其中對于工程區面積最大的兩種土地利用/土地覆蓋類型——草地和沙地,兩者的用戶者精度分別提高了14%和9%。綜合整體分類精度和典型地類單一識別精度結果發現,本研究提出的SCS模型優化方案能夠有效提高稀疏植被區域的植被與沙地信息提取精度,能夠滿足大區域尺度對稀疏植被區生態效應、環境變化和人類活動等深入研究的需要。

表3 2020年原有模型土地覆蓋分類結果混淆矩陣

表4 2020年優化模型土地覆蓋分類結果混淆矩陣

3.3 2020年工程區典型土地利用/土地覆蓋空間分布

采用優化的SCS模型,獲得了2020年工程區典型土地利用/土地覆蓋的空間分布情況(圖6)。結果表明,總體上,研究區多以草地、林地和沙地為主,其中草地面積最大,占研究區總面積的51.52%,廣布于研究區的大部,東北部的烏珠穆沁草原區域是最重要的分布范圍；林地范圍次之,喬木林和灌木林分別占研究區總面積的19.55%和7.23%,分布范圍上來說主要位于河北省北部和內蒙古錫林郭勒盟南部以及零星分布于沙地和草地周邊；沙地則主要集中于一個沙漠三個沙地區域(庫布齊沙漠、毛烏素沙地、渾善達克沙地和科爾沁沙地),總面積占研究區總面積的11.96%；耕地主要分布于黃河沿岸和北京市、天津市以及河北山西的平原上,約占研究區總面積的1.54%；戈壁則在研究區的西北部烏拉特后旗,約占研究區總面積的5.20%。稀疏植被覆蓋(草地、沙地、戈壁)區域面積占比約為68.68%,研究區大部分區域屬于脆弱生態區。

圖6 研究區2020年土地利用/土地覆蓋類型分布Fig.6 The distribution of land use/land cover types in 2020

3.4 近20年工程區典型土地利用/土地覆蓋類型空間變化格局

自2000年以來,研究區相繼實施了退耕還林還草工程、野生動植物保護及自然保護區建設工程、天然林保護工程、京津風沙源治理工程等生態保護重大工程[38],在森林面積增加,沙化土地減少,區域經濟發展等方面取得了豐碩的成果,充分發揮了生態修復和改善,維護國家生態安全等方面的重要作用。由于京津風沙源治理工程(一期)實施時間為(2001—2010年),二期工程實施時間為(2013—2022年),因此本研究選擇2000年和2010年作為對比時間節點,深入分析近20年京津風沙源治理二期工程區土地利用/土地覆蓋類型變化情況。圖7所示為2000年和2010年土地利用/土地覆蓋類型分布,圖8所示為近20年工程區典型土地利用/土地覆蓋類型的變化情況。結果表明,整體上,近20年來喬木林、灌木林和草地的面積增加最多,其中,喬木林占研究區的比例由2000年的10.17%增加到2010年的19.22%,再到2020年的19.55%,增加了9.38%；灌木林的面積則從1.43%(2000年)增加到6.65%(2010年)再到2020年的7.23%,增加了5.79%；草地的面積由46.05%(2000年)下降到44.26%(2010年)再增加到51.52%(2020年),增加了5.47%。與之相應的是,耕地和沙地的面積明顯減少,兩者占工程區總面積的比例分別由2000年的13.33%和16.54%,減少到2020年2.43%和11.96%,分別減少了10.89%和4.58%。

圖7 研究區2000年和2010年土地利用/土地覆蓋類型分布Fig.7 The distribution of land use/land cover types in 2000 and 2010

圖8 研究區2000—2020年土地利用/土地覆蓋類型變化統計Fig.8 The statistics of land use/land cover changes in the study area from 2000 to 2020

在自然條件下,土地覆蓋類型的變化是量變到質變的過程。植被逐步退化最終變為荒地或沙地,而沙化土地隨著植被的逐步改善,最終變為草原或灌木林[39]。因此,本研究進一步從喬灌草植被覆蓋區和沙化土地兩方面的相互轉化情況來反映近20年來工程區的實施效果。變化情況如圖9所示。從2000年到2020年,位于渾善達克沙地、毛烏素沙地和科爾沁沙地的腹地的沙化土地相對穩定,而沙地邊緣植被恢復明顯,特別是杭錦旗、達爾罕茂明安聯合旗等沙地邊緣地帶,沙地轉換區域中的86.45%轉換為喬灌草植被覆蓋區；由其他土地利用/土地覆蓋類型轉為沙地的區域主要集中在錫林郭勒盟的蘇尼特左旗和蘇尼特右旗等地,轉換為沙地的區域中76.14%來自于原有喬灌草植被類型覆蓋區。植被方面,喬灌草土地覆蓋類型不變區域主要位于研究區北部和東部大部分區域；喬灌草土地覆蓋類型減少區域主要位于渾善達克沙地和科爾沁沙地邊緣地帶,減少區域中的33.41%轉換成沙化土地；喬灌草土地覆蓋類型增加區域主要位于毛烏素沙地、庫布齊沙漠區域,增加區域的25.62%來源于沙化土地。總體來說,京津風沙源治理工程實施的十多年時間內,沙化土地面積大幅度降低,喬灌草總面積大幅度上升,土地退化的趨勢得到遏制,草地質量和生態環境改善明顯,表明研究區植被狀況不斷改善,實施的各項生態工程作用顯著,能夠更有效地服務于多維度生態功能的發揮和評價。

圖9 2000—2020年研究區植被類型(喬灌草) 和沙化土地變化分布Fig.9 The change distribution of vegetation types (arbor, shrub, grass) and desertification land in the study area from 2000 to 2020

4 討論

前人的研究結果[40—41]表明,京津風沙源治理工程區土地利用/土地覆蓋變化在2000—2010年主要表現為沙地、耕地和林地的增加,草地的減少,2010—2018年則表現為耕地、沙地的減少,草地和林地的增加,恢復區域中以鄂爾多斯高原、渾善達克沙地和科爾沁沙地為主,這些結論都基本與本研究獲得的結果一致。只是部分認為2018年研究區內草地的分布占比為55.25%,略高于本研究51.52%(2020年)而沙地的占比略低。可能的原因是本研究利用SAR數據(Sentinel- 1A)的穿透能力,將研究區內表層具有稀疏植被,但土壤為沙礫的土地覆蓋類型劃分為沙化土地,一定程度上減少了草地占比而增加了沙地占比,但本研究的提取結果能夠更準確地表征沙化土地,具有更強的科學性和實際意義。

全球尺度方面,根據與劉良云、宮鵬和陳軍等人研究的結果[10—13]對比來看,在各類別的區域分布上與本研究的結果基本一致,但因三人的研究成果均基于全球尺度,在區域尺度上并未細致研究各類別的特異性特征,因此也缺乏對京津風沙源區域復雜的植被土壤環境的深入研究,易出現植被和沙化土地錯分的情況,特別是在稀疏植被覆蓋的沙化土地區域,而本研究在該方面有較為明顯的改善。

由于大尺度、長時間序列工程措施生態效應評價的需要,本研究采用的超像素分割算法一方面可以增加像元分類中缺少的對象信息,另外一方面又減少面向對象分割中最優尺度等人工參與,實現稀疏植被區沙化土地類型識別和土地覆蓋分類的自動化、全流程的快速提取,以SAR數據表征沙化土地為基礎的機器學習方法能夠更加細致刻畫稀疏植被區中沙化土地和稀疏植被的邊界。但由于驗證樣本不足,本研究在部分區域的精度還需進一步進行論證,特別是本研究主要聚焦于稀疏植被和沙化土地等自然地表類型,對于建筑物、道路等人工地表類型的識別能力還需深入研究。此外,可以拓展研究該區域外非工程區土地覆蓋類型的變化,從而與本研究的結果進行對比分析,進一步論證生態工程措施的積極作用。下一步,通過融合更多時序信息、專題信息等對研究區的土地利用/土地覆蓋類型進行更深入的分析,進一步提高該區域的遙感識別精度。

5 結論

本研究以典型大尺度稀疏植被區——京津風沙源治理二期工程區為研究區,通過構建懲罰性機制優化超像素分割算法,研建了SCS模型方案,實現了研究區典型林草沙要素信息提取和主要土地利用/土地覆蓋識別。研究結果表明：

1)引入懲罰性機制優化后的SNIC分割算法,有效提升了稀疏植被區與沙地區的邊界區分度,有助于分類精度的提升；

2)基于改進SNIC-CNN-SVM模型方案的研究區總體分類精度達89.41%,較優化前提高了11.17%,特別是針對草地和沙地這兩種工程區核心類型的識別精度分別提高了14%和9%,表明該模型具有較好的沙化土地和稀疏植被識別能力,該優化方案在以研究區為代表的稀疏植被區域分類中具有較好的應用效果和推廣價值；

3)2020年工程區中以草地面積最大,占研究區總面積的51.52%；沙地主要集中于一個沙漠三個沙地區域(庫布齊沙漠、毛烏素沙地、渾善達克沙地和科爾沁沙地),總面積占工程區總面積的11.96%；稀疏植被覆蓋(草地、沙地、戈壁)區域面積占比約為68.68%,表明工程區處在林地-稀疏植被-沙地的過渡地帶,生態環境保護壓力與防沙治沙形勢依然嚴峻。

4)通過對比近20年來的典型土地利用/土地覆蓋類型空間分布格局發現,喬灌草等植被增加面積約占工程區20.64%,主要由沙化土地轉化,沙化土地減少面積約占研究區的4.58%,沙地減少的區域中86.45%轉換為喬灌草等植被類型。植被覆蓋范圍增加,質量提升,沙化土地面積降低,特別是在毛烏素沙地和庫布齊沙漠等區域改善效果突出,表明經過20余年的治理,研究區土地退化趨勢得到遏制,草地質量和生態環境改善明顯,實施的各項生態工程作用顯著,能夠更有效地服務于多維度生態服務功能的發揮。該研究以期能夠為京津風沙源二期工程區的生態系統演變規律研究及生態工程評價等工作提供重要科學支撐。