孫灝，阮琳，郭航，徐丹，孫中平

(1.中國礦業(yè)大學(北京) 地球科學與測繪工程學院，北京 100083；2.山東新巨龍能源有限責任公司，山東 菏澤 274918；3.生態(tài)環(huán)境部衛(wèi)星環(huán)境應用中心，北京 100094)

0 引言

生態(tài)用地比例是指評價區(qū)域綠地、水域濕地和耕地面積占評價區(qū)域總面積的比例，是衡量生態(tài)系統(tǒng)宏觀構(gòu)成性能的重要指標，也是《生態(tài)環(huán)境狀況評價技術(shù)規(guī)范》中評價生態(tài)環(huán)境質(zhì)量的重要參數(shù)[1]。寧夏沿黃城市帶是以地緣相近、交通便利、經(jīng)濟關(guān)聯(lián)度較高的銀川市為中心，石嘴山、吳忠、中衛(wèi)三個地級市為主干，形成大中小城市相結(jié)合呈帶狀分布的城鎮(zhèn)集合體，是寧夏城鎮(zhèn)分布最密集、生產(chǎn)力要素最活躍的區(qū)域，也是寧夏經(jīng)濟社會發(fā)展的精華地帶。然而，工業(yè)化和城市化帶來巨大經(jīng)濟效益的同時，也可能引發(fā)諸多生態(tài)環(huán)境問題[2-5]。定量評價與分析寧夏沿黃城市帶的生態(tài)用地比例變化，有助于科學地揭示沿黃城市帶已經(jīng)和可能引發(fā)的生態(tài)環(huán)境問題，為寧夏可持續(xù)發(fā)展和生態(tài)文明建設提供重要的輔助決策資料。

目前，遙感技術(shù)已成為獲取生態(tài)用地比例的重要手段[6-7]。其中，基于現(xiàn)有土地覆蓋類型遙感產(chǎn)品計算生態(tài)用地比例，無需進行初始的分類操作，產(chǎn)品已經(jīng)過嚴格的精度驗證，能避免主觀誤差，且具有較為完善的分類體系，為大范圍生態(tài)用地比例的獲取提供了有利條件。然而，土地覆蓋類型產(chǎn)品多式多樣，且數(shù)據(jù)來源、制作原理等有著明顯差異[8-10]，致使基于典型土地覆蓋類型產(chǎn)品計算生態(tài)用地比例時標準難以統(tǒng)一。同時，分類系統(tǒng)的不同、語義分辨率的差異，也會導致各土地覆蓋類型產(chǎn)品的不同土地覆蓋類型間存在混淆，缺乏詳細的一致性分析等問題[11-15]。另外，鮮有研究利用全球土地覆蓋類型產(chǎn)品，定量評價與分析寧夏沿黃城市帶的生態(tài)用地比例變化。

本文選取四種典型的全球土地覆蓋類型產(chǎn)品，即MODIS-MCD12Q1[16]、GLOBCOVER[17-18]、GLOBELAND30[19-20]，以及FROM-GLC(finer resolution observation and monitoring of global land cover)[21-22]，從生態(tài)用地角度將土地覆蓋類型進行重分類，在統(tǒng)一分類體系下設計生態(tài)用地比例計算方法。并針對寧夏沿黃城市帶，從單個土地覆蓋類別以及總體生態(tài)用地角度，定量分析了四種典型產(chǎn)品的相對一致性。在此基礎上，選擇時間連續(xù)性較好的MODIS-MCD12Q1產(chǎn)品，利用曼-肯德爾法非參檢驗(Mann-Kendall，M-K)、森氏斜率估計(Sen’s slope)、赫斯特指數(shù)(Hurts)等時間序列分析技術(shù)，分析寧夏沿黃城市帶近20年生態(tài)用地比例的演變過程和趨勢，以期為寧夏沿黃城市帶的可持續(xù)發(fā)展規(guī)劃提供科學依據(jù)，也為其他生態(tài)環(huán)境問題研究提供參考。

1 數(shù)據(jù)與研究方法

1.1 研究區(qū)概況

寧夏沿黃城市帶地處黃河中上游引黃灌區(qū)(104°22′E～107°50′E，36°43′N～39°22′N)，屬干旱半干旱氣候，年均氣溫5～9 ℃，年降水180～680 mm，是生態(tài)環(huán)境保護研究的典型區(qū)域[23]。圖1為寧夏沿黃城市帶各區(qū)縣空間分布圖，圖中背景為2010年MODIS土地覆蓋類型產(chǎn)品MCD12Q1數(shù)據(jù)信息，采用IGBP分類系統(tǒng)[24]，空間分辨率為500 m。

注：該圖基于寧夏地理信息公共服務平臺下載的審圖號為寧S[2019]第018號的標準地圖制作，底圖無修改。圖1 研究區(qū)域空間分布圖

1.2 典型土地覆蓋類型遙感產(chǎn)品

四種典型土地覆蓋類型遙感產(chǎn)品的基本信息如表1所示，其中MODIS-MCD12Q1產(chǎn)品提供了2001—2017年，每年一幅空間分辨率為500 m的影像，具有良好的時序性。它是基于決策樹的監(jiān)督分類方法，并借助先驗知識和輔助信息做后續(xù)處理完成的分類。GLOBCOVER是歐洲空間局通過全球合作，利用ENVISAT/MERIS數(shù)據(jù)，采用非監(jiān)督分類方法生產(chǎn)而成。GLOBELAND30是首批分辨率達到30 m的土地覆蓋產(chǎn)品之一，它是利用Landsat TM、ETM+數(shù)據(jù)和中國環(huán)境減災衛(wèi)星(HJ-1)多光譜影像，結(jié)合大量輔助信息(全球MODIS NDVI年序數(shù)據(jù)、全球DEM數(shù)據(jù)等)，基于“像元-對象-知識”POK方法制作而成，并委托第三方進行精度評估。GLOBELAND30-2010數(shù)據(jù)的總體精度為83.51%，Kappa系數(shù)為0.78，具有很高可靠性和可信度。FROM-GLC是利用TM和ETM+數(shù)據(jù)，采用最大似然分類器、J4.8決策樹分類器、隨機森林分類器和支持向量機分類器四種分類器，最終合成全球30 m分辨率土地覆蓋圖。

表1 四種土地覆蓋類型產(chǎn)品概要信息

1.3 生態(tài)用地比例計算

1)土地覆蓋類型產(chǎn)品預處理。使用MRT工具將MODIS-MCD12Q1產(chǎn)品在WGS-1984坐標系下按UTM Zone 48N投影成空間分辨率480 m的影像。以MODIS-MCD12Q1為準，將GLOBCOVER、GLOBELAND30、FROM-GLC轉(zhuǎn)換成相同坐標系。應用最近鄰采樣，將GLOBCOVER聚合到480 m，利用眾數(shù)聚合，將GLOBELAND30和FROM-GLC聚合到480 m。

為能在同一分類基準下比較四類數(shù)據(jù)產(chǎn)品，從生態(tài)用地角度出發(fā)，對地物類型進行歸并[25]。歸并后共六類，分別為：自然植被、耕地、自然植被與耕地鑲嵌體、濕地、水域以及非生態(tài)用地。表2為歸并后類別與歸并前類別編碼的對應關(guān)系。

表2 統(tǒng)一分類系統(tǒng)后分類編碼信息

2)一致性分析。MODIS-MCD12Q1產(chǎn)品具有良好的時序性，故本文以其為參考數(shù)據(jù)，對GLOBCOVER、GLOBELAND30、FROM-GLC產(chǎn)品從面積一致性、空間位置一致性兩個方面進行對比[26]。面積一致性方面，主要對比分析各類別的面積占比?？臻g一致性方面，將土地覆蓋類別進一步聚合為生態(tài)用地(E)和非生態(tài)用地(NE)，經(jīng)過疊加分析生成四種數(shù)據(jù)類別，即E-E、E-NE、NE-E、NE-NE，最后運用式(1)計算四類產(chǎn)品的空間相似系數(shù)[27]。

(1)

式中：O代表空間相似系數(shù)；A代表E-E像元總數(shù)；B代表E-NE像元總數(shù)；C代表NE-E像元總數(shù)。

3)生態(tài)用地比例計算。依據(jù)生態(tài)用地比例定義，可將自然植被、耕地、自然植被/耕地鑲嵌體、濕地、水域歸入生態(tài)用地，得到式(2)。

(2)

式中：E為生態(tài)用地比例；X為生態(tài)用地的像元數(shù)目總和；Y為研究區(qū)范圍內(nèi)有效像元數(shù)目的總數(shù)(不包含無數(shù)據(jù)的像元)。為分析生態(tài)用地比例的時空演變趨勢，研究使用了Mann-Kendall非參數(shù)檢驗、Sen’s 斜率估計以及Hurts指數(shù)分析[28-29]。

2 結(jié)果

2.1 不同土地覆蓋類型產(chǎn)品的面積一致性

圖2對各產(chǎn)品土地覆蓋類型面積占比進行展示。與同期MODIS-MCD12Q1產(chǎn)品相比，各產(chǎn)品在自然植被、耕地、自然植被/耕地鑲嵌體、濕地、水域、非生態(tài)用地的差異值為：GLOBCOVER-2005為-57.75%、-11.81%、23.26%、-0.07%、0.18%、46.19%；GLOBELAND30-2010為-22.49%、13.88%、0.0%、0.06%、1.17%、7.38%；FROM-GLC為-38.55%、1.99%、0%、-0.08%、2.31%、34.33%(2015年)，-44.93%、4.31%、0%、0.05%、0.9%、39.67%(2017年)。結(jié)果表明，其主要差異均體現(xiàn)在自然植被和非生態(tài)用地類型層面，MODIS-MCD12Q1產(chǎn)品中大多數(shù)被識別為自然植被類型(草地)的像元，在其他產(chǎn)品中被分為非生態(tài)用地類型(裸地)，是這一差異的主要原因。

圖2 各產(chǎn)品與MCD12Q1產(chǎn)品土地覆蓋類型面積占比

自然植被、耕地、濕地以及水體間的相互錯分并不會引起生態(tài)用地比例計算的誤差，故將土地類型進一步歸并為生態(tài)用地和非生態(tài)用地2種類型，從而更為直觀地反映各產(chǎn)品在生態(tài)用地比例計算方面的差異。圖3是基于每種產(chǎn)品生對態(tài)用地比例的統(tǒng)計，與同期MODIS-MCD12Q1產(chǎn)品相比，生態(tài)用地比例差異值：GLOBCOVER-2005為-46%；GLOBELAND30-2010為-7%，具有最高一致性；FROM-GLC為-34%(2015年)、-40%(2017年)。寧夏沿黃城市帶區(qū)域，稀疏草地分布較廣，各產(chǎn)品對裸地和草地判別標準的不同是產(chǎn)生差異的主要原因。

圖3 不同年份產(chǎn)品間生態(tài)用地面積占比

2.2 不同土地覆蓋類型產(chǎn)品的空間一致性

圖4為三種典型土地覆蓋產(chǎn)品與同期MODIS-MCD12Q1空間一致性空間分布圖，其中E-E代表類別一致，E-NE和NE-E代表兩類產(chǎn)品在該處類別不一致。該圖表明，MODIS-MCD12Q1與GLOBELAND30在統(tǒng)計生態(tài)用地和非生態(tài)用地方面的空間一致性較高。表3統(tǒng)計了三種土地覆蓋產(chǎn)品(GLOBECOVER、GLOBELAND30以及FROM-GLC)與同期MCD12Q1產(chǎn)品生態(tài)用地總體空間相似系數(shù)。結(jié)果表明：FROM-GLC產(chǎn)品在2015、2017年，相似系數(shù)分別為60.93%、56.40%；GLOBELAND30-2010數(shù)據(jù)產(chǎn)品與MODIS-MCD12Q1產(chǎn)品相比，生態(tài)用地總體相似系數(shù)高達85.07%，具有較高的契合度。GLOBECOVER與MODIS-MCD12Q1產(chǎn)品的相似系數(shù)最低，僅為44.57%。

注：該圖基于寧夏地理信息公共服務平臺下載的審圖號為寧S[2019]第018號的標準地圖制作，底圖無修改。圖4 空間一致性分析(E為生態(tài)用地;NE為非生態(tài)用地)

表3 三種全球土地覆蓋產(chǎn)品與同期MCD12Q1產(chǎn)品生態(tài)用地總體空間相似性 %

2.3 寧夏沿黃城市帶生態(tài)用地比例演變趨勢分析

通過不同產(chǎn)品間的一致性分析可知，MODIS-MCD12Q1與高分辨率產(chǎn)品GLOBELAND30-2010有較高的一致性。其中，GLOBELAND30-2010數(shù)據(jù)產(chǎn)品經(jīng)第三方評價，具有很好的精度保證，間接驗證了MODIS-MCD12Q1產(chǎn)品的精度和可靠性。同時，MODIS-MCD12Q1產(chǎn)品具有良好的時間連續(xù)性，適宜于長時間序列的研究分析。故本文基于MODIS-MCD12Q1土地覆蓋產(chǎn)品，以寧夏沿黃城市帶為例，探究生態(tài)用地比例的演變過程。

從圖5可以看出，寧夏沿黃城市帶生態(tài)用地比例主要分為三種趨勢類型，即增長型、波動型、縮減型。大部分地區(qū)生態(tài)用地比例呈明顯上升趨勢，屬增長型，如惠農(nóng)縣、大武口區(qū)；部分區(qū)縣生態(tài)用地比例無明顯變化趨勢，屬波動型，如中寧縣、賀蘭縣；金鳳區(qū)、西夏區(qū)兩個區(qū)生態(tài)用地比例呈持續(xù)減少趨勢，屬縮減型。但寧夏沿黃城市帶整體上呈明顯增長趨勢。

圖5 研究區(qū)域生態(tài)用地比例時間序列圖

對寧夏沿黃城市帶2001—2017年的生態(tài)用地比例進行時間上的趨勢顯著性分析和定量變化率估計，如圖6所示。結(jié)果表明：1)非參檢驗角度，除金鳳區(qū)和西夏區(qū)生態(tài)用地比例呈下降趨勢外，其余地區(qū)生態(tài)用地比例均呈上升趨勢。金鳳區(qū)和西夏區(qū)M-K趨勢統(tǒng)計量絕對值均大于2.32，表明其下降趨勢通過99%的顯著性檢驗，即呈極顯著減少趨勢；生態(tài)用地比例趨勢上升的所有地區(qū)中只有大武口區(qū)和利通區(qū)的M-K趨勢統(tǒng)計量不滿足99%的顯著性檢驗，但其絕對值也大于1.64，通過了95%的顯著性檢驗，即呈顯著增加趨勢。2)Sen’s斜率估計中，僅有金鳳區(qū)和西夏區(qū)的估計量小于零，其他區(qū)縣及研究區(qū)整體生態(tài)用地比例的估計量均大于零，即研究區(qū)內(nèi)大部分地區(qū)生態(tài)用地比例都呈上升趨勢，金鳳區(qū)和西夏區(qū)生態(tài)用地比例呈下降趨勢；惠農(nóng)區(qū)估計值最大，表明惠農(nóng)區(qū)生態(tài)用地比例增大趨勢最明顯。

圖6 研究區(qū)域生態(tài)用地比例時間序列分析結(jié)果

圖7為M-K統(tǒng)計量和Sen’s斜率估計空間分布圖，更直觀地展現(xiàn)寧夏沿黃城市帶核心區(qū)域生態(tài)用地比例近20年的變化趨勢及速率。

注：該圖基于寧夏地理信息公共服務平臺下載的審圖號為寧S[2019]第018號的標準地圖制作，底圖無修改。圖7 研究區(qū)域生態(tài)用地比例時間序列分析結(jié)果空間分布圖

3 討論

Hurst指數(shù)可以定量刻畫時間序列發(fā)展趨勢的持續(xù)性或反持續(xù)性強弱程度[30-31]，Hurst指數(shù)值為0.5表示序列為獨立同分布的隨機序列，屬于布朗運動，即現(xiàn)在不會影響未來；大于0.5表示時間序列數(shù)據(jù)趨勢具有持續(xù)性；小于0.5表示時間序列數(shù)據(jù)趨勢具有反持續(xù)性。本文選擇Hurst指數(shù)定量描述研究區(qū)生態(tài)用地比例的持續(xù)性或反持續(xù)性強度，結(jié)果如圖8所示。

注：圖(b)基于寧夏地理信息公共服務平臺下載的審圖號為寧S[2019]第018號的標準地圖制作，底圖無修改。圖8 研究區(qū)域生態(tài)用地比例Hurst指數(shù)及其空間分布圖

圖8表明，除西夏區(qū)外其余各區(qū)和研究區(qū)整體的生態(tài)用地比例序列Hurst指數(shù)均大于0.80，表明這些區(qū)縣生態(tài)用地比例時間序列變化趨勢持續(xù)性很強，即金鳳區(qū)將保持其下降趨勢，其他區(qū)縣將保持上升趨勢。西夏區(qū)Hurst指數(shù)也達到了0.75以上，表明其持續(xù)性也較強，也將繼續(xù)保持下降趨勢。

圖9為基于MODIS-MCD12Q1產(chǎn)品的土地覆蓋類型變化圖?？梢钥闯觯簩幭难攸S城市帶自然植被類型在2006年之前呈微弱下降、之后呈現(xiàn)增長趨勢；耕地先呈緩慢增加，后趨于變化穩(wěn)定的趨勢；裸地類型占比在2006年之后先迅速下降，后趨于穩(wěn)定；建設用地則逐年穩(wěn)步增加。該圖表明寧夏沿黃城市帶發(fā)展過程中建設用地的增加主要來源于裸地覆蓋類型的減少，并未以破壞生態(tài)用地為代價，因此研究區(qū)整體生態(tài)用地比例呈現(xiàn)良好發(fā)展態(tài)勢。

圖9 研究區(qū)域主要土地覆蓋類型變化分析圖

4 結(jié)束語

本文比較分析了MODIS-MCD12Q1、GLOBCOVER、GLOBELAND30和FROM-GLC四種典型全球土地覆蓋產(chǎn)品的相對一致性?；贛ODIS-MCD12Q1運用M-K非參檢驗、Sen’s斜率估計以及Hurst指數(shù)分析了寧夏沿黃城市帶近20年生態(tài)用地比例演變過程和趨勢，結(jié)論如下。

1)從生態(tài)用地角度，MODIS-MCD12Q1產(chǎn)品與GLOBELAND30-2010數(shù)據(jù)產(chǎn)品面積差異僅為-7.38%，空間相似系數(shù)為85.07%，一致性程度最高；MODIS-MCD12Q1與FROM-GLC產(chǎn)品的一致性次之，與GLOBCOVER產(chǎn)品的一致性最低。同時發(fā)現(xiàn)，各產(chǎn)品在寧夏沿黃城市帶的差異主要源于草地和裸地間的混淆。

2)M-K非參檢驗和Sen’s斜率估計結(jié)果表明，2001—2017年，寧夏沿黃城市帶除金鳳區(qū)和西夏區(qū)生態(tài)用地比例呈下降趨勢外，其余區(qū)縣均呈增長趨勢或在小范圍內(nèi)波動，研究區(qū)總體呈現(xiàn)顯著增長趨勢。Hurst指數(shù)表明，各區(qū)趨勢具有強、很強的持續(xù)性。通過對自然植被、耕地、建設用地以及裸地類型面積占比的趨勢分析，發(fā)現(xiàn)總體生態(tài)用地比例增長趨勢主要來源于自然植被和耕地的增加以及裸地類型的減少，建設用地的增加并未導致整體生態(tài)用地比例的減少。

雖然建設用地逐漸擴張有可能侵占生態(tài)用地，但通過生態(tài)環(huán)境綜合治理(如對裸地等非生態(tài)用地綠化改造)，可維持整體生態(tài)用地比例的穩(wěn)定，為寧夏沿黃城市帶可持續(xù)發(fā)展提供保障。