基于洪水淹沒(méi)模擬的水文水動(dòng)力模型研究

李選彧

(遼寧省遼陽(yáng)水文局，遼寧 遼陽(yáng) 111000)

1 水文水動(dòng)力耦合模型

1.1 水文模型

文章對(duì)山區(qū)產(chǎn)匯流的計(jì)算利用BTOPMC模型實(shí)現(xiàn)，該模型是一種分布式物理機(jī)制較強(qiáng)的模擬軟件。TOPMODEL蓄滿產(chǎn)流概化模式和馬斯京根-康奇法是產(chǎn)流過(guò)程、匯流過(guò)程的計(jì)算基礎(chǔ)，在下墊面蒸散發(fā)計(jì)算過(guò)程中土壤水蒸發(fā)能力與冠層截留蒸發(fā)能力，一般利用S-W雙源模型確定。產(chǎn)流、匯流以及地形子模型為組成BTOPMC模型的3個(gè)模塊，該模型主要有馬斯京根-康奇法河道匯流演算參數(shù)和土壤、土地覆蓋參數(shù)兩部分[2]。其中，曼寧糙率系數(shù)n、馬斯京根參數(shù)dt與dl屬于馬斯京根-康奇道匯流演算參數(shù)；而最大根區(qū)蓄水容量Smax、蒸發(fā)函數(shù)參數(shù)α、流量衰減系數(shù)m和飽和土壤水分傳導(dǎo)度T0屬于土壤、土地覆蓋參數(shù)。參數(shù)取值，見(jiàn)表1。

表1 參數(shù)取值

對(duì)于以地形為主要因素影響產(chǎn)流的山區(qū)BTOPMC模型具有較強(qiáng)適用性，水動(dòng)力模型的入流邊界條件選取建成區(qū)與山區(qū)交界處模型計(jì)算的斷面流量時(shí)間序列。

1.2 水動(dòng)力模型

美國(guó)聯(lián)邦應(yīng)管理署(FEMA)認(rèn)證、FLO-2D公司商業(yè)研發(fā)和科羅拉多大學(xué)研發(fā)而成的FLO-2D模型，現(xiàn)已廣泛應(yīng)用于我國(guó)臺(tái)灣地區(qū)，但缺乏其他地區(qū)的案例研究及應(yīng)用[3]。FLO-2D模型普遍應(yīng)用于涵洞與橋梁及不規(guī)則形狀河道、工程風(fēng)險(xiǎn)設(shè)計(jì)、洪水平原管理、城市洪水、泥石流等模擬演算，模型由出流、入流、降雨、水工建筑物、寬度與面積折減、街道、堤防、渠道、洪水平原等多種子模塊組成[4]。其中，建筑物阻礙洪水的做種使用折減系數(shù)表征，比較符合寬度與面積折減模塊的實(shí)際應(yīng)用。

實(shí)際使用過(guò)程中應(yīng)將FLO-2D作如下基本假設(shè)：①演算時(shí)段內(nèi)的流體符合靜水壓力分布，流態(tài)穩(wěn)定；②每個(gè)單元具有相同的糙率及河道形狀，糙率符合曼寧公式且為均勻穩(wěn)定紊流的阻抗。洪水演進(jìn)計(jì)算的理論基礎(chǔ)是簡(jiǎn)單的體積守恒，使水流障礙物能夠在網(wǎng)格系統(tǒng)內(nèi)排出，并以獨(dú)立緩變量流作為正八邊形網(wǎng)格系統(tǒng)單元的各邊水流。采用一維緩變量流的圣維南方程組，構(gòu)造完整動(dòng)力波方程與連續(xù)方程，其表達(dá)式為[5]：

(1)

(2)

式中：Sf、S0為摩擦坡度和河床坡度，m/m，其中Sf值取單位長(zhǎng)度上的總水頭損失；x、V為沿流動(dòng)方向的空間坐標(biāo)和平均流速，m、m/s；h、t、i為水流深度、時(shí)間、超額降雨強(qiáng)度，m、s、mm/h；g為重力加速度，m/s2。

采用中央有限差分顯式法計(jì)算完整動(dòng)力波方程式的數(shù)值，并將其按照運(yùn)算過(guò)程進(jìn)一步分為動(dòng)力波和擴(kuò)散波項(xiàng)目：二階牛頓-拉弗森法的起點(diǎn)選取擴(kuò)散方程式求解的數(shù)值，然后對(duì)二階牛頓-拉弗森法利用完整動(dòng)力波方程式實(shí)行3次迭代運(yùn)算，若3次迭代未能收斂則運(yùn)算結(jié)果取擴(kuò)散波方程式的數(shù)值解[6]。

以獨(dú)立一維水流實(shí)現(xiàn)八邊形網(wǎng)格單元的FLO-2D運(yùn)算，考慮洪水演進(jìn)受建筑物的影響給網(wǎng)格賦值寬度和面積折減系數(shù)。通過(guò)提取網(wǎng)格系統(tǒng)重疊信息與建筑物多邊形圖層，可以確定相應(yīng)的折減系數(shù)，寬度與面積折減系數(shù)反映了建筑物對(duì)水流流動(dòng)的阻礙作用、對(duì)水流存儲(chǔ)的影響程度。

1.3 耦合的模型

采用最簡(jiǎn)單、最有效的松散耦合法，即以模型的輸出與輸入條件轉(zhuǎn)換實(shí)現(xiàn)不同模型的聯(lián)合。該方法保證了模型組分的獨(dú)立完整性且無(wú)需對(duì)模型代碼修改，耦合過(guò)程中不影響各模型功能，但對(duì)于反饋機(jī)制特別是關(guān)鍵物理過(guò)程的描述存在一定不足。將城區(qū)水動(dòng)力模型入流邊界條件設(shè)定為水文模型求解的產(chǎn)區(qū)出水口總流量。耦合結(jié)構(gòu)，見(jiàn)圖1。

圖1 耦合結(jié)構(gòu)

2 實(shí)例應(yīng)用

2.1 流域概況

研究區(qū)為遼寧省遼陽(yáng)市境內(nèi)太子河流域，境內(nèi)河長(zhǎng)142.8km，其面積4000km2，各支流從東、南、北三向匯入構(gòu)成向心水系。太子河葠窩水庫(kù)以上地勢(shì)陡峭、山嶺連綿、森林茂密、植被良好，該河段占全流域60%，將其設(shè)定為山區(qū)；葠窩水庫(kù)至遼陽(yáng)段人口密集、植被覆蓋較少、地勢(shì)較緩，土壤侵蝕嚴(yán)重，將此區(qū)域設(shè)為低山丘陵區(qū)；而遼陽(yáng)以下坡度較緩、河流彎曲、斷面較淺，地貌形態(tài)以平原為主，地勢(shì)低洼，河流泄洪能力差易引起洪澇災(zāi)害，此區(qū)段設(shè)為建成區(qū)。此外，建成區(qū)面積遠(yuǎn)遠(yuǎn)低于山區(qū)匯水面積，并利用水動(dòng)力模型、水模型對(duì)建成區(qū)和山區(qū)建模，由于水位過(guò)高建成區(qū)暴雨洪水時(shí)常出現(xiàn)頂托作用，洪水無(wú)法排出而加劇城市內(nèi)澇。以2019年“利奇馬”號(hào)臺(tái)風(fēng)暴雨為情景，利用水文水動(dòng)力耦合模型模擬城區(qū)洪水淹沒(méi)及上游山區(qū)匯流情況，以期為區(qū)域洪水治理及防洪規(guī)劃提供支持。

2.2 數(shù)據(jù)來(lái)源

土壤分布、土地利用和DEM數(shù)據(jù)為建?；A(chǔ)的GIS數(shù)據(jù)類(lèi)型，其中獲取土壤分布數(shù)據(jù)的流程為：根據(jù)流域邊界條件切割1：100萬(wàn)中國(guó)土壤圖，將土壤圖利用SWAT軟件加載后重分類(lèi)，從而獲取所需的土壤圖。采用以下流程獲取土地利用數(shù)據(jù)：基于易康法、ERDAS軟件和TM7波段數(shù)據(jù)實(shí)現(xiàn)遙感解譯及幾何校正，通過(guò)對(duì)比實(shí)地考察數(shù)據(jù)合理確定土地利用類(lèi)型；然后合并轉(zhuǎn)化用地類(lèi)型圖，將Coverage格式轉(zhuǎn)變?yōu)?shp格式，并結(jié)合邊界條件完成切割，構(gòu)造Landuse查找表文件完成重分類(lèi)，在此基礎(chǔ)上可以輸出土地利用類(lèi)型圖。研究所用DEM數(shù)據(jù)包括1m、30m兩種分辨率數(shù)據(jù)，分別來(lái)源于能實(shí)時(shí)獲取地形表面三維影像與空間信息的航空遙感系統(tǒng)以及中科院計(jì)算機(jī)網(wǎng)絡(luò)信息中心地理空間數(shù)據(jù)云平臺(tái)。以Albers等積圓錐投影所有輸入的數(shù)據(jù)，并且30m分辨率的DEM數(shù)據(jù)與此對(duì)應(yīng)，利用同樣格柵大小的ESRI GRID格式轉(zhuǎn)化土壤空間分布與土地利用數(shù)據(jù)。

將GIMMS NDVI遙感數(shù)據(jù)作為S-W蒸散發(fā)模型數(shù)據(jù)來(lái)源，并采用CRU數(shù)據(jù)集作為月平均氣象數(shù)據(jù)。通過(guò)數(shù)字化處理LIDAR機(jī)載正射影像，可以獲取建筑物圖層以及概化后的建筑物，把生成的1527座建筑物用于后續(xù)寬度與面積折減系數(shù)的賦值。由遼陽(yáng)水文站提供河流水位以及降雨數(shù)據(jù)，并利用壓力式水位計(jì)結(jié)合實(shí)地調(diào)查數(shù)據(jù)獲取研究期間洪泛區(qū)的淹水深度，經(jīng)對(duì)比分析驗(yàn)證模型計(jì)算結(jié)果。

2.3 水文模型率定

考慮到流域內(nèi)水文資料匱乏且水文站點(diǎn)較少的實(shí)際，通過(guò)對(duì)比斷面設(shè)計(jì)流量實(shí)現(xiàn)模型的率定驗(yàn)證。2019年“利奇馬”號(hào)臺(tái)風(fēng)暴雨重現(xiàn)期按照歷史降雨數(shù)據(jù)推算，可以判定為20a一遇，若以20a一遇作為流量重現(xiàn)期，則山區(qū)和低山丘陵區(qū)斷面流量達(dá)到308m3/s、251m3/s。水文模型參數(shù)率定，見(jiàn)表2。

表2 水文模型參數(shù)率定

將2019年“利奇馬”號(hào)臺(tái)風(fēng)暴雨過(guò)程利用率定后的水文模型進(jìn)行模擬，設(shè)定模擬輸出步長(zhǎng)1h。

2.4 水動(dòng)力模型設(shè)置

以建成區(qū)覆蓋面積6.5km2完成水動(dòng)力模型的建模，從模擬精準(zhǔn)度和平衡模型運(yùn)行速度兩個(gè)角度選擇分辨率為10m×10m的網(wǎng)格，由此實(shí)現(xiàn)較高的模擬運(yùn)算速度以及街道、河道等細(xì)節(jié)的全面覆蓋，最終生成52160個(gè)正交網(wǎng)格。選擇2019年8月12日12時(shí)—2019年8月13日18時(shí)為計(jì)算時(shí)段。

研究區(qū)山區(qū)、低山丘陵區(qū)匯流通過(guò)太子河干流輸入系統(tǒng)，由此分為山區(qū)、低山丘陵區(qū)兩部分入流邊界條件。其中，建成區(qū)上游處為低山丘陵區(qū)的匯流點(diǎn)，流域出現(xiàn)暴雨時(shí)將家中城鎮(zhèn)洪澇災(zāi)害和漫堤風(fēng)險(xiǎn)；城鎮(zhèn)洪水淹沒(méi)受上游山區(qū)匯流的影響較弱，由于受到水位頂托作用也添加入流網(wǎng)格點(diǎn)。通過(guò)設(shè)置太子河流域總出流邊界條件約束外河出流水位時(shí)間序列，以保證水文模型與降雨數(shù)據(jù)相同。流量邊界條件，見(jiàn)圖2。

圖2 流量邊界條件

水體流動(dòng)受流體邊界的阻礙作用選用糙率系數(shù)表征，洪水演進(jìn)過(guò)程中的流體、流速取決于糙率的取值大小。流域內(nèi)各下墊面類(lèi)型下糙率n的取值參考FLO-2D手冊(cè)，糙率n的取值，見(jiàn)表3。

表3 糙率n的取值

可透水區(qū)的下滲計(jì)算利用FLO-2D系統(tǒng)自帶的Green-Ampt模塊實(shí)現(xiàn)，參考FLO-2D手冊(cè)及不同的下墊面類(lèi)型確定相關(guān)參數(shù)，選取的參數(shù)如表4?？紤]到城區(qū)土壤性質(zhì)變化不大且面積相對(duì)較小，以分布均勻的沙壤土概化建成區(qū)，其中降雨填洼為區(qū)別下墊面入滲參數(shù)的主要依據(jù)。由于建筑物對(duì)洪水演進(jìn)過(guò)程能夠產(chǎn)生較大的影響，必須概化處理間距可不略補(bǔ)給、密度較大的建筑群，最終獲取的多邊形圖層覆蓋1251座建筑物。寬度與面積折減系數(shù)利用網(wǎng)格重疊區(qū)、多邊形圖層計(jì)算確定，折減系數(shù)賦值的網(wǎng)格數(shù)達(dá)到19255個(gè)。

河道斷面形狀利用高精度DEM數(shù)據(jù)和ARCGIS軟件提取，充分考慮河道沿程變化以及建成區(qū)河道長(zhǎng)度提取18個(gè)斷面數(shù)據(jù)，河底坡度與河床沿程斷面形狀利用插值法及18個(gè)斷面數(shù)據(jù)合理確定。

3 模擬驗(yàn)證

3.1 模型驗(yàn)證

設(shè)模型的平均時(shí)間步長(zhǎng)0.895s、運(yùn)算時(shí)長(zhǎng)0.72h和模擬時(shí)長(zhǎng)20h，水量平衡相對(duì)誤差(5.71×10-6)基本可忽略不計(jì)，絕對(duì)誤差48.05m3。7處新聞圖片點(diǎn)與14處實(shí)地考察點(diǎn)的平均絕對(duì)誤差為0.216m、0.175m，最大絕對(duì)誤差為0.312m、0.450m，水深實(shí)測(cè)值與模擬結(jié)果總體具有較高的吻合度，并且新聞圖片點(diǎn)吻合程度明顯低于實(shí)地考察點(diǎn)。通過(guò)進(jìn)一步分析發(fā)現(xiàn)模型誤差以及實(shí)測(cè)誤差為誤差的主要來(lái)源，其中地形數(shù)據(jù)垂向誤差為引起模型誤差的關(guān)鍵因素，而淹沒(méi)水深信息轉(zhuǎn)換新聞普遍過(guò)程不夠準(zhǔn)確、居民描述存在偏差以及現(xiàn)場(chǎng)勘測(cè)測(cè)量不夠精準(zhǔn)等為實(shí)測(cè)誤差來(lái)源。水動(dòng)力模型驗(yàn)證，見(jiàn)圖3。

圖3 水動(dòng)力模型驗(yàn)證

3.2 洪水淹沒(méi)分析

采用水文水動(dòng)力耦合模型和2019年“利奇馬”號(hào)臺(tái)風(fēng)降雨數(shù)據(jù)，模擬計(jì)算遼陽(yáng)市太子河流域達(dá)最大水深時(shí)間和淹沒(méi)水深。結(jié)果顯示，淹沒(méi)最為嚴(yán)重的地區(qū)位于太子河附近，遼陽(yáng)市東部為淹沒(méi)集中區(qū)域，局部水深達(dá)到0.8m。部分農(nóng)田的淹沒(méi)水深超過(guò)1m，少數(shù)農(nóng)田達(dá)到1.5m，受河道洪水漫堤影響河道兩岸也受到較為嚴(yán)重的淹沒(méi)，河流兩岸20m范圍內(nèi)淹沒(méi)水深達(dá)到1.2m。

根據(jù)達(dá)最大淹水深度所需時(shí)間及其數(shù)據(jù)變化規(guī)律，將各網(wǎng)格數(shù)據(jù)按時(shí)間劃分為5個(gè)區(qū)間，峰現(xiàn)時(shí)間tpeak分布，見(jiàn)表5。由表5可知，基本上峰現(xiàn)時(shí)間≤8.5h的地區(qū)均符合地勢(shì)較高的特征，噶區(qū)域受洪水的影響較弱；峰現(xiàn)時(shí)間在8.5-9.1h的區(qū)域其淹沒(méi)水深基本<0.3m，所占比例達(dá)到最大為39.6%，此外單位小時(shí)降雨量在雨強(qiáng)達(dá)到9h時(shí)最大，該變化特征與峰現(xiàn)時(shí)間保持一致，該區(qū)域淹沒(méi)主要與當(dāng)?shù)亟涤晗嚓P(guān)；太子河流入主城區(qū)后發(fā)生漫堤處附近的峰現(xiàn)時(shí)間界于9.1-11.6h，洪水漫堤流量最大將直接導(dǎo)致其附近出現(xiàn)水深峰值；峰現(xiàn)時(shí)間在11.6-12.5h的區(qū)域其地勢(shì)低洼屬于淹水重災(zāi)區(qū)，淹沒(méi)水深普遍>0.5m，這主要與當(dāng)?shù)亟涤暝谙聣|面產(chǎn)生匯流以及漫堤洪水的演進(jìn)等因素有關(guān)；最后，峰現(xiàn)時(shí)間在12.5-20h之間的低洼區(qū)域其洪水蓄容量達(dá)到最大，淹沒(méi)水深≤1m屬于中等程度受災(zāi)，由于無(wú)法及時(shí)排出洪水致使過(guò)量洪水向四周蔓延。

表5 峰現(xiàn)時(shí)間tpeak分布

4 結(jié) 論

將高運(yùn)算效率的水文模型與具有更為清晰物理機(jī)制的水動(dòng)力模型耦合，以遼陽(yáng)市太子河流域?yàn)槔M城市洪水演進(jìn)過(guò)程。下游城區(qū)水動(dòng)力模型輸入條件設(shè)定為水文模型計(jì)算得到的山區(qū)出口斷面流量，在保證淹沒(méi)模擬精度的同時(shí)實(shí)現(xiàn)洪水模型的快速運(yùn)算。結(jié)果表明，淹沒(méi)水深驗(yàn)證數(shù)據(jù)與耦合模型計(jì)算結(jié)果高度一致，模擬成果可靠且運(yùn)算穩(wěn)定，淹沒(méi)水深平均誤差小可用于洪水預(yù)報(bào)及模擬。由于率定驗(yàn)證的水文數(shù)據(jù)比較缺乏，未來(lái)仍需要進(jìn)一步探討耦合模型優(yōu)化方式，如為了增強(qiáng)模型時(shí)效性、提高運(yùn)算效率采用緊密耦合的方式等。水動(dòng)力模型受地形數(shù)據(jù)的影響顯著，雖然獲取的地形數(shù)據(jù)現(xiàn)階段已達(dá)到較高的精度，還要不斷提升數(shù)據(jù)的準(zhǔn)確性。