基于截斷奇異值低秩矩陣恢復的Canny邊緣檢測算法*

郭 偉，董宏亮，趙德冀

(遼寧工程技術大學軟件學院，遼寧 葫蘆島 125105)

1 引言

邊緣特征作為圖像的底層特征，是圖像分割、模式識別、計算機視覺等領域的研究基礎，邊緣檢測結果的好壞直接影響著后續(xù)高層次的圖像處理和分析。近年來，越來越多的新技術引入到了邊緣檢測算法中，如水平集[1]、距離圖[2]、遺傳算法、神經網絡、Gabor濾波器[3]等。Canny是基于最小二乘法推導出來的多級邊緣檢測算子，因具有較好的信噪比和較高的檢測精度而被廣泛應用。由于經典的Canny邊緣檢測算法中使用的高斯濾波器是線性濾波器，在對圖像進行平滑處理時，只考慮到像素間的空間距離關系，而沒有考慮像素值之間的相似性，是從圖像的整體結構來修改圖像，使得圖像變得模糊，邊緣信息被平滑處理，因而定位精度和信噪比下降，在很大程度上影響了后續(xù)的工作。因此，許多研究者提出了改進型的邊緣檢測算法[4 - 14]，并取得了較為理想的效果。Canny算法雖然應用廣泛，但也存在著不足，尺度較小的各向同性高斯核可以提取圖像的細節(jié)變化，但對噪聲非常敏感，而尺度較大的高斯核有較好的魯棒性，但定位精度較低。由此，李德軍等人[4]提出了基于雙邊濾波的圖像邊緣檢測算法，不僅考慮了圖像空間上的鄰近關系，同時也考慮到了圖像灰度上的相似性。此算法對于受低密度高斯噪聲干擾的圖像有很好的平滑效果，而對于受到椒鹽噪聲干擾的圖像檢測效果卻不理想。考慮到脈沖信號對邊緣檢測的影響，戚曉偉等人[5]提出了改進的雙邊濾波算法，首先對圖像進行中值濾波，在一定程度上減小了椒鹽噪聲的影響，由于中值濾波的效果受濾波窗口大小的影響，因此采用固定大小的濾波窗口很難達到理想效果。張潔等人[6]采用各向異性擴散方程來進行濾波操作，該算法雖然有效地去除了偽邊緣，但對椒鹽噪聲的處理結果仍然不理想。馬蒙蒙等人[7]采用空間尺度自適應的各向同性高斯濾波算法，自適應地獲取空間尺度。凌鳳彩等人[8]設計了新的梯度算子以及非極大值抑制的條件。李敏花等人[9]在Canny算法的基礎上設計了自適應的高斯濾波器尺寸和滯后閾值，提高了邊緣檢測算法參數設定的靈活性，但該算法在含有脈沖信號噪聲干擾的情況下仍存在著不足。而在實際生產應用中，由于傳感器故障、光照等原因，會造成圖像中出現噪聲點或者離群點，這時傳統邊緣檢測算法就顯得無能為力。Thombare等人[10]結合分割所具有的廣泛適用性和抵制噪聲的優(yōu)點，且基于塊統計來計算閾值，提出了基于閾值分割的分布式Canny邊緣檢測算法。Canny算法敏感于噪聲，在濾除噪聲時容易失去弱邊緣信息，且固定參數降低了算法適應性。針對這些問題，Rong等人[11]用重力場強度的概念代替圖像梯度，并獲得了引力場強度算子，提出一種改進的Canny算法。馮珂等人[12]首先利用非線性擴散濾波器有效去除噪聲，保持圖像的邊緣信息；其次，在鄰域梯度幅度的計算中考慮了像素對角線方向的梯度計算，進一步抑制了噪聲的影響；最后，使用平均類間方差可以自適應地計算不同圖像的雙閾值。文獻[13]對低秩矩陣恢復進行了理論分析，提出將矩陣分解為低秩矩陣和稀疏矩陣，由于其具有較強的魯棒性，被廣泛地應用于人工智能、圖像處理、計算機視覺等領域。據此，牛發(fā)發(fā)等人[14]提出了基于魯棒主成分分析法RPCA(Robust Principal Component Analysis method)的Canny邊緣檢測算法。將邊緣檢測算法作用于低秩矩陣上，從而實現對圖像的邊緣檢測。RPCA應用于圖像時，可以認為其將圖像矩陣分解為無噪聲圖像矩陣和噪聲圖像矩陣之和，其中無噪聲圖像矩陣具有低秩性，噪聲圖像矩陣具有稀疏性。由于RPCA分解得到的主成分可以去除少量的椒鹽噪聲及離散噪聲，所以在含有低密度椒鹽噪聲的情況下此算法得到了較好的檢測效果。由于低秩矩陣數據具有較強的全局描述能力和抗干擾能力，能夠更充分地利用相似塊之間的非局部信息，更好地保護圖像原有信息，因此，利用矩陣的低秩性可以較好地恢復出原數據矩陣，以達到較好的檢測效果。

本文在研究了經典Canny算法以及低秩矩陣恢復相關知識的基礎上，結合二者各自的優(yōu)點，提出一種基于低秩矩陣恢復的Canny邊緣檢測算法，并提出一種新的凸優(yōu)化模型及求解方法，將邊緣檢測算法直接作用于矩陣的低秩部分，既可以保留Canny邊緣檢測算法的優(yōu)越性又能增強算法的抗噪抗干擾能力。

2 數學工具

本節(jié)簡單地介紹一下本文與矩陣相關的基本概念。首先介紹一種非常重要的分解模式，即奇異值分解SVD(Singular Value Decomposition)。

定義4[16,17]矩陣A∈Rm×n的奇異值分解為：A=USVT,S=diag({σi}1≤i≤min(m,n))，定義奇異值閾值算子：Θτ(A)=UΘτ(D)VT,Θτ(Σ)=diag(max{σi-τ,0})，其中σi是矩陣A中的第i個奇異值。

3 低秩矩陣恢復算法分析

低秩矩陣恢復又稱魯棒主成分分析或稀疏與低秩分解。從傳統的主成分分析PCA(Principal Component Analysis)方法的角度來看待矩陣恢復問題，即傳統的PCA方法是一種無監(jiān)督的、使用SVD對復雜數據進行降維的方法，同時也可以認為其是一種去相關的方法，這種方法可以有效地找出數據中最“主要”的元素和結構，去除噪聲及不相關部分，揭示隱藏在復雜數據背后的簡單結構。PCA方法的本質是高維數據在其低維線性子空間上的投影，其目標就是使用另一組基去重新描述得到的目標空間，希望能盡量揭示原有的數據間的關系?？梢员硎龀扇缦碌臄祵W模型：

min ‖S‖F,

s.t.rank(A)≤r,D=A+S

(1)

其中,D為輸入的數據矩陣，S為誤差項，‖S‖F表示矩陣的Frobenius范數，r為目標子空間的維數。通過上述的優(yōu)化約束問題可以找到數據矩陣D在一個最近的r維的線性子空間上的投影，也就是說要估計這個低維子空間的數學模型就是要找到一個低秩矩陣A，使得其與觀測矩陣D之間的差異最小化,稱為損失最小化。當S服從Gaussian分布并且獨立時，PCA方法可以通過一次SVD準確地找到最優(yōu)逼近解A，此時PCA方法具有較好的表現。而當A受到稀疏噪聲干擾時，PCA方法并不能很好地完成降維或者去相關操作，也就是說PCA方法對于離群點缺少魯棒性，估計出的解往往不是最優(yōu)解，并且PCA方法還需要預設目標子空間的維數r。于是Wright等人[18]提出了RPCA方法來解決矩陣中受到稀疏或者“污點”噪聲干擾的情況。RPCA方法考慮的是這樣一個問題：一般收集到的數據矩陣中包含結構信息的同時也包含噪聲信息。那么，我們可以將這個矩陣分解為兩個矩陣之和，一個是低秩的(由于內部有一定的結構信息，造成各行或列間是線性相關的)，另一個是稀疏的(存在的噪聲是稀疏的)。與經典PCA問題一樣，RPCA問題本質上也是尋找數據在低維空間上的最佳投影問題。對于低秩數據觀測矩陣D，假如D受到隨機(稀疏)噪聲的影響，則D的低秩性就會破壞，使D變成滿秩的。所以，就需要將D分解成包含其真實結構的低秩矩陣和稀疏噪聲矩陣之和，以消除稀疏噪聲的干擾。那么，RPCA可以描述成原始數據矩陣D由兩個部分組成，即低秩矩陣A和稀疏矩陣S，原始數據矩陣可以表示為D=A+S。并且只要誤差矩陣S是稀疏的，不論其大小，找到了低秩矩陣A，實際上就找到了數據的本質低維空間。那么,魯棒主成分分析可以描述成如下的優(yōu)化問題：

minrank(A)+λ‖S‖0,

s.t.D=A+S

(2)

其中,目標函數表示為矩陣A的秩函數及誤差矩陣S的零范數，λ>0表示平衡因子，我們希望在合適的平衡因子λ的作用下能夠精確地恢復出低秩矩陣A，然而，公式(2)的優(yōu)化問題是一個NP-Hard非凸問題，在理論和實踐過程中均只存在指數級復雜度的算法，并不能在多項式時間內直接求解。因此，需要尋找可以直接求解的優(yōu)化問題來近似原問題。Candès等人[19]從理論上證明了L1范數最小化近似于L0范數最小化解，這樣就可以將L0范數最小化問題松弛到L1范數最小化問題，而矩陣的rank()是非凸不連續(xù)函數，是奇異值的L0范數，而核范數是奇異值的L1范數，因此根據上述理論，可以采用凸松弛法，采用核范數來近似矩陣的rank()函數，L1范數近似L0范數。將式(2)轉化為如下凸優(yōu)化問題：

min‖A‖*+λ‖S‖1,

s.t.D=A+S

(3)

4 截斷奇異值模型的提出

min ‖A‖t+λ‖S‖1,

s.t.D=A+S

(4)

然而，截斷奇異值方法同樣是非凸的，因此式(4)同樣無法直接求解。根據矩陣相關知識，由定理1和引理1可以容易地推出如下結論：

(5)

因此根據式(5)，式(4)可以轉化為求解如下優(yōu)化問題：

s.t.D=Z+S

(6)

因此,本文根據截斷奇異值模型提出基于截斷奇異值的魯棒主成分分析方法RPCA-TSV(Robust Principal Component Analysis method with Truncated Singular Value)。

5 新模型的提出及算法描述

根據文獻[20]中相關介紹，用于求解魯棒性主成分分析的方法是基于這樣一個假設：原始數據矩陣由低秩矩陣部分和稀疏誤差矩陣兩部分組成。在現實中，圖像在獲取和傳輸的過程中往往會遭受到各種噪聲的污染，尤其是高斯噪聲，從而降低了圖像的質量，給后續(xù)的圖像處理帶來諸多的不利影響。也就是說現實中的數據矩陣常常同時受到高斯噪聲和稀疏噪聲的干擾，使得算法受到噪聲干擾而變得不準確。因此，本文提出一種新的凸優(yōu)化模型，把原始數據分為低秩部分、稀疏噪聲部分和高斯噪聲三部分，并添加高斯正則化項，構造雙噪聲凸優(yōu)化模型:

min ‖Z‖t+λ‖S‖1+β‖G‖F,

s.t.D=Z+S+G

(7)

其中,‖Z‖t表示截斷奇異值，‖G‖F表示矩陣的Frobenius范數，作為高斯噪聲的約束項，β為高斯噪聲與低秩矩陣之間的權衡參數。

5.1 優(yōu)化模型求解方法

求解凸優(yōu)化模型的最優(yōu)解是一個關鍵問題，實際上，各個模型均可以使用迭代閾值法[21]來求解，但其中存在著如何選取迭代步長和迭代次數的問題。針對迭代閾值法的缺陷，Lin等人[22]提出了加速近端梯度算法和對偶算法。加速近端梯度算法使用部分奇異值分解減輕了求解過程對奇異值分解的依賴，提高了求解魯棒主成分分析問題時的收斂速度和精度，但是當奇異值個數小于某一閾值時，其收斂速度甚至比使用完整的奇異值分解還要緩慢。而對偶算法在求解的過程中只依賴于最大的奇異值相關主奇異空間，然而目前并不存在一個十分有效的方法來計算該主奇異空間。隨后Lin等人[23]又提出了增廣拉格朗日乘子法來求解優(yōu)化問題，使得其迭代解能夠收斂到該優(yōu)化問題的精確解，并且耗時相對較少。

采用交替更新的方法，當更新矩陣Z時，固定矩陣S、Y和G，求解下面的優(yōu)化模型：

(8)

其中,Θμ-1(D)為伸縮奇異值算子。當更新S時，固定其他三項Z、Y和G，求解下面的優(yōu)化模型：

(9)

其中,δε(A)為伸縮絕對值算子。

當更新G時，固定其他三項Z、Y和S，求解下面的優(yōu)化模型：

(10)

當更新Y時，固定Z、S和G，此時Y的更新公式為：

Yk+1=Yk+μk(D-Zk+1-Sk+1-Gk+1)

(11)

更新μ:

(12)

其中,ρ為大于1的常數，ε為大于0且比較小的正整數。文獻[24]中給出ρ=1.5,ε=10-8。初始Lagrange懲罰因子μ=0.25/mean，其中mean是D中所有元素的平均值。

5.2 截斷奇異值魯棒主成分分析求解流程

交替迭代更新算法流程描述如下：

步驟1初始化數據矩陣以及迭代誤差:

設定初始數據矩陣M=D，懲罰因子μ，迭代收斂誤差eps=ε‖D‖F，終止誤差ε。

步驟2計算中間變量Ak和Bk:

將初始化矩陣M進行奇異值分解即：[Uk,Σ,Vk]=svd(M)，其中Uk=(u1,u2,u3,…,um)∈Rm×n,Vk=(v1,v2,v3,…,vn)∈Rm×n，從而得到Ak=(u1,u2,u3,…,ur)T,Bk=(v1,v2,v3,…,vr)T。

步驟3使用交替更新方法求解式(7)的優(yōu)化問題：

(1)用式(8)更新變量Z；

(2)用式(9)更新變量S；

(3)用式(10)更新變量G；

(4)用式(11)更新拉格朗日乘子Y；

(5)用式(12)更新懲罰因子μ。

直到內部算法收斂，即‖Zk+1-Zk‖F≤eps,‖Sk+1-Sk‖F≤eps,‖Gk+1-Gk‖F≤eps。

步驟4轉到步驟2重復循環(huán)更新變量，直到外部算法收斂，即‖D-Z-S-G‖<ε‖D‖F。

5.3 本文邊緣檢測算法流程

本文算法首先求取輸入圖像的主成分部分，即低秩部分，然后將經典的Canny邊緣檢測算法作用于主成分上，最后輸出邊緣檢測結果。算法流程如圖1所示：

Figure 1 Algorithm flow chart圖1 算法流程

6 實驗結果與評價分析

6.1 算法評價指標

為了從客觀的角度去評價一個算法的效果，本文采用以下評價標準。

6.1.1 峰值信噪比PSNR

采用峰值信噪比PSNR(Peak Signal to Noise Ratio)來評價圖像的去噪效果，其定義如下[25]：

(13)

其中,L表示灰度等級的最大取值范圍，本文中圖像的灰度級為0～255;M，N分別表示圖像的長度和寬度；i(x,y)表示原始圖像;i′(x,y)表示去噪后的圖像。PSNR越大表示去噪的效果越好。

Figure 2 Detection results of Lenna under salt noise of different densities 圖2 Lenna圖在不同密度椒鹽噪聲環(huán)境下的檢測結果

6.1.2 品質因數FOM

品質因數FOM(Figure Of Merit)用于衡量邊緣保存能力，其定義如下[26]：

(14)

其中,NID表示理想的邊緣像素數目，NDE表示圖像檢測后的像素數目，di表示檢測出的像素與離其最近的真實像素的距離，β是一個常量，一般取1/9。FOM是歸一化的結果，其取值為0～1，與邊緣的定位精度成正比，即FOM越大表示檢測的邊緣越準確。

6.1.3 邊緣保持指數EPI

邊緣保持指數EPI(Edge Preservation Index)用于評價邊緣保持能力[27],其定義為:

(15)

其中,x表示無噪聲原始圖像，y表示去除噪聲后的圖像，Δ表示Laplacian算子。EPI的大小與邊緣保持能力成正相關，其值越大表示邊緣保持得越好。

6.2 實驗結果分析

為驗證本文算法的優(yōu)越性，采用256×256像素的Lenna和Pepper標準圖作為輸入，由于各種算法在對無噪聲圖像進行檢測時效果類似，而在處理含有椒鹽噪聲圖像時，因較為敏感而沒有很好的效果。因此，本文采用文獻[5,14]以及本文算法對以下兩種情況進行驗證。

(1)僅受椒鹽噪聲干擾。

圖2和圖3分別顯示了Lenna圖和Pepper圖在受不同密度(0.05和0.1)椒鹽噪聲干擾時各個方法的檢測結果。圖2中文獻[5]的方法雖然可以很好地去除椒鹽噪聲的干擾，但是會丟失細節(jié)信息,其原因是文獻[5]首先對圖像采用固定大小窗口的中值濾波，然后再使用雙邊濾波進行處理，而造成了對圖像的過度平滑，使得圖像的檢測效果下降。使用RPCA方法雖然可以得到較好的圖像的邊緣檢測結果，同樣也會把噪聲當成邊緣而進行誤檢，產生偽邊緣。比如說圖2中Lenna帽子上的偽邊緣。而本文算法由于采用更加準確的表示方式來刻畫模型中的秩函數而得到更加理想的邊緣檢測效果。表1顯示了圖2在不同密度椒鹽噪聲干擾情況下對應的評價標準結果，無論是在客觀上還是在主觀上，都可以看出本文算法的準確性較高、抗干擾能力較強。

(2)椒鹽噪聲和高斯噪聲混合干擾。

為了驗證本文算法的魯棒性，針對椒鹽噪聲和高斯噪聲混合噪聲進行檢測處理。圖4和圖5顯示了受到椒鹽噪聲和高斯噪聲混合噪聲干擾時各個算法的邊緣檢測結果，檢測圖像中分別添加椒鹽噪聲(0.05)和高斯噪聲(0.002)的混合噪聲、椒鹽噪聲(0.1)和高斯噪聲(0.002)的混合噪聲。由于在魯棒主成分分析過程中采用了更加精確的截斷核范數形式來代替秩函數，因此本文算法的檢測結果與RPCA方法相比更加準確；同時，本文方法增加了高斯噪聲的約束項，使得在混合噪聲干擾的情況下具有較好的檢測效果。表2顯示了圖4所示的混合噪聲情況下各算法對應的評價標準，從檢測結果和評價標準可以看出，本文算法具有較好的檢測結果，體現了本文算法的魯棒性。

Table 1 PSNR,FOM,EPI of each algorithm under noise of different densities表1 不同密度噪聲下各算法對應的PSNR，FOM，EPI

Figure 3 Detection results of Pepper under salt noise of different densities 圖3 Pepper圖在不同密度椒鹽噪聲環(huán)境下的檢測結果

Figure 4 Detection results of Lenna in mixed noise environments with different densities圖4 Lenna圖在不同密度混合噪聲環(huán)境下的檢測結果

Figure 5 Detection results of Pepper in mixed noise environments with different densities圖5 Pepper圖在不同密度混合噪聲環(huán)境下的檢測結果

Table 2 PSNR,FOM,EPI of each algorithm under mixed noise environments with different densities表2 不同混合噪聲密度下各算法對應的PSNR，FOM，EPI

7 結束語

針對圖像邊緣檢測這一圖像的底層特征提取問題，對于噪聲圖像，傳統的邊緣檢測方法抗噪性能較差，往往不能很好地獲取其邊緣特征。本文通過改進的截斷核范數魯棒主成分分析方法將圖像分解成低秩矩陣、稀疏誤差矩陣和模擬的高斯噪聲3個部分，去除原始數據矩陣中的稀疏噪聲和冗余部分，起到增強原始數據矩陣的作用，有利于后續(xù)邊緣信息的提取及圖像處理工作。同時，結合經典Canny邊緣檢測算子的優(yōu)點，將其作用于分解后的低秩主成分上。因此，本文算法可以很好地消除椒鹽噪聲以及混合噪聲的影響，既能增加算法的魯棒性，又能提高算法的去噪性能，提高邊緣檢測的準確性。