CN103698763A - 基于硬阈值omp的线阵sar稀疏成像方法 - Google Patents

Abstract

本发明提供了一种基于硬阈值OMP的线阵SAR稀疏成像方法，它是针对线阵SAR观测场景目标空间中主散射目标在空间上稀疏的特征，通过建立线阵SAR原始回波信号与观测场景目标空间中散射系数的线性测量矩阵，利用目标散射系数最大最小对比度和目标散射系数变化率作为硬阈值OMP算法迭代处理的迭代终止条件，克服了传统正交匹配追踪算法在线阵SAR稀疏成像中对主散射目标个数的依赖，与基于传统正交匹配追踪算法的线阵SAR稀疏成像相比，它无需已知观测场景目标空间的主散射目标个数，更适用于实际情况中主散射目标个数未知时的线阵SAR稀疏成像；提高了线阵SAR的成像精度。本发明可以应用于合成孔径雷达成像和地球遥感等领域。

Description

基于硬阈值OMP的线阵SAR稀疏成像方法

技术领域

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

背景技术

传统合成孔径雷达（SAR）利用单天线的直线运动合成一维的虚拟线阵天线，获得方位向高分辨率，再利用脉冲压缩技术获得雷达视线方向高分辨率，从而实现观测场景的二维成像。但由于侧视成像时地形遮挡阴影效应以及对称模糊问题，二维SAR在城市、山地和峡谷等复杂起伏地形情况不能获得满意的成像结果。三维SAR基本原理是通过天线运动合成虚拟二维面阵天线，获得面阵平面内二维高分辨，再结合脉冲压缩技术获得雷达视线方向高分辨率，实现对观测场景目标三维成像，克服了二维SAR成像技术在复杂起伏地形区域的缺陷。三维SAR是SAR成像技术未来发展的必然趋势以及当前的研究热点。线阵合成孔径雷达（Linear array SAR，LASAR，简称线阵SAR）是利用与方位向和雷达视线方向垂直放置线阵天线，再结载荷平台的运动合成二维虚拟面阵的三维SAR成像技术。与传统合成孔径雷达系统相比，线阵SAR具有多模式工作能力，除了可工作于传统侧视模式，还可工作于下视模式以及前视模式，在成像应用中更加灵活。根据每个脉冲周期内线阵阵元个数不同，线阵SAR可分为全阵元(满阵)线阵SAR和稀疏阵元线阵SAR。全阵元线阵SAR需要接收、存储和处理全采样线阵天线的回波信号，系统硬件成本高、数据量及处理难度大。为了降低系统硬件与数据处理成本，实际中线阵SAR系统中通常采用稀疏阵元，即利用稀疏线阵阵列在每个脉冲重复周期内接收一个或少数个天线阵元回波数据，但采用稀疏线阵的代价是利用经典成像方法时分辨率降低、成像质量下降。

目前线阵SAR成像方法主要基于匹配滤波(Matched Filter,MF)理论，如三维距离-多普勒(RD)算法和三维后向投影(BP)算法，见参考文献“G.Fornaro,F.Serafino,and F.Soldovieri.Three-dimensionalFocusingwithMultipass SARData.IEEE Trans.Geosci.Remote Sens,vol.41,no.3,pp.507–517,Mar.2003.”和“ShiJun,Zhang Xiaoling,YangJianyu,Wang yinbo.Surface-Tracing-Based LASAR3-D Imaging Method via Multiresolution Approximation.IEEETrans.Geosci.Remote Sens,vol.46,no.11,pp.3719–3730,Nov.2008.”，该类算法在频域或时域上通过回波数据的相参积累获得观测场景目标的三维成像。虽然传统匹配滤波算法运算效率较高，但匹配滤波成像算法由于受线阵长度和分辨率瑞利准则限制，成像分辨率较低而且在稀疏线阵时存在较严重主瓣展宽和旁瓣干扰，成像质量不足。随着各应用领域对雷达成像分辨率要求越来越高，研究获取高分辨或者超分辨能力成像算法成为了当前线阵SAR成像技术的一个迫切需要。

线阵SAR成像是从原始回波信号中重构出目标散射系数的过程，该成像过程本质上是一个逆散射问题求解过程。如果能够建立线阵SAR回波信号和观测场景目标空间散射系数的线阵测量模型，线阵SAR成像问题就可等效为三维目标空间散射系数的线性方程逆求解问题。在线阵SAR成像的三维观测场景目标空间中，由于大多数区域不包含散射点（如，空气）或散射点被其他散射点遮挡而无法被入射波束照射（如，地下目标），线阵SAR三维图像往往表现出典型的空间稀疏特征。因此线阵SAR成像问题可以进一步转化为观测场景目标空间中稀疏目标的散射系数的估计与重构问题，在成像处理过程中只需要估计出包含稀疏目标的空间单元散射系数，并不需要估计观测场景目标空间所有单元的散射系数。基于线阵SAR观测场景目标的稀疏性，稀疏信号重构理论为克服经典线阵SAR成像算法的缺陷，提高线阵SAR成像质量提供了一种新的技术途径。正交匹配追踪(OMP)算法是稀疏信号重构理论中的经典算法，详见参考文献“J.A.Tropp,A.C.Gilbert.SignalRecovery fromRandomMeasurements viaOrthogonalMatchingPursuit.IEEE Transactions on Information Theory,vol.53,no.12,pp:4655-4666,2007.”。OMP算法基本思想是利用贪婪迭代追踪方法，在每一次迭代中选定测量矩阵中与重构残余误差的最大相干项作为索引原子，并对所选定原子矩阵进行正交化以保证迭代结果最优。因正交匹配追踪算法结构简单、计算复杂度低和运算时间快等优点，近两年已在图像处理、医学成像和无线通信等领域中相关稀疏信号重构得到了成功应用。但是，正交匹配追踪算法成功重构稀疏信号的前提是需要已知原始稀疏信号中非零元素的个数，如果原始稀疏信号中非零个数未知，该算法重构性能会严重下降。然而，在线阵SAR实际成像中，观测场景中主散射目标的个数通常不是已知的。虽然观测场景中主散射目标的个数可以利用其它算法进行近似估计，但是估计处理过程复杂且结果往往不够精确。为了在线阵SAR中成功应用正交匹配追踪算法对观测场景主散射目标进行稀疏成像，必须克服正交匹配追踪算法对已知原始稀疏信号中非零元素个数的约束。

为了提高稀疏线阵SAR成像精度，并且为了克服观测场景中主散射目标个数对正交匹配追踪算法的约束，本发明提出一种基于硬阈值OMP的线阵SAR稀疏成像方法。硬阈值OMP算法利用线阵SAR观测场景中目标最大最小散射系数对比度和目标散射系数变化率代替主散射目标的个数作为OMP算法迭代终止的判定条件，通过合理的设置目标最大最小散射系数对比度和目标散射系数变化率的值，硬阈值OMP在主散射目标个数未知时也能较精确地分离观测场景中的主散射目标和弱散射背景，因而可以有效应用于主散射目标个数未知时线阵SAR稀疏成像。根据本人了解，当前还没有出现基于硬阈值OMP的线阵SAR稀疏重构成像方法。

发明内容

为了提高线阵SAR的成像精度，本发明结合线阵SAR三维观测场景目标的稀疏特征，将稀疏重构理论应用于线阵SAR成像，提出了一种基于硬阈值OMP的线阵SAR稀疏成像方法，它是针对线阵SAR观测场景目标空间中主散射目标在空间上稀疏的特征，通过建立线阵SAR原始回波信号与观测场景目标空间中散射系数的线性测量矩阵，提出了基于硬阈值OMP的线阵SAR稀疏成像算法，该算法利用目标散射系数最大最小对比度和目标散射系数变化率作为硬阈值OMP算法迭代处理的迭代终止条件，克服了传统正交匹配追踪算法在线阵SAR稀疏成像中对主散射目标个数的依赖，更加适用于实际情况中未知主散射目标个数时的线阵SAR稀疏成像，并且相对传统匹配滤波方法提高了线阵SAR的成像精度。

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

定义1、稀疏信号

如果一个离散信号中非零值的个数远小于信号本身的长度，则该信号可认为是稀疏的。设X＝[x₁,x₂,…,x_N]^T为N个离散信号组成的列向量，其中x₁表示向量X中的第1个元素，x₂表示向量X中的第2个元素，x_N表示向量X中的第N个元素，右上角T为转置运算符号。如果向量X中仅有K₀个元素非零或远大于零，则向量X定义为K₀稀疏向量，

的值定义为信号向量X的稀疏度。详见文献“S.Mallat.AWaveletTour ofSignalProcessing:TheSparseWay.Access Online via Elsevier,2008.”。

定义2、范数

设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范数表达式为

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

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

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

定义4、正交匹配追踪(OMP)算法

正交匹配追踪算法(简称OMP算法)是已知信号线性测量模型中的测量信号以及相对应的测量矩阵，在每一次迭代处理过程中从测量矩阵中选择最匹配原子构建稀疏基矩阵，再对所选定基矩阵进行正交化，然后计算正交化后稀疏基矩阵下测量信号的残差，再从测量矩阵寻找与残差最匹配的原子，经过多次迭代后利用选定的稀疏基矩阵估计稀疏信号的信号重构算法。对于定义3中的信号线性测量模型Y＝AX，其中向量Y为测量信号，矩阵A为测量矩阵，向量X为被测量信号，正交匹配追踪算法本质上是求解以下L0范数最优化问题得到信号的稀疏估计解，稀疏估计解记为

$\hat{X}=\arg \underset{X}{\min }\left({\left|\right|X\left|\right|}_0\right)s.\mathrm{tY}=\mathrm{AX}$

其中，

表示求取满足括号里面函数最小值的自变量向量X最优值，||X||₀表示向量X的L0范数，s.t表示存在或使得数学符号。在正交匹配追踪算法中，需要原始信号的非零信号个数作为算法迭代处理的终止条件才能正确重构信号。正交匹配追踪算法详见参考文献“J.A.Tropp,A.C.Gilbert.Signal recovery from random measurements via orthogonal matchingpursuit.IEEE Transactions on Information Theory,vol.53,no.12,pp:4655-4666,2007.”。

定义5、线阵合成孔径雷达（Linear array SAR,LASAR,线阵SAR）

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

定义6、线阵SAR慢时刻与快时刻

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

定义7、线阵SAR观测场景目标空间

线阵SAR观测场景目标空间是指现实空间中所有待观测场景目标散射点的位置集合。观测场景目标空间在不同空间坐标系下有不同的表示，但一旦坐标系确立以后其表示是唯一的。一般情况下为了方便成像，线阵SAR观测场景目标空间取为地面坐标系。本发明中用以下数学关系表示场景目标空间，记为Ω：

其中

和

表示构成观测场景目标空间Ω的地表正交坐标基，分别表示水平横向、水平纵向和垂直地表的高度向，

为场景目标空间中一个分辨单元位置向量，x、y和z分别表示该分布单元的水平横向、水平纵向和高度向坐标，

表示实数域。

定义8、线阵SAR成像空间

线阵SAR成像空间是指将场景目标空间中的散射点投影到切航迹向－沿航迹向－距离向的三维空间坐标系，该空间由线阵SAR成像空间中的三个相互正交的坐标基确定。本发明中用以下数学关系表示成像空间，记为M：

其中和

表示构成线阵SAR成像空间M的正交坐标基，分别表示切航迹向、沿航迹向和距离向，

为成像空间中的待观测点位置向量，u、v和z分别表示该点的切航迹向、沿航迹向和距离向坐标，

表示实数域。

定义9、线阵SAR传统理论成像分辨率

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

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

其中D_a为天线在方位向的真实孔径；切航迹向的分辨率记为ρ_c，近似表达式为

其中λ为线阵SAR雷达载频波长，R₀为线阵SAR平台到观测场景中心的参考斜距，L为线阵天线长度。详见参考文献“Shi.J,Zhang.X.L,et al.,APCTrajectoryDesign for One-ActiveLinear-arrayThree-dimensionalImaging SAR,IEEE Transactions on Geoscience and Remote Sensing,Vol.48,No.3,pp:1470-1486,2010”。

定义10、目标最大最小散射系数对比度

目标最大最小散射系数对比度是指观测场景目标空间Ω中感兴趣目标区域分辨单元最大散射系数与最小散射系数的比值。本发明中目标最大最小散射系数对比度记为η_TMM，η_TMM用以下数学关系式表示：

${\mathrm{\η}}_{\mathrm{TMM}}=\frac{\max \left|{\mathrm{\α}}_T\right|}{\min \left|{\mathrm{\α}}_T\right|}$

其中Τ表示观测场景目标空间Ω中感兴趣目标区域，α_Τ表示感兴趣目标区域Τ分辨单元对应的目标散射系数集，max|·|表示求取向量中元素的最大绝对值，min|·|表示求取向量中元素的最小绝对值。

定义11、目标散射系数变化率

目标散射系数变化率定是指观测场景目标空间Ω中分辨单元的散射系数按大小排序，组成向量后相邻元素之间散射系数的变化率。本发明中用以下数学关系表示目标散射系数变化率记为β_T，β_T用以下数学关系表示：

其中，

表示观测场景目标空间Ω中分辨单元散射系数元素按大小排序后得到散射系数向量，表示向量

第i个元素值，

表示向量第i-1个元素值，|·|表示取绝对值符号。

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

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

本发明提供的一种基于硬阈值OMP的线阵SAR稀疏成像方法，它包括以下步骤（如图1和图2所示）：

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

初始化线阵SAR系统参数包括：平台速度矢量，记做

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

其中n为天线各阵元序号，为自然数，n＝1,2,...,N，N为线阵天线的阵元总数；线阵天线长度，记做L；雷达工作中心频率，记做f_c；雷达载频波长，记做λ；雷达发射基带信号的信号带宽，记做B_r；雷达发射信号脉冲宽度，记做T_P；雷达发射信号的调频斜率，记做f_dr；雷达接收波门持续宽度，记做T_o；雷达接收系统的采样频率，记做f_s；雷达发射系统的脉冲重复频率，记做PRF；雷达系统的脉冲重复时间，记为PRI；雷达接收系统接收波门相对于发射信号发散波门的延迟，记做T_D；天线在方位向的有效孔径长度，记做D_a；光在空气中的转播速度，记做C；距离向快时刻，记做t，t＝1,2,…,T，T为距离向快时刻总数；方位向慢时刻，记做l，l＝1,2,…,K，K为方位向慢时刻总数；上述参数均为线阵SAR系统标准参数，其中线阵天线的阵元总数N，线阵天线长度L，相邻天线阵元之间的间距d，雷达中心频率f_c，雷达载频波长λ，雷达发射基带信号的信号带宽B_r，雷达发射信号脉冲宽度T_P，雷达发射信号调频斜率f_dr，雷达接收波门持续宽度T_o，雷达接收系统的采样频率f_s,雷达系统的脉冲重复频率PRF，雷达系统的脉冲重复时间PRI和雷达接收系统接收波门相对于发射信号发散波门的延迟T_D，天线在方位向的有效孔径长度D_a在线阵SAR系统设计过程中已经确定；平台速度矢量

及线阵天线各阵元初始位置矢量在线阵SAR观测方案设计中已经确定。根据线阵SAR成像系统方案和观测方案，线阵SAR成像方法需要的初始化成像系统参数均为已知。

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

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

m表示观测场景目标空间Ω中第m个单元格，m＝1,2,…,M,M为观测场景目标空间Ω中的单元格总数；观测场景目标空间Ω中所有单元格的散射系数按位置顺序排列组成向量，记做α，向量α由M行1列组成。散射系数向量α中第m个元素的散射系数，记做α_m。观测场景目标空间Ω在线阵SAR成像方案设计中已经确定。

步骤3、建立线阵SAR原始回波信号与观测场景目标散射系数的线性测量矩阵：

根据步骤1中初始化的平台速度矢量线阵天线各阵元初始位置矢量

和雷达系统的脉冲重复频率PRF，采用公式 ${\stackrel{\‾}{P}}_n\left(l\right)={\stackrel{\‾}{P}}_n\left(0\right)+\stackrel{\‾}{V}\·l/\mathrm{PRF},n=\mathrm{1,2},\·\·\·,N,l=\mathrm{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中的向量L2范数，

为步骤2中初始化得到观测场景目标空间Ω中第m个单元格的坐标矢量，M为步骤2中初始化得到的场景目标空间Ω中单元格总数。采用公式 ${\mathrm{\τ}}_{\mathrm{nm}}\left(l\right)=2\·R({\stackrel{\‾}{P}}_n\left(l\right),{\stackrel{\‾}{P}}_m)/C,n=\mathrm{1,2},\·\·\·,N,l=\mathrm{1,2},\·\·\·,K,m=\mathrm{1,2},\·\·\·,M,$ 计算得到在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格到第n个线阵天线阵元的时间延时，记为τ_nm(l)，其中C为步骤1中初始化得到的光在空气中的传播速度。在第l个方位向慢时刻和第t个距离向快时刻中线阵SAR第n个线阵天线阵元的原始回波数据记为s(t,l,n)，t＝1,2,…,T，l＝1,2,…,K，n＝1,2,…,N，其中T为步骤1中初始化得到的距离向快时刻总数。在线阵SAR实际成像中，s(t,l,n)可由数据接收机提供；而在仿真过程中，s(t,l,n)为观测场景目标空间Ω中所有单元格回波的总数，采用定义12中传统的合成孔径雷达原始回波仿真方法产生提供，近似表示公式可表示为 $s(t,l,n)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\α}}_m\·\exp [-j\·2\·\mathrm{\π}\·f_c\·{\mathrm{\τ}}_{\mathrm{nm}}\left(l\right)]\exp \{j\·\mathrm{\π}\·f_{\mathrm{dr}}\·{[t-{\mathrm{\τ}}_{\mathrm{nm}}\left(l\right)]}^2\},$ 其中Σ(·)表示求和运算符号，exp(·)表示e指数运算符号，f_c为步骤1初始化得到的雷达工作中心频率，f_dr为步骤1初始化得到的发射信号调频斜率，α_m为步骤2初始化得到的场景目标空间中第m个单元格的散射系数，t为距离向的第t个快时刻，j为虚数单位(即-1的开根值)，π为圆周率。将所有线阵SAR原始回波信号s(t,l,n)按顺序排列组成向量，记为回波信号向量S，回波信号向量S由O行1列组成，其中O＝T·K·N,T为步骤1中初始化的距离向快时刻总数，K为步骤1初始化的方位向慢时刻总数，N为步骤1初始化的线阵天线的阵元总数。采用公式φ_i(m)＝exp[j·2·π·f_c·τ_nm(l)]exp{j·π·f_dr·[t-τ_nm(l)]²}，t＝1,2,…,T，l＝1,2,…,K，n＝1,2,…,N，m＝1,2,…,M，i＝1,2,…,O，计算得到观测场景目标空间Ω中第m个单元格在回波信号向量S第i个元素信号对应的时延函数，记为φ_i(m)。令矩阵A为线阵SAR原始回波信号向量S与观测场景目标空间Ω中所有单元格散射系数向量_α之间的测量矩阵，测量矩阵A由线阵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个元素信号对应的时延函数，φ_O(1)为观测场景目标空间Ω中第1个单元格在回波信号向量S第O个元素信号对应的时延函数，φ_O(2)为观测场景目标空间Ω中第2个单元格在回波信号向量S第O个元素信号对应的时延函数，φ_O(M)为观测场景目标空间Ω中第M个单元格在回波信号向量S第O个元素信号对应的时延函数，φ₁(1),φ₁(2),…,φ₁(M)分别为观测场景目标空间Ω中第1,2,…,M个单元格在回波信号向量S第1个元素信号对应的时延函数向量，φ₂(1),φ₂(2),…,φ₂(M)分别为观测场景目标空间Ω中第1,2,…,M个单元格在回波信号向量S第2个元素信号对应的时延函数向量，φ_O(1),φ_O(2),…,φ_O(M)分别为观测场景目标空间Ω中第1,2,…,M个单元格在回波信号向量S第O个元素信号对应的时延函数向量。线阵SAR的线性测量矩阵A为O行M列的二维矩阵。

步骤4、设定硬阈值正交匹配追踪(简称OMP)算法的初始参数：

初始化硬阈值OMP算法的参数包括：硬阈值OMP算法重构迭代处理的最大迭代次数，记做MaxIter；目标散射系数最大最小对比度门限，记做η₀；目标散射射系数变化率门限，记做β₀；重构残余误差门限，记做ε₀；观测场景目标散射系数向量α的初始迭代值记为α⁽⁰⁾，一般α⁽⁰⁾的值选择为α⁽⁰⁾＝0或者α⁽⁰⁾＝A^HS，其中A为步骤3中得到的线阵SAR测量矩阵，S为步骤3中得到的线阵SAR原始回波信号向量，上标H表示共轭转置运算符号；重构残余误差的初始迭代值记为r⁽⁰⁾，一般r⁽⁰⁾的值选择为r⁽⁰⁾＝S；索引集合的初始迭代值，记为Ξ⁽⁰⁾，一般Ξ⁽⁰⁾的值选择为

其中

表示空集；k表示硬阈值OMP算法中的第k迭代次数，k的初始值设置为k＝0，且k的取值范围是从0到MaxIter。

步骤5、寻找测量矩阵与重构残余误差的最大相干项：

采用公式

计算得到硬阈值OMP算法第k次迭代过程中测量矩阵与重构残余误差的最大相干项，记为I_m，其中

表示求取满足括号中最大值时对应自变量m的最优值，||·||₂为向量的L2范数，Α_m为测量矩阵Α中的第m列，A为步骤3中得到的线阵SAR测量矩阵，右上角T为转置运算符号，r^(k-1)为算法第k-1次迭代过程中得到的重构残余误差，k表示硬阈值OMP算法中的第k迭代次数。若k＝1，r^(k-1)的值为步骤4中得到的初始迭代值r⁽⁰⁾，否则r^(k-1)通过硬阈值OMP算法第k-1次迭代过程中步骤10的重构残余误差计算提供。

步骤6、更新索引集合：

采用公式Ξ^(k)＝Ξ^(k-1)∪I_m计算得到硬阈值OMP算法第k次迭代过程中的索引集合，记为Ξ^(k)，其中Ξ^(k-1)为硬阈值OMP算法第k-1次迭代过程中得到的索引集合，I_m为步骤5得到的测量矩阵与重构残余误差的最大相干项，∪表示并集运算符号，k表示硬阈值OMP算法中的第k迭代次数。若k＝1，Ξ^(k-1)的值为步骤4中得到的初始值Ξ⁽⁰⁾；k＝2时，Ξ^(k-1)的值为硬阈值OMP算法第1次迭代过程中得到的索引集合Ξ⁽¹⁾；k＝3时，Ξ^(k-1)的值为硬阈值OMP算法第2次迭代过程中得到的索引集合Ξ⁽²⁾；依次类推，在第k次迭代中且k＞1时，Ξ^(k-1)为硬阈值OMP算法第k-1次迭代过程中得到的索引集合。

步骤7、重构观测场景目标空间的散射系数：

采用公式和计算得到硬阈值OMP算法第k次迭代过程中观测场景目标空间的散射系数向量，记为α^(k)，其中

为观测场景目标空间散射系数向量α^(k)中由索引集合Ξ^(k)对应的元素向量，Ξ^(k)为步骤6中硬阈值OMP算法第k次迭代过程中得到的索引集合，

为观测场景目标空间散射系数向量α^(k)中索引集合

对应的元素组成的向量，

为索引集合Ξ^(k)在观测场景目标空间Ω的补集，为矩阵的伪逆矩阵，右上角为矩阵伪逆运算符号，

为测量矩阵A中由索引集合对应的列组成的矩阵，A为步骤3中得到的线阵SAR测量矩阵，S为步骤3中得到的线阵SAR回波信号向量，∪表示并集运算符号，k表示硬阈值OMP算法中的第k迭代次数。

步骤8、计算目标散射系数最大最小对比度和目标散射系数变化率：

采用公式计算得到硬阈值OMP算法第k次迭代过程中的目标散射系数最大最小对比度，记为

其中其中

为步骤7中得到的观测场景目标空间散射系数向量α^(k)中由索引集合Ξ^(k)对应的元素向量，α^(k)为步骤7中得到的硬阈值OMP算法第k次迭代得到的观测场景目标空间的散射系数向量，Ξ^(k)为步骤6中第k次迭代过程中得到的索引集合，k表示硬阈值OMP算法中的第k迭代次数，max|·|表示求取向量中元素的最大绝对值，min|·|表示求取向量中元素的最小绝对值。采用公式

计算得到第k次迭代过程中的目标散射系数变化率，记为

其中α^(k)为步骤7中得到的硬阈值OMP算法第k次迭代得到观测场景目标空间的散射系数向量，α^(k-1)为硬阈值OMP算法第k-1次迭代得到观测场景目标空间的散射系数向量，α^(k-2)为硬阈值OMP算法第k-2次迭代得到观测场景目标空间的散射系数向量，||·||₂为向量L2范数。

步骤9、迭代终止判定：

如果

且则执行步骤10，否则硬阈值OMP算法终止迭代，此刻硬阈值OMP算法第k次迭代得到的散射系数向量值α^(k)即为观测场景目标空间Ω最终的散射系数向量，其中

为步骤8中得到的硬阈值OMP算法第k次迭代目标散射系数最大最小对比度，

为步骤8中得到的硬阈值OMP算法第k次迭代目标散射系数变化率，η₀为步骤4中初始化得到的目标散射系数最大最小对比度门限；β₀为步骤4中初始化得到的目标散射射系数变化率门限，k为硬阈值OMP算法中的第k迭代次数。

步骤10、计算重构残余误差：

采用公式r^(k)＝S-Αα^(k)计算得到硬阈值OMP算法第k次迭代过程中的重构残余误差，记为r^(k)，其中S为步骤3中得到的线阵SAR回波信号向量，A为步骤3中得到的线阵SAR测量矩阵，α^(k)为步骤7中得到的硬阈值OMP算法第k次迭代的观测场景目标空间散射系数向量，k表示硬阈值OMP算法中的第k迭代次数。

步骤11、迭代终止判定：

如果若r^(k)≥ε₀且k≤MaxIter，则k←k+1，则执行步骤5至步骤9，否则终止迭代，此刻硬阈值OMP算法第k次迭代得到的散射系数向量值α^(k)即为线阵SAR观测场景目标空间Ω最终的散射系数向量，其中r^(k)为步骤10中得到的硬阈值OMP算法第k次迭代的重构残余误差，k表示硬阈值OMP算法中的第k迭代次数，MaxIter为步骤4中初始化得到的硬阈值OMP算法重构迭代处理的最大迭代次数。最后将观测场景目标空间散射系数向量α^(k)转换成三维矩阵形式，得到线阵SAR观测场景目标空间Ω的三维成像结果。

本发明方法的主要思路是：利用线阵SAR雷达系统参数、运动平台参数和观测场景目标的空间参数与原始回波信号的相互关系，建立线阵SAR原始回波信号与三维观测场景目标散射系数之间的线性测量模型，然后基于该信号线性测量模型，利用硬阈值OMP方法完成原子矩阵的选择以及观测场景目标散射系数重构，而且在每一次迭代过程中利用目标最大最小散射系数对比度和目标散射系数变化率作为算法迭代终止条件。该方法的特点是：1)基于线阵SAR回波信号线性测量模型进行成像；2)在成像中利用贪婪迭代匹配追踪方法寻找最优的原子矩阵，获得稀疏散射目标的成像结果；3)利用观测场景中目标最大最小散射系数对比度和目标散射系数变化率作为迭代终止条件，克服正交匹配追踪算法对主散射目标个数已知要求的依赖。

本发明的优点在于利用目标散射系数最大最小对比度和目标散射系数变化率作为硬阈值OMP算法迭代处理的迭代终止条件，适用于实际情况中观测场景目标空间主散射目标个数未知的线阵SAR稀疏成像，提高了正交匹配追踪算法在主散射目标个数未知条件下线阵SAR稀疏成像性能。本发明可以应用于合成孔径雷达成像和地球遥感等领域。

附图说明

图1为线阵SAR成像几何关系图

其中，长平行四边形框表示线阵，黑点表示线阵阵元，线阵天线长度为L，K为方位向慢时刻总数，l为方位向第l个慢时刻，PRI表示线阵SAR发射信号的脉冲重复时间，

为线阵天线中第n个阵元在方位向第l个慢时刻的位置矢量，x_n(l)、y_n(l)和z_n(l)分别表示线阵天线中第n个阵元在方位向第l个慢时刻的水平横向、水平纵向和高度向坐标；

表示观测场景目标空间中第m单元格的位置向量，

为在方位向第l个慢时刻时观测场景目标空间中第m个单元格到线阵天线第n阵元的距离，x、y和z分别表示观测场景目标空间中水平横向、水平纵向和高度向坐标，0表示观测场景目标空间中坐标原点。

图2为本发明所提供的基于硬阈值OMP的线阵SAR稀疏成像方法的处理流程示意框图

图3为本发明具体实施方式采用的线阵SAR系统仿真参数表

具体实施方式

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

步骤1、初始化仿真所需的线阵SAR系统参数：

初始化线阵SAR系统参数的值如图3所示，包括：运动平台速度矢量

线阵天线的阵元总数N＝201，线阵天线各个阵元的初始位置矢量

其中n为第n个线阵天线的阵元序号，n＝1,2,…,N，N＝201，线阵天线的长度L＝3m，线阵天线相邻阵元之间的间距d＝0.015m，雷达中心频率f_c＝10GHz,雷达发射机基带信号的信号宽度B_r＝300MHz，雷达发射信号脉冲宽度T_P＝10^-6s，雷达发射信号调频斜率f_dr＝3×10¹⁴Hz/s，雷达接收系统的采样频率f_s＝500MHz，雷达系统发射信号的脉冲重复频率PRF＝500Hz，发射信号脉冲重复时间

线阵天线在方位向的有效孔径长度D_a＝1m，光在空气中的传播速度C＝3×10⁸m/s，距离向快时刻总数T＝256，距离向快时刻序列t＝1,2,…,T，其中T＝256，方位向慢时刻总数K＝256，方位向慢时刻序列l＝1,2,…,K，其中K＝256。

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

以雷达波束照射场区域地平面和垂直于该地平面向上的单位向量所构成的空间直角坐标作为线阵SAR的观测场景目标空间Ω。初始化观测场景目标空间Ω的大小为128×128×1像素，观测场景目标空间Ω的中心坐标位置位于[0,0,0]，每一个单元网格在水平横向、水平纵向和高度向边长为d_x＝d_y＝d_z＝0.5m，计算得到场景目标空间的总单元格数M＝16384，观测场景目标空间Ω中每一个单元格的位置为

其中x′＝1,2,…,128，y′＝1,2,…,128，m＝(x′-1)·128+y′。

为观测场景目标空间Ω中第m个单元格的位置矢量，m表示观测场景目标空间Ω中第m个单元格，m＝1,2,…,M，M＝16384。在观测场景目标空间Ω中加入仿真点目标散射体，点目标散射体数个数为5个，它们散射系数值均为1，坐标位置分别为[0,0,0]、[20,20,0]、[20,-21,0]、[-20,20,0]、[-20,-21,0]，单位均为m；观测场景目标空间Ω中没有包含点目标单元格的散射系数设置为0。将观测场景目标空间Ω中所有单元格的目标散射系数按单元格位置顺序排列组成散射向量α。确定观测场景目标空间Ω所有单元散射系数后，散射系数向量α在线阵SAR三维成像观测仿真过程中就已经确定。场景目标散射系数向量α由M行1列组成，α_m为向量α中对应场景目标空间Ω中第m个单元格的散射系数值。在本仿真观测场景目标空间Ω中，只有包含点散射目标的5个单元格散射系数值α设置为1，其余单元格的散射系数都为0。利用传统的合成孔径雷达原始回波仿真方法产生线阵SAR的原始回波数据

步骤3、建立线阵SAR原始回波信号与场景目标散射系数的线性测量矩阵：

采用公式

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

其中n表示天线阵元序号n＝1,2,…,N，N＝201，l表示方位向慢时刻序号l＝1,2,…,K，K＝256，运动平台初始位置矢量

运动平台速度矢量

脉冲重复频率PRF＝500Hz。采用公式

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

其中n＝1,2,…,N，N＝201，l＝1,2,…,K，K＝256，m表示场景目标空间Ω中第m个单元格，m＝1,2,…,M，M＝16384，||·||₂表示向量L2范数，

为步骤2中初始化得到观测场景目标空间Ω中第m个单元格的坐标矢量。采用公式计算得到在第l个方位向慢时刻线阵SAR场景目标空间Ω中第m个单元格到第n个线阵天线阵元的时间延时τ_nm(l)，其中n＝1,2,…,N，l＝1,2,…,K，m＝1,2,…,M，C＝3×10⁸m/s。采用公式 $s(t,l,n)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\α}}_m\·\exp [-j\·2\·\mathrm{\π}\·f_c\·{\mathrm{\τ}}_{\mathrm{nm}}\left(l\right)]\exp \{j\·\mathrm{\π}\·f_{\mathrm{dr}}\·{[t-{\mathrm{\τ}}_{\mathrm{nm}}\left(l\right)]}^2\},$ 和合成孔径雷达原始回波仿真方法得到观测场景目标空间Ω的近似线阵SAR原始信号回波s(t,l,n)，其中t表示距离向快时刻序号，t＝1,2,…,T，T＝256，l＝1,2,…,K，n＝1,2,…,N，Σ(·)表示求和运算符号，exp(·)表示e指数运算符号，雷达中心频率f_c＝10GHz，雷达发射信号调频斜率f_dr＝3×10¹⁴Hz/s，α_m为步骤2初始化得到的场景目标空间中第m个单元格的散射系数，j为虚数单位(即-1的开根值)，π＝3.14159。将所有线阵SAR原始回波信号s(t,l,n)按顺序排列组成回波信号向量S，回波信号向量S由O行1列组成，其中O＝T·K·N＝13172736。采用公式φ_i(m)＝exp［-j·2·π·f_c·τ_nm(l)]exp{j·π·f_dr·[t-τ_nm(l)]²}，t＝1,2,…,T，T＝256，l＝1,2,…,K，K＝256，n＝1,2,…,N，N＝201，m＝1,2,…,M，M＝16384，i＝[(t-1)K+l-1]N+n，计算得到观测场景目标空间Ω中第m个单元格在回波信号向量S第i个元素信号对应的时延函数φ_i(m)，其中i的取值范围为i＝1,2,…,O，O＝13172736。

采用矩阵表达公式

计算得到线阵SAR原始回波信号与场景目标空间所有单元格的线性测量矩阵A，线性测量矩阵A为O行M列的二维矩阵。

步骤4、设定硬阈值OMP算法的初始参数：

初始化硬阈值OMP算法稀疏重构处理的最大迭代次数为

目标散射系数最大最小对比度门限值η₀＝0.25；目标散射射系数变化率门限值β₀＝0.2；重构的残余误差门限值ε₀＝0.001；线阵SAR观测场景目标空间散射系数向量α的初始值α⁽⁰⁾选择为α⁽⁰⁾＝A^HS，其中A为步骤3中得到的线阵SAR测量矩阵，S为步骤3中得到的线阵SAR原始回波信号向量，上标H表示共轭转置元算符号；重构残余误差的初始迭代值r⁽⁰⁾选择为r⁽⁰⁾＝S；索引集合的初始迭代值Ξ⁽⁰⁾选择为

其中

表示空集；k表示硬阈值OMP算法中的第k迭代次数，k的初始值设置为k＝0，且k的取值范围是从0到MaxIter，其中MaxIter＝2048。

步骤5、寻找测量矩阵与重构残余误差的最大相干项：

采用公式

计算得到硬阈值OMP算法第k次迭代过程中测量矩阵与重构残余误差的最大相干项，记为I_m，其中

表示求取满足括号中最大值时对应自变量m的最优值，||·||₂为向量L2范数，Α_m为测量矩阵Α中的第m列，A为步骤3中得到的线阵SAR测量矩阵，右上角T为转置运算符号，r^(k-1)为硬阈值OMP算法第k-1次迭代过程中得到的重构残余误差。若k＝1，r^(k-1)的值为步骤4中得到的初始迭代值r⁽⁰⁾，否则r^(k-1)通过硬阈值OMP算法第k-1次迭代过程中步骤10的重构残余误差计算提供，k表示硬阈值OMP算法中的第k迭代次数。

步骤6、更新索引集合：

采用公式Ξ^(k)＝Ξ^(k-1)∪I_m计算得到硬阈值OMP算法第k次迭代过程中的索引集合，记为Ξ^(k)，其中Ξ^(k-1)为硬阈值OMP算法第k-1次迭代过程中得到的索引集合，I_m为步骤5得到的测量矩阵与重构残余误差的最大相干项，∪表示并集符号。若k＝1，Ξ^(k-1)的值为步骤4中得到的初始值Ξ⁽⁰⁾；k＝2时，Ξ^(k-1)的值为硬阈值OMP算法第1次迭代过程中得到的索引集合Ξ⁽¹⁾；k＝3时，Ξ^(k-1)的值为硬阈值OMP算法第2次迭代过程中得到的索引集合Ξ⁽²⁾；依次类推，在第k次迭代中且k＞1时，Ξ^(k-1)为硬阈值OMP算法第k-1次迭代过程中得到的索引集合，k表示硬阈值OMP算法中的第k迭代次数。

步骤7、重构观测场景目标空间散射系数：

采用公式

和

计算得到硬阈值OMP算法第k次迭代过程中观测场景目标空间的散射系数向量，记为α^(k)，其中

为观测场景目标空间散射系数向量α^(k)中由索引集合Ξ^(k)对应的元素向量，Ξ^(k)为步骤6中硬阈值OMP算法第k次迭代过程中得到的索引集合，

为观测场景目标空间散射系数向量α^(k)中索引集合

对应的元素组成的向量，

为索引集合Ξ^(k)在观测场景目标空间Ω的补集，

为矩阵

的伪逆矩阵，右上角

为矩阵伪逆运算符号，

为测量矩阵A中由索引集合

对应的列组成的矩阵，A为步骤3中得到的线阵SAR测量矩阵，S为步骤3中得到的线阵SAR回波信号向量，∪表示并集运算符号，k表示硬阈值OMP算法中的第k迭代次数。

步骤8、计算目标散射系数最大最小对比度和目标散射系数变化率：

采用公式

计算得到硬阈值OMP算法第k次迭代过程中的目标散射系数最大最小对比度，记为

其中其中

为步骤7中得到的由索引集合Ξ^(k)对应的观测场景目标空间散射系数向量α^(k)元素，max|·|表示求取向量中元素的最大绝对值，min|·|表示求取向量中元素的最小绝对值。采用公式

计算得到硬阈值OMP算法第k次迭代过程中的目标散射系数变化率，记为

其中α^(k)为硬阈值OMP算法第k次迭代得到观测场景目标空间散射系数向量，α^(k-1)为硬阈值OMP算法第k-1次迭代得到观测场景目标空间散射系数向量，α^(k-2)为硬阈值OMP算法第k-2次迭代得到观测场景目标空间散射系数向量，k为硬阈值OMP算法中的第k迭代次数，||·||₂为向量L2范数。

步骤9、迭代终止判定：

如果若

且

则k←k+1，则执行步骤10，若不满足

和

任一条件，则硬阈值OMP算法停止迭代过程，此刻硬阈值OMP算法第k次迭代得到的散射系数向量值α^(k)即为线阵观测场景目标空间Ω最终的散射系数向量，将得到的观测场景目标空间散射系数向量α^(k)转换成128×128×1三维矩阵形式，得到线阵SAR观测场景目标空间Ω的三维成像结果，其中

为步骤8中得到第k次迭代目标散射系数最大最小对比度，为步骤8中得到第k次迭代目标散射系数变化率，η₀为步骤4中初始化得到的目标散射系数最大最小对比度门限η₀＝0.25，β₀为步骤4中初始化得到的目标散射射系数变化率门限β₀＝0.2，k为硬阈值OMP算法中的第k迭代次数。

步骤10、计算重构的残余误差：

采用公式r^(k)＝S-Αα^(k)计算得到硬阈值OMP算法第k次迭代过程中的重构残余误差，记为r^(k)，其中S为步骤3中得到的线阵SAR回波信号向量，A为步骤3中得到的线阵SAR测量矩阵，α^(k)为步骤7中得到第k次迭代的观测场景目标空间散射系数向量，k为硬阈值OMP算法中的第k迭代次数。

步骤11、迭代终止判定：

如果若r^(k)≥ε₀且k≤MaxIter，则k←k+1，则执行步骤5到步骤9，若不满足k≤MaxIter与r^(k)≥ε0任意一个条件，则硬阈值OMP算法停止迭代过程，此刻硬阈值OMP算法第k次迭代得到的散射系数向量值α^(k)即为线阵SAR观测场景目标空间Ω最终的散射系数向量，其中r^(k)为步骤10中得到第k次迭代重构残余误差，k为硬阈值OMP算法中的第k迭代次数，ε₀为步骤4中初始化得到的算法重构残余误差门限值ε₀＝0.001，MaxIter为步骤4中初始化得到的硬阈值OMP算法重构处理的最大迭代次数MaxIter=2048。最后将最终得到的观测场景目标空间散射系数向量α^(k)转换成128×128×1三维矩阵形式，得到线阵SAR观测场景目标空间Ω的三维成像结果。

通过本发明具体实施方式可以看出，本发明通过建立线阵SAR原始回波信号与场景目标空间散射系数的线性测量模型，并结合线阵SAR观测场景目标空间的稀疏特征，将线阵SAR成像处理过程转换成为正交匹配追踪算法的稀疏求解过程。本发明提供了基于硬阈值OMP的线阵SAR稀疏成像方法，该方法结合线阵SAR系统参数和目标场景参数，建立回波信号与稀疏目标散射系数的线性测量模型，利用目标散射系数最大最小对比度和目标散射系数变化率作为算法迭代处理的终止条件，与基于传统正交匹配追踪算法的线阵SAR稀疏成像相比，它无需已知观测场景目标空间的主散射目标个数，更适用于实际情况中主散射目标个数未知时的线阵SAR稀疏成像。

Claims (1)

1.一种基于硬阈值OMP的线阵SAR稀疏成像方法，其特征是它包括以下步骤：

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

初始化线阵SAR系统参数包括：平台速度矢量，记做

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

其中n为天线各阵元序号，为自然数，n＝1,2,...,N，N为线阵天线的阵元总数；线阵天线长度，记做L；雷达工作中心频率，记做f_c；雷达载频波长，记做λ；雷达发射基带信号的信号带宽，记做B_r；雷达发射信号脉冲宽度，记做T_P；雷达发射信号的调频斜率，记做f_dr；雷达接收波门持续宽度，记做T_o；雷达接收系统的采样频率，记做f_s；雷达发射系统的脉冲重复频率，记做PRF；雷达系统的脉冲重复时间，记为PRI；雷达接收系统接收波门相对于发射信号发散波门的延迟，记做T_D；天线在方位向的有效孔径长度，记做D_a；光在空气中的转播速度，记做C；距离向快时刻，记做t，t＝1,2,…,T，T为距离向快时刻总数；方位向慢时刻，记做l，l＝1,2,…,K，K为方位向慢时刻总数；上述参数均为线阵SAR系统标准参数，其中线阵天线的阵元总数N，线阵天线长度L，相邻天线阵元之间的间距d，雷达中心频率f_c，雷达载频波长λ，雷达发射基带信号的信号带宽B_r，雷达发射信号脉冲宽度T_P，雷达发射信号调频斜率f_dr，雷达接收波门持续宽度T_o，雷达接收系统的采样频率f_s,雷达系统的脉冲重复频率PRF，雷达系统的脉冲重复时间PRI和雷达接收系统接收波门相对于发射信号发散波门的延迟T_D，天线在方位向的有效孔径长度D_a在线阵SAR系统设计过程中已经确定；平台速度矢量

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

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

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

初始化线阵SAR的观测场景目标空间参数，包括：以雷达波束照射场区域地平面和垂直于该地平面向上的单位向量所构成的空间直角坐标作为线阵三维SAR的观测场景目标空间Ω；将观测场景目标空间Ω均匀划分成大小相等的立体单元格，单元网格在水平横向、水平纵向和高度向边长分别记为d_x、d_y和d_z，单元格大小一般选择为线阵SAR系统传统理论成像分辨率或该分辨率的二分之一；观测场景目标空间Ω中第m个单元格的坐标矢量，记做m表示观测场景目标空间Ω中第m个单元格，m＝1,2,…,M,M为观测场景目标空间Ω中的单元格总数；观测场景目标空间Ω中所有单元格的散射系数按位置顺序排列组成向量，记做α，向量α由M行1列组成；散射系数向量α中第m个元素的散射系数，记做α_m；观测场景目标空间Ω在线阵SAR成像方案设计中已经确定；

步骤3、建立线阵SAR原始回波信号与观测场景目标散射系数的线性测量矩阵：

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

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

和雷达系统的脉冲重复频率PRF，采用公式 ${\stackrel{\‾}{P}}_n\left(l\right)={\stackrel{\‾}{P}}_n\left(0\right)+\stackrel{\‾}{V}\·l/\mathrm{PRF},n=\mathrm{1,2},\·\·\·,N,l=\mathrm{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中的向量L2范数，

为步骤2中初始化得到观测场景目标空间Ω中第m个单元格的坐标矢量，M为步骤2中初始化的场景目标空间Ω中单元格总数；采用公式 ${\mathrm{\τ}}_{\mathrm{nm}}\left(l\right)=2\·R({\stackrel{\‾}{P}}_n\left(l\right),{\stackrel{\‾}{P}}_m)/C,n=\mathrm{1,2},\·\·\·,N,l=\mathrm{1,2},\·\·\·,K,m=\mathrm{1,2},\·\·\·,M,$ 计算得到在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格到第n个线阵天线阵元的时间延时，记为τ_nm(l)，其中C为步骤1中初始化得到的光在空气中的传播速度；在第l个方位向慢时刻和第t个距离向快时刻中线阵SAR第n个线阵天线阵元的原始回波数据记为s(t,l,n)，t＝1,2,…,T，l＝1,2,…,K，n＝1,2,…,N，其中T为步骤1中初始化的距离向快时刻总数；在线阵SAR实际成像中，s(t,l,n)可由数据接收机提供；而在仿真过程中，s(t,l,n)为观测场景目标空间Ω中所有单元格回波的总数，采用定义12中传统的合成孔径雷达原始回波仿真方法产生提供，近似表示公式可表示为 $s(t,l,n)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\α}}_m\·\exp [-j\·2\·\mathrm{\π}\·f_c\·{\mathrm{\τ}}_{\mathrm{nm}}\left(l\right)]\exp \{j\·\mathrm{\π}\·f_{\mathrm{dr}}\·{[t-{\mathrm{\τ}}_{\mathrm{nm}}\left(l\right)]}^2\},$ 其中Σ(·)表示求和运算符号，exp(·)表示e指数运算符号，f_c为步骤1初始化得到的雷达工作中心频率，f_dr为步骤1初始化得到的发射信号调频斜率，α_m为步骤2初始化得到的场景目标空间中第m个单元格的散射系数，t为距离向的第t个快时刻，j为虚数单位，即-1的开根值，π为圆周率；将所有线阵SAR原始回波信号s(t,l,n)按顺序排列组成向量，记为回波信号向量S，回波信号向量S由O行1列组成，其中O＝T·K·N,T为步骤1中初始化的距离向快时刻总数，K为步骤1初始化的方位向慢时刻总数，N为步骤1初始化的线阵天线的阵元总数；采用公式φ_i(m)＝exp[-j·2·π·f_c·τ_nm(l)]exp{j·π·f_dr·[t-τ_nm(l)]²}，t＝1,2,…,T，l＝1,2,…,K，n＝1,2,…,N，m＝1,2,…,M，i＝1,2,…,O，计算得到观测场景目标空间Ω中第m个单元格在回波信号向量S第i个元素信号对应的时延函数，记为φ_i(m)；令矩阵A为线阵SAR原始回波信号向量S与观测场景目标空间Ω中所有单元格散射系数向量α之间的测量矩阵，测量矩阵A由线阵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个元素信号对应的时延函数，φ_O(1)为观测场景目标空间Ω中第1个单元格在回波信号向量S第O个元素信号对应的时延函数，φ_O(2)为观测场景目标空间Ω中第2个单元格在回波信号向量S第O个元素信号对应的时延函数，φ_O(M)为观测场景目标空间Ω中第M个单元格在回波信号向量S第O个元素信号对应的时延函数，φ₁(1),φ₁(2),…,φ₁(M)分别为观测场景目标空间Ω中第1,2,…,M个单元格在回波信号向量S第1个元素信号对应的时延函数向量，φ₂(1),φ₂(2),…,φ₂(M)分别为观测场景目标空间Ω中第1,2,…,M个单元格在回波信号向量S第2个元素信号对应的时延函数向量，φ_O(1),φ_O(2),…,φ_O(M)分别为观测场景目标空间Ω中第1,2,…,M个单元格在回波信号向量S第O个元素信号对应的时延函数向量；线阵SAR的线性测量矩阵A为O行M列的二维矩阵；

步骤4、设定硬阈值正交匹配追踪(简称OMP)算法的初始参数：

初始化硬阈值OMP算法的参数包括：硬阈值OMP算法重构迭代处理的最大迭代次数，记做MaxIter；目标散射系数最大最小对比度门限，记做η₀；目标散射射系数变化率门限，记做β₀；重构残余误差门限，记做ε₀；观测场景目标散射系数向量α的初始迭代值记为α⁽⁰⁾，一般α⁽⁰⁾的值选择为α⁽⁰⁾＝0或者α⁽⁰⁾＝A^HS，其中A为步骤3中得到的线阵SAR测量矩阵，S为步骤3中得到的线阵SAR原始回波信号向量，上标H表示共轭转置运算符号；重构残余误差的初始迭代值记为r⁽⁰⁾，一般r⁽⁰⁾的值选择为r⁽⁰⁾＝S；索引集合的初始迭代值，记为Ξ⁽⁰⁾，一般Ξ⁽⁰⁾的值选择为

其中表示空集；k表示硬阈值OMP算法中的第k迭代次数，k的初始值设置为k＝0，且k的取值范围是从0到MaxIter；

步骤5、寻找测量矩阵与重构残余误差的最大相干项：

采用公式计算得到硬阈值OMP算法第k次迭代过程中测量矩阵与重构残余误差的最大相干项，记为I_m，其中

表示求取满足括号中最大值时对应自变量m的最优值，||·||₂为向量的L2范数，Α_m为测量矩阵Α中的第m列，A为步骤3中得到的线阵SAR测量矩阵，右上角T为转置运算符号，r^(k-1)为算法第k-1次迭代过程中得到的重构残余误差，k表示硬阈值OMP算法中的第k迭代次数；若k＝1，r^(k-1)的值为步骤4中得到的初始迭代值r⁽⁰⁾，否则r^(k-1)通过硬阈值OMP算法第k-1次迭代过程中步骤10的重构残余误差计算提供；

步骤6、更新索引集合：

采用公式Ξ^(k)＝Ξ^(k-1)∪I_m计算得到硬阈值OMP算法第k次迭代过程中的索引集合，记为Ξ^(k)，其中Ξ^(k-1)为硬阈值OMP算法第k-1次迭代过程中得到的索引集合，I_m为步骤5得到的测量矩阵与重构残余误差的最大相干项，∪表示并集运算符号，k表示硬阈值OMP算法中的第k迭代次数；若k＝1，Ξ^(k-1)的值为步骤4中得到的初始值Ξ⁽⁰⁾；k＝2时，Ξ^(k-1)的值为硬阈值OMP算法第₁次迭代过程中得到的索引集合Ξ⁽¹⁾；k＝3时，Ξ^(k-1)的值为硬阈值OMP算法第2次迭代过程中得到的索引集合Ξ⁽²⁾；依次类推，在第k次迭代中且k＞1时，Ξ^(k-1)为硬阈值OMP算法第k-1次迭代过程中得到的索引集合；

步骤7、重构观测场景目标空间的散射系数：

采用公式

和

计算得到硬阈值OMP算法第k次迭代过程中观测场景目标空间的散射系数向量，记为α^(k)，其中为观测场景目标空间散射系数向量α^(k)中由索引集合Ξ^(k)对应的元素向量，Ξ^(k)为步骤6中硬阈值OMP算法第k次迭代过程中得到的索引集合，

为观测场景目标空间散射系数向量α^(k)中索引集合

对应的元素组成的向量，

为索引集合Ξ^(k)在观测场景目标空间Ω的补集，

为矩阵

的伪逆矩阵，右上角

为矩阵伪逆运算符号，

为测量矩阵A中由索引集合

对应的列组成的矩阵，A为步骤3中得到的线阵SAR测量矩阵，S为步骤3中得到的线阵SAR回波信号向量，∪表示并集运算符号，k表示硬阈值OMP算法中的第k迭代次数；

步骤8、计算目标散射系数最大最小对比度和目标散射系数变化率：

采用公式

计算得到硬阈值OMP算法第k次迭代过程中的目标散射系数最大最小对比度，记为其中其中

为步骤7中得到的观测场景目标空间散射系数向量α^(k)中由索引集合Ξ^(k)对应的元素向量，α^(k)为硬阈值OMP算法第k次迭代得到的观测场景目标空间的散射系数向量，Ξ^(k)为步骤6中第k次迭代过程中得到的索引集合，k表示硬阈值OMP算法中的第k迭代次数，max|·|表示求取向量中元素的最大绝对值，min|·|表示求取向量中元素的最小绝对值；采用公式

计算得到第k次迭代过程中的目标散射系数变化率，记为其中α^(k)为硬阈值OMP算法第k次迭代得到观测场景目标空间的散射系数向量，α^(k-1)为硬阈值OMP算法第k-1次迭代得到观测场景目标空间的散射系数向量，α^(k-2)为硬阈值OMP算法第k-2次迭代得到观测场景目标空间的散射系数向量，||·||₂为向量L2范数；

步骤9、迭代终止判定：

如果

且

则执行步骤10，否则硬阈值OMP算法终止迭代，此刻硬阈值OMP算法第k次迭代得到的散射系数向量值α^(k)即为观测场景目标空间Ω最终的散射系数向量，其中

为步骤8中得到的硬阈值OMP算法第k次迭代目标散射系数最大最小对比度，

为步骤8中得到的硬阈值OMP算法第k次迭代目标散射系数变化率，η₀为步骤4中初始化得到的目标散射系数最大最小对比度门限；β₀为步骤4中初始化得到的目标散射射系数变化率门限，k为硬阈值OMP算法中的第k迭代次数；

步骤10、计算重构残余误差：

采用公式r^(k)＝S-Αα^(k)计算得到硬阈值OMP算法第k次迭代过程中的重构残余误差，记为r^(k)，其中S为步骤3中得到的线阵SAR回波信号向量，A为步骤3中得到的线阵SAR测量矩阵，α^(k)为步骤7中得到的硬阈值OMP算法第k次迭代的观测场景目标空间散射系数向量，k表示硬阈值OMP算法中的第k迭代次数；

步骤11、迭代终止判定：

如果若r^(k)≥ε₀且k≤MaxIter，则k←k+1，则执行步骤5至步骤9，否则终止迭代，此刻硬阈值OMP算法第k次迭代得到的散射系数向量值α^(k)即为线阵SAR观测场景目标空间Ω最终的散射系数向量，其中r^(k)为步骤10中得到的硬阈值OMP算法第k次迭代的重构残余误差，k表示硬阈值OMP算法中的第k迭代次数，MaxIter为步骤4中初始化得到的硬阈值OMP算法重构迭代处理的最大迭代次数；最后将观测场景目标空间散射系数向量α^(k)转换成三维矩阵形式，得到线阵SAR观测场景目标空间Ω的三维成像结果。

Images