CN105068071A - 一种基于反投影算子的快速成像方法 - Google Patents

Abstract

本发明提供了一种基于反投影算子的快速成像方法。该快速成像方法包括：步骤A：根据观测场景中有效非零像素点总数目的估计值，设定阈值μ的取值；步骤B：令观测场景后向散射系数估计值的初值，残差向量的初值其中，s_u为原始回波数据；步骤C：在复数域近似信息传递算法的算法框架下，利用基于反投影成像原理的距离方位解耦算子，通过循环迭代的方式计算观测场景的后向散射系数。本发明在保证与原稀疏重建算法拥有相同的成像质量前提下，基于反投影算子的快速成像方法具有较低的时间复杂度与空间复杂度。

Description

一种基于反投影算子的快速成像方法

技术领域

本发明涉及稀疏微波成像技术领域，尤其涉及一种基于反投影算子的快速成像方法。

背景技术

与传统微波成像方法相比，稀疏微波成像方法在降低成像所需数据采样率、抑制旁瓣与加性噪声、提升图像分辨率等方面具有较为明显的优势。而作为实现稀疏微波成像方法的主要途径，近些年来稀疏重建算法的改进与发展已得到了电子信息工程领域科研人员的广泛关注。从目前来看，在现有的稀疏重建算法中，大多存在与观测矩阵相关的矩阵向量运算。当利用此类算法处理从空间尺度较大的观测场景返回的原始回波数据时，便会产生巨大的计算负担与存储损耗，从而导致常规计算机很难满足微波成像所需的系统硬件要求。

目前，解决该问题的方法主要分为两大类，即原始回波数据预处理方法和基于传统成像技术的近似观测方法。原始回波数据预处理方法通过距离向脉冲压缩，使预处理后的回波数据中每个距离门之间相互独立，然后，根据每个距离门内的回波数据，单独构建低维度的观测矩阵，最后，再利用各种稀疏重建算法实现微波成像。尽管此类方法能够有效地重建出空间尺度较大的观测场景，并同时实现方位向所需采样率的降低，但是它们却无法利用回波数据在距离向上的冗余信息。为此，微波遥感与雷达成像研究领域的科研人员又进一步提出了基于传统成像技术的近似观测方法。近似观测方法的核心思想是，通过将原始回波数据在距离向和方位向上进行解耦，利用传统成像技术中的基本算法运算(快速傅里叶变换/逆快速傅里叶变换(FastFourierTransform/InverseFastFourierTransform，FFT/IFFT)、相位复乘、插值操作等)，替代在原稀疏重建算法中由与观测矩阵相关的矩阵向量运算所实现的功能，从而达到在保证成像质量的前提下显著地降低算法时间复杂度与空间复杂度的目的。

因为FFT算法的时间复杂度与空间复杂度分别为线性对数阶和线性阶，所以目前在近似观测方法中的距离方位解耦算子均是根据基于FFT的传统成像方法的基本原理推导而来的。虽然上述方法能够极大地降低稀疏重建算法的时间复杂度与空间复杂度，但是受制于基于FFT的传统成像方法推导过程中所用到的观测几何近似与特定应用假设，导致它们只适用于相应的雷达观测模式与平台运动方式，使其应用范围受到了很大的局限。相比基于FFT的传统成像方法，反投影成像方法的应用范围更为广泛。它不但可以对平台运动的任意轨迹进行精确的运动补偿，还可以对在任意带宽与合成孔径角度下获取的回波数据进行准确的聚焦，同时也具有并行性好、成像质量高等诸多优点。鉴于上述内容，依据反投影成像原理与逆成像回波仿真思想，我们推导出了基于反投影成像原理的距离方位解耦算子。该算子不仅可以降低稀疏重建算法的时间复杂度与空间复杂度，而且继承了反投影成像方法的全部优点，具有更为广泛的应用范围。

复数域近似信息传递算法(ComplexApproximateMessagePassingAlgorithm，CAMP)是根据图模型理论推导而来的高效稀疏重建算法，它具有收敛速度快、重建精度高、恢复信号所需采样率低等诸多优点。为了进一步提升稀疏微波成像方法的运算效率与成像性能，我们将基于反投影成像原理的距离方位解耦算子与复数域近似信息传递算法相结合，提出了一种基于反投影算子的快速成像方法(BP-CAMP)。理论分析与成像实验均已证明，基于反投影算子的快速成像方法具有较低的时间复杂度与空间复杂度，能够在低于奈奎斯特采样率的条件下完成高质量微波成像，同时相比基于其他解耦算子的稀疏重建算法，具有更为普遍的适用性。

发明内容

(一)要解决的技术问题

在现有的稀疏重建算法中，大多存在与观测矩阵相关的矩阵向量运算。当利用此类算法处理从空间尺度较大的观测场景返回的原始回波数据时，便会产生巨大的计算负担与存储损耗，从而导致常规计算机很难满足微波成像所需的系统硬件要求。尽管可以利用由基于FFT的传统成像方法推导的距离方位解耦算子替代上述矩阵向量运算，来降低稀疏重建算法的时间复杂度与空间复杂度，但是受制于基于FFT的传统成像方法推导过程中所用到的观测几何近似与特定应用假设，导致它们只适用于相应的雷达观测模式与平台运动方式，使其应用范围受到了很大的局限。鉴于上述技术问题，本发明提供了一种基于反投影算子的快速成像方法，以降低原稀疏重建算法的时间复杂度与空间复杂度，进而降低实现稀疏微波成像所需的计算负担和存储损耗。

(二)技术方案

根据本发明的一个方面，提供了一种基于反投影算子的快速成像方法。该快速成像方法包括：步骤A：根据观测场景中有效非零像素点总数目的估计值，设定阈值μ的取值；步骤B：令观测场景后向散射系数估计值的初值残差向量的初值其中，s_u为原始回波数据；步骤C：在复数域近似信息传递算法的算法框架下，利用基于反投影成像原理的距离方位解耦算子，通过循环迭代的方式计算观测场景的后向散射系数

(三)有益效果

从上述技术方案可以看出，本发明基于反投影算子的快速成像方法具有以下有益效果：

(1)在保证与原稀疏重建算法拥有相同的成像质量前提下，基于反投影算子的快速成像方法具有较低的时间复杂度与空间复杂度；

(2)当数据采样率低于奈奎斯特采样定理要求时，基于反投影算子的快速成像方法仍能完成高质量的微波成像；

(3)与由基于FFT的传统成像方法推导的距离方位解耦算子相比，基于反投影成像原理的距离方位解耦算子具有更为广泛的应用范围，且基于反投影算子的快速成像方法在微波成像领域也具有更为普遍的适用性。

附图说明

图1为基于反投影成像原理的距离方位解耦算子的原理框图；

图2为根据本发明实施例的基于反投影算子的快速成像方法流程图；

图3(a)～图3(d)为利用不同的成像方法处理实际数据得到的成像结果。

具体实施方式

本发明根据反投影成像原理与逆成像回波仿真思想，推导出基于反投影成像原理的距离方位解耦算子，通过将该算子引入到复数域近似信息传递算法之中，构成了基于反投影算子的快速成像方法。

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明做进一步详细说明。

图1为基于反投影成像原理的距离方位解耦算子的原理框图。如图1所示，基于反投影成像原理的距离方位解耦算子主要由微波成像过程与回波生成过程两个部分组成。

根据反投影成像原理，微波成像过程可以分为以下三个步骤，即距离向脉冲压缩、空间域sinc插值和子图像沿方位向相干累加。据此，可以获得基于反投影成像原理的成像算子的表达式为：

该成像算子用于将回波域数据s转变为与其对应的图像域数据。其中，⊙表示哈达马乘法运算；表示把按采样时间排序的一维向量重新排列成按距离方位采样获得的二维矩阵的功能算子，而则表示其逆过程；表示距离向脉冲压缩算子；表示空间域sinc插值算子；表示子图像沿方位向相干累加算子；ψ₁和分别表示sinc插值前和sinc插值后的附加相位补偿矢量。

假设观测场景在一定的合成孔径角度范围内具有各向同性的性质，那么，可以将的逆过程认为是雷达在不同方位向采样位置观测到的是同一幅微波图像。而依据距离压缩与sinc插值的基本原理，可以很容易地明白两者均是可逆的。由此根据逆成像回波仿真思想，便可获得基于反投影成像原理的逆成像算子的表达式为：

该逆成像算子用于将图像域数据r转变为与其对应的回波域数据。其中，上角标*表示取向量元素的共轭，上角标-1表示相应算子的逆过程。

综上可知，(1)式和(2)式表示的微波成像过程与回波生成过程共同构成了基于反投影成像原理的距离方位解耦算子。

下面我们以雷达发射信号为步进频信号作为例子，详细地说明基于反投影成像原理的距离方位解耦算子的基本原理。为方便解释说明，我们已将功能算子和所表示的向量与矩阵间的转换关系表示在相应元素的标号序列中。

步进频雷达系统在接收端获得的解调后基带信号可以表示为：

其中，f_k＝f₀+kΔf，k＝0，…，K-1为第k个脉冲的载波频率，K为距离向采样点总数目，Δf为频率的步进长度；R(u_n，v_n，θ_m)表示雷达方位向采样位置θ_m与观测场景中目标点坐标(u_n，v_n)之间的距离；r(u_n，v_n)表示目标点的后向散射系数；c为电磁波在自由空间的传播速度。

首先，利用逆离散傅里叶变换对解调后的基带信号进行距离向脉冲压缩操作，可以得到：

其中，为雷达在方位向采样位置θ_m处获得的回波数据；M为方位向采样点总数目；在这里表示距离向上的逆离散傅里叶变换。

向量中的元素可以表示为：

其中，χ＝l-2KΔfR(u_n，v_n，θ_m)/c，l＝O，…，K-1；t_l为与标号l对应的采样时刻；f_c为步进频信号的中心频率。

然后，利用sinc插值方法计算出观测区域内网格点上的信号样值，即

其中，为在sinc插值前的附加相位补偿矢量，它的元素可以表示为ψ₁(l)＝exp(-jπl(K-1)/K)。

向量中的元素可以表示为：

最后，将低分辨率微波图像沿方位向相干累加，获得最终的成像结果

其中，为在sinc插值后的附加相位补偿矢量，它的元素可以表示为

与上述根据微波成像原理获得成像算子的过程类似，根据(2)式所示的回波生成原理，同理可得在接收到的回波数据为步进频信号时的逆成像算子，其具体过程可以表示为：

首先，根据算子所要表达的含义可知，雷达系统在不同方位向采样位置θ_m处观测到的是同一幅微波图像即

然后，利用在上述微波成像过程中所用sinc插值方法的逆过程，便可得到

其中，为附加相位补偿矢量的共轭。

最后，利用离散傅里叶变换将时域信号变换为频域信号，即

便可获得雷达系统在方位向采样位置θ_m处接收到的步进频信号的近似值其中，为附加相位补偿矢量ψ₁的共轭；在这里表示距离向上的离散傅里叶变换。

以上采用步进频信号为例，说明基于反投影成像原理的距离方位解耦算子的具体操作流程。基于上述解释说明，本领域的技术人员可以清楚地了解与其他类型的雷达发射信号(例如：线性调频信号、高斯随机噪声信号、非线性调频信号、相位编码信号、时间-频率编码信号等)相对应的距离方位解耦算子的具体操作流程，此处不再一一说明。

相比其他常用的稀疏重建算法，复数域近似信息传递算法具有收敛速度快、重建精度高、恢复信号所需采样率低等诸多优点。为了进一步提升稀疏微波成像方法的运算效率与成像性能，我们将基于反投影成像原理的距离方位解耦算子引入到复数域近似信息传递算法之中，构成了基于反投影算子的快速成像方法。

本发明的一个示例性实施例中，提供了一种基于反投影算子的快速成像方法。图2为根据本发明实施例的基于反投影算子的快速成像方法流程图。如图2所示，本实施例基于反投影算子的快速成像方法包括：

步骤A：根据观测场景中有效非零像素点总数目的估计值，来设定阈值μ＞0的取值；

有效非零像素点总数目的估计值越大，阈值μ的取值就越小，反之，阈值μ的取值就越大。在一般情况下，阈值μ的取值为大于0的正实数。本实施例中，观测场景中像素点总数目为2.5×10⁵，其中有效非零像素点总数目的估计值为1×10⁵，该阈值μ取10。

步骤B：变量初始化，令观测场景后向散射系数估计值的初值残差向量的初值其中，s_u为雷达系统接收到的原始回波数据；

步骤C：在复数域近似信息传递算法的算法框架下，利用基于反投影成像原理的距离方位解耦算子，通过循环迭代的方式计算观测场景的后向散射系数

在本步骤中，设迭代变量为t，最大迭代次数为maxiter。其中，最大迭代次数maxiter由所需的收敛精度决定。一般情况下，最大迭代次数maxiter的取值介于30～1000之间。本实施例中，收敛精度设定为10^-6，最大迭代次数设定为50。

每一次循环迭代的过程包括：

子步骤C1：利用下式计算中间变量

其中，表示重采样过程；为基于反投影成像原理的成像算子。

子步骤C2：利用下式计算的标准差估计

其中，median(·)表示中值滤波过程。

子步骤C3：利用下式计算中间变量γ^(t)；

其中，δ＞0为降采样率；〈·〉表示求解向量元素均值的运算；η(x+jy；τ)表示以复数x+jy为输入、以τ＞0为阈值的软阈值函数；η^R和η^I分别为软阈值函数的实部和虚部；和分别为η^R关于复输入实部x的偏导数和η^I关于复输入虚部y的偏导数。

子步骤C4：利用下式计算残差向量

其中，表示降采样过程；为基于反投影成像原理的逆成像算子。

子步骤C5：利用下式计算观测场景后向散射系数的估计

至此，便实现了基于本发明的稀疏微波成像。

在实际应用中，本实施例可以按照以下方式执行：

1、输入：原始回波数据s_u，基于反投影成像原理的距离方位解耦算子和降采样算子和重采样算子指定的阈值μ。

2、初始化：

3、循环迭代：

fort＝1tomaxiter

endfor

4、输出：观测场景的后向散射系数

根据上述具体实施方式，下面将利用实际数据成像实验，对本实施例所述方法的应用优势进行验证。在实际数据成像实验中，数据采集所用雷达的系统参数如表1所示。而用来完成微波成像的不同方法是在硬件配置为8核2.4GHz中央处理器、32GB内存的服务器上，利用MATLABR2013a软件实现的。值得注意的是，在实际数据处理过程中，所用降采样数据是在距离向和方位向的降采样率均为70％的条件下，从满采样数据中随机抽取获得的，同时，基于反投影算子的快速成像方法和复数域近似信息传递方法在实验中的迭代次数均设为50次。

实际数据成像实验的处理结果如表2和图3所示。表2列出了基于反投影算子的快速成像方法和复数域近似信息传递方法的程序运行时间与所需内存大小的比较结果。如表2所示，由于复数域近似信息传递方法需要进行与观测矩阵相关的矩阵向量运算，因此它的程序运行时间和所需内存大小要多于基于反投影算子的快速成像方法。图3(a)～图3(d)为利用不同的成像方法处理实际数据得到的成像结果。其中，图3(a)是反投影成像方法在满采样条件下获得的成像结果；图3(b)是反投影成像方法在使用49％满采样数据条件下获得的成像结果；图3(c)是基于反投影算子的快速成像方法在使用49％满采样数据条件下获得的成像结果；图3(d)是复数域近似信息传递方法在使用49％满采样数据条件下获得的成像结果。通过比较图3(a)～图3(d)所示的不同方法的成像结果，可知当数据采样率低于奈奎斯特采样定理要求时，基于反投影算子的快速成像方法和复数域近似信息传递方法均能实现微波成像，且成像质量要明显好于后投影成像方法的成像结果(图3(b))。

因为实验所用的实际数据是由圆轨步进频合成孔径雷达系统获取的回波信号，所以可以验证：与由基于FFT的传统成像方法推导的距离方位解耦算子相比，基于反投影成像原理的距离方位解耦算子具有更为广泛的应用范围，且基于反投影算子的快速成像方法在微波成像领域也具有更为普遍的适用性。

表1为圆轨步进频合成孔径雷达系统的系统参数

表2为不同成像方法间程序运行时间与所需内存大小的比较结果

成像方法	程序运行时间	所需内存大小
			BP-CAMP	150s	20MB

CAMP

395s

9GB

至此，已经结合附图对本发明实施例进行了详细描述。依据以上描述，本领域技术人员应当对本发明基于反投影算子的快速成像方法有了清楚的认识。

需要说明的是，在附图或说明书正文中，未绘示或描述的实现方式，均为所属技术领域中普通技术人员所知的形式，并未进行详细说明。此外，上述对各元件和方法的定义并不仅限于实施例中提到的各种具体结构、形状或方式，本领域普通技术人员可对其进行简单地更改或替换，例如：本文可提供包含特定值的参数的示范，但这些参数无需确切等于相应的值，而是可在可接受的误差容限或设计约束内近似于相应值。

综上所述，本发明依据反投影成像原理与逆成像回波仿真思想，推导出基于反投影成像原理的距离方位解耦算子，通过将该算子引入到复数域近似信息传递方法之中，构成了基于反投影算子的快速成像方法。本发明基于反投影算子的快速成像方法不但具有较低的时间复杂度与空间复杂度，而且还继承了反投影成像方法的全部优点，在微波成像领域具有更为普遍的适用性。

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

Claims (10)

1.一种基于反投影算子的快速成像方法，其特征在于，包括：

步骤A：根据观测场景中有效非零像素点总数目的估计值，设定阈值μ的取值；

步骤B：令观测场景后向散射系数估计值的初值残差向量的初值z⁽⁰⁾＝s_u，其中，s_u为原始回波数据；

步骤C：在复数域近似信息传递算法的算法框架下，利用基于反投影成像原理的距离方位解耦算子，通过循环迭代的方式计算观测场景的后向散射系数

2.根据权利要求1所述的快速成像方法，其特征在于，所述步骤C中，设定迭代变量为t，最大迭代次数为maxiter，每次迭代过程包括：

子步骤C1：利用下式计算中间变量

其中，表示重采样过程；为基于反投影成像原理的成像算子；

子步骤C2：利用下式计算的标准差估计

${\widehat{\mathrm{\σ}}}_t=median\left(\right|{\tilde{r}}^{\left(t\right)}\left|\right)/\sqrt{\ln 2}$

其中，median(·)表示中值滤波过程；

子步骤C3：利用下式计算中间变量γ^(t)；

${\mathrm{\&gamma;}}^{\left(t\right)}=\frac{1}{2\mathrm{\&delta;}}(<\frac{\&part;{\mathrm{\&eta;}}^R}{\&part;x}({\tilde{r}}^{\left(t\right)},\mathrm{\&mu;}{\widehat{\mathrm{\&sigma;}}}_t)>+<\frac{\&part;{\mathrm{\&eta;}}^I}{\&part;y}({\tilde{r}}^{\left(t\right)},\mathrm{\&mu;}{\widehat{\mathrm{\&sigma;}}}_t)>)$

其中，δ＞0为降采样率；<·>表示求解向量元素均值的运算；η(x+jy；τ)表示以复数x+jy为输入、以τ＞0为阈值的软阈值函数；η^R和η^I分别为软阈值函数的实部和虚部；和分别为η^R关于复输入实部x的偏导数和η^I关于复输入虚部y的偏导数；

子步骤C4：利用下式计算残差向量z^(t)；

其中，表示降采样过程；为基于反投影成像原理的逆成像算子；

子步骤C5：利用下式计算观测场景后向散射系数的估计

${\hat{r}}^{\left(t\right)}=\mathrm{\&eta;}({\tilde{r}}^{\left(t\right)};\mathrm{\&mu;}{\widehat{\mathrm{\&sigma;}}}_t)$

至此，实现了基于反投影算子的快速成像。

3.根据权利要求2所述的快速成像方法，其特征在于，所述子步骤C1中，基于反投影成像原理的成像算子的表达式为：

该成像算子用于将回波域数据s转变为与其对应的图像域数据，其中，⊙表示哈达马乘法运算；表示把按采样时间排序的一维向量重新排列成按距离方位采样获得的二维矩阵的功能算子，而则表示其逆过程；表示距离向脉冲压缩算子；表示空间域sinc插值算子；表示子图像沿方位向相干累加算子；ψ₁和分别表示sinc插值前和sinc插值后的附加相位补偿矢量。

4.根据权利要求3所述的快速成像方法，其特征在于，所述子步骤C4中，基于反投影成像原理的逆成像算子的表达式为：

该逆成像算子用于将图像域数据r转变为与其对应的回波域数据，其中，上角标*表示取向量元素的共轭，上角标-1表示相应算子的逆过程。

5.根据权利要求4所述的快速成像方法，其特征在于，所述原始回波数据为通过以下雷达发射信号中的一种而获得的回波数据：步进频信号、线性调频信号、非线性调频信号、相位编码信号、时间-频率编码信号和高斯随机噪声信号。

6.根据权利要求5所述的快速成像方法，其特征在于，所述雷达发射信号为步进频信号，在接收端解调后的基带信号为：

其中，f_k＝f₀+kΔf，k＝0，…，K-1为第k个脉冲的载波频率，K为距离向采样点总数目，Δf为频率的步进长度；R(u_n，v_n，θ_m)表示雷达方位向采样位置θ_m与观测场景中目标点坐标(u_n，v_n)之间的距离；r(u_n，v_n)表示目标点的后向散射系数；c为电磁波在自由空间的传播速度。

7.根据权利要求6所述的快速成像方法，其特征在于，所述基于反投影成像原理的成像算子的计算过程包括：

首先，利用逆离散傅里叶变换对解调后的基带信号进行距离向脉冲压缩操作，得到：

其中，为雷达在方位向采样位置θ_m处获得的回波数据；M为方位向采样点总数目；在这里表示距离向上的逆离散傅里叶变换；

向量中的元素表示为：

其中，χ＝l-2KΔfR(u_n，v_n，θ_m)/c，l＝0，…，K-1；t_l为与标号l对应的采样时刻；f_c为步进频信号的中心频率；

然后，利用sinc插值方法计算观测区域内网格点上的信号样值：

其中，为在sinc插值前的附加相位补偿矢量，它的元素表示为ψ₁(l)＝exp(-jπl(K-1)/K)；

向量中的元素可以表示为：

$s_1(u_n,v_n,{\mathrm{\&theta;}}_m)=r(u_n,v_n)\exp (-\frac{j4{\mathrm{\&pi;}}_{c}R(u_n,v_n,{\mathrm{\&theta;}}_m)}{c})$

最后，将低分辨率微波图像沿方位向相干累加，获得最终的成像结果：

其中，为在sinc插值后的附加相位补偿矢量，它的元素表示为 ${\mathrm{\&psi;}}_{2,{\mathrm{\&theta;}}_m}(u_n,v_n)=\exp \left(j4{\mathrm{\&pi;f}}_{c}R\right(u_n,v_n,{\mathrm{\&theta;}}_m)/c).$

8.根据权利要求7所述的快速成像方法，其特征在于，所述基于反投影成像原理的逆成像算子的计算过程包括：

首先，确定雷达系统在不同方位向采样位置θ_m处观测到的是同一幅微波图像即：

然后，利用sinc插值方法的逆过程，得到

其中，为附加相位补偿矢量的共轭；

最后，利用离散傅里叶变换将时域信号变换为频域信号，即

获得雷达系统在方位向采样位置θ_m处接收到的步进频信号的近似值其中，为附加相位补偿矢量ψ₁的共轭；在这里表示距离向上的离散傅里叶变换。

9.根据权利要求1至8中任一项所述的快速成像方法，其特征在于，所述步骤C中，最大迭代次数maxiter的取值由所需收敛精度决定，其取值范围介于30～1000之间。

10.根据权利要求1至8中任一项所述的快速成像方法，其特征在于，所述步骤A中，阈值μ的取值为大于0的正实数。