一种快速微波成像方法与流程

本发明属于微波成像技术领域，提出了一种快速微波成像方法。

背景技术：

微波成像技术以微波在各种复杂媒质中的传播和散射的研究为基础，通过测量被测媒质外部的散射场数据，重建被测媒质内部的复介电常数图像。被测的散射场携带大量有关散射体的信息，利用关于散射目标的先验知识，经过适当的数学处理之后可以提取出散射体本身所具有的某些特性，如散射体的形状，介电常数的分布等。由于微波成像能实现对目标的无损检测，能同时对目标进行几何成像和物理成像，且具有较高分辨率，因此微波成像技术在遥感，医学成像，雷达目标识别等领域有广泛的用途。

微波成像技术一般包含两部分：一部分是正问题计算，即计算给定模型的电场分布；另一部分是逆问题的求解，即根据给定的测量场重构电场的分布。当未知目标的物理尺寸比拟，甚至大于入射波的波长，逆问题就变得非线性，病态，并且计算量巨大。这些困难在处理三维微波成像问题中变得尤其棘手。根据等效电流是否包含在内，反演方法可大致分为两类：场型反演方法和源型反演方法。场型反演方法主要包括变型波恩迭代法(dbim)和牛顿迭代法等。在dbim方法中，每次需要更新非均匀背景下的格林函数，并且采用了波恩近似用入射场来代替总场。源型反演方法主要包括对比源方法(csi)和子空间优化方法(som)等。在som方法中，等效感应电流被划分为确定性电流和非确定性电流。在som方法中，通过对确定性电流进行一定的处理，提出了一种基于子空间的dbim方法(s-dbim)来加速迭代收敛。在s-dbim的迭代过程中，需要更新背景非均匀格林函数和总场，此时总场不再由入射场代替，而是入射场和散射场的总和。对总场赋予更精确的物理含义，这使得s-dbim收敛更快。但是在s-dbim中，非均匀格林函数和总场都需要更新，s-dbim在处理三维反演问题时计算量变得很大。

自从周永祖教授首次提出dbim反演算法，人们提出了各种改进的dbim算法来解决不同的成像问题，比如乳腺成像、密度成像。由于dbim算法反演速度不够快，我们将其运用到实际问题上还是存在诸多的困难，特别是在三维的微波成像问题上。在dbim算法中，每次迭代我们都需要计算一次正问题。反演的速度在很大程度上取决于正问题的计算速度。因此人们将不同的高效正问题求解器集成到dbim反演算法中来加快正问题的求解速度。稳定双共轭梯度快速傅立叶变换(bcgs-fft)曾作为正问题求解器来反演多层介质问题。但是在这种方法中，均匀离散化的必要性使得这种方法处理复杂几何物体变得很吃力。在2009年，国外学者提出将多层快速多级子(mlfma)和dbim相结合来。然而mlfma这种高效正问题求解器在处理多层介质问题的反演时，变得很棘手。另一方面，之后报导的一种叫作多层格林函数插值方法(mlgfim)被用来解决各种电磁正问题。这种方法采用了bcgs-fft的插值思想和mlfma的多层结构。因此这个算法作为一个可行性高的正问题求解器mlgfim-slover是一个很好的选择。

技术实现要素：

本发明所要解决的技术问题是传统的三维微波成像方法反演速度慢的问题，提出了一种快速微波成像方法，这种方法将修正的ms-dbim,基于变形波恩近似(dbim)反演算法的基础上，集成了一个高效的正问题求解器，该求解器采用了一种与积分核无关的插值算法，即多层格林函数插值算法(mlgfim)，可以将这个求解器命名为mlgfim-solver.其优点是在dbim算法的每一次迭代过程中，通过调用这个高效的正问题求解器，我们能相对快速地计算正问题来提高反演速度。同时为了使整个反演算法稳定地进行迭代反演并且，采用了截断奇异值分解(tsvd)技术，这避免了大量的数值测试来确定正则化项。为了进一步的减少计算量，在反演算法的成本函数中，用入射场来代替总场。因此该发明能快速，稳定地进行反演，保证了图像的反演速度和质量。

本发明的技术方案：

本发明设计方法利用发射天线照射电磁波到探测区域，用接收天线获取到一定量的散射场数据后，在电磁场积分方程的基础上通过离散化探测区域构建相对应的矩阵方程和成本函数方程，用发明的ms-dbim-mlgfim方法来快速迭代求解成本函数的最小值，同时为了使整个算法在迭代过程能稳定地进行反演，我们采用了截断奇异值分解技术来代替繁琐的正则化参数选取，具体过程如下:

本发明快速微波成像方法，包括以下步骤：

步骤(1)、探测区域外设置ninc个发射天线和nr个接收天线，发射天线发射频率f的电磁波到有未知散射体的探测区域，接收天线接收对应的测量散射场

步骤(2)、首先将探测区域离散为m个足够小的离散方块，并根据已知的电磁波频率f、探测区域位置、发射天线和接收天线的位置和极化方式、测量散射场计算格林函数矩阵和离散后入射场并构建mlgfim-solve求解器。

格林函数矩阵为离散后的格林函数的积分算子；为离散后的入射场矩阵其中并矢格林函数为并矢格林函数,表示探测区域内点r′处的点源对其接收天线rs处的作用，是一个3×3的矩阵；标量格林函数g(rs,r′)＝exp(ik0|rs-r′|)/(4π|rs-r′|)，k0为背景介质的波数。是一个3×3的单位对角矩阵，表示旋度算子，i表示复数的虚部。

mlgfim-solve求解器是对(矩阵向量相乘操作进行加速；其中代表任意一个向量)为离散化的格林函数的积分算子；并矢格林函数表示的是探测区域内点r′对探测区域内另一点r的作用；

是一个3m×3m的矩阵，可由9个m×m的矩阵组成,即

当m≠n,(n,m＝1,2,…,m)

当m＝n，

其中w分别是电磁波的角频率；μ0是背景介质的磁导率；ρ(u-v)表示冲激函数，当u＝v时ρ(u-v)＝1，当u≠v时ρ(u-v)＝0；rn,m＝|rn-rm|，rn＝(r1:n,r2；n,r3；n)，rm＝(r1:m,r2；m,r3；m)分别是第n,m个离散方块的中心坐标,其中1，2，3分别代表三维坐标下的x,y,z方向上的分量。

当未知数m很大时,会变成非常大的密集矩阵，如果通过传统的矩量法去填充的每个元素来完成矩阵向量乘法操作计算量很大并且耗费的时间也比较长，因此我们通过构建mlgfim-solver求解器来完成操作。

步骤(3)、定义新的成本函数

在原始的dbim算法中，成本函数可以定义为：

其中是二次入射电场；是上一迭代后的最新离散后的格林函数的积分算子；δ是散射张量系数,为散射强度张量，为一个3m×3m的对角矩阵,可以表示为

当m≠m,2m,3m，则q＝mod(m,m)，mod表示的是取余符号；其他情况下q＝m；vq,∈r；q分别是第q个离散方块的体积和介电常数。γ为tikhonov正则化参数，主要的作用是为了保证反演算法在迭代过程中的稳定性，但是这个参数需要大量的数值测试来确定。

在本发明的ms-dbim-mlgfim算法中，采用了截断奇异值分解(tsvd)技术,可以代替正则化参数的作用。对进行奇异值分解：

其中表示对奇异值分解得到的左边矩阵的第j列，表示第j个电流基；表示的共轭；将奇异值σj的绝对值按照从大到小的顺序排，取前l个较大奇异值σj对应的电流基构成一个电流子空间，这个电流子空间构成了电流的主要组成部分。即l是一个临界值，它满足条件|σj|≥阈值，j≤l。上述电流子空间即为感应电流的确定性部分，可以表示为：

其中表示的共轭。利用tsvd技术，将归一化处理后的探测区域内的测量散射场和感应电流确定性部分的叠加后得到新的成本函数：

其中为每次更新得到的第p个入射场对应的感应电流确定性部分。

步骤(4)、初始化：设定迭代次数

步骤(5)、利用散射强度张量通过公式(8)-(9)分别更新二次入射电场和非均匀格林函数将步骤(2)mlgfim-solver嵌入到共轭梯度迭代法(cg)中来求解

其中是3m×3m的单位对角矩阵；

步骤(6)、对更新后的格林函数矩阵进行奇异值分解

则感应电流相对应的确定性部分可以通过下式得到

步骤(7)、然后将代入到公式(7)得到当前迭代的成本函数f,若f小于阈值，则停止迭代,跳转至步骤(9)，否则跳转至步骤8。

步骤(8)、重置n＝n+1通过求公式(7)的最小值来更新根据公式(12)得到新的散射强度张量；然后跳转至步骤5继续迭代：

步骤(9)、根据当前更新得到的散射强度张量数据，在计算机上编程可视化进行成像。

本发明的有益效果：在三维微波成像中，提出了一种快速成像方法，在dbim反演迭代过程中，每次计算正问题通过调用mlgfim-solver这个高效的正问题求解器。并且在反演算法ms-dbim中，采用了截断奇异值分解技术，避免了繁琐的正则化参数选取过程,为了进一步的减少计算量，在反演成本函数中，用入射场来代替总场,避免了散射场那部分的计算量。本发明所提出的成像方法有很高的反演速度，并且构造出的图像质量高。

附图说明

图1是本发明方法流程图。

图2是三维微波成像的实验测量装置图。

图3是mlgfim-solver中二维插值示意图。

图4是实施例1的实物图。

图5是ms-dbim-mlgfim对实施例1的反演结果，其中(a1)、(a2)分别为z＝0截面的实部、虚部成像图；(b1)、(b2)分别为x＝0截面的实部、虚部成像图；(c1)、(c2)分别为y＝-0.5λ截面的实部、虚部成像图。

图6是实施例2的实物图。

图7是ms-dbim-mlgfim对仿真例2的反演结果；其中(a1)、(a2)分别为x＝0截面的实部、虚部成像图；(b1)、(b2)分别为z＝-0.4λ截面的实部、虚部成像图；(c1)、(c2)分别为y＝-0.5λ截面的实部、虚部成像图。

具体实施方式

下面结合附图和仿真测试例子对本技术做进一步详细说明。

首先我们对三维电磁场正问题进行描述，如图2所示。在三维坐标系下，一个在长方体探测区域内d内的非磁性散射体被ninc个平面波照射。

假设三维坐标系下任何一点可以表示为r＝(r1,r2,r3)。对于每一个入射天线而言，有nr个接收天线位于r′n＝(r′1:l,r′2；l,r′3；l),l＝1,2,...,nr来测量三个方向上的矢量场。在r处的散射场满足下列积分方程：

其中为对比度函数，∈(r′)表示离散方块r′的介电常数，∈0表示空气中的介电常数，为三维并矢格林函数，表示一个探测区域r′处的点源对另一点r处的作用，是一个3×3的矩阵；标量格林函数g(r,r′)＝exp(ik0|r-r′|)/(4π|rs-r′|)，k0为电磁波的波数；e^tot(r′)表示了位于r′的总场。e^tot＝e^inc+e^sca代表了总场是入射场和散射场之和，所以可得下列方程：

考虑到计算的简便性和mlgfim的应用，可将探测区域离散为m个足够小的小立方体，这些小立方体的中心坐标分别为rm＝(r1；m,r2；m,r3；m)，m＝1,2,3,...,m。因此(14)可以转换成一个矩阵方程：

其中由组成，t表示对一个矩阵进行转置，这表示入射场在m个小立方体上k方向上的场，1,2,3分别代表三维坐标下的x,y,z方向上的分量。和有着相同的结构，式(15)中其他变量在发明步骤中已经详细阐明，这里就不做说明。

在探测区域外的nr个接收天线上的散射场可以通过下式得到：

其中由组成，这表示散射场在nr个接收天线上k方向上的场。

当未知数m很大时，将会变成一个很大的密集矩阵，所以我们构建了mlgfim-solver求解器来降低这个复杂度。mlgfim-solve采用了分层的思想，首先我们将正方体探测区域d作为mlgfim-solver第一层的最大立方体，然后把它划分为8个更小的立方体，每个子立方体又继续往下划分，直到划分后的最低层立方体满足条件。

mlgfim-solver用插值方法来快速地估计矩阵中的元素，如图3所示，立方体m有一源点r′，在立方体n中有一场点r，如果立方体m和立方体n的中心点距离大于其2倍的边长，白色的空心点为插值点。我们可以通过插值的方法来计算这两点的作用gd；uv(r,r′)，即：

其中wm,p，wn,q分别为场立方体m的第p个插值基函数和源立方体n的第q个插值基函数，rm,p,rn,q分别为场立方体m的第p个插值点和源立方体n的第q个插值点，k为插值点数。

通过(17)，场立方体m中所有点和源立方体n中所有点之间的作用可由矩阵表示，其具体表达式为：

其中是插值函数矩阵，由wm,p(r)组成；是作用函数矩阵，由gd；uv(rm,p,rn,q)组成。

假设立方体m中有m个点，立方体n中有n个点，那么用矩量法(mom)直接计算的时间复杂度为o(mn)，用插值公式计算的时间复杂度为o(mk+kk+kn)，其中k为插值点的个数，插值点的数目是有限的，而且很少的插值点就可以达到很高的精度，即k<<min(m,n),所以时间复杂度又可以写成o(mk+kn)。可见mlgfim-solver比mom效率要好得多。

具体的迭代流程如图1所示，本发明设计的一种快速微波成像方法具体实施方法包括以下步骤：

步骤1：探测区域外设置ninc个发射天线和nr个接收天线，发射天线发射频率f的电磁波到有未知散射体的探测区域，接收天线接收对应的测量散射场

步骤2：根据公式(1)-(2)计算格林函数矩阵入射场矩阵由平面波的数学表达式得到，根据公式(17)-(18)来构建mlgfim-solver正问题求解器。

步骤3：根据截断奇异值分解技术定义新的成本函数(7)，并且进行初始化：设定迭代次数n＝0，并且利用波恩近似得到

步骤4：更新二次入射电场和非均匀格林函数在共轭梯度迭代法中cg中嵌入构建好的mlgfim-solver求解器来求解方程(8)-(9)。

步骤5：对更新后的格林函数矩阵进行奇异值分解，然后通过公式(11)得到确定性电流的主导部分

步骤6：计算成本函数(7)的值，判断是否小于阈值，如果满足要求则停止迭代，如果不满足要求，通过求(7)的最小值更新通过(12)得到新的散射强度张量然后转到步骤4继续执行。

本发明可用于快速安检成像，将以往的x光成像换成微波成像，降低了用途成本，并且成像速度快，实验装置结构图如图2所示，实施例的探测区域d为一个边长为2λ，(λ为在自由空间介质中的入射波波长)中心为原点的立方体。用48个发射天线去照射探测区域，这48个发射天线分别位于xoy,yoz,xoz的三个平面圆上，每个圆半径为3λ并且均匀地分布16个发射天线。在xoy,yoz,xoz平面上的天线极化方式分别为z极化,x极化,y极化。接收天线的位置和个数和入射天线是一样的。为了得到实验仿真数据，对探测区域进行剖分为60×60×60个网格，我们采用偶极耦合法(cdm)正向求解器来得到散射场数据。

实施例1

如图4，实施例1是两个半径为0.5λ的球，两个球的相对介电常数为1.5，它们的圆心坐标分别为(0,0.5λ,0)，(0,0.5λ,0)。在ms-dbim-mlgfim的迭代过程中，将成像区域离散为32×32×32个网格，通过18次的迭代，图5为所发明方法的成像结果，图5中的(a),(b),(c)分别为反演结果的z＝0,x＝0,y＝-0.5λ的截面图。可以看出，所发明的反演算法ms-dbim-mlgfim能精确的反演出两个球的位置和尺寸，并且能估计它们的相对介电常数大约为1.5。

仿真例2

为了进一步验证所提出方法的有效性，用一个更复杂的设计目标来测试该方法。这个未知由两个球体和u形结构组成，如图6所示。两个球体的球心分别位于(0,-0.5λ,0.4λ),(0,0.5λ,0.4λ)。u型结构的最低点位于(0,0,-0.9λ)，整个结构关于x,y两个坐标轴对称。将成像区域均匀划分为48×48×48均匀小立方体，每个小立方体的边长为0.03125m，图6物体中的各参数r＝0.35λ,h＝0.5λ,w＝0.4λ,l＝0.1912m,s＝0.2λ。经过16次的迭代，图7为所发明方法的成像结果，图7中的(a),(b),(c)分别为反演结果的x＝0,z＝0.4λ,y＝-0.5λ的截面图。我们可以观察到，u形结构部分的实部介电常数略低于1.5，但散射体的形状和位置仍然得到了很好的重构。

上述两实例仅仅只是例证本发明方法，并非是对于本发明的限制，本发明也并非仅限于上述实例，只要符合本发明方法的要求，均属于本发明方法的保护范围。