CN112435267A - 一种高分辨率城市卫星立体图像的视差图计算方法 - Google Patents

Abstract

本发明属于高分辨率卫星影像三维重建技术领域，具体涉及一种高分辨率城市卫星立体图像的视差图计算方法。本发明解决了现有卫星影像稠密立体匹配方法在遮挡区域出现匹配漏洞以及在视差不连续区域出现误匹配的问题。本发明在Census代价和Patch代价加权和的基础上引入对于纹理特征敏感的PHOW特征计算卫星图像的匹配代价，提高了光照变化区域、弱纹理区域和建筑物犄角等区域的匹配代价的描述能力。本发明将匹配代价优化问题转化为区域标记问题，通过区域边界和区域模型来约束区域内匹配代价的优化，能够定位深度不连续的视差断层区域，填补视差图中的空洞，区域边界和模型是动态更新的，有效地提高了区域边界和区域模型的准确性。

Description

一种高分辨率城市卫星立体图像的视差图计算方法

技术领域

本发明属于高分辨率卫星影像三维重建技术领域，具体涉及一种高分辨率城市卫星立体图像的视差图计算方法。

背景技术

城市三维信息对于城市智能交通，城市智能规划和灾害防控等应用都起到了至关重要的作用，因此城市三维地理信息具有巨大的经济价值和国防战略意义。因此，通过高分辨率城市卫星影像进行城市三维重建的方法已经受到了很大的关注。该方法一般包括立体校正，稠密视差图计算和三角测量三个过程。作为三维重建的核心，已经提出的适用于卫星立体图像对的稠密视差图构建方法大致可以分成两类，一类是基于半全局匹配(SGM)的方法，另一类是基于区域的匹配方法。

原始半全局匹配(SGM)方法首先通过互信息熵计算匹配代价，然后利用代价聚合获得准确的视差图。存在利用Census变换替代互信息熵计算匹配代价的方法，进而获得视差结果。但是Census是通过构建中心像素点与周围像素点灰度值的大小关系来表示当前像素的特征的，利用大小关系计算不同像素之间的相似程度经常导致匹配代价丧失匹配的单一性。基于区域的匹配方法大多利用切块匹配(Patch Match，PM)获取匹配代价。PM计算匹配代价是以参考图像的像素为中心开辟支撑窗口，然后在第二图像中移位该支撑窗口将，最低代价像素点作为匹配点。但是，该过程存在隐含的假设，即支撑窗口内的所有像素具有恒定的视差，这导致支撑窗口内的其他像素视差被中心像素“归一化”。在视差不连续区域PM方法误差会很大。

此外，在匹配代价优化阶段，动态规划的方法被用于减少“归一化”的影响。也有研究者利用自适应支持权重的导向滤波器法对代价量进行滤波，目的是还原深度不连续区域的真实物体边缘。但是高分辨率城市卫星影像是卫星在高速运动条件下获取的，因此在卫星立体图像中存在更多倾斜的建筑，导致图像中存在大量遮挡区域，上述两类方法在在遮挡区域无法获取可靠的匹配信息，因此获取的视差图包含很多视差空洞。同时，在大幅面的城市卫星影像中存在大量的密集的建筑，不同高度的建筑物之间、建筑物屋顶和立墙之间以及立墙与地面之间都形成了视差不连续的视差断层。通过上述两类方法获得的视差图在深度不连续的视差断层区域都存在大量的误差。这些视差空洞和视差不连续区域处的误差都严重影响了城市三维重建的精度。

发明内容

本发明的目的在于提供能够有效的填充遮挡区域的视差空洞，并优化视差不连续区域的视差结果，为高分辨率城市卫星影像三维重建提供准确的稠密视差结果的一种高分辨率城市卫星立体图像的视差图计算方法。

本发明的目的通过如下技术方案来实现：包括以下步骤：

步骤1：通过超像素分割方法将立体校正的高分辨率城市卫星立体图像对应的参考图像划分成贴合实际目标轮廓的超像素区域，得到n个超像素区域集合

步骤2：使用基于PHOW特征的下采样匹配方法在每个超像素区域内获得分布均匀的匹配点，并利用多模型拟合的方法计算每个超像素区域对应的单应模型；

步骤3：利用三种相似性度量的联合加权的方法计算高分辨率城市卫星立体图像中同名点之间的匹配代价COST；

步骤4：使用双重传播动态模型优化方法，在超像素区域模型的约束下对原始匹配代价的进行迭代优化，实现区域模型的动态更新，以解决遮挡区域和视差不连续区域对匹配造成的影响，并通过最小化能量函数的方式获得优化后的城市卫星立体图像对应的视差图；

首先执行外部传播，通过置信传播BP构建外传播概率图模型，超像素块作为图的节点，超像素之间的相邻关系作为图的边；为了解决因BP结构中的“环结构”导致消息的循环传播，在每次迭代中会通过区域面积筛选出对应的“目标节点”，这些节点只作为消息的接收者，不会向外传播消息，“目标节点”融合成新的“非目标节点”之后会继续向外传播消息，具体传播步骤如下：

步骤4.1：计算外部传播中的消息；

其中，

表示节点j传到相邻节点i第t次迭代的消息；φ_i(x_i)为φ_i(x_i,y_i)的简写，表示似然函数；ψ_j→i(·)表示势函数；

表示排除节点x_i的x_j邻域系统

中所有传入节点x_j的消息的乘积；g表示超像素区域i中的像素索引；f_g表示在超像素区域i中的第g个像素的图像特征；δ_i和μ_i分别为在超像素区域i中的标准差和均值；DisG(·)和DisH(·)分别为超像素区域i与超像素区域j之间的灰度距离和单应模型距离；τ_g和τ_h分别对应于灰度距离和单应模型距离的截断值，被用于控制区域间视差差值在合理范围内；

和

分别为超像素区域i与超像素区域j内的像素点集合；I_p为像素点p的灰度值；Num(p)为超像素区域内像素点p的数量；

和

分别为超像素区域i与超像素区域j内下采样点集合，

||H_j·l^t,l^t′||₂表示下采样点t的坐标l^t在单应模型H_j作用下获得的坐标与匹配点坐标l^t′之间的欧氏距离，下采样点t和单应模型H_j分别属于超像素区域i和超像素区域j；

步骤4.2：计算外部传播的置信度；节点i的置信度b_i(x_i)是所有传入节点i的消息与节点i似然函数的乘积；

其中，k为归一化常数；

步骤4.3：外部传播迭代；“目标节点”与置信度最大的节点融合，形成新的“非目标节点”；如果结构中存在“目标节点”，则返回步骤4.1；否则，终止迭代，执行步骤4.4，进入内部传播；

步骤4.4：执行内部传播的惩罚传播；内部传播对外部传播形成的各个区域中的每个像素的匹配代价都沿着扫描线的路径进行迭代更新，内部传播包括八个方向的区域扫描线迭代；通过增加对区域边缘的惩罚来保护视差断层，同时也平滑区域内部代价；区域内沿扫描线的方向的惩罚传播计算公式如下：

其中，

为超像素区域S中迭代计算的匹配代价；C^s(p,d)为超像素区域S中当前点匹配代价；N^r(p)为超像素区域S中r路径上像素p的邻域；

为超像素区域S中上一轮计算的匹配代价，q为位于超像素区域S中路径r上像素点p的前一个像素点；V^s(D_p,D_q)为超像素区域S中的惩罚函数，计算方法如下：

其中，p₁，p₂和p₃是惩罚因子，并且满足p₂＞p₃＞p₁；p₁惩罚超像素区域内视差变化为1个像素的情况；p₂惩罚位于视差断层处的代价；p₃惩罚超像素区域内视差变化大于1个像素的情况；惩罚函数既支持超像素内相邻像素之间的相似的视差和轻微的视差变化，也允许在超像素边缘处出现大的视差变化；

步骤4.5：通过对每个方向上的代价求和获得超像素区域内优化后的能量S^s(p,d)；通过对所有区域内的代价求和获得整张图像的能量函数E(D)；

步骤4.6：通过最小化能量函数的方式获得优化后的城市卫星立体图像对应的视差图；最小化能量函数的计算公式为：

D＝arg min E(D)。

本发明还可以包括：

所述的步骤2中使用基于PHOW特征的下采样匹配方法在每个超像素区域内获得分布均匀的匹配点，并利用多模型拟合的方法计算每个超像素区域对应的单应模型的方法具体为：

步骤2.1：使用基于PHOW特征的下采样匹配方法在M个像素(x,y)处获取描述子PHOW(x,y)；

PHOW(x,y)＝SIFT₄(x,y)+SIFT₈(x,y)+SIFT₁₂(x,y)

其中，SIFT₄(x,y)、SIFT₈(x,y)和SIFT₁₂(x,y)分别是半径为4、8、12的SIFT特征；

步骤2.2：通过最小特征距离与次小特征距离之比的方式获得符合要求的匹配点，构成下采样匹配点集合

判断条件如下式所示：

其中，Min(·)和Sec(·)分别为最小值函数和次小值函数；PHOW(x_l,y_l)表示左视图中点(x_l,y_l)处的PHOW特征，并且

表示右视图中第i个点(x_r,y_r)处的PHOW特征，并且

Thr表示阈值；

步骤2.3：利用模型拟合的方法对所有超像素区域进行模型拟合，获得每个区域对应的单应模型；单应矩阵是3×3的矩阵H：

通过变换公式Hp₁＝p₂，把p₁＝[x₁,y₁,1]^T变成一个新的点p₂＝[x₂,y₂,1]^T；其中，p₁和p₂都是齐次坐标，并且在一个单应矩阵H中，令右下角的元素h₃₃＝1进行归一化，单应性矩阵只存在8个自由元素h₁₁,h₁₂,h₁₃,h₂₁,h₂₂,h₂₃,h₃₁,h₃₂，因此可以获得下式：

至少需要4对不在同一条直线上匹配点(x₁,y₁)和(x₂,y₂)使用SVD分解求解H的八个元素，如下所示：

所述的步骤3中利用三种相似性度量的联合加权的方法计算高分辨率城市卫星立体图像中同名点之间的匹配代价COST的方法具体为：

步骤3.1：计算Patch匹配代价COST_P；

其中，W_p表示以像素点p为中心的方形窗口，S_W是该窗口的尺寸；κ是用于平衡灰度距离与坐标距离之间的影响比重；||I_p-I_q||₁表示点p灰度与点q灰度的L1距离；||C_p-C_q||₂表示点p坐标与点q坐标的L2距离；ω是用于平衡灰度与梯度之间的比重；||I_q-I_q′||₁表示点q与点q′处的灰度L1距离，

表示点q与点q′处的梯度L1距离；q′为q在另一视图中的匹配点；参数τ_gray和参数τ_grad表示对于灰度距离与梯度距离的截断值，对于遮挡区域具有鲁棒性；

步骤3.2：计算Census匹配代价COST_C；

Census匹配代价是通过计算像素点之间bit字符串的汉明距离得到的，令左视图和右视图中支撑窗口中心像素点的Census变换值分别为Census(x_l,y_l)和Census(x_r,y_r)，则两像素点之间的Census代价为：

COST_C＝Dist_HAMMING(Census(x_l,y_l),Census(x_r,y_r))

其中，Dist_HAMMING(str1,str2)运算表示两个bit字符串str1和str2中对应位置不同bit值的个数之和，汉明距离越小说明两支撑窗口中心像素点的相似程度越高；基于左右视图的Census变换为：

其中，I_p(·)为待匹配像素点p(x,y)的灰度值；m×n为以p(x,y)为中心的矩形窗口的尺寸；

表示对经过δ(·)运算后所得的结果的前后连接运算；

步骤3.3：计算PHOW匹配代价COST_PH；

将像素点上不同半径圆形支撑块内的SIFT描述子作为该像素的PHOW特征，采用单一半径的SIFT特征，每个点由一个128维SIFT描述符表示，如下所示

PHOW(x,y)＝SIFT₈(x,y)

通过计算对应点PHOW描述符之间的欧几里得距离，作为对应点之间的PHOW匹配代价值COST_PH，计算公式如下：

COST_PH＝||PHOW(x_l,y_l)-PHOW(x_r,y_r)||₂

步骤3.4：将三种代价联合加权作为最终的高分辨率城市卫星立体图像中同名点之间的匹配代价COST；

COST＝α·COST_P+β·COST_C+γ·COST_PH

其中，α、β和γ分别表示代价权重。

本发明的有益效果在于：

本发明解决了现有卫星影像稠密立体匹配方法在遮挡区域出现匹配漏洞以及在视差不连续区域出现误匹配的问题。本发明在Census代价和Patch代价加权和的基础上引入对于纹理特征敏感的PHOW特征计算卫星图像的匹配代价，以提高光照变化区域、弱纹理区域和建筑物犄角等区域的匹配代价的描述能力。本发明将匹配代价优化问题转化为区域标记问题，通过区域边界和区域模型来约束区域内匹配代价的优化。区域边界能够定位深度不连续的视差断层区域；区域模型能够填补视差图中的空洞。而且区域边界和模型并非是一成不变的，而是动态更新的，这有效地提高了区域边界和区域模型的准确性。

附图说明

图1为本发明的流程图。

图2为本发明中卫星立体图像对应的局部区域的单应关系示意图。

图3为本发明中Census变换示意图。

图4为本发明中外部传播示意图。

图5为本发明中内部传播示意图。

图6(a)为高分辨率城市卫星立体图像对应的参考图像。

图6(b)为图6(a)中方框区域的放大图。

图6(c)为图6(b)对应的视差图。

图6(d)为图6(c)的3D显示图。

具体实施方式

下面结合附图对本发明做进一步描述。

本发明属于高分辨率卫星影像三维重建技术领域，具体涉及一种高分辨率城市卫星立体图像的视差图计算方法。本发明解决了现有卫星影像稠密立体匹配方法在遮挡区域出现匹配漏洞以及在视差不连续区域出现误匹配的问题。

城市三维信息对于城市智能交通，城市智能规划和灾害防控等应用都起到了至关重要的作用，因此城市三维地理信息具有巨大的经济价值和国防战略意义。因此，通过高分辨率城市卫星影像进行城市三维重建的方法已经受到了很大的关注。该方法一般包括立体校正，稠密视差图计算和三角测量三个过程。作为三维重建的核心，已经提出的适用于卫星立体图像对的稠密视差图构建方法大致可以分成两类，一类是基于半全局匹配(SGM)的方法，另一类是基于区域的匹配方法。

原始半全局匹配(SGM)方法首先通过互信息熵计算匹配代价，然后利用代价聚合获得准确的视差图。存在利用Census变换替代互信息熵计算匹配代价的方法，进而获得视差结果。但是Census是通过构建中心像素点与周围像素点灰度值的大小关系来表示当前像素的特征的，利用大小关系计算不同像素之间的相似程度经常导致匹配代价丧失匹配的单一性。基于区域的匹配方法大多利用切块匹配(Patch Match，PM)获取匹配代价。PM计算匹配代价是以参考图像的像素为中心开辟支撑窗口，然后在第二图像中移位该支撑窗口将，最低代价像素点作为匹配点。但是，该过程存在隐含的假设，即支撑窗口内的所有像素具有恒定的视差，这导致支撑窗口内的其他像素视差被中心像素“归一化”。在视差不连续区域PM方法误差会很大。本发明在Census代价和Patch代价加权和的基础上引入对于纹理特征敏感的PHOW特征计算卫星图像的匹配代价，以提高光照变化区域、弱纹理区域和建筑物犄角等区域的匹配代价的描述能力。

此外，在匹配代价优化阶段，动态规划的方法被用于减少“归一化”的影响。也有研究者利用自适应支持权重的导向滤波器法对代价量进行滤波，目的是还原深度不连续区域的真实物体边缘。但是高分辨率城市卫星影像是卫星在高速运动条件下获取的，因此在卫星立体图像中存在更多倾斜的建筑，导致图像中存在大量遮挡区域，上述两类方法在在遮挡区域无法获取可靠的匹配信息，因此获取的视差图包含很多视差空洞。同时，在大幅面的城市卫星影像中存在大量的密集的建筑，不同高度的建筑物之间、建筑物屋顶和立墙之间以及立墙与地面之间都形成了视差不连续的视差断层。通过上述两类方法获得的视差图在深度不连续的视差断层区域都存在大量的误差。这些视差空洞和视差不连续区域处的误差都严重影响了城市三维重建的精度。本发明则是将匹配代价优化问题转化为区域标记问题，通过区域边界和区域模型来约束区域内匹配代价的优化。区域边界能够定位深度不连续的视差断层区域；区域模型能够填补视差图中的空洞。而且区域边界和模型并非是一成不变的，而是动态更新的，这有效地提高了区域边界和区域模型的准确性。

本发明的目的在于提供能够有效的填充遮挡区域的视差空洞，并优化视差不连续区域的视差结果，为高分辨率城市卫星影像三维重建提供准确的稠密视差结果。

首先通过超像素分割方法将立体校正的高分辨率城市卫星立体图像对的参考图像划分成贴合实际目标轮廓的超像素区域。然后使用基于金字塔直方图描述子(PyramidHistogram Of visual Words，PHOW)特征的下采样匹配方法在每个超像素区域内获得分布均匀的匹配点，并利用多模型拟合的方法计算每个超像素区域对应的模型。之后利用基于切块(Patch)代价、Census代价和金字塔直方图描述子(Pyramid Histogram Of visualWords，PHOW)代价的联合匹配代价计算方法计算高分辨率城市卫星立体图像中同名点之间的匹配代价。最后使用双重传播动态模型优化方法，在超像素区域模型的约束下对原始匹配代价的进行迭代优化，实现区域模型的动态更新，以解决遮挡区域和视差不连续区域对匹配造成的影响，并通过最小化能量函数的方式获得优化后的城市卫星立体图像对的视差图。本发明的流程图如图1所示。

本发明的有益效果是：本发明针对传统立体匹配方法应用在城市卫星立体图像上出现的视差图漏洞和视差断层区域误差较大的问题，本发明结合城市卫星图像特点，采用区域模型约束的思想实现城市卫星影像视差图的高精度计算。对比传统卫星图像视差图计算方法，本发明具有以下优点：

第一，本发明利用基于PHOW特征的下采样匹配方法以获得均匀分布的匹配点对。该方法克服了传统特征匹配方法具有的匹配点分布以及匹配点数量均不受控制的问题。利用此方法，每个超像素区域都能获得匹配点对。因此每个超像素区域的模型拟合都能获得充足的匹配点。

第二，本发明利用基于Patch代价、Census代价和PHOW代价的联合加权匹配代价计算方法计算高分辨率城市卫星立体图像中同名点之间的相似程度即匹配代价。联合加权匹配代价能够克服城市卫星影像中的建筑物犄角、弱纹理、相似纹理和光照变化等因素对视差计算获取造成的影响。

第三，本发明使用双重传播动态模型优化方法实现区域内匹配代价、区域划分以及区域模型的迭代优化，为区域视差图计算提供更加准确的区域边界和区域模型的约束。在双重传播动态模型优化中，先利用外部传播方法实现区域划分及其对应的区域模型的迭代优化。然后利用内部传播方法，在超像素区域内实现以像素为单位的扫描线的惩罚传播，获得像素级代价的优化结果。在双重传播的作用下，遮挡区域能够根据区域模型对视差进行了合理的推理，使视差不连续区域能够利用区域边界的约束获得准确的视差断层，最终获得更加完整和准确的城市卫星视差图。

图2是本发明中卫星立体图像对的局部区域的单应关系示意图，其中O-XYZ为世界坐标系，O′-X′Y′Z′为参考图像的摄像机坐标系，I_r和I_t分别为参考图像和目标图像，I_r中红色多边形区域为某块超像素区域，蓝色虚线箭头表示该区域内下采样匹配点对，H为该区域对应的单应矩阵；

图3是本发明中Census变换示意图，其中窗口大小为5×5，红色位置为窗口中心像素，灰度值为100，深蓝色位置为灰度值小于100的像素点，浅蓝色位置为灰度值大于或等于100的像素点。由公式12对当前窗口进行Census变换，通过公式13求出窗口内所有非中心像素点的bit值，经过前后连接获得像素点p(x,y)的Census变换值；

图4是本发明中外部传播图，其中I图为包含7个超像素块的预分割结果示意图，实线为超像素分割线，中间区域为“目标区域”对应G图的“目标节点”。在G图中，空白节点表示隐含节点，黑色节点表示显式节点，虚线表示隐含节点与显式节点之间的似然函数φ_i(x_i,y_i)，实线表示存在一对势函数ψ_ij与ψ_ji，箭头表示节点与节点之间的消息传递方向；

图5是本发明中内部传播示意图，其中(a)-(h)为两个超像素区域按照八个方向进行内部传播的示意图，图中黑色曲线为超像素区域边界，红色点与蓝色点表示位于不同超像素区域上的像素点。

实施例1：

本发明首先通过SLIC超像素分割方法将立体校正的卫星影像划分成贴合实际目标轮廓的超像素区域。然后通过基于PHOW特征的下采样匹配和多模型拟合的方法初始化超像素区域模型。初始匹配代价的构建是通过联合三种匹配代价的加权和实现的。最后通过双重传播动态模型优化技术对初始匹配代价进行迭代优化，并通过最小化能量函数获得视差图。本发明流程图如图1所示。与现有方法相比，在卫星影像遮挡区域本发明方法能够根据区域模型对视差进行合理推理，在视差不连续区域能够依据分割边界获得准确的视差结果。本发明可以应用于高分辨率城市卫星影像三维重建技术领域。

步骤一、城市卫星参考图像超像素分割。本发明使用超像素分割算法对卫星影像进行超像素分割，以获得准确划分的不同视差区域块。本发明通过简单线性迭代聚类(SLIC)超像素分割方法对图像进行分割，获得n个超像素集合

其中i表示超像素的label。在本发明中超像素边缘能够实现对视差断层的“定位”，保证“相似视差”像素的出现在同一区域内。

步骤二、基于PHOW特征的下采样匹配和多模型拟合。首先，使用基于金字塔直方图描述子(Pyramid Histogram Of visual Words，PHOW)的下采样匹配方法能够在M个像素(x,y)处获取描述子PHOW(x,y)。

PHOW(x,y)＝SIFT₄(x,y)+SIFT₈(x,y)+SIFT₁₂(x,y) (1)

式中SIFT₄(x,y)，SIFT₈(x,y)和SIFT₁₂(x,y)分别是半径为4,8,12的SIFT特征。

然后，通过最小特征距离与次小特征距离之比的方式获得符合要求的匹配点，构成下采样匹配点集合

判断条件如下式所示：

Min(·)和Sec(·)分别为最小值函数和次小值函数；PHOW(x_l,y_l)表示左视图中点(x_l,y_l)处的PHOW特征，并且

表示右视图中第i个点(x_r,y_r)处的PHOW特征，并且

Thr表示阈值。

最后，利用模型拟合的方法对步骤一的超像素区域都进行模型拟合，以获得每个区域对应的单应模型。拟合卫星立体图像对应区域之间的几何模型，区域间单应关系示意图如图2所示。单应矩阵是3×3的矩阵H，

能够通过变换公式Hp₁＝p₂，把p₁＝[x₁,y₁,1]^T变成一个新的点p₂＝[x₂,y₂,1]^T。其中，p₁和p₂都是齐次坐标，并且在一个单应矩阵H中，令右下角的元素h₃₃＝1进行归一化，单应性矩阵只存在8个自由元素h₁₁,h₁₂,h₁₃,h₂₁,h₂₂,h₂₃,h₃₁,h₃₂。因此可以获得下式：

因此至少需要4对不在同一条直线上匹配点(x₁,y₁)和(x₂,y₂)使用SVD分解求解H的八个元素，如下所示：

步骤三、联合加权匹配代价计算。本发明利用三种相似性度量的联合加权的方法计算城市卫星立体图像对中同名点之间的匹配代价。首选，分别计算Patch代价、Census代价和PHOW代价，具体计算过程如下：

Step1：计算Patch匹配代价。Patch匹配代价包括灰度与梯度代价，其对于弱纹理和相似纹理区域具有鲁棒性，因此窗口中的灰度与梯度联合代价计算如下：

其中，W_p表示以像素点p为中心的方形窗口，S_W是该窗口的尺寸。权重函数W(p,q)被用于解决同一视图内边缘不连续和目标形变较大的问题。他通过计算窗口内的像素q与窗口中心像素p的颜色差异与欧几里德距离来度量窗口内像素所占比重，如果颜色越相似、与中心像素距离越近，则比重越大，权重计算如下：

其中，κ是用于平衡灰度距离与坐标距离之间的影响比重，||I_p-I_q||₁表示点p灰度与点q灰度的L1距离，||C_p-C_q||₂表示点p坐标与点q坐标的L2距离。

公式6中的函数ρ(q,q′)计算不同视图之间的像素差异：

其中，ω是用于平衡灰度与梯度之间的比重。||I_q-I_q′||₁表示点q与点q′处的灰度L1距离，

表示点q与点q′处的梯度L1距离，q′为q在另一视图中的匹配点。参数τ_gray和参数τ_grad表示对于灰度距离与梯度距离的截断值，对于遮挡区域具有鲁棒性。

Step2：计算Census匹配代价。Census匹配代价是基于左右视图的Census变换实现的，计算过程如下：

其中，I_p(·)为待匹配像素点p(x,y)的灰度值，m×n为以p(x,y)为中心的矩形窗口的尺寸。

表示对经过δ(·)运算后所得的结果的前后连接运算。公式9描述的一次Census变换如图3所示。从图中可知，Census变换中bit字符串的本质是对支撑窗口内中心像素点与周围像素点灰度级关系的直观描述，因此对于光照变化具有鲁棒性。Census代价是通过计算像素点之间bit字符串的汉明距离得到的。令左视图和右视图中支撑窗口中心像素点的Census变换值分别为Census(x_l,y_l)和Census(x_r,y_r)，则两像素点之间的Census代价为：

COST_C＝Dist_HAMMING(Census(x_l,y_l),Census(x_r,y_r)) (11)

其中，Dist_HAMMING(str1,str2)运算表示两个bit字符串str1和str2中对应位置不同bit值的个数之和。汉明距离越小说明两支撑窗口中心像素点的相似程度越高。

Step3：计算PHOW匹配代价。PHOW匹配代价将像素点上不同半径圆形支撑块内的SIFT描述子作为该像素的PHOW特征，因此PHOW代价对于建筑物犄角和建筑物边缘处的像素具有很强的描述能力。本发明将该描述子作为一种特殊的“匹配代价”来计算每个下采样像素点处的“PHOW代价”，为了加速运算，本专利采用单一半径的SIFT特征，因此每个点由一个128维SIFT描述符表示，如下所示：

PHOW(x,y)＝SIFT₈(x,y) (12)

本专利通过计算对应点PHOW描述符之间的欧几里得距离，作为对应点之间的代价值，计算公式如下：

COST_PH＝||PHOW(x_l,y_l)-PHOW(x_r,y_r)||₂ (13)

Step4：三种代价中间的联合加权作为最终的匹配代价，计算公式如下：

COST＝α·COST_P+β·COST_C+γ·COST_PH (14)

其中α、β和γ分别表示代价权重。COST_P、COST_C和COST_PH分别表示Patch代价、Census代价和PHOW代价。

步骤四、通过双重传播的动态优化方法，利用步骤二获得的区域模型对步骤三的初始匹配代价进行动态优化，最后通过最小化能量函数获得视差。

双重传播优化由外部传播和内部传播组成，外部传播通过区域间信息传播实现区域划分和区域模型的迭代优化。内部传播通过基于局部区域的扫描线惩罚传播实现像素级匹配代价优化，并最小化能量函数获得视差。

首先执行外部传播，通过置信传播(Belief Propagation，BP)构建外传播概率图模型，超像素块作为图的节点，超像素之间的相邻关系作为图的边。为了解决因BP结构中的“环结构”导致消息的循环传播，在每次迭代中会通过区域面积筛选出对应的“目标节点”，这些节点只作为消息的接收者，不会向外传播消息。“目标节点”融合成新的“非目标节点”之后会继续向外传播消息，包含“目标节点”的置信传播的图结构如图4所示。外部传播的具体传播步骤如下：

Step1：计算外部传播中的消息。外部传播中的消息计算公式如下：

其中

表示节点j传到相邻节点i第t次迭代的消息；φ_i(x_i)(φ_i(x_i,y_i)的简写)表示似然函数，ψ_j→i(·)表示势函数，

表示排除节点x_i的x_j邻域系统

中所有传入节点x_j的消息的乘积。

似然函数φ_i(x_i)计算已知节点融合希望的参数，即可能进行节点融合的反过程的参数。这在高斯形式的似然函数中体现为均值和方差。似然函数如下所示：

式中g表示超像素区域i中的像素索引，f_g表示超像素在区域i中的第g个像素的图像特征，δ_i和μ_i分别为在超像素区域i中的标准差和均值。

势函数ψ_j→i(·)计算如下：

ψ_j→i(x_i,x_j)＝u·min(Dis G(x_i,x_j),τ_g)+v·min(DisH(x_i,x_j),τ_h) (17)

其中DisG(·)和DisH(·)分别为超像素区域i与超像素区域j之间的灰度距离和单应模型距离。τ_g和τ_h分别对应于灰度距离和单应模型距离的截断值，被用于控制区域间视差差值在合理范围内。DisG(x_i,x_j)计算区域i和j之间灰度平均值的差值，具体形式如下：

其中，

和

分别为超像素区域i与超像素区域j内的像素点集合，I_p为像素点p的灰度值，Num(p)为超像素区域内像素点p的数量。DisH(x_i,x_j)计算区域i与j中下采样像素点到相邻区域单应模型之间的距离，目的是量化计算区域i与区域j之间的几何差异，具体形式如下：

其中，

和

分别为超像素区域i与超像素区域j内下采样点集合，因此有

||H_j·l^t,l^t′||₂表示下采样点t的坐标l^t在单应模型H_j作用下获得的坐标与匹配点坐标l^t′之间的欧氏距离，注意下采样点t和单应模型H_j分别属于超像素区域i和超像素区域j。另一项

具有相似的含义。

Step2：计算外部传播的置信度。节点i的置信度是所有传入节点i的消息与节点i似然函数的乘积。因此节点i的置信度b_i(x_i)的具体形式为：

其中k是一个归一化常数。

Step3：外部传播迭代。“目标节点”与置信度最大的节点融合，形成新的“非目标节点”。如果结构中存在“目标节点”，那么转到Step1；否则，终止迭代，进入内部传播。

Step4：执行内部传播的惩罚传播。内部传播对外部传播形成的各个区域中的每个像素的匹配代价都沿着扫描线的路径进行迭代更新，内部传播包括八个方向的区域扫描线迭代，示意图如图5所示。通过增加对区域边缘的惩罚来保护视差断层，同时也平滑区域内部代价。区域内沿扫描线的方向的惩罚传播计算公式如下：

其中，L_r(p,d)为S区域中迭代计算的匹配代价。C^s(p,d)为S区域中当前点匹配代价，N^r(p)为超像素区域s中r路径上像素p的邻域。

为S区域中上一轮计算的匹配代价，q为位于超像素区域s中路径r上像素点p的前一个像素点。V^s(D_p,D_q)S区域中惩罚函数，计算如下：

其中，p₁，p₂和p₃是惩罚因子，并且满足p₂＞p₃＞p₁。p₁惩罚超像素区域内视差变化为1个像素的情况，p₂惩罚位于视差断层处的代价，p₃惩罚超像素区域内视差变化大于1个像素的情况。惩罚函数既支持超像素内相邻像素之间的相似的视差和轻微的视差变化，因此也允许在超像素边缘处出现大的视差变化。

Step5：计算能量函数。通过对每个方向上的代价求和获得超像素区域内优化后的能量：

通过对所有区域内的代价求和获得整张图像的能量函数，计算公式如下：

Step6：通过最小化能量函数的方式获得视差图，如图6所示，最小化能量函数的计算公式如下：

D＝arg min E(D) (25)

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (3)

1.一种高分辨率城市卫星立体图像的视差图计算方法，其特征在于，包括以下步骤：

步骤1：通过超像素分割方法将立体校正的高分辨率城市卫星立体图像对应的参考图像划分成贴合实际目标轮廓的超像素区域，得到n个超像素区域集合

步骤2：使用基于PHOW特征的下采样匹配方法在每个超像素区域内获得分布均匀的匹配点，并利用多模型拟合的方法计算每个超像素区域对应的单应模型；

步骤3：利用三种相似性度量的联合加权的方法计算高分辨率城市卫星立体图像中同名点之间的匹配代价COST；

步骤4：使用双重传播动态模型优化方法，在超像素区域模型的约束下对原始匹配代价的进行迭代优化，实现区域模型的动态更新，以解决遮挡区域和视差不连续区域对匹配造成的影响，并通过最小化能量函数的方式获得优化后的城市卫星立体图像对应的视差图；

首先执行外部传播，通过置信传播BP构建外传播概率图模型，超像素块作为图的节点，超像素之间的相邻关系作为图的边；为了解决因BP结构中的“环结构”导致消息的循环传播，在每次迭代中会通过区域面积筛选出对应的“目标节点”，这些节点只作为消息的接收者，不会向外传播消息，“目标节点”融合成新的“非目标节点”之后会继续向外传播消息，具体传播步骤如下：

步骤4.1：计算外部传播中的消息；

ψ_j→i(x_i,x_j)＝u·min(DisG(x_i,x_j),τ_g)+v·min(DisH(x_i,x_j),τ_h)

其中，

表示节点j传到相邻节点i第t次迭代的消息；φ_i(x_i)为φ_i(x_i,y_i)的简写，表示似然函数；ψ_j→i(·)表示势函数；

表示排除节点x_i的x_j邻域系统

中所有传入节点x_j的消息的乘积；g表示超像素区域i中的像素索引；f_g表示在超像素区域i中的第g个像素的图像特征；δ_i和μ_i分别为在超像素区域i中的标准差和均值；DisG(·)和DisH(·)分别为超像素区域i与超像素区域j之间的灰度距离和单应模型距离；τ_g和τ_h分别对应于灰度距离和单应模型距离的截断值，被用于控制区域间视差差值在合理范围内；

和

分别为超像素区域i与超像素区域j内的像素点集合；I_p为像素点p的灰度值；Num(p)为超像素区域内像素点p的数量；

和

分别为超像素区域i与超像素区域j内下采样点集合，

||H_j·l^t,l^t′||₂表示下采样点t的坐标l^t在单应模型H_j作用下获得的坐标与匹配点坐标l^t′之间的欧氏距离，下采样点t和单应模型H_j分别属于超像素区域i和超像素区域j；

步骤4.2：计算外部传播的置信度；节点i的置信度b_i(x_i)是所有传入节点i的消息与节点i似然函数的乘积；

其中，k为归一化常数；

步骤4.3：外部传播迭代；“目标节点”与置信度最大的节点融合，形成新的“非目标节点”；如果结构中存在“目标节点”，则返回步骤4.1；否则，终止迭代，执行步骤4.4，进入内部传播；

步骤4.4：执行内部传播的惩罚传播；内部传播对外部传播形成的各个区域中的每个像素的匹配代价都沿着扫描线的路径进行迭代更新，内部传播包括八个方向的区域扫描线迭代；通过增加对区域边缘的惩罚来保护视差断层，同时也平滑区域内部代价；区域内沿扫描线的方向的惩罚传播计算公式如下：

其中，

为超像素区域S中迭代计算的匹配代价；C^s(p,d)为超像素区域S中当前点匹配代价；N^r(p)为超像素区域S中r路径上像素p的邻域；

为超像素区域S中上一轮计算的匹配代价，q为位于超像素区域S中路径r上像素点p的前一个像素点；V^s(D_p,D_q)为超像素区域S中的惩罚函数，计算方法如下：

其中，p₁，p₂和p₃是惩罚因子，并且满足p₂＞p₃＞p₁；p₁惩罚超像素区域内视差变化为1个像素的情况；p₂惩罚位于视差断层处的代价；p₃惩罚超像素区域内视差变化大于1个像素的情况；惩罚函数既支持超像素内相邻像素之间的相似的视差和轻微的视差变化，也允许在超像素边缘处出现大的视差变化；

步骤4.5：通过对每个方向上的代价求和获得超像素区域内优化后的能量S^s(p,d)；通过对所有区域内的代价求和获得整张图像的能量函数E(D)；

步骤4.6：通过最小化能量函数的方式获得优化后的城市卫星立体图像对应的视差图；最小化能量函数的计算公式为：

D＝arg min E(D)。

2.根据权利要求1所述的一种高分辨率城市卫星立体图像的视差图计算方法，其特征在于：所述的步骤2中使用基于PHOW特征的下采样匹配方法在每个超像素区域内获得分布均匀的匹配点，并利用多模型拟合的方法计算每个超像素区域对应的单应模型的方法具体为：

步骤2.1：使用基于PHOW特征的下采样匹配方法在M个像素(x,y)处获取描述子PHOW(x,y)；

PHOW(x,y)＝SIFT₄(x,y)+SIFT₈(x,y)+SIFT₁₂(x,y)

其中，SIFT₄(x,y)、SIFT₈(x,y)和SIFT₁₂(x,y)分别是半径为4、8、12的SIFT特征；

步骤2.2：通过最小特征距离与次小特征距离之比的方式获得符合要求的匹配点，构成下采样匹配点集合

判断条件如下式所示：

其中，Min(·)和Sec(·)分别为最小值函数和次小值函数；PHOW(x_l,y_l)表示左视图中点(x_l,y_l)处的PHOW特征，并且

表示右视图中第i个点(x_r,y_r)处的PHOW特征，并且

Thr表示阈值；

步骤2.3：利用模型拟合的方法对所有超像素区域进行模型拟合，获得每个区域对应的单应模型；单应矩阵是3×3的矩阵H：

通过变换公式Hp₁＝p₂，把p₁＝[x₁,y₁,1]^T变成一个新的点p₂＝[x₂,y₂,1]^T；其中，p₁和p₂都是齐次坐标，并且在一个单应矩阵H中，令右下角的元素h₃₃＝1进行归一化，单应性矩阵只存在8个自由元素h₁₁,h₁₂,h₁₃,h₂₁,h₂₂,h₂₃,h₃₁,h₃₂，因此可以获得下式：

至少需要4对不在同一条直线上匹配点(x₁,y₁)和(x₂,y₂)使用SVD分解求解H的八个元素，如下所示：

3.根据权利要求1或2所述的一种高分辨率城市卫星立体图像的视差图计算方法，其特征在于：所述的步骤3中利用三种相似性度量的联合加权的方法计算高分辨率城市卫星立体图像中同名点之间的匹配代价COST的方法具体为：

步骤3.1：计算Patch匹配代价COST_P；

ρ(q,q′)＝(1-ω)·min(||I_q-I_q′||₁,τ_gray)+ω·min(||▽I_q-▽I_q′||₁,τ_grad)

其中，W_p表示以像素点p为中心的方形窗口，S_W是该窗口的尺寸；κ是用于平衡灰度距离与坐标距离之间的影响比重；||I_p-I_q||₁表示点p灰度与点q灰度的L1距离；||C_p-C_q||₂表示点p坐标与点q坐标的L2距离；ω是用于平衡灰度与梯度之间的比重；||I_q-I_q′||₁表示点q与点q′处的灰度L1距离，||▽I_q-▽I_q′||₁表示点q与点q′处的梯度L1距离；q′为q在另一视图中的匹配点；参数τ_gray和参数τ_grad表示对于灰度距离与梯度距离的截断值，对于遮挡区域具有鲁棒性；

步骤3.2：计算Census匹配代价COST_C；

Census匹配代价是通过计算像素点之间bit字符串的汉明距离得到的，令左视图和右视图中支撑窗口中心像素点的Census变换值分别为Census(x_l,y_l)和Census(x_r,y_r)，则两像素点之间的Census代价为：

COST_C＝Dist_HAMMING(Census(x_l,y_l),Census(x_r,y_r))

其中，Dist_HAMMING(str1,str2)运算表示两个bit字符串str1和str2中对应位置不同bit值的个数之和，汉明距离越小说明两支撑窗口中心像素点的相似程度越高；基于左右视图的Census变换为：

其中，I_p(·)为待匹配像素点p(x,y)的灰度值；m×n为以p(x,y)为中心的矩形窗口的尺寸；

表示对经过δ(·)运算后所得的结果的前后连接运算；

步骤3.3：计算PHOW匹配代价COST_PH；

将像素点上不同半径圆形支撑块内的SIFT描述子作为该像素的PHOW特征，采用单一半径的SIFT特征，每个点由一个128维SIFT描述符表示，如下所示

PHOW(x,y)＝SIFT₈(x,y)

通过计算对应点PHOW描述符之间的欧几里得距离，作为对应点之间的PHOW匹配代价值COST_PH，计算公式如下：

COST_PH＝||PHOW(x_l,y_l)-PHOW(x_r,y_r)||₂

步骤3.4：将三种代价联合加权作为最终的高分辨率城市卫星立体图像中同名点之间的匹配代价COST；

COST＝α·COST_P+β·COST_C+γ·COST_PH

其中，α、β和γ分别表示代价权重。

Images