[发明专利]一种基于非停走模型的星载SAR回波信号仿真方法

权利要求书

1.一种基于非停走模型的星载SAR回波信号仿真方法，首先，建立星载SAR的空间几何关系建模，S点表示卫星位置；E表示星下点，即地心O与卫星S连线和地球表面的交点；P点表示卫星波束中心同地球表面的交点；P_N表示北极点；C点表示升交点；K表示近地点；A点表示过P_N和P的大圆同赤道的交点；B点表示过P_N和E的大圆同赤道的交点；DE和D'P分别表示星下点轨迹和波束瞄准点轨迹；F点表示和点D、D'、P_N在同一个大圆上且在赤道上的点，其中，∠POA＝L为波束指向点的纬度；∠PSO＝β为天线波束下视角；为天线波束对应的地心角；∠COK＝ω为近地点幅角；∠FOD＝i为轨道倾角；∠EOC＝u为升交点幅角；

其特征在于，仿真方法具体包含如下几个步骤：

步骤一：确定仿真中心时刻；

获取星载SAR回波仿真参数，包括：卫星半长轴a、轨道倾角i、轨道偏心率e、近地点俯角ω、过近心点时刻τ、目标区域纬度L₁、天线波束下视角β、地球平均半径R_e、地球引力场引力常数μ；

利用卫星位置S、天线波束中心线同地球表面的交点P、地心O确定平面，利用式(1)确定天线波束对应的地心角

结合地心角轨道倾角i、目标区域纬度L₁，设P'_NO与DO垂直，且P'_NO与地球表面交于点P'_N，D、D'、P_N、P'_N在同一个地球大圆上，O'D'垂直P'_NO于O'，则利用式(2)和式(3)分别获取∠P_NO'P和∠D'O'P：

基于所获取的角∠D'O'P，获取相应的升交点幅角u、真近心角f及偏心角E：

f=u-ω        (5)

E=2·atan(1-e1+e·tan(f2))---(6)]]>

根据偏心角，获取回波信号仿真的仿真中心时刻t_M：

tM=τ+(E-e·sinE)a3μ---(7)]]>

步骤二：获取多普勒参数；

将仿真中心时刻带入到开普勒轨道方程，得到仿真中心时刻卫星在不转动地球坐标系中的位置矢量

θ=tMa3μ---(8)]]>

Aov=cosΩ-sinΩ0sinΩcosΩ00011000cosi-sini0sinicosicosω-sinω0sinωcosω0001---(9)]]>

R→s=xosyoszos=Aovrcosθrsinθ0---(10)]]>

其中，x_os、y_os和z_os为不转动地球坐标系中卫星位置矢量的坐标，A_ov为轨道平面坐标系到不转动地心坐标系的变换矩阵，Ω为升交点赤经，i为轨道倾角，ω为近地点幅角，θ为中心时刻真近心角，r为极矢径；

获取该时刻不转动地球坐标系中卫星的速度矢量和加速度矢量

V→s=VsxVsyVsz=μa(1-e2)Aov-sinθe+cosθ0---(11)]]>

A→s=AsxAsyAsz=-μ(1+ecosθ)2a2(1-e2)2Aovcosθsinθ0---(12)]]>

其中，V_sx、V_sy和V_sz为不转动地球坐标系中速度矢量的坐标，A_sx、A_sy和A_sz为不转动地球坐标系中卫星加速度矢量的坐标；

同时，结合仿真中心时刻卫星在不转动地球坐标系中的位置矢量天线波束下视角β以及仿真中心时刻t_M，确定不转动地球坐标系中天线波束中心线同地球表面交点的位置矢量考虑到该交点位于天线波束指向方向上，如果以R表示仿真中心时刻天线相位中心与天线波束中心线同地球表面的交点之间的距离，相应的不转动地球坐标系中天线波束中心线同地球表面交点的位置矢量满足如下方程：

R→T=xTyTzT=AovAvrAreAea0R0+xosyoszos+AovAvrArexeyeze---(13)]]>

Avr=-sin(θ-γ)-cos(θ-γ)0cos(θ-γ)-sin(θ-γ)0001---(14)]]>

γ=atanesinθ1+ecosθ---(15)]]>

Aea=1000cosθL-sinθL0sinθLcosθL---(16)]]>

Are=cosθy0-sinθy010sinθy0cosθycosθp-sinθp0sinθpcosθp00011000cosθr-sinθr0sinθrcosθr---(17)]]>

其中：x_T、y_T和z_T为不转动地球坐标系中天线波束中心线同地球表面交点位置矢量的坐标，A_ov为轨道平面坐标系到不转动地心坐标系的变换矩阵；A_vr为卫星平台坐标系到轨道平面坐标系的变换矩阵；A_re为卫星星体坐标系到卫星平台坐标系的变换矩阵；A_ea为天线坐标系到卫星星体坐标系的变换矩阵；γ为航迹角，x_e、y_e和z_e表示天线相位中心在卫星星体坐标系中的位置坐标；

将式(13)代入到地球椭球模型，通过求解方程得到斜距变量R，并将斜距变量R代入到式(13)，确定不转动地球坐标系中天线波束中心线同地球表面的交点的位置

xT2+yT2Ea2+zT2Eb2=1---(18)]]>其中：x_T、y_T和z_T为不转动地球坐标系中天线波束中心线同地球表面交点位置矢量的坐标，E_a表示地球的长半轴，E_b表示地球的短半轴；

根据不转动地球坐标系中天线波束中心线同地球表面交点的位置矢量得到不转动地球坐标系中交点处的速度矢量和加速度矢量

V→T=VTxVTyVTz=ωe×R→T=-ωeyTωexT0---(19)]]>

A→T=ATxATyATz=ωe2·R→T=ωe2xTωe2yTωe2zT---(20)]]>

其中：V_Tx、V_Ty和V_Tz为不转动地球坐标系中天线波束中心线同地球表面交点的速度矢量的坐标，A_Tx、A_Ty和A_Tz为不转动地球坐标系中天线波束中心线同地球表面交点的加速度矢量的坐标，ω_e表示地球自转的角速度；

计算不转动地球坐标系中雷达天线相位中心与天线波束中心线同地球表面交点之间的相对位置矢量相对速度矢量及相对加速度矢量

R→l=xT-xosyT-yoszT-zos-AovAvr_iArexeyeze---(21)]]>

V→l=VTx-VsxVTy-VsyVTz-Vsz---(22)]]>

A→l=ATx-AsxATy-AsyATz-Asz---(23)]]>

将相对位置矢量相对速度矢量及相对加速度矢量代入到多普勒参数计算公式，得到脉冲发射时刻对应仿真回波数据的多普勒中心频率f_{d_l}、多普勒调频率f_{r_l}及等效速度V_l和等效斜视角多普勒参数计算公式具体为：

R=(xT-xos)2+(yT-yos)2+(zT-zos)2---(24)]]>

fd_l=2V→l·R→lλR---(25)]]>

fr_l=2λ[V→l·V→lR+A→l·R→lR+(V→l·R→l)2R3]---(26)]]>

Vl=λRfr_l2+(λfd_l2)2---(27)]]>

其中：λ表示发射信号波长；

根据计算得到的斜距变量R、等效速度V_l和等效斜视角利用式(29)确定脉冲接收时刻t_{M_r}，并以此为基础结合式(8)～(26)，计算脉冲接收时刻对应仿真回波数据的多普勒中心频率f_{d_r}和多普勒调频率f_{r_r}；

其中，c为光速；

在获取发射脉冲仿真中心时刻对应仿真回波数据的多普勒中心频率f_{d_l}和多普勒调频率f_{r_l}及接收脉冲仿真中心时刻对应仿真回波数据的多普勒中心频率f_{d_r}和多普勒调频率f_{r_r}，取两者的平均值作为对应仿真回波数据的多普勒参数；

fd=fd_l+fd_r2---(30)]]>

fr=fr_l+fr_r2---(31)]]>

其中，f_d为平均多普勒中心频率，f_r为平均多普勒调频率；

步骤三：获取各脉冲发射时刻天线相位中心与各地面目标之间的距离矢量；

根据雷达系统的脉冲重复频率，得到仿真过程中雷达系统的各脉冲的发射时刻t_i，如式(32)所示；

ti=tM+(i-Np2)fPRF---(32)]]>

其中：i表示发射脉冲序列的序号；N_p表示回波数据仿真的脉冲数目；f_PRF表示雷达系统的脉冲重复频率；

将式(32)所得到的各脉冲发射时刻t_i代入到式(8)，得到各脉冲发射时刻所对应的真近心角θ_i，并将该真近心角θ_i代入到式(10)，得到各时刻卫星在不转动坐标系中的位置矢量获取雷达天线相位中心在不转动坐标系中的位置矢量速度矢量和加速度矢量

R→s_li=xos_iyos_izos_i+AovAvr_iArexeyeze---(33)]]>

V→s_li=μa(1-e2)Aov-sinθie+cosθi0---(34)]]>

A→s_li=μ(1+ecosθi)2a2(1-e2)2Aovcosθisinθi0---(35)]]>

其中，x_{os_i}、y_{os_i}和z_{os_i}为t_i时刻卫星在不转动坐标系中的位置矢量的坐标，x_e、y_e和z_e表示天线相位中心在卫星星体坐标系中的位置坐标，A_{vr_i}为t_i时刻卫星平台坐标系到轨道平面坐标系的变换矩阵，可通过将θ_i代入到式(14)得到t_i时刻的航迹角γ_i，再将θ_i和γ_i代入到式(13)得到；

步骤四：获取各脉冲发射时刻各地面目标对应回波数据的瞬时多普勒参数；

根据各脉冲发射时刻雷达天线相位中心在不转动坐标系中的位置矢量速度矢量和加速度矢量同时结合各地面目标在不转动坐标系中的位置矢量速度矢量和加速度矢量得到在不转动坐标系中雷达天线相位中心与各地面目标之间的相对位置矢量速度矢量和加速度矢量

R→stp_i=xstp_iystp_izstp_i=xtpytpztp-xos_iyos_izos_i-AovAvr_iArexeyeze---(36)]]>

V→stp_i=Vstpx_iVstpy_iVstpz_i=VtxVtyVtz-μa(1-e2)Aov-sinθie+cosθi0---(37)]]>

A→stp_i=Astpx_iAstpy_iAstpz_i=AtxAtyAtz+μ(1+ecosθi)2a2(1-e2)2Aovcosθisinθi0---(38)]]>

其中，x_{stp_i}、y_{stp_i}和z_{stp_i}为不转动坐标系中雷达天线相位中心与各地面目标之间的相对位置矢量的坐标，x_tp、y_tp和z_tp为不转动坐标系中地面目标位置矢量的坐标，V_{stpx_i}、V_{stpy_i}和V_{stpz_i}为不转动坐标系中雷达天线相位中心与各地面目标之间的相对速度矢量的坐标，V_tx、V_ty和V_tz为不转动坐标系中地面目标速度矢量的坐标，A_{stpx_i}、A_{stpy_i}和A_{stpz_i}为不转动坐标系中雷达天线相位中心与各地面目标之间的相对加速度矢量的坐标，A_tx、A_ty和A_tz为不转动坐标系中地面目标加速度矢量的坐标；

依据式(36)所得到的不转动坐标系中雷达天线相位中心与地面目标之间的位置矢量通过坐标转换处理得到天线坐标系内中雷达天线相位中心与地面目标之间位置矢量进而计算雷达天线相位中心与地面目标之间的距离R_{stp_i}及视线夹角θ_{stp_i}：

Rstp_i=(xstp_i_a2+ystp_i_a2+zstp_i_a2)1/2---(39)]]>

θstp_i=sin-1[xstp_i_a2+ystp_i_a2xstp_i_a2+ystp_i_a2+zstp_i_a2]1/2---(40)]]>

利用不转动坐标系中雷达天线相位中心与地面目标之间的相对位置矢量速度矢量和加速度矢量获取对应该目标的瞬时多普勒中心频率f_{dr_pi}、瞬时多普勒调频率f_{rr_pi}及天线相位中心与地面目标之间的等效速度V_{l_pi}和等效斜视角

fdr_pi=2V→stp_i·R→stp_iλRstp_i---(41)]]>

frr_pi=2λ[V→stp_i·V→stp_iRstp_i+A→stp_i·R→stp_iRstp_i+(V→stp_i·R→stp_i)2Rstp_i3]---(42)]]>

Vl_pi=λRfr_pi2+(λfd_pi2)2---(43)]]>

步骤五：确定各脉冲的接收时刻；

根据天线相位中心与地面目标之间的距离R_{stp_i}、天线相位中心与地面目标之间的等效速度V_{l_pi}和天线相位中心与地面目标之间的等效斜视角通过式(45)确定各地面目标对应仿真回波数据的脉冲接收时刻t_{i_r}；

步骤六：获取各地面目标接收时刻对应回波数据的瞬时多普勒参数；

将各地面目标对应仿真回波数据的脉冲接收时刻t_{i_r}代入到式(8)，得到该时刻雷达卫星对应的真近心角θ_{i_r}，并将该真近心角θ_{i_r}代入到式(10)，得到该时刻卫星在不转动坐标系中的位置矢量获取不转动坐标系中雷达天线相位中心的位置矢量速度矢量和加速度矢量

R→s_lir=xos_iryos_irzos_ir+AovAvr_irArexeyeze---(46)]]>

V→s_lir=μa(1-e2)Aov-sinθi_re+cosθi_r0---(47)]]>

A→s_lir=μ(1+ecosθi_r)2a2(1-e2)2Aovcosθi_rsinθi_r0---(48)]]>

其中，x_{os_ir}、y_{os_ir}和z_{os_ir}为t_{i_r}时刻卫星在不转动坐标系中的位置矢量的坐标，A_{vr_ir}为t_{i_r}时刻卫星平台坐标系到轨道平面坐标系的变换矩阵，可通过将θ_{i_r}代入到式(14)得到t_{i_r}时刻的航迹角γ_{i_r}，再将θ_{i_r}和γ_{i_r}代入到式(13)得到；

结合不转动坐标系中各地面目标的位置矢量速度矢量和加速度矢量得到该时刻不转动坐标系中雷达天线相位中心与各地面目标之间的相对位置矢量速度矢量和加速度矢量

R→stp_ir=xstp_irystp_irzstp_ir=xtpytpztp-xos_iryos_irzos_ir-AovAvr_irArexeyeze---(49)]]>

V→stp_ir=Vstpx_irVstpy_irVstpz_ir=VtxVtyVtz-μa(1-e2)Aov-sinθi_re+cosθi_r0---(50)]]>

A→stp_ir=Astpx_irAstpy_irAstpz_ir=AtxAtyAtz+μ(1+ecosθi_r)2a2(1-e2)2Aovcosθi_rsinθi_r0---(51)]]>

其中，x_{stp_ir}、y_{stp_ir}和z_{stp_ir}为t_{i_r}时刻不转动坐标系中雷达天线相位中心与地面目标之间的相对位置矢量的坐标，V_{stpx_ir}、V_{stpy_ir}和V_{stpz_ir}为t_{i_r}时刻不转动坐标系中雷达天线相位中心与地面目标之间的相对速度矢量的坐标，A_{stpx_ir}、A_{stpy_ir}和A_{stpy_ir}为t_{i_r}时刻不转动坐标系中雷达天线相位中心与地面目标之间的相对加速度矢量的坐标；

依据式(49)所得到的不转动坐标系中雷达天线相位中心与地面目标之间的位置矢量通过坐标转换处理得到天线坐标系内中雷达天线相位中心与地面目标之间位置矢量进而计算雷达天线相位中心与地面目标之间的距离R_{stp_ir}及视线夹角θ_{stp_ir}：

Rstp_ir=(xstp_ir_a2+ystp_ir_a2+zstp_ir_a2)1/2---(52)]]>

θstp_ir=sin-1[xstp_ir_a2+ystp_ir_a2xstp_ir_a2+ystp_ir_a2+zstp_ir_a2]1/2---(53)]]>

利用不转动坐标系中雷达天线相位中心与地面目标之间的位置矢量速度矢量和加速度矢量获取对应该目标的瞬时多普勒中心频率f_{dr_pir}、瞬时多普勒调频率f_{rr_pir}及天线相位中心与地面目标之间的等效速度V_{l_pir}和等效斜视角

fdr_pir=2V→stp_ir·R→stp_irλRstp_ir---(54)]]>

frr_pir=2λ[V→stp_ir·V→stp_irRstp_ir+A→stp_ir·R→stp_irRstp_ir+(V→stp_ir·R→stp_ir)2Rstp_ir3]---(55)]]>

Vl_pir=λRfr_pir2+(λfd_pir2)2---(56)]]>

步骤七：获取仿真回波数据；

根据式(39)、式(40)、式(52)和式(53)所得到的各地面目标发射时刻天线相位中心与地面目标之间的距离R_{stp_i}及视线夹角θ_{stp_i}和接收时刻天线相位中心与地面目标之间的距离R_{stp_ir}及视线夹角θ_{stp_ir}，获取每个距离采样时刻的瞬时收发双程斜距R_p(t_i,τ)，以及发射时刻瞬时雷达天线方位向增益函数W_{a_pl}(t_i,τ)和接收时刻瞬时雷达天线方位向增益函数W_{a_pr}(t_i,τ)：

Rp(ti,τ)=Rstp_i+Rstp_ir2+λ(fd_pi+fd_pir)4τ+λ(fr_pi+fr_pir)8τ2---(58)]]>

Wa_pl(ti,τ)=sin(πθstp_iLac(fc+fτ))(πθstp_iLac(fc+fτ))=sin(πθstp_iLac(fc+bτ))(πθstp_iLac(fc+bτ))---(59)]]>

Wa_pr(ti,τ)=sin(πθstp_irLac(fc+fτ))(πθstp_irLac(fc+fτ))=sin(πθstp_irLac(fc+bτ))(πθstp_irLac(fc+bτ))---(60)]]>

其中：L_a表示方位向天线尺寸;f_c表示发射信号载频;f_τ表示距离向瞬时频率;b表示发射线性调频信号的调频率；τ表示回波数据距离向时刻；

在得到每个距离采样时刻的瞬时收发双程斜距R_p(t_i,τ)，以及发射时刻瞬时雷达天线方位向增益函数W_{a_pl}(t_i,τ)和接收时刻瞬时雷达天线方位向增益函数W_{a_pr}(t_i,τ)的基础上，将上述参数带入到回波信号的数学模型，得到最终的回波信号；

Secho(ti,τ)=σWa_pl(ti,τ)Wa_pr(ti,τ)·exp{-j4πλRp(ti,τ)}·exp{-jπb(τ+NTp-2Rp(ti,τ)c)2}---(61)]]>

当逐脉冲完成各个脉冲仿真处理后，得到最终的回波信号数据。