一种双基地ISAR稀疏孔径成像方法与流程

本发明涉及雷达信号处理技术领域，具体涉及一种双基地isar稀疏孔径成像方法。

背景技术：

双基地isar是发射站和接收站分置的isar成像系统，不仅具有双基地雷达的“四抗”特性，而且与单基地isar相比，没有几何成像盲区并能获得更多的目标信息，日益成为现代雷达技术研究的热点问题之一。

在对目标的观测过程中，通常需要切换雷达波束对多个目标进行多视角观测，容易形成稀疏孔径。若直接利用传统的rd算法实现双基地isar成像，不仅会产生强烈的副瓣和能量泄漏，而且在低信噪比条件下存在大量干扰噪声，严重影响成像质量。压缩感知理论可以在获得少量数据的情况下，利用重构算法高概率重构出原始信号，由于isar空域稀疏性，故可以将其应用到双基地isar稀疏孔径成像中，以提高成像质量。在重构算法中，基于l1范数约束的凸优化算法虽能有效改善成像质量和抗噪能力，但在求解时，不能进行参数自学习，需要先估算约束参数，且估算误差对成像结果影响较大，有时可能不能获得最优稀疏解，影响成像质量。另外，在求解时涉及回波数据矢量化和实数化操作，导致数据存储量和运算量较大，运算效率有待提高。

技术实现要素：

本发明所要解决的技术问题是提供了一种双基地isar稀疏孔径成像方法，能够在提高成像质量的同时减少运算时间，当数据缺失较多时仍能得到较好的成像结果。

本发明的技术方案如下：

一种双基地isar稀疏孔径成像方法，具体包括如下步骤：

步骤1、建立双基地isar成像回波模型，并对平动补偿后的回波进行多普勒位移补偿，得到全孔径回波数据；

步骤2、构造稀疏基矩阵得到全孔径回波数据的稀疏表示，构造有效数据选择矩阵得到稀疏孔径回波数据，建立基于压缩感知的双基地isar稀疏孔径成像模型；

步骤3、在双基地isar稀疏孔径成像模型的基础上，将整个二维回波数据进行分块处理，并设目标图像各像元服从高斯先验进而建立稀疏贝叶斯模型；所述二维回波数据为稀疏孔径回波数据；

步骤4、利用快速边缘似然函数最大化方法对每块回波数据的稀疏贝叶斯模型进行求解，实现对每块回波数据对应的目标图像重构；

步骤5、将所求的每块回波数据对应的目标图像合成整个二维图像，得到重构的目标图像。

进一步的，所述步骤1具体包括：

设雷达发射的线性调频信号，经包络对齐和相位校正后的双基地isar回波用如下式(1)表示：

其中，fc为载波中心频率，tp为发射信号脉冲宽度，μ为调频斜率，σp为散射点p的信号复幅度，xp和yp分别为散射点p的坐标，θ(tm)和β(tm)分别为成像期间内的旋转角度和双基地角，随慢时间tm变化；表示快时间；c表示波速；

为了避免双基地角时变引起越分辨单元徙动和图像畸变，构造相应的补偿相位进行相位补偿，构造的补偿项如下式(2)：

得到相位补偿后的一维距离像如下式(3)：

其中，σp为第p个散射点的信号复幅度。

进一步的，所述步骤2具体包括：

设全孔径回波信号中共包含l个脉冲视角，累积转角为δθ，构造稀疏基矩阵fall将二维成像场景离散化为n个距离单元和m个多普勒单元，其表达式(4)和表达式(5)如下：

其中，ω表示fall中的元素值，其上角标m的取值为0,1,…,m-1，其下角标l的取值为0,1,…,l-1；

考虑到实际噪声的存在，则双基地isar全孔径回波可稀疏表示为如下式(6)：

sall＝falla+ε0(6)

其中，sall表示经过运动补偿和相位补偿后的全孔径二维回波数据，ε0为噪声，a为需求的目标图像；设有效的回波脉冲为j个，构造有效数据选择矩阵t，得到稀疏孔径回波数据的矩阵表达式(7)为：

s＝tsall+ε＝tfalla+ε＝fa+ε(7)

其中，ε为稀疏孔径回波中的噪声。

进一步的，所述步骤3具体包括：

步骤3-a、对二维回波数据分块化：

将回波矩阵sk×n按距离单元分成小块进行处理，每块数据矩阵中包含n1个距离单元的回波数据，则n个距离单元的回波数据一共可分为h＝[n/n1]块，其中，[]表示向上取整；

步骤3-b、将分块回波数据按列矢量化并从复数形式转化为实数形式，得到每块数据的实数域矢量表达式(8)：

其中，为回波数据矢量；为对角块矩阵；为目标图像矢量；为噪声矢量；

步骤3-c、设噪声和目标图像矢量各像元均服从独立高斯先验建立稀疏先验模型，得到分块稀疏贝叶斯学习模型，推得目标矢量的后验概率也服从均值为μ、方差为∑的高斯分布，其中均值μ和方差∑分别如下表达式(9)和表达式(10)：

其中，表示对角元素由组成的对角矩阵，此时的均值μ即为对应的目标矢量估计值σ为信号复幅度。

进一步的，所述步骤4具体包括：

利用快速边缘似然函数最大化方法对分块稀疏贝叶斯模型进行求解，得到每块回波数据对应的目标图像矢量化重构结果。

进一步的，快速边缘似然函数最大化方法实现的具体步骤包括：

步骤4-a、输入回波矢量和稀疏基对角块矩阵初始化迭代次数n＝1；

步骤4-b、将矩阵中的每一列基作为一个原子，将每个原子与回波矢量做内积，选择内积值最大的原子作为被选中的第一个原子fi，利用式和进行初始化，其中表示的标准差；

步骤4-c、计算∑,μ，引入两个中间变量gi和qi，表达式为：

其中，

步骤4-d、令gi＝αigi/(αi-gi)，qi＝αiqi/(αi-gi)，计算出每个原子对应的gi和qi以及θi＝qi²-gi值，分三种情况进行处理：

①若θi＞0且αi＝+∞，则将fi添加到被选中的原子中；

②若θi＞0且αi＜+∞，则对超参数αi进行更新，有αi＝gi²/(qi²-gi)；

③若θi＜0，则将fi从被选中的原子中删除，并设置αi＝+∞；

步骤4-e、利用更新σ²，得到新的α和σ²后，继续更新∑,μ，并计算新的gi和qi；

步骤4-f、判断是否终止迭代，若满足迭代终止条件则停止迭代，此时的μ值即为目标矢量的估计值进一步的，所述步骤5具体包括：

将所求目标矢量估计值中的实部与相应的虚部相结合变为复数矢量信号，再将其转化为二维矩阵形式，得到分块回波数据对应的二维目标图像再将所求的每块回波数据所求的目标图像合成整个二维图像，即可得到重构的目标图像。

本发明的有益效果如下：

本发明采取了分块处理的思想，在每块图像重构时减少了数据存储量和计算量，另外，相比于传统的稀疏贝叶斯学习求解方法，本发明所涉及的快速算法能够在保证重构质量的同时进一步缩短了运算时间，当孔径缺失数量较多时，利用本发明方法仍能高质量地实现图像重构，得到清晰且聚焦良好的目标图像。

附图说明

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

图2是根据本发明的一个实施例的双基地isar成像几何模型。

图3是根据本发明的一个实施例的回波数据分块处理示意图。

图4是根据本发明的一个实施例的稀疏贝叶斯学习两层稀疏先验示意图。

图5是根据本发明的一个实施例的快速边缘似然函数最大化方法流程示意图。

图6-a是根据本发明的一个实施例的一维稀疏信号重构性能的时域仿真图。

图6-b是根据本发明的一个实施例的一维稀疏信号重构性能的频域仿真图。

图7-a是根据本发明的一个实施例的仿真场景图。

图7-b是根据本发明的一个实施例的目标散射点模型图。

图8-a是根据本发明的一个实施例在成像时间段内双基地角变化曲线图。

图8-b是根据本发明的一个实施例在成像时间段内累积转角变化曲线图。

图9-a是50％数据缺失时的一维距离像、rd成像结果以及本发明方法成像结果。

图9-b是65％数据缺失时的一维距离像、rd成像结果以及本发明方法成像结果。

图9-c是80％数据缺失时的一维距离像、rd成像结果以及本发明方法成像结果。

具体实施方式

为了能够更清楚地理解本发明的上述目的、特征和优点，下面结合附图和具体实施方式对本发明进行进一步的详细描述。需要说明的是，在不冲突的情况下，本申请的实施例及实施例中的特征可以相互结合。

在下面的描述中，阐述了很多具体细节以便于充分理解本发明，但是，本发明还可以采用其他不同于在此描述的其他方式来实施，因此，本发明的保护范围并不受下面公开的具体实施例的限制。

本实施例涉及一种双基地isar稀疏孔径成像方法，其基本思路为：步骤1、对双基地isar目标进行回波建模，并对平动补偿后的回波进一步进行多普勒位移补偿，得到全孔径回波数据，再构造稀疏基矩阵和有效数据选择矩阵；步骤2、建立基于压缩感知的稀疏孔径成像模型；步骤3、将整个二维回波数据进行分块处理，并建立稀疏贝叶斯模型；步骤4、通过快速边缘似然函数最大化方法进行求解；步骤5、将所求的每块回波对应的目标图像合成整个二维图像，得到重构的目标图像。

进一步的，所述步骤1中：

假设雷达发射线性调频信号，经包络对齐和相位校正后的双基地isar回波的可近似为：

其中，fc为载波中心频率，tp为发射信号脉冲宽度，μ为调频斜率，σp为散射点p的信号复幅度，xp和yp分别为散射点p的坐标，θ(tm)和β(tm)分别为成像期间内的旋转角度和双基地角，随慢时间tm变化。

为了避免双基地角时变引起越的分辨单元徙动和图像畸变，需要构造相应的补偿相位进行相位补偿，构造的补偿项表达式为：

得到相位补偿后的一维距离像表达式为：

假设在距离单元(2yp/c)cos(β(tm)/2)内有q个强散射点，则此单元的回波信号可表示为：

其中，aq为第q个散射点的信号复幅度。

进一步的，所述步骤2中：

假设全孔径回波信号中共包含l个脉冲视角，累积转角为δθ，构造稀疏基矩阵fall将二维成像场景离散化为n个距离单元和m个多普勒单元，其中：

考虑到实际噪声的存在，则双基地isar全孔径回波可稀疏表示为：

sall＝falla+ε0

其中，ε0为噪声，a为需求的目标图像。假设有效的回波脉冲为j个，构造观测矩阵t，得到稀疏孔径回波数据：

s＝tsall+ε＝tfalla+ε＝fa+ε

进一步的，所述步骤3中：

将整个二维回波数据进行分块处理，并假设目标图像各像元服从高斯先验建立稀疏贝叶斯模型，主要包括以下步骤：

步骤a，对二维回波数据分块化，即将回波矩阵sk×n按距离单元分成小块进行处理，每块数据矩阵中包含n1个距离单元的回波数据，则n个距离单元的回波数据一共可分为h＝[n/n1]块，[]表示向上取整。

步骤b，将分块回波数据按列矢量化，得到每块数据的矢量表达式：

式中sn、εn和an分别表示第“n”个距离单元对应的回波数据矢量、噪声矢量和目标图像矢量。

将矢量化后的数据由复数转化为实数形式，得到表达式：

其中，re()和im()分别表示实部和虚部。转化之后再利用快速稀疏贝叶斯学习算法求得包含实部和虚部信息的目标矢量步骤c，假设噪声和目标图像矢量各像元均服从独立高斯先验建立稀疏先验模型，得到分块稀疏贝叶斯学习模型。可以推得目标矢量的后验概率也服从均值为μ，方差为∑的高斯分布，其中：

其中，表示对角元素由组成的对角矩阵，此时的均值μ即为对应的目标矢量估计值进一步的，所述步骤4中：

利用快速边缘似然函数最大化方法求解，实现分块回波数据目标图像重构，主要包括以下步骤：

步骤a，输入回波矢量和稀疏基对角块矩阵初始化迭代次数n＝1；

步骤b，将矩阵中的每一列基作为一个原子，将每个原子与回波矢量做内积，选择内积值最大的原子作为被选中的第一个原子fi，利用式和进行初始化，其中表示的标准差；

步骤c，计算∑,μ，引入两个中间变量gi和qi，表达式为：

其中，

步骤d，令gi＝αigi/(αi-gi)，qi＝αiqi/(αi-gi)，计算出每个原子对应的gi和qi

以及θi＝qi²-gi值，分三种情况进行处理：

①若θi＞0且αi＝+∞，则将fi添加到被选中的原子中；

②若θi＞0且αi＜+∞，则对超参数αi进行更新，有αi＝gi²/(qi²-gi)；

③若θi＜0，则将fi从被选中的原子中删除，并设置αi＝+∞；

步骤e，利用更新σ²，得到新的α和σ²后，继续更新∑,μ，并计算新的gi和qi；

步骤f，判断是否终止迭代，若满足迭代终止条件则停止迭代，此时的μ值即为目标矢量的估计值

进一步的，所述步骤5中：

将所求目标矢量估计值中的实部与相应的虚部相结合变为复数矢量信号，再将其转化为二维矩阵形式，得到分块回波数据对应的二维目标图像再将所求的每块回波数据所求的目标图像合成整个二维图像，即可得到重构的目标图像。参照图1，本发明的具体实施步骤如下：

步骤1，建立双基地isar成像回波模型假设目标平稳运动，并且其运动轨迹与双基地雷达基线共面。双基地isar平稳运动目标成像的几何模型如图2所示，其中，发射站雷达为t，接收站雷达为r，双基地雷达基线长度为l，等效单基地雷达位置为e。假设目标以速度v匀速运动，成像起始时刻为t1，双基地角为β1，以目标质心o1为原点，双基地角平分线方向为y1轴建立右手直角坐标系x1o1y1，散射点p1在x1o1y1坐标系中的坐标为(xp,yp)，设o1p1长度为d，与x1轴夹角为α1，目标质心o1、散射点p1到发射站t、接收站r的距离分为记为rt1、rr1和rtp1、rrp1。假设在tm成像时刻，目标质心运动至om点，散射点由p1运动至pm，此时的双基地角为βm，等效单基地雷达视角变化为θm(可看作由坐标系x1o1y1旋转到坐标系xmomym的旋转角度)。同样以目标质心om为原点，双基地角平分线方向为ym轴建立右手直角坐标系xmomym，ompm与xm轴夹角为αm且αm＝θm+α1，此时目标质心om、散射点pm到发射站t、接收站r的距离分别记为rtm、rrm和rtpm、rrpm。

双基地isar成像时，由于满足远场条件，则散射点pm到发射站和接收站的距离rtpm、rrpm可表示为：

假设雷达发射线性调频信号，经包络对齐和相位校正后的双基地isar回波的可近似为：

其中，fc为载波中心频率，tp为发射信号脉冲宽度，μ为调频斜率，σp为散射点p的信号复幅度，θ(tm)和β(tm)分别为成像期间内的旋转角度和双基地角，随慢时间tm变化。

为了避免双基地角时变引起越的分辨单元徙动和图像畸变，需要构造相应的补偿相位进行相位补偿，构造的补偿项表达式为：

得到相位补偿后的一维距离像表达式为：

假设在距离单元(2yp/c)cos(β(tm)/2)内有q个强散射点，则此单元的回波信号可表示为：

其中，aq为第q个散射点的信号复幅度。

步骤2，建立基于压缩感知的双基地isar稀疏孔径成像模型假设全孔径回波信号中共包含l个脉冲视角，累积转角为δθ，构造稀疏基矩阵fall将二维成像场景离散化为n个距离单元和m个多普勒单元，则在方位向上坐标可表示为xp＝mδx，其中，δx＝(l/m)(λ/2δθ)表示每个多普勒单元格在方位向的代表的尺度大小，回波信号可离散化为：

其中，l＝1,2,…l。可将稀疏基矩阵fall构造为：

考虑到实际噪声的存在，则双基地isar全孔径回波可稀疏表示为：

sall＝falla+ε0

其中，ε0为噪声，a为需求解的目标图像。假设有效的回波脉冲为j个，构造观测矩阵t，得到稀疏孔径回波数据：

s＝tsall+ε＝tfalla+ε＝fa+ε

步骤3，建立分块快速稀疏贝叶斯成像模型将整个二维回波数据进行分块处理，并假设目标图像各像元服从高斯先验建立稀疏贝叶斯模型，主要包括以下步骤：

步骤a，对二维回波数据分块化，即将回波矩阵sk×n按距离单元分成小块进行处理，如图3所示，假设每块数据矩阵中包含n1个距离单元的回波数据，则n个距离单元的回波数据一共可分为h＝[n/n1]块([]表示向上取整)，回波数据可表示为s＝[s1s2…sh]。为方便后续的算法推导，此处假设回波数据刚好可分为h个相同的数据块(若最后一块数据不足n1个距离单元，也可与前面的数据块按同样的方法处理)。在求解时可先利用重构算法进行逐块数据重构，然后再将得到的成像结果拼接成目标二维图像。

步骤b，只取其中一个数据块的回波矩阵进行分析求解，将分块回波数据按列矢量化，得到每块数据的矢量表达式：

式中sn、εn和an分别表示第“n”个距离单元对应的回波数据矢量、噪声矢量和目标图像矢量。

由于雷达回波信号为复数信号，将矢量化后的数据由复数转化为实数形式，得到表达式：

其中，re()和im()分别表示实部和虚部。转化之后再利用快速稀疏贝叶斯学习算法求得包含实部和虚部信息的目标矢量

步骤c，假设噪声和目标图像矢量各像元均服从独立高斯先验建立稀疏先验模型，得到分块稀疏贝叶斯学习模型。

首先，假设目标图像各像元的实部和虚部ai,(i＝1,2,…,2mn1)均满足依赖于超参数的零均值gaussian分布，即ai～n(0,αi^-1)，假设各像元之间是独立分布的，则其条件概率密度函数为：

其中，gamma(α|a,b)＝γ(a)^-1b^aα^a-1e^-bα，可以看出，目标图像矢量实际上是通过超参数α实现控制的，而超参数α又由系数a,b约束，所以对目标图像的稀疏促进作用可以看作是一个两层稀疏先验模型，如图4所示。利用贝叶斯准则可以推得目标矢量的后验概率也服从均值为μ，方差为∑的高斯分布，其中：

式中，表示对角元素由组成的对角矩阵，此时的均值μ即为对应的目标矢量估计值

步骤4，利用快速边缘似然函数最大化法实现分块数据图像重构利用快速边缘似然函数最大化方法求解，实现分块回波数据目标图像重构，实现流程示意图如图5所示，主要包括以下步骤：

步骤a，输入回波矢量和稀疏基对角块矩阵初始化迭代次数n＝1；

步骤b，将矩阵中的每一列基作为一个原子，将每个原子与回波矢量做内积，选择内积值最大的原子作为被选中的第一个原子fi，利用式和进行初始化，其中表示的标准差；

步骤c，计算∑,μ，引入两个中间变量gi和qi，表达式为：

其中，

步骤d，令gi＝αigi/(αi-gi)，qi＝αiqi/(αi-gi)，计算出每个原子对应的gi和qi

以及θi＝qi²-gi值，分三种情况进行处理：

①若θi＞0且αi＝+∞，则将fi添加到被选中的原子中；

②若θi＞0且αi＜+∞，则对超参数αi进行更新，有αi＝gi²/(qi²-gi)；

③若θi＜0，则将fi从被选中的原子中删除，并设置αi＝+∞；

步骤e，利用更新σ²，得到新的α和σ²后，继续更新∑,μ，并计算新的gi和qi；

步骤f，判断是否终止迭代，若满足迭代终止条件则停止迭代，此时的μ值即为目标矢量的估计值

步骤5，利用分块重构结果合成整个二维图像

将所求目标矢量估计值中的实部与相应的虚部相结合变为复数矢量信号，再将其转化为二维矩阵形式，得到分块回波数据对应的二维目标图像再将所求的每块回波数据所求的目标图像合成整个二维图像，即可得到重构的目标图像

本发明的效果可以通过下述仿真实验加以说明：

仿真实验1，算法运算时间仿真：

为验证算法的运算时间，选取简单的线性调频信号作为原始信号进行重构。假设信号为x(t)＝cos(2π(1000t+10t²/2))，采样频率为400hz，信号长度为l＝500，对原始信号进行随机稀疏采样，稀疏基矩阵为标准傅里叶基矩阵，利用本发明的快速边缘似然函数最大化算法和典型的直接求导法对稀疏贝叶斯学习模型进行求解，总迭代次数均设为1000次，收敛条件均设为在观测点数为m＝250时利用快速算法的重构结果如图6-a和图6-b所示，其中图6-a为时域重构结果，图6-b为频域重构结果，可以看出利用sbl算法能够很好地重构出目标信号。改变观测数m，采用运算时间和均方根误差rmse衡量两种求解方法对信号的重构性能，其中，rmse定义为结果如表1所示。

表1两种方法性能对比

仿真实验2，双基地isar稀疏孔径成像性能仿真：

1.设定仿真模型本仿真采用图7-a所示的双基地isar仿真场景，假设双基地基线长度为400km，目标在300km的高度以3km/s的速度沿基线方向由发射站向接收站匀速运动，运动轨迹为图中粗横线所示，成像起点运动轨迹高度上距接收站右侧水平距离70km处。目标的散射点仿真模型如图7-b所示，该模型由100个散射点组成，成像的仿真参数设置如表2所示。

表2成像参数设置

在本发明实例中假设成像观测时间为10s，累积脉冲数为500个，在此成像段内双基地角和累积转角变化曲线如图8-a和图8-b所示，其中，图8-a为双基地角变化曲线，图8-b为累积转角变化曲线。

2.仿真结果及分析

为比较本发明方法在稀疏孔径条件下的成像性能，在snr为10db条件下，假设三种数据缺失比分别为50％、65％和80％，为减少数据处理量，截取其中200个距离单元的数据进行处理，分别利用rd算法和本发明算法实现成像，结果如图9-a、图9-b和图9-c所示，其中，图9-a表示50％数据缺失时的一维距离像、rd成像结果以及本发明方法成像结果，图9-b表示65％数据缺失时的一维距离像、rd成像结果以及本发明方法成像结果，图9-c表示80％数据缺失时的一维距离像、rd成像结果以及本发明方法成像结果。

从成像结果可以看出，在稀疏孔径条件下，采用rd算法直接成像会有大量的能量泄漏并产生散焦现象，无法实现图像重构，而采用本发明算法成像能有效提高成像质量，能量泄漏问题和图像散焦问题得到明显改善，当孔径缺失数量较多时，如80％数据缺失的情况，仍能高质量地实现图像重构，得到清晰且聚焦良好的目标图像，体现了算法的优越性。

以上结合附图详细说明了本发明的技术方案，本发明中的步骤可根据实际需求进行顺序调整、合并和删减。

尽管参考附图详地公开了本发明，但应理解的是，这些描述仅仅是示例性的，并非用来限制本发明的应用。本发明的保护范围由附加权利要求限定，并可包括在不脱离本发明保护范围和精神的情况下针对发明所作的各种变型、改型及等效方案。