吳宇凡，王江波，茍愛萍

（1.南京工業(yè)大學(xué) 建筑學(xué)院，南京 211816；2.上海應(yīng)用技術(shù)大學(xué) 生態(tài)技術(shù)與工程學(xué)院，上海 201418）

隨著城市化進程加快，人類活動對地表水體的影響程度不斷提高，特別是在城市化進程最快的長三角、珠三角地帶，該區(qū)域是江、河、湖、海等水體交互作用最為活躍，也是地表水體對人類活動最為敏感的地區(qū)（楊達源 等，2009）。粵港澳大灣區(qū)作為繼紐約、舊金山、東京灣區(qū)之后的世界第四大灣區(qū)，中國三大城市群之一，其快速城市化在帶來巨大經(jīng)濟效益的同時，也對地表水體造成較大程度的影響，部分區(qū)域出現(xiàn)大面積水體萎縮，同時引發(fā)水土流失、水質(zhì)污染等生態(tài)問題。因此，探究粵港澳大灣區(qū)地表水體演變特征及相關(guān)影響因素是促進地表水體保護與合理利用、維持社會可持續(xù)發(fā)展的重點。

遙感技術(shù)和地理信息技術(shù)的快速發(fā)展為開展水體監(jiān)測與研究提供了強力的技術(shù)支撐。目前地表水體研究主要關(guān)注水體時空演變、水質(zhì)變化監(jiān)測等方向。在水體時空演變領(lǐng)域，學(xué)者多側(cè)重于從宏觀與微觀空間尺度開展研究。宏觀方面，如嚴登華等（2020）在提出不同分區(qū)流域產(chǎn)水系數(shù)移植方法的基礎(chǔ)上，評估了1960—2016年全球陸地地表水資源量并分析其演變特征；莫貴芬等（2022）通過動態(tài)度法等方法研究了2001—2018年中亞干旱區(qū)地表水資源時空變化特征，認為中亞干旱區(qū)內(nèi)氣候因子對地表水體影響不顯著。微觀層面，當(dāng)前學(xué)者多圍繞地表流域、保護區(qū)或單個市縣開展研究，如針對喀斯特山區(qū)（董飛 等，2022）、遼河（閆曉寒 等，2022）、新疆開孔河（孟記超，2021）、甘肅黑河（閆柏琨 等，2022）、黃土高原莊浪縣（曹杰 等，2021）等，研究成果頗豐，但對于城市群視角下長時序水體演變及深層次驅(qū)動因素的研究尚未深入。

水體提取層面，多采用水體指數(shù)如NDWI（胡雪松 等，2019）、MNDWI（梅元勛 等，2018）結(jié)合人機交互提取本地遙感數(shù)據(jù)中的水體或岸線，如Tao 等（2015）借助Landsat 系列衛(wèi)星數(shù)據(jù)，采用NDWI 探索了蒙古高原湖泊面積和湖泊數(shù)量的變化；程彥培等（2010）利用1989、1999、2009年3個時期的MODIS 數(shù)據(jù)，依據(jù)光譜特性提取了中亞地區(qū)的地表水體數(shù)據(jù)。隨著谷歌地球引擎（Google Earth Engine）云計算平臺的研發(fā)與興起，GEE 強大的交互能力與高算力等優(yōu)點彌補了以往傳統(tǒng)遙感技術(shù)在處理長時序、大范圍空間數(shù)據(jù)的不足，已成為當(dāng)前水體研究領(lǐng)域的高效工具。如Zou 等（2017）使用GEE 平臺，通過調(diào)用可獲取的Landsat5、7、8 系列影像，獲取美國1984—2016 年地表水?dāng)?shù)據(jù)集并量化其時空變化程度；Pekel 等（2016）利用GEE 平臺，使用專家系統(tǒng)水體提取方法，生產(chǎn)出1984—2020 年間年尺度以及月尺度全球范圍內(nèi)30 m 空間分辨率的地表水覆蓋數(shù)據(jù)集，相較于該領(lǐng)域的美國NASA SRTM Water Body Data（SWBD）、加拿大Global Lakes and Wetlands Database （GLWD）、日本Global 3 Arc-Secend Water Body Map（G3WBM）等數(shù)據(jù)集，在精度、時間跨度、公開性上都有明顯優(yōu)勢，后續(xù)研究（于濤，2019；趙金玉，2020）基于該數(shù)據(jù)集也取得良好的結(jié)果。

對于驅(qū)動因素的探討多采用統(tǒng)計分析（Yang et al., 2021）、主成分分析（Hu et al., 2020）、相關(guān)性分析、Mann-Kendall檢驗、Morlet小波分析（歐文杰 等，2020）等數(shù)學(xué)方法來量化各因子的驅(qū)動程度。但眾多方法的評價對象多為單個自變量對因變量的驅(qū)動解釋力，忽略了驅(qū)動過程中的因子整體性。王勁峰等（2017）提出地理探測器模型可以更好量化地理空間中自變量對因變量的解釋力，該方法能精準計算單因子以及各因子交互作用對因變量空間分異的解釋力，目前在自然地理（裴志林 等，2019；丁倩 等，2021）、人文地理（劉彥隨 等，2017；李在軍 等，2020）、應(yīng)用生態(tài)分析（黃欣等，2022）等方面得到廣泛應(yīng)用。

因此，本文選取粵港澳大灣區(qū)作為研究區(qū)，在通過GEE 平臺調(diào)用1985—2020 年全球地表水?dāng)?shù)據(jù)集（GSW）計算灣區(qū)地表水體面積變化速率、格網(wǎng)水體密度基礎(chǔ)上，使用土地利用轉(zhuǎn)換矩陣和地理探測器模型，結(jié)合各時期城市發(fā)展策略與區(qū)域發(fā)展導(dǎo)向，從城市發(fā)展進程的角度探究1985—2020年粵港澳大灣區(qū)地表水體時空演變特征與相關(guān)驅(qū)動因素。以期為更好掌握大灣區(qū)地表水體時空變化情況提供數(shù)據(jù)支撐，為相關(guān)部門加強水體保護與管理、維護城鄉(xiāng)可持續(xù)發(fā)展提供科學(xué)依據(jù)。

1 研究區(qū)、數(shù)據(jù)與方法

1.1 研究區(qū)概況

粵港澳大灣區(qū)（The Guangdong-Hong Kong-Macao Greater Bay Area, GBA）共包括廣州、深圳、珠海、佛山、惠州、東莞、中山、江門、肇慶9個地級市和香港、澳門2個特別行政區(qū)。地處于珠江下游，瀕臨南海，總占地面積55 900 km2，屬南亞熱帶濕潤季風(fēng)氣候，四季分明，雨熱同期，水系眾多。大灣區(qū)以國際一流灣區(qū)和世界級城市群作為發(fā)展目標，目前大灣區(qū)經(jīng)濟體量已逾12 萬億人民幣（廣東省人民政府，2023）。

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

選取1985—2020 年全球地表水年度數(shù)據(jù)集（JRC Yearly Water Classification History, v1.3）作為主要數(shù)據(jù)源，該數(shù)據(jù)集是歐盟聯(lián)合研究中心（JRC）利 用1984—2020 年 的300 多 萬 景Landsat TM/ETM+/OLI衛(wèi)星遙感影像數(shù)據(jù)，基于GEE平臺計算獲得的1984—2020 年30 m 分辨率全球地表水年度分布數(shù)據(jù)。數(shù)據(jù)集將水體像元數(shù)據(jù)分為永久性水體和季節(jié)性水體兩類，劃分標準為觀測周期內(nèi)的水體復(fù)現(xiàn)頻率SWO（Surface Water Occurrence），在1年中對該像元的所有觀測均判定為水體（SWO=100%），即為永久性水體，如水庫、湖泊等；否則，即判定該像元為季節(jié)性水體（1%＜SWO＜100%），如洪澇水體等。Pekel等（2016）通過對GSW 數(shù)據(jù)集的精度驗證，得出該數(shù)據(jù)集具有較高的水體提取精度，其中全球總漏檢率＜5%，總誤檢率＜1%，滿足數(shù)據(jù)研究精度。年度數(shù)據(jù)基于Landsat 系列衛(wèi)星在該年12個月所有可獲取數(shù)據(jù)，相比月度數(shù)據(jù)極大降低了數(shù)據(jù)受云量過高、衛(wèi)星未過境等問題的影響，更全面地反映年度水體分布，因此，選取1985—2020年每5年的全球地表水年度數(shù)據(jù)作為主要數(shù)據(jù)源。

輔助數(shù)據(jù)主要為全國土地利用類型遙感監(jiān)測空間分布數(shù)據(jù)（LULC）、高程、坡度、年均氣溫、年均降水量、NDVI、植被類型、土壤類型、國民生產(chǎn)總值、人口總量、耕地密度、建設(shè)用地密度、交通路網(wǎng)數(shù)據(jù)。其中，ASTER GDEM 30 m 分辨率數(shù)字高程數(shù)據(jù)來自地理空間數(shù)據(jù)云①http://www.gscloud.cn/；坡度數(shù)據(jù)由高程數(shù)據(jù)計算得出；交通路網(wǎng)數(shù)據(jù)來自O(shè)pen Street Map網(wǎng)站；GDP、人口、年均氣溫等數(shù)據(jù)來自中國科學(xué)院資源環(huán)境科學(xué)與數(shù)據(jù)中心②https://www.resdc.cn/Default.aspx；粵港澳大灣區(qū)行政邊界來自國家基礎(chǔ)地理信息中心③https://www.webmap.cn/2021 年發(fā)布的1∶100 萬公眾版基礎(chǔ)地理信息數(shù)據(jù)集；統(tǒng)計數(shù)據(jù)如城市化率等來自國家統(tǒng)計局、廣東省統(tǒng)計局以及研究區(qū)各地市統(tǒng)計局，詳細數(shù)據(jù)參數(shù)見表1所示。

表1 地表水體及相關(guān)驅(qū)動因素選取數(shù)據(jù)參數(shù)Table 1 Selected data parameters for surface water bodies and relevant driving factors

1.3 指標選取與方法

選取水體面積變化速率、水體密度2類指標對大灣區(qū)地表水體規(guī)模進行動態(tài)量化分析，使用土地利用轉(zhuǎn)換矩陣和地理探測器分析地表水體演化過程并探測其相關(guān)驅(qū)動因素。

1.3.1 面積變化速率 通過地表水體變化面積（ΔS）、變化速率（v）指標量化灣區(qū)地表水體面積變化情況，數(shù)學(xué)表達式為：

式中：Sa和Sb分別為研究期間前期和后期的地表水面積；ΔT為時間間隔。

1.3.2 水體密度 為進一步探究地表水體的動態(tài)空間變化，使用水體密度指標（D）及空間可視化來量化和表現(xiàn)地表水體在長時間序列中的空間演化情況。數(shù)學(xué)表達式為：

式中：GSa和GSb分別為單個格網(wǎng)（Grid）內(nèi)地表水體面積和該格網(wǎng)總面積。通過在GEE平臺生成1 km2的格網(wǎng)覆蓋地表水體空間范圍，統(tǒng)計網(wǎng)格內(nèi)地表水體面積占網(wǎng)格面積比例并進行可視化處理，以分析地表水體區(qū)域密度的空間演化。

1.3.3 土地利用轉(zhuǎn)移矩陣 土地利用轉(zhuǎn)移矩陣是馬爾科夫模型在土地利用變化方面的應(yīng)用，可定量描述系統(tǒng)內(nèi)的狀態(tài)變化（劉紀遠 等，2000）。運用ArcGIS對相鄰年份土地利用數(shù)據(jù)進行疊加，得到反映各土地利用類型相互之間的轉(zhuǎn)化數(shù)據(jù)，得出土地利用轉(zhuǎn)換矩陣，用以反映各用地類型轉(zhuǎn)換的數(shù)量變化。轉(zhuǎn)換矩陣的數(shù)學(xué)原理為：

式中：Cij為相鄰年份首年的i種土地利用類型轉(zhuǎn)變?yōu)橄噜從攴菽┠甑膉種土地利用類型的面積。

1.3.4 地理探測器 空間分異是自然和社會經(jīng)濟過程的空間表現(xiàn)，地理探測器是探測空間分異性，以及揭示其背后驅(qū)動因子的一種統(tǒng)計學(xué)方法（王勁峰等，2017）。該模型由因子探測、生態(tài)探測、交互作用探測和風(fēng)險因子探測4種探測器構(gòu)成，在衡量自變量對因變量空間分布的解釋力上相互完善與支撐，本研究主要針對單因子探測和交互作用探測2種探測類型展開研究。

分異及因子探測：探測研究區(qū)地表水體面積（Y）的空間分異性以及探測各因子X在多大程度上解釋了屬性Y的空間分異。用q值度量，數(shù)學(xué)表達式為：

式中：h=1, …,L為變量Y或因子X的分層；Nh和N分別為層h和全區(qū)的單元數(shù)；和σ2分別是層h和全區(qū)的Y值的方差。SST和SSW分別為全區(qū)總方差和層內(nèi)方差之和。q的值域為[0, 1]，值越大說明Y的空間分異性越明顯；如果分層是由自變量X生成的，q值越大表示自變量X對屬性Y的解釋力越強，反之越弱。

交互作用探測：用于評估影響因子之間的交互作用，通過比較單一因子q值、雙因子之和q值以及雙因子交互作用后的q值，判斷兩因子在交互作用后增大還是減弱了對地表水體空間分布的影響，將交互類型分為非線性減弱、單因子非線性減弱、雙因子增強、獨立關(guān)系和非線性增強5類。

選取自然與社會經(jīng)濟2方面共12個指標探測地表水體面積與相關(guān)因子之間的關(guān)聯(lián)程度（表2），因子選取具有較明顯的自然與社會經(jīng)濟屬性，對于驅(qū)動因素探究來說具有較強的代表性。數(shù)據(jù)處理全程使用ArcGIS Pro 2.8.6 與Matlab 操作，首先生成1 km2的格網(wǎng)覆蓋研究區(qū)全域，接著對各因子統(tǒng)一采用自然間斷點分類方法進行離散化處理，將統(tǒng)計數(shù)據(jù)導(dǎo)入Excel 與Matlab 中進行數(shù)據(jù)清洗后計算Spearman相關(guān)系數(shù)進行結(jié)果驗證，其他年份在替換相應(yīng)年份數(shù)據(jù)后重復(fù)上述操作。

表2 本研究地理探測器構(gòu)建參數(shù)Table 2 Geodetector parameters of this study

2 結(jié)果與分析

2.1 演變階段劃分

在1985-2020年，粵港澳大灣區(qū)全域地表水體面積變化速率呈現(xiàn)“緩慢下降—加速下降—減速放緩—穩(wěn)步提升”的階段性特征（圖1-a、表3），據(jù)此將研究時段劃分為4個階段：1）“水體緩慢萎縮”階段（1985-2000年），該階段地表水體緩慢萎縮，變化速率為-30.69～-4.66 km2/a，地表水體面積共減少274.04 km2；2）“水體加速萎縮”階段（2000-2010 年），該階段地表水體變化速率由-40.1 km2/a加速降至-84.18 km2/a，在此10 年內(nèi)地表水體面積大幅度減少621.44 km2，減速在2010 年達到頂峰；3）“水體減速放緩”階段（2010-2015 年），在5年內(nèi)地表水體減速由上一階段的-84.18放緩至-8.65 km2/a，但水體面積仍在減少；4）“水體穩(wěn)步增長”階段（2015-2020 年），該階段地表水體面積轉(zhuǎn)減為增，2020較2015年增加了444.28 km2，面積變化速率首次達到88.86 km2/a。

表3 1985—2020年粵港澳大灣區(qū)地表水體面積變化速率Table 3 The rate of surface water change area in the GBA from 1985 to 2020

圖1 粵港澳大灣區(qū)地表水體面積變化Fig.1 Change of surface water area change in GBA

從大灣區(qū)內(nèi)各城市角度上，各區(qū)域總體變化趨勢相同，部分區(qū)域變化趨勢略有不同（圖1-b）。在水體緩慢萎縮階段（1985-2000 年），除惠州市水體面積保持增長狀態(tài)，其余區(qū)域均呈現(xiàn)面積萎縮狀態(tài)，其中，佛山、廣州兩地水體面積降速較高，分別為-26.75 和-8.35 km2/a，肇慶、江門、珠海、中山、深圳5 市降速較高，降速在-7.06～-3.76 km2/a；水體加速萎縮階段（2000-2010 年），各市與特別行政區(qū)水體面積均呈現(xiàn)快速下降特征，其中，佛山仍保持最高降速，由上一階段的-26.75 km2/a加速降至-31.6 km2/a，江門、肇慶、珠海、廣州、惠州、深圳、東莞7 市降速也分別加速降至-23.01、-13.28、-12.1、-10.08、-9.46、-9.36、-8.96 km2/a。水體減速放緩階段（2010-2015 年），各市與特別行政區(qū)水體面積仍在減少，但減速較上一階段大幅度放緩，且出現(xiàn)城市組團的增長勢頭，江門、佛山、肇慶3 市水體面積變化速率分別增至8.1、7.55、6.28 km2/a；水體穩(wěn)步增長階段（2015-2020年），全域各市與特別行政區(qū)水體面積變化速率均呈穩(wěn)步增長態(tài)勢，其中，江門、中山、佛山、廣州4市速率增至14.18、13.86、13.74、11.09 km2/a。從空間上看，地表水體的大幅變化一般以區(qū)域城市組團的空間形式呈現(xiàn)，如水體緩慢萎縮階段的佛山、廣州、肇慶組團、水體加速萎縮階段的佛山、江門、肇慶組團、水體減速放緩階段的江門、佛山、肇慶組團和水體穩(wěn)步增長階段的江門、中山、佛山、廣州組團。

從研究區(qū)全域以及各城市角度對永久性水體和季節(jié)性水體面積進行統(tǒng)計分析（圖1-c、d）得出，永久性水體面積變化速率呈“小周期交替增減”特征，而季節(jié)性水體面積變化速率則呈“大幅度先減后增”特征（圖1-e、f）。大灣區(qū)永久性水體面積在各研究時段均呈“小周期豐枯交替”特征，而大灣區(qū)季節(jié)性水體面積在2010年以前多處于萎縮狀態(tài)，在2010年后得到大幅度提升。受初始面積占總面積比例的影響，2010年前，大灣區(qū)地表水體面積變化受永久性水體面積變化影響較大，變化趨勢與永久性水體變化趨勢基本一致；2010年后，大灣區(qū)水體面積變化與季節(jié)性水體變化不斷趨同。研究期內(nèi)，永久性水體與季節(jié)性水體對總體面積變化呈“交替支配”特征，這主要是由于永久性水體面積變化幅度在各時段逐漸低于季節(jié)性水體面積的變化幅度，且永久性水體面積在交替增減中逐漸下降，而季節(jié)性水體面積在2010 年后保持大幅度增加勢頭所造成的。

2.2 水體密度分析

通過格網(wǎng)化分析水體密度的方法進一步探究地表水體的空間變化過程，對水體密度格網(wǎng)采用相同間隔（20%）分為高、較高、中、較低、低5 種密度分區(qū)，得出各年份粵港澳大灣區(qū)水體密度分布（圖2）。

圖2 粵港澳大灣區(qū)格網(wǎng)水體密度分區(qū)空間分布Fig.2 Distribution of grid water density in GBA

在“水體緩慢萎縮”階段（1985-2000 年），佛山中部和東南部水體密度大幅度降低，大量季節(jié)性水體萎縮或消失，廣州中部及南部也存在較大范圍季節(jié)性水體的消失，其余地區(qū)如肇慶、江門、珠海等市也都存在不同區(qū)域的低密度區(qū)水體消失，該現(xiàn)象與上述水體面積變化相互佐證；萎縮或消失水體在該階段多出現(xiàn)在中、較高、高水體密度區(qū)域。2000-2010年，隨著各地城鎮(zhèn)化進程不斷加速，水體面積加速萎縮，大灣區(qū)全域均存在大幅度水體萎縮與消失現(xiàn)象，廣佛肇、深莞惠、珠中江等多個區(qū)域城市群呈現(xiàn)同步加速萎縮現(xiàn)象，水體萎縮與消失區(qū)多為中、較低、低密度水體區(qū)域，原本地表水面積較少的區(qū)域被其他類型用地占用，大范圍的填湖、填海工程使高密度水體區(qū)進一步加速萎縮。

2010-2015年，大灣區(qū)地表水體面積減速大幅度放緩，沿海城市如珠海、深圳、廣州南端的沿海區(qū)域，高密度的水域存在密度降低或水體消失現(xiàn)象，這主要是該階段的填海工程、連島工程所導(dǎo)致的。但江門、佛山、肇慶等城市已出現(xiàn)水體密度增長勢頭，低、較低、中、較高水體密度區(qū)域呈現(xiàn)較大幅度的增長狀態(tài)。進入2015年后，大灣區(qū)水體面積轉(zhuǎn)減為增，各類型水體進入穩(wěn)步增長階段，大灣區(qū)各市與特別行政區(qū)各密度水域均出現(xiàn)較大幅度增長，如江門、珠海、中山、廣州、佛山等城市形成的城市組團，低、較低、中、較高密度水體均呈現(xiàn)同步增長態(tài)勢。

2.3 用地轉(zhuǎn)換分析

為了進一步探究地表水體的地類轉(zhuǎn)換情況，挖掘長時序下驅(qū)動灣區(qū)地表水體面積變化的相關(guān)因素，對粵港澳大灣區(qū)進行1980-2020年各時間段土地利用轉(zhuǎn)換矩陣分析（表4）。數(shù)據(jù)源選取中國科學(xué)院出版的全國土地利用類型遙感監(jiān)測空間分布數(shù)據(jù)（徐新良 等，2018）該數(shù)據(jù)土地利用一級類型綜合評價精度達到93%以上，二級類型分類綜合精度達90%以上，滿足研究精度要求，數(shù)據(jù)選取時間間隔為10 a。

表4 1980-2020年粵港澳大灣區(qū)水域轉(zhuǎn)入轉(zhuǎn)出Table 4 Statistical table of water transfer in and out of GBA from 1980 to 2020 km2

1980-1990 年，隨著改革開放的開始，深圳、珠海、汕頭、廈門4城設(shè)立經(jīng)濟特區(qū)，各城市開始城市化進程，大灣區(qū)逐漸進入地表水體緩慢萎縮階段，10 a 內(nèi)水域轉(zhuǎn)出556.95 km2，其中轉(zhuǎn)入耕地與城鄉(xiāng)建設(shè)用地居多，分別為244.02 與208.89 km2，而轉(zhuǎn)入水域共385.02 km2，其中耕地、海域分別轉(zhuǎn)入268.4和113.24 km2，水域面積凈減171.93 km2。

1990-2000年，大灣區(qū)處于地表水體緩慢萎縮階段，隨著城市化進程不斷加快，水域分別向城鄉(xiāng)建設(shè)用地、耕地、林地分別轉(zhuǎn)出327.26、177.71、112.58 km2，城鄉(xiāng)建設(shè)用地取代耕地成為水域轉(zhuǎn)出面積最高地類，10 a水域凈減316.08 km2。

2000-2010年，大灣區(qū)進入地表水體加速萎縮階段，水域共轉(zhuǎn)出1 380.04 km2，其中，城鄉(xiāng)建設(shè)用地、耕地、林地分別轉(zhuǎn)出631.23、449.39、244.74 km2，快速城市化帶來大量對城市建設(shè)用地的需求，促使各主體通過填湖填塘等手段獲取用地，此期間轉(zhuǎn)入水域面積僅724.69 km2，其中耕地（62%）為主要轉(zhuǎn)入來源，這主要是由于大灣區(qū)水系密集，大量耕地被當(dāng)?shù)鼐用裼米魉a(chǎn)養(yǎng)殖等用途，水域面積凈減555.71 km2。

2010-2020年，大灣區(qū)地表水體由減速放緩階段過渡到穩(wěn)步增長階段，各城市發(fā)展陸續(xù)進入穩(wěn)定階段，各城市對建設(shè)用地需求降低；同時環(huán)境保護力度不斷加大、城市藍線等各種城市建設(shè)管控條件更加嚴格，期間水域共轉(zhuǎn)出526.68 km2，轉(zhuǎn)入水域面積752.90 km2，水域面積凈增226.22 km2，轉(zhuǎn)入地類中耕地轉(zhuǎn)入404.28 km2、林地轉(zhuǎn)入215.86 km2、城鄉(xiāng)建設(shè)用地轉(zhuǎn)入116.77 km2。

2.4 探測結(jié)果分析

為了進一步探究不同階段驅(qū)動大灣區(qū)地表水體空間分異的影響因子，選擇水體變化階段之間的關(guān)鍵時間節(jié)點，即2000、2010、2015、2020年進行地理探測器分析，選取自然與社會經(jīng)濟2方面共12個指標進行各年份的單因子探測和交互因子探測，并結(jié)合Spearman系數(shù)進行數(shù)據(jù)驗證。

2.4.1 單因子探測 在自然因子中，年降水量、植被覆蓋率、年均溫、土壤類型、高程等自然因子在對地表水體空間分異擁有較高的解釋力，解釋力分別為40%、34%、32%、30%、23%（表5），說明在長時序角度下，自然條件仍是影響地表水體空間分異的重要因子。高程決定的地表形態(tài)直接驅(qū)動水體產(chǎn)生空間分異；降水能直接作用于地表徑流，匯入永久性水體后進入海洋，同時也是城市與農(nóng)村重要的用水來源；植被覆蓋率高的區(qū)域水土保持能力較強，能起較好的水源涵養(yǎng)作用；年均溫能顯著影響年內(nèi)植被、土壤、水體的蒸散發(fā)過程從而以降水形式補給區(qū)域水體；研究區(qū)水系密度高的區(qū)域往往存在大量潮土、水稻土、淺育水稻土，多為位于河網(wǎng)兩側(cè)的水田或?；~塘，方便灌溉的同時起水土涵養(yǎng)作用；可以得出，高程、降水、植被、溫度、土壤等要素共同影響著地表水體的空間分異。在社會經(jīng)濟因子中，耕地密度、人口對地表水體空間分異擁有較高的解釋力，q值分別為31%和30%，灣區(qū)存在大量水田及?；~塘，這也是季節(jié)性水體的主要組成部分，水田、魚塘多位于河道兩旁以方便灌溉與取水。社會經(jīng)濟指標中的建設(shè)用地密度、地區(qū)產(chǎn)值、交通路網(wǎng)長度的解釋力較低，分別為15%、13%、5%。

表5 各年份單因子探測結(jié)果Table 5 List of single factor detection results in each year

從各階段因子解釋力變化上看，在大灣區(qū)水體萎縮階段，人口、耕地密度因子解釋力分別由27%、30%上升至31%、31%，從側(cè)面說明快速城市化促使大灣區(qū)工業(yè)、農(nóng)業(yè)的快速發(fā)展；在地表水體減速放緩和穩(wěn)步增長階段，社會經(jīng)濟等指標如地區(qū)生產(chǎn)總值因子解釋力由15%降至12%、人口因子由31%降至29%、耕地密度、建設(shè)用地密度也均有下降，從側(cè)面反映灣區(qū)城市整體發(fā)展方向與政策的轉(zhuǎn)向，越來越注重人與自然的協(xié)調(diào)關(guān)系，也間接驅(qū)動大灣區(qū)地表水體面積的增長。

2.4.2 因子交互探測 圖3表明，各因子存在雙因子增強和非線性增強2種交互方式，即因子交互的解釋力均大于單個因子的解釋力。在自然因子中，年降水量（X4）和植被覆蓋率（X5）交互解釋力達到86%、年均溫（X3）和植被覆蓋率（X5）為81%、年降水量（X4）與耕地密度（X9）為79%。年降水量與除建設(shè)用地密度和交通路網(wǎng)長度之外的所有因子均呈現(xiàn)雙因子增強的交互關(guān)系，年降水量顯著增強各因子對于地表水體空間分異的影響，印證了年降水量對于大灣區(qū)地表水體空間分異影響較大。其次，年均溫（X3）與植被覆蓋率（X5）為67%、年均溫（X3）與高程（X1）為65%、土壤類型（X6）與植被覆蓋率（X5）為64%，都呈現(xiàn)較高的解釋力，說明植被覆蓋率、年均溫、高程、土壤類型也是影響地表水體空間分異的重要因素。

圖3 因子交互探測熱力矩陣Fig.3 The thermal matrix of factor interaction detection

在社會經(jīng)濟因子中，建設(shè)用地密度（X10）、人口因子（X8）在與其他因子交互過程中多呈現(xiàn)非線性增強，在一定程度上提高交互因子對地表水體時空分異的影響力，其中，建設(shè)用地密度與人口交互達到38%的解釋力，耕地密度和人口達到32%的解釋力。這說明建設(shè)用地密度、人口、耕地密度因子對于大灣區(qū)地表水體空間分異也產(chǎn)生一定程度的影響。

2.4.3 地表水體演變驅(qū)動機制淺析 大灣區(qū)水體緩慢萎縮階段（1985-2000 年）：該階段灣區(qū)各城市在改革開放的大背景下陸續(xù)進入快速城市化階段（圖4-a），第二、三產(chǎn)業(yè)發(fā)展迅猛，人口聚集程度提高，城市快速擴張，帶來的用地需求不斷增長，因此，通過填湖、填海等方式獲取大量建設(shè)用地以緩解人地矛盾。該階段大灣區(qū)共421.73 km2水域轉(zhuǎn)為耕地，536 km2轉(zhuǎn)為建設(shè)用地，城市建成區(qū)面積共計增加456.68 km2，約為原面積的2.5 倍（圖4-b），導(dǎo)致河網(wǎng)水系逐漸萎縮。

圖4 1985—2020年粵港澳大灣區(qū)各地市城市化率（a）與城市建成區(qū)面積變化（b）Fig.4 Urbanization rate(a) and urban built-up area(b) of cities in the GBA from 1985 to 2020

大灣區(qū)水體加速萎縮階段（2000-2010 年）：隨著城市化率進一步提升，城市規(guī)模進一步擴大。該階段大灣區(qū)共計631.23 km2水域轉(zhuǎn)為建設(shè)用地，耕地也同樣因城市開發(fā)被大量侵占（圖5-a），約1 700 km2耕地轉(zhuǎn)為建設(shè)用地，城市建成區(qū)面積增加至3 202.48 km2，約為原面積的4.12倍。城市對水域的過度侵占使灣區(qū)地表水體加速萎縮，帶來流域下墊面改變、水質(zhì)下降、水土流失等問題；至2010年，廣東省全域僅7個水庫水功能區(qū)水質(zhì)達標，62個水功能區(qū)監(jiān)測未達標（廣東省水利廳，2022a）。

圖5 往年廣東省耕地面積、工業(yè)企業(yè)（a）及水功能區(qū)（b）統(tǒng)計數(shù)據(jù)Fig.5 Arable land area, industry(a) and water function zone(b)statistics in Guangdong Province in previous years

大灣區(qū)水體萎縮減速放緩階段（2010-2015年）：2010 年大灣區(qū)內(nèi)已有6 個城市（包括香港、澳門）城市化率高于70%，在經(jīng)歷大規(guī)模城市擴張與內(nèi)部調(diào)整后，此時該類城市進入城市化穩(wěn)定階段，城市結(jié)構(gòu)趨于穩(wěn)定。本階段城市發(fā)展理念逐漸轉(zhuǎn)變，產(chǎn)業(yè)結(jié)構(gòu)也隨之調(diào)整，高能耗高污染產(chǎn)業(yè)逐漸被淘汰，政府部門陸續(xù)出臺一系列水體保護規(guī)劃與水體保持實施方案（表6），在一定程度上驅(qū)動大灣區(qū)地表水體面積降速放緩。

表6 水體保護與管控相關(guān)文件Table 6 List of documents related to water protection and management

大灣區(qū)水體穩(wěn)步增長階段（2015-2020 年）：2020年僅江門與肇慶的城市化率低于70%，其余城市均進入城市化穩(wěn)定階段。隨著城市產(chǎn)業(yè)結(jié)構(gòu)進一步轉(zhuǎn)型以及水體保護管理政策的不斷推進，各城市萬元GDP用水量指標穩(wěn)步下降，整體水質(zhì)得到明顯提升（圖5-b、圖6），水土流失面積逐漸減少（廣東省水利廳，2022b），河網(wǎng)管理制度逐漸完善，灣區(qū)地表水體面積呈現(xiàn)穩(wěn)步增長狀態(tài)。

綜上，城市化發(fā)展階段不同，自然因素和社會經(jīng)濟因素對于地表水體的驅(qū)動方向與方式各不相同，自然因素在早期城市化過程中對于驅(qū)動地表水體的空間分異起主導(dǎo)作用，但隨著城市化進程加速與城市不斷擴張，人類活動對于地表水體空間分異的影響力快速上升，逐漸成為驅(qū)動地表水體時空演變的重要因素；進入城市化穩(wěn)定階段，城市發(fā)展理念逐漸從經(jīng)濟導(dǎo)向轉(zhuǎn)為生態(tài)導(dǎo)向，在當(dāng)前中國“綠水青山就是金山銀山”生態(tài)發(fā)展理念指導(dǎo)下，對于地表水體保護力度與管理措施將不斷深化。

3 結(jié)論與討論

1）通過分析粵港澳大灣區(qū)地表水體面積變化速率，根據(jù)變化速率的階段性特征將研究時段劃分為地表水緩慢萎縮（1985-2000 年）、加速萎縮（2000-2010 年）、減速放緩（2010-2015 年）、穩(wěn)步增長（2015-2020 年）4 個階段。大灣區(qū)在水體萎縮階段地表水體面積共下降895.48 km2，降速逐漸提高；在減速放緩階段，雖然降速放緩，但面積仍下降約43.24 km2；穩(wěn)步增長階段地表水體面積轉(zhuǎn)降為增，面積增加約444.28 km2。

2）粵港澳大灣區(qū)內(nèi)永久性水體變化存在“小周期交替增減”特征，而季節(jié)性水體呈現(xiàn)“大幅度先減后增”特征。兩類水體對地表水總面積變化影響以2010年為時間節(jié)點呈“交替支配”特征，這是由于永久性水體面積變化幅度在各時間段內(nèi)逐漸低于季節(jié)性水體，而季節(jié)性水體面積在2010年后保持大幅度增加勢頭所造成的。

3）空間上，各階段地表水體的大幅度變化一般以城市組團的形式出現(xiàn)，如佛廣肇、深莞惠、珠中江等區(qū)域城市群，在長時序水體演變中呈同步加速萎縮或增長。地類轉(zhuǎn)換上，水域轉(zhuǎn)入轉(zhuǎn)出情況與經(jīng)濟發(fā)展理念和城市化進程緊密相關(guān)，在城市化初期為了獲取更多的建設(shè)用地用于城市發(fā)展，大量水域被建設(shè)用地、耕地侵占；隨著城市化率的提高，城市發(fā)展進入穩(wěn)定階段，耕地、林地轉(zhuǎn)入水域的面積不斷增多。

4）從研究期全時段看，自然因子相較社會經(jīng)濟因子對粵港澳大灣區(qū)地表水體的空間分異擁有更高解釋力，其中，年降水量和植被覆蓋率因子解釋力達到40%和34%，且與除建設(shè)用地密度和交通路網(wǎng)長度之外的所有因子均呈現(xiàn)雙因子增強的交互關(guān)系；社會經(jīng)濟因子中，耕地密度、人口因子解釋力為31%和30%，與多因子存在非線性增強的交互關(guān)系。

5）從城市發(fā)展階段看，自然因素和社會經(jīng)濟因素對于地表水體的影響程度在不同階段存在差異。自然因素在早期城市化過程中對地表水體時空分異起主導(dǎo)作用，隨著城市化發(fā)展進程加快并進入穩(wěn)定階段，人類活動對于地表水體的影響力快速上升，成為驅(qū)動地表水體時空分異的重要因素。

本文從長時序視角探究了1985-2020年粵港澳大灣區(qū)地表水體時空演變特征，同時量化了相關(guān)影響因素對于地表水體空間分異的解釋力。與現(xiàn)有研究的不同與創(chuàng)新之處在于，本研究從城市化發(fā)展進程的角度，結(jié)合城市發(fā)展政策與區(qū)域發(fā)展導(dǎo)向，解析灣區(qū)各階段地表水體演變的內(nèi)在規(guī)律與驅(qū)動機制。但仍存在不足之處，如受到數(shù)據(jù)年份限制，水體演變階段劃分過于粗略、地理探測因子選取不夠全面等問題，后續(xù)研究可以通過解譯高時空分辨率的遙感影像（如Sentinal-2影像）來獲取更高精度、更長時序的歷史水體數(shù)據(jù)。同時，就地表水體的管理與決策而言，后續(xù)可以通過相關(guān)模型算法（MCCA 模型、PLUS 模型等）模擬地表水體的動態(tài)演變，能更加有效的為地表水體的可持續(xù)利用提供科學(xué)依據(jù)。