基于RSEI的山東省不同土地覆蓋區(qū)生態(tài)環(huán)境質(zhì)量變化特征及驅(qū)動(dòng)因素

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

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

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

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

因此,本文基于Landsat數(shù)據(jù)計(jì)算NDVI(綠色度)、WET(濕度)、LST(熱量)和歸一化差異累積和裸土指數(shù)(NDBSI)(干燥度),然后基于PAC 法獲得1991—2021年山東省RSEI時(shí)間序列數(shù)據(jù),并基于RSEI時(shí)序數(shù)據(jù)探究山東省不同土地覆蓋區(qū)域生態(tài)環(huán)境質(zhì)量空間分布特征、變化趨勢及環(huán)境因素和土地利用轉(zhuǎn)移變化對RSEI的驅(qū)動(dòng)情況,以期為山東省生態(tài)環(huán)境保護(hù)戰(zhàn)略與經(jīng)濟(jì)可持續(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)性氣候,植被類型以溫帶闊葉林、溫帶針葉林和溫帶針葉闊葉混交林為主,全年平均降水量達(dá)到710 mm,年平均氣溫11～14℃,年均日照時(shí)長為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)濟(jì)圈的交匯處,連接了華東與華中地區(qū),經(jīng)濟(jì)發(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)平臺(tái)提供1991—2021年的30 m 空間分辨率的Landsat系列數(shù)據(jù)提供的地表反射率數(shù)據(jù)構(gòu)建RSEI指數(shù)需要的歸一化差異植被指數(shù)(NDVI)、纓帽變換的濕度分量(WET)、歸一化差異不透水表面指數(shù)(NDBSI)和地表溫度(LST)。GEE 平臺(tái)提供的Landsat的T2_SR 產(chǎn)品已經(jīng)過輻射定標(biāo)、大氣校正等預(yù)處理,在構(gòu)建RSEI時(shí)通過多次對比試驗(yàn)最終選擇每年的3—10月份的數(shù)據(jù)構(gòu)建每年的RSEI,因?yàn)樵撾A段能更好地體現(xiàn)RSEI指數(shù)的特性。

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

1.3 研究方法

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

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

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

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

式中:Grescale為特定于頻帶的重新縮放增益因子;Qcal為量化的校準(zhǔn)像素值(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)分別是校準(zhǔn)常數(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)濕度指標(biāo)。對地表濕度信息較為敏感的纓帽變換濕度分量(TCW)已被證明對植被、水和土壤濕度條件較為敏感[11]。因此,本文選擇TCW 作為濕度指數(shù)(WET),并且可以分別使用方程(5)和(6)來估計(jì)Landsat TM 和Landsat 8的濕度指數(shù)[12]。

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

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

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

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

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

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

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

為了便于指標(biāo)的度量和比較,同樣對RSEI0進(jìn)行歸一化:

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

本文基于文獻(xiàn)[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個(gè)等級(jí)。

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

式中:β為RSEI變化趨勢;Median為取中值函數(shù);當(dāng)β＞0時(shí),RSEI呈增加趨勢;β＜0時(shí),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類:低波動(dòng)變化(＜0.049)、較低波動(dòng)變化(0.049～0.1)、中波動(dòng)變化(0.1～0.185)、較高波動(dòng)變化(0.185～0.392)和高波動(dòng)變化(0.392～1)。

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

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

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

2 結(jié)果與分析

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

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

表1 各指標(biāo)第一主成分及貢獻(xiàn)率Table 1 First principal component and contribution rate of each indicator

2.2 各指標(biāo)時(shí)序變化特征

由圖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指標(biāo)呈增加趨勢和具有負(fù)效益的LST,NDBSI呈減小趨勢,但具有正效益的WET 卻呈減小趨勢,經(jīng)過主成分分析的4 個(gè)遙感指標(biāo)在抵消了部分生態(tài)效益的不同影響后,將導(dǎo)致RSEI指標(biāo)下降,且從各指標(biāo)曲線走勢來看,RSEI曲線走勢與4種指標(biāo)走勢均不一致,由此可見RSEI變化是4 種指標(biāo)共同變化引起。

圖2 4個(gè)指標(biāo)和RSEI時(shí)序曲線Fig.2 Four indicators and RSEI time series curve

2.3 RSEI的空間分布特征

從山東省RSEI的空間分布來看(圖3),1991年山東省威海市、煙臺(tái)市、濟(jì)寧與臨沂市交界處,濟(jì)寧市與秦安市交界處生態(tài)環(huán)境質(zhì)量表現(xiàn)為優(yōu)等級(jí),聊城南部與菏澤北部生態(tài)環(huán)境質(zhì)量表現(xiàn)為較差等級(jí)。1996年生態(tài)環(huán)境質(zhì)量較1991年整體有所改善,大部分地區(qū)表現(xiàn)良好等級(jí)和優(yōu)等級(jí),濟(jì)寧東南部表現(xiàn)為較差等級(jí)。2001年山東省大部分地區(qū)生態(tài)質(zhì)量環(huán)境表現(xiàn)較差,其中煙臺(tái)、威海、青島南部表現(xiàn)為差等級(jí)。2006年相較于2001年生態(tài)環(huán)境質(zhì)量有明顯提升,山東省中部、東部、北部表現(xiàn)為良好等級(jí),南部部分地區(qū)表現(xiàn)為較差等級(jí)。2011年相較于2006年,山東省整體生態(tài)環(huán)境質(zhì)量有惡化趨勢,其中德州、聊城、濟(jì)南、濱州等地區(qū)大部分表現(xiàn)為差等級(jí),菏澤、淄博、臨沂大部分地區(qū)表現(xiàn)為較差等級(jí)。2016年山東省整體生態(tài)環(huán)境質(zhì)量有較大改善,大部分地區(qū)表現(xiàn)為良好等級(jí),相較于2011年,生態(tài)環(huán)境質(zhì)量有顯著提升。2021年山東省生態(tài)環(huán)境質(zhì)量嚴(yán)重下滑,其中聊城、德州表現(xiàn)為差等級(jí),濟(jì)南北部、濱州中南部、淄博北部表現(xiàn)為較差等級(jí),煙臺(tái)、青島和濰坊東南部表現(xiàn)為較差等級(jí)。從生態(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)計(jì)結(jié)果,1991—2021年RSEI的結(jié)果值在0～1范圍內(nèi)均有分布,均值在1991—2001 年為下降趨勢,2006—2021 年呈“增—減—增—減”趨勢。從分布來看,2001年RSEI范圍分布集中在0.1～0.6,分布相對集中且低于其他年份,均值為7 個(gè)年份中最低年份(0.37);2006 年RSEI集中在0.4～0.8,分布集中且高于其他年份,均值為7個(gè)年份中最高(0.59)。從不同等級(jí)占比來看,良等級(jí)占比與RSEI均值變化特征一致,1991—2001年減少,2006—2021年呈“增—減—增—減”趨勢,但較差等級(jí)變化特征與良等級(jí)呈相反趨勢,中等級(jí)占比保持穩(wěn)定不變。

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

2.4 RSEI變化趨勢

由圖5,圖6A 可知,不同時(shí)間段山東省RSEI變化為穩(wěn)定狀態(tài)的區(qū)域占主體地位,其他趨勢具有時(shí)空差異。其中,1991—2000年14.33%的區(qū)域RSEI呈上升趨勢,集中在日照市、臨沂市和菏澤市南部地區(qū),19.66%的區(qū)域RSEI呈下降趨勢,集中在泰安市、濟(jì)南市、萊蕪市、淄博市。2001—2010年RSEI呈增加上升趨勢的區(qū)域較1991—2000年增加了30.69%,轉(zhuǎn)變區(qū)域主要集中在山東省西部大部分城市和東北部的威海市和煙臺(tái)市。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ì)量變異性計(jì)算結(jié)果和實(shí)際情況,將變異劃分為5個(gè)等級(jí)〔變化相對穩(wěn)定(≤0.2)、低波動(dòng)變化(0.2～0.25)、較低波動(dòng)變化(0.25～0.3)、較高波動(dòng)變化(0.3～0.4)、高波動(dòng)變化(＞0.4)〕(圖6B、圖7)。1991—2000年RSEI變異系數(shù)屬于變化相對穩(wěn)定、低波動(dòng)變化和較低波動(dòng)變化的區(qū)域占比為69.859%,較高波動(dòng)變化和高波動(dòng)變化占30.141%,較高波動(dòng)變化的區(qū)域主要分布在山東省北部和東部地區(qū)。2001—2010年RSEI呈變化相對穩(wěn)定、低波動(dòng)變化和較低波動(dòng)變化的區(qū)域占比為60.354%,較高波動(dòng)變化和高波動(dòng)變化占39.646%,其中呈高波動(dòng)變化的區(qū)域較1991—2000年增加9.248%,這些區(qū)域集中在東北部的威海市、煙臺(tái)市和東部的青島市和日照市和中部的泰安市,這些區(qū)域主要土地利用類型為草地、林地和建設(shè)用地。2011—2021年僅有25.582%的區(qū)域?yàn)檩^高和高波動(dòng)變化,集中在西部的聊城市和德州市,其他區(qū)域均以低波動(dòng)變化為主。

圖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均為較低波動(dòng)變化,且草地RSEI波動(dòng)性＞未利用地＞減少用地＞耕地＞林地。2001—2010年不同土地利用區(qū)域RSEI均為增加趨勢,且林地區(qū)域RSEI增加速率最大(0.158/a),其次為耕地區(qū)域(0.155/a),未利用地RSEI增加速率最小(0.053/a)。但從變異性來看,未利用地RSEI波動(dòng)性＞草地＞耕地＞建設(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變化較大,且波動(dòng)性也為3個(gè)時(shí)間段最大。

表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%。零星分布在山東省煙臺(tái)市和威海市,其他區(qū)域H均小于0.5,表明未來一段時(shí)期內(nèi)山東省RSEI的變化趨勢與1991—2021年變化趨勢為相反趨勢。其中反持續(xù)性顯著減少(0≤H＜0.4)占79.82%,反持續(xù)性減少不顯著(0.4≤H＜0.5)占19.41%。從不同土地類型區(qū)域RSEI的H來看,不同土地利用類型區(qū)域未來一段時(shí)間內(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ū)動(dòng)因素的響應(yīng)

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

表3 各階段山東省土地利用類型轉(zhuǎn)換像元統(tǒng)計(jì)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個(gè)生態(tài)指標(biāo)對RSEI交互作用探測結(jié)果 為了探究不同環(huán)境因子對山東省不同時(shí)間段內(nèi)RSEI空間分布的影響,采用地理探測器分析不同階段不同因子對RSEI的解釋力(表4)。從表4—5可知,1990—2000年、2001—2010年LUCC對RSEI空間分布變化的解釋力最強(qiáng),q值分別為0.184 2,0.184 0,但與2001—2010年相比,WET 和NDVI對RSEI的影響較大。2011—2021年RSEI空間分布變化受到WET 的影響最大,其次為LUCC,NDBSI和NDVI對其影響較小。

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

從交互探測結(jié)果可知(表5),不同階段環(huán)境因子的交互作用產(chǎn)生的q值均大于單獨(dú)因子的q值,即不同環(huán)境因子的交互作用會(huì)增加對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)濟(jì)強(qiáng)省和農(nóng)業(yè)大省,其生態(tài)的好壞直接關(guān)系到人民生活水平。但隨著城鎮(zhèn)化快速發(fā)展,城市用地的快速擴(kuò)張,導(dǎo)致其他土地利用覆蓋類型發(fā)生了顯著變化(尤其是生態(tài)用地)[3-6]。雖然山東省實(shí)施了一系列生態(tài)文明建設(shè)(如生態(tài)優(yōu)先、綠色低碳發(fā)展、黃河流域生態(tài)保護(hù)與高質(zhì)量發(fā)展策略等),但山東省生態(tài)環(huán)境質(zhì)量對地形地貌、氣候變化和人類活動(dòng)的響應(yīng)過程是一個(gè)復(fù)雜的動(dòng)態(tài)過程,不同地貌條件會(huì)導(dǎo)致熱量和水分的空間差異,進(jìn)而使得生態(tài)環(huán)境質(zhì)量變化存在明顯空間差異[7]。因此,亟需厘清生態(tài)工程實(shí)施和城鎮(zhèn)化快速發(fā)展背景下,山東省生態(tài)環(huán)境質(zhì)量時(shí)空演變及其主要驅(qū)動(dòng)因子。因此,本文利用1991—2021年的Landsat數(shù)據(jù)計(jì)算RSEI指數(shù),并輔以Sen、變異系數(shù)法、Hurst法、隨機(jī)森林對山東省生態(tài)環(huán)境質(zhì)量變化特征及其變化的驅(qū)動(dòng)因子進(jìn)行分析,主要結(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ū)呈改善趨勢,泰安市、濟(jì)南市、萊蕪市、淄博市為下降趨勢;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ì)量總體為下降趨勢;莊會(huì)波等[5]基于NDVI指數(shù)研究了山東省1982—2018年植被覆蓋度變化特征,發(fā)現(xiàn)近30 a山東省NDVI均呈增加趨勢。以上研究結(jié)論與本文結(jié)論相一致,這表明山東省生態(tài)環(huán)境質(zhì)量仍在惡化,需要大力推進(jìn)生態(tài)文明建設(shè)和城市化綠色發(fā)展理念從而改善山東生態(tài)環(huán)境質(zhì)量。

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

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

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

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