，

(長(zhǎng)江科學(xué)院 河流研究所，武漢 430010)

自然潰堤或人工分洪情況下，洪水在洪泛區(qū)內(nèi)的演進(jìn)將給洪泛區(qū)內(nèi)人民生命財(cái)產(chǎn)安全帶來巨大威脅，因此模擬洪水在洪泛區(qū)內(nèi)的演進(jìn)過程，將給防洪決策和防洪搶險(xiǎn)提供技術(shù)支撐。目前，許多學(xué)者對(duì)分蓄洪區(qū)洪水演進(jìn)開展了數(shù)值模擬研究[1-3]，然而分蓄洪區(qū)地形獲取的不便和模擬成果顯示不直觀等問題，仍然是制約數(shù)學(xué)模型發(fā)展的重要因素。近10多年來，地理信息技術(shù)得到快速發(fā)展，形成了許多的地理信息平臺(tái)。因此，如何將洪水演進(jìn)模擬結(jié)果顯示在地理信息平臺(tái)之上，是一件非常有意義的事情。

Google Earth軟件(以下簡(jiǎn)稱GE)便是其中應(yīng)用最為廣泛的通用的地理信息平臺(tái)之一，它把衛(wèi)星照片或航空照相和GIS信息布置在一個(gè)三維的地球模型上，而且免費(fèi)對(duì)大眾開放。本文基于Google Earth軟件平臺(tái)，建立了分蓄洪區(qū)水沙演進(jìn)數(shù)學(xué)模型，能快速完成模型數(shù)據(jù)前處理和模擬計(jì)算，并實(shí)現(xiàn)了計(jì)算結(jié)果在三維虛擬場(chǎng)景中的動(dòng)態(tài)演示。該模型系統(tǒng)適用于緊急條件下的快速計(jì)算和直觀展示，計(jì)算成果可為分蓄洪區(qū)防洪決策與防洪搶險(xiǎn)提供技術(shù)參考。

1 基于非結(jié)構(gòu)網(wǎng)格的數(shù)值求解方法

分蓄洪區(qū)一般不規(guī)則，為此采用非結(jié)構(gòu)網(wǎng)格[4]進(jìn)行分蓄洪區(qū)的網(wǎng)格劃分是非常合適的。非結(jié)構(gòu)網(wǎng)格可以是三角形網(wǎng)格、四邊形網(wǎng)格或三角形與四邊形混合網(wǎng)格。

模型基本方程采用守恒形式的水深平均平面二維水流基本方程[5]。采用基于同位網(wǎng)格的有限體積算法進(jìn)行方程的離散[6]，流速矢量基于直角坐標(biāo)系，節(jié)點(diǎn)布置采用單元中心法；采用SIMPLEC算法進(jìn)行離散方程的求解，該算法的優(yōu)點(diǎn)在于即便采用較粗的網(wǎng)格也能保持水量和動(dòng)量的守恒。

2 基于Google Earth的模擬演示系統(tǒng)

本文基于Google Earth軟件進(jìn)行二次開發(fā)，將數(shù)值模擬成果在由航拍地貌和三維建筑物構(gòu)成的虛擬地球上進(jìn)行顯示，有助于更好地分析計(jì)算成果的合理性，同時(shí)還可利用Google Earth軟件縮放、平移、鳥瞰等功能，更直觀地查看計(jì)算成果。

Google Earth軟件的二次開發(fā)有2種方式，一種是基于ComAPI；一種是基于KML(全稱為Keyhole Markup Language)?；?ComAPI的開發(fā)方式主要用來控制GE的視角,實(shí)現(xiàn)動(dòng)畫效果，而基于KML的開發(fā)方式主要用來生成地理要素(如點(diǎn)、線、面、多邊形、多面體以及模型等)，實(shí)現(xiàn)數(shù)據(jù)的動(dòng)態(tài)更新等。在實(shí)際應(yīng)用中，需要結(jié)合這2種開發(fā)方式，才能得到比較理想的效果。

筆者采用VB和Fortran語(yǔ)言混合編程方法，對(duì)Google Earth軟件進(jìn)行二次開發(fā)，研制了“基于虛擬地球的水沙水質(zhì)數(shù)值模擬與演示系統(tǒng)V2.0”，系統(tǒng)由前處理、模型庫(kù)、后處理及數(shù)據(jù)庫(kù)等部分構(gòu)成，模型可實(shí)現(xiàn)干河床的水沙演進(jìn)計(jì)算。該系統(tǒng)不僅實(shí)現(xiàn)了界面操作的可視化，而且實(shí)現(xiàn)了計(jì)算結(jié)果在虛擬地球上的動(dòng)態(tài)演示，系統(tǒng)主界面見圖1。

圖1 可視化界面圖

2.1 基于ComAPI的開發(fā)

GE 公布的API目前是1.0a版，數(shù)量比較少，所能實(shí)現(xiàn)的功能也不多，在實(shí)際開發(fā)中主要用來加載KML數(shù)據(jù)、控制視角和實(shí)現(xiàn)動(dòng)畫播放。加載KML數(shù)據(jù)命令有：

(1) IApplicationGE. OpenKmlFile(kmlFilePath)，加載指定的KML文件到GE中；

(2) IApplicationGE.LoadKmlData(kmlData)，加載KML文本流到GE中。

2.2 基于KML文件的開發(fā)

2.2.1 坐標(biāo)轉(zhuǎn)換

目前我國(guó)的河道地形圖大地基準(zhǔn)面一般采用54坐標(biāo)系或西安80坐標(biāo)系，而Google Earth軟件是基于WGS84坐標(biāo)系統(tǒng)，因此在生成KML文件之前需對(duì)計(jì)算結(jié)果文件進(jìn)行坐標(biāo)轉(zhuǎn)換。上述3類坐標(biāo)系的轉(zhuǎn)換在相關(guān)教科書中有詳細(xì)的說明，這里不再贅述。

2.2.2 流速矢量場(chǎng)

計(jì)算成果顯示包括流速矢量場(chǎng)和水位、水深、流速等物理量的標(biāo)量場(chǎng)，同時(shí)還需文字標(biāo)注、圖例。為進(jìn)行流場(chǎng)的動(dòng)態(tài)演示，需進(jìn)行每幀流場(chǎng)的時(shí)間控制。

圖2 流速箭頭示意圖

計(jì)算域內(nèi)每個(gè)點(diǎn)的流速矢量可用箭頭表示，箭頭長(zhǎng)短表示流速大小，箭頭方向代表流速方向。每個(gè)箭頭可由6個(gè)有序點(diǎn)形成的折線表示，如圖2所示。

在KML文件中，每個(gè)點(diǎn)的流速矢量可用LineString標(biāo)簽來定義，代碼如下：

112.138 084 0, 30.284 331 1, 0

112.140 711 6, 30.283 183 1, 0

......

其中字段表示定義標(biāo)簽；定義該箭頭名稱；定義流速箭頭顏色，需根據(jù)流速的大小預(yù)先定義；繪制折線；下面定義6個(gè)點(diǎn)的坐標(biāo)，注意這里采用的是經(jīng)轉(zhuǎn)換后的經(jīng)緯度坐標(biāo)。

2.2.3 等值線云圖

等值線的繪制與流速箭頭的繪制類似，只需繪制每條折線就行，云圖的繪制需將封閉折線構(gòu)成的區(qū)域進(jìn)行顏色填充，可用Polygon標(biāo)簽來定義，代碼格式如下：

112.141 431 5, 30.283 587 2, 0

112.140 788 8, 30.281 146 2, 0

112.144 288 4, 30.282 297 2, 0

112.141 431 5, 30.283 587 2, 0

112.141 431 5, 30.283 587 2, 0

2.2.4 文字標(biāo)注

流速值、水深值或其他文本的標(biāo)注可采用Google Earth中的Point地標(biāo)標(biāo)簽來定義，代碼如下：

2.2.5 流場(chǎng)圖例

完整的流場(chǎng)顯示需定義圖例，在Google Earth中圖例由ScreenOverlay標(biāo)簽控制，在動(dòng)態(tài)演示過程中保持不變，代碼如下：

2.2.6 時(shí)間控制

為實(shí)現(xiàn)流場(chǎng)的動(dòng)態(tài)演示，需定義每幀流場(chǎng)的起止時(shí)間，采用TimeSpan標(biāo)簽定義，代碼如下：

3 在荊江分蓄洪區(qū)的應(yīng)用

荊江分蓄洪區(qū)位于長(zhǎng)江南岸，是長(zhǎng)江中游防洪系統(tǒng)的主要組成部分，分洪區(qū)面積921 km2，有效蓄洪容量54億m3。荊江分洪工程包括進(jìn)洪閘、節(jié)制閘和分洪區(qū)圍堤。進(jìn)洪閘位于太平口東岸，設(shè)計(jì)進(jìn)洪流量8 000 m3/s。荊江分蓄洪區(qū)的功能在于分泄荊江河段的超額洪水，使洪水能夠安全通過，以保護(hù)荊江兩岸廣大農(nóng)田和人民生命財(cái)產(chǎn)安全。

1954年7月22日至8月1日荊江分洪區(qū)進(jìn)行了3次分洪運(yùn)用，分洪總量為122.6億m3。采用建立的模型系統(tǒng)，對(duì)1954年荊江分蓄洪區(qū)分洪過程進(jìn)行了模擬計(jì)算。

3.1 計(jì)算條件

計(jì)算范圍為整個(gè)荊江分蓄洪區(qū)，以北閘為進(jìn)流邊界，南閘為出流邊界。網(wǎng)格采用非結(jié)構(gòu)網(wǎng)格，見圖3。計(jì)算網(wǎng)格地形根據(jù)航拍的DEM地形數(shù)據(jù)插值得到。

圖3 分洪區(qū)計(jì)算網(wǎng)格布置

將荊江分蓄洪區(qū)分洪過程概化為潰壩過程，潰壩計(jì)算關(guān)鍵在于如何進(jìn)行干河床與動(dòng)邊界的處理，筆者提出了“干濕雙最小水深法[6]”，能較好反映潰壩水流在干河床上的演進(jìn)過程。

由于荊江分蓄洪區(qū)內(nèi)缺乏實(shí)測(cè)的洪水資料，因此，糙率參考相關(guān)資料確定[6]：樹林0.070，旱地0.065，水田0.050，水面0.025。如果某網(wǎng)格內(nèi)含有多種地形，則按照各種地形糙率的加權(quán)平均值確定該網(wǎng)格的糙率。

模擬計(jì)算從1954年7月22日8:00開始，考慮為全潰，直至蓄滿為止。進(jìn)口含沙量假定為零，床沙中值粒徑取0.2 mm。

3.2 計(jì)算成果分析

3.2.1 洪水演進(jìn)計(jì)算成果

計(jì)算結(jié)果表明，分洪閘下泄流量先緩慢增大，最大可達(dá)8 000 m3/s，后快速減小，直至蓄滿所需時(shí)間為11 d左右。

圖4給出不同時(shí)刻洪水在分蓄洪區(qū)的演進(jìn)情況(圖中箭頭表示流速矢量大小和方向，藍(lán)色云圖表示水深，藍(lán)色越渾則水深越大)。從圖4中可見，系統(tǒng)較直觀地反映了洪水在航拍地貌上的演進(jìn)過程，在潰口和分洪區(qū)縮窄段，流速明顯增大，最大流速為1.2 m/s，圖中彩色云圖能直觀反映不同時(shí)刻洪水淹沒范圍及淹沒水深。

圖4 不同時(shí)刻洪水演進(jìn)流場(chǎng)圖

3.2.2 泥沙輸移計(jì)算成果

圖5(a)給出了分洪67 h時(shí)分蓄洪區(qū)含沙量分布情況(圖中顏色越紅表示含沙量越大)。從圖4中可見，由于進(jìn)口未給定含沙量，泥沙主要來自于床面的沖刷，因此流速越大部位(如圖5(a)中深藍(lán)色區(qū)域)，含沙量也會(huì)越大，最大含沙量可達(dá)2 kg/m3左右，位于分洪口門下游附近；而在分洪區(qū)下游，由于流速不斷減小，含沙量也隨之減小。

圖5 分蓄洪區(qū)含沙量和沖淤厚度分布(T=67 h)

圖5(b)給出了分洪67 h后分蓄洪區(qū)沖淤厚度分布情況(圖中紅色表示淤積，藍(lán)色表示沖刷，顏色越深則沖淤幅度越大)。從圖中可見，一般大流速部位，床面易發(fā)生沖刷，最大沖深為1.7 m，其余大部分區(qū)域發(fā)生淤積，最大淤厚約1.5 m。

4 結(jié) 論

(1) 建立了基于地理信息技術(shù)的分蓄洪區(qū)平面二維水沙演進(jìn)數(shù)學(xué)模型。模型能直接在Google Earth軟件環(huán)境中確定計(jì)算域邊界，采用無結(jié)構(gòu)的三角形、四邊形網(wǎng)格進(jìn)行計(jì)算域的網(wǎng)格劃分，根據(jù)航拍的DEM地形插值得到計(jì)算網(wǎng)格地形高程，提高了前處理的效率。

(2) 對(duì)Google Earth軟件進(jìn)行二次開發(fā)，研制了計(jì)算結(jié)果后處理模塊，實(shí)現(xiàn)了分蓄洪區(qū)水沙演進(jìn)過程在三維航拍地貌中的動(dòng)態(tài)演示，并能與Google Earth同步實(shí)現(xiàn)縮放、鳥瞰、漫游等多種查看功能。

(3) 該模型系統(tǒng)在長(zhǎng)江荊江分洪區(qū)得到成功運(yùn)用，復(fù)演了1954年分洪后的水沙演進(jìn)過程，并在Google Earth地貌中進(jìn)行動(dòng)態(tài)演示，可直觀查看水深、流速、洪水淹沒范圍、沖深等要素，可為洪水風(fēng)險(xiǎn)評(píng)估和應(yīng)急預(yù)案的制定提供參考。

參考文獻(xiàn)：

[1] MACCHIONE F, MORELLI M. Practical Aspects in Comparing Shock-Capturing Schemes for Dam-break Problems [J]. Journal of Hydraulic Engineering, 2002, 129(3)：187-195.

[2] 余明輝，張小峰. 平面二維潰堤水流泥沙數(shù)值模擬[J]. 水科學(xué)進(jìn)展，2001， 12(3)：286-290.(YU Ming-hui, ZHANG Xiao-feng. Horizontal 2-D Uneven Sediment Mathematical Model in Dike Burst[J]. Advances in Water Science, 2001, 12(3): 286-290.(in Chinese))

[3] 王志剛，汪德爟，賴錫軍, 等. 下游為干河床的潰壩水流數(shù)值模擬[J]. 水利水運(yùn)工程學(xué)報(bào)，2003，(2):18-23.(WANG Zhi-gang, WANG De-huan, LAI Xi-jun,etal. Numerical Simulation of Flow on Dry Bed for Dam Break[J]. Hydro-science and Engineering, 2003, (2):18-23.(in Chinese))

[4] 張細(xì)兵.河道有限元網(wǎng)格自動(dòng)剖分方法研究[J]. 長(zhǎng)江科學(xué)院院報(bào)，2002，19(3):19-21.(ZHANG Xi-bing. Auto-generating Method of Triangulation Grids for River Reach[J]. Journal of Yangtze River Scientific Research Institute, 2002,19(3):19-21.(in Chinese))

[5] 張細(xì)兵，范北林. 潰壩洪水演進(jìn)平面二維數(shù)學(xué)模型初步研究[J]. 長(zhǎng)江科學(xué)院院報(bào), 2006，23(6):19-21.(ZHANG Xi-bing, FAN Bei-lin. Numerical Simulation for Two-Dimensional Dam-break Flood Flow[J]. Journal of Yangtze River Scientific Research Institute, 2006, 23(6): 19-21.(in Chinese))

[6] 張細(xì)兵，歐治華，崔占峰，等. 基于非結(jié)構(gòu)網(wǎng)格的分蓄洪區(qū)水沙演進(jìn)數(shù)學(xué)模型研究[J].長(zhǎng)江科學(xué)院院報(bào)，2011，28(4):75-79.(ZHANG Xi-bing, OU Zhi-hua, CUI Zhan-feng,etal. Unstructured Grid Model of Flow and Sediment Evolution in Flood Diversion and Storage Area[J]. Journal of Yangtze River Scientific Research Institute, 2011, 28(4):75-79.(in Chinese))