基于模糊局部信息C均值的太赫茲圖像目標檢測

江 天，沈會良，楊冬曉*，劉建軍，鄒 哲

(1.浙江大學信息與電子工程學系，杭州310027;2.中國計量學院太赫茲技術與應用研究所，杭州310018)

引 言

太赫茲(terahertz，THz)波段的頻率范圍為0.1THz～10THz。太赫茲輻射作為一種波源，和其它波段的輻射一樣，可以作為物體成像的信號源。太赫茲波對于很多非極性電介質材料(如塑料、陶瓷、紙箱、布料等)有很好的穿透性。太赫茲成像是近年發展快速的一種成像技術，已成為X射線成像、毫米波成像等成像技術的有力補充，在無損檢測、安全檢查、病變組織檢測等領域體現出了特有的實用價值。

根據太赫茲波的形式，大體上可以分為太赫茲脈沖成像和太赫茲連續波(continuous wave，CW)成像。太赫茲脈沖成像不僅能夠獲得物體的形狀，還可以得到物體的光譜信息，但由于需要進行時間掃描與空間掃描，成像速度緩慢，因此常用于物體的波譜檢測。太赫茲連續波成像雖然只能獲得物體反射或透射的太赫茲波強度信息，但只需要進行空間掃描，成像速度比太赫茲脈沖成像快得多。由于成像系統結構簡單、速度快、成本較低、易于小型化和商業化，因此在安檢、無損探傷等領域具有更大的應用價值。

目前太赫茲連續波成像的圖像分辨率與信噪比普遍較低，現有的研究大都集中在去除圖像噪聲［1-3］和增強圖像對比度［4］等方面。圖像分割是實現目標識別的方法之一，分割出目標后可以進一步對其進行特征提取和模式識別。因此，研究太赫茲圖像分割算法，對于提高太赫茲成像安檢應用的疑似目標檢測效果具有重要意義。

在圖像分割領域，KRINIDIS和CHATZIS提出了一種模糊局部信息C均值聚類(fuzzy local information C-means，FLICM)分割算法［5］，它對經典模糊 C 均值聚類算法(fuzzy C-means，FCM)算法［6］加以改進，加入與局部空間信息和灰度信息相關的模糊因子，從而大幅提升了對噪聲的魯棒性。本文中根據太赫茲圖像的性質，增加了先驗信息，改變FLICM算法隸屬度函數的計算方式，以提高該算法在太赫茲圖像處理領域的適用性，實現太赫茲連續波圖像的目標檢測。

1 基于返波管的太赫茲連續波成像

1．1 基于返波管的太赫茲連續波實驗系統

大功率、高效率的THz發射源是THz診斷和成像技術實現的關鍵。目前這些波源大致可以按基于光學方法的太赫茲輻射源和電子學方法的太赫茲輻射源分類。基于電子學方法的太赫茲輻射源包括太赫茲自由電子激光器［7］、納米行波管、返波管(backward wave oscillator，BWO)等。本文實驗中采用Microtech公司的返波管BWO_OV30作為太赫茲輻射源，返波管輸出連續、單頻、線偏振、可調諧的THz波，輸出功率最大值約15mW。通過軟件或電源儀表控制陰極電壓，可在231GHz～375GHz頻段上實現頻率調節，338.288GHz頻率處的太赫茲波輸出功率最大，以該頻率的太赫茲波作為成像系統的波源。

THz探測器是另一項關鍵部件。由于目前THz輻射源普遍功率較低，因此高靈敏度、高信噪比的探測技術顯得很重要。本文中采用Corehent公司的熱釋電探測器(P4-45)。

太赫茲連續波成像系統分為透射式和反射式兩種不同的模式。本文中采用結構簡單、易于搭建的透射式系統，系統主要由返波管、聚乙烯透鏡、斬波器、2維平移臺、熱釋電探測器、鎖相放大器和計算機組成，如圖1所示。太赫茲波由透鏡準直并聚焦于樣品處，并在與太赫茲波束垂直的平面內進行x-y方向的2維掃描。從物體每一點透射的太赫茲波，經另兩個聚乙烯透鏡準直并聚焦到室溫工作的焦熱電探測器上。太赫茲波由斬波器調制，斬波器與探測器的信號同步進入鎖相放大器，得到表征樣品透射特性的電壓信號，經過A/D轉換輸入計算機［8］，計算機控制平移臺移動，逐點獲取每一個像素點的電壓值，再將所有像素點的電壓值量化到0～255的灰度空間，最終得到THz灰度圖像。

Fig.1 Diagram of THz-CW transmission imaging system

1．2 太赫茲成像結果分析

本實驗系統對多種物體進行了成像，選取粘貼在白紙背面的剃須刀片、粘貼在聚丙烯薄片背面的手術刀片、信封里白紙上的鉛筆字、塑料盒里的金屬雙叉等特點較鮮明的成像結果加以分析，如圖2～圖5所示。

Fig.2 Imaging result of the blade by THz-CWa—visible image b—THz image of the blade c—histogram of Fig.2b

Fig.3 Imaging result of the surgical blade by THz-CWa—visible image b—THz image c—histogram of Fig.3b

從圖2到圖5可以較為清晰地辨識出各個物體，但從圖像中也可以看到，各個樣品的太赫茲透射成像都會有不規則的條紋干擾的存在。對于條紋的成因，尚未有明確的定論，相關研究者給出了不同的解釋，大多數研究者認為這是由于干涉的結果［9］，因為激光器本身以及探測器的透鏡都有可能提供相干光束，塑料和白卡紙本身都可能存在一定的反射。

圖2～圖5也顯示了所成圖像的直方圖，可以發現很難從這些直方圖上找到各種物體共有的全局閾值分割方法，尤其是圖4c，直方圖只有一個明顯的峰值。因此，本文中提出了基于FLICM的圖像聚類算法，用于將圖像中的目標提取出來。

Fig.4 Imaging result of the words in pencil by THz-CWa—visible image b—THz image c—histogram of Fig.4b

Fig.5 Imaging result of the double forks in the envelope by THz-CWa—visible image b—THz image c—histogram of Fig.5b

為客觀比較圖像分割精度，對實驗圖像進行手工目標提取，得出參考的分割基準集，如圖6所示。

Fig.6 Ground truths of THz images by manual markinga—ground truth of Fig.2b b—ground truth of Fig.3b c—ground truth ofFig.4b d—ground truth of Fig.5b

2 基于FLICM的太赫茲圖像目標檢測算法

2．1 模糊局部信息C均值聚類算法

模糊聚類是指模糊集里每一個樣本點屬于各個分組的隸屬度是一個取值在［0，1］區間內的數，一個樣本點同時屬于所有的類，但是通過隸屬度的大小來區分其差異。FCM算法對無噪聲圖像的分割比較有效，但由于并未考慮任何空間信息，因此對野值、圖像偽影等噪聲比較敏感。FLICM算法采用的局部窗口模糊因子Gki［5］可以同時將第i個像素為中心的局部窗口內的空間信息和像素特征考慮在一起:

式中，Ni是i的鄰域集合，dij是像素i到像素j的歐氏距離，ukj為像素j相對于聚類k的隸屬度，m是模糊權重因子，xj是像素j的值，vk是第k聚類的聚類中心。Gki表示Ni內像素與聚類中心vk的歐氏距離的加權和。如果Ni內像素屬于第k類的隸屬度較小，則Gki的取值較大。因此，目標函數和隸屬度計算中都加入了Gki。

FLICM 的目標函數［5］為:

式中，N表示圖像像素總個數，c表示聚類的總數目。隸屬度函數［5］為:

聚類中心［5］為:

FLICM的算法步驟為:(1)確定聚類數目c、模糊權重因子m、迭代停止閾值ε和最大迭代次數Tmax;(2)隨機初始化聚類中心，計算隸屬度矩陣U，迭代計數器l=0;(3)使用(4)式計算并更新聚類中心;(4)使用(3)式計算并更新隸屬度矩陣 U;(5)如果＜ε或l＞Tmax，轉至下一步驟，否則轉至步驟(3);(6)根據每個像素的模糊隸屬度完成圖像分割。

這里模糊權重因子m設定為2。在針對前述太赫茲成像實驗結果進行目標檢測時，由于要區分的是物體目標和背景兩類，因此c設定為2，ε設定為0．0001，Tmax設定為 50。

2．2 隸屬度函數的改進

隸屬度函數的設定是運用模糊集合理論解決實際問題的基礎。在FLICM算法中，uki的計算僅考慮樣本值xi到第k個聚類中心的歐氏距離，雖然考慮了局部窗口中像素j對像素i的影響，具備一定的抗噪性能，但針對前述太赫茲圖像時，仍然發現其分割準確性有待提高。主要原因是圖2～圖5所示的太赫茲圖像的目標大都處在直方圖中偏暗的波峰附近，條紋干擾占據了直方圖的很大部分，FLICM的聚類分割沒有結合相關的太赫茲圖像特征加以考慮。

因此，作者設計了先驗權重函數H(xi)。首先，需要找到一個全局二值化閾值作為先驗信息，因此使用最大類間方差法(Otsu法)［10］對圖像進行二值化的分割。Otsu法是一種自適應的閾值確定的方法，它遍歷圖像的灰度范圍，將圖像分成目標和背景兩部分，背景和目標之間的類間方差越大，說明構成圖像的兩部分的差別越大。因此，使類間方差最大的分割點也就是Otsu法給出的全局閾值點。

對于1幅圖像，設當目標與背景的分割灰度閾值為a時，目標類點占圖像概率為p0(a)，均值為u0(a)，背景類點占圖像概率為p1(a)，均值為u1(a)，ut為圖像的總計均值，則建立目標函數式:

式中，Sb(a)就是當分割閾值為a時的類間方差表達式。Otsu算法使得Sb(a)取得全局最大值，當Sb(a)為最大時所對應的a′就是Otsu法得出的二值化全局閾值，如下式所示:

式中，L是圖像的灰度級，一般為256。

但是，僅僅使用Otsu法分割太赫茲圖像的結果也不十分理想。將這個閾值點的值作為先驗信息，當前灰度值較閾值小時，其屬于背景的可能性比較小，在計算其對亮聚類中心的隸屬度時，將其調低更合理;當前灰度值較閾值大時，其屬于目標的可能性也較小，在計算其對暗聚類中心的隸屬度時，將其調低也更合理。另外，灰度值離該閾值點越近的點，被誤分割的概率也越大，反之則越小。

依據上述分析，以 sigmoid函數為基礎，設計H(xi)，得出新的隸屬度函數uki′的計算公式:

其中，

式中，u1i為像素i與相對亮聚類中心v1的隸屬度;u2i為像素i與相對暗聚類中心v2的隸屬度，a′是Otsu法對圖像分割得出的二值化全局閾值，通過q可以調節H(xi)函數曲線的變化率。

算法的目標函數、聚類中心的計算方式和算法迭代步驟與FLICM算法一致。

通過改變q取值，可以調整H(x)曲線的斜率。通過觀察H(xi)與H(2a′-xi)的函數曲線，可以發現q越大，x=a′軸附近的H(x)曲線越陡峭，也就意味著對閾值點附近的像素做更大幅度的隸屬度調整，如圖7所示。

Fig.7 Diagram of H(xi)and H(2a′-xi)under different parameter q when a′=0．2

應用改進的FLICM算法，在不同q取值下對圖3b進行圖像分割，結果如圖8所示。

Fig.8 Segmentation result of the blade by our method under different parameter q

圖9 是q=5時，H(xi)和H(2a′-xi)疊加于圖2b直方圖上的示意圖。

Fig.9 Histogram of Fig.2b with H(xi)and H(2a′-xi)when q=5

觀察圖8中的4個分割結果，發現目標提取結果對q不敏感，不同的q取值都可以較明顯地分割出刀片的輪廓。

接下來，采用概率蘭德指數(probabilistic Rand index，PRI)［11］來對新算法的分割結果進行分割質量評價。PRI算法是一種有監督的圖像分割評價算法，將測試分割結果Stest與標準分割集合進行比較，得到一個數值作為分割結果的評測值。PRI的值反映了實際分割結果與理想分割圖像之間的屬性共生一致性。PRI值范圍在0到1之間，PRI值越接近1表示分割效果越好，越接近0表示分割效果越差。本文中，將圖6所示的人工分割基準作為參考的標準分割集合。

將圖8與圖6a進行比較，得出不同q取值下PRI算法的定量評價結果，如表1所示。

Table 1 PRI quantitative assessment results under different parameter q values shown in Fig.8

從表1可以看出，客觀評價的結果分割結果和觀察結果相似，q由1增長到10，PRI值有一定增加，但總體而言，PRI值對q的變化并不敏感，不同q的分割結果都和人工分割基準比較接近，q=10時，PRI值最接近1，表示分割效果最接近人工分割基準。盡管人工分割基準本身具有一定的主觀性，仍然可以選取q=10作為新算法參量q的經驗值。

另外，分析了q=10時新算法的收斂性，如圖10所示。正如CELIK和LEE所指出的，由于目標函數中包含聚類能量和局部信息能量兩個部分，目標函數并不會一直隨著迭代次數的增加而增加，而是趨向于穩定［12］。本文中提出的算法在迭代次數10次左右目標函數值就已基本穩定。

Fig.10 Convergence of objective function on segmentation of Fig.2b when q=10 by our method

3 算法比較分析

下面將本文中算法與傳統算法對比，選取Otsu［9］，，FLICM 種算法進行比較，結果如圖11～圖14所示。

Fig.11 Comparative result on segmentation of the blade by different methoda—Otsu b—K-means c—FLICM d—our method

Fig.12 Comparative result on segmentation of the surgical blade by different methoda—Otsu b—K-means c—FLICM d—our method

Fig.13 Comparative result on segmentation of the words in pencil by different methoda—Otsu b—K-means c—FLICM d—our method

Fig.14 Comparative result on segmentation of the double forks by different methoda—Otsu b—K-means c—FLICM d—our method

對于以上不同算法的目標檢測比較結果可以看出，作者提出的算法較其它算法能更為清晰地提取出目標的形狀和邊緣，基本濾去了干擾條紋。相比之下，Otsu，K-means，FLICM 3種算法或多或少地都把一些干擾條紋誤判為目標。

對圖11～圖14的各分割結果應用PRI算法進行客觀評價，結果如表2所示。

Table 2 PRI quantitative assessment results on various terahertz images

由表2可知，針對不同太赫茲圖像，本文中算法的目標分割結果都較其它算法更接近人工分割基準，也說明客觀評價與主觀評價基本是一致的。可見本文中算法比幾種經典算法具有更強的太赫茲連續波圖像的目標檢測能力。

在太赫茲波段內，包括本文實驗中使用的樣品在內的大部分常見物品對太赫茲波的色散都不大，因此使用其它頻率的太赫茲發射源進行太赫茲成像時，圖像特性是基本一致的，因此都可以適用于本文中的圖像分割算法。

4 結論

部分物體在基于返波管連續波太赫茲成像系統的成像結果對比度低，存在較強條紋的干擾。本文中針對太赫茲圖像的特點，對模糊局部信息C均值算法做了改進。對成像結果的圖像處理實驗表明，改進后的聚類分割算法較傳統算法可以更好地分割出物體的輪廓，避開條紋干擾，達到了目標檢測和提取的目的。

［1］ HOUSHMAND K，TIZHOOSH H R.Filtering and fusion of THz images for defect detection in composite materials［C］//Proceeding of IEEE International Conference on Fuzzy Systems.New York，USA:IEEE，2008:1301-1305.

［2］ ZOU Y Y，GE Q P，HAN Y，et al.Stripe noise of THz image processing based on frequency filtering［J］.Computer Engineering and Applications，2009，45(17):241-243(in Chinese).

［3］ LI Q，XIA Z W，DING S H，et al.Image denoising of CW THz images by use of non-local mean ［J］.Infrared and Laser Engineering，2012，41(2):517-522(in Chinese).

［4］ XU Y ，HONG Z.Study of multi-scale enhancement algorithm for THz images combining wavelet denoising［J］.Chinese Journal of Sensors and Actuators，2011，24(3):398-401(in Chinese).

［5］ KRINIDIS S，CHATZIS V.A robust fuzzy local information C-means clustering algorithm［J］.IEEE Transactions on Image Processing，2010，19(5):1328-1337.

［6］ BEZDEK J C.Pattern recognition with fuzzy objective function algorithms［M］.New York，USA:Plenum Press，1981:1-13.

［7］ WANG Q，CHI X，LI Q.Imaging theory and development of THz free electron lasers［J］.Laser Technology，2006，30(6):643-646(in Chinese).

［8］ WANG L F，LI J，HONG Z.CW terahertz image system for detecting defects in phenolic foam［J］.Optical Instruments，2010，32(1):13-15(in Chinese).

［9］ YIN Q G.Applied research of markov random field in continuouswave terahrtz image processing［D］.Haerbin:Harbin Institute of Technology，2009:23-25(in Chinese).

［10］ OTSU N.A threshold selection method from gray-level histograms［J］.IEEE Transactions on Systems，Man and Cybernetics，1979，9(1):62-66.

［11］ RANJITH U，CAROLINE P，MARTIAL H.Toward objective evaluation of image segmentation algorithms［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，2007，29(6):929-944.

［12］ CELIK T，LEE H K.Comments on“A robust fuzzy local information C-means clustering algorithm”［J］.IEEE Transactions on Image Processing，2013，22(3):1258-1261.

［13］ MAcQUEEN J.Some methods for classification and analysis of multivariate observations［C］//Procedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability.Berkeley，USA:University of California Press，1967:281-297.