樂榮武,李 巍, *,周思楊,宋南奇

1 北京師范大學(xué)環(huán)境學(xué)院環(huán)境模擬與污染控制國家重點(diǎn)聯(lián)合實(shí)驗(yàn)室,北京 100875

2 交通運(yùn)輸部水運(yùn)科學(xué)研究院,北京 100088

生態(tài)系統(tǒng)服務(wù)(ES)被廣泛定義為直接或間接促進(jìn)可持續(xù)人類福祉的生態(tài)特征、功能或過程[1—2]。生態(tài)系統(tǒng)服務(wù)價(jià)值(ESV)是人類直接或間接地從生態(tài)系統(tǒng)功能中獲得的有形或無形的利益的客觀體現(xiàn), 是區(qū)域生態(tài)文明與可持續(xù)發(fā)展的重要表征[3]。ESV 的空間格局和演變受到自然和社會(huì)經(jīng)濟(jì)等多方面驅(qū)動(dòng)因素的影響[4],而且各驅(qū)動(dòng)因素通過直接或間接途徑產(chǎn)生復(fù)雜的相互作用關(guān)系。因此,識(shí)別ESV多個(gè)驅(qū)動(dòng)因素的交互效應(yīng),對(duì)于區(qū)域人類活動(dòng)協(xié)同管控和系統(tǒng)提升區(qū)域生態(tài)系統(tǒng)服務(wù)功能具有重要意義[5]。

當(dāng)前,對(duì)ESV的研究逐漸從ESV的評(píng)價(jià)和分析其時(shí)空變化轉(zhuǎn)向到分析驅(qū)動(dòng)機(jī)制。要摸清ESV的驅(qū)動(dòng)機(jī)制,必須識(shí)別ESV的驅(qū)動(dòng)因素并探究其交互作用關(guān)系[6—7]。探究ESV驅(qū)動(dòng)因素常用的方法包括冗余分析[8—9]、相關(guān)性分析[10]、主成分分析[11—12]、逐步回歸[13]和Logistic回歸[14]等傳統(tǒng)統(tǒng)計(jì)分析方法,這些方法僅能識(shí)別驅(qū)動(dòng)因素的作用大小和方向,無法識(shí)別多個(gè)驅(qū)動(dòng)因素之間的交互作用。有較多的研究使用地理探測(cè)器模型探究ESV的主導(dǎo)驅(qū)動(dòng)因素及各驅(qū)動(dòng)因素的交互作用[15—18],雖然地理探測(cè)器可以對(duì)ESV驅(qū)動(dòng)因素進(jìn)行非線性歸因且能識(shí)別驅(qū)動(dòng)因素之間的交互作用,但不能判斷驅(qū)動(dòng)因素的作用方向,也無法刻畫ESV對(duì)各驅(qū)動(dòng)因素的變化響應(yīng)特征[19—20]。地理探測(cè)器對(duì)驅(qū)動(dòng)因素的數(shù)據(jù)要求轉(zhuǎn)化為類型數(shù)據(jù),該處理方法涉及可變面元問題增加了研究結(jié)果的不確定性[21—22]。在交互路徑識(shí)別方面,已有學(xué)者使用結(jié)構(gòu)方程模型對(duì)驅(qū)動(dòng)因素間的作用路徑進(jìn)行定量測(cè)度[23—27],但這類研究仍然較少,有待進(jìn)一步探究?？傊?當(dāng)前的研究方法尚未系統(tǒng)全面分析ESV驅(qū)動(dòng)因素的重要性、作用方向、交互作用和刻畫ESV對(duì)驅(qū)動(dòng)因素的變化響應(yīng)特征。如何系統(tǒng)性揭示ESV變化的驅(qū)動(dòng)機(jī)制,是當(dāng)前亟需解決的問題。近年來,機(jī)器學(xué)習(xí)模型因其可以很好地表征驅(qū)動(dòng)因素與ESV之間的非線性關(guān)系而被大量用于ESV的驅(qū)動(dòng)分析[28—30]。然而,機(jī)器學(xué)習(xí)模型的大多數(shù)算法都被視為黑箱模型,并不能充分解釋驅(qū)動(dòng)作用機(jī)制。沙普利加性解釋(SHAP)方法是近期開發(fā)的機(jī)器學(xué)習(xí)模型解釋工具,不僅可以衡量驅(qū)動(dòng)因素的貢獻(xiàn)大小,還能識(shí)別驅(qū)動(dòng)因素的影響方向及交互效應(yīng)[31—32],但將其用于ESV驅(qū)動(dòng)作用機(jī)制的研究尚未見到。因此,將SHAP方法用于探索ESV驅(qū)動(dòng)因素的交互作用機(jī)制方面具有較大的應(yīng)用潛力。

呼包鄂榆城市群作為我國典型的快速城市化地區(qū),是黃河流域人口、生產(chǎn)力布局的主要載體。城市群地處農(nóng)牧交錯(cuò)帶,生態(tài)本底脆弱,資源開采、城市用地?cái)U(kuò)張等導(dǎo)致生態(tài)用地破碎化,生態(tài)系統(tǒng)極易受到破壞。本研究重點(diǎn)關(guān)注 ESV 變化與土地利用類型轉(zhuǎn)移模式之間的關(guān)系,探索ESV驅(qū)動(dòng)因素交互作用的潛在機(jī)制。具體研究內(nèi)容為:(1)分析研究期間土地利用和ESV的動(dòng)態(tài)演變;(2) 耦合隨機(jī)森林模型與SHAP 方法分析ESV驅(qū)動(dòng)因素的重要程度、變化響應(yīng)特征及其交互效應(yīng);(3)采用偏最小二乘路徑模型識(shí)別驅(qū)動(dòng)因素的交互路徑。本研究系統(tǒng)性探究了呼包鄂榆城市群ESV驅(qū)動(dòng)因素的交互作用機(jī)制,以期為黃河流域生態(tài)保護(hù)和高質(zhì)量發(fā)展提供科學(xué)參考。

1 研究區(qū)概況與數(shù)據(jù)來源

1.1 研究區(qū)概況

呼包鄂榆城市群是以內(nèi)蒙古自治區(qū)呼和浩特市、包頭市、鄂爾多斯市和陜西省榆林市組成的國家級(jí)城市群[33]。該地區(qū)國土面積約17.5萬km2,位于我國中西部的干旱半干旱區(qū),土地利用以草地為主,約占總面積的54%(圖1)。2020年總?cè)丝诩s1197萬人,平均海拔約1300 m,年均氣溫約8℃,年均降水約320 mm。該區(qū)域是我國重要的能源、煤化工基地、農(nóng)牧產(chǎn)品加工業(yè)基地和稀土新材料產(chǎn)業(yè)基地,也是黃河流域經(jīng)濟(jì)發(fā)展的重要增長極。

圖1 研究區(qū)土地利用類型分布

1.2 數(shù)據(jù)來源

研究使用的主要數(shù)據(jù)和來源包括:①土地利用類型數(shù)據(jù)、行政邊界矢量數(shù)據(jù)、氣溫、降水、潛在蒸散發(fā)、國內(nèi)生產(chǎn)總值(GDP)格網(wǎng)和人口密度格網(wǎng)數(shù)據(jù)來源于中國科學(xué)院資源與環(huán)境科學(xué)數(shù)據(jù)中心(https://www.resdc.cn/),空間分辨率均為1 km,將土地利用類型重分類分為6個(gè)一級(jí)類, 分別為耕地、林地、草地、水域、建設(shè)用地、未利用地;②高程數(shù)據(jù)從地理空間數(shù)據(jù)云(https://www.gscloud/)獲取,空間分辨率為90 m,在ArcGIS中使用分析工具計(jì)算坡度;③植被凈初級(jí)生產(chǎn)力(NPP)數(shù)據(jù)從美國國家航空航天局(https://www.nasa.gov/)獲取, 空間分辨率為500 m;④土壤屬性數(shù)據(jù)從國家地球系統(tǒng)科學(xué)數(shù)據(jù)共享服務(wù)平臺(tái)(http://www.geodata.cn/)獲取, 空間分辨率為1 km,土壤保持量數(shù)據(jù)通過通用土壤流失方程(ULSE)進(jìn)行計(jì)算;⑤其他社會(huì)經(jīng)濟(jì)數(shù)據(jù)來自《全國農(nóng)產(chǎn)品成本收益資料匯編》《呼和浩特市統(tǒng)計(jì)年鑒》《包頭市統(tǒng)計(jì)年鑒》《鄂爾多斯市統(tǒng)計(jì)年鑒》《榆林市統(tǒng)計(jì)年鑒》等。將所有空間數(shù)據(jù)統(tǒng)一到同一坐標(biāo)系,并重采樣至1 km。本研究的ESV時(shí)空分析部分均在1 km的柵格尺度上進(jìn)行,綜合權(quán)衡數(shù)據(jù)分析精度和模型計(jì)算效率,以5 km×5 km 的格網(wǎng)為基本單元,將研究區(qū)劃分為6817個(gè)格網(wǎng),在每個(gè)格網(wǎng)內(nèi)統(tǒng)計(jì)所有變量的均值,用于本研究的驅(qū)動(dòng)因素分析部分。

2 研究方法

2.1 生態(tài)系統(tǒng)服務(wù)價(jià)值核算

基于謝高地等的研究[34],1個(gè)標(biāo)準(zhǔn)單位ESV當(dāng)量因子是指1 hm2當(dāng)年平均產(chǎn)量的農(nóng)田自然糧食產(chǎn)量經(jīng)濟(jì)價(jià)值的1/7。糧食產(chǎn)量價(jià)值主要依據(jù)稻谷、小麥和玉米進(jìn)行計(jì)算,本文以研究區(qū)2015年的物價(jià)水平為基準(zhǔn),計(jì)算出研究區(qū)主要農(nóng)作物平均單產(chǎn)價(jià)值,最終得到研究區(qū)1個(gè)標(biāo)準(zhǔn)當(dāng)量經(jīng)修正后的價(jià)值量為1463.69元/hm2(表1)。

表1 不同土地利用類型單位面積價(jià)值系數(shù)/(元/hm2)

由于研究區(qū)不同年份、不同區(qū)域的生態(tài)系統(tǒng)的基本情況是變化的,ESV也相應(yīng)發(fā)生動(dòng)態(tài)變化。本文參考謝高地等的研究[34],認(rèn)為生態(tài)系統(tǒng)食物生產(chǎn)、原材料生產(chǎn)、氣體調(diào)節(jié)、氣候調(diào)節(jié)、凈化環(huán)境、維持養(yǎng)分循環(huán)、生物多樣性和美學(xué)景觀功能與生物量在總體上呈正相關(guān),水資源供給和水文調(diào)節(jié)與降水變化相關(guān),而土壤保持與降水、地形坡度、土壤性質(zhì)和植被蓋度密切相關(guān)。故選取NPP、降水量和土壤保持量三項(xiàng)因子對(duì)當(dāng)量進(jìn)行動(dòng)態(tài)調(diào)節(jié),以此構(gòu)建生態(tài)服務(wù)時(shí)空動(dòng)態(tài)價(jià)值當(dāng)量表,計(jì)算公式為:

(1)

式中,Fnij指某種生態(tài)系統(tǒng)在第i年第j地區(qū)第n類生態(tài)服務(wù)功能的單位面積價(jià)值當(dāng)量因子;Fn指該類生態(tài)系統(tǒng)的第n種生態(tài)服務(wù)價(jià)值當(dāng)量因子;n1表示與NPP相關(guān)的服務(wù)功能;n2表示與降水相關(guān)的服務(wù)功能;n3指土壤保持服務(wù)功能,Pij指NPP時(shí)空調(diào)節(jié)系數(shù),Rij指降水時(shí)空調(diào)節(jié)系數(shù),Sij指土壤保持時(shí)空調(diào)節(jié)系數(shù)。

研究區(qū)ESV計(jì)算公式為:

(2)

(3)

式中,c為第c種生態(tài)系統(tǒng)服務(wù)功能;Ec為第c種生態(tài)系統(tǒng)服務(wù)功能價(jià)值;Fnij表示某種生態(tài)系統(tǒng)在第i年第j地區(qū)第n類生態(tài)服務(wù)功能的單位面積價(jià)值當(dāng)量因子;D為多年的1個(gè)標(biāo)準(zhǔn)當(dāng)量因子的生態(tài)系統(tǒng)服務(wù)平均價(jià)值量(元/hm2),此處為1463.69元/hm2;Aij為第i年第j地區(qū)的面積。

2.2 生態(tài)系統(tǒng)服務(wù)價(jià)值變化趨勢(shì)分析

運(yùn)用一元線性回歸方程擬合1990—2020 年ESV的時(shí)空演變特征,IESV表征每個(gè)柵格ESV線性變化斜率值,若IESV>0,則表明ESV有改善趨勢(shì); 反之則表示存在退化趨勢(shì)。采用F檢驗(yàn)ESV變化趨勢(shì)的顯著性。結(jié)合IESV及F檢驗(yàn)結(jié)果,ESV變化趨勢(shì)可分為5類:大幅改善(IESV>0,P<0.01),一般改善(IESV>0,P<0.05),變化不顯著(P≥0.05),一般退化(IESV<0,P<0.05)和大幅退化(IESV<0,P<0.01)[35]。

2.3 隨機(jī)森林模型與SHAP方法

本研究使用機(jī)器學(xué)習(xí)中的隨機(jī)森林模型分析ESV的驅(qū)動(dòng)因素[36],耦合SHAP方法對(duì)隨機(jī)森林模型進(jìn)行解釋。隨機(jī)森林相對(duì)于傳統(tǒng)統(tǒng)計(jì)模型在進(jìn)行預(yù)測(cè)時(shí)往往有更好的精度,但是同時(shí)也失去了統(tǒng)計(jì)模型的可解釋性,所以隨機(jī)森林通常被認(rèn)為是黑箱模型。針對(duì)基于樹的機(jī)器學(xué)習(xí)模型,2017年Lundberg和Lee提出了SHAP方法用來解釋各種黑箱模型[37]。與以往研究中多采用線性擬合相關(guān)的方法或地理探測(cè)器模型量化ESV的驅(qū)動(dòng)因素且無法排除驅(qū)動(dòng)因素間的相互干擾相比,耦合隨機(jī)森林模型與SHAP 方法不僅可以識(shí)別ESV和驅(qū)動(dòng)因素的非線性關(guān)系,還能夠分離出每個(gè)驅(qū)動(dòng)因素對(duì)ESV的獨(dú)立影響以及各個(gè)因素間的交互影響。

將GDP密度、人口密度、降水、氣溫、潛在蒸散發(fā)、高程、坡度、土壤容重、黏土含量9種空間變量(圖2)以及6種土地類型面積比例共15個(gè)驅(qū)動(dòng)因素作為隨機(jī)森林模型的輸入變量,將ESV作為模型的輸出變量,模型樣本量為6817個(gè)。

圖2 驅(qū)動(dòng)因素的空間分布

2.4 偏最小二乘路徑模型

偏最小二乘路徑模型(PLSPM)是一種分析多變量間復(fù)雜因果關(guān)系的綜合分析模型,屬于結(jié)構(gòu)方程模型的一種方法。該模型不僅能解決變量間存在的共線性問題,對(duì)變量的分布狀態(tài)無要求,還可以計(jì)算不同變量對(duì)響應(yīng)變量的直接效應(yīng)和間接效應(yīng)[38]。模型中路徑系數(shù)反映了變量之間關(guān)系的方向和強(qiáng)度, 而預(yù)測(cè)變量和響變量之間相乘的路徑系數(shù)則顯示出間接影響的強(qiáng)度。擬合度指數(shù)(GOF)用于評(píng)估模型的預(yù)測(cè)性能, 數(shù)值越大, 模型預(yù)測(cè)效果越好。本研究使用偏最小二乘路徑模型識(shí)別氣象、地形、土壤、土地利用和社會(huì)經(jīng)濟(jì)因素對(duì)ESV的交互影響作用路徑。

3 結(jié)果與分析

3.1 土地利用類型演變

在1990—2020年間,耕地面積整體呈現(xiàn)下降趨勢(shì),在1995年面積達(dá)到峰值后逐漸下降,損失的耕地主要轉(zhuǎn)化為草地(33.95%)和林地(5.12%)。林地面積在1995年達(dá)到最低值后逐漸上升,林地增加的主要原因是來自草地(35.03%)和耕地(19.21%)的轉(zhuǎn)化。草地面積變化劇烈,草地與耕地和未利用地存在較大比例地互相轉(zhuǎn)化;水域呈現(xiàn)先減后增的趨勢(shì);建設(shè)用地持續(xù)增加,主要來自于草地(40.11%)和耕地(28.89%)的轉(zhuǎn)化(圖3)。

圖3 1990—2020年研究區(qū)土地利用轉(zhuǎn)移變化

3.2 ESV的時(shí)空演變

3.2.1ESV的時(shí)間變化

1990—2020年研究區(qū)ESV呈先下降后上升趨勢(shì),其中在1990—1995年ESV下降的最多,在此期間林地面積和水域面積減少幅度均為最大,分別減少了7.64%、14.14%,損失的林地主要轉(zhuǎn)變成耕地和草地,損失的水域主要轉(zhuǎn)變成草地。在2000—2005年ESV上升的最多,在此期間林地面積增幅最大,達(dá)到了11.89%,增長的林地主要來源于耕地和草地。ESV在研究期內(nèi)總體上增加了62.28億元,主要在于由草地和耕地轉(zhuǎn)化而來的林地面積增加了21.18%。(表2和圖4)。

表2 1990—2020年研究區(qū)各地類ESV/億元

3.2.2ESV的空間變化

在空間格局上,研究區(qū)ESV空間分異明顯(圖5)。研究區(qū)ESV總體呈現(xiàn)東高西低的分布格局,ESV高值分布以河流水系為中心,ESV較高值集中分布在呼和浩特市大青山保護(hù)區(qū)的林地,低值以未利用地類型為主要分布區(qū)。該分布趨勢(shì)與土地利用分布基本吻合,研究區(qū)西北部干旱少雨,人類活動(dòng)較少,土地利用以未利用地為主,植被覆蓋率低,ESV較低。而東南部地勢(shì)相對(duì)平坦,海拔較低降水較多,植被覆蓋率高,ESV較高;河流成為高值中心主要與水域的單位面積ESV高有關(guān)。

圖5 1990—2020年研究區(qū)ESV密度分布

總體來看,在1990—2020年間73.52%研究區(qū)的ESV變化不顯著,ESV的退化區(qū)面積占比(9.76%)高于改善區(qū)(6.38%)(圖6)。大幅退化和一般退化面積占比分別為1.62%、8.14%,大幅改善和一般改善面積占比分別為1.47%、4.91%。退化區(qū)主要集中分布在包頭市北部和零星分布在鄂爾多斯市大部分地區(qū),改善區(qū)主要分布在呼和浩特市和榆林市的東南部地區(qū)。雖然ESV的退化面積大于改善面積,但由于改善區(qū)主要位于土地類型由未利用地、草地向林地、水域轉(zhuǎn)化的區(qū)域,單位面積林地和水域的ESV更高,因此整體上研究區(qū)的ESV呈增加趨勢(shì)。

本文進(jìn)一步通過熱點(diǎn)分析揭示ESV的空間集聚特征及其演變規(guī)律。研究區(qū)ESV 冷熱點(diǎn)總體上呈“西冷東熱”的空間分布格局(圖7),研究區(qū)約84%地區(qū)的ESV空間集聚特征不顯著,ESV高值集聚和低值集聚所占區(qū)域面積比例分別約為10%、6%。ESV 熱點(diǎn)區(qū)集中分布在研究區(qū)東北部呼和浩特市北部的武川縣境內(nèi),該區(qū)域林地分布廣泛,ESV 次熱點(diǎn)區(qū)圍繞熱點(diǎn)區(qū)四周分布。 ESV 冷點(diǎn)區(qū)主要分布在研究區(qū)西北部的鄂爾多斯市杭錦旗境內(nèi),ESV 次冷點(diǎn)區(qū)聚集分布在包頭市的最北部境內(nèi),該區(qū)域分布較廣的是未利用地。從各時(shí)期階段來看,熱點(diǎn)區(qū)和冷點(diǎn)區(qū)面積呈波動(dòng)下降趨勢(shì),次熱點(diǎn)區(qū)和次冷點(diǎn)區(qū)面積呈波動(dòng)上升趨勢(shì),表明研究區(qū)ESV 高值和低值聚類均逐漸弱化。

圖7 1990—2020年研究區(qū)ESV冷熱點(diǎn)的空間分布

3.3 ESV驅(qū)動(dòng)因素分析

3.3.1驅(qū)動(dòng)因素的重要性識(shí)別

耦合隨機(jī)森林模型和 SHAP方法分析研究區(qū)ESV的驅(qū)動(dòng)因素作用方向及各驅(qū)動(dòng)因素的相對(duì)重要性(圖8)。ESV的驅(qū)動(dòng)因素重要性從大到小依次為水域比例、未利用地比例、林地比例、降水、坡度、草地比例、土壤容重、人口密度、潛在蒸散發(fā)、GDP密度、氣溫、黏土含量、高程、建設(shè)用地比例和耕地比例,表明土地類型中水域比例是影響ESV最重要的驅(qū)動(dòng)因素,其次是未利用地比例,耕地比例在研究區(qū)的重要性最低。從驅(qū)動(dòng)因素的作用方向來看,水域比例、林地比例、降水、坡度、草地、土壤容重、黏土含量和高程對(duì)研究區(qū)ESV有促進(jìn)作用,未利用地比例、人口密度、潛在蒸散發(fā)、GDP密度、氣溫、建設(shè)用地比例和耕地比例對(duì)研究區(qū)ESV有抑制作用。從驅(qū)動(dòng)因素的維度來看,土地利用類型的貢獻(xiàn)度為61.24%,成為最重要的驅(qū)動(dòng)力;其次是地形和氣象的貢獻(xiàn)度分別為17.59%和17.05%,土壤和社會(huì)經(jīng)濟(jì)的貢獻(xiàn)度分別為2.39%和1.73%。

圖8 驅(qū)動(dòng)因素重要性排序圖和模型SHAP摘要圖

3.3.2ESV對(duì)驅(qū)動(dòng)因素的響應(yīng)特征

通常,ESV驅(qū)動(dòng)因素間存在不同強(qiáng)度的相關(guān)性,這會(huì)干擾單一驅(qū)動(dòng)因素對(duì)ESV的影響分析,SHAP 方法可以較好地排除其他因素的干擾,剝離出ESV隨單一驅(qū)動(dòng)因素的變化趨勢(shì)(圖9)。

圖9 ESV對(duì)驅(qū)動(dòng)因素的變化響應(yīng)特征

總體來看,ESV對(duì)驅(qū)動(dòng)因素的響應(yīng)呈現(xiàn)出非線性變化特征。在社會(huì)經(jīng)濟(jì)因素中,當(dāng)GDP密度和人口密度較低時(shí),GDP和人口的SHAP值大于0,顯示正貢獻(xiàn),表現(xiàn)出對(duì)ESV的促進(jìn)作用;隨著經(jīng)濟(jì)的發(fā)展和人口的增長開始對(duì)ESV產(chǎn)生抑制作用。在氣象因素中,降水和氣溫對(duì)ESV的影響均存在明顯的閾值效應(yīng),降水在低于和高于250 mm時(shí)分別對(duì)ESV起抑制和促進(jìn)作用,氣溫在低于和高于7.5℃時(shí)分別對(duì)ESV起促進(jìn)和抑制作用;在潛在蒸散發(fā)低于1000 mm時(shí),潛在蒸散發(fā)對(duì)ESV的促進(jìn)作用隨潛在蒸散發(fā)的增加持續(xù)減弱,當(dāng)潛在蒸散發(fā)高于1000 mm時(shí),響應(yīng)曲線在SHAP為0值的附近維持較為平穩(wěn),表明對(duì)ESV的影響較小。高程在低于和高于1300 m時(shí)分別對(duì)ESV起抑制和促進(jìn)作用,隨著高程的增加促進(jìn)作用不再增強(qiáng),坡度對(duì)ESV始終起促進(jìn)作用。土壤容重和黏土含量在低值時(shí)對(duì)ESV無顯著作用,在高值時(shí)具有促進(jìn)作用。耕地比例較低時(shí)對(duì)ESV無顯著影響,較高時(shí)開始呈現(xiàn)抑制作用。未利用地比例在低于和高于25%時(shí)對(duì)ESV分別起促進(jìn)和抑制作用,建設(shè)用地比例對(duì)ESV始終起抑制作用。林地和水域比例始終對(duì)ESV起促進(jìn)作用,草地比例在低于和高于50%時(shí)分別起抑制和促進(jìn)作用。

3.3.3驅(qū)動(dòng)因素的交互效應(yīng)

ESV的變化受驅(qū)動(dòng)因素的綜合影響,不同驅(qū)動(dòng)因素之間還存在相互影響,SHAP 方法可捕捉交互性最強(qiáng)的成對(duì)的驅(qū)動(dòng)因素對(duì)ESV的相互作用效果(圖10)。

圖10 驅(qū)動(dòng)因素的交互效應(yīng)

GDP密度和人口密度與其他驅(qū)動(dòng)因素?zé)o明顯交互作用,這可能是由于研究區(qū)的社會(huì)經(jīng)濟(jì)因素對(duì)ESV影響的貢獻(xiàn)率較低。在以水域比例為交互項(xiàng)中,在水域比例較高的地區(qū),降水與水域比例的交互作用對(duì)ESV呈現(xiàn)出由抑制作用減弱到促進(jìn)作用增強(qiáng)的變化趨勢(shì),而氣溫、潛在蒸散發(fā)與水域比例的交互作用對(duì)ESV均呈現(xiàn)出由促進(jìn)作用減弱到抑制作用增強(qiáng)的變化趨勢(shì)。在水域比例較低的地區(qū),降水與水域比例的交互作用對(duì)ESV呈現(xiàn)出由抑制作用減弱到促進(jìn)作用增強(qiáng)的變化趨勢(shì),但交互作用程度相對(duì)較低;氣溫、潛在蒸散發(fā)與水域比例的交互作用均不明顯,未利用地比例對(duì)ESV具有更強(qiáng)的抑制作用。在降水量較多的地區(qū),高程與降水無明顯交互作用;但在降水量較少的地區(qū),高程與降水的交互作用對(duì)ESV呈現(xiàn)出由抑制作用減弱到促進(jìn)作用增強(qiáng)后又趨于穩(wěn)定的變化趨勢(shì)。在以林地比例為交互項(xiàng)中,在林地比例較高的地區(qū),坡度、土壤容重和黏土含量的增加均對(duì)ESV呈現(xiàn)抑制減弱到促進(jìn)增強(qiáng)的變化趨勢(shì),耕地比例的增加對(duì)ESV具有抑制作用;在林地比例較低的地區(qū),坡度、土壤容重、黏土含量和耕地比例與林地比例無明顯交互作用。林地比例在坡度大的地區(qū)較在坡度小的地區(qū)對(duì)ESV有更強(qiáng)的促進(jìn)作用。在以耕地比例為交互項(xiàng)中,在耕地比例較大(較小)的地區(qū),草地比例的增加對(duì)ESV的抑制(促進(jìn))作用減弱(增強(qiáng)),建設(shè)用地比例在耕地比例較大的地區(qū)對(duì)ESV具有更強(qiáng)的抑制作用。水域比例在降水較多的地區(qū)對(duì)ESV具有更強(qiáng)的促進(jìn)作用。

隨機(jī)森林模型耦合SHAP方法只能衡量兩個(gè)驅(qū)動(dòng)因素間的交互作用大小、方向及變化趨勢(shì),多種驅(qū)動(dòng)因素間的交互作用路徑如何需要進(jìn)一步探討。將驅(qū)動(dòng)因素分為社會(huì)經(jīng)濟(jì)(GDP、人口)、氣象(降水、氣溫、潛在蒸散發(fā))、地形(高程、坡度)、土壤(黏土含量、土壤容重)和土地利用類型共5個(gè)潛變量,采用偏最小二乘路徑模型探究各驅(qū)動(dòng)因素潛變量間的交互路徑(圖11),模型擬合度指數(shù)(GOF)為0.64,表明模擬結(jié)果符合精度要求。

圖11 驅(qū)動(dòng)因素的交互作用路徑和影響程度

土地利用類型對(duì)ESV的影響系數(shù)為0.696(圖11),再次表明土地利用對(duì)ESV的影響起著決定作用。社會(huì)經(jīng)濟(jì)受到地形因素的制約影響,社會(huì)經(jīng)濟(jì)對(duì)土地利用和ESV的影響系數(shù)分別為0.222和-0.186,表明社會(huì)經(jīng)濟(jì)對(duì)ESV的直接影響是負(fù)效應(yīng),社會(huì)經(jīng)濟(jì)主要通過直接影響土地利用進(jìn)而間接影響ESV。氣象對(duì)ESV的直接作用較弱,主要是通過氣象-土壤-土地利用或氣象-土地利用作用路徑影響 ESV。地形對(duì)ESV的影響路徑有多條,地形對(duì)社會(huì)經(jīng)濟(jì)、氣象、土壤和土地利用均有直接影響,其中對(duì)土壤的直接影響力最強(qiáng)。土壤受到地形、社會(huì)經(jīng)濟(jì)和氣象的影響,土壤對(duì)ESV的影響主要是通過影響土地利用產(chǎn)生的間接影響。除土地利用以外,其他因素對(duì)ESV的間接影響均大于直接影響。

4 討論

在ESV的變化響應(yīng)特征方面,以前的研究尚未有ESV對(duì)驅(qū)動(dòng)因素變化響應(yīng)特征的探討;在交互作用分析方面,當(dāng)前大量的研究采用地理探測(cè)器方法,但該方法未能識(shí)別交互作用的方向性;耦合隨機(jī)森林模型和SHAP的方法能識(shí)別驅(qū)動(dòng)因素在不同值域內(nèi)的交互作用大小、方向和非線性響應(yīng)特征;在交互路徑分析方面,當(dāng)前研究主要采用結(jié)構(gòu)方程模型識(shí)別交互路徑,本文選擇結(jié)構(gòu)方程模型中的偏最小二乘路徑模型,該方法對(duì)樣本量和樣本分布要求低[23,38],模擬結(jié)果精度較好,表現(xiàn)出了較好的適用性。綜上,本研究耦合隨機(jī)森林模型和SHAP方法并結(jié)合偏最小二乘路徑模型,可以較為系統(tǒng)全面地分析ESV的驅(qū)動(dòng)因素及其交互效應(yīng)。

本研究表明土地利用的變化強(qiáng)烈影響 ESV,這與之前的研究一致[23,39]。ESV在2000年前后分別呈下降和上升的趨勢(shì),ESV的增加主要來自于耕地和草地向林地轉(zhuǎn)換的結(jié)果。林地面積在2000年之后逐漸增加,這得益于當(dāng)?shù)亍巴烁€林還草”政策和“三北”防護(hù)林工程的實(shí)施。尤其是在ESV大幅改善的區(qū)域,例如,位于榆林北部的毛烏素沙地,由于近二十年的持續(xù)造林治沙活動(dòng),ESV得到了大幅提升。人類社會(huì)經(jīng)濟(jì)活動(dòng)對(duì)生態(tài)系統(tǒng)服務(wù)的需求是其影響ESV變化的根本原因[24],本研究中GDP和人口的增長對(duì)ESV的影響由促進(jìn)作用轉(zhuǎn)變?yōu)橐种谱饔?表明較低程度的人類活動(dòng)對(duì)生態(tài)環(huán)境具有改善作用,但在高收入地區(qū)大量人口的高需求和有限的供應(yīng)水平將導(dǎo)致生態(tài)系統(tǒng)服務(wù)稀缺,從而對(duì)ESV產(chǎn)生抑制作用。偏最小二乘路徑模型進(jìn)一步揭示了社會(huì)經(jīng)濟(jì)因素對(duì)ESV的直接和間接影響是相反的,社會(huì)經(jīng)濟(jì)通過氣象、土壤和土地利用對(duì) ESV 產(chǎn)生間接的積極影響(圖11),表明經(jīng)濟(jì)發(fā)展可能會(huì)在自然變化過程中改變生態(tài)系統(tǒng)服務(wù)演化的方向[40]。以往的研究表明,城市擴(kuò)張和經(jīng)濟(jì)發(fā)展會(huì)直接加劇生態(tài)系統(tǒng)服務(wù)供需失衡[41]。此外,更多生態(tài)環(huán)境政策的實(shí)施也會(huì)間接促進(jìn)了生態(tài)系統(tǒng)服務(wù)功能的恢復(fù)和提升[42]。因此,結(jié)合經(jīng)濟(jì)發(fā)展和自然條件的空間差異,優(yōu)化土地利用配置,對(duì)提高ESV具有重要意義。考慮到水域是最重要的驅(qū)動(dòng)因素,草地對(duì)研究區(qū)的ESV貢獻(xiàn)最大。未來的生態(tài)恢復(fù)實(shí)踐應(yīng)優(yōu)先考慮濕地和湖泊生態(tài)系統(tǒng)提供的生態(tài)價(jià)值和效益,通過跨省聯(lián)動(dòng)、多方補(bǔ)水、關(guān)閉周圍煤礦等保護(hù)措施,加強(qiáng)紅堿淖濕地國家自然保護(hù)區(qū)的生態(tài)保護(hù);在沿黃河濕地建設(shè)沿黃生態(tài)廊道,嚴(yán)格保護(hù)基本林地和草原。

本研究發(fā)現(xiàn)降水和氣溫對(duì)ESV的影響具有閾值效應(yīng),在年降水量低于250 mm或年均氣溫高于7.5℃的地區(qū),降水或氣溫成為該地區(qū)ESV提升的限制因子,即該地區(qū)降水和氣溫對(duì)ESV具有負(fù)貢獻(xiàn),隨著降水的增加或氣溫的降低這種負(fù)貢獻(xiàn)逐漸減小直至變?yōu)檎暙I(xiàn)。在整個(gè)研究區(qū)內(nèi),降水的增加或氣溫的降低始終對(duì)ESV具有提升作用,原因在于研究區(qū)位于干旱半干旱區(qū),熱量充足而降水量小,有限的降水還未被植被充分吸收利用就已經(jīng)蒸散發(fā)了。在交互分析中,水域比例是參與最多的交互項(xiàng),主要在于ESV當(dāng)量表中水域的值最大,水域比例成為最重要的驅(qū)動(dòng)因素,這也導(dǎo)致降水和氣溫在水域比例較高的地區(qū)對(duì)ESV具有更大的影響程度。

本研究仍存在一定的局限性。首先,基于當(dāng)量因子法的ESV的評(píng)價(jià)嚴(yán)重依賴于土地利用數(shù)據(jù),在一定程度上可能夸大了土地利用類型對(duì)ESV 的影響,后期應(yīng)嘗試采用其他ESV 評(píng)價(jià)方法進(jìn)一步驗(yàn)證不同驅(qū)動(dòng)因素的重要性。其次,本研究ESV的驅(qū)動(dòng)因素分析是基于分辨率為5 km的格網(wǎng)尺度,但是不同空間尺度上的驅(qū)動(dòng)影響和交互關(guān)系尚不清楚,尺度效應(yīng)增加了研究結(jié)論的不確定性,后續(xù)應(yīng)該在不同空間尺度探究ESV的驅(qū)動(dòng)因素。此外,本研究尚未考慮政策因素對(duì)ESV的影響,因此未來的工作應(yīng)該將政策因素充分考慮進(jìn)來。

5 結(jié)論

本研究核算了1990—2020年呼包鄂榆城市群的ESV并分析其時(shí)空變化特征,耦合隨機(jī)森林模型和SHAP方法探究了ESV驅(qū)動(dòng)因素的貢獻(xiàn)程度、作用方向、變化響應(yīng)特征及其交互效應(yīng),進(jìn)一步采用偏最小二乘路徑模型識(shí)別了交互路徑。結(jié)論如下:(1)草地對(duì)研究區(qū)的ESV貢獻(xiàn)最大,其次是水域;研究區(qū)的ESV在1990—2020年呈波動(dòng)上升趨勢(shì),總體增加了62.28億元,ESV的增長主要源自于單位面積ESV較低的草地和耕地向單位面積ESV較高的林地的轉(zhuǎn)化。(2)ESV對(duì)驅(qū)動(dòng)因素呈非線性變化響應(yīng)特征;土地利用類型是最重要的驅(qū)動(dòng)因素,其貢獻(xiàn)度達(dá)到61.24%,地形和氣象的貢獻(xiàn)度分別為17.59%和17.05%,土壤和社會(huì)經(jīng)濟(jì)的貢獻(xiàn)度較低,分別為2.39%和1.73%。(3)水域比例是最重要的交互項(xiàng),其次是林地比例;不同因素間交互作用在因素處于不同范圍內(nèi)表現(xiàn)出不同的交互效應(yīng)。(4)土地利用類型直接作用于ESV,其他因素主要通過直接影響土地利用類型進(jìn)而間接影響ESV,且對(duì)ESV的間接影響均大于直接影響。