CN110727923A - 一种短期潮位预报方法 - Google Patents

Abstract

本发明涉及一种短期潮位预报方法，该方法基于数据融合技术，利用扩展卡尔曼滤波(EKF)，将通用的潮汐模型计算水位融合到短序列潮位直接调和预报模型之中，不依靠额外高精度数据或计算关系，即可生成比两种源数据精度都更高的融合后预报数据，数据测验表明，对于3天的验潮数据，EKF方法至少对向后5天的预报有效，融合后数据较两种源数据的平均改善度分别为33％和60％；对于7天的验潮数据，EKF方法在向后20的预报中几乎都产生了改善效果，对应两种源数据的平均改善度分别为17％和38％。

Description

一种短期潮位预报方法

技术领域

本发明涉及潮位预报技术领域，特别是一种近岸海域的短期潮位预报方法。

背景技术

近岸海域的潮位观测、分析和预报，对主要集中于此区域的各类海洋工程施工、交通运输、渔业生产等经济活动有直接的影响。尤其是海洋测绘工作，其工作质量相当程度建立在潮位资料的精度之上。我国大部分地区的长期潮位观测起步较晚、总体数量少且在各地区的分布也不均匀。以江苏省为例，原海洋局系统进行长期观测并提供预报服务的，仅连云港、吕四等8个站位。除数量少之外，这些长期验潮站主要设置在海岸线附近，向海的控制距离有限。实际工作中，大量的岸外潮位主要还是通过基于不同的工作任务设立的短期临时验潮站获得的。

海上临时验潮站虽然数量多，但受工程工期、航次安排、仪器设备、海况条件等因素影响，观测持续时间都不长，一般在几天到几十天不等。即使是同一任务所设的多个站点，资料的长度、起止时间、数据是否连续等也很难做到统一化。对这类数据作准确的分析和向后预报是工程潮位处理的难点所在。

为提高短期验潮资料的分析精度，前人也研究提出了不同的解决方案。最基本的方案是引入差比关系来解决短期资料调和分析时可能产生的分潮混淆现象，以更准确的获得基本分潮的调和常数。余水位统计法也是近年来较多研究的方案，其基本思想是认为可以通过分离出长期验潮站实测水位中的余水位部分(即非天文潮位部分)，用于改正对应短期验潮站的预报水位。这一思路的两个认识前提是：(1)余水位在一定的区域范围内是稳定的。(2)短期资料能准确预报天文潮位部分，即原有的不准是因为不能准确预测非天文潮位。基于类似的认识，也有学者提出了较新的基于支持向量机(SVM)的模块化潮汐实时预报模型。特别是在对非天文潮部分的预测中，将SVM与灰色模型相结合，并利用粒子群优化算法对SVM的参数进行优化以提高预报精度。以上这些方案的使用，都使短期验潮资料的准确分析成为可能，预报精度得以提高。但是，其中仍存在一些问题：(1)仍依赖于对应的长期验潮站资料或已知的区域差比系数，搜集这些资料现实中有很大难度；(2)部分方案的精度提高指标建立在调和常数精度对比之上，和按逐时实际水位比较仍有区别；(3)部分方案原理及计算过程相对复杂，不能快速用于实际生产。

发明内容

本发明的目的是以较为简捷的方式改善短期验潮数据预报的准确度，从而提供一种短期潮位预报方法。

为实现上述目的，本发明提供的一种短期潮位预报方法，包括：

步骤(一)，对短期验潮资料使用调和分析法进行基本处理，获取主要分潮的调和常数，用于向后预报；记水位方程为：

式1中，ζ表示观测到的实际水位；H和g分别表示分潮的振幅和迟角，即最终要求取的分潮调和常数；

△即水位中未能调和分解的成分，主要包括平均海面水位及非天文潮位；σ是分潮角速度，对不同的分潮而言，该值是一个固定的常数；V₀是分潮天文初相角，由天文变量计算得到；f、u分别为分潮交点因子和交点订正角，实际计算过程中，以fcosu、fsinu的形式给出，同样可以由天文变量计算得到，以上变量中，除了ζ值外，V₀、f、u同样是时间的变量；

步骤(二)，通过扩展卡尔曼滤波(EKF)，将通用的潮汐模型计算水位融合到短序列潮位直接调和预报模型中；

记状态水位为d，观测水位为z，则有EKF下状态方程及观测方程：

d_(t)＝f(d_(t-1))+w_(t) (式2)

z_(t)＝h(d_(t))+v_(t) (式3)

式2、式3中，d(t)是t时刻下的系统状态(真实值)，用t-1时刻下潮位预报方程预报值表示；f(d_(t-1))为已知分潮调和常数后按式(1)计算出的水位值，w(t)为系统过程噪声即预报误差；本文将模型值融入预报值，故记模型在t时刻的计算水位为观测值z(t)，v(t)是模型误差，同时取w(t)、v(t)的方差为q和r；

对式2、式3用泰勒级数展开以线性化，取一阶导数，有：

d_(t)＝f(d_(t-1t-1))+F_(t-1)(d_(t-1)-d_(t-1|t-1))+w_(t) (式4)

z_(t)＝h(d_(t|t-1))+H_(t)(d_(t)-d_(t|t-1))+v_(t) (式5)

其中，F_(t-1)和H_(t)分别表示函数f(d)、h(d)在d_(t-1t-1)、d_(tt-1)处的雅克比矩阵；

基于以上公式，EKF的预测过程为：

d_(tt-1)＝f(d_(t-1t-1))(式6)

估算过程为：

K_(t)＝P_(t|t-1)H_(t) ^T(H_(t)P_(t|t-1)H_(t) ^T+r)^-1 (式8)

d_(t|t)＝d_(t|t-1)+K_(t)(z_(t)-h(d_(t|t-1))) (式9)

P_(t|t)＝(1-K_(t)H_(t))P_(t|t-1) (式10)

式中，P为均方差，K为卡尔曼增益，d_(t|t)即为两种数据融合后所得的最优估计水位值。

优选地，采用NAO.99jb模型值作为融入数据。

本发明利用扩展卡尔曼滤波(EKF)，将通用的潮汐模型计算水位融合到短序列潮位直接调和预报模型之中，不依靠额外高精度数据或计算关系，即可生成比两种源数据精度都更高的融合后预报数据，数据测验表明，对于3天的验潮数据，EKF方法至少对向后5天的预报有效，融合后数据较两种源数据的平均改善度分别为33％和60％；对于7天的验潮数据，EKF方法在向后20的预报中几乎都产生了改善效果，对应两种源数据的平均改善度分别为17％和38％。

本发明以较简捷的方式改善一定的时间范围内短期验潮数据预报的准确度，在不依赖于长期资料或额外的已知关系，尽可能准确的分析短期数据并向后预报，从而解决了工程运用层面的实际需求。

附图说明

图1是实施例中7天数据对比图(2018.01.19～2018.01.25，2016点)；

图2是实施例中20天数据对比图(2018.01.21～2018.02.09，5760点)。

具体实施方式

针对海上临时验潮站验潮时间短、分析难度大、预报精度低的问题。本发明一种典型的实施方式提了一种短期潮位预报方法，利用扩展卡尔曼滤波(EKF)，将通用的潮汐模型计算水位融合到短序列潮位直接调和预报模型之中，以改善短期验潮数据预报的准确度。所述的短期潮位预报方法，包括以下步骤。

步骤(一)，对短期验潮资料使用调和分析法进行基本处理，获取主要分潮的调和常数，用于向后预报；记水位方程为：

式1中，ζ表示观测到的实际水位；H和g分别表示分潮的振幅和迟角，即最终要求取的分潮调和常数；

△即水位中未能调和分解的成分，主要包括平均海面水位及非天文潮位；σ是分潮角速度，对不同的分潮而言，该值是一个固定的常数；V₀是分潮天文初相角，由天文变量计算得到；f、u分别为分潮交点因子和交点订正角，实际计算过程中，以fcosu、fsinu的形式给出，同样可以由天文变量计算得到，以上变量中，除了ζ值外，V₀、f、u同样是时间的变量；

步骤(二)，通过扩展卡尔曼滤波(EKF)，将通用的潮汐模型计算水位融合到短序列潮位直接调和预报模型中；

记状态水位为d，观测水位为z，则有EKF下状态方程及观测方程：

d_(t)＝f(d_(t-1))+w_(t) (式2)

z_(t)＝h(d_(t))+v_(t) (式3)

式2、式3中，d(t)是t时刻下的系统状态(真实值)，用t-1时刻下潮位预报方程预报值表示；f(d_(t-1))为已知分潮调和常数后按式(1)计算出的水位值，w(t)为系统过程噪声即预报误差；本文将模型值融入预报值，故记模型在t时刻的计算水位为观测值z(t)，v(t)是模型误差，同时取w(t)、v(t)的方差为q和r；

对式2、式3用泰勒级数展开以线性化，取一阶导数，有：

d_(t)＝f(d_(t-1|t-1))+F_(t-1)(d_(t-1)-d_(t-1|t-1))+w_(t) (式4)

z_(t)＝h(d_(t|t-1))+H_(t)(d_(t)-d_(t|t-1))+v_(t) (式5)

其中，F_(t-1)和H_(t)分别表示函数f(d)、h(d)在d_(t-1|t-1)、d_(t|t-1)处的雅克比矩阵；

基于以上公式，EKF的预测过程为：

d_(tt-1)＝f(d_(t-1t-1)) (式6)

估算过程为：

K_(t)＝P_(t|t-1)H_(t) ^T(H_(t)P_(t|t-1)H_(t) ^T+r)^-1 (式8)

d_(t|t)＝d_(t|t-1)+K_(t)(z_(t)-h(d_(t|t-1))) (式9)

P_(t|t)＝(1-K_(t)H_(t))P_(t|t-1) (式10)

式中，P为均方差，K为卡尔曼增益，d_(tt)即为两种数据融合后所得的最优估计水位值。

对短期潮位数据进行直接的调和分析，获取最主要几个分潮的调和常数并进行递推计算，是一种不需引入任何附加资料的基本的潮位预报方法。如在江苏沿海半日潮海区，利用几天的数据(如7天之内)分解出最主要的M₂、S₂、O₁、K₁(时长过短对后两个有影响)分潮是可以做到的，其它分潮则较难得到准确、稳定的计算结果。事实上，由于非天文潮位等易变因素的影响，即便使用长时间观测数据，调和分析的最好精度也只在20～30cm之间。

另一种可单独获得预报水位的途径是通过潮汐模型进行理论计算。潮汐模型一般都是建立在水动力学方程之上，通过同化潮汐实测资料或卫星测高资料来提高其计算精度。全球型潮汐模型共通的特点是在外海精度较高，在近海精度普遍偏低。有学者利用验潮站数据对NAO.99b、CSR4.0、TPX07.2等三种当前主流的海潮模型进行了水位的直接对比分析，结果表明，在黄海和东海海域，表现最好的模型数据与实测数据的标准差分别是43.11cm和21.24cm。

以上两种独立方法的分析总体精度近似但有各自的特点。直接调和分析所得的调和常数用于预报，对潮时的定位较准确，但潮高部分受所取的分析资料时间段影响较大。同时，随着预报时长的增加，准确度逐渐降低。潮汐模型方法预报值则没有后一点的困扰，其准确度和预报时间点并无直接关联，但由于近岸地形(底边界)复杂等因素对水动力方程计算的影响，导致在潮时的定位上常与实际有出入。

本发明的整体构思是利用数据融合技术将这两种独立的数据结合起来，优势互补，期望得到较任一单一数据源精度都高的融合后数据，以这种较简捷的方式改善一定时间范围内短期验潮数据预报的准确度。

卡尔曼滤波(KF)是一种利用线性系统状态方程，通过系统输入输出观测数据，对系统状态进行最优估计的算法。但由于潮波状态方程不是线性的，故本发明使用扩展卡尔曼滤波(EKF)方案进行数据融合。扩展卡尔曼滤波是标准卡尔曼滤波在非线性情形下的一种扩展形式，其基本思想是利用泰勒级数展开将非线性系统线性化，然后采用卡尔曼滤波框架对信号进行滤波。

一般而言，分辨率高的区域性潮汐模型计算结果优于全球性潮汐模型。如已掌握了地形边界等各类边界条件，自建区域模型也是好的选择。本发明考虑到通用性，采用NAO.99jb模型值作为融入数据。NAO.99jb是日本国立天文观测台(National AstronomicalObservatory)开发的区域潮汐模型。范围为纬度20°～65°N，经度110°～155°E，分辨率为1/12°×1/12°。事实上，NAO.99jb在动力学模型中同化了卫星测高数据及部分沿岸验潮站数据。

下面结合具体实施例对本发明的技术方案及技术效果作进一步清楚、完整的说明。以下实施例使用扩展卡尔曼滤波(EKF)方案，将潮汐模型数据融入到基于短期验潮数据的调和预报模型中，并将融合前后数据分别与实测数据相对比，分析融合方案的有效程度。

南黄海辐射沙脊海域地形多变、水动力条件复杂，长时间设立验潮站的难度大、潮汐模型的预报精度低。本实施例以位于盐城岸外辐射沙脊海域的短期验潮站为期28天的实测数据作为分析对象，按上述数据融合方案进行实验测试，以期获得‘改善’的合成后预报数据。同时，重点分析这种精度改善随时间的变化规律。考虑到海上临时验潮时长常在一周内，故分别以3天及7天数据进行分析、预报、融合、对比研究。

供分析的验潮站坐标：112.443°E、33.104°N，为半日潮海区。验潮时长2018.01.13～2018.02.09(格林威治时间)，约28天。水位取样间隔5分钟。

(一)三天数据

在实测水位序列中，取3天(2018.01.16～2018.01.18)连续水位数据，调和分析得到M₂、S₂分潮的调和常数并建立预报方程(式1反向计算)。

知道了各个分潮的调和常数(H和g)，等式右半部分的其它参数都是已知的，代入后即可以求出水位ζ。向后预报7天(2018.01.19～2018.01.25)的水位值，作为调和预报值，作对比使用。

采用NAO.99jb潮汐模型值作为融入数据。使用EKF方案在预报方程中融入同时段潮汐模型值(7天)进行数据合成。

将合成后的7天水位值分别与调和预报值、潮汐模型值、实测水位作对比统计分析，共计对比5分钟间隔水位2016点。EKF方案融合时，取过程和测量的标准差为0.3m及0.15m，对应的方差q、r即为0.09及0.0225。

表1：7天数据对比对比统计表(单位：米)

对于3天实测数据的数据融合预报结果，结合图1、表1的分析，认识如下：

(1)在本实施例中,潮汐模型值较调和预报值精度高。这与调和分析采用3天时长水位资料，仅分析M₂、S₂两个分潮而导致本身预报精度偏低有关。

(2)使用EKF方法融合后数据，在预报的前5天较两种原始数据，均值及STD值普遍得到改善。均值较模型平均改善6.3cm改善度33％(本实施例改善度计算为(0.19-0.127)/0.19×100％，下类同)、较调和预报平均改善18.7cm改善度60％；STD值较模型平均改善2.9cm改善度25％、较调和预报平均改善13.1cm改善度60％。由此可见，使用EKF方法的融合后数据，至少在本实施例的前5天内，比使用任何单一方案(潮汐模型或调和预报)精度都明显提高。并且，这种提高仅是数据融合的结果，而不是另外引入了更高精度的数据或其它计算关系。

(3)从第6天开始，融合后数据在均值指标上较模型数据精度低，主要原因是作为源数据的调和预报数据精度过低引起的(使用三天大潮期数据直接计算2个分潮的调和常数用于小潮期的预报时精度较低)。

(4)由上一条可知，源数据本身的精度对数据融合效果有明显的影响。结合一系列测试认为，当两种源数据本身的精度近似时，融合后数据精度会更理想。

(5)另外，使用EKF方案进行数据融合时，过程及测量标准差的取值对结果影响很大，需要在广泛测试的基础上选取合适的值。

(二)七天数据

类似方法，在实测水位序列中，取7天(2018.01.14～2018.01.20)连续水位数据，调和分析得到M₂、S₂、O₁、K₁分潮的调和常数并建立预报方程，同时向后预报20天(2018.01.21～2018.02.09)的水位值(调和预报值)作对比使用。

采用NAO.99jb潮汐模型值作为融入数据。使用EKF方案在预报方程中融入同时段潮汐模型值(20天)进行数据合成，将合成后的20天水位值分别与调和预报值、潮汐模型值、实测水位作对比统计分析，共计对比5分钟间隔水位5760点。EKF方案融合时，取过程和测量的标准差为0.3m及0.15m，对应的方差q、r即为0.09及0.0225。

表2：20天数据对比对比统计表

由图2、表2可知，随着用于调和分析的水位资料延长为7天，解算分潮扩展为4个，调和预报本身的精度得到提高(小潮期仍影响较大)。故而，和采用3天资料不同，在持续20天的预报中，融合后数据精度较两种源数据普遍是提高的。其中，均值较模型平均改善6.3cm改善度17％、较调和预报平均改善18.5cm改善度38％；STD值较模型平均改善4.4cm改善度19％、较调和预报平均改善12.5cm改善度40％。

总的来说，EKF方法能有效合成出比任何一种源数据精度更高的融合数据。且随着源数据本身精度的提高，这种‘有效合成’的持续时间更长。

小结

本发明使用扩展卡尔曼数据融合方法，将NAO.99jb潮汐模型数据和短序列水位调和预报数据融合在一起，形成新的融合后预报数据。通过在不同环境下对该数据的分析对比试验，形成主要认识如下：

(1)无论使用三天或七天序列的水位数据进行分析、融合，融合后数据至少在一定的持续时间内，能比任何一种源数据有明显的精度提高，证明EKF数据融合方案是行之有效的。

(2)源数据本身的精度情况对融合后数据精度有重要影响，直接的表现就是有效融合的持续时间不同。理想的状况是两种源数据精度高且接近。

(3)EKF数据融合过程中，过程及测量标准差等参数的取值同样对融合后数据精度影响很大，实际使用时需对各种参数组合进行数据试验。

总之，本发明提供的数据处理方案，使用的源资料易于获得，不需要额外长期资料数据或某些已知的关系，且计算方法简洁便于程序实现，试验效果理想。可以在资料实际较少的情况下测试使用，以提高数据精度，降低工作成本。

Claims (2)

1.一种短期潮位预报方法，其特征在于，包括：

步骤(一)，对短期验潮资料使用调和分析法进行基本处理，获取主要分潮的调和常数，用于向后预报；记水位方程为：

式1中，ζ表示观测到的实际水位；H和g分别表示分潮的振幅和迟角，即最终要求取的分潮调和常数；

△即水位中未能调和分解的成分，主要包括平均海面水位及非天文潮位；σ是分潮角速度，对不同的分潮而言，该值是一个固定的常数；V₀是分潮天文初相角，由天文变量计算得到；f、u分别为分潮交点因子和交点订正角，实际计算过程中，以fcosu、fsinu的形式给出，同样可以由天文变量计算得到，以上变量中，除了ζ值外，V₀、f、u同样是时间的变量；

步骤(二)，通过扩展卡尔曼滤波(EKF)，将通用的潮汐模型计算水位融合到短序列潮位直接调和预报模型中；

记状态水位为d，观测水位为z，则有EKF下状态方程及观测方程：

d_(t)＝f(d_(t-1))+w_(t) (式2)

z_(t)＝h(d_(t))+v_(t) (式3)

式2、式3中，d(t)是t时刻下的系统状态(真实值)，用t-1时刻下潮位预报方程预报值表示；f(d_(t-1))为已知分潮调和常数后按式(1)计算出的水位值，w(t)为系统过程噪声即预报误差；本文将模型值融入预报值，故记模型在t时刻的计算水位为观测值z(t)，v(t)是模型误差，同时取w(t)、v(t)的方差为q和r；

对式2、式3用泰勒级数展开以线性化，取一阶导数，有：

d_(t)＝f(d_(t-1|t-1))+F_(t-1)(d_(t-1)-d_(t-1|t-1))+w_(t) (式4)

z_(t)＝h(d_(t|t-1))+H_(t)(d_(t)-d_(t|t-1))+v_(t) (式5)

其中，F_(t-1)和H_(t)分别表示函数f(d)、h(d)在d_(t-1|t-1)、d_(t|t-1)处的雅克比矩阵；

基于以上公式，EKF的预测过程为：

d_(t|t-1)＝f(d_(t-1|t-1)) (式6)

估算过程为：

K_(t)＝P_(t|t-1)H_(t) ^T(H_(t)P_(t|t-1)H_(t) ^T+r)^-1 (式8)

d_(t|t)＝d_(t|t-1)+K_(t)(z_(t)-h(d_(t|t-1))) (式9)

P_(t|t)＝(1-K_(t)H_(t))P_(t|t-1) (式10)

式中，P为均方差，K为卡尔曼增益，d_(t|t)即为两种数据融合后所得的最优估计水位值。

2.根据权利要求1所述的方法，其特征在于：采用NAO.99jb模型值作为融入数据。

Images