傅里叶功率谱在尘肺阴影密集度判读中的应用

罗海峰 翟荣存

（1.中国科学技术大学，安徽 合肥 230026；2.安徽工业职业技术学院，安徽 铜陵 244000；3.铜陵市立医院，安徽 铜陵 244000）

一、引言

尘肺病是我国目前最严重的职业病之一，它是由于长期工作在生成性粉尘的环境下，使人的肺组织呈现纤维化的疾病。在检查方面，目前各大医院对尘肺病的检查基本都采用直接数字X线摄影DR系统。国家卫生部2009年3月发布的《尘肺病诊断标准GBZ70-2009》将尘肺的严重程度划分为五期，分别是无尘肺、观察对象以及壹期、贰期和叁期尘肺[1]。为充分发挥数字化图像优势，提高诊断工作效率，同时，尽量避免因人为因素所造成的错判，因此，实现对尘肺病的智能判读分期具有重要的实用价值。

DR胸片的各个肺区尘肺阴影密集度是尘肺病分期判读的重要依据，而胸片的纹理特征又是尘肺阴影密集度判读的最重要特征。经观察发现，不同尘肺阴影密集度的胸片其纹理结构、分布走向、光点粗细有一定差别。当采用傅里叶变换将图像从空间域变换到频域，其频谱自然也会随着纹理结构的变化而发生变化，也会导致其傅里叶功率谱发生变化。因此，本文首先邀请医学影像专家对不同期别尘肺的DR胸片样本进行分期判读，并对各分区阴影密集度进行记录。再将DR胸片影像进行傅里叶变换，获得不同频段的傅里叶功率谱P（u,v)=|F(u,v)|2占总的功率谱能量的百分比特征值，然后选取区分度较大的3个特征值组成描述纹理的特征值向量并运用BP神经网络进行分类测试[2]，得出结果与影像专家判读吻合度最高达83.3%，效果良好。

二、尘肺阴影密集度及分级判读

DR胸片尘肺阴影密集度的划分是根据各个分区内小阴影的大小、形状和数量，共分为4级，分别为0、1、2、3级，并对其进行字母和数字标识。其中每个级别又分为3个小级，如果通过观察发现阴影与相应的标准片基本一致，则在“/”的上方和下方记录该级别，如判定为1级密集度，则标识为：“1/1”，但如果发现密集度比1级低一些，也可以考虑将其判定为“0”级，则标识为“1/0”；如果比 1级高一些，也可以考虑将其判定为2级，则标识为“1/2”。这样一来，4个级别的阴影密集度可细分为12个小的级别：分别表示为：0/-、0/0、0/1；1/0、1/1、1/2；2/1、2/2、2/3；3/2、3/3、3/+等。

三、离散傅里叶变换及傅里叶功率谱

离散傅里叶变换可以将二维离散图像从空间域变换到频域，生成二维的离散傅里叶频谱，其变换表示式为：

式中，u=0,1,2…，M-1，v=0,1,2…，N-1

其中f(x,y)是数字图像，x,y是空间域中图像的横、纵坐标轴；F(u,y)是频谱，通常F(u,y)是两个实频率变量u和v的复值函数，u、v频率分别对应x轴和y轴。傅里叶变换后的幅谱、相位谱和能量谱表达式分别为：

其中 R(u,v)，I（u,v)分别表示 F(u,v)的实部和虚部。P(u,v)=|F(u,v)|2是傅里叶功率谱。通过对一般纹理图像的实验，容易发现：纹理粗细不同的图像，在傅里叶功率谱分布上是不同的，粗纹理在低频段的功率分布相对较大，而细纹理在低频段的功率分布相对较小[3]。

四、BP神经网络

BP神经网络的拓扑结构如图1所示，它是一种多层的前馈神经网络,该网络对信号处理的主要特点是信号从输入层开始传递，经隐含层处理，逐层向前传递，一直传递到输出层。对于输出层来说，如果没有获得指定的期望值输出，则将期望值与实际输出之间形成的误差信号通过网络反向传播，并将误差信号作为网络权值和阈值实时调整的参数，通过不断循环迭代运算，使得网络的实际输出不断逼近期望输出，直到达到迭代运算次数或满足误差条件。

图1 BP神经网络拓扑结构图

图中X1,X2……Xi是BP神经网络的输入，y1,y2……yo是BP神经网络的实际输出,wij是输入层I与隐含层J之间权值系数，wjo是隐含层J与输出层O之间的权值系数。根据BP神经网络的结构特点，在网络训练过程中，当数据从输入层I输入，经过网络运算到输出层O，将其输出结果与期望值比较，如果其误差值超过指定值，则修正网络的权值系数[4]，如此循环，直到误差值达到指定值或循环次数达到指定次数。

五、实验结果与分析

（一）DR胸片影像预处理

因为通过X线DR系统产生的胸片影像为256阶灰度模式，为增强正常胸片纹理在受到尘肺病变影响后的纹理效果，增大不同阴影密集度影像通过傅里叶变换后在不同频段的傅里叶功率谱强度分布区分度，可以首先对影像进行预处理，即：灰度值降阶。其效果如图2、图3、图4所示。

图2 密集度为0的原始影像和降阶后影像

图3 密集度为1的原始影像和降阶后影像

图4 密集度为2的原始影像和降阶后影像

（二）傅里叶变换及傅里叶功率谱特征值提取

对降阶后的灰度影像进行傅里叶变换，分别得到的傅里叶频谱如图5所示。将整个频段从低频到高频分为10个分频段，统计各个不同频段的傅里叶功率谱占整副图功率谱能量的百分比，其分布情况如图6所示。从图中可以看出，不同阴影密集度的傅里叶功率谱在1、2、3频段的百分比区分度最大，因此，可考虑将此作为纹理的特征值进行提取[5]。

图5 0、1、2密集度影像的傅里叶频谱图像

图6 不同频段傅里叶功率谱强度分布

安徽铜陵地区是尘肺病发生的典型区域。本文提取了97份有效样本，抽取21副经过医学影像专家判读的不同期别尘肺胸片，对其右上部分进行如上方法特征值提取，分别得到傅里叶功率谱前3频段的数据分布如表1所示。

表1 不同阴影密集度傅里叶功率谱能量百分比

（三）BP网络训练与测试

设置BP神经网络相关参数，将输入层节点数设置为3，输出层节点数设置为3，隐含层节点数设置为6，最小误差值设置0.08，迭代最大循环次数为50000次，权值修正的速率设置为0.3，用小的随机数初始化权值W和偏差B，运用MATLAB进行仿真实验。

训练与测试方式一：将表1中1-5、8-12、15-19号样本的傅里叶功率谱分布百分比值和输出期望值作为训练样本，得到的训练误差函数值变化情况如图7所示。

图7 BP神经网络误差函数

图7横坐标显示迭代次数，纵坐标显示网络输出值与期望值之间的误差值，黑色横线为设置的最小误差值，黑色曲线随着迭代次数增加不断逼近最小误差值，从图中可以看出，训练过程在迭代次数时结束未能达到最小误差值。

然后将 6-7、13-14、20-21号样本作为测试样本，并对输出值按式（5）函数进行处理（其中Y为网络输出值，i为样本序号），得出各测试样本的密集度值，结果如表2所示。从表2中可以看出只有14号样本未能与预期值保持一致，正确率为83.3%，对于该样本未能获得预期结果，其主要原因是：该样本在医生判读时是介于1密集度与2密集度之间，生成的纹理功率谱特征值与2密集度纹理特征更加接近，因此网络将其判读为2[6]。该问题需要后期进一步加强研究予以解决。

表2 测试样本的输出值与分类结果

训练与测试方式二：为增强傅里叶功率谱特征值的可信度，采用Bootstrap法随机从小样本集中抽取样本进行训练和测试，具体做法是：从表1中不同级别阴影密集度的7个样本分别抽取5组和剩余2组数据作为训练和测试样本，并采用上述BP神经网络进行训练和测试，共进行20次。最终实验结果如表3所示，其中最高吻合度为83.33%，最低吻合度为50.00%，平均吻合度为64.17%，方差为12.10%，充分体现特征值提取和数据分类的有效性。

表3 采用bootstrap法测试结果

六、结语

本文通过对不同尘肺期别和阴影密集度DR胸片样本进行试验，并采用bootstrap法进行测试，从其结果可以看出：通过傅里叶变换生成的傅里叶功率谱在不同频段的强度分布所组成的3个特征值能够对不同密集度纹理特征进行有效描述，并采用神经网络进行分类，能够实现对密集度的有效判读，且效果良好。