CN107037429B - 基于门限梯度追踪算法的线阵sar三维成像方法 - Google Patents

Abstract

本发明提供了一种门限梯度追踪算法的线阵SAR三维成像方法，它是利用线阵SAR雷达系统参数、运动平台参数和观测场景目标的空间参数与原始回波信号的相互关系，建立线阵SAR原始回波信号与三维观测场景目标散射系数之间的线性测量模型，然后基于该信号线性测量模型，利用TBGP方法完成观测场景目标散射系数重构，并且在每一次迭代过程中利用目标最大最小散射系数对比度和目标散射系数变化率作为算法迭代终止条件，提高了GP算法在观测场景稀疏度未知条件下线阵SAR稀疏成像性能，相对OMP算法提高了运算效率和空间存储效率，可以应用于合成孔径雷达成像和地球遥感等领域。

Description

基于门限梯度追踪算法的线阵SAR三维成像方法

技术领域

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

背景技术

传统合成孔径雷达(SAR)利用单天线的直线运动合成一维的虚拟线阵天线，获得方位向高分辨率，再利用脉冲压缩技术获得雷达视线方向高分辨率，从而实现观测场景的二维成像。但传统二维SAR主要工作于侧视成像模式，而侧视成像时存在地形遮挡、阴影效应及顶底倒置问题，故传统二维SAR在城市、山地和峡谷等复杂起伏地形情况不能获得满意的成像结果。三维SAR基本原理是通过天线运动合成虚拟二维面阵天线，获得面阵平面内二维高分辨，再结合脉冲压缩技术获得雷达视线方向高分辨率，实现对观测目标三维成像，克服了传统二维SAR成像技术在复杂起伏地形区域的缺陷。三维SAR是SAR成像技术未来发展的必然趋势以及当前的研究热点。目前，三维SAR成像技术主要包括曲线SAR、圆周SAR、层析SAR和线阵SAR。线阵SAR(Linear array SAR，LASAR)是利用与方位向垂直放置线阵天线，再结合载荷平台的运动合成二维虚拟面阵的三维SAR成像技术。与其它三维SAR成像技术相比，线阵SAR具有多模式工作能力，除了可工作于传统侧视模式，还可工作于下视模式以及前视模式，在地形观测成像应用中更加灵活。

然而，目前线阵三维SAR成像技术主要基于传统奈奎斯特(Nyquist)采样定理和经典信号处理理论。在传统成像体制和数据处理方法前提下，若要实现高分辨对地观测，线阵三维SAR成像技术面临以下诸多挑战：(1)载荷条件受限：对于ka波段机载平台，在3km高度要实现0.05m分辨率观测，则要求线阵长度达到300m以上，现有机载平台尺寸难以满足；(2)线阵阵元数过多：为了抑制旁瓣和栅瓣，通常要采用全采样满阵线阵天线，少则需要采用上千个阵元，多则上万个阵元，硬件系统实现复杂、成本高。(3)成像精度不足：目前线阵三维SAR数据处理主要采用后向投影(Back Projection，BP)算法，BP算法受瑞利准则约束，图像存在主瓣旁瓣模糊，弱散射目标容易被邻近强散射目标旁瓣掩盖，且受线阵长度限制切航迹分辨率往往远低于其它两维，制约线阵三维SAR应用；(4)回波数据量过大：线阵三维SAR所有阵元在每一个脉冲重复时间内都采集回波，回波数据量往往是传统SAR上千倍，数据采样、存储、传输和处理难度增大、成本也大大增加。

因此，针对传统线阵三维SAR成像体制和数据处理方法的缺陷，开展具有高分辨能力的线阵三维SAR成像新机制和新处理方法研究成为亟待解决的课题。近年来，一种新兴的压缩传感(Compressed Sensing，CS)理论为稀疏信号精确重构技术带来了革命性的突破，并在压缩成像、信道编码、图像处理、医学成像、模式识别、无线网络和雷达技术等不同应用领域引起高度关注。CS理论主要采用非自适应线性投影来保持原始信号的结构信息，再通过求解方程最优解重构出原始信号。CS理论指出只要原始信号存在稀疏性或可压缩性，就可用远低于Nyquist采样率的采样信号恢复出原始信号，并且信号稀疏性越强，稀疏重构所需的观测数据越少。一般情况下，线阵三维SAR目标散射体在三维成像空间中具有很强的空间稀疏性。因此，CS新理论的出现为克服传统线阵三维SAR成像技术和数据处理方法的缺陷提供了契机。线阵三维SAR成像技术可以与CS新理论良好地结合，用低于Nyquist采样率进行回波欠采样，如切航迹向采用稀疏线阵、沿航迹向稀疏发射脉冲及距离向稀疏采样等，从而通过稀疏重构方法实现稀疏目标高精度三维成像。

近年来，基于压缩感知理论的一些算法在三维SAR成像中具有广泛的应用。OMP算法作为一种常见的贪婪算法在SAR三维成像中也有重要应用，但是对于大规模问题，OMP算法在运行时间和空间存储上效率比较低。Thomas Blumensath等人提出了GP算法，该算法在运行时间和空间存储上效率比OMP算法高，但是同OMP算法一样，GP算法也需要预设场景的稀疏度，而在SAR三维成像中一般无法获取场景的真实稀疏度。因而限制了它们的应用。

发明内容：

为了提高线阵SAR的成像精度，结合线阵SAR三维观测场景目标的稀疏特征，将稀疏重构理论应用于线阵SAR成像，本发明提供了一种门限梯度追踪算法的线阵SAR三维成像方法，利用线阵SAR雷达系统参数、运动平台参数和观测场景目标的空间参数与原始回波信号的相互关系，建立线阵SAR原始回波信号与三维观测场景目标散射系数之间的线性测量模型，然后基于该信号线性测量模型，利用TBGP方法完成观测场景目标散射系数重构，而且在每一次迭代过程中利用目标最大最小散射系数对比度和目标散射系数变化率作为算法迭代终止条件。该方法的特点是：1)基于线阵SAR回波信号线性测量模型进行成像；2)在成像中利用门限梯度追踪方法获得稀疏散射目标的成像结果；3)利用观测场景中目标最大最小散射系数对比度和目标散射系数变化率作为迭代终止条件，克服梯度追踪算法需要预设场景的稀疏度。

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

定义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.A wavelet tour of signal processing:the sparseway.Access Online via Elsevier,2008.”。

定义2、合成孔径雷达距离向压缩方法

合成孔径雷达距离向压缩方法是指利用合成孔径雷达发射参数进行匹配滤波的过程，主要包括：采用距离向参考信号对合成孔径雷达距离向回波信号进行匹配滤波压缩和距离走动补偿得到合成孔径雷达的成像。详细内容可参考文献“雷达成像技术”，保铮等编著，电子工业出版社出版。

定义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、信号线性测量模型

对于一个数字信号测量系统，假设定义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的测量矩阵。

定义5、梯度追踪算法

梯度追踪算法(简称GP算法)是已知信号线性测量模型中的测量信号以及相对应的测量矩阵，在每一次迭代处理过程中从测量矩阵中选择最大相干的原子构建稀疏基矩阵，更新方向、步长和估计信号，然后计算残差，再从测量矩阵中寻找与残差最大相干的原子，经过多次迭代后利用选定的稀疏基矩阵估计稀疏信号的信号重构算法。

定义6、线阵合成孔径雷达(Linear array SAR,LASAR,线阵SAR)

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

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

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

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

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

其中

和

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

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

表示实数域。

定义9、线阵SAR成像空间

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

其中

和

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

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

表示实数域。

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

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

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

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

其中λ为线阵SAR雷达载频波长，R₀为线阵SAR平台到观测场景中心的参考斜距，L为线阵天线长度。详见参考文献“Shi,J.,etal.,APC trajectory design for One-Active linear-array three-dimensionalimaging SAR,IEEE Transactions on Geoscience and Remote Sensing,Vol.48,No.3,pp:1470-1486,2010”。

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

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

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

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

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

其中，

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

表示向量

第i个元素值，

表示向量

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

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

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

本发明提供的基于门限梯度追踪算法的线阵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，单元网格在水平横向、水平纵向和高度向单元格数分别为M_x、M_y和M_z，单元格大小选择为线阵SAR系统传统理论成像分辨率或该分辨率的二分之一；水平横向和水平纵向构成阵列平面维成像空间，在阵列平面维成像空间上第m个单元格的坐标矢量，记做

m表示阵列平面维成像空间第m个单元格，m＝1,2,…,M,M为阵列平面维成像空间的单元格总数,M＝M_x·M_y；阵列平面维成像空间中所有单元格的散射系数按位置顺序排列组成向量，记做α，向量α由M行1列组成；散射系数向量α中第m个元素的散射系数，记做α_m；观测场景目标空间Ω在线阵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个线阵天线阵元的距离，记为

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

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

采用公式

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

其中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_AC(t,l,n)，近似表示公式可表示为

t=1,2,…,T，l=1,2,…,K，n=1,2,…,N，m=1,2,…,M,其中

表示第t个距离单元中阵列平面维成像空间上第m个分辨单元对应的散射系数值，

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

采用公式

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、设定基于门限梯度追踪算法(Threshold-Based Gradient Pursuit,简称TBGP)的初始参数：

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

其中

表示空集；k表示TBGP算法中的第k迭代次数，k的初始值设置为k＝0；

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

采用公式

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

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

步骤6、更新索引集合：

采用公式Ξ^(k)＝Ξ^(k-1)∪I_m计算得到TBGP算法第k次迭代过程中的索引集合，记为Ξ^(k)，其中Ξ^(k-1)为TBGP算法第k-1次迭代过程中得到的索引集合，I_m为步骤5得到的测量矩阵与重构残余误差的最大相干项，∪表示并集符号，k表示TBGP算法中的第k迭代次数；

若k＝1，Ξ^(k-1)的值为步骤4中得到的初始值Ξ⁽⁰⁾；k＝2时，Ξ^(k-1)的值为TBGP算法第1次迭代过程中得到的索引集合Ξ⁽¹⁾；k＝3时，Ξ^(k-1)的值为TBGP算法第2次迭代过程中得到的索引集合Ξ⁽²⁾；依次类推，在第k次迭代中且k＞1时，Ξ^(k-1)为TBGP算法第k-1次迭代过程中得到的索引集合；

步骤7、计算更新方向：

采用公式

计算得到TBGP算法第k次迭代过程中的更新方向，记为

其中

为测量矩阵Α中由索引集合Ξ^(k)对应的列组成的矩阵，A为步骤3中得到的线阵SAR测量矩阵，右上角T为转置运算符号，Ξ^(k)为步骤6中TBGP算法第k次迭代过程中得到的索引集合，S为步骤3中得到的线阵SAR阵列平面维信号向量，

为观测场景目标散射系数向量α^(k-1)中由索引集合Ξ^(k)对应的元素向量，α^(k-1)为步骤9中TBGP算法第k-1次迭代过程中的阵列平面维成像空间的散射系数估计，k表示TBGP算法中的第k迭代次数，若k＝1，α^(k-1)的值为步骤4中得到的初始迭代值α⁽⁰⁾，否则α^(k-1)通过TBGP算法第k-1次迭代过程中步骤9的阵列平面维成像空间的散射系数估计计算；

步骤8、计算步长：

采用公式

计算得到TBGP算法第k次迭代过程中的步长，记为a^(k)，其中r^(k)为TBGP算法第k次迭代过程中得到的重构残余误差，

为测量矩阵Α中由索引集合Ξ^(k)对应的列组成的矩阵，Ξ^(k)为步骤6中TBGP算法第k次迭代过程中得到的索引集合，

为步骤7中得到的TBGP算法第k次迭代过程中的更新方向，<,>为向量内积，

为向量的L2范数的平方，k表示TBGP算法中的第k迭代次数；

步骤9、更新阵列平面维成像空间的散射系数估计：

采用公式

和

计算得到TBGP算法第k次迭代过程中的阵列平面维成像空间的散射系数估计，记为α^(k)，其中

为阵列平面维成像空间的散射系数向量α^(k)中由索引集合Ξ^(k)对应的元素组成的向量，Ξ^(k)为步骤6中TBGP算法第k次迭代过程中得到的索引集合，

为阵列平面维成像空间散射系数向量α^(k)中索引集合

对应的元素组成的向量，

为索引集合Ξ^(k)在阵列平面维成像空间Ω的补集，

为阵列平面维成像空间的散射系数向量α^(k-1)中由索引集合Ξ^(k)对应的元素组成的向量，α^(k-1)为TBGP算法第k-1次迭代过程中的阵列平面维成像空间的散射系数估计，a^(k)为步骤8中TBGP算法第k次迭代过程中的步长，

为步骤7中TBGP算法第k次迭代过程中的更新方向，k表示TBGP算法中的第k迭代次数；

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

采用公式

计算得到TBGP算法第k次迭代过程中的目标散射系数最大最小对比度，记为η^(k)，其中

为步骤9中得到的阵列平面维成像空间散射系数向量α^(k)中由索引集合Ξ^(k)对应的元素向量，α^(k)为TBGP算法第k次迭代得到的阵列平面维成像空间的散射系数向量，Ξ^(k)为步骤6中第k次迭代过程中得到的索引集合，k表示TBGP算法中的第k迭代次数，max|·|表示求取向量中元素的最大绝对值，min|·|表示求取向量中元素的最小绝对值；采用公式

计算得到第k次迭代过程中的目标散射系数变化率，记为β^(k)，其中α^(k)为TBGP算法第k次迭代得到阵列平面维成像空间的散射系数向量，α^(k-1)为TBGP算法第k-1次迭代得到阵列平面维成像空间的散射系数向量，α^(k-2)为TBGP算法第k-2次迭代得到阵列平面维成像空间的散射系数向量，||·||₂为向量L2范数；

步骤11、基于目标散射系数最大最小对比度门限和目标散射系数变化率门限的迭代终止判定：

如果η^(k)≤η₀且β^(k)≤β₀，则执行步骤12，否则TBGP算法终止迭代，此时TBGP算法第k次迭代得到的散射系数向量值α^(k)即为阵列平面维成像空间最终的散射系数向量，其中η^(k)为步骤10中得到的TBGP算法第k次迭代目标散射系数最大最小对比度，β^(k)为步骤10中得到的TBGP算法第k次迭代目标散射系数变化率，η₀为步骤4中初始化得到的目标散射系数最大最小对比度门限；β₀为步骤4中初始化得到的目标散射系数变化率门限，k为TBGP算法中的第k迭代次数；

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

采用公式r^(k)＝r^(k-1)-a^(k)c^(k)计算得到TBGP算法第k次迭代过程中的重构残余误差，记为r^(k)，其中S为步骤4中得到的线阵SAR回波信号向量，A为步骤4中得到的线阵SAR测量矩阵，α^(k)为步骤11中得到的TBGP算法第k次迭代的阵列平面维成像空间散射系数向量，k表示TBGP算法中的第k迭代次数；

步骤13、迭代终止判定：

如果若r^(k)≥ε₀且k≤MaxIter，则k←k+1，则执行步骤5至11，否则终止迭代，此刻TBGP算法第k次迭代得到第t个距离单元所对应的散射系数向量值α^(k)即为线阵SAR阵列平面维成像空间最终的散射系数向量，其中r^(k)为步骤12中得到的TBGP算法第k次迭代的重构残余误差，ε₀为步骤4中初始化得到的算法重构残余误差门限值，k表示TBGP算法中的第k迭代次数，MaxIter为步骤4中初始化得到的TBGP算法重构迭代处理的最大迭代次数；

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

如果t≤T，则执行步骤3至步骤13，得到所有距离单元所对应的阵列平面维成像空间散射系数向量α^(k)，并将所有距离单元所对应的阵列平面维成像空间散射系数向量α^(k)转换成三维矩阵形式，最终得到线阵SAR观测场景目标空间Ω的三维成像结果，t为距离向快时刻，t＝1,2,…,T，T为距离向快时刻总数。

本发明的实质是在GP算法的基础上提出了一种基于门限的梯度追踪(TBGP)算法。该算法既保留了GP算法相对于OMP算法在运行时间和存储空间上的优势，又克服了GP算法需要预设稀疏度的缺点。TBGP算法用散射系数对比度、散射系数变化率代替GP算法的预设稀疏度来控制迭代的次数。根据本人了解，当前还没有出现基于门限梯度追踪算法的线阵SAR三维成像方法。

本发明的创新点在于针对线阵SAR观测场景目标空间中主散射目标在空间上稀疏的特征，建立了线阵SAR原始回波信号与观测场景目标空间中散射系数的线性测量矩阵，并提供了基于TBGP的线阵SAR稀疏成像算法，该算法利用目标散射系数最大最小对比度和目标散射系数变化率作为TBGP算法迭代处理的迭代终止条件，克服了传统GP算法在线阵SAR稀疏成像中对观测场景稀疏度的依赖，更加适用于实际情况中未知观测场景稀疏度时的线阵SAR稀疏成像，并且相对OMP算法提高了运算效率和空间存储效率。

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

附图说明

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

图1中长平行四边形框表示线阵，黑点表示线阵阵元，线阵天线长度为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为本发明所提供的基于TBGP的线阵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×128像素，观测场景目标空间Ω的中心坐标位置位于[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；在本仿真观测场景目标空间Ω中，只有包含点散射目标的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范数，

为中初始化得到阵列平面维成像空间中第m个单元格的坐标矢量；采用公式

计算得到在第l个方位向慢时刻线阵SAR阵列平面维成像空间中第m个单元格到第n个线阵天线阵元的时间延时

其中n＝1,2,…,N，l＝1,2,…,K，m＝1,2,…,M，C＝3×10⁸m/s；

采用合成孔径雷达原始回波仿真方法得到观测场景目标空间Ω的近似线阵SAR原始信号回波s(t,l,n)；采用合成孔径雷达距离向压缩方法对每个接收天线阵元的距离向回波数据进行距离向二维压缩，得到距离向压缩后的线阵合成孔径雷达数据，记做s_AC(t,l,n)，近似表示公式可表示为

t＝1,2,…,T，l＝1,2,…,K，n＝1,2,…,N，m＝1,2,…,M,其中

表示第t个距离单元中阵列平面维成像空间上第m个分辨单元对应的散射系数值，

为阵列平面中第m个分辨单元在慢时间l到第n个线阵阵元的时延时间；将所有线阵SAR第t个距离单元的阵列平面维回波信号s_AC(t,l,n)按顺序排列组成向量，记为回波信号向量S，回波信号向量S由O行1列组成，其中O＝K·N＝51456；采用公式

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＝51456；

采用矩阵表达公式

计算得到线阵SAR原始回波信号与阵列平面维成像空间所有单元格的线性测量矩阵A，线性测量矩阵A为O行M列的二维矩阵；

步骤4、设定TBGP算法的初始参数：

初始化TBGP算法稀疏重构处理的最大迭代次数为

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

其中

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

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

采用公式

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

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

步骤6、更新索引集合：

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

步骤7、计算更新方向：

采用公式

计算得到TBGP算法第k次迭代过程中的更新方向，记为

其中

为测量矩阵Α中由索引集合Ξ^(k)对应的列组成的矩阵，A为步骤3中得到的线阵SAR测量矩阵，右上角T为转置运算符号，Ξ^(k)为步骤6中TBGP算法第k次迭代过程中得到的索引集合，S为步骤3中得到的线阵SAR原始回波信号向量，

为观测场景目标散射系数向量α^(k-1)中由索引集合Ξ^(k)对应的元素向量，α^(k-1)为步骤9中TBGP算法第k-1次迭代过程中的阵列平面维成像空间的散射系数估计，k表示TBGP算法中的第k迭代次数，若k＝1，α^(k-1)的值为步骤4中得到的初始迭代值α⁽⁰⁾，否则α^(k-1)通过TBGP算法第k-1次迭代过程中步骤9的阵列平面维成像空间的散射系数估计计算；

步骤8、计算步长：

采用公式

计算得到TBGP算法第k次迭代过程中的步长，记为a^(k)，其中r^(k)为TBGP算法第k次迭代过程中得到的重构残余误差，

为测量矩阵Α中由索引集合Ξ^(k)对应的列组成的矩阵，Ξ^(k)为步骤6中TBGP算法第k次迭代过程中得到的索引集合，

为步骤7中得到的TBGP算法第k次迭代过程中的更新方向，<,>为向量内积，

为向量的L2范数的平方，k表示TBGP算法中的第k迭代次数；

步骤9、更新阵列平面维成像空间的散射系数估计：

采用公式

和

计算得到TBGP算法第k次迭代过程中的阵列平面维成像空间的散射系数估计，记为α^(k)，其中

为阵列平面维成像空间的散射系数向量α^(k)中由索引集合Ξ^(k)对应的元素组成的向量，Ξ^(k)为步骤6中TBGP算法第k次迭代过程中得到的索引集合，

为阵列平面维成像空间散射系数向量α^(k)中索引集合

对应的元素组成的向量，

为索引集合Ξ^(k)在阵列平面维成像空间的补集，

为阵列平面维成像空间的散射系数向量α^(k-1)中由索引集合Ξ^(k)对应的元素组成的向量，α^(k-1)为TBGP算法第k-1次迭代过程中的阵列平面维成像空间的散射系数估计，a^(k)为步骤8中TBGP算法第k次迭代过程中的步长，

为步骤7中TBGP算法第k次迭代过程中的更新方向，k表示TBGP算法中的第k迭代次数；

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

采用公式

计算得到TBGP算法第k次迭代过程中的目标散射系数最大最小对比度，记为η^(k)，其中其中

为步骤9中得到的阵列平面维成像空间散射系数向量α^(k)中由索引集合Ξ^(k)对应的元素向量，α^(k)为TBGP算法第k次迭代得到的阵列平面维成像空间的散射系数向量，Ξ^(k)为步骤6中第k次迭代过程中得到的索引集合，k表示TBGP算法中的第k迭代次数，max|·|表示求取向量中元素的最大绝对值，min|·|表示求取向量中元素的最小绝对值；采用公式

计算得到第k次迭代过程中的目标散射系数变化率，记为β^(k)，其中α^(k)为TBGP算法第k次迭代得到阵列平面维成像空间的散射系数向量，α^(k-1)为TBGP算法第k-1次迭代得到阵列平面维成像空间的散射系数向量，α^(k-2)为TBGP算法第k-2次迭代得到阵列平面维成像空间的散射系数向量，||·||₂为向量L2范数；

步骤11、基于目标散射系数最大最小对比度门限和目标散射系数变化率门限的迭代终止判定：

如果η^(k)≤η₀且β^(k)≤β₀，则k←k+1，执行步骤10，若不满足

和

任一条件，则TBGP算法终止迭代，此刻TBGP算法第k次迭代得到的散射系数向量值α^(k)即为阵列平面维成像空间最终的散射系数向量，其中η^(k)为步骤10中得到的TBGP算法第k次迭代目标散射系数最大最小对比度，β^(k)为步骤10中得到的TBGP算法第k次迭代目标散射系数变化率，η₀为步骤4中初始化得到的目标散射系数最大最小对比度门限η₀=0.25；β₀为步骤4中初始化得到的目标散射系数变化率门限β₀=0.2，k为TBGP算法中的第k迭代次数；

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

采用公式r^(k)=r^(k-1)-a^(k)c^(k)计算得到TBGP算法第k次迭代过程中的重构残余误差，记为r^(k)，其中S为步骤3中得到的线阵SAR回波信号向量，A为步骤3中得到的线阵SAR测量矩阵，α^(k)为步骤10中得到的TBGP算法第k次迭代的阵列平面维成像空间散射系数向量，k表示TBGP算法中的第k迭代次数；

步骤13、迭代终止判定：

如果若r^(k)≥ε₀且k≤MaxIter，则k←k+1，则执行步骤5至11，若不满足k≤MaxIter与r^(k)≥ε₀任意一个条件，则TBGP算法停止迭代过程，此刻TBGP算法第k次迭代得到的散射系数向量值α^(k)即为线阵SAR阵列平面维成像空间最终的散射系数向量，其中r^(k)为步骤10中得到的TBGP算法第k次迭代的重构残余误差，ε₀为步骤4中初始化得到的算法重构残余误差门限值ε₀＝0.001，k表示TBGP算法中的第k迭代次数，MaxIter为步骤四中初始化得到的TBGP算法重构迭代处理的最大迭代次数MaxIter＝2048；

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

执如果t≤T，则执行步骤3至步骤10，得到所有距离单元所对应的阵列平面维成像空间散射系数向量α^(k)，并将所有距离单元所对应的阵列平面维成像空间散射系数向量α^(k)转换成128×128×128三维矩阵形式，最终得到线阵SAR观测场景目标空间Ω的三维成像结果。

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

Claims (1)

1.一种基于门限梯度追踪算法的线阵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，单元网格在水平横向、水平纵向和高度向单元格数分别为M_x、M_y和M_z，单元格大小选择为线阵SAR系统传统理论成像分辨率或该分辨率的二分之一；水平横向和水平纵向构成阵列平面维成像空间，在阵列平面维成像空间上第m个单元格的坐标矢量，记做

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

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

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

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

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

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

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

采用公式

计算得到在第l个方位向慢时刻线阵SAR阵列平面维成像空间中第m个单元格到第n个线阵天线阵元的距离，记为

||·||₂表示向量的L2范数，

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

采用公式

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

其中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_AC(t,l,n)，近似表示公式可表示为

其中

表示第t个距离单元中阵列平面维成像空间上第m个分辨单元对应的散射系数值，

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

采用公式

计算得到阵列平面维成像空间中第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、设定基于门限梯度追踪TBGP算法的初始参数：

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

其中

表示空集；k表示TBGP算法中的第k迭代次数，k的初始值设置为k=0；

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

采用公式

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

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

步骤6、更新索引集合：

采用公式Ξ^(k)＝Ξ^(k-1)∪I_m计算得到TBGP算法第k次迭代过程中的索引集合，记为Ξ^(k)，其中Ξ^(k-1)为TBGP算法第k-1次迭代过程中得到的索引集合，I_m为步骤5得到的测量矩阵与重构残余误差的最大相干项，∪表示并集符号，k表示TBGP算法中的第k迭代次数；

若k＝1，Ξ^(k-1)的值为步骤4中得到的初始值Ξ⁽⁰⁾；k=2时，Ξ^(k-1)的值为TBGP算法第1次迭代过程中得到的索引集合Ξ⁽¹⁾；k=3时，Ξ^(k-1)的值为TBGP算法第2次迭代过程中得到的索引集合Ξ⁽²⁾；依次类推，在第k次迭代中且k＞1时，Ξ^(k-1)为TBGP算法第k-1次迭代过程中得到的索引集合；

步骤7、计算更新方向：

采用公式

计算得到TBGP算法第k次迭代过程中的更新方向，记为

其中

为测量矩阵Α中由索引集合Ξ^(k)对应的列组成的矩阵，A为步骤3中得到的线阵SAR测量矩阵，右上角T为转置运算符号，Ξ^(k)为步骤6中TBGP算法第k次迭代过程中得到的索引集合，S为步骤3中得到的线阵SAR阵列平面维信号向量，

为观测场景目标散射系数向量α^(k-1)中由索引集合Ξ^(k)对应的元素向量，α^(k-1)为步骤9中TBGP算法第k-1次迭代过程中的阵列平面维成像空间的散射系数估计，k表示TBGP算法中的第k迭代次数，若k＝1，α^(k-1)的值为步骤4中得到的初始迭代值α⁽⁰⁾，否则α^(k-1)通过TBGP算法第k-1次迭代过程中步骤9的阵列平面维成像空间的散射系数估计计算；

步骤8、计算步长：

采用公式

计算得到TBGP算法第k次迭代过程中的步长，记为a^(k)，其中r^(k)为TBGP算法第k次迭代过程中得到的重构残余误差，

为测量矩阵Α中由索引集合Ξ^(k)对应的列组成的矩阵，Ξ^(k)为步骤6中TBGP算法第k次迭代过程中得到的索引集合，

为步骤7中得到的TBGP算法第k次迭代过程中的更新方向，<,>为向量内积，

为向量的L2范数的平方，k表示TBGP算法中的第k迭代次数；

步骤9、更新阵列平面维成像空间的散射系数估计：

采用公式

和

计算得到TBGP算法第k次迭代过程中的阵列平面维成像空间的散射系数估计，记为α^(k)，其中

为阵列平面维成像空间的散射系数向量α^(k)中由索引集合Ξ^(k)对应的元素组成的向量，Ξ^(k)为步骤6中TBGP算法第k次迭代过程中得到的索引集合，

为阵列平面维成像空间散射系数向量α^(k)中索引集合Ξ^(k)对应的元素组成的向量，Ξ^(k)为索引集合Ξ^(k)在阵列平面维成像空间Ω的补集，

为阵列平面维成像空间的散射系数向量α^(k-1)中由索引集合Ξ^(k)对应的元素组成的向量，α^(k-1)为TBGP算法第k-1次迭代过程中的阵列平面维成像空间的散射系数估计，a^(k)为步骤8中TBGP算法第k次迭代过程中的步长，

为步骤7中TBGP算法第k次迭代过程中的更新方向，k表示TBGP算法中的第k迭代次数；

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

采用公式

计算得到TBGP算法第k次迭代过程中的目标散射系数最大最小对比度，记为η^(k)，其中

为步骤9中得到的阵列平面维成像空间散射系数向量α^(k)中由索引集合Ξ^(k)对应的元素向量，α^(k)为TBGP算法第k次迭代得到的阵列平面维成像空间的散射系数向量，Ξ^(k)为步骤6中第k次迭代过程中得到的索引集合，k表示TBGP算法中的第k迭代次数，max|·|表示求取向量中元素的最大绝对值，min|·|表示求取向量中元素的最小绝对值；采用公式

计算得到第k次迭代过程中的目标散射系数变化率，记为β^(k)，其中α^(k)为TBGP算法第k次迭代得到阵列平面维成像空间的散射系数向量，α^(k-1)为TBGP算法第k-1次迭代得到阵列平面维成像空间的散射系数向量，α^(k-2)为TBGP算法第k-2次迭代得到阵列平面维成像空间的散射系数向量，||·||₂为向量L2范数；

步骤11、基于目标散射系数最大最小对比度门限和目标散射系数变化率门限的迭代终止判定：

如果η^(k)≤η₀且β^(k)≤β₀，则执行步骤12，否则TBGP算法终止迭代，此时TBGP算法第k次迭代得到的散射系数向量值α^(k)即为阵列平面维成像空间最终的散射系数向量，其中η^(k)为步骤10中得到的TBGP算法第k次迭代目标散射系数最大最小对比度，β^(k)为步骤10中得到的TBGP算法第k次迭代目标散射系数变化率，η₀为步骤4中初始化得到的目标散射系数最大最小对比度门限；β₀为步骤4中初始化得到的目标散射系数变化率门限，k为TBGP算法中的第k迭代次数；

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

采用公式r^(k)=r^(k-1)-a^(k)c^(k)计算得到TBGP算法第k次迭代过程中的重构残余误差，记为r^(k)，其中S为步骤4中得到的线阵SAR回波信号向量，A为步骤4中得到的线阵SAR测量矩阵，α^(k)为步骤11中得到的TBGP算法第k次迭代的阵列平面维成像空间散射系数向量，k表示TBGP算法中的第k迭代次数；

步骤13、迭代终止判定：

如果若r^(k)≥ε₀且k≤MaxIter，则k←k+1，则执行步骤5至11，否则终止迭代，此刻TBGP算法第k次迭代得到第t个距离单元所对应的散射系数向量值α^(k)即为线阵SAR阵列平面维成像空间最终的散射系数向量，其中r^(k)为步骤12中得到的TBGP算法第k次迭代的重构残余误差，ε₀为步骤4中初始化得到的算法重构残余误差门限值，k表示TBGP算法中的第k迭代次数，MaxIter为步骤4中初始化得到的TBGP算法重构迭代处理的最大迭代次数；

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

如果t≤T，则执行步骤3至步骤13，得到所有距离单元所对应的阵列平面维成像空间散射系数向量α^(k)，并将所有距离单元所对应的阵列平面维成像空间散射系数向量α^(k)转换成三维矩阵形式，最终得到线阵SAR观测场景目标空间Ω的三维成像结果，t为距离向快时刻，t＝1,2,…,T，T为距离向快时刻总数。

Images