本发明属于电阻抗断层成像(Electrical Impedance Tomography，EIT)领域，特别涉及一种颅脑EIT局部最优正则化参数选取的方法。
背景技术：
电阻抗断层成像技术(Electrical Impedance Tomography，EIT)通过安装在物体表面的电极对物体有规律地施加激励，若物体内部存在阻抗变化，则会引起表面测量电极电位的变化，基于测量电极电位变化，结合相应重建算法可以获得物体内部阻抗变化图像。中国专利申请(专利号：ZL 99115885.5)，公开了名称为一种电阻抗断层成像方法，对EIT成像技术方案进行了详细披露。
在EIT中，合理选取正则化参数对于图像重建至关重要。若选取的正则化参数过小，会导致重建图像伪影过多，甚至不能得到正确的重建结果；若选取的正则化参数过大，会降低重建图像分辨率，使得重建目标产生严重的变形。目前，常用的正则化参数优选方法包括L曲线法和梯度下降法等。上述方法均是根据边界测量数据选取一个最优的正则化参数。然而，边界测量数据中包含了脑实质、颅骨层和头皮层区域的阻抗变化信息，这样选取的正则化参数为全局最优正则化参数。对于颅脑EIT，已知病灶如脑出血一般发生在脑实质区域，故我们更关注发生在脑实质的阻抗变化，头皮层与颅骨层的阻抗变化可以认为是噪声或者干扰引起的。基于上述分析，亟需寻找一定的方法将仅包含脑实质区域信息的数据从原来数据中剥离出来，然后基于处理后的数据再次优选正则化参数，来获得一个局部最优的正则化参数来提高颅脑EIT成像质量。
技术实现要素：
针对颅脑电阻抗图像重建合理选取正则化参数的问题，本发明目的在于提供一种颅脑EIT局部最优正则化参数选取的方法，通过该方法不仅可以提高信号的信噪比，同时可以获取更优的正则化参数，提高重建图像质量。
为了实现上述目的，本发明采取如下技术解决方案：
一种颅脑EIT局部最优正则化参数选取的方法，包括如下步骤：
步骤1，获取颅脑电阻抗初始重建数据f0，根据f0选取全局最优正则化参数λglobal并重建图像，获得初始重建图像，其中，颅脑电阻抗初始重建数据f0包含脑实质的阻抗变化信息、颅骨的阻抗变化信息与头皮层的阻抗变化信息；
步骤2，提取初始重建图像背景区域的阻抗变化ρBGD，初始重建图像背景区域为颅骨与头皮区域；
步骤3，计算背景区域的阻抗变化引起的边界电位变化fBGD；
步骤4，从颅脑电阻抗初始重建数据f0中去除背景区域的阻抗变化引起的边界电位变化fBGD，获得仅包含脑实质区域阻抗变化的数据fbrain；
步骤5，基于仅包含脑实质区域阻抗变化的数据fbrain选择局部正则化参数λlocal，以局部正则化参数λlocal作为颅脑EIT局部最优正则化参数。
所述步骤1中，根据f0基于L曲线法选取全局最优正则化参数λglobal。
所述步骤1中，基于阻尼最小二乘算法获取初始阻抗变化ρ0，ρ0＝(JTJ+λglobalWTW)-1JTf0，其中ρ0对应的图像为初始重建图像，J为敏感系数矩阵，W＝diag(J)为正则化矩阵。
所述步骤2中，初始重建图像背景区域的阻抗变化ρBGD＝Dρ0，其中，D为背景提取矩阵，背景提取矩阵D为一对角矩阵，其对角元素dii满足m表示单元编号，Ei表示第i个单元，Ω头皮颅骨表示头皮颅骨区域包含的单元集合，Ω脑实质表示脑实质区域包含的单元集合。
所述步骤3中，背景区域的阻抗变化引起的边界电位变化fBGD＝JρBGD，J为敏感系数矩阵。
所述步骤4中，仅包含脑实质区域阻抗变化的数据fbrain＝f0-fBGD。
所述步骤5中，基于仅包含脑实质区域阻抗变化的数据fbrain再次利用L曲线法选择局部正则化参数λlocal。
与现有技术相比，本发明具有如下有益效果：
本发明的颅脑EIT局部最优正则化参数选取的方法是基于颅脑脑损伤监护中，病灶通常发生在脑实质区域的先验信息，将颅骨层与头皮的阻抗变化从重建数据中的剔除，然后根据处理后的数据选择仅考虑了脑实质阻抗变化的局部最优正则化参数。本方法一方面可以利用空间分布信息抑制噪声与干扰，提高数据的信噪比，尤其是当头皮层存在强干扰的时，该方法可有效恢复正确的EIT图像。通过本方法选取的颅脑EIT局部最优正则化参数进行图像重建，能够避免正则化参数选取不合理对EIT图像重建的影响，改善重建的准确性，提高成像质量。
附图说明
图1是本发明颅脑电阻抗成像局部最优正则化参数选取方法的流程图；
图2是利用comsol建立的颅脑分层重构模型图；
图3是基于颅脑仿真模型图(仿真例一)；
图4是基于颅脑仿真模型图(仿真例二)；
图5是基于传统方法选取全局正则化参数对仿真例一的重建过程图；
图6是基于传统方法选取全局正则化参数对仿真例二的重建过程图；
图7是基于本发明方法选取局部正则化参数对仿真例一的重建过程图；
图8是基于本发明方法选取局部正则化参数对仿真例二的重建过程图。
具体实施方式
下面结合附图和实施例来对本发明作进一步的说明。
如图1所示，本发明颅脑电阻抗成像局部最优正则化参数选取方法采取如下技术方案：
1)获取颅脑电阻抗初始重建数据为f0，f0中包含了脑实质、颅骨与头皮层的阻抗变化信息。根据f0基于L曲线法选取全局最优正则化参数λglobal，再基于阻尼最小二乘算法获取初始阻抗变化ρ0：ρ0＝(JTJ+λglobalWTW)-1JTf0，其中ρ0对应的图像为初始重建图像，J为敏感系数矩阵，W＝diag(J)为正则化矩阵。
2)提取初始重建图像背景区域(颅骨与头皮区域)的阻抗变化ρBGD，ρBGD＝Dρ0。其中背景提取矩阵D为一对角矩阵，其对角元素dii满足m表示单元编号，Ei表示第i个单元，Ω头皮颅骨表示头皮颅骨区域包含的单元集合，Ω脑实质表示脑实质区域包含的单元集合。
3)计算出(2)中背景阻抗变化引起的边界电位变化fBGD，fBGD＝JρBGD；
4)从f0中去除背景阻抗变化引起的边界电位变化fBGD，获取仅包含脑实质区域阻抗变化的数据fbrain，fbrain＝f0-fBGD
5)基于仅包含脑实质区域阻抗变化的数据fbrain再次利用L曲线法选择局部正则化参数λlocal，以局部正则化参数λlocal作为颅脑EIT局部最优正则化参数。
参见图1，本实施例给出一种颅脑EIT局部最优正则化参数选取方法。其核心是提取初始重建图像中背景(头皮层与颅骨层)阻抗变化，然后将背景区域对应的边界电位变化从原始信号中去除，基于处理后的重建数据，再次选取局部最优正则化参数。该方法可以去除背景区域的噪声与干扰，提高颅脑EIT数据质量，同时为重建算法提供更为客观的最优正则化参数，减小重加误差，获取更为准确的重建图像。
在仿真施例中，按照如下步骤实施：
(1)如图2，为利用comsol软件建立二维颅脑分层仿真模型图，并对头皮层、颅骨层、脑实质层电导率进行设置：头皮层电导率设为0.44S/m，颅骨层电导率设为0.01259/m，脑实质电导率设为0.1499S/m；
(2)参照图3和图4，均为基于颅脑仿真模型图，模拟EIT监护脑出血，并设计两个仿真例(如图3和图4)来比较本文所提方法与传统方法所选正则化参数对EIT图像重建的影响。如图3所示，在脑实质区域设置一个圆形目标表示出血病灶，其电导率设置为0.7S/m，与血液相同；如图4所示，在头皮层设置导率为1.79S/m的长方形小目标，模拟出汗干扰，并于脑实质区域设置与仿真例一中相同的出血目标。利用comsol软件计算出血前后边界电位变化，并在电位变化中添加1％的高斯噪声，作为初始重建数据f0。
(3)利用传统方式选取正则化参数并完成仿真例一、仿真例二的图像重建。即根据仿真一、仿真例二中获取的f0，基于L曲线法选取一个全局最优正则化参数，然后根据阻尼最小二乘算法计算阻抗变化ρ0：ρ0＝(JTJ+λglobalWTW)-1JTf0，ρ0即为本发明方法中的初始重建图像。图5和图6为基于传统方法选取全局正则化参数，对仿真例一、仿真例二的重建结果。
(4)根据颅脑分层结构信息，构建背景区域(颅骨与头皮层区域)提取矩阵D，D为一对角矩阵，其每一个对角元素都对应着重建图像中的一个三角单元。设D的第i个对角元素为dii，m表示单元编号，Ei表示第i个单元，Ω头皮颅骨表示头皮颅骨区域包含的单元集合，Ω颅骨以内表示颅骨以内区域包含的单元集合。从初始图像ρ0中提取背景阻抗变化ρBGD，ρBGD＝Dρ0。
(5)计算出(4)中背景阻抗变化引起的边界电位变化fBGD，fBGD＝JρBGD，再从初始重建数据f0中去除背景阻抗变化对应的边界电位变化，获取仅包含脑实质区域阻抗变化的数据fbrain，fbrain＝f0-fBGD。
(6)基于(5)中fbrain再次利用L曲线法选择局部最优正则化参数λlocal，λlocal即为本发明方法所得的最优正则化参数。
(7)基于新方法选取的局部最优正则化参数λlocal，使用阻尼最小二乘算法完成图像重建：ρ＝(JTJ+λlocalWTW)-1JTfbrain。
如图7和图8所示，为基于本发明方法选取局部最优正则化参数，分别对仿真例一和仿真例二的重建结果。对于仿真例一，由图7结果与图5结果对比可知，在噪声存在的情况下，本发明方法对图像质量有较大的改善，其主要表现在两个方面：其一，新方法通过先验信息将背景区域的阻抗变化去除，故图7最终重建结果中该区域不存在伪影；第二，在去除背景层干扰后，信号质量有了提升，使得选取的局部最优正则化参数小于全局最优正则化参数，这样可以减小正则化对图像的模糊作用，由图像可知，相比图5，图7最终得到的重建结果形状更贴切原来的目标。对于仿真例二，在头皮层存在较强干扰的情况下，由图6可知，传统的方法仅能重建出头皮层的干扰，而如图8所示，本发明所提的方法可以重建出正确的图像。
需要说明的是，以上的实施例仅用于本领域的技术人员进一步理解本发明，本发明并不限于该实施例，凡是由本领域技术人员根据发明的技术方案做出的等效替换和增加，同样属于本发明保护的范围。