虫蛀麦粒的太赫兹成像检测方法研究

琚新刚 廉飞宇 张 元1, 葛宏义 蒋玉英

(河南粮食作物协同创新中心1, 郑州 450046) (河南省粮食光电探测与控制重点实验室; 河南工业大学2, 郑州 450001) (河南教育学院电路与系统重点学科;河南教育学院3，郑州 450046)

太赫兹成像是当前太赫兹技术应用的一个热点。作为一个最新被利用的波段，通过使用太赫兹波对物体进行辐射,可以获取被测物对太赫兹波的吸收率,进而反映出物体内部空间结构，并且还可以通过使用相位测量来得到被测物的折射率的空间分布阵列。利用这一特点,并与断层扫描原理相结合,可以实现层析成像,并获得被测物的二维或三维结构的信息。与其他波段的成像技术相比,太赫兹成像所得到的图像的分辨率和景深都有明显的增加。近年来，太赫兹成像逐渐被应用到诸多的领域[1]，如材料检测[2]，医学诊断[3-4]，文物保护[5-6]，安检及反恐[7-8]等，成为目前检测领域中的一个研究热点[9]。

太赫兹光子辐射能量很低，只有X光的千分之一，因而特别适合做生物组织的成像，而不会对生物组织造成电离损害[10]。将太赫兹成像技术应用于食品质量安全检测，具有潜在的市场前景[11]。粮食作为食品的主要原料，在储藏期间易受到虫害威胁，特别是蛀食性害虫，由于在粮粒内部生长发育，目前还没有非破坏性的有效检测手段。史卫亚等[12]基于生物光子学理论研究了粮粒内隐蔽性害虫的非破坏性检测方法，取得了较好的效果，但需要专业的隔音室，不能满足实际应用的需要。Guo等[13]使用TeraView TPS3000 观察小麦种子的状态和结构，通过设置不同的回波延时，清晰地观测到了小麦胚部位的状态，但没有涉及麦粒内部虫害状况的观测；Jiang等[14]采用太赫兹成像和化学计量学方法研究了发芽小麦早期的状态，发现小麦在储存的第6小时是阿尔法淀粉酶活跃的高峰期，但没有考虑有虫害小麦的状态变化情况。由于蛀食性害虫是储藏小麦的主要威胁之一[15]，其早期检测十分必要，而太赫兹层析成像可实现物体内部三维成像[16]，因而可对蛀食性害虫提供非破坏性的成像检测，可有效解决蛀食性害虫早期难以检测的问题，本文正是基于这一原理，开展了利用太赫兹成像对蛀食性害虫进行早期检测的研究。

1 设备与材料

实验设备采用的是美国ZOmega公司的Z3型THz-TDS系统。Z3系统以反射模式获得太赫兹成像。THz成像系统包括一个脉冲飞秒太蓝宝石激光器，脉冲宽度和中心波长为100 fs和800 nm。光束经分束镜分为较强的泵浦光和较弱的探测光。泵浦光被频率为1.1 kHz斩波器调制、经透镜聚焦后射向光电导天线砷化镓(GaAs)晶体上，通过光整流效应产生频率范围约为0.2～3 THz的太赫兹脉冲。在试验过程中，系统的信噪比是1 000，谱分辨率好于40 GHz。为了成像，样品被固定在一个X-Y驱动的平台，其移动垂直于进入的太赫兹波束。在测量期间，温度被保持在室温，室内湿度仿真实际应用的条件。

把不同损害程度的麦粒被放在太赫兹成像系统的移动平台上，该平台有着最大的扫描区域50 mm×50 mm。以0.05 s的速率点对点扫描，分辨率为0.1 mm。总的图像获取过程被太赫兹分析和控制软件控制和执行。

一旦扫描过程完成，每个麦粒的太赫兹图像带着一个三维的形式，不仅包括空间信息，而且包括光谱信息的图像被产生、记录和存储。而一个完全的太赫兹波形包含512个时域点，覆盖了一个30 ps的时间范围，相应的频率范围是0～3.5 T，能够被获取在每一个像素位置。这些被获取的太赫兹图像和光谱经由MATLAB和Origin8.5软件分析为特征提取，并且构建模型用于识别麦粒内部的状态。

对于内部有空洞或虫卵的麦粒，将在太赫兹灰度图上呈现出明或暗的区域，但存在着低对比对度、高背景的问题，给分割目标区域带来一定的困难。

2 分析方法

2.1 Otsu法获取最佳分割阈值

阈值法是实现图像目标区域分割的一种常用方法[17]，根据选定的灰度阈值，将图像分割成不同的区域。可以有多种阈值选定准则，其中最大类间方差法(Otsu)[18]是目前比较流行的一种，它以目标区域类与背景区域类之间的方差的最大值作为分割准则，在分割效果上较其他阈值法有明显的优势。Otsu获取最佳阈值的过程如下：

设灰度图的灰度级个数为L，ni为灰度值为i的像素点的数目，则总像素数为：

(1)

像素点的灰度值为i的概率为：

pi=ni/N

(2)

设Th为分割阈值，它将图像分割成R1和R2两个区域，则R1出现的概率为：

(3)

而R2出现的概率为：P2=1-P1

区域R1和R2的灰度均值分别为：

(4)

(5)

整个图像的灰度均值为：

(6)

区域R1和R2的的类间方差d(Th)定义为：

d(Th)=P1(av1-av)2+P2(av2-av)2

(7)

将式(6)带入式(7)，可得

d(Th)=P1P2(av1-av2)2

(8)

最佳阈值Thopt即为d(Th)值最大的阈值，即

d(Thopt)=max[d(Th)]

(9)

利用Z3系统获取样品的太赫兹成像时，将对每一个像素产生一个时域谱和经傅里叶变换后得到的频域谱。利用像素的时域谱成像时，可以有两种方法：一是采用所有记录到的时间轴上的点对应的幅度的平均值；另一种方法是采用峰值，但无论采用哪种方法，得到的幅度都是一个用4个字节表示的浮点数。Z3系统本身自带了一个成像软件，它可以将以浮点数表示的4字节原始数据转换为256个灰度级，从而呈现一个只有256个灰度值的灰度图。对于本文研究的图像，很难直接从256个灰度级的灰度图上分割出目标区域，所以，可以采用从原始数据上应用Otsu方法进行分割，即在原始数据中找到一个最佳分割阈值。对原始数据进行区域分割，最后映射到256级灰度图上。由于原始数据对应的每个像素值是一个4字节的浮点数，因而可认为是对一个灰度级很大的图像进行分割，而Otsu的最佳阈值的求解是一个在整个灰度级的取值范围内的线性求解，这就使得Otsu的最佳阈值求解法的效率很低，甚至是无法完成的。而遗传算法(GA)能够以非线性的方式快速求解最优解，并具有全局寻优能力。将GA算法引入Otsu中的最大类间方差和最佳阈值的查找中，可以大大提高本文太赫兹成像目标区域的分割效率。

2.2 结合GA算法的最佳阈值获取

遗传算法是一种通过模拟自然进化过程搜索最优解的方法[19]。在使用GA算法进行寻优的过程中，一般包括解空间编码、初始化种群、选定适应度函数，交叉和变异，繁衍等操作。本研究中，以类间方差作为适应度函数。在使用遗传算法求解最佳阈值时，适应度函数对算法的收敛性以收敛速度的影响较大。适应度函数的选取与目标函数有关，但在传统遗传算法中，没有考虑目标函数的变化趋势，致使遗传过程中优选的点(染色体)在一个时期内只具有较优的目标数值，但可能都是在较平坦的区域内徘徊，经多代遗传后，目标函数值仍然得不到较明显的优势， 这种情况有可能导致“早熟”现象。在本文的太赫兹成像中，由于小麦本身状态(霉变、牙变、虫蛀等)的复杂性，使得异常区域的识别难以采用恒定的最大类间方差作为适应度函数，因而采用GA算法搜索最佳阈值时容易产生“早熟”现象。为了避免少数个体霸占整个群体而导致的“早熟”现象以及因平均适应度、最佳适应度而使优化过程趋于随机漫游的现象，本文提出了利用目标函数梯度(即变化率)的改进遗传算法。

考虑连续函数的无约束极小化问题,f(x)的适应度函数为：

f(x)=minf(x1,x2,…,xn),xi∈[ai,bi],i=1,2,…,n

(10)

如果需要,还要选取相应的适应值标定方式,引进适应数是为了衡量遗传种群个体的优劣程序，对连续函数求极小问题而言,自然应该要求各个个体适应度函数值越小越好,但是,在一些函数求极值的问题中,为了求函数的极小(即下山),往往不会选择函数值小的搜索点, 而且还必须保证该搜索点的函数值变化率足够大,即只有合适地选择函数值小而且函数值变化率大的点作为下山过程的搜索点,才能使下山过程迈大步,否则,就使下山过程在一块“平坦的高原”上来回漫游,而很难到欲到达的“山谷”,现在把这种思想借到遗传算法中来,在选取适应度函数时,充分考虑目标函数的变化率,对可微函数而言就是利用适应度,提出选用适应度函数:

(11)

(12)

(13)

其中λ∈[0,1]称为控制因子，用以体现适应度函数值和函数值变化率对欲解问题的重要程度。fit(x)为原始适应度函数，fitmin和fitmax分别表示当前代群体中个体适应度的最小和最大值。(vk)i表示向量vk的第i个分量，‖·‖表示某种范数，fit′(x)表示新的适应度函数，ps表示种类规模。

为了实现上述算法，只需在一般的二进制数编码的遗传算法中，在计算应适度函数值时加如下步骤：1)计算群体中fitmin,fitmax,fmin和fmax；2)对所有染色体vk,k=1,2,…,ps，计算函数f(vk)的梯度f(vk)，若fitmax-fitmin<ε或fitmax-fitmin<ε，则fit′(vk)=fit(vk)，否则按照式(11)计算fit′(vk)，其中ε为一小正数。

借鉴上述适应度函数调整的方法，适当调整权值λ，就可以构造出即考虑函数值又充分考虑了函数值变换率的适应度函数。

3 结果与分析

实验结果如图1所示。图1展示了以峰对峰切片模式获取的麦粒2D太赫兹图像。由于回波信号的强弱不同，麦粒的2D太赫兹图像显示出灰度的差别，图中亮度高的地方表示回波信号的强度大。设置不同的回波延时，可以得到不同深度的的回波[13]。回波延时设置时间较短( 如6～11 ps)，只能看到麦粒表层的一些信息，回波延时设置为13～16 ps，可以看到胚的结构，本文将回波延时设置为30～33 ps， 以试图看到麦粒内部的结构。

对两幅小麦粒太赫兹图像(如图1所示)应用本文的方法进行分割，遗传算法的相关参数设置如表1所示。

图1 小麦粒的太赫兹图像

参数名称参数值染色体长度L5种群规模N20交叉概率pc0.7变异概率pm0,01终止进化代数T200

基于表1参数设置，利用遗传算法通过寻优计算得到基于最大类间方差的阈值分割结果如图2所示，其最佳适应度值计划曲线和每一代的最佳阈值进化曲线分别如图3所示。

图2 基于最大类间方差的阈值分割结果

图3 进化曲线

如果采用人工阈值法对图2图像进行分割处理，结果如图4所示。由于阈值是图像分割效果的关键，而基于图像直方图的人工阈值分割法，直方图缺乏典型的双峰结构，使得阈值的选取具有一定的盲目性，如当分割阈值Th=112时，由于雷达图中噪声、干扰的影响，一些不属于目标区域的地方被划分成了目标区域，使得划分出的目标区域比实际区域要大，当分割阈值Th=213时，虽然比较好地抑制了噪声和干扰的影响，但也对目标区域产生了过划分，使得划分出的目标区域比实际区域要小。当Th=168时，划分出的目标区域，与采用本文的方法基本相同，而且与实际目标区域的大小也十分吻合。为了寻找最佳的分割阈值，需要针对每一阈值进行最大类间方差计算，这在图像像素较多、图像灰度级较大的情况下，算法具有很低的效率，而且难以保证分割性能。本文基于改进遗传算法通过对最佳阈值进行寻优，提高了计算最佳阈值的算法效率，取得了较好的分割效果，也为后续的识别和处理工作提供了有力的支持。

图4 基于人工阈值法的分割效果

为了进一步验证本文算法的分割性能，从算法的执行时间和算法对目标区域分割的定量识别两个方面与常规的人工阈值分割方法做了对比。表2是算法执行时间的对比结果。

表2 算法执行时间对比/s

注：运行平台为：Intel Pentium CPU G3240 @3.10GHz/2G/Windows 7。

由表2可知，基于常规GA的Otsu算法，获取最佳阈值的时间是2.321 s，是三种算法中运行速度最快的，其次是本文改进GA的Otsu算法，获取最佳阈值时间是2.342 s，比基于常规GA的Otsu算法稍慢一些，而常规Otsu算法获取最佳阈值的时间最长，为4.794 s，几乎达到了本文算法的2倍以上。

表3 定量识别结果

我们采用链码方法对三种算法分割的目标区域进行了定量识别。链码方法可以对目标区域的面积进行定量计算，其基本思想是用一些首位相连的具有一定方向和长度的线段来表示区域的边界，本文采用8邻域链码(FC8)，用起始点的坐标和8邻域链码值来表示图像中分割出的目标区域的位置和形状。用边界内像素点的数目表示目标区域的面积。通过用链码方法计算三种算法分割出的目标区域的面积与目标区域的实际面积进行对比，评判算法对目标区域的定量分析性能。三种算法对图像中目标区域的定量识别结果如表3和表4所示。

表4 定量识别结果

由表3可以看出，使用基于改进GA的Otsu算法和基于常规GA的Otsu算法，与使用常规Otsu算法，都取得了较好的定量识别效果，采用链码方法获得的目标区域分割的面积与实际面积的误差都不超过10%，且基于GA的Otsu算法误差与实际值相比均不超过5%。而基于改进GA的Otsu算法，其分割精度要稍高于基于常规GA是Otsu算法，且两者的分割精度均高于常规Otsu方法。从识别目标区域的个数来看，如表4所示，基于改进GA的Otsu算法，识别出的目标区域个数与实际个数完全相符，识别率达到100%，基于常规GA的Otsu算法有一个识别错误的目标区域，而常规Otsu算法目标区域的识别错误数达到3个，表明本文的基于GA的Otsu算法不仅在识别时间和识别精度上优于常规Otsu算法，而且在识别稳定性和可靠性上也优于常规Otsu算法。

4 结论

利用太赫兹反射成像技术探测麦粒内部结构，通过采用最大类间方差法，对图像中的空洞区域进行了有效地分割。为了提高最佳分割阈值的寻优效率，引入了遗传算法，同时为了解决由于没有考虑目标函数的变化趋势导致的“早熟”现象，提出了利用目标函数梯度(即变化率)的改进遗传算法。实验结果表明，利用改进遗传算法对最佳分割阈值进行寻优，在保证定量识别精度与传统人工阈值法相比不下降的前提下，寻优的时间大为减少，算法的效率得到了明显提高。需要指出的是，本文对麦粒虫蛀程度的成像与选取的回波延时时间有关，而本文对太赫兹回波延时的设置采用的是经验值，已经可以通过计算分割出的虫蛀空洞图像像素数确定虫蛀的程度。下一步的研究工作主要包含两个方面：一是通过设置不同的回波延时，得到多层的太赫兹成像切片，进行虫蛀空洞的空间构象，进而更精确的判定麦粒虫蛀的程度，二是关注太赫兹成像器件技术的发展，采用可在线化的太赫兹时域光谱技术，实现堆积麦粒虫蛀损害的在线监测。