珠江口風(fēng)暴潮數(shù)值模擬系統(tǒng)的設(shè)計(jì)與實(shí)現(xiàn)

葉榮輝，宋志堯，沈 正，孔 俊，張 卓，金光球

(1.珠江水利委員會(huì)珠江水利科學(xué)研究院，廣東廣州 510611;2.河海大學(xué)水文水資源及水利工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，江蘇南京 210098;3.南京師范大學(xué)虛擬地理環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室，江蘇南京 210046)

珠江水系經(jīng)八大口門注入南海，河口地區(qū)水流受徑流和潮流的共同作用，在臺(tái)風(fēng)的侵襲下，珠江口各口門均會(huì)發(fā)生程度不同的臺(tái)風(fēng)風(fēng)暴潮災(zāi)害，水動(dòng)力環(huán)境較為復(fù)雜，因此建立一套能夠準(zhǔn)確反映珠江口復(fù)雜水動(dòng)力特征的數(shù)值模擬系統(tǒng)對于臺(tái)風(fēng)風(fēng)暴潮的預(yù)警預(yù)報(bào)、河口綜合治理規(guī)劃的制定具有現(xiàn)實(shí)意義。目前經(jīng)典的河口、海岸數(shù)學(xué)模型如 POM［1］，ECOM［2］，F(xiàn)VCOM［3］，ELCIRC［4］等大都采用 Fortran 語言開發(fā)，沒有提供人機(jī)交互界面，模型參數(shù)調(diào)試、邊界條件設(shè)置過程繁瑣，特別是對于存在多條開邊界的復(fù)雜研究區(qū)域的風(fēng)暴潮數(shù)值模擬，模型參數(shù)更多，邊界條件給定更為復(fù)雜。當(dāng)計(jì)算時(shí)段、模型時(shí)間步長等發(fā)生改變時(shí)，對邊界條件的手動(dòng)設(shè)置需要耗費(fèi)大量時(shí)間。同時(shí)，這些模型的計(jì)算結(jié)果為一系列文本文件，不清晰直觀，對結(jié)果的后處理通常需要借助于第三方軟件如 Excel，Surfer，Tecplot等［5-7］，而在運(yùn)用 Surfer，Tecplot等軟件對流場進(jìn)行可視化操作時(shí)，還需要熟悉這些商業(yè)軟件復(fù)雜的數(shù)據(jù)格式，工作量大、效率低。

本文針對珠江口地區(qū)風(fēng)暴潮數(shù)值模擬過程復(fù)雜繁瑣的特點(diǎn)，運(yùn)用混合編程技術(shù)，充分利用現(xiàn)有數(shù)學(xué)模型代碼，結(jié)合Fortran語言較強(qiáng)的數(shù)值計(jì)算能力及地理信息系統(tǒng)(GIS)在可視化方面的優(yōu)勢，設(shè)計(jì)并建立了珠江口風(fēng)暴潮數(shù)值模擬系統(tǒng)。經(jīng)過對200814號臺(tái)風(fēng)風(fēng)暴潮的實(shí)際模擬應(yīng)用，結(jié)果表明，該系統(tǒng)界面友好、性能穩(wěn)定、模擬精度較高。

1 系統(tǒng)設(shè)計(jì)

1.1 功能結(jié)構(gòu)

珠江口風(fēng)暴潮數(shù)值模擬系統(tǒng)分為四大模塊。(a)地圖操作模塊:實(shí)現(xiàn)對地圖的放大、縮小、漫游、圖層控制等功能。(b)臺(tái)風(fēng)信息模塊:提供臺(tái)風(fēng)路徑可視化、臺(tái)風(fēng)屬性信息動(dòng)態(tài)查詢等功能。(c)風(fēng)暴潮模擬計(jì)算模塊:集成臺(tái)風(fēng)模型、珠江口風(fēng)暴潮數(shù)學(xué)模型(簡稱珠江口模型)、南中國海風(fēng)暴潮數(shù)學(xué)模型(簡稱南海模型)等專業(yè)數(shù)學(xué)模型，提供人機(jī)交互界面，實(shí)現(xiàn)模型參數(shù)設(shè)置、邊界條件設(shè)置、風(fēng)暴潮數(shù)值模擬等功能。(d)結(jié)果展示模塊:提供模型計(jì)算結(jié)果的后處理功能，主要用于實(shí)現(xiàn)重要站點(diǎn)潮位過程線、流場的繪制。

1.2 層次結(jié)構(gòu)

珠江口風(fēng)暴潮數(shù)值模擬系統(tǒng)主要分為3個(gè)層次，如圖1所示。(a)數(shù)據(jù)資源層:由空間及非空間數(shù)據(jù)庫組成，主要為該系統(tǒng)提供數(shù)據(jù)支持。空間數(shù)據(jù)庫主要包含基礎(chǔ)地理信息數(shù)據(jù)及屬性信息;非空間數(shù)據(jù)庫包含臺(tái)風(fēng)數(shù)據(jù)庫、水文數(shù)據(jù)庫、模型方法庫等。(b)專業(yè)模型層:在整合臺(tái)風(fēng)、水文等數(shù)據(jù)的基礎(chǔ)上，依次調(diào)用臺(tái)風(fēng)模型計(jì)算臺(tái)風(fēng)風(fēng)場，為風(fēng)暴潮數(shù)值模擬提供驅(qū)動(dòng)風(fēng)場;調(diào)用南海模型為珠江口風(fēng)暴潮數(shù)值計(jì)算提供開邊界;最后調(diào)用珠江口模型進(jìn)行水力計(jì)算。(c)服務(wù)層:為用戶提供模型參數(shù)設(shè)置、臺(tái)風(fēng)信息查詢以及計(jì)算結(jié)果可視化等功能。

2 系統(tǒng)實(shí)現(xiàn)

2.1 開發(fā)方式

珠江口風(fēng)暴潮數(shù)值模擬系統(tǒng)前臺(tái)集成開發(fā)環(huán)境選擇Visual Studio 2008，開發(fā)語言 C#，數(shù)據(jù)庫平臺(tái) MS SQL Server 2008，各專業(yè)數(shù)學(xué)模型均采用Fortran語言開發(fā)，編譯器采用Visual Fortran 6.5。系統(tǒng)集成基于組件式GIS技術(shù)，采用混合編程思路實(shí)現(xiàn)，GIS二次開發(fā)組件選用ArcGISEngine 9.3。

2.2 模型開發(fā)

珠江口風(fēng)暴潮數(shù)值模擬系統(tǒng)中風(fēng)暴潮模擬計(jì)算模塊是核心。由于珠江口地區(qū)河網(wǎng)密布，水道縱橫，分八大口門入海，因此模式計(jì)算范圍覆蓋整個(gè)珠江三角洲河網(wǎng)區(qū)及口外海區(qū)［8］(圖2)。模型上游邊界分別取在石咀、高要、石角、博羅和老鴉崗;模型的下邊界取在珠江口外南海30 m等深線附近。模型采用非結(jié)構(gòu)網(wǎng)格，網(wǎng)格數(shù)為94334，網(wǎng)格節(jié)點(diǎn)為92221，網(wǎng)格最小步長為20 m。

2.2.1 模型理論

珠江口模型基于ELCIRC模型建立，控制方程如下:

圖1 系統(tǒng)層次結(jié)構(gòu)Fig.1 System framework structure

其中

式中:u，v——x，y 向流速;t——時(shí)間;η——水位;D——全水深;f——柯氏力因子;g——重力加速度;ρ0——水的密度;τwx，τwy——x，y 向 10 m 高空處的風(fēng)應(yīng)力;ρa(bǔ)——空氣密度，CDs——風(fēng)應(yīng)力系數(shù)，wx，wy——x，y 向 10 m 高空處的風(fēng)速，采用 Myers公式［9］計(jì)算獲得;τbx，τby——x，y 向底部摩擦應(yīng)力;n——糙率系數(shù)。

2.2.2 定解條件

定解條件包括初始條件和邊界條件:潮位和流速初始值均設(shè)為零。閉邊界滿足流體不可入條件。上游開邊界給定流量過程，由系統(tǒng)從水文數(shù)據(jù)庫中讀取歷史數(shù)據(jù)獲得;外海開邊界給定潮位過程，由系統(tǒng)調(diào)用南海模型［10］計(jì)算獲得。

2.3 混合編程

由于各專業(yè)數(shù)學(xué)模型均采用Fortran語言開發(fā)，而系統(tǒng)界面采用C#開發(fā)。為充分利用兩種語言各自的特點(diǎn)，將各專業(yè)數(shù)學(xué)模型編譯為動(dòng)態(tài)鏈接庫(DLL)，以C#為主程序來調(diào)用風(fēng)暴潮計(jì)算內(nèi)核程序。

用混合語言編程存在由于各自語言特點(diǎn)不同而相互不匹配的問題，在C#與Fortran語言混合編程的過程中，應(yīng)該注意不同語言之間的命名約定、調(diào)用約定以及參數(shù)傳遞約定三方面的問題［11-12］，解決這些問題是實(shí)現(xiàn)C#與Fortran語言混合編程的關(guān)鍵。

圖2 模型計(jì)算范圍Fig.2 Model’s calculation area

2.4 可視化功能

珠江口風(fēng)暴潮數(shù)值模擬系統(tǒng)采用組件式GIS技術(shù)［13-14］實(shí)現(xiàn)臺(tái)風(fēng)路徑、潮位過程線、流場等的可視化功能，對數(shù)據(jù)庫的訪問采用ADO.NET數(shù)據(jù)訪問技術(shù)，空間信息與屬性數(shù)據(jù)通過ID號進(jìn)行關(guān)聯(lián)。

臺(tái)風(fēng)路徑數(shù)據(jù)存儲(chǔ)在臺(tái)風(fēng)數(shù)據(jù)庫中，獲取數(shù)據(jù)以后，先生成路徑定位點(diǎn)并賦上空間坐標(biāo)等屬性信息，而后根據(jù)每個(gè)定位點(diǎn)的坐標(biāo)生成臺(tái)風(fēng)路徑。

潮位數(shù)據(jù)包括歷史潮位數(shù)據(jù)及模擬潮位數(shù)據(jù)。歷史潮位數(shù)據(jù)存儲(chǔ)在數(shù)據(jù)庫中，而模擬潮位數(shù)據(jù)以文本形式存儲(chǔ)，系統(tǒng)直接從文本中讀取。潮位過程線的繪制通過Microsoft Chart Control 6.0(SP4)控件完成。

流場數(shù)據(jù)是矢量數(shù)據(jù)，同時(shí)具有數(shù)值和方向的特征，不能直接可視化［15］。珠江口風(fēng)暴潮數(shù)值模擬系統(tǒng)首先從數(shù)值計(jì)算結(jié)果文本文件中讀取流場信息數(shù)據(jù)，包括網(wǎng)格空間點(diǎn)坐標(biāo)、流速、流向等，然后進(jìn)行矢量數(shù)據(jù)的映射及繪制和顯示。矢量數(shù)據(jù)的映射采用幾何形狀的點(diǎn)圖標(biāo)法，用矢量箭頭表示水流的方向，用其長度表示流速的大小。矢量箭頭的起始點(diǎn)A的坐標(biāo)為該網(wǎng)格中心點(diǎn)坐標(biāo)，其余3個(gè)頂點(diǎn)B，C，D坐標(biāo)可參考A點(diǎn)坐標(biāo)通過簡單的幾何計(jì)算獲得(圖3)。在程序中矢量箭頭被定義為一個(gè)Vector類，在讀取數(shù)據(jù)的時(shí)候，把所有的數(shù)據(jù)以一個(gè)個(gè)Vector實(shí)例的形式進(jìn)行編輯并存儲(chǔ)到泛型集合List＜Vector＞Varr中。最后通過ArcEngine中幾何對象以要素的形式進(jìn)行可視化表達(dá)，并將同一類型的幾何對象以圖形圖層方式進(jìn)行管理。

圖3 箭頭結(jié)構(gòu)示意Fig.3 Schematic diagram of arrow structure

3 系統(tǒng)應(yīng)用

200814 號強(qiáng)臺(tái)風(fēng)引發(fā)了珠江口特大風(fēng)暴潮，珠江口地區(qū)燈籠山、黃埔和南沙站等潮位均超過歷史極值。本文以該次強(qiáng)臺(tái)風(fēng)為例，運(yùn)用珠江口風(fēng)暴潮數(shù)值模擬系統(tǒng)對該次臺(tái)風(fēng)風(fēng)暴潮進(jìn)行模擬。

3.1 臺(tái)風(fēng)信息查詢

在界面中選擇需要查看的年份，系統(tǒng)會(huì)根據(jù)所選年份查找符合條件的數(shù)據(jù)，并以列表的形式顯示。圖4為200814號臺(tái)風(fēng)路徑，從圖4可以看出，該臺(tái)風(fēng)屬西進(jìn)型臺(tái)風(fēng)，珠江口地區(qū)位于臺(tái)風(fēng)路徑前進(jìn)方向的右半圓區(qū)為風(fēng)暴增水區(qū)。通過鼠標(biāo)點(diǎn)擊該場臺(tái)風(fēng)路徑上離珠江口直線距離較近的點(diǎn)，獲取這些時(shí)刻臺(tái)風(fēng)強(qiáng)度屬性信息發(fā)現(xiàn):在離珠江口地區(qū)較近的區(qū)域，該臺(tái)風(fēng)中心氣壓保持在935 hPa左右，中心最大風(fēng)速達(dá)到了50 m/s。結(jié)合其路徑，可以從宏觀上判斷該臺(tái)風(fēng)可引起珠江口地區(qū)較大的風(fēng)暴潮增水。

3.2 風(fēng)暴潮數(shù)值模擬

在進(jìn)行珠江口風(fēng)暴潮數(shù)值模擬前，需要對模型參數(shù)進(jìn)行設(shè)置，包括起止時(shí)間、時(shí)間步長、風(fēng)應(yīng)力系數(shù)、糙率等。在本次應(yīng)用中，通過模型參數(shù)設(shè)置界面，模型計(jì)算時(shí)間段設(shè)為2008年9月22日8:00至2008年9月24日8:00，時(shí)間步長設(shè)為180s，CDs取為0.0026，n設(shè)定為0.023。

圖4 臺(tái)風(fēng)路徑Fig.4 Tyhoon track

本文建立的珠江口模型包含5條上游開邊界以及3條外海開邊界，若采用手動(dòng)設(shè)置邊界條件的方法，設(shè)置過程繁瑣，工作量大。而在珠江口風(fēng)暴潮數(shù)值模擬系統(tǒng)中，通過界面上的邊界計(jì)算按鈕，即可根據(jù)計(jì)算時(shí)段、時(shí)間步長等相關(guān)設(shè)置自動(dòng)完成指定時(shí)間段的上邊界、外海開邊界文件的生成。風(fēng)暴潮模擬計(jì)算所需的風(fēng)壓場文件通過調(diào)用臺(tái)風(fēng)模型生成。在完成對珠江口模型運(yùn)行所需的文本設(shè)置后，運(yùn)行該模型，模型將自動(dòng)調(diào)用地形文件、邊界文本文件、風(fēng)壓場文本文件等進(jìn)行模擬計(jì)算。數(shù)值計(jì)算結(jié)果通過GIS平臺(tái)統(tǒng)一展現(xiàn)。數(shù)值模擬過程如圖5所示。

3.3 結(jié)果可視化

模擬結(jié)果輸出時(shí)段為2008年9月23日8:00至2008年9月24日8:00，圖6為利用珠江口風(fēng)暴潮數(shù)值模擬系統(tǒng)計(jì)算得出的200814號臺(tái)風(fēng)期間燈籠山站潮位過程線計(jì)算值。從圖6可以看出，前期由于距離臺(tái)風(fēng)中心較遠(yuǎn)，燈籠山站未發(fā)生明顯風(fēng)暴增水，從9月23日16:00開始出現(xiàn)明顯的風(fēng)暴增水，并于9月24日3:30達(dá)到整個(gè)風(fēng)暴潮的最高潮位2.66 m。而實(shí)測此次風(fēng)暴潮的最高潮位為2.73 m，出現(xiàn)在9月24日3:00［16］。與實(shí)測值相比，模型計(jì)算得出的最高潮位誤差僅為0.07 m，相位差為30 min，表明該系統(tǒng)模型具有較高的精度。

圖7為相應(yīng)計(jì)算時(shí)段內(nèi)某時(shí)刻全局流場圖，用戶可以通過對圖層的縮放、平移等功能查看感興趣區(qū)域的流場等內(nèi)容。

圖5 風(fēng)暴潮數(shù)值模擬過程Fig.5 Storm surge numerical simulation process

4 結(jié) 論

針對珠江河口地區(qū)風(fēng)暴潮數(shù)值模擬過程復(fù)雜繁瑣的特點(diǎn)，運(yùn)用混合編程技術(shù)，結(jié)合Fortran語言較強(qiáng)的數(shù)值計(jì)算能力及GIS在可視化方面的優(yōu)勢，設(shè)計(jì)并建立了珠江口風(fēng)暴潮數(shù)值模擬系統(tǒng)。該系統(tǒng)集成了前處理、計(jì)算模擬和后處理等功能，可實(shí)現(xiàn)臺(tái)風(fēng)路徑動(dòng)態(tài)顯示、潮位過程線繪制及流場可視化，方便用戶對計(jì)算結(jié)果的分析。通過對200814號臺(tái)風(fēng)風(fēng)暴潮的實(shí)際模擬測試，表明該系統(tǒng)的建立能有效提高工作效率，用戶可借助該系統(tǒng)從宏觀與微觀尺度上分析、評價(jià)臺(tái)風(fēng)對珠江口地區(qū)水動(dòng)力環(huán)境的影響。

圖6 潮位過程線Fig.6 Time series of tide

圖7 流場可視化Fig.7 Visualization of current field

［1］BLUMBERG A F，MELLOR G L.A description of a three-dimensional coastal ocean circulation model［M］//HEAPS N.Three-Dimentional Coastal Ocean Models.Washington:American Geophys Union，1987.

［2］BLUMBERG A F，GOODRICH D M.A simulation of wind-induced destratification in Chesapeake Bay［J］.Estuaries，1990，13:1236-1249.

［3］CHEN C，LIU H，BEARDSLEY R C.An unstructured，finite-volume，three-dimensional，primitive e-quation ocean model:application to coastal ocean and estuaries［J］.Journal of Atmospheric and Oceanic Technology，2003，20:159-186.

［4］ZHANG Y，BAPTISTA A M.A cross-scale model for 3D baroclinic circulation in estuary-plume-shelf systems:I.formulation and skill assessment［J］.Continental Shelf Research，2005，25:935-972.

［5］陳可鋒，喻國華，徐永輝.水流數(shù)值模擬中的混合編程技術(shù)［J］.水運(yùn)工程，2007(2):14-18.(CHEN Kefeng，YU Guohua，XU Yonghui.Mixed-language programming of flow numerical simulation［J］.Port＆ Waterway Engineerin，2007(2):14-18.(in Chinese))

［6］胡艷進(jìn).平面二維水流數(shù)學(xué)模型可視化系統(tǒng)的研究及應(yīng)用［D］.重慶:重慶交通大學(xué)，2009.

［7］蘇奮振，周成虎，楊曉梅，等.海洋地理信息系統(tǒng):原理、技術(shù)與應(yīng)用［M］.北京:海洋出版社，2005.

［8］楊明遠(yuǎn)，嚴(yán)以新，孔俊，等.珠江口水流泥沙運(yùn)動(dòng)模擬研究［M］.北京:海洋出版社，2008.

［9］MYERS V A.Characteristics of United States hurricanes pertinent to levee design for Lake Okechobeem［R］.FL Hydromet Report 32.Washington D C:Government Printing Office，1954.

［10］葉榮輝.基于GIS的珠江口風(fēng)暴潮數(shù)值模擬系統(tǒng)研究［D］.南京:河海大學(xué)，2009.

［11］歐陽永忠，王瑞，陸秀平，等.VC、VB與 FORTRAN的混合編程技術(shù)及其實(shí)現(xiàn)［J］.海洋測繪，2004，24(1):54-59.(OUYANG Yongzhong，WANG Rui，LU Xiupin，et al.On mixed-language programming and realization with VC，VB and FORTRAN［J］.Hydrographic Surveying and Chartin，2004，24(1):54-59.(in Chinese))

［12］朱泰山，王一一，馮國泰.基于Fortran與C#混編數(shù)值仿真軟件系統(tǒng)的實(shí)現(xiàn)［J］.電腦編程技巧與維護(hù)，2008(11):12-15.(ZHU Taishan，WANG Yiyi，F(xiàn)ENG Guotai.Realization of numerical simulation system based on fortran-C#mixed-language programming［J］.Computer Programming Skills ＆ Maintenance，2008(11):12-15.(in Chinese))

［13］宋關(guān)福，鐘耳順.組件式地理信息系統(tǒng)研究與開發(fā)［J］.中國圖象圖形學(xué)報(bào)，1998，3(4):313-317.(SONG Guanfu，ZHONG Ershun.Research and development of components geographic information systems［J］.Journal of Image and Graphics，1998，3(4):313-317.(in Chinese))

［14］楊旭，黃家柱，許建軍，等.基于組件式GIS的地下水動(dòng)態(tài)管理系統(tǒng)設(shè)計(jì)與開發(fā)［J］.地理與地理信息科學(xué)，2004，20(1):47-50.(YANG Xu，HUANG Jiazhu，XU Jianju，et al.Design and development of groundwater dynamic management system based on ComGIS ［J］.Geography and Geo-Information Science，2004，20(1):47-50.(in Chinese))

［15］楊峰，杜云艷，蘇奮振，等.基于Web服務(wù)的海洋矢量場遠(yuǎn)程可視化研究［J］.地球信息科學(xué)，2008，10(6):749-755.(YANG Feng，DU Yunyan，SU Fenzhen，et al.A study on remote visualization of“ocean vector field data”based on web service［J］.Geo-Information Science，2008，10(6):749-755.(in Chinese))

［16］李杰，于福江，李洋，等.珠江口地區(qū)臺(tái)風(fēng)風(fēng)暴潮的數(shù)值模擬試驗(yàn)［J］.海洋預(yù)報(bào)，2009，26(2):1-6.(LI Jie，YU Fujiang，LI Yang，et al.Numerical simulation of typhoon storm surge in Zhujiang River Estuary［J］.Marine Forecasts，2009，26(2):1-6.(in Chinese))