Track segment association algorithm based on statistical binary thresholds

1.   引言

合成孔径雷达(Synthetic Aperture Radar, SAR)作为一种全天时，全天候的遥感手段，在海洋观测，灾害监制，资源勘察，精细农业，地质测绘，政府公共决策，环境保护等方面有着越来越广泛的应用^[1]。但随着技术的发展与社会的进步，由于传统SAR只具备方位向和地距向两维的分辨能力，已然不能满足我们实际应用中的需求。干涉合成孔径雷达(Interference SAR, InSAR)具有大面积生成地表高度模型的能力，近些年来获得了广泛的应用。然而，在进行陡峭地形测绘时，由于叠掩和阴影等现象，干涉相位图中经常出现干涉系数比较低，不能进行相位解缠的区域，这些区域称为困难区域。阵列干涉SAR技术是InSAR技术的延伸，通过在交轨向布放多个阵元形成天线阵列，具备了高程向分辨能力，同时提高了多个通道之间的相参性，保证相干测量精度，从而实现观测场景的高精度3维成像。

利用3维点云数据进行目标重建的技术在激光雷达技术(Light Detection And Ranging, LiDAR)中广泛使用。在建筑类目标重建中，通过使用3D Hough变换或随机样本一致算法选择共面点，对分段建筑区域的表面拟合完成3维物体表面重建^[2]。大量无监督聚类方法^[3]、基于区域增长算法^[4]、图匹配技术^[5] 等用于建筑屋顶分割和重建。然而，阵列干涉SAR成像机理和LiDAR存在很大区别，常规LiDAR点云处理方法不太适用于阵列干涉SAR点云处理。

不过由于微波的散射特性和阵列干涉SAR的侧视模式，SAR 3维点云可以获取更多建筑立面的散射点信息和散射特性，不仅反映了建筑物的空间位置和结构形状信息，还反映了特定散射结构和散射材料特性。传统建筑物提取方法是结合SAR 3维点云建筑立面的散射点信息较丰富的特性，根据方位斜距投影密度和表面法向量来提取建筑立面^[6]，进而根据建筑物立面信息，使用密度聚类^[7]、区域增长^[8]的手段找到建筑物屋顶信息，从而完成建筑物点云的提取。但是由于阵列干涉SAR存在较多杂散点，而且场景中植被信息较多，根据表面法向量分辨立面的方法存在一定的困难，因此在进行建筑物提取之前，需要对3维点云数据进行相应的滤波处理。

经典LiDAR点云滤波方法有形态学滤波算法^[9]，基于坡度滤波算法^[10]，迭代加密三角网滤波算法^[11]等。形态学滤波和基于坡度滤波算法需要利用精细的高度信息，迭代加密三角网算法在迭代加密的过程中对误差更是十分敏感，因此经典LiDAR点云滤波方法不太适用于高程向误差较大的阵列干涉SAR。本文提出的滤波算法，不需要精细的高度信息，也没有迭代加密过程，对误差不敏感，主要利用建筑物区域点云密集且规整的性质，结合空间聚类种子生长手段就可实现点云滤波，从而为后续的建筑物提取重建提供便利，使得提取重建效果更加精确。

本文安排如下，本文的第2节主要介绍实验处理所需的数据集以及建筑物模型假设，第3节主要介绍本文提出的算法原理和处理流程，并在第4节使用真实的数据进行实验验证，在第5节中，通过和传统LiDAR滤波算法比较，验证本文算法的性能，并在第6节得到相应的结论。

2.   数据集以及建筑物模型假设

2.1   阵列干涉SAR点云数据

对阵列干涉SAR多个通道的2维复图像进行图像配准、幅度校正和相位校正，从而满足观测区域3维成像的基本要求。基于压缩感知的超分辨算法对观测区域进行高程向超分辨后，进行坐标转换，得到观测场景的原始3维点云。本文的数据采用的是中国科学院电子学研究所国内首次机载阵列干涉SAR实验的实际数据，也即图1所示数据。

图  1  原始3维点云

Figure  1.  The original three-dimensional point cloud

下载: 全尺寸图片 幻灯片

2.2   建筑物模型假设

在大多数情况下，建筑轮廓被假定为一个矩形多边形，而且边界跟踪算法通常确定规范的边界为直线段，这样就可以用一个多边形的形状代表建筑轮廓。在本文中，我们假设建筑物轮廓为矩形多边形且建筑物立面垂直于地面，如图2所示。虽然存在桥、塔以及球形等复杂建筑的情况，但矩形多边形的模型适用于大多数的建筑。

图  2  简单建筑物模型

Figure  2.  Simple building model

下载: 全尺寸图片 幻灯片

3.   基于空间聚类种子生长算法的建筑物点云滤波算法

3.1   算法原理

在阵列干涉SAR信号处理流程中对于从多幅2维图像中配准校正以及超分辨得到观测区域的原始3维点云分布的相关研究较多，可以参考文献[12–14]的算法进行处理。而对于3维点云的后处理研究较少，还处于人机交互阶段。由于原始3维点云中存在着大量的噪点以及虚假信息，在建筑物提取之前进行点云滤波显得极为重要，这也是本文的研究重点。

本文点云滤波算法原理主要如下：

(1) 对象的高程向高度以及邻域内点数

此场景下我们关注的建筑物高度是比植被要高的，该信息在进行滤波的时候可以很好被利用，而且植被点形成的表面较建筑物表面更为粗糙，所以在每一个建筑物点周围应该有相似高程的其他建筑物点，而植被则不一定。通过设定MinPoints可以去除一部分无效的杂散点，然后根据先验知识设置一个高程阈值就可以一定程度上区分建筑物点和低小植被点及其他地物。

(2) 对象的聚集程度和面积

由于建筑物的形状相较于植被来说具有一定的规则性，植被点云相对零散，聚集程度不如建筑物。同时建筑物由于窗户、阳台、建筑物楼顶边缘的凸起保护措施以及屋顶的不严格平面都会使得建筑物区域的点云聚集程度大大超过植被区域，因此对象的聚集程度可以作为分辨建筑物和植被的很重要的依据。

另外对象的面积是一个重要指标，经过高程-密度阈值滤波后，剩下的植被相对分布零散。可以设置一个面积的阈值，认为建筑物的面积理应在指定范围内。该面积大小的设定只需要依据场景的情况估算即可，需要部分先验知识，难度不大，而本文采用的面积阈值为50。

根据建筑物区域点云相对其他杂点高度较高以及密集且规整的性质，本文采用一种基于空间聚类种子生长算法进行点云滤波，最终达到去除杂点的目的。

3.2   算法流程

本文基于密度-高程图像的空间聚类种子生长算法能够去除植被以及地面其他地物的影响，为后续的建筑物提取重建提供便利，使得提取重建效果更加精确。其主要算法流程如下：

步骤1 数据分块(等网格分割)，基于密度与高程双重阈值生成密度-高程图像，去除一部分杂散点和低小地物。

首先将3维点云在地距-方位平面上投影，将不规则的SAR点云数据进行数据分块处理。现有的数据分块方法有：均匀条带分割、等格网分割、四叉树分割和自适应分块等^[15]，由于本步骤需要去除的仅仅只有离散杂点，所以比较适合数据分块的方式是等网格分割。然后通过向地距方位向做投影得到每个网格的密度值。将小于窗口阈值的格网认为离散噪点，直接去除就可以。而对于那些高于阈值的网格不做处理，虽然它们的密度值达到了要求，可是依然包含着一些小型的植被等其他点云数据。对于它们来说有一个比较普遍的现象就是它们在地面的高度上不会太高，因此我们选择使用每个网格的散射点的高程的平均值作为我们每个网格的高程值。利用高程阈值，去除一些小型地面目标以及低矮植被等点云数据，遍历所有格网，如果该网格存在散射点将该网格标记为1，而如果该网格为空，即网格内的散射点已经被去除或者根本不存在散射点，将该网格标记为0。这样我们就能得到原点云数据所对应的密度-高程图像。

步骤2 基于图像处理手段去除较为离散的散射点，填补建筑物空洞。

首先遍历整幅图像，利用像素周围的其他像素特性将不满足情况的离散点去除，由于建筑物比较规整，而杂点相对比较离散所以利用图像处理的手段可以将比较离散的一部分杂点去除，将一部分建筑物空洞利用形态学处理手段进行补充，从而得到较为完整的建筑物密度-高程图。

参考形态学图像处理具体方法如下：

(1) 对于灰度值为1的点云，只要它的上下左右4个临近方向的灰度值为0，那我们将其灰度值重新设置为0。对于灰度值为0的点云，如果它的4个方向灰度值都为1，那我们重新设置该点灰度值为1；

(2) 对于灰度值为1的点云，只要它的左上，右上，左下，右下4个方向的灰度值均为0，那我们将其点的灰度值重新设置为0。对于灰度值为0的点云，如果它的4个方向灰度值都为1，那我们重新设置该点灰度值为1；

(3) 基于形态学的内部空洞填充手段对密度-高程图像进行空洞填充，消除建筑物空洞。

步骤3 利用空间聚类种子生长算法，去除小面积区域，得到较为完整的建筑物区域。

相对于完整的建筑来说，植被相对零散，形成的对象面积较小，因此利用空间聚类种子生长算法去除小面积图斑的形式去除植被信息。

空间聚类种子生长的思想是在生长的过程中不关注会生成区域的形状，而是以目标元素内部的任意一个元素作为种子，只要符合我们的生长条件，就将其加入到种子队列中，这样不断地扩大种子区域，最终将整个区域的符合条件的种子区域全部找到。在生长的过程中对新增加的种子个数进行统计，每增加1个种子数，对应的面积数加1，最后得到该种子区域的总面积。该算法还要额外说明一下连通像素的定义。常规的情况下对于连通的搜索判断方式可以分为两种：四邻域搜索和八邻域搜索。图3展示了两种搜索方式的示意图，其中图3(a)是四邻域搜索示意，图3(b)是八邻域搜索示意，其中X表示不连通单元。

图  3  搜索中邻域的编码

Figure  3.  The encoding of the neighborhood in the search

下载: 全尺寸图片 幻灯片

四邻域搜索认为，只有以种子元素为中心的上下左右4个方向的相邻像素是与种子像素相连通的。八邻域搜索则认为，以种子元素为中心的周围8个方向相邻的元素均与种子像素是相连通的。如果使用数字代码来描述，那么四邻域中连通代码为0, 2, 4, 6，而八邻域中连通代码为0～7，本文采用八邻域搜索方式。

种子生长算法流程如下：

(1) 建立种子栈，将所有种子放入栈中；

(2) 从种子栈中取出一个种子，进行8方向搜索，将搜索到的种子点从种子栈中清除并放入新栈，每搜索到一个种子点对应的sum值加1；

(3) 从新栈中取出一个种子点，重复步骤2的操作，一直到新栈为空，此时获得的sum就是该区域的面积；

(4) 判断面积是否符合建筑物面积标准，如果大于阈值，保留该区域，否则去除；

(5) 判断种子栈是否为空，若为空，跳出循环，否则重复步骤2–步骤5，直到种子栈为空。

具体流程图如图4所示。

图  4  种子生长流程图

Figure  4.  Flowchart of seed growth

下载: 全尺寸图片 幻灯片

4.   实验结果和验证

4.1   利用密度高程双重阈值生成密度-高程图像

根据文献[6]的思路我们设置将网格划分为0.3 m×0.3 m，设置密度阈值为200 points/m²，然后根据先验知识，采用的高度阈值为5 m。

图5(a)为做完地距方位向投影后的点云数据，可以看出除了建筑物点云之外，点云中存在着较多的杂点。图5(b)为经密度高程双重阈值滤波后生成的密度-高程图像。

图  5  点云处理前后示意图

Figure  5.  Schematic diagram before and after cloud processing

下载: 全尺寸图片 幻灯片

通过对比图5(a)和图5(b)，可以看出经过密度高程双重阈值滤波之后，有一部分杂点已经被去除，但是从图6的放大图可以看出，依然存在一些比较离散的点，同时建筑物区域内部出现了空洞，所以需要进一步去除离散点并且填补建筑物内部的空洞。

图  6  处理后局部放大图

Figure  6.  Local enlarged image after processing

下载: 全尺寸图片 幻灯片

4.2   基于图像处理的离散点去除以及空洞填补

对于密度-高程图像中存在的离散点和建筑物空洞，可以利用形态学图像处理手段进行处理。图7为参考形态学图像处理前后的对比情况。

图  7  图像处理前后的密度-高程图像

Figure  7.  Density-elevation image before and after image processing

下载: 全尺寸图片 幻灯片

从图7(a)和图7(b)的对比可以看出，经过图像处理后，有一些离散点已经被很好地去除掉，而且从图8(a)和图8(b)的局部放大图中可以看出建筑物空洞得到了很好地填充。但是仍然会有一些小面积图斑存在，如图8(b)的红色圆圈处，这些可能是比较密集的植被所形成，需要下一步的去除。

图  8  图像处理前后的局部放大图

Figure  8.  Local enlargement before and after image processing

下载: 全尺寸图片 幻灯片

4.3   基于空间聚类种子生长算法去除小面积图斑

对于小面积图斑的处理，可以利用空间聚类种子生长算法计算每个种子对应的面积，设置一个面积的阈值去除小面积图斑，本文采用的面积阈值为50。处理效果如图9所示。

图  9  小面积区域去除前后的密度-高程图像

Figure  9.  Density-elevation images before and after removal of small area

下载: 全尺寸图片 幻灯片

从图9(a)和图9(b)对比可以看出经过空间聚类种子生长算法处理后，小面积图斑被去除，没有了杂点的影响，建筑物区域被很好地提取出来。图10为滤波后的建筑物区域与实验区域光学图像中的建筑物区域的对照图，从图中可以看出本文算法可以有效地实现点云滤波，去除场景点云中的杂点，从而证实了算法的可行性。

图  10  滤波后结果与光学图像对照图

Figure  10.  Comparison of filter results with optical image

下载: 全尺寸图片 幻灯片

5.   质量评估

为了评估建筑物提取的质量，本文采取了一种逐点比较的方法，滤波后得到的建筑物散射点和参考建筑物进行对比，完成性能评估。由于实际实验中我们并没有真正代表建筑物区域的准确数据，所以在滤波有效性的精确评估是无法做到的。我们假设距离建筑轮廓线中心2个像素内的点为有效的散射点^[16]。原始3维点云中被分为两类：建筑物散射点或者非建筑物散射点。滤波后的散射点属于建筑物的定义为T_p，不属于建筑物的定义为F_p。类似属于建筑物的散射点被检测为非建筑物散射点的定义为F_N。滤波算法的性能可以参考如下的评估准则^[6]：

完整度Completeness代表着所有杂点去除的比例，完整度越高，说明在点云数据中杂点去除得越干净；正确率Correctness则表示为杂点去除正确的概率，正确率越高，说明杂点识别得越清楚，去除杂点越有效；而质量评价Quality代表我们整体算法的完成质量，即融合了完整度和准确率两者信息，质量评价参数越高则代表杂点去除的精确度和完整度越好，因此我们使用质量评价参数Quality作为对滤波算法性能分析的指标。

利用国内首次机载阵列干涉SAR实验的实际数据，该区域包含的点云数目为104601个。目前对于阵列干涉SAR点云处理还停留在人机交互的阶段，因此采用人机交互完成的滤波结果作为标准用来与各个算法得到的结果作一个比较，并结合上述评价准则对算法性能做一个检验。图11(a)是传统的人机交互完成的滤波结果，图11(c)，图11(e)，图11(g)分别为本文算法、形态学滤波算法、坡度滤波算法的处理结果图。

下载: 全尺寸图片 幻灯片

通过表1的对比可以看出本文算法相比于传统的LiDAR滤波算法能够更好地去除杂点，完成质量达到94.81%，达到了滤波的目的。

表  1  滤波算法质量评价表

Table  1.  Filter algorithm quality evaluation table

滤波算法	After filtering	Tp	Fp	FN	Completeness (100%)	Correctness (100%)	Quality (100%)
本文算法	64858	63753	1105	2381	96.40	98.30	94.81
形态学滤波	81443	64032	17441	2102	96.82	78.62	76.64
坡度滤波	84039	64395	19944	1739	97.37	76.63	75.07

下载: 导出CSV 

| 显示表格

6.   结论

由于阵列干涉SAR存在杂散点较多、高程向误差较大的问题，传统LiDAR的点云滤波方法效果并不太理想。针对上述问题，本文提出了基于空间聚类种子生长算法的阵列干涉SAR点云滤波算法。国内首次机载阵列干涉SAR实验数据的处理结果表明本文算法具有良好的滤波效果，杂点去除质量达到了94.81%，验证了本文算法的有效性。