2. 航迹几何特征提取与分段

为什么要在清洗后立即做几何特征提取

例1输出的干净航迹只是一串按时间排列的经纬度点，它本身不携带“船在往哪走、走多快、是否在转弯”的信息。后续的分段优化、约束建模都依赖这些几何量，因此本例的任务是：对清洗后的点序列逐点计算航向 COG、航速 SOG、曲率 κ 与转角率 ω，再依据转角率与航向变化识别转向点，把整条航迹切成若干单调、近似直线的航段。

四个特征的定义与差分方式必须统一，否则分段结果会随实现漂移。设第 i 个点为 P_i=(λ_i, φ_i, t_i)，经度 λ、纬度 φ 单位为度，t 为秒：

• 航向：θ_i = atan2(Δλ·cos φ̄, Δφ)，其中 φ̄ 为两点平均纬度，结果归一到 [0°, 360°)。
• 航速：v_i = d_i / Δt_i，d_i 用 haversine 公式计算的大圆距离（米），v 单位 m/s，展示时换算为节（1 kn = 0.5144 m/s）。
• 转角率：ω_i = Δθ_i / Δt_i，Δθ_i 必须做角度环绕处理（把差值折算到 (−180°, 180°]），否则 359°→1° 会被误判为 −358° 的大转向。
• 曲率：κ_i = |Δθ_i| / d_i，单位 rad/m，反映单位航程上的方向改变量，与航速无关，适合做几何分段判据。

特征计算与转向点识别的实现

下面是一段可直接运行的 Python 实现，输入为例1输出的干净航迹（CSV，含 lon、lat、timestamp 三列，timestamp 为 Unix 秒），输出带特征的 DataFrame 与分段结果。示例参数取 AIS 常见上报间隔 10 s、转向判定阈值 ω_th=0.5 °/s、航向累计变化阈值 Δθ_th=15°，这些是教学示例值，实际应按船型与采样率调整。

import numpy as np
import pandas as pd

R_EARTH = 6371000.0  # 米

def haversine(lon1, lat1, lon2, lat2):
    lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])
    dlon, dlat = lon2 - lon1, lat2 - lat1
    a = np.sin(dlat/2)**2 + np.cos(lat1)*np.cos(lat2)*np.sin(dlon/2)**2
    return 2 * R_EARTH * np.arcsin(np.sqrt(a))

def wrap180(deg):
    """把角度差折算到 (-180, 180]"""
    return (deg + 180.0) % 360.0 - 180.0

def extract_features(df, omega_th=0.5, dtheta_th=15.0):
    df = df.sort_values("timestamp").reset_index(drop=True)
    lon, lat, t = df["lon"].values, df["lat"].values, df["timestamp"].values

    # 相邻点差分
    dlon = np.diff(lon); dlat = np.diff(lat); dt = np.diff(t)
    lat_mid = np.radians((lat[:-1] + lat[1:]) / 2.0)

    # 航向（度），归一到 [0,360)
    cog = np.degrees(np.arctan2(dlon * np.cos(lat_mid), dlat)) % 360.0
    # 大圆距离与航速
    dist = haversine(lon[:-1], lat[:-1], lon[1:], lat[1:])
    sog = dist / np.maximum(dt, 1e-6)          # m/s
    # 转角率与曲率
    dtheta = wrap180(np.diff(cog))              # 相邻航向差，度
    omega = dtheta / np.maximum(dt[1:], 1e-6)   # °/s
    kappa = np.abs(np.radians(dtheta)) / np.maximum(dist[1:], 1e-6)  # rad/m

    # 对齐到点序列（首点无前向差分，置 NaN）
    df["cog"] = np.r_[np.nan, cog]
    df["sog_kn"] = np.r_[np.nan, sog / 0.5144]
    df["omega"] = np.r_[np.nan, np.nan, omega]
    df["kappa"] = np.r_[np.nan, np.nan, kappa]

    # 转向点：转角率超阈值，或窗口内航向累计变化超阈值
    turn = np.zeros(len(df), dtype=bool)
    turn[2:] = np.abs(omega) > omega_th
    # 合并密集转向点，保留窗口内航向变化最大的点
    idx = np.where(turn)[0]
    if len(idx) > 1:
        keep = [idx[0]]
        for j in idx[1:]:
            if j - keep[-1] <= 3:               # 3 点窗口内视为同一次转向
                if abs(omega[j-2]) > abs(omega[keep[-1]-2]):
                    keep[-1] = j
            else:
                keep.append(j)
        turn[:] = False
        turn[keep] = True

    # 分段：转向点作为航段边界
    seg_id = np.cumsum(turn)
    df["is_turn"] = turn
    df["seg_id"] = seg_id
    return df

# 使用示例
# df = pd.read_csv("clean_track.csv")
# out = extract_features(df)
# print(out[["timestamp","cog","sog_kn","omega","kappa","is_turn","seg_id"]].head(20))
# print("航段数:", out["seg_id"].nunique())

关键点有三处：其一，wrap180 必须在求转角率之前调用，否则跨 0°/360° 的航向会被算成巨大跳变；其二，np.diff 使特征序列比点序列短 1~2 个元素，用 np.r_ 前置 NaN 对齐，保证每行点都能查到自己的特征；其三，转向点合并窗口（示例取 3 点）用于抑制 AIS 抖动造成的连续误判，窗口内只保留转角率最大的点作为该次转向的代表。

特征与分段的判读关系

下图把一条航迹的四个特征沿航程展开，并标出转向点如何把航迹切成航段。可以看到：直航段内 κ 接近 0、ω 在 0 附近小幅波动；进入转向时 ω 与 κ 同步抬升，越过阈值的位置即为分段边界；航速在转向段通常略有下降，可作为辅助佐证但不能单独作为分段依据。

分段结果的验证与输出

分段是否合理，不能只看航段数量，要回到几何量上核对。建议对每个航段统计三项指标：航段内航向标准差（应小于约 5°，说明该段近似直线）、航段长度（过短的航段往往是抖动误判，可合并到相邻段）、以及航段首尾航向差（应接近该次转向的实际转角）。下面给出一个可直接套用的航段汇总与校验片段：

def summarize_segments(df):
    rows = []
    for sid, g in df.groupby("seg_id"):
        g = g.dropna(subset=["cog"])
        if len(g) < 2:
            continue
        # 航向标准差（环绕处理后再算）
        c = np.radians(g["cog"].values)
        mean_c = np.arctan2(np.sin(c).mean(), np.cos(c).mean())
        std_deg = np.degrees(np.sqrt(-2*np.log(np.hypot(np.sin(c).mean(), np.cos(c).mean()))))
        length = haversine(g["lon"].iloc[0], g["lat"].iloc[0],
                           g["lon"].iloc[-1], g["lat"].iloc[-1])
        rows.append({
            "seg_id": sid,
            "点数": len(g),
            "航段长度_m": round(length, 1),
            "航向标准差_deg": round(std_deg, 2),
            "首航向": round(g["cog"].iloc[0], 1),
            "末航向": round(g["cog"].iloc[-1], 1),
            "平均航速_kn": round(g["sog_kn"].mean(), 2),
        })
    return pd.DataFrame(rows)

# seg = summarize_segments(out)
# print(seg)
# 校验：航向标准差 > 5° 的航段需复查是否漏判转向
# 校验：航段长度 < 200 m 的航段考虑与相邻段合并

通过条件建议设为：所有航段航向标准差 ≤ 5°、无长度小于 200 m 的碎片航段、航段数在合理范围（示例航迹 4 段）。若某航段航向标准差明显偏大，说明该段内还藏着未被识别的转向，应下调 ω 阈值或缩小合并窗口后重跑；若碎片航段过多，则说明阈值过敏感，应上调 ω 阈值或扩大合并窗口。最终把带 seg_id 与四个特征的 DataFrame 落盘为 CSV，即为例3约束建模与后续分段优化的输入。