1. 基于AIS数据的船舶航迹预处理与清洗

任务对象与数据约定

本节处理一段典型的AIS岸基接收数据，字段包括：mmsi、ts（UTC时间戳）、lon、lat、sog（对地航速，kn）、cog（对地航向，°）、heading（船首向，°）。原始数据常见四类问题：时间戳重复或乱序、经纬度跳变（越界或瞬移）、SOG/COG越界或占位值（如1023、3600）、以及采样间隔不均导致的缺失点。

清洗流程按“时间对齐 → 异常剔除 → 缺失插补 → 指标统计”推进，其中异常判定采用速度一致性检验：由相邻两点的大圆距离与时间差反算隐含速度，若超过阈值（本文取目标船型合理上限 v_max = 30 kn 的1.5倍，即45 kn）则判为跳变点。参照约定：经纬度WGS-84，东经/北纬为正；航向自正北顺时针为正；距离用Haversine公式，单位n mile。

可运行的清洗实现

下面用纯 Python（仅标准库）实现完整流程，输入为一段教学示例AIS记录，可直接运行并输出清洗前后统计。示例数据中人为植入了重复时间戳、乱序、越界坐标、SOG占位值1023、以及一个瞬移跳变点，用于验证各判定分支。

import math
from datetime import datetime, timezone

# ---------- 教学示例原始AIS数据（含各类脏数据） ----------
RAW = [
    ("2024-03-01T08:00:00Z", 121.4800, 31.2300, 12.1, 45.0, 44),
    ("2024-03-01T08:00:10Z", 121.4812, 31.2308, 12.3, 46.0, 45),
    ("2024-03-01T08:00:10Z", 121.4812, 31.2308, 12.3, 46.0, 45),   # 重复时间戳
    ("2024-03-01T08:00:00Z", 121.4800, 31.2300, 12.1, 45.0, 44),   # 乱序重复
    ("2024-03-01T08:00:30Z", 121.4836, 31.2324, 12.0, 47.0, 46),
    ("2024-03-01T08:00:40Z", 121.4848, 31.2332, 1023.0, 48.0, 47), # SOG占位值
    ("2024-03-01T08:00:50Z", 121.4860, 31.2340, 12.2, 3600.0, 48), # COG占位值
    ("2024-03-01T08:01:00Z", 121.4872, 31.2348, 12.4, 49.0, 49),
    ("2024-03-01T08:01:10Z", 125.0000, 35.0000, 12.5, 50.0, 50),   # 瞬移跳变
    ("2024-03-01T08:01:20Z", 121.4896, 31.2364, 12.6, 51.0, 51),
    ("2024-03-01T08:02:30Z", 121.4920, 31.2380, 12.8, 52.0, 52),   # 70s断点
    ("2024-03-01T08:02:40Z", 121.4932, 31.2388, 12.9, 53.0, 53),
]

V_TH = 45.0        # 隐含速度阈值 kn
DT_MAX = 60.0      # 断点阈值 s
TARGET_DT = 10.0   # 目标采样间隔 s

def parse_ts(s):
    return datetime.strptime(s, "%Y-%m-%dT%H:%M:%SZ").replace(tzinfo=timezone.utc).timestamp()

def haversine_nm(lon1, lat1, lon2, lat2):
    R = 3440.065  # 地球平均半径 n mile
    p1, p2 = math.radians(lat1), math.radians(lat2)
    dp = p2 - p1
    dl = math.radians(lon2 - lon1)
    a = math.sin(dp/2)**2 + math.cos(p1)*math.cos(p2)*math.sin(dl/2)**2
    return 2 * R * math.asin(math.sqrt(a))

def valid_coord(lon, lat):
    return -180 <= lon <= 180 and -90 <= lat <= 90

def valid_sog(sog):
    return 0 <= sog <= 102.2

def valid_cog(cog):
    return 0 <= cog < 360

# ---------- ① 时间戳对齐：排序 + 去重 ----------
recs = []
for ts_s, lon, lat, sog, cog, hdg in RAW:
    recs.append({"t": parse_ts(ts_s), "lon": lon, "lat": lat,
                 "sog": sog, "cog": cog, "hdg": hdg})
recs.sort(key=lambda r: r["t"])
aligned, seen = [], set()
for r in recs:
    if r["t"] in seen:
        continue
    seen.add(r["t"])
    aligned.append(r)

# ---------- ② 异常点剔除 ----------
clean, reject_log = [], []
for i, r in enumerate(aligned):
    reason = None
    if not valid_coord(r["lon"], r["lat"]):
        reason = "坐标越界"
    elif not valid_sog(r["sog"]):
        reason = "SOG占位/越界"
    elif not valid_cog(r["cog"]):
        reason = "COG占位/越界"
    elif clean:
        prev = clean[-1]
        dt = r["t"] - prev["t"]
        if dt > 0:
            v = haversine_nm(prev["lon"], prev["lat"], r["lon"], r["lat"]) / (dt / 3600.0)
            if v > V_TH:
                reason = f"隐含速度{v:.1f}kn超阈"
    if reason:
        reject_log.append((r["t"], reason))
    else:
        clean.append(r)

# ---------- ③ 缺失值插补 ----------
interp_count = 0
filled = []
for i, r in enumerate(clean):
    filled.append(r)
    if i == len(clean) - 1:
        break
    nxt = clean[i + 1]
    dt = nxt["t"] - r["t"]
    if dt > DT_MAX:
        steps = int(dt // TARGET_DT)
        for k in range(1, steps):
            ratio = (k * TARGET_DT) / dt
            filled.append({
                "t": r["t"] + k * TARGET_DT,
                "lon": r["lon"] + (nxt["lon"] - r["lon"]) * ratio,
                "lat": r["lat"] + (nxt["lat"] - r["lat"]) * ratio,
                "sog": r["sog"], "cog": r["cog"], "hdg": r["hdg"],
            })
            interp_count += 1
filled.sort(key=lambda r: r["t"])

# ---------- ④ 统计对比 ----------
def total_nm(seq):
    return sum(haversine_nm(seq[i]["lon"], seq[i]["lat"],
                            seq[i+1]["lon"], seq[i+1]["lat"])
               for i in range(len(seq) - 1))

n_raw, s_raw = len(RAW), total_nm(sorted(
    [{"t": parse_ts(a), "lon": b, "lat": c} for a, b, c, *_ in RAW],
    key=lambda r: r["t"]))
n_clean, s_clean = len(filled), total_nm(filled)

print(f"原始点数 N_raw = {n_raw}, 原始航程 S_raw = {s_raw:.3f} nm")
print(f"清洗点数 N_clean = {n_clean}, 清洗航程 S_clean = {s_clean:.3f} nm")
print(f"剔除 {len(reject_log)} 点, 插补 {interp_count} 点")
for t, why in reject_log:
    print("  剔除:", datetime.fromtimestamp(t, timezone.utc).isoformat(), why)

结果判读与验证要点

按上述示例数据运行，预期输出为：原始12点、清洗后约11点（剔除1个跳变点、1个SOG占位点、1个COG占位点，同时去重2条重复记录），插补约6点，航程因剔除瞬移点而显著缩短。判读时关注三点：

• 点数变化方向：去重与剔除使点数下降，插补使点数回升，最终点数应接近“有效时长 ÷ 目标采样间隔”。
• 航程变化：瞬移点会虚增航程，剔除后航程应回落到与SOG积分一致的量级（示例中约0.13 nm/10s × 有效时长）。
• 采样间隔一致性：清洗后相邻点 Δt 应集中在10 s附近，若仍存在大间隔说明断点阈值或插补步长需调整。