ISSN 
-
SAR图像具有全天时、全天候以及信息量大等特点，因此被广泛用于监测建筑物损毁、滑坡、冰雪灾害和洪涝灾害等领域。目前在道路损毁信息提取中，主要采用光学遥感图像，并获得了较多的研究成果。但光学遥感图像在特殊条件下，易受时间和气候等因素的影响，图像获取的实时性和成像质量无法保证，因此研究基于SAR图像的道路损毁信息提取具有研究价值和实用价值。目前与SAR图像道路损毁相关的研究包括道路损毁、建筑物垮塌和滑坡等信息的提取[1]。
基于SAR图像的道路提取方法已有大量研究，通常分为道路局部检测和道路全局连接两个步骤。局部检测的方法包括边缘检测和FCM分割等，道路全局连接的方法包括MRF模型、遗传算法、粒子滤波和动态规划等[2-4]。这些研究主要集中在城区或乡村环境中，取得了一些成果，但通常处理过程较为复杂，提取效果容易受到前期处理结果的影响。目前针对山区环境的道路提取研究较少，由于山区环境中阴影和地形因素等的干扰，进行边缘检测时目标误判和断裂情况较多，增加了全局连接的处理难度。因此本文选择粒子滤波方法进行道路提取，并根据山区环境进行了图像预处理、边缘特征方向获取和粒子权重设计等方面的改进，通过局部最优的求解方法，降低了道路提取的难度，实现了山区环境中道路的准确提取。
道路损毁的检测方法主要包括：视觉判读法、边缘检测法、光谱纹理分析法、面向对象法和复合分析法[5]。其中，视觉判读法为人工辅助方法，准确度最高，但效率较低；面向对象法利用了多尺度分割的思想，参数设置及算法复杂度较高；复合分析法采用多种数据进行联合分析，具有较高的实用价值[6]。目前该方面的研究仍主要集中于光学遥感[7-8]，SAR图像领域中的相关研究鲜见报道。根据现有的建筑物损毁提取和滑坡提取的成果及研究方法，文献[9]通过在道路区域内利用变化检测进行了道路损毁提取。
本文在多时相SAR图像中实现了道路损毁信息提取，其提取过程包括两个步骤：1) 提取损毁前的基准SAR图像区域内的道路中心线；2) 以步骤1)的道路中心线为基础，构建道路损毁模型，在配准后的待检测SAR图像中进行道路损毁检测，得到道路损毁的坐标、长度和损毁程度等信息。
HTML
-
在高分辨率SAR图像中，道路通常表现为具有一定长度和方向的低灰度暗黑条带，与道路背景区具有一定对比度，方向变化缓慢，多条道路相互连接形成道路网络。但是在实际成像条件下，SAR图像中存在严重的乘性相干斑噪声，具体表现为同质区域的像元强度存在涨落，从而导致了道路与背景区域对比度不明显、边缘模糊等现象。
本文的研究范围主要集中在山区，除以上由SAR图像成像特点所导致的复杂性外，山区中地物的分布对道路提取造成了较大的影响，如由于地形影响导致SAR图像异质区域较多，当进行边缘特征提取时，恒虚警分割无法直接有效地区分道路边缘与非道路边缘。此外一些与道路特征类似的线性结构的存在，也给道路连接造成了较大的干扰。因此在传统道路提取方法的基础上，山区环境中的道路提取需要充分利用道路的灰度、边缘和方向信息，依靠信息之间的互补性完成道路连接过程。道路中心线提取流程如图 1所示。
-
对于SAR图像中的道路提取，通常要进行一定的滤波处理以降低相干斑噪声对后续操作的影响。一种滤波思路认为道路是在长度方向上具有一定排列的低灰度值像素点，因此可利用该特征进行滤波器设计以便得到较好的结果。其中，文献[10]以标准FROST滤波器(SFF)为基础提出了一种改进方案，该方案综合利用了方向和轮廓等线特征信息，因此被称为方向性改进FROST滤波器(DMFF)。DMFF不仅能有效地提高道路对比度，平滑非道路区域，还能增加道路与河流及阴影的区分能力。
DMFF的核函数采用一维指数函数(GEF)描述，它具有更好的灵活性。它可以表示为高斯函数、指数函数和矩形函数等多种形式，因此能更好地描述空域上局部道路的特征。其定义为：
式中，参数a与GEF的宽度有关，a越小，则轮廓宽度越大；参数λ控制GEF的平坦度，λ越大，GEF的形状越平坦，当λ=1时，GEF则为高斯函数。
在式(1)的基础上，选择二维平面上的旋转GEF作为DMFF的核函数，其中方向k上的核函数以fk(x, y)表示，其参数a、λ要参考道路的宽度等属性进行选择。以7*7的滤波窗口为例，对当前像素点(x, y)，经过滤波后其像元值定义为：
-
文献[3]针对SAR图像的乘性噪声模型提出了一种具有恒虚警的边缘检测方法，它主要利用异质区域中像素点的灰度均值比进行检验。与传统的检测方法相比，该方法对SAR图像中的边缘检测具有更好的适应性和稳定性。
Duda算子[11]为ROA算子的一种改进方法，其检测模型利用道路中心灰度低于两边背景灰度的特点，结合了区域内部的均匀度进行道路边缘特征检测，因此具有更好的检测效果。它的检测窗口分为0°、45°、90°度和135° 4个方向，对每一个当前像素点计算出4个检测窗口对应的响应值，取它们的最大值作为该点的边缘强度值，并进行阈值分割。其计算过程如式(3)～式(7)所示：
式中，R1和R2分别为中心区域与左右背景区域的灰度均值比；Fi(x)和G(x)分别代表边缘对比度分值函数和中心区域均匀度函数；T1为对比度范围阈值，可以根据图像中道路区域对比度进行选择；H(x)为对应坐标点及相应方向上的响应值。本文直接以边缘强度值作为下一步粒子滤波的输入。
-
边缘方向表征了道路的走向以及道路像素之间的上下文联系，是粒子滤波权重设置中的重要参数。边缘方向可以根据Duda算子的最大响应值对应方向进行确定，但该方法离散化程度较低，不能够满足道路连接的需要。因此本文先对边缘响应图进行二值化处理，再对该二值图采用八方向模板进行方向计算，最终计算得到灰度范围为0～8的边缘方向图，如图 2c所示。该边缘方向图包含了更高分辨能力的边缘方向信息，可以满足后续道路连接的要求。
式中，I(x, y)、dk分别代表二值化边缘图以及k方向上的滤波窗口函数。
-
粒子滤波(particle filter)是一种基于贝叶斯理论和重要性采样的蒙特卡洛方法，被广泛应用于控制、跟踪及模式识别等领域。与卡尔曼滤波相比，粒子滤波更适合于非高斯和非线性系统。在乘性噪声影响的SAR图像中，对具有低信噪比和强杂波背景的目标检测，粒子滤波可以有效降低虚警，增大发现概率。
-
1) 道路状态描述及预测方程
首先建立道路的状态矢量St={xt, yt, θt}，其中，{xt, yt}代表当前道路片段的中心位置，θt为当前道路片段的走向，对道路上任一点的粒子样本，进行道路候选点状态量的进一步估计，由St-1到St粒子状态转移方程为：
式中，d为设定的预测步长；u1、u2、u3分别为用于修正预测状态的随机噪声。
2) 粒子权重
粒子权重的选择和更新是粒子滤波中最重要的一个步骤，它描述了备选粒子样本属于道路区域的确信度。本文在边缘检测响应、图像强度以及道路片段方向的基础上得到了粒子权重值。
首先，检测得到边缘强度，道路的边缘响应值一般在0.3～0.8之间，边缘响应越高，其属于道路区域的可能性越大。在进行粒子权重计算时，采用高斯分布函数对边缘响应值进行归一化，有：
式中，σ为一个常量，它反映了粒子强度的散度，在本文中采用0.3；ri为第i个粒子的边缘响应值。
其次，道路的灰度值通常处于一个较为稳定低灰度范围内，相邻区域内趋向于保持相同的灰度区间。因此在经过方向性滤波后的SAR图像中，结合阈值分割方法，计算灰度权值为：
式中，ThGrayL、ThGrayH分别为道路的低灰度门限和高灰度门限，低门限用于降低同处于低灰度区内阴影及水体的干扰，高门限用于降低强目标的干扰；gt-1、gti分别表示当前位置处和粒子i在预测位置处的灰度值。由于道路具有上下文联系，因此局部区域内灰度越接近，当前粒子的道路置信度越高。
另外，局部区域内道路片段的方向存在一致性。根据边缘特征提取的结果，对二值化处理后求得方向信息，然后代入式(11)求得粒子的方向权值为：
最后将这3个权重参数融合为粒子权重，得到：
式中，α、β分别为灰度权值和方向权值。在山区的SAR图像中，由于边缘特征检测结果受非道路干扰较多，因此灰度权值的比重应加大，增大灰度权值有利于道路检测。步长d以及方向权值的设定视道路方向变化速度以及道路的分辨率而定。
3) 估计
在这一阶段，粒子滤波的后验概率由每个粒子的粒子权重计算得到，据此，道路状态的估计为：
4) 重采样
粒子在不断传播的过程中会产生粒子退化现象，进而产生粒子贫化，因此需要进行重采样，文中采用的是重要性重采样方法[14]。
-
SNAKE模型[15]又被称为活动轮廓模型，是一个充分利用高层信息自上而下的处理过程，在文中该方法被用于对粒子滤波结果进行精确校正，以得到符合真实道路的道路中心线。其定义为：
式中，$v(s)=(x(s),y(s)),s\in [0,1]$；α、β、γ为控制参数；${{E}_{\operatorname{int}}}(v(s))$为曲线形状的内能量，控制模板的光滑与连续；${{E}_{\text{ext}}}(v(s))$为曲线形状外能量，使曲线轮廓向着目标边缘移动；${{E}_{\text{con}}}$是约束力，根据具体使用环境设定。
将粒子滤波的结果作为SNAKE模型初始控制点，在内能量和外能量的优化约束下进行变换得到精确位置。外能量通常采用图像的边缘信息，如光学图像中采用梯度信息作为外能量。文中采用边缘强度作为外能量函数，有：
式中，${{\gamma }_{\text{edge}}}$为边缘强度。采用贪婪算法寻找总能量的最优解，可以得到精确的道路中心线位置。
1.1. 方向性滤波
1.2. 道路边缘特征检测
1.3. 边缘方向信息获取
1.4. 道路连接
1.4.1. 粒子滤波
1.4.2. 道路追踪模型
1.4.3. SNAKE方法
-
如前文所述，道路在高分辨率SAR图像中表现为具有一定长度和方向的双边缘暗黑条带，在单幅SAR图像中道路损毁发生时通常会表现为道路区域的灰度变化以及与背景区域对比度的损失，因此利用道路的灰度特征和边缘特征进行建模，在确信道路区域内分别提取边缘损毁参数和灰度损毁参数，然后综合二者进行分析并得到道路损毁信息。
-
均值比算子(ROA)通常被用于边缘检测中，基于ROA模型及第二部分所提取出的道路中心线，通过将中心线沿水平方向和垂直方向分别进行多次平移，得到道路损毁检测模型的中心区域及背景区域，从而得到待检测点的两个方向上的局部均值比，利用式(15)计算得到最终的边缘比值参数为：
式中，${{\gamma }_{vij}}{{\gamma }_{hij}}$分别为垂直方向和水平方向上的道路均值比，有：
由参数${{\gamma }_{ij}}$可以得到中心线区域的边缘损毁参数为：
式中，Te为分割道路损毁区域的边缘损毁参数阈值。
-
由于道路在图像中的灰度值较低，因此可以通过对灰度特征的分割将道路分为可能的道路区域和明显非道路区域，利用式(18)获得道路线区域的灰度损毁参数为：
式中，Tg为灰度损毁参数阈值；G(i, j)为点(i, j)邻近位置处一定长度的道路点序列的灰度均值。
最后，将Ig和Ie融合得到最终得到损毁信息提取结果为：
2.1. 边缘损毁参数
2.2. 灰度损毁参数
-
将本文的方法应用于SAR图像中，结果如图 3所示，其中图 3a的区域选择为都汶高速老虎嘴区域，经纬度为31.092 942°N，103.485 025°E，地理环境为山区；图 3b中可知经过方向性滤波后道路区域的对比度和均匀度增加；图 3c中的粒子滤波结果整体符合道路的位置及走向，但在局部范围内仍不准确，存在一定的偏差，图中道路的断裂部分为图像中阴影部分，该区域无法得到正确的粒子滤波结果，因此舍弃该部分，重新选择起始点进行道路追踪；图 3d为SNAKE方法校正后得到的道路中心线，可以看出与图像中道路的真实位置具有更高的符合程度，能满足后续道路损毁信息提取的要求。综上所述，在山区环境中该道路提取方法能有效提取出道路信息。
-
将本文的道路损毁提取方法应用于5.12汶川地震前后的两幅同区域Radarsat2图像中，震后该区域道路损毁严重，同期历史数据表明，地震后都汶高速、国道213线等重要干线的全线损毁率约为80%[16]。在该损毁信息提取实验中，道路损毁阈值确定的范围为Te(0.5～0.7)、Tg(0.0.13～0.015)。为避免人工选择阈值的复杂性，可以采用自动阈值方法进行阈值确定，具体实施方法为根据震前图像的边缘损毁参数和灰度损毁参数的直方分布图，按照一定的恒虚警率(0.05)计算得到两个阈值分割参数。如图 4所示，图 4a为震后SAR图像，对应的震前图像如图 3a所示，两幅图像以图 3a为主图像进行配准；图 4b为震后同时期的航拍光学图像，该图中可以看出该区域由于山体的大面积崩塌和堰塞湖的存在，导致道路被掩埋或损毁的程度极为严重。图 4c、图 4d所示为震后待检测图像和震前基准图像的道路损毁信息提取结果图，其中白色为损毁部分，黑色为未损毁部分，该结果图中的阈值分割方法采用了自动阈值方法，从而计算得到震后与震前的损毁比率分别为91.81%和3.63%。震后待检测图像中的道路损毁信息提取结果与图 4b的分析结论以及同期历史数据相符合，震前基准图像图 3a的道路损毁信息提取结果在虚警条件下也满足震前道路完好的假设。综上，该方法可以有效地区分完好道路与损毁道路，识别出道路的损毁区域。
3.1. 道路信息提取实验
3.2. 道路损毁信息提取实验
-
本文提出了一种基于SAR图像的道路损毁信息提取方法，研究背景主要为山区。该方法包括两部分：1) 基于粒子滤波方法的道路提取，由于山区环境的道路提取具有地形起伏的复杂性，采用方向性滤波增强道路特征，同时增加阴影区域、道路区域以及其他中高灰度均匀区域的可区分性，道路的粒子权重设置以降低阴影以及其他低灰度线形结构的影响，同时道路的上下文联系对道路提取具有重要意义；2) 基于已知道路中心线进行道路损毁模板构建检测损毁。实验结果验证了本文方法能有效地提取道路损毁信息，具有较高的准确度。随着SAR成像技术的发展，在高分辨率和多极化SAR图像中，道路特征将更明显，后续研究中将采用全极化和多波段融合等方法进行道路信息提取。
DownLoad: