一維管網與二維地表雙向耦合的城市暴雨內澇模擬

鄭茂輝， 姚 帥， 周念清， 劉俊兵

（1. 同濟大學 上海防災救災研究所，上海 200092；2. 同濟大學 城市安全風險監測預警應急管理部重點實驗室，上海 200092；3. 同濟大學 水利工程系，上海 200092）

受人類活動和氣候變化雙重影響，國內外極端強降雨頻發，城市內澇問題加?。?-2］，特別是部分沿海城市受海平面上升和地面沉降影響面臨著額外的內澇災害風險［3］。探索準確、高效的城市內澇模擬方法對于預防和減輕城市內澇災害具有重要意義。

城市內澇通常與排水系統能力不足有關［4］。歐美發達國家從20 世紀60 年代開始研制多種城市雨洪模型用來評估排水管網的水力特征，代表性模型有：美國環境保護署的暴雨洪水管理模型（storm water management model，SWMM）［5］、丹麥水力學研究所的城市水文模型（Mike Urban）［6］和城市排水系統模型（MOUSE）［7］、英國Wallingford公司的城市綜合流域排水模型（InfoWorks ICM）［8］等。其中，SWMM 模型免費開源，對管網匯流計算表現良好，應用最為廣泛。一維管網模型不能模擬二維地表積水擴散，結合水位-體積曲線方法能夠確定積水深度及范圍，但對于管網和上覆地面之間復雜的水力交互仍缺乏有效模擬手段。

為準確模擬城市區域暴雨內澇淹沒過程，需要將一維管網模型與二維淹沒模型進行結合［9-10］。一些商業軟件，如MIKE Urban、XP-SWMM 和PCSWMM等，雖然實現了一維/二維水動力耦合模擬，但軟件價格昂貴，且耦合機制未公開，一定程度上限制了其使用。一些學者［11-13］嘗試結合SWMM和開源或半開源二維水動力模型，實現城市暴雨內澇過程模擬。例如，曾照洋等［12］耦合SWMM和半開源二維水動力模型（LISFLOOD-FP）［14］對東莞市某區域進行了洪澇模擬，結果表明耦合模型在研究區具有較好適用性；李鵬等［13］通過連接SWMM 和LISFLOOD-FP模擬了濟南市某流域實測場次降水過程，取得了較好的模擬效果，認為耦合模型可用于城市流域暴雨內澇的模擬。以上研究通過SWMM管網模型和LISFLOOD-FP 淹沒模型的單向連接，將管網節點溢流作為點源邊界條件驅動LISFLOOD-FP 模型，實現過程相對簡單、易用，但忽略了管網節點的地表積水回流，內澇模擬結果可能趨于嚴重。

本文提供了一種基于開源暴雨洪水管理模型接口（PySWMM）［15-16］實現一維管網模型與二維淹沒模型雙向耦合的方法，構建SWMM/LISFLOODFP耦合模型，實現排水管網及其上覆地表的雙向水量交換，并利用上海市浦東外高橋地區的兩次降水序列對耦合模型進行校驗，同時對比單向和雙向2種耦合方式的模擬結果，驗證模型的適用性和精確性。

1 模型方法

所構建的SWMM/LISFLOOD-FP 雙向耦合模型中，SWMM 模型用于產匯流和一維管網模擬計算，LISFLOOD-FP模型則用于地表二維淹沒模擬。一維和二維模型以管網的檢查井為連接點，通過節點溢流和回流實現管網和地表的雙向流量交互，其中節點溢流量通過PySWMM 提取，回流量采用堰流和孔口流量公式確定；一維模型采用固定時間步長，二維模型采用可變自適應步長，二者通過時間同步實現步長級數據交換。

1.1 SWMM管網模型

SWMM 模型是一個降雨徑流水量和水質分析的綜合性計算機模型，由水文、水力及水質三個主要模塊組成，可以實現地表產匯流、管網一維水動力模擬以及水質模擬。其中，地表產流可選擇霍頓（Horton）下滲模型、格林-安普特（Green-Ampt）下滲模型和徑流曲線法（soil conservation service curve number method，SCS-CN）模擬下滲過程，地表匯流采用非線性水庫法，管網一維水動力可選擇運動波或動力波模擬。SWMM 能夠模擬匯水區域、檢查井、管道等水文和水力要素的時空分布［17］，具有適用性良好、開源和應用簡便等優點，在城市排水系統模擬中得到了廣泛應用［18］。本文重點考慮暴雨、徑流、內澇的關系，未將水質模擬納入到此項研究之中。

雨水管網匯流計算采用一維非恒定流模型，通過連續性方程和動量方程聯立求解：

式（1）—（2）中：Q為管道中的流量，m3·s－1；A為過流斷面的面積，m2；S0為管道坡度，P為濕周長，m；n為曼寧摩擦系數；h為水流深度，m。

1.2 LISFLOOD-FP淹沒模型

LISFLOOD-FP 模型是英國布里斯托大學Bates 等人開發的一種基于柵格計算的二維水動力模型［14，19］，該模型可以和高精度的數字高程模型（digital elevation model， DEM）結合，在低海拔洪泛平原區適用性較好，是一種簡單可靠的二維洪水淹沒模型［20］。經過不斷的升級更新，近年來LISFLOOD-FP模型在城市內澇模擬中也得到應用和推廣。

采用LISFLOOD-FP的洪泛區求解器模擬地表匯流，將二維地表水流運動離散到正方形柵格上，用連續性方程和動量方程描述：

式（3）—（4）中：t為時間，s；Qup、Qdown、Qleft、Qright分別為來自上游、下游、左邊和右邊柵格單元的流量，m3·s－1；當前時間步長某柵格單元的水量變化等價于其4 個鄰接柵格的流量累加；Qij為相鄰柵格單元i，j交界處的流量，m3·s－1；Aij為相鄰柵格單元i，j交界處的截面積，m2；Rij為i，j柵格交界處的水力半徑，m；Sij為相鄰柵格i和j之間的水面坡度；n為曼寧系數。

1.3 模型耦合

1.3.1 模型耦合原理

一、二維水動力模型的耦合可分為單向耦合與雙向耦合。單向耦合將一維管網模型（SWMM）的節點溢流過程作為點源邊界條件驅動二維淹沒模型（LISFLOOD-FP）（圖1a），該耦合方式不考慮管網節點回流，數據單向傳遞。雙向耦合則是將一維、二維模型分步模擬，互為驅動（圖1b），該耦合方式不光考慮管網節點溢流，地表積水也可通過節點回流到地下管網，數據雙向交互。

圖1 模型耦合原理及過程Fig. 1 Coupling principle and process of models

美國環境保護署的SWMM 模型不具備在運行過程中讀取模擬結果和交互數據的功能，通常需要修改模型源代碼方能實現其同二維水動力模型的數據互操作，本文選擇基于SWMM5 開發的第三方開源庫PySWMM，實現任意時間步長模擬結果的讀取，通過開發控制算法對節點流量交換實現步長級控制。

雙向耦合模型通過發生溢流或回流的節點連接一維管網模型和二維淹沒模型，實現雙向數據交換。當降雨強度超過管網排水能力時，管網中的水量通過超載節點溢流到上覆地表；相反，當管網排水能力充足時，地面積水通過對應的井點回流到地下管網中。雙向耦合模型主要由4個模塊組成，（1）SWMM一維管網模擬：SWMM模型每完成一個步長模擬，均記錄發生溢流或回流的井點坐標、水深和流量（溢流為正，回流為負），作為點源邊界驅動二維淹沒模型；（2）LISFLOODFP二維淹沒模擬：接收一維管網模型的節點流量，模擬二維地表漫流，更新地表積水深度；（3）一維和二維模型的數據交互：兩個模型通過節點溢流和回流進行雙向流量交換，溢流量通過PySWMM讀取，回流量采用堰流或孔流公式計算；（4）一維和二維模型的時間同步：一維、二維模型以串行方式連接，耦合模型每完成一個時間步長的模擬即進行數據交互，必須保證將兩個模型的時間同步。

1.3.2 雙向流量交互

一維、二維模型通過檢查井連接地下管網和地表柵格單元，進行雙向水量交換。假設管網節點水頭為hw，其對應的地表柵格單元水位為h2d，則有如下垂向流量交換：①hw＞h2d，管網中水流由節點溢出到地表；②hw＜h2d，地表水由節點回流到地下排水管網；③hw=h2d，或地表無積水，則地下管網和地表不存在垂向水量交換。

目前關于地下管網和地表的垂向水流交換的基礎理論尚不成熟，計算方法也較為有限［21］。本文通過PySWMM 編程獲取每個時間步長溢流節點的流量，而節點回流量則采用堰流公式和孔口流量公式計算。假定井口地表高程為hz，概化井口過水面積A=BL，其中B表示過水斷面寬度，L表示過水斷面長度。節點回流計算方法如下：

（1） 當h2d＞hz＞hw時，節點回流量采用自由堰流公式計算：

式中：Q為當前時間步長的節點回流量，m3·s－1；cw為堰流流量系數，取值［0，1］；g為重力加速度，m·s－2。

（3） 當h2d＞hw＞hz且（h2d-hz） ＞L時，節點回流量采用孔口流量公式計算：

式中：co為孔口流量系數，取值［0，1］。

1.3.3 時間同步

一維、二維模型的時間同步是數據雙向交互的重要前提。本研究中管網匯流和地表淹沒模擬采用獨立的動力學模型，其中一維管網模型采用固定的時間步長，二維淹沒模型采用的是可變自適應時間步長，因此需要將一維、二維模型的時間步長進行同步匹配。

如圖2 所示，將耦合模型總模擬時長T0劃分成多個等時長間隔，每個等時長間隔即為一維模型的一個固定時間步長Δt1D，同時對應二維模型的一個模擬批次的模擬時長ΔT2D，每個二維模擬批次包含若干自適應時間步長（ΔT2D=ΣΔt2D）。二維模型運行過程中通過設置“檢查點”屬性，使得下一個批次的模擬計算可以繼承上一個批次模擬結束時的水力狀態，由此實現一維模型固定時間步長和二維模型可變步長的時間同步。

圖2 一維和二維模型時間同步Fig. 2 Time synchronization between 1D and 2D models

式（8）—（9）中：T0表示耦合模型的總模擬時長，s；Δt1D表示一維模型的固定時間步長，s；ΔT2D表示二維模型一個批次的模擬時長，s；Δt2D表示二維模型的自適應時間步長，s。為確保數值模擬的穩定性，自適應時間步長可基于克朗 （Courant-Friedrichs-Lewy，CFL）條件進行估算［22］。

2 應用案例

2.1 研究區與數據

受亞熱帶季風氣候影響，上海雨量充沛，年平均降水量約為1 200 mm，且60 %的雨量集中在5 至9月。本文研究區位于上海市浦東新區北部，長江入?？谀蟼?，整體地形平緩，海拔4～5 m，區域面積為5.56 km2。北側和東側為外環運河，西側和南側分別以主干道楊高北路、洲海路為界（圖3a）。該區為上海外高橋保稅區核心區域，建筑物密集，地面硬化程度高，遭遇強降雨時易因排水能力不足造成積水內澇。

（二）分析電路圖，找到問題與已知條件間的關系。在做電學題時，要在讀清題之后，仔細的分析電路圖，先要正確的判斷電路的串并聯情況。教師可以用去掉一個用電器的方法，教學生正確的判斷電路。若去掉用電器后，電路相互影響則為串聯，若不會相互影響，則為并聯。其次要找清楚電表測量的對象，在看電壓表時，可以觀察他并聯在誰的兩端，就是測誰的電壓。在判斷電流表時，可以去掉電流表，看哪個用電器被影響到就是測誰的電流。最后根據歐姆定律，運用所學的電學公式解題。

圖3 研究區位置與排水設施分布Fig. 3 Location and distribution of drainage facilities in the study area

本文采用研究區2015 年的土地利用數據和排水管網數據， DEM 的分辨率為2 m×2 m?？紤]到建筑物對地面積水的阻擋作用，采用設定高程法，在ARCGIS（地理信息軟件）中將建筑物矢量圖層轉化為2 m×2 m柵格數據，疊加到DEM上，并根據建筑物分布修正地面高程。研究區排水管網設計暴雨重現期為1 年，自2010 年以來未有過系統性更新和改造。管網以航津路為界分為南、北兩個部分，互不相連，汛時主要依靠東部邊界2 個雨水泵站抽排至外環運河，匯入長江口。圖3 b 給出管網概化圖，共包括檢查井925 口，管道987 根、泵站排水口2 個。根據管網連通狀況，將研究區劃分上下2 個一級子匯水區；利用泰森多邊形進一步劃分二級子匯水區，確保每個二級子匯水區對應一個檢查井節點。采用2013 年8 月4 日（20130804）和2013 年9 月13 日（20130913）兩次短歷時暴雨過程對模型進行校準和驗證，降雨數據源自研究區高東雨量站，積水數據采用圖3 a中4個位置調研的積水深度數據。

2.2 模型校驗

模型經驗參數取值參考SWMM用戶手冊［23］，采用“20130804”場次降雨過程以及實測積水數據對管道曼寧系數、霍頓下滲參數、孔口系數、堰流系數等敏感性參數進行率定。利用“20130913”場次降雨驗證模型，該場次降雨強度較大，降雨量主要集中在15：30～17：30，2 h 降雨量占全天總降雨量的97 %。大暴雨導致浦東中北部地區內河水位急劇上漲，其中外高橋內河水位達3.18 m，多條馬路短時積水20～50 cm。

采用雙向耦合模型對“20130913”場次降雨進行模擬，模擬時長為6 h。圖4a給出降雨過程和兩個出水口的流量過程線。與降雨雙峰曲線不同，2個出水口流量過程基本上均呈單峰特征，漲洪段較陡，流量峰值出現在第2次雨量峰值后30～60 min，在降雨停止后流量持續走弱。提取圖4b中4個監控點位的最大水深數據，與調研水深對比以驗證模擬結果。由表1可知：模擬水深與調研結果較為接近，說明模型在該片區具有較好的適應性，模擬結果較為可靠。

表1 模擬積水深度與實測結果對比Tab. 1 Comparison of simulated ponding depth with measured depth

圖4 研究區“20130913”場次降雨模擬結果Fig. 4 Measured rainfall process and simulated water outlet flow process

2.3 模擬結果分析

分別采用單向、雙向耦合方式對“20130913”場次降雨進行模擬，對比2 種耦合方式的積水面積和水深分布，并對2 種耦合方式的模擬計算效率進行了比較。

2.3.1 最大積水分布范圍對比

圖5a、圖5b給出單向和雙向耦合的最大積水分布圖，兩者最大積水面積分別為27.21×104m2和21.07×104m2，交集區域（圖5c）面積達17.67 ×104m2，2種耦合方式對內澇淹沒范圍的預測重疊率較高。由表2不同積水深度的面積統計結果發現，對于占比80 %以上、水深0.2 m以下的輕度積水區域，單、雙向耦合的模擬積水面積相近，前者約為后者的1.21倍，但對于中等（0.2～0.3 m）和嚴重（＞0.3 m）積水區域，該比值分別達到1.88和2.1?？梢姡鄬τ陔p向耦合，單向耦合的模擬積水范圍，特別是對于水深0.2 m以上積水區域的預測結果趨于嚴重。此外，盡管單向耦合的模擬積水面積為雙向耦合的1.3倍，但前者出現溢流的節點卻略少于后者。這應與單向耦合未考慮管網節點回流有關，即降低了排水管網超載的壓力。

圖5 單向和雙向耦合最大積水分布圖Fig. 5 Maximum water accumulation area

2.3.2 積水面積與水深動態變化

為定量分析2種耦合方式模擬的積水面積和水深動態變化，圖6 給出了積水面積和最大水深隨時間的變化曲線。結果表明：

圖6 積水面積與最大水深隨時間變化曲線Fig. 6 Time-varying curves of stagnant water area and maximum water depth

（1）在2次雨峰之后即17：00左右，模擬積水面積達到峰值；其后，單向耦合的積水面積基本保持穩定，雙向耦合則更如實反映了地表積水的消退過程（圖6a）。

（2）不同于積水面積過程曲線，水深變化曲線具有明顯的雙峰特征，與降雨過程相符（圖6b）。單向和雙向耦合方式模擬最大水深分別為0.79 m 和0.61 m，水深峰值時刻出現在雨峰后的20～30 min，且單向略滯后于雙向（表3）。

表3 單向和雙向耦合模型結果對比Tab. 3 Comparison of result between unidirectional and bidirectional coupling models

（3） 17：30 后單向耦合最大模擬水深基本保持在0.3 m，雙向耦合最大水深則在18：00 后降至0.2 m以下，并逐漸消減。對照2種方式模擬結果 ，雙向耦合模擬的積水范圍與水深變化與研究區9.13 暴雨內澇報道情形基本一致，能夠更好揭示城市內澇積水以及擴散、消退的全過程。

這一差異增加了雙向耦合的一維模擬計算量。不過考慮到雙向耦合數據交互耗時與計算總時長及步長大小有關，對耦合模型總體運行效率的影響尚需進一步探討。

2.3.3 單、雙向耦合計算效率對比

2種耦合方式固定時間步長均設為30 s、模擬時長為6 h，在同一工作站配置計算環境下對模型運行耗時進行對比。由表4結果可知：在一維模擬階段，雙向耦合的計算耗時是單向耦合的1.8 倍；但在二維模擬階段，雙向耦合計算耗時要略低于單向耦合；總模擬耗時雙向耦合稍大于單向耦合，約為后者的1.1 倍。2 種耦合方式的計算耗時差異原因在于模型運行原理不同：單向耦合模式下一、二維模型可獨立運行，一維模型中管網節點溢流驅動二維水力計算；雙向耦合模式中一維和二維模型需要進行數據雙向交互，互為驅動。

表4 耦合模型計算時間對比Tab. 4 Comparison of calculation time between two coupling models

3 結語

城市下墊面建筑物和基礎設施密布、排水系統立體式分層結構顯著，使得城市的產匯流機制遠比天然流域復雜。為快速、準確模擬城市區域暴雨產匯流過程，本文提供了一種一維管網與二維地表的動態水力交互方法，構建了SWMM/LISFLOODFP雙向耦合模型。以上海外高橋地區為例，采用兩場次降雨過程對耦合模型進行了校準和驗證，并對比了單向耦合和雙向耦合的模擬結果，主要結論如下：

（1）所構建的雙向耦合模型通過時間同步實現了一維/二維模型之間步長級的流量交換控制，能夠較好地模擬地下排水管網與上覆地表之間復雜的水力交互。對實測場次降雨過程的模擬結果表明，該耦合模型具有較好的精度，在研究區適用性良好。

（2）比較單向耦合和雙向耦合模擬結果發現，對于占淹沒區域80 %以上的輕度（＜0.2 m）積水區二者模擬積水面積接近，對于水深0.2 m 以上的積水區，單向耦合模擬結果則相對趨于嚴重，約為后者的1.9倍。降雨結束后，雙向耦合的積水面積和水深隨著積水消退而逐漸減小，更好揭示暴雨積澇的全過程。在計算效率上，雙向耦合增加了數據交互耗時，但二維模擬階段用時略有下降，后續需結合時間步長設置進一步討論，并通過圖形處理器（GPU）高性能計算進一步提升運算效率。

（3）值得補充說明的是，本文采用開源PySWMM 和半開源LISFLOOD-FP 構建了雙向耦合模型，該方法同樣可選擇其他水動力模型進行雙向耦合，耦合模型及成果可用于城市區域暴雨內澇數值模擬及推演，為城市防洪和內澇治理措施的制定提供科學依據。

作者貢獻聲明：

鄭茂輝：耦合模型總體設計，論文修改與定稿。

姚 帥：算法實現，論文初稿撰寫。

周念清：指導論文框架，協助論文修改。

劉俊兵：協助數據處理與分析。