严志讯，林 晓，童慧鑫，蒋林华

（上海理工大学 光电信息与计算机工程学院，上海 200093）

1 太赫兹图像处理的国内外研究现状

目前市场上主流的安检仪为X射线扫描仪，但是，由于X射线电子能量较高，容易对被检测物质造成电离性伤害，因此，并不适用对人体进行安检扫描。太赫兹是指频率在0.1～10 THz的电磁波，如图1（见下页）所示。类似于X射线，太赫兹波能够穿透各种各样的材料，如塑料、纸张及服装等非金属材料[1-2]。此外，太赫兹波是一种电子能量低的电磁波，这意味着太赫兹辐射将不会对人体组织造成伤害[3-4]。考虑到太赫兹辐射的优点，一种基于太赫兹的成像安检设备正受到越来越多研究人员的关注。

图1 太赫兹频率范围Fig.1 Frequency range of terahertz wave

许多公司和机构已经成功地研制了不同类型的太赫兹成像系统，例如，上海理工大学成功研制了94 GHz被动式太赫兹成像系统[5-9]。太赫兹成像是一种较新的技术，其检测结果极易受到检测灵敏度和噪声的影响。低信噪比的太赫兹图像会给识别危险品造成极大难度。本文提出了一种自动检测太赫兹图像中可疑物体的算法，解决了上述缺点，提高了太赫兹成像系统的实用性。

近年来，许多学者提出了自动分割太赫兹人体图像的算法。Otsu[10]使用Ostu’s阈值方法分割太赫兹图像中藏匿的武器，然而，此方法仅仅考虑图像的直方图并使用单一阈值，很难准确地分割出可疑物。Yeom等[11]采用EM算法和多级检测技术来解决图像信噪比低、小区域对比度较小的问题，导致不能有效地检测藏匿的物品。Shen等[12]提出了一种基于混合高斯密度模型的多级阈值方法，分割太赫兹图像的可疑物。太赫兹图像存在较低的信噪比，导致图像分辨率低，并且直方图中没有明显的波峰或者波谷，所以，上述提到的方法很难直接应用到本文所研究的太赫兹图像之中。因此，减少噪声，借此提升图像的质量是至关重要的。目前有许多优秀的降噪算法可以减少噪声，增强图像质量，如非局部均值滤波器（NL-means filters）[13]、 高斯滤波、形态学滤波器[14]、各项异性滤波[15]等。考虑到太赫兹图像的特征，故采用一种改进的NL-means滤波器来移除太赫兹图像中的噪声并减少NL-means算法的时间复杂度。

如图2所示，由于图中可疑物品和背景有着较为接近的灰度值，导致可疑物体和背景相连通，很难直接分割出可疑物体。因此，提出了一种两级分割模型来解决这一难题。

图3为本文方法的主要步骤。在第一级中，首先通过改进的NL-means降低噪声，并利用主动形状模型（ASM）[16]提取封闭的人体轮廓图，解决可疑物体和背景相连通的情况[17]。在第二级中，通过改进的k-means聚类算法在分割出的人体轮廓中标定可疑物体的位置。为了保护人们的隐私，最终的输出图像仅显示人体的封闭轮廓和可疑物体在人体轮廓中的对应位置与大致形状。

图2 太赫兹安检仪图像Fig.2 Terahertz images showing different dangerous goods in human body area

图3 本文算法的主要步骤Fig.3 Key steps of the proposed model

2 太赫兹图像特性

太赫兹能够穿透衣服并发现隐藏在衣服下面的诸如小刀、手枪等危险品。此外，它还能够有效地识别非金属材料的物体（陶瓷刀、塑料手枪等），然而，这些物品难以被X射线或者金属探测器检测到。

由于本文所研究的太赫兹图像噪声高，即使是人眼也很难准确地分辨出人是否携带隐藏物体。图2展示了太赫兹图像中不同类型的危险品。图4（b）为图4（a）的直方图，分别从3个区域来描述，R1区域是背景区域，R2区域是隐藏物体区域，R3区域是人体区域。根据太赫兹图像的直方图，可以发现其中并没有明显的波峰和波谷，无法通过单一阈值分割可疑区域。

图4 太赫兹图像和太赫兹图像的直方图Fig.4 Terahertz image and the histogram of terahertz image

3 基于优化的NL-means滤波算法

如图5所示（见下页），被动式太赫兹成像系统产生的太赫兹图像有比较高的噪声，本文采用改进的NL-means算法提升太赫兹图像的质量。

NL-means算法2005年首次由Baudes等提出，其基本思想是利用图像的冗余信息降噪，原理为：每个像素的值由图像中与它具有相似邻域结构的像素加权平均得到。有别于双线性滤波器和中值滤波器，NL-means算法利用的是全局信息，它将一个图像块作为一个单元，搜寻相似区域，并计算相关区域的平均值，以此减少图像中的噪声。假设含有噪声的图像为V，去噪后的图像为NL（v）。对于图像 V中的每一个像素v（j），其NL-means算法可以表示为

式中： w (i,j)是一个非负权重，表示像素i和j间的相似度； Z (i)为归一化系数；h为平滑参数，控制高斯函数的衰减程度；d为搜索的邻域半径。

图5 原始太赫兹图像和基于改进的NL-means算法处理后的太赫兹图像Fig.5 Original terahertz images and processed images by the improved NL-means algorithm

采用NL-means算法对自然图像处理时能得到很好的效果[13]。然而，这个算法最大的问题在于时间复杂度较高，以至于不能满足系统实时处理太赫兹图像的要求。因此，本文从两方面改进NL-means算法：首先考虑到非零均值加性噪声，相似邻域应该有着相似的平均灰度值。因此，除了使用固定的预设数量的图像块之外，计算w（i, j）时，还需要考虑像素i和j邻域中的平均灰度值的比 率 。 η1＜(¯(i)/¯(j))＜ η2，(i)和¯ (j)是 像 素i和j邻域中的平均灰度值，η1＜1 和 η2＞1，并且是接近1的常数。另外，可以利用平均梯度来预估邻域间的相似度。像素i邻域的平均梯度可以表示为

对于无噪声图像块，有着不同的平均梯度；对于零均值加性噪声，梯度方向应该较相似。像素i和j的平均梯度的夹角

综上所述，式（2）所述滤波器更新为

式中： σ 是梯度幅值的阈值；σθ是一个阈值[18]。

通过上述优化算法可以大幅度减少NL-means算法的计算量，减小了NL-means算法的时间复杂度。在同一台电脑（Intel I7-4790 3.6 Ghz）上运行改进后的NL-means算法，其平均运行时间从3.3 s减少到0.1 s。

4 从太赫兹图像中提取人体轮廓

如图4（a）所示，部分太赫兹图像中可疑物体的灰度值与图像背景过于相似，使得可疑物和背景区域连接，增加了分割可疑物体的难度。因此，基于得到的降噪图像，需要从太赫兹图像中提取人体区域，以便后续处理。

考虑到人体太赫兹图像的边界不能通过图像的强度信息清晰地显示出来。因此，本文采用ASM[16]提取太赫兹图像中的人体轮廓，而不选用传统的基于单一灰度值的边缘检测方法。

4.1 构建人体模型

一个物体的形状可以通过带有n个关键特征点的形状向量xi描述[19]。

为了构建可变的人体模型，首先选取一些具有代表性的太赫兹图像，并手动标定人体轮廓的关键特征点，如图6所示。通过初步实验，确定最佳关键特征点的数量。若关键特征点的数量太少，则不能够很好地表示人体的姿势形变；若关键点的数量过多，会导致ASM计算量增加，使得时间复杂度增高。本文经过大量的实验，最终选取33个关键特征点绘制人体轮廓。这些点的坐标在形状空间中构成一组样本数据，采用主成分分析（PCA）获得形状空间中所有样本数据的点分布模型（PDM）。

训练集中的任何形状都可以通过训练的平均形状和变化的加权点近似表示。

图6 使用红点标记太赫兹图像中的人体轮廓Fig.6 Human body outline in the terahertz image marked with landmarks

根据式（9），每个元素的姿势与平均形状的差异可以分解为一个主要变化模式的组合。如图7所示，改变第一个形状参数b，手臂逐渐地由外向内收敛。

图7 通过调整参数b, 手臂轮廓由外向内逐渐收敛Fig.7 Arms gradually moving from outside to inside by adjusting parameter b

4.2 应用ASM模型检测人体轮廓

在开始匹配之前，人体轮廓的模板初始化至太赫兹图像中央。然后通过平移、旋转和缩放训练集中给出的形状，经由迭代搜索得出最佳匹配的人体轮廓，这个过程可以由式（10）给出。

这个步骤和检测人体轮廓的过程类似。从太赫兹图像中提取出人体轮廓如图8所示。

图8 通过ASM模型提取的人体轮廓Fig.8 Human body contour extracted by the ASM model

5 基于TSED的k-means聚类算法标定可疑物体

经过上述处理，可以从模糊的太赫兹图像得到清晰的人体轮廓。现在的主要目的是标定可疑物体。

分析观察太赫兹图像，可以发现可疑物体和人体区域之间的差异，主要来自灰度值和纹理信息。在太赫兹图像中，人体区域通常为明亮和平坦的区域，可疑物体呈现出黑暗和不规则的纹理。受文献[20-21]的启发，引入基于TSED（total square euclidean distance）的k-means聚类算法，准确地标记出可疑物区域，增强了算法的鲁棒性。

5.1 基于TSED的k-means的改进的k-means算法

给定一个数据集 X ={x1,x2,···,xn}，k-means聚类目的是将n个数据点划分到k个集合C={c1,c2,···,ck}之中，每个数据点属于距离聚类中心最近的一类。其原理如下：

式中，ui表示聚类中心的均值。

在式（11）中，欧氏距离用于将数据点划分到不同的聚类中心。然而，欧氏距离只适用于凸形或者球形分布的数据集。本文采用基于Total Squared Euclidean Distance（TSED）[20]的改进的 kmeans聚类，以提升其处理复杂分布数据的鲁棒性。TSED的表达式可以描述为δ（x1，x2）。

与欧式距离的平方 d (x1,x2)= ∥x1-x2∥2相比，TSED距离计算得到的是一个不会随坐标轴旋转而变换的正交距离，所以，δ（x1，x2）有着更好的鲁棒性。Liu等[21]提出了一种新方法来估计聚类中心，称为T-center。T-center是权重均值，这种权重对噪声和异常值不敏感。受其启发，式（11）可以改写为

其中

式中， wj为第j个像素跟第i个像素的高斯权重。

5.2 基于G-features，W-Features和S-Features检测可疑物体

由上述分析的太赫兹图像特征可知，可疑物体与正常人体区域的主要区别在于灰度值和纹理。本文采用Gabor滤波器提取纹理信息（W-features），同时结合灰度信息和位置信息，应用于改进的k-means聚类算法。Gabor滤波器在提取纹理特征方面具有很大的优势。对于一幅二维图像I（x，y），Gabor滤波器可以方便地提取图像在各个尺度和方向上的纹理信息，其表达式为

式中： h (x,y)表 示图像的纹理信息 δu和 δv高斯核两个方向上的尺度；u，v分别表示频域变量。

选取5个尺度和8个角度(θ=0，15，30，90，120，150，180 °)。W-features通过 40 个 Gabor滤波器获得，其卷积结果如图9所示。

图9 对图4（b）使用不同Gabor滤波器的结果Fig.9 Convolving results of Fig. 4(b) using different Gabor filters

Gabor滤波处理之后，再用一个简单的低通滤波器平滑Gabor的幅度信息。Jain等[22]利用非线性sigmoidal函数从滤波结果中获取Gabor特征。

将特征用于k-means，k-means聚类算法可改写为

式中：Gt为灰度特征； α 为常量；I为图像；g(σu,σv)是标准差 σu和 σu的二维高斯核。

构建一个特征空间K，它包含了Gabor特征Wi,j，灰度特征Gt，位置信息S。特征空间可以描述为

图10展示了从人体区域检测可疑物体的结果。通过引入TSED优化的k-means算法，可以很好地检测可疑物体。

图10 不同的聚类算法比较Fig.10 Comparison of different clustering algorithms

6 实验与分析

6.1 可疑物体检测与分割

在图10中，第一行中的原始太赫兹图像展示了人体携带的各种各样的危险品，第二行是手工标定的真值图，第三行和第四行分别为现有传统的聚类算法得到的分割结果和本文算法得到的效果图。

6.2 定量分析分割的准确性

本文采用准确率、回归率及F均值这3个指标来判断算法的有效性，如图11所示。

图11 Pr-R曲线与F均值柱状图Fig.11 Precision-recall curve and F-measure histogram

准确率Pr和回归率R的定义为

式中：M表示分割后图像前景的像素个数；G表示真值图中前景的像素个数。

F均值Fβ本质上是准确率以及回归率的平均数。

式中，β2一般设定为0.3。

7 结 论

综上所述，考虑到原始太赫兹图像具有较低的信噪比，采用改进的NL-means算法改善图像质量。然后采用基于TSED的k-means算法，来自动标定可疑物体。实验结果表明，本文的算法可以准确地检测出可疑物体，大大提升了太赫兹安检系统的实用性。