CN103149561A - 一种基于场景块稀疏的稀疏微波成像方法 - Google Patents

Abstract

本发明公开了一种基于场景块稀疏的稀疏微波成像方法，涉及雷达成像技术，利用了实际成像场景中目标的几何尺寸远大于高分辨成微波像雷达的空间分辨率，并且通常覆盖多个雷达分辨单元这一事实，采用压缩采样的方式对块稀疏场景进行稀疏观测，基于场景的块稀疏特性进行高效的稀疏图像重构。本发明方法与基于匹配滤波成像的方法不同，稀疏微波成像雷达以稀疏信号处理为基础，通过求解一个优化问题恢复被观测的稀疏或者变换域稀疏场景。本发明方法在同样采样条件下，实现对稀疏成像场景比常规稀疏微波成像更为精确的重构，也可以在保证同样重构性能条件下，实现对稀疏场景更为稀疏的观测，从而进一步降低采样的数据率。

Description

一种基于场景块稀疏的稀疏微波成像方法

技术领域

本发明涉及雷达成像技术领域，具体涉及一种基于场景块稀疏的稀疏微波成像方法，适用于高分辨稀疏微波成像，为稀疏微波成像雷达图像的高精度重构提供技术保障。

背景技术

合成孔径雷达(Synthetic Aperture Radar，SAR)作为一种可全天时、全天候工作的主动式高分辨率微波遥感设备，它已经被广泛应用于军事和民用对地观测等相关领域中，如地球资源环境勘察与高精度测绘、自然灾害监测与突发事件应对、宽幅海洋资源和海面舰船监视、侦察与预警等。随着半导体器件的迅猛发展以及微波成像理论和方法的不断发展，实现米级甚至厘米级的分辨率以及上千公里的测绘带成像能力的微波成像雷达在技术上已经成为可能。

常规微波成像技术基于匹配滤波技术实现雷达场景的二位反卷积，它以奈奎斯特采样定理为基础保证雷达场景的无模糊重建。奈奎斯特采样定理要求无失真恢复有限带宽信号所需的采样率不得低于信号带宽，因而奈奎斯特定理要求的采样率决定了微波成像系统的规模和复杂度。由于传统微波成像技术需要发射/接收宽带距离信号以获得较高的距离分辨率，同时获得高方位分辨率需要较大的多普带宽信号，对距离信号和方位信号进行耐奎斯特采样不可避免的增加了系统的数据采集率、数据存储量以及对下行数据链路传输能力的要求，因而增加了系统的复杂度。

几年来随着诸如宽幅海洋监视的宽测绘度、高分辨率微波成像需求的不断增长，人们寻求一种可以对宽幅场景中的稀疏目标高分辨成像的稀疏微波成像新体制。稀疏微波成像是指将稀疏信号处理理论引入微波成像并有机结合形成的微波成像新理论、新体制和新方法，它以压缩采样与恢复理论为基础，以稀疏信号处理为技术手段，通过寻找被观测对象的稀疏表征域(如空间、时间、极化域)，在此稀疏表征域内以二维降速率采样为观测条件，采用l_q(0＜q≤1)稀疏恢复算法获取被观测对象的空间位置、散射特征和运动特性等几何与物理特征。稀疏微波成像雷达利用了微波场景的稀疏性或者变换稀疏性这一特点，以远低于耐奎斯特速率要求的采样率对稀疏微波场景进行采样，大大减少了数据采集、存储以及传输量，降低了系统复杂度。

稀疏微波成像适用于微波场景具有稀疏特征的雷达成像领域，如海面舰船监视、空中目标ISAR成像等，具有较强后向散射系数目标区域的几何、物理特征的重构性能为稀疏微波成像系统的首要侧重点。由于稀疏微波成像技术以稀疏场景为先验信息，采用较少的采样实现强散射目标的检测、位置估计以及后向散射系数重构。但是在实际应用中，由于杂波的存在使得任意的成像场景不可能是理想稀疏的，这将明显抑制了稀疏微波成像的性能，在实际操作中通常设定一个阈值，当估计值小于某个阈值时即认为该值为零进而使其强制性的满足稀疏性条件。微波图像重构的性能严重依赖于该阈值的选取，阈值选取较大导致低检测概率，阈值选取过小将导致高虚警概率；另一方面，由于成像雷达的高分辨，稀疏微波成像场景中的感兴趣的目标往往覆盖了多个分辨单元，即如果某个分辨单元为强散射点，那么它周围的若干个分辨单元往往也是强散射点，即除非人为设定的角反射器外不存在孤立的强散射点目标。现实成像场景的这一先验信息可以用来明显改善稀疏微波成像的性能，它也可以在保证同等图像重构的前提下进一步降低对微波场景的观测数据量。本专利提出基于稀疏微波成像雷达的稀疏场景的块稀疏特性的高精度微波图像重构方法。以下列出本发明的参考文献，通过引用将它们并入于此，如同在本说明书中作了详尽描述。

参考文献：

[1]D.L.Donoho，“Compressed sensing，”IEEE Trans.Inform.Theory，

vol.52，no.4，pp.1289-1306，Apr.2006.

[2]E.J.Cand`es，M.B.Wakin，and S.P.Boyd，“Enhancing sparsityby reweighted `1 minimization.”The Journal of Fourier Analysis andApplications (Preprint).

[3]Richard G.Baraniuk and Volkan Cevher，Volkan Cevher andChinmay Hegde，“Model-Based Compressive Sensing，”IEEE Trans.Inform.Theory，vol.56，no.4，pp.1982-2001，Apr.2006.

[4]Mihailo Stojnic，Farzad Parvaresh and Babak Hassibi，“On theReconstruction of Block-Sparse Signals With an Optimal Number ofMeasurements，”IEEE Trans.Signal Process.，vol.57，no.8，pp.3075-3085，Aug.2009.

[5]Vishal M.Patel，Glenn R.Easley，Dennis M.Healy，Healy，Jr.，and Rama Chellappa，“Compressed Synthetic Aperture Radar，”IEEE J.Sel.Topics Signal Process.，vol.4，no.2 pp.244-254，Apr.2010.

[6]M.Herman and T.Strohmer，“High resolution radar viacompressed sensing，”IEEE Trans.Signal Process.，vol.57，no.6，pp.2275-2284，Jun.2009.

[7]Lin Yueguan，Zhang Bingchen，Hong Wen and Wu Yirong.Multi-channel SAR imaging based on Distributed Compressive Sensing.Science in China Series：F：Information Sciences.

发明内容

本发明的目的是公开一种基于场景块稀疏的稀疏微波成像方法，利用稀疏目标的分块稀疏特性，及块稀疏先验信息进行图像重构，以克服现有稀疏重构算法的高虚警概率或低检测概率的缺点。

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

一种基于块稀疏的稀疏微波成像方法，其包括以下步骤：

(1)利用成像场景的稀疏性，对场景进行稀疏观测，得到降采样率的原始回波数据；

(2)利用稀疏优化的方法，基于稀疏观测到的原始回波数据重构微波成像场景；

(3)基于场景的块稀疏性，实现快稀疏场景的稀疏重构：利用实际成像场景的分块稀疏特性，将稀疏微波成像场景在距离向和方位向分割为许多小的矩形块，分别计算这些矩形块内的所有分辨单元的后向散射系数的某种平均结果：算术平均、几何平均、以及由它们构成向量的2-范数

并以此计算迭代场景重建算法中阈值。

所述的稀疏微波成像方法，其所述步骤(1)中的稀疏观测，包括：

假设稀疏的成像场景按照分辨单元大小被离散化表示为x_l，k，其中l，k分别代表距离和方位下标，雷达平台向成像场景发射宽带脉冲信号s(t)；给定雷达脉冲信号的发射时刻(慢时间)、及对每个脉冲的采样时刻(慢时间)，含有噪声n的原始离散回波数据y和成像场景x满足关系式：y＝Φ·x+n，其中Φ定义为雷达观测矩阵；

s(t)是线性调频信号，或是宽带波形编码信号，且雷达多次发射的宽带脉冲信号相同或不同；

以最小化观测矩阵Φ的所有不同两列的互相关系数的最大值为准则确定发射时刻、采样时刻、及选择波形s(t)。

所述的稀疏微波成像方法，其所述步骤(2)中的场景重构方法，包括：

3a、给定雷达观测矩阵Φ、观测噪声的功率α的任意初始估计值

对角元素均为非负数的对角矩阵Q的任意初始估计值

观测数据y以及一个充分小的正数ε(如ε＝10^-6)；

3b、令n＝0，1，2，3，…，执行如下迭代步骤：

3b1、计算： ${\hat{x}}^{(n+1)}={({\mathrm{\Φ}}^H\mathrm{\Φ}+{\widehat{\mathrm{\α}}}^{\left(n\right)}{\hat{Q}}^{\left(n\right)})}^{-1}{\mathrm{\Φ}}^{H}y;$

3b2、更新： ${\widehat{\mathrm{\α}}}^{(n+1)}={\left|\right|y-\mathrm{\Φ}{\hat{x}}^{(n+1)}\left|\right|}_2^2,$ 定义 ${\hat{q}}_i^{(n+1)}={\left|{\hat{x}}_i^{(n+1)}\right|}^2+\mathrm{\ϵ},$

更新 ${\hat{Q}}^{(n+1)}=\mathrm{diag}(1/{\hat{q}}_1^{(n+1)},\·\·\·,1/{\hat{q}}_N^{(n+1)});$

3c、判断：如果连续两次迭代估计值

小于预先给定的常数或者迭代次数达到人为设定的上限，迭代结束；否则n值自动加1，重复上述计算、更新两个步骤。

所述的稀疏微波成像方法，其所述步骤3b中场景重构的迭代方法，在每次迭代计算的过程中利用场景的一维或两维分块稀疏特性修改迭代步骤。

所述的稀疏微波成像方法，其所述使用场景的一维分块稀疏特性使用如下迭代步骤，包括以下步骤：

(1)给个观测矩阵Φ、α的初始估计值

Q的初始估计值

一个假定的正常数ρ、观测数据y并令n＝0；

(2)计算 ${\hat{x}}^{(n+1)}={({\mathrm{\Φ}}^H\mathrm{\Φ}+{\widehat{\mathrm{\α}}}^{\left(n\right)}{\hat{Q}}^{\left(n\right)})}^{-1}{\mathrm{\Φ}}^{H}y;$

(3)更新 ${\widehat{\mathrm{\α}}}^{(n+1)}={\left|\right|y-\mathrm{\Φ}{\hat{x}}^{(n+1)}\left|\right|}_2^2,$

${\mathrm{\δ}}_d^{(n+1)}=\left(\frac{1}{D}\underset{l=(d-1)N/D+1}{\overset{\mathrm{dN}/D}{\mathrm{\Σ}}}{\left|{\hat{x}}_l^{(n+1)}\right|}^2\right),$

${\hat{q}}_{(d-1)N/D+1}^{(n+1)}=\&CenterDot;\&CenterDot;\&CenterDot;={\hat{q}}_{\mathrm{dN}/D}^{(n+1)}=\left\{\begin{array}{cc}{\mathrm{\&delta;}}_d^{(n+1)}-\mathrm{\&rho;}+\mathrm{\&epsiv;},& \mathrm{if}{\mathrm{\&delta;}}_d^{(n+1)}\&GreaterEqual;\mathrm{\&rho;}\\ \mathrm{\&epsiv;},& \mathrm{if}{\mathrm{\&delta;}}_d^{(n+1)}<\mathrm{\&rho;}\end{array}\right.,$

${\hat{Q}}^{(n+1)}=\mathrm{diag}(1/{\hat{q}}_1^{(n+1)},\&CenterDot;\&CenterDot;\&CenterDot;,1/{\hat{q}}_N^{(n+1)});$

(4)如果连续两次迭代估计值

小于预先给定的常数或者迭代次数达到人为设定的上限，迭代结束；否则n值自动加1，重复步骤(2)和(3)。

所述的稀疏微波成像方法，其采用基于场景块稀疏特性的阈值迭代类算法、或贪心类的重构方法，包括OMP算法、CoSaMP算法、TWIST算法、FPC算法、或Lq正则化方法其中之一，及它们的组合。

本发明的一种基于块稀疏的稀疏微波成像方法，是指将稀疏信号处理理论引入微波成像并有机结合形成的微波成像的新理论、新体制和新方法。与基于匹配滤波成像的方法不同，稀疏微波成像雷达以稀疏信号处理为基础，通过求解一个优化问题恢复被观测的稀疏或者变换域稀疏场景。本发明方法在同样采样条件下，实现对稀疏成像场景比常规稀疏微波成像更为精确的重构，也可以在保证同样重构性能条件下，实现对稀疏场景更为稀疏的观测，从而进一步降低采样的数据率。

附图说明

图1本发明的一种基于块稀疏的稀疏微波成像方法流程图；

图2具有分块稀疏目标的微波成像场景示意图；

图3稀疏微波成像雷达成像流程图；

图4基于块稀疏的信号平均方法示意图；其中：

图4a为一维稀疏结构示意图；

图4b为二维稀疏结构示意图；

图5软门限限幅加权示意图；

图6基于块稀疏特性的稀疏微波重构图像示意图。

具体实施方式

本发明的一种基于块稀疏的稀疏微波成像方法，根据雷达系统参数和平台几何关系建立的稀疏微波成像模型，包括：稀疏观测矩阵、成像雷达观测矩阵、稀疏化变换矩阵、稀疏系数矢量矩阵。

y＝Θ·H·x+n□Φ·x+n

其中，

表示稀疏微波成像雷达对稀疏场景进行观测得到的离散回波采样数据；

表示典型稀疏的离散成像场景的后向散射系数矢量(x中只有少数元素明显非零，其他元素近似为零。)；

是稀疏微波成像雷达观测矩阵，

是稀疏采样矩阵，

是成像雷观测矩阵，

表示系统热噪声。

对于上式中的由稀疏场景的各个分辨单元后向散射系数构成的稀疏变量

注意到x的元素的幅度之间具有一定的相关性，即如果x的某个元素为非零的强散射点，那么这一点周围的元素往往也是非零的强点，而且明显非零元素的出项大多是成块出现的，使得场景表现出分块稀疏的特点。

如成像场景是二维的，如距离+方位，此时向量x由二维矩阵

代替(N和L是二维场景的大小)，类似一维块稀疏特性，二维场景在两个维度均表现为块稀疏结构，即二维场景中的一个分辨单元

为强散射点，那么与它相邻分辨单元X_i，j，(|i-i₀|+|j-j₀|＜d，d为正数)通常也是明显非零的的强点。对二维回波采样数据采用按列展开的形式，仍将其表示为一维信号形式，如上的稀疏微波成像输入输出模型依然成立。

利用成像场景x的稀疏性这一特点，稀疏微波成像雷达采用L-1优化理论重构成像场景x，

$\hat{x}=\arg \min {\left|x\right|}_1$

s.t.‖y-Φ·x‖₂＜ε_n

这里的参数ε_n控制对观测噪声的容忍程度，L-1范数定义为|x|₁＝∑_k|x_k|。

如上的优化问题为一个二次非光滑问题，需要将其转化为内点法求解。在稀疏微波成像系统中，成像场景中除了具有较大幅度的目标外，其余的杂波点的幅度并非严格为零，这些小的杂波点将导致如上凸优化求解得到的整个图像重构性能的明显恶化。另一方面，为了克服这些非理想稀疏特性带来的影响，同时降低如上优化问题的计算复杂度，基于阈值迭代的L-1正则化方法受到了广泛关注，它的核心思想是采用迭代方法得到如下的正则化问题的解：

$\hat{x}=\arg \min \{{\left|\right|y-\mathrm{\&Phi;}\&CenterDot;x\left|\right|}_2^2+\mathrm{\&lambda;}{\left|x\right|}_1\}$

迭代公式可以统一写成：

${\hat{x}}^{(n+1)}=H_s({\hat{x}}^{\left(n\right)}+{\mathrm{\&mu;\&Phi;}}^H(y-\mathrm{\&Phi;}\&CenterDot;x)).$

以上两式中的λ、s分别为正则化参数以及阈值，阈值迭代函数H_s(□)是一个非线性算子，它设置所有幅度小于阈值s的元素为一个充分小的正常数。经过若干次迭代后，算法满足收敛条件即可以得到稀疏信号x的重构结果。由于在稀疏微波成像雷达体制下，对稀疏场景的采样通常采用降速率采样方式，即观测量y的维数小于未知变量x的维数，即观测方程为欠定方程。另一方面由于未知场景变量x在实际中并不是严格稀疏的，常规的稀疏微波成像方法具有较高虚警概率或者较低的检测概率。

为了克服以上缺点，利用微波场景的后向散射系数的高斯分布特性，即x_i□N(0，σ_i)这一先验信息(σ_i为场景中第i个分辨单元的后向散射系数的标准差)，得到未知场景变量x的贝叶斯估计：

$x=\max P\left(x\right|y)$

$=\max P\left(y\right|x)\&CenterDot;P\left(x\right)$

$\&Proportional;\max \{\log P\left(y\right|x)+\log P\left(x\right)\}$

$\&Proportional;\max \{{\left|\right|y-\mathrm{\&Phi;x}\left|\right|}_2^2+\mathrm{\&alpha;}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\left|\frac{x_i}{{\mathrm{\&sigma;}}_i}\right|}^2\}$

$=\max \{{\left|\right|y-\mathrm{\&Phi;x}\left|\right|}_2^2+\mathrm{\&alpha;}{x}^H\mathrm{Qx}\}$

式中α为已知常数，

此时给定σ₁，…，σ_N，可以得到x的最大后验估计为

$\hat{x}={({\mathrm{\&Phi;}}^H\mathrm{\&Phi;}+\mathrm{\&alpha;Q})}^{-1}{\mathrm{\&Phi;}}^{H}y$

这里σ_i为x₁，…，x_N的标准差，其理论估计值为

α的估计值为

给定任意的初值：

结合这两点可以得到x的迭代计算公式，

$\left\{\begin{array}{c}{\hat{x}}^{(n+1)}={({\mathrm{\&Phi;}}^H\mathrm{\&Phi;}+{\widehat{\mathrm{\&alpha;}}}^{\left(n\right)}{Q}^{\left(n\right)})}^{-1}{\mathrm{\&Phi;}}^{H}y\\ {\widehat{\mathrm{\&alpha;}}}^{(n+1)}={\left|\right|y-\mathrm{\&Phi;}{\hat{x}}^{(n+1)}\left|\right|}_2^2\\ {\hat{q}}_i^{(n+1)}={\left|{\hat{x}}_i^{(n+1)}\right|}^2\\ {\hat{Q}}^{(n+1)}=\mathrm{diag}(1/{\hat{q}}_1^{(n+1)},\&CenterDot;\&CenterDot;\&CenterDot;,1/{\hat{q}}_N^{(n+1)})\end{array}\right.$

考虑到x的块稀疏特点，将上一次估计得到的分成D等分，计算：

${\mathrm{\&delta;}}_d^{(n+1)}=\frac{1}{D}\underset{l=(d-1)N/D+1}{\overset{\mathrm{dN}/D}{\mathrm{\&Sigma;}}}{\left|{\hat{x}}_l^{(n+1)}\right|}^2$

另一方面，根据

的计算公式，如果某个

值趋于零，那么的幅度也将趋近于零。为了明显消除杂波的影响，设置一个阈值ρ使得

小于ρ的区间内的所有都被赋予一个很小的正常数ε。因而以上的迭代稀疏重构方法修改为：

$\left\{\begin{array}{c}{\hat{x}}^{(n+1)}={({\mathrm{\&Phi;}}^H\mathrm{\&Phi;}+{\widehat{\mathrm{\&alpha;}}}^{(n+1)}{Q}^{\left(n\right)})}^{-1}{\mathrm{\&Phi;}}^{H}y\\ {\widehat{\mathrm{\&alpha;}}}^{(n+1)}={\left|\right|y-\mathrm{\&Phi;}{\hat{x}}^{(n+1)}\left|\right|}_2^2\\ {\mathrm{\&delta;}}_d^{(n+1)}={\frac{1}{D}\underset{/=(d-1)N/D+1}{\overset{\mathrm{dN}/D}{\mathrm{\&Sigma;}}}\left|{\hat{x}}_l^{(n+1)}\right|}^2\\ {\hat{q}}_{(d-1)N/D+1}^{(n+1)}=\&CenterDot;\&CenterDot;\&CenterDot;={\hat{q}}_{\mathrm{dN}/D}^{(n+1)}=\left\{\begin{array}{cc}{\mathrm{\&delta;}}_d^{(n+1)}-\mathrm{\&rho;}+\mathrm{\&epsiv;},& \mathrm{if}{\mathrm{\&delta;}}_d^{(n+1)}\&GreaterEqual;\mathrm{\&rho;}\\ \mathrm{\&epsiv;},& \mathrm{if}{\mathrm{\&delta;}}_d^{(n+1)}<\mathrm{\&rho;}\end{array}\right.\\ {\hat{Q}}^{(n+1)}=\mathrm{diag}(1/{\hat{q}}_1^{(n+1)},\&CenterDot;\&CenterDot;\&CenterDot;,1/{\hat{q}}_N^{(n+1)})\end{array}\right.$

ε为一个充分小的常数。如上的迭代图像重构方法考虑了x的一维块稀疏性，也可以将二维块稀疏性质应用到稀疏重构中以进一步改善重构性能，降低目标的虚警概率和漏警概率。

下面参照附图对本发明一种基于场景块稀疏的稀疏微波成像方法的具体实施方式进行详细的说明。为了清楚和简明起见，在说明书中并未描述实际实施方式的所有特征。

本发明方法主要集中于稀疏微波场景的成像，图2给出了一幅稀疏微波场景的后向散射系数(RCS)的示意图，在一片均匀的海面上存在6条船只，这些船只和海面背景相比具有较大的后向散射系数，它们的幅度明显非零，海面对雷达照射信号的反射能量较小，它们的后向散射实际虽然非零但近似为零。稀疏微波成像雷达的工作原理是：根据特定的参数设置，随着雷达平台的运动断续的向成像场景发射宽带脉冲信号，对待成像场景进行稀疏微波观测，得到离散时间采样原始回波数据，利用这些回波数据重构场景的RCS。假设在第η时刻雷达发射脉冲为s(t)，接收来自一个二维成像区域Ω的回波信号经过混频、采样(方位、距离采样时刻假设为

)，二维回波数据可以表述成：

$s_R(\mathrm{\&eta;},t)=\underset{\mathrm{\&Omega;}}{\&Integral;}{\mathrm{\&sigma;}}_{\mathrm{\&omega;}}{w}_a^{i,j}\left(\mathrm{\&eta;}\right)\&CenterDot;\mathrm{rect}\left(\frac{t-2R_{i,j}\left(\mathrm{\&eta;}\right)/c}{T_p}\right)\&CenterDot;\exp \{-j\frac{4\mathrm{\&pi;}{f}_0{R}_{i,j}\left(\mathrm{\&eta;}\right)}{c}\}\&CenterDot;s(t-\frac{2R_{i,j}\left(\mathrm{\&eta;}\right)}{c})\mathrm{d\&omega;}$

$\&ap;\underset{i,j}{\mathrm{\&Sigma;}}{\mathrm{\&sigma;}}_{i,j}{w}_a^{i,j}\left(\mathrm{\&eta;}\right)\&CenterDot;\mathrm{rect}\left(\frac{t-2R_{i,j}\left(\mathrm{\&eta;}\right)/c}{T_p}\right)\&CenterDot;\exp \{-j\frac{4\mathrm{\&pi;}{f}_0{R}_{i,j}\left(\mathrm{\&eta;}\right)}{c}\}\&CenterDot;s(t-\frac{2R_{i,j}\left(\mathrm{\&eta;}\right)}{c})$

其中

为雷达天线在η时刻照射到第(i，j)个分辨单元的方位向双程天线方向图，R_i，j(η)为η时刻雷达天线到第(i，j)个分辨单元的距离，c为光速，f₀为载频，T_p为脉冲信号的时间长度，rect(□)为矩形窗函数。

将二维回波数据以及二维场景的RCS按列展开，可以得到修改的观测信号形式：

y＝Φx+n

其中y和x为别表示二维回波数据以及二维场景的RCS按列展开，n为服从复高斯分布n□N(0，σI)的观测噪声。压缩感知理论表明对一个二维稀疏场景进行成像并不需要对其进行Nyquist采样，因而允许在如上的测量方程中，未知数的个数可能远大于测量方程个数。稀疏微波成像的主要目的就是利用场景的稀疏性，使用很少的观测值重构稀疏场景的RCS。

基于压缩感知的稀疏微波成像方法采用最小化x的p-范数(0＜p＜1)的优化问题重构稀疏微波场景，在未知场景稀疏度(明显非零元素的个数)的情况下可以通过以下优化问题的求解来实现：

$\hat{x}=\arg \min \{{\left|\right|y-\mathrm{\&Phi;}\&CenterDot;x\left|\right|}_2^2+\mathrm{\&lambda;}{\left|x\right|}_p^p\}$

在稀疏度已知时，可以采用正交匹配追踪(OMP)，压缩采样匹配追踪(CoSAMP)等算法实现。

以上两类稀疏重构方法需要事先知道信号的稀疏度或者正则化参数λ，如果设置较大的稀疏度或较小正则化参数将会导致较高的虚警概率，反之会导致较高的漏警概率。在稀疏微波成像系统设计时，设计者并不能精确得到这些先验信息，而且如上的重构算法中没有利用除了场景稀疏性之外的任何信息，因此常规的稀疏微波成像雷达很难得到令人满意的性能。

对于回波信号中来自成像场景中某个特定分辨单元的信号形成分，它是由包含了大量离散散射体回波的叠加，当该分辨单元与微波相互作用时，每个散射体都产生具有特定相位和幅度的后向散射波，因而来自该分辨单元目标的总回波为这些后向散射波的叠加，最终将该分辨单元目标的后向散射系数等效为一个复高斯随机变量。即成像场景x的每个元素均服从具有不同方差的零均值复高斯分布(x_i□N(0，σ_i)，σ_i为x_i的标准差)。结合这一先验信息和雷达测量值，稀疏场景稀疏的后向散射稀疏向量x的最大后验估计为：

$\hat{x}=\max \{{\left|\right|y-\mathrm{\&Phi;x}\left|\right|}_2^2+\mathrm{\&alpha;}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{{\left|\right|x_i\left|\right|}^2}{{\mathrm{\&sigma;}}_i^2}\}$

α为由观测噪声功率决定。定义对角矩阵 $Q=\mathrm{diag}(1/{\mathrm{\&sigma;}}_1^2,\&CenterDot;\&CenterDot;\&CenterDot;,1/{\mathrm{\&sigma;}}_N^2),$ 如果α和σ_i均为已知，则：

$\hat{x}={({\mathrm{\&Phi;}}^H\mathrm{\&Phi;}+\mathrm{\&alpha;Q})}^{-1}{\mathrm{\&Phi;}}^{H}y$

在给定x的估计值后，需要更新α、Q，α为观测噪声功率σ_n的估计值，Q的第i个对角线元素为x_i的方差的估计值

的倒数。

$\mathrm{\&alpha;}=\widehat{\mathrm{\&sigma;}}={\left|\right|y-\mathrm{\&Phi;}\hat{x}\left|\right|}_2^2$

${\widehat{\mathrm{\&sigma;}}}_i^2=\left|\right|{\hat{x}}_i^2\left|\right|$

基于以上分析的重构算法的迭代形式如下(参见1和图3)：

输入：观测矩阵Φ、x的初始估计值

α的初始估计值

的初始估计值

充分小的正常数ε，观测数据y，

令n＝0，1，2，3，…，执行如下迭代步骤

1.计算： ${\hat{x}}^{(n+1)}={({\mathrm{\&Phi;}}^H\mathrm{\&Phi;}+{\widehat{\mathrm{\&alpha;}}}^{\left(n\right)}{\hat{Q}}^{\left(n\right)})}^{-1}{\mathrm{\&Phi;}}^{H}y;$

2.更新： ${\widehat{\mathrm{\&alpha;}}}^{(n+1)}={\left|\right|y-\mathrm{\&Phi;}{\hat{x}}^{(n+1)}\left|\right|}_2^2,$ ${\widehat{\mathrm{\&sigma;}}}_i^{(n+1)}={\left|{\hat{x}}_i^{(n+1)}\right|}^2+\mathrm{\&epsiv;},$

${\hat{Q}}^{(n+1)}=\mathrm{diag}(1/{\widehat{\mathrm{\&sigma;}}}_1^{2(n+1)},\&CenterDot;\&CenterDot;\&CenterDot;,1/{\widehat{\mathrm{\&sigma;}}}_N^{2(n+1)});$

判断：如果连续两次迭代估计值

小于预先给定的常数或者迭代次数达到人为设定的上限，迭代结束，输出x的估计值

否则n值自动加1，重复上述计算、更新两个步骤。以上方法适用于理想的稀疏场景重构。但是由于微波场景的实际非理想稀疏特性，如上的成像方法在实际应用中的效果并不理想。重构考虑到稀疏成像场景的x的块稀疏性，即在高分辨雷达成像中，实际目标的大小通常覆盖了多个分辨单元，因而使得整个成像场景表现为分块稀疏的，即如果某个分辨单元的散射强度较大，那么与它相邻的若干个距离以及方位分辨单元的也是大概率的强散射点。基于这一先验信息，假设信号长度N可以被整数D整除，将上一次估计得到的分成D等分，参加图4(a)。计算每个等分区间内的信号的平均功率：

${\mathrm{\&delta;}}_d^{\left(n\right)}=\frac{1}{D}\underset{l=(d-1)N/D+1}{\overset{\mathrm{dN}/D}{\mathrm{\&Sigma;}}}{\left|{\hat{x}}_l^{\left(n\right)}\right|}^2$

对于第d个区间段内的

采用如下的阈值迭代：

给定观测矩阵Φ、x的初始估计值α的初始估计值

的初始估计值

正数d、充分小的正常数ε，观测数据y，

$\left\{\begin{array}{c}{\hat{x}}^{(n+1)}={({\mathrm{\&Phi;}}^H\mathrm{\&Phi;}+\mathrm{\&alpha;}{Q}^{\left(n\right)})}^{-1}{\mathrm{\&Phi;}}^{H}y\\ {\widehat{\mathrm{\&alpha;}}}^{(n+1)}={\left|\right|y-\mathrm{\&Phi;}{\hat{x}}^{(n+1)}\left|\right|}_2^2\\ {\mathrm{\&delta;}}_d^{(n+1)}={\frac{1}{D}\underset{/=(d-1)N/D+1}{\overset{\mathrm{dN}/D}{\mathrm{\&Sigma;}}}\left|{\hat{x}}_l^{(n+1)}\right|}^2\\ {\hat{q}}_{(d-1)N/D+1}^{(n+1)}=\&CenterDot;\&CenterDot;\&CenterDot;={\hat{q}}_{\mathrm{dN}/D}^{(n+1)}=\left\{\begin{array}{cc}{\mathrm{\&delta;}}_d^{(n+1)}-\mathrm{\&rho;}+\mathrm{\&epsiv;},& \mathrm{if}{\mathrm{\&delta;}}_d^{(n+1)}\&GreaterEqual;\mathrm{\&rho;}\\ \mathrm{\&epsiv;},& \mathrm{if}{\mathrm{\&delta;}}_d^{(n+1)}<\mathrm{\&rho;}\end{array}\right.\\ {\hat{Q}}^{(n+1)}=\mathrm{diag}(1/{\hat{q}}_1^{(n+1)},\&CenterDot;\&CenterDot;\&CenterDot;,1/{\hat{q}}_N^{(n+1)})\end{array}\right.$

其中ε为一个很小的常数，其作用时避免病态对角矩阵

的出项；ρ为迭代阈值，通常可以设置为

式中q采用如下的软阈值计算方法，如图5所示。

$q=\left\{\begin{array}{cc}f-\mathrm{\&rho;}+\mathrm{\&epsiv;},& \mathrm{if\; f}\&GreaterEqual;\mathrm{\&rho;}\\ \mathrm{\&epsiv;},& \mathrm{if\; f}<\mathrm{\&rho;}\end{array}\right.$

如上的迭代图像重构方法考虑了x的一维块稀疏性，如果结合成像场景的二维稀疏性，在每一次迭代时，平均一个小矩形平面内的若干个分辨单元的散射强度(参见图4(b))，可以得到改进的阈值迭代重构算法，以进一步降低虚警概率和漏警概率。假设图4(b)中的场景重构是逐行进行的，在第n次迭代后第k行信号的第i个元素的估计值为

那么包含该位置信号矩形块(宽度为D，高度为D₁)内的信号平均幅度定义为：

${\mathrm{\&delta;}}_d^{\left(n\right)}=\frac{1}{D{D}_1}\underset{l=(d-1)N/D+1}{\overset{\mathrm{dN}/D}{\mathrm{\&Sigma;}}}{\underset{k=1}{\overset{D_1}{\mathrm{\&Sigma;}}}\left|{\hat{x}}_{l,k}^{\left(n\right)}\right|}^2$

相邻的第d₁(1，2，…，D₁)行的场景重构方法中迭代运算采用以下方法实施：

$\left\{\begin{array}{c}{{\hat{x}}_{:,d_1}}^{(n+1)}={({\mathrm{\&Phi;}}^H\mathrm{\&Phi;}+{\widehat{\mathrm{\&alpha;}}}^{(n+1)}{Q}^{\left(n\right)})}^{-1}{\mathrm{\&Phi;}}^H{y}_{:,d_1}\\ {\widehat{\mathrm{\&alpha;}}}^{(n+1)}={\left|\right|y-\mathrm{\&Phi;}{\hat{x}}^{(n+1)}\left|\right|}_2^2\\ {\mathrm{\&delta;}}_d^{(n+1)}={\frac{1}{D{D}_1}\underset{/=(d-1)N/D+1}{\overset{\mathrm{dN}/D}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{D_1}{\mathrm{\&Sigma;}}}\left|{\hat{x}}_{l,k}^{(n+1)}\right|}^2\\ {\hat{q}}_{(d-1)N/D+1}^{(n+1)}=\&CenterDot;\&CenterDot;\&CenterDot;={\hat{q}}_{\mathrm{dN}/D}^{(n+1)}=\left\{\begin{array}{cc}{\mathrm{\&delta;}}_d^{(n+1)}-\mathrm{\&rho;}+\mathrm{\&epsiv;},& \mathrm{if}{\mathrm{\&delta;}}_d^{(n+1)}\&GreaterEqual;\mathrm{\&rho;}\\ \mathrm{\&epsiv;},& \mathrm{if}{\mathrm{\&delta;}}_d^{(n+1)}<\mathrm{\&rho;}\end{array}\right.\\ {\hat{Q}}^{(n+1)}=\mathrm{diag}(1/{\hat{q}}_1^{(n+1)},\&CenterDot;\&CenterDot;\&CenterDot;,1/{\hat{q}}_N^{(n+1)})\end{array}\right.$

为了验证本发明方法的有效性，对如图2的稀疏场景进行降采样稀疏观测，观测数据量为Nyquist采样数据量的25％，图6给出了采用本专利方法得到的稀疏场景图2的重构结果，对比图2和图6，特别是其中用矩形标注的区域1和2，图6中这两块区域内目标(船只)的方位(纵向)副瓣较低，另外对比还可以看出图6中重构图像中的目标轮廓较为清晰，目标与杂波的对比度较高，虚假点较少，基本满足雷达图像理解和目标识别的要求。

以上为本发明方法所提的基于块稀疏的稀疏微波成像方法的主要内容，以上内容并不局限于使用所提出的阈值迭代方法，也可适用于其他阈值迭代方法。

Claims (6)

1.一种基于块稀疏的稀疏微波成像方法，其特征在于，包括以下步骤：

(1)利用成像场景的稀疏性，对场景进行稀疏观测，得到降采样率的原始回波数据；

(2)利用稀疏优化的方法，基于稀疏观测到的原始回波数据重构微波成像场景；

(3)基于场景的块稀疏性，实现快稀疏场景的稀疏重构：利用实际成像场景的分块稀疏特性，将稀疏微波成像场景在距离向和方位向分割为许多小的矩形块，分别计算这些矩形块内的所有分辨单元的后向散射系数的某种平均结果：算术平均、几何平均、以及由它们构成向量的2-范数并以此计算迭代场景重建算法中阈值。

2.根据权利要求1所述的稀疏微波成像方法，其特征在于，所述步骤(1)中的稀疏观测，包括：

假设稀疏的成像场景按照分辨单元大小被离散化表示为x_l，k，其中l，k分别代表距离和方位下标，雷达平台向成像场景发射宽带脉冲信号s(t)；给定雷达脉冲信号的发射时刻(慢时间)、及对每个脉冲的采样时刻(慢时间)，含有噪声n的原始离散回波数据y和成像场景x满足关系式：y＝Φ·x+n，其中Φ定义为雷达观测矩阵；

s(t)是线性调频信号，或是宽带波形编码信号，且雷达多次发射的宽带脉冲信号相同或不同；

以最小化观测矩阵Φ的所有不同两列的互相关系数的最大值为准则确定发射时刻、采样时刻、及选择波形s(t)。

3.根据权利要求1所述的稀疏微波成像方法，其特征在于，所述步骤(2)中的场景重构方法，包括：

3a、给定雷达观测矩阵Φ、观测噪声的功率α的任意初始估计值

对角元素均为非负数的对角矩阵Q的任意初始估计值

观测数据y以及一个充分小的正数ε(如ε＝10^-6)；

3b、令n＝0，1，2，3，…，执行如下迭代步骤：

3b1、计算： ${\hat{x}}^{(n+1)}={({\mathrm{\&Phi;}}^H\mathrm{\&Phi;}+{\widehat{\mathrm{\&alpha;}}}^{\left(n\right)}{\hat{Q}}^{\left(n\right)})}^{-1}{\mathrm{\&Phi;}}^{H}y;$

3b2、更新： ${\widehat{\mathrm{\&alpha;}}}^{(n+1)}={\left|\right|y-\mathrm{\&Phi;}{\hat{x}}^{(n+1)}\left|\right|}_2^2,$ 定义 ${\hat{q}}_i^{(n+1)}={\left|{\hat{x}}_i^{(n+1)}\right|}^2+\mathrm{\&epsiv;},$

更新 ${\hat{Q}}^{(n+1)}=\mathrm{diag}(1/{\hat{q}}_1^{(n+1)},\&CenterDot;\&CenterDot;\&CenterDot;,1/{\hat{q}}_N^{(n+1)});$

3c、判断：如果连续两次迭代估计值

小于预先给定的常数或者迭代次数达到人为设定的上限，迭代结束；否则n值自动加1，重复上述计算3b1、更新3b2两个步骤。

4.根据权利要求1或3所述的稀疏微波成像方法，其特征在于，所述步骤3b中场景重构的迭代方法，在每次迭代计算的过程中利用场景的一维或两维分块稀疏特性修改迭代步骤。

5.根据权利要求4所述的稀疏微波成像方法，其特征在于，所述使用场景的一维分块稀疏特性使用如下迭代步骤，包括以下步骤：

(1)给个观测矩阵Φ、α的初始估计值

Q的初始估计值

一个假定的正常数ρ、观测数据y并令n＝0；

(2)计算 ${\hat{x}}^{(n+1)}={({\mathrm{\&Phi;}}^H\mathrm{\&Phi;}+{\widehat{\mathrm{\&alpha;}}}^{\left(n\right)}{\hat{Q}}^{\left(n\right)})}^{-1}{\mathrm{\&Phi;}}^{H}y;$

(3)更新 ${\widehat{\mathrm{\&alpha;}}}^{(n+1)}={\left|\right|y-\mathrm{\&Phi;}{\hat{x}}^{(n+1)}\left|\right|}_2^2,$

${\mathrm{\&delta;}}_d^{(n+1)}=\left(\frac{1}{D}\underset{l=(d-1)N/D+1}{\overset{\mathrm{dN}/D}{\mathrm{\&Sigma;}}}{\left|{\hat{x}}_l^{(n+1)}\right|}^2\right),$

${\hat{q}}_{(d-1)N/D+1}^{(n+1)}=\&CenterDot;\&CenterDot;\&CenterDot;={\hat{q}}_{\mathrm{dN}/D}^{(n+1)}=\left\{\begin{array}{cc}{\mathrm{\&delta;}}_d^{(n+1)}-\mathrm{\&rho;}+\mathrm{\&epsiv;},& \mathrm{if}{\mathrm{\&delta;}}_d^{(n+1)}\&GreaterEqual;\mathrm{\&rho;}\\ \mathrm{\&epsiv;},& \mathrm{if}{\mathrm{\&delta;}}_d^{(n+1)}<\mathrm{\&rho;}\end{array}\right.,$

${\hat{Q}}^{(n+1)}=\mathrm{diag}(1/{\hat{q}}_1^{(n+1)},\&CenterDot;\&CenterDot;\&CenterDot;,1/{\hat{q}}_N^{(n+1)});$

(4)如果连续两次迭代估计值

小于预先给定的常数或者迭代次数达到人为设定的上限，迭代结束；否则n值自动加1，重复步骤(2)和(3)。

6.根据权利要求1或3所述的稀疏微波成像方法，其特征在于，采用基于场景块稀疏特性的阈值迭代类算法、或贪心类的重构方法，包括OMP算法、CoSaMP算法、TWIST算法、FPC算法、或Lq正则化方法其中之一，及它们的组合。

Images