CN103837874A - 用于地球同步轨道sar成像的二维非线性变调频方法 - Google Patents

Abstract

本发明公开了一种用于地球同步轨道SAR成像的二维非线性变调频方法。使用本发明能够有效地完成GEO SAR大场景回波的成像处理。本发明首先建立GEO SAR空变模型，在该模型中考虑了SAR成像处理中的聚焦参数即雷达的斜距历史4阶泰勒级数展开中的二阶系数在距离向、方位向上的二次变化，以及三阶系数在方位向上的线性变化，然后根据该空变模型进行距离向NCS成像处理和方位向NCS成像处理。与已有技术相比，本发明能够有效地完成地球同步轨道SAR大场景的回波成像处理，消除了参数两维空变性带来的方位向相位误差，能够完成精确的方位向聚焦，实现GEO SAR大场景（400km×400km）成像。

Description

用于地球同步轨道SAR成像的二维非线性变调频方法

技术领域

本发明涉及合成孔径雷达成像技术领域，具体涉及一种用于地球同步轨道SAR成像的二维非线性变调频方法。

背景技术

地球同步轨道合成孔径雷达（GEO SAR）是运行在距地面高度36000km的地球同步轨道卫星上的合成孔径雷达（SAR）系统。由于GEO SAR可以覆盖超过1000km以上的区域，同时对于重点区域的重访周期可以小于1天，GEO SAR在陆地探测、土壤水分检测、大气水汽测绘、地表形变检测等方面有巨大的实用价值。

由于GEO SAR轨道高度高，离心率不为零，卫星速度较小，合成孔径时间较长，观测范围较大，成像聚焦参数不仅具有距离向空变性，还具有方位向的空变性，因此需要有效的成像处理算法。在传统GEO SAR的成像算法中，非线性变调频（NCS）算法是比较精确的成像算法，但是由于忽略了GEO SAR成像参数的方位向空变性，成像场景和成像分辨率都受到了制约。

发明内容

有鉴于此，本发明提供了一种用于地球同步轨道SAR成像的二维非线性变调频方法，能够有效地完成GEO SAR大场景回波的成像处理。

本发明的用于地球同步轨道SAR成像的二维非线性变调频方法，包括如下步骤：

步骤1、建立GEO SAR空变模型：

SAR成像处理中聚焦参数k₁，k₂，k₃，k₄的距离向空变模型为：

k₁(R)＝k₁₀+k_r11·(R-R₀)

k₂(R)＝k₂₀+k_r21·(R-R₀)+k_r22·(R-R₀)²

k₃(R)＝k₃₀+k_r31·(R-R₀)

k₄(R)＝k₄₀+k_r41·(R-R₀)

其中，k₁(R)、k₂(R)、k₃(R)、k₄(R)分别为目标距雷达的距离为R时的、雷达斜距历史的4阶泰勒展开式

中的一阶、二阶、三阶、四阶系数；R表示孔径中心照射时刻目标与雷达的距离；R₀为孔径中心时刻场景中心点与雷达的距离；t_a为方位向时间；R(t_a)表示t_a时刻目标与雷达的距离，k₁₀、k₂₀、k₃₀、k₄₀为场景中心点的各阶参数，k_r11，k_r21，k_r22，k_r31，k_r41为各阶参数沿距离向变化的变化系数；

聚焦参数的方位向空变模型为：

$k_2(R,t_p)=k_{a20}\left(R\right)+k_{a21}\left(R\right)\·t_p+k_{a22}\left(R\right)\·t_p^2$

k₃(R,t_p)＝k_a30(R)+k_a31(R)·t_p

其中，t_p为孔径中心照射目标的方位时刻；k_a20(R)=k₂(R)，k_a30(R)=k₃(R)，k_a21(R)、k_a22(R)、k_a31(R)表示位于同一距离R处的不同方位向目标参数的变化系数；

步骤2，距离向NCS成像聚焦处理：

采用基于级数展开的距离向NCS算法进行距离向NCS成像聚焦处理，距离向聚焦之后的信号为：

$S_1(t_r,f_a)={\sin c}_r[t_r-\frac{2R}{c}]\·W_a\left(f_a\right)\exp \left[j2{\mathrm{\π\φ}}_{\mathrm{az}}(f_a,R)\right]\exp \left[j{\mathrm{\φ}}_{\mathrm{RP}}\left(R\right)\right]$

其中S₁(t_r，f_a)表示GEO SAR距离多普勒域信号，sinc_r(·)表示回波信号距离向聚焦后包络,t_r表示距离向时间，c表示光速，W_a(f_a)表示方位向频域包络，f_a表示方位向频率，φ_RP(R)为距离向NCS残留相位，φ_az(f_a,R)表示方位调制相位；

${\mathrm{\φ}}_{\mathrm{RP}}\left(R\right)=\frac{k_1^2\left(R\right)}{2\·\mathrm{\λ}\·k_2\left(R\right)}+\frac{k_1^3\left(R\right)\·k_3\left(R\right)}{4\·\mathrm{\λ}\·k_2^3\left(R\right)}+\frac{k_1^4\left(R\right)\·(9\·k_3^2\left(R\right)-4\·k_2\left(R\right)\·k_4\left(R\right))}{32\·\mathrm{\λ}\·k_2^5\left(R\right)}-\frac{2\·R}{\mathrm{\λ}}$

${\mathrm{\φ}}_{\mathrm{az}}(f_a;R)=\mathrm{\π}\·\left\{\begin{array}{c}(\frac{k_1\left(R\right)}{k_2\left(R\right)}+\frac{{3k}_3\left(R\right){k_1}^2\left(R\right)}{4{k_2}^3\left(R\right)}+\frac{{k_1}^3\left(R\right)({{9k}_3}^2\left(R\right)-{4k}_2\left(R\right)k_4\left(R\right))}{{{8k}_2}^5\left(R\right)})f_a\\ +(\frac{\mathrm{\λ}}{{4k}_2\left(R\right)}+\frac{{3k}_3\left(R\right)k_1\left(R\right)\mathrm{\λ}}{{{8k}_2}^3\left(R\right)}+\frac{{{3\mathrm{\λk}}_1}^2\left(R\right)({{9k}_3}^2\left(R\right)-{4k}_2\left(R\right)k_4\left(R\right))}{{{32k}_2}^5\left(R\right)}){f_a}^2\\ +(\frac{k_3\left(R\right){\mathrm{\λ}}^2}{{{16k}_2}^3\left(R\right)}+\frac{k_1\left(R\right){\mathrm{\λ}}^2({{9k}_3}^2\left(R\right)-{4k}_2\left(R\right)k_4\left(R\right))}{{{32k}_2}^5\left(R\right)}){f_a}^3\\ +\left(\frac{{\mathrm{\λ}}^3({{9k}_3}^2\left(R\right){\mathrm{\λ}}^3-4k_4\left(R\right)k_2\left(R\right))}{{{256k}_2}^5\left(R\right)}\right){f_a}^4\end{array}\right\}$

对残留相位φ_RP(R)进行补偿，补偿因子H₁为

H₁＝exp[-j2πφ_RP(R)]

补偿后的信号S₂(t_r，f_a)为

$S_2(t_r,f_a)=\sin c(t_r-\frac{2R}{c})W_a\left(f_a\right)\exp \left[j2{\mathrm{\π\φ}}_{\mathrm{az}}(f_a,R)\right]$

步骤3、方位向NCS成像处理

步骤3.1、补偿S₂(t_r，f_a)的φ_az(f_a;R)项中没有方位向空变性的相位项，补偿因子为H₂、H₃，同时加入一个非线性滤波补偿项H₄：

$H_2(f_a;R)=\begin{array}{c}(\frac{k_1\left(R\right)}{k_2\left(R\right)}+\frac{{3k}_3\left(R\right)k_1^2\left(R\right)}{{4k}_2^3\left(R\right)}+\frac{k_1^3\left(R\right)\·({9k}_3^2\left(R\right)-{4k}_2\left(R\right)k_4\left(R\right))}{{8k}_2^5\left(R\right)})f_a\\ +(\frac{{3k}_3\left(R\right)k_1\left(R\right)\mathrm{\λ}}{{{8k}_2}^3\left(R\right)}+\frac{{{3\mathrm{\λk}}_1}^2\left(R\right)({{9k}_3}^2\left(R\right)-{4k}_2\left(R\right)k_4\left(R\right))}{{32k}_2^5\left(R\right)}){f_a}^2\\ +\left(\frac{k_1\left(R\right){\mathrm{\λ}}^2({{9k}_3}^2\left(R\right)-{4k}_2\left(R\right)k_4\left(R\right))}{{32k}_2^5\left(R\right)}\right){f_a}^3\end{array}$

$H_3(f_a;R)=\exp [j\·\mathrm{\π}\·(-\frac{{\mathrm{\λ}}^2{k}_{a30}\left(R\right)f_a^3}{{32k}_{a20}^4\left(R\right)}-\frac{{\mathrm{\λ}}^3({9k}_{a30}^2\left(R\right)-{4k}_4\left(R\right)k_{a20}\left(R\right)){f_a}^4}{{256k}_{a20}^5\left(R\right)})]$

H₄(f_a)＝exp[jπ(p₃f_a ³+p₄f_a ⁴)]

其中，

$\left\{\begin{array}{c}{p}_3=\frac{f_{\mathrm{dr}1}(4\mathrm{\α}-1)}{{3f}_{\mathrm{dr}0}^3(2\mathrm{\α}-1)}\\ p_4=\frac{(4\mathrm{\α}-1)({3a}_{\mathrm{rt}}{f}_{\mathrm{dr}0}^4+{5f}_{\mathrm{dr}1}^2)-{4f}_{\mathrm{dr}0}{f}_{\mathrm{dr}2}\mathrm{\α}}{12(2\mathrm{\α}-1)f_{\mathrm{dr}0}^5}\end{array}\right.$

$f_{\mathrm{dr}0}\left(R\right)=-\frac{{4k}_{a20}\left(R\right)}{\mathrm{\λ}}$

$f_{\mathrm{dr}1}\left(R\right)=-\frac{{4k}_{a21}\left(R\right)}{\mathrm{\λ}}$

$f_{\mathrm{dr}2}\left(R\right)=-\frac{{4k}_{a22}\left(R\right)}{\mathrm{\λ}}$

$a_{\mathrm{rt}}=\frac{{\mathrm{\λ}}^2}{16}\·(-\frac{{3k}_{a30}\·k_{a21}}{k_{a20}^4}+\frac{k_{a31}}{k_{a20}^3})$

其中，α为方位向变标因子，为一常数，且α≠0.5；

补偿后的信号S₃(t_r，f_a)为：

S₃(t_r,f_a)＝S₂(t_r,f_a)·H₂(f_a;R)·H₃(f_a;R)·H₄(f_a)

步骤3.2、通过方位向逆傅里叶变换将信号S₃(t_r，f_a)转换到方位时域，方位时域信号为S₄(t_r，t_a)；在方位时域对信号S₄(t_r，t_a)进行方位向非线性变标，得到信号S₅(t_r，t_a)，其中变标因子H₅为

H₅(t_a)＝exp(jπ(q₂t_a ²+q₃t_a ³+q₄t_a ⁴))

其中，

$\left\{\begin{array}{c}{q}_2=(2\mathrm{\α}-1)f_{\mathrm{dr}0}\\ q_3=\frac{1}{3}(2\mathrm{\α}-1)f_{\mathrm{dr}1}\\ q_4=\frac{{14f}_{\mathrm{dr}1}^2\mathrm{\α}-{5f}_{\mathrm{dr}1}^2-{4f}_{\mathrm{dr}0}{f}_{\mathrm{dr}2}\mathrm{\α}+{6f}_{\mathrm{dr}0}^4{a}_{\mathrm{rt}}\mathrm{\α}-{3a}_{\mathrm{rt}}{f}_{\mathrm{dr}0}^4}{{12f}_{\mathrm{dr}0}}\end{array}\right.$

步骤3.3、利用傅里叶变化，将信号S₅(t_r，t_a)变换到方位频域，得到的信号S₆(t_r，f_a)：

$S_6(t_r,f_a)={\sin c}_r(t_r-\frac{2R}{c})W_a\left[\frac{f_a-q_2{t}_p}{(q_2+f_{\mathrm{dr}0})T_s}\right]\exp \left(\begin{array}{c}-j2\mathrm{\π}\frac{t_p}{2\mathrm{\α}}{f}_a+\mathrm{j\π}\left(\begin{array}{c}-\frac{{f_a}^2}{q_2+f_{\mathrm{dr}0}}\\ +\frac{(f_{\mathrm{dr}0}^3{p}_3+q_3){f_a}^3}{{(q_2+f_{\mathrm{dr}0})}^3}\\ +\frac{(f_{\mathrm{dr}0}^4{p}_4+q_4){f_a}^4}{{(q_2+f_{\mathrm{dr}0})}^4}\end{array}\right)\end{array}\right)$

其中，T_s为合成孔径时间；

对S₆(t_r，f_a)进行方位向压缩，方位向压缩函数H₆(f_a)为

$H_6\left(f_a\right)=\exp (-\mathrm{j\π}(-\frac{{f_a}^2}{q_2+f_{\mathrm{dr}0}}+\frac{({f_{\mathrm{dr}0}}^3{p}_3+q_3){f_a}^3}{{(q_2+f_{\mathrm{dr}0})}^3}+\frac{({f_{\mathrm{dr}0}}^4{p}_4+q_4){f_a}^4}{{(q_2+f_{\mathrm{dr}0})}^4}))$

步骤3.4、将步骤3.3的压缩结果利用方位向逆FFT变化转换到方位时域，得到最终的GEO SAR二维时域图像信号s(t_r，t_a)。

有益效果：

本发明方法，对比已有技术，能够有效地完成地球同步轨道SAR大场景的回波成像处理，具有以下效果：

（1）建立了精确的GEO SAR参数空变模型，保证了成像算法的精确性；

（2）消除了参数两维空变性带来的方位向相位误差，能够完成精确的方位向聚焦，实现GEO SAR大场景（400km×400km）成像。

附图说明

图1为本发明流程图。

图2为GEO SAR400km×400km成像结果图。

图3（a）为传统距离向NCS算法成像结果；图3（b）为本发明两维NCS算法成像结果。

具体实施方式

下面结合附图并举实施例，对本发明进行详细描述。

本发明提供了一种用于地球同步轨道SAR成像的二维非线性变调频方法，本发明首先建立GEO SAR空变模型，在该模型中考虑了SAR成像处理中的聚焦参数即雷达的斜距历史4阶泰勒级数展开中的二阶系数在距离向、方位向上的二次变化，以及三阶系数在方位向上的线性变化，然后根据该空变模型进行距离向NCS成像处理和方位向NCS成像处理，如图1所示，具体包括如下步骤：

步骤1、建立GEO SAR空变模型

SAR的工作模式是一种“Stop-and-go”的模式。设雷达的不同时刻下的位置坐标为（x_s(t_a),y_s(t_a),z_s(t_a)），其中t_a为方位向时间，目标位置坐标为（x,y,z），则目标相对于雷达的斜距历史可以表示为：

$R\left(t_a\right)=\sqrt{{(x_s\left(t_a\right)-x)}^2+{(y_s\left(t_a\right)-y)}^2+{(z_s\left(t_a\right)-z)}^2}---\left(1\right)$

其中，R(t_a)表示t_a时刻目标与雷达的距离。

为了精确描述目标的斜距历史，同时便于成像处理，GEO SAR成像处理一般将式（1）关于t_a做4阶泰勒级数展开，得到：

$R\left(t_a\right)=R+k_1{t}_a+k_2{t}_a^2+k_3{t}_a^3+k_4{t}_a^4---\left(2\right)$

其中，R表示孔径中心照射时刻（t_a=0）目标与雷达的距离，k₁，k₂，k₃，k₄分别表示目标斜距历史关于方位时间t_a的各次项系数，可以通过对目标的斜距历史（式（1））进行四阶多项式拟合得到。

该模型不仅能够精确描述GEO SAR各目标的斜距历程，同时能够描述GEOSAR远地点附近特殊的“近-远-近”斜距历程。

但是该模型是针对场景中心点建立的，随着成像场景的增大，场景中不同目标的k₁、k₂、k₃、k₄是不同的，即式（2）中各次项系数将具有空变性。由于各次项系数是SAR成像处理中的聚焦参数，这种空变的聚焦参数会引入空变的距离徙动校正误差和聚焦相位误差，影响最终的SAR图像聚焦。要获得精确的成像结果，成像处理中必须消除参数空变性。传统的距离向NCS算法只考虑了聚焦参数k₁、k₂、k₃、k₄在距离向上的线性变化，由于忽略了二阶系数k₂距离向、上的二次变化、二阶系数k₂在方位向上的二次变化以及三阶系数k₃在方位向上的线性变化，传统的距离向NCS算法成像场景和分辨率受到了很大限制。为了实现GEO SAR大场景高精度聚焦成像，本发明首先建立式（2）中各参数的两维的空变模型。

为了实现GEO SAR大场景精确聚焦成像，对现有GEO SAR参数距离向空变模型进行修正，将k2的距离向空变模型选定为二阶。故参数距离向空变模型为：

$\begin{array}{c}{k}_1\left(R\right)=k_{10}+k_{r11}\·\left(R{-R}_0\right)\\ k_2\left(R\right)=k_{20}+k_{r21}\·(R-R_0)+k_{r22}\·{(R-R_0)}^2\\ k_3\left(R\right)=k_{30}+k_{r31}\·(R-R_0)\\ k_4\left(R\right)=k_{40}+k_{r41}\·(R-R_0)\end{array}---\left(3\right)$

其中，k₁(R)，k₂(R)，k₃(R)，k₄(R)表示距离为R处的参数，R₀为孔径中心时刻场景中心点与雷达的距离。k₁₀，k₂₀，k₃₀，k₄₀为参考点（场景中心点）各阶参数，k_r11，k_r21，k_r22，k_r31，k_r41为各阶参数沿距离向变化的变化系数。

经过仿真实验，模型选定最多二阶即可满足精度要求，故参数方位向空变模型为：

$\begin{array}{c}{k}_2(R,t_p)=k_{a20}\left(R\right)+k_{a21}\left(R\right)\·t_p+k_{a22}\left(R\right)\·t_p^2\\ k_3(R,t_p)=k_{a30}\left(R\right)+k_{a31}\left(R\right)\·t_p\end{array}---\left(4\right)$

其中t_p为孔径中心照射目标的方位时刻，其中，k_a20(R)=k₂(R)，k_a30(R)=k₃(R)，k_a21(R)，k_a22(R)，k_a31(R)表示位于同一距离R处的不同方位向目标参数的变化系数，可以利用系统参数和成像几何关系拟合得到。

现有技术仅考虑了聚焦参数k₂在距离向上的一阶线性变化，当场景变大之后，不能实现精确成像，本发明为了实现大场景高精度的成像，考虑了二阶参数k₂的2阶变化。并且，现有技术没有考虑方位向参数空变性，本发明经过仿真试验，确定将方位向参数k₂的空变性模型建为如式（4）所示的二阶模型；k₃的空变模型如式（4）所示的一阶模型。

其中，距离向空变模型的各参数和方位向空变模型的各参数可以按照以下步骤计算得到：

（1）利用GEO SAR系统参数计算得到卫星在成像场景坐标系下不同时刻的三维位置坐标；

（2）选择场景坐标系原点作为参考点，该点三维位置坐标为（0,0,0）；

（3）在场景中沿距离向和方位向放置一个二维点阵，根据卫星场景坐标系下的位置坐标和各点目标的位置坐标计算得到各点目标的斜距历史；

（4）对各点目标的斜距历史进行四阶多项式拟合，可以得到各点目标的各阶参数k₁、k₂、k₃、k₄；其中参考点的各阶参数记做k₁₀、k₂₀、k₃₀、k₄₀；

（5）将各点目标的各阶参数分别沿距离向进行多项式拟合，便可以得到距离向空变模型式（3）中各阶参数沿距离向的变化系数k_r22、k_r11、k_r21、k_r22、k_r31、k_r41；

（6）将各点目标的各阶参数分别沿方位向进行多项式拟合，便可以得到方位向空变模型式（4）中各阶参数沿方位向的变化系数k_a21(R)、k_a22(R)、k_a31(R)。

然后分别进行距离向NCS成像处理和方位向NCS成像处理。

步骤2、距离向NCS成像处理

基于级数展开的距离向NCS算法可以精确的校正目标的距离徙动，实现大场景数据的距离向聚焦，因此本发明直接采用该方法实现GEO SAR回波数据的距离向聚焦。

距离向聚焦之后的信号在距离多普勒可以表示为：

$S_1(t_r,f_a)={\sin c}_r[t_r-\frac{2R}{c}]\·W_a\left(f_a\right)\exp \left[j2{\mathrm{\π\φ}}_{\mathrm{az}}(f_a,R)\right]\exp \left[j{\mathrm{\φ}}_{\mathrm{RP}}\left(R\right)\right]---\left(5\right)$

其中S₁(t_r，f_a)表示GEO SAR距离多普勒域信号，sinc_r(·)表示回波信号距离向聚焦后包络,t_r表示距离向时间，c表示光速，W_a(f_a)表示方位向频域包络，f_a表示方位向频率，φ_RP(R)为距离向NCS残留相位，φ_az(f_a,R)表示方位调制相位，可以表示为：

${\mathrm{\φ}}_{\mathrm{az}}(f_a;R)=\mathrm{\π}\·\left\{\begin{array}{c}(\frac{k_1\left(R\right)}{k_2\left(R\right)}+\frac{{3k}_3\left(R\right){k_1}^2\left(R\right)}{4{k_2}^3\left(R\right)}+\frac{{k_1}^3\left(R\right)({{9k}_3}^2\left(R\right)-{4k}_2\left(R\right)k_4\left(R\right))}{{{8k}_2}^5\left(R\right)})f_a\\ +(\frac{\mathrm{\λ}}{{4k}_2\left(R\right)}+\frac{{3k}_3\left(R\right)k_1\left(R\right)\mathrm{\λ}}{{{8k}_2}^3\left(R\right)}+\frac{{{3\mathrm{\λk}}_1}^2\left(R\right)({{9k}_3}^2\left(R\right)-{4k}_2\left(R\right)k_4\left(R\right))}{{{32k}_2}^5\left(R\right)}){f_a}^2\\ +(\frac{k_3\left(R\right){\mathrm{\λ}}^2}{{{16k}_2}^3\left(R\right)}+\frac{k_1\left(R\right){\mathrm{\λ}}^2({{9k}_3}^2\left(R\right)-{4k}_2\left(R\right)k_4\left(R\right))}{{{32k}_2}^5\left(R\right)}){f_a}^3\\ +\left(\frac{{\mathrm{\λ}}^3({{9k}_3}^2\left(R\right){\mathrm{\λ}}^3-4k_4\left(R\right)k_2\left(R\right))}{{{256k}_2}^5\left(R\right)}\right){f_a}^4\end{array}\right\}---\left(6\right)$

${\mathrm{\φ}}_{\mathrm{RP}}\left(R\right)=\frac{k_1^2\left(R\right)}{2\·\mathrm{\λ}\·k_2\left(R\right)}+\frac{k_1^3\left(R\right)\·k_3\left(R\right)}{4\·\mathrm{\λ}\·k_2^3\left(R\right)}+\frac{k_1^4\left(R\right)\·(9\·k_3^2\left(R\right)-4\·k_2\left(R\right)\·k_4\left(R\right))}{32\·\mathrm{\λ}\·k_2^5\left(R\right)}-\frac{2\·R}{\mathrm{\λ}}$

其中λ表示信号波长。

在进行方位向NCS成像处理之前，需要先补偿距离向NCS操作的残留相位φ_RP(R)，补偿因子H₁为：

H₁＝exp[-j2πφ_RP(R)]      （7）

补偿因子H₁与基于级数展开的距离向NCS算法中的φ_RP(R)的补偿因子一致。

补偿了残留相位φ_RP(R)之后的距离多普勒域信号S₂(t_r，f_a)为

$S_2(t_r,f_a)=\sin c(t_r-\frac{2R}{c})W_a\left(f_a\right)\exp \left[j2{\mathrm{\π\φ}}_{\mathrm{az}}(f_a,R)\right]---\left(8\right)$

步骤3、方位向NCS成像处理

首先，补偿φ_az(f_a;R)中没有方位向空变性的相位项，补偿因子为H₂、H₃，其中，

$H_2(f_a;R)=\begin{array}{c}(\frac{k_1\left(R\right)}{k_2\left(R\right)}+\frac{{3k}_3\left(R\right)k_1^2\left(R\right)}{{4k}_2^3\left(R\right)}+\frac{k_1^3\left(R\right)\·({9k}_3^2\left(R\right)-{4k}_2\left(R\right)k_4\left(R\right))}{{8k}_2^5\left(R\right)})f_a\\ +(\frac{{3k}_3\left(R\right)k_1\left(R\right)\mathrm{\λ}}{{{8k}_2}^3\left(R\right)}+\frac{{{3\mathrm{\λk}}_1}^2\left(R\right)({{9k}_3}^2\left(R\right)-{4k}_2\left(R\right)k_4\left(R\right))}{{32k}_2^5\left(R\right)}){f_a}^2\\ +\left(\frac{k_1\left(R\right){\mathrm{\λ}}^2({{9k}_3}^2\left(R\right)-{4k}_2\left(R\right)k_4\left(R\right))}{{32k}_2^5\left(R\right)}\right){f_a}^3\end{array}---\left(9\right)$

$H_3(f_a;R)=\exp [j\·\mathrm{\π}\·(-\frac{{\mathrm{\λ}}^2{k}_{a30}\left(R\right)f_a^3}{{32k}_{a20}^4\left(R\right)}-\frac{{\mathrm{\λ}}^3({9k}_{a30}^2\left(R\right)-{4k}_4\left(R\right)k_{a20}\left(R\right)){f_a}^4}{{256k}_{a20}^5\left(R\right)})]---\left(10\right)$

其中，由于φ_az(f_a;R)中的项 $(\frac{k_1\left(R\right)}{k_2\left(R\right)}+\frac{{3k}_3\left(R\right)k_1^2\left(R\right)}{{4k}_2^3\left(R\right)}+\frac{k_1^3\left(R\right)\·({9k}_3^2\left(R\right)-{4k}_2\left(R\right)k_4\left(R\right))}{{8k}_2^5\left(R\right)})f_a,$ $(\frac{{3k}_3\left(R\right)k_1\left(R\right)\mathrm{\λ}}{{{8k}_2}^3\left(R\right)}+\frac{{{3\mathrm{\λk}}_1}^2\left(R\right)({{9k}_3}^2\left(R\right)-{4k}_2\left(R\right)k_4\left(R\right))}{{32k}_2^5\left(R\right)}){f_a}^2,$ $\left(\frac{k_1\left(R\right){\mathrm{\λ}}^2({{9k}_3}^2\left(R\right)-{4k}_2\left(R\right)k_4\left(R\right))}{{32k}_2^5\left(R\right)}\right){f_a}^3$ 没有空变性，所以直接采用补偿因子H₂进行补偿；

由于φ_az(f_a;R)中

项包含空变性的项，也包含没有空变性的项，故将

进行泰勒展开后，采用补偿因子H₃中的对其没有空变性的项进行补偿；由于φ_az(f_a;R)中是4阶的相位项，可忽略其方位向空变性，直接采用补偿因子H3中的

进行补偿。

同时为了便于之后NCS操作，加入一个非线性滤波补偿项H₄。考虑到方位向空变模型中参数k₂为二阶，需构造一个四阶的滤波补偿项，k₃为一阶，需构造一个三阶的滤波补偿项，故非线性滤波补偿项H₄为

H₄(f_a)＝exp[jπ(p₃f_a ³+p₄f_a ⁴)]     （11）

其中，p₃、p₄因子为

$\left\{\begin{array}{c}{p}_3=\frac{f_{\mathrm{dr}1}(4\mathrm{\α}-1)}{{3f}_{\mathrm{dr}0}^3(2\mathrm{\α}-1)}\\ p_4=\frac{(4\mathrm{\α}-1)({3a}_{\mathrm{rt}}{f}_{\mathrm{dr}0}^4+{5f}_{\mathrm{dr}1}^2)-{4f}_{\mathrm{dr}0}{f}_{\mathrm{dr}2}\mathrm{\α}}{12(2\mathrm{\α}-1)f_{\mathrm{dr}0}^5}\end{array}\right.---\left(12\right)$

其中

$\begin{array}{c}{f}_{\mathrm{dr}0}\left(R\right)=-\frac{4k_{a20}\left(R\right)}{\mathrm{\λ}}\\ f_{\mathrm{dr}1}\left(R\right)=-\frac{{4k}_{a21}\left(R\right)}{\mathrm{\λ}}\\ f_{\mathrm{dr}2}\left(R\right)=-\frac{{4k}_{a22}\left(R\right)}{\mathrm{\λ}}\end{array}---\left(13\right)$

其中，α为一常数，且α≠0.5（p₄的分母不为0），定义为方位向变标因子。

补偿之后得到的信号为：

S₃(t_r,f_a)＝S₂(t_r,f_a)·H₂(f_a;R)·H₃(f_a;R)·H₄(f_a)       （14）

其中S₃(t_r，f_a)为方位向频域处理之后得到的信号，通过方位向逆傅里叶变换将信号S₃(t_r，f_a)转换到方位时域，方位时域信号为S₄(t_r，t_a)，在方位时域对信号S₄(t_r，t_a)进行方位向非线性变标操作，即乘以一个非线性的变标因子H₅，由于方位向空变模型中参数k₂为二阶，需构造一个四阶的变标因子，因此，H₅为

H₅(t_a)＝exp(jπ(q₂t_a ²+q₃t_a ³+q₄t_a ⁴))      （15）

其中，

$\left\{\begin{array}{c}{q}_2=(2\mathrm{\α}-1)f_{\mathrm{dr}0}\\ q_3=\frac{1}{3}(2\mathrm{\α}-1)f_{\mathrm{dr}1}\\ q_4=\frac{{14f}_{\mathrm{dr}1}^2\mathrm{\α}-{5f}_{\mathrm{dr}1}^2-{4f}_{\mathrm{dr}0}{f}_{\mathrm{dr}2}\mathrm{\α}+{6f}_{\mathrm{dr}0}^4{a}_{\mathrm{rt}}\mathrm{\α}-{3a}_{\mathrm{rt}}{f}_{\mathrm{dr}0}^4}{{12f}_{\mathrm{dr}0}}\end{array}\right.---\left(16\right)$

其中

$a_{\mathrm{rt}}=\frac{{\mathrm{\λ}}^2}{16}\·(-\frac{{3k}_{a30}\·k_{a21}}{k_{a20}^4}+\frac{k_{a31}}{k_{a20}^3})$

方位向信号非线性变标之后的信号为：

S₅(t_r,t_a)＝S₄(t_r,t_a)·H₅(t_a)（17）

利用傅里叶变化，将信号S₅(t_r，t_a)变换到方位频域，得到的信号为：

$S_6(t_r,f_a)={\sin c}_r(t_r-\frac{2R}{c})W_a\left[\frac{f_a-q_2{t}_p}{(q_2+f_{\mathrm{dr}0})T_s}\right]\exp \left(\begin{array}{c}-j2\mathrm{\π}\frac{t_p}{2\mathrm{\α}}{f}_a+\mathrm{j\π}\left(\begin{array}{c}-\frac{{f_a}^2}{q_2+f_{\mathrm{dr}0}}\\ +\frac{(f_{\mathrm{dr}0}^3{p}_3+q_3){f_a}^3}{{(q_2+f_{\mathrm{dr}0})}^3}\\ +\frac{(f_{\mathrm{dr}0}^4{p}_4+q_4){f_a}^4}{{(q_2+f_{\mathrm{dr}0})}^4}\end{array}\right)\end{array}\right)---\left(18\right)$

S₆(t_r，f_a)表示方位向NCS操作之后的方位向频域信号，其中，T_s为合成孔径时间。可以看出，此时信号中的各个参数均没有空变性，即聚焦参数的方位向空变性已经被消除，因此可以使用统一的聚焦参数进行方位向压缩，将S₆(t_r,f_a)中关于方位向频率f_a的2、3、4次项进行补偿，方位向压缩函数H₆(f_a)为：

$H_6\left(f_a\right)=\exp (-\mathrm{j\π}(-\frac{{f_a}^2}{q_2+f_{\mathrm{dr}0}}+\frac{({f_{\mathrm{dr}0}}^3{p}_3+q_3){f_a}^3}{{(q_2+f_{\mathrm{dr}0})}^3}+\frac{({f_{\mathrm{dr}0}}^4{p}_4+q_4){f_a}^4}{{(q_2+f_{\mathrm{dr}0})}^4}))---\left(19\right)$

式（18）与式（19）相乘，便完成了方位频域压缩，将压缩结果利用方位向逆FFT变化转换到方位时域，便可以得到最终聚焦的GEO SAR图像：

$s(t_r,t_a)={\sin c}_r(t_r-\frac{2R}{c})\·{\sin c}_a(t_a-\frac{t_p}{2\mathrm{\α}})---\left(20\right)$

其中，s(t_r，t_a)为最终聚焦的GEO SAR二维时域图像信号，sinc_a(·)表示图像方位向聚焦后包络。

自此，就实现了一种用于地球同步轨道SAR成像的两维NCS的方法。

为了验证本发明给出的两维NCS算法，进行了地球同步轨道SAR大场景（400km×400km）的回波仿真，信号波长为0.24m，信号带宽为8MHz。使用本发明方法处理的成像结果如图2所示，从图中可以看出，场景中各个目标都已经精确聚焦；同时，分别使用传统的距离向NCS算法和本发明方法中的两维NCS算法对场景边缘目标（200km,200km）进行成像处理，得到的处理结果如图3所示。由图可以明显观察出，传统的NCS算法的聚焦结果质量很差，方位向旁瓣展宽并且出现不对称旁瓣，而两维NCS算法聚焦结果良好。因此可以证明，两维NCS算法能够有效地对GEO SAR回波数据进行良好的成像处理。

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (1)

1.一种用于地球同步轨道SAR成像的二维非线性变调频方法，其特征在于，包括如下步骤：

步骤1、建立GEO SAR空变模型：

SAR成像处理中聚焦参数k₁，k₂，k₃，k₄的距离向空变模型为：

k₁(R)＝k₁₀+k_r11·(R-R₀)

k₂(R)＝k₂₀+k_r21·(R-R₀)+k_r22·(R-R₀)²

k₃(R)＝k₃₀+k_r31·(R-R₀)

k₄(R)＝k₄₀+k_r41·(R-R₀)

其中，k₁(R)、k₂(R)、k₃(R)、k₄(R)分别为目标距雷达的距离为R时的、雷达斜距历史的4阶泰勒展开式

中的一阶、二阶、三阶、四阶系数；R表示孔径中心照射时刻目标与雷达的距离；R₀为孔径中心时刻场景中心点与雷达的距离；t_a为方位向时间；R(t_a)表示t_a时刻目标与雷达的距离，k₁₀、k₂₀、k₃₀、k₄₀为场景中心点的各阶参数，k_r11，k_r21，k_r22，k_r31，k_r41为各阶参数沿距离向变化的变化系数；

聚焦参数的方位向空变模型为：

$k_2(R,t_p)=k_{a20}\left(R\right)+k_{a21}\left(R\right)\·t_p+k_{a22}\left(R\right)\·t_p^2$

k₃(R,t_p)＝k_a30(R)+k_a31(R)·t_p

其中，t_p为孔径中心照射目标的方位时刻；k_a20(R)=k₂(R)，k_a30(R)=k₃(R)，k_a21(R)、k_a22(R)、k_a31(R)表示位于同一距离R处的不同方位向目标参数的变化系数；

步骤2，距离向NCS成像聚焦处理：

采用基于级数展开的距离向NCS算法进行距离向NCS成像聚焦处理，距离向聚焦之后的信号为：

$S_1(t_r,f_a)={\sin c}_r[t_r-\frac{2R}{c}]\·W_a\left(f_a\right)\exp \left[j2{\mathrm{\π\φ}}_{\mathrm{az}}(f_a,R)\right]\exp \left[j{\mathrm{\φ}}_{\mathrm{RP}}\left(R\right)\right]$

其中S₁(t_r，f_a)表示GEO SAR距离多普勒域信号，sinc_r(·)表示回波信号距离向聚焦后包络,t_r表示距离向时间，c表示光速，W_a(f_a)表示方位向频域包络，f_a表示方位向频率，φ_RP(R)为距离向NCS残留相位，φ_az(f_a,R)表示方位调制相位；

${\mathrm{\φ}}_{\mathrm{RP}}\left(R\right)=\frac{k_1^2\left(R\right)}{2\·\mathrm{\λ}\·k_2\left(R\right)}+\frac{k_1^3\left(R\right)\·k_3\left(R\right)}{4\·\mathrm{\λ}\·k_2^3\left(R\right)}+\frac{k_1^4\left(R\right)\·(9\·k_3^2\left(R\right)-4\·k_2\left(R\right)\·k_4\left(R\right))}{32\·\mathrm{\λ}\·k_2^5\left(R\right)}-\frac{2\·R}{\mathrm{\λ}}$

${\mathrm{\φ}}_{\mathrm{az}}(f_a;R)=\mathrm{\π}\·\left\{\begin{array}{c}(\frac{k_1\left(R\right)}{k_2\left(R\right)}+\frac{{3k}_3\left(R\right){k_1}^2\left(R\right)}{4{k_2}^3\left(R\right)}+\frac{{k_1}^3\left(R\right)({{9k}_3}^2\left(R\right)-{4k}_2\left(R\right)k_4\left(R\right))}{{{8k}_2}^5\left(R\right)})f_a\\ +(\frac{\mathrm{\λ}}{{4k}_2\left(R\right)}+\frac{{3k}_3\left(R\right)k_1\left(R\right)\mathrm{\λ}}{{{8k}_2}^3\left(R\right)}+\frac{{{3\mathrm{\λk}}_1}^2\left(R\right)({{9k}_3}^2\left(R\right)-{4k}_2\left(R\right)k_4\left(R\right))}{{{32k}_2}^5\left(R\right)}){f_a}^2\\ +(\frac{k_3\left(R\right){\mathrm{\λ}}^2}{{{16k}_2}^3\left(R\right)}+\frac{k_1\left(R\right){\mathrm{\λ}}^2({{9k}_3}^2\left(R\right)-{4k}_2\left(R\right)k_4\left(R\right))}{{{32k}_2}^5\left(R\right)}){f_a}^3\\ +\left(\frac{{\mathrm{\λ}}^3({{9k}_3}^2\left(R\right){\mathrm{\λ}}^3-4k_4\left(R\right)k_2\left(R\right))}{{{256k}_2}^5\left(R\right)}\right){f_a}^4\end{array}\right\}$

对残留相位φ_RP(R)进行补偿，补偿因子H₁为

H₁＝exp[-j2πφ_RP(R)]

补偿后的信号S₂(t_r，f_a)为

$S_2(t_r,f_a)=\sin c(t_r-\frac{2R}{c})W_a\left(f_a\right)\exp \left[j2{\mathrm{\π\φ}}_{\mathrm{az}}(f_a,R)\right]$

步骤3、方位向NCS成像处理

步骤3.1、补偿S₂(t_r，f_a)的φ_az(f_a;R)项中没有方位向空变性的相位项，补偿因子为H₂、H₃，同时加入一个非线性滤波补偿项H₄：

$H_2(f_a;R)=\begin{array}{c}(\frac{k_1\left(R\right)}{k_2\left(R\right)}+\frac{{3k}_3\left(R\right)k_1^2\left(R\right)}{{4k}_2^3\left(R\right)}+\frac{k_1^3\left(R\right)\·({9k}_3^2\left(R\right)-{4k}_2\left(R\right)k_4\left(R\right))}{{8k}_2^5\left(R\right)})f_a\\ +(\frac{{3k}_3\left(R\right)k_1\left(R\right)\mathrm{\λ}}{{{8k}_2}^3\left(R\right)}+\frac{{{3\mathrm{\λk}}_1}^2\left(R\right)({{9k}_3}^2\left(R\right)-{4k}_2\left(R\right)k_4\left(R\right))}{{32k}_2^5\left(R\right)}){f_a}^2\\ +\left(\frac{k_1\left(R\right){\mathrm{\λ}}^2({{9k}_3}^2\left(R\right)-{4k}_2\left(R\right)k_4\left(R\right))}{{32k}_2^5\left(R\right)}\right){f_a}^3\end{array}$

$H_3(f_a;R)=\exp [j\·\mathrm{\π}\·(-\frac{{\mathrm{\λ}}^2{k}_{a30}\left(R\right)f_a^3}{{32k}_{a20}^4\left(R\right)}-\frac{{\mathrm{\λ}}^3({9k}_{a30}^2\left(R\right)-{4k}_4\left(R\right)k_{a20}\left(R\right)){f_a}^4}{{256k}_{a20}^5\left(R\right)})]$

H₄(f_a)=exp[jπ(p₃f_a ³+p₁f_a ⁴)]

其中，

$\left\{\begin{array}{c}{p}_3=\frac{f_{\mathrm{dr}1}(4\mathrm{\α}-1)}{{3f}_{\mathrm{dr}0}^3(2\mathrm{\α}-1)}\\ p_4=\frac{(4\mathrm{\α}-1)({3a}_{\mathrm{rt}}{f}_{\mathrm{dr}0}^4+{5f}_{\mathrm{dr}1}^2)-{4f}_{\mathrm{dr}0}{f}_{\mathrm{dr}2}\mathrm{\α}}{12(2\mathrm{\α}-1)f_{\mathrm{dr}0}^5}\end{array}\right.$

$f_{\mathrm{dr}0}\left(R\right)=-\frac{{4k}_{a20}\left(R\right)}{\mathrm{\λ}}$

$f_{\mathrm{dr}1}\left(R\right)=-\frac{{4k}_{a21}\left(R\right)}{\mathrm{\λ}}$

$f_{\mathrm{dr}2}\left(R\right)=-\frac{{4k}_{a22}\left(R\right)}{\mathrm{\λ}}$

$a_{\mathrm{rt}}=\frac{{\mathrm{\λ}}^2}{16}\·(-\frac{{3k}_{a30}\·k_{a21}}{k_{a20}^4}+\frac{k_{a31}}{k_{a20}^3})$

其中，α为方位向变标因子，为一常数，且α≠0.5；

补偿后的信号S₃(t_r，f_a)为：

S₃(t_r,f_a)＝S₂(t_r,f_a)·H₂(f_a;R)·H₃(f_a;R)·H₄(f_a)

步骤3.2、通过方位向逆傅里叶变换将信号S₃(t_r，f_a)转换到方位时域，方位时域信号为S₄(t_r，t_a)；在方位时域对信号S₄(t_r，t_a)进行方位向非线性变标，得到信号S₅(t_r，t_a)，其中变标因子H₅为

H₅(t_a)＝exp(jπ(q₂t_a ²+q₃t_a ³+q₄t_a ⁴))

其中，

$\left\{\begin{array}{c}{q}_2=(2\mathrm{\α}-1)f_{\mathrm{dr}0}\\ q_3=\frac{1}{3}(2\mathrm{\α}-1)f_{\mathrm{dr}1}\\ q_4=\frac{{14f}_{\mathrm{dr}1}^2\mathrm{\α}-{5f}_{\mathrm{dr}1}^2-{4f}_{\mathrm{dr}0}{f}_{\mathrm{dr}2}\mathrm{\α}+{6f}_{\mathrm{dr}0}^4{a}_{\mathrm{rt}}\mathrm{\α}-{3a}_{\mathrm{rt}}{f}_{\mathrm{dr}0}^4}{{12f}_{\mathrm{dr}0}}\end{array}\right.$

步骤3.3、利用傅里叶变化，将信号S₅(t_r，t_a)变换到方位频域，得到的信号S₆(t_r，f_a)：

$S_6(t_r,f_a)={\sin c}_r(t_r-\frac{2R}{c})W_a\left[\frac{f_a-q_2{t}_p}{(q_2+f_{\mathrm{dr}0})T_s}\right]\exp \left(\begin{array}{c}-j2\mathrm{\π}\frac{t_p}{2\mathrm{\α}}{f}_a+\mathrm{j\π}\left(\begin{array}{c}-\frac{{f_a}^2}{q_2+f_{\mathrm{dr}0}}\\ +\frac{(f_{\mathrm{dr}0}^3{p}_3+q_3){f_a}^3}{{(q_2+f_{\mathrm{dr}0})}^3}\\ +\frac{(f_{\mathrm{dr}0}^4{p}_4+q_4){f_a}^4}{{(q_2+f_{\mathrm{dr}0})}^4}\end{array}\right)\end{array}\right)$

其中，T_s为合成孔径时间；

对S₆(t_r，f_a)进行方位向压缩，方位向压缩函数H₆(f_a)为

$H_6\left(f_a\right)=\exp (-\mathrm{j\π}(-\frac{{f_a}^2}{q_2+f_{\mathrm{dr}0}}+\frac{({f_{\mathrm{dr}0}}^3{p}_3+q_3){f_a}^3}{{(q_2+f_{\mathrm{dr}0})}^3}+\frac{({f_{\mathrm{dr}0}}^4{p}_4+q_4){f_a}^4}{{(q_2+f_{\mathrm{dr}0})}^4}))$

步骤3.4、将步骤3.3的压缩结果利用方位向逆FFT变化转换到方位时域，得到最终的GEO SAR二维时域图像信号s(t_r，t_a)。

Images