CN103543452B - 一种基于多普勒频率展开的双基地合成孔径雷达成像方法 - Google Patents

Abstract

本发明公开了一种基于多普勒频率展开的双基地合成孔径雷达系统成像方法，它是通过将精确的双基SAR二维频谱表达式关于多普勒频率进行泰勒展开，简化了求解二维频谱表达式；利用omega-k算法的思想，根据二维频谱的表达式，得到系统的二维STOLT关系，最后通过利用二维STOLT关系，采用二维非均匀傅里叶变换来近似二维STOLT插值，完成系统的成像。与传统方法相比，由于它是基于简洁且精确的二维频谱，使得点目标聚焦良好，另外，充分考虑多普勒频率的高次项对成像结果的影响，使得成像算法能够对场景精确成像，并且由于只考虑多普勒频率，在很大程度上减少双基地成像算法的复杂度。

Description

一种基于多普勒频率展开的双基地合成孔径雷达成像方法

技术领域

本技术发明属于雷达技术领域，它涉及了双基地合成孔径(BiSAR)雷达成像技术领域。

背景技术

双基地合成孔径雷达(Bistatic synthetic aperture radar，简写为BiSAR)是指收发天线分置于两个不同运动平台的雷达系统。与单基地SAR相比，双基SAR具有构造特殊，安全性高，隐蔽性好，抗干扰能力强，低成本和灵活性强的优点，并且能够实现一些单基地SAR所无法实现的特殊模式，如前视成像。双基成像是一种非常具有应用价值的成像模式，可应用于导弹导航、恶劣天气下的飞机导航及着陆等方面。鉴于双基地SAR的多种优势，对双基地SAR成像技术的研究具有重要意义。

由于二维频谱的表达式是双基频域成像算法的研究的基础，但是双基地SAR斜距史表达式复杂，导致驻定相位时间的求解变得复杂，成为二维频谱的表达式求解的最严重的障碍。目前，已出现一些对双基地二维频谱求解的近似表达式，如Extending Lofeld's bistatic formula(ELBF)，Method of Series Reversion(MSR)和the Method of air-phase(AP)，详见R.Wang.et al，“Extending lofeld’sbistatic formula for the general bistatic sar confguration”，Yew Lam Neo.et al，“A Two-DimensionalSpectrum fr Bistatic SAR Processing Using Series Reversion”和Liu.z.et al，“Study on spaceborne/airborne hybrid bistatic sar image foration in frequency domain”。这些方法能够近似得到二维频谱表达式，但是表达式过于复杂，不利于后续成像处理的应用。

发明内容

本发明的目的是为了解决双基地斜距史的双根号问题及克服二维频谱求解复杂的缺点，提供一种基于多普勒频率展开的双基地合成孔径雷达系统成像方法。该方法一方面基于简洁且精确的二维频谱，使得点目标聚焦良好，另一方面充分考虑多普勒频率的高次项对成像结果的影响，使得成像算法能够对场景精确成像，并且由于只考虑多普勒频率，在很大程度上减少双基地成像算法的复杂度。

为了方便描述本发明的内容，首先作以下术语定义：

定义1、慢时间和快时间

慢时间就是方位向时间，记为t_a，指收发平台飞过一个飞行孔径所需要的时间，由于雷 达以一定的周期T_r发射接收脉冲，慢时间可以表示为一个离散化的时间变量t_a＝nT_r，其中n＝1，…N，n为脉冲重复周期序号，N为一个合成孔径内慢时间的离散个数，T_r为脉冲重复周期。

快时间就是距离向时间，指电磁波从平台传播到目标之间的时间，记为t。

方位向多普勒频率是指慢时间傅里叶变换之后的域。(详见“保铮，邢孟道，王彤，《雷达成像技术》，24-30页”)

距离向频率是指快时间傅里叶变换之后的域。(详见“保铮，邢孟道，王彤，《雷达成像技术》，24-30页”)

定义2、零多普勒时刻

零多普勒时刻是指多普勒频率为零时的慢时间，记为t_ac。

定义3、双基SAR系统相关参数描述

发射机平台斜距史记为R_s(t_a)， $R_s\left(t_a\right)=\sqrt{R_{s0}^2+V_s^2{(t_a-T_{s0})}^2}$

其中，t_a为慢时间，T_s0为发射平台相对于目标点的零多普勒时刻：

R_s0为发射平台在零多普勒时刻距目标的最近斜距：V_s是发射平台相对目标的运动速度大小，下标s表示发射机，x表示观测场景点的横轴坐标，y表示观测场景点的纵轴坐标，X_s0，Y_s0，H_s表示发射平台的初始三维坐标，所有坐标都是三维笛卡尔直角坐标系中的坐标；

接收机平台斜距史记为R_r(t_a)， $R_r\left(t_a\right)=\sqrt{R_{r0}^2+V_r^2{(t_a-T_{r0})}^2}$

其中，t_a为慢时间，T_r0为接收平台相对于目标点的零多普勒时刻：

R_r0为接收平台在零多普勒时刻距目标的最近斜距：V_r是接收平台相对目标的运动速度大小，下标r表示接收机，x表示观测场景点的横轴坐标，y表示观测场景点的纵轴坐标，X_r0，Y_r0，H_r表示接收平台的初始三维坐标，所有坐标都是三维笛卡尔直角坐标系 中的坐标；

双基SAR系统斜距史记为R(t_a)，R(t_a)＝R_s(t_a)+R_r(t_a) 

波数域频率： $k=2\mathrm{\π}\frac{f+F_0}{C}$

其中，f为距离向频率，F₀为发射信号中心频率，C为光速大小； 

发射机平台相位史记为θ_s(t_a)，θ_s(t_a)＝kR_s(t_a) 

接收机平台相位史记为θ_r(t_a)，θ_r(t_a)＝kR_r(t_a) 

其中，k为波数域频率，R_s(t_a)和R_r(t_a)分别为发射平台和接收机平台的斜距史；

双基SAR系统相位史记为，θ(t_a)＝θ_s(t_a)+θ_r(t_a)+2πf_at_a

其中，t_a为慢时间，f_a为方位向多普勒频率；

驻定相位时间点t_a1满足θ′(t_a1)＝0；

其中θ′(t_a1)表示双基SAR系统相位史θ(t_a)在时刻t_a1的一阶导数。

定义4、泰勒展开式 

泰勒展开式是一种数学上的函数等效方法。如果函数足够光滑的话，在已知函数在某一点各阶导数值的情况之下，可以用泰勒展开式的值做系数，构建一个多项式来近似已知函数在这一点邻域的值。泰勒展开式可以用若干项连加式来表示一个函数，这些相加的项由函数在某一点的导数求得。

定义5、双基SAR二维频谱模型 

双基SAR二维频谱： 

其中，表示二维频谱的相位项信息，

f为距离向频率，f_a为方位向多普勒频率，exp(·)表示以自然常数为底的指数函数，f_c表示载波频率，C表示光速的大小，j表示-1的平方根，R_M是等效单基的最短斜距史： T_M是等效单基零多普勒时刻：V_M是等效单基运动速度：a₃为误差修正项：其中，k_i为双基SAR斜距史R(t_a)关于方位向慢时间t_a在零多普勒时刻t_ac的第i阶泰勒展开式系数， 表示函数求导。

定义6、基于多普勒频率展开双基SAR二维频谱

基于多普勒频率展开双基SAR二维频谱： 

其中，定示二维频谱相位信息，exp(·)表示以自然常数为底的指数函数，关于多普勒频率f_a在多普勒中心频率f_ac处的泰勒展开式：

其中，M₁为点目标的距离向位置信息：M₂为点目标的方位向位置信息：M₃，M₄和M₅表示点目标位置信息的高次耦合项，其中 $M_3=K_2-k_1{K}_3-k_1^2{K}_4,$ M₄＝K₃-k₁K₄，M₅＝M₄，而 $K_l=\frac{{\∂}^{l}K}{\∂f_a^l}{|}_{f_a=f_{\mathrm{ac}}},l=\mathrm{0,1,2,3,4},$ K为二维频谱相位信息的一部分，k₁为双基SAR斜距史R(t_a)关于方位向慢时间t_a在零多普勒时刻t_ac的第1阶泰勒展开式系数，f_ac为多普勒中心频率：f为距离向频率，f_c表示载波频率，C表示光速的大小，R_M是等效单基的最短斜距史，T_M是等效单基零多普勒时刻，V_M是等效单基运动速度，a₃为误差修正项。j表示-1的平方根。

定义7、合成孔径雷达标准距离压缩方法

合成孔径雷达标准距离压缩方法是指利用合成孔径雷达发射信号参数，采用匹配滤波技术对合成孔径雷达的距离向信号进行滤波的过程。详见文献“雷达成像技术”，保铮等编著，电子工业出版社出版。

定义8、二维非均匀快速傅里叶变换

信号S(f，f_a)的二维非均匀快速傅里叶变化NUFFT(S(f，f_a))为

NUFFT(S(f，f_a))＝∫∫S(f，f_a)exp(-j2πft-j2πf_at_a)dfdf_a

其中，f和f_a分别表示距离向频率和方位向多普勒频率，是非均匀分布的，t和t_a分别表示快时间和慢时间，exp(·)表示以自然常数为底的指数函数，j表示-1的平方根。

定义9、双基雷达系统成像坐标系和观测场景坐标系

双基雷达成像坐标系：[ΔM₁，ΔM₁]

观测场景坐标系：[x，y]

其中，ΔM₁为非参考点与参考点的距离向位置信息之差：ΔM₁＝M_1no-M_1ref，ΔM₂为非参考点与参考点的方位向位置信息之差：ΔM₂＝M_2no-M_2ref，M_1no是非参考点距离向位置信息，M_1ref是参考点距离向位置信息，M_2no是非参考点方位向位置信息，M_2ref是参考点方位向位置信息，参考点是观测场景中的中心点，非参考点是观测场景中除了中心点以外的点，x表示观测场景点的横轴坐标，y表示观测场景点的纵轴坐标。

定义10、多元线性拟合

多元线性拟合的目的是将复合函数拟合为两个或者两个以上自变量的线性函数关系。在数值分析中，线性拟合是用解析表达式逼近复合函数。详见文献“数值计算引论”，白峰杉编著，高等教育出版社。

本发明提供一种基于多普勒频率展开的双基地合成孔径雷达成像方法，该方法的步骤如下：

步骤1、初始化双基地回波信号

双基地前视合成孔径雷达系统参数如下：在三维笛卡尔坐标系中，收发平台初始位置，分别记做P_s0(X_s0，Y_s0，H_s)和P_r0(X_r0，Y_r0，H_r)，其中X_s0，Y_s0，H_s分别表示发射平台的初始三维坐 标，X_r0，Y_r0，H_r，分别表接收平台的初始三维坐标；收发平台速度矢量分别为V_s(0，v_s，0)和V_r(0，v_r，0)，其中v_s表示发射平台y方向的速度大小，v_r表示接收平台y方向的速度大小；雷达发射线性调频信号，发射信号的中心频率为F₀，脉冲重复周期为PRF，发射脉冲的时宽为T，发射脉冲的调频斜率为K，发射脉冲的带宽为B，回波方位向采样点数Nplus，回波距离向上的采样点数N，其中Nplus和N均为正整数，距离向上的采样频率为F，观测场景的距离向总长度为R米，方位向总长度为Z米。

回波信号数据矩阵s(t，t_a)，它是一个Nplus行和N列的矩阵，回波数据矩阵s(t，t_a)的每行数据是快时间的回波信号采样数据，每列数据是慢时间的回波采样数据；以上双基地前视合成孔径雷达系统参数均为已知；

步骤2、回波信号距离向脉冲压缩

将步骤1中的回波信号数据矩阵s(t，t_a)在快时间上做快速傅里叶变换，并采用脉冲压缩处理方法进行标准脉冲压缩处理，得到距离向脉冲压缩后的距离频域信号SS(f，t_a)，其中，t为快时间，t_a为慢时间，f为距离向频率；

步骤3、回波信号的方位向傅里叶变换

对步骤2中得到距离向脉冲压缩后的距离频域信号SS(f，t_a)在慢时间上做快速傅里叶变换，得到回波信号的二维频谱S₁(f，f_a)，其中，t_a为慢时间，f_a为方位向多普勒频率，f为距离向频率；

步骤4、参考点相位补偿

采用多普勒频率展开方法，得到参考点二维频谱S₀(f，f_a)，再对参考点二维频谱S₀(f，f_a)取复共轭，得到将步骤3中得到的回波信号的二维频谱S₁(f，f_a)与参考点二维频谱S₀(f，f_a)的复共轭相乘，得到参考点相位补偿之后的回波信号的二维频谱S₂(f，f_a)，如下式所示

$S_2(f,f_a)=S_1(f,f_a)\×S_0^{*}(f,f_a)$

其中，参考点是观测场景中目标中心点，可以用下标ref表示。

根据计算公式： $M_{1\mathrm{ref}}-K_{0\mathrm{ref}}+k_{1\mathrm{ref}}^2{K}_{2\mathrm{ref}}+k_{1\mathrm{ref}}^3{K}_{3\mathrm{ref}}+k_{1\mathrm{ref}}^4{K}_{4\mathrm{ref}},$ $M_{2\mathrm{ref}}-K_{1\mathrm{ref}}-k_{1\mathrm{ref}}{K}_{2\mathrm{ref}}-k_{1\mathrm{ref}}^2{K}_{3\mathrm{ref}}-k_{1\mathrm{ref}}^3{K}_{4\mathrm{ref}},$ M_4ref＝K_3ref-k_1refK_4ref，M₅ref＝M_4ref，得到参考点的距离向位置信息M_1ref，参考点的方位向位置信息M_2ref以及参考点位置信息的高次耦合项M_3ref，M_4ref和M_5ref。其中， K_ref为参考点二维频谱相位信息的一部分，f_ec为多普勒中心频率，f_ac可以由下式计算：k_1ref为参考点斜距史R(t_a)关于方位向慢时间t_a在零多普勒时刻t_ac的第1阶泰勒展开式系数，R_Mref是参考点等效单基的最短斜距史，T_Mref是参考点等效单基零多普勒时刻，V_Mref是参考点等效单基运动速度，a_3ref为参考点误差修正项，R_Mref，T_Mref，V_Mref，a_3ref和k_iref可由下式计算： $R_{\mathrm{Mref}}=\sqrt{\frac{k_{\mathrm{oref}}^3}{9k_{2\mathrm{ref}}(9k_{2\mathrm{ref}}^3{k}_{\mathrm{oref}}+k_{3\mathrm{ref}}^2{k}_{\mathrm{oref}}^2)}},T_{\mathrm{Mref}}=\frac{3k_{3\mathrm{ref}}{k}_{2\mathrm{ref}}{k}_{\mathrm{oref}}^2}{9k_{2\mathrm{ref}}^3{k}_{\mathrm{oref}}+k_{3\mathrm{ref}}^2{k}_{\mathrm{oref}}^2},$ $V_{\mathrm{Mref}}=\sqrt{\frac{9k_{2\mathrm{ref}}^3{k}_{\mathrm{oref}}+k_{3\mathrm{ref}}^2{k}_{\mathrm{oref}}^2}{3k_{2\mathrm{ref}}}},a_{3\mathrm{ref}}=k_{1\mathrm{ref}}+\frac{k_{3\mathrm{ref}}{k}_{\mathrm{oref}}}{3k_{2\mathrm{ref}}}$ 和 $k_{\mathrm{iref}}=\frac{{\∂}^{i}R\left(t_a\right)}{\∂t_a^i}{|}_{t_a=t_{\mathrm{ac}}},i=\mathrm{0,1,2,3}$ 得到，其中，R(t_a)为参考点系统斜距史，t_ac为零多普勒时刻。

将参考点相位补偿之后的回波信号的二维频谱

S₂(f，f_a)表示为其中exp(·)表示以自然常数为底的指数函数，参考点相位补偿之后的回波信号的二维频谱S₂(f，f_a)的相位信息

f为距离向频率，f_a是方位向多普勒频率，f_c为载波频率，C表示光速大小，，j表示-1的平 方根，M_1ref为参考点的距离向位置信息，M_2ref代表参考点的方位向位置信息，M_3ref，M_4ref和M_5ref表示的是参考点位置信息的高次耦合项，M₁为点目标的距离向位置信息，M₂为点目标的方位向位置信息，M₃，M₄和M₅为点目标位置信息的高次耦合项，可以由多普勒频率展开方法得到。

步骤5、位置信息高次耦合项线性拟合

把步骤4中得到的参考点相位补偿之后的回波信号的二维频谱S₂(f，f_a)中位置信息高次项(M₃-M_3ref)，(M₄-M_4ref)和(M₅-M_5ref)按照多元线性拟合方法，线性拟合为(M₁-M_1ref)和(M₂-M_2ref)的复合函数，即  $\left\{\begin{array}{c}(M_3-M_{3\mathrm{ref}})\≈{\mathrm{\α}}_1(M_1-M_{1\mathrm{ref}})+{\mathrm{\β}}_1(M_2-M_{2\mathrm{rref}})\\ (M_4-M_{4\mathrm{ref}})\≈{\mathrm{\α}}_2(M_1-M_{1\mathrm{ref}})+{\mathrm{\β}}_2(M_2-M_{2\mathrm{rref}})\\ (M_5-M_{5\mathrm{ref}})\≈{\mathrm{\α}}_3(M_1-M_{1\mathrm{ref}})+{\mathrm{\β}}_4(M_2-M_{2\mathrm{rref}})\end{array}\right.,$ 其中，α_m，β_m分别为线性拟合系数， 则将参考点相位补偿之后回波信号的二维频谱S₂(f，f_a)中位置信息高次项(M₃-M_3ref)，(M₄-M_4ref)和(M₅-M_5ref)代入步骤4中得到的参考点相位补偿之后的回波信号二维频谱S₂(f，f_a)相位信息中，得到因为

而

利用 $\left\{\begin{array}{c}(M_3-M_{3\mathrm{ref}})\≈{\mathrm{\α}}_1(M_1-M_{1\mathrm{ref}})+{\mathrm{\β}}_1(M_2-M_{2\mathrm{rref}})\\ (M_4-M_{4\mathrm{ref}})\≈{\mathrm{\α}}_2(M_1-M_{1\mathrm{ref}})+{\mathrm{\β}}_2(M_2-M_{2\mathrm{rref}})\\ (M_5-M_{5\mathrm{ref}})\≈{\mathrm{\α}}_3(M_1-M_{1\mathrm{ref}})+{\mathrm{\β}}_4(M_2-M_{2\mathrm{rref}})\end{array}\right.,$

则

将此时的二维频谱记为S₃(f，f_a)，用如下表示

其中，

f为距离向频率，f_c为载波频率，C表示光速大小，f_a为方位向多普勒频率，M_1no是非参考点距离向位置信息，M_1ref是参考点距离向位置信息，M_2no是非参考点方位向位置信息，M_2ref是参考点方位向位置信息，参考点是观测场景中的目标中心点，非参考点是观测场景中除目标中心点的其他目标点，非参考点距离向位置信息M_1no和非参考点方位向位置信息M_2no可以由多普勒频率展开方法得到，下标no表示非参考点，exp(·)表示以自然常数为底的指数函数，j表示-1的平方根。

步骤6、多普勒频谱耦合项等效

对步骤5中得到二维频谱S₃(f，f_a)，将相位项中的频率耦合项等效到新的频率项，即  $\left\{\begin{array}{c}\left(\frac{f+f_C}{C}\right)-{\mathrm{\α}}_1\left(\frac{f+f_C}{C}\right)f_a^2-{\mathrm{\α}}_2{\left(\frac{f+f_C}{C}\right)}^2{f}_a^3-{\mathrm{\α}}_3{\left(\frac{f+f_C}{C}\right)}^3{f}_a^4=f^{\′}\\ f_a-{\mathrm{\β}}_1\left(\frac{f+f_C}{C}\right)f_a^2-{\mathrm{\β}}_2\frac{f+{f_C}^2}{C}{f}_a^3-{\mathrm{\β}}_3{\left(\frac{f+f_C}{C}\right)}^3{f}_a^4=f_a^{\′}\end{array}\right.,$ 其中，f为距离向频率，f_c为载波频率，C表示光速大小，f_a为方位向多普勒频率，f′为等效的距离向频率，f_a′为等效的方位向多普勒频率，α_m，β_m，m＝1，2，3表示位置信息耦合项的线性拟合系数；

等效频率的二维频谱S₄(f′，f_a′)，如下所示：

s₄(f′，f_a′)＝exp(-j2π{(M_1no-M_1ref)f′+(M_2no-M_2ref)f_a′})

其中，M_1no是非参考点距离向位置信息，M_1ref是参考点距离向位置信息，M_2no是非参考点方位向位置信息，M_2ref是参考点方位向位置信息，下标no表示飞参考点，下标ref表示参考点，exp(·)表示以自然常数为底的指数函数，j表示-1的平方根。

步骤7、二维非均匀快速傅里叶变换

对步骤6中得到的等效频率的二维频谱S₄(f′，f_a′)做非均匀傅里叶变换NUFFT(S₄(f′，f_a′))＝∫∫S₄(f′，f_a′)exp(-j2πf′t-j2πf_a′t_a)df′df_a′，实现将信号变换到斜距史图像域-方位向图像域，其中，f′为等效的距离向频率，f_a′为等效的方位向多普勒频率，t_a为慢时间，t为快时间。

经过上述步骤处理，就可以从双基地合成孔径雷达系统接收到的观测区域回波数据s(t，t_a)中获取具有较高分辨率的目标成像结果。

本发明的基本原理在于通过将精确的双基SAR二维频谱表达式关于多普勒频率进行泰勒展开，整理后得到简洁且意义更加明确的二维频谱表达式，并利用omega-k算法的思想，根据二维频谱的表达式，得到系统的二维STOLT关系，最后通过研究二维STOLT关系，利用二维非均匀傅里叶变换来近似二维STOLT插值，以完成对该系统的成像。

本发明的创新点是将双基SAR二维频谱关于多普勒频率展开，并利用omega-k算法的思想，根据二维频谱的表达式，得到系统的二维STOLT关系，最后通过研究二维STOLT关系，利用二维非均匀傅里叶变换来近似二维STOLT插值，以完成对该系统的成像。

本发明主要是针对精确的双基地SAR二维频谱表达式的复杂度，将二维频谱关于多普勒频率展开，从而简化了求解二维频谱表达式，并有效地利用了单基地成像算法的思想；通过二维频谱的解析表达式可以分析得到STOLT映射关系及其近似处理办法，避免了二维STOLT插值带来的运算量大和复杂度高的问题。

附图说明

图1为本发明具体实施方式采用的双基地合成孔径雷达飞行几何关系图

其中，P_t表达发射平台的飞行轨迹位置，P_r表达接收平台的飞行轨迹位置，表示发射平台的运动速度矢量，表示接收平台的运动速度矢量；0表示观测场景中的参考点；X，Y和Z表 示场景坐标轴。表示发射平台的雷达波束指向向量，表示接收平台的雷达波束指向向量，A₀表示观测场景中的参考点，A₁，A₂，A₃，A₄表示观测场景中的4个非参考点。

图2是本发明的工作流程框图

图3为双基地SAR回波数据经过本发明步骤一至步骤七处理后的结果图

其中，横轴表示快时间，纵轴表示慢时间，A₀表示观测场景中的参考点，A₁，A₂，A₃，A₄表示观测场景中的4个非参考点。

图4为本实验实施中采用的双基地SAR系统平台参数

具体实施方式

本发明主要采用仿真实验的方法进行验证该方案的可行性，所有步骤、结论都在MATLAB7.0上验证正确。

本例采用收发平台平行同向飞行，两平台的天线波束速度分别为2100(米每秒)和700(米每秒)。发射信号中心频率为9.65GHz，脉冲重复度为2500Hz，发射信号时宽为2(微秒)，系统平台其他参数如图4所示。观测区域的X轴总长度为1000米，Y轴总长度为200米，整个观测区域分布着5个反射点，相邻两个点的距离向差距为250米，方位向差距为100米。参考点A的坐标为(0,700)。

具体实施步骤如下：

步骤一、初始化双基地回波信号

使用图4所示的双基地SAR系统参数，仿真得到一个以1799行1024列数值矩阵存放的双基地SAR回波信号数据s(t，t_a)；每行数据存放的是对快时间回波信号的采样数据，每列数据存放的是对慢时间回波信号的采样数据；

步骤二、回波信号距离向压缩

将步骤一中的回波信号s(t，t_a)在快时间上做快速傅里叶变换，并进行标准脉冲压缩处理，得到距离向脉冲压缩后的距离频域信号SS(f，t_a)；

步骤三、回波信号的方位向傅里叶变换

对步骤二中得到距离向脉冲压缩后的距离频域信号SS(f，t_a)在慢时间上做快速傅里叶变换，则可以得到回波信号的二维频谱S₁(f，f_a)；

步骤四、参考点相位补偿

选择成像场景中心点为参考点，利用如下公式 $T_{\mathrm{Mref}}=\frac{3k_{3\mathrm{ref}}{k}_{2\mathrm{ref}}{k}_{\mathrm{oref}}^2}{9k_{2\mathrm{ref}}^3{k}_{\mathrm{oref}}+k_{3\mathrm{ref}}^2{k}_{\mathrm{oref}}^2},V_{\mathrm{Mref}}=\sqrt{\frac{9k_{2\mathrm{ref}}^3{k}_{\mathrm{oref}}+k_{3\mathrm{ref}}^2{k}_{\mathrm{oref}}^2}{3k_{2\mathrm{ref}}}},a_{3\mathrm{ref}}=k_{1\mathrm{ref}}+\frac{k_{3\mathrm{ref}}{k}_{\mathrm{oref}}}{3k_{2\mathrm{ref}}}$ 和 计算得到参考点等效的单基斜距史R_Mref＝8407.7(米)，参考点等效的单基零多普勒时间T_Mref＝20.0763(秒)，参考点等效的单基运动速度V_Mref＝239.5887(米每秒)和误差修正项系数a_3ref＝-12.09(米)，根据公式：

$M_{1\mathrm{ref}}-K_{0\mathrm{ref}}+k_{1\mathrm{ref}}^2{K}_{2\mathrm{ref}}+k_{1\mathrm{ref}}^3{K}_{3\mathrm{ref}}+k_{1\mathrm{ref}}^4{K}_{4\mathrm{ref}},$

$M_{2\mathrm{ref}}-K_{1\mathrm{ref}}-k_{1\mathrm{ref}}{K}_{2\mathrm{ref}}-k_{1\mathrm{ref}}^2{K}_{3\mathrm{ref}}-k_{1\mathrm{ref}}^3{K}_{4\mathrm{ref}},$

M_4ref＝K_3ref-k_1refK_4ref，M₅ref＝M_4ref，得到参考点的距离向位置信息M_1ref＝7.9527e3(秒)，参考点方位向位置信息M_2ref＝28.8495(秒)，位置信息的高次耦合项分别是M_3ref＝0.1668(秒)，M_4ref＝-5.87e-4(秒)和M_5ref＝1.7e-6(秒)，从而根据公式：得到参考点的基于多普勒频率展开的二维频谱S₀(f，f_a)，将其复共轭与回波信号的二维频谱S₁(f，f_a)相乘得到参考点相位补偿后的二维频谱S₂(f，f_a)；

步骤五、位置信息高次耦合项线性拟合

将参考点相位补偿后的二维频谱S₂(f，f_a)中的位置信息耦合项(M₃-M_3ref)，(M₄-M_4ref)和(M₅-M_5ref)按照线性拟合方法拟合为(M₁-M_1ref)和(M₂-M_2ref)的复合函数，根据公式 ${\mathrm{\α}}_m=\frac{d{M}_{m+2}}{d{M}_1}{{|}_{\begin{array}{c}x=0\\ y=0\end{array}},{\mathrm{\β}}_m=\frac{d{M}_{m+2}}{d{M}_2}{|}_{\begin{array}{c}x=0\\ y=0\end{array}},}_{}m=\mathrm{1,2,3},$ 得到各自的等效系数  $\begin{array}{c}{\mathrm{\α}}_1=-2.1e-5,{\mathrm{\β}}_1=0.008,{\mathrm{\α}}_2=2.2e-7,{\mathrm{\β}}_2=-4.46e-5,{\mathrm{\α}}_3=-3.4e-\\ 10,{\mathrm{\β}}_3=9.4e-8\end{array},$ 得到等效后 的二维频谱S₃(f，f_a)；

步骤六、多普勒频谱耦合项等效

对步骤五中得到等效后的二维频谱S₃(f，f_a)，将相位项中的频率耦合项等效为新的频率项，即 $\left\{\begin{array}{c}\left(\frac{f+f_C}{C}\right)-{\mathrm{\α}}_1\left(\frac{f+f_C}{C}\right)f_a^2-{\mathrm{\α}}_2{\left(\frac{f+f_C}{C}\right)}^2{f}_a^3-{\mathrm{\α}}_3{\left(\frac{f+f_C}{C}\right)}^3{f}_a^4=f^{\′}\\ f_a-{\mathrm{\β}}_1\left(\frac{f+f_C}{C}\right)f_a^2-{\mathrm{\β}}_2\frac{f+{f_C}^2}{C}{f}_a^3-{\mathrm{\β}}_3{\left(\frac{f+f_C}{C}\right)}^3{f}_a^4=f_a^{\′}\end{array}\right.,$ 因此，得到等效频率的二维频谱S₄(f′，f_a′)；

步骤七、二维非均匀快速傅里叶变换及成像结果显示

我们将步骤六中得到的等效频率的二维频谱S₄(f′，f_a′)利用非均匀傅里叶变换

NUFFT(S₄(f′，f_a′))＝∫∫S₄(f′，f_a′)exp(-j2πf′t-j2πf_a′t_a)df′df_a′，得到σ((M₁-M_1ref)，(M₂-M_2ref))，则实现将信号变换到斜距史图像域-方位向图像域，通过传统的Matlab绘图方法得到二维图形。本发明应用MATLAB的“contour”函数为获得等高线的二维成像结果。

经过上述步骤的处理，就可以从双基地SAR回波数据中获取既有较高分辨率的场景图像。

通过本发明具体实施方案的仿真及测试，本发明所提供的成像方案能够实现双基SAR成像，不仅简化双基地SAR二维频谱的求解表达式，而且充分利用了已成熟的单基地成像算法的思想，同时还获得了高质量的成像结果。

Claims (1)

1.一种基于多普勒频率展开的双基地合成孔径雷达成像方法，其特征是它包括如下步骤：

步骤1、初始化双基地回波信号

双基地前视合成孔径雷达系统参数如下：在三维笛卡尔坐标系中，收发平台初始位置，分别记做P_s0(X_s0，Y_s0，H_s)和P_r0(X_r0，Y_r0，H_r)，其中X_s0，Y_s0，H_s分别表示发射平台的初始三维坐标，X_r0，Y_r0，H_r分别表接收平台的初始三维坐标；收发平台速度矢量分别为V_s(0，v_s，0)和V_r(0，v_r，0)，其中v_s表示发射平台y方向的速度大小，v_r表示接收平台y方向的速度大小；雷达发射线性调频信号，发射信号的中心频率为F₀，脉冲重复周期为PRF，发射脉冲的时宽为T，发射脉冲的调频斜率为K，发射脉冲的带宽为B，回波方位向采样点数Nplus，回波距离向上的采样点数N，其中Nplus和N均为正整数，距离向上的采样频率为F，观测场景的距离向总长度为R米，方位向总长度为Z米；

回波信号数据矩阵s(t，t_a)，它是一个Nplus行和N列的矩阵，回波数据矩阵s(t，t_a)的每行数据是快时间的回波信号采样数据，每列数据是慢时间的回波采样数据；以上双基地前视合成孔径雷达系统参数均为已知；

步骤2、回波信号距离向脉冲压缩

将步骤1中的回波信号数据矩阵s(t，t_a)在快时间上做快速傅里叶变换，并采用脉冲压缩处理方法进行标准脉冲压缩处理，得到距离向脉冲压缩后的距离频域信号SS(f，t_a)，其中，t为快时间，t_a为慢时间，f为距离向频率；

步骤3、回波信号的方位向傅里叶变换

对步骤2中得到距离向脉冲压缩后的距离频域信号SS(f，t_a)在慢时间上做快速傅里叶变换，得到回波信号的二维频谱S₁(f，f_a)，其中，t_a为慢时间，f_a为方位向多普勒频率，f为距离向频率；

步骤4、参考点相位补偿

采用多普勒频率展开方法，得到参考点二维频谱S₀(f，f_a)，再对参考点二维频谱S₀(f，f_a)取复共轭，得到将步骤3中得到的回波信号的二维频谱S₁(f，f_a)与参考点二维频谱S₀(f，f_a)的复共轭相乘，得到参考点相位补偿之后的回波信号的二维频谱S₂(f，f_a)，如下式所示

其中，参考点是观测场景中目标中心点，可以用下标ref表示；

根据计算公式： 得到参考点的距离向位置信息M_1ref，参考点的方位向位置信息M_2ref以及参考点位置信息的高次耦合项M_3ref，M_4ref和M_5ref；其中， K_ref为参考点二维频谱相位信息的一部分，f_ac为多普勒中心频率，f_ac可以由下式计算：k_1ref为参考点斜距史R(t_a)关于方位向慢时间t_a在零多普勒时刻t_ac的第1阶泰勒展开式系数，R_Mref是参考点等效单基的最短斜距史，T_Mref是参考点等效单基零多普勒时刻，V_Mref是参考点等效单基运动速度，a_3ref为参考点误差修正项，R_Mref，T_Mref，V_mref，a_3ref和k_5ref可由下式计算： 和得到，其中，R(t_a)为参考点系统斜距史，t_ac为零多普勒时刻；

将参考点相位补偿之后的回波信号的二维频谱

S₂(f，f_a)表示为其中exp(·)表示以自然常数为底的指数函数，参考点相位补偿之后的回波信号的二维频谱S₂(f，f_a)的相位信息

f为距离向频率，f_a是方位向多普勒频率，f_c为载波频率，C表示光速大小，，j表示-1的平方根，M_1ref为参考点的距离向位置信息，M_2ref代表参考点的方位向位置信息，M_3ref，M_4ref和M_5ref表示的是参考点位置信息的高次耦合项，M₁为点目标的距离向位置信息，M₂为点目标的方位向位置信息，M₃，M₄和M₅为点目标位置信息的高次耦合项，可以由多普勒频率展开方法得到；

步骤5、位置信息高次耦合项线性拟合

把步骤4中得到的参考点相位补偿之后的回波信号的二维频谱S₂(f，f_a)中位置信息高次项(M₃-M_3ref)，(M₄-M_4ref)和(M₅-M_5ref)按照多元线性拟合方法，线性拟合为(M₁-M_1ref)和(M₂-M_2ref)的复合函数，即 其中，α_m，β_m分别为线性拟合系数， 则将参考点相位补偿之后回波信号的二维频谱S₂(f，f_a)中位置信息高次项(M₃-M_3ref)，(M₄-M_4ref)和(M₅-M_5ref)代入步骤4中得到的参考点相位补偿之后的回波信号二维频谱S₂(f，f_a)相位信息中，得到

因为

而

利用

则

将此时的二维频谱记为S₃(f，f_a)，用如下表示

其中，

f为距离向频率，f_e为载波频率，C表示光速大小，f_a为方位向多普勒频率，M_1no是非参考点距离向位置信息，M_1ref是参考点距离向位置信息，M_2no是非参考点方位向位置信息，M_2ref是参考点方位向位置信息，参考点是观测场景中的目标中心点，非参考点是观测场景中除目标中心点的其他目标点，非参考点距离向位置信息M_1no和非参考点方位向位置信息M_2no可以由多普勒频率展开方法得到，下标no表示非参考点，_exp(·)表示以自然常数为底的指数函数，j表示-1的平方根；

步骤6、多普勒频谱耦合项等效

对步骤5中得到二维频谱S₃(f，f_a)，将相位项中的频率耦合项等效到新的频率项，即 其中，f为距离向频率，f_c为载波频率，C表示光速大小，f_a为方位向多普勒频率，f′为等效的距离向频率，f′_a为等效的方位向多普勒频率， α_m，β_m，m＝1，2，3表示位置信息耦合项的线性拟合系数；

等效频率的二维频谱S₄(f′，f′_a)，如下所示：

S₄（f′，f′₂)＝exp(-j2π{(M_1no-M_1ref)f′+(M_2no-M_2ref)f′_a})

其中，M_1no是非参考点距离向位置信息，M_1ref是参考点距离向位置信息，M_2no是非参考点方位向位置信息，M_2ref是参考点方位向位置信息，下标no表示飞参考点，下标ref表示参考点，exp(·)表示以自然常数为底的指数函数，j表示-1的平方根；

步骤7、二维非均匀快速傅里叶变换

对步骤6中得到的等效频率的二维频谱S₄(f′，f′_a)做非均匀傅里叶变换NUFFT(S₄(f′，f′_a))＝∫∫S₄(f′，f′_a)exp(-j2πf′t-j2πf′_at_a)df′df′_a，实现将信号变换到斜距史图像域-方位向图像域，其中，f′为等效的距离向频率，f′_a为等效的方位向多普勒频率，t_a为慢时间，t为快时间；

经过上述步骤处理，从双基地合成孔径雷达系统接收到的观测区域回波数据s(t，t_a)中获取具有较高分辨率的目标成像结果。