基于GWR模型的東江水質(zhì)空間分異與水生態(tài)功能分區(qū)驗證

和克儉,黃曉霞,*,丁 佼,劉 琦,江 源

1 云南大學(xué)資源環(huán)境與地球科學(xué)學(xué)院,昆明 650091 2 云南省水利廳,昆明 650021 3 北京師范大學(xué)地理科學(xué)學(xué)部,北京 100875

水質(zhì)表征水環(huán)境質(zhì)量狀態(tài),水質(zhì)變化是流域地貌、水文、生物過程和人類活動等多因素綜合影響的結(jié)果[1]。水生態(tài)區(qū)是淡水生態(tài)系統(tǒng)或生物體及其與環(huán)境相互關(guān)系相對一致的土地單元[2- 3]。已有研究指出以水生態(tài)區(qū)為基礎(chǔ)的水污染防治體系,比單純以工程技術(shù)為主的污染防治體系更有前景[4]。2008年我國啟動國家水體污染控制與治理重大科技專項課題(水專項),系統(tǒng)地開展了流域水生態(tài)功能分區(qū)研究工作[5- 9]。水生態(tài)功能分區(qū)強調(diào)基于水陸生態(tài)系統(tǒng)的耦合關(guān)系研究,通過辨識影響水生態(tài)系統(tǒng)結(jié)構(gòu)與功能的關(guān)鍵陸域或水環(huán)境因子,將這些因子的空間差異作為分區(qū)指標(biāo),揭示流域水生態(tài)系統(tǒng)格局與功能的空間異質(zhì)性特征[10- 13]。盡管水生態(tài)功能分區(qū)以生態(tài)學(xué)理論為基礎(chǔ),但分區(qū)結(jié)果的科學(xué)性和合理性還有待驗證,如何驗證是當(dāng)前亟待解決的問題。

目前,流域特征與地表水質(zhì)關(guān)系的定量研究主要基于全局統(tǒng)計方法,比如相關(guān)分析、主成分分析、冗余分析、普通線性回歸等。全局統(tǒng)計方法一般假定水質(zhì)指標(biāo)與流域影響因素之間的關(guān)系在研究區(qū)域內(nèi)穩(wěn)定,不隨空間位置改變而變化[14]。而實際上水質(zhì)與流域影響因素之間的關(guān)系可能存在局部差異,具有空間非平穩(wěn)性。比如流域內(nèi)不同區(qū)域,由于主導(dǎo)的土地利用方式不同,可能導(dǎo)致影響區(qū)域水質(zhì)的關(guān)鍵因素發(fā)生變化。地理加權(quán)回歸模型(Geographically weighted regression,GWR)[15]是一種局域空間分析方法,該方法將數(shù)據(jù)的地理位置嵌入回歸模型中,使回歸系數(shù)成為與空間位置有關(guān)的函數(shù),其分析結(jié)果是在每個樣點都會得到一組局部的參數(shù)估計。因此,GWR模型對于空間差異化關(guān)系具有強大的分析能力。作為一門新興技術(shù),近年來GWR模型逐漸成為土地利用對水質(zhì)影響研究領(lǐng)域的新熱點[16- 20],但目前相關(guān)研究側(cè)重于定量描述GWR結(jié)果,少有研究將GWR結(jié)果與水生態(tài)區(qū)劃結(jié)果相結(jié)合來探討分區(qū)合理性。

東江流域是我國流域水生態(tài)功能分區(qū)研究的重點流域之一,目前已完成的一二級水生態(tài)功能分區(qū)結(jié)果主要表達了溫度、降水、地形、地表覆被等大尺度環(huán)境因子影響下的區(qū)域水資源量和水質(zhì)空間異質(zhì)性特征[21]。本文以東江流域為例,基于大范圍采樣數(shù)據(jù)和GWR模型,探討水質(zhì)空間規(guī)律及其與流域因素的空間關(guān)系與水生態(tài)功能分區(qū)結(jié)果是否一致,以檢驗分區(qū)結(jié)果的合理性,為分區(qū)管理和水生態(tài)系統(tǒng)保護提供科學(xué)依據(jù)；并對比GWR模型與普通最小二乘(OLS)模型的性能,探討GWR模型在分區(qū)結(jié)果驗證中的應(yīng)用價值與不足,為其他流域水生態(tài)功能分區(qū)研究提供借鑒。

1 材料與方法

1.1 研究區(qū)概況

東江(113°04′—115°50′ E,22°21′—25°12′ N)是珠江三大水系之一,發(fā)源于江西省尋鄔縣椏髻缽山,自東北向西南流經(jīng)廣東省,于東莞注入獅子洋(圖1)。東江流域總面積35340 km2,干流全長562 km,地處亞熱帶季風(fēng)氣候區(qū),多年平均氣溫21°C,多年平均降水量1800 mm,雨季(4—9月)降水占全年降水的80%以上[22]。

圖1 東江流域采樣點及水生態(tài)功能分區(qū)Fig.1 Sampling sites and aquatic ecoregions of the Dongjiang River basinR：河流型流域 River；F：東江流域編碼；"I、II、III…"為一級分區(qū)編碼；下標(biāo)"1、2、3…"為二級分區(qū)編碼

流域內(nèi)主要土地利用類型為林地、農(nóng)業(yè)用地(包括耕地和園地)和城鎮(zhèn)用地,林地主要分布在北部區(qū)域,城鎮(zhèn)用地和農(nóng)業(yè)用地主要分布在中部和南部；流域地勢北高南低,北部多為中、低山地,中部以低山丘陵為主,南部為沿江平原和河口三角洲；流域內(nèi)坡地(坡度>7°)主要分布在北部和中部,低平地(坡度<7°)主要分布在南部；流域內(nèi)土壤以赤紅壤、紅壤和水稻土為主(圖2)。基于自然背景區(qū)域差異,東江流域被劃分為3個水生態(tài)功能一級區(qū)和9個水生態(tài)功能二級分區(qū)(圖1)[21]。

1.2 水質(zhì)數(shù)據(jù)采樣與測定

采用快照采樣(Snapshot sampling)方法于2010年7月穩(wěn)定流量狀況(48小時內(nèi)降雨小于10 mm)下共采集102個樣點水樣[22],樣點覆蓋干流及全流域主要支流(圖1)。實驗設(shè)計基于如下考慮：1)2010年東江流域氣候特征與流域長期氣候特征一致,降雨量為1787 mm,年月均溫為21.3°,雨季降雨量占全年降雨量的82%[23]；2)前期研究表明營養(yǎng)鹽和有機污染等指標(biāo)的季節(jié)差異不明顯[24]；3)大范圍空間采樣能為揭示水質(zhì)空間差異提供穩(wěn)健結(jié)果[25-26]。

1.3 空間統(tǒng)計數(shù)據(jù)庫

基于水質(zhì)影響因素研究相關(guān)成果[18,28- 30],本文從土地利用和自然背景兩方面構(gòu)建流域影響因素指標(biāo)集(表1)。

土地利用數(shù)據(jù)基于2009年Landsat TM/ETM影像(分辨率 30 m)進行解譯[31],基于流域特征及研究目的,將土地利用類型重分類為耕地、園地、水域、林地、草地、城鎮(zhèn)用地及其他用地類型(圖2)。土地利用指標(biāo)(斑塊面積百分比PLAND、斑塊密度PD)用FRAGSTATS 4.0軟件[32]進行計算。海拔(ELEVA)和坡度(SLOPE)指標(biāo)基于分辨率30 m的ASTER GDEM(http://glovis.usgs.gov/)數(shù)據(jù)進行計算。

土壤類型數(shù)據(jù)使用東江流域1∶100萬的土壤類型圖(中國科學(xué)院南京土壤研究所)。根據(jù)美國農(nóng)業(yè)部水土保持局(Soil Conservation Server,SCS)開發(fā)的流域水文模型土壤水文組(hydrologic soil group,HSG)的劃分標(biāo)準(zhǔn)[33],參考東江流域土壤的水文物理特征[34],按照產(chǎn)流能力從低到高將土壤類型分為4組(表 2),水域不計入土壤水文組類別(NA)(圖2)。

表1 流域環(huán)境因素指標(biāo)集及基本統(tǒng)計特征

PLAND：斑塊面積百分比 Percentage of land use；PD：斑塊密度 Patch density；ELEVA：海拔 Elevation；SLOPE：坡度 Slope；HSG：土壤水文組 Hydrologic soil group；SOLAR：太陽輻射量 Solar radiation；RAIN：降雨量 Precipitation

圖2 東江流域土地利用、海拔、坡度、土壤水文組、2010年流域雨季前期太陽輻射量、2010年流域雨季前期降雨量空間分布圖Fig.2 The spatial distribution of land use, elevation, slope, hydrological soil group, solar radiation and rainfall during the early stage of rainy season in 2010 in the Dongjiang River basinHSGA：土壤水文組A面積百分比 Hydrologic soil group A；HSGB：土壤水文組B面積百分比 Hydrologic soil group B；HSGC：土壤水文組C面積百分比 Hydrologic soil group C；HSGD：土壤水文組D面積百分比 Hydrologic soil group D

土壤水文組Hydrologic soil group土壤質(zhì)地Soil texture穩(wěn)定下滲率/(mm/h)Steady infiltration rate滲透能力Infiltration abilityA砂土、壤質(zhì)砂土、砂質(zhì)壤土>7.26強B壤土、粉砂壤土3.81—7.26較強C砂粘壤土1.27—3.81中等D粘壤土、粉砂壤土、砂粘土、粉砂粘土、粘土0.00—1.27弱

SCS為美國農(nóng)業(yè)部水土保持局(Soil Conservation Server,SCS)開發(fā)的流域水文模型

由于夏季樣點間氣溫差異小,因此本研究舍棄氣溫指標(biāo),選取太陽輻射量以表征熱量空間差異(圖2)。太陽輻射的觀測密度通常較小,采用簡單的空間插值技術(shù)不能合理地揭示太陽輻射空間分布特征[35]。因此本研究采用ArcGIS 9.3軟件中Solar radiation 工具模擬雨季開始至采樣前(2010年4月—6月)累積的太陽總輻射量[36](天氣狀況：晴朗；透射率：0.5,散射比例：0.3)。并根據(jù)相鄰氣象站輻射日值觀測數(shù)據(jù),選取每5天中的最大日總輻射量作為晴朗天氣下太陽總輻射量,驗證模型模擬的精度,驗證結(jié)果表明模型模擬效果較好(圖3)。2010年4月—6月降雨量數(shù)據(jù)(http://cdc.cma.gov.cn/)以氣象站高程作為協(xié)方差進行cokriging空間插值(圖2)。

圖3 日太陽總輻射量實測值與模擬值比較 Fig.3 The comparison between measured and simulated daily solar radiation

為保證樣點環(huán)境數(shù)據(jù)不受鄰近樣點位置的影響,分別以每個樣點為出水口提取其上游集水區(qū)[24],并基于每個樣點上游集水區(qū)范圍分別計算其土地利用和自然背景指標(biāo),共計提取612個圖層進行分層統(tǒng)計。

1.4 統(tǒng)計方法

水質(zhì)區(qū)域差異分析采用Kruskal-Wallis非參數(shù)方差分析,并進行Dunn多重比較(posthoctest)分析,顯著性水平為P<0.05。方差分析在SPSS 21.0軟件中實現(xiàn)；均值圖在Origin 9.0軟件中實現(xiàn)。

1.4.1普通線性回歸模型OLS

采用普通最小二乘法(Ordinary Least Squares,OLS)線性回歸模型,評估單個水質(zhì)指標(biāo)與流域影響變量之間的全局關(guān)系,計算公式[32]如下：

(1)

式中,y為因變量；xi為解釋變量；p為預(yù)測變量個數(shù)；β0為截距；βi為回歸系數(shù)；ε為均值為0、方差為σ2的誤差項。水質(zhì)指標(biāo)預(yù)先進行Ln(X+1)轉(zhuǎn)換,解釋變量采用Z標(biāo)準(zhǔn)化。采用逐步多元回歸挑選顯著的解釋變量進入模型。根據(jù)如下標(biāo)準(zhǔn)挑選各水質(zhì)指標(biāo)的最佳OLS模型：1)模型校正R2(adjustedR2)最高；2)模型參數(shù)及解釋變量參數(shù)顯著(P≤0.05)；3)解釋變量方差膨脹因子(Variance Inflation Factor,VIF)小于3。以上操作在SPSS 21.0實現(xiàn)。

1.4.2地理加權(quán)回歸模型GWR

地理加權(quán)回歸模型(Geographically weighted regression,GWR)是對普通線性回歸模型的擴展,其公式[30]如下：

(2)

式中,(uj,vj)為第j個采樣點的坐標(biāo)(公里網(wǎng)坐標(biāo))；βi(uj,vj)為第j個采樣點上的第i個回歸系數(shù),即地理位置的函數(shù)；εj為誤差項。根據(jù)加權(quán)最小二乘法,βi(uj,vj)可通過使

(3)

達到最小來進行估計,式中,wjk為回歸點j與其他已知觀測點k之間的距離衰減函數(shù),其基本假設(shè)為距離j點越近的觀測點對求解局部回歸系數(shù)的重要性越大,越遠的觀測點重要性越小[13]。wjk為空間權(quán)重矩陣,其計算方法包括距離閾值法、距離反比法、Gauss函數(shù)法及bi-square函數(shù)法等[14- 15]。本研究選用Gauss函數(shù)計算空間權(quán)重矩陣,公式[37]如下：

(4)

式中,djk為觀測點j、k之間距離；b為帶寬。本研究以AICC(修正的Akaike信息準(zhǔn)則)為衡量標(biāo)準(zhǔn)的自適應(yīng)帶寬形式,識別最佳自適應(yīng)鄰近點個數(shù),即在數(shù)據(jù)密集地方采用較小的帶寬而在數(shù)據(jù)稀疏的地方采用較大的帶寬。為避免數(shù)據(jù)共線性影響,采用最佳逐步OLS模型挑選出的顯著解釋變量輸入GWR模型。

1.4.3OLS模型與GWR模型評價

采用穩(wěn)態(tài)評估Koenker (BP) 統(tǒng)計量[31]對OLS模型進行空間平穩(wěn)性檢驗,若K(BP)-Prob<0.05表示模型具有顯著的非穩(wěn)態(tài),即模型變量間關(guān)系存在空間差異,說明模型適合進行地理加權(quán)回歸(GWR)分析。

通過比較OLS和GWR模型的adjustedR2和殘差空間自相關(guān)指數(shù)(Moran′s I),分析GWR模型在預(yù)測精度和處理空間自相關(guān)過程中是否優(yōu)于OLS模型。采用空間自相關(guān)工具計算模型標(biāo)準(zhǔn)化回歸殘差的全局 Moran′s I 指數(shù),其計算公式[26]如下：

(5)

式中,Xi、Xj分別為要素i和j的值；n為要素個數(shù)；wij為要素i和j之間空間權(quán)重,定義為兩者之間距離的倒數(shù)。如果 Moran′s I 指數(shù)值為正則指示聚集趨勢,其值為負則指示離散趨勢,其值接近零則指示隨機模式。

GWR模型估計以及空間自相關(guān)分析均在ArcGIS 9.3軟件中進行。

2 結(jié)果與分析

2.1 水生態(tài)功能分區(qū)水質(zhì)差異特征

圖4的Kruskal-Wallis檢驗結(jié)果顯示：各水質(zhì)指標(biāo)(除NO3-N外)在水生態(tài)一級區(qū)RFⅠ與RFⅢ、RFⅡ與RFⅢ之間差異顯著(P<0.05)；DO、CODMn、TN和NO3-N在RFⅠ內(nèi)二級區(qū)之間差異顯著(P<0.05)；TEM、EC、TN和NO3-N在RFⅡ內(nèi)二級區(qū)之間差異顯著(P<0.05)；EC、TP、TN、NO3-N和Chl-a在RFⅢ內(nèi)二級區(qū)之間差異顯著(P<0.05)。說明一二級水生態(tài)功能分區(qū)結(jié)果能較好地反映流域水質(zhì)空間差異。此外,TN在一級區(qū)之間,以及一級區(qū)內(nèi)二級區(qū)之間,均存在顯著差異,這表明TN是理想的分區(qū)驗證指標(biāo)。

圖4 一二級水生態(tài)區(qū)水質(zhì)指標(biāo)差異(均值±標(biāo)準(zhǔn)誤)Fig.4 The differences of water quality parameters among ecoregions

2.2 基于GWR模型評估水質(zhì)與影響因素關(guān)系

通過逐步多元回歸篩選最優(yōu)OLS模型(表 3),結(jié)果顯示除了TEM和Chl-a模型解釋率較低外(adjusted

LnTEM：水溫 Water temperature；LnDO：溶解氧 Dissolved oxygen；LnEC：電導(dǎo)率 Electrical conductivity；LnCODMn：高錳酸鹽指數(shù) Permanganate index；LnTP：總磷 Total phosphorus；LnTN：總氮 Total nitrogen；LnNH3-N：氨氮 Ammonia nitrogen；LnNO3-N：硝氮 Nitrate nitrogen；LnChl-a：葉綠素a Chlorophyll a

R2<0.4),其他指標(biāo)模型解釋率均較高(adjustedR2>0.4),說明流域影響因素能解釋大部分水質(zhì)差異。GWR模型局部解釋率(LocalR2)空間差異說明流域環(huán)境因素對水質(zhì)的影響存在明顯的空間差異,且這種空間差異與水生態(tài)功能分區(qū)結(jié)果具有高度一致性(圖5和圖6)。各水質(zhì)指標(biāo)的Local R2在三個水生態(tài)一級區(qū)(RFⅠ、RFⅡ和RFⅢ)之間差異顯著(P<0.05)；除EC在RFⅠ1和RFⅠ2之間差異不顯著外,各指標(biāo)Local R2在一級區(qū)內(nèi)二級區(qū)間也存在顯著差異(P<0.05)。這說明基于流域特征的水生態(tài)功能分區(qū)結(jié)果與流域環(huán)境因素對水質(zhì)影響關(guān)系的空間差異狀況一致。采用GWR模型評估流域特征與水質(zhì)關(guān)系的空間差異,能簡單有效地驗證水生態(tài)功能分區(qū)的合理性。

2.3 OLS和GWR模型性能評價

Koenker(BP)穩(wěn)態(tài)檢驗結(jié)果(表 4)表明EC、TP、NH3-N、NO3-N、Chl-a的OLS回歸模型存在空間非穩(wěn)態(tài)(P<0.05),說明這些水質(zhì)指標(biāo)更適合運行GWR模型。通過比較模型adjustedR2和Moran′s I值,對于營養(yǎng)鹽和有機污染指標(biāo)而言,GWR模型比OLS模型具有更高的解釋率,殘差空間自相關(guān)性更低。表明GWR模型能改善OLS模型殘差空間依賴性問題,具有更好的擬合效果。

值得注意的是,DO不僅OLS模型具有空間平穩(wěn)性(P>0.05),其OLS的adjustedR2高于GWR模型,且OLS殘差自相關(guān)程度弱于GWR模型,因此更適合進行OLS模型分析。此外,盡管TP模型存在顯著的空間非平穩(wěn)性(P<0.05),但其OLS模型和GWR模型回歸殘差均呈聚集狀態(tài),說明該模型缺失了某個關(guān)鍵解釋變量,有可能是未考慮人為點源排放的影響。因此,DO和TP可能不適合作為分區(qū)驗證指標(biāo)。

圖5 地理加權(quán)回歸(GWR)模型的local R2空間分布圖Fig.5 The spatial distribution of local R2 from Geographically Weighted Regression models

圖6 一二級水生態(tài)區(qū)GWR模型local R2差異(均值±標(biāo)準(zhǔn)誤)圖Fig.6 The differences of local R2 from GWR models among ecoregions

水質(zhì)參數(shù)Water quality parameters穩(wěn)態(tài)評估K(BP)-Prob校正R2 Adjusted R2殘差空間自相關(guān)指數(shù) Moran′s IOLSGWROLSGWR水溫 LnTEM0.1790.3670.3720.140.12溶解氧 LnDO0.0760.5060.498-0.21-0.22電導(dǎo)率 LnEC0.040?0.6090.6890.310.00高錳酸鹽指數(shù) LnCODMn0.0900.5160.5210.140.12總磷 LnTP<0.001??0.5890.6011.431.43總氮 LnTN0.1570.6530.6650.070.05氨氮 LnNH3-N<0.001??0.6490.6530.200.18硝氮 LnNO3-N0.019?0.4030.4500.130.00葉綠素a LnChl-a0.004??0.2200.2170.460.44

P<0.05； **P<0.01

隨著區(qū)劃理論和技術(shù)的發(fā)展,未來的生態(tài)區(qū)劃將更多的從自然、社會經(jīng)濟等不同的影響因素出發(fā)進行指標(biāo)選取,并且指標(biāo)選取的依據(jù)將更加量化。GWR模型作為一種空間明確的分析方法,不僅估計結(jié)果有明確的解析表示,而且得到的參數(shù)估計還能進行統(tǒng)計檢驗,在空間數(shù)據(jù)分析和變量關(guān)系的空間制圖方面具有明顯的優(yōu)勢,能提供傳統(tǒng)全局統(tǒng)計方法無法給出的信息。因此GWR方法能夠更好的用于識別區(qū)劃目標(biāo)要素與自然、社會經(jīng)濟因子之間耦合關(guān)系的空間非平穩(wěn)性,為分區(qū)指標(biāo)的選取、分區(qū)結(jié)果的驗證、不同區(qū)劃方案的直接比較等方面提供新的思路和方法。

然而GWR模型尚存在一些不足之處,比如GWR模型雖然在一定程度上考慮了數(shù)據(jù)的空間自相關(guān),但它并沒有明確地解決空間自相關(guān)問題；GWR模型基于歐式距離測度距離衰減效應(yīng),而對于河流系統(tǒng)而言,還應(yīng)考慮沿徑流網(wǎng)絡(luò)的距離衰減效應(yīng),因此基于歐式距離的測度方式不能反映水質(zhì)采樣點間真實的空間接近程度。如何降低數(shù)據(jù)空間自相關(guān)帶來的影響,以及使用更合適的距離測度方法,是未來GWR模型應(yīng)用于水生態(tài)系統(tǒng)研究的難點問題。

3 結(jié)論

東江流域各水質(zhì)指標(biāo)(除NO3-N外)在水生態(tài)一級區(qū)之間差異顯著；DO、CODMn、TN和NO3-N在RFⅠ內(nèi)二級區(qū)之間差異顯著；TEM、EC、TN和NO3-N在RFⅡ內(nèi)二級區(qū)之間差異顯著；EC、TP、TN、NO3-N和Chl-a在RFⅢ內(nèi)二級區(qū)之間差異顯著。說明東江水生態(tài)功能分區(qū)結(jié)果能夠較好地反映流域水質(zhì)空間差異?；贕WR模型,各水質(zhì)指標(biāo)的Local R2在水生態(tài)一級區(qū)之間差異顯著；除EC外,各指標(biāo)Local R2在一級區(qū)內(nèi)二級區(qū)間也存在顯著差異,說明流域環(huán)境因素對水質(zhì)的影響存在明顯的空間差異,且這種空間差異與水生態(tài)功能分區(qū)結(jié)果具有高度一致性,表明東江水生態(tài)功能分區(qū)結(jié)果能合理反映水陸耦合關(guān)系。此外,TN在一級區(qū)之間,以及一級區(qū)內(nèi)二級區(qū)之間均存在顯著差異,可以作為理想的分區(qū)驗證指標(biāo)。而DO由于模型變量關(guān)系表現(xiàn)出空間平穩(wěn)性,TP由于OLS模型和GWR模型回歸殘差均呈聚集狀態(tài),不適合作為分區(qū)驗證指標(biāo)。GWR模型相比OLS模型具有更好的模型擬合效果,對空間數(shù)據(jù)分析和空間制圖具有明顯的優(yōu)勢,在分區(qū)研究中具有廣泛應(yīng)用前景。降低數(shù)據(jù)空間自相關(guān)影響及改善距離測度方法是未來GWR模型研究的難點問題。