[发明专利]基于星历数据的星载SAR高阶多普勒参数估算方法

权利要求书

1.基于星历数据的星载SAR高阶多普勒参数估算方法，包括如下步骤：

步骤1：根据回波数据存放格式，读取回波数据所对应的辅助数据，辅助数据包括：波前斜距R_min、信号采样率f_s、脉冲重复频率f_PRF、发射信号脉冲宽度τ、卫星的位置矢量R_os＝(R_o,sx,R_o,sy,R_o,sz)'、卫星的速度矢量V_os＝(V_sx,V_sy,V_sz)'、卫星的偏航角θ_y、卫星的俯仰角θ_p、卫星的横滚角θ_r、雷达天线工作视角θ_L；

步骤2：结合卫星的位置矢量R_os＝(R_o,sx,R_o,sy,R_o,sz)'、卫星的速度矢量V_os＝(V_sx,V_sy,V_sz)'、卫星的偏航角θ_y、卫星的俯仰角θ_p、卫星的横滚角θ_r，计算不转动的地心坐标系到卫星平台坐标系之间的转换矩阵A_or、卫星平台坐标系到卫星星体坐标系之间的转换矩阵A_re、卫星星体坐标系到天线坐标系之间的转换矩阵A_ea；

Aor=1000y1/r1x1/r10-x1/r1y1/r1vsx/v1-vsy/v10vsy/v1vsx/v10001v1/v0vsz/v010-vsz/v0v1/v---(14)]]>其中：r1=v1/v0-vsz/v010vsz/v0v1/vvsx/v1vsy/v10-vsy/v1vsx/v10001·(Ros×Vos)=(0,x1,y1);]]>

r₁＝|r₁|；

v＝|V_os|；

v1=vsx2+vsy2;]]>

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

Aea=1000cosθLsinθL0-sinθLcosθL]]>

步骤3：在依据步骤2计算各坐标系之间的相互转换矩阵后，进一步计算不转动地心坐标系中卫星与地面波束指向点之间的距离矢量，结合卫星的位置矢量和卫星与地面波束指向点之间的距离矢量，通过矢量相加处理得到地面波束指向点在不转动地心坐标系中的位置矢量；考虑到天线波束指向点在天线坐标系的Y轴上，满足如下向量[0,R_p,0]，通过矢量叠加处理获得地面波束指向点在不转动地球坐标系中的坐标位置：

Rop=Ros+Ro_sp]]>

=Ro,sxRo,syRo,sz+AorAreAea0Rp0]]>

=Ro,sxRo,syRo,sz+Aoe_11Aoe_12Aoe_13Aoe_12Aoe_22Aoe_23Aoe_13Aoe_23Aoe_330RpcosθL-RpsinθL---(15)]]>

=Ro,sxRo,syRo,sz+Aoe_12·RpcosθL-Aoe_13·RpsinθLAoe_22·RpcosθL-Aoe_23·RpsinθLAoe_32·RpcosθL-Aoe_33·RpsinθL=Ro,nxRo,nyRo,nz]]>

其中：R_{o_sp}表示不转动地心坐标系内卫星与地面波束指向点之间的距离矢量，N_r表示距离向采样点数；c表示光速；

步骤4：确定卫星的轨道根数，包涵轨道半长轴a，轨道倾角i，升交点赤经Ω，偏心率e，近地点幅角ω，真近心角θ，结合坐标转换矩阵之间的相互关系，有：

Aor=1000y1/r1x1/r10-x1/r1y1/r1vsx/v1-vsy/v10vsy/v1vsx/v10001v1/v0vsz/v010-vsz/v0v1/v=aor11aor21aor31aor12aor22aor32aor13aor23aor33---(16)]]>

=cosΩ-sinΩ0sinΩcosΩ00011000cosi-sini0sinicosicosω-sinω0sinωcosω0001-sin(θ-γ)-cos(θ-γ)0cos(θ-γ)-sin(θ-γ)0001]]>

其中：A_or中的各系数a_orij通过卫星位置矢量及速度矢量，i表明该系数为矩阵中的第i列系数，j表明该系数为矩阵中的第j行系数，1≤i≤3,1≤j≤3；γ表示航迹角，由下式计算获得：

tanγ=esinθ1+ecosθ---(17)]]>

结合上述转换矩阵之间的关系得到如下方程：

cosi=aor33sini·sinΩ=aor31---(18)]]>

通过求解方程（18），确定卫星轨道的轨道倾角i和升交点赤经Ω后，方程（16）修正为：

cosω-sinω0sinωcosω0001-sin(θ-γ)-cos(θ-γ)0cos(θ-γ)-sin(θ-γ)0001=-sin(ω+θ-γ)-cos(ω+θ-γ)0cos(ω+θ-γ)-sin(ω+θ-γ)0001]]>

=1000cosi-sini0sinicosi-1cosΩ-sinΩ0sinΩcosΩ0001-1aor11aor21aor31aor12aor22aor32aor13aor23aor33---(19)]]>

=bor11bor21bor31bor12bor22bor32bor13bor23bor33]]>

其中，上式中的右侧矩阵为已知矩阵，相应的计算得到：

cos(ω+θ-γ)=bor12sin(ω+θ-γ)=-bor11---(20)]]>

通过求解上式可计算得到角度(ω+θ-γ)；

结合卫星在不转动地球坐标系中的位置坐标R_os＝(R_o,sx,R_o,sy,R_o,sz)'，计算得到矢半径r为：

r=Ro,sx2+Ro,sy2+Ro,sz2---(21)]]>

同时位置坐标满足如下方程：

rcosΩ-sinΩ0sinΩcosΩ00011000cosi-sini0sinicosicosω-sinω0sinωcosω0001cos(θ)sin(θ)0=Ro,sxRo,syRo,sz---(22)]]>

将矢半径、轨道倾角及升交点赤经带入上述方程后有：

cosω-sinω0sinωcosω0001cos(θ)sin(θ)0=1000cosi-sini0sinicosi-1cosΩ-sinΩ0sinΩcosΩ0001-1xsyszs]]>   (23)

=cor1cor2cor3=cos(ω+θ)sin(ω+θ)0]]>

依据上式得到如下方程：

cos(ω+θ)=cor1sin(ω+θ)=cor2---(24)]]>通过求解上式可计算得到角度(ω+θ)；

综合式（20）和式（24），得到雷达系统的航迹角γ；

考虑到卫星在不转动地球坐标系中的速度矢量V_os＝(V_sx,V_sy,V_sz)'，满足如下方程关系：

VsxVsyVsz=μa(1-e2)Aov-sinθe+cosθ0---(25)]]>相应的得到如下方程：

Vsx2+Vsy2+Vsz2=μa(1-e2)[1+e2+2ecosθ]]]>

=μa(1-e2)[2(1+ecosθ)+e2-1]]]>           (26)

=2μ(1+ecosθ)a(1-e2)+μa(1-e2)[e2-1]]]>

=2μr-μa]]>

利用方程式（26），可以得到轨道的半长轴a；

结合开普勒雷达方程和航迹角计算方程，得到如下方程组（27）：

a(1-e2)=r(1+ecosθ)tanγ=esinθ1+ecosθ---(27)]]>对式（27）进行化简后有：

tan2γ·[1+2a(1-e2)-rr+(a(1-e2)-rr)2]=e2-(a(1-e2)-rr)2---(28)]]>通过求解上述四次方程组获得轨道的偏心率e，同时将轨道偏心率带入到航迹角的计算公式获得雷达卫星的瞬时真近心角θ，结合式（24）得到卫星轨道的近地点幅角ω；

步骤5：计算轨道平面坐标系与不转动地心坐标系之间的转换矩阵Α_ov，将轨道倾角、升交点赤经和近地点幅角带入到轨道平面坐标系与不转动地心坐标系之间的转换矩阵计算公式，计算得到轨道平面坐标系与不转动地心坐标系之间的转换矩阵A_ov；

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

步骤6：计算不转动地心坐标系中卫星的瞬时加速度矢量A_os及加速度变化矢量

Aos=μ(1+ecos(θ))2a2(1-e2)2Aov-cos(θ)-sin(θ)0---(29)]]>

A·os=μ(1+ecos(θ))3a3(1-e2)3μa(1-e2)Aovsin(θ)+3ecos(θ)sin(θ)e(2sin2(θ)-cos2(θ))-cos(θ)0---(30)]]>

步骤7：地算不转动地心坐标系中，地面波束指向点的速度矢量V_op、加速度矢量A_op及加速度变化矢量

Vop=Vo_pxVo_pyVo_pz=-ωeRo,nyωeRo,nx0---(31)]]>

Aop=Ao_pxAo_pyAo_pz=-ωe2Ro,nx-ωe2Ro,ny0---(32)]]>

A·op=A·o_pxA·o_pyA·o_pz=ωe3Ro,ny-ωe3Ro,nx0---(33)]]>

步骤8：计算相对位置矢量R_{o_sp}、相对速度矢量V_{o_sp}、相对加速度矢量A_{o_sp}和相对加速度变化率矢量

R_{o_sp}＝R_os+R_op                    （34）

V_{o_sp}＝V_os+V_op                       （35）

A_{o_sp}＝A_os+A_op                       （36）

A·o_sp=A·os+A·op---(37)]]>

步骤9：结合相对位置矢量、相对速度矢量、相对加速度矢量、相对加速度变化率矢量，估算回波信号多普勒参数；

f1=2λVo_spxRo_spx+Vo_spyRo_spy+Vo_spzRo_spzRo_spx2+Ro_spy2+Ro_spz2---(38)]]>

f2=2λ(Vo_spxVo_spx+Vo_spyVo_spy+Vo_spzVo_spzRo_spx2+Ro_spy2+Ro_spz2+Ao_spxRo_spx+Ao_spyRo_spy+Ao_spzRo_spzRo_spx2+Ro_spy2+Ro_spz2]]>    （39）

-(Vo_spxRo_spx+Vo_spyRo_spy+Vo_spzRo_spz)2(Ro_spx2+Ro_spy2+Ro_spz2)3/2)]]>

f3=2λ(A·o_spxRo_spx+A·o_spyRo_spy+A·o_spzRo_spz6Ro_spx2+Ro_spy2+Ro_spz2+Ao_spxVo_spx+Ao_spyVo_spy+Ao_spzVo_spz2Ro_spx2+Ro_spy2+Ro_spz2]]>

-(Ao_spxRo_spx+Ao_spyRo_spy+Ao_spzRo_spz)(Vo_spxRo_spx+Vo_spyRo_spy+Vo_spzRo_spz)2(Ro_spx2+Ro_spy2+Ro_spz2)3/2]]>       （40）

-(Vo_spxVo_spx+Vo_spyVo_spy+Vo_spzVo_spz)(Vo_spxRo_spx+Vo_spyRo_spy+Vo_spzRo_spz)2(Ro_spx2+Ro_spy2+Ro_spz2)3/2]]>

-(Vo_spxRo_spx+Vo_spyRo_spy+Vo_spzRo_spz)32(Ro_spx2+Ro_spy2+Ro_spz2)5/2)]]>

f4=2λ(A·o_spxVo_spx+A·o_spyVo_spy+A·o_spzVo_spz6Ro_spx2+Ro_spy2+Ro_spz2+Ao_spxAo_spx+Ao_spyAo_spy+Ao_spzAo_spz8Ro_spx2+Ro_spy2+Ro_spz2]]>   （41）

-(λf2)232(Ro_spx2+Ro_spy2+Ro_spz2)1/2-(λf1f3)224(Ro_spx2+Ro_spy2+Ro_spz2)1/2)]]>

综合式（38）～（41）获得回波信号的四阶多普勒参数，进而得到星载SAR高阶多普勒参数，应用于星载SAR成像处理。