珠江口风暴潮数值模拟系统的设计与实现

叶荣辉，宋志尧，沈 正，孔 俊，张 卓，金光球

(1.珠江水利委员会珠江水利科学研究院，广东广州 510611;2.河海大学水文水资源及水利工程科学国家重点实验室，江苏南京 210098;3.南京师范大学虚拟地理环境教育部重点实验室，江苏南京 210046)

珠江水系经八大口门注入南海，河口地区水流受径流和潮流的共同作用，在台风的侵袭下，珠江口各口门均会发生程度不同的台风风暴潮灾害，水动力环境较为复杂，因此建立一套能够准确反映珠江口复杂水动力特征的数值模拟系统对于台风风暴潮的预警预报、河口综合治理规划的制定具有现实意义。目前经典的河口、海岸数学模型如 POM［1］，ECOM［2］，FVCOM［3］，ELCIRC［4］等大都采用 Fortran 语言开发，没有提供人机交互界面，模型参数调试、边界条件设置过程繁琐，特别是对于存在多条开边界的复杂研究区域的风暴潮数值模拟，模型参数更多，边界条件给定更为复杂。当计算时段、模型时间步长等发生改变时，对边界条件的手动设置需要耗费大量时间。同时，这些模型的计算结果为一系列文本文件，不清晰直观，对结果的后处理通常需要借助于第三方软件如 Excel，Surfer，Tecplot等［5-7］，而在运用 Surfer，Tecplot等软件对流场进行可视化操作时，还需要熟悉这些商业软件复杂的数据格式，工作量大、效率低。

本文针对珠江口地区风暴潮数值模拟过程复杂繁琐的特点，运用混合编程技术，充分利用现有数学模型代码，结合Fortran语言较强的数值计算能力及地理信息系统(GIS)在可视化方面的优势，设计并建立了珠江口风暴潮数值模拟系统。经过对200814号台风风暴潮的实际模拟应用，结果表明，该系统界面友好、性能稳定、模拟精度较高。

1 系统设计

1.1 功能结构

珠江口风暴潮数值模拟系统分为四大模块。(a)地图操作模块:实现对地图的放大、缩小、漫游、图层控制等功能。(b)台风信息模块:提供台风路径可视化、台风属性信息动态查询等功能。(c)风暴潮模拟计算模块:集成台风模型、珠江口风暴潮数学模型(简称珠江口模型)、南中国海风暴潮数学模型(简称南海模型)等专业数学模型，提供人机交互界面，实现模型参数设置、边界条件设置、风暴潮数值模拟等功能。(d)结果展示模块:提供模型计算结果的后处理功能，主要用于实现重要站点潮位过程线、流场的绘制。

1.2 层次结构

珠江口风暴潮数值模拟系统主要分为3个层次，如图1所示。(a)数据资源层:由空间及非空间数据库组成，主要为该系统提供数据支持。空间数据库主要包含基础地理信息数据及属性信息;非空间数据库包含台风数据库、水文数据库、模型方法库等。(b)专业模型层:在整合台风、水文等数据的基础上，依次调用台风模型计算台风风场，为风暴潮数值模拟提供驱动风场;调用南海模型为珠江口风暴潮数值计算提供开边界;最后调用珠江口模型进行水力计算。(c)服务层:为用户提供模型参数设置、台风信息查询以及计算结果可视化等功能。

2 系统实现

2.1 开发方式

珠江口风暴潮数值模拟系统前台集成开发环境选择Visual Studio 2008，开发语言 C#，数据库平台 MS SQL Server 2008，各专业数学模型均采用Fortran语言开发，编译器采用Visual Fortran 6.5。系统集成基于组件式GIS技术，采用混合编程思路实现，GIS二次开发组件选用ArcGISEngine 9.3。

2.2 模型开发

珠江口风暴潮数值模拟系统中风暴潮模拟计算模块是核心。由于珠江口地区河网密布，水道纵横，分八大口门入海，因此模式计算范围覆盖整个珠江三角洲河网区及口外海区［8］(图2)。模型上游边界分别取在石咀、高要、石角、博罗和老鸦岗;模型的下边界取在珠江口外南海30 m等深线附近。模型采用非结构网格，网格数为94334，网格节点为92221，网格最小步长为20 m。

2.2.1 模型理论

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

图1 系统层次结构Fig.1 System framework structure

其中

式中:u，v——x，y 向流速;t——时间;η——水位;D——全水深;f——柯氏力因子;g——重力加速度;ρ0——水的密度;τwx，τwy——x，y 向 10 m 高空处的风应力;ρa——空气密度，CDs——风应力系数，wx，wy——x，y 向 10 m 高空处的风速，采用 Myers公式［9］计算获得;τbx，τby——x，y 向底部摩擦应力;n——糙率系数。

2.2.2 定解条件

定解条件包括初始条件和边界条件:潮位和流速初始值均设为零。闭边界满足流体不可入条件。上游开边界给定流量过程，由系统从水文数据库中读取历史数据获得;外海开边界给定潮位过程，由系统调用南海模型［10］计算获得。

2.3 混合编程

由于各专业数学模型均采用Fortran语言开发，而系统界面采用C#开发。为充分利用两种语言各自的特点，将各专业数学模型编译为动态链接库(DLL)，以C#为主程序来调用风暴潮计算内核程序。

用混合语言编程存在由于各自语言特点不同而相互不匹配的问题，在C#与Fortran语言混合编程的过程中，应该注意不同语言之间的命名约定、调用约定以及参数传递约定三方面的问题［11-12］，解决这些问题是实现C#与Fortran语言混合编程的关键。

图2 模型计算范围Fig.2 Model’s calculation area

2.4 可视化功能

珠江口风暴潮数值模拟系统采用组件式GIS技术［13-14］实现台风路径、潮位过程线、流场等的可视化功能，对数据库的访问采用ADO.NET数据访问技术，空间信息与属性数据通过ID号进行关联。

台风路径数据存储在台风数据库中，获取数据以后，先生成路径定位点并赋上空间坐标等属性信息，而后根据每个定位点的坐标生成台风路径。

潮位数据包括历史潮位数据及模拟潮位数据。历史潮位数据存储在数据库中，而模拟潮位数据以文本形式存储，系统直接从文本中读取。潮位过程线的绘制通过Microsoft Chart Control 6.0(SP4)控件完成。

流场数据是矢量数据，同时具有数值和方向的特征，不能直接可视化［15］。珠江口风暴潮数值模拟系统首先从数值计算结果文本文件中读取流场信息数据，包括网格空间点坐标、流速、流向等，然后进行矢量数据的映射及绘制和显示。矢量数据的映射采用几何形状的点图标法，用矢量箭头表示水流的方向，用其长度表示流速的大小。矢量箭头的起始点A的坐标为该网格中心点坐标，其余3个顶点B，C，D坐标可参考A点坐标通过简单的几何计算获得(图3)。在程序中矢量箭头被定义为一个Vector类，在读取数据的时候，把所有的数据以一个个Vector实例的形式进行编辑并存储到泛型集合List＜Vector＞Varr中。最后通过ArcEngine中几何对象以要素的形式进行可视化表达，并将同一类型的几何对象以图形图层方式进行管理。

图3 箭头结构示意Fig.3 Schematic diagram of arrow structure

3 系统应用

200814 号强台风引发了珠江口特大风暴潮，珠江口地区灯笼山、黄埔和南沙站等潮位均超过历史极值。本文以该次强台风为例，运用珠江口风暴潮数值模拟系统对该次台风风暴潮进行模拟。

3.1 台风信息查询

在界面中选择需要查看的年份，系统会根据所选年份查找符合条件的数据，并以列表的形式显示。图4为200814号台风路径，从图4可以看出，该台风属西进型台风，珠江口地区位于台风路径前进方向的右半圆区为风暴增水区。通过鼠标点击该场台风路径上离珠江口直线距离较近的点，获取这些时刻台风强度属性信息发现:在离珠江口地区较近的区域，该台风中心气压保持在935 hPa左右，中心最大风速达到了50 m/s。结合其路径，可以从宏观上判断该台风可引起珠江口地区较大的风暴潮增水。

3.2 风暴潮数值模拟

在进行珠江口风暴潮数值模拟前，需要对模型参数进行设置，包括起止时间、时间步长、风应力系数、糙率等。在本次应用中，通过模型参数设置界面，模型计算时间段设为2008年9月22日8:00至2008年9月24日8:00，时间步长设为180s，CDs取为0.0026，n设定为0.023。

图4 台风路径Fig.4 Tyhoon track

本文建立的珠江口模型包含5条上游开边界以及3条外海开边界，若采用手动设置边界条件的方法，设置过程繁琐，工作量大。而在珠江口风暴潮数值模拟系统中，通过界面上的边界计算按钮，即可根据计算时段、时间步长等相关设置自动完成指定时间段的上边界、外海开边界文件的生成。风暴潮模拟计算所需的风压场文件通过调用台风模型生成。在完成对珠江口模型运行所需的文本设置后，运行该模型，模型将自动调用地形文件、边界文本文件、风压场文本文件等进行模拟计算。数值计算结果通过GIS平台统一展现。数值模拟过程如图5所示。

3.3 结果可视化

模拟结果输出时段为2008年9月23日8:00至2008年9月24日8:00，图6为利用珠江口风暴潮数值模拟系统计算得出的200814号台风期间灯笼山站潮位过程线计算值。从图6可以看出，前期由于距离台风中心较远，灯笼山站未发生明显风暴增水，从9月23日16:00开始出现明显的风暴增水，并于9月24日3:30达到整个风暴潮的最高潮位2.66 m。而实测此次风暴潮的最高潮位为2.73 m，出现在9月24日3:00［16］。与实测值相比，模型计算得出的最高潮位误差仅为0.07 m，相位差为30 min，表明该系统模型具有较高的精度。

图7为相应计算时段内某时刻全局流场图，用户可以通过对图层的缩放、平移等功能查看感兴趣区域的流场等内容。

图5 风暴潮数值模拟过程Fig.5 Storm surge numerical simulation process

4 结 论

针对珠江河口地区风暴潮数值模拟过程复杂繁琐的特点，运用混合编程技术，结合Fortran语言较强的数值计算能力及GIS在可视化方面的优势，设计并建立了珠江口风暴潮数值模拟系统。该系统集成了前处理、计算模拟和后处理等功能，可实现台风路径动态显示、潮位过程线绘制及流场可视化，方便用户对计算结果的分析。通过对200814号台风风暴潮的实际模拟测试，表明该系统的建立能有效提高工作效率，用户可借助该系统从宏观与微观尺度上分析、评价台风对珠江口地区水动力环境的影响。

图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］陈可锋，喻国华，徐永辉.水流数值模拟中的混合编程技术［J］.水运工程，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］胡艳进.平面二维水流数学模型可视化系统的研究及应用［D］.重庆:重庆交通大学，2009.

［7］苏奋振，周成虎，杨晓梅，等.海洋地理信息系统:原理、技术与应用［M］.北京:海洋出版社，2005.

［8］杨明远，严以新，孔俊，等.珠江口水流泥沙运动模拟研究［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的珠江口风暴潮数值模拟系统研究［D］.南京:河海大学，2009.

［11］欧阳永忠，王瑞，陆秀平，等.VC、VB与 FORTRAN的混合编程技术及其实现［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#混编数值仿真软件系统的实现［J］.电脑编程技巧与维护，2008(11):12-15.(ZHU Taishan，WANG Yiyi，FENG 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］宋关福，钟耳顺.组件式地理信息系统研究与开发［J］.中国图象图形学报，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的地下水动态管理系统设计与开发［J］.地理与地理信息科学，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服务的海洋矢量场远程可视化研究［J］.地球信息科学，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］李杰，于福江，李洋，等.珠江口地区台风风暴潮的数值模拟试验［J］.海洋预报，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))