CN109655830A - 一种针对星上去走动与二维空变校正的地球同步轨道sar斜视成像处理装置 - Google Patents

Abstract

本发明公开了针对星上去走动与二维空变校正的地球同步轨道SAR斜视成像处理装置，本装置顺次采用星上卫星调整接收窗去除信号走动、去走动回波链路下传、星地几何关系构建空变模型、距离向残余走动校正、距离徙动校正与压缩、方位向部分去斜降低时间跨度、方位向残余走动去除、方位向非线性变换与压缩、最后经几何校正与成像结果输出。与传统算法比较，本发明利用GEO SAR下重频极低的特点，通过星上去走动降低距离和方位向耦合的同时提高了星上数据存储的有效性；对于接收到的回波通过部分对方位和距离向进行精细修正，提高了算法的聚焦精度。

Description

一种针对星上去走动与二维空变校正的地球同步轨道SAR斜 视成像处理装置

技术领域

本发明涉及GEO SAR的斜视成像。具体来说是利用GEO SAR重频极低的特点进行星上去走动，并采用基于去走动后的二维空变模型进行二维空变性校正，从而实现GEO SAR斜视处理的方法，属于遥感成像处理领域。

背景技术

合成孔径雷达(Synthetic Aperture Radar，SAR)是一种全天时、全天候的主动式信息获取系统，合成孔径雷达的全链路系统可以分为三个部分，即雷达载体平台上的回波数据获取部分、空间段数据传输部分以及地面段数据处理部分。合成孔径雷达系统通过发射和接收脉冲信号，获取合成孔径雷达原始回波数据。

地球同步轨道合成孔径雷达(GEO SAR，Geosynchronous orbit SyntheticAperture Radar)卫星运行在高度约35,786千米的倾斜或非圆的地球静止轨道上的SAR卫星。它的轨道周期与地球自转周期相同，使得卫星相对于地球的理论运动轨迹线是闭合的，这意味着每隔24小时卫星都会回归到同一位置上空。这种“经度锁定”特性确保了GEO SAR能够对感兴趣区域进行长时间观测。闭合轨迹的这种形式当辅以斜视模式时即可构成GEOSAR所特有的“凝视”工作模式：即通过对天线波束的间断性调整实现对区域进行全轨持续观测。这种特性在地质灾害，如洪涝，地震等区域进行持续监测方面具有重要的意义。

由于卫星轨道极高，卫星的运行速度相较于低轨SAR下降到和地球自转速度相同的量级，这将产生两个主要的影响：一是累积与低轨相同分辨率的孔径时间变得极长，使得弯曲轨道特性不可忽略，传统的双曲模型不再适用；二是地球自转效应必须加以考虑，地面不可再认为是平面。综合这两方面，传统的平地直线模型不再成立，而高效率低轨成像处理所依赖的同距离单元多普勒历程相同的假设不再成立，因此回波信号出现明显的二维空变性。回波信号的二维空变性是影响GEO SAR高精度成像的根本原因，因此GEO SAR成像处理的研究方向都集中在空变性校正的方面。

常规情况下GEO SAR的数据量已然很大，当卫星处于斜视观测时数据量就会急剧上升，然而观测区域场景大小并没改变，这对于星上存储是极大地浪费。同时，斜视观测会加重距离和方位向的耦合，对后续的处理增加了难度。时域成像算法虽然简便易行，但是信号处理的效率难以保证。因而GEO SAR成像的研究方向还主要集中在频域成像处理上。在GEO SAR信号频域处理方面，可以将现有研究成果大致分为基于RD算法的成像算法、基于Chirp变标的成像算法以及基于波数域算法的成像算法。在以RD算法为基础的算法方面，胡程在2014年发表了《An Improved Frequency Domain Focusing Method inGeosynchronous SAR》，文中称所提方法可以实现40公里×40公里、10米分辨率成像。除了此文所提算法之外，其余算法都是基于双曲模型进行算法研究，主要解决方案是分块处理，理论贡献较弱。在以波数域算法为基础的算法方面，2015年，胡斌利用多项式斜距模型和级数反演法得出的二维频谱表达式，提出一种基于多项式斜距模型的STOLT变标方法，在线性空变性的假设下修正了低轨波数域算法，具体请参见《Generalized Omega-K Algorithmfor Geosynchronous SAR Image Formation》。由于基于Chirp变标的算法不需要类如RD和Omega-K一样的插值操作，因而计算效率高，这使得目前GEO SAR成像算法的研究大多集中在以Chirp变标为基础的方向上。GEO SAR回波信号的二维空变性使得高性能的成像算法必须要考虑二维空变的校正。最早引入二维空变性校正的是孙光才(《A 2-D Space-VariantChirp Scaling Algorithm Based on the RCM Equalization and Subband Synthesisto Process Geosynchronous SAR Data》)，但他提出的这种方法需要方位向子带合成，会产生寄生旁瓣。而后有代表性的则是《Modeling and Processing of Two-DimensionalSpatial-Variant Geosynchronous SAR Data》，这里融合了Omega-K算法与CS算法，具有比较好的成像性能。但是这些算法均是针对正侧视模式进行设计的。

发明内容

本发明针对GEO SAR斜视观测情况下距离走动量大和二维空变性严重等问题，利用同步轨道SAR重频极低的特点，提出回波信号星上去走动接收策略。针对接收的回波存在严重的二维空变的特点，建立相应的二维空变信号模型。对于距离向处理，先降低去除距离向残余走动，然后采用基于场景二维空变模型的非线性变标进行徙动和距离脉冲压缩；对于方位向，则采用部分去斜降低信号的时域跨度，然后进行方位向残余走动去除，再通过时频域联合的非线性变换进行方位向均衡处理，从而实现方位向脉冲压缩，最后进过几何校正形成最终图像。本发明采用星上去走动方式降低回波信号数据量，提高星上数据存储的同时简化后续成像处理压力。采用的部分去斜技术和时频联合非线性变标技术提高信号一次和二次空变性的校正能力，图像聚焦效果良好。

本发明利用GEO SAR重频极低的特点，通过星上去走动使得回波信号的数据量变小，降低后续的处理压力。同时，引入去走动后的二维空变模型，利用时频混合变标和部分去斜技术去除方位向和距离向的一次和二次空变，从而实现斜视观测的二维精确聚焦。

本发明的一种针对星上去走动与二维空变校正的GEO SAR斜视成像处理装置，该装置包括有雷达参数与空变性获取模块、回波去走动模块、方位时域变标模块、徙动校正与脉冲压缩模块和方位向处理模块；

回波去走动模块，该模块在星上完成，通过对接收回波按照后续给出的接收窗调整策略去场景中心走动，然后频谱搬移将方位向信号从射频带搬移到基带，获得去走动后的回波信号；将该回波作为后续处理操作的输入；

方位时域变标模块，该模块对上述回波进行方位时间变标，用来校正回波信号的方位向一次空变性；该模块接收雷达参数输出模块的全部输出参数，并将方位时间变标后的信号作为该模块的输出；

徙动校正与脉冲压缩模块，该模块接收雷达参数输出模块以及方位时间变标模块的全部输出参数，进行距离单元徙动校正及距离脉冲压缩，并将距离单元徙动校正及距离压缩后的回波信号作为输出参数；

方位向处理模块，该模块利用部分去斜处理技术将方位向的时域跨度缩短，然后方位向划分子孔径去除多普勒中心频率，之后利用时频混合变标去除信号的方位一次和二次空变性，以实现方位向信号聚焦；最后对信号进行几何校正，输出最终图像。

本发明GEO SAR斜视成像处理装置的优点在于：

①本发明GEO SAR斜视成像处理装置利用星上回波去走动，降低了信号的回波复杂度和数据量，降低后续成像过程的计算负担，可以用于对目标多角度持续观测。

②本发明GEO SAR斜视成像处理装置的斜距模型及信号频谱模型都采用精度极高的高阶多项式模型，保证了后续成像处理精度。

③本发明GEO SAR斜视成像处理装置利用时频混合变标的方式校正了回波信号的二维高阶空变性，提高了空变性校正的能力及精度。

④本发明GEO SAR斜视成像处理装置在方位向分块去除多普勒中心频率，便于并行计算。

⑤本发明GEO SAR斜视成像处理装置在信号处理当中考虑了信号频谱幅度非矩形特性并对此特性进行了补偿，二维剖面聚焦效果更好。

⑥本发明GEO SAR斜视成像处理装置在频域实现成像处理效果高。

附图说明

图1是本发明针对星上去走动与二维空变校正的GEO SAR斜视成像处理装置的结构框图。

图2是本发明针对星上去走动与二维空变校正的GEO SAR斜视成像处理装置进行星上与地面划分的结构框图。

图3是传统开窗时间与本发明开窗时间对比图。

图4是本发明针对星上去走动与二维空变校正的GEO SAR斜视成像处理装置流程图。

图5是仿真目标点布局示意图。

图5A1是目标T1的点扩展函数轮廓图。 图5A2是目标T1的方位向剖面图。

图5A3是目标T1的距离向剖面图。 图5B1是目标T2的点扩展函数轮廓图。

图5B2是目标T2的方位向剖面图。 图5B3是目标T2的距离向剖面图。

图5C1是目标T3的点扩展函数轮廓图。 图5C2是目标T3的方位向剖面图。

图5C3是目标T3的距离向剖面图。 图5D1是目标T4的点扩展函数轮廓图。

图5D2是目标T4的方位向剖面图。 图5D3是目标T4的距离向剖面图。

图5E1是目标T5的点扩展函数轮廓图。 图5E2是目标T5的方位向剖面图。

图5E3是目标T5的距离向剖面图。 图5F1是目标T6的点扩展函数轮廓图。

图5F2是目标T6的方位向剖面图。 图5F3是目标T6的距离向剖面图。

图5G1是目标T7的点扩展函数轮廓图。 图5G2是目标T7的方位向剖面图。

图5G3是目标T7的距离向剖面图。 图5H1是目标T8的点扩展函数轮廓图。

图5H2是目标T8的方位向剖面图。 图5H3是目标T8的距离向剖面图。

图5I1是目标T9的点扩展函数轮廓图。 图5I2是目标T9的方位向剖面图。

图5I3是目标T9的距离向剖面图。

具体实施方式

下面将结合附图和实施例对本发明做进一步的详细说明。

本发明的一种针对星上去走动与二维空变校正的GEO SAR斜视成像处理装置，本装置对信号的处理是顺次采用星上卫星调整接收窗去除信号走动→去走动回波链路下传→星地几何关系构建空变模型→距离向残余走动校正→距离徙动校正与压缩→方位向部分去斜降低时间跨度→方位向残余走动去除→方位向非线性变换与压缩→几何校正与成像结果输出。与传统算法比较，本发明利用GEO SAR下重频极低的特点，通过星上去走动降低距离和方位向耦合的同时提高了星上数据存储的有效性；对于接收到的回波通过部分对方位和距离向进行精细修正，提高了算法的聚焦精度。

(一)雷达参数与空变性获取模块

在本发明中，雷达参数与空变性获取模块用于为整个成像处理系统输入参数，该模块的输出参数分为两类，一类是雷达工作参数f₁＝{Na,Nr,PRF,Fs,Br,Tp,Kr,λ,T_{s_general},t_center,ta,fa,tr,fr}，一类是回波各像素点对应地面目标点的斜距系数的空变参数f₂＝{P_n,Q_n,0≤n≤5}。

Na是发射脉冲总数。Nr是距离采样点总数。PRF是脉冲重复频率。Fs是距离向采样率。Br是雷达工作带宽。Tp是发射脉冲宽度。Kr是距离向调频率。λ是工作波长。T_{s_general}为采样延时，也称为开窗时间。t_center为仿真中心时刻。

ta为方位向时间序列，即ta＝[η₁,η₂,...,η_m,...,η_Na]^T，上角标T为转置运算符。η₁为雷达发射的第1个脉冲时刻，η₂为雷达发射的第2个脉冲时刻，η_m为雷达发射的第m个脉冲时刻，且η_Na为雷达发射的最后一个脉冲时刻。m也表示雷达发射的脉冲时刻标识号，m＝1,2,…,Na。

fa为方位向频率序列，即fa＝[f_η,1,f_η,2,...,f_η,m′,...,f_η,Na]^T，f_η,1表示方位频域的第1个频点，f_η,2表示方位频域的第2个频点，f_η,m′表示方位频域的第m′个频点，且f_η,Na表示方位频域的最后一个频点。m′也表示方位频域的频点标识号，m′＝1,2,…,Na。

tr为距离向时间序列，且tr＝[τ₁,τ₂,…,τ_k,…,τ_Nr]^T，τ₁表示从接收窗开启时间起的第1个距离向采样时刻，τ₂表示从接收窗开启时间起的第2个距离向采样时刻，τ_k表示从接收窗开启时间起的第k个距离向采样时刻，且τ_Nr表示从接收窗开启时间起的最后一个距离向采样时刻。k也表示从接收窗开启时间起的距离向采样时刻标识号，k＝1,2,...,Nr。

fr为距离向频率序列，且fr＝[f_τ,1,f_τ,2,...,f_τ,k′,...,f_τ,Nr]^T，f_τ,1表示距离向频域的第1个频点，f_τ,2表示距离向频域的第2个频点，f_τ,k′表示距离向频域的第k′个频点，且f_τ,Nr表示距离向频域的最后一个频点。k′也表示距离向频域频点标识号，k′＝1,2,...,Nr。

P_n是第n阶斜距系数的距离向空变参数矩阵，尺寸为1×5。

Q_n是第n阶斜距系数的方位空变参数矩阵，尺寸为Nr×5。

在GEO SAR中，根据星地几何关系计算出的回波像素点地面对应点的第n阶斜距系数分布是个二维矩阵，记为其中r_n,m,k表示穿越时刻为中心斜距为的目标点的第n阶斜距系数，也就是斜距系数矩阵当中第m行、第k列的元素。

在本发明中，所述也是去除多普勒中心频率后的斜距系数。另外：

也可以将所述表示为元素行向量的列向量

也可以将所述表示为元素列向量的行向量

其中表示的第m行行－列－向量；

其中表示的第k列列－行－向量。

在本发明中，斜距系数的空变性建模方法如下：

步骤1-1，从中提取出第Na/2+1行的行－列－向量利用MATLAB当中的Polyfit函数，将所述对距离时间tr-2R_ref/c进行多项式拟合，得到关于tr-2R_ref/c的各阶多项式系数P_n＝[p_n,4p_n,3p_n,2p_n,1p_n,0]，即是第n阶斜距系数的距离空变系数；

步骤1-2，依次提取出的每一列的列－行－向量利用MATLAB当中的Polyfit函数，将所述对方位时间ta进行多项式拟合，得到关于ta的各阶多项式系数[q_n,4,k q_n,3,k q_n,2,k q_n,1,k q_n,0,k]，即是第n阶第k位置斜距系数的方位空变系数；

步骤1-3，将[q_n,4,k q_n,3,k q_n,2,k q_n,1,k q_n,0,k]存入尺寸为Nr×5的空变系数矩阵Q_n的第k行当中；

步骤1-4，从1至Nr遍历第k列并重复执行步骤1-1至步骤1-3，即可得到Q_n的所有元素；

步骤1-5，使n从1至5遍历，并重复执行步骤1-4，即可得到各阶斜距系数的距离空变系数矩阵{P_n,0≤n≤5}以及方位空变系数{Q_n,0≤n≤5}。

为斜距系数中第1个行向量的列向量。为斜距系数中第m个行向量的列向量，简称为行－列－向量。为斜距系数中第Na个行向量的列向量。为斜距系数中第1个列向量的行向量。为斜距系数中第k个列向量的行向量，简称为列－行－向量。为斜距系数中第Nr个列向量的行向量。r_n,m,1为第n阶斜距系数中第m行第1列的元素。r_n,m,k为第n阶斜距系数中第m行第k列的元素。为第n阶斜距系数中第m行第Nr列的元素。r_n,1,k为第n阶斜距系数中第1行第k列的元素。为第n阶斜距系数中第Na行第k列的元素。p_n，4表示将拟合成tr-2R_ref/c的多项式的第4阶系数。p_n,3表示将拟合成tr-2R_ref/c的多项式的第3阶系数。p_n,2表示将拟合成tr-2R_ref/c的多项式的第2阶系数。p_n,1表示将拟合成tr-2R_ref/c的多项式的第1阶系数。p_n,0表示将拟合成tr-2R_ref/c的多项式的第0阶系数。R_ref表示参考中心的斜距。q_n,4,k表示将拟合成ta的多项式的第4阶系数。q_n,3,k表示将拟合成ta的多项式的第3阶系数。q_n,2,k表示将拟合成ta的多项式的第2阶系数。q_n,1,k表示将拟合成ta的多项式的第1阶系数。q_n,0,k表示将拟合成ta的多项式的第0阶系数。

(二)回波去走动模块

参见图1、图2、图3、图4所示，回波去走动模块在星上完成，通过对接收回波按照“接收窗调整策略”去场景中心走动，然后采用频谱搬移将方位向信号从射频带搬移到基带，获得去走动后的回波信号。

相比较于低轨SAR卫星，GEO SAR卫星的重频极低，因此可以通过调整接收窗时间的方式进行星上去走动。传统SAR的开窗时间T_{s_general}是固定不变的(如图3中虚线所示)，此时信号的走动量极大，造成观测同样大小场景产生的无用数据非常多。当动态地改变开窗时间(如图3中实线所示)后，此时全场景的绝大部分走动量被去除，回波信号则转化为类似正侧视模式。因此接收窗调整策略应满足公式(1)。

按照公式(1)处理后，接收回波信号变成类正侧视模式，并将数据通过链路下传到地面接收站，记此时的回波信号为S₁。在图3中，按照实线动态地改变接收窗位置，可以将原本的斜视回波变成图中所示的类正侧视回波形式。

T_{s_RWC}为新的变化的开窗时间。R_{ref_RWC}(t_a)为去走动后的参考斜距。R_ref(t_a)为参考中心的斜距历程。c为光速。为t_a的n次幂。n为斜距拟合次数。n！表示n的阶乘。为r_{n_ref}的n次幂。r_{0_ref}为参考目标点的0阶斜距系数。r_{1_ref}为参考目标点的一阶斜距系数。r_{n_ref}(n＝0～5)为参考目标点的斜距系数，在“雷达参数与空变性获取模块”已给出详细的获取方法。

(三)方位时域变标模块

参见图1、图2、图4所示，方位时域变标模块对所述回波S₁进行方位时间变标，用来校正回波信号的方位向一次空变性。在本发明中的方位时域变标的方法如下：

步骤3-1，将回波信号S₁做距离向傅里叶变换；

在本发明中，提取S₁的第一行，将此行矩阵的前个元素与后个元素对换位置，得到交换位置后的行向量g_交换；然后对所述g_交换进行快速傅里叶变换FFT，得到频谱G_交换；最后对所述G_交换进行前个元素与后个元素对换位置，即完成对回波信号S₁第一行数据的距离向傅里叶变换；

对所述S₁的第二行至Na行采用上述相同方法处理，从而完成对回波信号S₁进行距离向傅里叶变换的操作，记为距离向傅里叶变换后的信号为S_{rng_fft}。

由于将回波信号S₁做距离向的傅里叶变换与快速傅里叶变换(FFT)有所不同，快速傅里叶变换(FFT)的方法参见郑君里等编著的《信号与系统》下册9.6节，2000年5月第2版。

步骤3-2，构造方位时间变标滤波器：

其中H_AS1表示方位时间变标滤波器，j表示虚部，f₀表示载波频率，r_pt(ta)表示方位依赖的调制项，且A_n表示驻定项相位点系数，所述A_n为公式(3)。

步骤3-3，将距离傅里叶变换后的回波信号S_{rng_fft}与H_AS1逐点相乘，即可完成方位时间变标，将变标后的信号记为S₂作为此模块的输出之一。

在本发明中，计算方位时间变标后的斜距系数空变参数，计算方法如下：

步骤(31)方位时间变标后的距离空变系数P′_n＝P_n+[0 0 0 0 A_n]，0≤n≤5；

步骤(32)方位时间变标后的方位空变系数计算方法为：

步骤(33)将{P′_n,0≤n≤5}、{Q′_n,0≤n≤5}以及参数{A_n,1≤n≤7}作为方位时域变标模块的第二组输出参数，记为f₃＝{P′_n,Q′_n,A_n}。

本发明的方位时域变标模块采用了五阶方位依赖调制项做方位时域均衡，从而使得的徙动校正和方位聚焦范围扩大，结果精确。

Q₀表示第0阶斜距系数的方位空变参数。Q′₀表示时域变标后的第0阶斜距系数的方位空变参数。Q₁表示第1阶斜距系数的方位空变参数。Q′₁表示时域变标后的第1阶斜距系数的方位空变参数。Q₂表示第2阶斜距系数的方位空变参数。Q′₂表示时域变标后的第2阶斜距系数的方位空变参数。Q₃表示第3阶斜距系数的方位空变参数。Q′₃表示时域变标后的第3阶斜距系数的方位空变参数。Q₄表示第4阶斜距系数的方位空变参数。Q′₄表示时域变标后的第4阶斜距系数的方位空变参数。Q₅表示第5阶斜距系数的方位空变参数。Q′₅表示时域变标后的第5阶斜距系数的方位空变参数。A₁表示第1阶驻定项相位点系数。A₂表示第2阶驻定项相位点系数。A₃表示第3阶驻定项相位点系数。A₄表示第4阶驻定项相位点系数。A₅表示第5阶驻定项相位点系数。A₆表示第6阶驻定项相位点系数。A₇表示第7阶驻定项相位点系数。

(四)徙动校正与脉冲压缩模块

参见图1、图2、图4所示，徙动校正与脉冲压缩模块接收方位时间变标后的回波信号S₂以及斜距系数空变参数f₃＝{P′_n,Q′_n,A_n}，完成对S₂进行徙动校正及距离压缩的操作，具体操作方法如下：

步骤4-1，计算任意一个目标点的频谱相位；

假设任意一个目标点其穿越时刻及穿越时刻的中心距离分别为η_m,cross＝ta(m)、中心斜距为则这个目标点的各阶斜距系数表达式如下：

在本发明中，依据公式(5)其频谱相位表达式为：

θ_2df,m,k表示二维频谱相位。

Φ_u,m,k表示二维频谱距离频率展开式中第u阶系数，u表示二维频谱距离频率的标识号，且u≥3。u≥3，表示在l+u-2个数中选取u个数的个数，即组合数操作，表示频谱相位中第l阶斜距系数中第m行第k列的元素，表示参考中心频率f₀的u-1次幂。在本发明中，l≥2，表示在v+1个数中选取u个数的个数。

当u＝0时，r_0,m,k为第0阶斜距系数中第m行第k列的元素，表示频谱相位中第0阶斜距系数中第m行第k列的元素，λ是工作波长，B_l-1,m,k表示将二维频谱距离频率展开式中的第m行第k列的第l－1项系数展开成方位向频率fa多项式的系数，l表示属于fa的序号。在本发明中，v是求和标识号，v＝1,2,…,9，B_v,m,k表示将二维频谱距离频率展开式中的第m行第k列的第v项系数展开成方位向频率fa多项式的系数，r_1,m,k表示由公式(5)计算得到的第1阶斜距系数。

当u＝1时，

当u＝2时， 表示在l个数中选取u个数的个数，即组合数操作，Kr是距离向调频率。

在本发明中，其中B_v,m,k，1≤v≤9表示目标点的驻定相位点表达式，每一个B_v,m,k的表达式都是各阶斜距系数的四则运算组合，所述四则运算组合参见《CorrectingSpatial Variance of RCM for GEOSAR Imaging based on Time-Frequency Scaling》，2016年7月。

步骤4-2，计算频谱相位系数的空变参数；

将目标点的各阶斜距系数表达式(5)代入到B_v,m,k的表达式当中并整理出关于穿越时刻η_m,cross以及穿越时刻中心斜距ΔR_k的多项式：

再将斜距一次项系数即公式(7)的B_v,m,k代入到P_l,m,k的表达式中并整理出关于穿越时刻η_m,cross以及穿越时刻中心斜距ΔR_k的多项式：

再将公式(8)的P_l,m,k代入到Φ_u,m,k中，并整理出关于穿越时刻以及穿越时刻中心斜距的多项式：

步骤4-3，补偿高次相位；

首先将方位时间变标后的信号S₂进行方位向傅里叶变换，方位向傅里叶变换只需将距离向傅里叶变换当中所有对行矩阵的操作变成对列矩阵的操作即可。进行方位向傅里叶变换后，即获得方位时间变标后的二维频域回波信号S_2,2df。构造高阶相位补偿滤波器：

在本发明中，当中所有的矩阵乘法都是逐元素点乘。将S_2,2df与H_{4plus_Y}逐点相乘即可得到补偿高次相位后的二维频谱，记为S_4plus,2df。

步骤4-4，非线性Chirp变标；

将补偿高次相位后的二维频谱S_4plus,2df距离向傅里叶反变换，即可得到补偿高次相位后的RD域频谱，记为S_4plus,rd，距离向傅里叶反变换的方法与距离向傅里叶变换的步骤相同，只是将过程当中所用的FFT函数替换成IFFT函数即可。

构造非线性变标方程，表达式如下：

QC₂表示距离多普勒域下的二阶变标系数，且QC₂＝αK_mref。

QC₃表示距离多普勒域下的三阶变标系数，且

对公式(11)所有的矩阵乘除法为逐点运算，将S_4plus,rd与H_NCS进行点乘操作，并对点乘操作后的二维矩阵记为S_ncs,rd。

步骤4-5，距离压缩及徙动校正；

将非线性变标后的距离多普勒域信号S_ncs,rd做距离向傅里叶变换，变换至二维频域，操作方法参见前文，记距离向傅里叶变换后的信号为S_ncs,2df。构造距离压缩及徙动校正滤波器如下：

表示行矩阵Φ_1,ref的第个元素。将S_ncs,2df与H_{rpc_rcmc}逐点相乘，并将相乘后的二维矩阵进行距离向傅里叶反变换即完成徙动校正及脉冲压缩，记变换回距离多普勒域的信号记为S₃，S₃即是此模块的输出参量。

p′_n,0表示p_n,0经方位时间变标后的距离向第0阶空变系数。p′_n,1表示p_n,1经方位时间变标后的距离向第1阶空变系数。p′_n,2表示p_n,2经方位时间变标后的距离向第2阶空变系数。p′_n,3表示p_n,3经方位时间变标后的距离向第3阶空变系数。p′_n,4表示p_n,4经方位时间变标后的距离向第4阶空变系数。q′_n,1,k表示q_n,1,k经方位时间变标后的方位向第1阶空变系数。q′_n,2,k表示q_n,2,k经方位时间变标后的方位向第2阶空变系数。q′_n,3,k表示q_n,3,k经方位时间变标后的方位向第3阶空变系数。q′_n,4,k表示q_n,4,k经方位时间变标后的方位向第4阶空变系数。CO_v,0表示第v次沿距离向拟合的第0次系数。CO_v,1表示第v次沿距离向拟合的第1次系数。CO_v,2表示第v次沿距离向拟合的第2次系数。CO_v,3表示第v次沿距离向拟合的第3次系数。CO_v,4表示第v次沿距离向拟合的第4次系数。DO_v,1,k表示第v次沿方位向拟合的第1行第k列系数。DO_v,2,k表示第v次沿方位向拟合的第2行第k列系数。DO_v,3,k表示第v次沿方位向拟合的第3行第k列系数。DO_v,4,k表示第v次沿方位向拟合的第4行第k列系数。CO_l,0表示第l次沿距离向拟合的第0次系数。CO_l,1表示第l次沿距离向拟合的第1次系数。CO_l,2表示第l次沿距离向拟合的第2次系数。CO_l,3表示第l次沿距离向拟合的第3次系数。CO_l,4表示第l次沿距离向拟合的第4次系数。DO_l,1,k表示第l次沿方位向拟合的第1行第k列系数。DO_l,2,k表示第l次沿方位向拟合的第2行第k列系数。DO_l,3,k表示第l次沿方位向拟合的第3行第k列系数。DO_l,4,k表示第l次沿方位向拟合的第4行第k列系数。Φ_0,ref表示参考中心位置第0阶相位展开系数。Φ_1,ref表示参考中心位置第1阶相位展开系数。Φ_2,ref表示参考中心位置第2阶相位展开系数。Φ_u,ref表示参考中心位置第u阶相位展开系数。M_1,r表示目标点第0阶相位展开系数的第1阶距离依赖系数。r表示距离向标识号。M_2,r表示目标点第0阶相位展开系数的第2阶距离依赖系数。M_3,r表示目标点第0阶相位展开系数的第3阶距离依赖系数。M_4,r表示目标点第0阶相位展开系数的第4阶距离依赖系数。M_1,a,k表示目标点第0阶相位展开系数的第1阶方位向第k列依赖系数。a表示方位向标识号。M_2,a,k表示目标点第0阶相位展开系数的第2阶方位向第k列依赖系数。M_3,a,k表示目标点第0阶相位展开系数的第3阶方位向第k列依赖系数。M_4,a,k表示目标点第0阶相位展开系数的第4阶方位向第k列依赖系数。L_1,r表示目标点第1阶相位展开系数的第1阶距离依赖系数。L_2,r表示目标点第1阶相位展开系数的第2阶距离依赖系数。L_3,r表示目标点第1阶相位展开系数的第3阶距离依赖系数。L_4,r表示目标点第1阶相位展开系数的第4阶距离依赖系数。L_1,a,k表示目标点第1阶相位展开系数的第1阶方位向第k列依赖系数。L_2,a,k表示目标点第1阶相位展开系数的第2阶方位向第k列依赖系数。L_3,a,k表示目标点第1阶相位展开系数的第3阶方位向第k列依赖系数。L_4,a,k表示目标点第1阶相位展开系数的第4阶方位向第k列依赖系数。J_1,r表示目标点第2阶相位展开系数的第1阶距离依赖系数。J_2,r表示目标点第2阶相位展开系数的第2阶距离依赖系数。J_3,r表示目标点第2阶相位展开系数的第3阶距离依赖系数。J_4,r表示目标点第2阶相位展开系数的第4阶距离依赖系数。J_1,a,k表示目标点第2阶相位展开系数的第1阶方位向第k列依赖系数。J_2,a,k表示目标点第2阶相位展开系数的第2阶方位向第k列依赖系数。J_3,a,k表示目标点第2阶相位展开系数的第3阶方位向第k列依赖系数。J_4,a,k表示目标点第2阶相位展开系数的第4阶方位向第k列依赖系数。Y表示距离向变标扰动系数，且K_mref表示参考位置距离多普勒域的调制系数，表示K_mref的二次幂，表示K_mref的三次幂，且α表示一阶变标系数，且L_1,r,ref表示L_1,r在参考频率处的取值；β表示二阶变标系数，且L_2,r,ref表示L_2,r在参考频率处的取值；表示将目标点距离多普勒域的调制系数展开成距离向多项式的一阶系数，且Φ_3,ref表示参考中心位置第3阶相位展开系数。

(五)方位向处理模块

参见图1、图2、图4所示，方位向处理模块接收徙动校正及距离压缩模块的输出参量S₃、雷达参数输出模块输出参量，先对信号采用部分去斜去除方位向依赖的多普勒中心频率，逐列进行时频混合变标及方位压缩操作，完成方位聚焦并形成最终图像。方位向处理模块的操作为：

步骤5-1，构造滤波器恢复原始信号的斜距模型；

其中r_pt(ta)在方位时间变标模块已有定义。将输入参量S₃进行方位向傅里叶反变换，记为S_rngproc,td，将H_-AS1与S_rngproc,td逐点相乘，记相乘后的二维矩阵进行方位向傅里叶变换，记为S_acs2。然后，计算r_pt(ta)的最大值以及最小值，分别记为r_pt,max与r_pt,min，令N_pt为最小的大于r_pt,max-r_pt,min的2的整数次幂。另外，构造一个Na×Nr的空矩阵，记为S_image。

依次对位于S_acs2每一列的目标点进行方位聚焦。下面以S_acs2的第k列为例给出此模块的操作方法，k可以表示至范围内的任意一个数。

在本发明中，提取出S_acs2的第列至第列，将此列矩阵记为S_acs2,k，尺寸为N_a×N_pt，与此同时提取出原始各阶斜距系数的方位空变系数Q_n的第k列以及P_n，记Q_n的第k列为Q_n,k，其中Q_n,k以及P_n均为尺寸为1×5。

步骤5-2，构造线性校正函数；

q_1,1,k表示将1阶斜距系数拟合成ta的多项式的第k列的1阶系数。

将该函数与S_acs2,k相乘，相乘之后修改空变模型系数为公式(15)，其余各个空变系数保持不变，记此时的信号为S_acs3,k；

步骤5-3，构造时频混合变标—频域变标滤波器；

Y_3,k表示方位向时域变标三次调制系数，且

Y_4,k表示方位向时域变标四次调制系数，且

将S_acs3,k变换到方位频域，并与H_fs,k逐点相乘，完成方位频率变标。将频率变标后的矩阵进行方位傅里叶反变换，记此时的二维信号矩阵为S_acs3td,k。

步骤5-4，计算频率变标后的斜距系数空变系数；

经过方位频率变标后，该目标点的各阶斜距系数变为如下形式：

将Y_3,k以及Y_4,k代入到上式当中，即可得到r″_n,m,k关于穿越时刻η_m,cross以及ΔR_k的多项式：

步骤5-5，构造时频混合变标—时域变标滤波器；

构造方位时域变标滤波器：

其中r′_pt(ta)表示方位频率变标后的方位依赖的调制项，且A′表示方位频率变标后的驻定相位点系数，并且

将S_acs3td,k与H_AS2逐点相乘，完成方位频率变标后的方位时间变标，记时间变标后波信号为S′_acs3td，k，另外，方位时间变标后的斜距系数计算公式如下：

步骤5-6，频域部分去斜；

将上述变标后的信号S′_acs3td,k变换到二维频域，从而构建部分去斜和高阶相位补偿滤波器：

其中

各个的计算公式如下：

表示频域部分去斜后第v阶驻定相位点系数。

表示从v+1中选取i个的个数，即组合数操作。

在本发明中，将经过该操作后的信号变换到时域，并记该信号为S_pda,k，以作后部分输入。

步骤5-7，分子场景去除方位向一次空变性；

计算孔径时域跨度缩短后的时长为(1-χ)T_s，T_s为观测孔径时间；依据此获得方位向可分割最小跨度，即对每一块方位数据做去除方位向一次空变性处理，从而构建所述每一块去除方位一次空变性滤波器；

其中为当前子块的方位向时刻，即波束穿越时刻。将S_pda,k划分后的每一块与上述滤波器进行点乘，获得去方位一次空变性后的结果。将每一块均这样处理，并按原位置拼接起来，获得输出信号S_{pda_out,k}。

步骤5-8，方位向一致压缩；构造方位一致压缩滤波器：

其中r_{2_PD}为经过部分去斜后的斜距二次项系数，为r_{2_PD}＝r″_2,m,k/(1-χ)。将上述信号变换到频域，并与上式进行相乘，并进行方位向傅里叶反变换，即完成方位一致压缩，此时该距离单元内的目标点已经被精确聚焦，但是由于进行第一次方位时间变标造成了距离向位置偏移，因此聚焦之后还要进行几何校正。

步骤5-9，几何校正；

在本发明中，记一致压缩并变至进行方位向傅里叶反变换后的矩阵为S_ar,k，对其进行距离向傅里叶变换，记为S_geocor,k，构造几何校正滤波器:

其中将H_geoCorrect与S_geocor,k逐点相乘，并进行距离向傅里叶反变换，提取出第列，存储至S_image的第k列，即完成了原始回波矩阵当中第k列信号的成像。

frr₁表示距离向分割子块的第1个频点。

frr_y表示距离向分割子块的第y个频点，y表示子块频点标识号。

表示距离向分割子块的最后一个频点。

步骤5-10，遍历所有的y，即可完成全场景的精确聚焦。

P₀″表示线性校正之后第0阶斜距系数的距离空变参数。P₁″表示线性校正之后第1阶斜距系数的距离空变参数。P₂″表示线性校正之后第2阶斜距系数的距离空变参数。Q₀″表示线性校正之后的第0阶斜距系数的方位空变参数。Q₁″表示线性校正之后的第1阶斜距系数的方位空变参数。Q₂″表示线性校正之后的第2阶斜距系数的方位空变参数。r″_0,m,k表示经过方位频率变标后的第0阶第m行第k列的阶斜距系数。r″_1,m,k表示经过方位频率变标后的第1阶第m行第k列的阶斜距系数。r″_2,m,k表示经过方位频率变标后的第2阶第m行第k列的阶斜距系数。r″_3,m,k表示经过方位频率变标后的第3阶第m行第k列的阶斜距系数。r″_4,m,k表示经过方位频率变标后的第4阶第m行第k列的阶斜距系数。r″_5,m,k表示经过方位频率变标后的第5阶第m行第k列的阶斜距系数。p″_n,0表示属于P_n″中的第0阶系数。p″_n,1表示属于P″_n中的第1阶系数。p″_n,2表示属于P_n″中的第2阶系数。p″_n,3表示属于P_n″中的第3阶系数。p″_n,4表示属于P_n″中的第4阶系数。q″_n,1,k表示属于Q_n″中的第1阶第k列系数。q″_n,2,k表示属于Q_n″中的第2阶第k列系数。q″_n,3,k表示属于Q_n″中的第3阶第k列系数。q″_n,4,k表示属于Q_n″中的第4阶第k列系数。Q″_0,k表示所属Q₀″的第k列方位向空变参数。Q″_0,k表示经过方位时间变标后的第0阶第k列方位向空变参数。Q″_1,k表示所属Q₁″的第k列方位向空变参数。Q″′_1，k表示经过方位时间变标后的第1阶第k列方位向空变参数。Q″_2,k表示所属Q₂″的第k列方位向空变参数。Q″′_2,k表示经过方位时间变标后的第2阶第k列方位向空变参数。Q″_3,k表示所属Q₃″的第k列方位向空变参数。Q″′_3,k表示经过方位时间变标后的第3阶第k列方位向空变参数。Q″_4,k表示所属Q₄″的第k列方位向空变参数。Q″′_4,k表示经过方位时间变标后的第4阶第k列方位向空变参数。Q″_5,k表示所属Q₅″的第k列方位向空变参数。Q″′_5,k表示经过方位时间变标后的第5阶第k列方位向空变参数。A₁′表示经过方位时间变标后的第1阶驻定项相位点系数。A₂′表示经过方位时间变标后的第2阶驻定项相位点系数。A₃′表示经过方位时间变标后的第3阶驻定项相位点系数。A₄′表示经过方位时间变标后的第4阶驻定项相位点系数。A₅′表示经过方位时间变标后的第5阶驻定项相位点系数。A₆′表示经过方位时间变标后的第6阶驻定项相位点系数。A₇′表示经过方位时间变标后的第7阶驻定项相位点系数。AMP表示二维频谱幅度补偿系数。表示第i阶波束穿越时刻的相位系数。表示第0阶波束穿越时刻的相位系数。表示第1阶波束穿越时刻的相位系数。表示第2阶波束穿越时刻的相位系数。χ表示部分去斜因子，具体应根据方位向时宽带宽积计算而得，典型值为0.95。表示方位向波束穿越时刻的相位。

实施例1

仿真实施参数以表1所示，场景中心设置在东经103.4°，北纬31°，为我国汶川所在地。本方法用于同步轨道斜视观测，主要目的是对感兴趣目标进行长时间监控，如地质灾害地，因而对场景大小的要求不大，为此仿真场景设置为30km(方位向)×30km(距离向)。另外需要注意的是，由于RCM随着位置而发生变化，本方法能够实现的幅宽大小也是随着观测位置而变的。例如，斜视角较小时，场景中心设置在赤道处，则观测幅宽可达到150km(方位向)×150km(距离向)。

表1雷达参数

参数	取值
		轨道倾角(度)	60
轨道偏心率	0
		轨道半长轴(米)	42164431.345
近心点角距(度)	0
		升交点赤经(度)	100
过近心点时刻(秒)	0
		仿真中心时刻(秒)	9663.3421
波长λ(m)	0.24
		方位向采样点数	36000
重频PRF(赫兹)	80
		发射信号脉宽(微秒)	8
发射信号带宽(兆赫兹)	80
		天线下视角(度)	2.5882
天线斜视角(度)	1.6642

由于存在复杂的二维空变性，全场景的回波严格的讲不是完全对称的，因此选择场景中心及边缘八个点作为考核位置，这九个目标点都能精确聚焦才能说明成像算法性能良好。为此，场景按照图5所示的方式进行布设，其中T₅位于场景中心。

根据图5的场景布设方式，对T₁至T₉这九个目标点分别成像，其成像结果如图5A1～图5I3所示，其方位向剖面和距离向剖面评估结果如表2所示。

表2目标点成像质量评估结果

从图中可以看出，场景中心及八个边缘点的两个方向的剖面十分接近理想的sinc函数，由此可见此成像算法能够对全场景进行精确聚焦。另外，从全场景的质量评估表中可以看出，全场景目标点的方位分辨率波动在±0.21dB以内，方位峰值旁瓣比波动在±0.035dB以内，方位积分旁瓣比波动在±0.173dB以内，距离分辨率几乎无波动，距离峰值旁瓣比波动在±0.020dB以内，距离积分旁瓣比波动在±0.097dB以内。所有指标的全场景波动都很小，说明此算法具有良好的全场景一致性。

Claims (6)

1.一种针对星上去走动与二维空变校正的地球同步轨道SAR斜视成像处理装置，其特征在于：该装置包括有雷达参数与空变性获取模块、回波去走动模块、方位时域变标模块、徙动校正与脉冲压缩模块和方位向处理模块；

雷达参数与空变性获取模块，用于为整个成像处理系统输入参数，该模块的输出参数分为两类，一类是雷达工作参数f₁＝{Na,Nr,PRF,Fs,Br,Tp,Kr,λ,td,t_center,ta,fa,tr,fr}，，一类是回波各像素点对应地面目标点的斜距系数的空变参数，f₂＝{P_n,Q_n,0≤n≤5}；

回波去走动模块，该模块在星上完成，通过对接收回波按照后续给出的接收窗调整策略去场景中心走动，然后频谱搬移将方位向信号从射频带搬移到基带，获得去走动后的回波信号；将该回波作为后续处理操作的输入；

方位时域变标模块，该模块对上述回波进行方位时间变标，用来校正回波信号的方位向一次空变性；该模块接收雷达参数输出模块的全部输出参数，并将方位时间变标后的信号作为该模块的输出；

徙动校正与脉冲压缩模块，该模块接收雷达参数输出模块以及方位时间变标模块的全部输出参数，进行距离单元徙动校正及距离脉冲压缩，并将距离单元徙动校正及距离压缩后的回波信号作为输出参数；

方位向处理模块，该模块利用部分去斜处理技术将方位向的时域跨度缩短，然后方位向划分子孔径去除多普勒中心频率，之后利用时频混合变标去除信号的方位一次和二次空变性，以实现方位向信号聚焦；最后对信号进行几何校正，输出最终图像。

2.根据权利要求1所述的一种针对星上去走动与二维空变校正的地球同步轨道SAR斜视成像处理装置，其特征在于斜距系数的空变性建模方法如下：

步骤1-1，从中提取出第Na/2+1行的行－列－向量利用MATLAB当中的Polyfit函数，将所述对距离时间tr-2R_ref/c进行多项式拟合，得到关于tr-2R_ref/c的各阶多项式系数P_n＝[p_n,4 p_n,3 p_n,2 p_n,1 p_n,0]，即是第n阶斜距系数的距离空变系数；

步骤1-2，依次提取出的每一列的列－行－向量利用MATLAB当中的Polyfit函数，将所述对方位时间ta进行多项式拟合，得到关于ta的各阶多项式系数[q_n,4,kq_n,3,k q_n,2,k q_n,1,k q_n,0,k]，即是第n阶第k位置斜距系数的方位空变系数；

步骤1-3，将[q_n,4,k q_n,3,k q_n,2,k q_n,1,k q_n,0,k]存入尺寸为Nr×5的空变系数矩阵Q_n的第k行当中；

步骤1-4，从1至Nr遍历第k列并重复执行步骤1-1至步骤1-3，即可得到Q_n的所有元素；

步骤1-5，使n从1至5遍历，并重复执行步骤1-4，即可得到各阶斜距系数的距离空变系数矩阵{P_n,0≤n≤5}以及方位空变系数{Q_n,0≤n≤5}。

3.根据权利要求1所述的一种针对星上去走动与二维空变校正的地球同步轨道SAR斜视成像处理装置，其特征在于：接收窗调整策略应满足公式

4.根据权利要求1所述的一种针对星上去走动与二维空变校正的地球同步轨道SAR斜视成像处理装置，其特征在于方位时域变标的方法如下：

步骤3-1，将回波信号S₁做距离向傅里叶变换；

提取S₁的第一行，将此行矩阵的前个元素与后个元素对换位置，得到交换位置后的行向量g_交换；然后对所述g_交换进行快速傅里叶变换FFT，得到频谱G_交换；最后对所述G_交换进行前个元素与后个元素对换位置，即完成对回波信号S₁第一行数据的距离向傅里叶变换；

步骤3-2，构造方位时间变标滤波器

其中A_n表示驻定项相位点系数，所述A_n为

步骤3-3，将距离傅里叶变换后的回波信号S_{rng_fft}与H_AS1逐点相乘，即可完成方位时间变标，将变标后的信号第一组输出参数，记为S₂；

方位时间变标后的距离空变系数P′_n＝P_n+[0 0 0 0 A_n]，0≤n≤5；

方位时间变标后的方位空变系数

将{P′_n,0≤n≤5}、{Q′_n,0≤n≤5}以及参数{A_n,1≤n≤7}作为方位时域变标模块的第二组输出参数，记为f₃＝{P′_n,Q′_n,A_n}。

5.根据权利要求1所述的一种针对星上去走动与二维空变校正的地球同步轨道SAR斜视成像处理装置，其特征在于徙动校正及距离压缩的操作如下：

步骤4-1，计算任意一个目标点的频谱相位；

步骤4-2，计算频谱相位系数的空变参数；

步骤4-3，补偿高次相位；

步骤4-4，非线性Chirp变标；所述非线性变标方程为

步骤4-5，距离压缩及徙动校正；距离压缩及徙动校正滤波器为

6.根据权利要求1所述的一种针对星上去走动与二维空变校正的地球同步轨道SAR斜视成像处理装置，其特征在于方位向处理模块的操作为：；

步骤5-1，构造滤波器恢复原始信号的斜距模型；

步骤5-2，构造线性校正函数；

步骤5-3，构造时频混合变标—频域变标滤波器

步骤5-4，计算频率变标后的斜距系数空变系数；

步骤5-5，构造时频混合变标—时域变标滤波器

步骤5-6，频域部分去斜，构建部分去斜和高阶相位补偿滤波器

步骤5-7，分子场景去除方位向一次空变性；

步骤5-8，方位向一致压缩；构造方位一致压缩滤波器：

步骤5-9，几何校正；

步骤5-10，遍历所有的y，即可完成全场景的精确聚焦。