CN104504935A - 航海交通管制方法 - Google Patents

Abstract

本发明涉及一种航海交通管制方法，包括如下几个步骤，首先通过海面雷达获得船舶的实时和历史位置信息；然后在每一采样时刻，依据船舶的实时和历史位置信息滚动推测未来时段内船舶的轨迹；再而基于船舶当前的运行状态和历史位置观察序列，获取海域风场变量的数值；再而基于各船舶的运行状态和设定的船舶在海域内运行时需满足的安全规则集，对船舶的动态行为实施监控并为海上交通控制中心提供及时的告警信息；当告警信息出现时，在满足船舶物理性能和海域交通规则的前提下，通过设定优化指标函数以及融入风场变量数值，采用自适应控制理论方法对船舶避撞轨迹进行滚动规划，并将规划结果传输给各船舶执行。本发明实时预测规划轨迹，安全性较好。

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]；

②在每一采样时刻，依据步骤①得到的船舶的实时和历史位置信息滚动推测未来时段内船舶的轨迹，其具体过程如下：

2.1)船舶轨迹数据预处理，依据所获取的船舶原始离散二维位置序列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)；

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

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

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

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

③在每一采样时刻，基于船舶当前的运行状态和历史位置观察序列，获取海域风场变量的数值；

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

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

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

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

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

进一步的，所述步骤①中，通过应用小波变换理论对原始离散二维位置序列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)d{x}^{\′}\\ {\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′]，也采用上述方法进行去噪处理。

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

2.3.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\≤i\≤N),$ 由此得到λ₀＝(π₀，A₀，B₀)，其中o_k表示某一显观测值，π₀、A₀和B₀分别是由元素和构成的矩阵，令参数l＝0，o＝(o_t-T′+1，...，o_t-1，o_t)为当前时刻t之前的T′个历史位置观测值；

2.3.2)执行E-M算法：

2.3.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表示某一隐状态；

2.3.2.2)M-步骤：运用 ${\stackrel{\‾}{\mathrm{\π}}}_i={\mathrm{\γ}}_i\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；

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

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

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

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

2.4.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表示隐状态的集合；

2.4.2)递推过程： ${\mathrm{\δ}}_g\left(s_i\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];$

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

2.4.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],$ 转到步骤2.4.5)；

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

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

2.4.5.2)后向递推：

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

进一步的，所述步骤②中，聚类个数M′的值为4，隐状态数目N的值为3，参数更新时段τ′为30秒，T′为10，预测时域W为300秒。

进一步的，所述步骤③获取海域风场变量的数值的具体过程如下：

3.1)设定船舶的停靠位置为轨迹参考坐标原点并在水平面上建立横坐标轴和纵坐标轴；

3.2)在船舶处于直线运行状态和匀速转弯运行状态时，构建海域风场线性滤波模型x₁(t+Δt)＝F(t)x₁(t)+w(t)和z(t)＝H(t)x₁(t)+v(t)获取风场变量数值，其中Δt表示采样间隔，x₁(t)表示t时刻的状态向量，z(t)表示t时刻的观测向量，且x₁(t)＝[x(t)，y(t)，v_x(t)，v_y(t)，w_x(t)，w_y(t)]^T，其中x(t)和y(t)分别表示t时刻船舶位置在横坐标轴和纵坐标轴上的分量，v_x(t)和v_y(t)分别表示t时刻船舶速度在横坐标轴和纵坐标轴上的分量，w_x(t)和w_y(t)分别表示t时刻风场数值在横坐标轴和纵坐标轴上的分量，F(t)和H(t)分别表示状态转移矩阵和输出测量矩阵，w(t)和v(t)分别表示系统噪声向量和测量噪声向量：

$F\left(t\right)=\left[\begin{array}{cccccc}1& 0& \frac{\sin \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a\left(t\right)}& \frac{1-\cos \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a\left(t\right)}& \mathrm{\Δt}& 0\\ 0& 1& \frac{\cos \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)-1}{{\mathrm{\ω}}_a\left(t\right)}& \frac{\sin \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a\left(t\right)}& 0& \mathrm{\Δt}\\ 0& 0& \cos \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)& \sin \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)& 0& 0\\ 0& 0& -\sin \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)& \cos \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right]$

$H\left(k\right)=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\end{array}\right];$

在船舶处于变速转弯运行状态时，构建海域风场非线性滤波模型x₁(t+Δt)＝Ψ(t，x₁(t)，u(t))+w(t)、z(t)＝Ω(t，x₁(t))+v(t)和u(t)＝[ω_a(t)，γ_a(t)]^T，其中Ψ(·)和Ω(·)分别表示状态转移矩阵和输出测量矩阵，ω_a(t)和γ_a(t)分别表示转弯率和加速率：

$\mathrm{\Ψ}=\left[\begin{array}{c}x\left(t\right)+v_x\left(t\right)(\frac{\sin \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a\left(t\right)}+\frac{{\mathrm{\γ}}_a\left(t\right)\mathrm{\Δt}}{\sqrt{v_x^2\left(t\right)+v_y^2\left(t\right)}}{C}_5)+v_y\left(t\right)(\frac{1-\cos \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a\left(t\right)}+\frac{{\mathrm{\γ}}_a\left(t\right)\mathrm{\Δt}}{\sqrt{v_x^2\left(t\right)+v_y^2\left(t\right)}}{C}_6)+w_x\left(t\right)\\ y\left(t\right)-v_x\left(t\right)(\frac{1-\cos \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a\left(t\right)}+\frac{{\mathrm{\γ}}_a\left(t\right)\mathrm{\Δt}}{\sqrt{v_x^2\left(t\right)+v_y^2\left(t\right)}}{C}_6)+v_y\left(t\right)(\frac{\sin \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a\left(t\right)}+\frac{{\mathrm{\γ}}_a\left(t\right)\mathrm{\Δt}}{\sqrt{v_x^2\left(t\right)+v_y^2\left(t\right)}}{C}_5)+w_y\left(t\right)\\ \left((1+\frac{{\mathrm{\γ}}_a\left(t\right)\mathrm{\Δt}}{\sqrt{v_x^2\left(t\right)+v_y^2\left(t\right)}})(v_x\left(t\right)\cos \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)+v_y\left(t\right)\sin \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right))\right)\\ ((1+\frac{{\mathrm{\γ}}_a\left(t\right)\mathrm{\Δt}}{\sqrt{v_x^2\left(t\right)+v_y^2\left(t\right)}})(v_y\left(t\right)\cos \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)-v_x\left(t\right)\sin \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)))\\ w_x\left(t\right)\\ w_y\left(t\right)\end{array}\right],$

其中：Δt表示采样时间间隔，

$C_5=(\frac{\sin \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a\left(t\right)}-\frac{1-\cos \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a^2\left(t\right)\mathrm{\Δt}}),$

$C_6=(\frac{\sin \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a^2\left(t\right)\mathrm{\Δt}}-\frac{\cos \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a\left(t\right)});$

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

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

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

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

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

进一步的，步骤⑤中，步骤5.2)的具体过程是：令

$d_{\mathrm{Rt}}^2={\left|\right|P_R\left(t\right)-P_R^f\left|\right|}_2^2={(x_{\mathrm{Rt}}-x_R^f)}^2+{(y_{\mathrm{Rt}}-y_R^f)}^2,$

其中表示t时刻船舶R当前所在位置和下一航道点间的距离的平方，P_R(t)＝(x_Rt，y_Rt)，那么t时刻船舶R的优先级指数可设定为：

$L_{\mathrm{Rt}}=100\frac{d_{\mathrm{Rt}}^{-2}}{\underset{R=1}{\overset{Z_t}{\mathrm{\Σ}}}{d}_{\mathrm{Rt}}^{-2}},$

其中z_t表示t时刻海域内存在冲突的船舶数目，由优先级指数的含义可知，船舶距离其下一航道点越近，其优先级越高；

设定优化指标

$\begin{array}{c}{\mathrm{\Φ}}^{*}(u_1^{\left(t\right)},u_1^{(t+\mathrm{\Δt})},...,u_1^{(t+\mathrm{p\Δt})},...,u_{Z_t}^{\left(t\right)},u_{Z_t}^{(t+\mathrm{\Δt})},...,u_{Z_t}^{(t+\mathrm{p\Δt})})=\underset{h=1}{\overset{p}{\mathrm{\Σ}}}\underset{R=1}{\overset{Z_t}{\mathrm{\Σ}}}{L}_{\mathrm{Rt}}{\left|\right|P_R(t+\mathrm{h\Δt})-P_R^f\left|\right|}_2^2\\ =\underset{h=1}{\overset{p}{\mathrm{\Σ}}}\underset{R=1}{\overset{Z_t}{\mathrm{\Σ}}}{(P_R(t+\mathrm{h\Δt})-P_R^f)}^T{Q}_{\mathrm{Rt}}(P_R(t+\mathrm{h\Δt})-P_R^f)\end{array}$

，其中R∈I(t)表示船舶代码且I(t)＝{1，2，...，Zt}，P_R(t+hΔt)表示船舶在时刻(t+hΔt)的位置向量，表示船舶R的解脱终止点，u_R表示待优化的船舶R的最优控制序列，Q_Rt为正定对角矩阵，其对角元素为船舶R在t时刻的优先级指数L_Rt，并且 $Q_{\mathrm{Rt}}=\left[\begin{array}{cc}{L}_{\mathrm{Rt}}& 0\\ 0& L_{\mathrm{Rt}}\end{array}\right].$

进一步的，所述步骤⑤中终止参考点位置P设定为船舶运行的下一个航道点，避撞策略控制时域Θ为300秒；轨迹预测时域W为300秒。

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

(2)本发明在船舶冲突解脱过程中，融入了海域内风场的影响，所采用的滚动解脱轨迹规划方案能够根据海域内风场的变化及时调整解脱轨迹，提高了船舶冲突解脱的鲁棒性。

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

附图说明

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

图2为本发明中的风场滤波方法流程示意图；

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

图4为本发明中的船舶避撞轨迹优化方法流程示意图。

具体实施方式

(实施例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)d{x}^{\′}\\ {\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′]，也采用上述方法进行去噪处理。

②在每一采样时刻，依据步骤①得到的船舶的实时和历史位置信息滚动推测未来时段内船舶的轨迹，见图1，其具体过程如下：

2.1)船舶轨迹数据预处理，依据所获取的船舶原始离散二维位置序列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)；

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

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

2.3.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\≤i\≤N),$ 由此得到λ₀＝(π₀，A₀，B₀)，其中o_k表示某一显观测值，π₀、A₀和B₀分别是由元素和构成的矩阵，令参数l＝0，o＝(o_t-T′+1，...，o_t-1，o_t)为当前时刻t之前的T′个历史位置观测值；

2.3.2)执行E-M算法：

2.3.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表示某一隐状态；

2.3.2.2)M-步骤：运用 ${\stackrel{\‾}{\mathrm{\π}}}_i={\mathrm{\γ}}_i\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；

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

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

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

2.4)依据隐马尔科夫模型参数，采用Viterbi算法获取当前时刻观测值所对应的隐状态q；确定船舶航迹最佳隐状态序列的迭代过程如下：

2.4.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表示隐状态的集合；

2.4.2)递推过程： ${\mathrm{\δ}}_g\left(s_i\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];$

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

2.4.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],$ 转到步骤2.4.5)；

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

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

2.4.5.2)后向递推：

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

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

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

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

3.1)设定船舶的停靠位置为轨迹参考坐标原点并在水平面上建立横坐标轴和纵坐标轴；

3.2)在船舶处于直线运行状态和匀速转弯运行状态时，构建海域风场线性滤波模型x₁(t+Δt)＝F(t)x₁(t)+w(t)和z(t)＝H(t)x₁(t)+v(t)获取风场变量数值，其中Δt表示采样间隔，x₁(t)表示t时刻的状态向量，z(t)表示t时刻的观测向量，且x₁(t)＝[x(t)，y(t)，v_x(t)，v_y(t)，w_x(t)，w_y(t)]^T，其中x(t)和y(t)分别表示t时刻船舶位置在横坐标轴和纵坐标轴上的分量，v_x(t)和v_y(t)分别表示t时刻船舶速度在横坐标轴和纵坐标轴上的分量，w_x(t)和w_y(t)分别表示t时刻风场数值在横坐标轴和纵坐标轴上的分量，F(t)和H(t)分别表示状态转移矩阵和输出测量矩阵，w(t)和v(t)分别表示系统噪声向量和测量噪声向量：

$F\left(t\right)=\left[\begin{array}{cccccc}1& 0& \frac{\sin \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a\left(t\right)}& \frac{1-\cos \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a\left(t\right)}& \mathrm{\Δt}& 0\\ 0& 1& \frac{\cos \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)-1}{{\mathrm{\ω}}_a\left(t\right)}& \frac{\sin \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a\left(t\right)}& 0& \mathrm{\Δt}\\ 0& 0& \cos \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)& \sin \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)& 0& 0\\ 0& 0& -\sin \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)& \cos \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right]$

$H\left(k\right)=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\end{array}\right];$

在船舶处于变速转弯运行状态时，构建海域风场非线性滤波模型x₁(t+Δt)＝Ψ(t，x₁(t)，u(t))+w(t)、z(t)＝Ω(t，x₁(t))+v(t)和u(t)＝[ω_a(t)，γ_a(t)]^T，其中Ψ(·)和Ω(·)分别表示状态转移矩阵和输出测量矩阵，ω_a(t)和γ_a(t)分别表示转弯率和加速率：

$\mathrm{\Ψ}=\left[\begin{array}{c}x\left(t\right)+v_x\left(t\right)(\frac{\sin \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a\left(t\right)}+\frac{{\mathrm{\γ}}_a\left(t\right)\mathrm{\Δt}}{\sqrt{v_x^2\left(t\right)+v_y^2\left(t\right)}}{C}_5)+v_y\left(t\right)(\frac{1-\cos \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a\left(t\right)}+\frac{{\mathrm{\γ}}_a\left(t\right)\mathrm{\Δt}}{\sqrt{v_x^2\left(t\right)+v_y^2\left(t\right)}}{C}_6)+w_x\left(t\right)\\ y\left(t\right)-v_x\left(t\right)(\frac{1-\cos \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a\left(t\right)}+\frac{{\mathrm{\γ}}_a\left(t\right)\mathrm{\Δt}}{\sqrt{v_x^2\left(t\right)+v_y^2\left(t\right)}}{C}_6)+v_y\left(t\right)(\frac{\sin \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a\left(t\right)}+\frac{{\mathrm{\γ}}_a\left(t\right)\mathrm{\Δt}}{\sqrt{v_x^2\left(t\right)+v_y^2\left(t\right)}}{C}_5)+w_y\left(t\right)\\ \left((1+\frac{{\mathrm{\γ}}_a\left(t\right)\mathrm{\Δt}}{\sqrt{v_x^2\left(t\right)+v_y^2\left(t\right)}})(v_x\left(t\right)\cos \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)+v_y\left(t\right)\sin \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right))\right)\\ ((1+\frac{{\mathrm{\γ}}_a\left(t\right)\mathrm{\Δt}}{\sqrt{v_x^2\left(t\right)+v_y^2\left(t\right)}})(v_y\left(t\right)\cos \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)-v_x\left(t\right)\sin \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)))\\ w_x\left(t\right)\\ w_y\left(t\right)\end{array}\right],$

其中：Δt表示采样时间间隔，

$C_5=(\frac{\sin \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a\left(t\right)}-\frac{1-\cos \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a^2\left(t\right)\mathrm{\Δt}}),$

$C_6=(\frac{\sin \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a^2\left(t\right)\mathrm{\Δt}}-\frac{\cos \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a\left(t\right)});$

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

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

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

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

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

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

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

5.2)设定在给定优化指标函数的前提下，基于合作式避撞轨迹规划思想，通过给各个船舶赋予不同的权重以及融入实时风场变量滤波数值，得到各个船舶的避撞轨迹和避撞控制策略并将规划结果传输给各船舶执行，且各船舶在滚动规划间隔内仅实施其第一个优化控制策略：令 $d_{\mathrm{Rt}}^2={\left|\right|P_R\left(t\right)-P_R^f\left|\right|}_2^2={(x_{\mathrm{Rt}}-x_R^f)}^2+{(y_{\mathrm{Rt}}-y_R^f)}^2,$

其中表示t时刻船舶R当前所在位置和下一航道点间的距离的平方，P_R(t)＝(x_Rt，y_Rt)，那么t时刻船舶R的优先级指数可设定为：

$L_{\mathrm{Rt}}=100\frac{d_{\mathrm{Rt}}^{-2}}{\underset{R=1}{\overset{Z_t}{\mathrm{\Σ}}}{d}_{\mathrm{Rt}}^{-2}},$

其中z_t表示t时刻海域内存在冲突的船舶数目，由优先级指数的含义可知，船舶距离其下一航道点越近，其优先级越高；

设定优化指标

$\begin{array}{c}{\mathrm{\Φ}}^{*}(u_1^{\left(t\right)},u_1^{(t+\mathrm{\Δt})},...,u_1^{(t+\mathrm{p\Δt})},...,u_{Z_t}^{\left(t\right)},u_{Z_t}^{(t+\mathrm{\Δt})},...,u_{Z_t}^{(t+\mathrm{p\Δt})})=\underset{h=1}{\overset{p}{\mathrm{\Σ}}}\underset{R=1}{\overset{Z_t}{\mathrm{\Σ}}}{L}_{\mathrm{Rt}}{\left|\right|P_R(t+\mathrm{h\Δt})-P_R^f\left|\right|}_2^2\\ =\underset{h=1}{\overset{p}{\mathrm{\Σ}}}\underset{R=1}{\overset{Z_t}{\mathrm{\Σ}}}{(P_R(t+\mathrm{h\Δt})-P_R^f)}^T{Q}_{\mathrm{Rt}}(P_R(t+\mathrm{h\Δt})-P_R^f)\end{array}$

，其中R∈I(t)表示船舶代码且I(t)＝{1，2，...，Z_t}，P_R(t+hΔt)表示船舶在时刻(t+hΔt)的位置向量，表示船舶R的解脱终止点，u_R表示待优化的船舶R的最优控制序列，Q_Rt为正定对角矩阵，其对角元素为船舶R在t时刻的优先级指数L_Rt，并且

$Q_{\mathrm{Rt}}=\left[\begin{array}{cc}{L}_{\mathrm{Rt}}& 0\\ 0& L_{\mathrm{Rt}}\end{array}\right];$

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

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

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

Claims (9)

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]；

②在每一采样时刻，依据步骤①得到的船舶的实时和历史位置信息滚动推测未来时段内船舶的轨迹，其具体过程如下：

2.1)船舶轨迹数据预处理，依据所获取的船舶原始离散二维位置序列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)；

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

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

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

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

③在每一采样时刻，基于船舶当前的运行状态和历史位置观察序列，获取海域风场变量的数值；

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

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

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

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

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

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{\∀K}{\underset{\∀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′]，也采用上述方法进行去噪处理。

3.根据权利要求1或2所述的航海交通管制方法，其特征在于：所述步骤②中，步骤2.3)中确定航迹隐马尔科夫模型参数λ′＝(π，A，B)的过程如下：

2.3.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′个历史位置观测值；

2.3.2)执行E-M算法：

2.3.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表示某一隐状态；

2.3.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；

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

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

4.根据权利要求1至3之一所述的航海交通管制方法，其特征在于：所述步骤②中，步骤2.4)确定船舶航迹最佳隐状态序列的迭代过程如下：

2.4.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表示隐状态的集合；

2.4.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];$

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

2.4.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],$ 转到步骤2.4.5)；

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

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

2.4.5.2)后向递推：

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

5.根据权利要求1至4之一所述的航海交通管制方法，其特征在于：所述步骤②中，聚类个数M′的值为4，隐状态数目N的值为3，参数更新时段τ′为30秒，T′为10，预测时域W为300秒。

6.根据权利要求1至5之一所述的航海交通管制方法，其特征在于：所述步骤③获取海域风场变量的数值的具体过程如下：

3.1)设定船舶的停靠位置为轨迹参考坐标原点并在水平面上建立横坐标轴和纵坐标轴；

3.2)在船舶处于直线运行状态和匀速转弯运行状态时，构建海域风场线性滤波模型x₁(t+Δt)＝F(t)x₁(t)+w(t)和z(t)＝H(t)x₁(t)+v(t)获取风场变量数值，其中Δt表示采样间隔，x₁(t)表示t时刻的状态向量，z(t)表示t时刻的观测向量，且x₁(t)＝[x(t)，y(t)，v_x(t)，v_y(t)，w_x(t)，w_y(t)]^T，其中x(t)和y(t)分别表示t时刻船舶位置在横坐标轴和纵坐标轴上的分量，v_x(t)和v_y(t)分别表示t时刻船舶速度在横坐标轴和纵坐标轴上的分量，w_x(t)和w_y(t)分别表示t时刻风场数值在横坐标轴和纵坐标轴上的分量，F(t)和H(t)分别表示状态转移矩阵和输出测量矩阵，w(t)和v(t)分别表示系统噪声向量和测量噪声向量：

$F\left(t\right)=\left[\begin{array}{cccccc}1& 0& \frac{\sin \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a\left(t\right)}& \frac{1-\cos \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a\left(t\right)}& \mathrm{\Δt}& 0\\ 0& 1& \frac{\cos \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)-1}{{\mathrm{\ω}}_a\left(t\right)}& \frac{\sin \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a\left(t\right)}& 0& \mathrm{\Δt}\\ 0& 0& \cos \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)& \sin \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)& 0& 0\\ 0& 0& -\sin \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)& \cos \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right]$

$H\left(k\right)=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\end{array}\right];$

在船舶处于变速转弯运行状态时，构建海域风场非线性滤波模型x₁(t+Δt)＝Ψ(t，x₁(t)，u(t))+w(t)、z(t)＝Ω(t，x₁(t))+v(t)和u(t)＝[ω_a(t)，γ_a(t)]T，其中Ψ(·)和Ω(·)分别表示状态转移矩阵和输出测量矩阵，ω_a(t)和γ_a(t)分别表示转弯率和加速率：

$\mathrm{\Ψ}=\left[\begin{array}{c}x\left(t\right)+v_x\left(t\right)(\frac{\sin \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a\left(t\right)}+\frac{{\mathrm{\γ}}_a\left(t\right)\mathrm{\Δt}}{\sqrt{v_x^2\left(t\right)+v_y^2\left(t\right)}}{C}_5)+v_y\left(t\right)(\frac{1-\cos \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a\left(t\right)}+\frac{{\mathrm{\γ}}_a\left(t\right)\mathrm{\Δt}}{\sqrt{v_x^2\left(t\right)+v_y^2\left(t\right)}}{C}_6)+w_x\left(t\right)\\ y\left(t\right)-v_x\left(t\right)(\frac{1-\cos \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a\left(t\right)}+\frac{{\mathrm{\γ}}_a\left(t\right)\mathrm{\Δt}}{\sqrt{v_x^2\left(t\right)+v_y^2\left(t\right)}}{C}_6)+v_y\left(t\right)(\frac{\sin \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a\left(t\right)}+\frac{{\mathrm{\γ}}_a\left(t\right)\mathrm{\Δt}}{\sqrt{v_x^2\left(t\right)+v_y^2\left(t\right)}}{C}_5)+w_y\left(t\right)\\ \left((1+\frac{{\mathrm{\γ}}_a\left(t\right)\mathrm{\Δt}}{\sqrt{v_x^2\left(t\right)+v_y^2\left(t\right)}})(v_x\left(t\right)\cos \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)+v_y\left(t\right)\sin \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right))\right)\\ \left((1+\frac{{\mathrm{\γ}}_a\left(t\right)\mathrm{\Δt}}{\sqrt{v_x^2\left(t\right)+v_y^2\left(t\right)}})(v_y\left(t\right)\cos \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)-v_x\left(t\right)\sin \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right))\right)\\ w_x\left(t\right)\\ w_y\left(t\right)\end{array}\right],$

其中：Δt表示采样时间间隔，

$C_5=(\frac{\sin \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a\left(t\right)}-\frac{1-\cos \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a^2\left(t\right)\mathrm{\Δt}}),$

$C_6=(\frac{\sin \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a^2\left(t\right)\mathrm{\Δt}}-\frac{\cos \left({\mathrm{\ω}}_a\left(t\right)\mathrm{\Δt}\right)}{{\mathrm{\ω}}_a\left(t\right)});$

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

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

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

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

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

8.根据权利要求1至7之一所述的航海交通管制方法，其特征在于：步骤⑤中，步骤5.2)的具体过程是：令 $d_{\mathrm{Rt}}^2={\left|\right|P_R\left(t\right)-P_R^f\left|\right|}_2^2={(x_{\mathrm{Rt}}-x_R^f)}^2+{(y_{\mathrm{Rt}}-y_R^f)}^2,$

其中表示t时刻船舶R当前所在位置和下一航道点间的距离的平方，P_R(t)＝(x_Rt，y_Rt)，那么t时刻船舶R的优先级指数可设定为：

$L_{\mathrm{Rt}}=100\frac{d_{\mathrm{Rt}}^2}{\underset{R=1}{\overset{Z_t}{\mathrm{\Σ}}}{d}_{\mathrm{Rt}}^2},$

其中z_t表示t时刻海域内存在冲突的船舶数目，由优先级指数的含义可知，船舶距离其下一航道点越近，其优先级越高；

设定优化指标

$\begin{array}{c}{\mathrm{\Φ}}^{*}(u_1^{\left(t\right)},u_1^{(t+\mathrm{\Δt})},...,u_1^{(t+\mathrm{p\Δt})},...,u_{Z_t}^{\left(t\right)},u_{Z_t}^{(t+\mathrm{\Δt})},...,u_{Z_t}^{(t+\mathrm{p\Δt})})=\underset{h=1}{\overset{p}{\mathrm{\Σ}}}\underset{R=1}{\overset{Z_t}{\mathrm{\Σ}}}{L}_{\mathrm{Rt}}{\left|\right|P_R(t+\mathrm{h\Δt})-P_R^f\left|\right|}_2^2\\ =\underset{h=1}{\overset{p}{\mathrm{\Σ}}}\underset{R=1}{\overset{Z_t}{\mathrm{\Σ}}}{(P_R(t+\mathrm{h\Δt})-P_R^f)}^T{Q}_{\mathrm{Rt}}(P_R(t+\mathrm{h\Δt})-P_R^f)\end{array}$

，其中R∈I(t)表示船舶代码且I(t)＝{1，2，...，Z_t}，P_R(t+hΔt)表示船舶在时刻(t+hΔt)的位置向量，表示船舶R的解脱终止点，u_R表示待优化的船舶R的最优控制序列，Q_Rt为正定对角矩阵，其对角元素为船舶R在t时刻的优先级指数L_Rt，并且 $Q_{\mathrm{Rt}}=\left[\begin{array}{cc}{L}_{\mathrm{Rt}}& 0\\ 0& L_{\mathrm{Rt}}\end{array}\right].$

9.根据权利要求1至9之一所述的航海交通管制方法，其特征在于：所述步骤⑤中终止参考点位置P设定为船舶运行的下一个航道点，避撞策略控制时域Θ为300秒；轨迹预测时域W为300秒。