基于數字流域系統的漳衛河系洪水演進模型建立

邢 斌，朱曉春，李永慶

（海河水利委員會科技咨詢中心，天津 300170）

洪水的運動規律復雜多變，威脅人們正常的生產、生活，帶來嚴重的生命財產損失和安全隱患。利用現代計算機技術，洪水演進模型能夠根據上游干支流的洪水過程及區間來水，模擬計算下游任意斷面的洪峰流量、出現時間、最高水位及洪水過程，同時能夠模擬水流在蓄滯洪區內的運動規律，為流域防洪規劃、洪水預報、控制工程調度等提供科學的決策依據。

目前，洪水演進模擬方法主要有水文學法、歷史洪災調查法、一維或二維水力學法等。其中，二維水力學模擬方法能夠分別從時間與空間角度模擬洪水演進過程，提供更為直觀的洪水淹沒范圍變化，特別是一維、二維、零維耦合的水力學模型能在大流域尺度內更好地為防洪調度的決策提供服務。

隨著水力學模型技術日趨成熟，建模平臺可以使用丹麥的MIKE 系列、荷蘭的De1ft3D 等商業化軟件，也可以采用我國自主開發的如天津大學、河海大學、清華大學等高校研制的洪水演進模型。其中，河海大學王船海開發的數字流域系統應用開發軟件，因其可視化的建模、需求方案的可定制、計算成果的在線查詢分析、動態顯示等功能，在太湖、長江、淮河等流域有很多成功應用案例，如太湖流域洪水預報系統、三峽水庫動庫容預報調度系統、江蘇省里下河河網水動力模型系統等。筆者采用數字流域系統建立海河流域漳衛河系洪水演進模型，依據漳衛河系的現行調度規劃模擬了漳衛河系發生50年一遇洪水的演進過程。

1 流域概況

漳衛河系是海河流域最南部的防洪骨干水系，流經山西、河南、河北、山東、天津四省一市，流域面積3.76 萬km2，西以太岳山為界，南接黃河、徒駭馬頰河，北界滏陽河，東達渤海。流域地形總的趨勢西高東低，地面坡度山區丘陵區為10‰～0.5‰、平原為0.1‰～0.3‰。

漳衛河系由漳河、衛河、衛運河、漳衛新河、南運河組成。流域上游有漳河和衛河兩大支流。漳、衛兩河于徐萬倉匯合后稱衛運河，至四女寺樞紐，以下為漳衛新河和南運河。漳衛新河是漳衛河水系的入海尾閭河道。主要的蓄滯洪區有良相坡、長虹渠、白寺坡等。

2 模型建立及參數率定

漳衛河流域面積較大，河流眾多，蓄滯洪區分布其中。根據流域實際情況，洪水演進模型需要建立一維、二維及零維模型。

2.1 模型原理

2.1.1 一維模擬

河道洪水波傳播屬不恒定流，其演變規律符合圣維南方程組，并采用偏微分方程的數值解法進行水流的數值模擬，公式如下：

式中：q 為旁側入流（m3/s）；Q，A，B，Z 分別為河道斷面流量（m3/s），過水面積（m2），河寬（m）和水位（m）；Vx為旁側入流流速在水流方向上的分量，一般可以近似為0；K為流量模數，反映河道的實際過流能力；α為動量校正系數，是反映河道斷面流速分布均勻性的系數；g 為重力加速度（m/s2）；t 為時間（s）；x 為距上一斷面的距離（m）。

2.1.2 二維模擬

蓄滯洪區的水流采用二維淺水波方程來描述：

式中：Z為水位（m）；u，v分別為x與y方向上的流速（m/s）；U，V分別為x與y方向上的單寬流量（m3/s）為單寬流量的矢量，為它的模，；q為考慮降雨等因素的源項（m/s）；g 為重力加速度（m/s2）；c為謝才系數；f為柯氏力系數；wxt，wyt分別為風應力沿x和y方向的分量；t為時間（s）；h為水深（m）；ρ為空氣密度（kg/m3）；x，y為自變量坐標。

對上述二維淺水波方程直接求解有一定的困難。因此，采用破開算子法將該方程分成如下兩分步方程組，然后分別對其采用合適的方法進行求解。

第一分步為：

第二分步為：

對上面兩分步方程組的數值求解，采用直角坐標系下非均勻矩形規則網格的控制體積法，具體方法如圖1所示。

圖1 二維差分網格示意

蓄滯洪區邊界條件一般有堤，四周封閉，在大水情況下通過口門、漫堤、閘堰等與河道干流進行水量交換，其交換水量對行洪區講是邊界條件，對整個模型計算講是內邊界條件。

由單元I流進單元J的流量為QX=δX(ZI-ZJ)+βX，單元K流到單元J的流量為QY=δY(ZK-ZJ)+βY，對式（4）離散可得：

式中：A 為單元J 的面積（m2）；∑Qi表示包括降雨在內單位時間內流進單元J水量的代數和（m3/s）；Z為單元J的水位（m）；Δt為計算步長（s）。

2.1.3 零維模擬

計算軟件對無地形資料的蓄洪區、湖泊等對象處理為零維區域，主要考慮水流在該單元區域內的調蓄作用，至于其中的運動規律則不考慮，僅僅考慮區域內的水位變化。水位的變化規律遵循水量平衡原理，流入區域的凈水量等于區域內的蓄量增量，即進出蓄洪區內的水量等于其內的蓄量變化。其計算公式為：

式中： )(zA為蓄洪區內水面面積（m2），一般為水位的某種函數關系；∑Q為蓄變量（m3/s）；Z為t時刻的水位（m）；t為時間（s）。

對式（6）離散后為：

式中：∑Q為蓄變量（m3/s）； )(zA為蓄洪區內水面面積（m2）；Z0為初始水位（m）；Z 為t 時刻的水位（m）；t為時間（s）。

2.1.4 模型耦合

流域洪水運動模擬由零維、一維、二維模擬所組成，各部分模擬必須耦合聯立才能求解，而各部分模擬的耦合是通過“聯系”來實現的。“聯系”要素是指流域中控制水流運動的堰、閘及口門等，其過流流量一般以寬頂堰流量公式計算。

寬頂堰上的水流可分為自由出流、淹沒出流兩種流態，不同流態采用不同的計算公式。當出流為自由出流時，其計算公式為：

當出流為淹沒出流時其計算公式為：

式中：B 為堰寬（m）；Z1為堰上水位（m）；Z2為堰下水位（m）；H0為堰上總水頭（m），H0=Z1-Zd；hs為堰頂水深（m），hs=Z2-Zd；g為重力加速度（m/s2）；m為自由出流系數，一般取m=0.325～0.385；?m為淹沒出流系數，一般取?m＝1.0～1.18；Zd為堰頂高程（m）。

對式（8）—（9）離散后，自由出流流量Q=δZ1Z1+βZ1，淹沒出流流量Q=δZ2(Z1-Z2)。其中，δZ1，δZ2，1Zβ為與Z1，Z2有關的系數，一般常采用時段初水位來計算，有時為了提高計算精度，可采用迭代法計算。

2.2 模型建立

模型的建立主要包括地形及網格剖分、人工控制工程概化、控制條件設置、邊界條件確立等。

2.2.1 模擬范圍

本次涉及的模擬范圍主要包括漳河、衛河出山口以下的河道及蓄滯洪區，即漳河（京廣鐵路橋下—徐萬倉）、衛河（淇門—徐萬倉）、衛運河（徐萬倉—四女寺）、漳衛新河（四女寺—入?？冢? 條主要河道，還包括淇河、湯河、安陽河3 條衛河的主要支流；大名泛區、良相坡、長虹渠（含柳圍坡）、白寺坡、共渠西行洪區、小灘坡、任固坡7個蓄滯洪區。計算范圍概化示意，如圖2所示。

圖2 洪水演進模型建模范圍示意

2.2.2 數字地形與網格剖分

一維模型的地形數據是以1∶2 000 河道實測大斷面資料為基礎，將兩堤坐標、起點距、高程導入斷面編輯器，參與模型計算。斷面編輯器中可以對河道進行主槽、灘地分界，設置糙率。河道斷面編輯器界面，如圖3所示。

圖3 河道斷面編輯器界面

良相坡、長虹渠等蓄滯洪區地形數據是以1∶10 000地形圖提取的散點高程為基礎，以蓄滯洪區圍堤為固定邊界，將x，y，z文件導入軟件，生成不規則三角網，得到模擬區域的地形文件，結果如圖4所示。

本次二維計算采用單元中心的有限體積法進行數值模擬，計算區域采用直角坐標系下非均勻矩形網格進行地形的剖分。網格間距根據模擬區域面積大小、地形變化情況、模擬精度等因素綜合分析后設置，當關注某一區域的水文要素時可以設置網格加密，分為單獨橫向加密和單獨縱向加密以及局部單元格加密。模型模擬區域網格剖分結果，如圖5所示。

小灘坡和任固坡無近期地形資料，利用水位－蓄量關系建立零維模型，以水位節點的形式參與計算，忽略其內部流場。零維要素編輯器界面，如圖6所示。

圖4 蓄滯洪區（以大名泛區為例）地形高程渲染

圖5 蓄滯洪區（以大名泛區為例）地形剖分

圖6 零維要素編輯器界面

2.2.3 工程模擬

工程設施主要是指流域中控制水流運動的堰、水閘及蓄滯洪區分洪口門等。這些設施的過流滿足水力學上的堰流公式，其原理與模型耦合一致。模型中將工程概化為“聯系”。河道上相鄰的斷面之間、不同河道的斷面之間、河道斷面與二維網格之間、不同的二維網格之間均可根據水流方向建立“聯系”，分別模擬河道攔河閘、穿堤涵洞、分洪（閘）口門等水利工程。模擬這些工程時，首先根據本身實際坐標找到在模型中的位置，然后將閘底板高程、閘孔寬度、閘孔高度輸入模型，如果是寬頂堰型，則只需輸入堰頂高程和寬度。為了提高計算穩定性，根據剖分網格尺寸大小考慮將較寬的口門分為若干“聯系”進行模擬，如網格邊長200 m、口門寬度700 m可用3個寬度為200 m、1個寬度為100 m的“聯系”組合模擬。

2.2.4 控制條件設置

水閘、口門等工程的具體控制調度方法在模型中是以控制條件的形式實現的。本次計算的控制條件依據《漳衛河洪水調度方案》編寫。臨時破堤型口門，當水位或者流量達到調度條件時被扒開。由于破堤口門被炸開分洪時并不是瞬間達到最大寬度，模型中考慮了口門寬度發展速度控制，即模擬口門自扒開至最大口門寬度這個過程隨時間變化的實際情況。

2.2.5 邊界條件確立

模型中共有5個洪水入口，分別是漳河入流、良相坡入流、淇河入流、湯河入流和安陽河入流，其設計洪水流量過程作為模型的上邊界條件；1 個洪水出口，即漳衛新河河口，其設計潮位過程作為模型的下邊界條件。模型邊界條件示意，如圖7所示。

圖7 模型邊界條件示意

2.2.6 露灘問題處理

對于水流還沒有流到的露灘區域，需進行特殊處理。目前，比較合理的方法是采用動邊界方法確定出計算域中有水和無水區域的界限進行模擬，但是在程序編譯處理過程中比較復雜困難。本模型作近似簡化處理，通過“凍結法”和最小水深假設的方法把露灘問題處理為固定邊界。

2.3 參數率定

計算模型的主要參數包括用于計算河道流量模數（K）的曼寧糙率系數（n）、計算分洪口門流量的堰流系數（m）。其中，分洪口門一般按照自由出流考慮，堰流系數一般取0.325～0.385，考慮到分洪口門沒有控制閘門，計算中一般按低限考慮，可不設定。因此，模型計算時需要設置的主要參數是糙率。

天然河道的糙率一般與斷面形狀、河床特性等因素有關，變化范圍較大。由于沒有完整的實測資料來率定模型糙率參數，目前采用漳衛河系防洪規劃成果中的設計糙率。模型糙率取值，見表1。

表1 模型糙率（n）取值

3 模擬計算

運行模型之前需要給定初始條件，對于流域模型而言，上游至下游主要洪水流路必須通暢，有條件控制的分洪口門及蓄滯洪區不參加初始條件計算。給定一個高于最高地面高程的水位值進行初始條件計算，當整個模型的水位變幅為0 即初始條件計算達到恒定時開始正式計算。

選取50年一遇設計洪水過程，利用之前制作的地形及模型參數，進行漳衛河系洪水演進模擬，不僅可以獲得蓄滯洪區洪水的淹沒要素，還可以輸出河道水面線、控制站水位（流量）過程、蓄滯洪區控制點水位、蓄滯洪區最大入流、出流等。主要計算結果，見表2。

表2 50年一遇洪水模擬主要成果

洪水演進過程支持實時查詢，即在地圖上點擊任意位置，就可獲得該時刻該位置的水位、水深、流速、水量、地形高程等信息。以洪水入流后160 h為例，蓄滯洪區對應的淹沒范圍、河道水面線、流量過程、水位過程，如圖8—10所示。

圖8 長虹渠淹沒范圍

圖9 衛運河水面線

圖10 四女寺樞紐流量、水位過程（當發生50年一遇設計洪水時）

經過對比，本次研究洪水演進的成果與防洪規劃成果有所不同，主要不同點見表3—4。

表3 河道計算成果對比水位：m；流量：m3/s

表4 蓄滯洪區計算成果對比水位：m；流量：m3/s

本次研究洪水演進成果與防洪規劃成果之所以有所不同，是由于兩次洪水演進采用的口門規模、運用規則不同。本次研究按照現行調度方案編寫控制條件。如，規劃中王灣口門寬度采用800 m、淇門160 m、圈里口門200 m，本次研究分別采用370、380、140 m，底高也不盡相同（見表5）；本次研究長虹渠僅通過淇門入流，而規劃先運用宋村進行分洪；本次研究白寺坡王灣入流大于規劃，是由于共渠西上部李橋水位高于63.471 m 時即通過邢固扒口分洪入白寺坡，白寺坡最大入流為邢固、王灣合計最大，而規劃中邢固口門未啟用?？傮w上講，本次研究洪水演進的結果較為合理，與規劃報告的趨勢保持一致。

表5 規劃與本次研究采用口門規模及控制條件對比

4 結論

筆者以漳衛河系洪水演進為例研究了數字流域系統模擬一維、二維、零維水力學的建模過程，得到以下主要結論。

（1）數字流域系統中的水動力學模塊能夠較好地應用于大型流域河網一維、蓄滯洪區二維及零維洪水演進數值模擬。本次研究以漳衛河山區以下發生50年一遇洪水時按照《漳衛河洪水調度方案》進行調度為模擬對象，與以往成果進行對比驗證，結果顯示模型的模擬精度滿足工作需要，可用于平原河網地區、大型流域的洪水演進分析。

（2）二維模型采用隱格式有限體積法計算，在滿足模型穩定性和精度要求的條件下計算步長可達600 s以上，計算效率較高。較快的計算速度使模型能夠應用于洪水預報、調度方案調整、大型流域防洪規劃等領域。

（3）由于海河流域缺乏完整的實測洪水資料，本次計算僅使用以往規劃成果所采用的模型參數進行模型調試，再與規劃成果進行對比，分析各項條件的異同，判斷計算結果及其趨勢是否合理，因此對模型驗證結果的可靠性支撐略顯欠缺。未來，隨著流域洪水記錄的積累，可對模型進行更為全面的驗證，以進一步提高模擬精度。

[1]田景環，李芳芳.河道洪水演進淺析及一維數學模型的建立[J].中國水運，2007，（8）：94-95.

[2]魏凱，梁忠民，王軍.基于MIKE21 的濛洼蓄滯洪區洪水演算模擬[J].南水北調與水利科技，2013，（6）：16-19.

[3]王船海，李光熾.行蓄洪區型流域洪水模擬[J].成都科技大學學報，1995，（2）：6-14.

[4]周潮洪，杜津媛，常守權.蓄滯洪區二維洪水演進系統的建立[J].海河水利，2002，（6）：33-34.

[5]海河水利委員會.漳衛河系防洪規劃[R].天津：海河水利委員會，2008.