[0001] 本
申请是申请号为：2014108415648，
发明创造名称为《一种船舶轨迹实时预测方法》，申请日为：2014年12月30日的发明
专利申请的分案申请。
技术领域
[0002] 本发明涉及一种海域交通管制方法，尤其涉及一种基于滚动规划策略的船舶轨迹预测方法。
背景技术
[0003] 随着全球航运业的快速发展，部分繁忙海域内的交通愈加拥挤。在船舶交通流密集复杂海域，针对船舶间的冲突情形仍然采用航行计划结合人工间隔调配的管制方式已不能适应航运业的快速发展。为保证船舶间的安全间隔，实施有效的冲突调配就成为海域交通管制工作的重点。船舶冲突解脱是航海领域中的一项关键技术，安全高效的解脱方案对于增加海域船舶流量以及确保海运安全具有重大意义。
[0004] 为了提高船舶的航行效率，船用雷达自动标绘仪目前已经被广泛应用到船舶监控和避碰中，该设备通过提取船舶相关信息为船舶间冲突情形的判定提供参考依据。尽管此类设备极大降低了人工监控的负荷，但它并不具备船舶自动冲突解脱功能。而船舶冲突解脱是基于对船舶轨迹的预测的
基础上，在船舶实际航行中，受气象条件、导航设备以及驾驶员操作等各种因素的影响，它的运行状态往往不完全属于某一特定的运动状态，在船舶轨迹预测过程中需要考虑各种随机因素的影响，通过获取各类随机因素的最新特性对其未来轨迹实施滚动预测并增强其轨迹预测的鲁棒性。
发明内容
[0005] 本发明要解决的技术问题是提供一种鲁棒性较好的船舶轨迹预测方法，该方法的船舶轨迹预测
精度较高。
[0006] 实现本发明目的的技术方案是提供一种船舶轨迹预测方法，包括如下几个步骤：
[0007] ①通过海面雷达获得船舶的实时和历史
位置信息，各船舶的位置信息为离散二维位置序列x'＝[x1',x2',...,xn']和y'＝[y1',y2',...,yn']，通过应用
小波变换理论对原始离散二维位置序列x'＝[x1',x2',...,xn']和y'＝[y1',y2',...,yn']进行初步处理，从而获取船舶的去噪离散二维位置序列x＝[x1,x2,...,xn]和y＝[y1,y2,...,yn]；
[0008] ②在每一
采样时刻对船舶轨迹数据预处理，依据所获取的船舶原始离散二维位置序列x＝[x1,x2,...,xn]和y＝[y1,y2,...,yn]，采用一阶差分方法对其进行处理获取新的船舶离散位置序列Δx＝[Δx1,Δx2,...,Δxn-1]和Δy＝[Δy1,Δy2,...,Δyn-1]，其中Δxi＝xi+1-xi,Δyi＝yi+1-yi(i＝1,2,...,n-1)；
[0009] ③在每一采样时刻对船舶轨迹数据聚类，对处理后新的船舶离散二维位置序列Δx和Δy，通过设定聚类个数M'，采用K-means聚类
算法分别对其进行聚类；
[0010] ④在每一采样时刻对船舶轨迹数据利用隐
马尔科夫模型进行参数训练，通过将处理后的船舶运行轨迹数据Δx和Δy视为隐马尔科夫过程的显观测值，通过设定隐状态数目N和参数更新时段τ'，依据最近的T'个位置观测值并采用B-W算法滚动获取最新隐马尔科夫模型参数λ'；
[0011] ⑤在每一采样时刻依据隐马尔科夫模型参数，采用Viterbi算法获取当前时刻观测值所对应的隐状态q；
[0012] ⑥在每一采样时刻，通过设定预测时域W，基于船舶当前时刻的隐状态q，获取未来时段船舶的位置预测值O，从而在每一采样时刻滚动推测到未来时段内船舶的轨迹。
[0013] 进一步的，所述步骤①中，通过应用小波变换理论对原始离散二维位置序列x'＝[x1',x2',...,xn']和y'＝[y1',y2',...,yn']进行初步处理，从而获取船舶的去噪离散二维位置序列x＝[x1,x2,...,xn]和y＝[y1,y2,...,yn]：对于给定的原始二维序列数据x'＝[x1',x2',...,xn']，利用如下形式的线性表达式分别对其进行近似：
[0014] 其中：
[0015] f'(x')表示对数据平滑处理后得到的函数表达式，ψ(x')表示母波，δ、J和K均为小波变换常数，ψJ,K(x')表示母波的转换形式，cJ,K表示由小波变换过程得到的函数系数，它体现了子波ψJ,K(x')对整个函数近似的权重大小，若此系数很小，那么它意味着子波ψJ,K(x')的权重也较小，因而可以在不影响函数主要特性的前提下，从函数近似过程中将子波ψJ,K(x')除去；在实际
数据处理过程中，通过设定
阈值χ来实施“阈值转换”，当cJ,K＜χ时，设定cJ,K＝0；阈值函数的选取采用如下两种方式：
[0016] 和
[0017] 对于y'＝[y1',y2',...,yn']，也采用上述方法进行去噪处理。
[0018] 进一步的，所述步骤④中确定航迹隐马尔科夫模型参数λ'＝(π,A,B)的过程如下：
[0019] 4.1)变量赋初值：应用均匀分布给变量πi，aij和bj(ok)赋初值πi0， 和 并使其满足约束条件： 和 由此得到λ0＝(π0,A0,B0)，其中ok表示某一显观测值，π0、A0和B0分别是由元素 和 构成的矩阵，令参数l＝0，o＝(ot-T'+1,...,ot-1,ot)为当前时刻t之前的T'个历史位置观测值；
[0020] 4.2)执行E-M算法：
[0021] 4.2.1)E-步骤：由λl计算ξe(i,j)和γe(si)；
[0022] 变量 那么
[0023] 其中s表示某一隐状态；
[0024] 4.2.2)M-步骤：运用 分别估计πi，aij和bj(ok)并由此得到λl+1；
[0025] 4.2.3)循环：l＝l+1，重复执行E-步骤和M-步骤，直至πi、aij和bj(ok)收敛，即[0026] |P(o|λl+1)-P(o|λl)|＜ε，其中参数ε＝0.00001，返回步骤4.2.4)；
[0027] 4.2.4)：令λ'＝λl+1，算法结束。
[0028] 进一步的，所述步骤⑤中确定船舶航迹最佳隐状态序列的
迭代过程如下：
[0029] 5.1)变量赋初值：令g＝2，βT'(si)＝1(si∈S)，δ1(si)＝πibi(o1)，ψ1(si)＝0，其中，[0030]
[0031] ，其中变量ψg(sj)表示使变量δg-1(si)aij取最大值的船舶航迹隐状态si，参数S表示隐状态的集合；
[0032] 5.2)递推过程：
[0033] 5.3)时刻更新：令g＝g+1，若g≤T'，返回步骤5.2)，否则迭代终止并转到步骤5.4)；
[0034] 5.4) 转到步骤5.5)；
[0035] 5.5)最优隐状态序列获取：
[0036] 5.5.1)变量赋初值：令g＝T'-1；
[0037] 5.5.2)后向递推：
[0038] 5.5.3)时刻更新：令g＝g-1，若g≥1，返回步骤5.5.2)，否则终止。
[0039] 进一步的，所述步骤③中，聚类个数M'的值为4。
[0040] 进一步的，所述步骤④中，状态数目N的值为3，参数更新时段τ'为30秒，T'为10。
[0041] 进一步的，所述步骤⑥中，预测时域W为300秒。
[0042] 本发明具有积极的效果：(1)本发明在船舶轨迹实时预测的过程中，融入了随机因素的影响，所采用的滚动轨迹预测方案能够及时提取外界随机因素的变化状况，提高了船舶轨迹预测的准确性。
[0043] (2)本发明基于不同性能指标，其船舶轨迹实时预测结果可以为存在冲突的多个船舶提供解
脱轨迹规划方案，提高船舶运行的经济性和海域资源的利用率。
附图说明
[0044] 图1为本发明中的船舶运行短期轨迹生成流程示意图。
具体实施方式
[0046] 见图1，本实施例的一种船舶轨迹预测方法包括如下几个步骤：
[0047] ①通过海面雷达获得船舶的实时和历史位置信息，各船舶的位置信息为离散二维位置序列x'＝[x1',x2',...,xn']和y'＝[y1',y2',...,yn']，通过应用小波变换理论对原始离散二维位置序列x'＝[x1',x2',...,xn']和y'＝[y1',y2',...,yn']进行初步处理，从而获取船舶的去噪离散二维位置序列x＝[x1,x2,...,xn]和y＝[y1,y2,...,yn]：y＝[y1,y2,...,yn]：对于给定的原始二维序列数据x'＝[x1',x2',...,xn']，利用如下形式的线性表达式分别对其进行近似：
[0048] 其中：
[0049] f'(x')表示对数据平滑处理后得到的函数表达式，ψ(x')表示母波，δ、J和K均为小波变换常数，ψJ,K(x')表示母波的转换形式，cJ,K表示由小波变换过程得到的函数系数，它体现了子波ψJ,K(x')对整个函数近似的权重大小，若此系数很小，那么它意味着子波ψJ,K(x')的权重也较小，因而可以在不影响函数主要特性的前提下，从函数近似过程中将子波ψJ,K(x')除去；在实际数据处理过程中，通过设定阈值χ来实施“阈值转换”，当cJ,K＜χ时，设定cJ,K＝0；阈值函数的选取采用如下两种方式：
[0050] 和
[0051] 对于y'＝[y1',y2',...,yn']，也采用上述方法进行去噪处理；
[0052] ②在每一采样时刻对船舶轨迹数据预处理，依据所获取的船舶原始离散二维位置序列x＝[x1,x2,...,xn]和y＝[y1,y2,...,yn]，采用一阶差分方法对其进行处理获取新的船舶离散位置序列Δx＝[Δx1,Δx2,...,Δxn-1]和Δy＝[Δy1,Δy2,...,Δyn-1]，其中Δxi＝xi+1-xi,Δyi＝yi+1-yi(i＝1,2,...,n-1)；
[0053] ③在每一采样时刻对船舶轨迹数据聚类，对处理后新的船舶离散二维位置序列Δx和Δy，通过设定聚类个数M'，采用K-means
聚类算法分别对其进行聚类；
[0054] ④在每一采样时刻对船舶轨迹数据利用隐马尔科夫模型进行参数训练，通过将处理后的船舶运行轨迹数据Δx和Δy视为隐马尔科夫过程的显观测值，通过设定隐状态数目N和参数更新时段τ'，依据最近的T'个位置观测值并采用B-W算法滚动获取最新隐马尔科夫模型参数λ'；确定航迹隐马尔科夫模型参数λ'＝(π,A,B)的过程如下：
[0055] 4.1)变量赋初值：应用均匀分布给变量πi，aij和bj(ok)赋初值 和 并使其满足约束条件： 和 由此得
到λ0＝(π0,A0,B0)，其中ok表示某一显观测值，π0、A0和B0分别是由元素 和 构成的矩阵，令参数l＝0，o＝(ot-T'+1,...,ot-1,ot)为当前时刻t之前的T'个历史位置观测值；
[0056] 4.2)执行E-M算法：
[0057] 4.2.1)E-步骤：由λl计算ξe(i,j)和γe(si)；
[0058] 变量 那么
[0059] 其中s表示某一隐状态；
[0060] 4.2.2)M-步骤：运用 分别估计πi，aij和bj(ok)并由此得到λl+1；
[0061] 4.2.3)循环：l＝l+1，重复执行E-步骤和M-步骤，直至πi、aij和bj(ok)收敛，即[0062] |P(o|λl+1)-P(o|λl)|＜ε，其中参数ε＝0.00001，返回步骤4.2.4)；
[0063] 4.2.4)：令λ'＝λl+1，算法结束。
[0064] ⑤在每一采样时刻依据隐马尔科夫模型参数，采用Viterbi算法获取当前时刻观测值所对应的隐状态q：
[0065] 5.1)变量赋初值：令g＝2，βT'(si)＝1(si∈S)，δ1(si)＝πibi(o1)，ψ1(si)＝0，其中，[0066]
[0067] ，其中变量ψg(sj)表示使变量δg-1(si)aij取最大值的船舶航迹隐状态si，参数S表示隐状态的集合；
[0068] 5.2)递推过程：
[0069] 5.3)时刻更新：令g＝g+1，若g≤T'，返回步骤5.2)，否则迭代终止并转到步骤5.4)；
[0070] 5.4) 转到步骤5.5)；
[0071] 5.5)最优隐状态序列获取：
[0072] 5.5.1)变量赋初值：令g＝T'-1；
[0073] 5.5.2)后向递推：
[0074] 5.5.3)时刻更新：令g＝g-1，若g≥1，返回步骤5.5.2)，否则终止。。
[0075] ⑥在每一采样时刻，通过设定预测时域W，基于船舶当前时刻的隐状态q，获取未来时段船舶的位置预测值O。
[0076] 上述聚类个数M'的值为4，，状态数目N的值为3，参数更新时段τ'为30秒，T'为10，预测时域W为300秒。
[0078] 本实施例的航海交通管制方法包括如下几个步骤：
[0079] 步骤A、根据实施例1得到的船舶轨迹预测方法获得船舶在每一采样时刻推测到的未来时段内船舶的轨迹；
[0080] 步骤B、在每一采样时刻，基于船舶当前的运行状态和历史位置观察序列，获取海域
风场变量的数值，其具体过程如下：
[0081] B.1)设定船舶的停靠位置为轨迹参考坐标原点；
[0082] B.2)在船舶处于直线运行状态和匀速转弯运行状态时，构建海域风场线性滤波模型；
[0083] B.3)根据所构建的滤波模型获取风场变量的数值。
[0084] 步骤C、在每一采样时刻，基于各船舶的运行状态和设定的船舶在海域内运行时需满足的安全规则集，当船舶间有可能出现违反安全规则的状况时，对其动态行为实施监控并为海上交通控制中心提供及时的告警信息；
[0085] 步骤D、当告警信息出现时，在满足船舶物理性能和海域交通规则的前提下，通过设定优化指标函数以及融入风场变量数值，采用模型预测控制理论方法对船舶避撞轨迹进行滚动规划，并将规划结果传输给各船舶执行，其具体过程如下：
[0086] D.1)设定船舶避撞轨迹规划的终止参考点位置P、避撞策略控制时域Θ、轨迹预测时域γ；
[0087] D.2)设定在给定优化指标函数的前提下，基于合作式避撞轨迹规划思想，通过给各个船舶赋予不同的权重以及融入实时风场变量滤
波数值，得到各个船舶的避撞轨迹和避撞控制策略并将规划结果传输给各船舶执行，且各船舶在滚动规划间隔内仅实施其第一个优化控制策略；
[0088] D.3)在下一采样时刻，重复步骤5.2)直至各船舶均到达其解脱终点。
[0089] 上述终止参考点位置P设定为船舶位置冲突点的下一个航道点，避撞策略控制时域Θ为300秒；轨迹预测时域γ为300秒。
[0090] 显然，上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。而这些属于本发明的精神所引伸出的显而易见的变化或变动仍处于本发明的保护范围之中。
