基于改進遙感生態指數的寧夏沿黃平原區生態環境質量評價

董春媛,喬榮榮,楊智程,羅立輝,常學禮,*

1 魯東大學資源與環境工程學院, 煙臺 264025 2 南京大學生命科學學院, 南京 210008 3 中國科學院西北生態環境資源研究院, 蘭州 730000

生態環境質量是指一定時間范圍和區域內對影響其社會發展與人類活動的部分或全部生態因素進行定性分析與評判的區域性生態環境優劣程度[1]。區域生態環境質量評價是國民經濟建設與可持續性發展規劃的基礎,為此國家環境保護部在2006年以行業標準的形式頒布了《生態環境狀況評價技術規范》[2],其后在全國一些省、縣行政單元進行了生態環境評價[3—6]。在應用中由于規范中的核心指標生態指數(Ecological Index,EI)不具備空間表達且對主觀因素(權重)依賴較大,存在無法精準落實到應用規劃圖中的短板[7]。針對這一問題徐涵秋在2013年提出了一種基于遙感影像信息提取的遙感生態指數(Remote Sensing Based Ecological Index,RSEI),該指數可以很好的彌補EI指數不能進行空間表達的缺憾,并借助遙感技術數據獲取容易、計算程序化等優點迅速在諸多地區生態環境質量變化與評價中得到應用[7—10]。RSEI是以遙感數據計算的綠度(用歸一化植被指數(Normalized Difference Vegetation Index, NDVI)表示)、濕度(Wet Index,WI)、干度(用歸一化裸土指數和建筑指數(Normalized Difference Building-Soil Index,NDBSI)平均值表示)和熱度(用地表溫度(Land Surface Temperature, LST)表示)四個指標為輸入量,經過主成分分析(Principal Component Analysis,PCA)獲得分量,然后用各因子在分量軸上載荷向量值判斷其代表性和對軸的影響力[11]。與EI相比,RSEI具有客觀、易獲得和二維空間表達的優勢,對后續區域發展規劃和生態環境質量評價具有重要支撐意義[12—14]。但是,RSEI評價因子NDBSI與LST之間存在的高度正相關引起導入因子(圖層)生態學意義重復帶來的偏差沒有得到重視;同時也沒有考慮基于柵格計算的RSEI如何概括表達相鄰柵格性質差異所帶來的影響,如某一柵格只與另一類柵格相鄰或與多種類型柵格相鄰所產生的RSEI不同無法表達。那么如何修正上述RSEI計算中的漏洞?本研究擬采用景觀多樣性指數(Landscape Diversity Index,LDI)替換LST對普遍采用的RSEI指數進行改進,稱其為改進遙感生態指數(Modified Remote Sensing Ecological Index,MRSEI),主要依據包括:(1)與其它三個指標相比LST反映的是短時間尺度地表溫度狀態,而LDI因子與其他因子屬性相似都是相對較長時間尺度地表屬性累積狀態表達;(2)LDI計算考慮到了基本分析單元周邊的組成類型數和其面積比例信息量。

在MRSEI應用中,若用LDI替代LST首先需要解決的問題是如何確定適合的LDI分析尺度,即:確定LDI尺度依賴特征。已有的景觀尺度效應研究表明,在水域分布為特點的巢湖和南四湖分別采用3000 m×3000 m和1000 m×1000 m分析網格(基本分析單元)進行研究獲得了較為合理的計算結果[15—16];而在山西平朔礦山復墾生態恢復區的研究則以500 m×500 m網格獲得了較為可信的結果[17]。由此可以看出不論是在景觀相似或相異區域其尺度依賴規律都存在差異[15—18]。因此識別研究區LDI尺度依賴特征不僅是MRSEI計算中首先要確定的問題,同時也是景觀變化研究中恰當的分析尺度應用的基礎。

寧夏沿黃平原處于干旱、半干旱交錯帶屬于典型人工綠洲。近幾十年高速的非均質化的城鎮和農田擴張導致域內生態環境質量空間異質性發生了極大變化[19—20]。為了深入解析寧夏沿黃平原區生態環境質量空間分布特征,本文擬采用景觀多樣性指數LDI替代LST的RSEI計算方法對研究區生態環境質量空間分布格局進行分析。擬解決的科學問題是:(1)在干旱半干旱人工綠洲生態環境質量評價中MRSEI計算的合理性與生態學意義解釋;(2)寧夏沿黃平原生態環境質量空間異質性與總體生態環境現狀。

1 研究區概況

寧夏沿黃平原位于寧夏回族自治區北部,范圍介于北緯 37°20′—39°20′,東經 105°0′—107°0′,面積約10831.3 km2(圖1)。該區北起石嘴山,南止黃土高原,東到鄂爾多斯臺地,西接賀蘭山。寧夏沿黃平原核心區是由黃河串聯起來的衛寧灌域和銀川灌域組成的人工綠洲區域,綠洲所占比例60%以上[20]。研究區內湖泊和濕地密布,人工灌渠縱橫,主要有東干渠、西干渠、漢延渠、唐徠渠等。該區土壤以隱域性灌淤土和草甸土為主,天然植被以沿黃河分布的沙棗林和零散分布的灌叢濕地植被為主。主要植物種有沙棗(ElaeagnusangustifoliaLinn.)、枸杞(LyciumchinenseMiller)、檉柳(TamarixchinensisLour.)和蘆葦(Phragmitesaustralis(Cav.) Trin. ex Steud.)等為主。該區是國家級沿黃經濟區的核心地區,也是國家生態功能區劃中的重點區域。

圖1 研究區位置Fig.1 Location of the study area 研究區涉及到各級行政邊界來源于全國地理信息資源目錄服務系統

2 數據與方法

2.1 數據

為了獲得準確的景觀多樣性尺度依賴特征,分析過程采用研究區多年土地利用/覆蓋數據平均值進行。其中2000年、2010年和2020年數據來源于中國向全世界免費分享的GLOBELAND30數據(http://www.globallandcover.com/)。1975年和1987年由課題組下載的Landsat MSS和TM遙感數據,在遵循上述產品分類原則并參考中國科學院資源環境科學數據中心上世紀80年代土地利用/覆蓋矢量數據,通過野外和歷史訪問調查在寧夏林科院專家指導下解譯完成。研究區土地利用/覆蓋類型為農田、喬木林地、灌木林地、草地、濕地、水體、人工地表和裸地等8類。在Google Earth Engine(GEE)平臺上對2020年7月19日和28日的Landsat OLI數據進行計算(網址:https://developers.google.com/earth-engine/datasets/catalog/Landsat_LC08_C01_T1_SR),研究區涉及到的軌道號分別為:129—33、129—34和130—34。

2.2 方法

2.2.1景觀多樣性指數與閾值

在已有的研究中RSEI在計算流程和評價因子選擇上基本固定[21—22],在本文研究中為了避免存在評價因子生態學含義重復表達的問題(NDBSI和LST),采用分析確定后最佳分析尺度的LDI替換LST。本文的NDVI、WI和NDBSI計算與已有文獻一致[7],LDI和RSEI的計算如下:

景觀多樣性指數(LDI):借用植物多樣性中Shannon-Weiner index進行計算,計算過程在ArcGIS的鄰域分析中完成。

(1)

式中,Pi為景觀中第i類土地利用類型(或類型分級)占分析單元面積比例。

由于本文采用基礎土地利用數據分辨率為30 m,而且涉及到8類土地利用類型,故文中不同尺度LDI計算從3×3個柵格(90 m×90 m)起始,確保最小分析尺度也有包含8種類型的能力。同時,為了減少計算量,分析單位按照柵格倍數逐級增大,從90 m×90 m起,其后依次為300 m×300 m、600 m×600 m、900 m×900 m、1200 m×1200 m、1500 m×1500 m、3000 m×3000 m、4500 m×4500 m、6000 m×6000 m 共9個梯度。最終LDI尺度依賴性確定采用1975年、1987年、2000年、2010年和2020年共5年的土地利用/覆蓋數據計算結果平均值判定。

2.2.2評價因子歸一化與MRSEI計算

針對上述4個輸入因子在賦值單位和數值變化幅度上存在差異,需要進行標準化處理使指標值統一到0—1之間。經極差歸一化處理后獲得輸入因子圖層(圖2)。

圖2 PCA分析中的評價層Fig.2 Evaluation layer in PCAPCA: 主成分分析 Principal component analysis

(2)

式中, MRSEI為改進遙感生態指數,MRSEI值越大,生態環境越好;反之亦然。n表示主成分特征根累積達到90%以上的分量數,在MRSEI分析中因為輸入因子為4,故1≤n≤4(i為整數)。EVi和ETi分別為某個評價因子的特征值和特征向量。其它因子定義同前。全域可視化表達是建立在采用歸一化差異水體指數剔除研究區內水體信息后的基礎上。此外,文中MRSEI空間異質性分析是建立在分級基礎上,分級原則采用均值標準差法[21],將分級間距分為4級(表1中第一和第二列)。

表1 MRSEI分級

最后,整個研究區MRSEI水平(評價)可通過公式(3)計算結果與表1中分級閾值比較確定研究區整體生態環境質量級別。

(3)

式中,MRSEI為時間t的遙感生態指數,Mi為評價區域MRSEI第i類級別中值,PAi為評價區域MRSEI第i類級別相對面積。文中統計分析在SPSS和Excel中完成,回歸方程顯著性檢查采用F檢查(P<0.001為極顯著),關聯系數顯著性檢查采用R顯著性檢查。閾值分別為R4,0.05=0.811,R4,0.01=0.917和R4,0.001=0.991。需要補充說明的是本文中的MRSEI計算考慮到了PC軸數目累計貢獻率要超過90%,很顯然參與計算PC軸數大于等于1,與其它研究在方法上存在一些不同[7—10]。

3 結果

3.1 LDI的尺度依賴性

從研究采用的5個時段以土地利用/覆蓋類型為基礎計算的LDI平均值的尺度依賴特征來看,存在顯著一元二次方程變化規律(P<0.001,R=0.988)。在分析矩形邊長6000 m范圍內可以捕捉到LDI變化拐點(在3000 m處,圖3中的三角點)。當分析尺度小于拐點時,LDI增加非常陡峭,大于拐點時增加趨于平緩。因此文中LDI指數采用3000 m×3000 m的基本單元計算,計算完成后用ArcMap數據管理中的重采樣功能,用近鄰法重采樣至30 m×30 m與其它三個因子圖層在分辨率上保持一致。

3.2 MRSEI主成分分析特征

從PCA輸出評價因子與MRSEI的關聯分析結果來看(表2),在所有評價因子中NDVI與MRSEI關聯程度最大為0.7878,但仍未通過最低統計學顯著性檢查(R4,0.05=0.811),說明采用任何一個因子都無法獲得通過統計學顯著性檢查的MRSEI狀況。同時,從各評價因子與MRSEI的關聯關系來看特點非常明顯,WI和NDBSI與MRSEI呈負相關,NDVI和LDI與MRSEI呈正相關。此外,PCA結果顯示PC1軸占特征值信息的68.98%,PC2軸占28.76%,二軸合計承載了特征值97.74%。可以確定,采用此2軸依據公式(2)可獲得所有參與評價因子絕大部分信息。從各評價因子對PC1軸信息集成的貢獻大小(特征向量大小和正負關系)來看(表3),PC1主要貢獻者是NDVI,特征向量為0.8901;其它3個評價因子特征向量為負值,變化在-0.4146—-0.0317之間。PC2主要貢獻者是LDI,特征向量為0.9100,其次是NDVI特征向量為0.4056;WI和NDBSI特征向量為負且最大不超過-0.0200。

表2 評價因子與MRSEI 和RSEI的相關系數矩陣

表3 不同評價因子PCA分析特征向量與特征值

從LDI替換LST對PCA結果的影響來看(表2),在MRSEI中LDI與MRSEI正相關與RSEI中LST與RSEI關系相反,說明LDI引入避免了RSEI中LST與NDBSI對結果的影響高度一致重復表達的不足(二者關聯系數分別為-0.9617和-0.9330)。同時,從統計學顯著性檢查結果來看,在RSEI中4個評價因子關聯系數顯著性都超過了0.05水平,說明在人工綠洲區采用任一因子單獨評價區域生態環境質量也可獲得具有統計學意義的結果,這明顯偏離了PCA分析應用的主旨。而反觀MRSEI與評價因子的關聯分析結果有效的避免了這種現象發生(表2)。此外,從PCA向量特征來看RSEI分析中的PC1累計貢獻率就在91%以上,從多維投影向量合成角度來看,RSEI分析采用的評價因子在結果貢獻上具有相似作用(聚集于PC1);而在MRSEI中PC1貢獻率不足70%,引入的LDI因子在PC2中占優勢并使特征向量貢獻率達到了28%以上(表3)。這一結果從細節上輔證了LDI替代LST在人工綠洲區遙感生態環境質量評價中的合理性。

從研究區總體環境質量評價角度來看,在MRSEI評價中由于PC1和PC2軸累積貢獻率超過90%以上,故依據表3中第一、二列數據采用公式(2)在ArcMap柵格計算器中對PCA輸出的第一、二圖層分別計算獲得了圖4中的研究區2020年MRSEI分布現狀圖,其中白色區域是采用歸一化差異水體指數剔除的研究區內的水體信息。其后通過對MRSEI分布現狀圖的屬性表進行分析獲得研究區MRSEI平均值和標準差并依據表1的原則進行分類獲得寧夏沿黃平原區MRSEI生態環境質量分級圖。最后逐級進行MRSEI中值計算(表1第4列)。

圖4 寧夏沿黃平原區2020年MRSEI分布現狀與分級Fig.4 Distribution status and classification of MRSEI of Ningxia Plain along the Yellow River in 2020

3.3 MRSEI空間異質性特征

從圖4 的MRSEI分布現狀來看,MRSEI低值區主要分布在研究區周邊,即賀蘭山東麓向寧夏沿黃平原過渡區和黃河東岸鄂爾多斯臺地與黃土高原向黃河階地過渡區。從圖4 的MRSEI生態環境質量分級結果來看,“差”和“較差”級別在研究區周邊呈連續分布格局,其中西側MRSEI以“差”級別為主,東側以“較差”級別為主。從不同MRSEI級別分布格局的數量化特征來看(表4),在分布面積方面,以“較差”級別面積最大為3528.7 km2,占研究區的31.9%;其次為“較好”級別,面積為3370.7 km2,占31.9%;其后依次為級別“好”和“差”,面積分別為1855.4 km2和1804.8 km2,所占比例分別為17.6%和17.1%。在斑塊密度方面,從“差”到“好”級別梯度上,斑塊密度為減少趨勢,由“差”級別的8.3個/km2減少到5.9個/km2。說明隨MRSEI增加其分布趨向集中分布,這一點從平均斑塊面積變化特點中可以得到印證,在同一梯度方向平均斑塊面積由0.120 km2增加到0.169 km2。在同類斑塊空間距離方面,“差”和“好”級別斑塊間相隔相對較遠,分別為253.4 m和236.2 m,相互之間的連通性相對較低;“較差”和“較好”級別分別為210.3 m和206.4 m。相互之間連通性相對較高。

表4 不同MRSEI級別空間格局數量特征

最后,研究區總體MRSEI水平可依據表4中第一列數據與研究區總面積比和表1中分級計算獲取的中值采用公式(3)完成計算。研究結果表明研究區整體MRSEI為0.0117,查表1中分級閾值可知,研究區生態環境質量剛好達到較好級別(研究區MRSEI值為0.0117大于較好級別最低閾限值0.0116)。

4 討論

從2013年RSEI在中國首次提出并應用于福建省長汀縣生態環境質量評價以來[7],RSEI的計算方法被許多研究者應用來分析區域環境質量[23—28]。這些研究都延續采用NDVI、WI、NDBSI和LST四個指標,這不僅忽視了遙感手段獲取的LST是表征的短時間尺度地表屬性量[29],而其它三個指標是相對較長時間地表屬性累積的反映,也忽視了在許多區域類型中NDBSI和LST在生態學意義解釋上高度相似的性質。因為大量遙感指數反演研究表明,地表溫度與人工地表和自然裸地有極高的正關聯[30—31]。此外,在RSEI分析中把研究基本單元(柵格)設定成相互獨立狀態,忽視了柵格與同類或異類或多種異類相鄰在生態環境質量評價中的意義,因此用能夠表達相鄰柵格異質性的生態環境評價指標替代LST成為提高RSEI在生態環境評價應用中的可靠性必須面對的問題。

LDI是景觀生態學研究中使用頻率最高的一個指數,它是基本分析單元(或研究區)中組成類型多樣性的綜合表達,其大小與組成類型數和面積比例信息均勻程度相關。LDI在不同研究區都表現出隨分析單元大小而變化的尺度依賴特征[32]。從寧夏沿黃平原LDI尺度依賴特征結果來看,當分析尺度(正方形邊長)達到3000 m時LDI變化趨于平緩(圖3),在尺度分別增加1.5倍(4500 m)和2倍(6000 m)時,LDI僅分別增了0.069和0.092;在尺度減小0.5倍(1500 m)時,LDI減少了0.123。因此可以推斷,在干旱半干旱人工綠洲區進行景觀多樣性研究中最佳閾值在3000 m。此外,需要指出的是寧夏沿黃平原面積為10831.3 km2,考慮到尺度效應不僅受分析單元的影響,而且還受研究區范圍大小的影響[33],所以其它類似地區應用時要考慮到研究區規模大小,可選擇3000 m×3000 m左右若干梯度做簡約分析獲得適宜分析尺度。

從MRSEI評價過程中的PCA結果來看,用LDI替代LST計算獲得各因子與MRSEI的關聯系數都未通過統計學最低顯著性檢查(R4, 0.05=0.811),說明用LDI替代LST后保持了各評價因子對MRSEI的相互獨立狀態(表2),滿足PCA分析n維向量投影條件。從PCA結果來看,PC1軸占特征值信息的68.98%,其中NDVI具有最大向量投射(0.8901)是決定因子,是影響PC1最大的因子,LDI在PC1軸上向量投射居第二(-0.4146)。其與NDVI比值為-0.4770,作用方向相反,是影響PC1軸性質的次要因子。因此,PC1軸在MRSEI中的生態學解釋是NDVI和LDI是影響生態環境質量的主要特征分量,其值越大生態環境質量越好,這種作用在PC1特征值中居主導作用;相反LDI越大生態環境質量越差,但這種作用在PC1特征值中居次要作用。對現狀的解釋為在NDVI高值區,LDI往往處于低水平。其生態學意義為大面積高NDVI農田和喬木林地分布區對應地是LDI低值區,這樣的區域在MRSEI計算中分值較高。PC2軸占特征值信息的28.76%,其中LDI具有最大向量投射(0.9100)是決定因子,LDI是影響PC2最大的因子,NDVI在PC2軸上向量投射居第二(0.4056)與LDI比值為0.4457,作用與LDI方向一致是影響PC2軸性質的次要因子。因此,PC2軸在MRSEI中的生態學解釋就是NDVI和LDI是影響生態環境質量的次要特征分量,NDVI和LDI越大生態環境質量越好,反之亦然。對現狀的解釋為在LDI高值區,NDVI往往處于低水平,這樣的區域在MRSEI計算中分值較低。其生態學意義為若任一分析單元相鄰有多種土地利用/覆蓋類型,其植被生產力對MRSEI貢獻被弱化。結合研究區整體MRSEI水平為0.0117略超過“較好”級別低限(0.0116),說明在研究區中NDVI和LDI高值區重疊現象發生概率較低。從MRSEI空間異質性特點來看,一是以“差”和“較差”級別多分布在賀蘭山東麓—寧夏平原過渡帶、鄂爾多斯臺地—寧夏平原過渡帶和黃土高原—寧夏平原過渡帶三個區域,呈環狀圍繞整個研究區。二是格局指數斑塊密度和平均斑塊面積與MRSEI變化具有較好的規律性,主要表現為在從“差”到“好”梯度上,斑塊密度為減少趨勢,而平均斑塊面積呈增加趨勢(表4)。總的來看,LDI替代LST進行的MRSEI分析首先保證了與原有三個因子(NDVI、WI和NDBSI)的低關聯度,其中與NDVI最高僅為-0.3221(表2),沒有發生對MRSEI造成重復貢獻表達的結果。同時,基于MRSEI分級中值閾值評價寧夏沿黃平原區生態環境質量總體上處于較好水平。從具體應用角度來看,在ArcMap環境中對研究區內任一指定范圍(或行政區)的總體生態環境質量水平評估可通過表2中值閾值和選定范圍內各級別相對面積用公式(3)完成計算。

在MRSEI分析中研究區自然特征必須予以高度重視。在大多數情況下采用LDI替代LST在評價因子屬性代表性和生態學解釋上是可行的,因為其有效地避免了RSEI分析中評價因子NDBSI和LST之間存在的生態學意義重復表達和多因子向量投影中的高度聚集。但是在特殊地理單元如青藏高原等常年凍土分布區和不同沙丘(固定、半固定等)類型為主的地貌類型區域,LST作為極其重要的指標是否適用需要慎重考慮。但是LST所反映的短時間尺度屬性必須通過恰當方法(多數據統計等)改進,增強其與其它評價因子在時間尺度代表性上相對匹配。

5 結論

在寧夏沿黃平原區基于土地利用/覆蓋分類的景觀多樣性指數具有顯著的尺度依賴特征(P<0.001),閾值出現在3000 m×3000 m。PCA解釋了研究區生態環境質量主要受到NDVI和LDI影響,其中NDVI是PC1(特征值貢獻率68.98%)的決定因子,特征向量為0.8901;LDI為次要決定因子,特征向量為-0.4146,該類型區在MRSEI計算中分值較高。LDI是PC2(特征值貢獻率28.76%)的決定因子,特征向量為0.9100;NDVI為次要決定因子,特征向量為0.4056。研究區空間異質性主要以“差”和“較差”級別分布在不同土地利用/覆蓋類型交錯區且呈環繞研究區為主要特點。在MRSEI由“差”到“好”梯度上,數量特征主要表現為斑塊密度為減少趨勢由8.3 個/km2減少到5.9 個/km2,而平均斑塊面積呈增加趨勢由0.120 km2增加到0.169 km2。綜合來看寧夏沿黃平原區生態環境質量總體MRSEI值為0.0117處在較好水平。