结合统计模型与动力模型的乘潮水位预报方法

摘要

结合统计模型与动力模型的乘潮水位预报方法。该方法包括：获得该地潮汐调和常数，预报该地的预测天文潮位和增水水位，得出基于统计模型的预测乘潮水位Ls；根据该大气预报模型计算并得出预测的风场结果与潮汐的开边界条件，以此作为强迫场，使用海洋动力学模型计算得出包含该地天文潮位与增水水位的预测乘潮水位Lm；基于对两个模型在各种时间尺度和气象条件下预测结果的准确度分析，构造最佳的乘潮水位预报公式，给出最终乘潮水位L的预报结果，并通过可视化技术将该预报结果实时的展示出来。该方法有机的耦合气象与潮位的监测数据、潮位值的统计预测方法、大气与海洋动力学模型以及实时的可视化系统，使其达到准确和快速预报乘潮水位的要求。

说明书

技术领域

本发明涉及海洋科学与海洋工程领域，具体地说，是海港或海上乘潮水位的预报方法。

背景技术

我国自改革开放之后，进出口贸易快速增加。特别是21世纪以来，船舶的进出港压力愈加凸显。另一方面，由于全球气候变化所带来的具体天气异常情况越来越多，港口和航道区域也不例外，其中由于气象与外海信号等因素造成的增减水为港口和航道的通航安全提出了严峻的挑战。乘潮水位是指在一定时间间隔且可用以通航的某一高潮水位，在海港工程的设计以及航道的选取时常常要求参考此值。虽然这个值多是出现在高潮前后，但因为诸多如天气等不确定的因素，科学地掌握它仍然需要进行若干复杂的计算和经验，否则就可能因为通航水位设计过浅而影响通航能力，或者由于疏浚过深而带来较大的浪费。更重要的是，对乘潮水位掌握不准会影响船舶通航的调度，并在异常情况下带来诸多的不安全因素。因此，准确的掌握乘潮水位在船舶运输、港口工程等诸多方面有着重要的实用价值，尤其是在水深较浅和潮差较大的港口与航道，其作用更加显而易见。但目前由于气候异常事件日渐增多，传统的乘潮水位估计方法的准确性和及时性已无法满足当前对乘潮水位预报的需求。主要原因在于一方面经验方法所基于的物理基础较为薄弱，难以准确了解大区域天气和海洋环境的变化对当地乘潮水位的影响，而另一方面，纯动力学数值模型的计算受当今动力学基础研究的限制，在短时间内尚不能替代经验公式估算的准确性。因此，有机结合经验统计方法与动力学模型计算方法的优势，在各种时间尺度上给出乘潮水位更为准确和及时的预报，不仅可以为保护我国近海港口的运行效率提供技术支持，同时也将促进我国航运事业的安全发展。

发明内容

本发明的目的是提供一种在各种时间尺度上准确和及时地预报乘潮水位的方法。

本发明提供的乘潮水位预报方法，有机的耦合气象与潮位的监测数据、潮位值的统计预测方法、大气与海洋动力学模型以及实时的可视化系统，使其达到准确和快速预报乘潮水位的要求，其具体方法和步骤如下(参见图1)：

第1、整理历史气象观测数据与历史潮位观测数据，以及当地的地形数据，分别建立历史气象数据库、历史潮位数据库以及当地的地形数据库；

第2、通过第1步建立的历史潮位数据库中的历史观测潮位数据获得该地潮汐调和常数，并通过该调和常数进行调和分析并预报出该地的预测天文潮位；

第3、实时获取气象监测数据与潮位监测数据，使用增减水的统计模型预测出该地的增水水位；

第4、加和第2步中预测的天文潮位与第3步中预测的增水水位，得出基于统计模型的预测乘潮水位L_s；

第5、使用数据同化方法同化第3步中的实时气象监测数据至大气预报模型(如WRF)，根据该大气预报模型计算并得出预测的风场结果与潮汐的开边界条件；

第6、由第5步中的大气预报模型预测的风场与潮汐的开边界条件作为强迫场，使用海洋动力学模型计算得出包含该地天文潮位与其他因素增水水位的预测乘潮水位L_m；

第7、根据第1步建立的历史潮位数据库中的历史观测潮位数据序列，分别使用统计模型与海洋动力学模型计算出相应时间段的预测潮位时间序列，并基于对两个模型在各种时间尺度和气象条件下预测结果的准确度分析，构造最佳的乘潮水位预报公式；

第8、结合第4步中统计模型预测的乘潮水位L_s与第6步中海洋动力学模型预测的乘潮水位L_m结果，使用第7步中构造的乘潮水位预报公式给出最终乘潮水位L的预报结果，并通过可视化技术将该预报结果实时的展示出来。

其中，第2步所述的潮汐调和常数的计算方法是，根据第1步中建立的历史潮位数据库中历史观测数据的时间长度，使用最小二乘法分析出该地的调和常数，其中调和常数的振幅与迟角分别用H与g表示；

第2步所述的天文潮位是指由引潮力直接或间接引起的潮位，基于以上得出的该地调和常数，根据以下表达式计算预测天文潮位：

$\zeta \left(t\right)=a_0+\underset{j=1}{\overset{m}{\Sigma }}{f}_j{H}_j\cos [{\sigma }_{j}t+{(v_0+u)}_j-g_j]$

式中a₀为平均海面在潮高基准面上的高度，是适用于航海使用的起算面；根据需要也可取平均海面在验潮站水尺零点上的高度；j为分潮的个数；f_j、H_j、g_j、σ_j和u_j分别表示分潮交点因数、分潮振幅、分潮迟角、分潮角速率和分潮交点订正角；f_j、H_j、g_j、σ_j和(V₀+u)_j分别对于不同分潮有不同的量值；H_j、g_j对于固定地点是常数；f_j、(V₀+u)_j是时间的函数；f_j、u的变化周期为18.61年，因此，可取每月或每年中间日期的时刻作计算，视为恒量；V₀为所取日期零时的量值，也是恒量。

第3步所述的增水水位的计算方法是，使用拉格朗日插值或者三次样条函数方法，基于临近时间的潮位曲线，计算出该地的增水水位ζ’，并得出基于统计模型的预测乘潮水位L_s＝ζ+ζ’。

第5步所述的数据同化方法以及预测风场结果与潮汐的开边界条件的方法是：

首先对气象观测数据进行格式转换与质量控制，然后利用3DVAR方法进行数据同化，为预报所使用的海洋动力学模型提供初始场与时变边界条件；采用并配置WRF-ARW大气动力学模型为所关注海区提供风场以及相关气象条件未来48小时的预报；该模型采用完全可压缩非静力欧拉方程组，水平网格采用Arakawa C网格，垂直坐标采用目前国际上广泛应用的基于质量的地形追随η坐标。

第6步所述的使用海洋动力学模型计算得出包含该地天文潮位与增水水位的预测乘潮水位L_m的方法是：

基于第5步计算所得的风场结果以及来自于NAO99全球潮汐模型(http://www.miz.nao.ac.jp/staffs/nao99/index_En.html)所提供的8个分潮M2、S2、N2、K2、K1、O1、P1、Q1的调和常数作为FVCOM海洋动力学模型的强迫条件与潮汐的开边界条件，所述的海洋动力学模型为FVCOM(An Unstructured Grid，Finite-Volume Coastal Ocean Model)海洋动力学模型，简称FVCOM模型，该模型是一个自由海表面、流体静力学和Boussinesq近似的、原始方程组的海洋模型；

在第5步中所述WRF-ARW模型运算完毕后，SHELL脚本自动运行FVCOM模型，并计算获取该地未来48小时的预测乘潮水位L_m。

第7步所述构造最佳乘潮水位L的预报公式的方法是：

根据第一步建立的数据库中针对当地的具体历史观测潮位数据，分别使用上述统计模型方法与海洋动力学模型方法进行预测乘潮水位并计算二者的平均误差，这里我们定义：

a＝|E_m|/(|E_m|+|E_s|)

其中E_m为海洋动力学模型预测乘潮水位L_m的平均误差，而E_s为统计模型预测乘潮水位L_s的误差；同时，我们发现统计模型在5个小时内的预测结果相对海洋动力学模型的预测结果更为准确，且时间越靠近预报时间准确度越高，故我们初步构造当地的乘潮水位L的预报公式如下：

$L=\frac{(1-a)t^2+25a}{25+t^2}(a\frac{25}{25+t^2}{L}_s+(1-a)\frac{t^2}{25+t^2}{L}_m)$

其中，t为预测提前的时间。

这里需要说明的是，系数a可以对不同风速与气压条件下分别统计，并在实际应用时对应相应的气象条件来选用系数a。

第8步所述通过可视化技术将该预报结果实时的展示出来的方法包括：

使用绘图技术实时的将48小时乘潮水位L的预测结果，连同历史的观测与预测潮位曲线显示出来，并使用web技术在互联网上发布，同时通过相关通讯手段第一时间将信息提供于港口调度人员与海上船舶相关负责人，为港口与航运安全、高效的进行，以及异常潮位变化时的决策提供有力的技术支持。

本发明提供的上述预报方法路线的计算与可视化过程均由计算机自动调度执行。

本发明所涉及技术术语的含义是：

潮位：文中所述潮位皆指总的潮位。

天文潮位：由引潮力直接或间接引起的潮位，在本文中称为天文潮位。

调和常数：将实测潮位资料分解出许多分潮，所求出每个分潮的平均振幅和迟角值。

增水水位：由于气象因素与外海非潮流海流传入信号影响所造成的水位变化值。

乘潮水位：船舶在通过航道(包括进港航道)的局部浅段时，由于水深不足，常利用一定的高潮位以增加航深使船舶通过。这种使船舶能在一定时间内，乘一定的较大潮位通过航道浅段的水位称为乘潮水位。

最小二乘法：最小二乘法(又称最小平方法)是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。

本发明的优点和有益效果：

本发明方法有机的结合了乘潮水位预报中统计模型与海洋动力学模型的优点，极大的提高了乘潮水位预测的准确性，本发明使用计算机脚本自动完成乘潮水位预报模型的运算过程，并将预报结果进行实时的可视化展示，为港口与海上航运安全提供了有力的技术保障。

附图说明

图1是乘潮水位预报方法的技术路线图。

图2是使用本发明方法对天津港2010年1月1日逐时水位的预报结果与实际观测结果的对比图。

具体实施方式

如图1所示，本发明给出了乘潮水位预报技术生成的详细步骤，同时为方便描述和理解，对这一技术方法的构造过程实例描述如下：

(1)建立该地历史观测的气象数据与潮位数据库，并使用Linux/UNIX-SHELL脚本调用lftp、wget等下载工具定时自动收集由气象监测站与潮位观测站实时传回的数据，且添加这些数据至已有数据库中。

(2)根据潮位历史观测资料的时间长度，使用最小二乘法分析出该地的调和常数。

(3)基于该调和常数，根据以下表达式计算预测天文潮位：

$\zeta \left(t\right)=a_0+\underset{j=1}{\overset{m}{\Sigma }}{f}_j{H}_j\cos [{\sigma }_{j}t+{(v_0+u)}_j-g_j]$

式中a₀为平均海面在潮高基准面上的高度，是适用于航海使用的起算面。根据需要也可取平均海面在验潮站水尺零点上的高度；j为分潮的个数；f_j、H_j、g_j、σ_j和(V₀+u)_j分别对于不同分潮有不同的量值。H_j、g_j对于固定地点是常数。f_j、(V₀+u)_j是时间的函数。f_j、u的变化周期为18.61年，因此，可取每月或每年中间日期的时刻作计算，视为恒量。V₀为所取日期零时的量值，也是恒量。

(4)使用拉格朗日插值或者三次样条函数方法，基于临近时间的潮位曲线，计算出该地的增水水位ζ’，并得出统计模型的预测乘潮水位L_s＝ζ+ζ’；

(5)首先对气象观测数据进行格式转换与质量控制，然后利用3DVAR方法进行数据同化，为预报所使用的动力学模型提供初始场与时变边界条件。采用并配置WRF-ARW大气动力学模型为所关注海区提供风场以及相关气象条件未来48小时的预报。该模型采用完全可压缩非静力欧拉方程组，水平网格采用Arakawa C网格，垂直坐标采用目前国际上广泛应用的基于质量的地形追随η坐标，η层可根据需要改变。

(6)FVCOM(An Unstructured Grid，Finite-Volume Coastal Ocean Model)是一个自由海表面、流体静力学和Boussinesq近似的、原始方程组的海洋模型。基于以上步骤计算所得的风场以及来自于NAO99全球潮汐模型(http://www.miz.nao.ac.jp/staffs/nao99/index_En.html)所提供的8个分潮(M2、S2、N2、K2、K1、O1、P1、Q1)的调和常数作为FVCOM海洋模型的的强迫条件与开边界条件，在(5)中所述WRF-ARW模型运算完毕后，SHELL脚本自动运行FVCOM模型，并计算获取该地未来48小时的预测乘潮水位L_m。

(7)针对某地的具体历史观测数据，我们分别使用上述统计方法与动力学数值模型方法进行预测并计算二者的平均误差，这里我们定义：

a＝|E_m|/(|E_m|+|E_s|)

其中E_m为动力学模型预测乘潮水位的平均误差，而E_s为统计模型预测乘潮水位的误差。同时，我们发现统计模型在5个小时内的预测相对动力学模型更为准确，且时间越靠近准确度越高，故我们初步构造该地的乘潮水位L预报公式如下：

$L=\frac{(1-a)t^2+25a}{25+t^2}(a\frac{25}{25+t^2}{L}_s+(1-a)\frac{t^2}{25+t^2}{L}_m)$

其中t为预测提前的时间。这里需要说明的是，系数a可以对不同风速与气压条件下分别统计，并在实际应用时对应相应的气象条件来选用系数a。

(8)基于上步的计算结果L，使用绘图技术实时的将48小时潮位预测结果，连同历史的观测与预测潮位曲线显示出来，并使用web技术在互联网上发布，同时通过相关通讯手段第一时间将信息提供于港口调度人员，船长等相关负责人，为港口与航运安全、高效的进行，以及异常潮位变化时的决策提供有力的技术支持。

预报结果实例见图2。

图2为使用本发明方法对天津港2010年1月1日逐时水位的预报结果与实际观测结果的对比图。图中细实线为水位观测结果，粗实线为使用本方法的预测结果(具体数值见下表)，可以看出平均预报误差已经小于5％，达到了非常理想的水平。

表乘潮水位预报值(单位：厘米)

  时间  实测值  预测值  统计模型  动力模型  0  186  216  236  256  1  263  289  309  329  2  322  345  355  375  3  346  360  380  390  4  345  354  374  394  5  324  338  348  358  6  295  304  314  334  7  261  266  276  296  8  228  238  241  261  9  214  217  217  227  10  227  222  219  219  11  266  262  252  242  12  320  311  301  291

  13  367  365  345  325  14  397  386  366  346  15  399  391  361  341  16  380  372  332  312  17  344  334  288  268  18  293  283  233  213  19  237  224  174  154  20  188  168  120  108  21  155  143  83  63  22  156  148  70  60  23  196  186  96  66