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

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

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

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

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

1 材料與方法

1.1 研究區(qū)概況

山東省地處我國東部沿海,黃河下游,瀕臨渤海與黃海(114.79°—122.7°,34.38°—38.4°),省內(nèi)地勢中部多山地,西南、西北地勢低洼平坦,東部丘陵起伏,總面積約155 800 km2。氣候?qū)儆谂瘻貛駶櫚霛駶櫦撅L(fēng)性氣候,植被類型以溫帶闊葉林、溫帶針葉林和溫帶針葉闊葉混交林為主,全年平均降水量達到710 mm,年平均氣溫11～14℃,年均日照時長為2 290～2 890 h,降水集中,雨熱同季。山東省作為我國主要的糧食產(chǎn)區(qū)之一,全省農(nóng)業(yè)總產(chǎn)值連續(xù)位居全國第一位,截至2022年底,全省人口1.02億,生產(chǎn)總值為87 435億元[8-9]。山東省地處三大經(jīng)濟圈的交匯處,連接了華東與華中地區(qū),經(jīng)濟發(fā)展迅速,土地利用變化對區(qū)域生態(tài)環(huán)境質(zhì)量產(chǎn)生了顯著的影響(圖1)。

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

1.2 數(shù)據(jù)源及預(yù)處理

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

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

1.3 研究方法

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

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

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

(2)熱度指標。熱度指標由地表溫度(LST)表示,雖然Landsat 8 TIRS傳感器中有兩個熱紅外波段,但由于波段11的輻射校準問題,本文僅選擇了波段10來反演LST[21]。當然,對于TM,我們?nèi)匀坏玫搅嘶诓ǘ?的LST。首先,將數(shù)字(DN)值轉(zhuǎn)換為傳感器孔徑處的光譜輻射亮度(Lλ),計算公式見(2):

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

式中: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)分別是校準常數(shù)1,2。最后,使用Tλ獲得LST。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2 結(jié)果與分析

2.1 生態(tài)環(huán)境整體分析

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

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

2.2 各指標時序變化特征

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

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

2.3 RSEI的空間分布特征

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

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

圖4顯示了每隔5 a的RSEI統(tǒng)計結(jié)果,1991—2021年RSEI的結(jié)果值在0～1范圍內(nèi)均有分布,均值在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年呈“增—減—增—減”趨勢,但較差等級變化特征與良等級呈相反趨勢,中等級占比保持穩(wěn)定不變。

圖4 RSEI的空間統(tǒng)計Fig.4 Spatial statistics of RSEI

2.4 RSEI變化趨勢

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

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

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

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

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

從不同土地利用覆蓋區(qū)域RSEI變化速率和變異性可知(表2),1991—2000年除草地區(qū)域RSEI未變化,其他區(qū)域RSEI均為下降趨勢,且建設(shè)用地區(qū)域RSEI下降速率最快(0.091/a),其次為林地(0.070/a)。但從RSEI變異性來看,不同土地利用區(qū)域RSEI均為較低波動變化,且草地RSEI波動性＞未利用地＞減少用地＞耕地＞林地。2001—2010年不同土地利用區(qū)域RSEI均為增加趨勢,且林地區(qū)域RSEI增加速率最大(0.158/a),其次為耕地區(qū)域(0.155/a),未利用地RSEI增加速率最小(0.053/a)。但從變異性來看,未利用地RSEI波動性＞草地＞耕地＞建設(shè)用地＞林地。2011—2021年林地和耕地區(qū)域均為增加趨勢,增加速率分別為0.021/a,0.019/a,其他土地利用類型均為減小趨勢,且未利用地減小速率最大(0.080/a)。從變異性來看,耕地RSEI變異性最大(變異系數(shù)為0.284),其次為林地(變異系數(shù)為0.270),草地變異性最小(變異系數(shù)為0.236),其次為未利用地(變異系數(shù)為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指數(shù)(H)分布以0.3～0.4區(qū)間為主,H＞0.5的區(qū)域僅占0.77%。零星分布在山東省煙臺市和威海市,其他區(qū)域H均小于0.5,表明未來一段時期內(nèi)山東省RSEI的變化趨勢與1991—2021年變化趨勢為相反趨勢。其中反持續(xù)性顯著減少(0≤H＜0.4)占79.82%,反持續(xù)性減少不顯著(0.4≤H＜0.5)占19.41%。從不同土地類型區(qū)域RSEI的H來看,不同土地利用類型區(qū)域未來一段時間內(nèi)變化特征均與過去呈相反趨勢,且不同土地類型區(qū)域的RSEI均小于0.4,這表明不同土地類型區(qū)域的反持續(xù)性均較顯著。

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

2.6 RSEI對驅(qū)動因素的響應(yīng)

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

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

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

2.6.2 不同階段4個生態(tài)指標對RSEI交互作用探測結(jié)果 為了探究不同環(huán)境因子對山東省不同時間段內(nèi)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 驅(qū)動力值Table 4 q driving force values of different factors on RSEI changes at different stages

從交互探測結(jié)果可知(表5),不同階段環(huán)境因子的交互作用產(chǎn)生的q值均大于單獨因子的q值,即不同環(huán)境因子的交互作用會增加對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 討論與結(jié)論

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

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

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

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

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

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