技术领域
[0001] 本
发明涉及海洋学模型参数率定领域，具体涉及一种多个海洋学模型参数的同时率定方法。
背景技术
[0002] 海洋学模型通过求解Navier-Stokes控制方程，获得海洋环境变量的
时空分布，具有定量、预报性强的特色，在
基础学科研究和资源开发活动中起着越来越重要的作用。此外，精确的海洋数值预报结果是开展海上
溢油漂移预测和海上遇险目标漂移预测的基础。
[0003] 海洋学模型中，存在着众多需要参数化的过程。海洋上边界处的
风场是海洋
能量的重要来源，海面风场对于海洋的影响常用依赖风速的风应
力计算公式进行参数化，其中风
应力拖曳系数的选取对于风应力的计算至关重要。底部摩擦效应对于海洋能量耗散的影响，是一个经典的海洋学问题，底摩擦效应常用二次形式的依赖于流速的计算公式进行参数化，其
中底摩擦系数的取值很大程度上影响能量耗散的计算结果。可以看出，风应力拖曳系数和底摩擦系数的合理取值及相互匹配，是决定海洋数值模拟结果准确与否的重要因素。
[0004] 众多学者基于实验室实验或者海上观测实验，单独给出了风应力拖曳系数或者底摩擦系数的经验公式。一般情况下，经验公式是根据一定的物理意义给出的，但其中的经验系数需要根据实验数据进行率定。而经验系数在不同的实验条件、不同的海域、不同的海洋环境下，取值往往不同，因此，不假思索的直接采用文献中的经验系数，可能会造成数值模拟结果与实际情况相差较大。此外，在利用海洋学模型开展数值模拟时，传统的试错法很难实现对于多个模型参数进行高效、同步率定。
[0005] 在此背景下，如何利用实验数据或者观测数据，对研究海域的风应力拖曳系数和底摩擦系数经验公式中的经验系数同时进行自动化率定，是进一步提高海洋数值模拟
精度和预报精度的重要问题，也是本发明主要解决的技术问题。
发明内容
[0006] 有鉴于此，本发明提供了一种多个海洋学模型参数的同时率定方法，解决了上述问题中的至少一个。
[0007] 为了达到上述目的，本发明的技术方案是这样实现的：
[0008] 一种多个海洋学模型参数的同时率定方法，包括：
[0009] 构建地波雷达观测系统，获取表
层流场观测数据；
[0010] 获取研究海域
海水层结和水深，由表层流场推导出深度平均的流场；
[0011] 构建球
坐标系下考虑风应力和海表面气压的二维
正压海洋伴随同化模型；
[0012] 给定二维正压海洋伴随同化模型中风应力拖曳系数和底摩擦系数的经验公式，给出经验公式中经验系数调整策略；
[0013] 以二维正压海洋伴随同化模型为基础，同化深度平均的流场数据，根据给出的经验系数调整策略，同时对风应力拖曳系数和底摩擦系数经验公式中的经验系数进行优化调整，从而实现对于风应力拖曳系数和底摩擦系数的同时率定。
[0014] 可选地，获取研究海域
海水层结和水深，由表层流场推导出深度平均的流场；包括：
[0015] 根据研究海域海水层结和水深，求解垂向本征方程得到垂向本征函数；
[0016] 利用模态分解理论和最小二乘理论，由表层流场构造全水深的流场，进而得到深度平均的流场。
[0017] 可选地，构建球坐标系下考虑风应力和海表面气压的二维正压海洋伴随同化模型；包括：
[0018] 所述二维正压海洋伴随同化模型包含正向模
块、反向模块和梯度模块。
[0019] 可选地，所述正向模块的控制方程为：
[0020]
[0021]
[0022]
[0023] 其中，t代表时间，λ和φ分别是经度和纬度，z为深度，自静止海表面向上为正，u和v分别表示λ和φ方向上深度平均流速分量，p为
大气压强，ρ为
水体密度，a＝Rcosφ，R是地球半径，g是重力
加速度，Ω为地球自转
角速度，f＝2Ωsinφ是Coriolis参数，Ah是水平涡动粘性系数，ρa为空气密度，Cs为风应力拖曳系数，W为海面10米高度的风速，h为静水深，ζ为海表面起伏， 为引潮势，Cd为底摩擦系数，Δ是拉普拉斯算子。
[0024] 可选地，反向模块的控制方程是针对于同化流速数据而推导出来的，具体为：
[0025]
[0026]
[0027]
[0028] 其中，t代表时间，λ和φ分别是经度和纬度，z为深度，自静止海表面向上为正，u和v分别表示λ和φ方向上深度平均流速分量，p为大气压强，ρ为水体密度，a＝Rcosφ，R是地球半径，g是
重力加速度，Ω为地球自转角速度，f＝2Ωsinφ是Coriolis参数，Ah是水平涡动粘性系数，ρa为空气密度，Cs为风应力拖曳系数，W为海面10米高度的风速，h为静水深，ζ为海表面起伏，为引潮势，Cd为底摩擦系数，Ku和Kv是权重矩阵，理论上是观测误差矩阵的逆矩阵，u和 为深度平均流速u和v的观测值，ζa、ua和va分别是ζ、u和v的伴随变量，Δ是拉普拉斯算子。
[0029] 可选地，给定二维正压海洋伴随同化模型中风应力拖曳系数和底摩擦系数的经验公式，给出经验公式中经验系数调整策略；包括：
[0030] 给定的风应力拖曳系数的经验公式为：
[0031]
[0032] 其中，Cs为风应力拖曳系数，W为10米高度风速，C0为风速小于11米时的常数分布的风应力拖曳经验系数，k、b、A、B、C为经验系数。
[0033] 可选地，给定二维正压海洋伴随同化模型中风应力拖曳系数和底摩擦系数的经验公式，给出经验公式中经验系数调整策略；包括：
[0034] 给定的底摩擦系数的经验公式为：
[0035]
[0036] 其中，Cd为底摩擦系数，g为重力加速度，Cchezy为Chezy系数，n为曼宁经验系数，α为经验系数。
[0037] 可选地，给定二维正压海洋伴随同化模型中风应力拖曳系数和底摩擦系数的经验公式，给出经验公式中经验系数调整策略；包括：
[0038] 以所述梯度模块为基础，经验公式中经验系数调整策略如下：
[0039]
[0040] 其中，pc代表任一经验系数，pcold为优化调整前的pc值，pcnew为优化调整后的pc值，ss为调整步长， 为代价函数关于经验系数的梯度。
[0041] 可选地，所述代价函数关于经验系数的梯度是根据所述梯度模块和分步求导公式求得的。
[0042] 可选地，以二维正压海洋伴随同化模型为基础，同化深度平均的流场数据，根据给出的经验系数调整策略，同时对风应力拖曳系数和底摩擦系数经验公式中的经验系数进行优化调整，从而实现对于风应力拖曳系数和底摩擦系数的同时率定，包括：
[0043] 风应力拖曳系数和底摩擦系数的同时率定是指，对风应力拖曳系数和底摩擦系数经验公式中的经验系数进行同时率定，从而实现对具有一定物理意义的风应力拖曳系数和底摩擦系数的同时率定。
[0044] 与
现有技术相比，本发明的优点是：本发明利用伴随同化方法同化地波雷达观测数据，将参数率定问题转化为数值
迭代过程，能够自适应地同时率定多个海洋学模型参数；此外，通过引入风应力拖曳系数和底摩擦系数的经验公式，对公式中的经验系数进行率定，而不是直接对风应力拖曳系数和底摩擦系数进行率定，避免了传统伴随同化方法率定参数时不考虑物理意义的弊端。基于地波雷达观测数据的多个海洋学模型参数的同时率定，为进一步提高海洋数值预报精度奠定了坚实的基础。
附图说明
[0045] 构成本发明的一部分的附图用来提供对本发明的进一步理解，本发明的示意性
实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：
[0047] 图2是本发明对风应力拖曳系数和底摩擦系数经验公式中的经验系数进行优化调整的流程示意图。
具体实施方式
[0048] 下面结合附图对本发明的较佳实施例进行详细阐述，以使本发明的优点和特征能更易于被本领域技术人员理解，从而对本发明的保护范围做出更为清楚明确的界定。
[0049] 如图1所示，一种多个海洋学模型参数的同时率定方法，包括以下步骤：
[0050] S100：构建地波雷达观测系统，获取表层流场观测数据；
[0051] 具体地，地波雷达是一种新兴的海洋监测技术，具有超视距、大范围、全天候以及低成本等优点，能够从雷达回波中提取风场、浪场、流场等海况信息；本步骤中，可以构造多站地波雷达组网的观测系统，从而提高表层流场的观测精度。
[0052] S200：获取研究海域海水层结和水深，由表层流场推导出深度平均的流场。
[0053] 具体地，包括：
[0054] S201：根据研究海域海水层结和水深，求解垂向本征方程得到垂向本征函数。
[0055] 其中，垂向本征方程如下所述：
[0056]
[0057] 其中，ψ为垂向本征函数，z为垂向坐标(从海表面起，向上为正)，N为表征海水层结的浮
频率，h为水深，c为
波速。由公式(1)可知，垂向本征函数具有
正交性。
[0058] 利用WOA数据集或者HYCOM模式实时预报系统提供的
温度、
盐度数据，可以计算出研究海域的浮频率。由浅水波的频散关系可求得正压模态的波速；利用WKB近似，可以得到斜压模态的波速。将浮频率和波速带入垂向本征方程(1)，可以求解出不同模态的垂向本征函数。
[0059] S202：利用模态分解理论和最小二乘理论，由表层流场构造全水深的流场，进而得到深度平均的流场。
[0060] 其中，所述模态分解理论为：
[0061]
[0062] 其中，u和v分别为东西向和南北向流速，t为时间，um和vm是第m模态(m＝0为正压模态，m>0为斜压模态)的流速分量，ψm为第m模态的垂向本征函数(即式(1)的解)，Um和Vm分别是um和vm对应的时间变化系数。
[0063] 最小二乘理论是使得地波雷达观测系统所得表层流速观测值和利用前M个垂向模态构造的表层流速值之间的误差平方达到最小(以流速分量u为例)：
[0064]
[0065] 其中，z1代表表层，M为所选取的模态数(包含第0模态的正压模态和前M个斜压模态)。
[0066] 通过将已经求解得到的前M个垂向本征函数带入式(3)，可在最小二乘意义下，求得Um的值；进而利用模态分解公式(2)，可以得到全水深的流场。将全水深的流场进行垂向加权平均，可以得到垂向平均的流场。
[0067] S300：构建球坐标系下考虑风应力和海表面气压的二维正压海洋伴随同化模型；
[0068] 具体地，二维正压海洋伴随同化模型包含正向模块、反向模块和梯度模块。
[0069] 其中，正向模块的控制方程为：
[0070]
[0071]
[0072]
[0073] 其中，t代表时间，λ和φ分别是经度和纬度，z为深度，自静止海表面向上为正，u和v分别表示λ和φ方向上深度平均流速分量，p为大气压强，ρ为水体密度(1025kg·m-3)，a＝Rcosφ，R是地球半径，g是重力加速度，Ω为地球自转角速度，f＝2Ωsinφ是Coriolis参数，Ah是水平涡动粘性系数，ρa为空气密度(1.27kg·m-3)，Cs为风应力拖曳系数，W为海面10米高度的风速，h为静水深，ζ为海表面起伏， 为引潮势，Cd为底摩擦系数，Δ是拉普拉斯算子:
[0074]
[0075]
[0076] 正向模块的初始条件采用HYCOM模式实时预报系统提供的流场数据，上边界条件采用NCEP
气候预测系统提供的再分析风场、气压场，开边界条件采用OSU全球潮汐模型提供的M2、S2、K1、O1、N2、K2、P1和Q1八个主要分潮的水位数据，闭边界条件为法向流速为0，水深由ETOPO5数据插值得到。
[0077] 为构建反向模块和梯度模块，特构造表征深度平均流场模拟结果和观测结果之间误差的代价函数：
[0078]
[0079] 其中，J(u,v)为代价函数，Ku和Kv是权重矩阵，理论上是观测误差矩阵的逆矩阵，u和 为深度平均流速u和v的观测值，Σ是正向模块的时空计算区域。
[0080] 进而根据拉格朗日乘子法，构造拉格朗日函数：
[0081]
[0082] 其中，ζa、ua和va分别是ζ、u和v的伴随变量。
[0083] 令拉格朗日函数关于伴随变量ζa、ua和va的偏导数为0，可以得到正向模块的控制方程(式(4))。令拉格朗日函数关于模型变量ζ、u和v的偏导数为0，可以得到反向模块(伴随模型)的控制方程，具体如下：
[0084]
[0085]
[0086]
[0087] 反向模块的初始条件为伴随变量ζa、ua和va的值为0，侧边界条件为伴随变量ua和va的法向分量为0。
[0088] 令拉格朗日函数关于风应力拖曳系数和底摩擦系数的偏导数为0，可以求出梯度模块(代价函数关于模型参数的梯度表达式)，具体如下：
[0089]
[0090]
[0091] S400：给定二维正压海洋伴随同化模型中风应力拖曳系数和底摩擦系数的经验公式，给出经验公式中经验系数调整策略；
[0092] 本步骤中，所述风应力拖曳系数经验公式采用依赖于风速的分段函数形式，具体如下：
[0093]
[0094] 其中，Cs为风应力拖曳系数，W为10米高度风速，C0为风速小于11米时的常数分布的风应力拖曳经验系数，k、b、A、B、C为经验系数。上述经验系数C0、k、b、A、B、C均假定为常数、需要进行率定。
[0095] 所述底摩擦系数的经验公式采用依赖于水深的形式，具体如下：
[0096]
[0097] 其中，Cd为底摩擦系数，g为重力加速度，Cchezy为Chezy系数，n为曼宁经验系数，α为经验系数，假定上述经验系数均为常数且需要进行率定。
[0098] 本步骤中，根据梯度模块和分步求导公式，可知代价函数关于经验公式(式(10)和式(11))中经验系数的梯度表达式为：
[0099]
[0100]
[0101]
[0102]
[0103]
[0104]
[0105]
[0106]
[0107] 根据最速下降法，给出经验公式中经验系数调整策略，如下：
[0108]
[0109] 其中，pc代表任一经验系数，pcold为优化调整前的pc值，pcnew为优化调整后的pc值，ss为调整步长， 为代价函数关于经验系数的梯度(即：式(12))。
[0110] S500：以二维正压海洋伴随同化模型为基础，同化深度平均的流场数据，根据给出的经验系数调整策略，同时对风应力拖曳系数和底摩擦系数经验公式中的经验系数进行优化调整，从而实现对于风应力拖曳系数和底摩擦系数的同时率定。
[0111] 具体地，风应力拖曳系数和底摩擦系数的同时率定是指，不是对其中的一个参数进行单独率定，也不是直接假定参数是常数或者空间分布的形式直接进行率定，而是对风应力拖曳系数和底摩擦系数经验公式中的经验系数进行同时率定，从而实现对具有一定物理意义的风应力拖曳系数和底摩擦系数的同时率定。
[0112] 如图2所示，对风应力拖曳系数和底摩擦系数经验公式中的经验系数进行优化调整的流程可细分为如下步骤：
[0113] S501：根据参考文献，给定风应力拖曳系数经验公式(式(10))和底摩擦系数经验公式(式(11))中经验系数C0、k、b、A、B、C、n、α的初始猜测值；
[0114] S502：对正向模块(式(4))进行数值离散，在时间上正向积分正向模块，然后在时间上反向积分反向模块(式(8))，进而计算出步骤S300中代价函数值(式(6))和步骤S400中代价函数关于经验系数的梯度值(式(12))；
[0115] S503：根据步骤S400中经验系数调整策略(式(13))，同时对风应力拖曳系数和底摩擦系数经验公式中的经验系数进行优化调整；
[0116] S504：重复步骤S502、S503，直至相邻两次的代价函数的差小于1×10-4或者重复次数达到100次；将所得经验系数C0、k、b、A、B、C、n、α的值带入步骤S400中风应力拖曳系数经验公式(式10)和底摩擦系数经验公式(式(11))，即可得到最终基于地波雷达观测数据所率定出的风应力拖曳系数和底摩擦系数。
[0117] 以上所述仅为本发明的实施例，并非因此限制本发明的
专利范围，凡是利用本发明
说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其他相关的技术领域，均同理包括在本发明的专利保护范围内。
