CN106960450A - 基于块的影像匹配数字表面模型的全局高程优化方法 - Google Patents

Abstract

本发明涉及基于块的影像匹配数字表面模型的全局高程优化方法，步骤为：1、输入原始的IMDSM；2、将输入的原始IMDSM分割成一系列规则的块；3、构建全局能量函数中的数据项；4、构建全局能量函数中的平滑项；5、根据数据项和平滑项，建立全局能量函数，并解算函数；6、执行羽化算法，消除块之间的拼接缝；本发明能够有效消除目前IMDSM中普遍存在的“高程阶梯”问题，优化后的IMDSM表面连续光滑，三维可视效果好。本发明将IMDSM的优化问题转化为矩阵函数的极值计算问题，无需迭代，可直接计算全局最优解，时间复杂度低，适合大范围DSM的全局优化。

Description

基于块的影像匹配数字表面模型的全局高程优化方法

技术领域

本发明涉及影像匹配数字表面模型的高程优化技术领域，具体地指一种基于块的影像匹配数字表面模型的全局高程优化方法。

技术背景

影像匹配数字表面模型(Image Matching Based Digital Surface Model,IMDSM)指通过影像密集匹配技术，获取大范围自然地表和人工地物表面的数字三维模型。生产IMDSM所需的影像包括卫星影像、航拍影像、无人机影像等等。IMDSM生成技术是计算机视觉、数字摄影测量等领域的核心技术之一，在数字地球、导航定位、巡航避障、文物重建、抢险救灾、灾后重建、空间数据库更新等应用中，发挥着重要的作用。

传统的IMDSM生产流程，往往采用局部匹配算法、半全局匹配算法、基于图像引导的非局部匹配算法、全局匹配算法等简单高效的一维编号(1D label)匹配算法。1D label匹配算法总是假设影像中相邻像素之间要尽可能满足视差一致的约束。实际上，在地表斜面区域，相邻像素的视差是不可能一致的。因此，虽然1D label算法的时间复杂度低，但是1D label算法生成的IMDSM在斜面区域会存在明显的“高程阶梯”问题，影响IMDSM的高程精度和三维可视效果。

发明内容

本发明的目的就是要提供一种基于块的影像匹配数字表面模型的全局高程优化方法，该方法将初始的影像匹配数字表面模型分割成一系列相互邻接的块，采用全局优化的方法和羽化算法，获取表面连续光滑的影像匹配数字表面模型。本发明能够有效解决目前影像匹配数字表面模型中普遍存在的“高程阶梯”问题，使得影像匹配数字表面模型表面连续光滑，增强影像匹配数字表面模型的三维可视效果。

为实现此目的，本发明所设计的一种基于块的影像匹配数字表面模型的全局高程优化方法，其特征在于，它包括如下步骤：

步骤1：将原始影像匹配数字表面模型输入到计算机中；

步骤2：计算机采用SLIC超像素分割的方法，将输入的原始影像匹配数字表面模型分割成一系列规则的块，用S_i表示第i个块；

步骤3：构建全局能量函数中的数据项，用E_data表示数据项，保证优化后的影像匹配数字表面模型的高程，与原始影像匹配数字表面模型的高程保持一致；

每个所述块都用一个高程平面方程来描述，即：

d(t_i)＝a_i·t_ix+b_i·t_iy+c_i；t_i∈S_i

其中，a_i、b_i、c_i表示块S_i的高程平面方程参数；t_i＝(t_ix,t_iy)^T表示块S_i内的一个三维点的平面坐标；d表示三维点的高程坐标；

令表示块S_i所在平面的高程平面方程系数，表示优化影像匹配数字表面模型的全局能量函数的变量，q表示块的数目，i∈1,...q，将数据项E_data表示为矩阵计算的形式，即：

式中，G_data表示数据项E_data的二次项系数矩阵；H_data表示数据项E_data的一次项的系数矩阵，E_data为优化影像匹配数字表面模型的全局能量函数的数据项；l_data表示数据项E_data的常数项，上述系数矩阵和常数项具体表达为：

其中，g_i、h_i和l_i分别表示G_data、H_data和l_data中的块矩阵；t_ix、t_iy表示块S_i内每个点的平面坐标；|S_i|表示块S_i内三维点的数目；Diag表示对角线矩阵，t_i表示原始影像匹配数字表面模型上的块S_i内的一个点，d₀(t_i)表示点t_i在原始影像匹配数字表面模型上的高程；

步骤4：构建优化影像匹配数字表面模型的全局能量函数中的平滑项，用E_smooth表示优化影像匹配数字表面模型的全局能量函数的平滑项，保证优化后的三维点云表面分片光滑连续；同样令表示优化影像匹配数字表面模型的全局能量函数的变量，将平滑项E_smooth表示为矩阵计算的形式：

式中，G_s表示平滑项的二次项系数矩阵，G_s可以表达为：

其中，q表示块的数目；S_j表示S_i的邻接块；N(S_i)表示块S_i的相邻块集合；E(S_i,S_j)表示块S_i内，与块S_j相邻接的像素，c_i＝(c_ix,c_iy)^T表示块S_i的重心；P(i,j)表示根据块S_i与块S_j的邻接关系计算出来的惩罚系数，t表示原始影像匹配数字表面模型上的一个点，该点位于集合E(S_i,S_j)内；g_sr(i,j,t)表示块S_i与块S_j之间的相关矩阵。

其中，g_sr代表块S_i与块S_j之间的邻接矩阵；σ₁表示分块矩阵，T为转置符号；0_3×3表示3×3的零矩阵；

步骤5：根据数据项和平滑项，构建优化影像匹配数字表面模型的全局能量函数，其中，优化影像匹配数字表面模型的全局能量函数的极值解，即为影像匹配数字表面模型的全局优化结果；

定义D表示影像匹配数字表面模型，E(D)表示优化影像匹配数字表面模型的全局能量函数，那么，将优化影像匹配数字表面模型的全局能量函数定义为：

其中，E_data表示优化影像匹配数字表面模型的全局能量函数的数据项，数据项E_data保证优化的影像匹配数字表面模型的高程与原始影像匹配数字表面模型的高程尽可能的保持一致；E_smooth表示优化影像匹配数字表面模型的全局能量函数的平滑项，平滑项E_smooth保证优化后的三维点云表面分片光滑连续；

求公式的最小值，等价于求采用Cholesky矩阵分解方法直接计算获得全局最优的影像匹配数字表面模型；

步骤6：采用羽化算法，修正全局最优的影像匹配数字表面模型中块与块之间交界处的高程，消除全局最优的影像匹配数字表面模型中块与块之间的拼接缝，使得全局最优的影像匹配数字表面模型中块与块之间平滑过渡，获得平滑的影像匹配数字表面模型。

本发明的有益效果为：

本发明能够有效消除目前影像匹配数字表面模型中普遍存在的高程阶梯问题，采用本发明的方法优化后的影像匹配数字表面模型表面连续光滑，三维可视效果好。本发明将影像匹配数字表面模型的优化问题转化为矩阵函数的极值计算问题，无需迭代，可直接计算全局最优解，时间复杂度低，适合大范围DSM(Digital Surface Model，数字表面模型)的全局优化，能够为数字摄影测量、计算机视觉、数字地球、巡航避障、虚拟显示等学科和应用提供技术服务。

附图说明

图1为“高程阶梯”问题示意图；

图2为本发明的流程图；

图3为SLIC超像素分割示意图；

图4为全局优化后的影像匹配数字表面模型示意图；

图5为羽化示意图。

具体实施方式

以下结合附图和具体实施例对本发明作进一步的详细说明：

本发明针对目前实际工程应用中生产的影像匹配数字表面模型，普遍存在的高程阶梯问题，如图1所示，提出一种基于块的影像匹配数字表面模型全局高程优化方法。该方法无需迭代，能够直接快速获得优化后的影像匹配数字表面模型，能够有效解决高程阶梯问题，使得优化后的影像匹配数字表面模型表面连续光滑，三维可视效果好。图1中(a)表示原始影像，图1中(b)表示根据影像匹配出来的影像匹配数字表面模型，可以比较分析发现，匹配出来的原始影像匹配数字表面模型，存在明显的高程阶梯问题。本发明的工作流程如图2所示，它包括如下步骤：

步骤1：将原始影像匹配数字表面模型输入到计算机中；

原始影像匹配数字表面模型即根据多视影像直接匹配出来的数字表面模型，影像可以采用卫星影像、航拍影像、无人机影像等，匹配的算法可以采用局部匹配算法、基于图像引导的非局部匹配算法、半全局匹配算法，全局匹配算法等；

步骤2：计算机采用SLIC(simple linear iterative cluster)超像素分割的方法，将输入的原始影像匹配数字表面模型分割成一系列规则的块，用S_i表示第i个块，认为每个块在三维空间近似为一个平面，可以采用高程平面方程来描述，超像素分割代码参见网址http://ivrl.epfl.ch/supplementary_material/RK_SLICSuperpixels/index.html；

步骤3：构建全局能量函数中的数据项，用E_data表示数据项，保证优化的影像匹配数字表面模型的高程，与原始影像匹配数字表面模型的高程保持一致；

每个所述块都用一个高程平面方程来描述，即：

d(t_i)＝a_i·t_ix+b_i·t_iy+c_i；t_i∈S_i

其中，a_i、b_i、c_i表示块S_i的高程平面方程参数；t_i＝(t_ix,t_iy)^T表示块S_i内的一个三维点的平面坐标；d表示三维点的高程坐标；

令表示块S_i所在平面的高程平面方程系数，表示优化影像匹配数字表面模型的全局能量函数的变量，q表示块的数目，i∈1,...q，将数据项E_data表示为矩阵计算的形式，即：

式中，G_data表示数据项E_data的二次项系数矩阵；H_data表示数据项E_data的一次项的系数矩阵，E_data为优化影像匹配数字表面模型的全局能量函数的数据项；l_data表示数据项E_data的常数项，上述系数矩阵和常数项具体表达为：

其中，g_i、h_i和l_i分别表示G_data、H_data和l_data中的块矩阵；t_ix、t_iy表示块S_i内每个点的平面坐标；|S_i|表示块S_i内三维点的数目；Diag表示对角线矩阵，t_i表示原始影像匹配数字表面模型上块S_i内的的一个点，d₀(t_i)表示点t_i在原始影像匹配数字表面模型上的高程，T为转置符号；

步骤4：构建优化影像匹配数字表面模型的全局能量函数中的平滑项，用E_smooth表示优化影像匹配数字表面模型的全局能量函数的平滑项，保证优化后的三维点云表面分片光滑连续；同样令表示优化影像匹配数字表面模型的全局能量函数的变量，将平滑项E_smooth表示为矩阵计算的形式：

式中，G_s表示平滑项的二次项系数矩阵，G_s可以表达为：

其中，q表示块的数目；S_j表示S_i的邻接块；N(S_i)表示块S_i的相邻块集合；E(S_i,S_j)表示块S_i内，与块S_j相邻接的像素，c_i＝(c_ix,c_iy)^T表示块S_i的重心；P(i,j)表示根据块S_i与块S_j的邻接关系计算出来的惩罚系数，t表示原始影像匹配数字表面模型上的一个点，该点位于集合E(S_i,S_j)内；g_sr(i,j,t)表示块S_i与块S_j之间的相关矩阵；

其中，g_sr代表块S_i与块S_j之间的邻接矩阵；σ₁表示分块矩阵；0_3×3表示3×3的零矩阵；

步骤5：根据数据项和平滑项，构建优化影像匹配数字表面模型的全局能量函数，其中，优化影像匹配数字表面模型的全局能量函数的极值解，即为影像匹配数字表面模型的全局优化结果；

定义D表示影像匹配数字表面模型，E(D)表示优化影像匹配数字表面模型的全局能量函数，那么，将优化影像匹配数字表面模型的全局能量函数定义为：

其中，E_data表示优化影像匹配数字表面模型的全局能量函数的数据项，数据项E_data保证优化后的影像匹配数字表面模型的高程与原始影像匹配数字表面模型的高程尽可能的保持一致；E_smooth表示优化影像匹配数字表面模型的全局能量函数的平滑项，平滑项E_smooth保证优化后的三维点云表面分片光滑连续；

求公式的最小值，等价于求采用Cholesky矩阵分解方法直接计算具体解算线性方程组的代码可以参见开源的eigen库：http://eigen.tuxfamily.org/index.php？title＝Main_Page，获得全局最优的影像匹配数字表面模型，解决了“高程阶梯”问题，如图4所示；

步骤6：采用羽化算法，修正全局最优的影像匹配数字表面模型中块与块之间交界处的高程，消除全局最优的影像匹配数字表面模型中块与块之间的拼接缝，使得全局最优的影像匹配数字表面模型中块与块之间平滑过渡，获得平滑的影像匹配数字表面模型。

上述技术方案中，所述步骤4中，根据块S_i与块S_j的邻接关系计算出来的惩罚系数P(i,j)有：

其中，P为人工给定的惩罚值阈值(根据经验，可设为16)，exp表示以自然指数e为底的指数函数；分别表示在块S_i与块S_j之间的邻接像素中平均的高程值；num(i,j)表示块S_i与块S_j之间的邻接像素数目；σ_d是优化影像匹配数字表面模型的高程平滑因子(一般设为2)、σ_n为优化影像匹配数字表面模型的连接性平滑因子(一般设为5)。

上述技术方案中，所述步骤5中，解算后，即获得每个块的平面方程参数(a_i,b_i，c_i)^T，然后根据每个块内的点云平面坐标和平面方程参数，计算每个点优化后的高程，如下式所示：

d(t_i)＝a_i·t_ix+b_i·t_iy+c_i；t_i∈S_i.

其中，a_i、b_i、c_i表示块S_i的高程平面方程参数；t_i＝(t_ix,t_iy)^T表示块S_i内的一个三维点的平面坐标；d表示三维点优化后的高程坐标，依次计算每个块内所有三维点的高程坐标，获得全局最优的影像匹配数字表面模型。

上述技术方案中，步骤6：采用羽化算法修正全局最优的影像匹配数字表面模型中块与块之间交界处的高程(在全局优化后的影像匹配数字表面模型的曲面区域中，块与块之间会存在明显的拼接缝，影响三维重建的可视效果)，首先，以全局最优的影像匹配数字表面模型的块与块之间的拼接缝为中心，定义一个半径为l的缓冲区，如图5中(a)所示，缓冲区内的所有点参与羽化计算，像素p为缓冲区内的一个像素，距离拼接缝的距离为l'，点p的第一高程由块S_i的平面方程参数确定，设为d_i；点p的第二高程由块S_j的平面方程参数确定，设为d_j.点p羽化后的深度，由d_i和d_j联合确定，羽化后的结果如图5中(b)所示，如下式所示：

d′(p)＝w·d_i+(1-w)·d_j.

式中，d′(p)表示像素p羽化后的深度；w表示权值，根据像素p距离拼接缝的距离为l'来确定：w＝0.5+l′/_2l°

本发明能够有效解决目前影像匹配数字表面模型中普遍存在的“高程”阶梯问题，所获取的影像匹配数字表面模型表面连续光滑，三维可视效果好。本发明所设计的算法时间复杂度低，适合大范围影像匹配数字表面模型的高程优化。

本说明书未作详细描述的内容属于本领域专业技术人员公知的现有技术。

Claims (4)

1.一种基于块的影像匹配数字表面模型的全局高程优化方法，其特征在于，它包括如下步骤：

步骤1：将原始影像匹配数字表面模型输入到计算机中；

步骤2：计算机采用SLIC超像素分割的方法，将输入的原始影像匹配数字表面模型分割成一系列规则的块，用S_i表示第i个块；

步骤3：构建全局能量函数中的数据项，用E_data表示数据项，保证优化后的影像匹配数字表面模型的高程，与原始影像匹配数字表面模型的高程保持一致；

每个所述块都用一个高程平面方程来描述，即：

d(t_i)＝a_i·t_ix+b_i·t_iy+c_i；t_i∈S_i

其中，a_i、b_i、c_i表示块S_i的高程平面方程参数；t_i＝(t_ix,t_iy)^T表示块S_i内的一个三维点的平面坐标；d表示三维点的高程坐标；

令表示块S_i所在平面的高程平面方程系数，表示优化影像匹配数字表面模型的全局能量函数的变量，q表示块的数目，i∈1,...q，将数据项E_data表示为矩阵计算的形式，即：

$E_{data}={\tilde{x}}^T{G}_{data}\tilde{x}-2H_{data}^T\tilde{x}+l_{data}$

式中，G_data表示数据项E_data的二次项系数矩阵；H_data表示数据项E_data的一次项的系数矩阵，E_data为优化影像匹配数字表面模型的全局能量函数的数据项；l_data表示数据项E_data的常数项，T为转置符号，上述系数矩阵和常数项具体表达为：

$G_{data}=Diag\left(g_i\right);H_{data}={\left(\begin{array}{cccc}{h}_0^T& h_1^T& ...& h_q^T\end{array}\right)}^T;$

$l_{data}=\underset{i=1}{\overset{q}{\Σ}}{l}_i;$

$g_i=\left(\begin{array}{ccc}\underset{t_i\∈S_i}{\mathrm{\Σ}}{t_{ix}}^2& \underset{t_i\∈S_i}{\mathrm{\Σ}}{t}_{ix}\·t_{iy}& \underset{t_i\∈S_i}{\mathrm{\Σ}}{t}_{ix}\\ \underset{t_i\∈S_i}{\mathrm{\Σ}}{t}_{ix}\·t_{iy}& \underset{t_i\∈S_i}{\mathrm{\Σ}}{t_{iy}}^2& \underset{t_i\∈S_i}{\mathrm{\Σ}}{t}_{iy}\\ \underset{t_i\∈S_i}{\mathrm{\Σ}}{t}_{ix}& \underset{t_i\∈S_i}{\mathrm{\Σ}}{t}_{iy}& \left|S_i\right|\end{array}\right);$

$h_i={\left(\begin{array}{ccc}\underset{t_i\∈S_i}{\mathrm{\Σ}}{t}_{ix}\·d_0\left(t_i\right)& \underset{t_i\∈S_i}{\mathrm{\Σ}}{t}_{iy}\·d_0\left(t_i\right)& \underset{t_i\∈S_i}{\mathrm{\Σ}}{d}_0\left(t_i\right)\end{array}\right)}^T;$

$l_i=\underset{t_i\∈S_i}{\Σ}{d}_0{\left(t_i\right)}^2$

其中，g_i、h_i和l_i分别表示G_data、H_data和l_data中的块矩阵；t_ix、t_iy表示块S_i内每个点的平面坐标；|S_i|表示块S_i内三维点的数目；Diag表示对角线矩阵，t_i表示原始影像匹配数字表面模型上块S_i内的一个点，d₀(t_i)表示点t在原始影像匹配数字表面模型上的高程；

步骤4：构建优化影像匹配数字表面模型的全局能量函数中的平滑项，用E_smooth表示优化影像匹配数字表面模型的全局能量函数的平滑项，保证优化后的三维点云表面分片光滑连续；同样令表示优化影像匹配数字表面模型的全局能量函数的变量，将平滑项E_smooth表示为矩阵计算的形式：

$E_{smooth}={\tilde{x}}^T{G}_S\tilde{x}$

式中，G_s表示平滑项的二次项系数矩阵，G_s可以表达为：

$G_s=\underset{i=1}{\overset{q}{\Σ}}\left(\underset{S_j\∈N\left(S_i\right)}{\Σ}P\right(i,j)\·\underset{t\∈\left(E\right(S_i,S_j)\∪c_i)}{\Σ}{g}_{sr}\left(i,j,t\right))$

其中，q表示块的数目；S_j表示S_i的邻接块；N(S_i)表示块S_i的相邻块集合；E(S_i,S_j)表示块S_i内，与块S_j相邻接的像素，c_i＝(c_ix,c_iy)^T表示块S_i的重心；P(i,j)表示根据块S_i与块S_j的邻接关系计算出来的惩罚系数，t表示原始影像匹配数字表面模型上的一个点，该点位于集合E(S_i,S_j)内；g_sr(i,j,t)表示块S_i与块S_j之间的相关矩阵；

${\mathrm{\σ}}_1(i,j,t)=\left(\begin{array}{ccc}{t_{ix}}^2& t_{ix}\·t_{iy}& t_{ix}\\ t_{ix}\·t_{iy}& {t_{iy}}^2& t_{iy}\\ t_{ix}& t_{iy}& 1\end{array}\right).$

其中，g_sr代表块S_i与块S_j之间的邻接矩阵；σ₁表示分块矩阵；0_3×3表示3×3的零矩阵；

步骤5：根据数据项和平滑项，构建优化影像匹配数字表面模型的全局能量函数，其中，优化影像匹配数字表面模型的全局能量函数的极值解，即为影像匹配数字表面模型的全局优化结果；

定义D表示影像匹配数字表面模型，E(D)表示优化影像匹配数字表面模型的全局能量函数，那么，将优化影像匹配数字表面模型的全局能量函数定义为：

$\begin{array}{c}E\left(D\right)=E_{data}+E_{smooth}\\ ={\tilde{x}}^T\left(G_{data}+G_s\right)\tilde{x}-2H_{data}^T\tilde{x}+l_{data}\end{array}$

其中，E_data表示优化影像匹配数字表面模型的全局能量函数的数据项；E_smooth表示优化影像匹配数字表面模型的全局能量函数的平滑项，平滑项E_smooth保证优化后的三维点云表面分片光滑连续；

求公式的最小值，等价于求采用Cholesky矩阵分解方法直接计算获得全局最优的影像匹配数字表面模型；

步骤6：采用羽化算法，修正全局最优的影像匹配数字表面模型中块与块之间交界处的高程，消除全局最优的影像匹配数字表面模型中块与块之间的拼接缝，使得全局最优的影像匹配数字表面模型中块与块之间平滑过渡，获得平滑的影像匹配数字表面模型。

2.根据权利要求1所述的基于块的影像匹配数字表面模型的全局高程优化方法，其特征在于：所述根据块S_i与块S_j的邻接关系计算出来的惩罚系数P(i,j)有：

$P(i,j)=P\·\exp (-|\stackrel{\‾}{{\mathrm{nd}}_i}-\stackrel{\‾}{{\mathrm{nd}}_j}|/{\mathrm{\σ}}_d)\·(1-\exp (-num\left(i,j\right)/{\mathrm{\σ}}_n\left)\right);$

其中，P为人工给定的惩罚值阈值，exp表示以自然指数e为底的指数函数；分别表示在块S_i与块S_j之间的邻接像素中平均的高程值；num(i,j)表示块S_i与块S_j之间的邻接像素数目；σ_d是优化影像匹配数字表面模型的高程平滑因子、σ_n为优化影像匹配数字表面模型的连接性平滑因子。

3.根据权利要求1所述的基于块的影像匹配数字表面模型的全局高程优化方法，其特征在于：所述步骤5中，解算后，即获得每个块的平面方程参数(a_i,b_i,c_i)^T，然后根据每个块内的点云平面坐标和平面方程参数，计算每个点优化后的高程，如下式所示：

d(t_i)＝a_i·t_ix+b_i·t_iy+c_i；t_i∈S_i.

其中，a_i、b_i、c_i表示块S_i的高程平面方程参数；t_i＝(t_ix,t_iy)^T表示块S_i内的一个三维点的平面坐标；d表示三维点优化后的高程坐标，依次计算每个块内所有三维点的高程坐标，获得全局最优的影像匹配数字表面模型。

4.根据权利要求1所述的基于块的影像匹配数字表面模型的全局高程优化方法，其特征在于：步骤6：在采用羽化算法修正全局最优的影像匹配数字表面模型中，首先，以全局最优的影像匹配数字表面模型的块与块之间的拼接缝为中心，定义一个半径为l的缓冲区，缓冲区内的所有点参与羽化计算，像素p为缓冲区内的一个像素，距离拼接缝的距离为l'，点p的第一高程由块S_i的平面方程参数确定，设为d_i；点p的第二高程由块S_j的平面方程参数确定，设为d_j.点p羽化后的深度，由d_i和d_j联合确定，如下式所示：

d′(p)＝w·d_i+(1-w)·d_j.

式中，d′(p)表示像素p羽化后的深度；w表示权值，根据像素p距离拼接缝的距离为l'来确定：w＝0.5+l′/2l。