线特征约束的建筑物密集匹配边缘全局优化方法

巩丹超，韩轶龙，黄 旭

1. 西安测绘研究所,陕西 西安 710054； 2. 地理信息工程国家重点实验室, 陕西 西安 710054； 3. 武汉大学遥感信息工程学院，湖北 武汉 430079； 4. 武汉市工程科学技术研究院, 湖北 武汉 430019

立体影像密集匹配是在立体影像中，逐像素或者接近逐像素地寻找两张影像之间的同名点，并根据同名光线对对相交的原理，生成三维点云的过程[1]。立体影像密集匹配具有点云密度大、匹配结果稳定、成本低等突出优势，因此可以广泛应用于一些智能3D场景，如：测绘制图、变化监测、国防军事、智慧城市、无人驾驶、虚拟现实、游戏动画、抢险救灾等[2-6]。

立体影像密集匹配常采用固定匹配窗口来计算和比较同名点之间的特征相似性，其基本假设是局部匹配窗口内，所有像素的视差均一致。该假设在大部分场景下适用，但是在建筑物边缘区域，由于存在明显的视差/高程阶跃变化，该假设无法成立。在建筑物边缘区域，固定匹配窗口会带来较大的匹配歧义性，从而降低建筑物边缘的匹配精度，具体表现为：建筑物边缘不规则，建筑物边缘外扩等问题，如图1所示。图1(a)、(c)表示原始影像。图1(b)、(d)表示对应的立体影像密集匹配结果。图1(a)、(c)中红色的线，以及图1(b)、(d)中白色的线是对应的建筑物直线边缘特征。其中，密集匹配算法采用著名的半全局密集匹配算法(semi-global matching,SGM)[7-8]，匹配特征算子采用Census[9]，匹配窗口大小为9×9像素。从图1可以看出，SGM密集匹配算法的匹配结果，在建筑物边缘区域存在明显的外扩问题，从而影响后续建筑物的三维重建精度。除了SGM算法以外，其他基于固定匹配窗口的全局/半全局密集匹配算法，也存在同样的问题[10-13]。除了固定匹配窗口因素以外，匹配常用的视差平滑约束，同样也是边缘外扩的因素之一。

图1 密集匹配算法在建筑物边缘区域的匹配结果Fig.1 The dense matching result in the edge area of the buildings

为了获取高精度的建筑物边缘，目前主要有两种解决方案。第1种方法是在密集匹配过程中加入线特征约束或者面特征约束，削弱建筑物边缘处的匹配歧义性，从而提高建筑物边缘处的匹配精度[14-15]。这类方法往往可以取得比传统SGM更高的匹配精度，但是这些方法的线、面特征约束是根据自身算法的实际情况来定制设计的，不一定适用于其他密集匹配方法，比如，文献[14]采用动态规划思想来优化密集匹配结果，并采用影像线特征来自适应调整动态规划方程中的惩罚项系数，以提高边缘处的密集匹配精度；但是对于不含惩罚项系数的密集匹配方法(如：基于最小生成树的密集匹配方法[16]，基于深度学习的密集匹配方法[17])，是无法使用这类线特征约束。此外，目前测绘单位和地理信息平台存储大量数字表面模型产品(digital surface model,DSM)，第1种方法只能用于生产新的DSM，无法对现有DSM产品进行边缘优化。第2种方法是通过后处理的方式，直接对密集匹配视差图或者DSM中的建筑物边缘进行优化。这类后处理的方法，要么考虑灰度特征分布和实际视差/高程分布之间的相关性，修正所有像素的视差/高程，使其与对应的影像灰度分布近似，从而达到边缘优化的目的；要么从影像中提取的线特征或者面特征，并修正视差/高程边缘，使其与影像线特征或面特征尽量接近。因此，通过后处理来优化边缘的方法可以分为两类：①基于灰度特征分布的边缘锐化算法；②基于线面特征的优化算法。

大多数基于灰度特征分布的边缘锐化算法通过局部窗口内的灰度特征，来改善视差/高程，从而也称为局部边缘锐化算法[18-24]。该类算法假设在局部窗口内，灰度相似的像素，其视差/高程也是一致的。由于建筑物边缘往往也是灰度阶跃的边缘，因此可以采用这类局部边缘锐化算子进行滤波，通过灰度分布来重新调整建筑物边缘，从而获得更加精细的建筑物边缘。但是受限于窗口大小，该类局部锐化算法很难改正建筑物边缘处一些较大的误匹配。为了解决边缘处较大误匹配的问题，需要从更加全局的角度来修正建筑物边缘，也称为全局边缘锐化算法[25-26]。该类算法假设整张影像中，相邻且灰度相似的像素，其视差/高程应该尽量接近，从而将视差/高程优化问题转化为全局能量函数的最优解计算问题。全局边缘锐化算法能够有效优化建筑物边缘，但是由于影像纹理(非边缘区域)也具有灰度阶跃特性，因此在影像纹理区域容易出现“锐化”假象。

基于线面特征的优化算法首先从影像中提取线、面等几何特征，然后改正建筑物边缘的视差/高度，使得边缘特征与影像的几何特征相符。文献[7]通过影像分割算法，提取影像中的面特征，并将每个面特征作为一个视差平面进行优化，从而达到改善边缘的目的。但是，过大的面特征不一定满足视差平面假设，可能导致匹配精度下降。文献[27]进一步采用直线特征提取算子(line segment detector,LSD)来提取建筑物的边缘，并假设边缘特征附近的局部区域满足平面约束，从而实现边缘优化。该方法计算效率高，但是由于只考虑直线附近区域的视差/高程优化，容易导致与直线区域以外的像素之间的视差/高程断裂，并且容易将直线附近区域的地形地貌强行抹平。文献[28]通过平移边缘附近的像素，来改善建筑物边缘，其本质是舍弃直线边缘附近的所有视差/高程，然后以直线边缘为中心，分别利用边缘以外的建筑物视差/高程信息和地面视差/高程信息，向中心边缘进行内插。但是，该方法同样会将建筑物边缘附近的地物抹平。文献[12]计算边缘像素的视差和匹配代价，通过取最小代价的视差，来优化边缘。文献[29]通过超像素分割来确定物体边缘，并通过加权中值滤波来优化物体边缘。

为了解决局部边缘锐化算法难以解决较大误匹配的问题，并解决基于线面特征优化算法强行抹平附近地物的问题，本文提出了一种基于直线特征约束的建筑物边缘全局优化方法，能够在解决较大的建筑物边缘误差的同时，保留建筑物附近的地物。该方法的贡献在于，将建筑物边缘优化问题，转化为一个新的全局能量函数的最优解计算问题，并采用图割算法进行全局最优计算。全局能量函数由代价项和正则项两部分组成。考虑到保留建筑物边缘附近的实际地物，代价项的设计主要根据估计视差/高程与原始视差/高程之间的距离，以及该像素的视差/高程置信度来共同决定的。正则项的设计，则是依据灰度相近的相邻像素，其视差/高程也是相似的原则来设计。本文所提出的建筑物边缘优化方法，能够有效改正建筑物边缘的误匹配问题，解决建筑物边缘外扩的问题，最终重建高精度的建筑物三维边缘。

1 基于线特征的建筑物密集匹配边缘的全局优化方法

本文重点解决密集匹配中，建筑物边缘精度较低的问题。具体流程如图2所示。算法的输入是原始核线影像和视差图(如图2(a)所示)，或者正射影像和对应的DSM，其中，图2(a)中视差图的视差范围是34像素至96像素。首先，采用直线特征提取算子LSD[30]，来检测并提取建筑物的边缘，如图2(b)所示，红色、白色的线分别表示采用LSD算法检测和提取的建筑物边缘。两条边缘线均位于屋顶边缘。然后，分别优化每一条直线特征，将直线特征的优化问题转化为全局能量函数最优解的计算问题。全局能量函数共有代价项和正则项两部分构成，如图2(c)所示。其中，代价项的计算是根据估计视差/高程与原始视差/高程之间的距离Δd，当前像素与建筑物边缘之间的距离s等因素，来共同决定的。Δd越小、s越小，则代价越小，反之则越大；正则项的计算，则是依据灰度相近的相邻像素，其视差/高程也是相似的原则来设计。全局能量函数的具体形式，如1.2节式(1)所示。图2(c)中的3个蓝色的圆形表示边缘附近的3个像素，黑色的连接线表示全局优化中的视差平滑约束，线越粗表示约束越强。灰度相近的像素之间的约束较强，从而在后续全局优化中鼓励他们的视差/高程一致，反之，则灰度差异较大像素之间的约束较小。最后采用图割算法求解整型最优解，并采用基于图像引导滤波的方法获得连续最优解。

图2 本文算法流程Fig.2 The pipeline of the algorithm proposed in this paper

1.1 建筑物边缘检测

本文采用LSD算子来检测和提取影像的直线特征。但是，受地物纹理、阴影等因素的影响，采用LSD算子提取的直线特征，大部分不属于建筑物的边缘。考虑到建筑物边缘的视差/高程阶跃特性，本文采用线缓冲区方法[27]来从LSD直线特征中，进一步提取建筑物的边缘。具体流程为：①首先，采用LSD算子检测影像的直线特征，如图3(b)所示；②以每一条直线为中心，建立缓冲区，统计直线两边的视差/高程情况；③如果直线位于视差/高程的阶跃区域，则作为建筑物的边缘，提取出来，如图3(c)所示。本文算法实质以直线边缘为基本单元，优化直线边缘附近的高程/视差。因此，除了文献[27]提出的基于缓冲区的直线边缘提取方法之外，其他能够提取直线边缘的方法，同样适用于本文算法，如手动提取直线。

图3 建筑物边缘检测结果Fig.3 Building edge detection results

1.2 全局能量函数设计

提取出建筑物边缘的直线特征之后，本文分别对每一条边缘直线特征进行优化。本文算法的核心思想是改变边缘直线附近像素的视差/高程，从而使得改正后的视差边缘/高程边缘与影像的边缘直线尽量接近。因此，首先需要寻找边缘直线附近的像素。本文采用“内/外缓冲区”的思路，来统计边缘直线附近的像素，如图4所示。白色直线表示建筑物的边缘直线。红色方框表示内缓冲区。本文算法通过改变内缓冲区内所有像素的视差/高程，使得视差/高程边缘与影像边缘尽量接近。蓝色方框表示外缓冲区。外缓冲区内像素的视差/高程是不做任何变化的。外缓冲区的作用有两个：①本文采用固定窗口来计算能量函数的代价项，因此，在计算代价项的过程中，会有一部分代价窗口超出内缓冲区的范围。超出的部分可以采用外缓冲区来计算。具体代价计算详见1.2.1小节。②传统的图割算法只能提供整型最优解，因此，在内缓冲区和外缓冲区之间，存在明显的视差/高程断裂。在后续连续最优解计算过程中，将考虑外缓冲区的像素，使得内外缓冲区之间的边界平滑过渡。外缓冲区的半径rout设置为代价窗口的一半。考虑到立体影像密集匹配窗口通常采用9×9窗口[31-32]，本文将内缓冲区半径rin设置为10。

图4 边缘直线的内/外缓冲区Fig.4 Inner/outer buffer zone with straight edges

在统计边缘直线附近的像素之后，本文通过改变内缓冲区中像素的视差/高程，来达到边缘优化的目的。总体来说，本文的边缘优化准则包含两个方面：①优化后的视差/高程，与原始视差/高程尽量接近；②相邻且灰度相近的像素，其视差也是相似的。根据上述两个准则，可以将视差边缘优化问题，转化为全局能量函数的最优解计算问题，如式(1)所示

minE(D)=∑p∈BinSc(p)·Cost(p,d)+

∑p,q∈NP·wc(p,q)·we(p,q)·

min(|dp-dq|,σ)

(1)

式中，D表示内缓冲区中所有像素的估计视差/高程的集合；E(D)表示全局能量函数的目标函数值；Bin表示内缓冲区所有像素的集合；p表示内缓冲区中的任意一个像素；Cost(p,d)表示像素p对应于视差/高程d的代价，是根据估计视差/高程与原始视差/高程之间的距离来计算，详见1.2.1小节介绍；Sc(p)表示像素p的置信度，表示像素p的代价计算的可靠性，用于降低误匹配点的代价权重，具体详见1.2.1小节介绍；N表示所有邻域像素的集合；q表示p的邻域像素；dp、dq表示像素p、q的估计视差/高程；σ表示截断值，用于定义视差阶跃条件的阈值；P表示惩罚项系数；wc(p,q)表示根据像素p、q的灰度差，计算出来的权重，一般灰度差越小，权越大；we(p,q)表示根据像素p、q的位置，定义的权重。如果p、q位于直线边缘的同一侧，则权we为1，否则，权we为0。we用于禁止直线两边的优化结果相互影响，从而能够获得更加锐利的边缘。本文将能量函数式(1)中的第1项称为代价项，将第2项称为正则项。代价项和正则项的具体介绍，详见1.2.1小节和1.2.2小节。

1.2.1 代价项计算

大多数算法[7,27-28]采用平面模型来改变边缘直线附近的视差/高程，因此，会强行删掉边缘直线附近的实际地物。为了解决这个问题，本文将估计视差/高程与原始视差/高程之间的距离，作为代价计算的一个依据。当距离越小，则代价越小；反之，距离越大，则代价越大，从而尽可能地保留直线附近的实际地物，如式(2)所示

(2)

但是，原始视差/高程同样存在一定的误匹配问题，单纯采用上述距离测度计算的代价，有可能在优化中将误匹配点也保留下来了。因此，本文定义了视差/高程置信度准则，用于衡量内缓冲区中每个像素是正确匹配点的概率。置信度越高则表示正确匹配点的概率越大，则在后续全局优化中的权重越大；反之，则权重越低。本文定义了两种置信度准则：第1种是像素距离边缘直线的距离，考虑到匹配边缘的外扩问题，距离越近则是正确匹配点的概率越小，距离越远则是正确匹配点的概率越大；第2种是像素与周围灰度相近像素之间的视差/高程差异，差异越小则是正确匹配点的概率越大，差异越大则是正确匹配点的概率越小。因此，综上两个准则，置信度的计算如式(3)所示

(3)

1.2.2 正则项计算

能量函数的正则项用于对缓冲区内的所有像素引入视差/高程平滑约束。相邻像素的灰度越接近，则视差/高程平滑约束越强，从而在能够改正边缘直线附近的误匹配点的同时，尽可能多地保留缓冲区内的细节信息。总的来说，本文定义的正则项与3个要素有关：①相邻像素的灰度差；②相邻像素是否位于直线的同一侧；③相邻像素的视差/高程差异。

相邻像素灰度越接近，则正则项的视差/高程平滑约束越强。因此，本文采用高斯核函数来定义相邻像素之间的灰度接近程度，如式(4)所示

(4)

式中，Ip、Iq分别表示像素p、q的灰度；σc表示与灰度相关的高斯核函数平滑因子，本文取10。当像素p、q之间的灰度越接近，灰度权wc越接近1。

考虑到有些情况下，边缘直线两侧的灰度差异较小。在这种情况下，灰度权wc会较大，从而导致边缘直线两侧的视差/高程平滑过渡，降低边缘直线的锐利程度。因此，为了保证优化后的视差/高程边缘，与影像边缘直线特征尽量接近，本文禁止直线两侧像素之间的正则平滑约束，如图5所示。其中，白色直线表示建筑物直线边缘；蓝色圆圈代表缓冲区内的像素；黄色线表示相邻像素之间的正则项，线的粗细表示正则项平滑约束的大小；红色线框表示内缓冲区。在本文算法中，只有位于同一侧的相邻像素才有正则项的平滑约束。

如果相邻像素属于同一侧缓冲区，则需要引入正则项平滑约束，将正则项中的同侧权值we(p,q)定义为1；反之，则定义为0，如式(5)所示

(5)

1.3 连续最优解计算

本文采用图割算法[34]，获取全局能量函数(式(1))的最优解。但是，传统的图割算法只能获得整型最优解，从而导致缓冲区内的像素和缓冲区外的像素出现轻微的视差/高程断裂，如图6(a)所示。其中，白色的线表示影像上的边缘直线特征；红色的线是剖分图轨迹，对应于右侧的折线剖分图。为了获得连续解，需要综合考虑外缓冲区像素的视差/高程，对内缓冲区像素进行平滑滤波。为了保证滤波后的视差/高程边缘的几何精度不受影响，本文采用保边滤波算子(如基于图像引导的滤波[22])，来对缓冲区内的像素进行平滑滤波。考虑到一些房子的屋顶与地面的灰度差异不是特别明显，即使采用保边滤波算子，也无法在滤波过程中保持建筑物边缘的尖锐性。因此，本文以边缘直线特征为界限，分别对边缘直线两侧进行滤波。当有滤波窗口中的像素超出边缘直线时，超出部分的像素不参与滤波计算，如图6(a)所示。图6(a)中的绿色方格表示当前需要滤波的像素；蓝色方格表示周围的8邻域像素。其中，最右侧的3个方格超出边缘直线部分，则这3个方格不参与滤波计算。

图6 连续最优解计算示意Fig.6 Continuous optimal solution calculation

基于图像引导的滤波[22]是一种高效的保边平滑算子，其核心思想是将局部窗口内的像素视差/高程进行加权平均。窗口内每个像素的权值大小取决于局部视差和灰度之间的线性模型。由于视差/高程边缘往往也是灰度边缘，因此基于图像引导的滤波具有很好的保边性质。本文采用基于图像引导的滤波算子，来进一步改善内缓冲区中的整型视差/高程。在滤波过程中，外缓冲区中的像素同样参与计算，从而能够为内缓冲区视差/高程的优化提供边界约束。基于图像引导的滤波算子的定义如式(6)所示

(6)

2 试验与结果分析

本文分别采用航拍立体像对数据集和卫星立体像对数据集，从定量和定性两个角度来分析本文算法的正确性和有效性。其中，航拍立体像对数据集包括International Society for Photogram-metry and Remote Sensing(ISPRS)提供的Vaihingen航拍数据集和Toronto航拍数据集；卫星立体像对数据集包括约翰斯·霍普金斯大学和IAPRA提供的杰克逊维尔地区的WorldView-3卫星影像数据集和对应的激光点云(LiDAR)DSM。

2.1 航拍立体像对数据集

该部分试验主要测试视差图中建筑物边缘的优化效果。其中，Toronto数据集由Microsoft Vexcel’s UltraCam-D(UCD)相机拍摄，影像大小为2000×2000像素，地面分辨率为15 cm；Vaihingen数据集由Intergraph/ZI DMC相机拍摄，影像大小为2000×2000像素，地面分辨率为8 cm。两个数据集共由9个立体像对组成。每个立体像对均重采样成核线影像，并采用SGM半全局密集匹配算法[8]计算对应的视差图，并采用左右一致性检测剔除误匹配点，然后再进行视差内插。具体测试数据集和对应的原始视差图如图7所示。图7的上半部分表示原始核线影像，下半部分表示SGM密集匹配结果。所有立体像对均包含较密集的房屋，而传统的密集匹配算法在这些建筑物的边缘区域的匹配精度较差。

图7 航拍数据集与SGM密集匹配视差Fig.7 Aerial data and dense matching disparity using SGM

采用本文算法，最常用的双边滤波算法和加权中值滤波算法，分别对SGM密集匹配结果进行边缘优化，并对各个算法的边缘优化结果进行性能评估、对比和分析。其中，双边滤波算法[33]和加权中值滤波算法[24]均假设距离相近且灰度相近的像素，其视差应该是一致的。由于建筑物边缘往往也是灰度阶跃的影像边缘，因此，两种方法均能对建筑物边缘的匹配结果进行优化，并且已经广泛应用于一些立体影像密集匹配算法中[13,35]。为了与本文算法的参数尽量接近(缓冲区半径为10)，双边滤波和加权中值滤波两种方法均采用11×11的滤波窗口。由于本文算法仅仅优化建筑物的边缘，因此，本文仅在提取的建筑物边缘区域(1.1节介绍的方法)分别评估3种方法的边缘优化结果。由于建筑物边缘的误匹配，绝大多数是由于建筑物外扩导致的，因此，本文根据直线缓冲区对应地面那一侧内的误匹配点数目，来评估各种算法的优化性能，如图8所示。其中，白色直线表示建筑物边缘；黄色像素表示误匹配点；蓝色像素表示正确匹配点；红色方框表示直线的缓冲区。本文根据图8中的黄色像素和蓝色像素数目，来综合评价边缘优化算法的性能。

perbad=|Serror|/|Serror+Scorrect|

(7)

式中，|·|表示集合中所有元素的数目；perbad表示误匹配像素百分比。perbad越小，表示边缘优化效果越好

σlow=

(8)

采用本文优化算法，双边滤波和加权中值滤波，分别对图7中9个立体像对的SGM密集匹配结果进行边缘优化，并采用式(7)和式(8)两个测度，分别对原始匹配结果、本文优化算法结果、双边滤波结果和加权中值滤波结果进行性能评定。评定结果如图9(a)和(b)所示。

图9 航拍数据集不同方法的性能对比结果Fig.9 Performance comparison results using different edge-refine methods on aerial dataset

由图9可以看出，在两种性能评定测度中，本文算法、双边滤波方法和加权中值方法均能够提高建筑物边缘的匹配效果。在误差像素百分比测度perbad中，原始视差图的误匹配点百分比为25.53%，加权中值方法平均能够将建筑物边缘的误匹配像素百分比减少到20.04%，双边滤波方法平均能够将误匹配像素百分比减少到14.21%，而本文算法平均能将误匹配像素百分比减少到9.12%。在方差测度σlow中，原始视差图建筑物边缘区域的方差平均为5.86像素，加权中值方法平均能够将方差减少到5.48像素，双边滤波方法平均能够将方差减少到4.34像素，而本文算法平均能将方差减少到3.63像素。总体来说，本文算法在两种测度中的整体表现均是最好的。具体到每个立体像对的边缘优化结果评定中，本文算法在9个像对中的表现，同样也都是最好的。这是由于双边滤波方法和加权中值方法，均是采用局部窗口，用于优化的信息较少；而本文方法用到了边缘直线缓冲区内的所有信息，并进行全局优化，从而优化效果更好。为了更全面地评估本文算法的性能，本文以第1对核线立体影像密集匹配的边缘优化结果为例，列出一些优化前后的建筑物边缘匹配结果的对比，以及对应的建筑物剖面图，如图10所示。

图10 本文算法的边缘优化效果Fig.10 Building edge refined results using our proposed algorithm

在图(a)展示了原始影像、原始视差图和本文算法优化后视差图的缩略图结果。由于建筑物边缘的误匹配像素往往只有几个像素的误差，因此，在缩略图中看不出太明显的变化。为了进一步展示本文算法的效果，本文挑选了3个红框所示的区域，并放大显示，如图10中(b)图—(d)图所示。在图10(b)图—(d)图的所有原始影像中，红色直线表示提取的影像直线边缘；所有原始视差图中，白色直线表示提取的影像直线边缘，红色直线表示建筑物剖面图轨迹，对应最右侧的剖面图，用于检查优化前后的建筑物边缘变化。剖面图中，红色的线表示原始视差图边缘；绿色的线表示经过本文算法修正后的边缘；蓝色的线表示影像直线边缘。影像直线边缘与地面之间的视差差异，主要是根据影像边缘直线的位置，及其缓冲区两侧平均视差所共同决定的。原始SGM密集匹配结果，其建筑物边缘距离影像直线边缘会外扩若干个像素，从而降低建筑物边缘的三维重建精度。本文算法能够高精度地将建筑物边缘修正到影像直线边缘附近0～1个像素。在图10(a)中，本文算法修正后的边缘，与影像直线边缘完全重合，而在图10(b)、(c)中，本文算法修正后的边缘，仅仅与影像直线边缘差1个像素，从而充分说明本文算法的正确性和有效性。

2.2 卫星立体像对数据集

该部分试验主要测试DSM中建筑物边缘的优化效果。其中，WorldView-3卫星数据集总共包含4个立体像对，空间分辨率为0.3 m，像对重叠区域约为0.5 km2。激光点云DSM的空间分辨率为0.5 m，范围约为0.1 km2。首先，将每个立体像对重采样成核线影像，采用SGM半全局密集匹配算法[8]计算对应的视差图，并采用左右一致性检测剔除误匹配点；然后根据视差图生产对应的DSM，并进行DSM内插。为了便于利用激光点云DSM来评估精度，每个卫星立体像对DSM和正射影像均切割成与激光点云DSM同样的范围，并将两者的空间分辨率重采样成0.5 m。具体测试数据集、对应的原始DSM和激光点云DSM如图11所示。图11的上半部分表示卫星正射影像，下半部分表示立体像对DSM和激光点云DSM。

采用本文算法和基于缓冲区平面模型的建筑物边缘修正方法(plane based boundary refinement,PBR)[27]，分别对图11中的4组原始DSM进行边缘优化，并以激光点云DSM作为高程真值，对各个算法的边缘优化结果进行性能评估、对比和分析。其中，基于缓冲区平面模型的建筑物边缘修正方法(PBR)[27]假设缓冲区内的高程/视差满足平面约束，通过计算缓冲区内的平面方程，来优化建筑物的边缘。由于两种算法均仅仅优化建筑物的边缘，因此，本文仅仅在提取的建筑物边缘区域(1.1节介绍的方法)评估两种方法的边缘优化精度。为了客观评价两种算法的优化效果，以激光点云DSM作为真值，定量分析优化后DSM与原始DSM相比，精度的提升情况。总的来说，本文通过如下指标来评价各个算法的性能：①平均误差变化，建筑物边缘区域的平均误差在优化前后的变化(即优化后的平均误差减去原始平均误差)；②小于2 m误差百分比变化(bad-2)，统计误差小于2 m的像素所占百分比在优化前后的变化(即优化后的百分比减去原始的百分比)；③小于3 m误差百分比变化(bad-3)，统计误差小于3 m的像素所占百分比在优化前后的变化。采用这3种精度评价准则，分别对PBR和本文算法进行性能评定，结果如图12所示。在图12中，平均误差变化越小表示优化后的精度提升越大；而bad-2，bad-3越大，表示优化后正确点的百分比提升越大。

图11 卫星数据集与对应的数字表面模型Fig.11 Satellite dataset and the corresponding DSM

由图12可以看出，本文算法的优化效果明显优于PBR算法(除了(c)图中的像对2)。在平均误差变化测度中，PBR算法在部分像对的平均误差变化大于0，说明PBR算法反而降低了原始DSM的精度，这是因为建筑物边缘附近经常存在其他地物(如树木)，不符合平面约束条件，因此基于平面约束的PBR算法在这种情况下会降低DSM精度。而本文算法采用了更加灵活的约束，依靠相邻且灰度相似的像素来优化建筑物边缘，从而能保留建筑物附近的不相似地物，最终提高DSM的精度。在大部分像对中，PBR算法的bad-2和bad-3两个测度小于0，说明PBR算法反而减少了建筑物边缘附近的正确点数目。而本文算法在大部分像对中的bad-2测度和bad-3测度大于0，说明经过本文算法优化，建筑物边缘附近的正确点数目得到提升。因此，综合高程误差和正确点百分比两个方面，均证明本文算法能够提高卫星DSM中的建筑物边缘精度。为了更全面地评估本文算法的性能，本文列出一些优化前后的建筑物边缘对比，如图13所示。

图12 卫星数据集的性能对比结果Fig.12 Performance comparison results on satellite dataset

在图13中，正射影像中的红色直线表示提取的边缘直线，与DSM中的白色直线是一一对应的关系。原始DSM中的红色像素表示因为误匹配而导致外扩的边缘，而被红色曲线所包围的像素(如图(c)所示)，表示因为遮挡区域高程内插而导致的边缘外扩问题。PBR算法通过强行将边缘附近的区域拟合为平面，获得与正射影像直线特征完全吻合的建筑物边缘，但是这类过强的约束，同时也会破坏建筑物附近的地形地貌，反而降低了DSM的精度，如图13(b)和(c)所示。本文算法能够在优化建筑物边缘的同时，有效保留建筑物附近的不相似地物，如图13(b)所示，并对遮挡区域的建筑物边缘，有着明显的改善效果(如图13(c)中的红色曲线区域所示)。但是，本文算法的效果依赖于图像的质量，当边缘误匹配像素的灰度信息，与周围正确匹配像素的灰度信息存在较大差异时，本文算法无法对这类边缘进行优化，如图13(a)中下方的建筑物边缘所示。此外，尽管采用图割算法进行全局优化，本文算法的时间复杂度较低。以WorldView-3卫星数据集为例，在单核CPU @2.60 GHz的配置下，采用本文算法处理一张512×512像素的影像，平均处理时间为0.699 s，平均每条边缘的处理时间为0.012 s。这是因为本文算法仅仅处理建筑物边缘附近的像素，而不是处理整张影像，因此，时间复杂度较低。

图13中的边缘优化对比是以直线段为基本单元。为了更加全面地评估本文算法的边缘优化效果，本文以图11中的立体像对4作为测试对象，测试整栋建筑物的直线边缘优化效果。测试中的立体像对DSM和正射影像采用切割前的产品，范围约为0.5 km2(2221×2223像素)，采用文献[27]的方法对建筑物进行直线边缘检测，部分没检测出来的直线，采用人工方法提取，从而形成完整的建筑物轮廓，总计提取出46栋建筑物的边缘。然后，采用本文算法进行建筑物边缘的全局优化，优化结果如图14所示。图14中，图(a)表示优化前后的全局DSM；图14(b)—(d)表示局部若干栋建筑物的边缘在优化前/后的对比。图14(b)—(d)的具体地理位置，如图(a)中的红框所示。由于本文算法只优化建筑物边缘，因此从整体上看，优化前后的两组全局DSM区别不大，如图(a)所示。但是，在局部范围内，建筑物边缘的优化效果较为明显，如图14(b)—(d)所示。其中，优化前建筑物边缘有一定程度的外扩。在建筑物和地面之间存在平滑过渡区域，从而导致边缘模糊化。经过本文算法优化后，建筑物在DSM中的边缘具有较明显的直线特征，从而提高了建筑物边缘重建的精度。在单核CPU @2.60 GHz的配置下，采用本文算法处理该区域的时间为6.23 s，平均每栋建筑物的边缘优化时间为0.14 s。本文算法以直线边缘为基本单元进行优化，即一个个的直线进行优化，而非DSM/视差图的整体处理。因此，本文优化方法的时间效率取决于建筑物直线边缘的提取数目。考虑到建筑物边缘及其附近的像素数目，在整张DSM中所占百分比较小，因此尽管采用了图割解法，本文优化算法的时间效率仍然较高，每栋建筑物的优化效率可以达到亚秒级。考虑到每条直线边缘单独处理，因此本文优化算法非常容易用于并行框架，从而能够进一步提高处理效率。

图13 两种算法的边缘优化效果对比Fig.13 The boundary refinement comparisons of the two algorithms

图14 整栋建筑物优化前后的边缘对比效果Fig.14 The comparisons of building boundaries before/after the optimization

3 结 论

本文提出了一种基于直线特征的建筑物边缘精化方法，用于提高建筑物边缘的匹配精度或者高程精度。采用航拍影像数据集Vaihingen和Toronto，以及WorldView-3卫星数据集来验证本文算法的正确性和有效性。试验结果表明：本文算法明显优于常用的边缘锐化算子(双边滤波和加权中值)和一种新的基于缓冲区平面的边缘优化方法，能够有效改善建筑物边缘的精度，减少建筑物边缘的误匹配像素。因此，本文算法在一些建筑物三维重建场景(如：虚拟现实、智慧城市、城市管理等)，均能得到应用。在未来，计划结合深度学习技术，研究更精确的影像边缘直线提取算法，从而进一步提升建筑物边缘精化的效果。

致谢：感谢ISPRS提供的Vaihingen和Toronto航拍数据集；感谢约翰霍普金斯大学应用物理实验室，IARPA以及IEEE GRSS图像分析和数据融合技术委员会提供的融合竞赛数据。