基于大氣-水動(dòng)力集成模型的城市內(nèi)澇模擬預(yù)報(bào)研究

陳光照，侯精明，同 玉，周 聶，高徐軍，蘇 鋒，李繼成，呂 鵬，楊 霄，張宏芳

(1.西安理工大學(xué) 西北旱區(qū)生態(tài)水利國家重點(diǎn)實(shí)驗(yàn)室， 陜西 西安 710048；2.中國電力建設(shè)集團(tuán)有限公司 西北勘測設(shè)計(jì)研究院有限公司， 陜西 西安 710065；3.陜西省氣象局 陜西省氣象服務(wù)中心， 陜西 西安 710014)

近年來，隨著城市化率的提高和氣候變化的加劇，城市內(nèi)澇災(zāi)害愈加頻繁[1]。在下墊面硬化、排水能力不足、短歷時(shí)強(qiáng)降水事件頻發(fā)等因素的共同作用下，許多城市將遭受內(nèi)澇災(zāi)害并造成嚴(yán)重?fù)p失[2]。政府間氣候變化專門委員會(IPCC)的研究表明，1981年至2000年間，500 mm極端降雨的發(fā)生概率約為1%，至2100年這一概率將高達(dá)18%[3]。因此，開展城市內(nèi)澇預(yù)報(bào)方法研究，對防澇減災(zāi)具有重要意義。

城市內(nèi)澇預(yù)報(bào)是降低內(nèi)澇災(zāi)害風(fēng)險(xiǎn)、減輕災(zāi)害損失的重要措施，但目前仍未開發(fā)出完全滿足防災(zāi)減災(zāi)要求的預(yù)報(bào)系統(tǒng)[4]。影響預(yù)報(bào)效果的問題主要為預(yù)見期短及無法精確描述內(nèi)澇變化過程。為了延長預(yù)見期，內(nèi)澇預(yù)報(bào)可采用數(shù)值天氣預(yù)報(bào)(NWP)模型生成的預(yù)報(bào)降雨產(chǎn)品作為地表內(nèi)澇模型的輸入數(shù)據(jù)[5]。例如，Li等[5]集成了柳溪河水文模型與WRF氣象模型，將其用于南方流域大尺度洪澇預(yù)報(bào)。該方法延長了內(nèi)澇預(yù)報(bào)的預(yù)見期，但只能提供網(wǎng)格分辨率為20km×20km的預(yù)報(bào)產(chǎn)品。Tian等[6]開發(fā)了耦合WRF模型、河北水文模型及3DVar數(shù)據(jù)同化模塊，構(gòu)建了大氣-水文預(yù)報(bào)系統(tǒng)。結(jié)果表明，該預(yù)報(bào)系統(tǒng)對空間分布均勻的降雨有較好的預(yù)報(bào)效果，但對短歷時(shí)強(qiáng)降雨預(yù)報(bào)誤差較大。Habibi等[7]采用了由X波段氣象雷達(dá)數(shù)據(jù)驅(qū)動(dòng)的NWS水文分布模型。結(jié)果表明，該系統(tǒng)能較好地預(yù)測水位和河流的變化，但其時(shí)間誤差較為嚴(yán)重。

常規(guī)水文模型無法精確模擬內(nèi)澇變化過程，而求解淺水方程(SWEs)的全水動(dòng)力模型數(shù)值則可解決這一問題[8]?；趧?dòng)力波方法的水動(dòng)力模型能夠模擬水力要素在復(fù)雜地形上的變化過程，生成準(zhǔn)確的內(nèi)澇變化計(jì)算結(jié)果。本文提出集成數(shù)值氣象模型和水動(dòng)力模型構(gòu)建內(nèi)澇預(yù)報(bào)模型，并以此模擬預(yù)報(bào)灃西新城內(nèi)澇淹沒過程，研究結(jié)果將為提升防災(zāi)減災(zāi)能力提供參考借鑒。

1 研究區(qū)域

研究區(qū)位于西咸新區(qū)灃西新城，總面積22.5 km2。全年降雨量主要集中在7～9月。排水管網(wǎng)局部淤堵及結(jié)構(gòu)設(shè)計(jì)不合理，限制了該地區(qū)的排水能力，當(dāng)遇到暴雨時(shí)，內(nèi)澇頻發(fā)。該地區(qū)布設(shè)有微型氣象站，可以氣象站監(jiān)測數(shù)據(jù)作為地面實(shí)況，對GRAPES_MESO預(yù)報(bào)數(shù)據(jù)進(jìn)行評價(jià)。

本文采用無人機(jī)航測技術(shù)對研究區(qū)地形數(shù)據(jù)進(jìn)行了測量。數(shù)字高程模型(DEM)和數(shù)字正射影像圖(DOM)的網(wǎng)格分辨率為2m，約560萬個(gè)網(wǎng)格單元，如圖1、圖2所示。根據(jù)正射影像圖，采用最大似然分類法將網(wǎng)格單元?jiǎng)澐譃槲宸N土地利用類型，如表1所示。曼寧系數(shù)參考相關(guān)標(biāo)準(zhǔn)和城市排水文獻(xiàn)確定[9-10]。

圖1 DEM數(shù)據(jù)Fig.1 DEM data

圖2 土地利用數(shù)據(jù)Fig.2 Land use data

表1 土地利用參數(shù)Tab.1 Properties of the different land uses

2 大氣-水動(dòng)力集成模擬預(yù)報(bào)模型

本文集成GRAPES_MESO大氣模型、二維水動(dòng)力模型和預(yù)報(bào)降雨重構(gòu)方法建立大氣-水動(dòng)力內(nèi)澇預(yù)報(bào)模型，圖3為模型工作流程圖。集成模型每隔12 h將滾動(dòng)更新的預(yù)報(bào)降雨數(shù)據(jù)輸入基于回歸分析方法的重構(gòu)模塊進(jìn)行改進(jìn)，而后將重構(gòu)的預(yù)報(bào)降雨數(shù)據(jù)作為水動(dòng)力模型的輸入數(shù)據(jù)。為了評估重構(gòu)法對預(yù)報(bào)降雨和內(nèi)澇數(shù)據(jù)的改善效果，本文還利用氣象站的實(shí)測降雨數(shù)據(jù)和原始預(yù)報(bào)數(shù)據(jù)驅(qū)動(dòng)內(nèi)澇模型，進(jìn)行內(nèi)澇淹沒過程模擬計(jì)算，并以淹沒面積和水深預(yù)報(bào)結(jié)果分析模型性能。

圖3 城市內(nèi)澇預(yù)報(bào)工作流程Fig.3 Flowchart for the urban inundation forecasting

2.1 GRAPES_MESO氣象模型

全球/區(qū)域同化預(yù)報(bào)系統(tǒng)(GRAPES)是中國氣象局自主開發(fā)的采用多級數(shù)據(jù)同化的數(shù)值天氣預(yù)報(bào)模型。根據(jù)不同分辨率要求，選擇區(qū)域或全球中期物理過程軟件包，可形成全球中期天氣預(yù)報(bào)系統(tǒng)(GRAPES_GFS)或區(qū)域中尺度天氣預(yù)報(bào)系統(tǒng)(GRAPES_MESO)。城市內(nèi)澇預(yù)報(bào)模型采用GRAPES_MESO 4.0版本預(yù)報(bào)降水?dāng)?shù)據(jù)。預(yù)報(bào)系統(tǒng)每天00:00和12:00對中國及周邊地區(qū)的未來72h降雨數(shù)據(jù)進(jìn)行滾動(dòng)更新，時(shí)空分辨率為3h、10km。

2.2 預(yù)報(bào)降雨數(shù)據(jù)重構(gòu)方法

該預(yù)報(bào)降雨數(shù)據(jù)重構(gòu)方法以GRAPES_MESO預(yù)報(bào)數(shù)據(jù)(數(shù)據(jù)來自國家氣象科學(xué)數(shù)據(jù)中心http://data.cma.cn/)為自變量，氣象站實(shí)測數(shù)據(jù)為因變量。利用SPSS統(tǒng)計(jì)分析軟件的回歸分析函數(shù)擬合修正公式。由于GRAPES_MESO的數(shù)值預(yù)報(bào)產(chǎn)品于2015年12月29日00:00首次發(fā)布，因此該方法選取的分析樣本時(shí)間為2016年1月1日00:00至2019年1月1日24:00。在SPSS回歸分析中，在篩去降雨量小于10 mm的場次降雨后，選取剩余降雨數(shù)據(jù)中的37個(gè)典型致澇降雨數(shù)據(jù)作為樣本數(shù)據(jù)。為增強(qiáng)公式的適用性，各樣本數(shù)據(jù)的降雨歷時(shí)、雨強(qiáng)、累計(jì)降雨量都有所不同。此外，選擇4個(gè)不同特征的降雨作為測試數(shù)據(jù)。

研究區(qū)第一和第三季度預(yù)報(bào)雨量小于實(shí)測數(shù)據(jù)，而第二季度預(yù)報(bào)雨量較大。由于第四季度的暴雨次數(shù)較少，因此不進(jìn)行分析。本文選取線性方程、二次方程、三次方程、對數(shù)方程、指數(shù)方程、冪方程6種常用模型擬合公式。去除過度擬合和低相關(guān)性的函數(shù)，得到預(yù)測數(shù)據(jù)與實(shí)測數(shù)據(jù)之間的函數(shù)關(guān)系，如表2所示。第一季度符合對數(shù)函數(shù)，相關(guān)系數(shù)R2為0.81；第二和第三季度符合三次函數(shù)，R2分別為0.91和0.93。其中，xR為原始預(yù)測數(shù)據(jù)；xH為預(yù)測相對濕度；y表示重構(gòu)值。

表2 修正公式Tab.2 Correction formulas

由于預(yù)報(bào)降雨數(shù)據(jù)的時(shí)間分辨率為3 h，當(dāng)降雨持續(xù)時(shí)間小于3 h時(shí)，雨峰將變得不明顯，致使降雨過程變得平緩，進(jìn)而造成較大的誤差。為了使重構(gòu)后的降雨數(shù)據(jù)更符合實(shí)際情況，對采用修正公式(見表2)重構(gòu)的降雨再進(jìn)行重新分布。基于多年降雨數(shù)據(jù)的芝加哥雨型公式可以描述該地區(qū)的短期降雨過程，如式(1)所示。降雨將按芝加哥雨型每小時(shí)的雨量分配比例進(jìn)行分配。根據(jù)公式計(jì)算結(jié)果，在第1小時(shí)分配降雨總量的10%，第2小時(shí)分配81%，第3小時(shí)分配9%。采用這種分配方法，可以在預(yù)報(bào)降雨時(shí)間小于3 h的情況下構(gòu)建雨峰。

(1)

式中：q為暴雨強(qiáng)度；p為重現(xiàn)期；m為暴雨歷時(shí)。

納什效率系數(shù)(NSE)用于評估預(yù)測數(shù)據(jù)與實(shí)際情況的相似程度，ηNSE值越接近于1，則模擬結(jié)果越可靠。

(2)

2.3 二維水動(dòng)力模型

2.3.1控制方程

模型控制方程為耦合水文過程的二維淺水方程(SWEs)，忽略運(yùn)動(dòng)黏性項(xiàng)、科氏力、風(fēng)應(yīng)力及紊流黏性項(xiàng)，其對應(yīng)的二維非線性淺水方程守恒格式的矢量形式為：

(3)

其中，

式中：q為變量矢量，包括x和y方向的單寬流量qx和qy；u、v表示x和y方向的流速；t為時(shí)間；F和G分別為x、y方向的通量矢量；i為凈雨率；S為源項(xiàng)矢量；zb為河床底面高程；Cf=gn2/h1/3為河床糙率系數(shù)，n為曼寧系數(shù)，g為重力加速度，h為水深。

2.3.2數(shù)值方法與GPU高性能計(jì)算技術(shù)

模型采用基于結(jié)構(gòu)化網(wǎng)格的Godunov格式有限體積法求解SWEs。通過二階顯式Runge-Kutta法構(gòu)造具有二階時(shí)空精度的MUSCL型格式，從而確保物質(zhì)守恒并有效解決不連續(xù)問題。針對模擬過程中可能產(chǎn)生急變流及非連續(xù)等復(fù)雜問題，模型選用HLLC近似黎曼求解器對單元界面上的質(zhì)量以及動(dòng)量的通量進(jìn)行求解。底坡源項(xiàng)通過底坡通量法計(jì)算[11]。摩阻源項(xiàng)通過改進(jìn)的顯式方法計(jì)算[12]。模型代碼使用C++和CUDA進(jìn)行編程，從而能夠在圖形處理單元(GPU)實(shí)現(xiàn)并行計(jì)算，大大加快了計(jì)算速度。使用矩陣的方式將等待計(jì)算的單元分配給每個(gè)線程。每個(gè)GPU線程按順序計(jì)算網(wǎng)格邊界通量和源項(xiàng)。

3 結(jié)果及討論

3.1 模型驗(yàn)證

以2016年8月25日00:36至8月25日12:27的實(shí)測降雨資料來驗(yàn)證模型的準(zhǔn)確性。該場次降雨為雙峰型，總歷時(shí)為12h，降雨總量為97.2 mm，如圖4所示。根據(jù)灃西新城管委會提供的監(jiān)測資料，對其中4個(gè)內(nèi)澇點(diǎn)(A、B、C、D)的模擬淹沒面積進(jìn)行分析，如圖5所示。

圖4 降雨過程線Fig.4 Rainfall process line

圖5 內(nèi)澇點(diǎn)分布圖Fig.5 Inundation location map

表3為實(shí)測結(jié)果與模擬結(jié)果的對比。從表3中可以看出，模擬內(nèi)澇面積與實(shí)測數(shù)據(jù)相吻合，相對誤差最大為-4.2%，表明該模型精度較高，可用于大尺度復(fù)雜地形下的城市內(nèi)澇過程模擬。

表3 實(shí)測結(jié)果與模擬結(jié)果比較(t=4 h)Tab.3 Comparison of measured results and simulated results (t=4 h)

3.2 預(yù)報(bào)降雨性能評估

以4場典型降雨(Ⅰ、Ⅱ、Ⅲ和Ⅳ)為例，對模型預(yù)報(bào)降雨性能進(jìn)行評估。詳細(xì)的降雨數(shù)據(jù)信息已在表4中列出。通過圖6所示的原始預(yù)報(bào)降雨、重構(gòu)預(yù)報(bào)降雨與氣象站實(shí)測數(shù)據(jù)對比結(jié)果可以得出，重構(gòu)數(shù)據(jù)的ηNSE值較原始預(yù)報(bào)數(shù)據(jù)的ηNSE值有了顯著提高，說明該重構(gòu)方法是合理的。但在原始預(yù)報(bào)誤差較大的情況下，重構(gòu)方法的修正能力也受到限制，存在著一定的不確定性，如Ⅰ、Ⅲ號降雨。

表4 典型降雨修正結(jié)果Tab.4 Corrected results of the typical rainfall events

圖6 預(yù)報(bào)降雨過程線對比Fig.6 Comparison of forecast rainfall processes

3.3 淹沒深度及面積預(yù)報(bào)性能評估

以2017年8月7日10:18發(fā)生的降雨致澇事件為例，該次模擬所用計(jì)算平臺的CPU為IntelRCoreTMi7-8700，GPU為NVIDIA RTX2080。集成內(nèi)澇預(yù)報(bào)模型耗時(shí)2.45小時(shí)完成了7小時(shí)內(nèi)澇過程的模擬。由于GRAPES_MESO的預(yù)報(bào)降雨產(chǎn)品于2017年8月7日00:00發(fā)布，故此次內(nèi)澇預(yù)報(bào)信息可以提前約9小時(shí)生成。內(nèi)澇點(diǎn)模擬預(yù)報(bào)與現(xiàn)場實(shí)況的對比如圖7所示，結(jié)果表明，預(yù)報(bào)模型可以準(zhǔn)確預(yù)報(bào)出內(nèi)澇位置。

以秦皇島與開元路十字路口(圖7中B點(diǎn))為例，分析該區(qū)域淹沒面積與最大水深的預(yù)報(bào)結(jié)果。圖8、圖9為采用三類降雨數(shù)據(jù)(實(shí)測、原始預(yù)報(bào)和重構(gòu)預(yù)報(bào)數(shù)據(jù))模擬預(yù)報(bào)的內(nèi)澇過程。從圖中可以看出，三類降雨數(shù)據(jù)模擬的面積和水深其變化趨勢是一致的，但具體數(shù)值存在差異。由于GRAPES_MESO生成的原始預(yù)報(bào)數(shù)值高估了降雨量和持續(xù)時(shí)間，故用原始預(yù)報(bào)降雨量驅(qū)動(dòng)預(yù)報(bào)系統(tǒng)得到的結(jié)果會稍大于用實(shí)測降雨模擬得到的預(yù)報(bào)結(jié)果，且峰現(xiàn)時(shí)間延遲了1小時(shí)；而利用重構(gòu)數(shù)據(jù)模擬的預(yù)測結(jié)果更接近于實(shí)測結(jié)果，峰現(xiàn)時(shí)間也與實(shí)測降雨模擬結(jié)果一致，但1小時(shí)至2小時(shí)的模擬數(shù)值超過了使用原始預(yù)報(bào)降雨得出的數(shù)值，這說明重構(gòu)方法仍存在一定的不確定性。

圖7 2017年8月7日現(xiàn)場圖片與模擬結(jié)果對比(t=6 h)Fig.7 Comparison of field pictures and simulation results on August 7, 2017 (t=6 h)

圖8 不同輸入降雨條件下的淹沒面積變化過程Fig.8 Simulated inundation area by using different input data

圖9 不同輸入降雨條件下的最大水深變化過程Fig.9 Simulated max water depth by using different input data

圖8、圖9所示結(jié)果表明，使用重構(gòu)降雨得出的內(nèi)澇預(yù)報(bào)結(jié)果的相對誤差在大多數(shù)時(shí)刻都較低。其預(yù)報(bào)面積和水深的平均相對誤差分別為11.24%和3.75%；預(yù)報(bào)面積和水深的ηNSE分別為0.89和0.94。重構(gòu)方法能夠提升模型的預(yù)報(bào)性能，但重構(gòu)方法改進(jìn)的降雨數(shù)據(jù)存在一定的不確定性，可能會使預(yù)報(bào)結(jié)果出現(xiàn)偏差。同時(shí)，原始預(yù)測數(shù)據(jù)的時(shí)空分辨率不足也是影響預(yù)報(bào)準(zhǔn)確性的另一因素。由于地形平坦，研究區(qū)域內(nèi)的氣候變化不大，因此，具有中尺度時(shí)空分辨率的預(yù)測數(shù)據(jù)可以在該區(qū)域使用。但是致澇強(qiáng)降雨具有局部分布的特點(diǎn)，應(yīng)盡可能使用具有精細(xì)時(shí)空分辨率的預(yù)報(bào)降雨數(shù)據(jù)。

4 結(jié) 論

本文通過集成GRAPES_MESO數(shù)值天氣預(yù)報(bào)模型與水動(dòng)力模型，構(gòu)建了城市內(nèi)澇預(yù)報(bào)模型，并將該模型應(yīng)用于灃西新城區(qū)域。

1) 針對GRAPES_MESO數(shù)值天氣預(yù)報(bào)模型的預(yù)報(bào)產(chǎn)品存在不確定性的問題，基于回歸分析方法，提出了一種重構(gòu)預(yù)報(bào)降雨數(shù)據(jù)的方法。

2) 集成模型能夠準(zhǔn)確預(yù)報(bào)內(nèi)澇點(diǎn)位置，預(yù)報(bào)淹沒面積和水深的平均相對誤差分別為11.24%和3.75%，ηNSE分別為0.89和0.94。

3) 在Ⅱ號典型降雨算例中，集成內(nèi)澇預(yù)報(bào)模型耗時(shí)2.45小時(shí)完成了7小時(shí)內(nèi)澇過程的模擬，內(nèi)澇預(yù)報(bào)信息可于內(nèi)澇發(fā)生前約9小時(shí)發(fā)布。

本文提出的城市內(nèi)澇預(yù)報(bào)模型預(yù)見期較長、準(zhǔn)確性較高，在城市防澇工作中具有巨大的應(yīng)用潛力。內(nèi)澇預(yù)報(bào)的準(zhǔn)確性很大程度上取決于預(yù)報(bào)降雨的質(zhì)量，因此，改進(jìn)和優(yōu)化數(shù)值天氣預(yù)報(bào)模型仍然是提高內(nèi)澇預(yù)報(bào)準(zhǔn)確性最主要和最直接的手段。在下一步研究中，將開發(fā)天氣預(yù)報(bào)模型降尺度方法，以便獲取更高分辨率的降雨數(shù)據(jù)。