CN109239771B - 一种基于非均匀背景介质的弹性波成像方法 - Google Patents

Abstract

本发明公开一种基于非均匀背景介质的弹性波成像方法，以重建嵌入非均匀背景介质中的未知目标。将非均匀背景，即在有限域内将其作为已知散射体的壁面处理，其优点是避免了非均匀背景格林函数的耗时计算。在此方案下，结合改进的Levenberg Marquardt(LM)算法，提出了一种基于差分积分方程模型的新型方法(ME‑LM)以实现图像重建，另外，使用改进的广义交叉验证(GCV)正则化技术来自适应地选择迭代中的正则化参数，提高了迭代过程的稳定性。结果表明，所提出的方法不仅运算速度快，收敛性好，而且能够减轻非线性，能够重构出背景对比度高的未知散射体。仿真测试和实际测试均验证了该方法的有效性。

Description

一种基于非均匀背景介质的弹性波成像方法

技术领域

本发明属于弹性波成像技术领域，提出了一种基于非均匀背景介质的弹性波成像方法。

背景技术

在生产实践中，我们人类往往是通过自身的感应器官，如眼睛、鼻子、耳朵等，来直观地感受，从而获取外界的信息，比如说通过眼睛获取的图像，包含了巨大的信息量。然而对于某些领域，比如微观领域，海洋领域等，由于人类感官的局限性，单单通过人体感官，不足以获取较多有用的信息，仅仅通过感官获取的信息非常有限，于是人类就发明了各种工具来扩展人体感官的范围，其中的一些发明工具的方法是用来代替眼睛的，这一类方法就属于成像方法。成像方法的基本原理就是利用相应的发射装置或者设备，发射出特定的波比如弹性波，发射波通过传播媒介进行传播，达到探测目标后，会发生一定的散射和透射，然后我们再利用一定的接收装置或者设备，将这些散射波和透射波接收起来，从这些波中可以得到特定的信息，然后我们利用成像算法，通过计算机技术对这些信息进行加工处理，就可以形成图像。

在弹性波成像方法的发展过程中，在早期，主要采用的是线性的算法，但是它的局限性比较大，只能得到一些简单的信息，比如合成孔径雷达技术，只能得到物体的大概位置和形状，而无法获得物体材料参数等相关信息。成像的关键在于解决逆问题，然而，在复杂的情况下，解决逆问题的主要困难就是它的病态性和非线性，所谓病态性，是指方程无解或者解不唯一的情况。为了较好的解决这些困难，我们考虑采用非线性算法。本发明是基于改进增强的LevenbergMarquardt算法实现非均匀背景介质弹性波成像的方法(简称ME-LM)，能够显著提高成像的质量，可广泛应用于雷达成像和医学生物成像等领域。

在弹性波的成像方法中，根据等效感应流是否存在进行分类，方法分为场类型和源类型两种方法。波恩迭代法和变形波恩迭代法是典型的场型成像方法，ME-LM方法也是一种基于高斯-牛顿迭代法的场型成像方法。对比源反演方法和基于子空间优化的方法都是典型的源类型优化方法。由于弹性波的物理特性，弹性波成像已经被广泛应用到无损检测、医学造影、地质勘探等实际工程领域。本发明主要是针对弹性波成像中场类型方法设计的一种，现有的非均匀背景下的介质成像没有将已知背景和背景里面的未知散射体分开的方法，这样在迭代的过程中，容易使已知的非均匀背景受到噪声的污染而有误差产生，而与现有的非均匀背景下的成像方法相比，本发明方法采用将背景与背景里面的未知散射体分开来的方法，通过在迭代过程中仅仅只需反演未知散射体，而不会影响已知的非均匀背景，使迭代中的误差更小，从而可以更好地改善图像的质量，这样通过图像就可以获取更多有用的信息，比如可以获得物体的位置、形状、材料信息等。

发明内容

本发明的目的是针对目前成像算法对于非均匀背景介质成像的局限性，提出了一种基于二维非均匀背景介质弹性波的成像方法(ME-LM)。将非均匀背景，即在有限域内将其作为已知散射体的壁面处理，其优点是避免了非均匀背景格林函数的耗时计算。在此方案下，结合改进的Levenberg Marquardt(LM)算法，提出了一种基于差分积分方程模型的新型方法(ME-LM)以实现图像重建，另外，在该方法中，改进的广义交叉验证方法(GCV)并结合截断奇异值分解方法被用于自适应地选择正则化参数，使整个算法稳定地进行迭代反演，能够较好地反演出未知散射体的相关参数，改善重建图像的质量。

本发明的技术方案：

本发明设计方法利用接收装置获取到弹性波散射场数据后，利用弹性波场积分方程构建非均匀背景介质下相关目标函数方程和成本函数方程，通过发明的ME-LM方法通过迭代的求解方式，求出成本函数最小值，同时在迭代的过程中，为了防止迭代的方向出现大的偏差，我们采用改进的广义交叉验证方法来自适应地选择正则化参数，以提高算法的优化性能，具体如下：

1、本发明一种基于非均匀背景介质弹性波成像方法，弹性波为电磁波，包括以下：

步骤(1)、根据离散的网格位置和发射装置、接收装置的位置，计算格林函数数

和

并根据非均匀背景的对比度

和场积分方程(4)-(6)计算出相应的非均匀背景总场场强

散射场场强

和对比源

其中

为离散后的格林函数G(r_s,r′)的积分算子；

是离散后的格林函数G(r,r′)的积分算子；

表示一个位于空间r_s处的点源对其周围空间某一点r′所产生的场；

表示一个位于空间r处的点源对其周围空间某一点r′所产生的场；

为第一类零阶汉克尔函数，i表示虚数，k₀是弹性波的波数；具体是：

总场积分方程：

其中表示位于r处的入射场场强；χ(r′)＝(∈(r′)-∈₀)/∈₀，它为∈_r的对比度函数，∈₀表示弹性波穿过的介质的某种物理特性；L为发射装置的个数；

散射场积分方程：

表示的是位于r_s处的接收装置所接收到的散射体产生的散射场数据，

M为接收装置的个数；

对比源为对比度和总场的乘积，定义为：

I_l(r)＝χ(r)E_l(r) (3)

将公式(1)-(3)离散化得到：

总场场强：

散射场场强：

对比源：

其中(m,n)代表离散网格的中心坐标，

为第(m,n)个离散网格的感应电流，

为各个离散网格感应电流

的集合，

是离散后的格林函数G(r,r′)的积分算子，

为离散后的格林函数G(r_s,r′)的积分算子，

是χ(r′)的对角矩阵形式；

由于非均匀背景介质的对比度

已知,因此进一步可以得到：

非均匀背景总场场强

非均匀背景散射场场强

非均匀背景对比源

步骤(2)、初始化未知散射体对比度参数

和正则化参数α，并对它们赋值初值为0，同时设定迭代次数p＝0；

步骤(3)、将未知散射体对比度

代入到公式(10)和(11)中，得到未知散射体的对比源

和未知散射场理论值F(Δχ)，并求取此次迭代获得的未知散射场理论值F(Δχ)和未知散射体对比度

之间的雅可比矩阵D，并进行奇异值分解；具体是：

由于非均匀背景介质的对比度

已知，因此可将探测区域的对比度

和总场

以及对比源

分为已知的非均匀背景介质和背景介质里面的未知散射体这两部分：

其中

分别表示非均匀背景的对比度、总场、对比源,分别表示未知散射体的对比度、总场、对比源；

将公式(7)、(8)、(9)代入公式(4)、(6)则得到关于未知散射体对比源

的对比源积分方程：

将公式(10)代入到公式(5)中，则可以得到关于未知散射体的散射场场强理论值，也即目标函数：

其中⊙的定义为两个矩阵的乘法，vec{}定义为向量化张量的操作，公式(11)表明目标函数仅仅只是关于未知散射体的对比度

的函数；

根据

e^s代表表示已知非均匀背景的散射场场强数据和未知散射体的散射场测量数据之和，即接收装置直接接收到的数据，

代表的是表示已知非均匀背景的散射场场强数据，由于

也是已知，故而得到背景里面的未知散射体的散射场实测数据：

根据上述构建未知散射体的对比度

的成本函数：

f(Δχ)＝||F(Δχ)-Δe^s||² _min (12)

由于公式(12)的病态性，利用ME-LM方法进行优化，可以得到以下方程：

[D^*D+αI]Δ(Δχ_α)＝D^*δe^s (13)

其中α表示正则化参数，I是单位矩阵；δe^s＝Δe^s-F(Δχ_c)表示测量的未知散射体的散射场测量数据与利用公式(11)计算的未知散射体的散射场场强理论值之间的差值；其中Δχ_c表示本次迭代的数值结果；

进而可以得到：

Δ(Δχ_α)＝[D^*D+αI]^-1D^*δe^s (14)

由于公式(14)计算量比较大，并且不容易滤除数据噪声，故而采用奇异值分解的方法来求Δ(Δχ_α)；所述的奇异值分解具体是:

采用雅可比矩阵D＝U∑V^*求解方程(14)，得到前后相邻两次迭代未知散射体的对比度误差：

进一步获取下一次迭代的对比度参数：

Δχ_p+1＝Δχ_p+Δ(Δχ_α)_p (16)

其中V是大小为N×N的单位矩阵，表示雅可比矩阵D的右奇异矩阵；Σ为LM×N的对角矩阵，[∑]_k＝σ_k，k＝1,2,...,min(LM,N)，σ_k表示的是Σ第k行k列上的对角矩阵元素；

步骤(4)、利用公式(11)重新计算下一次迭次的未知散射体的散射场数据F(Δχ_p+1)，并求出未知散射体的散射场场强的理论值和测量值之间的差值δe^s，然后判断是否满足迭代停止条件(一般设置为||δe^s||＜0.001)，如果满足条件则结束并输出未知散射体最佳的对比度值

进而重建相应的图像，如果不满足则进行步骤(5)，继续进行迭代优化；

步骤5：更新迭代次数p＝p+1,利用公式(17)函数V(α)最小值时所对应的α值计算新的正则化参数α_p，返回步骤(3)继续进行优化；

改进的广义交叉函数如下所示：

其中，U是大小为LM×LM的单位矩阵，表示雅可比矩阵D的左奇异矩阵；n是矩阵U的行数，q是矩阵U的列数；K为雅可比矩阵D奇异值截取的个数。

2、本发明一种基于非均匀背景介质弹性波成像方法，弹性波为声波，包括以下：

步骤1：根据离散的网格位置和发射装置、接收装置的位置，计算格林函数

和

并根据非均匀背景的对比度

和场积分方程(20)-(21)计算出相应的非均匀背景总场场强P^b，散射场场强P^b,sca；

其中

为离散后的格林函数g(r_s,r′)的积分算子，

是离散后的格林函数g(r,r′)的积分算子，

表示一个位于空间r_s处的点源对其周围空间某一点r′所产生的场；

表示一个位于空间r处的点源对其周围空间某一点r′所产生的场；

为第一类零阶汉克尔函数，i表示虚数，k₀是弹性波的波数；具体是：

总场积分方程：

散射场积分方程：

将以上两个公式离散化之后可得：

其中

均为密度的对比度函数值，以

的形式表示χ的离散化形式，其它公式字母也均以此方式表示，P^inc、P、P^sca分别表示入射场、总场和散射场，由于非均匀背景介质的对比度

已知,因此进一步可以得到：

非均匀背景总场场强P^b：

非均匀背景散射场场强：P^b,sca：

步骤2：初始化未知散射体对比度参数Δχ₁，Δχ₂和正则化参数α，并对它们赋值初值为0，同时设定迭代次数p＝0；

步骤3：将Δχ₁，Δχ₂代入到公式(27)和(28)中，得到F(Δχ₁,Δχ₂)，并求取此次迭代获得的未知散射场数据

和未知散射体对比度

之间的雅可比矩阵D，并进行奇异值分解；具体是：

由于非均匀背景介质的对比度

已知，因此可将探测区域的对比度

和总场P分为已知的非均匀背景介质和背景里面的未知散射体这两部分：

P＝P^b+ΔP (26)

其中

表示非均匀背景的对比度，P^b表示非均匀背景的总场；

表示非均匀背景里面的未知散射体的对比度，ΔP表示非均匀背景里面的未知散射体的总场；将公式(24),(25),(26)代入公式(20),(22)中，并相减，得到关于

的场差分方程：

消去中间变量ΔP，并代入公式(21)得到关于未知散射体散射场的理论值，也即目标函数：

公式(28)表明目标函数仅仅只是关于未知散射体的对比度

的函数，

根据Δe^s＝e^s-P^b,sca，e^s表示已知的非均匀背景的散射场数据和未知散射体的散射场的测量数据之和，也即接收装置直接接收到的数据，P^b,sca表示已知的非均匀背景散射体的散射场测量数据，由于P^b,sca也是已知，故而得到背景里面的未知散射体的实测数据：

根据上述构建未知散射体的对比度的成本函数，令Δχ为Δχ₁,Δχ₂的集合，则：Δχ＝[Δχ₁,Δχ₂]，则F(Δχ₁,Δχ₂)＝F(Δχ)，构建成本函数为：

f(Δχ)＝||F(Δχ)-Δe^s||² _min (29)

由于公式(29)的病态性，因此利用ME-LM方法来优化，可以得到以下方程：

[D^*D+αI]Δ(Δχ_α)＝D^*δe^s (30)

进而可以得到：

Δ(Δχ_α)＝[D^*D+αI]^-1D^*δe^s (31)

由于公式(31)计算量比较大，并且不容易滤除数据噪声，故而采用奇异值分解的方法来求Δ(Δχ_α)，所述的奇异值分解具体是:

采用雅可比矩阵D＝U∑V^*求解方程(30)，得到前后相邻两次迭代未知散射体的对比度误差：

进一步获取下一次迭代的对比度参数：

Δχ_p+1＝Δχ_p+Δ(Δχ_α)_p (33)

其中U是大小为LM×LM的单位矩阵，它表示雅可比矩阵D的左奇异矩阵；V是大小为N×N的单位矩阵，它表示雅可比矩阵D的右奇异矩阵；Σ为LM×N的对角矩阵，[∑]_k＝σ_k，k＝1,2,...,min(LM,N)，σ_k表示的是Σ第k行k列上的对角矩阵元素，

步骤4：利用公式(28)重新计算下一次迭次的未知散射体的散射场数据F(Δχ_p+1)，并求出散射场的理论值和测量值之间的差值δe^s，判断是否满足迭代停止条件(一般设置为||δe^s||＜0.001)，如果满足条件则结束并输出未知散射体最佳的对比度值

进而重建相应的图像，如果不满足，则进行步骤5，继续进行迭代优化；

步骤5：更新迭代次数p＝p+1,利用公式(34)计算新的正则化参数α_p，返回步骤3继续进行优化；

改进的广义交叉函数如下所示：

其中，n是矩阵U的行数，q是矩阵U的列数；K为雅可比矩阵D奇异值截取的个数；当函数V(α)取得最小值时的α值，即为所要选取的正则化参数。

本发明的有益效果是：

本发明主要设计了在非均匀背景介质下基于ME-LM方法实现对二维弹性波散射数据成像的一种新型方法。自适应选择参数的GCV正则化技术可以确保方法在迭代过程中解的稳定性和精确性，并且有利于提高迭代的速度，使所构建出来的图像的质量更高，在实际验证过程中，通过对实测实验数据的重构结果表明：所提出的成像方法在弹性波成像中具有很好的有效性和准确性。

附图说明

图1是所提出的弹性波成像方法的迭代优化流程图；

图2是弹性波成像的实验测量装置结构图；

图3是弹性波成像方法对探测目标仿真数据的重建结果；

图4是弹性波成像方法对探测目标实测数据的重建结果。

具体实施方式

弹性波包括电磁波、声波等波，下面以电磁波为例，结合附图对本发明的成像方法作进一步说明。

弹性波都包含有两个重要的积分方程，总场积分方程和散射场积分方程，则电磁波的两个积分方程为：

总场积分方程：

其中G(r,r′)是二维自由空间格林函数

它的物理意义表示一个位于空间r处的点源对其周围空间某一点r′所产生的场，括号里面的r′代表场点的位置，r代表的是源点的位置,

是数学上的一类特殊函数，为第一类零阶汉克尔函数，i表示的是虚数的含义，

代表了位于r处的入射场数据，k₀是弹性波的波数，χ(r′)＝(∈(r′)-∈₀)/∈₀为对比度函数，∈₀表示弹性波穿过的介质的某种物理特性，如密度，介电常数等，本例中代表介质的介电常数，可见χ(r′)是介电常数的对比度函数值。L为发射装置的个数。

散射场积分方程：

表示的是位于r_s处的接收装置或者设备所接受到的散射体产生的散射场数据，

也为自由空间的格林函数，其物理意义与前面的G(r,r′)的物理意义相同，只是具体的空间位置不同。M为接收装置的个数。

对比源为对比度和总场的乘积，定义为：

I_l(r)＝χ(r)E_l(r) (3)

为了便于后续计算，将以上公式整合并离散化可得：

(m,n)代表离散网格的中心坐标，

为第(m,n)个离散网格的感应电流，

为各个离散网格感应电流

的集合，

是离散后的格林函数G(r,r′)的积分算子，

为离散后的格林函数G(r_s,r′)的积分算子，

是χ(r′)的对角矩阵形式。离散后的网格总数为N，由于非均匀背景介质的介电常数已知，因此可将探测区域的对比度

和总场

以及对比源

分为已知的非均匀背景介质和背景里面的未知散射体这两部分：

其中

表示非均匀背景的对比度，总场和对比源,

表示非均匀背景里面的未知散射体的对比度，总场，对比源。将公式(7)，(8)，(9)代入公式(4)，(5)则我们可以得到关于

的对比源积分方程：

将公式(10)代入到公式(6)中，则可以得到关于未知散射体散射场的计算值，也即目标函数：

其中⊙的定义为两个矩阵的乘法，vec{}定义为向量化张量的操作，公式(11)表明目标函数仅仅只是关于未知散射体的对比度

的函数，

因为

e^s代表已知的非均匀背景的散射场数据和未知散射体的散射场的测量数据之和，也即接收装置或者设备直接接收到的数据，

代表的是已知的非均匀背景散射体的散射场测量数据，

也是已知的，因此Δe^s也是已知的，故而得到背景里面的未知散射体的实测数据

根据上述构建未知散射体的对比度

的成本函数：

f(Δχ)＝||F(Δχ)-Δe^s||² _min (12)

由于公式(12)的病态性，也即解的不唯一性或者无解性，比较严重，因此我们利用ME-LM方法来优化解决方案，可以得到以下方程：

[D^*D+αI]Δ(Δχ_α)＝D^*δe^s (13)

进而可以得到：

Δ(Δχ_α)＝[D^*D+αI]^-1D^*δe^s (14)

其中D为雅可比矩阵，表示矩阵D的第i行第j列的元素取值，δe^s＝Δe^s-F(Δχ_c)表示测量的未知散射体的散射场数据与利用公式(11)计算的未知散射体的散射场数据理论值之间的差值,Δ(Δχ)＝Δχ₊-Δχ_c表示迭代过程中，下一次迭代求出的未知散射体的对比度Δχ₊与本次迭代求出的未知散射体的对比度Δχ_c的差值,α是一个正则化参数，它是利用广义交叉验证算法自适应地选择，I是单位矩阵，为了便于后续处理噪声影响问题，采用奇异值分解的方法D＝U∑V^*来求方程(14)，其中U为大小是LM×LM的单位矩阵，它是雅可比矩阵D的左奇异矩阵，V为N×N的单位矩阵它是雅可比矩阵D的右奇异矩阵，Σ为LM×N的对角矩阵，[∑]_k＝σ_k，k＝1,2,...,min(LM,N)，σ_k表示的是Σ第k行k列上的对角矩阵元素，α是正则化参数，则我们可得:

对于第p次迭代，可以得到更新的对比度参数：

Δχ_p+1＝Δχ_p+Δ(Δχ_α)_p (16)

具体的迭代流程如图1所示，本发明设计的基于非均匀背景介质的弹性波成像方法的具体实施方案包括以下步骤：

步骤1：根据离散的网格位置和发射、接收装置或者设备的位置，计算出格林函数

和

并根据非均匀背景的对比度

和场积分方程(4)，(5)，(6)计算出相应的非均匀背景场强

和对比源

步骤2：初始化未知散射体对比度参数

和正则化参数α，对它们赋值初值为0，设定p＝0；

步骤3：将代入到公式(10)和(11)中，得到

和F(Δχ)，并求取此次迭代获得的未知散射场数据F(Δχ)和未知散射体对比度

之间的雅可比矩阵D，并进行奇异值分解。

步骤4：p＝p+1,利用公式(17)计算新的正则化参数α_p，并代入公式(14)得到新的Δ(Δχ)_p，代入公式(16)即可得到第p+1次的对比度Δχ_p+1。

步骤5：利用公式(11)重新计算新的未知散射体的散射场理论值F(Δχ_p)，并求出理论值和测量值之间的差值δe^s，判断是否满足迭代停止条件，如果满足停止迭代，如果不满足，则返回步骤3继续进行迭代优化。优化停止后我们将会得到未知散射体最佳的对比度值

进而重建相应的图像。

在以上迭代地过程中，我们还运用正则化技术来求得正则化参数α，以提高迭代的稳定性和解的精确性，在实际过程中，通过最小化改进的广义交叉验证函数来自适应地选择正则化参数，改进的广义交叉验证函数如下所示：

其中，n是矩阵U的行数，q是矩阵U的列数。K为雅可比矩阵D奇异值截取的个数。当函数V(α)取得最小值时的α值，即为所要选取的正则化参数。

接下来考虑当弹性波为声波的情况，列出声波的场积分方程

总场积分方程：

散射场积分方程：

将以上两个公式离散化之后可得：

其中

均为密度的对比度函数值，在整篇文章中我们以

的形式表示χ的离散化形式，其它公式字母也均以此方式表示，P^inc，P^sca分别表示入射场和散射场，方程中其它各个参数的含义与电磁波的情况均一致，接下来我们由以上方程可推导关于非均匀背景介质和未知散射体之间的差分方程。

令

当感兴趣区域只存在非均匀背景介质,假设其对比度函数值为

且已知，则有关的场积分方程为：

如果感兴趣区域不确定是否存在未知散射体，那么将(20)式与(22)式相减可得：

则由以上公式消去中间变量，并代入公式(21)即可得到目标函数：

从方程中不难发现，声波的场积分方程和电磁波的积分方程很相似，因此接下来的整个迭代优化过程也和电磁波的推导过程一致，这里就不具体展开推导。接下来，我们以具体实施例子来论证所提出的技术。

实施例1.

本发明设计采用的实验装置结构图如图2所示，本例采用实验仿真数据验证所提出的弹性波成像方法，仿真例子是一个由矩形框架和圆环组成的轮廓。该轮廓位于边长为2λ的方形区域，其中λ是入射波的波长。不均匀的背景由矩形框架和自由空间组成。环和矩形框的中心都位于(0，0)处。环的外半径是0.4λ，内半径是0.2λ。矩形框的外壁长度为1.6λ，内壁长度为1.2λ。作为非均匀背景，矩形框的相对介电常数为2.0是已知的，环的相对介电常数对比度为4.0。该域被离散化为30×30网格，并且使用均匀分布在2π立体角上的16个入射点作为发射器，并且以感兴趣区域中心为圆心均匀分布有16个接收器用于测量散射数据。在所有测量数据中，均加入了20dB的加性高斯白噪声。重建效果如图3所示，从图中可以看出：ME-LM方法可以比其他两种方法更好地重建图像。传统的LM方法重建图像的效果很差。对于传统的DBIM方法，尽管其图像重建效果优于LM方法，但其效果仍不如ME-LM方法。从重构后的图像中可以看出，采用ME-LM方法重建的图像虽然背景中仍然存在一些波动，但与其他两种方法相比，目标区域内的矩形框和圆环都可以很好地重建。结果令人满意。

实施例2.

为了验证本发明设计的弹性波成像方法对实测数据的成像效果，本例2依然采用例1中的实验装置结构图，并搭建了相应的实际实验测量装置，我们对两种模型进行了测量，并使用校准后的实验数据进行成像。两种模型的矩形框墙壁都是由特氟隆材料制成的。墙是方形结构。它的边长为21cm，厚度为1cm。第一个模型是由两个U型的散射体背对背构成的K型散射体，U型散射体具有大约为3的相对介电常数，下部U型散射体具有7cm的边长，上部的U型结构具有10cm的边长。它们的厚度是1cm。第二个模型是圆柱体散射体，圆柱体的半径约为5cm，相对介电常数约为3，位置在(-2.5，-7.5)cm处，成像结果如图4所示。从我们可以看到，对于K型散射体，算法对中间一段介电常数的估计比对其他位置的估计准确，因为中间部分是两个U形散射体的接触部分，其实际厚度相当于2cm，这表明，当厚度增加时，算法对散射体的介电常数的估计变得更准确；对于圆柱结构，该方法不仅可以实现较好的壁面成像，而且可以清楚地区分圆柱形散射体与壁面之间的间隙。这表明我们提出的成像方法具有很高的可行性。

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

Claims (2)

1.一种基于非均匀背景介质弹性波成像方法，弹性波为电磁波，其特征在于包括以下：

步骤(1)、根据离散的网格位置和发射装置、接收装置的位置，计算

和并根据非均匀背景的对比度和场积分方程(4)-(6)计算出相应的非均匀背景总场场强

散射场场强和对比源

其中

为离散后的格林函数G(r_s，r′)的积分算子；

是离散后的格林函数G(r，r′)的积分算子；

表示一个位于空间r_s处的点源对其周围空间某一点r′所产生的场；

表示一个位于空间r处的点源对其周围空间某一点r′所产生的场；

为第一类零阶汉克尔函数，i表示虚数，k₀是弹性波的波数；具体是：

总场积分方程：

其中表示位于r处的入射场场强；χ(r′)＝(∈(r′)-∈₀)/∈₀，为∈_r的对比度函数，∈₀表示弹性波穿过的介质的某种物理特性；L为发射装置的个数；

散射场积分方程：

表示的是位于r_s处的接收装置所接收到的散射体产生的散射场数据，

M为接收装置的个数；

对比源为对比度和总场的乘积，定义为：

I_l(r)＝χ(r)E_l(r) (3)

将公式(1)-(3)离散化得到：

第(m，n)个离散网格的总场场强：

散射场场强：

对比源：

其中(m，n)代表离散网格的中心坐标，

为第(m，n)个离散网格的感应电流，

为各个离散网格感应电流

的集合，

是离散后的格林函数G(r，r′)的积分算子，为离散后的格林函数G(r_s，r′)的积分算子，

是χ(r′)的对角矩阵形式；

由于非均匀背景介质的对比度

已知，因此进一步可以得到：

非均匀背景总场场强

非均匀背景散射场场强

非均匀背景对比源

步骤(2)、初始化未知散射体对比度

和正则化参数α，并对它们赋值初值为0，同时设定迭代次数p＝0；

步骤(3)、将未知散射体对比度

代入到公式(10)和(11)中，得到未知散射体的对比源

和未知散射场理论值F(Δχ)，并求取此次迭代获得的未知散射场理论值F(Δχ)和未知散射体对比度

之间的雅可比矩阵D，并进行奇异值分解；具体是：

由于非均匀背景介质的对比度

已知，因此可将探测区域的对比度

和总场

以及对比源

分为已知的非均匀背景介质和背景介质里面的未知散射体这两部分：

其中

分别表示非均匀背景的对比度、总场、对比源，

分别表示未知散射体的对比度、总场、对比源；

将公式(7)、(8)、(9)代入公式(4)、(6)则得到关于未知散射体对比源

的对比源积分方程：

将公式(10)代入到公式(5)中，则可以得到关于未知散射体的散射场场强理论值，也即目标函数：

其中⊙的定义为两个矩阵的乘法，vec{}定义为向量化张量的操作，公式(11)表明目标函数仅仅只是关于未知散射体的对比度的函数；

根据

e^s表示已知非均匀背景的散射场场强数据和未知散射体的散射场测量数据之和，即接收装置直接接收到的数据，表示已知非均匀背景的散射场场强数据，由于

也是已知，故而得到背景里面的未知散射体的散射场实测数据：

根据上述构建未知散射体的对比度

的成本函数：

f(Δχ)＝||F(Δχ)-Δe^s||² _min (12)

由于公式(12)的病态性，利用改进增强的Levenberg-Marquardt方法进行优化，可以得到以下方程：

[D*D+αI]Δ(Δχ_α)＝D*δe^s (13)

其中α表示正则化参数，I是单位矩阵；δe^s＝Δe^s-F(Δχ_c)表示测量的未知散射体的散射场测量数据与利用公式(11)计算的未知散射体的散射场场强理论值之间的差值；其中Δχ_c表示本次迭代的数值结果；

进而可以得到：

Δ(Δχ_α)＝[D*D+αI]^-1D*δe^s (14)

由于公式(14)计算量比较大，并且不容易滤除数据噪声，故而采用奇异值分解的方法来求Δ(Δχ_α)；所述的奇异值分解具体是：

采用雅可比矩阵D＝U∑V^*求解方程(14)，得到前后相邻两次迭代未知散射体的对比度误差：

进一步获取下一次迭代的对比度参数：

Δχ_p+1＝Δχ_p+Δ(Δχ_α)_p (16)

其中V是大小为N×N的单位矩阵，表示雅可比矩阵D的右奇异矩阵；∑为LM×N的对角矩阵，[∑]_k＝σ_k，k＝1，2，...，min(LM，N)，σ_k表示的是∑第k行k列上的对角矩阵元素；

步骤(4)、利用公式(11)重新计算下一次迭次的未知散射体的散射场数据F(Δχ_p+1)，并求出未知散射体的散射场场强的理论值和测量值之间的差值δe^s，然后判断是否满足迭代停止条件，如果满足条件则结束并输出未知散射体最佳的对比度值

进而重建相应的图像，如果不满足则进行步骤(5)，继续进行迭代优化；

步骤(5)、更新迭代次数p＝p+1，利用公式(17)函数V(α)最小值时所对应的α值计算新的正则化参数α_p，返回步骤(3)继续进行优化；

改进的广义交叉函数如下所示：

其中，U是大小为LM×LM的单位矩阵，表示雅可比矩阵D的左奇异矩阵；n是矩阵U的行数，q是矩阵U的列数；K为雅可比矩阵D奇异值截取的个数。

2.一种基于非均匀背景介质弹性波成像方法，弹性波为声波，其特征在于包括以下：

步骤(1)、根据离散的网格位置和发射装置、接收装置的位置，计算

和并根据非均匀背景的对比度

和场积分方程(20)-(21)计算出相应的非均匀背景总场场强P^b，散射场场强P^b，sca；

其中

为离散后的格林函数g(r_s，r′)的积分算子，

是离散后的格林函数g(r，r′)的积分算子，

表示一个位于空间r_s处的点源对其周围空间某一点r′所产生的场；

表示一个位于空间r处的点源对其周围空间某一点r′所产生的场；

为第一类零阶汉克尔函数，i表示虚数，k₀是弹性波的波数；具体是：

总场积分方程：

散射场积分方程：

将以上两个公式离散化之后可得：

其中

χ₂均为密度的对比度函数值，以

的形式表示χ的离散化形式，其它公式字母也均以此方式表示，P^inc、P、P^sca分别表示入射场、总场和散射场，由于非均匀背景介质的对比度已知，因此进一步可以得到：

非均匀背景总场场强P^b：

非均匀背景散射场场强：P^b，sca：

步骤(2)、初始化未知散射体对比度Δχ₁，Δχ₂和正则化参数α，并对它们赋值初值为0，同时设定迭代次数p＝0；

步骤(3)、将Δχ₁，Δχ₂代入到公式(27)和(28)中，得到F(Δχ₁，Δχ₂)，并求取此次迭代获得的未知散射场数据

和未知散射体对比度

之间的雅可比矩阵D，并进行奇异值分解；具体是：

由于非均匀背景介质的对比度

已知，因此可将探测区域的对比度和总场P分为已知的非均匀背景介质和背景里面的未知散射体这两部分：

P＝P^b+ΔP (26)

其中

表示非均匀背景的对比度，P^b表示非均匀背景的总场；

表示非均匀背景里面的未知散射体的对比度，ΔP表示非均匀背景里面的未知散射体的总场；将公式(24)，(25)，(26)代入公式(20)，(22)中，并相减，得到关于

的场差分方程：

消去中间变量ΔP，并代入公式(21)得到关于未知散射体散射场的理论值，也即目标函数：

公式(28)表明目标函数仅仅只是关于未知散射体的对比度

的函数，

根据Δe^s＝e^s-P^b，sca，e^s表示已知的非均匀背景的散射场数据和未知散射体的散射场的测量数据之和，也即接收装置直接接收到的数据，P^b，sca表示已知的非均匀背景散射体的散射场测量数据，由于P^b，sca也是已知，故而得到背景里面的未知散射体的实测数据：

根据上述构建未知散射体的对比度

的成本函数，令Δχ为Δχ₁，Δχ₂的集合，则：Δχ＝[Δχ₁，Δχ₂]，则F(Δχ₁，Δχ₂)＝F(Δχ)，构建成本函数为：

f(Δχ)＝||F(Δχ)-Δe^s||² _min (29)

由于公式(29)的病态性，因此利用改进增强的Levenberg-Marquardt方法来优化，得到以下方程：

[D*D+αI]Δ(Δχ_α)＝D*δe^s (30)

进而可以得到：

Δ(Δχ_α)＝[D*D+αI]^-1D*δe^s (31)

由于公式(31)计算量比较大，并且不容易滤除数据噪声，故而采用奇异值分解的方法来求Δ(Δχ_α)，所述的奇异值分解具体是：

采用雅可比矩阵D＝U∑V^*求解方程(30)，得到前后相邻两次迭代未知散射体的对比度误差：

进一步获取下一次迭代的对比度参数：

Δχ_p+1＝Δχ_p+Δ(Δχ_α)_p (33)

其中U是大小为LM×LM的单位矩阵，它表示雅可比矩阵D的左奇异矩阵；V是大小为N×N的单位矩阵，它表示雅可比矩阵D的右奇异矩阵；∑为LM×N的对角矩阵，[∑]_k＝σ_k，k＝1，2，...，min(LM，N)，σ_k表示的是∑第k行k列上的对角矩阵元素，

步骤(4)、利用公式(28)重新计算下一次迭次的未知散射体的散射场数据F(Δχ_p+1)，并求出散射场的理论值和测量值之间的差值δe^s，判断是否满足迭代停止条件，如果满足条件则结束并输出未知散射体最佳的对比度值

进而重建相应的图像，如果不满足，则进行步骤(5)，继续进行迭代优化；

步骤(5)、更新迭代次数p＝p+1，利用公式(34)计算新的正则化参数α_p，返回步骤(3)继续进行优化；

改进的广义交叉函数如下所示：

其中，n是矩阵U的行数，q是矩阵U的列数；K为雅可比矩阵D奇异值截取的个数；当函数V(α)取得最小值时的α值，即为所要选取的正则化参数。

Images