AIS数据修复-三次样条插值法（Cubic spline interpolation）

"前面我们讨论了关于船舶AIS的压缩(这里主要指船舶的轨迹)，压缩方法其中包括：基于时间比率的算法、基于时间-速度-航向的算法和基于改进的DP算法。"

在本次博客中，我们对船舶轨迹的修复方法做一个总结。由于船舶在航行的过程中会受到各种异常问题的影响，如 “AIS的位置传感器导致的异常、AIS信号传输过程导致的轨迹异常和AIS网络通信阻塞导致的轨迹异常等” ，进而导致船舶轨迹点不连续。现有的船舶轨迹修复方法有，基于几何的方法，基于人工智能的方法等。本此我们主要研究基于几何的方法——三次样条插值法。

三次样条插值法涉及到《数值分析》，我们对修复原理不再阐述，可参考《舶舶AIS轨迹异常的自动检测与修复算法》一文。下面针对船舶轨迹的部分段缺失，编写了修复算法——三次样条插值法。

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

注意 “pycubicspline” 是自定义的.py文件，这个文件参考了Github上一篇大牛的博客。

# -*- 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）和参考论文，在我的博客“资源”界面下，需要下载~

欢迎留言，欢迎指正~