技术领域
[0001] 本
发明属于海洋技术领域，涉及一种基于人工智能和数值模型的
风暴潮计算方法。
背景技术
[0002] 由于海洋状况和
水深地形的影响，现有的风暴潮模型对于每个地区的适用度差异很大，对不同海域的波浪模拟效果也大不相同，对于复杂海况海域的风暴潮模拟效果并不好。近年来，由于
气候变化的影响，极端恶劣天气状况的增加使得台风、极端海浪等要素出现的
频率增加，台风、波浪和风暴潮的研究计算已经成为海洋工程领域最为重要的课题之一。我国沿海海岸线长达1.8万公里，尤其是东南沿海地区经常遭受台风和海浪的侵害，对人民的生产生活都造成了极大的威胁。自上世纪五十年代计算机技术飞速发展以来，有限元计算方法自出现至今经历了极大的改进。然而复杂的海况和天气导致人们对于基于数值方法的风暴潮模拟到达了一定的
瓶颈期，模拟效果都没有达到对每个海域和每种海况都理想的效果。
发明内容
[0003] 本发明的目的在于提供一种基于人工智能和数值模型的风暴潮计算方法，本发明的有益效果是能实现根据海洋状况进行风暴潮模型自动选择、校核风暴潮相关特性参数和有限元方法的计算参数，在此
基础上进行风暴潮模拟，达到风暴潮模型对所有地区都能相对适应的效果，更加准确的模拟复杂海况的海域潮位变化。
[0004] 本发明所采用的技术方案是假定地球的半径远大于
海水的厚度，球
坐标系下的连续方程和动量方程表达形式为：
[0005]
[0006]
[0007]
[0008] 式中，地球平均半径R＝6.378×106m，λ和φ为经度和纬度；U和V为海水深度平均处的水平流速(U为东西方向，V为南北方向)， vdz为垂向积分的速度分量，uv为垂向变化的速度分量，总水深H≡ξ+h，其中h为海底至平均海平面的距离，ζ为自由海平面至平均海平面的距离；t为时间，
[0009] 科氏
力参数f＝2Ωsinφ，Ω＝7.292×10-5rads-1；Ps为海面上的
大气压强，g为重力
加速度，ρ0为海水
密度，η为
牛顿等价潮汐势；
[0010]
[0011]
[0012] 将以(λ,φ)为单位的球坐标系通过标准圆柱形投影(Carte Parallel ogrammatique Projection,CPP)向笛卡尔坐标系(λ0,φ0)进行转换，控制方程变为：
[0013]
[0014]
[0015]
[0016] 式中，(λ0,φ0)为计算区域中心点的经度和纬度，x＝R(λ-λ0)cosφ0，y＝Rφ；球坐标修正因子SP＝cosφ0/cosφ；
[0017]
[0018]
[0019] 通过量级分析可以发现，当研究区域不接近极地区域时，式(4)～(6)的最后一项忽略：
[0020]
[0021]
[0022]
[0023] 除了x坐标轴的导数要乘以球坐标修正因子SP，上述方程中的各项等价于笛卡尔坐标系中的各项。故式(7)～(9)是ADCIRC模型的计算通式，在SP＝1时等于笛卡尔坐标系下的方程组。
[0024] 球坐标系下的通用
波动连续性方程表达式为：
[0025]
[0026] 式中，τ0为随空间变化的时间权重系数，分项Ax、Ay的表达式如下：
[0027]
[0028]
[0029] 底
摩擦系数和风拖曳系数表达式分别为：
[0030]
[0031] Cd＝μ(0.75+0.667W) (14)
[0032] 其中Cfmin和Hbreak是常数，分别为线性底摩擦系数和
破碎波高。a0和b0均为随环境条件变化的常数。μ为摩擦系数，W是海面上10m处的风速。
附图说明
[0033] 图1是方法流程示意图。
具体实施方式
[0034] 下面结合具体实施方式对本发明进行详细说明。
[0035] 一、如图1所示本发明方法流程，建立控制方程
[0036] 应用静水压力和Bousinessq近似，并假定地球的半径远大于海水的厚度，球坐标系下的连续方程和动量方程表达形式为：
[0037]
[0038]
[0039]
[0040] 式中，地球平均半径R＝6.378×106m，λ和φ为经度和纬度； vdz为垂向积分的速度分量，总水深H≡ζ+h，其中h为海底至平均海平面的距离，ζ为自由海平面至平均海平面的距离；t为时间， 科氏力参数f＝2Ωsinφ，Ω＝7.292×10-5rads-1；Ps为海面上的大气压强，g为
重力加速度，ρ0为海水密度，η为牛顿等价潮汐势；
[0041]
[0042]
[0043] 将以(λ,φ)为单位的球坐标系通过标准圆柱形投影(Carte Parallel ogrammatiqueProjection,CPP)向笛卡尔坐标系(λ0,φ0)进行转换，控制方程变为：
[0044]
[0045]
[0046]
[0047] 式中，x＝R(λ–λ0)cosφ0，y＝Rφ；球坐标修正因子SP＝cosφ0/cosφ；
[0048]
[0049]
[0050] 通过量级分析可以发现，当研究区域不接近极地区域时，式(4)～(6)的最后一项可以忽略：
[0051]
[0052]
[0053]
[0054] 除了x坐标轴的导数要乘以球坐标修正因子SP，上述方程中的各项等价于笛卡尔坐标系中的各项。故式(7)～(9)是ADCIRC模型的计算通式，在SP＝1时等于笛卡尔坐标系下的方程组。
[0055] 为了避免方程原始的Galerkin有限元形式会产生诸多振荡和不守恒等问题，ADCIRC模型采用通用波动连读方程GWCE。Kinnmark发展的GWCE方程的基本思想是：初始深度积分的连续方程求取时间微分后，与求取空间微分的深度平均的初始动量方程结合，然后乘以权重系数。联立GWCE方程与初始动量方程，转换后的方程能够准确求解有限差分格式下的水位和流速。
[0056] 根据上述GWCE方程的转化思想，球坐标系下的通用波动连续性方程表达式为：
[0057]
[0058] 式中，τ0为随空间变化的时间权重系数，分项Ax、Ay的表达式如下：
[0059]
[0060]
[0061] 与原有连续方程相比，GWCE方程的计算效率更高，计算
精度也得以提高，同时归功于自然解耦式的求解，方程的
质量矩阵非定常，使得模型的
算法简单，计算耗时有了缩短，计算出错的几率也有了很大改善。
[0062] 最后，ADCIRC模型对控制方程联立求解。GWCE方程可以通过计算一致或集中质量矩阵(编译标志)、显式或隐式的时间推进方案(时间权重系数)进行求解。如果是集中的全显式公式，则无需求解矩阵。
[0063] 同时，ADCIRC模型在实际运行过程中，底摩擦系数和风拖曳系数也是重要[0064] 参数，表达式分别为：
[0065]
[0066] Cd＝μ(0.75+0.667V) (14)
[0067] 上述两个参数虽然有明确表达式，但在实际海况中是很难计算其精确值的，比如底摩擦系数与线性底摩擦系数Cfmin和破碎波高Hbreak有关，但这两个常数却很难测量。
[0068] 因此，选取了模型方程中的四个参数，随空间变化的时间权重系数τ0、底摩擦系数Cd、风拖曳系数Cf和阻尼时间T0来进行智能化处理。
[0069] 二、人工智能模型的结合：
[0070]
遗传算法(GeneticAlgorithms)以达尔文的“物竞天择，适者生存”的
生物进化论为原理，并将其引入优化参数形成的编码
串联群体中，按照适应度函数评估每个个体的适应度，并通过遗传中的选择、交叉和变异操作对个体进行筛选，使适应度值较好的个体被保留，适应度较差的个体逐渐被淘汰，新产生的群体在保留上一代的信息的基础上又要优于上一代。如此循环往复，直至满足要求。遗传算法基本的操作分为三个部分：即选择操作、交叉操作及变异操作。
[0071] 1、个体编码及初始化种群：每一个体均为一个由4个系数(即随空间变化的时间权重系数τ0、底摩擦系数Cd、风拖曳系数Cf和阻尼时间T0)组成的取值组合，即个体长度为4；根据个体编码方式、长度及种群的规模对种群进行初始化操作，即生产一S×4的二维随机数组储存每个网络的权值、
阈值作为初始种群，其中S为种群个数；
[0072] (2)计算每一个体的适应度：根据每一个体的初始权值和阈值，用
训练数据训练波浪模型，以预测输出和期望输出之间的误差平方和的倒数作为个体适应度值，计算公式如式(15)所示：
[0073]
[0074] 式中，Fi为第i个个体的适应度值；Ei为第i个个体的误差平方和；yj、oj分别为第j组波浪数据的期望输出、预测输出；m为波浪数据个数；k为系数，取0.5。
[0075] (3)选择操作：选择操作是指从父代群体中以一定概率选择个体到子代群体中，个体被选中的概率跟适应度值有关，适应度值大的个体被选中的概率也越大。为确保性能优良的个体能够遗传到下一代中，本文使用轮盘赌法作为选择算法，其公式如式(16)所示：
[0076]
[0077] 式中，pi为每个个体被选中的概率；其余符号含义与上文相同。
[0078] (4)交叉操作：交叉操作为遗传算法中最重要的操作。种群中个体间通过交叉产生新的个体，可逐渐扩展搜索空间，提高全局搜索的能力。本文交叉操作算法采用实数交叉法，该方法可确保父代的优良基因模式基本不发生破坏，提高种群的平均适应度。
[0079]
[0080] 式中，akj为k个个体的j位的实数编码(基因)；alj为l个个体的j位的实数编码(基因)；b为[0,1]之间的随机数。
[0081] (5)变异操作：变异操作是遗传变化的主要根源，是摆脱遗传算法的局部收敛的最有效方法。本文按式(5)进行变异操作，该算法随着遗传代数的增加而逐渐减少变异操作，可保护优良的基因模式，利于趋向最优解。
[0082]
[0083] 式中，aij为第i个个体amax为基因aij的上界；amin为基因aij的下界；f(g)＝r2(1-g/Gmax)；r2为一个随机数；g为当前遗传代数；Gmax为最大遗传代数；r为[0,1]之间随机数。
[0084] 以上所述仅是对本发明的较佳实施方式而已，并非对本发明作任何形式上的限制，凡是依据本发明的技术实质对以上实施方式所做的任何简单
修改，等同变化与修饰，均属于本发明技术方案的范围内。
