CN103823210A - 一种非合作式星地双基地sar时频同步方法 - Google Patents

Abstract

本发明公开了一种非合作星地双基地SAR时频同步方法，将天线接收的信号进行二维划分、脉冲压缩操作后进行数据处理，能够将接收信号中旁瓣接收的直达波信号数据提取出来，并根据短合成孔径时间内，卫星轨道近似直线的先验信息，利用斜距历史的平方具有关于时间的三次项系数为0的特点，抑制了噪声对估计精度的影响，从而实现利用天线背瓣中接收的直达波信号信噪比弱的情况下，实现星地双基地SAR系统的同步误差估计。

Description

一种非合作式星地双基地SAR时频同步方法

技术领域

本发明属于雷达信号处理技术领域，涉及一种非合作式星地双基地SAR时频同步方法。

背景技术

星地双基地SAR是利用卫星为发射平台，近地面放置接收机组成的双基地SAR系统，除了SAR本身具备的全天时、全天候观测能力，双基地SAR由于其灵活的配置，在抗干扰、前视成像、干涉处理方面有着独特的优势，因而得到了越来越多的关注。

基于星地双基地SAR系统由于收发平台分置，与收发端有相同的时间轴和本振信号的单基地SAR信号相比，双基地SAR信号存在时间和频率同步误差，会导致未知的起始采样距离、徙动曲线线性偏移和多普勒相位历史误差，从而引起最终的图像严重散焦。

对于收发同步的方法，最早是通过增加专门的通信链路来进行的。对于非合作式双基地SAR，由于不能采用通信链路，只能采用基于卫星直达波信息的方法。为了获取直达波信息，目前广泛采用的方案是，增加一个额外的专用天线指向卫星用于接收来自卫星的直达波数据，同时另一个天线指向场景，用于接收场景的回波信号。这种方案不但增加了硬件复杂度，而且存在两通道之间的一致性问题。实际上，来自卫星的直达波信号能够从回波天线旁瓣中进入，但是，由于天线旁瓣的增益小，直达波信号信噪比弱，目前还没有人利用其进行同步误差估计。

发明内容

本发明提出了一种非合作星地双基地SAR时频同步方法，能够利用天线背瓣中接收的直达波信号信噪比弱的情况下，实现同步误差的高精度估计。

一种非合作式星地双基地SAR时频同步方法，其特征在于，包括下列步骤：

步骤一、将天线的主瓣指向场景，接收场景回波，旁瓣接收来自卫星的直达波信号，当卫星波束扫过场景时，接收机连续采集数据并存储；

步骤二：对步骤一获得的数据进行脉冲压缩，获得相邻峰值之间的时间间隔PRT_R，按照PRT_R对步骤一获得的数据进行二维划分，获得以快时间τ和慢时间η表征的二维数据；

步骤三、对步骤二获得的二维数据进行距离向脉压，获得直达波的波束中心穿越时刻η0、视在斜距历史R_deo(η)和峰值位置处随慢时间变化的相位历史

步骤四、利用星历数据计算得到的卫星与天线之间的最短斜距R₂₀的估计值

步骤五：对步骤三获得的获取的视在斜距历史R_deo(η)进行三次拟合获得多普勒中心频率f_dc与划分误差ΔPRT的关系；所述划分误差ΔPRT为真实脉冲重复时间和估计的脉冲重复时间PRT_R之间的误差；

步骤六、根据步骤三获得的峰值位置处随慢时间变化的相位历史

获取多普勒中心频率f_dc与收发载频的固定频差Δf的关系；

步骤七、将ΔPRT在(-PRT_R/2,PRT_R/2)范围，

在(R₂₀-100km,R₂₀+100km)范围进行遍历，获得多普勒中心频率f_dc与起始采样距离R_m的关系；ΔPRT为划分误差ΔPRT的估计值，

为起始采样距离R_m的估计值；

步骤八、在[-1/PRT_R,1/PRT_R]范围内对f_dc遍历，得到f_dc和R₂₀的关系；

步骤九、根据步骤四获得的和步骤八得到f_dc和R₂₀的关系，获得f_dc的估计值

步骤十、根据步骤九获得的估计值

和步骤五获得的多普勒中心频率f_dc与划分误差ΔPRT的关系，获得多普勒中划分误差ΔPRT的估计值ΔPRT；根据步骤九获得的估计值和步骤六获得的多普勒中心频率f_dc与收发载频的固定频差Δf的关系，获得固定频差Δf的估计值

步骤十一、根据步骤十获得的划分误差估计值ΔPRT，固定频差估计值

和步骤二获得的时间间隔PRT_R，获得直达波信号或场景回波信号的同步误差补偿数据，具体为：

1）时间同步：对第n个脉冲，左移n·ΔPRT

2）频率同步：对第n个脉冲，乘以相位

3）根据步骤十获得的估计值

和步骤七获得的多普勒中心频率f_dc与起始采样距离R_m的关系获得起始采样距离R_m的估计值

进一步的，步骤五中获得多普勒中心频率f_dc与划分误差ΔPRT的关系具体方法为：

对步骤三获得的视在斜距历史R_deo(η)进行三次拟合获得一次项、二次项和三次项系数分别记为a₁、a₂、a₃；

将视在斜距距离向脉压后的表达式R_deo(η)在波束中心穿越时刻η₀进行泰勒展开，获得一次项、二次项和三次项系数的表达式；

根据a₁、a₂、a₃与泰勒展开的一次项、二次项、三次项系数表达式获得多普勒中心频率f_dc与划分误差ΔPRT的关系以及多普勒调频率f_dr：

$\mathrm{\ΔPRT}=(a_1+{\mathrm{\λf}}_{\mathrm{dc}}+{2a}_c{\mathrm{\η}}_0+{3a}_3{{\mathrm{\η}}_0}^2)\frac{{\mathrm{PRT}}_R}{c}---\left(34\right)$

$f_{\mathrm{dr}}=-\frac{{2a}_1+a_3{\mathrm{\η}}_0}{\mathrm{\λ}}---\left(35\right)$

c为光速，λ为信号载波波长，所述划分误差ΔPRT为真实脉冲重复时间和估计的脉冲重复时间PRT_R之间的误差；

进一步的，步骤六中获取多普勒中心频率f_dc与收发载频的固定频差Δf的关系具体为：

根据步骤三获得的峰值位置处随慢时间变化的相位历史获取视在多普勒中心频率f_{dc_view}的值，进而获得收发载频的固定频差Δf与多普勒中心频率f_dc的关系：

Δf＝f_{dc_view}-f_dc    （37）

进一步的，步骤七获得多普勒中心频率f_dc与起始采样距离R_m的关系具体方法为：

将ΔPRT在(-PRT_R/2,PRT_R/2)范围，在(R₂₀-100km,R₂₀+100km)范围进行遍历，并在每个遍历点对视在斜距历史R_deo(η)按照公式（42）进行补偿，获得每个遍历点对应的补偿斜距历史R_dec(η)，并对获得的补偿斜距历史的平方计算三次项系数α₃；其中，ΔPRT为划分误差ΔPRT的估计值，

为起始采样距离R_m的估计值；

$R_{\mathrm{dec}}\left(\mathrm{\η}\right)=R_{\mathrm{deo}}\left(\mathrm{\η}\right)-\frac{\mathrm{\ΔPRT}}{{\mathrm{PRT}}_R}\·c\·\mathrm{\η}+{\hat{R}}_m---\left(42\right)$

取（α₃，ΔPRT，

）构成的平面与（0，ΔPRT，

）构成的平面的交线的方程，即为ΔPRT与

的关系：

${\hat{R}}_m=r_0+r_1\mathrm{\ΔPRT}---\left(28\right)$

再根据步骤三获得的f_dc与ΔPRT的关系获得

与f_dc的关系为：

${\hat{R}}_m=q_0+q_1{\hat{f}}_{\mathrm{dc}}---\left(39\right)$

其中，r₁、q₁为一次项系数，r₀、q₀为常数项；

进一步的，步骤八中得到f_dc和R₂₀的关系具体方法为：

根据步骤四获得的多普勒中心频率f_dc与划分误差ΔPRT的关系和步骤六获得的多普勒中心频率f_dc与起始采样距离

的关系，将补偿后的斜距表示为：

R_dec(η;f_dc)＝R_deo(η)+k₀+k₁·f_dc    （55）

k₀、k₁为常数项和一次系数项，为已知量；

在[-1/PRT_R,1/PRT_R]范围内对f_dc遍历，并在每个遍历点按照公式（55）计算

并进行二次拟合，获得二次项系数b₀、一次项系数b₁、常数项b₂；

将真实斜距历史表达式的平方展开获得二次项系数、一次项系数和常数项的表达式；

根据b₀、b₁、b₂和真实斜距历史平方展开获得的二次项系数、一次项系数和常数项的表达式，获得f_dc与R₂₀的关系为：

$R_{20}=\sqrt{b_0-{({\mathrm{\η}}_0-\frac{f_{\mathrm{dc}}}{f_{\mathrm{dr}}})}^2\·b_2}---\left(56\right)$

有益效果

1）本发明利用天线旁瓣接收的直达波信号进行同步误差的估计和补偿，不需要单独的直达波天线，从而降低了硬件复杂度。

2）本发明将天线接收的信号进行二维划分、脉冲压缩操作后进行数据处理，能够将接收信号中旁瓣接收的直达波信号数据提取出来，从而实现星地双基地SAR系统的同步误差估计。

3）本发明利用短合成孔径时间内，卫星轨道近似直线的先验信息，该条件下，斜距历史的平方具有关于时间的三次项系数为0的特点，以该特点为约束条件，获取了起始采样距离和划分误差ΔPRT之间的关系，从而通过对先验信息的充分利用进一步抑制了噪声对估计精度的影响，为同步误差的高精度估计提供了基础。

附图说明

图1为星地双基地的系统示意图。

具体实施方式

下面结合附图对本发明进行详细描述。

本发明提供了一种非合作一站固定式星地双基地SAR的时频同步方法,具体包括下列步骤：

步骤一、将天线的主瓣指向场景，接收场景回波，旁瓣接收来自卫星的直达波信号，当卫星波束扫过场景时，接收机连续采集数据并存储；

如图1所示，建立一个地固坐标系，其中卫星以速度V运动，同时发射信号照射场景区域，T为场景中的任意点，坐标为(x,y)，A为天线。

首先，天线接收到的信号表示为：

S_r(t)＝S_d(t)+S_echo(t)+w′(t)    （1）

$\begin{array}{c}{S}_d\left(t\right)=\underset{n}{\mathrm{\Σ}}{W}_d\left(t\right)\mathrm{rect}\left(\frac{t+{\mathrm{\τ}}_m-n\·{\mathrm{PRT}}_T-{\mathrm{\τ}}_d\left(t\right)}{T_p}\right)\\ \·\exp \left({\mathrm{j\πk}}_r{(t+{\mathrm{\τ}}_m-n\·{\mathrm{PRT}}_T-{\mathrm{\τ}}_d\left(t\right))}^2\right)\\ \·\exp \left(j2\mathrm{\π}{\∫}_0^{t+{\mathrm{\τ}}_m-{\mathrm{\τ}}_d\left(t\right)}{f}_T\left(t^{\′}\right){\mathrm{dt}}^{\′}\right)\end{array}---\left(2\right)$

$S_{\mathrm{echo}}\left(t\right)\underset{x,y}{\∫\∫}\mathrm{\σ}(x,y)\left[\begin{array}{c}\underset{n}{\mathrm{\Σ}}{W}_{\mathrm{echo}}\left(t\right)\mathrm{rect}\left(\frac{t+{\mathrm{\τ}}_m-n\·{\mathrm{PRT}}_T-{\mathrm{\τ}}_{\mathrm{echo}}(t;x,y)}{T_p}\right)\\ \exp \left({\mathrm{j\πk}}_r{(t+{\mathrm{\τ}}_m-n\·{\mathrm{PRT}}_T-{\mathrm{\τ}}_{\mathrm{echo}}(t;x,y))}^2\right)\\ \·\exp \left(j2\mathrm{\π}{\∫}_0^{t+{\mathrm{\τ}}_m-{\mathrm{\τ}}_{\mathrm{echo}}(t;x,y)}\stackrel{}{}{f}_T\left(t^{\′}\right){\mathrm{dt}}^{\′}\right)\end{array}\right]\mathrm{dxdy}---\left(3\right)$

其中，t表示时间，n表示脉冲号，t′为积分变量，w′(t)为信号中的噪声，S_d(t)为直达波信号，S_echo(t)为回波信号，W_d(t)和W_echo(t)分别表示直达波信号和回波信号中收发天线的天线方向图、传播损耗造成的加权因子。σ(x,y)为位于(x,y)位置处目标的后向散射系数。PRT_T为发射脉冲的脉冲重复周期，T_p表示脉冲宽度，k_r表示发射脉冲的调频率，f_T(t′)为发射机本振的时变频率误差。

τ_m为起始采样时间：

τ_m＝R_m/c    （4）

其中，R_m表示起始采样距离，c为光速。

rect(·)为门函数，定义为：

$\mathrm{rect}\left(v\right)=\left\{\begin{array}{cc}1,& \left|v\right|\≤1\\ 0,& \left|v\right|>1\end{array}\right.---\left(5\right)$

τ_d(t)为直达波的时延:

τ_d(t)＝R_d(t)/c    （6）

其中，R_d(t)表示收发平台间的真实斜距历史。

τ_echo(t;x,y)为场景回波的时延：

${\mathrm{\τ}}_{\mathrm{echo}}(t;x,y)=\frac{R_{\mathrm{ST}}(t;x,y)+R_{\mathrm{RT}}(t;x,y)}{c}---\left(7\right)$

其中，R_ST(t;x,y)和R_RT(t;x,y)分别表示t时刻，卫星到场景目标T(x,y)以及场景目标和回波天线之间的距离。

对发射和接收本振频率分别建模，在较短的合成孔径时间下，可认为收发本振频率误差主要表现为一个恒定的偏差，则：

f_T(t)＝f_c+Δf_T+n_T(t)    （8）

f_R(t)＝f_c+Δf_R+n_R(t)    （9）

其中，f_c为本振的标称频率，Δf_T与Δf_R分别为发射和接收本振频率与标称频率的平均偏差，n_T(t)与n_R(t)是发射和接收本振频率的零均值的随机抖动。

其次，将天线接收到的信号解调到基带，得到：

S_rb(t)＝S_db(t)+S_eb(t)+w(t)    （10）

$\begin{array}{c}{S}_{\mathrm{db}}=\underset{n}{\mathrm{\Σ}}{W}_d\left(t\right)\mathrm{rect}\left(\frac{t+{\mathrm{\τ}}_m-n\·{\mathrm{PRT}}_T-{\mathrm{\τ}}_d\left(n\right)}{T_p}\right)\\ \×\exp \left({\mathrm{j\πk}}_r{(t-n\·{\mathrm{PRT}}_T-{\mathrm{\τ}}_d\left(n\right))}^2\right)\\ \×\exp [-j2\mathrm{\π}(f_c+{\mathrm{\Δf}}_T){\mathrm{\τ}}_d\left(n\right)]\exp (j2\mathrm{\π\Δf}\·t)\\ \×\exp \left(j2\mathrm{\π}{\∫}_0^t\right[n_T\left(t^{\′}\right)\left]{\mathrm{dt}}^{\′}\right)\exp \left(j2\mathrm{\π}{\∫}_t^{t+{\mathrm{\τ}}_m-{\mathrm{\τ}}_d\left(n\right)}{n}_T\left(t^{\′}\right){\mathrm{dt}}^{\′}\right)\end{array}---\left(11\right)$

$S_{\mathrm{eb}}\left(t\right)=\underset{x,y}{\∫\∫}\mathrm{\σ}(x,y)\left[\begin{array}{c}\underset{n}{\mathrm{\Σ}}{W}_{\mathrm{echo}}\left(t\right)\mathrm{rect}\left(\frac{t+{\mathrm{\τ}}_m-n\·{\mathrm{PRT}}_T-{\mathrm{\τ}}_{\mathrm{echo}}(n;x,y)}{T_p}\right)\\ \×\exp \left({\mathrm{j\πk}}_r{(t-n\·{\mathrm{PRT}}_T-{\mathrm{\τ}}_{\mathrm{echo}}(n;x,y))}^2\right)\\ \×\exp [-j2\mathrm{\π}(f_c+{\mathrm{\Δf}}_T){\mathrm{\τ}}_{\mathrm{echo}}(n;x,y)]\exp (j2\mathrm{\π\Δf}\·t)\\ \×\exp \left(j2\mathrm{\π}{\∫}_0^t\right[n_T\left(t^{\′}\right)-n_R\left(t^{\′}\right)\left]{\mathrm{dt}}^{\′}\right)\\ \×\exp \left(j2\mathrm{\π}{\∫}_t^{t+{\mathrm{\τ}}_m-{\mathrm{\τ}}_{\mathrm{echo}}(n;x,y)}{n}_T\left(t^{\′}\right){\mathrm{dt}}^{\′}\right)\end{array}\right]---\left(12\right)$

其中，S_db(t)为直达波的基带输出，S_eb(t)为回波信号的基带输出，w(t)为基带输出的噪声，Δf＝Δf_T-Δf_R，表示收发载频的固定频差。

由式（11）、（12）可知，直达波和场景回波信号中存在共同的始采样时间τ_m以及频偏Δf因子，即直达波和场景回波信号中的未知起始采样时间τ_m以及频偏Δf是一样的，因此，通过直达波信号估计得到的同步误差补偿到回波中就能完成同步处理。

另一方面，直达波和回波信号存在时延差：

$\mathrm{\Δ\τ}(t;x,y)=\frac{R_{\mathrm{ST}}(t;x,y)+R_{\mathrm{RT}}(t;x,y)-R_d\left(t\right)}{c}---\left(13\right)$

一般而言，场景距离天线的最近距离为几百米，而距离向分辨率为几米的量级，因此，通过时延差Δτ(t;x,y)能够将场景回波和直达波明显地分开，因此能够将直达波信号从天线接收的信号分离出来，本实施例在理论论证时只对直达波信号进行说明。

考虑w(t)，直达波信号的基带输出表示为：

$\begin{array}{c}{S}_{\mathrm{db}}^{\′}\left(t\right)=\underset{n}{\mathrm{\Σ}}{W}_d\left(t\right)\mathrm{rect}\left(\frac{t+{\mathrm{\τ}}_m-n\·{\mathrm{PRT}}_T-{\mathrm{\τ}}_d\left(t\right)}{T_p}\right)\\ \×\exp \left({\mathrm{j\πl}}_r{(t-n\·{\mathrm{PRT}}_T-{\mathrm{\τ}}_d\left(t\right))}^2\right)\\ \×\exp [-j2\mathrm{\π}(f_c+{\mathrm{\Δf}}_T){\mathrm{\τ}}_d\left(t\right)]\exp (j2\mathrm{\π\Δf}\·t)\\ \×\exp \left(j2\mathrm{\π}{\∫}_0^t\right[n_T\left(t^{\′}\right)-n_R\left(t^{\′}\right)\left]{\mathrm{dt}}^{\′}\right)\exp \left(j2\mathrm{\π}{\∫}_t^{t+{\mathrm{\τ}}_m-{\mathrm{\τ}}_d\left(t\right)}{n}_T\left(t^{\′}\right){\mathrm{dt}}^{\′}\right)+w\left(t\right)\end{array}---\left(14\right)$

由于合成孔径时间较短，此式相位误差中，固定频偏造成的线性相位占相位同步误差的绝大部分，此时，可忽略n_T(t′)与n_R(t′)对直达波相位的影响，则直达波信号基带输出为：

$\begin{array}{c}{S}_{\mathrm{db}}^{\′}\left(t\right)=\underset{n}{\mathrm{\Σ}}{W}_d\left(t\right)\mathrm{rect}\left(\frac{t+{\mathrm{\τ}}_m-n\·{\mathrm{PRT}}_T-{\mathrm{\τ}}_d\left(t\right)}{T_p}\right)\\ \×\exp \left({\mathrm{j\πk}}_r{(t-n\·{\mathrm{PRT}}_T-{\mathrm{\τ}}_d\left(t\right))}^3\right)\\ \×\exp (j2\mathrm{\π\Δf}\·t)\exp [-j2\mathrm{\π}(f_c+{\mathrm{\Δf}}_T){\mathrm{\τ}}_d\left(t\right)]+w\left(t\right)\end{array}---\left(15\right)$

步骤二：对步骤一获得的数据进行脉冲压缩，获得相邻峰值之间的时间间隔PRT_R，按照PRT_R对步骤一获得的数据进行二维划分，获得以快时间τ和慢时间η表征的二维数据；

首先，对基带信号进行脉冲压缩，获得相邻峰值之间的时间间隔PRT_R，按照PRT_R对步骤一获得的数据进行二维划分。

其次，通过对一维数据进行脉压获得相邻脉冲峰值之间的位置，从而获得相邻峰值之间的时间间隔PRT_R，具体的获得方法为，取距离向点数为：

N_r＝round(PRT_R·f_s)    （16）

其中，f_s为采样率，round(·)表示邻近取整操作。

由于发射平台在运动，相邻峰值间的间距会发生变化，这可能导致利用式（16）估计得到的N_r对原始数据进行划分时，出现最后某个脉冲被截断为两部分的现象，因此，需要根据划分后的二维信号幅度的分布，对N_r进行微调，如果慢时间η越大，脉冲所在区域往右边偏移，则应增大N_r，否则应减小，直到不发生脉冲的截断。一般而言，微调量不会超过2。微调完成后，获得新的距离向点数，进而更新PRT_R：

PRT_R＝N_r/f_s    （17）

本实施例仅介绍直达波部分，按照PRT_R对直达波信号基带输出S′_db(t)进行划分。设快时间为τ，则全时间为：

t＝n·PRT_R+τ    （18）

设划分误差为ΔPRT，即真实脉冲重复时间和估计的脉冲重复时间之间的误差为：

ΔPRT＝PRT_T-PRT_R    （19）

将（18）、（19）代入（15），则直达波信号基带输出表示为：

$\begin{array}{c}{S}_{\mathrm{db}}^{\′}(n,\mathrm{\τ})=W_d\left(n\right)\mathrm{rect}\left(\frac{\mathrm{\τ}+{\mathrm{\τ}}_m-n\·\mathrm{\ΔPRT}-{\mathrm{\τ}}_d\left(n\right)}{T_p}\right)\\ \×\exp \left({\mathrm{j\πk}}_r{(\mathrm{\τ}-n\·\mathrm{\ΔPRT}-{\mathrm{\τ}}_d\left(n\right))}^2\right)\\ \×\exp (j2\mathrm{\π\Δf}\·(n\·{\mathrm{PRT}}_R+\mathrm{\τ}))\exp [-j2\mathrm{\π}(f_c+\mathrm{\Δ}{f}_T){\mathrm{\τ}}_d\left(n\right)]+w(n,\mathrm{\τ})\end{array}---\left(20\right)$

可见，划分误差将导致脉冲信号在每行存在一个线性偏移，该偏移与信号的真实距离徙动曲线耦合在一起，影响距离徙动校正。

步骤三：对步骤二获得的二维数据进行距离向脉压，获得直达波的波束中心穿越时刻估计η₀、视在斜距历史R_deo(η)和峰值位置处随慢时间变化的相位历史

距离向脉压参考函数为：

$S_{\mathrm{ref}}\left(\mathrm{\τ}\right)=\mathrm{rect}\left(\frac{\mathrm{\τ}}{T_p}\right)\exp \left({\mathrm{j\πk}}_r{\left(\mathrm{\τ}\right)}^2\right)---\left(21\right)$

本实施例仅介绍直达波部分，因此，对（20）进行脉压，得到：

$\begin{array}{c}{s}_{\mathrm{db},\mathrm{rc}}(n,\mathrm{\τ})=A\left(n\right)\sin c\left[(B-\mathrm{\Δf})\right[\mathrm{\τ}+{\mathrm{\τ}}_m-(n\·\mathrm{\ΔPRT}+{\mathrm{\τ}}_d\left(n\right)-\frac{\mathrm{\Δf}}{k_r})\left]\right]\\ \×\exp (-\mathrm{j\π\Δf}[\mathrm{\τ}+{\mathrm{\τ}}_m-(n\·\mathrm{\ΔPRT}+{\mathrm{\τ}}_d\left(n\right)-\frac{\mathrm{\Δf}}{k_r})\left]\right)\\ \×\exp (-j2\mathrm{\π}(f_c+{\mathrm{\Δf}}_T)\·{\mathrm{\τ}}_d\left(n\right))\\ \×\exp (-\mathrm{j\π}[\frac{{\mathrm{\Δf}}^2}{k_r}-2(n\·\mathrm{\ΔPRT}+{\mathrm{\τ}}_d\left(n\right))\mathrm{\Δf}\left]\right)\\ \×\exp (j2\mathrm{\π\Δf}\·n\·{\mathrm{PRT}}_R)+w^{\′\′}(n,\mathrm{\τ})\end{array}---\left(22\right)$

其中，w′′(n,τ)为脉压后的w(n,τ)。

由于ΔfB,Δf²

k_r，Δf_T

f_c因此，（22）可近似为：

$\begin{array}{c}{s}_{\mathrm{db},\mathrm{rc}}(n,\mathrm{\τ})=A\left(n\right)\sin c\left[B\right[\mathrm{\τ}+{\mathrm{\τ}}_m-(n\·\mathrm{\ΔPRT}+{\mathrm{\τ}}_d\left(n\right))\left]\right]\\ \×\exp (-\mathrm{j\π\Δf}[\mathrm{\τ}+{\mathrm{\τ}}_m-(n\·\mathrm{\ΔPRT}+{\mathrm{\τ}}_d\left(n\right))\left]\right)\\ \×\exp (j2\mathrm{\π\Δf}\·n\·{\mathrm{PRT}}_R)\exp (-j2\mathrm{\π}{f}_c\·{\mathrm{\τ}}_d\left(n\right))+w^{\′\′}(n,\mathrm{\τ})\end{array}---\left(23\right)$

直达波信号经过距离向脉压后，峰值位置对应的时间延迟历史表达式为：

τ_peak(n)＝n·ΔPRT+τ_d(n)-τ_m    （24）

则视在斜距历史R_deo(η)表达式为：

$R_{\mathrm{deo}}\left(\mathrm{\η}\right)=c\·\·{\mathrm{\τ}}_{\mathrm{peak}}\left(n\right)=\mathrm{\η}\·\frac{\mathrm{\ΔPRT}}{{\mathrm{PRT}}_R}\·c+R_d\left(\mathrm{\η}\right)-R_m---\left(25\right)$

直达波信号经过距离向脉压后，峰值位置处随慢时间变化的相位历史表达式为：

取出直达波信号的每个脉冲压缩后的峰值位置处的数据，并按照方位向排列成矢量，找到该矢量的幅度最大处对应的方位向时间，即为波束中心穿越时刻，记为η₀。

步骤四：利用星历数据获得卫星与天线之间的最短斜距R₂₀的估计值

步骤五、对步骤三获得的视在斜距历史R_deo(η)进行三次拟合获得一次项、二次项和三次项系数分别记为a₁、a₂、a₃；

对视在斜距历史表达式R_deo(η)在波束中心穿越时刻η₀进行泰勒展开，获得一次项、二次项和三次项系数的表达式：

$\begin{array}{c}{R}_{\mathrm{deo}}\left(\mathrm{\η}\right)\≈R_d\left({\mathrm{\η}}_0\right)-R_m+\mathrm{\η}\·\frac{\mathrm{\ΔPRT}}{{\mathrm{PRT}}_R}\·c{R}_d^{\′}\left({\mathrm{\η}}_0\right)\left(\mathrm{\η}-{\mathrm{\η}}_0\right)\\ +\frac{{R_d}^{\left(2\right)}}{2}{(\mathrm{\η}-{\mathrm{\η}}_0)}^2+\frac{{R_d}^{\left(3\right)}\left({\mathrm{\η}}_0\right)}{6}{\left({\mathrm{\η}-\mathrm{\η}}_0\right)}^3\end{array}---\left(27\right)$

令：

$a_1=\frac{\mathrm{\ΔPRT}}{{\mathrm{PRT}}_R}\·c+R_d^{\′}\left({\mathrm{\η}}_0\right)-R_d^{\left(2\right)}\left({\mathrm{\η}}_0\right){\mathrm{\η}}_0+\frac{R_d^{\left(3\right)}\left({\mathrm{\η}}_0\right)}{2}{\mathrm{\η}}_0^2---\left(28\right)$

$a_2=\frac{R_d^{\left(2\right)}\left({\mathrm{\η}}_0\right)}{2}-\frac{R_d^{\left(3\right)}\left({\mathrm{\η}}_0\right)}{2}{\mathrm{\η}}_0---\left(29\right)$

$a_3=\frac{R_d^{\left(3\right)}\left({\mathrm{\η}}_0\right)}{6}---\left(30\right)$

由于多普勒中心频率为：

$f_{\mathrm{dc}}=-\frac{R_d^{\′}\left({\mathrm{\η}}_0\right)}{\mathrm{\λ}}---\left(31\right)$

多普勒调频率为:

$f_{\mathrm{dr}}=-\frac{{R_d}^{\left(2\right)}\left({\mathrm{\η}}_0\right)}{\mathrm{\λ}}---\left(32\right)$

λ为信号载波波长。

由（29）-（31）代入到（28）得：

$a_1=\frac{\mathrm{\ΔPRT}}{{\mathrm{PRT}}_R}\·c-{\mathrm{\λf}}_{\mathrm{dc}}-{2a}_2{\mathrm{\η}}_0-{3a}_2{\mathrm{\η}}_0^2---\left(33\right)$

即获得多普勒中心频率f_dc与划分误差ΔPRT的关系：:

$\mathrm{\ΔPRT}=(a_1+{\mathrm{\λf}}_{\mathrm{dc}}+{2a}_2{\mathrm{\η}}_0+{3a}_3{\mathrm{\η}}_0^2)\frac{{\mathrm{PRT}}_R}{c}---\left(34\right)$

将（29）（30）代入到（32），获得多普勒调频率f_dr：

$f_{\mathrm{dr}}=-\frac{{2a}_1+a_3{\mathrm{\η}}_0}{\mathrm{\λ}}---\left(35\right)$

上式表明了，ΔPRT与f_dc呈线性关系，且各项系数可由对峰值位置历史进行拟合得到。也可看出，同步误差并不影响调频率f_dr，这是因为在短合成孔径时间下，同步误差主要表现为距离徙动随方位向的线性变化和频率的固定偏移。

实际中，脉压峰值的估计不如峰值处相位的估计准确，为提高估计精度，在实际操作中，对于峰值位置历史二阶及以上的系数，由相位进行估计，并按照比例因子换算。

步骤六、根据步骤三获得的峰值位置处随慢时间变化的相位历史

获取视在多普勒中心频率f_{dc_view}的值。

具体为：取峰值位置处随慢时间变化的相位历史

在波束中心穿越时刻的导数值并除以2π，即获得视在多普勒频率的值f_{dc_view}；

设收发载频存在固定频差Δf，则视在多普勒频率：

f_{dc_view}＝Δf+f_dc    （36）

即多普勒中心频率f_dc与收发载频的固定频差Δf满足关系：

Δf＝f_{dc_view}-f_dc    （37）

在没有频率同步误差的情况下，即Δf＝0时，视在多普勒中心频率与真实多普勒中心频率相等，而当Δf非零时，传统的f_dc估计方法并不能估计出真实的f_dc。

步骤七、将ΔPRT在(-PRT_R/2,PRT_R/2)范围，在(R₂₀-100km,R₂₀+100km)范围进行遍历，并在每个遍历点对视在斜距历史R_deo(η)按照公式（42）进行补偿，获得每个遍历点对应的补偿斜距历史R_dec(η)，并对获得的补偿斜距历史的平方计算三次项系数α₃；其中，ΔPRT为划分误差ΔPRT的估计值，

为起始采样距离Rm的估计值；

取（α₃、ΔPRT、

）构成的平面与（0、ΔPRT、

）构成的平面的交线的方程，即为ΔPRT与

的关系：

${\hat{R}}_m=r_0+r_1\mathrm{\ΔPRT}---\left(38\right)$

再根据步骤三获得的f_dc与ΔPRT的关系获得

与f_dc的关系为：

${\hat{R}}_m=q_0+q_1{\hat{f}}_{\mathrm{dc}}---\left(39\right)$

其中，r₁、q₁为一次项系数，r₀、q₀为常数项；

补偿斜距历史R_dec(η)及取（α₃，ΔPRT，

）构成的平面与（0，ΔPRT、

）构成的平面的交线的方程，即为ΔPRT与

的关系的理论依据为：

当卫星与天线距离最小时，对应的多普勒频率为0，称该时刻为零多普勒时刻。因此，从零多普勒时刻到波束中心穿越时刻对应的时间为

在短合成孔径时间的条件下，卫星的飞行轨迹近似直线，因此可将真实斜距历史近似写为：

$R_d\left(\mathrm{\η}\right)=\sqrt{R_{20}^2+{({\mathrm{\η}}_0-\frac{f_{\mathrm{dc}}}{f_{\mathrm{dr}}}-\mathrm{\η})}^2\·V^2}{-}^{}--\left(40\right)$

V为卫星的速度；

视在斜距历史表达式为：

$R_{\mathrm{deo}}\left(\mathrm{\η}\right)=R_d\left(\mathrm{\η}\right)+\frac{\mathrm{\ΔPRT}}{{\mathrm{PRT}}_R}\·c\·\mathrm{\η}-R_m---\left(25\right)$

设ΔPRT和R_m的估计量分别为ΔPRT和那么，按照估计值对视在斜距进行补偿后，补偿斜距历史为：

$R_{\mathrm{deo}}\left(\mathrm{\η}\right)=R_d\left(\mathrm{\η}\right)+\frac{\mathrm{\ΔPRT}-\mathrm{\ΔPRT}}{{\mathrm{PRT}}_R}\·c\·\mathrm{\η}-R_m+{\hat{R}}_m---\left(41\right)$

即 $R_{\mathrm{dec}}\left(\mathrm{\η}\right)=R_{\mathrm{deo}}\left(\mathrm{\η}\right)-\frac{\mathrm{\ΔPRT}}{{\mathrm{PRT}}_R}\·c\·\mathrm{\η}+{\hat{R}}_m---\left(42\right)$

设

$\mathrm{\Δv}=\frac{\mathrm{\ΔPRT}-\mathrm{\ΔPRT}}{{\mathrm{PRT}}_R}\·c---\left(43\right)$

$\mathrm{\ΔR}=-R_m+{\hat{R}}_m---\left(44\right)$

那么

R_dec(η)＝R_d(η)+Δv·η+ΔR    （45）

将R_d(η)展开为泰勒级数：

$R_d\left(\mathrm{\η}\right)=d_0+d_1\mathrm{\η}+d_2{\mathrm{\η}}^2+d_3{\mathrm{\η}}^3+\underset{m=4}{\overset{\∞}{\mathrm{\Σ}}}{d}_m{\mathrm{\η}}^m---\left(46\right)$

则（45）式写为：

$R_d\left(\mathrm{\η}\right)=d_0+d_1\mathrm{\η}+d_2{\mathrm{\η}}^2+d_3{\mathrm{\η}}^3+\underset{m=4}{\overset{\∞}{\mathrm{\Σ}}}{d}_m{\mathrm{\η}}^m+-\mathrm{\Δv}\·\mathrm{\η}+\mathrm{\ΔR}--\left(47\right)$

从（47）式可知，

的三次项系数为：

α₃＝2d₃·ΔR+2d₂·Δv+(2d₀d₃+2d₁d₂)    （48）

可见，点(ΔR,Δv,α₃)构成了一平面。

当ΔR＝Δv＝0时，补偿后的斜距历史即为真实的斜距历史，即：

$R_{\mathrm{dec}}\left(\mathrm{\η}\right)=R_d\sqrt{R_{20}^2+{({\mathrm{\η}}_0-\frac{f_{\mathrm{dc}}}{f_{\mathrm{dr}}}-\mathrm{\η})}^2\·V^2}{}^{}$

则 $R_{\mathrm{dec}}^2\left(\mathrm{\η}\right)=R_{20}^2+{({\mathrm{\η}}_0-\frac{f_{\mathrm{dc}}}{f_{\mathrm{dr}}}-\mathrm{\η})}^2\·V^2---\left(49\right)$

由（49）式可知三次项系数为0，因此(ΔR,Δv,α₃)构成不仅为一个平面，而且过原点，则有：

α₃＝2d₃·ΔR+2d₂·Δv    （50）

将（43）（44）代入（50），整理得：

${\mathrm{\α}}_3={2d}_3\·{\hat{R}}_m-{2d}_2\·\frac{\mathrm{\ΔPRT}}{{\mathrm{PRT}}_R}\·c+({2d}_2\·\frac{\mathrm{\ΔPRT}}{\mathrm{PRT}}\·c-{2d}_3{R}_m)---\left(51\right)$

即(α₃,

ΔPRT)构成平面，且该平面与α₃＝0的平面的交线为：

${\hat{R}}_m\frac{d_2}{d_3}+\frac{\mathrm{\ΔPRT}}{{\mathrm{PRT}}_R}\·c-(\frac{d_2}{d_3}\·\frac{\mathrm{\ΔPRT}}{\mathrm{PRT}}\·c-R_m)---\left(52\right)$

将（34）展开到无穷阶，有：

$\frac{\mathrm{\ΔPRT}}{{\mathrm{PRT}}_R}\·c=\mathrm{\λ}{\hat{f}}_{\mathrm{dc}}+\underset{m=0}{\overset{\∞}{\mathrm{\Σ}}}(m+1)a_m{\mathrm{\η}}_0^m---\left(53\right)$

将（53）代入（52）中得到：

${\hat{R}}_m=\frac{d_2}{d_3}\·\mathrm{\λ}{\hat{f}}_{\mathrm{dc}}+[\frac{d_2}{d_3}\·(\underset{m=0}{\overset{\∞}{\mathrm{\Σ}}}(m+1)a_m{\mathrm{\η}}_0^m-\frac{\mathrm{\ΔPRT}}{{\mathrm{PRT}}_R}\·c)+R_m]---\left(54\right){}_{}$

可见，

与

呈线性关系，因此，对每个遍历点的视在斜距历史Rdeo(η)按照公式（42）进行补偿后获得的ΔPRT与

的关系是意义的。

步骤八、根据步骤四获得的多普勒中心频率f_dc与划分误差ΔPRT的关系和步骤六获得的多普勒中心频率f_dc与起始采样距离R_m的关系，将补偿后的斜距写成：

R_dec(η;f_dc)＝R_deo(η)+k₀+k₁·f_dc    （55）

即将公式（34）和（39）代入到（42）中，k₀、k₁为整理后的常数项和一次项系数。

在[-1/PRT_R,1/PRT_R]范围内对f_dc遍历，并在每个遍历点按照公式（55）计算补偿后的斜距

并进行二次拟合，获得二次项系数b0、一次项系数b₁、常数项b₂；

将补偿后的斜距公式（49）展开获得二次项系数、一次项系数和常数项的表达式；

从而获得f_dc与R₂₀的关系为：

$R_{20}=\sqrt{b_0-{({\mathrm{\η}}_0-\frac{f_{\mathrm{dc}}}{f_{\mathrm{dr}}})}^2\·b_2}---\left(56\right)$

步骤九：将步骤四获得的

代入到公式（56）中的R20，获得fdc的估计值

步骤十，将步骤九获得的估计值

代入到公式（34）中，获得的多普勒中划分误差ΔPRT的估计值ΔPRT；将步骤九获得的估计值

代入到公式（37）中，获得固定频差Δf的估计值

步骤十、将步骤九获得的估计值代入到公式（34）中，获得的多普勒中划分误差ΔPRT的估计值ΔPRT；将步骤九获得的估计值

代入到公式（37）中，获得固定频差Δf的估计值

步骤十一：根据步骤九获得的划分误差估计值ΔPRT，固定频差估计值

和步骤二获得的时间间隔PRT_R，获得直达波信号或场景回波信号的同步误差补偿数据，具体为：

1）时间同步：对第n个脉冲，左移_n?ΔPRT

2）频率同步：对第n个脉冲，乘以相位

3）将步骤九中的估计值

代入公式（39）获得起始采样距离Rm的估计值

本发明公开了一种非合作式星地双基地SAR的时频同步方法，属于雷达信号处理领域。该方法的基础是接收天线接收场景回波的同时，来自卫星的直达波信号也会从天线旁瓣进入，由于直达波和场景回波的传播路径不同，所以它们可以从时间上区分，那么直达波信号就能被提取出来而不受场景回波信号的干扰。

但是由于天线旁瓣增益小，导致所获取的直达波信号的信噪比较弱，利用已有的方法可能导致较大的估计误差，使得同步处理失效。针对能够适应低信噪比条件的同步误差高精度估计方法进行了研究。建立了含同步误差情况下的信号模型，推导得到了关于时频同步误差未知量的估计方程组，其中利用短合成孔径时间内卫星轨道近似直线的特点来约束噪声的影响。

以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种非合作式星地双基地SAR时频同步方法，其特征在于，包括下列步骤： 

步骤一、将天线的主瓣指向场景，接收场景回波，旁瓣接收来自卫星的直达波信号，当卫星波束扫过场景时，接收机连续采集数据并存储； 

步骤二：对步骤一获得的数据进行脉冲压缩，获得相邻峰值之间的时间间隔PRT_R，按照PRT_R对步骤一获得的数据进行二维划分，获得以快时间τ和慢时间η表征的二维数据； 

步骤三、对步骤二获得的二维数据进行距离向脉压，获得直达波的波束中心穿越时刻η₀、视在斜距历史R_deo(η)和峰值位置处随慢时间变化的相位历史 

步骤四、利用星历数据计算得到的卫星与天线之间的最短斜距R₂₀的估计值 

步骤五：对步骤三获得的获取的视在斜距历史R_deo(η)进行三次拟合获得多普勒中心频率f_dc与划分误差ΔPRT的关系；所述划分误差ΔPRT为真实脉冲重复时间和估计的脉冲重复时间PRT_R之间的误差； 

步骤六、根据步骤三获得的峰值位置处随慢时间变化的相位历史

获取多普勒中心频率f_dc与收发载频的固定频差Δf的关系； 

步骤七、将ΔPRT在(-PRT_R/2,PRT_R/2)范围，

在(R₂₀-100km,R₂₀+100km)范围进行遍历，获得多普勒中心频率f_dc与起始采样距离R_m的关系；ΔPRT为划分误差ΔPRT的估计值，为起始采样距离R_m的估计值； 

步骤八、在[-1/PRT_R,1/PRT_R]范围内对f_dc遍历，得到f_dc和R₂₀的关系； 

步骤九、根据步骤四获得的和步骤八得到f_dc和R₂₀的关系，获得f_dc的估计值

步骤十、根据步骤九获得的估计值

和步骤五获得的多普勒中心频率f_dc与划分误差ΔPRT的关系，获得多普勒中划分误差ΔPRT的估计值ΔPRT；根据步骤九获得的估计值和步骤六获得的多普勒中心频率f_dc与收发载频的固定频差Δf的关系，获得固定频差Δf的估计值

步骤十一、根据步骤十获得的划分误差估计值ΔPRT，固定频差估计值

和步骤二获得的时间间隔PRT_R，获得直达波信号或场景回波信号的同步误差补偿数据，具体为： 

1）时间同步：对第n个脉冲，左移n·ΔPRT 

2）频率同步：对第n个脉冲，乘以相位

3）根据步骤九获得的估计值

和步骤七获得的多普勒中心频率f_dc与起始采样距离R_m的关系获得起始采样距离R_m的估计值

2.如权利要求1所述的一种非合作式星地双基地SAR时频同步方法，其特征在于，步骤五中获得多普勒中心频率f_dc与划分误差ΔPRT的关系具体方法为： 

对步骤三获得的视在斜距历史R_deo(η)进行三次拟合获得一次项、二次项和三次项系数分别记为a₁、a₂、a₃； 

将视在斜距距离向脉压后的表达式R_deo(η)在波束中心穿越时刻η₀进行泰勒展开，获得一次项、二次项和三次项系数的表达式； 

根据a₁、a₂、a₃与泰勒展开的一次项、二次项、三次项系数表达式获得多普勒中心频率f_dc与划分误差ΔPRT的关系以及多普勒调频率f_dr： 

c为光速，λ为信号载波波长，所述划分误差ΔPRT为真实脉冲重复时间和 估计的脉冲重复时间PRT_R之间的误差。 

3.如权利要求1所述的一种非合作式星地双基地SAR时频同步方法，其特征在于，步骤六中获取多普勒中心频率f_dc与收发载频的固定频差Δf的关系具体方法为： 

根据步骤三获得的峰值位置处随慢时间变化的相位历史

获取视在多普勒中心频率f_{dc_view}的值，进而获得收发载频的固定频差Δf与多普勒中心频率f_dc的关系： 

Δf＝f_{dc_view}-f_dc    （37）。 

4.如权利要求1所述的一种非合作式星地双基地SAR时频同步方法，其特征在于，步骤七获得多普勒中心频率f_dc与起始采样距离R_m的关系具体方法为： 

将ΔPRT在(-PRT_R/2,PRT_R/2)范围，

在(R₂₀-100km,R₂₀+100km)范围进行遍历，并在每个遍历点对视在斜距历史R_deo(η)按照公式（42）进行补偿，获得每个遍历点对应的补偿斜距历史R_dec(η)，并对获得的补偿斜距历史的平方计算三次项系数α₃；其中，ΔPRT为划分误差ΔPRT的估计值，

为起始采样距离R_m的估计值； 

取（α₃，ΔPRT，

）构成的平面与（0，ΔPRT，

）构成的平面的交线的方程，即为ΔPRT与

的关系： 

再根据步骤三获得的f_dc与ΔPRT的关系获得

与f_dc的关系为： 

其中，r₁、q₁为一次项系数，r₀、q₀为常数项。 

5.如权利要求1所述的一种非合作式星地双基地SAR时频同步方法，其 特征在于，步骤八中得到f_dc和R₂₀的关系具体方法为： 

根据步骤四获得的多普勒中心频率f_dc与划分误差ΔPRT的关系和步骤六获得的多普勒中心频率f_dc与起始采样距离

的关系，将补偿后的斜距表示为： 

R_dec(η;f_dc)＝R_deo(η)+k₀+k₁·f_dc    （55） 

k₀、k₁为常数项和一次系数项，为已知量； 

在[-1/PRT_R,1/PRT_R]范围内对f_dc遍历，并在每个遍历点按照公式（55）计算

并进行二次拟合，获得二次项系数b₀、一次项系数b₁、常数项b₂； 

将真实斜距历史表达式的平方展开获得二次项系数、一次项系数和常数项的表达式； 

根据b₀、b₁、b₂和真实斜距历史平方展开获得的二次项系数、一次项系数和常数项的表达式，获得f_dc与R₂₀的关系为： 

。

Images