CN111766575A - 一种穿墙雷达自聚焦稀疏成像方法及计算机设备 - Google Patents

Abstract

本发明适用于穿墙雷达成像领域，提供了一种穿墙雷达自聚焦稀疏成像方法及计算机设备。方法包括：获取扩展目标的回波信号，并构建以墙体厚度和相对介电常数作为参数的参数化字典信号模型；将扩展目标的稀疏特性作为先验信息，基于全变分约束的最大后验概率估计和参数化字典信号模型，构建包含参数化字典的目标函数；根据目标函数交替迭代更新稀疏反射系数和墙体参数，其中，用哈希表和线性卷积对更新公式中包含的字典矩阵的相关运算进行替代；达到迭代终止条件时，输出外循环终止时对应的墙体厚度和相对介电常数作为墙体参数估计值，对应的稀疏反射系数用于成像。本发明有效地保留了扩展目标的边缘特性，在消除墙体参数未知引起的目标位置偏移的同时使成像结果更清晰，且有效的减少运算时间，并在空间复杂度方面得到了改善。

Description

一种穿墙雷达自聚焦稀疏成像方法及计算机设备

技术领域

本发明属于穿墙雷达成像领域，尤其涉及一种穿墙雷达自聚焦稀疏成像方法及计算机设备。

背景技术

穿墙雷达成像中，由于墙体的存在，当电磁波经过墙体时，会在墙体前后表面发生透射、反射以及信号衰减等现象。电磁波在墙体内的传播速度和路径与墙体参数有关，即墙体厚度、墙体相对介电常数。墙体参数的不确定性导致了时延补偿的不确定性，进而影响最终的成像质量。墙体时延补偿不当会引起目标定位出现偏差以及图像模糊等问题。现有的穿墙雷达成像算法大都是在墙体参数信息已知的前提下进行成像，而在实际的观测场景中，墙体参数信息无法提前准确获取。因此，如何在墙体参数未知的前提下实现目标精确重构是在穿墙成像中所需研究的问题之一。

针对上述问题，目前已经有大量研究致力于解决墙体参数未知导致的目标位置出现偏差的问题。例如，在CS(压缩感知)的架构下进行自聚焦稀疏重构，同时采用非线性共轭梯度法更新墙体参数，但是该方法只考虑了目标具有的稀疏性，而忽略了现实场景中目标多为扩展目标的事实。为了充分利用目标具有的结构稀疏特性，目前也有采用引入针板先验，然后通过交替迭代来求解墙体参数和目标反射系数，在更新墙体参数时，通过最小二乘估计墙体参数增量，取得了不错的效果。然而，在CS的架构下进行自聚焦稀疏重构时，不可避免的需要建立字典矩阵，并且当每次采用不同的墙体参数进行时延补偿实现成像时，都需要更新字典矩阵，这无疑会导致大的存储压力，同时多次字典求逆和乘法运算导致计算复杂度高。

发明内容

本发明的目的在于提供一种穿墙雷达自聚焦稀疏成像方法、计算机可读存储介质及计算机设备，旨在解决在CS的架构下进行自聚焦稀疏重构时，需要建立字典矩阵，并且随着墙体参数的更新，字典矩阵会发生变换，导致大的存储压力，同时多次字典求逆和乘法运算导致计算复杂度高的问题。

第一方面，本发明提供了一种穿墙雷达自聚焦稀疏成像方法，所述方法包括：

S101、获取扩展目标的回波信号，并构建以墙体厚度和相对介电常数作为参数的参数化字典信号模型；

S102、将扩展目标的稀疏特性作为先验信息，基于全变分约束的最大后验概率估计和所述参数化字典信号模型，构建包含参数化字典的目标函数；

S103、根据所述包含参数化字典的目标函数交替迭代更新稀疏反射系数和墙体参数，其中，用哈希表和线性卷积对更新公式中包含的字典矩阵的相关运算进行替代；

S104、达到迭代终止条件时，输出外循环终止时对应的墙体厚度和相对介电常数作为墙体参数估计值，对应的稀疏反射系数用于成像。

第二方面，本发明提供了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现如所述的穿墙雷达自聚焦稀疏成像方法的步骤。第三方面，本发明提供了一种计算机设备，包括：

一个或多个处理器；

存储器；以及

一个或多个计算机程序，所述处理器和所述存储器通过总线连接，其中所述一个或多个计算机程序被存储在所述存储器中，并且被配置成由所述一个或多个处理器执行，所述处理器执行所述计算机程序时实现如所述的穿墙雷达自聚焦稀疏成像方法的步骤。

在本发明中，由于穿墙雷达自聚焦稀疏成像过程中，将扩展目标的稀疏特性作为先验信息，基于全变分约束的最大后验概率估计和所述参数化字典信号模型，构建包含参数化字典的目标函数，然后通过交替迭代更新稀疏反射系数和墙体参数。从而在进行墙体参数估计的同时，有效地保留了扩展目标的边缘特性，在消除墙体参数未知引起的目标位置偏移的同时使成像结果更清晰。在此基础上，通过用哈希表和线性卷积对更新公式中包含的字典矩阵的相关运算进行替代，从而在整个自聚焦过程中无需构建、存储和计算字典矩阵，将迭代循环中繁杂的计算过程进行简化，有效的减少运算时间，并在空间复杂度方面得到了改善。

附图说明

图1是本发明实施例一提供的穿墙雷达自聚焦稀疏成像方法的流程图。

图2是仿真场景示意图，墙后一个扩展目标长0.5m，宽0.2m。

图3(a)是忽略墙体时采用延时求和算法直接进行成像的结果。

图3(b)是采用延时求和算法，在真实仿真参数d＝0.2,ε＝4.5下进行成像的结果。

图4(a)为k＝1时的成像结果，此时的墙体参数为d＝0.1734,ε＝5.4896。

图4(b)为迭代结束的成像结果，对应的墙体参数估计值为d＝0.2002,ε＝4.5136。

图5是本发明实施例三提供的计算机设备的具体结构框图。

具体实施方式

为了使本发明的目的、技术方案及有益效果更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

为了说明本发明所述的技术方案，下面通过具体实施例来进行说明。

实施例一：

请参阅图1，本发明实施例一提供的穿墙雷达自聚焦稀疏成像方法包括以下步骤：

S101、获取扩展目标的回波信号，并构建以墙体厚度和相对介电常数作为参数的参数化字典信号模型。

在本发明实施例一中，S101具体可以包括：

将成像区域离散化为N个像素网格点，堆叠对成像区域进行探测的L个收发共置天线的数据，并构建参数化字典信号模型，可表示为：

y＝A(θ)x+e (1)

式中，x是稀疏反射系数，x＝[x₁,x₂,...,x_N]^T，x∈C^N×1(C^N×1表示N×1维的向量)，N为总像素网格点数；y为回波数据矢量，y∈C^ML×1，M为每个收发共置天线处接收回波的采样点数；e为噪声矢量，e∈C^ML×1；A(θ)为参数化字典矩阵，θ＝[d,ε]，d是墙体厚度，ε是墙体相对介电常数，A(θ)∈C^ML×N，且

A_l(θ)是用第l个收发共置天线构建的子矩阵；A_l(θ)的第m行第i列的元素为：

A_l(m,i)＝s(mT_s-τ_li) (2)

其中，s(mT_s)是发射信号序列，m＝1,2,...,M；若记第l个天线位置坐标为(x_p,y_p)，第i个像素点位置坐标为(x_a,y_a)，电磁波在自由空间的传播速度为c，第l个收发共置天线和第i个像素点之间的传输时延为τ_li，则有

其中，

S102、将扩展目标的稀疏特性作为先验信息，基于全变分(Total Variation，TV)约束的最大后验概率(Maximum a Posteriori，MAP)估计和所述参数化字典信号模型，构建包含参数化字典的目标函数。

因为实际场景中目标大多为扩展目标，因此采用包含TV约束项的MAP成像模型来重构目标反射系数矢量x，既保留了扩展目标边缘特性，又避免了L1范数正则化问题中惩罚参数需事先人工设定参数的问题。

在本发明实施例一中，S102具体可以包括：

假设噪声矢量e服从均值为0，噪声的协方差矩阵βI的复高斯分布，I为单位矩阵，通过将墙体参数引入到基于TV约束的MAP稀疏成像方法当中，构建包含参数化字典矩阵的目标函数，可表述为：

其中，D为一阶差分算子，

β是噪声参数。

S103、根据所述包含参数化字典的目标函数交替迭代更新稀疏反射系数和墙体参数(墙体厚度和墙体相对介电常数)，其中，用哈希表和线性卷积对更新公式中包含的字典矩阵的相关运算进行替代。

在本发明实施例一中，S103具体可以包括：

S1031、固定所述包含参数化字典的目标函数中的墙体参数，即d^(k),ε^(k)，θ^(k)＝[d^(k),ε^(k)]，更新稀疏反射系数x^(k+1)，并用哈希表和线性卷积对更新公式中包含字典矩阵的相关运算进行替代；其中k为当前外循环迭代次数，令外循环最大循环次数为K，初始状态下k＝0，墙体厚度为d⁽⁰⁾，墙体相对介电常数为ε⁽⁰⁾。

在本发明实施例一中，S1031具体可以包括：

固定θ^(k)时，A(θ^(k))为确定性矩阵，记A＝A(θ^(k))，则包含参数化字典的目标函数(3)为：

此时，通过交替迭代更新β和x，求解使得目标函数J₁最小的最优x。

具体过程如下：

S10311、令内循环最大循环次数为T，当前迭代次数为t，初始状态下t＝0，采用DAS(延时求和)得到稀疏反射系数的初值x⁽⁰⁾。

S10312、对式(4)中的目标函数J₁(x,β)进行优化最小化(Majorize-Minimize，MM)替代，然后对目标函数J₁(x,β)的替代函数分别求解关于稀疏反射系数x和噪声参数β的偏导数，并令偏导数为0，得到噪声参数β和稀疏反射系数x的第i个元素x_i的表达式，即：

和

式中，

n_j是字典矩阵A(θ)的第j行中非零元素的数量；

S10313、交替迭代更新噪声参数β和x，当内循环次数达到最大迭代次数T或者||x^{(t +1)}-x^(t)||²/||x^(t)||²＜ξ，ξ为内循环终止阈值时(ξ是很小的正实数，本发明实施例一中ξ可以取值为1e-5)，终止迭代，输出内循环最终迭代结果的稀疏反射系数作为此次外循环所要求的稀疏反射系数x^(k+1)，否则，返回S10312中式(5)和式(6)继续迭代。

S1032、固定所述包含参数化字典的目标函数中的x^(k)，更新墙体参数θ^(k+1)，即更新d^(k+1),ε^(k+1)，并用哈希表和线性卷积对更新公式中包含字典矩阵的相关运算进行替代。

在本发明实施例一中，S1032具体可以包括：

固定x^(k)时，目标函数(3)变为：

求使得目标函数J₂(θ)最小的θ(d,ε)，即：

采用一阶泰勒级数对参数化字典矩阵A(d,ε)展开处理，得：

式中，

由最小二乘得到墙体厚度估计的变化量Δd和墙体相对介电常数估计变化量Δε：

则得到墙体参数更新公式:

d^(k+1)＝d^(k)+Δd,ε^(k+1)＝ε^(k)+Δε (11)

当外循环达到最大迭代次数K或者||x^(k+1)-x^(k)||²/||x^(k)||²＜ζ，ζ为内循环终止阈值时(ξ是很小的正实数，本发明实施例一中ξ可以取值为1e-5)，终止迭代，否则，返回S1031开始继续迭代。

经分析，在交替迭代过程中，由式(5)、(6)和式(10)可以看出，在计算β、x(

和H_i)和墙体参数变化量Δd,Δε的过程中，都需要字典的构造、存储和计算，因此需要的存储空间大，计算复杂度也高。为解决此问题，采用哈希表和线性卷积对更新公式中包含字典矩阵的相关运算进行替代，使得对任意超宽带发射波形都可以在无需构建字典的情况下进行求解。

包含的字典矩阵的相关运算进行替代具体包括：

(1)β^(t+1)的计算

在式(5)中，Ax^(t)＝[(A₁x^(t))^T,(A₂x^(t))^T,...,(A_Lx^(t))^T]^T，只要求出A_lx^(t),l＝1,2,...L就可得到Ax^(t)。令

m_li(i＝1,2,...,N)，并假设m_li,i＝1,2,...,N的最小值和最大值分别为m_min和m_max，则有A_lx^(t)的第m个元素为：

记：

g_l[m]@A_lx^(t)[m] (13)

s(mT_s)@u[m] (14)

利用冲击函数的抽样特性，则式(12)变为：

令

则：

可以看出，式(12)的运算转换为式(16)的p_l和u的卷积计算，避免了A的构建、存储和计算。而且，p_l[n]通过将延时序列中满足m_li＝n,i＝1,2,...,N的第i个像素点所对应的

的值进行累加得到，其调用哈希表(即MATLAB的accumarray函数)就可以快速计算。

根据式(5)和式(16)，可得：

(2)

的计算

将式(2)和式(13)代入

可得：

令υ[M+1-m]@s(mT_s)，则：

可以看出，

的运算转换为υ和(y_l-g_l)的卷积计算，也避免了A的构建、存储和计算。

(3)H_i的计算

由式(2)代入H_i，可得：

令f[M+1-m]＝s²(mT_s)，则：

其中，n_l∈C^M×1，n_l[m]表示对应的第l个子矩阵A_l∈C^M×N的第m行元素的非零元素的个数。由以上分析可知，同样，H_i的运算转换为n_l和f的卷积计算，也避免了A的构建、存储和计算。

值得一提的是，如果要快速得到H_i，那么n_l的计算是至关重要的。因为发射信号是超宽带信号，也就是时域为短时信号，所以对于长度为M的发射序列u(m)(即s(mT_s))而言，存在一个整数Q，当满足m≤Q时，使得u(m)不为零。考虑到矩阵A_l的第i列为A_l＝[s(T_s-τ_li),s(2T_s-τ_li),...,s(MT_s-τ_li)]^T，则元素A_l(m,i)非零当m_li满足：

1≤m-m_li≤Q (22)

且离散时延m_li满足：

m_min≤m_li≤m_max (23)

所以有：

max(m_min,m-Q)≤m_li≤min(m_max,m-1) (24)

即：

式中，

表示第l个天线关于N个像素点的时延m_li,i＝1,2,...,N中满足m_li＝n,i＝1,2,...,N的m_li的个数，可通过哈希表快速求得。

(4)Δd,Δε的计算

在式(10)中，包含字典项有A(d^(k),ε^(k))x^(k)和B。

对于A(d^(k),ε^(k))x^(k)的计算，因其与上述g_l的计算方法相似，所以采用g_l的快速计算方法。前面已经讨论过，此处不再累述，以下是B的快速计算。

对于B的运算，令

则B＝[B₁,B₂]。下面给出B₁,B₂的快速计算方法。关于B₁的计算如下：

其中，

用第l个通道子字典矩阵A_l关于墙体厚度的偏导矩阵，

的第m行第i列的元素为：

记向量

的第m个元素为：

对第l个通道而言，利用冲击函数的抽样特性，得到：

其中

与式(16)中p_l[n]类似，q_l[n]同样可以用哈希表求得。

令：

则式(29)转化为：

将式(30)代入式(26)可知：

可以看出，B₁的运算转换为q_l和w的卷积计算，也避免了字典的构建、存储和计算。

关于B₂的计算如下：

其中，

用第l个通道子字典矩阵A_l关于墙体相对介电常数的偏导矩阵，

的第m行第i列的元素为：

记向量

的第m个元素为：

对第l个通道而言，利用冲击函数的抽样特性，得到：

其中

与(16)中p_l[n]类似，r_l[n]同样可以用哈希表求得。

将(35)代入(32)可知：

可以看出，B₂的运算转换为r_l和w的卷积计算，也避免了字典的构建、存储和计算。

S104、达到迭代终止条件时，输出外循环终止时对应的墙体厚度和相对介电常数作为墙体参数估计值，对应的稀疏反射系数用于成像。

本发明实施例一中，利用仿真模型对墙体厚度为0.2m、墙体相对介电常数为4.5的墙体后一个长0.5m，宽0.2m的扩展目标进行成像。预设墙体参数d＝0.15m,ε＝6.8，所用场景示意图如图2所示，使用L＝20个收发共置天线合成线性阵列，对SOI进行探测，阵列平行于墙体，放置在墙前1m处。仿真中采用发射波中心频率为1GHz的窄高斯脉冲信号。

穿墙雷达墙体参数估计结果和估计值的相对误差如表所示。

参数	墙体厚度d(m)	墙体相对介电常数ε
			真实值	0.2	4.5
估计值	0.2002	4.5136
			误差	0.0002	0.0136

由上表可以看出，墙体参数估计结果的误差很小，验证了本发明的正确性和有效性。所以本发明能够比较精确地估计出墙体参数，同时在无需构建字典矩阵的情况下实现了墙后目标的自聚焦成像，减少了计算复杂度和计算机的存储压力。

图3(a)是忽略墙体时采用延时求和算法直接进行成像的结果，图3(b)是采用延时求和算法，在真实仿真参数d＝0.2,ε＝4.5下进行成像的结果。可以发现，加上墙体补偿后，传统的成像方法也可以矫正目标位置，但所得到的成像结果中栅旁瓣严重，图像不够清晰。图4(a)为k＝1时的成像结果，此时的墙体参数为d＝0.1734,ε＝5.4896。图4(b)中最终自聚焦成像结果较好的消除了位置偏移，同时，与图3(b)相比，最终成像结果也更加清晰。

实施例二：

本发明实施例二提供了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现如本发明实施例一提供的穿墙雷达自聚焦稀疏成像方法的步骤。

实施例三：

图5示出了本发明实施例三提供的计算机设备的具体结构框图，一种计算机设备100包括：一个或多个处理器101、存储器102、以及一个或多个计算机程序，其中所述处理器101和所述存储器102通过总线连接，所述一个或多个计算机程序被存储在所述存储器102中，并且被配置成由所述一个或多个处理器101执行，所述处理器101执行所述计算机程序时实现如本发明实施例一提供的穿墙雷达自聚焦稀疏成像方法的步骤。

在本发明中，由于穿墙雷达自聚焦稀疏成像过程中，将扩展目标的稀疏特性作为先验信息，基于全变分约束的最大后验概率估计和所述参数化字典信号模型，构建包含参数化字典的目标函数，然后通过交替迭代更新稀疏反射系数和墙体参数。从而在进行墙体参数估计的同时，有效地保留了扩展目标的边缘特性，在消除墙体参数未知引起的目标位置偏移的同时使成像结果更清晰。在此基础上，通过用哈希表和线性卷积对更新公式中包含的字典矩阵的相关运算进行替代，从而在整个自聚焦过程中无需构建、存储和计算字典矩阵，将迭代循环中繁杂的计算过程进行简化，有效的减少运算时间，并在空间复杂度方面得到了改善。

本领域普通技术人员可以理解上述实施例的各种方法中的全部或部分步骤是可以通过程序来指令相关的硬件来完成，该程序可以存储于一计算机可读存储介质中，存储介质可以包括：只读存储器(ROM，Read Only Memory)、随机存取记忆体(RAM，RandomAccess Memory)、磁盘或光盘等。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (9)

1.一种穿墙雷达自聚焦稀疏成像方法，其特征在于，所述方法包括：

S101、获取扩展目标的回波信号，并构建以墙体厚度和相对介电常数作为参数的参数化字典信号模型；

S102、将扩展目标的稀疏特性作为先验信息，基于全变分约束的最大后验概率估计和所述参数化字典信号模型，构建包含参数化字典的目标函数；

S103、根据所述包含参数化字典的目标函数交替迭代更新稀疏反射系数和墙体参数，其中，用哈希表和线性卷积对更新公式中包含的字典矩阵的相关运算进行替代；

S104、达到迭代终止条件时，输出外循环终止时对应的墙体厚度和相对介电常数作为墙体参数估计值，对应的稀疏反射系数用于成像。

2.如权利要求1所述的方法，其特征在于，所述S101具体包括：

将成像区域离散化为N个像素网格点，堆叠对成像区域进行探测的L个收发共置天线的数据，并构建参数化字典信号模型，可表示为：

式中，x是稀疏反射系数，x＝[x₁,x₂,...,x_N]^T，x∈C^N×1，C^N×1表示N×1维的向量，N为总像素网格点数；y为回波数据矢量，y∈C^ML×1，M为每个收发共置天线处接收回波的采样点数；e为噪声矢量，e∈C^ML×1；

为参数化字典矩阵，

d是墙体厚度，ε是墙体相对介电常数，

且

是用第l个收发共置天线构建的子矩阵；

的第m行第i列的元素为：

A_l(m,i)＝s(mT_s-τ_li) (2)

其中，s(mT_s)是发射信号序列，m＝1,2,...,M；若记第l个天线位置坐标为(x_p,y_p)，第i个像素点位置坐标为(x_a,y_a)，电磁波在自由空间的传播速度为c，第l个收发共置天线和第i个像素点之间的传输时延为τ_li，则有

其中，

3.如权利要求2所述的方法，其特征在于，所述S102具体包括：

假设噪声矢量e服从均值为0，协方差矩阵为βI的复高斯分布，I为单位矩阵，通过将墙体参数引入到基于全变分约束的最大后验概率估计的稀疏成像方法当中，构建包含参数化字典矩阵的目标函数，可表述为：

其中，D为一阶差分算子，

β是噪声参数。

4.如权利要求3所述的方法，其特征在于，所述S103具体包括：

S1031、固定所述包含参数化字典的目标函数中的墙体参数，即d^(k),ε^(k)，

更新稀疏反射系数x^(k+1)，并用哈希表和线性卷积对更新公式中包含字典矩阵的相关运算进行替代；其中k为当前外循环迭代次数，令外循环最大循环次数为K，初始状态下k＝0，墙体厚度为d⁽⁰⁾，墙体相对介电常数为ε⁽⁰⁾；

S1032、固定所述包含参数化字典的目标函数中的x^(k)，更新墙体参数

即更新d^{(k +1)},ε^(k+1)，并用哈希表和线性卷积对更新公式中包含字典矩阵的相关运算进行替代。

5.如权利要求4所述的方法，其特征在于，所述S1031具体包括：

固定

时，

为确定性矩阵，记

则包含参数化字典矩阵的目标函数(3)为：

此时，通过交替迭代更新噪声参数β和x，求解使得目标函数J₁最小的最优x；

具体过程如下：

S10311、令内循环最大循环次数为T，当前迭代次数为t，初始状态下t＝0，采用延时求和得到稀疏反射系数的初值x⁽⁰⁾；

S10312、对式(4)中的目标函数J₁(x,β)进行优化最小化替代，然后对目标函数J₁(x,β)的替代函数分别求解关于稀疏反射系数x和噪声参数β的偏导数，并令偏导数为0，得到噪声参数β和稀疏反射系数x的第i个元素x_i的表达式，即：

和

式中，

n_j是字典矩阵

的第j行中非零元素的数量；

S10313、交替迭代更新噪声参数β和x，当内循环次数达到最大迭代次数T或者||x^(t+1)-x^(t)||²/||x^(t)||²＜ξ时，终止迭代，输出内循环最终迭代结果的稀疏反射系数作为此次外循环所要求的稀疏反射系数x^(k+1)，否则，返回S10312中式(5)和式(6)继续迭代，其中ξ为内循环终止阈值。

6.如权利要求5所述的方法，其特征在于，所述S1032具体包括：

固定x^(k)时，目标函数(3)变为：

求使得目标函数

最小的

即：

采用一阶泰勒级数对参数化字典矩阵A(d,ε)展开处理，得：

式中，

由最小二乘得到墙体厚度估计的变化量Δd和墙体相对介电常数估计变化量Δε：

则得到墙体参数更新公式:

d^(k+1)＝d^(k)+Δd,ε^(k+1)＝ε^(k)+Δε (11)

当外循环达到最大迭代次数K或者||x^(k+1)-x^(k)||²/||x^(k)||²＜ζ时，终止迭代，否则，返回S1031开始继续迭代。

7.如权利要求6所述的方法，其特征在于，包含的字典矩阵的相关运算进行替代具体包括：

关于β^(t+1)的计算如下：

在式(5)中，Ax^(t)＝[(A₁x^(t))^T,(A₂x^(t))^T,...,(A_Lx^(t))^T]^T，求出A_lx^(t) _,l＝1,2,...L就可得到Ax^(t)；令

并假设m_li,i＝1,2,...,N的最小值和最大值分别为m_min和m_max，则有A_lx^(t)的第m个元素为：

记：

g_l[m]@A_lx^(t)[m] (13)

s(mT_s)@u[m] (14)

利用冲击函数的抽样特性，则式(12)变为：

令

则：

根据式(5)和式(16)，得到：

关于

的计算如下：

将式(2)和式(13)代入

得：

令υ[M+1-m]@s(mT_s)，则：

关于H_i的计算如下：

由式(2)代入H_i，得：

令f[M+1-m]＝s²(mT_s)，则：

其中，n_l∈C^M×1，n_l[m]表示对应的第l个子矩阵A_l∈C^M×N的第m行元素的非零元素的个数；

对于长度为M的发射序列u(m)，即s(mT_s)而言，存在一个整数Q，当满足m≤Q时，使得u(m)不为零；考虑到矩阵A_l的第i列为A_l＝[s(T_s-τ_li),s(2T_s-τ_li),...,s(MT_s-τ_li)]^T，则元素A_l(m,i)非零当m_li满足：

1≤m-m_li≤Q (22)

且离散时延m_li满足：

m_min≤m_li≤m_max (23)

所以有：

max(m_min,m-Q)≤m_li≤min(m_max,m-1) (24)

即：

式中，

表示第l个天线关于N个像素点的时延m_li,i＝1,2,...,N中满足m_li＝n,i＝1,2,...,N的m_li的个数，通过哈希表求得。

关于Δd,Δε的计算如下：

在式(10)中，包含字典项有A(d^(k),ε^(k))x^(k)和B；对于A(d^(k),ε^(k))x^(k)的计算，采用g_l的计算方法；

以下是B的计算：

对于B的运算，令

则B＝[B₁,B₂]，关于B₁的计算如下：

其中，

用第l个通道子字典矩阵A_l关于墙体厚度的偏导矩阵，

的第m行第i列的元素为：

记向量

的第m个元素为：

对第l个通道而言，利用冲击函数的抽样特性，得到：

其中

与式(16)中p_l[n]类似，q_l[n]用哈希表求得；

令：

则式(29)转化为：

将式(30)代入式(26)可知：

关于B₂的计算如下：

其中，

用第l个通道子字典矩阵A_l关于墙体相对介电常数的偏导矩阵，

的第m行第i列的元素为：

记向量

的第m个元素为：

对第l个通道而言，利用冲击函数的抽样特性，得到：

其中

与(16)中p_l[n]类似，r_l[n]用哈希表求得；

将(35)代入(32)可知：

8.一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现如权利要求1至7任一项所述的穿墙雷达自聚焦稀疏成像方法的步骤。

9.一种计算机设备，包括：

一个或多个处理器；

存储器；以及

一个或多个计算机程序，所述处理器和所述存储器通过总线连接，其中所述一个或多个计算机程序被存储在所述存储器中，并且被配置成由所述一个或多个处理器执行，其特征在于，所述处理器执行所述计算机程序时实现如权利要求1至7任一项所述的穿墙雷达自聚焦稀疏成像方法的步骤。

Images