CN113176570A - 一种斜视sar时域成像自聚焦方法 - Google Patents
一种斜视sar时域成像自聚焦方法 Download PDFInfo
- Publication number
- CN113176570A CN113176570A CN202110432393.3A CN202110432393A CN113176570A CN 113176570 A CN113176570 A CN 113176570A CN 202110432393 A CN202110432393 A CN 202110432393A CN 113176570 A CN113176570 A CN 113176570A
- Authority
- CN
- China
- Prior art keywords
- matrix
- size
- data
- azinum
- rannum
- 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.)
- Granted
Links
- 238000003384 imaging method Methods 0.000 title claims abstract description 74
- 238000000034 method Methods 0.000 title claims abstract description 33
- 208000004350 Strabismus Diseases 0.000 title claims abstract description 24
- 239000011159 matrix material Substances 0.000 claims abstract description 122
- 238000005070 sampling Methods 0.000 claims abstract description 33
- 238000007906 compression Methods 0.000 claims abstract description 15
- 238000012545 processing Methods 0.000 claims abstract description 11
- 230000006835 compression Effects 0.000 claims abstract description 9
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 3
- 230000009191 jumping Effects 0.000 claims description 3
- 230000009466 transformation Effects 0.000 claims description 3
- 230000006872 improvement Effects 0.000 description 9
- 238000010586 diagram Methods 0.000 description 5
- 230000000694 effects Effects 0.000 description 3
- 230000004048 modification Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 238000012952 Resampling Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000015556 catabolic process Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 238000006731 degradation reaction Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 238000013441 quality evaluation Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
- G01S13/88—Radar or analogous systems specially adapted for specific applications
- G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
- G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
- G01S13/904—SAR modes
- G01S13/9041—Squint mode
-
- 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
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Remote Sensing (AREA)
- Radar, Positioning & Navigation (AREA)
- Data Mining & Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Theoretical Computer Science (AREA)
- General Engineering & Computer Science (AREA)
- Software Systems (AREA)
- Databases & Information Systems (AREA)
- Computer Networks & Wireless Communication (AREA)
- Algebra (AREA)
- Computing Systems (AREA)
- Electromagnetism (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
本发明公开了一种斜视SAR时域成像自聚焦方法,包括:获取成像参数和回波数据S0;对回波数据S0进行距离向脉冲压缩,得到距离向聚焦信号Src;对距离向聚焦信号Src进行距离向插值,得到插值后的信号S3;建立旋转的地面成像网络;输入信号S3和地面成像网络坐标矩阵,逐像素点补偿相位因子,得到Sjj;遍历所有方位向采样点数Na,将各方位采样点处得到的矩阵Sjj相加得到最终的BP聚焦数据SBP;输入BP聚焦数据SBP进行PGA自聚焦处理。本发明利用地面网格与散焦的点目标方位旁瓣方向的正交使PGA估计的二次误差相位更准确,利用此方法估计的误差相位对散焦图像进行补偿,可实现散焦图像的精确聚焦。
Description
技术领域
本发明涉及信号处理技术领域,具体涉及一种斜视SAR时域成像自聚焦方法。
背景技术
在合成孔径雷达(SAR)成像中,回波多普勒信号的相位误差会造成图像的散焦,影响成像质量。随着SAR成像分辨率的提高,回波相位误差对成像质量的影响也愈显严重。虽然依靠惯性导航运动补偿系统能够补偿一定的相位误差,但是残留相位误差仍然有可能造成图像质量的严重下降。相比传统的自聚焦算法,如Map-Drift,Phase Difference等,相位梯度自聚焦算法(PGA)不需要场景中具有强反射点,并且对低阶、高阶以及随机误差都能够较好的进行补偿。因此,PGA算法自提出以来已经在SAR成像领域有了十分广泛的应用。
原始PGA算法主要分为四个部分,分别是圆移、加窗、相位梯度估计和迭代相位校正处理。其中,相位梯度估计是最重要的步骤,PGA算法通过方位向傅里叶逆变换将图像域数据返回数据域,根据数据域强散射目标的相位变化规律,计算得到成像中存在的相位误差,不同距离门的叠加效果将提升误差估计的精度。由此可见算法需要沿二维数据域中的某一个维度进行估计,即需要方位向为水平或垂直方向,从而才能进行估计。但对于斜视情况下的后向投影(BP)算法成像,最终目标的扩展函数方向与斜视角在地面的投影角度有关,使得方位向不再是严格的正交方向,从而导致斜视情况下SAR系统自聚焦后成像质量不能够明显提升。
发明内容
为解决斜视情况下采用BP算法进行目标成像后,由于点目标扩展函数二维不正交导致的传统PGA方法失效的问题,本发明提供一种斜视SAR时域成像自聚焦方法,能够完成斜视情况下时域成像的精确聚焦处理。
本发明公开了一种斜视SAR时域成像自聚焦方法,包括:
步骤1、获取成像参数和回波数据S0;
步骤2、对所述回波数据S0进行距离向脉冲压缩,得到距离向聚焦信号Src;
步骤3、对所述距离向聚焦信号Src进行距离向插值,得到插值后的信号S3;
步骤4、建立旋转的地面成像网络;
步骤5、输入信号S3和地面成像网络坐标矩阵,逐像素点补偿相位因子,得到Sjj;
步骤6、遍历所有方位向采样点数Na,将各方位采样点处得到的矩阵Sjj相加得到最终的BP聚焦数据SBP;
步骤7、输入BP聚焦数据SBP进行PGA自聚焦处理。
作为本发明的进一步改进,在所述步骤1中,
所述成像参数包括距离采样率Fr、调频率Kr、脉冲重复频率PRF、距离分辨率和方位分辨率指标ρr、ρa、雷达俯仰角β、场景中心点的位置坐标xt、yt、zt、回波数据距离向采样点数Nr,回波数据方位向采样点数Na,其中Na/2为整数;雷达三维位置坐标XT、YT、ZT以及距离门开启时刻的距离R0;
回波数据S0是大小为Na×Nr的复数矩阵。
作为本发明的进一步改进,所述步骤2,具体包括:
步骤201、调用MATLAB自带的FFT函数,输入回波信号S0,对每一行数据进行Nr点傅里叶变换,得到距离向频域的回波数据S1;
步骤202、根据公式(1)构造H1,H1是大小为1×Nr的行向量;fr是大小为1×Nr的行向量,第i个点的表达式如下:
fr=i×Fr/Nr
其中,j表示虚数单位,i=1,2,…,Nr表示距离向采样点序列;
步骤203、调用MATLAB中自带的复写函数repmat,输入步骤202中得到的大小为1×Nr的H1,设置复写行数参数Na,输出大小为Na×Nr的矩阵;与步骤201中得到的信号S1点乘;然后调用MATLAB中自带的IFFT函数,输入每一行的点乘结果,利用距离向快速傅里叶逆变换将点乘的结果变换到距离向时域,输出已完成距离向聚焦的信号Src。
作为本发明的进一步改进,所述步骤3,具体包括:
步骤301、调用MATLAB自带的FFT函数,输入距离脉冲压缩后的距离向聚焦信号Src,对输入信号的每一行进行Nr点距离向快速傅里叶变换得到信号SFC,SFC是大小为Na×Nr的矩阵;
步骤302、设定距离插值倍数M1,在步骤301得到的信号SFC两端分别补行数为Na,列数为Nr×(M1-1)/2个零,构成大小为Na×(Nr×M1)的矩阵S2;
步骤303、调用MATLAB自带的IFFT函数,输入补零后的数据S2,逐行进行Nr点快速傅里叶逆变换,完成对距离脉冲压缩后的距离向聚焦信号Src的距离向插值,得到插值后的信号S3,S3是大小为Na×(Nr×M1)的矩阵。
作为本发明的进一步改进,所述步骤4,具体包括:
步骤401、设定地面成像网格尺寸:
根据步骤1中载入的分辨率指标ρr、ρa,设定在X方向地面网格的尺寸为Xp,Y方向地面网格的尺寸为Yp,X方向网格像素间隔为Δx,Y方向网格像素间隔为Δy,地面成像网格点数为AziNum×RanNum,具体由如下公式得到:
Xp=AziNum·Δx
Yp=RanNum·Δy (3)
步骤402、对X方向和Y方向的网格矩阵化:
初始设置X是大小为1×AziNum的矩阵,Y是大小为1×RanNum的矩阵,以场景中心点(xt,yt)作为网格的中心,m=1,2,…,AziNum,n=1,2,…,RanNum,根据公式(4)对网格进行离散化表示为:
X(1,m)=xt-XP+(2·Xp)×(m-1)/(RanNum-1)
Y(1,n)=yt-YP+(2·Yp)×(n-1)/(AziNum-1) (4)
步骤403、调用MATLAB中自带的repmat函数,输入步骤402生成的X网格,首先将X网格转置成为AziNum×1的矩阵,然后设置复写行数参数RanNum,输出大小为AziNum×RanNum的矩阵X(m,n);同理,调用repmat函数,输入步骤402生成的Y网格,设置复写行数AziNum,输出大小为AziNum×RanNum的矩阵Y(m,n);
步骤404、计算旋转角度θs,构造旋转矩阵H:
输入步骤1中加载的俯仰角参数β,其中β是大小为1×Na的向量、雷达三维位置坐标XT、YT、ZT,其中XT、YT、ZT均为大小为1×Na的向量以及公式(2)计算得到的斜距R1,利用俯仰角向量β的第Na/2点取值根据公式(5)计算得到中心时刻斜距在地面上的投影长度RL:
输入雷达平台在Y方向坐标YT的第Na/2点的取值YT(Na/2)以及中公式(5)得到的投影长度RL,根据正弦定理,计算旋转矩阵的旋转角度θs为:
根据输入的旋转角度θs计算旋转矩阵H为:
步骤405、输入步骤403中构造的成像网格X(m,n)、Y(m,n),以及旋转矩阵H,将X(m,n)按列拆分后纵向排列成一个列向量X1,此时X1的大小为(AziNum×RanNum)×1,同样将Y(m,n)按列拆分后纵向排列成一个列向量Y1,此时Y1的大小为(AziNum×RanNum)×1;将X1和Y1合并成一个新的矩阵并转置,得到转置后的矩阵C1,C1是大小为2×(AziNum×RanNum)的矩阵:
步骤406、首先将C1的每一列与(xt,yt)的转置相减得到矩阵C2,然后将旋转矩阵H与C2相乘,最后对矩阵相乘结果的每一列与(xt,yt)的转置相加,输出为旋转后的成像网格坐标C3,此时矩阵C3的大小为2×(AziNum×RanNum);
步骤407、调用MATLAB中reshape函数,输入参数分别为C3矩阵的第一行以及重排后的矩阵大小AziNum×RanNum,输出为重排后的网格矩阵坐标X2;同理调用MATLAB中reshape函数,输入参数分别为C3矩阵的第二行以及重排后的矩阵大小AziNum×RanNum,输出为重排后的网格矩阵坐标Y2:
作为本发明的进一步改进,所述步骤5,具体包括:
步骤501、计算该采样点的雷达位置与成像网格中每一个像素点之间的斜距RT,RT为大小为AziNum×RanNum的矩阵;其中,
输入该采样时刻雷达的三维位置坐标XT、YT、ZT以及由步骤407得到的大小分别为AziNum×RanNum的地面成像网格坐标X2、Y2,根据公式(10)计算雷达平台与每一个网格的斜距,遍历所有网格点,得到每一个方位采样时刻雷达与所有网格点的斜距RT,RT是大小为AziNum×RanNum的矩阵:
步骤502、首先输入步骤501得到的斜距矩阵RT、光速c以及步骤302中的距离插值倍数M1,根据公式(11)、(12)得到tij和dtr;然后调用MATLAB中自带的round函数,根据公式(13)输入tij、dtr、c和R0,输出为Loc,Loc是大小为AziNum×RanNum的矩阵;
tij=2RT/c (11)
dtr=1/M1/Fr (12)
Loc=round((tij-2R0/c)/dtr+1) (13)
步骤503、取步骤303得到的插值后的距离向脉冲压缩信号第jj行数据中第Loc个点处的复像素值σmn与点乘得到Sjj,jj=1,2,…Na;σmn是大小为AziNum×RanNum的矩阵,Sjj是大小为AziNum×RanNum的矩阵。
作为本发明的进一步改进,在所述步骤6中,
作为本发明的进一步改进,所述步骤7,具体包括:
步骤701、输入BP聚焦数据SBP,调用MATLAB中的max求最大值函数,输出SBP中每一个列的最大值以及该最大值在每一列中的位置;
步骤702、输入SBP每一列数据以及该列数据中最大值的位置,调用MATLAB中自带的circshift函数,将每一列数据中的最大值移动到该列的中心位置,遍历所有列,得到移位后的数据S4,S4是大小为AziNum×RanNum的矩阵;
步骤703、输入步骤702中移位后的数据S4,将每一列数据以第AziNum/2个数据为中心,保留长度为W1的数据,其余数据置零,得到加窗处理后的数据S5:
步骤704、调用MATLAB自带的IFFT函数,输入加窗后的数据S5,对S5的每一列进行逆傅里叶变换得到S6;
步骤705、输入S6,遍历S6的AziNum行,对每一行数据取共轭,然后与下一行数据点乘得到S7,S7是大小为(AziNum-1)×RanNum的矩阵,然后将S7每一行中RanNum个数值相加得到S8,S8是大小为(AziNum-1)×1的矩阵;然后调用MATLAB中自带的angle函数,输入S8中每一个值,输出为相邻两行数据之间的相位差是大小为(AziNum-1)×1的矩阵;
步骤706、调用MATLAB中自带的cumsum函数,输入步骤705中得到的输出为各行数据与第一行数据之间的相位差,是大小为(AziNum-1)×1的矩阵;设置是大小为AziNum×1的矩阵,其中第一行的元素为0,其余(AziNum-1)行的元素为中的数值;
步骤707、判断如果中第AziNum个值大于等于π/4,则用S6的每一行数据点乘该行对应的相位差再调用MATLAB中自带的FFT函数对每一行数据进行傅里叶变换得到新的BP聚焦数据,然后重复执行步骤701~707;如果相位差小于π/4,则跳出循环,得到回波数据S0的最终聚焦数据I2。
与现有技术相比,本发明的有益效果为:
本发明的自聚焦方法利用地面网格与散焦的点目标方位旁瓣方向的正交使PGA估计的二次误差相位更准确,利用此方法估计的误差相位对散焦图像进行补偿,可实现散焦图像的精确聚焦,是一种高效、准确的针对斜视情况下目标扩展函数二维不正交的自聚焦方法。
附图说明
图1为本发明一种实施例公开的斜视SAR时域成像自聚焦方法的流程图;
图2为本发明一种实施例公开的机载斜视SAR数据获取的示意图;
图3为本发明一种实施例公开的BP算法成像网格的示意图;
图4为传统自聚焦算法对点目标聚焦成像的结果的示意图;
图5为本发明对仿真点目标聚焦成像的结果的示意图;
图6为本发明中点目标聚焦质量评估结果的示意图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚,下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例是本发明的一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例,都属于本发明保护的范围。
下面结合附图对本发明做进一步的详细描述:
本发明提供一种斜视SAR时域成像自聚焦方法,包括:获取成像参数和回波数据S0;对回波数据S0进行距离向脉冲压缩,得到距离向聚焦信号Src;对距离向聚焦信号Src进行距离向插值,得到插值后的信号S3;建立旋转的地面成像网络;输入信号S3和地面成像网络坐标矩阵,逐像素点补偿相位因子,得到Sjj;遍历所有方位向采样点数Na,将各方位采样点处得到的矩阵Sjj相加得到最终的BP聚焦数据SBP;输入BP聚焦数据SBP进行PGA自聚焦处理。本发明利用地面网格与散焦的点目标方位旁瓣方向的正交使PGA估计的二次误差相位更准确,利用此方法估计的误差相位对散焦图像进行补偿,可实现散焦图像的精确聚焦。
具体的:
如图1所示,本发明提供一种斜视SAR时域成像自聚焦方法,包括:
步骤1、获取成像参数和回波数据S0;
如图2所示,具体包括:
步骤101、读入成像参数,在本实施例中:包括载频f0=5.3GHz/s、脉冲宽度Tr=3μs、距离采样率Fr=108MHz、调频率Kr=30000GHz/s、脉冲重复频率PRF=300Hz、距离分辨率和方位分辨率指标ρr=1.7m、ρa=1.7m、回波数据方位向采样点数Na=512、回波数据距离向采样点数Nr=5000,雷达俯仰角β,β是大小为1×Na的向量,第Na/2个点取值为β=21.5°、雷达三维位置坐标XT、YT、ZT、其中XT、YT、ZT均为1×Na的向量,第Na/2个点取值为XT=-10Km、YT=7.82Km、ZT=5Km、场景中心点的位置坐标xt=0、yt=0、zt=0以及距离门开启时刻的距离R0=7.58Km;成像网格参数Xp=200m、Yp=200m、Δx=1m、Δy=1m、RanNum=200、AziNum=200;自聚焦参数Num=6、W1=50。
步骤102、根据回波数据方位向采样点数Na和回波数据距离向采样点数Nr读入SAR的回波数据S0,S0是大小为Na×Nr的复数矩阵。
步骤2、对回波数据S0进行距离向脉冲压缩,得到距离向聚焦信号Src;其中,
本发明采用传统匹配滤波的脉冲压缩方式对SAR回波信号进行距离向聚焦,具体包括:
步骤201、调用MATLAB自带的FFT函数,输入回波信号S0,对每一行数据进行Nr点傅里叶变换,得到距离向频域的回波数据S1;
步骤202、根据公式(1)构造H1,H1是大小为1×Nr的行向量;fr是大小为1×Nr的行向量,第i个点的表达式如下:
fr=i×Fr/Nr
其中,j表示虚数单位,i=1,2,…,Nr表示距离向采样点序列;
步骤203、调用MATLAB中自带的复写函数repmat,输入步骤202中得到的大小为1×Nr的H1,设置复写行数参数Na,输出大小为Na×Nr的矩阵;与步骤201中得到的信号S1点乘;然后调用MATLAB中自带的IFFT函数,输入每一行的点乘结果,利用距离向快速傅里叶逆变换将点乘的结果变换到距离向时域,输出已完成距离向聚焦的信号Src。
步骤3、对距离向聚焦信号Src进行距离向插值,得到插值后的信号S3;其中,
本发明采用频域插值的方式对距离向聚焦信号Src进行重采样,具体包括:
步骤301、调用MATLAB自带的FFT函数,输入距离脉冲压缩后的距离向聚焦信号Src,对输入信号的每一行进行Nr点距离向快速傅里叶变换得到信号SFC,SFC是大小为Na×Nr的矩阵;
步骤302、设定距离插值倍数M1,在步骤301得到的信号SFC两端分别补行数为Na,列数为Nr×(M1-1)/2个零,构成大小为Na×(Nr×M1)的矩阵S2;
步骤303、调用MATLAB自带的IFFT函数,输入补零后的数据S2,逐行进行Nr点快速傅里叶逆变换,完成对距离脉冲压缩后的距离向聚焦信号Src的距离向插值,得到插值后的信号S3,S3是大小为Na×(Nr×M1)的矩阵。
步骤4、建立旋转的地面成像网络;其中,
首先根据读入成像参数中的雷达轨迹XT、YT、ZT在地面沿X-Y方向建立成像网格,网格像素间隔按照ρr、ρa设定。然后计算雷达轨迹XT、YT、ZT第Na/2点处的取值与场景中心点间的斜距R1;
地面成像网格的建立过程,具体包括:
步骤401、设定地面成像网格尺寸:
根据步骤1中载入的分辨率指标ρr、ρa,设定在X方向地面网格的尺寸为Xp,Y方向地面网格的尺寸为Yp,X方向网格像素间隔为Δx,Y方向网格像素间隔为Δy,地面成像网格点数为AziNum×RanNum,具体由如下公式得到:
Xp=AziNum·Δx
Yp=RanNum·Δy (3)
步骤402、对X方向和Y方向的网格矩阵化:
初始设置X是大小为1×AziNum的矩阵,Y是大小为1×RanNum的矩阵,以场景中心点(xt,yt)作为网格的中心,m=1,2,…,AziNum,n=1,2,…,RanNum,根据公式(4)对网格进行离散化表示为:
X(1,m)=xt-XP+(2·Xp)×(m-1)/(RanNum-1)
Y(1,n)=yt-YP+(2·Yp)×(n-1)/(AziNum-1) (4)
步骤403、调用MATLAB中自带的repmat函数,输入步骤402生成的X网格,首先将X网格转置成为AziNum×1的矩阵,然后设置复写行数参数RanNum,输出大小为AziNum×RanNum的矩阵X(m,n);同理,调用repmat函数,输入步骤402生成的Y网格,设置复写行数AziNum,输出大小为AziNum×RanNum的矩阵Y(m,n);
步骤404、计算旋转角度θs,构造旋转矩阵H:
输入步骤1中加载的俯仰角参数β,其中β是大小为1×Na的向量、雷达三维位置坐标XT、YT、ZT,其中XT、YT、ZT均为大小为1×Na的向量以及公式(2)计算得到的斜距R1,利用俯仰角向量β的第Na/2点取值根据公式(5)计算得到中心时刻斜距在地面上的投影长度RL:
输入雷达平台在Y方向坐标YT的第Na/2点的取值YT(Na/2)以及中公式(5)得到的投影长度RL,根据正弦定理,计算旋转矩阵的旋转角度θs为:
根据输入的旋转角度θs计算旋转矩阵H为:
步骤405、输入步骤403中构造的成像网格X(m,n)、Y(m,n),以及旋转矩阵H,将X(m,n)按列拆分后纵向排列成一个列向量X1,此时X1的大小为(AziNum×RanNum)×1,同样将Y(m,n)按列拆分后纵向排列成一个列向量Y1,此时Y1的大小为(AziNum×RanNum)×1;将X1和Y1合并成一个新的矩阵并转置,得到转置后的矩阵C1,C1是大小为2×(AziNum×RanNum)的矩阵:
步骤406、首先将C1的每一列与(xt,yt)的转置相减得到矩阵C2,然后将旋转矩阵H与C2相乘,最后对矩阵相乘结果的每一列与(xt,yt)的转置相加,输出为旋转后的成像网格坐标C3,此时矩阵C3的大小为2×(AziNum×RanNum),旋转后的成像网格如图3所示。
步骤407、调用MATLAB中reshape函数,输入参数分别为C3矩阵的第一行以及重排后的矩阵大小AziNum×RanNum,输出为重排后的网格矩阵坐标X2;同理调用MATLAB中reshape函数,输入参数分别为C3矩阵的第二行以及重排后的矩阵大小AziNum×RanNum,输出为重排后的网格矩阵坐标Y2:
步骤5、输入信号S3和地面成像网络坐标矩阵,逐像素点补偿相位因子,得到Sjj;
具体包括:
步骤501、计算该采样点的雷达位置与成像网格中每一个像素点之间的斜距RT,RT为大小为AziNum×RanNum的矩阵;其中,
输入该采样时刻雷达的三维位置坐标XT、YT、ZT以及由步骤407得到的大小分别为AziNum×RanNum的地面成像网格坐标X2、Y2,根据公式(10)计算雷达平台与每一个网格的斜距,遍历所有网格点,得到每一个方位采样时刻雷达与所有网格点的斜距RT,RT是大小为AziNum×RanNum的矩阵:
步骤502、首先输入步骤501得到的斜距矩阵RT、光速c以及步骤302中的距离插值倍数M1,根据公式(11)、(12)得到tij和dtr;然后调用MATLAB中自带的round函数,根据公式(13)输入tij、dtr、c和R0,输出为Loc,Loc是大小为AziNum×RanNum的矩阵;
tij=2RT/c (11)
dtr=1/M1/Fr (12)
Loc=round((tij-2R0/c)/dtr+1) (13)
步骤503、取步骤303得到的插值后的距离向脉冲压缩信号第jj行数据中第Loc个点处的复像素值σmn与点乘得到Sjj,σmn是大小为AziNum×RanNum的矩阵,Sjj是大小为AziNum×RanNum的矩阵。
步骤6、遍历所有方位向采样点数Na,将各方位采样点处得到的矩阵Sjj相加得到最终的BP聚焦数据SBP;其中,
步骤7、输入BP聚焦数据SBP进行PGA自聚焦处理;
具体包括:
步骤701、输入BP聚焦数据SBP,调用MATLAB中的max求最大值函数,输出SBP中每一个列的最大值以及该最大值在每一列中的位置;
步骤702、输入SBP每一列数据以及该列数据中最大值的位置,调用MATLAB中自带的circshift函数,将每一列数据中的最大值移动到该列的中心位置,遍历所有列,得到移位后的数据S4,S4是大小为AziNum×RanNum的矩阵;
步骤703、输入步骤702中移位后的数据S4,将每一列数据以第AziNum/2个数据为中心,保留长度为W1的数据,其余数据置零,得到加窗处理后的数据S5:
步骤704、调用MATLAB自带的IFFT函数,输入加窗后的数据S5,对S5的每一列进行逆傅里叶变换得到S6;
步骤705、输入S6,遍历S6的AziNum行,对每一行数据取共轭,然后与下一行数据点乘得到S7,S7是大小为(AziNum-1)×RanNum的矩阵,然后将S7每一行中RanNum个数值相加得到S8,S8是大小为(AziNum-1)×1的矩阵;然后调用MATLAB中自带的angle函数,输入S8中每一个值,输出为相邻两行数据之间的相位差是大小为(AziNum-1)×1的矩阵;
步骤706、调用MATLAB中自带的cumsum函数,输入步骤705中得到的输出为各行数据与第一行数据之间的相位差,是大小为(AziNum-1)×1的矩阵;设置是大小为AziNum×1的矩阵,其中第一行的元素为0,其余(AziNum-1)行的元素为中的数值;
步骤707、判断如果中第AziNum个值大于等于π/4,则用S6的每一行数据点乘该行对应的相位差再调用MATLAB中自带的FFT函数对每一行数据进行傅里叶变换得到新的BP聚焦数据,然后重复执行步骤701~707;如果相位差小于π/4,则跳出循环,得到回波数据S0的最终聚焦数据I2。
其中,本发明的自聚焦结果如图5所示。将其与图4利用传统PGA自聚焦的结果进行对比,可看出聚焦效果有明显改善。图6中点目标的评估结果可看出本发明方法聚焦效果的有效性。
本发明适用于斜视情况下由于BP算法成像导致点目标扩展函数(PSF)二维不正交时PGA算法失效时的自聚焦方法,通过旋转BP成像网格,使点目标方位旁瓣与成像网格中的一维正交,再逐距离门估计二次误差相位进行误差补偿,是一种高效、准确的自聚焦方法。
以上仅为本发明的优选实施例而已,并不用于限制本发明,对于本领域的技术人员来说,本发明可以有各种更改和变化。凡在本发明的精神和原则之内,所作的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。
Claims (8)
1.一种斜视SAR时域成像自聚焦方法,其特征在于,包括:
步骤1、获取成像参数和回波数据S0;
步骤2、对所述回波数据S0进行距离向脉冲压缩,得到距离向聚焦信号Src;
步骤3、对所述距离向聚焦信号Src进行距离向插值,得到插值后的信号S3;
步骤4、建立旋转的地面成像网络;
步骤5、输入信号S3和地面成像网络坐标矩阵,逐像素点补偿相位因子,得到Sjj;
步骤6、遍历所有方位向采样点数Na,将各方位采样点处得到的矩阵Sjj相加得到最终的BP聚焦数据SBP;
步骤7、输入BP聚焦数据SBP进行PGA自聚焦处理。
2.如权利要求1所述的斜视SAR时域成像自聚焦方法,其特征在于,在所述步骤1中,
所述成像参数包括距离采样率Fr、调频率Kr、脉冲重复频率PRF、距离分辨率和方位分辨率指标ρr、ρa、雷达俯仰角β、场景中心点的位置坐标xt、yt、zt、回波数据距离向采样点数Nr,回波数据方位向采样点数Na,其中Na/2为整数;雷达三维位置坐标XT、YT、ZT以及距离门开启时刻的距离R0;
回波数据S0是大小为Na×Nr的复数矩阵。
3.如权利要求2所述的斜视SAR时域成像自聚焦方法,其特征在于,所述步骤2,具体包括:
步骤201、调用MATLAB自带的FFT函数,输入回波信号S0,对每一行数据进行Nr点傅里叶变换,得到距离向频域的回波数据S1;
步骤202、根据公式(1)构造H1,H1是大小为1×Nr的行向量;fr是大小为1×Nr的行向量,第i个点的表达式如下:
fr=i×Fr/Nr
其中,j表示虚数单位,i=1,2,…,Nr表示距离向采样点序列;
步骤203、调用MATLAB中自带的复写函数repmat,输入步骤202中得到的大小为1×Nr的H1,设置复写行数参数Na,输出大小为Na×Nr的矩阵;与步骤201中得到的信号S1点乘;然后调用MATLAB中自带的IFFT函数,输入每一行的点乘结果,利用距离向快速傅里叶逆变换将点乘的结果变换到距离向时域,输出已完成距离向聚焦的信号Src。
4.如权利要求3所述的斜视SAR时域成像自聚焦方法,其特征在于,所述步骤3,具体包括:
步骤301、调用MATLAB自带的FFT函数,输入距离脉冲压缩后的距离向聚焦信号Src,对输入信号的每一行进行Nr点距离向快速傅里叶变换得到信号SFC,SFC是大小为Na×Nr的矩阵;
步骤302、设定距离插值倍数M1,在步骤301得到的信号SFC两端分别补行数为Na,列数为Nr×(M1-1)/2个零,构成大小为Na×(Nr×M1)的矩阵S2;
步骤303、调用MATLAB自带的IFFT函数,输入补零后的数据S2,逐行进行Nr点快速傅里叶逆变换,完成对距离脉冲压缩后的距离向聚焦信号Src的距离向插值,得到插值后的信号S3,S3是大小为Na×(Nr×M1)的矩阵。
5.如权利要求4所述的斜视SAR时域成像自聚焦方法,其特征在于,所述步骤4,具体包括:
步骤401、设定地面成像网格尺寸:
根据步骤1中载入的分辨率指标ρr、ρa,设定在X方向地面网格的尺寸为Xp,Y方向地面网格的尺寸为Yp,X方向网格像素间隔为Δx,Y方向网格像素间隔为Δy,地面成像网格点数为AziNum×RanNum,具体由如下公式得到:
Xp=AziNum·Δx
Yp=RanNum·Δy (3)
步骤402、对X方向和Y方向的网格矩阵化:
初始设置X是大小为1×AziNum的矩阵,Y是大小为1×RanNum的矩阵,以场景中心点(xt,yt)作为网格的中心,m=1,2,…,AziNum,n=1,2,…,RanNum,根据公式(4)对网格进行离散化表示为:
X(1,m)=xt-XP+(2·Xp)×(m-1)/(RanNum-1)
Y(1,n)=yt-YP+(2·Yp)×(n-1)/(AziNum-1) (4)
步骤403、调用MATLAB中自带的repmat函数,输入步骤402生成的X网格,首先将X网格转置成为AziNum×1的矩阵,然后设置复写行数参数RanNum,输出大小为AziNum×RanNum的矩阵X(m,n);同理,调用repmat函数,输入步骤402生成的Y网格,设置复写行数AziNum,输出大小为AziNum×RanNum的矩阵Y(m,n);
步骤404、计算旋转角度θs,构造旋转矩阵H:
输入步骤1中加载的俯仰角参数β,其中β是大小为1×Na的向量、雷达三维位置坐标XT、YT、ZT,其中XT、YT、ZT均为大小为1×Na的向量以及公式(2)计算得到的斜距R1,利用俯仰角向量β的第Na/2点取值根据公式(5)计算得到中心时刻斜距在地面上的投影长度RL:
输入雷达平台在Y方向坐标YT的第Na/2点的取值YT(Na/2)以及中公式(5)得到的投影长度RL,根据正弦定理,计算旋转矩阵的旋转角度θs为:
根据输入的旋转角度θs计算旋转矩阵H为:
步骤405、输入步骤403中构造的成像网格X(m,n)、Y(m,n),以及旋转矩阵H,将X(m,n)按列拆分后纵向排列成一个列向量X1,此时X1的大小为(AziNum×RanNum)×1,同样将Y(m,n)按列拆分后纵向排列成一个列向量Y1,此时Y1的大小为(AziNum×RanNum)×1;将X1和Y1合并成一个新的矩阵并转置,得到转置后的矩阵C1,C1是大小为2×(AziNum×RanNum)的矩阵:
步骤406、首先将C1的每一列与(xt,yt)的转置相减得到矩阵C2,然后将旋转矩阵H与C2相乘,最后对矩阵相乘结果的每一列与(xt,yt)的转置相加,输出为旋转后的成像网格坐标C3,此时矩阵C3的大小为2×(AziNum×RanNum);
步骤407、调用MATLAB中reshape函数,输入参数分别为C3矩阵的第一行以及重排后的矩阵大小AziNum×RanNum,输出为重排后的网格矩阵坐标X2;同理调用MATLAB中reshape函数,输入参数分别为C3矩阵的第二行以及重排后的矩阵大小AziNum×RanNum,输出为重排后的网格矩阵坐标Y2:
6.如权利要求5所述的斜视SAR时域成像自聚焦方法,其特征在于,所述步骤5,具体包括:
步骤501、计算该采样点的雷达位置与成像网格中每一个像素点之间的斜距RT,RT为大小为AziNum×RanNum的矩阵;其中,
输入该采样时刻雷达的三维位置坐标XT、YT、ZT以及由步骤407得到的大小分别为AziNum×RanNum的地面成像网格坐标X2、Y2,根据公式(10)计算雷达平台与每一个网格的斜距,遍历所有网格点,得到每一个方位采样时刻雷达与所有网格点的斜距RT,RT是大小为AziNum×RanNum的矩阵:
步骤502、首先输入步骤501得到的斜距矩阵RT、光速c以及步骤302中的距离插值倍数M1,根据公式(11)、(12)得到tij和dtr;然后调用MATLAB中自带的round函数,根据公式(13)输入tij、dtr、c和R0,输出为Loc,Loc是大小为AziNum×RanNum的矩阵;
tij=2RT/c (11)
dtr=1/M1/Fr (12)
Loc=round((tij-2R0/c)/dtr+1) (13)
7.如权利要求6所述的斜视SAR时域成像自聚焦方法,其特征在于,在所述步骤6中,
SBP=S1+S2+...SNa (14)。
8.如权利要求6所述的斜视SAR时域成像自聚焦方法,其特征在于,所述步骤7,具体包括:
步骤701、输入BP聚焦数据SBP,调用MATLAB中的max求最大值函数,输出SBP中每一个列的最大值以及该最大值在每一列中的位置;
步骤702、输入SBP每一列数据以及该列数据中最大值的位置,调用MATLAB中自带的circshift函数,将每一列数据中的最大值移动到该列的中心位置,遍历所有列,得到移位后的数据S4,S4是大小为AziNum×RanNum的矩阵;
步骤703、输入步骤702中移位后的数据S4,将每一列数据以第AziNum/2个数据为中心,保留长度为W1的数据,其余数据置零,得到加窗处理后的数据S5:
步骤704、调用MATLAB自带的IFFT函数,输入加窗后的数据S5,对S5的每一列进行逆傅里叶变换得到S6;
步骤705、输入S6,遍历S6的AziNum行,对每一行数据取共轭,然后与下一行数据点乘得到S7,S7是大小为(AziNum-1)×RanNum的矩阵,然后将S7每一行中RanNum个数值相加得到S8,S8是大小为(AziNum-1)×1的矩阵;然后调用MATLAB中自带的angle函数,输入S8中每一个值,输出为相邻两行数据之间的相位差是大小为(AziNum-1)×1的矩阵;
步骤706、调用MATLAB中自带的cumsum函数,输入步骤705中得到的输出为各行数据与第一行数据之间的相位差,是大小为(AziNum-1)×1的矩阵;设置是大小为AziNum×1的矩阵,其中第一行的元素为0,其余(AziNum-1)行的元素为中的数值;
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110432393.3A CN113176570B (zh) | 2021-04-21 | 2021-04-21 | 一种斜视sar时域成像自聚焦方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110432393.3A CN113176570B (zh) | 2021-04-21 | 2021-04-21 | 一种斜视sar时域成像自聚焦方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN113176570A true CN113176570A (zh) | 2021-07-27 |
CN113176570B CN113176570B (zh) | 2022-11-15 |
Family
ID=76924039
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202110432393.3A Active CN113176570B (zh) | 2021-04-21 | 2021-04-21 | 一种斜视sar时域成像自聚焦方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN113176570B (zh) |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN118131167A (zh) * | 2024-05-08 | 2024-06-04 | 成都玖锦科技有限公司 | 一种近场sar成像中地面杂波的抑制方法 |
Citations (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103913741A (zh) * | 2014-03-18 | 2014-07-09 | 电子科技大学 | 一种合成孔径雷达高效自聚焦后向投影bp方法 |
CN104251990A (zh) * | 2014-09-15 | 2014-12-31 | 电子科技大学 | 合成孔径雷达自聚焦方法 |
CN104316923A (zh) * | 2014-10-14 | 2015-01-28 | 南京航空航天大学 | 针对合成孔径雷达bp成像的自聚焦方法 |
CN104931967A (zh) * | 2015-06-12 | 2015-09-23 | 西安电子科技大学 | 一种改进的高分辨率sar成像自聚焦方法 |
CN105842694A (zh) * | 2016-03-23 | 2016-08-10 | 中国电子科技集团公司第三十八研究所 | 一种基于ffbp sar成像的自聚焦方法 |
US20180196136A1 (en) * | 2017-01-11 | 2018-07-12 | Institute Of Electronics, Chinese Academy Of Sciences | Method and device for imaging by bistatic synthetic aperture radar |
CN110531338A (zh) * | 2019-10-12 | 2019-12-03 | 南京航空航天大学 | 基于fpga的多模式sar自聚焦快速处理方法及系统 |
CN110554385A (zh) * | 2019-07-02 | 2019-12-10 | 中国航空工业集团公司雷华电子技术研究所 | 机动轨迹合成孔径雷达自聚焦成像方法、装置及雷达系统 |
CN110632594A (zh) * | 2019-09-18 | 2019-12-31 | 北京航空航天大学 | 一种长波长星载sar成像方法 |
-
2021
- 2021-04-21 CN CN202110432393.3A patent/CN113176570B/zh active Active
Patent Citations (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103913741A (zh) * | 2014-03-18 | 2014-07-09 | 电子科技大学 | 一种合成孔径雷达高效自聚焦后向投影bp方法 |
CN104251990A (zh) * | 2014-09-15 | 2014-12-31 | 电子科技大学 | 合成孔径雷达自聚焦方法 |
CN104316923A (zh) * | 2014-10-14 | 2015-01-28 | 南京航空航天大学 | 针对合成孔径雷达bp成像的自聚焦方法 |
CN104931967A (zh) * | 2015-06-12 | 2015-09-23 | 西安电子科技大学 | 一种改进的高分辨率sar成像自聚焦方法 |
CN105842694A (zh) * | 2016-03-23 | 2016-08-10 | 中国电子科技集团公司第三十八研究所 | 一种基于ffbp sar成像的自聚焦方法 |
US20180196136A1 (en) * | 2017-01-11 | 2018-07-12 | Institute Of Electronics, Chinese Academy Of Sciences | Method and device for imaging by bistatic synthetic aperture radar |
CN110554385A (zh) * | 2019-07-02 | 2019-12-10 | 中国航空工业集团公司雷华电子技术研究所 | 机动轨迹合成孔径雷达自聚焦成像方法、装置及雷达系统 |
CN110632594A (zh) * | 2019-09-18 | 2019-12-31 | 北京航空航天大学 | 一种长波长星载sar成像方法 |
CN110531338A (zh) * | 2019-10-12 | 2019-12-03 | 南京航空航天大学 | 基于fpga的多模式sar自聚焦快速处理方法及系统 |
Non-Patent Citations (3)
Title |
---|
李根 等: "基于Keystone变换和扰动重采样的机动平台大斜视SAR成像方法", 《电子与信息学报》 * |
王昕 等: "一种改进的SAR反投影图像相位梯度自聚焦方法", 《现代雷达》 * |
王昕 等: "机载SAR反投影图像自聚焦处理方法", 《航空学报》 * |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN118131167A (zh) * | 2024-05-08 | 2024-06-04 | 成都玖锦科技有限公司 | 一种近场sar成像中地面杂波的抑制方法 |
Also Published As
Publication number | Publication date |
---|---|
CN113176570B (zh) | 2022-11-15 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
US7397418B1 (en) | SAR image formation with azimuth interpolation after azimuth transform | |
CN107229048A (zh) | 一种高分宽幅sar动目标速度估计与成像方法 | |
CN110865344B (zh) | 一种脉冲多普勒雷达体制下副瓣快速抑制方法 | |
CN105093223A (zh) | 双站前视sar的快速时域成像方法 | |
CN110673143A (zh) | 一种子孔径大斜视sar俯冲成像的两步处理方法 | |
CN106405552A (zh) | 基于wvd—pga算法的sar雷达目标聚焦方法 | |
CN111487614B (zh) | 基于子孔径的曲线航迹弹载sar波前重建成像方法及系统 | |
CN110954899B (zh) | 高海况下海面舰船目标成像方法及装置 | |
CN113030964B (zh) | 基于复Laplace先验的双基地ISAR稀孔径高分辨成像方法 | |
CN113176570B (zh) | 一种斜视sar时域成像自聚焦方法 | |
CN111551934A (zh) | 一种用于无人机载sar成像的运动补偿自聚焦方法与装置 | |
CN112433210A (zh) | 一种双站前视探地雷达快速时域成像方法 | |
CN113759372A (zh) | 弹载大斜视小孔径多通道sar的成像方法 | |
CN105044720A (zh) | 一种基于直角坐标系的后向投影成像方法 | |
CN113484862A (zh) | 一种自适应的高分宽幅sar清晰重构成像方法 | |
CN117556605A (zh) | 一种多体制雷达仿真系统及其控制方法 | |
CN110618409B (zh) | 顾及叠掩及阴影的多通道InSAR干涉图仿真方法及系统 | |
CN113030963B (zh) | 联合残余相位消除的双基地isar稀疏高分辨成像方法 | |
CN115792906A (zh) | 一种星载大斜视滑动聚束sar成像处理方法 | |
CN110297242A (zh) | 基于压缩感知的合成孔径雷达层析三维成像方法及装置 | |
CN115267706A (zh) | 合成孔径雷达的距离空变相位误差估计方法、装置及介质 | |
CN113885026A (zh) | 运动目标的sar稀疏成像方法、装置、电子设备及存储介质 | |
Li et al. | High squint multichannel SAR imaging algorithm for high speed maneuvering platforms with small-aperture | |
Garren et al. | Image preconditioning for a SAR image reconstruction algorithm for multipath scattering | |
CN115291211B (zh) | 一种无定标器的圆迹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 | ||
EE01 | Entry into force of recordation of patent licensing contract | ||
EE01 | Entry into force of recordation of patent licensing contract |
Application publication date: 20210727 Assignee: SPACETY Co.,Ltd. (CHANGSHA) Assignor: BEIHANG University Contract record no.: X2023990000783 Denomination of invention: A self focusing method for time-domain imaging of squint SAR Granted publication date: 20221115 License type: Common License Record date: 20230829 |