王世東 馮正英 余 洋,2 張合兵
(1.河南理工大學測繪與國土信息工程學院, 焦作 454000; 2.中國地質環境監測院, 北京 100081)
近年來,國內外學者較多利用轉移矩陣[1-2]和土地利用動態變化模型[3-5]對土地利用/覆被動態變化進行研究,分析土地利用和覆被的數量變化、程度變化及空間分布。上述研究大多是基于單期或相鄰兩期數據的對比揭示在多時間序列上土地利用變化的時空規律。構建土地利用變化軌跡模型不僅可以得到其時空變化特征,還可以通過各種軌跡分析方法[6-10]分析其轉移和流向[11],從軌跡過程中尋求變化規律,揭示時間序列上土地利用/覆被變化的時空動態演變特征[12]。
土地利用/覆被變化研究是一個長期動態過程,以上研究較多只是刻畫了一個相對“連續”時間尺度[13-14]的變化,若兩期間隔時間較長,勢必會割裂其連續變化,破壞變化規律中的過程完整性。因此,在遙感影像的支持下,研究基于連續時間尺度的土地利用/覆被變化軌跡,對于精細的土地利用動態監測及土地資源的可持續利用具有重要意義。當時間序列數據較多時,土地利用/覆被變化軌跡更加繁雜,軌跡提取和時空演變分析愈加困難,因此需要一種針對多時間序列數據的軌跡提取和分析方法。目前,穩定映射變化軌跡法已廣泛應用于各領域的軌跡變化中,但該方法在土地利用變化方面的研究較少,尤其是針對多時間序列土地利用/覆被變化的研究。本文結合ArcGIS技術,基于長時序遙感數據和改進穩定映射法的土地利用/覆被變化(Land use/cover change,LUCC)軌跡分析方法,研究時相數與穩定映射變化軌跡STD指標的對應關系,為長時序連續時間尺度的土地利用變化軌跡提取和變化軌跡類型劃分提供借鑒。
丹江流域水資源豐富,位于河南省淅川縣的丹江口水庫是南水北調中線工程的主要水源地[15],自2003年南水北調工程實施以來,不僅改變了丹江口庫區水資源的自然地理分布,同時受人類活動影響,丹江口庫區的人口密度增大、土地開發利用強度增加、環境資源過度消耗等也使得丹江口庫區的土地利用/覆被變化更加顯著[16]。鑒于此,本文以丹江流域(河南段)為例,選用2002—2017年16期Landsat TM/OLI及HJ-1A CCD遙感影像數據,提取連續時間尺度的土地利用變化圖譜,建立改進穩定映射法的LUCC軌跡分析方法,提取變化軌跡,并對軌跡類型進行劃分,同時結合土地利用動態度等定量模型和景觀指數,分析研究區16年間土地利用和覆被的轉移流向總體特征和時空變化規律。
丹江流域(河南段)位于河南省西南部,地理坐標為東經110°35′~111°58′,北緯32°30′~34°1′,南接湖北省,西臨陜西省,總面積約為8 419 km2。丹江干流自淅川縣荊紫關入境,匯水總面積占整個研究區的3.11%,流經河南省洛陽市、南陽市和三門峽市,包含三門峽市的盧氏縣,洛陽市的欒川縣部分區域,南陽市的西峽縣、淅川縣、內鄉縣以及鄧州市部分區域,其地理位置如圖1所示。

圖1 研究區地理位置圖Fig.1 Geographical location map of study area
研究區具有典型的流域特色,整體地形起伏較大,地勢由西北向東南傾斜,西北部環繞伏牛山脈,東南部依次為崗地和沖擊平原區,屬亞熱帶向暖溫帶過渡的季風性氣候。研究區景觀分布呈規律性,林地主要分布于研究區西部和北部的山地丘陵地帶,大部分為塊狀山林,且受河流切割使得地貌起伏跌宕;耕地主要分布于東南部的平原地帶;水域主要位于西南部淅川縣境內的部分丹江口庫區,是南水北調中線工程的主要水源地。2003年南水北調中線工程正式開工,2005年丹江口大壩開始實施加高加固,正常蓄水水位由157 m提升至170 m,庫容由1.745×1010m3擴大至2.905×1010m3,至2014年底南水北調中線工程正式通水。南水北調中線工程實施及建成后,優化了我國水資源配置,有力支撐了水資源的社會經濟發展,但對于丹江口庫區的生態質量、土地資源可持續利用及社會經濟發展都產生了重要影響。因此,以丹江流域(河南段)為例,基于長時序遙感數據,研究土地利用/覆被變化軌跡,揭示土地利用時空演變規律,為研究區土地可持續利用、生態環境建設以及水源地保護提供科學依據。
采用單一遙感數據源往往難以支持長時間序列對土地變化軌跡分析的需求,為獲得長時序多時相土地利用/覆被信息,選用2002—2017年共16期遙感影像數據,分別為Landsat TM(2002—2011年)、HJ-1A CCD(2012年)、Landsat OLI(2013—2017年),其中選用Landsat TM影像的1~5波段和第7波段、HJ-1A CCD相機波段和Landsat OLI影像的1~7波段,從而保證所用影像的分辨率均為30 m,其他數據包括研究區2015年土地利用變更數據庫和道路水系等。
為保證長時序遙感影像對空間疊加及土地利用分類的準確性,需要對影像進行輻射定標、校正、波段組合等預處理[17],從而減少光照、大氣等對不同時相影像造成的誤差,同時也可以提高土地利用/覆被信息的解譯精度。本研究在對長時序遙感影像進行預處理后,采用統計特征分析、主成分分析和相關分析等方法,利用相關系數矩陣、特征向量和特征值,按照最佳波段組合原則[18],得到適合研究區土地利用分類的最佳波段組合:Landsat TM(5、4、3波段)、HJ-1A CCD(1、2、3波段)、Landsat OLI(7、5、4波段)。
根據土地利用一級類型劃分體系和丹江流域(河南段)實際土地利用特征,將土地利用類型劃分為水域、裸地、耕地、林地、草地、建設用地6類。結合土地利用變更數據庫、Google歷史影像及實地勘探,每期影像均選擇兩套樣本點,分別用于分類和精度驗證[19]。采用支持向量機分類和人工目視解譯法,獲取研究區16年間土地利用/覆被分類結果,如圖2所示。經驗證,分類總體精度不小于90.43%,Kappa系數達到了0.84,滿足分類精度的要求。影像數據分類完成后,對地類屬性做重分類處理,將水域、裸地、耕地、林地、草地、建設用地依次重編碼為1、2、3、4、5、6,為后續研究區多時相土地利用變化圖譜的提取提供數據準備。

圖2 2002—2017年研究區土地利用分類圖Fig.2 Land use status maps of study area
土地利用變化軌跡的建立等同于多時相遙感數據土地空間單元變化圖譜的提取,史培軍等[1]提出了土地利用變化空間表達方法,變化圖譜提取公式為
Ci×j=10nAki×j+10n-1A(k+1)i×j+10n-2A(k+2)i×j+…+ 10n-nA(k+n)i×j
(1)
式中Aki×j——k時刻研究區土地利用圖譜
i——研究區行數
j——研究區列數
n——時相數
Ci×j——n+1個時相的土地利用變化圖譜
變化圖譜表現了由k時刻到k+n時刻的土地利用變化類型和空間分布,該方法適用于土地利用類型小于10的情況。
較多學者均采用土地利用變化空間表達方法,并結合ArcGIS技術的柵格運算功能,構建長時序相對“連續”的土地利用/覆被變化軌跡,時相數均為4~6個[20-22]。在本研究中,需要構建丹江流域(河南段)2002—2017年共16個時相的土地利用變化軌跡,變化圖譜在影像柵格中表示為16位的屬性值,例如某空間單元初始時相的土地利用類型為水域,并從始至終未發生變化,則其變化圖譜的柵格屬性值為1 111 111 111 111 111。但由于柵格運算功能只能使用系統字段進行計算,該字段值類型為長整型,取值范圍為-2 147 483 648~2 147 483 648,最多只能滿足10個時相的土地利用變化圖譜的表達。因此,在本研究中將16個時相的土地利用分類結果通過轉換工具轉換為矢量,使得每一個土地空間單元對應一個矢量點,采用空間連接和字段計算器功能,代入式(1),再將矢量點轉換為柵格,此時柵格數據中參與運算的字段類型為字符串,字段值即為研究區16個時相的土地利用變化圖譜信息。
AMY等[23]基于6時相23個隨機地點的土地利用變化結果,通過計算研究區長時序多時相變化圖譜中的STD指標,即研究時段內各土地空間單元所經歷演變過程的變化次數(Turnover,T)、相似性(Similarity,S)和多樣性(Diversity,D)來綜合反映土地利用軌跡變化過程,從而提出了土地利用/覆被長時序多時相穩定映射變化軌跡分析方法。其中,S表示某土地空間單元在研究時段內經歷的相同土地利用類型的數量;T表示研究時段內該土地空間單元在相鄰時間段發生變化的次數;D表示某土地空間單元在研究時段內經歷的不同土地利用類型的數量。
STD指標是穩定映射變化軌跡模型的基礎,是劃分研究區土地利用變化軌跡類型的關鍵。參考AMY等[23]提出的土地利用變化軌跡的STD判定原則,本研究結合長時序遙感數據變化圖譜信息,推導計算出穩定映射STD指標與時相數的關系式,并根據研究區的實際情況,提出適合研究區的基于改進穩定映射法的LUCC變化軌跡分析方法。
2.3.1STD指標計算
土地利用變化軌跡的種類受研究時相數和土地利用類型數量的影響,令研究時相數為n,土地利用類型數為m,則可能出現的變化軌跡有mn種。本研究參考土地利用類型一級劃分體系,將研究區劃分為6種土地利用類型,選用連續的16個時相的土地利用結果,則可能出現的變化軌跡有616種,在應用穩定映射變化軌跡STD判定方法時,各土地空間單元演變過程中經歷變化次數、多樣性和相似性3個指標的計算更加復雜。由穩定映射STD指標的原理可知,長時序土地利用變化軌跡中,研究時段內某土地空間單元在相鄰時間段發生變化的次數決定著該土地空間單元所經歷的相同和不同土地利用類型的數量,即其所經歷演變過程的變化次數T決定著該條變化軌跡上土地利用類型的相似性和多樣性。而隨著研究時段內時相數的增多,其變化軌跡中所經歷演變過程的變化次數也越來越多,則土地利用類型就具有更高的相似性和多樣性。本研究發現,時相數量的變化與變化次數、多樣性、相似性3個指標的計算具有一定的對應關系,由此推導出在土地利用類型為6的情況下,時相數與STD指標的計算關系式(表1),并通過隨機抽樣驗證了其正確性。表1中,Ceil為向右取整函數。

表1 長時序土地利用變化軌跡的STD指標Tab.1 STD index for long time series land use change track
將本研究中n=16代入表1的范圍中,可以計算出研究區土地利用變化過程中,經歷不同變化次數的土地利用變化軌跡所包含的地類多樣性和相似性,從而得到研究區16個時相的土地利用穩定映射STD指標(表2),為下一步土地利用變化軌跡類型劃分提供數據基礎。

表2 研究區16個時相的土地利用變化STD指標Tab.2 Changes in soil utilization STD index at 16 in study area
2.3.2改進的土地利用變化軌跡分析方法
結合丹江流域(河南段)實際土地利用變化情況,本研究提出了一種基于長時序遙感的土地利用變化軌跡的STD判定方法(表3)。將一級變化軌跡類型劃分為穩定型、漸變型、非連續漸變型、循環型和波動型5類。由于研究時序較長,研究區穩定土地空間單元會存在因土地利用規劃占地或分類誤差等影響而在某一年產生突變,但卻長期保持該種土地利用類型的情況,因此本研究將穩定型土地利用變化軌跡劃分為持續穩定型和長期保持型兩類二級過程;同時,針對漸變型土地利用變化軌跡,又根據變化發生時期的早晚將其劃分為早期變化后期穩定、中期變化和早期穩定后期變化3類二級過程;其次,根據研究區實際土地演變情況和STD指數,將循環型土地利用類型又劃分為AB類循環和ABC類循環。最后,根據土地利用類型,將持續穩定型和長期保持型分別劃分為持續穩定林地、持續穩定草地、持續穩定水域、持續穩定耕地、持續穩定裸地、持續穩定建設用地和長期保持林地、長期保持草地、長期保持水域、長期保持耕地、長期保持裸地、長期保持建設用地6類三級過程。
傳統的穩定映射變化軌跡分析方法是根據研究時段內區域土地利用變化的軌跡類型,分析其不同軌跡類型的面積變化和分布特征。近年來,國內外學者在分析其軌跡類型的基礎上,結合定量分析模型和景觀指數進行了更深層次的土地利用變化軌跡的演變特征分析。定量分析模型可以反映研究時段內區域的單一土地利用類型和綜合土地利用類型的變化量、變化幅度和變化狀態等特征,景觀指數可以計算景觀單元時空變化并反映景觀格局結構等空間分布特征,其中景觀水平上的景觀格局分析能夠反映出整個研究區的景觀格局的演變規律[24-26]。本文將研究區土地利用變化軌跡劃分為5種軌跡類型,結合其土地利用軌跡類型特征,選用單一土地利用動態度和綜合動態度兩個定量模型以及景觀水平上的面積(CA)、面積占比(PLAND)、斑塊密度(PD)、最大斑塊指數(LPI)、散布與并列指數(IJI)和聚集度(AI)6種常用的景觀指數,將變化軌跡相同的斑塊看作同一種景觀類型,利用ENVI軟件和Fragstats軟件,對研究區土地利用變化軌跡類型從動態變化總體特征和時空變化規律兩方面進行系統分析。

表3 研究區土地利用變化軌跡的STD判定方法Tab.3 STD determining method of change track of soil use in study area
在2002—2017年間,研究區各土地利用類型的組成結構和面積均發生了顯著變化。從土地利用類型的組成結構來看,研究區地勢從西北到東南依次為山地、丘陵、壟崗和平原,研究區的土地利用狀況主要以林地和耕地為主;從土地利用類型的面積變化來看(表4),隨著研究區可持續發展進程加快以及南水北調中線工程的實施,2002—2017年間,研究區土地利用變化總體趨勢主要以耕地向建設用地轉化以及水域面積不斷增加為主。具體表現為:①耕地面積在2004—2011年間持續減少,但2004年耕地面積較2003年增長約453.40 km2,主要原因為1999—2003年間,城鎮、市區的迅速擴張和鄉鎮企業的突出發展占用了大量耕地,國家在該階段相繼出臺了一系列保護耕地的法律法規,通過土地開發、復墾、整理等措施增加耕地面積[27]。2004年后,耕地面積逐年穩定減少。②林草地面積變化趨勢較為復雜,2002—2004年呈現下降趨勢,但在2006年后林草地面積有小幅度上升,主要因為該階段“退耕還林還草”政策對其有較大的影響。 ③建設用地面積穩步增加,由2002年的106.35 km2增加至2017年的416.68 km2,共增加了310.33 km2。④水域面積總體呈上升趨勢,原因在于2003年底南水北調中線工程開始實施,為使調水工程順利進行,2005年開始對丹江口大壩進行加高加固處理,至2014年庫區面積有了較為明顯的增加。在2010—2012年期間,水域面積有所減少,與該時期丹江上游季節性降水和丹江口的人工調蓄有一定關系[28-29]。

表4 研究區土地利用類型面積統計Tab.4 Statistics of land use type area in study area km2
對比連續16年間研究區各土地利用類型的年變化率和綜合土地利用動態度(圖3)可以看出:①2002—2008年各土地利用類型的年變化率較大,這期間內研究區的土地利用類型之間的轉變尤為頻繁。②建設用地在研究時段內基本保持正向增長趨勢,最大年增長率為51.00%,而耕地基本保持負向減少。③年變化率的計算受前一年度土地利用類型面積的影響,因此可以看到2005—2006年草地的年增長率為127.00%,主要原因為2003—2005年出臺的一系列“復墾耕地”政策導致草地面積減少。④綜合土地利用空間動態度曲線定量描述了研究區土地利用類型的變化幅度和活躍程度。由圖3可以看出,研究初期丹江流域(河南段)各土地利用類型轉換頻繁,其中2008年土地利用變化最為劇烈;2009—2015年研究區綜合土地利用動態度均在8%以下,其土地利用類型變化幅度較小;2016年后綜合土地利用動態度有明顯增加,表明這兩年研究區土地利用類型變化較為活躍。

圖3 研究區2002—2017年土地利用/覆被定量 分析折線圖Fig.3 Quantitative analysis of land use/covers in study area in 2002—2017
基于改進穩定映射法的變化軌跡分析方法,按照土地利用變化軌跡動態過程將研究區變化軌跡劃分為穩定型、漸變型、非連續漸變型、循環型和波動型5個一級類,并根據研究區實際土地利用情況劃分二級和三級軌跡類型,然后結合6個常用景觀指數,從空間構型上對研究區土地利用變化軌跡作深入分析(表5)。其構成特征如下:
(1)穩定型
此類軌跡類型表示自研究初期至研究末期,研究區土地利用類型從未發生過變化,或在研究期內長期保持某種地類的圖譜單元。本研究中穩定型面積占比為61.80%。此種軌跡類型可以細分為持續穩定型和長期保持型。
持續穩定過程表示在連續16年間一直保持同一種土地利用類型的演變過程。結合圖4a可以看出,在該種軌跡類型中,持續穩定林地面積最大,占研究區總面積的47.85%,主要分布于研究區北部和西北部的山地丘陵地帶,其次為持續穩定耕地和持續穩定水域,面積占比分別為4.85%和2.03%,分別位于研究區東南部的平原地帶和南部的丹江口庫區。同時,分析最大斑塊指數(LPI)和聚集度(AI)可知,持續穩定林地的最大斑塊景觀占比最高,為14.01%,聚集度為94.02%,說明持續穩定林地的景觀破碎度較低,斑塊聚集程度較高,林地的土地利用狀況較為穩定;草地的最大斑塊占比最小,趨近于0,聚集度為32.45%,說明該種地類景觀破碎度高,主要以細碎斑塊為主,土地利用變化活躍。裸地為過渡性土地利用類型,其面積占比很小,持續穩定裸地和長期保持裸地類型占比趨近于0,因此不對其進行景觀指數分析。
長期保持過程是指某土地空間單元在研究期內土地利用類型原則上未發生變化,但由于研究時相較多,可能因為土地利用規劃占地或分類結果存在誤差等原因導致某一時刻地類發生一次變化,但仍長期保持該種地類的演變過程。結合圖4b,此類軌跡斑塊面積占研究區總面積的6.47%,主要集中在持續穩定型軌跡斑塊附近,并在持續穩定型斑塊周圍離散分布。因此,長期保持型軌跡的AI值和整體水平AI值都偏低。
(2)漸變型
漸變型土地利用軌跡變化過程表示在研究區內土地利用類型只發生一次變化,即從一類轉換為另一類的過程。該類型圖譜單元變化次數為1,多樣性為2,占研究區總面積的7.32%。根據變化時期的早晚,又將2002—2007年期間發生的變化歸為早期變化后期穩定,2008—2011年期間發生的變化歸為中期變化,將2012—2017年期間發生的變化歸為早期穩定后期變化。結合圖4c和表6可以看出,早期變化后期穩定的面積占比為3.86%,主要分布在研究區中部、北部和西南部,轉化為林地、草地的面積較多;中期變化的面積占比為0.45%,主要轉化為草地和建設用地;早期穩定后期變化的面積占比為3.01%,主要分布在研究區東南部平原地帶和丹江流域支流的流經區域,對應的土地利用流向分別為耕地向草地、建設用地和耕地向水域的轉化。

表5 研究區土地利用變化軌跡景觀指數統計Tab.5 Land-use change track landscape index statistics

圖4 研究區土地利用變化軌跡的空間分布Fig.4 Spatial distribution maps of land-use change trajectories

表6 漸變型土地利用類型轉化面積統計Tab.6 Land use type conversion area statistics km2
(3)非連續漸變型
非連續漸變型是表示2種或3種土地利用類型的周期性轉化,但趨勢不明顯的土地利用軌跡過程,占研究區總面積的8.86%。本研究提取出面積占比超過0.30%的變化軌跡如圖4d所示,其中非連續耕地林地漸變和非連續耕地草地漸變共占此類型面積的87.33%。結合表5中非連續漸變型AI值為18.66%,說明該種軌跡類型由較多細碎小斑塊組成,各個地類之間轉換頻繁。
(4)循環型
循環型軌跡變化過程是描述2種或3種土地利用類型間的相互轉化,由于研究時序較長,時相較多,結合研究區實際土地利用狀況,將該軌跡類型細分為兩種地類間的循環轉化和3種地類間的循環轉化。經統計發現,耕地草地循環、耕地水域循環和林地草地循環的面積在該類循環中面積占比為83.50%,因此主要針對這3種循環轉化進行分析(圖4e)。其中耕地草地循環面積占比為3.79%,從空間分布來看,主要分布在穩定型耕地和穩定型草地的交錯區域;耕地水域循環主要分布在淅川縣丹江口庫區支流的徑流區域,面積占比為0.49%;林地草地循環面積占比為2.06%,受季節和區域影響,主要分布在研究區山地的山脊、山谷以及丘陵和平原相接的區域。
(5)波動型
波動型土地利用軌跡變化過程具有更高的多樣性,表達不同地類間的頻繁轉化,占研究區總面積的14.43%。同時,因為波動型軌跡類型較多,本研究中僅將面積占比大于0.03%的軌跡類型選取出來。在此類軌跡變化過程中,圖譜單元多樣性為3的軌跡面積占該類總面積的64.20%,即在波動型的軌跡變化過程中,主要是3種地類間的頻繁轉化。其中,涉及到耕地和草地的變化軌跡居多,說明在研究區連續16年的土地利用變化中,耕地草地的變化是較為活躍的。
綜上所述,在2002—2017年期間,研究區土地利用流向主要是耕地轉化為建設用地;林草轉化和耕草轉化較為頻繁;同時,水域面積明顯增大,與南水北調中線工程的實施和丹江口庫區大壩的增高緊密相關。從研究區土地利用變化總體特征和時空變化規律可以看出,16年來土地利用變化總體趨勢主要以耕地向建設用地轉化和水域面積不斷增加為主,其中耕地面積由2002年的2 474.54 km2減少至2017年的1 562.60 km2,建設用地面積和水域面積分別由2002年的106.35 km2和218.60 km2增加至2017年的416.68 km2和400.31 km2,林草地在16年間面積變化較為復雜,呈現先減少后上升的趨勢,主要受該階段“退耕還林還草”政策的影響。綜合分析16年間土地利用變化率可以得出,建設用地在研究時段內基本保持正向增長,這與16年間社會經濟的快速發展和城鎮、市區的迅速擴張有密切關系,同時2008年研究區綜合土地利用動態度呈增加趨勢,表明該年研究區土地利用類型的變化最為劇烈。研究區2002—2017年間穩定型軌跡面積占比最大,主要分布在研究區北部山區、丘陵地帶林地、丹江口庫區和東南部耕地等區域,達到61.80%,說明研究區總體土地利用狀況穩定,符合研究區土地利用實際情況。而地類轉化頻繁的區域主要位于丘陵與平原的過渡地帶,林地、耕地和草地出現了顯著的相互轉化現象。
提取土地利用/覆被變化軌跡,以研究土地利用類型的轉移和流向,揭示土地利用變化規律,從而分析時間序列中土地利用/覆被的時空演變特征[30-31]。但隨著研究時間序列數據的增多,土地利用/覆被變化軌跡更加繁雜,軌跡提取更加困難。本研究以丹江流域(河南段)為研究區,基于長時序遙感數據,提取土地利用變化圖譜信息,并以改進的穩定映射變化軌跡分析方法將變化軌跡重分類,結合土地利用動態度等定量模型和景觀格局指數,揭示研究區2002—2017年間土地利用類型的變化軌跡和時空演變規律,研究結果為土地資源可持續利用和水源地生態環境保護提供科學依據,同時也為南水北調工程后期進一步擴大引漢規模,并預計在2030年完成年均調水量達130億m3的目標提供數據參考與決策支持。
(1)在2002—2017年間,研究區穩定型軌跡面積占比最大,達到61.80%,說明研究區總體土地利用狀況穩定。
(2)隨著城市化進程的加快和丹江口庫區的建設,16年間研究區土地利用變化總體趨勢以耕地向建設用地轉化和水域面積不斷增加為主。由于丹江口庫區建設和南水北調中線工程的實施,水域面積由218.60 km2增長到400.31 km2。
(3)受地形影響,研究區北部山區和丘陵地帶林地穩定性較好,而丘陵與平原的過渡地帶林地和耕地、草地和耕地間的交替性變化頻繁,占研究區總面積的5.85%。隨著南水北調中線工程的實施,研究流域內土地利用變化軌跡對于制定區域土地利用政策及水源地保護意義重大。