CN109001732A - 一种优化的压缩感知步进频sar成像恢复重建方法 - Google Patents

Abstract

本发明提供了一种优化的压缩感知步进频SAR成像恢复重建方法，涉及雷达成像技术领域，利用克罗内克积的运算对应关系，计算列向量后更新传感矩阵中的列集合，再利用最小二乘法估计在当前列集合Φ下的目标信号，更新残余信号矢量，直至满足终止条件，则按照建立的索引位置对目标场景散射分布赋值，即可得散射率分布二维目标像。本发明由于根据步进频雷达回波信号的稀疏表示模型，利用感知矩阵和稀疏字典的二维可分离性，优化了在传感矩阵中需找相关性最强列的方法，降低了存储量需求，提高了计算效率。

Description

一种优化的压缩感知步进频SAR成像恢复重建方法

技术领域

本发明涉及雷达成像技术领域，尤其是一种基于压缩感知的SAR成像方法。

背景技术

高分辨雷达成像系统在军用和民用领域有着广泛的应用，但受奈奎斯特采样定理的限制，传统雷达在提高分辨率和满足实时性要求时面临采样率过高、数据量过大、快速处理困难等问题的挑战。具体来讲，主要有宽带雷达中A/D转换技术的限制，匹配滤波对分辨率的影响，不确定原理对时频分辨的约束，雷达数据本身不健全及数据存储量和传输量过大等问题。高频区雷达目标散射特性可由少数重要散射中心来刻画，其回波具有潜在的稀疏性，因而压缩感知对雷达成像处理具有极大的吸引力。基于压缩感知的雷达成像技术能够有效缓解高分辨雷达系统的数据采集和存储压力，但是付出的代价是信号重构算法的软件成本。正交匹配追踪(OMP)是最早的恢复重建贪婪类算法之一，从现有的研究来看，对OMP计算量的改善通常以损失些许重构性能为代价，而性能的改善则以较高的运算量或存储量为成本。为进一步提高性能，Fu Ning等人提出了一种基于相关系数的正交匹配追踪算法(Amodified orthogonal matching algorithm using correlation coefficient forcompressed sensing.I2MTC,May 10-12,2011:1-5)利用感知矩阵和测量矢量间的相关系数取代内积估计信号的支撑域，虽然提高了重构成功的概率，但运算量约为经典OMP的3倍，计算复杂度过高。

发明内容

为了克服现有技术的不足，依据雷达稀疏回波字典及测量矩阵的二维可分离特性，本发明提出了一种优化的压缩感知步进频SAR成像恢复重建方法，较大幅度提高了计算效率。

本发明解决其技术问题所采用的技术方案的详细步骤如下：

步骤一：输入初始值以及循环终止条件，即如下步骤：

1)对原始雷达回波数据做压缩采样所得的数据矩阵Y；

2)根据已知的包括脉冲时宽、距离向采样率、方位向采样率、成像场景区域大小和网格划分数的雷达参数先验知识构造距离向传感矩阵Θ_r和方位向传感矩阵Θ_d；

3)输入如下循环终止条件的初始值：

①输入目标的等效散射中心数目，即目标散射率分布的稀疏度I；

②输入残余信号的能量阈值ε；

步骤二：循环程序执行前，参数进行初始化：

1)将数据矩阵Y拉成列向量y＝vec(Y)，其中vec表示将矩阵按列整合为一个列向量，y即表示观测值；

2)残余信号矢量r赋初值，r＝y，即残余信号矢量初始值赋为观测值；

3)列向量索引集合Λ赋初值为空集；

4)列向量集合Φ赋初值为空集；

5)设循环标识k＝1，在循环中循环标识将逐次加1，直至满足循环终止条件；

步骤三：寻找传感矩阵Θ中与残余信号矢量r相关性最大的列对应的索引λ_k，此处的相关性通过传感矩阵Θ的列向量与残余信号矢量r的内积值来表征，计算Θ中所有列向量分别和r的内积，并记录最大内积值对应的Θ中的列向量，即数学表达式为其中φ_j表示矩阵Θ的第j列，argmax(f(x))表示使得f(x)取得最大值的的自变量x；

所述内积值的计算即计算Θ^Hr得到内积向量，再挑选内积向量中绝对值的最大值，具体实现方法为：

首先将残余信号矢量r按列重排为K×L的矩阵R，其中K为距离向压缩采样点数，L为方位向压缩采样点数，利用如下公式计算内积向量：

其中表示克罗内克积，得到传感矩阵Θ各个列向量和r的内积向量；

步骤四：利用克罗内克积的运算对应关系，即由索引关系式λ_k＝(q-1)×P+p，得出q和p的值，其中P和Q是成像场景划分为P×Q的二维网格后的网格数，其中，p＝1,2,3,...,P，q＝1,2,3,...,Q，计算可得Θ中对应的列向量更新索引集Λ＝Λ∪{λ_k}，将元素λ_k并入索引集合Λ，更新传感矩阵中的列集合Φ＝[Φ,φ_k]，将列向量φ_k并入列集合Φ中，即列向量集合Φ是由索引集合Λ所指示的Θ中的列向量构成的矩阵；

步骤五：利用最小二乘法估计在当前列集合Φ下的目标信号其中Φ^◇表示矩阵Φ的伪逆，计算公式为

步骤六：更新残余信号矢量令循环标识k加1，进入步骤七；

步骤七：判断是否满足k＞I或||r||₂＜ε，其中，ε为残余信号的能量阈值，若不满足上述条件之一，则跳转到步骤三，再次执行步骤三至步骤六；若至少满足以上两个条件之一则跳出循环，进入步骤八；

步骤八：按照建立的索引位置对赋值，其中，是重建得到的目标场景散射分布，它表示场景中不同网格点处的散射强度；Λ是步骤四中更新得到的索引集合，因此表示由索引集合Λ所指定的中的位置；是步骤五中计算得出的最小二乘解，然后将按列重排为大小P×Q的矩阵即为所求的目标散射率分布二维目标像。

本发明的有益效果在于由于根据步进频雷达回波信号的稀疏表示模型，利用感知矩阵和稀疏字典的二维可分离性，优化了在传感矩阵中需找相关性最强列的方法，降低了存储量需求，提高了计算效率。例如，设距离向和方位向压缩采样的点数均为100，离散网格数均为200，按照double数据类型计算，存储传感矩阵Θ所需要的内存容量就达到了2.98Gbyte，而按照存储Θ_d和Θ_r计算仅需内存容量0.305Mbyte。

附图说明

图1是本发明的整体操作流程图。

图2是本发明的飞机点散射仿真模型。

图3是本发明的原始正交匹配追踪方法的成像结果。

图4是本发明的采用本发明方法的成像结果。

具体实施方式

下面结合附图和实施例对本发明进一步说明。

图1为本发明所述方法的整体操作流程图，以此为依据，利用图2所示的飞机点散射模型生成的回波数据进行成像实验。系统参数设置为：雷达中心频率为10GHz，带宽为1GHz，距离向采样点数N为256个，目标稀疏度I为46，方位积累角度为8度，方位向采样点数M为256个。压缩采样比为4，距离向压缩采样点数K为64个，方位向压缩采样点数L为64个。距离向离散网格划分数P为128个，方位向离散网格划分数Q为128个。

步骤一：输入初始值以及循环终止条件，即如下步骤：

1)对原始雷达回波数据做压缩采样所得的数据矩阵Y；

2)根据已知的包括脉冲时宽、距离向采样率、方位向采样率、成像场景区域大小和网格划分数的雷达参数先验知识构造距离向传感矩阵Θ_r和方位向传感矩阵Θ_d；

3)输入如下循环终止条件的初始值：

①输入目标的等效散射中心数目，即目标散射率分布的稀疏度I；

②输入残余信号的能量阈值ε；

步骤二：循环程序执行前，参数进行初始化：

1)将数据矩阵Y拉成列向量y＝vec(Y)，其中vec表示将矩阵按列整合为一个列向量，y即表示观测值；

2)残余信号矢量r赋初值，r＝y，即残余信号矢量初始值赋为观测值；

3)列向量索引集合Λ赋初值为空集；

4)列向量集合Φ赋初值为空集；

5)设循环标识k＝1，在循环中循环标识将逐次加1，直至满足循环终止条件；

步骤三：寻找传感矩阵Θ中与残余信号矢量r相关性最大的列对应的索引λ_k，此处的相关性通过传感矩阵Θ的列向量与残余信号矢量r的内积值来表征，计算Θ中所有列向量分别和r的内积，并记录最大内积值对应的Θ中的列向量，即数学表达式为其中φ_j表示矩阵Θ的第j列，argmax(f(x))表示使得f(x)取得最大值的的自变量x；

所述内积值的计算即计算Θ^Hr得到内积向量，再挑选内积向量中绝对值的最大值，具体实现方法为：

首先将残余信号矢量r按列重排为K×L的矩阵R，其中K为距离向压缩采样点数，L为方位向压缩采样点数，利用如下公式计算内积向量：

其中表示克罗内克积，得到传感矩阵Θ各个列向量和r的内积向量；

步骤四：利用克罗内克积的运算对应关系(譬如对而言，A有2列B有10列,Φ的各列依次为：),即由索引关系式λ_k＝(q-1)×P+p，得出q和p的值，其中P和Q是成像场景划分为P×Q的二维网格后的网格数，其中，p＝1,2,3,...,P，q＝1,2,3,...,Q，计算可得Θ中对应的列向量更新索引集Λ＝Λ∪{λ_k}，将元素λ_k并入索引集合Λ，更新传感矩阵中的列集合Φ＝[Φ,φ_k]，将列向量φ_k并入列集合Φ中，即列向量集合Φ是由索引集合Λ所指示的Θ中的列向量构成的矩阵；

步骤五：利用最小二乘法估计在当前列集合Φ下的目标信号其中Φ◇表示矩阵Φ的伪逆，计算公式为

步骤六：更新残余信号矢量令循环标识k加1，进入步骤七；

步骤七：判断是否满足k＞I或||r||₂＜ε，其中，ε为残余信号的能量阈值，若不满足上述条件之一，则跳转到步骤三，再次执行步骤三至步骤六；若至少满足以上两个条件之一则跳出循环，进入步骤八；

步骤八：按照建立的索引位置对赋值，其中，是重建得到的目标场景散射分布，它表示场景中不同网格点处的散射强度；Λ是步骤四中更新得到的索引集合，因此表示由索引集合Λ所指定的中的位置；是步骤五中计算得出的最小二乘解，然后将按列重排为大小P×Q的矩阵即为所求的目标散射率分布二维目标像。

实施例：

步骤一：按照以上参数设置，通过仿真得到完整的回波数据，并按照4倍压缩比得到压缩采样的回波数据矩阵Y，其大小为64×64。目标散射率分布的稀疏度K＝46，以及残余信号的能量阈值ε。传感矩阵其中表示克罗内克积，构造Θ_d，其大小为64×128，Θ_d，其大小为64×128。

具体理论如下：第l个视角下压缩测量模型为其中l＝1,2,…,L，对应的是采样的视角顺序，表示距离向采样矩阵，Ψ^l表示第l个视角下的稀疏字典。令y＝[y¹,y²,…,y^l]，Φ＝[Φ¹,Φ²,…,Φ^l]，Ψ＝[Ψ¹,Ψ²,…,Ψ^l]，则联合压缩采样模型为y＝ΦΨσ＝Θσ，当矩阵Θ满足RIP准则或不相干条件时，通过非线性优化可重构σ。

尽管在每个观测视角下采用不同的测量矩阵有潜力获得更多目标信息，但在每个视角下采用相同的测量矩阵可以降低存储要求且更易于实现。进一步压缩测量模型可写为，其中Θ_d＝Φ_dΨ_d，Θ_r＝Φ_rΨ_r。

步骤二：系统循环迭代初始化，将压缩采样的数据矩阵按列堆叠转换为一维列向量，y＝vec(Y)，残余信号初值即为y，r＝y，表示从传感矩阵选择的列序号的索引集初始化为空，Λ＝空集，传感矩阵中选择出的列集合Φ，循环标识k＝1；

贪婪类算法通过每次迭代时选择一个局部最优解来逐步逼近原始信号。而传感矩阵Θ中列数PQ和列向量的长度KL均很大，因此在其中搜索与残余信号相关性最强的列向量时计算复杂度较高。所以本发明的关键在于步骤三和步骤四，避免了大型传感矩阵Θ的存储和计算，只需存储Θ_r和Θ_d，从而显著降低了存储量，使得算法效率大大提高。

步骤三：寻找传感矩阵Θ中与残余信号矢量r相关性最大的列对应的索引λ_k，即其中φ_j表示矩阵Θ的第j列。

具体方法是将残差向量r按列重排为K×L的矩阵R，通过关系式计算；

步骤四：由索引关系λ_k＝(q-1)×P+p，得到Θ中对应的列向量更新索引集Λ＝Λ∪{λ_k}和传感矩阵中的列集合Φ＝[Φ,φ_k]，即Φ是由Λ所指示的Θ中的列向量构成的矩阵；

以下步骤与经典正交匹配追踪算法处理方法相同，按照最小二乘法估计当前解，并以此更新残差信号，在本次验证中按照已知稀疏度I＝46设置终止条件。

步骤五：利用最小二乘法估计在当前列集合Φ下的目标信号其中Φ^◇表示矩阵Φ的伪逆。具体方法为

步骤六：更新残余信号分量令k＝k+1；

步骤七：判断是否满足k＞46，若满足则停止循环；若不满足，则执行步骤三；

步骤八：按照建立的索引位置对赋值，然后将按列重排为大小P×Q的矩阵即为所求的目标散射率分布二维目标像。

图2为实验生成回波所用的飞机散射模型，其稀疏度为46。图3为原始正交匹配追踪算法的成像结果，图4为本发明方法的成像结果，从成像效果上来看，本发明的方法具有和原始正交匹配追踪方法相似质量的成像结果。而本发明由于优化了在传感矩阵中需找最相关列向量的步骤，有效提高了步进频雷达系统的数据利用率。

从存储空间上来看，二维联合模型下直接采用原始正交匹配追踪方法需要的存储量为O(KLPQ)，而本发明仅需要O(KP)+O(LQ)大小的存储空间，有效节约了内存需求。

从计算复杂度上来看，贪婪类算法的计算负担主要来自残余信号与传感矩阵列向量间的内积计算和伪逆的求解过程。原始正交匹配追踪方法通过矩阵矢量积Θ^Tr计算内存最大值，每次迭代的时间是O(KLPQ)，总的运算复杂度为O(KLPQI)，而本发明方法每次迭代的时间为O(KLP+QLP)，总的运算复杂度为O(ILP(K+Q))。

综上，本发明的方法相比原始方法节约了大量的存储空间和运算开销。

Claims (1)

1.一种优化的压缩感知步进频SAR成像恢复重建方法，其特征在于包括下述步骤：步骤一：输入初始值以及循环终止条件，即如下步骤：

1)对原始雷达回波数据做压缩采样所得的数据矩阵Y；

2)根据已知的包括脉冲时宽、距离向采样率、方位向采样率、成像场景区域大小和网格划分数的雷达参数先验知识构造距离向传感矩阵Θ_r和方位向传感矩阵Θ_d；

3)输入如下循环终止条件的初始值：

①输入目标的等效散射中心数目，即目标散射率分布的稀疏度I；

②输入残余信号的能量阈值ε；

步骤二：循环程序执行前，参数进行初始化：

1)将数据矩阵Y拉成列向量y＝vec(Y)，其中vec表示将矩阵按列整合为一个列向量，y即表示观测值；

2)残余信号矢量r赋初值，r＝y，即残余信号矢量初始值赋为观测值；

3)列向量索引集合Λ赋初值为空集；

4)列向量集合Φ赋初值为空集；

5)设循环标识k＝1，在循环中循环标识将逐次加1，直至满足循环终止条件；

步骤三：寻找传感矩阵Θ中与残余信号矢量r相关性最大的列对应的索引λ_k，此处的相关性通过传感矩阵Θ的列向量与残余信号矢量r的内积值来表征，计算Θ中所有列向量分别和r的内积，并记录最大内积值对应的Θ中的列向量，即数学表达式为其中φ_j表示矩阵Θ的第j列，argmax(f(x))表示使得f(x)取得最大值的的自变量x；

所述内积值的计算即计算Θ^Hr得到内积向量，再挑选内积向量中绝对值的最大值，具体实现方法为：

首先将残余信号矢量r按列重排为K×L的矩阵R，其中K为距离向压缩采样点数，L为方位向压缩采样点数，利用如下公式计算内积向量：

其中表示克罗内克积，得到传感矩阵Θ各个列向量和r的内积向量；

步骤四：利用克罗内克积的运算对应关系，即由索引关系式λ_k＝(q-1)×P+p，得出q和p的值，其中P和Q是成像场景划分为P×Q的二维网格后的网格数，其中，p＝1,2,3,...,P，q＝1,2,3,...,Q，计算可得Θ中对应的列向量更新索引集Λ＝Λ∪{λ_k}，将元素λ_k并入索引集合Λ，更新传感矩阵中的列集合Φ＝[Φ,φ_k]，将列向量φ_k并入列集合Φ中，即列向量集合Φ是由索引集合Λ所指示的Θ中的列向量构成的矩阵；

步骤五：利用最小二乘法估计在当前列集合Φ下的目标信号其中Φ^◇表示矩阵Φ的伪逆，计算公式为

步骤六：更新残余信号矢量令循环标识k加1，进入步骤七；

步骤七：判断是否满足k＞I或||r||₂＜ε，其中，ε为残余信号的能量阈值，若不满足上述条件之一，则跳转到步骤三，再次执行步骤三至步骤六；若至少满足以上两个条件之一则跳出循环，进入步骤八；

步骤八：按照建立的索引位置对赋值，其中，是重建得到的目标场景散射分布，它表示场景中不同网格点处的散射强度；Λ是步骤四中更新得到的索引集合，因此表示由索引集合Λ所指定的中的位置；是步骤五中计算得出的最小二乘解，然后将按列重排为大小P×Q的矩阵即为所求的目标散射率分布二维目标像。