CN110109107B - 一种合成孔径雷达频域bp算法的运动误差补偿方法 - Google Patents

一种合成孔径雷达频域bp算法的运动误差补偿方法 Download PDF

Info

Publication number
CN110109107B
CN110109107B CN201910334110.4A CN201910334110A CN110109107B CN 110109107 B CN110109107 B CN 110109107B CN 201910334110 A CN201910334110 A CN 201910334110A CN 110109107 B CN110109107 B CN 110109107B
Authority
CN
China
Prior art keywords
vector
pixel points
pixel
grid
calculated
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
CN201910334110.4A
Other languages
English (en)
Other versions
CN110109107A (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 CN201910334110.4A priority Critical patent/CN110109107B/zh
Publication of CN110109107A publication Critical patent/CN110109107A/zh
Application granted granted Critical
Publication of CN110109107B publication Critical patent/CN110109107B/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
    • 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

Landscapes

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

Abstract

本发明公开了一种合成孔径雷达频域BP算法的运动误差补偿方法，它是基于图像锐度最优准则通过估计等效相位中心误差进行运动补偿，首先采用雷达平台匀速直线运动轨迹进行并利用频域BP算法进行成像，获得通过频域BP算法成像得到的未聚焦的SAR图像，然后以SAR图像锐度作为目标函数，利用最优化技术估计出天线相位中心位置误差，最后将估计出的APC位置误差加到匀速直线运动轨迹上，获得APC绝对位置，并用频域BP算法进行最终的高精度成像。本发明由于首先在频域进行成像处理，且无需要计算散射点到阵元的距离史，所需运算时间比时域的自聚焦BP算法大大减少，因此更适用于大场景、长孔径、高精度SAR成像。

Description

一种合成孔径雷达频域BP算法的运动误差补偿方法
技术领域
本发明属于合成孔径雷达(Synthetic Aperture Radar,SAR)高分辨率成像的技术领域，它特别涉及到了SAR高精度运动补偿的技术领域。
背景技术
合成孔径雷达(SAR)是一种高分辨率微波成像雷达，具有全天时和全天候工作的优点，已被广泛应用各个领域，如地形测绘、制导、环境遥感和资源勘探等。SAR应用的重要前提和信号处理的主要目标是通过成像算法获取高分辨、高精度的微波图像。但是诸多环境因素(如风场、湍流等)将导致雷达载体平台的运动轨迹偏离设计的理想轨迹，从而严重影响SAR图像的质量(包括聚焦深度、对比度等)。因此，运动补偿技术成为了SAR成像过程中的关键技术。
频域后向投影(BP)算法首先将合成孔径雷达原始数据沿距离向进行距离压缩(脉冲压缩)，从每个孔径的原始回波频谱中获得图像频谱的反投影值，并乘以相位补偿因子。然后，对每个孔径的图像频谱进行累积，重建出最终的图像频谱。最后，通过图像频谱的反傅立叶变换，计算出在时域的图像。由于在已知平台轨迹的前提下，频域BP算法可以有效的补偿运动误差，并且成像效率高，因此成为一种具有应用前景的算法。详见“Zhe Li,JianWang,Qing Huo Liu,Frequency-Domain Backprojection Algorithm for SyntheticAperture Radar Imaging[J].IEEE Geoscience and Remote Sensing Letters.2015,12(4):905-909”。然而，当运动轨迹出现偏差，测量设备精度不能满足要求时，不能精确补偿相位，将导致成像结果散焦。因此研究运动误差补偿算法对频域BP高精度成像具有重要意义。
自聚焦算法有基于相位误差函数的自聚焦算法，如相位梯度自聚焦、子视图相关算法等，和基于雷达成像质量的自聚焦算法，如基于最小熵算法、基于最大锐度和基于最大对比度算法等。详见“皮亦鸣,杨建宇,付毓生,et al.合成孔径雷达成像原理[M].电子科技大学出版社,2007”。其中，基于锐度的自聚焦算法能够实现较好的相位误差补偿，目前有大量基于图像锐度的自聚焦算法。如，基于图像锐度采用自然几何优化的自聚焦BP算法(详见“J.N.Ash,An autofocus method for backprojection imagery in synthetic apertureradar[J].IEEE Geoscience and Remote Sensing Letters.2012,9(1):104-108”)，基于半正定规划最大锐度自聚焦算法(详见“Wei S J,Zhang X L,Hu K B,et al.“LASARautofocus imaging using maximum sharpness back projection via semidefiniteprogramming,IEEE Radar Conference,1311-1315,2015.”)等。但是，这些自聚焦算法都面临着计算量大或无法对整个场景进行精确的相位补偿的问题。且这些方法不适用于频域BP算法，难以对频域BP算法成像进行高精度误差补偿。
发明内容
针对频域后向投影算法，本发明提出了一种合成孔径雷达频域BP算法的运动误差补偿方法，该方法基于图像锐度最优准则通过估计等效相位中心误差进行运动补偿，具有高效和高精度的优点。首先采用雷达平台匀速直线运动轨迹进行并利用频域BP算法进行成像，获得通过频域BP算法成像得到的未聚焦的SAR图像，然后以SAR图像锐度作为目标函数，利用最优化技术估计出天线相位中心(Antenna Phase Center,APC)位置误差，最后将估计出的APC位置误差加到匀速直线运动轨迹上，获得APC绝对位置，并用频域BP算法进行最终的高精度成像。
为了方便描述本发明的内容，首先作以下术语定义：
定义1、脉冲压缩
脉冲压缩是一种现代雷达信号处理技术，简单来说就是雷达发射宽脉冲，然后再接收端“压缩”为窄脉冲，从而改善雷达的两种性能：作用距离和距离分辨率。详见“皮亦鸣,杨建宇,付毓生,杨晓波.合成孔径雷达成像原理.第一版.电子科技大学出版社.2007.3”。
定义2、升采样
升采样是一种在离散信号域提高信号采样率的方法，有时域升采样和频域升采样两种实现方式。
定义3、快速傅里叶变换
计算离散傅里叶变换的一种快速算法，简称FFT。采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少，特别是被变换的抽样点数越多，FFT算法计算量的节省就越显著。FFT的逆变换叫做逆傅里叶变换，简称IFFT。详见“程乾生.数字信号处理.北京大学出版社，北京，2003”。
定义4、频域后向投影算法
频域后向投影算法，简称频域BP算法。频域BP算法首先将原始回波数据在距离向变换到频域获得原始回波频谱，然后计算图像频谱每个像素对应的索引值，根据该索引值可以从每个孔径的原始回波频谱中获得图像频谱的反投影值，并乘以相位补偿因子。然后，对每个孔径的图像频谱进行累积，重建出最终的图像频谱。最后，通过图像频谱的反傅立叶变换，计算出在时域的图像。详见“Zhe Li,Jian Wang,Qing Huo Liu,Frequency-DomainBackprojection Algorithm for Synthetic Aperture Radar Imaging[J].IEEEGeoscience and Remote Sensing Letters.2015,12(4):905-909”。
定义5、方位向、距离向
将雷达平台运动的方向叫做方位向，将垂直于方位向的方向叫做距离向。
定义6、快时间、慢时间、慢时刻
快时间是距离向采样的时间，慢时间是方位向采样的时间，将离散的慢时间从1开始编号，每一个编号叫做一个慢时刻。
定义7、天线相位中心
天线相位中心，简称APC，是雷达天线发射信号的位置，精确的天线相位中心是频域BP算法能够精确成像的前提。
定义8、图像锐度
图像强度是指一幅复图像中每个像素点幅度平方之和，图像锐度则是图像强度的平方。可以用来表征SAR图像聚焦好坏，SAR图像聚焦越好，图像锐度越大。
定义9、共轭梯度法
共轭梯度法是一种最优化方法，用迭代点处的负梯度向量为基础产生一组共轭方向。共轭梯度法的收敛速度快，且不用求矩阵的逆，在使用计算机求解时，所需存储量较小。详见“何坚勇.最优化方法.第一版.清华大学出版社.北京.2007.1”。
定义10、波数矢量
波矢是波的矢量表示方法。波矢是一个矢量，其大小表示波数，其方向表示波传播的方向。
定义11、Armijo算法
Armijo算法是一种一维搜索算法，可以保证目标函数在搜索方向上充分下降。详见“ARMIJO.Minimization of functions having Lipschitz continuous first partialderivatives[J].Pacific Journal of Mathematics.vol.16,no.1,pp.1-3.1966”。
本发明提供了一种合成孔径雷达高精度运动补偿方法，它包括如下步骤(如附图1所示)：
步骤1、相关参数的初始化
初始化的参数均为已知，且初始化的参数如下：光速为C；雷达发射线性调频信号，载频为ω；雷达发射脉冲的带宽为B；雷达发射脉冲的时宽为Tp；雷达脉冲重复周期为T；雷达回波距离向采样频率为fs；雷达回波数据矩阵为SK×L；雷达回波数据矩阵SK×L的方位向点数和距离向点数分别为K和L(K和L均为正整数)，K也称为慢时刻数；慢时刻向量为ts＝[-K/2,1-K/2,…,K/2-1]×T；频域数据矩阵为QSK×L，QSK×L的大小为K行、L列，K为雷达回波数据矩阵的方位向点数，L为雷达回波数据的距离向点数；升采样数据矩阵为SSK×P，SSK×P的大小为K行、8×L列，K为雷达回波数据矩阵的方位向点数，L为雷达回波数据的距离向点数；几何坐标系为三维笛卡尔坐标系O-XY；雷达平台速度向量为V，V的大小为1行3列；雷达平台在零时刻的位置向量为Pt0，Pt0的大小为1行3列；雷达参考距离为R0；O-XY平面内的矩形场景为Θ；将Θ离散化为像素点网格，记为ΩM×N；像素点网格ΩM×N中X方向和Y方向的像素点点数分别为M和N；像素点网格ΩM×N中X方向和Y方向的像素点间隔分别为dx和dy；场景中心位置向量为Pc，Pc的大小为1行3列；共轭梯度算法最大迭代次数为Q，Q为正整数；初始APC误差向量为
Figure GDA0003534881160000041
Figure GDA0003534881160000042
01×2K为1行、2×K列的零矩阵，K为雷达回波数据矩阵SK×L方位向点数；G＝[gx,gy]为波数矢量，gx为距离向波数，gy为方位向波数，离散化形式为gi＝g1+(i-1)Δg,i＝1,...,M，g1＝2πf1/c，Δg＝2πΔf/c，Δg为波数间隔，f1和Δf分别是信号的起始频率和频率间隔。
步骤2、对雷达回波数据矩阵的每一行进行脉冲压缩
取出步骤1中所有雷达回波数据SK×L，采用脉冲压缩方法对SK×L的每一行进行脉冲压缩，得到脉冲压缩后的数据矩阵PSK×L。
步骤3、将脉冲压缩后的数据矩阵转至频域
对步骤2得到的脉冲压缩后的数据矩阵PSK×L做快速傅里叶变换(FFT)：
步骤3.1、取出步骤2中脉冲压缩后的数据矩阵PSK×L的第k行向量，记为ppk，k＝1,2,…,K，K为步骤1定义的慢时刻数。
步骤3.2、对向量ppk做快速傅里叶变换(FFT)，得到向量zpk。
步骤3.3、将向量zpk存放到矩阵QSK×L的第k行，QSK×L为步骤1定义的频域数据矩阵。
步骤4、对频域的数据矩阵的每一行进行频域升采样
对步骤3得到的数据矩阵QSK×L的每一行统一做如下8倍频域升采样处理：
步骤4.1、取出步骤3中的频域数据矩阵QSK×L的第k行向量，记为fk，k＝1,2,…,K，K为步骤1定义的慢时刻数。
步骤4.2、从向量fk的L/2+1位置开始插入7×L个零，得到向量zk，zk＝[fk(1),fk(2),…,fk(L/2),01×7L,fk(L/2+1),…,fk(L)]，fk(1)为向量fk中的第1个元素，fk(2)为向量fk中的第2个元素，fk(L/2)为向量fk中的第L/2个元素，01×7L为1行、7×L列的零向量，fk(L/2+1)为向量fk中的第L/2+1个元素，fk(L)为向量fk中的第L个元素，L为步骤1定义的雷达回波数据矩阵距离向点数。
步骤4.3、将向量zk存放到矩阵SSK×P的第k行，SSK×P为步骤1定义的升采样数据矩阵。
步骤5、计算初始搜索方向
步骤5.1、采用公式
Figure GDA0003534881160000051
计算第k个匀速直线APC，记为
Figure GDA0003534881160000052
Pt0为步骤1定义的雷达平台在零时刻的位置向量，V为步骤1定义的雷达平台速度向量，ts(k)为向量ts的第k个元素，k＝1,2,…,K，K为步骤1定义的慢时刻数，ts为步骤1定义的慢时刻向量。
步骤5.2、采用公式
Figure GDA0003534881160000053
计算像素点网格ΩM×N中第m行n列像素点的位置向量，记为
Figure GDA0003534881160000054
Pc为步骤1定义的场景中心位置向量，m＝1,2,…,M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1,2,…,N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数，dx为步骤1定义的像素点网格ΩM×N中X方向的像素点间隔，dy为步骤1定义的像素点网格ΩM×N中Y方向的像素点间隔。
步骤5.3、采用公式
Figure GDA0003534881160000055
计算第k个APC误差向量，记为
Figure GDA0003534881160000056
Figure GDA0003534881160000057
为步骤1定义的初始APC误差向量，
Figure GDA0003534881160000058
为
Figure GDA0003534881160000059
中第2(k-1)+1个元素，
Figure GDA00035348811600000510
为
Figure GDA00035348811600000511
中第2(k-1)+2个元素，k＝1,2,…,K，K为步骤1定义的慢时刻数；
Figure GDA0003534881160000061
为步骤5.1计算出的第k个匀速直线APC，
Figure GDA0003534881160000062
为步骤5.2计算出的像素点网格ΩM×N中第m行n列像素点的位置向量，m＝1,2,…,M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1,2,…,N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数。
步骤5.4、采用公式
Figure GDA0003534881160000063
计算第k个APC对应的在像素点网格ΩM×N中第m行n列像素点回波的索引值，记为
Figure GDA0003534881160000064
gx为距离向波数，gy为方位向波数，Δg为波数间隔，m＝1,2,…,M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1,2,…,N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数。
步骤5.5、利用
Figure GDA0003534881160000065
在SSK×P第k行数据里找到
Figure GDA0003534881160000066
对应的回波信号值，记为
Figure GDA0003534881160000067
Figure GDA0003534881160000068
为步骤5.4计算出的第k个APC对应的在像素点网格ΩM×N中第m行n列像素点回波的索引值，SSK×P为步骤3计算出的升采样数据矩阵，m＝1,2,…,M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1,2,…,N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数，k＝1,2,…,K，K为步骤1定义的慢时刻数，C为步骤1定义的光速。
步骤5.6、采用公式
Figure GDA0003534881160000069
计算像素点网格ΩM×N中第m行n列像素点在频域的后向投影值，
Figure GDA00035348811600000610
为步骤5.5计算出的
Figure GDA00035348811600000611
对应的回波信号值，
Figure GDA00035348811600000612
为步骤5.4计算出的第k个APC对应的在像素点网格ΩM×N中第m行n列像素点回波的索引值，G为步骤1定义的波数矢量，m＝1,2,…,M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1,2,…,N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数，k＝1,2,…,K，K为步骤1定义的慢时刻数，
Figure GDA00035348811600000613
为步骤5.3计算的APC误差向量，
Figure GDA00035348811600000614
为步骤5.1计算的匀速APC向量。
步骤5.7、采用公式
Figure GDA0003534881160000071
计算像素点网格ΩM×N中第m行n列像素点在时域的后向投影值，
Figure GDA0003534881160000072
为步骤5.6计算出的像素点网格ΩM×N中第m行n列像素点在频域的后向投影值，G为步骤1定义的波数矢量，
Figure GDA0003534881160000073
为步骤5.3计算的APC误差向量，
Figure GDA0003534881160000074
为步骤5.2中计算出的第m行n列像素点的位置向量，m＝1,2,…,M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1,2,…,N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数。
步骤5.8、采用公式
Figure GDA0003534881160000075
计算
Figure GDA0003534881160000076
关于
Figure GDA0003534881160000077
的偏导数，记为
Figure GDA0003534881160000078
Figure GDA0003534881160000079
为步骤5.5计算出的
Figure GDA00035348811600000710
对应的回波信号值，
Figure GDA00035348811600000711
为步骤5.4计算出的第k个APC对应的在像素点网格ΩM×N中第m行n列像素点回波的索引值，
Figure GDA00035348811600000712
为步骤5.3计算的APC误差向量，
Figure GDA00035348811600000713
为步骤5.1计算的匀速APC向量，G为步骤1定义的波数矢量，k＝1,2,…,K，K为步骤1定义的慢时刻数，Ro为步骤1定义的参考距离，gx为距离向波数，gy为距离向波数，m＝1,2,…,M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1,2,…,N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数。
步骤5.9、采用公式
Figure GDA00035348811600000714
计算第k个慢时刻对应的中间向量，记为
Figure GDA00035348811600000715
Figure GDA00035348811600000716
为步骤5.8计算的
Figure GDA00035348811600000717
关于
Figure GDA00035348811600000718
的偏导数，
Figure GDA00035348811600000719
为步骤5.7计算的网格ΩM×N中第m行n列像素点像素点在时域的后向投影值，k＝1,2,…,K，K为步骤1定义的慢时刻数，m＝1,2,…,M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1,2,…,N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数，(·)*表示向量共轭运算，Re[·]表示取实数部分。
步骤5.10、重复步骤5.1到步骤5.9，对所有的k＝1,2,…,K，K为步骤1定义的慢时刻数，计算出
Figure GDA0003534881160000081
Figure GDA0003534881160000082
为第1个慢时刻对应的中间向量，
Figure GDA0003534881160000083
为第2个慢时刻对应的中间向量，
Figure GDA0003534881160000084
为第K个慢时刻对应的中间向量。
步骤5.11、采用公式
Figure GDA0003534881160000085
计算初始梯度向量，记为
Figure GDA0003534881160000086
Figure GDA0003534881160000087
为步骤5.10计算出的第1个慢时刻对应的中间向量，
Figure GDA0003534881160000088
为步骤5.10计算出的第2个慢时刻对应的中间向量，
Figure GDA0003534881160000089
为步骤5.10计算出的第K个慢时刻对应的中间向量，(·)T表示向量转置运算。
步骤5.12、采用公式
Figure GDA00035348811600000810
计算初始搜索方向，记为d0，
Figure GDA00035348811600000811
为步骤5.11计算出的初始梯度向量。
步骤6、初始化迭代变量
步骤6.1、定义当前迭代次数为q，q＝0,1,2,…,Q，Q为步骤1定义的共轭梯度算法最大迭代次数，初始化q＝0。
步骤6.2、定义第q次迭代APC误差向量为
Figure GDA00035348811600000812
并初始化为
Figure GDA00035348811600000813
Figure GDA00035348811600000814
为步骤1定义的初始APC误差向量，q为步骤6.1定义的当前迭代次数。
步骤6.3、定义第q次迭代搜索方向为dq，并初始化为dq＝d0，d0为步骤5.12计算出的初始搜索方向，q为步骤6.1定义的当前迭代次数。
步骤6.4、定义第q次迭代梯度向量为
Figure GDA00035348811600000815
并初始化为
Figure GDA00035348811600000816
Figure GDA00035348811600000817
为步骤5.11计算出的初始梯度向量，q为步骤6.1定义的当前迭代次数。
步骤7、判断迭代是否结束
如果当前迭代次数q满足q≥Q，Q为步骤1定义的共轭梯度算法最大迭代次数，则结束迭代，输出
Figure GDA00035348811600000818
Figure GDA00035348811600000819
为步骤6.2定义的第q次迭代APC误差向量；如果当前迭代次数q满足q<Q，则继续执行步骤7。
步骤8、计算最佳搜索步长
利用Armijo算法计算第q次迭代最佳搜索步长，记为λq，λq的大小为1行、2×K列。
步骤9、计算第q+1次迭代APC误差向量
采用公式
Figure GDA0003534881160000091
计算第q+1次迭代APC误差向量，记为
Figure GDA0003534881160000092
λq为步骤8计算出的第q次迭代最佳搜索步长，dq为步骤6.3计算出的第q次迭代搜索方向，
Figure GDA0003534881160000093
为步骤6.2定义的第q次迭代APC误差向量，q为步骤6.1定义的当前迭代次数。
步骤10、计算第q+1次迭代搜索方向
步骤10.1、计算第q+1次迭代梯度向量，具体方法是：
步骤10.1.1、采用公式
Figure GDA0003534881160000094
计算第q+1次迭代、第k个匀速直线APC，记为
Figure GDA0003534881160000095
Pt0为步骤1定义的雷达平台在零时刻的位置向量，V为步骤1定义的雷达平台速度向量，ts(k)为向量ts的第k个元素，k＝1,2,…,K，K为步骤1定义的慢时刻数，ts为步骤1定义的慢时刻向量，q为步骤6.1定义的当前迭代次数。
步骤10.1.2、采用公式
Figure GDA0003534881160000096
计算第q+1次迭代、像素点网格ΩM×N中第m行n列像素点的位置向量，记为
Figure GDA0003534881160000097
Pc为步骤1定义的场景中心位置向量，m＝1,2,…,M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1,2,…,N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数，dx为步骤1定义的像素点网格ΩM×N中X方向的像素点间隔，dy为步骤1定义的像素点网格ΩM×N中Y方向的像素点间隔，q为步骤6.1定义的当前迭代次数。
步骤10.1.3、采用公式
Figure GDA0003534881160000098
Figure GDA0003534881160000099
计算第q+1次迭代、第k个APC误差向量，记为
Figure GDA00035348811600000910
Figure GDA00035348811600000911
为步骤9计算出的第q+1次迭代APC误差向量，
Figure GDA00035348811600000912
为
Figure GDA00035348811600000913
中第2(k-1)+1个元素，
Figure GDA00035348811600000914
为
Figure GDA00035348811600000915
中第2(k-1)+2个元素，k＝1,2,…,K，K为步骤1定义的慢时刻数；
Figure GDA00035348811600000916
为步骤10.1.1计算出的第q+1次迭代、第k个匀速直线APC，
Figure GDA00035348811600000917
为步骤10.1.2计算出的第q+1次迭代、像素点网格ΩM×N中第m行n列像素点的位置向量，m＝1,2,…,M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1,2,…,N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数，q为步骤6.1定义的当前迭代次数。
步骤10.1.4、采用公式
Figure GDA0003534881160000101
计算第q+1次迭代、计算第k个APC对应的在像素点网格ΩM×N中第m行n列像素点回波的索引值，记为
Figure GDA0003534881160000102
gx为距离向波数，gy为方位向波数，Δg为波数间隔，m＝1,2,…,M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1,2,…,N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数，q为步骤6.1定义的当前迭代次数。
步骤10.1.5、利用第q+1次迭代、第k个APC与像素点网格ΩM×N中第m行n列像素点的回波的索引值
Figure GDA0003534881160000103
在SSK×P第k行数据里找到
Figure GDA0003534881160000104
对应的回波信号值，记为
Figure GDA0003534881160000105
Figure GDA0003534881160000106
为步骤10.1.4计算出的第q+1次迭代、第k个APC对应的在像素点网格ΩM×N中第m行n列像素点回波的索引值，SSK×P为步骤3计算出的升采样数据矩阵，m＝1,2,…,M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1,2,…,N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数，k＝1,2,…,K，K为步骤1定义的慢时刻数，C为步骤1定义的光速，q为步骤6.1定义的当前迭代次数。
步骤10.1.6、采用公式
Figure GDA0003534881160000107
计算像素点网格ΩM×N中第m行n列像素点在频域的后向投影值，
Figure GDA0003534881160000108
为步骤10.1.5计算出的
Figure GDA0003534881160000109
对应的回波信号值，
Figure GDA00035348811600001010
为步骤10.1.4计算出的第k个APC对应的在像素点网格ΩM×N中第m行n列像素点回波的索引值，G为步骤1定义的波数矢量，m＝1,2,…,M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1,2,…,N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数，k＝1,2,…,K，K为步骤1定义的慢时刻数，
Figure GDA00035348811600001011
为步骤10.1.3计算的APC误差向量，
Figure GDA0003534881160000111
为步骤10.1.1计算的匀速APC向量，q为步骤6.1定义的当前迭代次数。
步骤10.1.7、采用公式
Figure GDA0003534881160000112
计算像素点网格ΩM×N中第m行n列像素点在时域的后向投影值，
Figure GDA0003534881160000113
为步骤10.1.6计算出的像素点网格ΩM×N中第m行n列像素点在频域的后向投影值，
Figure GDA0003534881160000114
为步骤10.1.3计算的APC误差向量，
Figure GDA0003534881160000115
为步骤10.1.2中计算出的第m行n列像素点的位置向量，m＝1,2,…,M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1,2,…,N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数，G为步骤1定义的波数矢量，q为步骤6.1定义的当前迭代次数。
步骤10.1.8、采用公式
Figure GDA0003534881160000116
计算
Figure GDA0003534881160000117
关于
Figure GDA0003534881160000118
的偏导数，记为
Figure GDA0003534881160000119
Figure GDA00035348811600001110
为步骤10.1.5计算出的
Figure GDA00035348811600001111
对应的回波信号值，
Figure GDA00035348811600001112
为步骤10.1.4计算出的第k个APC对应的在像素点网格ΩM×N中第m行n列像素点回波的索引值，
Figure GDA00035348811600001113
为步骤10.1.3计算的APC误差向量，
Figure GDA00035348811600001114
为步骤10.1.1计算的匀速APC向量，G为步骤1定义的波数矢量，k＝1,2,…,K，K为步骤1定义的慢时刻数，Ro为步骤1定义的参考距离，gx为距离向波数，gy为方位向波数，m＝1,2,…,M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1,2,…,N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数，q为步骤6.1定义的当前迭代次数。
步骤10.1.9、
Figure GDA00035348811600001115
计算第k个慢时刻对应的中间向量，记为
Figure GDA00035348811600001116
Figure GDA00035348811600001117
为步骤10.1.8计算的
Figure GDA00035348811600001118
关于
Figure GDA00035348811600001119
的偏导数，
Figure GDA00035348811600001120
为步骤10.1.7计算的网格ΩM×N中第m行n列像素点像素点在时域的后向投影值，k＝1,2,…,K，K为步骤1定义的慢时刻数，m＝1,2,…,M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1,2,…,N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数，(·)*表示向量共轭运算，Re[·]表示取实数部分，q为步骤6.1定义的当前迭代次数。
步骤10.1.10、重复步骤10.1.1到步骤10.1.9，对所有的k＝1,2,…,K，K为步骤1定义的慢时刻数，计算出
Figure GDA0003534881160000121
Figure GDA0003534881160000122
为第1个慢时刻对应的第q+1次迭代中间向量，
Figure GDA0003534881160000123
为第2个慢时刻对应的第q+1次迭代中间向量，
Figure GDA0003534881160000124
为第K个慢时刻对应的第q+1次迭代中间向量，q为步骤6.1定义的当前迭代次数。
步骤10.1.11、采用公式
Figure GDA0003534881160000125
计算第q+1次迭代梯度向量，记为
Figure GDA0003534881160000126
Figure GDA0003534881160000127
为步骤10.1.10计算出的第1个慢时刻对应的第q+1次迭代中间向量，
Figure GDA0003534881160000128
为步骤10.1.10计算出的第2个慢时刻对应的第q+1次迭代中间向量，
Figure GDA0003534881160000129
为步骤10.1.10计算出的第K个慢时刻对应的第q+1次迭代中间向量，(·)T表示向量转置运算，q为步骤6.1定义的当前迭代次数。
步骤10.2、采用公式
Figure GDA00035348811600001210
计算中间参数，
Figure GDA00035348811600001211
为步骤10.1.11计算的第q+1次迭代梯度向量，
Figure GDA00035348811600001212
为步骤6.4定义的第q次迭代梯度向量，q为步骤6.1定义的当前迭代次数，‖·‖表示求向量二范数运算。
步骤10.3、采用公式
Figure GDA00035348811600001213
计算第q+1次迭代搜索方向，记为dq+1，
Figure GDA00035348811600001214
为步骤10.1.11计算出的第q+1次迭代梯度向量，βq为步骤10.2计算出的中间参数，dq为步骤6.3定义的第q次迭代搜索方向，q为步骤6.1定义的当前迭代次数。
步骤11、更新迭代次数，进入下一次迭代
采用公式q←q+1更新下一次迭代次数，运算符←表示赋值操作，重复步骤7～步骤11，直到迭代结束。
本发明的创新点在于：基于频域图像锐度最优进行合成孔径雷达天线相位中心误差估计，从而进行高精度运动误差补偿。本发明首先使用匀速直线天线相位中心采用频域后向投影成像算法进行粗聚焦成像，然后以图像锐度最优为准则，并利用频域成像算法的高效成像特点，迭代估计天线相位中心误差，迭代结束后获得天线相位中心误差估计，最后利用匀速直线天线相位中心和估计出的天线相位中心误差进行高精度后向投影成像。本发明方法与现有的时域自聚焦BP算法相比，简化了自聚焦算法的公式推导过程，并且由于频域BP算法的特点，大大地提高了成像速度和精度，且降低了内存消耗。
本发明的优点：与时域自聚焦BP算法相比，本发明由于首先在频域进行成像处理，且无需要计算散射点到阵元的距离史，本发明所需运算时间比时域的自聚焦BP算法大大减少，因此更适用于大场景、长孔径、高精度SAR成像。除此之外，在成像过程中未采用任何近似，考虑了相位误差的空变性，通过有效估计天线相位中心误差，对场景中每个像素点都能进行高精度的相位补偿。
附图说明
图1是本发明流程图。
具体实施方式
本发明主要采用计算机仿真的方法进行验证，所有步骤、结论都在cup为InterI7-8700K 3.7GHz，MATLAB版本为R2017b上验证正确。具体实施步骤如下：
步骤1、相关参数的初始化
初始化的参数均为已知，且初始化的参数如下：光速为C＝3×108m/s；雷达发射线性调频信号，载频为ω＝6×π×109rad/s；雷达发射脉冲的带宽为B＝300×106Hz；雷达发射脉冲的时宽为Tp＝1×10-6s；雷达脉冲重复周期为T＝0.002s；雷达回波距离向采样频率为fs＝390×106Hz；雷达回波数据矩阵为SK×L；雷达回波数据矩阵SK×L的方位向点数和距离向点数分别为K＝500和L＝2048(K和L均为正整数)，K也称为慢时刻数；慢时刻向量为ts＝[-500/2,1-500/2,…,500/2-1]×0.002s；频域数据矩阵为QSK×L，QSK×L的大小为K行、L列，K为雷达回波数据矩阵的方位向点数，L为雷达回波数据的距离向点数；升采样数据矩阵为SSK×P，SSK×P的大小为K行、8×L列，K为雷达回波数据矩阵的方位向点数，L为雷达回波数据的距离向点数；几何坐标系为三维笛卡尔坐标系O-XY；雷达平台速度向量为V＝[0,50,0]m/s，V的大小为1行3列；雷达平台在零时刻的位置向量为Pt0＝[0,0,4000]m，Pt0的大小为1行3列；雷达参考距离为R0；O-XY平面内的矩形场景为Θ；将Θ离散化为像素点网格，记为ΩM×N；像素点网格ΩM×N中X方向和Y方向的像素点点数分别为M＝25和N＝1000；像素点网格ΩM×N中X方向和Y方向的像素点间隔分别为dx＝0.2m和dy＝0.2m；场景中心位置向量为Pc＝[3000,0,0]m，Pc的大小为1行3列；共轭梯度算法最大迭代次数为Q＝100，Q为正整数；初始APC误差向量为
Figure GDA0003534881160000141
Figure GDA0003534881160000142
01×2K为1行、2×K列的零矩阵，K为雷达回波数据矩阵SK×L方位向点数；G＝[gx,gy]为波数矢量，gx为距离向波数，gy为方位向波数，离散化形式为gi＝g1+(i-1)Δg,i＝1,...,M，g1＝2πf1/c，Δg＝2πΔf/c，Δg为波数间隔，f1和Δf分别是信号的起始频率和频率间隔。
步骤2、对雷达回波数据矩阵的每一行进行脉冲压缩
取出步骤1中所有雷达回波数据S500×2048，采用脉冲压缩方法对S500×2048的每一行进行脉冲压缩，得到脉冲压缩后的数据矩阵PS500×2048。
步骤3、将脉冲压缩后的数据矩阵转至频域
对步骤2得到的脉冲压缩后的数据矩阵PS500×2048做快速傅里叶变换(FFT)：
步骤3.1、取出步骤2中脉冲压缩后的数据矩阵PS500×2048的第k行向量，记为ppk，k＝1,2,…,500，其为步骤1定义的慢时刻数。
步骤3.2、对向量ppk做快速傅里叶变换(FFT)，得到向量zpk。
步骤3.3、将向量zpk存放到矩阵QS500×2048的第k行，QS500×2048为步骤1定义的频域数据矩阵。
步骤4、对频域的数据矩阵的每一行进行频域升采样
对步骤3得到的数据矩阵QS500×2048的每一行统一做如下8倍频域升采样处理：
步骤4.1、取出步骤3中的频域数据矩阵QS500×2048的第k行向量，记为fk，k＝1,2,…,500，其为步骤1定义的慢时刻数。
步骤4.2、从向量fk的2048/2+1位置开始插入7×2048个零，得到向量zk，zk＝[fk(1),fk(2),…,fk(2048/2),01×16384,fk(2048/2+1),…,fk(2048)]，fk(1)为向量fk中的第1个元素，fk(2)为向量fk中的第2个元素，fk(2048/2)为向量fk中的第2048/2个元素，01×16384为1行、16384列的零向量，fk(2048/2+1)为向量fk中的第2048/2+1个元素，fk(2048)为向量fk中的第2048个元素，2048为步骤1定义的雷达回波数据矩阵距离向点数。
步骤4.3、将向量zk存放到矩阵SS500×16384的第k行，SS500×16384为步骤1定义的升采样数据矩阵。
步骤5、计算初始搜索方向
步骤5.1、采用公式
Figure GDA0003534881160000151
计算第k个匀速直线APC，记为
Figure GDA0003534881160000152
Pt0＝[0,0,4000]m为步骤1定义的雷达平台在零时刻的位置向量，V＝[0,50,0]m/s为步骤1定义的雷达平台速度向量，ts(k)为向量ts＝[-500/2,1-500/2,…,500/2-1]×0.002s的第k个元素，k＝1,2,…,500，其为步骤1定义的慢时刻数，ts为步骤1定义的慢时刻向量。
步骤5.2、采用公式
Figure GDA0003534881160000153
计算像素点网格Ω25×1000中第m行n列像素点的位置向量，记为
Figure GDA0003534881160000154
Pc＝[3000,0,0]m，其为步骤1定义的场景中心位置向量，m＝1,2,…,25为步骤1定义的像素点网格Ω25×1000中X方向的像素点点数，n＝1,2,…,1000为步骤1定义的像素点网格Ω25×1000中Y方向的像素点点数，dx＝0.2m为步骤1定义的像素点网格Ω25×1000中X方向的像素点间隔，dy＝0.2m为步骤1定义的像素点网格Ω25×1000中Y方向的像素点间隔。
步骤5.3、采用公式
Figure GDA0003534881160000155
计算第k个APC误差向量，记为
Figure GDA0003534881160000156
Figure GDA0003534881160000157
为步骤1定义的初始APC误差向量，
Figure GDA0003534881160000158
为
Figure GDA0003534881160000159
中第2(k-1)+1个元素，
Figure GDA00035348811600001510
为
Figure GDA00035348811600001511
中第2(k-1)+2个元素，k＝1,2,…,500为步骤1定义的慢时刻数；
Figure GDA00035348811600001512
为步骤5.1计算出的第k个匀速直线APC，
Figure GDA00035348811600001513
为步骤5.2计算出的像素点网格Ω25×1000中第m行n列像素点的位置向量，m＝1,2,…,25为步骤1定义的像素点网格Ω25×1000中X方向的像素点点数，n＝1,2,…,1000为步骤1定义的像素点网格Ω25×1000中Y方向的像素点点数。
步骤5.4、采用公式
Figure GDA0003534881160000161
计算第k个APC对应的在像素点网格Ω25×1000中第m行n列像素点回波的索引值，记为
Figure GDA0003534881160000162
gx为距离向波数，gy为方位向波数，Δg为波数间隔，m＝1,2,…,25为步骤1定义的像素点网格Ω25×1000中X方向的像素点点数，n＝1,2,…,1000为步骤1定义的像素点网格Ω25×1000中Y方向的像素点点数。
步骤5.5、利用
Figure GDA0003534881160000163
在SS500×16384第k行数据里找到
Figure GDA0003534881160000164
对应的回波信号值，记为
Figure GDA0003534881160000165
Figure GDA0003534881160000166
为步骤5.4计算出的第k个APC对应的在像素点网格Ω25×1000中第m行n列像素点回波的索引值，SS500×16384为步骤3计算出的升采样数据矩阵，m＝1,2,…,25为步骤1定义的像素点网格Ω25×1000中X方向的像素点点数，n＝1,2,…,1000为步骤1定义的像素点网格Ω25×1000中Y方向的像素点点数，k＝1,2,…,500为步骤1定义的慢时刻数，C＝3×108m/s为步骤1定义的光速。
步骤5.6、采用公式
Figure GDA0003534881160000167
计算像素点网格Ω25×1000中第m行n列像素点在频域的后向投影值，
Figure GDA0003534881160000168
为步骤5.5计算出的
Figure GDA0003534881160000169
对应的回波信号值，
Figure GDA00035348811600001610
为步骤5.4计算出的第k个APC对应的在像素点网格Ω25×1000中第m行n列像素点回波的索引值，G为步骤1定义的波数矢量，m＝1,2,…,25为步骤1定义的像素点网格Ω25×1000中X方向的像素点点数，n＝1,2,…,1000为步骤1定义的像素点网格Ω25×1000中Y方向的像素点点数，k＝1,2,…,500为步骤1定义的慢时刻数，
Figure GDA00035348811600001611
为步骤5.3计算的APC误差向量，
Figure GDA00035348811600001612
为步骤5.1计算的匀速APC向量。
步骤5.7、采用公式
Figure GDA00035348811600001613
计算像素点网格Ω25×1000中第m行n列像素点在时域的后向投影值，
Figure GDA00035348811600001614
为步骤5.6计算出的像素点网格Ω25×1000中第m行n列像素点在频域的后向投影值，G为步骤1定义的波数矢量，
Figure GDA00035348811600001615
为步骤5.3计算的APC误差向量，
Figure GDA00035348811600001616
为步骤5.2中计算出的第m行n列像素点的位置向量，m＝1,2,…,25为步骤1定义的像素点网格Ω25×1000中X方向的像素点点数，n＝1,2,…,1000为步骤1定义的像素点网格Ω25×1000中Y方向的像素点点数。
步骤5.8、采用公式
Figure GDA0003534881160000171
计算
Figure GDA0003534881160000172
关于
Figure GDA0003534881160000173
的偏导数，记为
Figure GDA0003534881160000174
Figure GDA0003534881160000175
为步骤5.5计算出的
Figure GDA0003534881160000176
对应的回波信号值，
Figure GDA0003534881160000177
为步骤5.4计算出的第k个APC对应的在像素点网格Ω25×1000中第m行n列像素点回波的索引值，
Figure GDA0003534881160000178
为步骤5.3计算的APC误差向量，
Figure GDA0003534881160000179
为步骤5.1计算的匀速APC向量，G为步骤1定义的波数矢量，k＝1,2,…,500为步骤1定义的慢时刻数，Ro为步骤1定义的参考距离，gx为距离向波数，gy为距离向波数，m＝1,2,…,25为步骤1定义的像素点网格Ω25×1000中X方向的像素点点数，n＝1,2,…,1000为步骤1定义的像素点网格Ω25×1000中Y方向的像素点点数。
步骤5.9、采用公式
Figure GDA00035348811600001710
计算第k个慢时刻对应的中间向量，记为
Figure GDA00035348811600001711
Figure GDA00035348811600001712
为步骤5.8计算的
Figure GDA00035348811600001713
关于
Figure GDA00035348811600001714
的偏导数，
Figure GDA00035348811600001715
为步骤5.7计算的网格Ω25×1000中第m行n列像素点像素点在时域的后向投影值，k＝1,2,…,500为步骤1定义的慢时刻数，m＝1,2,…,25为步骤1定义的像素点网格Ω25×1000中X方向的像素点点数，n＝1,2,…,1000为步骤1定义的像素点网格Ω25×1000中Y方向的像素点点数，(·)*表示向量共轭运算，Re[·]表示取实数部分。
步骤5.10、重复步骤5.1到步骤5.9，对所有的k＝1,2,…,500，其为步骤1定义的慢时刻数，计算出
Figure GDA00035348811600001716
Figure GDA00035348811600001717
为第1个慢时刻对应的中间向量，
Figure GDA00035348811600001718
为第2个慢时刻对应的中间向量，
Figure GDA00035348811600001719
为第500个慢时刻对应的中间向量。
步骤5.11、采用公式
Figure GDA0003534881160000181
计算初始梯度向量，记为
Figure GDA0003534881160000182
Figure GDA0003534881160000183
为步骤5.10计算出的第1个慢时刻对应的中间向量，
Figure GDA0003534881160000184
为步骤5.10计算出的第2个慢时刻对应的中间向量，
Figure GDA0003534881160000185
为步骤5.10计算出的第500个慢时刻对应的中间向量，(·)T表示向量转置运算。
步骤5.12、采用公式
Figure GDA0003534881160000186
计算初始搜索方向，记为d0，
Figure GDA0003534881160000187
为步骤5.11计算出的初始梯度向量。
步骤6、初始化迭代变量
步骤6.1、定义当前迭代次数为q，q＝0,1,2,…,100，100为步骤1定义的共轭梯度算法最大迭代次数，初始化q＝0。
步骤6.2、定义第q次迭代APC误差向量为
Figure GDA0003534881160000188
并初始化为
Figure GDA0003534881160000189
Figure GDA00035348811600001810
为步骤1定义的初始APC误差向量，q为步骤6.1定义的当前迭代次数。
步骤6.3、定义第q次迭代搜索方向为dq，并初始化为dq＝d0，d0为步骤5.12计算出的初始搜索方向，q为步骤6.1定义的当前迭代次数。
步骤6.4、定义第q次迭代梯度向量为
Figure GDA00035348811600001811
并初始化为
Figure GDA00035348811600001812
Figure GDA00035348811600001813
为步骤5.11计算出的初始梯度向量，q为步骤6.1定义的当前迭代次数。
步骤7、判断迭代是否结束
如果当前迭代次数q满足q≥100，则结束迭代，输出
Figure GDA00035348811600001814
Figure GDA00035348811600001815
为步骤6.2定义的第q次迭代APC误差向量；如果当前迭代次数q满足q<100，则继续执行步骤7。
步骤8、计算最佳搜索步长
利用Armijo算法计算第q次迭代最佳搜索步长，记为λq，λq的大小为1行、2×500列。
步骤9、计算第q+1次迭代APC误差向量
采用公式
Figure GDA00035348811600001816
计算第q+1次迭代APC误差向量，记为
Figure GDA00035348811600001817
λq为步骤8计算出的第q次迭代最佳搜索步长，dq为步骤6.3计算出的第q次迭代搜索方向，
Figure GDA00035348811600001818
为步骤6.2定义的第q次迭代APC误差向量，q为步骤6.1定义的当前迭代次数。
步骤10、计算第q+1次迭代搜索方向
步骤10.1、计算第q+1次迭代梯度向量，具体方法是：
步骤10.1.1、采用公式
Figure GDA0003534881160000191
计算第q+1次迭代、第k个匀速直线APC，记为
Figure GDA0003534881160000192
Pt0＝[0,0,4000]m，其为步骤1定义的雷达平台在零时刻的位置向量，V＝[0,50,0]m/s为步骤1定义的雷达平台速度向量，ts(k)为向量ts＝[-500/2,1-500/2,…,500/2-1]×0.002s的第k个元素，k＝1,2,…,500，其为步骤1定义的慢时刻数，ts为步骤1定义的慢时刻向量，q为步骤6.1定义的当前迭代次数。
步骤10.1.2、采用公式
Figure GDA0003534881160000193
Figure GDA0003534881160000194
计算第q+1次迭代、像素点网格Ω25×1000中第m行n列像素点的位置向量，记为
Figure GDA0003534881160000195
Pc＝[3000,0,0]m为步骤1定义的场景中心位置向量，m＝1,2,…,25，其为步骤1定义的像素点网格Ω25×1000中X方向的像素点点数，n＝1,2,…,1000，其为步骤1定义的像素点网格Ω25×1000中Y方向的像素点点数，dx＝0.2m为步骤1定义的像素点网格Ω25×1000中X方向的像素点间隔，dy＝0.2m为步骤1定义的像素点网格Ω25×1000中Y方向的像素点间隔，q为步骤6.1定义的当前迭代次数。
步骤10.1.3、采用公式
Figure GDA0003534881160000196
Figure GDA0003534881160000197
计算第q+1次迭代、第k个APC误差向量，记为
Figure GDA0003534881160000198
Figure GDA0003534881160000199
为步骤9计算出的第q+1次迭代APC误差向量，
Figure GDA00035348811600001910
为
Figure GDA00035348811600001911
中第2(k-1)+1个元素，
Figure GDA00035348811600001912
为
Figure GDA00035348811600001913
中第2(k-1)+2个元素，k＝1,2,…,500，其为步骤1定义的慢时刻数；
Figure GDA00035348811600001914
为步骤10.1.1计算出的第q+1次迭代、第k个匀速直线APC，
Figure GDA00035348811600001915
为步骤10.1.2计算出的第q+1次迭代、像素点网格Ω25×1000中第m行n列像素点的位置向量，m＝1,2,…,25，其为步骤1定义的像素点网格Ω25×1000中X方向的像素点点数，n＝1,2,…,1000，其为步骤1定义的像素点网格Ω25×1000中Y方向的像素点点数，q为步骤6.1定义的当前迭代次数。
步骤10.1.4、采用公式
Figure GDA0003534881160000201
计算第q+1次迭代、计算第k个APC对应的在像素点网格Ω25×1000中第m行n列像素点回波的索引值，记为
Figure GDA0003534881160000202
gx为距离向波数，gy为方位向波数，Δg为波数间隔，m＝1,2,…,25，其为步骤1定义的像素点网格Ω25×1000中X方向的像素点点数，n＝1,2,…,1000，其为步骤1定义的像素点网格Ω25×1000中Y方向的像素点点数，q为步骤6.1定义的当前迭代次数。
步骤10.1.5、利用第q+1次迭代、第k个APC与像素点网格Ω25×1000中第m行n列像素点的回波的索引值
Figure GDA0003534881160000203
在SS500×16384第k行数据里找到
Figure GDA0003534881160000204
对应的回波信号值，记为
Figure GDA0003534881160000205
Figure GDA0003534881160000206
为步骤10.1.4计算出的第q+1次迭代、第k个APC对应的在像素点网格Ω25×1000中第m行n列像素点回波的索引值，SS500×16384为步骤3计算出的升采样数据矩阵，m＝1,2,…,25，其为步骤1定义的像素点网格Ω25×1000中X方向的像素点点数，n＝1,2,…,1000，其为步骤1定义的像素点网格Ω25×1000中Y方向的像素点点数，k＝1,2,…,500，其为步骤1定义的慢时刻数，C＝3×108m/s为步骤1定义的光速，q为步骤6.1定义的当前迭代次数。
步骤10.1.6、采用公式
Figure GDA0003534881160000207
计算像素点网格Ω25×1000中第m行n列像素点在频域的后向投影值，
Figure GDA0003534881160000208
为步骤10.1.5计算出的
Figure GDA0003534881160000209
对应的回波信号值，
Figure GDA00035348811600002010
为步骤10.1.4计算出的第k个APC对应的在像素点网格Ω25×1000中第m行n列像素点回波的索引值，G为步骤1定义的波数矢量，m＝1,2,…,25，其为步骤1定义的像素点网格Ω25×1000中X方向的像素点点数，n＝1,2,…,1000，N为步骤1定义的像素点网格Ω25×1000中Y方向的像素点点数，k＝1,2,…,500，其为步骤1定义的慢时刻数，
Figure GDA00035348811600002011
为步骤10.1.3计算的APC误差向量，
Figure GDA00035348811600002012
为步骤10.1.1计算的匀速APC向量，q为步骤6.1定义的当前迭代次数。
步骤10.1.7、采用公式
Figure GDA0003534881160000211
计算像素点网格Ω25×1000中第m行n列像素点在时域的后向投影值，
Figure GDA0003534881160000212
为步骤10.1.6计算出的像素点网格Ω25×1000中第m行n列像素点在频域的后向投影值，
Figure GDA0003534881160000213
为步骤10.1.3计算的APC误差向量，
Figure GDA0003534881160000214
为步骤10.1.2中计算出的第m行n列像素点的位置向量，m＝1,2,…,25，其为步骤1定义的像素点网格Ω25×1000中X方向的像素点点数，n＝1,2,…,1000，其为步骤1定义的像素点网格Ω25×1000中Y方向的像素点点数，G为步骤1定义的波数矢量，q为步骤6.1定义的当前迭代次数。
步骤10.1.8、采用公式
Figure GDA0003534881160000215
计算
Figure GDA0003534881160000216
关于
Figure GDA0003534881160000217
的偏导数，记为
Figure GDA0003534881160000218
Figure GDA0003534881160000219
为步骤10.1.5计算出的
Figure GDA00035348811600002110
对应的回波信号值，
Figure GDA00035348811600002111
为步骤10.1.4计算出的第k个APC对应的在像素点网格Ω25×1000中第m行n列像素点回波的索引值，
Figure GDA00035348811600002112
为步骤10.1.3计算的APC误差向量，
Figure GDA00035348811600002113
为步骤10.1.1计算的匀速APC向量，G为步骤1定义的波数矢量，k＝1,2,…,500，其为步骤1定义的慢时刻数，Ro为步骤1定义的参考距离，gx为距离向波数，gy为方位向波数，m＝1,2,…,25，其为步骤1定义的像素点网格Ω25×1000中X方向的像素点点数，n＝1,2,…,1000，其为步骤1定义的像素点网格Ω25×1000中Y方向的像素点点数，q为步骤6.1定义的当前迭代次数。
步骤10.1.9、
Figure GDA00035348811600002114
计算第k个慢时刻对应的中间向量，记为
Figure GDA00035348811600002115
Figure GDA00035348811600002116
为步骤10.1.8计算的
Figure GDA00035348811600002117
关于
Figure GDA00035348811600002118
的偏导数，
Figure GDA00035348811600002119
为步骤10.1.7计算的网格Ω25×1000中第m行n列像素点像素点在时域的后向投影值，k＝1,2,…,500，其为步骤1定义的慢时刻数，m＝1,2,…,25，其为步骤1定义的像素点网格Ω25×1000中X方向的像素点点数，n＝1,2,…,1000，其为步骤1定义的像素点网格Ω25×1000中Y方向的像素点点数，(·)*表示向量共轭运算，Re[·]表示取实数部分，q为步骤6.1定义的当前迭代次数。
步骤10.1.10、重复步骤10.1.1到步骤10.1.9，对所有的k＝1,2,…,500，其为步骤1定义的慢时刻数，计算出
Figure GDA0003534881160000221
Figure GDA0003534881160000222
为第1个慢时刻对应的第q+1次迭代中间向量，
Figure GDA0003534881160000223
为第2个慢时刻对应的第q+1次迭代中间向量，
Figure GDA0003534881160000224
为第500个慢时刻对应的第q+1次迭代中间向量，q为步骤6.1定义的当前迭代次数。
步骤10.1.11、采用公式
Figure GDA0003534881160000225
计算第q+1次迭代梯度向量，记为
Figure GDA0003534881160000226
Figure GDA0003534881160000227
为步骤10.1.10计算出的第1个慢时刻对应的第q+1次迭代中间向量，
Figure GDA0003534881160000228
为步骤10.1.10计算出的第2个慢时刻对应的第q+1次迭代中间向量，
Figure GDA0003534881160000229
为步骤10.1.10计算出的第500个慢时刻对应的第q+1次迭代中间向量，(·)T表示向量转置运算，q为步骤6.1定义的当前迭代次数。
步骤10.2、采用公式
Figure GDA00035348811600002210
计算中间参数，
Figure GDA00035348811600002211
为步骤10.1.11计算的第q+1次迭代梯度向量，
Figure GDA00035348811600002212
为步骤6.4定义的第q次迭代梯度向量，q为步骤6.1定义的当前迭代次数，‖·‖表示求向量二范数运算。
步骤10.3、采用公式
Figure GDA00035348811600002213
计算第q+1次迭代搜索方向，记为dq+1，
Figure GDA00035348811600002214
为步骤10.1.11计算出的第q+1次迭代梯度向量，βq为步骤10.2计算出的中间参数，dq为步骤6.3定义的第q次迭代搜索方向，q为步骤6.1定义的当前迭代次数。
步骤11、更新迭代次数，进入下一次迭代
采用公式q←q+1更新下一次迭代次数，运算符←表示赋值操作，重复步骤7～步骤11，直到迭代结束。
通过本发明的具体实施可以看出，相比于时域的自聚焦BP算法，本发明在成像过程中未采用任何近似，具有成像精度高的优点，且同样适用于任意成像模式和任意轨迹，实现了对场景中每个像素点都进行高精度运动补偿，从而大大提高了成像精度；由于本发明首先在频域进行成像处理，且无需要计算散射点到阵元的距离史，所以所需时间大大降低。

Claims (1)

1.一种合成孔径雷达频域BP算法的运动误差补偿方法，其特征是它包括以下步骤：
步骤1、相关参数的初始化
初始化的参数均为已知，且初始化的参数如下：光速为C；雷达发射线性调频信号，载频为ω；雷达发射脉冲的带宽为B；雷达发射脉冲的时宽为Tp；雷达脉冲重复周期为T；雷达回波距离向采样频率为fs；雷达回波数据矩阵为SK×L；雷达回波数据矩阵SK×L的方位向点数和距离向点数分别为K和L(K和L均为正整数)，K也称为慢时刻数；慢时刻向量为ts＝[-K/2，1-K/2，...，K/2-1]×T；频域数据矩阵为QSK×L，QSK×L的大小为K行、L列，K为雷达回波数据矩阵的方位向点数，L为雷达回波数据的距离向点数；升采样数据矩阵为SSK×P，SSK×P的大小为K行、8×L列，K为雷达回波数据矩阵的方位向点数，L为雷达回波数据的距离向点数；几何坐标系为三维笛卡尔坐标系O-XY；雷达平台速度向量为V，V的大小为1行3列；雷达平台在零时刻的位置向量为Pt0，Pt0的大小为1行3列；雷达参考距离为R0；O-XY平面内的矩形场景为Θ；将Θ离散化为像素点网格，记为ΩM×N；像素点网格ΩM×N中X方向和Y方向的像素点点数分别为M和N；像素点网格ΩM×N中X方向和Y方向的像素点间隔分别为dx和dy；场景中心位置向量为Pc，Pc的大小为1行3列；共轭梯度算法最大迭代次数为Q，Q为正整数；初始APC误差向量为
Figure FDA0003534881150000011
Figure FDA0003534881150000012
01×2K为1行、2×K列的零矩阵，K为雷达回波数据矩阵SK×L方位向点数；G＝[gx，gy]为波数矢量，gx为距离向波数，gy为方位向波数，离散化形式为gi＝g1+(i-1)Δg，i＝1，...，M，g1＝2πf1/c，Δg＝2πΔf/c，Δg为波数间隔，f1和Δf分别是信号的起始频率和频率间隔；
步骤2、对雷达回波数据矩阵的每一行进行脉冲压缩
取出步骤1中所有雷达回波数据SK×L，采用脉冲压缩方法对SK×L的每一行进行脉冲压缩，得到脉冲压缩后的数据矩阵PSK×L；
步骤3、将脉冲压缩后的数据矩阵转至频域
对步骤2得到的脉冲压缩后的数据矩阵PSK×L做快速傅里叶变换(FFT)：
步骤3.1、取出步骤2中脉冲压缩后的数据矩阵PSK×L的第k行向量，记为ppk，k＝1，2，...，K，K为步骤1定义的慢时刻数；
步骤3.2、对向量ppk做快速傅里叶变换(FFT)，得到向量zpk；
步骤3.3、将向量zpk存放到矩阵QSK×L的第k行，QSK×L为步骤1定义的频域数据矩阵；
步骤4、对频域的数据矩阵的每一行进行频域升采样
对步骤3得到的数据矩阵QSK×L的每一行统一做如下8倍频域升采样处理：
步骤4.1、取出步骤3中的频域数据矩阵QSK×L的第k行向量，记为fk，k＝1，2，...，K，K为步骤1定义的慢时刻数；
步骤4.2、从向量fk的L/2+1位置开始插入7×L个零，得到向量zk，zk＝[fk(1)，fk(2)，...，fk(L/2)，01×7L，fk(L/2+1)，...，fk(L)]，fk(1)为向量fk中的第1个元素，fk(2)为向量fk中的第2个元素，fk(L/2)为向量fk中的第L/2个元素，01×7L为1行、7×L列的零向量，fk(L/2+1)为向量fk中的第L/2+1个元素，fk(L)为向量fk中的第L个元素，L为步骤1定义的雷达回波数据矩阵距离向点数；
步骤4.3、将向量zk存放到矩阵SSK×P的第k行，SSK×P为步骤1定义的升采样数据矩阵；
步骤5、计算初始搜索方向
步骤5.1、采用公式
Figure FDA0003534881150000021
计算第k个匀速直线APC，记为
Figure FDA0003534881150000022
Pt0为步骤1定义的雷达平台在零时刻的位置向量，V为步骤1定义的雷达平台速度向量，ts(k)为向量ts的第k个元素，k＝1，2，…，K，K为步骤1定义的慢时刻数，ts为步骤1定义的慢时刻向量；
步骤5.2、采用公式
Figure FDA0003534881150000023
计算像素点网格ΩM×N中第m行n列像素点的位置向量，记为
Figure FDA0003534881150000024
Pc为步骤1定义的场景中心位置向量，m＝1，2，…，M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1，2，…，N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数，dx为步骤1定义的像素点网格ΩM×N中X方向的像素点间隔，dy为步骤1定义的像素点网格ΩM×N中Y方向的像素点间隔；
步骤5.3、采用公式
Figure FDA0003534881150000031
计算第k个APC误差向量，记为
Figure FDA0003534881150000032
Figure FDA0003534881150000033
为步骤1定义的初始APC误差向量，
Figure FDA0003534881150000034
为
Figure FDA0003534881150000035
中第2(k-1)+1个元素，
Figure FDA0003534881150000036
为
Figure FDA0003534881150000037
中第2(k-1)+2个元素，k＝1，2，...，K，K为步骤1定义的慢时刻数；
Figure FDA0003534881150000038
为步骤5.1计算出的第k个匀速直线APC，
Figure FDA0003534881150000039
为步骤5.2计算出的像素点网格ΩM×N中第m行n列像素点的位置向量，m＝1，2，…，M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1，2，…，N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数；
步骤5.4、采用公式
Figure FDA00035348811500000310
计算第k个APC对应的在像素点网格ΩM×N中第m行n列像素点回波的索引值，记为
Figure FDA00035348811500000311
gx为距离向波数，gy为方位向波数，Δg为波数间隔，m＝1，2，…，M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1，2，…，N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数；
步骤5.5、利用
Figure FDA00035348811500000312
在SSK×P第k行数据里找到
Figure FDA00035348811500000313
对应的回波信号值，记为
Figure FDA00035348811500000314
Figure FDA00035348811500000315
为步骤5.4计算出的第k个APC对应的在像素点网格ΩM×N中第m行n列像素点回波的索引值，SSK×P为步骤3计算出的升采样数据矩阵，m＝1，2，…，M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1，2，…，N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数，k＝1，2，…，K，K为步骤1定义的慢时刻数，C为步骤1定义的光速；
步骤5.6、采用公式
Figure FDA00035348811500000316
计算像素点网格ΩM×N中第m行n列像素点在频域的后向投影值，
Figure FDA00035348811500000317
为步骤5.5计算出的
Figure FDA00035348811500000318
对应的回波信号值，
Figure FDA00035348811500000319
为步骤5.4计算出的第k个APC对应的在像素点网格ΩM×N中第m行n列像素点回波的索引值，G为步骤1定义的波数矢量，m＝1，2，…，M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1，2，…，N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数，k＝1，2，…，K，K为步骤1定义的慢时刻数，
Figure FDA0003534881150000041
为步骤5.3计算的APC误差向量，
Figure FDA0003534881150000042
为步骤5.1计算的匀速APC向量；
步骤5.7、采用公式
Figure FDA0003534881150000043
计算像素点网格ΩM×N中第m行n列像素点在时域的后向投影值，
Figure FDA0003534881150000044
为步骤5.6计算出的像素点网格ΩM×N中第m行n列像素点在频域的后向投影值，G为步骤1定义的波数矢量，
Figure FDA0003534881150000045
为步骤5.3计算的APC误差向量，
Figure FDA0003534881150000046
为步骤5.2中计算出的第m行n列像素点的位置向量，m＝1，2，…，M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1，2，…，N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数；
步骤5.8、采用公式
Figure FDA0003534881150000047
计算
Figure FDA0003534881150000048
关于
Figure FDA0003534881150000049
的偏导数，记为
Figure FDA00035348811500000410
Figure FDA00035348811500000411
为步骤5.5计算出的
Figure FDA00035348811500000412
对应的回波信号值，
Figure FDA00035348811500000413
为步骤5.4计算出的第k个APC对应的在像素点网格ΩM×N中第m行n列像素点回波的索引值，
Figure FDA00035348811500000414
为步骤5.3计算的APC误差向量，
Figure FDA00035348811500000415
为步骤5.1计算的匀速APC向量，G为步骤1定义的波数矢量，k＝1，2，…，K，K为步骤1定义的慢时刻数，Ro为步骤1定义的参考距离，gx为距离向波数，gy为距离向波数，m＝1，2，…，M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1，2，…，N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数；
步骤5.9、采用公式
Figure FDA0003534881150000051
计算第k个慢时刻对应的中间向量，记为
Figure FDA0003534881150000052
Figure FDA0003534881150000053
为步骤5.8计算的
Figure FDA0003534881150000054
关于
Figure FDA0003534881150000055
的偏导数，
Figure FDA0003534881150000056
为步骤5.7计算的网格ΩM×N中第m行n列像素点像素点在时域的后向投影值，k＝1，2，…，K，K为步骤1定义的慢时刻数，m＝1，2，…，M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1，2，…，N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数，(·)*表示向量共轭运算，Re[·]表示取实数部分；
步骤5.10、重复步骤5.1到步骤5.9，对所有的k＝1，2，…，K，K为步骤1定义的慢时刻数，计算出
Figure FDA0003534881150000057
Figure FDA0003534881150000058
为第1个慢时刻对应的中间向量，
Figure FDA0003534881150000059
为第2个慢时刻对应的中间向量，
Figure FDA00035348811500000510
为第K个慢时刻对应的中间向量；
步骤5.11、采用公式
Figure FDA00035348811500000511
计算初始梯度向量，记为
Figure FDA00035348811500000512
Figure FDA00035348811500000513
为步骤5.10计算出的第1个慢时刻对应的中间向量，
Figure FDA00035348811500000514
为步骤5.10计算出的第2个慢时刻对应的中间向量，
Figure FDA00035348811500000515
为步骤5.10计算出的第K个慢时刻对应的中间向量，(·)T表示向量转置运算；
步骤5.12、采用公式
Figure FDA00035348811500000516
计算初始搜索方向，记为d0，
Figure FDA00035348811500000517
为步骤5.11计算出的初始梯度向量；
步骤6、初始化迭代变量
步骤6.1、定义当前迭代次数为q，q＝0，1，2，…，Q，Q为步骤1定义的共轭梯度算法最大迭代次数，初始化q＝0；
步骤6.2、定义第q次迭代APC误差向量为
Figure FDA00035348811500000518
并初始化为
Figure FDA00035348811500000519
Figure FDA00035348811500000520
为步骤1定义的初始APC误差向量，q为步骤6.1定义的当前迭代次数；
步骤6.3、定义第q次迭代搜索方向为dq，并初始化为dq＝d0，d0为步骤5.12计算出的初始搜索方向，q为步骤6.1定义的当前迭代次数；
步骤6.4、定义第q次迭代梯度向量为
Figure FDA0003534881150000061
并初始化为
Figure FDA0003534881150000062
Figure FDA0003534881150000063
为步骤5.11计算出的初始梯度向量，q为步骤6.1定义的当前迭代次数；
步骤7、判断迭代是否结束
如果当前迭代次数q满足q≥Q，Q为步骤1定义的共轭梯度算法最大迭代次数，则结束迭代，输出
Figure FDA0003534881150000064
Figure FDA0003534881150000065
为步骤6.2定义的第q次迭代APC误差向量；如果当前迭代次数q满足q＜Q，则继续执行步骤7；
步骤8、计算最佳搜索步长
利用Armijo算法计算第q次迭代最佳搜索步长，记为λq，λq的大小为1行、2×K列；
步骤9、计算第q+1次迭代APC误差向量
采用公式
Figure FDA0003534881150000066
计算第q+1次迭代APC误差向量，记为
Figure FDA0003534881150000067
λq为步骤8计算出的第q次迭代最佳搜索步长，dq为步骤6.3计算出的第q次迭代搜索方向，
Figure FDA0003534881150000068
为步骤6.2定义的第q次迭代APC误差向量，q为步骤6.1定义的当前迭代次数；
步骤10、计算第q+1次迭代搜索方向
步骤10.1、计算第q+1次迭代梯度向量，具体方法是：
步骤10.1.1、采用公式
Figure FDA0003534881150000069
计算第q+1次迭代、第k个匀速直线APC，记为
Figure FDA00035348811500000610
Pt0为步骤1定义的雷达平台在零时刻的位置向量，V为步骤1定义的雷达平台速度向量，ts(k)为向量ts的第k个元素，k＝1，2，…，K，K为步骤1定义的慢时刻数，ts为步骤1定义的慢时刻向量，q为步骤6.1定义的当前迭代次数；
步骤10.1.2、采用公式
Figure FDA00035348811500000611
计算第q+1次迭代、像素点网格ΩM×N中第m行n列像素点的位置向量，记为
Figure FDA00035348811500000612
Pc为步骤1定义的场景中心位置向量，m＝1，2，…，M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1，2，…，N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数，dx为步骤1定义的像素点网格ΩM×N中X方向的像素点间隔，dy为步骤1定义的像素点网格ΩM×N中Y方向的像素点间隔，q为步骤6.1定义的当前迭代次数；
步骤10.1.3、采用公式
Figure FDA0003534881150000071
计算第q+1次迭代、第k个APC误差向量，记为
Figure FDA0003534881150000072
Figure FDA0003534881150000073
为步骤9计算出的第q+1次迭代APC误差向量，
Figure FDA0003534881150000074
为
Figure FDA0003534881150000075
中第2(k-1)+1个元素，
Figure FDA0003534881150000076
为
Figure FDA0003534881150000077
中第2(k-1)+2个元素，k＝1，2，...，K，K为步骤1定义的慢时刻数；
Figure FDA0003534881150000078
为步骤10.1.1计算出的第q+1次迭代、第k个匀速直线APC，
Figure FDA0003534881150000079
为步骤10.1.2计算出的第q+1次迭代、像素点网格ΩM×N中第m行n列像素点的位置向量，m＝1，2，…，M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1，2，…，N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数，q为步骤6.1定义的当前迭代次数；
步骤10.1.4、采用公式
Figure FDA00035348811500000710
计算第q+1次迭代、计算第k个APC对应的在像素点网格ΩM×N中第m行n列像素点回波的索引值，记为
Figure FDA00035348811500000711
gx为距离向波数，gy为方位向波数，Δg为波数间隔，m＝1，2，…，M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1，2，…，N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数，q为步骤6.1定义的当前迭代次数；
步骤10.1.5、利用第q+1次迭代、第k个APC与像素点网格ΩM×N中第m行n列像素点的回波的索引值
Figure FDA00035348811500000712
在SSK×P第k行数据里找到
Figure FDA00035348811500000713
对应的回波信号值，记为
Figure FDA00035348811500000714
Figure FDA00035348811500000715
为步骤10.1.4计算出的第q+1次迭代、第k个APC对应的在像素点网格ΩM×N中第m行n列像素点回波的索引值，SSK×P为步骤3计算出的升采样数据矩阵，m＝1，2，…，M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1，2，…，N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数，k＝1，2，…，K，K为步骤1定义的慢时刻数，C为步骤1定义的光速，q为步骤6.1定义的当前迭代次数；
步骤10.1.6、采用公式
Figure FDA0003534881150000081
计算像素点网格ΩM×N中第m行n列像素点在频域的后向投影值，
Figure FDA0003534881150000082
为步骤10.1.5计算出的
Figure FDA0003534881150000083
对应的回波信号值，
Figure FDA0003534881150000084
为步骤10.1.4计算出的第k个APC对应的在像素点网格ΩM×N中第m行n列像素点回波的索引值，G为步骤1定义的波数矢量，m＝1，2，…，M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1，2，…，N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数，k＝1，2，…，K，K为步骤1定义的慢时刻数，
Figure FDA0003534881150000085
为步骤10.1.3计算的APC误差向量，
Figure FDA0003534881150000086
为步骤10.1.1计算的匀速APC向量，q为步骤6.1定义的当前迭代次数；
步骤10.1.7、采用公式
Figure FDA0003534881150000087
计算像素点网格ΩM×N中第m行n列像素点在时域的后向投影值，
Figure FDA0003534881150000088
为步骤10.1.6计算出的像素点网格ΩM×N中第m行n列像素点在频域的后向投影值，
Figure FDA0003534881150000089
为步骤10.1.3计算的APC误差向量，
Figure FDA00035348811500000810
为步骤10.1.2中计算出的第m行n列像素点的位置向量，m＝1，2，…，M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1，2，…，N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数，G为步骤1定义的波数矢量，q为步骤6.1定义的当前迭代次数；
步骤10.1.8、采用公式
Figure FDA0003534881150000091
计算
Figure FDA0003534881150000092
关于
Figure FDA0003534881150000093
的偏导数，记为
Figure FDA0003534881150000094
Figure FDA0003534881150000095
为步骤10.1.5计算出的
Figure FDA0003534881150000096
对应的回波信号值，
Figure FDA0003534881150000097
为步骤10.1.4计算出的第k个APC对应的在像素点网格ΩM×N中第m行n列像素点回波的索引值，
Figure FDA0003534881150000098
为步骤10.1.3计算的APC误差向量，
Figure FDA0003534881150000099
为步骤10.1.1计算的匀速APC向量，G为步骤1定义的波数矢量，k＝1，2，…，K，K为步骤1定义的慢时刻数，Ro为步骤1定义的参考距离，gx为距离向波数，gy为方位向波数，m＝1，2，…，M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1，2，…，N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数，q为步骤6.1定义的当前迭代次数；
步骤10.1.9、
Figure FDA00035348811500000910
计算第k个慢时刻对应的中间向量，记为
Figure FDA00035348811500000911
Figure FDA00035348811500000912
为步骤10.1.8计算的
Figure FDA00035348811500000913
关于
Figure FDA00035348811500000914
的偏导数，
Figure FDA00035348811500000915
为步骤10.1.7计算的网格ΩM×N中第m行n列像素点像素点在时域的后向投影值，k＝1，2，…，K，K为步骤1定义的慢时刻数，m＝1，2，…，M，M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数，n＝1，2，…，N，N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数，(·)*表示向量共轭运算，Re[·]表示取实数部分，q为步骤6.1定义的当前迭代次数；
步骤10.1.10、重复步骤10.1.1到步骤10.1.9，对所有的k＝1，2，…，K，K为步骤1定义的慢时刻数，计算出
Figure FDA00035348811500000916
Figure FDA00035348811500000917
为第1个慢时刻对应的第q+1次迭代中间向量，
Figure FDA0003534881150000101
为第2个慢时刻对应的第q+1次迭代中间向量，
Figure FDA0003534881150000102
为第K个慢时刻对应的第q+1次迭代中间向量，q为步骤6.1定义的当前迭代次数；
步骤10.1.11、采用公式
Figure FDA0003534881150000103
计算第q+1次迭代梯度向量，记为
Figure FDA0003534881150000104
Figure FDA0003534881150000105
为步骤10.1.10计算出的第1个慢时刻对应的第q+1次迭代中间向量，
Figure FDA0003534881150000106
为步骤10.1.10计算出的第2个慢时刻对应的第q+1次迭代中间向量，
Figure FDA0003534881150000107
为步骤10.1.10计算出的第K个慢时刻对应的第q+1次迭代中间向量，(·)T表示向量转置运算，q为步骤6.1定义的当前迭代次数；
步骤10.2、采用公式
Figure FDA0003534881150000108
计算中间参数，
Figure FDA0003534881150000109
为步骤10.1.11计算的第q+1次迭代梯度向量，
Figure FDA00035348811500001010
为步骤6.4定义的第q次迭代梯度向量，q为步骤6.1定义的当前迭代次数，||·||表示求向量二范数运算；
步骤10.3、采用公式
Figure FDA00035348811500001011
计算第q+1次迭代搜索方向，记为dq+1，
Figure FDA00035348811500001012
为步骤10.1.11计算出的第q+1次迭代梯度向量，βq为步骤10.2计算出的中间参数，dq为步骤6.3定义的第q次迭代搜索方向，q为步骤6.1定义的当前迭代次数；
步骤11、更新迭代次数，进入下一次迭代
采用公式q←q+1更新下一次迭代次数，运算符←表示赋值操作，重复步骤7～步骤11，直到迭代结束。
CN201910334110.4A 2019-04-24 2019-04-24 一种合成孔径雷达频域bp算法的运动误差补偿方法 Active CN110109107B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910334110.4A CN110109107B (zh) 2019-04-24 2019-04-24 一种合成孔径雷达频域bp算法的运动误差补偿方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910334110.4A CN110109107B (zh) 2019-04-24 2019-04-24 一种合成孔径雷达频域bp算法的运动误差补偿方法

Publications (2)

Publication Number Publication Date
CN110109107A CN110109107A (zh) 2019-08-09
CN110109107B true CN110109107B (zh) 2022-05-31

Family

ID=67486490

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910334110.4A Active CN110109107B (zh) 2019-04-24 2019-04-24 一种合成孔径雷达频域bp算法的运动误差补偿方法

Country Status (1)

Country Link
CN (1) CN110109107B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110568410B (zh) * 2019-10-09 2021-08-31 上海无线电设备研究所 一种空间频率色散的微波雷达超分辨方法
CN113126057B (zh) * 2021-04-20 2022-09-16 哈尔滨工业大学 一种基于调频率估计的sar运动补偿方法
CN113484862B (zh) * 2021-08-04 2023-10-17 电子科技大学 一种自适应的高分宽幅sar清晰重构成像方法
CN113805176B (zh) * 2021-09-18 2022-04-26 哈尔滨工业大学 基于锐度分析和成像投影平面选取的最优成像时间段选取方法

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR930010562A (ko) * 1991-11-27 1993-06-22 완다 케이. 덴슨-로우 합성 어레이 레이다에서의 포커스 에러의 자동 보정 방법
US7456780B1 (en) * 2006-07-26 2008-11-25 Science Applications International Corporation Method and system for developing and using an image reconstruction algorithm for detecting and imaging moving targets
CN103913741A (zh) * 2014-03-18 2014-07-09 电子科技大学 一种合成孔径雷达高效自聚焦后向投影bp方法
CN104181514A (zh) * 2014-08-18 2014-12-03 电子科技大学 一种合成孔径雷达高精度运动补偿方法
CN104316923A (zh) * 2014-10-14 2015-01-28 南京航空航天大学 针对合成孔径雷达bp成像的自聚焦方法
CN104730520A (zh) * 2015-03-27 2015-06-24 电子科技大学 基于子孔径合成的圆周sar后向投影自聚焦方法
CN107015225A (zh) * 2017-03-22 2017-08-04 电子科技大学 一种基于自聚焦的sar平台初始高度误差估计方法
CN107748362A (zh) * 2017-10-10 2018-03-02 电子科技大学 一种基于最大锐度的线阵sar快速自聚焦成像方法
CN109031222A (zh) * 2018-07-09 2018-12-18 中国科学院电子学研究所 重航过阵列合成孔径雷达三维成像运动误差补偿方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8816896B2 (en) * 2012-05-11 2014-08-26 Raytheon Company On-board INS quadratic correction method using maximum likelihood motion estimation of ground scatterers from radar data

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR930010562A (ko) * 1991-11-27 1993-06-22 완다 케이. 덴슨-로우 합성 어레이 레이다에서의 포커스 에러의 자동 보정 방법
US7456780B1 (en) * 2006-07-26 2008-11-25 Science Applications International Corporation Method and system for developing and using an image reconstruction algorithm for detecting and imaging moving targets
CN103913741A (zh) * 2014-03-18 2014-07-09 电子科技大学 一种合成孔径雷达高效自聚焦后向投影bp方法
CN104181514A (zh) * 2014-08-18 2014-12-03 电子科技大学 一种合成孔径雷达高精度运动补偿方法
CN104316923A (zh) * 2014-10-14 2015-01-28 南京航空航天大学 针对合成孔径雷达bp成像的自聚焦方法
CN104730520A (zh) * 2015-03-27 2015-06-24 电子科技大学 基于子孔径合成的圆周sar后向投影自聚焦方法
CN107015225A (zh) * 2017-03-22 2017-08-04 电子科技大学 一种基于自聚焦的sar平台初始高度误差估计方法
CN107748362A (zh) * 2017-10-10 2018-03-02 电子科技大学 一种基于最大锐度的线阵sar快速自聚焦成像方法
CN109031222A (zh) * 2018-07-09 2018-12-18 中国科学院电子学研究所 重航过阵列合成孔径雷达三维成像运动误差补偿方法

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
A Less-Memory and High-Efficiency Autofocus Back Projection Algorithm for SAR Imaging;Kebin Hu等;《 IEEE Geoscience and Remote Sensing Letters》;20141120;第12卷(第4期);全文 *
Autofocus for Correcting Three Dimensional Trajectory Deviations in Synthetic Aperture Radar;Ran, L等;《2016 CIE INTERNATIONAL CONFERENCE ON RADAR (RADAR)》;20161013;全文 *
Fast Back-Projection Autofocus for Linear Array SAR 3-D imaging via Maximum Sharpness;Wei, SJ等;《2018 IEEE RADAR CONFERENCE (RADARCONF18)》;20181231;全文 *
基于图像强度最优的SAR高精度运动补偿方法;胡克彬等;《雷达学报》;20150228;第4卷(第1期);全文 *
快速后向投影合成孔径雷达成像的自聚焦方法;张磊等;《西安电子科技大学学报（自然科学版）》;20140314;第41卷(第1期);全文 *
机载SAR BP算法成像的运动补偿及GPU并行化实现研究;刘斌;《中国优秀博硕士学位论文全文数据库(硕士)信息科技辑》;20140115;全文 *

Also Published As

Publication number Publication date
CN110109107A (zh) 2019-08-09

Similar Documents

Publication Publication Date Title
CN110109107B (zh) 一种合成孔径雷达频域bp算法的运动误差补偿方法
Zhang et al. A robust motion compensation approach for UAV SAR imagery
Ash An autofocus method for backprojection imagery in synthetic aperture radar
Liu et al. Superresolution ISAR imaging based on sparse Bayesian learning
Rao et al. Adaptive sparse recovery by parametric weighted L $ _ {1} $ minimization for ISAR imaging of uniformly rotating targets
Xu et al. High-resolution inverse synthetic aperture radar imaging and scaling with sparse aperture
CN105842694B (zh) 一种基于ffbp sar成像的自聚焦方法
Rao et al. Parametric sparse representation method for ISAR imaging of rotating targets
Li et al. An ISAR imaging algorithm for maneuvering targets with low SNR based on parameter estimation of multicomponent quadratic FM signals and nonuniform FFT
CN109597075B (zh) 一种基于稀疏阵列的成像方法及成像装置
Andersson et al. Fast Fourier methods for synthetic aperture radar imaging
CN105699969A (zh) 基于广义高斯约束的最大后验估计角超分辨成像方法
CN110095787B (zh) 基于MEA和deramp的SAL全孔径成像方法
CN104251990A (zh) 合成孔径雷达自聚焦方法
CN113671492B (zh) 一种面向机动平台前视成像的samp重构方法
Güngör et al. Autofocused compressive SAR imaging based on the alternating direction method of multipliers
CN103809180B (zh) 用于InSAR地形测量的方位向预滤波处理方法
CN110879391B (zh) 基于电磁仿真和弹载回波仿真的雷达图像数据集制作方法
Zeng et al. Two‐dimensional autofocus technique for high‐resolution spotlight synthetic aperture radar
CN103076608B (zh) 轮廓增强的聚束式合成孔径雷达成像方法
CN113608218B (zh) 一种基于后向投影原理的频域干涉相位稀疏重构方法
Chen et al. Iterative minimum entropy algorithm for refocusing of moving targets in SAR images
CN111722227A (zh) 基于近似观测矩阵的聚束sar压缩感知成像方法
CN103792534B (zh) 一种基于先验相位结构知识的sar两维自聚焦方法
CN107229050B (zh) 一种基于极坐标格式的雷达成像优化方法

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