融合時序遙感分析的國土空間生態保護修復關鍵區識別

陳元鵬 周 旭 陳 妍 劉巖濤 蘇香燕 張成鵬

(1.自然資源部國土整治中心, 北京 100035; 2.中國地質大學(北京)土地科學技術學院, 北京 100083;3.南京大學地理與海洋科學學院, 南京 210023)

0 引言

國土空間生態保護修復是指遵循生態系統演替規律和內在機理,對生態功能退化、生態系統受損、空間格局失衡、自然資源開發利用不合理的生態、農業、城鎮國土空間,統籌和科學開展山水林田湖草沙一體化保護修復的活動[1-3]。國土空間生態保護修復是落實國家生態文明建設的重要措施,也是構建國家生態安全格局和統籌山水林田湖草沙系統治理的重要舉措[4-5]。作為提升生態系統結構和功能完整性的有效方案,國土空間生態保護修復已上升為國家戰略工程[6]。

國土空間生態保護修復關鍵區識別至關重要,是國土空間生態保護修復規劃編制、生態保護修復工程規劃和時空布局、生態保護修復工程措施選取等系列工作的前置條件和基礎保障,是國土空間生態保護修復工作的重點和難點[7-9]。

目前國土空間保護修復關鍵區的識別主要依賴于最小累積阻力模型[10-12]或電路理論[3, 13-14]構建的生態安全格局。“識別生態源地-構建阻力面-提取生態廊道”的框架模式成為當前生態安全格局構建的基本范式[15],其中生態源地的提取是重要基礎[16]。生態源地是物種或生態流向外擴散的起點,可以促進生態過程、維持系統完整性和提供高質量的生態系統服務,是保障區域生態安全的關鍵區域[17]。識別生態源地的方法可以分為直接法[18]與綜合評價法[19]。現在研究和采用較多的方法,是基于多指標生態系統服務并考慮人為因素識別生態源地。然而,目前的方案雖然充分考慮了變量間的數值關系,但對變量間的生態關系考慮不足,較少考慮不同生態系統服務之間的權衡協同場景[20-22]。

在目前生態安全格局構建的范式下,國土空間生態保護修復關鍵區識別的方法通常從單一時間節點出發,在識別生態源地的基礎上提取生態廊道與生態戰略點進而實現目標,但缺少對生態系統時間序列變化的考慮[23]。因此,本文在生態系統健康評估中考慮生態系統服務權衡協同的場景,更加科學合理地提取生態源地;另一方面在生態源地提取的基礎上,增加對長時間序列植被指數變化趨勢的分析,融合國土空間生態安全格局的動態和靜態因素,提出一種更全面、準確、合理的生態保護修復關鍵區域的識別技術方法,為國土空間生態保護修復本底調查、問題識別、空間規劃以及工程布局等提供技術支撐。

1 研究區與數據源

1.1 研究區概況

研究區為河南省某生態保護修復工程所在地,位于秦嶺山脈東段北部、河南省西部,地理坐標為33°38′33″～35°04′50″N、110°21′50″～113°36′57″E。該區域地形地貌變化劇烈,由西向東形成秦嶺中低山地、黃土丘陵臺塬、伊洛山間盆地、黃河灘區濕地的過渡帶地貌,第1級地貌為西部小秦嶺、伏牛山等中山區,平均海拔1 000～2 000 m,部分山峰海拔超過2 000 m;第2級地貌為中部崤山、嵩山等低山丘陵一帶,平均海拔200～1 000 m;第3級地貌為東北部伊洛盆地和黃河灘區,平均海拔200 m以下。研究區地形地貌如圖1所示。

圖1 研究區地形地貌圖Fig.1 Topographic of study area

特有的區域氣候和地貌特征奠定了實施區森林、濕地、河湖、農田、城鎮等5類陸地生態系統發育與演變的自然基礎,以及社會經濟發展的空間格局。

研究區域主要為暖溫帶季風型氣候,四季分明,年平均氣溫14.7℃,1月最冷,月平均氣溫0.3℃,7月最熱,月平均氣溫27.5℃。年平均降水量606.9 mm,降水主要集中在6—9月,占全年降水量的62.4%。多年平均蒸發量為1 829.10 mm。

1.2 數據來源及處理

本文所用數據包括氣象數據、遙感數據、土地覆被數據等。氣象數據來源于中國科學院資源環境科學與數據中心網站,中國氣象要素年度空間插值數據集,包括2020年氣溫、降水量、蒸發量等,空間分辨率為1 km。土壤數據來源于中國科學院資源環境科學與數據中心網站1∶1 000 000土壤數據庫。地形數據為地面數字高程模型數據(DEM),來源于地理空間數據云SRTMDEM,空間分辨率為90 m。遙感數據為基于MOD13A1構建的時序(2011—2020年)歸一化植被指數(NDVI)數據集,共計230景,空間分辨率為500 m。土地覆被數據為MODIS數據產品MCD12Q1土地覆被類型(Land cover),選擇其中2010、2015、2020年3期,空間分辨率為 500 m。 植被凈初級生產力(NPP)來源于MODIS數據產品MOD17A3HGF數據集,時間序列為2011—2020年,空間分辨率為500 m。其中,時序NDVI、NPP、土地覆被類型數據均通過Google Earth Engine(GEE)云平臺獲取并批量處理計算。

GEE是當前世界上較為先進的處理衛星影像等地理空間觀測數據的云計算平臺。GEE云端數據庫中集成了近40年的Landsat系列衛星的歷史存檔數據,給個人用戶提供了強大的算力和云存儲空間,同時提供了方便快捷的JavaScript語言API接口進行數據處理、算法實現以及結果分析[24-25]。本研究中,通過GEE云計算平臺對研究所需的遙感、土地覆被、NPP數據進行了批量處理,減少了前期大量數據的準備工作,有效提升了工作效率,降低了數據處理與算法實現過程中對本地硬件設備的依賴。

2 研究方法

本文研究方法框架主要包括3部分:首先,基于GEE云平臺、MOD13A1數據,采用Mann-Kendall方法對研究區2011—2020年時序NDVI開展趨勢分析,形成時序分析結果;其次,基于氣象、地形、土壤、NPP等數據,開展基于生態系統服務協同和權衡計算的生態系統健康分析評估,選擇生態源地;最后,基于時序分析結果和生態源地選擇結果,采用疊置分析,開展生態保護修復關鍵區域識別。具體技術路線如圖2所示。

圖2 技術流程圖Fig.2 Workflow of this study

2.1 Mann-Kendall趨勢分析(時序遙感趨勢分析)

Theil-Sen Median(TSM)趨勢分析[26]和Mann-Kendall(MK)非參數檢驗法[27]常用來分析NDVI空間演變特征。本文采用 Theil-Sen Median 趨勢分析方法并結合MK檢驗方法判斷Sen趨勢的顯著性,與一元線性回歸趨勢分析方法相比,其可以規避時間序列數據分布和數據缺失對分析結果的影響[28],同時可以降低異常值對結果影響[29]。TSM計算公式為

(1)

式中i、j——時間序列數

NDVIi、NDVIj——第i、j時間序列的NDVI值

Median——多年NDVI中位數函數

n——時間序列長度

MK是一種非參數統計檢驗方法,用來判斷趨勢的顯著性,其計算公式為

(2)

(3)

本文使用顯著性水平α=0.05,對應的u1-α/2為1.96[30]。

2.2 基于生態系統健康的生態源地識別

本文在生態系統服務權衡與協同框架下開展生態系統健康評估,進而開展生態源地識別。

2.2.1生態系統健康評估體系

生態源地是對區域生態過程和功能起決定性作用的生境斑塊,提供必要的生態系統服務,對區域生態系統健康安全具有重要意義。生態系統健康(Ecosystem health,EH)評估是生態源地識別的基礎[31],本文從生態系統活力(Ecosystem vigor)、生態系統結構(Ecosystem organization)、生態系統彈性(Ecosystem resilience)、生態系統服務(Ecosystem service,ES)4個維度構建評估框架體系[32]。計算公式分別為

(4)

(5)

式中EH——生態系統健康水平

PH——生態系統物理健康水平

ES——生態系統服務

V——生態系統活力

O——生態系統結構

R——生態系統彈性

ES、V、O和R計算結果均進行歸一化處理。

V可用NPP表征,因為NPP能夠很好地反映植物固定和轉化光合產物的效率,并表征可供人類使用的物質和能量。

O表示區域生態系統結構的穩定性,主要包括景觀異質性和連通性。為進一步描述景觀異質性和連通性,選擇景觀格局指數(包括面積加權平均斑塊分維數(AWMPFD)、香農多樣性指數(SHDI)、景觀聚集度指數(CONT)、景觀破碎化指數(FN))用以表征O,公式為

O=0.3SHDI+0.2AWMPFD+
0.25FN+0.25CONT

(6)

式(6)中,景觀格局指數采用Fragstats軟件計算。

R表示生態系統在面對人為擾動時保持自身結構穩定性的能力,反映一個區域在生態系統過程中抵御和適應外部干擾的能力,因此包括兩方面,即抵抗力和恢復力[33],本研究中為抵抗力和恢復力分別賦值0.6和0.4[34],為反映同一土地覆被類型不同生態系統的彈性,采用NDVI進行校正,公式為

(7)

式中Ri——第i個像元總彈性系數

NDVIm——第m個像元的NDVI值

NDVImeanj——第j類土地覆被類型的NDVI均值

Resilj、Resisj——第j類土地覆被類型的恢復力系數和抵抗力系數(表1)

表1 生態系統彈性系數Tab.1 Ecological resilience coefficient of land use types

生態系統服務是自然系統和社會經濟系統之間的橋梁[35],反映了人與自然之間的耦合關系,是生態系統健康的重要組成部分。針對研究區的自然特征和生態環境狀況,本文選取水源涵養、水土保持、防風固沙3類生態系統服務作為區域關鍵生態系統服務功能進行定量評估。3類生態系統服務按照凈初級生產力(NPP)定量指標評估法計算,計算公式分別為

WR=NPPmeanFsicFpre(1-Fslo)

(8)

Spro=NPPmean(1-K)(1-Fslo)

(9)

Sws=NPPmeanKFqD

(10)

其中

式中WR——生態系統水源涵養服務能力指數

NPPmean——多年植被凈初級生產力平均值

Fsic——土壤滲流因子

Fpre——多年平均降水量因子

Fslo——坡度因子

Spro——生態系統水土保持服務能力指數

K——土壤可蝕性因子

Sws——生態系統防風固沙服務能力指數

Fq——多年平均氣候侵蝕力

u——高度2 m的月平均風速

ETPi——月潛在蒸發量

Pi——月降水量

d——當月天數

Ti——月平均氣溫

ri——月平均相對濕度

D——地表粗糙度因子

θ——坡度

2.2.2生態系統服務權衡與協同

在生態系統服務的背景下,生態源地的選擇過程必須充分考慮生態系統服務的權衡和協同關系。目前,有序加權平均模型(Ordered weighted average model,OWA)已有效用于量化生態系統服務之間的權衡和協同關系[36],因此基于OWA可量化不同場景下的生態系統服務供給,反映不同生態系統服務之間的權衡和協同效應[34, 37],其公式為

(11)

式中axj——第x柵格數據中第j個像元的屬性值

Sxj——重新賦權后的柵格數據

ωx——Sxj的權重

N——參與計算的生態系統服務柵格數據數量

此外,采用風險值(risk)和權衡值(trade-off)兩個參數作為OWA的約束和優化條件,公式為

(12)

(13)

(14)

通過計算風險值和權衡值,可以得出不同場景模式下的最優化的生態系統服務權重系數組合。計算最大權衡值和相應的權重,共設定11個場景。

2.2.3生態源地識別

為有效集約地選擇生態源地,對于不同的場景模式,選擇生態系統健康值前60%的區域作為優先保護區域,并采用生態系統保護效率公式在優先保護區域中進一步篩選生態源地。公式為[38]

(15)

2.3 生態保護修復關鍵區識別

以時序遙感趨勢分析中NDVI的下降趨勢、保持穩定和輕微上升趨勢,以及生態系統健康評估和生態源地識別結果為基礎,采用疊置分析,進一步識別生態保護修復關鍵區[39]。

3 結果與分析

3.1 時序遙感趨勢

基于GEE云平臺得到了研究區域2011—2021年的NDVI時間序列數據,共計230景,從統計描述看,230景NDVI數據的均值在0.16～0.76范圍內波動,波動周期符合年度季相特征。標準差在0.06～0.23范圍內波動,穩定性較好,數據質量較好,時序NDVI數據構建結果符合預期。研究區數據統計見圖3。

圖3 時序NDVI數據的均值和方差Fig.3 Mean and variance of time-series NDVI data

基于TSM趨勢分析,結合Mann-Kendall非參數檢驗法,在顯著性水平α=0.05下,形成時序NDVI趨勢變化統計數據和空間分布(圖4)。首先,從NDVI趨勢變化統計數據結果看,研究區域內時序NDVI的斜率變化范圍為-0.001 7～0.002 1,均值為0.000 3,其中斜率為負值的像元6 136個,占比5.33%,說明2011—2020年間,研究區域內少數局部區域NDVI呈下降趨勢,但整體NDVI變化呈上升趨勢,一定程度反映出研究區域內的生態系統質量穩中向好。其次,疊加土地覆被類型(圖5)看NDVI趨勢變化空間分布可知,2011—2020年間研究區域內NDVI呈下降趨勢的區域主要分布在中部洛河、伊河沿線的不透水面集中地區及周邊緩沖區域,東北側黃河沿線的不透水面集中地區及周邊緩沖區域,以及西北角小秦嶺北側的山前平原,其余區域還有些零星分布。NDVI呈上升趨勢主要分布在中部黃土丘陵和東部山間盆地區域,主要土地覆被類型為耕地、草地和少量灌木林地等。NDVI變化趨勢不顯著的區域主要分布在秦嶺中低山地區域,包括崤山、熊耳山、伏牛山等山脈,其土地覆被類型以林地為主。

圖5 2020年土地覆被類型Fig.5 Spatial distribution of land cover types in 2020

3.2 生態系統健康評估結果

3.2.1生態系統物理健康(PH)

根據式(5)～(7),分別計算生態系統活力(V)、生態系統結構(O)、生態系統彈性(R),并經歸一化處理,進一步計算得到生態系統物理健康(PH)評估結果及其空間分布(圖6)。如圖6所示,生態系統物理健康均值為0.592,研究區域內高值區主要分布在西側秦嶺中低山地區域,如小秦嶺、熊耳山、伏牛山等,中低山地區域植被總體生長良好,生態系統健康程度高。由西向東,生態系統物理健康值逐漸下降,低值區域主要分布在伊洛山間盆地、河流灘區等土壤侵蝕嚴重、水土流失嚴重的區域。此外,城市化程度較高的地區,由于人類活動的干擾較大,生態系統物理健康水平較低。

圖6 生態系統物理健康(PH)評估結果和空間分布Fig.6 Evaluation result and spatial distribution of ecosystem physical health

3.2.2生態系統服務(ES)

根據式(8)～(10)凈初級生產力(NPP)定量指標評估算法,并經歸一化處理,計算得到生態系統服務(ES)評估結果及其空間分布(圖7)。不同生態系統服務能力之間存在空間異質性,反映了生態系統服務能力之間的權衡性,也揭示了不同生態過程對于研究區整體生態安全格局的影響和作用。如圖7所示,研究區域內水源涵養服務能力由西南向東北逐步減弱,水源涵養能力較強的區域主要分布在中低山區,但并非區域內所有中低山區都具有較強的水源涵養服務能力,較強區域主要集中在伏牛山-熊耳山-外方山區域。較于水源涵養能力,水土保持能力分布相對均衡,重要區域主要分布在中低山、黃土丘陵區以及陸渾水庫周邊區域。防風固沙能力與水源涵養能力存在一定協同性,較強區域主要集中在中低山區,其中最強的區域主要分布在小秦嶺、熊耳山、伏牛山區域,次強區域為崤山、青要山、外方山和嵩山區域。從單一生態系統服務空間格局看,研究區域內生態系統服務能力較強的區域主要分布在中低山區、黃土丘陵區等海拔較高且植被覆蓋較強的區域。

圖7 生態系統服務(ES)評估結果和空間分布Fig.7 Evaluation result and spatial distribution of different ecosystem services

3.2.3生態系統服務協同與權衡

根據式(11)～(14)生態系統服務權衡與協同算法,通過改變決策風險和協同權衡值,可以生成一系列場景。結合工作實際,以0為起始、0.1為間隔設置風險值,計算最大協同權衡值及其相應的權重,經計算共得到11個場景(表2)。

表2 不同風險場景及權衡水平下的最優系數Tab.2 Optimal order weights for different decision risk scenarios and trade-off levels

3.2.4生態系統健康(EH)評估和生態源地識別結果

根據生態系統物理健康(PH)和生態系統服務(ES)計算結果,結合生態系統服務權衡與協同的11個場景,按照式(4)計算得出11個生態系統健康(EH)評估結果及其空間分布特征(圖8)。

圖8 生態系統健康(EH)評估結果及其空間分布Fig.8 Evaluation result and spatial distribution of ecosystem health in different scenarios

如圖8所示,生態系統健康(EH)的高值和低值空間分布呈一致趨勢,但局部的高值區和低值區也存在一定差異,這是由于不同的場景、不同的生態系統服務權衡與協同水平下,權重等次不同,一定程度上說明了簡單的加權求和不能充分反映區域生態系統健康(EH)的空間分布,因此,有必要將生態系統服務協同權衡效應納入生態系統健康評估。

基于11個場景的生態系統健康評估結果,按照高值的前60%篩選了不同場景下的優先保護區域,并根據式(15)測量了其保護效率(表3)。如表3所示,場景9的保護效率最高,為3.02,因此選擇場景9的高值前60%區域作為生態源地,其空間分布如 圖9 所示。研究區域內生態源地約1 386.11 km2,主要分布在西南部的伏牛山-熊耳山-外方山區域,少量分布在小秦嶺區域,涉及的土地覆被類型主要為林地和草地。

表3 不同場景下的優先保護區域保護效率EiTab.3 Protection efficiency of priority protected areas in different scenarios

圖9 生態源地選擇結果及其空間分布Fig.9 Select result and spatial distribution of ecological sources

3.3 生態保護修復關鍵區

以時序遙感趨勢分析結果(圖4)為基礎,選取NDVI下降趨勢(斜率小于0)、NDVI保持穩定和輕微上升趨勢(斜率在0～0.000 2之間)的區域,與生態源地選擇結果(圖9)進行疊置分析。

將NDVI保持穩定和輕微上升的區域和生態源地合并為生態保護關鍵區,用以表征研究區域內需生態保護的關鍵區域,其中:①將NDVI保持穩定和輕微上升區與生態源地重合的區域劃定為一級生態保護區,因為該區域既是保持生態系統健康的重要區域,但從時序上看,近10年間該區域可以表征生態系統的NDVI指數并未形成顯著變化趨勢,一直屬臨界狀態,如不作為重點保護區域則可能出現生態系統退化或破壞情況,因此該區域有必要作為一級生態保護區,以逐步形成生態系統穩定性較強、NDVI指數穩中有升的區域。②將NDVI上升區與生態源地重合的區域劃定為二級生態保護區,二級生態保護區同樣是保持生態系統健康的重要區域,但與一級生態保護區相比,近10年間區域內NDVI已形成上升趨勢,生態系統已穩中向好,故此作為二級生態保護區。③將未與生態源地重合的NDVI保持穩定和輕微上升區劃定為三級生態保護區,與一、二級生態保護區相比較,三級生態保護區未包含生態源地,其對于保持區域內生態系統完整、提供高質量生態系統服務的貢獻度相對不強,但近10年間該區域內生態系統未遭到明顯破壞,因此有必要進行保護以維持現狀或將其逐步轉化成為生態系統穩定性較強、NDVI指數穩中有升的區域。④將NDVI下降區劃定為生態修復區,采取適當的措施實施生態修復。形成的生態保護修復關鍵區識別結果及其空間分布如圖10所示。

圖10 生態保護修復關鍵區識別結果及其空間分布Fig.10 Spatial distribution of ecological protection and restoration key areas identification results

由圖10可知,一、二級生態保護區主要分布在西部和西南部小秦嶺、熊耳山和伏牛山等區域,主要為中山森林生態系統,承擔了生物多樣性保護、固碳釋氧等生態系統服務功能;三級生態保護區主要分布在邙山、崤山區域,在伊洛盆地和黃河灘區也有零星分布,主要為低山丘陵森林生態系統、山間盆地農田生態系統、河流生態系統,承擔了水源涵養、水土保持、固碳釋氧、礦產資源提供、農林產品提供生態系統服務功能;生態修復區主要分布在中部洛河、伊河沿線和下游不透水面集中地區及周邊緩沖區域,東北側黃河沿線的不透水面集中地區及周邊緩沖區域,以及西北角小秦嶺北側的山前平原,主要為盆地區城鎮、農田、河流生態系統,承擔了農產品提供、城鎮生態涵養生態系統服務功能。

4 結束語

基于時序遙感分析、生態系統服務協同權衡計算、生態系統健康評估、生態源地識別和空間疊置分析等方法,對研究區的生態保護修復關鍵區進行了有效識別。結果表明,提出的“基于生態系統服務協同權衡計算的生態系統健康評估——時序遙感趨勢分析”的研究框架,對于識別構建區域生態保護修復關鍵區尤為重要。結合生態系統服務協同權衡的生態系統健康評估能夠更有效地選擇關鍵生態源地,在此基礎上創新性地結合時序遙感分析進一步識別生態保護修復關鍵區,不僅顧及了研究區域內“靜態”的生態系統服務、生態系統健康屬性,同時衡量了“動態”的生態系統變化趨勢,更全面、準確地識別了研究區的生態保護修復關鍵區域。研究成果可為國土空間生態保護修復本底調查、問題識別、規劃和工程布局等提供技術支撐。