月尺度天然徑流量還原計(jì)算模型研究

劉開(kāi)磊 王敬磊 祝得領(lǐng)

（淮河水利委員會(huì)水文局（信息中心） 蚌埠 233001 安徽淮河水資源科技有限公司 蚌埠 233000山東水之源水利規(guī)劃設(shè)計(jì)有限公司 濟(jì)南 250013）

水資源是與經(jīng)濟(jì)、社會(huì)發(fā)展密切相關(guān)的重要戰(zhàn)略資源，水資源的開(kāi)發(fā)利用程度反映著地區(qū)生產(chǎn)力發(fā)展水平。由于區(qū)域水資源量十分有限，水資源可利用量在相當(dāng)大程度上制約著區(qū)域社會(huì)經(jīng)濟(jì)的發(fā)展。因此科學(xué)開(kāi)展區(qū)域水資源量的調(diào)查研究，對(duì)于水資源可持續(xù)利用，以支撐區(qū)域社會(huì)、經(jīng)濟(jì)、建設(shè)的健康可持續(xù)發(fā)展，具有著舉足輕重的重要意義。

我國(guó)在20世紀(jì)便已把水資源領(lǐng)域業(yè)務(wù)提到國(guó)家戰(zhàn)略的高度。1980年3月水利部下達(dá)《全國(guó)水資源調(diào)查評(píng)價(jià)工作重點(diǎn)》，部署開(kāi)展全國(guó)第一次水資源調(diào)查評(píng)價(jià)工作；2002年水利部、國(guó)家計(jì)委頒發(fā)《關(guān)于開(kāi)展全國(guó)水資源綜合規(guī)劃編制工作的通知》，要求在全國(guó)范圍內(nèi)開(kāi)展水資源調(diào)查評(píng)價(jià)工作；2017年4月，水利部召開(kāi)第三次全國(guó)水資源調(diào)查評(píng)價(jià)工作啟動(dòng)會(huì)議，依照水利部、國(guó)家發(fā)展改革委的《關(guān)于開(kāi)展第三次全國(guó)水資源調(diào)查評(píng)價(jià)工作的通知》，全面啟動(dòng)和部署第三次全國(guó)水資源調(diào)查評(píng)價(jià)工作。

在以往的水資源調(diào)查評(píng)價(jià)工作中，1985年鄭濯清提出的分項(xiàng)調(diào)查法是其中應(yīng)用較為廣泛的一種，張佑民（1991）、陸中央（2000）等均曾采用該方法對(duì)天然徑流量進(jìn)行還原計(jì)算。然而，在實(shí)際工作中，該方法所依賴的調(diào)查數(shù)據(jù)量較大、參數(shù)多，徑流還原實(shí)現(xiàn)難度較大。由于人類活動(dòng)對(duì)降雨徑流過(guò)程影響復(fù)雜，各地采用的徑流還原方法均存在一定局限。考慮到傳統(tǒng)的降雨徑流相關(guān)法在描述流域降雨徑流相關(guān)關(guān)系時(shí)，在全國(guó)范圍內(nèi)均具有良好的應(yīng)用效果，本研究嘗試借鑒該方法的產(chǎn)流計(jì)算原理，依據(jù)水資源工作的要求進(jìn)行改造，以構(gòu)建適合于徑流還原計(jì)算的模型。

1 模型介紹

1.1 模型概述

研究中提出的月尺度天然徑流還原計(jì)算模型（簡(jiǎn)稱月徑流模型）以傳統(tǒng)的經(jīng)驗(yàn)水文模型（API：Antecedent Precipitation IndexMethod）為基礎(chǔ)，參考流域特點(diǎn)，擬在考慮流域滯蓄以及下墊面水源涵養(yǎng)等因素的前提下，提出改進(jìn)的月尺度天然徑流量模擬方法，以用于逐月水文過(guò)程及水資源量模擬。該方法認(rèn)為，上一月的產(chǎn)流量一部分作為快速?gòu)搅髟诋?dāng)月全部流出；慢速部分（河網(wǎng)地下水）通過(guò)流域內(nèi)水庫(kù)等水利工程設(shè)施調(diào)蓄以及下墊面蓄滯作用影響下，在下一個(gè)月全部流出。

以月為計(jì)算步長(zhǎng)，依據(jù)上一個(gè)月降水量所在區(qū)間，分別進(jìn)行月徑流模擬計(jì)算。由于計(jì)算過(guò)程不考慮蒸發(fā)、水庫(kù)塘壩調(diào)蓄、工農(nóng)業(yè)取用水等因素，因此模型所提供徑流模擬結(jié)果即可視作是天然徑流還原量。月徑流模型以傳統(tǒng)降雨徑流相關(guān)法（API）為基礎(chǔ)，在模型構(gòu)建過(guò)程中需要考慮時(shí)間尺度變化所帶來(lái)的描述對(duì)象改變，同時(shí)也要保證實(shí)際使用中的可操作性，因此模型在繼承API模型的基本原理同時(shí)，在如下方面進(jìn)行了改造：

（1）以上一個(gè)月降水量而非前期影響雨量作為衡量流域濕潤(rùn)程度的指標(biāo)。API模型一般要求以日為計(jì)算步長(zhǎng)，提前15日開(kāi)始遞推估算前期雨量值。API運(yùn)算過(guò)程中所要求的計(jì)算步長(zhǎng)與本模型差別較大，且逐日降雨數(shù)據(jù)相對(duì)逐月降雨更難以獲得。因此本模型考慮采用上一個(gè)月的降水量代替API中的前期雨量值，并根據(jù)降水量所在區(qū)間范圍適配不同的降雨徑流曲線。

（2）僅考慮兩種不同濕潤(rùn)程度下的降雨徑流曲線。降雨徑流曲線本身即是基于統(tǒng)計(jì)方法的，原始數(shù)據(jù)的可靠性程度越高，該曲線所可以反映越多的水文物理過(guò)程的細(xì)節(jié)；反之，當(dāng)原始數(shù)據(jù)可靠性程度較低時(shí)，該曲線并不能夠從細(xì)節(jié)上描述對(duì)應(yīng)的產(chǎn)流過(guò)程。例如在天然徑流還原過(guò)程中，待分析的天然徑流量往往缺乏實(shí)測(cè)數(shù)據(jù)驗(yàn)證，因此，在研究中我們僅考慮干旱、濕潤(rùn)兩種土壤狀態(tài)下的降雨徑流關(guān)系。

（3）以兩段式的分段曲線描述降雨徑流關(guān)系。區(qū)分未蓄滿、蓄滿產(chǎn)流兩部分進(jìn)行降雨徑流模擬，在未蓄滿狀態(tài)下，降雨量越大，徑流深增加的幅度變大，P-R曲線斜率越大；當(dāng)流域蓄滿時(shí)，所有降雨可以轉(zhuǎn)換成流域徑流深，此時(shí)降雨徑流關(guān)系曲線的斜率穩(wěn)定為1。

（4）不考慮匯流過(guò)程。API模型中，產(chǎn)流與匯流過(guò)程是分離開(kāi)的，在計(jì)算出產(chǎn)流量之后需要配合采用單位線模擬流域匯流過(guò)程，以獲得流域出口斷面位置流量。在月尺度天然徑流還原工作中，我們所關(guān)注的并非流域出口位置流量，兼且由于流域內(nèi)閘壩調(diào)蓄、工農(nóng)業(yè)取水等數(shù)據(jù)難以精確獲得，因此只需關(guān)注降雨及對(duì)應(yīng)產(chǎn)流深的關(guān)系即可。

1.2 模型的構(gòu)建

搜集歷史記錄中月降雨量、還原徑流量序列，依據(jù)上月降雨量從小到大排序，認(rèn)為上月降雨的概率分布小于、不小于10%數(shù)據(jù)分別對(duì)應(yīng)干旱、濕潤(rùn)兩種土壤狀態(tài)，據(jù)此構(gòu)建兩種土壤濕潤(rùn)程度下的降雨、徑流數(shù)據(jù)集合。分別以兩段式的函數(shù)，以擬合精度最高為目標(biāo)，對(duì)以上兩組數(shù)據(jù)進(jìn)行擬合適線。本研究中，對(duì)于任意濕潤(rùn)條件下的降雨徑流曲線，都嘗試以兩段式的分段函數(shù)進(jìn)行描述，根據(jù)產(chǎn)流過(guò)程進(jìn)行劃分，各分段函數(shù)分別對(duì)應(yīng)未蓄滿、蓄滿兩個(gè)產(chǎn)流階段。

假設(shè)試驗(yàn)流域有N個(gè)雨量站，依據(jù)泰森多邊形算法將試驗(yàn)流域劃分為N個(gè)子流域，每個(gè)子流域與唯一的雨量站對(duì)應(yīng)，且權(quán)重值和為1。

式中，p、pi為月降水量（mm）；wi為權(quán)重項(xiàng)，取值范圍[0，1]；i用于指示雨量站或子流域序號(hào)。在已知降雨量p之后，便需要以擬合精度最高為優(yōu)化目標(biāo)，構(gòu)建合適的降雨徑流曲線模擬產(chǎn)流深。基于前人研究成果，未蓄滿狀態(tài)下各子流域還原徑流量可以用以下形式的方程估算。

式中，R為產(chǎn)流深，S為還原徑流量，t為權(quán)重值，代表滯留于流域內(nèi)的上月徑流量所占比例。將上月的慢速?gòu)搅髁考爱?dāng)月快速?gòu)搅髁恐鸵曌鳟?dāng)月還原凈流量S。在確定計(jì)算公式的形式之后，需要通過(guò)參數(shù)優(yōu)化或參照參數(shù)的物理意義進(jìn)行參數(shù)求解，本研究基于各參數(shù)物理、統(tǒng)計(jì)意義，采用SCE-UA算法進(jìn)行參數(shù)優(yōu)化，尋找各參數(shù)的最優(yōu)解。對(duì)于蓄滿狀態(tài)下流域產(chǎn)流量的計(jì)算，需要考慮與未蓄滿狀態(tài)下降雨徑流關(guān)系曲線光滑連接，即1）拐點(diǎn)處斜率均為1；2）拐點(diǎn)處產(chǎn)流量相等。

1.3 參數(shù)特征討論

在參數(shù)尋優(yōu)過(guò)程中，首先需要確定各參數(shù)的取值范圍。d與R的量級(jí)相近，當(dāng)月降水為零時(shí)，流域內(nèi)徑流深僅包括以往月份來(lái)得及排出的部分，因此d應(yīng)當(dāng)為正值。由于R-P顯然呈正相關(guān)關(guān)系，因此我們?nèi)藶樵O(shè)定p的系數(shù)a、b均大于零。一般來(lái)講降雨量越多，產(chǎn)流深越大，因此上式中R對(duì)p的一階偏導(dǎo)應(yīng)當(dāng)始終不小于零，即：

圖1 月尺度降雨徑流模型流程圖

圖2 月徑流模型原理示意圖

參照流域內(nèi)其他相關(guān)工作可知，a值一般在（0，1），b值一般不大于4，c值與R'量級(jí)相近。另，根據(jù)蓄滿產(chǎn)流原理，剛開(kāi)始產(chǎn)流時(shí)降雨與其所產(chǎn)徑流深的比值較小，并隨著降雨的進(jìn)行而逐漸增大，并且僅當(dāng)流域完全蓄滿之后，降雨能夠全部轉(zhuǎn)化為徑流，R'接近最大值1。因此，R'應(yīng)當(dāng)滿足1）在未蓄滿時(shí)，R'＜=1，由于函數(shù)第一、二段函數(shù)之間需要銜接平滑的特性及第二段函數(shù)的單調(diào)特性，這一項(xiàng)對(duì)于R是始終滿足的；2）R'的值應(yīng)當(dāng)隨著p值增大而增大，R對(duì)p的二階偏導(dǎo)數(shù)應(yīng)當(dāng)滿足如下條件：

上式中，p＞0，則p的乘冪必定滿足不小于零的條件；a、b均大于零，則b＞1。根據(jù)以上，考慮到流域內(nèi)還原徑流量的實(shí)際情況，可知各參數(shù)滿足如下條件：

在已知參數(shù)取值范圍之后，即可以依據(jù)SCE-UA工具進(jìn)行參數(shù)尋優(yōu)，獲取兩段式函數(shù)的參數(shù)值。SCE-UA算法是一種全局優(yōu)化算法，這種方法以信息共享和自然界生物演化規(guī)律的概念為基礎(chǔ)，是段青云等人在亞利桑那州大學(xué)發(fā)展的一種基于非線性單純形法的混合方法。SCE-UA算法被認(rèn)為是流域水文模型參數(shù)優(yōu)選中最有效的方法，在流域水文模型參數(shù)優(yōu)選中應(yīng)用十分廣泛。

2 試驗(yàn)流域介紹

在日照市22個(gè)主要河流中，選取典型試驗(yàn)流域進(jìn)行模擬試驗(yàn)。青峰嶺水庫(kù)位于沭河干流上游、日照市莒縣城西北30km處，壩址在棋山、洛河、茅埠三鄉(xiāng)交界處的棋山鄉(xiāng)蘆家岔河村。流域面積770km2，現(xiàn)狀總庫(kù)容4.101億m3，興利庫(kù)容2.687億m3。青峰嶺水庫(kù)流域?qū)倥瘻貛Ъ撅L(fēng)區(qū)大陸性氣候。多年平均氣溫12.1℃，多年平均年降水量754.2mm，其中汛期（6～9月）多年平均降水量557.1mm，占多年平均年降水量的74%。

表1 SCE-UA算法的參數(shù)設(shè)定表

表2 月降雨徑流模型參數(shù)表

青峰嶺水庫(kù)流域內(nèi)的雨量站點(diǎn)既有汛期站點(diǎn)也有常設(shè)站點(diǎn)，且各站點(diǎn)雨量觀測(cè)年限長(zhǎng)度差異較大。為保證建模所用資料的一致性，本研究剔除汛期雨量站，僅選用輝泉、馬站、沙溝水庫(kù)、陳家莊、青峰嶺水庫(kù)五處常設(shè)雨量站的資料作為降雨輸入，資料年限為1978～2015年。

還原徑流量成果基于分項(xiàng)調(diào)查評(píng)價(jià)法獲得，由于分項(xiàng)調(diào)查法已經(jīng)在全國(guó)范圍的水資源調(diào)查評(píng)價(jià)業(yè)務(wù)中獲得了廣泛應(yīng)用，其計(jì)算結(jié)果是合理的；且某地區(qū)的天然徑流量難以獲得實(shí)測(cè)觀測(cè)數(shù)據(jù)，因此，本研究將該還原徑流量視作參考依據(jù)，用于模型及其參數(shù)的可靠性評(píng)價(jià)。為區(qū)分采用分項(xiàng)調(diào)查法及本研究所用模型計(jì)算獲得的還原徑流量，將分項(xiàng)調(diào)查法結(jié)果稱為還原徑流量，本研究的計(jì)算結(jié)果成為計(jì)算徑流量。

3 模擬結(jié)果討論

3.1 參數(shù)優(yōu)化結(jié)果

應(yīng)用SCE-UA算法，分別對(duì)前期干旱、濕潤(rùn)兩種情況的月降雨—徑流數(shù)據(jù)進(jìn)行模型參數(shù)擬合。優(yōu)化所選目標(biāo)函數(shù)為確定性系數(shù)（NSE），NSE取值范圍在（-∞,1]，其取值越大則表明模擬精度越高。以青峰嶺流域濕潤(rùn)狀態(tài)下月降雨～徑流模型的參數(shù)率定過(guò)程為例，SCE-UA算法的參數(shù)設(shè)定及收斂過(guò)程示意如表1、圖3。

本次模擬試驗(yàn)，在分析1978～2015年數(shù)據(jù)的基礎(chǔ)上，基于SCE-UA算法優(yōu)化公式（2）中各項(xiàng)參數(shù)，獲得如下率定結(jié)果列于表2。

3.2 模擬結(jié)果分析

以相關(guān)圖展示還原與計(jì)算徑流深之間的差異性（圖4），各點(diǎn)距離45°線越遠(yuǎn)、點(diǎn)據(jù)越分散，表明模擬精度越差，模型、參數(shù)越不可靠；反之，當(dāng)圖中各點(diǎn)集中分布在45°線兩側(cè)時(shí)，模擬結(jié)果越好。

圖3 （歸一化后）濕潤(rùn)狀態(tài)下月降雨徑流模型參數(shù)優(yōu)化過(guò)程示意圖

圖4 還原—計(jì)算徑流量相關(guān)圖

分析圖中干旱狀態(tài)下，流域計(jì)算與還原徑流量的關(guān)系可知，干旱狀態(tài)土壤條件下，流域計(jì)算徑流量幾乎為水平直線，流域計(jì)算徑流量變化范圍很窄，說(shuō)明所采用模型基本上不能夠準(zhǔn)確反映出干旱狀態(tài)流域降雨—徑流關(guān)系，或者干旱狀態(tài)下，流域降雨徑流關(guān)系難以用單一的降雨徑流曲線進(jìn)行精確描述。濕潤(rùn)狀態(tài)下，還原—計(jì)算徑流量點(diǎn)據(jù)相對(duì)均勻的分布在45°線的兩側(cè)，反映出月模型在濕潤(rùn)土壤狀態(tài)下能夠較為準(zhǔn)確地描述試驗(yàn)流域的產(chǎn)流過(guò)程。兩種前期土壤狀態(tài)下的降雨徑流相關(guān)關(guān)系如圖5所示。

圖5 降雨徑流關(guān)系圖

圖6 2005～2015年逐月降雨、還原及計(jì)算徑流過(guò)程圖

圖中拐點(diǎn)位置以上均為45°直線，即當(dāng)月降雨超過(guò)拐點(diǎn)時(shí)，降雨可全部轉(zhuǎn)化為徑流深。對(duì)于青峰嶺區(qū)間流域來(lái)說(shuō)，所用資料中，各年份的月降雨量均未達(dá)到拐點(diǎn)處降雨量，如圖6，將逐年降雨及對(duì)應(yīng)的還原、模擬徑流量繪制在同一幅圖中，前期土壤狀態(tài)為干旱、濕潤(rùn)的兩種情況均包含在內(nèi)。限于篇幅，我們僅摘取2005～2015年期間降雨、徑流過(guò)程進(jìn)行展示：

如圖6所示，月降雨與還原、模擬徑流量的相關(guān)性強(qiáng)，月徑流模型模擬結(jié)果與還原徑流量結(jié)果較為一致，說(shuō)明以公式（2）為核心的月徑流模型能夠較為準(zhǔn)確地反映流域降雨、徑流關(guān)系，對(duì)于天然徑流量的還原計(jì)算具備參考價(jià)值。

4 結(jié)論

根據(jù)以上模擬結(jié)果可知：（1）不同的土壤狀態(tài)下，流域降雨徑流模擬結(jié)果差異明顯，因此根據(jù)流域前一月的干旱狀態(tài)，區(qū)分干旱、濕潤(rùn)兩種情況，分別采用不同的降雨徑流關(guān)系曲線進(jìn)行天然徑流模擬計(jì)算的做法，是合理的；（2）基于傳統(tǒng)API模型改進(jìn)而來(lái)的月尺度天然徑流還原模型，其計(jì)算成果與常規(guī)的分項(xiàng)調(diào)查法所得到計(jì)算結(jié)果基本一致，該模型適用于青峰嶺水庫(kù)流域的天然徑流量還原的計(jì)算。

本研究所提出的月尺度天然徑流還原計(jì)算模型，無(wú)需探究水資源利用過(guò)程中各分項(xiàng)水量，依賴于已經(jīng)搜集到的大量水文數(shù)據(jù)進(jìn)行建模計(jì)算，即可獲得具備較高可靠性模擬結(jié)果。區(qū)分不同前期土壤濕潤(rùn)程度分別進(jìn)行徑流還原、區(qū)分不同速度的徑流成分進(jìn)行合成運(yùn)算的思路，在本研究中已被證明是可靠的，以上徑流還原計(jì)算的思路具備參考價(jià)值，值得在更廣泛的流域范圍內(nèi)推廣使用■