CN109870687B - 一种用于互质采样星载sar的成像处理方法 - Google Patents

Abstract

本发明公开了一种用于互质采样星载SAR的成像处理方法，属于信号处理领域。本发明方法在获得成像参数、回波数据及互质采样矩阵之后，首先对方位互质采样SAR回波信号进行距离向脉冲压缩，然后从第一个距离门开始，依据各距离门的多普勒参数计算该距离门回波信号跨越的距离门个数以及构造相应的稀疏字典，并以此为基准截取二维观测信号，最后用改进的二维观测信号稀疏度自适应的稀疏重构算法对二维观测信号进行重构，由此得到场景目标的后向散射信息。本发明补偿了成像参数随距离门变化对稀疏重构造成的影响，可实现全场景精确聚焦，能够实现互质采样工作模式下的成像处理，并且成像精度高、实用性强。

Description

一种用于互质采样星载SAR的成像处理方法

技术领域

本发明属于信号处理领域，涉及一种星载SAR成像处理方法，具体涉及一种用于互质采样星载SAR的成像处理方法。

背景技术

合成孔径雷达(Synthetic Aperture Radar，SAR)是一种主动式微波成像传感器，其通过发射宽带信号，结合合成孔径技术，能够在距离向和方位向上同时获得高分辨率图像，实现对目标区域的远距离主动成像。与传统光学遥感和高光谱遥感相比，SAR具备全天候、全天时、高分辨对地观测的能力，还具有一定的穿透性，是实现军事侦察、自然资源普查、自然灾害监测等的重要技术手段。

高分辨率和宽测绘带对地观测是星载SAR的发展趋势之一。对于传统的星载SAR系统，方位空间分辨率与距离向测绘幅宽对SAR系统的脉冲重复频率(Pluse RepeatedFrequency，PRF)提出了相互矛盾的约束要求：一方面，方位高分辨率要求SAR系统的PRF足够高，以保证方位采样满足奈奎斯特采样定理；另一方面，宽测绘带要求SAR系统的PRF足够低，以保证测绘带回波在两次脉冲发射间隔时间内被完整接收。受固定PRF的限制，传统星载SAR不能同时实现高分辨率和宽测绘带对地观测。

互质阵列和互质采样是近些年在阵列信号领域中被提出的一种新颖的稀疏非均匀阵列及采样方式。互质阵列由两个均匀线性阵列(Uniform Linear Array,ULA)组成，如图1所示，第一个ULA包含M个阵元，阵元间隔为Nλ/2，第二个ULA包含N个阵元，阵元间隔为Mλ/2，其中M和N为互质数，λ为工作波长。由于两个子阵共用第一个阵元，一个(M,N)-互质阵列仅有M+N-1个阵元，但其可提供的自由度达到MN，远大于其物理阵元数。此外，互质阵列还有阵元间互耦影响小及硬件上易于实现的优点，因此在空间谱估计，也称波达方向(Direction-Of-Arrival，DOA)估计，波束形成等领域得到广泛的关注和应用。

互质采样星载SAR将互质采样技术引入星载SAR系统，其利用互质阵列可由低物理阵元数提供高自由度的特点，在采样率远低于奈奎斯特采样定律的条件下实现具有稀疏特性场景的精确重构，有效减小雷达系统所需要的方位采样数目。此时，若进一步采用时分复用技术，如图2所示，在非互质采样时刻将波束切换至其他测绘带区域，可有效增加SAR系统的测绘带宽度。

由于在方位向采用了互质稀疏采样，如何实现互质采样星载SAR回波信号成像处理变得十分困难。目前对于互质采样SAR成像方法的研究中，均假设在SAR回波的距离方位耦合和二维空变特性可忽略不计的情况下，引入压缩感知(CS)理论，构造由方位调频信号组成的稀疏字典，采用压缩感知理论中的稀疏重构算法对各方位信号进行重构。参考文献1[Tao Y,Zhang G,Li D.Coprime sampling with deterministic digital filters incompressive sensing radar[C]//Cie International Conference on Radar.IEEE,2017:1-4.]提出一种采用确定性滤波器实现的互质采样方式用于CS雷达采样。参考文献2[史洪印,贾宝京,齐兆龙.基于压缩感知的非均匀脉冲SAR欺骗性干扰抑制方法[J].仪器仪表学报,2016,37(3):525-532.]将SAR信号方位向互质采样方式与嵌套稀疏采样方式同CS算法相结合，不但可以抑制欺骗式干扰，而且能够在远低于Nyquist采样率的条件下实现目标的高分辨成像。参考文献3[Han B,Wang P,Fang Y,et al.The Research of SARImaging Based on Co-Prime Arrays Sampling[C]//Advances in Materials,Machinery,Electrical Engineering.2017.]在距离向和方位向均进行互质采样，并采用CS算法实现SAR信号的二维精确聚焦。上述几篇文献是从不同角度探讨互质采样在SAR中的应用，有不同的侧重点，但对于互质采样SAR成像方面，他们主要侧重分析用压缩感知成像的可行性，因此都是在忽略SAR回波信号二维耦合和空变特性情况下，直接对SAR回波的一维方位向信号进行压缩感知重构，采用的重构算法也基本是压缩感知中最普遍的OMP算法，并且采用这种重构算法需要已知信号稀疏度，不能实现稀疏度自适应。

所以针对目前由于完全忽略了SAR回波信号的二维耦合及二维空变特性，无法实现回波信号距离方位二维耦合补偿处理及成像参数空变补偿处理的情况，本发明实现了一种用于互质采样星载SAR的成像处理方法，主要侧重于互质采样SAR成像算法的精确性及实用性。

发明内容

本发明的目的是为了解决互质采样模式下星载SAR的成像问题，提出了一种用于互质采样星载SAR的成像处理方法，能够完成距离方位二维耦合不可忽略情况下的精确聚焦处理。

本发明的用于互质采样星载SAR的成像处理方法，实现步骤包括：

步骤一：读入成像参数，包括互质数co₁、co₂以及回波数据方位向采样点数M、回波数据距离向采样点数N_r；根据回波数据方位向采样点数M和回波数据距离向采样点数N_r读入互质采样星载SAR的回波数据S_cos。

步骤二：根据互质数co₁、co₂以及回波数据方位向采样点数M构造互质采样矩阵Φ。

步骤三：对回波数据S_cos进行距离向脉冲压缩；设距离向脉冲压缩后的回波数据表示为S_{cos_rc}(τ,η)，τ表示距离向时间，η表示方位向时间。

针对回波数据S_{cos_rc}(τ,η)，按照距离门次序，从它的第一个距离门开始，依次对每个距离门执行步骤四～步骤六。

步骤四：二维观测信号截取；

利用雷达的等效斜视距离模型，获取距离门k对应的雷达波束中心照射到点目标的斜距R_k，并得到相应的波束照射范围内的最短斜距与斜距R_k之差跨越的距离门数Δk₁、波束照射范围内的最大斜距与斜距R_k之差跨越的距离门数Δk₂；

判断k＞Δk₁及k+Δk₂＜N_r是否成立，若是，在S_{cos_rc}(τ,η)上截取第k-Δk₁～k+Δk₂个距离门作为二维观测数据y；否则，停止对当前距离门k的操作，k自增1，继续对距离门k执行步骤四。

步骤五：对步骤四截取的二维观测数据y构造稀疏字典。

首先，计算距离向脉冲压缩包络p_r(τ_l)，τ_l,l＝1,…,N_r为距离向采样时刻序列；

然后，构造二维参考信号s_{k_rc}(τ_l,η_i)如下：

s_{k_rc}(τ_l,η_i)＝p_r[τ_l-2R(η_i)/c]exp{-j4πf₀R(η_i)/c},τ_l∈[τ_a,τ_b]

其中，τ_a＝2min[R(R(η_i)＞0)]/c，τ_b＝2max[R(η_i)]/c；c为光速；R(R(η_i)＞0)表示R(η_i)中所有非零值；min(·)表示求最小值函数，max(·)表示求最大值函数；η_i,i＝1,…,N_a为方位时间序列，R(η_i)是在距离门k对应的波束照射范围内的斜距；f₀为雷达工作频率；

最后，将二维参考信号s_{k_rc}(τ_l,η_i)在方位向进行循环移位，舍去由于循环移位导致有效回波位置不连续的原子，得到对应于观测数据y的稀疏字典

步骤六：利用互质采样矩阵Φ和稀疏字典Ψ对观测数据y进行稀疏重构。在稀疏重构时，首先获得传感矩阵A＝ΦΨ，然后将A中所有二维原子均拉成一维列向量，将二维观测信号y也拉成一个一维列向量，然后对观测数据y进行稀疏重构。

所述的步骤六中，对观测数据y进行稀疏重构的过程包括：

步骤601、由下式获得大小为M×(Δk×N)传感矩阵A＝ΦΨ；

步骤602、将传感矩阵A中所有二维原子均拉成一维列向量，此时A的大小为(M×Δk)×N，用a_j表示拉长后的各原子，j＝1,2,…,N表示各原子在A中的序列号。将二维观测信号y也拉成一个一维列向量，此时y的大小为(M×Δk)×1；

步骤603、进行初始化，包括：设迭代次数为t，初始t＝1，t≤M；设重构残差为r_t，设置初始残差r₀＝y；设第t次迭代后的传感矩阵A中原子的序号集合为Λ_t，初始的序列集合Λ₀为空集；设L_y表示选择与二维观测信号y最匹配的原子步长，L_r表示选择与残差r_t最匹配的原子步长，初始设置L_y＝l，L_r＝l，l表示初始步长；设

表示二维观测信号y的稀疏表示系数估计，其大小为N×1，初始

步骤604、预测试一，获得原子初选集P，分三步进行：

第一步，筛选A中有效回波位置完全包含于二维观测数据y中有效回波位置的原子，这些原子所在A中的序号构成原子预备集R；第二步，设A_R＝{a_j}(forallj∈R)，对A_R中各原子与y进行相关性测试，计算相关度u_{R_y}＝abs[A_R ^Ty]，其中，abs[·]表示求模运算；第三步，选择u_{R_y}中所有大于阈值Th＝0.1×max(u_{R_y})的值所对应的A中原子的序号构成原子初选集P，设A_P＝{a_j}(forallj∈P)；

步骤605、预测试二，获得原子候选集C_t，分三步进行：

第一步，选择原子初选集P中前L_y个与y之间的相关度最大值对应的A中原子的序号构成集合S_yt；第二步，对A_P中各原子与残差r_t-1进行相关性测试，相关度u_{P_r}＝abs[A_P ^Tr_t-1]，选择所有u_{P_r}中前L_r个最大值对应的A中原子的序号构成集合S_rt；第三步，令原子候选集C_t＝Λ_t-1∪S_yt∪S_rt，其中，∪表示集合并运算；设

步骤606、最终测试，获得原子支撑集F_t，包括：对

中各原子按使残差最小原则测试，也就是求

的最小二乘解，θ_t表示在第t次迭代的y的稀疏表示系数，得到θ_t的估计值

选择θ_t中绝对值最大的前L_y项记为

对应

中的L_y个原子记为A_tL，这L_y个原子对应在A中的序号记为Λ_tL，令原子支撑集F_t＝Λ_tL；

步骤607、计算残差r_new＝y-A_tL(A_tL ^TA_tL)^-1A_tL ^Ty，并计算该残差能量对于初始残差能量的减少率rate_0＝(||r₀||₂-||r_new||₂)/||r₀||₂，以及相邻两次迭代残差能量的减少量对初始残差能量的比率rate_1＝(||r_t-1||₂-||r_new||₂)/||r₀||₂，其中||·||₂表示求向量的二范数运算；

步骤608、进行迭代条件判断，包括：

首先对rate_0进行判断，如果rate_0＞1-ε0，则令最终原子支撑集F_final＝F_t，对应最终稀疏表示系数

并停止迭代进入步骤609；否则，继续对rate_1进行判断；

如果rate_1＜0，则更新步长L_y＝L_y+l，L_r＝L_r+l，返回步骤605继续迭代；否则，更新步长L_y＝L_y+l，令L_r＝l，Λ_t＝F_t，r_t＝r_new，t自增1，并进行以下判断：若rate_1＞ε1，则令最终原子支撑集F_final＝F_t，对应最终稀疏表示系数

并返回步骤605继续迭代；若rate_1≤ε1，则停止迭代进入步骤609；ε0和ε1是设置的参数，取值接近于0；

步骤609、根据最终原子支撑集F_final和对应的最终稀疏表示系数

得到大小为N×1的稀疏表示向量

最后，在向量

两端分别补V/2个0，得到y最终的稀疏表示系数估计向量

中的非零元素及位置表示二维观测数据y所包含地面目标的后向散射信息，且为最终压缩成像数据矩阵

的第k列；

在对回波数据S_{cos_rc}(τ,η)的所有距离门完成从步骤四到步骤六的操作后，把每个距离门对应的稀疏表示系数估计向量按距离门序号进行排列，得到回波数据S_cos的最终压缩成像数据矩阵

本发明与现有技术相比，具有以下优势：

(1)精度高。本发明根据距离单元徙动数截取二维观测信号，并通过所改进的稀疏度自适应的二维信号稀疏重构算法，能够完成距离方位二维耦合不可忽略情况下的精确聚焦处理。

(2)全场景聚焦。本发明根据各距离门的多普勒参数计算截取观测信号的距离门个数以及构造稀疏字典，补偿了成像参数随距离门变化对稀疏重构造成的影响，由此可实现全场景精确聚焦。

(3)实用性强。本发明能够实现互质采样工作模式下的成像处理，相比于目前忽略了距离方位耦合补偿和成像参数空变补偿的成像处理方法，本发明具有更好的实用性。

附图说明

图1是互质阵列结构示意图；

图2是互质采样星载SAR工作方式示意图；

图3是本发明互质采样星载SAR成像处理方法的流程图；

图4是雷达的等效斜视距离模型的示意图；

图5是本发明中二维信号稀疏度自适应的稀疏重构算法的流程图；

图6是本发明多点目标成像仿真处理结果示意图，(a)为多点目标成像仿真结果；(b)为(a)中右下角点目标的水平剖面示意图；(c)为在距离向脉冲压缩时加入Kaiser窗的多点目标成像仿真结果；(d)为(c)中右下角点目标的水平剖面示意图。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图对本发明作进一步的详细和深入描述。

本发明的一种用于互质采样星载SAR的成像处理方法，在获得成像参数、回波数据及互质采样矩阵之后，首先对方位互质采样SAR回波信号进行距离向脉冲压缩，然后从第一个距离门开始，依据各距离门的多普勒参数计算该距离门回波信号跨越的距离门个数以及构造相应的稀疏字典，并以此为基准截取二维观测数据，最后用改进的二维观测信号稀疏度自适应的稀疏重构算法对二维观测信号进行重构，由此得到场景目标的后向散射信息。

本发明方法整体上主要包括计算互质采样矩阵、距离向脉冲压缩、二维观测信号截取、稀疏字典构造和二维观测信号稀疏重构五个部分。本发明实施例中，设置9个单位幅度点目标，显然单个点目标在方位向的稀疏度为1。根据所设置的多点目标进行互质采样星载SAR回波仿真，得到回波数据，进行成像处理，本发明方法流程如图3所示，下面说明各步骤。

步骤一：读入相关的成像参数和需要进行成像处理的互质采样星载SAR的回波数据，具体包括步骤101和102。

步骤101、读入成像参数，包括：雷达工作频率f₀、发射脉冲时宽T_r、发射脉冲调频率K_r、距离采样率f_s、脉冲重复频率PRF、互质数co₁和co₂、方位向波束宽度θ_bw、参考斜距R_ref、雷达波束照射中心时刻η_c、回波数据方位向采样点数M、回波数据距离向采样点数N_r、各个距离门k的多普勒中心频率f_{D_k}、各个距离门k的多普勒调频率f_{R_k}，k为距离门编号，k∈{1,2,…,N_r}。

步骤102、根据回波数据方位向采样点数M和回波数据距离向采样点数N_r，读入需要进行成像处理的互质采样星载SAR的回波数据S_cos，其大小为M×N_r。

其中，本实施例中具体成像参数为：f₀＝10GHz，T_r＝40μs，K_r＝1.5MHz/μs，f_s＝72MHz，PRF＝5680Hz，co₁＝4，co₂＝5，θ_bw＝0.0089，R_ref＝850km，η_c＝-0.8358s，M＝2913，N_r＝4096，f_{D_k}和f_{R_k}与距离门k有关，如k＝2048时，f_{D_k}＝3306.8Hz，f_{R_k}＝-3956.3Hz/s。

步骤二：根据互质数co₁、co₂以及回波数据方位向采样点数M构造互质采样矩阵，具体包括如下子步骤201～203。

步骤201、计算原始回波数据S_o的方位向长度N_a，如下：

rate＝(co₁+co₂-1)/co₁/co₂ (1)

其中，rate表示互质采样的降采样率，

表示向上取整。原始回波数据S_o表示未进行方位互质采样时的回波数据，S_o是未知的，大小为N_a×N_r，在实际中直接获得互质采样后的回波数据S_cos，此处目的是为了计算S_o的方位向长度N_a。

在本实施例中，计算得到rate＝0.4，N_a＝7283。

步骤202、根据互质数co₁、co₂计算互质采样时刻序列H，如下：

其中，H中元素表示互质采样时刻，H中元素按升序排列且无重复元素，

表示向下取整。

步骤203、构造互质采样矩阵Φ如下：

其中，m＝1,2,…,M，H_m表示互质采样时刻序列H中第m个元素。向量

是一个分段序列，表示为{X_n}，当n＝P_m+1时，X_n为1，当n＝1,2,…,N_a,且n≠H_m+1时，X_n为0。

由步骤二构造的互质采样矩阵Φ大小为M×N_a，满足S_cos＝ΦS_o。

步骤三：对互质采样星载SAR的回波数据S_cos进行距离向脉冲压缩。本发明采用传统匹配滤波的脉冲压缩方式对互质采样星载SAR回波信号进行距离向聚焦，包括如下子步骤301～303。

步骤301、互质采样回波数据S_cos在距离向时间域的表示为S_cos(τ,η)，τ表示距离向时间，η表示方位向时间，对在距离向时域的回波数据S_cos(τ,η)进行距离向快速傅里叶变换得到距离向频域的回波数据S_cos(f_r,η)如下：

S_cos(f_r,η)＝FFT_τ{S_cos(τ,η)} (6)

其中，f_r表示距离向频率，FFT_τ{·}表示进行距离向快速傅立叶变换。

步骤302、在距离向频域构造如(7)式所示的滤波器H(f_r)：

H(f_r)＝exp{jπf_r ²/K_r} (7)

其中，j表示虚数单位。

步骤303、将滤波器H(f_r)与变换到距离向频域的回波数据S_cos(f_r,η)相乘，再利用快速傅里叶逆变换将相乘的结果变换到距离向时域，得到已完成距离向聚焦的信号S_{cos_rc}(τ,η)：

S_{cos_rc}(τ,η)＝IFFT_τ{S_cos(f_r,η)H(f_r)} (8)

其中，IFFT_τ{·}表示进行距离向快速傅里叶逆变换。

下面将从已完成距离向聚焦的信号S_{cos_rc}的第一个距离门开始，按照距离门次序，逐次对每个距离门完成步骤四～步骤六。

步骤四：截取二维观测信号。

如图4所示，X轴表示雷达平台飞行方向，Y表示垂直于地面方向，Z表示地面垂直于雷达平台飞行方向，波束中心时刻η_c表示雷达波束中心正好照射到目标Tar的时刻。对于第k个距离门而言，波束中心时刻雷达平台到地面目标间的斜距为R_k，对应的等效速度和等效斜视角分别为V_k和

R(η_i)表示在η_i时刻雷达平台到地面目标间的斜距，η_i,i＝1,2,…N_a代表方位时间序列。设Δk₁表示波束照射范围内的最短斜距与斜距R_k之差跨越的距离门数，Δk₂就表示波束照射范围内的最大斜距与斜距R_k之差跨越的距离门数，设Δk＝Δk₁+Δk₂。如果k＞Δk₁，且k+Δk₂＜N_r，则在已完成距离向聚焦的信号S_{cos_rc}上截取第k-Δk₁～k+Δk₂个距离门信号作为二维观测信号y，y大小为M×Δk；否则，距离门k＝k+1，对下一个距离门执行步骤四～六。

其中，计算Δk₁和Δk₂的步骤如下：

步骤401、计算距离门k对应的斜距R_k：

R_k＝R_ref+(k-N_r/2)bin_r (9)

bin_r＝c/2/f_s (10)

其中，bin_r表示距离采样间隔对应的斜距变化量，c为光速，取3×10⁸m/s。

步骤402、根据距离门k对应的多普勒中心频率f_{D_k}和多普勒调频率f_{R_k}计算等效速度V_k和等效斜视角

其中，λ表示波长，V_k表示距离门k对应的等效速度，

表示距离门k对应的等效斜视角。

步骤403、计算方位时间序列η_i公式如下：

η_i＝(i-N_a/2)/PRF+η_c,i＝1,…,N_a (13)

在不同的方位时间η_i，对应一个等效斜距R(η_i)。

步骤404、计算方位向包络w_a(η_i)公式如下：

其中，abs(·)表示求绝对值函数，

表示雷达在运动过程中在方位时间η_i时、波束中心偏离目标的方位离轴角在斜距平面上的投影。

步骤405、计算在波束照射范围内的斜距R(η_i)公式如下：

步骤406、计算Δk₁和Δk₂，公式如下：

Δk₁＝{R_k-min[R(R(η_i)＞0)]}/bin_r (17)

Δk₂＝{max[R(η_i)]-R_k}/bin_r (18)

其中，R(R(η_i)＞0)表示R(η_i)中所有非零值，min(·)表示求最小值函数，max(·)表示求最大值函数。

本发明截取的观测信号为二维信号，可将SAR回波的二维耦合作为一个整体进行重构，避免了在方位稀疏非均匀采样情况下难以实现的二维解耦操作；而且截取观测信号的距离门数根据回波各距离门的多普勒参数进行更新，使得回波时频分布特性随距离门变化对重构算法产生的影响得以消除。

步骤五：根据步骤四中截取的距离门k的二维观测信号y，构造稀疏字典，具体如下三个子步骤501～503。本发明方法中通过等效斜距模型构造参考信号，使得稀疏字典的精确度提高，从而提高重构精度；在下面构造稀疏字典时舍去了有效回波位置不连续的原子，使得在重构过程中可以避免非有效原子的干扰，同时还可以避免额外的计算开销。

步骤501、计算距离向采样时刻序列τ_l公式如下：

τ_l＝(l-N_r/2)/f_s+2R_ref/c,l＝1,…,N_r (19)

其中，l表示采样步长。

计算距离向时域包络w_r(τ_l)公式如下：

计算距离向频域包络W_r(f_τ)公式如下：

W_r(f_τ)＝w_r(f_τ/K_r) (21)

其中，f_τ表示距离向频率。

计算距离向脉冲压缩包络p_r(τ_l)公式如下：

其中，sinc是表示辛格函数。

步骤502、构造二维参考信号s_{k_rc}(τ_l,η_i)如下：

s_{k_rc}(τ_l,η_i)＝p_r[τ_l-2R(η_i)/c]exp{-j4πf₀R(η_i)/c},τ_l∈[τ_a,τ_b] (23)

其中，τ_a＝2min[R(R(η_i)＞0)]/c，τ_b＝2max[R(η_i)]/c；R(η_i)为波束照射范围内的斜距，由(16)式计算。

步骤503、将二维参考信号s_{k_rc}(τ_l,η_i)在方位向进行循环移位，得到对应于二维观测信号y的稀疏字典。若单点目标原始回波中有效回波的方位向采样点数为V，V是互质采样SAR回波数据的原始回波数据(方位长度为N_a)中有效回波部分的采样点数，V＜N_a，舍去由于循环移位导致有效回波位置不连续的原子，得到稀疏字典Ψ，如公式(24)所示，其大小为N_a×(Δk×N)，其中包含的原子个数N＝N_a-V，各原子大小为N_a×Δk。

步骤六：根据改进的稀疏度自适应匹配追踪算法对二维观测信号进行稀疏重构。

本发明将稀疏度自适应匹配追踪算法(Sparsity Adaptive Matching Pursuit，SAMP)(参考文献[4]：Do T T,Lu G,Nguyen N,et al.Sparsity adaptive matchingpursuit algorithm for practical compressed sensing[C]//Signals,Systems andComputers,2008,Asilomar Conference on.IEEE,2009:581-587.)进行改进，使其适用于步骤四中所得二维观测信号y的稀疏重构。

本步骤对y进行稀疏重构的流程如图5所示，包括如下子步骤601～609。

步骤601、由公式(25)获得大小为M×(Δk×N)的传感矩阵A：

A＝ΦΨ (25)

步骤602、将传感矩阵A中所有二维原子均拉成一维列向量，此时A的大小为(M×Δk)×N，用a_j表示拉长后的各原子，j＝1,2,…,N表示各原子在A中的序列号。将二维观测信号y也拉成一个一维列向量，此时y的大小为(M×Δk)×1。

本发明方法首先将二维观测信号y、传感矩阵A的所有二维原子均拉长成一维列向量，通过这样简单的拉长操作，使一维信号的稀疏重构算法适用于二维信号的稀疏重构。

步骤603、初始化：设重构残差为r_t，设置初始残差r₀＝y；设迭代次数为t，初始t＝1，t≤M；设置第t次迭代后的传感矩阵A中原子的序号集合为Λ_t，第一次迭代前初始的序列集合

表示空集；设L_y表示选择与二维观测信号y最匹配的原子步长，L_r表示选择与残差r_t最匹配的原子步长，初始设置L_y＝l，L_r＝l，l表示初始步长，本实施例中，l＝1；设

表示二维观测信号y的稀疏表示系数估计，其大小为N×1，初始

步骤604、预测试一，获得原子初选集P：

分三小步进行，第一步，筛选A中有效回波位置完全包含于二维观测信号y中有效回波位置的原子，这些原子所在A中的序号构成原子预备集R；第二步，设A_R＝{a_j}(forallj∈R)，对A_R中各原子与二维观测信号y进行相关性测试，即计算u_{R_y}＝abs[A_R ^Ty]，其中，abs[·]表示求模运算，u_{R_y}表示A_R中各原子与观测信号y之间的相关度；第三步，选择u_{R_y}中所有大于阈值Th＝0.1×max(u_{R_y})的值，这些值对应A中原子的序号构成原子初选集P，这些值就是原子初选集P中元素与y之间的相关度u_{P_y}，设A_P＝{a_j}(forallj∈P)。后续的迭代过程将在A_P中进行。

通过原子初选集的筛选，可以排除有效回波位置不匹配的原子，大大减小了重构算法不必要的计算开销，提高重构效率。

步骤605、预测试二，获得原子候选集C_t：

也分三小步进行，第一步，选择原子初选集P中前L_y个与y的相关度最大值u_{P_y}对应A中原子的序号构成集合S_yt；第二步，对A_P中各原子与残差r_t-1进行相关性测试，即计算u_{P_r}＝abs[A_P ^Tr_t-1]，选择所有u_{P_r}中前L_r个最大值对应的A中原子的序号构成集合S_rt；第三步，令原子候选集C_t＝Λ_t-1∪S_yt∪S_rt，其中，符号∪表示集合并运算。

原子候选集的选取中，通过分别测试原子初选集与观测信号的相关度以及原子初选集与残差的相关度进行筛选，且选择与观测信号最匹配原子的步长和选择支撑集的步长一致，并随迭代次数递增，但选择与残差最匹配原子的步长只随迭代阶段递增。

步骤606、最终测试，获得原子支撑集F_t：

对

中各原子按使残差最小原则测试，即求

的最小二乘解：

选择θ_t中绝对值最大的前L_y项记为

对应

中的L_y个原子记为A_tL，这L_y个原子对应在A中的序号记为Λ_tL，并令原子支撑集F_t＝Λ_tL；θ_t表示二维观测信号y的在第t次迭代的稀疏表示系数，θ_t表示在第t次迭代的稀疏表示系数估计值。

步骤607、计算残差r_new＝y-A_tL(A_tL ^TA_tL)^-1A_tL ^Ty，并计算此残差能量对于初始残差能量的减少率rate_0＝(||r₀||₂-||r_new||₂)/||r₀||₂，以及相邻两次迭代残差能量的减少量对初始残差能量的比率rate_1＝(||r_t-1||₂-||r_new||₂)/||r₀||₂，其中||·||₂表示求向量的二范数运算。

步骤608、判断迭代条件：

首先对rate_0进行判断：

如果rate_0＞1-ε0，则令最终原子支撑集F_final＝F_t，对应最终稀疏表示系数

并停止迭代进入步骤609；

否则，则对rate_1进行判断：

如果rate_1＜0，则更新步长L_y＝L_y+l，L_r＝L_r+l，返回步骤605继续迭代；

否则，则更新步长L_y＝L_y+l，令L_r＝l，Λ_t＝F_t，r_t＝r_new，t＝t+1。并且进行以下判断：若rate_1＞ε1，则令最终原子支撑集F_final＝F_t，对应最终稀疏表示系数

并返回步骤605继续迭代；若rate_1≤ε1，则停止迭代进入步骤609。

在本发明实施例中，取ε0＝1e-4，ε1＝1e-8。

需要说明，ε0、ε1均取接近于0的值。其中，ε0为判断点目标重构的迭代停止条件参数。点目标重构成功时残差能量接近于0，但对面目标来说，所截取的观测信号为目标所有散射中心回波在所截取距离门范围内的叠加结果，在迭代过程中残差能量不会减小至0，在本发明中一个稳妥的方法是取ε0＝1e-4。ε1为判断面目标重构的迭代停止条件参数，即相邻两次迭代的残差减小量对初始残差的比率足够小时停止迭代。ε1取值越接近于0，重构的精度就越高，但不会超过本发明中稀疏重构算法的理论重构精度，即方位向采样间隔。而另一方面，如果ε1过小，将导致不必要的算法计算开销，综合考虑以及进行实验，取ε1＝1e-8最为合适。

由于本发明中二维观测信号的截取以及稀疏字典的构造方式，使得面目标的二维观测信号在迭代过程中残差能量不能减小至零或者无限接近噪声功率， 因此本发明在判断残差能量减小率的基础上，引入对连续两次迭代的残差能量减小量对初始残差能量的比率的判断条件，即使连续两次迭代之间的残差能量减小量在小于某个阈值时停止迭代。

步骤609、根据最终原子支撑集F_final和对应的最终稀疏表示系数

得到大小为N×1的稀疏表示向量

最后，在向量

两端分别补V/2个0，得到二维观测信号y最终的稀疏表示系数估计向量

稀疏表示系数估计向量

中的非零元素及位置即表示二维观测信号y所包含地面目标的后向散射信息，且为最终压缩成像数据矩阵

的第k列。当对S_{cos_rc}所有距离门都完成从步骤四到步骤六的所有操作后，把每个距离门对应的稀疏表示系数估计向量按距离门序号进行排列，就可以得到回波数据S_cos的最终压缩成像数据矩阵

本实施例的成像仿真结果如图6所示，其中图6(a)为本实施例中的多点目标成像仿真结果，可以看出仅利用原始回波40％的数据量就完成了方位向的精确聚焦，且没有方位向旁瓣。图6(b)为(a)中右下角点目标成像结果在距离向插值后的旁瓣情况，发现距离向峰值旁瓣比PSLR＜-13dB，这说明，本发明中的稀疏重构方法没有影响传统SAR距离多普勒成像方法在距离向的分辨率。此外，图6(c)显示了在距离向脉冲压缩时加入Kaiser窗以抑制距离旁瓣的多点目标成像仿真结果，图6(d)同样为(c)中右下角点目标成像结果在距离向插值后的旁瓣情况，可发现在加入Kaiser窗以后距离向峰值旁瓣比PSLR＜-60dB。

Claims (9)

1.一种用于互质采样星载SAR的成像处理方法，其特征在于，包括如下步骤：

步骤一：读入成像参数，包括互质数co₁、co₂以及回波数据方位向采样点数M、回波数据距离向采样点数N_r；根据回波数据方位向采样点数M和回波数据距离向采样点数N_r读入互质采样星载SAR的回波数据S_cos；

步骤二：根据互质数co₁、co₂以及回波数据方位向采样点数M构造互质采样矩阵Φ，并计算原始回波数据的方位向长度N_a；

步骤三：对回波数据S_cos进行距离向脉冲压缩，设距离向脉冲压缩后的回波数据表示为S_{cos_rc}(τ,η)，τ表示距离向时间，η表示方位向时间；

针对回波数据S_{cos_rc}(τ,η)，按照距离门次序，从它的第一个距离门开始，依次对每个距离门执行步骤四～步骤六；设距离门编号为k,k∈{1,2,…,N_r}；

步骤四：截取二维观测数据；

利用雷达的等效斜视距离模型，获取距离门k对应的雷达波束中心照射到点目标的斜距R_k，并得到相应的波束照射范围内的最短斜距与斜距R_k之差跨越的距离门数Δk₁、波束照射范围内的最大斜距与斜距R_k之差跨越的距离门数Δk₂；

判断k＞Δk₁及k+Δk₂＜N_r是否成立，若是，在S_{cos_rc}(τ,η)上截取第k-Δk₁～k+Δk₂个距离门作为二维观测数据y；否则，停止对当前距离门k的操作，k自增1，继续对下一个距离门执行步骤四；

步骤五：对步骤四截取的二维观测数据y构造稀疏字典；

首先，计算距离向脉冲压缩包络p_r(τ_l)；其中τ_l,l＝1,…,N_r为距离向采样时刻序列；

然后，构造二维参考信号s_{k_rc}(τ_l,η_i)如下：

s_{k_rc}(τ_l,η_i)＝p_r[τ_l-2R(η_i)/c]exp{-j4πf₀R(η_i)/c},τ_l∈[τ_a,τ_b]

其中，τ_a＝2min[R(R(η_i)＞0)]/c，τ_b＝2max[R(η_i)]/c，c为光速，R(R(η_i)＞0)表示R(η_i)中所有非零值，η_i,i＝1,…,N_a为方位时间序列，R(η_i)是在距离门k对应的波束照射范围内的斜距；min(·)表示求最小值函数，max(·)表示求最大值函数；f₀为雷达工作频率；

最后，将二维参考信号s_{k_rc}(τ_l,η_i)在方位向进行循环移位，舍去由于循环移位导致有效回波位置不连续的原子，得到对应于观测数据y的稀疏字典Ψ；

步骤六：利用互质采样矩阵Φ和稀疏字典Ψ对观测数据y进行稀疏重构；

在稀疏重构时，首先获得传感矩阵A＝ΦΨ，然后将A中所有二维原子均拉成一维列向量，将二维观测信号y也拉成一个一维列向量，然后对观测数据y进行稀疏重构。

2.根据权利要求1所述的方法，其特征在于，所述的步骤一中，读入的成像参数还包括：雷达工作频率f₀、发射脉冲时宽T_r、发射脉冲调频率K_r、距离采样率f_s、脉冲重复频率PRF、方位向波束宽度θ_bw、参考斜距R_ref、雷达波束照射中心时刻η_c、各个距离门k的多普勒中心频率f_{D_k}、各个距离门k的多普勒调频率f_{R_k}。

3.根据权利要求1或2所述的方法，其特征在于，所述的步骤二中构造互质采样矩阵，包括：

步骤201、计算原始回波数据S_o的方位向长度N_a，如下：

rate＝(co₁+co₂-1)/co₁/co₂

其中，原始回波数据S_o表示未进行方位互质采样前的回波数据，大小为N_a×N_r，rate表示互质采样的降采样率，

表示向上取整；

步骤202、计算互质采样时刻序列H，如下：

其中，H中元素表示互质采样时刻，H中元素按升序排列且无重复元素，

表示向下取整；

步骤203、构造互质采样矩阵Φ，如下：

且n≠H_m+1,m＝1,2,…,M

其中，H_m表示互质采样时刻序列H中第m个元素，向量

是一个分段序列，表示为{X_n}。

4.根据权利要求1所述的方法，其特征在于，所述的步骤三采用匹配滤波的脉冲压缩方式对回波数据S_cos进行距离向压缩，包括：

步骤301、对距离向时间域的回波数据S_cos(τ,η)进行距离向快速傅里叶变换得到距离向频域的回波数据S_cos(f_r,η)＝FFT_τ{S_cos(τ,η)}；其中，f_r表示距离向频率，FFT_τ{·}为进行距离向快速傅立叶变换；

步骤302、在距离向频域构造滤波器H(f_r)＝exp{jπf_r ²/K_r}；其中，j表示虚数单位，K_r为发射脉冲调频率；

步骤303、将滤波器H(f_r)与S_cos(f_r,η)相乘，再利用快速傅里叶逆变换将相乘的结果变换到距离向时域，完成距离向脉冲压缩；

距离向脉冲压缩后的回波数据S_{cos_rc}(τ,η)＝IFFT_τ{S_cos(f_r,η)H(f_r)}；其中，IFFT_τ{·}表示进行距离向快速傅里叶逆变换。

5.根据权利要求1所述的方法，其特征在于，所述的步骤四中，计算Δk₁和Δk₂的步骤包括：

步骤401、计算距离门k对应的斜距R_k如下：

R_k＝R_ref+(k-N_r/2)bin_r

bin_r＝c/2/f_s

其中，R_ref为参考斜距，bin_r表示距离采样间隔对应的斜距变化量，f_s为距离采样率；

步骤402、根据距离门k对应的多普勒中心频率f_{D_k}和多普勒调频率f_{R_k}计算等效速度V_k和等效斜视角

如下：

其中，λ表示波长；

步骤403、方位时间序列η_i＝(i-N_a/2)/PRF+η_c,i＝1,…,N_a；其中，N_a为原始回波数据S_o的方位向长度，PRF为脉冲重复频率，η_c为雷达波束照射中心时刻；

步骤404、计算方位向包络w_a(η_i)如下：

其中，abs(·)表示求绝对值函数，

表示雷达运动过程中在方位时间η_i波束中心偏离目标的方位离轴角在斜距平面上的投影；θ_bw为方位向波束宽度；

步骤405、计算在波束照射范围内在方位时间η_i的斜距R(η_i)如下：

步骤406、计算Δk₁和Δk₂，公式如下：

Δk₁＝{R_k-min[R(R(η_i)＞0)]}/bin_r；

Δk₂＝{max[R(η_i)]-R_k}/bin_r。

6.根据权利要求1所述的方法，其特征在于，所述的步骤五中，计算距离向脉冲压缩包络p_r(τ_l)如下：

距离向采样时刻序列τ_l＝(l-N_r/2)/f_s+2R_ref/c,l＝1,…,N_r；f_s为距离采样率，R_ref为参考斜距，c为光速；

计算距离向时域包络

其中，T_r为发射脉冲时宽；

则距离向频域包络W_r(f_τ)＝w_r(f_τ/K_r)；其中，f_τ表示距离向频率，K_r为发射脉冲调频率；

对距离向频域包络进行快速傅里叶逆变换，根据下式计算距离向脉冲压缩包络p_r(τ_l)：

p_r(τ_l)＝IFFT_τ{W_r(f_τ)}

≈T_rsinc[K_rT_r(τ_l-2R_ref/c)]

其中，sinc是表示辛格函数，IFFT_τ{·}表示进行快速傅里叶逆变换。

7.根据权利要求1或6所述的方法，其特征在于，所述的步骤五中，对二维观测数据y构造的稀疏字典Ψ为：

其中，V是原始回波数据中有效回波部分的采样点数，V＜N_a；

稀疏字典Ψ的大小为N_a×(Δk*N)，其中包含的原子个数N＝N_a-V，原子大小为N_a×Δk，Δk＝Δk₁+Δk₂。

8.根据权利要求1所述的方法，其特征在于，所述的步骤六中，对观测数据y进行稀疏重构的过程包括如下：

步骤601、获得传感矩阵A，大小为M×(Δk*N)；

步骤602、将A中所有二维原子均拉成一维列向量，此时A的大小为(M*Δk)×N，用a_j表示拉长后的各原子，j＝1,2,…,N表示各原子在A中的序列号；将二维观测信号y也拉成一个一维列向量，此时y的大小为(M*Δk)×1；

步骤603、进行初始化，包括：设迭代次数为t，初始t＝1，t≤M；设重构残差为r_t，设置初始残差r₀＝y；设第t次迭代后的传感矩阵A中原子的序号集合为Λ_t，初始的序列集合Λ₀为空集；设L_y表示选择与二维观测信号y最匹配的原子步长，L_r表示选择与残差r_t最匹配的原子步长，初始设置L_y＝l，L_r＝l，l表示初始步长；设

表示二维观测信号y的稀疏表示系数估计，其大小为N×1，初始

步骤604、预测试一，获得原子初选集P，分三步进行：

第一步，筛选A中有效回波位置完全包含于二维观测数据y中有效回波位置的原子，这些原子所在A中的序号构成原子预备集R；第二步，设A_R＝{a_j},j∈R，对A_R中各原子与y进行相关性测试，计算相关度u_{R_y}＝abs[A_R ^Ty]，其中，abs[·]表示求模运算；第三步，选择u_{R_y}中所有大于阈值Th＝0.1×max(u_{R_y})的值所对应的A中原子的序号构成原子初选集P，设A_P＝{a_j},j∈P；

步骤605、预测试二，获得原子候选集C_t，分三步进行：

第一步，选择原子初选集P中前L_y个与y之间的相关度最大值对应的A中原子的序号构成集合S_yt；第二步，对A_P中各原子与残差r_t-1进行相关性测试，相关度u_{P_r}＝abs[A_P ^Tr_t-1]，选择所有u_{P_r}中前L_r个最大值对应的A中原子的序号构成集合S_rt；第三步，令原子候选集C_t＝Λ_t-1∪S_yt∪S_rt，其中，∪表示集合并运算；设

j∈C_t；

步骤606、最终测试，获得原子支撑集F_t，包括：对

中各原子按使残差最小原则测试，也就是求

的最小二乘解，θ_t表示在第t次迭代的y的稀疏表示系数，得到θ_t的估计值

选择θ_t中绝对值最大的前L_y项记为

对应

中的L_y个原子记为A_tL，这L_y个原子对应在A中的序号记为Λ_tL，令原子支撑集F_t＝Λ_tL；

步骤607、计算残差r_new＝y-A_tL(A_tL ^TA_tL)^-1A_tL ^Ty，并计算该残差能量对于初始残差能量的减少率rate_0＝(||r₀||₂-||r_new||₂)/||r₀||₂，以及相邻两次迭代残差能量的减少量对初始残差能量的比率rate_1＝(||r_t-1||₂-||r_new||₂)/||r₀||₂，其中||·||₂表示求向量的二范数运算；

步骤608、进行迭代条件判断，包括：

首先对rate_0进行判断，如果rate_0＞1-ε0，则令最终原子支撑集F_final＝F_t，对应最终稀疏表示系数

并停止迭代进入步骤609；否则，继续对rate_1进行判断；

如果rate_1＜0，则更新步长L_y的值增加l，L_r的值增加l，返回步骤605继续迭代；否则，更新步长L_y的值增加l，令L_r＝l，Λ_t＝F_t，r_t＝r_new，t自增1，并进行以下判断：若rate_1＞ε1，则令最终原子支撑集F_final＝F_t，对应最终稀疏表示系数

并返回步骤605继续迭代；若rate_1≤ε1，则停止迭代进入步骤609；ε0和ε1是设置的参数，取值接近于0；

步骤609、根据最终原子支撑集F_final和对应的最终稀疏表示系数

得到大小为N×1的稀疏表示向量

最后，在向量

两端分别补V/2个0，得到y最终的稀疏表示系数估计向量

中的非零元素及位置表示二维观测数据y所包含地面目标的后向散射信息，且为最终压缩成像数据矩阵

的第k列；

在对回波数据S_cos__rc(τ,η)的所有距离门完成从步骤四到步骤六的操作后，把每个距离门对应的稀疏表示系数估计向量按距离门序号进行排列，得到回波数据S_cos的最终压缩成像数据矩阵

9.根据权利要求8所述的方法，其特征在于，所述的步骤608中，设置参数ε0和ε1分别为：ε0＝1e-4，ε1＝1e-8。

Images