一阶相关免疫布尔函数的计数

田海建; 杨义先; 王建宇

一阶相关免疫布尔函数的计数

田海建^{①;杨义先},
杨义先,
王建宇

计量
- 文章访问数: 2030
- HTML全文浏览量: 90
- PDF下载量: 406
- 被引次数: 1
出版历程
- 收稿日期: 1995-11-06
- 修回日期: 1996-07-22
- 刊出日期: 1997-09-19

ENUMERATING CORRELATION-IMMUNE FUNCTIONS OF ORDER ONE

摘要

摘要: 本文研究线性结构布尔函数的相关免疫性,得到了一大类满足一阶相关免疫的线性结构函数,并大大改进了一阶相关免疫函数的计数结果,得到了目前为止的最好下界.
- 布尔函数; 序列密码; 相关免疫函数
Abstract: A large class of linear structure functions satisfying correlation immunity of order one is found by the studying of linear structural Boolean functions. The known enumeration bounds of correlation-immune Boolean functions are greatly improved. The best updated lower bounds are found.

HTML全文

1. 引言

脑出血 (IntraCerebral Hemorrhage, ICH)是由于脑内血管破裂或渗透引发的出血性病变，占脑卒中的15%～20%，其死亡率高达50%，是一种严重危及生命的疾病^[1]，对于脑出血疾病，早期、准确的检测及诊断对挽救患者生命至关重要，后面持续开展对脑出血病变的监测对于有效治疗也具有重大意义^[2]。

传统的医学影像检测技术有计算机断层成像 (Computed Tomography, CT)、磁共振成像 (Magnetic Resonance Imaging, MRI)和光学相干层析成像 (Optical Coherence Tomography, OCT)^[3,4]。CT重建图像的分辨率较高但病灶和背景间对比度低，且辐射剂量较大，MRI能对软组织成像且无放射性损害，但伪影也较CT多^[5]，OCT多用于眼底疾病的成像检测，也被证明其用于脑出血检测的可行性，但其穿透深度较低无法检测深层的颅脑疾病^[6]。此外，CT和MRI设备庞大，成本高、成像时间长，在紧急情况下应用受到限制，无法给患者提供快速的检测和诊断及实时和持续的监测^[3]。磁感应断层成像 (Magnetic Induction Tomography, MIT)是一种基于电磁检测原理对人体的被动电磁特性进行成像的功能性成像技术，具有非接触、安全无创、低成本、便携、无需耦合剂、可持续监测等特点^[7]。MIT因其产生的激励磁场能够穿透低电导率的头皮和颅骨，在医学成像领域有着巨大的潜在优势，能实现对脑出血的早期诊断和床旁实时、长期、无创的监测^[8]。因此，将MIT用于脑出血成像对脑出血的治疗、监测具有重大意义。

近年来，很多学者对脑出血MIT开展了研究。由于颅脑组织的复杂结构及其各向异性使正问题的求解变得十分困难，目前在脑出血MIT建模与仿真计算的研究中多使用简化的脑模型，如2维同心圆、3维同心球及圆柱模型等^[9-12]，与实际头部几何结构相差较大，在求解正问题时不可避免会带来相应的计算误差，导致后续图像重建结果误差增大；在MIT图像重建的研究中，柯丽等人^[13,14]建立包含头皮、颅骨、脑实质、脑出血的颅脑模型，并在其基础上进行了脑组织的阻抗分析，利用改进反投影算法进行图像重建，文献[8,15,16]对多频成像方法及线圈阵列方式进行了研究，并采用Tikhonov算法进行图像重建，Zhang等人^[17]提出了一种自适应阈值的分裂布雷格曼 (Split Bregman, SB)算法，降低了重建图像的噪声和面积误差，但并未分析重建图像的相关系数且最低只能对4 ml的脑出血进行成像，Chen等人^[18]提出了一种自适应阈值优化的Tikhonov算法，实现了不同出血量脑出血的成像，但重建图像的相关系数有待提高，周曦等人^[19]从初始值和正则化因子选取、灵敏度矩阵更新策略以及图像阈值优化等4个方面对牛顿-拉夫逊 (Newton Raphson, NR)算法进行改进，实现了不同噪声的图像重建，但重建图像质量还有待提高，Han等人^[12]在NR算法目标函数中加入加权矩阵和L1范数惩罚项提高了对脑出血图像重建的精度，但增大了计算量，算法收敛较慢。从上述研究可以看出，线性反投影 (Linear Back Projection, LBP)算法和Tikhonov正则化算法等单步图像重建算法可直接求解场域内的电磁特性分布，成像速度较快，但图像重建质量相对较低；Landweber算法、NR算法和共轭梯度最小二乘法 (Conjugate Gradients Least Squares, CGLS)等迭代图像重建算法需不断迭代求解场域内的电磁特性分布，成像速度比单步成像算法慢，但图像重建质量相对较高，在迭代图像重建算法中，NR算法应用最为普遍，其成像质量相对较高，收敛性相对其他迭代图像重建算法较好，但在脑出血MIT中由于组织的电导率极小，发生脑出血时测量数据较小的变化会对重建结果产生较大的影响，导致算法的稳定性和收敛性较差，目标与背景间伪影较大，图像重建质量也有待提高，因此有必要对其进行改进，提高其稳定性、收敛性和重建质量，同时高质量的重建图像也为后续的医学图像处理及临床中进行辅助诊断和评估患者预后提供基础^[20]。

针对上述问题，为提高解的稳定性和重建图像的质量，本文提出一种改进的NR算法，在目标函数中加入自适应加速正则化惩罚项和L2范数正则化惩罚项，并引入投影算子P施加物理意义上的约束，提高算法每一步迭代的效率，减少重建图像的伪影，改善成像质量；其次为验证算法在脑出血MIT中的实用性以及提高脑模型的真实性，构建了包含头皮、颅骨、脑脊液和脑实质4种组织的真实3维颅脑模型，并开展了3个位置处出血量各为24 ml, 14 ml, 2 ml的脑出血仿真实验，通过模型计算得到图像重建所需的相位差检测值和灵敏度矩阵，然后利用改进NR算法与Tikhonov, Landweber, NR, CGLS、优化迭代策略的NR等算法分别对设置的脑出血进行磁感应断层成像，验证所提算法的有效性。

2. 脑出血MIT检测原理

脑出血MIT检测的基本原理：采用电磁感应原理，向激励线圈中通入正弦交流电，激励线圈在被测区域内产生初级交变磁场 ${{{\boldsymbol{B}}}_0}$ ，初级交变磁场 ${{{\boldsymbol{B}}}_0}$ 作用于人脑时，会在脑组织中感应出涡流，从而产生次级磁场 $\Delta {{\boldsymbol{B}}}$ ，在不同病理状态下次级磁场 $\Delta {{\boldsymbol{B}}}$ 的值会发生改变，可通过检测线圈得到初级磁场 ${{{\boldsymbol{B}}}_0}$ 和叠加磁场 ${{{\boldsymbol{B}}}_0} + \Delta {{\boldsymbol{B}}}$ 之间的相位差 $\Delta {{\boldsymbol{\varphi }}}$ ，再通过图像重建算法得到成像域内的电导率分布图像，进而通过电导率分布图像诊断脑出血的严重程度。MIT的基本原理及初、次磁场关系如图1所示。

图 1 MIT的基本原理及初、次磁场间关系

下载: 全尺寸图片幻灯片

对于较小的电导率扰动 $\Delta {{\boldsymbol{\sigma}} }$ ， $\Delta {{\boldsymbol{\varphi}} }$ 和 $\Delta {{\boldsymbol{\sigma}} }$ 间的关系可以简化为近似离散线性方程^[21]

$\Delta {{{\boldsymbol{\varphi}} }_{M \times 1}}{ = }{{{\boldsymbol{S}}}_{M \times N}} \cdot \Delta {{{\boldsymbol{\sigma}} }_{N \times 1}} + {{\boldsymbol{v}}}$

(1)

其中， $\Delta {{{\boldsymbol{\varphi }}}_{M \times 1}} \in {\boldsymbol{R}}$ 是线圈测量相位差， ${{{\boldsymbol{S}}}_{M \times N}} \in {\boldsymbol{R}}$ 是灵敏度矩阵或雅可比矩阵，表示成像区域内电导率分布变化引起的线圈检测相位差变化率^[21]， $\Delta {{{\boldsymbol{\sigma}} }_{N \times 1}} \in {\boldsymbol{R}}$ 是被测区域内的电导率分布变化， $M$ 是激励-检测线圈组合数量， $N$ 是被测域内的网格数量， ${\boldsymbol{{v}}}$ 是噪声。

脑出血MIT是将通过模型或实验获得的检测值 $\Delta {{\boldsymbol{\varphi}} }$ 和成像域内的灵敏度矩阵 ${{\boldsymbol{S}}}$ 求解未知的电导率特性分布 $\Delta {{\boldsymbol{\sigma}} }$ ，但是由于MIT中线圈数量有限，因此线圈检测值 $\Delta {{\boldsymbol{\varphi}} }$ 的数量远远小于待求解的电导率特性分布 $\Delta {{\boldsymbol{\sigma }}}$ 的数量，同时 $\Delta {{\boldsymbol{\varphi }}}$ 的微小变化将对 $\Delta {{\boldsymbol{\sigma}} }$ 的值产生较大影响，因此式(1)是一个欠定、病态的方程，需要通过一定的图像重建算法进行求解，将求解得到的 $\Delta {{\boldsymbol{\sigma}} }$ 按幅值大小可绘制成反映脑出血情况的图像。

3. 改进的Newton-Raphson图像重建算法

脑出血MIT图像重建是根据检测信号和场域内的灵敏度分布重建头部不同组织中的电导率特性分布。NR算法的基本思想是不断迭代更新电导率分布值，寻找电导率分布的最优解，使得目标函数的值最小^[22]。采用NR算法求解MIT图像重建问题时通常以相位差测量值和计算值的误差范数平方作目标函数

$\Delta \mathop {{\boldsymbol{\sigma}} }\limits^ \sim = \arg \min \left( {\frac{1}{2}\left\| {{\boldsymbol{F}}(\Delta {{\boldsymbol{\sigma}} }) - \Delta {{\boldsymbol{\varphi}} }} \right\|_2^2} \right)$

(2)

其中， ${\boldsymbol{F}}(\Delta {{\boldsymbol{\sigma}} })$ 和 $\Delta {{\boldsymbol{\varphi}} }$ 分别是相位差的计算值和测量值，NR算法的迭代公式为

$\Delta {{\boldsymbol{\sigma}}}_{{k+}1}=\Delta {{\boldsymbol{\sigma}}}_{{k}}-{\boldsymbol{H}}_{k}^{-1}\nabla {f}(\Delta {{\boldsymbol{\sigma}}}_{{k}})$

(3)

其中， $\Delta {{\boldsymbol{{\sigma}} }_{{k + }1}}$ 是待求解的成像区域内电导率分布变化， $\nabla {f}(\Delta {{\boldsymbol{\sigma}}}_{{k}})$ 是梯度方向， ${{{\boldsymbol{H}}}_k}$ 是Hessian矩阵。

NR算法每一步迭代都需要求解目标函数的2阶偏导数，计算复杂度较大，同时Hesse矩阵具有较严重的病态性，尤其在脑出血MIT时生物组织的电导率极小，脑出血测量数据较小的变化会对重建结果产生较大的影响，导致算法的稳定性和收敛性较差，重建图像中背景与脑出血间的伪影较大。

采用优化方法求解式(2)的过程是不稳定的，通过在目标函数中增加L2范数正则化惩罚项 $\alpha \left\| {{{\boldsymbol{R}}}\Delta {{{\boldsymbol{\sigma}} }_k}} \right\|_2^2$ 来降低解的不稳定性，同时为了提高NR算法中每步迭代的计算效率，提高图像重建质量，在目标函数中加入与前一刻计算结果相关的自适应加速正则化惩罚项 ${\beta _k}\left\| {\Delta {{{\boldsymbol{\sigma}} }_k} - \Delta {{{\boldsymbol{\sigma}} }_{k - 1}}} \right\|_2^2$ 进行约束，加快算法的收敛速度，该惩罚项在目标函数中的权重随着迭代次数的增加而自适应改变，改进算法的目标函数为

$\begin{split} \Delta \mathop {{\boldsymbol{\sigma}} }\limits^ \sim =& {\text{arg}}\min \frac{1}{2}\left\{ {\left\| {{\boldsymbol{F}}(\Delta {{\boldsymbol{\sigma}} }) - \Delta {{\boldsymbol{\varphi }}}} \right\|_2^2} \right. \\ & + \left. {{\beta _k}\left\| {\Delta {{{\boldsymbol{\sigma}} }_k} - \Delta {{{\boldsymbol{\sigma}} }_{k - 1}}} \right\|_2^2 + \alpha \left\| {{{\boldsymbol{R}}}\Delta {{{\boldsymbol{\sigma}} }_k}} \right\|_2^2} \right\} \end{split}$

(4)

其中， ${\boldsymbol{F}}(\Delta {{\boldsymbol{\sigma}} })$ 是相位差计算值， $\Delta {{\boldsymbol{\varphi}} }$ 是相位差检测值， ${\beta _k}$ 是自适应加权系数，随着迭代次数的增加而自适应改变该惩罚项的权重， $\Delta {{{\boldsymbol{\sigma}} }_k}$ 是当前时刻电导率分布变化， $\Delta {{{\boldsymbol{\sigma}} }_{k - 1}}$ 是前一时刻电导率分布变化， $\alpha > 0$ ，是用于控制解平滑性和连续性的正则化参数，从理论上确定其值较为困难，一般根据经验进行设定， ${{\boldsymbol{R}}}$ 是正则化矩阵，一般为单位矩阵。将式(4)展开并求其下降梯度 $\nabla {f}(\Delta {\boldsymbol{\sigma}}_{{k}})$ 和Hesse矩阵 ${{{\boldsymbol{H}}}_k}$

$\begin{split} \nabla {f}\left( {\Delta {{\boldsymbol{\sigma}}_{k}}} \right) =& {\left( {{\boldsymbol{F}}'\left( {\Delta {\boldsymbol{\sigma}}} \right)} \right)^{\text{T}}}\left( {{\boldsymbol{F}}(\Delta {\boldsymbol{\sigma}}) - \Delta {{\boldsymbol{\varphi}} }} \right) \\ & + {\beta _k}\left( {\Delta {{\boldsymbol{\sigma}}_k} - \Delta {{\boldsymbol{\sigma}}_{k - 1}}} \right) + \alpha {{{\boldsymbol{R}}}^{\text{T}}}{{\boldsymbol{R}}}\Delta {{\boldsymbol{\sigma}}_k} \end{split}$

(5)

${{{\boldsymbol{H}}}_k} = {\left( {{\boldsymbol{F}}'\left( {\Delta {{\boldsymbol{\sigma}} }} \right)} \right)^{\text{T}}}{\boldsymbol{F}}'\left( {\Delta {{\boldsymbol{\sigma}} }} \right) + \alpha {{{\boldsymbol{R}}}^{\text{T}}}{{\boldsymbol{R}}}$

(6)

其中， ${\boldsymbol{F}}'\left( {\Delta {{\boldsymbol{\sigma}} }} \right)$ 是灵敏度矩阵 ${{\boldsymbol{S}}}$ ，将 ${\boldsymbol{F}}'\left( {\Delta {{\boldsymbol{\sigma}} }} \right)$ 用 ${{\boldsymbol{S}}}$ 表示，并将式(5)、式(6)代入式(3)后得到

$\begin{split} \Delta {{\boldsymbol{\sigma} }_{{k + }1}} =& \Delta {{\boldsymbol{\sigma} }_{k}} - {\left( {{{{\boldsymbol{S}}}^{\text{T}}}{{\boldsymbol{S}} + }\alpha {{{\boldsymbol{R}}}^{\text{T}}}{{\boldsymbol{R}}}} \right)^{ - 1}}\\ & \cdot\left[ {{{\boldsymbol{S}}}^{\text{T}}}\left( {{{\boldsymbol{S}}}\Delta {{\boldsymbol{\sigma} }_{k}} - \Delta {\boldsymbol{{\varphi}} }} \right) + {\beta _k}\left( {\Delta {{\boldsymbol{\sigma} }_k} - \Delta {{\boldsymbol{\sigma} }_{k - 1}}} \right)\right. \\ & \left.+ \alpha {{{\boldsymbol{R}}}^{\text{T}}}{{\boldsymbol{R}}}\Delta {{\boldsymbol{\sigma} }_k} \right]\\[-10pt] \end{split}$

(7)

其中， ${\beta _k}$ 除了可以根据经验选取为固定值，还可以利用v方法求解其在每次迭代中的最优值，v方法是由Brakhage提出的一种用于求解不适定问题或条件数很差的适定问题的计算方法，与其他用于求解此类问题的方法相比，采用v方法计算目标函数误差项的加权系数，仅需较少的计算步骤，便可得到较小的最小误差范数^[23]。计算公式为

$\left. \begin{aligned} &{{\beta _0} = 0} \\ & {{\beta _k} = \frac{{\left( {2k + 2v + 1} \right)\left( {2k - 1} \right)k}}{{\left( {2k + 4v + 1} \right)\left( {2k + 2v - 1} \right)\left( {k + 2v} \right)}},}\;\;{k \ge 1} \end{aligned} \right\}$

(8)

为了进一步加快算法的收敛，将LBP算法计算得到的电导率分布值作为改进算法的迭代初值，同时根据脑出血MIT检测的实际情况在改进算法中引入投影算子 ${P}$ 确保解收敛。综上，所提改进NR算法的迭代公式如式(9)，其流程图如图2所示。

图 2 改进NR算法流程图

下载: 全尺寸图片幻灯片

$\left.\begin{aligned} & \Delta {{\boldsymbol{\sigma}} }_{-1}=0 \\ & \Delta {{\boldsymbol{\sigma}} }_{0}=\Delta {{\boldsymbol{\sigma}} }_{\text{LBP}} \\ & \Delta {{\boldsymbol{\sigma}} }_{{k+}1}={P}\left[\Delta {{\boldsymbol{\sigma}} }_{{k}}-{\left({\boldsymbol{S}}^{\text{T}}{{\boldsymbol{S}}+}\alpha {\boldsymbol{R}}^{\text{T}}\boldsymbol{R}\right)}^{-1}\left[{\boldsymbol{S}}^{\text{T}}\left(\boldsymbol{S}\Delta {{\boldsymbol{\sigma}} }_{{k}}-\Delta {\boldsymbol{\varphi}}\right)+{\beta }_{k}\left(\Delta {{\boldsymbol{\sigma}} }_{k}-\Delta {{\boldsymbol{\sigma}} }_{k-1}\right)+\alpha {\boldsymbol{R}}^{\text{T}}{\boldsymbol{R}}\Delta {{\boldsymbol{\sigma}} }_{k}\right]\right] \\ & P\left[{f}(\Delta {\boldsymbol{\sigma}} )\right]=\left\{\begin{array}{llllll} 0,& \text{ }{f}(\Delta {\boldsymbol{\sigma}} )\le 0\\ {f}(\Delta {\boldsymbol{\sigma}} ),& \text{ }{f}(\Delta {\boldsymbol{\sigma}} ) > 0\end{array}\right. \end{aligned}\right\}$

(9)

其中，LBP算法的计算结果作为改进NR算法的迭代初值； $P$ 是非负凸集上的投影算子，通过投影算子 $P$ 引入物理意义上的约束，确保每次迭代的解都非负有界且收敛于一个凸集上； ${\beta _k}$ 使用v方法确定其最优值，v的值一般通过经验法选取^[23]，本文v取8。在迭代过程中设置迭代终止条件为 $\left| \Delta {{{\boldsymbol{\sigma}} }_{k + 1}} - \Delta {{{\boldsymbol{\sigma}} }_k} \right| < \varepsilon$ 或 $k > {k_{\max }}$ ，当迭代过程中两个相邻解间的差小于设定的较小阈值 $\varepsilon = 1 \times {10^{ - 12}}$ 时停止迭代，同时设置了最大迭代次数 ${k_{\max }} = 1\;000$ ，防止算法陷入半收敛或无限迭代，进行重复的无效计算，占用计算资源。满足迭代终止条件时输出电导率分布 $\Delta {{{\boldsymbol{\sigma}} }_{{k + }1}}$ ，并按其幅值大小绘制得到重建图像。

4. 建模与仿真

4.1 脑出血MIT有限元模型的建立

构建具有真实几何形状的颅脑模型，能提高对模型涡流场、磁场的求解精度，在此基础上开展线圈检测值和场域内灵敏度的仿真计算，为后续的图像重建提供可靠的仿真实验结果。建立的仿真模型共包含3个域：真实颅脑模型、空气域、线圈域。颅脑模型的几何结构采用电气与电子工程师协会 (Institute of Electrical and Electronics Engineers, IEEE)、国际电工委员会 (International Electrotechnical Commission, IEC)和欧洲电工标准化委员会 (European Committee for Electrotechnical Standardization, CENELEC)提供的标准比吸收率 (Specific Absorption Rate, SAR)测量规范中相同的几何结构，利用COMSOL Multiphysics对MRI数据进行放样、转换、旋转、缩放等操作后得到3维头部轮廓；将具有类似介电特性的组织进行合并，保留头皮、颅骨、脑脊液和脑实质等4种主要头部组织^[24]，并将得到的3维头部轮廓按各组织1:0.92:0.87:0.84的比例进行缩放^[25]，最终得到具有4种组织的真实3维颅脑模型，空间中笛卡尔坐标系的原点位于建立的脑模型的几何中心。

为模拟实际检测环境，设置一个足够大的球型空气域，半径为70 cm；根据线圈检测特性和颅脑模型的结构，采用8线圈均匀布置的方式，依次编号1～8，每个线圈大小、形状完全一致，均匀分布于脑模型中与头皮距离为12 cm的圆上，线圈外径为20 mm，内径为15 mm，高为30 mm，匝数为100，线圈导线电导率为 $6 \times {10^7}{\text{ S/m}}$ ，导线的截线面积定义为 ${10^{ - 7}}{\text{ }}{{\text{m}}^{\text{2}}}$ ，建立的脑出血MIT有限元模型如图3所示。

图 3 脑出血MIT有限元模型

下载: 全尺寸图片幻灯片

将模型中线圈的材料属性设置为铜，球形空气域的材料属性为空气；在脑出血MIT中常用的激励频率在1～10 MHz^[16]，将激励频率设置为1 MHz，线圈中激励电流为1 A，并按表1设置模型中各头部组织的电磁特性参数(数据来源文献[24])。

表 1 1 MHz下的脑组织电磁特性

	脑组织
	头皮	颅骨	脑脊液	脑实质	脑出血
电导率(S/m)	0.044	0.024	2.000	0.102	0.822
介电常数(F/m)	50.8	145	109	480	3030

下载: 导出CSV

| 显示表格

在脑模型中分别设置了3个不同位置的脑出血，分别位于空间位置A(0, 0.03, 0.075)处，位置B(–0.026, –0.04, 0.075)处及位置C(0.02, 0.02, 0.075)处，在每个位置分别设置出血量为24 ml, 14 ml及2 ml的脑出血。设置的脑出血分布情况如图4所示。

图 4 脑出血分布情况

下载: 全尺寸图片幻灯片

4.2 仿真计算

4.2.1 相位差仿真计算

采用循环激励，循环检测的模式，即其中1个线圈为激励，其余7个进行检测，当遍历8个线圈进行激励后共得到56个检测相位值。仿真计算图4中位置A, B, C处3个出血量的相位差如图5所示。

图 5 位置A,B,C处3个出血量的相位差

下载: 全尺寸图片幻灯片

图5(a)—图5(c)分别表示位置A, B, C处3个出血量的相位差，每个图中黑色、蓝色、红色虚线分别表示在该位置24 ml, 14 ml, 2 ml出血量情况下的仿真计算得到相位差测量值，由图5(a)可以看出，检测相位差数量级很小；脑出血体积的大小对线圈测量值影响很大，出血量越大检测相位差的幅值变化越大，出血量越小检测相位差的幅值变化越小，位置B, C处不同出血量的测量相位差规律与此类似。

4.2.2 灵敏度仿真计算

灵敏度矩阵是图像重建的先验信息^[15]，常见的灵敏度矩阵求解方法有扰动法和互易定理法，采用基于互易定理的计算方法能大量减少计算时间^[16]，灵敏度矩阵通过图3中建立的脑出血MIT有限元模型仿真计算得到，假设有激励-检测线圈对a和b，分别对a和b施加激励，另外一个线圈作检测，此时可计算场域内各网格上电导率发生变化时的灵敏度，计算公式为

${S_{({a}{{-}}1){N_{S}} + {b}{\text{,}}{i}}} = {{\rm{Im}}} \left( {\int\limits_{{{\varOmega }_i}} {\frac{{{\boldsymbol{E}}_i^a \cdot {\boldsymbol{E}}_i^b}}{{{{\left( {{U_0}} \right)}_{{a}{\text{,}}{b}}}}}} {\text{d}}{\varOmega }} \right)$

(10)

其中， ${N_S}$ 是检测线圈数量； ${{\rm{Im}}} ( \cdot )$ 代表虚部； ${{\varOmega }_i}$ 是成像区域内网格单元； ${\boldsymbol{E}}_i^a$ 和 ${\boldsymbol{E}}_i^b$ 分别是线圈a和线圈b施加激励时检测区域内的电场分布； ${\left( {{U_0}} \right)_{{a}{\text{,}}{b}}}$ 是空场时由激励-检测线圈ab计算得到的感应电压。建立的模型是8线圈传感器，因此检测线圈数量 ${N_S} = 7$ 。在参考频率为1 MHz的情况下，通过将整个检测区域设置为1 S/m的均匀电导率来计算灵敏度矩阵 ${{\boldsymbol{S}}}$ 。通过有限元仿真计算相关场量并根据式(10)计算灵敏度矩阵，计算6种典型线圈组合下的灵敏度如图6所示。

图 6 不同激励-检测组合的灵敏度图

下载: 全尺寸图片幻灯片

图6(a)—图6(f)分别表示线圈1-线圈5、线圈4-线圈8、线圈2-线圈3、线圈5-线圈6、线圈5-线圈8、线圈3-线圈6组合时的灵敏度图，每两组分别代表相对、相邻、相隔线圈激励检测情况下的典型灵敏度分布情况，其中不同的颜色反映该区域对电导率变化的敏感度不同，颜色对应数值的绝对值越大表示该区域灵敏度越高，颜色对应数值的绝对值越小表示该区域灵敏度越低；从图6可以看出，在单线圈激励单线圈检测的情况下，越靠近激励或检测线圈的区域，灵敏度的绝对值越大，远离该区域灵敏度绝对值逐渐减小，如图6(a)中，线圈1-线圈5组合分别作激励或检测线圈时，越靠近线圈1和线圈5区域的灵敏度绝对值越大，该区域发生脑出血引起颅脑内电导率分布发生变化时，检测线圈测量到的信号变化也越大，因此脑出血在靠近线圈的区域时成像效果更好，精度也更高。

5. 实验及结果分析

为了验证所提改进NR算法的有效性，建立了真实颅脑和线圈的有限元模型，并设置了图3中9种情况的脑出血，通过模型仿真计算得到9种脑出血情况下的线圈检测相位差和灵敏度矩阵后，应用所改进的NR算法进行图像重建实验，并与NR算法^[26]、Tikhonov算法^[8]、CGLS算法^[27]、Landweber算法^[28]和优化迭代策略的NR算法^[19]进行比较，采用改进的NR算法和5种对比算法分别对z=0.075平面进行断层成像，图像重建结果如图7所示。实验在处理器Intel i7-10875H，内存16 GB，显卡Nvidia GeForce RTX 2060, 6 GB显存的电脑上进行。

下载: 全尺寸图片幻灯片

在图7(a)—图7(c)、图7(d)—图7(f)、图7(g)—图7(i)分别对应位置A处、B处及C处出血量分别为24 ml, 14 ml, 2 ml共9种情况的脑出血图像重建结果，其中的黑色虚线表示脑出血的实际轮廓。从图7可以看出，6种成像算法均能重建出不同出血情况的脑出血位置，在同一位置下，随着出血量的减小，与之对应区域的重建电导率幅值也减小，这是因为出血量减小，检测的相位差信号也相对减小；Tikhonov, Landweber, CGLS, NR、优化迭代策略的NR等5种算法的重建图像，出血的范围都大于实际范围，而使用所提改进NR算法的重建图像，出血范围与实际范围更吻合，重建图像中脑出血周围的伪影比Tikhonov, Landweber, CGLS, NR、优化迭代策略的NR等5种算法少，且脑出血与背景的边界更清晰，对出血量的判断也更准确，可有效实现对2 ml脑出血的图像重建。

图 7 图像重建结果

下载: 全尺寸图片幻灯片

为了评价不同算法的图像重建性能，采用图像相关系数(Correlation Coefficient, CC)、图像重建误差(Image Error, IE)、归一化均方距离(Normalization Mean Square Distance criterion, NMSD)、图像重建时间和迭代次数作为评价指标。其中图像重建时间是算法开始到重建图像所花费的时间，迭代次数是目标函数值小于设定阈值 $\varepsilon$ 时已进行迭代计算的次数，各评价指标的计算公式为

$\qquad\qquad\qquad\qquad {\text{CC}} = \frac{{\displaystyle\sum\limits_{i = 1}^N {\left( {\Delta {{\boldsymbol{\sigma} }_{{\text{true}}}}\left( i \right) - \Delta {{\overline {\boldsymbol{\sigma} } }_{{\text{true}}}}} \right)} \left( {\Delta {{\boldsymbol{\sigma} }_{{\text{rec}}}}\left( i \right) - \Delta {{\overline {\boldsymbol{\sigma} } }_{{\text{rec}}}}} \right)}}{{\sqrt {{{\displaystyle\sum\limits_{i = 1}^N {\left( {\Delta {{\boldsymbol{\sigma} }_{{\text{true}}}}\left( i \right) - \Delta {{\overline {\boldsymbol{\sigma} } }_{{\text{true}}}}} \right)^2} }}\displaystyle\sum\limits_{i = 1}^N {{{\left( {\Delta {{\boldsymbol{\sigma} }_{{\text{rec}}}}\left( i \right) - \Delta {{\overline {\boldsymbol{\sigma} } }_{{\text{rec}}}}} \right)}^2}} } }}$

(11)

$\qquad\qquad\qquad\qquad {\text{IE}} = \frac{{\left\| {\Delta {{\boldsymbol{\sigma} }_{{\text{true}}}} - \Delta {{\boldsymbol{\sigma} }_{{\text{rec}}}}} \right\|}}{{\left\| {\Delta {{\boldsymbol{\sigma} }_{{\text{true}}}}} \right\|}}$

(12)

${\text{NMSD = }}\sqrt {\frac{{\displaystyle\sum\limits_{i = 1}^N {{{\left( {\Delta {{\boldsymbol{\sigma} }_{{\text{true}}}}\left( i \right) - \Delta {{\boldsymbol{\sigma} }_{{\text{rec}}}}\left( i \right)} \right)}^2}} }}{{\displaystyle\sum\limits_{i = 1}^N {{{\left( {\Delta {{\boldsymbol{\sigma} }_{{\text{true}}}}\left( i \right) - \Delta {{\overline {\boldsymbol{\sigma} } }_{{\text{true}}}}} \right)}^2}} }}}$

(13)

其中， $\Delta {{\boldsymbol{\sigma} }_{{\text{true}}}}$ 是真实电导率分布， $\Delta {{\boldsymbol{\sigma} }_{{\text{rec}}}}$ 是重建的电导率分布， $\Delta {\overline {\boldsymbol{\sigma} } _{{\text{true}}}}$ 和 $\Delta {\overline {\boldsymbol{\sigma} } _{{\text{rec}}}}$ 分别是 $\Delta {{\boldsymbol{\sigma} }_{{\text{true}}}}$ 和 $\Delta {{\boldsymbol{\sigma} }_{{\text{rec}}}}$ 的平均值。利用CC和NMSD评价重建图像的分辨率^[11]，CC越大，NMSD越小，表示重建图像分辨率更优；CC越大，IE和NMSD越小，表示图像重建质量越高，分别计算9种出血情况下6种算法重建图像的相关系数、重建误差、归一化均方距离、图像重建时间和迭代次数，并分析了所提改进NR算法中投影算子 ${P}$ 对收敛速度提高和成像质量改善的影响，得到的结果如表2、表3和表4所示。

表 2 相关系数

脑出血分布	Tikhonov	Landweber	CGLS	NR	NR (优化迭代策略)	改进NR (无投影算子)	改进NR
位置A：24 ml	0.505	0.355	0.490	0.483	0.512	0.513	0.676
位置A：14 ml	0.437	0.306	0.422	0.416	0.443	0.441	0.595
位置A： 2 ml	0.269	0.209	0.255	0.253	0.275	0.267	0.377
位置B：24 ml	0.455	0.311	0.449	0.416	0.468	0.448	0.590
位置B：14 ml	0.378	0.282	0.380	0.342	0.388	0.367	0.401
位置B： 2 ml	0.207	0.171	0.199	0.189	0.213	0.202	0.286
位置C：24 ml	0.481	0.364	0.476	0.446	0.491	0.474	0.619
位置C：14 ml	0.407	0.308	0.404	0.375	0.417	0.397	0.528
位置C： 2 ml	0.244	0.180	0.243	0.227	0.250	0.240	0.333

下载: 导出CSV

| 显示表格

表2、、分别是重建图像的相关系数、图像重建误差和归一化均方距离、图像重建时间和迭代次数，其中表的第1列是9种脑出血分布情况，第2列—第8列是9种出血情况下6种算法重建图像以及无投影算子的改进NR算法的相关系数、重建误差和归一化均方距离、图像重建时间和迭代次数；对比3个表中有无投影算子的改进NR算法，可以看出施加投影算子 ${P}$ 后，算法的相关系数提高，重建误差和归一化均方距离均减小，且图像重建时间和迭代次数均减少，平均为NR算法的1/2，表明通过投影算子 ${P}$ 施加物理意义上的约束后提高了重建图像的成像质量和收敛速度；由表2和表3可以看出，对于9种脑出血情况，采用所提改进NR算法重建图像的相关系数均大于其他5种算法重建图像的相关系数，归一化均方距离均小于其他算法的归一化均方距离，重建图像的分辨率更优，同时其重建误差也均小于其他算法的图像重建误差，表明所提改进NR算法重建图像的质量更高；在同一出血位置，随着出血量的减小，各算法重建图像的相关系数减小，重建误差增大，这是因为当出血量较小时，检测到的相位信号变化也较小，图像重建难度增大，同时对于2 ml微小出血的重建误差较大，这是由于出血量较小检测信号微弱以及检测信号由各个组织及其相互间耦合的信号构成，可通过研究更有效的图像重建算法、在图像重建过程中引入脑组织的结构先验信息以及脑出血MIT信号的有效分解来进一步降低重建误差。结合图7所给出的图像重建结果，同一位置处出血量减少时，所提改进NR算法得到的脑出血范围更符合实际脑出血范围，其余5种算法得到的脑出血范围虽然有所减小，但依然大于实际范围，同时所提改进NR算法的相关系数更高、重建误差和归一化均方距离更低，并且能实现2 ml脑出血的图像重建，表明所提改进NR算法求解得到的颅脑截面电导率分布及其数值解更符合实际，求解过程更加稳定。

表 3 图像误差/归一化均方距离

脑出血分布	Tikhonov	Landweber	CGLS	NR	NR (优化迭代策略)	改进NR (无投影算子)	改进NR
位置A：24 ml	0.962/1.182	0.993/1.005	0.956/0.968	0.957/0.969	0.961/0.975	0.959/0.970	0.910/0.949
位置A：14 ml	1.022/1.030	1.510/1.527	1.163/1.174	1.184/0.987	0.991/0.977	0.975/0.983	0.947/0.876
位置A： 2 ml	1.867/1.872	2.709/2.610	2.505/2.514	2.517/2.525	1.662/1.666	1.860/1.865	1.544/1.550
位置B：24 ml	0.970/0.983	0.994/1.007	0.967/0.980	0.969/0.983	0.970/0.984	0.971/0.984	0.938/0.967
位置B：14 ml	1.013/1.022	1.525/1.540	1.043/1.052	1.135/1.145	0.997/0.993	1.144/1.153	0.968/0.920
位置B： 2 ml	1.825/1.779	2.226/2.247	2.098/2.111	2.256/2.262	1.716/1.675	1.887/1.891	1.501/1.505
位置C：24 ml	0.965/0.978	0.993/1.006	0.961/0.974	0.964/0.977	0.966/0.980	0.965/0.978	0.933/0.961
位置C：14 ml	0.990/0.989	1.540/1.567	1.038/1.046	1.101/1.111	0.981/0.967	1.107/1.117	0.962/0.910
位置C： 2 ml	1.663/1.622	1.877/1.765	1.855/1.860	2.026/2.038	1.575/1.510	1.705/1.710	1.354/1.358

下载: 导出CSV

| 显示表格

表 4 图像重建时间(s)/迭代次数

脑出血分布	Tikhonov	Landweber	CGLS	NR	NR (优化迭代策略)	改进NR (无投影算子)	改进NR
位置A：24 ml	0.064/1	0.027/300	0.002/40	0.750/12	1.949/33	0.411/8	0.223/5
位置A：14 ml	0.078/1	0.029/303	0.002/45	0.801/13	1.841/30	0.406/8	0.217/5
位置A： 2 ml	0.071/1	0.041/460	0.002/52	1.134/17	2.133/35	0.410/8	0.287/6
位置B：24 ml	0.061/1	0.026/300	0.001/38	0.706/11	1.534/25	0.405/8	0.285/6
位置B：14 ml	0.068/1	0.032/332	0.002/40	0775/12	1.797/30	0.391/8	0.235/5
位置B： 2 ml	0.063/1	0.042/510	0.002/46	0.983/15	1.910/32	0.393/8	0.293/6
位置C：24 ml	0.061/1	0.026/302	0.002/39	0.765/12	1.527/26	0.414/8	0.253/5
位置C：14 ml	0.070/1	0.030/303	0.002/40	0.770/12	1.678/28	0.413/8	0.297/6
位置C： 2 ml	0.064/1	0.030/310	0.002/44	1.053/15	1.873/31	0.396/8	0.295/6

下载: 导出CSV

| 显示表格

由表4可以看出，对于9种脑出血情况，改进NR 算法重建图像的时间比Tikhonov, Landweber, CGLS等3种算法较多，这是因为Tikhonov仅需一步成像，Landweber和CGLS在迭代过程中只需计算目标函数的1阶导数^[29,30]，而NR算法每步迭代都需要计算其2阶导数，导致图像重建时间较多。但所提改进NR算法相比NR和优化迭代策略的NR算法所需的重建时间大幅减少，成像时间平均只需NR算法的1/3。结合表2和表3分析，与Landweber, CGLS, NR和优化迭代策略的NR算法4种迭代图像重建算法相比，所提改进NR算法迭代次数最少，重建图像质量最高，表明所提改进NR算法每一步迭代的效率更高，收敛性更好且能够得到更高质量的重建图像。

6. 结束语

本文提出一种用于脑出血MIT检测的改进NR图像重建算法，首先建立了一个包含头皮、颅骨、脑脊液和脑实质4种组织的真实3维颅脑仿真模型，并仿真计算了9种脑出血情况下的线圈检测值和灵敏度矩阵，采用改进的NR算法和其他5种图像重建算法开展了图像重建实验。实验结果表明，对于不同脑出血分布情况，所提改进NR算法的重建图像相比于其他5种图像重建算法在相关系数上提高了0.013～0.321，在重建误差上减小了0.019～1.165，在归一化均方距离上减小了0.013～1.060，成像时间平均只需NR算法的1/3，改进NR算法能实现较少的迭代次数和时间重建出更高质量的图像，且能实现2 ml微小脑出血的图像重建，表明所提改进NR算法的迭代效率更高、稳定性更好、重建图像质量更高，由此可见，所提改进NR算法用于脑出血MIT是可行且有效的，有力推动了MIT技术的发展。

在本研究中，考虑到仿真模拟所需的计算资源，简化了脑组织仅保留了主要的头部组织，而真实脑组织更为复杂，同时改进NR算法重建图像中脑出血和背景之间仍存在一定的伪影，下一步研究将构建更为精细的脑组织模型，并搭建脑出血MIT实验平台，利用琼脂块或3D打印技术制作不同脑出血情况下的真实3维颅脑模型，采集相关实验数据，进一步验证所提改进NR算法有效性，并开展应用研究。

参考文献(1)

Mitchell C. Enumerating Boolean function of cryptographic significance. J. Cryptology,1990,2(3): 155-170.[2]杨义先,胡正名.用于序列密码的布尔函数计数问题.通信学报,1992,13(4): 18-24.[3]郭宝安.非线性序列的分析与综合:[博士论文].北京:北京邮电大学,1993.6.[4]Yang Yi Xian, Guo Bao An. Further enumeration of Boolean functions of cryptographic significance. J.Cryptology, 1996, 8(1):115-122.[5]OConnor L, Klapper A. Algebraic nonlinearity and its applications to cryptography. J. Cryptology, 1994, 7(4): 213-227.[6]杨义先,林须端,胡正名.编码密码学.北京:人民邮电出版社,1992年,第15章,538-549.

施引文献

期刊类型引用(0)
其他类型引用(1)

资源附件(0)

访问统计

计量

文章访问数: 2030
HTML全文浏览量: 90
PDF下载量: 406
被引次数: 1

1. 引言
2. 脑出血MIT检测原理
3. 改进的Newton-Raphson图像重建算法
4. 建模与仿真
4.1 脑出血MIT有限元模型的建立
4.2 仿真计算
5. 实验及结果分析
6. 结束语

姓名
邮箱
手机号码
标题
留言内容
验证码

留言板

一阶相关免疫布尔函数的计数

计量