一种飞行器SAR实时成像方法与流程

一种飞行器sar实时成像方法
技术领域
1.本发明属于雷达信号处理技术领域，特别涉及一种飞行器sar实时成像方法。


背景技术：

2.与其他探测感知系统相比，sar成像系统具有全天候、全天时、高分辨等诸多优势，可通过图像匹配实现精确制导与飞行器定位。不同于机载、星载平台，机动性较高的飞行器如无人机、导弹等，其应用场景决定了它高速高机动、非线性航迹等突出特点，导致雷达录取回波存在几何形变、距离徙动大、多普勒参数变化大等问题，这对成像算法提出了很高的误差补偿要求与实时性要求。
3.目前，常采用子孔径方法解决机动平台sar成像难题，核心思想是通过减小合成孔径时长来降低数据处理的复杂度，近似认为飞行器在子孔径内做匀加速运动。常用的子孔径成像算法主要有specan(spectral analysis)算法、ecs(extended chirp scaling)算法等。specan算法的优点是成像效率高、易于硬件实现，但该算法基于简单模型，无法校正距离徙动的空变量，成像分辨率低；ecs算法，能实现空变距离徙动的精确校正，但方位向调频变标的补零点数多，硬件实现困难、实时性差。因此，适用于机动平台实时匹配导航的简洁、快速、高效sar成像算法极具实际研究价值。


技术实现要素：

4.本发明的目的是解决现有技术中存在的问题，提供一种飞行器sar实时成像方法。
5.为了实现以上目的，本发明通过以下技术方案实现：
6.一种飞行器sar实时成像方法，采用间歇burst成像模式，包括：
7.步骤s1，结合惯导信息，对单个burst回波数据做多普勒中心频率估计；
8.步骤s2，结合该burst回波数据及其多普勒中心频率估计值，进行多普勒调频率估计，输出飞行器的真实速度与径向加速度；
9.步骤s3，利用输出的径向加速度a
r
生成相位补偿函数h
acc
，将该burst回波数据s0(t
r
,t
a
)乘上h
acc
完成飞行器加速度误差补偿，得到信号s1(t
r
,t
a
)；
10.步骤s4，基于信号s1(t
r
,t
a
)，采用改进ecs算法完成回波数据的sar成像预处理；
11.步骤s5，沿方位向划分子孔径数据，分别做方位fft得到多幅子孔径图像，经多视处理得到burst图像；
12.步骤s6，对多幅burst图像进行sar定位、配准、拼接，得到大场景sar图像。
13.优选的，步骤s1进一步包括：
14.s11，利用单个burst内的惯导数据计算飞行器速度大小均值v
g
和波束中心斜视角依据几何法计算多普勒中心频率f
d1
；
15.s12，利用基于数据的相关函数法计算多普勒中心频率基带值f
d2
；
16.s13，结合f
d1
、f
d2
进行多普勒中心解模糊，获得精确的绝对多普勒中心频率估计值f
dc
。
17.优选的，步骤s2进一步包括：
18.s21，完成方位频谱基带搬移与方位加窗，方位频谱搬移在方位时域进行，将基带解调的回波信号s0(t
r
,t
a
)乘上方位频谱基带搬移与方位加窗函数之乘积h0，得到方位基带回波信号s
01
(t
r
,t
a
)；
19.s22，对方位基带回波信号s
01
(t
r
,t
a
)做距离向fft，再与距离向频域匹配滤波函数h
r
相乘，得到信号s
02
(f
r
,t
a
)；
20.s23，对信号s
02
(f
r
,t
a
)做方位向fft，进行多普勒带宽的截取得到信号s
03
(f
r
,f
a
)；
21.s24，结合burst内的惯导加速度矢量，计算飞行器加速度在波束指向上的径向加速度分量a
r0
，并选取距离中心临近的δn个距离单元，计算各距离单元多普勒调频率初值k
a0
(n
r
)；
22.s25，利用中心距离单元的多普勒调频率计算飞行器等效飞行速度v
eq
，再生成距离徙动校正函数h
rcmc
，将信号s
03
(f
r
,f
a
)乘上h
rcmc
完成距离徙动校正，并回到两维时域得到信号s
04
(t
r
,t
a
)；
23.s26，利用多普勒调频率构造二次相位函数s
ref
对两维时域信号s
04
(t
r
,t
a
)进行方位dechirp，得到信号s
05
(t
r
,t
a
)，然后在方位时域划分两视子孔径数据，各自做方位fft生成两视子孔径图像s
sub1
和s
sub2
，通过最大相关法估计两视图像的相对移动量δ，进而计算得到多普勒调频率误差δk
a
；
24.s27，更新飞行器速度，判定速度估计收敛条件：若收敛，则输出飞行器速度、斜视角、径向加速度，否则更新中心距离单元的多普勒调频率，返回执行s25及之后的所有步骤，直至满足迭代终止条件，输出速度、斜视角、径向加速度估计结果v、a
r
。
25.优选的，在步骤s24中，第n
r
个距离单元多普勒调频率初值k
a0
(n
r
)的计算表达式为：
[0026][0027]
式中，v
g
为由该burst内的惯导数据计算的飞行器速度大小均值，表示雷达波束中心斜视角初值，λ为雷达工作波长，r(n
r
)表示第n
r
个距离单元对应的斜距大小，计算式为
[0028][0029]
r
min
＝cτ0/2为雷达回波最短延时τ0对应的斜距大小，δr＝c/(2f
s
)为距离采样间隔，f
s
为距离采样频率，t
p
为发射信号脉冲宽度，n
r
为回波距离向采样点数，c为光速。
[0030]
优选的，飞行器等效飞行速度v
eq
的计算表达式如下：
[0031][0032]
其中，k
a0
(n
r
/2+1)表示中心距离单元的多普勒调频率，r
s
＝r
min
+n
r
·
δr/2
‑
ct
p
/4。
[0033]
优选的，距离徙动校正函数h
rcmc
的计算表达式如下：
[0034][0035]
其中，f
r
为距离频率，为最短距离；分别表示徙动参数和参考徙动参数，式中的v初次迭代采用v
g
计算，之后均采用前一次迭代更新的等效飞行速度v
eq
计算，f
a
为方位频率，f
dc
为多普勒中心频率估计值。
[0036]
优选的，二次相位函数的计算表达式如下：
[0037][0038]
其中，δk
a
为多普勒调频率误差，初始值为0，t
a
为方位向慢时间。
[0039]
优选的，多普勒调频率误差δk
a
的计算表达式为：
[0040][0041]
式中，prf表示雷达的脉冲重复频率，n
a
表示多普勒中心频率估计所用回波数据的方位单元数，δ为通过最大相关法估计的两视子孔径图像的相对移动量。
[0042]
优选的，在步骤s3中，相位补偿函数为
[0043]
优选的，步骤s4进一步包括：
[0044]
s41，将加速度误差补偿完的回波数据s1(t
r
,t
a
)乘以相位函数h0进行方位频谱搬移与加窗，得到信号s2(t
r
,t
a
)；
[0045]
s42，将信号s2(t
r
,t
a
)经方位向fft得到距离多普勒信号s2(t
r
,f
a
)，再乘以线性调频变标方程h1(t
r
,f
a
)完成线性调频变标；
[0046]
s43，将信号s2(t
r
,f
a
)做距离向fft，乘以相位补偿函数h2(f
r
,f
a
)完成距离脉压与距离徙动校正，经距离向ifft得到距离多普勒信号s3(t
r
,f
a
)；
[0047]
s44，将s3(t
r
,f
a
)乘以相位补偿函数h3(t
r
,f
a
)，以补偿线性调频变标引入的附加相位，再乘以方位向变标的相位补偿函数h4(t
r
,f
a
)，将不同斜距处目标的方位向调频率校正为一致的，再沿距离向乘以扇贝效应校正函数h5(n
a
)完成方位天线图导致的扇贝效应；
[0048]
s45，对上一步的信号做方位向ifft回到两维时域，乘以线性相位函数h6(t
a
)完成方位频谱搬移复原，再乘以相位补偿函数h7(t
r
,t
a
)完成方位向dechirp。
[0049]
与现有技术相比，本发明具有如下优点：
[0050]
(1)针对间歇burst成像体制，避免了传统ecs成像算法方位调频变标处理的大量补零操作，对平台存储资源的需求更低，实时性好；
[0051]
(2)针对飞行器匹配导航的应用需求与条件限制，兼顾算法效率与成像精度，保相性好。
附图说明
[0052]
为了更清楚地说明本发明的技术方案，下面将对描述中所需要使用的附图作简单
地介绍，显而易见地，下面描述中的附图是本发明的一个实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图：
[0053]
图1是本发明方法的间歇bust工作体制示意图；
[0054]
图2是sar成像几何关系图；
[0055]
图3是实时sar成像算法的流程框图；
具体实施方式
[0056]
以下结合附图和具体实施方式对本发明提出的方案作进一步详细说明。根据下面说明，本发明的优点和特征将更清楚。需要说明的是，附图采用非常简化的形式且均使用非精准的比例，仅用以方便、明晰地辅助说明本发明实施方式的目的。为了使本发明的目的、特征和优点能够更加明显易懂，请参阅附图。须知，本说明书所附图式所绘示的结构、比例、大小等，均仅用以配合说明书所揭示的内容，以供熟悉此技术的人士了解与阅读，并非用以限定本发明实施的限定条件，故不具技术上的实质意义，任何结构的修饰、比例关系的改变或大小的调整，在不影响本发明所能产生的功效及所能达成的目的下，均应仍落在本发明所揭示的技术内容能涵盖的范围内。
[0057]
参照图1，本发明方法采用间歇burst体制的sar成像模式，系统工作模式为雷达开机工作一段时间后(小于合成孔径时间)待机，待机期间进行信号处理，然后雷达再次开机，如此循环多次，对录取的多个burst数据单独进行sar成像处理，最后对多个burst图像进行拼接获取满足指标要求的sar图像。
[0058]
图2给出了sar成像几何关系，o
‑
xyz为匹配中心航向坐标系，原点o为dem匹配区域中心点，x轴水平指向预定航向方向，y轴垂直向天(地表法线方向)指天为正，z轴由右手系确定；h表示飞行器高度，v表示飞行器速度，r
s
表示雷达斜距，r0表示最短距离，θ表示雷达下视角，θ0表示零多普勒面下视角，表示雷达波束中心斜视角，φ表示雷达波束中心方位角。假设burst时间内，飞行器近似做匀加速运动，即飞行器加速度可以描述为恒定加速度a，则多普勒调频率的计算式为
[0059][0060]
其中，λ为雷达工作波长，a
r
为矢量a在雷达视线反方向上的投影分量。可见，恒加速模型下，sar回波多普勒调频率会多出与径向加速度有关的一项。
[0061]
据此，本发明所提出的一种飞行器sar实时成像方法包括以下步骤：
[0062]
步骤s1，对单个burst的回波数据做多普勒中心估计，具体如下：
[0063]
s11斜视情况下，多普勒中心频率会发生模糊，首先利用惯导提供的飞行器速度和斜视角，依据几何法计算多普勒中心频率
[0064][0065]
其中，v
g
为由该burst内的惯导数据计算的飞行器速度大小均值，为惯导提供的波束中心斜视角。
[0066]
s12采用基于数据的相关函数法计算多普勒中心频率基带值
[0067][0068]
其中，n表示多普勒中心估计所用回波数据的距离单元数，n
a
表示其方位单元数，prt表示雷达的脉冲重复周期，arg(
·
)表示取相角函数，s0(n
r
,n
a
)表示回波中距离、方位单元索引分别为n
r
、n
a
的复信号，上标*表示求共轭运算。
[0069]
s13结合两种方法的估计结果f
d1
、f
d2
进行多普勒中心解模糊，得到精确的绝对多普勒中心频率估计值
[0070]
f
dc
＝round[(f
d1
‑
f
d2
)/prf]
×
prf+f
d2
[0071]
其中，round(
·
)表示四舍五入取整函数，prf＝1/prt为脉冲重复频率。
[0072]
步骤s2，结合该burst回波数据与多普勒中心频率估计值，估计其多普勒调频率，输出飞行器的真实速度与径向加速度，具体如下：
[0073]
(2a)利用多普勒中心频率估计值将回波方位频谱搬移到基带，并完成距离脉压；
[0074]
单个点目标的基带回波信号表达式为
[0075][0076]
其中，t
r
为距离向快时间，t
a
为方位向慢时间，a0为一个复常数，后续推导予以省略，ω
r
(
·
)为发射脉冲包络，r
s
(t
a
)为雷达至目标的瞬时斜距，c为光速，ω
a
(
·
)为方位向双程天线方向图，t
ac
为雷达波束中心穿越时刻，f
c
为雷达载频，k
r
为发射脉冲线性调频率(距离调频率)。
[0077]
方位频谱搬移在方位时域进行，得到方位基带回波信号
[0078]
s
01
(t
r
,t
a
)＝s0(t
r
,t
a
)
·
h0[0079]
其中，h0＝h
00
·
h
01
，h
00
＝exp(
‑
j2πf
dc
t
a
)，h
00
为方位频谱搬移函数，h
01
＝(1
‑
α)
‑
α
·
cos(2π(n
a
‑
1)/n
a
),1≤n
a
≤n
a
，h
01
为hamming窗函数(锐化窗)，窗系数α通常取0.3。
[0080]
完成方位频谱基带搬移后，需在距离频域完成脉冲压缩处理
[0081]
s
02
(f
r
,t
a
)＝fft
r
[s
01
(t
r
,t
a
)]
·
h
r
[0082]
其中，f
r
为距离频率，fft
r
(
·
)表示距离向fft；h
r
＝rect{f
r
/b
r
}exp(jπf
r2
/k
r
)为距离频域匹配滤波函数，rect(
·
)为矩形窗函数，b
r
为发射信号带宽，exp(
·
)表示e指数函数。
[0083]
(2b)方位频谱截取，以减小斜距模型近似误差；
[0084]
由于map drift(md)算法默认的斜距模型为抛物线，与实际的双曲线相比存在近似误差，且该误差随着方位波束宽度的增大而增大，根据实际方位波束宽度情况进行加窗截取
[0085][0086]
其中，f
a
为方位频率，fft
a
(
·
)表示方位向fft，|
·
|表示求绝对值；为波束宽度θ
limit
对应的多普勒带宽，是采用多普勒中心频率估计值f
dc
和v
g
计算的斜视角初值，θ
limit
＝θ
a
/m为缩小后的方位波束宽度(θ
a
表示方位波束原始宽度，m为波束宽度缩小倍数，取正整数)。
[0087]
(2c)选取距离中心临近的δn个距离单元，计算相应的多普勒调频率初值，同时利用中心距离单元的多普勒调频率初值计算飞行器等效速度，完成距离徙动校正；
[0088]
调频率初值由惯组提供的速度、加速度信息，结合斜视角初值、多普勒中心频率估计模块输出的波束中心指向向量信息获得。第n
r
个距离单元的多普勒调频率为
[0089][0090]
其中，r(n
r
)表示第n
r
个距离单元对应的斜距大小，计算式为(n
r
‑
δn)/2+1≤n
r
≤(n
r
+δn)/2，式中，r
min
＝cτ0/2为雷达回波最短延时τ0对应的斜距大小，δr＝c/(2f
s
)为距离采样间隔(f
s
为距离采样频率)，t
p
为发射信号脉冲宽度，n
r
为回波距离向采样点数。加速度在波束中心指向上的投影量a
r0
初始计算式为
[0091]
a
r0
＝a
·
d
beam
[0092]
其中，a＝[a
xg
,a
yg
,a
zg
]为惯导加速度矢量，a
xg
,a
yg
,a
zg
分别为该burst内由惯导给出的飞行器在三个坐标轴方向的各自加速度均值，是根据几何关系计算的波束中心指向反方向单位向量的初始值。
[0093]
利用距离单元中心处的多普勒调频率计算飞行器等效速度
[0094][0095]
其中，r
s
＝r
min
+n
r
·
δr/2
‑
ct
p
/4。
[0096]
利用上面计算的等效速度生成距离徙动校正函数
[0097][0098]
其中，为最短距离；
分别表示徙动参数和参考徙动参数，式中的v初次迭代采用v
g
计算，之后均采用前一次迭代更新的等效速度v
eq
计算。
[0099]
对方位谱截取信号s
03
(f
r
,f
a
)完成距离徙动校正，并回到两维时域
[0100]
s
04
(t
r
,t
a
)＝ifft2[s
03
(f
r
,f
a
)
·
h
rcmc
]
[0101]
其中，ifft2(
·
)为两维逆fft操作。
[0102]
(2d)通过md算法估计当前多普勒调频率的误差；
[0103]
利用多普勒调频率构造二次相位函数对两维时域信号s
04
(t
r
,t
a
)进行方位dechirp，二次相位函数表达式为
[0104][0105]
其中，δk
a
为多普勒调频率误差，初始值为0，即s
ref
初次采用计算的多普勒调频率初值生成，后续采用前一次迭代的多普勒调频率更新值生成。
[0106]
s
05
(t
r
,t
a
)＝s
04
(t
r
,t
a
)
·
s
ref
[0107]
去斜后，在时域划分两视子孔径数据，经方位fft生成两视子孔径图像s
sub1
和s
sub2
[0108]
s
sub1
＝|fft
a
[s
05
(t
r
,
‑
t
a
/2≤t
a
≤0)]|
[0109]
s
sub2
＝|fft
a
[s
05
(t
r
,0≤t
a
≤t
a
/2)]|
[0110]
其中，t
a
是该burst时长。
[0111]
通过最大相关法估计两视图像的相对移动量δ，得到多普勒调频率误差量
[0112][0113]
(2e)更新飞行器飞行速度，判定速度估计收敛条件：若收敛，则输出飞行器速度、斜视角、径向加速度，否则更新多普勒调频率，回到(2c)重新往下执行；
[0114]
先更新多普勒调频率
[0115]
k
a
＝k
a0
(n
r
/2+1)+δk
a
[0116]
从多普勒调频率结果中去除飞行器加速度对多普勒调频率影响的部分，计算飞行器实际飞行速度
[0117][0118]
同时更新斜视角和径向加速度
[0119][0120][0121]
若速度估计值已收敛，则输出速度、斜视角、径向加速度估计结果；若速度估计值未收敛，则回到步骤(2c)，重新计算飞行器等效速度
[0122][0123]
重新执行步骤(2c)
‑
(2e)，直至达到迭代终止条件，输出速度、斜视角、径向加速度估计结果v、a
r
。
[0124]
步骤s3，利用多普勒参数估计结果，对该burst回波数据做加速度误差补偿；
[0125]
加速度误差补偿的目的是通过相位相乘补偿由于飞行器平台恒加速运动导致的相位误差
[0126]
s1(t
r
,t
a
)＝s0(t
r
,t
a
)
·
h
acc
[0127]
其中，相位补偿函数为
[0128]
步骤s4，采用改进的ecs算法完成回波数据的成像预处理；
[0129]
(4a)首先，同步骤(2a)，将加速度误差补偿完的回波数据乘以h0进行方位频谱搬移与方位加窗
[0130]
s2(t
r
,t
a
)＝s1(t
r
,t
a
)
·
h0[0131]
(4b)方位向fft，进行线性调频变标
[0132]
对s2(t
r
,t
a
)进行方位向fft，得到的距离多普勒域信号为
[0133][0134]
其中，w
a
(
·
)为方位向双程天线方向图的频谱，k
m
为距离/方位耦合后的距离向调频率
[0135][0136]
式中的f
c
为雷达载频。
[0137]
线性调频变标方程式为
[0138][0139]
式中，参考方位频率f
aref
＝f
dc
，r0(n
r
,n
a
)为第n
r
个距离单元、第n
a
个方位单元对应
信号的最短距离(不同于常数r0)，计算式为1≤n
r
≤n
r
,1≤n
a
≤n
a
，对于每一个采样点，均需要将其对应的斜距r(n
r
)和多普勒频率f
a
(n
a
)代入上式计算得到r0(n
r
,n
a
)。
[0140]
参考距离r
ref
由雷达斜距r
s
计算，表达式为
[0141][0142]
此处距离向尚未做脉压，斜距计算均需补偿半个脉宽对应的斜距量。
[0143]
将s2(t
r
,f
a
)与h1(t
r
,f
a
)相乘完成线性调频变标操作。
[0144]
(4c)在距离频域完成距离脉压和距离徙动校正
[0145]
对线性调频变标后的信号做距离向fft变换到两维频域，与相位补偿函数相乘，完成距离压缩、二次距离压缩和一致rcmc，相位补偿函数表达式为
[0146][0147]
再经距离向ifft后，得到距离多普勒域信号
[0148][0149]
其中，sinc(
·
)表示辛克函数。
[0150]
(4d)附加相位补偿
[0151]
该步骤用于补偿线性调频变标引入的附加相位，相位补偿函数为
[0152][0153]
由于距离脉压后，脉冲峰值位于脉冲中心，因此距离向每个采样点的斜距计算式更新成
[0154]
r
′0＝r(n
r
)
·
d(f
aref
,v),1≤n
r
≤n
r
[0155]
r
′
ref
＝r
s
·
d(f
aref
,v)
[0156]
(4e)方位向变标
[0157]
方位向变标目的是通过相位相乘将不同斜距处目标的方位向调频率校正为一致，
且为线性调频率信号。相位补偿函数为
[0158][0159]
其中，为参考多普勒调频率，t
aref
(r
′
ref
)＝f
dc
/k
aref
；t
a
(r0′
)＝
‑
f
dc
λr
′0/(2v2d(f
dc
,v))。
[0160]
(4f)校正扇贝效应
[0161]
该步校正由于方位向天线方向图导致的扇贝效应。通过对不同方位向的采样点乘以一个天线方向图平方的倒数来实现。扇贝效应校正函数为
[0162][0163]
其中，
[0164][0165]
然后，经方位向ifft回到方位时域
[0166][0167]
(4g)方位频谱搬移复原，进行方位向dechirp
[0168]
方位向频谱搬移的目的是避免方位向成像后图像混叠，该操作通过线性相位相乘实现。线性相位函数表达式为
[0169]
h6(t
a
)＝exp{j2πf
dc
′
t
a
}
[0170]
其中，为多普勒中心频率更新值。
[0171]
方位向dechirp通过相位相乘补偿方位向信号的二次相位，相位补偿函数为
[0172][0173]
补偿后的信号为
[0174]
s5(t
r
,t
a
)＝s4(t
r
,t
a
)
·
h6·
h7。
[0175]
步骤s5，沿方位向划分子孔径数据，分别做方位fft得到多幅子孔径图像，经多视处理得到burst图像；
[0176][0177]
其中，m
s
为子孔径图像数，i表示子孔径图像编号，s
5,i
(t
r
,t
a
)表示第i幅子孔径图像。
[0178]
步骤s6，对多幅burst图像进行sar定位、配准、拼接，得到大场景sar图像，可供飞行器匹配导航使用。
[0179]
尽管本发明的内容已经通过上述优选实施例作了详细介绍，但应当认识到上述的描述不应被认为是对本发明的限制。在本领域技术人员阅读了上述内容后，对于本发明的多种修改和替代都将是显而易见的。因此，本发明的保护范围应由所附的权利要求来限定。