CN114144109A

CN114144109A - 卒中监测

Info

Publication number: CN114144109A
Application number: CN202080018447.5A
Authority: CN
Inventors: 阿明·阿布什; 阿尔曼·阿夫萨里
Original assignee: Meishi Medical Instrument Co ltd
Current assignee: Meishi Medical Instrument Co ltd
Priority date: 2019-03-04
Filing date: 2020-03-04
Publication date: 2022-03-04
Also published as: EP3934531C0; JP7465888B2; AU2020230859A1; EP3934531A4; WO2020176940A1; US20220079443A1; EP3934531A1; JP2022522894A; EP3934531B1; CA3131076A1; KR20220112176A

Abstract

一种用于在关键康复期期间连续监测脑卒中的计算机实现的方法，该方法包括以下步骤：(i)访问初始图像数据，该初始图像数据表示包含卒中区域的受试者的脑的初始图像；(ii)访问散射参数数据，该散射参数数据表示由受试者的脑散射并源自围绕受试者的脑设置的多个天线的微波；以及(iii)使用无梯度优化方法处理散射参数数据和初始图像数据，以生成受试者的脑内的卒中区域的空间维度的估计，其中，受试者的脑的初始图像被用作先验信息以提高所生成的估计的准确度，并且卒中区域的空间维度被生成作为无梯度优化方法的全局参数。

Description

卒中监测

技术领域

本发明涉及医学成像，特别地涉及用于在关键康复期期间连续监测脑卒中(stroke)的装置和方法。

背景技术

世界卫生组织将卒中定义为“持续超过24小时或在24小时内因死亡而中断的脑血管原因的神经功能缺损”。这表明需要快速行动，使得能够施用治疗来降低卒中的严重程度，要知道在症状发作之后约6小时后，由此导致的脑功能丧失变得不可逆转。然而，这样的快速行动目前是在部分医学盲区的情况下进行的，原因是目前没有卒中监测工具来有效地测量卒中“几何形状”及其对治疗的响应。

卒中发生后的前六小时通常被称为“关键康复期”(“CRP”)，其被定义为“由卒中导致的功能丧失可以恢复的时段”。因此，一旦诊断出卒中，就立即进行卒中的医学治疗。然而，特别是在CRP期间，医学治疗需要包括基本上连续的监测，以评估药物或其他治疗对卒中的“几何尺寸”的影响，并且如果需要，则提供相应的医疗调整(例如，根据卒中的行为，立即改变为不同的治疗)。然而，目前可用的成像模式目前并没有促进这样的关键需求，对于这，在以下文献中描述了示例：P.D.Schellinger等，“Monitoring intravenousrecombinant tissue plasminogen activator thrombolysis for acute ischemicstroke with diffusion and perfusion MRI”，卒中，第31卷，第6期，第1318至1328页，2000年3月(“Schellinger”)。

卒中在CRP期间的几何变化

为了演示在CRP期间的卒中行为，图1A示出了在症状发作后3小时的MRI得出的卒中模型，类似于Schellinger中的医学情况。该模型是从在紧急情况下使用的“卒中MRI”成像模式所生成的图像得出的。这种成像模式的主要目的仅仅是揭示卒中的严重程度，因此与需要较长计算时间来生成的常规MRI图像相比，这种成像模式通常具有更低的特异性。

图1B和图1C表示在症状发作后5个半小时在CRP期间图1A的卒中的形状的两种可能变化。图1B中所示的卒中的扩大通常被称为“低灌注(hypoperfusion)”，并且其特征在于由于“凝块”的延伸导致对脑组织供氧不足，从而将受影响区域的介电特性降低至低于周围组织的介电特性10％至15％的值。这种情况发生在溶栓治疗(通过向受影响区域注射或直接输送(使用导管)溶解凝块的药物来分解“凝块”的方法)进行得相当晚时——通常是自症状发作的4个半小时之后，并且相反地作用使得卒中区域生长，而不是按预期使卒中区域收缩。

另一方面，图1C演示了通常在卒中发作后的前4小时内及时进行溶栓治疗时的凝块的分解。图1B和图1C的几何变化是对卒中治疗的不同关键响应，但由于缺乏有效的卒中监测方法，目前无法检测到图1B和图1C的几何变化。

在卒中期间，当脑动脉破裂时，涌入周围组织的血液约4分钟后形成凝块。因此，在患者到达时，脑的受影响区域已经被凝块占据。

实时监测的障碍

缺乏连续监测技术源于现有医学成像模式(主要是MRI和X射线)的一些固有限制。例如，MRI的庞大且静态的结构使其无法用作连续监测工具，因为患有严重卒中的患者通常处于紧急医疗阶段，这妨碍在重症监护病房(ICU)与用于监测的成像单元之间反复转运他们。此外，X射线的电离效应和MRI对人体的遗传毒性影响阻止了这些成像模式被用来将身体持续暴露于可能不安全剂量的辐射。

为了解决这些困难，已经引入电磁断层扫描(EMT)作为可能的监测模式。EMT涉及根据通常以本领域中被称为“S矩阵”或“散射参数”的形式的测量数据来重建未知受试者的图像。这种成像模式受益于EMT硬件重量轻，因此可在每个医疗单元(包括ICU)中携带和安装EMT硬件，并且这种成像模式还受益于EMT硬件是安全的，因为EMT硬件不会引起MRI和X射线成像模式的电离和遗传毒性效应。

然而，与这些成熟的技术相比，EMT存在一些局限性。特别地，在通常进行医学EMT的UHF(300MHz至3GHz)和S(2GHz至4GHz)频带中，人头部组织的尖角(特别是如图1A所示的脑和脑脊髓液之间的圆角)处的高度局部化倏逝波(evanescent wave)的存在和衍射效应非常显著。因此，EMT图像在几何上是粗糙的，并且无法准确地表示正在被成像的受试者(例如，卒中区域)的实际形状，但是这对于CRP监测是期望的。

例如，图1D是示出图1A的示例的相对介电常数(ε_r)的空间分布的EMT图像，如使用在1.5GHz下工作的八个成像偶极子天线所生成的。甚至与图1A进行粗略比较表明：计算出的图像没有准确地表示卒中的实际空间分布或介电特性(用σ表示的检索的电导率被排除在演示之外，原因是关于电导率σ的相应的检索准确度通常比关于相对介电常数ε_r的相应的检索准确度更差)。此外，通过EMT生成医学图像所需的相应计算时间相当高。例如，使用具有3.6GHz下的Intel Core^TM i7-4790 CPU和16GB RAM的计算机平台，在15dB信噪比(SNR)下生成图1D的图像的时间约为76分钟。因此，UHF频带和S频带中的倏逝波的存在、衍射效应以及所需的过高计算时间都排除了将常规EMT用于独立卒中监测。如Schellinger所提出的：“仍然需要一种快速的卒中成像工具，这种卒中成像工具对于在最初6小时内检测脑内出血(ICH)具有足够高的灵敏度，并且如果存在这种卒中成像工具的话，这种卒中成像工具能够标识处于危险中的组织”。

因此，期望克服或减轻现有技术的一个或更多个困难或者至少提供有用的替选方案。

发明内容

根据本发明的一些实施方式，提供了一种用于在关键康复期期间连续监测脑卒中的计算机实现的方法，所述方法包括以下步骤：

(i)访问初始图像数据，所述初始图像数据表示包含卒中区域的受试者的脑的初始图像；

(ii)访问散射参数数据，所述散射参数数据表示由受试者的脑散射并源自围绕受试者的脑设置的多个天线的微波；以及

(iii)使用无梯度优化方法处理散射参数数据和初始图像数据，以生成受试者的脑内的卒中区域的空间维度的估计，其中，受试者的脑的初始图像被用作先验信息以提高所生成的估计的准确度，并且卒中区域的空间维度被生成作为无梯度优化方法的全局参数。

在一些实施方式中，通过优化受试者的脑的卒中区域的第一预定介电常数值的空间维度和受试者的脑的非卒中区域的第二预定介电常数值的空间维度来最初确定卒中区域的空间维度。

在一些实施方式中，卒中区域的空间维度和相对介电常数被生成作为无梯度优化方法的全局参数。

在一些实施方式中，通过二维平面中的交叠椭圆来近似卒中区域的形状，并且通过确定交叠椭圆的空间维度来确定卒中区域的空间维度。

在一些实施方式中，交叠椭圆具有短轴，所述短轴具有固定空间维度，并且交叠椭圆的空间维度被确定为与所述交叠椭圆的长轴对应的两个参数。在一些其他实施方式中，通过确定四个几何参数来确定卒中区域的空间维度。

在一些实施方式中，该方法包括在连续时间内重复步骤(ii)和步骤(iii)以监测卒中区域随着时间的增长或收缩。

无梯度优化方法可以是Nelder-Mead无梯度优化方法。

可以通过磁共振成像或x射线成像或电磁断层扫描成像来生成受试者的脑的初始图像。

根据本发明的一些实施方式，提供了一种用于在关键康复期期间连续监测脑卒中的装置，所述装置包括：

存储器；

至少一个处理器；以及

至少一种计算机可读存储介质，所述至少一种计算机可读存储介质上存储有指令，所述指令当由所述至少一个处理器执行时，使所述至少一个处理器执行以下步骤：

(iii)处理散射参数数据和初始图像数据，以估计受试者的脑内的卒中区域的空间维度，其中，受试者的脑的初始图像被用作先验信息以提高确定的准确度，并且卒中区域的空间维度被确定为无梯度优化方法的全局参数。

在一些实施方式中，交叠椭圆中的每一个的空间维度被确定为两个几何参数。

在一些实施方式中，通过确定四个几何参数来确定卒中区域的空间维度。

在一些实施方式中，该装置包括在连续时间内重复步骤(ii)和步骤(iii)以监测卒中区域随着时间的增长或收缩。

无梯度优化方法可以是Nelder-Mead无梯度优化方法。

根据本发明的一些实施方式，提供了至少一种计算机可读存储介质，所述至少一种计算机可读存储介质上存储有指令，所述指令当由脑监测装置的至少一个处理器执行时，使所述至少一个处理器执行上述方法中的任何方法的步骤。

附图说明

在下文中参照附图仅通过示例的方式描述了本发明的一些实施方式，在附图中：

图1A是从包含卒中区域的受试者的头部的截面平面图MRI图像得出的，并且示出了用来发起EMT的卒中模型信息，其中在1.5GHz下的相应组织介电特性(ε_r，σ)取自Gabriel；

图1B和图1C示出了根据Schellinger图1A的卒中在CRP期间的两个相应的可能关键响应；

图1D是受试者的头部的图像，该图像由基于鲁棒的基于梯度的优化的现有技术EMT方法生成并且没有由卒中MRI图像发起，而且该图像显然不能准确地标识卒中区域；

图1E和图1F是由与针对图1D的方法相同的方法生成的相应的EMT图像，但是该图像由图1A的卒中MRI图像发起，图1E准确地标识卒中区域，图1F示出了在每个像素处其对药物治疗(溶栓)的后续响应。然而，这种改进的EMT检索方法花费多于一个小时，从而使得该方法对于在CRP期间进行卒中监测而言不切实际；

图2示出了被表示为针对三个卒中参数a₀、b₀和ε_r(几何参数a₀和b₀在图1A中被示出)的3D参数空间中的几何操作的Nelder-Mead(“NM”)优化的不同优化操作；

图3是根据本发明的所描述的实施方式的NM优化方法的流程图；在构造顶点之后，该方法首先检查“最好”顶点是否满足截断条件。如果“最好”顶点不满足截断条件，则执行第一误差减少操作，即“反射”。因此，如果反射顶点相对于最好顶点具有更低的误差，则执行扩展操作以希望找到甚至更好的顶点。如果反射顶点相对于最好顶点不具有更低的误差，则根据反射顶点相对于最差顶点和第二最差顶点的误差水平，执行向前/向后紧缩和收缩以最终找到新顶点，该新顶点在每次迭代中最起码比“最差”顶点更好(具有较低的误差)；

图4A至图4C是表示本文中描述的NM优化方法的检索方法和收敛行为的相应图；

图5A至图5D是示出在具有8个天线并使用两个卒中形状参数a₀、b₀的情况下在CRP期间的卒中监测中一些影响因素对NM优化的准确度和计算时间的影响的相应图；图5A示出了当SNR＝15dB时天线总数的影响，图5B示出了SNR的影响，图5C示出了当SNR＝15dB时卒中尺寸的影响，图5D示出了当SNR＝15dB时每个单独组织与Gabriel中的数据库之间的差异的影响；

图6：3D MRI得出的卒中模式的NM优化的检索方法，其中：

图6A是具有围绕患者的头部设置的圆形天线阵列的三维EMT几何结构的表示；

图6B示出了通过本文中描述的三维EMT方法成像的患者头部内的卒中区域的特写图；

图6C示出了具有定义三维卒中区域的四个几何参数的成像的卒中区域；

图6D示出了在凝块分解期间成像的三维卒中区域；

图6E是通过该方法确定的四个几何卒中区域参数和相对介电常数随着迭代次数变化的图；

图6F是示出方法的收敛(以dB为单位)随着迭代次数变化的图；

图7是根据本发明的实施方式的卒中监测装置的示意图；以及

图8是根据本发明的实施方式的卒中监测方法的流程图。

具体实施方式

EMT作为监测模式的可能性促使发明人开发新的EMT方法，该新的EMT方法适用于在CRP期间进行卒中监测。实际上，尽管上述EMT的一些固有限制——即，在UHF频带和S频带处的倏逝波的存在和衍射效应——是不可避免的，但发明人确定现有技术的EMT方法的过长的计算时间是由于在每个现有技术EMT系统中使用的基于梯度的优化的数值公式而导致的。由于这些优化在所得到的图像的每个像素处寻找“变量”的最佳值，因此通常需要高计算时间来寻找这些最佳值。换言之，如果脑内卒中区域的空间分布和介电特性是未知的，则它们被认为是变量，它们的最佳值必须在图像的每个“像素”处被检索(例如图1A中所示的那样)。由于具有临床上可接受的空间分辨率需要相对高数目的像素(即使当考虑EMT的有限的最大可实现的空间分辨率时)，因此现有技术的EMT计算方法对于CRP期间的卒中监测是非常低效的。此外，如从图1D中明显的是，由这些低效方法计算出的得到的图像在任何情况下不够准确。

考虑到这一点，发明人认识到，由于患者在到达时立即被转移到成像单元，因此在将患者转移到ICU之前(使用卒中MRI成像模式生成卒中区域的初始图像花费不到一个小时)，使用诸如MRI(或X射线或甚至EMT)的高分辨率成像模式从早期诊断阶段获得卒中区域的初始空间分布作为先验信息，如图1A所示。虽然MRI(或X射线或EMT)成像由于上述原因(安全性、成本和计算时间)实际上不能用于连续监测，但是如在以下文献中描述的那样可以对准或“配准”在患者到达时生成的卒中MRI(或X射线或EMT)图像：G.Boverman,C.E.L.Davis，S.D.Geimer和P.M.Meaney，“Image registration for microwavetomography of the breast using priors from nonsimultaneous previous magneticresonance images”，医学和生物学中的电磁学、射频和微波的IEEE期刊(IEEE Journal ofElectromagnetics,RF and Microwaves in Medicine and Biology)，第2卷，第1期，第2至9页，2018年3月；以及R.L.Leijsen，W.M.Brink，C.A.T.van den Berg，A.G.Webb和R.F.Remis，“Three-dimensional contrast source inversion-electrical propertiestomography”，IEEE医学成像汇刊(IEEE Trans.Med.Imag.)，第37卷，第9期，2018年9月(“Leijsen”)，其中，Leijsen给出了不同组织的介电特性。

如图1E和图1F所示，(在该示例中图1A的)配准图像可以在ICU中安装的另外的常规EMT装置中被采用(实际上，MRI仅被执行一次)，然后被用作先验信息以生成卒中响应的图像(在图1B和图1C中示意性的被示出)，表明了检索准确度的显著提高。提高的原因在于：配准的卒中MRI图像被用于限制由基于梯度的优化生成的可能的可检索图像的空间，从而使EMT收敛到更准确的图像。然而，由于常规的(并因此基于梯度的)EMT的基于像素的优化特性，该方法的计算时间仍然过高，对于该示例，该方法的计算时间约为一小时。与没有MRI图像补充EMT(76分钟)的图1D的示例相比，这仅仅是时间上的轻微减少。因此，如果基于梯度的EMT由MRI、X射线或CT扫描图像发起，则基于梯度的EMT可能仅用于非紧急情况情形，例如肿瘤筛查(如Leijsen所述，对于3D问题，在标准计算机上，报告的计算时间约为11小时)。

本文描述了用于在关键康复期(CRP)期间连续监测脑卒中的装置和方法，为了方便起见，所述装置和方法在本文中也被称为卒中监测装置和方法。在本文中描述的卒中监测方法和装置中，MRI(或X射线或EMT)配准的图像被实现为先验信息以提供卒中区域的初始几何形状和介电特性。然而，代替计算成像区域的每个像素处的变量的低效的现有技术的基于梯度的EMT方法，计算卒中区域的形状、维度和介电特性作为全局参数。因此，如下所述，在稍后的监测时间通过无梯度优化方法连续地更新这些参数的全局值。

例如，在所描述的实施方式中，参数是将卒中区域的形状和维度定义为具有1cm的固定短轴的两个交叠椭圆的半长轴的两个几何参数——如图1A所示，即a₀和b₀(初始值为a₀＝1.5cm，b₀＝2cm)——以及定义卒中区域的相对介电常数即ε_r的第三参数。将卒中的形状和介电特性视为全局参数将未知量的数目从例如图1A中图像的约1000个基于像素的变量减少到仅3个全局参数。这种重新组成使无梯度优化能够被用于EMT作为用于参数优化的高效技术(因为参数是全局优化的，而不是在每个像素处局部优化的)。显然，这种无梯度优化不能提供基于像素的图像，例如图1D至图1F的那些图像。然而，对于CRP期间的医疗需要，每两分钟快速且准确地跟踪卒中参数(例如，a₀、b₀和ε_r)比通过基于梯度的优化以每小时时间为基础获得卒中的粗略图像在医学上更具相关性。

开发了无梯度优化方法来解决基于梯度的优化不适用的问题，特别是在优化期间要最小化的函数不可微分或不平滑的情况下。

已经针对不同的电磁和天线应用开发了各种无梯度优化方法，包括Nelder-Mead(“NM”)、遗传算法(“GA”)和粒子群优化(“PSO”)方法。在计算时间是优化方法中除了准确度之外的主要考虑因素的情况下，NM优化方法通常是最快的无梯度优化方法。在CRP期间的卒中监测的情况下，时间就是生命，发明人认为NM优化最适合监测要求(尽管在其他实施方式中可以使用其他无梯度优化方法)。由于以前未将NM优化用于EMT应用，因此下面描述用于该特定应用的NM优化方法的概述。在以下文献中描述了一般的NM方法：N.Pham，A.Malinowski和T.Bartczak，“Comparative study of derivative free optimizationalgorithms”，IEEE工业信息汇刊(IEEE Trans.Industr.Inform.)，第7卷，第4期，第592页至600页，2011年11月(“Pham”)。

在所描述的实施方式中，如图7所示，由卒中监测装置执行所描述的方法。在使用中，该装置与矢量网络分析器(VNA)或收发器701进行通信，该矢量网络分析器或收发器701又连接到微波天线阵列705。

如所示的，微波天线阵列701被布置成接收其脑要被成像的患者的头部704，使得阵列中的每个天线可以被选择性地激励，以将电磁波或微波频率的信号辐射到受试者的头部中并通过受试者的头部以被散射，并且相应的散射的信号由所有天线阵列(包括发射相应信号的天线)检测。

对于本领域技术人员明显的是，矢量网络分析器(VNA)701如上所述激励天线，并将来自天线的相应信号记录为数据(在本文中被称为“散射”数据)，该数据以在本领域中被称为“散射参数”或“S参数”的形式表示散射的微波的幅度和相位。VNA 701将该数据发送到用于处理的装置以生成关于成像的对象的内部特征(例如，脑凝块、出血部位和其他特征)的信息。在所描述的实施方式中，具有大于700dB的大动态范围和低于-700dBm的本底噪声(noise floor)的VNA可以被用来激活天线，以在0.5GHz至4GHz的频带上发射电磁信号并从那些天线接收散射的信号。

尽管所描述的实施方式的装置是计算机的形式，但是在其他实施方式中也不必是这种情况。如图7所示，所描述的实施方式的卒中监测装置是64位Intel架构计算机系统，并且由卒中监测装置执行的卒中监测方法被实现为存储在与计算机系统相关联的非易失性(例如，硬盘或固态驱动器)存储装置704上的一个或更多个软件模块702的编程指令。然而，明显的是，这些方法的至少一部分可以替选地以一种或更多种其他形式被实现，例如被实现为现场可编程门阵列(FPGA)的配置数据，或被实现为诸如专用集成电路(ASIC)的一种或更多种专用硬件部件，或被实现为这样的形式的任何组合。

卒中监测装置包括随机存取存储器(RAM)706、至少一个处理器708和外部接口710、712、713、714，所有这些都通过总线716互连。外部接口包括：网络接口连接器(NIC)712，其将卒中监测装置连接到诸如因特网720的通信网络；通用串行总线(USB)接口710，通用串行总线(USB)接口710中的至少一个可以连接到键盘718和诸如鼠标719的定点设备；以及显示适配器714，其可以连接到诸如LCD面板显示器722的显示设备。

卒中监测装置还包括诸如Linux或微软视窗软件(Microsoft Windows)的操作系统724，并且在一些实施方式中包括附加软件模块726至730，包括：web服务器软件726，例如Apache，可在http://www.apache.org处获得；脚本语言支持728，例如PHP，可在http://www.php.net处获得，或Microsoft ASP；以及结构化查询语言(SQL)支持730，例如MySQL，可从http://www.mysql.com获得，其允许数据被存储在SQL数据库732中以及从SQL数据库732检索数据。

web服务器726、脚本语言模块728和SQL模块730一起向卒中监测装置提供一般能力，以允许具有配备有标准web浏览器软件的标准计算设备的远程用户访问卒中监测装置，尤其是在CRP期间监测卒中的进展。A.EMT的NM优化：初始化

对于图1B和图1C(由图1A发起)中所示的EMT问题，八个成像偶极子天线以圆形配置围绕受试者的头部被布置，使得受试者的头部的2D MRI模型穿过偶极子天线的端口(因为在优化方法中实现的S矩阵在偶极子端口处被记录)。

天线以1.5GHz的正弦电磁波照射头部，比如说症状发作之后的5个半小时，此时已经施用了溶栓。该计时示例取自Schellinger描述的卒中病例，其中受试者的到达时间为症状发作后约3小时，并且在到达之后0.75小时准备2D单切片卒中MRI。然后，天线记录以S矩阵形式的相应散射的场(下面讨论天线数目对准确度和检索时间的影响)。然后该S矩阵被实现成更新三个卒中参数a0、b0、ε_r，三个卒中参数在图1A中的初始值为X₀＝(a₀＝1.5cm，b₀＝2cm，ε_r＝39)，其可以被认为定义3D空间(因为参数的数目为3)中的点(或如下所述的顶点)。

出于以下原因而选择这些参数。由于通过受试者的脑的卒中区域的相对介电常数ε_r的变化来反映施用至卒中的任何药物的效果，因此相对介电常数ε_r可以被认为是监测卒中并在稍后时间收敛至与S矩阵匹配的值的可靠参数。如上所述，由于与相对介电常数ε_r相比，卒中区域的电导率σ的检索准确度通常较差，因此从检索中排除电导率。此外，几何参数或形状参数a₀、b₀被选择成最好地反映卒中区域的几何变化(延伸或收缩)。对于图1A所示的平滑的卒中形状，a₀、b₀是两个椭圆的半长轴，它们被认为相对于可能的替选几何参数最准确地表示卒中区域的几何变化(下面讨论更复杂的几何形状)。发明人已经发现，这些形状参数也连续不断地收敛至最适合相应S矩阵的值。图1A中的2D物理问题域与图2A中构造的参数的3D数学空间不同。

除了X₀之外，三个参数的变化范围限定了可用的参数空间，假设遭受低灌注的卒中区域(图1B)可以延伸到脑的整个左半球，并且成功的凝块分解(图1C)可以完全去除凝块(即，使卒中区域收缩为零)。这个大的变化范围对应于0≤a₀≤4cm、0≤b₀≤7cm且39≤ε_r≤43.5的参数空间，表明即使没有关于对药物的卒中响应的可靠预测，该优化方法在寻找最佳顶点时也不需要非常精确的参数范围。如图2A所示，为了优化这3个参数，NM优化首先在参数的3D空间内构造具有3+1个顶点的等长单纯形(广义三角形)。实际上，顶点之一是由初始MRI图像提供的初始猜测(即，起始点)X₀。如下通过将以下矢量添加到该初始猜测来得出具有相等距离c的所有其他顶点：

其中，矢量分量为：

并且其中，N是参数的数目(即，在所描述的实施方式中为3)。通常，c＝1允许该方法在初始步骤处在足够大的体积中进行搜索。小的c值通常需要长的计算时间，并且可能误导该方法以使该方法仅找到局部最小值。此外，因为噪声的存在可能导致稍微不同的顶点产生相同的频率响应(S矩阵)，所以使顶点彼此位于足够远处(相距c≥1)使得方法在早期迭代时对噪声具有鲁棒性。

通过构造单纯形，下一步骤是评估每个顶点X₀、X₁、X₂、X₃处的目标函数(即，通过NM优化要被最小化的函数)。在图1B或图1C的EMT问题中，该目标函数被如下定义为测量的(“meas.”)S矩阵与检索的(“retr.”)S矩阵之间的L₂范数失配：

其中，Ω表示成像域(在柱面坐标

中)，χ是被定义为

的人头部的介电特性的对比，其中，

是复介电常数，

是自由空间介电常数，ω＝2π×1.5GHz是角频率。在所描述的具有八个天线的示例中，S矩阵的大小因此是8×8；因此失配是S^meas.与S^retr.的相应矩阵元素之间的差。如从以下文献中所描述的得出的，E是图1B或图1C的成像域上的总电场：W.C.Chew，“Waves and Fields inInhomogeneous Media”，IEEE出版社，纽约，1995(“Chew”)的第9章以及A.Afsari，A.Abbosh和Y.Rahmat-Samii，“Modified Born iterative method in medical electromagnetictomography using magnetic field fluctuation contrast source operator”，IEEE微波理论技术汇刊，DOI：10.1109/TMTT.2018.2876228，并且E由以下式子给出：

其中，E^inc是在Ω中不存在任何对象时的入射电场。最后，G是Chew的第1章中给出的并矢格林函数(dyadic Green’s function)，并且

是自由空间波数，其中，μ₀＝4π×10^-7是自由空间磁导率。

如图2B所示，在评估所有顶点的目标函数之后，三个顶点在NM技术中具有特殊重要性。给出目标函数的最小值的顶点(在所描述的实施方式中由等式(3)给出)被称为“最好”顶点并且被表示为X_b。具有最高L₂范数失配的顶点被称为“最差”顶点X_w。为了将该方法引导到总是减小等式(3)中的失配的方向，该方法确定等式(3)具有第二最高值的辅助点。该顶点被称为“第二最差”顶点X_sw。然后除X_w之外的所有顶点的平均值被计算为X_a。X_a与X_w之间的线段(L_XwXa)总是向下的，并且包含一些有用的点，其中等式(3)在优化方法期间可能具有比X_w更低的失配。根据该向下方向搜索最佳点，该方法也被称为“下坡”优化。

在以上初始化之后，该方法执行至少两个且至多五个“错误减少”操作。在这些步骤的每一个中，X_w的旧值从计算机存储器中被移除(即，不被存储用于下一次迭代)，并且所有其他顶点被重新布置以提供X_b、X_w、X_sw的新值。在下文中，介绍以下操作：

反射：NM方法中的第一优化步骤是如下以相同的长度在L上反射最差的顶点X：

X_r＝2X_a-X_w (5)

如图2c所示。该操作是检查在

方向上的移动是否应当由该方法继续还是另一方向将该方法引导至最佳顶点。如果等式(3)在X_r处相对于X_b具有较低的失配，即F(X_r)＜F(X_b)，则在替代之前：

该方法评估找到甚至更好的顶点的机会(其中，编程约定A←B表示旧值A被新值B代替)。为此，总是通过在同一方向

上进一步移动来执行扩展操作。

扩展：根据图2d，该方法沿着

以相同步长进一步移动，即

X_e＝2X_r-X_a (7)

然后，在这个扩展顶点处也评估等式(3)。如果它的值低于X_b(即使它比X_r更差)，即F(X_e)＜F(X_b)，则该方法进行如下替代：

并且(迭代地)返回到反射步骤。尽管它是在其他顶点中找到的最好顶点但该方法不立即接受X_b←X_r的原因来自以下事实：因为该顶点位于由X_e形成的新单纯形的内部，因此该顶点被该方法保留。因此，通过执行扩展，X_r的邻域仅被保护为子域，其中可能存在一些其他好的顶点或甚至更好的顶点以使等式(3)最小化。然而，如果等式(3)在X_e处不低于X_b，则执行等式(6)的代替，并且该方法迭代地返回到第一操作(即，反射)。

向前紧缩：等式(6)或等式(8)假设F(X_r)＜F(X_b)。如果这没有被实现，但F(X_sw)＜F(X_r)＜F(X_w)，则该方法已经沿

方向过度地被移动，并且更好的顶点可能位于比X_r更近的距离处。因此，如图2E所示，通过从X_r向X_a返回半个步长(通常为平衡的步长)来执行沿

的向前紧缩，即

X_fc＝0.5X_a-0.5X_w (9)

如果F(X_fc)＜F(X_r)，则通过返回至初始化步骤并将这些顶点从最差顶点重新布置到最好顶点，在顶点X、X_b、X_sw、X_fc上形成新的单纯形(simplex)。

向后紧缩：如图2f所示，如果F(X_w)＜F(X_r)，则

可能仍然包含可以改进等式(3)的失配的一些顶点，并且如下以与向前紧缩相同的方式，通过从X_a向X_w向后移动半个步长来执行向后紧缩：

X_bc＝0.5X_a+0.5X_w (10)

如果F(X_bc)＜F(X_w)，则通过返回至初始化步骤并将顶点从最差重新布置到最好，在X、X_b、X_sw、X_bc上构造如图2f所示的新单纯形。

收缩：尽管如此，如果上述条件都没有发生，则找到朝向最佳顶点的更好方向的最后步骤是使单纯形收缩。为此，仅保持最好顶点X_b，并且对于其他顶点，如下执行收缩操作(对于每个第i顶点)：

X_i(新)＝0.5X_b+0.5X_i(旧) (11)

然后，该方法返回到初始化步骤以重新布置在图2g所示的收缩步骤中形成的新顶点。继续NM优化的该迭代过程，直到F(X_b)满足目标函数的截断条件，在所描述的实施方式中截断条件为：

F(X_b)＜10^-7 (12)

将满足截断条件的X_b的值存储为最终结果。等式(12)中的截断条件被选择得非常小，以确保满足检索参数所需的准确度。截断条件的较大值不会导致非常准确的参数检索。为了一次演示所有这些步骤，图3是用于EMT的NM优化方法的流程图，并且下面给出了相应的伪代码。

I.实际中的NM无梯度优化：2D检索

为了有效地监测被称为低灌注且在图1B中示出的不同卒中行为或者图1C中的凝块分解，利用上述NM优化方法来优化三个参数a₀、b₀和ε_r。仅选择形状参数a₀和b₀进行优化并排除诸如图1B或图1C中的卒中的中心的坐标的位置参数的原因来自以下事实：在低灌注或凝块分解过程中，卒中位置不改变。相反，卒中仅围绕其中心延伸或收缩。

A.NM优化性能

对于低灌注，给定示例中的参数的实际值为a₀＝2.5cm、b₀＝3cm、ε_r＝39，而检索到的值——在图4A中示出了其在连续迭代中的演变——为a₀＝2.68cm、b₀＝3.32cm、ε_r＝39.08。如图5A所示，当单纯形在迭代期间变小时，NM方法变化逐渐变得稳定。对于图1C的凝块分解，NM优化方法的实际值和检索值分别为a₀＝1cm、b₀＝1cm、ε_r＝40以及a₀＝0.997cm、b₀＝1.002cm、ε_r＝40.07(图4B)。应当注意的是，在凝块分解期间，受影响区域的相对介电常数在恢复时再次接近脑的介电常数。在这两种情况下，尽管模拟环境是高噪声的(S矩阵的SNR被设置为15dB)，但是参数仍然被准确地更新。为了提供对这种准确度的定量分析，“检索误差”被定义为：

在每种情况下，检索误差都很低：对于低灌注结果，检索误差为-18.86dB；而对于凝块分解结果，检索误差为-27.95dB。

然而，关键点是检索时间。根据图4C，一旦满足等式(12)的终止条件，则方法终止。在所描述的示例中，低灌注和凝块分解检索分别需要94次迭代和69次迭代，并且每次迭代花费1.25秒。因此，在每两分钟时间范围内，卒中被连续监测以评估治疗方法并在需要时立即准备开颅手术(手术移除受试者的颅骨的一部分以直接进入受影响区域)。

图4A和图4B提供了NM方法分别仅在前70次迭代和前40次迭代中的性能的图形说明(在进一步迭代中检索到的参数的变化太小而不能被有效地绘制在同一图中)。与凝块分解过程相比，由低灌注针对初始拓扑(即，图1A)引起的更大的几何变化导致与凝块分解检索相比更高的检索误差和更长的收敛时间。

A.对准确度和计算时间的影响因素

在影响基于NM优化的EMT的准确度(检索误差)和计算时间的不同因素中，主要贡献因素是：成像天线的总数、SNR、形状参数a₀、b₀以及每个受试者的头部组织的介电特性相对于在以下文献中描述的数据库的微小差异，这些因素被用来配准初始卒中MRI图像；所述文献如下：C.Gabriel，S.Gabriel和E.Corthout，“The dielectric properties ofbiological tissues:I.Literature survey”，物理医学生物学(Phys.Med.Biol.)，第41卷，第1期，第2231页至2249页，1996年(“Gabriel”)。为了描述这些影响，由于低灌注的检索更加具有挑战性(具有更高的检索误差)，因此针对这种医学情境来描述天线数目、SNR和Gabriel的准确度水平的影响。

图5A示出了：以增加每次迭代的计算时间为代价，将天线的数目增加到特定水平显著地提高准确度。这种提高来自以下事实：具有附加天线相当于具有更多信息。由于相同的原因，迭代的总数从94次减少到78次，原因是可能的解的空间现在被更多地限制并且该方法需要更少的迭代来满足截断条件。然而，由于目标函数的大小的增加，每次迭代的计算时间急剧增加，目标函数的大小对应于S矩阵的大小。实际上，具有M个成像天线会产生M×MS矩阵。在该文献中，8个被建议作为以可接受的准确度提供头部的医学图像的最少天线数目。

要考虑的第二因素是SNR。当测量的数据被噪声高度污染时，彼此接近的顶点(在图2A中具有接近的参数值)具有最可能被噪声主导的非常类似的频率响应。因此，该方法的准确度会降低，原因是它不再能够区分接近的顶点。如图5B所示，如果医学成像的域与噪声高度隔离，则改进了检索准确度和所需的迭代次数二者。这是因为等式(3)中的测量的S矩阵包含用来检索卒中响应(低灌注)的更可靠的信息，因此使等式(3)中的目标函数最好地最小化的顶点更接近于实际的顶点。每次迭代的计算时间保持不变，即1.25秒，原因是该时间是成像天线的数目和检索参数的函数。

要考虑的下一个因素是可以准确地被检索的卒中区域的尺寸范围。为此，根据经验，引入截面因数a₀×b₀来表示卒中尺寸。图5c示出了以下：如果a₀×b₀＜0.64cm²或者a₀×b₀＞9cm²，则对于图1B中的示例，检索方法变得不准确。如上所述，衍射效应和倏逝波的存在在UHF频带和S频带中是显著的，并且当卒中区域的尺寸变得小于四分之一波长时，这些效应变得极具破坏性。当在1.5GHz下的波长在凝块内为3.2cm时(其中ε_r＝39)，四分之一波长为0.8cm，得到0.64cm²的截面因数，低于该截面因数，检索准确度会显著降低。另一方面，对于非常高的截面因数，例如，a₀×b₀＞9cm²，卒中形状的变化不能被准确地跟踪，原因是该变化不可避免地需要附加的形状参数(除a₀、b₀之外)以对卒中的不规则性和非对称性更准确地建模，特别是在卒中通过颅骨从一侧变形但仍然可以从另一侧延伸的情况下。从计算时间的角度来看，如果对于大尺寸的卒中区域仅采用两个形状参数，则方法针对图1B中的卒中需要更多的迭代以使卒中的初始形状部分地符合延伸的形状，如图5c所示。这是由以下事实造成的：添加每个形状参数及其相应的变化范围将更加限制可能的解的空间。仅处理两个形状参数确实需要更多的迭代以部分地匹配延伸的卒中的不规则性和非对称性(这在下面被说明)。然而，要注意的是，增加卒中参数的数目作为改进检索误差和减少总迭代的解决方案并不一定会减少计算时间，原因是向问题增加每个参数会向参数空间增加另一个维度。因此，每次迭代的计算时间呈指数增长。

最后一个因素是Gabriel中的数据库与每个单独的人的头部的介电特性之间的一致性水平的影响。由于无梯度优化方法仅集中于检索全局卒中参数，因此确保Gabriel中的介电特性(通过其将卒中MRI配准到EMT装置中)高度准确以避免与期望值的显著检索偏差是非常重要的。幸运的是，这种一致性通常非常高，原因是构造图1A中的组织的材料对于不同的个体是相同的并且这些材料与Gabriel中的数据库之间的微小差异主要来自部分统计随机性。考虑到这个因素，图5d对图1A中的每个组织应用了高达5％的随机差异——相对于在Gabriel中报告的每个组织的相应值——以证明该方法针对不同头部组织的介电特性的微小差异的鲁棒性。如所看到的，计算速度和检索误差二者略有下降，但仍然落在CRP的可接受的要求内。

通过研究影响因素对医学EMT应用中的NM优化方法的准确度和计算时间的影响，该方法在下面被应用于更复杂的问题，其中，根据Leijsen的图4，3D MRI得出的卒中模型补充NM方法，如图6所示。在3D中实现四个形状参数以监测更复杂的凝块分解过程。

I.实际中的NM无梯度优化：3D检索

对于真实的3D EMT问题，如果使用诸如有限元建模(FEM)的发展良好的数值方法来模拟域内的总电场或所检索的S矩阵，则可以进一步提高检索到的参数的准确度。这可以通过直接求解波动方程并考虑成像天线的整个三维物理结构(参见图6A)而不是使用上面给出的基于点源的等式(3)和(4)来实现。实际上，实现上面给出的等式(3)的基于点源的目标函数的原因在于：在EMT中，在现有技术的基于梯度的优化方法中使用相同的目标函数。因此，这样的点源近似为比较这两种方法提供了更好的基础。然而，由于图6所示的问题的真实天线结构的额外复杂性，使用商业COMSOL软件包进行模拟。如图6B所示，整个头部模型内的卒中区域被确定为由三角形网格元素限定的复杂形状，并且卒中参数在图6C中被示出为具有初始值X′₀＝(a₁＝0.41cm，b₁＝0.4cm，c₁＝0.41cm，d₁＝0.42cm，ε_r＝39)，这些值是从卒中区域的中心点到该区域的周界上的四个最远点的距离。

如在图6d中看到的，当及时进行溶栓治疗时，凝块分解方法开始。该方法会将卒中参数轻微地或显著地更改为X₁＝(a₁＝0.4cm，b₁＝0.395cm，c₁＝0.39cm，d₁＝0.41cm，ε_r＝42)。

根据图6E的图，检索到的卒中区域形状参数为X₁＝(a₁＝0.401cm，b₁＝0.395cm，c₁＝0.392cm，d₁＝0.411cm)并且检索到的卒中区域的介电常数为ε_r＝41.43，与它们的实际值非常一致。如图6D所示，由于凝块分解引起的卒中区域的空间维度的变化针对形状参数c₁、d₁最大，这就是为什么它们需要额外的迭代来接近它们的实际值的原因。图6F是NM优化方法的收敛图。由于具有足够的形状参数，因此需要少于50次的迭代来达到截断条件。在安装在相同PC上的COMSOL多重物理量(Multiphysics)中，由于3D问题域的大尺寸，每次迭代需要15秒才能完成。总之，可以每11分钟监测卒中的行为(即，所描述的装置可以每11分钟产生如图6E的一组新图)，这是在CRP期间进行3D监测的合理时间。将Pham中描述的结果视为基本事实，当使用GA或PSO无梯度方法来代替NM方法时，需要多至少10倍的计算时间，例如，对于该特定问题，计算时间约为110分钟。该时间接近现有技术的基于梯度的优化方法的计算时间，并且通常超过CRP时间范围。因此，发明人认为，在无梯度方法和在编写时可用的通用计算机硬件中，仅基于NM无梯度优化的EMT系统可以在CRP期间被有效地用作2D或3D卒中监测工具。

本文中描述并基于Nelder-Mead无梯度优化的EMT监测方法和装置提供了在CRP期间监测卒中的扩展或紧缩的能力，因此可以潜在地被认为是增加卒中存活机会的转化医学进展。本文中描述的结果表明：该方法在通用计算机平台上在每2分钟内检索2D卒中响应或者在每11分钟内检索3D卒中响应是非常高效的，而诸如GA或PSO的其他无梯度方法可以产生相同的输出，但是代价是更长的计算时间。所描述的方法可以由从早期诊断获得的卒中MRI数据发起。然后，卒中区域的形状和介电特性(介电常数的实部)被定义为全局参数。随后，可以使用少量的成像天线通过如本文中描述的便携式EMT系统对患者进行连续成像，并且利用由这些天线在每个成像步骤处记录的S矩阵来更新全局参数，从而标识卒中区域的扩展或紧缩，特别是响应于一种或更多种治疗。因此，所描述的方法和装置可以改善治疗过程，并因此改善卒中受害者的存活机会。

在不脱离本发明的范围的情况下，许多修改对于本领域技术人员将是明显的。

Claims

1.一种用于在关键康复期期间连续监测脑卒中的计算机实现的方法，所述方法包括以下步骤：

(ii)访问散射参数数据，所述散射参数数据表示由所述受试者的脑散射并源自围绕所述受试者的脑设置的多个天线的微波；以及

(iii)使用无梯度优化方法处理所述散射参数数据和所述初始图像数据，以生成所述受试者的脑内的卒中区域的空间维度的估计，其中，所述受试者的脑的初始图像被用作先验信息以提高所生成的估计的准确度，并且所述卒中区域的空间维度被生成作为所述无梯度优化方法的全局参数。

2.根据权利要求1所述的方法，其中，通过优化所述受试者的脑的卒中区域的第一预定介电常数值的空间维度和所述受试者的脑的非卒中区域的第二预定介电常数值的空间维度来最初确定所述卒中区域的空间维度。

3.根据权利要求1或2所述的方法，其中，通过二维平面中的交叠椭圆来近似所述卒中区域的形状，并且通过确定所述交叠椭圆的空间维度来确定所述卒中区域的空间维度。

4.根据权利要求3所述的方法，其中，所述交叠椭圆具有短轴，所述短轴具有固定空间维度，并且所述交叠椭圆的空间维度被确定为与所述交叠椭圆的长轴对应的两个参数。

5.根据权利要求1或2所述的方法，其中，通过确定四个几何参数来确定所述卒中区域的空间维度。

6.根据权利要求1至5中任一项所述的方法，包括在连续时间内重复步骤(ii)和步骤(iii)以监测所述卒中区域随着时间的增长或收缩。

7.根据权利要求1至6中任一项所述的方法，其中，所述无梯度优化方法是Nelder-Mead无梯度优化方法。

8.根据权利要求1至7中任一项所述的方法，其中，所述卒中区域的空间维度和相对介电常数被生成作为所述无梯度优化方法的全局参数。

9.根据权利要求1至8中任一项所述的方法，其中，通过磁共振成像或x射线成像或电磁断层扫描成像来生成所述受试者的脑的初始图像。

10.一种用于在关键康复期期间连续监测脑卒中的装置，所述装置包括：

存储器；

至少一个处理器；以及

至少一种计算机可读存储介质，所述计算机可读存储介质上存储有指令，所述指令当由所述至少一个处理器执行时，使所述至少一个处理器执行以下步骤：

(iii)处理所述散射参数数据和所述初始图像数据，以估计所述受试者的脑内的卒中区域的空间维度，其中，所述受试者的脑的初始图像被用作先验信息以提高确定的准确度，并且所述卒中区域的空间维度被确定为无梯度优化方法的全局参数。

11.根据权利要求10所述的装置，其中，通过优化所述受试者的脑的卒中区域的第一预定介电常数值的空间维度和所述受试者的脑的非卒中区域的第二预定介电常数值的空间维度来最初确定所述卒中区域的空间维度。

12.根据权利要求10或11所述的装置，其中，通过二维平面中的交叠椭圆来近似所述卒中区域的形状，并且通过确定所述交叠椭圆的空间维度来确定所述卒中区域的空间维度。

13.根据权利要求12所述的装置，其中，所述交叠椭圆中的每一个的空间维度被确定为两个几何参数。

14.根据权利要求10或11所述的装置，其中，通过确定四个几何参数来确定所述卒中区域的空间维度。

15.根据权利要求10至14中任一项所述的装置，包括在连续时间内重复步骤(ii)和步骤(iii)以监测所述卒中区域随着时间的增长或收缩。

16.根据权利要求10至15中任一项所述的装置，其中，所述无梯度优化方法是Nelder-Mead无梯度优化方法。

17.根据权利要求10至16中任一项所述的装置，其中，所述卒中区域的空间维度和相对介电常数被生成作为所述无梯度优化方法的全局参数。

18.根据权利要求10至16中任一项所述的装置，其中，通过磁共振成像或x射线成像或电磁断层扫描成像来生成所述受试者的脑的初始图像。

19.至少一种计算机可读存储介质，所述至少一种计算机可读存储介质上存储有指令，所述指令当由脑监测装置的至少一个处理器执行时，使所述至少一个处理器执行根据权利要求1至9中任一项所述的步骤。