基于RSEI的山東省不同土地覆蓋區生態環境質量變化特征及驅動因素

馬林,梁勇,張雅芹,李云鵬

(1.濟南市土地儲備中心,濟南 250099;2.濟南市勘察測繪研究院,濟南 250101)

我國工業化和城市化進程的不斷推進,給人類帶來了很多福祉同時也使很多區域面臨著各種各樣的生態問題,對區域經濟及社會的可持續發展也產生了顯著影響,目前已成為政府及學術界關注的焦點[1-2]。山東省作為我國經濟強省和全球變化較為敏感的生態系統地帶,在過去30 a的發展過程中土地覆蓋類型發生了巨大變化[3-4]。因此,定量描述和評價城鎮化和氣候雙重影響下山東省生態變化特征及其變化的主要驅動因子對指導當前生態文明建設工作具有重要作用。

目前國內外學者對山東省生態環境監測與評估的大多數研究均使用單一生態因子[3-5]。如孫匯穎等[3]基于土地利用類型指標發現2000—2018年山東省生境環境質量主要受到土地利用變化的影響;劉玉紅等[4]基于植被初級生產力(NPP)指數發現山東省2000—2015年NPP呈增加趨勢,降水量是影響植被NPP變化的主導因子;莊會波等[5]基于歸一化植被指數(NDVI)指標發現山東省1998—2018年植被覆蓋度呈增加趨勢,但不同地形區域存在差異;也有很多學者采用地表溫度(LST)、歸一化差異不透水表面指數(NDBSI)、裸土指數(BI)分別探究了山東省城市熱島效應、不透水表面面積及裸土區域等變化情況從而反映山東省生態環境質量等[6-9]。以上研究雖可以在一定程度上反映山東省生態環境質量情況,但由于山東省生態系統的復雜性和多樣性,僅采用一個或兩個生態指標來評估生態系統的狀況存在片面性。因此有學者研究發現通過主成分分析(PCA)方法整合綠色、濕度、干燥和熱量生成的完全基于遙感的生態指數(RSEI)可以實現快速、高精度的監測區域生態環境質量及其變化情況[10-14]。RSEI不僅可以反映區域植被覆蓋、土壤濕度狀況,還可以反映地表干旱程度和地表熱度情況,能對區域生態環境狀況進行客觀定量評價,成為區域生態環境質量監測和評估的新方法[15-19]。如張靜等[10]基于RSEI對西安市生態環境質量監測發現RSEI可以很好反映西安市生態環境質量;羅镕基等[11]基于改進的RSEI對甘肅省古浪縣生態環境質量監測發現RSEI能較好地綜合各分量指標信息全面準確反映古浪縣生態質量狀況;Yang等[20]也基于RSEI對長江流域生態環境變化特征及其驅動因素發現長江流域整體生態質量有所改善,RSEI能較好地反映區域生態環境質量狀況。

因此,本文基于Landsat數據計算NDVI(綠色度)、WET(濕度)、LST(熱量)和歸一化差異累積和裸土指數(NDBSI)(干燥度),然后基于PAC 法獲得1991—2021年山東省RSEI時間序列數據,并基于RSEI時序數據探究山東省不同土地覆蓋區域生態環境質量空間分布特征、變化趨勢及環境因素和土地利用轉移變化對RSEI的驅動情況,以期為山東省生態環境保護戰略與經濟可持續發展戰略的協調發展提供理論依據。

1 材料與方法

1.1 研究區概況

山東省地處我國東部沿海,黃河下游,瀕臨渤海與黃海(114.79°—122.7°,34.38°—38.4°),省內地勢中部多山地,西南、西北地勢低洼平坦,東部丘陵起伏,總面積約155 800 km2。氣候屬于暖溫帶濕潤半濕潤季風性氣候,植被類型以溫帶闊葉林、溫帶針葉林和溫帶針葉闊葉混交林為主,全年平均降水量達到710 mm,年平均氣溫11～14℃,年均日照時長為2 290～2 890 h,降水集中,雨熱同季。山東省作為我國主要的糧食產區之一,全省農業總產值連續位居全國第一位,截至2022年底,全省人口1.02億,生產總值為87 435億元[8-9]。山東省地處三大經濟圈的交匯處,連接了華東與華中地區,經濟發展迅速,土地利用變化對區域生態環境質量產生了顯著的影響(圖1)。

圖1 山東省土地類型空間分布Fig.1 Spatial distribution of land types in Shandong Province

1.2 數據源及預處理

(1)Landsat數據。本文使用Google Earth Engine(GEE)平臺提供1991—2021年的30 m 空間分辨率的Landsat系列數據提供的地表反射率數據構建RSEI指數需要的歸一化差異植被指數(NDVI)、纓帽變換的濕度分量(WET)、歸一化差異不透水表面指數(NDBSI)和地表溫度(LST)。GEE 平臺提供的Landsat的T2_SR 產品已經過輻射定標、大氣校正等預處理,在構建RSEI時通過多次對比試驗最終選擇每年的3—10月份的數據構建每年的RSEI,因為該階段能更好地體現RSEI指數的特性。

(2)土地利用數據(LUCC)。本文土地利用數據來源于中國科學院資源環境科學數據平臺(http:∥www.resdc.cn),該數據集是由多個研究所基于Landsat系列的遙感影像結合人機交互目視解譯方法獲取,其二級分類包括草地、林地、耕地、建設用地、水域和未利用地,分辨率為30 m。其中耕地和建設用地平均分類精度達到85%以上,其他土地利用類型平均分類精度達到75%以上,可滿足本文研究需求。本文采用1990年、2010年、2020年的土地利用數據。

1.3 研究方法

1.3.1 遙感生態指數 在反映生態質量的諸多自然因素中,綠度、濕度、熱度、干度可謂是與人類生存息息相關的4個重要指標,也是人類直觀感覺生態條件優劣的重要因素,因此常被用于評價生態系統。因此徐涵秋等[21]提出的遙感生態指數(RSEI)從反映生態質量的眾多自然因素中選擇與人類生存相關的4個因素作為生態評價因素,該指數利用主成分分析耦合了表示地表植被覆蓋和生長情況的綠色指數(NDVI)、表示土壤的濕度狀況的濕度(WET)指數、表示地表干旱程度的干燥度指數(NDBSI)和反映地表熱度情況的熱量指數(LST)[16-18]。其中,NDVI和NDBSI指標可以反映生態環境對人類活動引起的土地覆蓋變化的生態響應,LST 和WET 揭示了生態質量對代表氣候變化的響應,4個生態指標綜合起來能夠直觀地反映區域生態環境惡劣情況[21],4個指標的計算如下:

(1)綠度指標。NDVI是應用最廣泛的植被指數之一,許多生態學研究均證明了其有用性[12-13]。因此,本文選擇NDVI作為綠度指數,計算公式為:

式中:ρnir為遙感圖像的近紅外波段;ρred為紅光波段。

(2)熱度指標。熱度指標由地表溫度(LST)表示,雖然Landsat 8 TIRS傳感器中有兩個熱紅外波段,但由于波段11的輻射校準問題,本文僅選擇了波段10來反演LST[21]。當然,對于TM,我們仍然得到了基于波段6的LST。首先,將數字(DN)值轉換為傳感器孔徑處的光譜輻射亮度(Lλ),計算公式見(2):

式中:Grescale為特定于頻帶的重新縮放增益因子;Qcal為量化的校準像素值(DN);Brescale為不同于頻帶的重縮放偏置因子。傳感器處光譜輻射亮度(Lλ)可使用公式(3)轉換為傳感器亮度溫度下的有效輻射亮度。

式中:Tλ為傳感器亮度溫度下的有效值,K1〔TM 波段6為607.76 W/(m2·sr·μm),TIR波段10為774.89 W/(m2·sr·μm)〕和K2(TM 波段4為260.56 K,TIR波段10為1 321.08 K)分別是校準常數1,2。最后,使用Tλ獲得LST。

式中:λ為發射輻射的波長(TM 波段6 為11.435 μm,TIR 波段10 為10.9μm);ρ為常數(1.438×10-2m K);ε為表面發射率,可通過使用Sobrino模型的NDVI估算[21]。

(3)濕度指標。對地表濕度信息較為敏感的纓帽變換濕度分量(TCW)已被證明對植被、水和土壤濕度條件較為敏感[11]。因此,本文選擇TCW 作為濕度指數(WET),并且可以分別使用方程(5)和(6)來估計Landsat TM 和Landsat 8的濕度指數[12]。

式中:ρi(i=1,2,3,…,6)分別為TM 和OLI傳感器對應的藍光波動、綠光波段、紅光波段、近紅外波段、短波紅外波段1和短波紅外波段2。

(4)干度指標。由于建設用地和裸露土壤會對區域生態環境造成嚴重危害。因此,本文選擇了歸一化差異不透水表面指數(NDBSI)結合基于指數的建成指數(IBI)和土壤指數(BSI)表征研究區干燥度[21]。

式中:S1,ρred,ρgreen,ρnir,ρbule分別為TM 和OLI對應的短波紅外、紅光、綠光、近紅外和藍光波段。

(5)構建RSEI。通過主成分變化(PCA)對以上4個指標空間坐標軸旋轉從而去掉各指標間的相關性,把4個指標的主要信息集中到1～2個主成分上,從而構建RSEI。這種方法最大優點是綜合指標的權重不是人為確定的,而是根據每個指標對PC1 的貢獻自動客觀地確定。因此,在計算中避免了由于權重設置而導致的結果的可能偏差,該權重設置因個人和方法而異,這大大提高了結果的客觀性和可信度。在預處理后需要對4個指標進行歸一化以實現量綱統一,然后采用PCA 方法將NDVI,LST,WET,NDBSI的信息集中到第一、二主成分上[21]。

指標歸一化公式見(10):

式中:NIi為各指標歸一化結果;Ii為各指標在第i像元的值;Imin和Imax分別為各指標最小值和最大值。

為了使PC1大的數值代表好的生態條件,進一步用1減去PC1,獲取初始RSEI,RSEI計算公式見公式(11)—(12):

為了便于指標的度量和比較,同樣對RSEI0進行歸一化:

式中:RSEI0max,RSEI0min分別為RSEI0最大值和最小值。RSEI值越接近1,生態越好,反之,生態越差。

本文基于文獻[14]將RSEI劃分為RSEI＜0.2(差),0.2≤RSEI＜0.4(較差),0.4≤RSEI＜0.6(中等),0.6≤RSEI＜0.8(良好)和RSEI≥0.8(優)5個等級。

1.3.2 RSEI變化趨勢分析 采用Sen趨勢法探究1991—2021年山東省RSEI變化趨勢和變化速率。Sen斜率的計算公式如下[22]:

式中:β為RSEI變化趨勢;Median為取中值函數;當β＞0時,RSEI呈增加趨勢;β＜0時,RSEI呈減小趨勢。

1.3.3 穩定性分析 研究使用變異系數Cv分析RSEI近31 a的持續的穩定性,Cv越大表明RSEI的變化越不穩定[20]。

式中:Cv為變異系數;Xi為第i年對應的RSEI值;Xmean為基于1991—2021年的RSEI均值。本文根據自然斷裂法將RSEI的穩定性分為5類:低波動變化(＜0.049)、較低波動變化(0.049～0.1)、中波動變化(0.1～0.185)、較高波動變化(0.185～0.392)和高波動變化(0.392～1)。

1.3.4 未來趨勢分析 Hurst指數是根據RSEI的長時間序列數據預測未來數據的變化趨勢,基于重標極差(R/S)的Hurst指數是一種定量表述時間序列信息長期依賴性的有效方法[23]。本文基于Hurst指數研究RSEI未來的趨勢,若0≤Hurst＜0.5,說明RSEI將來總體趨勢和過去趨勢相反;若Hurst=0.5,表明RSEI互相獨立且沒有依賴性,現在不會影響未來;0.5＜Hurst≤1,表明未來RSEI總體趨勢與過去總體趨勢相同,越接近1,持續性越強[23]。

1.3.5 地理探測器 本文利用地理探測器度量RSEI的空間分布特性,以及探測不同因子與RSEI的相關性程度,用q值解釋這個程度[24]。公式如下:

式中:q為影響因子對RSEI時空變化的解釋力;h為不同驅動因子的分類或分區數據;L為影響因子的樣本數量;Nh和N分別為h和整個區域的單元數;δ2h和δ2為h和整個區域的方差。q值越大,代表該因子對RSEI影響程度越大。并使用交互作用探測識別不同因子之間的相互作用,并評估他們的組合效應以觀測任何一對因素共同作用是否會增加或降低對RSEI空間分布的解釋力[24]。

2 結果與分析

2.1 生態環境整體分析

由山東省RSEI主成分分析結果(表1)可知,1991—2021年的第一主成分貢獻率均達到80%以上,說明第一主成分集中了NDVI,LST,WET 和NDBSI指標的大部分特征信息,能夠較好地展示該區域內的生態環境特征情況。從各指標的載荷來看,NDVI和WEI的載荷值成正比,表示綠度和濕度越高,土壤水分充足,地表植被覆蓋越高,生態環境越好。綠度和濕度越低,土壤水分不足,地表植被覆蓋降低,生態環境惡劣。LST 與NDBSI的載荷值成正比,干度和熱度越高,表示土壤沙化、荒漠化、巖石裸露等問題越嚴重,干度和熱度越低,表示土壤沙化、荒漠化、巖石裸露等問題得到緩解。

表1 各指標第一主成分及貢獻率Table 1 First principal component and contribution rate of each indicator

2.2 各指標時序變化特征

由圖2可知,RSEI總體以每年0.000 5的速率減小,即生態環境質量輕微退化。NDVI每年以0.003 1的速率增加,WET,LST,NDBSI分別以0.002/a,0.008/a,0.008 6/a的速率下降。雖然具有正效益的NDVI指標呈增加趨勢和具有負效益的LST,NDBSI呈減小趨勢,但具有正效益的WET 卻呈減小趨勢,經過主成分分析的4 個遙感指標在抵消了部分生態效益的不同影響后,將導致RSEI指標下降,且從各指標曲線走勢來看,RSEI曲線走勢與4種指標走勢均不一致,由此可見RSEI變化是4 種指標共同變化引起。

圖2 4個指標和RSEI時序曲線Fig.2 Four indicators and RSEI time series curve

2.3 RSEI的空間分布特征

從山東省RSEI的空間分布來看(圖3),1991年山東省威海市、煙臺市、濟寧與臨沂市交界處,濟寧市與秦安市交界處生態環境質量表現為優等級,聊城南部與菏澤北部生態環境質量表現為較差等級。1996年生態環境質量較1991年整體有所改善,大部分地區表現良好等級和優等級,濟寧東南部表現為較差等級。2001年山東省大部分地區生態質量環境表現較差,其中煙臺、威海、青島南部表現為差等級。2006年相較于2001年生態環境質量有明顯提升,山東省中部、東部、北部表現為良好等級,南部部分地區表現為較差等級。2011年相較于2006年,山東省整體生態環境質量有惡化趨勢,其中德州、聊城、濟南、濱州等地區大部分表現為差等級,菏澤、淄博、臨沂大部分地區表現為較差等級。2016年山東省整體生態環境質量有較大改善,大部分地區表現為良好等級,相較于2011年,生態環境質量有顯著提升。2021年山東省生態環境質量嚴重下滑,其中聊城、德州表現為差等級,濟南北部、濱州中南部、淄博北部表現為較差等級,煙臺、青島和濰坊東南部表現為較差等級。從生態環境特征總體來看,1991—2020年山東省生態環境質量整體表現為下降趨勢。

圖3 RSEI的空間分布Fig.3 Spatial distribution map of RSEI

圖4顯示了每隔5 a的RSEI統計結果,1991—2021年RSEI的結果值在0～1范圍內均有分布,均值在1991—2001 年為下降趨勢,2006—2021 年呈“增—減—增—減”趨勢。從分布來看,2001年RSEI范圍分布集中在0.1～0.6,分布相對集中且低于其他年份,均值為7 個年份中最低年份(0.37);2006 年RSEI集中在0.4～0.8,分布集中且高于其他年份,均值為7個年份中最高(0.59)。從不同等級占比來看,良等級占比與RSEI均值變化特征一致,1991—2001年減少,2006—2021年呈“增—減—增—減”趨勢,但較差等級變化特征與良等級呈相反趨勢,中等級占比保持穩定不變。

圖4 RSEI的空間統計Fig.4 Spatial statistics of RSEI

2.4 RSEI變化趨勢

由圖5,圖6A 可知,不同時間段山東省RSEI變化為穩定狀態的區域占主體地位,其他趨勢具有時空差異。其中,1991—2000年14.33%的區域RSEI呈上升趨勢,集中在日照市、臨沂市和菏澤市南部地區,19.66%的區域RSEI呈下降趨勢,集中在泰安市、濟南市、萊蕪市、淄博市。2001—2010年RSEI呈增加上升趨勢的區域較1991—2000年增加了30.69%,轉變區域主要集中在山東省西部大部分城市和東北部的威海市和煙臺市。2011—2021 年73.6%的區域RSEI變化呈穩定狀態,僅有13.58%和12.82%的區域RSEI表現為上升和下降趨勢,呈上升趨勢的區域零星分布在中部部分城市,呈下降趨勢的區域零星分布在東部和東北部部分城市。

圖5 各階段RSEI的趨勢分布Fig.5 Trend distribution of RSEI in each stage

圖6 各階段RSEI變化趨勢和穩定性所占像元百分比Fig.6 Statistical chart of the changes trends of RSEI and stability at each stage

基于1991—2000年、2001—2010年、2011—2021年山東省生態環境質量變異性計算結果和實際情況,將變異劃分為5個等級〔變化相對穩定(≤0.2)、低波動變化(0.2～0.25)、較低波動變化(0.25～0.3)、較高波動變化(0.3～0.4)、高波動變化(＞0.4)〕(圖6B、圖7)。1991—2000年RSEI變異系數屬于變化相對穩定、低波動變化和較低波動變化的區域占比為69.859%,較高波動變化和高波動變化占30.141%,較高波動變化的區域主要分布在山東省北部和東部地區。2001—2010年RSEI呈變化相對穩定、低波動變化和較低波動變化的區域占比為60.354%,較高波動變化和高波動變化占39.646%,其中呈高波動變化的區域較1991—2000年增加9.248%,這些區域集中在東北部的威海市、煙臺市和東部的青島市和日照市和中部的泰安市,這些區域主要土地利用類型為草地、林地和建設用地。2011—2021年僅有25.582%的區域為較高和高波動變化,集中在西部的聊城市和德州市,其他區域均以低波動變化為主。

圖7 各階段RSEI的變化穩定性分布Fig.7 Stability distribution of RSEI changes in different stages

從不同土地利用覆蓋區域RSEI變化速率和變異性可知(表2),1991—2000年除草地區域RSEI未變化,其他區域RSEI均為下降趨勢,且建設用地區域RSEI下降速率最快(0.091/a),其次為林地(0.070/a)。但從RSEI變異性來看,不同土地利用區域RSEI均為較低波動變化,且草地RSEI波動性＞未利用地＞減少用地＞耕地＞林地。2001—2010年不同土地利用區域RSEI均為增加趨勢,且林地區域RSEI增加速率最大(0.158/a),其次為耕地區域(0.155/a),未利用地RSEI增加速率最小(0.053/a)。但從變異性來看,未利用地RSEI波動性＞草地＞耕地＞建設用地＞林地。2011—2021年林地和耕地區域均為增加趨勢,增加速率分別為0.021/a,0.019/a,其他土地利用類型均為減小趨勢,且未利用地減小速率最大(0.080/a)。從變異性來看,耕地RSEI變異性最大(變異系數為0.284),其次為林地(變異系數為0.270),草地變異性最小(變異系數為0.236),其次為未利用地(變異系數為0.237)?？傮w來看,2001—2010年RSEI變化較大,且波動性也為3個時間段最大。

表2 不同土地利用類型下各階段RSEI變化速率和變異性Table 2 Statistical chart of RSEI rate of change and variability in different land use types at different stages

2.5 RSEI未來趨勢分析

從圖8可知,山東省Hurst指數(H)分布以0.3～0.4區間為主,H＞0.5的區域僅占0.77%。零星分布在山東省煙臺市和威海市,其他區域H均小于0.5,表明未來一段時期內山東省RSEI的變化趨勢與1991—2021年變化趨勢為相反趨勢。其中反持續性顯著減少(0≤H＜0.4)占79.82%,反持續性減少不顯著(0.4≤H＜0.5)占19.41%。從不同土地類型區域RSEI的H來看,不同土地利用類型區域未來一段時間內變化特征均與過去呈相反趨勢,且不同土地類型區域的RSEI均小于0.4,這表明不同土地類型區域的反持續性均較顯著。

圖8 RSEI的Hurst指數分布Fig.8 Distribution map of hurst index for RSEI

2.6 RSEI對驅動因素的響應

2.6.1 土地覆蓋變化 從不同階段山東省土地利用轉移矩陣可知(圖9,表3):1990—2000年耕地和未利用面積以減少為主,分別減少了1 111,941 個像元,而水域和建設用地面積以增加為主,分別增加了461,318個像元,發生轉移的區域零星分布在整個山東省。2000—2010年耕地和未利用地仍為大面積減少,分別減少了1 808,523個像元,建設用地面積大幅度增長,增加了2 791,發生土地利用轉移的區域也零星分布在山東省整個區域,且在濱州市最北部建設用地顯著增加。2010—2020年不同土地利用類型均發生了較大的轉移,其中以草地和未利用地大面積減少、建設用地大面積增加為主,其中在濱州市北部、濰坊市北部等水域區域的建設用地轉移為水域,中部的濟南市、淄博市和棗莊市濟寧市交匯處和東北部的煙臺市新增草地面積大于減少草地面積。

表3 各階段山東省土地利用類型轉換像元統計Table 3 Pixel statistics of land use type conversion in Shandong Province at different stages

圖9 各階段山東省土地利用類型的轉換變化Fig.9 Transformation change of land use types in Shandong Province at different stages

2.6.2 不同階段4個生態指標對RSEI交互作用探測結果 為了探究不同環境因子對山東省不同時間段內RSEI空間分布的影響,采用地理探測器分析不同階段不同因子對RSEI的解釋力(表4)。從表4—5可知,1990—2000年、2001—2010年LUCC對RSEI空間分布變化的解釋力最強,q值分別為0.184 2,0.184 0,但與2001—2010年相比,WET 和NDVI對RSEI的影響較大。2011—2021年RSEI空間分布變化受到WET 的影響最大,其次為LUCC,NDBSI和NDVI對其影響較小。

表4 不同因子對不同階段RSEI變化的q 驅動力值Table 4 q driving force values of different factors on RSEI changes at different stages

從交互探測結果可知(表5),不同階段環境因子的交互作用產生的q值均大于單獨因子的q值,即不同環境因子的交互作用會增加對RSEI空間分布變化的解釋。其中1990—2000 年LUCC 與其他因子之間的交互作用對RSEI的影響均大于其他因子彼此間的交互作用,解釋率均在30%以上,LST 與其他因子間對RSEI的解釋力均在15%以上。2001—2010年LUCC 與其他因子之間的交互作用對RSEI也均大于其他因子之間的交互作用,解釋力均在30%以上。2011—2021年RSEI的空間分布變化主要受到LST∩LUCC,WET∩NDBSI,NDVI∩NDBSI和WET∩LUCC 的交互作用,其對RSEI空間分布變化的解釋力均在10%以上,其他因子之間的交互作用對RSEI的影響較小。

表5 不同因子對不同階段RSEI變化的交互作用Table 5 Interaction diagram of different factors on RSEI changes at different stages

3 討論與結論

山東省作為我國經濟強省和農業大省,其生態的好壞直接關系到人民生活水平。但隨著城鎮化快速發展,城市用地的快速擴張,導致其他土地利用覆蓋類型發生了顯著變化(尤其是生態用地)[3-6]。雖然山東省實施了一系列生態文明建設(如生態優先、綠色低碳發展、黃河流域生態保護與高質量發展策略等),但山東省生態環境質量對地形地貌、氣候變化和人類活動的響應過程是一個復雜的動態過程,不同地貌條件會導致熱量和水分的空間差異,進而使得生態環境質量變化存在明顯空間差異[7]。因此,亟需厘清生態工程實施和城鎮化快速發展背景下,山東省生態環境質量時空演變及其主要驅動因子。因此,本文利用1991—2021年的Landsat數據計算RSEI指數,并輔以Sen、變異系數法、Hurst法、隨機森林對山東省生態環境質量變化特征及其變化的驅動因子進行分析,主要結論如下:

(1)1991—2021年山東省生態環境質量總體以下降趨勢為主,且以每年0.000 5的速率減小,NDVI每年以0.003 1 的速率增加,WET,LST,NDBSI分別以0.002/a,0.008/a,0.008 6/a 的速率下降。1991—2000年在日照市、臨沂市和菏澤市南部地區呈改善趨勢,泰安市、濟南市、萊蕪市、淄博市為下降趨勢;2001—2010年RSEI呈增加趨勢的區集中在山東省西部和東北部;2011—2021年RSEI變化大部分區域呈穩定狀態。孫匯穎等[3]基于In Vest對山東省2000—2018年生境質量研究同樣也發現近20 a山東省生境環境質量總體為下降趨勢;莊會波等[5]基于NDVI指數研究了山東省1982—2018年植被覆蓋度變化特征,發現近30 a山東省NDVI均呈增加趨勢。以上研究結論與本文結論相一致,這表明山東省生態環境質量仍在惡化,需要大力推進生態文明建設和城市化綠色發展理念從而改善山東生態環境質量。

(2)1991—2000年山東省生態環境質量變化以穩定趨勢為主,變化以低波動變化為主,且除草地區域生態環境質量為上升趨勢,其他區域均為下降趨勢,但草地區域生態環境質量變化波動性＞未利用地＞減少用地＞耕地＞林地。2001—2010年生態環境質量仍以穩定變化趨勢為主,但在山東省西部大部分城市和東北部的威海市和煙臺市以增加趨勢為主,且不同土地利用區域生態環境均為增加趨勢,未利用地波動性＞草地＞耕地＞建設用地＞林地。2011—2021年生態環境質量變化呈穩定狀態,但林地和耕地區域生態環境質量為增加趨勢,其他土地利用類型均為減小趨勢,耕地波動性＞林地＞減少用地＞未利用地＞草地。

(3)山東省生態環境質量未來一段時期內的變化趨勢以增加趨勢為主,且不同土地利用類型區域未來一段時間內變化特征均與過去呈相反趨勢。

(4)1990—2000 年、2001—2010 年土地利用變化對山東省生態環境質量空間分布變化的影響最大,2010—2021年生態環境質量空間分布變化受到濕度變化和土地利用變化的影響較大。不同階段不同環境因子的交互作用會增加對生態環境質量的影響。1990—2000年、2001—2010年土地利用與其他因子之間的交互作用對生態環境質量的影響均較大,2011—2021年生態環境質量的空間分布變化主要受到地表溫度和土地利用、濕地和干度、綠度和干度以及土地利用和濕度的交互作用。

20世紀以來,隨著山東省新舊動能轉換重大工程實施,城市用地面積快速擴張,青連高速、魯南高速等交通干線相繼開通,山東省土地利用空間格局發生了巨大變化[25]。土地利用變化是人類活動最直接的反映,本研究發現土地利用是驅動1991—2010年山東省生態環境質量發生變化的主要驅動因子,同樣本文也發現這一階段山東省建設用地面積快速擴張,而作為特殊的半人工自然生態系統的耕地面積快速收縮,因此造成這一階段生態環境遭到破壞。這一發現已在吳健生[26]、李勝鵬[27]、孫匯穎等[3]的研究結果中得以證明。近年來,隨著人類我國生態文明建設戰略的實施以及人們對生態意識的增強,生態環境質量較前20 a有所改善,且2010—2021年之后土地利用變化對山東省生態環境質量的影響較小[8]。但目前山東省生態環境質量仍以中等及以下等級為主,因此仍需要進一步加強生態文明建設,在城市群中增加能夠提升生態環境質量的綠色植物,如提高交通道路邊、小區內部植被覆蓋度等,從而增加城市綠地率,構建國土空間安全格局,改善山東省生態環境質量。