CN106802416B - 一种快速因式分解后向投影sar自聚焦方法 - Google Patents

Abstract

本发明公开了一种快速因式分解后向投影SAR自聚焦方法，结合快速因式分解后向投影成像算法(FFBP)建立相位误差优化模型，并利用坐标下降方法和割线方法来求解该优化模型，有效的解决了相位误差补偿问题。与现有技术相比，本发明能够更加精确地估计相位误差以及得到良好聚焦的SAR图像，解决了现有FFBP自聚焦方法的运算量高、精度低、适应小的问题，从而实现结合FFBP成像的精确运动误差估计和良好聚焦。

Description

一种快速因式分解后向投影SAR自聚焦方法

技术领域

本发明属于雷达技术领域，具体涉及一种基于最大图像锐度的快速因式分解后向投影SAR自聚焦方法的设计。

背景技术

合成孔径雷达(Synthetic Aperture Radar，SAR)是一种具有高分辨率的成像雷达，与光学传感器相比，SAR具有全天时全天候工作能力的独特优点。随着SAR技术的发展和提高，其分辨率越来越高，目前已接近或超光学成像的分辨率，因而被广泛的应用于地球遥感、海洋研究、资源勘探、灾情预报和军事侦察等领域。

基于时域的快速因式分解后向投影(FFBP)成像方法不存在近似处理，具有高精度成像能力和良好的相位保持特性，几乎适用于所有合成孔径雷达系统，包括双/多基地合成孔径雷达等具有复杂成像几何的合成孔径雷达系统；并且解决了传统后向投影(BP)成像方法的大的计算复杂度问题。

对于SAR的高分辨成像，自聚焦是关键步骤之一。FFBP成像算法的应用已经逐渐趋于成熟，但是对于FFBP的自聚焦还有待进一步的研究。传统的自聚焦方法，如相位梯度自聚焦、子视图相关等都是建立在图像域和距离压缩相位历史域的傅里叶变换关系上；然而这种变换关系对于FFBP却是不成立的。为了解决这个问题，文献Zhang Lei,Li Haolin,XingMengdao,Bao Zheng，“Integrating Autofocus Techniques With Fast FactorizedBack-projection For High-Resolution Spotlight SAR Imaging,”IEEE Geoscienceand remote sensingpp:1394-1398，2013引入虚拟极坐标系，使二者的傅里叶变换关系近似成立，同时对误差函数进行多项式拟合，提出子孔径图像偏移自聚焦，但是却受限于实际大场景且对高频误差进行估计时性能下降；文献Jan Torgrimsson,Patrik Dammert,HansHellsten,Lars M.H.Ulander，“Factorized Geometrical Autofocus for SyntheticApertuer Radar Processing,”IEEE Transactions on geoscience and remotesensing，2014和文献Jan Torgrimsson,Lars M.H Ulander,Patrik Dammert,HansHellsten，“Factorized geometrical autofocus:On the geometry search,”IEEETransactions on geoscience and remote sensing，2016从几何参数出发，建立评价图像聚焦性能的代价函数，利用优化求解的方法对几何参数进行估计，但是涉及多维参数优化求解，运算量大，速度慢，且当涉及复杂运动时估计精度降低，适用性小。

发明内容

本发明的目的是为了解决现有技术中针对快速因式分解后向投影成像方法中的运动误差补偿问题，提出了一种基于最大图像锐度的快速因式分解后向投影SAR自聚焦方法。

本发明的技术方案为：一种快速因式分解后向投影SAR自聚焦方法，包括以下步骤：

S0、系统参数初始化；

S1、快速因式分解后向投影成像；

S2、相位误差估计。

进一步地，步骤S0具体为：

以场景参考目标点O＝[0 0 0]^T为原点建立坐标系，场景中任意像素点位置记为P＝[x y 0]^T，雷达平台理想位置记为P_A＝[0 vt h]^T；

其中，v为载机沿着y轴的理想飞行速度，h为载机沿着y轴的理想飞行高度，t为方位慢时间；

任意像素点距雷达平台的理想瞬时距离记为R(t)＝||P-P_A||；当飞行轨迹存在偏差时，雷达平台实际位置记为

其中[Δx(t) Δy(t) Δz(t)]^T表示载机位置偏差；

任意像素点距雷达平台的实际瞬时距离记为

即：

将公式(1)展开可以得到：

其中ΔR表示距离误差；得到对应相位误差表达式：

式中λ表示雷达发射信号对应波长。

进一步地，步骤S1具体为：

根据因式分解的原理确定最佳的初始孔径长度l₀以及每次合并的子图像个数I，即分解因子；

以各个初始子孔径的中心为原点建立极坐标系，划分出子图像的像素点坐标(r,θ)的取值范围,有

其中，R_srn为近端距离，R_srf为远端距离，r表示像素点到子孔径坐标原点的距离，θ表示孔径矢量和距离矢量r之间的夹角，Θ表示积累角；极坐标系下的瞬时距离表示为：

其中，d表示雷达平台到对应坐标系原点的距离；

确定初始子图像的角度向分辨率 同时根据信号带宽确定距离向分辨率ρr；方 位向采样点个数为N，初始成像子孔径个数为K；通过公式(5)(6)对成像场景中的每一点(r, θ)进行初始成像：

其中，w(θ)为天线方向函数，s(x,τ)为获取的回波，δ(·)为单位冲激函数，x为雷达在方位向的位置，τ为对应的快时间，下标q表示方位向采样点，k表示第k个子孔径，且有q∈[1,N],k∈[1,K]，上标(1)表示第一级(初始)成像；

考虑确定雷达位置对应确定的相位误差，记飞行轨迹偏差导致的方位向相位误差为φ＝{φ₁,φ₂.....φ_N}，代入公式(6)，得到

其中，j表示为虚数单位；则初始成像更新为：

对于i+1级某个成像网格中的某点(r,θ)，计算其在第i级第k幅子图像中的位置(r′,θ′)，将第k幅图像在该点的结果采用插值的方式累加至新的i+1级图像f(r,θ)⁽ⁱ⁺¹⁾中，根据公式(8)-(10)对子孔径成像结果进行逐级合并：

rcos(θ)-r′cosθ′＝d_k (8)

此时，d_k表示为当前子孔径对应的上一级合并时第k个子孔径距离该子孔径几何中心的距离；

每进行一级合并，将I幅子图像生成次一级子图像，第i+1幅子图像分辨率和第i级子图像分辨率存在以下的关系：

其中，ρ_θm为最终分辨率；得到最终图像：

其中，

表示多次合并插值之后的值，

表示第q个方位向采样点相位误差的估计值，用于补偿相位误差，则f^(F)为最终补偿后的图像。

进一步地，步骤S2具体为：

记图像锐度为：

其中i表示第i个像素点，记最终得到的图像f^(F)→f，对应的共轭表示为f^*，则v_i＝f_if_i ^*；建立带有未知相位误差的最大图像锐度的优化模型来求解相位误差

即是对φ的估计，最终使得图像良好聚焦，表示为：

使用坐标下降的迭代方法来求解最优解；在坐标下降的算法中，每一次迭代都使得优化变量按序得以更新，而保持其他变量固定不变；记

为第i次迭代的第n个相位误差校正量，则第i+1次迭代表示为：

最终图像表示如下：

其中，x表示除了第n次采样点的其他所有采样点的后向投影之和，y表示对当前采样点校正后的后向投影；将公式(16)代入(13)得到：

其中，x_i和y_i为第i个像素点在公式(16)下的对应表示；x_i ^*和y_i ^*为x_i和y_i对应的共轭表示，

因此在坐标下降处理下将多维优化求解简化为单变量求极值问题

由公式(13)得到关于φ的导数：

将公式(19)(20)代入公式(18)最终得到：

其中：

利用割线法对单变量优化问题进行求解，具体步骤如下：

T1、选择初始点(φ₀)¹，参数ε＞0,γ＞0，令(φ₁)¹＝(φ₀)¹-γφ′((φ₀)¹)，k:＝0；

T2、若|φ′((φ_k)ⁱ)|＞ε，则利用公式(23)和(24)更新α_k，否则停止运算；其中，上标i表示第i个方位向采样点；

T3、通过公式(25)和(26)计算(φ_k+1)ⁱ,k＝k+1，返回步骤T2；

公式(23)-(26)中的下标k表示对于第i个方位向采样点的第k次迭代求解；

通过上述步骤，实现对相位误差进行精确估计并实现良好成像。

本发明的有益效果是：本发明结合快速因式分解后向投影成像算法(FFBP)建立相位误差优化模型，并利用坐标下降方法和割线方法来求解该优化模型，有效的解决了相位误差补偿问题。与现有技术相比，本发明能够更加精确地估计相位误差以及得到良好聚焦的SAR图像，解决了现有FFBP自聚焦方法的运算量高、精度低、适应小的问题，从而实现结合FFBP成像的精确运动误差估计和良好聚焦。

附图说明

图1为本发明提供的一种快速因式分解后向投影SAR自聚焦方法流程图。

图2为本发明实施例采用的快速因式分解后向投影成像具体实现结构示意图。

图3为本发明实施例采用的雷达系统结构示意图。

图4为本发明实施例的场景点目标分布图。

图5为本发明实施例的成像散焦结果图。

图6为本发明实施例的相位误差估计图。

图7为本发明实施例的成像聚焦结果图。

具体实施方式

下面结合附图对本发明的实施例作进一步的说明。

本发明提供了一种基于最大图像锐度的快速因式分解后向投影SAR自聚焦方法，如图1所示，包括以下步骤：

S0、系统参数初始化。

以场景参考目标点O＝[0 0 0]^T为原点建立坐标系，场景中任意像素点位置记为P＝[x y 0]^T，雷达平台理想位置记为P_A＝[0 vt h]^T。

其中，v为载机沿着y轴的理想飞行速度，h为载机沿着y轴的理想飞行高度，t为方位慢时间。

任意像素点距雷达平台的理想瞬时距离记为R(t)＝||P-P_A||；当飞行轨迹存在偏差时，雷达平台实际位置记为

其中[Δx(t) Δy(t) Δz(t)]^T表示载机位置偏差。

任意像素点距雷达平台的实际瞬时距离记为

即：

将公式(1)展开可以得到：

其中ΔR表示距离误差；得到对应相位误差表达式：

式中λ表示雷达发射信号对应波长。

S1、快速因式分解后向投影成像。

根据因式分解的原理确定最佳的初始孔径长度l₀以及每次合并的子图像个数I，即分解因子。

以各个初始子孔径的中心为原点建立极坐标系，划分出子图像的像素点坐标(r,θ)的取值范围,有

其中，R_srn为近端距离，R_srf为远端距离，r表示像素点到子孔径坐标原点的距离，θ表示孔径矢量和距离矢量r之间的夹角，Θ表示积累角；极坐标系下的瞬时距离表示为：

其中，d表示雷达平台到对应坐标系原点的距离。

确定初始子图像的角度向分辨率

同时根据信号带宽确定距离向分辨率ρ_r；方位向采样点个数为N，初始成像子孔径个数为K；通过公式(5)(6)对成像场景中的每一点(r,θ)进行初始成像：

其中，w(θ)为天线方向函数，s(x,τ)为获取的回波，δ(·)为单位冲激函数，x为雷达在方位向的位置，τ为对应的快时间，下标q表示方位向采样点，k表示第k个子孔径，且有q∈[1,N],k∈[1,K]，上标(1)表示第一级(初始)成像。

考虑确定雷达位置对应确定的相位误差，记飞行轨迹偏差导致的方位向相位误差为φ＝{φ₁,φ₂.....φ_N}，代入公式(6)，得到

其中，j表示为虚数单位；则初始成像更新为：

对于i+1级某个成像网格中的某点(r,θ)，计算其在第i级第k幅子图像中的位置(r′,θ′)，将第k幅图像在该点的结果采用插值的方式累加至新的i+1级图像f(r,θ)⁽ⁱ⁺¹⁾中，从图2所示的几何关系中可以得到公式(8)-(10)：

rcos(θ)-r′cosθ′＝d_k (8)

此时，d_k表示为当前子孔径对应的上一级合并时第k个子孔径距离该子孔径几何中心的距离。

根据公式(8)-(10)对子孔径成像结果进行逐级合并，每进行一级合并，将I幅子图像生成次一级子图像，第i+1幅子图像分辨率和第i级子图像分辨率存在以下的关系：

其中，ρ_θm为最终分辨率；得到最终图像：

其中，

表示多次合并插值之后的值，

表示第q个方位向采样点相位误差的估计值，用于补偿相位误差，则f^(F)为最终补偿后的图像。

S2、相位误差估计。

记图像锐度为：

其中i表示第i个像素点，记最终得到的图像f^(F)→f，对应的共轭表示为f^*，则v_i＝f_if_i ^*；建立带有未知相位误差的最大图像锐度的优化模型来求解相位误差

即是对φ的估计，最终使得图像良好聚焦，表示为：

使用坐标下降的迭代方法来求解最优解；在坐标下降的算法中，每一次迭代都使得优化变量按序得以更新，而保持其他变量固定不变；记

为第i次迭代的第n个相位误差校正量，则第i+1次迭代表示为：

最终图像表示如下：

其中，x表示除了第n次采样点的其他所有采样点的后向投影之和，y表示对当前采样点校正后的后向投影；将公式(16)代入(13)得到：

其中，x_i和y_i为第i个像素点在公式(16)下的对应表示；x_i ^*和y_i ^*为x_i和y_i对应的共轭表示，

因此在坐标下降处理下将多维优化求解简化为单变量求极值问题

由公式(13)得到关于φ的导数：

将公式(19)(20)代入公式(18)最终得到：

其中：

利用割线法对单变量优化问题进行求解，具体步骤如下：

T1、选择初始点(φ₀)¹，参数ε＞0,γ＞0，令(φ₁)¹＝(φ₀)¹-γφ′((φ₀)¹)，k:＝0。

T2、若|φ′((φ_k)ⁱ)|＞ε，则利用公式(23)和(24)更新α_k，否则停止运算；其中，上标i表示第i个方位向采样点；

T3、通过公式(25)和(26)计算(φ_k+1)ⁱ,k＝k+1，返回步骤T2。

公式(23)-(26)中的下标k表示对于第i个方位向采样点的第k次迭代求解；

通过上述步骤，即可实现对相位误差进行精确估计并实现良好成像。

下面以一个具体实施例对本发明提供的一种基于最大图像锐度的快速因式分解后向投影SAR自聚焦方法作进一步描述：

S0、系统参数初始化。

本发明实施例采用的SAR几何结构如图3所示，系统仿真参数如下表所示：

仿真参数	取值
		载频	9.6e9Hz
带宽	200e6Hz
		脉冲重复频率	800Hz
时间宽度	5e-6s
		相位误差	Δφ＝9×(t)<sup>2</sup>+8×(t)<sup>3</sup>+7×(t)<sup>4</sup>+6×(t)<sup>5</sup>rad

本发明实施例采用的目标场景如图4所示，图中的黑色圆点为布置于地面上的9个点目标，这9个点分布于x轴和y轴上，沿x方向(切航迹)间隔10m，沿y方向(沿航迹)间隔10m，平台沿y轴运动，且速度v＝30m/s。

波束中心位于场景坐标原点处时记为零时刻，零时刻的雷达平台位置坐标为P_A＝(0 1000 1000)^T，场景中心坐标为(0 0 0)^T；场景中任一点目标的位置坐标为P[x y 0]^T；在MATLAB平台上对成像区域内的任意目标计算其距离历程R，仿真出雷达回波数据，记为s₀(t,τ)。

S1、对回波数据进行快速因式分解后向投影成像。

对回波数据s₀(t,τ)按照下式进行距离压缩：

s′₀(t,τ)＝IFFT(FFT(s₀(t,τ))×FFT(f_τ(τ))) (27)

其中

K_r为距离向调频斜率，取K_r＝4.0e+13，f_c为载频。

在距离压缩后的回波数据中加入相位误差Δφ，误差具体表达式如上表所示，记加入相位误差的回波数据为s(t,τ)。

取分解因子I为2，初始孔径长度为l₀为0.0375m，确定初始子图像的角度向分辨率

为0.0768rad，同时根据信号带宽确定距离向分辨率ρ_r为1.5m。在极坐标系下按照公式(4)计算每一个像素点距雷达平台的双程距离和，按照公式(5)(6)进行初始子图像成像。所有子孔径成像完成后按照公式(8)-(10)对子孔径成像结果进行逐级合并，下一级成像分辨率计算如公式(11)所示，直至得到最终成像结果，带有相位误差的成像结果如图5所示。

S2、相位误差估计。

根据公式(13)计算图像锐度，得到

设置初始误差

以及参数ε＝0.001，γ＝0.001。将

ε，γ代入步骤T1中计算出

将其代入公式(18)计算

根据公式(23)和(24)更新参数(α₀)¹，代入(25)和(26)计算下一个

重复以上步骤计算下一轮迭代的

最终得到所有N个误差估计值

相位误差估计结果如图6所示，误差补偿后成像结果如图7所示。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (3)

1.一种快速因式分解后向投影SAR自聚焦方法，其特征在于，包括以下步骤：

S0、系统参数初始化；

S1、快速因式分解后向投影成像；步骤S1具体为：

根据因式分解的原理确定最佳的初始孔径长度l₀以及每次合并的子图像个数I，即分解因子；

以各个初始子孔径的中心为原点建立极坐标系，划分出子图像的像素点坐标(r,θ)的取值范围,有

其中，R_srn为近端距离，R_srf为远端距离，r表示像素点到子孔径坐标原点的距离，θ表示孔径矢量和距离矢量r之间的夹角，Θ表示积累角；极坐标系下的瞬时距离表示为：

其中，d表示雷达平台到对应坐标系原点的距离；

确定初始子图像的角度向分辨率

同时根据信号带宽确定距离向分辨率ρ_r；方位向采样点个数为N，初始成像子孔径个数为K；通过公式(5)(6)对成像场景中的每一点(r,θ)进行初始成像：

其中，w(θ)为天线方向函数，s(x,τ)为获取的回波，δ(·)为单位冲激函数，x为雷达在方位向的位置，τ为对应的快时间，下标q表示方位向采样点，k表示第k个子孔径，且有q∈[1,N],k∈[1,K]，上标(1)表示第一级成像；

考虑确定雷达位置对应确定的相位误差，记飞行轨迹偏差导致的方位向相位误差为φ＝{φ₁,φ₂.....φ_N}，代入公式(6)，得到

其中，j表示为虚数单位；则初始成像更新为：

对于i+1级某个成像网格中的某点(r,θ)，计算其在第i级第k幅子图像中的位置(r′,θ′)，将第k幅图像在该点的结果采用插值的方式累加至新的i+1级图像f(r,θ)⁽ⁱ⁺¹⁾中，根据公式(8)-(10)对子孔径成像结果进行逐级合并：

rcos(θ)-r′cosθ′＝d_k (8)

此时，d_k表示为当前子孔径对应的上一级合并时第k个子孔径距离该子孔径几何中心的距离；

每进行一级合并，将I幅子图像生成次一级子图像，第i+1幅子图像分辨率和第i级子图像分辨率存在以下的关系：

其中，ρ_θm为最终分辨率；得到最终图像：

其中，

表示多次合并插值之后的值；

表示第q个方位向采样点相位误差的估计值，用于补偿相位误差，则f^(F)为最终补偿后的图像；

S2、相位误差估计。

2.根据权利要求1所述的快速因式分解后向投影SAR自聚焦方法，其特征在于，所述步骤S0具体为：

以场景参考目标点O＝[0 0 0]^T为原点建立坐标系，场景中任意像素点位置记为P＝[xy 0]^T，雷达平台理想位置记为P_A＝[0 vt h]^T；

其中，v为载机沿着y轴的理想飞行速度，h为载机沿着y轴的理想飞行高度，t为方位慢时间；

任意像素点距雷达平台的理想瞬时距离记为R(t)＝||P-P_A||；当飞行轨迹存在偏差时，雷达平台实际位置记为

其中[Δx(t) Δy(t) Δz(t)]^T表示载机位置偏差；

任意像素点距雷达平台的实际瞬时距离记为

即：

将公式(1)展开可以得到：

其中ΔR表示距离误差；得到对应相位误差表达式：

式中λ表示雷达发射信号对应波长。

3.根据权利要求1所述的快速因式分解后向投影SAR自聚焦方法，其特征在于，所述步骤S2具体为：

记图像锐度为：

其中i表示第i个像素点，记最终得到的图像f^(F)→f，对应的共轭表示为f^*，则v_i＝f_if_i ^*；建立带有未知相位误差的最大图像锐度的优化模型来求解相位误差

即是对φ的估计，最终使得图像良好聚焦，表示为：

使用坐标下降的迭代方法来求解最优解；在坐标下降的算法中，每一次迭代都使得优化变量按序得以更新，而保持其他变量固定不变；记

为第i次迭代的第n个相位误差校正量，则第i+1次迭代表示为：

最终图像表示如下：

其中，x表示除了第n次采样点的其他所有采样点的后向投影之和，y表示对当前采样点校正后的后向投影；将公式(16)代入(13)得到：

其中，x_i和y_i为第i个像素点在公式(16)下的对应表示；x_i ^*和y_i ^*为x_i和y_i对应的共轭表示，

因此在坐标下降处理下将多维优化求解简化为单变量求极值问题

由公式(13)得到关于φ的导数：

将公式(19)(20)代入公式(18)最终得到：

其中：

利用割线法对单变量优化问题进行求解，具体步骤如下：

T1、选择初始点(φ₀)¹，参数ε＞0,γ＞0，令(φ₁)¹＝(φ₀)¹-γφ′((φ₀)¹)，k:＝0；

T2、若|φ′((φ_k)ⁱ)|＞ε，则利用公式(23)和(24)更新α_k，否则停止运算；其中，上标

i表示第i个方位向采样点；

T3、通过公式(25)和(26)计算(φ_k+1)ⁱ,k＝k+1，返回步骤T2；

公式(23)-(26)中的下标k表示对于第i个方位向采样点的第k次迭代求解；

通过上述步骤，实现对相位误差进行精确估计并实现良好成像。

Images