CN102866393A - 一种基于pos与dem数据的sar多普勒参数估计方法 - Google Patents

Abstract

本发明涉及一种基于POS与DEM数据的SAR多普勒参数估计方法，该方法充分利用机载GPS导航数据、IMU姿态根据和数字高程数据，根据距离方程、共面方程和考虑地面高程数据的地球椭球数据来求解反射信号的目标位置，根据传感器状态和目标位置信息估计波束中心的多普勒信息。本发明充分能够对非正侧视、非平稳飞行以及地面地形复杂条件下机载SAR波束中心多普勒参数进行高精度估计，在SAR信号处理和成像领域有着重要的作用。

Description

一种基于POS与DEM数据的SAR多普勒参数估计方法

技术领域

本发明属于雷达数据处理领域，特别涉及到雷达波束中心的多普勒参数估计。

背景技术

由于飞机在飞行的过程中受到大气和其它干扰因素，复杂地形环境下非匀速直线飞行对多普勒参数的估计精度有很大的影响，而SAR多普勒参数估计精度会影响到SAR成像聚焦，导致噪声和模糊度的增高，严重时还会影响到图像的质量。目前SAR多普勒中心频率估计方法虽然有很多种，但大多基于一个相同的基本原理：利用SAR***的大时宽带宽积特性，多普勒频率和方位时间有着确定的对应关系，回波某个多普勒频率上的能量必然来自于雷达波束中某个特定方向上的目标，每个点目标回波在方位上都是一个线性调频信号，其中心频率f_DC被天线方向图调制。SAR目标由大量散射点组成，如果他们有相同的散射截面，那么回波的方位向功率谱密度就与天线功率方向图有相同的形状。通常天线方向图都是相对于波束中心对称的，所以回波方位向功率谱密度的峰值即为f_DC。

在整个成像过程中有三个地方的数据可以用来估计f_DC。一种早期的杂波锁定方法是在距离压缩之前完成；另一种方法是使用距离压缩之后、方位压缩之前的数据，把数据方位频谱的峰值位置作为多普勒中心频率的估计值；还可以利用功率谱与自相关函数，是一对傅里叶变换，从时域来估计多普勒中心频率，即自相关函数法。估计多普勒调频斜率的方法一直是研究的热点，目前大多数SAR处理器都使用图像偏移法来估计多普勒调频斜率，其估计精度也是与目标特性有密切关系。

SAR遥感影像的测量主要以几何构像模型为基础展开的。目前，雷达影像严密构像方程主要有基于传感器状态矢量和多普勒参数的模型，以及基于传感器状态矢量和姿态参数(外方位元素)的模型。Leberl F与R-D模型是侧视雷达影像严密定位使用得最广泛的几何模型；以外方位元素作为定向参数的模型除雷达共线方程外，距离-共面(R-Cp)方程是一种简洁严密的几何构像模型，其根据距离条件和波束中心共面条件来构建。距离-共面条件是指一行影像对应的所有地面点均在该行影像摄影时刻天线发射的雷达波束中心面内，这个波束中心面由传感器状态矢量和姿态确定；传感器至地面目标点的空间距离与雷达波测量的距离相等。距离-共面方程根据如下公式建立：

$\left\{\begin{array}{c}R=|\stackrel{\→}{\mathrm{OP}}-\stackrel{\→}{\mathrm{OS}}|\\ \stackrel{\→}{i}\·(\stackrel{\→}{\mathrm{OP}}-\stackrel{\→}{\mathrm{OS}})=0\end{array}\right.---\left(2\right)$

分别为波束中心面的法向量、地面点位置向量和传感器位置向量。R为雷达波测量斜距。

姿态转角***的距离-共面方程的展开式为：

上式中，M_r为距离向分辨率，R₀为最近端斜距，y为像元在影像上的列坐标。

距离多普勒模型不仅用于SAR信号的成像，也用于SAR影像的对地定位。距离多普勒模型根据传感器与地面点间的距离条件和被观测目标的多普勒频率条件构建，即点目标的空间定位由两方面来决定:一是根据目标回波时间的长短来确定点目标到雷达的距离,二是根据目标回波的多普勒特性。根据距离向上回波信号的时延和方位向上的多普勒频移这两个信息来确定点目标位置的方法，即距离-多普勒(R-D)定位方法。R-D方程表达式如下：

$\left\{\begin{array}{c}\sqrt{{(X-\mathrm{Xs})}^2+{(Y-\mathrm{Ys})}^2+{(Z-\mathrm{Zs})}^2}-R_j=0\\ 2[(X-\mathrm{Xs})V_X+(Y-\mathrm{Ys})V_Y+(Z-\mathrm{Zs})V_Z]/\mathrm{\λ}{R}_j-f_{\mathrm{DC}}=0\end{array}\right.$

其中(X，Y，Z)为地面点坐标、(Xs，Ys，Zs，V_S，V_Y，V_Z)为SAR天线中心状态矢量，λ、f_DC、R_j为SAR波长、多普勒频率及雷达波实测距离。

尽管对SAR处理进行了相当的改进，但是许多实用化SAR处理***仍饱受大量景成像中不可靠多普勒中心频率的困扰。距离－多普勒模型和距离－共面模型都可用于SAR影像的几何处理，表明两者之间有一定的联系。本发明将距离共面方程与距离－多普勒方程结合起来，实现多普勒参数的高精度估计，来提高SAR成像处理的精确度。

发明内容

本发明在于解决目前SAR天线存在较大斜视角、非平稳飞行及复杂地形条件下，SAR波束中心的多普勒参数估计精度较低的不足的问题，给出一种利用POS测量数据以及DEM数据估算SAR信号的多普勒参数，提高非正侧视、载机非平稳飞行和复杂地形条件下的参数的估计精度，对复杂条件下的机载SAR成像有重要的意义。

本发明的技术方案：充分利用机载GPS导航数据、IMU姿态根据和数字高程数据，根据距离方程、共面方程和考虑地面高程数据的地球椭球数据来求解反射信号的目标位置，根据传感器状态和目标位置信息估计波束中心的多普勒信息，即通过计算机进行以下步骤：

步骤1，给定方位向时间，以及POS测量值，插值获得该时刻GPS天线位置、速度、加速度以及传感器平台的姿态，并计算相应时刻SAR传感器天线的位置，具体的说，确定需要求解多普勒参数的相应时刻t，根据载机上GPS和IMU在测量值插值获得该时刻GPS天线位置、速度、加速度以及平台的姿态参数；根据插值获得的GPS天线位置[X_GPS,Y_GPS,Z_GPS]、传感器平台姿态以及在传感器平台坐标系中，SAR天线中心相对于GPS天线中心的偏移量eXs，eYs，eZs，计算相应时刻SAR天线的位置：

$\left[\begin{array}{c}\mathrm{Xs}\\ \mathrm{Ys}\\ \mathrm{Zs}\end{array}\right]=\left[\begin{array}{c}{X}_{\mathrm{GPS}}\\ X_{\mathrm{GPS}}\\ Z_{\mathrm{GPS}}\end{array}\right]+{R_n}^E{R_b}^n\left[\begin{array}{c}\mathrm{eXs}\\ \mathrm{eYs}\\ \mathrm{eZs}\end{array}\right]$

其中，[Xs，Ys，Zs]为计算得到的SAR传感器天线中心位置，R_b ⁿ为平台姿态旋转矩阵，R_n ^E为导航坐标系到地心只直角坐标系的旋转矩阵，该旋转矩阵容易通过传感器所在大地坐标经纬坐标计算得到；

姿态根据相应时刻线性内插POS姿态测量值得到，传感器速度和加速度可以直接采用GPS在相应方位向时间的测量值，也可以通过该时间附近的位置二次多项式拟合函数对时间的一阶和二阶导数求出；

步骤2，给定距离向快时间，建立波束面法向量矢量，构建距离方程、波束中心扫描面方程以及考虑地面高程数据的地球椭球方程。建立扫描波束面法向量矢量

$\stackrel{\→}{i}={R_n}^E{R_b}^n{\left[\begin{array}{ccc}1& 0& 0\end{array}\right]}^T={R_b}^E{\left[\begin{array}{ccc}1& 0& 0\end{array}\right]}^T=\left[\begin{array}{ccc}{a}_1& a_2& a_3\\ b_1& b_2& b_3\\ c_1& c_2& c_3\end{array}\right]{\left[\begin{array}{ccc}1& 0& 0\end{array}\right]}^T={\left[\begin{array}{ccc}{a}_1& b_1& c_1\end{array}\right]}^T$

式中a₁、a₂、a₃、b₁、b₂、b₃、c₁、c₂、c₃为平台坐标系到地心坐标系的旋转矩阵R_b ^E中的元素，R_b ⁿ为传感器平台的姿态旋转矩阵，根据步骤(1)中计算得到的SAR传感器姿态计算，R_n ^E为导航坐标系到地心直角坐标系的旋转矩阵，根据SAR传感器天线中心位置的经纬坐标计算；

距离方程、波束扫描面方程和考虑地面高程数据的地球椭球方程为：

$\left\{\begin{array}{c}{\mathrm{\α}}_1.(X-\mathrm{Xs})+b_1.(Y-\mathrm{Ys})+c_1.(Z-\mathrm{Zs})=0\\ \sqrt{{(X-\mathrm{Xs})}^2+{(Y-\mathrm{Ys})}^2+{(Z-\mathrm{Zs})}^2}=t_r\·c/2\\ \frac{X^2+Y^2}{{(A+H)}^2}+\frac{Z^2}{{(B+H)}^2}=1\end{array}\right.$

式中X，Y，Z为待求的信号反射体目标坐标，Xs，Ys,Zs为步骤(1)中计算得到的SAR传感器天线中心位置坐标，c为光速，t_r为给定的SAR传感器接受信号距离向快时间，H为反射体目标的地面高程，A为地球椭球的长半轴，B为地球椭球的短半轴，取值为A＝6378137.0m，B＝6356752.314m；

步骤3，对距离方程、波束面方程以及考虑地面高程数据的地球椭球方程线性化，构建误差方程组，并利用DEM数据通过迭代计算SAR信号反射体目标的位置坐标，即通过迭代计算相应快时间反射信号的地面目标点坐标的计算；

步骤3.1，其中对距离方程、波束面方程以及考虑地面高程数据的地球椭球方程线性化是通过相应方程对SAR信号反射体目标坐标（X，Y，Z）一阶求导实现，组成的误差方程组为：

$\left\{\begin{array}{c}{v}_c=f_{\mathrm{cX}}\mathrm{dX}+f_{\mathrm{cY}}\mathrm{dY}+f_{\mathrm{cZ}}\mathrm{dZ}-l_c\\ v_r=f_{\mathrm{rX}}\mathrm{dX}+f_{\mathrm{rY}}\mathrm{dY}+f_{\mathrm{rZ}}\mathrm{dZ}-l_r\\ v_e=f_{\mathrm{eX}}\mathrm{dX}+f_{\mathrm{eY}}\mathrm{dY}+f_{\mathrm{eZ}}\mathrm{dZ}-l_e\end{array}\right.$

式中v_c，v_r，v_e，分别为扫描面方程、距离方程以及考虑地面高程数据的地球椭球方程对应的误差增量，f_cX、f_cY、f_cZ，f_rX、f_rY、f_rZ，f_eX、f_eY、f_eZ分别为三个方程对反射体目标坐标的线性化系数，dX，dY，dZ为SAR信号反射体目标坐标未知数增量，l_c，l_r，l_e分别为三个误差方程的常数项；

由于步骤（2）中建立的方程为非线性形式，X、Y、Z坐标的计算通过迭代的方法通过最小二乘求解本步骤中的dX，dY，dZ计算；

步骤3.2，通过迭代计算方位向时间t_a、距离向快时间为t_r的SAR信号反射体目标坐标流程为：

1）首先根据SAR信号接受时的方位向时间t_a，按步骤1中方法计算传感器的位置、速度、加速度及传感器姿态,给定临时高程参数H_i并赋值H_i＝H₀，H₀为测区平均高程；

2）给定步骤2中的高程值H＝H_i，根据接受信号的距离向快时间t_r按步骤3.1建立误差方程组迭代计算反射体坐标增量dX，dY，dZ，进而求解步骤2中的坐标X，Y，Z；

3）将反射体坐标X，Y，Z转换到DEM（数字高程模型）数据相同的坐标系中，得到DEM所属坐标系中的平面坐标（x，y），并根据相应的坐标，获取DEM数据中该点的高程值H_i+1；

4）比较H_i和H_i+1的值，若H_i和H_i+1之差的绝对值小于给定的阈值δ，即|H_i+1-H_i|<δ，流程2)中求得的（X，Y，Z）即为SAR信号反射体目标的坐标，转至步骤4，否则将高程参数H_i重新赋值为H_i＝H_i+1，转至流程2)；

步骤4，根据SAR传感器天线的位置、速度、加速度及SAR信号反射体目标的坐标，计算SAR波束中心信号的多普勒参数，具体的说是根据载机的位置、速度、加速度及相应信号的地面点坐标，计算相应信号的多普勒参数；即根据载机的位置、速度、加速度、地面反射信号点的坐标以及快时间t_r值，获得该时刻SAR波束中心多普勒频率f_DC和调频率f_DR分别为:

f_DC=2[(X-Xs)V_X+(Y-Ys)V_Y+(Z-Zs)V_Z]/(λ·t_r·c)

f_DR=2[(X-Xs)a_X+(Y-Ys)a_Y+(Z-Zs)a_Z]/(λ·t_r·c)(V_X，V_Y，V_Z)为传感器的速度，(a_X，a_Y，a_Z)为传感器的加速度，c为光速，λ为SAR的波长。

步骤5，将SAR数据获取历程的方位向时间区间和距离向快时间区间分别按一定时间间隔构建时间格网，根据(1)－(4)步骤中方法计算每个格网点相应时间SAR波束中心信号的多普勒参数，其它任意时间SAR波束中心信号的多普勒参数根据格网点上的值双线性插值获得。

以方位向时间和距离向时间为轴建立二维坐标系，方位向起始时间t_a0和结束时间t_an内按时间间隔Δt_a进行均匀划分，获取n+1个方位向时间节点，在距离向近距端对应的快时间t_r0和远距端对应的快时间t_rm范围内按时间间隔Δt_r进行均匀划分，获取m+1个时间节点，其中Δt_a=(t_an-t_a0)/n,Δt_r＝(t_rm-t_r0)/m，m、n为整数，大小根据多普勒参数估计精度需求选定；根据方位向和距离向时间节点构建时间格网，方位向第i个时间节点对应的时间t_ai和距离向第j个节点对应的快时间t_rj可以组合成（n+1）×（m+1）个坐标为(t_ai，t_rj)(i＝0,1,2...n;j＝0,1,2...m)的点，称为时间格网点；

根据格网点对应的方位向时间和距离向快时间，按步骤(1)-(4)计算所有格网点时间的多普勒参数，其它距离向和方位向时间（t_a,t_r）处的多普勒参数，根据该时间坐标点周边最近的四个格网点时间的多普勒参数值按双线性插值计算得到。

本发明也可以通过将扫描面的平面方程替换为锥面方程，使发明内容适用于锥面扫描形式的真实孔径雷达或SAR传感器的波束多普勒参数计算。本发明实现了利用POS（即定位定姿***，由GPS和惯性测量装置IMU组成）数据和地形数据来提高SAR信号多普勒参数的估计精度，进而提高SAR成像质量。

附图说明

图1是本发明提供的基于POS与DEM数据的SAR多普勒参数估计方法示意图。

具体实施方式

本发明通过GPS、IMU（惯性测量装置）和地形数据求解地面目标位置，进而求解多普勒参数，具体通过计算机来实施，实施流程参见附图1，多普勒参数估算通过以下步骤实现：

步骤1，确定需要求解多普勒参数的相应时刻t，根据载机上GPS和IMU在测量值插值获得该时刻GPS天线位置、速度、加速度以及平台的姿态参数，以及根据插值获得的GPS天线位置[X_GPS,Y_GPS,Z_GPS]、传感器平台姿态，以及在传感器平台坐标系中，SAR（合成孔径雷达）天线中心相对于GPS天线中心的偏移量[eXs，eYs，eZs]，计算相应时刻SAR天线的位置：

$\left[\begin{array}{c}\mathrm{Xs}\\ \mathrm{Ys}\\ \mathrm{Zs}\end{array}\right]=\left[\begin{array}{c}{X}_{\mathrm{GPS}}\\ X_{\mathrm{GPS}}\\ Z_{\mathrm{GPS}}\end{array}\right]+{R_n}^E{R_b}^n\left[\begin{array}{c}\mathrm{eXs}\\ \mathrm{eYs}\\ \mathrm{eZs}\end{array}\right]$

其中，[Xs，Ys，Zs]为计算得到的SAR传感器天线中心位置，R_b ⁿ为平台姿态旋转矩阵，R_n ^E为导航坐标系到地固地心直角坐标系的旋转矩阵，该旋转矩阵容易通过传感器所在大地坐标经纬坐标计算得到。

平台坐标系以运动载体的几何中心作为原点，运动载体的几何轴作为载体坐标系的坐标轴。通常情况下平台坐标系与传感器坐标系之间存在固定的相对关系。航天传感器的传感器平台坐标系常称作卫星本体坐标系。

导航坐标系也称当地水平坐标系，是以地球椭球面、摄站点处的法线为基准面和基准线建立的局部空间直角坐标系。原点位于载体中心，以X轴取朝北向，Y轴朝东，Z轴取参考椭球外法线反方向(朝向地心)，又称北东地坐标系。

地心坐标系(E-XYZ)与地球固联，随地球一起旋转，也称作地球坐标系。原点O与地球质心重合，Z轴指向北极(由于地极的移动，一般指向国际协议原点CIO，此时也称协议地球坐标系)，X轴指向协议赤道面与格林尼治子午线的交点，Y轴按右手规则确定。

传感器平台坐标系到坐标系到地心直角坐标系的转换矩阵，由IMU测量得到的Ψ,Θ,Ф传感器姿态值来计算，即：

${R_b}^n=R_Z\left(\mathrm{\&Psi;}\right)R_Y\left(\mathrm{\&Theta;}\right)R_X\left(\mathrm{\&Phi;}\right)=\left[\begin{array}{ccc}\cos \mathrm{\&Psi;}& -\sin \mathrm{\&Psi;}& 0\\ \sin \mathrm{\&Psi;}& \cos \mathrm{\&Psi;}& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}\cos \mathrm{\&Theta;}& 0& \sin \mathrm{\&Theta;}\\ 0& 1& 0\\ -\sin \mathrm{\&Theta;}& 0& \cos \mathrm{\&Theta;}\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \mathrm{\&Phi;}& -\sin \mathrm{\&Phi;}\\ 0& \sin \mathrm{\&Phi;}& \cos \mathrm{\&Phi;}\end{array}\right]$

$=\left[\begin{array}{ccc}\cos \mathrm{\&Psi;}\cos \mathrm{\&Theta;}& \cos \mathrm{\&Psi;}\sin \mathrm{\&Theta;}\sin \mathrm{\&Phi;}-\sin \mathrm{\&Psi;}\cos \mathrm{\&Phi;}& \cos \mathrm{\&Psi;}\sin \mathrm{\&Theta;}\cos \mathrm{\&Phi;}+\sin \mathrm{\&Psi;}\sin \mathrm{\&Phi;}\\ \sin \mathrm{\&Psi;}\cos \mathrm{\&Theta;}& \sin \mathrm{\&Psi;}\sin \mathrm{\&Theta;}\sin \mathrm{\&Phi;}+\cos \mathrm{\&Psi;}\cos \mathrm{\&Phi;}& \sin \mathrm{\&Psi;}\sin \mathrm{\&Theta;}\cos \mathrm{\&Phi;}-\cos \mathrm{\&Psi;}\sin \mathrm{\&Phi;}\\ -\sin \mathrm{\&Theta;}& \cos \mathrm{\&Theta;}\sin \mathrm{\&Phi;}& \cos \mathrm{\&Theta;}\cos \mathrm{\&Phi;}\end{array}\right]$

导航坐标系到地心直角坐标系的转换矩阵R_n ^E由传感器的大地经纬度来计算：

${R_n}^E=\left[\begin{array}{ccc}-\sin L& -\cos L\sin B& \cos L\cos B\\ \cos L& -\sin L\sin B& \sin L\cos B\\ 0& \cos B& \sin B\end{array}\right]$

姿态根据相应时刻线性内插POS姿态测量值得到，传感器速度和加速度可以直接采用GPS在相应方位向时间的测量值，也可以通过该时间附近的位置二次多项式拟合函数对时间的一阶和二阶导数求出；

步骤2，建立扫描波束面法向量矢量

$\stackrel{\&RightArrow;}{i}={R_n}^E{R_b}^n{\left[\begin{array}{ccc}1& 0& 0\end{array}\right]}^T={R_b}^E{\left[\begin{array}{ccc}1& 0& 0\end{array}\right]}^T=\left[\begin{array}{ccc}{a}_1& a_2& a_3\\ b_1& b_2& b_3\\ c_1& c_2& c_3\end{array}\right]{\left[\begin{array}{ccc}1& 0& 0\end{array}\right]}^T={\left[\begin{array}{ccc}{a}_1& b_1& c_1\end{array}\right]}^T$

式中a₁、a₂、a₃、b₁、b₂、b₃、c₁、c₂、c₃为平台坐标系到地心坐标系的旋转矩阵R_b ^E中的元素，R_b ⁿ为传感器平台的姿态旋转矩阵，根据步骤(1)中计算得到的SAR传感器姿态计算，R_n ^E为导航坐标系到地心直角坐标系的旋转矩阵，根据SAR传感器天线中心位置的经纬坐标计算；

距离方程、波束扫描面方程和考虑地面高程数据的地球椭球方程为：

$\left\{\begin{array}{c}{a}_1.(X-\mathrm{Xs})+b_1.(Y-\mathrm{Ys})+c_1.(Z-\mathrm{Zs})=0\\ \sqrt{{(X-\mathrm{Xs})}^2+{(Y-\mathrm{Ys})}^2+{(Z-\mathrm{Zs})}^2}=t_r\&CenterDot;c/2\\ \frac{X^2+Y^2}{{(A+H)}^2}+\frac{Z^2}{{(B+H)}^2}=1\end{array}\right.$

式中X，Y，Z为待求的信号反射体目标坐标，Xs,Ys,Zs为步骤(1)中计算得到的SAR传感器天线中心位置坐标，c为光速，t_r为给定的SAR传感器接受信号距离向快时间，H为反射体目标的地面高程，A为地球椭球的长半轴，B为地球椭球的短半轴，取值为A＝6378137.0m，B＝6356752.314m；

步骤3，对距离方程、波束面方程以及考虑地面高程数据的地球椭球方程线性化，构建误差方程组，并通过迭代计算相应快时间反射信号的地面目标点坐标的计算；

而且，步骤3.1，其中对距离方程、波束面方程以及考虑地面高程数据的地球椭球方程线性化是通过相应方程对SAR信号反射体目标坐标（X，Y，Z）一阶求导实现，组成的误差方程组为：

$\left\{\begin{array}{c}{v}_c=f_{\mathrm{cX}}\mathrm{dX}+f_{\mathrm{cY}}\mathrm{dY}+f_{\mathrm{cZ}}\mathrm{dZ}-l_c\\ v_r=f_{\mathrm{rX}}\mathrm{dX}+f_{\mathrm{rY}}\mathrm{dY}+f_{\mathrm{rZ}}\mathrm{dZ}-l_r\\ v_e=f_{\mathrm{eX}}\mathrm{dX}+f_{\mathrm{eY}}\mathrm{dY}+f_{\mathrm{eZ}}\mathrm{dZ}-l_e\end{array}\right.$

式中v_c，v_r，v_e，分别为扫描面方程、距离方程以及考虑地面高程数据的地球椭球方程对应的误差增量，f_cX、f_cY、f_cZ，f_rX、f_rY、f_rZ，f_eX、f_eY、f_eZ分别为三个方程对反射体目标坐标的线性化系数，dX，dY，dZ为SAR信号反射体目标坐标未知数增量，l_c，l_r，l_e分别为三个误差方程的常数项；

由于步骤（2）中建立的方程为非线性形式，X、Y、Z坐标的计算通过迭代的方法通过求解dX，dY，dZ计算；

步骤3.2，通过迭代计算方位向时间t_a、距离向快时间为t_r的SAR信号反射体目标坐标流程为：

1）首先根据SAR信号接受时的方位向时间t_a，按步骤1中方法计算传感器的位置、速度、加速度及传感器姿态,给定临时高程参数H_i并赋值H_i＝H₀，H₀为测区平均高程；

2）给定步骤2中的高程值H＝H_i，根据接受信号的距离向快时间t_r按步骤3.1建立误差方程组迭代计算反射体坐标增量dX，dY，dZ，进而求解步骤2中的坐标X，Y，Z；

3）将反射体坐标X，Y，Z转换到DEM数据相同的坐标系中，得到DEM所属坐标系中的平面坐标（x，y），并根据相应的坐标，获取DEM数据中该点的高程值H_i+1；

4）比较H_i和H_i+1的值，若H_i和H_i+1之差的绝对值小于给定的阈值δ，即|H_i+1-H_i|<δ，流程2)中求得的（X，Y，Z）即为SAR信号反射体目标的坐标，转至步骤4，否则将高程参数H_i重新赋值为H_i＝H_i+1，转至流程2)；

步骤4，根据载机的位置、速度、加速度及SAR信号反射体坐标，计算相应信号的多普勒参数；即根据载机的位置、速度、加速度、SAR信号反射体目标的坐标以及SAR接受信号距离向快时间t_r值，按下式计算SAR信号波束中心多普勒频率f_DC和调频率f_DR:

f_DC=2[(X-Xs)V_X+(Y-Ys)V_Y+(Z-Zs)V_Z]/(λ·t_r·c)

f_DR=2[(X-Xs)a_X+(Y-Ys)a_Y+(Z-Zs)a_Z]/(λ·t_r·c)

(V_X，V_Y，V_Z)、(a_X，a_Y，a_Z)分别为传感器的速度和加速度，c为光速，λ为SAR的波长；

步骤5，以方位向时间和距离向时间为轴建立二维坐标系，方位向起始时间t_a0和结束时间t_an内按时间间隔Δt_a进行均匀划分，获取n+1个方位向时间节点，在距离向近距端对应的快时间t_r0和远距端对应的快时间t_rm范围内按时间间隔Δt_r进行均匀划分，获取m+1个时间节点，其中Δt_a＝(t_an-t_a0)/n,Δt_r=(t_rm-t_r0)/m，m、n为整数，大小根据多普勒参数估计精度需求选定；根据方位向和距离向时间节点构建时间格网，方位向第i个时间节点对应的时间t_ai和距离向第j个节点对应的快时间t_rj可以组合成（n+1）×（m+1）个坐标为(t_ai，t_rj)(i＝0,1,2...n;j＝0,1,2...m)的点，称为时间格网点；

根据格网点对应的方位向时间和距离向快时间，按步骤(1)-(4)计算所有格网点时间的多普勒参数，其它距离向和方位向时间（t_a,t_r）处的多普勒参数，根据该时间坐标点周边最近的四个格网点时间的多普勒参数值按双线性插值计算得到，即：

$\begin{array}{c}f(x,y)\&ap;\frac{f\left(Q_{11}\right)}{(x_2-x_1)(y_2-y_1)}(x_2-x)(y_2-y)+\frac{f\left(Q_{21}\right)}{(x_2-x_1)(y_2-y_1)}(x-x_1)(y_2-y)\\ +\frac{f\left(Q_{12}\right)}{(x_2-x_1)(y_2-y_1)}(x_2-x)(y-y_1)+\frac{f\left(Q_{22}\right)}{(x_2-x_1)(y_2-y_1)}(x-x_1)(y-y_1).\end{array}$

其中，(x₁，y₁)，(x₁，y₂)(x₂，y₁)(x₂，y₂)为待插值点周边最邻近的四个格网点方位和距离时间坐标值。f(Q₁₁)，f(Q₁₂),f(Q₂₁)，f(Q₂₂)为四个格网点的多普勒频率或调频率参数值，f(x，y)则为时间坐标为(x，y)处的插值结果。

本发明实现了利用POS（即定位定姿***，由GPS和惯性测量装置IMU组成）数据和地形数据来提高SAR信号多普勒参数的估计精度，进而提高SAR成像质量。

需要说明的是，本发明在有SAR传感器精确几何定标参数的情况下，有更好的效果。

Claims (7)

1.一种基于POS测量数据与DEM数据的SAR多普勒参数估计方法，其特征在于：利用POS测量数据以及DEM数据估算SAR信号的多普勒参数，提高非正侧视、载机非平稳飞行和复杂地形条件下的参数的估计精度，包括如下步骤：（1）给定方位向时间，以及POS测量值，插值获得该时刻GPS天线位置、速度、加速度以及传感器平台的姿态，并计算相应时刻SAR传感器天线的位置；（2）给定距离向快时间，建立波束面法向量矢量，构建距离方程、波束中心扫描面方程以及考虑地面高程数据的地球椭球方程;（3）对距离方程、波束扫描面方程以及考虑地面高程数据的地球椭球方程线性化，构建误差方程组，并利用DEM数据通过迭代计算SAR信号反射体目标的位置坐标；（4）根据SAR传感器天线的位置、速度、加速度及SAR信号反射体目标的坐标，计算SAR波束中心信号的多普勒参数；(5)将SAR数据获取历程的方位向时间区间和距离向快时间区间分别按一定时间间隔构建时间格网，根据(1)－(4)步骤中方法计算每个格网点相应时间SAR波束中心信号的多普勒参数，其它任意时间SAR波束中心信号的多普勒参数根据格网点上的值双线性插值获得。

2.根据权利要求1所述的基于POS与DEM数据的SAR多普勒参数估计方法，其特征在于：所述步骤（1）中，SAR天线位置[Xs,Ys,Zs]的计算根据插值获得的GPS天线位置坐标[X_GPS，Y_GPS，Z_GPS]、传感器平台姿态以及在传感器平台坐标系中，SAR天线中心相对于GPS天线中心的偏移量向量[eXs,eYs,eZs]^T来计算：

其中，Xs，Ys，Zs为计算得到的SAR传感器天线中心位置坐标，R_b ⁿ为传感器平台姿态的旋转矩阵，R_n ^E为导航坐标系到地心直角坐标系的旋转矩阵，该旋转矩阵通过传感器所在大地坐标经纬坐标计算得到； 

姿态根据相应时刻线性内插POS姿态测量值得到，传感器速度和加速度可以直接采用GPS在相应方位向时间的测量值，也可以通过该时间附近的位置二次多项式拟合函数对时间的一阶和二阶导数求出；

该步获得的传感器运动状态参数用于步骤（2）（3）中SAR信号反射体目标（即反射SAR信号的地面目标）坐标的计算。

3.根据权利要求1所述的基于POS与DEM数据的SAR多普勒参数估计方法，其特征在于：所述步骤（2）中，扫描波束面法向量矢量 为：

式中a₁、a₂、a₃、b₁、b₂、b₃、c₁、c₂、c₃为平台坐标系到地心坐标系的旋转矩阵R_b ^E中的元素，R_b ⁿ为传感器平台的姿态旋转矩阵，根据步骤(1)中计算得到的SAR传感器姿态计算，R_n ^E为导航坐标系到地心直角坐标系的旋转矩阵，根据SAR传感器天线中心位置的经纬坐标计算；

距离方程、波束扫描面方程和考虑地面高程数据的地球椭球方程为：

式中X，Y，Z为待求的信号反射体目标（即地面目标）坐标，Xs,Ys,Zs为步骤(1)中计算得到的SAR传感器天线中心位置坐标，c为光速，t_r为给定的SAR传感器接受信号距离向快时间，H为给定的地面高程，A为地球椭球的长半轴，B为地球椭球的短半轴，取值为A＝6378137.0m，B＝6356752.314m。

4.根据权利要求1所述的基于POS与DEM数据的SAR多普勒参数估计方法，其特征在于：所述步骤（3）中，对距离方程、波束面方程以及考虑地面高程数据的地球椭球方程线性化是通过相应方程对SAR信号反射体目标坐标（X，Y，Z）一阶求导实现，组成的误差方程组为： 

式中v_c，v_r，v_e，分别为扫描面方程、距离方程以及考虑地面高程数据的地球椭球方程对应的误差增量，f_cX、f_cY、f_cZ，f_rX、f_rY、f_rZ，f_eX、f_eY、f_eZ分别为三个方程对反射体目标坐标的线性化系数，dX，dY，dZ为SAR信号反射体目标坐标未知数增量，l_c，l_r，l_e分别为三个误差方程的常数项；

由于步骤（2）中建立的方程为非线性形式，X、Y、Z坐标的计算通过迭代的方法通过求解dX，dY，dZ计算。

5.根据权利要求1所述的基于POS与DEM数据的SAR多普勒参数估计方法，其特征在于：所述步骤（3）中，SAR信号反射体的位置坐标通过下述流程计算：

1），首先根据SAR信号接受时的方位向时间t_a，按步骤（1）中方法计算传感器的位置、速度、加速度及传感器姿态,给定临时高程参数H_i并赋值H_i＝H₀，H₀为测区平均高程；

2），给定步骤（2）中的高程值H＝H_i，根据接受信号的距离向快时间t_r按步骤（3）建立误差方程组迭代计算反射体坐标增量dX，dY，dZ，进而求解步骤（2）中的坐标X，Y，Z；

3），将反射体坐标X，Y，Z转换到DEM数据相同的坐标系中，得到DEM所属坐标系中的平面坐标（x，y），并根据相应的坐标，获取DEM数据中该点的高程值H_i+1；

4），比较H_i和H_i+1的值，若H_i和H_i+1之差的绝对值小于给定的阈值δ，即|H_i+1-H_i|<δ，流程2)中求得的地面点坐标（X，Y，Z）即为SAR信号反射体目标的坐标，转至步骤（4），否则将高程参数H_i重新赋值为H_i＝H_i+1，转至流程2)；

获取的信号反射体坐标用于步骤（4）多普勒参数的计算。

6.根据权利要求1所述的基于POS与DEM数据的SAR多普勒参数估计方法，其特征在于：所述步骤（4）中，根据载机的位置、速度、加速度、SAR信号反射体的坐标以及距离向快时间t_r值，按下式计算该时刻SAR波束中心多普勒频率f_DC和调频率f_DR: 

f_DC=2[(X-Xs)V_X+(Y-Ys)V_Y+(Z-Zs)V_Z]/(λ·t_r·c)

f_DR=2[(X-Xs)a_X+(Y-Ys)a_Y+(Z-Zs)a_Z]/(λ·t_r·c)(V_X，V_Y，V_Z)为传感器的速度，(a_X，a_Y，a_Z)为传感器的加速度，c为光速，λ为SAR的波长。

7.根据权利要求1所述的基于POS与DEM数据的SAR多普勒参数估计方法，其特征在于，所述步骤(5)中，时间格网的建立和格网点时间多普勒参数计算方法为：以方位向时间和距离向时间为轴建立二维坐标系，方位向起始时间t_a0和结束时间t_an内按时间间隔Δt_a进行均匀划分，获取n+1个方位向时间节点，在距离向近距端对应的快时间t_r0和远距端对应的快时间t_rm范围内按时间间隔Δt_r进行均匀划分，获取m+1个时间节点，其中Δt_a＝(t_an-t_a0)/n,Δt_r＝(t_rm-t_r0)/m，m、n为整数，m和n的大小根据多普勒参数估计精度需求选定，值越大，后面的内插结果精度越高，反之越低；根据方位向和距离向时间节点构建时间格网，方位向第i个时间节点对应的时间t_ai和距离向第j个节点对应的快时间t_rj可以组合成（n+1）×（m+1）个坐标为(t_ai,t_rj)(i＝0,1,2...n;j＝0,1,2...m)的点，称为时间格网点；

根据格网点对应的方位向时间和距离向快时间，按步骤(1)-(4)计算所有格网点时间的多普勒参数，其它距离向和方位向时间（t_a,t_r）处的多普勒参数，根据该时间坐标点周边最近的四个格网点时间的多普勒参数值按双线性插值计算得到。 

Images