CN105891826B - 机载雷达快速最大后验成像方法 - Google Patents

Abstract

本发明公开了一种机载雷达快速最大后验成像方法，包括以下步骤：S1、对二维回波数据进行距离向脉冲压缩；S2、获取天线方向图信息，构造卷积矩阵；S3、导出乘性迭代求解公式；S4、初始化迭代；S5、初始化向量外推预测；S6、向量外推预测迭代，利用相邻的前三次的迭代结果进行进一步向量预测，并将预测结果带入乘性迭代求解公式得到下一次迭代结果；S7、判断相邻两次迭代结果是否满足迭代终止条件，若满足则进行步骤S8；否则返回步骤S6；S8、对回波的每行数据按照步骤S4～S7进行处理，处理完成后输出成像结果。本发明在传统迭代步骤之前利用向量外推方法增加一步预测来估计目标，可以在保证锐化效果的前提下，有效减少迭代所需次数，提高运算效率。

Description

机载雷达快速最大后验成像方法

技术领域

本发明属于雷达成像技术领域，特别涉及一种机载雷达快速最大后验成像方法。

背景技术

雷达成像技术凭借其全天时、全天候和高成像质量的优点，已成为当今探测领域不可取代的技术手段，在民用、军事领域都发挥着越来越重要的作用。由于常规的单站SAR和DBS技术无法实现前视成像，实波束机载扫描雷达波束锐化技术成为实现前视成像的一种重要技术手段。

文献“Iterative noncoherent angular superresolution”(Richards M.RadarConference,1988.,Proceedings of the 1988IEEE National.IEEE,1988:100-105.)采用了一种非线性限制迭代逆滤波技术对实波束进行锐化，但此方法没有利用噪声和目标的先验信息，在低信噪比情况下性能较差。文献“Superresolution of coherent sources inreal-beam data”(Uttam S,Goodman N.Aerospace and Electronic Systems,IEEETransactions on,2010,46(3):1557-1566.)利用扫描回波的二阶统计特性对回波进行锐化，但此方法需要大量快拍数，无法应用在在单次或快拍数较少的情况。文献“BayesianDeconvolution for Angular Super-Resolution in Forward-Looking Scanning Radar”(Zha Y,Huang Y,Sun Z,et al..Sensors,2015,15(3):6924-6946)在贝叶斯框架下导出最大后验算法，利用噪声与目标的先验信息对目标散射系数进行估计，取得了鲁棒性较强的波束锐化结果，但由于采用的乘性迭代策略收敛速度慢限制了其在实时处理中的应用。

发明内容

本发明的目的是针对背景技术存在的缺陷，提供一种利用最大后验方法相邻迭代向量方向的相似性，在传统迭代步骤之前利用向量外推方法增加一步预测来估计目标，解决了最大后验迭代方法里存在的收敛速度慢和无法实时成像问题的机载雷达快速最大后验成像方法。

本发明的目的是通过以下技术方案来实现的：机载雷达快速最大后验成像方法，包括以下步骤：

S1、对二维回波数据进行距离向脉冲压缩；

S2、获取天线方向图信息，根据系统参数对天线方向图离散化，得到离散化的天线方向图，由离散化的天线方向图构造卷积矩阵；

S3、根据步骤S2得到的卷积矩阵导出乘性迭代求解公式；

S4、初始化迭代，对二维回波信号的一行数据取幅值，并将其设置为迭代初始值，利用乘性迭代求解公式得到第一次迭代结果，再利用第一次迭代结果和乘性迭代求解公式得到第二次迭代结果；

S5、初始化向量外推预测，利用第二次迭代结果、第一次迭代结果和迭代初始值进行第一步向量预测，并将得到的预测结果带入乘性迭代求解公式得到第三次迭代结果；

S6、向量外推预测迭代，利用相邻的前三次的迭代结果进行进一步向量预测，并将得到的预测结果带入乘性迭代求解公式得到下一次迭代结果；

S7、判断相邻两次迭代结果是否满足迭代终止条件，若满足则进行步骤S8；否则返回步骤S6，继续进行迭代计算直至满足迭代终止条件；

S8、输出成像结果，对回波的每行数据按照步骤S4～S7进行处理，处理完成后输出成像结果。

进一步地，步骤S1具体实现方法为：针对获取的二维回波数据Echo(τ,t)，利用雷达发射信号参数构造距离向脉冲压缩参考函数refer(τ)＝exp(iπk_rτ²)，其中k_r＝10MHz/μs，τ＝-1μs～1μs；将二维回波数据Echo(τ,t)沿距离向做傅里叶变换，将傅里叶变换得到的结果与距离向脉冲压缩参考函数的傅里叶变换的共轭refer^*(f)相乘，再将相乘后结果变换到时域，完成对二维回波数据Echo(τ,t)的距离向脉冲压缩；记压缩后的回波数据为s(τ,t)，其大小为M×L，其中M为回波数据的距离采样点数，L为方位采样点数。

进一步地，步骤S2具体实现方法为：获取天线方向图信息，根据系统参数对天线方向图离散化，得到离散化的天线方向图h＝[h₀h₁...h_m-1]，其长度由此天线方向图构造卷积矩阵H：

其中，θ为天线方向图波束宽度，ω为扫描速度，prf为脉冲重复频率。

进一步地，步骤S3具体实现方法为：根据以下成像模型导出乘性迭代求解公式：

s＝Hx+n (1)

其中，x为待求解的目标函数，s为方位回波，n为噪声；

将噪声和目标函数分别建模为泊松分布和拉普拉斯分布，得到方位回波的后验概率密度：

其中，σ为分布参数，则目标函数为此后验概率密度的自然对数：

将上式中与变量x无关的项略去，得到新的目标函数为：

||x||₁为向量x的一范数；对此目标函数求导，得到其梯度为：

其中ε是为解决一范数的不可导性而引入的一个很小的常量；diag{·}为对角矩阵；

最后结合目标函数导数由固定点求解策略得到乘性迭代公式：

其中，x_k和x_k+1为相邻两次迭代结果。

进一步地，步骤S4具体实现方法为：对二维回波信号s(τ,t)的一行数据取幅值，即为待处理的原始方位回波s，并将其设置为迭代初始值x₀；利用公式(6)，得到迭代结果x₁；再利用x₁和迭代公式(6)得到第二次迭代结果x₂。

进一步地，步骤S5具体实现方法为：利用x₂、x₁、x₀，根据式(7)进行第一步向量预测：

y₁＝x₂+α₂h₂ (7)

其中，

h₂＝x₂-x₁ (8)

g₂＝x₂-x₁ (10)

g₁＝x₁-x₀ (11)

得到预测结果y₁，并将y₁赋值给x_k代入迭代公式(6)，得到第三次迭代结果x₃，并将得到的x₃、x₂、y₁、g₂分别赋值给x_k、x_k-1，y_k，g_k-1，为下一次迭代做准备。

进一步地，步骤S6具体实现方法为：将x_k、x_k-1，y_k，g_k-1，代入式(12)进行一步预测

y_k+1＝x_k+α_kh_k (12)

其中，

h_k＝x_k-x_-k1 (13)

g_k＝x_k-y_k (15)

得到预测结果y_k+1，并将其赋值给x_k，代入公式(6)，得到迭代结果x_k+1。

进一步地，步骤S7的具体实现方法为：判断相邻两次迭代结果是否满足下列迭代终止条件：

||x_k+1-x_k||²＜δ (16)

其中，δ为用户设定的阈值；若相邻两次迭代结果x_k+1、x_k满足式(16)，则进入步骤S8；若不满足则将x_k+1、x_k、y_k+1、g_k分别赋值给x_k、x_k-1、y_k、和g_k-1，返回步骤S6，直至满足迭代终止条件。

本发明的有益效果是：本发明采用的快速迭代求解策略利用了最大后验方法相邻迭代向量方向的相似性，在传统迭代步骤之前利用向量外推方法增加一步预测来估计目标，可以在保证锐化效果的前提下，有效减少迭代所需次数，由于预测步骤计算量小，不会显著增加每次迭代所需时间，因此能大幅度减少算法的最终运算时间，提高运算效率，解决了最大后验算法无法实时成像的问题。

附图说明

图1是本发明具体实施例中采用的实波束前视扫描雷达(RBSR)成像示意图；

图2是本发明具体实施例中采用的目标场景布置图；图中的黑色圆点为布置于地面上的3个点目标，这3个点沿方位方向间隔0.2°；

图3是本发明方法的流程框图；

图4是本发明具体实施例中距离向脉冲压缩后的数据；

图5是本发明具体实施例中距离向脉冲压缩后目标所在距离单元的剖面图；

图6是本发明具体实施例中迭代1125次后结果；

图7是图6结果中目标所在距离单元的剖面图。

具体实施方式

本发明主要采用仿真实验的方法进行验证，所以步骤、结论都在Matlab2012上验证正确。下面结合附图和具体实施例对本发明方法做进一步的阐述。

图1是本实施例所采用的实波束前视扫描雷达(RBSR)成像示意图，雷达边扫描边发射线性调频信号，其系统参数与扫描参数如表一所示，则方位采样点数L＝334。图2展示了仿真目标场景，其主要目标为斜距3km处相邻的三个等幅度的点，方位间距为0.2°。以采样率f_r＝100MHz对距离向回波进行采样，得到回波矩阵，记为Echo(τ,t)，其大小为440×334。

表一雷达系统参数与扫描参数

参数	符号	数值
			载频	f0	30.75GHz
发射信号带宽	Br	40MHz
			发射信号时宽	Tr	2μs
发射信号调频率	kr	10MHz/μs
			距离采样率	fr	100MHz
脉冲采样频率	prf	2000Hz
			天线方向图长度	θ	1.95°
扫描速度	ω	600/s
			扫描范围	θmin～θmax	-50～50

如图3所示，本发明的机载雷达快速最大后验成像方法，包括以下步骤：

S1、对二维回波数据进行距离向脉冲压缩；其具体实现方法为：针对获取的二维回波数据Echo(τ,t)，利用雷达发射信号参数构造距离向脉冲压缩参考函数refer(τ)＝exp(iπk_rτ²)，其中k_r＝10MHz/μs，τ＝-1μs～1μs；将二维回波数据Echo(τ,t)沿距离向做傅里叶变换，将傅里叶变换得到的结果与距离向脉冲压缩参考函数的傅里叶变换的共轭refer^*(f)相乘，再将相乘后结果变换到时域，完成对二维回波数据Echo(τ,t)的距离向脉冲压缩；记压缩后的回波数据为s(τ,t)，其大小保持440×334不变，如图4所示。

S2、获取天线方向图信息，根据天线方向图宽度θ、脉冲重复频率prf和扫描速度ω对扫描天线方向图离散化，由可得到离散化的天线方向图点数m＝66，则h＝[h₀h₁...h₆₅]，由此天线方向图构造卷积矩阵H：

S3、根据步骤S2得到的卷积矩阵导出乘性迭代求解公式；具体实现方法为：根据以下成像模型导出乘性迭代求解公式：

s＝Hx+n (1)

其中，x为待求解的目标函数，s为方位回波，n为噪声；

将噪声和目标函数分别建模为泊松分布和拉普拉斯分布，得到方位回波的后验概率密度：

其中，σ为分布参数，则目标函数为此后验概率密度的自然对数：

将上式中与变量x无关的项略去，得到新的目标函数为：

||x||₁为向量x的一范数；对此目标函数求导，得到其梯度为：

其中，ε是为解决一范数的不可导性而引入的一个很小的常量，diag{·}为对角矩阵；本实施例中σ取9.8，ε取10e(-8)。

最后结合目标函数导数由固定点求解策略得到乘性迭代公式：

其中，x_k和x_k+1为相邻两次迭代结果。

S4、初始化迭代，对二维回波信号的一行数据取幅值，并将其设置为迭代初始值，利用乘性迭代求解公式得到第一次迭代结果，再利用第一次迭代结果和乘性迭代求解公式得到第二次迭代结果；具体实现方法为：对二维回波信号s(τ,t)的一行数据取幅值，即为待处理的原始方位回波s，并将其设置为迭代初始值x₀；利用公式(6)，得到迭代结果x₁；再利用x₁和迭代公式(6)得到第二次迭代结果x₂。

S5、初始化向量外推预测，利用第二次迭代结果、第一次迭代结果和迭代初始值进行第一步向量预测，并将得到的预测结果带入乘性迭代求解公式得到第三次迭代结果；具体实现方法为：利用x₂、x₁、x₀，根据式(7)进行第一步向量预测：

y₁＝x₂+α₂h₂ (7)

其中，

h₂＝x₂-x₁ (8)

g₂＝x₂-x₁ (10)

g₁＝x₁-x₀ (11)

得到预测结果y₁，并将y₁赋值给x_k代入迭代公式(6)，得到第三次迭代结果x₃，并将得到的x₃、x₂、y₁、g₂分别赋值给x_k、x_k-1，y_k，g_k-1，为下一次迭代做准备。

S6、向量外推预测迭代，利用相邻的前三次的迭代结果进行进一步向量预测，并将得到的预测结果带入乘性迭代求解公式得到下一次迭代结果；具体实现方法为：将x_k、x_k-1，y_k，g_k-1，代入式(12)进行一步预测

y_k+1＝x_k+α_kh_k (12)

其中，

h_k＝x_k-x_-k1 (13)

g_k＝x_k-y_k (15)

得到预测结果y_k+1，并将其赋值给x_k，代入公式(6)，得到迭代结果x_k+1。

S7、判断相邻两次迭代结果是否满足迭代终止条件，若满足则进行步骤S8；否则返回步骤S6，继续进行迭代计算直至满足迭代终止条件；具体实现方法为：判断相邻两次迭代结果是否满足下列迭代终止条件：

||x_k+1-x_k||²＜δ (16)

其中，δ为用户设定的阈值，本实施例δ取为0.001；若相邻两次迭代结果x_k+1、x_k满足式(16)，则进入步骤S8；若不满足则将x_k+1、x_k、y_k+1、g_k分别赋值给x_k、x_k-1、y_k、和g_k-1，返回步骤S6，直至满足迭代终止条件。

S8、输出成像结果，对回波的440行数据按照步骤S4～S7进行处理，处理完成后输出成像结果。

目标回波经脉冲压缩后的仿真结果s(τ,t)及目标点所在行的剖面图如图4、图5所示，由于三个点目标方位间隔小于天线波束宽度，故在实波束回波中无法分别三点，图6、图7为本发明的快速最大后验方法迭代1125次后结果，从结果可以看出三个目标点被明显分开，且幅度差异不大；表二为背景技术最大后验方法锐化结果的均方值误差达到图7所示结果的均方值误差水平时所需迭代次数和计算时间。从表中可以看出要得到快速迭代方法1125次时的锐化结果，原始最大后验方法需要迭代14816次，即收敛速度提高了近13.2倍；通过对比两种方法所需时间可知运算效率提高了近12.6倍，说明本发明的算法相较于原始方法确实可以更高效的对实波束进行锐化。

表二迭代次数与计算时间表

	最大后验	加速最大后验
			迭代次数	14816	1125
计算时间	5.15(s)	0.41(s)

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (8)

1.机载雷达快速最大后验成像方法，其特征在于，包括以下步骤：

S1、对二维回波数据进行距离向脉冲压缩；

S2、获取天线方向图信息，根据系统参数对天线方向图离散化，得到离散化的天线方向图，由离散化的天线方向图构造卷积矩阵；

S3、根据步骤S2得到的卷积矩阵导出乘性迭代求解公式；

S4、初始化迭代，对二维回波信号的一行数据取幅值，并将其设置为迭代初始值，利用乘性迭代求解公式得到第一次迭代结果，再利用第一次迭代结果和乘性迭代求解公式得到第二次迭代结果；

S5、初始化向量外推预测，利用第二次迭代结果、第一次迭代结果和迭代初始值进行第一步向量预测，并将得到的预测结果带入乘性迭代求解公式得到第三次迭代结果；

S6、向量外推预测迭代，利用相邻的前三次的迭代结果进行进一步向量预测，并将得到的预测结果带入乘性迭代求解公式得到下一次迭代结果；

S7、判断相邻两次迭代结果是否满足迭代终止条件，若满足则进行步骤S8；否则返回步骤S6，继续进行迭代计算直至满足迭代终止条件；

S8、输出成像结果，对回波的每行数据按照步骤S4～S7进行处理，处理完成后输出成像结果。

2.根据权利要求1所述的机载雷达快速最大后验成像方法，其特征在于，所述的步骤S1具体实现方法为：针对获取的二维回波数据Echo(τ,t)，利用雷达发射信号参数构造距离向脉冲压缩参考函数refer(τ)＝exp(iπk_rτ²)，其中k_r＝10MHz/μs，τ＝-1μs～1μs；将二维回波数据Echo(τ,t)沿距离向做傅里叶变换，将傅里叶变换得到的结果与距离向脉冲压缩参考函数的傅里叶变换的共轭refer^*(f)相乘，再将相乘后结果变换到时域，完成对二维回波数据Echo(τ,t)的距离向脉冲压缩；记压缩后的回波数据为s(τ,t)，其大小为M×L，其中M为回波数据的距离采样点数，L为方位采样点数。

3.根据权利要求2所述的机载雷达快速最大后验成像方法，其特征在于，所述的步骤S2具体实现方法为：获取天线方向图信息，根据系统参数对天线方向图离散化，得到离散化的天线方向图h＝[h₀h₁...h_m-1]，其长度由此天线方向图构造卷积矩阵H：

其中，θ为天线方向图波束宽度，ω为扫描速度，prf为脉冲重复频率。

4.根据权利要求3所述的机载雷达快速最大后验成像方法，其特征在于，所述的步骤S3具体实现方法为：根据以下成像模型导出乘性迭代求解公式：

s＝Hx+n (1)

其中，x为待求解的目标函数，s为方位回波，n为噪声；

将噪声和目标函数分别建模为泊松分布和拉普拉斯分布，得到方位回波的后验概率密度：

其中，σ为分布参数，则目标函数为此后验概率密度的自然对数：

将上式中与变量x无关的项略去，得到新的目标函数为：

x^T为x的转置；

||x||₁为向量x的一范数；对此目标函数求导，得到其梯度为：

其中ε是为解决一范数的不可导性而引入的一个很小的常量；diag{·}为对角矩阵；

最后结合目标函数导数由固定点求解策略得到乘性迭代公式：

其中，x_k和x_k+1为相邻两次迭代结果。

5.根据权利要求4所述的机载雷达快速最大后验成像方法，其特征在于，所述的步骤S4具体实现方法为：对二维回波信号s(τ,t)的一行数据取幅值，即为待处理的原始方位回波s，并将其设置为迭代初始值x₀；利用公式(6)，得到迭代结果x₁；再利用x₁和迭代公式(6)得到第二次迭代结果x₂。

6.根据权利要求5所述的机载雷达快速最大后验成像方法，其特征在于，所述的步骤S5具体实现方法为：利用x₂、x₁、x₀，根据式(7)进行第一步向量预测：

y₁＝x₂+α₂h₂ (7)

其中，

h₂＝x₂-x₁ (8)

g₂＝x₂-x₁ (10)

g₁＝x₁-x₀ (11)

得到预测结果y₁，并将y₁赋值给x_k代入迭代公式(6)，得到第三次迭代结果x₃，并将得到的x₃、x₂、y₁、g₂分别赋值给x_k、x_k-1，y_k，g_k-1，为下一次迭代做准备。

7.根据权利要求6所述的机载雷达快速最大后验成像方法，其特征在于，所述的步骤S6具体实现方法为：将x_k、x_k-1，y_k，g_k-1，代入式(12)进行一步预测

y_k+1＝x_k+α_kh_k (12)

其中，

h_k＝x_k-x_-k1 (13)

g_k＝x_k-y_k (15)

得到预测结果y_k+1，并将其赋值给x_k，代入公式(6)，得到迭代结果x_k+1。

8.根据权利要求7所述的机载雷达快速最大后验成像方法，其特征在于，所述的步骤S7的具体实现方法为：判断相邻两次迭代结果是否满足下列迭代终止条件：

||x_k+1-x_k||²＜δ (16)

其中，δ为用户设定的阈值；若相邻两次迭代结果x_k+1、x_k满足式(16)，则进入步骤S8；若不满足则将x_k+1、x_k、y_k+1、g_k分别赋值给x_k、x_k-1、y_k、和g_k-1，返回步骤S6，直至满足迭代终止条件。