馮嫻慧 ，曾芝琳

1.華南理工大學(xué)建筑學(xué)院，廣東 廣州 510641；2.亞熱帶建筑科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，廣東 廣州 510641；3.長江大學(xué)醫(yī)學(xué)部，湖北 荊州 434023

植被是覆蓋地表的植物群落的總稱，是陸地生態(tài)系統(tǒng)中重要的組成部分，承擔(dān)著物質(zhì)循環(huán)、能量轉(zhuǎn)換等重要功能，在保持水土、調(diào)節(jié)氣候以及實(shí)現(xiàn)全球碳平衡等方面都具有不可替代的作用（Piao et al.，2011；Peng et al.，2012）。植被覆蓋度（Fractional Vegetation Cover，F(xiàn)VC）是衡量地表植被覆蓋狀況的一個(gè)綜合量化指標(biāo)，指單位面積內(nèi)植被地上部分（包括葉、莖、枝）在地面的垂直投影面積占統(tǒng)計(jì)總面積的百分比（Gitelson et al.，2002）。植被覆蓋度對環(huán)境變化較敏感，能夠反映區(qū)域范圍的生態(tài)環(huán)境特征，是區(qū)域乃至全球生態(tài)、水文、土壤等學(xué)科研究領(lǐng)域的重要參數(shù)和基礎(chǔ)數(shù)據(jù)（賈坤等，2013）。獲取植被覆蓋度的途徑主要可分為地表實(shí)測和遙感監(jiān)測兩類。由于遙感監(jiān)測具有高效、精確的顯著優(yōu)勢，已成為獲取區(qū)域植被覆蓋度的重要手段（程紅芳等，2008）。植被指數(shù)（Vegetation Index，VI）是通過遙感方法表征植被覆蓋度的指標(biāo)，常用植被指數(shù)包括歸一化植被指數(shù)（Normalized Difference Vegetation Index，NDVI）、增強(qiáng)植被指數(shù)（Enhanced Vegetation Index，EVI）、比值植被指數(shù)（Ratio Vegetation Index，RVI）、差值植被指數(shù)（Difference Vegetation Index，DVI）等（Xiao et al.，2005）。

國內(nèi)外學(xué)者利用遙感方法，在不同規(guī)模尺度上采用植被指數(shù)進(jìn)行植被覆蓋研究。在全球尺度上，Myneni et al.（1997）、Tucker et al.（2001）使用NDVI數(shù)據(jù)分析北半球1981—1991年和1982—1999年的植被覆蓋動(dòng)態(tài)變化；Fensholt et al.（2012）分析全球半干旱地區(qū) 1981—2007年的植被覆蓋變化和驅(qū)動(dòng)因素；Piao et al.（2011）研究歐亞大陸1982—2006年間的植被覆蓋變化。在區(qū)域尺度，甘春英等（2011）、穆少杰等（2012）、劉憲鋒等（2013）分別使用Landsat TM影像、MODIS-NDVI數(shù)據(jù)研究連江流域、內(nèi)蒙古、黃土高原地區(qū)等地10余年間的植被覆蓋變化。在城市群尺度，董晨煒等（2017）基于 MODIS NDVI數(shù)據(jù)分析環(huán)杭州灣地區(qū) 2000—2013年的植被變化；賀忠華等（2020）使用EVI指數(shù)分析浙江省20年植被變化；李卓等（2017）基于NDVI數(shù)據(jù)分析京津冀城市群 2005—2015年的植被覆蓋度變化；陳丹等（2006）利用NDVI數(shù)據(jù)分析廣東省2000—2005年的植被覆蓋狀況。近3年，鄧玉嬌等（2021）、王思等（2022）基于NDVI數(shù)據(jù)分析廣東省2000—2018和2000—2020年以來的植被覆蓋變化。何全軍（2019）、Hu et al.（2019）利用NDVI分別研究了珠三角 2001—2017年和 2000—2016年的植被覆蓋變化以及氣候響應(yīng)。馮嫻慧等（2022）、趙桔超等（2022）分別使用NDVI指數(shù)分析粵港澳大灣區(qū)近20年的植被覆蓋變化特征。NDVI是最常用的植被覆蓋度指數(shù)，但粵港澳灣區(qū)屬于南亞熱帶—熱帶季風(fēng)氣候區(qū)，植被生長茂盛，NDVI易達(dá)到飽和，增強(qiáng)型植被指數(shù)（EVI）在減少大氣作用以及飽和問題上優(yōu)于 NDVI；在相同的空間分辨率下，EVI能更好地反映植被覆蓋的空間差異。因此，除NDVI指數(shù)研究外，有必要進(jìn)一步采用EVI指數(shù)進(jìn)行粵港澳灣區(qū)的植被覆蓋變化研究，目前未見采用EVI在大灣區(qū)尺度進(jìn)行植被覆蓋的相關(guān)研究。

在植被覆蓋變化的驅(qū)動(dòng)影響因子研究中，以氣候類因子中的氣溫和降水因子的影響研究為主。陳超男等（2019）、張凱選等（2019）分別采用相關(guān)系數(shù)分析秦巴山區(qū)、西南喀斯特地區(qū)植被覆蓋變化與氣溫、降水的關(guān)系；穆少杰等（2012）采用偏相關(guān)系數(shù)探究內(nèi)蒙古植被覆蓋度對氣溫和降水的響應(yīng)；晏利斌等（2011）、徐勇等（2020）分別采用奇異值分解法（SVD）和相關(guān)系數(shù)法分析京津冀地區(qū)植被覆蓋變化與氣溫、降水的關(guān)系；陳云浩等（2001）通過偏相關(guān)和復(fù)相關(guān)系數(shù)建立氣候因子驅(qū)動(dòng)的植被分區(qū)。何全軍（2019）、賀忠華等（2020）、鄧玉嬌等（2021）采用相關(guān)系數(shù)分別分析了浙江省、廣東省、珠三角地區(qū)植被覆蓋變化對氣溫、降水、日照時(shí)數(shù)的響應(yīng)變化。影響植被的自然因素比較復(fù)雜，雖已有大量關(guān)于氣候因素對植被影響研究，但更多的自然因子，尤其是非氣候類因子，如土壤類型、高程、坡度等對植被覆蓋變化的影響以及驅(qū)動(dòng)力貢獻(xiàn)的研究亟待開展，但非氣候類自然因子對植被覆蓋變化影響的相關(guān)研究很少。

綜上所述，基于空間分辨率為 250 m的MOD13Q1-EVI數(shù)據(jù)，采用EVI指數(shù)，使用Sen趨勢分析和Mann-Kendall模型檢驗(yàn)，Pearson相關(guān)系數(shù)分析和t檢驗(yàn)，地理探測器驅(qū)動(dòng)力分析模型等方法研究粵港澳大灣區(qū)20年期間（2001—2020）植被覆蓋變化以及氣候類因子和非氣候類因子等共9個(gè)自然因子對植被變化的影響及其驅(qū)動(dòng)力貢獻(xiàn)，研究結(jié)果預(yù)期對監(jiān)測灣區(qū)植被動(dòng)態(tài)，了解自然要素對植被的影響，以及制定灣區(qū)生態(tài)系統(tǒng)保護(hù)的統(tǒng)籌建設(shè)和優(yōu)化管理具有重要意義。

1 數(shù)據(jù)與方法

1.1 研究區(qū)概況

本研究區(qū)域?yàn)榛浉郯拇鬄硡^(qū)（The Guangdong-Hong Kong-Macao Greater Bay Area，GBA），位于廣東省中部（21°26′—24°28′N，111°14′—115°24′E），總面積5.6×104km2，包括珠三角九市（廣州、深圳、珠海、佛山、惠州、東莞、中山、江門、肇慶）和香港、澳門兩個(gè)特別行政區(qū)?；浉郯拇鬄硡^(qū)是中國重點(diǎn)發(fā)展建設(shè)的世界級城市群及都市連綿區(qū)，是中國經(jīng)濟(jì)增長的重要引擎和參與全球競爭的戰(zhàn)略性空間載體。在過去20年間，快速的經(jīng)濟(jì)增長和城市大幅擴(kuò)張使得大灣區(qū)生態(tài)環(huán)境面臨嚴(yán)峻考驗(yàn)。對此，《粵港澳大灣區(qū)發(fā)展規(guī)劃綱要》特別指出，推進(jìn)生態(tài)文明建設(shè)、加強(qiáng)生態(tài)環(huán)境保護(hù)是進(jìn)一步提升粵港澳大灣區(qū)發(fā)展質(zhì)量的重要支撐。

大灣區(qū)屬南亞熱帶—熱帶季風(fēng)氣候，年均氣溫21—23 ℃，年均降水量1300—2500 mm。大灣區(qū)東、西、北部三面山地丘陵圍繞，中部為平原，南部瀕臨南海。2020年末，常住人口7264.92萬人。大灣區(qū)主要土地覆蓋類型為林地、耕地與居住及建設(shè)用地，林地主要分布于東、西、北部山地丘陵地帶，如肇慶、江門、惠州及廣州北部，居住及建設(shè)用地主要分布于中部平原城市集群區(qū)，如廣州南部、佛山、東莞、深圳、香港、中山、珠海、澳門等地。研究區(qū)域的地帶性植被以南亞熱帶季風(fēng)常綠闊葉林、中亞熱帶典型常綠闊葉林和海岸紅樹林為主。

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

1.2.1 植被覆蓋數(shù)據(jù)

研究區(qū)域植被覆蓋生長茂密，為減少植被指數(shù)飽和對研究結(jié)果的影響，本研究采用增強(qiáng)型植被指數(shù)（EVI）數(shù)據(jù)進(jìn)行植被趨勢變化研究。數(shù)據(jù)來源于美國國家航空航天局（National Aeronautics and Space Administration，NASA）發(fā)布的 2001—2020年MOD13Q1產(chǎn)品，時(shí)間分辨率為16 d，空間分辨率為250 m，每年有23幅影像。數(shù)據(jù)記錄始于2000年2月，為保證數(shù)據(jù)序列的完整性，所選取的影像時(shí)間范圍為 2001—2020年。使用 MRT（MODIS Reprojection Tools）工具進(jìn)行投影轉(zhuǎn)換，投影坐標(biāo)系選擇Albers Equal Area，借用研究區(qū)矢量邊界進(jìn)行掩膜裁剪。運(yùn)用最大值合成法（MVC）將每年全部影像進(jìn)行合成，得到 2001—2020年大灣區(qū)逐年的年最大EVI值。采用S-G濾波法（Savitzky-Golay filter）對年最大EVI值進(jìn)行時(shí)間序列重構(gòu)，再經(jīng)計(jì)算得到大灣區(qū)EVI多年均值，并按等間距法分為以下 5 類：低（EVI＜0.2）、較低（0.2≤EVI＜0.4）、中（0.4≤EVI＜0.6）、較高（0.6≤EVI＜0.8）、高（EVI≥0.8）。

植被類型數(shù)據(jù)來源于中國科學(xué)院資源環(huán)境科學(xué)數(shù)據(jù)中心的《1∶1000000中國植被集》，分辨率為1 km，按照植被大類分為7類。土地利用數(shù)據(jù)來源于中國科學(xué)院資源環(huán)境科學(xué)數(shù)據(jù)中心所提供的2020年中國土地利用數(shù)據(jù)，分辨率為1 km，包括耕地、林地、草地、水域、居民地和未利用土地6個(gè)一級類型。以上數(shù)據(jù)經(jīng)投影變換、研究區(qū)邊界掩膜、重采樣等處理后與EVI數(shù)據(jù)投影、分辨率一致。

1.2.2 氣候自然因子數(shù)據(jù)來源與預(yù)處理

本研究使用同期（2001—2020年），時(shí)間跨度為 20年的氣象數(shù)據(jù)，數(shù)據(jù)源為中國科學(xué)院資源環(huán)境科學(xué)數(shù)據(jù)中心所提供的大灣區(qū)及周邊國家地面氣象觀測站和香港基準(zhǔn)天氣站共96個(gè)站點(diǎn)的20年逐日（日值）氣象數(shù)據(jù)，氣象站點(diǎn)位置分布參見圖1。氣象數(shù)據(jù)包括氣溫、降水、風(fēng)速、相對濕度、日照時(shí)數(shù)等共5個(gè)因子。日值數(shù)據(jù)采用線性內(nèi)插法對各站點(diǎn)數(shù)據(jù)進(jìn)行填補(bǔ)缺失值處理，并合成為站點(diǎn)年值數(shù)據(jù)。采用年值氣象數(shù)據(jù)，使用薄板樣條法和克里金法（Kriging）進(jìn)行插值。其中氣溫、降水受高程影響較大，采用薄板樣條法，以高程為協(xié)變量對氣溫、降水進(jìn)行插值；風(fēng)速、相對濕度、日照時(shí)數(shù)氣象數(shù)據(jù)使用常規(guī)的克里金法進(jìn)行插值，輸出插值結(jié)果分辨率均為250 m。

圖1 氣象站點(diǎn)分布Figure 1 Distribution of Meteorological station

1.2.3 非氣候類自然因子數(shù)據(jù)來源與預(yù)處理

非氣候類自然因子包括高程、坡度、坡向、土壤類型共4個(gè)自然因子。高程數(shù)據(jù)來源于中國科學(xué)院計(jì)算機(jī)網(wǎng)絡(luò)信息中心地理空間數(shù)據(jù)云平臺(tái)所提供的GDEMV3數(shù)據(jù)，分辨率為30 m，按自然間斷點(diǎn)法分為5類。坡度、坡向數(shù)據(jù)在ArcGIS軟件中基于GDEMV3高程數(shù)據(jù)計(jì)算得出，分辨率為30 m，坡度按自然間斷點(diǎn)法分為5類，坡向等分為9類。土壤類型數(shù)據(jù)來源于中國科學(xué)院資源環(huán)境科學(xué)數(shù)據(jù)中心提供的《1∶100萬中華人民共和國土壤圖》中的數(shù)字化柵格圖，分辨率為1 km，按照土壤發(fā)生的分類系統(tǒng)，將土類分類為18類。

1.3 研究方法

1.3.1 Sen趨勢分析和Mann-Kendall檢驗(yàn)

Theil-Sen Median方法又被稱為Sen趨勢分析，是一種穩(wěn)健的非參數(shù)統(tǒng)計(jì)的趨勢計(jì)算方法。該方法對于測量誤差和離群數(shù)據(jù)不敏感，常被用于長時(shí)間序列數(shù)據(jù)的趨勢分析，與線性回歸法要求時(shí)間序列數(shù)據(jù)符合正態(tài)分布的要求相比較，Sen趨勢分析是經(jīng)過計(jì)算序列中的中值，更適用于估計(jì)植被長期變化的趨勢值，通常與Mann-Kendall非參數(shù)檢驗(yàn)結(jié)合使用。Mann-Kendall是一種非參數(shù)統(tǒng)計(jì)檢驗(yàn)方法，其計(jì)算不需要測量值服從正態(tài)分布，也不要求趨勢是線性的，且不受缺失值和異常值的影響，在長時(shí)間序列數(shù)據(jù)的趨勢顯著檢驗(yàn)中被廣泛應(yīng)用。Sen趨勢分析法和Mann-Kendall檢驗(yàn)均為非參數(shù)分析方法。

本研究采用Theil-Sen法計(jì)算大灣區(qū)近20年來逐像元的EVI變化趨勢，并使用Mann-Kendall法進(jìn)行顯著性檢驗(yàn)。Theil-Sen斜率計(jì)算公式如下：

式中：

β——EVI的年際變化斜率，β＞0表示 EVI呈上升趨勢，β＜0表示EVI呈下降趨勢；

xj和xi——時(shí)間序列數(shù)據(jù)，分別為第i年和第j年的EVI值；

n——時(shí)間序列長度。

Man-Kendall趨勢檢驗(yàn)的統(tǒng)計(jì)量S計(jì)算公式如下：

式中：

sgn(xj-xi)——符號函數(shù)，其定義見式（3）：

檢驗(yàn)統(tǒng)計(jì)量S由公式（2）計(jì)算，當(dāng)n＜10時(shí)，如果S＞0，則認(rèn)為序列存在上升趨勢；S=0，無趨勢；S＜0認(rèn)為序列存在下降趨勢。當(dāng)n≥10時(shí)，統(tǒng)計(jì)量S大致服從正態(tài)分布，Var(S)是S統(tǒng)計(jì)量的方差，其計(jì)算見公式（4），檢驗(yàn)統(tǒng)計(jì)量Z由公式（5）計(jì)算。本研究時(shí)間序列長度為20（2001—2020年），采用檢驗(yàn)統(tǒng)計(jì)量Z進(jìn)行趨勢檢驗(yàn)。

標(biāo)準(zhǔn)化統(tǒng)計(jì)量Z的計(jì)算公式見式（5）：

給定顯著性水平α下，當(dāng)|Z|＞u_(1-α/2)時(shí)，表明所研究序列在α水平上存在顯著變化。由于Theil-Sen趨勢分析中基本不存在β=0的區(qū)域，故將-0.0005≤β≤0.0005的區(qū)域視為穩(wěn)定不變的區(qū)域，將β＞0.0005的區(qū)域劃分為改善區(qū)，β＜-0.0005的區(qū)域劃分為退化區(qū)。本研究取顯著性水平0.01和0.05進(jìn)行EVI檢驗(yàn)，將Mann-Kendall檢驗(yàn)結(jié)果依據(jù)α=0.01和α=0.05的置信水平劃分為極顯著變化（|Z|＞2.58）、顯著變化（1.96＜|Z|≤2.58）和不顯著變化（|Z|≤1.96），即表示趨勢分別通過了置信度為95%和99%的顯著性檢驗(yàn)。將β與Z值結(jié)合進(jìn)行分類，從而將研究區(qū)的EVI變化趨勢劃分為極顯著改善、顯著改善、不顯著改善、穩(wěn)定不變、極顯著退化、顯著退化、不顯著退化7種類型。

1.3.2 皮爾森相關(guān)系數(shù)和t檢驗(yàn)

本研究采用相關(guān)系數(shù)法分析植被覆蓋變化與氣候自然因子的相關(guān)性。Pearson相關(guān)系數(shù)是衡量兩個(gè)變量線性相關(guān)性強(qiáng)弱的指標(biāo)，其計(jì)算公式如下：

式中：

r——變量X與Y的相關(guān)系數(shù)，分別為EVI序列和氣候因子序列；

n——樣本數(shù)；

Xi——序列第i個(gè)數(shù)據(jù)的EVI值；

Y——?dú)夂蛞蜃又笖?shù)的第i個(gè)數(shù)據(jù)值。當(dāng)r＞0時(shí)，代表兩變量呈線性正相關(guān)；當(dāng)r＜0時(shí)，代表兩變量呈線性負(fù)相關(guān)。r的絕對值越大，說明兩變量的線性相關(guān)的程度越強(qiáng)。通過 MATLAB軟件，逐像元計(jì)算2001—2020年EVI值與氣溫、降水量、風(fēng)速、相對濕度、日照時(shí)數(shù)年均值之間的相關(guān)系數(shù)，并采用t檢驗(yàn)法進(jìn)行顯著性檢驗(yàn)。結(jié)合置信度P值將相關(guān)系數(shù)的顯著性劃分為4個(gè)等級：顯著負(fù)相關(guān)（r＜0，P＜0.05）、不顯著負(fù)相關(guān)（r＜0，P≥0.05）、不顯著正相關(guān)（r＞0，P≥0.05）、顯著正相關(guān)（r＞0，P＜0.05）。

1.3.3 驅(qū)動(dòng)力分析模型

本研究采用地理探測器模型（geographical detectors）進(jìn)行驅(qū)動(dòng)力分析，地理探測器是一種探測空間分異性并揭示其背后驅(qū)動(dòng)因子的統(tǒng)計(jì)學(xué)模型，也是進(jìn)行驅(qū)動(dòng)力分析的工具之一（王勁峰等，2017）。其理論核心為基于空間分異理論，結(jié)合GIS空間疊加技術(shù)，通過檢驗(yàn)兩個(gè)變量空間分布的一致性，來探測兩者間可能的因果關(guān)系（劉彥隨等，2017）。模型能夠定量分析影響因子對因變量的貢獻(xiàn)程度，具有無線性假設(shè)、對多自變量共線性免疫等優(yōu)勢，被應(yīng)用于分析地理要素格局演變和地域空間分異等方面，包括因子探測、交互作用探測等算法。因子探測算法主要用于探測因變量Y在空間分布上的異質(zhì)性，并量化自變量因子X對Y值的空間分異的解釋力大小，以q值來度量，其計(jì)算式如下：

式中：

σh2和σ2——層h和全區(qū)的Y值的方差，h=1,2, 3…；

q——某因子X解釋了100×q%的Y，q值值域?yàn)閇0, 1]，q值越大表示自變量X對Y的解釋力越強(qiáng)；

L——因變量Y或自變量X的分層，即分類或分區(qū)；

Nh和N——層h和全區(qū)的單元數(shù)。

交互作用探測主要用于識(shí)別不同X因子間的交互作用，即評估因子X1和X2共同作用時(shí)是否會(huì)增加或減弱對因變量Y的解釋力，兩個(gè)因子間包括非線性減弱、非線性增強(qiáng)、雙因子增強(qiáng)和獨(dú)立等交互關(guān)系分析。

2 結(jié)果與分析

2.1 粵港澳大灣區(qū)植被覆蓋特征與變化趨勢

基于 EVI逐年最大值空間分布，逐像元進(jìn)行Sen趨勢分析和Mann-Kendall顯著性檢驗(yàn)，可以有效反映近20年大灣區(qū)EVI的空間變化趨勢。將βEVI在-0.0005和0.0005之間的區(qū)域劃分為穩(wěn)定區(qū)域，高于0.0005的區(qū)域劃分為改善區(qū)域，低于-0.0005的區(qū)域劃分為退化區(qū)域。取置信度水平0.01和0.05，相應(yīng)Z值分別為2.58和1.96，通過Sen趨勢分析和Mann-Kendall顯著性檢驗(yàn)，植被變化趨勢可劃分為7個(gè)等級，大灣區(qū)植被EVI指數(shù)的均值分布特征和變化趨勢如圖2、表1所示。近20年，大灣區(qū)植被覆蓋整體呈改善趨勢，植被改善區(qū)的面積大于退化區(qū)面積，極顯著改善的植被區(qū)域占比66.98%，極顯著退化區(qū)面積占比在5.70%。EVI變化趨勢的空間分布具有明顯的差異性和集聚性（圖2b）。顯著改善的植被區(qū)域具有高度集聚性，主要分布于大灣區(qū)北部、西部、東部大部分區(qū)域；顯著退化區(qū)域分布比較零散，空間差異性顯著，其中相對集中分布區(qū)域主要是廣州南部、佛山西部、珠海南部、中山北部和東莞中部。整體上看，退化區(qū)主要分布在城鎮(zhèn)發(fā)展的周邊區(qū)域，極顯著改善區(qū)主要分布在大灣區(qū)東、西、北部等外圍區(qū)域，與中部城市密集建設(shè)區(qū)有一定的距離，整體受人類活動(dòng)的干擾程度較小。

圖2 2001—2020粵港澳大灣區(qū)年均EVI分布（a）和EVI變化趨勢（b）Figure 2 Multi-year mean EVI (a) and EVI trends (b) in The Guangdong-Hong Kong-Macao Greater Bay Area(GBA) from 2001 to 2020

表1 EVI變化趨勢統(tǒng)計(jì)Table 1 Statistic of EVI trends

據(jù)EVI年均值分布圖（圖2a）和不同土地利用類型、不同植被類型的EVI均值統(tǒng)計(jì)表（表2）分析，大灣區(qū)植被覆蓋特征表現(xiàn)為，整體覆蓋狀況良好，EVI均值在0.4—0.6之間的區(qū)域占比為37.43%，0.6—0.8之間占比40.66%。在空間分布上，植被覆蓋度呈現(xiàn)“四周高、中間低”的特征，中部珠三角城市密集區(qū)的城鄉(xiāng)建設(shè)空間植被覆蓋度低。在不同土地利用類型中，城鄉(xiāng)建設(shè)用地以及未利用土地EVI均值在0.2—0.4區(qū)間，植被類型以栽培植被為主；EVI 均值0.4—0.6區(qū)間的土地類型主要是耕地，其次是草地，仍以栽培植被為主要類型；EVI均值0.6—0.8區(qū)間的主要是林地，以闊葉林為主，其次為針葉林。

表2 土地利用類型和不同植被上EVI的統(tǒng)計(jì)Table 2 Statistic of EVI in different land use and vegetation type

將大灣區(qū)土地利用類型和植被類型空間與EVI變化趨勢空間分布疊加統(tǒng)計(jì)，得到在不同土地利用類型下的EVI變化趨勢以及不同類型植被的EVI變化趨勢（圖 3）。在不同土地利用類型上，植被變化趨勢均表現(xiàn)為極顯著改善。其中，林地、草地、未利用土地的植被均未顯示出退化和極顯著退化；城鄉(xiāng)建設(shè)用地、水域、耕地的植被變化趨勢表現(xiàn)為以改善為主，但改善和退化等7種不同程度的變化趨勢共存，城鄉(xiāng)建設(shè)用地的極顯著退化趨勢最高，占比面積不足18.59%，水域次之，占比12.79%。不同類型植被的EVI均存在顯著改善的趨勢；其中，栽培植被、高山植被的發(fā)展趨勢具有明顯的改善、退化等7種不同程度的變化共存。

圖3 2001—2020大灣區(qū)不同植被類型（a）和不同土地利用類型（b）的EVI變化趨勢Figure 3 EVI trends of vegetation types (a) and EVI trends of land use type (b) in GBA from 2001 to 2020

2.2 氣候因子對植被變化趨勢的影響分析

2.2.1 氣候因子年際變化特征

統(tǒng)計(jì)大灣區(qū)及周邊國家地面氣象觀測站和香港基準(zhǔn)天氣站共96個(gè)站點(diǎn)的20年氣象數(shù)據(jù)，得到氣溫、降水、風(fēng)速、濕度、日照等氣候因子的總體變化特征。

近20年間，氣溫、相對濕度呈上升趨勢，大灣區(qū)地處濕熱多雨的夏熱冬暖地區(qū)，從氣候變化來看，變化趨勢日益表現(xiàn)出高溫多濕的特征。氣溫呈現(xiàn)緩慢上升的趨勢，但氣溫變化幅度很小，氣溫增速約 0.046 ℃·(10 a)-1，r2=0.0078，平均氣溫22.17 ℃，最高平均氣溫出現(xiàn)在 2019年，為22.66 ℃，最低值出現(xiàn)在2011年，為21.56 ℃。相對濕度呈先下降后上升的趨勢，總體呈現(xiàn)出增加的趨勢，濕度增速約1.972 ℃·(10 a)-1，R2=0.2322，平均值為77.33%，最高值出現(xiàn)在2016年，相對濕度81.52%，最低值出現(xiàn)在2011年，為72.39%。

降水量、風(fēng)速、日照時(shí)數(shù)呈現(xiàn)出穩(wěn)定的趨勢。降水量的變化沒有明顯總體年際上升或下降趨勢，處于波動(dòng)狀態(tài)，在1500—2500 mm之間，平均值為1913 mm，最高值出現(xiàn)在2016年，降水量2466 mm，最低值出現(xiàn)在2011年，降水量1426 mm。平均風(fēng)速趨于穩(wěn)定，年際間有波動(dòng)，平均值為 2.07 m·s-1，最高值出現(xiàn)在2011年，風(fēng)速值2.16 m·s-1，最低值出現(xiàn)在2001年，為1.85 m·s-1。日照時(shí)數(shù)趨于穩(wěn)定，年際間波動(dòng)，介于1494—1989之間，平均日照時(shí)數(shù)1696 h，最高值出現(xiàn)在2003年，最低值出現(xiàn)在2012年。

2.2.2 氣候因子與植被覆蓋變化趨勢的相關(guān)分析

采用相關(guān)系數(shù)法分析大灣區(qū)EVI均值與氣溫、降水、風(fēng)速、相對濕度和日照時(shí)數(shù)等5個(gè)氣候自然因子的相關(guān)性。通過逐像元計(jì)算 2001—2020年間EVI均值與氣候因子的相關(guān)系數(shù)，進(jìn)行顯著性檢驗(yàn)，得到 EVI變化趨勢與氣候因子的相關(guān)空間分布圖（圖4）和相關(guān)變化統(tǒng)計(jì)表（表3）。

表3 氣候因子影響EVI變化趨勢相關(guān)性面積占比統(tǒng)計(jì)Table 3 Correlation area statistics of climate factors influence the proportion of EVI trends

氣溫因子正相關(guān)影響區(qū)域占比49.16%，但顯著正相關(guān)的面積僅有7.15%。在空間分布上，氣溫顯著相關(guān)變化的植被區(qū)域呈極分散狀，多分布于城市建設(shè)的密集區(qū)域（圖4a）。相對濕度與EVI變化趨勢表現(xiàn)出顯著的正相關(guān)性，影響區(qū)域約80.23%，其中顯著正相關(guān)面積占比37.74%（r＞0，P＜0.05），影響的空間區(qū)域分布較廣，聚集分布在大灣區(qū)北部和東西兩側(cè)（圖4d）。降水量與EVI變化趨勢呈現(xiàn)出弱相關(guān)性，不顯著相關(guān)區(qū)域占比超過90%（圖4b）。EVI變化趨勢與風(fēng)速呈現(xiàn)出比較顯著的正向影響力，顯著相關(guān)影響區(qū)域32.99%（r＞0，P＜0.05），空間分布呈集聚性，集中于灣區(qū)東北、西北兩翼，多為林地，遠(yuǎn)離城鄉(xiāng)建設(shè)密集區(qū)域（圖 4c）。日照時(shí)數(shù)與EVI變化趨勢呈現(xiàn)出弱相關(guān)性，與EVI變化無顯著相關(guān)，其中不顯著負(fù)相關(guān)面積占比 68.22%（r＜0，P≥0.05）（圖 4e）。

圖4 EVI變化趨勢與氣候因子的相關(guān)性空間分析Figure 4 Correlation and significance between EVI trends and climate factors

2.3 非氣候自然因子對植被特征與變化趨勢的影響分析

2.3.1 高程、坡度、坡向與植被特征、變化趨勢的相關(guān)分析

采用 ArcGIS分級分析法研究高程和坡度等非氣候自然因子與EVI變化的相關(guān)性。在ArcGIS中利用自然間斷點(diǎn)法將高程分為＜82 m、82—217 m、217—394 m、394—653 m、＞653 m共5個(gè)梯度等級；依據(jù)1984年中國農(nóng)業(yè)區(qū)劃委員會(huì)頒發(fā)的《土地利用現(xiàn)狀調(diào)查技術(shù)規(guī)程》，將研究區(qū)坡度劃分為≤2°、2°—6°、6°—15°、15°—25°和＞25°共 5 個(gè)坡度等級，將研究區(qū)坡向均等劃分為北、東北、東、東南、南、西南、西、西北8個(gè)方向，各坡向所占面積比重較一致。分別將不同高程、不同坡度、不同坡向與2001—2020年逐年 EVI等級分布進(jìn)行空間耦合疊加，統(tǒng)計(jì)不同等級EVI均值隨不同高程梯度、不同坡度、不同坡向的分布情況（表4，圖5a）。利用逐像元計(jì)算2001—2020年間EVI與高程梯度、坡向、坡度的相關(guān)系數(shù)，通過 Sen趨勢分析和 Mann-Kendall顯著性檢驗(yàn)，結(jié)果劃分為7個(gè)等級，得到不同高程梯度、坡向、坡度的EVI均值特征和變化趨勢統(tǒng)計(jì)表（表5，圖5b）。

表4 高程梯度、坡度的EVI均值等級分布統(tǒng)計(jì)Table 4 Statistic of EVI values over different elevation gradients and different slope

表5 高程梯度、坡度影響EVI變化趨勢相關(guān)性面積占比統(tǒng)計(jì)Table 5 Elevation gradient and slope affect the proportion of EVI trends correlation area

研究結(jié)果表明，在植被分布特征上，伴隨高程的增加，EVI均值等級越高。高程≤82 m，EVI等級多在0.2—0.6之間；高程82—217 m，EVI等級在0.4—0.8之間，占比近79%；高程217 m以上，EVI等級集中在0.6—0.8之間，占比在85%以上。伴隨坡度的逐漸增大，EVI的等級均值越高。坡度≤2°時(shí)，EVI等級多在0.2—0.4之間。坡向因素對植被特征和變化趨勢的影響較弱，8個(gè)坡向的EVI分布情況及變化趨勢較一致（圖5）。

圖5 EVI分布特征、變化趨勢與坡向相關(guān)性分析Figure 5 Correlation between EVI Distribution (a), EVI trends (b) and different slope factors

在植被的變化趨勢上，不同高程梯度、不同坡度的EVI均表現(xiàn)為極顯著改善，其中改善程度最高的區(qū)域集中在82—394 m高程之間；伴隨坡度的增大，改善程度增大。在高程較低和坡度較緩地區(qū)，EVI等級均值和改善程度均低于高程較高和坡度較陡的區(qū)域，推測主要受人類活動(dòng)的影響，平地的EVI等級均值較低，改善趨勢較弱。

2.3.2 土壤類型與植被特征、變化趨勢的相關(guān)分析

土壤類型數(shù)據(jù)包含18個(gè)土類，赤紅壤、水稻土和紅壤是大灣區(qū)最主要的3種土壤類型，面積占比分別為49.31%、32.31%和10.62%。將土壤類型與逐年EVI均值等級分布進(jìn)行空間耦合疊加，統(tǒng)計(jì)不同等級EVI均值在不同土壤類型中的分布情況。借用逐像元計(jì)算EVI均值與土壤類型的相關(guān)系數(shù)，通過Sen趨勢分析和Mann-Kendall顯著性檢驗(yàn)，得到不同土壤類型中EVI變化的7種趨勢統(tǒng)計(jì)（表6）。

表6 不同土壤類型EVI分布與EVI變化趨勢Table 6 Distribution and trends statistic of EVI on the different type of soil

研究結(jié)果顯示，在植被特征上，城區(qū)土壤的EVI等級多為Ⅰ和Ⅱ類，分別占41.4%和51.19%；沙洲、島嶼等土壤類型的植被也多為Ⅰ和Ⅱ類，分別占33.33%和49.48%；赤紅壤、黃壤和紅壤的EVI等級均值高，赤紅壤EVI等級均值集中在0.4—0.8之間，黃壤和紅壤的EVI等級集中在0.6—0.8之間。在植被的變化趨勢上，各種土壤類型中對應(yīng)的植被生長狀況均表現(xiàn)為極顯著改善或顯著改善。濱海養(yǎng)殖場、城區(qū)、濱海鹽土、水稻土中的植被極顯著退化比例較高，均大于10%，分別為26.67%、14.58%、14.29%、12.11%；其他土壤類型中的極顯著退化和退化面積均小于或等于10%。

2.4 自然因子對植被覆蓋特征的驅(qū)動(dòng)力分析

采用地理探測器模型對氣候類和非氣候類自然因子影響植被空間格局特征的驅(qū)動(dòng)力進(jìn)行分析。以大灣區(qū)EVI均值作為因變量Y，以9個(gè)自然因子作為自變量因子X1,X2, ……,X9，如表7所示。在EVI多年均值水平上利用模型研究單個(gè)自然因子以及雙因子交互作用對EVI空間分布的驅(qū)動(dòng)力，探測結(jié)果如圖6、7所示。

圖6 單自然因子對EVI驅(qū)動(dòng)力分析Figure 6 Analysis of driving force of EVI by single natural factor

圖7 雙自然因子對EVI驅(qū)動(dòng)力的交互作用探測Figure 7 Detection of interaction of two natural factors on EVI driving force

表7 自變量編碼Table 7 Independent variable coding

單因子的驅(qū)動(dòng)力分析表明，9個(gè)單因子驅(qū)動(dòng)力影響排序?yàn)椋焊叱蹋練鉁兀鞠鄬穸龋就寥李愋停酒露龋撅L(fēng)速＞降水量＞日照時(shí)數(shù)＞坡向。4個(gè)非氣候類自然因子中，高程、土壤類型、坡度是影響植被覆蓋特征的主要因子。其中，高程對EVI空間格局的解釋力超過50%；土壤類型解釋力接近40%；坡度解釋力約34%；坡向解釋力則在10%及以下，EVI受坡向的影響很弱。5個(gè)氣候因子中，年均氣溫、年均相對濕度、風(fēng)速是影響植被覆蓋特征的主要因子，年均氣溫、年均相對濕度對EVI空間格局的解釋力均在40%以上；風(fēng)速解釋力在30%以上，降水解釋力低于20%，日照因子解釋力低于10%。

在雙因子協(xié)同驅(qū)動(dòng)作用力上，大部分因子之間的交互作用表現(xiàn)為線性關(guān)系，雙因子增強(qiáng)，但增強(qiáng)的程度存在差異。唯一的非線性關(guān)系主要是降水和日照的協(xié)同關(guān)系，但影響力極弱。其中，非氣候類自然因子中，高程（X6）與其他8個(gè)自然因子的交互作用均表現(xiàn)為驅(qū)動(dòng)力協(xié)同增強(qiáng)，雙因子增強(qiáng)q值均大于0.5，最大值0.6669出現(xiàn)在高程和相對濕度的協(xié)同驅(qū)動(dòng)，雙因子交互驅(qū)動(dòng)的解釋力接近70%，表明高程和相對濕度的交互驅(qū)動(dòng)作用對EVI空間格局影響的解釋力最強(qiáng)。其次是土壤因子（X9），與其他8個(gè)自然因子的交互作用也表現(xiàn)為驅(qū)動(dòng)力協(xié)同增強(qiáng)，雙因子增強(qiáng)q值均大于0.4，最大值0.6002出現(xiàn)在土壤類型和相對濕度的協(xié)同驅(qū)動(dòng)。坡向和其他各因子的協(xié)同關(guān)系均比較弱。在氣候因子方面，氣溫（X1）、風(fēng)速（X3）、相對濕度（X4）與其他自然因子都表現(xiàn)出協(xié)同增強(qiáng)的驅(qū)動(dòng)力，其中，氣溫、相對濕度與各因子的雙因子驅(qū)動(dòng)解釋力均在50%左右，風(fēng)速與相對濕度、高程、土壤類型、坡度的雙因子驅(qū)動(dòng)解釋力也均在50%左右。

3 討論

粵港澳大灣區(qū)是中國開放程度最高、經(jīng)濟(jì)活力最強(qiáng)的區(qū)域之一，在國家發(fā)展大局中具有重要戰(zhàn)略地位。在20年發(fā)展歷程中，始終面臨生態(tài)環(huán)境保護(hù)和開發(fā)建設(shè)的雙重需求和壓力。從EVI的逐年均值變化看，大灣區(qū)的植被覆蓋年際變化總體趨勢是波動(dòng)上升，這一結(jié)果與珠三角植被覆蓋變化趨勢一致（何全軍，2019；王思等，2022）。雖植被覆蓋以改善為主，但同時(shí)存在退化。由于所采用的指數(shù)和具體方法不同，退化的面積占比存在差異，對顯著退化面積的研究結(jié)果包括大灣區(qū)顯著退化3.72%（趙桔超等，2022）、珠三角顯著退化 9.16%（何全軍2019）、珠三角九市顯著退化17.6%（王思等，2022），大灣區(qū)顯著退化20.95%（馮嫻慧等，2022）。結(jié)合同時(shí)期（2000—2020）大灣區(qū)土地利用變化研究，大灣區(qū)的用地類型仍以林地為主，但居住和建設(shè)用地增加9%，林地、草地、耕地減少面積約8%（高照忠等，2021）。2015年，林地在用地類型中占比約53.99%（張杰等，2021），而到2020年，占比約為49%（高照忠等，2021）?？梢酝茢?，大灣區(qū)的植被退化面積在10%—20%區(qū)間，極顯著退化的面積在10%以下。經(jīng)濟(jì)發(fā)展對于環(huán)境改善有一定促進(jìn)作用，城市土地開發(fā)配套30%綠地硬性指標(biāo)、新建城市公園、廣東省綠道建設(shè)、綠化廣東大行動(dòng)等，均有效改善了植被退化的程度。從空間分布看，大灣區(qū)年最大EVI的均值分布具有明顯的空間差異性和集聚性，呈現(xiàn)出“外高中低”的分布特點(diǎn)，與前人研究結(jié)果一致（趙桔超等，2022；何全軍2019）。

氣候類自然因子對植被變化的驅(qū)動(dòng)影響力中，已有的研究集中于氣溫、降水、日照時(shí)數(shù)3個(gè)因子。不同氣候區(qū)域的氣候因子對植被的影響程度存在差異，氣溫、降水因子對植被具有顯著正向影響的地區(qū)有內(nèi)蒙古（穆少杰等，2012）、京津冀地區(qū)（晏利斌等，2011；徐勇等，2020）和黃河流域（解晗等，2022）。張顧萍等（2021）認(rèn)為半干旱區(qū)的植被覆蓋變化主要受降水影響。而在浙江的研究表明，氣溫是引起植被覆蓋變化的主要因子（賀忠華等，2020）；在廣東省以及大灣區(qū)，研究表明氣溫對植被生長的影響最顯著（鄧玉嬌等，2021；趙桔超等，2022）。本研究發(fā)現(xiàn)，在高溫多雨的大灣區(qū)，除氣溫因子外，相對濕度、風(fēng)速也是影響植被的重要因子，其影響程度超過了降水、日照因子。由此推斷，高溫多雨地區(qū)的植被生長特征與變化受到降雨因子的影響遠(yuǎn)低于半干旱地區(qū)，但空氣濕度、風(fēng)速對植被生長有更顯著的影響。

本研究雖未涉及人為影響因子，但從自然因子的分析看，高程、土壤類型、坡度對土地可供開發(fā)利用具有重要影響。在土壤類型中，EVI發(fā)生了顯著變化的土地利用類型主要是城區(qū)、濱海土類、水稻土，可印證20年來，大灣區(qū)土地利用類型變化以耕地向城鄉(xiāng)用地轉(zhuǎn)變?yōu)橹鞯陌l(fā)展?fàn)顩r，這與杭州灣、浙江省的研究結(jié)果具有一致性（董晨煒等，2017；賀忠華等，2020），進(jìn)一步推斷，海岸帶的開發(fā)和建設(shè)，也是造成灣區(qū)植被化的原因之一。

4 結(jié)論

本研究基于2001—2020年的MODIS NDVI遙感數(shù)據(jù)，采用增強(qiáng)植被指數(shù)（EVI），結(jié)合Sen趨勢分析和Mann-Kendall檢驗(yàn)、Pearson相關(guān)系數(shù)法、分級分析法、地理探測器驅(qū)動(dòng)力分析模型等方法，分析 20年長時(shí)間序列的大灣區(qū)植被覆蓋特征與變化趨勢，并進(jìn)一步結(jié)合氣象數(shù)據(jù)、高程、坡度、坡向、土壤類型等數(shù)據(jù)，研究9個(gè)氣候與非氣候自然因子對植被覆蓋特征與變化趨勢的影響和驅(qū)動(dòng)力。研究發(fā)現(xiàn)，（1）大灣區(qū)植被覆蓋整體呈改善趨勢，植被改善區(qū)的面積大于退化區(qū)面積。在不同土地利用類型上，植被變化趨勢均表現(xiàn)為極顯著改善。其中，城鄉(xiāng)建設(shè)用地、水域、耕地的植被變化趨勢表現(xiàn)為以改善為主，但改善和退化等7種不同程度的變化趨勢共存，城鄉(xiāng)建設(shè)用地的極顯著退化趨勢最高。在植被變化趨勢上，城鄉(xiāng)建設(shè)用地、栽培植被的改善、退化趨勢比較復(fù)雜。（2）大灣區(qū)的氣候變化主要變現(xiàn)在氣溫、相對濕度的緩慢上升，相對濕度、風(fēng)速、氣溫對植被覆蓋變化表現(xiàn)出顯著影響力；但在空間分布上，氣溫影響區(qū)域與城市建設(shè)密集區(qū)域緊密相關(guān)。（3）在非氣候因子中，高程、坡度與植被分布特征緊密相關(guān)。在植被的變化趨勢上，改善程度最高的區(qū)域集中在82—394 m高程之間，退化區(qū)域的高程≤82 m，不同坡度的EVI改善程度伴隨坡度增大而增大。植被覆蓋發(fā)生顯著退化占比面積較大的有城區(qū)、濱海土類、水稻土。高程、土壤類型、坡度是導(dǎo)致植被分布變化的主要驅(qū)動(dòng)因子群。（4）通過驅(qū)動(dòng)力模型分析，5個(gè)氣候因子中，年均氣溫、年均相對濕度、風(fēng)速是影響植被空間特征的主要因子；4個(gè)非氣候因子中，高程、土壤類型是影響植被空間特征的主要因子。在雙因子交互驗(yàn)證中，氣溫、相對濕度與各因子的雙因子驅(qū)動(dòng)解釋力均在 50%左右，表明高溫高濕的氣候?qū)Υ鬄硡^(qū)的植被生長具有正向驅(qū)動(dòng)作用。高程、土壤因子與其他自然因子的交互作用均表現(xiàn)為驅(qū)動(dòng)力協(xié)同增強(qiáng)。