山丘區小流域單位線峰值經驗公式推求研究

2019-06-04 06:00:34姬榮彬

中國農村水利水電 2019年5期

姬榮彬，葉磊，吳劍，郭良，史可，張弛

(1. 大連理工大學水利工程學院，遼寧大連 116024；2. 中國水利水電科學研究院，北京 100038)

0 引言

我國山丘區面積占國土面積的67%，人口占全國總人口的56%[1]，廣大山丘區山高坡陡、溪河密集，局部短歷時強降雨頻發，洪水陡漲陡落，常易爆發山洪災害，給人們的生命財產安全帶來巨大損失[2]。短歷時強降雨引發的山丘區小流域洪水，在陡峭地形影響下，匯流速度快，時間短，致災迅速，給山丘區小流域的暴雨洪水分析計算帶來了很大的挑戰。山丘區水文氣象站點稀少，水文資料稀缺，是典型的無資料或缺資料地區[3]；同時山丘區小流域暴雨所引發的水文響應更為直接，短歷時的強降雨常常是山洪致災的動力因子[4]，雨強及下墊面異質性的影響更為凸顯。因此，探索山丘區小流域的暴雨洪水計算方法對山洪防災減災具有重要意義。

無資料條件下，由地貌學的角度尋求水文響應的合理解釋是解決小流域暴雨洪水分析計算的重要途徑。地貌單位線[5,6]從地貌成因上揭示了流域水文響應，即“單位線等價于呈弱粒子相互作用下的水滴匯流時間的概率密度函數”，通常稱之為降雨徑流形成的“粒子說”。這一理論的建立，使得推求單位線的過程轉變為計算水滴匯流時間概率分布的過程，Maidment等[7]、芮孝芳等[8]基于該理論，進一步提出了由DEM推求分布式地貌單位線的方法；孔凡哲等[9]針對坡地匯流不應被忽略的事實，通過在坡地和河網分別采用不同的公式計算流速，提出了基于空間分布流速場的地貌單位線；郭良等[10]將雨強因子引入流速計算公式，考慮了雨強對流域水文響應的影響，并依托山洪災害調查評價平臺，提取了全國范圍內不同時段長、不同雨強下的分布式地貌單位線用于山洪預警指標分析，取得了良好的效果。分布式地貌單位線考慮了雨強、流域下墊面異質性，與小流域匯流非線性響應突出的特點相符，在山丘區小流域洪水匯流計算中具有較大應用價值，但其推求過程繁瑣，且需要一定精度的地形地貌、土地利用等資料，因此目前在無資料山丘區小流域應用較少。

建立流域水文過程與下墊面地形地貌特征的直接定量關系，對于山丘區小流域暴雨洪水分析計算具有現實意義。考慮山洪預警多以臨界雨量或臨界流量為預警指標[11]，更多關注于實際暴雨洪水過程中洪峰流量的大小。因此，本文基于考慮雨強、流域下墊面異質性的分布式地貌單位線，開展單位線峰值與下墊面因子的經驗關系研究，首先通過成因分析結合逐步回歸識別對小流域單位線峰值影響顯著的下墊面因子，進而構建兩者間的經驗公式。研究成果可為山洪防災減災工作提供實際參考。

1 研究流域及數據

1.1 研究流域

我國山洪災害防治基礎數據及開發應用中，首次劃分了全國范圍內53萬個山丘區小流域，基于數字高程模型和高分辨率遙感影像提取了小流域面積、坡降等基礎屬性，構建了小流域基礎屬性庫；并以降雨徑流形成的“粒子說”[5,6]為基礎，提取了各小流域不同時段長(10，30，60 min)、不同降雨強度(5，10，30 mm)條件下的分布式地貌單位線，形成了分布式地貌單位線庫[10]。依托山洪災害防治數據中小流域基礎屬性庫和分布式地貌單位線庫，考慮不同地區的產匯流機制不同導致流域水文響應差異顯著，本文選取地處不同省份的伊河(河南省)、龍河(重慶市)、涇河(甘肅省)、拒馬河(河北省)、牛欄江(云南省)共5個典型流域進行研究，5個流域整體上位于我國地形的二級階梯及一二級、二三級階梯的過渡帶，流域位置分布如圖1所示。

圖1 研究流域位置示意Fig.1 Sketch map of the location of study region

表1列出了本文所選研究流域的主要水文特征及小流域劃分情況。所選5個研究流域多位于不同的氣候類型區，水文特征差異明顯。伊河流域屬溫帶大陸性季風氣候，多年平均降雨量為791 mm，徑流量受降水影響明顯，每年6-8月為汛期；龍河流域屬亞熱帶季風濕潤氣候區，立體氣候明顯，春早冷暖多變，夏熱多雨多伏旱，秋涼多綿雨，冬季干冷，龍河支流眾多，為典型的山地冷水急流型河流；涇河流域地處副熱帶氣候區，多年平均徑流量21.4 億m3，徑流年內分配不均勻，夏季大于秋季，冬季最小；拒馬河地處溫帶大陸性季風氣候區，具有春季干旱多風沙、夏季炎熱多雨水、秋季晴朗涼爽、冬季寒冷少雨的特征，降水量在年內和年際的分布都很不均勻，汛期6-9月，豐水年汛期產生的徑流量一般占全年徑流量65%～75%；牛欄江流域內多數地區屬暖溫帶高原季風氣候，年溫差小、日溫差大，干濕季節分明，多年平均徑流量38.1 億m3，水力資源十分豐富，干流理論蘊藏量占全水系的90.7%。

表1 研究流域主要水文特征及小流域劃分情況Tab.1 Key hydrologic characteristics and small basin division in study region

此外，5個流域都屬于地形起伏較大的山丘區地帶，屬于山洪災害易發區，山洪形勢極為嚴峻。例如，伊河流域每年汛期(6-8月)至少出現2次及數次上游洪水致使山洪暴發；龍河則流淌于崇山峻嶺之間，集雨面積大，河床小，每逢下雨就極易形成山洪，據石柱縣防辦數據顯示，新中國成立以來，該縣縣城4次被龍河洪水淹沒，上游中益、橋頭兩鄉鎮被淹沒的次數則更多。涇河洪水同樣猛烈，是渭河及黃河洪水主要源地之一，涇河一場洪水即可對下游構成威脅。拒馬河流域和牛欄江流域同樣如此，山區洪水陡漲陡落，極易引發山洪災害。所選5個流域均屬于山洪災害易發和頻發的高危區域，且氣候水文特征不同，因此，認為所選5個流域對于本文的研究內容而言具有一定代表性。5個流域內共劃分了小流域803個，其中，牛欄江流域內7個坡度為0的小流域予以剔除，最終選取了796個小流域開展研究。

1.2 小流域分布式地貌單位線

單位線是流域匯流計算的常用方法，傳統單位線是基于線性假定的，而實際情況中，即使降雨空間分布均勻，對同一種徑流成分來說，流域單位線也并非固定不變，而是隨著降雨強度變化而變化[12]。受降雨強度和流域內下墊面時空變異性的影響，山丘區小流域匯流過程呈強非線性特征，給流域峰值計算帶來很大難度。分布式地貌單位線以“粒子說”[5,6]理論為依據，充分考慮了雨強對單位線的影響，同時還考慮了流域下墊面局部地形、植被覆蓋等的空間分布特征，因此是一種物理機制較為明確且更符合實際情況的單位線。

本文收集到796個小流域不同時段長(10、30、60 min)、不同降雨強度(5、10、30 mm/h)下的分布式地貌單位線。某小流域時段長為10 min對應不同雨強的單位線如圖2所示，可見，雨強越大，單位線越尖瘦，峰值越大，峰現時間越早，非線性現象明顯。不同時段的單位線可通過S曲線相互轉換，為研究方便，本文以時段長為10 min、雨強為5 mm/h的分布式地貌單位線為例，提取單位線峰值，分析其與小流域下墊面因子間的經驗關系。各典型流域的小流域單位線峰值分布箱線圖如圖3所示(箱的寬度與小流域個數呈正比)，各單位線峰值來自不同小流域樣本，異常值未剔除。

圖2 某小流域時段長為10 min對應不同雨強條件下的分布式地貌單位線Fig.2 Distributed geomorphological unit hydrograph under 10 min and diverse rainfall intensity in a small basin

圖3 各典型流域時段長為10 min、雨強為5 mm/h的分布式地貌單位線峰值分布箱線圖Fig.3 Box plot of peak values of distributed geomorpho-logical unit hydrograph under 10 min and 5 mm/h in each typical watershed

1.3 小流域下墊面因子

流域下墊面地物地貌特征深刻影響著流域水文響應，而流域下墊面因子是衡量下墊面地物地貌特征的重要參數。本文提取了流域下墊面諸多特征因子，各因子的名稱、符號表示、釋義及單位信息詳見表2。所收集到的5個典型流域的各下墊面因子數據見表3(限于篇幅，未展示所有數據記錄)，部分數據來源于全國山洪災害防治小流域基礎屬性庫。

表2 小流域下墊面因子及釋義Tab.2 Introduction and definition of underlying surface factors of small basins

表3 研究區域內各小流域下墊面因子數據Tab.3 Underlying surface factors values of each small basin in study region

2 小流域下墊面因子選取

分布式地貌單位線是流域水文響應的直觀表現，綜合反映了流域內各下墊面因子對水流匯集過程的影響。流域下墊面因子識別和選取是山丘區小流域洪峰流量計算的關鍵。為了充分識別影響小流域單位線峰值的關鍵因子，以下首先從成因上進行定性分析，接著利用逐步回歸法進行定量分析，兩種方法相結合以保證因子選取的正確性和可靠性。

2.1 基于成因的定性分析

水文學家普遍認為，流域下墊面對于流域洪峰形成起著至關重要的作用，當降雨、蒸散發等氣象過程相同時，流域出口斷面流量過程的形成將取決于其下墊面條件[13]。但是，對于哪些下墊面因子能夠顯著影響洪峰過程，并無統一的規律可循，多數學者曾進行過相關研究。Snyder[14]根據統計規律認為，洪峰流量與流域面積、主河道長度、流域出口到形心的距離有關，因此構建了相應經驗公式；陳家琦等[15]認為，洪峰流量與流域面積、主河道長度、主河道坡降有關，因此提出了推理公式；其他學者[16,17]同樣也提出了一些計算公式，所選取的下墊面因子多圍繞流域面積、主河道長、坡降、彎曲率等一些常用屬性展開。若從定性分析來看，當其他條件均保持一致時，流域面積作為表征流域大小的常用因子，面積越大，往往意味著水量大，進而所形成的洪峰大；坡降反映流域地形起伏，坡降越大，水流流速越快，峰值必然大；主河道是輸移洪水的主要通道，主河道越長且越彎曲，水流到達出口前所經過的路徑就越長，峰值就會衰減；因此，利用這些常用的下墊面因子構建洪峰經驗公式是比較合理的。

2.2 基于逐步回歸的定量分析

在定性分析的基礎上，為了更為準確的分析各下墊面因子對單位線峰值的影響，接著采用逐步回歸法進行定量分析。逐步回歸法是一種逐步篩選自變量的多元線性回歸方法，由于其實質是建立最優回歸模型，且能夠在一定程度上克服因子間的多重共線性[18]，常被用于分析水文因子的響應關系[16-19]。

考慮流域匯流的非線性響應及不同流域匯流特性差異，本文以分布式地貌單位線峰值的自然對數值作為回歸因變量，以小流域各下墊面因子的自然對數值作為回歸自變量，分別在伊河、龍河、涇河、拒馬河、牛欄江5個典型流域構建逐步回歸模型，取定常用進入和移除水平分別為0.05和0.1，經逐步引入與剔除，得到5個流域的分布式地貌單位線峰值逐步回歸結果，見表4(以因子選取為目的，未列出回歸常量)。有關逐步回歸的實施細節可參見文獻[18]。

表4 不同典型流域的分布式地貌單位線峰值逐步回歸結果Tab.4 The stepwise results of peak values of distributed geomorphological unit hydrograph in each typical watershed

注：**表示通過水平為0.05的顯著性檢驗。

表4中所有回歸方程及回歸系數都通過了顯著性檢驗，決定系數R2接近于1，表明由下墊面因子所構成的回歸方程對小流域單位線峰值具有較強的解釋作用。其中，伊河、龍河、牛欄江三個流域優選出的下墊面因子一致，均為F、J、Lm以及La；涇河流域的優選結果中含S，不含J；拒馬河流域中J、S均包含在內，也即不同流域優選出的下墊面因子存在差別，實際中每個流域的地形地貌、匯流特性均有差異，因此認為優選結果的不一致與實際是相符的。總體上，F、Lm、La三種下墊面因子在五個流域內均被優選為顯著性因子，表明這三種因子對單位線峰值影響最為顯著，Lm、La均為表征流域內水流匯流路徑的因子，但前者較后者更為顯著(前者的回歸系數絕對值明顯大于后者)；J、S同為表征流域下墊面地形起伏的因子，但前者顯著于后者(前者被優選為顯著性因子4次，后者2次)；其他下墊面因子Lmax、Lv、R未被優選出，表明對單位線峰值的影響不顯著。結合定性分析結果，最終確定優選結果為流域面積F、最長匯流路徑長度Lm及最長匯流路徑坡降J。這與我國設計洪水規范中計算小流域設計洪水的推理公式法[15]所采用的流域下墊面因子(流域面積、主河道長度、主河道坡降)基本一致。

3 分布式地貌單位線經驗公式推求

3.1 經驗公式擬合

對于構建何種類型的經驗公式，多數學者[13-17,20]經過深入研究，認為峰值經驗公式應為冪函數型比較合理，故根據前述小流域下墊面因子選取結果擬定了如下經驗公式，

(1)

式中：Qm為分布式地貌單位線峰值，m3/s；J、F、Lm意義同前；k，a，b，c為待定參數。

為確定式(1)在5個典型流域的峰值計算效果，以單位線洪峰實際值與擬合值的均方根誤差RMSE最小作為目標函數，采用遺傳算法優化公式中參數，分別得到式(1)在5個典型流域的擬合情況(表5)。為分析需要，同時以相對誤差不超過20%作為小流域洪峰合格標準，統計合格率QR。

(2)

(3)

式中：RMSE為單位線實際峰值Qm,obs與擬合峰值Qm,sim的均方根誤差，m3/s；n為小流域個數。

3.2 結果分析討論

由表5可知，從擬合RMSE來看，伊河流域的擬合效果最好(RMSE=0.114)，其次為涇河流域(RMSE=0.131)；從擬合合格率QR來看，涇河流域的擬合效果最好(QR=94.2%)，其次是龍河流域(QR=93.4%)。總體上，式(1)在五個流域的擬合效果都非常好(RMSE均比較小，QR皆大于90%)。參數a，b，c在不同流域的不同取值，意味著不同流域內下墊面因子對單位線峰值的非線性響應不同，但總體變幅不大；參數k作為量綱轉換系數，不同流域取值不同，但必為正值；參數a、b皆大于0，表明流域面積F及最長匯流路徑坡降J與單位線峰值呈正相關，參數c小于0，表明最長匯流路徑長度Lm與單位線峰值呈負相關。也即流域面積F越大，單位線峰值相應越大，這符合前文定性分析結果，與文獻[14-17,21]的研究結果相一致；最長匯流路徑坡降J表征流域的地形起伏，坡降大導致水流流速快，所形成峰值必然大，文獻[21]認為最長匯流路徑坡降與洪峰流量呈負相關，這與本文所得出的結果不一致，原因是文獻[21]利用兩者的相關系數判定其相關關系，實際上洪峰流量受多種因素影響，僅用兩者的相關系數判定具有較大的不確定性，本文采用的逐步回歸法剔除了流域面積等其他因素影響，所得結果與實際情況較為一致；最長匯流路徑長度Lm表征實際流域內的主河道長度，主河道是水流匯至流域出口斷面所必經的通道，其他條件一定時，最長匯流路徑越長，單位線峰值相應越小(洪水波坦化作用)。

表5 式(1)在五個典型流域的擬合結果Tab.5 Fitting results by Eq(1) in each typical watershed

為直觀地比較經驗公式的擬合效果，繪制5個流域的擬合效果圖(圖4)，絕大多數擬合點均分布于合格范圍以內，緊緊圍繞在45°線周圍，表明由顯著性因子流域面積F、最長匯流路徑長度Lm及坡降J所構建的經驗公式擬合效果十分理想。其中，伊河、涇河流域擬合點距與45°線更為貼近，與從RMSE分析的結果相符；此外，5個流域中位于合格范圍以外的點大多集中于小洪峰段，即峰值較小的小流域更難滿足合格條件，分析原因是由于采用相對誤差指標，峰值較小時，20%的合格域很小，多數洪峰較小的小流域難于滿足合格條件。

為了更好地應用本文所推求的小流域分布式地貌單位線峰值經驗公式，本文通過分析合格小流域的下墊面因子特征進一步確定式(1)的適用條件和邊界條件。提取五個典型流域中所有計算合格的小流域的流域面積F、最長匯流路徑長度Lm及坡降J數據，結果如圖5所示。經相關性分析發現F與Lm有較好的相關性(相關系數為0.79)，而與J相關性較差(相關系數為-0.04)，因此只需分析F和J的相互關系以確定式(1)的適用范圍。

圖4 式(1)在五個典型流域的擬合散點圖Fig.4 Scatter plot of fitting results by Eq(1) in each typical watershed

圖5 合格小流域的F、Lm、J數據Fig.5 F、Lm and J data of qualified small basins注：圖中F、Lm、J的單位依次為km2、km、%。

以所有合格小流域的下墊面因子F和J繪制F-J邊際圖(邊際圖是在x軸和y軸邊際中包含箱線圖的散點圖)，如圖6所示。圖6的主體是合格小流域的F和J的散點圖，散點圖右側和上側的箱線圖分別表征F和J的數據分布，從箱線圖中可直觀看出數據的異常值范圍，以兩個箱線圖中正常值與異常值間的臨界點(F=31.7 km2，J=9.8%)分別作水平線和豎直線，將整個散點圖區域劃為4個區，分別是Ⅰ區(F≤31.7 km2,J≤9.8%)、Ⅱ區(F≤31.7 km2,J>9.8%)、Ⅲ區(F>31.7 km2,J≤9.8%)和Ⅳ區(F>31.7 km2,J>9.8%)。從圖6各個區中的合格小流域點距的分布情況可明顯看出，Ⅰ區是式(1)的主要適用區域，大量點距集聚于此；Ⅱ區和Ⅲ區內點距稀疏，屬于適用式(1)的邊界地帶，若要采用式(1)計算峰值，建議首先做適用性分析，實地驗證式(1)的試用效果；Ⅳ區屬于式(1)的非適用區域，不建議采用本文所提出的方法。