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

Shize Zhang

已于 2022-09-28 10:46:03 修改

阅读量1.8k

点赞数 11

分类专栏： python AIS 船舶自动识别系统文章标签：算法机器学习人工智能

于 2022-09-27 23:04:04 首次发布

本文链接：https://blog.csdn.net/john_zsz/article/details/127074531

版权

python 同时被 3 个专栏收录

5 篇文章 2 订阅

订阅专栏

AIS

4 篇文章 1 订阅

订阅专栏

船舶自动识别系统

4 篇文章 2 订阅

订阅专栏

在上两篇博客中，对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)

三、实验结果