本發(fā)明涉及水動力模擬和水文預(yù)報(bào)技術(shù)領(lǐng)域，特別涉及一種基于水動力數(shù)值模擬的感潮河網(wǎng)潮位預(yù)報(bào)方法。

背景技術(shù)：

隨著河口及三角洲經(jīng)濟(jì)建設(shè)的不斷發(fā)展，對水文預(yù)報(bào)工作提出了更高要求，如何快速準(zhǔn)確地進(jìn)行潮位預(yù)報(bào)，成為當(dāng)前急需解決的問題。防汛決策、水利樞紐設(shè)計(jì)參考、閘泵群啟閉方案優(yōu)化等對感潮網(wǎng)河區(qū)的潮位預(yù)報(bào)需求日益增加。感潮河網(wǎng)的水位變化過程受上游徑流、下游潮汐、河汊分流、區(qū)間入流及降水等多種因素共同影響，其水文特性既有內(nèi)陸河流的特性，表現(xiàn)出洪、中、枯水期的季節(jié)變化，又有海洋的某些特征，表現(xiàn)出大、中、小潮的變化和潮汐引起水位漲落與流速流向的變化，此外，還受陸地氣候的影響，上游徑流不僅年內(nèi)有明顯的變化，而且在年際間變化也很不相同；受河道地形的影響也較大，情況十分復(fù)雜，感潮河網(wǎng)內(nèi)潮位預(yù)報(bào)是長期以來河口及三角洲治理開發(fā)中的一個(gè)難點(diǎn)問題。

傳統(tǒng)的潮位預(yù)報(bào)多為高水位預(yù)報(bào)，一般采用相應(yīng)水位法，即由某站水位、傳輸距離推求另一站水位，其預(yù)報(bào)有效期受傳輸時(shí)間左右，無法滿足長期、中期預(yù)報(bào)需求；另有分離過程法將預(yù)報(bào)站水位分離成上、下游水位相關(guān)的兩部分分別進(jìn)行預(yù)報(bào)，其不足之處與相關(guān)法相同，主要是預(yù)報(bào)時(shí)間較為受限；神經(jīng)網(wǎng)絡(luò)法預(yù)報(bào)水位則主要存在可能由于局部極值難收斂、振幅不大時(shí)較精準(zhǔn)、高水位預(yù)報(bào)誤差大的問題。為滿足三角洲對潮位預(yù)報(bào)需求，發(fā)展一種基于數(shù)值模擬的能實(shí)現(xiàn)感潮河網(wǎng)內(nèi)任意水域潮位快速預(yù)報(bào)的方法是十分必要的。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明的目的在于克服河網(wǎng)地區(qū)現(xiàn)有潮位預(yù)報(bào)方法誤差較大、適用范圍小的不足，提供一種基于水動力數(shù)學(xué)模型的感潮河網(wǎng)任意水域潮位快速預(yù)報(bào)及實(shí)時(shí)校正方法，提高感潮河網(wǎng)潮位預(yù)報(bào)的適用范圍和精度。

本發(fā)明的目的通過下述技術(shù)方案實(shí)現(xiàn)：

一種基于水動力數(shù)值模擬的感潮河網(wǎng)潮位預(yù)報(bào)方法，包括以下步驟：

步驟S1、根據(jù)預(yù)報(bào)位置選擇合適的研究范圍，確定河口下邊界及上游控制站點(diǎn)與邊界；

步驟S2、對研究范圍內(nèi)河道進(jìn)行斷面剖分和地形插值；

步驟S3、建立上游控制站至下游河口邊界之間的網(wǎng)河區(qū)的一維非恒定流水動力數(shù)學(xué)模型；

步驟S4、采用實(shí)測水文數(shù)據(jù)對一維非恒定流水動力數(shù)學(xué)模型進(jìn)行率定和驗(yàn)證；

步驟S5、根據(jù)河口控制站點(diǎn)調(diào)和參數(shù)進(jìn)行潮位過程預(yù)報(bào)計(jì)算：

步驟S6、根據(jù)實(shí)測數(shù)據(jù)對河口潮位邊界進(jìn)行校正；

步驟S7、采用加權(quán)移動平均方法推求上游控制站點(diǎn)流量過程；

步驟S8、根據(jù)實(shí)時(shí)實(shí)測數(shù)據(jù)對上游流量邊界進(jìn)行校正；

步驟S9、輸入預(yù)報(bào)上游邊界流量和下游河口邊界潮位過程，進(jìn)行網(wǎng)河區(qū)的一維非恒定流水動力數(shù)學(xué)模型計(jì)算，求得網(wǎng)河區(qū)內(nèi)各斷面預(yù)潮位過程。

進(jìn)一步地，所述步驟S3、建立上游控制站至下游河口邊界之間的網(wǎng)河區(qū)的一維非恒定流水動力數(shù)學(xué)模型具體過程如下：

步驟S31、采用圣維南方程組作為一維潮流數(shù)學(xué)模型，方程如下：

連續(xù)方程：

動量方程：

將相關(guān)支流匯入或流出點(diǎn)作為網(wǎng)河區(qū)汊點(diǎn)，汊口點(diǎn)水流按以下連接條件控制：

流量連接條件：

水位連接條件：Z_i,j＝Z_m,n＝……＝Z_l,k；

式中，Z為斷面平均水位；Q，A，R分別為斷面流量、過水面積、水力 半徑；B為水面寬度；X,T為距離和時(shí)間；G為重力加速度；謝才系數(shù)C＝R^1/6/n，n為糙率；

步驟S32、采用四點(diǎn)加權(quán)Preissmann固定網(wǎng)格隱式差分格式將方程離散，采用一維河網(wǎng)三級聯(lián)解算法求解方程，具體求解方法為：進(jìn)行網(wǎng)格劃分，以S代表流量Q和水位Z，則S在Δx河段、Δt時(shí)段內(nèi)的加權(quán)平均量及相應(yīng)偏導(dǎo)數(shù)可分別表示為：

式中：θ為加權(quán)系數(shù)，θ一般取0.5～1.0。

進(jìn)一步地，所述步驟S5、根據(jù)河口控制站點(diǎn)調(diào)和參數(shù)進(jìn)行潮位過程預(yù)報(bào)計(jì)算具體過程如下：

步驟S51、有分潮參數(shù)的直接根據(jù)預(yù)報(bào)時(shí)間計(jì)算潮位過程；

步驟S52、無分潮參數(shù)，但有長于一年的實(shí)測歷史數(shù)據(jù)的，調(diào)和分析各分潮參數(shù)，然后重復(fù)步驟S51；

步驟S53、分潮參數(shù)和歷史數(shù)據(jù)均無則參考臨近站點(diǎn)簡化處理。

進(jìn)一步地，所述步驟S51、有分潮參數(shù)的直接根據(jù)預(yù)報(bào)時(shí)間計(jì)算潮位過程的具體計(jì)算公式如下：

式中：為t時(shí)刻的潮位預(yù)報(bào)值，m；A₀為從某基準(zhǔn)面算起的平均海平面高度，m；f_i為第i個(gè)分潮的交點(diǎn)因子；H_i為第i個(gè)分潮的振幅，m；σ_i為分潮角速度；(v0+u)_i為分潮的天文初相角；g_i為分潮遲角；i為分潮編號；N為分潮總數(shù)；H_i、g_i亦稱為調(diào)和常數(shù)。

進(jìn)一步地，所述步驟S6、根據(jù)實(shí)測數(shù)據(jù)對河口潮位邊界進(jìn)行修正具體過程如下：

步驟S61、在[t0-26,t0]范圍內(nèi)搜索判別出高高潮位Z_H^PRE、低低潮位Z_L^PRE，計(jì)算潮幅Z_AM^PRE，Z_AM^PRE＝Z_H^PRE+Z_L^PRE；

步驟S62、同樣在[t0-26,t0]范圍內(nèi)搜索出實(shí)測潮位過程中的高高潮位Z_H^MESU、低低潮位Z_L^MESU，計(jì)算潮幅Z_AM^MESU，Z_AM^MESU＝Z_H^MESU+Z_L^MESU；

步驟S63、根據(jù)步驟S61及步驟S62結(jié)果對預(yù)報(bào)潮位中軸線進(jìn)行調(diào)整，使得Z_AM^PRE＝Z_AM^MESU；

步驟S64、計(jì)算中軸線偏移量DIF_Z_ZERO，DIF_Z_ZERO＝(Z_AM^PRE-Z_AM^MESU)/2.0；

步驟S65、對于預(yù)報(bào)潮位過程值Z^PRE逐一進(jìn)行中軸線偏移校正，

Z^PRE＝Z^PRE-DIF_Z_ZERO；

步驟S66、計(jì)算潮幅偏差DIF_Z_AM，DIF_Z_AM＝Z_AM^PRE/Z_AM^MESU；

進(jìn)一步地，所述步驟S7中加權(quán)移動平均方法的公式如下：

$Q_{pree}^t=Q_A^t*2+Q_{pree}^{t-1}$

式中：為t時(shí)刻預(yù)報(bào)流量，m³/s；為第t時(shí)刻的流量趨勢率的修正量，m³/s；為第t-1時(shí)刻的預(yù)報(bào)流量，m³/s。

進(jìn)一步地，所述流量趨勢率的修正量的求解步驟如下：

步驟S71、對于每一段預(yù)報(bào)期(t0～tn)內(nèi)的首、末趨勢率：式中：分別為預(yù)報(bào)期內(nèi)首時(shí)刻t0、末時(shí)刻tn對應(yīng)流量趨勢率；Q_T0-1和Q_T0-2分別為預(yù)報(bào)前1、前2小時(shí)的實(shí)測流量；

步驟S72、按線性插值法計(jì)算每一段預(yù)報(bào)期內(nèi)首末時(shí)刻以外其它各時(shí)刻趨勢率，

式中：為預(yù)報(bào)期內(nèi)首末時(shí)刻以外其它各時(shí)刻趨勢率；T為預(yù)報(bào)期總時(shí)長。

本發(fā)明相對于現(xiàn)有技術(shù)具有如下的優(yōu)點(diǎn)及效果：

本發(fā)明公開了一種基于水動力數(shù)學(xué)模型的感潮河網(wǎng)任意水域潮位快速預(yù)報(bào)及實(shí)時(shí)校正方法，綜合運(yùn)用了計(jì)算流體力學(xué)(CFD)方法和水位預(yù)報(bào)方法，解決了感潮河網(wǎng)內(nèi)水位受多汊道分流及外海潮汐影響而難以準(zhǔn)確預(yù)報(bào)的問題，提高感潮河網(wǎng)潮位預(yù)報(bào)的適用范圍和精度。

附圖說明

圖1是本發(fā)明公開的的一種基于水動力數(shù)值模擬的感潮河網(wǎng)潮位預(yù)報(bào)方法 的流程步驟圖；

圖2是一維河網(wǎng)水動力模型斷面剖分示意圖；

圖3(a)是本發(fā)明在中順大圍調(diào)度系統(tǒng)中一維河網(wǎng)水動力模型研究范圍圖；

圖3(b)是本發(fā)明在中順大圍調(diào)度系統(tǒng)中中順大圍區(qū)域的研究范圍圖；

圖4(a)是98.6洪水期三沙口站水位過程驗(yàn)證成果圖；

圖4(b)是98.6洪水期勒竹站水位過程驗(yàn)證成果圖；

圖4(c)是05.1枯水期小欖站水位過程驗(yàn)證成果圖；

圖4(d)是05.1枯水期三沙口站水位過程驗(yàn)證成果圖；

圖4(e)是05.1枯水期竹銀站水位過程驗(yàn)證成果圖；

圖4(f)是05.1枯水期亭角站水位過程驗(yàn)證成果圖；

圖5為河口控制點(diǎn)實(shí)測潮位、預(yù)報(bào)潮位、預(yù)報(bào)修正潮位對比圖；

圖6為上游梧州站實(shí)測流量、預(yù)報(bào)流量、預(yù)報(bào)修正流量對比圖；

圖7(a)東河水閘實(shí)測潮位與預(yù)報(bào)潮位過程對比圖；

圖7(b)西河水閘實(shí)測潮位與預(yù)報(bào)潮位過程對比圖。

具體實(shí)施方式

為使本發(fā)明的目的、技術(shù)方案及優(yōu)點(diǎn)更加清楚、明確，以下參照附圖并舉實(shí)施例對本發(fā)明進(jìn)一步詳細(xì)說明。應(yīng)當(dāng)理解，此處所描述的具體實(shí)施例僅僅用以解釋本發(fā)明，并不用于限定本發(fā)明。

實(shí)施例

請參見圖1，圖1是本實(shí)施例中公開的一種基于水動力數(shù)值模擬的感潮河網(wǎng)潮位預(yù)報(bào)方法的流程步驟圖。圖1所示的一種基于水動力數(shù)值模擬的感潮河網(wǎng)潮位預(yù)報(bào)方法，具體包括以下步驟：

步驟S1、根據(jù)預(yù)報(bào)位置選擇合適的研究范圍，確定河口下邊界及上游控制站點(diǎn)與邊界；

步驟S2、對研究范圍內(nèi)河道進(jìn)行斷面剖分和地形插值；

步驟S3、建立上游控制站至下游河口邊界之間的網(wǎng)河區(qū)的一維非恒定流水動力數(shù)學(xué)模型；具體步驟如下：

步驟S31、一維潮流數(shù)學(xué)模型采用圣維南方程組，方程如下：

連續(xù)方程：

動量方程：

相關(guān)支流匯入或流出點(diǎn)作為網(wǎng)河區(qū)汊點(diǎn)，汊口點(diǎn)水流按以下連接條件控制：

流量連接條件：

水位連接條件：Z_i,j＝Z_m,n＝……＝Z_l,k；

式中，Z為斷面平均水位；Q，A，R分別為斷面流量、過水面積、水力半徑；B為水面寬度；X,T為距離和時(shí)間；G為重力加速度；謝才系數(shù)C＝R^1/6/n，n為糙率。

步驟S32、方程離散采用四點(diǎn)加權(quán)Preissmann固定網(wǎng)格隱式差分格式，采用目前應(yīng)用廣泛的一維河網(wǎng)三級聯(lián)解算法，網(wǎng)格布置圖2所示。

具體求解方法為：對于圖2所示的網(wǎng)格劃分，以S代表流量Q和水位Z，則S在Δx河段、Δt時(shí)段內(nèi)的加權(quán)平均量及相應(yīng)偏導(dǎo)數(shù)可分別表示為：

式中：θ為加權(quán)系數(shù)，θ一般取0.5～1.0。

步驟S4、采用實(shí)測水文數(shù)據(jù)對一維非恒定流水動力數(shù)學(xué)模型進(jìn)行率定和驗(yàn)證；

步驟S5、根據(jù)河口控制站點(diǎn)調(diào)和參數(shù)進(jìn)行潮位過程預(yù)報(bào)計(jì)算，具體過程如下：

步驟S51、有分潮參數(shù)的直接根據(jù)預(yù)報(bào)時(shí)間計(jì)算潮位過程；

步驟S52、無分潮參數(shù)，但有長于一年的實(shí)測歷史數(shù)據(jù)的，調(diào)和分析各分潮參數(shù)，然后重復(fù)步驟S51；

步驟S53、分潮參數(shù)和歷史數(shù)據(jù)均無則參考臨近站點(diǎn)簡化處理。

式中：為t時(shí)刻的潮位預(yù)報(bào)值，m；A₀為從某基準(zhǔn)面算起的平均海平 面高度，m；f_i為第i個(gè)分潮的交點(diǎn)因子；H_i為第i個(gè)分潮的振幅，m；σ_i為分潮角速度；(v0+u)_i為分潮的天文初相角；g_i為分潮遲角；i為分潮編號；N為分潮總數(shù)；H_i、g_i亦稱為調(diào)和常數(shù)。

步驟S6、根據(jù)實(shí)測數(shù)據(jù)對河口潮位邊界進(jìn)行修正；

步驟S7、采用加權(quán)移動平均方法推求上游控制站點(diǎn)流量過程，并同時(shí)進(jìn)行修正；

式中：為t時(shí)刻預(yù)報(bào)流量，m³/s；為第t時(shí)刻的流量趨勢率(修正量)，m³/s；為第t-1時(shí)刻的預(yù)報(bào)流量，m³/s。

步驟S8、輸入預(yù)報(bào)上游邊界流量和下游河口邊界潮位過程，進(jìn)行網(wǎng)河區(qū)的一維非恒定流水動力數(shù)學(xué)模型計(jì)算，求得網(wǎng)河區(qū)內(nèi)各斷面預(yù)潮位過程。

其中，步驟S6所述潮位修正步驟如下：

步驟S61、在[t0-26,t0](預(yù)報(bào)起點(diǎn)前一個(gè)潮期，約26h)范圍內(nèi)搜索判別出高高潮位Z_H^PRE、低低潮位Z_L^PRE，計(jì)算潮幅Z_AM^PRE；

步驟S62、同樣在[t0-26,t0]范圍內(nèi)搜索出實(shí)測潮位過程中的高高潮位Z_H^MESU、低低潮位Z_L^MESU，計(jì)算潮幅Z_AM^MESU；

步驟S63、根據(jù)步驟S61及步驟S62結(jié)果對預(yù)報(bào)潮位中軸線進(jìn)行調(diào)整，使得(Z_H^PRE+Z_L^PRE)＝(Z_H^MESU+Z_L^MESU)；

步驟S64、計(jì)算中軸線偏移量DIF_Z_ZERO，

DIF_Z_ZERO＝((Z_H^PRE+Z_L^PRE)-(Z_H^MESU+Z_L^MESU))/2.0；

步驟S65、對于預(yù)報(bào)潮位過程值Z^PRE逐一進(jìn)行中軸線偏移校正，

Z^PRE＝Z^PRE-DIF_Z_ZERO；

步驟S66、計(jì)算潮幅偏差DIF_Z_AM，DIF_Z_AM＝Z_AM^PRE/Z_AM^MESU；

其中，步驟S7所述流量修正中的求解步驟如下：

步驟S71、對于每一段預(yù)報(bào)期(t0～tn)內(nèi)的首、末趨勢率：式中：分別為預(yù)報(bào)期內(nèi)首時(shí)刻t0、末時(shí)刻tn對應(yīng)流量趨勢率；Q_T0-1和Q_T0-2分別為預(yù)報(bào)前1、前2小時(shí)的實(shí)測流量。

步驟S72、按線性插值法計(jì)算每一段預(yù)報(bào)期內(nèi)首末時(shí)刻以外其它各時(shí)刻趨勢率，

式中：為預(yù)報(bào)期內(nèi)首末時(shí)刻以外其它各時(shí)刻趨勢率；T為預(yù)報(bào)期總時(shí)長。

實(shí)施例二

以本發(fā)明在珠江三角洲河網(wǎng)中順大圍調(diào)度系統(tǒng)外江潮位預(yù)報(bào)實(shí)例及附圖對本發(fā)明作進(jìn)一步詳細(xì)的描述。

(1)根據(jù)中順大圍的位置特性，確定潮位預(yù)報(bào)模型研究范圍——包了西、北江三角洲、東江三角洲及廣州出海水道、潭江水道；上邊界取以下控制水文(位)站：西江梧州、北江石角、老鴉崗(流溪河)、麒麟咀(增江)、東江博羅、潭江石咀站；下邊界取至八大口門控制站點(diǎn)：大虎、南沙、馮馬廟、橫門和磨刀門燈籠山、黃金、西炮臺、官沖，見附圖3(a)和圖3(b)；

(2)對河道進(jìn)行斷面剖分和地形插值，共設(shè)4140個(gè)斷面，采用2005年實(shí)測地形；

(3)建立包括珠江三角洲網(wǎng)河區(qū)的一維非恒定流水動力數(shù)學(xué)模型；

(4)采用1998年6月、1999年7月、2005年1月實(shí)測水文資料對珠江三角洲網(wǎng)河區(qū)一維非恒定流水動力數(shù)學(xué)模型進(jìn)行率定和驗(yàn)證，部分成果見圖4(a)～圖4(f)；

(5)根據(jù)河口控制站點(diǎn)調(diào)和參數(shù)進(jìn)行潮位過程預(yù)報(bào)計(jì)算，8個(gè)下邊界控制站點(diǎn)無分潮參數(shù)資料，其分潮參數(shù)均由1～2年實(shí)測歷史數(shù)據(jù)調(diào)和分析而得；

(6)根據(jù)實(shí)測數(shù)據(jù)對河口潮位邊界進(jìn)行校正，以燈籠山為例，實(shí)測值、預(yù)測值、預(yù)測校正值如圖5所示，所有控制站點(diǎn)潮位預(yù)報(bào)誤差統(tǒng)計(jì)見表1、表2所示；

(7)采用加權(quán)移動平均方法推求上游控制站梧州站流量過程，流量預(yù)報(bào)采用15小時(shí)滾動預(yù)報(bào)；

(8)根據(jù)實(shí)測數(shù)據(jù)對上游梧州站流量邊界進(jìn)行預(yù)報(bào)及校正，實(shí)測流量、預(yù)報(bào)流量及經(jīng)校正后的預(yù)報(bào)流量過程如圖6所示。

表1控制站點(diǎn)潮位預(yù)報(bào)誤差統(tǒng)計(jì)(2014.05.01.0:00～5.31.23:00一個(gè)月)

表2控制站點(diǎn)2～4個(gè)月預(yù)報(bào)期潮位預(yù)報(bào)誤差統(tǒng)計(jì)

(9)輸入預(yù)報(bào)的上游梧州站流量過程和河口八大口門潮位過程，進(jìn)行一維非恒定流水動力數(shù)學(xué)模型計(jì)算，求得中順大圍內(nèi)石岐河與外江相接的兩個(gè)主要控制閘——東河水閘、西河水閘預(yù)報(bào)潮位過程，如圖7(a)、圖7(b)所示。

綜上所述，本發(fā)明公開的一種基于水動力數(shù)值模擬的感潮河網(wǎng)潮位預(yù)報(bào)方法，首先，建立感潮河網(wǎng)區(qū)及上游河道一維水動力數(shù)學(xué)模型；其次，利用調(diào)和分析預(yù)測河口控制站潮位過程并利用預(yù)報(bào)前一日實(shí)測數(shù)據(jù)進(jìn)行預(yù)報(bào)值校正；然后采用加權(quán)移動平均方法推求上游控制站點(diǎn)流量過程；最后利用河網(wǎng)區(qū)及上游干流區(qū)一維水動力數(shù)學(xué)模型求解感潮河網(wǎng)區(qū)各斷面潮位過程，實(shí)現(xiàn)河網(wǎng)區(qū)任意斷面潮位預(yù)報(bào)。

上述實(shí)施例為本發(fā)明較佳的實(shí)施方式，但本發(fā)明的實(shí)施方式并不受上述實(shí)施例的限制，其他的任何未背離本發(fā)明的精神實(shí)質(zhì)與原理下所作的改變、修飾、替代、組合、簡化，均應(yīng)為等效的置換方式，都包含在本發(fā)明的保護(hù)范圍之內(nèi)。