PostGIS轨迹分析——AIS数据删除异常点

点位连成线

select st_makeline(t.the_geom)
from (select the_geom
      from points
      where name = 'xxx'
        and point_time > '2023-09-01'
        and point_time < '2023-09-03'
      order by point_time ) as t;

异常点

红框区域明显断开了，原因是有一个“飞点”

解决思路

异常点一般距离前一个点的距离会大于某一个阈值（比如50000米），当然这个阈值视情况而定，以下是解决此问题的步骤：

1. 计算距离下一个点的大圆距离
2. 当前点与下一个点的距离，如果距离大于阈值，则认为下一个点是异常点
3. 找出多个异常点，这两个异常点中间的点都可以认为是异常点，都要移除

1.计算距离下一个点的大圆距离

select id,
       the_geom,
       point_time,
       ---计算距离下一个点的大圆距离，单位：米
       st_distance(the_geom::geography, lead(the_geom) over (order by point_time)::geography) as d
from points
where name = 'xxx'
  and point_time > '2023-09-01'
  and point_time < '2023-09-03'
order by point_time;

2.当前点与下一个点的距离，如果距离大于阈值，则认为下一个点是异常点

with agg as (select id,
                    the_geom,
                    point_time,
                    ---计算距离下一个点的大圆距离，单位：米
                    st_distance(the_geom::geography, lead(the_geom) over (order by point_time)::geography) as d
             from points
             where name = 'xxx'
               and point_time> '2023-09-01'
               and point_time< '2023-09-03'
             order by point_time),
     ---超过阈值的无效点
     invalid as (select id,
                        the_geom,
                        point_time,
                        d,
                        row_number() over (order by point_time) as number
                 from agg
                 ---50000为阈值，单位：米
                 where d > 50000
                 order by point_time)
select *
from invalid;

3.找出多个异常点，这两个异常点中间的点都可以认为是异常点，都要移除

with agg as (select id,
                    the_geom,
                    point_time,
                    ---计算距离下一个点的大圆距离，单位：米
                    st_distance(the_geom::geography, lead(the_geom) over (order by point_time)::geography) as d
             from points
             where name = 'xxx'
               and point_time> '2023-09-01'
               and point_time< '2023-09-03'
             order by point_time),
     ---超过阈值的无效点
     invalid as (select id,
                        the_geom,
                        point_time,
                        d,
                        row_number() over (order by point_time) as number
                 from agg
                 ---50000为阈值，单位：米
                 where d > 50000
                 order by point_time),
     ---过滤无效值后的点
     valid as (select id,
                      the_geom,
                      point_time,
                      d
               from agg
               ---条件invalid.number = 1和invalid.number = 2代表一组异常点，根据实际情况添加条件
               ---这里添加了三组异常点，当然多添加几个条件并不会影响结果
               where not (
                           point_time > (select point_time from invalid where invalid.number = 1 limit 1)
                       and point_time <= (select point_time from invalid where invalid.number = 2 limit 1)
                   )
                  or not (
                           point_time > (select point_time from invalid where invalid.number = 3 limit 1)
                       and point_time <= (select point_time from invalid where invalid.number = 4 limit 1)
                   )
                  or not (
                           point_time > (select point_time from invalid where invalid.number = 5 limit 1)
                       and point_time <= (select point_time from invalid where invalid.number = 6 limit 1)
                   )
               order by point_time)
select st_makeline(valid.the_geom)
from valid;

最后可以使用平滑函数来平滑一下轨迹线

...
select st_simplify(st_makeline(valid.the_geom),0.002)
from valid;