CN105891827A - 一种机载mimo-sar下视三维稀疏成像方法 - Google Patents

Abstract

本发明提供一种机载MIMO‑SAR下视三维稀疏成像方法，包括：第一步：分析MIMO‑SAR成像处理过程，构建了MIMO‑SAR回波数据的稀疏表示模型，并对回波数据进行降采样处理；第二步：依据稀疏模型对数据进行重构，获得经过高度向与方位向补偿和匹配滤波的结果，并对结果进行俯仰向成像处理，实现场景的三维成像。该方法能够利用少量数据进行重构，实现高分辨率三维成像，有效地降低数据量，为解决MIMO‑SAR数据量过大问题提供一种有力的方法。

Description

一种机载MIMO-SAR下视三维稀疏成像方法

技术领域

本发明涉及信号与信息处理技术，具体涉及一种机载MIMO-SAR下视三维稀疏成像方法。

背景技术

机载下视多输入多输出合成孔径雷达(Multiple Input Multiple OutputSynthetic Aperture Radar,MIMO-SAR)通过采用多个发射天线阵元和多个接收天线阵元形成虚拟阵元的方法，有效地减少了物理阵元的数量，并且采用下视三维成像技术，实现了场景的三维成像。高度向上，机载下视MIMO-SAR采用脉冲压缩技术实现高度向上的高分辨率成像；方位向上，利用载机运动形成的合成孔径实现方位向上的高分辨率成像；俯仰向上，通过虚拟阵元形成的合成孔径实现俯仰向的高分辨率成像。随着MIMO-SAR分辨率的不断提高，采用奈奎斯特(Nyquist)采样定理所生成的SAR数据给数据存储和传输带来了巨大的挑战。压缩感知(Compressive Sensing,CS)理论的提出使得信号处理与获取等方面发生了革命性变化，并为解决高分辨率雷达数据量过大的问题提供了思路，CS理论指出，当信号具有稀疏性时，可以用远低于Nyquist采样定理所要求的采样速率来采集信号，并通过重构算法准确重构出原始信号。

现有采用压缩感知对机载下视MIMO-SAR进行成像处理的技术，如基于压缩感知的下视三维SAR成像新方法(参见李学仕的《基于压缩感知的下视三维SAR成像新方法》，发表在《电子与信息学报》，2012，第34期)；基于贝叶斯压缩感知的噪声MIMO雷达目标成像方法(参见王超宇的《基于贝叶斯压缩感知的噪声MIMO雷达目标成像》，发表在《南京理工大学学报》，2013，第37期)，这些成像算法，将CS理论应用于俯仰向的成像处理，致力于提高成像质量，但对于高分辨率SAR数据量过大的问题，这些方法都未能有效地解决。

由于现有成像方法无法解决高分辨率SAR数据量过大的问题，将会带来以下不足：对于载机平台而言，数据的存储和传输能力都是有限的，SAR数据量过大将会增加载机平台的负荷；对于稀疏性场景，运用压缩感知理论，只需少量的数据就能实现场景的三维成像，若直接存储和传输原始数据将会造成资源的浪费。

发明内容

本发明的目的在于克服上述现有技术中的不足之处，提出一种机载MIMO-SAR下视三维稀疏成像方法。

本发明是通过如下方式实现的：

步骤一：构建MIMO-SAR回波数据稀疏表示模型，并进行降采样处理；

步骤二：依据稀疏模型对数据进行重构，获得高度向与方位向的成像处理结果，并对结果进行俯仰向成像处理，实现场景的三维成像。

所述的步骤一具体包括：

采集的回波数据为T_i表示第i个发射阵元，i＝1,2,K,N_t，共有N_t个发射阵元，阵元间隔为2d；R_j表示第j个发射阵元，j＝1,2,K,N_r，共有N_r个接收阵元；

对回波数据进行成像处理得其中Θ表示高度向与方位向的成像处理结果，Ξ(·)表示高度向与方位向成像算子；

对回波数据进行降采样处理，其中Φ＝{φ_c,d}是观测矩阵，N_r为距离向上的采样点数，且φ_c,d满足：

${\mathrm{\φ}}_{c,d}=\left\{\begin{array}{cc}1,& c=m_d\\ 0,& others\end{array}\right.,d=1,2,\mathrm{...},{\hat{N}}_r;m_d=1,2,\mathrm{...},N_r;.$

所述的步骤二具体包括下述步骤：

Step 1)初始化：

迭代次数计数器k＝0，Θ₀＝0，残差p₀＝s_com，设置最大迭代次数K_max，观测矩阵Φ与迭代门限ξ；

Step 2)计算ΔΘ_k＝Ξ(p_kΦ^H)，并根据公式和δ_k＝|Θ_k+η_kΔΘ_k|_I/η_k计算η_k与δ_k，其中|·|_I表示第I大的值；

Step 3)计算Θ_k+1＝soft{Θ_k+η_kΔΘ_k,δ_k}，soft(x,δ)＝sign(x)·max(|x|-δ,0)；

Step 4)计算s_k+1＝Ξ^-1(Θ_k+1)，并更新残差p_k+1＝s_com-s_k+1Φ；

Step 5)如果k<K_max且|Θ_k+1-Θ_k|₂/|Θ_k|₂>ξ，令k＝k+1并转到Step2)，否则算法结束，输出重构结果Θ_k；

Step 6)对Θ_k进行俯仰向成像处理，实现场景的三维成像。

进一步，所述Step 1)中的最大迭代次数K_max的优选为100-10000，迭代门限的优选为0-0.5。

进一步，所述Step 1)中的最大迭代次数K_max的优选为300-5000，迭代门限的优选为0-0.2。

进一步，所述Step 1)中的最大迭代次数K_max的优选为500-2000，迭代门限的优选为0-0.1。

本发明的有益效果在于：针对现有的机载MIMO-SAR下视三维成像算法无法解决数据量过大的问题，提出的机载MIMO-SAR下视三维稀疏成像方法，利用场景的稀疏性，该方法只需少量的数据即可实现场景的高分辨率三维成像，有效地降低了存储与传输的数据量，减轻载机平台的负荷，减少了资源的浪费。

附图说明

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

图2为MIMO-SAR信号发射接收模型；

图3为散射点模型图；

图4(a)为原始信号的散射点三维像，图4(b)为欠采样比例为1/2的散射点三维像，图4(c)为欠采样比例为1/3的散射点三维像，图4(d)为欠采样比例为1/4的散射点三维像；

图5(a)为高度200m处的原始信号成像结果的剖面图，图5(b)为高度200m处欠采样比例为1/2成像结果的剖面图；

图6(a)为方位0m处原始信号成像结果的剖面图，图6(b)为方位0m处欠采样比例为1/2成像结果的剖面图；

图7(a)为俯仰0m处原始信号成像结果的剖面图，图7(b)为俯仰0m处欠采样比例为1/2成像结果的剖面图；

具体实施方式

下面结合附图和本发明的实例，对本发明作进一步的描述。

如图1所示，本发明通过下列步骤实现：构建回波数据的稀疏表示模型，并对其进行降采样处理；对降采样后的数据进行重构，在重构过程中融入高度向与方位向的成像处理，获得经过高度向与方位向处理后的重构结果，并对重构结果进行俯仰向处理，实现场景的三维成像。具体说明如下：

步骤一：构建回波数据的稀疏表示模型与降采样处理

如图2所示，x轴表示方位向，y轴表示俯仰向，z轴表示高度向。机载下视MIMO-SAR沿俯仰向布置一个线性天线阵列，采用两端阵元发射、中间阵元接收的方式；天线阵列平行于俯仰向，沿方位向以速度v移动，阵列高度为H，目标场景位于阵列正下方；假设天线阵列由N_t个发射天线阵元和N_r个接收天线阵元组成，发射天线阵元用T_i,i＝1,2,K,N_t表示，阵元间隔为2d；接收天线阵元用R_j,j＝1,2,K,N_r表示，阵元间隔为N_td。成像场景中共有N个散射点，其中各个散射点坐标表示为P_n＝(x_n,y_n,z_n),n＝1,2,...,N，此时由发射天线阵元T_i发射线性调频信号：

$s_{T_i}\left(t\right)=rect\left(\frac{t}{T_P}\right)\exp \left(j2\mathrm{\&pi;}\left(f_{c}t+\frac{1}{2}\mathrm{\&mu;}{\hat{t}}^2\right)\right)---\left(1\right)$

其中t为快时间，T_P为脉冲宽度，f_c为起始载频，μ调频斜率，对于成像场景中的N个散射点，接收阵元R_j接收到的信号为：

$\begin{array}{c}{s}_{T_i{R}_j}\left(t\right)=\underset{n=1}{\overset{N}{\&Sigma;}}{\mathrm{\&sigma;}}_{n}rect\left(\frac{t-R_{T_i{R}_j}/c}{T_P}\right)rect\left(\frac{x-x_n}{L_{sar}}\right)\\ \&times;\exp \left(j2\mathrm{\&pi;}\left(f_c\left(t-\frac{R_{T_i{R}_j}}{c}\right)+\frac{1}{2}\mathrm{\&mu;}{\left(t-\frac{R_{T_i{R}_j}}{c}\right)}^2\right)\right)\end{array}---\left(2\right)$

其中σ_n表示第n个散射点的散射系数，表示发射阵元T_i发射信号经过第n个散射点反射后接收阵元R_j接收到的信号所经历的距离，L_sar表示合成孔径长度，c表示电磁波的速度。

对于式(2)，回波信号不具有稀疏性，而对于稀疏场景，高度向与方位向构成的二维平面是稀疏的，因此可将CS理论运用在高度向与方位向的成像处理中；设Θ表示高度向与方位向二维成像处理结果(Θ是稀疏的)，对于信号二维成像处理可表示为：

$\mathrm{\&Theta;}=\mathrm{\&Xi;}\left(s_{T_i{R}_j}\right)---\left(6\right)$

其中Ξ(·)表示二维成像算子，由于无法Θ直接获得，将式(6)转化为：

$s_{T_i{R}_j}={\mathrm{\&Xi;}}^{-1}\left(\mathrm{\&Theta;}\right)---\left(7\right)$

其中Ξ^-1(·)是二维Ξ(·)逆运算，对信号进行降采样处理，可得：

$s_{com}=s_{T_i{R}_j}\mathrm{\&Phi;}={\mathrm{\&Xi;}}^{-1}\left(\mathrm{\&Theta;}\right)\mathrm{\&Phi;}---\left(8\right)$

其中Φ＝{φ_c,d}是观测矩阵，N_r为距离向上的采样点数，观测矩阵中的元素φ_c,d满足：

${\mathrm{\&phi;}}_{c,d}=\left\{\begin{array}{cc}1,& c=m_d\\ 0,& others\end{array}\right.,d=1,2,\mathrm{...},{\hat{N}}_r;m_d=1,2,\mathrm{...},N_r;---\left(9\right)$

其中m_d是个随机数，已知s_com和Φ，为了求解Θ，我们还需得到二维成像算子Ξ(·)与二维成像算子的逆运算Ξ^-1(·)；对于回波数据，在进行高度向成像处理时，先对等效相位中心误差进行补偿，补偿函数为：

$H_{EPC,ij}=\exp \left(j\mathrm{\&pi;}\frac{{(y_{T_i}-y_{R_j})}^2}{2{\mathrm{\&lambda;R}}_s}\right)---\left(10\right)$

其中λ表示波长，雷达的信号发射模式为交替发射，随着载机平台的移动，导致等效阵元在方位向上产生位移，需对其进行补偿。第N_t个PRT的补偿函数为：

$H_{EPC,N_t}=\exp \left(j\mathrm{\&pi;}\frac{{\left(v(N_t-1)PRT\right)}^2}{2{\mathrm{\&lambda;R}}_s}\right)---\left(11\right)$

其中v表示载机速度，经过相位补偿后，对回波数据进行匹配滤波处理，将其变换至快时间频域进行匹配滤波处理，匹配滤波函数为：

$H_t=\exp \left(-j\mathrm{\&pi;}\frac{{f_r}^2}{\mathrm{\&mu;}}\right)---\left(12\right)$

至此完成高度向的成像处理过程，得到：

$s_{T_i{R}_j}\left(f_r\right)=\underset{n=1}{\overset{N}{\&Sigma;}}{\mathrm{\&sigma;}}_{n}rect\left(\frac{f_r}{{\mathrm{\&mu;T}}_P}\right)rect\left(\frac{x-x_n}{L_{sar}}\right)a\left(y_{ij}\right)\exp \left(-j4{\mathrm{\&pi;R}}_{ij}\frac{f_r+f_c}{c}\right)---\left(13\right)$

接下来对方位向进行成像处理。根据泰勒展开公式可得：

$R_{ij}=R_B+\frac{{({\mathrm{vt}}_m-x_n)}^2}{2R_B}+\frac{{y_{ij}}^2-2y_{ij}{y}_n}{2R_B}---\left(14\right)$

其中t_m表示慢时间，将式(14)代入式(15)，并变换到二维时间域得：

$\begin{array}{c}{s}_{T_i{R}_j}\left(f_r,t_m\right)=\underset{n=1}{\overset{N}{\&Sigma;}}{\mathrm{\&sigma;}}_{n}rect\left(\frac{x-x_n}{L_{sar}}\right)a\left(y_{ij}\right)\sin c\left(B_r\left(t-\frac{2R_{ij}}{c}\right)\right)\exp \left(-j2\mathrm{\&pi;}\frac{\left({\mathrm{vt}}_m-x_n\right)}{{\mathrm{\&lambda;R}}_B}\right)\\ \&times;\exp \left(-j2\mathrm{\&pi;}\frac{{y_{ij}}^2-2y_{ij}{y}_n}{{\mathrm{\&lambda;R}}_B}\right)\exp \left(-j\frac{4\mathrm{\&lambda;}}{{\mathrm{\&lambda;R}}_B}\right)\end{array}---\left(15\right)$

其中B_r表示线性调频脉冲的带宽，B_r＝μ·T_P；sinc函数项反应高度信息，第一个反应方位向信息，是方位向成像处理的基础。首先进行方位向的距离徙动校正，校正函数为：

$H_{RCMC1,}=\exp \left(j2\mathrm{\&pi;}\frac{{({\mathrm{vt}}_m-x_n)}^2}{{\mathrm{\&lambda;R}}_B}\right)---\left(16\right)$

徙动校正完成后对其进行方位向的匹配滤波，将回波数据变换至多普勒频域，使用下式进行匹配滤波：

$H_{t_m}=\exp \left(-j\mathrm{\&pi;}\frac{{f_a}^2}{K_a}\right)---\left(17\right)$

其中K_a＝-2v²/λR_s。

总结高度向与方位向成像处理过程可得：

$\mathrm{\&Theta;}=\mathrm{\&Xi;}\left(s_{T_i{R}_j}\right)={\widehat{\mathrm{\&omega;}}}_a\left(\left({\mathrm{\&omega;}}_a\left(\left(\left(\left(\left(s_{T_i{R}_j}\&CenterDot;H_{EPC,ij}\&CenterDot;H_{EPC,N_t}\right)\&CenterDot;{\mathrm{\&omega;}}_r\right)\&CenterDot;H_t\right)\&CenterDot;{\widehat{\mathrm{\&omega;}}}_r\right)\&CenterDot;H_{RCMC,1}\right)\right)\&CenterDot;H_{t_m}\right)---\left(18\right)$

其中(·)表示Hadamard运算，ω_r和表示快时间域的傅里叶变换和傅里叶逆变换，ω_a和表示多普勒域的傅里叶变换和傅里叶逆变换，由式(18)可得其逆运算为：

$s_{T_i{R}_j}={\mathrm{\&Xi;}}^{-1}\left(\mathrm{\&Theta;}\right)=\left(\left(\left(\left(\left({\widehat{\mathrm{\&omega;}}}_a\&CenterDot;\left(\left({\mathrm{\&omega;}}_a\&CenterDot;\mathrm{\&Theta;}\right)\&CenterDot;H_{t_m}^{*}\right)\right)\&CenterDot;H_{RCMC,1}^{*}\right)\&CenterDot;{\mathrm{\&omega;}}_r\right)\&CenterDot;H_t^{*}\right)\&CenterDot;{\widehat{\mathrm{\&omega;}}}_r\right)\&CenterDot;H_{EPC,N_t}^{*}\&CenterDot;H_{EPC,ij}^{*}---\left(19\right)$

其中(·)^*表示共轭运算。

至此，完成了回波数据稀疏表示模型的构建与数据的降采样处理。

步骤二：利用重构算法完成二维成像处理结果的重构，并利用该结果进行三维成像；

本发明对欠采样的数据进行重构，并对其进行三维成像；

具体重构成像算法如下：

输入：压缩数据s_com，降采样矩阵Φ，二维成像算子Ξ(·)与二维成像算子的逆运算Ξ^-1(·)；

输出：重构的二维成像处理结果Θ以及三维成像结果；

Step 1)初始化：

迭代次数计数器k＝0，Θ₀＝0，残差p₀＝s_com，设置最大迭代次数K_max，观测矩阵Φ与迭代门限ξ；

Step 2)计算ΔΘ_k＝Ξ(p_kΦ^H)，并根据公式和δ_k＝|Θ_k+η_kΔΘ_k|_I/η_k计算η_k与δ_k，其中|·|_I表示第I大的值；

Step 3)计算Θ_k+1＝soft{Θ_k+η_kΔΘ_k,δ_k}，soft(x,δ)＝sign(x)·max(|x|-δ,0)；

Step 4)计算s_k+1＝Ξ^-1(Θ_k+1)，并更新残差p_k+1＝s_com-s_k+1Φ；

Step 5)如果k<K_max且|Θ_k+1-Θ_k|₂/|Θ_k|₂>ξ，令k＝k+1并转到Step2)，否则算法结束，输出重构结果Θ_k；

Step 6)对Θ_k进行俯仰向成像处理，实现场景的三维成像。

由此完成对机载MIMO-SAR下视三维稀疏成像。

仿真实验：为了验证本发明所提算法的有效性，我们进行机载MIMO-SAR下视三维稀疏成像计算机仿真实验，雷达发射线性调频信号，采用下视的方式发射信号，数据模拟所需参数设定见表1。

表1仿真参数设置

Table 1 Values of parameters for simulation

仿真1：为了验证算法的有效性，现进行如下仿真实验。采用六个散射点的模型，六个散射点坐标分别为P₁＝(0,0,200)，P₂＝(0,0,180)，P₃＝(15,0,200)，P₄＝(15,0,180)，P₅＝(-15,0,200)，P₆＝(-15,0,180)，散射系数σ_i＝1,i＝1,2,3,4,5,6，具体如图3所示。

对上述散射点模型进行成像处理，成像结果如图4所示，其中图4(a)表示采用原始数据的三维成像结果，成像结果可以清晰地分辨出散射点，采用本发明方法结果如图4(b)、4(c)、4(d)所示，其中图4(b)表示欠采样比例为1/2的三维成像结果；图4(c)表示欠采样比例为1/3的三维成像结果；图4(d)表示欠采样比例为1/4的三维成像结果；成像结果也能清晰地分辨出各散射点，进而验证了算法的有效性。

为了便于分析，给出原始成像方法与欠采样比例1/2的各维度的剖面图。图5给出了高度200m出的俯仰向-方位向的剖面图。图5(a)表示原始成像算法的剖面图，图5(b)表示欠采样比例1/2剖面图，从图中可以清晰地分辨出散射点P₁、P₃、P₅。图6给出了方位0m处的俯仰向-高度向的剖面图，图6(a)表示原始成像算法的剖面图，图6(b)表示欠采样比例1/2的剖面图，从图中可以清晰地分辨出散射点P₃、P₄。图7给出了俯仰0m处的方位向-高度向的剖面图，图7(a)表示原始成像算法的剖面图，图7(b)表示欠采样比例1/2的剖面图，从图中可以清晰地分辨出各散射点，这表明对于稀疏性场景，该方法能够较好地重构数据，利用少量的数据实现三维成像，大幅度的降低了数据量，为解决MIMO-SAR数据量过大问题提供一种有力的方法。

Claims (6)

1.一种机载MIMO-SAR下视三维稀疏成像方法，包括下列步骤：

步骤一：构建MIMO-SAR回波数据稀疏表示模型，并进行降采样处理；

步骤二：依据稀疏模型对数据进行重构，获得高度向与方位向的成像处理结果，并对结果进行俯仰向成像处理，实现场景的三维成像。

2.根据权利要求1所述的机载MIMO-SAR下视三维稀疏成像方法，其特征在于：所述的步骤一具体包括下述步骤：

采集的回波数据为T_i表示第i个发射阵元，i＝1,2,K,N_t，共有N_t个发射阵元，阵元间隔为2d；R_j表示第j个发射阵元，j＝1,2,K,N_r，共有N_r个接收阵元；

对回波数据进行成像处理得其中Θ表示高度向与方位向的成像处理结果，Ξ(·)表示高度向与方位向成像算子；

对回波数据进行降采样处理，其中Φ＝{φ_c,d}是观测矩阵，N_r为距离向上的采样点数，且φ_c,d满足：

${\mathrm{\&phi;}}_{c,d}=\left\{\begin{array}{cc}1,& c=m_d\\ 0,& others\end{array}\right.,d=1,2,\mathrm{...},{\hat{N}}_r;m_d=1,2,\mathrm{...},N_r;.$

3.根据权利要求1所述的机载MIMO-SAR下视三维稀疏成像方法，其特征在于：所述的步骤二具体包括下述步骤：

Step 1)初始化：

迭代次数计数器k＝0，Θ₀＝0，残差p₀＝s_com，设置最大迭代次数K_max，观测矩阵Φ与迭代门限ξ；

Step 2)计算ΔΘ_k＝Ξ(p_kΦ^H)，并根据公式和δ_k＝|Θ_k+η_kΔΘ_k|_I/η_k计算η_k与δ_k，其中|·|_I表示第I大的值；

Step 3)计算Θ_k+1＝soft{Θ_k+η_kΔΘ_k,δ_k}，soft(x,δ)＝sign(x)·max(|x|-δ,0)；

Step 4)计算s_k+1＝Ξ^-1(Θ_k+1)，并更新残差p_k+1＝s_com-s_k+1Φ；

Step 5)如果k<K_max且|Θ_k+1-Θ_k|₂/|Θ_k|₂>ξ，令k＝k+1并转到Step2)，否则算法结束，输出重构结果Θ_k；

Step 6)对Θ_k进行俯仰向成像处理，实现场景的三维成像。

4.根据权利要求3所述的机载MIMO-SAR下视三维稀疏成像方法，其特征在于：所述Step1)中的最大迭代次数K_max的优选为100-10000，迭代门限的优选为0-0.5。

5.根据权利要求4所述的机载MIMO-SAR下视三维稀疏成像方法，其特征在于：所述Step1)中的最大迭代次数K_max的优选为300-5000，迭代门限的优选为0-0.2。

6.根据权利要求5所述的机载MIMO-SAR下视三维稀疏成像方法，其特征在于：所述Step1)中的最大迭代次数K_max的优选为500-2000，迭代门限的优选为0-0.1。