CN113273988B - 基于电流量的电阻抗成像方法、装置及存储介质 - Google Patents

基于电流量的电阻抗成像方法、装置及存储介质 Download PDF

Info

Publication number
CN113273988B
CN113273988B CN202010079013.8A CN202010079013A CN113273988B CN 113273988 B CN113273988 B CN 113273988B CN 202010079013 A CN202010079013 A CN 202010079013A CN 113273988 B CN113273988 B CN 113273988B
Authority
CN
China
Prior art keywords
distribution
grid
current
imaging
electric field
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010079013.8A
Other languages
English (en)
Other versions
CN113273988A (zh
Inventor
马而昉
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xian Jiaotong Liverpool University
Original Assignee
Xian Jiaotong Liverpool University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xian Jiaotong Liverpool University filed Critical Xian Jiaotong Liverpool University
Priority to CN202010079013.8A priority Critical patent/CN113273988B/zh
Publication of CN113273988A publication Critical patent/CN113273988A/zh
Application granted granted Critical
Publication of CN113273988B publication Critical patent/CN113273988B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/05—Detecting, measuring or recording for diagnosis by means of electric currents or magnetic fields; Measuring using microwaves or radio waves 
    • A61B5/053—Measuring electrical impedance or conductance of a portion of the body
    • A61B5/0536—Impedance imaging, e.g. by tomography

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Biomedical Technology (AREA)
  • Medical Informatics (AREA)
  • Biophysics (AREA)
  • Pathology (AREA)
  • Engineering & Computer Science (AREA)
  • Radiology & Medical Imaging (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Physics & Mathematics (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • General Health & Medical Sciences (AREA)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • Investigating Or Analyzing Materials By The Use Of Electric Means (AREA)

Abstract

本发明公开了一种基于电流量的电阻抗成像方法、装置及存储介质。成像方法包括：对成像区域进行网格化，并获取网格参数；在电流密度在每个元素中处处相同的条件下，基于网格参数、成像区域中电流密度分布的无源性、电场分布的无旋性以及成像区域边界电流量条件，确定当成像区域中电导率服从某个分布时流过每个元素边界的电流量；基于每个元素边界的电流量确定成像区域的电流密度分布以及电场强度分布；基于电场强度分布确定成像区域的电势分布；基于电势分布进行图像重构。在同一网格条件下，本发明提高了计算电阻抗成像正问题的精度，从而提高了电阻抗成像的精度。

Description

基于电流量的电阻抗成像方法、装置及存储介质
技术领域
本发明实施例涉及电阻抗成像技术领域，尤其涉及一种基于电流量的电阻抗成像方法、装置及存储介质。
背景技术
电阻抗成像(Electrical Impedance Tomography,EIT)技术基于生物组织的电学特性，是一种可视化的、响应快的、无损的、无辐射的以及价格低廉的理想成像技术。
目前，电阻抗成像的一种常用方法是利用电阻抗成像正问题在不同电导率分布下的解来估计成像区域中的电导率分布。而计算电阻抗成像的正问题通常即是求解以成像区域内的电压分布为基本变量的边界值问题。由于在实际应用中成像区域是不规则的，该边界值问题通常需用数值方法求解。目前常用的数值求解该边界值问题的方法是有限元方法。在有限元方法中，首先，成像区域被网格化成许多个小元素，比如二维成像区域被划分成多个小三角形，三维成像区域被划分成多个小四面体。然后，通过假设在每个小元素中电压的分布满足一定的条件(比如是线性变化的)，来离散化边界值问题，进而得到一个以网格中每个点上电压为未知变量的稀疏线性方程组。最后，通过求解稀疏线性方程组可以得到成像区域中电压分布的一个估计。特别地，还将得到附在成像区域边界上的电极上的电压估计。理论上已经证明，对满足一定条件的成像区域及其网格化，当网格中所有元素的尺寸都趋向于零时，有限元方法得到的近似解将收敛到边界值问题的真实解。特别地，有限元方法得到的电极电压的估计将收敛到其真实值。
在另一方面，为使其解能达到一定的精度，传统的有限元方法通常需要使用很密的网格来划分成像区域。而这会增加有限元方法中稀疏线性方程组的尺寸以及非零元的个数，这将大大加大求解该方程组所需的计算量和计算时间，延长整个求解电阻抗成像正问题的计算时间，进而导致电阻抗成像的效率较低。
发明内容
有鉴于此，本发明的目的是提出一种基于电流量的电阻抗成像方法、装置及存储介质，以解决在给定网格条件下电阻抗成像精度低的问题。
为实现上述目的，本发明采用如下技术方案：
第一方面，本发明实施例提供了一种基于电流量的电阻抗成像方法，包括：
对成像区域进行网格化，并获取网格参数，其中，所述网格参数包括网格节点信息、构成网格的元素的元素信息、元素间的邻近关系信息和成像区域的边界信息；
基于所述网格参数，确定所述成像区域的边界电流量条件以及电导率分布的估计；
在电流密度在每个所述元素中恒定不变且处处相同的条件下，基于所述网格参数、所述边界电流量条件、所述成像区域中电流密度分布的无源性、电场分布的无旋性以及所述成像区域中的电导率分布的估计，确定流过每个所述元素边界的电流量；
基于每个所述元素边界的电流量与所述网格参数，确定所述成像区域网格中的电流密度分布；
基于所述电流密度分布与所述电导率分布的估计，确定所述成像区域网格中的电场强度分布；
基于所述电场强度分布确定所述成像区域网格中的电势分布以及电极电压；
反复更新对电导率分布的估计直至在一电导率分布的估计下计算得到的电极电压与实际测量得到的电极电压的相对误差不超过0.1％；
基于最终更新的电导率分布的估计进行图像重构。
第二方面，本发明实施例提供了一种基于电流量的电阻抗成像装置，包括：
网格化模块，用于对成像区域进行网格化，并获取网格参数，其中，所述网格参数包括网格节点信息、构成网格的元素的元素信息、元素间的邻近关系信息和成像区域的边界信息；
电导率估计确定模块，用于基于所述网格参数，确定所述成像区域的边界电流量条件以及电导率分布的估计；
电流量确定模块，用于在电流密度在每个所述元素中恒定不变且处处相同的条件下，基于所述网格参数、所述边界电流量条件、所述成像区域中电流密度分布的无源性、电场分布的无旋性以及所述成像区域中的电导率分布的估计，确定流出每个所述元素边界的电流量；
电流密度确定模块，用于基于每个所述元素边界的电流量与所述网格参数，确定所述成像区域网格中的电流密度分布；
电场强度确定模块，用于基于所述电流密度分布与所述电导率分布的估计，确定所述成像区域网格中的电场强度分布；
电势确定模块，用于基于所述电场强度分布确定所述成像区域网格中的电势分布以及电极电压；
更新模块，用于反复更新对电导率分布的估计直至在一电导率分布的估计下计算得到的电极电压与实际测量得到的电极电压的相对误差不超过0.1％；
图像重构模块，用于基于最终更新的电导率分布的估计进行图像重构。
第三方面，本发明实施例提供了一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现本发明任一实施例提供的基于电流量的电阻抗成像方法。
本发明的有益效果是：本发明提供的基于电流量的电阻抗成像方法，通过对成像区域进行网格化，在电流密度在每个元素中恒定不变且处处相同的条件下，基于网格参数、边界电流量条件、成像区域中电流密度分布的无源性、电场分布的无旋性以及成像区域中的电导率分布的估计，确定流过每个元素边界的电流量；然后基于每个元素边界的电流量得到成像区域的电势分布；最后基于该电势分布进行图像重构。在同一网格基础上，本发明技术能比传统的使用线性元素的有限元方法更精确地估计电极上的电压，进而提高了计算电阻抗成像正问题的精度，以及电阻抗成像的精度。
附图说明
下面将通过参照附图详细描述本发明的示例性实施例，使本领域的普通技术人员更清楚本发明的上述及其他特征和优点，附图中：
图1是本发明实施例提供的基于电流量的电阻抗成像方法的流程示意图；
图2是本发明实施例提供的网格中的元素为三角形的示意图；
图3是本发明实施例提供的网格中的元素为四面体的示意图；
图4是本发明实施例提供的一个三角形元素中电流密度的方向示意图；
图5是本发明实施例提供的围绕三角形网格中某个内部节点的有向闭合折线段的示意图；
图6是本发明实施例提供的二维电阻抗模型图；
图7是本发明实施例提供的成像区域上的二维网格与非均匀分布的电导率的示意图；
图8是传统方法得到的成像区域中的电流密度分布图；
图9是本发明实施例提供的成像区域中的电流密度分布图；
图10-图13是本发明实施例提供的在不同视角下的三维电阻抗模型图；
图14是本发明实施例提供的成像区域上非均匀分布的电导率示意图
图15是传统有限元方法得到的圆柱下部区域中电流密度分布示意图；
图16是本发明技术方案得到的圆柱下部区域中电流密度分布示意图；
图17是传统有限元方法得到的圆柱中部区域中电流密度分布示意图；
图18是本发明技术方案得到的圆柱中部区域中电流密度分布示意图；
图19是传统有限元方法得到的圆柱上部区域中电流密度分布示意图；
图20是本发明技术方案得到的圆柱上部区域中电流密度分布示意图；
图21是本发明实施例提供的基于电流量的电阻抗成像装置的结构框图。
具体实施方式
下面结合附图并通过具体实施方式来进一步说明本发明的技术方案。可以理解的是，此处所描述的具体实施例仅仅用于解释本发明，而非对本发明的限定。另外还需要说明的是，为了便于描述，附图中仅示出了与本发明相关的部分而非全部结构。
图1是本发明实施例提供的基于电流量的电阻抗成像方法的流程示意图。该方法可以由基于电流量的电阻抗成像装置来执行。基于电流量的电阻抗成像装置可以由软件和/或硬件的方式来实现。如图1所示，该方法包括：
步骤110、对成像区域进行网格化，并获取网格参数。
其中，网格参数包括网格节点信息、构成网格的元素的元素信息、元素间的邻近关系信息和成像区域的边界信息。
具体地，可采用网格化程序对成像区域进行网格化，每一个小网格即为一个元素，成像区域可以为二维区域，也可以为三维区域；当成像区域为二维区域时，每个元素均为三角形，当成像区域可以为三维区域时，每个元素均为四面体。
本实施例中，网格节点信息包括节点矩阵，在节点矩阵中，第i行是第i个节点的笛卡尔坐标,i＝1,2,…,N，其中N是网格中的总节点数；元素信息包括元素矩阵和元素的电导率常数(电导率常数在一个元素中处处相同)，其中，元素矩阵的第t行是第t个元素所有顶点的序号。例如，若第三个元素由网格中第1,4,5,8个节点构成，则元素矩阵的第三行是(1,4,5,8)，每个元素的顶点按图2或图3中的相对位置排序。元素间的邻近关系信息包括描述元素间邻近关系的邻近关系矩阵；对于由三角形组成的二维网格，如果第j个元素的第i个边也是第k个元素的某个边，则邻近关系矩阵中位于第i 行、第j列的数是k，其中元素的第i个边是与其第i个顶点相对的边。对于由四面体组成的三维网格，如果第j个元素的第i个面也是第k个元素的某个面，则邻近关系矩阵中位于第i行、第j列的数是k，其中元素的第i个面是与其第i个顶点相对的面。另外，成像区域的边界信息包括一组描述成像区域边界信息的向量；对于由三角形组成的二维网格，第一个向量给出了所有在成像区域边界上但不属于任何一个电极的边的序号。第l个向量给出了所有属于第l-1个电极的边的序号，l＝2,…,L+1，其中L表示总的电极数，其中网格中的第(i-1)·3+j条边是其第i个元素的第j条边。例如，若第一个电极是由网格中第1,3,6,9条边构成，则第2个边界向量是(1,3,6,9)。对于由四面体组成的三维网格，第一个向量给出了所有在成像区域边界上但不属于任何一个电极的面的序号。第l个向量给出了所有属于第l-1个电极的面的序号,l＝2,…,L+1，其中网格中的第(i-1)·4+j个面是其第i个元素的第j个面。
步骤120、基于网格参数，确定成像区域的边界电流量条件以及电导率分布的估计。
步骤130、在电流密度在每个元素中恒定不变且处处相同的条件下，基于网格参数、边界电流量条件、成像区域中电流密度分布的无源性、电场分布的无旋性以及成像区域中的电导率分布的估计，确定流过每个元素边界的电流量。
示例性地，基于网格参数、成像区域内电流密度分布的无源性、成像区域边界电流量条件以及成像区域的电场分布的无旋性，确定当成像区域中电导率服从某个分布时，每个元素边界的电流量，可包括：
A、基于网格参数、成像区域中电流密度分布的无源性、电场分布的无旋性，以及成像区域边界电流量条件，建立关于电流量的第一线性方程组。
B、基于稀疏QR分解的直接法求第一线性方程组的最小二乘解，得到每个元素边界的电流量。
具体地，由于成像区域中电流密度分布是无源的，因此，对于三角形二维网格，通过元素边界的电流量满足以下关系：
It,1+It,2+It,3＝0 (2.1.1)
对于四面体三维网格，通过元素边界的电流量满足以下关系：
It,1+It,2+It,3+It,4＝0 (2.1.2)
元素边界的电流量可以用It,v表示，It,v表示通过第t个元素的第v条边流出该元素的电流量，t＝1,2,…T，T为元素的个数，v＝1,2,3或v＝1,2,3,4。
由于电流不通过成像区域边界上电极以外的区域，因此，对于三角形二维网格，通过在网格的边界但不属于任何一个电极的元素边界的电流量满足以下关系：
Figure GDA0003782851430000071
Figure GDA0003782851430000072
Figure GDA0003782851430000073
表示在网格的边界但不属于任何一个电极的边的条数，其中的第b条边是第
Figure GDA0003782851430000074
个元素的第
Figure GDA0003782851430000075
条边。
对于四面体三维网格，通过在网格的边界但不属于任何一个电极的元素边界的电流量满足以下关系：
Figure GDA0003782851430000076
Figure GDA0003782851430000077
Figure GDA0003782851430000078
表示在网格的边界但不属于任何一个电极的面的个数，其中的第b个面是第
Figure GDA0003782851430000079
个元素的第
Figure GDA00037828514300000710
个面。
在电极附近，电流密度应垂直于电极所在的每个面以保证每个电极上的电压处处相同。因此，对于三角形二维网格，元素中的电流密度具有如下关系：
Figure GDA00037828514300000711
Figure GDA00037828514300000712
L表示总的电极数，Sl表示第l个电极的边的条数，第l个电极的第s条边是第tl,s个元素的第vl,s条边。例如，参考图4，电极的某条边是第t个三角形元素的第三条边，且电流通过该条边所属的电极流出成像区域，则电流密度垂直于该三角形元素的第三条边，且
Figure GDA00037828514300000713
Figure GDA00037828514300000714
对于四面体三维网格，元素中的电流密度具有如下关系：
Figure GDA0003782851430000081
Figure GDA0003782851430000082
L表示总的电极数，Sl表示第l个电极的面的个数，第l个电极的第s个面是第tl,s个元素的第vl,s个面。例如，若构成电极的某个面是第t个四面体的第三个面，则
Figure GDA0003782851430000083
Figure GDA0003782851430000084
Figure GDA0003782851430000085
其中θk3表示元素第k个面与第三个面之间的夹角,k＝1,2,4。
本实施例中，通过每个电极的电流量是已知的：
Figure GDA0003782851430000086
Il表示通过第l个电极流出成像区域的电流量。需要说明的是，当Il<0时，电流通过该电极流入成像区域。
由于电流在成像区域内是连续流动的，因此，对于二维三角形网格，通过每个元素边界的电流量还具有以下关系：
Figure GDA0003782851430000087
Figure GDA0003782851430000088
EI表示网格内边的个数，第i条边既是第
Figure GDA0003782851430000089
个元素的第
Figure GDA00037828514300000810
条边又是第
Figure GDA00037828514300000811
个元素的第
Figure GDA00037828514300000812
条边。
对于三维四面体网格，通过每个元素边界的电流量还具有以下关系：
Figure GDA00037828514300000813
Figure GDA00037828514300000814
FI表示网格内面的个数，第i个面既是第
Figure GDA00037828514300000815
个元素的第
Figure GDA00037828514300000816
个面又是第
Figure GDA0003782851430000091
个元素的第
Figure GDA0003782851430000092
个面。
重要的是，由于成像区域中电场分布是无旋的，因此，对于二维三角形网格，电场强度沿着围绕每个内部节点的一个闭合折线段的线积分为零，例如，参考图5，闭合折线段由包含该内部节点的每个三角形的一条边组成，闭合折线段中的每条边都是该内部节点在包含它的每个三角形中的一条对立边。电场强度的具体关系如下：
Figure GDA0003782851430000093
其中，
Figure GDA0003782851430000094
Figure GDA0003782851430000095
Figure GDA0003782851430000096
k＝1,2,…,NI，NI表示除边界外的网格内部的节点的个数，第k个节点属Nk个不同的元素，包含第k个节点的第n个元素是网格中第
Figure GDA0003782851430000097
个元素，在第
Figure GDA0003782851430000098
个元素中，这第 k个节点是其第
Figure GDA0003782851430000099
个顶点，
Figure GDA00037828514300000910
表示第
Figure GDA00037828514300000911
个元素中处处相同的电场强度；
Figure GDA00037828514300000912
表示第t个元素中处处相同的电场强度，vt,v是和第t个三角形中的第v个顶点在该三角形中对立的边有关的向量，vt,v的大小是这条边的长度，方向是这条边的外法向方向，σt表示第t个元素中的电导率常数，At表示第t个元素的面积。
上述关系式(2.6.2-4)是由于：
Figure GDA00037828514300000913
Figure GDA00037828514300000914
并且，
Figure GDA00037828514300000915
Figure GDA00037828514300000916
Figure GDA00037828514300000917
其中
Figure GDA0003782851430000101
是vt,v逆时针旋转90度后得到的向量。
对于三维四面体网格，电场强度沿着围绕网格内每条边的一条闭合折线段的线积分为零。围绕网格内第k条边的闭合折线段由包括第k条边的每个四面体中的一条边组成。这个闭合折线中的每条边分别是第k条边在包含它的每个四面体中对立的边。电场强度的具体关系如下：
Figure GDA0003782851430000102
其中，
Figure GDA0003782851430000103
Figure GDA0003782851430000104
Figure GDA0003782851430000105
Figure GDA0003782851430000106
Figure GDA0003782851430000107
Figure GDA0003782851430000108
k＝1,2,…,EI，EI表示除边界外的网格内部的边的条数，第k条边属Ek个不同的元素，第n个元素是网格中第
Figure GDA0003782851430000109
个元素，在第
Figure GDA00037828514300001010
个元素中，第k条边是其第
Figure GDA00037828514300001011
面与第
Figure GDA00037828514300001012
面相交得到的，右手四指顺着从第
Figure GDA00037828514300001013
面到第
Figure GDA00037828514300001014
面的方向弯曲时，右手大拇指指向第 k条边的方向，
Figure GDA00037828514300001015
表示第
Figure GDA00037828514300001016
个元素中处处相同的电场强度，
Figure GDA00037828514300001017
表示第
Figure GDA00037828514300001018
个元素中从其顶点
Figure GDA00037828514300001019
到其顶点
Figure GDA00037828514300001020
的向量；
Figure GDA00037828514300001021
表示第t个元素中处处相同的电场强度， vt,(i,j)表示第t个元素中从其顶点i到其顶点j的向量，
Figure GDA00037828514300001022
表示第t个元素中从第e个顶点到第f个顶点的向量，e,f＝1,2,3,4；e≠f，σt表示第t个元素中的电导率常数，Vt表示第t个元素的体积。
上述关系式(2.6.11-16)是由于(2.6.6)和
It,v＝jt·vt,v (2.6.17)
其中v＝1,2,3,4，并且在第t个四面体中，
Figure GDA0003782851430000111
Figure GDA0003782851430000112
Figure GDA0003782851430000113
Figure GDA0003782851430000114
Figure GDA0003782851430000115
Figure GDA0003782851430000116
示例性地，当成像区域为二维的，且网格元素为三角形时，通过联立上述关系式(2.1.1,2.2.1,2.3.1,2.4.1,2.5.1,2.6.1)，得到第一线性方程组：
Figure GDA0003782851430000117
其中，
Figure GDA0003782851430000121
通过以下关系式中的一个来表示：
Figure GDA0003782851430000122
Figure GDA0003782851430000123
Figure GDA0003782851430000124
当成像区域为三维的，且网格元素为四面体时，通过联立上述关系式(2.1.2,2.2.2,2.3.4,2.4.1,2.5.2,2.6.10)，得到第一线性方程组：
Figure GDA0003782851430000125
其中，
Figure GDA0003782851430000126
通过以下关系式中的一个来表示：
Figure GDA0003782851430000127
Figure GDA0003782851430000128
Figure GDA0003782851430000129
Figure GDA00037828514300001210
Figure GDA00037828514300001211
Figure GDA00037828514300001212
有上述分析可以将第一线性方程组表示为：
KI＝F (2.7.1)
对于三角形元素，I是一个有3·T个分量的列向量，其第3·(t-1)+v个分量是It,v。K的列数是3·T，行数是T+EB+L+EI+NI。对于四面体元素，I是一个有4·T个分量的列向量，其第4·(t-1)+v个分量是It,v。K的列数是4·T，行数是T+FB+L+ EI+FI。由此，K是稀疏的并且其行数大于列数。因此，可基于稀疏QR分解的直接法求第一线性方程组的最小二乘解，得到每个元素边界的电流量。
步骤140、基于每个元素边界的电流量与网格参数，确定成像区域网格中的电流密度分布。
具体地，基于上述步骤130，对于二维三角形网格，基于以下公式确定所述成像区域网格中的电流密度分布
Figure GDA0003782851430000131
其中，jt表示第t个元素中的电流密度，
Figure GDA0003782851430000132
表示第t个元素中从第e个顶点到第f个顶点的向量；
基于各元素的电流密度以及电导率常数确定各元素的电场强度。
对于三维四面体网格，基于以下公式确定所述成像区域网格中的电流密度分布
Figure GDA0003782851430000133
其中，jt表示第t个元素中的电流密度；
基于各元素的电流密度以及电导率常数确定各元素的电场强度。
另外，对于上述两种情况的电流密度公式，每个电流密度公式中，各等式得到的电流密度的值略有不同，因此，对于每一情况的电流密度公式，可将各等式得到的电流密度取平均后作为最终的电流密度jt。
步骤150、基于电流密度分布与电导率分布的估计，确定成像区域网格中的电场强度分布。
对于二维三角形网格，基于以下公式确定成像区域网格中的电场强度分布：
Figure GDA0003782851430000141
其中，
Figure GDA0003782851430000142
表示第t个元素中的电场强度。
对于三维四面体网格，基于以下公式确定成像区域网格中的电场强度分布：
Figure GDA0003782851430000143
其中，
Figure GDA0003782851430000144
表示第t个元素中的电场强度。
步骤160、基于电场强度分布确定成像区域网格中的电势分布以及电极电压。
具体地，基于上述技术方案，基于电场强度分布确定成像区域网格中的电势分布以及电极电压，包括：
基于电场强度分布建立关于网格中每个节点处电压的第二线性方程组；
利用最小二乘法求解第二线性方程组，得到网格中所有节点处的电压。
可选地，第二线性方程组为：
Figure GDA0003782851430000145
其中，ut,e-ut,f表示第t个元素的第e个顶点与第f个顶点处的电势差。
步骤170、反复更新对电导率分布的估计直至在一电导率分布的估计下计算得到的电极电压与实际测量得到的电极电压的相对误差不超过0.1％。
步骤180、基于最终更新的电导率分布的估计进行图像重构。
基于上述技术方案，本发明实施例分别对二维电阻抗模型以及三维电阻抗模型进行了仿真验证。
具体地，针对二维电阻抗模型，参考图6，考虑一个半径为一米的圆形区域的电导率分布。在图6所述成像区域的边界安置两个尺寸相同的电极(电极I和电极II)。每个电极与圆形区域的圆心张成的角都是60度。1安培的电流通过底部的电极I输入圆形区域，然后通过顶部的电极II从圆形区域中输出。
当电导率在该区域中处处是1西门子每米时，电极I处的电压减去电极II处的电压的正确值是
Figure GDA0003782851430000151
这可通过复变函数中的保角变换方法得到。同时，我们使用传统的有限元方法和本发明技术方案分别计算这个两个电极间的电势差。具体地，我们首先使用如参考图7所示的网格对圆形区域进行划分。该网格共有4423个三角形元素。在使用传统的有限元方法时，我们使用线性元素进行近似，然后得到两个电极间的电势差的估计是1.2651 伏特。在使用本发明技术方案时，我们使用电极上靠近其边缘的点处的电压作为电极电压的估计，然后得到电极间电势差的估计是1.2749伏特。显然，本发明技术得到的结果比传统有限元方法得到的结果更加地接近真实值。在同一网格条件下，本发明技术方案能比使用线性元素的传统有限元方法能更精确地估计电极上的电压。
当电导率在圆形区域中按图7所示的方式非均匀分布(黑色、灰色区域的电导率分别代表电导率是0.1、0.3西门子每米，其他区域中的电导率是1西门子每米)时，使用线性元素的传统的有限元方法与本发明技术方案得到的成像区域中的电流密度分布的估计分别如图8和9所示。从图中可以看出，本发明技术方案得到电流密度分布与传统方法得到的结果很接近。同时，我们知道电流密度在每个电极的两个端点附近应是剧烈变化的，而本发明技术方案能比传统有限元方法更精确地模拟这一特点。
另外，针对三维电阻抗模型，参考图10和图11，成像区域是一个近似为圆柱的柱形，其高3米，底面上的所有边界点都在半径是1米的圆上，两个大小相同的电极安置在柱形的侧面，两个电极关于一个垂直于地面的平面对称。两个电极的边界点在其关于对称的平面上的投影都在一个半径是0.75米、圆心的高是1.5米的圆上。1安培的电流从一侧的电极进入区域，然后从另一侧的电极离开区域。参考图12和图13，成像区域上生成四面体网格，网格中共有12849个四面体元素。成像区域中的电导率分布如图14所示。在成像区域的中心处和两个电极的附近，电导率是0.1西门子每米，其他处的电导率分布都是1西门子每米。图15和图16显示了圆柱中高度在0.1米至0.5米之间的元素中电流密度分布的估计，图17和图18显示了圆柱中高度在1.3米至1.7米之间的元素中电流密度分布的估计，图19和图20显示了圆柱中高度在2.5 米与2.9米之间的元素中电流密度分布的估计。由图15-图20可以看出，本发明技术方案与使用线性元素的传统有限元方法得到的对应区域中电流密度分布的估计在大部分处都是非常接近的。和二维情形类似，我们知道电流密度在靠近电极边界处应是剧烈变化的。而本发明技术能比传统有限元方法更精确地模拟这一特点。
综上，本实施例提供的基于电流量的电阻抗成像方法，通过对成像区域进行网格化，在电流密度在每个元素中恒定不变且处处相同的条件下，基于网格参数、成像区域边界的电流量条件、成像区域中电流密度分布的无源性、以及电场分布的无旋性，确定当成像区域中电导率服从某个分布时，流过每个元素边界的电流量；然后基于流过每个元素边界的电流量得到成像区域的电势分布；最后通过比较某个电导率分布下的电极电压与实际测量得到的电极电压进行图像重构。在给定网格条件下，本方法能比使用线性元素的有限元方法更精确地计算电阻抗成像正问题的解，从而提高了电阻抗成像的精度。
另外，本发明实施例还提供了一种基于电流量的电阻抗成像装置。图21是本发明实施例提供的基于电流量的电阻抗成像装置的结构框图，如图21所述，基于电流量的电阻抗成像装置包括：
网格化模块10，用于对成像区域进行网格化，并获取网格参数，其中，网格参数包括网格节点信息、构成网格的元素的元素信息、元素间的邻近关系信息和成像区域的边界信息；
电导率估计确定模块20，用于基于网格参数，确定成像区域的边界电流量条件以及电导率分布的估计；
电流量确定模块30，用于在电流密度在每个元素中恒定不变且处处相同的条件下，基于网格参数、边界电流量条件、成像区域中电流密度分布的无源性、电场分布的无旋性以及成像区域中的电导率分布的估计，确定流出每个元素边界的电流量；
电流密度确定模块40，用于基于每个元素边界的电流量与网格参数，确定成像区域网格中的电流密度分布；
电场强度确定模块50，用于基于电流密度分布与电导率分布的估计，确定成像区域网格中的电场强度分布；
电势确定模块60，用于基于电场强度分布确定成像区域网格中的电势分布以及电极电压；
更新模块70，用于反复更新对电导率分布的估计直至在一电导率分布的估计下计算得到的电极电压与实际测量得到的电极电压的相对误差不超过0.1％；
图像重构模块80，用于基于最终更新的电导率分布的估计进行图像重构。
本实施例提供的基于电流量的电阻抗成像装置，与本发明任意实施例所提供的基于电流量的电阻抗成像方法属于同一发明构思，可执行本发明任意实施例所提供的基于电流量的电阻抗成像方法，具备相应的功能和有益效果。未在本实施例中详尽描述的技术细节，可参见本发明任意实施例提供的基于电流量的电阻抗成像方法。
本发明实施例还提供了一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现如本发明实施例所提供的基于电流量的电阻抗成像方法，该方法包括：
对成像区域进行网格化，并获取网格参数，其中，网格参数包括网格节点信息、构成网格的元素的元素信息、元素间的邻近关系信息和成像区域的边界信息；
基于网格参数，确定成像区域的边界电流量条件以及电导率分布的估计；
在电流密度在每个元素中恒定不变且处处相同的条件下，基于网格参数、边界电流量条件、成像区域中电流密度分布的无源性、电场分布的无旋性以及成像区域中的电导率分布的估计，确定流过每个元素边界的电流量；
基于每个元素边界的电流量与网格参数，确定成像区域网格中的电流密度分布；
基于电流密度分布与电导率分布的估计，确定成像区域网格中的电场强度分布；
基于电场强度分布确定成像区域网格中的电势分布以及电极电压；
反复更新对电导率分布的估计直至在一电导率分布的估计下计算得到的电极电压与实际测量得到的电极电压的相对误差不超过0.1％；
基于最终更新的电导率分布的估计进行图像重构。
本发明实施例的计算机存储介质，可以采用一个或多个计算机可读的介质的任意组合。计算机可读介质可以是计算机可读信号介质或者计算机可读存储介质。计算机可读存储介质例如可以是——但不限于——电、磁、光、电磁、红外线、或半导体的系统、装置或器件，或者任意以上的组合。计算机可读存储介质的更具体的例子(非穷举的列表)包括：具有一个或多个导线的电连接、便携式计算机磁盘、硬盘、随机存取存储器(RAM)、只读存储器(ROM)、可擦式可编程只读存储器(EPROM或闪存)、光纤、便携式紧凑磁盘只读存储器(CD-ROM)、光存储器件、磁存储器件、或者上述的任意合适的组合。在本文件中，计算机可读存储介质可以是任何包含或存储程序的有形介质，该程序可以被指令执行系统、装置或者器件使用或者与其结合使用。
计算机可读的信号介质可以包括在基带中或者作为载波一部分传播的数据信号，其中承载了计算机可读的程序代码。这种传播的数据信号可以采用多种形式，包括但不限于电磁信号、光信号或上述的任意合适的组合。计算机可读的信号介质还可以是计算机可读存储介质以外的任何计算机可读介质，该计算机可读介质可以发送、传播或者传输用于由指令执行系统、装置或者器件使用或者与其结合使用的程序。
计算机可读介质上包含的程序代码可以用任何适当的介质传输，包括——但不限于无线、电线、光缆、RF等等，或者上述的任意合适的组合。
可以以一种或多种程序设计语言或其组合来编写用于执行本发明操作的计算机程序代码，所述程序设计语言包括面向对象的程序设计语言—诸如Java、Smalltalk、 C++，还包括常规的过程式程序设计语言—诸如“C”语言或类似的程序设计语言。程序代码可以完全地在用户计算机上执行、部分地在用户计算机上执行、作为一个独立的软件包执行、部分在用户计算机上部分在远程计算机上执行、或者完全在远程计算机或终端上执行。在涉及远程计算机的情形中，远程计算机可以通过任意种类的网络——包括局域网(LAN)或广域网(WAN)—连接到用户计算机，或者，可以连接到外部计算机(例如利用因特网服务提供商来通过因特网连接)。
注意，上述仅为本发明的较佳实施例及所运用技术原理。本领域技术人员会理解，本发明不限于这里所述的特定实施例，对本领域技术人员来说能够进行各种明显的变化、重新调整、相互结合和替代而不会脱离本发明的保护范围。因此，虽然通过以上实施例对本发明进行了较为详细的说明，但是本发明不仅仅限于以上实施例，在不脱离本发明构思的情况下，还可以包括更多其他等效实施例，而本发明的范围由所附的权利要求范围决定。

Claims (10)

1.一种基于电流量的电阻抗成像方法，其特征在于，包括：
对成像区域进行网格化，并获取网格参数，其中，所述网格参数包括网格节点信息、构成网格的元素的元素信息、元素间的邻近关系信息和成像区域的边界信息；
基于所述网格参数，确定所述成像区域的边界电流量条件以及电导率分布的估计；
在电流密度在每个所述元素中恒定不变且处处相同的条件下，基于所述网格参数、所述边界电流量条件、所述成像区域中电流密度分布的无源性、电场分布的无旋性以及所述成像区域中的电导率分布的估计，确定流过每个所述元素边界的电流量；
基于每个所述元素边界的电流量与所述网格参数，确定所述成像区域网格中的电流密度分布；
基于所述电流密度分布与所述电导率分布的估计，确定所述成像区域网格中的电场强度分布；
基于所述电场强度分布确定所述成像区域网格中的电势分布以及电极电压；
反复更新对电导率分布的估计直至在一电导率分布的估计下计算得到的电极电压与实际测量得到的电极电压的相对误差不超过0.1％；
基于最终更新的电导率分布的估计进行图像重构。
2.根据权利要求1所述的电阻抗成像方法，其特征在于，基于所述网格参数、所述边界电流量条件、所述成像区域中电流密度分布的无源性、电场分布的无旋性以及所述成像区域中的电导率分布，确定流过每个所述元素边界的电流量，包括：
基于所述网格参数、所述边界电流量条件、所述成像区域中电流密度分布的无源性、电场分布的无旋性以及所述成像区域中的电导率分布，建立关于电流量的第一线性方程组；
基于稀疏QR分解的直接法求所述第一线性方程组的最小二乘解，得到流过每个所述元素边界的电流量。
3.根据权利要求2所述的电阻抗成像方法，其特征在于，所述成像区域为二维的，且所述元素为三角形时，所述第一线性方程组为：
Figure FDA0003782851420000021
其中，It,v表示通过第t个元素的第v条边流出该元素的电流量，t＝1,2,…T，T为元素的个数，v＝1,2,3；
Figure FDA0003782851420000022
Figure FDA0003782851420000023
表示在网格的边界但不属于任何一个电极的边的条数，其中的第b条边是第
Figure FDA0003782851420000024
个元素的第
Figure FDA0003782851420000025
条边；l＝1,2,…,L，s＝1,2,…,Sl，L表示总的电极数，Sl表示第l个电极的边的条数，第l个电极的第s条边是第tl,s个元素的第vl,s条边；Il表示通过第l个电极流出成像区域的电流量；
Figure FDA0003782851420000026
Figure FDA0003782851420000027
EI表示网格内边的个数，第i条边既是第
Figure FDA0003782851420000028
个元素的第
Figure FDA0003782851420000029
条边又是第
Figure FDA00037828514200000210
个元素的第
Figure FDA00037828514200000211
条边；k＝1,2,…,NI，NI表示除边界外的网格内部的节点的个数，第k个节点属Nk个不同的元素，包含第k个节点的第n个元素是网格中第
Figure FDA00037828514200000212
个元素，在第
Figure FDA00037828514200000213
个元素中，这第k个节点是其第
Figure FDA00037828514200000214
个顶点，
Figure FDA00037828514200000215
表示第
Figure FDA00037828514200000216
个元素中处处相同的电场强度；其中，
Figure FDA00037828514200000217
通过以下关系式中的一个来表示：
Figure FDA00037828514200000218
Figure FDA00037828514200000219
Figure FDA00037828514200000220
其中，
Figure FDA0003782851420000031
表示第t个元素中处处相同的电场强度，vt,v是和第t个三角形中的第v个顶点在该三角形中对立的边有关的向量，vt,v的大小是这条边的长度，方向是这条边的外法向方向，σt表示第t个元素中的电导率常数，At表示第t个元素的面积。
4.根据权利要求3所述的电阻抗成像方法，其特征在于，基于以下公式确定所述成像区域网格中的电流密度分布：
Figure FDA0003782851420000032
其中，jt表示第t个元素中的电流密度，
Figure FDA0003782851420000033
表示第t个元素中从第e个顶点到第f个顶点的向量；
基于以下公式确定成像区域网格中的电场强度分布：
Figure FDA0003782851420000034
其中，
Figure FDA0003782851420000035
表示第t个元素中的电场强度。
5.根据权利要求2所述的电阻抗成像方法，其特征在于，所述成像区域为三维的，且所述元素为四面体时，所述第一线性方程组为：
Figure FDA0003782851420000036
其中，It,v表示通过第t个元素的第v个面流出该元素的电流量，t＝1,2,…T，T为元素的个数，v＝1,2,3,4；
Figure FDA0003782851420000041
Figure FDA0003782851420000042
表示在网格的边界但不属于任何一个电极的面的个数，其中的第b个面是第
Figure FDA0003782851420000043
个元素的第
Figure FDA0003782851420000044
个面；l＝1,2,…,L，s＝1,2,…,Sl，L表示总的电极数，Sl表示第l个电极中面的个数，第l个电极的第s个面是第tl,s个元素的第vl,s个面；Il表示通过第l个电极流出成像区域的电流量；
Figure FDA0003782851420000045
Figure FDA0003782851420000046
FI表示网格内面的个数，第i个面既是第
Figure FDA0003782851420000047
个元素的第
Figure FDA0003782851420000048
个面又是第
Figure FDA0003782851420000049
个元素的第
Figure FDA00037828514200000410
个面；k＝1,2,…,EI，EI表示除边界外的网格内部的边的条数，第k条边属Ek个不同的元素，第n个元素是网格中第
Figure FDA00037828514200000411
个元素，在第
Figure FDA00037828514200000412
个元素中，第k条边是其第
Figure FDA00037828514200000413
面与第
Figure FDA00037828514200000414
面相交得到的，右手四指顺着从第
Figure FDA00037828514200000415
面到第
Figure FDA00037828514200000416
面的方向弯曲时，右手大拇指指向第k条边的方向，
Figure FDA00037828514200000417
表示第
Figure FDA00037828514200000418
个元素中处处相同的电场强度，
Figure FDA00037828514200000419
表示第
Figure FDA00037828514200000420
个元素中从其顶点
Figure FDA00037828514200000421
到其顶点
Figure FDA00037828514200000422
的向量；其中，
Figure FDA00037828514200000423
通过以下关系式中的一个来表示：
Figure FDA00037828514200000424
Figure FDA00037828514200000425
Figure FDA00037828514200000426
Figure FDA00037828514200000427
Figure FDA00037828514200000428
Figure FDA00037828514200000429
其中，
Figure FDA00037828514200000430
表示第t个元素中处处相同的电场强度，
Figure FDA00037828514200000431
表示第t个元素中从第e个顶点到第f个顶点的向量，e,f＝1,2,3,4；e≠f，σt表示第t个元素中的电导率常数，Vt表示第t个元素的体积。
6.根据权利要求5所述的电阻抗成像方法，其特征在于，基于以下公式确定所述成像区域的电流密度分布：
Figure FDA0003782851420000051
其中，jt表示第t个元素中的电流密度；
基于以下公式确定成像区域网格中的电场强度分布：
Figure FDA0003782851420000052
其中，
Figure FDA0003782851420000053
表示第t个元素中的电场强度。
7.根据权利要求4或6所述的电阻抗成像方法，其特征在于，基于所述电场强度分布确定所述成像区域网格中的电势分布以及电极电压，包括：
基于所述电场强度分布建立关于网格中每个节点处电压的第二线性方程组；
利用最小二乘法求解所述第二线性方程组，得到网格中所有节点处的电压。
8.根据权利要求7所述的电阻抗成像方法，其特征在于，所述第二线性方程组为：
Figure FDA0003782851420000054
其中，ut,e-ut,f表示第t个元素的第e个顶点与第f个顶点处的电势差。
9.一种基于电流量的电阻抗成像装置，其特征在于，包括：
网格化模块，用于对成像区域进行网格化，并获取网格参数，其中，所述网格参数包括网格节点信息、构成网格的元素的元素信息、元素间的邻近关系信息和成像区域的边界信息；
电导率估计确定模块，用于基于所述网格参数，确定所述成像区域的边界电流量条件以及电导率分布的估计；
电流量确定模块，用于在电流密度在每个所述元素中恒定不变且处处相同的条件下，基于所述网格参数、所述边界电流量条件、所述成像区域中电流密度分布的无源性、电场分布的无旋性以及所述成像区域中的电导率分布的估计，确定流出每个所述元素边界的电流量；
电流密度确定模块，用于基于每个所述元素边界的电流量与所述网格参数，确定所述成像区域网格中的电流密度分布；
电场强度确定模块，用于基于所述电流密度分布与所述电导率分布的估计，确定所述成像区域网格中的电场强度分布；
电势确定模块，用于基于所述电场强度分布确定所述成像区域网格中的电势分布以及电极电压；
更新模块，用于反复更新对电导率分布的估计直至在一电导率分布的估计下计算得到的电极电压与实际测量得到的电极电压的相对误差不超过0.1％；
图像重构模块，用于基于最终更新的电导率分布的估计进行图像重构。
10.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，该程序被处理器执行时实现如权利要求1-8中任一所述的基于电流量的电阻抗成像方法。
CN202010079013.8A 2020-02-03 2020-02-03 基于电流量的电阻抗成像方法、装置及存储介质 Active CN113273988B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010079013.8A CN113273988B (zh) 2020-02-03 2020-02-03 基于电流量的电阻抗成像方法、装置及存储介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010079013.8A CN113273988B (zh) 2020-02-03 2020-02-03 基于电流量的电阻抗成像方法、装置及存储介质

Publications (2)

Publication Number Publication Date
CN113273988A CN113273988A (zh) 2021-08-20
CN113273988B true CN113273988B (zh) 2022-10-11

Family

ID=77274840

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010079013.8A Active CN113273988B (zh) 2020-02-03 2020-02-03 基于电流量的电阻抗成像方法、装置及存储介质

Country Status (1)

Country Link
CN (1) CN113273988B (zh)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101794453A (zh) * 2010-03-25 2010-08-04 河北工业大学 基于回归分析的节点映射图像重构方法
CN102599907A (zh) * 2012-04-09 2012-07-25 成都晨德科技有限公司 基于网格位移模型的电阻抗断层成像方法
JP2015160004A (ja) * 2014-02-27 2015-09-07 日本光電工業株式会社 電気的インピーダンス測定装置

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101794453A (zh) * 2010-03-25 2010-08-04 河北工业大学 基于回归分析的节点映射图像重构方法
CN102599907A (zh) * 2012-04-09 2012-07-25 成都晨德科技有限公司 基于网格位移模型的电阻抗断层成像方法
JP2015160004A (ja) * 2014-02-27 2015-09-07 日本光電工業株式会社 電気的インピーダンス測定装置

Also Published As

Publication number Publication date
CN113273988A (zh) 2021-08-20

Similar Documents

Publication Publication Date Title
Chen et al. On choosing the location of the sources in the MFS
Christen et al. Markov chain Monte Carlo using an approximation
Ren et al. Inclusion boundary reconstruction and sensitivity analysis in electrical impedance tomography
Liu et al. Eses: Software for e ulerian solvent excluded surface
Fox et al. Sampling conductivity images via MCMC
Duraiswami et al. Boundary element techniques for efficient 2-D and 3-D electrical impedance tomography
Yao et al. Enhanced supervised descent learning technique for electromagnetic inverse scattering problems by the deep convolutional neural networks
Li et al. A three-dimensional model-based inversion algorithm using radial basis functions for microwave data
Achitouv et al. Improving reconstruction of the baryon acoustic peak: The effect of local environment
Pollok et al. Magnetic field prediction using generative adversarial networks
Cao et al. Direct image reconstruction for 3-D electrical resistance tomography by using the factorization method and electrodes on a single plane
Kriváchy et al. High-speed batch processing of semidefinite programs with feedforward neural networks
Sritharan et al. Computing the Riemannian curvature of image patch and single-cell RNA sequencing data manifolds using extrinsic differential geometry
Mücke et al. Markov chain generative adversarial neural networks for solving Bayesian inverse problems in physics applications
Dimas et al. Advances in Electrical Impedance Tomography Inverse Problem Solution Methods: From Traditional Regularization to Deep Learning
CN113273988B (zh) 基于电流量的电阻抗成像方法、装置及存储介质
Elhag et al. Manifold diffusion fields
Hon et al. A cell based particle method for modeling dynamic interfaces
Rashid et al. A dynamic oppositional biogeography-based optimization approach for time-varying electrical impedance tomography
Arif et al. State estimation approach to dual-modal imaging of two-phase flow based on electromagnetic flow tomography and electrical tomography
Carpio et al. Hybrid Topological Derivative‐Gradient Based Methods for Nondestructive Testing
Hoffmann et al. Robustness of topological defects in discrete domains
Zhang et al. HybridDenseU-Net: learning a multi-scale convolution and dense connectivity CNN for inverse imaging problems
CN113204905A (zh) 一种接触式激发极化法有限单元数值模拟方法
Rapetti Comments on a high-order Whitney complex for simplices

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant