CN104849713B - 一种基于slim算法的sar成像实现方法 - Google Patents

Abstract

本发明公开了一种基于SLIM算法的SAR成像实现方法，包括步骤一：读入相关参数和回波数据，步骤二：根据SAR回波模型，计算出场景中每个点在每个脉冲发射时刻的斜距、方位向包络和每个点的零多普勒时刻，以及每个点在每个脉冲发射时刻、每个距离向采样时刻的距离向包络，步骤三：建立SAR回波数据估计矩阵A，并将回波数据列向量化，步骤四：根据估计矩阵A，运用SLIM估计算法，对列向量化后的回波数据进行迭代运算直至收敛，输出场景目标RCS估计结果。本发明能够获得地面目标后向散射系数的精确值，同时有效抑制旁瓣能量，提升SAR图像质量，便于SAR图像的判读及后续的进一步应用处理。

Description

一种基于SLIM算法的SAR成像实现方法

技术领域

本发明属于图像处理领域，涉及一种SAR成像实现方法，特别涉及运用稀疏学习迭代最小化(SLIM)算法的SAR成像实现方法。

背景技术

合成孔径雷达(SAR)卫星由于不受天气、地理、时间等因素的限制，能够实现全天时、全天候的对地观测，且具有一定的穿透力，因而被广泛的应用于军事侦察、地形测绘、资源探测、海洋观测、生态监测、自然灾害监测、快速救援等方面。

和光学雷达不同，SAR系统是通过发射线性调频信号，并接收地面目标反射的回波信号，来完成对地观测。通过后续的成像处理获得对应的雷达图像。因此，与光学雷达相比，SAR系统多一个处理步骤，需要将二维回波数据通过一定的处理算法转换成地面目标的后向散射系数图。

随着SAR技术的发展，各种成像方法被不断提出。目前常用的成像算法主要依托于匹配滤波处理来展开，典型的成像处理算法包括距离-多普勒(R-D)算法、Chirp Scaling(CS)算法、ωK算法等等。

1、距离-多普勒(R-D)算法

距离-多普勒(R-D)算法通过距离和方位上的频域操作，达到了高效的模块化处理要求，同时又具有一维操作的简便性。该算法根据距离和方位上的大尺度时间差异，在两个一维操作之间使用距离徙动校正(RCMC)，对距离和方位进行了近似的分离处理。

2、Chirp Scaling(CS)算法

CS成像处理算法是一种精确的成像算法，它通过频域相位补偿完成成像处理，是一种具有良好相位保持特性的成像处理方法。CS算法从原始回波信号入手，精确推导回波信号在距离-多普勒域的表达式，利用CS原理来完成不同距离门距离徙动曲线的修正，通过二维频域内的相位补偿来完成距离徙动校正、距离压缩处理及二次距离脉压，最后通过在方位多普勒域内补偿方位相位和方位向傅立叶逆变换来获得最终的成像处理结果

3、ωK算法

ωK算法是一种将电磁波传播方式等效成地震波传播方式而导出的一种合成孔径雷达成像算法。在算法的推导过程中没有进行任何近似，是一种优秀的成像算法，具有很高的成像处理精度。由于处理算法在进行变量替换处理过程中需要采用插值处理，这不仅增加了成像处理的运算量，也影响成像处理的精度

然而不管是距离-多普勒(R-D)算法、Chirp Scaling算法还是ωK算法，均是依托于匹配滤波器来完成图像聚焦。虽然能从回波信号中反演出场景目标的散射系数，但是具有较高的旁瓣。过高的旁瓣能量会淹没图像中的弱目标，降低图像质量，进而影响图像判读。

发明内容

本发明提出了一种基于SLIM算法的SAR成像实现方法，该成像方法以SAR回波模型为基础，根据SAR雷达与地面场景目标的空间几何关系，建立SAR回波数据估计矩阵，并通过SLIM算法进行迭代计算，完成对SAR回波数据的成像处理，获得地面目标后向散射系数的精细估计结果，并有效抑制旁瓣能量。

一种基于SLIM算法的SAR成像实现方法，包括以下几个步骤：

步骤一：读入相关参数和回波数据，具体步骤又分为：

A、读入成像参数，包括斜距参数、多普勒参数和雷达信号参数：参考斜距R_ref、各个距离门y的多普勒中心频率f_d,y、各个距离门y的多普勒调频率f_r,y、距离向调频率K_r、波长λ、信号采样率f_s、脉冲重复频率f_prf、脉冲持续时间T_r、雷达速度v、雷达照射中心时刻t₀、天线方位向长度L_a；

B、读入SAR回波信号数据Echo和SAR回波大小X×Y(表示Echo为一个X行Y列的数据矩阵)；

步骤二：根据SAR回波模型及相关参数，计算出场景中每个点(x,y)在每个脉冲发射时刻η_a的斜距R_x,y(η_a)、方位向包络w_a,x,y(η_a)和每个点(x,y)的零多普勒时刻η_c,x,y，以及每个点(x,y)在每个脉冲发射时刻η_a、每个距离向采样时刻τ_r的距离向包络w_r,x,y(η_a,τ_r)，其中x表示目标位于第x行(方位向)，y表示目标位于第y列(距离向)。

A、根据脉冲重复频率f_prf、雷达照射中心时刻t₀、信号采样率f_s、参考斜距R_ref，计算第a个脉冲发射时刻η_a和第r个距离门采样时刻τ_r；

η_a＝(a-X/2)/f_prf+t₀,a＝1,…,X (1a)

τ_r＝(r-Y/2)/f_s+2R_ref/c,r＝1,…,Y (1b)

其中，c为光速，取3×10⁸m/s。

B、根据点(x,y)与场景中心位置点的关系以及脉冲重复频率f_prf、雷达照射中心时刻t₀，计算出点(x,y)的零多普勒时刻η_c,x,y；

η_c,x,y＝(x-X/2)/f_prf+t₀ (2)

C、根据点(x,y)与场景中心位置点的关系、参考斜距R_ref、波长λ、信号采样率f_s、其所在距离门y的多普勒中心频率f_d,y和f_r,y以及其零多普勒时刻η_c,x,y，计算出点(x,y)在每个脉冲发射时刻η_a的斜距R_x,y(η_a)；

D、根据雷达速度v、点(x,y)的零多普勒时刻η_c,x,y、其在每个脉冲发射时刻η_a的斜距R_x,y(η_a)、波长λ、天线方位向长度L_a，计算出点(x,y)在每个脉冲发射时刻η_a的方位向包络w_a,x,y(η_a)；

w_{a,x,y_sita}(η_a)＝arc sin[v×(η_a-η_c,x,y)/R_x,y(η_a)] (4a)

其中，arcsin()为反正弦函数，sinc()表示sinc函数。

E、根据点(x,y)在每个脉冲发射时刻η_a的斜距R_x,y(η_a)、脉冲持续时间T_r，计算出点(x,y)在每个脉冲发射时刻η_a、每个距离采样时刻τ_r的距离向包络w_r,x,y(η_a,τ_r)。

步骤三：建立回波数据估计矩阵A，并将回波数据列向量化，具体处理过程分别为：

A、根据距离向包络w_a,x,y(η_a)、方位向包络w_r,x,y(η_a,τ_r)、斜距R_x,y(η_a)、波长λ、距离向调频率K_r、距离采样时刻τ_r，计算估计矩阵A；

矩阵A是一个XY×XY的复数矩阵，矩阵A中的元素A(m,n)方法如下，其中m是行数，n是列数：

m＝a+(r-1)×X,a＝1,…,X r＝1,…,Y (6a)

n＝x+(y-1)×X,x＝1,…,X y＝1,…,Y (6b)

A(m,n)＝w_a,x,y(η_a)w_r,x,y(η_a,τ_r)

(6c)

×exp{-j4πR_x,y(η_a)/λ}exp{jπK_r(τ_r-2R_x,y(η_a)/c)²}

B、将回波数据Echo列向量化，变为echo：

回波数据Echo是一个X×Y的数据矩阵，Echo(x,y)是其第x行、第y列的元素；echo是一个XY×1的数据矩阵，echo(m)是其第m行的元素，则

m＝x+(y-1)×X,x＝1,…,X y＝1,…,Y (7a)

echo(m)＝Echo(x,y) (7b)

步骤四：根据估计矩阵A，运用SLIM估计算法，对列向量化后的回波数据echo进行迭代运算直至收敛，输出场景目标的后向散射系数估计结果。

A、根据估计矩阵A和回波数据echo计算场景目标的初始后向散射系数；

场景目标的初始后向散射系数α_c,0是一个XY×1的复数矩阵，α_c,0(m)是其第m行的元素，则

α_c,0(m)＝A(:,m)^Hecho/[A(:,m)^HA(:,m)] (8)

其中A(:,m)是估计矩阵A第m列的所有元素，(·)^H表示求共轭矩阵。

B、将上一次得到的场景目标的后向散射系数α_c,i-1、回波数据echo和估计矩阵A代入SLIM估计算法方程组，计算估算误差值γ_i，若本次估算误差与上一次迭代中计算得到误差值γ_i-1之比大于0.9，表明收敛，停止迭代；否则重复步骤B；

SLIM估计算法方程组如下：

其中||·||₂表示矩阵的二范数，diag{·}表示矩阵对角化，I为XY×XY的单位矩阵。

式(9a)中第二个方程的矩阵对角化具体过程如下：对角矩阵P是一个XY×XY的矩阵，P(m,n)是其第m行、第n列的元素，则

m＝x+(y-1)×X,x＝1,…,X y＝1,…,Y (9b)

先根据上一次得到的场景目标的后向散射系数α_c,i-1、估计矩阵A和列向量化后的回波数据echo计算收敛值γ_i，然后计算对角矩阵P，最后利用估计矩阵A、对角矩阵P、收敛值γ_i以及列向量化后的回波数据echo计算场景目标新的后向散射系数α_c,i。如果收敛值γ_i不满足要求，则重复以上步骤；如果收敛值γ_i满足要求，获得场景目标最终的后向散射系数α_c，则跳至步骤C。

C、将步骤B中得到的场景目标的后向散射系数α_c去列向量化，得到矩阵化的场景目标的后向散射系数α_T；

步骤B中得到的场景目标的后向散射系数α_c是一个XY×1的复数矩阵，α_c(m)是其第m行的元素；矩阵化的场景目标的后向散射系数α_T是一个X×Y的复数矩阵，α_T(x,y)是其第x行、第y列的元素，则

m＝x+(y-1)×X,x＝1,…,X y＝1,…,Y (10a)

α_T(x,y)＝α_c(m) (10b)

D、将步骤C中得到的场景目标后向散射系数α_T取绝对值，得到SAR图像的灰度值矩阵α_g，输出SAR灰度图像。

SAR图像灰度值矩阵α_g是一个X×Y的实数矩阵，α_g(x,y)是其第x行、第y列的元素，则

α_g(x,y)＝|α_T(x,y)| (11)

其中|·|表示取绝对值。

本发明的优点在于：

本方法能够获得地面目标后向散射系数的精确值，同时能够有效抑制旁瓣能量，提升SAR图像质量，便于SAR图像的判读及后续的进一步应用处理。

附图说明

图1是本发明的方法流程图；

图2是本发明的计算场景中每个点的零多普勒时刻、斜距、方位向包络、距离向包络的流程图；

图3是本发明的计算求解SLIM估计算法方程组的流程图；

图4是本发明的针对单点目标采用匹配滤波成像算法的成像结果；

图5是本发明的针对单点目标采用SLIM成像算法的成像结果。

图6是本发明的针对单点目标采用SLIM成像算法与匹配滤波成像算法的距离向剖面对比图。

图7是本发明的针对单点目标采用SLIM成像算法与匹配滤波成像算法的方位向剖面对比图。

图8是本发明的场景强弱目标分布图。

图9是本发明的针对强弱目标采用匹配滤波成像算法的成像结果。

图10是本发明的针对强弱目标采用SLIM成像算法的成像结果。

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

本发明的一种基于SLIM算法的SAR成像实现方法，流程如图1所示，包括以下几个步骤：

步骤一：读入成像参数和回波信号数据，具体步骤为：

A、读入成像参数，包括斜距参数、多普勒参数和雷达信号参数：参考斜距R_ref、各个距离门y的多普勒中心频率f_d,y、各个距离门y的多普勒调频率f_r,y、距离向调频率K_r、波长λ、信号采样率f_s、脉冲重复频率f_prf、脉冲持续时间T_r、雷达速度v、雷达照射中心时刻t₀、天线方位向长度L_a；

B、读入SAR回波信号数据Echo和SAR回波大小X×Y(表示Echo为一个X行Y列的数据矩阵)；

步骤二：根据SAR回波模型及相关参数，计算出场景中每个点(x,y)在每个脉冲发射时刻η_a的斜距R_x,y(η_a)、方位向包络w_a,x,y(η_a)和每个点(x,y)的零多普勒时刻η_c,x,y，以及每个点(x,y)在每个脉冲发射时刻η_a、每个距离向采样时刻τ_r的距离向包络w_r,x,y(η_a,τ_r)，其中x表示目标位于第x行(方位向)，y表示目标位于第y列(距离向)。

流程如图2所示，具体步骤为：

A、根据脉冲重复频率f_prf、雷达照射中心时刻t₀、信号采样率f_s、参考斜距R_ref，计算第a个脉冲发射时刻η_a和第r个距离门采样时刻τ_r；

η_a＝(a-X/2)/f_prf+t₀,a＝1,…,X (1a)

τ_r＝(r-Y/2)/f_s+2R_ref/c,r＝1,…,Y (1b)

其中，c为光速，取3×10⁸m/s。

B、根据点(x,y)与场景中心位置点的关系以及脉冲重复频率f_prf、雷达照射中心时刻t₀，计算出点(x,y)的零多普勒时刻η_c,x,y；

η_c,x,y＝(x-X/2)/f_prf+t₀ (2)

C、根据点(x,y)与场景中心位置点的关系、参考斜距R_ref、波长λ、信号采样率f_s、其所在距离门y的多普勒中心频率f_d,y和f_r,y以及其零多普勒时刻η_c,x,y，计算出点(x,y)在每个脉冲发射时刻η_a的斜距R_x,y(η_a)；

D、根据雷达速度v、点(x,y)的零多普勒时刻η_c,x,y、其在每个脉冲发射时刻η_a的斜距R_x,y(η_a)、波长λ、天线方位向长度L_a，计算出点(x,y)在每个脉冲发射时刻η_a的方位向包络w_a,x,y(η_a)；

w_{a,x,y_sita}(η_a)＝arc sin[v×(η_a-η_c,x,y)/R_x,y(η_a)] (4a)

其中，arcsin()为反正弦函数，sinc()表示sinc函数。

E、根据点(x,y)在每个脉冲发射时刻η_a的斜距R_x,y(η_a)、脉冲持续时间T_r，计算出点(x,y)在每个脉冲发射时刻η_a、每个距离采样时刻τ_r的距离向包络w_r,x,y(η_a,τ_r)。

步骤三：建立回波数据估计矩阵A，并将回波数据列向量化，具体处理过程为：

A、根据距离向包络w_a,x,y(η_a)、方位向包络w_r,x,y(η_a,τ_r)、斜距R_x,y(η_a)、波长λ、距离向调频率K_r、距离采样时刻τ_r，计算估计矩阵A；

矩阵A是一个XY×XY的复数矩阵，矩阵A中的元素A(m,n)方法如下，其中m是行数，n是列数：

m＝a+(r-1)×X,a＝1,…,X r＝1,…,Y (6a)

n＝x+(y-1)×X,x＝1,…,X y＝1,…,Y (6b)

A(m,n)＝w_a,x,y(η_a)w_r,x,y(η_a,τ_r)

(6c)

×exp{-j4πR_x,y(η_a)/λ}exp{jπK_r(τ_r-2R_x,y(η_a)/c)²}

B、将回波数据Echo列向量化，变为echo；

回波数据Echo是一个X×Y的数据矩阵，Echo(x,y)是其第x行、第y列的元素；echo是一个XY×1的数据矩阵，echo(m)是其第m行的元素，则

m＝x+(y-1)×X,x＝1,…,X y＝1,…,Y (7a)

echo(m)＝Echo(x,y) (7b)

步骤四：根据估计矩阵A，运用SLIM估计算法，对列向量化后的回波数据echo进行迭代运算直至收敛，输出场景目标的后向散射系数估计结果，最终得到SAR灰度图像。

流程如图3所示，具体步骤又分为：

A、根据估计矩阵A和回波数据echo计算场景目标的初始后向散射系数；

场景目标的初始后向散射系数α_c,0是一个XY×1的复数矩阵，α_c,0(m)是其第m行的元素，则

α_c,0(m)＝A(:,m)^Hecho/[A(:,m)^HA(:,m)] (8)

其中A(:,m)是估计矩阵A第m列的所有元素，(·)^H表示求共轭矩阵。

B、将上一次得到的场景目标的后向散射系数α_c,i-1、回波数据echo和估计矩阵A代入SLIM估计算法方程组，计算估算误差值γ_i，若本次估算误差与上一次迭代中计算得到误差值γ_i-1之比大于0.9，表明收敛，停止迭代；否则重复步骤B；

SLIM估计算法方程组如下：

其中||·||₂表示矩阵的二范数，diag{·}表示矩阵对角化，I为XY×XY的单位矩阵。

式(9a)中第二个方程的矩阵对角化具体过程如下：对角矩阵P是一个XY×XY的矩阵，P(m,n)是其第m行、第n列的元素，则

m＝x+(y-1)×X,x＝1,…,X y＝1,…,Y (9b)

先根据上一次得到的场景目标的后向散射系数α_c,i-1、估计矩阵A和列向量化后的回波数据echo计算收敛值γ_i，然后计算对角矩阵P，最后利用估计矩阵A、对角矩阵P、收敛值γ_i以及列向量化后的回波数据echo计算场景目标新的后向散射系数α_c,i。如果收敛值γ_i不满足要求，则重复以上步骤；如果收敛值γ_i满足要求，获得场景目标最终的后向散射系数α_c，则跳至步骤C。

C、将步骤B中得到的场景目标的后向散射系数α_c去列向量化，得到矩阵化的场景目标的后向散射系数α_T；

步骤B中得到的场景目标的后向散射系数α_c是一个XY×1的复数矩阵，α_c(m)是其第m行的元素；矩阵化的场景目标的后向散射系数α_T是一个X×Y的复数矩阵，α_T(x,y)是其第x行、第y列的元素，则

m＝x+(y-1)×X,x＝1,…，X y＝1,…,Y (10a)

α_T(x,y)＝α_c(m) (10b)

D、将步骤C中得到的场景目标后向散射系数α_T取绝对值，得到SAR图像的灰度值矩阵α_g，输出SAR灰度图像。

SAR图像灰度值矩阵α_g是一个X×Y的实数矩阵，α_g(x,y)是其第x行、第y列的元素，则

α_g(x,y)＝|α_T(x,y)| (11)

其中|·|表示取绝对值。

实施例：以一幅大小为128×128的SAR图像为例，具体流程如图1所示，包括以下步骤：

步骤一：读入相关参数和回波数据，具体步骤又分为：

A、读入成像参数，包括斜距参数、多普勒参数和雷达信号参数：参考斜距R_ref、各个距离门y的多普勒中心频率f_d,y、各个距离门y的多普勒调频率f_r,y、距离向调频率K_r、波长λ、信号采样率f_s、脉冲重复频率f_prf、脉冲持续时间T_r、雷达速度v、雷达照射中心时刻t₀、天线方位向长度L_a；

B、读入SAR回波信号数据Echo和SAR回波大小X×Y(表示Echo为一个X行Y列的数据矩阵)；

其中，本实施例中具体参数为：R_ref＝5.5×10³m，K_r＝50×10¹³Hz，λ＝0.03m，f_s＝170×10⁶Hz，f_prf＝210Hz，T_r＝3×10^-7s，v＝150m/s，t₀＝0s，L_a＝1.6m，X＝128，Y＝128，f_d,y与f_r,y是与距离门y有关的参数，f_d,y＝0，f_r,y＝2v²/[λ(R_ref+1.5×10⁸×(y-Y/2)/f_s)]，y＝1,…,Y。

步骤二：根据SAR回波模型及相关参数，计算出场景中每个点(x,y)在每个脉冲发射时刻η_a的斜距R_x,y(η_a)、方位向包络w_a,x,y(η_a)和每个点(x,y)的零多普勒时刻η_c,x,y，以及每个点(x,y)在每个脉冲发射时刻η_a、每个距离向采样时刻τ_r的距离向包络w_r,x,y(η_a,τ_r)，其中x表示目标位于第x行(方位向)，y表示目标位于第y列(距离向)。

流程如图2所示，具体步骤又分为：

A、根据脉冲重复频率f_prf、雷达照射中心时刻t₀、信号采样率f_s、参考斜距R_ref，计算第a个脉冲发射时刻η_a和第r个距离门采样时刻τ_r；

η_a＝(a-X/2)/f_prf+t₀,a＝1,…,X (1a)

τ_r＝(r-Y/2)/f_s+2R_ref/c,r＝1,…,Y (1b)

其中，本实施例中具体参数为：f_prf＝210Hz，t₀＝0s，X＝128，R_ref＝5.5×10³m，f_s＝170×10⁶Hz，Y＝128，c＝3×10⁸m/s，根据a的不同取值得到η_a，根据r的不同取值得到τ_r。

B、根据点(x,y)与场景中心位置点的关系以及脉冲重复频率f_prf、雷达照射中心时刻t₀，计算出点(x,y)的零多普勒时刻η_c,x,y；

η_c,x,y＝(x-X/2)/f_prf+t₀ (2)

其中，本实施例中具体参数为：f_prf＝210Hz，t₀＝0s，X＝128，根据x的不同取值得到η_c,x,y。

C、根据点(x,y)与场景中心位置点的关系、参考斜距R_ref、波长λ、信号采样率f_s、其所在距离门y的多普勒中心频率f_d,y和f_r,y以及其零多普勒时刻η_c,x,y，计算出点(x,y)在每个脉冲发射时刻η_a的斜距R_x,y(η_a)；

其中，本实施例中具体参数为：R_ref＝5.5×10³m，c＝3×10⁸m/s，Y＝128，f_s＝170×10⁶Hz，λ＝0.03m，f_d,y与f_r,y由步骤一中读入参数得到，η_a按公式(1a)获得，η_c,x,y按公式(2)获得，根据x，y，a的不同取值得到R_x,y(η_a)，x＝1,…,X，y＝1,…,Y，a＝1,…,X。

D、根据雷达速度_v、点(x,y)的零多普勒时刻η_c,x,y、其在每个脉冲发射时刻η_a的斜距R_x,y(η_a)、波长λ、天线方位向长度L_a，计算出点(x,y)在每个脉冲发射时刻η_a的方位向包络w_a,x,y(η_a)；

w_{a,x,y_sita}(η_a)＝arc sin[v×(η_a-η_c,x,y)/R_x,y(η_a)] (4a)

其中，arcsin()为反正弦函数，sinc()表示sinc函数。

其中，本实施例中具体参数为：v＝150m/s，λ＝0.03m，L_a＝1.6m，η_a按公式(1a)获得，η_c,x,y按公式(2)获得，R_x,y(η_a)按公式(3)获得，根据x，y，a的不同取值得到w_a,x,y(η_a)，x＝1,…,X，y＝1,…,Y，a＝1,…,X。

E、根据点(x,y)在每个脉冲发射时刻η_a的斜距R_x,y(η_a)、脉冲持续时间T_r，计算出点(x,y)在每个脉冲发射时刻η_a、每个距离采样时刻τ_r的距离向包络w_r,x,y(η_a,τ_r)。

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

其中，本实施例中具体参数为：T_r＝3×10^-7s，τ_r按公式(1b)获得，R_x,y(η_a)按公式(3)获得，根据x，y，a，_r的不同取值得到w_r,x,y(η_a,τ_r)，x＝1,…,X，y＝1,…,Y，a＝1,…,X，r＝1,…,Y。

步骤三：建立回波数据估计矩阵A，并将回波数据列向量化，具体处理过程分别为：

A、根据距离向包络w_a,x,y(η_a)、方位向包络w_r,x,y(η_a,τ_r)、斜距R_x,y(η_a)、波长λ、距离向调频率K_r、距离采样时刻τ_r，计算估计矩阵A；

矩阵A是一个XY×XY的复数矩阵，矩阵A中的元素A(m,n)方法如下，其中m是行数，n是列数：

m＝a+(r-1)×X,a＝1,…,X r＝1,…,Y (6a)

n＝x+(y-1)×X,x＝1,…,X y＝1,…,Y (6b)

A(m,n)＝w_a,x,y(η_a)w_r,x,y(η_a,τ_r)

(6c)

×exp{-j4πR_x,y(η_a)/λ}exp{jπK_r(τ_r-2R_x,y(η_a)/c)²}

其中，本实施例中具体参数为：f₀＝c/λ＝10¹⁰Hz，K_r＝50×10¹³Hz，c＝3×10⁸m/s，X＝128，Y＝128，wa_x,y(η_a)按公式(4a～4b)获得，wr_x,y(η_a,τ_r)按公式(5)获得，R_x,y(η_a)按公式(3)获得，τ_r按公式(1b)获得，根据m，n的不同取值得到A(m,n)。

B、将回波数据Echo列向量化，变为echo；

回波数据Echo是一个X×Y的数据矩阵，Echo(x,y)是其第x行、第y列的元素；echo是一个XY×1的数据矩阵，echo(m)是其第m行的元素，则

m＝x+(y-1)×X,x＝1,…,X y＝1,…,Y (7a)

echo(m)＝Echo(x,y) (7b)

其中，本实施例中具体参数为：X＝128，Y＝128，Echo是步骤一中读入的回波数据。

步骤四：根据估计矩阵A，运用SLIM估计算法，对列向量化后的回波数据echo进行迭代运算直至收敛，输出场景目标的后向散射系数估计结果，最终得到SAR灰度图像。

流程如图3所示，具体步骤又分为：

A、根据估计矩阵A和回波数据echo计算场景目标的初始后向散射系数；

场景目标的初始后向散射系数α_c,0是一个XY×1的复数矩阵，α_c,0(m)是其第m行的元素，则

α_c,0(m)＝A(:,m)^Hecho/[A(:,m)^HA(:,m)] (8)

其中A(:,m)是估计矩阵A第m列的所有元素，(·)^H表示求共轭矩阵。

其中，本实施例中具体参数为：A按公式(6a～6c)获得，echo按公式(7a～7b)获得。

B、将上一次得到的场景目标的后向散射系数α_c,i-1、回波数据echo和估计矩阵A代入SLIM估计算法方程组，计算估算误差值γ_i，若本次估算误差与上一次迭代中计算得到误差值γ_i-1之比大于0.9，表明收敛，停止迭代；否则重复步骤B；

SLIM估计算法方程组如下：

其中||·||₂表示矩阵的二范数，diag{·}表示矩阵对角化，I为XY×XY的单位矩阵。

式(9a)中第二个方程的矩阵对角化具体过程如下：对角矩阵P是一个XY×XY的矩阵，P(m,n)是其第m行、第n列的元素，则

m＝x+(y-1)×X,x＝1,…,X y＝1,…,Y (9b)

先根据上一次得到的场景目标的后向散射系数α_c,i-1、估计矩阵A和列向量化后的回波数据echo计算收敛值γ_i，然后计算对角矩阵P，最后利用估计矩阵A、对角矩阵P、收敛值γ_i以及列向量化后的回波数据echo计算场景目标新的后向散射系数α_c,i。如果收敛值γ_i不满足要求，则重复以上步骤；如果收敛值γ_i满足要求，获得场景目标最终的后向散射系数α_c，则跳至步骤C。

其中，本实施例中具体参数为：A按公式(6a～6c)获得，echo按公式(7a～7b)获得。

C、将B中得到的场景目标的后向散射系数α_c去列向量化，得到矩阵化的场景目标的后向散射系数α_T；

B中得到的场景目标的后向散射系数α_c是一个XY×1的复数矩阵，α_c(m)是其第m行的元素；矩阵化的场景目标的后向散射系数α_T是一个X×Y的复数矩阵，α_T(x,y)是其第x行、第y列的元素，则

m＝x+(y-1)×X,x＝1,…,X y＝1,…,Y (10a)

α_T(x,y)＝α_c(m) (10b)

其中，本实施例中具体参数为：X＝128，Y＝128，α_c按公式(9)获得。

D、将步骤C中得到的场景目标后向散射系数α_T取绝对值，得到SAR图像的灰度值矩阵α_g，输出SAR灰度图像；

SAR图像灰度值矩阵α_g是一个X×Y的实数矩阵，α_g(x,y)是其第x行、第y列的元素，则

α_g(x,y)＝|α_T(x,y)| (11)

其中|·|表示取绝对值。

为了说明该方法的有效性，进行如下仿真试验，仿真参数与步骤一中读入参数相同。针对单点目标进行仿真，单点目标后向散射系数设为1，匹配滤波成像算法结果如图4所示，本发明SLIM成像算法结果如图5所示，距离向剖面对比如图6所示，方位向剖面对比如图7所示，距离向和方位向峰值旁瓣比如表1所示；针对多点目标进行仿真，在场景中布置5个点目标如图8所示，四个弱目标后向散射系数设为0.1，一个强目标后向散射系数设为1，四个弱目标分别在距离向和方位向与强目标相间隔一个像素，匹配滤波成像算法结果如图9所示，本发明SLIM成像算法结果如图10所示。

表1峰值旁瓣比

峰值旁瓣比(dB)	匹配滤波成像算法	本发明SLIM成像算法
			方位向	-13.18	<-60
距离向	-13.26	<-60

从图4、图5、图6、图7可以看出：传统匹配滤波算法具有较高的旁瓣能量，其峰值旁瓣比在-13.2dB左右。相比较而言，本发明方法所获得处理结果明显优于传统匹配滤波算法，其峰值旁瓣比远小于匹配滤波算法，旁瓣能量近似可以忽略。

从图8、图9、图10可以看出：由于传统匹配滤波处理算法具有较高的旁瓣能量，导致弱目标被强目标的旁瓣淹没，影响图像的判读。相比较而言，本发明方法的旁瓣能量极低，不会出现强弱目标之间的相互干扰，有效提升图像质量。

Claims (3)

1.一种基于SLIM算法的SAR成像实现方法，包括以下几个步骤：

步骤一：读入成像参数和回波信号数据，具体步骤为：

A、读入成像参数，包括斜距参数、多普勒参数和雷达信号参数：参考斜距R_ref、各个距离门y的多普勒中心频率f_d,y、各个距离门y的多普勒调频率f_r,y、距离向调频率K_r、波长λ、信号采样率f_s、脉冲重复频率f_prf、脉冲持续时间T_r、雷达速度v、雷达照射中心时刻t₀、天线方位向长度L_a；

B、读入SAR回波信号数据Echo和SAR回波大小X×Y；

步骤二：根据SAR回波模型及相关参数，计算出场景中每个点(x,y)在每个脉冲发射时刻η_a的斜距R_x,y(η_a)、方位向包络w_a,x,y(η_a)和每个点(x,y)的零多普勒时刻η_c,x,y，以及每个点(x,y)在每个脉冲发射时刻η_a、每个距离向采样时刻τ_r的距离向包络w_r,x,y(η_a,τ_r)，其中x表示目标位于第x行，y表示目标位于第y列；

步骤三：建立回波数据估计矩阵A，并将回波数据列向量化，具体处理过程为：

A、根据方位向包络w_a,x,y(η_a)、距离向包络w_r,x,y(η_a,τ_r)、斜距R_x,y(η_a)、波长λ、距离向调频率K_r、距离采样时刻τ_r，计算估计矩阵A；

矩阵A为XY×XY的复数矩阵，矩阵A中的元素A(m,n)方法如下，其中m是行数，n是列数：

m＝a+(r-1)×X,a＝1,…,X r＝1,…,Y (6a)

n＝x+(y-1)×X,x＝1,…,X y＝1,…,Y (6b)

$\begin{array}{c}A(m,n)=w_{a,x,y}\left({\mathrm{\η}}_a\right)w_{r,x,y}({\mathrm{\η}}_a,{\mathrm{\τ}}_r)\\ \×\exp \{-j4{\mathrm{\πR}}_{x,y}\left({\mathrm{\η}}_a\right)/\mathrm{\λ}\}\exp \left\{{\mathrm{j\πK}}_r{({\mathrm{\τ}}_r-2R_{x,y}({\mathrm{\η}}_a)/c)}^2\right\}\end{array}---\left(6c\right)$

B、将回波数据Echo列向量化，变为echo；

回波数据Echo为X×Y的数据矩阵，Echo(x,y)是其第x行、第y列的元素；echo是一个XY×1的数据矩阵，echo(m)是其第m行的元素，则

m＝x+(y-1)×X,x＝1,…,X y＝1,…,Y (7a)

echo(m)＝Echo(x,y) (7b)

步骤四：根据估计矩阵A，运用SLIM估计算法，对列向量化后的回波数据echo进行迭代运算直至收敛，输出场景目标的后向散射系数估计结果，最终得到SAR灰度图像。

2.根据权利要求1所述的一种基于SLIM算法的SAR成像实现方法，所述步骤二具体步骤为：

A、根据脉冲重复频率f_prf、雷达照射中心时刻t₀、信号采样率f_s、参考斜距R_ref，计算第a个脉冲发射时刻η_a和第r个距离门采样时刻τ_r；

η_a＝(a-X/2)/f_prf+t₀,a＝1,…,X (1a)

τ_r＝(r-Y/2)/f_s+2R_ref/c,r＝1,…,Y (1b)

其中，c为光速，取3×10⁸m/s；

B、根据点(x,y)与场景中心位置点的关系以及脉冲重复频率f_prf、雷达照射中心时刻t₀，计算出点(x,y)的零多普勒时刻η_c,x,y；

η_c,x,y＝(x-X/2)/f_prf+t₀ (2)

C、根据点(x,y)与场景中心位置点的关系、参考斜距R_ref、波长λ、信号采样率f_s、其所在距离门y的多普勒中心频率f_d,y和f_r,y以及其零多普勒时刻η_c,x,y，计算出点(x,y)在每个脉冲发射时刻η_a的斜距R_x,y(η_a)；

$R_{x,y}\left({\mathrm{\η}}_a\right)=R_{ref}+\frac{c}{2}\×(y-\frac{Y}{2})/f_s+\frac{\mathrm{\λ}}{2}\[f_{d,y}({\mathrm{\η}}_a-{\mathrm{\η}}_{c,x,y})+\frac{1}{2}{f}_{r,y}{({\mathrm{\η}}_a-{\mathrm{\η}}_{c,x,y})}^2\]---\left(3\right)$

D、根据雷达速度v、点(x,y)的零多普勒时刻η_c,x,y、其在每个脉冲发射时刻η_a的斜距R_x,y(η_a)、波长λ、天线方位向长度L_a，计算出点(x,y)在每个脉冲发射时刻η_a的方位向包络w_a,x,y(η_a)；

w_{a,x,y_sita}(η_a)＝arc sin[v×(η_a-η_c,x,y)/R_x,y(η_a)] (4a)

$w_{a,x,y}\left({\mathrm{\η}}_a\right)=\sin c\left(\frac{w_{a,x,y\_sita}\left({\mathrm{\η}}_a\right)L_a}{\mathrm{\λ}}\right)---\left(4b\right)$

其中，arc sin()为反正弦函数，sin c()表示sinc函数；

E、根据点(x,y)在每个脉冲发射时刻η_a的斜距R_x,y(η_a)、脉冲持续时间T_r，计算出点(x,y)在每个脉冲发射时刻η_a、每个距离采样时刻τ_r的距离向包络w_r,x,y(η_a,τ_r)；

$w_{r,x,y}({\mathrm{\&eta;}}_a,{\mathrm{\&tau;}}_r)=\left\{\begin{array}{cc}1,& abs\&lsqb;({\mathrm{\&tau;}}_r-2R_{x,y}({\mathrm{\&eta;}}_a)/2)\&rsqb;<T_r/2\\ 0,& abs\&lsqb;({\mathrm{\&tau;}}_r-2R_{x,y}({\mathrm{\&eta;}}_a)/2)\&rsqb;\&GreaterEqual;T_r/2\end{array}\right.---\left(5\right).$

3.根据权利要求1所述的一种基于SLIM算法的SAR成像实现方法，所述步骤四具体步骤为：

A、根据估计矩阵A和回波数据echo计算场景目标的初始后向散射系数；

场景目标的初始后向散射系数α_c,0为XY×1的复数矩阵，α_c,0(m)是其第m行的元素，则

α_c,0(m)＝A(:,m)^Hecho/[A(:,m)^HA(:,m)] (8)

其中A(:,m)是估计矩阵A第m列的所有元素，(·)^H表示求共轭矩阵；

B、将上一次得到的场景目标的后向散射系数α_c,i-1、回波数据echo和估计矩阵A代入SLIM估计算法方程组，计算估算误差值γ_i，若本次估算误差与上一次迭代中计算得到误差值γ_i-1之比大于0.9，表明收敛，停止迭代；否则重复步骤B；

SLIM估计算法方程组如下：

$\left\{\begin{array}{c}{\mathrm{\&gamma;}}_i=\frac{1}{XY}\left|\right|echo-{\mathrm{A\&alpha;}}_{c,i-1}|{|}_2^2\\ P=diag\left\{\right|{\mathrm{\&alpha;}}_{c,i-1}\left(1\right){|}^2,|{\mathrm{\&alpha;}}_{c,i-1}\left(2\right){|}^2,...,|{\mathrm{\&alpha;}}_{c,i-1}\left(XY\right){|}^2\}\\ {\mathrm{\&alpha;}}_{c,i}={\mathrm{PA}}^H{({\mathrm{APA}}^H+{\mathrm{\&gamma;}}_{i}I)}^{-1}echo\end{array}\right.---\left(9a\right)$

其中||·||₂表示矩阵的二范数，diag{·}表示矩阵对角化，I为XY×XY的单位矩阵；

式(9a)中第二个方程的矩阵对角化具体过程如下：对角矩阵P为XY×XY的矩阵，P(m,n)是其第m行、第n列的元素，则

m＝x+(y-1)×X,x＝1,…,X y＝1,…,Y (9b)

$P(m,n)=\left\{\begin{array}{cc}|{\mathrm{\&alpha;}}_{c,i-1}\left(xy\right){|}^2,& m=n\\ 0,& m\&NotEqual;n\end{array}\right.---\left(9c\right)$

先根据上一次得到的场景目标的后向散射系数α_c,i-1、估计矩阵A和列向量化后的回波数据echo计算收敛值γ_i，然后计算对角矩阵P，最后利用估计矩阵A、对角矩阵P、收敛值γ_i以及列向量化后的回波数据echo计算场景目标新的后向散射系数α_c,i；如果收敛值γ_i不满足要求，则重复以上步骤；如果收敛值γ_i满足要求，获得场景目标最终的后向散射系数α_c，则跳至步骤C；

C、将步骤B中得到的场景目标的后向散射系数α_c去列向量化，得到矩阵化的场景目标的后向散射系数α_T；

步骤B中得到的场景目标的后向散射系数α_c为XY×1的复数矩阵，α_c(m)是其第m行的元素；矩阵化的场景目标的后向散射系数α_T是一个X×Y的复数矩阵，α_T(x,y)是其第x行、第y列的元素，则

m＝x+(y-1)×X,x＝1,…,X y＝1,…,Y (10a)

α_T(x,y)＝α_c(m) (10b)

D、将步骤C中得到的场景目标后向散射系数α_T取绝对值，得到SAR图像的灰度值矩阵α_g，输出SAR灰度图像；

SAR图像灰度值矩阵α_g是一个X×Y的实数矩阵，α_g(x,y)是其第x行、第y列的元素，则

α_g(x,y)＝|α_T(x,y)| (11)

其中|·|表示取绝对值。