CN102393517B - 基于密度加权的模拟退火稀疏合成孔径雷达成像方法 - Google Patents

Abstract

本发明公开了一种基于密度加权的模拟退火稀疏合成孔径雷达成像方法，主要解决现有技术的雷达回传数据传输率高和运算量大，均匀稀疏采样成像存在“鬼影”场景的问题。其实现过程为：采用密度加权获得中间采样孔径；利用模拟退火优化算法以主瓣和副瓣峰值比最大为代价函数，优化剩余两边对称采样孔径，得到方位稀疏采样位置；根据方位稀疏采样位置实现稀疏回波数据采样以降低数据量；回波数据方位向空缺位置补零后完成成像。本发明在付出较小成像质量下降的代价下实现降低前端回传数据传输率，减少后端距离向脉压的计算量，提高了成像的处理速度并消除均匀稀疏采样成像的“鬼影”场景，可用于机/星载雷达系统对地观测成像。

Description

基于密度加权的模拟退火稀疏合成孔径雷达成像方法

技术领域

本发明属于雷达成像技术领域，具体地说是一种稀疏合成孔径雷达的成像方法，该方法可用于机/星载雷达系统对地观测成像。

背景技术

合成孔径雷达SAR安置于机/星载平台，是一种主动式微波遥感设备，可以全天候、全天时、远距离对场景进行观测。根据奈奎斯特采样定理，方位向雷达发射重复频率必须大于方位向合成孔径带宽2倍才能保证频谱无混叠。由于机/星载合成孔径雷达往往平台速度快和天线孔径小，使得方位向带宽大，这就要求雷达重复频率高，导致方位向脉冲数多从而观测数据量大，不利于实时数据回传。

传统完全采样成像，需要采集并回传大量回波数据，导致成像计算量大，处理速度慢。稀疏采样可以降低回波采集数据量，从而减少成像计算量，提高处理速度。于是近些年出现了稀疏合成孔径雷达成像。

最简单的稀疏方式是均匀稀疏采样成像，由于这种成像方法不满足奈奎斯特采样定理，方位向匹配滤波后会出现成对回波即栅瓣，导致成像结果存在主场景的多次复制即鬼影。当雷达工作于条带模式时，“鬼影”场景会叠加到对应其成像中心的主场景上，造成成像模糊。

常规优化稀疏采样成像，通过设计好的代价函数，应用全局优化算法，如遗传算法GA、粒子群优化PSO算法、模拟退火SA算法等，优化出稀疏采样位置或保留优化位置的采样回波，能消除栅瓣和压低副瓣。但常规的优化稀疏采样成像方法存在运算量大、优化时间长、副瓣不够低等问题。

发明内容

本发明的目的在于克服上述已有技术的不足，提出一种基于密度加权的模拟退火的稀疏合成孔径雷达成像方法，以减小运算量和优化时间，进一步降低副瓣，提高成像质量。

为实现上述目的，本发明的技术方案包括如下步骤：

(1)采用密度加权获得中间采样孔径：

(1a)设定全孔径采样数为D、稀疏孔径采样总数为M、初始中间孔径采样数为N₀，初始化1×D维采样标记矢量x＝0，采样标记矢量x元素为1表示采样，元素为0表示不采样，D、M、N₀均取正整数，且满足N₀＜M＜D；

(1b)由全孔径采样数D，获得D×1维泰勒窗矢量h＝taylorwin(D)，对窗矢量的最大值归一化，得到D×1维归一化泰勒窗矢量h₀＝h/max(h)，其中taylorwin和max分别表示泰勒窗函数和求最大值函数；

(1c)给循环变量赋初值：i＝fix(D/2)-fix(N₀/2)，这里i表示中间孔径采样序号，fix表示朝零取整；

(1d)随机产生一个服从均匀分布[0，1]区间的数：r＝rand(1)，以归一化泰勒窗矢量中第i个元素为参考，判断[h₀]_i＞r是否成立，若成立则[x]_i＝1，[x]_D+1-i＝1，否则直接执行步骤(1e)，其中rand表示产生一服从均匀分布的随机数，[x]_i表示矢量x中的第i个元素；

(1e)令i＝i+1，判断i≤ceil(D/2)是否成立，若成立，则返回步骤(1d)继续循环，否则执行步骤(1f)，这里ceil表示朝正无穷大取整；

(1f)循环结束后，统计中间采样孔径实际所含采样数N＝sum(x)，sum表示求和函数；

(2)利用模拟退火智能优化算法以主瓣和副瓣峰值比最大为代价函数，优化除中间采样孔径外的剩余两边对称采样孔径，得到方位向稀疏采样位置：

(2a)设定初始温度T₀，模拟退火总的迭代次数Q；

(2b)在采样标记矢量x剩余左边孔径中，随机填(M-N)/2个数值1，根据对称性在右边孔径填(M-N)/2个数值1，并记录1的位置，得到扰动前的位置矢量I_p＝find(x＝1)，find表示寻找满足条件元素的位置；

(2c)设计代价函数为：f(I_p)＝FindLocalPeak2(P(I_p))，其中FindLocalPeak2是返回参数矢量的第二个极大值的函数，P(I_p)为合成阵列的归一化方向图，且 $P\left(I_p\right)=\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\exp (j2\mathrm{\π}{\left[I_p\right]}_m\left(\sin \mathrm{\θ}{-\sin \mathrm{\θ}}_0\right)/\mathrm{\λ}),$ 其中θ为波束指向，θ₀为导向方向，λ为波长，j为虚数单位；

(2d)根据设计好的代价函数利用模拟退火算法优化除中间采样孔径外的剩余两边对称采样孔径，得到方位向稀疏采样位置；

(3)雷达根据方位向稀疏采样位置发射并接收回波，实现回波数据的稀疏采样；

(4)对稀疏采样得到的稀疏回波数据在方位向空缺位置补零，然后利用RD算法完成成像处理。

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

1)与传统完全采样成像相比，本发明在付出较小成像质量下降的代价下实现前端数据采集量的显著减少，可降低实时回传数据传输率；由于方位稀疏采样，相比于传统完全采样处理，可减少后端距离向脉冲压缩的计算量，提高了成像的处理速度。

2)与均匀稀疏采样成像相比，本发明成像结果只有主场景，没有鬼影场景；

3)本发明由于采用的密度加权获得中间采样孔径，通过模拟退火对称优化剩余两边采样孔径，使优化运算量可降至常规模拟退火算法的(M-N)/M，同时优化出来的方向图副瓣更低。

附图说明

图1是本发明的流程图；

图2是本发明采用密度加权获得中间孔径的子流程图；

图3是本发明优化除中间采样孔径外的剩余两边对称采样孔径的子流程图；

图4是本发明进行成像处理的子流程图；

图5是本发明优化稀疏稀疏采样和现有疏稀疏采样的仿真方向图；

图6是本发明优化稀疏稀疏采样和现有疏稀疏采样的成像仿真结果图。

具体实施方式

参照图1，本发明的具体实现步骤如下：

步骤1，采用密度加权获得中间采样孔径；

参照图2，本步骤的具体实现如下：

1.1)设定全孔径采样数为D、稀疏孔径采样总数为M、初始中间孔径采样数为N₀，初始化1×D维采样标记矢量x＝0，采样标记矢量x元素为1表示采样，元素为0表示不采样，D、M、N₀均取正整数，且满足N₀＜M＜D；

1.2)由全孔径采样数D，获得D×1维泰勒窗矢量h＝taylorwin(D)，对窗矢量的最大值归一化，得到D×1维归一化泰勒窗矢量h₀＝h/max(h)，其中taylorwin和max分别表示泰勒窗函数和求最大值函数；

1.3)给循环变量赋初值：i＝fix(D/2)-fix(N₀/2)，这里i表示中间孔径采样序号，fix表示朝零取整；

1.4)随机产生一个服从均匀分布[0，1]区间的数：r＝rand(1)，以归一化泰勒窗矢量中第i个元素为参考，判断[h₀]_i＞r是否成立，若成立则[x]_i＝1，[x]_D+1-i＝1，否则直接执行步骤1.5)，其中rand表示产生一服从均匀分布的随机数，[x]_i表示矢量x中的第i个元素；

1.5)令i＝i+1，判断i≤ceil(D/2)是否成立，若成立，则返回步骤1.4)继续循环，否则执行步骤1.6)，这里ceil表示朝正无穷大取整；

1.6)循环结束后，统计中间采样孔径实际所含采样数N＝sum(x)，sum表示求和函数。

步骤2，利用模拟退火智能优化算法以主瓣和副瓣峰值比最大为代价函数，优化除中间采样孔径外的剩余两边对称采样孔径，得到方位向稀疏采样位置：

参照图3，本步骤的具体实现如下：

2.1)设定初始温度T₀，模拟退火总的迭代次数Q；

2.2)在采样标记矢量x剩余左边孔径中，随机填(M-N)/2个数值1，根据对称性在右边孔径填(M-N)/2个数值1，并记录1的位置，得到扰动前的位置矢量I_p＝find(x＝1)，find表示寻找满足条件元素的位置；

2.3)设计代价函数为：f(I_p)＝FindLocalPeak2(P(I_p))，其中FindLocalPeak2是返回参数矢量的第二个极大值的函数，P(I_p)为合成阵列的归一化方向图，且 $P\left(I_p\right)=\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\exp (j2\mathrm{\π}{\left[I_p\right]}_m\left(\sin \mathrm{\θ}{-\sin \mathrm{\θ}}_0\right)/\mathrm{\λ}),$ 其中θ为波束指向，θ₀为导向方向，λ为波长，j为虚数单位；

2.4)根据设计好的代价函数利用模拟退火算法优化除中间采样孔径外的剩余两边对称采样孔径，得到方位向稀疏采样位置：

2.4a)设外循环变量k＝1，k表示第k次模拟退火迭代；

2.4b)设第k次的退火温度值为T_k＝T_k-1/ln(k+1)；

2.4c)初始化内循环变量n＝2，并令[I_c]^k＝[I_p]^k，其中n表示剩余孔径采样序号，I_c表示扰动后的位置矢量，[I_p]^k表示第k次迭代扰动前的采样位置矢量，[I_c]^k表示第k次迭代扰动后的采样位置矢量；

2.4d)将第n个采样位置在第n-1和n+1采样位置之间进行扰动，该扰动表示为 ${\left[I_c\right]}_n^k={\left[I_c\right]}_{n-1}^k+\mathrm{round}[({\left[I_p\right]}_{n+1}^k-{\left[I_c\right]]}_{n-1}^k)\×\mathrm{rand}\left(1\right)],$ 其中round表示四舍五入，

表示第k次迭代第n个序号扰动前的采样位置矢量，

表示第k次迭代第n个序号扰动后的采样位置矢量；

2.4e)为满足对称性，令

其中D为全孔径采样数，M为稀疏孔径采样总数；

2.4f)计算扰动前后代价函数的变化，即ΔE＝f([I_c]^k)-f([I_p]^k)：

2.4g)判断ΔE＜0是否成立，若不成立，则执行步骤2.4h)，否则跳转到步骤2.4i)；

2.4h)随机产生一个服从均匀分布[0，1]区间的数r，r＝rand(1)，然后判断exp(-ΔE/K₀T_k)＞r是否满足，其中K₀为玻尔兹曼常数，若满足，则执行步骤2.4i)，否则跳转到步骤2.4j)；

2.4i)接受扰动后的状态，即[I_p]^k＝[I_c]^k；

2.4j)更新内循环变量，n＝n+1；并判断n≤(M-N)/2是否成立，若不成立，则跳转到步骤2.4b)，否则执行步骤2.4k)；

2.4k)更新外循环变量：k＝k+1，并判断k≤Q是否成立，若成立则返回到步骤2.4b)，否则迭代结束，迭代结束后，I_p即为优化得到的方位向稀疏采样位置矢量。

步骤3，雷达根据方位向稀疏采样位置发射并接收回波，实现回波数据的稀疏采样。

步骤4，对稀疏采样得到的稀疏回波数据在方位向空缺位置补零，然后利用RD算法完成成像处理。

参照图4，本步骤所述的利用RD算法完成成像处理，其具体实现如下：

4.1)利用相关函数法，对补零后的稀疏回波数据的多普勒中心频率进行估计；

4.2)根据估计出来的多普勒中心频率，进行距离徙动校正；

4.3)距离徙动校正后，对方位稀疏采样位置对应的距离向数据进行频域脉冲压缩处理；

4.4)距离脉冲压缩后，利用图像偏置法进行多普勒调频率估计；

4.5)根据估计出来的多普勒中心频率和多普勒调频率，进行方位频域脉冲压缩处理后，得到SAR图像。

本发明的优点可以通过以下仿真实验进一步说明。

仿真实验一：稀疏采样优化仿真

1a.仿真条件

以稀疏度M/D＝75％为例，设置稀疏采样优化仿真参数如表1所示：

表1  本发明稀疏采样优化仿真参数表

迭代次数Q	2000	初始温度T0(K)	2000
				全孔径采样数D	513	稀疏孔径采样总数M	385
初始中间孔径采样数N0	285

1b.仿真实例

根据表1所设参数，按本发明所述步骤进行稀疏采样优化仿真，得到本发明优化稀疏稀疏采样的采样位置，并根据优化出来的采样位置和现有疏稀疏采样方式得到的采样位置，画出本发明优化稀疏稀疏采样和现有疏稀疏采样的仿真方向图，如图(5)所示。其中图5(a)中的点划线为采用传统完全采样方式得到的方向图，图5(a)中实线为采用均匀稀疏采样方式得到的方向图，图5(b)中的点划线为采用常规模拟退火优化稀疏采样方式得到的方向图，图5(b)中的实线为采用本发明优化稀疏采样方式得到的方向图。

由图5(a)可知，传统完全采样方向图的最高副瓣为-13.4dB，均匀稀疏采样方向图有-4dB的栅瓣；由图5(b)可知，常规模拟退火优化稀疏采样方向图的最高副瓣为-15.96dB，本发明优化稀疏采样方向图的最高副瓣为-20dB。因此，采用本发明优化稀疏采样具有最好的副瓣性能，有利于后续的稀疏成像。

仿真实验二：成像仿真

2a.仿真条件

以正侧视条带式合成孔径雷达成像为例，仿真参数如表2所示：

表2  正侧视条带式雷达仿真参数表

波长(m)	0.03	场景中心斜距(Km)	30
				载频(GHz)	10	天线长度(m)	1.2
平台高度(Km)	4	发射信号带宽(MHz)	75
				平台速度(m/s)	250	发射信号时宽(us)	10
采样频率(MHz)	150	脉冲重复频率(Hz)	500

2b.仿真实例

根据仿真一得到的四种方位向采样位置及表2参数，在提取实际雷达成像场景反射系数后，高保真的生成雷达回波数据，并分别利用RD算法完成成像处理后得到的SAR图像如图6所示。其中图6(a)为采用传统完全采样方式的成像结果，图6(b)为采用均匀稀疏采样方式的成像结果，图6(c)为采用常规模拟退火优化稀疏采样方式的成像结果，图6(d)为采用本发明优化稀疏采样方式的成像结果。

由图6(b)可知，均匀稀疏采样成像存在“鬼影”场景，这与图5(a)中该采样的栅瓣相对应，由于是条带式成像，“鬼影”场景会叠加到对应其成像中心的主场景上，造成成像模糊；对比图6(c)和图6(d)可知，与常规模拟退火优化稀疏采样成像结果相比，本发明成像结果主场景成像质量更清晰，副瓣场景不明显，这与图5(b)中本发明副瓣更低有关；对比图6(a)和图6(d)可知，与完全采样成像结果相比，本发明主场景成像质量的下降不明显，是可以接受的。

Claims (3)

1.一种基于密度加权的模拟退火稀疏合成孔径雷达成像方法，包括如下步骤：

(1)采用密度加权获得中间采样孔径：

(1a)设定全孔径采样数为D、稀疏孔径采样总数为M、初始中间孔径采样数为N₀，初始化1×D维采样标记矢量x＝0，采样标记矢量x元素为1表示采样，元素为0表示不采样，D、M、N₀均取正整数，且满足N₀＜M＜D；

(1b)由全孔径采样数D，获得D×1维泰勒窗矢量h＝taylorwin(D)，对窗矢量的最大值归一化，得到D×1维归一化泰勒窗矢量h₀=h/max(h)，其中taylorwin和max分别表示泰勒窗函数和求最大值函数；

(1c)给循环变量赋初值：i=fix(D/2)-fix(N₀/2)，这里i表示中间孔径采样序号，fix表示朝零取整；

(1d)随机产生一个服从均匀分布[0,1]区间的数：r＝rand(1)，以归一化泰勒窗矢量中第i个元素为参考，判断[h₀]_i>r是否成立，若成立则[x]_i＝1，[x]_D+1-i＝1，否则直接执行步骤(1e)，其中rand表示产生一服从均匀分布的随机数，[x]_i表示矢量x中的第i个元素；

(1e)令i＝i+1,判断i≤ceil(D/2)是否成立，若成立，则返回步骤(1d)继续循环，否则执行步骤(1f)，这里ceil表示朝正无穷大取整；

(1f)循环结束后，统计中间采样孔径实际所含采样数N＝sum(x)，sum表示求和函数；

(2)利用模拟退火智能优化算法以主瓣和副瓣峰值比最大为代价函数，优化除中间采样孔径外的剩余两边对称采样孔径，得到方位向稀疏采样位置：

(2a)设定初始温度T₀，模拟退火总的迭代次数Q；

(2b)在采样标记矢量x剩余左边孔径中，随机填(M-N)/2个数值1，根据对称性在右边孔径填(M-N)/2个数值1，并记录1的位置，得到扰动前的位置矢量I_p＝find(x＝1)，find表示寻找满足条件元素的位置；

(2c)设计代价函数为：f(I_p)＝FindLocalPeak2(P(I_p))，其中FindLocalPeak2是返回参数矢量的第二个极大值的函数，P(I_p)为合成阵列的归一化方向图，且 $P\left(I_p\right)=\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\exp (j2\mathrm{\π}{\left[I_p\right]}_m(\sin \mathrm{\θ}-{\sin \mathrm{\θ}}_0)/\mathrm{\λ}),$ 其中θ为波束指向，θ₀为导向方向，λ为波长，j为虚数单位；

(2d)根据设计好的代价函数利用模拟退火算法优化除中间采样孔径外的剩余两边对称采样孔径，得到方位向稀疏采样位置；

(3)雷达根据方位向稀疏采样位置发射并接收回波，实现回波数据的稀疏采样；

(4)对稀疏采样得到的稀疏回波数据在方位向空缺位置补零，然后利用RD算法完成成像处理。

2.根据权利要求1所述的方法，其特征在于步骤(2d)，按如下步骤进行：

2d1)设外循环变量k＝1，k表示第k次模拟退火迭代；

2d2)设第k次的退火温度值为T_k＝T_k-1/ln(k+1)；

2d3)初始化内循环变量n＝2，并令[I_c]^k＝[I_p]^k，其中n表示剩余孔径采样序号，I_c表示扰动后的位置矢量，[I_p]^k表示第k次迭代扰动前的采样位置矢量，[I_c]^k表示第k次迭代扰动后的采样位置矢量；

2d4)将第n个采样位置在第n-1和n+1采样位置之间进行扰动，该扰动表示为 ${\left[I_c\right]}_n^k={\left[I_c\right]}_{n-1}^k+\mathrm{round}[({\left[I_p\right]}_{n+1}^k-{\left[I_c\right]}_{n-1}^k)\×\mathrm{rand}\left(1\right)],$ 其中round表示四舍五入，

表示第k次迭代第n个序号扰动前的采样位置矢量，

表示第k次迭代第n个序号扰动后的采样位置矢量；

2d5)为满足对称性，令

其中D为全孔径采样数，M为稀疏孔径采样总数；

2d6)计算扰动前后代价函数的变化，即ΔE＝f([I_c]^k)-f([I_p]^k)；

2d7)判断ΔE＜0是否成立，若不成立，则执行步骤2d8)，否则跳转到步骤2d9)；

2d8)随机产生一个服从均匀分布[0,1]区间的数r，r＝rand(1)，然后判断exp(-ΔE/K₀T_k)>r是否满足，其中K₀为玻尔兹曼常数，若满足，则执行步骤2d9)，否则跳转到步骤2d10)；

2d9)接受扰动后的状态，即[I_p]^k＝[I_c]^k；

2d10)更新内循环变量，n＝n+1；并判断n≤(M-N)/2是否成立，若不成立，则跳转到步骤2d2)，否则执行步骤2d11)；

2d11)更新外循环变量：k＝k+1，并判断k≤Q是否成立，若成立则返回到步骤2d2)，否则迭代结束，迭代结束后，I_p即为优化得到的方位向稀疏采样位置矢量。

3.根据权利要求1所述的方法，其特征在于步骤(4)所述的利用RD算法完成成像处理，按如下步骤进行：

1)利用相关函数法，对补零后的稀疏回波数据的多普勒中心频率进行估计；

2)根据估计出来的多普勒中心频率，进行距离徙动校正；

3)距离徙动校正后，对方位向稀疏采样位置对应的距离向数据进行频域脉冲压缩处理；

4)距离脉冲压缩后，利用图像偏置法进行多普勒调频率估计；

5)根据估计出来的多普勒中心频率和多普勒调频率，进行方位频域脉冲压缩处理后，得到SAR图像。

Images