CN104391295A - 一种图像熵最优的压缩传感sar稀疏自聚焦成像方法 - Google Patents

Abstract

本发明公开了一种图像熵最优的压缩传感SAR稀疏自聚焦成像方法，它是针对SAR回波信号测量模型中方位向相位误差对压缩传感SAR成像的影响，以及压缩传感SAR成像模型中的未知相位误差估计与补偿问题，利用成像模型中方位向相位误差特性和稀疏目标特征，采用压缩传感SAR图像熵作为评价准则，在每一次迭代处理过程中利用SAR方位向回波与观测目标的关系估计方位向相位误差，然后对压缩传感成像模型进行相位误差补偿，接着再进行压缩传感SAR成像，并采用逐次迭代方法使得压缩传感SAR成像的图像熵最优，从而提高了压缩传感SAR成像质量。

Description

一种图像熵最优的压缩传感SAR稀疏自聚焦成像方法

技术领域：

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

背景技术：

由于具有全天时、全天候和大场景观测等优势，合成孔径雷达(SAR)已成为当前大面积地形测绘的一项重要遥感技术，在地形测绘、自然灾害监测和自然资源调查等领域发挥越来越大的作用。压缩传感稀疏重构作为一种近几年新提出的信号处理理论，突破了传统Nyquist采样定理约束，可利用远低于Nyquist采样率精确重构原始稀疏信号(详见参考文献“D.L.Donoho.Compressed sensing.IEEE Transactions on Information Theory,2006,52(4):1289-1306”)，在降低SAR***采样率和提高成像质量等方面有着巨大的应用潜力。因此，压缩传感SAR成像已经成为了SAR领域的新兴热点课题。压缩传感重构方法对信号测量模型的精确度要求很高，如果信号测量模型存在误差或者不精确，压缩传感重构方法的重构精确度就会严重退化，有时甚至会出现错误的结果。对于压缩传感SAR实际数据成像，影响其成像效果的一个重要因素是运动误差估计与补偿(详见参考文献“Jihua Tian,Sun Jinping,Han Xiao,and Zhang Bingchen.Motion compensation for Compressive Sensing SAR imaging with autofocus.IEEE Conference on Industrial Electronics and Applications(ICIEA),1564-1567,2011”)。在压缩传感SAR成像时需要利用每个方位向平台位置建立SAR信号观测模型。但是，由于平台运动等因素存在测量误差，使得每个方位向回波信号中混入相位误差，导致压缩传感SAR成像质量下降甚至无法成像。因此在压缩传感SAR实际成像过程中，必须先对SAR数据每个方位向的未知相位误差进行有效估计与补偿，才可对SAR原始回波数据进行高精度压缩传感成像。 

熵作为描述信息不确定性的尺度，在信号处理领域有着广泛的应用。对于SAR成像，散焦图像对比度较低，表明图像各像素取值概率接近，不确定性较大，而聚焦图像对比度较高，说明各像素取值的确定性增加。因此，图像熵可作为衡量SAR成像图像质量优劣的重要尺度。对基于传统成像方法的SAR数据处理中，图像熵优化已被用于实现自聚焦成像的一个重要手段(详见参考文献“X.Li,G.S.Liu,and J.L.Ni.Autofocusing of ISAR images based on entropy minimization.IEEE Transactions on Aerospace and Electronic Systems,1999,35(4):1240-1252”)。然而，基于压缩传感的SAR成像方法与传统SAR成像方法在信号处理上有本质的差别，基于压缩传感的SAR成像方法可从欠采样回波数据中精确重构观测目标，从而实现高精度成像， 但是传统SAR成像方法很难实现欠采样回波数据高精度成像。目前，基于图像熵最优的传统自聚焦方法是针对全采样条件下的SAR回波数据进行处理，在欠采样情况下这类方法自聚焦性能会严重下降甚至失效。因此，为了在压缩传感SAR成像中利用图像熵进行自聚焦成像，不能直接利用传统图像熵自聚焦成像方法，必须要结合压缩传感SAR信号重构模型和重构结果进行处理。 

发明内容：

为了解决压缩传感SAR稀疏成像过程中未知方位向相位误差的估计与补偿问题，本发明根据压缩传感SAR成像模型中方位向相位误差特性和稀疏目标特征，结合压缩传感稀疏重构方法以及SAR图像熵，提出了一种图像熵最优的压缩传感SAR稀疏自聚焦成像方法，本发明利用成像模型中方位向相位误差特性和稀疏目标特征，利用压缩传感SAR图像熵作为评价准则，提出了一种图像熵最优的压缩传感SAR稀疏自聚焦成像方法，该方法在每一次迭代处理过程中利用SAR方位向回波与观测目标的关系估计方位向相位误差，然后对压缩传感成像模型进行相位误差补偿，接着再进行压缩传感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₀稀疏向量。详见文献“S.Mallat.A Wavelet Tour of Signal Processing:The Sparse Way.Access Online via Elsevier,2008”。 

定义2、范数 

设X是数域上线性空间，表示复数域，若它满足如下性质：||X||≥0，且||X||＝0仅有X＝0，||aX||＝|a|||X||，a为任意常数，||X₁+X₂||≤||X₁||+||X₂||，则称||X||为X空间上的范数，||·||表示范数符号，其中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。详见文献“矩阵理论”，黄廷祝等编著，高等教育出版社出版。 

定义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、合成孔径雷达慢时刻和快时刻 

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

定义5、合成孔径雷达标准距离压缩 

合成孔径雷达标准距离压缩是指利用合成孔径雷达发射信号参数，采用匹配滤波技术对合成孔径雷达的距离向信号进行滤波的过程。详见文献“雷达成像技术”，保铮等编著，电子工业出版社出版。 

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

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

定义7、相位误差向量和相位误差矩阵 

相位误差向量是指由SAR测量回波数据中每一个方位向相位误差值所组成的向量，相位误差矩阵是指由SAR测量数据中每一个方位向相位误差值所组成的对角矩阵。 

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

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

定义9、SAR成像空间 

SAR成像空间是指将场景目标空间中的散射点投影到方位向－距离向的空间坐标系，该空间由SAR成像空间中两个相互正交的坐标基确定。 

定义10、压缩传感稀疏重构方法 

压缩传感主要是将高维原始信号进行非自适应线性投影到低维空间以保持信号的结构信息，再通过求解线性最优解重构出原始信号的理论，该理论主要包括信号稀疏表示、稀疏测量和稀疏重构三个方面。压缩传感稀疏重构方法的基本思想为求解特定约束条件下的最优解或次最优解，主要方法有贪婪追踪算法和凸松弛算法等。详细内容可参考文献“Donoho D L.Compressed sensing.IEEE Transactions on Information Theory,2006,52(4):1289-1306.”。 

定义11、对角矩阵与单位矩阵 

对角矩阵首先是行数和列数相等的方阵，同时矩阵的主对角元素不全为零，矩阵非主对角元素全为零。单位矩阵指的是主对角元素全为1，非主对角元素全为0的对角矩阵。 

定义12、图像熵 

图像熵表示为图像灰度级集合的比特平均数，也描述了图像信源的平均信息量，图像模糊度越大，图像熵越大。详细内容可以参考文献“Brink A D.Using spatial information as an aid to maximum entropy image threshold selection.Pattern Recognition Letters,1996,17(1):29-36.” 

定义13、施密特正交化方法 

施密特正交化方法的基本思想是利用投影原理在已有正交基的基础上构造一个新的正交基，把一组线性无关的向量变成一单位正交向量组的方法。详细内容可以参考“线性代数”，居余马等编著，清华大学出版社出版。 

本发明提供的一种图像熵最优的SAR稀疏自聚焦成像方法，它包括以下步骤： 

步骤1、初始化SAR***参数： 

初始化SAR***参数包括：平台速度矢量，记做V；天线初始位置矢量，即方位向0慢时刻天线位置，记做P(0)；雷达工作中心频率，记做f_c；雷达载频波长，记做λ；雷达发射基带信号的信号带宽，记做B_r；雷达发射信号脉冲宽度，记做T_P；雷达发射信号的调频斜率，记做f_dr；雷达脉冲重复频率，记做PRF；方位向等效天线长度，记做D_a；雷达接收***的采样频率，记做f_s；光在空气中的传播速度，记做C；距离向快时刻序列，记做t，t＝1,2,…,N_R， N_R为距离向快时刻总数；方位向慢时刻序列，记做l，l＝1,2,…,N_A，N_A为方位向慢时刻总数。上述参数均为SAR***标准参数，在SAR***设计和观测过程中已经确定。根据SAR成像***方案和观测方案，SAR成像方法需要的初始化***参数均为已知。 

步骤2、初始化SAR成像空间参数以及获取原始回波信号： 

初始化SAR的成像空间参数，包括：以雷达波束照射场区域地平面所构成的空间坐标作为SAR的成像空间，该成像空间记为Ω；平台到SAR成像空间中心的参考斜距，记做R_ref；将成像空间Ω均匀划分成大小相等的平面单元格，也称为分辨单元，单元网格在水平横向、水平纵向边长分别记为d_x和d_y，单元网格大小选择为SAR***传统理论成像分辨率或SAR***传统理论成像分辨率的二分之一；观测场景目标空间Ω中第m个单元格的坐标矢量，记做P_m，m表示SAR成像空间Ω中第m个单元格，m＝1,2,…,M，M为成像空间Ω中的单元格总数；SAR成像空间Ω中所有单元格的散射系数按位置顺序排列组成向量，记做α，向量α由M行1列组成；散射系数向量α中第m个元素，记做α_m，m＝1,2,…,M。 

根据步骤1中初始化的平台速度矢量V，天线初始位置矢量P(0)和雷达***的脉冲重复频率PRF，采用公式P(l)＝P(0)+V·l/PRF，l＝1,2,…,N_A，计算得到天线在第l个方位向慢时刻的位置矢量，记为P(l)，N_A为步骤1的方位向慢时刻总数。采用公式R(P(l),P_m)＝||P(l)-P_m||₂，l＝1,2,…,N_A，m＝1,2,…,M，计算得到在第l个方位向慢时刻SAR成像空间Ω中第m个单元格到天线的距离，记为R(P(l),P_m)，其中||·||₂表示定义2中的向量L2范数，P_m为初始化得到成像空间Ω中第m个单元格的坐标矢量，M为初始化的成像空间Ω中单元格总数。采用公式τ_m(l)＝2·R(P(l),P_m)/C，l＝1,2,…,N_A，m＝1,2,…,M，计算得到在第l个方位向慢时刻SAR成像空间Ω中第m个单元格到天线的时间延时，记为τ_m(l)，其中C为步骤1中初始化得到的光在空气中的传播速度。 

第l个方位向慢时刻和第t个距离向快时刻中SAR天线的原始回波数据记为s(t,l)，t＝1,2,…,N_R，l＝1,2,…,N_A，其中N_R为步骤1中初始化的距离向快时刻总数。在SAR实际成像中，回波数据s(t,l)，t＝1,2,…,N_R，l＝1,2,…,N_A，可由SAR***中雷达数据接收机提供。 

步骤3、建立SAR回波信号的线性测量模型： 

将步骤2中获取所有SAR原始回波信号s(t,l)按顺序排列组成向量，记为回波信号向量S，回波信号向量S由D行1列组成，其中D＝N_A·N_R，N_R为步骤1中初始化的距离向快时刻总数，N_A为步骤1中初始化的方位向慢时刻总数。 

采用公式φ_d(m)＝exp(-j·2·π·f_c·τ_m(l))·exp(j·π·f_dr·[t-τ_m(l)]²)，t＝1,2,…,N_R，l＝1,2,…,N_A，m＝1,2,…,M，d＝1,2,…,D，计算得到成像空间Ω中第m个单元格在回波信号向量S第d个元素信号对应的时延函数，记为φ_d(m)，其中exp(·)表示e指数运算符号，f_c为步骤1初始化得到的雷达工作中心频率，f_dr为步骤1初始化得到的发射信号调频斜率，τ_m(l)为步骤2得到的在第l个方位向慢时刻SAR成像空间Ω中第m个单元格到天线的时间延时，t为距离向的第t个快时刻，j为虚数单位(即-1的开根值)，π为圆周率。 

SAR原始回波信号向量S与成像空间Ω中所有单元格散射系数向量α之间的测量矩阵，记为A。测量矩阵A由SAR成像空间Ω中所有单元格对应的时延函数构成，A为D行M列的二维矩阵，具体表达式为 

其中，φ₁(1)为成像空间Ω中第1个单元格在回波信号向量S第1个元素信号对应的时延函数，φ₁(2)为成像空间Ω中第2个单元格在回波信号向量S第1个元素信号对应的时延函数，φ₁(M)为成像空间Ω中第M个单元格在回波信号向量S第1个元素信号对应的时延函数，φ₂(1)为成像空间Ω中第1个单元格在回波信号向量S第2个元素信号对应的时延函数，φ₂(2)为成像空间Ω中第2个单元格在回波信号向量S第2个元素信号对应的时延函数，φ₂(M)为成像空间Ω中第M个单元格在回波信号向量S第2个元素信号对应的时延函数，φ_D(1)为成像空间Ω中第1个单元格在回波信号向量S第D个元素信号对应的时延函数，φ_D(2)为成像空间Ω中第2个单元格在回波信号向量S第D个元素信号对应的时延函数，φ_D(M)为成像空间Ω中第M个单元格在回波信号向量S第D个元素信号对应的时延函数，φ₁(1),φ₁(2),…,φ₁(M)分别为成像空间Ω中第1,2,…,M个单元格在回波信号向量S第1个元素信号对应的时延函数向 量，φ₂(1),φ₂(2),…,φ₂(M)分别为成像空间Ω中第1,2,…,M个单元格在回波信号向量S第2个元素信号对应的时延函数向量，φ_D(1),φ_D(2),…,φ_D(M)分别为成像空间Ω中第1,2,…,M个单元格在回波信号向量S第D个元素信号对应的时延函数向量。 

步骤4、对SAR成像空间的散射系数向量进行初步估计： 

采用公式计算得到SAR成像空间Ω中散射系数向量的初始估计值，记做其中A^H为矩阵A的共轭转置，A为步骤3中得到的SAR测量矩阵，上标H为共轭转置运算符号，S为步骤3中得到的SAR原始回波信号向量。 

步骤5、初始化稀疏自聚焦成像算法所需的参数： 

初始化稀疏自聚焦成像算法所需参数，包括：成像处理迭代的最大迭代次数，记做Maxiter；成像处理第i次迭代过程中SAR成像模型中的相位误差向量，记为相位误差向量Φ⁽ⁱ⁾为D行1列的向量，i为自然数，表示为成像处理的第i次迭代，i＝1,2,…,Maxite，r元素为向量Φ⁽ⁱ⁾中的第1个元素，元素为向量Φ⁽ⁱ⁾中的第2个元素，元素为向量Φ⁽ⁱ⁾中的第D个元素，右上角正体符号T为转置运算符号；成像处理第i次迭代过程中第l个方位向慢时刻的相位误差值，记为l＝1,2,…,N_A；成像处理迭代的阈值，记做τ；采用公式W⁽ⁱ⁾＝diag(exp(j·Φ⁽ⁱ⁾))，i＝1,2,…,Maxiter，计算得到成像处理第i次迭代过程中SAR测量模型的初始相位误差矩阵值，记做W⁽ⁱ⁾，i＝1,2,…,Maxiter，其中相位误差矩阵W⁽ⁱ⁾为D行D列的对角矩阵，diag(·)为将向量元素作为对角矩阵中对角元素的运算符号；成像处理第i次迭代过程中SAR成像图像的图像熵，记为L⁽ⁱ⁾，i＝1,2,…,Maxiter。在成像迭代处理之前，成像处理第0次迭代中SAR测量模型相位误差矩阵的初值，记为W⁽⁰⁾，成像处理第0次迭代中SAR成像图像熵的初值，记为L⁽⁰⁾。 

步骤6、估计SAR方位向的相位误差向量： 

包括以下步骤： 

步骤6.1、获取目标回波的响应向量： 

采用合成孔径雷达标准距离压缩方法对SAR原始回波数据s(t,l)，t＝1,2,…,N_R，l＝1,2,…,N_A，进行距离压缩处理，得到距离压缩后的SAR回波数据，记做s_r(t,l)， t＝1，2，…，N_R，l＝1，2，…，N_A，其中s(t,l)为步骤2得到的第l个方位向慢时刻和第t个距离向快时刻中SAR天线的原始回波数据。 

采用公式 $\mathrm{\β}(l,m)=s_r(\mathrm{ceil}\left(\frac{R(P\left(l\right),P_m)-R_{\mathrm{ref}}}{2\·C\·f_s}\right),l)\·\exp (j\·\frac{4\·\mathrm{\π}\·R(P\left(l\right),P_m)}{\mathrm{\λ}}),$ l＝1，2，…，N_A，m＝1，2，…，M，和公式s_l＝[β(l，1)，β（l，2)，…，β（l，M）]，l＝1，2，…，N_A，计算得到第l个方位向目标回波的响应向量，记为s_l，其中β(l,1)为第l个方位向慢时刻SAR成像空间Ω中第1个单元格对应的β(l,m)，β(l,2)为第l个方位向慢时刻SAR成像空间Ω中第2个单元格对应的β(l,m)，β(l,M)为第^l个方位向慢时刻SAR成像空间Ω中第M个单元格对应的β(l,m)，ceil(·)表示上取整运算符号，R_ref为步骤2初始化得到的平台到SAR成像空间中心参考斜距，C为步骤1初始化得到的光在空气中的传播速度，f_s为步骤1初始化得到的雷达接收***的采样频率，R(P(l),P_m)为步骤2中得到的第l个方位向慢时刻SAR成像空间Ω中第m个单元格到天线的距离，λ为步骤1中初始化得到的波长。 

令q为自然数，q取值范围为q＝1,2,…,N_A。当q＝1时，采用公式 计算得到第q个方位向慢时刻的数据向量，记为X。当q＝N_A时，采用公式计算得到第q个方位向慢时刻的数据向量，记为X；当1＜q＜N_A时，采用公式 计算得到第q个方位向慢时刻的数据向量，记为X，向量X的第m个元素记为x_m，其中为成像处理第i次迭代过程中第l个方位向慢时刻的相位误差值，l＝1,2,…,N_A，s_l为第l个方位向目标回波的响应向量。 

采用公式Y＝s_q计算得到第q个方位向慢时刻的数据向量，记为Y，向量Y的第m个元素记为y_m，其中s_q为l等于q时第l个方位向目标回波的响应向量，即s_q＝s_l。 

采用公式m＝1,2,…,M，计算得到元素f_m；再采用公式F＝[f₁,f₂,…,f_M]计算得到的数据向量，记为F，其中|·|²为绝对值平方运算符号，f₁为m＝1时 对应的元素f_m，f₂为m＝2时对应的元素f_m，f_M为m＝M时对应的元素f_m，x_m为向量X的第m个元素，y_m为向量Y的第m个元素。 

步骤6.2、利用施密特正交化方法获取单位正交向量： 

采用公式和m＝1,2,…,M，计算得到元素a_m和b_m，其中x_m为步骤6.1得到的向量X中第m个元素，y_m为步骤6.1得到的向量Y中第m个元素；再采用a＝[a₁,a₂,…,a_M]、b＝[b₁,b₂,…,b_M]，m＝1,2,…,M，计算得到向量a和向量b，a₁为m＝1时对应的元素a_m，a₂为m＝2时对应的元素a_m，a_M为m＝M时对应的元素a_m，b₁为m＝1时对应的元素b_m，b₂为m＝2时对应的元素b_m，b_M为m＝M时对应的元素b_m。 

利用施密特正交化方法对向量a和向量b进行标准正交化，得到向量a和向量b构成的平面对应的单位正交向量其中为向量的第1个元素，为向量的第2个元素，为向量的第1个元素，为向量的第2个元素。 

步骤6.3、计算二次型矩阵及参数估计： 

采用公式 $r_1=\frac{{{\tilde{a}}_2}^2+{{\tilde{b}}_2}^2}{{\tilde{a}}_2\·{\tilde{b}}_1-{\tilde{a}}_1\·{\tilde{b}}_2},r_2=\frac{{{\tilde{a}}_1}^2+{{\tilde{b}}_1}^2}{{\tilde{a}}_2\·{\tilde{b}}_1-{\tilde{a}}_1\·{\tilde{b}}_2},r_3=-\frac{{\tilde{a}}_1\·{\tilde{a}}_2+{\tilde{b}}_1\·{\tilde{b}}_2}{{\tilde{a}}_2\·{\tilde{b}}_1-{\tilde{a}}_1\·{\tilde{b}}_2}$ 和 $R=\left[\begin{array}{cc}{r}_1& r_3\\ r_3& r_2\end{array}\right],$ 计算得到二次型矩阵R，其中为向量的第1个元素，为向量的第2个元素，为向量的第1个元素，为向量的第2个元素，和分别为步骤6.2得到的单位正交向量。 

采用公式Γ·Λ·Γ^T＝R对二次型矩阵R进行特征值分解，得到二次型矩阵R的特征向量矩阵记为Γ，Γ^T为矩阵Γ的转置，得到二次型矩阵R的特征值矩阵记为Λ。 

采用公式η＝(Γ·x₀-I)·R^-1计算得到参数的估计值，记为η，其中x₀为向量F在向量a和向量b所张成的二维平面的垂足点，I为单位矩阵，R^-1为二次型矩阵R的逆矩阵。 

步骤6.4、估计方位向相位误差： 

采用公式σ＝[a,b]^-1·(η·R+I)^-1·x₀计算得到的向量，记为σ，其中向量σ维数为2×1，向量σ中第1个元素记为σ₁，向量σ中第2个元素记为σ₂，a和b为步骤6.2得到的向量，η为步骤6.3得到的估计值向量，R为步骤6.3得到的二次型矩阵，x₀为向量f在向量a和b所张成的二维平面Ω内的垂足点，I为单位矩阵，上标^-1表示矩阵求逆运算符号。 

采用公式计算更新成像处理第i次迭代过程中第q个方位向慢时刻的相 位误差向量，其中atan(·)为求解正切函数反函数运算符号。将相位误差向量Φ⁽ⁱ⁾中第(q-1)·N_R+1至q·N_R个元素的值全部赋值为其中N_R为距离向快时刻总数，为相位误差向量Φ⁽ⁱ⁾中的第(q-1)·N_R+1个元素，为相位误差向量Φ⁽ⁱ⁾中的第q·N_R个元素。 

步骤6.5、方位向相位误差值逐个估计： 

对于所有方位向序号q，q＝1,2,…,N_A，采用步骤6.1到步骤6.4逐个方位向直到估计所有方位向相位误差，最终得到方位向的相位误差向量Φ⁽ⁱ⁾。 

步骤7、压缩传感稀疏成像： 

采用公式W⁽ⁱ⁾＝diag(exp(j·Φ⁽ⁱ⁾))计算得到成像处理第i次迭代过程中的相位误差矩阵，记做W⁽ⁱ⁾，其中Φ⁽ⁱ⁾为步骤6.5得到的成像处理第i次迭代过程中的相位误差向量，diag(·)为将向量元素作为对角矩阵对角元素的运算符号。 

采用公式和标准压缩传感稀疏重构方法计算得到成像处理第i次迭代过程中SAR目标成像空间的散射系数向量，记做其中表示求取满足括号中最小值时对应自变量α的最优值，向量S为步骤3得到的SAR数据回波信号向量，表示向量L2范数的平方运算符号，||·||₁表示向量L1范数运算符号。 

步骤8、计算SAR图像熵： 

采用公式m＝1,2,…,M，和公式计算得到成像处理第i次迭代过程中SAR成像图像熵，记做L⁽ⁱ⁾，i＝1,2,…,Maxiter，其中为步骤8中得到的散射系数向量的第m个元素，i为成像处理的第i次迭代，i＝1,2,…,Maxiter，M为步骤2初始化得到的SAR成像空间Ω中划分的单元格总数，为向量L2范数的平方运算符号，Σ(·)表示向量元素求和运算符号，log₂(·)表示底数为2的对数运算符号。 

步骤9、成像算法迭代条件判定： 

若满足条件||L⁽ⁱ⁾-L^(i-1)||₂≥τ且i≤Maxiter，则令成像处理迭代次数i的值加1，然后重复执行步骤6到步骤9，其中L⁽ⁱ⁾为成像处理第i次迭代过程中SAR成像图像熵，L^(i-1)为成像处 理第i-1次迭代过程中SAR成像图像熵，τ为步骤5中初始化得到的算法迭代阈值，Maxiter为步骤5中初始化得到的成像处理最大迭代次数，||·||₂为向量L2范数运算符号；若满足条件||L⁽ⁱ⁾-L^(i-1)||₂＜τ或者i＞Maxiter，则成像处理第i次迭代过程得到的散射系数向量和相位误差向量Φ⁽ⁱ⁾为最终的SAR成像散射系数向量和相位误差估计结果。 

本发明的创新点：为了获得良好的压缩传感SAR成像质量，本发明考虑了SAR回波信号测量模型中方位向相位误差对压缩传感SAR成像的影响。针对压缩传感SAR成像模型中的未知相位误差估计与补偿问题，本发明利用成像模型中方位向相位误差特性和稀疏目标特征，利用压缩传感SAR图像熵作为评价准则，提出了一种图像熵最优的压缩传感SAR稀疏自聚焦成像方法，该方法在每一次迭代处理过程中利用SAR方位向回波与观测目标的关系估计方位向相位误差，然后对压缩传感成像模型进行相位误差补偿，接着再进行压缩传感SAR成像，并采用逐次迭代方法使得压缩传感SAR成像的图像熵最优，从而提高了压缩传感SAR成像质量。 

本发明的优点在于利用图像熵最优准则实现了压缩传感SAR成像模型中未知方位向相位误差的估计与补偿问题，提高了压缩传感SAR成像质量及稳健性。 

附图说明：

图1为本发明所提供的一种图像熵最优的压缩传感SAR稀疏自聚焦成像方法流程示意图 

图2为本发明中估计SAR方位向相位误差向量方法的流程示意图 

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

具体实施方式

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

步骤1、初始化SAR***参数： 

初始化SAR***参数包括：平台速度矢量V＝[0,150,0]m/s；天线初始位置矢量P(0)＝[0,0,6000]m；雷达工作中心频率f_c＝10×10⁹Hz；雷达载频波长λ＝0.03m；雷达发射基带信号的信号带宽B_r＝1.5×10⁸Hz；雷达发射信号脉冲宽度T_P＝1×10^-6s；雷达发射信号的调频斜率f_dr＝1.5×10¹⁴Hz/s；雷达接收***的采样频率f_s＝3×10⁸Hz；雷达脉冲重复频率PRF＝500Hz；方位向等效天线长度为D_a＝2m；光在空气中的传播速度C＝3×10⁸m/s；距 离向快时刻总数N_R＝128，距离向快时刻序列t＝1,2,…,N_R；方位向慢时刻总数N_A＝128，方位向慢时刻序列l＝1,2,…,N_A。上述参数均为SAR***标准参数，在SAR***设计和观测过程中已经确定。根据SAR成像***方案和观测方案，SAR成像方法需要的初始化***参数均为已知。 

步骤2、初始化SAR成像空间参数以及获取原始回波信号： 

初始化SAR的成像空间参数，包括：以雷达波束照射场区域地平面所构成的空间坐标作为SAR的成像空间，该成像空间记为Ω；初始化SAR成像空间Ω的大小为128×128像素，SAR成像空间Ω的中心坐标位置位于[8000,0,0]m，每一个单元网格在水平横向和水平纵向边长为d_x＝d_y＝0.5m，计算得到SAR成像空间的总单元格数M＝16384，SAR成像空间Ω中每一个单元格的位置为P_m＝[(x′-64)·0.5,(y′-64)·0.5,0]m，其中x′＝1,…,2，y′＝1,2,…,128，m＝(x′-1)·128+y′。场景中心参考斜距R_ref＝10000m。P_m为SAR成像空间Ω中第m个单元格的位置矢量，m表示SAR成像空间Ω中第m个单元格，m＝1,2,…,M，M＝16384。在SAR成像空间Ω中加入仿真点目标散射体，点目标散射体数个数为5个，它们散射系数值均为1，坐标位置分别为[0,0,]、[20,20,0]、[20,-21,0]、[-20,20,0]、[-20,-21,0]，单位均为m；SAR成像空间Ω中没有包含点目标单元格的散射系数设置为0。将SAR成像空间Ω中所有单元格的目标散射系数按单元格位置顺序排列组成散射向量α。确定SAR成像空间Ω所有单元散射系数后，散射系数向量α在SAR成像观测仿真过程中就已经确定。场景目标散射系数向量α由M行1列组成，α_m为向量α中对应SAR成像空间Ω中第m个单元格的散射系数值。在本仿真成像空间Ω中，只有包含点散射目标的5个单元格散射系数值α设置为1，其余单元格的散射系数都为0。利用传统的合成孔径雷达原始回波仿真方法产生SAR的原始回波数据。 

根据步骤1中初始化得到的平台速度矢量V＝[0,150,0]m/s，天线的初始位置矢量P(0)＝[0,0,6000]m和脉冲重复频率PRF＝500Hz，采用公式P(l)＝P(0)+V·l/PRF，l＝1,2,…,N_A，计算得到天线在第l个方位向慢时刻的位置矢量，记为P(l)，N_A为步骤1的方位向慢时刻总数N_A＝128。采用公式R(P(l),P_m)＝||P(l)-P_m||₂，l＝1,2,…,N_A，m＝1,2,…,M，M＝16384，计算得到在第l个方位向慢时刻SAR成像空间Ω中第m个单元 格到天线的距离，记为R(P(l),P_m)，其中||·||₂表示定义2中的向量L2范数，P_m为初始化得到成像空间Ω中第m个单元格的坐标矢量，M为初始化的成像空间Ω中单元格总数M＝16384。采用公式τ_m(l)＝2·R(P(l),P_m)/C，l＝1,2,…,N_A，m＝1,2,…,M，N_A＝128，M＝16384，计算得到在第l个方位向慢时刻SAR成像空间Ω中第m个单元格到天线的时间延时，记为τ_m(l)，其中C＝3×10⁸m/s。 

在仿真过程或者实际成像过程中获取SAR原始回波数据，在第l个方位向慢时刻和第t个距离向快时刻中SAR天线的原始回波数据记为s(t,l)，t＝1,2,…,N_R，l＝1,2,…,N_A，其中N_R为步骤1中初始化的距离向快时刻总数N_R＝128，N_A为步骤1中初始化的方位向慢时刻总数N_A＝128。在SAR实际成像中，s(t,l)可由SAR***中雷达数据接收机提供；而在仿真过程中，s(t,l)采用定义6中传统合成孔径雷达原始回波仿真方法产生提供。 

步骤3、建立SAR回波信号的线性测量模型： 

将步骤2中获取的SAR原始回波信号s(t,l)按顺序排列组成回波信号向量S，回波信号向量S由D行1列组成，其中D＝N_A·N_R＝16384，N_R为步骤1中初始化的距离向快时刻总数N_R＝128，N_A为步骤1中初始化的方位向慢时刻总数N_A＝128。 

采用公式φ_d(m)＝exp(-j·2·π·f_c·τ_m(l))·exp(j·π·f_dr·[t-τ_m(l)]2)，t＝1，2，...，N_R，l＝1,2,…,N_A，m＝1,2,…,M，d＝1,2,…,D，N_R＝128，N_A＝128，M＝16384，D＝16384，计算得到成像空间Ω中第m个单元格在回波信号向量S第d个元素信号对应的时延函数，记为φ_d(m)，其中exp(·)表示e指数运算符号，f_c为步骤1初始化得到的雷达工作中心频率，f_dr为步骤1初始化得到的发射信号调频斜率，τ_m(l)为步骤2得到的在第l个方位向慢时刻SAR成像空间Ω中第m个单元格到天线的时间延时，t为距离向的第t个快时刻，j为虚数单位(即-1的开根值)，π为圆周率，约为π＝3.1416。 

采用矩阵表达公式 

计算得到SAR原始回波信号与SAR成像空间Ω所有单元格的测量矩阵A，测量矩阵A为D行M列的二维矩阵。 

步骤4、对SAR成像空间的散射系数向量进行初步估计： 

采用公式计算得到SAR成像空间Ω的散射系数向量的初始估计值，记做其中A^H为矩阵A的共轭转置，A为步骤3中得到的SAR测量矩阵，上标H为共轭转置运算符号，S为步骤3中得到的SAR原始回波信号向量。 

步骤5、初始化稀疏自聚焦成像算法所需的参数： 

初始化稀疏自聚焦成像算法所需参数，包括：成像迭代处理的最大迭代次数Maxiter＝15；成像处理第i次迭代过程中SAR成像模型中的相位误差向量的初值为零向量，相位误差向量Φ⁽ⁱ⁾为D行1列的向量，i为自然数，表示为成像处理的第i次迭代，i＝1,2,…,Maxiter，元素为向量Φ⁽ⁱ⁾中的第1个元素，元素为向量Φ⁽ⁱ⁾中的第2个元素，元素为向量Φ⁽ⁱ⁾中的第D个元素，右上角正体符号T为转置运算符号；成像处理第i次迭代过程中第l个方位向慢时刻的相位误差值，记为l＝1,2,…,N_A，N_A＝128；成像处理迭代的阈值τ＝0.01；采用公式W⁽ⁱ⁾＝diag(exp(j·Φ⁽ⁱ⁾))，i＝1,2,…,Maxiter，计算得到成像处理第i次迭代过程中SAR测量模型的初始相位误差矩阵值，记做W⁽ⁱ⁾，i＝1,2,…,Maxiter，Maxiter＝15，其中相位误差矩阵W⁽ⁱ⁾为D行D列的对角矩阵，diag(·)为将向量元素作为对角矩阵中对角元素的运算符号；成像处理第i次迭代过程中SAR成像图像的图像熵，记为L⁽ⁱ⁾，i＝1,2,…,Maxiter，Maxiter＝15。在成像处理迭代之前，成像处理第0次迭代中SAR测量模型相位误差矩阵W⁽⁰⁾的初值为零矩阵，成像处理第0次迭代中SAR成像图像熵的初值L⁽⁰⁾＝0。 

步骤6、估计SAR方位向的相位误差向量： 

包括以下步骤： 

步骤6.1、获取目标回波的响应向量： 

采用标准合成孔径雷达标准距离压缩方法对SAR原始回波数据s(t,l)，t＝1,2,…,N_R，l＝1,2,…,N_A，进行距离压缩处理，得到距离压缩后的SAR回波数据，记做s_r(t,l)，t＝1,2,…,N_R，l＝1,2,…,N_A，N_R＝128，N_A＝128，其中s(t,l)为步骤2得到的第l个方位向慢时刻和第t个距离向快时刻中SAR天线的原始回波数据。 

采用公式 $\mathrm{\β}(l,m)=s_r(\mathrm{ceil}\left(\frac{R(P\left(l\right),P_m)-R_{\mathrm{ref}}}{2\·C\·f_s}\right),l)\·\exp (j\·\frac{4\·\mathrm{\π}\·R(P\left(l\right),P_m)}{\mathrm{\λ}}),$ l＝1,2,…,N_A，m＝1,2,…,M，和公式s_l＝[β(l,1),β(l,2),…,β(l,M)]，l＝1,2,…,N_A，计算得到第l个方位向目标回波的响应向量，记为s_l，其中β(l,1)为第l个方位向慢时刻SAR成像空间Ω中第1个单元格对应的β(l,m)，β(l,2)为第l个方位向慢时刻SAR成像空间Ω中第2个单元格对应的β(l,m)，β(l,M)为第l个方位向慢时刻SAR成像空间Ω中第M个单元格对应的β(l,m)，ceil(·)表示上取整运算符号，R_ref为步骤2初始化得到的平台到SAR成像空间中心参考斜距，C为步骤1初始化得到的光在空气中的传播速度，f_s为步骤1初始化得到的雷达接收***的采样频率，R(P(l),P_m)为步骤2中得到的第l个方位向慢时刻SAR成像空间Ω中第m个单元格到天线的距离，λ为步骤1中初始化得到的波长。 

令q为自然数，q取值范围为q＝1,2,…,N_A，N_A＝128。当q＝1时，采用公式 计算得到第q个方位向慢时刻的数据向量，记为X。当q＝N_A时，采用公式计算得到第q个方位向慢时刻的数据向量，记为X；当1＜q＜N_A时，采用公式 计算得到第q个方位向慢时刻的数据向量，记为X，向量X的第m个元素记为x_m，其中为成像处理第i次迭代过程中第l个方位向慢时刻的相位误差值，l＝1,2,…,N_A，s_l为第l个方位向目标回波的响应向量。 

采用公式Y＝s_q计算得到第q个方位向慢时刻的数据向量，记为Y，向量Y的第m个元素记为ym，其中s_q为l等于q时第l个方位向目标回波的响应向量，即s_q＝s_l。 

采用公式m＝1,2,…,M，计算得到元素f_m；再采用公式F＝[f₁,f₂,…,f_M]计算得到的数据向量，记为F，其中|·|²为绝对值平方运算符号，f₁为m＝1时对应的元素f_m，f₂为m＝2时对应的元素f_m，f_M为m＝M时对应的元素f_m，x_m为向量X的第m个元素，y_m为向量Y的第m个元素。 

步骤6.2、利用施密特正交化方法获取单位正交向量： 

采用公式和m＝1,2,…,M，M＝16384，计算得到元素a_m和b_m，其中x_m为步骤6.1得到的向量X中第m个元素，y_m为步骤6.1得到的向量Y中第m个元素；再采用a＝[a₁,a₂,…,a_M]、b＝[b₁,b₂,…,b_M]，m＝1,2,…,M，M＝16384，计算得到向量a和向量b，a₁为m＝1时对应的元素a_m，a₂为m＝2时对应的元素a_m，a_M为m＝M时对应的元素a_m，b₁为m＝1时对应的元素b_m，b₂为m＝2时对应的元素b_m，b_M为m＝M时对应的元素b_m。 

利用施密特正交化方法对向量a和向量b进行标准正交化，得到向量a和向量b构成的平面对应的单位正交向量其中为向量的第1个元素，为向量的第2个元素，为向量的第1个元素，为向量的第2个元素。 

步骤6.3、计算二次型矩阵及参数估计： 

采用公式 $r_1=\frac{{{\tilde{a}}_2}^2+{{\tilde{b}}_2}^2}{{\tilde{a}}_2\·{\tilde{b}}_1-{\tilde{a}}_1\·{\tilde{b}}_2},r_2=\frac{{{\tilde{a}}_1}^2+{{\tilde{b}}_1}^2}{{\tilde{a}}_2\·{\tilde{b}}_1-{\tilde{a}}_1\·{\tilde{b}}_2},r_3=-\frac{{\tilde{a}}_1\·{\tilde{a}}_2+{\tilde{b}}_1\·{\tilde{b}}_2}{{\tilde{a}}_2\·{\tilde{b}}_1-{\tilde{a}}_1\·{\tilde{b}}_2}$ 和 $R=\left[\begin{array}{cc}{r}_1& r_3\\ r_3& r_2\end{array}\right],$ 计算得到二次型矩阵R，其中为向量的第1个元素，为向量的第2个元素，为向量的第1个元素，为向量的第2个元素，和分别为步骤6.2得到的单位正交向量。 

采用公式Γ·Λ·Γ^T＝R对二次型矩阵R进行特征值分解，得到二次型矩阵R的特征向量矩阵记为Γ，Γ^T为矩阵Γ的转置，得到二次型矩阵R的特征值矩阵记为Λ。 

采用公式η＝(Γ·x₀-I)·R^-1计算得到参数的估计值，记为η，其中x₀为向量F在向量a和向量b所张成的二维平面的垂足点，I为单位矩阵，R^-1为二次型矩阵R的逆矩阵。 

步骤6.4、估计方位向相位误差： 

采用公式σ＝[a,b]^-1·(η·R+I)^-1·x₀计算得到的向量，记为σ，其中向量σ维数为2×1，向量σ中第1个元素记为σ₁，向量σ中第2个元素记为σ₂，a和b为步骤6.2得到的向量，η为步骤6.3得到的估计值向量，R为步骤6.3得到的二次型矩阵，x₀为向量f在向量a和b所张成的二维平面Ω内的垂足点，I为单位矩阵，上标^-1表示矩阵求逆运算符号。 

采用公式计算更新成像处理第i次迭代过程中第q个方位向慢时刻的相位误差向量，其中atan(·)为求解正切函数反函数运算符号。将相位误差向量Φ⁽ⁱ⁾中第 (q-1)·N_R+1至q·N_R个元素的值全部赋值为其中N_R为距离向快时刻总数，为相位误差向量Φ⁽ⁱ⁾中的第(q-1)·N_R+1个元素，为相位误差向量Φ⁽ⁱ⁾中的第q·N_R个元素。 

步骤6.5、方位向相位误差值逐个估计： 

对于所有方位向序号q，q＝1,2,…,N_A，N_A＝128，采用步骤6.1到步骤6.4逐个方位向直到估计所有方位向相位误差，最终得到方位向的相位误差向量Φ⁽ⁱ⁾。 

步骤7、压缩传感稀疏成像： 

采用公式W⁽ⁱ⁾＝diag(exp(j·Φ⁽ⁱ⁾))计算得到成像处理第i次迭代过程中的相位误差矩阵，记做W⁽ⁱ⁾，其中Φ⁽ⁱ⁾为步骤6.5得到的成像处理第i次迭代过程中的相位误差向量，diag(·)为将向量元素作为对角矩阵对角元素的运算符号。 

采用公式和标准压缩传感稀疏重构方法计算得到成像处理第i次迭代过程中SAR目标成像空间的散射系数向量，记做其中表示求取满足括号中最小值时对应自变量α的最优值，向量S为步骤3得到的SAR数据回波信号向量，表示向量L2范数的平方运算符号，||·||₁表示向量L1范数运算符号。 

步骤8、计算SAR图像熵： 

采用公式m＝1,2,…,M，M＝16384，和公式计算得到成像处理第i次迭代过程中SAR成像图像熵，记做L⁽ⁱ⁾，i＝1,2,…,Maxiter，Maxiter＝15，其中为步骤8中得到的散射系数向量的第m个元素，i为算法的第i次迭代，i＝1,2,…,15，M为步骤2初始化得到的SAR成像空间Ω中划分的单元格总数M＝16384，为向量L2范数的平方运算符号，Σ(·)表示向量元素求和运算符号，log₂(·)表示底数为2的对数运算符号。 

步骤9、成像处理迭代条件判定： 

若满足条件||L⁽ⁱ⁾-L^(i-1)||₂≥τ且i≤Maxiter，则令成像处理迭代次数i的值加1，然后重复执行步骤6到步骤9，其中L⁽ⁱ⁾为成像处理第i次迭代过程中SAR成像图像熵，L^(i-1)为成像处 理第i-1次迭代过程中SAR成像图像熵，τ为步骤5中初始化得到的算法迭代阈值τ＝0.01，Maxiter为步骤5中初始化得到的成像处理最大迭代次数Maxiter＝15，||·||₂为向量L2范数运算符号；若满足条件||L⁽ⁱ⁾-L^(i-1)||₂＜τ或者i＞Maxiter，则成像处理第i次迭代过程得到的散射系数向量和相位误差向量Φ⁽ⁱ⁾为最终的SAR成像散射系数向量和相位误差估计结果。 

Claims (1)

1.一种图像熵最优的SAR稀疏自聚焦成像方法，其特征是它包括如下步骤：

步骤1、初始化SAR***参数：

初始化SAR***参数包括：平台速度矢量，记做V；天线初始位置矢量，即方位向0慢时刻天线位置，记做P(0)；雷达工作中心频率，记做f_c；雷达载频波长，记做λ；雷达发射基带信号的信号带宽，记做B_r；雷达发射信号脉冲宽度，记做T_P；雷达发射信号的调频斜率，记做f_dr；雷达脉冲重复频率，记做PRF；方位向等效天线长度，记做D_a；雷达接收***的采样频率，记做f_s；光在空气中的传播速度，记做C；距离向快时刻序列，记做t，t＝1,2,…,N_R，N_R为距离向快时刻总数；方位向慢时刻序列，记做l，l＝1,2,…,N_A，N_A为方位向慢时刻总数；上述参数均为SAR***标准参数，在SAR***设计和观测过程中已经确定；根据SAR成像***方案和观测方案，SAR成像方法需要的初始化***参数均为已知；

步骤2、初始化SAR成像空间参数以及获取原始回波信号：

初始化SAR的成像空间参数，包括：以雷达波束照射场区域地平面所构成的空间坐标作为SAR的成像空间，该成像空间记为Ω；平台到SAR成像空间中心的参考斜距，记做R_ref；将成像空间Ω均匀划分成大小相等的平面单元格，也称为分辨单元，单元网格在水平横向、水平纵向边长分别记为d_x和d_y，单元网格大小选择为SAR***传统理论成像分辨率或SAR***传统理论成像分辨率的二分之一；观测场景目标空间Ω中第m个单元格的坐标矢量，记做P_m，m表示SAR成像空间Ω中第m个单元格，m＝1,2,…,M，M为成像空间Ω中的单元格总数；SAR成像空间Ω中所有单元格的散射系数按位置顺序排列组成向量，记做α，向量α由M行1列组成；散射系数向量α中第m个元素，记做α_m，m＝1,2,…,M；

根据步骤1中初始化的平台速度矢量V，天线初始位置矢量P(0)和雷达***的脉冲重复频率PRF，采用公式P(l)＝P(0)+V·l/PRF，l＝1,2,…,N_A，计算得到天线在第l个方位向慢时刻的位置矢量，记为P(l)，N_A为步骤1的方位向慢时刻总数；采用公式R(P(l),P_m)＝||P(l)-P_m||₂，l＝1,2,…,N_A，m＝1,2,…,M，计算得到在第l个方位向慢时刻SAR成像空间Ω中第m个单元格到天线的距离，记为R(P(l),P_m)，其中||·||₂表示定义2中的向量L2范数，P_m为初始化得到成像空间Ω中第m个单元格的坐标矢量，M为初始化的成像空间Ω中单元格总数；采用公式τ_m(l)＝2·R(P(l),P_m)/C，l＝1,2,…,N_A，m＝1,2,…,M，计算得到在第l个方位向慢时刻SAR成像空间Ω中第m个单元格到天线的时间延时，记为τ_m(l)，其中C为步骤1中初始化得到的光在空气中的传播速度；

第l个方位向慢时刻和第t个距离向快时刻中SAR天线的原始回波数据记为s(t,l)，t＝1,2,…,N_R，l＝1,2,…,N_A，其中N_R为步骤1中初始化的距离向快时刻总数；在SAR实际成像中，回波数据s(t,l)，t＝1,2,…,N_R，l＝1,2,…,N_A，可由SAR***中雷达数据接收机提供；

步骤3、建立SAR回波信号的线性测量模型：

将步骤2中获取所有SAR原始回波信号s(t,l)按顺序排列组成向量，记为回波信号向量S，回波信号向量S由D行1列组成，其中D＝N_A·N_R，N_R为步骤1中初始化的距离向快时刻总数，N_A为步骤1中初始化的方位向慢时刻总数；

采用公式t＝1,2,…,N_R，l＝1,2,…,N_A，m＝1,2,…,M，d＝1,2,…,D，计算得到成像空间Ω中第m个单元格在回波信号向量S第d个元素信号对应的时延函数，记为φ_d(m)，其中exp(·)表示e指数运算符号，f_c为步骤1初始化得到的雷达工作中心频率，f_dr为步骤1初始化得到的发射信号调频斜率，τ_m(l)为步骤2得到的在第l个方位向慢时刻SAR成像空间Ω中第m个单元格到天线的时间延时，t为距离向的第t个快时刻，j为虚数单位(即-1的开根值)，π为圆周率；

SAR原始回波信号向量S与成像空间Ω中所有单元格散射系数向量α之间的测量矩阵，记为A；测量矩阵A由SAR成像空间Ω中所有单元格对应的时延函数构成，A为D行M列的二维矩阵，具体表达式为

其中，φ₁(1)为成像空间Ω中第1个单元格在回波信号向量S第1个元素信号对应的时延函数，φ₁(2)为成像空间Ω中第2个单元格在回波信号向量S第1个元素信号对应的时延函数，φ₁(M)为成像空间Ω中第M个单元格在回波信号向量S第1个元素信号对应的时延函数，φ₂(1)为成像空间Ω中第1个单元格在回波信号向量S第2个元素信号对应的时延函数，φ₂(2)为成像空间Ω中第2个单元格在回波信号向量S第2个元素信号对应的时延函数，φ₂(M)为成像空间Ω中第M个单元格在回波信号向量S第2个元素信号对应的时延函数，φ_D(1)为成像空间Ω中第1个单元格在回波信号向量S第D个元素信号对应的时延函数，φ_D(2)为成像空间Ω中第2个单元格在回波信号向量S第D个元素信号对应的时延函数，φ_D(M)为成像空间Ω中第M个单元格在回波信号向量S第D个元素信号对应的时延函数，φ₁(1),φ₁(2),…,φ₁(M)分别为成像空间Ω中第1,2,…,M个单元格在回波信号向量S第1个元素信号对应的时延函数向量，φ₂(1),φ₂(2),…,φ₂(M)分别为成像空间Ω中第1,2,…,M个单元格在回波信号向量S第2个元素信号对应的时延函数向量，φ_D(1),φ_D(2),…,φ_D(M)分别为成像空间Ω中第1,2,…,M个单元格在回波信号向量S第D个元素信号对应的时延函数向量；

步骤4、对SAR成像空间的散射系数向量进行初步估计：

采用公式计算得到SAR成像空间Ω中散射系数向量的初始估计值，记做其中A^H为矩阵A的共轭转置，A为步骤3中得到的SAR测量矩阵，上标H为共轭转置运算符号，S为步骤3中得到的SAR原始回波信号向量；

步骤5、初始化稀疏自聚焦成像算法所需的参数：

初始化稀疏自聚焦成像算法所需参数，包括：成像处理迭代的最大迭代次数，记做Maxiter；成像处理第i次迭代过程中SAR成像模型中的相位误差向量，记为相位误差向量Φ⁽ⁱ⁾为D行1列的向量，i为自然数，表示为成像处理的第i次迭代，i＝1,2,…,Maxiter，元素为向量Φ⁽ⁱ⁾中的第1个元素，元素为向量Φ⁽ⁱ⁾中的第2个元素，元素为向量Φ⁽ⁱ⁾中的第D个元素，右上角正体符号T为转置运算符号；成像处理第i次迭代过程中第l个方位向慢时刻的相位误差值，记为l＝1,2,…,N_A；成像处理迭代的阈值，记做τ；采用公式W⁽ⁱ⁾＝diag(exp(j·Φ⁽ⁱ⁾))，i＝1,2,…,Maxiter，计算得到成像处理第i次迭代过程中SAR测量模型的初始相位误差矩阵值，记做W⁽ⁱ⁾，i＝1,2,…,Maxiter，其中相位误差矩阵W⁽ⁱ⁾为D行D列的对角矩阵，diag(·)为将向量元素作为对角矩阵中对角元素的运算符号；成像处理第i次迭代过程中SAR成像图像的图像熵，记为L⁽ⁱ⁾，i＝1,2,…,Maxiter；在成像迭代处理之前，成像处理第0次迭代中SAR测量模型相位误差矩阵的初值，记为W⁽⁰⁾，成像处理第0次迭代中SAR成像图像熵的初值，记为L⁽⁰⁾；

步骤6、估计SAR方位向的相位误差向量：

包括以下步骤：

步骤6.1、获取目标回波的响应向量：

采用合成孔径雷达标准距离压缩方法对SAR原始回波数据s(t,l)，t＝1,2,…,N_R，l＝1,2,…,N_A，进行距离压缩处理，得到距离压缩后的SAR回波数据，记做s_r(t,l)，t＝1,2,…,N_R，l＝1,2,…,N_A，其中s(t,l)为步骤2得到的第l个方位向慢时刻和第t个距离向快时刻中SAR天线的原始回波数据；

采用公式 $\mathrm{\β}(l,m)=s_r(\mathrm{ceil}\left(\frac{R(P\left(l\right),P_m)-R_{\mathrm{ref}}}{2\·C\·f_s}\right),l)\·\exp (j\·\frac{4\·\mathrm{\π}\·R(P\left(l\right),P_m)}{\mathrm{\λ}}),$ l＝1,2,…,N_A，m＝1,2,…,M，和公式s_l＝[β(l,1),β(l,2),…,β(l,M)]，l＝1,2,…,N_A，计算得到第l个方位向目标回波的响应向量，记为s_l，其中β(l,1)为第l个方位向慢时刻SAR成像空间Ω中第1个单元格对应的β(l,m)，β(l,2)为第l个方位向慢时刻SAR成像空间Ω中第2个单元格对应的β(l,m)，β(l,M)为第l个方位向慢时刻SAR成像空间Ω中第M个单元格对应的β(l,m)，ceil(·)表示上取整运算符号，R_ref为步骤2初始化得到的平台到SAR成像空间中心参考斜距，C为步骤1初始化得到的光在空气中的传播速度，f_s为步骤1初始化得到的雷达接收***的采样频率，R(P(l),P_m)为步骤2中得到的第l个方位向慢时刻SAR成像空间Ω中第m个单元格到天线的距离，λ为步骤1中初始化得到的波长；

令q为自然数，q取值范围为q＝1,2,…,N_A；当q＝1时，采用公式计算得到第q个方位向慢时刻的数据向量，记为X；当q＝N_A时，采用公式计算得到第q个方位向慢时刻的数据向量，记为X；当1＜q＜N_A时，采用公式计算得到第q个方位向慢时刻的数据向量，记为X，向量X的第m个元素记为x_m，其中为成像处理第i次迭代过程中第l个方位向慢时刻的相位误差值，l＝1,2,…,N_A，s_l为第l个方位向目标回波的响应向量；

采用公式Y＝s_q计算得到第q个方位向慢时刻的数据向量，记为Y，向量Y的第m个元素记为y_m，其中s_q为l等于q时第l个方位向目标回波的响应向量，即s_q＝s_l；

采用公式m＝1,2,…,M，计算得到元素f_m；再采用公式F＝[f₁,f₂,…,f_M]计算得到的数据向量，记为F，其中|·|²为绝对值平方运算符号，f₁为m＝1时对应的元素f_m，f₂为m＝2时对应的元素f_m，f_M为m＝M时对应的元素f_m，x_m为向量X的第m个元素，y_m为向量Y的第m个元素；

步骤6.2、利用施密特正交化方法获取单位正交向量：

采用公式和m＝1,2,…,M，计算得到元素a_m和b_m，其中x_m为步骤6.1得到的向量X中第m个元素，y_m为步骤6.1得到的向量Y中第m个元素；再采用a＝[a₁,a₂,…,a_M]、b＝[b₁,b₂,…，b_M]，m＝1,2,…,M，计算得到向量a和向量b，a₁为m＝1时对应的元素a_m，a₂为m＝2时对应的元素a_m，a_M为m＝M时对应的元素a_m，b₁为m＝1时对应的元素b_m，b₂为m＝2时对应的元素b_m，b_M为m＝M时对应的元素b_m；

利用施密特正交化方法对向量a和向量b进行标准正交化，得到向量a和向量b构成的平面对应的单位正交向量其中为向量的第1个元素，为向量的第2个元素，为向量的第1个元素，为向量的第2个元素；

步骤6.3、计算二次型矩阵及参数估计：

采用公式 $r_1=\frac{{{\tilde{a}}_2}^2+{{\tilde{b}}_2}^2}{{\tilde{a}}_2\·{\tilde{b}}_1-{\tilde{a}}_1\·{\tilde{b}}_2},r_2=\frac{{{\tilde{a}}_1}^2+{{\tilde{b}}_1}^2}{{\tilde{a}}_2\·{\tilde{b}}_1-{\tilde{a}}_1\·{\tilde{b}}_2},r_3=-\frac{{\tilde{a}}_1\·{\tilde{a}}_2+{\tilde{b}}_1\·{\tilde{b}}_2}{{\tilde{a}}_2\·{\tilde{b}}_1-{\tilde{a}}_1\·{\tilde{b}}_2}$ 和 $R=\left[\begin{array}{cc}{r}_1& r_3\\ r_3& r_2\end{array}\right],$ 计算得到二次型矩阵R，其中为向量的第1个元素，为向量的第2个元素，为向量的第1个元素，为向量的第2个元素，和分别为步骤6.2得到的单位正交向量；

采用公式Γ·Λ·Γ^T＝R对二次型矩阵R进行特征值分解，得到二次型矩阵R的特征向量矩阵记为Γ，Γ^T为矩阵Γ的转置，得到二次型矩阵R的特征值矩阵记为Λ；

采用公式η＝(Γ·x₀-I)·R^-1计算得到参数的估计值，记为η，其中x₀为向量F在向量a和向量b所张成的二维平面的垂足点，I为单位矩阵，R^-1为二次型矩阵R的逆矩阵；

步骤6.4、估计方位向相位误差：

采用公式σ＝[a,b]^-1·(η·R+I)^-1·x₀计算得到的向量，记为σ，其中向量σ维数为2×1，向量σ中第1个元素记为σ₁，向量σ中第2个元素记为σ₂，a和b为步骤6.2得到的向量，η为步骤6.3得到的估计值向量，R为步骤6.3得到的二次型矩阵，x₀为向量f在向量a和b所张成的二维平面Ω内的垂足点，I为单位矩阵，上标^-1表示矩阵求逆运算符号；

采用公式计算更新成像处理第i次迭代过程中第q个方位向慢时刻的相位误差向量，其中atan(·)为求解正切函数反函数运算符号；将相位误差向量Φ⁽ⁱ⁾中第(q-1)·N_R+1至q·N_R个元素的值全部赋值为其中N_R为距离向快时刻总数，为相位误差向量Φ⁽ⁱ⁾中的第(q-1)·N_R+1个元素，为相位误差向量Φ⁽ⁱ⁾中的第q·N_R个元素；

步骤6.5、方位向相位误差值逐个估计：

对于所有方位向序号q，q＝1,2,…,N_A，采用步骤6.1到步骤6.4逐个方位向直到估计所有方位向相位误差，最终得到方位向的相位误差向量Φ⁽ⁱ⁾；

步骤7、压缩传感稀疏成像：

采用公式W⁽ⁱ⁾＝diag(exp(j·Φ⁽ⁱ⁾))计算得到成像处理第i次迭代过程中的相位误差矩阵，记做W⁽ⁱ⁾，其中Φ⁽ⁱ⁾为步骤6.5得到的成像处理第i次迭代过程中的相位误差向量，diag(·)为将向量元素作为对角矩阵对角元素的运算符号；

采用公式和标准压缩传感稀疏重构方法计算得到成像处理第i次迭代过程中SAR目标成像空间的散射系数向量，记做其中表示求取满足括号中最小值时对应自变量α的最优值，向量S为步骤3得到的SAR数据回波信号向量，表示向量L2范数的平方运算符号，||·||₁表示向量L1范数运算符号；

步骤8、计算SAR图像熵：

采用公式m＝1,2,…,M，和公式计算得到成像处理第i次迭代过程中SAR成像图像熵，记做L⁽ⁱ⁾，i＝1,2,…,Maxiter，其中为步骤8中得到的散射系数向量的第m个元素，i为成像处理的第i次迭代，i＝1,2,…,Maxiter，M为步骤2初始化得到的SAR成像空间Ω中划分的单元格总数，为向量L2范数的平方运算符号，Σ(·)表示向量元素求和运算符号，log₂(·)表示底数为2的对数运算符号；

步骤9、成像算法迭代条件判定：

若满足条件||L⁽ⁱ⁾-L^(i-1)||₂≥τ且i≤Maxiter，则令成像处理迭代次数i的值加1，然后重复执行步骤6到步骤9，其中L⁽ⁱ⁾为成像处理第i次迭代过程中SAR成像图像熵，L^(i-1)为成像处理第i-1次迭代过程中SAR成像图像熵，τ为步骤5中初始化得到的算法迭代阈值，Maxiter为步骤5中初始化得到的成像处理最大迭代次数，||·||₂为向量L2范数运算符号；

若满足条件||L⁽ⁱ⁾-L^(i-1)||₂＜τ或者i＞Maxiter，则成像处理第i次迭代过程得到的散射系数向量和相位误差向量Φ⁽ⁱ⁾为最终的SAR成像散射系数向量和相位误差估计结果。