2005-2018年黃河三角洲景觀格局脆弱性的時空變化規律

宋雨桐, 張子璇, 牛蓓蓓, 李新舉

(1.山東農業大學 資源與環境學院, 山東 泰安 271018; 2.土肥資源高效利用國家工程實驗室, 山東 泰安 271018)

景觀格局變化是指由于土地利用方式不同導致的景觀類型和景觀格局的時空變化[1]。早期景觀格局的主要研究方法包括傳統的定性描述法、景觀生態圖疊置分析方法等。自20世紀90年代以來，隨著經濟社會的不斷發展，人類活動強度增大，土地開發與土地利用方式的轉變都對景觀格局研究提出了新要求，同時3S技術的快速發展和應用，使得景觀生態學的研究重點開始轉向景觀格局脆弱性及其動態變化的研究。景觀格局脆弱性研究起源于生態脆弱性研究，1905年由美國學者Clements提出，隨著相關研究的進一步深入，生態脆弱性的內涵也發生了變化，從單一研究方向轉向多元化，由此衍生出景觀格局脆弱性的研究[2]。景觀格局脆弱性反映景觀生態系統所受到的干擾及其脆弱程度[3]。國內外學者從多個方面對不同區域景觀格局脆弱性展開了研究。主要從景觀脆弱度指數[4]、土地利用類型變化[5]、地形因子[6]、脆弱性評價體系構建等[7]角度開展。目前針對黃河三角洲景觀格局脆弱性的國內外研究較少，國外方面主要從3S技術[8-9]角度展開。國內方面主要運用景觀指數分析法、構建生態環境評價體系等[10-12]方法對黃河三角洲景觀格局的演變趨勢及脆弱性展開研究，對探析景觀格局演變規律，開展生態環境脆弱區的保護研究做出了貢獻。景觀格局分析主要以柵格數據作為數據源進行研究[13]。因此，選擇適宜的粒度和尺度對景觀格局的分析十分重要。Benson等[14]、秦嶺等[15]、王朋沖等[16]通過不同粒度大小的柵格數據比較得出研究區景觀格局研究的最適宜分析粒度。王璐等[17]、張金茜等[18]通過變異系數、最佳研究尺度模型的構建得出研究區最適宜分析尺度。但現有研究大都從粒度或尺度出發，對于粒度尺度的綜合分析研究較少。

黃河三角洲地理位置特殊，近年來受人類活動、黃河水沙變化、海水侵蝕等影響，生態環境較為脆弱，景觀類型變化復雜、劇烈，景觀演變速度逐漸加快[19]。因此，開展景觀格局及其脆弱性的研究能一定程度掌握區域景觀格局變化機制及脆弱性的區域特征，為區域可持續發展和生態修復提供科學依據。

鑒于此，本研究以Landsat遙感影像為數據源，以監督分類得到的景觀類型數據為基礎，通過計算景觀格局指數，運用變異系數和地統計方法，在確定最佳空間分析粒度和尺度基礎上，探究黃河三角洲地區景觀格局脆弱性時空分異特征及其空間關聯性，以期可為生態脆弱區景觀格局的研究提供一定的科學依據。

1 研究區概況

黃河三角洲是中國三大河口三角洲之一，位于萊州灣西岸和渤海灣南岸，是黃河攜帶大量泥沙淤積所形成的沖積平原，是中國最年輕的土地。本文選擇黃河三角洲所處的東營市東營區、墾利區、河口區和利津縣為研究區(37°16′—38°0′N，118°06′—119°18′E)，總面積為5 962.64 km2。黃河三角洲屬暖溫帶大陸性季風氣候，四季分明，雨熱同期，平均年降水量530～630 mm，多集中于夏季[20]。土壤類型受人為和自然因素影響存在較大差異，以潮土、鹽土為主。黃河入海口處還成立了國家級自然保護區，該區具有中國最完整、最廣闊的河口新生濕地，在全球濕地生態系統和珍稀物種保護方面都有不可替代的作用。

2 數據來源與研究方法

2.1 數據來源

本研究所使用的2005，2012，2018年黃河三角洲Landsat遙感影像數據均下載自地理空間數據云網站(http:∥www.gscloud.cn)。為保證分類的準確性，選擇5月成像質量好、地物層次豐富、無云的影像。后在ENVI 5.3中對三期影像進行輻射定標、大氣校正、幾何精矯正，裁剪等預處理，參考黃河三角洲的區域特征和已有分類成果，采用最大似然法對預處理后的影像進行監督分類，將研究區劃分為耕地、林地、草地、未利用地(包括鹽堿地)、水域、建設用地、灘涂、鹽田和養殖池8種景觀類型。對不確定的圖斑進行實地踏勘并修正，3期影像分類結果的kappa系數均在85%以上，符合精度要求，表明分類結果能夠為本次研究提供數據支持。

2.2 研究方法

2.2.1 景觀格局變化分析

(1) 構建景觀類型轉移矩陣。景觀類型轉移矩陣可以在時間上和空間上展示出景觀類型的轉移方向和變化幅度。將得到的景觀類型分布圖在ENVI中運用change detection statistics分別獲取2005—2012年，2012—2018年景觀類型轉移矩陣，利用轉移矩陣對黃河三角洲景觀類型變化過程進行定量分析。

(2) 計算景觀類型動態度。動態度常用來表示特定時間內某種土地利用類型的數量變化情況，可以定量描述其變化幅度與變化速度[21]。動態度計算公式如下：

(1)

式中：i為第i類景觀類型;t1,t2為研究時間點,Ki為t1到t2時段內第i類景觀類型的動態度;Sit1,Sit2為t1,t2時間第i類景觀類型的面積。

2.2.2 景觀格局最佳分析粒度的確定

(1) 變異系數法選取敏感性指數。變異系數是衡量序列觀測值離散程度的相對統計量[22]。本研究運用變異系數法判定各景觀指數對空間粒度變化的敏感程度，隨著景觀空間粒度的增加，若某一景觀指數的變異系數隨之增加，則說明該景觀指數對空間粒度的敏感性較高，且變異系數越大，敏感性越高，反之亦然。本研究首先在ArcGIS中將研究區各年份景觀類型矢量數據柵格化，每個年份均得到28幅不同空間粒度大小(30～300 m)的景觀類型數據。然后從景觀水平和類型水平兩個層面，從面積—邊緣指數、形狀指數、多樣性指數、聚集度指數4個方面選取具有代表性的景觀格局指數進行變異系數的計算并選出敏感指數。其中景觀水平選取了景觀面積(TA)、斑塊個數(NP)、斑塊密度(PD)、最大斑塊指數(LPI)、邊緣密度(ED)、斑塊相似度指數(LSI)、平均斑塊面積(SHAPE-MN)、面積周長分維數(PAFRAC)、蔓延度指數(CONTAG)、相似臨近比(PLADJ)、散布與并列指數(IJI)、斑塊凝聚度(COHESION)、景觀分裂指數(DIVISION)、分離度指數(SPLIT)、香農多樣性指數(SHDI)、香農均勻度指數(SHEI)、聚合度指數(AI)17個景觀指數；類型水平上選取了斑塊類型面積(CA)、斑塊個數(NP)、斑塊密度(PD)、最大斑塊指數(LPI)、邊緣密度(ED)、斑塊相似度指數(LSI)、平均斑塊面積(SHAPE-MN)、面積周長分維數(PAFRAC)、相似臨近比(PLANDJ)、散布與并列指數(IJI)、斑塊凝聚度(COHESION)、景觀分裂指數(DIVISION)、分離度指數(SPLIT)、聚合度指數(AI)14個景觀指數。

(2) 拐點識別法。拐點識別法是通過分析各景觀指數隨粒度變化的趨勢圖來確定景觀指數的適宜粒度的一種方法[23]。通過綜合分析敏感景觀指數隨粒度變化的趨勢圖確定適宜粒度區間，根據趨勢拐點來確定適宜的景觀粒度。

2.2.3 景觀格局脆弱性研究

(1) 半變異函數。半變異函數是指區域化變量考慮系統屬性在所有分離距離上任意兩樣本間的差異[24]。通常認為區域化變量的空間異質性是由空間自相關和隨機誤差所引起的。本研究分別計算3個研究年份4個不同網格大小下的景觀適應性指數LAI，利用GS+軟件分別計算其半變異函數，根據半變異函數擬合模型參數確定黃河三角洲地區景觀脆弱性分析的最適宜尺度。

(2)

式中：h為步長;N(h)是分隔距離為h時的樣本總數,當h=0時,γ(h)=C0;C0稱為塊金值,表示隨機因素產生的空間異質性;A0為變程用于表示空間相關性的作用范圍,當h增大到A0時,γ(h)從非零值達到一個相對穩定的常數;C是結構方差,表示因空間自相關而產生的空間異質性;C0+C為基臺值,表示自相關部分的空間異質性,其值大小與空間差異性成正比。

(2) 景觀脆弱度計算。本研究通過構建景觀敏感度指數LSI和景觀適應性指數LAI得到景觀脆弱度指數LVI[1]，其計算公式見(3)—(5)。利用Fragstats和ArcGIS軟件計算處理后得到不同研究年份的黃河三角洲景觀脆弱度分布圖。

(3)

LAI=PRD×SHDI×SHEI

(4)

LVI=LSI×(1-LAI)

(5)

式中：n為景觀類型數目;i為景觀類型;Ui為景觀干擾度指數;Vi為景觀類型易損度; PRD，SHDI，SHEI分別為斑塊豐富度指數、香農多樣性指數、香農豐富度指數。

(3) 空間自相關分析。空間自相關分析用來描述事物在空間上的依賴關系[25]。本研究借助GeoDa軟件，來進行景觀空間自相關分析。首先構建空間權重矩陣，再進行全局自相關和局部自相關分析。通過Moran’sI指數和局部自相關LISA顯著圖[26]，探究不同研究年份研究區各網格單元景觀脆弱度指數之間的空間相關性和某一網格單元與鄰近網格單元景觀脆弱度指數的空間關聯程度。

3 結果與分析

3.1 景觀格局動態變化

3.1.1 景觀類型轉移特征 黃河三角洲2005—2012年和2012—2018年的景觀類型轉移矩陣見表1—2。2005—2012年，未利用地和灘涂的轉出面積，分別是555.90 km2，269.04 km2。占到了各自總面積的67.84%和57.56%，主要轉移方向是耕地、建設用地、鹽田和養殖池。耕地轉出面積為416.93 km2，占總變化面積的25.61%，主要轉移方向是建設用地。建設用地面積增加最多，7 a間增加321.32 km2，其次是鹽田和養殖池，增加了310.05 km2。2012—2018年，依舊是未利用地和灘涂的轉出面積最大，分別是236.13 km2，189.166 km2，鹽田和養殖池、耕地、建設用地分別占未利用地轉出的31.59%，40.05%，16.73%，而灘涂的主要轉移方向是鹽田和養殖池，占其總轉出面積的59.74%。面積增加最多的仍然是建設用地，增加了283.45 km2，占建設用地總面積的26.71%。鹽田和養殖池增加的面積僅次于建設用地為252.32 km2。綜合來看，建設用地的面積增加了2.33倍，這與近年來黃河三角洲城市化進程加快以及經濟發展的需求有關。2005—2018年黃河三角洲未利用地和灘涂的轉出幅度最大，兩者面積分別較2005年減少20.08%和41.51%，土地利用率提高。林地和草地的總面積基本保持不變，但轉入轉出較為頻繁。水域面積基本穩定。

表1 2005-2012年黃河三角洲景觀類型轉移矩陣 km2

表2 2012-2018年黃河三角洲景觀類型轉移矩陣 km2

3.1.2 景觀動態度分析 黃河三角洲景觀類型在研究時段內呈現出不同的變化特點，景觀類型動態度見圖1。由圖1可知建設用地、鹽田和養殖池面積大幅增加。建設用地在兩個研究時段內的動態度分別為10.06%，10.01%，呈現持續增長態勢，且增長速度最快。2018年鹽田和養殖池面積是2005年的2倍，增長速度先快后慢，其中灘涂和未利用地持續減少，減少速度加快，在2012—2018年，動態度分別達到了-6.84%和-9.40%。這表明黃河三角洲地區經濟的快速發展使得城鎮化進程持續推進，土地利用效率提高。林、草地的總面積基本保持不變，草地面積的增加基本來自于林地的轉化，原因是研究區土地鹽漬化比較嚴重樹木相對難以存活往往轉化成荒草地或者沼澤。

圖1 2005-2018年黃河三角洲單一景觀類型動態度

耕地面積先增加后減少，動態度為0.74%和-1.4%，這與一開始的擴大耕地規模以滿足生產需求為主要發展模式轉變為人與自然和諧發展的政策有關。水域動態度較小，為-1.70%和0.24%，研究時段初期水體轉出方向主要為耕地，但到后期耕地轉換為水體基本達到了一個動態的平衡。

3.2 景觀格局最佳空間分析粒度

3.2.1 景觀水平上最佳空間分析粒度確定 本研究通過變異系數法判斷各景觀指數對空間粒度的敏感程度，根據景觀指數的變異系數大小，將敏感度分為5個等級，分別為極高敏感度(變異系數>10%)，高敏感度(變異系數7%～10%)，中敏感度(變異系數4%～7%)，低敏感度(1%～4%)，不敏感(變異系數<1%)。變異系數計算結果見表3。

由表3可知，NP，PD，LPI，ED，LSI，DIVISION，SPLIT這7個指數是景觀水平上的極高敏感指數。且3個年份中，2012年各個景觀指數的變異系數整體較高，因此選用2012年的數據進行研究區最佳空間分析粒度的趨勢分析。

表3 2005-2018年黃河三角洲景觀水平景觀指數隨粒度變化變異系數 %

2012年7個敏感指數隨粒度的變化趨勢見圖2。可以看出，隨著空間粒度的增加，景觀指數的變化主要分為3類：第一類特征是隨著空間粒度的增加，景觀指數先快速降低然后緩慢降低，包括NP，PD；第2類是隨著空間粒度的增加基本保持同一斜率緩慢下降，包括ED，LSI；第3類是先基本保持平穩然后隨著空間粒度的增加上下波動，達到到某一空間粒度時又恢復到穩定狀態，包括LPI，DIVISION，SPLIT。

圖2 2012年黃河三角州景觀水平上敏感景觀指數隨粒度的變化趨勢

由圖2可知，NP，PD無明顯的“拐點”，但達到200 m粒度后下降的速率明顯降低；ED，LSI也沒有明顯的拐點；DIVISION，SPLIT，LPI在110～220 m粒度區間內波動明顯。因此，從景觀水平上綜合分析景觀指數的拐點認為，黃河三角洲整體的適宜粒度區間為220～300 m。

3.2.2 類型水平上最佳空間分析粒度確定 為了保證研究的準確性，類型水平上最佳分析粒度的確定同樣重要。與景觀水平的研究思路一致，以2012年為例，選取了CA，NP，PD，LPI，ED，LSI，SHAPE-MN，PAFRAC，PLANDJ，IJI，SPLIT，COHESION，DIVISION，AI這14個景觀指數進行不同空間粒度大小下變異系數的計算。其變異系數計算結果見表4。

由表4可知，類型水平上不同景觀類型的景觀格局指數表現出不同的敏感度，所有景觀類型均表現出極高敏感的指數有NP，PD，ED，LSI。同樣，將各景觀類型變異系數與空間粒度進行趨勢分析，發現鹽田和養殖池、水體的適宜粒度區間為190～300 m；灘涂、建設用地、未利用地、林地和草地的適宜粒度區間為200～300 m；耕地的適宜粒度區間為220～300 m。

表4 2005-2018年黃河三角洲不同類型水平景觀指數隨粒度變化變異系數 %

綜合景觀水平和類型水平的分析結果，黃河三角洲地區景觀分析的最佳粒度范圍為220～300 m，最終本研究選取220 m作為黃河三角洲最佳景觀分析粒度。

3.3 黃河三角洲景觀格局脆弱性時空演變規律

3.3.1 最優尺度的確定 景觀格局脆弱性具有尺度效應，已有研究通常將網格劃分尺度設置為是平均斑塊面積的2～5倍[27]，因此本研究將研究區劃分為1.5 km×1.5 km，2.0 km×2.0 km，2.5 km×2.5 km，3.0 km×3.0 km 4種尺度的單元網格，通過分析不同網格尺度下景觀適應性指數(LAI)的半變異函數參數特征，篩選景觀格局脆弱性的最優分析尺度。3個年份不同網格尺度下景觀適應性指數的半變異函擬合參數見表5。

表5 2005-2018年黃河三角洲不同網格尺度下景觀適應性指數的半變異函數擬合模型

整體來看，3個研究年份不同尺度下的最佳模型均為指數模型。2.5 km和3.0 km尺度下的塊金值較小，說明隨機部分造成的空間異質性較少。由表5可知，2.5 km，3.0 km尺度下的空間變異程度較為合理，既能反映一定的差異性，又不會使差異性過于顯著。C是結構方差，表示因空間自相關而產生的空間異質性。C/(C0+C)是結構方差與基臺值的比例，C/(C0+C)越大，表示因自相關造成的空間異質性越多。由表5可知，各個尺度下的C/(C0+C)均大于75%，說明變量具有顯著的空間相關性。A0表示變程，只有當網格大小小于變程時，才能保證研究區景觀格局脆弱性的高度相關關系。由此可知2.5 km和3.0 km尺度更適宜探討景觀格局的空間異質性及關聯性。相比而言，3 km尺度下函數復相關系數R2更高，綜合殘差RSS更小，因此3 km×3 km的網格是黃河三角洲景觀格局脆弱性分析的最適宜尺度。

3.3.2 景觀脆弱性時空特征分析 圖3展示了基于3 km×3 km尺度的黃河三角洲景觀格局脆弱度空間分布圖，各脆弱度等級所占比例見表6。由圖3可知，在研究時段內研究區景觀脆弱度變化幅度不大，高度脆弱區大都分布在北部和東部沿海灘涂處，由于鹽田和養殖池連片式分布，且景觀類型以灘涂為主，使得景觀優勢度較低，從而導致景觀敏感性較高。較高脆弱度分布在未利用地和建設用地上，其中以未利用地為主的較高脆弱區逐漸減少，而以城鎮為主的地區呈逐年擴大趨勢。中等脆弱度地區分布較為均勻，較低脆弱度和低脆弱度地區主要分布在研究區的中部、西部和南部。

圖3 2005-2018年黃河三角洲景觀格局脆弱度及其指數(LVI)分布特征

由表6可知，低脆弱度地區呈逐年減少趨勢，研究時段內比例從14.38%下降到6.99%。這是由于為滿足農業和城市化的發展，土地利用開發程度增加，景觀破碎化增加。較低脆弱區、中等脆弱區呈增加趨勢，所占比例分別增加了2.0%和3.0%。此外，2005—2018年黃河三角洲景觀脆弱度呈現不同的空間變化特征，研究區北部的景觀格局脆弱度有變好的趨勢，但南部有惡化的趨勢，這是因為北部未利用地土地利用結構不斷優化，景觀類型逐漸豐富，這對景觀格局脆弱度有一定的緩和作用；而研究區中部和南部的墾利縣、東營區和利津縣，近年來建設用地逐漸擴張，破壞了原有的板塊連通性，因此導致景觀脆弱程度的提高。

表6 2005-2018年黃河三角洲不同等級景觀格局脆弱度網格所占比例 %

3.3.3 景觀格局脆弱性的空間關聯分析 2005，2012，2018年黃河三角洲的景觀脆弱度指數全局Moran’sI結果如表7所示，3個年份的全局Moran’sI系數分別為0.354，0.365，0.399，且顯著性水平均達到0.001，說明黃河三角洲景觀脆弱度呈現顯著的空間正相關關系，且相關關系持續增強，空間集聚現象日趨明顯。

表7 2005-2018年黃河三角州景觀格局脆弱度全局Moran’s I指數

選擇局部空間自相關系數LISA作為度量指標，探究黃河三角洲不同方向上景觀格局脆弱度的空間集聚模式(圖4)。由圖4可知，研究區以高—高集聚區、低—低集聚區為主。景觀格局脆弱度高—高區主要集中在研究區北部和東部，這是由于北部和東部對應區域的中心單元均為高、較高脆弱區，景觀類型主要以灘涂、鹽田和養殖池為主，相鄰單元正相關性強，對周邊單元的輻射作用較強。高—高區分布范圍先縮小后增加，表明景觀格局脆弱性先有所緩解但又略微抬升。低—低脆弱區分布在中部和西南部，這與此區域景觀類型分布較為均勻有關。同時低—低脆弱區的減少，主要是因為經濟的快速發展、人為擾動增強，導致景觀破碎度增加。

圖4 2005-2018年黃河三角洲地區景觀脆弱度LISA集聚特征

由圖5 LISA顯著性水平圖分析可知，景觀格局脆弱度高—高對應的顯著性較強，P大都在0.001，0.01水平上，中部和西南部的低—低區顯著性水平大都為0.5。

圖5 2005-2018年黃河三角洲景觀脆弱度局部自相關LISA顯著特征

4 討論與結論

4.1 討 論

本研究以遙感影像為數據源，綜合變異系數法、網格法、地統計等方法在確定最佳分析粒度的基礎上選取最佳分析尺度，對黃河三角洲景觀格局脆弱性展開研究，并進行空間自相關分析，從而揭示黃河三角洲景觀格局在時間上和空間上的演變規律。目前關于黃河三角洲景觀格局脆弱性的研究較少，本文能在一定程度上提高黃河三角洲景觀格局研究的準確性，為景觀格局脆弱性研究提供了新的思路。

數據源、空間粒度、尺度的不同都會對景觀格局的研究結果產生影響。因此，本研究選取了2005，2012，2018年相同月份的Landsat遙感影像，避免了數據源不同引起的差異。合理確定研究區最佳分析粒度和尺度是保證景觀格局脆弱性結果準確的前提條件。相較于其他文章僅從最佳粒度或最佳尺度一個方面展開研究，本研究綜合確定了黃河三角洲分析的最佳粒度和尺度，提高了研究的準確性。很多學者采用拐點法來選擇最佳分析粒度[18]，也不乏有學者以信息損失最小為原則，運用信息損失評價法確定最佳分析粒度[28]。然而在最佳分析尺度的選取方面，很多學者通過綜合分析研究區大小以及工作量的大小直接選取。同時，最佳粒度和尺度有很強的地域性，鑒于此，如何選取最佳空間粒度和空間幅度的方法仍需進一步研究。

影響景觀格局脆弱性的因素有很多，例如氣候條件、土壤類型、地形地貌、人類活動等。本研究受限于數據的可獲得性，僅以遙感影像為數據源來探討景觀格局的脆弱性，缺乏景觀脆弱性的影響因子方面的研究，研究結果尚不夠全面。在后續的工作中，將從氣候、地形、土壤類型、人類活動等多角度展開綜合研究，從而更加全面的揭示研究區景觀生態的演化特點和規律。

4.2 結 論

(1) 2005—2018年，黃河三角洲景觀類型發生了較大幅度的改變，建設用地面積增加2.33倍，未利用地和灘涂面積大幅減少，分別減少了各自總面積的20.08%和41.51%。

(2) 景觀水平和類型水平上，高敏感性的景觀指數在粒度達到220 m后基本相對穩定，確定220 m是研究區最適宜的空間分析粒度。根據不同網格尺度下的景觀適應性指數的半變異函數參數對比，確定3 km×3 km尺度是研究區最適宜的景觀脆弱度分析尺度。

(3) 2005—2018年，黃河三角洲地區景觀格局脆弱度逐漸增強，脆弱地區大都分布在沿海地區，主要是由于沿海地區景觀系統主要為灘涂、未利用地、鹽田和養殖池，景觀集中連片分布導致景觀適應性較低。內陸地區景觀脆弱度也有所增強。

(4) 2005，2012，2018年景觀格局全局Moran’sI逐年上升，說明3個時期景觀格局脆弱度存在正相關現象。在總體空間分布上，2005—2018年黃河三角洲景觀格局脆弱度均以高—高脆弱區與低—低脆弱區的聚集為主，且顯著性水平相對較高。