知識數(shù)據(jù)雙驅(qū)動的感潮河網(wǎng)水動力智能模擬方法

關(guān)鍵詞：水動力模擬；感潮河網(wǎng)；智能模擬；知識驅(qū)動；數(shù)據(jù)驅(qū)動

感潮河網(wǎng)地區(qū)地勢平緩，水動力弱，且受到下游潮汐頂托，水體常出現(xiàn)滯流和往復(fù)流現(xiàn)象，自凈能力差，河網(wǎng)水質(zhì)易惡化[1]。通過大量水閘、泵站工程的聯(lián)合調(diào)度提升水動力是提高河網(wǎng)水體自凈能力和改善水質(zhì)的重要手段，而感潮河網(wǎng)水動力過程的精準(zhǔn)快速模擬是重要前提。目前感潮河網(wǎng)的水動力模擬以水流運動方程計算和神經(jīng)網(wǎng)絡(luò)方法為主，如上海河網(wǎng)[2-3]。前者基于描述淺水水流運動的圣維南方程組，通過對方程進(jìn)行離散求解。這種方法通常需要將連續(xù)的水動力過程離散到有限數(shù)量的網(wǎng)格點，為了保證計算精度，需要劃分比較精細(xì)的網(wǎng)格，對于較大河網(wǎng)系統(tǒng)，計算效率低。近年來，神經(jīng)網(wǎng)絡(luò)等機器學(xué)習(xí)方法因其逼近非線性系統(tǒng)能力強、計算效率高等特點而越來越多地應(yīng)用于水動力模擬和水文預(yù)報。這類方法使用數(shù)據(jù)驅(qū)動的人工神經(jīng)網(wǎng)絡(luò)代替圣維南方程組，避開了方程求解過程，通過大量水位、流量等數(shù)據(jù)進(jìn)行訓(xùn)練學(xué)習(xí)，實現(xiàn)快速地模擬和預(yù)測，且通過更新神經(jīng)網(wǎng)絡(luò)算法可以進(jìn)一步提升模型的計算速度[4-6]。然而，神經(jīng)網(wǎng)絡(luò)模型往往對訓(xùn)練數(shù)據(jù)的需求大、精準(zhǔn)度要求高[7-8]，而且缺乏物理可解釋性。Raissi等[9]在2019年提出了一種集成了數(shù)據(jù)約束和物理約束的計算方法，可以彌補上述2種方法的不足。這種方法所需的訓(xùn)練數(shù)據(jù)量少，緩解了數(shù)據(jù)稀缺對機器學(xué)習(xí)模型的限制，而且因為加入物理約束而增加了可解釋性，已被成功應(yīng)用于流體力學(xué)、彈性力學(xué)、傳熱系統(tǒng)和量子化學(xué)等領(lǐng)域[10-13]。

本文提出一種知識數(shù)據(jù)雙驅(qū)動的感潮河網(wǎng)水動力智能模擬方法，并將其應(yīng)用于上海蕰南片感潮河網(wǎng)，實現(xiàn)感潮河網(wǎng)水動力的快速、精準(zhǔn)模擬，對感潮河網(wǎng)洪水快速預(yù)報、水動力時空格局優(yōu)化提供有力支持。

1模型原理與方法

1.1模型框架

傳統(tǒng)的人工神經(jīng)網(wǎng)絡(luò)模型將輸入層的數(shù)據(jù)經(jīng)過多個隱藏層計算后得到輸出層的預(yù)測值，這些預(yù)測值與實際值的差異為損失函數(shù)，該函數(shù)指導(dǎo)優(yōu)化算法調(diào)整網(wǎng)絡(luò)權(quán)重，從而完成人工神經(jīng)網(wǎng)絡(luò)的訓(xùn)練過程[14]。但是，對于沒有提供訓(xùn)練數(shù)據(jù)的位置，神經(jīng)網(wǎng)絡(luò)便無法依此過程完成訓(xùn)練，而且訓(xùn)練結(jié)果是否滿足客觀物理規(guī)律不得而知。若將物理方程嵌入損失函數(shù)中，計算每個輸入樣本對應(yīng)輸出值的方程殘差，以殘差作為損失函數(shù)的一部分，指導(dǎo)神經(jīng)網(wǎng)絡(luò)訓(xùn)練，既可以指導(dǎo)無訓(xùn)練數(shù)據(jù)位置的神經(jīng)網(wǎng)絡(luò)訓(xùn)練，訓(xùn)練得到的結(jié)果也滿足給定的物理方程，具有較好的可解釋性。據(jù)此，圖1給出了知識數(shù)據(jù)雙驅(qū)動（簡稱KDM）的感潮河網(wǎng)智能模擬方法框架，主干部分為圖1（a）部分所示的人工神經(jīng)網(wǎng)絡(luò)，在計算損失函數(shù)時嵌入如圖1（b）部分所示的物理方程，對人工神經(jīng)網(wǎng)絡(luò)訓(xùn)練結(jié)果進(jìn)行約束，其中，ε為人工神經(jīng)網(wǎng)絡(luò)的閾值，當(dāng)損失函數(shù)小于此閾值時表示訓(xùn)練完成，本文取值為0.0001；Jf為阻力項。

該方法中，人工神經(jīng)網(wǎng)絡(luò)作為主體框架，負(fù)責(zé)處理和預(yù)測數(shù)據(jù)。神經(jīng)網(wǎng)絡(luò)主要由輸入層、隱藏層和輸出層構(gòu)成?？紤]到預(yù)測水位和流量預(yù)測中所需捕捉的特征可能不同，這里分別采用2個相同尺寸的深度神經(jīng)網(wǎng)絡(luò)來實現(xiàn)對水位（h）和流量（Q）的預(yù)測，神經(jīng)網(wǎng)絡(luò)的輸入層X為

式中：f（）為用于提供非線性關(guān)系學(xué)習(xí)能力的激活函數(shù)；Whi和Bhi分別為第i層隱藏層的權(quán)重和偏置值；Zh（i-1）為第層隱藏層的輸出值。最后人工神經(jīng)網(wǎng)絡(luò)預(yù)測的水位可由最后一層隱藏層l的權(quán)重（Whl）、偏置（Bhl）和輸出值（Zhl）依照式（2），以相同方法計算。

1.2知識數(shù)據(jù)雙約束損失函數(shù)

數(shù)據(jù)約束部分的損失函數(shù)（L1（ ））由神經(jīng)網(wǎng)絡(luò)預(yù)測值與給定初始條件、邊界條件和水文測站觀測值的均方根差值構(gòu)成：

式中：為包含水位預(yù)測神經(jīng)網(wǎng)絡(luò)和流量預(yù)測神經(jīng)網(wǎng)絡(luò)全部權(quán)重和偏置的權(quán)重集；、和分別為初始條件、邊界條件和水文測站觀測值對應(yīng)的誤差；ic、bc和ob分別為初始條件、邊界條件和水文測站觀測值對應(yīng)誤差的權(quán)系數(shù)。、和可分別由式（4）—式（6）計算。

考慮到真實感潮河網(wǎng)中絕大部分時空位置的實測數(shù)據(jù)難以獲取，無法依照此思路實現(xiàn)神經(jīng)網(wǎng)絡(luò)的訓(xùn)練過程，且訓(xùn)練得到的結(jié)果是否滿足河網(wǎng)物理規(guī)律也不得而知。因此，引入一維圣維南方程組作為河道水流控制方程；引入節(jié)點能量平衡方程作為河網(wǎng)分匯流節(jié)點控制方程，計算各樣本輸出值對應(yīng)的方程殘差，作為損失函數(shù)的一部分，指導(dǎo)神經(jīng)網(wǎng)絡(luò)完成訓(xùn)練過程，同時也為預(yù)測結(jié)果提供了物理依據(jù)。

在河網(wǎng)模擬計算中，除了需要對單條河道進(jìn)行水力計算外，還需要分匯流節(jié)點等內(nèi)邊界條件處的控制方程。在分匯流節(jié)點處，圣維南方程組不再適用，由以下質(zhì)量守恒方程和節(jié)點水位方程控制[15]{Σ：

式中：ΣQi為河網(wǎng)節(jié)點上游各河段的入流流量之和；ΣQj為河網(wǎng)節(jié)點下游各河段的出流流量之和；hi為與河網(wǎng)節(jié)點處相連的第i條河段的水位，i=1，2，…，n。水流質(zhì)量守恒方程假設(shè)入流流量總和等于出流流量總和；節(jié)點水位方程假設(shè)河網(wǎng)節(jié)點處不存在水位突變，水位相等。

式中：Npde為方程殘差所用于計算的時空點數(shù)；rj分別為式（7）和式（8）中所示4個方程的殘差；Ωpde為方程殘差計算中的時空域。此時，本方法的損失函數(shù)（L（））由人工神經(jīng)網(wǎng)絡(luò)部分的損失函數(shù)（L1（））和方程對應(yīng)損失函數(shù)（（））組合而成：

1.3權(quán)系數(shù)自適應(yīng)調(diào)整

權(quán)系數(shù)！ic、！bc和！ob可以用來平衡不同損失項的數(shù)量級，避免某些損失項因數(shù)量級較大而在訓(xùn)練過程中主導(dǎo)了其他損失項，從而影響優(yōu)化效果。而權(quán)系數(shù)的優(yōu)選是隨著不同求解問題和神經(jīng)網(wǎng)絡(luò)優(yōu)化過程不斷變化的，人為調(diào)整權(quán)系數(shù)費時且難以獲取權(quán)系數(shù)的最佳取值。本文中采用動態(tài)權(quán)系數(shù)法解決此問題，該方法使權(quán)系數(shù)在神經(jīng)網(wǎng)絡(luò)權(quán)重更新過程中動態(tài)調(diào)整以平衡損失函數(shù)的梯度，神經(jīng)網(wǎng)絡(luò)的權(quán)重更新過程可表示為

由于輸入數(shù)據(jù)之間的尺度差異較大，在誤差反向傳播時可能會出現(xiàn)梯度爆炸的現(xiàn)象，本文采用以下量綱一化方法將輸入數(shù)據(jù)縮放到[?1，1]區(qū)間內(nèi)：

式中：Xre為量綱一化的輸入層。

1.4參數(shù)設(shè)置

神經(jīng)網(wǎng)絡(luò)尺寸的選擇對模擬結(jié)果有著至關(guān)重要的影響。若選擇的神經(jīng)網(wǎng)絡(luò)尺寸過小，將難以捕捉到輸入和輸出之間的復(fù)雜關(guān)系，導(dǎo)致神經(jīng)網(wǎng)絡(luò)欠擬合；若神經(jīng)網(wǎng)絡(luò)尺寸過大，會過多地捕捉訓(xùn)練數(shù)據(jù)的細(xì)節(jié)，導(dǎo)致過擬合問題，同時增加計算的負(fù)擔(dān)。在預(yù)試驗中，本文測試了增加隱藏層和增加神經(jīng)元的案例，結(jié)果顯示增加的隱藏層和神經(jīng)元會顯著提升訓(xùn)練時間，而對模擬精度幾乎沒有提升，甚至過擬合導(dǎo)致精度降低。為此，本文在模型驗證中使用2個隱藏層，每層30個神經(jīng)元的全連接結(jié)構(gòu)構(gòu)造本方法中的神經(jīng)網(wǎng)絡(luò)部分，加入系數(shù)為0.01的正則化方法防止過擬合現(xiàn)象，提高神經(jīng)網(wǎng)絡(luò)的魯棒性。神經(jīng)網(wǎng)絡(luò)的激活函數(shù)f使用tanh激活函數(shù)，神經(jīng)網(wǎng)絡(luò)權(quán)重的初始化方法選擇適用于對稱激活函數(shù)tanh的xavier初始化，用于尋找最小損失的優(yōu)化器使用Adam，神經(jīng)網(wǎng)絡(luò)的學(xué)習(xí)率為0.005，訓(xùn)練迭代次數(shù)為30000次，權(quán)系數(shù)ic、bc和ob初始取值均為10。

2模型驗證

2.1案例設(shè)置

本文設(shè)置了概化河網(wǎng)案例，用以探討本方法模擬感潮河網(wǎng)水動力過程的可行性（圖2）。河網(wǎng)由8條河道和4個交匯分汊節(jié)點組成，擁有2個入口和2個出口邊界。入口邊界給定非恒定流量過程，出口邊界給定潮位，河道斷面均為矩形，河長為1000m，糙率為0.028m?1/3/s，河道其他物理信息如表1所示?？紤]到現(xiàn)實河流中存在水文測站可以提供當(dāng)?shù)氐乃弧⒘髁康葧r間序列數(shù)據(jù)，這里也提供若干斷面（即水文測站）的時間序列水位和流量數(shù)據(jù)以提高模擬精度，這些數(shù)據(jù)由一維水動力模擬軟件HEC-RAS計算值提供。最后，將本方法預(yù)測結(jié)果與HEC-RAS軟件計算結(jié)果對比，分析本方法的可行性。

模擬案例參數(shù)設(shè)置如下：使用HEC-RAS求解時，入口邊界為1-8和1-7河道上游，分別設(shè)置為50～130m3/s和20～50m3/s的流量過程（圖3）；出口邊界為1-1和1-2河道下游，設(shè)置為潮位邊界，水位變動幅度為20cm（圖3）；計算時長為24h，時間步長為15min，河道每隔100m設(shè)置1個斷面，共88個斷面。使用KDM求解時，訓(xùn)練的初始條件為t=0s時全部河道各斷面的水位和流量，邊界條件為1-8和1-7河道x=0m以及1-2和1-1河道x=1000m處的水位、流量時間序列數(shù)據(jù)，水文測站設(shè)在8條河道x=500m處共8個，計算中的Nic、Nbc、Nob和Npde分別為88、194、776和8536。由于案例中提供的數(shù)據(jù)點比較稀疏，為防止過擬合，設(shè)置正則化方法系數(shù)為0.1。

本文采用L2相對誤差衡量KDM模擬結(jié)果與HEC-RAS參考解之間的差異，評估KDM的模擬效果。通過計算兩者之間的歐氏距離，并進(jìn)行量綱一化處理，以便比較不同規(guī)模的數(shù)據(jù)。對水位和流量的L2相對誤差可由下式計算：

2.2案例模擬結(jié)果分析

圖4為概化河網(wǎng)模擬中各河道水位KDM解和HEC-RAS參考解的結(jié)果對比圖。顯然，KDM方法顯現(xiàn)出很高的精確性，水位和流量的L2相對誤差分別為3.83%和7.49%，潮汐邊界對水位的影響主要集中在河道1-7和1-8下游500m范圍內(nèi)，這種影響能被KDM很好地捕捉到。其中，水位的最大誤差在0.2m左右，這些位置集中在分匯流節(jié)點附近，誤差成因主要是交匯口的壅水和水位突變。分匯流節(jié)點處時常伴隨著復(fù)雜的水流結(jié)構(gòu)[16-19]，節(jié)點水位方程和質(zhì)量守恒無法精確捕捉分匯流節(jié)點處的水動力過程。若以更加精確的方程實現(xiàn)對河網(wǎng)分匯流節(jié)點的控制，本方法的預(yù)測精度將進(jìn)一步提升。同時，節(jié)點附近缺乏測量數(shù)據(jù)也是導(dǎo)致誤差的原因之一，可以通過增設(shè)交匯/分汊口附近的測站以提高水動力的預(yù)測精度。在模擬過程中，權(quán)系數(shù)ic、bc和ob由初始值10上升至10000附近，最后穩(wěn)定在9000左右。這表明，訓(xùn)練初期L1（）的極速收斂導(dǎo)致中期下降速度遠(yuǎn)小于，并且一直保持這樣的趨勢直至訓(xùn)練完成。除了計算精度，本方法的計算效率很高，本案例中HEC-RAS的計算時間為47s，而訓(xùn)練完成之后的KDM為1s。

2.3魯棒性分析

感潮河網(wǎng)實際測量數(shù)據(jù)往往會伴隨有一定的噪聲。為檢驗本方法的抗噪聲能力，在參考解中加入不同水平的噪聲（0.5%、1%、2%、5%、10%）用于KDM訓(xùn)練，對應(yīng)誤差如表2所示。隨著噪聲水平從0.5%增加到2%，KDM解仍然具有優(yōu)異的表現(xiàn)，水位誤差控制在6%以內(nèi)，流量誤差略微升高到19%左右；當(dāng)噪聲水平從5%增加到10%時，水位誤差從6%左右增加到10%左右，流量誤差升高到26%左右，模擬結(jié)果受影響較大。KDM在5%內(nèi)的噪聲影響下仍可以精確模擬感潮河網(wǎng)水動力過程，表明本方法具有良好的魯棒性。

3模型應(yīng)用

3.1模擬案例介紹

本節(jié)通過模擬真實感潮河網(wǎng)來進(jìn)一步驗證方法的可行性。上海蕰南片河網(wǎng)由楊樹浦港水系（包括虬江、東走馬塘和楊樹浦港）、虹口港水系（包括南泗塘、西泗塘、沙涇港、俞涇浦、虹口港）、彭越浦-木瀆港水系（包括東茭涇、彭越浦、桃浦河、木瀆港）和新江灣城水系組成，全片總面積為173.3km2，片內(nèi)河道總長159.51km，河道總面積為3.69km2，河網(wǎng)密度為0.92km/km2。該區(qū)域水系較復(fù)雜，河網(wǎng)水動力不足、水流往復(fù)現(xiàn)象突出[6]，通過河網(wǎng)水閘調(diào)度提升河網(wǎng)水動力是改善水環(huán)境狀況的主要措施。本案例模擬了蕰藻浜以南和蘇州河以北的整個河網(wǎng)共33條河道28個節(jié)點的水動力過程，模擬情景為2022年5月1日全天，選擇往復(fù)流現(xiàn)象突出的虹口港水系進(jìn)行水動力模擬結(jié)果分析。虹口港水系存在7條河道（河道編號為2-1—2-7）、3個分匯流節(jié)點以及3個水文測站，如圖5所示。

考慮到感潮河網(wǎng)的復(fù)雜性，本案例神經(jīng)網(wǎng)絡(luò)設(shè)置3層隱藏層，每層30個神經(jīng)元。模擬過程中，不提供初始條件，給定河網(wǎng)南側(cè)4個閘泵（虹口港水閘、彭越浦泵閘、真如港泵閘、木瀆港泵閘）以及北側(cè)5個閘泵（郝橋港泵閘、西泗塘泵閘、東茭涇泵閘、西彌浦水閘、桃浦河泵閘）的閘內(nèi)當(dāng)日水位、流量數(shù)據(jù)作為邊界條件。這里為了計算簡便，未考慮閘門調(diào)度，直接給定閘內(nèi)數(shù)據(jù)，數(shù)據(jù)的時間間隔為5min。河網(wǎng)內(nèi)存在江灣、涼城新村和虎林公園3個水文測站。以下分別討論虎林公園作為訓(xùn)練、江灣及涼城新村作為驗證，虎林公園及江灣作為訓(xùn)練、涼城新村作為驗證，以及3個測站同時作為訓(xùn)練的3種工況下3個水文測站真實測量數(shù)據(jù)與模擬結(jié)果的吻合程度，并將3個水文測站均用于訓(xùn)練，分析上海蕰南片河網(wǎng)虹口港水系的水動力模擬結(jié)果。模擬過程中，時間步長為15min，整個模擬河網(wǎng)選取187個斷面參與殘差計算。

3.2模擬結(jié)果驗證

以虎林公園水文測站作為訓(xùn)練數(shù)據(jù)，其余兩測站作為驗證，模擬水位、流量結(jié)果如圖6（a）、圖6（b）所示。由于缺乏初始水動力數(shù)據(jù)，且俞涇浦、沙涇港處均不存在訓(xùn)練數(shù)據(jù)，兩河的分流比以及水面坡降都難以確定，位于兩河上的涼城新村和江灣水文測站的水位模擬結(jié)果在關(guān)閘時刻基本與實測值吻合，但是在開閘時刻的漲落趨勢和幅度難以預(yù)測。以虎林公園、江灣水文測站作為訓(xùn)練數(shù)據(jù)，涼城新村作為驗證的模擬結(jié)果如圖6（c）、圖6（d）所示?？梢钥吹?，在給定了江灣的水位、流量數(shù)據(jù)之后，涼城新村的水位時間序列數(shù)據(jù)已經(jīng)基本與實測數(shù)據(jù)吻合；流量模擬結(jié)果也得到了很大的改善，但仍存在部分誤差。將3個測站的數(shù)據(jù)均作為訓(xùn)練數(shù)據(jù)的模擬結(jié)果如圖6（e）、圖6（f）所示。從3個測站的實測水位流量可以看出，主要的開閘時間在0～4h和12～18h，在這2個時間段中，3個測站的水位受到潮汐影響而產(chǎn)生波動，流量也隨之發(fā)生劇烈震蕩，在這種影響下，3個測站的水位、流量時間序列數(shù)據(jù)也都得到了很好的模擬，且變化趨勢也更加平滑，部分波動數(shù)據(jù)被模型認(rèn)為是測量誤差而沒有進(jìn)行學(xué)習(xí)，這表示潮汐變化過程能夠被捕捉到，也再次表明本方法具有良好的魯棒性。

3.3虹口港水動力模擬結(jié)果分析

虹口港水系各河道的水位、流量模擬結(jié)果如圖7所示。當(dāng)日0～4h及12～16h左右為虹口港開閘時間，虹口港閘內(nèi)水位提升了0.5m左右，流量提升了30m3/s左右，河道2-1的水動力也因此得到了明顯的提升；河道2-2與2-3位于河道2-1下游，且擁有相似的分流比，兩河道水動力相當(dāng)，都能很好地排出河道積水；河道2-5和河道2-6分別與河道2-2和河道2-3下游相連，且4條河道流向基本一致，具有大致相同的水動力特征；河道2-4和河道2-7作為東西走向的河道，幾乎與其余河道正交，且南北側(cè)閘門啟閉時間接近，因此，水動力較強的南北走向河道中的水流難以流入這2條河道，導(dǎo)致其水動力較差，整體流量較小，水流往復(fù)現(xiàn)象嚴(yán)重。這種往復(fù)水流往往不利于水體自凈和活水提質(zhì)，是河網(wǎng)局部水質(zhì)變差的重要原因[3]。

4結(jié)論

本文提出了知識數(shù)據(jù)雙驅(qū)動的感潮河網(wǎng)水動力智能模擬方法，應(yīng)用于感潮河網(wǎng)案例的水動力模擬，并與河網(wǎng)水動力模擬軟件的計算結(jié)果進(jìn)行驗證分析，同時應(yīng)用于上海蕰南片感潮河網(wǎng)水動力模擬。主要結(jié)論如下：

（1）本方法以河網(wǎng)水流控制方程約束人工神經(jīng)網(wǎng)絡(luò)，在保證計算效率的基礎(chǔ)上，提高了物理可解釋性，而且需要的訓(xùn)練數(shù)據(jù)大幅度減少，還可以得到?jīng)]有訓(xùn)練數(shù)據(jù)斷面的水位和流量。

（2）在感潮河網(wǎng)案例模擬中，該方法能夠捕捉潮汐對水位的影響，且具有較高的模擬精度和良好的魯棒性。由于分匯節(jié)點處水流復(fù)雜、水位變化大，水動力模擬較現(xiàn)實河網(wǎng)存在一定偏差。為此，在控制方程的選取方面仍然有待優(yōu)化，可以將節(jié)點水位方程和質(zhì)量守恒方程替換為更加準(zhǔn)確的節(jié)點控制方程，實現(xiàn)對分匯流節(jié)點處水動力過程的精準(zhǔn)捕捉以增加模擬精度。此外，針對水位梯度在分匯節(jié)點處的突變特性，可以通過在這些關(guān)鍵位置增設(shè)水文測站來有效降低預(yù)測偏差。上海蕰南片感潮河網(wǎng)水動力模擬結(jié)果顯示，虹口港水系中東西走向的河道水動力較弱，難以排出往復(fù)積水，導(dǎo)致河網(wǎng)水質(zhì)變差?？梢酝ㄟ^控制南北閘門啟閉及分匯流節(jié)點處的工程措施控制分流比，提升東西走向河道的水動力。

（3）本方法未來可以將水閘、泵站的調(diào)度規(guī)則納入模型中，進(jìn)一步耦合河網(wǎng)水環(huán)境、水生態(tài)相關(guān)控制方程，構(gòu)建感潮河網(wǎng)水動力-水質(zhì)-水生態(tài)耦合智能模型，服務(wù)于感潮河網(wǎng)地區(qū)水問題綜合治理，為流域數(shù)字孿生、國家水網(wǎng)建設(shè)提供支撐。