航空摄影测量技术因非接触、成本低等优点被广泛应用，可以从二维航摄影像自动得到三维空间的几何信息、辐射信息和语义信息，进而生产4D产品。影像密集匹配是在获得影像间的相对位置关系之后为重叠区域内的每个像素确定同名点的过程，在计算机视觉和数字摄影测量领域中相关研究较多。密集匹配点云因其具有颜色和纹理信息，已成为摄影测量数字表面模型生成、正射影像纠正以及三维重建工作的重要数据来源^{[1, 2]}。

摄影测量与计算机视觉都是研究物体与影像之间关系的，随着计算机的发展和进步，摄影测量技术的发展借鉴了许多计算机视觉方面的研究成果^[3]。密集匹配算法一般主要包括4个步骤^[4]：匹配代价计算、代价聚合、视差计算与优化、视差后处理。依据所采用的数学模型及最优化理论可将影像密集匹配分为局部密集匹配和全局密集匹配。局部密集匹配算法主要是基于窗口进行代价聚合，并且在视差优化方面直接采用WTA（winner takes all）策略获得视差图，虽然随着自适应支持权重窗口算法^[5]以及后续的各种改进算法的提出，匹配精度越来越高，但其在纹理缺乏区域、深度不连续处匹配效果仍不佳；全局密集匹配方法给出了显式平滑性假设，主要是优化由整个图像的数据项和平滑项组成的全局能量方程，常用动态规划^[6]、图割^[7]和置信度传播^[8]等方法求解能量方程。相对于局部匹配算法，全局匹配算法的准确率高，但其计算复杂度高、速度慢，不易于实时应用。半全局匹配（semi ⁃ global matching，SGM）算法^[9]以其高精度、高效率的特点，已被广泛用于商用摄影测量软件^[10]，是用多个一维方向的聚合能量近似代替全局能量来求解最优视差。

基于对实际场景的分段平滑假设，SGM算法采用双参数P₁和P₂来建立相邻像素间的视差平滑约束，P₁和P₂的大小决定了匹配代价和平滑项的权重，影响着视差图的平滑程度。在实际应用中，由于航摄影像场景复杂，存在高楼林立、阴影交错的城区等易出现深度不连续和遮挡的区域，也存在农田、森林、水域等重复纹理或弱纹理区域，这给密集匹配带来了一定困难和挑战，若对所有像素使用固定的P₁和P₂难以顾及复杂的地形条件。原始的SGM方法试图使用梯度来克服这个问题，即在高梯度区域使用低惩罚而在低梯度区使用高惩罚。但梯度大的地方并不完全对应着深度不连续，易受噪声影响。文献[11]中的SURE方法采用Canny算子检测边缘，并在边缘处减小参数P₂，从而达到保持边缘的效果，但是Canny算子仅利用梯度信息获取物体边缘，仍对噪声敏感，且检测到的边缘由单像素连接而成，而边缘实际上很难保证准确定位^[12]。朱庆等^[1]采用局部窗口的梯度和标准差来提取纹理信息，并根据纹理信息来调整惩罚参数，在一定程度上减小了灰度噪声的影响，但是仅利用局部信息定量表达纹理特征仍存在局限性。Chuang等^[13]采用两步匹配策略：第1步使用U⁃SURF特征匹配边缘像素，使用自适应支持窗口方法匹配非边缘像素，获得初始视差图；第2步利用初始匹配代价调整惩罚参数，使用SGM方法获得最终视差图，该方法不仅保留了清晰的边缘又不会引入冗余的噪声，但没有考虑效率。

1 边缘约束的SGM方法

本文提出了一种边缘约束的航摄立体影像密集匹配方法，整体流程如图 1所示。由于核线关系将同名点约束在同名核线上，将同名点的搜索范围从二维降至了一维，极大的提高了匹配效率，因此首先生成核线影像对，在核线影像上进行密集匹配。

1.1 基于点互信息的边缘提取

点互信息^[14]在图像处理领域被用来衡量像素特征之间联系的紧密程度，利用属于同一物体中的像素比不属于同一物体的像素的统计相关性更高的特性，判断从属于两个特征的两像素是否属于同一物体。若特征间的关联度高，则从属于这两个特征的像素点间的相似度高，应属于同一物体；反之，若特征间的关联度低，从属于这两个特征的像素点间的相似度低，应属于不同物体。像素点互信息是对影像全局特征的描述，具有更强的稳定性和抗噪性，通过判断特征间的点互信息值可将看似特征完全不同但实属同一目标的像素点融合，从而获得影像中目标物体的真实边缘信息。

根据上述点互信息的特点，本文采用基于点互信息的边缘提取方法，在一张影像中对相邻像素进行采样，随机采样提取N个点对，利用得到的N个特征对和N个点互信息值生成随机森林决策树，以获取影像中任意特征对的点互信息值，即对整张影像构建相似性矩阵，然后计算边缘概率，获取边缘。

为影像中的每个像素定义一个特征向量$f$, 对于一张影像中的一个点对$(i, j)$, 点对中的点可看作一个随机变量, 点对映射到特征域后可提取一个特征对$\left(f_{i}^{1}, f_{j}^{1}\right)$, 该特征对的点互信息可由式(1)得到:

$ \operatorname{PMI}_{\rho}\left(f_{i}^{1}, f_{j}^{1}\right)=\lg \frac{P\left(f_{i}^{1}, f_{j}^{1}\right)^{\rho}}{P\left(f_{i}^{1}\right) P\left(f_{j}^{1}\right)} $

(1)

式中, $\rho=1.25; P\left(f_{i}{ }^{1}, f_{j}{ }^{1}\right)$为特征$f_{i}{ }^{1}$和$f_{j}{ }^{1}$的联合概率$; P\left(f_{i}{ }^{1}\right)、P\left(f_{j}{ }^{1}\right)$分别为特征$f_{i}{ }^{1}$和$f_{j}{ }^{1}$的独立概率。采用非参数核密度估计^[15]来推导影像特征间的联合分布概率和独立分布概率。

为了相似性度量的稳定性和准确性，本文采用像素的颜色特征（L×a×b颜色空间）及3×3窗口内像素RGB（red⁃green⁃blue）颜色的对角协方差矩阵作为特征集，因此每个像素有两个三维的特征向量。则采用如式（2）所示的相似度矩阵对像素i，j间的相似性进行定量描述。

$ \boldsymbol{W}_{i, j}=e^{\sum\limits_{k=1}^{M} \operatorname{PMI}_{\rho}\left(f_{i}^{k}, f_{j}^{k}\right)} $

(2)

式中, $\boldsymbol{W}$为相似度矩阵; $f_{i}$与$f_{j}$分别为点$i$和$j$的特征向量; $k$表示特征; $\boldsymbol{M}$表示特征向量的个数。获取相似度矩阵$\boldsymbol{W}$后, 采用文献[15]中的方法计算得到每个像素的边缘概率。

1.2 边缘约束的立体影像密集匹配

SGM方法用多个方向的聚合能量来近似代替全局能量，然后采用WTA（winner-take-all）方法获取视差。通常，能量函数由数据项和平滑项两部分构成，表达式为：

$ \begin{aligned} E(d)=& {\sum\nolimits_p} C\left(p, d_{p}\right)+\\ & \sum_{q \in N_{p}} P_{1} T\left[\left|d_{p}-d_{q}\right|=1\right]+\\ & \sum_{q \in N_{p}} P_{2} T\left[\left|d_{p}-d_{q}\right|>1\right] \end{aligned} $

(3)

式中, 第1项为数据项, 计算像素$p$在视差$d_{p}$时的匹配代价; 第2项为平滑项, 惩罚相邻像素$(p, q)$的深度不连续性; $N_{p}$为像素$p$的邻域点集, $T[\bullet]$为Potts模型, $T[x]=\left\{\begin{array}{ll}0 & x \text { 为假} \\ 1 & x \text { 为真}\end{array}\right.$ 。$\mathrm{SGM}$中的平滑项涉及两个惩罚参数$P_{1}$和$P_{2}\left(P_{2}>P_{1}\right)$, 若相邻像素的视差差异等于1, 则使用$P_{1}$对其进行惩罚, 以允许倾斜平面或曲面的视差微小变化; 若相邻像素的视差差异大于1, 则使用$P_{2}$对其进行惩罚, 以惩罚深度不连续。

立体影像密集匹配的核心工作是在左右影像中找到同一物体的对应同名点。匹配代价是判断两个像素是否相似的度量函数，是估算视差的依据。Scharstein等^[4]对当前常用的匹配代价函数进行了研究和评估，经过实验发现基于Census变换的匹配代价^[16]对图像辐射变化不敏感且可以抵抗较小的局部不规则噪声，具有最稳健的匹配性能。但Census变换是根据两个像素间的灰度大小关系来计算匹配代价的，没有考虑特定的像素灰度值，在纹理简单的区域具有重复的结构和较低的区分度；另外，Census变换对中心像素的依赖性很强，当中心像素出现噪声时，匹配结果是不可预测的^[17]。针对上述问题，本文利用中心像素的支持窗口内的邻域均值代替中心像素的灰度来进行Census变换，以避免中心像素受噪声干扰生成错误的字符串；另外，考虑到梯度可以抵抗部分噪声，且使用了邻近像素灰度值之间的连续性约束，以及梯度对深度不连续与一定的识别能力，本文引入梯度算子与Census变换融合构建新的匹配代价。匹配代价计算公式为：

$ C_{C e n}\left(p, d_{P}\right)=\mathrm{H}\left(T_{L}(p), T_{R}\left(p^{\prime}\left(p, d_{p}\right)\right)\right) $

(4)

$ C_{G r a}\left(p, d_{p}\right)=\left|G_{L}(p, h(p))-G_{R}\left(p^{\prime}\left(p, d_{p}\right)\right)\right| $

(5)

$ C\left(p, d_{p}\right)=\rho\left(C_{C e n}\left(p, d_{p}\right), \lambda_{\text {Cen }}\right)+ \\ \rho\left(C_{G r a}\left(p, d_{p}\right), \lambda_{G r a}\right) $

(6)

$ \rho(c, \lambda)=1-\exp \left(-\frac{c}{\lambda}\right) $

(7)

式中, $p^{\prime}\left(p, d_{p}\right)$表示当视差为$d_{p}$时, 像素$p$为在右影像中的待匹配像素; $T(\bullet)$表示以邻域均值作为判断阈值的Census变换; $H(\bullet)$表示Hamming距离; $G(\bullet)$表示影像的一阶方向导数; 下标$L$和$R$分别表示左右影像; $h(q)$表示通过像素$p$的核线, 归一化函数$\rho(c, \lambda)$避免了匹配代价倾向于某一代价; $\lambda$用来平衡两个匹配代价的权重, 有益于结合Census变换和梯度的优点, 更好地适应灰度失真和噪声, 减少深度不连续处的匹配误差。

由于单一的固定惩罚参数无法适应航摄影像中复杂的场景，因此本文在能量函数模型中增加边缘约束，顾及像素的边缘属性对每个像素选择不同的惩罚参数P 1和P 2：若当前像素的边缘概率低，应施加大的惩罚参数保证视差平滑；若当前像素的边缘概率高，则施加小的惩罚参数允许视差跳跃，防止过度平滑。一种简单的调整参数的方式是根据像素的边缘概率逐像素计算惩罚参数，但在实际应用中，该方法存在诸多问题：一是需要额外的内存和存取操作；二是同一物体上的像素的边缘概率可能不同，这可能会造成错误。因此，本文根据每个像素的边缘概率，设定一个阈值，采用二值化处理方法，令单个像素对应一个字节标识边缘属性，见式（8），此时边缘不再是单一像素连接而成，而是由高边缘概率像素组成的邻近真实边缘像素附近的几个像素，这就避免了单像素边缘定位不准确的问题。

$ \rho(p)=\left\{\begin{array}{l} 255, \text { 边缘像素 } \\ 0, \text { 其他 } \end{array}\right. $

(8)

在进行某一方向的代价聚合时, 根据每个像素的边缘属性自适应选择惩罚参数$P_{1}$和$P_{2}$, 具体来说, 若$\rho(p)=255$, 即当前像素$p$的属性为边缘像素, 则对相邻像素的视差变化选取小的惩罚参数允许视差突变; 若$\rho(p)=0$, 即当前像素$p$的属性为非边缘像素, 则对相邻像素的视差变化选取大的惩罚参数保证视差平滑。改进后的能量函数为：

$ \begin{aligned} &E(d)=\sum\nolimits_{p} C\left(p, d_{p}\right)+ \\ &\sum\limits_{q \in N_{p}} P_{1} T\left[\left|d_{p}-d_{q}\right|=1 \cap \rho(p)=255\right]+ \\ &\sum\limits_{q \in N_{p}} P_{2} T\left[\left|d_{p}-d_{q}\right|>1 \cap \rho(p)=255\right]+ \\ &\sum\limits_{q \in N_{p}} P_{1}{ }^{\prime} T\left[\left|d_{p}-d_{q}\right|=1 \cap \rho(p)=0\right]+ \\ &\sum\limits_{q \in N_{p}} P_{2}{ }^{\prime} T\left[\left|d_{p}-d_{q}\right|>1 \cap \rho(p)=0\right] \end{aligned} $

(9)

其中, $ P_{2}^{\prime}>P_{2}, P_{1}^{\prime}>P_{1}, P_{2}>P_{1} $。

2 试验结果及其分析 2.1 试验设计

本文使用来自ISPRS/EuroSDR的“Benchmark on High Density Aerial Image Matching”项目所提供的Vaihingen公开数据集，采用UltraCamX相机获取3条航带共36张像片，像幅为9 420×14 430像素，相机焦距为100 mm，航高为2 900 m，设计航向重叠度为和旁向重叠度均为60%，地面分辨率为20 cm。受计算机运行内存的限制，原始大小的影像视差范围过大不能直接进行匹配，本文采用影像分块的策略，将原始核线影像裁剪为360×955像素大小的子图像，视差范围为256像素。

试验环境为：Windows10 64bit操作系统，处理器为Intel（R）Core（TM）i5⁃6300HQCPU@2. 30 GHz，8 GB内存。

试验的参数设置为：SGM方法按照8个方向进行代价累加，采用7×7窗口内的Census变换计算海明距离，λ_Cen=10，λ_Gra=50，惩罚参数设置如表 1所示。

2.2 试验结果与分析

为了验证边缘约束的密集匹配算法的有效性，将本文方法与SGM方法和SURE软件的点云结果进行对比。如图 2所示，本文截取试验区域中包含房屋、平地区域的3张影像块a、b、c，并选取一些检查线和检查面，其中，检查线对应房屋边缘，检查面对应屋顶、地面等平坦区域。

首先，本文采用M3C2计算不同方法生成的点云的Z坐标与参考Z坐标之间的差距以定量评价点云精度。其中参考Z坐标是通过将DSM（digital sur⁃ face model）栅格化获取的^[18]。表 2为本文提出的边缘约束的匹配方法与SGM方法在各个影像块中生成点云的Z坐标的平均误差对比。可以看出本文方法在平坦区域和深度不连续处的匹配精度相较于SGM方法有较大提高。

图 3定性比较了使用本文方法、SGM和SURE重建的点云图。可以看出，SURE与SGM相比，虽然在某些房屋边缘的表现优于后者，如图 3中的实线椭圆，但是也导致某些房屋的边角特征缺失，如图 3中的虚线椭圆。而本文提出方法不仅改善了SGM方法得到的点云边缘缺失、前景视差延伸到背景的问题，并且没有因为加入了边缘约束导致更多的边角特征缺失，如图 3中矩形区域。因此，从整体来看，本文方法优于SGM方法和文献[11]的SURE方法，不仅提高了匹配精度，且生成的点云在平坦区域保持平滑，在房屋边缘处更加完整且清晰。

3 结束语

本文在传统的SGM方法的基础上，提出了一种边缘约束的航摄立体影像密集匹配算法：首先通过计算颜色特征和协方差矩阵特征集下的图像像素点互信息构建整张影像的相似性矩阵，获取影像的边缘信息；然后根据像素的边缘属性自适应选择逐像素的惩罚参数。试验表明，本文方法与原始SGM相比，生成的点云精度有所提高，重建点云不仅在道路、屋顶等平坦区域保持平滑，且更有效地保持了房屋的边角特征。