CN102955159A - 一种基于压缩感知的电磁逆散射成像方法 - Google Patents

Abstract

本发明公开了一种基于压缩感知的电磁逆散射成像方法，涉及微波成像技术，包括：步骤一、微波激励与测量阶段，在目标外用微波对目标进行穿透照射，测量目标的散射回波，从而去感知目标内部结构和材质；步骤二、目标建模阶段，建立目标介电常数和散射回波间的非线性观测模型，建立描述目标内部结构稀疏性的表征模型；步骤三、计算成像阶段，使用目标模型和压缩感知处理方法对目标的介电常数分布成像，由于不同介电常数对应目标的不同组成材料，故成像结果作为目标内部结构的表征。相对于传统微波成像方法，本发明使用压缩感知处理技术，显著降低回波观测数据量，同时显著提高目标内部结构的成像清晰度。

Description

一种基于压缩感知的电磁逆散射成像方法

技术领域

本发明涉及压缩感知和基于电磁逆散射的微波成像技术，特别是一种通过微波照射目标得到电磁波回波，进而利用回波数据反演目标的介电常数分布。本发明不仅利用到了回波数据与介电常数分布间的非线性关系，还利用到介电常数分布自身的稀疏特性。相对于传统微波成像方法，它可以提高目标内部图像的轮廓分辨清晰度。

背景技术

定量的电磁逆散射方法，旨在重建目标的电磁参数分布，如介电常数、电导率等，在重建目标电磁参数分布的同时也得到了目标的几何图像。它由严格的麦克斯韦方程出发，充分考虑目标之间的多次散射及目标与背景媒质之间的相互作用，能够高精度地重建目标几何参数和电磁参数，得到目标的高分辨率图像。

目前发展的高维定量的电磁逆散射方法可分为线性逆散射方法和非线性逆散射方法两类。线性逆散射方法工作在目标对比度低的情况下，其测量系统模型为线性模型，典型的反演方法为衍射层析成像(Diffraction Tomography，DT)。非线性逆散射方法适用于高对比度目标反演，典型的包括波恩迭代方法(BIM)和变形波恩迭代方法(DBIM)，对比源方法(Contrast Source Inversion，CSI)等等。

2006年，由Donoho、Candès及华裔科学家Tao等人提出了一种新的信息获取理论，即，压缩感知(Compressive Sensing，CompressedSensing，CS)理论。该理论指出：对可压缩信号可通过远低于Nyquist采样率的测量数据精确地恢复。新理论的优势是能凭借被测信号的稀疏性大幅度压缩提取该信号所必需的测量数。该理论一经提出，即在多个信号处理领域如图像处理、医疗成像、模式识别、地质勘探、光学/雷达成像、无线通信等得到应用，并被美国科技评论评为2007年度十大科技进展。尤其在各种成像应用中，利用CS技术可以以极少的线性测量得到无模糊的目标图像，降低了成像系统的测量消耗和系统复杂度。

实际应用中待探测和成像的目标具有一定的规则性，如在总体变分或者小波变换下表示成稀疏的或可压缩的形式，因此可以将CS理论与实现技术应用于线性电磁逆散射成像中。如前所述，CS技术一般只适用于线性测量系统，目前未有见CS技术与非线性电磁逆散射结合的电磁成像方法出现，这限制了CS技术在电磁逆散射中的应用潜力。基于电磁逆散射的微波成像需要实现目标的非侵入式成像，因此它必须将目标的内部信息映射到目标外部周界上，这种测量方式天然构成一组压缩式的测量。传统的逆散射处理方法通过引入一些人为的先验信息来约束介电常数分布的恢复过程，例如约束信号能量和信号的符号特性，所得到的图像会出现一定程度的模糊，这种模糊在内部跳变边缘体现得最明显。相对于传统的处理方法，基于CS的电磁逆散射以目标内部不连续跳变边缘的稀疏性作为先验约束，能得到轮廓清晰度更高的目标图像。

参考资料：

[1]A.Tikhonov.On the solution of incorrectly formulatedproblems and the regularization method.Seviet Math，Doklady，1963，4：1035-1038.

[2]黄卡玛，赵翔.电磁场中的逆问题及应用.北京：科学出版社，2005.

[3]丁科，宋守根，谢忠球.逆散射理论的发展及应用前景.地球物理学进展，2005，20(3)：661-666.

[4]S.X.Pan，A.C.Kak.A computat ionalstudy ofreconstruction algorithms for diffraction tomography：Interpolation versus filtered backpropagation.IEEE Trans.Acoust.，Speech，SignalProcess.，1983，31：1262-1275.

[5]R.F.Remis，P.M.van den Berg.On the equivalence ofthe Newton-Kantorovich and distorted Born methods.InverseProblems，2000，16(1)：L1-L4.

[6]R.E.Kleinman，P.M.Vandenberg.A ModifiedGradient-Method for 2-Dimensional Problems in Tomography.Journal of Computational and Applied Mathematics，1992，42：17-35.

[7]P.M.van den Berg，R.E.Kleinman.A Contrast SourceInversion Method.Inverse Problems，1997，13：1607-1620.

[8]D.Donoho.Compressed sens ing.IEEE Trans.Inform.，2006，52(4)：5406-5425.

[9]E.Candès，M.Wakin.An  Introduction To CompressiveSampling.IEEE Signal Processing Magazine，2008，25(2)：21-30.

发明内容

本发明的目的是公开一种基于压缩感知的电磁逆散射成像方法，通过目标外部散射场和目标对比度函数之间的非线性函数关系构建压缩测量模型，将待恢复信号描述为均匀缓变部分和不连续跳变部分这两部分的叠加，约束其总体变分的稀疏性，即跳变部分的稀疏性。该方法结合多输入多输出、单点频或扫频的测量，以及待恢复信号的稀疏性，利用压缩感知的处理手段恢复原始介电常数分布，有效保持内部结构边缘的成像清晰度。在具体处理中，使之满足压缩感知的一般处理形式；同时，避免每迭代步骤中出现复杂的优化计算。

为了达到上述目的，本发明的技术解决方案是：

一种基于压缩感知的电磁逆散射成像方法，其包括：

步骤一、微波激励与测量阶段，在该阶段由多个发射机从不同照射角对目标实施单一频率的微波照射，或在设计频段内实施扫频照射，并由多个接收机从不同观测角对目标散射回波的幅度和相位信息，进行对应的单频或扫频测量，发射机和接收机一一对应，构成多发多收测量；

步骤二、目标建模阶段，建立目标介电常数和散射回波间的非线性观测模型，建立描述目标内部结构稀疏性的表征模型；

步骤三、计算成像阶段，使用目标模型和压缩感知处理方法对目标的介电常数分布成像，成像结果作为目标内部结构的表征。

所述的电磁逆散射成像方法，其特征在于，所述步骤一中，在微波激励与测量阶段获得多发多收、单频或多频的目标回波测量数据；包括：

1a、首先由发射天线在环绕目标的围线上，或者目标一侧的多个发射位置上分别发射特定频率或频段的电磁波，并由多个接收天线接收响应电场E^tot(r_r，r_i，ω)，其中，r_i为发射机位置坐标，r_i为接收机位置坐标，ω为工作角频率；

1b、然后将目标移出被测区域，在发射源不变的前提下测量相同接收位置上的直达电场E^inc(r_i，r_i，ω)；

1c、最后将两次测量数据相减，得到目标散射回波数据E^scatt(r_r，r_i，ω)，即：

E^scatt(r_r，r_i，ω)＝E^tot(r_r，r_i，ω)-E^inc(r_r，r_i，ω)          (1)

所述的电磁逆散射成像方法，其所述步骤二中，建立目标介电常数分布ε(r)与散射回波E^scatt(r_r，r_i，ω)之间的非线性观测模型，其中，r为目标场景位置坐标，r_r为接收机位置坐标，r_i为发射机位置坐标，ω为观测角频率；

根据探测微波穿透目标表面后，诱发目标内部不同部分间多次散射的作用机理，得到描述该关联关系的目标内部场效应模型和目标外部场效应模型，其中，内部场效应模型描述为：

$E^{\mathrm{tot}}(r,r_i,\mathrm{\ω})=E^{\mathrm{inc}}(r,r_i,\mathrm{\ω})+$

$\underset{D_{\mathrm{inv}}}{\∫}G(r,r^{\′};k_b,\mathrm{\ω})k_b^2(x^{\′},\mathrm{\ω})\mathrm{\χ}(r^{\′},\mathrm{\ω})E^{\mathrm{tot}}(r,r_i,\mathrm{\ω})\·{\mathrm{dr}}^{\′},r\∈D_{\mathrm{inv}}---\left(2\right)$

其中，D_inv为目标场景，r，r′为目标场景坐标，E^inc(r，r_i，ω)为D_inv内的激励电场或称入射场，E^tot(r，r_i，ω)为D_inv内的应激电场或称总场，G(r，r′；k_b，ω)是背景波数为k_b情况下的电磁波传递函数，称为格林函数，x(r′，ω)为目标介电常数定义为

其中，k为目标的波数；

外部场效应模型描述为：

$E^{\mathrm{scatt}}(r,r_i,\mathrm{\ω})=\underset{D_{\mathrm{inv}}}{\∫}G(r,r^{\′};k_b,\mathrm{\ω})k_b^2(x^{\′},\mathrm{\ω})\mathrm{\χ}(r^{\′},\mathrm{\ω})E^{\mathrm{tot}}(r,r_i,\mathrm{\ω})\·{\mathrm{dr}}^{\′}---\left(3\right)$

结合两组模型得到离散化的非线性观测模型为：

E^scatt(ω)＝A(x(ω))＝G_d(ω)X(ω)(I-G_s(ω)X(ω))^-1E^inc(ω)(4)

其中，X(ω)为由分布x(r′，ω)为对角线元素的对角矩阵，I为合适规模的单位矩阵，G_s，G_d分别表示格林函数在内、外部模型中的离散矩阵形式，不考虑频率ω，观测模型可简记为E^scatt＝A(x)，A是一个非线性函数。

所述的电磁逆散射成像方法，其所述步骤二中，建立描述目标内部结构稀疏性的表征模型；该表征模型由目标内部相同材质区域介电常数缓变的部分，与不同材质交界处介电常数跳变的部分这两部分刻画，用对比度函数x(r′，ω)方向导数的绝对值之和约束：

$D\left(\mathrm{\χ}\right)=\underset{D_{\mathrm{inv}}}{\∫}|\▿\mathrm{\χ}(r^{\′},\mathrm{\ω})|\·{\mathrm{dx}}^{\′}---\left(5\right)$

其中，

为x(r′，ω)的方向导数，表征目标内部不同结构的轮廓，相对整幅目标内部图像是稀疏的。

所述的电磁逆散射成像方法，其所述步骤三中，在计算成像阶段实现对目标内部的计算成像，即由下述非线性约束优化过程：

$\mathrm{\χ}=\underset{\mathrm{\χ}}{\arg }\min D\left(\mathrm{\χ}\right)$ s.t.E^scatt＝A(x)            (6)

其中，

是最小化计算式，s.t.是“subject to”的缩写，表示同时满足，得到介电常数分布，作为目标内部的恢复图像。

所述的电磁逆散射成像方法，其所述步骤一中，对于单一频率测量，在步骤三优化过程(6)式的求解分为初值假设步骤，迭代优化步骤，和终止步骤：

1a、初值假设步骤，以已知背景介电常数分布作为目标分布初值，即在该步骤假设目标对比度为零；

1b、迭代优化步骤，设第n迭代步的计算结果为x_n，n＞1，为了简化(6)式的计算，(6)式中的D(x)用以下形式代替，即：

$D\left(\mathrm{\χ}\right)\←{\∫}_D\frac{{|\▿\mathrm{\χ}|}^2}{{({|\▿{\mathrm{\χ}}_{n-1}|}^2+\mathrm{\δ})}^{\frac{1}{2}}}{\mathrm{dr}}^{\′}---\left(7\right)$

其中，δ为人为引入的正值常数，取其为D_inv离散网格面积与工作波长平方的比值；(6)式中的A(x)用A(x)在x_n-1处的一阶泰勒展开代替，即：

$A\left(\mathrm{\χ}\right)\←A\left({\mathrm{\χ}}_{n-1}\right)+\▿A\·(\mathrm{\χ}-{\mathrm{\χ}}_{n-1})---\left(8\right)$

得到：

${\mathrm{\χ}}_n=\underset{\mathrm{\χ}}{\arg }\min \left|\right|\mathrm{W\χ}\left|\right|$ s.t. $y_n=\▿A\·\mathrm{\χ}---\left(9\right)$

其中，算子 $W=\frac{1}{\sqrt{{\left|{\▿\mathrm{\χ}}_{n-1}\right|}^2+\mathrm{\δ}}}\▿,$ 向量 $y_n=E^{\mathrm{scatt}}-A\left({\mathrm{\χ}}_{n-1}\right)+{\▿A}_{n-1}\·{\mathrm{\χ}}_{n-1},$

(8)式的解为：

${\mathrm{\χ}}_n={\left(W^{*}W\right)}^{-1}{(\▿A)}^{*}{(\▿A{\left(W^{*}W\right)}^{-1}{(\▿A)}^{*})}^{-1}{y}_n---\left(10\right)$

(9)式为迭代步骤的迭代表达式；

1c、终止步骤，当

时终止迭代，其中，η为人为设置的终止门限，取为η＝10^-3。

所述的电磁逆散射成像方法，其所述步骤一中，对于扫频测量，在步骤二中设目标介电常数在扫频带宽内缓慢变化，扫频测量以提高分辨率；

测量时，若有校准后的矢量网络分析仪，利用网络分析仪发射扫频连续波，直接得到扫频测量数据；若有校准后的多个发射、接收机组件，利用多个发射机发射扫频带宽内的不同正交编码脉冲信号，接收机接收、解码，通过时-频傅里叶变换得到扫频测量数据；

计算时，利用测得的扫频数据，首先计算低频率下的目标对比度分布，再以此分布为初值，计算高频率下的目标对比度分布，最后以高频率下的计算结果作为最终成像结果。

所述的电磁逆散射成像方法，其所述接收天线接收响应电场，包括幅度与相位信息。

所述的电磁逆散射成像方法，其适用于不包含金属目标的场景成像：电感应测井和生物体软组织切片成像的应用。

相对于传统的电磁逆散射成像方法，本发明方法的优势包括以下几个方面：

1.本发明将CS技术应用于变形波恩迭代方法中(DBIM)，相较于传统DBIM，可以有效保持重建目标图像的边缘轮廓清晰图，实现高对比度目标的非线性无模糊成像。这是由于在传统DBIM迭代过程中对重建目标施加了可稀疏化的先验约束条件，限制重建结果必需在总体变分或小波变换下有最稀疏的表示。

2.本发明实现CS应用于非线性的电磁逆散射成像，扩展了CS技术在电磁逆散射中的使用范围。具体实现时，通过变形波恩迭代将原始的非线性测量模型近似为分段线性模型，在每个线性分段中使用CS的处理手段得到稀疏约束下的成像结果，再通过逐步迭代收敛到最终的目标图像。整个过程既保持了变形波恩迭代收敛速度快、迭代步数少的优点，又使用到CS处理技术在重建稀疏信号/目标时的强大能力，从而实现高对比度目标的有效重建。

3.本发明在使用CS处理手段同时，保持新成像方法与传统变形波恩迭代的计算复杂度基本相同，能实现对测量数据地实时处理与成像。这是由于，在变形波恩迭代的每个线性分段中，没有完整求解一个CS问题，而借助于前一步的反演结果松弛了原始的CS问题。松弛后的CS问题能够通过矩阵运算直接求解，相较于传统方法，其计算量没有显著提升。

4.本发明使用CS技术，可以缩减成像所必需的测量数据量。例如，在有限的发射角和有限的观测角中进行目标的散射场测量。而在通常电磁逆散射中，限制发射和观测角度会造成目标信息的丢失，进而导致目标成像结果变模糊或者缺失，出现伪影，甚至不可接受的结果。而在应用CS技术后，目标在已知变换下可稀疏化的先验条件有可能补充丢失的信息量，从而保持重建目标图像的正确性。

本发明使用的方法在获得上述优势的同时，也继承了变形波恩迭代的优势，例如1)迭代步数少，收敛快；2)单步迭代可用快速算法加速其中的矩阵运算；3)计算过程可并行化。

本发明可应用于如下几个重要领域：1)人体组织切片成像，例如对人体脑腔、胸腔等部位的成像，检测早期癌变。人体的切片图像通常可被稀疏化，利用本发明的方法可以快速有效地实现病灶的定位与定性。2)地底埋藏物和地下矿产的探测。3)电离层监控。4)无损检测。等等。

附图说明

图1是本发明一种基于压缩感知的电磁逆散射成像方法的微波激励与测量系统示意图；

图2是本发明一种基于压缩感知的电磁逆散射成像方法的扫频测量系统示意图；其中：

图2(a)为基于矢量网络分析仪-7扫频测量；

图2(b)为基于多个发射接收机组件-8组成多通道测量；

图3是本发明方法的实施流程图；

图4是本发明方法中非线性测量模型示意图，非线性测量模型(图4上)分解为内部场效应模型(图4左下)和外部场效应模型(图4右下)；

图5是本发明方法中一个成像目标及其稀疏表征模型实例，其中：

图5(a)为成像目标；

图5(b)为稀疏表征模型；

图6是在工作频率4GHz下实验数据成像结果；其中：

图6(a)为成像目标；

图6(b)为传统微波逆散射方法成像结果；

图6(c)为本发明使用压缩感知处理的成像结果；

图6(d)为目标对比度恢复误差随迭代步骤的变化示意图。

图中标号：

1-发射天线          2-接收天线           3-目标

4-成像区域          5-成像转台           6-接收天线转轴

7-矢量网络分析仪    8-发射/接收机组件    9-等效辐射源

具体实施方式

下面结合实施范例和附图说明本发明的具体实施方式。应当指出，设计该施例的主要目的是通过设置最简单的电磁逆散射问题模型，集中体现本发明的基本处理流程，它对本发明的具体应用不起任何限定作用。

如图1所示，是微波激励与测量的系统几何图，系统采用发射天线1从多个位置依次发射单频或扫频电磁波，另有接收天线2从多个位置接收回波信号，目标3放置于固定的成像区域4中。

图2是图1所示测量的两种可选的实施方式。实施方式一如图2a)所示，由矢量网络分析仪7配合两副天线实施测量，利用矢量网络分析仪7的端口1接入发射天线发射1，端口2接入接收天线2接收回波信号，以1-2端口间的微波网络散射参数S21作为测量数据，目标位于测量轨道中心的成像转台5上。由于在测量中只使用两副天线，测量时固定发射天线1位置不动，旋转接收天线转轴6使接收天线2在多个位置接收。完成一轮数据接收后多次转动成像转台5，重复同样的接收过程。实施方式二如图2b)所示，由多个发射/接收机组件8同时实施测量，利用发射机向目标发射正交编码信号，接收机接收后解码，得到多发多收的测量数据。与前一测量实施方式相比较，通过增加测量硬件节省了测量时间消耗。

结合图3，基于压缩感知的电磁逆散射成像方法的具体处理步骤为

步骤S1，微波激励与测量阶段，按照图1、2所示的测量系统结构和实施方式布置测量：

1a、首先由发射天线1在环绕目标3的围线上，或者目标3一侧的多个发射位置上分别发射特定频率或频段的电磁波，并由接收天线2接收响应电场(包括幅度与相位信息)E^tot(r_r，r_i，ω)，其中，r_i为发射天线1位置坐标，r_i为接收天线2位置坐标，ω为工作角频率；

1b、然后将目标3移出被测区域，在发射不变的前提下测量相同接收位置上的直达电场E^inc(r_r，r_i，ω)；

1c、最后将两次测量数据相减，得到目标3散射回波数据E^scatt(r_r，r_i，ω)，即：

E^scatt(r_r，r_i，ω)＝E^tot(r_r，r_i，ω)-E^inc(r_r，r_i，ω)

步骤S2，目标建模阶段：

2a、建立目标介电常数分布ε(r)与散射回波E^scatt(r_r，r_i，ω)之间的非线性观测模型。如图4所示，目标3的电磁波散射过程可等效的分解为，目标3受电磁波E^inc(r)激励产生内建电场E^tot(r)，r∈D_inv，以及等效辐射源9x(r)E^tot(r)，再由该等效源发射散射场E^scatt(r)＝E^tot(r)-E^inc(r)，r∈D_obs这两个过程。记测量线上的坐标集合为D_obs，依据空间坐标r从属于D_inv还是D_obs，可将散射方程分别表示成内部场效应方程(状态方程)和外部场效应方程(数据方程)。

设目标3由非磁性、各向同性非均匀媒质组成，分布于给定的背景域内，称该区域为目标区域，记为D_inv。以D_inv的中心为坐标原点，其所在面为坐标平面，建立直角坐标系。记背景介电常数分布为ε_b＝ε_b(r)，对一般媒质，磁导率μ_b≈μ₀近似为恒定值，μ₀为真空中的磁导率。记目标3复介电常数分布

$\widetilde{\mathrm{\ϵ}}\left(r\right)=\left\{\begin{array}{cc}\mathrm{\ϵ}\left(r\right)+\frac{\mathrm{\σ}\left(r\right)}{\mathrm{i\ω}}& r\∈D_{\mathrm{inv}}\\ 0,& r\∉D_{\mathrm{inv}}\end{array}\right.$

其中r为空间坐标，为广义复介电常数，ε(r)为目标3介电常数，当有介质损耗时为复数，σ(r)为目标3电导率，ω为探测波角频率，对单一频率，ω为给定常数，目标磁导率设为μ＝μ₀。

记坐标r上的总场及入射场分别为E^tot(r)、E^inc(r)，目标和背景波数为k，k_b，并且

$k^2\left(r\right)={\mathrm{\ω}}^2\widetilde{\mathrm{\ϵ}}\left(r\right){\mathrm{\μ}}_0$

$k_b^2\left(r\right)={\mathrm{\ω}}^2{\mathrm{\ϵ}}_b\left(r\right){\mathrm{\μ}}_0$

定义目标3的对比度函数为，

$\mathrm{\χ}\left(r\right)=\frac{k^2}{k_b^2}-1=\frac{\widetilde{\mathrm{\ϵ}}\left(r\right)}{{\mathrm{\ϵ}}_b\left(r\right)}-1$

它表示目标3参数与背景参数之间的相对差值。在频率域，目标3的散射方程能够简单的表示成格林函数积分的形式。记给定背景下的格林函数为G(r，r′；ε_b)，对自由背景：

$G(r,r^{\′};\mathrm{\ω})=-\frac{i}{4}{H}_0^{\left(2\right)}\left(k_b\right|r-r^{\′}\left|\right)$

其中，

为第二类零阶贝塞尔函数。则由麦克斯韦方程可得到散射积分方程：

$E^{\mathrm{tot}}\left(r\right)=E^{\mathrm{inc}}\left(r\right)+$

$\underset{D_{\mathrm{inv}}}{\∫}G(r,r^{\′};{\mathrm{\ϵ}}_b)k_b^2\left(r^{\′}\right)\mathrm{\χ}\left(r^{\′}\right)E^{\mathrm{tot}}\left(r^{\′}\right)\·{\mathrm{dr}}^{\′}$

上述方程对所有的空间坐标r都成立。

当r∈D_inv时，E^tot(r)未知，因此内部场效应方程写成：

$E^{\mathrm{scatt}}(r,r_i,\mathrm{\ω})=\underset{D_{\mathrm{inv}}}{\∫}G(r,r^{\′};k_b,\mathrm{\ω})k_b^2(x^{\′},\mathrm{\ω})\mathrm{\χ}(r^{\′},\mathrm{\ω})E^{\mathrm{tot}}(r,r_i,\mathrm{\ω})\·{\mathrm{dr}}^{\′}$

两者结合得到目标介电常数分布ε(r)与散射回波E^scatt(r_r，r_i，ω)之间的非线性观测模型：

$E^{\mathrm{scatt}}\left(r\right)=\underset{D_{\mathrm{inv}}}{\∫}G(r,r^{\′};{\mathrm{\ϵ}}_b)k_b^2\left(r^{\′}\right)\mathrm{\χ}\left(r^{\′}\right)E^{\mathrm{tot}}\left(r^{\′}\right)\·{\mathrm{dr}}^{\′},r\∈D_{\mathrm{obs}}$

2b、建立描述目标内部结构稀疏性的表征模型。如图5所示，目标包括内部相同材质区域介电常数缓变的部分，与不同材质交界处介电常数跳变的部分：图5a)为原始的目标图像，表示为x(r′，ω)，其中缓变部分占主要成分；图5b)为目标图像的跳变边缘，由x(r′，ω)的方向导数刻画，表示为

体现目标结构的稀疏性，总体上用

的绝对值之和描述目标内部结构稀疏性的表征模型：

$D\left(\mathrm{\χ}\right)=\underset{D_{\mathrm{inv}}}{\∫}|\▿\mathrm{\χ}(r^{\′},\mathrm{\ω})|\·{\mathrm{dx}}^{\′}$

步骤S3，计算成像阶段

1a、初值假设步骤，以已知背景介电常数分布作为目标分布初值，即在该步骤假设目标对比度为零；

1b、迭代优化步骤，设第n迭代步的计算结果为x_n，n＞1，上式中的D(x)用以下形式代替，即：

$D\left(\mathrm{\χ}\right)\←{\∫}_D\frac{{|\▿\mathrm{\χ}|}^2}{{({|\▿{\mathrm{\χ}}_{n-1}|}^2+\mathrm{\δ})}^{\frac{1}{2}}}{\mathrm{dr}}^{\′}$

其中，δ为人为引入的正值常数，取其为D_inv离散网格面积与工作波长平方的比值；(6)式中的A(x)用A(x)在x_n-1处的一阶泰勒展开代替，即：

$A\left(\mathrm{\χ}\right)\←A\left({\mathrm{\χ}}_{n-1}\right)+\▿A\·(\mathrm{\χ}-{\mathrm{\χ}}_{n-1})$

得到：

${\mathrm{\χ}}_n=\underset{\mathrm{\χ}}{\arg }\min \left|\right|\mathrm{W\χ}\left|\right|$ s.t. $y_n=\▿A\·\mathrm{\χ}$

其中，算子 $W=\frac{1}{\sqrt{{\left|{\▿\mathrm{\χ}}_{n-1}\right|}^2+\mathrm{\δ}}}\▿,$ 向量 $y_n=E^{\mathrm{scatt}}-A\left({\mathrm{\χ}}_{n-1}\right)+{\▿A}_{n-1}\·{\mathrm{\χ}}_{n-1},$

上式的解为：

${\mathrm{\χ}}_n={\left(W^{*}W\right)}^{-1}{(\▿A)}^{*}{(\▿A{\left(W^{*}W\right)}^{-1}{(\▿A)}^{*})}^{-1}{y}_n$

为迭代步骤的迭代表达式；

1c、终止步骤，当

时终止迭代，其中，η为人为设置的终止门限，取为η＝10^-3。

图6是利用实测数据的成像结果。该实验采样图2a)的测量结构，发射机固定不动，接收机在θ_r＝60°～300°范围内转动，接收机在θ_r＝60°～300°范围内转动，这样就获得不同角度的测量数据。完成一轮测量后转动目标转台，每间隔45°重复一轮测量。测量系统采样的工作频率为4GHz，成像区域为16cm×16cm的方形区域，包含如图6a)所示的介质圆柱，半径为1.5cm，相对介电常数为ε_r＝3±0.3，对称与坐标原点沿x轴放置，间隔9cm。发射机距离目标区域中心0.72m，接收机为0.76m。

图6是在工作频率4GHz下实验数据的成像结果，图6a)成像目标，图6b)传统逆散射方法成像结果，图6c)本发明使用压缩感知处理的成像结果，图6d)目标对比度恢复误差随迭代步骤的变化。

Claims (9)

1.一种基于压缩感知的电磁逆散射成像方法，其特征在于，包括：

步骤一、微波激励与测量阶段，在该阶段由多个发射机从不同照射角对目标实施单一频率的微波照射，或在设计频段内实施扫频照射，并由多个接收机从不同观测角对目标散射回波的幅度和相位信息，进行对应的单频或扫频测量，发射机和接收机一一对应，构成多发多收测量；

步骤二、目标建模阶段，建立目标介电常数和散射回波间的非线性观测模型，建立描述目标内部结构稀疏性的表征模型；

步骤三、计算成像阶段，使用目标模型和压缩感知处理方法对目标的介电常数分布成像，成像结果作为目标内部结构的表征。

2.根据权利要求1所述的电磁逆散射成像方法，其特征在于，所述步骤一中，在微波激励与测量阶段获得多发多收、单频或多频的目标回波测量数据；包括：

1a、首先由发射天线在环绕目标的围线上，或者目标一侧的多个发射位置上分别发射特定频率或频段的电磁波，并由多个接收天线接收响应电场E^tot(r_r，r_i，ω)，其中，r_i为发射机位置坐标，r_i为接收机位置坐标，ω为工作角频率；

1b、然后将目标移出被测区域，在发射源不变的前提下测量相同接收位置上的直达电场E^inc(r_r，r_i，ω)；

1c、最后将两次测量数据相减，得到目标散射回波数据E^scatt(r_r，r_i，ω)，即：

E^scatt(r_r，r_i，ω)＝E^tot(r_r，r_i，ω)-E^inc(r_r，r_i，ω)        (1)

3.根据权利要求1所述的电磁逆散射成像方法，其特征在于：所述步骤二中，建立目标介电常数分布ε(r)与散射回波E^scatt(r_r，r_i，ω)之间的非线性观测模型，其中，r为目标场景位置坐标，r_r为接收机位置坐标，r_i为发射机位置坐标，ω为观测角频率；

根据探测微波穿透目标表面后，诱发目标内部不同部分间多次散射的作用机理，得到描述该关联关系的目标内部场效应模型和目标外部场效应模型，其中，内部场效应模型描述为：

$E^{\mathrm{tot}}(r,r_i,\mathrm{\ω})=E^{\mathrm{inc}}(r,r_i,\mathrm{\ω})+$

$\underset{D_{\mathrm{inv}}}{\∫}G(r,r^{\′};k_b,\mathrm{\ω})k_b^2(x^{\′},\mathrm{\ω})\mathrm{\χ}(r^{\′},\mathrm{\ω})E^{\mathrm{tot}}(r,r_i,\mathrm{\ω})\·{\mathrm{dr}}^{\′},r\∈D_{\mathrm{inv}}---\left(2\right)$

其中，D_inv为目标场景，r，r′为目标场景坐标，E^inc(r，r_i，ω)为D_inv内的激励电场或称入射场，E^tot(r，r_i，ω)为D_inv内的应激电场或称总场，G(r，r′；k_b，ω)是背景波数为k_b情况下的电磁波传递函数，称为格林函数，x(r′，ω)为目标介电常数定义为

$\mathrm{\χ}(x^{\′},\mathrm{\ω})=\frac{k^2(x^{\′},\mathrm{\ω})}{k_b^2(x^{\′},\mathrm{\ω})}-1,$ 其中，k为目标的波数；

外部场效应模型描述为：

$E^{\mathrm{scatt}}(r,r_i,\mathrm{\ω})=\underset{D_{\mathrm{inv}}}{\∫}G(r,r^{\′};k_b,\mathrm{\ω})k_b^2(x^{\′},\mathrm{\ω})\mathrm{\χ}(r^{\′},\mathrm{\ω})E^{\mathrm{tot}}(r,r_i,\mathrm{\ω})\·{\mathrm{dr}}^{\′}---\left(3\right)$

结合两组模型得到离散化的非线性观测模型为：

E^scatt(ω)＝A(x(ω))＝G_d(ω)X(ω)(I-G_s(ω)X(ω))^-1E^inc(ω)(4)

其中，X(ω)为由分布x(r′，ω)为对角线元素的对角矩阵，I为合适规模的单位矩阵，G_s，G_d分别表示格林函数在内、外部模型中的离散矩阵形式，不考虑频率ω，观测模型可简记为E^scatt＝A(x)，A是一个非线性函数。

4.根据权利要求1所述的电磁逆散射成像方法，其特征在于：所述步骤二中，建立描述目标内部结构稀疏性的表征模型；该表征模型由目标内部相同材质区域介电常数缓变的部分，与不同材质交界处介电常数跳变的部分这两部分刻画，用对比度函数x(r′，ω)方向导数的绝对值之和约束：

$D\left(\mathrm{\χ}\right)=\underset{D_{\mathrm{inv}}}{\∫}|\▿\mathrm{\χ}(r^{\′},\mathrm{\ω})|\·{\mathrm{dx}}^{\′}---\left(5\right)$

其中，

为x(r′，ω)的方向导数，表征目标内部不同结构的轮廓，相对整幅目标内部图像是稀疏的。

5.根据权利要求1所述的电磁逆散射成像方法，其特征在于：所述步骤三中，在计算成像阶段实现对目标内部的计算成像，即由下述非线性约束优化过程：

$\mathrm{\χ}=\underset{\mathrm{\χ}}{\arg }\min D\left(\mathrm{\χ}\right)$ s.t.E^scatt＝A(x)                            (6)

其中，

是最小化计算式，s.t.是“subject to”的缩写，表示同时满足，得到介电常数分布，作为目标内部的恢复图像。

6.根据权利要求1或5所述的电磁逆散射成像方法，其特征在于：所述步骤一中，对于单一频率测量，在步骤三优化过程(6)式的求解分为初值假设步骤，迭代优化步骤，和终止步骤：

1a、初值假设步骤，以已知背景介电常数分布作为目标分布初值，即在该步骤假设目标对比度为零；

1b、迭代优化步骤，设第n迭代步的计算结果为x_n，n＞1，为了简化(6)式的计算，(6)式中的D(x)用以下形式代替，即：

$D\left(\mathrm{\χ}\right)\←{\∫}_D\frac{{|\▿\mathrm{\χ}|}^2}{{({|\▿{\mathrm{\χ}}_{n-1}|}^2+\mathrm{\δ})}^{\frac{1}{2}}}{\mathrm{dr}}^{\′}---\left(7\right)$

其中，δ为人为引入的正值常数，取其为D_inv离散网格面积与工作波长平方的比值；(6)式中的A(x)用A(x)在x_n-1处的一阶泰勒展开代替，即：

$A\left(\mathrm{\χ}\right)\←A\left({\mathrm{\χ}}_{n-1}\right)+\▿A\·(\mathrm{\χ}-{\mathrm{\χ}}_{n-1})$

得到：

${\mathrm{\χ}}_n=\underset{\mathrm{\χ}}{\arg }\min \left|\right|\mathrm{W\χ}\left|\right|$ s.t. $y_n=\▿A\·\mathrm{\χ}---\left(9\right)$

其中，算子 $W=\frac{1}{\sqrt{{\left|{\▿\mathrm{\χ}}_{n-1}\right|}^2+\mathrm{\δ}}}\▿,$ 向量 $y_n=E^{\mathrm{scatt}}-A\left({\mathrm{\χ}}_{n-1}\right)+{\▿A}_{n-1}\·{\mathrm{\χ}}_{n-1},$ (8)式的解为：

${\mathrm{\χ}}_n={\left(W^{*}W\right)}^{-1}{(\▿A)}^{*}{(\▿A{\left(W^{*}W\right)}^{-1}{(\▿A)}^{*})}^{-1}{y}_n---\left(10\right)$

(9)式为迭代步骤的迭代表达式；

1c、终止步骤，当

时终止迭代，其中，η为人为设置的终止门限，取为η＝10^-3。

7.根据权利要求1所述的电磁逆散射成像方法，其特征在于：所述步骤一中，对于扫频测量，在步骤二中设目标介电常数在扫频带宽内缓慢变化，扫频测量以提高分辨率；

测量时，若有校准后的矢量网络分析仪，利用网络分析仪发射扫频连续波，直接得到扫频测量数据；若有校准后的多个发射、接收机组件，利用多个发射机发射扫频带宽内的不同正交编码脉冲信号，接收机接收、解码，通过时-频傅里叶变换得到扫频测量数据；

计算时，利用测得的扫频数据，首先计算低频率下的目标对比度分布，再以此分布为初值，计算高频率下的目标对比度分布，最后以高频率下的计算结果作为最终成像结果。

8.根据权利要求2所述的电磁逆散射成像方法，其特征在于：所述接收天线接收响应电场，包括幅度与相位信息。

9.根据权利要求1所述的电磁逆散射成像方法，其特征在于：本发明方法适用于不包含金属目标的场景成像：电感应测井和生物体软组织切片成像的应用。

Images