裘國(guó)永+張力

摘 要：由于在成像過(guò)程中出現(xiàn)遮擋現(xiàn)象，圖像矩陣的元素有缺失。在正投影相機(jī)模型下，提出一種基于列約束的低秩矩陣恢復(fù)方法。該方法利用圖像矩陣是一個(gè)低秩矩陣從而圖像序列具有冗余性的特性，利用奇異值分解由圖像矩陣的列空間構(gòu)造出一個(gè)投影矩陣，得到圖像矩陣的列所滿足的約束條件，將缺失元素的恢復(fù)轉(zhuǎn)化為迭代求解二次型的極值問(wèn)題，利用它恢復(fù)出圖像矩陣的缺失元素。該方法從理論上能夠保證收斂到全局最小值。仿真實(shí)驗(yàn)表明，此方法具有收斂速度快，恢復(fù)精度高等優(yōu)點(diǎn)。

關(guān)鍵詞：遮擋點(diǎn)；低秩矩陣；奇異值分解；列約束

中圖分類號(hào)：TP391 文獻(xiàn)標(biāo)志碼：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 引 言

特征點(diǎn)跟蹤是計(jì)算機(jī)視覺(jué)研究領(lǐng)域中三維重建、運(yùn)動(dòng)分析等深度研究的前提和基礎(chǔ)[1]。由于在攝像機(jī)成像過(guò)程中視場(chǎng)改變、遮擋出現(xiàn)、光照變化等諸多原因，在特征點(diǎn)的跟蹤過(guò)程中，必然會(huì)導(dǎo)致部分特征點(diǎn)跟蹤丟失，從而圖像矩陣的相應(yīng)位置沒(méi)有值。如何恢復(fù)這些缺失的矩陣元素，是計(jì)算機(jī)視覺(jué)的熱門(mén)研究問(wèn)題之一[2-3]。

由于圖像矩陣構(gòu)造的特殊性，圖像矩陣在沒(méi)有任何誤差的情形下是一個(gè)秩不超過(guò)3的低秩矩陣[4]，這也表明圖像序列具有冗余性，且圖像矩陣的列向量組是線性相關(guān)的，列空間可以用不同的列向量部分組生成。利用圖像矩陣的低秩性，有的學(xué)者利用子矩陣法對(duì)低秩矩陣進(jìn)行恢復(fù)[4-5]，但是該方法的結(jié)果和所選擇的子矩陣有關(guān)。為了克服該方法的缺點(diǎn)，不少研究人員利用圖像矩陣可以分解為兩個(gè)矩陣相乘的特性，其中一個(gè)表示相機(jī)的旋轉(zhuǎn)，一個(gè)表示目標(biāo)物的形狀，交替固定一個(gè)求另外一個(gè)，進(jìn)行循環(huán)迭代求解[6-9]。但是該類方法適用于遮擋率比較低的情形，當(dāng)遮擋率較高時(shí)，算法就無(wú)法收斂。同時(shí)，也有研究人員通過(guò)重投影誤差構(gòu)造出目標(biāo)函數(shù)，采用非線性優(yōu)化

方法進(jìn)行迭代求解[10-11]，但由于該目標(biāo)函數(shù)是一個(gè)非凸函數(shù)，導(dǎo)致若初始值遠(yuǎn)離真實(shí)解時(shí)，不能求到全局最優(yōu)解。為了解決這個(gè)問(wèn)題，在假設(shè)相機(jī)為針孔模型下，Martinec等利用三線性約束來(lái)求解遮擋點(diǎn)的真實(shí)位置[12]，但該方法的缺點(diǎn)是要求三線性約束關(guān)系，而在求三線性約束關(guān)系過(guò)程中，不能把所有的圖像平等看待，因此必然導(dǎo)致誤差較大。

為了克服上述缺點(diǎn)，本文提出了一種基于列約束的低秩矩陣恢復(fù)方法，利用圖像矩陣是一個(gè)低秩矩陣從而圖像序列具有冗余性，通過(guò)奇異值分解由圖像矩陣的列空間構(gòu)造出一個(gè)投影矩陣，由此投影矩陣得到圖像矩陣的列所滿足的約束條件，根據(jù)這些約束條件將遮擋點(diǎn)的求解轉(zhuǎn)化為迭代求解一個(gè)半正定二次型的極值問(wèn)題，利用它恢復(fù)出圖像矩陣的缺失元素。由于半正定二次型對(duì)應(yīng)優(yōu)化問(wèn)題的目標(biāo)函數(shù)是凸的、一步就可以求解出函數(shù)的極值等，因此本文方法具有良好的數(shù)值求解特性。本文方法從理論上證明了算法的收斂性。為了驗(yàn)證方法的有效性，本文進(jìn)行了仿真實(shí)驗(yàn)，與Martinec方法進(jìn)行了對(duì)比，結(jié)果表明，本文算法具有收斂速度快，恢復(fù)精度高等優(yōu)點(diǎn)。

2 相機(jī)模型和圖像矩陣

為了表示圖像矩陣，首先需要確定相機(jī)模型。在最常用的針孔模型中，對(duì)應(yīng)的透視投影是非線性的，導(dǎo)致計(jì)算量大且在透視效果不明顯時(shí)表現(xiàn)出病態(tài)性，即對(duì)躁聲敏感。在某些實(shí)際情形，比如攝像機(jī)的視場(chǎng)小，且目標(biāo)物體與攝像機(jī)有一定的距離，則可以用線性相機(jī)模型—正投影相機(jī)模型近似代替針孔相機(jī)模型，這樣不僅可以簡(jiǎn)化公式的表示而且能減少計(jì)算量。因此本文取相機(jī)模型為正投影模型。endprint

設(shè)給定的圖像序列中有f幅二維圖像及s個(gè)三維空間特征點(diǎn)，則有

mi，j=Rixj（1）

式中mi，j=ui，jvi，jT表示二維圖像點(diǎn)，xj=xjyjzjT表示三維空間點(diǎn)，Pi為相機(jī)的投影矩陣，i和j分別表示第i幅圖像和第j個(gè)圖像點(diǎn)。

為了得到圖像矩陣，現(xiàn)按分塊矩陣的方法將所有特征點(diǎn)的圖像坐標(biāo)合并成一個(gè)2f×s圖像矩陣，每幅圖像對(duì)應(yīng)的投影矩陣也合并成一個(gè)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。

由此可見(jiàn)，圖像矩陣可以表示成兩個(gè)矩陣相乘，其中一個(gè)與相機(jī)的旋轉(zhuǎn)有關(guān)，另一個(gè)與目標(biāo)物的形狀有關(guān)。

由于R只有3列，S只有3行，所以根據(jù)線性數(shù)相關(guān)結(jié)論得到圖像矩陣M2f×s的秩為3。且M的列向量組和R的列向量組張成的是R2f的同一個(gè)三維子空間。

3 基于列約束的低秩矩陣恢復(fù)方法

由于圖像矩陣M2f×s的秩為3，所以它的列空間是R2f的三維子空間。因此先找出它的一組基底。為此，對(duì)M2f×s進(jìn)行奇異值分解（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階對(duì)角陣。

從上式可以看出，M2f×s中的列向量組和S2f×3的列向量組張成R2f的同一三維線性子空間。同時(shí)，R2f到S2f×3列向量組張成的子空間的投影矩陣為：

P2f×2f=S2f×3ST2f×3S2f×3-1ST2f×3（4）

R2f到其正交補(bǔ)空間上的投影矩陣為： P⊥2f×2f=I2f-S2f×3ST2f×3S2f×3-1ST2f×3 （5）

式中I2f為單位陣。

由于S2f×3的列向量組兩兩正交，且模都為1，因此S2f×3滿足

ST2f×3S2f×3=I3×3 （6）

從而式（5）可以簡(jiǎn)化為

P⊥2f×2f=I2f×2f-S2f×3ST2f×3（7）

為了方便描述，任取M2f×s中的一個(gè)列向量ci。由于M2f×s中的列向量和S2f×3的列向量張成R2f的同一線性子空間，因此，ci投影到S2f×3列向量生成的子空間的正交補(bǔ)空間上的投影為02m×1，即

P⊥2f×2fci=02f×1（8）

這就是M2f×s的元素要滿足的列約束條件。但通常情況下，躁聲使圖像矩陣的元素不可能準(zhǔn)確滿足（11）。因此，可以將上式的求解轉(zhuǎn)化為求余差ei的極小值，即

ei=P⊥2f×2fciTP⊥2f×2fci（9）

由于P⊥2f×2f為投影矩陣，因此上式可以化簡(jiǎn)成

ei=cTiP⊥2f×2fci（10）

為了表示方便，令

P⊥2f×2f=p1，1p1，2…p1，2fp2，1p2，2…p2，2fp2f，1p2f，2…p2f，2f（11）

同時(shí)，假設(shè)由于遮擋的原因，ci含有兩個(gè)缺失的元素，即

ci=c1c2…c-k…c-l…c2fT（12）

式中c-表示為遮擋點(diǎn)。

將式（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列構(gòu)成的子矩陣，

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為常數(shù)項(xiàng)，c′i為向量ci中將兩個(gè)缺失元素去除后得到的向量，Pc′2f-2×2f-2為矩陣P⊥2f×2f中刪除了第k、l行和第k、l列的矩陣。

式（13）是一個(gè)半正定二次型，對(duì)應(yīng)的極小值問(wèn)題的目標(biāo)函數(shù)是一個(gè)凸函數(shù)，因此它的最優(yōu)解可以直接求出，即為

xi=A-1ibi（14）

在上述推導(dǎo)過(guò)程中，我們假設(shè)ci僅含有兩個(gè)缺失元素，若含有n個(gè)缺失元素，則Ai為一個(gè)n×n的P⊥2f×2f子矩陣，其元素取自P⊥2f×2f對(duì)應(yīng)于未知元素所在的行和列；bi為一個(gè)n×1的向量，其元素為ci中已知元素對(duì)應(yīng)的P⊥2f×2f行中的元素和ci中已知元素相乘之和。

同時(shí)，在上述推導(dǎo)中，我們僅針對(duì)M2f×s矩陣中的第i列ci，對(duì)于所有的其他列，求解都是一樣的，則總余差為：

e=∑si=1ei（15）

在上面的推導(dǎo)過(guò)程中，我們假設(shè)投影矩陣P⊥2f×2f已知，而事實(shí)上P⊥2f×2f由M2f×s矩陣的奇異值分解得到。但在我們所研究的問(wèn)題中，由于出現(xiàn)了遮擋現(xiàn)象，低秩矩陣M2f×s中產(chǎn)生了缺失元素，因此，無(wú)法直接進(jìn)行奇異值分解。因此算法開(kāi)始時(shí)，可以假設(shè)缺失元素的初值由其他已知元素的值的均值確定，然后利用列約束條件，求出缺失元素，此時(shí)求出的缺失元素比開(kāi)始時(shí)給定的初值要更加接近真實(shí)值，由此可以構(gòu)造一個(gè)迭代算法，恢復(fù)出低秩圖像矩陣的缺失元素。endprint

4 算法總結(jié)與計(jì)算復(fù)雜度

現(xiàn)給出基于列約束的低秩矩陣恢復(fù)算法：

（1）假設(shè)所有的遮擋點(diǎn)m-i，j=1k∑kj=1mi，j（即第i幅圖像中所有已知分量的均值），迭代次數(shù)l=1，令ε為任意小的一個(gè)正數(shù)；

（2）利用式（3），對(duì)M2f×s進(jìn)行奇異值分解；

（3）利用式（7）求投影矩陣P⊥2f×2f；

（4）對(duì)于圖像矩陣的每一列，利用式 （14）求出該列的缺失元素xli，并將求出的元素xli代替M2f×s中的缺失元素；

（5）若‖xli-xl-1i‖

SymbolcB@ ε則停止；否則l=l+1，轉(zhuǎn)第（2）步。

由于式（13） 是一個(gè)半正定二次型，從而有唯一的極小值點(diǎn)，而且極小值可以通過(guò)式 （14）直接求出。又由于它是一個(gè)凸函數(shù)，所以迭代誤差滿足

∑si=1e（1）i≥∑si=1e（2）i≥…≥∑si=1e（l）i （15）

因此，理論上保證算法能夠收斂到全局最小值。

本算法的每一步迭代需要進(jìn)行一次奇異值分解，以及計(jì)算P⊥2f×2f，Ai，bi，xi。由于M2f×s奇異值分解的時(shí)間復(fù)雜度是O（f2s），P⊥2f×2f的計(jì)算復(fù)雜度是O（f2）。另外由于Ai是P⊥2f×2f的二階子矩陣，bi是二維向量，其計(jì)算復(fù)雜度為O（f）。而xi=A-1ibi的計(jì)算復(fù)雜度是一個(gè)常數(shù)，與f和s均無(wú)關(guān)。所以本算法的計(jì)算復(fù)雜度為O（f2sl），其中l(wèi)為算法收斂所需要的迭代次數(shù)。

5 仿真實(shí)驗(yàn)

為了檢驗(yàn)本文所提方法的收斂速度和恢復(fù)效果，現(xiàn)在對(duì)含有大量遮擋點(diǎn)的圖像序列進(jìn)行仿真實(shí)驗(yàn)。實(shí)驗(yàn)的實(shí)驗(yàn)對(duì)象是仿真圖像。本文首先用計(jì)算機(jī)模擬的方法得到100個(gè)三維空間點(diǎn)作為特征點(diǎn)，即S為3×100矩陣。另外為了得到矩陣R，模擬相機(jī)的位置變化，得到一個(gè)包括50幅大小640×480的圖像序列，即R為100×100矩陣。隨機(jī)選擇50%的特征點(diǎn)進(jìn)行遮擋，并在圖像中依次加上均值都為0，方差為0.5、1、1.5、2、2.5、3個(gè)像素的高斯噪聲。利用這些模擬得到的圖像點(diǎn)用本文算法進(jìn)行圖像矩陣的恢復(fù)，算法的迭代收斂圖如圖1所示。

同時(shí)，為了研究本文算法的收斂性與圖像點(diǎn)的遮擋率之間的關(guān)系，在圖像中加入2.0個(gè)像素的高斯噪聲，并分別遮擋20%，30%，40%，50%及60%的圖像點(diǎn)，實(shí)驗(yàn)結(jié)果如圖2所示。

從圖1和圖2可以看出，本文所提出的方法迭代20次以內(nèi)就能夠達(dá)到收斂，具有良好的收斂性能，原因是由于本文算法將遮擋點(diǎn)的求解轉(zhuǎn)化為二次型的迭代求解，而二次型的極值點(diǎn)可以一步求解。同時(shí)，從圖2還可以看出，遮擋率越高，算法收斂速度越慢，原因是由于遮擋率越高，圖像序列的冗余性也少，約束越少，因此求解的未知數(shù)就越多，因此計(jì)算量就更大，算法的收斂速度就越慢。

為了比較本文算法和Martinec方法的收斂速度和恢復(fù)誤差，用上述方法產(chǎn)生50幅圖像，遮擋率取為20%，并在圖像中依次加入均值為0，方差從0變化到3的高斯噪聲，分別用本文算法和Martinec方法進(jìn)行圖像矩陣恢復(fù)。在每種噪聲下各重復(fù)100次實(shí)驗(yàn)，然后求取平均誤差，實(shí)驗(yàn)結(jié)果如圖3所示。

從圖3可以看出，相比Martinec方法本文算法具有更高的恢復(fù)精度，原因主要有三方面，一方面是由于本文算法不需要計(jì)算基礎(chǔ)矩陣，基礎(chǔ)矩陣的計(jì)算對(duì)圖像噪聲非常敏感；其二方面是由于本文算法將遮擋點(diǎn)的求解轉(zhuǎn)換為二次型的迭代求解，從而對(duì)應(yīng)優(yōu)化問(wèn)題的目標(biāo)函數(shù)是一個(gè)凸函數(shù)，它可以一步求解到全局最優(yōu)值，而且算法能夠保證收斂；其三方面是由于本文算法把所有的圖像平等地看待，并沒(méi)有倚重某些圖像。因此，本文算法具有較小的誤差。

6 結(jié)束語(yǔ)

本文在正投影相機(jī)模型下，提出一種基于列約束的低秩矩陣恢復(fù)方法。該方法利用圖像矩陣是一個(gè)低秩矩陣從而圖像序列具有冗余性，利用奇異值分解由圖像矩陣的列空間構(gòu)造出一個(gè)投影矩陣，將遮擋點(diǎn)的求解轉(zhuǎn)化為迭代求解二次型的極值問(wèn)題，利用它恢復(fù)出遮擋點(diǎn)，補(bǔ)上圖像矩陣的缺失元素。由于該方法的關(guān)鍵在于迭代求解二次型的極值問(wèn)題，而二次型是一個(gè)凸函數(shù)，它可以一步求解到全局最優(yōu)值，因此具有良好的收斂性能。仿真實(shí)驗(yàn)結(jié)果表明，本文方法相比Martinec方法具有收斂速度快，恢復(fù)精度高等優(yōu)點(diǎn)。由于在此算法中我們只考慮了基于列向量的約束，進(jìn)一步的研究方向?qū)⑹前鸦谛邢蛄亢土邢蛄康募s束同時(shí)考慮，以提高算法的效率。

參考文獻(xiàn)

[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，F(xiàn)RAHM 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