基于模糊局部信息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.