CN103247075A - 基于变分机制的室内环境三维重建方法 - Google Patents

Abstract

本发明属于计算机视觉与智能机器人的交叉领域，公开了一种基于变分机制的大范围室内场景的重建方法，包括：步骤一，获取相机的标定参数，并建立畸变矫正模型；步骤二，建立相机位姿描述与相机投影模型；步骤三，利用基于SFM的单目SLAM算法实现相机位姿估计；步骤四，建立基于变分机制的深度图估计模型，并求解该模型；步骤五，建立关键帧选取机制，实现三维场景的更新。本发明采用RGB相机获取环境数据，针对利用高精度单目定位算法，提出了一种基于变分机制的深度图生成方法，实现了大范围的快速室内三维场景重建，有效地解决了三维重建算法成本与实时性问题。

Description

基于变分机制的室内环境三维重建方法

技术领域

本发明属于计算机视觉与智能机器人的交叉领域，涉及一种室内环境三维重建技术，尤其涉及一种基于变分机制的大范围室内场景的重建方法。

技术背景

随着同时定位与地图创建（Simultaneous Localization And Mapping，SLAM）研究的不断深入，环境三维立体化建模已逐步成为该领域研究热点，引起众多学者的关注。G.Klein等于2007年在增强现实（AR）领域首先提出同时定位与地图创建（Parallel Tracking and Mapping,PTAM）的概念，以解决环境实时建模问题。PTAM将摄像机跟踪与地图生成划分为两个独立线程，利用FastCorner方法更新检测特征点的同时，采用最优的局部与全局光束平差法(Bundle Adjustment,BA)，不断实现相机位姿与三维特征点地图的更新。该方法基于稀疏点云建立了环境三维地图，但该地图缺乏对环境的直观三维描述。Pollefeys等人通过多传感器融合实现了大型室外场景的三维重建。但该方法存在计算的高复杂性以及对噪音敏感等缺点。目前在实时跟踪和稠密环境模型重构方面也有了一些尝试性的进步，但是仅仅局限于一些简单物体的重构，并且只能在特定约束条件下可以获得较高的精度。Richard A.Newcombe等人，利用基于SFM（Structure from Moving）的SLAM算法获取空间稀疏特征点云，采用多尺度径向基插值，运用图形图像学中隐式曲面多边形化方法，构造三维空间初始化网格地图，并结合场景流约束与高精度TV-L1光流算法更新网格顶点坐标，以达到逼近真实场景的目的。该算法能获取高精度的环境模型，但由于其算法复杂度较高，在两个图形硬件处理器（GPU）加速情况下，处理一帧图像仍需花费几秒钟的时间。

发明内容

针对现有技术中存在的上述问题，本发明提供了一种基于变分机制的快速三维重建方法，以实现在室内复杂环境下的三维建模。该方法保证环境信息的同时降低了所需处理数据量，能实现大范围的快速室内三维场景重建。有效地解决了三维重建算法成本与实时性问题，提高了重建精度。

本发明采用的技术方案如下：

利用PTAM算法作为相机位姿估计手段，并在关键帧处选取适当图像序列构造基于变分模式的深度图估计能量函数，运用原始对偶算法优化上述能量函数，实现在当前关键帧处环境深度图的获取。由于该算法利用邻近帧信息构造能量函数，且有效利用了特定视角坐标***间的关联性，以及摄像机透视投影变换关系，使得数据项蕴含了多视成像约束，降低了算法模型求解的计算复杂度。在统一计算架构下，本发明利用图形加速硬件实现了算法的并行优化，有效提高了算法实时性。

一种基于变分机制的室内环境三维重建的方法，其特征在于包括以下步骤：

步骤一，获取相机的标定参数，并建立畸变矫正模型。

在计算机视觉应用中，通过相机成像的几何模型，有效建立图像中像素点与空间三维点之间的映射关系。构成相机模型的几何参数须通过实验与计算才能得到，求解上述参数的过程就称之为相机标定。在本发明中相机参数的标定是非常关键的环节，标定参数的精度直接影响最终结果三维地图的准确性。

相机标定的具体过程为：

（1）打印一张棋盘模板。本发明采用一张A4纸，棋盘的间隔为0.25cm。

（2）从多个角度拍摄棋盘。拍摄时，应尽量让棋盘占满屏幕，并保证棋盘的每一个角都在屏幕中，一共拍摄6张模板图片。

（3）检测出图像中的特征点，即棋盘的每一个黑色交叉点。

（4）求取相机的内部参数，方法如下：

RGB相机标定参数主要为相机内参。相机的内参矩阵K为：

$K=\left[\begin{array}{ccc}{f}_u& 0& u_0\\ 0& f_v& v_0\\ 0& 0& 1\end{array}\right]$

式中，u、v为相机平面坐标轴，(u₀,v₀)是相机像平面中心坐标，(f_u,f_v)是相机的焦距。

根据标定参数，RGB图像中点与三维空间点的映射关系如下：RGB图像中点p＝(u,v)在相机坐标系下的坐标P_3D＝(x,y,z)表示为：

$\left\{\begin{array}{c}x=(u-u_0)*z/f_u\\ y=(v-v_0)*z/f_v\\ z=d\end{array}\right.$

式中，d表示深度图像中点p的深度值。

本发明中相机坐标系如图2所示，向下为y轴正方向，向前为z轴正方向，向右为x正方向。将相机的起始点位置设定为世界坐标系原点，世界坐标系的X、Y、Z方向与相机的定义相同。

FOV（Field of Viewer）相机矫正模型为：

$u_d=\left[\begin{array}{c}{u}_0\\ v_0\end{array}\right]+\left[\begin{array}{cc}{f}_u& 0\\ 0& f_v\end{array}\right]\frac{r_d}{r_u}{x}_u$

$r_d=\frac{1}{\mathrm{\ω}}\arctan \left({2r}_u\tan \frac{\mathrm{\ω}}{2}\right)$

$r_u=\frac{\tan \left(r_d\mathrm{\ω}\right)}{2\tan \frac{\mathrm{\ω}}{2}}$

式中，x_u为z=1面的像素坐标，u_d为原始图像中像素坐标，ω为FOV相机畸变系数。

步骤二，建立相机位姿描述与相机投影模型。

在已建立起的世界坐标系下，相机位姿可以表示为如下矩阵：

$T_{\mathrm{cw}}=\left[\begin{array}{cc}{R}_{\mathrm{cw}}& t_{\mathrm{cw}}\\ 0& 1\end{array}\right]$

式中，“cw”表示从世界坐标系到当前相机坐标系，T_cw∈SE(3)，SE(3)为刚体的旋转平移变换空间。T_cw可由如下六元组μ＝(μ₁,μ₂,μ₃,μ₄,μ₅,μ₆)表示，即：

$T_{\mathrm{cw}}=\exp \left(\widehat{\mathrm{\μ}}\right)$

$\hat{u}=\left[\begin{array}{cccc}0& {\mathrm{\μ}}_6& -{\mathrm{\μ}}_5& {\mathrm{\μ}}_1\\ {\mathrm{\μ}}_6& 0& {\mathrm{\μ}}_4& {\mathrm{\μ}}_2\\ {\mathrm{\μ}}_5& -{\mathrm{\μ}}_4& 0& {\mathrm{\μ}}_3\\ 0& 0& 0& 0\end{array}\right]$

式中，μ₁,μ₂,μ₃分别为Kinect在全局坐标系下的平移量，μ₄,μ₅,μ₆表示局部坐标系下坐标轴的旋转量。

相机的位姿T_cw建立了当前坐标系下空间点云坐标p_c到世界坐标p_w的变换关系，即：

p_c＝T_cwp_w

在当前标系下，三维空间点云到z=1平面上投影定义为：

π(p)＝(xz,yz)^T

式中，p∈R³为三维空间点，x，y，z为该点的坐标值。根据当前坐标点深度值d，利用逆向投影法确定当前空间三维点坐标p，其坐标关系可表示为：

π^-1(u,d)＝dK^-1u

步骤三，利用基于SFM的单目SLAM算法实现相机位姿估计。

目前，单目视觉SLAM算法主要包括基于滤波与SFM(Structure from Moving)的SLAM算法。本发明采用PTAM算法实现对相机的定位。该算法是一种基于SFM的单目视觉SLAM方法，通过将***划分为相机跟踪与地图创建两个独立的线程。在相机跟踪线程，***利用相机获取当前环境纹理信息，并构建四层高斯图像金字塔，运用FAST-10角点检测算法提取当前图像中特征信息，采用块匹配的方式建立角点特征间的数据关联。在此基础上，根据当前投影误差，建立位姿估计模型实现相机的精确定位，并结合特征匹配信息与三角测量算法生成当前三维点云地图。相机位姿估计的具体过程为：

（1）稀疏地图的初始化

PTAM算法利用标准立体相机算法模型建立当前环境初始化地图，并在此基础上结合新增加关键帧不断更新三维地图。在地图的初始化过程中，通过人为选择两个独立关键帧，利用图像中FAST角点匹配关系，采用基于随机采样一致性(Random Sample Consensus,RANSAC)的五点法实现上述关键帧间重要矩阵F估计，并计算当前特征点处的三维坐标，同时，结合RANSAC算法选取适当空间点建立当前一致性平面，以确定全局世界坐标系，实现地图的初始化。

（2）相机位姿估计

***利用相机获取当前环境纹理信息，并构建四层高斯图像金字塔，运用FAST-10角点检测算法提取当前图像中特征信息，采用块匹配的方式建立角点特征间数据关联。在此基础上，根据当前投影误差，建立位姿估计模型，其数学描述如下：

$\mathrm{\ξ}=\underset{\mathrm{\ξ}}{\arg \min }\mathrm{\ΣObj}(\frac{\left|e_j\right|}{{\mathrm{\σ}}_j},{\mathrm{\σ}}_T)$

$e_j=\left(\begin{array}{c}{u}_i\\ v_i\end{array}\right)-\mathrm{K\π}\left(\exp \left(\widehat{\mathrm{\ξ}}\right)p\right)$

式中，e_j是投影误差，∑Obj(·,σ_T)为Tukey双权目标函数，σ_T为特征点的匹配标准差的无偏估计值，ξ为当前位姿6元组表示，

为由ξ组成的反对称矩阵。

根据上述位姿估计模型，选取位于图像金字塔顶层的50个特征匹配点，实现对相机的初始化位姿估计。更进一步，该算法结合相机初始位姿，采用极线收索的方式，建立图像金字塔中角点特征亚像素精度匹配关系，并将上述匹配对带入位姿估计模型，实现相机的精确重定位。

（3）相机位姿优化

***经初始化后，地图创建线程将等待新的关键帧进入。若相机与当前关键帧间图像帧数超出阈值条件，且相机跟踪效果最佳时，将自动执行添加关键帧过程。此时，***将会对新增加关键帧中所有FAST角点进行Shi-Tomas评估，以获取当前具有显著特征的角点信息，并选取与之最近的关键帧利用极线收索与块匹配方法建立特征点映射关系，结合位姿估计模型实现相机精确重定位，同时将匹配点投影到空间，生成当前全局环境三维地图。

为了实现全局地图的维护，在地图创建线程等待新关键帧进入的过程中，***将利用局部与全局的Levenberg-Marquardt集束调整算法实现当前地图的一致性优化。该集束调整算法的数学描述为：

$\left\{\right\{={\mathrm{\ξ}}_1...{\mathrm{\ξ}}_N\},\{p_1...p_M\left\}\right\}=\underset{\left\{\right\{\mathrm{\μ}\},\{p\left\}\right\}}{\arg \min }\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j\∈s_i}{\mathrm{\Σ}}\mathrm{Obj}(\frac{\left|e_{\mathrm{ji}}\right|}{{\mathrm{\σ}}_{\mathrm{ji}}},{\mathrm{\σ}}_T)$

式中，σ_ji为在第i个关键帧中，FAST特征点的匹配标准差的无偏估计，ξ_i表示第i个关键帧位姿的6元组表示，p_i为全局地图中的点。

步骤四，建立基于变分机制的深度图估计模型，并求解该模型。

在PTAM精确位姿估计前提下，本发明基于多视重建方法，利用变分机制建立深度求解模型。该方法基于光照不变性与深度图平滑性假设，建立L1型数据惩罚项与变分规则项，该模型由在光照不变性假设的前提下建立数据惩罚项，并利用数据惩罚项保证当前深度图的平滑性，其数学模型如下：

E_d＝∫_Ω(E_data+λE_reg)dx

式中，λ为数据惩罚项E_data与变分规则项E_reg间的权重系数，

为深度图取值范围。

通过选取当前关键帧为深度图估计算法的参考帧I_r，利用其相邻像序列I＝{I₁,I₂,...,I_n}，结合投影模型建立数据惩罚项E_data，其数学描述为：

$E_{\mathrm{data}}=\frac{1}{\left|I\left(r\right)\right|}\underset{I_i\∈I}{\mathrm{\Σ}}|I_r\left(x\right)-I_i\left(x^{\′}\right)|$

式中，|I(r)|为当前图像序列中与参考帧具有重合信息的图像帧数量，x′为在当深度d下参考帧x处在I_i处的投影坐标，即：

$x^{\′}={\mathrm{\π}}^{-1}\left({\mathrm{KT}}_r^i\mathrm{\π}(x,d)\right)$

在深度图平滑性假设前提下，为了确保在图像中边界处的不连续性，引入加权Huber算子构建变分规则项，其数学描述为：

E_reg＝g(u)||▽d(u)||_α

式中，▽d为深度图的梯度，g(u)为像素梯度权重系数，且g(u)＝exp(-a||▽I_r(u)||)

Huber算子||x||_α的数学描述为：

${\left|\right|x\left|\right|}_{\mathrm{\α}}=\left\{\begin{array}{c}\frac{{\left|\right|x\left|\right|}^2}{2\mathrm{\α}},\left|\right|x\left|\right|\≤\mathrm{\α}\\ \left|\right|x\left|\right|-\frac{\mathrm{\α}}{2},\mathrm{others}\end{array}\right.$

式中，α为常量。

根据Legendre-Fenchel变换，能量函数可表示为： $g{\left|\right|\&dtri;d\left|\right|}_{\mathrm{\&alpha;}}=<g\&dtri;d,q>-\mathrm{\&delta;}\left(q\right)-\frac{\mathrm{\&alpha;}}{2}{\left|\right|q\left|\right|}^2$

式中， $\mathrm{\&delta;}\left(q\right)=\left\{\begin{array}{cc}\frac{\mathrm{\&alpha;}}{2}& \mathrm{\&alpha;}<\left|\right|q\left|\right|\&le;1\\ \&infin;& \mathrm{others}\end{array}\right.$

上述Huber算子的引入为三维重建过程提供了光滑性保证，同时也为确保深度图中存在非连续边界，提高了三维地图创建质量。

针对上述数学模型求解复杂度高、计算量大的问题，引入辅助变量建立凸优化模型，采用交替下降法实现对上述模型的优化，其具体过程如下：

(1)固定h，求解：

$\underset{q}{\arg \max }\left\{\underset{d}{\arg \min }{E}_{d,q}\right\}$

$E_{d,q}={\&Integral;}_{\mathrm{\&Omega;}}(<g\&dtri;d,q>+\frac{1}{2\mathrm{\&theta;}}{(d-h)}^2-\mathrm{\&delta;}\left(q\right)-\frac{\mathrm{\&alpha;}}{2}{\left|\right|q\left|\right|}^2)\mathrm{dx}$

式中，θ为二次项常系数，g为变分规则项中梯度权重系数。

根据拉格朗日极值法，上述能量函数达到极值的条件为：

$\frac{{\&PartialD;E}_{d,q}}{\&PartialD;q}=g\&dtri;d-\mathrm{\&alpha;q}=0$

$\frac{{\&PartialD;E}_{d,q}}{\&PartialD;d}=g\mathrm{div}q+\frac{1}{\mathrm{\&theta;}}(d-h)=0$

式中，divq为q的散度。

结合偏导数离散化描述，上述极值条件可表示为：

$\frac{q^{n+1}-q^n}{{\mathrm{\&epsiv;}}_q}=g\&dtri;d-{\mathrm{\&alpha;q}}^{n+1}$

$\frac{d^{n+1}-d^n}{{\mathrm{\&epsiv;}}_d}=g\mathrm{div}p+\frac{1}{\mathrm{\&theta;}}(d^{n+1}-h)$

此时可采用原始对偶算法实现能量函数的迭代优化，即：

$p^{n+1}=\frac{(p^n+{\mathrm{\&epsiv;}}_{q}g\&dtri;d^n)/(1+{\mathrm{\&epsiv;}}_q\mathrm{\&alpha;})}{\max (1,(p^n+{\mathrm{\&epsiv;}}_{q}g\&dtri;d^n)/(1+{\mathrm{\&epsiv;}}_q\mathrm{\&alpha;}))}$

$d^{n+1}=\frac{d^n+{\mathrm{\&epsiv;}}_d({g\mathrm{div}q}^{n+1}+h^n/\mathrm{\&theta;})}{(1+{\mathrm{\&epsiv;}}_d/\mathrm{\&theta;})}$

式中，ε_q、ε_d为常数，分别表示最大化与最小化梯度描述系数。

(2)固定d，求解：

$\underset{h}{\arg \min }{E}_h$

$E_h={\&Integral;}_{\mathrm{\&Omega;}}(\frac{\mathrm{\&theta;}}{2}{(d-h)}^2+\frac{\mathrm{\&lambda;}}{\left|I\left(r\right)\right|}\underset{i=0}{\overset{n}{\mathrm{\&Sigma;}}}|I_i\left(x\right)-I_{\mathrm{ref}}(x,h)\left|\right)\mathrm{dx}$

在上述能量函数求解过程中，为了有效减少算法的复杂度，同时保证重建过程中的部分细节信息。本发明将深度取值范围[d_min,d_max]划分为S个采样平面，采用穷举的方式获取当前能量函数的最优解。其中步长的选择为：

$d_{\mathrm{inc}}^k=\frac{{\mathrm{Sd}}_{\min }{d}_{\max }}{(S-k)d_{\min }+d_{\max }}$

式中，为k与k-1采样平面间隔。

步骤五，建立关键帧选取机制，实现三维场景的更新。

考虑***冗余信息的消除，为了提高重建结果的清晰度以及实时性，减少***在计算负担，本发明只在关键帧处实现对三维场景的估计，并更新和维护所生成的三维场景。当新增一帧KeyFrame数据后，根据式

将当前新增KeyFrame数据转换到世界坐标系中，完成场景数据的更新。

利用深度模型中数据惩罚项，建立当前帧与关键帧间的信息重合程度评估函数，即：

$N=\underset{x\&Element;R^2}{\mathrm{\&Sigma;}}c\left(x\right)$

$c\left(x\right)=\left\{\begin{array}{cc}1,& |I_r\left(x\right)-I_i\left(x^{\&prime;}\right)|<\mathrm{\&zeta;}\\ 0,& \mathrm{others}\end{array}\right.$

式中，ζ为常数。

若此时N小于图像大小的0.7时，即确定当前帧为新关键帧。

本发明的有益效果是：本发明采用RGB相机获取环境数据。针对利用高精度单目定位算法，提出一种基于变分机制的深度图生成方法，实现了大范围的快速室内三维场景重建，有效地解决了三维重建算法成本与实时性问题。

附图说明

图1为基于变分模型的室内三维场景重建方法流程图；

图2为相机坐标系示意图；

图3为本发明应用实例的三维重建实验结果。

具体实施方式

图1是基于变分模型的室内三维场景重建方法流程图，包括以下步骤：

步骤一，获取相机的标定参数，并建立畸变矫正模型。

步骤二，建立相机位姿描述与相机投影模型。

步骤三，利用基于SFM的单目SLAM算法实现相机位姿估计。

步骤四，建立基于变分机制的深度图估计模型，并求解该模型。

步骤五，建立关键帧选取机制，实现三维场景的更新。

下面给出本发明的一个应用实例。

本实例采用的RGB相机为Point Grey Flea2，图像辨率为640×480，最高帧频为30fps，水平视场角为65°，焦距大约为3.5mm。所使用的PC机配备有GTS450GPU和i5四核CPU。

在实验过程中，通过彩色相机获取环境深度信息，结合相机位姿估计算法实现对自身精确定位。当进入关键帧后，选择关键帧周围20帧图像作为本文深度估计算法的输入。在深度估计算法执行过程中，令d⁰＝h⁰且q⁰＝0，计算

以获取当前深度图的初始化输入，并迭代优化E_d,q与E_h直到收敛。同时，该在算法迭代过程中不断减小θ值，增加二次函数在算法执行过程中的权重，有效提高了算法收敛速度。最终实验结果如图3所示，实验表明该方法能有效实现环境的稠密三维重建，并进一步验证了该方法的可行性。

Claims (3)

1.一种基于变分机制的室内环境三维重建的方法，其特征在于包括以下步骤：

步骤一，获取相机的标定参数，并建立畸变矫正模型；

相机标定的具体过程为：

（1）打印一张棋盘模板；

（2）从多个角度拍摄棋盘，应尽量让棋盘占满屏幕，并保证棋盘的每一个角都在屏幕中，一共拍摄6张模板图片；

（3）检测出图像中的特征点，即棋盘的每一个黑色交叉点；

（4）求取的内部参数，方法如下：

RGB相机标定参数主要为相机内参，相机的内参矩阵K为：

$K=\left[\begin{array}{ccc}{f}_u& 0& u_0\\ 0& f_v& v_0\\ 0& 0& 1\end{array}\right]$

式中，u、v为相机平面坐标轴，(u₀,v₀)是相机像平面中心坐标，(f_u,f_v)是相机的焦距；

根据标定参数，RGB图像中点与三维空间点的映射关系如下：RGB图像中点p＝(u,v)在相机坐标系下的坐标P_3D＝(x,y,z)表示为：

$\left\{\begin{array}{c}x=(u-u_0)*z/f_u\\ y=(v-v_0)*z/f_v\\ z=d\end{array}\right.$

式中，d表示深度图像中点p的深度值；

相机坐标系向下为y轴正方向，向前为z轴正方向，向右为x正方向；将相机的起始点位置设定为世界坐标系原点，世界坐标系的X、Y、Z方向与相机的定义相同；

FOV相机矫正模型为：

$u_d=\left[\begin{array}{c}{u}_0\\ v_0\end{array}\right]+\left[\begin{array}{cc}{f}_u& 0\\ 0& f_v\end{array}\right]\frac{r_d}{r_u}{x}_u$

$r_d=\frac{1}{\mathrm{\&omega;}}\arctan \left({2r}_u\tan \frac{\mathrm{\&omega;}}{2}\right)$

$r_u=\frac{\tan \left(r_d\mathrm{\&omega;}\right)}{2\tan \frac{\mathrm{\&omega;}}{2}}$

式中，x_u为z=1面的像素坐标，u_d为原始图像中像素坐标，ω为FOV相机畸变系数；

步骤二，建立相机位姿描述与相机投影模型，方向如下：

在已建立的世界坐标系下，相机位姿可以表示为如下矩阵：

$T_{\mathrm{cw}}=\left[\begin{array}{cc}{R}_{\mathrm{cw}}& t_{\mathrm{cw}}\\ 0& 1\end{array}\right]$

式中，cw表示从世界坐标系到当前相机坐标系，T_cw∈SE(3)，SE(3)为刚体的旋转平移变换空间；T_cw可由如下六元组μ＝(μ₁,μ₂,μ₃,μ₄,μ₅,μ₆)表示，即：

$T_{\mathrm{cw}}=\exp \left(\widehat{\mathrm{\&mu;}}\right)$

$\widehat{\mathrm{\&mu;}}=\left[\begin{array}{cccc}0& {\mathrm{\&mu;}}_6& -{\mathrm{\&mu;}}_5& {\mathrm{\&mu;}}_1\\ {\mathrm{\&mu;}}_6& 0& {\mathrm{\&mu;}}_4& {\mathrm{\&mu;}}_2\\ {\mathrm{\&mu;}}_5& -{\mathrm{\&mu;}}_4& 0& {\mathrm{\&mu;}}_3\\ 0& 0& 0& 0\end{array}\right]$

式中，μ₁,μ₂,μ₃分别为Kinect在全局坐标系下的平移量，μ₄,μ₅,μ₆表示局部坐标系下坐标轴的旋转量；

相机的位姿T_cw建立了当前坐标系下空间点云坐标p_c到世界坐标p_w的变换关系，即：

p_c＝T_cwp_w

在当前标系下，三维空间点云到z=1平面上投影定义为：

π(p)＝(xz,yz)^T

式中，p∈R³为三维空间点，x，y，z为该点的坐标值；根据当前坐标点深度值d，利用逆向投影法确定当前空间三维点坐标p，其坐标关系可表示为：

π^-1(u,d)＝dK^-1u

步骤三，利用基于SFM的单目SLAM算法实现相机位姿估计；

步骤四，建立基于变分机制的深度图估计模型，并求解该模型；

步骤五，建立关键帧选取机制，实现三维场景的更新，方法如下：

在关键帧处实现对三维场景的估计，并更新和维护所生成的三维场景；当新增一帧KeyFrame数据后，根据式

将当前新增KeyFrame数据转换到世界坐标系中，完成场景数据的更新；

利用深度模型中数据惩罚项，建立当前帧与关键帧间的信息重合程度评估函数，即：

$N=\underset{x\&Element;R^2}{\mathrm{\&Sigma;}}c\left(x\right)$

$c\left(x\right)=\left[\begin{array}{cc}1,& |I_r\left(x\right)-I_i\left(x^{\&prime;}\right)|<\mathrm{\&zeta;}\\ 0,& \mathrm{others}\end{array}\right]$

式中，ζ为常数；

若此时N小于图像大小的0.7时，即确定当前帧为新关键帧。

2.根据权利要求1所述的一种基于变分机制的室内环境三维重建的方法，其特征在于，步骤三利用基于SFM的单目SLAM算法实现相机位姿估计的方法还包括以下步骤：

（1）稀疏地图的初始化

PTAM算法利用标准立体相机算法模型建立当前环境初始化地图，并在此基础上结合新增加关键帧不断更新三维地图；在地图的初始化过程中，通过人为选择两个独立关键帧，利用图像中FAST角点匹配关系，采用基于随机采样一致性的五点法实现上述关键帧间的重要矩阵F的估计，并计算当前特征点处的三维坐标，同时，结合RANSAC算法选取适当空间点建立当前一致性平面，以确定全局世界坐标系，实现地图的初始化；

（2）相机位姿估计

***利用相机获取当前环境纹理信息，并构建四层高斯图像金字塔，运用FAST-10角点检测算法提取当前图像中特征信息，采用块匹配的方式建立角点特征间的数据关联；在此基础上，根据当前投影误差，建立位姿估计模型，其数学描述如下：

$\mathrm{\&xi;}=\underset{\mathrm{\&xi;}}{\arg \min }\mathrm{\&Sigma;Obj}(\frac{\left|e_j\right|}{{\mathrm{\&sigma;}}_j},{\mathrm{\&sigma;}}_T)$

$e_j=\left(\begin{array}{c}{u}_i\\ v_i\end{array}\right)-\mathrm{K\&pi;}\left(\exp \left(\mathrm{\&xi;}\right)p\right)$

式中，e_j是投影误差，ΣObj(·,σ_T)为Tukey双权目标函数函数，σ_T为特征点的匹配标准差的无偏估计值，ξ为当前位姿6元组表示，

为由ξ组成的反对称矩阵；

根据上述位姿估计模型，选取位于图像金字塔顶层的50个特征匹配点，实现对相机的初始化位姿估计；更进一步，该算法结合相机初始位姿，采用极线收索的方式，建立图像金字塔中角点特征亚像素精度匹配关系，并将上述匹配对带入位姿估计模型，实现相机的精确重定位；

（3）相机位姿优化

***经初始化后，地图创建线程等待新的关键帧进入；若相机与当前关键帧间图像帧数超出阈值条件，且相机跟踪效果最佳时，将自动执行添加关键帧过程；此时，***将会对新增加的关键帧中所有FAST角点进行Shi-Tomas评估，以获取当前具有显著特征的角点特征信息，并选取与之最近的关键帧利用极线收索与块匹配方法建立特征点映射关系，结合位姿估计模型实现相机精确重定位，同时将匹配点投影到空间，生成当前全局环境三维地图；

为了实现全局地图的维护，在地图创建线程等待新的关键帧进入的过程中，***利用局部与全局的Levenberg-Marquardt集束调整算法实现当前地图的全局一致性优化；该集束调整算法的数学描述为：

$\left\{\right\{{\mathrm{\&xi;}}_2..{\mathrm{\&xi;}}_N\},\{p_1..p_M\left\}\right\}=\underset{\left\{\right\{\mathrm{\&mu;}\},\{p\left\}\right\}}{\arg \min }\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j\&Element;S_i}{\mathrm{\&Sigma;}}\mathrm{Obj}(\frac{\left|e_{\mathrm{ji}}\right|}{{\mathrm{\&sigma;}}_{\mathrm{ji}}},{\mathrm{\&sigma;}}_T)$

式中，σ_ji为在第i个关键帧中，FAST特征点的匹配标准差的无偏估计，ξ_i表示第i个关键帧位姿的6元组表示，p_i为全局地图中的点。

3.根据权利要求1所述的一种基于变分机制的室内环境三维重建的方法，其特征在于，步骤四建立并求解基于变分机制的深度图估计模型的方法如下：

基于变分机制的深度图估计模型，在光照不变性假设的前提下建立数据惩罚项，并利用数据惩罚项保证当前深度图的平滑性，其数学模型如下：

E_d＝∫_Ω(E_data+λE_reg)dx

式中，λ为数据惩罚项E_data与变分规则项E_reg间的权重系数，

为深度图取值范围；

通过选取当前关键帧为深度图估计算法的参考帧I_r，利用其相邻像序列I＝{I₁,I₂,...,I_n}，结合投影模型建立数据惩罚项E_data，其数学描述为：

$E_{\mathrm{data}}=\frac{1}{\left|I\left(r\right)\right|}\underset{I_i\&Element;I}{\mathrm{\&Sigma;}}|I_r\left(x\right)-I_i\left(x^{\&prime;}\right)|$

式中，|I(r)|为当前图像序列中与参考帧具有重合信息的图像帧数量，x′为在当深度d下参考帧x处在I_i处的投影坐标，即：

$x^{\&prime;}={\mathrm{\&pi;}}^{-1}\left({\mathrm{KT}}_r^i\mathrm{\&pi;}(x,d)\right)$

在深度图平滑性假设前提下，为了确保在图像中边界处的不连续性，引入加权Huber算子构建变分规则项，其数学描述为：

E_reg＝g(u)||▽d(u)||_α

式中，▽d为深度图的梯度，g(u)为像素梯度权重系数，g(u)＝exp(-a||▽I_r(u)||)

Huber算子||x||_α的数学描述为：

${\left|\right|x\left|\right|}_{\mathrm{\&alpha;}}=\left\{\begin{array}{c}\frac{{\left|\right|x\left|\right|}^2}{2\mathrm{\&alpha;}},\left|\right|x\left|\right|\&le;\mathrm{\&alpha;}\\ \left|\right|x\left|\right|-\frac{\mathrm{\&alpha;}}{2},\mathrm{others}\end{array}\right.$

式中，α为常量；

根据Legendre-Fenchel变换，能量函数变换为：

$g{\left|\right|\&dtri;d\left|\right|}_{\mathrm{\&alpha;}}=<g\&dtri;d,q>-\mathrm{\&delta;}\left(q\right)-\frac{\mathrm{\&alpha;}}{2}{\left|\right|q\left|\right|}^2$

式中， $\mathrm{\&delta;}\left(q\right)=\left\{\begin{array}{cc}\frac{\mathrm{\&alpha;}}{2}& \mathrm{\&alpha;}<\left|\right|q\left|\right|\&le;1\\ \&infin;& \mathrm{other}\end{array}\right.$

鉴于上述数学模型求解复杂度高、计算量大，引入辅助变量建立凸优化模型，采用交替下降法实现对上述模型的优化，具体过程如下：

(1)固定h，求解：

$\underset{q}{\arg \max }\left\{\underset{d}{\arg \min }{E}_{d,q}\right\}$

$E_{d,q}={\&Integral;}_{\mathrm{\&Omega;}}(<g\&dtri;d,q>+\frac{1}{2\mathrm{\&theta;}}{(d-h)}^2-\mathrm{\&delta;}\left(q\right)-\frac{\mathrm{\&alpha;}}{2}{\left|\right|q\left|\right|}^2)\mathrm{dx}$

式中，g为变分规则项中梯度权重系数，θ为二次项常系数；

根据拉格朗日极值法，上述能量函数达到极值的条件为：

$\frac{\&PartialD;E_{d,q}}{\&PartialD;q}=g\&dtri;d-\mathrm{\&alpha;q}=0$

$\frac{{\&PartialD;E}_{d,q}}{\&PartialD;d}=g\mathrm{div}q+\frac{1}{\mathrm{\&theta;}}(d-h)=0$

式中，divq为q的散度；

结合偏导数离散化描述，上述极值条件可表示为：

$\frac{q^{n+1}-q^n}{{\mathrm{\&epsiv;}}_q}=g\&dtri;d-{\mathrm{\&alpha;q}}^{n+1}$

$\frac{d^{n+1}-d^n}{{\mathrm{\&epsiv;}}_d}=g\mathrm{div}p+\frac{1}{\mathrm{\&theta;}}(d^{n+1}-h)$

采用原始对偶算法实现能量函数的迭代优化，即：

$p^{n+1}=\frac{(p^n+{\mathrm{\&epsiv;}}_{q}g\&dtri;d^n)/(1+{\mathrm{\&epsiv;}}_q\mathrm{\&alpha;})}{\max (1,(p^n+{\mathrm{\&epsiv;}}_{q}g\&dtri;d^n)/(1+{\mathrm{\&epsiv;}}_q\mathrm{\&alpha;}))}$

$d^{n+1}=\frac{d^n+{\mathrm{\&epsiv;}}_d({g\mathrm{div}q}^{n+1}+h^n/\mathrm{\&theta;})}{(1+{\mathrm{\&epsiv;}}_d/\mathrm{\&theta;})}$

式中，ε_q、ε_d为常数，分别表示最大化与最小化梯度描述系数；

(2)固定d，求解：

$\underset{h}{\arg \min E_h}$

$E_h={\&Integral;}_{\mathrm{\&Omega;}}(\frac{\mathrm{\&theta;}}{2}{(d-h)}^2+\frac{\mathrm{\&lambda;}}{\left|I\left(r\right)\right|}\underset{i=0}{\overset{n}{\mathrm{\&Sigma;}}}|I_i\left(x\right)-I_{\mathrm{ref}}(x,h)|)\mathrm{dx}$

在上述能量函数求解过程中，为了有效减少算法的复杂度，同时保证重建过程中的部分细节信息，将深度取值范围[d_min,d_max]划分为S个采样平面，采用穷举的方式获取当前能量函数的最优解；其中步长的选择为：

$d_{\mathrm{inc}}^k=\frac{{\mathrm{Sd}}_{\min }{d}_{\max }}{(S-k)d_{\min }+d_{\max }}$

式中，

为k与k-1采样平面间隔。

Images