技术领域
[0001] 本
发明属于
信号处理领域，涉及一种星载SAR成像处理方法，具体涉及一种用于互质采样星载SAR的成像处理方法。
背景技术
[0002]
合成孔径雷达(Synthetic Aperture Radar，SAR)是一种主动式
微波成像
传感器，其通过发射宽带信号，结合合成孔径技术，能够在距离向和方位向上同时获得高
分辨率图像，实现对目标区域的远距离主动成像。与传统光学遥感和高
光谱遥感相比，SAR具备全天候、全天时、高分辨对地观测的能
力，还具有一定的穿透性，是实现军事侦察、自然资源普查、
自然灾害监测等的重要技术手段。
[0003] 高分辨率和宽测绘带对地观测是星载SAR的发展趋势之一。对于传统的星载SAR系统，方位空间分辨率与距离向测绘幅宽对SAR系统的脉冲重复
频率(Pluse Repeated Frequency，PRF)提出了相互矛盾的约束要求：一方面，方位高分辨率要求SAR系统的PRF足够高，以保证方位采样满足奈奎斯特采样定理；另一方面，宽测绘带要求SAR系统的PRF足够低，以保证测绘带回波在两次脉冲发射间隔时间内被完整接收。受固定PRF的限制，传统星载SAR不能同时实现高分辨率和宽测绘带对地观测。
[0004] 互质阵列和互质采样是近些年在阵列信号领域中被提出的一种新颖的稀疏非均匀阵列及采样方式。互质阵列由两个均匀线性阵列(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)估计，波束形成等领域得到广泛的关注和应用。
[0005] 互质采样星载SAR将互质采样技术引入星载SAR系统，其利用互质阵列可由低物理阵元数提供高自由度的特点，在采样率远低于奈奎斯特采样定律的条件下实现具有稀疏特性场景的精确重构，有效减小雷达系统所需要的方位采样数目。此时，若进一步采用时分复用技术，如图2所示，在非互质采样时刻将波束切换至其他测绘带区域，可有效增加SAR系统的测绘带宽度。
[0006] 由于在方位向采用了互质稀疏采样，如何实现互质采样星载SAR回波信号成像处理变得十分困难。目前对于互质采样SAR成像方法的研究中，均假设在SAR回波的距离方位耦合和二维空变特性可忽略不计的情况下，引入
压缩感知(CS)理论，构造由方位调频信号组成的稀疏字典，采用压缩感知理论中的稀疏重构
算法对各方位信号进行重构。参考文献1[Tao Y,Zhang G,Li D.Coprime sampling with deterministic digital filters in compressive 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 SAR Imaging Based on Co-Prime Arrays Sampling[C]//Advances in Materials,
Machinery,Electrical Engineering.2017.]在距离向和方位向均进行互质采样，并采用CS算法实现SAR信号的二维精确聚焦。上述几篇文献是从不同
角度探讨互质采样在SAR中的应用，有不同的侧重点，但对于互质采样SAR成像方面，他们主要侧重分析用压缩感知成像的可行性，因此都是在忽略SAR回波信号二维耦合和空变特性情况下，直接对SAR回波的一维方位向信号进行压缩感知重构，采用的重构算法也基本是压缩感知中最普遍的OMP算法，并且采用这种重构算法需要已知信号稀疏度，不能实现稀疏度自适应。
[0007] 所以针对目前由于完全忽略了SAR回波信号的二维耦合及二维空变特性，无法实现回波信号距离方位二维耦合补偿处理及成像参数空变补偿处理的情况，本发明实现了一种用于互质采样星载SAR的成像处理方法，主要侧重于互质采样SAR成像算法的精确性及实用性。
发明内容
[0008] 本发明的目的是为了解决互质采样模式下星载SAR的成像问题，提出了一种用于互质采样星载SAR的成像处理方法，能够完成距离方位二维耦合不可忽略情况下的精确聚焦处理。
[0009] 本发明的用于互质采样星载SAR的成像处理方法，实现步骤包括：
[0010] 步骤一：读入成像参数，包括互质数co1、co2以及回
波数据方位向采样点数M、回波数据距离向采样点数Nr；根据回波数据方位向采样点数M和回波数据距离向采样点数Nr读入互质采样星载SAR的回波数据Scos。
[0011] 步骤二：根据互质数co1、co2以及回波数据方位向采样点数M构造互质采样矩阵Φ。
[0012] 步骤三：对回波数据Scos进行距离向
脉冲压缩；设距离向脉冲压缩后的回波数据表示为Scos_rc(τ,η)，τ表示距离向时间，η表示方位向时间。
[0013] 针对回波数据Scos_rc(τ,η)，按照距离
门次序，从它的第一个
距离门开始，依次对每个距离门执行步骤四～步骤六。
[0014] 步骤四：二维观测信号截取；
[0015] 利用雷达的等效斜视距离模型，获取距离门k对应的雷达波束中心照射到点目标的斜距Rk，并得到相应的波束照射范围内的最短斜距与斜距Rk之差跨越的距离门数Δk1、波束照射范围内的最大斜距与斜距Rk之差跨越的距离门数Δk2；
[0016] 判断k＞Δk1及k+Δk2＜Nr是否成立，若是，在Scos_rc(τ,η)上截取第k-Δk1～k+Δk2个距离门作为二维观测数据y；否则，停止对当前距离门k的操作，k自增1，继续对距离门k执行步骤四。
[0017] 步骤五：对步骤四截取的二维观测数据y构造稀疏字典。
[0018] 首先，计算距离向脉冲压缩包络pr(τl)，τl,l＝1,…,Nr为距离向采样时刻序列；
[0019] 然后，构造二维参考信号sk_rc(τl,ηi)如下：
[0020] sk_rc(τl,ηi)＝pr[τl-2R(ηi)/c]exp{-j4πf0R(ηi)/c},τl∈[τa,τb][0021] 其中，τa＝2min[R(R(ηi)＞0)]/c，τb＝2max[R(ηi)]/c；c为光速；R(R(ηi)＞0)表示R(ηi)中所有非零值；min(·)表示求最小值函数，max(·)表示求最大值函数；ηi,i＝1,…,Na为方位时间序列，R(ηi)是在距离门k对应的波束照射范围内的斜距；f0为雷达工作频率；
[0022] 最后，将二维参考信号sk_rc(τl,ηi)在方位向进行循环移位，舍去由于循环移位导致有效回波
位置不连续的
原子，得到对应于观测数据y的稀疏字典
[0023] 步骤六：利用互质采样矩阵Φ和稀疏字典Ψ对观测数据y进行稀疏重构。在稀疏重构时，首先获得传感矩阵A＝ΦΨ，然后将A中所有二维原子均拉成一维列向量，将二维观测信号y也拉成一个一维列向量，然后对观测数据y进行稀疏重构。
[0024] 所述的步骤六中，对观测数据y进行稀疏重构的过程包括：
[0025] 步骤601、由下式获得大小为M×(Δk×N)传感矩阵A＝ΦΨ；
[0026] 步骤602、将传感矩阵A中所有二维原子均拉成一维列向量，此时A的大小为(M×Δk)×N，用aj表示拉长后的各原子，j＝1,2,…,N表示各原子在A中的序列号。将二维观测信号y也拉成一个一维列向量，此时y的大小为(M×Δk)×1；
[0027] 步骤603、进行初始化，包括：设
迭代次数为t，初始t＝1，t≤M；设重构残差为rt，设置初始残差r0＝y；设第t次迭代后的传感矩阵A中原子的序号集合为Λt，初始的序列集合Λ0为空集；设Ly表示选择与二维观测信号y最匹配的原子步长，Lr表示选择与残差rt最匹配的原子步长，初始设置Ly＝l，Lr＝l，l表示初始步长；设 表示二维观测信号y的稀疏表示系数估计，其大小为N×1，初始
[0028] 步骤604、预测试一，获得原子初选集P，分三步进行：
[0029] 第一步，筛选A中有效回波位置完全包含于二维观测数据y中有效回波位置的原子，这些原子所在A中的序号构成原子预备集R；第二步，设AR＝{aj}(forallj∈R)，对AR中各原子与y进行相关性测试，计算相关度uR_y＝abs[ARTy]，其中，abs[·]表示求模运算；第三步，选择uR_y中所有大于
阈值Th＝0.1×max(uR_y)的值所对应的A中原子的序号构成原子初选集P，设AP＝{aj}(forallj∈P)；
[0030] 步骤605、预测试二，获得原子候选集Ct，分三步进行：
[0031] 第一步，选择原子初选集P中前Ly个与y之间的相关度最大值对应的A中原子的序号构成集合Syt；第二步，对AP中各原子与残差rt-1进行相关性测试，相关度uP_r＝abs[APTrt-1]，选择所有uP_r中前Lr个最大值对应的A中原子的序号构成集合Srt；第三步，令原子候选集Ct＝Λt-1∪Syt∪Srt，其中，∪表示集合并运算；设
[0032] 步骤606、最终测试，获得原子
支撑集Ft，包括：对 中各原子按使残差最小原则测试，也就是求 的最小二乘解，θt表示在第t次迭代的y的稀疏表示系数，得到θt的估计值 选择θt中绝对值最大的前Ly项记为 对应 中的Ly个原子记为AtL，这Ly个原子对应在A中的序号记为ΛtL，令原子支撑集Ft＝ΛtL；
[0033] 步骤607、计算残差rnew＝y-AtL(AtLTAtL)-1AtLTy，并计算该残差
能量对于初始残差能量的减少率rate_0＝(||r0||2-||rnew||2)/||r0||2，以及相邻两次迭代残差能量的减少量对初始残差能量的比率rate_1＝(||rt-1||2-||rnew||2)/||r0||2，其中||·||2表示求向量的二范数运算；
[0034] 步骤608、进行迭代条件判断，包括：
[0035] 首先对rate_0进行判断，如果rate_0＞1-ε0，则令最终原子支撑集Ffinal＝Ft，对应最终稀疏表示系数 并停止迭代进入步骤609；否则，继续对rate_1进行判断；
[0036] 如果rate_1＜0，则更新步长Ly＝Ly+l，Lr＝Lr+l，返回步骤605继续迭代；否则，更新步长Ly＝Ly+l，令Lr＝l，Λt＝Ft，rt＝rnew，t自增1，并进行以下判断：若rate_1＞ε1，则令最终原子支撑集Ffinal＝Ft，对应最终稀疏表示系数 并返回步骤605继续迭代；若rate_1≤ε1，则停止迭代进入步骤609；ε0和ε1是设置的参数，取值接近于0；
[0037] 步骤609、根据最终原子支撑集Ffinal和对应的最终稀疏表示系数 得到大小为N×1的稀疏表示向量 最后，在向量 两端分别补V/2个0，得到y最终的稀疏表示系数估计向量 中的非零元素及位置表示二维观测数据y所包含地面目标
的后向散射信息，且为最终压缩成像数据矩阵 的第k列；
[0038] 在对回波数据Scos_rc(τ,η)的所有距离门完成从步骤四到步骤六的操作后，把每个距离门对应的稀疏表示系数估计向量按距离门序号进行排列，得到回波数据Scos的最终压缩成像数据矩阵
[0039] 本发明与
现有技术相比，具有以下优势：
[0040] (1)
精度高。本发明根据距离单元徙动数截取二维观测信号，并通过所改进的稀疏度自适应的二维信号稀疏重构算法，能够完成距离方位二维耦合不可忽略情况下的精确聚焦处理。
[0041] (2)全场景聚焦。本发明根据各距离门的多普勒参数计算截取观测信号的距离门个数以及构造稀疏字典，补偿了成像参数随距离门变化对稀疏重构造成的影响，由此可实现全场景精确聚焦。
[0042] (3)实用性强。本发明能够实现互质采样工作模式下的成像处理，相比于目前忽略了距离方位耦合补偿和成像参数空变补偿的成像处理方法，本发明具有更好的实用性。
附图说明
[0043] 图1是互质阵列结构示意图；
[0044] 图2是互质采样星载SAR工作方式示意图；
[0045] 图3是本发明互质采样星载SAR成像处理方法的
流程图；
[0046] 图4是雷达的等效斜视距离模型的示意图；
[0047] 图5是本发明中二维信号稀疏度自适应的稀疏重构算法的流程图；
[0048] 图6是本发明多点目标成像仿真处理结果示意图，(a)为多点目标成像仿真结果；(b)为(a)中右下角点目标的
水平剖面示意图；(c)为在距离向脉冲压缩时加入Kaiser窗的多点目标成像仿真结果；(d)为(c)中右下角点目标的水平剖面示意图。
具体实施方式
[0049] 为了便于本领域普通技术人员理解和实施本发明，下面结合附图对本发明作进一步的详细和深入描述。
[0050] 本发明的一种用于互质采样星载SAR的成像处理方法，在获得成像参数、回波数据及互质采样矩阵之后，首先对方位互质采样SAR回波信号进行距离向脉冲压缩，然后从第一个距离门开始，依据各距离门的多普勒参数计算该距离门回波信号跨越的距离门个数以及构造相应的稀疏字典，并以此为基准截取二维观测数据，最后用改进的二维观测信号稀疏度自适应的稀疏重构算法对二维观测信号进行重构，由此得到场景目标的后向散射信息。
[0051] 本发明方法整体上主要包括计算互质采样矩阵、距离向脉冲压缩、二维观测信号截取、稀疏字典构造和二维观测信号稀疏重构五个部分。本发明
实施例中，设置9个单位幅度点目标，显然单个点目标在方位向的稀疏度为1。根据所设置的多点目标进行互质采样星载SAR回波仿真，得到回波数据，进行成像处理，本发明方法流程如图3所示，下面说明各步骤。
[0052] 步骤一：读入相关的成像参数和需要进行成像处理的互质采样星载SAR的回波数据，具体包括步骤101和102。
[0053] 步骤101、读入成像参数，包括：雷达工作频率f0、发射脉冲时宽Tr、发射脉冲调频率Kr、距离采样率fs、脉冲重复频率PRF、互质数co1和co2、方位向波束宽度θbw、参考斜距Rref、雷达波束照射中心时刻ηc、回波数据方位向采样点数M、回波数据距离向采样点数Nr、各个距离门k的多普勒中心频率fD_k、各个距离门k的多普勒调频率fR_k，k为距离门编号，k∈{1,2,…,Nr}。
[0054] 步骤102、根据回波数据方位向采样点数M和回波数据距离向采样点数Nr，读入需要进行成像处理的互质采样星载SAR的回波数据Scos，其大小为M×Nr。
[0055] 其中，本实施例中具体成像参数为：f0＝10GHz，Tr＝40μs，Kr＝1.5MHz/μs，fs＝72MHz，PRF＝5680Hz，co1＝4，co2＝5，θbw＝0.0089，Rref＝850km，ηc＝-0.8358s，M＝2913，Nr＝4096，fD_k和fR_k与距离门k有关，如k＝2048时，fD_k＝3306.8Hz，fR_k＝-3956.3Hz/s。
[0056] 步骤二：根据互质数co1、co2以及回波数据方位向采样点数M构造互质采样矩阵，具体包括如下子步骤201～203。
[0057] 步骤201、计算原始回波数据So的方位向长度Na，如下：
[0058] rate＝(co1+co2-1)/co1/co2 (1)
[0059]
[0060] 其中，rate表示互质采样的降采样率， 表示向上取整。原始回波数据So表示未进行方位互质采样时的回波数据，So是未知的，大小为Na×Nr，在实际中直接获得互质采样后的回波数据Scos，此处目的是为了计算So的方位向长度Na。
[0061] 在本实施例中，计算得到rate＝0.4，Na＝7283。
[0062] 步骤202、根据互质数co1、co2计算互质采样时刻序列H，如下：
[0063]
[0064] 其中，H中元素表示互质采样时刻，H中元素按升序排列且无重复元素， 表示向下取整。
[0065] 步骤203、构造互质采样矩阵Φ如下：
[0066]
[0067]
[0068] 其中，m＝1,2,…,M，Hm表示互质采样时刻序列H中第m个元素。向量 是一个分段序列，表示为{Xn}，当n＝Pm+1时，Xn为1，当n＝1,2,…,Na,且n≠Hm+1时，Xn为0。
[0069] 由步骤二构造的互质采样矩阵Φ大小为M×Na，满足Scos＝ΦSo。
[0070] 步骤三：对互质采样星载SAR的回波数据Scos进行距离向脉冲压缩。本发明采用传统匹配滤波的脉冲压缩方式对互质采样星载SAR回波信号进行距离向聚焦，包括如下子步骤301～303。
[0071] 步骤301、互质采样回波数据Scos在距离向时间域的表示为Scos(τ,η)，τ表示距离向时间，η表示方位向时间，对在距离向时域的回波数据Scos(τ,η)进行距离向快速傅里叶变换得到距离向频域的回波数据Scos(fr,η)如下：
[0072] Scos(fr,η)＝FFTτ{Scos(τ,η)} (6)
[0073] 其中，fr表示距离向频率，FFTτ{·}表示进行距离向快速
傅立叶变换。
[0074] 步骤302、在距离向频域构造如(7)式所示的滤波器H(fr)：
[0075] H(fr)＝exp{jπfr2/Kr} (7)
[0076] 其中，j表示虚数单位。
[0077] 步骤303、将滤波器H(fr)与变换到距离向频域的回波数据Scos(fr,η)相乘，再利用快速傅里叶逆变换将相乘的结果变换到距离向时域，得到已完成距离向聚焦的信号Scos_rc(τ,η)：
[0078] Scos_rc(τ,η)＝IFFTτ{Scos(fr,η)H(fr)} (8)
[0079] 其中，IFFTτ{·}表示进行距离向快速傅里叶逆变换。
[0080] 下面将从已完成距离向聚焦的信号Scos_rc的第一个距离门开始，按照距离门次序，逐次对每个距离门完成步骤四～步骤六。
[0081] 步骤四：截取二维观测信号。
[0082] 如图4所示，X轴表示雷达平台飞行方向，Y表示垂直于地面方向，Z表示地面垂直于雷达平台飞行方向，波束中心时刻ηc表示雷达波束中心正好照射到目标Tar的时刻。对于第k个距离门而言，波束中心时刻雷达平台到地面目标间的斜距为Rk，对应的等效速度和等效斜视角分别为Vk和 R(ηi)表示在ηi时刻雷达平台到地面目标间的斜距，ηi,i＝1,2,…Na代表方位时间序列。设Δk1表示波束照射范围内的最短斜距与斜距Rk之差跨越的距离门数，Δk2就表示波束照射范围内的最大斜距与斜距Rk之差跨越的距离门数，设Δk＝Δk1+Δk2。如果k＞Δk1，且k+Δk2＜Nr，则在已完成距离向聚焦的信号Scos_rc上截取第k-Δk1～k+Δk2个距离门信号作为二维观测信号y，y大小为M×Δk；否则，距离门k＝k+1，对下一个距离门执行步骤四～六。
[0083] 其中，计算Δk1和Δk2的步骤如下：
[0084] 步骤401、计算距离门k对应的斜距Rk：
[0085] Rk＝Rref+(k-Nr/2)bin_r (9)
[0086] bin_r＝c/2/fs (10)
[0087] 其中，bin_r表示距离采样间隔对应的斜距变化量，c为光速，取3×108m/s。
[0088] 步骤402、根据距离门k对应的多普勒中心频率fD_k和多普勒调频率fR_k计算等效速度Vk和等效斜视角
[0089]
[0090]
[0091] 其中，λ表示波长，Vk表示距离门k对应的等效速度， 表示距离门k对应的等效斜视角。
[0092] 步骤403、计算方位时间序列ηi公式如下：
[0093] ηi＝(i-Na/2)/PRF+ηc,i＝1,…,Na (13)
[0094] 在不同的方位时间ηi，对应一个等效斜距R(ηi)。
[0095] 步骤404、计算方位向包络wa(ηi)公式如下：
[0096]
[0097]
[0098] 其中，abs(·)表示求绝对值函数， 表示雷达在运动过程中在方位时间ηi时、波束中心偏离目标的方位离轴角在斜距平面上的投影。
[0099] 步骤405、计算在波束照射范围内的斜距R(ηi)公式如下：
[0100]
[0101] 步骤406、计算Δk1和Δk2，公式如下：
[0102] Δk1＝{Rk-min[R(R(ηi)＞0)]}/bin_r (17)
[0103] Δk2＝{max[R(ηi)]-Rk}/bin_r (18)
[0104] 其中，R(R(ηi)＞0)表示R(ηi)中所有非零值，min(·)表示求最小值函数，max(·)表示求最大值函数。
[0105] 本发明截取的观测信号为二维信号，可将SAR回波的二维耦合作为一个整体进行重构，避免了在方位稀疏非均匀采样情况下难以实现的二维解耦操作；而且截取观测信号的距离门数根据回波各距离门的多普勒参数进行更新，使得回波时频分布特性随距离门变化对重构算法产生的影响得以消除。
[0106] 步骤五：根据步骤四中截取的距离门k的二维观测信号y，构造稀疏字典，具体如下三个子步骤501～503。本发明方法中通过等效斜距模型构造参考信号，使得稀疏字典的精确度提高，从而提高重构精度；在下面构造稀疏字典时舍去了有效回波位置不连续的原子，使得在重构过程中可以避免非有效原子的干扰，同时还可以避免额外的计算开销。
[0107] 步骤501、计算距离向采样时刻序列τl公式如下：
[0108] τl＝(l-Nr/2)/fs+2Rref/c,l＝1,…,Nr (19)
[0109] 其中，l表示采样步长。
[0110] 计算距离向时域包络wr(τl)公式如下：
[0111]
[0112] 计算距离向频域包络Wr(fτ)公式如下：
[0113] Wr(fτ)＝wr(fτ/Kr) (21)
[0114] 其中，fτ表示距离向频率。
[0115] 计算距离向脉冲压缩包络pr(τl)公式如下：
[0116]
[0117] 其中，sinc是表示辛格函数。
[0118] 步骤502、构造二维参考信号sk_rc(τl,ηi)如下：
[0119] sk_rc(τl,ηi)＝pr[τl-2R(ηi)/c]exp{-j4πf0R(ηi)/c},τl∈[τa,τb] (23)[0120] 其中，τa＝2min[R(R(ηi)＞0)]/c，τb＝2max[R(ηi)]/c；R(ηi)为波束照射范围内的斜距，由(16)式计算。
[0121] 步骤503、将二维参考信号sk_rc(τl,ηi)在方位向进行循环移位，得到对应于二维观测信号y的稀疏字典。若单点目标原始回波中有效回波的方位向采样点数为V，V是互质采样SAR回波数据的原始回波数据(方位长度为Na)中有效回波部分的采样点数，V＜Na，舍去由于循环移位导致有效回波位置不连续的原子，得到稀疏字典Ψ，如公式(24)所示，其大小为Na×(Δk×N)，其中包含的原子个数N＝Na-V，各原子大小为Na×Δk。
[0122]
[0123] 步骤六：根据改进的稀疏度自适应
匹配追踪算法对二维观测信号进行稀疏重构。
[0124] 本发明将稀疏度自适应匹配追踪算法(Sparsity Adaptive Matching Pursuit，SAMP)(参考文献[4]：Do T T,Lu G,Nguyen N,et al.Sparsity adaptive matching pursuit algorithm for practical compressed sensing[C]//Signals,Systems and Computers,2008,Asilomar Conference on.IEEE,2009:581-587.)进行改进，使其适用于步骤四中所得二维观测信号y的稀疏重构。
[0125] 本步骤对y进行稀疏重构的流程如图5所示，包括如下子步骤601～609。
[0126] 步骤601、由公式(25)获得大小为M×(Δk×N)的传感矩阵A：
[0127] A＝ΦΨ (25)
[0128] 步骤602、将传感矩阵A中所有二维原子均拉成一维列向量，此时A的大小为(M×Δk)×N，用aj表示拉长后的各原子，j＝1,2,…,N表示各原子在A中的序列号。将二维观测信号y也拉成一个一维列向量，此时y的大小为(M×Δk)×1。
[0129] 本发明方法首先将二维观测信号y、传感矩阵A的所有二维原子均拉长成一维列向量，通过这样简单的拉长操作，使一维信号的稀疏重构算法适用于二维信号的稀疏重构。
[0130] 步骤603、初始化：设重构残差为rt，设置初始残差r0＝y；设迭代次数为t，初始t＝1，t≤M；设置第t次迭代后的传感矩阵A中原子的序号集合为Λt，第一次迭代前初始的序列集合 表示空集；设Ly表示选择与二维观测信号y最匹配的原子步长，Lr表示选择与残差rt最匹配的原子步长，初始设置Ly＝l，Lr＝l，l表示初始步长，本实施例中，l＝1；设 表示二维观测信号y的稀疏表示系数估计，其大小为N×1，初始
[0131] 步骤604、预测试一，获得原子初选集P：
[0132] 分三小步进行，第一步，筛选A中有效回波位置完全包含于二维观测信号y中有效回波位置的原子，这些原子所在A中的序号构成原子预备集R；第二步，设AR＝{aj}(forallj∈R)，对AR中各原子与二维观测信号y进行相关性测试，即计算uR_y＝abs[ARTy]，其中，abs[·]表示求模运算，uR_y表示AR中各原子与观测信号y之间的相关度；第三步，选择uR_y中所有大于阈值Th＝0.1×max(uR_y)的值，这些值对应A中原子的序号构成原子初选集P，这些值就是原子初选集P中元素与y之间的相关度uP_y，设AP＝{aj}(forallj∈P)。后续的迭代过程将在AP中进行。
[0133] 通过原子初选集的筛选，可以排除有效回波位置不匹配的原子，大大减小了重构算法不必要的计算开销，提高重构效率。
[0134] 步骤605、预测试二，获得原子候选集Ct：
[0135] 也分三小步进行，第一步，选择原子初选集P中前Ly个与y的相关度最大值uP_y对应A中原子的序号构成集合Syt；第二步，对AP中各原子与残差rt-1进行相关性测试，即计算uP_r＝abs[APTrt-1]，选择所有uP_r中前Lr个最大值对应的A中原子的序号构成集合Srt；第三步，令原子候选集Ct＝Λt-1∪Syt∪Srt，其中，符号∪表示集合并运算。
[0136] 原子候选集的选取中，通过分别测试原子初选集与观测信号的相关度以及原子初选集与残差的相关度进行筛选，且选择与观测信号最匹配原子的步长和选择支撑集的步长一致，并随迭代次数递增，但选择与残差最匹配原子的步长只随迭代阶段递增。
[0137] 步骤606、最终测试，获得原子支撑集Ft：
[0138] 对 中各原子按使残差最小原则测试，即求 的最小二乘解： 选择θt中绝对值最大的前Ly项记为 对应 中的Ly个原
子记为AtL，这Ly个原子对应在A中的序号记为ΛtL，并令原子支撑集Ft＝ΛtL；θt表示二维观测信号y的在第t次迭代的稀疏表示系数，θt表示在第t次迭代的稀疏表示系数估计值。
[0139] 步骤607、计算残差rnew＝y-AtL(AtLTAtL)-1AtLTy，并计算此残差能量对于初始残差能量的减少率rate_0＝(||r0||2-||rnew||2)/||r0||2，以及相邻两次迭代残差能量的减少量对初始残差能量的比率rate_1＝(||rt-1||2-||rnew||2)/||r0||2，其中||·||2表示求向量的二范数运算。
[0140] 步骤608、判断迭代条件：
[0141] 首先对rate_0进行判断：
[0142] 如果rate_0＞1-ε0，则令最终原子支撑集Ffinal＝Ft，对应最终稀疏表示系数并停止迭代进入步骤609；
[0143] 否则，则对rate_1进行判断：
[0144] 如果rate_1＜0，则更新步长Ly＝Ly+l，Lr＝Lr+l，返回步骤605继续迭代；
[0145] 否则，则更新步长Ly＝Ly+l，令Lr＝l，Λt＝Ft，rt＝rnew，t＝t+1。并且进行以下判断：若rate_1＞ε1，则令最终原子支撑集Ffinal＝Ft，对应最终稀疏表示系数 并返回步骤605继续迭代；若rate_1≤ε1，则停止迭代进入步骤609。
[0146] 在本发明实施例中，取ε0＝1e-4，ε1＝1e-8。
[0147] 需要说明，ε0、ε1均取接近于0的值。其中，ε0为判断点目标重构的迭代停止条件参数。点目标重构成功时残差能量接近于0，但对面目标来说，所截取的观测信号为目标所有散射中心回波在所截取距离门范围内的
叠加结果，在迭代过程中残差能量不会减小至0，在本发明中一个稳妥的方法是取ε0＝1e-4。ε1为判断面目标重构的迭代停止条件参数，即相邻两次迭代的残差减小量对初始残差的比率足够小时停止迭代。ε1取值越接近于0，重构的精度就越高，但不会超过本发明中稀疏重构算法的理论重构精度，即方位向采样间隔。而另一方面，如果ε1过小，将导致不必要的算法计算开销，综合考虑以及进行实验，取ε1＝1e-8最为合适。
[0148] 由于本发明中二维观测信号的截取以及稀疏字典的构造方式，使得面目标的二维观测信号在迭代过程中残差能量不能减小至零或者无限接近噪声功率， 因此本发明在判断残差能量减小率的
基础上，引入对连续两次迭代的残差能量减小量对初始残差能量的比率的判断条件，即使连续两次迭代之间的残差能量减小量在小于某个阈值时停止迭代。
[0149] 步骤609、根据最终原子支撑集Ffinal和对应的最终稀疏表示系数 得到大小为N×1的稀疏表示向量 最后，在向量 两端分别补V/2个0，得到二维观测信号y最终的稀疏表示系数估计向量
[0150] 稀疏表示系数估计向量 中的非零元素及位置即表示二维观测信号y所包含地面目标的后向散射信息，且为最终压缩成像数据矩阵 的第k列。当对Scos_rc所有距离门都完成从步骤四到步骤六的所有操作后，把每个距离门对应的稀疏表示系数估计向量按距离门序号进行排列，就可以得到回波数据Scos的最终压缩成像数据矩阵
[0151] 本实施例的成像仿真结果如图6所示，其中图6(a)为本实施例中的多点目标成像仿真结果，可以看出仅利用原始回波40％的数据量就完成了方位向的精确聚焦，且没有方位向旁瓣。图6(b)为(a)中右下角点目标成像结果在距离向插值后的旁瓣情况，发现距离向峰值旁瓣比PSLR＜-13dB，这说明，本发明中的稀疏重构方法没有影响传统SAR距离
多普勒成像方法在距离向的分辨率。此外，图6(c)显示了在距离向脉冲压缩时加入Kaiser窗以抑制距离旁瓣的多点目标成像仿真结果，图6(d)同样为(c)中右下角点目标成像结果在距离向插值后的旁瓣情况，可发现在加入Kaiser窗以后距离向峰值旁瓣比PSLR＜-60dB。
