技术领域
[0001] 本
发明涉及海洋科学及海洋工程技术领域，具体涉及一种基于GIS的
风暴潮漫滩及淹没模型。
背景技术
[0002] 风暴潮是来自海上的一种巨大的
自然灾害现象，系指由于强烈的大气扰动——如强风和气压骤变所招致的海面异常升高的现象。风暴潮常伴随天文潮(正常潮位)、短周期(数秒)的海浪而来，使潮位暴涨，甚至
海水漫溢，酿成灾害。
[0003] 风暴潮可分为两类，即台风风暴潮和温带风暴潮。我国是一个自然灾害发生较为频繁的国家，尤其是海洋灾害，而风暴潮灾害居海洋灾害首位，常发生在经济较为发达的沿海地区。我国东部濒临渤海、黄海、东海，南部为南海，海岸线长达18000km，沿岸地区经常遭受风暴潮侵袭，往往浅海湾更容易发生风暴潮，作用时间通常会超过30个小时，是世界上风暴潮灾害最为严重的国家之一。
[0004] 因此，基于并行数值计算和GIS技术针对滨海城市的精细化风暴潮模拟，再采用GIS技术对风暴潮灾害漫滩过程中可能影响范围和程度进行预测，能够为灾害预警、预报、救灾等提供有效的技术支持。
[0005] 但是，传统的对风暴潮的模拟较为粗糙，未能达到理想的
精度，传统方式也没有利用可能最大风暴潮水位和海堤高程计算风暴潮漫堤或溃堤的过水量，因此对风暴潮灾害的漫滩过程的影响范围和程度预测的精度也并不高，有较大提升空间。
发明内容
[0006] 本发明要解决的技术问题在于，针对
现有技术的
缺陷，提供一种基于GIS技术建立风暴潮漫滩和淹没分析的方法及系统。
[0007] 本发明解决其技术问题所采用的技术方案是：构造一种基于GIS技术建立风暴潮漫滩和淹没分析的方法，具体包括以下步骤：
[0008] S1、构建地理信息、历史台风
数据库；
[0009] S2、利用数据库中的历史台风数据，采用台风风场经验模型，计算历史
热带气旋风场和气压场；
[0010] S3、基于ADCIRC海洋模式，利用步骤S2计算所得的历史热带气旋风场和气压场，在研究区域范围内，建立风暴潮精细化数值模型；
[0011] S4、在研究区域内，构造极端热带气旋，并利用步骤S3建立的建立风暴潮精细化数值模型，计算得到当前环境下该研究区域的模拟风暴潮水位；
[0012] S5、根据步骤S4中计算的模拟风暴潮水位，以及海堤高程，计算风暴潮漫堤和溃堤时的过水量；
[0013] S6、利用数据库中的地理信息数据，根据步骤S5计算所得的风暴潮漫堤和溃堤时的过水量，采用体积法计算风暴潮漫堤和溃堤时的淹没区水位，并根据计算所得的淹没区水位，利用GIS技术直观展现漫滩和淹没模拟结果。
[0014] 进一步的，在执行步骤S4之前，还需结合步骤S2计算所得的历史热带气旋风场和气压场，对步骤S3构建的风暴潮精细化数值模型的计算精度进行调整；其中：
[0015] 将历史热带气旋风场A和气压场B作为风暴模式的气象驱动场，将参数A和B带入到风暴潮精细化数值模型中，计算研究区域内，不同强度的热带气旋所引起的风暴潮水位；将所得的风暴潮水位与历史台风数据库中的数据进行比对，在有误差的情况下，返回到步骤S2，重新计算历史热带气旋风场和气压场，直到得到一个精准的风暴潮水位计算模型的时候，执行步骤S4。
[0016] 进一步的，步骤S5中，需在确定淹没系数ks和流量修正系数cv后，进一步求得淹没区过水量Qb，其中，淹没系数ks、流量修正系数cv的计算公式分别为：
[0017]
[0018]
[0019] 其中，h为堤坝一侧的潮水高度，h1是淹没区的水位高度，hb为溃口处高程，而其中当水流以自由出流的方式入流时，ks＝1.0；Bd为坝址处的坝宽，hbm为最终的溃口底高程，且在hbm＝hb时，淹没区过水量Qb的计算公式为：
[0020]
[0021] 其中，z为溃口边坡系数，bi为瞬间溃口宽。
[0022] 进一步的，步骤S6中，采用体积法计算风暴潮漫堤和溃堤时的淹没区水位，具体为使用固定水量模型，求解洪水水面高程EW；所述固定水量模型的数学表达式为：
[0023]
[0024] 其中，EW为洪水水面高程；V为风暴潮造成的淹没区域洪
水体积；将整个风暴潮洪水淹没区A分割为若干个小方
块，则Eg(i)为第i个方块的水面高程；Δσ为每个小方块的面积；N为小方块的总个数，N和A的关系为 N＝f(Ew，Eg)；i＝1，2...，N。
[0025] 进一步的，步骤S6中，采用体积法计算洪水水面高程EW，首先，需输入步骤S5计算所得的淹没区过水量Qb，以及研究区域的高程数据到固定水量模型；然后，将固定水量模型转换为：
[0026]
[0027] 最后针对上述计算公式，利用二分法求解洪水水面高程Ew。
[0028] 进一步的，根据计算所得的洪水水面高程Ew，一方面
结合水面高程Eg，计算洪水淹没水深El；另一方面通过Ew在反求到淹没区的分割区域总数N后，利用N反求淹没范围即所述整个风暴潮洪水淹没区A；
[0029] 其中，根据计算所得的洪水淹没水深El和淹没范围A，利用GIS技术的三维
可视化展现功能，使用渐变色标识表征，来评估区域内不同等级强度风暴潮情景下被淹没范围及水深分布；通过绘制不同等级强度台风引起的风暴潮淹没范围及水深分布图，对漫滩和淹没模拟结果进行直观显示。
[0030] 本发明公开的一种基于GIS技术建立风暴潮漫滩和淹没分析的系统，包括以下模块：
[0031] 数据构建模块，用于构建地理信息、历史台风数据库；
[0032] 风场和气压场计算模块，用于利用数据构建模块中构建的历史台风数据，采用台风风场经验模型，计算历史热带气旋风场和气压场；
[0033] 风暴潮精细化数值模型建立模块，用于基于ADCIRC海洋模式，以及风场和气压场计算模块计算所得的历史热带气旋风场和气压场，在研究区域范围内，建立风暴潮精细化数值模型；
[0034] 风暴潮水位计算模块，用于在研究区域内，构造极端的热带气旋，利用风暴潮精细化数值模型建立模块建立的风暴潮精细化数值模型，计算得到该研究区域的模拟风暴潮水位；
[0035] 过水量计算模块，根据风暴潮水位计算模块中计算的模拟风暴潮水位，以及海堤高程，计算风暴潮漫堤和溃堤时的过水量；
[0036] 结果展示模块，用于利用数据构建模块中构建的地理信息数据，根据过水量计算模块计算所得的风暴潮漫堤和溃堤时的过水量，采用体积法计算风暴潮漫堤和溃堤时的淹没区水位，并根据计算所得的淹没区水位，利用GIS技术直观展现漫滩和淹没模拟结果。
[0037] 在本发明所述的一种基于GIS技术建立风暴潮漫滩和淹没分析的方法及系统中，通过设计极端风暴潮，计算海堤边界处风暴潮水位，进一步计算漫堤过水量，结合高精度的DEM和地物类型信息，将计算得到的过水量利用GIS技术，建立漫滩和淹没分析模型，模拟出风暴潮漫滩和淹没过程及程度，为将来应对未知风暴潮可能产生的漫滩过程和范围提供了可靠的技术支持。
[0038] 实施本发明的一种基于GIS技术建立风暴潮漫滩和淹没分析的方法及系统，具有以下有益效果为，其不同于传统的风暴潮漫滩及淹没模型，本发明创新性的计算可能最大风暴潮水位和海底高程计算风暴潮漫堤或溃堤的过水量，并且结合高精度的DEM和地物类型信息，采用固定水量模型，即“体积法”，再利用GIS技术建立漫滩和淹没分析模型，通过网格
染色方法，直观展现出漫滩和淹没的模拟结果。能够为灾害预警、预报等提供有效的技术支持。
附图说明
[0039] 下面将结合附图及
实施例对本发明作进一步说明，附图中：
[0040] 图1是本发明公开的风暴潮漫滩和淹没分析方法的
流程图；
[0041] 图2是本发明公开的利用“体积法”计算洪水淹没水深和淹没范围的方法流程图；
[0042] 图3是本发明公开的利用GIS技术展示输出结果的效果图；
[0043] 图4是本发明公开的建立风暴潮漫滩和淹没分析的系统结构图。
具体实施方式
[0044] 为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本发明的具体实施方式。
[0045] 请参考图1，其为本发明公开的风暴潮漫滩和淹没分析方法的流程图，该方法具体包括以下步骤：
[0046] S1、确定研究区域，并
整理该区域内的
基础地理信息数据和历史热带气旋数据，构建地理信息、历史台风数据库；其中，研究区域的地理信息数据库包括高精度DEM、海域水深分布、沿岸提防数据、社区或村人口分布、土地利用现状二级分类空间分布数据、重要承载体数据；研究区域的历史台风数据库内容包括台风等级、潮位站历史观测数据、历史典型风暴潮淹没范围数据。
[0047] S2、利用数据库中的历史台风数据，采用台风风场经验模型，对研究区域历史热带气旋数据进行重构，具体为计算历史热带气旋风场和气压场；
[0048] 由于在风暴潮的计算过程中，对于气旋的风场W和气压场Pa的计算是一个重要环节；其中，台风域中的风场由两个矢量场
叠加而成，其一是相对台风中心对称的风场，其二是基本风场。本实施例中，选用长期在业务化预报中使用的Jelesnianski模型作为台风风场经验模型，其公式如下：
[0049]
[0050]
[0051] 其中，A＝-[(x-xc)sinθ+(y-yc)cosθ]，B＝[(x-xc)cosθ-(y-yc)sinθ]，W为台风最大风速；R为台风最大风速半径；r为计算点到台风中心的距离；V0为台风移动速度；(x,y)、(xc,yc)分别为计算点坐标和台风中心坐标；θ为入流
角；P0为台风中心气压，P∞为无穷远处气压。
[0052] 上述公式中，台风的最大风速半径R，综合考虑台风探测资料和多种经验计算公式确定，选取被普遍认可的一个最大风速半径经验统计公式：
[0053] R＝Rk-0.4×(P0-900)+0.01×(P0-900)2；
[0054] 其中Rk为经验常数，通常取[30,60]。
[0055] 上述公式中，台风中心气压P0的计算采用概率论法计算，利用中国气象局(CMA)1949-2015年热带气旋最佳路径集(Best-track)数据，选取出研究区域400km范围内历年路经本区域的台风，以这些台风的最小P0值作样本。采用极值I型分布计算1000年一遇的P0值作为可能最大台风的中心气压。
[0056] S3、基于ADCIRC海洋模式，利用步骤S2计算所得的历史热带气旋风场和气压场，在研究区域范围内，建立风暴潮精细化数值模型；具体为：
[0057] 确定研究区域的计算范围，配置计算网格和开边界条件，选取ADCIRC海洋模式建立研究区域风暴潮数值计算模型，由于ADCIRC海洋模式可模拟海洋、近岸与河口的水位、流场等，它基于有限元方法，采用可任意局部灵活加密的无结构网格，且ADCIRC模式的计算速度相对较快。
[0058] ADCIRC海洋模式连续方程为：
[0059]
[0060] 其动量方程为：
[0061]
[0062]
[0063] 其中，t为时间；(x,y)为水平笛卡尔坐标；(λ,φ)为经度和纬度；(λ0,φ0)为计算网格中心点的经度和纬度；H＝ζ+h为
海水水柱的总水深，单位为m；ζ为从平均海面起算的自由表面高度；h(x,y)为未扰动的海洋水深，即平均海平面至海底的距离，单位为m；R为地球半径，单位为m，取6378135m；(U,V)为深度平均的海水水平流速，单位为ms-1；f＝2Ωsinφ为科式
力参数，单位为s-1；Ω为地球自转
角速度；g为重力
加速度，单位为ms-2；ρ0为海水
密度，当前模式取ρ0为1025kgm-3；Ps为海水自由表面处的
大气压强，单位为Nm-2；η为
牛顿引潮势，单位为m；τbx,τby为海底
摩擦力的x和y方向的分量；Dx,Dy为动量方程的水平扩散项。
[0064] 作为一个优选的实施例，在执行步骤S4之前，还需结合步骤S2计算所得的历史热带气旋风场和气压场，对步骤S3构建的风暴潮精细化数值模型的计算精度进行调整；其中：
[0065] 将历史热带气旋风场A和气压场B作为风暴模式的气象驱动场，将参数A和B带入到风暴潮精细化数值模型中，计算研究区域内，不同强度的热带气旋所引起的风暴潮水位；将所得的风暴潮水位与历史台风数据库中的数据进行比对，在有误差的情况下，返回到步骤S2，重新计算历史热带气旋风场和气压场，直到得到一个精准的风暴潮水位计算模型的时候，执行步骤S4。
[0066] S4、在研究区域内，构造极端的热带气旋，利用步骤S3建立的建立风暴潮精细化数值模型，计算得到该研究区域的模拟风暴潮水位；其中，模拟风暴潮水位的计算需要考虑最大增水域天文高潮位(10％超越概率天文高潮位值)的叠加效应。
[0067] 本实施例依据研究区域潮位站的潮位调和常数，分别预报多年的逐时天文潮潮位，逐月挑选天文最高潮，然后计算高潮位的累积
频率，累积频率曲线。从中获取最高天文潮位值、历年天文潮平均高潮位值、10％超越概率天文高潮位值。在天文潮预报潮位序列中找到10％超越概率天文高潮位值相接近的潮位，并假定该时刻与最大增水发生的时刻相同，根据台风路径反向计算模拟起始时间，进行风暴潮和天文潮的耦合计算。
[0068] S5、根据步骤S4中计算的模拟风暴潮水位，以及海堤高程，计算风暴潮漫堤和溃堤时的过水量；其中，潮水漫堤过水量的计算公式为：
[0069]
[0070] 式中：ε为侧收缩系数，取1.0；b为漫堤长度，单位为m；H为堰顶
水头，单位为m，其中，H＝G(潮高)-P(堤顶高)；m为流量系数，由表1查算所得：
[0071] 表1 流量系数表
[0072]
[0073] 在无人为干扰的情况下模拟溃坝，计算淹没区的过水量。采用Fread提出的宽顶堰流公式计算风暴潮造成的决口处的过水量，这个模型的前提假设条件有两个，一个是堤坝属于逐渐溃堤，另一个是决堤口的断面为梯形。溃口处过水量Qb的计算过程如下，首先Qb的初始数学表达式为：
[0074] Qb＝cvks[3.1bi(h-hb)]1.5+2.45z(h-hb)2.5；
[0075] 其中，cv为流量修正系数；z为溃口边坡系数，一般取值0-2之间，取值大小主要与堤坝材料及材料的密实度等有关；ks为淹没系数；bi为瞬间溃口宽，单位为m；h为堤坝一侧的潮水高度，单位为m；hb为溃口处高程；其中淹没系数ks的计算公式如下：
[0076]
[0077] 其中，ks适用于 的情况，否则ks＝1.0。这是因为在实际情况中，当淹没区水位不断涨高时，溃口处过流能力减弱，流量减小，这种情况为淹没出流，反之为自由出流。故当 时判定为淹没出流，否则以自由出流方式入流，即ks＝1.0。式中h1是淹没区的水位高度。
[0078] 流量修正系数的计算公式为：
[0079]
[0080] 其中，Bd为坝址处坝宽，单位为m；hbm为最终溃口底高程，单位为m。
[0081] 假设hbm最终溃口底高程与hb相等，故由Qb和cv的计算公式，得到最终的流量计算公式为：
[0082]
[0083] S6、利用数据库中的地理信息数据，根据步骤S5计算所得的风暴潮漫堤和溃堤时的过水量，采用体积法计算风暴潮漫堤和溃堤时的淹没区水位，并根据计算所得的淹没区水位，利用GIS技术直观展现漫滩和淹没模拟结果。
[0084] “体积法”的基本原理为设定淹没区
位置为(x,y)，Ew(x,y)表示(x,y)处洪水水面高程，Eg(x,y)表示(x,y)处地面高程，H(x,y)表示淹没处(x,y)的水深。则得到以下的关系：
[0085] H(x,y)＝Ew(x,y)-Eg(x,y)
[0086] 采用体积法计算风暴潮漫堤和溃堤时的淹没区水位，具体为使用固定水量模型，求解洪水水面高程EW；所述固定水量模型的数学表达式为：
[0087]
[0088] 其中，EW为洪水水面高程；V为风暴潮造成的淹没区域洪水体积；将整个风暴潮洪水淹没区A分割为若干个小方块，则Eg(i)为第i个方块的水面高程；Δσ为每个小方块的面积；N为小方块的总个数，N和A的关系为 N＝f(Ew，Eg)；i＝1，...，N。
[0089] 采用体积法计算洪水水面高程EW的过程中，首先，需输入步骤S5计算所得的淹没区过水量Qb，以及研究区域的高程数据到固定水量模型；然后，将固定水量模型转换为：
[0090]
[0091] 最后针对上述计算公式，利用二分法求解洪水水面高程Ew。其中，二分法是指对于函数f(x)，如果存在实数c，使得f(x)＝0，则x＝c为函数的零点，假定零点位于区间(a,b)之间，每次取区间中间值对应的函数值与零值进行比较，不断缩小范围从而求得零值或者零值的近似值。
[0092] 根据以上分析，已知洪水淹没量及固定体积V，需要将固定体积与淹没区洪水体积进行判定，故设定求解函数：
[0093]
[0094] 由于该函数为单调递减函数，又已知的f(Ew0)＝V，Ew0为入口单元处的高程，现求取一个E'w，使f(Ew)趋近于0，当f(Ew)收敛于某收敛容差时，此时的E'w即为待求的风暴潮洪水水面高程；其中，在二分法的求解中，需要先求得一个水面高程Ew1，使f(Ew1)＜0，然后利用二分法在(Ew0,Ew1)的区间范围内寻找使f(Ew)趋近于0的E'w。
[0095] 最后，根据计算所得的洪水水面高程Ew，一方面结合水面高程Eg，计算洪水淹没水深El；另一方面通过Ew在反求到淹没区的分割区域总数N后，利用N反求淹没范围即所述整个风暴潮洪水淹没区A；
[0096] 其中，根据计算所得的洪水淹没水深El和淹没范围A，利用GIS技术的三维可视化展现功能，使用渐变色标识表征，来评估区域内不同等级强度风暴潮情景下被淹没范围及水深分布；通过绘制不同等级强度台风引起的风暴潮淹没范围及水深分布图，对漫滩和淹没模拟结果进行直观显示。
[0097] 故，请参考图2，其为本发明公开的利用“体积法”计算洪水淹没水深和淹没范围的方法流程图，其计算步骤如下：
[0098] S61、输入淹没水量数值及研究区的高程数据到固定水量淹没模型；
[0099] S62、针对以下计算公式 利用二分法求解洪水水面高程Ew；
[0100] S63、根据所求的洪水水面高程Ew，反推得到洪水淹没水深及淹没范围；
[0101] S64、利用渐变色标识表征评估区域内不同等级强度风暴潮情景下被淹没范围及水深分布，绘制不同等级强度台风引起的风暴潮淹没范围及水深分布图(如图3所示)。
[0102] 请参考图4，其为本发明公开的建立风暴潮漫滩和淹没分析的系统结构图，该系统包括数据构建模块L1、风场和气压场计算模块L2、风暴潮精细化数值模型建立模块L3、风暴潮水位计算模块L4、过水量计算模块L5、结果展示模块L6：
[0103] 数据构建模块L1用于构建地理信息、历史台风数据库；
[0104] 风场和气压场计算模块L2用于利用数据构建模块中构建的历史台风数据，采用台风风场经验模型，计算历史热带气旋风场和气压场；
[0105] 风暴潮精细化数值模型建立模块L3用于基于ADCIRC海洋模式，以及风场和气压场计算模块计算所得的历史热带气旋风场和气压场，在研究区域范围内，建立风暴潮精细化数值模型；
[0106] 风暴潮水位计算模块L4用于在研究区域内，构造极端的热带气旋，利用风暴潮精细化数值模型建立模块建立的风暴潮精细化数值模型，计算得到该研究区域的模拟风暴潮水位；
[0107] 过水量计算模块L5根据风暴潮水位计算模块中计算的模拟风暴潮水位，以及海堤高程，计算风暴潮漫堤和溃堤时的过水量；
[0108] 结果展示模块L6用于利用数据构建模块中构建的地理信息数据，根据过水量计算模块计算所得的风暴潮漫堤和溃堤时的过水量，采用体积法计算风暴潮漫堤和溃堤时的淹没区水位，并根据计算所得的淹没区水位，利用GIS技术直观展现漫滩和淹没模拟结果。
[0109] 上面结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和
权利要求所保护的范围情况下，还可做出很多形式，这些均属于本发明的保护之内。
