基于BP與PSO的河口雙邊界條件逆向推求方法

易 梓 楊，酈 建 鋒，錢 進，陸 卞 和，張 語 航，李 豐 鐸

(1.河海大學(xué) 環(huán)境學(xué)院，江蘇 南京 210098； 2.河海大學(xué) 淺水湖泊綜合治理與資源開發(fā)教育部重點實驗室，江蘇 南京 210098； 3.中國電建集團華東勘測設(shè)計研究院有限公司，浙江 杭州 311122； 4.浙江省華東生態(tài)環(huán)境工程研究院，浙江 杭州 311122)

0 引 言

臺風(fēng)期的強降雨與天文潮會給中國東部沿海河口地區(qū)帶來嚴(yán)重的洪澇災(zāi)害，造成大量經(jīng)濟損失[1-2]。為保證河口地區(qū)的防洪安全，水利部門經(jīng)常會采用數(shù)學(xué)模型來預(yù)測河口地區(qū)的水文過程。MIKE 11模型以其簡單干凈的數(shù)據(jù)處理界面、可靠的邊界條件輸入、模擬板塊的多樣化等特點在水文模擬中被廣泛應(yīng)用[3-5]。在河口地區(qū)水文模擬的過程中，一般是以上游流量和下游潮位為出發(fā)點，探究雙邊界條件對河口地區(qū)的影響[6-9]，而以限制河口計算區(qū)間某一斷面水位為出發(fā)點，逆向推求雙邊界(上游流量、下游潮位)條件的研究較少。然而在逆向推求河口雙邊界條件時，基于MIKE 11軟件的傳統(tǒng)試算法是通過手動調(diào)節(jié)雙邊界條件來完成該任務(wù)，搜索隨機性、效率、時間成本成為傳統(tǒng)試算法亟待解決的問題。

反饋式神經(jīng)網(wǎng)絡(luò)模型，即BP(Back Propagation)神經(jīng)網(wǎng)絡(luò)模型是最成熟及應(yīng)用范圍最廣的人工神經(jīng)網(wǎng)絡(luò)模型，適合解決模擬軟件計算效率低的問題[10]。粒子群優(yōu)化(Particle Swarm Optimizer，PSO)算法是一種新型的全局隨機優(yōu)化算法，其優(yōu)點是結(jié)構(gòu)簡單、全局搜索能力強、魯棒性好、易于實現(xiàn)和優(yōu)化效率高[11]。BP神經(jīng)網(wǎng)絡(luò)和PSO算法以其網(wǎng)絡(luò)結(jié)構(gòu)優(yōu)良和搜索效率高在各種復(fù)雜的優(yōu)化問題中得到廣泛應(yīng)用。例如，趙文剛等[12]基于BP神經(jīng)網(wǎng)絡(luò)與PSO算法通過相關(guān)性分析和因子貢獻率分析，建立西洞庭湖南咀水文站月平均徑流預(yù)報模型，實現(xiàn)了在保障擬合精度的前提下確定影響徑流預(yù)報因子最少個數(shù)，相比于傳統(tǒng)試算法，提高了計算效率；何勝男等[13]以調(diào)蓄池截污總量最大為優(yōu)化目標(biāo)，以污水處理廠處理規(guī)模和可用征地為約束條件，耦合BP神經(jīng)網(wǎng)絡(luò)與PSO算法推求各調(diào)蓄池的優(yōu)化截污率，實現(xiàn)了滿足占地面積和污水處理廠處理規(guī)模下的截污效率最大化，極大節(jié)約了計算時間。

本文基于BP神經(jīng)網(wǎng)絡(luò)與PSO算法，提出一種在河口某一斷面水位確定條件下逆向推求雙邊界(上游流量和下游潮位)條件的優(yōu)化方法，并以晉江河口為例，將該方法的計算結(jié)果與傳統(tǒng)試算法的計算結(jié)果進行對比分析。

1 傳統(tǒng)的河口雙邊界條件逆向推求過程

MIKE 11水動力模型由于特有的數(shù)據(jù)輸入格式與數(shù)據(jù)集成平臺的限制，只能根據(jù)輸入邊界得到輸出結(jié)果，為得到河道水位首先需要將上、下游邊界的數(shù)據(jù)修改成MIKE 11可識別的dfs0時間序列文件，然后將時間序列文件導(dǎo)入MIEK 11的邊界文件，最后運行模型得到河口沿程水面線計算結(jié)果文件[14-16]，計算效率較低。

因而，以河口某一斷面限制水位為約束條件，采用傳統(tǒng)試算法逆向推求雙邊界條件，首先需確定河口下游高潮位，接著根據(jù)下游高潮位對上游最大流量不斷地進行調(diào)整，并且在MIKE 11模型中一一計算，直到求得滿足該斷面水位條件的上游最大流量，最后輸出上游流量、下游潮位雙邊界條件，具體計算流程如圖1所示。由于下游高潮位是一個變化的區(qū)間，因此試算法需要根據(jù)下游高潮位的變化反復(fù)調(diào)整上游最大流量，搜索效率低，靈活性較差。

圖1 傳統(tǒng)試算法逆向推求河口雙邊界條件流程Fig.1 Flow chart of traditional trial algorithm for inverse derivation of estuarine double boundary conditions

2 基于BP神經(jīng)網(wǎng)絡(luò)和粒子群優(yōu)化算法的河口雙邊界條件逆向推求過程

MIKE 11水動力模型從輸入雙邊界條件到輸出某一斷面水位結(jié)果是一個繁瑣的過程，因此采用BP神經(jīng)網(wǎng)絡(luò)建立輸入雙邊界條件與某一斷面水位的非線性映射關(guān)系。非線性映射關(guān)系建立后，從邊界條件輸入到水位結(jié)果輸出的時間大大減少。同時為了改進試算法手動修改邊界、手動計算不同工況等重復(fù)性工作，引入粒子群優(yōu)化算法，以某一斷面水位為優(yōu)化目標(biāo)，根據(jù)下游高潮位的變化自動搜索上游最大流量進而提高搜索效率。

2.1 BP神經(jīng)網(wǎng)絡(luò)建立輸入與輸出的非線性映射關(guān)系

BP神經(jīng)網(wǎng)絡(luò)是人工神經(jīng)網(wǎng)絡(luò)的一種，主要由輸入層、隱含層和輸出層組成，具有多層前饋網(wǎng)絡(luò)和誤差反向傳播的學(xué)習(xí)機制[17-18]，圖2為建立的河口雙邊界與某一斷面水位的BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)圖。BP神經(jīng)網(wǎng)絡(luò)具有高度的非線性映射能力和較強的泛化能力，其基于梯度下降法原理，預(yù)測精度高[19-20]。只要隱含層和隱節(jié)點足夠多，可以以任意精度逼近一個非線性映射關(guān)系[21]。

圖2 BP神經(jīng)網(wǎng)絡(luò)模型結(jié)構(gòu)Fig.2 Diagram structure of BP neural network model

BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)的擬合精度受樣本數(shù)量和樣本代表性的影響。樣本數(shù)量根據(jù)實際情況確定。為提高樣本數(shù)據(jù)的代表性，對雙邊界條件采用拉丁超立方抽樣的方法，防止抽取數(shù)據(jù)過于集中[22]；接著將抽取的樣本用MIKE batch模塊并行計算；最后利用MATLAB提取指定斷面的水位結(jié)果，進而得到雙邊界條件和指定斷面水位的離線數(shù)據(jù)庫。

BP神經(jīng)網(wǎng)絡(luò)的訓(xùn)練函數(shù)選取的是L-M(Levenberg-Marquardt)算法的training函數(shù)，設(shè)置BP神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)率μ，訓(xùn)練誤差為δ。將離線數(shù)據(jù)庫分為訓(xùn)練集和驗證集。BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)的建立分為兩部分：① BP神經(jīng)網(wǎng)絡(luò)的訓(xùn)練。通過BP神經(jīng)網(wǎng)絡(luò)各層之間的權(quán)值與閾值，不斷調(diào)整訓(xùn)練集中預(yù)測輸出與期望輸出的偏差，直到誤差小于預(yù)設(shè)值δ，訓(xùn)練結(jié)束。② BP神經(jīng)網(wǎng)絡(luò)的驗證。將驗證集的雙邊界條件帶入BP神經(jīng)網(wǎng)絡(luò)得到預(yù)測輸出，對比預(yù)測輸出與期望輸出的差異以評價網(wǎng)絡(luò)結(jié)構(gòu)的優(yōu)良性。采取MSE(均方誤差)、R2(相關(guān)系數(shù))來評判BP神經(jīng)網(wǎng)絡(luò)的擬合效果。

2.2 粒子群優(yōu)化算法實現(xiàn)河口雙邊界條件的逆向推求

根據(jù)BP神經(jīng)網(wǎng)絡(luò)建立的河口雙邊界條件與指定斷面水位的映射關(guān)系，引入粒子群優(yōu)化算法，實現(xiàn)對邊界條件的自動搜尋，提高搜索效率[23-24]。PSO是科學(xué)家研究鳥類捕食行為而得出來的一種算法結(jié)構(gòu)，其原理為：PSO初始化為一組隨機粒子(隨機解)，在每次搜尋過程中粒子通過兩個極值(Pbest，Gbest)來更新自己的位置[25]，各個粒子根據(jù)對比自身最優(yōu)解與群體最優(yōu)解來調(diào)整自己的位置和速度，直至到達預(yù)設(shè)的精度或者預(yù)設(shè)的算法迭代次數(shù)[26]。

通過對PSO變量中初始種群數(shù)、學(xué)習(xí)因子、迭代次數(shù)以及對河口某一斷面水位(優(yōu)化目標(biāo))的設(shè)定，構(gòu)建PSO算法結(jié)構(gòu)。設(shè)置河口某一斷面水位為優(yōu)化目標(biāo)，初始化的粒子隨機分布在樣本空間中，粒子不斷調(diào)整自身空間位置以搜索滿足優(yōu)化目標(biāo)的最優(yōu)解。

2.3 優(yōu)化方法與傳統(tǒng)試算法對比

圖3為優(yōu)化方法與傳統(tǒng)試算法的計算流程對比圖。優(yōu)化方法相比于傳統(tǒng)試算法，主要從兩個方面進行了改進：① 基于BP神經(jīng)網(wǎng)絡(luò)構(gòu)建河口雙邊界條件與某一斷面水位的非線性映射關(guān)系，從邊界條件的輸入到斷面水位結(jié)果的輸出只需要幾秒，提高了傳統(tǒng)試算法的計算效率；② 基于粒子群優(yōu)化算法，通過優(yōu)化目標(biāo)的設(shè)定，實現(xiàn)雙邊界條件的自動搜索，提高了傳統(tǒng)試算法的搜索效率。

圖3 優(yōu)化法與傳統(tǒng)試算法計算流程對比Fig.3 Comparison of the calculation processes of the optimization method and the traditional trial method

3 實例應(yīng)用

3.1 研究區(qū)域概況

晉江位于福建省東南沿海，是中國臺風(fēng)登陸的主要路線之一，臺風(fēng)活動極為頻繁，多發(fā)生在7～9月。如遭遇天文大潮，其危害性更大。晉江作為主要的泄洪通道，在雨水頻繁期，承載著巨大的泄洪壓力。查閱晉江水文資料，金雞水閘100 a一遇的洪峰流量為11 100 m3/s；前浦站100 a一遇的最高潮位為5.16 m，多年平均最高潮位為4.00 m。本次研究范圍從金雞水閘到前浦站，共19.30 km，具體如圖4所示。

圖4 研究區(qū)域概況Fig.4 Overview of the research area

為了保證晉江河口兩岸人民的安全，需要限制防洪安全空間不足地區(qū)的最高水位。根據(jù)福建省水利水電勘測設(shè)計研究院2020年編制的《晉江下游防洪標(biāo)準(zhǔn)復(fù)核報告》可知，晉江河口地區(qū)民俗主題文化公園處的防洪措施最為薄弱。在100 a一遇洪水和100 a一遇潮位的情形下，民俗主題文化公園處計算水位為7.832 m。為了給民俗主題文化公園預(yù)留足夠的防洪空間，本文以民俗主題文化公園水位6.832 m為優(yōu)化目標(biāo)，對雙邊界條件進行逆向推求。

3.2 水動力模型構(gòu)建

以MIKE 11為計算平臺，對研究區(qū)域建立水動力模型，概化河長為19.30 km，從金雞水閘到前浦站，共概化斷面55個。為保證模型計算時的穩(wěn)定性，根據(jù)河道水位的實測值，設(shè)置初始水位為2.00 m。

泉州大橋(橋下)水位站為國家基本水文站，有逐日水位觀測數(shù)據(jù)。晉江下游沿程設(shè)有水尺，自動化記錄水位數(shù)據(jù)。2016年9月15日為晉江下游近年來最大洪水，2019年6月3日為當(dāng)年最大洪水，根據(jù)水文站及水尺的記錄數(shù)據(jù)，采用2016年9月14～16日測量水位數(shù)據(jù)對建立的水動力模型進行率定，選擇2019年6月2～5日測量水位數(shù)據(jù)對率定的水動力模型進行驗證。模型率定的結(jié)果表明：河道主槽糙率在 0.036～0.040之間，灘地糙率在0.042～0.050之間。模型驗證的結(jié)果表明：河道沿程水位的實測值與計算值擬合較為良好，水位絕對誤差小于5 cm(見表1)；泉州大橋(橋下)水位站實測水位與計算水位的絕對誤差小于8 cm(見圖5)。水利部頒布的SL 483-2017《洪水風(fēng)險圖編制導(dǎo)則》[27]規(guī)定，水位允許絕對誤差不應(yīng)超過20 cm。因此，本文建立的水動力模型的精度滿足實際應(yīng)用需要。

表1 沿程水位計算值與實測值比較Tab.1 Comparison of calculated and measured values of water levels along the course m

圖5 泉州大橋(橋下)水位站水動力模型驗證結(jié)果Fig.5 Hydrodynamic model validation results of Quanzhou Bridge (under bridge)water level station

3.3 計算結(jié)果

3.3.1傳統(tǒng)試算法計算結(jié)果

為了保證民俗主題文化公園最高水位小于6.832 m，通過已率定驗證后的MIKE 11水動力模型逆向推求上游流量和下游潮位邊界條件。選取前浦站100 a一遇最高潮位5.16 m、50 a一遇最高潮位5.04 m和多年平均高潮位4.00 m為典型情況逆向推求上游最大下泄流量。逆向推求過程及結(jié)果如表2所列。為保證民俗主題文化公園處于防洪安全位置，下游最高潮位為5.16 m時，上游最大流量為8 160 m3/s；下游最高潮位為5.04 m時，上游最大流量為8 300 m3/s；下游最高潮位為4.00 m時，上游最大流量為9 210 m3/s。由表2可以看出，試算法計算精度有限，隨機性大，花費的時間長，搜索效率較低。

表2 傳統(tǒng)試算法逆向推求雙邊界條件結(jié)果Tab.2 Results of traditional trial algorithm for inverse derivation of double boundary conditions

3.3.2優(yōu)化方法計算結(jié)果

根據(jù)晉江水文資料，選取100 a一遇最高潮位與多年平均潮位作為潮位樣本空間最大值與最小值，選取100 a一遇最大洪峰流量作為流量樣本空間的最大值，流量樣本空間在不影響計算效率的同時應(yīng)盡可能廣，以對應(yīng)下游潮位的變化，根據(jù)計算設(shè)定流量樣本空間的最小值為7 000 m3/s。因此設(shè)置樣本空間為上游最大洪峰流量7 000～11 100 m3/s，下游最高潮位為4.00～5.16 m。利用拉丁超立方抽樣在區(qū)間范圍內(nèi)抽取1 000 組上、下游邊界，然后采用MIKE batch模塊將得到的1 000組邊界數(shù)據(jù)進行并行計算，利用MATLAB批量提取民俗主題文化公園水位結(jié)果，進而建立雙邊界條件與民俗主題文化公園水位的離線數(shù)據(jù)庫?；贛ATLAB平臺搭建BP神經(jīng)網(wǎng)絡(luò)，具體參數(shù)設(shè)置為：輸入層2個，隱含層9個，輸出層1個，迭代次數(shù)1 000 次，學(xué)習(xí)率0.001，訓(xùn)練誤差達到10-6級。

選取750組數(shù)據(jù)作為訓(xùn)練集，250組數(shù)據(jù)作為驗證集。BP神經(jīng)網(wǎng)絡(luò)預(yù)測水位輸出結(jié)果與期望水位輸出結(jié)果如圖6所示，MSE=1.36×10-6、R2為0.996。由此判斷BP神經(jīng)網(wǎng)絡(luò)擬合效果較好，預(yù)測精度高。

圖6 BP神經(jīng)網(wǎng)絡(luò)預(yù)測輸出與期望輸出水位對比Fig.6 Predicted output versus desired output water level by BP neural network

粒子群優(yōu)化過程中，慣性權(quán)重C1=0.8，自我學(xué)習(xí)因子C2=2.0，群體學(xué)習(xí)因子C3=2.0，迭代次數(shù)為400次，種群規(guī)模為30。搜索空間為上游最大流量7 000～11 100 m3/s，下游最高潮位4.00～5.16 m。潮位搜索間隔為0.04 m，總共29種工況。

以民俗主題文化公園最高水位6.832 m為優(yōu)化目標(biāo)，逆向推求滿足優(yōu)化目標(biāo)的雙邊界條件映射關(guān)系。計算得到的上游最大下泄流量與下游高潮位之間映射關(guān)系如圖7所示。根據(jù)PSO算法計算結(jié)果，在100 a一遇最高潮位為5.16 m時，上游最大下泄流量為8 181.00 m3/s；在50 a一遇最高潮位為5.04 m時，上游最大下泄流量為8 332.51 m3/s；在多年平均最高潮位為4.00 m時，上游最大下泄流量為9 222.72 m3/s。對上游最大下泄流量隨下游高潮位的變化關(guān)系選擇二次函數(shù)進行擬合，x為下游高潮位，y為上游最大下泄流量。最終擬合公式為y=6195.1091+2028.2287x-318.2523x2，R2計算值為0.999 87，擬合效果較好。擬合公式可以很好地反映上游最大下泄流量隨下游最高潮位的變化情況，可為晉江金雞水閘的調(diào)度提供依據(jù)。根據(jù)下游潮位的變化，金雞水閘可相應(yīng)調(diào)整最大下泄流量，以保證民俗文化主題公園處于防洪安全空間，實用參考性強。

圖7 上游最大下泄流量隨下游最高潮位變化Fig.7 Variation of upstream maximum flow with downstream maximum tide level

3.4 對比分析

選取下游最高潮位為100 a一遇、50 a一遇和多年平均3種典型情形進行計算，傳統(tǒng)試算法與優(yōu)化方法計算結(jié)果對比如表3所列。在3種情形下，傳統(tǒng)試算法計算民俗主題文化公園水位與目標(biāo)水位均存在一定偏差，但優(yōu)化法計算水位與目標(biāo)水位一致，在計算精度上優(yōu)化法略大于傳統(tǒng)試算法；傳統(tǒng)試算法推求過程中實驗次數(shù)隨機性強，需要人工操作時間長，逆向推求一組滿足優(yōu)化目標(biāo)的雙邊界條件所需時間約為10 min，優(yōu)化方法逆向推求一組滿足優(yōu)化目標(biāo)的雙邊界的時間約為1 min，在計算效率上優(yōu)化法相比于傳統(tǒng)試算法提升了約10倍。因此，優(yōu)化方法計算精度和計算效率均優(yōu)于傳統(tǒng)試算法，更適合解決河口地區(qū)雙邊界的逆向推求。

表3 傳統(tǒng)試算法與優(yōu)化法計算結(jié)果對比Tab.3 Comparison of the results of the traditional trial method and the optimization method

4 結(jié) 語

本文為解決傳統(tǒng)試算法逆向推求河口雙邊界條件(上游流量、下游潮位)計算效率與搜索效率低、時間成本高的問題，提出基于BP神經(jīng)網(wǎng)絡(luò)與PSO算法逆向推求邊界條件的優(yōu)化方法。研究表明：基于BP神經(jīng)網(wǎng)絡(luò)與PSO算法的優(yōu)化方法相比于傳統(tǒng)試算法，在保障計算精度的前提下，計算時長約減少至原來的十分之一。優(yōu)化方法科學(xué)、高效且適用范圍廣，更適合河口地區(qū)雙邊界條件的逆向推求，同時可為河口地區(qū)防洪調(diào)度提供技術(shù)支持。