CN104714230B - 稀疏合成孔径雷达多视成像方法 - Google Patents

Abstract

本发明提供了一种稀疏合成孔径雷达多视成像的方法。该方法采用基于频域运算的回波模拟算子，建立多视观测与欠采样回波数据之间的观测模型，实现了多视技术和压缩感知技术的融合，实现欠采样下稀疏目标场景的多视成像。

Description

稀疏合成孔径雷达多视成像方法

技术领域

本发明涉及电子行业雷达技术领域，尤其涉及一种基于回波模拟算子及结构稀疏正则化的稀疏合成孔径雷达(Synthetic Aperture Radar，简称SAR)多视成像方法。

背景技术

合成孔径雷达作为一种主动式微波成像系统，具有全天时、全天候和高分辨率成像等特点。然而，在SAR中，由于每一分辨单元内的各散射点之间的相对相位与雷达视角紧密相关，其相干叠加是随机的，因而会产生斑点噪声。这种噪声表现为一种复杂的乘性噪声，不仅影响图像的性噪比，更极大的影响了目标的识别。多视处理技术是现有斑点噪声抑制中最为常用的技术之一，它通过对信号频谱的不同部分进行成像处理，并对幅度进行非相干叠加，抑制斑点噪声，提高辐射分辨率。

最近兴起的稀疏SAR成像方法旨在利用压缩感知理论与方法，通过远小于奈奎斯特采样定律要求的采样数(即欠采样)下恢复稀疏信号。一般来说，欠采样是在空域上实现的。因而在此情况下，回波频谱将会混叠且难以分割，导致多视技术无法直接应用压缩感知SAR成像方法中。

发明内容

(一)要解决的技术问题

为解决上述的一个或多个问题，本发明提供了一种基于回波模拟算子及结构稀疏正则化的稀疏合成孔径雷达多视成像，以实现欠采样下稀疏目标场景的多视成像。

(二)技术方案

本发明提供了一种稀疏合成孔径雷达多视成像的方法。该方法包括：步骤A，接收合成孔径雷达回波数据Y_s，构建相应的稀疏合成孔径的雷达多视观测方程：

Y_s＝Θ·M(X)

其中，X是经过子视分割的雷达子视图像的集合，M(·)是多视观测算子，Θ为采样矩阵；步骤B，根据多视观测方程，建立基于结构稀疏正则化的稀疏合成孔径雷达多视成像模型；步骤C，采用阈值迭代算法求解稀疏合成孔径雷达多视成像模型，重建获得雷达子视图像集合X；以及步骤D，对雷达子视图像集合X的幅度进行非相干叠加，获得多视输出图像Z。

(三)有益效果

从上述技术方案可以看出，本发明稀疏合成孔径雷达多视成像具有以下有益效果：

(1)采用基于频域运算的回波模拟算子，建立多视观测与欠采样回波数据之间的观测模型，实现了多视技术和压缩感知技术的融合，实现欠采样下稀疏目标场景的多视成像；

(2)基于多视输出图像特点，构造基于幅度非相干叠加的结构稀疏正则化恢复模型，可以在低于奈奎斯特采样率下恢复实现稀疏目标场景的多视成像。

附图说明

图1为本发明实施例稀疏合成孔径雷达多视成像的流程图；

图2为采用本实施例稀疏SAR多视成像方法实现不同视数、不同采样率下的SAR数据二维仿真结果。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。需要说明的是，在附图或说明书描述中，相似或相同的部分都使用相同的图号。附图中未绘示或描述的实现方式，为所属技术领域中普通技术人员所知的形式。

本发明稀疏合成孔径雷达多视成像中，输入的数据为SAR回波数据，最终输出的数据为观测场景后向散射幅度系数矩阵。本发明稀疏合成孔径雷达多视成像，解决了多视技术与压缩感知技术融合的问题，实现欠采样下稀疏目标场景的多视成像。

在本发明的一个示例性实施例中，提供了一种稀疏合成孔径雷达多视成像方法。如图1所示，该稀疏合成孔径雷达多视成像包括：

步骤A，接收合成孔径雷达回波数据Y_s，构建相应的稀疏合成孔径雷达多视观测方程：

Y_s＝ΘM(X)

其中，X是包含所有子视图像的集合，Θ为采样矩阵。

具体而言，该步骤A又可以包括：

子步骤A1，对子视图像X_il(i＝1，…，I；l＝1，2，…，L)，按照以下公式获得信号

其中，矩阵F_η第k_a行第n_a列上的元素为π为圆周率，k_a＝1，2，…，N_a/I，n_a＝1，2，…，N_a/I；N_a为方位向采样数；为全频谱多普勒频率；为方位解压算子，λ为波长，v为雷达平台相对于目标场景的运动速度，R₀(τ)为位于距离时间τ处雷达与目标的最短斜距；

子步骤A2，对进行方位向/距离向二维耦合，获得信号

其中，为二维耦合后的距离向时间，为距离徙动量，sinc(x)＝sin(x)/x，sin(x)为正弦函数；

子步骤A3，对按照以下公式获得信号S₃：

其中，F_τ第n_r行第k_r列上的元素k_r＝1，2，…，N_r/L，n_r＝1，2，…，N_r/L；N_r为距离向采样数；为距离向逆傅里叶变换矩阵，为方位向逆傅里叶变换矩阵，conj(·)表示取共轭；

子步骤A4，对信号S₃进行稀疏采样，完成多视雷达观测方程：

Y_s＝ΘS₃

具体来讲，上述方程为：

其中，η，η_l代表子视图像和数据方位向时间，τ，τ_l代表子视图像和数据距离向时间，F和F^-1分别表示傅里叶变换矩阵及其逆矩阵，P为方位向频域滤波器，Q为距离向频域滤波器，“ο”代表矩阵间哈达马乘积的运算符，Λ(·)为距离徙动算子，I，L分别为方位向和距离向视数。

步骤B，根据所述多视观测方程，建立基于结构稀疏正则化的稀疏合成孔径雷达多视成像模型：

其中，ω为正则化参数，‖X‖_G为结构稀疏范数，arg min是最小化计算式。

步骤C，采用阈值迭代算法求解多视SAR成像模型，重建获得雷达子视图像集合X；

具体而言，该步骤C又可以包括：

子步骤C1：初始化SAR回波数据对应场景的目标散射强度X⁰，正则化参数ω和迭代终止准则ε，令k＝0；

子步骤C2：更新梯度下降序列B^k，由下式实现：

其中，X^k为第k次迭代的子视图像的集合，Y_s表示真实回波数据，μ_k表示梯度下降的步长，M(·)为步骤2所定义的回波模拟算子，G_ij(·)：

其中，分别为方位向距离向频域滤波器，为P，Q的复共轭，φ_j选取矩阵的第[(i-1)N_a/I+1：i·N_a/I]行或者第[(j-1)N_r/L+1：j·N_r/L]列，Γ(·)为距离徙动校正算子：

其中，分别为距离徙动校正前后的数据，为距离徙动校正前的方位频率和距离向时间；

子步骤C3：更新目标场景散射强度X_k+1：X^k+1＝H(B^k)，其中，H(·)为阈值算子使得：

X^k+1(m，n)＝h(X^k(m，n))

其中，h(·)为阈值函数，X^k(m，n)为所有组成的向量，并且对任意

其中，sign(·)为取相位函数；

子步骤C4：若‖X^k+1-X^k‖_F/‖X^k‖_F＜ε，其中‖·‖_F为F范数，则迭代终止，算法输出为观测场景的后向散射系数矩阵X^*＝X^k+1；否则，令X^k＝X^k+1，k＝k+1，转回步骤C2。

步骤D，对子视图像幅度进行非相干叠加，获得多视输出图像Z：

申请人运用本实施例的方法对不同视数、不同采样率下的SAR数据二维仿真实验，在该仿真实验中，主要雷达参数如下：景中心距离R₀＝8000m，飞行器速度V＝150m/s，信号带宽B_r＝50MHz，脉冲持续时间T_r＝2μs，载频f₀＝3GHz，脉冲发射频率PRF＝75Hz，合成孔径时间T_s＝1.333s，信噪比SNR＝5dB。

图2为采用本实施例稀疏SAR多视成像方法实现不同视数、不同采样率下的SAR数据二维仿真结果。图2中(a)为RDA在满采样下一视(左)和三视(右)成像结果；图2中(b)为RDA在30％采样下一视(左)和三视(右)成像结果；图2中(c)为所提出的多视稀疏SAR成像方法在30％采样下一视(左)和三视(右)成像结果。

至此，已经结合附图对本实施例进行了详细描述。依据以上描述，本领域技术人员应当对本发明稀疏合成孔径雷达多视成像方法有了清楚的认识。

此外，上述对各元件的定义并不仅限于实施方式中提到的各种具体结构或形状，本领域的普通技术人员可对其进行简单地熟知地替换。

综上所述，本发明采用基于频域运算的回波模拟算子，可以建立多视观测与欠采样回波数据之间的观测模型，实现了多视技术和压缩感知技术的融合，实现欠采样下稀疏目标场景的多视成像，可广泛应用于雷达观测的各个领域。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种稀疏合成孔径雷达多视成像的方法，其特征在于，包括：

步骤A，接收合成孔径雷达回波数据Y_s，构建相应的稀疏合成孔径的雷达多视观测方程：

Y_s＝Θ·M(X)

其中，X是经过子视分割的雷达子视图像的集合，M(·)是多视观测算子，Θ为采样矩阵；

步骤B，根据所述多视观测方程，建立基于结构稀疏正则化的稀疏合成孔径雷达多视成像模型；

步骤C，采用阈值迭代算法求解所述稀疏合成孔径雷达多视成像模型，重建获得雷达子视图像集合X；以及

步骤D，对雷达子视图像集合X的幅度进行非相干叠加，获得多视输出图像Z；

其中，步骤A构建的多视观测方程为：

其中，X_il为X中方位向第i块、距离向第l块雷达子视图像，各子视图像对应到回波数据域的频谱互不重叠；F和F^-1分别表示傅里叶变换矩阵及其逆变换矩阵，其下标η和τ代表子视图像方位向和距离向时间，η_s和τ_s代表回波数据方位向和距离向时间；P为方位向频域滤波器，Q为距离向频域滤波器，符号代表矩阵哈达玛乘积，Λ(·)为方位向/距离向二维耦合算子，I和L分别为方位向和距离向多视数。

2.根据权利要求1所述的方法，其特征在于，所述步骤A包括：

子步骤A1，对子视图像X_il(i＝1,…,I；l＝1,2,…,L)，按照以下公式获得信号

$P(f_{{\mathrm{\η}}_s},\mathrm{\τ})=e^{{\mathrm{j\π\λR}}_0\left(\mathrm{\τ}\right)f_{{\mathrm{\η}}_s}^2/2v^2}$

其中，矩阵F_η第k_a行第n_a列上的元素为π为圆周率，k_a＝1,2,…,N_a/I，n_a＝1,2,…,N_a/I；N_a为方位向采样数；为全频谱多普勒频率；为方位解压算子，λ为波长，v为雷达平台相对于目标场景的运动速度，R₀(τ)为位于距离时间τ处雷达与目标的最短斜距；

子步骤A2，对进行方位向/距离向二维耦合，获得信号

$S_2^l(f_{{\mathrm{\η}}_s},\widetilde{\mathrm{\τ}})=\mathrm{\Λ}\left(S_1^l\right(f_{{\mathrm{\η}}_s},\mathrm{\τ}\left)\right)=\underset{\mathrm{\τ}}{\Σ}{S}_1^l\left(f_{{\mathrm{\η}}_s},\mathrm{\τ}\right)\sin c(\widetilde{\mathrm{\τ}}-(\mathrm{\τ}+2\mathrm{\Δ}R\left(f_{{\mathrm{\η}}_s},\mathrm{\τ}\right)/c\left)\right)$

其中，为二维耦合后的距离向时间，为距离徙动量，sinc(x)＝sin(x)/x，sin(x)为正弦函数；

子步骤A3，对按照以下公式获得信号S₃：

$Q({\mathrm{\η}}_s,f_{{\mathrm{\τ}}_s})=e^{{\mathrm{j\πf}}_{{\mathrm{\τ}}_s}^2/K_r}$

其中，F_τ第n_r行第k_r列上的元素n_r＝1,2,…,N_r/L；N_r为距离向采样数；为距离向逆傅里叶变换矩阵，为方位向逆傅里叶变换矩阵,conj(·)表示取共轭；

子步骤A4，对信号S₃进行稀疏采样，完成多视雷达观测方程：

Y_s＝ΘS₃。

3.根据权利要求1所述的方法，其特征在于，所述步骤B建立基于稀疏正则化的稀疏合成孔径雷达多视成像模型为：

$\underset{X}{argmin}\left(\right||Y_s-\mathrm{\Θ}M(X)|{|}_2^2+\mathrm{\ω}|\left|X\right|{|}_G)$

$\left|\right|X|{|}_G=\underset{m=1}{\overset{N_a/I}{\Σ}}\underset{n=1}{\overset{N_r/L}{\Σ}}{\left(\underset{i=1}{\overset{I}{\Σ}}\underset{l=1}{\overset{L}{\Σ}}\right|X_{il}(m,n){|}^2)}^{1/2}$

其中，ω为正则化参数；arg min是最小化计算式；X_il(m,n)为子视图像X_il中的元素。

4.根据权利要求1所述的方法，其中步骤C采用阈值迭代算法求解所述多视SAR成像模型，重建获得雷达子视图像集合X的具体步骤包括：

子步骤C1：初始化所述SAR回波数据对应场景的目标散射强度X⁰、正则化参数ω和迭代终止准则ε，令迭代数k＝0；

子步骤C2：更新梯度下降序列由下式实现：

$B_{il}^k=X_{il}^k+{\mathrm{\μ}}_k{G}_{il}\left({\mathrm{\Θ}}^T\right(Y_s-\mathrm{\Θ}M\left(X^k\right)\left)\right)$

其中，为第k次迭代的方位向第i块、距离向第l块雷达子视图像，μ_k表示梯度下降的步长，和分别为方位向距离向频域滤波器，为P和Q的复共轭，分别表示选取矩阵的第[(i-1)N_a/I+1:i·N_a/I]行或者第[(j-1)N_r/L+1:j·N_r/L]列，Γ(·)为距离徙动校正算子：

$\mathrm{\Γ}\left(\tilde{S}\right({\tilde{f}}_l,\widetilde{\mathrm{\τ}}\left)\right)=\underset{\widetilde{\mathrm{\τ}}}{\Σ}\tilde{S}({\tilde{f}}_{{\mathrm{\η}}_l},\widetilde{\mathrm{\τ}})\sin c(\widetilde{\mathrm{\τ}}-(\mathrm{\τ}+2\mathrm{\Δ}R\left(f_{{\mathrm{\η}}_l},\mathrm{\τ}\right)/c)$

其中，分别为距离徙动校正前后的数据，为距离徙动校正前的方位频率和距离向时间；

子步骤C3：更新目标场景散射强度其中，H(·)为阈值算子使得：

$X_{il}^{k+1}(m,n)=h\left(X_{il}^k\right(m,n\left)\right)$

其中，h(·)为阈值函数，并且对任意

$h\left(x\right)=\left\{\begin{array}{cc}sign\left(x\right)\left(\right|x|-{\mathrm{\ω\μ}}_k),& \left|x\right|>{\mathrm{\ω\μ}}_k\\ 0,& \left|x\right|\≤{\mathrm{\ω\μ}}_k\end{array}\right.$

其中，sign(·)为取相位函数；再由所有组成

子步骤C4：若||X^k+1-X^k||_F/||X^k||_F<ε，其中||·||_F为F范数，则迭代终止，算法输出为雷达子视图像集合X＝X^k+1；否则，令X^k＝X^k+1，k＝k+1，转回步骤C2。

5.根据权利要求2所述的方法，其中，所述步骤D对雷达子视图像集合X进行幅度非相干叠加，输出多视雷达图像：

$Z=\sqrt{\underset{i=1}{\overset{I}{\&Sigma;}}\underset{l=1}{\overset{L}{\&Sigma;}}|X_{il}{|}^2}.$