CN103344252B - 一种航空高光谱成像系统定位误差分析方法 - Google Patents

Abstract

一种航空高光谱成像系统定位误差分析方法，其步骤如下：（1）高光谱图像信息、姿态位置测量系统数据读入；（2）根据几何成像物理关系建立共线条件方程；（3）分析姿态位置测量系统测角误差、全球定位系统定位误差和成像光谱仪标定误差，对每个误差源进行建模；（4）利用步骤（3）得到的误差源分析模型，对步骤（2）所建立的共线条件方程进行一阶泰勒展开，建立误差传递方程；（5）计算步骤（4）所建立的误差传递方程的外方位元素误差传递矩阵；（6）基于协方差传播律，建立由步骤（4）和步骤（5）确定的姿态位置参数协方差矩阵；（7）根据步骤（6）建立的协方差矩阵计算定位误差。

Description

一种航空高光谱成像系统定位误差分析方法

技术领域

本发明涉及一种航空高光谱成像系统定位误差分析方法，属于高光谱数据处理方法与应用技术领域，适用于航空高光谱遥感几何校正的理论方法和应用技术研究。

背景技术

在航空高光谱成像系统中，成像光谱仪由于受地形起伏、自身技术参数、工作模式以及载荷平台姿态位置信息等因素的影响，所采集的航空高光谱数据会产生较大的几何畸变。而带有几何畸变的遥感数据无法直接进行后续应用处理，因此航空高光谱成像系统的定位误差已成为评价高光谱数据获取质量的一个重要指标。

目前针对航空高光谱成像系统定位误差分析的方法主要包括三大类：基于几何成像关系的方法、基于共线条件方程全微分展开的方法和基于共线条件方程一阶泰勒展开的方法。基于几何成像关系的方法是通过构建地辅坐标系下像点与地面点的相似三角形关系，通过成像光谱仪切向和径向的标定误差求解地面点定位误差，该方法虽然直观反映了成像光谱仪技术参数对定位误差的影响，但是没有考虑姿态位置测量系统和载荷实时姿态位置信息对定位误差的影响；基于共线条件方程全微分展开的方法是通过对共线条件方程全微分展开，由像片三个内方位元素的精度、六个外方位元素的精度以及高程精度求解地面点定位误差，该方法较完整地反映了各误差源对定位误差的作用，但是未考虑各误差源相互作用下对定位误差的影响，也未能根据姿态位置测量系统与载荷的实际技术参数对误差源进行分析；基于共线条件方程一阶泰勒展开的方法是对共线条件方程一阶泰勒展开，应用统计学的方法由像片外方位元素求解定位误差，该方法依然未能根据姿态位置测量系统与载荷的实际技术参数对误差源进行分析。目前航空高光谱成像系统定位误差分析存在的问题是：未能根据姿态位置测量系统和载荷的实际技术参数对误差源进行分析，未能很好地反映姿态位置测量系统和载荷实时姿态位置信息对定位误差的影响。

发明内容

本发明的目的是针对现有定位误差分析方法未能根据POS（PositionandOrientationSystem，姿态位置测量系统）和载荷的实际技术参数对误差源进行分析以及未能很好地反映姿态位置测量系统和载荷实时姿态位置信息对定位误差的影响等不足，提出一种航空高光谱成像系统定位误差分析方法。

本发明的技术解决方案是：一种利用已知的POS技术参数、成像光谱仪技术参数，应用协方差传播律求解航空高光谱成像系统定位误差的方法。该方法主要是对共线条件方程一阶泰勒展开，通过计算各误差源的误差传递系数，利用协方差传播律分析定位误差。

本发明一种航空高光谱成像系统定位误差分析方法，其步骤如下：

（1）高光谱图像信息、姿态位置测量系统数据读入；

（2）根据几何成像物理关系建立共线条件方程；

（3）分析姿态位置测量系统测角误差、全球定位系统定位误差和成像光谱仪标定误差，对每个误差源进行建模；

（4）利用步骤（3）得到的误差源分析模型，对步骤（2）所建立的共线条件方程进行一阶泰勒展开，建立误差传递方程；

（5）计算步骤（4）所建立的误差传递方程的外方位元素误差传递矩阵；

（6）基于协方差传播律，建立由步骤（4）和步骤（5）确定的姿态位置参数协方差矩阵；

（7）根据步骤（6）建立的协方差矩阵计算定位误差。

其中，步骤（2）中所述的“根据几何成像物理过程建立共线条件方程”，其计算方法如下：

设投影中心S与地面点P在地辅坐标系D-X_tpY_tpZ_tp中坐标分别为(X_S,Y_S,Z_S)和(X_P,Y_P,Z_P)，则地面点P在像空间辅助坐标系D-XYZ中的坐标为(X_P-X_S,Y_P-Y_S,Z_P-Z_S)，而相应像点p在像空间辅助坐标系中的坐标为(X,Y,Z)。由于投影中心S、像点p和地面点P三点共线，根据相似三角形关系可得共线方程：

$\left.\begin{array}{c}x=-f\·\frac{a_1(X_P-X_S)+b_1(Y_P-Y_S)+c_1(Z_P-Z_S)}{a_3(X_P-X_S)+b_3(Y_P-Y_S)+c_3(Z_P-Z_S)}\\ y=-f\·\frac{a_2(X_P-X_S)+b_2(Y_P-Y_S)+c_2(Z_P-Z_S)}{a_3(X_P-X_S)+b_3(Y_P-Y_S)+c_3(Z_P-Z_S)}\end{array}\right\}$

将其逆算式表示为矩阵形式：

$\left[\begin{array}{c}{X}_p\\ Y_p\\ Z_p\end{array}\right]=\left[\begin{array}{c}{X}_s\\ Y_s\\ Z_s\end{array}\right]+\mathrm{\λ}\·R_i^m\·\left[\begin{array}{c}u\\ v\\ w\end{array}\right]$

式中，X_p,Y_p,Z_p表示待求的地面点在地辅坐标系中的位置参数；X_s,Y_s,Z_s表示当前扫描行中心像元在地辅坐标系中的位置参数；λ表示在地辅坐标系中像平面内任一点与地面对应点的比例关系；u,v,w表示探测器上当前待求像元在像空间坐标系中的位置参数；表示像空间坐标系（i系）到地辅坐标系（m系）的旋转矩阵，可以表示为其中表示地心坐标系（E系）到地辅坐标系（m系）的旋转矩阵，表示导航坐标系（g系）到地心坐标系（E系）的旋转矩阵，表示惯性测量单元坐标系（b系）到导航坐标系（g系）的旋转矩阵，表示传感器坐标系（c系）到惯性测量单元坐标系（b系）的旋转矩阵，表示像空间坐标系（i系）到传感器坐标系（c系）的旋转矩阵。

其中，步骤（3）所述的“分析姿态位置测量系统测角误差、全球定位系统定位误差和成像光谱仪标定误差，对每个误差源进行建模”，其计算方法如下：

姿态位置测量系统的测角误差由陀螺仪所固有的漂移误差产生，其模型为ε＝ε_b+ε_r+ε_g，ε为姿态位置测量系统的测角误差，ε_b＝[ε_bxε_byε_bz]^T为测角误差中的常值误差，ε_r＝[ε_rxε_ryε_rz]^T为测角误差中的漂移误差，ε_g为测角误差中的白噪声。在实际应用中，姿态位置测量系统在技术参数中给出测角精度，该精度的数值在计算过程中已经考虑了陀螺仪漂移误差的影响，因此可以直接应用测角精度进行计算，即：

式中，表示滚转轴、俯仰轴、航向轴的测角误差，表示三个轴所对应的常值误差，表示三个轴所对应的漂移误差，t表示姿态位置测量系统所记录的航摄时间；

GPS（GlobalPositionSystem，全球定位系统）的定位误差在短时间的航摄飞行中随航摄时间t线性变化，该定位误差可以采用下述线性模型来拟合：

$\left[\begin{array}{c}\mathrm{\Δ}{X}_S\\ \mathrm{\Δ}{Y}_S\\ \mathrm{\Δ}{Z}_S\end{array}\right]=\left[\begin{array}{c}{a}_X\\ a_Y\\ a_Z\end{array}\right]+\left[\begin{array}{c}{b}_X\\ b_Y\\ b_Z\end{array}\right]\·L$

式中，[ΔX_SΔY_SΔZ_S]^T表示地面点在地辅坐标系内由GPS接收机测量而产生的定位误差，[a_Xa_Ya_Z]^T表示GPS定位误差中的常值误差，[b_Xb_Yb_Z]^T表示GPS定位误差中的漂移误差，L表示t时间内所飞过的路程，设t时刻的飞行速度为v_t，Δt为前后两个时刻的间隔，t₀为航摄飞行的初始时刻，T为航摄飞行的终止时刻，则有

成像光谱仪的标定误差由实验室定标获得，该误差可以表示为：

ε_c＝[ΔuΔvΔw]^T

式中，Δu、Δv表示成像光谱仪探测器中当前像元的切向标定误差，Δw表示成像光谱仪探测器中当前像元的径向标定误差。

其中，步骤（4）中所述的“利用步骤（3）得到的误差源分析模型，对步骤（2）所建立的共线条件方程进行一阶泰勒展开，建立误差传递方程”，其计算方法如下：

对共线方程按照泰勒级数展开至一次项，建立误差传递方程：

式中， ${\left[\begin{array}{ccc}{V}_{X_p}& V_{Y_p}& V_{Z_p}\end{array}\right]}^T$ 表示待求地面点的定位误差， $\left[\begin{array}{c}\mathrm{\Δ}{X}_s\\ \mathrm{\Δ}{Y}_s\\ \mathrm{\Δ}{Z}_s\end{array}\right]=\left[\begin{array}{c}{a}_X\\ a_Y\\ a_Z\end{array}\right]+\left[\begin{array}{c}{b}_X\\ b_Y\\ b_Z\end{array}\right]\·L$ 表示当前扫描行中心像元的定位误差，由定位误差中的常值误差[a_Xa_Ya_Z]^T、漂移误差[b_Xb_Yb_Z]^T和飞行距离L确定；[ΔuΔvΔw]^T表示探测器上当前待求像元在像空间坐标系中的位置误差，即成像光谱仪标定误差；

表示POS测角误差，由POS技术参数直接给出；D₁表示[ΔX_sΔY_sΔZ_s]^T的误差传递系数矩阵，λ表示在地辅坐标系下像平面内任一点与地面对应点的比例，W₁表示[ΔuΔvΔw]^T的误差传递系数矩阵，W₂表示的误差传递系数矩阵。

其中，步骤（5）中所述的“计算步骤（4）所建立的误差传递方程的外方位元素误差传递矩阵”，其计算方法如下：

$D_1=E=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right];$

式中A₁,A₂,A₃,B₁,B₂,B₃,C₁,C₂,C₃分别表示矩阵W₁中各对应位置的元素；

式中ω、κ分别表示由姿态位置测量系统给出的载荷相对于惯性测量单元坐标系的滚转偏角、俯仰偏角以及航向偏角，表示导航坐标系（g系）到地辅坐标系（m系）的旋转矩阵，可表示为表示惯性测量单元坐标系（b系）到地辅坐标系（m系）的旋转矩阵，可表示为 表示像空间坐标系（i系）到导航坐标系（g系）的旋转矩阵，可表示为a₁,a₂,a₃,b₁,b₂,b₃,c₁,c₂,c₃分别表示矩阵W₂中各对应位置的元素；

式中，H为飞机与地面高程差，u,v,w表示探测器上当前待求像元在像空间坐标系中的位置参数，C₁,C₂,C₃由W₁求得。

其中，步骤（6）中所述的“基于协方差传播律，建立由步骤（4）和步骤（5）确定的位置参数协方差矩阵”，其计算方法如下：

按照协方差传播律，位置参数的协方差矩阵D_XYZ应满足以下条件：

D_XYZ＝K'·D_m·K'^T

式中，参数向量D_m为参数向量m的协方差矩阵 $D_m=\left[\begin{array}{ccc}{D}_P& 0_{3\×3}& 0_{3\×3}\\ 0_{3\×3}& D_G& 0_{3\×3}\\ 0_{3\×3}& 0_{3\×3}& D_T\end{array}\right],$ 其中 $D_G=\left[\begin{array}{ccc}{\mathrm{\δ}}_{X_s}^2& 0& 0\\ 0& {\mathrm{\δ}}_{Y_s}^2& 0\\ 0& 0& {\mathrm{\δ}}_{Z_s}^2\end{array}\right],$ $D_T=\left[\begin{array}{ccc}{\mathrm{\δ}}_u^2& 0& 0\\ 0& {\mathrm{\δ}}_v^2& 0\\ 0& 0& {\mathrm{\δ}}_w^2\end{array}\right];$ 可直接由POS和成像光谱仪技术参数求得；可由GPS定位误差求得，即 $\left[\begin{array}{c}{{\mathrm{\δ}}_X}_p\\ {{\mathrm{\δ}}_Y}_p\\ {{\mathrm{\δ}}_Z}_p\end{array}\right]=\left[\begin{array}{c}{a}_X\\ a_Y\\ a_Z\end{array}\right]+\left[\begin{array}{c}{b}_X\\ b_Y\\ b_Z\end{array}\right]\·L;$ K'为像片位置参数对向量m各参数的一阶偏导数矩阵，K'^T为其转置矩阵，则有：

$K^{\′}=\left[\begin{array}{ccc}\mathrm{\λ}\·W_2& E& \mathrm{\λ}\·W_2\end{array}\right]=\left[\begin{array}{ccccccccc}\mathrm{\λ}\·a_1& \mathrm{\λ}\·a_2& \mathrm{\λ}\·a_3& 1& 0& 0& \mathrm{\λ}\·A_1& \mathrm{\λ}\·A_2& \mathrm{\λ}\·A_3\\ \mathrm{\λ}\·b_1& \mathrm{\λ}\·b_2& \mathrm{\λ}\·b_3& 0& 1& 0& \mathrm{\λ}\·B_1& \mathrm{\λ}\·B_2& \mathrm{\λ}\·B_3\\ \mathrm{\λ}\·c_1& \mathrm{\λ}\·c_2& \mathrm{\λ}\·c_3& 0& 0& 1& \mathrm{\λ}\·C_1& \mathrm{\λ}\·C_2& \mathrm{\λ}\·C_3\end{array}\right]$

$K^{\′T}={\left[\begin{array}{ccc}\mathrm{\λ}\·W_2& E& \mathrm{\λ}\·W_2\end{array}\right]}^T={\left[\begin{array}{ccccccccc}\mathrm{\λ}\·a_1& \mathrm{\λ}\·a_2& \mathrm{\λ}\·a_3& 1& 0& 0& \mathrm{\λ}\·A_1& \mathrm{\λ}\·A_2& \mathrm{\λ}\·A_3\\ \mathrm{\λ}\·b_1& \mathrm{\λ}\·b_2& \mathrm{\λ}\·b_3& 0& 1& 0& \mathrm{\λ}\·B_1& \mathrm{\λ}\·B_2& \mathrm{\λ}\·B_3\\ \mathrm{\λ}\·c_1& \mathrm{\λ}\·c_2& \mathrm{\λ}\·c_3& 0& 0& 1& \mathrm{\λ}\·C_1& \mathrm{\λ}\·C_2& \mathrm{\λ}\·C_3\end{array}\right]}^T$

将上述参数带入到等式中计算可得：

$D_{\mathrm{XYZ}}=K^{\′}\·D_m\·K^{\′T}=\left[\begin{array}{ccc}{D}_{\mathrm{XX}}& D_{\mathrm{XY}}& D_{\mathrm{XZ}}\\ D_{\mathrm{YX}}& D_{\mathrm{YY}}& D_{\mathrm{YZ}}\\ D_{\mathrm{ZX}}& D_{\mathrm{ZY}}& D_{\mathrm{ZZ}}\end{array}\right]$

式中，D_XX、D_YY、D_ZZ分别为X方向、Y方向、Z方向单独作用下的定位误差，D_XY、D_YX为X方向与Y方向耦合作用下的定位误差，D_XZ、D_ZX为X方向与Z方向耦合作用下的定位误差，D_YZ、D_ZY为Y方向与Z方向耦合作用下的定位误差。

其中，步骤（7）中所述的“根据步骤（6）建立的协方差矩阵计算定位误差。”，其计算方法如下：

地辅坐标系下，地面点在XY平面内的定位误差为测区在XY平面内的定位中误差可表示为：

$\mathrm{\δ}=\sqrt{\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\Δ}}_i{\mathrm{\Δ}}_i}{n}}$

式中，δ为定位中误差，n为测区采样点个数。

本发明与现有技术相比的优点在于：克服了现有航空高光谱成像系统定位误差分析方法无法反映姿态位置测量系统和载荷实时姿态信息对定位误差的影响以及不能根据姿态位置测量系统与载荷的实际技术参数对误差源进行建模等不足，本方法对共线条件方程一阶泰勒展开，通过计算各误差源的误差传递系数，利用协方差传播律求解定位点误差。它具有以下优点：（1）在研究姿态位置测量系统和成像光谱仪技术参数特点的基础上，对航空高光谱成像系统定位误差源进行分析；（2）通过对共线条件方程一阶泰勒展开，应用协方差传播律进行定位误差分析，反映了姿态位置测量系统和载荷实时姿态位置信息对定位误差的影响。

（四）附图说明

图1为本发明涉及的一种航空高光谱成像系统定位误差分析方法的实现流程；

（五）具体实施方式

为了更好的说明本发明涉及的一种航空高光谱成像系统定位误差分析方法，利用山东遥感测区的PHI（PushbroomHyperspectralImagery）航空高光谱数据对定位误差进行分析。本发明一种航空高光谱成像系统定位误差分析方法，实现流程如图1所示，具体实现步骤如下：

（1）高光谱图像信息、姿态位置测量系统数据读入；

山东遥感测区的高光谱数据原始大小为512×10000，飞机飞行时POS姿态位置数据原始大小17×10000，即每行记录17个姿态位置信息参数，共记录了10000行。取其中1000行作为实验数据，读入数据大小为17×1000；

（2）根据几何成像物理过程建立共线条件方程；

设投影中心S与地面点P在地辅坐标系D-X_tpY_tpZ_tp中坐标分别为(X_S,Y_S,Z_S)和(X_P,Y_P,Z_P)，则地面点P在像空间辅助坐标系D-XYZ中的坐标为(X_P-X_S,Y_P-Y_S,Z_P-Z_S)，而相应像点p在像空间辅助坐标系中的坐标为(X,Y,Z)。由于投影中心S、像点p和地面点P三点共线，根据相似三角形关系可得共线方程：

$\left.\begin{array}{c}x=-f\·\frac{a_1(X_P-X_S)+b_1(Y_P-Y_S)+c_1(Z_P-Z_S)}{a_3(X_P-X_S)+b_3(Y_P-Y_S)+c_3(Z_P-Z_S)}\\ y=-f\·\frac{a_2(X_P-X_S)+b_2(Y_P-Y_S)+c_2(Z_P-Z_S)}{a_3(X_P-X_S)+b_3(Y_P-Y_S)+c_3(Z_P-Z_S)}\end{array}\right\}$

将其逆算式表示为矩阵形式：

$\left[\begin{array}{c}{X}_p\\ Y_p\\ Z_p\end{array}\right]=\left[\begin{array}{c}{X}_s\\ Y_s\\ Z_s\end{array}\right]+\mathrm{\λ}\·R_i^m\·\left[\begin{array}{c}u\\ v\\ w\end{array}\right]$

式中，X_p,Y_p,Z_p表示待求的地面点在地辅坐标系中的位置参数；X_s,Y_s,Z_s表示当前扫描行中心像元在地辅坐标系中的位置参数；λ表示在地辅坐标系中像平面内任一点与地面对应点的比例关系；u,v,w表示探测器上当前待求像元在像空间坐标系中的位置参数；表示像空间坐标系（i系）到地辅坐标系（m系）的旋转矩阵，可以表示为其中表示地心坐标系（E系）到地辅坐标系（m系）的旋转矩阵，表示导航坐标系（g系）到地心坐标系（E系）的旋转矩阵，表示惯性测量单元坐标系（b系）到导航坐标系（g系）的旋转矩阵，表示传感器坐标系（c系）到惯性测量单元坐标系（b系）的旋转矩阵，表示像空间坐标系（i系）到传感器坐标系（c系）的旋转矩阵。

（3）分析姿态位置测量系统测角误差、全球定位系统定位误差和成像光谱仪标定误差，对每个误差源进行建模；

姿态位置测量系统的测角误差由陀螺仪所固有的漂移误差产生，其模型为ε＝ε_b+ε_r+ε_g，ε为姿态位置测量系统的测角误差，ε_b＝[ε_bxε_byε_bz]^T为测角误差的常值误差，ε_r＝[ε_rxε_ryε_rz]^T为测角误差的漂移误差，ε_g为测角误差的白噪声。在实际应用中，POS在技术参数中给出测角精度，该精度的数值在计算过程中已经考虑了陀螺仪漂移误差的影响，因此可以直接应用测角精度进行计算，即：

式中，表示滚转轴、俯仰轴、航向轴的测角误差，表示三个轴所对应的常值误差，表示三个轴所对应的漂移误差；t表示姿态位置测量系统所记录的航摄时间；

GPS的定位误差在短时间的航摄飞行中随航摄时间t线性变化，该定位误差可以采用下述线性模型来拟合：

$\left[\begin{array}{c}\mathrm{\Δ}{X}_S\\ \mathrm{\Δ}{Y}_S\\ \mathrm{\Δ}{Z}_S\end{array}\right]=\left[\begin{array}{c}{a}_X\\ a_Y\\ a_Z\end{array}\right]+\left[\begin{array}{c}{b}_X\\ b_Y\\ b_Z\end{array}\right]\·L$

式中，[ΔX_SΔY_SΔZ_S]^T表示地面点在地辅坐标系内由GPS接收机测量而产生的定位误差，[a_Xa_Ya_Z]^T表示GPS定位误差中的常值误差，[b_Xb_Yb_Z]^T表示GPS定位误差中的漂移误差，L表示t时间内所飞过的路程，设t时刻的飞行速度为v_t，Δt为前后两个时刻的间隔，t₀为航摄飞行的初始时刻，T为航摄飞行的终止时刻，则有

成像光谱仪的标定误差由实验室定标获得，该误差可以表示为：

ε_c＝[ΔuΔvΔw]^T

式中，Δu、Δv表示成像光谱仪探测器中当前像元的切向标定误差，Δw表示成像光谱仪探测器中当前像元的径向标定误差；

（4）在步骤（3）的基础上，对步骤（2）所建立的共线条件方程进行一阶泰勒展开，建立误差传递方程；

对共线方程按照泰勒级数展开至一次项，建立误差传递方程：

式中， ${\left[\begin{array}{ccc}{V}_{X_p}& V_{Y_p}& V_{Z_p}\end{array}\right]}^T$ 表示待求地面点的定位误差， $\left[\begin{array}{c}\mathrm{\Δ}{X}_s\\ \mathrm{\Δ}{Y}_s\\ \mathrm{\Δ}{Z}_s\end{array}\right]=\left[\begin{array}{c}{a}_X\\ a_Y\\ a_Z\end{array}\right]+\left[\begin{array}{c}{b}_X\\ b_Y\\ b_Z\end{array}\right]\·L$ 表示当前扫描行中心像元的定位误差，由GPS接收机的常值误差[a_Xa_Ya_Z]^T、漂移误差[b_Xb_Yb_Z]^T和飞行距离L确定；[ΔuΔvΔw]^T表示探测器上当前待求像元在像空间坐标系中的位置误差，即成像光谱仪标定误差；

表示POS测角误差，由POS技术参数直接给出；D₁表示[ΔX_sΔY_sΔZ_s]^T的误差传递系数矩阵，λ表示在地辅坐标系下像平面内任一点与地面对应点的比例，W₁表示[ΔuΔvΔw]^T的误差传递系数矩阵，W₂表示的误差传递系数矩阵；

（5）计算步骤（4）所建立的误差传递方程的外方位元素误差传递矩阵；

$D_1=E=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right];$

式中A₁,A₂,A₃,B₁,B₂,B₃,C₁,C₂,C₃分别表示矩阵W₁中各对应位置的元素；

式中ω、κ分别表示由姿态位置测量系统给出的载荷相对于惯性测量单元坐标系的滚转偏角、俯仰偏角以及航向偏角，表示导航坐标系（g系）到地辅坐标系（m系）的旋转矩阵，可表示为表示惯性测量单元坐标系（b系）到地辅坐标系（m系）的旋转矩阵，可表示为表示像空间坐标系（i系）到导航坐标系（g系）的旋转矩阵，可表示为a₁,a₂,a₃,b₁,b₂,b₃,c₁,c₂,c₃分别表示矩阵W₂中各对应位置的元素；

式中，H为飞机与地面高程差，u,v,w表示探测器上当前待求像元在像空间坐标系中的位置参数，C₁,C₂,C₃由W₁求得；

（6）基于协方差传播律，建立由步骤（4）和步骤（5）确定的姿态位置参数协方差矩阵；

按照协方差传播律，位置参数的协方差矩阵D_XYZ应满足以下条件：

D_XYZ＝K'·D_m·K'^T

式中，参数向量D_m为参数向量m的协方差矩阵 $D_m=\left[\begin{array}{ccc}{D}_P& 0_{3\×3}& 0_{3\×3}\\ 0_{3\×3}& D_G& 0_{3\×3}\\ 0_{3\×3}& 0_{3\×3}& D_T\end{array}\right],$ 其中 $D_G=\left[\begin{array}{ccc}{\mathrm{\δ}}_{X_s}^2& 0& 0\\ 0& {\mathrm{\δ}}_{Y_s}^2& 0\\ 0& 0& {\mathrm{\δ}}_{Z_s}^2\end{array}\right],$

$D_T=\left[\begin{array}{ccc}{\mathrm{\δ}}_u^2& 0& 0\\ 0& {\mathrm{\δ}}_v^2& 0\\ 0& 0& {\mathrm{\δ}}_w^2\end{array}\right];$ 可直接由POS和成像光谱仪技术参数求得；可由GPS定位误差求得，即 $\left[\begin{array}{c}{{\mathrm{\δ}}_X}_p\\ {{\mathrm{\δ}}_Y}_p\\ {{\mathrm{\δ}}_Z}_p\end{array}\right]=\left[\begin{array}{c}{a}_X\\ a_Y\\ a_Z\end{array}\right]+\left[\begin{array}{c}{b}_X\\ b_Y\\ b_Z\end{array}\right]\·L;$ K'为像片位置参数对向量m各参数的一阶偏导数矩阵，K'^T为其转置矩阵，则有：

$K^{\′}=\left[\begin{array}{ccc}\mathrm{\λ}\·W_2& E& \mathrm{\λ}\·W_2\end{array}\right]=\left[\begin{array}{ccccccccc}\mathrm{\λ}\·a_1& \mathrm{\λ}\·a_2& \mathrm{\λ}\·a_3& 1& 0& 0& \mathrm{\λ}\·A_1& \mathrm{\λ}\·A_2& \mathrm{\λ}\·A_3\\ \mathrm{\λ}\·b_1& \mathrm{\λ}\·b_2& \mathrm{\λ}\·b_3& 0& 1& 0& \mathrm{\λ}\·B_1& \mathrm{\λ}\·B_2& \mathrm{\λ}\·B_3\\ \mathrm{\λ}\·c_1& \mathrm{\λ}\·c_2& \mathrm{\λ}\·c_3& 0& 0& 1& \mathrm{\λ}\·C_1& \mathrm{\λ}\·C_2& \mathrm{\λ}\·C_3\end{array}\right]$

$K^{\′T}={\left[\begin{array}{ccc}\mathrm{\λ}\·W_2& E& \mathrm{\λ}\·W_2\end{array}\right]}^T={\left[\begin{array}{ccccccccc}\mathrm{\λ}\·a_1& \mathrm{\λ}\·a_2& \mathrm{\λ}\·a_3& 1& 0& 0& \mathrm{\λ}\·A_1& \mathrm{\λ}\·A_2& \mathrm{\λ}\·A_3\\ \mathrm{\λ}\·b_1& \mathrm{\λ}\·b_2& \mathrm{\λ}\·b_3& 0& 1& 0& \mathrm{\λ}\·B_1& \mathrm{\λ}\·B_2& \mathrm{\λ}\·B_3\\ \mathrm{\λ}\·c_1& \mathrm{\λ}\·c_2& \mathrm{\λ}\·c_3& 0& 0& 1& \mathrm{\λ}\·C_1& \mathrm{\λ}\·C_2& \mathrm{\λ}\·C_3\end{array}\right]}^T$

将上述参数带入到等式中计算可得：

$D_{\mathrm{XYZ}}=K^{\′}\·D_m\·K^{\′T}=\left[\begin{array}{ccc}{D}_{\mathrm{XX}}& D_{\mathrm{XY}}& D_{\mathrm{XZ}}\\ D_{\mathrm{YX}}& D_{\mathrm{YY}}& D_{\mathrm{YZ}}\\ D_{\mathrm{ZX}}& D_{\mathrm{ZY}}& D_{\mathrm{ZZ}}\end{array}\right]$

式中，D_XX、D_YY、D_ZZ分别为X方向、Y方向、Z方向单独作用下的定位误差，D_XY、D_YX为X方向与Y方向耦合作用下的定位误差，D_XZ、D_ZX为X方向与Z方向耦合作用下的定位误差，D_YZ、D_ZY为Y方向与Z方向耦合作用下的定位误差；

（7）根据步骤（6）建立的协方差矩阵计算定位误差；

地辅坐标系下，地面点在XY平面内的定位误差为测区在XY平面内的定位中误差可表示为：

$\mathrm{\δ}=\sqrt{\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\Δ}}_i{\mathrm{\Δ}}_i}{n}}$

式中，δ为定位中误差，n=512×1000。

Claims (5)

1.一种航空高光谱成像系统定位误差分析方法，其特征在于：它包含以下步骤：

(1)高光谱图像信息、姿态位置测量系统数据读入；

(2)根据几何成像物理关系建立共线条件方程；

(3)分析姿态位置测量系统测角误差、全球定位系统定位误差和成像光谱仪标定误差，综合考虑陀螺仪的漂移误差、光学成像系统的畸变误差，根据物理成像过程、原理以及误差产生机制，对每个误差源进行建模，其方法如下：

姿态位置测量系统的测角误差由陀螺仪所固有的漂移误差产生，其模型为ε＝ε_b+ε_r+ε_g，ε为姿态位置测量系统的测角误差，ε_b＝[ε_bxε_byε_bz]^T为测角误差中的常值误差，ε_r＝[ε_rxε_ryε_rz]^T为测角误差中的漂移误差，ε_g为测角误差中的白噪声；在实际应用中，姿态位置测量系统在技术参数中给出测角精度，该精度的数值在计算过程中已经考虑了陀螺仪漂移误差的影响，因此直接应用测角精度进行计算，即：

式中，表示滚转轴、俯仰轴、航向轴的测角误差，表示三个轴所对应的常值误差，表示三个轴所对应的漂移误差，t表示姿态位置测量系统所记录的航摄时间；

全球定位系统的定位误差在短时间的航摄飞行中随航摄时间t线性变化，该定位误差采用下述线性模型进行拟合：

$\left[\begin{array}{c}{\mathrm{\ΔX}}_S\\ {\mathrm{\ΔY}}_S\\ {\mathrm{\ΔZ}}_S\end{array}\right]=\left[\begin{array}{c}{a}_X\\ a_Y\\ a_Z\end{array}\right]+\left[\begin{array}{c}{b}_X\\ b_Y\\ b_Z\end{array}\right]\·L$

式中，[ΔX_SΔY_SΔZ_S]^T表示地面点在地辅坐标系内由全球定位系统测量而产生的定位误差，[a_Xa_Ya_Z]^T表示定位误差中的常值误差，[b_Xb_Yb_Z]^T表示定位误差中的漂移误差，L表示t时间内所飞过的路程，设t_i时刻的飞行速度为v_ti，Δt为前后两个时刻的间隔，t₀为航摄飞行的初始时刻，T为航摄飞行的终止时刻，则有

成像光谱仪的标定误差由实验室定标获得，该误差表示为：

ε_c＝[ΔuΔvΔw]^T

式中，ε_c表示成像光谱仪标定误差，Δu、Δv表示成像光谱仪探测器中当前像元的切向标定误差，Δw表示成像光谱仪探测器中当前像元的径向标定误差；

(4)利用步骤(3)得到的误差源分析模型，对步骤(2)所建立的共线条件方程进行一阶泰勒展开，建立误差传递方程；

(5)计算步骤(4)所建立的误差传递方程的外方位元素误差传递矩阵；

(6)基于协方差传播律，建立由步骤(4)和步骤(5)确定的姿态位置参数协方差矩阵；

(7)根据步骤(6)建立的协方差矩阵计算定位误差，其具体方法如下：

地辅坐标系下，地面点在XY平面内的定位误差为其中D_XX、D_YY分别为X方向、Y方向单独作用下的定位误差，D_XY为X方向与Y方向耦合作用下的定位误差，测区在XY平面内的定位中误差可表示为：

$\mathrm{\δ}=\sqrt{\frac{\underset{i=1}{\overset{n}{\Σ}}{\mathrm{\Δ}}_i{\mathrm{\Δ}}_i}{n}}$

式中，δ为定位中误差，n为测区采样点个数。

2.根据权利要求1所述的一种航空高光谱成像系统定位误差分析方法，其特征在于：步骤(2)中所述的根据几何成像物理过程建立的共线条件方程为：

$\left[\begin{array}{c}{X}_p\\ Y_p\\ Z_p\end{array}\right]=\left[\begin{array}{c}{X}_s\\ Y_s\\ Z_s\end{array}\right]+\mathrm{\λ}\·R_i^m.\left[\begin{array}{c}u\\ v\\ w\end{array}\right]$

式中，X_p,Y_p,Z_p表示待求地面点在地辅坐标系中的位置参数；X_s,Y_s,Z_s表示当前扫描行中心像元在地辅坐标系中的位置参数；λ表示在地辅坐标系中像平面内任一点与地面对应点的比例关系；u,v,w表示探测器上当前待求像元在像空间坐标系中的位置参数；表示像空间坐标系(i系)到地辅坐标系(m系)的旋转矩阵，可以表示为其中表示地心坐标系(E系)到地辅坐标系(m系)的旋转矩阵，表示导航坐标系(g系)到地心坐标系(E系)的旋转矩阵，表示惯性测量单元坐标系(b系)到导航坐标系(g系)的旋转矩阵，表示传感器坐标系(c系)到惯性测量单元坐标系(b系)的旋转矩阵，表示像空间坐标系(i系)到传感器坐标系(c系)的旋转矩阵。

3.根据权利要求1所述的一种航空高光谱成像系统定位误差分析方法，其特征在于：其步骤(4)中所述的利用步骤(3)得到的误差源分析模型，对步骤(2)所建立的共线条件方程进行一阶泰勒展开，建立误差传递方程的方法如下：

式中， ${\left[\begin{array}{ccc}{V}_{X_p}& V_{Y_p}& V_{Z_p}\end{array}\right]}^T$ 表示待求地面点的定位误差， $\left[\begin{array}{c}\mathrm{\Δ}{X}_S\\ {\mathrm{\ΔY}}_S\\ {\mathrm{\ΔZ}}_S\end{array}\right]=\left[\begin{array}{c}{a}_X\\ a_Y\\ a_Z\end{array}\right]+\left[\begin{array}{c}{b}_X\\ b_Y\\ b_Z\end{array}\right]\·L$ 表示地面点在地辅坐标系内由全球定位系统测量而产生的定位误差，由定位误差中的常值误差[a_Xa_Ya_Z]^T、漂移误差[b_Xb_Yb_Z]^T和飞行路程L确定；[ΔuΔvΔw]^T表示探测器上当前待求像元在像空间坐标系中的位置误差，即成像光谱仪标定误差；表示姿态位置测量系统的测角误差，由姿态位置测量系统技术参数直接给出；D₁表示[ΔX_sΔY_sΔZ_s]^T的误差传递系数矩阵，λ表示在地辅坐标系下像平面内任一点与地面对应点的比例，W₁表示[ΔuΔvΔw]^T的误差传递系数矩阵，W₂表示的误差传递系数矩阵。

4.根据权利要求1所述的一种航空高光谱成像系统定位误差分析方法，其特征在于：其步骤(5)中所述的计算步骤(4)所建立的误差传递方程的外方位元素误差传递矩阵，其方法如下：

D₁表示[ΔX_sΔY_sΔZ_s]^T的误差传递系数矩阵，其中[ΔX_sΔY_sΔZ_s]^T为地面点在地辅坐标系内由全球定位系统测量而产生的定位误差；λ表示在地辅坐标系下像平面内任一点与地面对应点的比例；W₁表示[ΔuΔvΔw]^T的误差传递系数矩阵，其中[ΔuΔvΔw]^T为探测器上当前待求像元在像空间坐标系中的位置误差；W₂表示的误差传递系数矩阵，其中为姿态位置测量系统的测角误差；则D₁、W₁、W₂、λ的计算方法如下：

$D_1=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right];$

$W_1=R_i^m=R_E^m\·R_g^E\·R_b^g\·R_c^b\·R_i^c=\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],$ 式中A₁,A₂,A₃,B₁,B₂,B₃,C₁,C₂,C₃分别表示矩阵W₁中各对应位置的元素；

式中ω、κ分别表示由姿态位置测量系统给出的载荷相对于惯性测量单元坐标系的滚转偏角、俯仰偏角以及航向偏角，表示导航坐标系(g系)到地辅坐标系(m系)的旋转矩阵，可表示为表示惯性测量单元坐标系(b系)到地辅坐标系(m系)的旋转矩阵，可表示为表示像空间坐标系(i系)到导航坐标系(g系)的旋转矩阵，可表示为a₁,a₂,a₃,b₁,b₂,b₃,c₁,c₂,c₃分别表示矩阵W₂中各对应位置的元素；

式中H为飞机与地面高程差，u,v,w表示探测器上当前待求像元在像空间坐标系中的位置参数，C₁,C₂,C₃由W₁求得。

5.根据权利要求1所述的一种航空高光谱成像系统定位误差分析方法，其特征在于：步骤(6)中所述的基于协方差传播律，建立由步骤(4)和步骤(5)确定的位置参数协方差矩阵，其过程如下：

按照协方差传播律，位置参数的协方差矩阵D_XYZ应满足以下条件

D_XYZ＝K'·D_m·K^'T

其中，参数向量D_m为参数向量m的协方差矩阵，K'为像片位置参数对向量m各参数的一阶偏导数矩阵，K^'T为其转置矩阵，则有：

$D_{XYZ}=K^{\′}\·D_m\·K^{\′T}=\left[\begin{array}{ccc}{D}_{XY}& D_{XY}& D_{XZ}\\ D_{YX}& D_{YY}& D_{YZ}\\ D_{ZX}& D_{ZY}& D_{ZZ}\end{array}\right]$

式中，D_XX、D_YY、D_ZZ分别为X方向、Y方向、Z方向单独作用下的定位误差，D_XY、D_YX为X方向与Y方向耦合作用下的定位误差，D_XZ、D_ZX为X方向与Z方向耦合作用下的定位误差，D_YZ、D_ZY为Y方向与Z方向耦合作用下的定位误差。