CN109061642B - 一种贝叶斯迭代重加权稀疏自聚焦阵列sar成像方法 - Google Patents

一种贝叶斯迭代重加权稀疏自聚焦阵列sar成像方法 Download PDF

Info

Publication number
CN109061642B
CN109061642B CN201810767253.XA CN201810767253A CN109061642B CN 109061642 B CN109061642 B CN 109061642B CN 201810767253 A CN201810767253 A CN 201810767253A CN 109061642 B CN109061642 B CN 109061642B
Authority
CN
China
Prior art keywords
array
equidistant
initialization
vector
total number
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.)
Active
Application number
CN201810767253.XA
Other languages
English (en)
Other versions
CN109061642A (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 CN201810767253.XA priority Critical patent/CN109061642B/zh
Publication of CN109061642A publication Critical patent/CN109061642A/zh
Application granted granted Critical
Publication of CN109061642B publication Critical patent/CN109061642B/zh
Active 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
    • 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
    • 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
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • G01S7/41—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section

Abstract

本发明公开了一种贝叶斯迭代重加权稀疏自聚焦阵列SAR成像方法，它是针对阵列SAR回波信号中存在的相位误差对成像结果的影响，基于传统的贝叶斯迭代最小化自聚焦稀疏成像(SAFBRIM)算法的基础上，通过建立阵列SAR原始回波信号与观测场景目标空间中散射系数的线性测量矩阵，对算法中代价函数中的范数项进行迭代自适应重加权处理，对距离向进行脉冲压缩、划分等距离面，然后再对每一个等距离的二维平面进行估计。本发明对每一个范数项赋予了不同的加权系数，然后对图像进行重构，能获得更高质量的阵列SAR成像结果。本发明具有重构精度高、有效降低相位误差的优势，可适用于阵列合成孔径雷达成像等领域。

Description

一种贝叶斯迭代重加权稀疏自聚焦阵列SAR成像方法
技术领域
本发明属于雷达技术领域，它特别涉及了合成孔径雷达(SAR)成像技术领域。
背景技术
作为一种工作在微波波段的有源雷达，合成孔径雷达(Synthetic ApertureRadar,SAR)具有全天时、全天候的成像能力，即无论是白天或黑夜、晴天还是雷雨风雪天气，都可以随时随地成像，克服了光学和红外系统不能在晚上和复杂天气条件进行成像的缺点。传统的SAR成像一般只具有二维成像分辨率，在一些起伏比较大的地方比如陡峭的山峰、峡谷以及城市中矗立挺拔的高楼时，传统SAR成像存在的失真(阴影遮挡效应、空间模糊、顶底倒置等)导致空间的一些重要信息(比如高度)丢失，所以能对目标进行三维成像是非常有必要的，为了适应这种需求，目前常见的三维成像技术有圆周SAR(Circular SAR)三维成像、层析SAR(Tomography SAR)三维成像、阵列SAR(Array SAR，ASAR)三维成像。
阵列SAR三维成像的基本原理是在切行迹向添加阵列天线，通过沿航迹向平台的飞行形成虚拟的面阵进而获得二维分辨率，距离向再通过脉冲压缩技术获得第三维的分辨率。相比于圆周SAR三维成像，阵列SAR三维成像不需要圆周运动的轨迹；相比于层析SAR三维成像需要航过多次，阵列SAR三维成像只需一次航过，所以阵列SAR三维成像相对于层析SAR和圆周SAR三维成像有更强的灵活性。目前阵列SAR三维成像技术在地形测绘、城市测绘、灾难救援、军事探测等领域发挥着重要的作用。
传统基于匹配滤波的SAR成像方法的分辨率受到限制，具体来说就是距离向的分辨率受信号带宽的影响，沿航迹分辨率受合成孔径长度的影响，切航迹的分辨率受阵列天线的影响。尤其是切航迹的分辨率，如果按照传统的方法很难提高。如果一个信号是稀疏的或者是可压缩的，那么这个信号就能以低于Nyquist采样定理要求的采样率精确的重构出该信号，这就是压缩感知(Compressed Sensing，CS)的基本思想。针对压缩感知理论用于SAR成像，目前的重构算法大概可以分为以下几类：贪婪追踪算法、凸松弛算法、贝叶斯框架算法、组合算法。
在实际阵列三维SAR回波数据获取过程中，由于系统外部和内部测量的不确定性，如平台受气流扰动和GPS/IMU导航定位精度的影响，即使使用代价昂贵的测量设备，运动平台与天线的位置测量仍不可避免存在误差，而测量参数精度值往往达不到稀疏成像的要求。为了实现高分辨阵列三维SAR稀疏成像，除了利用外部测量参数进行补偿，还必须校正测量数据中残余的相位误差。然而，目前大多数CS稀疏成像算法很少考虑测量矩阵中的误差或不确定影响。因此，研究阵列三维SAR回波信号的相位误差表示模型，实现高精度稀疏自聚焦成像，是阵列三维SAR稀疏成像中的迫切问题。
为了降低回波信号中相位误差对高精度稀疏成像的影响，本发明提出了一种基于贝叶斯迭代自适应重加权范数最小化的稀疏自聚焦(An Iterative Adaptive ReweightedNorm Minimization Sparsity Autofocus Algorithm via Bayesian Recovery，IARNSABR)阵列SAR成像算法。
发明内容
为了提高阵列SAR成像质量，降低回波信号中存在的相位误差对阵列SAR成像的影响，本发明提出的基于贝叶斯迭代自适应重加权范数最小化的一种稀疏自聚焦阵列SAR成像算法，该方法先对距离向进行脉冲压缩、划分等距离面，对代价函数中的范数项进行自适应重加权处理后，利用最大似然准则进行估计，该算法对每一个范数项赋予了不同的加权系数，有效降低了相位误差对于高精度稀疏成像的影响，能获得更高质量的阵列SAR成像结果。
为了方便描述本发明的内容，首先作以下术语定义：
定义1、合成孔径雷达(SAR)
合成孔径雷达是将雷达固定于载荷运动平台上，结合平台的运动以合成等效阵列以实现阵列向的分辨率，再利用雷达波束向回波延时实现距离一维成像，从而实现对观测目标二维成像的一种合成孔径雷达技术。
定义2、标准合成孔径雷达回波数据距离向脉冲压缩
标准合成孔径雷达回波数据距离向脉冲压缩是指利用合成孔径雷达发射信号参数，采用匹配滤波技术对合成孔径雷达的距离向信号进行信号聚焦成像的过程。详见文献“雷达成像技术”，保铮，邢孟道，王彤，电子工业出版社，2005。
定义3、范数
设X是复数域
Figure GDA0003389215180000023
上线性空间，其中
Figure GDA0003389215180000024
表示复数域，若它满足如下性质：||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范数表达式为
Figure GDA0003389215180000021
其中xi为向量X的第i个元素，∑|·|表示绝对值求和运算符号，向量X的L1范数表达式为
Figure GDA0003389215180000022
向量X的L2范数表达式为
Figure GDA0003389215180000031
向量X的L0范数表达式为
Figure GDA0003389215180000032
且xi≠0。详见文献“矩阵理论”，黄廷祝等编著，高等教育出版社出版。
定义4、方位向、距离向
将雷达平台运动的方向叫做方位向，将垂直于方位向的方向叫做距离向。
定义5、压缩感知稀疏重构理论
如果一个信号是稀疏的或可压缩的，那么该信号就可以用远低于奈奎斯特采样定理所要求的采样率来无失真的重构出该信号。如果信号稀疏，并且测量矩阵满足不相干和RIP属性，使用压缩感知恢复的信号稀疏重建可以通过解决以下最优化问题来实现：
Figure GDA0003389215180000033
其中，α是估计信号，y是测量信号，Θ是测量矩阵，ε是噪声门限。详见文献“阵列三维合成孔径雷达稀疏成像技术研究”韦顺军，2013。
定义6、基于贝叶斯迭代最小化自聚焦稀疏成像(SAFBRIM)算法
基于贝叶斯迭代最小化自聚焦稀疏成像算法(Sparse Autofocus BayesianRecovery via Iterative Minimum)由电子科技大学的韦顺军副教授于2013年提出,详见文献“阵列三维合成孔径雷达稀疏成像技术研究”韦顺军，2013
定义7、合成孔径雷达原始回波仿真方法
合成孔径雷达原始回波仿真方法是指基于合成孔径雷达成像原理仿真出一定系统参数条件下具有合成孔径雷达回波信号特性的原始信号的方法，详见文献“张朋，合成孔径雷达回波信号仿真研究，西北工业大学博士论文，2004”。
定义8、阵列SAR的快时刻和慢时刻
阵列SAR运动平台飞过一个方位向合成孔径长度所需要的时间称为慢时间，雷达系统以一定时间长度的重复周期发射接收脉冲，因此慢时间可以表示为一个以脉冲重复周期为步长的离散化时间变量，其中每一个脉冲重复周期离散时间变量值为一个慢时刻。快时刻是指在一个脉冲重复周期内，距离向采样回波信号的时间间隔变量。详见文献“合成孔径雷达成像原理”，皮一鸣等编著，电子科技大学出版社出版。
定义9、信号线性测量模型
对于一个数字信号测量系统，假设N×1维离散信号向量X＝[x1,x2,…,xN]T为该测量系统需要测量的信号，向量Y＝[y1,y2,…,yM]T为该测量系统输出的M维离散信号向量，其中T为转置运算符号，y1为向量Y中的第一个元素，y2表示向量Y中的第二个元素，yM表示向量Y中的第M个元素，信号的线性测量模型是指测量信号Y和被测量信号X的关系可以表示为Y＝AX，其中A为M×N矩阵，矩阵A为线性测量模型中信号X的测量矩阵。详见文献“阵列三维合成孔径雷达稀疏成像技术研究，韦顺军，2013”。
定义10传统理论成像分辨率
阵列SAR成像传统理论分辨率是指利用经典匹配滤波理论成像算法得到阵列SAR系统在距离向、方位向和切航迹向的成像分辨率。对于收发共用天线，阵列SAR距离向的分辨率记为ρr，近似表达式为
Figure GDA0003389215180000041
其中C为电磁波在空气中传播的速度，Br为阵列SAR发射信号带宽；方位向分辨率记为ρa，近似表达为
Figure GDA0003389215180000042
其中Da为天线在方位向的真实孔径；切航迹向的分辨率记为
Figure GDA0003389215180000043
其中λ为阵列SAR雷达载波波长，R0为阵列SAR平台到观测场景中心的参考斜距，L为阵列天线长度。
本发明提出的基于贝叶斯迭代自适应重加权范数最小化的一种稀疏自聚焦阵列SAR成像方法，它包括如下步骤：
步骤1、初始化SAR系统参数：
初始化SAR系统参数包括：平台速度矢量记为
Figure GDA0003389215180000044
阵列天线各阵元初始位置矢量，记做
Figure GDA0003389215180000045
其中n为天线各阵元序号，N为阵列天线的阵元总数；阵列天线长度，记做L；雷达发射信号载频为fc；雷达发射信号的调频斜率为fdr；脉冲重复时间记为PRI；雷达系统的脉冲重复频率为PRF；雷达发射信号带宽记做Br；电磁波在空气中的传播速度记做C；距离向快时间记做t，t＝1,2...T，T为距离向快时刻总数，方位向慢时刻记做l，l＝1,2,...K，K为方位向慢时刻总数；上述参数均为SAR系统标准参数，其中雷达信号载频fc，雷达发射信号的调频斜率fdr，脉冲重复时间PRI，雷达系统的脉冲重复频率PRF，雷达发射信号带宽Br，阵列天线的阵元总数N，阵列天线长度L在阵列SAR系统设计过程中已经确定；平台速度矢量记为
Figure GDA0003389215180000051
阵列天线各阵元初始位置矢量
Figure GDA0003389215180000052
在SAR观测方案设计中已经确定；根据SAR成像系统方案和观测方案，SAR成像方法需要的初始化成像系统参数均为已知；
步骤2、初始化SAR的观测场景目标空间参数：
初始化SAR的观测场景目标空间参数包括：以雷达波束照射场区域地平面和垂直于该地平面向上的单位向量所构成的空间直角坐标系作为阵列SAR的观测场景目标空间Ω；将观测场景目标空间Ω均匀划分为大小相等的立体单元格，称为分辨单元，单元网格在水平横向、水平纵向和高度向边长分别记为dx，dy和dz，观测场景空间在水平横向、水平纵向和高度向单元格数分别为Mx,My和Mz，单元格大小为阵列SAR系统传统理论成像分辨率；水平横向和水平纵向构成阵列维成像空间，在阵列平面维成像空间上第t个等距离单元格第m个元素的位置，记做
Figure GDA0003389215180000053
其中m＝(my-1)Mx+mx＝1,…,M，mx＝1,...,Mx，my＝1,...,My，t＝1,...,T；将观测场景目标空间中第t个等距离单元格第m个元素的散射系数记为
Figure GDA0003389215180000054
根据公式
Figure GDA0003389215180000055
计算得到散射系数矩阵，记做δ，散射系数矩阵δ由M行T列组成，其中T为步骤1中初始化得到的距离向快时刻总数，其中M＝Mx·My为阵列平面维成像空间的第t个等距离单元格阵列向单元格总数；初始化SAR的观测场景目标空间参数在SAR成像方案设计中已经确定
步骤3、建立阵列SAR的线性观测矩阵：
步骤3.1、根据公式
Figure GDA0003389215180000056
计算得到第n个阵列天线在第l个方位向慢时刻的位置矢量，记为
Figure GDA0003389215180000061
其中N为步骤1中初始化得到的阵列天线阵元总数，其中K为步骤1中初始化得到的方位向慢时刻总数，其中
Figure GDA0003389215180000062
为步骤1中初始化得到的阵列天线各阵元初始位置，其中
Figure GDA0003389215180000063
为步骤1中初始化得到的平台速度，其中PRF为步骤1中初始化得到的雷达系统的脉冲重复频率；
步骤3.2、采用公式
Figure GDA0003389215180000064
Figure GDA0003389215180000065
计算得到第l个方位向慢时刻阵列SAR观测场景目标空间Ω中第t个等距离单元格到第n个天线阵元的距离，记为
Figure GDA0003389215180000066
其中M为步骤2中初始化得到的阵列平面维成像空间的第t个等距离单元格阵列向单元格总数，其中||·||2表示定义3中向量L2范数，其中
Figure GDA0003389215180000067
为步骤2中初始化得到的阵列平面维成像空间中第t个等距离单元格中第m个元素的位置，
Figure GDA0003389215180000068
为步骤3.1中得到的第n个阵列天线在第l个方位向慢时刻的位置矢量，其中T为步骤1中初始化得到的距离向快时刻总数，其中K为步骤1中初始化得到的方位向慢时刻总数，其中N为步骤1中初始化得到的天线阵元总数；
步骤3.3、采用公式
Figure GDA0003389215180000069
计算得到在第l个方位向慢时刻阵列SAR观测场景目标空间Ω中第t个等距离单元格到第n个阵列阵元的时间延时，记为
Figure GDA00033892151800000610
其中C为步骤1中初始化得到的电磁波在空气中的传播速度，其中T为步骤1中初始化得到的距离向快时刻总数，其中
Figure GDA00033892151800000611
为步骤3.2中得到的第l个方位向慢时刻阵列SAR观测场景目标空间Ω中第t个等距离单元格到第n个天线阵元的距离，其中K为步骤1中初始化得到的方位向慢时刻总数，其中N为步骤1中初始化得到的阵列天线阵元总数；
步骤3.4、在第l个方位向慢时刻和第t个距离向快时刻中阵列SAR第n个天线阵元的原始回波数据记做s(t,l,n)，t＝1,2,...T，l＝1,2,...K，n＝1,2,...,N，其中T为步骤1中初始化得到的距离向快时刻总数，其中K为步骤1中初始化得到的方位向慢时刻总数，其中N为步骤1中初始化得到的阵列天线阵元总数；在阵列SAR实际成像中，原始回波数据s(t,l,n)由数据接收机提供；
步骤3.5、采用标准合成孔径雷达距离向脉冲压缩方法对s(t,l,n)进行距离向脉冲压缩，得到距离向压缩后的阵列合成孔径雷达数据，记做sAC(t,l,n)，其中为s(t,l,n)步骤3.4中得到的原始回波数据；
根据公式St＝sAC(t,l,n)，t＝1,2,...T，l＝1,2,...K，n＝1,2,...,N计算得到第t个等距离单元格回波信号向量，记为St，St由W＝K·N行1列组成，其中K是步骤1中初始化得到的慢时刻总数，其中N为步骤1中初始化得到的阵列天线的阵元总数，其中T为步骤1中初始化得到的距离向快时刻总数；
步骤3.6、采用公式
Figure GDA0003389215180000071
Figure GDA0003389215180000072
计算得到阵列平面中第m个单元格在慢时间l到回波信号向量St中第i个元素信号对应的时延函数，记为Φi(m)，其中
Figure GDA0003389215180000073
为步骤3.3中得到的在第l个方位向慢时刻阵列SAR观测场景目标空间Ω中第t个等距离单元格到第n个阵列阵元的时间延时；
根据公式Ψ＝Φi(m)，m＝1,2,...M，i＝1,2,...W，计算得到回波信号向量St与散射系数矩阵δ之间的测量矩阵，记做Ψ，其中T为步骤1中初始化得到的距离向快时刻总数，其中K是步骤1中初始化得到的慢时刻总数，其中δ为步骤2中初始化得到的散射系数矩阵，其中St为步骤3.5中得到的第t个等距离单元格回波信号向量，其中M为步骤2中初始化得到的阵列平面维成像空间的第t个等距离单元格阵列向单元格总数，其中W是步骤3.5中得到的第t个等距离单元格回波信号向量St的行数；
步骤4、初始化基于贝叶斯迭代自适应重加权范数最小化的稀疏自聚焦阵列SAR成像算法的初始参数：
步骤4.1、初始化算法的最大迭代次数Nmax，初始化误差迭代终止门限为ε0，初始化加权输入参数为ε1，初始化范数项系数为p，初始化算法迭代次数，记为gen；
步骤4.2、初始化相位误差：
初始化相位误差向量，记做
Figure GDA0003389215180000081
根据公式
Figure GDA0003389215180000082
计算得到第gen次迭代的相位误差矩阵，记做
Figure GDA0003389215180000083
其中W是步骤3.5中得到的原始回波信号向量的行数，其中Nmax为步骤4.1中初始化得到的最大迭代次数；
步骤4.3、初始化散射系数：
采用公式
Figure GDA0003389215180000084
初始化阵列SAR平面维第t个等距离子平面空间散射系数向量，记为
Figure GDA0003389215180000085
其中T为步骤1中初始化得到的距离向快时间总数，其中M为步骤2中初始化得到的阵列平面维成像空间的第t个等距离单元格阵列向单元格总数，其中
Figure GDA0003389215180000086
为步骤4.2中初始化得到的相位误差矩阵，其中Ψ为步骤3.6中得到的回波信号与散射系数之间的测量矩阵，其中Nmax为步骤4.1中初始化得到的最大迭代次数，其中St为步骤3.5中得到的第t个等距离单元格回波信号向量，其中St为步骤3.5中得到的第t个等距离单元格回波信号向量；
步骤4.4、初始化噪声方差：
采用公式
Figure GDA0003389215180000087
初始化系统噪声方差，记做
Figure GDA0003389215180000088
其中Ψ为步骤3.6中得到的测量矩阵，其中
Figure GDA0003389215180000089
为步骤4.2中初始化得到的相位误差矩阵，其中
Figure GDA00033892151800000810
为步骤4.3中初始化得到的第t个等距离子平面空间散射系数向量，其中St为步骤3.5中得到的第t个等距离单元格回波信号向量，其中W是步骤3.5中得到的原始回波信号向量S的行数，其中Ψ为步骤3.6中得到的回波信号与散射系数之间的测量矩阵，其中Nmax为步骤4.1中初始化得到的最大迭代次数，其中T为步骤1中初始化得到的距离向快时间总数；
步骤5、估计散射系数向量、系统噪声方差和相位误差：
步骤5.1、根据噪声方差、相位误差，估计散射系数向量：
在第gen次迭代中，若gen＝0，则散射系数向量为
Figure GDA00033892151800000811
噪声方差为
Figure GDA0003389215180000091
相位误差为
Figure GDA0003389215180000092
其中T为步骤1中初始化得到的距离向快时间总数，其中W是步骤3.5中得到的第个等距离单元格回波信号向量St的行数；
若gen≥1时，根据公式
Figure GDA0003389215180000093
计算得到的第n次迭代的相位误差矩阵，其中
Figure GDA0003389215180000094
为第gen-1次迭代得到的相位误差；
根据公式
Figure GDA0003389215180000095
计算得到第t个等距离子平面空间第gen-1次迭代的对角矩阵，记为
Figure GDA0003389215180000096
其中
Figure GDA0003389215180000097
为第gen-1次迭代后得到的第t个等距离子平面空间散射系数向量，其中ε1为步骤4.1中初始化得到的加权输入参数，其中p为步骤4.1中初始化得到的范数项系数；
采用公式
Figure GDA0003389215180000098
计算得到第gen次迭代后的第t个等距离子平面空间散射系数向量，记为
Figure GDA0003389215180000099
其中
Figure GDA00033892151800000910
为第gen-1次迭代后得到的噪声方差，其中Ψ为步骤3.6中得到的回波信号与散射系数之间的测量矩阵，其中St为步骤3.5中得到的第t个等距离单元格回波信号向量，其中M为步骤2中初始化得到的阵列平面维成像空间的第t个等距离单元格阵列向单元格总数，其中T为步骤1中初始化得到的距离向快时间总数，其中W是步骤3.5中得到的第t个等距离单元格回波信号向量St的行数，其中Nmax为步骤4.1中初始化得到的最大迭代次数；
步骤5.2、根据散射系数、相位误差，估计噪声方差：
采用公式
Figure GDA00033892151800000911
计算得到第gen次迭代后得到噪声方差估计值，记为
Figure GDA00033892151800000912
其中
Figure GDA00033892151800000913
为步骤5.1中得到的第gen次第t个等距离子平面空间散射系数向量，其中Ψ为步骤3.6中得到的回波信号与散射系数之间的测量矩阵，其中St为步骤3.5中得到的第t个等距离单元格回波信号向量，其中W是步骤3.5中得到的原始回波信号向量S的行数；其中
Figure GDA00033892151800000914
为步骤5.1中得到的相位误差矩阵，其中Nmax为步骤4.1中初始化得到的最大迭代次数,其中T为步骤1中初始化得到的距离向快时间总数；
步骤5.3、根据散射系数、噪声方差，估计相位误差向量：
采用公式
Figure GDA0003389215180000101
计算得到第gen次迭代后得到的相位误差估计向量，记做
Figure GDA0003389215180000102
其中S为步骤3.5中得到的回波信号向量，其中Ψ为步骤3.6得到的回波信号与散射系数之间的测量矩阵，其中
Figure GDA0003389215180000103
为步骤5.1中第gen次迭代后得到的第t个等距离子平面空间散射系数向量，其中Nmax为步骤4.1中初始化得到的最大迭代次数，其中∠表示求角度运算符号,其中T为步骤1中初始化得到的距离向快时间总数；
步骤6、迭代终止判断：
如果
Figure GDA0003389215180000104
并且gen≤Nmax，m＝1,2,…,M，t＝1,2,...,T，则继续执行步骤5～6，且gen＝gen+1，其中M为步骤2中初始化得到的阵列平面维成像空间的第t个等距离单元格阵列向单元格总数，其中T为步骤1中初始化的距离向快时间总数，其中ε0为步骤4.1中初始化得到的误差迭代终止门限，其中gen为步骤4.1初始化得到的IARNSABR算法迭代次数；
若不满足
Figure GDA0003389215180000105
和gen≤Nmax任一条件，算法迭代终止，则输出
Figure GDA0003389215180000106
Figure GDA0003389215180000107
得到的IARNSABR算法第gen次迭代得到的散射系数向量值δt即为阵列SAR平面维成像空间最终的第t个等距离子平面空间散射系数向量，其中
Figure GDA0003389215180000108
为步骤5.1中的得到第gen次迭代得到的第t个等距离子平面空间散射系数向量，其中
Figure GDA0003389215180000109
为步骤5.3中得到的第gen次迭代的相位误差，其中Nmax为步骤4.1中初始化得到的最大迭代次数；
步骤7、全场景三维成像：
采用公式AA＝[δ1，…，δT]，将各等距离子平面空间散射系数向量δt，t＝1,2,...,T排列为三维矩阵形式，得到三维SAR观测场景目标区间的三维成像结果，记为AA，其中T为步骤1中初始化得到的距离向快时间总数，其中δt为步骤6中得到的第t个等距离子平面空间散射系数向量；
至此，我们全场景阵列SAR的三维成像结果，整个重构方法结束。
本发明的创新点是：针对阵列SAR回波信号中存在的相位误差对成像结果的影响，本发明在定义6中的SAFBRIM算法的基础上，通过建立阵列SAR原始回波信号与观测场景目标空间中散射系数的线性测量矩阵，对算法中代价函数中的范数项进行迭代自适应重加权处理，对距离向进行脉冲压缩、划分等距离面，然后再对每一个等距离的二维平面进行估计。
本发明优点在于算法可以针对不同的信号进行不同的先验概率建模，在定义6中所定义的SAFBRIM算法的基础上对算法中代价函数中散射系数的范数项进行迭代自适应重加权处理，在重构过程中不断利用重构结果优化范数项自身加权系数，然后对图像进行重构，本算法具有重构精度高、有效降低相位误差的优势，本发明可适用于阵列合成孔径雷达成像等领域。
附图说明
图1为本发明流程图；
图2为系统参数表；
具体实施方式
本发明主要采用计算机仿真的方法进行验证，所有步骤、结论都在MATLAB-2017b上验证正确。具体实施步骤如下：
步骤1、初始化SAR系统参数：
初始化SAR系统参数包括：平台速度矢量记为
Figure GDA0003389215180000111
阵列天线各阵元初始位置矢量，记做
Figure GDA0003389215180000112
其中n为天线各阵元序号，N＝64为阵列天线的阵元总数；阵列天线长度，记做L＝3m；雷达发射信号载频为fc＝37.5GHz；雷达发射信号的调频斜率为fdr＝4×1014Hz/s；脉冲重复时间记为PRI＝2μs；雷达系统的脉冲重复频率为PRF＝0.5MHz；雷达发射信号带宽记做Br＝0.8GHz；电磁波在空气中的传播速度记做C＝3×108m/s；距离向快时间记做t，t＝1,2...T，T＝256为距离向快时刻总数，方位向慢时刻记做l，l＝1,2,...K，K＝64为方位向慢时刻总数；上述参数均为SAR系统标准参数，其中雷达信号载频fc，雷达发射信号的调频斜率fdr，脉冲重复时间PRI，雷达系统的脉冲重复频率PRF，雷达发射信号带宽Br，阵列天线的阵元总数N，阵列天线长度L在阵列SAR系统设计过程中已经确定；平台速度矢量记为
Figure GDA0003389215180000121
阵列天线各阵元初始位置矢量
Figure GDA0003389215180000122
在SAR观测方案设计中已经确定；根据SAR成像系统方案和观测方案，SAR成像方法需要的初始化成像系统参数均为已知；
步骤2、初始化SAR的观测场景目标空间参数：
初始化SAR的观测场景目标空间参数包括：以雷达波束照射场区域地平面和垂直于该地平面向上的单位向量所构成的空间直角坐标系作为阵列SAR的观测场景目标空间Ω，Ω为51×51×256像素；将观测场景目标空间Ω均匀划分为大小相等的立体单元格，称为分辨单元，单元网格在水平横向、水平纵向和高度向边长分别记为dx＝1，dy＝1和dz＝0.12，观测场景空间在水平横向、水平纵向和高度向单元格数分别为Mx＝51，My＝51和Mz＝256；水平横向和水平纵向构成阵列维成像空间，在阵列平面维成像空间上第t个等距离单元格第m个元素的位置，记做
Figure GDA0003389215180000123
其中m＝51(my-1)+mx＝1,…,M，mx＝1,…,51，my＝1,…,51，t＝1,…,T；将观测场景目标空间中第t个等距离单元格第m个元素的散射系数记为
Figure GDA0003389215180000124
Figure GDA0003389215180000125
根据公式
Figure GDA0003389215180000126
得到散射系数矩阵，记做δ，散射系数矩阵δ由M行T列组成，其中T＝256为步骤1中初始化得到的距离向快时刻总数，其中M＝Mx·My＝2601为阵列平面维成像空间的第t个等距离单元格阵列向单元格总数；初始化SAR的观测场景目标空间参数在SAR成像方案设计中已经确定；
步骤3、建立阵列SAR的线性观测矩阵：
步骤3.1、根据公式
Figure GDA0003389215180000127
计算得到第n个阵列天线在第l个方位向慢时刻的位置矢量，记为
Figure GDA0003389215180000128
其中N＝64为步骤1中初始化得到的阵列天线阵元总数，其中K＝64为步骤1中初始化得到的方位向慢时刻总数，其中
Figure GDA0003389215180000131
为步骤1中初始化得到的各阵元初始位置，其中
Figure GDA0003389215180000132
为步骤1中初始化得到的平台速度，其中PRF＝0.5MHz为步骤1中初始化得到的雷达系统的脉冲重复频率；
步骤3.2、采用公式
Figure GDA0003389215180000133
Figure GDA0003389215180000134
计算得到第l个方位向慢时刻阵列SAR观测场景目标空间Ω中第t个等距离单元格到第n个天线阵元的距离，记为
Figure GDA0003389215180000135
其中M＝2601为步骤2中初始化得到的阵列平面维成像空间中第t个等距离单元格阵列向单元格总数，其中||·||2表示定义3中向量L2范数，其中
Figure GDA0003389215180000136
为步骤2中初始化得到的阵列平面维成像空间中第t个等距离单元格中第m个元素的位置，
Figure GDA0003389215180000137
为步骤3.1中得到的第n个阵列天线在第l个方位向慢时刻的位置矢量，其中T＝256为步骤1中初始化得到的距离向快时刻总数，其中K＝64为步骤1中初始化得到的方位向慢时刻总数，其中N＝64为步骤1中初始化得到的阵列天线的阵元总数；
步骤3.3、采用公式
Figure GDA0003389215180000138
计算得到在第l个方位向慢时刻阵列SAR观测场景目标空间Ω中第t个等距离单元格到第n个阵列阵元的时间延时，记为
Figure GDA0003389215180000139
其中C＝3×108m/s为步骤1中初始化得到的电磁波在空气中的传播速度，其中T＝256为步骤1中初始化得到的距离向快时刻总数，其中
Figure GDA00033892151800001310
为步骤3.2中得到的第l个方位向慢时刻阵列SAR观测场景目标空间Ω中第t个等距离单元格到第n个天线阵元的距离，其中K＝64为步骤1中初始化得到的方位向慢时刻总数，其中N＝64为步骤1中初始化得到的阵列天线的阵元总数；
步骤3.4、在第l个方位向慢时刻和第t个距离向快时刻中阵列SAR第n个天线阵元的原始回波数据记做s(t,l,n)，t＝1,2,...T，l＝1,2,...K，n＝1,2,…,N，其中T＝256为步骤1中初始化得到的距离向快时刻总数，其中K＝64为步骤1中初始化得到的方位向慢时刻总数，其中N＝64为步骤1中初始化得到的阵列天线的阵元总数，在阵列SAR实际成像中，s(t,l,n)由数据接收机提供；
步骤3.5、采用定义2中的标准合成孔径雷达距离向脉冲压缩方法对s(t,l,n)进行距离向脉冲压缩后，得到距离向压缩后的阵列合成孔径雷达数据，记做sAC(t,l,n)，其中为s(t,l,n)步骤3.4中得到的原始回波数据；
根据公式St＝sAC(t,l,n)，t＝1,2,…T，l＝1,2,…K，n＝1,2,...,N计算得到第t个等距离单元格回波信号向量，记为St，St由W＝K·N＝4096行1列组成，其中K＝64是步骤1中初始化得到的慢时刻总数，其中N＝64为步骤1中初始化得到的阵列天线的阵元总数，其中T＝256为步骤1中初始化得到的距离向快时刻总数；
步骤3.6、采用公式
Figure GDA0003389215180000141
Figure GDA0003389215180000142
得到阵列平面中第m个单元格在慢时间l到回波信号向量St中第i个元素信号对应的时延函数Φi(m)，其中
Figure GDA0003389215180000143
为在第l个方位向慢时刻阵列SAR观测场景目标空间Ω中第t个等距离单元格到第n个阵列阵元的时间延时；
根据公式Ψ＝Φi(m)，m＝1,2,…M，i＝1,2,…W计算得到回波信号向量St与散射系数向量δ之间的测量矩阵，记做Ψ，其中T＝256为步骤1中初始化得到的距离向快时刻总数，其中K＝64是步骤1中初始化得到的慢时刻总数，其中δ为步骤2中初始化得到的散射系数矩阵，其中St为步骤3.5中得到的第t个等距离单元格回波信号向量，其中M＝2601为步骤2中初始化得到的阵列平面维成像空间中第t个等距离单元格阵列向单元格总数，其中W＝4096是步骤3.5中得到的第t个等距离单元格回波信号向量St的行数，其中N＝64为步骤1中初始化得到的阵列天线的阵元总数；
步骤4、初始化基于贝叶斯迭代自适应重加权范数最小化的稀疏自聚焦阵列SAR成像算法的初始参数：
步骤4.1、初始化算法的最大迭代次数Nmax＝200，初始化误差迭代终止门限为ε0＝10-18，初始化加权输入参数为ε1＝10-2，初始化范数项系数为p＝1，初始化算法迭代次数，记为gen＝0；
步骤4.2、初始化相位误差：
初始化相位误差向量，记做
Figure GDA0003389215180000151
根据公式
Figure GDA0003389215180000152
计算得到第n次迭代的相位误差矩阵，记做
Figure GDA0003389215180000153
其中W＝4096是步骤3.5中得到的第t个等距离单元格回波信号向量St的行数，其中Nmax＝256为步骤4.1中初始化得到的最大迭代次数；
步骤4.3、初始化散射系数：
采用公式
Figure GDA0003389215180000154
初始化阵列SAR平面维第t个等距离子平面空间散射系数向量，记为
Figure GDA0003389215180000155
其中T＝256为步骤1中初始化得到的距离向快时间总数，其中M＝2601为步骤2中初始化得到的阵列平面维成像空间中第t个等距离单元格阵列向单元格总数，其中
Figure GDA0003389215180000156
为步骤4.2中初始化得到的相位误差矩阵，其中Ψ为步骤3.6中得到的回波信号与散射系数之间的测量矩阵，其中Nmax＝200为步骤4.1中初始化得到的最大迭代次数，其中St为步骤3.5中得到的第t个等距离单元格回波信号向量；
步骤4.4、初始化噪声方差：
采用公式
Figure GDA0003389215180000157
初始化系统噪声方差，记做
Figure GDA0003389215180000158
其中Ψ为步骤3.6中得到的测量矩阵，其中
Figure GDA0003389215180000159
为步骤4.2中初始化得到的相位误差矩阵，其中
Figure GDA00033892151800001510
为步骤4.3中初始化得到的第t个等距离子平面空间散射系数向量，其中St为步骤3.5中得到的第t个等距离单元格回波信号向量，其中W＝4096是步骤3.5中得到的第t个等距离单元格回波信号向量St的行数，其中Ψ为步骤3.6中得到的回波信号与散射系数之间的测量矩阵，其中Nmax＝200为步骤4.1中初始化得到的最大迭代次数，其中T＝256为步骤1中初始化得到的距离向快时间总数；
步骤5、估计散射系数向量、系统噪声方差和相位误差：
步骤5.1、根据噪声方差、相位误差，估计散射系数向量：
在第gen次迭代中，若gen＝0，则散射系数向量为
Figure GDA00033892151800001511
噪声方差为
Figure GDA0003389215180000161
相位误差为
Figure GDA0003389215180000162
其中T＝256为步骤1中初始化得到的距离向快时间总数，其中W＝4096是步骤3.5中得到的第t个等距离单元格回波信号向量St的行数,其中St为步骤3.5中得到的第t个等距离单元格回波信号向量；若gen≥1时，根据公式
Figure GDA0003389215180000163
计算得到的第n次迭代的相位误差矩阵，其中
Figure GDA0003389215180000164
为第gen-1次迭代得到的相位误差；
根据公式
Figure GDA0003389215180000165
计算得到第t个等距离子平面空间第gen-1次迭代的对角矩阵，记为
Figure GDA0003389215180000166
其中
Figure GDA0003389215180000167
为第gen-1次迭代后得到的第t个等距离子平面空间散射系数向量，其中ε1＝10-2为步骤4.1中初始化得到的加权输入参数，其中p＝1为步骤4.1中初始化得到的范数项系数；
利用公式
Figure GDA0003389215180000168
计算得到第gen次迭代后的第t个等距离子平面空间散射系数向量，记为
Figure GDA0003389215180000169
其中
Figure GDA00033892151800001610
为第gen-1次迭代后得到的噪声方差，其中Ψ为步骤3.6中得到的回波信号与散射系数之间的测量矩阵，其中St为步骤3.5中得到的第t个等距离单元格回波信号向量，其中M＝2601为步骤2中初始化得到的阵列平面维成像空间中第t个等距离单元格阵列向单元格总数，其中T＝256为步骤1中初始化得到的距离向快时间总数，其中W是步骤3.5中得到的第t个等距离单元格回波信号向量St的行数，其中Nmax＝200为步骤4.1中初始化得到的最大迭代次数；
步骤5.2、根据散射系数、相位误差，估计噪声方差：
根据公式
Figure GDA00033892151800001611
计算得到第gen次迭代后得到噪声方差估计值，记为
Figure GDA00033892151800001612
其中
Figure GDA00033892151800001613
为步骤5.1中得到的第gen次第t个等距离子平面空间散射系数向量，其中Ψ为步骤3.6中得到的回波信号与散射系数之间的测量矩阵，其中St为步骤3.5中得到的第t个等距离单元格回波信号向量，其中W＝4096是步骤3.5中得到的第t个等距离单元格回波信号向量St的行数，；其中
Figure GDA00033892151800001614
为步骤5.1中得到的相位误差矩阵，其中Nmax＝200为步骤4.1中初始化得到的最大迭代次数,其中T＝256为步骤1中初始化得到的距离向快时间总数；
步骤5.3、根据散射系数、噪声方差，估计相位误差向量：
根据公式
Figure GDA0003389215180000171
计算得到第gen次迭代后得到的相位误差估计向量，记做
Figure GDA0003389215180000172
其中St为步骤3.5中得到的第t个等距离单元格回波信号向量，其中Ψ为步骤3.6得到的回波信号与散射系数之间的测量矩阵，其中
Figure GDA0003389215180000173
为步骤5.1中第gen次迭代后得到的第t个等距离子平面空间散射系数向量，其中Nmax＝200为步骤4.1中初始化得到的最大迭代次数，其中∠表示求角度运算符号,其中T＝256为步骤1中初始化得到的距离向快时间总数；
步骤6、迭代终止判断：
如果
Figure GDA0003389215180000174
并且gen≤Nmax，m＝1,2,…,M，t＝1,2,…,T，则继续执行步骤5～6，且gen＝gen+1，其中M＝2601为步骤2中初始化得到的阵列平面维成像空间中第t个等距离单元格阵列向单元格总数，其中T＝256为步骤1中初始化的距离向快时间总数，其中ε0＝10-18为步骤4.1中初始化得到的误差迭代终止门限，其中gen＝0为步骤4.1初始化得到的IARNSABR算法迭代次数；
若不满足
Figure GDA0003389215180000175
和gen≤Nmax任一条件，算法迭代终止，则输出
Figure GDA0003389215180000176
Figure GDA0003389215180000177
得到的IARNSABR算法第gen次迭代得到的散射系数向量值δt即为阵列SAR平面维成像空间最终的第t个等距离子平面空间散射系数向量，其中
Figure GDA0003389215180000178
为步骤5.1中的得到第gen次迭代得到的第t个等距离子平面空间散射系数向量，其中
Figure GDA0003389215180000179
为步骤5.3中得到的第gen次迭代的相位误差，其中Nmax＝200为步骤4.1中初始化得到的最大迭代次数；
步骤7、全场景三维成像：
根据公式AA＝[δ1，…，δT]将各等距离子平面空间散射系数向量δt，t＝1,2,…,T排列为三维矩阵形式，得到三维SAR观测场景目标区间的三维成像结果，记为AA，其中T＝256为步骤1中初始化得到的距离向快时间总数，其中δt为步骤6中得到的第t个等距离子平面空间散射系数向量；
至此，我们全场景阵列SAR的三维成像结果，整个重构方法结束。
经过计算机仿真及实测数据结果证明，本发明通过对代价函数中的范数项进行迭代自适应重加权处理，对每个距离子空间的散射系数赋予不同的加权系数，本发明能够更好的校正相位误差对于高精度阵列SAR成像的影响，获得了具有更高质量的成像结果。

Claims (1)

1.一种贝叶斯迭代重加权稀疏自聚焦阵列SAR成像方法，其特征是它包括如下步骤：
步骤1、初始化SAR系统参数：
初始化SAR系统参数包括：平台速度矢量记为
Figure FDA0003478944120000011
阵列天线各阵元初始位置矢量，记做
Figure FDA0003478944120000012
其中n为天线各阵元序号，N为阵列天线的阵元总数；阵列天线长度，记做L；雷达发射信号载频为fc；雷达发射信号的调频斜率为fdr；脉冲重复时间记为PRI；雷达系统的脉冲重复频率为PRF；雷达发射信号带宽记做Br；电磁波在空气中的传播速度记做C；距离向快时间记做t，t＝1,2...T，T为距离向快时刻总数，方位向慢时刻记做l，l＝1,2,...K，K为方位向慢时刻总数；上述参数均为SAR系统标准参数，其中雷达信号载频fc，雷达发射信号的调频斜率fdr，脉冲重复时间PRI，雷达系统的脉冲重复频率PRF，雷达发射信号带宽Br，阵列天线的阵元总数N，阵列天线长度L在阵列SAR系统设计过程中已经确定；平台速度矢量记为
Figure FDA0003478944120000013
阵列天线各阵元初始位置矢量
Figure FDA0003478944120000014
在SAR观测方案设计中已经确定；根据SAR成像系统方案和观测方案，SAR成像方法需要的初始化成像系统参数均为已知；
步骤2、初始化SAR的观测场景目标空间参数：
初始化SAR的观测场景目标空间参数包括：以雷达波束照射场区域地平面和垂直于该地平面向上的单位向量所构成的空间直角坐标系作为阵列SAR的观测场景目标空间Ω；将观测场景目标空间Ω均匀划分为大小相等的立体单元格，称为分辨单元，单元网格在水平横向、水平纵向和高度向边长分别记为dx，dy和dz，观测场景空间在水平横向、水平纵向和高度向单元格数分别为Mx,My和Mz，单元格大小为阵列SAR系统传统理论成像分辨率；水平横向和水平纵向构成阵列维成像空间，在阵列平面维成像空间上第t个等距离单元格第m个元素的位置，记做
Figure FDA0003478944120000015
其中m＝(my-1)Mx+mx＝1,…,M，mx＝1,…,Mx，my＝1,…,My，t＝1,…,T；将观测场景目标空间中第t个等距离单元格第m个元素的散射系数记为
Figure FDA0003478944120000016
根据公式
Figure FDA0003478944120000021
计算得到散射系数矩阵，记做δ，散射系数矩阵δ由M行T列组成，其中T为步骤1中初始化得到的距离向快时刻总数，其中M＝Mx·My为阵列平面维成像空间的第t个等距离单元格阵列向单元格总数；初始化SAR的观测场景目标空间参数在SAR成像方案设计中已经确定
步骤3、建立阵列SAR的线性观测矩阵：
步骤3.1、根据公式
Figure FDA0003478944120000022
计算得到第n个阵列天线在第l个方位向慢时刻的位置矢量，记为
Figure FDA0003478944120000023
其中N为步骤1中初始化得到的阵列天线阵元总数，其中K为步骤1中初始化得到的方位向慢时刻总数，其中
Figure FDA0003478944120000024
为步骤1中初始化得到的阵列天线各阵元初始位置，其中
Figure FDA0003478944120000025
为步骤1中初始化得到的平台速度，其中PRF为步骤1中初始化得到的雷达系统的脉冲重复频率；
步骤3.2、采用公式
Figure FDA0003478944120000026
Figure FDA0003478944120000027
计算得到第l个方位向慢时刻阵列SAR观测场景目标空间Ω中第t个等距离单元格到第n个天线阵元的距离，记为
Figure FDA0003478944120000028
其中M为步骤2中初始化得到的阵列平面维成像空间的第t个等距离单元格阵列向单元格总数，其中||·||2表示向量L2范数，其中
Figure FDA0003478944120000029
为步骤2中初始化得到的阵列平面维成像空间中第t个等距离单元格中第m个元素的位置，
Figure FDA00034789441200000210
为步骤3.1中得到的第n个阵列天线在第l个方位向慢时刻的位置矢量，其中T为步骤1中初始化得到的距离向快时刻总数，其中K为步骤1中初始化得到的方位向慢时刻总数，其中N为步骤1中初始化得到的天线阵元总数；
步骤3.3、采用公式
Figure FDA00034789441200000211
计算得到在第l个方位向慢时刻阵列SAR观测场景目标空间Ω中第t个等距离单元格到第n个阵列阵元的时间延时，记为
Figure FDA00034789441200000212
其中C为步骤1中初始化得到的电磁波在空气中的传播速度，其中T为步骤1中初始化得到的距离向快时刻总数，其中
Figure FDA00034789441200000213
为步骤3.2中得到的第l个方位向慢时刻阵列SAR观测场景目标空间Ω中第t个等距离单元格到第n个天线阵元的距离，其中K为步骤1中初始化得到的方位向慢时刻总数，其中N为步骤1中初始化得到的阵列天线阵元总数；
步骤3.4、在第l个方位向慢时刻和第t个距离向快时刻中阵列SAR第n个天线阵元的原始回波数据记做s(t,l,n)，t＝1,2,...T，l＝1,2,...K，n＝1,2,...,N，其中T为步骤1中初始化得到的距离向快时刻总数，其中K为步骤1中初始化得到的方位向慢时刻总数，其中N为步骤1中初始化得到的阵列天线阵元总数；在阵列SAR实际成像中，原始回波数据s(t,l,n)由数据接收机提供；
步骤3.5、采用标准合成孔径雷达距离向脉冲压缩方法对s(t,l,n)进行距离向脉冲压缩，得到距离向压缩后的阵列合成孔径雷达数据，记做sAC(t,l,n)，其中为s(t,l,n)步骤3.4中得到的原始回波数据；
根据公式St＝sAC(t,l,n)，t＝1,2,...T，l＝1,2,...K，n＝1,2,...,N计算得到第t个等距离单元格回波信号向量，记为St，St由W＝K·N行1列组成，其中K是步骤1中初始化得到的慢时刻总数，其中N为步骤1中初始化得到的阵列天线的阵元总数，其中T为步骤1中初始化得到的距离向快时刻总数；
步骤3.6、采用公式
Figure FDA0003478944120000031
Figure FDA0003478944120000032
计算得到阵列平面中第m个单元格在慢时间l到回波信号向量St中第i个元素信号对应的时延函数，记为Φi(m)，其中
Figure FDA0003478944120000033
为步骤3.3中得到的在第l个方位向慢时刻阵列SAR观测场景目标空间Ω中第t个等距离单元格到第n个阵列阵元的时间延时；
根据公式Ψ＝Φi(m)，m＝1,2,...M，i＝1,2,...W，计算得到回波信号向量St与散射系数矩阵δ之间的测量矩阵，记做Ψ，其中T为步骤1中初始化得到的距离向快时刻总数，其中K是步骤1中初始化得到的慢时刻总数，其中δ为步骤2中初始化得到的散射系数矩阵，其中St为步骤3.5中得到的第t个等距离单元格回波信号向量，其中M为步骤2中初始化得到的阵列平面维成像空间的第t个等距离单元格阵列向单元格总数，其中W是步骤3.5中得到的第t个等距离单元格回波信号向量St的行数；
步骤4、初始化基于贝叶斯迭代自适应重加权范数最小化的稀疏自聚焦阵列SAR成像算法的初始参数：
步骤4.1、初始化算法的最大迭代次数Nmax，初始化误差迭代终止门限为ε0，初始化加权输入参数为ε1，初始化范数项系数为p，初始化算法迭代次数，记为gen；
步骤4.2、初始化相位误差：
初始化相位误差向量，记做
Figure FDA0003478944120000041
根据公式
Figure FDA0003478944120000042
计算得到第gen次迭代的相位误差矩阵，记做
Figure FDA0003478944120000043
其中W是步骤3.5中得到的原始回波信号向量的行数，其中Nmax为步骤4.1中初始化得到的最大迭代次数；
步骤4.3、初始化散射系数：
采用公式
Figure FDA0003478944120000044
初始化阵列SAR平面维第t个等距离子平面空间散射系数向量，记为
Figure FDA0003478944120000045
其中T为步骤1中初始化得到的距离向快时间总数，其中M为步骤2中初始化得到的阵列平面维成像空间的第t个等距离单元格阵列向单元格总数，其中
Figure FDA0003478944120000046
为步骤4.2中初始化得到的相位误差矩阵，其中Ψ为步骤3.6中得到的回波信号与散射系数之间的测量矩阵，其中Nmax为步骤4.1中初始化得到的最大迭代次数，其中St为步骤3.5中得到的第t个等距离单元格回波信号向量，其中St为步骤3.5中得到的第t个等距离单元格回波信号向量；
步骤4.4、初始化噪声方差：
采用公式
Figure FDA0003478944120000047
初始化系统噪声方差，记做
Figure FDA0003478944120000048
其中Ψ为步骤3.6中得到的测量矩阵，其中
Figure FDA0003478944120000049
为步骤4.2中初始化得到的相位误差矩阵，其中
Figure FDA00034789441200000410
为步骤4.3中初始化得到的第t个等距离子平面空间散射系数向量，其中St为步骤3.5中得到的第t个等距离单元格回波信号向量，其中W是步骤3.5中得到的原始回波信号向量S的行数，其中Ψ为步骤3.6中得到的回波信号与散射系数之间的测量矩阵，其中Nmax为步骤4.1中初始化得到的最大迭代次数，其中T为步骤1中初始化得到的距离向快时间总数；
步骤5、估计散射系数向量、系统噪声方差和相位误差：
步骤5.1、根据噪声方差、相位误差，估计散射系数向量：
在第gen次迭代中，若gen＝0，则散射系数向量为
Figure FDA0003478944120000051
噪声方差为
Figure FDA0003478944120000052
相位误差为
Figure FDA0003478944120000053
其中T为步骤1中初始化得到的距离向快时间总数，其中W是步骤3.5中得到的第个等距离单元格回波信号向量St的行数；
若gen≥1时，根据公式
Figure FDA0003478944120000054
计算得到的第n次迭代的相位误差矩阵，其中
Figure FDA0003478944120000055
为第gen-1次迭代得到的相位误差；
根据公式
Figure FDA0003478944120000056
计算得到第t个等距离子平面空间第gen-1次迭代的对角矩阵，记为
Figure FDA0003478944120000057
其中
Figure FDA0003478944120000058
为第gen-1次迭代后得到的第t个等距离子平面空间散射系数向量，其中ε1为步骤4.1中初始化得到的加权输入参数，其中p为步骤4.1中初始化得到的范数项系数；
采用公式
Figure FDA0003478944120000059
计算得到第gen次迭代后的第t个等距离子平面空间散射系数向量，记为
Figure FDA00034789441200000510
其中
Figure FDA00034789441200000511
为第gen-1次迭代后得到的噪声方差，其中Ψ为步骤3.6中得到的回波信号与散射系数之间的测量矩阵，其中St为步骤3.5中得到的第t个等距离单元格回波信号向量，其中M为步骤2中初始化得到的阵列平面维成像空间的第t个等距离单元格阵列向单元格总数，其中T为步骤1中初始化得到的距离向快时间总数，其中W是步骤3.5中得到的第t个等距离单元格回波信号向量St的行数，其中Nmax为步骤4.1中初始化得到的最大迭代次数；
步骤5.2、根据散射系数、相位误差，估计噪声方差：
采用公式
Figure FDA0003478944120000061
计算得到第gen次迭代后得到噪声方差估计值，记为
Figure FDA0003478944120000062
其中
Figure FDA0003478944120000063
为步骤5.1中得到的第gen次第t个等距离子平面空间散射系数向量，其中Ψ为步骤3.6中得到的回波信号与散射系数之间的测量矩阵，其中St为步骤3.5中得到的第t个等距离单元格回波信号向量，其中W是步骤3.5中得到的原始回波信号向量S的行数；其中
Figure FDA0003478944120000064
为步骤5.1中得到的相位误差矩阵，其中Nmax为步骤4.1中初始化得到的最大迭代次数,其中T为步骤1中初始化得到的距离向快时间总数；
步骤5.3、根据散射系数、噪声方差，估计相位误差向量：
采用公式
Figure FDA0003478944120000065
计算得到第gen次迭代后得到的相位误差估计向量，记做
Figure FDA0003478944120000066
其中S为步骤3.5中得到的回波信号向量，其中Ψ为步骤3.6得到的回波信号与散射系数之间的测量矩阵，其中
Figure FDA0003478944120000067
为步骤5.1中第gen次迭代后得到的第t个等距离子平面空间散射系数向量，其中Nmax为步骤4.1中初始化得到的最大迭代次数，其中∠表示求角度运算符号,其中T为步骤1中初始化得到的距离向快时间总数；
步骤6、迭代终止判断：
如果
Figure FDA0003478944120000068
并且gen≤Nmax，m＝1,2,...,M，t＝1,2,...,T，则继续执行步骤5～6，且gen＝gen+1，其中M为步骤2中初始化得到的阵列平面维成像空间的第t个等距离单元格阵列向单元格总数，其中T为步骤1中初始化的距离向快时间总数，其中ε0为步骤4.1中初始化得到的误差迭代终止门限，其中gen为步骤4.1初始化得到的IARNSABR算法迭代次数；
若不满足
Figure FDA0003478944120000069
和gen≤Nmax任一条件，算法迭代终止，则输出
Figure FDA00034789441200000610
Figure FDA0003478944120000071
得到的IARNSABR算法第gen次迭代得到的散射系数向量值δt即为阵列SAR平面维成像空间最终的第t个等距离子平面空间散射系数向量，其中
Figure FDA0003478944120000072
为步骤5.1中的得到第gen次迭代得到的第t个等距离子平面空间散射系数向量，其中
Figure FDA0003478944120000073
为步骤5.3中得到的第gen次迭代的相位误差，其中Nmax为步骤4.1中初始化得到的最大迭代次数；
步骤7、全场景三维成像：
采用公式AA＝[δ1，...，δT]，将各等距离子平面空间散射系数向量δt，t＝1,2,...,T排列为三维矩阵形式，得到三维SAR观测场景目标区间的三维成像结果，记为AA，其中T为步骤1中初始化得到的距离向快时间总数，其中δt为步骤6中得到的第t个等距离子平面空间散射系数向量。
CN201810767253.XA 2018-07-13 2018-07-13 一种贝叶斯迭代重加权稀疏自聚焦阵列sar成像方法 Active CN109061642B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810767253.XA CN109061642B (zh) 2018-07-13 2018-07-13 一种贝叶斯迭代重加权稀疏自聚焦阵列sar成像方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810767253.XA CN109061642B (zh) 2018-07-13 2018-07-13 一种贝叶斯迭代重加权稀疏自聚焦阵列sar成像方法

Publications (2)

Publication Number Publication Date
CN109061642A CN109061642A (zh) 2018-12-21
CN109061642B true CN109061642B (zh) 2022-03-15

Family

ID=64816334

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810767253.XA Active CN109061642B (zh) 2018-07-13 2018-07-13 一种贝叶斯迭代重加权稀疏自聚焦阵列sar成像方法

Country Status (1)

Country Link
CN (1) CN109061642B (zh)

Families Citing this family (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109633646B (zh) * 2019-01-21 2022-05-06 中国人民解放军陆军工程大学 一种基于加权l1范数约束的双基地isar成像方法
CN110109101A (zh) * 2019-04-04 2019-08-09 电子科技大学 一种基于自适应阈值的压缩感知三维sar成像方法
CN110133651B (zh) * 2019-05-24 2021-04-06 中国科学院电子学研究所 一种稀疏sar成像自适应稀疏度估计方法、装置
CN110133656B (zh) * 2019-06-06 2022-05-03 电子科技大学 一种基于互质阵列的分解与融合的三维sar稀疏成像方法
CN110764086B (zh) * 2019-09-29 2022-09-09 西安电子科技大学 一种基于扰动矩阵估计的贝叶斯雷达关联成像方法
CN111766575B (zh) * 2020-06-08 2023-04-21 桂林电子科技大学 一种穿墙雷达自聚焦稀疏成像方法及计算机设备
CN113204022B (zh) * 2021-04-30 2022-07-29 电子科技大学 基于相关向量机线性阵列sar三维成像快速贝叶斯压缩感知方法
CN113608218B (zh) * 2021-07-19 2023-05-26 电子科技大学 一种基于后向投影原理的频域干涉相位稀疏重构方法
CN113484862B (zh) * 2021-08-04 2023-10-17 电子科技大学 一种自适应的高分宽幅sar清晰重构成像方法
CN115421115A (zh) * 2022-05-23 2022-12-02 中国人民解放军空军预警学院 一种用于联合相位校正与isar成像的重赋权交替方向乘子法
CN116702514B (zh) * 2023-08-02 2023-09-29 南京纳特通信电子有限公司 基于近电场优化的天线阵列优化方法、装置、介质和设备

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2015008310A1 (en) * 2013-07-19 2015-01-22 Consiglio Nazionale Delle Ricerche Method for filtering of interferometric data acquired by synthetic aperture radar (sar)
CN107037429A (zh) * 2017-04-17 2017-08-11 电子科技大学 基于门限梯度追踪算法的线阵sar三维成像方法
CN108008385A (zh) * 2017-11-20 2018-05-08 西安电子科技大学 基于稀疏贝叶斯学习的干扰环境isar高分辨成像方法
CN108226928A (zh) * 2017-12-18 2018-06-29 西安电子科技大学 基于期望传播算法的逆合成孔径雷达成像方法
CN108226927A (zh) * 2017-12-14 2018-06-29 电子科技大学 基于加权迭代最小稀疏贝叶斯重构算法的sar成像方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7969345B2 (en) * 2009-04-13 2011-06-28 Raytheon Company Fast implementation of a maximum likelihood algorithm for the estimation of target motion parameters

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2015008310A1 (en) * 2013-07-19 2015-01-22 Consiglio Nazionale Delle Ricerche Method for filtering of interferometric data acquired by synthetic aperture radar (sar)
CN107037429A (zh) * 2017-04-17 2017-08-11 电子科技大学 基于门限梯度追踪算法的线阵sar三维成像方法
CN108008385A (zh) * 2017-11-20 2018-05-08 西安电子科技大学 基于稀疏贝叶斯学习的干扰环境isar高分辨成像方法
CN108226927A (zh) * 2017-12-14 2018-06-29 电子科技大学 基于加权迭代最小稀疏贝叶斯重构算法的sar成像方法
CN108226928A (zh) * 2017-12-18 2018-06-29 西安电子科技大学 基于期望传播算法的逆合成孔径雷达成像方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
"Estimation Method for Pol-InSAR Multi-interferometric Phase Based on the MUSIC Method";Zhang Xiao-ling et al.;《Journal of University of Electronic Science and Technology of China》;20110930;652-657 *
"Multi-Static Passive SAR Imaging Based on Bayesian Compressive Sensing";Qisong Wu et al.;《COMPRESSIVE SENSING III》;20141231;1-9 *
"基于压缩感知的连续场景稀疏阵列SAR 三维成像";李烈辰 等;《电子与信息学报》;20140930;第36卷(第9期);2166-2171 *
"基于稀疏贝叶斯正则化的LASAR 高分辨成像算法";韦顺军 等;《第四届高分辨率对地观测学术年会论文集》;20170917;1-20 *

Also Published As

Publication number Publication date
CN109061642A (zh) 2018-12-21

Similar Documents

Publication Publication Date Title
CN109061642B (zh) 一种贝叶斯迭代重加权稀疏自聚焦阵列sar成像方法
CN108226927B (zh) 基于加权迭代最小稀疏贝叶斯重构算法的sar成像方法
CN107037429B (zh) 基于门限梯度追踪算法的线阵sar三维成像方法
CN107193003B (zh) 一种稀疏奇异值分解扫描雷达前视成像方法
CN111679277B (zh) 一种基于sbrim算法的多基线层析sar三维成像方法
CN103149561B (zh) 一种基于场景块稀疏的稀疏微波成像方法
CN106680817B (zh) 一种实现前视雷达高分辨成像的方法
CN111145337B (zh) 基于分辨率逼近的快速稀疏重构的线阵sar三维成像方法
Zhang et al. Superresolution downward-looking linear array three-dimensional SAR imaging based on two-dimensional compressive sensing
CN104536000A (zh) 一种实波束扫描雷达角超分辨方法
CN110109101A (zh) 一种基于自适应阈值的压缩感知三维sar成像方法
CN108226891B (zh) 一种扫描雷达回波计算方法
CN110346794B (zh) 一种资源优化配置的分布式雷达成像方法
CN105699969A (zh) 基于广义高斯约束的最大后验估计角超分辨成像方法
CN104391295A (zh) 一种图像熵最优的压缩传感sar稀疏自聚焦成像方法
CN107843875A (zh) 基于奇异值分解降噪的贝叶斯压缩感知雷达数据融合方法
Ren et al. 3D imaging algorithm for down-looking MIMO array SAR based on Bayesian compressive sensing
CN103728617B (zh) 双基地合成孔径雷达时域快速成像方法
CN113608218B (zh) 一种基于后向投影原理的频域干涉相位稀疏重构方法
Huan et al. Low elevation angle estimation with range super-resolution in wideband radar
CN110133656B (zh) 一种基于互质阵列的分解与融合的三维sar稀疏成像方法
Kang et al. Downward-looking linear array three-dimensional SAR imaging based on the two-dimensional mismatch compensation
CN113835090B (zh) 一种基于多通道sar系统的高精度干涉相位获取方法
CN113204022B (zh) 基于相关向量机线性阵列sar三维成像快速贝叶斯压缩感知方法
CN113238229A (zh) 一种geo星机双基sar无模糊成像方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant