張金萍，王宇昊

（1.鄭州大學(xué)水利科學(xué)與工程學(xué)院，鄭州450001；2.鄭州大學(xué)黃河生態(tài)保護與區(qū)域協(xié)調(diào)發(fā)展研究院，鄭州450001）

0 引言

城市河道徑流是城市水文過程的重要組成部分，也是城市用水和水資源管理的重要內(nèi)容。與自然河道徑流相比，城市河道徑流水源組成更為復(fù)雜，采用適宜的方法對城市河道徑流進行合理分割，不僅有助于對城市河道徑流水源組成的科學(xué)認識，而且還將提升城市水文模擬精度與城市用水安全［1］。近年來，許多學(xué)者針對徑流分割做出了相應(yīng)的研究［2-4］，例如張泳華等采用數(shù)字濾波法和平滑最小值法對東江流域控制性水文站博羅站的實測日徑流量進行徑流分割，探究了降水和人類活動對基流的影響［5］；馬曉婧等利用改進退水常數(shù)的數(shù)字濾波法對拒馬河徑流過程進行分割，分析了拒馬河豐、平、枯年的基流變化特征［6］。這些研究中徑流分割研究的對象主要是河川基流，而城市河道中存在大量的生活污水、工業(yè)廢水，徑流的成分發(fā)生了顯著的改變，特別是對于上游天然來水被全部攔截，無上游來水的城市河道，更需要對其進行徑流分割，從而進一步探究徑流的組成部分及其演變過程。

下墊面是影響降雨-徑流關(guān)系的重要因素［7，8］，伴隨著城市化的發(fā)展，城市下墊面對徑流的影響得到了廣泛研究，例如李慧等運用水文特征參數(shù)時間序列法探究了西安市下墊面變化對城市河流灞河徑流的影響［9］；趙彥軍等利用SWMM 模型研究了下墊面城市化對濟南小清河流域產(chǎn)匯流的影響［10］。然而城市河道尤其是無上游來水的城市河道，生活污水、工業(yè)廢水的排放主要受人口、社會經(jīng)濟等條件的影響，導(dǎo)致總徑流不能準確地反映城市真實的集水特性，因此通過徑流分割得出的直接徑流能更準確地反映下墊面對降雨-徑流關(guān)系的影響。

本文以穿鄭州市而過的賈魯河為例，引入數(shù)字濾波法、滑動最小值法和時間步長法，對鄭州市中牟水文站2003-2016年的逐日徑流進行徑流分割，篩選出最優(yōu)分割方法后，采用年直接徑流系數(shù)表征鄭州市區(qū)的集水特性，結(jié)合鄭州市區(qū)5 期土地利用/覆蓋數(shù)據(jù)，構(gòu)建年直接徑流系數(shù)和各類下墊面面積的多元回歸模型，并對2017年和2018年的年直接徑流系數(shù)進行預(yù)測，從而更準確地反映下墊面城市化對降雨-徑流關(guān)系的影響。

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

1.1 研究區(qū)域概況

賈魯河屬淮河水系，是淮河支流沙潁河的支流，發(fā)源于新密市（隸屬鄭州市），向東北流經(jīng)鄭州市區(qū)，至市區(qū)北郊折向東流，流經(jīng)中牟、開封后，在周口市入沙潁河，最后流入淮河。賈魯河全長246 km，流域面積5 896 km2，主要支流有金水河、索須河、熊兒河、七里河、東風渠等［11］。

本文研究區(qū)域為賈魯河流域鄭州市區(qū)段，鄭州市區(qū)境內(nèi)賈魯河長68 km，流域面積993 km2，約占鄭州市區(qū)總面積的94%。由于賈魯河上游存在常莊水庫，尖崗水庫，河王水庫等水利樞紐，因此除賈魯河主河道外，其他支流均屬小河渠，基本無天然水源，已經(jīng)成為城市排污、農(nóng)灌退水及泄洪排水的渠道。此外，鄭州市區(qū)雨水的排放，主要是就近排入河道，市區(qū)內(nèi)的雨水分別從賈魯河、東風渠、金水河、熊兒河、七里河等泄洪河道匯入賈魯河后，排出市外［12］，河流水系圖見圖1。日徑流數(shù)據(jù)來源于中牟水文站實測資料，市區(qū)年降雨數(shù)據(jù)來源于鄭州市水資源公報。

圖1 鄭州市區(qū)賈魯河水系Fig.1 Drainage map of the Jialu River in urban area of Zhengzhou

1.2 數(shù)字濾波法

數(shù)字濾波法的基本原理是將流量過程轉(zhuǎn)化為數(shù)字信號，通過數(shù)字濾波器將信號分解為高頻和低頻，對應(yīng)地將徑流過程劃分為地表徑流和基流［13］。目前應(yīng)用廣泛的數(shù)字濾波法有：Lyne-Hollick 濾波法（F1）、Chapman 濾波法（F2）、Chapman-Maxwell 濾波（F3）和Boughton-Chapman濾波法（F4）。

F1法的徑流分割方程為：

式中：Qi為第i時刻的徑流量；Qd（i）為第i時刻的地表徑流量；Qb（i）為第i時刻的基流量；f1為濾波系數(shù)，一般取0.8～1.0；f2一般取0.5［14］。

由F1法改進的F2法基流分割方程［15］為：

F3法的基流分割方程為：

式中：k為退水系數(shù)，一般取0.95［16］。

F4法的基流分割方程［16］為：

式中：參數(shù)C一般取0.15［17］。

1.3 滑動最小值法

滑動最小值法是一種徑流快速分割方法，基本原理是通過確定時間步長N將逐年的日徑流過程分成若干個時間段，若某一時間段的最小流量值與拐點檢驗因子f的乘積小于等于相鄰時間段的最小流量值，則該時段的最小流量值即為流量過程線上的一個拐點。重復(fù)以上步驟，確定所有拐點后進行連接得到基流過程線。滑動最小值法目前主要有標準BFI（f）法和改進BFI（k）法兩種，相應(yīng)的拐點檢驗因子f和k分別取0.9 和0.979 15［18］。

1.4 時間步長法

時間步長法又稱為HYSEP 法，包括固定時間間隔法（FI）、滑動時間間隔法（SI）和局部最小值法（LM）等3 種流量分割方法［19］。運用該方法進行徑流分割，首先要計算的是直接徑流所持續(xù)的時間，目前計算該時間采用的是經(jīng)驗公式：

式中：T為直接徑流持續(xù)天數(shù)，d；A為流域面積，km2。

確定直接徑流持續(xù)天數(shù)后，取與2T最為接近且介于3～11 d的奇數(shù)作為分割的時間間隔參數(shù)N進行流量分割［20］。

1.5 多元線性回歸模型

多元線性回歸模型的核心在于建立因變量和兩個或兩個以上自變量之間的線性關(guān)系，模型建立前，需要分析各自變量對因變量影響的顯著性并選擇對因變量影響較大的自變量，從而合理構(gòu)建模型。假設(shè)因變量為y，自變量為x1、x2、…、xn，則多元線性回歸模型如下：

式中：β0為模型的常數(shù)項；β1～βn為偏回歸系數(shù)；ε為去除自變量對因變量影響后的隨機誤差。常數(shù)項和偏回歸系數(shù)的值可根據(jù)最小二乘法確定［21］。

2 流量過程線的自動分割

2.1 徑流分割的可行性論證

進入21世紀以來，鄭州市城市化程度不斷提高，伴隨著人口城市化和下墊面城市化程度的提高，城市河道的徑流量和徑流成分發(fā)生了較大的變化。根據(jù)中牟水文站2003-2016年逐日徑流資料，計算168 個月的月徑流量和逐月最小日徑流量見圖2。

圖2 中牟水文站逐月降水量-月徑流量-月最小日徑流量Fig.2 Monthly precipitation-monthly runoff-monthly minimum daily runoff of Zhongmu hydrological station

由圖2 可得非汛期（10-4月）的月徑流量和月最小日徑流量均呈現(xiàn)一定的增長趨勢，徑流的年內(nèi)變化進一步減小。結(jié)合實際情況，隨著地下水位下降，賈魯河及其支流的地下水已經(jīng)少之又少，加上上游尖崗、常莊等水庫的截留，導(dǎo)致賈魯河基流嚴重匱乏［22］；另一方面，伴隨著人口的增長和社會經(jīng)濟的發(fā)展，市區(qū)排放的生活污水和工業(yè)廢水顯著增長，而這一部分廢污水排放在賈魯河徑流中占據(jù)了較大的比例。

在擬采用的3 種流量過程線自動分割方法中，數(shù)字濾波法實質(zhì)上是將快速響應(yīng)的直接徑流視為高頻信號，將慢速響應(yīng)的地下徑流視為低頻信號，從而進行流量分割，而市區(qū)排放的生活污水和工業(yè)廢水等徑流成分在年內(nèi)較為穩(wěn)定，可以視為低頻信號；滑動最小值法實質(zhì)上是由總徑流曲線上的一系列轉(zhuǎn)折點所形成的折線，從而實現(xiàn)對流量的自動分割；時間步長法則是根據(jù)直接徑流的持續(xù)天數(shù)進行流量分割。因此，這3 種自動流量分割方法在無天然來水的城市河道中均具有一定的適用性，可以將賈魯河的徑流劃分為直接徑流和基流，其中直接徑流指的是當?shù)亟涤晷纬傻牡乇韽搅髁浚鲃t包括生活污水和工業(yè)廢水等。

許多研究表明：對于F1，當濾波次數(shù)N取3 即正向-反向-正向濾波，且濾波參數(shù)f1控制在0.90～0.95 間能夠得到較好的分割結(jié)果，同時增大濾波參數(shù)減少次數(shù)可以得到相同的結(jié)果［23-24］，在此基礎(chǔ)上，f1分別取0.9，0.925，0.95進行流量分割；對于F2，F(xiàn)3和F4 法均選擇默認的參數(shù)；對于滑動最小值法，由于研究區(qū)域面積較小，時間步長N分別取3，4，5；對于時間步長法，根據(jù)流域面積計算得時間間隔參數(shù)N=9。

2.2 徑流分割結(jié)果分析

通過上述3 類9 種方法對實測日徑流量進行徑流分割，分割后的結(jié)果用年基流指數(shù)BFI表示（圖3），其中BFI為當年基流量和總徑流量的比值。

圖3 徑流分割結(jié)果Fig.3 The results of runoff division

由圖3 可得：數(shù)字濾波法中，F(xiàn)1 法的參數(shù)f1越大，BFI值越?。籉2 和F3 法對應(yīng)的BFI基本一致且顯著偏小，在年際間波動很小；F4 的分割結(jié)果與F1 相近，但又存在一定差異，且相應(yīng)的BFI值在年際間波動較小；滑動最小值法中，無論是標準BFI（f）法還是改進BFI（k）法，N值越大，BFI越??；時間步長法中，F(xiàn)I和SI的分割結(jié)果基本一致，而LM法對應(yīng)的BFI值小于二者。

為了更直觀地了解日基流量分割情況，選擇2011年為典型豐水年，2008年為典型平水年，2012年為典型枯水年，對汛期（5月至9月）日徑流的分割結(jié)果分別見圖4-6。

圖4 2011年汛期基流過程線Fig.4 Base flow process line in flood season of 2011

圖5 2008年汛期基流過程線Fig.5 Base flow process line in flood season of 2008

由圖4、5 和圖6 可得，不同分割方法所得的基流過程線差異較大。無論是豐水年、平水年還是枯水年，數(shù)字濾波法中，F(xiàn)1法的基流過程線波動較小且隨洪峰變化不大；F2 法和F3 法的基流過程線基本重合，與F4的波動趨勢基本一致，但F4的基流過程線顯著高于F2 和F3，且隨著洪峰變化有明顯反映；滑動最小值法中，當N相同時，標準BFI（f）法和改進BFI（k）法的基流過程線基本一致，且兩種方法的基流過程線均有明顯的起伏波動，會出現(xiàn)基流量等于徑流量的情況。時間步長法中，F(xiàn)I 和SI法的基流過程線呈現(xiàn)明顯的階梯型波動，有較多的明顯拐點，不符合無上游來水的城市河道中的基流產(chǎn)生機理；LM 法的基流過程線較FI和SI略為平緩，基流變化趨勢隨著徑流過程變化趨勢變化，且滯后于洪峰過程。

圖6 2012年汛期基流過程線Fig.6 Base flow process line in flood season of 2012

2.3 最優(yōu)分割方法的確定

枯水指數(shù)（Q90/Q50）是反映地下水補給河川徑流特性的重要指標，其中Q90和Q50分別代表時段內(nèi)出現(xiàn)頻率大于等于90%和50%時的徑流量，即枯水指數(shù)在一定程度上可以代表河道中常年存在的徑流量占總徑流量的比例［25］，因此本文徑流分割的合理性可以用枯水指數(shù)來衡量。

將枯水指數(shù)與年總徑流量的乘積視為河道中常年存在的徑流量，用于和分割出的年基流量進行對比，將決定系數(shù)R2、納什效率系數(shù)NSE和平均相對誤差Re作為評價標準，驗證結(jié)果見表1。評價標準的公式如下：

式中：Qk為第k年觀測的年基流量；Sk為第k年估算的年基流量；為觀測的年基流量均值；為估算的年基流量均值；n為樣本個數(shù)。

由表1 得：F2 和F3 法估算的年基流量值偏小，NSE為負數(shù)且Re高達34%，說明這2 種方法在本研究區(qū)域的年基流量估算中不理想；FI 和SI 法的NSE介于0.7 和0.8 之間，但Re均超過10%，說明這兩種方法在本研究區(qū)域的年基流估算中效果一般，具有一定的模擬效果；同理，無論是BFI（f）還是BFI（k），N=5 時的基流估算效果最好，N=3 或4 時效果一般。F1、F4、BFI（N=5）、LM 法的NSE均超過了0.8，且Re控制在10%以下，說明這些方法對基流的模擬效果很好。其中，F(xiàn)1 法（f1=0.925）的NSE最高且Re最小，估算的基流量與實際觀測值的模擬效果最好，結(jié)合2.2中對典型年內(nèi)的基流過程線分析，選擇F1法（f1=0.925/N=3）的分割結(jié)果作為實際的徑流組成。

表1 各徑流分割方法的驗證結(jié)果Tab.1 Verification of runoff division methods

3 下墊面對直接徑流的影響

3.1 徑流分割成果

根據(jù)F1 法（f=0.925/N=3）分割的成果，得到中牟水文站2003-2016年徑流分割結(jié)果見圖7。

圖7 2003-2016年徑流量-年基流量-年直接徑流量Fig.7 Annual runoff-annual base flow-annual direct flow from 2003 to 2016

由圖7可得鄭州市賈魯河年徑流量和年基流量的變化趨勢基本一致，且年基流量顯著大于年直接徑流量，因此賈魯河年徑流量增大的主要原因是伴隨著鄭州市人口的增長和社會經(jīng)濟的發(fā)展，市區(qū)的生活污水、工業(yè)廢水等排放量增大，這一部分徑流對總徑流增長的貢獻率為83.96%，而年直接徑流對總徑流增長的貢獻率為16.04%，年直接徑流量和市區(qū)年降雨量的變化見圖8。

圖8 2003-2016年直接徑流量和年降雨量變化Fig.8 The change of annual direct runoff and annual rainfall from 2003 to 2016

由圖8 可得，2003-2016年降雨量整體上呈減少趨勢，而年直接徑流量呈增長趨勢，且同一年降水量級下的年徑流量顯著增大，其中2013年降雨量最小，但年直接徑流量大于2012 和2014年直接徑流量，這是由于2013年賈魯河單位生態(tài)補水量為10.16 m3/s，顯著大于其他年份，因此在徑流分割中一部分生態(tài)補水量被計入了直接徑流中，在對該年徑流分割結(jié)果進行相應(yīng)的調(diào)整后計算多年直接徑流系數(shù)見圖9，其中年直接徑流系數(shù)為年直接徑流深與年降雨量的比值，年直接徑流深由年直接徑流量和流域面積計算所得。

圖9 2003-2016年直接徑流系數(shù)Fig.9 Annual direct runoff coefficient from 2003 to 2016

3.2 城市下墊面變化分析

根據(jù)2000，2005，2010，2015，2018年五期精度為30 m 的鄭州市區(qū)土地利用數(shù)據(jù)（圖10），將其重分類為六類：耕地、林地、草地、城鄉(xiāng)工礦居民用地、水域和未利用地，并分別計算各時期鄭州市區(qū)的耕地、林草地，城鄉(xiāng)工礦居民用地和水域的面積見表2。

圖10 鄭州市區(qū)2000，2005，2010，2015，2018年土地利用/覆蓋分類Fig.10 Land use/cover classification of Zhengzhou City in 2000，2005，2010，2015 and 2018

結(jié)合圖10和表2可得2000-2018年鄭州市區(qū)的耕地面積持續(xù)減少，城鄉(xiāng)工礦居民用地面積持續(xù)增長，二者的變化量基本相等，即減少的耕地大多轉(zhuǎn)化成了不透水地面，且二者在2005-2010 這一時段內(nèi)變化最為顯著，在2010-2015 這一時段內(nèi)變化不大；林草地面積在2005-2010年顯著減少，在其他時段變化不大；水域面積在2000-2015年持續(xù)增加，但變化不大，在2015-2018年呈顯著減少趨勢。

表2 不同時期各類土地利用/覆蓋面積 km2Tab.2 Area of different land use/cover in different periods

假設(shè)2000-2005，2005-2010，2010-2015，2015-2018 等時段內(nèi)各類土地利用面積呈線性均勻變化，生成2003-2016年各類土地利用面積的時間序列，結(jié)合年直接徑流系數(shù)序列作相關(guān)分析見表3。

由表3 可得：2003-2016年，城鄉(xiāng)工礦居民用地面積和年直接徑流系數(shù)的相關(guān)性最好且呈正相關(guān)，這是由于隨著市區(qū)不透水面積的增大，下滲水量減小，從而導(dǎo)致降雨產(chǎn)生的直接徑流量增大；耕地面積、林草地面積和年直接徑流系數(shù)相關(guān)性較好且呈負相關(guān)；水域面積和年直接系數(shù)呈正相關(guān)且相關(guān)性較差，這與水域面積和徑流量呈負相關(guān)的實際情況不符，主要原因是水域面積變化較小，它對年直接徑流系數(shù)的真實影響被其他下墊面的變化所掩蓋。

表3 年直接徑流系數(shù)和各類下墊面面積的相關(guān)系數(shù)Tab.3 The correlation coefficient between annual direct runoff coefficient and various underlying surface areas

3.3 多元線性回歸模型構(gòu)建

由于研究區(qū)域面積一定且未利用地面積可忽略不計，因此城鄉(xiāng)工礦居民用地、林草地、水域、耕地等4 類土地利用的面積之和大致為一定值，即4個變量間存在確定的線性關(guān)系，因此為避免自變量出現(xiàn)顯著的共線性關(guān)系，在進行多元回歸分析時，至多選取3個變量作為自變量。由于減少的耕地面積基本上全部轉(zhuǎn)變?yōu)槌青l(xiāng)工礦居民用地，因此選取城鄉(xiāng)工礦居民用地面積（S1）、林草地面積（S2）、水域面積（S3）作為自變量，年直接徑流系數(shù)為因變量y，根據(jù)最小二乘法構(gòu)建多元回歸模型為：

利用線性關(guān)系的顯著檢驗方法（F檢驗法）對該模型進行檢驗，計算得F=27.261，大于在0.01 水平上的臨界值6.552，即該模型在0.01水平上顯著。模型模擬結(jié)果見圖11。

圖11 年直接徑流系數(shù)實測值和模擬值Fig.11 Measured and simulated annual direct runoff coefficient

由圖11 可得：模型模擬值與實測值相比，納什效率系數(shù)NSE、決定系數(shù)R2、平均相對誤差Re分別為0.891，0.891 和6.72%，且各年年直接徑流系數(shù)的相對誤差均小于20%。將2017年和2018年的城鄉(xiāng)工礦居民用地面積、林草地面積、水域面積輸入該多元線性回歸模型，預(yù)測這兩年的年直接徑流系數(shù)分別為0.349 和0.394。

利用數(shù)字濾波F1 法（f1=0.925/N=3）對中牟水文站2017年和2018年的逐日徑流量進行分割，可得年直接徑流量分別為1.695億m3和2.114億m3，結(jié)合年降水量計算年直接徑流系數(shù)分別為0.333 和0.359。該多元線性回歸模型對2017、2018年直接徑流系數(shù)預(yù)測的相對誤差分別為4.80%和9.75%，因此該多元線性回歸模型具有良好的模擬效果，在一定程度上量化了下墊面對年直接徑流的影響，更加清晰地反映了鄭州市區(qū)下墊面城市化對降雨-徑流關(guān)系的影響。

4 結(jié)論

將不同的徑流分割方法應(yīng)用于賈魯河這一無上游來水的城市河道，分析典型年內(nèi)的基流過程線，并將枯水指數(shù)作為衡量分割合理性的標準，從而確定最優(yōu)分割方法，結(jié)果表明：數(shù)字濾波F1 法（f1=0.925/N=3）在賈魯河的徑流分割效果最好，有效地將徑流分割為直接徑流和基流；2003-2016年，基流和直接徑流對年總徑流增長的貢獻率分別為83.96%和16.04%。

基于徑流分割成果，利用年直接徑流系數(shù)表征鄭州市區(qū)的集水特性，結(jié)合5 期土地利用數(shù)據(jù)構(gòu)建2003-2016年直接徑流系數(shù)和各類下墊面面積的多元線性回歸模型，并對2017年和2018年的年直接徑流系數(shù)進行預(yù)測，得出以下結(jié)論：年直接徑流系數(shù)總體上呈增長趨勢，且與城鄉(xiāng)工礦居民用地的面積顯著正相關(guān)；該多元回歸模型在0.01 水平上顯著，且具有良好的模擬效果，NSE、R2和Re分別為0.891，0.891 和6.72%。此外，該模型對2017、2018年直接徑流系數(shù)預(yù)測的相對誤差分別為4.80%和9.75%。

對于無上游來水的城市河道，利用流量過程線的自動分割方法從總徑流中分割出直接徑流，剔除與產(chǎn)匯流無關(guān)的生活污水、工業(yè)廢水等徑流成分，能更清晰地反映下墊面對降雨徑流關(guān)系的影響，對流域水資源的規(guī)劃與管理具有重要的理論與現(xiàn)實意義。□