利用双窗口检测和粒子滤波提取高分辨率SAR图像道路中心点

程江华 关永峰 库锡树 孙即祥

(国防科技大学电子科学与工程学院 长沙 410073)

1 引言

在各类星载或机载成像传感器中，SAR具有卓越的空间信息获取能力，能够全天候工作，在国防、环境等方面具有重要的战略意义，自诞生之日起就受到了普遍关注。特别是近几年来，随着成像技术的发展，从不同星载或机载SAR传感器中获得了大量的高分辨率 SAR 图像，这些图像描述的地表信息更丰富，地物目标的细节更清晰。道路作为典型的地物目标，是地理信息系统(GIS)的重要组成部分，在城市规划、交通控制、GIS数据库更新、应急响应等诸多领域有着极其广泛的应用。

高分辨率道路提取方法主要分为自动和半自动两大类。自动提取方法将道路的先验知识融合进算法中，大致分为线特征提取和连接两个步骤，是未来的发展方向及目标。虽然分辨率图像表现的细节更加清晰，但是各种噪声干扰对提取的影响也加剧；另外，算法需要结合具体某类图像道路的先验知识，设置大量的参数，算法的稳健性、适用性不强；一般情况下，处理结果需要后期人工检验和调整。就目前的图像处理及模式识别技术水平，实现高分辨率SAR图像道路的自动提取非常困难。半自动提取方法能够将机器的快速计算和人的解译技巧进行有效结合[1]，是当前研究的热点和主流。

半自动提取方法一般分为局部道路特征提取和全局道路跟踪两个步骤[1-4]。局部道路特征信息主要有局部道路区域的方向、宽度和中心点坐标等；全局道路跟踪是根据局部道路特征提取计算获得的方向、宽度和中心点坐标等信息估计下一个道路中心点，并将获得的中心点连成道路网络。此类方法中，文献[1]根据相邻矩形模板间均方误差最小的准则，实现了主干道中心点的跟踪提取；文献[2]根据目标模板间距离的远近给参考模板组设定加权系数，提出当前模板与参考模板组相匹配的策略，克服道路中噪声对模板匹配的影响；文献[4]交替使用扩展卡尔曼滤波(正常情况)和粒子滤波(遇到障碍或路口)进行道路中心点跟踪连接。该类方法存在的普遍问题是：(1)为确定道路方向、中心点坐标和道路宽度等初始信息，需手动输入两个以上初始点，且后续跟踪效果易受初始点影响；(2)模板类型及尺寸固定，不能适应不同类型(高速公路、城区主干道、次干道、一般街道等)道路跟踪；(3)一般针对加性噪声模型的光学遥感图像，不能满足低信噪比SAR图像乘性噪声抗干扰要求；(4)方向、中心点坐标、宽度等观测值和估计值计算存在相关关系，一旦跟踪失败，不能自动纠正；(5)采用固定的步进长度进行跟踪，遇到护栏、绿化带、车辆、天桥等地物干扰时，不能跳过，需要人工引导，增加了人机交互次数。

为减少初始点的数目对后续跟踪的影响，提高算法的抗干扰性及适应性能力，本文提出了针对高分辨率SAR图像道路半自动提取新的框架。该方法分为局部道路检测和全局粒子滤波道路跟踪这两个步骤。本文结构安排为：第2节详细阐述局部道路检测，第3节介绍全局粒子滤波道路跟踪，第4节给出实验结果及分析，第5节得出结论。

图1 高分辨率SAR图像道路边缘情况

2 局部道路检测

在高分辨率SAR图像中，道路不再表现为线特征，而是呈现出由双边缘包围的较长的暗区域(图1(a))。因此，可通过寻找双边缘平行线或搜索狭长的暗同态区域检测道路。然而，由于强乘性噪声及护栏、绿化带、车辆、天桥等地物的干扰，在部分区域道路和周围地物区分不明显，有时人眼也难以分辨(图1(b))；道路与周围建筑物粘连，使得原本连续的道路双边缘变成单边缘(图1(c))，路旁建筑物、树木等遮挡使得部分路段甚至没有边缘(图1(d))。因此，采用常规的边缘检测方法很难有效解决道路中的干扰及遮挡问题。

本文提出采用双窗口模型，在外窗口区域采用非线性结构张量[5]提取局部道路区域的方向，在内窗口采用最佳搜索法计算出道路宽度和中心点坐标，如图2所示。

2.1 双窗口模型

在传统道路局部区域计算中，通常采用长宽固定，方向可变矩形窗口模型[1]，算法简单实用，但不能适应不同宽度情况道路提取；或采用长度固定、宽度和方向可变的矩形窗口模型[2]，能够根据道路的不同情况进行自适应窗口大小，但是在没有边缘的道路遮挡区域，窗口的宽度会无序扩大(图1(d))，脱离实际的道路宽度。

图2 本文提出的SAR图像道路网提取流程图

本文结合固定窗和可变窗的优点，提出采用双窗口模型(图3(a))。外窗口形状为方形，大小固定，便于快速计算道路方向。内窗口转向为外窗口计算得到的方向，形状为矩形，长度固定，宽度可变，规定内窗口不超过外窗口区域，避免在无边缘的情况下窗口的宽度无约束扩大。

图3 双窗口模型示意图

2.2 方向计算

外窗口大小位置确定后，在窗内计算道路的主方向。传统局部主方向计算通常采用边缘点梯度方向直方图统计求极值的方法。使用较为广泛的有：边缘方向直方图法[5](EDH)，以及其改进形式 对称边缘方向直方图[6](SEDH)等。这类方法的计算效果易受边缘检测性能的影响，特别是在乘性相干斑噪声干扰的低信噪比SAR图像中，边缘方向直方图的极值不明显，可能出现多个极值现象。而线性结构张量[7]能保持图像中的重要结构信息，可用于道路主方向提取。但该方法在高斯滤波的同时，也降低了图像边缘的对比度。为克服其不足，本文针对SAR图像乘性噪声干扰情况，提出采用非线性结构张量[8]计算局部主方向。

定义图像I在点(i,j)处的梯度为：∇Ii,j=((Ii,j)x,(Ii,j)y)T，则在该点结构张量表示为

G(J)的特征值为

假设大的特征值λ1对应特性向量ω1，小的特征值λ对应特性向量ω，反之亦然。则

那么点(i,j)处的方向：θi,j=arctan(2u'12/在外窗口区域内，局部道路主方向的计算流程如下：

(1)先利用式(2)计算每点的结构张量，然后依次计算对应的特征值和特征向量，通过特征向量，计算得到该点的方向；

(2)统计各点的方向值，形成方向直方图；

(3)求方向直方图中的峰值，该峰值对应外窗口区域内的主方向。

2.3 宽度和中心点计算

道路的宽度信息决定了中心点的坐标位置。以往半自动道路提取方法中，道路的宽度信息在跟踪的初始阶段由人工初始输入三点来确定[3]，存在跟踪过程易受初始点影响的问题；或是在跟踪开始时，根据初始输入的两点，结合道路的双边缘特征，自动计算道路宽度[2]。该方法需要满足双边缘条件，不能应用于受遮挡的单边缘或无边缘情况。本文采用内窗口先平移后扩大方式，在外窗口范围内自动搜索道路的宽度，具体流程如下：

(1)初始化 以外窗口计算得到的局部道路主方向为内窗口的方向，以外窗口的中心点为内窗口的初始中心点，形状为矩形；

(2)平移 以原始中心点为中心，垂直于道路方向，固定步进长度，在外窗口范围内进行平移，求解系列平移窗口的方差均值；

(3)扩大 以步骤(2)平移中求得系列平移窗口方差均值的最小值窗口为扩大的初始窗口，垂直于道路方向，固定步进长度，在外窗口范围内进行内窗口双向扩大，求解系列扩大窗口的方差均值；

(4)计算 以步骤(3)平移中求得的系列扩大窗口方差均值的最小值窗口为该区域道路的最佳搜索匹配窗口，以该窗口的宽度为道路的宽度，以该窗口的中心为道路的中心。

本文方法无需人工输入多点，无需满足双边缘条件。由于采用了具有一定面积的矩形区域，因此可降低道路上干扰物对搜索结果的影响。另外，给内窗口的搜索设定了范围，这样不至于在道路被遮挡处内窗口无限制扩大。

3 全局粒子滤波跟踪

将道路中心点的跟踪看成一个时间序列过程，在输入初始点后，跟踪算法根据初始状态预测下一状态，到达下一状态后，根据观测值更新预测值。这样周而复始的循环，直到满足停止条件，跟踪过程结束。

假定动态时变系统预测方程和观测方程分别为

其中sk表示第k个道路中心点的状态值，uk表示过程噪声；zk表示第k个道路中心点的观测值，vk表示测量噪声。

定义状态向量sk为3维向量，按照图4所示，预测方程具体形式为

其中()为中心点坐标，θk表示中心点方向，dk为中心点步进间距，(i=1,2,3)为过程转换噪声。观测方程具体形式为

图4 道路中心点状态转换示意图

根据贝叶斯跟踪理论，如果已知初始概率密度函数p(s0/z0)=p(s0)，则状态预测方程和状态更新方程分别为

粒子滤波是一种非线性、非高斯动态系统的实时推理算法，采用了蒙特卡洛思想近似公式(9)，式(10)，其核心思想是利用一系列随机样本的加权和表示后验概率密度，通过求和来近似积分操作。定义为粒子，则后验概率密度方程可近似为

假定观测值呈正态分布，则观测值的似然函数为

其中li为粒子与观测值的欧氏距离。则全局粒子滤波道路中心点跟踪的流程如下：

(1)初始化 以输入的一个道路中心点为中心，形成内外窗口。外窗口计算局部方向，内窗口根据计算方向值进行旋转，并估计道路宽度。在道路宽度范围内，以输入点为中心，随机产生N个粒子，并初始化权值=1 /N,i=1,…,N。

(2)预测 由于重要性密度函数q(sk/z1:k)选择根据式(7)预测下一状态1,2,…。

(3)更新 根据k时刻的测量值zk(由第 2节计算得到)，更新每个粒子的权值，并归一化为。

(4)输出 根据式(11)和步骤(3)得到的权值，求得期望，作为当前道路中心点的滤波输出。并计算有效粒子个数=效粒子数门限)，则转入步骤(5)进行重采样；否则，转入步骤(2)，进行下一中心点的预测。

(5)重采样 保持粒子数N不变，复制权值大的粒子，代替权值小的粒子。转入步骤(2)，进行下一中心点的预测。

在跟踪的过程中，计算测量值zk时，存在由于遮挡或路口道路宽度太宽，造成外窗口内无边缘方向。本文采用可变步进间距dk克服其影响。无遮挡时，dk为常量；遇到遮挡时，保持原来的方向，dk增大跳过遮挡物；当连续跳跃4次以上时，则认为跟踪失败，将控制权交给用户，重新初始化起点，进行跟踪。另外，当跟踪器到达图像边缘时，也将控制权交给用户，重新跟踪。

4 实验结果及分析

为了测试本文方法的有效性，选择3幅分辨率为1 m的机载高分辨率SAR图像切片进行实验。图5(a)外窗口中，存在明显的双边缘道路区域。图5(b)为采用边缘方向直方图法计算得到的结果，图5(c)为采用本文提出的结构张量直方图得到的结果。从图5(b)中可以看出，在乘性噪声干扰的低信噪比SAR图像中，直接计算图像的边缘点梯度方向得到的直方图，存在多个极值，且极值不明显。而结构张量由于在计算梯度的同时保存了边缘信息，可有效提取道路方向。图5(c)中存在明显的双峰，且峰值之间相差约 180°的相角，进一步证明外窗口区域中存在双边缘。

图6切片1中部分区域道路上带有车辆、铁护栏等高亮噪声干扰(窗2，窗3)，部分区域具有一定曲率(窗4，窗5)。实验结果表明，本文方法中外窗口能够有效检测出道路的方向，内窗口能够根据方向进行旋转，并根据方差值进行平移和扩大。图7切片2中有部分区域(窗3，窗4)被道路旁的高大建筑物遮挡，造成窗3计算得到的结果与前一窗2的结果相差甚远，窗4中无主方向，这些判断为方向提取失败。本文根据道路方向变换比较平缓的特性，采用可变步长的方法。在有边缘方向区域，固定步长为200；遇到遮挡等情况，步长每次增加40，方向保持不变，当步进到窗口5位置时，重新捕获到了道路区域。图6，图7中外窗口大小均为120×120，内窗口长度为80，宽度根据道路方差均值自动调整。表1为切片1，切片2道路中心点跟踪过程数据。

图6 切片1道路中心点跟踪图

另外，传统跟踪方法[1-4]输入需要 2-3个初始点，用于确定道路的起点、方向及宽度等信息。而本文只需输入1个起始点，道路的方向和宽度自动计算。且在后续的跟踪过程中，道路方向及中心点坐标的观测值与估计值相独立，这样一旦跟踪失败，能够自动纠正。由于 SAR图像道路上各种干扰存在，其灰度值变换较大，依靠灰度模板匹配的传统方法[1-4]实现跟踪很难奏效，且需要不断地比对、更新模板[2]。本文采用了双窗口方法，无需依靠模板匹配。传统跟踪[1,3]方法遇到障碍物遮挡时，需要人工引导重新指定方向点，交互次数多。本文方法能够跳过遮挡，重新实现自动跟踪，交互次数少。

图7 切片2道路中心点跟踪图

5 结论

道路上各种干扰物的存在以及道路周围各种地物遮挡的影响，使得传统的基于模板匹配的遥感图像道路提取方法在高分辨率 SAR图像道路提取中不再适用。本文根据高分辨率SAR图像的特点及道路干扰物具有方向性的特征，提出了双窗口局部检测和粒子滤波全局跟踪相结合的道路提取方法。实验结果表明，本文方法可有效降低路面高亮地物对道路局部检测的干扰，能够跳过遮挡物的影响，无需进行模板匹配，有效实现道路中心点的跟踪。

表1 切片1，切片2道路中心点跟踪数据

[1]Kim Tae-jung,Park Seung-Ran,Soo Jeong,et al..Tracking road centerlines from high resolution remote sensing images by least squares correlation matching[J].Photogrammetric Engineering＆Remote Sensing,2004,70(12): 1417-1422.

[2]Zhou Jun,Bischof W F,and Caelli T.Road tracking in aerial images based on human-computer interaction and Bayesian filtering[J].Photogrammetry＆Remote Sensing,2006,61(2):108-124.

[3]Lin Xiang-guo,Liu Zheng-jun,Zhang Ji-xian,et al..Combining multiple algorithms for road network tracking from multiple source remotely sensed imagery: a practical system and performance evaluation[J].Sonsors,2009,9(2):1237-1257.

[4]Movaghati S,Moghaddamjoo A,and Tavakoli A.Road extraction from satellite images using particle filtering and extended kalman filtering[J].IEEE Transactions on Geoscience and Remote Sensing,2010,48(7): 2807-2817.

[5]Jain K A and Vailaya A.Image retrieval using color and shape[J].Pattern Recognition,1996,29(8): 1233-1244.

[6]唐亮,谢维信,黄建军,等.一种新的道路描述子: 对称边缘方向直方图[J].电子学报,2005,33(1): 7-11.Tang Liang,Xie Wei-xin,Huang Jian-jun,et al..A new method for main roads description: symmetrical edge orientation histogram[J].Acta Electronica Sinica,2005,33(1):7-11.

[7]张道兵,张慧,张正,等.一种稳健的道路主方向提取算法[J].光子学报,2007,36(S1): 326-330.Zhang Dao-bing,Zhang Hui,Zhang Zheng,et al..A Robust main direction extraction method for road segment[J].Acta Photonica Sinica,2007,36(S1): 326-330.

[8]Brox T,Weickert J,Burgeth B,et al..Nonlinear structure tensors[J].Image and Vision Computing,2006,24(1): 41-55.

[9]鉴福升,徐跃民,阴泽杰.多模型粒子滤波跟踪算法研究[J].电子与信息学报,2010,32(6): 1271-1276.Jian Fu-sheng,Xu Yue-min,and Yin Ze-jie.Research of multiple model particle filter tracking algorithms[J].Journal of Electronics＆Information Technology,2010,32(6):1271-1276.