CN108334897A - 一种基于自适应高斯混合模型的海上漂浮物轨迹预测方法 - Google Patents

Abstract

本发明涉及机器学习领域，提出了一种海洋轨迹聚类及预测的方法。为了准确预测未来的轨迹点需要先进行轨迹聚类。本发明轨迹聚类方法主要是将海上复杂多变性和强波动性的轨迹点进行相似性度量，挖掘其中潜在的数据信息。该方法将高斯混合模型GP与狄利克雷过程DP相结合，利用DP的非参贝叶斯框架，确定聚类数目提高聚类自适应性。算法采用在DP基础上加入中国餐馆过程，采用坍塌的Gibbs sampling方法进行模型求解，实现了由有限混合模型到无限混合模型的无监督分类，可自动获取聚类数目,对已聚类的轨迹通过高斯过程回归预测方法预测未来的轨迹点。本发明避免了人工指定聚类数目和参数估计时局部最大化的弊端，在保证自适应聚类的前提下，提高预测的准确率。

Description

一种基于自适应高斯混合模型的海上漂浮物轨迹预测方法

技术领域：

本发明涉及机器学习与海洋轨迹领域，具体地说，是一种通过浮标实时获取漂浮物的轨迹点进行聚类及预测的方法。

背景技术：

随着海洋运输业的快速发展，通过轨迹分析是对海上移动对象最常用的方法，而海上复杂多变和波动性强的轨迹给数据挖掘工作带来了很大的挑战。由于在对轨迹数据进行分析之前并不知道彼此之间的信息，所以轨迹聚类方法非常适合对移动对象轨迹数据进行挖掘。聚类就是将相似运动模式的轨迹，通过分析轨迹自身所具有的结构特征进而确定轨迹间的相似程度，而后将具有相似度较高的轨迹归为一类。传统的聚类算法都需要事先确定聚类数目，如k均值聚类算法、高斯混合模型等。考虑到海上移动对象的高度动态化的轨迹特性，非参贝叶斯方法就适于处理复杂环境的数据。而传统高斯混合模型具有人工指定聚类数目的随意性和参数估计产生局部最大化的弊端，所以我们建议将移动模式建模为高斯过程即GP与狄利克雷过程即DP的结合。GP为每个移动模式提供灵活的表示，而DP则将观察到的轨迹分配给特定的移动模式。GP和DP与已存在的模型区别是两者可根据可用的数据动态调整模型的复杂性，且以原则性和强大的框架来处理不同长度的轨迹，而不仅仅是分离一个数据点。我们将此模型应用到海上移动对象的轨迹聚类研究中，在有噪音情景下分析轨迹的运动模式，使轨迹平滑化进而获取自适应轨迹聚类数目。最后对聚类的数据利用高斯过程回归方法进行轨迹预测，以对海上搜救和异常点检测等方面提供重要研究方向。

发明内容：

为了解决人工指定聚类数目的随意性，本发明的目的是弥补传统高斯混合模型所面临的不足，提出了一种基于自适应高斯混合模型的海上漂浮物轨迹预测方法。

本发明采取的技术方案是：将高斯混合模型与狄利克雷过程即DP相结合，利用DP的非参贝叶斯框架，确定聚类数目提高聚类自适应性。算法采用在高斯混合模型基础上加入中国餐馆过程思想，采用坍塌的Gibbs sampling方法进行模型近似求解，实现了由有限混合模型到无限混合模型的无监督分类，可自动获取聚类数目。在保证自适应聚类的前提下提高预测准确率，该方法包括步骤1.1和步骤1.2：

步骤1.1：对海上漂浮物的轨迹将高斯混合模型和狄利克雷混合模型相结合进行自适应聚类；

步骤1.2：对聚类好的轨迹利用高斯过程回归方法进行预测；

步骤1.1包括以下步骤：

步骤11：建立模型

通过高斯混合模型建立所需模型，采用经度x和纬度y的数据集,由于聚类数目k未知，(x_i,y_i)属于无监督学习，参数分布未知，则需要一个分布的分布即狄利克雷过程(DP)，其最重要的作用是作为数据所属分布的先验存在。其中(x₁,x₂,...x_N；y₁,y₂,...y_N)～Dir(α₁,α₂,...α_k)，Dir表示狄利克雷分布，是一种边缘分布，将其推广至无限维就是DP，则使得θ～DP(α,H)，α是集中度调节参数，H是基分布。

步骤12：估计参数

设浮漂物的数据集为X＝{x₁,x₂,...x_n；y₁,y₂,...y_n}，其中它的两个参数为Z＝{z₁,z₂,...z_n},n∈{1,2,...K},θ＝{θ₁,θ₂,...θ_k}其中Z为隐变量，表示为聚类样本的标签，z_i＝k代表当前第i个类有k个对应成员，θ是该模型的各个类的成员参数，根据贝叶斯理论可知：p(θ,z|X)∝p_X(θ)p_X(z)p(X|θ,z)，则参数θ的后验分布可以通过计算其似然函数及先验分布来得以实现，在此基础上计算出θ的后验分布，并通过Gibbs Sampling的方法来不断更新参数θ。

为了更简便的求解似然函数p(X|θ,z)，因现有数据均服从高斯分布特点，θ即为期望值μ，则p(x_i|θ,z)～N(μ_zi,∑)，对于每一类的先验p₀(z)，则采用CRP模型中的方法产生：

其中，将CRP模型中的每张桌子类比成一个类，其中z_-i初了i之外其他的成员，n_k代表当前在第k个类的轨迹点，α表示狄利克雷的参数集中度。

步骤13：求解模型

Gibbs sampling为各个值分配类标签z_i,然后更新各部件的分布参数θ_k和类别数目K。对于每一个z_i,具体的采样方法如下：

选择坐在已有类即第k个类的概率：

选择新开一个类即第k+1个类的概率：

当轨迹点增加时，此类的参数θ_k就要更新，其方法如下：

步骤1.2包括以下步骤：

步骤21：假设训练集为D＝(x,y),其中，输入经度值为x,输出纬度值为y，输入测试数据D＝(x^*,y^*),通过对较近的两个类簇进行高斯过程回归，则y*关于y的回归函数，即y*的预测值其公式如下：

对应的方差为：

其中

步骤22：从最近的历史轨迹中获取未来d+1时刻的X和Y方向上的预测增量Δx_d+1,Δy_d+1,即能得到该位置点的预测值：

s_d+1＝(x_d+1,y_d+1,t_d+1)＝((x_d+Δx_d+1),y_d+Δ_yd+1,t_d+1)

另外，在已知历史轨迹d+1时刻位置的真实值s_d+1＝(x_d+1,y_d+1,t_d+1)的基础上，可定量求得Δs_d+1为：

其中Δs_d+1为d+1时刻位置的预测误差。

该方法充分利用相邻轨迹之间没有意义相近的说法，根据不同特征轨迹间的相似度对轨迹进行自适应聚类和预测，合理的弥补了高斯混合模型的人工指定聚类数目的随意性和产生局部最大化的不足。

附图说明：

图1狄利克雷过程混合模型的图模型表示

图2CRP模型构造过程表示

图3自适应轨迹聚类及预测流程图

图4误差平方和对比图

图5预测准确率对比图

具体实施方式：

为了使本发明实现的技术手段、创作特性、达成目的与功效易于明白了解，下面结合一种实施方案中NOAA的实时获取的x方向经度值和y方向纬度值的浮标轨迹集，进一步阐述本发明，包括步骤1.1和步骤1.2，其流程图如图3所示：步骤1.1：对海上漂浮物的轨迹将高斯混合模型和狄利克雷混合模型相结合进行自适应聚类；

步骤1.2：对聚类好的轨迹利用高斯过程回归方法进行预测；

步骤1.1包括以下步骤：

步骤11：建立模型

通过高斯混合模型建立所需模型，采用经度值x方向和纬度值y方向,由于聚类数目k未知，(x_i,y_i)属于无监督学习，参数分布未知，则需要一个分布的分布即狄利克雷过程，表示为(x₁,x₂,...x_N；y₁,y₂,...y_N)～Dir(α₁,α₂,...α_k)，Dir表示狄利克雷分布，是一种边缘分布，将其推广至无限维就是DP，则使得θ～DP(α,H)，α是集中度调节参数，H是基分布。

步骤12：估计参数

设浮漂物的数据集为X＝{x₁,x₂,...x_n；y₁,y₂,...y_n}，其中它的两个参数为Z＝{z₁,z₂,...z_n},n∈{1,2,...K},θ＝{θ₁,θ₂,...θ_k}其中Z为隐变量，表示为聚类样本的标签，z_i＝k代表当前第i个类有k个对应成员，θ是该模型的各个类的成员参数，根据贝叶斯理论可知：p(θ,z|X)∝p_X(θ)p_X(z)p(X|θ,z)，则参数θ的后验分布可以通过计算其似然函数及先验分布来得以实现，在此基础上计算出θ的后验分布，并通过Gibbs Sampling的方法来不断更新参数θ。

为了求解似然函数p(X|θ,z)，因现有数据均服从高斯分布特点，θ即为期望值μ，则p(x_i|θ,z)～N(μ_zi,Σ)，对于每一类的先验p₀(z)，则采用CRP模型中的方法产生：

其中，将CRP模型中的每张桌子类比成一个类，其中z_-i初了i之外其他的成员，n_k代表当前坐在第k个类的其他轨迹点，α表示狄利克雷的参数集中度。。

步骤13：求解模型

Gibbs sampling为各个值分配类标签z_i,然后更新各部件的分布参数θ_k和类别数目K。对于每一个z_i,具体的采样方法如下：

选择坐在已有类即第k个类的概率：

选择新开一个类即第k+1个类的概率：

当轨迹点增加时，此类的参数θ_k就要更新，其方法如下：

步骤1.2包括以下步骤：

步骤21：假设训练集为D＝(x,y),其中，输入经度值为x,输出纬度值为y，输入测试数据D＝(x^*,y^*),通过对较近的两个类簇进行高斯过程回归，则y*关于y的回归函数，即y*的预测值其公式如下：

对应的方差为：

其中

步骤22：从最近的历史轨迹中获取未来d+1时刻的X和Y方向上的预测增量Δx_d+1,Δy_d+1,即能得到该位置点的预测值：

s_d+1＝(x_d+1,y_d+1,t_d+1)＝((x_d+Δx_d+1),y_d+Δ_yd+1,t_d+1)

另外，在已知历史轨迹d+1时刻位置的真实值s_d+1＝(x_d+1,y_d+1,t_d+1)的基础上，可定量求得Δs_d+1为：

其中Δs_d+1为d+1时刻位置的预测误差。

控制对象：

如图1为狄利克雷混合过程(DPMM)构造模型：

θ_i|G～G

x_i|θ_i～F(θ_i)

G|α₀，G₀:DP(α₀,G₀)

其中，θ_i服从参数为F(θ_i)的分布，参数θ_i服从概率测度G,而G可通过狄利克雷过程构造。θ_i既可以是单个参数，也可以是多个参数构成的向量，参数α是一个正实数。

如图2所示为中国餐馆模型(CRP)的构建模型：

其中，i表示餐馆里的桌子；n_i为第i张桌子上的顾客数；α为狄利克雷过程的先验参数，是一个正实数；Z为隐变量，表示为聚类样本的标签(label)，z_i＝k代表当前第i个类有k个对应成员；Z_-i为除了i之外的所有成员。假设一个中国餐馆有无限的桌子，第一个顾客到来之后坐在第一张桌子上。第二个顾客来到可以选择坐在第一张桌子上，也可以选择坐在一张新的桌子上，假设第n+1个顾客到来的时候，已经有k张桌子上有顾客了，分别坐了n₁,n₂,...n_k个顾客，那么第n+1个顾客可以以概率为坐在第i张桌子上，n_i为第i张桌子上的顾客数；同时有概率为选取一张新的桌子坐下。那么在n个顾客坐定之后，很显然CRP把这n个顾客分为了K个堆，即K个clusters。

仿真实验：

为了验证本发明方法的有效性，将本发明控制方法和现有的大数据环境下移动对象自适应轨迹预测模型、传统高斯混合模型进行聚类性能对比。

仿真的对象为NOAA的The GDPDrift Data Assembly Center(DAC)的HourlyData，其中每一天包含200个浮标的24个小时实时轨迹，连续记录180天，由于轨迹数量较大，本文就从200个浮标任意抽取其中的5个，然后抽取其中三个月的轨迹每隔72小时进行切分，选取波动性较大，复杂多变的轨迹组成了数据集{D1,D2,D3,D4,D5}。在仿真实验中，海上移动对象的轨迹收到天气和自然灾害的影响呈现为不规则，振动性强的轨迹，这是仿真实验的挑战。

为了更好地评价算法的性能优劣和衡量聚类准确率，其选取的是k＝4,α＝1.0条件下，将本文方法(记为A-GMM方法)与大数据环境下移动对象自适应轨迹预测模型(简记为HMM方法)、传统高斯混合模型(GMM)的误差平方和SSE、纯度和F1值分别见图4、图5、表1和表2。

表1GMM方法、HMM方法和A-GMM方法纯度值对比

表2GMM方法、HMM方法和A-GMM方法的F1值对比

仿真结果分析：

(1)通过观察表1可知，A-GMM纯度值平均高于GMM约为10％左右，高于HMM约为9％，说明了A-GMM不会因轨迹的复杂多样性而受到太大影响，聚类效果较佳。而HMM和GMM对特征变化较明显或者波动性强的轨迹聚类效果不理想；

(2)F1值是查全率和准确率的加权平均，其值越高，说明聚类性能较好。通过观表2可知，GMM和HMM的F1值相近，但均低于A-GMM约6％和5％左右，说明A-GMM的聚类可信度较高。

(3)通过观察图4的纵轴的误差平方和可知，GMM和HMM的误差平方和均高于A-GMM，而SSE越高说明类间相似性不高、差别较大，致使前两者聚类效果大大降低，而A-GMM每组的SSE低于前两者的数值约为10左右，很大程度上减少了误差，聚类可靠性较高；

(4)通过观察图5，在同一条件下对不同数量轨迹进行预测，A-GMM的预测准确率明显高于HMM和GMM，说明此发明对海上浮漂物的轨迹预测具有较好的效果，可以规避不测时间的发生。

综上所述，本发明所提方法能在自适应聚类的同时，有效提高了聚类准确率，在实际航途中可以规避不测时间的发生，减少人员和财产损失具有实际意义。

Claims (1)

1.一种基于自适应高斯混合模型的海上漂浮物轨迹预测方法，其特征在于包括步骤1.1、步骤1.2：

步骤1.1：对海上漂浮物的轨迹将高斯混合模型和狄利克雷混合模型相结合进行自适应聚类；

步骤1.2：对聚类好的轨迹利用高斯过程回归方法进行预测；

步骤1.1包括以下步骤：

步骤11：建立模型

通过高斯混合模型来建立所需模型，给定漂浮物轨迹点(x₁,x₂,...x_n；y₁,y₂,...y_n)，其中x表示经度，y表示纬度，已知聚类数目k，每个类簇都有一个对应均值μ和方差Σ，给出一个初始分布，则第i个数据(x_i,y_i)的似然函数为p(x_i,y_i|θ)，其中θ是关于(μ，Σ)的分布，但是现在聚类数目k未知，(x_i,y_i)属于无监督学习，参数分布未知，所以无法使用EM进行参数估计，则需要一个分布的分布即狄利克雷过程，其最重要的作用是作为数据所属分布的先验存在，其中(x₁,x₂,...x_N；y₁,y₂,...y_N)～Dir(α₁,α₂,...α_k)，Dir表示狄利克雷分布，是一种边缘分布，将其推广至无限维就是狄利克雷混合模型DPMM，则使得θ～DP(α,H)，α是集中度调节参数，H是基分布；

步骤12：估计参数

浮漂物的数据集为X＝{x₁,x₂,...x_n；y₁,y₂,...y_n}，它的两个参数为Z＝{z₁,z₂,...z_n},n∈{1,2,...K},θ＝{θ₁,θ₂,...θ_k}其中Z为隐变量，表示为聚类样本的标签，z_i＝k代表当前第i个类有k个对应成员，θ是该模型的各个类的成员参数，根据贝叶斯理论可知：p(θ,z|X)∝p_X(θ)p_X(z)p(X|θ,z)，则参数θ的后验分布可以通过计算其似然函数及先验分布来得以实现，在此基础上计算出θ的后验分布，并通过Gibbs Sampling的方法来不断更新参数θ；

为了更简便的求解似然函数p(X|θ,z)，因现有数据均服从高斯分布特点，θ即为期望值μ，则p(x_i|θ,z)～N(μ_zi,Σ)，对于先验p₀(θ)，可通过其共轭先验p₀(θ_i)～N(0,σ²I)，对于每一类的先验p₀(z)，则采用CRP模型中的方法产生：

其中，将CRP模型中的每张桌子类比成一个类，其中z_-i初了i之外其他的成员，n_k代表当前坐在第k个类的其他轨迹点，α表示狄利克雷的参数集中度；

步骤13：求解模型

Gibbs sampling为各个值分配类标签z_i,然后更新各部件的分布参数θ_k和类别数目K，对于每一个z_i,具体的采样方法如下：

选择坐在已有类即第k个类的概率：

选择新开一个类即第k+1个类的概率：

当轨迹点增加时，此类的参数θ_k就要如下更新：

步骤1.2包括以下步骤：

步骤21：假设训练集为D＝(x,y)其中，输入经度值为x,输出纬度值为y，输入测试数据D＝(x^*,y^*),通过对较近的两个类簇进行高斯过程回归，则y*关于y的回归函数，即y*的预测值其公式如下：

对应的方差为：

其中

步骤22：从最近的历史轨迹中获取未来d+1时刻的X和Y方向上的预测增量Δx_d+1,Δy_d+1,即能得到该位置点的预测值：

s_d+1＝(x_d+1,y_d+1,t_d+1)＝((x_d+Δx_d+1),y_d+Δ_yd+1,t_d+1)

另外，在已知历史轨迹d+1时刻位置的真实值s_d+1＝(x_d+1,y_d+1,t_d+1)的基础上，可定量求得Δs_d+1为：

其中Δs_d+1为d+1时刻位置的预测误差。