基于高阶扩展FMM和低秩矩阵恢复的信号重建研究.pdf
裘国永+张力
摘 要：由于在成像过程中出现遮挡现象，图像矩阵的元素有缺失。在正投影相机模型下，提出一种基于列约束的低秩矩阵恢复方法。该方法利用图像矩阵是一个低秩矩阵从而图像序列具有冗余性的特性，利用奇异值分解由图像矩阵的列空间构造出一个投影矩阵，得到图像矩阵的列所满足的约束条件，将缺失元素的恢复转化为迭代求解二次型的极值问题，利用它恢复出图像矩阵的缺失元素。该方法从理论上能够保证收敛到全局最小值。仿真实验表明，此方法具有收敛速度快，恢复精度高等优点。
关键词：遮挡点；低秩矩阵；奇异值分解；列约束
中图分类号：TP391 文献标志码：A
Low-rank Matrix Recovery Method Based on Column Metric Constraints
QIU Guo-yong，ZHANG Li
（School of Computer Science，Shaanxi Normal University，Xian，Shaanxi 710119，China）
Abstract：To recover the position of occlusion，the image matrix has some missing elements because of occluding in the imaging process.An low-rank matrix recovery method based on column metric constraints under orthographic projection is presented.Utilizing the property that the image matrix is of low rank and the redundancy that exists in an image sequence，a projective matrix is constructed via singular value decomposition to get the image matrixs column metric constraints.The method iteratively solves the minimum of a quadratic function to recover the missing elements of the image matrix.The method can guarantee to converge to the global optimal solution theoretically.The simulated experiments show that the proposed method has the advantages of fast convergence speed and small error.
Key words：occlusion；low-rank matrix；singular value decomposition；column metric constraints
1 引 言
特征点跟踪是计算机视觉研究领域中三维重建、运动分析等深度研究的前提和基础[1]。由于在摄像机成像过程中视场改变、遮挡出现、光照变化等诸多原因，在特征点的跟踪过程中，必然会导致部分特征点跟踪丢失，从而图像矩阵的相应位置没有值。如何恢复这些缺失的矩阵元素，是计算机视觉的热门研究问题之一[2-3]。
由于图像矩阵构造的特殊性，图像矩阵在没有任何误差的情形下是一个秩不超过3的低秩矩阵[4]，这也表明图像序列具有冗余性，且图像矩阵的列向量组是线性相关的，列空间可以用不同的列向量部分组生成。利用图像矩阵的低秩性，有的学者利用子矩阵法对低秩矩阵进行恢复[4-5]，但是该方法的结果和所选择的子矩阵有关。为了克服该方法的缺点，不少研究人员利用图像矩阵可以分解为两个矩阵相乘的特性，其中一个表示相机的旋转，一个表示目标物的形状，交替固定一个求另外一个，进行循环迭代求解[6-9]。但是该类方法适用于遮挡率比较低的情形，当遮挡率较高时，算法就无法收敛。同时，也有研究人员通过重投影误差构造出目标函数，采用非线性优化
方法进行迭代求解[10-11]，但由于该目标函数是一个非凸函数，导致若初始值远离真实解时，不能求到全局最优解。为了解决这个问题，在假设相机为针孔模型下，Martinec等利用三线性约束来求解遮挡点的真实位置[12]，但该方法的缺点是要求三线性约束关系，而在求三线性约束关系过程中，不能把所有的图像平等看待，因此必然导致误差较大。
为了克服上述缺点，本文提出了一种基于列约束的低秩矩阵恢复方法，利用图像矩阵是一个低秩矩阵从而图像序列具有冗余性，通过奇异值分解由图像矩阵的列空间构造出一个投影矩阵，由此投影矩阵得到图像矩阵的列所满足的约束条件，根据这些约束条件将遮挡点的求解转化为迭代求解一个半正定二次型的极值问题，利用它恢复出图像矩阵的缺失元素。由于半正定二次型对应优化问题的目标函数是凸的、一步就可以求解出函数的极值等，因此本文方法具有良好的數值求解特性。本文方法从理论上证明了算法的收敛性。为了验证方法的有效性，本文进行了仿真实验，与Martinec方法进行了对比，结果表明，本文算法具有收敛速度快，恢复精度高等优点。
2 相机模型和图像矩阵
为了表示图像矩阵，首先需要确定相机模型。在最常用的针孔模型中，对应的透视投影是非线性的，导致计算量大且在透视效果不明显时表现出病态性，即对躁声敏感。在某些实际情形，比如摄像机的视场小，且目标物体与摄像机有一定的距离，则可以用线性相机模型—正投影相机模型近似代替针孔相机模型，这样不仅可以简化公式的表示而且能减少计算量。因此本文取相机模型为正投影模型。endprint
设给定的图像序列中有f幅二维图像及s个三维空间特征点，则有
mi，j=Rixj（1）
式中mi，j=ui，jvi，jT表示二维图像点，xj=xjyjzjT表示三维空间点，Pi为相机的投影矩阵，i和j分别表示第i幅图像和第j个图像点。
为了得到图像矩阵，现按分块矩阵的方法将所有特征点的图像坐标合并成一个2f×s图像矩阵，每幅图像对应的投影矩陣也合并成一个2f×3矩阵。则由式（1）可得到
M2f×s=R2f×3S3×s （2）
式中M2f×s=u1，1u1，2…u1，sv1，1v1，2…v1，suf，1uf，2…uf，svf，1vf，2…vf，s2f×s，R2f×3=R1R2Rf2f×3，S3×s=x1x2…xs。
由此可见，图像矩阵可以表示成两个矩阵相乘，其中一个与相机的旋转有关，另一个与目标物的形状有关。
由于R只有3列，S只有3行，所以根据线性数相关结论得到图像矩阵M2f×s的秩为3。且M的列向量组和R的列向量组张成的是R2f的同一个三维子空间。
3 基于列约束的低秩矩阵恢复方法
由于图像矩阵M2f×s的秩为3，所以它的列空间是R2f的三维子空间。因此先找出它的一组基底。为此，对M2f×s进行奇异值分解（SVD），得到
M2f×s=S′2f×2fV′2f×sD′s×s=S2f×3V3×3D3×s（3）
其中S2f×3为S′2f×2f的前3列组成的矩阵，D3×s为D′s×s的前3行组成的矩阵，V3×3为V′2f×s的左上角3×3子矩阵，即3阶对角阵。
从上式可以看出，M2f×s中的列向量组和S2f×3的列向量组张成R2f的同一三维线性子空间。同时，R2f到S2f×3列向量组张成的子空间的投影矩阵为：
P2f×2f=S2f×3ST2f×3S2f×3-1ST2f×3（4）
R2f到其正交补空间上的投影矩阵为： P⊥2f×2f=I2f-S2f×3ST2f×3S2f×3-1ST2f×3 （5）
式中I2f为单位阵。
由于S2f×3的列向量组两两正交，且模都为1，因此S2f×3满足
ST2f×3S2f×3=I3×3 （6）
从而式（5）可以简化为
P⊥2f×2f=I2f×2f-S2f×3ST2f×3（7）
为了方便描述，任取M2f×s中的一个列向量ci。由于M2f×s中的列向量和S2f×3的列向量张成R2f的同一线性子空间，因此，ci投影到S2f×3列向量生成的子空间的正交补空间上的投影为02m×1，即
P⊥2f×2fci=02f×1（8）
这就是M2f×s的元素要满足的列约束条件。但通常情况下，躁声使图像矩阵的元素不可能准确满足（11）。因此，可以将上式的求解转化为求余差ei的极小值，即
ei=P⊥2f×2fciTP⊥2f×2fci（9）
由于P⊥2f×2f为投影矩阵，因此上式可以化简成
ei=cTiP⊥2f×2fci（10）
为了表示方便，令
P⊥2f×2f=p1，1p1，2…p1，2fp2，1p2，2…p2，2fp2f，1p2f，2…p2f，2f（11）
同时，假设由于遮挡的原因，ci含有两个缺失的元素，即
ci=c1c2…c-k…c-l…c2fT（12）
式中c-表示为遮挡点。
将式（11）和（12）代入式（10）中，则有
ei=ck—cl—pk，kpk，lpl，kpl，lck—cl—+∑2fm=1，m≠k，lpk，mcm+∑2fm=1，m≠k，lpm，kcm∑2fm=1，m≠k，lpl，mcm+∑2fm=1，m≠k，lpm，lcmTck—cl—+c′TiPc′2f-2×2f-2c′i
即
ei=xTiAixi+2bTixi+εi（13）
式中xi=c-kc-lT为未知向量，Ai=pk，kpk，lpl，kpl，l为矩阵P⊥2f×2f的第k、l行和第k、l列构成的子矩阵，
bi=12∑2fm=1，m≠k，lpk，mcm+∑2fm=1，m≠k，lpm，kcm∑2fm=1，m≠k，lpl，mcm+∑2fm=1，m≠k，lpm，lcm
=∑2fm=1，m≠k，lpk，mcm∑2fm=1，m≠k，lpl，mcm，εi=c′iTPc′（2m-2）×（2m-2）c′i为常数项，c′i为向量ci中将两个缺失元素去除后得到的向量，Pc′2f-2×2f-2为矩阵P⊥2f×2f中删除了第k、l行和第k、l列的矩阵。
式（13）是一个半正定二次型，对应的极小值问题的目标函数是一个凸函数，因此它的最优解可以直接求出，即为
xi=A-1ibi（14）
在上述推导过程中，我们假设ci仅含有两个缺失元素，若含有n个缺失元素，则Ai为一个n×n的P⊥2f×2f子矩阵，其元素取自P⊥2f×2f对应于未知元素所在的行和列；bi为一个n×1的向量，其元素为ci中已知元素对应的P⊥2f×2f行中的元素和ci中已知元素相乘之和。
同时，在上述推导中，我们仅针对M2f×s矩阵中的第i列ci，对于所有的其他列，求解都是一样的，则总余差为：
e=∑si=1ei（15）
在上面的推导过程中，我们假设投影矩阵P⊥2f×2f已知，而事实上P⊥2f×2f由M2f×s矩阵的奇异值分解得到。但在我们所研究的问题中，由于出现了遮挡现象，低秩矩阵M2f×s中产生了缺失元素，因此，无法直接进行奇异值分解。因此算法开始时，可以假设缺失元素的初值由其他已知元素的值的均值确定，然后利用列约束条件，求出缺失元素，此时求出的缺失元素比开始时给定的初值要更加接近真实值，由此可以构造一个迭代算法，恢复出低秩图像矩阵的缺失元素。endprint
4 算法总结与计算复杂度
现给出基于列约束的低秩矩阵恢复算法：
（1）假设所有的遮挡点m-i，j=1k∑kj=1mi，j（即第i幅图像中所有已知分量的均值），迭代次数l=1，令ε为任意小的一个正数；
（2）利用式（3），对M2f×s进行奇异值分解；
（3）利用式（7）求投影矩阵P⊥2f×2f；
（4）对于图像矩阵的每一列，利用式 （14）求出该列的缺失元素xli，并将求出的元素xli代替M2f×s中的缺失元素；
（5）若‖xli-xl-1i‖
SymbolcB@ ε则停止；否则l=l+1，转第（2）步。
由于式（13） 是一个半正定二次型，从而有唯一的极小值点，而且极小值可以通过式 （14）直接求出。又由于它是一个凸函数，所以迭代误差满足
∑si=1e（1）i≥∑si=1e（2）i≥…≥∑si=1e（l）i （15）
因此，理论上保证算法能够收敛到全局最小值。
本算法的每一步迭代需要进行一次奇异值分解，以及计算P⊥2f×2f，Ai，bi，xi。由于M2f×s奇异值分解的时间复杂度是O（f2s），P⊥2f×2f的计算复杂度是O（f2）。另外由于Ai是P⊥2f×2f的二阶子矩阵，bi是二维向量，其计算复杂度为O（f）。而xi=A-1ibi的计算复杂度是一个常数，与f和s均无关。所以本算法的计算复杂度为O（f2sl），其中l为算法收敛所需要的迭代次数。
5 仿真实验
为了检验本文所提方法的收敛速度和恢复效果，现在对含有大量遮挡点的图像序列进行仿真实验。实验的实验对象是仿真图像。本文首先用计算机模拟的方法得到100个三维空间点作为特征点，即S为3×100矩阵。另外为了得到矩阵R，模拟相机的位置变化，得到一个包括50幅大小640×480的图像序列，即R为100×100矩阵。随机选择50%的特征点进行遮挡，并在图像中依次加上均值都为0，方差为0.5、1、1.5、2、2.5、3个像素的高斯噪声。利用这些模拟得到的图像点用本文算法进行图像矩阵的恢复，算法的迭代收敛图如图1所示。
同时，为了研究本文算法的收敛性与图像点的遮挡率之间的关系，在图像中加入2.0个像素的高斯噪声，并分别遮挡20%，30%，40%，50%及60%的图像点，实验结果如图2所示。
从图1和图2可以看出，本文所提出的方法迭代20次以内就能够达到收敛，具有良好的收斂性能，原因是由于本文算法将遮挡点的求解转化为二次型的迭代求解，而二次型的极值点可以一步求解。同时，从图2还可以看出，遮挡率越高，算法收敛速度越慢，原因是由于遮挡率越高，图像序列的冗余性也少，约束越少，因此求解的未知数就越多，因此计算量就更大，算法的收敛速度就越慢。
为了比较本文算法和Martinec方法的收敛速度和恢复误差，用上述方法产生50幅图像，遮挡率取为20%，并在图像中依次加入均值为0，方差从0变化到3的高斯噪声，分别用本文算法和Martinec方法进行图像矩阵恢复。在每种噪声下各重复100次实验，然后求取平均误差，实验结果如图3所示。
从图3可以看出，相比Martinec方法本文算法具有更高的恢复精度，原因主要有三方面，一方面是由于本文算法不需要计算基础矩阵，基础矩阵的计算对图像噪声非常敏感；其二方面是由于本文算法将遮挡点的求解转换为二次型的迭代求解，从而对应优化问题的目标函数是一个凸函数，它可以一步求解到全局最优值，而且算法能够保证收敛；其三方面是由于本文算法把所有的图像平等地看待，并没有倚重某些图像。因此，本文算法具有较小的误差。
6 结束语
本文在正投影相机模型下，提出一种基于列约束的低秩矩阵恢复方法。该方法利用图像矩阵是一个低秩矩阵从而图像序列具有冗余性，利用奇异值分解由图像矩阵的列空间构造出一个投影矩阵，将遮挡点的求解转化为迭代求解二次型的极值问题，利用它恢复出遮挡点，补上图像矩阵的缺失元素。由于该方法的关键在于迭代求解二次型的极值问题，而二次型是一个凸函数，它可以一步求解到全局最优值，因此具有良好的收敛性能。仿真实验结果表明，本文方法相比Martinec方法具有收敛速度快，恢复精度高等优点。由于在此算法中我们只考虑了基于列向量的约束，进一步的研究方向将是把基于行向量和列向量的约束同时考虑，以提高算法的效率。
参考文献
[1] KENNEDY R，BALZANO L，WRIGHT S，et al.Online algorithms for factorization-based structure from motion[J].Computer Vision and Image Understanding，2016，150： 139-152.
[2] KIM H，LEUTENEGGER S，DAVISON A.Real-Time 3D Reconstruction and 6-DoF Tracking with an Event Camera[C].14th European Conference on Computer Vision，11-14 October 2016，Amsterdam，Netherlands，349-364.
[3] COHEN A，SCHONBERGER J，SPECIALE P，SATTLER T，FRAHM J，POLLEFEYS M.Indoor-Outdoor 3D Reconstruction Alignment[C].14th European Conference on Computer Vision，11-14 October 2016，Amsterdam，Netherlands，285-300.endprint
[4] TOMASI C，KANADE T.Shape and Motion from Image Streams under Orthography： a Factorization Method[J].International Journal of Computer Vision，1992，9（2）：137-154.
[5] JACOBS D.Linear fitting with missing data for structure-from-motion[J].Computer Vision and Image Understanding，2008，82（1）：57-81.
[6] STRELOW D.General and nested Wiberg minimization[C].IEEE Conference on Computer Vision and Pattern Recognition，Rhode Island，16-21 June 2012，1584-1591.
[7] OKATANI T，DEGUCHI K.On the Wiberg algorithm for matrix factorization in the presence of missing components[J].International Journal of Computer Vision，2007，72（3）：329-337
[8] ZHENG Y，SUGIMOTO S，YAN S，et al.Generalizing Wiberg algorithm for rigid and nonrigid factorizations with missing components and metric constraints[C].IEEE Conference on Computer Vision and Pattern Recognition，16-21 June 2012，Rhode Island，2010-2017.
[9] STRELOW D，WANG Q，SI L，et al.General，Nested，and Constrained Wiberg Minimization[J].IEEE Transactions on Pattern Analysis and Machine Intelligence.2016，38（9） ：1803-1815.
[10] ALBL C，SUGIMOTO A，PAJDLA T.Degeneracies in Rolling Shutter SfM[C].14th European Conference on Computer Vision，11-14 October 2016，Amsterdam，Netherlands，36-51.
[11] AGUDO A，NOGUER F，CALVO B，et al.Sequential Non-Rigid Structure from Motion Using Physical Priors[J].IEEE Transactions on Pattern Analysis and Machine Intelligence，2016，38（5）： 979-994.
[12] MARTINEC D，PAJDLA T.Outlier detection for factorization-based reconstruction from perspective images with occlusions[A].In Proceedings of the Photogrammetric Computer Vision[C].TU-Graz，2002，B：161-164.endprint
