AIS数据压缩-改进的DP算法（Improved DP algorithm）

在上两篇博客中，对AIS数据进行压缩用了两种方法：

1.AIS数据压缩-时间比率算法（Time-Ratio-algorithm）

2.AIS数据压缩-时间比率_速度_航向算法（Time-Speed-Heading-Dased）

《A Comparison of Trajectory Compression Algorithms Over AIS Data》一文中，对比了各类AIS数据压缩算法的性能，如下表所示：

从表中得知：

（1）在压缩比率得分对比中，DP算法的性能是最好的，TR(Time-Ratio-algorithm)算法其次，TSH(Time-Speed-Heading-Dased)算法较差。

（2）在执行时间得分对比中，DP算法也表现出优异的性能，TR(Time-Ratio-algorithm)算法折中，TSH(Time-Speed-Heading-Dased)算法运行时间最长。

（3）在与原始轨迹对比的相似度得分中，DP算法，TR(Time-Ratio-algorithm)，TSH(Time-Speed-Heading-Dased)表现不佳；反而基于航向的HD(Heading-Based)算法，和基于速度的SP(Speeding-Based)算法，在相似得分方面，表现出优异性能。

在以上分析中，可以看到，DP算法除了在原AIS船舶轨迹的相似度方面表现较差外，其他指标方面都很优异。因此，考虑船舶AIS数据的特点，对经典的DP算法进行改进成为必然。在最新的一篇文章中《考虑对地航速和航向的船舶典型轨迹提取方法》中提出了改进DP算法的原理，并给出了算法的执行流程。

一、改进的DP算法，原理和算法流程如下：

“本部分来自《考虑对地航速和航向的船舶典型轨迹提取方法》一文”

DP 算法思路是：假设 P1-P2-P3-P4-P5-P6-P7-P8 为船舶航行轨迹的位置点图，第一步把清洗后的船舶起始点连线，然后计算所有中间航迹点到这条直线的距离，依次比较这些距离，得到最大距离dmax 。第二步把距离最大值dmax 和初始阈值 dth 比较，若 dmax > dth ，那么把中间点全部舍去；若 dmax ≤ dth ，则把这条曲线分成两部分。第三步对这两部分轨迹分别执行上述操作，依次迭代，直到压缩完成 P1-P4-P6-P7-P8，如图 1 所示。

DP 算法是根据点到轨迹段的距离进行划分的，在压缩轨迹的同时保存了轨迹的形状特征。但是当处理移动对象时，传统 DP 算法缺点是没有考虑船舶其他重要维度。因此，本节对传统的 DP 算法进行改进，把航速和航向加入到 DP 算法的压缩过程中，改进的 DP 算法流程图如图 2 所示。

图2 改进的 DP 算法流程图

二、根据以上改机DP算法的原理和程序执行流程图，我们对此程序进行了编写：

1、首先导入所用到的包：

# -*- coding: utf-8 -*-
"""Description :   改进的DP算法Author :        张尺date :          2022/09/23
"""
from   __future__ import division
from   math       import sqrt, pow
import matplotlib.pyplot as plt
import pandas as pd

2、第二部分是主函数代码的编写，这里主要包含：

（1）给出形状阈值 SHAPE_THRESHOLD

（2）从EXCEL表中获取原始AIS数据，保存到列表compression_data中

（3）调用 " 方法：improved_dp_main "

（4）分别画出原始轨迹和采用改进DP算法压缩后的轨迹图

if __name__ == '__main__':SHAPE_THRESHOLD = 0.000009 * 100  # 形状阈值AIS_MMIS = []                     # 储存各个船舶的MMSI号AIS_MMIS_NUM = 1                  # MISS号的个数，手动输入AIS_MMIS_NUM_LIST = []            # MISS号每个船舶，对应的经纬度个数shipdata = pd.read_excel("MMSI-413813194-原始数据.xlsx")mmsi =  shipdata['MMSI']year_month_day     =  shipdata['年月日']hour_minute_second =  shipdata['时分秒']lon = shipdata['LON']lat = shipdata['LAT']sog = shipdata['SOG']cog = shipdata['COG']Temporary_variable_miss = mmsi[0]AIS_MMIS.append(mmsi[0])                      # 先赋值第一个船舶的MMSI号compression_data = list()                     # 创建要压缩的数据的AIS_MMIS_NUM个列表for i in range(AIS_MMIS_NUM):compression_data.append([])'''------统计各个MMSI号的AIS数据个数，这里要求数据必须是对齐的------'''i = 0j = 1while i < len(lon):if mmsi[i] == Temporary_variable_miss:j = j + 1Temporary_variable_miss = mmsi[i]else:AIS_MMIS_NUM_LIST.append(j)AIS_MMIS.append(mmsi[i])j = 1Temporary_variable_miss = mmsi[i]i = i + 1AIS_MMIS_NUM_LIST.append(j)AIS_MMIS_NUM_LIST[0] = AIS_MMIS_NUM_LIST[0] - 1'''------不同MMSI号对应着不同的AIS数据个数，创建这个列表------'''i = 0j = 0while i < AIS_MMIS_NUM:for j in range(AIS_MMIS_NUM_LIST[i]):compression_data[i].append([])i = i + 1'''------往compression_data赋值------'''i = 0j = 0k = 0for i in range(AIS_MMIS_NUM):for j in range(AIS_MMIS_NUM_LIST[i]):compression_data[i][j].append(year_month_day[k])compression_data[i][j].append(hour_minute_second[k])compression_data[i][j].append(lon[k])compression_data[i][j].append(lat[k])compression_data[i][j].append(sog[k])compression_data[i][j].append(cog[k])k = k + 1i = i + 1# print(compression_data)'''------创建压缩之后的数据的列表------'''compression_complete_data = list()for i in range(AIS_MMIS_NUM):compression_complete_data.append([])# print(compression_data[0])'''------开始压缩------'''d = DouglasPeuker()for i in range(AIS_MMIS_NUM):d.improved_dp_main(compression_data[i])compression_complete_data[i] = d.qualify_list# ----变量清零---- #d.shape_threshold = SHAPE_THRESHOLD    # 形状阈值d.speed_threshold = 1                  # 速度阈值d.heading_threshold = 30               # 航向阈值d.qualify_list    = list()d.disqualify_list = list()print('各船舶MMSI号:')print(AIS_MMIS)print('压缩完成各船舶经纬度点:')print(compression_complete_data)'''------画图------'''fig = plt.figure()show_original_line_x = []show_original_line_y = []for show_point in compression_data[0]:show_original_line_x.append(show_point[2])show_original_line_y.append(show_point[3])show_compression_complete_line_x = []show_compression_complete_line_y = []for show_point in compression_complete_data[0]:show_compression_complete_line_x.append(show_point[2])show_compression_complete_line_y.append(show_point[3])ax1 = fig.add_subplot(211)  # 211，指的就是将这块画布分为2×1，然后1对应的就是1号区plt.plot(show_original_line_x, show_original_line_y, color='green', linestyle='--')plt.scatter(show_original_line_x, show_original_line_y, color='green')plt.title('Original_data—413XXXXXX')plt.sca(ax1)ax2 = fig.add_subplot(111)  # 211，指的就是将这块画布分为2×1，然后2对应的就是2号区plt.plot(show_compression_complete_line_x, show_compression_complete_line_y, color='green', linestyle='--')plt.scatter(show_compression_complete_line_x, show_compression_complete_line_y, color='red')plt.title('Compression_data—413XXXXXX, Speed_threshold:0.15and0.2')plt.sca(ax2)plt.show()

3、第三部分是创建：Improved_DP_algorithm 类：

（1）定义“方法 __init__”中的形参speed_threshold 的速度阈值和形参heading_threshold 的转向阈值。

（2）根据 “图2 改进的 DP 算法流程图” 找到 “与首尾两点连线距离的最大值 ”，“找到速度差值均值的最大值“，“找到转向的最大值”。

（3）把上述这些最大值与各自的阈值进行比较，若大于阈值，则分割曲线；若小于，则对下一个属性的最大值进行判断。

def point2LineDistance(point_a, point_b, point_c):"""-------计算点a到点b c所在直线的距离-------"""# 首先计算 b c 所在直线的斜率和截距if point_b[2] == point_c[2]:return 9999999slope = (point_b[3] - point_c[3]) / (point_b[2] - point_c[2])  # 计算斜率intercept = point_b[3] - slope * point_b[2]  # 计算截距# 计算点a到 b c所在直线的距离distance = abs(slope * point_a[2] - point_a[3] + intercept) / sqrt(1 + pow(slope, 2))return distanceclass Improved_DP_algorithm(object):def __init__(self):self.shape_threshold = SHAPE_THRESHOLD          # 形状阈值self.speed_threshold = 0.1                      # 速度阈值self.heading_threshold = 30                     # 航向阈值self.qualify_list = list()self.disqualify_list = list()def method_shape_speed_heading(self, point_list):if len(point_list) < 2:self.qualify_list.extend(point_list[::-1])else:#---------1.找到与首尾两点连线距离最大的点---------#max_distance_index, max_distance = 0, 0for index, point in enumerate(point_list):if index in [0, len(point_list) - 1]:continuedistance = point2LineDistance(point, point_list[0], point_list[-1])if distance > max_distance:max_distance_index = indexmax_distance = distance# ---------------2.找到速度差值均值的最大值---------------#previous_velocity_difference = 0latter_velocity_difference = 0ave_velocity_difference = 0for index, point in enumerate(point_list):if index in [0, len(point_list) - 1]:continueprevious_velocity_difference = point[4] - point_list[index - 1][4]latter_velocity_difference = point_list[index + 1][4] - point[4]if abs(previous_velocity_difference) >= self.speed_threshold and abs(latter_velocity_difference) \>= self.speed_threshold:if (abs(previous_velocity_difference)+abs(latter_velocity_difference))/2 > self.speed_threshold:ave_velocity_difference = (abs(previous_velocity_difference)+abs(latter_velocity_difference))/2max_speed_index = index# ---------------3.找到航向的最大值---------------#heading_max = 0for index, point in enumerate(point_list):if index in [0, len(point_list) - 1]:continueif abs(point_list[index+1][5] - point[5]) > heading_max:heading_max = abs(point_list[index+1][5] - point[5])max_heading_index = index#---------若最大距离小于阈值，则去掉所有中间点。 反之，则将曲线按最大距离点分割---------#if max_distance > self.shape_threshold:# ---------1.将曲线按最大距离的点分割成两段---------#sequence_a = point_list[:max_distance_index]sequence_b = point_list[max_distance_index:]for sequence in [sequence_a, sequence_b]:if len(sequence) < 2 and sequence == sequence_b:self.qualify_list.extend(sequence[::-1])else:self.disqualify_list.append(sequence)elif ave_velocity_difference > self.speed_threshold:# ---------2.将曲线按最大平均速度差值的点分割成两段---------#sequence_a = point_list[:max_speed_index]sequence_b = point_list[max_speed_index:]for sequence in [sequence_a, sequence_b]:if len(sequence) < 2 and sequence == sequence_b:self.qualify_list.extend(sequence[::-1])else:self.disqualify_list.append(sequence)elif heading_max > self.heading_threshold:# ---------3.将曲线按最大航向角的点分割成两段---------#sequence_a = point_list[:max_heading_index]sequence_b = point_list[max_heading_index:]for sequence in [sequence_a, sequence_b]:if len(sequence) < 2 and sequence == sequence_b:self.qualify_list.extend(sequence[::-1])else:self.disqualify_list.append(sequence)else:# ---------都不满足, 去掉所有中间点---------#self.qualify_list.append(point_list[-1])self.qualify_list.append(point_list[0])def improved_dp_main(self, point_list):self.method_shape_speed_heading(point_list)while len(self.disqualify_list) > 0:self.method_shape_speed_heading(self.disqualify_list.pop())# print('压缩之后的数据长度:%s' % len(self.qualify_list))# print('压缩之后的数据:\n %s' %      self.qualify_list)

三、实验结果

原始AIS数据个数：3485个

压缩完成后AIS数据个数: 114个

参考文献：

[1] Makris A , Kontopoulos I , Alimisis P , et al. A comparison of trajectory compression algorithms over AIS data[J]. IEEE Access, 2021, PP(99):1-1.

[2] 刘畅,张仕泽,李倍莹,李波.考虑对地航速和航向的船舶典型轨迹提取方法[J].交通运输系统工程与信息:1-14[2022-09-27].

完整工程（数据集+Python）和参考论文，在我的博客“资源”界面下，需要下载~

欢迎留言，欢迎指正~

AIS数据压缩-改进的DP算法（Improved DP algorithm）相关推荐

模式识别读书报告---关于DP算法的…
我喜欢玩星际争霸,当我操作我的部队用快捷键A进攻敌人时,部队就会自动避开障碍物attack目标.当时我就对为什么部队能自动选择最短的路进攻产生了兴趣,它是用什么算法实现的.而且当你的对手是电脑时,是用 ...
boost::geometry模块多边形DP算法简化示例
boost::geometry模块多边形DP算法简化示例实现功能 C++实现代码实现功能 boost::geometry模块多边形DP算法简化示例 C++实现代码 #include <boo ...
dp笔记：关于DP算法和滚动数组优化的思考
从网上总结了一些dp的套路以及对滚动数组的一些思考,现记录如下,希望以后回顾此类算法时会有所帮助. 目录 1.DP算法经验 1.DP算法核心: 2.DP算法类别以及例题例1:三步问题例2:最小路径 ...
c语言dp算法,C++动态规划dp算法题
问题1:找硬币,换钱的方法输入: penny数组代表所有货币的面值,正数不重复 aim小于等于1000,代表要找的钱输出: 换钱的方法总数解法1:经典dp,空间复杂度O(n*aim) class ...
Leetcode 64. 最小路径和 -- DP算法
Time: 20190831 题目描述给定一个包含非负整数的 m x n 网格,请找出一条从左上角到右下角的路径,使得路径上的数字总和为最小. 说明:每次只能向下或者向右移动一步. 示例: 输入: ...
DP算法：动态规划算法
步骤 (1)确定初始状态 (2)确定转移矩阵,得到每个阶段的状态,由上一阶段推到出来 (3)确定边界条件. 例题蓝桥杯--印章(python实现) 使用dp记录状态,dp[i][j]表示买i张印章, ...
灰度图像压缩 DP算法位运算详解
作者码字不易,白天敲代码,晚上熬夜赶报告,要转载请注明出处哦,程序猿的辛酸泪目录位运算回顾压缩过程解压过程关于一个莫得感情的小bug 实用小工具的下载地址完整版代码位运算回顾若 a = ...
第十四周DP算法总结
这周自己主要再看DP算法的博客,感觉DP这一部分内容确实比之前的都要麻烦一些,最后攻克这一部分难题还是挺好的. 这周自己总结了一些题型,以及一些方法思路,最后再把动态规划和之前的分治和贪心做一下比较. ...
C++入门算法1——dp(动态规划)
dp(动态规划)是十分重要的一个算法,一般来说这种算法会比dfs(深度优先搜索)快很多. 首先先来看一道例题题目链接:P1048 [NOIP2005 普及组] 采药 - 洛谷 | 计算机科学教育新生 ...

AIS数据压缩-改进的DP算法（Improved DP algorithm）

AIS数据压缩-改进的DP算法（Improved DP algorithm）相关推荐

最新文章

热门文章