一种基于遗传算法的X射线CT多相流检测方法

摘要

一种基于遗传算法的X射线CT多相流检测方法，该方法通过结合被检测多相流体系的相分布或介质分布规律构造具体遗传算法，可以实现对有限角度X射线投影数据的流场重构，其具体步骤包括：1)对被检测体系进行简化处理；2)对被检测区域进行几何剖分；3)X射线投影数据的采集；4)X射线投影方程组的离散化；5)构造具体遗传算法应用于图像重构过程；6)将数值解转化为被检测区域的相/介质分布图。本发明一方面能更好地完成多相流快速成像过程，另外一方面，结合体系物理特征构造具体遗传算法，实现了约束优化的求解过程，具有较强的抗测量噪声干扰的能力，有效控制了图像重构过程的不适定性，明显减弱了测量误差对重构结果的影响。

说明书

技术领域

本发明涉及一种基于遗传算法的X射线CT多相流检测方法，属于多相流测量技术领域。

背景技术

医学CT技术已经成熟地应用于人体的断面和立体成像研究及病情诊断。至今，医学CT已经进入以电子束CT机为代表的第五代设备。五代医学CT机的单投影和断层扫描时间如下表所示。

五代医学CT机的单投影和断层扫描时间

 扫描时间  第一代  第二代  第三代  第四代  第五代 单投影(s)  1  1  ～0.001  ～0.25  ＜0.0004 断层扫描(s)  ～200  ～18  5  1  0.01

尽管医学CT在快速检测方面已经取得重要的进步，但是医学CT主要针对人体，设备庞大，造价极其昂贵，要求大量的扫描角度以达到完备投影数据从而确保精确的图像重构。这些因素都限制了医学CT直接用于工业过程测量和实验室的科研应用。另一方面，除了第五代医学CT具有较快的断层扫描速度，其它的CT技术一般用于获得时均的密度场信息。

X射线CT在理论上已有相对成熟和完备的图像重构方法，其图像重构算法一般要求很高的投影数据完备性，以保证对被检测物体内部信息的准确重现，这一点也同时限制了现有的医学CT的时间分辨率。欲实现快速CT测量，则需要昂贵的设备，提供足够多的X射线源和检测器。

随着医学CT原理越来越多地应用于工业工程，医学CT的时间分辨率的限制也逐渐地突出。

在石油、化工、冶金、动力、能源、医药、食品等领域中，多相流动过程普遍存在，如气-液、气-固、液-固、气-液-固等体系。由于存在界面效应和相间作用，多相流动过程极其复杂。要建立过程模型，并进行过程的预测、设计和控制，首先要解决多相流动过程的测量技术问题。科学研究和生产技术的迅速发展要求开发和更新具有高时间分辨率、高空间分辨率和适宜价位的多相流检测设备，以应用于多相流场中二维截面或者三维空间内的瞬态温度场分布、速度场分布和相组成分布等信息的获取，甚至达到原位在线检测的效果(也即流场信息瞬态“冻结”并且快速重构)。

过程层析成像(PT)技术是医学CT技术应用于工业过程而产生的。由于其对流场无干扰、可获取过程参数的二维/三维分布等特性而备受关注，目前主要应用来测量各相相含率分布和流型判断等。PT技术的基本原理为：在被检测流场四周按一定规律排列若干个传感器，传感器阵列用非接触或非侵入方式获取被检测对象的流场信息，然后利用图像重建算法计算出断面流体分布并显示。PT系统有多种获取信息的形式，按传感器的检测原理，可以分为电磁辐射式、声学式和电测式三大类，分别以X射线CT(X-CT)、超声CT(UCT)和电容CT(ECT)为典型代表。目前的PT技术都还在发展过程当中，有各自的局限性。

X-CT的主要优点是在投影数据较完备的情况下图像重建运算较为简单，成像精度较高(决定于射线源和投影数，最高可达100纳米)；主要缺点是如果考虑经济性因素而采用第三、四代CT扫描方式，每断层射线源旋转过程耗时达到1秒甚至更高的量级，严重影响了瞬态信息的获取效果。X-CT的高空间分辨率的特点对实现多相流的流场精细测量意义重大，保证其优异的空间分辨率又同时提高其时间分辨率是目前X-CT发展的重点。

同时，CT重构问题可以归结为基于第一类积分算子的不适定反问题的求解问题，其不适定性体现在：重建结果敏感于测量投影的微小变化，难以得到稳定解。对于ECT，检测电极间电容的微小变化，会不可避免地引入测量噪声，这些噪声严重影响了重建图像的稳定性，造成了成像系统的失真，降低了重建结果的可靠性。陈德运等(陈德运，郑贵滨，于晓洋，孙立镌.基于遗传算法电容层析成像图像重建算法的研究，电机与控制学报，2003，7(3)：207-211.)引入遗传算法，尝试结合流型特点，以遗传算法固有的高随机性、自适应性和复杂非线性问题全局最优解搜索高鲁棒性来克服图像重建不适定性，取得了较好的效果。对于X-CT，当投影数据不完备(如投影角度较少)时，一般将投影方程组离散化为线性方程组之后采用迭代法求解，也存在类似于ECT的不适定性的问题。

遗传算法是一种基于生物自然选择与遗传机理的随机搜索算法，它相对于传统的基于目标函数最小化的搜索算法相比，具有以下特点：遗传算法的搜索始于解的一个种群，而不是传统优化方法的单个解；遗传算法只使用适值函数，而不使用目标函数的梯度或高阶导数等辅助条件；遗传算法采用随机的，而不是确定的状态转移规则。以上特点使得遗传算法能够非常有效地进行概率意义下的全局搜索，从而有能力收敛于问题的最优解。遗传算法的操作通常包括：编码与初始种群的生成；适应度评估及判断收敛；选择操作，交叉操作和变异操作。

遗传算法从任一初始种群出发，通过选择(使优秀个体有更多机会传给子代)，交叉(体现优秀个体间的信息交换)，变异(引入新的个体，保持种群的多样性)操作种群一代一代地进化到搜索空间中最优点附近，直至收敛到最优解点。

另外，考虑到多相流体系的物理特征，如气-液体系是典型的局部浓度为0-1的体系，气浆体系为近0-1体系(浆态相视为拟均相)、气固两相鼓泡过程也可近似为0-1体系(即乳相和泡相)；同时，多相流中的大尺度结构的动态行为(如气液鼓泡动力学等)往往对流动和反应器行为起控制作用。针对二值或近似二值的多相流体系，充分结合体系的物理特征，简化X-CT技术在多相流测量应用过程中的硬件和软件要求，对化工多相流动的深层次研究有重要的科学意义和工程价值。

发明内容

本发明的目的是提供一种基于遗传算法的X射线CT多相流检测方法。用这种方法通过结合被检测多相流体系的相分布或介质分布规律构造具体遗传算法，可以实现对有限角度(如3～24个角度的扇形束)X射线投影数据的流场重构，相比于传统CT技术的完备投影数据的采集要求(数百个角度的扇形束)能更好地完成快速成像过程。

因此，为了达到上述目的，本发明提供了一种基于遗传算法的X射线CT多相流检测方法，该方法按如下步骤进行：

1)对被检测体系进行简化处理

结合被检测多相流体系的相分布或介质分布规律，使之简化为二值或近似二值体系，整个空间的相/介质分布非0即1；

2)对被检测区域进行几何剖分

选择多相流体系的某一待测二维截面作为被检测区域，将被检测区域剖分成若干个单元子区域，设定各单元子区域的内部介质分布均匀，则所有单元子区域的介质分布赋值构成一个图像向量；

3)X射线投影数据的采集

利用X射线源、X射线检测器和数据采集控制设备采集有限角度下的扇形束X射线投影数据，不同角度的投影数据采集通过单个X射线源的旋转-发射交替进行，或者多个X射线源由电子开关控制依次发射完成；

4)X射线投影方程组的离散化

计算表征各条X射线与图像向量各分量之间相互关系的加权系数，将投影方程组离散化为线性方程组Ax＝b，示意图见图1，其中：A表示关系系数矩阵，x表示图像向量，b表示投影向量，同时选择图像重构最优准则；所述的加权系数基于以下四种射线定义方式之一确定：

①忽略射线宽度，当射线穿过一单元子区域时，取该单元子区域对应的加权系数为1，否则为0；

②忽略射线宽度，当射线穿过一单元子区域时，取该单元子区域截得射线的长度作为表征其投影贡献的加权系数；

③把射线看成具有一定宽度的粗线，取该射线覆盖各单元子区域的面积作为表征其投影贡献的加权系数；

④把射线看成具有一定宽度的粗线，认为单元子区域的赋值集中于单元子区域的中心，射线经过该中心时该单元子区域对应的加权系数为1，否则为0；

5)构造具体遗传算法应用于图像重构过程

结合被检测多相流体系的相分布或介质分布规律构造具体遗传算法，设计基于流动物理模型的选择算子、交叉算子和变异算子，并引入多种群同时进化的计算方法，实现对图像向量的有效重构计算，所述遗传算法的具体构造步骤如下：

①选择“排序选择”之选择算子和“单点交叉”之交叉算子；

②引入基于节点的变异算子：执行每代遗传操作时等概率随机取空间若干节点，逐一统计节点所涉及区域内单元子区域0-1分布情况，记N(0)为0值单元子区域数目，N(1)为1值单元子区域数目；当N(0)、N(1)任一为0时，所涉及区域内各单元子区域变异率设定为0；当N(0)、N(1)均不为0时，所涉及区域内0、1值各单元子区域的变异率分别设定为：

$>>p>>(>0>)>>=>>>N>>(>1>)>>>>>(>N>>(>0>)>>+>N>>(>1>)>>)>>·>N>>(>0>)>>>>>s>$

$>>p>>(>1>)>>=>>>N>>(>0>)>>>>>(>N>>(>0>)>>+>N>>(>1>)>>)>>·>N>>(>1>)>>>>;>>s>$

6)将数值解转化为被检测区域的相/介质分布图

将所求得的图像向量数值解转化为被检测截面的相/介质分布图并输出。

在上述技术方案中，本发明的技术特征还在于：所述步骤2)中的被检测区域的几何剖分所采用的网格包括方形网格和有限元三角网格；步骤3)中的所述的有限角度包括3～36个角度；步骤4)中所述的最优准则包括以下四种：最小二乘方准则，最大均匀性准则及平滑准则，最大熵准则，贝叶斯准则；其中最小二乘方准则为主则，另外三者为附则，根据具体问题需要选择某些附则与主则结合使用；步骤5)中所述的遗传算法包括其他方式的改进遗传算法和模糊遗传算法；步骤6)中所述的相/介质分布包括相含率分布、介质分布和流型。

本发明于现有技术相比，具有以下实质性特点和显著效果。一方面，继承了X射线CT技术空间分辨能力高的优点，通过引入体系先验认识以降低X射线投影数据完备性要求，进而提高了时间分辨能力，可达到快速采集投影数据的效果，拓展了X射线CT技术的适应性和应用范围。另外一方面，结合体系物理特征构造具体遗传算法，实现了约束优化的求解过程，具有较强的抗测量噪声干扰的能力，有效控制了图像重构过程的不适定性，明显减弱了测量误差对重构结果的影响。

附图说明

图1为投影方程离散化示意图。

图2为投影数据采集示意图。

图3为被检测区域有限元三角网格剖分示意图。

图4为数值仿真实验中不同投影角度数下本发明的图像重构结果。

图5为数值仿真实验中不同投影角度数下FBP法的图像重构结果。

图6为静态实验中不同投影角度数下本发明的图像重构结果。

图7为静态实验中不同投影角度数下FBP法的图像重构结果。

图8为静态实验X射线源投影值在某三个投影角度下与理论值的比较。

图9为不同网格剖分下复杂多泡体系应用本发明的图像重构结果。

具体实施方式

以下结合实例，对本发明的具体实施作进一步的描述。

被检测体系如图2所示，待测截面为一直径为90mm的圆面，也即最大的圆所在区域，以其圆心为原点，x轴正向水平朝右，y轴正向竖直向上，建立直角坐标系，单位为mm。截面由空气-水两相组成：空气相分布于四个圆内，圆心位置分别为(-24，10)、(22，15)、(2，-2)和(-2，-29)，对应的直径分别为20、16、16和20mm；其余区域均为水。

1)对被检测体系进行简化处理

将空气-水两相体系简化为0-1二值体系，0代表空气，1代表水。

2)对被检测区域进行几何剖分

进行方形网格剖分如图1所示，将被检测区域放置到圆面的外切矩形中去，矩形区域被均匀划分为若干个方形单元子区域。进行有限元网格剖分，如采用三角网格时，网格划分的粗细决定了空间分辨能力的高低，越精细的网格越能准确描述截面相/介质分布。在算法实现过程中，根据空间分辨率的要求，决定是否对网格进行“变焦操作”，也即对每个单元子区域进行再剖分，如连接三角形的三边中点将其剖分为四个全等的三角形。

如图3所示，对待测区域进行一次三角网格剖分，共产生128个单元子区域，调整节点位置，使三角形尽量接近正三角形。在上述剖分下的重构结果基础之上进行两次“变焦操作”，则可得到2048个单元子区域下的重构结果。

3)X射线投影数据的采集

利用X射线源、X射线检测器和数据采集控制设备采集有限角度(如3～36个角度的扇形束)下的扇形束X射线投影数据，不同角度的投影数据采集通过单个X射线源的旋转-发射交替进行，或者多个X射线源由电子开关控制依次发射完成。所需射线源的能量和数量，需要由事先对该体系进行估算而确定。

投影数据采集分两种形式：数值仿真实验和静态实验。前者由已知截面相/介质分布精确计算出投影向量；后者采用一台工作电压为150kV的闪光X光机和一个X射线平板检测器，通过旋转静态模型，分步实现不同角度的投影数据采集。根据数值仿真实验或静态实验等间距安置若干组X射线源和检测器，而且二者与原点距离保留固定值(分别为180mm和130mm)。每个方向发射出的40度张角的X射线扇形面上共有196条射线，X射线宽度为1mm，在探测器上的投影距离为等距。在数值仿真实验和静态实验过程中均取了24个角度(每隔15度完成一次投影)，重构过程中分别取其中的3、4、6、8、12个投影角度或全部投影角度数据进行图像重构。

4)X射线投影方程组的离散化

迭代重建算法和反投影重建算法最大的区别在于前者一开始就将截面图像f(r，θ)离散化，如图1所示，每个单元子区域内部的离散值x(k)＝f(r_k，θ_k)＝pixel(i，j)，而且为常数，该离散值表示灰度或密度，其中x(k)的表达方式适用于有限元三角网格剖分方式，pixel(i，j)适用于方格剖分方式。图像向量x为K维列向量。

第p条射线的总投影值为 $>>b>>(>p>)>>=>>\Sigma >k>{{}_{}}_{}$_pk为第k个单元子区域对总投影值贡献的加权系数。若投影数为P，则得到P个离散化的投影方程，再记A＝{r_pk}_P×K，进一步可得到线性方程组Ax＝b。其中A在设置数据采集系统和选择投影射线的时候就已经确定，b也就是对应的投影值向量。

在本实施例中，图像向量x的每个分量只能取0或1，投影向量的每个分量代表每条对应射线所穿透的所有水的厚度。图像向量维度K＝128或2048；射线的总数P＝196*NA，其中NA为扇形束投影的数目。

由于该方程具有较强的不适定性，解向量不在连续空间域内，使得问题求解难度很大，我们将之转化为约束优化问题 $>>>min>x>{}_{}$

5)构造具体遗传算法应用于图像重构过程

对于这一问题，遗传算法的作用是不断搜索、取舍x，使之越来越能满足缩小Ax和b之间距离的要求，从而避开了直接求解Ax＝b过程中产生的病态效应。

通用遗传算法对这一特定问题在收敛性、有效性和效率等方面都有局限，必须结合被检测多相流体系的相分布或介质分布规律构造具体遗传算法，设计基于流动物理模型的选择算子、交叉算子和变异算子，并引入多种群同时进化的计算方法，实现对图像向量的有效重构计算，所述遗传算法的具体构造步骤如下：

①选择“排序选择”之选择算子和“单点交叉”之交叉算子，其中“单点交叉”的交叉率设定为0.8；

②引入基于节点的变异算子：执行每代遗传操作时等概率随机取空间若干节点，逐一统计节点所涉及区域内单元子区域0-1分布情况，记N(0)为0值单元子区域数目，N(1)为1值单元子区域数目；当N(0)、N(1)任一为0时，所涉及区域内各单元子区域变异率设定为0；当N(0)、N(1)均不为0时，所涉及区域内0、1值各单元子区域的变异率分别设定为：

$>>p>>(>0>)>>=>>>N>>(>1>)>>>>>(>N>>(>0>)>>+>N>>(>1>)>>)>>·>N>>(>0>)>>>>>s>$

$>>p>>(>1>)>>=>>>N>>(>0>)>>>>>(>N>>(>0>)>>+>N>>(>1>)>>)>>·>N>>(>1>)>>>>.>>s>$

实施“变焦操作”时，将粗网格下的重构结果适当处理，作为细网格下的迭代初值，有利于重构过程更可靠、更有效地收敛。

6)将数值解转化为被检测区域的相/介质分布图

遗传算法有两个终止准则，也即繁衍代数和代间适值偏差，前者取决于问题本身性质以及选择、交叉、变异算子的设置，后者取决于适值函数的设置。为兼顾解的可靠性和运算的效率，加上本身为随机优化过程，必须进行多次重构，确定更优的操作设置和终止准则，得到真正稳定的可靠解。之后，将所求得的图像向量数值解转化为被检测截面的相/介质分布图并输出。

在本例中，根据重构过程的复杂程度设置繁衍代数为500～10000，代间适值偏差取10^-6，保证重构过程得到稳定的可靠解。再将图像向量转化为被检测截面的相/介质分布图并输出，如图4、6系列所示(K＝2048)，图中还给出了气泡的期望位置以作参照，a、b、c、d、e和f分别表示3、4、6、8、12和24个角度投影数据的重构结果。

经过数值仿真实验和静态实验验证，可以发现经过改进的遗传算法的确能有效地控制图像重构的不适定性，在较大的噪声干扰下，还能得到稳定的可靠解。为说明本发明的优越性，使用经典CT重构算法之滤波反投影重构算法(FBP)计算同一算例，重构结果如图5、7系列所示，a、b、c、d、e和f分别表示3、4、6、8、12和24个角度投影数据的重构结果，其中图7中a、b二图由于缺乏物理意义，略去。

如图4，由数值仿真实验的结果可知，从3到24个投影角度，本发明都可以重构出气泡(空气相)的位置和大致形状，并且随着投影角度数增加，相边界的描述不断得到改善，重构结果相当理想。如图5，考察FBP法的重构结果，在8个投影角度下，气泡的位置和形状仍无法识别判断，至少需要12个投影角度配合人脑智能才能识别出气泡的位置和形状，而重构图像在被检测区域内外的伪影均相当严重。对比二者的结果可知，FBP法对于有限角度数据重构是无能为力的。本发明的算法耦合了体系的物理特征，显示出了解决二值体系下有限角度数据重构问题所具有的较大优势。

数值仿真实验说明了算法的潜在优势，而实际测量过程中所采集的投影数据将伴有系统误差和随机误差，本方法的可行性必须经过物理实验的验证。受限于实验条件，我们的静态实验设备的精度不高，远不如医学CT装置的投影精度，某三个角度扇形束投影的实测值和理论值如图8所示。虽然存在不小的偏差，但重构结果仍然非常引人注目。如图6，由静态实验的结果可知，即使在三个角度下，本发明也能给出气泡大致的位置和形状，并且随着投影角度数增加，相边界的描述也能不断得到改善，重构结果趋于理想；同时可以注意到，较多投影数下重构结果和期望结果存在固定的偏差，这应该归咎于实验系统的非理想性。如图7，考察FBP法的重构结果，即使在24个投影角度下，重构结果还是不理想，气泡的边界处有较大偏差；而在12个投影角度下，重构结果已经非常模糊。对比二者的结果可知，本发明在解决二值体系下有限角度数据重构问题方面具有较大的优越性。

针对上述体系，本发明在数值仿真实验和静态实验中都得到了比较理想的重构效果，但一般体系中气泡的数量可能远远多于4个，而且气泡的尺寸也会比较小。将本发明应用于更复杂的包含10个气泡的体系，采用4个投影角度，在不同的网格剖分下数值仿真结果如图9所示。由结果可知，本发明成功地确定了所有气泡的位置和大致形状。

本发明的另外一个独特之处在于：只要计算机硬件条件许可，“变焦操作”可以不断进行，重构结果可以得到不断的改进，更趋理想。

因此，本发明提出的基于遗传算法的X射线CT多相流检测方法是确实可行的，具有其独到的优越性。