背景

在处理GPS数据时，记录中往往会有很多时间戳和距离都接近的数据，一方面浪费了较多的存储空间，另一方面造成所要表达的图形不光滑或不符合标准。因此要通过某种规则，在保证矢量曲线形状不变的情况下， 最大限度地减少数据点个数，这个过程称为抽稀。

常用的抽稀算法，包括但不限于（这部分内容来源于网络）

• 步长压缩算法：沿连续曲线每隔一定的步长抽取一点，其余点全部压缩掉，然后在相邻抽取点间用直线连续或曲线拟合逼近。
• 线段过滤算法：当某一段的长度小于某一过滤值时，就以该段的中点代替该段，如同此段的两端退化到中点一样。
• 道格拉斯-普克（Douglas-Peuker)算法：亦称为拉默-道格拉斯-普克算法、迭代适应点算法、分裂与合并算法)是将曲线近似表示为一系列点，并减少点的数量的一种算法。

PostGIS内置了部分常用算法，我们进行一下简单的测试。

使用PostGIS函数

使用到的PostGIS函数如下

• ST_Simplify
• ST_SimplifyVW
• ST_SimplifyPreserveTopology

轨迹抽稀

函数ST_Simplify()基于Douglas-Peucker算法进行轨迹抽稀，具体原理这里不解释。但是有一点要注意的是，当使用WGS84时，算法的容差参数单位是度。

经纬度在不同地区，每度距离差是不同的，如果假定地球是完美的球体（这样假设误差不是很大）的话，纬度为 B 的地区：

• 纬度变化一度，球面南北方向距离变化：πR/180 ........111.7km
• 经度变化一度，球面东西方向距离变化：πR/180cosB ....111.7cosB

比如北京 B = 40、cosB = 0.766，经度变化1度，则东西方向距离变化 85.567km。

为了简便计算，我们假设中国范围内，1度约等于100公里。

Douglas-Peucker 算法

抽稀后，因为拓扑形状变化，轨迹距离会有损失。我们希望的结果是，空间比例尽量小，距离比例尽量大。

测试，当容差距离分别为1000米，100米，10米的情况下，抽稀结果的对比。

with simplified_trajectory as 
(
	-- 1 degree = 100 km
	select 
		tr.tid, tr.dt,
		st_npoints(traj) as np_orgin,
		(st_length(traj::geography)/1000) ::int as km_orgin,
		-- 0.01 degree = 1000 m
		st_npoints( ST_Simplify(traj,0.01,true) ) as np_01,
		(st_length(ST_Simplify(traj,0.01,true)::geography)/1000) ::int as km_01,
		-- 0.001 degree = 100 m
		st_npoints( ST_Simplify(traj,0.001,true) ) as np_001,
		(st_length(ST_Simplify(traj,0.001,true)::geography)/1000) ::int as km_001,
		-- 0.0001 degree = 10 m
		st_npoints( ST_Simplify(traj,0.0001,true) ) as np_0001,
		(st_length(ST_Simplify(traj,0.0001,true)::geography)/1000) ::int as km_0001,
		1 as endflag
	from demo.t_taxi_trajectory tr
)
select st.tid, st.dt, st.np_orgin, st.km_orgin,
	(st.np_01 / st.np_orgin::float)::decimal(4,2) as np_01_pct,
	(st.km_01 / st.km_orgin::float)::decimal(4,2) as km_01_pct,
	(st.np_001 / st.np_orgin::float)::decimal(4,2) as np_001_pct,
	(st.km_001 / st.km_orgin::float)::decimal(4,2) as km_001_pct,
	(st.np_0001 / st.np_orgin::float)::decimal(4,2) as np_0001_pct,
	(st.km_0001 / st.km_orgin::float)::decimal(4,2) as km_0001_pct
from simplified_trajectory st
where 1=1	
	and st.dt = '2008-02-06'
	and st.tid  in (28,10012)
	order by st.tid, st.dt
;

Visvalingam-Whyatt 算法

ST_SimplifyVW 内置的另一个抽稀算法，没有DP那么有名。使用方法一样，这里不写语句了，有兴趣看官方文档。

测试中发现，VW算法中的容差和DP不太一样，原因未知。

步长压缩算法

对原始轨迹点使用逢十抽一的步长算法，即每10个原始点中抽一个。

-- 常规做法应该从原始GPS点数据开始，但是为了演示和测试，选择把轨迹还原成GPS然后再抽稀
with dump_pt as
(
	select tid, 
		(st_dumppoints(tr.traj)).path[1] as sn, 
		(st_dumppoints(tr.traj)).geom as pt 
	from demo.t_taxi_trajectory tr
),
-- 每10个原始点划分为一个组
rank_pt as
(
	select 
		tid, to_char(to_timestamp(st_m(pt)), 'yyyy-mm-dd') as ts ,
		st_setsrid( st_makepointm( st_x(pt), st_y(pt), st_m(pt) ), 4326) pt,
		row_number() over (partition by tid, sn/10 order by sn) as rn
	from dump_pt 
)
select rpt.tid, rpt.ts,
	st_makeline( array_agg(rpt.pt order by rpt.ts )::geometry[] ) as traj
from rank_pt rpt
where rpt.rn = 1
group by rpt.tid, rpt.ts
;

DP算法中容差阈值的影响

  tid  |     dt     | np_orgin | km_orgin | np_01_pct | km_01_pct | np_001_pct | km_001_pct | np_0001_pct | km_0001_pct
-------+------------+----------+----------+-----------+-----------+------------+------------+-------------+-------------
    28 | 2008-02-06 |      947 |      180 |      0.01 |      0.94 |       0.03 |       0.97 |        0.15 |        0.98
 10012 | 2008-02-06 |      687 |        5 |      0.00 |      0.00 |       0.00 |       0.00 |        0.14 |        0.40
(2 rows)

从这个样本数据中看出

• 长原始距离（947km）
  • 容差1000米时，坐标点个数减少为1%，距离减少为94%；
  • 容差100米时，坐标点个数减少为3%，距离减少为97%；
  • 容差10米时，坐标点个数减少为15%，距离减少为98%；
• 短原始距离（5km）
  • 容差1000米时，坐标点个数减少为0%，距离减少为0% （偏差严重）；
  • 容差100米时，坐标点个数减少为0%，距离减少为0% （偏差严重）；
  • 容差10米时，坐标点个数减少为15%，距离减少为40%；

所以，当原始轨迹距离比较长的情况下，用100米容差的效果很好。如果原始距离很短，需要更小的容差。

轨迹抽稀效果对比

把一条样本轨迹用不同抽稀算法处理，然后比较容量和准确度。

drop table if exists demo.t_02;
create table if not exists demo.t_02 (id text, traj geometry(linestringm, 4326));


insert into demo.t_02 
select 'src', traj from demo.t_taxi_trajectory where tid=28 and dt='2008-02-02'
union all
select 'dp', ST_Simplify(traj,0.001,true)  from demo.t_taxi_trajectory where tid=28 and dt='2008-02-02'
union all
select 'vw', ST_SimplifyVW(traj,0.00001)  from demo.t_taxi_trajectory where tid=28 and dt='2008-02-02'
;

--
with dump_pt as
(
	select tid, 
		(st_dumppoints(tr.traj)).path[1] as sn, 
		(st_dumppoints(tr.traj)).geom as pt 
	from demo.t_taxi_trajectory tr
	where tr.tid in (28) 
		and tr.dt = '2008-02-02' 
),
--
rank_pt as
(
	select 
		tid, to_char(to_timestamp(st_m(pt)), 'yyyy-mm-dd') as ts ,
		st_setsrid( st_makepointm( st_x(pt), st_y(pt), st_m(pt) ), 4326) pt,
		row_number() over (partition by tid, sn/10 order by sn) as rn
	from dump_pt 
)
--
insert into demo.t_02 
select '10-1',
	st_makeline( array_agg(rpt.pt order by rpt.ts )::geometry[] ) as traj
from rank_pt rpt
where rpt.rn = 1
group by rpt.tid, rpt.ts
;

算法	轨迹点数量	轨迹公里数
原始轨迹	387	190
DP 0.001	130	188
VW 0.00001	107	174
步长10抽1	39	140

效果对比图如下（结论仅对本条数据有效），可见DP的准确性最好，步长10抽1最差。

拓扑抽稀

Douglas-Peucker 算法

和ST_Simplify类似，函数ST_SimplifyPreserveTopology也使用DP算法进行抽稀。但是ST_SimplifyPreserveTopology可以保证抽稀后的数据是一个合法的拓扑图形。

我不是特别明白，但是我的理解

• 不支持M值，所以不能用来抽稀基于LinestringM的轨迹数据。
• 是Polygon对象抽稀，最好使用这个函数。（别问为什么，因为我也不明白）

-- 容差0.001的DP拓扑抽稀
select 
	st_multi(st_simplifypreservetopology(fence_polygon, 0.001)) as fence_polygon
from demo.t_area_fence where id = 510107;

拓扑抽稀效果对比

把一条样本围栏用不同抽稀算法处理，然后比较容量和准确度。

警告 对国家领土边界进行数据抽稀，仅仅作为验证技术可行性的研究行为，不能用于其他任何场景！！！

drop table if exists demo.t_03;
create table if not exists demo.t_03 (id text, fence_polygon geometry(multipolygon, 4326));


insert into demo.t_03 
select 'src', fence_polygon from demo.t_area_fence where id = 510107
union all
select 'dp', st_multi(ST_Simplify(fence_polygon,0.001,true)) from demo.t_area_fence where id = 510107
union all
select 'dppt', st_multi(st_simplifypreservetopology(fence_polygon, 0.001)) as fence_polygon from demo.t_area_fence where id = 510107
;

算法	围栏点数量
原始围栏	830
DP 0.001	126
DP拓扑 0.001	127

效果对比图如下（结论仅对本条数据有效），红色方框部分有细微差异，技术研究性质的GIS业务场景下可考虑使用。

警告 对国家领土边界进行数据抽稀，仅仅作为验证技术可行性的研究行为，不能用于其他任何场景！！！