CN102914773B - 一种多航过圆周sar三维成像方法 - Google Patents

一种多航过圆周sar三维成像方法 Download PDF

Info

Publication number
CN102914773B
CN102914773B CN201210333250.8A CN201210333250A CN102914773B CN 102914773 B CN102914773 B CN 102914773B CN 201210333250 A CN201210333250 A CN 201210333250A CN 102914773 B CN102914773 B CN 102914773B
Authority
CN
China
Prior art keywords
data
echo
vector
time
note
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
CN201210333250.8A
Other languages
English (en)
Other versions
CN102914773A (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 CN201210333250.8A priority Critical patent/CN102914773B/zh
Publication of CN102914773A publication Critical patent/CN102914773A/zh
Application granted granted Critical
Publication of CN102914773B publication Critical patent/CN102914773B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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/904—SAR modes
    • G01S13/9088—Circular SAR [CSAR, C-SAR]
    • 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/9011—SAR image acquisition techniques with frequency domain processing of the SAR signals in azimuth

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三维成像方法，它是通过对距离向数据采用距离压缩处理，并进行一维快速傅里叶变换（FFT）将数据转换到波数域，采用三维线性插得到直角坐标系下均匀分布的数据格式，构造出雷达回波信号在波数域的线性观测模型，采用加权的L1范数构造最优化目标函数，利用特征增强图像形成方法中的迭代算法求解出最优化解向量，从而得到目标场景的三维成像结果，实现了三维成像，有效的降低了旁瓣和提高了分辨率，获得良好的图像聚焦效果。

Description

一种多航过圆周SAR三维成像方法
技术领域
本发明属于雷达成像技术领域，它特别涉及了合成孔径雷达(SAR)成像技术领域。
背景技术
合成孔径雷达(SAR)是一种具有全天候、全天时工作能力的微波成像系统，它能够在环境恶劣的情况下得到较为优秀的成像结果，因此它在民用和军用方面的地位是不可替代的，是进行地形测绘和军事侦察等方面应用的有效手段。它利用雷达平台沿着航迹向运动形成虚拟天线阵列合成孔径来获得航迹向的分辨率，利用脉冲压缩技术获得斜距向分辨率。
三维合成孔径雷达(SAR)是一种新型合成孔径雷达技术。三维合成孔径雷达成像系统的基本原理是通过形成虚拟天线面阵以获得二维分辨力，再结合距离压缩技术得到距离向的聚焦能力，从而得到对成像区域的三维聚焦结果。SAR三维成像技术是SAR成像系统区别于其他遥测遥感系统的重要特征，由于其测绘范围广、具有三维成像能力等特点，在地形测绘、环境检测、灾害预测等方面具有广阔的应用场景。当前三维SAR成像技术主要有干涉SAR(InSAR)技术、线阵SAR(LASAR)技术和曲线SAR技术(CLSAR)。
圆周合成孔径雷达(CSAR)是在航迹向和切航迹向形成圆周孔径的一种特殊的曲线合成孔径雷达模式。该模式是利用天线平台绕着场景中心旋转飞行形成圆周轨迹，同时照射成像场景，再通过对收集到的回波数据进行处理来得到成像结果。由于圆周合成孔径雷达可以对场景进行360度观测，获得成像场景中散射点的全角度散射信息，该模式主要应用于对隐蔽目标的侦测以及特定区域的全方位侦查。文献Mehrdad Soumekh“Reconnaissance with Slant Plane Circular SARImaging”及Akira Ishimaru，Tsz-king Chan and Yasuo Kuga“An ImagingTechnique Using Confocal Circular Synthetic Aperture Radar”都对圆周SAR成像系统进行了分析，并提出了相应的成像算法。在实际的观测中，圆周合成孔径雷达方位向的相关角度无法达到360度的范围，根据实验可以得出圆周合成孔径雷达的方位向相关角度不超过20度，因此圆周合成孔径雷达三维成像能力较差。为了提高圆周合成孔径雷达的三维分辨率，Matthew Ferrara等人在文献“Enhancement of Multi-pass 3D Circular SAR Images Using SparseReconstruction Techniques”中对多航过圆周合成孔径雷达模式进行了分析。多航过圆周合成孔径雷达是雷达运动平台沿不同的高度平面进行多次圆周飞行，在高度向形成合成孔径以获得高度向的分辨率力的一种雷达运动模式。由于传统的圆周合成孔径雷达成像算法相干角度无法达到360度，划分子孔径的办法将会导致图像出现高旁瓣。
发明内容
本发明的目的是针对现有的圆周SAR中场景散射点无法在360度的观测角度内保持相关性，导致三维分辨率下降，并且划分子孔径的办法将会引起图像出现高旁瓣的缺点，提出了基于特征增强图像形成(Feature-Enhanced ImageFormation)的一种多航过圆周SAR三维成像方法，采用本发明的方法可以有效的降低旁瓣，提高分辨率。
为了方便描述本发明的内容，首先做以下术语定义：
定义1、圆周三维成像合成孔径雷达
圆周合成孔径雷达是利用运动平台绕着场景中心旋转飞行形成圆周轨迹，在航迹向和切航迹向形成圆周孔径，可以对观测区域进行三维成像的合成孔径雷达系统。
定义2、三维合成孔径雷达成像空间
三维合成孔径雷达成像空间是指运用成像算法由三维合成孔径雷达数据空间所得到的三维合成孔径雷达的图像空间。
定义3、三维合成孔径雷达数据空间
三维合成孔径雷达数据空间是指三维合成孔径雷达回波数据所构成的回波信号空间。
定义4、快时间与慢时间
慢时间是指收发平台飞过一个飞行孔径所需要的时间，由于雷达以一定的脉冲重复周期Tr发射接收脉冲，在第σ个发射脉冲时刻，慢时间可以表示为一个离散化的时间变量ts(σ)=σTr，σ=1，2，…，Na，Na表示为一个合成孔径内慢时间的离散个数。
快时间是指雷达发射接收脉冲的一个周期的时间，由于雷达接收回波是以采样率fs进行采样，在第ξ个采样时刻，快时间可以表示为一个离散化的时间变量tf(ξ)=ξ／fs，m=1，2，…，Nr，Nr表示为一个脉冲重复周期Tr内快时间的离散个数。
定义5、合成孔径雷达标准距离压缩方法
合成孔径雷达标准距离压缩方法是指利用合成孔径雷达发射信号参数，采用匹配滤波技术对合成孔径雷达的距离向信号进行滤波的过程。详见文献“雷达成像技术”，保铮等编著，电子工业出版社出版。
定义6、去调频方法
去调频方法又称为STRETCH处理法、宽带压缩法或时频转换法，是针对线性调频信号提出的一种信号处理方法，在接收端设置一个参考函数，以场景中心为参考点，将参考函数与接收信号在时域共轭相乘，实现去调频。
定义7、合成孔径雷达发射机
合成孔径雷达发射机是指目前合成孔径雷达采用的向观测区域发射电磁信号的系统，主要包括信号发生器、混频器、放大器等模块。
定义8、合成孔径雷达接收机
合成孔径雷达接收机是指目前合成孔径雷达采用的接收观测区域回波的系统，主要包括混频器、放大器、模／数转换器、存储设备等。
定义9、三维合成孔径雷达成像空间
三维合成孔径雷达成像空间是指运用成像算法由三维合成孔径雷达数据空间所得到的三维合成孔径雷达的图像空间。
定义10、diag(·)与ln(·)
diag(a1，a2，…，aM)表示对角元素为a1，a2，…，aM的M×M维方阵，其中a1，a2，…，aM均为数域C上的数。
定义11、向量范数
设映射||·||：Cn→R满足：
(1)正定条件，||x||≥0，当且仅当x=0时||x||=0；
(2)齐次条件，||λx||=|λ|||x||，λ∈C，x∈Cn；
(3)三角不等式， | | x + y | | ≤ | | x | | + | | y | | , ∀ x , y ∈ C n ,
则称映射||·||为Cn上向量x的范数。定义L1范数和L2范数为
Figure GDA0000428831100000042
| | x | | 2 = ( Σ i = 1 n | x i | 2 ) 1 / 2 .
定义12、直角坐标系和极坐标系
三维直角坐标系通常由三个相互垂直的坐标轴组成，通常三个坐标轴称为x轴、y轴和z轴，三个坐标轴交于一点称为原点。
极坐标系的坐标用
Figure GDA0000428831100000044
表示，其中ρ便是距离原点的距离，θ表示距离x轴的角度，表示距离z轴的角度。
极坐标系与直角坐标系的坐标对应关系为
Figure GDA0000428831100000047
定义13、快速傅里叶变换(FFT)
快速傅里叶变换是离散傅里叶变换的快速算法，离散傅里叶变换的公式为
Figure GDA0000428831100000048
详细内容请参考何子述等编写的“现代数字信号处理及其应用”
定义14、波数和波数域
波数是指单位长度内电磁波变化的相位。定义为ka=2πfc／v，其中fc为电磁波传播频率，v为电磁波传播速度。详细内容请参考皮亦鸣等编写的“合成孔径雷达成像原理”。
沿某方向传播的电磁波波数与该方向距离单位为一对傅里叶变换对。将雷达观测场景在空间直角坐标系的散射系数进行三维快速傅里叶变换(FFT)，可以得到散射场景在波数域的数据。
定义15、三维线性插值
线性插值为数学、计算机图像学等领域广泛应用的一种简单插值方法，该方法已知两点，求解两点连线上某点的值。详见Grant D.Martin等编写文献“SARPolar Format Implementation with MATLAB”。
定义16、转置和共轭转置
设矩阵A∈CN×M，其中C为复数域，N与M为正整数，对矩阵A的转置取共轭即可得到矩阵A的共轭转置矩阵AH，将H记做共轭转置运算符。对矩阵A取转置可得到矩阵A的转置矩阵，记做AT，T为转置运算符。
定义17、复指数函数
设x是复数域C上的复数exp(x)=ex，其中e为自然对数的底数。
定义18，纯虚数和π
纯虚数
Figure GDA0000428831100000051
π为圆周率取3.1415926。
定义19，向量的梯度
设一列向量w=[w1，w2，…wp]，其中wi=αi+jβi，J=J(w)是向量w的一个标量函数。则J(W)关于w的梯度定义为：
▿ J ( w ) = [ ∂ J ∂ α 1 + j ∂ J ∂ β 1 , ∂ J ∂ α 2 + j ∂ J ∂ β 2 , . . . , ∂ J ∂ α p + j ∂ J ∂ β p ] T , 其中为函数J对变量α的偏导数。具体定义可参考同济大学数学系主编的“高等数学”。
定义20、Matlab
Matlab是矩阵实验室(Matrix Laboratory)的简称，美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境。具体用法详见文献“MATLAB 5手册”，EvaPart-Enander等编著，机械工业出版社出版。
本发明提出了一种多航过圆周SAR三维成像方法，它包括以下几个步骤，如附图2所示：
1、一种多航过圆周SAR三维成像方法，其特征包括以下步骤：
步骤1、初始化圆周三维成像合成孔径雷达成像系统参数以及特征增强图像形成方法参数：
初始化成像系统参数包括：雷达发射电磁波的载频，记做fc；雷达发射电磁波的带宽，记做BW；雷达运动平台的飞行半径，记做R：雷达运动平台的方位角，记做θ：雷达运动平台的俯仰角，记做
Figure GDA0000428831100000061
电磁波的传播速度，记做C：雷达在方位维的观测范围，记做θw；合成孔径雷达距离向回波数据，记做
Figure GDA0000428831100000062
雷达快时间采样点数，记做K：雷达慢时间采样点数，记做NP；多航过雷达飞行模式航过次数，记做L：观测场景在距离向的观测范围，记做Wx；观测场景在方位向的观测范围，记做Wy；观测场景在高度向的观测范围，记做Wz；回波数据在波数域插值界限，记做B：观测场景像素划分间隔，记做Δw；观测场景的中心点，记做O；其中，雷达发射电磁波的载频，雷达发射电磁波的带宽，电磁波的传播速度在雷达系统的设计过程中均已确定，雷达运动平台的飞行半径，雷达运动平台的方位角，雷达运动平台的俯仰角，合成孔径雷达在方位维的观测范围，观测场景在距离向、方位向和高度向的观测范围，多航过雷达飞行模式航过次数，回波数据在波数域插值界限，观测场景像素划分间隔在预先确定观测方案中为均为已知；
初始化的特征增强图像形成方法参数包括：拉格朗日乘数，记做λ：特征增强图像形成算法误差门限，记做Δ：共轭梯度算法求解误差门限，记做ε：共轭梯度算法最大迭代次数，记做cgnum；特征增强图像形成方法最大迭代次数，记做itenum；初始解向量，记做f(0)；L1范数逼近常数，记做η；其中，拉格朗日常数特征增强图像形成算法误差门限Δ，共轭梯度算法求解误差门限ε.共轭梯度算法最大迭代次数cgnum，特征增强图像形成方法最大迭代次数itenum，L1范数逼近常数η，初始解向量f(0)均为已知；
步骤2、构建多航过圆周三维成像合成孔径雷达的运动模型
以场景中心点O为原点，以地面为xy平面，以高度向为z轴，构建空间直角坐标系，雷达运动平台以[0，0，z1]，[0,0，z2]，…，[0，0，zτ]，…，[0,0，zL]为圆心做半径为R的圆周运动，其中zτ为雷达运动平台飞行高度，合成孔径雷达运动平台不同高度平面编号为τ，取值范围是：τ=1，2，…，L，其中z1为第1次时合成孔径雷达雷达运动平台飞行高度，z2为第2次合成孔径雷达运动平台飞行高度，…，zτ为第τ次合成孔径雷达运动平台飞行高度，…，zL为第L次合成孔径雷达运动平台飞行高度；
步骤3、圆周三维成像合成孔径雷达原始数据进行距离压缩
采用合成孔径雷达标准距离压缩方法对步骤1中提供合成孔径雷达距离向回波数据进行压缩，得到距离压缩后的圆周三维成像合成孔径雷达数据，记做
Figure GDA0000428831100000072
距离压缩后的圆周三维成像合成孔径雷达数据
Figure GDA0000428831100000073
是一个三维数组，包括距离维、沿航迹维、高度维；
步骤4、获取圆周三维成像合成孔径雷达波数域的回波数据
对步骤3得到的距离压缩后的圆周三维成像合成孔径雷达数据中的每个距离维的数据进行一维快速傅里叶变换(FFT)，得到圆周三维成像合成孔径雷达波数域的回波数据，该数据是一个三维数组，记做
Figure GDA0000428831100000075
表示圆周三维成像合成孔径雷达波数域的回波数据，包括波数维、方位角维和俯仰角维，ρk表示合成孔径雷达回波数据在波数域中波数数据，θi表示合成孔径雷达回波数据在波数域中方位角数据，
Figure GDA0000428831100000076
表示合成孔径雷达回波数据在波数域中俯仰角数据，其中k表示波数域中波数方向采样编号，取值范围是：k=1，2，…，K，i表示波数域中方位角方向采样编号，取值范围是：i=1，2，…，NP，l表示波数域中俯仰角方向采样编号，取值范围是：l=1，2，…，L；
步骤5、通过插值获得直角坐标系下波数域回波数据，包括步骤5.1、步骤5.2；
步骤5.1
取步骤4中得到的圆周三维成像合成孔径雷达波数域的回波数据
Figure GDA0000428831100000077
的中心点，记做
Figure GDA0000428831100000078
其中ρ0=(ρ1+ρK)／2，
Figure GDA0000428831100000079
ρ1为步骤4中k=1时的波数，ρk为步骤4中k=K时的波数，θ1为步骤4中i=1时的方位角，
Figure GDA00004288311000000710
为步骤4中i=NP时的方位角，
Figure GDA0000428831100000081
为步骤4中i=1时的俯仰角，为步骤4中l=L时的俯仰角；
采用坐标系变换方法进行坐标系变换：利用直角坐标系和极坐标系的对应关系构建三维空间波数域的直角坐标系，得到空间直角坐标系；其中的坐标轴记做kx、ky和kz；以
Figure GDA0000428831100000083
为中心点，在kx轴方向以2π／Wx为间隔，在ky轴方向以2π／Wy为间隔，在kz轴方向以2π／Wz为间隔，在空间直角坐标系下将三维空间波数域离散为M×N×P个单元格，其中，M为kx轴方向单元格数目，N为ky轴方向单元格数目，P为kz轴方向单元格数目，其中M，N，P均为自然数，并满足 M ≤ BW x 2 π , N ≤ BW y 2 π , P ≤ BW z 2 π ;
步骤5.2
利用三维线性插值方法，将步骤4得到的回波数据
Figure GDA0000428831100000085
插值到步骤5.1得到的直角坐标系下M×N×P个单元格内，得到直角坐标系下的波数域回波数据，记做
步骤6、构造回波数据的线性测量矩阵和散射场景的线性观测模型，包括步骤6.1、步骤6.2、步骤6.3、步骤6.4；
步骤6.1
将步骤5得到的直角坐标系下的波数域回波数据用
Figure GDA0000428831100000088
表示直角坐标系下波数域回波数据在坐标(kx，m，ky，n，kz，μ)处的数据，其中，m表示步骤5得到的空间直角坐标系坐标轴kx方向采样编号，取值范围是：m=1.2…，M，n表示步骤5得到的空间直角坐标系坐标轴ky方向采样编号，取值范围是：n=1,2…，N，μ表示步骤5得到的空间直角坐标系坐标轴kz方向采样编号，取值范围是：μ=1,2…，P；
步骤6.2
将步骤5得到的直角坐标系下的波数域回波数据
Figure GDA0000428831100000089
定义为一个回波向量，记做
Figure GDA0000428831100000091
其中
Figure GDA0000428831100000092
为回波向量
Figure GDA0000428831100000093
的第1行的回波数据，
Figure GDA0000428831100000094
为回波向量
Figure GDA0000428831100000095
的第2行的回波数据，…，
Figure GDA00004288311000000937
为回波向量
Figure GDA0000428831100000096
的第μ行的回波数据，…，
Figure GDA0000428831100000097
为回波向量的第P行的回波数据；
回波向量
Figure GDA0000428831100000099
第μ行回波数据
Figure GDA00004288311000000910
其中
Figure GDA00004288311000000911
为回波向量
Figure GDA00004288311000000912
第1行的数据，
Figure GDA00004288311000000913
为回波向量
Figure GDA00004288311000000914
的第2行数据，…，
Figure GDA00004288311000000915
为回波向量
Figure GDA00004288311000000916
第n行数据，…，为回波向量
Figure GDA00004288311000000918
第N行的数据；
Figure GDA00004288311000000919
的第n行数据 y ‾ μ , n = [ y μ , n , 1 , y μ , n , 2 , . . . , y μ , n , m , . . . , y μ , N , M ] T , 其中yμ，n，1为回波向量
Figure GDA00004288311000000922
在波数域回波数据
Figure GDA00004288311000000923
中坐标(kx，1，ky，nkz，u)处的数据，uμ，n，2为回波向量在波数域回波数据中坐标(kx，2，ky，nkz，u)处的数据，…，yμ，n，m为回波向量
Figure GDA00004288311000000926
在波数域回波数据
Figure GDA00004288311000000927
中坐标(kx，m，ky，nkz，u)处的数据，…，yμ，N，M为回波向量
Figure GDA00004288311000000928
在波数域回波数据
Figure GDA00004288311000000929
中坐标(kx，M，ky，nkz，u)处的数据；
将yμ，n，m记做yq，则回波向量
Figure GDA00004288311000000930
表示为其中y1为回波向量
Figure GDA00004288311000000932
第1行的回波数据，y2为回波向量
Figure GDA00004288311000000933
第2行的回波数据，…，yq为回波向里
Figure GDA00004288311000000934
第q行的回波数据，…，yQ为回波向量
Figure GDA00004288311000000938
第Q行的回波数据，其中q为回波向量
Figure GDA00004288311000000935
的行标号，取值范围是：q=1，2，…，Q，Q=M×N×P；yq在波数域对应的位置坐标记做
Figure GDA00004288311000000936
步骤6.3
根据步骤2中构建的空间直角坐标系，将观测场景以场景中心O为中心点，在x轴方向以Δw为间隔，将x轴方向的观测场景范围Wx等分为x份；在y轴方向以Δw为间隔，将y轴方向的观测场景范围Wy等分为Y份；在z轴方向以Δw为间隔，将z轴方向的观测场景范围Wz等分为Z份；通过上述方法观测场景被划分为X×Y×Z个单元格，其中X为距离向单元格数目，Y为方位向单元格数目，Z为高度向单元格数目，且满足
Figure GDA0000428831100000101
每个单元格在直角坐标系下的位置矢量记做
Figure GDA0000428831100000102
其中，f表示空间直角坐标系坐标轴x轴方向采样编号，取值范围是：f=1,2，…X；g表示空间直角坐标系坐标轴y轴方向采样编号，取值范围是：g=1.2，…，Y；h表示空间直角坐标系坐标轴z轴方向采样编号，取值范围是：h=1,2，…，z；
构造测量矩阵为其中，
Figure GDA0000428831100000104
为测量矩阵
Figure GDA0000428831100000105
的第1行数据，
Figure GDA0000428831100000106
为测量矩阵
Figure GDA0000428831100000107
的第2行数据，…，
Figure GDA0000428831100000108
为测量矩阵
Figure GDA0000428831100000109
的第t行数据，…，
Figure GDA00004288311000001010
为测量矩阵
Figure GDA00004288311000001011
的第Q行数据，其中t为测量矩阵
Figure GDA00004288311000001012
的行标号，t=1,2，…，Q，Q=M×N×P；
测量矩阵
Figure GDA00004288311000001013
的第t行数据为 α ‾ t = [ α t ( 1,1,1 ) , α t ( 1,1,2 ) , . . . , α t ( f , g , h ) , . . . , α t ( X , Y , Z ) ] , 其中αt(1，1,1)为测量矩阵
Figure GDA00004288311000001015
的第t行第一列数据，αt(1，1，2)为为测量矩阵
Figure GDA00004288311000001016
的第t行第二列数据，αt(f，g，h)=exp{j(kx，fxf+ky，tyg+kz，tzh)}，(kx，f，ky，t，kz，t)为步骤6.3中回波向量
Figure GDA00004288311000001017
第t行数据对应的波数域位置坐标，(xf，yg，zh)为观测场景的单元格位置坐标，αt(X，Y，Z)为f=X，g=Y，h=Z时的值；
步骤6.4
随机抽取步骤6.2中回波向量中的第N1行作为抽取后的回波向量y的第1行，抽取回波向量
Figure GDA0000428831100000112
中的第N2行作为抽取后的回波向量y的第2行，…，抽取回波向量中的第Nr行作为抽取后的回波向量y的第r行，…，同理，抽取回波向量
Figure GDA0000428831100000114
中的第Nnum行作为抽取后的回波向量y的第nun行，N1为第1次抽取回波向量
Figure GDA0000428831100000115
的行标，N2为第2次抽取回波向量
Figure GDA0000428831100000116
的行标，…，Nr为第r次抽取回波向量的行标，…，Nnum为第num次抽取回波向量
Figure GDA0000428831100000118
的行标，其中r为抽取后的回波向量y的行标，取值范围是：r=1，2，…，num，num为抽取行的总数，并且N1，N2，…，Nnum，num均为不大于Q的自然数；
同理，抽取步骤6.3中测量矩阵
Figure GDA0000428831100000119
中的第N1行作为抽取后的测量矩阵A的第1行，抽取测量矩阵中的第N2行作为抽取后的测量矩阵A的第2行，…，抽取测量矩阵
Figure GDA00004288311000001111
中的第Nr行作为抽取后的测量矩阵A的第r行，…，同理，抽取测量矩阵
Figure GDA00004288311000001112
中的第Nnum行作为抽取后的测量矩阵A的第mum行，N1为第1次抽取测量矩阵的行标，N2为第2次抽取测量矩阵
Figure GDA00004288311000001114
的行标，…，Nr为第r次抽取测量矩阵
Figure GDA00004288311000001115
的行标，…，Nnum为第num次抽取测量矩阵
Figure GDA00004288311000001116
的行标；
步骤7、改进特征增强图像形成方法迭代求解
第1次求解，将步骤6中抽取后的回波向量y和抽取后的测量矩阵A带入公式(1)中求解：
f(1)=(2AHA+λΛ(f(0)))-12AHy                      (1)
公式(1)中： Λ ( f ( 0 ) ) = diag { 1 ( | f 1 ( 0 ) + η | ) 2 , 1 ( | f 2 ( 0 ) + η | ) 2 , . . . , 1 ( | f Q ( 0 ) + η | ) 2 } , f1 (0)为步骤1中初始解向量f(0)的第1行数据，f2 (0)为步骤1中初始解向量f(0)的第2行数据，fQ (0)为步骤1中初始解向量f(0)的第Q行数据，Λ(f(0))为第1次求解使用的加权L1范数的梯度，f(1)为第1次求解得出的解向量；
判断
Figure GDA0000428831100000121
是否成立，若成立则求解结束；若不成立，则判断求解次数是否达到第itenum次，若本次求解为第itenum求解则求解结束，否则继续进行第2次求解；
第2次求解，将步骤6中抽取后的回波向量y和抽取后的测量矩阵A带入公式(2)中求解：
f(2)=(2AHA+λΛ(f(1)))-12AHy                     (2)
公式(2)中： Λ ( f ( 1 ) ) = diag { 1 ( | f 1 ( 1 ) + η | ) 2 , 1 ( | f 2 ( 1 ) + η | ) 2 , . . . , 1 ( | f Q ( 1 ) + η | ) 2 } , f1 (1)为第1次求解得出的解向量f(1)的第1行数据，f2 (1)为第1次求解得出的解向量f(1)的第2行数据，fQ (1)为第1次求解得出的解向量f(1)的第Q行数据，Λ(f(1))为第2次求解使用的加权l1范数的梯度，f(2)为第2次求解得出的解向量；
判断
Figure GDA0000428831100000124
是否成立，若成立则结束求解；若不成立，判断求解次数是否达到第itrnum次，若本次求解为第itrnum求解则求解结束，否则进行第3次求解；
同理，第J次求解，将步骤6中抽取后的回波向量y和抽取后的测量矩阵A带入公式(J)中求解：
f(J)=(2AHA+λΛ(f(J-1)))-12AHy                 (J)
公式(J)中： Λ ( f ( J - 1 ) ) = diag { 1 ( | f 1 ( J - 1 ) + η | ) 2 , 1 ( | f 2 ( J - 1 ) + η | ) 2 , . . . , 1 ( | f Q ( J - 1 ) + η | ) 2 } , f1 (J-1)为第J-1次求解得出的解向量f(J-1)的第1行数据，f2 (J-1)为第J-1次求解得出的解向量f(J-1)的第2行数据，fQ (J-1)为第J-1次求解得出的解向量f(J-1)的第Q行数据，Λ(f(J-1))为第J次求解使用的加权L1范数的梯度，f(J)为第J次求解得出的解向量，J表示求解编号，且满足1≤j≤itenum；
判断
Figure GDA0000428831100000131
是否成立，若成立则结束求解；若不成立，判断求解次数是否达到第itenum次，若本次求解为第itenum求解则求解结束，否则进行第J+1次求解；
经过上述步骤，得到圆周三维成像合成孔径雷达图像数据。
需要说明的是，本发明方法中的利用公式(1)、(2)…(J)求解出的解向量即为合成孔径雷达图像数据。
本发明的创新点：针对现有圆周三维成像合成孔径雷达成像处理算法在成像过程中无法保证散射点在方位观测角度范围内保持相关性，并且划分子孔径的做法会导致图像出现旁瓣串扰降低成像的分辨率，而在目标侦查等方面，对测量精度提出了很高的要求，本发明提出了一种基于特征增强图像形成的圆周三维成像合成孔径雷达成像方法。该方法将成像问题转化为一个最优化问题，并对原算法中最优化函数提出改进，增加了加权矩阵和L1范数近似方法，从而获得更好的图像聚焦效果，有效的降低了旁瓣提高了分辨率。
本发明的基本原理：本方法对距离向数据采用距离压缩处理，并对其进行一维快速傅里叶变换(FFT)将数据转换到波数域，通过三维线性插得到直角坐标系下均匀分布的数据格式，构造出雷达回波信号在波数域的线性观测模型，通过采用加权的L1范数构造最优化目标函数，利用特征增强图像形成方法的迭代算法求解出最优化解向量，从而得到目标场景的三维成像结果，实现了三维成像，可以有效的降低旁瓣，提高分辨率。
本发明优点：针对圆周三维成像合成孔径雷达在全观测角度内无法保持场景散射点的相关性，而且划分子孔径的做法将导致图像出现高旁瓣降低分辨率，提出了基于特征图像增强图像形成方法，该方法将原始算法中的L1范数用加权L1范数代替，获得了更加稀疏的解向量，从而降低了旁瓣，提高了分辨率。与传统的基于FFT的成像算法(BP算法、PFA算法等)，克服了划分子孔径会导致图像出现高旁瓣的问题，从而提高了分辨率，改善了成像质量。
附图说明
图1多航过圆周三维成像合成孔径雷达模型
其中R为飞行半径，z1为第一次飞行高度，z2为第二次飞行高度，zL为第L次飞行高度，O为场景中心点，x，y，z为空间直角坐标系坐标轴标号
图2本发明的结构流程图
图3多航过圆周三维成像合成孔径雷达系统参数及特征化图像增强方法初始化参数
具体实施方式
本发明主要采用仿真实验的方式进行验证该系统模型的可行性，所有步骤、结论都在MATLAB7.0上验证正确。具体实施步骤如下：
步骤1：系统参数设定及初始化特征图像增强方法的参数
本具体实施方式所采用的系统参数详见图3。
步骤2、构建多航过圆周三维合成孔径雷达的运动模型
以场景中心点O为原点，以地面为xy平面，以高度向为z轴，构建空间直角坐标系，雷达运动平台以[0，0，z1]，[0，0，z2]…[0，0，z6]为圆心做半径为R=1000m的圆周运动，其中z1=1414m，z2=1466m，z3=1524m，z4=1589m，z5=1662m，z6=1743m。
步骤3、圆周三维成像合成孔径雷达原始数据进行距离压缩
采用传统的合成孔径雷达距离压缩方法对合成孔径雷达距离向回波数据进行压缩，得到距离压缩后的圆周三维成像合成孔径雷达数据，记做
Figure GDA0000428831100000142
距离压缩后的圆周三维成像合成孔径雷达数据是一个三维数组，包括距离维、沿航迹维、高度维；
步骤4、获取圆周三维成像合成孔径雷达波数域的回波数据
将步骤3得到的距离压缩后的三维回波数组
Figure GDA0000428831100000151
对每个距离维的数据进行传统的一维快速傅里叶变换(FFT)，得到圆周三维成像合成孔径雷达波数域的回波数据，该数据是一个三维数组，记做
Figure GDA0000428831100000152
表示极坐标下三维空间波数域的数据，包括波数维、方位角维和俯仰角维。用
Figure GDA0000428831100000153
表示合成孔径雷达回波数据在波数域中波数ρk，方位角θi，俯仰角
Figure GDA0000428831100000154
处的数据，其中k=1，2，…K，i=1，2，…Np，l=1，2，…L，K=512为步骤1中快时间采样点数，Np=128为步骤1中慢时间采样点数，L=6为步骤1中多航过雷达飞行模式航过次数。步骤5、通过插值获得直角坐标系下波数域回波数据，包括步骤5.1，步骤5.2。
步骤5.1
取步骤4中得到的极坐标下波数域回波数据
Figure GDA0000428831100000155
的中心点，记做
Figure GDA0000428831100000156
其中ρ0=(ρ1+ρ512)／2，θ0=(θ1+θ128)／2，
Figure GDA0000428831100000157
采用传统坐标系变换方法进行坐标系变换，利用直角坐标系和极坐标系的对应关系构建三维空间波数域的直角坐标系，得到空间直角坐标系，其中的坐标轴记做kx、ky和kz。
以
Figure GDA0000428831100000158
为中心点，在kx轴方向以2π／Wx为间隔，在ky轴方向以2π／Wy为间隔，在kz轴方向以2π／Wz为间隔，在直角坐标系下将三维空间波数域离散为M×N×P个单元格，其中，Wx=2m为步骤1中初始化参数中观测场景在距离向的观测范围，Wy=2m为步骤1中初始化参数中观测场景在方位向的观测范围，Wz=1m为步骤1中初始化参数中观测场景在高度向的观测范围，M=20为kx轴方向单元格数目，N=20为ky轴方向单元格数目，P=10为kz轴方向单元格数目。
步骤5.2
利用传统的三维线性插值方法，将步骤4得到的回波数据
Figure GDA0000428831100000159
插值到步骤5.1得到的直角坐标系下20×20×10个单元格内，得到直角坐标系下的波数域回波数据，记做
Figure GDA0000428831100000161
步骤6、构造回波数据的线性测量矩阵和散射场景的线性观测模型，包括步骤6.1，步骤6.2，步骤6.3，步骤6.4。
步骤6.1
将步骤5得到的直角坐标系下回波数据
Figure GDA0000428831100000162
用
Figure GDA0000428831100000163
表示直角坐标系下波数域回波数据在坐标(kx，m，k，y，nk，处的数据，其中
Figure GDA0000428831100000164
步骤6.2
将步骤5得到的空间直角坐标系下的波数域回波数据定义为一个回波向量，记做
Figure GDA0000428831100000166
其中
Figure GDA0000428831100000167
为回波向量的第1行的回波数据，
Figure GDA0000428831100000169
为回波向量
Figure GDA00004288311000001610
的第2行的回波数据，…，为回波向量
Figure GDA00004288311000001612
的第μ行的回波数据，…，
Figure GDA00004288311000001613
为回波向量
Figure GDA00004288311000001614
的第1O行的回波数据，μ表示步骤5得到的空间直角坐标系坐标轴kz方向采样编号，取值范围是：μ=1，2，…，10；
回波向量
Figure GDA00004288311000001615
第μ行回波数据
Figure GDA00004288311000001616
其中
Figure GDA00004288311000001617
为回波向量
Figure GDA00004288311000001618
第1行的数据，
Figure GDA00004288311000001619
为回波向量
Figure GDA00004288311000001620
的第2行数据，…，为回波向量
Figure GDA00004288311000001622
第n行数据，…，
Figure GDA00004288311000001623
为回波向量
Figure GDA00004288311000001624
第20行的数据，n表示步骤5得到的空间直角坐标系坐标轴ky方向采样编号，取值范围是：n=1，2，…，20；
Figure GDA00004288311000001625
的第n行数据
Figure GDA00004288311000001626
其中
Figure GDA00004288311000001627
Figure GDA00004288311000001628
为步骤6.1中波数域回波数据
Figure GDA00004288311000001629
在坐标(kx，m，ky，n，kz，μ)处的数据，yμ，n，1为回波向量
Figure GDA00004288311000001630
在波数域回波数据
Figure GDA00004288311000001631
中坐标(kx，m，ky，n，kz1)处的数据，yμ，n，2为回波向量
Figure GDA00004288311000001632
在波数域回波数据中坐标(kx，m，ky，n，kz2)处的数据，…，yμ，n，m为回波向量
Figure GDA00004288311000001634
在波数域回波数据
Figure GDA00004288311000001635
中坐标(kx，m，ky，n，kz，μ)处的数据，…，yp，n，20为回波向量
Figure GDA00004288311000001636
在波数域回波数据
Figure GDA00004288311000001637
中坐标(kx，m，ky，n，kz，20)处的数据，m表示步骤5得到的空间直角坐标系坐标轴kx方向采样编号，取值范围是：m=1，2，…，20。
将yμ，n，m记做yq，则回波向量
Figure GDA0000428831100000171
可以表示为
Figure GDA0000428831100000172
其中y1为回波向量
Figure GDA0000428831100000173
第1行的回波数据，y2为回波向量第2行的回波数据，…，yq为回波向量
Figure GDA0000428831100000175
第q行的回波数据，…，y4000为回波向量
Figure GDA0000428831100000176
第4000行的回波数据，其中q为回波向量的行标号，取值范围是：q=1，2，…4000。yq在波数域对应的位置坐标记做 s ‾ q = ( k x , q , k y , q , k z , q ) .
步骤6.3
根据步骤2中构建的空间直角坐标系，将观测场景以场景中心O为中心点，在x轴方向以0.1m为间隔，将x轴方向的观测场景范围等分为20份；在y轴方向以0.1m为间隔，将y轴方向的观测场景范围等分为20份；在z轴方向以0.1m为间隔，将z轴方向的观测场景范围等分为10份；通过上述方法观测场景被划分为20×20×10个单元格；每个单元格在直角坐标系下的位置矢量记做
Figure GDA0000428831100000179
其中f=1，2，…20,g=1，2，…20,h=1，2，…10。
构造测量矩阵为测量矩阵
Figure GDA00004288311000001711
的第t行数据为
Figure GDA00004288311000001712
(20，2(，其中αt(f=，g，，hx)+tj，(kx，t，ky，t，kz，t)为步骤6.3中回波向量
Figure GDA00004288311000001713
第t行数据对应的波数域位置坐标，(xf，yg，zh)为观测场景的单元格位置坐标，αt(1，1，1)为f=1，g=2，h=1时的值，αt(1，1，2)为f=1，g=1，h=2时的值，…，αt(20，20,10)为f=20,g=20,h=10时的值，f表示空间直角坐标系坐标轴x轴方向采样编号，取值范围是：f=1，2，…20；g表示空间直角坐标系坐标轴y轴方向采样编号，取值范围是：g=1，2，…20；h表示空间直角坐标系坐标轴z轴方向采样编号，取值范围是：h=1，2，…10。
步骤6.4
随机抽取步骤6.2中回波向量
Figure GDA0000428831100000181
中的第N1行作为抽取向量y的第1行，抽取
Figure GDA0000428831100000182
中的第N2行作为抽取向量y的第2行，…，同理，抽取
Figure GDA0000428831100000183
中的第N800行作为抽取向量y的第800行，其中y为抽取后的回波向量。
同理，随机抽取步骤6.3中测量矩阵
Figure GDA0000428831100000184
的第N1行作为抽取测量矩阵A的第1行，抽取
Figure GDA0000428831100000185
中的第N2行作为抽取测量矩阵A的第2行，…，抽取
Figure GDA0000428831100000186
中的第N800行作为抽取测量矩阵A的第800行，A为抽取后的测量矩阵。
步骤7、特征增强图像形成方法迭代求解
第1次求解，将步骤6中抽取后的回波向量y和抽取后的测量矩阵A带入公式(1)中求解：
Figure GDA0000428831100000187
公式(1)中：λ为步骤1中初始化的拉格朗日乘数， Λ ( f ( 0 ) ) = diag { 1 ( | f 1 ( 0 ) | + 0.1 ) 2 , 1 ( | f 2 ( 0 ) | + 0.1 ) 2 , . . . , 1 ( | f 4000 ( 0 ) | + 0.1 ) 2 } , f1 (0)为步骤1中初始解向量f(0)的第1行数据，f2 (0)为步骤1中初始解向量f(0)的第2行数据，为步骤1中初始解向量f(0)的第4000行数据，Λ(f(0))为第1次求解使用的加权L1范数的梯度，f(1)为第1次求解得出的解向量。
判断
Figure GDA0000428831100000189
是否成立，若成立则求解结束；若不成立，则判断求解次数是否达到第1000次，若本次求解为第1000求解则跳出求解，否则继续进行第2次求解。
第2次求解，将步骤6中抽取后的回波向量y和抽取后的测量矩阵A带入公式(2)中求解：
f(2)=(2AHA+λΛ(f(1)))-12AHy                  (2)
公式(2)中：λ为步骤1中初始化的拉格朗日乘数， Λ ( f ( 1 ) ) = diag { 1 ( | f 1 ( 1 ) | + 0.1 ) 2 , 1 ( | f 2 ( 1 ) | + 0.1 ) 2 , . . . , 1 ( | f 4000 ( 1 ) | + 0.1 ) 2 } , f1 (1)为第1次求解得出的解向量f(1)的第1行数据，f2 (1)为第1次求解得出的解向量f(1)的第2行数据，f4000 (1)为第1次求解得出的解向量f(1)的第4000行数据，Λ(f(1))为第2次求解使用的加权L1范数的梯度，f(2)为第2次求解得出的解向量。
判断是否成立，若成立则结束求解；若不成立，判断求解次数是否达到第1000次，若本次求解为第1000求解则跳出求解，否则进行第3次求解。
……
同理，第J次求解，将步骤6中抽取后的回波向量y和抽取后的测量矩阵A带入公式(J)中求解：
f(J)=(2AHA+λΛ(f(J-1)))-12AHy             (J)
公式(J)中：λ为步骤1中初始化的拉格朗日乘数， Λ ( f ( J ) ) = diag { 1 ( | f 1 ( J - 1 ) | + 0.1 ) 2 , 1 ( | f 2 ( J - 1 ) | + 0.1 ) 2 , . . . , 1 ( | f 4000 ( J - 1 ) | + 0.1 ) 2 } , f1 (J-1)为第J-1次求解得出的解向量f(J-1)的第1行数据，f2 (J-1)为第J-1次求解得出的解向量f(J-1)的第2行数据，
Figure GDA0000428831100000195
为第J-1次求解得出的解向量f(J-1)的第4000行数据，Λ(f(J-1))为第J次求解使用的加权L1范数的梯度，f(J)为第J次求解得出的解向量，J表示求解编号，且满足1≤J≤1000。
判断
Figure GDA0000428831100000194
是否成立，若成立则结束求解；若不成立，判断求解次数是否达到第1000次，若本次求解为第1000求解则跳出求解，否则进行第J+1次求解。
经过上述操作，即可得到基于特征增强图像形成方法圆周三维成像合成孔径雷达图像数据。
通过本发明具体实施方式的仿真及测试，本发明所提出的基于特征增强图像形成方法圆周三维成像合成孔径雷达成像方法，与现有的三维合成孔径雷达成像图像相比，它实现降低旁瓣和主瓣宽度从而提高分辨率的目的。

Claims (1)

1.一种多航过圆周SAR三维成像方法，其特征包括以下步骤：
步骤1、初始化圆周三维成像合成孔径雷达成像系统参数以及特征增强图像形成方法参数：
初始化成像系统参数包括：雷达发射电磁波的载频，记做fC；雷达发射电磁波的带宽，记做BW；雷达运动平台的飞行半径，记做R；雷达运动平台的方位角，记做θ；雷达运动平台的俯仰角，记做
Figure FDA0000428831090000012
；电磁波的传播速度，记做C；雷达在方位维的观测范围，记做θw；合成孔径雷达距离向回波数据，记做
Figure FDA0000428831090000011
雷达快时间采样点数，记做K；雷达慢时间采样点数，记做NP；多航过雷达飞行模式航过次数，记做L；观测场景在距离向的观测范围，记做Wx；观测场景在方位向的观测范围，记做Wy；观测场景在高度向的观测范围，记做Wz；回波数据在波数域插值界限，记做B；观测场景像素划分间隔，记做Δw；观测场景的中心点，记做O；其中，雷达发射电磁波的载频，雷达发射电磁波的带宽，电磁波的传播速度在雷达系统的设计过程中均已确定，雷达运动平台的飞行半径，雷达运动平台的方位角，雷达运动平台的俯仰角，合成孔径雷达在方位维的观测范围，观测场景在距离向、方位向和高度向的观测范围，多航过雷达飞行模式航过次数，回波数据在波数域插值界限，观测场景像素划分间隔在预先确定观测方案中为均为已知；
初始化的特征增强图像形成方法参数包括：拉格朗日乘数，记做λ；特征增强图像形成算法误差门限，记做Δ；共轭梯度算法求解误差门限，记做ε；共轭梯度算法最大迭代次数，记做cgnum；特征增强图像形成方法最大迭代次数，记做itenum；初始解向量，记做f(0)；L1范数逼近常数，记做η；其中，拉格朗日常数特征增强图像形成算法误差门限Δ，共轭梯度算法求解误差门限ε，共轭梯度算法最大迭代次数cgnum，特征增强图像形成方法最大迭代次数itenum，L1范数逼近常数η，初始解向量f(0)均为已知；
步骤2、构建多航过圆周三维成像合成孔径雷达的运动模型
以场景中心点O为原点，以地面为xy平面，以高度向为z轴，构建空间直角坐标系，雷达运动平台以[0，0，z1]，[0，0，z2]，…，[0，0，zτ]，…，[0，0，zL]为圆心做半径为R的圆周运动，其中zτ为雷达运动平台飞行高度，合成孔径雷达运动平台不同高度平面编号为τ，取值范围是：τ=1，2，…，L，其中z1为第1次时合成孔径雷达雷达运动平台飞行高度，z2为第2次合成孔径雷达运动平台飞行高度，…，zτ为第τ次合成孔径雷达运动平台飞行高度，…，zL为第L次合成孔径雷达运动平台飞行高度；
步骤3、圆周三维成像合成孔径雷达原始数据进行距离压缩
采用合成孔径雷达标准距离压缩方法对步骤1中提供合成孔径雷达距离向回波数据进行压缩，得到距离压缩后的圆周三维成像合成孔径雷达数据，记做距离压缩后的圆周三维成像合成孔径雷达数据是一个三维数组，包括距离维、沿航迹维、高度维；
步骤4、获取圆周三维成像合成孔径雷达波数域的回波数据
对步骤3得到的距离压缩后的圆周三维成像合成孔径雷达数据
Figure FDA0000428831090000024
中的每个距离维的数据进行一维快速傅里叶变换(FFT)，得到圆周三维成像合成孔径雷达波数域的回波数据，该数据是一个三维数组，记做
Figure FDA0000428831090000025
表示圆周三维成像合成孔径雷达波数域的回波数据，包括波数维、方位角维和俯仰角维，ρk表示合成孔径雷达回波数据在波数域中波数数据，θi表示合成孔径雷达回波数据在波数域中方位角数据，
Figure FDA0000428831090000026
表示合成孔径雷达回波数据在波数域中俯仰角数据，其中k表示波数域中波数方向采样编号，取值范围是：k=1，2，…，K，i表示波数域中方位角方向采样编号，取值范围是：i=1，2，…，NP，l表示波数域中俯仰角方向采样编号，取值范围是：l=1，2，…，L；
步骤5、通过插值获得直角坐标系下波数域回波数据，包括步骤5.1、步骤5.2：
步骤5.1
取步骤4中得到的圆周三维成像合成孔径雷达波数域的回波数据
Figure FDA0000428831090000027
的中心点，记做其中ρ0=(ρ1+ρK)／2，
Figure FDA00004288310900000213
Figure FDA0000428831090000029
为步骤4中k=1时的波数，ρk为步骤4中k=K时的波数，θ1为步骤4中i=1时的方位角，
Figure FDA00004288310900000210
为步骤4中i=NP时的方位角，
Figure FDA00004288310900000211
为步骤4中l=1时的俯仰角，
Figure FDA00004288310900000212
为步骤4中l=L时的俯仰角；
采用坐标系变换方法进行坐标系变换：利用直角坐标系和极坐标系的对应关系构建三维空间波数域的直角坐标系，得到空间直角坐标系；其中的坐标轴记做kx、ky和kz；以
Figure FDA0000428831090000031
为中心点，在kx轴方向以2π／Wx为间隔，在ky轴方向以2π／Wy为间隔，在kz轴方向以2π／Wz为间隔，在空间直角坐标系下将三维空间波数域离散为M×N×P个单元格，其中，M为kx轴方向单元格数目，N为ky轴方向单元格数目，P为kz轴方向单元格数目，其中M，N，P均为自然数，并满足 M ≤ BW x 2 π , N ≤ BW y 2 π , P ≤ BW z 2 π ;
步骤5.2
利用三维线性插值方法，将步骤4得到的回波数据
Figure FDA0000428831090000033
插值到步骤5.1得到的直角坐标系下M×N×P个单元格内，得到直角坐标系下的波数域回波数据，记做
Figure FDA00004288310900000323
步骤6、构造回波数据的线性测量矩阵和散射场景的线性观测模型，包括步骤6.1、步骤6.2、步骤6.3、步骤6.4；
步骤6.1
将步骤5得到的直角坐标系下的波数域回波数据
Figure FDA0000428831090000034
用表示直角坐标系下波数域回波数据在坐标(kx，m，ky，n，kz，μ)处的数据，其中，m表示步骤5得到的空间直角坐标系坐标轴kx方向采样编号，取值范围是：m=1，2…，M，n表示步骤5得到的空间直角坐标系坐标轴ky方向采样编号，取值范围是：n=1，2…，N，μ表示步骤5得到的空间直角坐标系坐标轴kz方向采样编号，取值范围是：μ=1，2…，P；
步骤6.2
将步骤5得到的直角坐标系下的波数域回波数据
Figure FDA0000428831090000036
定义为一个回波向量，记做其中
Figure FDA0000428831090000038
为回波向量
Figure FDA00004288310900000324
的第1行的回波数据，
Figure FDA0000428831090000039
为回波向量
Figure FDA00004288310900000310
的第2行的回波数据，…，
Figure FDA00004288310900000311
为回波向量
Figure FDA00004288310900000325
的第μ行的回波数据，…，
Figure FDA00004288310900000312
为回波向量
Figure FDA00004288310900000313
的第P行的回波数据；
回波向量
Figure FDA00004288310900000314
第μ行回波数据 y ‾ μ = [ y ‾ μ , 1 , y ‾ μ , 2 , . . . , y ‾ μ , n , . . . , y ‾ μ , N ] T , 其中
Figure FDA00004288310900000316
为回波向量
Figure FDA00004288310900000317
第1行的数据，为回波向量的第2行数据，…，为回波向量
Figure FDA00004288310900000320
第n行数据，…，
Figure FDA00004288310900000321
为回波向量
Figure FDA00004288310900000322
第N行的数据；
Figure FDA0000428831090000041
的第n行数据 y ‾ μ , n = [ y μ , n , 1 , y μ , n , 2 , . . . , y μ , n , m , . . . , y μ , N , M ] T , 其中 yμ，n，1为回波向量在波数域回波数据
Figure FDA0000428831090000046
中坐标(kx，1，ky，nkz，u)处的数据，yμ，n，2为回波向量
Figure FDA0000428831090000047
在波数域回波数据中坐标(kx，2，ky，nkz，u)处的数据，…，yμ，n，m为回波向量
Figure FDA0000428831090000049
在波数域回波数据
Figure FDA00004288310900000410
中坐标(kx，m，ky，nkz，u)处的数据，…，yμ，N，M为回波向量
Figure FDA00004288310900000411
在波数域回波数据
Figure FDA00004288310900000412
中坐标(kx，M，ky，nkz，u)处的数据；
将yμ，n，m记做yq，则回波向量
Figure FDA00004288310900000413
表示为其中y1为回波向量
Figure FDA00004288310900000415
第1行的回波数据，y2为回波向量
Figure FDA00004288310900000416
第2行的回波数据，…，yq为回波向里
Figure FDA00004288310900000417
第q行的回波数据，…，yQ为回波向量
Figure FDA00004288310900000418
第Q行的回波数据，其中q为回波向量的行标号，取值范围是：q=1,2，…，Q，Q=M×N×P；yq在波数域对应的位置坐标记做
Figure FDA00004288310900000420
步骤6.3
根据步骤2中构建的空间直角坐标系，将观测场景以场景中心O为中心点，在x轴方向以Δw为间隔，将x轴方向的观测场景范围Wx等分为X份；在y轴方向以Δw为间隔，将y轴方向的观测场景范围Wy等分为Y份；在z轴方向以Δw为间隔，将z轴方向的观测场景范围Wz等分为Z份；通过上述方法观测场景被划分为X×Y×Z个单元格，其中X为距离向单元格数目，Y为方位向单元格数目，Z为高度向单元格数目，且满足
Figure FDA00004288310900000421
每个单元格在直角坐标系下的位置矢量记做
Figure FDA00004288310900000422
其中，f表示空间直角坐标系坐标轴x轴方向采样编号，取值范围是：f=1,2，…X；g表示空间直角坐标系坐标轴y轴方向采样编号，取值范围是：g=1,2，…，Y；h表示空间直角坐标系坐标轴z轴方向采样编号，取值范围是：h=1,2，…，Z；
构造测量矩阵为
Figure FDA00004288310900000423
其中，
Figure FDA00004288310900000424
为测量矩阵
Figure FDA00004288310900000425
的第1行数据，
Figure FDA00004288310900000426
为测量矩阵
Figure FDA00004288310900000427
的第2行数据，…，
Figure FDA00004288310900000428
为测量矩阵的第t行数据，…，
Figure FDA00004288310900000430
为测量矩阵
Figure FDA00004288310900000431
的第Q行数据，其中t为测量矩阵
Figure FDA00004288310900000432
的行标号，t=1,2，…，Q，Q=M×N×P；
测量矩阵
Figure FDA0000428831090000051
的第t行数据为 α ‾ t = [ α t ( 1,1,1 ) , α t ( 1,1,2 ) , . . . , α t ( f , g , h ) , . . . , α t ( X , Y , Z ) ] , 其中αt(1，1，1)为测量矩阵
Figure FDA0000428831090000053
的第t行第一列数据，αt(1，1，2)为为测量矩阵
Figure FDA0000428831090000054
的第t行第二列数据，αt(f，g，h)=exp{j(kx，fxf+ky，tyg+kz，tzh)}，(kx，f，ky，t，kz，t)为步骤6.3中回波向量
Figure FDA0000428831090000055
第t行数据对应的波数域位置坐标，(xf，yg，zh)为观测场景的单元格位置坐标，αt(X，Y，Z)为f=X，g=Y，h=Z时的值；
步骤6.4
随机抽取步骤6.2中回波向量
Figure FDA0000428831090000056
中的第N1行作为抽取后的回波向量y的第1行，抽取回波向量
Figure FDA0000428831090000057
中的第N2行作为抽取后的回波向量y的第2行，…，抽取回波向量
Figure FDA0000428831090000058
中的第Nr行作为抽取后的回波向量y的第r行，…，同理，抽取回波向量
Figure FDA0000428831090000059
中的第Nnum行作为抽取后的回波向量y的第num行，N1为第1次抽取回波向量
Figure FDA00004288310900000510
的行标，N2为第2次抽取回波向量
Figure FDA00004288310900000511
的行标，…，Nr为第r次抽取回波向量
Figure FDA00004288310900000512
的行标，…，Nnum为第num次抽取回波向量的行标，其中r为抽取后的回波向量y的行标，取值范围是：r=1，2，…，num，num为抽取行的总数，并且N1，N2，…，Nnum，num均为不大于Q的自然数；
同理，抽取步骤6.3中测量矩阵
Figure FDA00004288310900000514
中的第N1行作为抽取后的测量矩阵A的第1行，抽取测量矩阵
Figure FDA00004288310900000515
中的第N2行作为抽取后的测量矩阵A的第2行，…，抽取测量矩阵
Figure FDA00004288310900000516
中的第Nr行作为抽取后的测量矩阵A的第r行，…，同理，抽取测量矩阵
Figure FDA00004288310900000517
中的第Nnum行作为抽取后的测量矩阵A的第num行，N1为第1次抽取测量矩阵
Figure FDA00004288310900000518
的行标，N2为第2次抽取测量矩阵
Figure FDA00004288310900000519
的行标，…，Nr为第r次抽取测量矩阵
Figure FDA00004288310900000520
的行标，…，Nnum为第num次抽取测量矩阵
Figure FDA00004288310900000521
的行标；
步骤7、改进特征增强图像形成方法迭代求解
第1次求解，将步骤6中抽取后的回波向量y和抽取后的测量矩阵A带入公式(1)中求解：
f(1)=(2AHA+λΛ(f(0)))-12AHy                   (1)
公式(1)中： Λ ( f ( 0 ) ) = diag { 1 ( | f 1 ( 0 ) + η | ) 2 , 1 ( | f 2 ( 0 ) + η | ) 2 , . . . , 1 ( | f Q ( 0 ) + η | ) 2 } , f1 (0)为步骤1中初始解向量f(0)的第1行数据，f2 (0)为步骤1中初始解向量f(0)的第2行数据，fQ (0)为步骤1中初始解向量f(0)的第Q行数据，Λ(f(0))为第1次求解使用的加权L1范数的梯度，f(1)为第1次求解得出的解向量；
判断
Figure FDA0000428831090000061
是否成立，若成立则求解结束；若不成立，则判断求解次数是否达到第itenum次，若本次求解为第itenum求解则求解结束，否则继续进行第2次求解；
第2次求解，将步骤6中抽取后的回波向量y和抽取后的测量矩阵A带入公式(2)中求解：
f(2)=(2AHA+λΛ(f(1)))-12AHy                   (2)
公式(2)中： Λ ( f ( 1 ) ) = diag { 1 ( | f 1 ( 1 ) + η | ) 2 , 1 ( | f 2 ( 1 ) + η | ) 2 , . . . , 1 ( | f Q ( 1 ) + η | ) 2 } , f1 (1)为第1次求解得出的解向量f(1)的第1行数据，f2 (1)为第1次求解得出的解向量f(1)的第2行数据，fQ (1)为第1次求解得出的解向量f(1)的第Q行数据，Λ(f(1))为第2次求解使用的加权L1范数的梯度，f(2)为第2次求解得出的解向量；
判断
Figure FDA0000428831090000064
是否成立，若成立则结束求解；若不成立，判断求解次数是否达到第itenum次，若本次求解为第itenum求解则求解结束，否则进行第3次求解；
同理，第J次求解，将步骤6中抽取后的回波向量y和抽取后的测量矩阵A带入公式(J)中求解：
f(J)=(2AHA+λΛ(f(J-1)))-12AHy       (J)
公式(J)中： Λ ( f ( J - 1 ) ) = diag { 1 ( | f 1 ( J - 1 ) + η | ) 2 , 1 ( | f 2 ( J - 1 ) + η | ) 2 , . . . , 1 ( | f Q ( J - 1 ) + η | ) 2 } , f1 (J-1)为第J-1次求解得出的解向量f(J-1)的第1行数据，f2 (J-1)为第J-1次求解得出的解向量f(J-1)的第2行数据，fQ (J-1)为第J-1次求解得出的解向量f(J-1)的第Q行数据，Λ(f(J-1))为第J次求解使用的加权L1范数的梯度，f(J)为第J次求解得出的解向量，J表示求解编号，且满足1≤J≤itenum；
判断
Figure FDA0000428831090000071
是否成立，若成立则结束求解；若不成立，判断求解次数是否达到第itenum次，若本次求解为第itenum求解则求解结束，否则进行第J+1次求解；
经过上述步骤，得到圆周三维成像合成孔径雷达图像数据。
CN201210333250.8A 2012-09-11 2012-09-11 一种多航过圆周sar三维成像方法 Expired - Fee Related CN102914773B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201210333250.8A CN102914773B (zh) 2012-09-11 2012-09-11 一种多航过圆周sar三维成像方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201210333250.8A CN102914773B (zh) 2012-09-11 2012-09-11 一种多航过圆周sar三维成像方法

Publications (2)

Publication Number Publication Date
CN102914773A CN102914773A (zh) 2013-02-06
CN102914773B true CN102914773B (zh) 2014-04-09

Family

ID=47613218

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201210333250.8A Expired - Fee Related CN102914773B (zh) 2012-09-11 2012-09-11 一种多航过圆周sar三维成像方法

Country Status (1)

Country Link
CN (1) CN102914773B (zh)

Families Citing this family (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103675813B (zh) * 2013-09-22 2017-07-25 中国科学院电子学研究所 一种基于标志点相位梯度提取的圆迹sar轨迹重建方法
CN109683148A (zh) * 2017-10-19 2019-04-26 深圳市新益技术有限公司 目标散射特性测试方法
CN109471073B (zh) * 2018-10-31 2020-08-28 中国科学院电子学研究所 基于增广拉格朗日粒子群算法的nlfm信号生成方法及装置
CN111220979B (zh) * 2020-01-16 2022-05-13 电子科技大学 一种用于曲线合成孔径雷达的成像方法
CN112379371B (zh) * 2020-10-29 2022-09-09 中国科学技术大学 基于优化理论的无线电信号三维成像方法及系统
CN112330560B (zh) * 2020-11-05 2024-02-20 中国科学院国家空间科学中心 一种合成孔径雷达数据图像可视化增强方法及系统
CN113050087B (zh) * 2021-03-17 2022-08-05 电子科技大学 一种用于曲线合成孔径雷达的子孔径划分方法
CN113391309B (zh) * 2021-06-15 2022-09-09 电子科技大学 一种火星探测器雷达径向下视成像方法
CN115657032B (zh) * 2022-12-27 2023-03-10 中国人民解放军国防科技大学 一种极化csar车辆目标三维重建方法及装置
CN117289277B (zh) * 2023-11-27 2024-01-30 中山大学 一种基于子带分割合成的多频雷达三维成像方法及系统

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2001083243A (ja) * 1999-09-13 2001-03-30 Mitsubishi Electric Corp 干渉型合成開口レーダによる地形の３次元情報抽出装置
CN101900812B (zh) * 2009-05-25 2012-11-14 中国科学院电子学研究所 一种圆迹合成孔径雷达的大场景极坐标格式三维成像方法
CN102313888A (zh) * 2010-06-29 2012-01-11 电子科技大学 一种基于压缩传感的线阵sar三维成像方法

Also Published As

Publication number Publication date
CN102914773A (zh) 2013-02-06

Similar Documents

Publication Publication Date Title
CN102914773B (zh) 一种多航过圆周sar三维成像方法
CN102323583B (zh) 一种超分辨线阵三维合成孔径雷达成像方法
US8193967B2 (en) Method and system for forming very low noise imagery using pixel classification
US11506776B2 (en) Method and device with improved radar resolution
CN103941257B (zh) 一种基于波数能量谱的导航雷达图像反演海面风向的方法
CN104635230B (zh) 一种用于mimo‑sar近场测量成像方位向旁瓣抑制的方法
CN104898119B (zh) 一种基于相关函数的动目标参数估计方法
CN110568434B (zh) 一种多通道匀加速sar动目标二维速度估计方法
CN102313888A (zh) 一种基于压缩传感的线阵sar三维成像方法
CN103869311A (zh) 实波束扫描雷达超分辨成像方法
CN104950305A (zh) 一种基于稀疏约束的实波束扫描雷达角超分辨成像方法
CN103163523A (zh) 基于压缩感知的低空风切变风速估计方法
CN112415515B (zh) 一种机载圆迹sar对不同高度目标分离的方法
CN113534065B (zh) 一种雷达目标微动特征提取与智能分类方法及系统
Setsu et al. Super-Resolution Doppler Velocity Estimation by Kernel-Based Range–$\tau $ Point Conversions for UWB Short-Range Radars
CN104280566A (zh) 基于空时幅相估计的低空风切变风速估计方法
CN104166129A (zh) 一种实波束雷达迭代最小均方误差角超分辨方法
CN103630903B (zh) 基于顺轨干涉sar测量海面流场径向速度的方法
CN103630899B (zh) 地面运动目标高分辨雷达压缩感知成像的方法
CN105447867A (zh) 基于isar图像的空间目标姿态估计方法
CN103605118A (zh) 一种利用极地探冰雷达提取极地冰层位方法
CN106646466A (zh) 一种基于主成分分析的加权后向投影算法的成像方法
Garry et al. Passive ISAR part I: framework and considerations
CN103954962B (zh) 一种基于压缩感知的isar成像脉冲估计算法
CN110967677B (zh) 时频域分级多目标分辨与测量方法

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20140409