基于1990—2019年多時相影像的干旱區(qū)綠洲景觀格局分析

宋 奇， 史 舟， 馮春暉， 馬自強， 紀文君，彭 杰， 高 琪， 蔣學瑋

（1.塔里木大學園藝與林學學院，新疆 阿拉爾 843300；2.浙江大學環(huán)境與資源學院，浙江 杭州 310058；3.塔里木大學農(nóng)學院，新疆 阿拉爾 843300；4.北京大學地球與空間科學學院遙感與地理信息系統(tǒng)研究所，北京 100871；5.中國農(nóng)業(yè)大學土地科學與技術(shù)學院，北京 100193；6.新疆昌吉州地質(zhì)環(huán)境監(jiān)測站，新疆 昌吉 831100）

全球變暖、冰川融化、生物多樣性降低等生態(tài)環(huán)境問題日益突出，加之頻繁的人類活動，打破了正常的生態(tài)平衡，區(qū)域生態(tài)環(huán)境逐漸惡化［1-3］，對世界景觀格局造成不小的影響［4-6］。為探明全球變暖環(huán)境下生態(tài)環(huán)境所受的影響，有必要對區(qū)域景觀格局進行分析，探索人類活動對景觀生態(tài)系統(tǒng)的影響［7-9］。對景觀格局信息的提取和分析不僅能夠表征空間的復(fù)雜性而且也決定著生態(tài)系統(tǒng)中資源環(huán)境的分布情況。對景觀格局進行研究，能在無序的景觀斑塊鑲嵌中發(fā)現(xiàn)潛在的生態(tài)學意義和景觀規(guī)律［10-11］。其中，區(qū)域景觀分類是分析景觀格局變化的基礎(chǔ)［12］，傳統(tǒng)地面調(diào)查方法不僅費時費力而且適用范圍較小。隨著遙感技術(shù)的應(yīng)用，遙感影像數(shù)據(jù)以其低成本、獲取方式便捷、覆蓋范圍廣、具有長時間序列等諸多優(yōu)勢，成為景觀格局變化分析中首選的數(shù)據(jù)源［13-15］。

已有相關(guān)研究［16-20］表明，目前景觀格局研究的遙感數(shù)據(jù)源多為典型的時間斷面數(shù)據(jù)，以基于像元的分類方法對長時序的區(qū)域景觀格局變化進行分析。然而，景觀格局的變化很多情況下并不是在關(guān)鍵年份發(fā)生的，典型的時間斷面數(shù)據(jù)很容易遺漏關(guān)鍵節(jié)點信息，難以把握準確的景觀格局變化情況，對連續(xù)性信息的把握不足［21］。此外，基于像元的分類方法（例如支持向量機分類方法）容易出現(xiàn)碎斑，分類后處理相對繁瑣，容易導(dǎo)致分類后精度偏低，最終影響區(qū)域景觀格局的分析。

阿拉爾墾區(qū)從地理方位來看，南抵塔克拉瑪干沙漠，北靠天山山脈，具有不可替代的景觀生態(tài)屏障功能。墾區(qū)受人類活動影響較大，區(qū)域內(nèi)的植被以人工植被為主，自然植被很少，景觀結(jié)構(gòu)單一，生態(tài)環(huán)境較為脆弱，景觀格局變化較大。因此，本研究選擇西北地區(qū)典型人工綠洲阿拉爾墾區(qū)為例，對長時序區(qū)域景觀格局變化信息進行分析，能為西北干旱地區(qū)相關(guān)研究提供一定的參考價值。

本研究獲取阿拉爾墾區(qū)近30 a 所有可用的Landsat 遙感影像，共226 景。對其逐年進行NDVI最大值合成，由此得到連續(xù)30 a 長時間序列的數(shù)據(jù)源，并對比多時相面向?qū)ο?、多時相基于像元和單時相面向?qū)ο? 種分類方法的精度，選取本區(qū)域精度最高的分類方法進行景觀分類，進而得到1990—2019年間阿拉爾墾區(qū)的景觀分類結(jié)果，以此分類結(jié)果分別從斑塊類型、景觀水平和空間分布3 方面對研究區(qū)景觀格局信息進行深入分析，并通過累積距平法探明區(qū)域景觀格局的突變信息，并檢驗了其突變點的真?zhèn)?。本文以連續(xù)長時間序列和NDVI最大值合成為創(chuàng)新點進行區(qū)域景觀格局分析，能夠從數(shù)據(jù)源開始就提高景觀格局信息的精確度，以此保證后續(xù)景觀格局分析的可靠性。

1 研究區(qū)域與數(shù)據(jù)來源

1.1 研究區(qū)概況

阿拉爾墾區(qū)位于阿克蘇河、和田河及葉爾羌河的三河交匯處，地理位置為40°22′～40°57′N，80°30′～81°58′E（圖1），總面積4105.92 km2，地勢呈西北高東南低走向，平均海拔1012 m，生態(tài)系統(tǒng)屬于荒漠-綠洲型，土壤類型主要為砂壤土，水資源分布主要以塔里木河及勝利、多浪和上游三大平原水庫為主。研究區(qū)地處歐亞大陸腹地，塔克拉瑪干沙漠北緣，屬暖溫帶極端大陸性干旱荒漠氣候，氣象多變，常有沙暴、揚塵和霜凍等自然災(zāi)害，全年降水少而蒸發(fā)大。阿拉爾墾區(qū)以棉花和紅棗種植為主，是中國最大的長絨棉種植基地［22］。此外，有小面積種植香梨、蘋果、葡萄、玉米、小麥、高粱等作物。

圖1 阿拉爾墾區(qū)位置示意圖Fig.1 Location map of Alar Reclamation Area

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

獲取1990—2019年南疆阿拉爾墾區(qū)（軌道號為146—32）空間分辨率為30 m 的Landsat 5、7 和8 所有遙感影像，為提高后期遙感影像分類的判讀精度，從所有影像中篩選出226 景云量低于40%的遙感影像作為本文基礎(chǔ)數(shù)據(jù)，其中各年遙感影像的云覆蓋率情況及影像數(shù)如圖2所示。

圖2 遙感數(shù)據(jù)Fig.2 Remote sensing data

2 研究方法

2.1 NDVI最大值合成

首先獲取1990—2019 年所有Landsat 遙感影像，從中篩選出226 景可用影像，在ENVI 軟件中將本文獲取的所有遙感影像（共226景）逐一進行輻射定標、大氣校正、圖像配準、裁剪等影像預(yù)處理操作，進而逐一提取歸一化植被指數(shù)（Normalized Difference Vegetation Index,NDVI），其計算公式為：

式中：NIR為Landsat近紅外波段反射率；R為紅光波段反射率。

相關(guān)合成的流程如圖3 所示，得到30 a 間連續(xù)時間序列的NDVI合成圖（圖4）。從NDVI合成圖中可看出，多時相遙感影像合成圖不僅消除了云量的影響，而且圖像中的植被覆蓋度更高，有利于后期景觀分類。

圖3 技術(shù)流程圖Fig.3 Technique flow chart

圖4 NDVI最大值合成圖Fig.4 Maximum NDVI composite

將1990—2019年各年份的NDVI年際變化信息（圖5）進行匯總，可以明顯發(fā)現(xiàn)近30 a 間，阿拉爾墾區(qū)的NDVI 處于增加的趨勢，從1990 年的0.1661增加到2019 年的0.4319，植被覆蓋度有明顯的增加趨勢。

圖5 1990—2019年阿拉爾墾區(qū)NDVI年際變化Fig.5 Annual variation of NDVI in Alar Reclamation Area from 1990 to 2019

2.2 分類方法

參考《新疆生產(chǎn)建設(shè)兵團統(tǒng)計年鑒》（1990—2019）中的土地利用情況、研究區(qū)地表分布特點和新疆已有案例［23-24］中的分類方法，選用面向?qū)ο蠓诸惙椒ㄟM行景觀分類。該分類方法是綜合考慮各景觀類型的光譜、形狀和灰度值等特征，利用eCognition軟件對多時相遙感影像進行多尺度分割，經(jīng)過多次測試，通過ESP（Estimation of Scale Parameter）評價工具找到最佳分割尺度，選用面向?qū)ο蟮碾S機森林（Random Forest）分類算法對其進行信息提取。本區(qū)域的耕地、林草地、園地和非植被的NDVI值不同，由此對其進行分類；此外，水體、建設(shè)用地和未利用地的幾何和紋理特征區(qū)別較大，以此作為這3類用地的分類依據(jù)，構(gòu)建合理的影像分類規(guī)則，進行景觀分類。

2.3 景觀格局變化分析方法

參考已有案例［25］，分別從斑塊類型水平和景觀水平對阿拉爾墾區(qū)景觀格局變化特征進行定量分析。選用常用的景觀指數(shù)［26］進行分析，分別從景觀破碎度、優(yōu)勢度、復(fù)雜度、聚合度、連通性和多樣性方面體現(xiàn)研究區(qū)的景觀格局變化情況。由于本研究的時間序列較長，在長時序的景觀格局變化中很可能會出現(xiàn)突變現(xiàn)象，我們選用累積距平法［27］探明研究區(qū)的景觀格局突變情況，計算公式如下：

式中：LRi為第i年的距平累積值；Ri為第i年的景觀指數(shù)值；為景觀指數(shù)序列的多年平均值。

為了驗證所得突變點的可靠性，筆者對累積距平法所得的突變點進行驗證［28］，保證所得結(jié)果的可靠性；最后對景觀指數(shù)進行可視化分析，選用不同大小的移動窗口，所得結(jié)果不同，因此在對比不同的移動窗口后，篩選出最佳窗口來反映景觀空間格局的空間變化特征。最終從景觀破碎度、優(yōu)勢度、復(fù)雜度、聚合度、連通性、多樣性、突變點和空間分布等方面綜合分析阿拉爾墾區(qū)景觀格局變化特征。

3 結(jié)果與分析

3.1 不同分類方法的分類結(jié)果及其精度評價

以2019年為例，參考高分遙感影像并對研究區(qū)進行實地調(diào)查后，在基于像元的分類方法中共選取723951 個分類樣本，其中耕地有173940 個，林草地有131591 個，園地有135138 個，水體有139291 個，建設(shè)用地有17488 個，未利用地有126503 個，進而得到相對應(yīng)的分類后結(jié)果，最后對3 種分類方法所得結(jié)果進行對比（圖6）。從圖中可以明顯看出，多時相面向?qū)ο蠓诸惙椒ǖ姆诸惤Y(jié)果最佳，而多時相基于像元和單時相面向?qū)ο蟮姆诸惤Y(jié)果較為破碎，容易出現(xiàn)“錯分”和“漏分”現(xiàn)象。

圖6 3種分類方法的分類后結(jié)果比較Fig.6 The results of three classification methods were compared

隨機選取阿拉爾墾區(qū)351173個驗證樣本，選用常用的混淆矩陣方法對分類后結(jié)果進行精度評價。從精度評價結(jié)果中得出多時相面向?qū)ο蠓诸惙椒ǖ目傮w精度和Kappa 系數(shù)最高，分別為96.57%和0.95，說明利用多時相面向?qū)ο蠓诸惙椒ㄟM行景觀分類能夠有效區(qū)分不同的景觀類型，分類結(jié)果可靠。此外，最優(yōu)分類方法多時相面向?qū)ο蠓诸惙椒ǖ目傮w精度相較于多時相基于像元分類方法提高了7%，Kappa系數(shù)相較于多時相基于像元分類方法提高了0.11，各景觀類型的生產(chǎn)者精度均得到不同程度的提高；相較于單時相面向?qū)ο蠓诸惙椒ㄌ岣吡?2.22%，Kappa系數(shù)提高了0.25。

多時相面向?qū)ο蠓椒ǖ脑敿毦仍u價情況如表1 所示。從表中可以看出，該最優(yōu)分類方法的總體精度OA達到了98.79%。同時各類型的用戶精度UA 和制圖精度PA 均在90%以上，說明該分類方法在本研究區(qū)是可行的。但水體和未利用地之間容易混淆，這是因為阿拉爾墾區(qū)水體附近存在小面積濕地，很容易被錯分成為未利用地，但對整個研究區(qū)而言，被錯分的濕地面積較小，不影響整體分類結(jié)果。同理，其他歷史年份均按此方法重新建立樣本和規(guī)則，逐年進行分類，即1990—2019 年共建立30 次分類樣本進行景觀分類。參考新疆現(xiàn)有案例［29-30］中的分類情況對本文的分類結(jié)果進行精度對比，對每年的分類結(jié)果重新建立精度評價混淆矩陣，保證每年分類的用戶精度、制圖精度和總體分類精度均在85%以上。

表1 阿拉爾墾區(qū)用地類型精度評價Tab.1 Evaluation of land type accuracy in Alar Reclamation Area

3.2 斑塊類型水平上的景觀格局分析

圖7為30 a間阿拉爾墾區(qū)在類型水平上的景觀格局指數(shù)變化情況。（1）斑塊密度（PD）：耕地、林草地和園地的斑塊密度較大。30 a 間，耕地類型的斑塊密度數(shù)值呈現(xiàn)出增加的趨勢，說明耕地類型的景觀趨于破碎化；園地類型的斑塊密度數(shù)值在2005年發(fā)生突變，數(shù)值突然增大，這很可能是與2005 年的政策導(dǎo)向有關(guān)，紅棗單價上漲，致使部分耕地轉(zhuǎn)變?yōu)閳@地［31-32］，而水體、建設(shè)用地和未利用地的斑塊密度數(shù)值變化較小，30 a 間沒有明顯的變化趨勢。（2）最大斑塊指數(shù)（LPI）：從最大斑塊指數(shù)的數(shù)值變化情況來看，未利用地類型的最大斑塊指數(shù)明顯大于其他景觀類型，但整體呈現(xiàn)出下降的趨勢；此外，園地類型的數(shù)值整體呈現(xiàn)出上升的趨勢。30 a 間，林草地類型在1994年出現(xiàn)了有史以來的最大值，而2017年的數(shù)值最小，出現(xiàn)最大值和最小值的可能原因是墾區(qū)最開始不斷開荒擴地，后面未利用地的開墾速度降低，耕地的開源主要是林草地的轉(zhuǎn)化。相比之下，水體和建設(shè)用地類型的最大斑塊指數(shù)基本沒有變化。（3）面積加權(quán)平均分維數(shù)（FRAC_AM）：林草地和園地的分維數(shù)明顯高于其他類型，表明林草地和園地的景觀斑塊形狀最復(fù)雜。30 a間，耕地、建設(shè)用地和未利用地的分維數(shù)都有所增加，而水體的分維數(shù)呈現(xiàn)出減少的趨勢。（4）聚集度指數(shù)（COHESION）：園地類型在1990—1999 年的數(shù)值呈現(xiàn)出增加的趨勢，說明此期間園地的聚合性增強。

圖7 阿拉爾墾區(qū)在斑塊類型水平上的景觀格局指數(shù)變化Fig.7 Index change of landscape pattern at patch type level in Alar Reclamation Area

3.3 景觀水平上的景觀格局分析

30 a 間，研究區(qū)在景觀水平上的指數(shù)變化情況如圖8 所示。斑塊個數(shù)（NP）呈現(xiàn)出波動式上升的趨勢，說明區(qū)域內(nèi)的景觀碎斑有所增加；景觀形狀指數(shù)（LSI）同斑塊個數(shù)呈現(xiàn)出相似的變化趨勢，說明研究區(qū)景觀斑塊呈現(xiàn)出一定的復(fù)雜性；聚集度指數(shù)（CONTAG）呈現(xiàn)出明顯的減小趨勢，說明研究區(qū)近年來出現(xiàn)了較多的小斑塊，而大型斑塊的數(shù)量在減少，整體空間分布較為破碎；香農(nóng)多樣性指數(shù)（SHDI）整體上處于增加的趨勢，說明研究區(qū)異質(zhì)性有所增強。

圖8 阿拉爾墾區(qū)在景觀水平上的景觀格局指數(shù)變化Fig.8 Index change of landscape pattern at landscape level in Alar Reclamation Area

圖9為利用累積距平法進行突變點檢驗所得的結(jié)果，從圖9a、圖9b、圖9d中可以看出，以2005年為界，累積距平呈現(xiàn)出先降后增的趨勢；由此說明阿拉爾墾區(qū)在1990—2019年間的景觀破碎度、復(fù)雜度和多樣性呈現(xiàn)出了先降低后增加的趨勢；圖9c的曲線變化情況剛好相反，同樣以2005 年為界，在2005年之前曲線呈現(xiàn)出一定的增加趨勢，在2005年之后呈現(xiàn)出下降的趨勢；由此可以看出，2005 年為阿拉爾墾區(qū)的突變年份，區(qū)域景觀連通性呈先增后降的趨勢。研究區(qū)從2004 年下半年開始大面積荒地被開墾為農(nóng)田，2005 年開始種植農(nóng)作物，農(nóng)作物種植面積急劇增加，從而植被覆蓋度增加，致使阿拉爾墾區(qū)的景觀破碎度、復(fù)雜度和多樣性顯著增加，最終導(dǎo)致阿拉爾墾區(qū)的景觀格局發(fā)生了明顯的突變現(xiàn)象。本研究使用累積距平法對該突變點進行檢驗，驗證其真?zhèn)危员ＷC研究的可靠性。我們將1990—2019年的連續(xù)時間序列數(shù)據(jù)以2005年為界，分為兩個子時間序列，結(jié)合SPSS軟件中的K-S兩樣本檢驗法對兩子時間序列進行差異性檢驗，所得結(jié)果如表2所示，從表中可以看出，各指數(shù)的雙尾顯著性檢驗（sig）均低于0.05，說明兩個子時間序列的差異性顯著，所得結(jié)果可靠。

圖9 阿拉爾墾區(qū)在景觀水平上的突變分析Fig.9 Analysis of abrupt change at landscape level in Alar Reclamation Area

表2 雙樣本K-S檢驗結(jié)果Tab.2 Result of K-S test for two-independent samples

3.4 基于移動窗口的景觀格局分析

在對景觀指數(shù)進行可視化分析中，移動窗口大小的選擇很重要，選用不同的移動窗口所得的景觀格局圖差異顯著，經(jīng)過不斷的調(diào)試，對比在邊長為300 m、600 m、900 m、1200 m、1500 m 和1800 m，共計6 種正方形移動窗下的景觀格局結(jié)果，最終選用1200 m進行可視化分析，所得結(jié)果如圖10所示。

圖10 基于移動窗口的景觀格局分布Fig.10 Landscape pattern distribution based on moving windows

過度的人類活動造成斑塊破碎化，斑塊密度能夠反映出研究區(qū)景觀破碎化程度，對斑塊密度進行分析能夠有效反映人類活動強度和土地利用程度。由于研究區(qū)中部有塔里木河流，使得靠近水源的區(qū)域農(nóng)業(yè)發(fā)展更好，景觀結(jié)構(gòu)更為穩(wěn)定，同時人類活動也更劇烈，從圖10a可以看出，中間塔里木河區(qū)域的斑塊密度更高，這很可能和近年來人類活動不斷增強，工業(yè)用地和農(nóng)田面積增加，然而阿拉爾墾區(qū)的工業(yè)廢水和農(nóng)業(yè)排水多是往偏遠地區(qū)排放，導(dǎo)致偏遠地區(qū)生長出零星的植物組團，植物多樣性增加，這些植被并不是連片分布，而是呈現(xiàn)出離散的空間分布趨勢，最終偏遠地方的斑塊數(shù)量增加。

同理對香農(nóng)多樣性指數(shù)進行空間可視化分析，所得結(jié)果如圖10b 所示，在中部塔里木河區(qū)域香農(nóng)多樣性指數(shù)的值會更高一些，整體呈現(xiàn)出以塔里木河為中心向四周不斷減弱的趨勢。這和中部人類活動更頻繁，土地利用結(jié)構(gòu)更穩(wěn)固有關(guān)，周邊區(qū)域多為未利用地，景觀結(jié)構(gòu)單一，開發(fā)空間充足。

4 討論

大多數(shù)研究使用的遙感影像數(shù)據(jù)往往只有3～6景，這種分析只能體現(xiàn)景觀格局的大致演變趨勢，不能在長時間尺度下得到詳細的景觀格局演變特征，容易遺漏連續(xù)長時間序列中的關(guān)鍵“拐點”信息［33-36］。本文對研究區(qū)近30 a 的所有Landsat 遙感影像數(shù)據(jù)進行了無差異下載，共使用226 景遙感數(shù)據(jù)參與研究區(qū)景觀格局演變分析，并將同一年中各月份影像進行NDVI最大值合成，得到周年NDVI最大值合成的多時相影像。這種基于連續(xù)長時間序列的多時相遙感影像進行的景觀格局分析，不僅能有效的彌補了以往基于3～6 景影像數(shù)據(jù)分析的不足，充分反映出研究區(qū)近30 a來景觀格局的細部變化特征，而且能夠避免由于不同作物生長旺季不同造成的“錯分”現(xiàn)象［37-42］。

此外，在人為干擾影響較大的年份，如2005年，阿拉爾墾區(qū)景觀格局發(fā)生了明顯的突變現(xiàn)象，如果僅用少量的遙感影像數(shù)據(jù)進行景觀格局分析，則不能有效捕捉研究區(qū)景觀格局的細部變化特征。因此，在進行景觀格局變化分析時，影像數(shù)據(jù)的時相選擇至關(guān)重要，直接關(guān)系到研究結(jié)果的準確性和可靠性。

在景觀格局變化的驅(qū)動機制方面：（1）滴灌技術(shù)的普及。2005 年以后，滴灌技術(shù)逐漸普及，不僅使得原先無法墾殖、地勢較高的土地可以墾殖了，而且地勢較低的、鹽堿重的區(qū)域也能墾殖，致使綠洲耕地面積不斷增加，未利用地面積不斷減少，景觀豐富度和異質(zhì)性逐漸增加。（2）種植業(yè)結(jié)構(gòu)的調(diào)整。如棉田向林果業(yè)調(diào)整對綠洲景觀產(chǎn)生了一定的影響，使得園地的斑塊密度、最大斑塊指數(shù)、分維數(shù)、聚集度指數(shù)增加，致使綠洲景觀的空間格局分布不均勻，景觀聚合性增強。（3）農(nóng)業(yè)技術(shù)進步對景觀格局的影響。技術(shù)的進步（滴灌普及）使得原本少雨、高溫、缺水的干旱區(qū)域大大提高了地下水的利用率，進而影響到土壤和植被，致使研究區(qū)景觀豐富度和異質(zhì)性增加，對綠洲景觀的景觀格局造成了一定的影響。

景觀指數(shù)分析的客觀性：結(jié)合研究區(qū)的實際情況以及研究重點，例如，本文對景觀斑塊數(shù)及其相關(guān)的斑塊密度進行分析時，著重對研究區(qū)所占面積比例較大的耕地、林草地和園地進行了景觀格局分析，并對其三者之間的相互轉(zhuǎn)化進行了解釋說明，而研究區(qū)中面積較小的水體和建設(shè)用地進行了簡要介紹。因此，景觀指數(shù)分析存在一定的客觀性，不同的研究者分析，其結(jié)果差異顯著。

5 結(jié)論

本研究對阿拉爾墾區(qū)近30 a間的景觀格局變化進行分析研究，主要結(jié)論如下：

（1）對多時相面向?qū)ο?、多時相基于像元和單時相面向?qū)ο? 種分類方法的比較分析表明，多時相面向?qū)ο蠓椒ǖ姆诸惥茸罡?，總體精度為96.57%，Kappa 系數(shù)為0.95，較后兩者分別提高了7%、0.11和12.22%、0.25，有效避免了“椒鹽”現(xiàn)象以及植被物候差異造成的“漏分”和“錯分”。

（2）從斑塊類型水平上看，1990—2019 年阿拉爾墾區(qū)耕地景觀趨于破碎化，未利用地優(yōu)勢度降低，園地和林草地的景觀形狀趨于復(fù)雜化，水體和建設(shè)用地的景觀趨于均衡化；從景觀水平上看，研究區(qū)景觀形狀趨于復(fù)雜化和破碎化，景觀連通性降低，景觀豐富度和異質(zhì)性增加；墾區(qū)景觀格局于2005 年發(fā)生突變，景觀破碎度、復(fù)雜度和多樣性顯著增加，景觀連通性降低。

（3）從空間分布方面看，塔里木河沿岸區(qū)域景觀指數(shù)空間分布較為均勻，人類活動頻繁，景觀結(jié)構(gòu)完整，景觀格局呈現(xiàn)出以塔里木河沿岸區(qū)域為核心向偏遠區(qū)域蔓延的趨勢，其中向東北和西南方向的擴張趨勢更為明顯。