CN108226927B - 基于加权迭代最小稀疏贝叶斯重构算法的sar成像方法 - Google Patents

Abstract

本发明公开一种基于加权迭代最小稀疏贝叶斯重构算法的SAR成像方法，它是针对线阵SAR观测场景目标空间中主散射目标在空间上稀疏的特征，通过建立线阵SAR原始回波信号与观测场景目标空间中散射系数的线性测量矩阵，本发明在迭代最小稀疏贝叶斯重构(SBRIM)算法的基础上，对代价函数中L1范数进行加权，对距离向进行脉冲压缩、划分等距离面，然后再对每一个等距离的二维平面进行稀疏重构。与传统方法相比，本发明具有重构精度高、运算效率高的特点，本发明可适用于合成孔径雷达成像和地球遥感等领域。

Description

基于加权迭代最小稀疏贝叶斯重构算法的SAR成像方法

技术领域

本技术发明属于雷达技术领域，它特别涉及了合成孔径雷达(SAR)成像技术领域。

背景技术

作为一种工作在微波波段的有源雷达，合成孔径雷达(Synthetic ApertureRadar,SAR)具有全天时、全天候的成像能力，即无论是白天或黑夜、晴天还是雷雨风雪天气，都可以随时随地成像，克服了光学和红外系统不能在晚上和复杂天气条件进行成像的缺点。传统的SAR成像一般只具有二维分辨率，在一些起伏比较大的地方比如陡峭的山峰、峡谷以及城市中矗立挺拔的高楼时，传统SAR成像存在的失真(阴影遮挡效应、空间模糊、顶底倒置等)导致空间的一些重要信息(比如高度)丢失，所以能对目标进行三维成像是非常有必要的，为了适应这种需求，目前常见的三维成像技术有圆周SAR(Circular SAR)三维成像、层析SAR(Tomography SAR)三维成像、线阵SAR(LASAR)三维成像。

圆周SAR三维成像的基本思想是载有雷达的平台沿着空间一定的轨迹做圆周运动形成一个曲面线阵得到这个平面的二维分辨率，再结合脉压技术获得距离向的分辨率，实现三维成像。而层析SAR三维成像技术是在干涉SAR技术的基础上发展而来的。它主要是利用多次平行航过的二维数据获得在层析向的第三维分辨率。但是由于圆周SAR三维成像需要雷达平台做圆周运动，层析SAR三维成像不仅要求雷达平台多次航过，而且要求多次航过的轨迹必须平行。这不仅对实际的雷达成像实验中增加了难度，而且增加了飞行的成本，从而限制了它们的发展。

线阵SAR三维成像的基本原理是在切行迹向添加线阵天线，通过沿航迹向平台的飞行形成虚拟的面阵进而获得二维分辨率，距离向再通过脉压技术获得第三维的分辨率，从而具有三维的成像能力。相比于圆周SAR三维成像，线阵SAR三维成像不需要圆周运动的轨迹；相比于层析SAR三维成像需要航过多次，线阵SAR三维成像只需一次航过，所以线阵SAR三维成像相对于层析SAR和圆周SAR三维成像有更强的灵活性。目前线阵SAR三维成像技术在地形测绘、城市测绘、灾难救援、军事探测等领域发挥着中重要的作用。

传统基于匹配滤波的SAR成像方法的分辨率受到限制，具体来说就是距离向的分辨率受信号带宽的影响，沿航迹分辨率受合成孔径长度的影响，切航迹的分辨率受线阵天线的影响。尤其是切航迹的分辨率，如果按照传统的方法很难提高。如果一个信号是稀疏的或者是可压缩的，那么这个信号就能以低于Nyquist采样定理要求的采样率精确的重构出该信号，这就是压缩感知的基本思想。针对压缩感知理论用于SAR成像，目前的重构算法大概可以分为以下几类：贪婪追踪算法、凸松弛算法、贝叶斯框架算法、组合算法。

基于稀疏贝叶斯理论的算法，它的主要原理是对测量信号进行合理的先验概率建模，然后利用构造的似然函数进行相关的参数估计。详见“Wei S J,Zhang X L,ShiJ.Sparse reconstruction for linear array SAR 3-D imaging based on Bayesianestimation[C]//Radar(Radar),2011 IEEE CIE International Conference on.IEEE,2011,2:1522-1525”。

基于贝叶斯理论的算法的一个优点是可以针对不同的信号进行不同的先验概率建模，具体问题，具体建模。从而该类算法具有灵活性强、重构精度高的特点。迭代最小稀疏贝叶斯重构(SBRIM)算法相对于贝叶斯压缩感知(BCS)算法而言具有参数选择简单的特点。(详见Wei S J,Zhang X L,Shi J.Sparse reconstruction for linear array SAR 3-Dimaging based on Bayesian estimation[C]//Radar(Radar),2011 IEEE CIEInternational Conference on.IEEE,2011,2:1522-1525)。基于迭代最小稀疏贝叶斯理论的重构算法，本文提出一种基于加权迭代最小稀疏贝叶斯重构(Weighted SparsityBayesian Recovery via Iterative Minimum,WSBRIM)算法。

发明内容

为了提高算法重构精度，本发明提出一种基于加权迭代最小稀疏贝叶斯重构算法的SAR成像方法，该方法在迭代最小稀疏贝叶斯重构(SBRIM)算法的基础上，对1范数进行加权，先对距离向进行脉冲压缩、划分等距离面，然后再对每一个等距离的二维平面进行稀疏重构，有效的提高了运算效率和重构精度。

为了方便描述本发明的内容，首先作以下术语定义：

定义1、合成孔径雷达(SAR)

合成孔径雷达是将雷达是将雷达固定于载荷运动平台上，结合运动平台的运动以合成线性阵列以达到运动向的分辨率，再利用雷达波束向回波延时实现距离一维成像，从而实现观测目标二维成像的一种合成孔径雷达技术。

定义2、标准合成孔径雷达回波数据距离向脉冲压缩

标准合成孔径雷达回波数据距离向脉冲压缩是指利用合成孔径雷达发射信号参数，采用匹配滤波技术对合成孔径雷达的距离向信号进行信号聚焦成像过程。详见文献“雷达成像技术”，保铮，邢孟道，王彤，电子工业出版社，2005。

定义3、范数

设X是数域

上线性空间，其中

表示复数域，若它满足如下性质：||X||≥0，且||X||＝0仅有X＝0；||aX||＝|a|||X||，a为任意常数；||X₁+X₂||≤||X₁||+||X₂||，则称||X||为X空间上的范数(norm)，其中X₁和X₂为X空间上的任意两个值。对于定义1中的N×1维离散信号向量X＝[x₁,x₂,…,x_N]^T，向量X的LP范数表达式为

其中x_i为向量X的第i个元素，∑|·|表示绝对值求和运算符号，向量X的L1范数表达式为

向量X的L2范数表达式为

向量X的L0范数表达式为

且x_i≠0。详见文献“矩阵理论”，黄廷祝等编著，高等教育出版社出版。

定义4、方位向、距离向

将雷达平台运动的方向叫做方位向，将垂直于方位向的方向叫做距离向。

定义5、压缩感知稀疏重构理论

如果一个信号是稀疏的或可压缩的，那么该信号就可以用远低于奈奎斯特采样定理所要求的采样率来无失真的重构出该信号。如果信号稀疏，并且测量矩阵满足不相干和RIP属性，使用压缩感知恢复的信号稀疏重建可以通过解决以下最优方程来实现：

α是估计信号，y是测量信号，Θ是测量矩阵，ε是噪声门限。详见文献“线阵三维合成孔径雷达稀疏成像技术研究”韦顺军，2013。

定义6、迭代最小稀疏贝叶斯重构(SBRIM)算法

迭代最小稀疏贝叶斯重构(Sparsity Bayesian Recovery via IterativeMinimum,SBRIM)是由电子科技大学的韦顺军在2011年提出的一个算法,详见Wei S J,Zhang X L,Shi J.Sparse reconstruction for linear array SAR 3-D imaging basedon Bayesian estimation[C]//Radar(Radar),2011 IEEE CIE InternationalConference on.IEEE,2011,2:1522-1525

定义7、合成孔径雷达原始回波仿真方法

合成孔径雷达原始回波仿真方法是指基于合成孔径雷达成像原理仿真出一定系统参数条件下具有合成孔径雷达回波信号特性的原始信号的方法，详细内容可参考文献：“张朋，合成孔径雷达回波信号仿真研究，西北工业大学博士论文，2004”。

定义8、线阵SAR的快时刻和慢时刻

线阵SAR运动平台飞过一个方位向合成孔径长度所需要的时间称为慢时间，雷达系统以一定时间长度的重复周期发射接收脉冲，因此慢时间可以表示为一个以脉冲重复周期为步长的离散化时间变量，其中每一个脉冲重复周期离散时间变量值为一个慢时刻。快时刻是指在一个脉冲重复周期内，距离向采样回波信号的时间间隔变量。详见文献“合成孔径雷达成像原理”，皮一鸣等编著，电子科技大学出版社出版。

定义9、信号线性测量模型

对于一个数字信号测量系统，假设N×1维离散信号向量X＝[x₁,x₂,…,x_N]^T为该测量系统需要测量的信号，向量Y＝[y₁,y₂,…,y_M]^T为该测量系统输出的维离散信号向量，其中右上角T为转置运算符号，y₁为向量Y中的第一个元素，y₂表示向量Y中的第二个元素，y_M表示向量Y中的第M个元素，信号的线性测量模型是指测量信号Y和被测量信号X的关系可以表示为Y＝AX，其中A为M×N矩阵，矩阵A成为测量系统中信号X的测量矩阵。

定义10传统理论成像分辨率

线阵SAR成像传统理论分辨率是指利用经典匹配滤波理论成像算法得到线阵SAR系统在距离向、方位向和切航迹向的成像分辨率。对于收发共用天线，线阵SAR距离向的分辨率记为ρ_r，近似表达式为

其中C为电磁波在空气中传播的速度，B_r为线阵SAR发射信号带宽；方位向分辨率记为ρ_a，近似表达为

其中D_a为天线在方位向的真实孔径；切航迹向的分辨率记为

其中λ为线阵SAR雷达载波波长，R₀为线阵SAR平台到观测场景中心的参考斜距，L为线阵天线长度。

本发明提出的一种基于加权迭代最小稀疏贝叶斯重构算法的SAR成像方法，它包括如下步骤：

步骤1、初始化SAR系统参数：

初始化SAR系统参数包括：平台速度矢量记为

线阵天线各阵元初始位置矢量，记做

其中n为天线各阵元序号为自然数，n＝1,2,...,N，N为线阵天线的阵元总数，线阵天线长度，记做L；雷达发射信号载频f_c；雷达发射信号的调频斜率f_dr；脉冲重复时间记为PRI；雷达系统的脉冲重复频率PRF；雷达发射信号带宽B_r；电磁波在空气中的传播速度C；距离向快时间记做t，t＝1,2...T，T为距离向快时刻总数，方位向慢时刻记做l，l＝1,2,...K，K为方位向慢时刻总数；上述参数均为SAR系统标准参数，其中雷达发射信号载频f_c，雷达发射信号的调频斜率f_dr，脉冲重复时间PRI，雷达系统的脉冲重复频率PRF，雷达发射信号带宽B_r，线阵天线的阵元总数N，线阵天线长度L在线阵SAR系统设计过程中已经确定；平台速度矢量记为

线阵天线各阵元初始位置矢量，记做

在SAR观测方案设计中已经确定；根据SAR成像系统方案和观测方案，SAR成像方法需要的初始化成像系统参数均为已知；

步骤2、初始化SAR的观测场景目标空间参数：

初始化SAR的观测场景目标空间参数包括：以雷达波束照射场区域地平面和垂直于该地平面向上的单位向量所构成的空间直角坐标系作为线阵SAR的观测场景目标空间Ω；将观测场景目标空间Ω均匀划分为大小相等的立体单元格，也可称为分辨单元，单元网格在水平横向、水平纵向和高度向边长分别记为d_x,d_y和d_z，单元网格在水平横向、水平纵向和高度向单元格数分别为M_x,M_y和M_z，单元格大小选择为线阵SAR系统传统理论成像分辨率；水平横向和水平纵向构成阵列维成像空间，在阵列平面维成像空间上第m个单元格的坐标矢量，记做

m表示阵列平面维成像空间第m个单元格，m＝1,2...M，M为阵列平面维成像空间的单元格总数，M＝M_x·M_y；阵列平面维成像空间中所有单元格的散射系数按位置顺序排列成向量，记做δ，向量δ由M行1列组成；散射系数向量δ中第t个距离单元中第m个元素的散射系数，记做

初始化SAR的观测场景目标空间参数在SAR成像方案设计中已经确定；

步骤3、建立线阵SAR的线性观测矩阵：

根据步骤1中初始化得到的平台速度矢量

线阵天线各阵元初始位置矢量

和雷达系统的脉冲重复频率PRF，采用公式

n＝1,2,...,N，l＝1,2,...K，计算得到第n个线阵天线在第l个方位向慢时刻的位置矢量，记为

其中N为步骤1中线阵天线阵元总数，K为步骤1的方位向慢时刻总数；

采用公式

n＝1,2,...,N，l＝1,2,...K，m＝1,2,...M，计算得到在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格到第n个天线阵元的距离，记为

其中||·||₂表示范数，

为步骤2中初始化得到的阵列平面维成像空间中第m个单元格的坐标矢量，M为步骤2中初始化的阵列平面维成像空间中单元格总数；

采用公式

n＝1,2,...,N，l＝1,2,...K，m＝1,2,...M，t＝1,2,...T计算得到在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格到第n个线阵阵元的时间延时，记为

其中C为步骤1中初始化得到的电磁波在空气中的传播速度；

在第l个方位向慢时刻和第t个距离向快时刻中线阵SAR第n个线阵天线阵元的原始回波数据记做s(t,l,n)，t＝1,2,...T，T为步骤1初始化得到的距离向快时刻总数；在线阵SAR实际成像中，s(t,l,n)由数据接收机提供；

采用标准合成孔径雷达回波数据距离向脉冲压缩方法对原始回波数据进行距离向脉冲压缩后，得到距离向压缩后的线阵合成孔径雷达数据；记做s_AC(t,l,n)；

将所有线阵SAR第t个距离单元阵列平面维回波信号s_AC(t,l,n)按顺序排列组成向量，记为原始回波信号向量S，S由W行1列组成，其中W＝KN，K是步骤1初始化的慢时刻总数，N为步骤1初始化的线阵天线的阵元总数；

采用公式

得到阵列平面中第m个单元格在慢时间l到原始回波信号向量S中第i个元素信号对应的时延函数Φ_i(m)，其中，

为在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格到第n个线阵阵元的时间延时，t＝1,2,...T，l＝1,2,...K，n＝1,2,...,N，m＝1,2,...M，i＝1,2,...W；

令矩阵Ψ为原始回波信号向量S与散射系数向量δ之间的测量矩阵，测量矩阵由线阵SAR阵列平面维成像空间所有单元格对应的时延函数构成，具体表达式为：

其中，Φ₁(1)表示阵列平面中第1个单元格在原始回波信号向量S中第1个元素信号对应的时延函数，Φ₁(2)表示阵列平面中第2个单元格在原始回波信号向量S中第1个元素信号对应的时延函数，Φ₁(M)表示阵列平面中第M个单元格在原始回波信号向量S中第1个元素信号对应的时延函数，Φ₂(1)表示阵列平面中第1个单元格在原始回波信号向量S中第2个元素信号对应的时延函数，Φ₂(2)表示阵列平面中第2个单元格在原始回波信号向量S中第2个元素信号对应的时延函数，Φ₂(M)表示阵列平面中第M个单元格在原始回波信号向量S中第2个元素信号对应的时延函数，Φ_W(1)表示阵列平面中第1个单元格在原始回波信号向量S中第W个元素信号对应的时延函数，Φ_W(2)表示阵列平面中第2个单元格在原始回波信号向量S中第W个元素信号对应的时延函数，Φ_W(M)表示阵列平面中第M个单元格在原始回波信号向量S中第W个元素信号对应的时延函数；

步骤4、设定基于加权迭代最小稀疏贝叶斯重构算法的初始参数：

初始化加权迭代最小稀疏贝叶斯重构算法的迭代次数n＝0，采用公式

初始化线阵SAR平面维成像空间散射系数向量

其中t＝1,2,...T，m＝1,2...M；

采用公式

初始化系统噪声方差值β⁽⁰⁾，其中T为步骤1中初始化的快时间总数，Ψ为步骤3中的测量矩阵，S为步骤3中得到的原始回波信号向量，H为转置运算符号，W是步骤3中计算得到的原始回波信号向量S的行数；

定义：初始化指数分布中的参数为λ，初始化迭代误差门限为ε₀，初始化迭代次数门限为I_iter；

步骤5、估计散射系数向量和系统噪声方差：

当n＝1时，

其中

和β⁽⁰⁾由步骤4中初始化得到，采用公式

计算得到Λ⁽¹⁾；

根据Λ⁽¹⁾和

采用公式

计算得到第t个距离单元上第m个分辨单元第1次迭代后散射系数向量

再根据

采用公式

计算得到β⁽¹⁾；其中，λ是步骤4中初始化的指数分布中的参数，Ψ是步骤3中计算得到的测量矩阵，S是步骤3中得到的原始回波信号向量，H是矩阵转置符号，η是根据基于加权迭代最小稀疏贝叶斯重构算法的稳健性来选择，η选择为10^-6，N是步骤1初始化得到的线阵天线的阵元总数，W是步骤3中计算得到的原始回波信号向量S的行数；

当n＝k,k∈(1,N)，

采用公式

计算得到Λ^(k)；

根据Λ^(k)和

采用公式

计算得到第t个距离单元上第m个分辨单元第k次迭代后散射系数向量

再根据

采用公式

计算得到β^(k)；其中，λ是步骤4中初始化的指数分布中的参数，Ψ是步骤3中计算得到的测量矩阵，S是步骤3中得到的原始回波信号向量，H是矩阵转置符号，η是根据基于加权迭代最小稀疏贝叶斯重构算法算法的稳健性来选择，η选择为10^-6，N是步骤1初始化得到的线阵天线的阵元总数，W是步骤3中计算得到的原始回波信号向量S的行数；

当n＝N，

采用公式

计算得到Λ^(N)；

根据Λ^(N)和

采用公式

计算得到第t个距离单元上第m个分辨单元第N次迭代后散射系数向量

再根据

采用公式

计算得到β^(N)；其中，λ是步骤4中初始化的指数分布中的参数，Ψ是步骤3中计算得到的测量矩阵，S是步骤3中得到的原始回波信号向量，H是矩阵转置符号，η是根据基于加权迭代最小稀疏贝叶斯重构算法的稳健性来选择，η选择为10^-6，N是步骤1初始化得到的线阵天线的阵元总数，W是步骤3中计算得到的原始回波信号向量S的行数；

步骤6、迭代终止判断：

如果

并且n≤I_iter，则继续执行步骤5；

若不满足

和n≤I_iter任一条件，则基于加权迭代最小稀疏贝叶斯重构算法的迭代终止，输出

得到的基于加权迭代最小稀疏贝叶斯重构算法第n次迭代得到的散射系数向量值δ_t即为线阵SAR平面维成像空间最终的散射系数向量，其中，m＝1,2...M,t＝1,2,...,T，其中ε₀、I_iter是步骤4中初始化得到的；

步骤7、构建三维目标全场景目标成像空间：

如果t≤T，则执行步骤5到步骤6，得到所有距离单元所对应的阵列平面维成像空间散射系数向量A＝[δ₁,δ₂...δ_T]，其中T是步骤1初始化得到的距离向快时间总数；

采用公式B＝t×t×A将所有距离单元所对应的阵列平面维成像空间散射散射系数向量转换成三维矩阵形式，记为B，最终得到线阵SAR观测场景目标空间的三维成像结果，t为距离向快时间，T为距离向快时间总数。

至此，我们得到散射系数的最优估计解，整个重构方法结束。

本发明的创新点是：针对线阵SAR观测场景目标空间中主散射目标在空间上稀疏的特征，建立了线阵SAR原始回波信号与观测场景目标空间中散射系数的线性测量矩阵，并提供了基于加权迭代的最小稀疏贝叶斯重构(WSBRIM)算法，该算法在定义6中提到的迭代最小稀疏贝叶斯重构(Sparsity Bayesian Recovery via Iterative Minimum，SBRIM)算法的基础上，对L1范数进行加权，提高了重构精度和运算效率。

本发明优点在于算法可以针对不同的信号进行不同的先验概率建模，在迭代最小稀疏贝叶斯重构(SBRIM)算法的基础上对代价函数中L1范数进行加权，然后对图像进行重构，具有重构精度高、运算效率高的优势，本发明可适用于合成孔径雷达成像和地球遥感等领域。

附图说明

图1为本发明流程图；

图2为系统参数表；

具体实施方式

本发明主要采用计算机仿真的方法进行验证，所有步骤、结论都在MATLAB-R2014b上验证正确。具体实施步骤如下：

步骤1、初始化SAR系统参数：

初始化SAR系统参数包括：平台速度矢量记为

线阵天线各阵元初始位置矢量，记做

其中n为天线各阵元序号为自然数，n＝1,2,...,N，N＝4096为线阵天线的阵元总数，线阵天线长度，记做L＝3m；雷达发射信号载频f_c＝30GHz；雷达发射信号的调频斜率f_dr＝3×10¹⁴Hz/s；脉冲重复时间记为PRI＝2ms；雷达系统的脉冲重复频率PRF＝500Hz；雷达发射信号带宽B_r＝1.5×10⁸Hz，电磁波在空气中的传播速度C＝3×10⁸m/s；距离向快时间记做t，t＝1,2...T，T＝256为距离向快时刻总数，方位向慢时刻记做l，l＝1,2,...K，K＝256为方位向慢时刻总数；

步骤2、初始化SAR的观测场景目标空间参数：

初始化SAR的观测场景目标空间参数包括：以雷达波束照射场区域地平面和垂直于该地平面向上的单位向量所构成的空间直角坐标系作为线阵SAR的观测场景目标空间Ω；观测场景目标空间Ω大小为128×64×64像素，每一个单元网格在水平横向、水平纵向和高度向边长分别记为d_x＝1.25,d_y＝0.25和d_z＝0.25，单元网格在水平横向、水平纵向和高度向单元格数分别为M_x＝40,M_y＝80和M_z＝40；水平横向和水平纵向构成阵列维成像空间，在阵列平面维成像空间上每一个单元格的位置

其中x'＝1,2,...128，y'＝1,2,...64，m＝(x'-1)128+y'，

为阵列平面成像空间中第m个单元格，m表示阵列平面维成像空间第m个单元格，m＝1,2...M，M为阵列平面维成像空间的单元格总数，M＝M_x·M_y＝8192；在观测场景目标空间Ω中加入仿真点目标散射体，点目标散射体个数为5个，散射系数值均为1，坐标位置分别为[32,32,0],[32,24,0],[24,32,0],[40,32,0],单位均为m，观测场景Ω中没有包含散射点的单元格散射系数值设置为0；利用传统的合成孔径雷达原始回波仿真方法产生线阵SAR的原始回波数据；

步骤3、建立线阵SAR的线性观测矩阵：

根据步骤1中初始化得到的平台速度矢量

线阵天线各阵元初始位置矢量

和雷达系统的脉冲重复频率PRF＝500Hz，采用公式

n＝1,2,...,N，l＝1,2,...K，计算得到第n个线阵天线在第l个方位向慢时刻的位置矢量，记为

其中N＝4096为步骤1中线阵天线阵元总数，K＝256为步骤1的方位向慢时刻总数；

采用公式

n＝1,2,...,N，l＝1,2,...K，m＝1,2,...M，计算得到在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格到第n个天线阵元的距离，记为

其中||·||₂表示定义3中的向量L2范数，

为步骤2中初始化得到的阵列平面维成像空间中第m个单元格的坐标矢量，M＝8192为步骤2中初始化的阵列平面维成像空间中单元格总数；

采用公式

n＝1,2,...,N，l＝1,2,...K，m＝1,2,...M计算得到在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格到第n个线阵阵元的时间延时，记为

其中C＝3×10⁸m/s为步骤1中初始化得到的电磁波在空气中的传播速度；在第l个方位向慢时刻和第t个距离向快时刻中线阵SAR第n个线阵天线阵元的原始回波数据记做s(t,l,n)，t＝1,2,...T，T＝256为步骤1初始化得到的距离向快时刻总数；在线阵SAR实际成像中，s(t,l,n)由数据接收机提供；

采用定义2中的标准合成孔径雷达回波数据距离向脉冲压缩方法对雷达回波信号进行距离向脉冲压缩后，得到距离向压缩后的线阵合成孔径雷达数据。记做s_AC(t,l,n)；

将所有线阵SAR第t个距离单元阵列平面维回波信号s_AC(t,l,n)按顺序排列组成向量，记为原始回波信号向量S，S由W行1列组成，其中K＝256是步骤1初始化的慢时刻总数，N＝4096为步骤1初始化的线阵天线的阵元总数,W＝K·N＝1048576；

采用公式

计算得到阵列平面中第m个单元格在慢时间l到原始回波信号向量S中第i个元素信号对应的时延函数，其中，

为阵列平面中第m个分辨单元在慢时间l到第n个线阵阵元的时延时间，t＝1,2,...T，l＝1,2,...K，n＝1,2,...,N，m＝1,2,...M，i＝1,2,...W；

令矩阵Ψ为原始回波信号向量S与散射系数向量δ之间的测量矩阵，测量矩阵由线阵SAR阵列平面维成像空间所有单元格对应的时延函数构成，具体表达式为：

其中，Φ₁(1)表示阵列平面中第1个单元格在原始回波信号向量S中第1个元素信号对应的时延函数，Φ₁(2)表示阵列平面中第2个单元格在原始回波信号向量S中第1个元素信号对应的时延函数，Φ₁(M)表示阵列平面中第M个单元格在原始回波信号向量S中第1个元素信号对应的时延函数，Φ₂(1)表示阵列平面中第1个单元格在原始回波信号向量S中第2个元素信号对应的时延函数，Φ₂(2)表示阵列平面中第2个单元格在原始回波信号向量S中第2个元素信号对应的时延函数，Φ₂(M)表示阵列平面中第M个单元格在原始回波信号向量S中第2个元素信号对应的时延函数，Φ_W(1)表示阵列平面中第1个单元格在原始回波信号向量S中第W个元素信号对应的时延函数，Φ_W(2)表示阵列平面中第2个单元格在原始回波信号向量S中第W个元素信号对应的时延函数，Φ_W(M)表示阵列平面中第M个单元格在原始回波信号向量S中第W个元素信号对应的时延函数；

步骤4、设定基于加权迭代最小稀疏贝叶斯重构算法的初始参数

初始化基于加权迭代最小稀疏贝叶斯重构算法的迭代次数n＝0，利用公式

初始化线阵SAR平面维成像空间散射系数向量

采用公式

初始化系统噪声方差值β⁽⁰⁾，t＝1,2,...T，m＝1,2...M，其中T为步骤1中初始化的快时间总数，Ψ为步骤3中的测量矩阵，S为步骤3中得到的原始回波信号向量，H为转置运算符号；初始化指数分布中的参数为λ＝1，初始化迭代误差门限为ε₀＝1e-13，初始化迭代次数门限为I_iter＝20；

步骤5、更新散射系数向量和系统噪声方差：

当n＝1时，

其中

和β⁽⁰⁾由步骤4中初始化得到，采用公式

计算得到Λ⁽¹⁾。根据Λ⁽¹⁾和

采用公式

计算得到第t个距离单元上第m个分辨单元第1次迭代后散射系数向量

再根据

采用公式

计算得到β⁽¹⁾。其中，λ是步骤4中初始化的指数分布中的参数，Ψ是步骤3中计算得到的测量矩阵，S是步骤3中得到的原始回波信号向量，H是矩阵转置符号，η是根据基于加权迭代最小稀疏贝叶斯重构算法的稳健性来选择，η选择为10^-6，N＝4096是步骤1初始化得到的线阵天线的阵元总数，W＝1048576是步骤3中原始回波信号向量S的行数；

当n＝k,k∈(1,N)，

采用公式

计算得到Λ^(k)，根据Λ^(k)和

采用公式

计算得到第t个距离单元上第m个分辨单元第k次迭代后散射系数向量

再根据

采用公式

计算得到β^(k)。其中，λ是步骤4中初始化的指数分布中的参数，Ψ是步骤3中计算得到的测量矩阵，S是步骤3中得到的原始回波信号向量，H是矩阵转置符号，η是根据基于加权迭代最小稀疏贝叶斯重构算法的稳健性来选择，η选择为10^-6，N＝4096是步骤1初始化得到的线阵天线的阵元总数，W＝1048576是步骤3中原始回波信号向量S的行数；

当n＝N，

采用公式

计算得到Λ^(N)，根据Λ^(N)和

采用公式

计算得到第t个距离单元上第m个分辨单元第N次迭代后散射系数向量

再根据

采用公式

计算得到β^(N)。其中，λ是步骤4中初始化的指数分布中的参数，Ψ是步骤3中计算得到的测量矩阵，S是步骤3中得到的原始回波信号向量，H是矩阵转置符号，η是根据基于加权迭代最小稀疏贝叶斯重构算法的稳健性来选择，一般情况下η选择在10^-6附近，N＝4096是步骤1初始化得到的线阵天线的阵元总数，W＝1048576是步骤3中原始回波信号向量S的行数；

步骤6、迭代终止判断：

如果

并且n≤I_iter，则继续执行步骤5，若不满足

和n≤I_iter任一条件，则基于加权迭代最小稀疏贝叶斯重构算法的迭代终止，输出

此刻基于加权迭代最小稀疏贝叶斯重构算法第n次迭代得到的散射系数向量值δ_t即为线阵SAR平面维成像空间最终的散射系数向量，其中，m＝1,2,…,M,t＝1,2,...,T，ε₀＝1e-13、I_iter＝20；

步骤7、构建三维目标全场景目标成像空间：

如果t≤T，则执行步骤5到步骤6，得到所有距离单元所对应的阵列平面维成像空间散射系数向量A＝[δ₁,δ₂...δ_T]，其中T是步骤1初始化得到的距离向快时间总数，并采用公式B＝t×t×A将所有距离单元所对应的阵列平面维成像空间散射散射系数向量转换成三维矩阵形式，记为B，最终得到线阵SAR观测场景目标空间的三维成像结果，t＝1,2,...T为距离向快时间，T＝256为距离向快时间总数。

Claims (1)

1.一种基于加权迭代最小稀疏贝叶斯重构算法的SAR成像方法，其特征是它包括以下步骤：

步骤1、初始化SAR系统参数：

初始化SAR系统参数包括：平台速度矢量记为

线阵天线各阵元初始位置矢量，记做

其中n为天线各阵元序号为自然数，n＝1,2,...,N，N为线阵天线的阵元总数，线阵天线长度，记做L；雷达发射信号载频f_c；雷达发射信号的调频斜率f_dr；脉冲重复时间记为PRI；雷达系统的脉冲重复频率PRF；雷达发射信号带宽B_r；电磁波在空气中的传播速度C；距离向快时间记做t，t＝1,2...T，T为距离向快时刻总数，方位向慢时刻记做l，l＝1,2,...K，K为方位向慢时刻总数；上述参数均为SAR系统标准参数，其中雷达发射信号载频f_c，雷达发射信号的调频斜率f_dr，脉冲重复时间PRI，雷达系统的脉冲重复频率PRF，雷达发射信号带宽B_r，线阵天线的阵元总数N，线阵天线长度L在线阵SAR系统设计过程中已经确定；平台速度矢量记为

线阵天线各阵元初始位置矢量，记做

在SAR观测方案设计中已经确定；根据SAR成像系统方案和观测方案，SAR成像方法需要的初始化成像系统参数均为已知；

步骤2、初始化SAR的观测场景目标空间参数：

初始化SAR的观测场景目标空间参数包括：以雷达波束照射场区域地平面和垂直于该地平面向上的单位向量所构成的空间直角坐标系作为线阵SAR的观测场景目标空间Ω；将观测场景目标空间Ω均匀划分为大小相等的立体单元格，也可称为分辨单元，单元网格在水平横向、水平纵向和高度向边长分别记为d_x,d_y和d_z，单元网格在水平横向、水平纵向和高度向单元格数分别为M_x,M_y和M_z，单元格大小选择为线阵SAR系统传统理论成像分辨率；水平横向和水平纵向构成阵列维成像空间，在阵列平面维成像空间上第m个单元格的坐标矢量，记做

m表示阵列平面维成像空间第m个单元格，m＝1,2...M，M为阵列平面维成像空间的单元格总数，M＝M_x·M_y；阵列平面维成像空间中所有单元格的散射系数按位置顺序排列成向量，记做δ，向量δ由M行1列组成；散射系数向量δ中第t个距离单元中第m个元素的散射系数，记做

初始化SAR的观测场景目标空间参数在SAR成像方案设计中已经确定；

步骤3、建立线阵SAR的线性观测矩阵：

根据步骤1中初始化得到的平台速度矢量

线阵天线各阵元初始位置矢量

和雷达系统的脉冲重复频率PRF，采用公式

计算得到第n个线阵天线在第l个方位向慢时刻的位置矢量，记为

其中N为步骤1中线阵天线阵元总数，K为步骤1的方位向慢时刻总数；

采用公式

计算得到在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格到第n个天线阵元的距离，记为

其中||·||₂表示范数，

为步骤2中初始化得到的阵列平面维成像空间中第m个单元格的坐标矢量，M为步骤2中初始化的阵列平面维成像空间中单元格总数；

采用公式

计算得到在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格到第n个线阵阵元的时间延时，记为

其中C为步骤1中初始化得到的电磁波在空气中的传播速度；

在第l个方位向慢时刻和第t个距离向快时刻中线阵SAR第n个线阵天线阵元的原始回波数据记做s(t,l,n)，t＝1,2,...T，T为步骤1初始化得到的距离向快时刻总数；在线阵SAR实际成像中，s(t,l,n)由数据接收机提供；

采用标准合成孔径雷达回波数据距离向脉冲压缩方法对原始回波数据进行距离向脉冲压缩后，得到距离向压缩后的线阵合成孔径雷达数据；记做s_AC(t,l,n)；

将所有线阵SAR第t个距离单元阵列平面维回波信号s_AC(t,l,n)按顺序排列组成向量，记为原始回波信号向量S，S由W行1列组成，其中W＝K·N，K是步骤1初始化的慢时刻总数，N为步骤1初始化的线阵天线的阵元总数；

采用公式

得到阵列平面中第m个单元格在慢时间l到原始回波信号向量S中第i个元素信号对应的时延函数Φ_i(m)，其中，

为在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格到第n个线阵阵元的时间延时，t＝1,2,...T，l＝1,2,...K，n＝1,2,...,N，m＝1,2,...M，i＝1,2,...W；

令矩阵Ψ为原始回波信号向量S与散射系数向量δ之间的测量矩阵，测量矩阵由线阵SAR阵列平面维成像空间所有单元格对应的时延函数构成，具体表达式为：

其中，Φ₁(1)表示阵列平面中第1个单元格在原始回波信号向量S中第1个元素信号对应的时延函数，Φ₁(2)表示阵列平面中第2个单元格在原始回波信号向量S中第1个元素信号对应的时延函数，Φ₁(M)表示阵列平面中第M个单元格在原始回波信号向量S中第1个元素信号对应的时延函数，Φ₂(1)表示阵列平面中第1个单元格在原始回波信号向量S中第2个元素信号对应的时延函数，Φ₂(2)表示阵列平面中第2个单元格在原始回波信号向量S中第2个元素信号对应的时延函数，Φ₂(M)表示阵列平面中第M个单元格在原始回波信号向量S中第2个元素信号对应的时延函数，Φ_W(1)表示阵列平面中第1个单元格在原始回波信号向量S中第W个元素信号对应的时延函数，Φ_W(2)表示阵列平面中第2个单元格在原始回波信号向量S中第W个元素信号对应的时延函数，Φ_W(M)表示阵列平面中第M个单元格在原始回波信号向量S中第W个元素信号对应的时延函数；

步骤4、设定基于加权迭代最小稀疏贝叶斯重构算法的初始参数：

初始化加权迭代最小稀疏贝叶斯重构算法的迭代次数n＝0，采用公式

初始化线阵SAR平面维成像空间散射系数向量

其中t＝1,2,...T，m＝1,2…M；

采用公式

初始化系统噪声方差值β⁽⁰⁾，其中T为步骤1中初始化的快时间总数，Ψ为步骤3中的测量矩阵，S为步骤3中得到的原始回波信号向量，H为转置运算符号，W是步骤3中计算得到的原始回波信号向量S的行数；

定义：初始化指数分布中的参数为λ，初始化迭代误差门限为ε₀，初始化迭代次数门限为I_iter；

步骤5、估计散射系数向量和系统噪声方差：

当n＝1时，

其中

和β⁽⁰⁾由步骤4中初始化得到，采用公式

计算得到Λ⁽¹⁾；

根据Λ⁽¹⁾和

采用公式

计算得到第t个距离单元上第m个分辨单元第1次迭代后散射系数向量

再根据

采用公式

计算得到β⁽¹⁾；其中，λ是步骤4中初始化的指数分布中的参数，Ψ是步骤3中计算得到的测量矩阵，S是步骤3中得到的原始回波信号向量，H是矩阵转置符号，η是根据基于加权迭代最小稀疏贝叶斯重构算法的稳健性来选择，η选择为10^-6，N是步骤1初始化得到的线阵天线的阵元总数，W是步骤3中计算得到的原始回波信号向量S的行数；

当n＝k,k∈(1,N)，

采用公式

计算得到Λ^(k)；

根据Λ^(k)和

采用公式

计算得到第t个距离单元上第m个分辨单元第k次迭代后散射系数向量

再根据

采用公式

计算得到β^(k)；其中，λ是步骤4中初始化的指数分布中的参数，Ψ是步骤3中计算得到的测量矩阵，S是步骤3中得到的原始回波信号向量，H是矩阵转置符号，η是根据基于加权迭代最小稀疏贝叶斯重构算法算法的稳健性来选择，η选择为10^-6，N是步骤1初始化得到的线阵天线的阵元总数，W是步骤3中计算得到的原始回波信号向量S的行数；

当n＝N，

采用公式

计算得到Λ^(N)；

根据Λ^(N)和

采用公式

计算得到第t个距离单元上第m个分辨单元第N次迭代后散射系数向量

再根据

采用公式

计算得到β^(N)；其中，λ是步骤4中初始化的指数分布中的参数，Ψ是步骤3中计算得到的测量矩阵，S是步骤3中得到的原始回波信号向量，H是矩阵转置符号，η是根据基于加权迭代最小稀疏贝叶斯重构算法的稳健性来选择，η选择为10^-6，N是步骤1初始化得到的线阵天线的阵元总数，W是步骤3中计算得到的原始回波信号向量S的行数；

步骤6、迭代终止判断：

如果

并且n≤I_iter，则继续执行步骤5；

若不满足

和n≤I_iter任一条件，则基于加权迭代最小稀疏贝叶斯重构算法的迭代终止，输出

得到的基于加权迭代最小稀疏贝叶斯重构算法第n次迭代得到的散射系数向量值δ_t即为线阵SAR平面维成像空间最终的散射系数向量，其中，m＝1,2…M,t＝1,2,...,T，其中ε₀、I_iter是步骤4中初始化得到的；

步骤7、构建三维目标全场景目标成像空间：

如果t≤T，则执行步骤5到步骤6，得到所有距离单元所对应的阵列平面维成像空间散射系数向量A＝[δ₁,δ₂...δ_T]，其中T是步骤1初始化得到的距离向快时间总数；

采用公式B＝t×t×A将所有距离单元所对应的阵列平面维成像空间散射散射系数向量转换成三维矩阵形式，记为B，最终得到线阵SAR观测场景目标空间的三维成像结果，t为距离向快时间，T为距离向快时间总数。

Images