北運河北京段水生態承載力研究及關鍵控制要素識別

杜 芳, 彭劍峰, 王亞俊, 張 旺

1.中國礦業大學， 江蘇 徐州 221116

2.清華大學環境學院水質與水生態研究中心， 北京 100084

3.中國礦業大學(北京)， 北京 100083

北京市副中心承擔著疏解非首都功能、推進京津冀協同發展的功能，區域內主要水系屬于北運河水系[1]. 北運河北京段承載了北京市主要人口，多年來水質狀況較差，水生態現狀不容樂觀[2-3]. 隨著北京市副中心的快速建設和發展，區域內水資源-水環境-水生態面臨嚴峻的承載壓力. 生態系統作為城市生態系統的承載基底，水生態承載力制約著城市發展的速度和規模. 如何解析北運河北京段水生態承載力演變過程，并識別關鍵控制要素，是提升副中心地區生態承載力的前提要求[4-5].

水生態承載力(WECC)是近年新興的，基于承載力理論，綜合水資源、水環境等要素對水生態系統進行研究的一個綜合性的概念，目前學術界對其定義及研究方法仍在不斷探索中，尚未達成統一共識[6]. 水生態承載力指標體系應綜合水生態系統的水資源、水環境等要素. 國內外學者多采用供需平衡法、生態足跡法、綜合指標法等方法評估流域水生態承載狀況[7-8]，其中，綜合指標評價法可綜合考慮各類因素，結構體系清晰，便于厘清各因素間關系及確定關鍵控制要素[9]. 水生態系統是城市生態系統的重要組成部分，同時水生態承載力制約著城市發展的速度和規模. 北京市副中心所屬北運河北京段，承擔著北京市排水、灌溉等多項功能，當前關于水生態承載力的相關研究主要集中于水資源承載力、水環境容量以及水生生物結構功能完整性等內容[10-11]. 現有研究中仍缺乏綜合水資源、水環境、水生態指標長時間序列下的水生態承載力評估和關鍵要素識別[12-13].

該研究通過綜合指標法構建北運河北京段的水生態承載力評估體系，分別計算了1995年、2000年、2005年、2010年、2015年和2018年北運河北京段區域整體和子流域的水生態承載力，分析了近30年水生態承載力演變特征，并采用多層感知器神經網絡模型識別水生態承載力的關鍵控制要素，以期為副中心生態承載力提升和壓力減負提供科學參考.

1 材料與方法

1.1 研究區概況

北運河發源于北京市軍都山南麓，自西北向東南流入通州區，后流經河北省香河縣、天津市武清區，并在天津市大紅橋匯入海河. 北運河上游以山區丘陵為主，中下游為沖積平原，北運河北京段區域面積約 4 348 km2，主河道長89.4 km. 北運河流域屬于暖溫帶大陸性季風氣候，降水集中于6—9月，多年平均降水量為581.7 mm，多年平均天然徑流量為4.8×108m3，年內來水量變化明顯，非汛期以再生水排水為主. 北運河北京段包括沙河、藺溝河、清河、通惠河、溫榆河、涼水河、壩河、小中河、北運河、鳳河10個子流域(見圖1)，承載了北京市70%以上的人口，該區域經濟總量占北京市的80%以上，容納了北京市90%以上的污水排放量. 根據《北京市水資源公報(2019年)》，2019年北運河水質以GB 3838—2002《地表水環境質量標準》Ⅳ類和Ⅴ類為主，水環境質量較差. 此外，就流域水生態狀況而言，京津冀對北運河流域水資源利用率較高，流域季節性斷流頻發，河流自凈能力差，生態基流難以保障.

圖1 北運河北京段子流域示意

1.2 研究方法

1.2.1水生態承載力評估體系構建及權重確定

該研究從水資源、水環境、水生態3個層面收集了30多個指標，并進行了因子有效性分析[14]. 在SPSS 19.0軟件中對指標進行主成分分析并進行最大方差法因子旋轉，指標體系通過了KMO度量、Bartlett球形度檢驗(KMO>0.6，Sig.=0.000). 在此基礎上劃分出水資源稟賦、水資源利用、水環境納污、水環境凈化、水生生境、水生生物6個分指標，人均水資源量、年降雨量、水資源開發利用率、農業污染強度指數、工業和生活污染強度指數、水環境質量指數、河岸帶林草覆蓋率、生境連通性指數、生境整體性指數、生境聚集度指數和葉綠素含量11個子指標(見表1)，構建了北運河北京段的水生態承載力評估體系.

為避免人為賦值帶來的主觀性影響，對各項指標進行標準化處理后采用熵權法確定指標權重[15-16]. 熵權法計算權重過程中，如果指標信息熵越大，則該指標提供的信息越充分，權重賦值越大. 正向指標和負向指標標準化的計算分別如式(1)(2)所示，第i個樣點第j項指標權重的計算公式如式(3)～(5)所示.

(1)

(2)

(3)

(4)

(5)

通過指標標準化和權重的確定，北運河北京段水生態承載力評估體系及指標權重如表1所示.

表1 北運河北京段水生態承載力評估體系

1.2.2水生態承載力分值計算

水生態承載力指標體系中水資源、水環境和水生態等3個總指標的計算如式(6)所示. 水生態承載狀態分值T值由總指標加權獲得，計算公式如式(7)所示.

(6)

T=A×WA+B×WB+C×WC

(7)

式中，Z為水生態承載力承載總指標，Ij為分指標第j個子指標指標值，Wj為第j個子指標權重，Wk為該總指標第k個分指標權重，WA、WB、WC分別為水資源、水環境、水生態總指標，A、B、C分別為總指標WA、WB、WC的權重.

根據水生態承載力狀態分值進行承載力狀態等級劃分，結果如表2所示.

表2 水生態承載力等級劃分

1.2.3水生態承載力差異性識別

單因素方差分析僅考慮一個因素的影響而固定其他因素，也可用于檢驗樣本量為兩個以上時均數差別是否具有顯著性[17-18]. 為了對子流域間水生態承載力的差異性進行識別，首先應對各子流域近30年的水生態承載力綜合分值進行方差齊性檢驗和正態分布檢驗(Shapiro-Wilk檢驗和K-S檢驗). 在方差齊性檢驗和正態分布檢驗均通過后進行單因素方差分析，這一過程在R 4.0軟件中運用multcomp包完成.

1.2.4區域承載力關鍵控制要素識別方法

為識別承載力關鍵控制要素，在IBM SPSS 19.0軟件中對各項指標進行正態分布檢驗，檢驗方法為K-S(Kolmogorov-Smirnov)法，對滿足正態分布和不滿足正態分布的指標分別用Pearson相關系數和Spearman相關系數進行分析[19]. 然后在IBM SPSS 19.0軟件的神經網絡模塊中按照訓練集和測試集7∶3 比例建模，對各子流域1995年、2000年、2005年、2010年、2015年和2018年共6期的水生態承載力綜合指標和子指標進行多層感知器神經網絡模型分析，并通過標準化后的重要性分析得出區域承載力關鍵控制要素[20]. 綜合分析相關系數和神經網絡模型重要性，最終得出區域承載力關鍵控制要素.

1.3 數據來源

1.3.1指標數據來源

北運河北京段區域1995年、2000年、2005年、2010年、2015年和2018年6期各子指標量化方法和數據來源如表3所示. 其中，水資源相關數據來自1995—2018年《北京市統計年鑒》《北京市通州區統計年鑒》和《北京市水資源公報》. Landsat TM/ETM+/OLI遙感影像(分辨率為30 m)、土地利用類型數據、相應年份的年降雨量、人口、GDP等空間柵格數據(分辨率為1 km)來自中國科學院資源環境科學與數據中心(http://www.resdc.cn/Default.aspx)[21-22]. 水生生境指數中的河岸帶林草覆蓋率由ArcGIS(500 m)緩沖區分析得到，生境連通性指數、生境整體性指數、生境聚集度指數分別由Frgstats 4.2軟件中的CONNECT指標、COHESION指標和AI指標計算得到. 生境連通性指數用來表示景觀間的連接狀況，景觀間連接狀況越好則生境連通性越強；生境整體性指數表示景觀斑塊間的整體性，景觀間整體性越強則生境整體性指數值越高；生境聚集度指數表示景觀聚集程度，即景觀中不同斑塊類型相鄰的概率；生境聚集度指數隨著景觀聚集程度的增加而增大[23-24].

表3 北運河北京段水生態承載力數據來源及量化方法

生境連通性指數(CONNECT)、生境整體性指數(COHESION，%)和生境聚集度指數(AI，%)的計算公式分別如式(8)(9)(10)所示：

(8)

(9)

(10)

式中：dijk為一定臨界距離內，與斑塊類型i相關的斑塊類型j、k的連接狀況；ni為斑塊類型i的斑塊個數；pij為斑塊ij的周長，m；aij為斑塊ij的面積，hm2；N為柵格總個數；eii為基于單倍法的斑塊類型i像元之間的結點數；maxeii為基于單倍法的斑塊類型i像元之間的最大結點數；CONNECT、COHESION、AI分別為生境連通性指數、生境整體性指數和生境聚集度指數，取值范圍均為[0，100],除CONNECT無單位外，其余兩個單位均為%.

遙感水質反演技術較常規實地采樣水質監測方式在監測面積和監測時長上均具有顯著優勢[25]. 傳統的實地采樣水質監測方式的不足之處主要有三方面：一是所需人力成本和物質財力成本較高；二是采樣所得數據在時間范圍和空間范圍上不連續，無法滿足長時間序列下，區域大空間范圍內持續監測的實際需求；三是因特殊原因在部分年份或部分采樣點數據缺失時，單靠傳統采樣方式難以彌補部分數據缺失問題. 因此，傳統監測方法難以滿足流域管理與決策的實際需求. 近年來，為解決傳統采樣獲取數據時空連續性不足、部分區域水質歷史數據缺乏等實際問題，遙感水質反演技術被廣泛運用于內陸水質監測[26-27]. 由于北運河水質數據來源限制，部分年份存在數據缺失問題，現有數據多為監測斷面數據，缺乏空間和時間連續性，無法滿足水生態承載力評估需要. 因此，該研究通過遙感反演的方法獲取了6期水質及水生生物數據，反演指標包括COD、TN、TP、氨氮濃度以及葉綠素含量.

1.3.2水質指標遙感反演

已有較多研究運用實測水質數據和遙感數據進行水質遙感反演，主要反演的水質指標包括COD、TP、TN、BOD、氨氮濃度以及葉綠素含量等. 該研究通過2019年9月和2020年10月北運河的實測數據和部分歷史數據，運用遙感反演模型獲得1995年、2000年、2005年、2010年、2015年和2018年6期的COD、TP、TN、氨氮濃度以及葉綠素含量數據，從而完善水生態承載力評估體系. 2019年設置水樣采樣點34個，2020年加密到47個采樣點，并補充采集了樣點水體光譜信息. 光譜信息運用美國ASD公司(Analytical Spectral Devices,Inc.)生產的便攜式地物光譜儀(ASD Field Spec 4)收集[28]. 監測樣品的采集、保存、測定均按照相關技術規范處理[29-31].

在MATLAB 2017和Envi 5.3軟件中結合實測光譜數據建立遙感影像和水質數據的反演模型. 在MATLAB 2017軟件中對實測光譜數據進行平滑降噪、散射校正等預處理，并結合實驗室實測水質數據進行特征波段的選取，然后對實測光譜和實測指標值運用最小二乘法進行建模，得到基于實測光譜數據的反演模型，并運用反演模型得到水質指標的預測值[32]. Landsat遙感影像數據經過輻射定標、Flash大氣校正、裁剪等預處理，將所選擇的反演模型特征波段范圍與Landsat遙感影像數據波段范圍進行對比，選擇遙感影像中對應波段與實測光譜預測值運用線性回歸方法進行建模，從而完成實測光譜數據對遙感影像數據的校正；將反演模型運用于各年份研究區遙感影像，從而得到北運河北京段范圍內各年份水質數據[33-34]. 得到各年份水質數據后，參考常用的精度檢驗方法，選用決定系數(R2)和均方根誤差(RMSE)并結合已有數據對反演精度進行評估[35]. 將反演出的水質指標值與從北京市北運河管理處收集到的歷史水質資料進行對照，并且與《北京市水資源公報》中該區域的水質標準等級指示的水質指標范圍進行對比[36].

2 結果與分析

2.1 水質指標遙感反演精度評價

實測指標數據與光譜數據的相關系數在0.85以上，遙感影像與實測光譜預測值建模得到的相關系數在0.8以上，COD、TP、TN、氨氮濃度以及葉綠素含量的相關系數(R2)和均方根誤差(RMSE)如表4所示. 由表4可見，水質指標反演整體相關性較好，反演模型精度較高，因此可用于評價指標體系的計算.

表4 水質指標反演精度評價

《北京市水資源公報》顯示，1986—2018年，北運河流域水質基本為Ⅴ類或劣Ⅴ類，將流域水質范圍與各樣點歷年水質指標等級以及北京市北運河管理處部分水質等級歷史數據進行對比，判定反演的水質指標等級基本與公布水質等級符合，因此反演精度與流域實際相符. 水利部和原環境保護部對北運河水質考核要求為Ⅴ類水質標準，由于承載力評估體系中水質指標反映的是與北京市市控斷面等水質控制標準相比較的達標情況，通過水質超標(Ⅴ類水質標準)倍數來表征，因此反演精度與《北京市水資源公報》中區域水質等級基本一致即可滿足評價精度要求.

2.2 北運河北京段水生態承載力時空變化分析

通過構建的北運河北京段水生態承載力評估體系，計算得到1995年、2000年、2005年、2010年、2015年、2018年6期區域整體和子流域的水生態承載力，結果如圖2所示. 區域水生態承載力評估結果表明，北運河北京段水生態承載力在近30年間呈周期性變化. 北運河北京段水生態承載狀態變化大致可分為3個階段：1995—2000年、2000—2015年、2015—2018年. 1995年區域綜合承載力分值為51.60，2000年區域綜合承載力處于超載狀態，分值為30.53，因此1995—2000年區域綜合承載狀態下降，2000年左右區域綜合承載狀態為近25年間最差. 2005年起區域承載狀態有所恢復，區域承載力綜合分值為52.43，為臨界超載狀態；2010年和2015年區域承載力綜合分值分別為42.19和51.70，處于臨界超載狀態，2015年區域綜合承載狀態較2005年有所下降. 2015—2018年區域綜合承載狀態有所恢復，至2018年，區域綜合承載分值達到63.91，已基本處于安全承載狀態. 2018年區域整體水生態承載力較1995年和2015年高，說明近年來區域水生態承載力有上升趨勢.

圖2 北運河北京段總體水生態承載力時序圖

從北運河北京段子流域承載狀態(見圖3)來看，北運河上游承載狀態好于中下游. 近30年上游的沙河子流域、藺溝河子流域以及中游的壩河子流域承載狀態較好，均處于臨界超載狀態；而承載狀態較差的子流域主要有涼水河、清河、溫榆河，近30年間基本處于超載狀態. 就具體階段而言，藺溝河自2010年以來承載狀態下降趨勢明顯. 至2018年，北運河承載狀態最好的子流域為壩河、沙河，水生態承載力分值分別為53.72、48.12；承載力最差的子流域分別為涼水河、北運河、清河，水生態承載力分值分別為28.59、30.33、30.97.

圖3 北運河北京段子流域水生態承載力時序圖

對各子流域綜合承載力指標值作單因素方差分析以檢驗子流域間綜合承載力的差異顯著性，結果如圖4所示. 進行單因素方差分析前首先進行了正態分布檢驗，Shapiro-Wilk檢驗和K-S檢驗結果顯示，P值分別為 0.403 9 和0.956；方差齊性檢驗P值為0.894，在α=0.05的檢驗水平下，P>0.05，因此數據滿足方差齊性及正態分布，可進行單因素方差分析. 單因素方差分析結果顯示，P值為2.48×10-11，小于 0.000 1，表明各子流域水生態承載力差異性顯著. 沙河、壩河的水生態承載力顯著高于其他子流域，藺溝河和小中河水生態承載力顯著高于北運河、涼水河、溫榆河和清河，清河的水生態承載力分值最低，長期處于超載狀態. 沙河、藺溝河子流域處于北運河水系上游，早期受人為擾動較小，因此整體承載狀態好于其他子流域；壩河近年來實施了一系列的綜合整治措施，生態恢復效果明顯.

注： 不同的小寫字母表示差異性顯著.

2.3 區域水生態承載力關鍵要素識別

運用多層感知器神經網絡進行建模、結合Pearson相關系數和Spearman相關系數方法得到區域水生態承載力的影響因子重要性，識別出區域水生態承載力關鍵要素. 結果(見圖5)顯示，相關系數與神經網絡建模識別出的關鍵因子基本一致，重要性最大的影響因子為河岸帶林草覆蓋率. 多層感知器神經網絡得出河岸帶林草覆蓋率標準化的重要性為100%，其次為河流連通性、人均水資源量、COD超標倍數、生境整體性指數，標準化的重要性分別為62.4%、52.4%、38.1%、34.3%. 從總指標層和分指標層來看，區域水生態承載力主要控制因子為水生態總指標中的水生生境指標，其次為水環境總指標中的子指標，水資源總指標的重要性較水生態和水環境指標低.

圖5 水生態承載力影響因子重要性

3 討論

3.1 社會經濟因素對北運河水生態承載力的影響

該研究運用綜合指標法從水資源、水生態、水環境3個層面構建了北運河北京段區域水生態承載力評估體系，結果顯示，北運河北京段區域水生態承載力在近30年間呈周期性變化，即1995—2000年呈下降趨勢，2000—2005年呈上升趨勢，2005—2010年呈下降趨勢，2018年較2015年呈上升趨勢，并且2000—2010年水生態承載力的變化趨勢與京津冀地區水生態承載力的變化趨勢基本一致[37]. 已有研究[38-39]表明，水生態承載力反映的是生態系統對區域社會經濟系統的承載能力，區域所能承載的人口、GDP必須與當地的資源、環境相協調，因而歷年人口、社會經濟指標的變化趨勢也可以反映區域水生態承載力變化. 從區域人口變化來看，根據《北京市統計年鑒(1995—2018)》，2018年末北京市常住人口為 2 154.2×104人，較2015年的 2 170.5×104人有所下降，常住人口下降對緩解環境壓力起到了一定作用，人口帶來的環境壓力變化與北運河整體水生態承載力變化趨勢一致. 從區域經濟發展速度來看，GDP環比增速與區域整體水生態承載力變化趨勢基本一致，均在1995—2000年有所下降、2000—2005年有所上升、2005—2015年持續下降，但是2015—2018年GDP環比增速下降，區域整體水生態承載力有所提升，說明經濟增速水平在一定程度上反映了區域生態承載能力，但是經濟增長不一定就會造成區域水生態承載力的下降，這與關鍵控制要素識別結果相符. 因此，從環境、人口、社會經濟等各角度綜合分析，北運河北京段區域水生態承載力變化趨勢與區域實際狀況基本相符. 區域整體和部分子流域水生態綜合狀況的改善可能與北京市近年來采取的一系列水環境綜合整治措施有關. 根據北京市北運河管理處相關資料，發現北京市從2005年開始進行黑臭水體修復、污水設施升級改造以及施行水環境區域補償機制，因此近幾年北運河部分區域水質得到了明顯改善[36].

3.2 生態環境因素對北運河北京段水生態承載力的影響

沙河子流域和藺溝河子流域均處于北運河上游，有研究[40-41]表明，這兩個子流域地勢海拔相對較高，土地利用類型以林地、草地為主，受人類干擾較小，生態環境基礎狀況較好，生態系統的自我維持與凈化能力較強，水生態承載力較高，其中沙河子流域歷年植被覆蓋度均為各子流域中最高，除2000年以外均在60%以上. 據《北京市統計年鑒》數據顯示，沙河子流域和藺溝河子流域歷年地均GDP和地均人口規模在所有子流域中屬中等偏下，城市開發強度相對較低，說明沙河子流域和藺溝河子流域受人為擾動較小，生態環境狀況相對較好[40,42]. 但2010年后藺溝河子流域水生態承載力呈明顯下降趨勢，北京市北運河管理處已有資料顯示，2018年藺溝河日排污量在 4×104m3以上，明顯高于其他子流域，說明部分子流域可能存在治理不徹底的狀況，從而導致子流域水生態的惡化. 此外，2010—2015年北運河水系中游壩河子流域水生態承載力呈上升趨勢，2018年左右水系上游沙河子流域和中游壩河子流域水生態承載力分值明顯較其他子流域高，這與已有研究結果[40,43]基本一致. 2010年前北運河水系中下游子流域承載狀態較差，最差的子流域主要有清河子流域和溫榆河子流域，較差的有涼水河子流域和小中河子流域，這些子流域均處于超載狀態. 已有研究[44-45]表明，溫榆河存在大量未達標污水排水口，清河和涼水河是北京市區主要排水河道，清河承載著城市化程度較高的北京市海淀區的人口、社會經濟增長壓力，而小中河沿岸工業污染排放較大，這些因素給區域水生態綜合承載能力帶來了較大壓力. 2010年后除涼水河外，其他子區域水生態承載力均有不同程度改善，這一結論基本與北京市北運河管理處持續監測結果一致.

北運河北京段水生態承載力最主要的影響因素為水生態，如河岸帶林草覆蓋率. 《北京市通州區統計年鑒》顯示，2005—2018年，通州區綠化覆蓋率從34.87%升至45.48%，提升趨勢顯著. 已有研究[46-48]表明，生態環境的改善，尤其是區域植被覆蓋率的提升對于區域生態輸水能力具有重要影響. 也有研究[49-51]表明，區域植被覆蓋度能調節氣候，有效控制水土流失，有助于水質凈化，并能進行雨洪滯蓄、涵養水源. 北運河水生態承載力受水環境影響也較為顯著，如COD和TN濃度. 根據《通州區水環境質量限期達標規劃(2018—2020年)》，通州區北運河流域污染來源主要為城鎮生活污水、畜禽養殖污染、農村生活污水等，主要污染物為COD、氮和磷，各類污染物中COD排放量較為突出，這與關鍵控制要素中污染物排放規模等要素相符合. 此外，北運河流域水生態承載力也在一定程度上受到水資源變化的影響，主要影響因素為人均水資源量，其反映了流域水資源稟賦. 1995—2018年北京市水資源總量與北京市人均水資源量變化趨勢基本一致，2000年北京市水資源總量為16.86×108m3，較1995年的47.62×108m3顯著下降，2005—2015年呈周期性變化，2018年為35.46×108m3，較2015年的26.76×108m3有所增加. 2000年北運河北京段區域整體水生態承載力較1995年下降明顯，2005—2015年呈周期性變化，2018年水生態承載力較2015年有所提升(見圖3)，變化趨勢基本與北京市水資源總量和人均水資源量一致. 北運河北京段區域水資源稟賦不佳，水資源補給以降雨為主，同時水資源總量年內變化明顯，水資源十分寶貴；所以區域水資源綜合利用效率的改善也可能是區域水生態承載力提升的重要途徑[43,52]. 總體而言，北運河北京段區域仍需加強對河流沿岸濕地及植被的保護，并通過人工恢復與自然恢復相結合的方法進一步改善流域生態環境[52-53].

4 結論

a) 近30年來，北運河北京段水生態承載力呈周期性變化. 2000年區域整體水生態承載力分值為30.53，為歷年最差，到2018年區域整體水生態承載力分值為63.91，呈上升趨勢，但仍需加強水環境污染防治及水資源節約利用.

b) 北運河北京段各子流域間歷年水生態承載力差異顯著. 1995—2010年，水生態承載力最好的子流域為沙河子流域；2015年、2018年水生態承載力最好的均是壩河子流域. 除2015年水生態承載力最差的子流域為涼水河子流域外，其余評價年份水生態承載力最差的子流域均為清河子流域.

c) 通過對北運河北京段區域歷年承載狀態影響因素進行分析可知，水生態要素對區域承載狀態影響最為顯著，其次為水環境要素.