基于GWR模型的東江水質空間分異與水生態功能分區驗證

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

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

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

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

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

1 材料與方法

1.1 研究區概況

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

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

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

1.2 水質數據采樣與測定

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

1.3 空間統計數據庫

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

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

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

表1 流域環境因素指標集及基本統計特征

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土壤質地Soil texture穩定下滲率/(mm/h)Steady infiltration rate滲透能力Infiltration abilityA砂土、壤質砂土、砂質壤土>7.26強B壤土、粉砂壤土3.81—7.26較強C砂粘壤土1.27—3.81中等D粘壤土、粉砂壤土、砂粘土、粉砂粘土、粘土0.00—1.27弱

SCS為美國農業部水土保持局(Soil Conservation Server,SCS)開發的流域水文模型

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

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

為保證樣點環境數據不受鄰近樣點位置的影響,分別以每個樣點為出水口提取其上游集水區[24],并基于每個樣點上游集水區范圍分別計算其土地利用和自然背景指標,共計提取612個圖層進行分層統計。

1.4 統計方法

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

1.4.1普通線性回歸模型OLS

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

(1)

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

1.4.2地理加權回歸模型GWR

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

(2)

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

(3)

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

(4)

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

1.4.3OLS模型與GWR模型評價

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

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

(5)

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

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

2 結果與分析

2.1 水生態功能分區水質差異特征

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

圖4 一二級水生態區水質指標差異(均值±標準誤)Fig.4 The differences of water quality parameters among ecoregions

2.2 基于GWR模型評估水質與影響因素關系

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

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

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

2.3 OLS和GWR模型性能評價

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

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

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

圖6 一二級水生態區GWR模型local R2差異(均值±標準誤)圖Fig.6 The differences of local R2 from GWR models among ecoregions

水質參數Water quality parameters穩態評估K(BP)-Prob校正R2 Adjusted R2殘差空間自相關指數 Moran′s IOLSGWROLSGWR水溫 LnTEM0.1790.3670.3720.140.12溶解氧 LnDO0.0760.5060.498-0.21-0.22電導率 LnEC0.040?0.6090.6890.310.00高錳酸鹽指數 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

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

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

3 結論

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