CN102073036B - 分布式卫星sar系统非平行航迹补偿方法 - Google Patents

Abstract

本发明公开了一种分布式卫星SAR系统非平行航迹补偿方法，主要解决分布式卫星SAR系统中非平行航迹补偿及沿航迹基线的计算的问题。其实现过程为：采用总体最小二乘法对星历位置测量数据拟合，得到成像时间段内卫星航迹的时间参数方程；对航迹位置重采样，得到每一个方位脉冲发射时间的航迹位置坐标；构造平行于主星航迹的新辅星航迹的时间参数方程；对辅星回波数据补偿，得到等效为沿平行于主星航迹的新辅星航迹接收的回波；计算沿航迹基线长度。本发明对星历位置测量数据的拟合能够克服具有短时相关性的GPS测量误差的影响，能为GMTI杂波抑制提供相关性更高的图像对及为运动目标测速定位提供准确的沿航迹基线。

Description

分布式卫星SAR系统非平行航迹补偿方法

技术领域

本发明属于目标检测技术领域，更进一步涉及一种分布式卫星SAR系统非平行航迹补偿方法，可以应用于分布式卫星合成孔径雷达地面动目标检测SAR-GMTI，该方法可以通过对非平行航迹的补偿，提高分布式卫星合成孔径雷达地面动目标检测所需复图像对的相关性，还可以为运动目标测速定位提供准确的沿航迹基线。

背景技术

在进行分布式卫星SAR-GMTI处理时，需要对接收星的成像结果进行高程干涉相位补偿，使得各通道空域导向近似相同。为了克服分布式卫星雷达系统中，各卫星航迹不平行给高程干涉相位补偿带来的影响，提高各通道SAR图像的相关性；同时为了克服非平行航迹、混合基线和非严格正侧视给GMTI沿航迹基线精确计算带来的困难，需要对分布式卫星非平行航迹进行补偿，补偿主要进行非平行航迹的数据域补偿和沿航迹基线的计算。

对于分布式卫星非平行航迹的数据域补偿，王彤等人曾在文章“分布式小卫星干涉高程测量”(《系统工程和电子技术》2004年7月26卷第7期第859-862页)中提出了一种方法。该方法根据理想平行航迹与实际航迹对地观测波程差，对回波数据进行快时间平移变换，从而实现非平行航迹补偿。由于该方法在计算理想平行航迹和实际航迹对地观测波程差中，所采用的公式是在理想平行航迹处于过真实航迹垂直于地面的平面中的几何关系下推导出来的，而在进行分布式卫星SAR-GMTI处理时，往往无法出现公式推导中需要的理想状态下的几何关系，导致无法实现该方法的补偿效果。

对于沿航迹基线的计算，余慧等人曾在文章“分布式卫星ATI-SAR性能分析与基线估计”(《电子与信息学报》2009年6月31卷第6期第1301-1304页)中，提出了一种基于回波数据的基线估计方法，通过对接收回波信号的特征分解，得到两路信号的干涉相位，并根据干涉相位与多普勒频率的线性关系解算出星间基线长度。但是，该方法仅适用于分布式卫星构型只具有沿航迹基线情况，对同时存在沿航迹基线和垂直航迹基线的情况并不适用。

发明内容

本发明的目的在于克服现有技术的不足，从分布式卫星SAR系统的观测几何构型出发，提供一种对非平行航迹补偿方法，实现非平行航迹的数据域补偿，提高SAR图像对的相关性，同时获得精确的沿航迹基线计算结果。

本发明对工作在一发两收模式下的双星系统补偿的具体步骤包括如下：

(1)拟合航迹方程。处理系统对卫星下传的主、辅星GPS星历位置测量数据采用总体最小二乘法拟合主、辅星成像时间段内航迹时间参数方程。

(2)航迹位置重采样。将成像时间段内的每一个方位脉冲发射时间，代入步骤(1)获得的航迹时间参数方程中，得到以脉冲重复周期为采样周期的航迹位置坐标。

(3)构造新辅星航迹

3a)在步骤(1)获得的主星时间参数方程中，选取时间参数的一次及一次以上项的系数，将其确定为新辅星航迹时间参数方程中的对应项系数。

3b)在步骤(1)获得的辅星时间参数方程中，选取常数项，将其确定为新辅星航迹时间参数方程的常数项。

(4)补偿辅星回波数据。对成像时间段内每一方位脉冲发射时刻，逐一分别进行步骤4a)和步骤4b)操作，补偿回波数据，实现对整个成像段回波数据的非平行航迹数据域补偿。

4a)计算波束中心视线方向单位矢量。建立波束中心视线矢量的几何关系模型，构造t＝τ时刻波束中心视线矢量与其他矢量几何关系的方程组

其中：b(τ)为t＝τ时刻辅星波束中心视线方向单位矢量，

a(τ)为t＝τ时刻由辅星指向地心的矢量，

v(τ)为t＝τ时刻辅星的速度矢量，

θ为天线斜视角，

为波束中心的下视角，

求解方程组得到t＝τ时刻的辅星波束中心视线方向单位矢量b(τ)。

4b)计算观测波程差，补偿辅星回波数据。将t＝τ时刻新辅星航迹与原辅星航迹的位置差矢量，向步骤4a)中b(τ)投影，得到观测波程差δR(τ)，并对辅星回波数据进行快时间平移补偿。

(5)计算沿航迹基线。建立补偿为平行航迹的分布式卫星的观测几何模型，计算沿航迹基线长度。

本发明与现有技术相比具有以下优点：

第一，由于本发明采用总体最小二乘法对星历位置测量数据拟合，能够克服具有短时相关性的GPS测量误差的影响；

第二，由于本发明在步骤(4)中建立分布式卫星非平行航迹的观测几何模型时，并没有限定星间的相对构型，因此相比现有技术，本发明对非平行航迹的数据域补偿处理方法不受观测几何构型的限制，在任意GMTI工作的观测时刻都可以进行非平行航迹的补偿处理，有利于提高主、辅星SAR图像对的相关性，从而获得具有更好的GMTI性能；

第三，由于本发明在步骤(5)中建立计算沿航迹基线几何模型时，考虑到了存在混合基线和雷达非严格正侧视的情况，因此相比现有技术，本发明对沿航迹基线的计算，不受基线构型和雷达视角的影响，能为动目标测速定位提供精确的沿航迹基线参数，利于获得更准确的测速定位结果。

附图说明

图1是本发明的非平行航迹补偿方法流程图。

图2是本发明波束中心视线矢量的几何关系示意图。

图3是本发明航迹平行的分布式卫星的观测几何模型示意图。

图4是未采用本发明补偿的配准后主、辅星SAR图像干涉相位图。

图5是采用本发明补偿的配准后主、辅星SAR图像干涉相位图。

具体实施方式

下面结合附图对本发明做进一步的描述。

参照图1，本发明的具体实施步骤如下：

步骤1.拟合航迹方程。

主、辅星位置测量数据是地心惯性坐标系下卫星X、Y、Z三个位置坐标的GPS测量值。在很短的成像时间内，卫星的运动航迹可以采用以时间t为参数的低阶多项式模型近似。多项式各项系数的求解，可以通过处理系统对GPS位置测量数据采用常用的最小二乘法拟合得到。实际中考虑到GPS测量误差在短时间内具有相关性，因此采用总体最小二乘法能够获得更好的拟合性能。

以主星Y坐标时间参数方程的拟合为例。将成像的中间时刻取为t＝0时刻，选定拟合阶数为三阶，则参数方程为y_m(t)＝a_ymt³+b_ymt²+c_ymt+d_ym，其中y_m(t)为主星航迹上t时刻的Y位置坐标分量，a_ym，b_ym，c_ym，d_ym为t的多项式的各阶系数，是所需求解的值，令k＝[a_ym，b_ym，c_ym，d_ym]^T，其中[·]^T表示转置。建立方程组Ak＝y，其中 $A=\left[\begin{array}{cccc}{t}_1^3& t_1^2& t_1& 1\\ t_2^3& t_2^2& t_2& 1\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ t_N^3& t_N^2& t_N& 1\end{array}\right]$ 为时间参数矩阵，N表示成像时间段内GPS采样点数，t₁，t₂，…，t_N

表示每个采样时间；y＝[y_m1，y_m2，…y_mN]^T为主星Y坐标GPS测量值构成的矢量，其中y_m1，y_m2，…y_mN为每个采样时间主星Y坐标的测量值。对Ak＝y方程组采用总体最小二乘法求解向量k，得到主星Y坐标时间参数方程的各项系数a_ym，b_ym，c_ym，d_ym。

采用与拟合主星Y坐标时间参数方程相同的方法，分别获得主、辅星航迹各坐标分量的时间参数方程的各项系数，从而构成主星航迹各坐标的时间参数方程和辅星航迹各坐标的时间参数方程

$\left\{\begin{array}{c}{x}_m\left(t\right)=a_{\mathrm{xm}}{t}^3+b_{\mathrm{xm}}{t}^2+c_{\mathrm{xm}}t+d_{\mathrm{xm}}\\ y_m\left(t\right)=a_{\mathrm{ym}}{t}^3+b_{\mathrm{ym}}{t}^2+c_{\mathrm{ym}}t+d_{\mathrm{ym}}\\ z_m\left(t\right)=a_{\mathrm{zm}}{t}^3+b_{\mathrm{zm}}{t}^2+c_{\mathrm{zm}}t+d_{\mathrm{zm}}\end{array}\right.$

$\left\{\begin{array}{c}{x}_s\left(t\right)=a_{\mathrm{xs}}{t}^3+b_{\mathrm{xs}}{t}^2+c_{\mathrm{xs}}t+d_{\mathrm{xs}}\\ y_s\left(t\right)=a_{\mathrm{ys}}{t}^3+b_{\mathrm{ys}}{t}^2+c_{\mathrm{ys}}t+d_{\mathrm{ys}}\\ z_s\left(t\right)=a_{\mathrm{zs}}{t}^3+b_{\mathrm{zs}}{t}^2+c_{\mathrm{zs}}t+d_{\mathrm{zs}}\end{array}\right.$

其中：x_m(t)、y_m(t)、z_m(t)为主星航迹上t时刻X、Y、Z三个坐标分量，

a_xm，b_xm，c_xm，d_xm为主星X坐标时间参数方程的各阶系数，

a_ym，b_ym，c_my，d_my为主星Y坐标时间参数方程的各阶系数，

a_zm，b_zm，c_zm，d_zm为主星Z坐标时间参数方程的各阶系数，

x_s(t)，y_s(t)，z_s(t)为辅星航迹上t时刻X、Y、Z三个坐标分量，

a_xs，b_xs，c_xs，d_xs为辅星X坐标时间参数方程的各阶系数，

a_ys，b_ys，c_ys，d_ys为辅星Y坐标时间参数方程的各阶系数，

a_zs，b_zs，c_zs，d_zs为辅星Z坐标时间参数方程的各阶系数。

步骤2.航迹位置重采样。

非平行航迹数据域补偿时需要对每一个方位脉冲发射时间的回波分别进行处理，而GPS星历位置坐标采样频率远小于雷达发射脉冲重复频率，不能满足处理要求。为了得到每一个方位脉冲发射时间对应的辅星位置坐标，需要对成像段内时间重新取值。将成像的中间时刻取为t＝0时刻，然后以脉冲重复周期为采样周期对成像时间段内时间重采样，获得每一个方位脉冲发射时间，再代入到步骤1中得到的辅星航迹时间参数方程里，获得每个脉冲发射时刻的辅星位置坐标。

步骤3.构造新辅星航迹。

为了过成像过程的中点时刻(t＝0时刻)辅星所在的空间位置，构造平行于主星航迹的新辅星航迹时间参数方程，将步骤1获得的主星航迹时间参数方程中t的一次项和一次以上各项系数a_xm、b_xm、c_xm、a_ym、b_ym、c_ym、a_zm、b_zm、c_zm，作为新辅星航迹时间参数方程的对应项系数，将步骤1获得的辅星航迹时间参数方程中常数项d_xs、d_ys、d_zs，作为新辅星航迹时间参数方程的常数项。得到新辅星航迹的时间参数方程如下表示

$\left\{\begin{array}{c}{x}_{\mathrm{ps}}\left(t\right)=a_{\mathrm{xm}}{t}^3+b_{\mathrm{xm}}{t}^2+c_{\mathrm{xm}}t+d_{\mathrm{xs}}\\ y_{\mathrm{ps}}\left(t\right)=a_{\mathrm{ym}}{t}^3+b_{\mathrm{ym}}{t}^2+c_{\mathrm{ym}}t+d_{\mathrm{ys}}\\ z_{\mathrm{ps}}\left(t\right)=a_{\mathrm{zm}}{t}^3+b_{\mathrm{zm}}{t}^2+c_{\mathrm{zm}}t+d_{\mathrm{zs}}\end{array}\right.$

其中：x_ps(t)，y_ps(t)，z_ps(t)为新辅星航迹上t时刻X、Y、Z三个坐标分量。

在本发明中，根据上述的新辅星航迹时间参数方程，采用与步骤2相同的方法，即可获得每个脉冲发射时刻的新辅星位置坐标。

步骤4.补偿辅星回波数据。

以t＝τ时刻为例，对每个方位脉冲发射时刻回波数据进行如下补偿：

4a)计算波束中心视线方向单位矢量。

波束中心视线矢量的几何关系示意如图2所示，图中O点为坐标系原点(地心)，S点为辅星雷达位置，v为辅星速度方向，

为波束中心的下视角，θ为天线斜视角，b(b′)为待求的波束中心视线方向单位矢量。

将步骤2中获得的t＝τ时刻辅星位置S坐标(s_x(τ)，y_s(τ)，z_s(τ))与坐标系原点相连，得到辅星指向地心的矢量a＝(-x_s(τ)，-y_s(τ)，-z_s(τ))；将步骤1中获得的辅星时间参数方程对时间t求导，得到辅星速度方向矢量v(τ)＝(v_x(τ)，v_y(τ)，v_z(τ))，其中 $v_x\left(\mathrm{\τ}\right)={\frac{{\∂x}_s\left(t\right)}{\∂t}|}_{t=\mathrm{\τ}},$ $v_y\left(\mathrm{\τ}\right)={\frac{{\∂y}_s\left(t\right)}{\∂t}|}_{t=\mathrm{\τ}},$ $v_z\left(\mathrm{\τ}\right)={\frac{{\∂z}_s\left(t\right)}{\∂t}|}_{t=\mathrm{\τ}};$ θ、

和天线左右视类型可从雷达参数中读出。根据上述已知条件和图2的几何关系，列出求解波束中心视线单位矢量b(τ)＝(b_x(τ)，b_y(τ)，b_z(τ))方程组如下：

对于上述非线性方程组，由图2所示的几何关系可以看出，将会得到两个解b和b′，分别位于卫星速度方向的左右两侧，其中只有一个解是真实的波束视线方向单位矢量。为了判断真实的波束视线方向单位矢量，需要处理系统从雷达参数中的读出雷达左右视类型，若雷达工作在左视状态，则卫星速度方向左侧的解是波束视线方向单位矢量；若雷达工作在右视状态，则卫星速度方向右侧的解是波束视线方向单位矢量。

4b)计算观测波程差，补偿辅星回波数据

在步骤2中可以得到t＝τ时刻辅星位置的坐标S(τ)，在步骤3中可以得到t＝τ时刻新辅星位置的坐标S′(τ)。新辅星航迹相对原辅星航迹的观测波程差δR(τ)，就是新辅星航迹与原辅星航迹的位置差矢量

在波波束中心视线方向单位矢量b(τ)上的投影，即

其中

将该时刻辅星接收的回波乘以

的相位因子，即获得在新辅星航迹上收到的等效回波，其中λ表示雷达工作波长。

对成像段内每个方位脉冲发射时刻t₁，t₂，…，的辅星回波数据，首先按照步骤4a)计算波束中心视线矢量，再按照步骤4b)计算观测波程差，补偿回波数据，实现了对整个成像段回波数据的非平行航迹数据域补偿，其中Na为方位脉冲个数。

步骤5.计算沿航迹基线。

图3给出了本发明在分布式卫星构型下计算等效沿航迹基线的几何示意图。对于分布式卫星系统，逐一通过步骤1、步骤2、步骤3、步骤4，将各卫星的非平行航迹补偿为平行航迹。补偿后的主星航迹L1与辅星航迹L2的空间位置虽然平行，但不重合，如图3所示，不仅具有沿航迹基线，同时具有垂直航迹基线。两种基线的同时存在，使得基于回波数据的基线估计方法在分布式卫星系统中无法适用。

当雷达工作在非严格正侧视条件下，如果把沿航迹基线等同于星间距离矢量向卫星速度方向上的投影，即使斜视角很小，由于星间距离通常具有数百米的长度，也会有较大的误差。

为了在混合基线和非严格正侧视情况下，精确计算沿航迹基线，在本发明中，采用下式计算沿航迹基线D：

其中：S₁和S₂分别为同一时刻主、辅星位置坐标，表示它们连线的矢量，

为卫星速度方向的单位矢量，

θ为雷达斜视角，波束中心视线方向单位矢量b与

的夹角，

b_s为主星航迹L₁与主星波束中心视线方向单位矢量b构成的平面内，主星正侧视情况下波束中心视线方向单位矢量。

S₁、S₂和

可以从卫星星历参数中读出，θ可以从雷达参数中读出，b_s由如下方程计算获得：

$\left\{\begin{array}{c}<b_s,\stackrel{\&RightArrow;}{V_S}>=\mathrm{\&pi;}/2\\ <b_s,b>=|\mathrm{\&pi;}/2-\mathrm{\&theta;}|\\ \left|b_s\right|=1\end{array}\right.$

从图3几何关系可以看出此方程组有唯一解，其中b的计算，可以采用与步骤4b)中相同的方法。

本发明的效果可以通过下述仿真实验加以说明：

仿真条件

雷达工作波长λ＝0.03m，轨道高度H＝512km，场景中心斜距R₀＝680km，雷达平台运动速度v_s＝7600m/s，雷达脉冲重复频率3500Hz，双星在绕飞构型的轨道上飞行，处于非平行航迹的工作状态。场景为500m×500m的平坦地面，杂波等效后向散射系数满足复高斯分布，杂噪比10dB。

仿真结果

图4是未经过非平行航迹补偿，配准后的主、辅星图像干涉相位图。图5是经过非平行航迹补偿，配准后的主、辅星图像干涉相位图。图中横坐标表示图像的方位像素位置，纵坐标表示图像的距离像素位置，仿真结果的分析如下：

图5显示的经过本发明非平行航迹的补偿后，获得的配准的主、辅星成像结果的干涉相位图较图4更为清晰。

对配准、去平地相位后得到的图像对的相关性进行统计得到：未经过非平行航迹补偿获得的图像对相关系数为0.9152；经过非平行航迹补偿获得的图像对相关系数为0.9584。由此得出，经过本发明非平行航迹补偿后，可以提高SAR图像对的相关性。根据最大杂波抑制性能计算公式

其中ρ表示相关系数，可以求得未经过非平行航迹补偿最大杂波抑制性能为7.8939dB；经过非平行航迹补偿之后最大杂波抑制性能可达10.8901dB，相比未经过非平行航迹补偿结果，杂波抑制性能明显提高。

为进一步验证本发明的有效性，在上述场景中均匀放置的157个点目标，所有目标具有相同的径向速度9.7647m/s，信杂比5dB，实验结果如下：

未经过非平行航迹补偿的杂波抑制后，检测目标数为97个；而经过非平行航迹补偿的杂波抑制后，检测目标数为115个。

当估计得到运动目标干涉相位φ′₂₁后，可以根据动目标径向速度估计公式

和方位向定位偏移量公式得到目标的径向速度，并对其重新定位，其中：v_r为目标径向速度；v为平台速度；λ为雷达工作波长；D为沿航迹基线长度；Δa为方位向定位偏移量；r为目标斜距。采用本发明方法计算得到的沿航迹基线，对所有检测到的动目标径向速度估计的均值为9.7657m/s，均方根误差为0.1691m/s，方位向相对定位均方根误差为18.3303m；以星间相对距离矢量投影到卫星速度方向的长度作为沿航迹基线，对所有检测到的动目标径向速度估计的均值为10.6299m/s，均方根误差为0.8853m/s，方位向相对定位均方根误差为95.9658m。

实验结果可以看出，采用本发明方法进行非平行航迹数据域补偿后，可以获得更好的动目标检测性能；采用本发明方法计算得到的沿航迹基线作为测速定位的参数，可以得到更接近无偏的测速定位估计，利于获得高精度的测速定位性能。

Claims (4)

1.一种分布式卫星SAR系统非平行航迹补偿方法，包括如下步骤：

(1)拟合航迹方程：处理系统对卫星下传的主、辅星GPS星历位置测量数据采用总体最小二乘法拟合主、辅星成像时间段内航迹时间参数方程；

(2)航迹位置重采样：将成像时间段内的每一个方位脉冲发射时间，代入步骤(1)获得的航迹时间参数方程中，得到以脉冲重复周期为采样周期的航迹位置坐标；

(3)构造新辅星航迹

3a)在步骤(1)获得的主星航迹时间参数方程中，选取时间参数的一次及一次以上项的系数，将其确定为新辅星航迹时间参数方程中的对应项系数；

3b)在步骤(1)获得的辅星航迹时间参数方程中，选取常数项，将其确定为新辅星航迹时间参数方程的常数项；

(4)补偿辅星回波数据：对成像时间段内每一方位脉冲发射时刻，逐一分别进行步骤4a)和步骤4b)操作，补偿辅星回波数据，实现对整个成像时间段回波数据的非平行航迹数据域补偿；

4a)计算辅星波束中心视线方向单位矢量：建立辅星波束中心视线方向单位矢量的几何关系模型，构造t＝τ时刻辅星波束中心视线方向单位矢量与其他矢量几何关系的方程组

其中：b(τ)为t＝τ时刻辅星波束中心视线方向单位矢量，

a(τ)为t＝τ时刻由辅星指向地心的矢量，

v(τ)为t＝τ时刻辅星的速度矢量，

θ_s为天线斜视角，

为辅星波束中心的下视角，

求解方程组得到t＝τ时刻的辅星波束中心视线方向单位矢量b(τ)；

4b)计算观测波程差，补偿辅星回波数据：将t＝τ时刻新辅星航迹与原辅星航迹的位置差矢量，向步骤4a)中b(τ)投影，得到观测波程差δR(τ)，并对辅星回波数据进行快时间平移补偿；

(5)计算沿航迹基线：建立补偿为平行航迹的分布式卫星的观测几何模型，计算沿航迹基线长度。

2.根据权利要求1所述的分布式卫星SAR系统非平行航迹补偿方法，其特征在于：所述步骤(1)中时间参数方程的拟合是在对主、辅星每个坐标建立方程组Ak＝y的基础上，采用总体最小二乘法求解得到时间参数方程中各项系数，其中A为时间参数矩阵，y为由所取坐标GPS测量值构成的矢量，k为由时间参数方程中各项系数构成的矢量。

3.根据权利要求1所述的分布式卫星SAR系统非平行航迹补偿方法，其特征在于：所述步骤(2)中方位脉冲发射时间是指以脉冲重复周期为采样周期对成像时间段内时间重采样获得的各个时间。

4.根据权利要求1所述的分布式卫星SAR系统非平行航迹补偿方法，其特征在于：所述步骤(5)中沿航迹基线长度计算采用如下公式：

其中：D表示沿航迹基线长度，

S₁、S₂分别表示同一时刻主、辅星位置坐标，

表示它们连线的矢量，

为卫星速度方向的单位矢量，

θ为雷达斜视角，是主星波束中心视线方向单位矢量b与

的夹角，

b_s为主星航迹与主星波束中心视线方向单位矢量b构成的平面内，主星正侧视情况下主星波束中心视线方向单位矢量。

Images