基于曲率濾波和N-P準則的路面裂縫識別方法*

王墨川,何 莉,胡成雪,陶 健,張德津

(1.湖北工業大學電氣與電子工程學院，湖北 武漢 430068；2.深圳大學機電與控制工程學院，廣東 深圳 518060；3.深圳大學廣東省城市空間信息工程重點實驗室，廣東 深圳 518060)

1 引言

裂縫是瀝青路面的主要破損情況，準確高效地檢測裂縫對路面養護有著積極的意義。瀝青路面圖像紋理復雜，裂縫形狀多樣、尺度多變，特別是細小裂縫與路面背景對比度低，使瀝青路面裂縫自動檢測有一定困難。

Cheng等[1]基于背景差法消除光照不均勻，對于陰影等低頻噪聲的處理效果較好，但容易平滑掉裂縫信息。基于多尺度幾何的Ridgelet方法[2]實現了裂縫目標的有效增強，但模型參數選擇復雜，算法適應性一般。基于閾值化的圖像分割得到了廣泛研究，Jitprasithsiri等[3]使用不同路面裂縫圖像對比了不同閾值分割方法，發現回歸方法的效果最好。Oh等[4]針對回歸方法的缺點，提出了一種迭代的裁剪法，使用均值和方差去除噪聲像素，保留裂縫像素。許多研究人員將機器學習的方法用于路面裂縫識別[5,6]，將路面圖像劃分成許多小的子圖像，然后針對子圖像進行特征提取構建特征向量，進而通過樣本訓練實現對新樣本分類。Oliveira等[5]計算一系列統計量，如均值、方差，作為類別特征值。Nguyen等[6]將各向異性作為類別特征。對于路面裂縫的連接，鄒勤等[7,8]采用了分別基于最小生成樹和最短路徑的路面裂縫連接算法，此類算法根據距離因素，設置經驗閾值對目標點內、外連接，對相互鄰近、交叉的復雜裂縫(如龜裂)的檢測效果不理想。

本文提出一種基于曲率濾波和N-P準則的路面裂縫識別方法。通過改進曲率濾波算法平滑干擾并保留裂縫信息，采用二次分割及噪聲去除的方法對裂縫進行定位，基于N-P準則實現裂縫片段端點的連接，最終實現裂縫檢測。

2 圖像特征分析

裂縫圖像具有如下特征：(1)灰度特征：裂縫相對于紋理是一些灰度值較低的像素集合，裂縫主要分布在灰度直方圖左側的一小部分區域內，如圖1a所示。(2)幾何特征：裂縫的幾何形狀接近線狀和帶狀，有較大的長度和較大的長寬比，而紋理的長度較短，或呈點狀，如圖1b所示。(3)位置特征：復雜裂縫由多條簡單裂縫發展而來，即空間上距離較近的簡單裂縫，可能相互滲透和影響，發展為復雜裂縫，如圖1c所示。(4)方向特征：路面裂縫的局部方向與裂縫的整體方向基本一致，如圖1c所示。

一幅含有裂縫的圖像由3部分組成，設S(I)表示含有裂縫的路面圖像，則可表示為：

S(I)=Sb(I)+Sn(I)+Sc(I)

(1)

其中，Sb(I)是紋理噪聲信號，Sn(I)是隨機噪聲信號，Sc(I)是裂縫信號，I是像素坐標{x,y}。由于Sb(I)和Sn(I)的存在，使得Sc(I)淹沒其中，很難直接提取出Sc(I)。考慮干擾濾除和裂縫保持要求，本文通過結合最小矩形切平面和最小三角形切平面、修正正則化能量的曲率濾波，實現消除隨機噪聲、平滑紋理，并盡可能保留裂縫的目的。在此基礎上，本文參考裂縫灰度特征，以子塊為單元，基于直方圖擬合的方法粗分割Sc(I)與Sb(I)，去除不滿足要求的區域，獲得感興趣區域ROI(Region of Interest)，通過ROI子塊對非ROI子塊進行分類后，對所有子塊再次分割；根據幾何特性，去除與裂縫具有相似灰度特征的Sb(I)；最后融合位置和方向信息，連接斷裂的裂縫片段，實現裂縫完整檢測。本文方法過程如圖2所示。

3 路面裂縫定位

3.1 基于曲率濾波的圖像平滑

全變分曲率濾波TVF(Total Variation Filter)和高斯曲率濾波GCF(Gaussian Curvature Filter)[9]是常用的濾波方法。對比發現,前者能更有效地去除噪聲，后者能更好地保持圖像的邊緣和紋理細節。本文組合全變分曲率濾波、高斯曲率濾波的最小矩形切平面、最小三角切平面，通過鄰域均值修正正則能量項，在去除干擾的同時保留裂縫邊緣信息。曲率濾波[9]的思想是定義一能量泛函，這個能量泛函的值越小，變量就越接近于期望的結果。能量函數的能量公式如式(2)所示:

E(S)=E0(S,I)+E1(S)

(2)

其中，E(S)表示總能量，E0(S,I)表示擬合項能量，E1(S)表示正則化能量。能量泛函的演化過程如圖3所示，因為擬合項能量E0(S,I)在增加，正則化能量E1(S)在減少，且目標是使總能量E(S)減少，所以從曲率濾波的角度得出在優化過程中正則化能量是主導部分，先將曲率正則化能量降至極小值，再觀察總能量決定是收斂還是繼續迭代。

本文改進曲率濾波算法的具體步驟如下:對輸入圖像U(i,j)構造最小矩形切平面和最小三角形切平面，如圖4所示。

d1=(U(i-1,j-1)+U(i-1,j)+U(i-1,j+1)+U(i,j+1)+U(i,j-1))/5-U(i,j),

d2=(U(i-1,j+1)+U(i,j+1)+U(i+1,j+1)+U(i-1,j)+U(i+1,j))/5-U(i,j),

d3=(U(i+1,j-1)+U(i+1,j)+U(i+1,j+1)+U(i,j-1)+U(i,j+1))/5-U(i,j),

d4=(U(i-1,j-1)+U(i,j-1)+U(i+1,j-1)+U(i-1,j)+U(i+1,j))/5-U(i,j),

d5=(U(i-1,j)+U(i,j-1)+U(i-1,j-1))/3-U(i,j),

d6=(U(i-1,j)+U(i,j+1)+U(i-1,j+1))/3-U(i,j),

d7=(U(i,j-1)+U(i+1,j)+U(i+1,j-1))/3-U(i,j),

d8=(U(i,j+1)+U(i+1,j)+U(i+1,j+1))/3-U(i,j),

dm=min{|di|,i=1,…,8},

通過多次實驗發現，曲率濾波的正則能量項僅考慮了曲率能量的約束，對于路面裂縫圖像，圖像的裂縫信息被嚴重污染，幾何結構被嚴重破壞，若使用結合的投影算子對能量函數進行迭代，如圖5a所示，總能量和數據擬合項能量在第2次迭代就開始增加，影響結合2種投影算子后的模型的平滑能力，使得結合2種投影算子后的算法無法收斂到較好的效果。

根據裂縫的灰度特征，裂縫鄰域均值較小，而非裂縫區域鄰域均值較大。受雙邊濾波器[10]的啟發，本文在考慮鄰域微分幾何空間的同時，結合鄰域內相鄰像素的灰度值差距，修正正則化能量。且正則化能量在第2次開始時是增加的，通過曲率正則化能量與鄰域像素作差以適應裂縫圖像的特征，提出一個新的正則能量項，如式(3)所示:

(3)

其中，K(U)為曲率正則化能量，Ψ=(y,U)表示幾何曲面，y為空間坐標，|gm|為8鄰域灰度值均值。

修正后的能量函數曲線如圖5b所示，正則化能量和總能量在多次迭代后可收斂至較平穩的結果，證明了本文方法的有效性。改進的曲率濾波結果如圖6d所示。通過改進曲率濾波濾除了干擾，并且可以保留裂縫信息，可通過圖7b進一步得到驗證，表明本文方法可減少圖像紋理像素灰度值變化，且不會顯著影響裂縫像素的灰度值。

3.2 二次分割裂縫提取

裂縫是空間上聚集灰度相似的像素點的集合，經改進的曲率濾波平滑圖像后，可通過閾值分割提取路面裂縫。本文采用子塊分析方法，利用二次分割的策略提取疑似裂縫。

(1)基于直方圖擬合的粗分割。

粗分割大致定位裂縫的灰度區間及可能的分布區域。由裂縫的灰度特征，裂縫像素主要分布在灰度直方圖左側的一小部分區域內，可以利用直方圖擬合的方法粗分割。

首先，將經改進曲率濾波后的圖像劃分為互不重疊的256×256像素的子塊圖像。其次，以子塊為單元，利用迭代式閾值分割算法[11]計算子塊圖像的灰度閾值iimg，在直方圖上選取灰度值小于iimg的區間。然后,利用高斯函數G(μ,s2)擬合位于最大值左側的直方圖區間(其中,μ和s2分別表示最大值左側直方圖的均值和方差),計算閾值Th，使P(Th)分別為65.39%,95.4%,99.7%，其對應Th=μ-s,Th=μ-2s,Th=μ-3s。最后，將灰度值小于Th的像素標記為“1”(裂縫像素)，其余標記為“0”(非裂縫像素)，完成圖像所有子塊分割。通過實驗采用不同閾值對3 050幅裂縫圖像進行檢測，其漏檢率和誤檢率如表1所示。

Table 1 Influence of different values of threshold on detection results

考慮漏檢率應盡可能小，而誤檢率應盡可能大，本文選取Th=μ-2s。圖8c是圖8b一子塊灰度值小于iimg的區間，黑點為高斯函數均值μ(125)，均值左側曲線表示擬合的高斯函數，計算出的閾值為Th=91。閾值分割的結果如圖8d所示。去除離心率不滿足要求的連通域，得到ROI，如圖8e所示，后續根據ROI對子塊進行分類并計算子塊圖像的二次分割閾值，因此提取的ROI必須準確可靠。

(2)子塊圖像自適應分類及細分割。

細分割進一步提取疑似裂縫。本文對含有ROI的子塊采用直方圖估計[12]的方法計算閾值，按不同閾值將子塊分為6類，如圖8f所示。比較非ROI子塊與已知類別子塊灰度均值、方差的相似性，將其歸為最為接近的那一類，實現所有子塊圖像的分類，如圖8g所示。將子塊分成若干類別，同時獲得每一類別的分割閾值，按照此分割閾值對子塊二次分割得到最終提取的裂縫，提取結果如圖8h所示。經二次分割，裂縫區域顯著增多。

3.3 裂縫噪聲自適應去除

經過裂縫提取后，裂縫主體已被識別出來，但存在大量的塊狀或點狀噪聲。噪聲產生的主要原因是部分小尺度紋理與裂縫的灰度區間相似，且紋理像素的坐標分布不連續，而文中的粗分割方法是基于直方圖的，對目標的空間信息不敏感，會將這些噪聲錯誤地識別為裂縫，這種情況下可根據裂縫表現出的幾何特征將二者區分開。Tong等[13-15]分別利用連通域的長寬比、離心率和面積率，通過設定經驗閾值來區分裂縫與噪聲的方法，都取得了較好的效果，但閾值的確定存在一定的主觀性。本文采用K-means算法對目標的塊狀或點狀噪聲進行自適應分類，而不依賴于特征因子的閾值。特征因子連通域的長寬比、離心率和面積率的定義如下所示：

定義1長寬比L-Wratio是指最小外接矩形的長L與寬W的比值，如式(4)所示。該比值反映裂縫的線狀特征,其值越大，說明連通域是裂縫的可能性越大。

(4)

定義2離心率e是指與連通域具有相同標準二階中心矩的橢圓的2個焦點間的距離c和長軸長度a的比值，如式(5)所示。該值反映了裂縫的線狀特征，其值越大，表明該連通域越接近于直線，是裂縫的可能性越大。

(5)

定義3面積率Sratio是指連通域的面積bwarea與最小外接矩面積的比值，如式(6)所示。Sratio越小，表明該連通域的線性特征越明顯，是裂縫的可能性越大。

(6)

聚類是一種無監督學習的方法，其原理是把數據集合中相似的對象分成不同的組別或者更多的子集，使在同一個子集中的成員都具有相似的屬性。裂縫連通域包含裂縫和噪聲2類，設定K-means算法中的類別數K=2，特征向量由3個特征因子構成，即X=(L-Wratio,e,Sratio)。圖8h的聚類結果如圖9所示，其中第1類表示噪聲，第2類表示裂縫。去除噪聲的結果如圖10所示，至此實現了裂縫的定位。

4 裂縫片段連接

受裂縫壁脫落、積灰等因素影響，裂縫僅在宏觀上呈現連續性，去除裂縫噪聲后得到的是裂縫的零碎片段，是不完整、不連續的裂縫，需要將裂縫片段的端點連接起來，以獲得較完整、連續的裂縫。并且并非所有裂縫片段的端點都需要連接，因此在連接過程中需要對所有片段進行判斷。

傳統的判斷方法通常采用貝葉斯準則進行判斷，該準則以先驗概率和損失函數為已知條件。對于如雷達信號的檢測，雷達信號的先驗概率是變化的且先驗概率的獲取十分困難，也無法對各種判斷結果給定損失函數。為此，本文采用奈曼皮爾遜(Neyman-Pearson)準則[16]，簡稱N-P準則。N-P判斷規則如式(7)所示：

(7)

其中，λ為檢測門限。可以連接的全體連通域對lla,b組成的類別為c1類，不能連接的全體連通域對lla,b組成的類別為c2類。a,b∈{1，…，n}為連通域的索引，n為二進制圖像中連通域的數量。X表示某連通域對，p(X∣c1)表示連通域對X在可以連接的條件下的連接概率密度函數。

在裂縫片段連接過程中，實際裂縫片段不能連接但錯誤判斷為可以連接的概率為虛警概率PF。N-P準則在給定PF的條件下，使得漏檢概率PM最小，即使得實際裂縫片段連接，判斷結果為不連接的錯誤判斷概率最小。本文融合位置和方向特性，得到2個特征因子，基于N-P準則對所有裂縫片段進行判斷，在判斷為連接的基礎上進行端點連接。具體步驟如下所示：

步驟1對所有連通域進行評估，考慮連通域的空間位置和方向特征得到2個特征因子x1a,b和x2a,b。(1)x1a,b：表示連通域a和b之間的最短距離。x1a,b指區域的位置特征，依據相鄰相似的原則，連通域在空間位置上距離越近,屬于同一裂縫的可能性就越大。(2)x2a,b表示連通域a和b具有相同標準二階中心矩的橢圓的長軸方向之間的最小角度。本文將與連通域具有相同標準二階中心矩的橢圓的長軸方向作為連通域的方向，x2a,b越小，裂縫的相似度越高。此時將連通域對X表示為向量形式，如式(8)所示：

(8)

步驟2假設每個隨機變量x1a,b,x2a,b都分別服從高斯分布且相互獨立，即N1(u1,σ1)，N2(u2,σ2),其中，u1和u2分別表示隨機變量x1a,b和x2a,b的均值，σ1和σ2分別表示隨機變量x1a,b和x2a,b的方差，則隨機變量x1a,b與x2a,b的二維均值向量和協方差矩陣分別如式(9)所示：

(9)

其中,ω∈{1,2}是類別。

步驟3根據全概率公式、正態分布的概率密度函數公式，利用式(8)和式(9)對式(7)進行推導，得到式(10)。當式(10)成立時，可將2個連通域中最近的端點連接起來。

(10)

其中，m∈{1,2}是隨機變量的維數。

經過上述步驟進行連接判斷及端點連接，最終實現裂縫識別。圖11a和圖11c為經裂縫噪聲去除后的效果圖，圖11b和圖11d分別展示的是圖11a和圖11c經裂縫片段連接的效果，如圖11a矩形框中的a1～a6、圖11c矩形框中的c1～c4經連接后，分別對應圖11b矩形框中的b1～b6、圖11d矩形框中的d1～d4。可以看出，本文裂縫片段連接方法能有效連接斷裂的裂縫。

5 實驗與結果分析

5.1 有效性驗證

本節利用公開數據集[17]，對現有較好的裂縫檢測方法Kittler法[18]、大津法[19]、松弛因子法[20]以及本文方法進行對比實驗，從測試效果和裂縫識別的完整性2個方面對實驗結果進行對比分析。圖12a第1列從上至下分別代表連續性差的橫縫、低對比度的縱縫、高對比度的縱縫及光照不均勻的網狀裂縫圖像。從圖12可知，本文方法與現有其它方法相比，對于漏檢、局部不連續等狀況有明顯的改善，且檢測出的裂縫清晰完整，與人工描點的裂縫較吻合，其檢測結果明顯優于其它方法。

本文使用機器學習中常用的精確率P(Precision)、召回率R(Recall)以及F值[14]進行完整性評估，各指標的定義如式(11)所示:

(11)

其中，TP表示真陽性,FP表示假陽性,FN表示假陰性。通過式(11)可得誤檢率(1-P)和漏檢率(1-R)。將不同方法的識別結果與人工描點進行比對，得到裂縫完整性評估結果，如圖13所示，本文方法的P、R、F均高于其他3種方法的。

5.2 工程數據分析

本文還對道路檢測裝備采集的工程數據進行了實驗。圖14a第1列從上至下分別代表橫裂、縱裂、塊裂及龜裂縫圖像，各行表示對于不同裂縫圖像，本文方法的各步驟圖及人工描點圖。

具體分析如下：

(1)對第3.1節中的曲率濾波方法與改進的曲率濾波方法的平滑度進行測試，結果如圖15所示。改進的高斯曲率濾波方法的P、R指標均高于曲率濾波的，分別提高了5.6%,9.4%。

(2)對第3.2節中所示粗分割和二次分割結果進行評估，評估結果如圖16所示。經過第二次分割后，評估指標F值提高了8.7%。

(3)對第3.3節中的裂縫噪聲去除，分別利用長寬比[13]、離心率[14]和面積率[15]，通過設置經驗閾值的方法，與本文方法進行對比實驗，評估結果如圖17所示。采用本文方法對裂縫噪聲消除最好，誤檢率分別降低了44.1%，34.5%和61.2%，漏檢率分別降低了78.1%，76.5%和79.3%。

(4)裂縫片段連接的判斷過程中，需要設置λ，從而獲得虛警概率PF與檢測概率PD=1-PM。設λ的值在[3.5,5.5]，以0.1為步長計算PF，部分結果如表2所示。繪制受試工作特性曲線如圖18所示，計算坐標點(0,1)到曲線最近距離對應的坐標A(0.0320，0.9394)，即PF=0.0320，根據表2查出，此時λ為5.0，檢測概率PD為0.939 4。

Table 2 Correspondence between λ and PF

(5)利用工程數據中的裂縫圖像分別對第3.3節(連接前)與第4節(連接后)測試裂縫識別的完整性，連接前P、R和F值分別為81.1%,79.8%和80.4%，連接后P、R和F值分別為91.2%,90.2%和90.7%，P、R、F值分別提高了12.5%，13.1%和12.8%。相比本文提出的濾波、分割方法，本文提出的連接方法對裂縫的完整性檢測起到了至關重要的作用。

6 結束語

本文提出了一種基于曲率濾波與N-P準則的路面裂縫識別方法，通過組合最小矩形切平面和最小三角形切平面及修正正則化能量的改進曲率濾波算法，達到了消除隨機噪聲、平滑路面紋理噪聲和裂縫保留的目的。以子塊為單元，利用裂縫灰度特征采取二次分割的策略提取疑似裂縫目標集合，運用裂縫幾何特性自適應去除塊狀或點狀噪聲，進而獲得裂縫片段和實現裂縫的定位。在此基礎上，融合裂縫片段的位置和方向信息，基于N-P準則連接裂縫片段端點，實現裂縫的完整性檢測。與其它裂縫檢測方法比較，本文方法可以對不同特征、位置的裂縫目標進行有效檢測，具有更好的檢測效果和更高的檢測精度。

本文研究建立在裂縫與背景灰度形成差異的基礎上，對于裂縫與紋理灰度相似的裂縫的檢測具有局限性。未來將結合裂縫灰度及脊邊緣特征的不同優點，在裂縫多特征融合方面展開深入研究，進一步提高識別的準確性。