CN104462856A - 船舶冲突预警方法 - Google Patents

Abstract

本发明涉及一种船舶冲突预警方法，包括如下几个步骤，首先通过海面雷达获得船舶的实时和历史位置信息并做初步处理；然后在每一采样时刻对船舶轨迹数据预处理，然后在每一采样时刻对船舶轨迹数据聚类，再而在每一采样时刻对船舶轨迹数据利用隐马尔科夫模型进行参数训练，然后在每一采样时刻依据隐马尔科夫模型参数，采用Viterbi算法获取当前时刻观测值所对应的隐状态q，最后在每一采样时刻通过设定预测时域W，基于船舶当前时刻的隐状态q，获取未来时段船舶的位置预测值O，并通过建立安全规则集对船舶的动态行为实施监控并及时发出告警信息。本发明滚动实时对船舶轨迹进行预测，有效预警海域冲突，提高海上交通的安全性。

Description

船舶冲突预警方法

技术领域

本发明涉及一种海域交通管制方法，尤其涉及一种基于滚动规划策略的船舶冲突预警方法。

背景技术

随着全球航运业的快速发展，部分繁忙海域内的交通愈加拥挤。在船舶交通流密集复杂海域，针对船舶间的冲突情形仍然采用航行计划结合人工间隔调配的管制方式已不能适应航运业的快速发展。为保证船舶间的安全间隔，实施有效的冲突预警就成为海域交通管制工作的重点。船舶冲突预警是航海领域中的一项关键技术，安全高效的预警方案对于增加海域船舶流量以及确保海运安全具有重大意义。

为了提高船舶的航行效率，船用雷达自动标绘仪目前已经被广泛应用到船舶监控和避碰中，该设备通过提取船舶相关信息为船舶间冲突情形的判定提供参考依据。而船舶冲突预警是基于对船舶轨迹的预测的基础上，在船舶实际航行中，受气象条件、导航设备以及驾驶员操作等各种因素的影响，它的运行状态往往不完全属于某一特定的运动状态，故目前对船舶轨迹的预测及船舶冲突预警均没有较为准确的方案。

发明内容

本发明要解决的技术问题是提供一种鲁棒性较好的船舶冲突预警方法，该方法的船舶轨迹预测精度较高，船舶冲突预警的准确性及时效性均较好。

实现本发明目的的技术方案是提供一种船舶冲突预警方法，包括如下几个步骤：

①通过海面雷达获得船舶的实时和历史位置信息，各船舶的位置信息为离散二维位置序列x′＝[x₁′，x₂′，...，x_n′]和y′＝[y₁′，y₂′，...，y_n′]，通过应用小波变换理论对原始离散二维位置序列x′＝[x₁′，x₂′，...，x_n′]和y′＝[y₁′，y₂′，...，y_n′]进行初步处理，从而获取船舶的去噪离散二维位置序列x＝[x₁，x₂，...，x_n]和y＝[y₁，y₂，...，y_n]；

②在每一采样时刻对船舶轨迹数据预处理，依据所获取的船舶原始离散二维位置序列x＝[x₁，x₂，...，x_n]和y＝[y₁，y₂，...，y_n]，采用一阶差分方法对其进行处理获取新的船舶离散位置序列Δx＝[Δx₁，Δx₂，...，Δx_n-1]和Δy＝[Δy₁，Δy₂，...，Δy_n-1]，其中Δx_i＝x_i+1-x_i，Δ_yi＝y_i+1-y_i(i＝1，2，...，n-1)；

③在每一采样时刻对船舶轨迹数据聚类，对处理后新的船舶离散二维位置序列Δx和Δy，通过设定聚类个数M′，采用遗传聚类算法分别对其进行聚类；

④在每一采样时刻对船舶轨迹数据利用隐马尔科夫模型进行参数训练，通过将处理后的船舶运行轨迹数据Δx和Δy视为隐马尔科夫过程的显观测值，通过设定隐状态数目N和参数更新时段τ′，依据最近的T′个位置观测值并采用B-W算法滚动获取最新隐马尔科夫模型参数λ′；

⑤在每一采样时刻依据隐马尔科夫模型参数，采用Viterbi算法获取当前时刻观测值所对应的隐状态q；

⑥在每一采样时刻，通过设定预测时域W，基于船舶当前时刻的隐状态q，获取未来时段船舶的位置预测值O，从而在每一采样时刻滚动推测到未来时段内船舶的轨迹；

⑦在每一采样时刻，基于各船舶的运行状态和设定的船舶在海域内运行时需满足的安全规则集，当船舶间有可能出现违反安全规则的状况时，对其动态行为实施监控并为海上交通控制中心提供及时的告警信息。

进一步的，所述步骤①中，通过应用小波变换理论对原始离散二维位置序列x′＝[x₁′，x₂′，...，x_n′]和y′＝[y₁′，y₂′，...，y_n′]进行初步处理，从而获取船舶的去噪离散二维位置序列x＝[x₁，x₂，...，x_n]和y＝[y₁，y₂，...，y_n]：对于给定的原始二维序列数据x′＝[x₁′，x₂′，...，x_n′]，利用如下形式的线性表达式分别对其进行近似： $\tilde{f}\left(x^{\′}\right)=\underset{\underset{\∀K}{\∀J}}{\mathrm{\Σ}}{c}_{J,K}{\mathrm{\ψ}}_{J,K}\left(x^{\′}\right),$

其中： $\begin{array}{c}{c}_{J,K}={\∫}_{-\∞}^{\∞}{f}^{\′}\left(x^{\′}\right){\mathrm{\ψ}}_{J,K}\left(x^{\′}\right){\mathrm{dx}}^{\′}\\ {\mathrm{\ψ}}_{J,K}\left(x^{\′}\right)=\mathrm{\δ}\·\mathrm{\ψ}(2^J{x}^{\′}-K)\end{array},$

f′(x′)表示对数据平滑处理后得到的函数表达式，ψ(x′)表示母波，δ、J和K均为小波变换常数，ψ_J，K(x′)表示母波的转换形式，c_J，K表示由小波变换过程得到的函数系数，它体现了子波ψ_J，K(x′)对整个函数近似的权重大小，若此系数很小，那么它意味着子波ψ_J，K(x′)的权重也较小，因而可以在不影响函数主要特性的前提下，从函数近似过程中将子波ψ_J，K(x′)除去；在实际数据处理过程中，通过设定阈值χ来实施“阈值转换”，当c_J，K＜χ时，设定c_J，K＝0；阈值函数的选取采用如下两种方式：

${\mathrm{\ρ}}_1(d,\mathrm{\χ})=\left\{\begin{array}{cc}d& \mathrm{if}\left|d\right|>\mathrm{\χ}\\ 0& \mathrm{if}\left|d\right|\≤\mathrm{\χ}\end{array}\right.$ 和 ${\mathrm{\ρ}}_2(d,\mathrm{\χ})=\left\{\begin{array}{cc}d-\left(\frac{d}{\left|d\right|}\mathrm{\χ}\right)& \mathrm{if}\left|d\right|>\mathrm{\χ}\\ 0& \mathrm{if}\left|d\right|\≤\mathrm{\χ}\end{array}\right.;$

对于y′＝[y₁′，y₂′，...，y_n′]，也采用上述方法进行去噪处理。

进一步的，所述步骤④中确定航迹隐马尔科夫模型参数λ′＝(π，A，B)的过程如下：

4.1)变量赋初值：应用均匀分布给变量π_i，a_ij和b_i(o_k)赋初值和并使其满足约束条件： ${\mathrm{\Σ}}_{i=1}^N{\mathrm{\π}}_i^0=1,{\mathrm{\Σ}}_{j=1}^N{a}_{\mathrm{ij}}^0=1(1\≤i\≤N)$ 和 ${\mathrm{\Σ}}_{k=1}^{M^{\′}}{b}_j^0\left(o_k\right)=1(1\≤j\≤N),$ 由此得到λ₀＝(π₀，A₀，B₀)，其中o_k表示某一显观测值，π₀、A₀和B₀分别是由元素和构成的矩阵，令参数l＝0，o＝(o_t-T′+1，...，o_t-1，o_t)为当前时刻t之前的T′个历史位置观测值；

4.2)执行E-M算法：

4.2.1)E-步骤：由λ_l计算ξ_e(i，j)和γ_e(s_i)；

变量 ${\mathrm{\ξ}}_e(i,j)=\frac{P(q_e=s_i,q_{e+1}=s_j,o|{\mathrm{\λ}}_l)}{P\left(o\right|{\mathrm{\λ}}_l)},$ 那么 ${\mathrm{\γ}}_e\left(s_i\right)={\mathrm{\Σ}}_{j=1}^N{\mathrm{\ξ}}_e(i,j),$

其中s表示某一隐状态；

4.2.2)M-步骤：运用 ${\stackrel{\‾}{\mathrm{\π}}}_i={\mathrm{\γ}}_1\left(s_i\right),{\stackrel{\‾}{a}}_{\mathrm{ij}}=\frac{{\mathrm{\Σ}}_{e=1}^{T^{\′}-1}{\mathrm{\ξ}}_e(i,j)}{{\mathrm{\Σ}}_{e=1}^{T^{\′}-1}{\mathrm{\γ}}_e\left(s_i\right)},{\stackrel{\‾}{b}}_j\left(o_k\right)=\frac{\underset{o_e=o_k}{{\mathrm{\Σ}}_{e=1}^{T^{\′}}{\mathrm{\γ}}_e\left(s_j\right)}}{{\mathrm{\Σ}}_{e=1}^{T^{\′}}{\mathrm{\γ}}_e\left(s_j\right)}$ 分别估计π_i，a_ij和b_j(o_k)并由此得到λ_l+1；

4.2.3)循环：l＝l+1，重复执行E-步骤和M-步骤，直至π_i、a_ij和b_j(o_k)收敛，即

|P(o|λ_l+1)-P(o|λ_l)|＜ε，其中参数ε＝0.00001，返回步骤4.2.4)；

4.2.4)：令λ′＝λ_l+1，算法结束。

进一步的，所述步骤⑤中确定船舶航迹最佳隐状态序列的迭代过程如下：

5.1)变量赋初值：令g＝2，β_T′(s_i)＝1(s_i∈S)，δ₁(s_i)＝π_ib_i(o₁)，ψ₁(s_i)＝0，其中，

${\mathrm{\δ}}_g\left(s_i\right)=\underset{q_1,q_2,...,q_{g-1}}{\max }P(q_1,q_2,...,q_{g-1},q_g=s_i,o_1,o_2,...,o_g|{\mathrm{\λ}}^{\′})$

，其中变量ψ_g(s_j)表示使变量δ_g-1(s_i)a_ij取最大值的船舶航迹隐状态s_i，参数S表示隐状态的集合；

5.2)递推过程： ${\mathrm{\δ}}_g\left(s_j\right)=\underset{s_i\∈S}{\max }\left[{\mathrm{\δ}}_{g-1}\left(s_i\right)a_{\mathrm{ij}}\right]b_j\left(o_g\right),{\mathrm{\ψ}}_g\left(s_j\right)=\underset{s_i\∈S}{\arg \max }\left[{\mathrm{\δ}}_{g-1}\left(s_i\right)a_{\mathrm{ij}}\right];$

5.3)时刻更新：令g＝g+1，若g≤T′，返回步骤5.2)，否则迭代终止并转到步骤5.4)；

5.4) $p^{*}=\underset{s_i\∈S}{\max }\left[{\mathrm{\δ}}_{T^{\′}}\left(s_i\right)\right],q_{T^{\′}}^{*}=\underset{s_i\∈S}{\arg \max }\left[{\mathrm{\δ}}_{T^{\′}}\left(s_i\right)\right],$ 转到步骤5.5)；

5.5)最优隐状态序列获取：

5.5.1)变量赋初值：令g＝T′-1；

5.5.2)后向递推：

5.5.3)时刻更新：令g＝g-1，若g≥1，返回步骤5.5.2)，否则终止。

进一步的，所述步骤③中，聚类个数M′的值为4。

进一步的，所述步骤④中，状态数目N的值为3，参数更新时段τ′为30秒，T′为10。

进一步的，所述步骤⑥中，预测时域W为300秒。

进一步的，所述步骤⑦中对各船舶的动态行为实施监控并为海上交通控制中心提供及时的告警信息的具体过程如下：

7.1)构造船舶在海域内运行时需满足的安全规则集D_mr(t)≥D_min，其中D_mr(t)表示任意两个船舶m和船舶r在t时刻的距离，D_min表示船舶间的最小安全距离；

7.2)依据采样时间，建立由船舶连续运行状态至离散采样状态的观测器Λ：г→Ξ，其中г表示船舶的连续运行状态，Ξ表示船舶的离散采样状态；

7.3)当船舶m和r的观测器Λ_m和Λ_r的离散观测数值Ξ_m和Ξ_r在t时刻表明该向量不在安全规则集中时，即关系式D_mr(t)≥D_min不成立时，立刻向海上交通控制中心发出告警信息。

本发明具有积极的效果：(1)本发明在船舶轨迹实时预测的过程中，融入了随机因素的影响，所采用的滚动轨迹预测方案能够及时提取外界随机因素的变化状况，提高了船舶轨迹预测的准确性。

(2)本发明基于不同性能指标，其船舶轨迹实时预测结果可以为存在冲突的多个船舶提供解脱轨迹规划方案，提高船舶运行的经济性和海域资源的利用率。

(3)本发明对船舶冲突的预警效果较好，可有效、准确、实时地预测船舶的轨迹并预测船舶冲突，有效提高海域交通的安全性。

附图说明

图1为本发明中的船舶运行短期轨迹生成流程示意图；

图2为本发明中的船舶运行态势监控流程示意图。

具体实施方式

(实施例1)

见图1，本实施例的船舶冲突预警方法包括如下几个步骤：

①通过海面雷达获得船舶的实时和历史位置信息，各船舶的位置信息为离散二维位置序列x′＝[x₁′，x₂′，...，x_n′]和y′＝[y₁′，y₂′，...，y_n′]，通过应用小波变换理论对原始离散二维位置序列x′＝[x₁′，x₂′，...，x_n′]和y′＝[y₁′，y₂′，...，y_n′]进行初步处理，从而获取船舶的去噪离散二维位置序列x＝[x₁，x₂，...，x_n]和y＝[y₁，y₂，...，y_n]，对于给定的原始二维序列数据x′＝[x₁′，x₂′，...，x_n′]，利用如下形式的线性表达式分别对其进行近似： $\tilde{f}\left(x^{\′}\right)=\underset{\underset{\∀K}{\∀J}}{\mathrm{\Σ}}{c}_{J,K}{\mathrm{\ψ}}_{J,K}\left(x^{\′}\right),$

其中： $\begin{array}{c}{c}_{J,K}={\∫}_{-\∞}^{\∞}{f}^{\′}\left(x^{\′}\right){\mathrm{\ψ}}_{J,K}\left(x^{\′}\right){\mathrm{dx}}^{\′}\\ {\mathrm{\ψ}}_{J,K}\left(x^{\′}\right)=\mathrm{\δ}\·\mathrm{\ψ}(2^J{x}^{\′}-K)\end{array},$

f′(x′)表示对数据平滑处理后得到的函数表达式，ψ(x′)表示母波，δ、J和K均为小波变换常数，ψ_J，K(x′)表示母波的转换形式，c_J，K表示由小波变换过程得到的函数系数，它体现了子波ψ_J，K(x′)对整个函数近似的权重大小，若此系数很小，那么它意味着子波ψ_J，K(x′)的权重也较小，因而可以在不影响函数主要特性的前提下，从函数近似过程中将子波ψ_J，K(x′)除去；在实际数据处理过程中，通过设定阈值χ来实施“阈值转换”，当c_J，K＜χ时，设定c_J，K＝0；阈值函数的选取采用如下两种方式：

${\mathrm{\ρ}}_1(d,\mathrm{\χ})=\left\{\begin{array}{cc}d& \mathrm{if}\left|d\right|>\mathrm{\χ}\\ 0& \mathrm{if}\left|d\right|\≤\mathrm{\χ}\end{array}\right.$ 和 ${\mathrm{\ρ}}_2(d,\mathrm{\χ})=\left\{\begin{array}{cc}d-\left(\frac{d}{\left|d\right|}\mathrm{\χ}\right)& \mathrm{if}\left|d\right|>\mathrm{\χ}\\ 0& \mathrm{if}\left|d\right|\≤\mathrm{\χ}\end{array}\right.;$

对于y′＝[y₁′，y₂′，...，y_n′]，也采用上述方法进行去噪处理；

②在每一采样时刻对船舶轨迹数据预处理，依据所获取的船舶原始离散二维位置序列x＝[x₁，x₂，...，x_n]和y＝[y₁，y₂，...，y_n]，采用一阶差分方法对其进行处理获取新的船舶离散位置序列Δx＝[Δx₁，Δx₂，...，Δx_n-1]和Δy＝[Δy₁，Δy₂，...，Δy_n-1]，其中Δx_i＝x_i+1-x_i，Δy_i＝y_i+1-y_i(i＝1，2，...，n-1)；

③在每一采样时刻对船舶轨迹数据聚类，对处理后新的船舶离散二维位置序列Δx和Δy，通过设定聚类个数M′，采用遗传聚类算法分别对其进行聚类；

④在每一采样时刻对船舶轨迹数据利用隐马尔科夫模型进行参数训练，通过将处理后的船舶运行轨迹数据Δx和Δy视为隐马尔科夫过程的显观测值，通过设定隐状态数目N和参数更新时段τ′，依据最近的T′个位置观测值并采用B-W算法滚动获取最新隐马尔科夫模型参数λ′；确定航迹隐马尔科夫模型参数λ′＝(π，A，B)的过程如下：

4.1)变量赋初值：应用均匀分布给变量π_i，a_ij和b_j(o_k)赋初值和并使其满足约束条件： ${\mathrm{\Σ}}_{i=1}^N{\mathrm{\π}}_i^0=1,{\mathrm{\Σ}}_{j=1}^N{a}_{\mathrm{ij}}^0=1(1\≤i\≤N)$ 阳 ${\mathrm{\Σ}}_{k=1}^{M^{\′}}{b}_j^0\left(o_k\right)=1(1\≤j\≤N),$ 由此得到λ₀＝(π₀，A₀，B₀)，其中o_k表示某一显观测值，π₀、A₀和B₀分别是由元素阳构成的矩阵，令参数l＝0，o＝(o_t-T′+1，...，o_t-1，o_t)为当前时刻t之前的T′个历史位置观测值；

4.2)执行E-M算法：

4.2.1)E-步骤：由λ_l计算ξ_e(i，j)和γ_e(s_i)；

变量 ${\mathrm{\ξ}}_e(i,j)=\frac{P(q_e=s_i,q_{e+1}=s_j,o|{\mathrm{\λ}}_l)}{P\left(o\right|{\mathrm{\λ}}_l)},$ 那么 ${\mathrm{\γ}}_e\left(s_i\right)={\mathrm{\Σ}}_{j=1}^N{\mathrm{\ξ}}_e(i,j),$

其中s表示某一隐状态；

4.2.2)M-步骤：运用 ${\stackrel{\‾}{\mathrm{\π}}}_i={\mathrm{\γ}}_1\left(s_i\right),{\stackrel{\‾}{a}}_{\mathrm{ij}}=\frac{{\mathrm{\Σ}}_{e=1}^{T^{\′}-1}{\mathrm{\ξ}}_e(i,j)}{{\mathrm{\Σ}}_{e=1}^{T^{\′}-1}{\mathrm{\γ}}_e\left(s_i\right)},{\stackrel{\‾}{b}}_j\left(o_k\right)=\frac{\underset{o_e=o_k}{{\mathrm{\Σ}}_{e=1}^{T^{\′}}{\mathrm{\γ}}_e\left(s_j\right)}}{{\mathrm{\Σ}}_{e=1}^{T^{\′}}{\mathrm{\γ}}_e\left(s_j\right)}$ 分别估计π_i，a_ij和b_j(o_k)并由此得到λ_l+1；

4.2.3)循环：l＝l+1，重复执行E-步骤和M-步骤，直至π_i、a_ij和b_j(o_k)收敛，即

|P(o|λ_l+1)-P(o|λ_l)|＜ε，其中参数ε＝0.00001，返回步骤4.2.4)；

4.2.4)：令λ′＝λ_l+1，算法结束。

⑤在每一采样时刻依据隐马尔科夫模型参数，采用Viterbi算法获取当前时刻观测值所对应的隐状态q：

5.1)变量赋初值：令g＝2，β_T′(s_i)＝1(s_i∈S)，δ₁(s_i)＝π_ib_i(o₁)，ψ₁(s_i)＝0，其中，

${\mathrm{\δ}}_g\left(s_i\right)=\underset{q_1,q_2,...,q_{g-1}}{\max }P(q_1,q_2,...,q_{g-1},q_g=s_i,o_1,o_2,...,o_g|{\mathrm{\λ}}^{\′})$ ，其中变量ψ_g(s_j)表示使变量δ_g-1(s_i)a_ij取最大值的船舶航迹隐状态s_i，参数S表示隐状态的集合；

5.2)递推过程： ${\mathrm{\δ}}_g\left(s_j\right)=\underset{s_i\∈S}{\max }\left[{\mathrm{\δ}}_{g-1}\left(s_i\right)a_{\mathrm{ij}}\right]b_j\left(o_g\right),{\mathrm{\ψ}}_g\left(s_j\right)=\underset{s_i\∈S}{\arg \max }\left[{\mathrm{\δ}}_{g-1}\left(s_i\right)a_{\mathrm{ij}}\right];$

5.3)时刻更新：令g＝g+1，若g≤T′，返回步骤5.2)，否则迭代终止并转到步骤5.4)；

5.4) $p^{*}=\underset{s_i\∈S}{\max }\left[{\mathrm{\δ}}_{T^{\′}}\left(s_i\right)\right],q_{T^{\′}}^{*}=\underset{s_i\∈S}{\arg \max }\left[{\mathrm{\δ}}_{T^{\′}}\left(s_i\right)\right],$ 转到步骤5.5)；

5.5)最优隐状态序列获取：

5.5.1)变量赋初值：令g＝T′-1；

5.5.2)后向递推：

5.5.3)时刻更新：令g＝g-1，若g≥1，返回步骤5.5.2)，否则终止。。

⑥在每一采样时刻，通过设定预测时域W，基于船舶当前时刻的隐状态q，获取未来时段船舶的位置预测值O。

上述聚类个数M′的值为4，，状态数目N的值为3，参数更新时段τ′为30秒，T′为10，预测时域W为300秒。

(应用例、航海交通管制方法)

本实施例的航海交通管制方法包括如下几个步骤：

步骤A、根据实施例1得到的船舶冲突预警方法获得船舶在每一采样时刻推测到的未来时段内船舶的轨迹；

步骤B、在每一采样时刻，基于船舶当前的运行状态和历史位置观察序列，获取海域风场变量的数值，其具体过程如下：

B.1)设定船舶的停靠位置为轨迹参考坐标原点；

B.2)在船舶处于直线运行状态和匀速转弯运行状态时，构建海域风场线性滤波模型；

B.3)根据所构建的滤波模型获取风场变量的数值。

步骤C、在每一采样时刻，基于各船舶的运行状态和设定的船舶在海域内运行时需满足的安全规则集，当船舶间有可能出现违反安全规则的状况时，对其动态行为实施监控并为海上交通控制中心提供及时的告警信息；

步骤D、当告警信息出现时，在满足船舶物理性能和海域交通规则的前提下，通过设定优化指标函数以及融入风场变量数值，采用模型预测控制理论方法对船舶避撞轨迹进行滚动规划，并将规划结果传输给各船舶执行，其具体过程如下：

D.1)设定船舶避撞轨迹规划的终止参考点位置P、避撞策略控制时域Θ、轨迹预测时域γ；

D.2)设定在给定优化指标函数的前提下，基于合作式避撞轨迹规划思想，通过给各个船舶赋予不同的权重以及融入实时风场变量滤波数值，得到各个船舶的避撞轨迹和避撞控制策略并将规划结果传输给各船舶执行，且各船舶在滚动规划间隔内仅实施其第一个优化控制策略；

D.3)在下一采样时刻，重复步骤D.2直至各船舶均到达其解脱终点。

上述终止参考点位置P设定为船舶位置冲突点的下一个航道点，避撞策略控制时域Θ为300秒；轨迹预测时域γ为300秒。

⑦见图2，在每一采样时刻，基于各船舶的运行状态和设定的船舶在海域内运行时需满足的安全规则集，当船舶间有可能出现违反安全规则的状况时，对其动态行为实施监控并为海上交通控制中心提供及时的告警信息，其具体过程如下：

7.1)构造船舶在海域内运行时需满足的安全规则集D_mr(t)≥D_min，其中D_mr(t)表示任意两个船舶m和船舶r在t时刻的距离，D_min表示船舶间的最小安全距离；

7.2)依据采样时间，建立由船舶连续运行状态至离散采样状态的观测器Λ：г→Ξ，其中г表示船舶的连续运行状态，Ξ表示船舶的离散采样状态；

7.3)当船舶m和r的观测器Λ_m和Λ_r的离散观测数值Ξ_m和Ξ_r在t时刻表明该向量不在安全规则集中时，即关系式D_mr(t)≥D_min不成立时，立刻向海上交通控制中心发出告警信息。

显然，上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。而这些属于本发明的精神所引伸出的显而易见的变化或变动仍处于本发明的保护范围之中。

Claims (8)

1.一种船舶冲突预警方法，其特征在于包括如下几个步骤：

①通过海面雷达获得船舶的实时和历史位置信息，各船舶的位置信息为离散二维位置序列x′＝[x₁′，x₂′，...，x_n′]和y′＝[y₁′，y₂′，...，y_n′]，通过应用小波变换理论对原始离散二维位置序列x′＝[x₁′，x₂′，...，x_n′]和y′＝[y₁′，y₂′，...，y_n′]进行初步处理，从而获取船舶的去噪离散二维位置序列x＝[x₁，x₂，...，x_n]和y＝[y₁，y₂，...，y_n]；

②在每一采样时刻对船舶轨迹数据预处理，依据所获取的船舶原始离散二维位置序列x＝[x₁，x₂，...，x_n]和y＝[y₁，y₂，...，y_n]，采用一阶差分方法对其进行处理获取新的船舶离散位置序列Δx＝[Δx₁，Δx₂，...，Δx_n-1]和Δy＝[Δy₁，Δy₂，...，Δy_n-1]，其中Δx_i＝x_i+1-x_i，Δy_i＝y_i+1-y_i(i＝1，2，...，n-1)；

③在每一采样时刻对船舶轨迹数据聚类，对处理后新的船舶离散二维位置序列Δx和Δy，通过设定聚类个数M′，采用遗传聚类算法分别对其进行聚类；

④在每一采样时刻对船舶轨迹数据利用隐马尔科夫模型进行参数训练，通过将处理后的船舶运行轨迹数据Δx和Δy视为隐马尔科夫过程的显观测值，通过设定隐状态数目N和参数更新时段τ′，依据最近的T′个位置观测值并采用B-W算法滚动获取最新隐马尔科夫模型参数λ′；

⑤在每一采样时刻依据隐马尔科夫模型参数，采用Viterbi算法获取当前时刻观测值所对应的隐状态q；

⑥在每一采样时刻，通过设定预测时域W，基于船舶当前时刻的隐状态q，获取未来时段船舶的位置预测值O，从而在每一采样时刻滚动推测到未来时段内船舶的轨迹；

⑦在每一采样时刻，基于各船舶的运行状态和设定的船舶在海域内运行时需满足的安全规则集，当船舶间有可能出现违反安全规则的状况时，对其动态行为实施监控并为海上交通控制中心提供及时的告警信息。

2.根据权利要求1所述的船舶冲突预警方法，其特征在于：所述步骤①中，通过应用小波变换理论对原始离散二维位置序列x′＝[x₁′，x₂′，...，x_n′]和y′＝[y₁′，y₂′，...，y_n′]进行初步处理，从而获取船舶的去噪离散二维位置序列x＝[x₁，x₂，...，x_n]和y＝[y₁，y₂，...，y_n]：对于给定的原始二维序列数据x′＝[x₁′，x₂′，...，x_n′]，利用如下形式的线性表达式分别对其进行近似： $\tilde{f}\left(x^{\′}\right)=\underset{\underset{\∀K}{\∀J}}{\mathrm{\Σ}}{c}_{J,K}{\mathrm{\ψ}}_{J,K}\left(x^{\′}\right),$

其中： $c_{J,K}={\∫}_{-\∞}^{\∞}{f}^{\′}\left(x^{\′}\right){\mathrm{\ψ}}_{J,K}\left(x^{\′}\right){\mathrm{dx}}^{\′},$

ψ_J，K(x′)＝δ·ψ(2^Jx′-K)

f′(x′)表示对数据平滑处理后得到的函数表达式，ψ(x′)表示母波，δ、J和K均为小波变换常数，ψ_J，K(x′)表示母波的转换形式，c_J，K表示由小波变换过程得到的函数系数，它体现了子波ψ_J，K(x′)对整个函数近似的权重大小，若此系数很小，那么它意味着子波ψ_J，K(x′)的权重也较小，因而可以在不影响函数主要特性的前提下，从函数近似过程中将子波ψ_J，K(x′)除去；在实际数据处理过程中，通过设定阈值χ来实施“阈值转换”，当c_J，K＜χ时，设定c_J，K＝0；阈值函数的选取采用如下两种方式：

${\mathrm{\ρ}}_1(d,\mathrm{\χ})=\left\{\begin{array}{cc}d& \mathrm{if}\left|d\right|>\mathrm{\χ}\\ 0& \mathrm{if}\left|d\right|\≤\mathrm{\χ}\end{array}\right.$ 和 ${\mathrm{\ρ}}_2(d,\mathrm{\χ})=\left\{\begin{array}{cc}d-\left(\frac{d}{\left|d\right|}\mathrm{\χ}\right)& \mathrm{if}\left|d\right|>\mathrm{\χ}\\ 0& \mathrm{if}\left|d\right|\≤\mathrm{\χ}\end{array}\right.;$

对于y′＝[y₁′，y₂′，...，y_n′]，也采用上述方法进行去噪处理。

3.根据权利要求1或2所述的船舶冲突预警方法，其特征在于：所述步骤④中确定航迹隐马尔科夫模型参数λ′＝(π，A，B)的过程如下：

4.1)变量赋初值：应用均匀分布给变量π_i，a_ij和b_j(o_k)赋初值和并使其满足约束条件： (1≤i≤N)和(1≤j≤N)，由此得到λ₀＝(π₀，A₀，B₀)，其中o_k表示某一显观测值，π₀、A₀和B₀分别是由元素和构成的矩阵，令参数l＝0，o＝(o_t-T′+1，...，o_t-1，o_t)为当前时刻t之前的T′个历史位置观测值；

4.2)执行E-M算法：

4.2.1)E-步骤：由λ_l计算ξ_e(i，j)和γ_e(s_i)；

变量 ${\mathrm{\ξ}}_e(i,j)=\frac{P(q_e=s_i,q_{e+1}=s_j,o|{\mathrm{\λ}}_l)}{P\left(o\right|{\mathrm{\λ}}_l)},$ 那么 ${\mathrm{\γ}}_e\left(s_i\right)={\mathrm{\Σ}}_{j=1}^N{\mathrm{\ξ}}_e(i,j),$

其中s表示某一隐状态；

4.2.2)M-步骤：运用 ${\stackrel{\‾}{\mathrm{\π}}}_i={\mathrm{\γ}}_1\left(s_i\right),$ ${\stackrel{\‾}{a}}_{\mathrm{ij}}=\frac{{\mathrm{\Σ}}_{e=1}^{T^{\′}-1}{\mathrm{\ξ}}_e(i,j)}{{\mathrm{\Σ}}_{e=1}^{T^{\′}-1}{\mathrm{\γ}}_e\left(s_i\right)},$ ${\stackrel{\‾}{b}}_j\left(o_k\right)=\frac{\underset{o_e=o_k}{{\mathrm{\Σ}}_{e=1}^{T^{\′}}{\mathrm{\γ}}_e\left(s_j\right)}}{{\mathrm{\Σ}}_{e=1}^{T^{\′}}{\mathrm{\γ}}_e\left(s_j\right)}$ 分别估计π_i，a_ij和b_j(o_k)并由此得到λ_l+1；

4.2.3)循环：l＝l+1，重复执行E-步骤和M-步骤，直至π_i、a_ij和b_j(o_k)收敛，即其中参数ε＝0.00001，返回步骤4.2.4)；

4.2.4)：令λ′＝λ_l+1，算法结束。

4.根据权利要求1至3之一所述的船舶冲突预警方法，其特征在于：所述步骤⑤中确定船舶航迹最佳隐状态序列的迭代过程如下：

5.1)变量赋初值：令g＝2，β_T′(s_i)＝1(s_i∈S)，δ₁(s_i)＝π_ib_i(o₁)，ψ₁(s_i)＝0，其中， ${\mathrm{\δ}}_g\left(s_i\right)=\underset{q_1,q_2,...,q_{g-1}}{\max }P(q_1,q_2,...,q_{g-1},q_g=s_i,o_1,o_2,...,o_g|{\mathrm{\λ}}^{\′}),$

其中变量ψ_g(s_j)表示使变量δ_g-1(s_i)a_ij取最大值的船舶航迹隐状态s_i，参数S表示隐状态的集合；

5.2)递推过程： ${\mathrm{\δ}}_g\left(s_j\right)=\underset{s_i\∈S}{\max }\left[{\mathrm{\δ}}_{g-1}\left(s_i\right)a_{\mathrm{ij}}\right]b_j\left(o_g\right),$ ${\mathrm{\ψ}}_g\left(s_j\right)=\underset{s_i\∈S}{\arg \max }\left[{\mathrm{\δ}}_{g-1}\left(s_i\right)a_{\mathrm{ij}}\right];$

5.3)时刻更新：令g＝g+1，若g≤T′，返回步骤5.2)，否则迭代终止并转到步骤5.4)；

5.4) $p^{*}=\underset{s_i\∈S}{\max }\left[{\mathrm{\δ}}_{T^{\′}}\left(s_i\right)\right],$ $q_{T^{\′}}^{*}=\underset{s_i\∈S}{\arg \max }\left[{\mathrm{\δ}}_{T^{\′}}\left(s_i\right)\right],$ 转到步骤5.5)；

5.5)最优隐状态序列获取：

5.5.1)变量赋初值：令g＝T′-1；

5.5.2)后向递推：

5.5.3)时刻更新：令g＝g-1，若g≥1，返回步骤5.5.2)，否则终止。

5.根据权利要求1至4之一所述的船舶冲突预警方法，其特征在于：所述步骤③中，聚类个数M′的值为4。

6.根据权利要求1至5之一所述的船舶冲突预警方法，其特征在于：所述步骤④中，状态数目N的值为3，参数更新时段τ′为30秒，T′为10。

7.根据权利要求1至6之一所述的船舶冲突预警方法，其特征在于：所述步骤⑥中，预测时域W为300秒。

8.根据权利要求1至7之一所述的航海交通管制方法，其特征在于：所述步骤⑦中对各船舶的动态行为实施监控并为海上交通控制中心提供及时的告警信息的具体过程如下：

7.1)构造船舶在海域内运行时需满足的安全规则集D_mr(t)≥D_min，其中D_mr(t)表示任意两个船舶m和船舶r在t时刻的距离，D_min表示船舶间的最小安全距离；

7.2)依据采样时间，建立由船舶连续运行状态至离散采样状态的观测器Λ：Γ→Ξ，其中Γ表示船舶的连续运行状态，Ξ表示船舶的离散采样状态；

7.3)当船舶m和r的观测器Λ_m和Λ_r的离散观测数值Ξ_m和Ξ_r在t时刻表明该向量不在安全规则集中时，即关系式D_mr(t)≥D_min不成立时，立刻向海上交通控制中心发出告警信息。