在上两篇博客中,对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 = 0
j = 1
while i < len(lon):
if mmsi[i] == Temporary_variable_miss:
j = j + 1
Temporary_variable_miss = mmsi[i]
else:
AIS_MMIS_NUM_LIST.append(j)
AIS_MMIS.append(mmsi[i])
j = 1
Temporary_variable_miss = mmsi[i]
i = i + 1
AIS_MMIS_NUM_LIST.append(j)
AIS_MMIS_NUM_LIST[0] = AIS_MMIS_NUM_LIST[0] - 1
'''------不同MMSI号对应着不同的AIS数据个数,创建这个列表------'''
i = 0
j = 0
while i < AIS_MMIS_NUM:
for j in range(AIS_MMIS_NUM_LIST[i]):
compression_data[i].append([])
i = i + 1
'''------往compression_data赋值------'''
i = 0
j = 0
k = 0
for 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 + 1
i = 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 9999999
slope = (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 distance
class 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, 0
for index, point in enumerate(point_list):
if index in [0, len(point_list) - 1]:
continue
distance = point2LineDistance(point, point_list[0], point_list[-1])
if distance > max_distance:
max_distance_index = index
max_distance = distance
# ---------------2.找到速度差值均值的最大值---------------#
previous_velocity_difference = 0
latter_velocity_difference = 0
ave_velocity_difference = 0
for index, point in enumerate(point_list):
if index in [0, len(point_list) - 1]:
continue
previous_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))/2
max_speed_index = index
# ---------------3.找到航向的最大值---------------#
heading_max = 0
for index, point in enumerate(point_list):
if index in [0, len(point_list) - 1]:
continue
if 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)和参考论文,在我的博客“资源”界面下,需要下载~
欢迎留言,欢迎指正~