牛 帥，崔信民

(南京水利科學(xué)研究院水文水資源研究所，南京 210029)

隨著我國(guó)城鎮(zhèn)化進(jìn)程的加快，村鎮(zhèn)地區(qū)社會(huì)經(jīng)濟(jì)水平發(fā)展迅速，洪水災(zāi)害造成的損失日益嚴(yán)重。在工程與非工程措施并重的防洪減災(zāi)新思路下，對(duì)村鎮(zhèn)地區(qū)進(jìn)行動(dòng)態(tài)洪水風(fēng)險(xiǎn)分析，可以為村鎮(zhèn)地區(qū)的水利規(guī)劃、防洪減災(zāi)等提供關(guān)鍵信息參考，具有非常重要的現(xiàn)實(shí)意義。

近年來，國(guó)內(nèi)外學(xué)者對(duì)村鎮(zhèn)地區(qū)的防洪標(biāo)準(zhǔn)、洪水風(fēng)險(xiǎn)、洪災(zāi)損失等內(nèi)容進(jìn)行了廣泛的研究[1,2]，并取得了較大進(jìn)展。王高丹等[3]利用多源數(shù)據(jù)通過GIS方法、層次分析法進(jìn)行了海南島中部村鎮(zhèn)地區(qū)山洪災(zāi)害風(fēng)險(xiǎn)分析。申邵洪等[4]采用一維、二維耦合水動(dòng)力學(xué)模型開展了?？悼h城鎮(zhèn)區(qū)域的洪水影響及淹沒分析研究。練繼建等[5]采用二維水動(dòng)力模型分析了萬泉河下游農(nóng)村地區(qū)洪水淹沒情況,繪制了脆弱性分布圖。由于洪水來源的不確定性，傳統(tǒng)的靜態(tài)洪水風(fēng)險(xiǎn)圖不能滿足實(shí)際需要，在實(shí)際防汛中需要進(jìn)行動(dòng)態(tài)洪水風(fēng)險(xiǎn)分析[6,7]。吳濱濱等[8]研發(fā)了大陸澤、寧晉泊洪水風(fēng)險(xiǎn)實(shí)時(shí)分析系統(tǒng)?？芗维|等[9]結(jié)合水動(dòng)力學(xué)模型和WebGIS技術(shù),構(gòu)建了洪澤湖周邊滯洪區(qū)暴雨內(nèi)澇預(yù)報(bào)預(yù)警系統(tǒng)。

本文根據(jù)村鎮(zhèn)易澇區(qū)的暴雨內(nèi)澇特點(diǎn)，采用一、二維水動(dòng)力學(xué)模型構(gòu)建了易澇區(qū)暴雨內(nèi)澇耦合計(jì)算模型，模型可以根據(jù)未來預(yù)測(cè)降雨過程和預(yù)測(cè)潮位過程進(jìn)行易澇區(qū)內(nèi)澇風(fēng)險(xiǎn)的動(dòng)態(tài)分析，可以為村鎮(zhèn)地區(qū)的實(shí)際防汛工作提供及時(shí)有效的決策參考信息。

1 研究區(qū)域概況

萬城鎮(zhèn)易澇區(qū)位于海南省萬寧市，瀕臨南海，地理位置如圖1所示。萬城鎮(zhèn)轄居10個(gè)社區(qū)，32個(gè)村委會(huì)，2個(gè)鎮(zhèn)辦農(nóng)場(chǎng)，擁有東星工業(yè)開發(fā)區(qū)和烏場(chǎng)港經(jīng)濟(jì)開發(fā)區(qū)。易澇區(qū)處于海南省暴雨高發(fā)區(qū)，暴雨多、強(qiáng)度大。易澇區(qū)內(nèi)河網(wǎng)水系較多，河道下游連接小海，小海洪潮受潮汐通道束縛影響，水位易于高漲，形成頂托，排水不暢。根據(jù)近年來統(tǒng)計(jì)數(shù)據(jù)，年平均受淹造成不同程度損失達(dá)7～8次之多(下游農(nóng)村區(qū)受淹次數(shù)居多)。

圖1 易澇區(qū)位置示意圖Fig.1 Geographical position of waterlogged area

2 易澇區(qū)洪水分析模型

2.1 河網(wǎng)一維水流模型

采用一維水動(dòng)力學(xué)模型對(duì)易澇區(qū)內(nèi)河網(wǎng)的水流運(yùn)動(dòng)進(jìn)行模擬，描述河流一維水流運(yùn)動(dòng)的圣維南方程組包括連續(xù)性方程和動(dòng)量方程：

(1)

(2)

式中：Q為流量；Z為斷面水位；A為過水?dāng)嗝婷娣e；q為旁側(cè)入流量；g為重力加速度；K為流量模數(shù)；x為空間坐標(biāo)；t為時(shí)間坐標(biāo)。

采用廣泛使用的Pressman格式對(duì)方程組進(jìn)行數(shù)值離散求解[10]。該格式采用隱式離散格式，具有較好的穩(wěn)定性，能夠滿足模擬計(jì)算的要求。根據(jù)河道地形資料，對(duì)易澇區(qū)內(nèi)的銅鼓溪、面前溪、太陽河舊河等河道進(jìn)行一維建模。共構(gòu)建概化斷面750個(gè)，概化河段733個(gè)，河道汊點(diǎn)9個(gè)，河段平均間距約70 m。河網(wǎng)模型上游邊界取至太陽河左岸堤防處，考慮到易澇區(qū)內(nèi)河道下游受到小海潮位頂托，下游邊界采用小海潮位邊界條件。

2.2 地表二維水動(dòng)力學(xué)模型

易澇區(qū)內(nèi)地表的洪水演進(jìn)模擬采用二維水動(dòng)力學(xué)模型，可以計(jì)算得到洪水淹沒水深、到達(dá)時(shí)間、流速分布等信息。若不考慮科氏力和風(fēng)力的影響，引入靜水壓力假設(shè)，則二維淺水方程的守恒形式為：

(3)

(4)

式中：U為守恒變量；E、F分別為x和y方向通量；S為源項(xiàng)；S0x、S0y分別為x和y方向的底坡；Sfx、Sfy分別為x和y方向上的摩阻比降；h為水深；g為重力加速度；u、v分別為x和y方向的垂線平均流速。

基于非結(jié)構(gòu)網(wǎng)格采用有限體積法對(duì)二維淺水方程進(jìn)行數(shù)值離散求解[11,12]，采用Roe格式的近似Riemann解來求解界面通量，對(duì)底坡源項(xiàng)特征分解以保證和諧性。根據(jù)易澇區(qū)地形特征、水系分布等資料，確定萬城鎮(zhèn)易澇區(qū)建模范圍面積約59.5 km2。采用非結(jié)構(gòu)網(wǎng)格對(duì)易澇區(qū)區(qū)域進(jìn)行網(wǎng)格剖分，見圖2。網(wǎng)格單元共17 325 個(gè)，網(wǎng)格平均邊長(zhǎng)55 m。為了避免高程點(diǎn)數(shù)據(jù)失準(zhǔn)導(dǎo)致局部高程凹陷造成的淹沒水深的奇異解，對(duì)萬城鎮(zhèn)易澇區(qū)的數(shù)字高程模型進(jìn)行數(shù)據(jù)校核，采用距離加權(quán)法對(duì)網(wǎng)格單元進(jìn)行高程插值。根據(jù)萬城鎮(zhèn)下墊面的土地利用分類情況(水田、草地、林地等)資料，結(jié)合糙率參數(shù)取值表對(duì)網(wǎng)格單元的糙率參數(shù)進(jìn)行賦值。土壤下滲和蒸發(fā)分別按照下滲能力和蒸發(fā)能力在降雨中進(jìn)行扣除后得到凈雨過程，作為暴雨邊界條件，采用直接降雨法進(jìn)行地表二維模型積澇計(jì)算。

圖2 河網(wǎng)概化與計(jì)算網(wǎng)格Fig.2 Generalization of river network and computation grid

2.3 模型耦合

為了實(shí)現(xiàn)河道水流通過河道堤防漫溢進(jìn)入易澇區(qū)及易澇區(qū)內(nèi)洪水回流河道，對(duì)一、二維水動(dòng)力學(xué)模型進(jìn)行耦合計(jì)算[13]，采用下面的寬頂堰流公式計(jì)算堤防漫溢水量實(shí)現(xiàn)一、二維模型水量交換，同時(shí)滿足水量守恒條件。

(5)

h1=max(Zu,Zd)-Ab；h2=min(Zu,Zd)-Zb

(6)

式中：q為通過堤防的單寬流量；Zu、Zd分別為堤防處河道內(nèi)外的水位；Zb為河道堤防高程。

2.4 模型驗(yàn)證

采用2010年10月實(shí)測(cè)暴雨洪水資料對(duì)易澇區(qū)洪水風(fēng)險(xiǎn)分析模型進(jìn)行驗(yàn)證分析。2010年9月30日至10月7日受到臺(tái)風(fēng)影響，萬寧市出現(xiàn)持續(xù)強(qiáng)降雨，20個(gè)觀測(cè)站累計(jì)降雨量均超過800 mm，其中3個(gè)超過1 000 mm，過程降雨為歷史罕見。由于在實(shí)際中很難得到暴雨造成的整個(gè)淹沒區(qū)域內(nèi)所有積澇點(diǎn)的淹沒水深數(shù)據(jù)，考慮到暴雨內(nèi)澇時(shí)的積水一般積蓄在地勢(shì)低洼地帶，所以通常選擇一些代表性低洼處的最大淹沒水深作為模型計(jì)算的驗(yàn)證資料。通過實(shí)地調(diào)研分析收集到5個(gè)地勢(shì)低洼處的代表性積澇點(diǎn)的最大淹沒水深數(shù)據(jù)，調(diào)查點(diǎn)位置見圖3。

圖3 調(diào)查點(diǎn)地理位置圖Fig.3 Geographical position of investigation point

將本次實(shí)測(cè)降雨過程作為降雨輸入邊界條件，采用模型進(jìn)行模擬計(jì)算，將實(shí)際調(diào)查點(diǎn)的計(jì)算最大淹沒水深與實(shí)測(cè)最大淹沒水深進(jìn)行比較驗(yàn)證，見表1。從比較結(jié)果可知，計(jì)算淹沒水深與實(shí)測(cè)淹沒水深基本一致，滿足模型驗(yàn)證要求，說明模型是合理可靠的。

表1 實(shí)測(cè)水深與計(jì)算水深比較 m

3 積澇風(fēng)險(xiǎn)動(dòng)態(tài)分析

3.1 模型封裝與接口調(diào)用

為了實(shí)現(xiàn)對(duì)預(yù)測(cè)降雨、潮位等突發(fā)的洪水風(fēng)險(xiǎn)因子下的易澇區(qū)積澇風(fēng)險(xiǎn)動(dòng)態(tài)分析，將降雨邊界條件和下游潮位邊界條件允許外部訪問對(duì)兩個(gè)邊界條件過程進(jìn)行動(dòng)態(tài)修改，然后調(diào)用模型進(jìn)行洪水淹沒風(fēng)險(xiǎn)計(jì)算。既可以針對(duì)特定水文頻率制作靜態(tài)洪水風(fēng)險(xiǎn)圖，也可以針對(duì)未來預(yù)測(cè)降雨過程、潮位過程進(jìn)行易澇區(qū)積澇風(fēng)險(xiǎn)動(dòng)態(tài)分析。

采用計(jì)算機(jī)程序語言對(duì)計(jì)算模型進(jìn)行封裝處理。設(shè)計(jì)模型與外部系統(tǒng)交互的接口共3個(gè)：①模型輸入接口，用于供外部系統(tǒng)輸入任意的預(yù)測(cè)降雨過程、潮位過程等模型計(jì)算邊界條件，作為模型計(jì)算的動(dòng)態(tài)風(fēng)險(xiǎn)因子；②模型計(jì)算接口，用于供外部系統(tǒng)顯式調(diào)用模型進(jìn)行洪水淹沒方案的模擬計(jì)算；③模型輸出接口，用于模型計(jì)算結(jié)果(淹沒水深、流速等)輸出供外部系統(tǒng)展示洪水演進(jìn)、查詢淹沒水深等。各個(gè)模塊調(diào)用流程圖見圖4。

圖4 模塊調(diào)用流程圖Fig.4 Chart of module transfer flow

3.2 設(shè)計(jì)頻率靜態(tài)洪水風(fēng)險(xiǎn)圖

根據(jù)《廣東省暴雨參數(shù)等值線圖》及《廣東省暴雨徑流查算圖表》中與海南省相關(guān)的水文設(shè)計(jì)分析成果，推求計(jì)算萬城鎮(zhèn)易澇區(qū)的5年一遇、10年一遇、20年一遇、50年一遇24 h設(shè)計(jì)暴雨過程。模型計(jì)算邊界條件中，以設(shè)計(jì)暴雨過程經(jīng)過扣損后的凈雨過程作為二維模型的輸入邊界條件，以小海的多年平均高潮位作為河網(wǎng)一維模型的下游潮位邊界條件。對(duì)5年一遇、10年一遇、20年一遇、50年一遇設(shè)計(jì)頻率下的萬城鎮(zhèn)易澇區(qū)積澇風(fēng)險(xiǎn)進(jìn)行計(jì)算，統(tǒng)計(jì)得到各個(gè)設(shè)計(jì)頻率暴雨下的最大淹沒水深分布圖，如圖5所示。

圖5 不同設(shè)計(jì)降雨下的最大淹沒水深分布圖Fig.5 Maximum submerged water depth distribution map under different design rainfall

4 結(jié) 語

本文構(gòu)建了萬城鎮(zhèn)易澇區(qū)一維、二維耦合水動(dòng)力學(xué)模型，并采用歷史洪水資料進(jìn)行了模型驗(yàn)證，可以有效地模擬暴雨條件下的萬城鎮(zhèn)易澇區(qū)洪水風(fēng)險(xiǎn)。面對(duì)實(shí)際防汛需要，通過降雨、潮位邊界條件的動(dòng)態(tài)設(shè)定構(gòu)建了萬城鎮(zhèn)易澇區(qū)動(dòng)態(tài)洪水風(fēng)險(xiǎn)分析模型?？梢詾榇彐?zhèn)地區(qū)的洪水預(yù)測(cè)預(yù)警、防洪減災(zāi)等提供關(guān)鍵技術(shù)支撐，具有良好的實(shí)際應(yīng)用價(jià)值。