CN110109107B - 一种合成孔径雷达频域bp算法的运动误差补偿方法 - Google Patents
一种合成孔径雷达频域bp算法的运动误差补偿方法 Download PDFInfo
- 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
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
技术领域
本发明属于合成孔径雷达(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误差向量为 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、采用公式计算第k个匀速直线APC,记为Pt0为步骤1定义的雷达平台在零时刻的位置向量,V为步骤1定义的雷达平台速度向量,ts(k)为向量ts的第k个元素,k=1,2,…,K,K为步骤1定义的慢时刻数,ts为步骤1定义的慢时刻向量。
步骤5.2、采用公式计算像素点网格ΩM×N中第m行n列像素点的位置向量,记为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、采用公式计算第k个APC误差向量,记为 为步骤1定义的初始APC误差向量,为中第2(k-1)+1个元素,为中第2(k-1)+2个元素,k=1,2,…,K,K为步骤1定义的慢时刻数;为步骤5.1计算出的第k个匀速直线APC,为步骤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、采用公式计算第k个APC对应的在像素点网格ΩM×N中第m行n列像素点回波的索引值,记为gx为距离向波数,gy为方位向波数,Δg为波数间隔,m=1,2,…,M,M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数,n=1,2,…,N,N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数。
步骤5.5、利用在SSK×P第k行数据里找到对应的回波信号值,记为 为步骤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、采用公式计算像素点网格ΩM×N中第m行n列像素点在频域的后向投影值,为步骤5.5计算出的对应的回波信号值,为步骤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定义的慢时刻数,为步骤5.3计算的APC误差向量,为步骤5.1计算的匀速APC向量。
步骤5.7、采用公式计算像素点网格ΩM×N中第m行n列像素点在时域的后向投影值,为步骤5.6计算出的像素点网格ΩM×N中第m行n列像素点在频域的后向投影值,G为步骤1定义的波数矢量,为步骤5.3计算的APC误差向量,为步骤5.2中计算出的第m行n列像素点的位置向量,m=1,2,…,M,M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数,n=1,2,…,N,N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数。
步骤5.8、采用公式计算关于的偏导数,记为 为步骤5.5计算出的对应的回波信号值,为步骤5.4计算出的第k个APC对应的在像素点网格ΩM×N中第m行n列像素点回波的索引值,为步骤5.3计算的APC误差向量,为步骤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、采用公式计算第k个慢时刻对应的中间向量,记为 为步骤5.8计算的关于的偏导数,为步骤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.11、采用公式计算初始梯度向量,记为 为步骤5.10计算出的第1个慢时刻对应的中间向量,为步骤5.10计算出的第2个慢时刻对应的中间向量,为步骤5.10计算出的第K个慢时刻对应的中间向量,(·)T表示向量转置运算。
步骤6、初始化迭代变量
步骤6.1、定义当前迭代次数为q,q=0,1,2,…,Q,Q为步骤1定义的共轭梯度算法最大迭代次数,初始化q=0。
步骤6.3、定义第q次迭代搜索方向为dq,并初始化为dq=d0,d0为步骤5.12计算出的初始搜索方向,q为步骤6.1定义的当前迭代次数。
步骤7、判断迭代是否结束
步骤8、计算最佳搜索步长
利用Armijo算法计算第q次迭代最佳搜索步长,记为λq,λq的大小为1行、2×K列。
步骤9、计算第q+1次迭代APC误差向量
采用公式计算第q+1次迭代APC误差向量,记为λq为步骤8计算出的第q次迭代最佳搜索步长,dq为步骤6.3计算出的第q次迭代搜索方向,为步骤6.2定义的第q次迭代APC误差向量,q为步骤6.1定义的当前迭代次数。
步骤10、计算第q+1次迭代搜索方向
步骤10.1、计算第q+1次迭代梯度向量,具体方法是:
步骤10.1.1、采用公式计算第q+1次迭代、第k个匀速直线APC,记为Pt0为步骤1定义的雷达平台在零时刻的位置向量,V为步骤1定义的雷达平台速度向量,ts(k)为向量ts的第k个元素,k=1,2,…,K,K为步骤1定义的慢时刻数,ts为步骤1定义的慢时刻向量,q为步骤6.1定义的当前迭代次数。
步骤10.1.2、采用公式计算第q+1次迭代、像素点网格ΩM×N中第m行n列像素点的位置向量,记为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、采用公式 计算第q+1次迭代、第k个APC误差向量,记为 为步骤9计算出的第q+1次迭代APC误差向量,为中第2(k-1)+1个元素,为中第2(k-1)+2个元素,k=1,2,…,K,K为步骤1定义的慢时刻数;为步骤10.1.1计算出的第q+1次迭代、第k个匀速直线APC,为步骤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、采用公式计算第q+1次迭代、计算第k个APC对应的在像素点网格ΩM×N中第m行n列像素点回波的索引值,记为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列像素点的回波的索引值在SSK×P第k行数据里找到对应的回波信号值,记为 为步骤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、采用公式计算像素点网格ΩM×N中第m行n列像素点在频域的后向投影值,为步骤10.1.5计算出的对应的回波信号值,为步骤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定义的慢时刻数,为步骤10.1.3计算的APC误差向量,为步骤10.1.1计算的匀速APC向量,q为步骤6.1定义的当前迭代次数。
步骤10.1.7、采用公式计算像素点网格ΩM×N中第m行n列像素点在时域的后向投影值,为步骤10.1.6计算出的像素点网格ΩM×N中第m行n列像素点在频域的后向投影值,为步骤10.1.3计算的APC误差向量,为步骤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、采用公式计算关于的偏导数,记为 为步骤10.1.5计算出的对应的回波信号值,为步骤10.1.4计算出的第k个APC对应的在像素点网格ΩM×N中第m行n列像素点回波的索引值,为步骤10.1.3计算的APC误差向量,为步骤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、计算第k个慢时刻对应的中间向量,记为 为步骤10.1.8计算的关于的偏导数,为步骤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定义的慢时刻数,计算出 为第1个慢时刻对应的第q+1次迭代中间向量,为第2个慢时刻对应的第q+1次迭代中间向量,为第K个慢时刻对应的第q+1次迭代中间向量,q为步骤6.1定义的当前迭代次数。
步骤10.1.11、采用公式计算第q+1次迭代梯度向量,记为 为步骤10.1.10计算出的第1个慢时刻对应的第q+1次迭代中间向量,为步骤10.1.10计算出的第2个慢时刻对应的第q+1次迭代中间向量,为步骤10.1.10计算出的第K个慢时刻对应的第q+1次迭代中间向量,(·)T表示向量转置运算,q为步骤6.1定义的当前迭代次数。
步骤10.3、采用公式计算第q+1次迭代搜索方向,记为dq+1,为步骤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误差向量为 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、采用公式计算第k个匀速直线APC,记为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、采用公式计算像素点网格Ω25×1000中第m行n列像素点的位置向量,记为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、采用公式计算第k个APC误差向量,记为 为步骤1定义的初始APC误差向量,为中第2(k-1)+1个元素,为中第2(k-1)+2个元素,k=1,2,…,500为步骤1定义的慢时刻数;为步骤5.1计算出的第k个匀速直线APC,为步骤5.2计算出的像素点网格Ω25×1000中第m行n列像素点的位置向量,m=1,2,…,25为步骤1定义的像素点网格Ω25×1000中X方向的像素点点数,n=1,2,…,1000为步骤1定义的像素点网格Ω25×1000中Y方向的像素点点数。
步骤5.4、采用公式计算第k个APC对应的在像素点网格Ω25×1000中第m行n列像素点回波的索引值,记为gx为距离向波数,gy为方位向波数,Δg为波数间隔,m=1,2,…,25为步骤1定义的像素点网格Ω25×1000中X方向的像素点点数,n=1,2,…,1000为步骤1定义的像素点网格Ω25×1000中Y方向的像素点点数。
步骤5.5、利用在SS500×16384第k行数据里找到对应的回波信号值,记为 为步骤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、采用公式计算像素点网格Ω25×1000中第m行n列像素点在频域的后向投影值,为步骤5.5计算出的对应的回波信号值,为步骤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定义的慢时刻数,为步骤5.3计算的APC误差向量,为步骤5.1计算的匀速APC向量。
步骤5.7、采用公式计算像素点网格Ω25×1000中第m行n列像素点在时域的后向投影值,为步骤5.6计算出的像素点网格Ω25×1000中第m行n列像素点在频域的后向投影值,G为步骤1定义的波数矢量,为步骤5.3计算的APC误差向量,为步骤5.2中计算出的第m行n列像素点的位置向量,m=1,2,…,25为步骤1定义的像素点网格Ω25×1000中X方向的像素点点数,n=1,2,…,1000为步骤1定义的像素点网格Ω25×1000中Y方向的像素点点数。
步骤5.8、采用公式计算关于的偏导数,记为 为步骤5.5计算出的对应的回波信号值,为步骤5.4计算出的第k个APC对应的在像素点网格Ω25×1000中第m行n列像素点回波的索引值,为步骤5.3计算的APC误差向量,为步骤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、采用公式计算第k个慢时刻对应的中间向量,记为 为步骤5.8计算的关于的偏导数,为步骤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定义的慢时刻数,计算出 为第1个慢时刻对应的中间向量,为第2个慢时刻对应的中间向量,为第500个慢时刻对应的中间向量。
步骤5.11、采用公式计算初始梯度向量,记为 为步骤5.10计算出的第1个慢时刻对应的中间向量,为步骤5.10计算出的第2个慢时刻对应的中间向量,为步骤5.10计算出的第500个慢时刻对应的中间向量,(·)T表示向量转置运算。
步骤6、初始化迭代变量
步骤6.1、定义当前迭代次数为q,q=0,1,2,…,100,100为步骤1定义的共轭梯度算法最大迭代次数,初始化q=0。
步骤6.3、定义第q次迭代搜索方向为dq,并初始化为dq=d0,d0为步骤5.12计算出的初始搜索方向,q为步骤6.1定义的当前迭代次数。
步骤7、判断迭代是否结束
步骤8、计算最佳搜索步长
利用Armijo算法计算第q次迭代最佳搜索步长,记为λq,λq的大小为1行、2×500列。
步骤9、计算第q+1次迭代APC误差向量
采用公式计算第q+1次迭代APC误差向量,记为λq为步骤8计算出的第q次迭代最佳搜索步长,dq为步骤6.3计算出的第q次迭代搜索方向,为步骤6.2定义的第q次迭代APC误差向量,q为步骤6.1定义的当前迭代次数。
步骤10、计算第q+1次迭代搜索方向
步骤10.1、计算第q+1次迭代梯度向量,具体方法是:
步骤10.1.1、采用公式计算第q+1次迭代、第k个匀速直线APC,记为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、采用公式 计算第q+1次迭代、像素点网格Ω25×1000中第m行n列像素点的位置向量,记为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、采用公式 计算第q+1次迭代、第k个APC误差向量,记为 为步骤9计算出的第q+1次迭代APC误差向量,为中第2(k-1)+1个元素,为中第2(k-1)+2个元素,k=1,2,…,500,其为步骤1定义的慢时刻数;为步骤10.1.1计算出的第q+1次迭代、第k个匀速直线APC,为步骤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、采用公式计算第q+1次迭代、计算第k个APC对应的在像素点网格Ω25×1000中第m行n列像素点回波的索引值,记为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列像素点的回波的索引值在SS500×16384第k行数据里找到对应的回波信号值,记为 为步骤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、采用公式计算像素点网格Ω25×1000中第m行n列像素点在频域的后向投影值,为步骤10.1.5计算出的对应的回波信号值,为步骤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定义的慢时刻数,为步骤10.1.3计算的APC误差向量,为步骤10.1.1计算的匀速APC向量,q为步骤6.1定义的当前迭代次数。
步骤10.1.7、采用公式计算像素点网格Ω25×1000中第m行n列像素点在时域的后向投影值,为步骤10.1.6计算出的像素点网格Ω25×1000中第m行n列像素点在频域的后向投影值,为步骤10.1.3计算的APC误差向量,为步骤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、采用公式计算关于的偏导数,记为 为步骤10.1.5计算出的对应的回波信号值,为步骤10.1.4计算出的第k个APC对应的在像素点网格Ω25×1000中第m行n列像素点回波的索引值,为步骤10.1.3计算的APC误差向量,为步骤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、计算第k个慢时刻对应的中间向量,记为 为步骤10.1.8计算的关于的偏导数,为步骤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定义的慢时刻数,计算出 为第1个慢时刻对应的第q+1次迭代中间向量,为第2个慢时刻对应的第q+1次迭代中间向量,为第500个慢时刻对应的第q+1次迭代中间向量,q为步骤6.1定义的当前迭代次数。
步骤10.1.11、采用公式计算第q+1次迭代梯度向量,记为 为步骤10.1.10计算出的第1个慢时刻对应的第q+1次迭代中间向量,为步骤10.1.10计算出的第2个慢时刻对应的第q+1次迭代中间向量,为步骤10.1.10计算出的第500个慢时刻对应的第q+1次迭代中间向量,(·)T表示向量转置运算,q为步骤6.1定义的当前迭代次数。
步骤10.3、采用公式计算第q+1次迭代搜索方向,记为dq+1,为步骤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误差向量为 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、采用公式计算第k个匀速直线APC,记为Pt0为步骤1定义的雷达平台在零时刻的位置向量,V为步骤1定义的雷达平台速度向量,ts(k)为向量ts的第k个元素,k=1,2,…,K,K为步骤1定义的慢时刻数,ts为步骤1定义的慢时刻向量;
步骤5.2、采用公式计算像素点网格ΩM×N中第m行n列像素点的位置向量,记为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、采用公式计算第k个APC误差向量,记为 为步骤1定义的初始APC误差向量,为中第2(k-1)+1个元素,为中第2(k-1)+2个元素,k=1,2,...,K,K为步骤1定义的慢时刻数;为步骤5.1计算出的第k个匀速直线APC,为步骤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、采用公式计算第k个APC对应的在像素点网格ΩM×N中第m行n列像素点回波的索引值,记为gx为距离向波数,gy为方位向波数,Δg为波数间隔,m=1,2,…,M,M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数,n=1,2,…,N,N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数;
步骤5.5、利用在SSK×P第k行数据里找到对应的回波信号值,记为 为步骤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、采用公式计算像素点网格ΩM×N中第m行n列像素点在频域的后向投影值,为步骤5.5计算出的对应的回波信号值,为步骤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定义的慢时刻数,为步骤5.3计算的APC误差向量,为步骤5.1计算的匀速APC向量;
步骤5.7、采用公式计算像素点网格ΩM×N中第m行n列像素点在时域的后向投影值,为步骤5.6计算出的像素点网格ΩM×N中第m行n列像素点在频域的后向投影值,G为步骤1定义的波数矢量,为步骤5.3计算的APC误差向量,为步骤5.2中计算出的第m行n列像素点的位置向量,m=1,2,…,M,M为步骤1定义的像素点网格ΩM×N中X方向的像素点点数,n=1,2,…,N,N为步骤1定义的像素点网格ΩM×N中Y方向的像素点点数;
步骤5.8、采用公式计算关于的偏导数,记为 为步骤5.5计算出的对应的回波信号值,为步骤5.4计算出的第k个APC对应的在像素点网格ΩM×N中第m行n列像素点回波的索引值,为步骤5.3计算的APC误差向量,为步骤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、采用公式计算第k个慢时刻对应的中间向量,记为 为步骤5.8计算的关于的偏导数,为步骤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.11、采用公式计算初始梯度向量,记为 为步骤5.10计算出的第1个慢时刻对应的中间向量,为步骤5.10计算出的第2个慢时刻对应的中间向量,为步骤5.10计算出的第K个慢时刻对应的中间向量,(·)T表示向量转置运算;
步骤6、初始化迭代变量
步骤6.1、定义当前迭代次数为q,q=0,1,2,…,Q,Q为步骤1定义的共轭梯度算法最大迭代次数,初始化q=0;
步骤6.3、定义第q次迭代搜索方向为dq,并初始化为dq=d0,d0为步骤5.12计算出的初始搜索方向,q为步骤6.1定义的当前迭代次数;
步骤7、判断迭代是否结束
步骤8、计算最佳搜索步长
利用Armijo算法计算第q次迭代最佳搜索步长,记为λq,λq的大小为1行、2×K列;
步骤9、计算第q+1次迭代APC误差向量
采用公式计算第q+1次迭代APC误差向量,记为λq为步骤8计算出的第q次迭代最佳搜索步长,dq为步骤6.3计算出的第q次迭代搜索方向,为步骤6.2定义的第q次迭代APC误差向量,q为步骤6.1定义的当前迭代次数;
步骤10、计算第q+1次迭代搜索方向
步骤10.1、计算第q+1次迭代梯度向量,具体方法是:
步骤10.1.1、采用公式计算第q+1次迭代、第k个匀速直线APC,记为Pt0为步骤1定义的雷达平台在零时刻的位置向量,V为步骤1定义的雷达平台速度向量,ts(k)为向量ts的第k个元素,k=1,2,…,K,K为步骤1定义的慢时刻数,ts为步骤1定义的慢时刻向量,q为步骤6.1定义的当前迭代次数;
步骤10.1.2、采用公式计算第q+1次迭代、像素点网格ΩM×N中第m行n列像素点的位置向量,记为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、采用公式计算第q+1次迭代、第k个APC误差向量,记为 为步骤9计算出的第q+1次迭代APC误差向量,为中第2(k-1)+1个元素,为中第2(k-1)+2个元素,k=1,2,...,K,K为步骤1定义的慢时刻数;为步骤10.1.1计算出的第q+1次迭代、第k个匀速直线APC,为步骤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、采用公式计算第q+1次迭代、计算第k个APC对应的在像素点网格ΩM×N中第m行n列像素点回波的索引值,记为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列像素点的回波的索引值在SSK×P第k行数据里找到对应的回波信号值,记为 为步骤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、采用公式计算像素点网格ΩM×N中第m行n列像素点在频域的后向投影值,为步骤10.1.5计算出的对应的回波信号值,为步骤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定义的慢时刻数,为步骤10.1.3计算的APC误差向量,为步骤10.1.1计算的匀速APC向量,q为步骤6.1定义的当前迭代次数;
步骤10.1.7、采用公式计算像素点网格ΩM×N中第m行n列像素点在时域的后向投影值,为步骤10.1.6计算出的像素点网格ΩM×N中第m行n列像素点在频域的后向投影值,为步骤10.1.3计算的APC误差向量,为步骤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、采用公式计算关于的偏导数,记为 为步骤10.1.5计算出的对应的回波信号值,为步骤10.1.4计算出的第k个APC对应的在像素点网格ΩM×N中第m行n列像素点回波的索引值,为步骤10.1.3计算的APC误差向量,为步骤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、计算第k个慢时刻对应的中间向量,记为 为步骤10.1.8计算的关于的偏导数,为步骤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定义的慢时刻数,计算出 为第1个慢时刻对应的第q+1次迭代中间向量,为第2个慢时刻对应的第q+1次迭代中间向量,为第K个慢时刻对应的第q+1次迭代中间向量,q为步骤6.1定义的当前迭代次数;
步骤10.1.11、采用公式计算第q+1次迭代梯度向量,记为 为步骤10.1.10计算出的第1个慢时刻对应的第q+1次迭代中间向量,为步骤10.1.10计算出的第2个慢时刻对应的第q+1次迭代中间向量,为步骤10.1.10计算出的第K个慢时刻对应的第q+1次迭代中间向量,(·)T表示向量转置运算,q为步骤6.1定义的当前迭代次数;
步骤10.3、采用公式计算第q+1次迭代搜索方向,记为dq+1,为步骤10.1.11计算出的第q+1次迭代梯度向量,βq为步骤10.2计算出的中间参数,dq为步骤6.3定义的第q次迭代搜索方向,q为步骤6.1定义的当前迭代次数;
步骤11、更新迭代次数,进入下一次迭代
采用公式q←q+1更新下一次迭代次数,运算符←表示赋值操作,重复步骤7~步骤11,直到迭代结束。
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)
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)
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)
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 |
-
2019
- 2019-04-24 CN CN201910334110.4A patent/CN110109107B/zh active Active
Patent Citations (9)
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)
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算法的运动误差补偿方法 | |
CN105842694B (zh) | 一种基于ffbp sar成像的自聚焦方法 | |
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 | |
CN103901429B (zh) | 基于稀疏孔径的机动目标逆合成孔径雷达成像方法 | |
Rao et al. | Adaptive sparse recovery by parametric weighted L $ _ {1} $ minimization for ISAR imaging of uniformly rotating targets | |
Rao et al. | Parametric sparse representation method for ISAR imaging of rotating targets | |
CN103760558B (zh) | 一种太赫兹雷达isar成像方法 | |
CN105699969B (zh) | 基于广义高斯约束的最大后验估计角超分辨成像方法 | |
CN103454638B (zh) | 一种圆迹合成孔径雷达三维层析成像方法 | |
CN109597075B (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 | |
Jiu et al. | Joint ISAR imaging and cross-range scaling method based on compressive sensing with adaptive dictionary | |
CN110879391B (zh) | 基于电磁仿真和弹载回波仿真的雷达图像数据集制作方法 | |
Zeng et al. | Two‐dimensional autofocus technique for high‐resolution spotlight synthetic aperture radar | |
CN112147608A (zh) | 一种快速高斯网格化非均匀fft穿墙成像雷达bp方法 | |
CN113608218B (zh) | 一种基于后向投影原理的频域干涉相位稀疏重构方法 | |
CN103076608B (zh) | 轮廓增强的聚束式合成孔径雷达成像方法 | |
Chen et al. | Iterative minimum entropy algorithm for refocusing of moving targets in SAR images | |
CN107229050B (zh) | 一种基于极坐标格式的雷达成像优化方法 | |
CN104181514B (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 |