CN102854499A - 基于对比源反演算法的三维电磁成像方法 - Google Patents

Abstract

三维电磁成像在无损探伤、医学成像、地表探测、地球物理勘探等实际工程领域中的应用具有巨大的潜力。由于线性成像算法采用了近似处理，重建图像在精确度和准确度方面存在局限性；而非线性算法，如对比源反演(Contrast Source Inverison，CSI)算法，则可以在迭代过程中逐步改善图像，如兴趣目标的位置、形状及材料参数，最后重建与兴趣目标最近似的图像。本发明提出了对多频率、多收发设置的三维电磁实测数据的非线性成像方法一结合了频率跳变方法和正则化方法的三维对比源反演算法。由于无须正演计算，在迭代过程中计算三维并矢格林(Dyadic Green)函数算子及其共轭算子的计算量是可控的，其重建结果验证了该方法在三维电磁成像中应用的精确性和可行性，可广泛应用于要求高精度和高分辨率成像的实际应用领域。

Description

基于对比源反演算法的三维电磁成像方法

技术领域

三维微波成像在无损探伤、医学成像、地表探测、地球物理勘探等实际工程领域中的应用具有巨大的潜力。电磁成像是逆问题，具有非线性和病态特性，难于求解，因此迫切需要精确的成像算法。目前实际应用领域中采用的成熟算法多是线性算法，由于采用了近似处理，其重建图像的分辨率和精确度是非常有限的。而非线性算法则可以在迭代过程中逐步改善图像，如兴趣目标的位置、形状及材料物理参数，最后取得与兴趣目标最近似的图像。正在进行的国家863计划-“绕月探测工程科学数据应用与研究”就明确提出了采用非线性(优化)算法对月球表面探测数据进行处理的研究。本发明基于一种非线性算法-对比源反演(Contrast Source Inversion，CSI)算法对三维电磁数据的成像方法，能够提高重建图像的质量，可广泛应用于工业的无损检测，遥感，医学成像及地球资源探测等领域。

背景技术

成像技术的原理与人的眼睛的成像原理基本类似，由太阳(源)发射出的可见光照射在目标上，一部分光反射(散射)、一部分光透射，眼睛接收到这些反射、透射的可见光，把接收到的信息(数据)传到大脑里，经过数据加工处理后，就形成了图像。人的眼睛只可以接收可见光，可见光的波长范围大约是390nm-780nm。对于自然界存在的大多数固体材料来说，可见光是没法穿透的，也就是说人无法用肉眼观察到目标的内部情况。成像技术就是利用发射设备(天线，传感器等)发射出电磁波，超声波，机械波(弹性波、地震波)等，通过传播媒质(可以是气体，也可以是固体、液体，如空气、水、大地)，到达探测目标后，一部分波反射(散射)、一部分波透射，再利用接收设备接收到这些反射(散射)、透射波后，把接收到的信息(数据)传到计算机里，用成像算法把数据进行加工处理后，就形成了图像。

目前成像技术的研究主要包括发射和接收设备、传播媒质(材料)及成像算法三个部分，本发明是基于对比源反演算法的三维电磁成像方法，该算法是在域积分方程的基础上实现的，因此在本说明书中电磁波的波动方程将用积分形式来表示。如图2所示目标坐落在一个立方体区域D中，假设构成目标的材料是各向同性和非均质的，其材料特性通过与背景媒质的差异来定义，即由对比度值函数χ来表示。共有N个发射器被放置在域D周围的一系列的点上，形成球面区域T，该域半径为r_T，有L个接收器被放置在环形区域M上，该域半径为r_M。在整个说明书中我们使用的坐标系是直角坐标系，其中R表示三维的位置矢量。

电磁成像属于逆散射问题，而逆散射问题具有非线性和病态(ill-posed)特性，目前算法大致可以分为线性和非线性算法。对于线性算法，适当的近似被应用于波动方程中，如Bom近似、Kirchoff近似等。线性算法计算速度快，因此目前所采用的成熟算法多是线性算法，如合成孔径雷达技术(Synetic Aperture Radar，SAR)等。但只是在有限的条件下，这些线性算法才准确，而且从图像中获得的信息也有限，只能获得被测目标的位置和大致形状，而无法获得探测目标的介电常数、电导率和密度等物理参数。非线性算法是在不改变逆散射问题非线性的前提下对目标成像，这些算法具有迭代性质，在迭代过程中逐步的改善图像，最后取得与探测目标最接近的重建图像，从图像中可以获得更多、更准确的信息。Chew提出了能够在实际工程领域应用的非线性算法应当满足的基本条件：1.重建图像的分辨率和能见度要高；2.重建图像的对比度和尺寸要大；3.算法的速度要快；4.对于具有多维空间的物体有效；5.算法鲁棒性好并且稳定，可以经受噪音或数据不完整的影响；6.实践中散射检测模型要有效；7.具有通用性，可以应用于其他波动方程和散射现象；8.算法具有独立性，只需要少量或不需要关于被测目标的先验知识。目前没有一种非线性算法能够满足上述所有的条件。由于对比源反演算法无须正演计算的特点，在迭代过程中利用快速傅立叶(Fast Fourier)变换计算出三维并矢格林(Dyadic Green)函数算子及其共轭算子，确保了在三维成像中反演过程的效率及稳定性，而且通过其他扩展手段，如分域并行方法(正在研究中)，可以大大提高其计算速度。

发明内容

本发明提出了基于对比源反演算法对三维电磁数据的成像方法，在电磁波的波动方程的基础上建立的符合三维电磁成像的数学模型，利用对比源反演算法无须正演计算的特点，在迭代过程中利用快速傅立叶(Fast Fourier)变换计算出三维并矢格林(Dyadic Green)函数算子及其共轭算子，并同时采用正则化方法和频率跳变(Frequency Hopping，FH)方法扩展了该成像算法，提高重建图像的质量。

本发明的技术方案：

本发明是基于一种非线性算法-对比源反演算法对三维电磁数据的成像方法，该方法读取接受到的三维电磁数据，该数据也包括了与实测电磁数据信息(如背景媒质的电磁参数、成像区域的分辨率、发射器的数量及位置、接收器的数量及位置等)，利用数据构建三维的数据积分方程及初始化目标积分方程，用三维数据和目标方程构成代价泛函，然后利用改进的共轭梯度法求解代价泛函函数的最小值，并采用正则化方法和频率跳变方法提高算法在迭代计算过程中的性能，一旦达到预期条件，则停止迭代并输出重建结果，如对比源函数、对比度值函数及总场。

本发明的有益效果：

本发明设计了一种基于非线性对比源反演算法对三维电磁数据的成像方法。该方法避免了正演计算，有利于反演计算的效率及稳定。并且该方法结合了正则化方法及频率跳变方法等扩展手段，大大提高了重建图像的质量，使算法能够实现对具有(介电常数、电导率等参数)高对比度值的三维目标的精确重建，扩大了算法的应用范围，增强了该算法的实用性；

附图说明

图1是三维成像方法的示意图

图2是三维电磁成像的设置

图3是三维目标1“twocubes”的实际结构

图4是三维目标2“cubepheres”的实际结构

图5采用频率跳变(FH)方法的三维对比源反演算法(CSI)对目标“twocubes”重建结果的实部

图6采用频率跳变(FH)方法的三维对比源反演算法(CSI)对目标“twocubes”重建结果的虚部

图7采用频率跳变(FH)方法和正则化方法的三维对比源反演算法(CSI)对目标“twocubes”重建结果的实部

图8采用频率跳变(FH)方法和正则化方法的三维对比源反演算法(CSI)对目标“twocubes”重建结果的虚部

图9采用频率跳变(FH)方法的三维对比源反演算法(CSI)对目标“cubespheres”重建结果的实部

图10采用频率跳变(FH)方法的三维对比源反演算法(CSI)对目标“cubespheres”重建结果的虚部

图11采用频率跳变(FH)方法和正则化方法的三维对比源反演算法(CSI)对目标“cubespheres”重建结果的实部

图12采用频率跳变(FH)方法和正则化方法的三维对比源反演算法(CSI)对目标“cubespheres”重建结果的虚部

具体实施方式

以下结合附图和通过实施例对本发明的具体实施方式作进一步说明：

本发明设计了基于非线性对比源反演算法对三维电磁数据的成像方法，其特征在于：如图1所示，本方案包括以下步骤：

a.通过最小化代价泛函，计算出对比源；

b.通过把步骤a中得出的对比源的近似值代入到“目标”方程，计算出“兴趣区域”中总场

c.再通过最小化代价泛函，确定最佳的对比度值，即重建图像。

步骤a所述的对比源是指的对比度值函数和总场的积，即

w(R′)＝χ(R′)E(R′)，                (1)

式中χ为对比度值函数(Contrast Function)，E表示总场。

步骤a所述的代价泛函是由两个三维误差函数构成的，这两个误差函数有是在两个积分方程的基础上定义的，即“数据”方程和“目标”方程。如图2所示，在以均质、各向同性的媒质为背景的立方体“兴趣区域”D中存在一个非均质、各向同性的目标S，发射天线和发射天线被放置在D周围的一系列的点上，分别形成圆环形区域T和M。这样三维的电磁“数据”方程和“目标”方程的表达式以此为

$E^{\mathrm{sc}}\left(R\right)=k^2\underset{R^{\′}\∈D}{\∫\∫\∫}\mathrm{\χ}\left(R^{\′}\right)E\left(R^{\′}\right)G(R^{\′}-R)d^3{R}^{\′}+\▿\▿\·\underset{R^{\′}\∈D}{\∫\∫\∫}\mathrm{\χ}\left(R^{\′}\right)E\left(R^{\′}\right)G(R^{\′}-R)d^3{R}^{\′},R\∈M,---\left(2\right)$

$E\left(R\right){=E}^{\mathrm{in}}\left(R\right)=k^2\underset{R^{\′}\∈D}{\∫\∫\∫}\mathrm{\χ}\left(R^{\′}\right)E\left(R^{\′}\right)G(R^{\′}-R)d^3{R}^{\′}+\▿\▿\·\underset{R^{\′}\∈D}{\∫\∫\∫}\mathrm{\χ}\left(R^{\′}\right)E\left(R^{\′}\right)G(R^{\′}-R)d^3{R}^{\′},R\∈D---\left(3\right)$

式中，上标“sc”和“in”分别代表了散射和入射，G表示在真空中电磁波的三维格林函数，k为波数，因为本文中的探测目标由金属和电介质柱体构成，所以可被定义为介电常数和电导率的对比度值函数，即

$\mathrm{\χ}\left(R^{\′}\right)=\frac{{\mathrm{\ϵ}}_r\left(R^{\′}\right)}{{\mathrm{\ϵ}}_r}-1+j\frac{{\mathrm{\σ}}_e\left(R^{\′}\right)}{{\mathrm{\ω\ϵ}}_0{\mathrm{\ϵ}}_r},R^{\′}\∈D---\left(4\right)$

其中ε₀为真空的介电常数，ε_r(R′)和ε_r分别为目标和背景媒质的相对介电常数，σ_e(R′)为目标的电导率，所以可以看出实数部分Re(χ)代表了介电常数的对比度值，虚数部分Im(χ)与角频率ω有关，代表了电导率的对比度值。本说明书中的三维探测目标都是由介质材料构成的，因此重建的对比度值的实数部分应该远远大于其虚数部分。

根据上述的“数据”和“目标”方程，我们就可以计算在三维的“数据”误差函数

和“目标”误差函数

分别表示为

其中和

为积分算子，分别代表了数据方程和目标方程。这样代价泛函可以定义为

式中，

代表了“目标”方程的那一部分代价泛函，

代表了“数据”方程的那一部分代价泛函，

和为加权因子，如下表示为

可以看到代价泛函是对比源和对比度值的函数，假定当前迭代下的对比度值函数

不变，通过最小化

就可以计算出当前迭代下的对比源，这一计算过程是通过共轭梯度法实现的。

在步骤a中得出的对比源的近似值后，在步骤b中将得出的对比源代入到“目标”方程(3)中，就可以计算出当前迭代下的“兴趣区域”中总场。

在得到了当前迭代下的对比源和总场后，这时代价泛函就只是对比度值的函数，即

这样在步骤c中最小化代价泛函

可以确定最佳的对比度值，这一计算过程也是通过共轭梯度法实现的。

在上述过程中还结合了正则化方法及频率跳变方法，具体内容如下：

正则化方法：在电磁成像中，由于噪声的存在及逆散射所固有的病态特性，造成了解的不稳定。为了提高重建图像的质量，我们采用了以全变差方式为基础的乘法正则化(Multiplicative Regularized，MR)的概念，这样既加强了解的稳定性，同时又保持解的特性。与对比源反演算法相比较，乘法正则化的对比源反演算法在代价范函中多了一项全变差因子，即

式中，为全变差因子，可写为

其中V代表区域D的体积，δ²用来保证

的可微性，可表示为

其中Δ³为域D离散化后的子域体积的倒数，即

${\mathrm{\Δ}}^3=\frac{1}{\mathrm{\Δx\Δy\Δz}}.---\left(12\right)$

频率跳变方法：对于频域电磁数据来讲，由于数据包含了多个频率的数据，我们首先用算法对单个频率的数据进行重建成像，在这个频率下取得的重建结果作为算法对下一个较高频率数据重建的初始值，然后不断重复这一过程，依次对这些频率下的数据进行重建，直到算法完成对最高频率数据重建成像，我们称这种方法为频率跳变方法(Frequency HoppingApproach，FH)。

实施例

共有2个非均匀的电介质目标：“twocubes”和“cubespheres”。三维实测电磁数据都是在暗室中获得的，其中发射天线为抛物面天线，接收天线为喇叭天线，试验设置如图2所示，以球坐标为例，发射天线在距试验设置中心的距离为1.796米的球面上移动；方位角φ^s(azimuthal angle)由0°到320°变换，每次变换40°；极角θ^s(polar angle)由50°到150°变换，每次变换15°；这样就共有81个发射天线分布在一个半径为1.796米的球面上。接收天线距试验设置中心的距离也是1.796米；方位角φ^r＝90°，即位于z＝0的平面；极角θ^r表示如下：

θ^r＝θ^s+n×10°+40°，n＝1，2，.....，27，

即对于每一个发射天线，都有27个接收天线分布在z＝0的平面上的半径为1.796米的圆弧上与之对应。

如图3所示，该目标由两个介电常数为2.3(相对于背景媒质，其介电常数的对比度值χ＝1.3)，边长均为25mm的电介质立方体构成，两个立方体在z方向的位置分别为25mm和50mm，高于z＝0的平面。在这我们采用的成像区域D的体积为112×112×112mm³，中心位置为(0，0，50)mm，域D被离散为32×32×32个子域，即子域的边长为3.5mm，以实测数据最高频率8GHz为例，约为波长的十分之一。首先我们用对比源反演算法及频率跳变(CF)方法对目标“twocubes”进行重建，其中选用3GHz、4GHz、5GHz、6GHz、7GHz和8GHz共6个频率的数据，算法对每个频率数据的迭代次数均为100。图5显示了该方法重建结果的实部，图6显示了重建结果的虚部，从上述图中我们可以看到重建结果的实部(介电常数)的对比度值Re(χ)远远大于虚部(电导率)的对比度值Im(χ)，根据式(3)，可以判断出目标是电介质目标。从图中我们还可以看到重建的两个立方体在z方向的位置是不同的，两个立方体在xy平面的位置、大小都重建得非常好，重建的最大对比度值Re(χ)约为1.2，略低于实际的对比度值。对于z方向上不同位置，重建结果也是不同，图5(e)～(m)清晰的显示一个立方体的重建结果，在图5(e)中z＝27.25mm，在图5(m)中z＝51.75mm，两图在z方向的距离为24.5mm，与该立方体在z方向的实际边长基本相符；图5(m)～(t)显示了另一个立方体的重建结果，图5(m)～(t)间的距离为24.5mm，也与实际目标的大小基本相符。图7和图8显示了采用正则化方法和频率跳变方法的对比源反演算法的重建结果的实部和虚部，与图5比较，我们可以看到图7中的两个立方体的形状更加清晰，而且最大的介电常数的对比度值Re(χ)分布更广、更平均，只是重建的最大对比度值Re(χ)约为1.1，低于对比源反演算法。

如图4所示，该目标是由介电常数为2.6(对比度值χ＝1.6)，直径为15.9mm的27个电介质球体聚集在一起而构成的一个边长为47.7mm的立方体(目标在xy平面上的投影是一个由9个圆构成的正方形)，该立方体位于z＝0平面之上，在这我们采用的成像区域D的体积为80×80×80mm³，中心位置为(0，0，12.5)mm，域D被离散为32×32×32个子域。我们采用对比源反演算法及频率跳变方法对目标“cubespheres”进行重建，其中选用3、4、5、6、7和8GHz共6个频率的数据，算法对每个频率数据的迭代次数均为100。图9显示了该方法最终取得的重建结果的实部，图10显示了重建结果的虚部。与重建结果的实部(介电常数的对比度值)Re(χ)比较，虚部(电导率的对比度值)Im(χ)基本可以忽略不计，因此也可以判断出目标同样是电介质目标。随着频率的升高，重建目标的位置，形状变得越来越清楚，可以看到是由9个圆构成的正方形。重建的介电常数的对比度值也越来越准确。在图9中，在xy平面上9个圆构成的正方形的位置、大小都重建得非常好，重建的最大对比度值Re(χ)约为1.5，与实际的对比度值接近。从图9(f)～(g)中，我们可以清楚的看到，对于z方向上不同位置，重建的正方形的尺寸是基本不变的，在图9(f)显示了z＝-11.0935mm上的重建结果，图9(v)显示z＝36.71875mm的位置上，两个图的距离约为47.8mm，与目标在z方向的实际长度接近。图11和12显示了采用正则化方法和频率跳变方法的对比源反演算法所取得的重建结果的实部和虚部。与图9中的重建结果比较，图11中重建的最大对比度值约为0.7，低于实际的对比度值，重建的立方体的整体形状更加清晰，只是没有能够清楚的重建构成立方体的小球体。

Claims (7)

1.一种基于非线性对比源反演算法对三维电磁数据的成像方法，其特征在于：本方案包括以下步骤：

a.通过最小化代价泛函，计算出三维的对比源函数；

b.通过把步骤a中得出的三维的对比源函数的近似值代入到三维“目标”方程，计算出三维“兴趣区域”中总场

c.再通过最小化代价泛函，确定最佳的三维对比度值，即重建三维图像。

2.根据权利要求1所述的一种基于非线性对比源反演算法对三维电磁数据的成像方法，其特征在于：步骤a所述的对比源是指的三维对比度值函数和总场的积，即

w(R′)＝χ(R′)E(R′)

式中χ为三维对比度值函数(Contrast Function)，E表示三维总场；步骤a所述的代价泛函是由两个三维误差函数构成的，这两个误差函数有是在两个三维积分方程的基础上定义的，即“数据”方程和“目标”方程；三维“数据”方程和“目标”方程的表达式以此为

$E^{\mathrm{sc}}\left(R\right)=k^2\underset{R^{\′}\∈D}{\∫\∫\∫}\mathrm{\χ}\left(R^{\′}\right)E\left(R^{\′}\right)G(R^{\′}-R)d^3{R}^{\′}+\▿\▿\·\underset{R^{\′}\∈D}{\∫\∫\∫}\mathrm{\χ}\left(R^{\′}\right)E\left(R^{\′}\right)G(R^{\′}-R)d^3{R}^{\′},R\∈M,$

$E\left(R\right){=E}^{\mathrm{in}}\left(R\right)=k^2\underset{R^{\′}\∈D}{\∫\∫\∫}\mathrm{\χ}\left(R^{\′}\right)E\left(R^{\′}\right)G(R^{\′}-R)d^3{R}^{\′}+\▿\▿\·\underset{R^{\′}\∈D}{\∫\∫\∫}\mathrm{\χ}\left(R^{\′}\right)E\left(R^{\′}\right)G(R^{\′}-R)d^3{R}^{\′},R\∈D$

其中D为“兴趣区域”即成像区域，M为接收天线所在的区域，上标“sc”和“in”分别代表了散射和入射，G表示在真空中电磁波的三维格林函数，k为波数，因为本文中的探测目标由电介质构成，所以可被定义为介电常数和电导率的对比度值函数，即

$\mathrm{\χ}\left(R^{\′}\right)=\frac{{\mathrm{\ϵ}}_r\left(R^{\′}\right)}{{\mathrm{\ϵ}}_r}-1+j\frac{{\mathrm{\σ}}_e\left(R^{\′}\right)}{{\mathrm{\ω\ϵ}}_0{\mathrm{\ϵ}}_r},R^{\′}\∈D$

其中ε₀为真空的介电常数，ε_r(R′)和ε_r分别为目标和背景媒质的相对介电常数，σ_e(R′)为目标的电导率，所以可以看出实数部分Re(χ)代表了介电常数的对比度值，虚数部分Im(χ)与角频率ω有关，代表了电导率的对比度值。

3.根据权利要求2所述的三维目标和数据方程，其特征在于：根据上述的“数据”和“目标”方程，我们就可以计算三维“数据”误差函数

和“目标”误差函数

其中

和

为三维积分算子，分别代表了数据方程和目标方程；根据上述的三维误差方程，该代价泛函可以定义为

式中，

代表了三维“目标”方程的那一部分代价泛函，

代表了三维“数据”方程的那一部分代价泛函，

和

为加权因子，如下表示为

4.根据权利要求3所述的代价泛函，其特征在于：该代价泛函是三维对比源和三维对比度值的函数，假定当前迭代下的三维对比度值函数

不变，通过最小化就可以计算出当前迭代下的三维对比源，这一计算过程是通过共轭梯度法实现的。

5.根据权利要求1所述的一种基于非线性对比源反演算法对三维电磁数据的成像方法，其特征在于：在步骤a中得出的三维对比源的近似值后，在步骤b中将得出的对比源代入到三维“目标”方程中，就可以计算出当前迭代下“兴趣区域”中总场。

6.根据权利要求1所述的一种基于非线性对比源反演算法对三维电磁数据的成像方法，其特征在于：在得到了当前迭代下的三维的对比源和总场后，这时代价泛函就只是三维对比度值的函数，即这样在步骤c中最小化代价泛函可以确定最佳的三维对比度值，这一计算过程也是通过共轭梯度法实现的。

7.根据权利要求1所述的分为3个步骤a，b，c的整个迭代过程，其特征在于：整个过程中结合了频率跳变(FH)方法和正则化方法，具体内容如下：

频率跳变方法：对于频域电磁数据来讲，由于数据包含了多个频率的数据，我们首先用算法对单个频率的数据进行重建成像，在这个频率下取得的重建结果作为算法对下一个较高频率数据重建的初始值，然后不断重复这一过程，依次对这些频率下的数据进行重建，直到算法完成对最高频率数据重建成像，我们称这种方法为频率跳变方法(Frequency HoppingApproach，FH)；

正则化方法：在电磁成像问题中，由于噪声的存在及逆散射所固有的不适定性，造成了解的不稳定，这些都影响了重建图像的质量；受到图像处理算法中正则化方法的启发，我们在对比源反演算法中采用了以全变差(Total Variation)为基础的乘法正则化，这样既加强了解的稳定性，同时又保持解的特性；由于在算法实现过程中把正则化因子乘在了代价泛函上，因此扩展后的算法称为乘法正则化(Multiplicative Regularization，MR)的对比源反演算法；在乘法正则化的对比源反演算法中，代价泛函可如下表示为

式中，

为全变差因子，可写为

其中V代表区域D的体积，δ²用来保证的可微性，可表示为

其中Δ³为域D离散化后的子域面积的倒数，即

${\mathrm{\Δ}}^3=\frac{1}{\mathrm{\Δx\Δy\Δz}}.$

Images