一种适用于常规雷达低角跟踪的改进算法

沈福民

摘要: 本文主要讨论在常规雷达低角跟踪中运用波束空间域最大似然(BDML)估计的实现技术。为了改善目标方向估计的性能,文中介绍了频率分集(多个工作频率),波束空间延拓,利用先验信息,接收数据的预处理和估计数据的滤波处理等技术。最后还给出了雷达对目标跟踪的计算机模拟结果。

关键词:

波束空间域; 最大似然; 雷达低角跟踪

Abstract: This paper proposes a BDML estimator used for low angle tracking in conventional radar. In order to improve estimation performance for target direction, some signal processing methods are introduced, such as frequency diversity (multi-frequency processing), beam-space expansion, utilization of information a priari, snapshot data preprocessing and data filtering after BDML estimation etc. Finally, computer simulations are provided for target tracking.

1. 引言

合成孔径雷达(Synthetic Aperture Radar , SAR)是最早提出并投入实用的成像雷达。合成孔径雷达可实现全天时、全天候的区域监测成像，且对植被等介质具有穿透能力，在灾害评估、环境监测和军事等领域得到了广泛的应用^[1-3]。

由雷达分辨率理论和奈奎斯特采样定理可知，合成孔径雷达系统性能的提高通常伴随着采样数据量的显著增加，这给系统设计和实现带来困难。近年来，随着压缩感知^[4]和相关理论的发展，稀疏信号处理在SAR成像中得到广泛应用，这表明当回波信号具有稀疏性或可压缩性时，便能够以远低于奈奎斯特的采样频率，用较少的观测数据高概率地恢复出原信号。现有的压缩感知SAR成像模型中，大致分为两大类，其一，单独考虑距离向和方位向的稀疏模型，然后与匹配滤波相结合进行成像。例如文献[5]仅在距离向实现了压缩感知技术，文献[6]仅在方位向实现了压缩感知技术，上述方法都可以实现SAR成像，但也存在一些问题。首先，仍然需要处理较大的数据量且没有充分利用空间的稀疏性；其次，成像过程中用到的匹配滤波方式导致很大的局限性。其二，在距离向和方位向同时构建压缩感知模型进行求解^[7-9]，其中文献[7]将接收到的雷达回波和2维后向散射系数矩阵叠加成非常大的列向量，然后进行方位向距离向压缩感知处理，虽然比仅对距离向和方位向处理更进一步减少了采样数据量，但求解过程需要对大的列向量进行处理，大幅增加大数据运算量与信号的重构时间。针对出现的问题，文献[8,9]将接收到的回波与后向散射系数保持为2维矩阵，而不是堆叠成向量，因此可以在降低采样率的同时很好地重建目标信号。不管是求解第1类成像模型还是第2类成像模型，都需要用到重构算法，最终目的都是解决l₀范数最小化问题。传统的压缩感知重构算法包括凸优化算法、贪婪算法和统计优化算法3大类，例如文献[10]实现了基于随机孔径贝叶斯压缩感知理论的高分辨合成孔径雷达成像，这里用到了统计优化类算法，文献[11]将隶属于贪婪算法类的正交匹配追踪(Orthogonal Matching Pursuit, OMP)算法应用到合成孔径雷达成像中，文献[12]将求解凸优化问题的快速迭代收缩算法(Fast Iterative Shrinkage Thresholding Algorithm, FISTA)应用到压缩感知ISAR成像的过程中。另有一些学者将非凸函数族中的平滑l₀算法应用到合成孔径雷达方位向稀疏模型求解过程中。但这些重构算法都存在重构精度不高、抗噪性能差等问题。

针对上述问题，本文在文献[8,9]构建的压缩感知SAR成像模型的基础上，将非凸非平滑优化算法引入到压缩感知SAR成像模型的求解问题中，提出一种基于迭代近端投影的2维欠采样合成孔径雷达成像重建算法。首先将压缩感知SAR成像模型转化为近端函数^[13,14]稀疏信号模型，然后采用SCAD判罚函数^[15]来进一步求解该模型，进而重构出图像。仿真与实测数据处理结果表明，本文方法具有很好的重构性能，具有较好的工程应用前景。

2. SAR回波信号2维稀疏表示模型

图1为SAR成像几何模型，假设雷达平台沿 $y$ 轴以恒定的速度 $v$ 前行，若雷达发射线性调频(Linearly Frequency Modulated, LFM)信号，则理想点目标 $P(x,y)$ 的基带回波为

图 1 合成孔径雷达成像几何模型

下载: 全尺寸图片幻灯片

$\begin{split} s({t_{\rm{r}}},{t_{\rm{a}}};r) =& \iint\limits_D {{\boldsymbol{g}}\left( {x,y} \right){a_{\rm{a}}}({t_{\rm{a}}}){a_{\rm{r}}}\left({t_{\rm{r}}} - \frac{{2R({t_{\rm{a}}};r)}}{c}\right)} \\ & \cdot \exp\left[{\rm{j}}\pi \gamma {\left({t_{\rm{r}}} - \frac{{2R({t_{\rm{a}}};r)}}{c}\right)^2}\right]\\ & \cdot\exp \left[ - {\rm{j}}\frac{{4\pi }}{\lambda }R({t_{\text{a}}};r)\right]{\text{d}}x{\text{d}}y \\[-15pt] \end{split}$

(1)

其中， $D = \left\{ {\left( {x,y} \right)|{x_{\min }} \le x \le {x_{\max }},{y_{\min }} \le y \le {y_{\max }}} \right\}$ , $R({t_{\text{a}}};r)$ 为 ${t_{\text{a}}}$ 时刻雷达与目标间的距离，其表达式为

$R({t_{\rm{a}}};r){\text{ = }}\sqrt {{r^2} + {v^2} \cdot {{({t_{\rm{a}}} - {t_{{{\rm{a}}_0}}})}^2}} \approx r + \frac{{{v^2}}}{{2r}}{\left({t_{\rm{a}}} - {t_{{{\rm{a}}_0}}}\right)^2}$

(2)

其中， ${t_{\rm{r}}}$ 和 ${t_{\rm{a}}}$ 分别为距离向快时间和方位向慢时间； $r$ 为点目标 $P(x,y)$ 到雷达航线的垂直距离； ${\boldsymbol{g}}\left( {x,y} \right)$ 为目标的后向散射系数矩阵； $a{}_{\rm{a}}( \cdot )$ 为方位向窗函数； ${a_{\rm{r}}}({t_{\rm{r}}}) = {\rm{rect}}\left(\dfrac{{{t_{\rm{r}}}}}{{{T_{\rm{p}}}}}\right)$ 为发射信号包络， ${T_{\rm{p}}}$ 为脉冲持续时间； $c$ 为光速； $\gamma$ 为发射的LFM信号的线性调频率； $\lambda$ 为波长。根据文献[8]和文献[9]知，理想点目标的回波经采样后可以表示为

${\boldsymbol{S}} = {{\boldsymbol{\varPhi}} _{\rm{r}}}{\boldsymbol{G}}{{\boldsymbol{\varPhi}} _{\rm{a}}}$

(3)

其中， ${\boldsymbol{S}} \in {\mathbb{C}^{M \times N}}$ 表示均匀采样回波信号； ${\boldsymbol{G}} \in {\mathbb{C}^{M \times N}}$ 表示观测场景目标点的后向散射矩阵；矩阵 ${{\boldsymbol{\varPhi}} _{\text{r}}} \in {\mathbb{C}^{M \times M}}$ 为距离向的卷积矩阵，其列是距离向参考函数的共轭与快时间域中采样数的卷积，其中距离向参考函数可以表示为 ${s_{\rm{r}}}({t_{\rm{r}}}) = {a_{\rm{r}}}\left( {{t_{\rm{r}}}} \right)\exp ({\rm{j}}\pi \gamma t_r^2)$ ；矩阵 ${{\boldsymbol{\varPhi}} _{\text{a}}} \in {\mathbb{C}^{N \times N}}$ 为方位向的卷积矩阵，其行是方位向参考函数的共轭与慢时间域中采样数的卷积，其中方位向参考函数可以表示为 ${s_{\rm{a}}}({t_{\rm{a}}}) = {a_{\rm{a}}}\left( {{t_{\rm{a}}}} \right)\exp \left( {{\rm{j}}\pi {k_{\rm{a}}}t_{\rm{a}}^2} \right)$ ，其中 ${k_{\rm{a}}} = \dfrac{{2{v^2}r}}{\lambda }$ 为多普勒调频率； $\mathbb{C}$ 表示复数域； $M$ 为距离向采样点数； $N$ 为方位向采样点数。

对飞机、桥梁等具有空域稀疏性的场景，首先，构建稀疏随机测量矩阵来减少2维采样量，随机测量矩阵构建方式为随机选取方位向和距离向部分采样点信号，如图2所示，其中灰色部分代表采样点信号。此时均匀采样回波信号 ${\boldsymbol{S}} \in {\mathbb{C}^{M \times N}}$ 变为随机采样的回波信号 ${{\boldsymbol{S}}_{2{\rm{CS}}}} \in {\mathbb{C}^{P \times Q}}$ ， $P\left( {P \le M} \right)$ 为距离向随机采样点数， $Q\left( {Q \le N} \right)$ 为方位向随机采样点数。与此同时，距离向卷积矩阵 ${{\boldsymbol{\varPhi}} _{\rm{r}}}$ 变为 ${{\boldsymbol{\varPhi}} _{{\rm{RCS}}}} \in {\mathbb{C}^{P \times M}}$ ，方位向卷积矩阵 ${{\boldsymbol{\varPhi}} _{\rm{a}}}$ 变为 ${{\boldsymbol{\varPhi}} _{{\rm{ACS}}}} \in {\mathbb{C}^{N \times Q}}$ ，通过上述整理可得，2维随机降采样合成孔径雷达回波信号的欠定方程表示为

图 2 SAR回波数据随机采样

下载: 全尺寸图片幻灯片

${{\boldsymbol{S}}_{2{\rm{CS}}}}{\text{ = }}{{\boldsymbol{\varPhi}} _{{\rm{RCS}}}}{\boldsymbol{G}}{{\boldsymbol{\varPhi}} _{{\rm{ACS}}}}$

(4)

考虑到现实场景中加性噪声的影响，式(4)可以表示为

${{\boldsymbol{S}}_{2{\rm{ZCS}}}}{\text{ = }}{{\boldsymbol{\varPhi}} _{{\rm{RCS}}}}{\boldsymbol{G}}{{\boldsymbol{\varPhi}} _{{\rm{ACS}}}} + {{\boldsymbol{N}}_Z}$

(5)

其中， ${{\boldsymbol{S}}_{2{\rm{CS}}}} \in {\mathbb{C}^{P \times Q}},{{\boldsymbol{S}}_{2{\rm{ZCS}}}} \in {\mathbb{C}^{P \times Q}}$ 为随机采样的回波信号， ${{\boldsymbol{\varPhi}} _{{\rm{RCS}}}} \in {\mathbb{C}^{P \times M}}$ 为距离维观测矩阵， ${{\boldsymbol{\varPhi}} _{{\rm{ACS}}}} \in {\mathbb{C}^{N \times Q}}$ 为方位维观测矩阵， ${\boldsymbol{G}} \in {\mathbb{C}^{M \times N}}$ 为待重构的目标图像， ${{\boldsymbol{N}}_Z} \in {\mathbb{C}^{P \times Q}}$ 为复高斯白噪声信号。

3. 基于迭代近端投影的合成孔径雷达成像方法

由上述章节描述可知，式(5)所表示的信号模型是2维卷积模型，它可以用两个连续的1维方程表示。首先，令 ${\boldsymbol{X}} \in {\mathbb{C}^{M \times Q}}$ 作为方位维压缩数据。其表达式为

${\boldsymbol{X}} = {\boldsymbol{G}}{{\boldsymbol{\varPhi}} _{{\rm{ACS}}}}$

(6)

则式(5)可以写为

${{\boldsymbol{S}}_{{\text{2ZCS}}}} = {{\boldsymbol{\varPhi}} _{{\rm{RCS}}}}{\boldsymbol{X}} + {{\boldsymbol{N}}_Z}$

(7)

式(6)表示方位向压缩感知模型，式(7)表示距离向压缩感知模型。

求解式(7)和式(6)的过程如下，首先，求解欠定方程式(7)得到距离向压缩后的信号，其次，将接收到的回波 ${{\boldsymbol{S}}_{{\text{2ZCS}}}}$ 分解为列向量的形式，其表达式为

${{\boldsymbol{S}}_{{\text{2ZCS}}}} = \left[ {{s_1},{s_2}, \cdots ,{s_Q}} \right]$

(8)

同时， ${\boldsymbol{X}}$ 也可以分解为列向量形式

${\boldsymbol{X}} = \left[ {{{\boldsymbol{x}}_1},{{\boldsymbol{x}}_2}, \cdots ,{{\boldsymbol{x}}_Q}} \right]$

(9)

那么，式(7)可以表示为

${{\boldsymbol{s}}_q} = {{\boldsymbol{\varPhi}} _{{\rm{RCS}}}}{{\boldsymbol{x}}_q}{\kern 1pt} {\kern 1pt} + {{\boldsymbol{n}}_z},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} q = 1,2, \cdots ,Q$

(10)

其中， ${{\boldsymbol{s}}_q} \in {\mathbb{C}^P}$ , ${{\boldsymbol{x}}_q} \in {\mathbb{C}^M}$ , ${{\boldsymbol{n}}_z} \in {\mathbb{C}^P}$ , $Q$ 为方位向接收到的脉冲数。由于式(10)是欠定方程，其对应的优化问题为式(11)的拉格朗日(Lagrangian)形式^[16]，

$\mathop {\min }\limits_{{x_q}} {\left\| {{x_q}} \right\|_0}, {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\text{s}}{\text{.t}}{\text{.}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\left\| {{s_q} - {{\boldsymbol{\varPhi}} _{{\rm{RCS}}}}{x_q}} \right\|_2} \le \xi$

(11)

其中， ${\left\| \cdot \right\|_0}$ 是l₀范数，它表示向量 ${x_q}$ 中非零元素个数； ${\Vert ·\Vert }_{2}$ 是l₂范数； $\xi$ 是一个与噪声水平有关的参数。众所周知，求解式(11)最小化问题是一个非凸非平滑优化问题，传统基于压缩感知的SAR成像方法是在一定条件下将l₀范数最小化的非凸非平滑优化问题转化为凸或平滑问题进行求解，这会使稀疏模型不可避免地存在一定误差，进而导致SAR成像效果不佳。而本文将能够解决非凸非平滑问题的迭代近端投影的方法应用到合成孔径雷达成像中，具体过程如下。由文献[13,14]知，为求得式(10)稀疏解，可将其转化为近端函数优化问题，其稀疏信号表示模型为

$\mathop {\min }\limits_{{x_q},z} F\left( z \right) + {\delta _\mathcal{A}}_{_\varepsilon }\left( {{x_q}} \right),{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\text{s}}{\text{.t}}{\text{.}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} z = {x_q}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt}$

(12)

其中， $F\left( \cdot \right)$ 是一个类似l₀范数的非凸非平滑函数， $z$ 是一个列向量的辅助变量， ${\delta _{{\mathcal{A}_\varepsilon }}}$ 定义为 ${\mathcal{A}_\varepsilon } \triangleq \{ {x_q}:{\left\| {{s_q} - {{\boldsymbol{\varPhi}} _{{\rm{RCS}}}}{{\boldsymbol{x}}_q}} \right\|_2} \le \varepsilon \}$ 可行集的指示函数^[13]，即

${\delta _{{\mathcal{A}_\varepsilon }}} = \left\{ {\begin{array}{*{20}{c}} {0,}&{{{\boldsymbol{x}}_q} \in {\mathcal{A}_\varepsilon }} \\ { + \infty ,}&{{{\boldsymbol{x}}_q} \notin {\mathcal{A}_\varepsilon }} \end{array}} \right.$

(13)

其中， ${\Vert ·\Vert }_{2}$ 表示l₂范数，式(12)的稀疏信号表示模型中引入惩罚函数方法，即

$\mathop {\min }\limits_{{x_q},z} F\left( z \right) + {\delta _\mathcal{A}}_{_\varepsilon }\left( {{x_q}} \right){\kern 1pt} + \frac{1}{{2\alpha }}\left\| {{x_q} - z} \right\|{\kern 1pt} _2^2{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt}$

(14)

其中， $\alpha > 0$ 为惩罚因子，通过交替最小化方法^[17]解决式(14)所表示问题。表达式为

$\left. \begin{aligned} & {{{\boldsymbol{z}}_{k + 1}} = \arg {{\min }_z}\left\{ \alpha F\left( {\boldsymbol{z}} \right) + \frac{1}{2}\left\| {{\boldsymbol{z}} - {{\boldsymbol{x}}_{q,z}}} \right\|_2^2\right\} } \\ & {{{\boldsymbol{x}}_{q,z + 1}} = \arg {{\min }_{{x_q}}}\left\{ {\delta _{{\mathcal{A}_\varepsilon }}} + \frac{1}{2}\left\| {{{\boldsymbol{x}}_q} - {{\boldsymbol{z}}_{k + 1}}} \right\|_2^2\right\} } \end{aligned} \right\}$

(15)

根据文献[13]可知，式(15)可以进一步写为式(16)和式(17)

$\left. \begin{aligned} \quad & {{{\boldsymbol{z}}_{k + 1}} = {\text{pro}}{{\text{x}}_{\alpha ,F}}({{\boldsymbol{x}}_{q,k}})} \\ & {{{\boldsymbol{x}}_{q,k + 1}} = {\text{pro}}{{\text{x}}_{{\delta _{{\mathcal{A}_\varepsilon }}}}}({{\boldsymbol{z}}_{k + 1}})} \end{aligned} \right\}$

(16)

$\begin{split} \qquad\;\;\quad {{\boldsymbol{x}}_{q,k + 1}} & = {\text{pro}}{{\text{x}}_{{\delta _{{\mathcal{A}_\varepsilon }}}}}\left( {{\text{pro}}{{\text{x}}_{\alpha ,F}}({{\boldsymbol{x}}_{q,k}})} \right) \\ & = {\mathcal{P}_{{\mathcal{A}_\varepsilon }}}({\text{pro}}{{\text{x}}_{\alpha ,F}}({{\boldsymbol{x}}_{q,k}})) \end{split}$

(17)

其中， ${\text{pro}}{{\text{x}}_{\alpha ,F}}({{\boldsymbol{x}}_{q,k}}) = \mathop {\arg \min }\limits_{{x_q}} \Bigr\{ \dfrac{1}{2}\left\| {{\boldsymbol{z}} - {{\boldsymbol{x}}_{q,k}}} \right\|_2^2 + \alpha F\left( {\boldsymbol{z}} \right)\Bigr\}$ 是 $F\left( z \right)$ 的近端算子，即求解 ${{\boldsymbol{x}}_{q,k}}$ 在满足约束条件下的最小的变量值 $z$ ，其核心思想就是将凸或非凸函数的优化问题转化为易求解的近端映射函数，从而实现近似解。 ${\mathcal{P}_{{\mathcal{A}_\varepsilon }}}$ 代表可行集 ${\mathcal{A}_\varepsilon }$ 的投影。此外，为了提高算法的性能和更好地解决非凸问题存在多个局部极小值的问题，本文利用外推步骤来改善算法性能，即

$\left. \begin{aligned} & {{{\bar {\boldsymbol{x}}}_{q,k}} = {{\boldsymbol{x}}_{q,k}} + \omega \left( {{{\boldsymbol{x}}_{q,k}} - {{\boldsymbol{x}}_{q,k - 1}}} \right)} \\ & {{{\bar {\boldsymbol{x}}}_{q,k + 1}} = {\mathcal{P}_{{\mathcal{A}_\varepsilon }}}\left( {{\text{pro}}{{\text{x}}_{\alpha ,F}}({{\bar {\boldsymbol{x}}}_{q,k}})} \right)} \end{aligned} \right\}$

(18)

其中， $\omega \ge 0$ 为权重常数。由文献[13]可知，非平滑函数的近端算子 ${\text{pro}}{{\text{x}}_{\alpha ,F}}({\bar {\boldsymbol{x}}_{q,k}})$ 可由SCAD惩罚函数 $F\left( {{{\bar {\boldsymbol{x}}}_{q,k}}} \right)$ 产生相应的SCAD阈值函数 $T_{\lambda ,a}^{{\rm{scad}}}$ 来获得，该阈值函数可以进一步减小软、硬阈值收缩函数带来的偏差和促进解的稀疏性。因此，为进一步求得更优解，非平滑函数的近端算子 ${\text{pro}}{{\text{x}}_{\alpha ,F}}({\bar {\boldsymbol{x}}_{q,k}})$ 可以表示为

$\begin{split} {\text{pro}}{{\text{x}}_{\alpha ,F}}({\bar x_{q,k}}) &= \mathop {\arg \min }\limits_{{{\bar x}_q}} \left\{ {\frac{1}{2}\left\| {{{\bar x}_q} - {{\bar x}_{q,k}}} \right\|_2^2 + \alpha F\left( {{{\bar x}_q}} \right)} \right\} \\ &= T_\lambda ^{{\rm{scad}}}\left( {{{\bar x}_q}} \right)\\[-10pt] \end{split}$

(19)

其中，惩罚函数 $F\left( {{{\bar {\boldsymbol{x}}}_{q,k}}} \right)$ 的表达式为

$\begin{split} & F\left( {{{\bar {\boldsymbol{x}}}_{q,k}}} \right){\text{ = }}\\ & \left\{ \begin{aligned} & \lambda \left| {{{\bar {\boldsymbol{x}}}_{q,k}}} \right|,\qquad\qquad\qquad\qquad\qquad 0 \le \;\left| {{{\bar {\boldsymbol{x}}}_{q,k}}} \right| < \lambda \hfill \\ & - \frac{{\left( {{{\left| {{{\bar {\boldsymbol{x}}}_{q,k}}} \right|}^2} - 2a\lambda \left| {{{\bar {\boldsymbol{x}}}_{q,k}}} \right| + {\lambda ^2}} \right)}}{{2\left( {a - 1} \right)}},\; \lambda \le \;\left| {{{\bar {\boldsymbol{x}}}_{q,k}}} \right| < a\lambda \hfill \\ & \frac{{\left( {a + 1} \right){\lambda ^2}}}{2},\qquad\qquad\qquad\qquad\;\; \left| {{{\bar {\boldsymbol{x}}}_{q,k}}} \right| \ge a\lambda \end{aligned} \right. \end{split}$

(20)

SCAD函数的近端映射为

$\begin{split} & {T}_{\lambda ,a}^{{\rm{scad}}}\left({\overline{{\boldsymbol{x}}}}_{q}\right)=\\ & \left\{ \begin{aligned} & {\rm{sign}}\left({\overline{{\boldsymbol{x}}}}_{q}\right){\left(\left|{\overline{{\boldsymbol{x}}}}_{q}\right|-\lambda \right)}_+,\qquad \left|{\overline{{\boldsymbol{x}}}}_{q}\right|\le 2\lambda \\ & \frac{\left(a-1\right){\overline{{\boldsymbol{x}}}}_{q}-{\rm{sign}}\left({\overline{{\boldsymbol{x}}}}_{q}\right)a\lambda }{a-2},2\lambda < \left|{\overline{{\boldsymbol{x}}}}_{q}\right|\le a\lambda \\ & {\overline{{\boldsymbol{x}}}}_{q},\qquad \qquad \qquad \qquad \quad\;\; \left|{\overline{{\boldsymbol{x}}}}_{q}\right| > a\lambda \end{aligned}\right. \end{split}$

(21)

其中， $\lambda$ 为调整参数， $a$ 为常量且 $a > 2$ ， ${\rm{sign}}\left( \cdot \right)$ 为符号函数， ${\left(\beta \right)}_{\text{+}}\text{=}\mathrm{max}\left(\beta ，0\right)$ 。通过求解式(18)得到稀疏向量 ${\bar {\boldsymbol{x}}_q}$ ，至此，由 ${\bar{\boldsymbol{ x}}_q}$ 构成的 ${\boldsymbol{X}}$ 就是经过距离向压缩感知得到的1维距离像。然后，再对方位向压缩感知模型式(6)进行求解，为了更好地表述，将式(6)进行转置得

${{\boldsymbol{X}}^{\rm{T}}} = {\boldsymbol{\varPhi}} _{{\rm{ACS}}}^{\rm{T}}{{\boldsymbol{G}}^{\rm{T}}}$

(22)

将 ${{\boldsymbol{X}}^{\rm{T}}}$ 表示为列向量形式即 ${{\boldsymbol{X}}^{\rm{T}}} = \left[ {{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{{\boldsymbol{x}}} }_1},{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{{\boldsymbol{x}}} }_2} \cdot \cdot \cdot ,{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{{\boldsymbol{x}}} }_M}} \right]$ ，将 ${{\boldsymbol{G}}^{\rm{T}}}$ 也表示为列向量形式即 ${{\boldsymbol{G}}^{\rm{T}}} = \left[ {{{\boldsymbol{g}}_{_1}},{{\boldsymbol{g}}_2}, \cdots ,{{\boldsymbol{g}}_M}} \right]$ ，那么欠定方程式(22)可以进一步写为

${\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{{\boldsymbol{x}}} _m} = {\boldsymbol{\varPhi}} _{{\rm{ACS}}}^{\rm{T}}{{\boldsymbol{g}}_{m,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} }}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} m = 1,2, \cdots ,M$

(23)

其中， ${\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{{\boldsymbol{x}}} _m} \in {\mathbb{C}^Q}$ ， ${{\boldsymbol{g}}_m} \in {\mathbb{C}^N}$ ， $M$ 为距离向采样点。对方位向压缩感知模型也使用本文提出的迭代近端投影的方法求解，求解得到稀疏向量 ${{\boldsymbol{g}}_m}$ ，进一步得到2维的SAR图像。

4. 仿真结果与分析

为综合比较分析本文方法的性能，将正交匹配追踪(OMP)算法、平滑l₀范数(Smoothed l₀ norm, SL0)算法、贝叶斯压缩感知(Bayesian Compressive Sensing, BCS)算法与本文提出的迭代近端投影(Iterative Proximal Projection, IPP)算法作为对比；同时，为进一步定量分析成像结果，本文利用目标杂波比(Target to Clutter Ratio, TCR)^[18]和图像熵(Image Entropy, IE)^[18]来定量评价成像质量，用成像时间来评定算法的运算效率。同时，为验证算法在噪声情况下的性能，对回波加入高斯白噪声进行仿真实验。

4.1 不同采样率成像实验

为进一步验证算法的性能，对30个点目标仿真产生的数据在不同采样率和不同信噪比条件下进行仿真实验。主要仿真参数设置如表1所示。为方便分析，设置相同的距离采样率和方位采样率 $\delta$ 。图3给出了目标散射点模型与不同采样率下各算法成像结果图，其中图3(a)为散射点模型，图3(b-e)表示采样率为原采样率1/2时4种算法成像结果。图3(f-i)表示采样率为原采样率1/4时4种算法成像结果。

表 1 雷达仿真参数

参数	数值
雷达信号载频	3 GHz
雷达信号带宽	150 MHz
采样频率	300 MHz
雷达距目标区域中心点	4200 m

下载: 导出CSV

| 显示表格

图 3 目标散射点模型与不同采样率下各算法成像结果

下载: 全尺寸图片幻灯片

由图3可以看出，在采样率为原采样率1/2时，几种算法都可以清晰地实现目标点SAR成像；当采样率为1/4时，OMP算法的成像效果变差，其他算法依然能够清晰成像。

为了进一步比较上述算法的成像质量，利用TCR和IE作为衡量指标。图4(a)、图4(b)分别给出了不同采样率条件下几种算法的TCR及IE对比曲线，从图4(a)可以看出，随着采样率的下降，SL0算法、BCS算法与本文算法的TCR相对保持稳定，而OMP算法在采样率较低时产生较大波动。同时，可以看出，本文算法的TCR一直大于其他算法，说明本文算法成像结果聚焦性能更好，虚假散射点更少。从图4(b)可以看出，本文算法的IE一直低于其他算法，进一步说明本文算法具有很好的成像性能。

图 4 不同采样率下成像性能曲线

下载: 全尺寸图片幻灯片

为进一步验证算法的运算效率，表2给出了采样率同为原采样率的1/2时几种算法成像时间对比，从中进一步得出本文算法成像时间也相对较短。说明本文算法在保证成像质量的前提下能够实现快速成像。

表 2 采样率为原采样率1/2时各算法成像时间

算法	时间(s)
IPP	26.7381
OMP	18.7749
SL0	27.3091
BCS	178.2851

下载: 导出CSV

| 显示表格

4.2 抗噪声性能实验

首先在不同信噪比条件下进行数据仿真实验，将距离向采样率和方位向采样率同时设置为原采样率的1/2，即 $\delta = \dfrac{1}{2}$ ，并在回波信号处添加高斯白噪声使SNR分别为–10 dB和20 dB。成像结果如图5，其中图5(a)—图5(d)为SNR=20 dB时4种算法成像结果，图5(e)—图5(h)为SNR=–10 dB时4种算法成像结果。由图5可以看出，当信噪比为20 dB时，几种算法都可以清晰成像；当采样率为–10 dB时，OMP算法与BCS算法成像结果出现严重失真，SL0算法成像结果虽可以成像，但是出现较多虚假目标，只有本文算法受噪声影响较小，可以较为清晰地成像。

图 5 不同信噪比情况下各算法成像结果比较

下载: 全尺寸图片幻灯片

为了进一步验证算法的鲁棒性，以TCR, IE作为衡量指标。图6(a)、图6(b)分别为不同信噪比条件下各算法的TCR与IE曲线，由图6(a)可以看出，随着信噪比的上升，各算法的TCR逐渐上升，但在较低信噪比时本文算法依然高于其他算法，说明本文算法有一定的抗噪声能力。同时，由图6(b)可以得出，本文算法的IE随着信噪比的上升逐渐下降且一直低于其他算法，进一步说明本文算法具有良好的鲁棒性。综合实验对比可知，本文算法相较于其他算法具有更好的欠采样成像能力。

图 6 不同信噪比情况下成像性能曲线

下载: 全尺寸图片幻灯片

5. 实测数据处理

为了充分验证本文算法对实际场景处理的有效性，本文将采用来自加拿大温哥华RADARSAT-1精细模式2(加拿大航天局版权所有)的实测数据进行实验。主要参数设置如表3所示。图7给出各算法成像结果，其中图7(a)为基于距离多普勒算法的温哥华地区英吉利海湾附近区域的实测场景图。本文将选取图7(a)中红色矩形区域中6个目标点(即6艘货船)进行实验，其原因是这6个目标点相较于海平面来说是强散射点，而又相较于英吉利海湾是稀疏的。图7(b)—图7(e)是4种算法在距离采样率为0.5，方位采样率为0.3时的成像结果。成像结果表明，当采样数据量减少后，本文算法仍能重建出目标，且相较于其他算法能够消除图像中的模糊，进一步验证了本文算法的有效性。

表 3 温哥华场景RADARSAT-1参数

参数	数值
距离带宽	30.3 MHz
距离向采样频率	32.317 MHz
脉冲宽度	30.111 MHz
卫星轨道半径	7186029 m
雷达波长	0.05657 m

下载: 导出CSV

| 显示表格

图 7 各算法成像结果

下载: 全尺寸图片幻灯片

6. 结束语

针对RD等成像算法存在数据量大、采样率高及OMP、SL0和BCS等算法存在成像精度低及抗噪性能差等问题，本文提出一种迭代近端投影的2维欠采样SAR成像方法，首先对SAR回波信号进行分析，构建基于距离向与方位向的稀疏表示模型，在此基础上，利用近端函数优化模型来表示距离向与方位向稀疏模型，并采用SCAD函数获得近端算子来求解该模型。仿真和实际数据结果验证了该方法的有效性。结果表明，该文方法在低采样率以及低信噪比条件下成像结果优于OMP算法、SL0算法与BCS算法。

Zoltowski M, Lee T S. IEEE Trans. on SP, 1991, SP-39(3): 656-671.[2]高世伟,保铮.电子学报,1989,17(1): 1-7.[3]Ziskind I, Wax M. IEEE Trans. on ASSP, 1988, ASSP-36(10): 1553-1560.[4]Wang H, Kaveh M, IEEE Trans, on ASSP, 1985, ASSP-33(4): 823-831.[5]Shan T J, et al. IEEE Trans, on ASSP, 1985, ASSP-33(4): 806-811.

期刊类型引用(0)

其他类型引用(1)

姓名
邮箱
手机号码
标题
留言内容
验证码

留言板