輸水渠系水動(dòng)力數(shù)字孿生模型糙率估計(jì)方法

管光華 劉王嘉儀 陳曉楠 史良勝

摘要：為提高水動(dòng)力數(shù)字孿生模型校正環(huán)節(jié)中糙率估計(jì)的實(shí)時(shí)性和精細(xì)化，考慮糙率值在渠道縱向上的空間變異性，提出基于水力半徑變化和估計(jì)精度分段估計(jì)糙率的思路;基于渠道分段，提出獨(dú)立估計(jì)法和聯(lián)合估計(jì)法2種不同估計(jì)框架。基于有限的觀測(cè)水位，在框架內(nèi)應(yīng)用集合卡爾曼濾波算法，在線估計(jì)各渠段的糙率值。結(jié)果表明：相比未分段時(shí)，2種估計(jì)方法可提高模型精度20%～50%，獨(dú)立估計(jì)法誤差累積小，適合復(fù)雜渠系;而聯(lián)合估計(jì)法適用于觀測(cè)量缺失的簡(jiǎn)單渠道。研究成果可服務(wù)于水動(dòng)力數(shù)字孿生模型的參數(shù)估計(jì)和變量更新，為建設(shè)數(shù)字孿生水網(wǎng)提供參考。

關(guān)鍵詞：糙率;輸水渠系;數(shù)字孿生;南水北調(diào);集合卡爾曼濾波;分段估計(jì)

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

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

文章編號(hào)：1001-6791（2023）06-0901-12

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

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

基金項(xiàng)目：國(guó)家自然科學(xué)基金資助項(xiàng)目（51979202;51879199）

作者簡(jiǎn)介：管光華（1979—），男，江蘇阜寧人，副教授，博士，主要從事灌排自動(dòng)化及輸水調(diào)度管理研究。

Email：GGH@whu.edu.cn

隨著物聯(lián)網(wǎng)、大數(shù)據(jù)等信息技術(shù)的發(fā)展，數(shù)字孿生技術(shù)在水利領(lǐng)域的應(yīng)用為輸水渠系智能管理提供了新思路［1］。數(shù)字孿生技術(shù)可使得水動(dòng)力模型的狀態(tài)變量在線更新，從而對(duì)渠道系統(tǒng)的實(shí)時(shí)狀態(tài)進(jìn)行高精度的監(jiān)測(cè)和預(yù)測(cè)［2］。基于上述特性，研究指出，實(shí)時(shí)校正是數(shù)字孿生模型中必不可少的環(huán)節(jié)，尤其是模型參數(shù)的估計(jì)［3］。糙率是一維水動(dòng)力模型中的重要參數(shù)，所以構(gòu)建和校正水動(dòng)力數(shù)字孿生模型時(shí)，輸水明渠中糙率實(shí)時(shí)準(zhǔn)確的估計(jì)十分重要。

基于觀測(cè)手段和智能算法的發(fā)展，國(guó)內(nèi)外大量研究將實(shí)測(cè)數(shù)據(jù)與智能算法相結(jié)合以實(shí)時(shí)更新糙率值。李光熾等［4］用卡爾曼濾波求解糙率，提高仿真水位精度;Lai等［5］應(yīng)用集合卡爾曼濾波（EnKF）算法將糙率與水力變量都視作狀態(tài)量同步進(jìn)行更新，考慮觀測(cè)誤差的不確定性，進(jìn)一步提升模型精度。類似地，擴(kuò)展卡爾曼濾波、神經(jīng)網(wǎng)絡(luò)也被用于糙率的實(shí)時(shí)估計(jì)和校正［6-7］。徐興亞等［8］采用粒子濾波方法同化觀測(cè)水位，有效校正糙率以提高實(shí)時(shí)洪水預(yù)報(bào)的精度；王麟等［9］通過(guò)神經(jīng)網(wǎng)絡(luò)建立糙率與水力要素的關(guān)系，從而有效預(yù)測(cè)人工明渠的糙率變化。這些算法可實(shí)現(xiàn)在線校正，且結(jié)果準(zhǔn)確性較高，能克服傳統(tǒng)擬合和反演方法計(jì)算效率低、數(shù)據(jù)量需求大的問(wèn)題。然而，大部分研究均假定在2個(gè)觀測(cè)點(diǎn)位之間的渠段內(nèi)，糙率沿縱向是恒定值。這個(gè)假定偏離了糙率會(huì)因環(huán)境而變化的情況，從而導(dǎo)致水力模型的仿真誤差［10］。針對(duì)這個(gè)問(wèn)題，陳一帆等［11］以空間分布平滑性作為糙率反演的條件，可同時(shí)提高計(jì)算的穩(wěn)定性和精度;高學(xué)平等［12］對(duì)河道橫斷面進(jìn)行分區(qū)，并計(jì)算等效糙率;Attari等［13-14］將渠道隨機(jī)分為若干段，利用迭代尋優(yōu)的方法確定各渠段的長(zhǎng)度和糙率值，進(jìn)一步得到糙率的沿程分布。但是，上述研究多采用隨機(jī)或假定分段，并未結(jié)合影響糙率的水力因素。

基于以上問(wèn)題，本文針對(duì)實(shí)時(shí)的糙率分段估計(jì)，提出新的分段思路，同時(shí)考慮影響糙率的水力因素和估計(jì)結(jié)果的數(shù)值穩(wěn)定性;為實(shí)現(xiàn)數(shù)字孿生校正環(huán)節(jié)的實(shí)時(shí)性，高效同化觀測(cè)數(shù)據(jù)，提出2種采用EnKF的糙率估計(jì)框架;以不同方式結(jié)合未測(cè)量斷面的水位估計(jì)和各渠段的糙率估計(jì)，分別命名為獨(dú)立估計(jì)法和聯(lián)合估計(jì)法;以南水北調(diào)中線工程京石段為例，分析2種框架各自的效果和適用性。其中，聯(lián)合估計(jì)法參考文獻(xiàn)［15］，將水位與糙率的估計(jì)過(guò)程進(jìn)行嵌套。

1 渠道水動(dòng)力數(shù)字孿生模型

如圖1所示，渠道系統(tǒng)的水動(dòng)力數(shù)字孿生模型應(yīng)包含4個(gè)階段［16］，分別為數(shù)據(jù)預(yù)處理、仿真預(yù)熱、模型校正和仿真預(yù)測(cè)。本文針對(duì)模型校正階段，將水位和糙率都視為系統(tǒng)狀態(tài)量，提出2種新的分段糙率估計(jì)框架，提供準(zhǔn)確在線的參數(shù)估計(jì)方法；在仿真預(yù)測(cè)階段，對(duì)比水力變量的預(yù)測(cè)精度，檢驗(yàn)糙率分段估計(jì)方法的有效性。

1.1 數(shù)據(jù)預(yù)處理

數(shù)據(jù)預(yù)處理針對(duì)實(shí)時(shí)的觀測(cè)數(shù)據(jù)，包括識(shí)別并刪除異常值、插值填充缺失值和去除噪聲3個(gè)步驟?？山Y(jié)合數(shù)理統(tǒng)計(jì)和水力特征曲線擬合的方法檢測(cè)異常值。當(dāng)前成熟的去噪方法較多，類似有濾波器、時(shí)間序列分解等等，可根據(jù)實(shí)際情況進(jìn)行選擇。

1.2 仿真預(yù)熱

模型的初始狀態(tài)會(huì)影響后續(xù)階段的仿真結(jié)果。為了降低模型前期的不穩(wěn)定性，提高初始的仿真精度，需要在水動(dòng)力數(shù)字孿生模型中設(shè)置預(yù)熱階段。在預(yù)熱階段，模型與渠道系統(tǒng)同步運(yùn)行模擬，但不記錄仿真結(jié)果且不調(diào)整模型，提高后續(xù)階段仿真的可靠性。

1.3 模型校正

本文將糙率和水位視作狀態(tài)量，并采用EnKF進(jìn)行實(shí)時(shí)更新：水位估計(jì)的輸入量為觀測(cè)斷面的水位，輸出量是未測(cè)量斷面的水位;糙率估計(jì)的輸入量為斷面間水位差，輸出量為該段的糙率。EnKF每個(gè)時(shí)間步長(zhǎng)更新一次狀態(tài)量，時(shí)間步長(zhǎng)等于傳感器采集步長(zhǎng)［17］。水位和糙率估計(jì)對(duì)應(yīng)的關(guān)鍵變量和控制方程如表1所示。渠道系統(tǒng)的狀態(tài)空間方程由一維圣維南（Saint-Venant，SV）方程的離散格式推導(dǎo)而來(lái)，由四點(diǎn)Preissmann格式求解［18］，且忽略側(cè)向入流的影響。

基于渠段沿程（m+1）斷面的水位值，SV動(dòng)量方程可估計(jì)分段糙率。為便于建立狀態(tài)空間方程，取糙率值的平方作為狀態(tài)量。方程組包含m段渠段建立的m個(gè)動(dòng)量方程，如式（2）所示。

Nk+1=Nk+λ2，k

f（Zk+1）=Zk+1，j+1-Zk+1k+1，j=φk+1Nk+1+βk+1+ξ2，k+1（2）

式中：Nk+1=（nk+1）2，n為渠段糙率;βk+1，j=-ΔxjgΔteAk+1，j（Qk+1，j-Qk，j）-αQ2Ak+1，j+1-αQ2Ak+1，jgAk+1，j;φ為觀測(cè)矩陣，φk+1，j=-（|Q|Q）k+1，j（Ak+1，j）2（Rk+1，j）4/3Δxj;α為動(dòng)能系數(shù)，本文取α=1;A為斷面的過(guò)水面積，m2;g為重力加速度，m2/s;R為水力半徑，m;λ2為糙率狀態(tài)預(yù)測(cè)中添加的白噪聲;f（Z）為觀測(cè)變量，此處為相鄰兩斷面的水位差，m;β為觀測(cè)方程中的輸入項(xiàng);ξ2為觀測(cè)噪聲。

對(duì)于不同斷面處的流量輸入值，本文僅考慮渠段沿程取水的影響。已知上游來(lái)流量和沿程各分水口的取水流量，各斷面處的流量可由下式計(jì)算：

Qj+1=Qj-Qout（3）

式中：Qout為各斷面的分水流量，m3/s，如果斷面不處于分水口處，則Qout=0。

觀測(cè)數(shù)據(jù)較為充足時(shí)，分水流量值為分水口處的流量觀測(cè)值。在分水口處流量計(jì)不足的情況下，則根據(jù)渠段進(jìn)、出口流量假定取水流量。本文將渠段內(nèi)分水總流量設(shè)置為渠段上、下游閘門處過(guò)閘流量之差，過(guò)閘流量為閘門處流量計(jì)采集的觀測(cè)數(shù)據(jù)。當(dāng)渠段內(nèi)有多個(gè)分水口時(shí)，取水流量按分水口設(shè)計(jì)流量的比例進(jìn)行分配。

1.4 仿真預(yù)測(cè)

在模型仿真階段，取所有更新時(shí)刻糙率估計(jì)值的平均值作為參數(shù)輸入一維渠道水力模型，對(duì)比不同參數(shù)下的模型誤差。采用Preissmann四點(diǎn)差分隱格式的離散圣維南方程求解非恒定流。模型的上游邊界條件為流量—水位關(guān)系（Q=f（Z）），下游邊界條件為動(dòng)態(tài)水位邊界（Z（t）），其值均為對(duì)應(yīng)的觀測(cè)數(shù)據(jù)。該水力模型以各時(shí)刻渠池的下游水位為輸入值來(lái)校正模型仿真，預(yù)測(cè)上游的各斷面水位。模型的校正間隔等于觀測(cè)數(shù)據(jù)采集的時(shí)間間隔，且大于仿真時(shí)間步長(zhǎng)，運(yùn)算邏輯如下式所示：

（Zu，…，Zj，…，Zm-1）=F（Zd，q，n，Δtp，ΔT，Δs）（4）

式中：F（·）為自校正水力模型的運(yùn)行邏輯;Zu和Zd分別為各渠池的上、下游水位，m;Zj為渠池內(nèi)各斷面的水位，m;Zm-1為渠池內(nèi)最后一個(gè)待估計(jì)斷面的水位，m;q為分水口的取水流量，m3/s;Δtp為預(yù)測(cè)段的仿真時(shí)間步長(zhǎng)，s;ΔT為自校正時(shí)間間隔，s;Δs為仿真空間步長(zhǎng)，m。

2 糙率分段估計(jì)方法

2.1 渠道分段方法

糙率的變化體現(xiàn)在過(guò)流斷面水力半徑的變化，與流量和水深的變化相關(guān)［19］。在人工渠道中，糙率值可能在2種節(jié)點(diǎn)發(fā)生變化：斷面尺寸變化處和取水點(diǎn)。另外，分段估計(jì)糙率以斷面間的水位差作為觀測(cè)量，過(guò)小的水位差會(huì)導(dǎo)致數(shù)值計(jì)算不穩(wěn)定以及估計(jì)值失真;過(guò)大的水位差則表示兩斷面相隔較遠(yuǎn)，估計(jì)的糙率分布偏離實(shí)際情況。渠道中斷面水位差與斷面間距相關(guān)，所以需要探究合適的分段長(zhǎng)度。結(jié)合以上2點(diǎn)，本文提出同時(shí)考慮渠道節(jié)點(diǎn)和適宜分段長(zhǎng)度劃分渠段，整體流程如下。

首先根據(jù)糙率估計(jì)值的方差確定適合分段長(zhǎng)度。在已知上、下游觀測(cè)水位的條件下，用未分段時(shí)整個(gè)渠段的糙率值推算渠道各斷面的初始水位。整段糙率值可通過(guò)式（2）和EnKF算法得到。由于觀測(cè)量和輸入量是隨時(shí)間變化的，對(duì)應(yīng)的水位差和糙率值是時(shí)變量。通過(guò)下式計(jì)算參數(shù)估計(jì)時(shí)期內(nèi)糙率值的方差，以此確定分段長(zhǎng)度，合適的水位差對(duì)應(yīng)的糙率估計(jì)結(jié)果方差較?。?/p>

V=1Mn∑kTk=1（nk-n）2（5）

式中：Mn為校正段內(nèi)糙率估計(jì)值的總數(shù);kT為估計(jì)期的計(jì)算時(shí)刻總數(shù);nk為第k步的糙率估計(jì)值;n為估計(jì)期內(nèi)的糙率平均值。

尋優(yōu)過(guò)程中相鄰循環(huán)之間設(shè)置的水位差關(guān)系和循環(huán)的終止條件如式（6）—式（7）所示：

ΔZr=ΔZr-1+Δz（6）

（1） 若分段區(qū)間兩側(cè)任一斷面位于節(jié)制閘處，則刪除斷面變化和分水口處的斷面，保留基于適宜分段長(zhǎng)度而設(shè)置的斷面;

（2） 若分段區(qū)間兩側(cè)任一斷面因適宜長(zhǎng)度而設(shè)置，則刪除因適宜長(zhǎng)度而設(shè)的斷面，保留斷面變化或分水口處的分段斷面;

（3） 若分段區(qū)間兩側(cè)斷面均處于斷面變化或分水口處，則刪除上游側(cè)斷面，保留下游側(cè)斷面。

基于以上分段步驟，渠道可分為長(zhǎng)度不等的若干渠段。一般地，模型校正期應(yīng)處于同一工況下，故應(yīng)用以上分段方案時(shí)，在糙率估計(jì)的過(guò)程中分段結(jié)果保持不變。當(dāng)運(yùn)行工況發(fā)生明顯變化時(shí)，應(yīng)重新進(jìn)行模型校正，分段方案不變，由于觀測(cè)數(shù)據(jù)和糙率方差可能有明顯不同，分段結(jié)果也會(huì)發(fā)生相應(yīng)變化。

2.2 糙率估計(jì)框架

受傳感器數(shù)量的限制，大部分水位傳感器僅設(shè)置在每個(gè)渠段的段首和段尾斷面（節(jié)制閘附近），渠段中間斷面的水位缺少直觀測(cè)量。因此，在估計(jì)分段糙率值時(shí)，還需要估計(jì)未測(cè)量、處于中間斷面的水位值?；诖诵枨?，本文提出了2種不同的糙率估計(jì)框架，分別為獨(dú)立估計(jì)法和聯(lián)合估計(jì)法，2種方法的具體步驟如圖2所示。為體現(xiàn)糙率估計(jì)值的物理意義，在2種框架中均設(shè)置約束糙率變化范圍為0.012～0.025［20］。

2.2.1 獨(dú)立估計(jì)法

對(duì)于水位、糙率這2個(gè)狀態(tài)量的初始狀態(tài)，本文提出初值的計(jì)算方式。初始糙率值為通過(guò)渠段上游端、下游端觀測(cè)水位估計(jì)的整段糙率值，此時(shí)各段糙率值一致;初始水位則為由初始糙率值根據(jù)明渠恒定非均勻流方程推求的各斷面水位。得到初始狀態(tài)后，使用式（1）—式（3）和EnKF算法，估計(jì)各時(shí)刻的斷面水位，而后再估計(jì)各時(shí)刻的分段糙率。水位和糙率值的估計(jì)是分別完成的，水位估計(jì)時(shí)僅引入下游斷面水位觀測(cè)值對(duì)水位進(jìn)行更新估計(jì)，糙率更新值不參與水位計(jì)算。

2.2.2 聯(lián)合估計(jì)法

聯(lián)合估計(jì)法獲得初始狀態(tài)的步驟與獨(dú)立估計(jì)方法一致，且同樣基于式（1）—式（2）和EnKF更新狀態(tài)量。二者的不同在于后續(xù)步驟：每個(gè)時(shí)間步長(zhǎng)內(nèi)，在水位估計(jì)時(shí)用糙率推求的水面線代替觀測(cè)水位，由各分段斷面的水位推算各分段糙率。該時(shí)刻的糙率估計(jì)值繼而影響下一時(shí)刻的水位計(jì)算。水位與糙率估計(jì)嵌套循環(huán)，該過(guò)程中糙率的更新值參與水位計(jì)算。

2.3 評(píng)價(jià)指標(biāo)

本文將獨(dú)立估計(jì)法、聯(lián)合估計(jì)法與未分段時(shí)EnKF估計(jì)的糙率值用于水位仿真，比較3種糙率值對(duì)應(yīng)的模型仿真精度。本文選用均方根誤差（ERMS）、絕對(duì)平均誤差（EMA）和相對(duì)平均誤差（EMAP）衡量預(yù)測(cè)值與觀測(cè)值的偏離程度，這3個(gè)指標(biāo)常用來(lái)表示誤差的離散程度和平均水平，與未分段相比，以EMA的降幅作為本文方法改善模型精度的幅度［21］。另外，使用相關(guān)系數(shù)（R）衡量估計(jì)水位與觀測(cè)水位的一致性，R越接近1效果越好［21］。各指標(biāo)計(jì)算公式如下：

ERMS=1MZ∑MZi=1（Zi，o-Zi，s）2（10）

EMA=1MZ∑MZi=1|Zi，o-Zi，s|（11）

EMAP=1MZ∑MZi=1Zi，o-Zi，sZi，o×100%（12）

R=∑MZi=1（Zi，o-Zo）（Zi，s-Zs）∑MZi=1（Zi，o-Zo）2∑MZi=1（Zi，s-Zs）2（13）

式中：MZ為水位值數(shù)據(jù)點(diǎn)的總數(shù);Zi，o和Zi，s分別為第i個(gè)水位觀測(cè)值和模擬值，m;Zo和Zs分別為水位觀測(cè)值和模擬值的平均值，m。

3 南水北調(diào)中線工程應(yīng)用實(shí)例

3.1 案例概況

本文選取南水北調(diào)中線工程京石段內(nèi)2個(gè)渠池（漠道溝節(jié)制閘—唐河節(jié)制閘段、唐河節(jié)制閘—放水河節(jié)制閘段）為研究案例，下文稱為渠池1和渠池2，渠底坡降均為0.04‰。渠池1長(zhǎng)為8.419 km，斷面尺寸的沿程變化較少，起點(diǎn)、終點(diǎn)底高程分別為66.321、65.998 m;而渠池2長(zhǎng)為25.567 km，斷面尺寸變化較多，起點(diǎn)、終點(diǎn)底高程分別為65.985、64.937 m。2個(gè)渠池在長(zhǎng)度和復(fù)雜程度上相差較大，以說(shuō)明方法在不同規(guī)模渠段中的應(yīng)用效果。案例具有以下特征：過(guò)水?dāng)嗝鏋樘菪?，渠底坡度平緩，無(wú)逆坡段;各渠池的蓄量較大，水位波動(dòng)不顯著;各渠池水面線趨于水平，即呈現(xiàn)下游水位高于上游水位的現(xiàn)象。

渠池1僅在下游處設(shè)有1個(gè)分水口，渠池2在中游和下游處共2個(gè)分水口。由于缺乏水流量的實(shí)時(shí)記載，渠池2中各分水流量需要按照分水口的設(shè)計(jì)流量進(jìn)行分配。本文將局部水頭損失對(duì)水位的影響納入糙率值中，且選取的案例不含大型渡槽、隧道、倒虹吸等大型過(guò)水建筑物，以降低建筑物帶來(lái)的誤差影響。渠道路線和渠池示意見圖3，按上游至下游方向的斷面相關(guān)參數(shù)變化見表2。

3.2 觀測(cè)數(shù)據(jù)

分別選取兩渠段在2個(gè)不同時(shí)期內(nèi)連續(xù)10 d的水位、流量觀測(cè)數(shù)據(jù)作為模型的輸入值和對(duì)比值，傳感器的采樣間隔為2 h。將10 d的數(shù)據(jù)進(jìn)行如下劃分：第1天（0～24 h）為模型的預(yù)熱階段;第2～4天（24～96 h）為模型校正階段，使用2種框架估計(jì)糙率值;第5—10天（96～240 h）為仿真預(yù)測(cè)階段，以估計(jì)期內(nèi)的糙率平均值輸入模型，模擬水位并和觀測(cè)數(shù)據(jù)進(jìn)行對(duì)比。2個(gè)時(shí)期分別命名為時(shí)段A和時(shí)段B：時(shí)段A的數(shù)據(jù)采集于2018年4月，平均水深約為4.0 m，平均流量約為90 m3/s，可視為正常流量規(guī)模的輸水時(shí)期;時(shí)段B的數(shù)據(jù)采集于2019年10月，平均水深約為4.5 m，平均流量約為115 m3/s，為調(diào)水工程的高流量時(shí)期。以上調(diào)度過(guò)程中的觀測(cè)數(shù)據(jù)均由南水北調(diào)中線建管局總調(diào)中心提供。對(duì)比2份數(shù)據(jù)樣本，可說(shuō)明本文提出的糙率估計(jì)方法在不同工況下的適用性。

3.3 仿真設(shè)置

仿真通過(guò)基于Matlab語(yǔ)言設(shè)計(jì)的輸水渠道系統(tǒng)運(yùn)行仿真與控制軟件實(shí)現(xiàn)［22］，利用四點(diǎn)差分隱格式求解圣維南方程組仿真渠道一維非恒定流的水動(dòng)力過(guò)程。參數(shù)估計(jì)階段的更新時(shí)間步長(zhǎng)等于傳感器采集間隔，皆為2 h;模型仿真階段的時(shí)間步長(zhǎng)為20 min，下游水位數(shù)據(jù)輸入的時(shí)間間隔為2 h，總仿真時(shí)長(zhǎng)根據(jù)數(shù)據(jù)樣本設(shè)定，為240 h。2個(gè)階段的空間離散步長(zhǎng)皆為100 m。

4 結(jié)果分析及討論

根據(jù)2.1節(jié)的渠道分段方法，探究得到案例適合的分段長(zhǎng)度在3 000～4 000 m之間，該區(qū)間可適應(yīng)案例工程中的不同渠池。渠池1長(zhǎng)度較短，故選擇3 000 m為分段參考;渠池2較長(zhǎng)，則選擇4 000 m為分段參考。綜合考慮合適的分段長(zhǎng)度和水力半徑變化后，渠池1被分為3段，分別為3 000、3 000和2 419 m;渠池2被分為7段，分別為4 000、4 000、4 200、3 000、2 500、4 000和3 867 m。本節(jié)從流量數(shù)據(jù)影響、糙率估計(jì)值和模型仿真精度方面，對(duì)比2個(gè)框架的效果和適用性。

4.1 流量數(shù)據(jù)影響對(duì)比

流量數(shù)據(jù)作為水位、糙率估計(jì)的輸入項(xiàng)，會(huì)顯著影響估計(jì)結(jié)果。然而，現(xiàn)實(shí)中流量數(shù)據(jù)可能存在缺失或精度低的問(wèn)題，所以本節(jié)選取有中、下游分水口的渠池2在時(shí)段A內(nèi)的數(shù)據(jù)，分析2種估計(jì)方法受觀測(cè)數(shù)據(jù)影響的程度。

由于僅在渠池的上、下游斷面（節(jié)制閘附近）設(shè)有水位傳感器，故選取觀測(cè)斷面的相鄰分段斷面進(jìn)行分析。渠池2共8個(gè)斷面，本節(jié)選取上游觀測(cè)斷面的相鄰分段斷面，用2號(hào)斷面表示;下游觀測(cè)斷面的相鄰分段斷面，用7號(hào)斷面表示。從圖4可看出，獨(dú)立估計(jì)法得到的水位波動(dòng)呈現(xiàn)變幅較大且不穩(wěn)定的現(xiàn)象，水位估計(jì)值偏離真實(shí)狀態(tài)；而運(yùn)用聯(lián)合估計(jì)法時(shí)，估計(jì)水位的波動(dòng)更加平緩，變化趨勢(shì)貼近相鄰斷面的實(shí)測(cè)水位。聯(lián)合估計(jì)法的相關(guān)系數(shù)則明顯大于獨(dú)立估計(jì)法（R1表示獨(dú)立估計(jì)法的擬合度，R2表示聯(lián)合估計(jì)法的擬合度），且較接近1，說(shuō)明水位與觀測(cè)值變化更一致。水位估計(jì)的準(zhǔn)確性和穩(wěn)定性會(huì)直接影響糙率估計(jì)結(jié)果。由此可看出，獨(dú)立估計(jì)法受數(shù)據(jù)影響明顯更大，聯(lián)合估計(jì)法能夠在一定程度內(nèi)隔離流量數(shù)據(jù)的不利影響。

4.2 糙率估計(jì)結(jié)果對(duì)比

針對(duì)渠池2的各段糙率值進(jìn)行分析，從時(shí)間分布、空間分布的角度對(duì)比2種估計(jì)框架的糙率估計(jì)結(jié)果。中、下游2個(gè)分水口分別處于渠池2的第4段和第7段。

4.2.1 時(shí)間分布

表3展示了在估計(jì)期內(nèi)各渠段糙率估計(jì)值的方差，方差越大則表示各時(shí)刻的估計(jì)值相差較大，反之則表示估計(jì)值越穩(wěn)定。一般來(lái)說(shuō)，短期內(nèi)糙率值變化幅度較小，所以較穩(wěn)定的結(jié)果可信度較高。由表3可知，在時(shí)段A，獨(dú)立估計(jì)法的糙率方差基本小于聯(lián)合估計(jì)法；在時(shí)段B中，二者估計(jì)的糙率方差相差不大，僅在4號(hào)渠段相差較大，且獨(dú)立估計(jì)法的方差明顯小于聯(lián)合估計(jì)法。2個(gè)時(shí)期內(nèi)都出現(xiàn)4號(hào)渠段糙率方差較大的情況，其中時(shí)段A中2種方法估計(jì)得到的方差幾乎一致，原因可能是4號(hào)渠段內(nèi)的分水口進(jìn)行取水，渠段的上、下游斷面流量作為輸入項(xiàng)，其非恒定變化對(duì)估計(jì)結(jié)果影響較大。從糙率時(shí)間分布的方差來(lái)看，獨(dú)立估計(jì)法的估計(jì)值總體上方差更小，在穩(wěn)定性上更具優(yōu)勢(shì)。

4.2.2 空間分布

渠池2內(nèi)各渠段的糙率估計(jì)值見表3?？梢娫跁r(shí)段A內(nèi)，各渠段糙率值分布范圍為0.012～0.020;而在時(shí)段B內(nèi)，同一渠池內(nèi)各段糙率值的分布范圍更小，平均值更大，分布范圍為0.015～0.020。隨著時(shí)間的推移和環(huán)境、流態(tài)的變化，例如藻類生長(zhǎng)或泥沙淤積，同一區(qū)域可能會(huì)出現(xiàn)糙率增大的情況［23］。并且參考相關(guān)文獻(xiàn)［20］，糙率估計(jì)值在大型調(diào)水工程中糙率的合理分布范圍內(nèi)。從空間分布的范圍來(lái)看，2個(gè)方法在同一時(shí)期內(nèi)的估計(jì)結(jié)果分布較為接近，僅在個(gè)別含有分水口的渠段中有所差別，說(shuō)明2個(gè)方法皆較為合理，可使用在實(shí)際渠道中。

4.3 模型精度對(duì)比

表4展示了模型仿真的上游水位誤差。在2個(gè)時(shí)期的2個(gè)渠池中，本文提出2種方法的EMA為2～4 cm，EMAP大多在1%以內(nèi)。相較于無(wú)分段時(shí)EnKF的估計(jì)值，誤差的降幅大約為1～3 cm，說(shuō)明估計(jì)方法較為明顯地降低了參數(shù)不準(zhǔn)確帶來(lái)的模型誤差。另外，聯(lián)合估計(jì)法的ERMS普遍高于獨(dú)立估計(jì)法，可知獨(dú)立估計(jì)法誤差分布區(qū)間更小，預(yù)測(cè)更穩(wěn)定。表4對(duì)比了2種方法分別改善模型精度的幅度，總體達(dá)到20%～50%。在不同場(chǎng)景下，方法的性能表現(xiàn)并未出現(xiàn)明顯的優(yōu)劣傾向，可見本文提出的方法適用于南水北調(diào)中線工程中不同長(zhǎng)度規(guī)模、不同流量規(guī)模的渠池。在大部分工況下，獨(dú)立估計(jì)法的改善幅度高于聯(lián)合估計(jì)法，在渠池1、時(shí)段B和渠池2、時(shí)段A這2組高出10%左右，所以獨(dú)立估計(jì)法的仿真精度略高。原因可能是聯(lián)合估計(jì)法中，水位估計(jì)與糙率估計(jì)嵌套，二者的誤差相互影響和累計(jì)導(dǎo)致最終估計(jì)值對(duì)應(yīng)的模型精度降低。

4.4 敏感性分析

選取渠池1在時(shí)段A的數(shù)據(jù)樣本，保持下游觀測(cè)水位不變，人為修改上游觀測(cè)量，改變水位差并計(jì)算各分段糙率值的變化幅度。表5反映出2種估計(jì)方法在同一種水位變幅下，糙率值變化無(wú)明顯差異。當(dāng)水位變幅絕對(duì)值為0.02 m時(shí)，各段糙率值變化約為0.000 4～0.000 6;當(dāng)水位變幅絕對(duì)值為0.05 m時(shí)，糙率值變化約為0.001 2～0.001 5。水位值與實(shí)際觀測(cè)值偏離越大，糙率變化值越大，這符合水力規(guī)律。各段糙率值變化幅度較為一致，且范圍較小，說(shuō)明本文提出的2種估計(jì)方法對(duì)水位輸入值敏感性不高，表現(xiàn)較為穩(wěn)定。

5 結(jié)? 論

針對(duì)水動(dòng)力模型的實(shí)時(shí)校正問(wèn)題，本文提出了2種分段估計(jì)糙率的方法：獨(dú)立估計(jì)法和聯(lián)合估計(jì)法。以南水北調(diào)中線的長(zhǎng)距離輸水系統(tǒng)為例，進(jìn)行糙率分段估計(jì)的研究。主要結(jié)論如下：

（1） 本文提出的2個(gè)分段估計(jì)糙率方法可有效、實(shí)時(shí)地估計(jì)渠池中間未測(cè)量斷面的水位和各渠段的糙率值，使得模型同步貼合系統(tǒng)的實(shí)際狀態(tài)。處于合理范圍內(nèi)的糙率估計(jì)值體現(xiàn)該方法適用于類似南水北調(diào)中線的長(zhǎng)距離輸水系統(tǒng)。

（2） 相較于單一的經(jīng)驗(yàn)值，使用EnKF的分段估計(jì)方法可有效提高一維水力模型的精度，改善幅度在20%～50%之間。且對(duì)于不同長(zhǎng)度、不同流量的渠道，精度都有明顯的提升效果，說(shuō)明方法對(duì)南水北調(diào)中線等大型調(diào)水工程具有一定適用性。

（3） 本文對(duì)比所提2種估計(jì)方法，二者優(yōu)勢(shì)各有不同。獨(dú)立估計(jì)法計(jì)算過(guò)程更穩(wěn)定，誤差累積程度更小，適用于復(fù)雜多渠池系統(tǒng)并行估計(jì)的情況，適合大型輸水渠系;聯(lián)合估計(jì)法受流量數(shù)據(jù)影響較小，水位估計(jì)值更符合觀測(cè)水位的變化趨勢(shì)，適用于觀測(cè)數(shù)據(jù)缺失或質(zhì)量較差的簡(jiǎn)單渠道系統(tǒng)。

本文提出的分段估計(jì)糙率法是對(duì)模型精細(xì)化的初步探索，根據(jù)監(jiān)測(cè)數(shù)據(jù)實(shí)時(shí)估計(jì)參數(shù)、推斷水情，可為水動(dòng)力數(shù)字孿生模型的校正提供思路和方法。對(duì)于分段糙率估計(jì)問(wèn)題，仍然面臨著在有限觀測(cè)條件下保證真實(shí)性和準(zhǔn)確性的挑戰(zhàn)。

參考文獻(xiàn)：

［1］冶運(yùn)濤，蔣云鐘，梁犁麗，等.數(shù)字孿生流域：未來(lái)流域治理管理的新基建新范式［J］.水科學(xué)進(jìn)展，2022，33（5）：683-704.（YE Y T，JIANG Y Z，LIANG L L，et al.Digital twin watershed：new infrastructure and new paradigm of future watershed governance and management［J］.Advances in Water Science，2022，33（5）：683-704.（in Chinese））

［2］黃艷，喻杉，羅斌，等.面向流域水工程防災(zāi)聯(lián)合智能調(diào)度的數(shù)字孿生長(zhǎng)江探索［J］.水利學(xué)報(bào)，2022，53（3）：253-269.（HUANG Y，YU S，LUO B，et al.Development of the digital twin Changjiang River with the pilot system of joint and intelligent regulation of water projects for flood management［J］.Journal of Hydraulic Engineering，2022，53（3）：253-269.（in Chinese））

［3］PILLAI U P A，PINARDI N，ALESSANDRI J，et al.A Digital Twin modelling framework for the assessment of seagrass nature based solutions against storm surges［J］.Science of the Total Environment，2022，847：157603.

［4］李光熾，周晶晏，張貴壽.用卡爾曼濾波求解河道糙率參數(shù)反問(wèn)題［J］.河海大學(xué)學(xué)報(bào)（自然科學(xué)版），2003，31（5）：490-493.（LI G C，ZHOU J Y，ZHANG G S.Inverse analysis of river channel roughness parameter by use of Kalman filter［J］.Journal of Hohai University（Natural Sciences），2003，31（5）：490-493.（in Chinese））

［5］LAI R X，F(xiàn)ANG H W，HE G J，et al.Dual state-parameter optimal estimation of one-dimensional open channel model using ensemble Kalman filter［J］.Journal of Hydrodynamics：Ser B，2013，25（4）：564-571.

［6］陳一帆，程海洋，萬(wàn)曉麗，等.結(jié)合糙率校正的河網(wǎng)水情數(shù)據(jù)同化［J］.水科學(xué)進(jìn)展，2015，26（5）：731-738.（CHEN Y F，CHENG H Y，WAN X L，et al.Research on data-assimilation combined with roughness correction for dynamic river system［J］.Advances in Water Science，2015，26（5）：731-738.（in Chinese））

［7］REN T，LIU X F，NIU J W，et al.Real-time water level prediction of cascaded channels based on multilayer perception and recurrent neural network［J］.Journal of Hydrology，2020，585：124783.

［8］徐興亞，方紅衛(wèi)，張?jiān)婪?，?河道洪水實(shí)時(shí)概率預(yù)報(bào)模型與應(yīng)用［J］.水科學(xué)進(jìn)展，2015，26（3）：356-364.（XU X Y，F(xiàn)ANG H W，ZHANG Y F，et al.A real-time probabilistic channel flood forecasting model and application based on particle filters［J］.Advances in Water Science，2015，26（3）：356-364.（in Chinese））

［9］王麟，趙濤.基于BP神經(jīng)網(wǎng)絡(luò)的人工加糙明渠糙率預(yù)測(cè)模型［J］.水電能源科學(xué)，2021，39（9）：125-127，124.（WANG L，ZHAO T.Prediction model of roughness coefficient of artificially roughened open channel based on BP neural network［J］.Water Resources and Power，2021，39（9）：125-127，124.（in Chinese））

［10］TIMBADIYA P V，PATEL P L，POREY P D.Calibration of HEC-RAS model on prediction of flood for Lower Tapi River，India［J］.Journal of Water Resource and Protection，2011，3（11）：805-811.

［11］陳一帆，程偉平，蔣建群.一種穩(wěn)健的河流糙率反演方法［J］.浙江大學(xué)學(xué)報(bào)（工學(xué)版），2013，47（8）：1361-1365，1378.（CHEN Y F，CHENG W P，JIANG J Q.A robust river roughness calibration model［J］.Journal of Zhejiang University （Engineering Science），2013，47（8）：1361-1365，1378.（in Chinese））

［12］高學(xué)平，呂建璋，孫博聞，等.含植物河道等效床面阻力試驗(yàn)研究［J］.水利學(xué)報(bào)，2021，52（9）：1024-1035，1046.（GAO X P，LYU J Z，SUN B W，et al.Experimental study on equivalent bed resistance of river containing vegetation［J］.Journal of Hydraulic Engineering，2021，52（9）：1024-1035，1046.（in Chinese））

［13］ATTARI M，HOSSEINI S M.A simple innovative method for calibration of Manning′s roughness coefficient in rivers using a similarity concept［J］.Journal of Hydrology，2019，575：810-823.

［14］ATTARI M，TAHERIAN M，HOSSEINI S M，et al.A simple and robust method for identifying the distribution functions of Manning′s roughness coefficient along a natural river［J］.Journal of Hydrology，2021，595：125680.

［15］ZHOU X，XU W，XIN K，et al.Self-adaptive calibration of real-time semand and roughness of water distribution systems［J］.Water Resources Research，2018，54（8）：5536-5550.

［16］LIU W，GUAN G H，TIAN X，et al.Exploiting a real-time self-correcting digital twin model for the middle route of the south-to-north water diversion project of China［J］.Journal of Water Resources Planning and Management，2023，149（7）：04023023.

［17］EVENSEN G.The ensemble Kalman filter：theoretical formulation and practical implementation［J］.Ocean Dynamics，2003，53（4）：343-367.

［18］LITRICO X，F(xiàn)ROMION V.Analytical approximation of open-channel flow for controller design［J］.Applied Mathematical Modelling，2004，28（7）：677-695.

［19］何建京，王惠民.流動(dòng)型態(tài)對(duì)曼寧糙率系數(shù)的影響研究［J］.水文，2002，22（6）：22-24，53.（HE J J，WANG H M.The effects of types of flow on manning′s roughness coefficient［J］.Hydrology，2002，22（6）：22-24，53.（in Chinese））

［20］王開，魏加華，王光謙.大型渠道糙率系數(shù)設(shè)計(jì)取值的不確定性及影響分析［J］.應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報(bào)，2008，16（6）：870-878.（WANG K，WEI J H，WANG G Q.Uncertainty in design value of roughness coefficient for large scale channel and effects analysis［J］.Journal of Basic Science and Engineering，2008，16（6）：870-878.（in Chinese））

［21］NAJAFZADEH M，OLIVETO G.Riprap incipient motion for overtopping flows with machine learning models［J］.Journal of Hydroinformatics，2020，22（4）：749-767.

［22］武漢大學(xué).輸水渠道系統(tǒng)運(yùn)行仿真與控制軟件：中國(guó)，2011SR034392［P/OL］.2011-06-03［2023-04-08］.（Wuhan University.Simulation and control of canal system：China，2011SR034392［P/OL］.2011-06-03［2023-04-08］.（in Chinese））

［23］ERRICO A，PASQUINO V，MAXWALD M，et al.The effect of flexible vegetation on flow in drainage channels：estimation of roughness coefficients at the real scale［J］.Ecological Engineering，2018，120：411-421.

Roughness estimation methods of hydrodynamic digital

twin models for canal systems

The study is financially supported by the National Natural Science Foundation of China （No.51979202;No.51879199）.

GUAN Guanghua1，LIU-WANG jiayi1，CHEN Xiaonan2，SHI Liangsheng1

（1. State Key Laboratory of Water Resources Engineering and Management，Wuhan University，Wuhan 430072，China;

2. China South-to-North Water Diversion Middle Route Corporation Limited，Beijing 100038，China）

Abstract：A segmented estimation method of roughness based on the variation of hydraulic radius and estimation accuracy is proposed in order to realize the real-time and high-fidelity roughness estimation of hydrodynamic digital twin （DT） models.This method considers the spatial variability of the roughness value in the longitudinal direction of canals.Based on canal segmentation，two different estimation frameworks，the independent estimation method and joint estimation method，are proposed.The ensemble Kalman filter algorithm is applied to estimate the roughness of each canal segment online based on the limited observed water levels.The results show that the two estimation methods can improve the accuracy of the model by 20%—50%.In addition，the independent estimation method is suitable for a complex canal system with a small error accumulation，while the joint estimation method is suitable for simple canals with unavailable observations.The proposed method can be used for parameter estimation and variable updating of hydrodynamic DT models，and provide a reference for the construction of DT water networks.

Key words：roughness;water transfer systems;digital twin;South-to-North Water Diversion Project;ensemble Kalman filter;segmented estimation