CN103487803A - 迭代压缩模式下机载扫描雷达成像方法 - Google Patents

Abstract

本发明公开了一种迭代压缩模式下机载扫描雷达成像方法，本发明的方法在贝叶斯理论框架下建立方位向超分辨成像的后向模型，该理论下建立的模型能够有效的融合目标的先验信息。在建立超分辨成像后向模型时，遵循雷达成像基本思想对噪声统计特性的假设，采用高斯分布函数描述噪声统计特性。由于目标散射相对于成像背景具有稀疏性，模型中采用拉普拉斯分布函数描述目标散射的先验信息；将雷达超分辨成像精确转化为数学上最大后验概率问题，通过求解后验概率最大时对应的目标信息，然后对目标信息重构，实现雷达超分辨成像。

Description

迭代压缩模式下机载扫描雷达成像方法

技术领域

本发明属于雷达技术领域，具体涉及机载脉冲多普勒雷达前视扫描成像、直升机雷达环扫全方位成像以及航天器交会对接口扫描成像。

背景技术

雷达能够穿透云雨且不受太阳光线的制约，具有全天候、全天时对地面、海洋进行探测和搜索能力,在军事侦察、地形测绘、陆/海追踪与救援、以及海上缉私等领域，发挥不可或缺的重要作用。

机载扫描雷达，是在飞行器飞行过程中，通过机械转动的方式，使雷达波束能够在方位上均匀地扫描被探测区域。利用雷达天线掠过目标顺序的时间关系，通过对回波数据的有效处理，获取目标在空间方位上的反射信息，达到对地面或海洋探测和成像的目的。由于机载扫描雷达仅仅改变原有机载雷达波束对地照射方式，且能够通过信号处理的方法，实现方位高分辨成像。

传统机载扫描雷达成像，一般采样两种方式，其一如文献：Richards M A.Iterativenoncoherent angular superresolution[radar][C]Radar Conference,1988,Proceedings of the1988IEEE National.IEEE,1988:100-105，直接利用实波束回波，进行成像处理，获取探测区域图像。这种方法没有从本质上突破信号波长和天线孔径对雷达天线角分辨率的限制，实用领域相当狭窄；第二种方式是如文献：李悦丽，梁甸农，黄晓涛，一种单脉冲雷达多通道解卷积前视成像方法，信号处理，Vol.23,No.5,2007，具体利用天线和差波束分别获取目标域数据，它主要是利用和差通道幅度信息差异，采用积分累加的方式，在图像域重构目标信息。该方法由于角闪烁的影响，易产生虚假目标，且仅适用于均匀背景下的强反射目标（高于背景10dB以上）或者弱反射目标（低于背景10dB以上）。

在文献Lucy L B.An iterative technique for the rectification of observed distributions[J].Theastronomical journal,1974,79:745.中提出一种基于统计方法建模实现反卷积；该方法假设目标的先验信息在统计学上服从均匀分布，采用最大似然的方法实现反卷积。当该方法应用到实波束雷达成像上时，由于对成像场景中目标散射的先验信息的假设不合理，使得该方法在实波束雷达成像中仅对单点目标效果较好。当场景中存在多目标时，该方法会因为噪声放大，造成目标幅度失真、位置偏移、产生虚假目标等现象这些现象都影响目标的检测、定位和识别性能，给实际应用带来严重后果。

发明内容

本发明针对应用技术领域中重大需求和背景技术存在的缺陷，设计了一种迭代压缩模式下机载扫描雷达超分辨成像方法。

本发明的解决方案是将扫描雷达超分辨成像在数学上表征为使用凸优化方法实现超分辨成像后向模型在最大后验概率标准下的目标散射信息重构。首先，在贝叶斯理论框架下建立雷达超分辨成像的后向模型，该模型中采用高斯分布函数和拉普拉斯分布函数分别描述噪声的统计特性和目标散射的先验信息；然后在概率背景下，将已获取噪声分布函数和目标散射的先验信息精确的表征为目标与回波的后验概率，并对后验概率进行负对数变换，得到对应的目标函数，再通过求解目标函数的全局最小值实现雷达超分辨成像。因为使用拉普拉斯分布函数描述目标散射的先验信息，所以目标函数中含有不可微的l1范数，这使得基于梯度的算法无法实现目标函数全局最小值的实现。为此，本方法构造相应的压缩函数，再使用凸优化理论中的次梯度代替梯度，构造线性的迭代过程，使得迭代的结果在最大后验概率标准下收敛于目标函数全局最小值，最终实现了机载扫描雷达的高分辨成像。

为了方便描述本发明的内容，首先对下列术语进行解释：

术语1：机载扫描雷达

机载扫描雷达，是在飞行器飞行过程中，通过机械转动的方式，使雷达波束在方位上均匀地扫描地面、海洋或者空间目标的一种雷达。

术语2：先验信息

假设样本空间为Θ，在针对样本空间Θ进行试验之前获取样本Θ的经验或者历史资料称为样本空间Θ的先验信息。

术语3：似然函数

似然函数是一种关于统计模型中参数的函数，反映参数变化时函数值的变化。设样本空间Θ服从分布p(Θ,ζ)，ζ表示待估参数，X₁,X₂,…X_n是来自Θ的样本，x₁,x₂,…x_n是X₁,X₂,…X_n的观测值，则样本的联合分布称为似然函数。

术语4：后验概率

假设A,B分别表示样本空间Θ中的两个事件，p(A)表示没有对样本空间Θ训练之前事件A发生的概率，p(B)表示事件B的先验概率，则p(A|B)为给定事件B时，事件A发生的概率，称为事件A发生的后验概率。

本发明的技术方案为：一种迭代压缩模式下机载扫描雷达成像方法，具体包括如下步骤：

步骤一：成像系统参数初始化，

初始化的扫描雷达成像系统参数如下：雷达平台运动速度，记为V；雷达平台高度，记为H；雷达天线波束俯仰角，记为θ；雷达平台初始位置，记为(0,0,H)，脉冲重复时间，记为PRI；成像场景回波方位向采样点数，记为N_a；场景回波距离向采样点数，记为N_r；t时刻雷达平台与场景中目标距离，记为

其中，R₀为初始时刻天线与场景中目标的斜距，

为目标的方位角；

扫描雷达成像区域的方位时间向量记为T_a=[-PRI·N_a/2,-PRI·(N_a/2-1),…,PRI·(N_a/2-1)]；距离时间向量T_r=[-1/f_s·N_r/2,-1/f_s·(N_r/2-1),…,1/f_s·(N_r/2-1)]，其中，f_s为距离向采样率;

步骤二：回波数据沿距离向脉冲压缩，

雷达发射信号记为 $s\left(\mathrm{\τ}\right)=\mathrm{rect}\left(\frac{\mathrm{\τ}}{T_p}\right)\·\exp (j2\mathrm{\π}{f}_c\mathrm{\τ}+\mathrm{j\πk}{\mathrm{\τ}}^2),$ 其中，rect(·)表示矩形函数，T_p表示发射信号脉冲时宽，f_c表示载频，k表示调频斜率，τ表示快时间；

成像目标区域回波可以表示为噪声加上发射信号与目标的卷积结果，对于面目标而言，回波的解析表达式可以写成：

$g_1(\mathrm{\τ},\mathrm{\η})=\underset{(x,y)\∈\mathrm{\Ω}}{\∫\∫}f(x,y)\·{\mathrm{\ω}}_a\left(\frac{\mathrm{\η}-{\mathrm{\η}}_{a_0}}{T_{\mathrm{\β}}}\right)\·s(\mathrm{\τ}-\frac{2\·R(x,y,t)}{c})\mathrm{dxdy}+N_1(\mathrm{\τ},\mathrm{\η})---\left(1\right)$

其中，Ω表示成像区域；(x,y)表示场景中目标的位置；f(x,y)表示点(x,y)处目标散射函数；ω_a表示慢时间域窗函数，表示方位向天线方向图函数的调制；η表示方位向时间，

表示天线方位角初始时刻；T_β表示目标在3dB天线波束宽度驻留时间；c表示电磁波传播速度；N₁(τ,η)表示回波中的噪声；

为了处理需要，将面目标成像区域的回波（1）表示成下列离散形式：

$g_2(\mathrm{\τ},\mathrm{\η})=\underset{(x,y)\∈\mathrm{\Ω}}{\mathrm{\Σ}}f(x,y)\·{\mathrm{\ω}}_a\left(\frac{\mathrm{\η}-{\mathrm{\η}}_0}{T_{\mathrm{\β}}}\right)\·s(\mathrm{\τ}-\frac{2\·R(x,y,t)}{c})+N_2(\mathrm{\τ},\mathrm{\η})---\left(2\right)$

其中，∑表示求和运算，N₂(τ,η)表示（1）中N₁(τ,η)经离散处理后表达形式；

根据距离向参考时间τ_ref和发射信号调频斜率k，构造距离向脉压参考信号

再将s_ref与回波数据g₂(τ,η)进行最大自相关运算，实现回波信号在距离向脉冲压缩，脉冲压缩后的回波数据可表示为：

$g_3(\mathrm{\τ},\mathrm{\η})=\underset{(x,y)\∈\mathrm{\Ω}}{\mathrm{\Σ}}f(x,y)\·{\mathrm{\ω}}_a\left(\frac{\mathrm{\η}-{\mathrm{\η}}_0}{T_{\mathrm{\β}}}\right)\·\sin c\left\{B\right[\mathrm{\τ}-\frac{2\·R(x,y,t)}{c}\left]\right\}+N_3(\mathrm{\τ},\mathrm{\η})---\left(3\right)$

其中，B表示发射信号带宽，N₃(τ,η)表示数据g₂(τ,η)脉压后的噪声；

步骤三：距离徙动校正，

成像场景区域Ω中的点(x,y)在时刻t与雷达平台之间的距离为

对R(x,y,t)进行泰勒展开，并只保留一次项，又因为θ与

较小，所以

cosθ≈1，从而R(x,y,t)≈R₀-Vt。这里可以通过惯导设备测量雷达平台飞行速度V、回波时刻t。在此基础上，消除雷达平台运动造成的距离徙动。具体表现为：对式(3)中的R(x,y,t)进行尺度变换，消除时间t的影响。可得：

$\begin{array}{c}{g}_4(\mathrm{\τ},\mathrm{\η})=\underset{(x,y)\∈\mathrm{\Ω}}{\mathrm{\Σ}}f(x,y)\·{\mathrm{\ω}}_a\left(\frac{\mathrm{\η}-{\mathrm{\η}}_0}{T_{\mathrm{\β}}}\right)\·\exp \{-j4\mathrm{\π}{f}_c\frac{R(x,y,t)}{c}\}\\ \·\sin c\left\{B\right[\mathrm{\τ}-\frac{2\·R_0}{c}\left]\right\}+N_4(\mathrm{\τ},\mathrm{\η})\end{array}---\left(4\right)$

其中，N₄(τ,η)表示数据g₃(τ,η)经距离徙动校正后的噪声；

步骤四：雷达方位向回波信号建模，

方位向回波信号获取的前向模型表示方位向回波数据g₄(τ,η)与天线H、目标散射信息f以及噪声n的关系，记为：

g=Hf+n               （5）

其中，

g_i(i=1,2,…N_r)是长度为N_a的行向量，表示第i个距离单元内所有的方位向目标的回波数据。因此，面目标Ω的回波数据g的矩阵表现形式如下：

$g={\left[\begin{array}{cccc}{g}_{11}& g_{12}& \·\·\·& g_{{1N}_a}\\ \·& \·& & \·\\ \·& \·& \·\·\·& \·\\ \·& \·& & \·\\ g_{i1}& g_{i2}& \·\·\·& g_{{\mathrm{iN}}_a}\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ g_{N_{r}1}& g_{N_{r}2}& \·\·\·& g_{N_r\·N_a}\end{array}\right]}_{N_r\×N_a}---\left(6\right)$

其中，g_i,j(i=1,2,…N_r;j=1,2,…N_a)表示场景Ω中距离向第i个单元，方位向第j个单元的回波数据，H表示雷达天线信息矩阵，n表示噪声向量。

需要说明的是，在实际操作的过程中，针对（6）中的数据采用的是逐行进行超分辨处理。

本步骤中要完成目标信息f后向模型的建模，即已知回波数据g,天线信息矩阵H和噪声n的统计特征的基础上，通过后向模型的建立并求解实现目标信息f的重构。因为雷达天线具有低通效应使得图像高频信息的丢失，如果在频域对前向模型（5）求逆，会产生噪声无限放大，无法实现超分辨成像，因此使用贝叶斯方法建立超分辨成像的后向模型，该模型在考虑噪声统计特性的同时，还有效融合了目标的先验信息，避免噪声放大现象。

雷达成像思想中假设噪声的统计特性是服从零均值高斯分布，从方位向信号前向模型g=Hf+n出发，本方法采用零均值高斯分布描述噪声n的统计特性，并且假设各次接收到的方位向数据g_k(k=1,2,…N_r)中噪声n_k(k=1,2,…N_r)服从独立同分布；又因为前向模型（5）关于场景f和噪声n都是线性的，噪声n的分布特性直接反映了场景中目标f_i与回波数据g_i之间的似然关系，记为： $p\left(g\right|f)=\underset{i=1}{\overset{N_r\·N_a}{\mathrm{\Π}}}p\left(g_i\right|f_i)\∝\exp (-\frac{{\left|\right|g-\mathrm{Hf}\left|\right|}_2^2}{2\·{\mathrm{\σ}}_n^2});$ 其中，表示从1到N_r·N_a的连乘积运算，∝表示正比关系；

表示欧氏距离下的2范数；

表示噪声n的方差。

在贝叶斯理论下建立的超分辨成像后向模型建模时，还需要考虑场景中目标的先验信息。从信号采样数上分析，目标采样数目要远远小于对场景采样的总数，这种情况下，目标相对于回波数据具有稀疏性。拉普拉斯分布函数的图像具有拖尾长的特点，恰能描述稀疏目标的散射特性。因此，本发明使用拉普拉斯分布来描述目标散射的先验信息，并假设各目标f_k(k=1,2,…N_r·N_a)散射的统计特性服从独立同分布，记为： $p\left(f\right)=\underset{k=1}{\overset{N_r\·N_a}{\mathrm{\Π}}}p\left(f_k\right)=\underset{k=1}{\overset{N_r\·N_a}{\mathrm{\Π}}}\frac{1}{\sqrt{2}\mathrm{\β}}\exp (-\frac{\sqrt{2}\·\left|f_k\right|}{\mathrm{\β}});$ 其中，参数β表示尺度参数，决定分布函数拖尾长度，|·|表示取绝对值运算。

将p(g|f)和p(f)代入全概率公式p(f|g)∝p(g|f)·p(f)，并求后验概率最大时对应的f_MAP，经负对数变换后，对应的最大后验概率f_MAP问题转化为相应目标函数的全局最小解。即

$\begin{array}{c}{f}_{\mathrm{MAP}}=\underset{f}{\arg \min }\{-\ln p\left(g\right|f)-\ln p\left(f\right)\}\\ =\underset{f}{\arg \min }\{\frac{1}{2}{\left|\right|\mathrm{Hf}-g\left|\right|}_2^2+\frac{2\sqrt{2}\·{\mathrm{\σ}}_n^2}{\mathrm{\β}}\·{\left|\right|f\left|\right|}_1\}\end{array}---\left(7\right)$

其中，||·||₁表示欧式距离下的l₁范数问题；令

表示正则参数，其作用是用来调节先验信息与似然函数之间的权重。本发明在此之后，所有的λ_r取值均与此处相同。

步骤五：扫描雷达方位向超分辨处理，

在完成步骤四的基础上，本步骤中使用凸优化方法求解式（7）中f_MAP。

由于式（7）中含有不可微的||f||₁，因此传统的基于梯度的求解方法无法实现满足条件的f_MAP。本步骤从||f||₁不可微的特点出发，构造对应目标函数的压缩函数；由于压缩函数的使用，避免了因梯度的计算需要对||f||₁的近似处理，从而提高了超分辨的精度；在此基础上使用次梯度代替梯度构造迭代过程，并将迭代终值作为超分辨处理结果。

首先，引入辅助变量y，对目标函数 $f_{\mathrm{MAP}}=\underset{f}{\arg \min }\{\frac{1}{2}{\left|\right|\mathrm{Hf}-g\left|\right|}_2^2+{\mathrm{\λ}}_r\·{\left|\right|f\left|\right|}_1\}$ 构造如下压缩函数，记为：

${\mathrm{shrink}}_{{\mathrm{\λ}}_r}\left(y\right)=\underset{f}{\arg \min }\frac{1}{2}{\left|\right|f-y\left|\right|}_2^2+{\mathrm{\λ}}_r{\left|\right|f\left|\right|}_1---\left(8\right)$

经计算可得，

${\mathrm{shrink}}_{{\mathrm{\λ}}_r}\left(y\right)=\mathrm{sign}\left(y\right)\·\max \left(\right|y|-{\mathrm{\λ}}_r,0)---\left(9\right)$

其中，sign(·)表示符号函数，max(·,·)表示括号内两个元素取最大值,y是引入的一个辅助变量,在这里没有具体的物理意义，完全是数学上的处理。

为叙述简便，设φ(f)=||f||₁，最大后验f_MAP可以表示成下列形式：

$f_{\mathrm{MAP}}=\underset{f}{\arg \min }\{\frac{1}{2}{\left|\right|\mathrm{Hf}-g\left|\right|}_2^2+{\mathrm{\λ}}_r\mathrm{\φ}\left(f\right)\}---\left(10\right)$

在（10）的基础上，构造迭代的过程如下：

首先对（10）式右边

和φ(f)分别进行梯度和次梯度运算，将计算结果相加并令其等于零，可得：

$0={\left(H\right)}^T(\mathrm{Hf}-g)+{\mathrm{\λ}}_r\·\∂\mathrm{\φ}\left(f\right)---\left(11\right)$

对于任意α>0的常数，将（11）中的(H)^T(Hf-g)移到等式左边，并在等式两边同时加上αI，再乘以f_MAP，可得：

$[\mathrm{\αI}-{\left(H\right)}^T(\mathrm{Hf}-g)]f_{\mathrm{MAP}}=(\mathrm{\αI}+{\mathrm{\λ}}_r\∂\mathrm{\φ})f_{\mathrm{MAP}}---\left(12\right)$

在（12）的基础上，两边同时乘以

可得：

$f_{\mathrm{MAP}}={(\mathrm{\αI}+{\mathrm{\λ}}_r\∂\mathrm{\φ})}^{-1}\·\left\{\right[\mathrm{\αI}-{\left(H\right)}^T(\mathrm{Hf}-g)]\·f_{\mathrm{MAP}}\}---\left(13\right)$

其中，(·)^T表示对括号内元素进行转置操作，表示函数φ(f)的次梯度，(·)^-1表示对括号内元素进行求逆运算，I表示元素全为数字1的列向量；

联立（9）得到下列迭代过程：

$\begin{array}{c}{f}_{k+1}={\mathrm{shrink}}_{\frac{{\mathrm{\λ}}_r}{\mathrm{\α}}}[f_k-\frac{1}{\mathrm{\α}}{\left(H\right)}^T(H\·f_k-g)]\\ =\mathrm{sign}(f_k-\frac{1}{\mathrm{\α}}{\left(H\right)}^T({\mathrm{Hf}}_k-g))\·\max \left(\right|f_k-\frac{1}{\mathrm{\α}}{\left(H\right)}^T({\mathrm{Hf}}_k-g)|-\frac{{\mathrm{\λ}}_r}{\mathrm{\α}},0)\end{array}---\left(14\right)$

其中，f_k与f_k+1分别表示第k次和第k+1次迭代结果，当迭代次数等于预先设定的次数时，迭代结束，将迭代的结果作为最终的扫描雷达成像结果。

本发明的有益效果：本发明方法的特点在于通过数学方法处理扫描雷达方位向信号，实现超分辨成像。本发明在贝叶斯理论框架下建立方位向超分辨成像的后向模型，该理论下建立的模型能够有效的融合目标的先验信息。本发明在建立超分辨成像后向模型时，遵循雷达成像基本思想对噪声统计特性的假设，采用高斯分布函数描述噪声统计特性，由于目标散射相对于成像背景具有稀疏性，模型中采用拉普拉斯分布函数描述目标散射的先验信息；将雷达超分辨成像精确转化为数学上最大后验概率问题，通过求解后验概率最大时对应的目标信息，进而对目标信息重构，实现雷达超分辨成像。

由于采用拉普拉斯分布函数描述目标散射的先验信息，使得目标函数含有不可微的l1范数，使得基于梯度的算法无法实现目标函数全局最小值的求解，因此本发明的方法构造相应的压缩函数，采用凸优化理论中的次梯度代替梯度构造迭代过程，避免了对目标函数进行近似处理，从而提高了超分辨成像精度。本发明方法的特点是采用线性方法求解线性的超分辨成像后向模型，不仅突破了传统合成孔径雷达成像方法无法实现雷达平台正前方区域分辨成像的技术瓶颈，而且在理论上保证了卷积反演问题由病态向良态的转化，使得雷达在超分辨成像过程中对噪声有效压制，防止目标畸变现象的出现。

附图说明

图1是发明方法的流程示意图。

图2是本发明具体实施方式采用的扫描雷达成像系统结构图。

图3是本发明具体实施时采用的仿真目标场景布置图。

图4是距离向脉压后的数据加入SNR=0dB高斯白噪声图形。

图5是脉压后数据加入SNR=0dB高斯白噪声沿方位向的剖面图。

图6是实施本发明中方法对图3中10个点目标进行成像的结果示意图。

图7是对应图6的成像结果沿方位向的剖面图。

具体实施方式

本发明采用仿真实验来论证，所有的步骤、结论都在Matlab2012上验证正确，以下是对本发明的具体实施作进一步的详细阐述。

步骤一：成像系统参数初始化

本仿真验证采用的成像几何模式如图2所示，系统坐标系以雷达在正下方地面上投影为俯视图坐标原点，雷达平台以速度V沿X轴方向运动，雷达天线以角速度为ω对平台正前方区域进行扫描，Y轴为切航迹方向，Z为垂直地面方向。对应的雷达成像系统参数如表1所示。

表1

本示例采用的成像场景如图3所示，途中圆点为布置与地面上的2×5的点目标，沿Y轴正方向，幅度依次为0.5、0.8、1、0.8、0.5，间隔为2°，这样设置目标散射幅度，有助于仿真验证本发明方法在保持目标幅度上的优势；沿X轴方向间隔为500m，雷达平台初始时刻位置坐标为(0,0,5000m)，XOY平面内目标散射函数记为f(x,y)。t时刻XOY平面内的点(x,y)与雷达平台距离记为R(x,y,t)。图中圆点是布置于地面上的10个点目标。

构造方位时间向量T_a=[-PRI·N_a/2,-PRI·(N_a/2-1),…,PRI·(N_a/2-1)]；其中PRI为发射信号脉冲重复时间，N_a为目标回波方位向采样点数；距离时间向量T_r=[-1/f_s·N_r/2,-1/f_s·(N_r/2-1),…,1/f_s·(N_r/2-1)]，其中，f_s为距离向采样率，N_r为目标回波距离向采样点数。

步骤二：对雷达接收到的回波数据进行距离向脉冲压缩，

使用Matlab仿真平台下，生成成像区域Ω的回波g₂(τ,η)，再根据发射信号调频斜率k和距离向参考时间τ_ref，构造距离向脉压参考信号

与回波数据g₂(τ,η)进行最大自相关运算，从而完成距离向脉冲压缩，压缩后的回波数据表示为g₃(τ,η)。

步骤三：距离徙动校正，

步骤二处理得到的数据g₃(τ,η)中含有目标的斜距距离历史R(x,y,t)。对斜距距离历史R(x,y,t)进行泰勒级数展开，因为由于θ与

较小，所以

cosθ≈1。从而R(x,y,t)≈R₀-Vt，在此基础上，对数据g₃(τ,η)进行尺度变换和消除雷达平台运动造成的距离徙动可得g₄(τ,η)。

为了模拟实际处理中噪声的存在，在数据g₄(τ,η)中加入SNR=0dB的高斯白噪声，相应的结果如图4所示，沿方位向剖面如图5所示。

步骤四：方位向信号建模，

将步骤三处理的结果，表示成矩阵向量的形式：g=Hf+n。将似然函数p(g|f)和先验信息p(f)无近似的表示成场景与方位向回波之间的后验概率p(f|g)。再经反对数变换后，方位向信号处理问题转化为下列目标函数最优解的实现问题。即

$f_{\mathrm{MAP}}=\underset{f}{\arg \min }\{\frac{1}{2}{\left|\right|\mathrm{Hf}-g\left|\right|}_2^2+{\mathrm{\λ}}_r\·{\left|\right|f\left|\right|}_1\}.$

步骤五：扫描雷达方位向超分辨处理，

本步骤是用凸优化方法实现 $f_{\mathrm{MAP}}=\underset{f}{\arg \min }\{\frac{1}{2}{\left|\right|\mathrm{Hf}-g\left|\right|}_2^2+{\mathrm{\λ}}_r\·{\left|\right|f\left|\right|}_1\}.$ 首先构造和目标函数对应的压缩函数： ${\mathrm{shrink}}_{{\mathrm{\λ}}_r}\left(y\right)=\mathrm{sign}\left(y\right)\·\max \left(\right|y|-{\mathrm{\λ}}_r,0);$ 在此基础上可得 $f_{\mathrm{MAP}}={(\mathrm{\αI}+{\mathrm{\λ}}_r\·\∂\mathrm{\φ})}^{-1}\·[\mathrm{\αI}-{\left(H\right)}^T(\mathrm{Hf}-g)]\·f_{\mathrm{MAP}};$

从而构造下列迭代过程：

$\begin{array}{c}{f}_{k+1}={\mathrm{shrink}}_{\frac{{\mathrm{\λ}}_r}{\mathrm{\α}}}[f_k-\frac{1}{\mathrm{\α}}{\left(H\right)}^T(H\·f_k-g)]\\ =\mathrm{sign}(f_k-\frac{1}{\mathrm{\α}}{\left(H\right)}^T({\mathrm{Hf}}_k-g))\·\max \left(\right|f_k-\frac{1}{\mathrm{\α}}{\left(H\right)}^T({\mathrm{Hf}}_k-g)|-\frac{{\mathrm{\λ}}_r}{\mathrm{\α}},0)\end{array}$

并将其在Matlab仿真平台下迭代运算预先设定的次数，将迭代的结果作为最终的扫描雷达成像结果，处理结果如图6所示，对应沿方位向剖面图如图7所示。

通过本发明的示例可以看出，本发明可以实现对扫描雷达前视区域进行高分辨成像。成像结果对于目标的幅度和位置信息、图像分辨率等具有良好的效果。

本发明方法的优点是不仅克服了传统雷达无法对前视区域高分辨成像的问题，而且在保持雷达图像中目标散射强度、位置信息、虚假目标和噪声抑制上具有显著优势；本发明的方法可以应用于物资空投、地球遥感、空间探测等领域。本领域工程技术人员可根据本发明公开的技术做出相应的应用，相关的应用仍在本发明的保护范围之内。

Claims (1)

1.一种迭代压缩模式下机载扫描雷达成像方法，具体包括如下步骤：

步骤一：成像系统参数初始化，

初始化的扫描雷达成像系统参数如下：雷达平台运动速度，记为V；雷达平台高度，记为H；雷达天线波束俯仰角，记为θ；雷达平台初始位置，记为(0,0,H)，脉冲重复时间，记为PRI；成像场景回波方位向采样点数，记为N_a；场景回波距离向采样点数，记为N_r；t时刻雷达平台与场景中目标距离，记为

其中，R₀为初始时刻天线与场景中目标的斜距，

为目标的方位角；

扫描雷达成像区域的方位时间向量记为T_a=[-PRI·N_a/2,-PRI·(N_a/2-1),…,PRI·(N_a/2-1)]；距离时间向量T_r=[-1/f_s·N_r/2,-1/f_s·(N_r/2-1),…,1/f_s·(N_r/2-1)]，其中，f_s为距离向采样率；

步骤二：回波数据沿距离向脉冲压缩，

雷达发射信号记为 $s\left(\mathrm{\τ}\right)=\mathrm{rect}\left(\frac{\mathrm{\τ}}{T_p}\right)\·\exp (j2\mathrm{\π}{f}_c\mathrm{\τ}+\mathrm{j\πk}{\mathrm{\τ}}^2),$ 其中，rect(·)表示矩形函数，T_p表示发射信号脉冲时宽，f_c表示载频，k表示调频斜率，τ表示快时间；

对于面目标而言，回波的解析表达式可以写成：

$g_1(\mathrm{\τ},\mathrm{\η})=\underset{(x,y)\∈\mathrm{\Ω}}{\∫\∫}f(x,y)\·{\mathrm{\ω}}_a\left(\frac{\mathrm{\η}-{\mathrm{\η}}_{a_0}}{T_{\mathrm{\β}}}\right)\·s(\mathrm{\τ}-\frac{2\·R(x,y,t)}{c})\mathrm{dxdy}+N_1(\mathrm{\τ},\mathrm{\η})---\left(1\right)$

其中，Ω表示成像区域；(x,y)表示场景中目标的位置；f(x,y)表示点(x,y)处目标散射函数；ω_a表示慢时间域窗函数，表示方位向天线方向图函数的调制；η表示方位向时间，

表示天线方位角初始时刻；T_β表示目标在3dB天线波束宽度驻留时间；c表示电磁波传播速度；N₁(τ,η)表示回波中的噪声；

将面目标成像区域的回波（1）表示成下列离散形式：

$g_2(\mathrm{\τ},\mathrm{\η})=\underset{(x,y)\∈\mathrm{\Ω}}{\mathrm{\Σ}}f(x,y)\·{\mathrm{\ω}}_a\left(\frac{\mathrm{\η}-{\mathrm{\η}}_0}{T_{\mathrm{\β}}}\right)\·s(\mathrm{\τ}-\frac{2\·R(x,y,t)}{c})+N_2(\mathrm{\τ},\mathrm{\η})---\left(2\right)$

其中，∑表示求和运算，N₂(τ,η)表示（1）中N₁(τ,η)经离散处理后表达形式；

根据距离向参考时间τ_ref和发射信号调频斜率k，构造距离向脉压参考信号

再将s_ref与回波数据g₂(τ,η)进行最大自相关运算，实现回波信号在距离向脉冲压缩，脉冲压缩后的回波数据可表示为：

$g_3(\mathrm{\τ},\mathrm{\η})=\underset{(x,y)\∈\mathrm{\Ω}}{\mathrm{\Σ}}f(x,y)\·{\mathrm{\ω}}_a\left(\frac{\mathrm{\η}-{\mathrm{\η}}_0}{T_{\mathrm{\β}}}\right)\·\sin c\left\{B\right[\mathrm{\τ}-\frac{2\·R(x,y,t)}{c}\left]\right\}+N_3(\mathrm{\τ},\mathrm{\η})---\left(3\right)$

其中，B表示发射信号带宽，N₃(τ,η)表示数据g₂(τ,η)脉压后的噪声；

步骤三：距离徙动校正，

对R(x,y,t)进行泰勒展开，并只保留一次项，从而R(x,y,t)≈R₀-Vt，对(3)中函数中的R(x,y,t)进行尺度变换，消除时间t的影响，可得：

$\begin{array}{c}{g}_4(\mathrm{\τ},\mathrm{\η})=\underset{(x,y)\∈\mathrm{\Ω}}{\mathrm{\Σ}}f(x,y)\·{\mathrm{\ω}}_a\left(\frac{\mathrm{\η}-{\mathrm{\η}}_0}{T_{\mathrm{\β}}}\right)\·\exp \{-j4\mathrm{\π}{f}_c\frac{R(x,y,t)}{c}\}\\ \·\sin c\left\{B\right[\mathrm{\τ}-\frac{2\·R_0}{c}\left]\right\}+N_4(\mathrm{\τ},\mathrm{\η})\end{array}---\left(4\right)$

其中，N₄(τ,η)表示数据g₃(τ,η)经距离徙动校正后的噪声；

步骤四：雷达方位向回波信号建模，

方位向回波信号获取的前向模型表示方位向回波数据g₄(τ,η)与天线H、目标散射信息f以及噪声n的关系，记为：

g=Hf+n          （5）

其中，

g_i(i=1,2,…N_r)是长度为N_a的行向量，表示第i个距离单元内所有的方位向目标的回波数据，面目标Ω的回波数据g的矩阵表现形式如下：

$g={\left[\begin{array}{cccc}{g}_{11}& g_{12}& \·\·\·& g_{{1N}_a}\\ \·& \·& & \·\\ \·& \·& \·\·\·& \·\\ \·& \·& & \·\\ g_{i1}& g_{i2}& \·\·\·& g_{{\mathrm{iN}}_a}\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ g_{N_{r}1}& g_{N_{r}2}& \·\·\·& g_{N_r\·N_a}\end{array}\right]}_{N_r\×N_a}---\left(6\right)$

其中，g_i,j(i=1,2,…N_r;j=1,2,…N_a)表示场景Ω中距离向第i个单元，方位向第j个单元的回波数据，H表示雷达天线信息矩阵，n表示噪声向量；

假设噪声的统计特性服从零均值高斯分布，从方位向信号前向模型g=Hf+n出发，采用零均值高斯分布描述噪声n的统计特性，并且假设各次接收到的方位向数据g_k(k=1,2,…N_r)中噪声n_k(k=1,2,…N_r)服从独立同分布；又因为前向模型（5）关于场景f和噪声n都是线性的，噪声n的分布特性直接反映了场景中目标f_i与回波数据g_i之间的似然关系，记为： $p\left(g\right|f)=\underset{i=1}{\overset{N_r\·N_a}{\mathrm{\Π}}}p\left(g_i\right|f_i)\∝\exp (-\frac{{\left|\right|g-\mathrm{Hf}\left|\right|}_2^2}{2\·{\mathrm{\σ}}_n^2});$ 其中，

表示从1到N_r·N_a的连乘积运算，∝表示正比关系；表示欧氏距离下的2范数；表示噪声n的方差。

使用拉普拉斯分布来描述目标散射的先验信息，并假设各目标f_k(k=1,2,…N_r·N_a)散射的统计特性服从独立同分布，记为： $p\left(f\right)=\underset{k=1}{\overset{N_r\·N_a}{\mathrm{\Π}}}p\left(f_k\right)=\underset{k=1}{\overset{N_r\·N_a}{\mathrm{\Π}}}\frac{1}{\sqrt{2}\mathrm{\β}}\exp (-\frac{\sqrt{2}\·\left|f_k\right|}{\mathrm{\β}});$ 其中，参数β表示尺度参数，决定分布函数拖尾长度，|·|表示取绝对值运算。

将p(g|f)和p(f)代入全概率公式p(f|g)∝p(g|f)·p(f)，并求后验概率最大时对应的f_MAP，经负对数变换后，对应的最大后验概率f_MAP问题转化为相应目标函数的全局最小解，即，

$\begin{array}{c}{f}_{\mathrm{MAP}}=\underset{f}{\arg \min }\{-\ln p\left(g\right|f)-\ln p\left(f\right)\}\\ =\underset{f}{\arg \min }\{\frac{1}{2}{\left|\right|\mathrm{Hf}-g\left|\right|}_2^2+\frac{2\sqrt{2}\·{\mathrm{\σ}}_n^2}{\mathrm{\β}}\·{\left|\right|f\left|\right|}_1\}\end{array}---\left(7\right)$

其中，||·||₁表示欧式距离下的l₁范数问题；令

表示正则参数；

步骤五：扫描雷达方位向超分辨处理，

引入辅助变量y，对目标函数 $f_{\mathrm{MAP}}=\underset{f}{\arg \min }\{\frac{1}{2}{\left|\right|\mathrm{Hf}-g\left|\right|}_2^2+{\mathrm{\λ}}_r\·{\left|\right|f\left|\right|}_1\}$ 构造如下压缩函数，记为：

${\mathrm{shrink}}_{{\mathrm{\λ}}_r}\left(y\right)=\underset{f}{\arg \min }\frac{1}{2}{\left|\right|f-y\left|\right|}_2^2+{\mathrm{\λ}}_r{\left|\right|f\left|\right|}_1---\left(8\right)$

经计算可得，

${\mathrm{shrink}}_{{\mathrm{\λ}}_r}\left(y\right)=\mathrm{sign}\left(y\right)\·\max \left(\right|y|-{\mathrm{\λ}}_r,0)---\left(9\right)$

其中，sign(·)表示符号函数，max(·,·)表示括号内两个元素取最大值；

设φ(f)=||f||₁，最大后验f_MAP表示成下列形式：

$f_{\mathrm{MAP}}=\underset{f}{\arg \min }\{\frac{1}{2}{\left|\right|\mathrm{Hf}-g\left|\right|}_2^2+{\mathrm{\λ}}_r\mathrm{\φ}\left(f\right)\}---\left(10\right)$

在（10）的基础上，构造迭代的过程如下：

首先对（10）式右边

和φ(f)分别进行梯度和次梯度运算，将计算结果相加并令其等于零，可得：

$0={\left(H\right)}^T(\mathrm{Hf}-g)+{\mathrm{\λ}}_r\·\∂\mathrm{\φ}\left(f\right)---\left(11\right)$

对于任意α>0的常数，将（11）中的(H)^T(Hf-g)移到等式左边，并在等式两边同时加上αI，再乘以f_MAP，可得：

$[\mathrm{\αI}-{\left(H\right)}^T(\mathrm{Hf}-g)]f_{\mathrm{MAP}}=(\mathrm{\αI}+{\mathrm{\λ}}_r\∂\mathrm{\φ})f_{\mathrm{MAP}}---\left(12\right)$

在（12）的基础上，两边同时乘以

可得：

$f_{\mathrm{MAP}}={(\mathrm{\αI}+{\mathrm{\λ}}_r\∂\mathrm{\φ})}^{-1}\·\left\{\right[\mathrm{\αI}-{\left(H\right)}^T(\mathrm{Hf}-g)]\·f_{\mathrm{MAP}}\}---\left(13\right)$

其中，(·)^T表示对括号内元素进行转置操作，

表示函数φ(f)的次梯度，(·)^-1表示对括号内元素进行求逆运算，I表示元素全为数字1的列向量；

联立（9）得到下列迭代过程：

$\begin{array}{c}{f}_{k+1}={\mathrm{shrink}}_{\frac{{\mathrm{\λ}}_r}{\mathrm{\α}}}[f_k-\frac{1}{\mathrm{\α}}{\left(H\right)}^T(H\·f_k-g)]\\ =\mathrm{sign}(f_k-\frac{1}{\mathrm{\α}}{\left(H\right)}^T({\mathrm{Hf}}_k-g))\·\max \left(\right|f_k-\frac{1}{\mathrm{\α}}{\left(H\right)}^T({\mathrm{Hf}}_k-g)|-\frac{{\mathrm{\λ}}_r}{\mathrm{\α}},0)\end{array}---\left(14\right)$

其中，f_k与f_k+1分别表示第k次和第k+1次迭代结果，当迭代次数等于预先设定的次数时，迭代结束，将迭代的结果作为最终的扫描雷达成像结果。

Images