微分形式新安江模型

DOI：10.14042/j.cnki.32.1309.2024.03.002

摘要：現(xiàn)有新安江模型數(shù)學(xué)上是代數(shù)方程并限于一階差分方法求解，不可避免存在數(shù)值誤差，探索數(shù)值求解誤差控制新途徑，對(duì)于提高模型計(jì)算精度具有重要意義。基于微分系統(tǒng)框架，識(shí)別新安江模型的狀態(tài)變量和通量，推導(dǎo)其控制方程和本構(gòu)方程，提出微分形式的新安江模型（ODE-XAJ），并采用四階顯式Runge Kutta法求解。數(shù)值實(shí)驗(yàn)結(jié)果顯示，以解析解為基準(zhǔn)，ODE-XAJ絕對(duì)誤差處于或小于10-4 mm量級(jí)，可實(shí)現(xiàn)對(duì)解析解的高階近似;以O(shè)DE-XAJ結(jié)果為基準(zhǔn)，按歸一化平均絕對(duì)誤差評(píng)價(jià)，現(xiàn)有新安江模型的數(shù)值求解誤差約為8.7%。典型流域應(yīng)用結(jié)果顯示，ODE-XAJ確定性系數(shù)提升0.02，洪量相對(duì)誤差降低4.3%。研究表明，ODE-XAJ理論上分離了模型的數(shù)學(xué)方程與具體解法，可有效控制數(shù)值求解誤差，提升模型模擬精度。

關(guān)鍵詞：新安江模型;數(shù)值誤差;微分形式;Runge Kutta法;屯溪流域

中圖分類號(hào)：P338

文獻(xiàn)標(biāo)志碼：A

文章編號(hào)：1001-6791（2024）03-0374-13

收稿日期：2023-11-18;網(wǎng)絡(luò)出版日期：2024-05-06

網(wǎng)絡(luò)出版地址：https：∥link.cnki.net/urlid/32.1309.P.20240430.1549.002

基金項(xiàng)目：國(guó)家自然科學(xué)基金資助項(xiàng)目（52379007）;水利部重大科技項(xiàng)目（SKR-2022032）

作者簡(jiǎn)介：梁忠民（1962—），男，遼寧鳳城人，教授，博士，主要從事水文水資源研究。E-mail：zmliang@hhu.edu.cn

通信作者：趙建飛，E-mail：jianfei.zhao@hhu.edu.cn

水文模型是對(duì)復(fù)雜自然水文規(guī)律的合理概化，是水文預(yù)報(bào)的主要技術(shù)手段［1-2］，一般可分為黑箱模型、集總式模型和分布式模型［3］，其中，集總式模型因其結(jié)構(gòu)簡(jiǎn)單、易于實(shí)現(xiàn)等優(yōu)點(diǎn)，得到了廣泛的研究與應(yīng)用，如中國(guó)的新安江模型等［3-5］。一般而言，集總式水文模型的數(shù)學(xué)方程為基于時(shí)段的代數(shù)方程，是以離散方式刻畫連續(xù)的水文過程，在建立模型的數(shù)學(xué)方程時(shí)，包含了采用一階精度的特定數(shù)值解法（如使用時(shí)段初值代替時(shí)段均值），如此不可避免地會(huì)引入數(shù)值求解誤差［1］。這種誤差不僅會(huì)降低模型模擬或預(yù)報(bào)精度，而且給模型參數(shù)確定帶來困難［6-7］，同時(shí)，數(shù)值求解帶來的誤差會(huì)隨降水量的增大而增大［8］，更增加了模型應(yīng)用的不確定性。

新安江模型于1973年提出［1-2］，模型包含蒸散發(fā)、產(chǎn)流、分水源和匯流4個(gè)模塊。自提出以來，模型不斷得到改進(jìn)，其中，蒸散發(fā)模塊由1層發(fā)展到3層形式，以解決久旱后實(shí)際蒸散發(fā)計(jì)算值偏小的問題;使用地表徑流、壤中流和地下徑流的三水源形式代替了地表徑流、地下徑流的二水源形式，使得模型適用于壤中流豐富的流域，由此形成了現(xiàn)有的新安江模型［9］。在此基礎(chǔ)上，目前對(duì)新安江模型的輸入、結(jié)構(gòu)和參數(shù)誤差的研究較多［10-12］，但對(duì)數(shù)值求解引入的固有誤差研究相對(duì)不足。已有研究［1，13］通過縮短時(shí)間步長(zhǎng)方式（根據(jù)凈雨按5 mm劃分子時(shí)段）對(duì)數(shù)值誤差進(jìn)行控制，但由于解析解的缺失［14］，無法對(duì)模型求解數(shù)值誤差的大小做出有效評(píng)判。

基于微分方程構(gòu)建水文模型可實(shí)現(xiàn)自然水文過程的連續(xù)描述，對(duì)微分方程進(jìn)行時(shí)間域的積分，可直接得到一些簡(jiǎn)單的水文模型（如Collie水量平衡模型［15］）或模型的子過程（如線性水庫匯流過程）的解析解;對(duì)于復(fù)雜的水文模型，即使無法得到其解析解，但當(dāng)將其以微分形式表達(dá)后，可實(shí)現(xiàn)模型數(shù)學(xué)方程和具體數(shù)值解法的分離，通過引入更高階精度的數(shù)值解法（如Runge Kutta法等）得到其理論解析解的高階近似。因此，研究水文模型的微分形式表達(dá)是控制模型數(shù)值求解誤差的有效途徑［16-17］。

本文以新安江模型（XAJ）為對(duì)象，識(shí)別微分系統(tǒng)框架下的模型狀態(tài)變量和通量，推導(dǎo)相應(yīng)的控制方程和本構(gòu)方程，建立微分形式的新安江模型（ODE-XAJ）。通過數(shù)值實(shí)驗(yàn)和實(shí)際流域驗(yàn)證，對(duì)比分析了ODE-XAJ與現(xiàn)有XAJ之間的差異，驗(yàn)證了ODE-XAJ的優(yōu)勢(shì)。研究成果可為控制集總式模型數(shù)值求解誤差、提高模型精度提供一個(gè)新途徑。

1" 微分形式新安江模型推導(dǎo)

1.1" 狀態(tài)變量及模型通量的識(shí)別

本文中，新安江模型［1］的蒸散發(fā)模塊采用3層蒸散發(fā)公式，產(chǎn)流計(jì)算使用考慮不透水面積影響的張力水蓄水容量曲線，分水源模塊為三水源形式，匯流模塊包含坡面與河網(wǎng)匯流，分別采用線性水庫和Nash單位線［18］。新安江模型邏輯結(jié)構(gòu)見圖1。

構(gòu)建微分系統(tǒng)框架下的水文模型，包括狀態(tài)變量和通量的識(shí)別、狀態(tài)變量控制方程及通量本構(gòu)方程的推導(dǎo)等環(huán)節(jié)［8］。狀態(tài)變量是表征流域狀態(tài)的物理量，其量綱與時(shí)間無關(guān);通量用于描述流域內(nèi)部或流域與外部之間的水分交換量，可以表示為狀態(tài)變量或其他通量的函數(shù)，其量綱與時(shí)間有關(guān);控制方程用于描述流域系統(tǒng)狀態(tài)變量隨時(shí)間的演化;本構(gòu)方程則用于建立未知通量與狀態(tài)變量或其他已知通量之間的關(guān)系。因此，通過識(shí)別模型狀態(tài)變量和通量，推導(dǎo)其相應(yīng)數(shù)學(xué)表達(dá)方程，即可構(gòu)建微分形式的水文模型。

通過對(duì)新安江模型中物理量的量綱分析，共確定了9個(gè)狀態(tài)變量和22個(gè)通量（表1）。其中，9個(gè)狀態(tài)變量包括：蒸散發(fā)模塊的Wu、Wl、Wd（流域上層、下層、深層土壤平均張力水蓄量），分水源模塊的S0（流域平均自由水蓄量），坡面匯流子模塊的Oi、Og（壤中流和地下徑流線性水庫蓄量），河網(wǎng)匯流子模塊的F1、F2、F3（Nash單位線第1、第2、第3個(gè)水庫的蓄量），其單位均為mm。

1.2" 控制方程及本構(gòu)方程的推導(dǎo)

一般而言，控制方程的左側(cè)為狀態(tài)變量的微分，右側(cè)包含1個(gè)或多個(gè)通量。逐一確定右側(cè)通量的本構(gòu)方程，使得微分系統(tǒng)閉合，從而形成微分形式水文模型的完整表達(dá)。限于篇幅，本文僅以新安江模型的狀態(tài)變量Wu為例，推導(dǎo)其控制方程及對(duì)應(yīng)通量的本構(gòu)方程。

對(duì)于Wu而言，其接收凈雨輸入，供給上層土壤蒸散發(fā)并形成產(chǎn)流，當(dāng)上層土壤蓄滿之后，多余水分補(bǔ)給下層土壤，其控制方程為

dWudt=Pn-R-Eu-Iu（1）

Pobs超出蒸散發(fā)能力部分形成Pn，不足的部分為En，對(duì)應(yīng)公式為：

Pn=max（Pobs-KeEobs，0）（2）

En=max（KeEobs-Pobs，0）（3）

式中：Ke為潛在蒸散發(fā)折算系數(shù)。

根據(jù)蓄滿產(chǎn)流理論［1］，采用考慮不透水面積影響的流域張力水蓄水容量曲線計(jì)算R，對(duì)應(yīng)曲線公式為

fw=Aimp+（1-Aimp）［1-（1-W/Wmm）b］（4）

式中：fw為流域蓄滿面積比例（徑流系數(shù)）;Aimp為不透水面積比例;W為流域內(nèi)單點(diǎn)張力水蓄量，mm;b為流域張力水蓄水容量曲線指數(shù);Wmm為流域內(nèi)單點(diǎn)最大蓄水容量，mm。對(duì)式（4）沿縱坐標(biāo)軸方向進(jìn)行積分（如圖1產(chǎn)流模塊所示），得到公式：

W0=∫aw0（1-fw）dW=（1-Aimp）Wmm1+b［1-（1-aw/Wmm）1+b］（5）

式中：W0為流域平均張力水蓄量，W0=Wu+Wl+Wd，mm;aw為此時(shí)流域內(nèi)單點(diǎn)張力水蓄量的最大值，mm。將aw=Wmm代入式（5），可得：

Wm=（1-Aimp）Wmm/（1+b）（6）

式中：Wm為流域平均張力水容量，mm。根據(jù)式（5），可推導(dǎo)aw關(guān)于W0的公式，具體為

aw=Wmm［1-（1-W0/Wm）1/（1+b）］（7）

令W=aw，將式（7）代入式（4），得到W0對(duì)應(yīng)的fw，具體為

fw=Aimp+（1-Aimp）［1-（1-W0/Wm）b/（1+b）］（8）

最終得到產(chǎn)流強(qiáng)度公式為

R=Pnfw=Pn{Aimp+（1-Aimp）［1-（1-W0/Wm）b/（1+b）］}（9）

根據(jù)3層蒸散發(fā)公式，當(dāng)Wugt;0時(shí)，上層土壤按照凈蒸發(fā)強(qiáng)度進(jìn)行蒸散發(fā)，否則沒有蒸散發(fā);當(dāng)Wugt;Wum后，超出水量均補(bǔ)給下層土壤，對(duì)應(yīng)公式為：

Eu=EnWugt;0

0Wu≤0（10）

Iu=Pn-R-EuWu≥Wum

0Wult;Wum（11）

上述公式中，式（1）為Wu的控制方程，式（2）、（9）、（10）和（11）為式（1）右側(cè)通量Pn、R、Eu和Iu的本構(gòu)方程。

同理，也可以推導(dǎo)出新安江模型其他狀態(tài)變量的控制方程及通量本構(gòu)方程，最終得到微分形式新安江模型，其狀態(tài)變量控制方程組見式（12），各通量本構(gòu)方程見表1。ODE-XAJ參數(shù)和取值范圍與現(xiàn)有新安江模型保持一致［13，19］，具體見表2。

ddt［Wu，Wl，Wd，S0，Oi，Og，F(xiàn)1，F(xiàn)2，F(xiàn)3］T=［Pn-R-Eu-Iu，Iu-El-Il，Il-Ed，Pn-（Rs+Ri+Rg）/

（fw-Aimp），Ri-Qi，Rg-Qg，Qt-Q1，Q1-Q2，Q2-Q3］T（12）

2" ODE-XAJ驗(yàn)證實(shí)驗(yàn)方案

由新安江模型原理可知，其數(shù)學(xué)模型由簡(jiǎn)單線性函數(shù)和冪函數(shù)的串并聯(lián)組合構(gòu)成，整體表現(xiàn)出高度的非線性。本文雖然推導(dǎo)出了新安江模型微分形式的數(shù)學(xué)方程，但從表1中可以看出，ODE-XAJ仍然無法得到任意初值及參數(shù)組合條件下的全局解析解（只能得到特定條件下的局部解析解）。然而，由于ODE-XAJ分離了現(xiàn)有新安江模型的數(shù)學(xué)方程和具體解法，意味著在無法直接獲取全局解析解的情況下，可使用高階精度方法（而非原有一階精度差分方法）進(jìn)行數(shù)值求解，理論上可以獲得全局解析解的高階近似。因此，本文建議采用具有四階精度的顯式Runge Kutta法（RK4）求解ODE-XAJ。為說明該做法（ODE-XAJ+RK4）的有效性，本文設(shè)計(jì)了如下3種驗(yàn)證方案：

（1） 理想數(shù)值實(shí)驗(yàn)。設(shè)定特定初始值和參數(shù)條件，推導(dǎo)出蒸散發(fā)、產(chǎn)流、分水源、匯流4個(gè)獨(dú)立模塊的通量解析解;生成模型輸入樣本（降水、蒸發(fā)）及參數(shù)值，采用ODE-XAJ+RK4得到模型輸出的數(shù)值解;對(duì)比數(shù)值解與解析解，從而評(píng)定數(shù)值解的近似程度。

（2） 對(duì)比數(shù)值實(shí)驗(yàn)。采用相同初始值和模型參數(shù)，以實(shí)際流域的降水和蒸發(fā)數(shù)據(jù)為輸入，同時(shí)使用ODE-XAJ+RK4與現(xiàn)有新安江模型得到相應(yīng)輸出，對(duì)比2類模型的計(jì)算精度和效率。

（3） 實(shí)際流域驗(yàn)證。使用實(shí)際流域的降水、蒸發(fā)和流量數(shù)據(jù)，分別率定ODE-XAJ及現(xiàn)有新安江模型的最優(yōu)參數(shù)，并進(jìn)行流量過程模擬，對(duì)比2類模型對(duì)實(shí)測(cè)流量數(shù)據(jù)的擬合效果。

2.1" 理想數(shù)值實(shí)驗(yàn)

理想數(shù)值實(shí)驗(yàn)包含4個(gè)子實(shí)驗(yàn)，設(shè)置見表3。子實(shí)驗(yàn)均使用了相同的1 000組參數(shù)，其通過對(duì)稱拉丁超立方體方法（SLH）［20］從新安江模型參數(shù)及其范圍（表2）采樣1 000次得到。

子實(shí)驗(yàn)1使用20組恒定強(qiáng)度的Pobs和Eobs，其總量均分別為0和（1-c）Wlm，序列長(zhǎng)度為1～20。參數(shù)c、Wlm和Wdm來源于SLH方法采樣生成的1 000組參數(shù)，參數(shù)Ke固定為1。記T為時(shí)段的開始時(shí)間，Δt為時(shí)段長(zhǎng)。T～T+Δt時(shí)段內(nèi)，下層土壤蒸散發(fā)累積量（El，*）的解析解為

El，*（T）=Wlmexp（-EobsT/Wlm）-Wlmexp［-Eobs（T+Δt）/Wlm］（13）

子實(shí)驗(yàn)2使用5組恒定強(qiáng)度和15組變化強(qiáng)度的Pobs和Eobs，其總量均分別為Wmm和0，其中5組恒定強(qiáng)度的序列長(zhǎng)度分別為1、10、20、50和100，15種變化強(qiáng)度的輸入通過隨機(jī)數(shù)生成。為便于推導(dǎo)產(chǎn)流累積量（R*）、透水面積上地表徑流累積量（Rs，*）的解析解，參數(shù)Ki和Kg固定為0以保證只生成地表徑流，參數(shù)Sm設(shè)置為Wmm/（1+ex），其余參數(shù)來源于SLH方法采樣生成的1 000組參數(shù)。T～T+Δt時(shí)段內(nèi)，R*和Rs，*的解析解分別為：

R*（T）=Pobs（T）Δt+W0-Wm+Wm{1-［aw（T）+Pobs（T）Δt］/Wmm}1+b（14）

Rs，*（T）=-（1-Aimp）Sm1-aw（T）Sm（1+ex）1+ex-1-aw（T）+Pobs（T）ΔtSm（1+ex）1+ex

+（1-Aimp）Pobs（T）Δt-Wm1-aw（T）Wmm1+b-1-aw（T）+Pobs（T）ΔtWmm1+b

+Wmm（1-Aimp）1+ex+b1-aw（T）Wmm1+ex+b-1-aw（T）+Pobs（T）ΔtWmm1+ex+b（15）

子實(shí)驗(yàn)3使用20組恒定強(qiáng)度的Ri和Rg，其總量均為1 mm，輸入類型與子實(shí)驗(yàn)1類同。參數(shù)Ci和Cg來源于SLH方法采樣生成的1 000組參數(shù)。T～T+Δt時(shí)段內(nèi)，壤中流及地下徑流線性水庫出流時(shí)段累積量（Qi，*和Qg，*）的解析解分別為：

Qi，*（T）=RiΔt-RiΔt{exp［（TlnCi+ΔtlnCi）/Δt］-exp［（TlnCi）/Δt］}/lnCi（16）

Qg，*（T）=RgΔt-RgΔt{exp［（TlnCg+ΔtlnCg）/Δt］-exp［（TlnCg）/Δt］}/lnCg（17）

子實(shí)驗(yàn)4使用5組恒定強(qiáng)度和15組變化強(qiáng)度的Qt，其總量均為1 mm，輸入類型與子實(shí)驗(yàn)2類同。參數(shù)Kf來源于SLH方法采樣生成的1 000組參數(shù)。T時(shí)刻Q3的解析解為

Q3（T）=Qt（T）exp-TKfexp-ΔtKf∑3i=1［（T-Δt）/Kf］3-i（3-i）！-∑3i=1（T/Kf）3-i（3-i）?。?8）

式中：為卷積符號(hào)。

使用四階顯式Runge Kutta法對(duì)ODE-XAJ進(jìn)行數(shù)值求解，采用Shampine自適應(yīng)時(shí)間步長(zhǎng)技術(shù)［21］，以進(jìn)一步提升求解精度。該技術(shù)采用歸一化后的局部估計(jì)誤差動(dòng)態(tài)調(diào)整時(shí)間步長(zhǎng)，局部估計(jì)誤差歸一化公式為

escaled=|y∧n+1-yn+1|/［a+rmax（yn，yn-1）］（19）

式中：escaled為歸一化后的局部估計(jì)誤差;y為狀態(tài)變量向量，y=［Wu，Wl，Wd，S0，Oi，Og，F(xiàn)1，F(xiàn)2，F(xiàn)3］T;yn、yn+1和yn-1為y在tn、tn+1和tn-1的狀態(tài)變量值，tn和tn+1分別為計(jì)算時(shí)段的開始和結(jié)束時(shí)間，tn-1為上個(gè)計(jì)算時(shí)段的開始時(shí)間;y∧n+1為使用高階精度公式計(jì)算的y在tn+1時(shí)狀態(tài)變量值;|y∧n+1-yn+1|為局地誤差估計(jì)值;a為絕對(duì)誤差閾值，mm;r為相對(duì)誤差閾值，本文中a和r均取10-4。

選取平均絕對(duì)誤差對(duì)數(shù)值解法的精度進(jìn)行評(píng)估，對(duì)應(yīng)的方程為

EMA=1N∑Ni=1|z（i）-z∧（i）|（20）

式中：z為通量或狀態(tài)變量的計(jì)算值;z∧為通量或狀態(tài)變量的解析值（或解析解的高階近似值）;N為時(shí)間序列長(zhǎng)度。

2.2" 對(duì)比數(shù)值實(shí)驗(yàn)

基于SLH方法采樣生成的1 000組參數(shù)和相同的零初始條件，以屯溪流域的實(shí)測(cè)流域面平均日降水和日蒸發(fā)數(shù)據(jù)作為模型輸入，對(duì)比ODE-XAJ+RK4與現(xiàn)有新安江模型在計(jì)算精度和計(jì)算效率上的差異。

屯溪流域是典型濕潤(rùn)區(qū)流域，年平均降水量約1 750 mm，流域面積為2 670 km2。流域日尺度實(shí)測(cè)降水、屯溪水文站逐日流量資料來源于水文年鑒，流域面平均日降水通過算術(shù)平均法計(jì)算得到，實(shí)測(cè)蒸發(fā)資料的來源為中國(guó)地面氣候資料日值數(shù)據(jù)集（V3.0），資料年限為2007—2018年。屯溪流域地理位置及站點(diǎn)分布如圖2所示。

現(xiàn)有新安江模型進(jìn)一步分為基于或不基于縮短時(shí)間步長(zhǎng)法2類?？s短時(shí)間步長(zhǎng)法（STS）由趙人俊提出，用以控制新安江模型中由差分解法引入的數(shù)值誤差［1，13］。STS的核心為根據(jù)實(shí)測(cè)降水和蒸發(fā)的時(shí)段累積量將計(jì)算時(shí)段劃分為若干正整數(shù)（G）個(gè)子時(shí)段，對(duì)應(yīng)公式為

G=ceil（|Pobs，*-KeEobs，*|/M）+1（21）

式中：ceil為向下取整函數(shù);Pobs，*和Eobs，*分別為實(shí)測(cè)降水和蒸發(fā)的時(shí)段累積量，mm;M為子時(shí)段入流累積量的最大值，mm。由于實(shí)測(cè)降水和蒸發(fā)在不同計(jì)算時(shí)段存在區(qū)別，對(duì)應(yīng)G的具體取值在不同時(shí)段有所變化。子時(shí)段內(nèi)實(shí)測(cè)降水和蒸發(fā)的累積量為原計(jì)算時(shí)段的1/G，使用的公式與原計(jì)算時(shí)段基本保持一致。原計(jì)算時(shí)段的初值作為第1個(gè)子時(shí)段的初值，上個(gè)子時(shí)段的末值作為下個(gè)子時(shí)段的初值，最后1個(gè)子時(shí)段的末值作為原下個(gè)計(jì)算時(shí)段的初值。參數(shù)Ki、Kg、Ci和Cg由于時(shí)段長(zhǎng)度發(fā)生變化需要進(jìn)行轉(zhuǎn)換，子時(shí)段對(duì)應(yīng)參數(shù)值分別以Ki，G、Kg，G、Ci，G和Cg，G表示，其余模型參數(shù)保持不變，對(duì)應(yīng)轉(zhuǎn)換公式為：

Ki，G=Ki［1-（1-Ki+Kg）1/G］/（Ki+Kg）（22）

Kg，G=Kg［1-（1-Ki+Kg）1/G］/（Ki+Kg）（23）

Ci，G=（Ci）1/G（24）

Cg，G=（Cg）1/G（25）

對(duì)比數(shù)值實(shí)驗(yàn)使用ODE-XAJ+RK4、不基于縮短時(shí)間步長(zhǎng)法的新安江模型（XAJ）和基于縮短時(shí)間步長(zhǎng)的新安江模型（XAJ+STS）3類模型，所有模型的初始值均設(shè)置為0。STS的M分別取5、0.5、0.05和0.000 5 mm，用以分析M的具體取值對(duì)XAJ+STS計(jì)算精度和計(jì)算效率的影響。

基于相同的硬件平臺(tái)（CPU為AMD Ryzen 9 5950X）運(yùn)行模型，根據(jù)運(yùn)行時(shí)間評(píng)估模型的計(jì)算效率。使用EMA和歸一化平均絕對(duì)誤差（ENMA）評(píng)估數(shù)值求解誤差的整體情況，ENMA的公式為

ENMA=∑Ni=1|z（i）-z∧（i）|/∑Ni=1z∧（i）×100%（26）

2.3" 實(shí)際流域驗(yàn)證

采用XAJ和ODE-XAJ對(duì)屯溪流域日徑流進(jìn)行模擬。使用資料為日尺度的流域面平均降水、蒸發(fā)和屯溪水文站的逐日平均流量，資料年限為2007—2018年，其中預(yù)熱期為2007年、率定期為2008—2015年、驗(yàn)證期為2016—2018年。XAJ和ODE-XAJ的初始值均設(shè)置為0。根據(jù)確定性系數(shù)（CD）［22］和洪量相對(duì)誤差（EFVR）分別評(píng)定XAJ和ODE-XAJ的模擬精度。采用SCE-UA方法［23］分別率定XAJ和ODE-XAJ的模型參數(shù)，搜索次數(shù)為2 000次，目標(biāo)函數(shù)為|1-CD|最小，分別重復(fù)10次搜索過程取最優(yōu)參數(shù)。

3" 結(jié)果及分析

3.1" 理想數(shù)值實(shí)驗(yàn)結(jié)果及分析

采用理想的輸入、參數(shù)和初始條件組合，使用四階顯式Runge Kutta法數(shù)值求解ODE-XAJ，以輸出通量的解析解為基準(zhǔn)，計(jì)算輸出通量數(shù)值解的EMA。由于使用20種輸入類型和1 000組參數(shù)，模型對(duì)應(yīng)模擬20 000次，對(duì)每個(gè)輸出通量而言，可計(jì)算20 000個(gè)EMA，其分布具體見圖3。從圖3可看出，各輸出通量的EMA均值整體處于或小于10-4量級(jí);Rs，*的EMA最大值大于其他通量，但小于5.0×10-3 mm。理想數(shù)值實(shí)驗(yàn)結(jié)果表明，四階顯式Runge Kutta法在求解ODE-XAJ（ODE-XAJ+RK4）時(shí)有很高的精度，可以很好地逼近解析解。

3.2" 對(duì)比數(shù)值實(shí)驗(yàn)結(jié)果及分析

以屯溪流域的實(shí)際降水、蒸發(fā)為輸入，基于1 000組模型參數(shù)，使用ODE-XAJ、XAJ和XAJ+STS 3種模型進(jìn)行模擬。以O(shè)DE-XAJ結(jié)果為基準(zhǔn)，計(jì)算其余模型模擬的狀態(tài)變量和通量時(shí)段累積量的EMA。由于使用了1 000組模型參數(shù)，對(duì)XAJ和XAJ+STS的特定狀態(tài)變量或通量時(shí)段累積量，均可計(jì)算出1 000個(gè)EMA，具體見表4。

從表4可看出，隨M取值減小，XAJ+STS狀態(tài)變量和通量的EMA最大值和平均值均表現(xiàn)出下降的趨勢(shì);在M取值為0.000 5 mm時(shí)，XAJ+STS無論在狀態(tài)變量還是通量時(shí)段累積量，都與ODE-XAJ對(duì)應(yīng)結(jié)果十分接近，其EMA平均值處于或小于10-3量級(jí)。XAJ狀態(tài)變量和通量的EMA均較大，其中S0的EMA最大值和平均值分別可達(dá)到1.77和1.47 mm，河網(wǎng)出流時(shí)段累積量（Q*）的EMA最大值和平均值分別為44.30和9.52 m3/s，Q*的ENMA最大值和平均值分別為37.7%和8.7%。對(duì)Q*而言，采用M取值為5 mm的STS方法后，XAJ模型W0、Et，*（實(shí)際蒸散發(fā)強(qiáng)度時(shí)段累積量）、R*、S0、Rs，*、Ri，*（壤中流時(shí)段累積量）、Rg，*（地下徑流時(shí)段累積量）、Qi，*、Qg，*和Q*的EMA平均值分別下降27.4%、30.9%、26.2%、77.6%、53.9%、38.6%、38.4%、71.8%、77.0%和63.3%，表明STS方法可改善新安江模型的數(shù)值求解誤差，但EMA均值整體處于10-1的量級(jí)，Q*的ENMA均值為3.2%，其精度仍需進(jìn)一步的提升。

對(duì)比數(shù)值實(shí)驗(yàn)基于相同的1 000組參數(shù)，對(duì)每個(gè)模型運(yùn)行1 000次并記錄其運(yùn)行時(shí)間，使用運(yùn)行時(shí)間平均值來分析不同模型的計(jì)算效率，具體見表5。從表5可看出，ODE-XAJ的平均運(yùn)行時(shí)間與M取值為5 mm的XAJ+STS較為接近;在M取值分別為0.5和0.05 mm時(shí)，XAJ+STS的平均運(yùn)行時(shí)間分別是ODE-XAJ的1.6和17.0倍;在M取值為0.000 5 mm時(shí)，XAJ+STS的平均運(yùn)行時(shí)間遠(yuǎn)超ODE-XAJ。

綜合而言，在原有新安江模型的基礎(chǔ)上使用縮短時(shí)間步長(zhǎng)法，雖然可提升模型計(jì)算精度，但其精度的提升意味著模型計(jì)算效率的快速下降，這是由差分解法的一階精度特性決定的。ODE-XAJ在理論上實(shí)現(xiàn)了模型數(shù)學(xué)方程與具體解法的分離，而不限于原一階精度差分解法，因而可以通過引入高階精度的方法（如Runge Kutta）求解，達(dá)到計(jì)算精度和計(jì)算效率的有效平衡。

3.3" 實(shí)際流域驗(yàn)證結(jié)果與分析

使用屯溪流域的降水、蒸發(fā)和流量數(shù)據(jù)，分別率定了ODE-XAJ和XAJ的模型參數(shù)并進(jìn)行日徑流模擬，模型參數(shù)見表6，對(duì)應(yīng)模擬結(jié)果見表7。

從表7可以看出，XAJ除2018年洪量相對(duì)誤差EFVR為20.8%以外，其余年份EFVR均在±20%以內(nèi)，其EFVR的絕對(duì)值平均為8.7%（率定期為6.6%、驗(yàn)證期為14.1%）。在徑流過程評(píng)定方面，XAJ率定期和驗(yàn)證期的確定性系數(shù)CD平均值分別為0.93和0.86，所有年份的CD均大于0.75，率定期和驗(yàn)證期的CD最大分別可達(dá)0.95和0.90。上述結(jié)果表明，XAJ對(duì)屯溪流域日徑流有較好的模擬結(jié)果。

無論是率定期和驗(yàn)證期，ODE-XAJ的EFVR均在±20%以內(nèi)，EFVR的絕對(duì)值平均為4.4%（率定期為2.9%、驗(yàn)證期為8.5%），所有年份的EFVR都優(yōu)于XAJ。在徑流過程評(píng)定方面，ODE-XAJ率定期和驗(yàn)證期的CD平均值分別為0.94和0.92，所有年份的CD均大于0.85，率定期和驗(yàn)證期的CD最大分別可達(dá)0.96和0.95;ODE-XAJ的CD僅在2011年和2014年略低于XAJ，ODE-XAJ的CD平均值和最大值在率定期和驗(yàn)證期均高于XAJ。上述結(jié)果表明，在屯溪流域日徑流模擬中，ODE-XAJ優(yōu)于XAJ。作為示例，圖4提供了XAJ和ODE-XAJ的模型模擬結(jié)果對(duì)比。

4" 結(jié)" 論

本文提出了微分形式新安江模型，通過數(shù)值實(shí)驗(yàn)和典型流域應(yīng)用，對(duì)比分析了微分形式新安江模型和現(xiàn)有新安江模型在計(jì)算精度、效率和實(shí)際模擬中的差異。主要結(jié)論如下：

（1） 在微分系統(tǒng)框架下，識(shí)別新安江模型的狀態(tài)變量和模型通量，推導(dǎo)其控制方程和本構(gòu)方程，建立了微分形式新安江模型，豐富了新安江模型理論。

（2） 數(shù)值實(shí)驗(yàn)結(jié)果表明，四階顯式Runge-Kutta法可以高階近似微分形式新安江模型的解析解。與采用縮短時(shí)間步長(zhǎng)法的新安江模型相比，微分形式新安江模型可更好地平衡模型計(jì)算精度和效率。典型濕潤(rùn)流域的應(yīng)用結(jié)果表明，微分形式新安江模型可進(jìn)一步提升模擬精度。

（3） 本文使用的方法可直接推廣到其他水文模型，具有較好的普適性。微分系統(tǒng)能夠自然地描述具有不同演變速率、動(dòng)態(tài)變化且相互作用的多個(gè)過程，其意味著微分形式的水文模型更易與其他學(xué)科模型（如大氣、生態(tài)模型等）交叉耦合，因此具有廣闊的應(yīng)用前景。

參考文獻(xiàn)：

［1］趙人俊.流域水文模擬：新安江模型與陜北模型［M］.北京：水利電力出版社，1984.（ZHAO R J.Hydrological simulation of watershed：Xin′anjiang model and Northern Shaanxi model［M］.Beijing：Water Resources and Electric Power Press，1984.（in Chinese））

［2］張建云.中國(guó)水文預(yù)報(bào)技術(shù)發(fā)展的回顧與思考［J］.水科學(xué)進(jìn)展，2010，21（4）：435-443.（ZHANG J Y.Review and reflection on China′s hydrological forecasting techniques［J］.Advances in Water Science，2010，21（4）：435-443.（in Chinese））

［3］芮孝芳，凌哲，劉寧寧，等.新安江模型的起源及對(duì)其進(jìn)一步發(fā)展的建議［J］.水利水電科技進(jìn)展，2012，32（4）：1-5.（RUI X F，LING Z，LIU N N，et al.Origin of Xin′anjiang model and its further development［J］.Advances in Science and Technology of Water Resources，2012，32（4）：1-5.（in Chinese））

［4］劉金濤，宋慧卿，張行南，等.新安江模型理論研究的進(jìn)展與探討［J］.水文，2014，34（1）：1-6.（LIU J T，SONG H Q，ZHANG X N，et al.A discussion on advances in theories of Xin′anjiang model［J］.Journal of China Hydrology，2014，34（1）：1-6.（in Chinese））

［5］陸旻皎.新安江模型研究的回顧和展望［J］.水利學(xué)報(bào)，2021，52（4）：432-441.（LU M J.Recent and future studies of the Xin′anjiang model［J］.Journal of Hydraulic Engineering，2021，52（4）：432-441.（in Chinese））

［6］CLARK M P，KAVETSKI D.Ancient numerical daemons of conceptual hydrological modeling：1：fidelity and efficiency of time stepping schemes［J］.Water Resources Research，2010，46：W10510.

［7］WOLDEGIORGIS B T，BAULCH H，WHEATER H，et al.Impacts of uncontrolled operator splitting methods on parameter identification，prediction uncertainty，and subsurface flux representation in conceptual hydrological models［J］.Water Resources Research，2023，59（7）：e2022WR033250.

［8］la FOLLETTE P T，TEULING A J，ADDOR N，et al.Numerical daemons of hydrological models are summoned by extreme precipitation［J］.Hydrology and Earth System Sciences，2021，25（10）：5425-5446.

［9］王厥謀，張恭肅，李玉瑤，等.趙人俊水文預(yù)報(bào)文集［M］.北京：水利電力出版社，1994：110-115.（WANG J M，ZHANG G S，LI Y Y，et al.Proceedings of ZHAO Renjun′s hydrological forecasting［M］.Beijing：Water Resources and Electric Power Press，1994：110-115.（in Chinese））

［10］關(guān)鐵生，鮑振鑫，賀瑞敏，等.無資料地區(qū)水文模型參數(shù)移植不確定性分析［J］.水科學(xué)進(jìn)展，2023，34（5）：660-672.（GUAN T S，BAO Z X，HE R M，et al.Uncertainties of model parameters regionalization in ungauged basins［J］.Advances in Water Science，2023，34（5）：660-672.（in Chinese））

［11］WANG J，ZHUO L，HAN D W，et al.Hydrological model adaptability to rainfall inputs of varied quality［J］.Water Resources Research，2023，59（2）：e2022WR032484.

［12］KNOBEN W J M，F(xiàn)REER J E，F(xiàn)OWLER K J A，et al.Modular Assessment of Rainfall-Runoff Models Toolbox （MARRMoT） v1.2：an open-source，extendable framework providing implementations of 46 conceptual hydrologic models as continuous state-space formulations［J］.Geoscientific Model Development，2019，12（6）：2463-2480.

［13］包為民.水文預(yù)報(bào)［M］.5版.北京：中國(guó)水利水電出版社，2017.（BAO W M.Hydrologic forecast［M］.5th ed.Beijing：China Water amp; Power Press，2017.（in Chinese））

［14］BISHT G，RILEY W J.Development and verification of a numerical library for solving global terrestrial multiphysics problems［J］.Journal of Advances in Modeling Earth Systems，2019，11（6）：1516-1542.

［15］JOTHITYANGKOON C，SIVAPALAN M，F(xiàn)ARMER D L.Process controls of water balance variability in a large semi-arid catchment：downward approach to hydrological model development［J］.Journal of Hydrology，2001，254（1/2/3/4）：174-198.

［16］SANTOS L，THIREL G，PERRIN C.Continuous state-space representation of a bucket-type rainfall-runoff model：a case study with the GR4 model using state-space GR4 （version 1.0）［J］.Geoscientific Model Development，2018，11（4）：1591-1605.

［17］GOUDARZI S，MILLEDGE D，HOLDEN J.A generalized multistep dynamic （GMD） TOPMODEL［J］.Water Resources Research，2023，59（1）：e2022WR032198.

［18］SINGH V P.Estimation of parameters of a uniformly nonlinear surface runoff model［J］.Hydrology Research，1977，8（1）：33-46.

［19］陸桂華，酈建強(qiáng)，楊曉華.水文模型參數(shù)優(yōu)選遺傳算法的應(yīng)用［J］.水利學(xué)報(bào)，2004，35（2）：50-56.（LU G H，LI J Q，YANG X H.Application of genetic algorithms to parameter optimization of hydrology model［J］.Journal of Hydraulic Engineering，2004，35（2）：50-56.（in Chinese））

［20］GONG W，DUAN Q，LI J，et al.An intercomparison of sampling methods for uncertainty quantification of environmental dynamic models［J］.Journal of Environmental Informatics，2016，28（1）：11-24.

［21］SHAMPINE L F.Solving ODEs and DDEs with residual control［J］.Applied Numerical Mathematics，2005，52（1）：113-127.

［22］水文情報(bào)預(yù)報(bào)規(guī)范：GB/T 22482—2008［S］.北京：中國(guó)標(biāo)準(zhǔn)出版社，2009.（Standard for hydrological information and hydrological forecasting：GB/T 22482—2008［S］.Beijing：Standards Press of China，2009.（in Chinese））

［23］DUAN Q Y，SOROOSHIAN S，GUPTA V.Effective and efficient global optimization for conceptual rainfall-runoff models［J］.Water Resources Research，1992，28（4）：1015-1031.

Differential-form Xin′anjiang model

The study is financially supported by the National Natural Science Foundation of China （No.52379007） and the Major Science and Technology Projects of the Ministry of Water Resources，China （No.SKR-2022032）.

LIANG Zhongmin1，ZHAO Jianfei1，DUAN Ya′nan2，HUANG Jialu1，LI Binquan1，WANG Jun1，HU Yiming1

（1. College of Hydrology and Water Resources，Hohai University，Nanjing 210098，China;

2. Jiangsu Province Water

Conservancy Engineering Sci-tech Consulting Co.，Ltd.，Nanjing 210029，China）

Abstract：The existing Xin′anjiang model is mathematically represented by algebraic equations，and its solution is limited by first-order accuracy finite difference methods，inevitably resulting in numerical errors.Therefore，exploring new approaches to controlling numerical errors is crucial for increasing the computational accuracy of the model.Within the framework of differential systems，the state variables and fluxes of the Xin′anjiang model are identified，and the corresponding control equations and constitutive equations are then derived.A differential-form Xin′anjiang model （ODE-XAJ） is proposed and solved using the fourth-order explicit Runge-Kutta method.The results of the numerical experiments show that，when benchmarked against the analytical solution，the absolute error of ODE-XAJ is on the order of 10-4 mm or below，enabling a high-order approximation of the analytical solution.With the results of ODE-XAJ as the benchmark，the numerical error of the existing Xin′anjiang model is evaluated to be approximately 8.7% based on the normalized mean absolute error.The application results for a typical watershed show that the determination coefficient of ODE-XAJ increased by 0.02 and the relative error of flood volume decreased by 4.3%.It is concluded that ODE-XAJ theoretically separates the mathematical equations of the model from specific numerical methods，which can effectively control the numerical errors and improve the accuracy of model simulations.

Key words：Xin′anjiang model;numerical error;differential form;Runge-Kutta method;Tunxi watershed