均匀背景下的三维快速互相关对比源电磁反演方法

1.本技术涉及电磁反演成像技术领域，特别是涉及一种均匀背景下的三维快速互相关对比源电磁反演方法。


背景技术：

2.电磁反演技术已广泛应用于雷达成像等许多领域，反演技术和研究的重点是计算效率，在电磁反演的计算中，关于散射场和梯度场的两类线性求解方程占据了反演算法中大部分的计算复杂度，因此对于两类线性方程组的求解计算过程的简化成为简化电磁反演计算的关键。
3.对于反演技术的计算效率问题，提出了一种互相关对比度源反演(cross-correlated contrast source inversion，cc-csi)方法，cc-csi是一种非线性迭代反演方法，在cc-csi方法中，状态误差和数据误差是互相关的，通过最小化互相关误差来稳定反演过程，相较于传统的对比源反演(csi)方法和乘法正则化csi(mr-csi)方法，cc-csi方法具有更高的反演精度和更好的鲁棒性。
4.在三维反演中，通过cc-csi方法求解散射场的第一类多线性方程组和梯度场的第二类多线性方程组，传统的求解算法为lu矩阵方法和bicg算法，lu矩阵算法和bicg算法是基于包含二阶中心差分近似误差的刚度矩阵，由于刚度矩阵的误差，从而影响了反演的精度，此外，两种算法在反演区域边界需要额外的完美匹配层(pmls)，造成了计算资源的浪费，计算复杂度也大大增加。


技术实现要素：

5.基于此，有必要针对上述技术问题，提供一种能够在三维反演中保证反演精度的前提下，提高反演计算效率的均匀背景下的三维快速互相关对比源电磁反演方法。
6.一种均匀背景下的三维快速互相关对比源电磁反演方法，所述方法包括：
7.在均匀背景下，基于三维快速互相关对比源电磁反演方法，构建反演求解模型；反演求解模型中包括：根据对比源矩阵和刚度矩阵计算散射场的第一类多线性方程组以及根据残差矩阵和共轭转置刚度矩阵计算梯度场的第二类多线性方程组；
8.获取对比源矩阵，对对比源矩阵进行三维傅里叶变换得到三维对比源空间谱矩阵，构建刚度矩阵对应的第一类核函数矩阵，对第一类核函数矩阵进行三维傅里叶变换得到第一类核函数三维空间谱矩阵；
9.对三维对比源空间谱矩阵和第一类核函数三维空间谱矩阵进行计算得到三维散射场空间谱矩阵，对三维散射场空间谱矩阵进行三维逆傅里叶变换得到散射场，完成第一类多线性方程组的求解；
10.获取残差矩阵，对残差矩阵进行三维傅里叶变换得到三维残差空间谱矩阵，构建共轭转置刚度矩阵对应的第二类核函数矩阵，对第二类核函数矩阵进行三维傅里叶变换得到第二类核函数三维空间谱矩阵；
11.对三维残差空间谱矩阵和第二类核函数三维空间谱矩阵进行计算得到三维梯度场空间谱矩阵，对三维梯度场空间谱矩阵进行三维逆傅里叶变换得到梯度场，完成第二类多线性方程组的求解；
12.根据第一类多线性方程组和第二类多线性方程组的求解完成反演求解模型的计算。
13.在其中一个实施例中，在均匀背景下，基于三维快速互相关对比源电磁反演方法，构建反演求解模型，包括：
14.反演求解模型中的两类多线性方程组表示为
15.ae＝j
16.ahg＝s
17.其中，ae＝j表示第一类多线性方程组，ahg＝s表示第二类多线性方程组，a表示刚度矩阵，e表示散射场，j＝χe
tot
表示对比源矩阵，χ表示对比度，e
tot
表示总场，ah表示共轭转置刚度矩阵，g表示梯度场，s表示残差矩阵。
18.在其中一个实施例中，获取对比源矩阵，对对比源矩阵进行三维傅里叶变换得到三维对比源空间谱矩阵，包括：
19.获取对比源矩阵函数jm(x)，m∈[1,2,3]；
[0020]
对对比源矩阵函数jm(x)进行三维傅里叶变换，得到三维对比源空间谱矩阵其中，表示三维空间谱的频率向量，x＝(x1,x2,x3)表示三维空间位置坐标向量。
[0021]
在其中一个实施例中，构建刚度矩阵对应的第一类核函数矩阵，对第一类核函数矩阵进行三维傅里叶变换得到第一类核函数三维空间谱矩阵，包括：
[0022]
构建刚度矩阵对应的第一类核函数矩阵分别表示为
[0023][0024]
[0025][0026]
其中，n,m∈[1,2,3]表示不同分量，ω表示角频率，i2＝-1，k表示不同频率的波数，表示三维空间位置坐标向量x到原点的距离；
[0027]
对第一类核函数矩阵进行三维傅里叶变换，得到第一类核函数三维空间谱矩阵其中，在均匀背景中保持不变。
[0028]
在其中一个实施例中，对三维对比源空间谱矩阵和第一类核函数三维空间谱矩阵进行计算得到三维散射场空间谱矩阵，包括：
[0029]
根据逐点乘法对三维对比源空间谱矩阵和第一类核函数三维空间谱矩阵进行计算，得到三维散射场空间谱矩阵表示为
[0030][0031]
在其中一个实施例中，对三维散射场空间谱矩阵进行三维逆傅里叶变换得到散射场，完成第一类多线性方程组的求解，包括：
[0032]
对三维散射场空间谱矩阵进行三维逆傅里叶变换，得到散射场的空间分布en(x)，表示为
[0033][0034]
其中，表示三维位置坐标空间；
[0035]
根据en(x)完成第一类多线性方程组的求解。
[0036]
在其中一个实施例中，获取残差矩阵，对残差矩阵进行三维傅里叶变换得到三维残差空间谱矩阵，包括：
[0037]
获取残差矩阵函数sm(y)，对残差矩阵函数sm(y)进行三维傅里叶变换，得到三维残差空间谱矩阵其中，表示反演域的三维空间谱的频率向量，y＝(y1,y2,y3)表示反演域的三维空间位置坐标向量。
[0038]
在其中一个实施例中，构建共轭转置刚度矩阵对应的第二类核函数矩阵，对第二类核函数矩阵进行三维傅里叶变换得到第二类核函数三维空间谱矩阵，包括：
[0039]
构建共轭转置刚度矩阵对应的第二类核函数矩阵其中，
表示共轭运算；
[0040]
对第二类核函数矩阵进行三维傅里叶变换得到第二类核函数三维空间谱矩阵其中，在均匀背景中保持不变。
[0041]
在其中一个实施例中，对三维残差空间谱矩阵和第二类核函数三维空间谱矩阵进行计算得到三维梯度场空间谱矩阵，包括：
[0042]
根据逐点乘法对三维残差空间谱矩阵和第二类核函数三维空间谱矩阵进行计算，得到三维梯度场空间谱矩阵表示为
[0043][0044]
在其中一个实施例中，对三维梯度场空间谱矩阵进行三维逆傅里叶变换得到震源产生的梯度场，完成第二类多线性方程组的求解，包括：
[0045]
对三维梯度场空间谱矩阵进行三维逆傅里叶变换，得到梯度场的空间分布gn(y)，表示为
[0046][0047]
根据gn(y)完成第二类多线性方程组的求解。
[0048]
上述均匀背景下的三维快速互相关对比源电磁反演方法，通过在均匀背景下，基于三维快速互相关对比源电磁反演方法，构建包括第一类多线性方程组和第二类多线性方程组的反演求解模型；通过对对比源矩阵和刚度矩阵对应的第一类核函数矩阵进行计算，得到三维散射场空间谱矩阵，对三维散射场空间谱矩阵进行三维逆傅里叶变换得到散射场，完成第一类多线性方程组的求解；通过对残差矩阵和共轭转置刚度矩阵对应的第二类核函数矩阵进行计算，得到三维梯度场空间谱矩阵，对三维梯度场空间谱矩阵进行三维逆傅里叶变换得到梯度场，完成第二类多线性方程组的求解；根据第一类多线性方程组和第二类多线性方程组的求解完成反演求解模型的计算。与现有技术相比，本发明基于三维快速互相关对比源电磁反演方法，对电磁反演成像中的两类多线性方程组进行快速求解，来完成反演求解模型的计算并实现快速电磁反演成像，降低了电磁反演成像技术的计算复杂度，提高了电磁反演的计算精度和计算速度，从而有效提升电磁反演算法在实际问题中的可用性。
附图说明
[0049]
图1为一个实施例中均匀背景下的三维快速互相关对比源电磁反演方法的流程示意图；
[0050]
图2为为一个实施例中在twospheres_pp(at 4ghz)数据集中均匀背景下的三维快速互相关对比源电磁反演方法的反演结果示意图：(a)为目标的真实形状示意图；(b)为目
标的反演形状示意图；(c)为在z＝0mm处的反演得到的相对介电常数示意图；(d)为在z＝0mm处的反演得到的电导率示意图；(e)为在y＝0mm处反演得到的相对介电常数示意图；(f)为在y＝0mm处反演得到的电导率示意图；(g)为在x＝0mm处反演得到的相对介电常数示意图；(h)为在x＝0mm处反演得到的电导率示意图。
具体实施方式
[0051]
为了使本技术的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本技术进行进一步详细说明。应当理解，此处描述的具体实施例仅仅用以解释本技术，并不用于限定本技术。
[0052]
在一个实施例中，提供了一种均匀背景下的三维快速互相关对比源电磁反演方法，包括以下步骤：
[0053]
步骤102，在均匀背景下，基于三维快速互相关对比源电磁反演方法，构建反演求解模型；反演求解模型中包括：根据对比源矩阵和刚度矩阵计算散射场的第一类多线性方程组以及根据残差矩阵和共轭转置刚度矩阵计算梯度场的第二类多线性方程组。
[0054]
可以理解，均匀背景是指在均匀背景介质下；第一类多线性方程组和第二类多线性方程组并不是分别只有一组，而是各有多组线性方程组。
[0055]
步骤104，获取对比源矩阵，对对比源矩阵进行三维傅里叶变换得到三维对比源空间谱矩阵，构建刚度矩阵对应的第一类核函数矩阵，对第一类核函数矩阵进行三维傅里叶变换得到第一类核函数三维空间谱矩阵。
[0056]
可以理解，通过对对比源空间谱矩阵和第一类核函数三维空间谱矩阵分别进行点乘和三维傅里叶变换，将散射场的计算转换到三维空间谱域上，简化了散射场的计算过程。
[0057]
步骤106，对三维对比源空间谱矩阵和第一类核函数三维空间谱矩阵进行计算得到三维散射场空间谱矩阵，对三维散射场空间谱矩阵进行三维逆傅里叶变换得到散射场，完成第一类多线性方程组的求解。
[0058]
可以理解，通过对三维散射场空间谱矩阵进行三维逆傅里叶变换，将散射场空间谱还原到空间维度。
[0059]
步骤108，获取残差矩阵，对残差矩阵进行三维傅里叶变换得到三维残差空间谱矩阵，构建共轭转置刚度矩阵对应的第二类核函数矩阵，对第二类核函数矩阵进行三维傅里叶变换得到第二类核函数三维空间谱矩阵。
[0060]
可以理解，通过对残差空间谱矩阵和第二类核函数三维空间谱矩阵分别进行点乘和三维傅里叶变换，将梯度场的计算转换到三维空间谱域上，简化了梯度场的计算过程。
[0061]
步骤110，对三维残差空间谱矩阵和第二类核函数三维空间谱矩阵进行计算得到三维梯度场空间谱矩阵，对三维梯度场空间谱矩阵进行三维逆傅里叶变换得到梯度场，完成第二类多线性方程组的求解。
[0062]
可以理解，通过对三维梯度场空间谱矩阵进行三维逆傅里叶变换，将梯度场还原到原来的空间维度。
[0063]
步骤112，根据第一类多线性方程组和第二类多线性方程组的求解完成反演求解模型的计算。
[0064]
可以理解，通过对两类多线性方程组的求解，完成反演求解模型的计算，从而实现
快速电磁反演成像。
[0065]
上述均匀背景下的三维快速互相关对比源电磁反演方法，通过在均匀背景下，基于三维快速互相关对比源电磁反演方法，构建包括第一类多线性方程组和第二类多线性方程组的反演求解模型；通过对对比源矩阵和刚度矩阵对应的第一类核函数矩阵进行计算，得到三维散射场空间谱矩阵，对三维散射场空间谱矩阵进行三维逆傅里叶变换得到散射场，完成第一类多线性方程组的求解；通过对残差矩阵和共轭转置刚度矩阵对应的第二类核函数矩阵进行计算，得到三维梯度场空间谱矩阵，对三维梯度场空间谱矩阵进行三维逆傅里叶变换得到梯度场，完成第二类多线性方程组的求解；根据第一类多线性方程组和第二类多线性方程组的求解完成反演求解模型的计算。与现有技术相比，本发明基于三维快速互相关对比源电磁反演方法，对电磁反演成像中的两类多线性方程组进行快速求解，来完成反演求解模型的计算并实现快速电磁反演成像，降低了电磁反演成像技术的计算复杂度，提高了电磁反演的计算精度和计算速度，从而有效提升电磁反演算法在实际问题中的可用性。
[0066]
在其中一个实施例中，在均匀的背景介质下，基于三维快速互相关对比源电磁反演方法，构建反演求解模型，其中反演求解模型包括散射场的第一类多线性方程组和梯度场的第二类多线性方程组，模型表示如下
[0067]
ae＝j
[0068]ah
g＝s
[0069]
其中，ae＝j表示第一类多线性方程组，ahg＝s表示第二类多线性方程组，表示频域有限差分(fdfd)方法中的刚度矩阵，e表示散射场，表示对比源矩阵，χ表示对比度，e
tot
表示总场，该矩阵中的每列均为有限差分(fd)模式中的矢量形式，ah表示共轭转置刚度矩阵，g表示梯度场，s表示残差矩阵，n
src
表示激励源数，n表示表示反演在每个维度上划分的网格数，3n3×
3n3表示刚度矩阵的维度，3n3×nsrc
表示对比源矩阵的维度。
[0070]
在其中一个实施例中，获取在x1、x2和x3分量上产生的对比源矩阵函数jm(x)，对对比源矩阵函数jm(x)做2n
×
2n
×
2n点的三维傅里叶变换，得到三维对比源空间谱矩阵其中表示三维空间谱的频率向量，x＝(x1,x2,x3)表示三维空间位置坐标向量，m∈[1,2,3]。
[0071]
在其中一个实施例中，根据对比源矩阵函数jm(x)激发的散射场构建刚度矩阵对应的第一类核函数矩阵对第一类核函数矩阵做2n
×
2n
×
2n点的三维傅里叶变换，得到第一类核函数三维空间谱矩阵，其中，在均匀背景中保持不变，n,m∈[1,2,3]；
[0072]
具体地，根据对比源矩阵函数jm(x)激发的散射场构建刚度矩阵对应的第一类核函数矩阵包括：
[0073]
根据对比源矩阵函数j1(x)在x1、x2和x3分量上激发的散射场构建表示为
[0074][0075]
根据对比源矩阵函数j2(x)在x1、x2和x3分量上激发的散射场构建表示为
[0076][0077]
根据对比源矩阵函数j3(x)在x1、x2和x3分量上激发的散射场构建表示为
[0078][0079]
其中，ω表示角频率，i2＝-1，k表示不同频率的波数，表示三维空间位置坐标向量x到原点的距离。
[0080]
可以理解，通过构建核函数，可以使方程求解的精度不受网格划分大小的影响。
[0081]
根据逐点乘法对三维对比源空间谱矩阵和第一类核函数三维空间谱矩阵进行计算，得到三维散射场空间谱矩阵表示为
[0082][0083]
在其中一个实施例中，对三维散射场空间谱矩阵做2n
×
2n
×
2n点三维逆傅里叶变换，得到散射场的空间分布en(x)，表示为
[0084][0085]
其中，表示三维位置坐标空间，具体地，当n＝1,2,3时，en(x)分别表示散射场的第1,4,7,
…
,3n
3-2行、第2,5,8,
…
,3n
3-1行和第3,6,9,
…
,3n3行；
[0086]
根据散射场的空间分布en(x)完成第一类多线性方程组的求解。
[0087]
在其中一个实施例中，获取在x1、x2和x3分量上产生的残差矩阵函数sm(y)，对残差矩阵函数sm(y)做2n
×
2n
×
2n点三维傅里叶变换，得到三维残差空间谱矩阵其中，表示反演域的三维空间谱的频率向量，y＝(y1,y2,y3)表示反演域的三维空间位置坐标向量。
[0088]
在其中一个实施例中，构建共轭转置刚度矩阵对应的第二类核函数矩阵其中，对第二类核函数矩阵做2n
×
2n
×
2n点的三维傅里叶变换得到第二类核函数三维空间谱矩阵其中，在均匀背景中保持不变，可以预先计算和存储以供重用，表示共轭运算。
[0089]
在其中一个实施例中，根据逐点乘法对三维残差空间谱矩阵和第二类核函数三维空间谱矩阵进行计算，得到三维梯度场空间谱矩阵，表示为
[0090][0091]
在其中一个实施例中，对三维梯度场空间谱矩阵做2n
×
2n
×
2n点三维逆傅里叶变换，得到梯度场的空间分布gn(y)，表示为
[0092][0093]
其中，当n＝1,2,3时，gn(y)分别表示梯度场的第1,4,7,
…
,3n
3-2行、第2,5,8,
…
,3n
3-1行和第3,6,9,
…
,3n3行。
[0094]
根据梯度场的空间分布gn(y)完成第二类多线性方程组的求解。
[0095]
为进一步说明本发明所提出的均匀背景下的三维快速互相关对比源电磁反演方法的有益效果，在twospheres_pp(at 4ghz)数据集上进行了实验验证，其中，twospheres_
pp(at 4ghz)数据集的目标由两个直径为50毫米、介电常数为26.的介电球组成，反演区域设置为[-64,64；-64,64；-64,64]mm3，网格尺寸是32
×
32
×
32。数据集中的目标被从81个入射角照射，接收天线被限制在距离中心r＝1.796m的方位角平面内，并且由于技术原因，它们不能比50
°
更接近发射源顶点，方位角范围为0
°
到350
°
，步长为10
°
，测量复数数据进行反演，入射场在坐标系原点处模拟为振幅为1，相位为0的平面波。
[0096]
反演结果如图2所示，图2中的(a)为目标的真实形状示意图，(b)为目标的反演形状示意图，(c)为在z＝0mm处的反演得到的相对介电常数示意图，(d)为在z＝0mm处的反演得到的电导率示意图，(e)为在y＝0mm处反演得到的相对介电常数示意图，(f)为在y＝0mm处反演得到的电导率示意图，(g)为在x＝0mm处反演得到的相对介电常数示意图，(h)为在x＝0mm处反演得到的电导率示意图。由图2可知，反演得到的目标的反演形状和相对介电常数的精度都在可接受的范围内。
[0097]
在一个具体的实施例中，还将本发明提出的均匀背景下的三维快速互相关对比源电磁反演方法在twospheres_pp、twocubes pp、cylinder pp、cubespheres pp以及myster pp这五种数据集上的运行时间进行了比较，如表1所示，表1中iteration number表示迭代次数，total time表示完整运行时间，表示平均运行时间，n
iter
表示迭代数，nf表示频率数。
[0098]
表1均匀背景下的三维快速互相关对比源电磁反演方法在五种数据集中的运行时间
[0099][0100]
由表1可知，本发明所提出的反演方法在几分钟内完成了反演过程，而在目前常规的反演方法中，使用基于lu矩阵算法的反演方法在同样的反演情况下则需要几个小时才能完成反演，使用基于gpu加速的bicg算法的反演方法在同样的反演情况下虽然比lu矩阵算法快近十倍，但本发明的提出的反演方法仍然比基于gpu加速的bicg算法快二十倍以上。相比于现有技术，本发明所提出均匀背景下的三维快速互相关对比源电磁反演方法的计算速
度更快，计算效率和计算精度更高。并且相较于现有技术，本发明在反演的区域边界不需要额外的完美匹配层，实现了计算资源的有效利用。
[0101]
应该理解的是，虽然图1流程图中的各个步骤按照箭头的指示依次显示，但是这些步骤并不是必然按照箭头指示的顺序依次执行。除非本文中有明确的说明，这些步骤的执行并没有严格的顺序限制，这些步骤可以以其它的顺序执行。而且，图1中的至少一部分步骤可以包括多个子步骤或者多个阶段，这些子步骤或者阶段并不必然是在同一时刻执行完成，而是可以在不同的时刻执行，这些子步骤或者阶段的执行顺序也不必然是依次进行，而是可以与其它步骤或者其它步骤的子步骤或者阶段的至少一部分轮流或者交替地执行。
[0102]
以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。
[0103]
以上所述实施例仅表达了本技术的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本技术构思的前提下，还可以做出若干变形和改进，这些都属于本技术的保护范围。因此，本技术专利的保护范围应以所附权利要求为准。