"前面我们讨论了关于船舶AIS的压缩(这里主要指船舶的轨迹),压缩方法其中包括:基于时间比率的算法、基于时间-速度-航向的算法和基于改进的DP算法。"
在本次博客中,我们对船舶轨迹的修复方法做一个总结。由于船舶在航行的过程中会受到各种异常问题的影响,如 “AIS的位置传感器导致的异常、AIS信号传输过程导致的轨迹异常和AIS网络通信阻塞导致的轨迹异常等” ,进而导致船舶轨迹点不连续。现有的船舶轨迹修复方法有,基于几何的方法,基于人工智能的方法等。本此我们主要研究基于几何的方法——三次样条插值法。
三次样条插值法涉及到《数值分析》,我们对修复原理不再阐述,可参考《舶舶AIS轨迹异常的自动检测与修复算法》一文。下面针对船舶轨迹的部分段缺失,编写了修复算法——三次样条插值法。
1、首先导入所用到的包:
注意 “pycubicspline” 是自定义的.py文件,
# -*- coding: utf-8 -*-
"""
Description : 三次样条插值法
Author : 张尺
date : 2022/07/15
"""
import math
import warnings
import pycubicspline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
"""-----设置属性防止中文乱码及拦截异常信息-----"""
plt.rcParams['font.sans-serif'] = [u'SimHei']
plt.rcParams['axes.unicode_minus'] = False
warnings.filterwarnings('ignore', category=FutureWarning)
2、第二部分主要完成的是,把EXCEL中的AIS数据赋值到列表Interpolation_data中。
"""-----主函数-----"""
fig = plt.figure()
AIS_MMIS = [] # 储存各个船舶的MMSI号
AIS_MMIS_NUM = 15 # MISS号的个数,手动输入
AIS_MMIS_NUM_LIST = [] # MISS号每个船舶,对应的经纬度个数
shipdata = pd.read_excel("测试数据(15条).xlsx")
mmsi = shipdata['MMSI']
Lon = shipdata['Lon']
Lat = shipdata['Lat']
Temporary_variable_miss = mmsi[0]
# print(shipdata.index) # 获取行的索引名称
# print(shipdata.columns) # 获取列的索引名称
AIS_MMIS.append(mmsi[0]) # 先赋值第一个船舶的MMSI号
"""-------创建要插值的数据的AIS_MMIS_NUM个向量-------"""
Interpolation_data = list()
for i in range(AIS_MMIS_NUM):
Interpolation_data.append([])
"""-------统计各个MMSI号,经纬度个数-------"""
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
print(AIS_MMIS_NUM_LIST)
print(AIS_MMIS)
"""-------往各个MMSI号的列表中加入对应经纬度组的个数-------"""
i = 0
j = 0
while i < AIS_MMIS_NUM:
for j in range(AIS_MMIS_NUM_LIST[i]):
Interpolation_data[i].append([])
i = i + 1
"""-------往Interpolation_data赋值-------"""
i = 0
j = 0
k = 0
for i in range(AIS_MMIS_NUM):
for j in range(AIS_MMIS_NUM_LIST[i]):
Interpolation_data[i][j].append(Lon[k])
Interpolation_data[i][j].append(Lat[k])
k = k + 1
i = i + 1
3、△ 第三部分,也是最重要的一部分。
(1)插值船舶轨迹的数据,一是经度 lon,二是纬度 lat。我们需要把 lon 或者 lat 任意一个作为插值坐标系中的 x 轴和 y 轴。考虑到最后的可视化效果,我们以 lon 作为 x 轴、以 lat 作为 y 轴。
(2)我们想得到的是,本算法对任意船舶的轨迹都能进行插值。而在实验过程中,我发现如果碰到一些“特殊”的船舶轨迹,程序就跑不动了。经过了很长一段时间的思考和调试,我发现两个相邻的轨迹点 lon 不能相等,这就意味在 x轴上,两个相邻点的横坐标是不能相等的,这也是因为在三次样条插值法时,涉及到了求一个点的导数,如果两个点的横坐标相等,求得的这点导数是无穷的,就会导致程序终止。因此,我们写了一个算法,把即将要插值轨迹中两个相邻轨迹点 lon 变的不重复。算法如下(为了方便大家测试以 lat 作为x轴,也同时编写了相邻轨迹点 lat 不重复的代码):
"""-------去除重复数据-lon-------"""
i = 0
Interpolation_data_new = []
for i in range(AIS_MMIS_NUM):
Interpolation_data_new.append([])
i = 0
j = 1
while i < AIS_MMIS_NUM:
Interpolation_temp = Interpolation_data[i][0][0]
Interpolation_data_new[i].append(Interpolation_data[i][0])
while j < AIS_MMIS_NUM_LIST[i]:
if Interpolation_data[i][j][0] != Interpolation_temp:
Interpolation_data_new[i].append(Interpolation_data[i][j])
Interpolation_temp = Interpolation_data[i][j][0]
j = j + 1
j = 1
i = i + 1
"""-------去除重复数据-lat-------"""
"""
i = 0
Interpolation_data_new = []
for i in range(AIS_MMIS_NUM):
Interpolation_data_new.append([])
i = 0
j = 1
while i < AIS_MMIS_NUM:
Interpolation_temp = Interpolation_data[i][0][1]
Interpolation_data_new[i].append(Interpolation_data[i][0])
while j < AIS_MMIS_NUM_LIST[i]:
if Interpolation_data[i][j][1] != Interpolation_temp:
Interpolation_data_new[i].append(Interpolation_data[i][j])
Interpolation_temp = Interpolation_data[i][j][1]
j = j + 1
j = 1
i = i + 1
"""
4、开始插值,并画图。
"""-------开始插值,一个循环体-------"""
i = 0
while i < AIS_MMIS_NUM:
Interpolation_data_x = [row[0] for row in Interpolation_data_new[i]] # 原始数据——取第一列数据 lon
Interpolation_data_y = [row[1] for row in Interpolation_data_new[i]] # 原始数据——取第二列数据 lat
print("Spline ship: ", AIS_MMIS[i] ,' ' , i)
x, y, yaw, k, travel = pycubicspline.calc_2d_spline_interpolation\
(Interpolation_data_x, Interpolation_data_y,num=int(len(Interpolation_data_x)*1.3))
""" -------画图------- """
plt.clf() # 用于画多个图时, 清空画面
# -------画图——图1------- #
ax1 = fig.add_subplot(211) # 211,指的就是将这块画布分为2×1,然后1对应的就是1号区
plt.scatter(Interpolation_data_x, Interpolation_data_y, marker='o')
plt.title(u'原始数据分布')
# -------画图——图2------- #
ax2 = fig.add_subplot(212)
plt.scatter(x, y, marker='o')
plt.title('Cubic_spline')
plt.pause(2) # 显示秒数
5、插值结果如下图1-图3所示。
图1 修复轨迹图
图2 修复轨迹图
图3 修复轨迹图
可以看到,利用三次样条插值法对缺失的轨迹段进行修复,得到的船舶完整轨迹是很成功的。
参考文献:吴建华,吴琛,刘文,郭俊纬.舶舶AIS轨迹异常的自动检测与修复算法[J].中国航海,2017,40(01):8-12+101.
完整工程(数据集+Python)和参考论文,在我的博客“资源”界面下,需要下载~
欢迎留言,欢迎指正~