道路是国家经济和军事的动脉。地震灾害发生后，快速评估灾区的道路损毁状况对于应急救援工作具有重大意义。合成孔径雷达(synthetic aperture radar,SAR)能够克服天气和光照条件的影响，更适合在恶劣天气中用于灾情监测。但是，其特有的相干斑噪声增加了道路提取的困难，同时也为损毁信息带来了干扰。

根据在道路损毁信息提取过程中所用信息源，可将道路损毁提取方法分为三类^[1]：(1) 基于灾后单时相影像的方法；(2) 基于灾前-灾后多时相影像的方法^[2]；(3) 基于灾前先验信息-灾后影像多源信息融合的方法^{[3, 4]}。与基于灾后单时相影像检测方法相比，通过灾前-灾后不同时相数据利用变化检测方法进行对比分析的研究较多，且识别精度相对较高。道路对象，特别是主要道路，通常比较稳定，自然条件下变化较小，根据灾前GIS道路网数据判断灾后发生变化的道路具有较高的可靠性^[5]，且能克服灾前影像数据不完备的问题。目前利用变化检测进行损毁信息提取的方法直接将变化结果作为损毁结果，没有对变化结果作进一步分析。实际上道路背景本身的复杂性，如车辆、树木的干扰及SAR影像中的阴影、叠掩、相干斑噪声等，容易造成实际未变化道路被误检测为损毁信息，形成虚警。因此，在提取变化信息后，有必要进一步排除干扰，剔除误检测的疑似损毁区域。

基于以上分析，本文提出了一种基于GIS与贝叶斯网络的高分辨率SAR影像道路损毁提取方法。其特点在于：(1)在GIS数据引导下利用水平集分割与改进D1线检测融合的方法提取疑似道路损毁区。(2)将多损毁证据引入贝叶斯概率模型对疑似区域进行进一步分析并剔除干扰信息，最终提取出实际道路损毁区。

道路在高分辨率SAR影像上通常呈暗条带状区域，但噪声、道路与周围建筑物粘连、树木遮挡等使得部分道路双边缘变成单边缘，从而利用单一线检测或面分割的方法容易导致道路检测不完整^[6]。因此，本文提出在GIS数据引导下以水平集分割与改进D1检测融合的方法检测出道路，进而提取疑似损毁区。

水平集是一种将曲线(或曲面)隐藏在更高一维连续曲面的零水平集中隐式地完成曲线演化的方法，大量研究已证明利用水平集理论实现高分辨率遥感影像目标分割的可行性。但传统的水平集分割容易受到各种干扰信息的影响，因此越来越多的研究考虑加入先验知识，如利用目标的先验形状来约束水平集分割，获得更加精确的结果^[7]。

加入先验形状约束后的能量函数的计算方法为：

式中，φ表示水平集曲线；r、g分别为图像中目标和背景两个区域；α为控制权系数；E_total 为加入先验形状约束后的能量函数；E_C是水平集方法的能量项； E_prior是先验形状能量项。针对最简单的情况，即先验形状不进行缩放、平移、旋转等变化时，先验能量项可采取如下的形式：

式中，φ₀表示已知形状的水平集函数；Ω表示图像域；x、y则为曲线的坐标变量。在使用水平集方法进行图像分割时，需要式 (1) 的能量泛函进行最小化，一般使用变分法，从而式(1)可重写为：

对上述泛函求变分，可得到u₊和 u_-的形式。通过上述方程，对于水平集函数φ，其欧拉-拉格朗日(E-L)方程为：

式中,H(φ)为Heaviside函数；δ(φ)为数值化的狄拉克函数，即δ(φ)=H′(φ)； ψ表示先验形状的水平集函数；f为图像灰度值；u₊和u_-为两个随着水平集淡化不断更新的标量，分别代表了ψ≥0和ψ＜0的不同区域的灰度均值。加入先验约束的水平集方法可以很好地解决由于底层信息缺失造成的图像分割错误。

为了充分利用GIS数据作为先验知识，考虑对GIS道路数据以一定宽度构建长条型缓冲区域，获取闭合的初始道路边界形状作为初始水平集函数。缓冲区宽度d可以任意设置为一个小于实际道路宽度的值，通常d的取值为3~5像素。最终通过水平集分割出的道路目标宽度记为实际道路宽度W。将分割后的道路结果与矢量道路缓冲区进行叠加分析，分割道路上出现断裂而矢量道路完整的区域即为疑似损毁道路。

在低分辨率SAR影像上道路呈细长暗线，D1算子^[8]可以很好地提取出这类线特征。但在高分辨率SAR影像上，道路宽度所占的像素远大于D1算子中央区域的1~4像素，因此利用D1算子在高分辨率影像上检测道路不太合适。同时，由于D1算子是一种线检测算子，直接用于影像会提取出很多非道路信息。因此，借助于GIS数据的位置和方向信息，本文提出了一种基于D1算子的改进道路检测算子。

图 1 D1 算子和改进的D1 算子检测窗口 Fig. 1 Detection Windows of D1 Detector and the Improved D1 Detector

图选项

改进的道路检测算子与D1相似，都属于比率检测算子。设其中央区域的灰度均值为a₂,左右两侧区域的灰度均值为a₁、a₃。对于某一段道路，其线检测响应值(line detector response,LDR)的计算方法为：

在矢量数据与影像粗配准的基础上，本文在宽度为W的矢量数据缓冲区内对影像进行道路检测。如图 2所示，检测窗口的中央区域宽度为W，在每段缓冲区内，检测窗口从矢量出发，沿着垂于矢量方向移动，每次移动距离为W/4,记录所有大于检测阈值的检测窗口中心线为道路基元。这样最坏的情况是检测出来的道路线段基元和实际道路中心线的最大距离为W/8，即道路中心线定位精度为1/8个道路宽度。

图 2 改进的D1检测窗口在缓冲区内的检测方式 Fig. 2 Detection Mode of the Improved D1 Detector in GIS Buffer

图选项

本文提出的改进算子与D1算子的主要区别在于：

1) 本文算子以矢量路段为基元，判断线段的响应值，而D1算子是判断某点的响应值。

2) 中央区域宽度固定，其宽度来源于影像分割得到的实际道路宽度，这样避免了其他与道路目标相似的地物干扰。

3) 方向由8个可移动方向变成1个方向。在矢量数据的引导下，线检测的方向固定为该位置对应的道路矢量数据的方向。

经过道路基元的编组、合并和筛选之后，得到最终道路线段集合。如果SAR影像上道路在某处发生损毁，该损毁区必定将整段道路划分成为两个线段组集合，且损毁区域不存在线段组。将线段组集合与矢量数据结合即可找出道路损毁的位置。

图 3 基于改进D1算子的疑似损毁道路判断 Fig. 3 Judgment of Suspected Damaged Road Based on the Improved D1 Detector

图选项

为了增加损毁区检测的可信度、减少漏检，本文将基于水平集分割与基于改进D1检测的结果进行融合。对于由水平集分割提取的疑似损毁路段与由改进D1方法提取的疑似损毁路段，均认为其为最终的疑似损毁区。融合后，疑似损毁区的虚警率虽高，但能够保证较低的漏检率，同时虚警的疑似损毁区在后续过程中会进一步判断剔除。

贝叶斯网络模型是一种不定性因果关联模型，利用条件概率表达各个信息要素之间的相关关系，能在有限、不完整、不确定的信息条件下进行学习和推理^[9]。因此,本文采用贝叶斯网络模型对检测的断裂区作进一步分析判断。

根据实际情况，将造成道路疑似损毁的原因简要分为以下几类情况：

1) 地震引发道路两旁建筑物倒塌，导致道路堵塞。这类损毁与房屋到道路的距离(房屋轮廓数据A)相关。

2) 地震或其他原因造成的滑坡、泥石流等阻断道路。这类损毁与该区域的地质条件(滑坡隐患点数据B及泥石流隐患点数据C)相关。

3) 暴雨或者地震引起的堰塞湖等原因造成道路被洪水淹没。这主要与洪水淹没信息(洪水淹没数据D)相关。

4) 高大建筑物和地形起伏引起的叠掩与阴影。这主要跟DSM/DEM信息(DSM/DEM数据E)相关。

5) 其他干扰。

损毁区的实际属性依赖于以上5件证据，而不同原因造成的道路断裂在SAR影像上的观测值是不同的，其观测值以灰度和纹理来加以表示，本文采用能量、熵、惯性矩、相关性和局部平稳性等5种常用纹理统计量进行实验。因而可构造如下道路损毁的贝叶斯网络模型：

图 4 基于多证据判断的贝叶斯网络模型 Fig. 4 Bayes Network Based on Multi-evidence

图选项

5个先验证据信息变量A、B、C、D、E是可疑断裂区实际属性X的条件，可疑断裂区实际属性X是观测值F、G的条件。那么，在各种证据辅助下，某疑似断裂区X的实际后验概率P(X|A,B,C,D,E,F,G)可计算为：

其中，P(X|A,B,C,D,E)表示该断裂区属于某种损毁原因的先验概率；P(F|X)，P(G|X)表示当前断裂情况下的灰度和纹理概率分布。

为了验证本文提出方法的有效性，以汶川县灾后COSMO影像和开放街道地图(Open Street Map)矢量道路为数据源提取由地震引发的道路阻塞，其中影像获取日期为2012年5月13日，分辨率为1 m。矢量数据与影像通过手动选择3~4个同名点进行配准。

图 5显示了利用水平集分割与改进D1检测融合方法得到的疑似损毁区检测结果，共有26处疑似损毁路段。其中，由于地形的影响，SAR影像上存在着明显的叠掩现象，且桥梁、交叉路口等干扰信息被检测为断裂区域。右侧显示了其中三个疑似断裂区Ⅰ、Ⅱ、Ⅲ的放大效果，分别对应着光学影像和SAR影像，依次是由叠掩、滑坡、道路经过水体造成的。

图 5 基于水平集分割与改进D1检测融合的疑似损毁区检测 Fig. 5 Uspected Road Damage Detection Based on the Fusion of Level-Set Segmentation and the Improved D1 Line Detection

图选项

图 6显示了基于贝叶斯网络模型验证的实道路际损毁提取。由于辅助数据有限，本文只借用了OpenStreetMap矢量数据的水域信息和房屋轮廓信息、DEM数据和滑坡隐患点数据对断裂区进行判断。对于断裂区可能存在的5种状态，在本文实验中简化到了房屋倒塌导致的损毁、滑坡导致的损毁、叠掩与阴影导致的虚警和其他干扰信息导致的虚警4种。本文的训练样本来自5\512地震北川地区的Radarsat2影像、都江堰地区的Cosmo影像以及茂县地区的TerraSAR影像，其中由房屋倒塌、滑坡、叠掩与阴影及其他干扰产生的道路损毁正负样本各取5例。假定A表示房屋轮廓数据，其取值为{0,1}，1表示断裂区两旁不存在房屋，否则为0；B表示滑坡隐患点数据，取值为{1,0},1表示断裂区未处于滑坡隐患点范围内，否则为0；C表示DEM数据，取值为{1,0}，1表示断裂区未处于叠掩与阴影范围内，否则为0；D表示水域轮廓数据，取值为{1,0}，1表示断裂区处于水体上，否则为0。X表示断裂区域的实际属性，被简化为4个取值: 1表示房屋倒塌，2表示滑坡，3表示叠掩与阴影，4表示其他干扰。断裂区域实际属性之间的先验概率定义方式为:

上式共包含了X、A、B、C、D的64(4×2×2×2×2)种取值情况。P(G|X)、P(H|X)的概率密度分布通过样本学习获得。最终，基于各种辅助证据和观测值，疑似损毁道路的后验概率为:

不论X取何值其分母都是一致的，因此判断X的最大后验概率只需要计算分子即可，取后验概率概率最大的情况作为X的实际属性。

本次实验最终有5处被判定为损毁路段。与图 5比较可知，经过多证据和断裂区观测值判断后，叠掩引起的道路损毁、道路经过水体引起的虚警和一些其他干扰信息都被剔除。

本文以国家减灾委员会发布的汶川地震317和213国道汶川段堵塞情况遥感监测图作为参考，检验本文方法道路损毁提取的准确性，结果如表 1所示。

综合图 6与图 7可以看出：

图 6 基于贝叶斯网络模型验证的实际道路损毁提取 Fig. 6 Real Road Damage Detection Based on the Bayes Network

图选项

图 7 四川汶川地震317和213国道汶川段堵塞情况遥感监测图 Fig. 7 Remote Sensing Monitoring Image for 317 and 213 National Road in Wenchuan Earthquake

图选项

1) 利用本文算法从高分辨率SAR影像上检测出的5段滑坡堵塞区域中，B、D、E三处区域在6图中都有对应区域，属于真实的滑坡引起的道路堵塞；

2) 图 6检测出A区域在图 7的范围之外，因而未被标识，事实上通过对应区域的光学影像对照可知，A区域确实发生了滑坡。

3) C区域未发生滑坡引起的堵塞，由于其观测值和滑坡相近，被误检测为了损毁区。

4) 除了A、B、C、D、E区域之外，图 7所示的范围内还有4处滑坡区域未被检测出来，其中F、G、H三处是由于在SAR影像的叠掩区域而被判断为叠掩，1处是由于没有对应的道路矢量数据而未被检测出来。

从表 1中可以看出，在3个正确提取的损毁区中，道路损毁检测的正确率都达到了84%以上。可以看出，基于多证据辅助和断裂区域观测值构建的贝叶斯网络模型能够较好地剔除干扰信息。其中，漏检的主要原因在于损毁区域处于叠掩与阴影区域或者对应的区域缺乏道路矢量数据而无法判断，虚警主要是由于干扰信息和滑坡的观测值接近，因此需要更多训练样本和辅助数据更加精确地定义疑似损毁区各种情况下观测值的概率密度分布。

本文利用GIS数据作为引导，提出了基于水平集分割与改进D1检测融合的方法从高分辨率SAR影像中提取疑似道路损毁区，并利用贝叶斯网络模型对疑似区进行进一步判断。实验结果证明，该方法能够快速、准确地对损毁信息进行提取。然而，由于损毁信息的复杂性，本文在构建贝叶斯网络模型时需要进一步考虑到更多的影响因素及训练样本，同时由于单极化SAR影像在表现损毁信息上的局限性，如何有效描述灾害损毁在SAR影像上的特征并有效结合利用其他信息，如SAR极化信息，准确地提取损毁区域是今后进一步研究的问题。