本发明涉及海洋大气边界层领域，公开了一种海面空气动力学粗糙度的计算方法，包括如下步骤：获取参数数据；建立初始参数组合区间并进行采样，在不同组合区间内，使用敏感性分析方法得到各种参数化方案在不同参数组合区间内的敏感性；依据敏感性结果，重新建立新的参数组合区间；分析新的参数组合区间内各参数化方案的适用性，得到当前参数组合区间内适用性最高的方案；使用智能优化算法在不同的参数组合区间内寻找方案的最优系数，得到计算粗糙度的最优方法。本发明所公开的方法可以提高海面空气动力学粗糙度的计算精度，为海上电磁波传播路径诊断、通信链路信道分析、蒸发波导诊断、雷达目标探测和遥感信息获取等方面提供技术(56)对比文件韩佳彤;彭怀午;许昌;刘艺芳.海上风电场的海面粗糙度模型研究.内蒙古工业大学学报(自然科学版).2017,(05),第63-67页.Sittichoke Pookpunt 等.Design ofoptimal wind farm configuration using abinary particle swarm optimization atHuasai district, Southern Thailand.《Energy Conversion and Management》.2016,第160-180页.
1.一种海面空气动力学粗糙度的计算方法，其特征在于，包括如下步骤：步骤一，数据采集：
在海上观测平台安装风速风向仪、三维超声风、毫米波雷达获取海平面10m高度处的平均风速 、海平面10m高度处的三维风速( 、、)、波高 和波周期 参数数据，其中，u为纬向风速，v为经向风速，w为垂直风速；
依据海平面10m高度处的平均风速 、波高 和波周期 参数数据的取值范围，设置初始间隔、划分初始参数区间并组合，建立初始参数组合区间；在初始参数组合区间内采用蒙特卡洛采样法，对每个组合区间内的参数进行采样，获取所需参数样本；在不同组合区间内，使用敏感性分析方法利用参数样本计算海面空气动力学粗糙度参数化方案的一阶敏感性指数和全局敏感性指数，得到各种海面空气动力学粗糙度参数化方案在不同参数组合区间内，对于三种参数的敏感性；依据敏感性结果，拆分或合并初始组合区间，重新建立新的参数组合区间；海面空气动力学粗糙度参数化方案包括YT96、TY01、O02、GW06、PS07；
在新的参数组合区间内，依据实测海平面10m高度处的三维风速( 、、)，以涡动相关法计算的摩擦速度 和不同海面空气动力学粗糙度参数化方案结合CAORE3.0模型计算的摩擦速度 为依据，对比不同方案的均方根误差 ，分析新的参数组合区间内各海面空气动力学粗糙度参数化方案的适用性，得到当前参数组合区间内适用性最高的海面空气动力学粗糙度参数化方案；
其中， 为纬向风速脉动量，为经向风速脉动量， 为垂直方向风速脉动量， 表示纬向和经向脉动量在20秒内的平均值的乘积；脉动量计算公式如下所示：；
其中，为纬向风速，为经向风速，为垂直方向风速，、、为30分钟内 的平均值；
在不同数据区间内以 为标准，计算五种海面空气动力学粗糙度参数化方案的摩擦速度和涡动相关法计算的摩擦速度之间的均方根误差 、 、 、 ，每种海面空气动力学粗糙度参数化方案单独计算，均方根误差的计算方法如下所示：；
其中，n为样本总数， 为第m个样本使用涡动相关法计算的摩擦速度， 为第m个样本使用五种海面空气动力学粗糙度参数化方案中每个方案单独计算的摩擦速度；
选择均方根误差最小的方案作为当前参数组合区间内适用性最高的海面空气动力学粗糙度参数化方案；
使用智能优化算法在不同的参数组合区间内，依据适应度函数值逐步迭代，寻找当前参数组合区间内适用性最高的海面空气动力学粗糙度参数化方案的最优系数，得到一种在不同参数组合区间内使用拥有最优系数的不同海面空气动力学粗糙度参数化方案计算粗糙度的方法。
2.根据权利要求1所述的一种海面空气动力学粗糙度的计算方法，其特征在于，敏感性分析方法包括Sobol指数法和扩展傅里叶振幅灵敏度检验法。
3.根据权利要求1所述的一种海面空气动力学粗糙度的计算方法，其特征在于，敏感性参数包含海平面10m高度处的平均风速 、波高 和波周期 ，参数取值范围、 、 ， 、、 的单位分别为m/s、m、s； 、、 的间隔均设置为1。
4.根据权利要求1所述的一种海面空气动力学粗糙度的计算方法，其特征在于，步骤二中，依据各模型某个输入参数在同一区间的最大全局敏感性指数作为区间划分的判断依据，将参数敏感性指数较大的参数区间进行拆分，将敏感性指数较小的区间进行合并。
5.根据权利要求1所述的一种海面空气动力学粗糙度的计算方法，其特征在于，步骤四中，智能优化算法包括粒子群算法、遗传算法、蚁群算法、模拟退火算法中的一种。
6.根据权利要求5所述的一种海面空气动力学粗糙度的计算方法，其特征在于，步骤四中，使用粒子群算法进行优化的过程如下：S1：首先确定当前方案待优化的参数个数D，在粒子群搜索空间内，随机生成N个包含待优化系数的D维粒子作为初代粒子群；
S2：选择摩擦速度均方根误差的相反数作为适应度函数，当均方根误差越小时，适应度越大，粒子群向适应度大的解的方向移动，适应度函数计算如下所示：；
S3：计算第t代粒子群中粒子的适应度函数值 ， ，…， ；
S4：依据第t代粒子的适应度函数值确定粒子飞行速度和粒子移动方向，第t代粒子群中第i个粒子的位置 ，第i个粒子在第j个维度上的飞行速度；
S5：将第i个粒子以如下公式的进化方法从第t代向第 代进化；
其中，t为迭代次数，为惯性权重， 为个体学习因子、为群体学习因子，取值范围[0，
2]，和 为[0,1]内的均匀随机数； 为第i个粒子在前t代进化过程中的最优适应值对应的位置， 为第1代到第t代粒子群中适应度最高的粒子的位置， 代表第i个粒子在第j维度上第 代的速度， 代表第i个粒子在第j维度上的第 代位置；
S6：重复步骤S3‑S5直到进化次数G=100，得到第i个粒子在G=100代中搜索到的最优位置 为 ，群 体 在 G = 1 0 0 代 中 搜 索 到 的 最 优 位 置 为，其中粒子群更新100代后的最优位置 对应的系数a和b即为最优系数。
7.根据权利要求1所述的一种海面空气动力学粗糙度的计算方法，其特征在于，步骤四中，使用遗传算法进行优化的过程如下：步骤1，依据经验设置交叉概率、变异概率和迭代次数，依据待优化的海面空气动力学粗糙度参数化方案确定2个待优化系数a和b、维度为D=2，设置初始种群大小N=30，给定a的精度0.1、b的精度0.001，确定个体长度L，给定参数a、b数值解范围a [ ， ]，b [ ， ]；
步骤2，随机生成二进制形式的包含N个体的初代种群 ， ，…， ；
计算初始种群中个体的二进制编码对应的十进制在解区间对应的解 ，计算公式如下所示：；
其中， 为个体的二进制值对应的十进制值，A、B为参数a、b的边界，对于参数a，A=A1，B=B1，对于参数b，A=A2，B=B2；
步骤3，对初代种群中每个个体对应的二进制形式进行交叉运算，即对两个个体二进制数的第9 13位进行交换；
对初代种群中每个个体对应的二进制形式进行变异运算，即对每一个个体的最后一位进行取反；
使用步骤2的公式计算交叉变异后的种群中每个个体对应的十进制解；
步骤4，使用实测10m平均风速数据、波高数据和波周期数据，在第一组的数据组合区间内，依据海面空气动力学粗糙度参数化方案结合COARE3.0模型计算该组合区间内所有样本计算的摩擦速度 ；
在第一组数据组合区间内，使用实测10m三维风速脉动量，结合涡动相关法计算摩擦速度 ；
在第一组数据组合区间内，依据区间内所有样本数据 ，计算第二代种群中每个个体的适应度值 ， ，…， ，第s代的适应度计算函数 如下所示：；
其中，s表示代数 ，为该区间内所有实测数据样本的总数，适应度函数选择摩擦速度的均方根误差的相反数，当均方根误差越小时，适应度越大，其中 为第m个样本使用涡动相关法计算的摩擦速度， 为第m个样本使用海面空气动力学粗糙度参数化方案计算的摩擦速度；
步骤5，使用轮盘赌方法依据第二代每个个体的适应度值与总适应度值之比 ， k=1，2，…，N；将一个轮盘划分成多个扇区，通过旋转轮盘，随机选择个体作为下一代，每个个体被选择的概率与其适应度成正比，第s代的适应度值与总适应度值之比 计算如下所示：；
轮盘赌方法具体如下：将轮盘展开成区间，区间表示为每个个体适应度与所有个体适应度总和之比的求和，在区间[0， ]内随机产生N个数，对随机数所在区间对应的个体进行选择运算，组成新一代种群；
判断迭代次数G是否达到100，若未达到，继续进行迭代，若达到迭代次数，则输出最优解a、b即为最优系数。
8.根据权利要求1所述的一种海面空气动力学粗糙度的计算方法，其特征在于，还包括步骤五，实验验证：依据实测数据，对比分析使用涡动相关法计算摩擦速度 和步骤四得到的计算粗糙度的方法结合COARE3.0模型计算的摩擦速度 ，验证步骤四得到的计算粗糙度的方法是否满足实际应用需求。
一种海面空气动力学粗糙度的计算方法
技术领域
[0001] 本发明涉及海洋大气边界层领域，特别涉及一种海面空气动力学粗糙度的计算方法。
背景技术
[0002] 海面空气动力学粗糙度的概念是从陆面上对数风廓线理论中的粗糙度延伸应用到海面而得来的，定义为风速等于零的高度，通常用 来表示。海面空气动力学粗糙度是描述海面小尺度粗糙程度的一个量，它的大小体现了海面微尺度起伏的程度，变化规律在一定程度上描述了海洋与大气之间进行动量传输的主要特征。
[0003] 目前，海面空气动力学粗糙度的研究主要依据海面风速、波龄、风浪的谱能量等相关参数和海面粗糙度之间的关系建立粗糙度参数化方案，然后借助于微波散射仪、卫星高度计和合成孔径雷达等设备测量风场而得到的相关数据验证糙度参数化方案的准确性。海面空气动力学粗糙度获取的准确性程度，对海上电磁波传播路径诊断、通信链路信道分析、蒸发波导诊断、雷达目标探测和遥感信息获取等研究的有重要的影响，因此，对于海面空气动力学粗糙度的研究显得尤为重要。
[0004] 由于海面空气动力学粗糙度空间分布不均，导致大范围测量比较困难，所以传统的海洋调查无法实时、实地测量海面空气动力学粗糙度，而目前的研究是基于量纲分析的海面空气动力学粗糙度参数化方案，不同的参数化方案是基于不同海域的气象水文参数数据建立的经验方案，而经验方案的确定高度依赖其研究海域的气象水文条件，因此，在实际应用计算时有必要对不同范围内的气象水文参数数据选择性地使用不同的海面空气动力学粗糙度参数化方案，并且需要依据该研究海域的气象水文参数数据进一步优化改进其经验方案，得到适用于该海域的海面空气动力学粗糙度计算方法。
发明内容
[0005] 为解决上述技术问题，本发明提供了一种海面空气动力学粗糙度的计算方法，依据海上实测数据结合涡动相关法建立在不同参数组合区间内使用拥有最优系数的海面空气动力学粗糙度参数化方案，达到提高参数化方案在该海域适用性的目的。
[0006] 为达到上述目的，本发明的技术方案如下：
[0007] 一种海面空气动力学粗糙度的计算方法，包括如下步骤：
[0008] 步骤一，数据采集：
[0009] 在海上观测平台安装风速风向仪、三维超声风、毫米波雷达获取海平面10m高度处的平均风速 、海平面10m高度处的三维风速( 、、)、波高 和波周期 参数数据；
[0010] 步骤二，敏感性分析：
[0011] 依据海平面10m高度处的平均风速 、波高 和波周期 参数数据的取值范围，设置初始间隔、划分初始参数区间并组合，建立初始参数组合区间；在初始参数组合区间内采用蒙特卡洛采样法，对每个组合区间内的参数进行采样，获取所需参数样本；在不同组合区间内，使用敏感性分析方法利用参数样本计算海面空气动力学粗糙度参数化方案的一阶敏感性指数和全局敏感性指数，得到各种海面空气动力学粗糙度参数化方案在不同参数组合区间内，对于三种参数的敏感性；依据敏感性结果，拆分或合并初始组合区间，重新建立新的参数组合区间；
[0012] 步骤三，适用性分析：
[0013] 在新的参数组合区间内，依据实测海平面10m高度处的三维风速( 、、)，以涡动相关法计算的摩擦速度 和不同海面空气动力学粗糙度参数化方案结合CAORE3.0模型计算的摩擦速度 为依据，对比不同方案的均方根误差 ，分析新的参数组合区间内各海面空气动力学粗糙度参数化方案的适用性，得到当前参数组合区间内适用性最高的海面空气动力学粗糙度参数化方案；
[0014] 步骤四，系数优化与新方法建立：
[0015] 使用智能优化算法在不同的参数组合区间内，依据适应度函数值逐步迭代，寻找当前参数组合区间内适用性最高的海面空气动力学粗糙度参数化方案的最优系数，得到一种在不同参数组合区间内使用拥有最优系数的不同海面空气动力学粗糙度参数化方案计算粗糙度的方法。
[0016] 上述方案中，海面空气动力学粗糙度参数化方案包括YT96、TY01、O02、GW06、PS07；敏感性分析方法包括Sobol指数法和扩展傅里叶振幅灵敏度检验法。
[0017] 上述方案中，敏感性参数包含海平面10m高度处的平均风速 、波高 和波周期，参数取值范围 、 、 ，单位为m/s、m、s，间隔设置为1。
[0018] 上述方案中，步骤二中，依据各模型某个输入参数在同一区间的最大全局敏感性指数作为区间划分的判断依据，将参数敏感性指数较大的参数区间进行拆分，将敏感性指数较小的区间进行组合。
[0019] 上述方案中，涡动相关法计算如下所示：
[0020] ；
[0021] 其中， 为纬向风速脉动量， 为经向风速脉动量， 为垂直方向风速脉动量，表示纬向和经向脉动量在20秒内的平均值的乘积；脉动量计算公式如下所示：
[0022] ；
[0023] ；
[0024] ；
[0025] 其中， 为纬向风速， 为经向风速， 为垂直方向风速， 、 、 为30分钟内的平均值。
[0026] 上述方案中，步骤三中，在不同数据区间内以 为标准，计算五种海面空气动力学粗糙度参数化方案的摩擦速度和涡动相关法计算的摩擦速度之间的均方根误差、 、 、 ，每种海面空气动力学粗糙度参数化方案单独计算，均方根误差的计算方法如下所示：
[0027] ；
[0028] 其中，n为样本总数， 为第m个样本使用涡动相关法计算的摩擦速度，为第m个样本使用五种海面空气动力学粗糙度参数化方案中每个方案单独计算的摩擦速度。
[0029] 上述方案中，步骤四中，智能优化算法包括粒子群算法、遗传算法、蚁群算法、模拟退火算法中的一种。
[0030] 上述方案中，步骤四中，使用粒子群算法进行优化的过程如下：
[0031] S1：首先确定当前方案待优化的参数个数D，在粒子群搜索空间内，随机生成N个包含待优化系数的D维粒子作为初代粒子群；
[0032] S2：选择摩擦速度均方根误差的相反数作为适应度函数，当均方根误差越小时，适应度越大，粒子群向适应度大的解的方向移动，适应度函数计算如下所示：
[0033] ；
[0034] S3：计算第t代粒子群中粒子的适应度函数值 ， ，…， ；
[0035] S4：依据第t代粒子的适应度函数值确定粒子飞行速度和粒子移动方向，第t代粒子群中第i个粒子的位置 ，第i个粒子在第j个维度上的飞行速度；
[0036] S5：将第i个粒子以如下公式的进化方法从第t代向第 代进化；
[0037] ；
[0038] ；
[0039] 其中，t为迭代次数，为惯性权重，为个体学习因子、为群体学习因子，取值范围[0，2]， 和 为[0,1]内的均匀随机数； 为第i个粒子在前t代进化过程中的最优适应值对应的位置， 为第1代到第t代粒子群中适应度最高的粒子的位置， 代表第i个粒子在第j维度上第 代的速度， 代表第i个粒子在第j维度上的第
代位置；
[0040] S6：重复步骤S3‑S5直到进化次数G=100，得到第i个粒子在G=100代中搜索到的最优位置为 ，群体在G=100代中搜索到的最优位置为，其中粒子群更新100代后的最优位置 对应的系数a和b
即为最优系数。
[0041] 上述方案中，步骤四中，使用遗传算法进行优化的过程如下：
[0042] 步骤1，依据经验设置交叉概率、变异概率和迭代次数，依据待优化的海面空气动力学粗糙度参数化方案确定2个待优化系数a和b、维度为D=2，设置初始种群大小N=30，给定a的精度0.1、b的精度0.001，确定个体长度L，给定参数a、b数值解范围a [ ， ]，b [ ，]；
[0043] 精度计算公式如下所示：
[0044] ；
[0045] ；
[0046] 步骤2，随机生成二进制形式的包含N个体的初代种群 ， ，…，；计算初始种群中个体的二进制编码对应的十进制在解区间对应的解 ，计算
公式如下所示：
[0047] ；
[0048] 其中， 为个体的二进制值对应的十进制值，A、B为边界；
[0049] 步骤3，对初代种群中每个个体对应的二进制形式进行交叉运算，即对两个个体二进制数的第9 13位进行交换；~
[0050] 对初代种群中每个个体对应的二进制形式进行变异运算，即对每一个个体的最后一位进行取反；
[0051] 使用步骤2的公式计算交叉变异后的种群中每个个体对应的十进制解；
[0052] 步骤4，使用实测10m平均风速数据、波高数据和波周期数据，在第一组的数据组合区间内，依据海面空气动力学粗糙度参数化方案结合COARE3.0模型计算该组合区间内所有样本计算的摩擦速度 ；
[0053] 在第一组数据组合区间内，使用实测10m三维风速脉动量，结合涡动相关法计算摩擦速度 ；
[0054] 在第一组数据组合区间内，依据区间内所有样本数据，计算第二代种群中每个个体的适应度值 ， ，…， ，第s代的适应度计算函数 如下所示：
[0055] ；
[0056] 其中，s表示代数 ，为该区间内所有实测数据样本的总数，适应度函数选择摩擦速度的均方根误差的相反数，当均方根误差越小时，适应度越大，其中 为第m个样本使用涡动相关法计算的摩擦速度， 为第m个样本使用海面空气动力学粗糙度参数化方案计算的摩擦速度；
[0057] 步骤5，使用轮盘赌方法依据第二代每个个体的适应度值与总适应度值之比， k=1，2，…，N；将一个轮盘划分成多个扇区，通过旋转轮盘，随机选择个体作为下一代，每个个体被选择的概率与其适应度成正比，第s代的适应度值与总适应度值之比计算如下所示：
[0058] ；
[0059] 轮盘赌方法具体如下：将轮盘展开成区间，区间表示为每个个体适应度与所有个体适应度总和之比的求和，在区间[0， ]内随机产生N个数，对随机数所在区间对应的个体进行选择运算，组成新一代种群；
[0060] 判断迭代次数G是否达到100，若未达到，继续进行迭代，若达到迭代次数，则输出最优解a、b即为最优系数。
[0061] 上述方案中，还包括步骤五，实验验证：依据实测数据，对比分析使用涡动相关法计算摩擦速度 和步骤四得到的计算粗糙度的方法结合COARE3.0模型计算的摩擦速度，验证步骤四得到的计算粗糙度的方法是否满足实际应用需求。
[0062] 通过上述技术方案，本发明提供的一种海面空气动力学粗糙度的计算方法具有如下有益效果：
[0063] 本发明使用敏感性分析方法，得到了不同海面空气动力学粗糙度参数化方案对不同气象水文参数的敏感性区间，为海面空气动力学粗糙度分区间计算方法指明了方向。
[0064] 本发明使用传统海面空气动力学粗糙度参数化方案与人工智能相关智能优化算法相结合的方法，改进了各传统海面空气动力学粗糙度参数化方案的系数，提高了各传统海面空气动力学粗糙度参数化方案在该研究海域适用性。
[0065] 本发明依托传统海面空气动力学粗糙度参数化方案，建立组合形式的粗糙度参数化方案，并且，本发明依据实测海上数据，对组合形式中的各海面空气动力学粗糙度参数化方案进行优化，得到该海域的最优系数，因此、与单一海面空气动力学粗糙度参数化方案相比，使用本发明优化后的海面空气动力学粗糙度计算方法在理论上更为合理，计算精度更高。
附图说明
[0066] 为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍。
[0067] 图1本发明实施例中所公开的一种海面空气动力学粗糙度的计算方法流程图。
[0068] 图2为10m平均风速、波高和波周期参数的初始组合区间示例图。
[0069] 图3为10m平均风速、波高和波周期新的参数组合区间示例图。
[0070] 图4为新的参数组合区间内粗糙度参数化方案适用性分析流程图。
[0071] 图5为新的参数组合区间内适用性最高的粗糙度参数化方案分布示例图。
[0072] 图6为粒子群算法优化流程图。
[0073] 图7为遗传算法优化流程图。
[0074] 图8为使用轮盘赌方法选择下一代种群的示意图。
[0075] 图9为粗糙度参数化方案优化结果示例图。
具体实施方式
[0076] 下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。
[0077] 本发明提供了一种海面空气动力学粗糙度的计算方法，如图1所示，包括如下步骤：
[0078] 步骤一，数据采集：
[0079] 在海上观测平台安装风速风向仪、三维超声风、毫米波雷达获取海平面10m高度处的平均风速 、海平面10m高度处的三维风速( 、、)、波高 和波周期 参数数据。
[0080] 海上观测平台包括但不限于海上气象梯度塔、舰载气象梯度观测仪、海上浮标等能够满足测量仪器安装的平台。
[0081] 使用风速风向仪测量10m高度处的平均风速 、使用三维超声风测量三维风速(、、)、使用毫米波雷达测量波高 和波周期 。
[0082] 步骤二，敏感性分析：
[0083] 依据海平面10m高度处的平均风速 、波高 和波周期 参数数据的取值范围，设置初始间隔、划分初始参数区间并组合，建立初始参数组合区间；在初始参数组合区间内采用蒙特卡洛采样法，对每个组合区间内的参数进行采样，获取所需参数样本；在不同组合区间内，使用敏感性分析方法利用参数样本计算海面空气动力学粗糙度参数化方案的一阶敏感性指数和全局敏感性指数，得到各种海面空气动力学粗糙度参数化方案在不同参数组合区间内，对于三种参数的敏感性；依据敏感性结果，拆分或合并初始组合区间，重新建立新的参数组合区间。
[0084] 海面空气动力学粗糙度参数化方案包括YT96、TY01、O02、GW06、PS07；敏感性分析方法包括Sobol指数法和扩展傅里叶振幅灵敏度检验法。
[0085] 敏感性参数包含海平面10m高度处的平均风速 、波高 和波周期 ，参数取值范围 、 、 ，单位为m/s、m、s，间隔设置为1。
[0086] 依据参数取值范围和间隔划分并组合参数区间，建立初始参数组合区间。
[0087] 依据各模型某个输入参数在同一区间的最大全局敏感性指数作为区间划分的判断依据，将参数敏感性指数较大的参数区间进行拆分，将敏感性指数较小的区间进行合并。
[0088] 如表1所示，本发明以五种粗糙度参数化方案YT96、TY01、O02、PS07、GW06举例说明。
[0089] 表1 五种海面空气动力学粗糙度参数化方案：
[0090]
[0091] 10m平均风速 、波高 和波周期 参数的取值范围和初始间隔设置如表2所示。
[0092] 表2 10m平均风速、波高和波周期参数范围和初始间隔：
[0093]
[0094] 如图2所示，依据表2中参数取值范围和间隔划分区间并组合，得到的3750组初始组合区间。在3750组初始组合区间内使用蒙特卡洛采样法，在每个初始组合区间中，分别获取所需输入参数样本。以第1个初始组合区间 、 、 为例，在该组合区间中获取输入参数样本（10m平均风速、波高和波周期）n组，样本参数个数为E=3，则生成n*2E的样本矩阵M如下：
[0095] ；
[0096] 其中，第一列 和第四列 代表10m风速数据，第二列 和第五列 代表波高数据，第三列 和第六列 代表波周期数据，将矩阵的前E列设置为矩阵A，后E列设置为矩阵B，则矩阵A、B如下：
[0097] ；
[0098] 对于矩阵A，使用矩阵B中的第i列替换其第i列，得到 （i=1，2，…，E），依此建立n*5组输入样本数据，则 、 和 如下：
[0099] ；
[0100] 以TY01方案为例，将n*5组输入参数样本数据代入TY01粗糙度参数化方案公式中结合COARE3.0模型计算粗糙度，得到n*5组TY01方案的粗糙度数据，将n*5组输入参数样本数据和n*5 数据结合构建n*5组样本数据。
[0101] 敏感性分析以Sobol指数法为例，依据n*5组样本数据计算每个初始区间内参数的一阶响应指数 、 、 和全局响应指数 、 、 ，计算方法如式（1）、（2）：
[0102] （1）
[0103] （2）
[0104] 其中，i=1，2，…，E。
[0105] 同理计算其他初始组合区间内各粗糙度参数化方案的一阶响应指数和全局响应指数，不同粗糙度方案在不同参数组合区间内全局敏感性指数分布如表3所示。
[0106] 表 3 不同粗糙度方案在不同参数组合区间内全局敏感性指数分布表：
[0107]
[0108] 以各模型某个输入参数在同一区间的最大全局敏感性指数作为区间划分的判断依据，当参数在某部分组合区间内敏感性指数大于0.5时，将该部分组合区间进行密集拆分，得到更多组合区间，当参数在某部分组合区间内敏感性指数小于0.3时，则将数组合区间进行合并，得到一个组合区间，否则不变，以此方法建立新的参数组合区间。
[0109] 例如：若风速在10 11m/s、11 12m/s、12 13m/s…区间内的ST指数较大，则将区间~ ~ ~10 11m/s、11 12m/s…进行拆分，例如将区间10 11m/s拆分为区间10 10.5m/s和区间10.5~ ~ ~ ~ ~
11m/s两个区间，将区间11 12m/s拆分为区间11 11.5m/s和区间11.5 12m/s两个区间；若风~ ~ ~
速在0 1m/s区间和0 2m/s内ST指数较小，则将两个区间合并为一个区间0 2m/s，具体拆分~ ~ ~
合并方式以敏感性指数为依据。
[0110] 对拆分和合并后的区间重新进行组合的得到新的参数组合区间。组合方式与初始组合区间类似。如：风速在第1个区间，波高在第1个区间，波周期第1个区间为第一组组合区间；风速在第1个区间，波高在第1个区间，波周期第2个区间为第二组组合区间；依此类推，风速在第P个区间，波高在第M个区间，波周期第N个区间为第P*M*N组组合区间。新的组合区间示例如图3所示。
[0111] 步骤三，适用性分析：
[0112] 在新的参数组合区间内，依据实测海平面10m高度处的三维风速( 、、)，以涡动相关法计算的摩擦速度 和五种不同海面空气动力学粗糙度参数化方案结合CAORE3.0模型计算的摩擦速度 、 、 、 、 为依据，对比不同方案的
均方根误差 ，分析新的参数组合区间内各海面空气动力学粗糙度参数化方案的适用性，得到当前参数组合区间内适用性最高的海面空气动力学粗糙度参数化方案。适用性分析方法如图4所示。
[0113] 依据图3中新的组合区间规则，将实测10m平均风速、波高和波周期三种参数数据和实测10m三维风速数据划分到对应的区间中。
[0114] 在每个数据组合区间内，依据三种参数数据，使用五种粗糙度参数化方案结合COARE3.0模型，计算五种粗糙度参数化方案对应的摩擦速度 、 、 、 、。
[0115] 使用实测10m三维风速数据（其中 为纬向风速，为经向风速， 垂直方向风速）结合涡动相关法计算摩擦速度 。
[0116] 涡动相关法计算如式（3）所示：
[0117] （3）；
[0118] 其中， 为纬向风速脉动量， 为经向风速脉动量， 为垂直方向风速脉动量，表示纬向和经向脉动量在20秒内的平均值的乘积；
[0119] 脉动量计算公式如式（4）、（5）、（6）所示：
[0120] （4）；
[0121] （5）；
[0122] （6）；
[0123] 其中， 为纬向风速， 为经向风速， 为垂直方向风速， 、 、 为30分钟内的平均值。
[0124] 在不同数据区间内以 为标准，计算五种海面空气动力学粗糙度参数化方案的摩擦速度和涡动相关法计算的摩擦速度之间的均方根误差 、 、、 ，每种海面空气动力学粗糙度参数化方案单独计算，均方根误差的计算方法如式（7）所示：
[0125] （7）；
[0126] 其中，n为样本总数， 为第m个样本使用涡动相关法计算的摩擦速度，为第m个样本使用五种海面空气动力学粗糙度参数化方案中每个方案单独计算的摩擦速度。
[0127] 选择均方根误差最小的方案作为当前参数组合区间内适用性最高的海面空气动力学粗糙度参数化方案，如图5所示。
[0128] 步骤四，系数优化与新方法建立：
[0129] 使用智能优化算法在不同的参数组合区间内，依据适应度函数值逐步迭代，寻找当前参数组合区间内适用性最高的海面空气动力学粗糙度参数化方案的最优系数，得到一种在不同参数组合区间内使用拥有最优系数的不同海面空气动力学粗糙度参数化方案计算粗糙度的方法。
[0130] 智能优化算法包括但不限于粒子群算法、遗传算法、蚁群算法、模拟退火算法中的一种。
[0131] 以粒子群算法和遗传算法优化粗糙度参数化方案TY01为例，系数a和b为待优化的2个系数。待优化的粗糙度参数化方案TY01如表4所示。
[0132] 表4 待优化的粗糙度参数化方案TY01：
[0133]
[0134] 如图6所示，使用粒子群算法进行优化的过程如下：
[0135] S1：首先确定当前方案待优化的最大迭代次数G=100，待优化系数个数D=2，设置位置限制a∈[0,2000]、b∈[0,10]，设置速度限制 。在粒子群搜索空间内，随机生成N=30个D=2维包含参数a、b的粒子 作为初代种群，第i个粒子 是一个二维的向量形式如 ，i=1，2，…，N。
[0136] S2：选择摩擦速度均方根误差的相反数作为适应度函数，当均方根误差越小时，适应度越大，粒子群向适应度大的解的方向移动，适应度函数计算如公式（8）所示：
[0137] （8）；
[0138] S3：计算第t代粒子群中粒子的适应度函数值 ， ，…， ；
[0139] S4：依据第t代粒子的适应度函数值确定粒子飞行速度和粒子移动方向，第t代粒子群中第i个粒子的位置 ，第i个粒子在第j个维度上的飞行速度；
[0140] S5：将第i个粒子以公式（9）、（10）的进化方法从第t代向第 代进化；
[0141] （9）；
[0142] （10）
[0143] 其中，t为迭代次数，为惯性权重，为个体学习因子、为群体学习因子，取值范围[0，2]， 和 为[0,1]内的均匀随机数； 为第i个粒子在前t代进化过程中的最优适应值对应的位置， 为第1代到第t代粒子群中适应度最高的粒子的位置， 代表第i个粒子在第j维度上第 代的速度， 代表第i个粒子在第j维度上的第
代位置；
[0144] S6：重复步骤S3‑S5直到进化次数G=100，得到第i个粒子在G=100代中搜索到的最优位置为 ，群体在G=100代中搜索到的最优位置为，其中粒子群更新100代后的最优位置 对应的系数a和b
即为最优系数。
[0145] 如图7所示，使用遗传算法进行优化的过程如下：
[0146] 步骤1，依据经验设置交叉概率、变异概率和迭代次数，依据待优化的海面空气动力学粗糙度参数化方案确定2个待优化系数a和b、维度为D=2，设置初始种群大小N=30，给定a的精度0.1、b的精度0.001，确定个体长度L=20（二进制数的位数），给定参数a、b数值解范围a [0,2000]，b [0,10]；
[0147] 精度计算公式如（11）、（12）所示：
[0148] （11）；
[0149] （12）；
[0150] 步骤2，随机生成二进制形式的包含N=30个体的初代种群 ， ，…，，的二进制编码形式示例如下所示：
[0151] （10111010101011110111，11101010111011110101）；
[0152] （11101101010101110110，10101101011001111101）；
[0153] ……
[0154] （10001101110001110101，10101101110001110100）；
[0155] 计算初始种群中个体的二进制编码对应的十进制在解区间对应的解 ，计算公式如（13）所示：
[0156] （13）；
[0157] 其中， 为个体的二进制值对应的十进制值，A、B为边界；则初始种群的解为：， ，…， 。
[0158] 步骤3，对初代种群中每个个体对应的二进制形式进行交叉运算，即对两个个体二进制数的第9 13位进行交换，如 和 交叉得到 和 ；~
、 、 、 如下所示：
[0159] （10111010101011110111，11101010111011110101）；
[0160] （11101101010101110110，10101101011001111101）；
[0161] （10111010010101110111，11101010011001110101）；
[0162] （11101101101011110110，10101101111011111101）；
[0163] 对初代种群中每个个体对应的二进制形式进行变异运算，即对每一个个体的最后一位进行取反；如 和 变异得到 和 ， 、 、、 如下所示：
[0164] （10111010101011110110，11101010111011110100）；
[0165] （11101101010101110111，10101101011001111100）；
[0166] （10111010010101110110，11101010011001110100）；
[0167] （11101101101011110111，10101101111011111100）；
[0168] 使用公式（13）计算交叉变异后的种群中每个个体对应的十进制解 ，，…， 。
[0169] 步骤4，使用实测10m平均风速数据、波高数据和波周期数据，在第一组的数据组合区间内，依据海面空气动力学粗糙度参数化方案TY01结合COARE3.0模型计算该组合区间内所有样本计算的摩擦速度 ；
[0170] 在第一组数据组合区间内，使用实测10m三维风速脉动量，结合涡动相关法计算摩擦速度 ；具体公式如式（3）。
[0171] 在第一组数据组合区间内，依据区间内所有样本数据，计算第二代种群中每个个体的适应度值 ， ，…， ，第s代的适应度计算函数 如公式（14）所示：
[0172] （14）；
[0173] 其中，s表示代数 ，为该区间内所有实测数据样本的总数，适应度函数选择摩擦速度的均方根误差的相反数，当均方根误差越小时，适应度越大，其中 为第m个样本使用涡动相关法计算的摩擦速度， 为第m个样本使用TY01海面空气动力学粗糙度参数化方案计算的摩擦速度；
[0174] 步骤5，使用轮盘赌方法依据第二代每个个体的适应度值与总适应度值之比， k=1，2，…，N；将一个轮盘划分成多个扇区，通过旋转轮盘，随机选择个体作为下一代，每个个体被选择的概率与其适应度成正比，第s代的适应度值与总适应度值之比计算如公式（15）所示：
[0175] （15）；
[0176] 轮盘赌方法具体实现如图8所示：将轮盘展开成区间，区间表示为每个个体适应度与所有个体适应度总和之比的求和，在区间[0， ]内随机产生N个数，对随机数所在区间对应的个体进行选择运算，组成新一代种群。依据该方法既能保留适应度较高的个体，也能给与适应度较低的个体一定的生存机会。
[0177] 判断迭代次数G是否达到100，若未达到，继续进行迭代，若达到迭代次数，则输出最优解a、b即为最优系数。
[0178] 同理，在第二组数据区间至第P*M*N组区间内依据上述方法可以得到每个区间内的最优系数。
[0179] 对于不同的组合区间内，适用性最高的粗糙度参数化方案，使用不同的智能优化算法，优化方案后得到每个方案的最优系数，从而得到每个区间内的最优粗糙度参数化方案，在不同参数组合区间的使用拥有最优系数的不同海面空气动力学粗糙度参数化方案，即为本发明建立的新方法。不同算法的优化结果如图9所示，其中a、b、c、d…表示不同粗糙度参数化方案最优系数，TY01best，PS07best，GW06best，…，O02best表示不同组合区间内拥有最优系数的海面空气动力学粗糙度参数化方案。
[0180] 步骤五，实验验证：
[0181] 依据实测海域海平面10m高度处的三维风速( 、、)，使用涡动相关法计算摩擦速度 ，依据海平面10m高度处的平均风速数据、波高和波周期数据使用每个参数组合区间最优海面空气动力学粗糙度参数化方案结合COARE3.0模型计算摩擦速度 ；依据和 使用公式（7）计算 ，依据 值验证优化后的海面空气动力学粗糙度参数化方案是否满足实际应用需求。
[0182] 对于实施例中，依据实测数据，使用本发明提出的粗糙度参数化方案优化方法在不同参数组合区间，使用不同的海面空气动力学粗糙度参数化方案和其最优系数，理论上与单一的海面空气动力学粗糙度参数化方案相比，更加合理。
[0183] 对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。
