CN102426353B - 一种利用高精度POS离线获取机载InSAR运动误差的方法 - Google Patents

Abstract

一种利用高精度POS离线获取机载InSAR运动误差的方法，该方法利用高精度POS中SINS与GPS的组合测量结果获取SINS中心的二次及二次以下项位置信息，利用SINS的纯惯性测量结果获取SINS中心的二次以上项位置信息，将计算得到的二次及二次以下项位置信息和二次以上项位置信息相加并补偿杆臂效应误差，计算出机载InSAR两个天线相位中心的高精度平滑位置信息，再将高精度平滑位置信息减去其一次多项式拟合结果，最终得到机载InSAR两个天线的运动误差。本发明具有精度高、计算简单的特点，可用于离线计算机载高精度InSAR的运动误差，以实现机载InSAR的高精度离线成像。

Description

一种利用高精度POS离线获取机载InSAR运动误差的方法

技术领域

本发明涉及一种利用位置测量精度达到厘米级、姿态测量精度达到10″量级的高精度位置姿态测量系统(Position and Orientation System，POS)获取机载InSAR运动误差的方法。

背景技术

机载合成孔径雷达(Synthetic Aperture Radar，SAR)是一种先进的对地观测的手段。相对于光学对地观测设备，机载SAR具有全天候、全天时、作用距离远等特点，并对地面植被、人工伪装甚至地表土壤层具有一定的穿透性，日益成为当今世界最具有代表性的对地观测手段之一。尤其是机载高精度InSAR，能够进行立体测绘，在国家基础测绘、国土资源调查、地质资源勘探、农业农情监测、自然灾害监测以及军事侦察等领域具有重要用途，是国家十分重要的军民两用的高技术装备。

机载SAR成像的基本原理要求载机作匀速直线飞行，即要求在合成孔径时间内载机所受到的合力为零。但由于大气湍流、设备性能等因素的影响，载机沿飞行方向和侧向必然存在干扰力的作用，导致机载SAR的天线相位中心(Anten na Phase Center，APC)偏离理想的匀速直线平移运动。常见的情况是载机在自动驾驶仪的控制下围绕理想航迹来回摆动，从而产生运动误差。这对SAR成像分辨率具有不可忽视的影响，严重时甚至导致SAR无法成像。因此，在机载高分辨率SAR成像过程中必须进行运动误差补偿。

机载SAR运动误差补偿方法有两种，一种是基于雷达回波数据的自聚焦方法，另一种是基于POS的运动误差测量方法。目前，比较成熟的方法是自聚焦方法与POS相结合，利用自聚焦方法获取低次运动误差，利用POS获取高次运动误差。但是，对于机载InSAR来说，采用自聚焦方法会破坏干涉条纹，使得InSAR无法实现高程精度测量。因此，机载InSAR不能够采用自聚焦与POS相结合的运动误差补偿方法，迫切需要一种能够精确获取机载InSAR运动误差的方法。

发明内容

本发明的技术解决问题是：克服现有技术的不足，提供一种利用高精度POS离线获取机载InSAR运动误差的方法，该方法充分利用了SINS/GPS组合测量和SINS纯惯性测量各自的优点，提高运动误差的精度，为机载InSAR运动误差补偿提供了一条有效途径。

本发明的技术解决方案为：一种基于高精度POS的机载InSAR运动误差离线计算方法，利用高精度POS中SINS与GPS的组合测量结果获取SINS中心的二次及二次以下项位置信息，利用SINS的纯惯性测量结果获取SINS中心的二次以上项位置信息，将计算得到的二次及二次以下项位置信息和二次以上项位置信息相加并补偿杆臂效应误差，计算出机载InSAR两个天线相位中心的高精度平滑位置信息，再将高精度平滑位置信息减去其一次多项式拟合结果，最终得到机载InSAR两个天线的运动误差。

具体实现步骤如下：

(1)在获取到POS记录的SINS数据和GPS数据的前提下，采用“前向卡尔曼滤波+后向平滑”的方法对POS记录的SINS数据和GPS数据进行信息融合，得到SINS中心的组合测量结果，包含位置、速度和姿态；

(2)确定机载InSAR第i次成像的开始时间T_i0和结束时间T_i1；

(3)在T_i0时刻，以组合测量结果为初值，启动捷联解算程序B，到T_i1时刻结束捷联解算程序B，计算出[T_i0，T_i1]时间段内SINS中心的纯惯性测量结果，包含位置、速度和姿态；

(4)将[T_i0，T_i1]时间段内组合测量的位置进行二次多项式拟合，组合测量位置的二次拟合结果即为[T_i0，T_i1]时间段内SINS中心的二次及二次以下项位置信息；

(5)将[T_i0，T_i1]时间段内纯惯性测量位置进行二次多项式拟合，再利用纯惯性测量位置减去其二次拟合结果，得到[T_i0，T_i1]时间段内SINS中心的二次以上项位置信息；

(6)将SINS中心的二次以上项位置误差和二次及二次以下项位置误差求和，即得到SINS中心的精确位置信息；

(7)采用如下公式补偿SINS中心与两个InSAR天线相位中心之间的杆臂误差，计算得到机载InSAR两个天线相位中心的高精度平滑位置信息：

x_i＝x+Δx_i    (i＝1，2)

y_i＝y+Δy_i    (i＝1，2)

z_i＝z+Δz_i    (i＝1，2)

其中，x、y、z表示SINS中心的位置信息，x_i、y_i、z_i为第i个InSAR天线相位中心的位置信息，Δx_i、Δy_i、Δz_i为SINS中心到第i个InSAR天线相位中心的杆臂误差；

(8)将步骤七得到高精度平滑位置信息减去其一次多项式拟合结果，得到机载InSAR两个天线的运动误差。

本发明的原理是：利用高精度POS中SINS与GPS的组合测量结果获取SINS中心的二次及二次以下项位置信息，利用SINS的纯惯性测量结果获取SINS中心的二次以上项位置信息，将计算得到的二次及二次以下项位置信息和二次以上项位置信息相加并补偿杆臂效应误差，计算出机载InSAR两个天线相位中心的高精度平滑位置信息，再将高精度平滑位置信息减去其一次多项式拟合结果，最终得到机载InSAR两个天线的运动误差。

本发明与现有技术相比的优点在于：

(1)本发明充分发挥了SINS/GPS组合测量结果绝对精度高但有高次误差和SINS纯惯性测量结果误差积累但数据平滑的特点，通过组合测量结果与纯惯性测量结果的结合，获取了机载InSAR天线相位中心的高精度、平滑信息，再通过减去其一次拟合结果得到了机载InSAR天线相位中心的运动误差。

(2)本发明具有精度高(可达到毫米级相对位置误差)、计算简单的特点，可用于离线计算机载高精度InSAR的运动误差，以实现机载InSAR的高精度离线成像。

附图说明

图1为本发明的一种基于高精度POS的机载InSAR运动误差离线计算方法流程图；

图2为本发明中涉及的卡尔曼滤波解算流程图；

图3为本发明中涉及的纯惯性解算流程图。

具体实施方式

如图1所示，本发明的具体实施包括以下步骤：

1、在获取到POS记录的SINS数据和GPS数据的前提下，采用“前向卡尔曼滤波+后向平滑”的方法对POS记录的SINS数据和GPS数据进行信息融合，得到SINS中心的组合测量结果，包含位置、速度和姿态，所述的“前向卡尔曼滤波+后向平滑”数据处理方法的思路是：在时间T₀-T_n内，使用前向卡尔曼滤波并存储各个时刻的估计和方差阵(包含 P_k、P_k+1，k、K_k、Φ_k+1，k)，当前向卡尔曼滤波结束后，采用平滑方程并用存储的全部估计与方差阵信息反顺序得到平滑的和P_k/n；

①所述的前向卡尔曼滤波具体公式如下：

●系统状态方程：

$\stackrel{\·}{X}=\mathrm{FX}+\mathrm{GW}---\left(1\right)$

其中，X为系统状态矢量，W为系统噪声矢量，F为系统转移矩阵，G为噪声转换矩阵：

$X={\left[\begin{array}{ccccccccccccccc}{\mathrm{\φ}}_x& {\mathrm{\φ}}_y& {\mathrm{\φ}}_z& {\mathrm{\δv}}_x& {\mathrm{\δv}}_y& {\mathrm{\δv}}_z& \mathrm{\δL}& \mathrm{\δ\λ}& \mathrm{\δh}& {\mathrm{\ϵ}}_x& {\mathrm{\ϵ}}_y& {\mathrm{\ϵ}}_z& {\▿}_x& {\▿}_y& {\▿}_z\end{array}\right]}^T$

$F=\left[\begin{array}{cc}{F}_{\mathrm{INS}}& F_S\\ o_{6\×6}& F_M\end{array}\right],$ $F_S=\left[\begin{array}{cc}{C}_b^n& 0_{3\×3}\\ 0_{3\×3}& C_b^n\\ 0_{3\×3}& 0_{3\×3}\end{array}\right],$ F_M＝[0_6×15]， $G=\left[\begin{array}{cc}{C}_b^n& 0_{3\×3}\\ 0_{3\×3}& C_b^n\\ 0_{9\×3}& 0_{9\×3}\end{array}\right]$

其中，φ_x、φ_y和φ_z为平台失准角，δv_x、δv_y和δv_z为速度误差，δL、δλ和δh分别为纬度误差、经度误差和高度误差，ε_x、ε_y和ε_z为陀螺仪随机常值漂移误差，Δ_x、Δ_y和Δ_z为加速度计随机常值偏置误差，F_INS为惯性系统矩阵，

为系统姿态转换矩阵；

●系统的量测方程

Z＝HX+η(2)

其中：Z为观测矢量，H为观测矩阵，η为量测噪声，I为单位矩阵：

Z＝[δlδλδh δv_x δv_y δv_z]^T

$H=\left[\begin{array}{ccc}{0}_{3\×6}& I_{3\×3}& 0_{3\×6}\\ 0_{3\×3}& I_{3\×3}& 0_{3\×9}\end{array}\right]$

●卡尔曼滤波基本算法编排，该算法的流程图如图2所示：

状态一步预测方程

${\stackrel{\mathrm{\Λ}}{X}}_{k/k-1}={\mathrm{\φ}}_{k,k-1}{\stackrel{\mathrm{\Λ}}{X}}_{k-1}---\left(3\right)$

其中，

φ_k，k-1分别为k时刻系统状态一步预测值、k-1时刻系统状态估值、k-1时刻到k时刻的系统状态转移矩阵；

状态估值计算方程

${\stackrel{\mathrm{\Λ}}{X}}_k={\stackrel{\mathrm{\Λ}}{X}}_{k/k-1}+K_k(Z_k-H_k{\stackrel{\mathrm{\Λ}}{X}}_{k/k-1})---\left(4\right)$

其中，

K_k、Z_k、H_k分别为k时刻系统状态估值、系统增益矩阵、量测向量和量测矩阵；

滤波增量方程

$K_k=P_{k/k-1}{H}_k^T{(H_k{P}_{k/k-1}{H}_k^T+R_k)}^{-1}---\left(5\right)$

其中，P_k/k-1、R_k分别为k时刻系统协方差阵的一步预测、k时刻系统量测噪声矩阵；

一步预测均方误差方程

$P_{k/k-1}={\mathrm{\φ}}_{k,k-1}{P}_{k-1}{\mathrm{\φ}}_{k,k-1}^T+{\mathrm{\Γ}}_{k-1}{Q}_{k-1}{\mathrm{\Γ}}_{k-1}^T---\left(6\right)$

其中，P_k-1、Q_k-1、Г_k-1分别为k-1时刻系统协方差阵、k-1时刻系统噪声矩阵、k-1时刻系统噪声驱动矩阵；

估计均方误差方程

$P_k=(I-K_k{H}_k)P_{k/k-1}{(I-K_k{H}_k)}^T+K_k{R}_k{K}_k^T---\left(7\right)$

其中，P_k为k时刻系统状态协方差阵；

②所述的后向平滑公式编排：

${\hat{X}}_{k/n}={\hat{X}}_k+P_k{\mathrm{\φ}}_{k+1,k}^T{P}_{k,k-1}^{-1}({\hat{X}}_{k+1/n}-{\hat{X}}_{k+1,k})---\left(8\right)$

$P_{k/n}=P_k+A_k(P_{k+1/n}-P_{k+1,k})A_k^T---\left(9\right)$

$A_k=P_k{\mathrm{\φ}}_{k+1,k}^T{P}_{k+1,k}^{-1}---\left(10\right)$

其中，

P_k/n、P_k+1/n分别为利用T₀-T_n时刻所有前向滤波结果平滑得到的k时刻状态估计值、k+1时刻状态估计值、k时刻系统协方差阵和k+1时刻系统协方差阵，

P_k、P_k+1，k、P_k，k-1、Φ_k+1，k分别为采用前向卡尔曼滤波得到的k时刻状态估计值，k+1时刻状态一步预测值、k时刻系统状态协方差、k+1时刻系统状态协方差的一步预测、k时刻系统状态协方差的一步预测、k时刻到k+1时刻的系统状态转移矩阵，A_k为平滑时的增益矩阵。

2、确定机载InSAR第i次成像的开始时间T_i0和结束时间T_i1；

3、在T_i0时刻，以组合测量结果为初值，进行捷联解算，到T_i1时刻结束捷联解算，计算出[T_i0，T_i1]时间段内SINS中心的纯惯性测量结果，包含位置、速度和姿态，所述的捷联解算公式为：

①在WGS-84全球大地坐标系下，载体所处位置的重力加速度g的数值为：

$g=9.7803267714\×(1-2\frac{h}{R_{e}h})\×\frac{(1+0.00193185138639{\sin }^{2}L)}{{(1-0.00669437999013{\sin }^{2}L)}^{1/2}}---\left(11\right)$

其中，h为高度，R_e为地球半径，L为地理纬度；

②姿态矩阵

的计算公式为：

$C_n^b=\left[\begin{array}{ccc}\cos \mathrm{\γ}\cos \mathrm{\ψ}-\sin \mathrm{\γ}\sin \mathrm{\θ}\sin \mathrm{\ψ}& \cos \mathrm{\γ}\sin \mathrm{\ψ}+\sin \mathrm{\γ}\sin \mathrm{\θ}\cos \mathrm{\ψ}& -\sin \mathrm{\γ}\cos \mathrm{\θ}\\ -\cos \mathrm{\θ}\sin \mathrm{\ψ}& \cos \mathrm{\θ}\cos \mathrm{\ψ}& \sin \mathrm{\θ}\\ \sin \mathrm{\γ}\cos \mathrm{\ψ}+\cos \mathrm{\γ}\sin \mathrm{\θ}\sin \mathrm{\ψ}& \sin \mathrm{\γ}\sin \mathrm{\ψ}-\cos \mathrm{\γ}\sin \mathrm{\θ}\cos \mathrm{\ψ}& \cos \mathrm{\γ}\cos \mathrm{\θ}\end{array}\right]---\left(12\right)$

其中，ψ、θ、γ分别为航向角、俯仰角和横滚角；

③位置矩阵

的计算公式为：

$C_e^n=\left[\begin{array}{ccc}-\sin \mathrm{\λ}& \cos \mathrm{\λ}& 0\\ -\sin L\cos \mathrm{\λ}& -\sin L\sin \mathrm{\λ}& \cos L\\ \cos L\cos \mathrm{\λ}& \cos L\sin \mathrm{\λ}& \sin L\end{array}\right]---\left(13\right)$

其中，L和λ分别为地理纬度和经度；

④四元数q计算公式为

$q=\left[\begin{array}{c}{q}_0\\ q_1\\ q_2\\ q_3\end{array}\right]=\left[\begin{array}{c}\cos \frac{\mathrm{\ψ}}{2}\cos \frac{\mathrm{\θ}}{2}\cos \frac{\mathrm{\γ}}{2}-\sin \frac{\mathrm{\ψ}}{2}\sin \frac{\mathrm{\θ}}{2}\sin \frac{\mathrm{\γ}}{2}\\ \cos \frac{\mathrm{\ψ}}{2}\sin \frac{\mathrm{\θ}}{2}\cos \frac{\mathrm{\γ}}{2}-\sin \frac{\mathrm{\ψ}}{2}\cos \frac{\mathrm{\θ}}{2}\sin \frac{\mathrm{\γ}}{2}\\ \cos \frac{\mathrm{\ψ}}{2}\cos \frac{\mathrm{\θ}}{2}\sin \frac{\mathrm{\γ}}{2}+\sin \frac{\mathrm{\ψ}}{2}\sin \frac{\mathrm{\θ}}{2}\cos \frac{\mathrm{\γ}}{2}\\ \cos \frac{\mathrm{\ψ}}{2}\sin \frac{\mathrm{\θ}}{2}\sin \frac{\mathrm{\γ}}{2}+\sin \frac{\mathrm{\ψ}}{2}\cos \frac{\mathrm{\θ}}{2}\cos \frac{\mathrm{\γ}}{2}\end{array}\right]---\left(14\right)$

⑤使用角增量法计算姿态更新矩阵，四阶计算公式为：

$q(n+1)=\{(1-\frac{{\left(\mathrm{\Δ}{\mathrm{\θ}}_0\right)}^2}{8}+\frac{{\left({\mathrm{\Δ\θ}}_0\right)}^4}{384})I+(\frac{1}{2}-{\left(\frac{{\mathrm{\Δ\θ}}_0}{48}\right)}^2\left(\mathrm{\Δ\θ}\right))\}q\left(n\right)---\left(15\right)$

其中，

$\mathrm{\Δ\θ}=\left[\begin{array}{cccc}0& -{\mathrm{\Δ\θ}}_x& -{\mathrm{\Δ\θ}}_y& -\mathrm{\Δ}{\mathrm{\θ}}_z\\ {\mathrm{\Δ\θ}}_x& 0& {\mathrm{\Δ\θ}}_z& -{\mathrm{\Δ\θ}}_y\\ {\mathrm{\Δ\θ}}_y& -{\mathrm{\Δ\θ}}_z& 0& {\mathrm{\Δ\θ}}_x\\ {\mathrm{\Δ\θ}}_z& {\mathrm{\Δ\θ}}_y& -{\mathrm{\Δ\θ}}_x& 0\end{array}\right]$

${\mathrm{\Δ\θ}}_0^2={\mathrm{\Δ\θ}}_x^2+{\mathrm{\Δ\θ}}_y^2+{\mathrm{\Δ\θ}}_z^2$

⑥速度计算公式为：

$\left[\begin{array}{c}{\stackrel{\·}{v}}_x^n\\ {\stackrel{\·}{v}}_y^n\\ {\stackrel{\·}{v}}_z^n\end{array}\right]=\left[\begin{array}{c}{a}_{\mathrm{ibx}}^n\\ a_{\mathrm{iby}}^n\\ a_{\mathrm{ibz}}^n\end{array}\right]+\left[\begin{array}{ccc}0& {2\mathrm{\ω}}_{\mathrm{iez}}^n+{\mathrm{\ω}}_{\mathrm{enz}}^n& -({2\mathrm{\ω}}_{\mathrm{iey}}^n+{\mathrm{\ω}}_{\mathrm{eny}}^n)\\ -(2{\mathrm{\ω}}_{\mathrm{iez}}^n++{\mathrm{\ω}}_{\mathrm{enz}}^n)& 0& {2\mathrm{\ω}}_{\mathrm{iex}}^n+{\mathrm{\ω}}_{\mathrm{enx}}^n\\ 2{\mathrm{\ω}}_{\mathrm{iey}}^n+{\mathrm{\ω}}_{\mathrm{eny}}^n& -({2\mathrm{\ω}}_{\mathrm{iex}}^n+{\mathrm{\ω}}_{\mathrm{enx}}^n)& 0\end{array}\right]\left[\begin{array}{c}{v}_x^n\\ v_y^n\\ v_z^n\end{array}\right]-\left[\begin{array}{c}0\\ 0\\ g\end{array}\right]---\left(16\right)$

其中，

和

为东北天地理坐标系下的相对加速度，

和

为东北天地理坐标系下的非引力加速度，

和

为东北天地理坐标系下的速度，

和

为地球自转角速度在东北天坐标系的投影，

和

为由导航坐标系相对地球的旋转角速度在东北天坐标系的投影；

⑦导航位置矩阵

或

更新计算公式为：

${\stackrel{\·}{C}}_n^e=C_n^e{\mathrm{\Ω}}_{\mathrm{en}}^n---\left(17\right)$

其中：

${\mathrm{\Ω}}_{\mathrm{en}}^n=\left[\begin{array}{ccc}0& -{\mathrm{\ω}}_{\mathrm{enz}}^n& {\mathrm{\ω}}_{\mathrm{eny}}^n\\ {\mathrm{\ω}}_{\mathrm{enz}}^0& 0& -{\mathrm{\ω}}_{\mathrm{enx}}^n\\ -{\mathrm{\ω}}_{\mathrm{eny}}^n& {\mathrm{\ω}}_{\mathrm{enx}}^n& 0\end{array}\right]---\left(18\right)$

4、将[T_i0，T_i1]时间段内组合测量的位置进行二次多项式拟合，组合测量位置的二次拟合结果即为[T_i0，T_i1]时间段内SINS中心的二次及二次以下项位置信息(包括Lonz2、Latz2、Altz2)，具体步骤为：

①采用Matlab软件读入[T_i0，T_i1]时间段内组合测量的位置数据，并将经度、纬度、高度、时间分别用4个n×1维列向量Lonz、Latz、Altz和T表示，其中n为[T_i0，T_i1]时间组合测量位置数据的个数；

②利用Matlab软件中的polyfit函数得到[T_i0，T_i1]时间段内组合测量位置数据与时间的二次函数，即得到该二次函数中的各项系数，具体命令如下：

Aloz＝polyfit(T，Lonz，2)；

Alaz＝polyfit(T，Latz，2)；

Atz＝polyfit(T，Altz，2)；

其中，三个向量Aloz、Alaz和Atz分别表示经度、纬度和高度的时间二次函数的各项系数；

③利用上一步得到的二次多项式的各项系数，采用Matlab软件中的polyval函数得到[T_i0，T_i1]时间段内组合测量经度、纬度和高度的二次拟合结果，分别用Lonz2、Latz2和Altz2表示，具体命令如下：

Lonz2＝polyval(Aloz，T)；

Latz2＝polyval(Alaz，T)；

Altz2＝polyval(Atz，T)；

5、将[T_i0，T_i1]时间段内纯惯性测量位置进行二次多项式拟合，再利用纯惯性测量位置减去其二次拟合结果，得到[T_i0，T_i1]时间段内SINS中心的二次以上项位置信息，具体步骤为：

①采用Matlab软件读入[T_i0，T_i1]时间段内纯惯性测量的位置数据，并将经度、纬度、高度、时间分别用4个n×1维列向量Loni、Lati、Alti和T表示，其中n为[T_i0，T_i1]时间纯惯性测量位置数据的个数；

②利用Matlab软件中的polyfit函数得到[T_i0，T_i1]时间段内纯惯性测量位置数据(经度、纬度和高度)与时间的二次函数，即得到该二次函数中的各项系数，分别用Aloi、Alai和Ati表示，具体命令如下：

Aloi＝polyfit(T，Loni，2)；

Alai＝polyfit(T，Lati，2)；

Ati＝polyfit(T，Alti，2)；

③利用上一步得到的二次多项式的各项系数，采用Matlab软件中的polyval函数得到[T_i0，T_i1]时间段内纯惯性测量位置数据的二次拟合结果，即Loni2、Lati2和Alti2，具体命令如下：

Loni2＝polyval(Aloi，T)；

Lati2＝polyval(Alai，T)；

Alti2＝polyval(Ati，T)；

④利用纯惯性测量位置信息及其二次拟合结果，前者减去后者得到纯惯性测量的二次以上项位置信息Loni3、Lati3和Alti3，具体公式如下：

Loni3＝Loni-Loni2；

Lati3＝Lati-Lati2；

Alti3＝Alti-Alti2；

其中，Lon3、Lat3和Alt3为利用纯惯性测量数据获得的SINS中心二次以上项经度、纬度和高度，Loni2、Lati2和Alti2为利用纯惯性测量数据获得的SINS中心二次及二次以下项经度、纬度和高度，Loni、Lati和Alti为纯惯性测量得到的经度、纬度和高度。

6、将SINS中心的二次以上项位置信息和二次及二次以下项位置信息求和，即得到SINS中心的精确位置信息Lon、Lat和Alt，具体公式如下：

Lon＝Lonz2+Loni3；

Lat＝Latz2+Lati3；

Alt＝Altz2+Alti3；

其中，Lon、Lat和Alt为最终得到的SINS中心的经度、纬度和高度信息，Lon2、Lat2和Alt2为利用组合测量数据获得的SINS中心二次及二次以下项经度、纬度和高度，Lon3、Lat3和Alt3为利用纯惯性测量数据获得的SINS中心二次以上项经度、纬度和高度。

7、采用如下公式补偿SINS中心与两个InSAR天线相位中心之间的杆臂误差，计算得到机载InSAR两个天线相位中心的高精度平滑位置信息：

x_i＝x+Δx_i    (i＝1，2)(1)

y_i＝y+Δy_i    (i＝1，2)(2)

z_i＝z+Δz_i    (i＝1，2)(3)

其中，x、y、z表示SINS中心的位置信息，x_i、y_i、z_i为第i个InSAR天线相位中心的位置信息，Δx_i、Δy_i、Δz_i为SINS中心到第i个InSAR天线相位中心的杆臂误差；

8、将步骤7得到高精度平滑位置信息减去其一次多项式拟合结果，得到机载InSAR两个天线的运动误差，具体公式如下：

①利用Matlab软件中的polyfit函数得到[T_i0，T_i1]时间段内高精度平滑位置信息(包含经度Lon、纬度Lat和高度Alt)与时间的一次函数，即得到该一次函数中的各项系数，分别用Alo1、Ala1和At1表示，具体命令如下：

Alo1＝polyfit(T，Lon，1)；

Ala1＝polyfit(T，Lat，1)；

At1＝polyfit(T，Alt，1)；

②利用上一步得到的二次多项式的各项系数，采用Matlab软件中的polyval函数得到[T_i0，T_i1]时间段内高精度平滑位置信息(经度、纬度和高度)的一次拟合结果，即Lon1、Lat1和Alt1，具体命令如下：

Lon1＝polyval(Alo1，T)；

Lat1＝polyval(Ala1，T)；

Alt1＝polyval(At1，T)；

③利用高精度平滑位置信息及其一次拟合结果，前者减去后者得到I nSAR两个天线的运动误差ΔLon、ΔLat和ΔAlti，具体公式如下：

ΔLon＝Lon-Lon1；

ΔLat＝Lat-Lat1；

ΔAlt＝Alt-Alt1；

其中，ΔLon、ΔLat和ΔAlti分别表示InSAR两个天线的经度误差、纬度误差和高度误差，Lon、Lat和Alti为高精度平滑得到的InSAR两个天线的经度、纬度和高度，Lon1、Lat1和Alti1为利用Lon、Lat和Alti进行一次拟合得到的经度、纬度和高度。

本发明说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。

Claims (2)

1.一种利用高精度POS离线获取机载InSAR运动误差的方法，其特征在于实现步骤如下：

(1)在获取到POS记录的SINS数据和GPS数据的前提下，采用“前向卡尔曼滤波+后向平滑”的方法对POS记录的SINS数据和GPS数据进行信息融合，得到SINS中心的组合测量结果，包含位置、速度和姿态；

(2)确定机载InSAR第i次成像的开始时间T_i0和结束时间T_i1；

(3)在T_i0时刻，以组合测量结果为初值，进行捷联解算，到T_i1时刻结束捷联解算，计算出[T_i0，T_i1]时间段内SINS中心的纯惯性测量结果，包含位置、速度和姿态；

(4)将[T_i0，T_i1]时间段内组合测量的位置进行二次多项式拟合，组合测量位置的二次拟合结果即为[T_i0，T_i1]时间段内SINS中心的二次及二次以下项位置信息；

(5)将[T_i0，T_i1]时间段内纯惯性测量位置进行二次多项式拟合，再利用纯惯性测量位置减去其二次拟合结果，得到[T_i0，T_i1]时间段内SINS中心的二次以上项位置信息；

(6)将SINS中心的二次以上项位置信息和二次及二次以下项位置信息求和，即得到SINS中心的精确位置信息；

(7)补偿SINS中心与两个InSAR天线相位中心之间的杆臂误差，计算得到机载InSAR两个天线相位中心的高精度平滑位置信息；

(8)将步骤(7)得到高精度平滑位置信息减去其一次多项式拟合结果，得到机载InSAR两个天线的运动误差。

2.根据权利要求1所述的一种利用高精度POS离线获取机载InSAR运动误差的方法，其特征在于：所述步骤(7)中补偿SINS中心与两个InSAR天线相位中心之间的杆臂误差的公式如下：

x_i=x+Δx_i

y_i=y+Δy_i

z_i=z+Δz_i

其中，i＝1,2；x、y、z表示SINS中心的精确位置信息；x_i、y_i、z_i为第i个InSAR天线相位中心的高精度平滑位置信息；Δx_i、Δy_i、Δz_i为SINS中心到第i个InSAR天线相位中心的杆臂误差。

Images