CN104833973B - 基于半正定规划的线阵sar后向投影自聚焦成像方法 - Google Patents

基于半正定规划的线阵sar后向投影自聚焦成像方法 Download PDF

Info

Publication number
CN104833973B
CN104833973B CN201510231100.XA CN201510231100A CN104833973B CN 104833973 B CN104833973 B CN 104833973B CN 201510231100 A CN201510231100 A CN 201510231100A CN 104833973 B CN104833973 B CN 104833973B
Authority
CN
China
Prior art keywords
linear array
vector
designated
matrix
sar
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CN201510231100.XA
Other languages
English (en)
Other versions
CN104833973A (zh
Inventor
韦顺军
张晓玲
吴文俊
刘嘉宁
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
University of Electronic Science and Technology of China
Original Assignee
University of Electronic Science and Technology of China
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by University of Electronic Science and Technology of China filed Critical University of Electronic Science and Technology of China
Priority to CN201510231100.XA priority Critical patent/CN104833973B/zh
Publication of CN104833973A publication Critical patent/CN104833973A/zh
Application granted granted Critical
Publication of CN104833973B publication Critical patent/CN104833973B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9004—SAR image acquisition techniques
    • G01S13/9019—Auto-focussing of the SAR signals
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9004—SAR image acquisition techniques
    • G01S13/9017—SAR image acquisition techniques with time domain processing of the SAR signals in azimuth
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/904—SAR modes

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • Electromagnetism (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Signal Processing (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明公开了基于半正定规划的线阵SAR后向投影自聚焦成像方法，它通过构造线阵SAR后向投影积累基函数矩阵，采用凸优化理论中的半正定规划方法估计线阵SAR回波数据的相位误差矩阵，利用后向投影积累基函数矩阵与相位误差矩阵的乘积矩阵的最小秩，求解约束条件下的最优矩阵实现线阵SAR相位误差的估计和自聚焦成像，从而将线阵SAR后向投影算法成像过程转化为矩阵与向量相乘的形式，得到线阵SAR成像结果。

Description

基于半正定规划的线阵SAR后向投影自聚焦成像方法
技术领域：
本技术发明属于雷达技术领域，它特别涉及了合成孔径雷达(SAR)成像技术领域。
背景技术：
传统合成孔径雷达(SAR)利用单天线的直线运动合成一维的虚拟线阵天线，获得方位向高分辨率，再利用脉冲压缩技术获得雷达视线方向高分辨率，从而实现观测场景的二维成像。但传统二维SAR主要工作于侧视成像模式，而侧视成像时存在地形遮挡、阴影效应及顶底倒置问题，故传统二维SAR在城市、山地和峡谷等复杂起伏地形情况不能获得满意的成像结果。三维SAR基本原理是通过天线运动合成虚拟二维面阵天线，获得面阵平面内二维高分辨，再结合脉冲压缩技术获得雷达视线方向高分辨率，实现对观测目标三维成像，克服了传统二维SAR成像技术在复杂起伏地形区域的缺陷。三维SAR是SAR成像技术未来发展的必然趋势以及当前的研究热点。目前，三维SAR成像技术主要包括曲线SAR、圆周SAR、层析SAR和线阵SAR。线阵SAR(Linear array SAR，LASAR)是利用与方位向垂直放置线阵天线，再结载荷平台的运动合成二维虚拟面阵的三维SAR成像技术。与其它三维SAR成像技术相比，线阵SAR具有多模式工作能力，除了可工作于传统侧视模式，还可工作于下视模式以及前视模式，在地形观测成像应用中更加灵活。
按照回波信号处理域的不同，线阵SAR成像算法可分为时域成像算法和频域成像算法。典型的频域成像算法包括三维距离多普勒(RD)算法、三维Chirp Scaling(CS)算法等。然而由于风流、湍流等因素将导致线阵SAR载荷平台往往偏离设计的理想轨迹，这些偏离误差极大程度上影响频域类成像算法的成像质量，甚至不能聚焦成像。经典的时域成像算法是三维后向投影(BP)算法(详见参考文献“G.Fornaro,F.Serafino,andF.Soldovieri.Three-dimensional focusing with multipass SAR data.IEEETrans.Geosci.Remote Sens,vol.41,no.3,pp.507–517,Mar.2003.”和“Shi Jun,ZhangXiaoling,Yang Jianyu,Wang yinbo.Surface-Tracing-Based LASAR 3-D ImagingMethod via Multiresolution Approximation.IEEE Trans.Geosci.Remote Sens,vol.46,no.11,pp.3719–3730,Nov.2008.”)，三维BP算法是一种精确的线阵SAR时域成像算法，它首先将合成孔径雷达原始数据沿距离向进行距离压缩(脉冲压缩)，然后通过选择不同慢时间观测空间中任意像素点在距离压缩后SAR数据空间中的数据，根据天线相位中心位置补偿方位向多普勒相位，然后再进行原始回波相干积累，最终获得各像素点散射系数的成像算法。相对频域成像算法，三维BP算法具有任意平台运动轨迹和任意成像空间成像的优势。但是，三维BP算法高精度聚焦成像的前提是精确已知阵列天线相位中心的位置，并且可以利用天线相位中心的位置进行运动误差补偿，详见“师君.双基地SAR与线阵SAR原理及成像技术研究[D].电子科技大学博士论文.2009”。当线阵SAR平台远动轨迹的位置测量精度较低或者存在未知运动误差时，会导致回波数据中存在相位误差，从而造成三维BP算法成像精度大大降低。因此，针对未知线阵SAR天线相位中心精确位置的情况下，如何实现相位误差估计以及三维BP算法成像具有重要意义。
SAR自聚焦成像算法是一种基于SAR回波数据进行运动误差估计与补偿的成像方法，自聚焦在SAR成像处理中占据着非常重要的地位。现有的SAR自聚焦成像算法主要分为两类。第一类自聚焦成像算法是基于SAR图像特征显点的方法，其代表算法为相位梯度自聚焦(详见参考文献“D.E.Wahal,P.H.Eichel,D.C.Ghiglia,and C.V.Jakowatz.Phasegradient autofocus-A robust tool for high resolution SAR phasecorrection.IEEE Transactions on Aerospace and Electionic Systems,1994,30(3):827-835.”)，该算法利用多个强点目标的聚焦深度来进行相位误差估计。第二类自聚焦成像算法是基于SAR图像质量的方法，该类算法首先将SAR图像的整体信息作为评价准则，目前常用的评价准则包括：熵最小化(见参考文献“M.Liu,C.and X.Shi.A back-projectionfast autofocus algorithm based on minimum entropy for SAR imaging.3rdInternational Asia-Pacific Conference,2005:310-314”)、对比度最大化(详见参考文献“J.Kolman.PACE:An autofocus algorithm for SAR.IEEE International RadarConference,2005:310-314”)、图像锐度最大化(详见参考文献“N.A.Joshua.An autofocusmethod for back projection imagery in synthetic aperture radar.IEEEGeoscience and Remote Sensing Letters,2012,9(1):104-108”)，该类算法通过调节相位误差的估计值使得SAR图像的质量指标达到最优，以此来确定补偿的相位误差。相对于基于SAR图像特征显点的自聚焦方法，基于SAR图像质量的自聚焦方法约束条件更少，更是适合于各种目标场景的SAR自聚焦成像。BP自聚焦成像算法是一类基于SAR图像质量的自聚焦算法，其主要思想是根据SAR图像最优质量指标来估计方位向相位误差，当SAR图像质量指标达到最优时，即认为SAR图像的聚焦最好。
为了提高线阵SAR后向投影算法成像精度，本发明提出一种基于半正定规划的线阵SAR后向投影自聚焦成像方法。该方法通过构造后向投影成像积累的基函数矩阵，利用后向投影成像积累基函数矩阵与相位误差矩阵的乘积矩阵的最小秩作为最优化目标，通过凸优化理论中的半正定规划方法求解约束条件下的最优矩阵实现线阵SAR相位误差，从而实现线阵SAR的自聚焦成像。根据本人了解，当前还没有出现基于半正定规划的线阵SAR后向投影自聚焦成像方法。
发明内容：
为了提高线阵SAR的成像精度，本发明将凸优化理论中的半正定规划方法应用于线阵SAR自聚焦成像，提供了基于半正定规划的线阵SAR后向投影自聚焦成像方法。该方法的主要思路是：根据线阵SAR回波信号，构造线阵SAR后向投影成像积累的基函数矩阵，利用后向投影积累基函数矩阵与相位误差的乘积矩阵的最小秩作为最优化目标，通过半正定规划方法估计求解约束条件下的最优相位误差矩阵实现线阵SAR相位误差的估计和自聚焦成像。
为了方便描述本发明的内容，首先作以下术语定义：
定义1、范数
设X是数域上线性空间，表示复数域，若它满足如下性质：||X||≥0，且||X||＝0仅有X＝0；||aX||＝|a|||X||，a为任意常数；||X1+X2||≤||X1||+||X2||，则称||X||为X空间上的范数(norm)，其中X1和X2为X空间上的任意两个值。对于定义1中的N×1维离散信号向量X＝[x1,x2,…,xN]T，向量X的LP范数表达式为其中xi为向量X的第i个元素，..表示绝对值求和符号，向量X的L1范数表达式为向量X的L2范数表达式为向量X的L0范数表达式为且xi≠0。详见文献“矩阵理论”，黄廷祝等编著，高等教育出版社出版。
定义2、凸优化
凸优化是指目标函数为凸函数，且由约束条件得到的定义域为凸集的优化问题，典型的凸优化问题包括：线性规划、二次规划、半正定规划等。详见文献“S.Boyd,L.Vandenberghe,Convex optimization,Cambridge university press,2009.”。
定义3、标准半正定规划
标准半正定规划是一种凸优化问题。给定半正定矩阵和系数矩阵i＝1,…,P，测量向量i＝1,…,P，其中i为自然数，p表示为测量向量的总个数，标准半正定规划是为了求解以下最优化问题的最优解
s.t tr(AiX)＝bi,i＝1,…,P
X≥0
其中tr(·)为矩阵求秩运算符，为求解满足目标函数的值最小时对应的矩阵X，s.t表示为投影或受约束符号，表示为实数集，表示为n×n维矩阵实数集，表示为n×1维向量实数集，∈表示为属于符号。半正定问题求解可利用各种凸优化计算软件进行求解，如CVX软件工具箱等，详见文献“G.Michael,S.Boyd,and Yinyu Ye.cvxusers’guide.Technical Report Build 711,Citeseer,2009.”。
定义4、线阵合成孔径雷达(Linear array SAR,LASAR,线阵SAR)
线阵合成孔径雷达成像是将线性阵列天线固定于载荷运动平台上并与平台运动方向垂直，结合运动平台的运动以合成二维平面阵列以实现阵列平面二维成像，再利用雷达波束向回波延时实现距离一维成像，从而实现观测目标三维成像的一种合成孔径雷达技术。
定义5、合成孔径雷达回波数据距离压缩
标准合成孔径雷达距离压缩方法是指利用合成孔径雷达发射信号参数，采用匹配滤波技术对合成孔径雷达的距离向信号进行信号聚焦成像的过程。详见文献“雷达成像技术”，保铮，邢孟道，王彤，电子工业出版社，2005。
定义6、合成孔径雷达原始回波仿真方法
合成孔径雷达原始回波仿真方法是指基于合成孔径雷达成像原理仿真出一定系统参数条件下具有合成孔径雷达回波信号特性的原始信号的方法，详细内容可参考文献：“张朋，合成孔径雷达回波信号仿真研究，西北工业大学博士论文，2004”。
定义7、线阵SAR观测场景目标空间
线阵SAR观测场景目标空间是指现实空间中所有待观测场景目标散射点的位置集合。观测场景目标空间在不同空间坐标系下有不同的表示，但一旦坐标系确立以后其表示是唯一的。一般情况下为了方便成像，线阵SAR观测场景目标空间取为地面坐标系。本发明中用以下数学关系表示观测场景目标空间Ω：
其中ex、ey和ez表示构成观测场景目标空间Ω的地表正交坐标基，分别表示水平横向、水平纵向和垂直地表的高度向，Ps为场景目标空间中一个分辨单元位置向量，x、y和z分别表示该分布单元的水平横向、水平纵向和高度向坐标，表示实数集。
定义8、线阵SAR成像空间
线阵SAR成像空间是指将场景目标空间中的散射点投影到切航迹向－沿航迹向－距离向的三维空间坐标系，该空间由线阵SAR成像空间中的三个相互正交的坐标基确定。本发明中用以下数学关系表示成像空间Ξ：
其中和表示构成线阵SAR成像空间Ξ的正交坐标基，分别表示切航迹向、沿航迹向和距离向，Pω为成像空间中的待观测点位置向量，u、v和z分别表示该点的切航迹向、沿航迹向和距离向坐标，表示实数域。
定义9、合成孔径雷达后向投影成像算法
合成孔径雷达的后向投影成像算法，简称BP成像算法。BP成像算法首先利用雷达平台的轨迹信息求出雷达平台与场景像素点的距离历史，然后通过距离历史找出回波数据中对应的复数据，再进行相位补偿并相干累加，从而得到该像素点的复图像值。详见“师君.双基地SAR与线阵SAR原理及成像技术研究[D].电子科技大学博士论文.2009”。
定义10、矩阵的秩
在线性代数中，矩阵的秩是指该矩阵中线性无关的纵列及横行的极大数目。在本发明中，求解矩阵的秩用符号tr(·)表示。
定义11、合成孔径雷达的距离向快时刻和方位向慢时刻
合成孔径雷达的距离向快时刻是指在一个雷达系统工作的脉冲重复周期内，距离向回波信号采样过程中不同采样点的时间间隔变量。合成孔径雷达系统以一定时间长度的重复周期发射和接收脉冲信号，方位向慢时刻表示为一个以脉冲重复周期为步长的离散化时间变量，其中每一个脉冲重复周期离散时间变量值为一个方位向慢时刻。详见文献“合成孔径雷达成像原理”，皮亦鸣等编著，电子科技大学出版社出版。
本发明提供的基于半正定规划的线阵SAR后向投影自聚焦成像方法，它包括以下步骤：
步骤1、初始化线阵SAR系统参数和回波数据：
初始化线阵SAR系统参数包括：运动平台的速度矢量，记做V＝[vx,vy,vz]，其中vx为平台在水平面横轴向上的速度，vy为平台在水平面纵轴向上的速度，vz为平台在地面垂直高度向上的速度；线阵天线各阵元的初始位置矢量，记做Pn(0)＝[xn(0),yn(0),zn(0)]，其中n为天线阵元序号，为自然数，n＝1,2,...,N，N为线阵天线的阵元总数，xn(0)为第n个天线阵元在水平面横轴向上的初始位置，yn(0)为第n个天线阵元在水平面纵轴向上的初始位置，zn(0)为第n个天线阵元在地面垂直高度向上的初始位置；线阵天线的总长度，记做L；雷达工作的中心频率，记做fc；雷达工作的载频波长，记做λ；雷达发射基带信号的信号带宽，记做Br；雷达发射信号的脉冲宽度，记做TP；雷达信号的调频斜率，记做fdr；雷达接收系统的采样频率，记做fs；雷达系统的脉冲重复频率，记做PRF；雷达系统的脉冲重复时间，记为PRI；天线在方位向的有效孔径长度，记做Da；光在空气中的转播速度，记做c；距离向快时刻的总数，记为T，距离向快时刻序列记为t＝1,2,…,T，其中t为第t个距离向快时刻；方位向慢时刻的总数，记为K，方位向慢时刻序列记为l＝1,2,…,K，其中l为第l个方位向慢时刻；上述参数均为线阵SAR系统的标准参数，其中线阵天线的阵元总数N，线阵天线总长度L，雷达中心频率fc，雷达载频波长λ，雷达发射基带信号的信号带宽Br，雷达发射信号脉冲宽度TP，雷达发射信号调频斜率fdr，雷达接收系统的采样频率fs,雷达系统的脉冲重复频率PRF，天线在方位向的有效孔径长度Da在线阵SAR系统设计过程中已经确定；平台速度矢量V及线阵天线各阵元初始位置矢量Pn(0)在线阵SAR观测方案设计中已经确定；根据线阵SAR成像系统方案和观测方案，线阵SAR成像方法需要的初始化成像系统参数均为已知；
线阵SAR第n阵元在第l个方位向慢时刻和第t个距离向快时刻采样得到的原始回波数据，记为s(t,l,n)，t＝1,2,…,T，l＝1,2,…,K，n＝1,2,...,N；在实际线阵SAR系统中，原始回波数据s(t,l,n)可由线阵SAR系统数据接收机提供；在仿真线阵SAR成像过程中，原始回波数据s(t,l,n)可根据线阵SAR成像系统参数，采用标准合成孔径雷达原始回波仿真方法产生得到；在线阵SAR数据进行成像之前，原始回波数据s(t,l,n)均已知；
步骤2、初始化线阵SAR的观测场景目标空间参数：
初始化线阵SAR的观测场景目标空间参数，包括：以雷达波束照射场区域地平面和垂直于该地平面向上的单位向量所构成的三维空间直角坐标作为线阵SAR的观测场景目标空间，记为Ω；将观测场景目标空间Ω均匀划分成大小相等的立体单元格，也称为分辨单元，单元网格在水平横向、水平纵向和高度向边长分别记为dx、dy和dz，分辨单元在水平横向、水平纵向和高度向的总数分别记为Mx、My和Mz；观测场景目标空间Ω中第m个单元格的坐标矢量，记做qm，m表示观测场景目标空间Ω中第m个单元格，m为正整数，m＝1,2,…,M,其中M记为观测场景目标空间Ω中的单元格总数，且M＝Mx·My·Mz；
步骤3、线阵SAR原始回波信号距离压缩：
采用标准合成孔径雷达距离压缩方法对线阵SAR原始回波数据s(t,l,n)，t＝1,2,…,T，l＝1,2,…,K，n＝1,2,...,N进行距离压缩，得到距离压缩后的线阵SAR回波数据，记为sr(t,l,n)，t＝1,2,…,T，l＝1,2,…,K，n＝1,2,...,N，其中s(t,l,n)为步骤1中得到的线阵SAR原始回波数据，T为步骤1中得到的距离向快时刻总数，K为步骤1中得到的方位向慢时刻总数，N为步骤1中得到的线阵天线的阵元总数；
步骤4、计算观测场景单元对应的回波信号序号：
采用公式Pn(l)＝Pn(0)+V·l/PRF，n＝1,2,…,N，l＝1,2,…,K，计算得到第n个线阵天线阵元在第l个方位向慢时刻的位置矢量，记为Pn(l)，其中V为步骤1中初始化得到的平台速度矢量，Pn(0)为步骤1中初始化得到的线阵天线各阵元初始位置矢量，PRF为步骤1中初始化得到的雷达系统的脉冲重复频率，N为步骤1中得到的线阵天线阵元总数，K为步骤1中得到的方位向慢时刻总数；
采用公式R(Pn(l),qm)＝||Pn(l)-qm||2，n＝1,2,…,N，l＝1,2,…,K，m＝1,2,…,M，计算得到在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格到第n个线阵天线阵元的距离，记为R(Pn(l),qm)，其中||·||2表示向量的L2范数，qm为步骤2中初始化得到观测场景目标空间Ω中第m个单元格的坐标矢量，M为步骤2中初始化的场景目标空间Ω中单元格总数；
采用公式I(l,n,m)＝ceil(2·R(Pn(l),qm)·fs/c)，n＝1,2,…,N，l＝1,2,…,K，m＝1,2,…,M，计算得到在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格在第n个线阵天线阵元对应的原始回波序号，记为I(l,n,m)，其中c为步骤1中初始化得到的光在空气中的传播速度，fs为步骤1中初始化得到的雷达接收系统的采样频率，ceil(·)为上取整数运算符号；
采用公式Ir(l,n,m)＝2·R(Pn(l),qm)·fs/c-I(l,n,m)，n＝1,2,…,N，l＝1,2,…,K，m＝1,2,…,M，计算得到在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格在第n个线阵天线阵元对应的原始回波序号残余量，记为Ir(l,n,m)；
步骤5、利用sinc函数进行回波数据插值：
设定sinc插值函数的窗口长度，记为W0，其中W0为大于2的偶数；令w为整数，w的值为w＝-W0/2,-W0/2+1,…,W0/2-1,W0/2；
采用公式计算得到sinc函数插值后在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格在第n个线阵天线阵元对应的回波数据，记为sc(l,n,m)，其中sr(I(l,n,m)+w,l,n)为回波数据sr(t,l,n)中距离向快时刻t＝I(l,n,m)+w时对应的值，sr(t,l,n)为步骤3得到的距离压缩后线阵SAR回波数据，I(l,n,m)为步骤4得到的在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格在第n个线阵天线阵元对应的原始回波序号，Ir(l,n,m)为步骤4得到的在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格在第n个线阵天线阵元对应的原始回波序号残余量，sinc(·)为sinc函数运算符号，为元素w值从-W0/2至W0/2范围内的函数求和符号；
步骤6、构造线阵SAR后向投影积累基函数矩阵：
距离压缩后线阵SAR回波数据的后向投影积累基函数矩阵，记为B，矩阵B的维数大小M×(K·N)，其中M为步骤2中得到的观测场景目标空间Ω中的单元格总数，K为步骤1中得到的方位向慢时刻总数，N为步骤1中得到的线阵天线的阵元总数；矩阵B中第m行第g列的元素值，记为B(m,g)，其中m＝1,2,…,M，g为自然数，g的取值为g＝1,2,…,KN；
采用公式B(m,(l-1)·N+n)＝sc(l,n,m)·exp(j·4π·R(Pn(l),qm)/λ)，n＝1,2,…,N，l＝1,2,…,K，m＝1,2,…,M，计算得到矩阵B第m行第(l-1)·N+n列的元素值，记为B(m,(l-1)·N+n)，其中sc(l,n,m)为步骤5得到的sinc函数插值后在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格在第n个线阵天线阵元对应的回波数据，R(Pn(l),qm)为步骤4得到的第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格到第n个线阵天线阵元的距离，λ为步骤1初始化得到的雷达工作的载频波长，exp(·)为自然常数e为底的指数运算符号，j为虚数符号，π为圆周率；
步骤7、选择观测场景目标空间中的强散射单元：
设定选择线阵SAR观测场景目标空间Ω中强散射单元的幅度阈值，记为ε0；
采用公式f＝B·IK·N计算得到线阵SAR观测场景散射系数的后向投影成像向量，记为f，向量f的维数大小为M×1，向量f的第m个元素记为fm，m＝1,2,…,M，其中B为步骤6得到的后向投影积累基函数矩阵，向量IK·N记为(K·N)×1的列向量，且向量IK·N中所有元素都为1；
采用公式β＝{m||fm|≥ε0·max(f),m＝1,2,…,M}计算得到观测场景目标空间Ω中强散射单元所在的序号向量，记为β，向量β的第y个元素记为βy，其中y为自然数，y的取值为y＝1,2,…,Q，Q为向量β中元素的总个数，fm为向量f的第m个元素，|·|表示为取绝对值运算符号，m||f(m)|≥ε0·max(f)表示满足条件|f(m)|≥ε0·max(f)所对应的序号m，max(·)表示为取最大值运算符号；
在矩阵B中选择向量β中元素值所对应的行向量按照序号大小组成矩阵，得到观测场景目标空间Ω中强散射单元的后向投影积累基函数矩阵，记为BS，矩阵BS的维数大小为Q×(K·N)，其中B为步骤6得到的后向投影积累基函数矩阵；
步骤8、初始化相位误差迭代估计的参数：
初始化相位误差迭代估计的参数包括：迭代估计过程的最大迭代次数，记做MaxIter；k记为迭代估计过程的第k次迭代，k为自然数，k初始值设置为k＝0，并且k的取值范围为k＝0,1,2,…,MaxIter；第k次迭代估计得到的相位误差向量，记为φ(k)，k＝0,1,2,…,MaxIter，其中，φ(k)为列向量，φ(k)的维数大小为(K·N)×1；初始化第k次迭代估计得到的相位误差向量的值为φ(k)＝0，k＝0,1,2,…,MaxIter；第k次迭代估计得到的线阵SAR观测场景后向投影成像结果向量，记为F(k)，k＝0,1,2,…,MaxIter，F(k)为列向量，向量F(k)的维数大小为M×1；初始化第k次迭代估计得到的线阵SAR观测场景后向投影成像结果向量的值为F(k)＝0，k＝0,1,2,…,MaxIter；
步骤9、利用半正定规划求解相位误差矩阵：
在第k次迭代中，若k＝0时采用公式γ(0)＝exp(jφ(0))计算得到相位误差迭代估计过程中第0次迭代的相位误差指数向量，记为γ(0)，其中φ(0)为步骤8初始化得到的第0次迭代估计相位误差向量，exp(·)为自然常数e为底的指数运算符号；若k≥1时采用公式γ(k)＝exp(jφ(k-1))计算得到相位误差迭代估计过程中第k次迭代的相位误差指数向量，记为γ(k)，其中φ(k-1)为第k-1次迭代估计得到的相位误差向量；
采用公式计算得到第k次迭代中强散射单元对应的相位误差半正定矩阵，记为C(k)，k＝0,1,2,…,MaxIter，其中上标H表示为共轭转置运算符号，BS为步骤7得到的观测场景目标空间Ω中强散射单元的后向投影积累基函数矩阵，为矩阵BS的共轭转置矩阵，diag(·)表示将向量元素值赋予对角矩阵中的对角元素运算符号
采用标准半正定规划计算以下公式
s.t X(g,g)＝1,g＝1,…,K·N
X≥0
得到第k次迭代中的相位误差矩阵，记为X(k)，矩阵X(k)的维数大小为(K·N)×(K·N)，其中表示求取满足括号中最小值时对应的自变量X的最优值，tr(·)表示为矩阵求秩运算符，s.t表示为受约束符号，X(g,g)表示为矩阵X中第g行第g列的元素值，K为步骤1初始化得到的方位向慢时刻总数，N为步骤1初始化得到的线阵天线的阵元总数；
步骤10、估计相位误差向量：
对于步骤9得到的相位误差矩阵X(k)，如果矩阵X(k)满足tr(X(k))＝1，采用公式X(k)＝α·αH对矩阵X(k)进行向量分解，分解得到的列向量记为α，列向量α的维数大小为(K·N)×1，然后采用公式φ(k)＝∠α得到第k次迭代过程的相位误差向量，其中∠表示为求解角度运算符号，其中tr(·)表示为矩阵求秩运算符；
如果步骤9得到的矩阵X(k)满足tr(X(k))≠1，采用公式X(k)＝D·DH对矩阵X(k)进行平方矩阵分解，分解得到的矩阵记为D，矩阵D的维数大小为(K·N)×MS，矩阵D的第r个列向量记为dr，r为自然数，r的取值为r＝1,2,…,Ms；令Mr为正整数，Mr的值要大于Ms；随机产生Mr个均值为0、方差为IK·N的独立同分布复高斯向量，分别记为其中u1为第1个复高斯向量，u2为第2个复高斯向量，为第Mr个复高斯向量，IK·N为步骤7得到的元素都为1的向量；令w为自然数，w的取值为采用公式ηw＝exp(j∠D·uw)，w＝1,2,…,Mr，计算得到相位向量，记为ηw；采用公式计算得到最优的相位向量，记为ηopt，其中表示为在序号为1≤w≤Mr范围中求取满足中括号目标函数最小值时对应的向量序号，C(k)为步骤9得到的第k次迭代中强散射单元对应的相位误差半正定矩阵；然后，采用公式φ(k)＝∠ηopt得到第k次迭代过程的相位误差向量；
步骤11、线阵SAR观测场景目标空间后向投影成像
采用公式计算得到第k次迭代过程的相位误差指数补偿向量，记为向量为列向量，的维数大小为(K·N)×1，其中φ(k)为步骤10得到的第k次迭代过程的相位误差向量，exp(·)为自然常数e为底的指数运算符号；
采用公式计算得到第k次迭代过程中线阵SAR观测场景后向投影成像结果，其中B为步骤6得到的线阵SAR后向投影积累基函数矩阵；
步骤12、迭代终止判定：
如果且k≤MaxIter，则k←k+1，执行步骤9至步骤12，否则终止迭代，此刻第k次迭代得到的F(k)即为线阵SAR观测场景最终的后向投影成像结果，其中k表示迭代估计过程中的第k迭代次数，MaxIter为步骤8中初始化得到的算法重构处理的最大迭代次数，F(k)为步骤11得到的第k次迭代过程中观测场景的后向投影成像结果，F(k-1)为第k-1次迭代过程中观测场景的后向投影成像结果，||·||2表示向量的L2范数运算符号，表示向量的L2范数平方运算符号，←表示赋值运算符号；最后将观测场景目标空间散射系数向量F(k)转换成三维矩阵形式，得到线阵SAR观测场景的最终三维成像结果。
本发明的创新点在于将线阵SAR后向投影算法成像过程转换为矩阵与向量相乘形式，采用凸优化理论中的半正定规划方法对线阵SAR回波数据中相位误差进行估计，提供了基于半正定规划的线阵SAR后向投影自聚焦成像方法，该方法通过构造线阵SAR后向投影积累基函数矩阵，利用后向投影积累基函数矩阵与相位误差矩阵的乘积矩阵的最小秩，通过求解约束条件下的最优矩阵实现线阵SAR相位误差的估计和自聚焦成像。
本发明的优点在于无需已知相位误差分布的先验信息，可估计线阵SAR回波数据中存在的低频、高频以及随机分布的相位误差，适用于复杂运动轨迹情况下的线阵SAR自聚焦成像。本发明可以应用于合成孔径雷达成像和地球遥感等领域。
附图说明：
图1为本发明所提供的基于半正定规划的线阵SAR后向投影自聚焦成像方法的处理流程示意框图。
图2为本发明具体实施方式采用的线阵SAR系统仿真参数表。
具体实施方式
本发明主要采用仿真实验的方法进行验证，所有步骤、结论都在MATLABR2008b上验证正确。具体实施步骤如下：
步骤1、初始化线阵SAR系统参数和回波数据：
初始化线阵SAR系统参数包括：平台速度矢量V＝[0,100,0]m/s，平台在水平面横轴向上的速度vx＝0m/s，平台在水平面纵轴向上的速度vy＝100m/s，平台在地面垂直高度向上的速度vz＝0m/s；线阵天线各阵元初始位置矢量Pn(0)＝[0.015·(n-251),0,3000]m，其中n为天线阵元序号，为自然数，n＝1,2,...,N，N为线阵天线的阵元总数且N＝501，第n个天线阵元在水平面横轴向上的初始位置xn(0)＝0.015·(n-251)m，第n个天线阵元在水平面纵轴向上的初始位置yn(0)＝0m，第n个天线阵元在地面垂直高度向上的初始位置zn(0)＝3000m；线阵天线的总长度L＝7.5m；雷达工作的中心频率fc＝10GHz；雷达工作的载频波长λ＝0.03m；雷达发射基带信号的信号带宽Br＝300MHz；雷达发射信号的脉冲宽度TP＝10-6s；雷达发射信号的调频斜率fdr＝3×1014Hz/s；雷达接收系统的采样频率fs＝500MHz；雷达系统的脉冲重复频率PRF＝500Hz；雷达系统的脉冲重复时间PRI＝2×10-3s；天线在方位向的有效孔径长度Da＝1m；光在空气中的转播速度c＝3×108m/；s距离向快时刻总数T＝256，距离向快时刻序列t＝1,2,…,T，其中t为第t个距离向快时刻；方位向慢时刻总数K＝256，方位向慢时刻序列l＝1,2,…,K，其中l为第l个方位向慢时刻；
线阵SAR第n阵元在第l个方位向慢时刻和第t个距离向快时刻采样得到的原始回波数据，记为s(t,l,n)，t＝1,2,…,T，l＝1,2,…,K，n＝1,2,...,N，其中T＝256，K＝256，N＝501；在仿真线阵SAR成像过程中，仿真点目标散射体的个数为6个，它们的散射系数值均为1，坐标位置分别为[5,5,12]，[5,-5,12]，[-5,-5,12]，[5,5,-12]，[5,-5,-12]，[-5,-5,-12]，单位均为m，原始回波数据s(t,l,n)采用标准合成孔径雷达原始回波仿真方法产生得到；
步骤2、初始化线阵SAR的观测场景目标空间参数：
初始化线阵SAR的观测场景目标空间参数，包括：以雷达波束照射场区域地平面和垂直于该地平面向上的单位向量所构成的三维空间直角坐标作为线阵SAR的观测场景目标空间Ω；将观测场景目标空间Ω均匀划分成大小相等立体单元格，也称为分辨单元，分辨单元在水平横向、水平纵向和高度向的边长分别为dx＝0.5m、dy＝0.5m和dz＝0.5m，分辨单元在水平横向、水平纵向和高度向的总数分别为Mx＝128、My＝128和Mz＝128；观测场景目标空间Ω中第m个单元格的坐标矢量qm＝[(x′-64)·dx,(y′-64)·dy,(z′-64)·dz]，其中x′＝1,2,…,Mx，y′＝1,2,…,My，z′＝1,2,…,Mz，m＝((x′-1)·My+y′-1)Mz+z′，m为观测场景目标空间Ω中第m个单元格，m＝1,2,…,M,M为观测场景目标空间Ω中的单元格总数且M＝Mx·My·Mz＝2097152；
步骤3、线阵SAR原始回波信号距离压缩：
采用标准的合成孔径雷达距离压缩方法对线阵SAR原始回波数据s(t,l,n)，t＝1,2,…,T，l＝1,2,…,K，n＝1,2,...,N进行距离压缩，得到距离压缩后的线阵SAR回波数据，记为sr(t,l,n)，t＝1,2,…,T，l＝1,2,…,K，n＝1,2,...,N，其中s(t,l,n)为步骤1得到的线阵SAR原始回波数据，T为步骤1初始化得到的距离向快时刻总数T＝256，K为步骤1初始化得到的方位向慢时刻总数K＝256，N为步骤1初始化得到的线阵天线的阵元总数N＝501；
步骤4、计算观测场景单元对应的回波信号序号：
采用公式Pn(l)＝Pn(0)+V·lPRF，n＝1,2,…,N，l＝1,2,…,K，N＝501，K＝256，计算得到第n个线阵天线阵元在第l个方位向慢时刻的位置矢量，记为Pn(l)，其中V为步骤1中初始化得到的平台速度矢量V＝[0,100,0]m/s，Pn(0)为步骤1中初始化得到的线阵天线各阵元初始位置矢量Pn(0)＝[0.015·(n-251),0,3000]m，PRF步骤1中初始化得到的雷达系统的脉冲重复频率PRF＝500Hz，K为步骤1初始化得到的方位向慢时刻总数K＝256，N为步骤1初始化得到的线阵天线的阵元总数N＝501；
采用公式R(Pn(l),qm)＝||Pn(l)-qm||2，n＝1,2,…,N，l＝1,2,…,K，m＝1,2,…,M，N＝501，K＝256，M＝2097152，计算得到在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格到第n个线阵天线阵元的距离，记为R(Pn(l),qm)，其中||·||2向量的L2范数，qm为步骤2中初始化得到观测场景目标空间Ω中第m个单元格的坐标矢量，M为步骤2中初始化的场景目标空间Ω中单元格总数M＝2097152；
采用公式I(l,n,m)＝ceil(2·R(Pn(l),qm)·fs/c)，n＝1,2,…,N，l＝1,2,…,K，m＝1,2,…,M，N＝501，K＝256，M＝2097152，计算得到在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格在第n个线阵天线阵元对应的原始回波序号，记为I(l,n,m)，其中c为步骤1中初始化得到的光在空气中的传播速度c＝3×108m/s，fs为步骤1中初始化得到的雷达接收系统的采样频率fs＝500MHz，ceil(·)为上取整数运算符号；
采用公式Ir(l,n,m)＝2·R(Pn(l),qm)·fs/c-I(l,n,m)，n＝1,2,…,N，l＝1,2,…,K，m＝1,2,…,M，N＝501，K＝256，M＝2097152，计算得到在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格在第n个线阵天线阵元对应的原始回波序号残余量，记为Ir(l,n,m)；
步骤5、利用sinc函数进行回波数据插值：
设定sinc插值函数的窗口长度为W0＝8；令w为整数，w＝-W0/2,-W0/2+1,…,W0/2-1,W0/2，其中W0＝8；
采用公式n＝1,2,…,N，l＝1,2,…,K，m＝1,2,…,M，N＝501，K＝256，M＝2097152，计算得到sinc函数插值后在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格在第n个线阵天线阵元对应的回波信号，记为sc(l,n,m)，其中sr(I(l,n,m)+w,l,n)为回波数据sr(t,l,n)中距离向快时刻t＝I(l,n,m)+w时对应的值，sr(t,l,n)为步骤3得到的距离压缩后线阵SAR回波数据，I(l,n,m)为步骤4得到的在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格在第n个线阵天线阵元对应的原始回波序号，Ir(l,n,m)为步骤4得到的在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格在第n个线阵天线阵元对应的原始回波序号残余量，sinc(·)为sinc函数运算符号，为元素w值从-W0/2至W0/2范围内的函数求和符号，W0为sinc插值函数的窗口长度窗口长度W0＝8；
步骤6、构造线阵SAR后向投影积累基函数矩阵：
距离压缩后线阵SAR回波数据的后向投影积累基函数矩阵，记为B，矩阵B的维数大小M×(K·N)，其中M为步骤2中得到的观测场景目标空间Ω中单元格总数M＝2097152，K为步骤1中得到的方位向慢时刻总数K＝256，N为线阵天线的阵元总数N＝501；矩阵B中第m行第g列的元素值，记为B(m,g)，其中m＝1,2,…,M，M＝2097152，g为自然数，g的取值为g＝1,2,…,K·N，N＝501，K＝256；
采用公式B(m,(l-1)·N+n)＝sc(l,n,m)·exp(j·4π·R(Pn(l),qm)/λ)，n＝1,2,…,N，l＝1,2,…,K，m＝1,2,…,M，N＝501，K＝256，M＝2097152，计算得到矩阵B第m行第(l-1)·N+n列的元素值，其中sc(l,n,m)为步骤5得到的sinc函数插值后在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格在第n个线阵天线阵元对应的原始回波，R(Pn(l),qm)为步骤4得到的第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格到第n个线阵天线阵元的距离，λ为步骤1初始化得到的雷达工作载频波长λ＝0.03m，exp(·)为自然常数e为底的指数运算符号，j为虚数符号，圆周率π＝3.1415926；
步骤7、选择观测场景目标空间中的强散射单元：
设定选择线阵SAR观测场景目标空间Ω中强散射单元的幅度阈值ε0＝0.2；
采用公式f＝B·IK·N计算得到线阵SAR观测场景目标空间Ω中散射系数的后向投影成像向量，记为f，向量f的维数大小为M×1，其中M＝2097152，向量f的第m个元素记为fm，m＝1,2,…,M，其中B为步骤6得到的距离压缩后线阵SAR回波数据的后向投影积累基函数矩阵，向量IK·N记为(K·N)×1的列向量，且向量IK·N的所有元素都为1，其中N＝501，K＝256；
采用公式β＝{m||fm|≥ε0·max(|f|),m＝1,2,…,M}计算得到观测场景目标空间Ω中强散射单元所在的序号向量，记为β，向量β的第y个元素记为βy，其中y为自然数，y＝1,2,…,Q，Q为向量β中元素的总个数，fm为向量f的第m个元素，|·|表示为取绝对值运算符号，m||f(m)|≥ε0·max(|f|)表示满足条件|f(m)|≥ε0·max(|f|)所对应的序号m，max(·)表示为取最大值运算符号；
在矩阵B中选择向量β中元素值所对应的行向量按照序号大小组成矩阵，得到观测场景目标空间Ω中强散射单元对应的后向投影积累基函数矩阵，记为BS，其中B为步骤6得到的距离压缩后线阵SAR回波数据的后向投影积累基函数矩阵，β为步骤7得到的观测场景目标空间Ω中强散射单元所在的序号向量；
步骤8、初始化相位误差迭代估计的参数：
初始化相位误差迭代估计的参数包括：迭代估计过程的最大迭代次数MaxIter＝32；k记为迭代估计过程的第k次迭代，k为自然数，k初始值设置为k＝0，并且k的取值为k＝0,1,2,…,MaxIt，er其中MaxIter＝32；第k次迭代估计得到的相位误差向量，记为φ(k)，k＝0,1,2,…,MaxIter且MaxIter＝32，其中φ(k)为列向量，φ(k)维数大小为(K·N)×1，其中N＝501，K＝256；初始化第k次迭代估计得到的相位误差向量的值为φ(k)＝0，k＝0,1,2,…,MaxIter且MaxIter＝32；第k次迭代估计得到的线阵SAR观测场景后向投影成像结果向量，记为F(k)，k＝0,1,2,…,MaxIter且MaxIter＝32，F(k)为列向量，向量F(k)的维数大小为M×1，其中M＝2097152；初始化第k次迭代估计得到的线阵SAR观测场景后向投影成像结果向量的值为F(k)＝0，k＝0,1,2,…,MaxIter且MaxIter＝32；
步骤9、利用半正定规划求解相位误差矩阵：
在第k次迭代中，若k＝0时采用公式γ(0)＝exp(jφ(0))计算得到相位误差迭代估计过程中第0次迭代的相位误差指数向量，记为γ(0)，其中φ(0)为步骤8初始化得到的第0次迭代估计相位误差向量，exp(·)为自然常数e为底的指数运算符号；若k≥1时采用公式γ(k)＝exp(jφ(k-1))计算得到相位误差迭代估计过程中第k次迭代的相位误差指数向量，记为γ(k)，其中φ(k-1)为第k-1次迭代估计得到的相位误差向量；
采用公式计算得到第k次迭代中强散射单元对应的相位误差半正定矩阵，记为C(k)，k＝0,1,2,…,MaxIter，其中上标H表示为共轭转置运算符号，BS为步骤7得到的观测场景目标空间Ω中强散射单元的后向投影积累基函数矩阵，MaxIter为步骤8得到的迭代估计过程的最大迭代次数MaxIter＝32，diag(·)表示将向量元素值赋予对角矩阵中的对角元素运算符号；
采用标准半正定规划计算以下公式
s.t X(g,g)＝1,g＝1,…,K·N
X≥0
得到第k次迭代中的相位误差矩阵，记为X(k)，矩阵X(k)的维数大小为(K·N)×(K·N)，其中表示求取满足括号中最小值时对应的自变量X的最优值，tr(·)表示为矩阵求秩运算符，s.t表示为受约束符号，X(g,g)表示为矩阵X中第g行第g列的元素值，K为步骤1初始化得到的方位向慢时刻总数K＝256，N为步骤1初始化得到的线阵天线的阵元总数N＝501；
步骤10、估计相位误差向量：
对于步骤9得到的相位误差矩阵X(k)，如果矩阵X(k)满足tr(X(k))＝1，采用公式X(k)＝α·αH对矩阵X(k)进行向量分解，得到的列向量记为α，列向量α的维数大小为(K·N)×1，其中tr(·)表示为矩阵求秩运算符，K为步骤1初始化得到的方位向慢时刻总数K＝256，N为步骤1初始化得到的线阵天线的阵元总数N＝501；然后采用公式φ(k)＝∠α得到第k次迭代过程的相位误差向量，其中∠表示为求解角度运算符号；
如果步骤9得到的相位误矩阵X(k)满足tr(X(k))≠1，采用公式X(k)＝D·DH对矩阵X(k)进行平方矩阵分解，得到的矩阵记为D，矩阵D的维数大小为(K·N)×MS，Ms表示为矩阵D的列数，其中K＝256，N＝501，矩阵D的第r个列向量记为dr，r为自然数，r的取值为r＝1,2,…,Ms；令Mr为正整数，初始化Mr的值为Mr＝16384；随机产生Mr个均值为0、方差为IK·N的独立同分布复高斯向量，分别记为其中u1为第1个复高斯向量，u2为第2个复高斯向量，为第Mr个复高斯向量，其中IK·N为步骤7得到的元素都为1的向量；令w为自然数，w的取值为w＝1,2,…,Mr；采用公式ηw＝exp(j∠D·uw)，w＝1,2,…,Mr且Mr＝16384，计算得到的相位向量，记为ηw；采用公式计算得到最优的相位向量，记为ηopt，其中表示为在序号为1≤w≤Mr范围中求取满足括号中最小值时对应的向量，C(k)为步骤9得到的第k次迭代中强散射单元对应的相位误差半正定矩阵，上标H表示为共轭转置运算符号；然后，采用公式φ(k)＝∠ηopt得到第k次迭代过程的相位误差向量；
步骤11、线阵SAR观测场景目标空间后向投影成像：
采用公式计算得到第k次迭代过程的相位误差指数补偿向量，记为向量为列向量，的维数大小为(K·N)×1，其中φ(k)为步骤10得到的第k次迭代过程的相位误差向量，K为步骤1初始化得到的方位向慢时刻总数K＝256，N为步骤1初始化得到的线阵天线的阵元总数N＝501；
采用公式计算得到第k次迭代过程中线阵SAR观测场景目标空间的后向投影成像结果，其中B为步骤6得到的距离压缩后线阵SAR回波数据的后向投影积累基函数矩阵；
步骤12、迭代终止判定：
如果且k≤MaxIter，则k←k+1，执行步骤9至步骤12，否则终止迭代，此刻第k次迭代得到的F(k)即为线阵SAR观测场景最终的后向投影成像结果，其中k表示迭代估计过程中的第k迭代次数，MaxIter为步骤8得到的迭代估计过程的最大迭代次数MaxIter＝32，F(k)为步骤11得到的第k次迭代过程中线阵SAR观测场景目标空间的后向投影成像结果，F(k-1)为第k-1次迭代过程中线阵SAR观测场景目标空间的后向投影成像结果，||·||2表示向量的L2范数运算符号，表示向量的L2范数平方运算符号，←表示赋值运算符号；最后将观测场景目标空间散射系数向量F(k)转换成三维矩阵形式，得到线阵SAR观测场景目标空间Ω的三维成像结果。

Claims (1)

1.基于半正定规划的线阵SAR后向投影自聚焦成像方法，其特征是它包括如下步骤：
步骤1、初始化线阵SAR系统参数和回波数据：
初始化线阵SAR系统参数包括：运动平台的速度矢量，记做V＝[vx,vy,vz]，其中vx为平台在水平面横轴向上的速度，vy为平台在水平面纵轴向上的速度，vz为平台在地面垂直高度向上的速度；线阵天线各阵元的初始位置矢量，记做Pn(0)＝[xn(0),yn(0),zn(0)]，其中n为天线阵元序号，为自然数，n＝1,2,...,N，N为线阵天线的阵元总数，xn(0)为第n个天线阵元在水平面横轴向上的初始位置，yn(0)为第n个天线阵元在水平面纵轴向上的初始位置，zn(0)为第n个天线阵元在地面垂直高度向上的初始位置；线阵天线的总长度，记做L；雷达工作的中心频率，记做fc；雷达工作的载频波长，记做λ；雷达发射基带信号的信号带宽，记做Br；雷达发射信号的脉冲宽度，记做TP；雷达发射信号的调频斜率，记做fdr；雷达接收系统的采样频率，记做fs；雷达系统的脉冲重复频率，记做PRF；雷达系统的脉冲重复时间，记为PRI；天线在方位向的有效孔径长度，记做Da；光在空气中的转播速度，记做c；距离向快时刻的总数，记为T，距离向快时刻序列记为t＝1,2,…,T，其中t为第t个距离向快时刻；方位向慢时刻的总数，记为K，方位向慢时刻序列记为l＝1,2,…,K，其中l为第l个方位向慢时刻；上述参数均为线阵SAR系统的标准参数，其中线阵天线的阵元总数N，线阵天线总长度L，雷达中心频率fc，雷达载频波长λ，雷达发射基带信号的信号带宽Br，雷达发射信号脉冲宽度TP，雷达发射信号调频斜率fdr，雷达接收系统的采样频率fs,雷达系统的脉冲重复频率PRF，天线在方位向的有效孔径长度Da在线阵SAR系统设计过程中已经确定；平台速度矢量V及线阵天线各阵元初始位置矢量Pn(0)在线阵SAR观测方案设计中已经确定；根据线阵SAR成像系统方案和观测方案，线阵SAR成像方法需要的初始化成像系统参数均为已知；
线阵SAR第n阵元在第l个方位向慢时刻和第t个距离向快时刻采样得到的原始回波数据，记为s(t,l,n)，t＝1,2,…,T，l＝1,2,…,K，n＝1,2,...,N；在实际线阵SAR系统中，原始回波数据s(t,l,n)可由线阵SAR系统数据接收机提供；在仿真线阵SAR成像过程中，原始回波数据s(t,l,n)可根据线阵SAR成像系统参数，采用标准合成孔径雷达原始回波仿真方法产生得到；在线阵SAR数据进行成像之前，原始回波数据s(t,l,n)均已知；
步骤2、初始化线阵SAR的观测场景目标空间参数：
初始化线阵SAR的观测场景目标空间参数，包括：以雷达波束照射场区域地平面和垂直于该地平面向上的单位向量所构成的三维空间直角坐标作为线阵SAR的观测场景目标空间，记为Ω；将观测场景目标空间Ω均匀划分成大小相等的立体单元格，也称为分辨单元，单元网格在水平横向、水平纵向和高度向边长分别记为dx、dy和dz，分辨单元在水平横向、水平纵向和高度向的总数分别记为Mx、My和Mz；观测场景目标空间Ω中第m个单元格的坐标矢量，记做qm，m表示观测场景目标空间Ω中第m个单元格，m为正整数，m＝1,2,…,M,其中M记为观测场景目标空间Ω中的单元格总数，且M＝Mx·My·Mz；
步骤3、线阵SAR原始回波信号距离压缩：
采用标准合成孔径雷达距离压缩方法对线阵SAR原始回波数据s(t,l,n)，t＝1,2,…,T，l＝1,2,…,K，n＝1,2,...,N进行距离压缩，得到距离压缩后的线阵SAR回波数据，记为sr(t,l,n)，t＝1,2,…,T，l＝1,2,…,K，n＝1,2,...,N，其中s(t,l,n)为步骤1中得到的线阵SAR原始回波数据，T为步骤1中得到的距离向快时刻总数，K为步骤1中得到的方位向慢时刻总数，N为步骤1中得到的线阵天线的阵元总数；
步骤4、计算观测场景单元对应的回波信号序号：
采用公式Pn(l)＝Pn(0)+V·l/PRF，n＝1,2,…,N，l＝1,2,…,K，计算得到第n个线阵天线阵元在第l个方位向慢时刻的位置矢量，记为Pn(l)，其中V为步骤1中初始化得到的平台速度矢量，Pn(0)为步骤1中初始化得到的线阵天线各阵元初始位置矢量，PRF为步骤1中初始化得到的雷达系统的脉冲重复频率，N为步骤1中得到的线阵天线阵元总数，K为步骤1中得到的方位向慢时刻总数；
采用公式R(Pn(l),qm)＝||Pn(l)-qm||2，n＝1,2,…,N，l＝1,2,…,K，m＝1,2,…,M，计算得到在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格到第n个线阵天线阵元的距离，记为R(Pn(l),qm)，其中||·||2表示向量的L2范数，qm为步骤2中初始化得到观测场景目标空间Ω中第m个单元格的坐标矢量，M为步骤2中初始化的场景目标空间Ω中单元格总数；
采用公式I(l,n,m)＝ceil(2·R(Pn(l),qm)·fs/c)，n＝1,2,…,N，l＝1,2,…,K，m＝1,2,…,M，计算得到在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格在第n个线阵天线阵元对应的原始回波序号，记为I(l,n,m)，其中c为步骤1中初始化得到的光在空气中的传播速度，fs为步骤1中初始化得到的雷达接收系统的采样频率，ceil(·)为上取整数运算符号；
采用公式Ir(l,n,m)＝2·R(Pn(l),qm)·fs/c-I(l,n,m)，n＝1,2,…,N，l＝1,2,…,K，m＝1,2,…,M，计算得到在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格在第n个线阵天线阵元对应的原始回波序号残余量，记为Ir(l,n,m)；
步骤5、利用sinc函数进行回波数据插值：
设定sinc插值函数的窗口长度，记为W0，其中W0为大于2的偶数；令w为整数，w的值为w＝-W0/2,-W0/2+1,…,W0/2-1,W0/2；
采用公式计算得到sinc函数插值后在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格在第n个线阵天线阵元对应的回波数据，记为sc(l,n,m)，其中sr(I(l,n,m)+w,l,n)为回波数据sr(t,l,n)中距离向快时刻t＝I(l,n,m)+w时对应的值，sr(t,l,n)为步骤3得到的距离压缩后线阵SAR回波数据，I(l,n,m)为步骤4得到的在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格在第n个线阵天线阵元对应的原始回波序号，Ir(l,n,m)为步骤4得到的在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格在第n个线阵天线阵元对应的原始回波序号残余量，sinc(·)为sinc函数运算符号，为元素w值从-W0/2至W0/2范围内的函数求和符号；
步骤6、构造线阵SAR后向投影积累基函数矩阵：
距离压缩后线阵SAR回波数据的后向投影积累基函数矩阵，记为B，矩阵B的维数大小M×(K·N)，其中M为步骤2中得到的观测场景目标空间Ω中的单元格总数，K为步骤1中得到的方位向慢时刻总数，N为步骤1中得到的线阵天线的阵元总数；矩阵B中第m行第g列的元素值，记为B(m,g)，其中m＝1,2,…,M，g为自然数，g的取值为g＝1,2,…,KN；
采用公式B(m,(l-1)·N+n)＝sc(l,n,m)·exp(j·4π·R(Pn(l),qm)/λ)，n＝1,2,…,N，l＝1,2,…,K，m＝1,2,…,M，计算得到矩阵B第m行第(l-1)·N+n列的元素值，记为B(m,(l-1)·N+n)，其中sc(l,n,m)为步骤5得到的sinc函数插值后在第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格在第n个线阵天线阵元对应的回波数据，R(Pn(l),qm)为步骤4得到的第l个方位向慢时刻线阵SAR观测场景目标空间Ω中第m个单元格到第n个线阵天线阵元的距离，λ为步骤1初始化得到的雷达工作的载频波长，exp(·)为自然常数e为底的指数运算符号，j为虚数符号，π为圆周率；
步骤7、选择观测场景目标空间中的强散射单元：
设定选择线阵SAR观测场景目标空间Ω中强散射单元的幅度阈值，记为ε0；
采用公式f＝B·IK·N计算得到线阵SAR观测场景散射系数的后向投影成像向量，记为f，向量f的维数大小为M×1，向量f的第m个元素记为fm，m＝1,2,…,M，其中B为步骤6得到的后向投影积累基函数矩阵，向量IK·N记为(K·N)×1的列向量，且向量IK·N中所有元素都为1；
采用公式β＝{m||fm|≥ε0·max(|f|),m＝1,2,…,M}计算得到观测场景目标空间Ω中强散射单元所在的序号向量，记为β，向量β的第y个元素记为βy，其中y为自然数，y的取值为y＝1,2,…,Q，Q为向量β中元素的总个数，fm为向量f的第m个元素，|·|表示为取绝对值运算符号，m||f(m)|≥ε0·max(|f|)表示满足条件|f(m)|≥ε0·max(|f|)所对应的序号m，max(·)表示为取最大值运算符号；
在矩阵B中选择向量β中元素值所对应的行向量按照序号大小组成矩阵，得到观测场景目标空间Ω中强散射单元的后向投影积累基函数矩阵，记为BS，矩阵BS的维数大小为Q×(K·N)，其中B为步骤6得到的后向投影积累基函数矩阵；
步骤8、初始化相位误差迭代估计的参数：
初始化相位误差迭代估计的参数包括：迭代估计过程的最大迭代次数，记做MaxIter；k记为迭代估计过程的第k次迭代，k为自然数，k初始值设置为k＝0，并且k的取值范围为k＝0,1,2,…,MaxIter；第k次迭代估计得到的相位误差向量，记为φ(k)，k＝0,1,2,…,MaxIter，其中，φ(k)为列向量，φ(k)的维数大小为(K·N)×1；初始化第k次迭代估计得到的相位误差向量的值为φ(k)＝0，k＝0,1,2,…,MaxIter；第k次迭代估计得到的线阵SAR观测场景后向投影成像结果向量，记为F(k)，k＝0,1,2,…,MaxIter，F(k)为列向量，向量F(k)的维数大小为M×1；初始化第k次迭代估计得到的线阵SAR观测场景后向投影成像结果向量的值为F(k)＝0，k＝0,1,2,…,MaxIter；
步骤9、利用半正定规划求解相位误差矩阵：
在第k次迭代中，若k＝0时采用公式γ(0)＝exp(jφ(0))计算得到相位误差迭代估计过程中第0次迭代的相位误差指数向量，记为γ(0)，其中φ(0)为步骤8初始化得到的第0次迭代估计相位误差向量，exp(·)为自然常数e为底的指数运算符号；若k≥1时采用公式γ(k)＝exp(jφ(k-1))计算得到相位误差迭代估计过程中第k次迭代的相位误差指数向量，记为γ(k)，其中φ(k-1)为第k-1次迭代估计得到的相位误差向量；
采用公式计算得到第k次迭代中强散射单元对应的相位误差半正定矩阵，记为C(k)，k＝0,1,2,…,MaxIter，其中上标H表示为共轭转置运算符号，BS为步骤7得到的观测场景目标空间Ω中强散射单元的后向投影积累基函数矩阵，为矩阵BS的共轭转置矩阵，diag(·)表示将向量元素值赋予对角矩阵中的对角元素运算符号
采用标准半正定规划计算以下公式
X ( k ) = arg min X ( t r ( C ( k ) · X ) )
s.t X(g,g)＝1,g＝1,…,K·N
X≥0
得到第k次迭代中的相位误差矩阵，记为X(k)，矩阵X(k)的维数大小为(K·N)×(K·N)，其中表示求取满足括号中最小值时对应的自变量X的最优值，tr(·)表示为矩阵求秩运算符，s.t表示为受约束符号，X(g,g)表示为矩阵X中第g行第g列的元素值，K为步骤1初始化得到的方位向慢时刻总数，N为步骤1初始化得到的线阵天线的阵元总数；
步骤10、估计相位误差向量：
对于步骤9得到的相位误差矩阵X(k)，如果矩阵X(k)满足tr(X(k))＝1，采用公式X(k)＝α·αH对矩阵X(k)进行向量分解，分解得到的列向量记为α，列向量α的维数大小为(K·N)×1，然后采用公式φ(k)＝∠α得到第k次迭代过程的相位误差向量，其中∠表示为求解角度运算符号，其中tr(·)表示为矩阵求秩运算符；
如果步骤9得到的矩阵X(k)满足tr(X(k))≠1，采用公式X(k)＝D·DH对矩阵X(k)进行平方矩阵分解，分解得到的矩阵记为D，矩阵D的维数大小为(K·N)×MS，矩阵D的第r个列向量记为dr，r为自然数，r的取值为r＝1,2,…,Ms；令Mr为正整数，Mr的值要大于Ms；随机产生Mr个均值为0、方差为IK·N的独立同分布复高斯向量，分别记为其中u1为第1个复高斯向量，u2为第2个复高斯向量，为第Mr个复高斯向量，IK·N为步骤7得到的元素都为1的向量；令w为自然数，w的取值为采用公式ηw＝exp(j∠D·uw)，w＝1,2,…,Mr，计算得到相位向量，记为ηw；采用公式计算得到最优的相位向量，记为ηopt，其中表示为在序号为1≤w≤Mr范围中求取满足中括号目标函数最小值时对应的向量序号，C(k)为步骤9得到的第k次迭代中强散射单元对应的相位误差半正定矩阵；然后，采用公式φ(k)＝∠ηopt得到第k次迭代过程的相位误差向量；
步骤11、线阵SAR观测场景目标空间后向投影成像：
采用公式计算得到第k次迭代过程的相位误差指数补偿向量，记为向量为列向量，的维数大小为(K·N)×1，其中φ(k)为步骤10得到的第k次迭代过程的相位误差向量，exp(·)为自然常数e为底的指数运算符号；
采用公式计算得到第k次迭代过程中线阵SAR观测场景后向投影成像结果，其中B为步骤6得到的线阵SAR后向投影积累基函数矩阵；
步骤12、迭代终止判定：
如果且k≤MaxIter，则k←k+1，执行步骤9至步骤12，否则终止迭代，此刻第k次迭代得到的F(k)即为线阵SAR观测场景最终的后向投影成像结果，其中k表示迭代估计过程中的第k迭代次数，MaxIter为步骤8中初始化得到的算法重构处理的最大迭代次数，F(k)为步骤11得到的第k次迭代过程中观测场景的后向投影成像结果，F(k-1)为第k-1次迭代过程中观测场景的后向投影成像结果，||·||2表示向量的L2范数运算符号，表示向量的L2范数平方运算符号，←表示赋值运算符号；最后将观测场景目标空间散射系数向量F(k)转换成三维矩阵形式，得到线阵SAR观测场景的最终三维成像结果。
CN201510231100.XA 2015-05-08 2015-05-08 基于半正定规划的线阵sar后向投影自聚焦成像方法 Expired - Fee Related CN104833973B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201510231100.XA CN104833973B (zh) 2015-05-08 2015-05-08 基于半正定规划的线阵sar后向投影自聚焦成像方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201510231100.XA CN104833973B (zh) 2015-05-08 2015-05-08 基于半正定规划的线阵sar后向投影自聚焦成像方法

Publications (2)

Publication Number Publication Date
CN104833973A CN104833973A (zh) 2015-08-12
CN104833973B true CN104833973B (zh) 2017-05-10

Family

ID=53811973

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201510231100.XA Expired - Fee Related CN104833973B (zh) 2015-05-08 2015-05-08 基于半正定规划的线阵sar后向投影自聚焦成像方法

Country Status (1)

Country Link
CN (1) CN104833973B (zh)

Families Citing this family (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105093225A (zh) * 2015-08-25 2015-11-25 西安电子科技大学 基于双重稀疏约束的逆合成孔径雷达自聚焦成像方法
CN105223573B (zh) * 2015-10-10 2018-11-13 电子科技大学 宽域高分辨率多目标逆合成孔径雷达成像技术
CN105487052B (zh) * 2015-12-08 2017-10-17 电子科技大学 基于低相干性的压缩感知lasar稀布线阵优化方法
CN107015225B (zh) * 2017-03-22 2019-07-19 电子科技大学 一种基于自聚焦的sar平台初始高度误差估计方法
CN107037429B (zh) * 2017-04-17 2020-06-16 电子科技大学 基于门限梯度追踪算法的线阵sar三维成像方法
CN107728144B (zh) * 2017-10-10 2020-06-16 电子科技大学 一种基于前视双基模式的干涉sar成像方法
CN107748362A (zh) * 2017-10-10 2018-03-02 电子科技大学 一种基于最大锐度的线阵sar快速自聚焦成像方法
CN107589421B (zh) * 2017-10-31 2022-03-29 西安电子科技大学 一种阵列前视sar成像方法
CN108872981B (zh) * 2018-04-20 2020-07-17 中国人民解放军国防科技大学 一种mimo雷达正则增强成像方法
CN111537999B (zh) * 2020-03-04 2023-06-30 云南电网有限责任公司电力科学研究院 一种稳健高效的分解投影自动聚焦方法
CN111679277B (zh) * 2020-05-28 2022-05-03 电子科技大学 一种基于sbrim算法的多基线层析sar三维成像方法
CN112946599B (zh) * 2021-02-04 2022-05-06 哈尔滨工业大学（威海） 一种基于稀布阵列的雷达空间谱估计方法
CN113608218B (zh) * 2021-07-19 2023-05-26 电子科技大学 一种基于后向投影原理的频域干涉相位稀疏重构方法
CN115061138B (zh) * 2022-08-18 2022-10-28 中南大学 探地雷达自聚焦并行掩膜后向投影成像方法、设备及介质

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103439693A (zh) * 2013-08-16 2013-12-11 电子科技大学 一种线阵sar稀疏重构成像与相位误差校正方法
CN103913741A (zh) * 2014-03-18 2014-07-09 电子科技大学 一种合成孔径雷达高效自聚焦后向投影bp方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103439693A (zh) * 2013-08-16 2013-12-11 电子科技大学 一种线阵sar稀疏重构成像与相位误差校正方法
CN103913741A (zh) * 2014-03-18 2014-07-09 电子科技大学 一种合成孔径雷达高效自聚焦后向投影bp方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
"LINEAR ARRAY SAR IMAGING VIA COMPRESSED SENSING";S.-J. Wei et al.;《Progress In Electromagnetics Research》;20111231;第117卷;299-319 *
"SPARSE AUTOFOCUS RECOVERY FOR UNDER-SAMPLED LINEAR ARRAY SAR 3-D IMAGING";Shun-Jun Wei et al.;《Progress In Electromagnetics Research》;20131231;第140卷;43-62 *
"Synthetic Aperture Radar Autofocus via Semidefinite Relaxation";Kuang-Hung Liu et al.;《IEEE TRANSACTIONS ON IMAGE PROCESSING》;20130630;第22卷(第6期);2317-2326 *

Also Published As

Publication number Publication date
CN104833973A (zh) 2015-08-12

Similar Documents

Publication Publication Date Title
CN104833973B (zh) 基于半正定规划的线阵sar后向投影自聚焦成像方法
CN103713288B (zh) 基于迭代最小化稀疏贝叶斯重构线阵sar成像方法
CN103698763B (zh) 基于硬阈值正交匹配追踪的线阵sar稀疏成像方法
CN107037429B (zh) 基于门限梯度追踪算法的线阵sar三维成像方法
CN103941243B (zh) 一种基于sar三维成像的自旋式飞行器测高方法
CN103439693B (zh) 一种线阵sar稀疏重构成像与相位误差校正方法
CN105677942B (zh) 一种重复轨道星载自然场景sar复图像数据快速仿真方法
CN105699969B (zh) 基于广义高斯约束的最大后验估计角超分辨成像方法
CN104730520B (zh) 基于子孔径合成的圆周sar后向投影自聚焦方法
CN104950305B (zh) 一种基于稀疏约束的实波束扫描雷达角超分辨成像方法
CN104898118B (zh) 一种基于稀疏频点的三维全息成像的重建方法
CN104698457B (zh) 一种迭代曲面预测InSAR成像及高度估计方法
CN104391295A (zh) 一种图像熵最优的压缩传感sar稀疏自聚焦成像方法
Wang et al. TPSSI-Net: Fast and enhanced two-path iterative network for 3D SAR sparse imaging
CN106405552B (zh) 基于wvd—pga算法的sar雷达目标聚焦方法
CN107015225B (zh) 一种基于自聚焦的sar平台初始高度误差估计方法
CN107576961B (zh) 一种互质降采样间歇合成孔径雷达稀疏成像方法
CN107748362A (zh) 一种基于最大锐度的线阵sar快速自聚焦成像方法
CN107271993A (zh) 一种基于最大后验的扫描雷达角超分辨成像方法
CN104536000A (zh) 一种实波束扫描雷达角超分辨方法
CN108226927A (zh) 基于加权迭代最小稀疏贝叶斯重构算法的sar成像方法
CN104251991B (zh) 一种基于稀疏度估计的分维度阈值迭代稀疏微波成像方法
CN105137425A (zh) 基于卷积反演原理的扫描雷达前视角超分辨方法
CN102854505A (zh) 一种加权稀疏驱动自聚焦sar成像方法
CN107607945A (zh) 一种基于空间嵌入映射的扫描雷达前视成像方法

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
EXSB Decision made by sipo to initiate substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20170510

Termination date: 20200508

CF01 Termination of patent right due to non-payment of annual fee