CN105910607A - 基于地面控制的卫星长周期姿态误差修正方法 - Google Patents

Abstract

本发明公开了一种基于地面控制的卫星长周期姿态误差修正方法，主要包括三个步骤：步骤1，基于星敏感器光轴夹角计算长周期误差的频率和振幅初值；步骤2，利用地面控制点和外方位观测值计算姿态改正量；步骤3，利用步骤1计算的初值和步骤2的姿态改正量拟合出长周期误差。该方法通过对卫星长周期姿态的探测与修正可以进一步提高星敏陀螺联合处理的定姿精度，从而有助于提升测绘卫星的全球无控测图能力。

Description

基于地面控制的卫星长周期姿态误差修正方法

技术领域

本发明涉及卫星姿态后处理技术领域，特别涉及一种基于地面控制的卫星长周期误差修正方法，应用于高精度遥感测绘卫星绝对姿态确定。

背景技术

由于星敏感器与陀螺联合定姿具备精度高、可靠性好等优点，该组合模式已广泛应用于国内外高分辨率遥感测绘卫星中。组合模式中常用的星敏感器有德国Jena-Optronik的ASTRO系列、法国SODERN的SED系列等[Spot image，2005.Pre-processinglevels and location accuracy.Technical information.www.spotimage.com.；REDUCTION OF LOW FREQUENCY ERROR FOR SED36 AND APS BASED HYDRA STAR TRACKERS，Julien OUAKNINE(1)，Ludovic BLARRE(1)，Lionel ODDOS-MARCEL，Johan MONTEL(2)，Jean-Marc JULIO(2)，Proc.‘6th Internat.Conf.on Space Optics’，ESTEC，Noordwijk，The Netherlands，27-30June 2006(ESA SP-621，June 2006.]，其中，ASTRO 10主要应用于近地轨道卫星，包括SAR-Lupe，TerraSAR，DARPAs Orbital Express以及我国的HJ-1、FY-3、ZY-3等)。ASTRO 15较ASTRO 10视场基本不变，观星能力增强，单星精度提高。1997年开始研制的SED16应用于2001年5月首飞的SPOT5卫星。SED26是SED16的ITAR(国际军品贸易条例)自由版本，可提供三轴姿态和载体运动角速度信息。最新的SED36是专门为Pléiades卫星提供高姿态精度的星敏感器，对光学畸变进行了精确的校正，升级了星表，增加了导航星数目。

表1示出了几种典型星敏感器精度的比较。

表1

星敏感器的定姿误差分为三个部分，第一部分误差为安装偏差(Bias)，主要包括星敏感器内部变形误差、星敏感器相对安装误差(包括星敏之间、星敏陀螺之间)等，在卫星入轨后基本保持不变，这部分误差已经通过星敏感器在轨标定和星敏陀螺相互标定加以消除。[Ju，G.2001.Autonomous star sensing，pattern identification and attitudedetermination for spacecraft：An analytical and experimental study.Texas A&MUniversity：176-177.；Samaan，2003，Toward faster and more accurate star sensorsusing recursive centroiding and star identification，Texas A&M University，pp.25-29..郝雪涛，张广军，江洁，2005.星敏感器模型参数分析与校准方法研究[J]，光电工程，32(3)：5-8.；邢飞，董瑛，武延鹏，尤政，2005.星敏感器参数分析与自主校正[J]，清华大学学报(自然科学版).45(11)：1484-1488.；谢俊峰；龚健雅；江万寿，2009；一种改进的恒星相机在轨检校方法，测绘科学，v.34；No.158(02)122-124.；钱曾波，刘静宇，肖国超，1992.航天摄影测量[M]，解放军出版社：116～120.谢俊峰，江万寿，龚健雅.顾及星像点分布的恒星相机在轨检校[J]北京航空航天大学学报，2011，V37(10)：1271-1276.；申娟，张广军，魏新国，2010，基于卡尔曼滤波的星敏感器在轨校准方法，航空学报，31(6)：1221-1224.；K.Lai，J.Crassidis，and R.Harman，“In-space spacecraft alignmentcalibration using the unscented filter，”in Proceedings of the Guidance，Navigation，and Control Conference and Exhibit(AIAA’03)，Honolulu，Hawaii，USA，August 2003.]。

第二部分误差主要为测量噪声(Noise)，现在大多遥感卫星联合定姿采用星上实时扩展kalman滤波(Extended Kalman Filter，EKF)/(Unscented Kalman Filtering，UKF)或者事后地面离线处理加以抑制[Rtliter，A.H.J.de and C.J.Damaren，2002.ExtendedKalman Filtering and Nonlinear Predictive Filtering for Spacecratf AttiutdeDetermination.Canadian Aeronuatics and Space Journal 48(1)：13～23；MatthewC.Vandyke Jana L.Schwartz Christopher D.，2004.Hall Unscented Kalman Filteringfor Spacecraft Attitude State and Parameter Estimation，AAS-04-115.www.dept.aoe.vt.edu/～cdhall/papers/AAS04-115.pdf；Landsat 7 ImageAssessment System(IAS)Geometric Algorithm Theoretical Basis Document，Department of the Interior U.S.Geological Survey，LS-IAS-01，Version 1.0，http：//landsat.usgs.gov/documents/LS-IAS-01_Geometric_ATBD.pdf；TakanoriIwata，Tetsuo Kawahara，2007.Precision attitude determination for the advancedland observing satellite(ALOS)：Design，Verification，and on-orbitcalibration.]。这两类误差的抑制和消除方法已经成熟。

第三部分误差为低频误差(Low Fequency Error，LEF)，即本发明描述和需要消除的长周期误差，它主要是由于星敏感器在卫星在轨运行过程中，受轨道周期内温度变化影响，在较长一段时间内起伏变化较大，对星敏感器姿态测量精度产生周期性的影响。针对此类误差修正方法，国外一般从温控设计优化、镜头畸变标定、增加恒星星表数量等几个方面来改善，温控优化和完善星表等方法一般在软硬件设计阶段就已加以考虑，对改善此类误差能起到一定抑制作用。相机镜头畸变标定方法采用多项式拟合，仅能消除镜头畸变引起的大部分测量误差，无法完全消除由热环境引起的全部误差。国内熊凯等人所根据星敏感器特性，利用地面控制点计算已知参考矢量，加入到星敏陀螺联合kalman滤波模型中，从而标定星敏感器低频误差。该方法设计较为理想，忽略了卫星摄影测量过程中，由于存在地标像点人工量测或者自动匹配误差、相机内外标定误差、定轨误差等因素，地标点计算的参考矢量精度不高，反而给联合滤波引入干扰误差源，影响了姿态滤波精度。此外，国内有关学者针对星敏感器温度对姿态测量精度的影响，主要从星敏感器硬件上，分析了温度对器件精度影响，但是由于卫星及星敏感器热学特性的复杂性，没有给出温度与定姿精度的具体量化关系。

可见，现有温控设计优化、镜头畸变标定、增加恒星星表数量等方法改善后仍无法消除低频长周期姿态误差，同时现有的在星敏感器和陀螺卡尔曼滤波(KF)模型加入并消除低频长周期误差的方法由于存在地标像点人工量测或者自动匹配误差、相机内外标定误差、定轨误差等因素，地标点计算的参考矢量精度不高，反而给联合滤波引入干扰误差源，影响了姿态滤波精度，降低了卫星影像的有控制定位精度。

因此，需要一种更有效的卫星长周期姿态误差修正方法。

发明内容

为此，本发明提出一种基于地面控制的卫星长周期误差修正方法，可充分地消除由于现有技术的限制和缺陷导致的一个或多个问题。

本发明另外的优点、目的和特性，一部分将在下面的说明书中得到阐明，而另一部分对于本领域的普通技术人员通过对下面的说明的考察将是明显的或从本发明的实施中学到。通过在文字的说明书和权利要求书及附图中特别地指出的结构可实现和获得本发明目的和优点。

本发明提出了一种基于地面控制的卫星长周期姿态误差修正方法，其中，所述方法包括以下步骤：

步骤1，基于星敏感器光轴夹角计算长周期姿态误差的频率和振幅初值，其中，步骤1具体包括以下子步骤：

步骤1.1，随机抽取一段卫星下传的数据，得到一组离散的以归一化四元数形式表示的姿态数据，所述姿态数据为第一星敏感器和第二星敏感器测量的姿态数据，其中，所述第一星敏感器测量的姿态数据优选的用四元数表示为：

$q=\left[\begin{array}{c}\stackrel{\‾}{q}\\ q_4\end{array}\right]=\left[\begin{array}{c}{q}_1\\ q_2\\ q_3\\ q_4\end{array}\right],$

其中，为矢量，q₄为标量，且满足约束方程：

q₁ ²+q₂ ²+q₃ ²+q₄ ²＝1；

所述第二星敏感器测量的姿态数据用四元数表示为：

$q=\left[\begin{array}{c}{\stackrel{\‾}{q}}^{\′}\\ q_4^{\′}\end{array}\right]=\left[\begin{array}{c}{q}_1^{\′}\\ q_2^{\′}\\ q_3^{\′}\\ q_4^{\′}\end{array}\right],$

其中，为矢量，q′₄为标量，且满足约束方程：

$q_1^{\′2}+q_2^{\′2}+q_3^{\′2}+q_4^{\′2}=1.$

步骤1.2，计算所述第一星敏感器和第二星敏感器的光轴矢量，其中：

所述第一星敏感器在某一个时刻的光轴矢量A为：

A＝(2(q₁q₃+q₂q₄) 2(q₂q₃-q_1q4) -q₁ ²-q₂ ²+q₃ ²+q₄ ²)；

所述第二星敏感器在同一时刻的光轴矢量B为：

$B=\left(\begin{array}{ccc}2(q_1^{\′}{q}_3^{\′}+q_2^{\′}{q}_4^{\′})& 2(q_2^{\′}{q}_3^{\′}-q_1^{\′}{q}_4^{\′})& -q_1^{\′2}-q_2^{\′2}+q_3^{\′2}+q_4^{\′2}\end{array}\right);$

步骤1.3，基于所述第一星敏感器和第二星敏感器的光轴矢量，计算每一个时刻的所述第一星敏感器和第二星敏感器的光轴夹角实测值；

步骤1.4，基于所述第一星敏感器和第二星敏感器在卫星本体上的安装矢量，计算所述第一星敏感器和第二星敏感器的光轴夹角安装值；

步骤1.5，计算星敏感器的光轴夹角变化差值；

步骤1.6，计算星敏感器的光轴夹角变化差值的频率和振幅的初值，并将所述星敏感器的光轴夹角变化差值的频率和振幅的初值作为所述长周期误差的频率和振幅初值。

步骤2，利用地面控制点和外方位观测值计算姿态改正量；其中，步骤2具体包括以下子步骤：

步骤2.1，将相机坐标系转换到卫星本体坐标系；

步骤2.2，将卫星本体坐标系转换到J2000惯性坐标系；

步骤2.3，将J2000惯性坐标系转换到WGS-84坐标系；

步骤2.4，基于步骤2.1至步骤2.3的坐标系转换，构建附加了姿态改正量的严密成像模型；

步骤2.5，利用步骤2.4建立的严密成像模型，基于地面控制点和外方位观测值计算姿态改正量；

步骤3，利用步骤1计算的初值和步骤2计算的姿态改正量拟合出长周期误差，获得长周期误差修正模型；其中，步骤3具体包括以下步骤：

步骤3.1，针对步骤2解求的姿态改正量构建卫星姿态长周期误差模型：

其中，振幅a_{Δω_1}、a_{Δκ_1}，截距c_Δω、c_Δκ，频率ω_Δω、ω_Δκ为待求参数。

步骤3.2，基于最小二乘原理对振幅a_{Δω_1}、a_{Δκ_1}，截距c_Δω、c_Δκ，频率ω_Δω、Δ_Δκ进行解算，将解算得到的最终的参数值a_{Δω_1}，a_{Δκ_1}，ω_Δω，ω_Δκ，c_Δω，c_Δκ代入到所述卫星姿态长周期误差模型，即得到了长周期误差修正模型。

优选的，所述步骤1.3具体包括：

根据公式(2.1)计算光轴夹角实测值θ：

$\mathrm{\θ}=arctan\frac{A\×B}{A\·B}---\left(2.1\right)$

其中，A和B分别为所述第一星敏感器和第二星敏感器的光轴矢量；

根据公式(2.1)，依次计算每一个时刻的第一星敏感器和第二星敏感器的光轴夹角θ_t1，θ_t2，θ_t3，…，θ_tm，…，θ_tN，其中，tm表示第m个采样时刻，θ_tm表示计算的tm时刻的第一星敏感器和第二星敏感器的光轴夹角，tN表示第N个采样时刻，N表示光轴夹角的总数。

优选的，所述步骤1.4具体包括：

根据公式(3.1)计算所述第一星敏感器和第二星敏感器的光轴夹角安装值所述光轴夹角安装值为光轴夹角理论值；

$\stackrel{\‾}{\mathrm{\θ}}=arctan\frac{ST1\×ST2}{ST1\·ST2}---\left(1.3\right)$

其中，ST1和ST2分别为第一星敏感器和第二星敏感器在卫星本体上的安装矢量。

优选的，所述步骤1.5具体包括：

依次计算每一时刻的星敏感器光轴夹角实测值θ_t1，θ_t2，θ_t3，…，θ_tm，…，θ_tN与光轴安装值的差Δθ_t1，Δθ_t2，Δθ_t3，…，Δθ_tm，…，Δθ_tN：

${\mathrm{\Δ\θ}}_{tm}={\mathrm{\θ}}_{tm}-\stackrel{\‾}{\mathrm{\θ}}$

其中，tm表示第m个采样时刻，Δθ_tm表示计算的tm时刻的第一星敏感器和第二星敏感器的光轴夹角差值，N表示计算的光轴夹角差值的总数。优选的，所述步骤1.6具体包括：

将光轴夹角变化差值的主波形设为y＝f(t)＝L cos(w×(t-t₀))，其中L为振幅，w为频率，t为时间，t₀为初始时刻偏移；

计算星敏感器光轴夹角变化差值的主波形的振幅L的初值为：

其中，ξ＝max(Δθ_i)，为光轴夹角变化差值的最大值，为光轴夹角变化差值的最小值；

查找两个光轴夹角变化差值的最大值ξ₁和ξ₂的对应的时间间隔T₁和T₂，则主波的周期的初值为ΔT＝T₁-T₂，频率的初值为w＝2π/ΔT。

优选的，步骤2.4具体包括：

设待求量三轴姿态改正量为(Δω_t，Δκ_t)，则基于姿态改正角计算姿态改正矩阵R_C(t)为：

$R_C\left(t\right)=\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₃由Δω_t、和Δκ_t计算，公式如下：

将姿态改正矩阵R_C(t)代入到严密成像模型中，则可以得到附加了姿态改正量的严密成像模型：

${\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_{WGS-84}={\left[\begin{array}{c}{X}_{GPS}\\ Y_{GPS}\\ Z_{GPS}\end{array}\right]}_{WGS-84}+{\mathrm{mR}}_C\left(t\right)R_{J2000}^{WGS-84}{R}_{Star}^{J2000}{R}_{body}^{star}{R}_{camera}^{body}\left[\begin{array}{c}tan\left({\mathrm{\ψ}}_X\right)\\ tan\left({\mathrm{\ψ}}_Y\right)\\ -1\end{array}\right]$

式中，为WGS-84坐标系下地面点的坐标，其中X，Y，Z分别为地面点的横坐标、纵坐标和高程值，为星上GPS装置测得的卫星质心的位置，X_GPS、Y_GPS、Z_GPS分别表示由星上载荷GPS接收机测量的卫星质心的横坐标，纵坐标和高度值，m为比例因子，其由卫星飞行高度决定，R_C(t)为t时刻的姿态改正矩阵，为J2000惯性坐标系转换到WGS-84坐标系的旋转矩阵，为卫星本体坐标系转换到J2000惯性坐标系的旋转矩阵，为相机坐标系转换到卫星本体坐标系的旋转矩阵，为探元指向角。

优选的，步骤2.5具体包括以下步骤：

步骤2.5.1，设置定向片个数，并计算控制点范围内的定向片成像时刻和定向片成像时刻之间的间隔；

步骤2.5.2，将控制点的经纬度坐标(B，L，H)转换为地心直角坐标(X，Y，Z)；

步骤2.5.3，获得定向片时刻处姿态值，并设定定向片姿态改正数初值；

步骤2.5.4，计算各控制点在立体像对的姿态和轨道；

步骤2.5.5，获取控制点所在位置的CCD像元对应的指向角；

步骤2.5.6，将由卫星方提供的传感器在卫星本体上的安装矩阵代入到下述公式(16)中，并在计算中作为已知且不变，设定定向片姿态改正量初值为0：

${\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_{WGS-84}={\left[\begin{array}{c}{X}_{GPS}\\ Y_{GPS}\\ Z_{GPS}\end{array}\right]}_{WGS-84}+{\mathrm{mR}}_C\left(t\right)R_{J2000}^{WGS-84}{R}_{Star}^{J2000}{R}_{body}^{star}{R}_{camera}^{body}\left[\begin{array}{c}tan\left({\mathrm{\ψ}}_X\right)\\ tan\left({\mathrm{\ψ}}_Y\right)\\ -1\end{array}\right]$

步骤2.5.7，利用公式(16)，对定向片姿态改正参数进行线性化，得到像控点对应的姿态改正数的误差方程；

步骤2.5.8，利用最小二乘解算未知数，根据下视和前后视分辨率不同以及对平差结果的影响，给定不同的权，对误差方程进行法化，解算定向片的姿态改正量。

优选的，所述步骤3.2具体包括：

步骤(3.2)具体为：

根据步骤(3.1)中的卫星姿态长周期误差模型，列出误差方程形式：

其中，

f_{ω_i}，f_{κ_i}表示ω，k分别对应的第i个测量值，Δa_{Δω_1}，Δa_{Δθ_1}，Δω_{Δω_1}，Δω_{Δκ_1}，Δc_{Δω_1}，Δc_{Δκ_1}表示待求的未知参数，为振幅、频率和截距的初值，其中，振幅、频率的初值为在步骤1中计算获得的振幅和频率的初值，截距c_Δω、c_Δκ的初值设为0。

根据最小二乘原理，最小二乘的解为：

将各参数初值代入上式，根据最小二乘计算求得未知参数改正量Δa_{Δω_1}，Δa_{Δκ_1}，Δω_Δω，Δω_Δκ，Δc_Δω，Δc_Δκ，并计算出对应的残差V_Δω，V_Δκ，若结果残差中误差小于给定阈值K＝0.00001度，则计算结束；否则累加未知参数值Δa_{Δω_1}，Δa_{Δκ_1}，Δω_Δω，Δω_Δκ，Δc_Δω，Δc_Δκ得到新的参数a_{Δω_1}，a_{Δκ_1}，ω_Δω，ω_Δκ，c_Δω，c_Δκ，重新代入上式中解算。迭代结束后，得到最终的参数值a_{Δω_1}，a_{Δκ_1}，ω_Δω，ω_Δκ，c_Δω，c_Δκ，代入到卫星姿态长周期误差模型，即得到了长周期误差修正模型。

本发明可以提高卫星绝对姿态精度，拓展卫星全球无控制测图能力，在完成精密定轨、相机内外方位标定的基础上，采用地面控制实现对卫星长周期误差修正。由于姿态长周期系统误差在短时间内变化不大，将其上行注入到星上或者在地面定姿系统中进行补偿，即可以应用境内控制点的标定结果来提高卫星在境外的定姿精度，进而提高无地面控制点情况下的地面定位精轴夹角变化量随时间变化图。

附图说明

图1为根据本发明实施例的基于两个星敏感器光轴夹角变化探测长周期姿态误差的流程图；

图2为根据本发明实施例的传感器坐标系和卫星本体坐标系示意图；

图3为根据本发明实施例的空间固定参考坐标系J2000示意图；

图4为根据本发明实施例的地球固定地面参考系WGS84示意图；

图5为根据本发明实施例的两个固联安装的星敏感器光轴在卫星本体坐标系的指向；

图6为根据本发明实施例的基于地面控制的姿态角解算流程；

图7为根据本发明实施例的模型选取与参数估计方案；

图8为根据本发明实施例的沿CCD方向指向角内插示意图；

图9为根据本发明实施例的解算定向片的姿态改正量的示意图；

图10为根据本发明实施例的最小二乘拟合姿态迭代过程；

图11为根据本发明实施例的两个星敏光轴夹角随时间变化图；

图12为根据本发明实施例的两个星敏光轴夹角变化量随时间变化图；

图13为根据本发明实施例的两个星敏感器光轴夹角变化的初步拟合结果。

具体实施方式

下面参照附图对本发明进行更全面的描述，其中说明本发明的示例性实施例。

本发明提出的基于地面控制的卫星长周期误差修正方法包括以下步骤：

步骤1，基于星敏感器光轴夹角计算长周期误差的频率和振幅初值；

步骤2，利用地面控制点和外方位观测值计算姿态改正量；

步骤3，利用步骤1计算的初值和步骤2的姿态改正量拟合出长周期误差。

下面结合附图对每个步骤作详细说明。

如图1所示，步骤1包括以下子步骤：

步骤1.1，随机抽取一段卫星下传的数据，得到一组离散的以归一化四元数形式表示的姿态数据，所述姿态数据为第一星敏感器和第二星敏感器测量的姿态数据。

根据本发明的实施例，可以随机抽取卫星原始姿态数据；优选的，本发明的实施案例以资源三号卫星下传的原始姿态数据为例，但并不限于此。随机抽取一段卫星下传的原始姿态数据，该段数据获取时间为2012年1月14日，持续时长为181分钟左右，采样频率1Hz，离散的以四元数形式表示的姿态数据值。数据包含10906帧双星敏感器(即，星敏感器1和星敏感器2)的原始测量数据，双星敏感器在卫星平台上固联安装，两个星敏感器的光轴呈一定夹角，相对关系理论上认为不变。

第一星敏感器测量的姿态数据用四元数表示为：

$q=\left[\begin{array}{c}\stackrel{\‾}{q}\\ q_4\end{array}\right]=\left[\begin{array}{c}{q}_1\\ q_2\\ q_3\\ q_4\end{array}\right]---\left(1\right)$

其中，为矢量，一组q₁、q₂和q₃值描述了卫星的一种姿态，q₄为标量，且满足约束方程：

q₁ ²+q₂ ²+q₃ ²+q₄ ²＝1 (2)

同理，第二星敏感器测量的姿态数据用四元数表示为：

$q=\left[\begin{array}{c}{\stackrel{\‾}{q}}^{\′}\\ q_4^{\′}\end{array}\right]=\left[\begin{array}{c}{q}_1^{\′}\\ q_2^{\′}\\ q_3^{\′}\\ q_4^{\′}\end{array}\right],$

其中，为矢量，一组q′₁、q′₂、和q′₃、值描述了卫星的一种姿态，q′₄为标量，且满足约束方程：

$q_1^{\′2}+q_2^{\′2}+q_3^{\′2}+q_4^{\′2}=1.$

双星敏感器采样时间同步，以第一帧下传数据为例：(94638601.65112 -0.383423537 -0.470189244 -0.431782395 94638601.65112 0.123371840 -0.065405965-0.755955458)。其中，第1个数为星敏感器1的时标，第2-4个数为该时刻星敏感器1测量的J2000惯性坐标系下姿态q值3个矢量，第5个数为星敏感器2的时标，第6-8为星敏感器2测量的J2000惯性坐标系下姿态q值3个矢量。

从而，第一星敏感器的测量的姿态数据可以表示为：

$\stackrel{\‾}{q}=\[\begin{array}{ccc}{q}_1& q_2& q_3\end{array}\]={\left[\begin{array}{ccc}-0.383423537& -0.470189244& -0.431782395\end{array}\right]}^T,$ q₄＝0.66743721。

第二星敏感器的测量的姿态数据可以表示为：

${\stackrel{\‾}{q}}^{\′}={\left[\begin{array}{ccc}{q}_1^{\′}& q_2^{\′}& q_3^{\′}\end{array}\right]}^T={\left[\begin{array}{ccc}0.12337184& -0.065405965& -0.755955458\end{array}\right]}^T,$ q₄′＝0.63957798。

步骤1.2，计算所述第一星敏感器和第二星敏感器的光轴矢量；

将观测的姿态数据以四元数表示，在直角坐标系下，惯性系与星敏感器本体坐标系下姿态转换关系如下式：

$\left(\begin{array}{c}{x}_b\\ y_b\\ z_b\end{array}\right)=R\left(q\right)\left(\begin{array}{c}{x}_i\\ y_i\\ z_i\end{array}\right)---\left(3\right)$

其中，为惯性坐标系下姿态，为星敏感器本体坐标系下姿态。R(q)为旋转矩阵，且R(q)为：

$R\left(q\right)=\left[\begin{array}{ccc}{q_1}^2-{q_2}^2-{q_3}^2+{q_4}^2& 2(q_1{q}_2+q_3{q}_4)& 2(q_1{q}_3-q_2{q}_4)\\ 2(q_1{q}_2-q_3{q}_4)& -{q_1}^2+{q_2}^2-{q_3}^2+{q_4}^2& 2(q_2{q}_3+q_1{q}_4)\\ 2(q_1{q}_3+q_2{q}_4)& 2(q_2{q}_3-q_1{q}_4)& -{q_1}^2-{q_2}^2+{q_3}^2+{q_4}^2\end{array}\right].$

星敏感器本体坐标系的Zb轴即为星敏感器光轴，Zb轴在惯性坐标系中的坐标为：

(2(q₁q₃+q₂q₄) 2(q₂q₃-q₁q₄) -q₁ ²-q₂ ²+q₃ ²+q₄ ²)。

下面对第一星敏感器和第二星敏感器的光轴矢量分别说明：

第一星敏感器的旋转矩阵为：

$R\left(q\right)=\left[\begin{array}{ccc}{q_1}^2-{q_2}^2-{q_3}^2+{q_4}^2& 2(q_1{q}_2+q_3{q}_4)& 2(q_1{q}_3-q_2{q}_4)\\ 2(q_1{q}_2-q_3{q}_4)& -{q_1}^2+{q_2}^2-{q_3}^2+{q_4}^2& 2(q_2{q}_3+q_1{q}_4)\\ 2(q_1{q}_3+q_2{q}_4)& 2(q_2{q}_3-q_1{q}_4)& -{q_1}^2-{q_2}^2+{q_3}^2+{q_4}^2\end{array}\right]$

第二星敏感器的旋转矩阵为：

$R\left(q^{\′}\right)=\left[\begin{array}{ccc}{q}_1^{\′2}-q_2^{\′2}-q_3^{\′2}+q_4^{\′2}& 2(q_1^{\′}{q}_2^{\′}+q_3^{\′}{q}_4^{\′})& 2(q_1^{\′}{q}_3^{\′}-q_2^{\′}{q}_4^{\′})\\ 2(q_1^{\′}{q}_2^{\′}-q_3^{\′}{q}_4^{\′})& -q_1^{\′2}+q_2^{\′2}-q_3^{\′2}+q_4^{\′2}& 2(q_2^{\′}{q}_3^{\′}+q_1^{\′}{q}_4^{\′})\\ 2(q_1^{\′}{q}_3^{\′}-q_2^{\′}{q}_4^{\′})& 2(q_2^{\′}{q}_3^{\′}-q_1^{\′}{q}_4^{\′})& -q_1^{\′2}-q_2^{\′2}+q_3^{\′2}+q_4^{\′2}\end{array}\right]$

第一星敏感器在当前时刻的光轴指向(即光轴的方向)用矢量A表示：

A＝(2(q₁q₃+q₂q₄) 2(q₂q₃-q₁q₄) -q₁ ²-q₂ ²+q₃ ²+q₄ ²)＝(-0.296532528，0.917861147，0.263816932)；

第二星敏感器在同一时刻的光轴指向用矢量B表示：

$\begin{array}{c}B=\left(\begin{array}{ccc}2\left(q_1^{\′}{q}_3^{\′}+q_2^{\′}{q}_4^{\′}\right)& 2\left(q_2^{\′}{q}_3^{\′}-q_1^{\′}{q}_4^{\′}\right)& -q_1^{\′2}-q_2^{\′2}+q_3^{\′2}+q_4^{\′2}\end{array}\right)=(-0.270188880,\\ \begin{array}{cc}-0.058918585,& 0.961002897\end{array}).\end{array}$

(1.3)基于所述第一星敏感器和第二星敏感器的光轴矢量，计算所述第一星敏感器和第二星敏感器的光轴夹角实测值；

两星敏感器的安装关系如图5所示。

根据两个矢量的点积公式如下：

A·B＝||A|| ||B||cosθ (4)

其中θ为光轴矢量A，B的夹角。

根据两个矢量的叉积公式如下：

A×B＝||A|| ||B||sinθ (5)

对叉积取模然后除以点积：

$tan\mathrm{\θ}=\frac{A\×B}{A\·B}---\left(6\right)$

所以：

$\mathrm{\θ}=arctan\frac{A\×B}{A\·B}---\left(7\right)$

其中，A和B分别为所述第一星敏感器和第二星敏感器的光轴矢量，×和·分别表示向量的叉乘和点乘。

根据两个星敏感器时标，内插四元数到统一时刻t下，由于两个星敏为同时曝光，不存在时标统一问题。步骤1.2中的矢量A和B的夹角可以根据公式(7)计算：

重复步骤(1.2)、(1.3)，依次计算每一个时刻的第一星敏感器和第二星敏感器的光轴夹角实测值θ_t1，θ_t2，θ_t3，…，θ_tm，…，θ_tN，其中，tm表示第m个采样时刻，θ_tm表示计算的tm时刻的第一星敏感器和第二星敏感器的光轴夹角，tN表示第N个采样时刻，N表示光轴夹角的总数。

如图11所示，x轴为采样时刻，从2009年1月1日12：59：59为起始点累计到当前时间，单位为秒，y轴表示计算的两个星敏感器的光轴夹角实测值。

(1.4)基于所述第一星敏感器和第二星敏感器在卫星本体上的安装矢量，计算所述第一星敏感器和第二星敏感器的光轴夹角安装值θ，即光轴夹角理论值；

如前所述，两个星敏感器是固联安装，光轴的理论夹角值是不变的，本发明充分利用这一特点，以两个星敏感器为参考，计算并判断基于星敏感器光轴测量值计算的夹角变化情况，如果发生变化，则说明是由姿态测量误差引入的，从而判断姿态误差变化情况，有效解决了由于星上无静态参照物无法探测姿态系统误差变化规律的问题。

由于星敏感器的安装固连在卫星平台，两个星敏感器的实际光轴夹角认为不变，设第一星敏感器和第二星敏感器在卫星本体上的安装矢量分别为ST1、ST2，该值在卫星发射前由相关研制单位在地面实验室测量并提供给用户，如图7所示，星敏感光轴在卫星本体坐标系内的关系为：第一星敏感器(ASTRO 10)光轴ZST2与OXZ面夹角为30°，与+Y轴夹角为60°，在OXZ平面的投影与星体-X轴夹角为70°，与-Z轴的夹角为20°。第二星敏感器(ASTRO10)光轴ZST3与OXZ面夹角为25°，也+Y轴夹角为65°，在OYZ平面的投影与星体+X轴夹角为25°，与-Z轴的夹角为65°。

则ST1、ST2在卫星平台的表示为：

ST1＝[-cos(30°)cos(70°) sin(30°) -cos(30°)sin(70°)]；

ST2＝[cos(25°)cos25° sin(25°) -cos(25°)sin(25°)]；

则第一星敏感器和第二星敏感器的光轴安装值(即光轴夹角理论值)为：

$\stackrel{\‾}{\mathrm{\θ}}=arctan\frac{ST1\×ST2}{ST1\·ST2}---\left(8\right)$

其中，ST1和ST2分别为第一星敏感器和第二星敏感器在卫星本体上的安装矢量。

(1.5)计算星敏感器的光轴夹角变化差值；

具体的，依次计算每一时刻的星敏感器光轴夹角实测值θ_t1，θ_t2，θ_t3，…，θ_tm，…，θ_tN与光轴安装值(光轴夹角理论值)的差Δθ_t1，Δθ_t2，Δθ_t3，…，Δθ_tm，…，Δθ_tN：

${\mathrm{\Δ\θ}}_{tm}={\mathrm{\θ}}_{tm}-\stackrel{\‾}{\mathrm{\θ}}---\left(9\right)$

其中，tm表示第m个采样时刻，Δθ_tm表示计算的tm时刻的第一星敏感器和第二星敏感器的光轴夹角差值，N表示计算的光轴夹角差值的总数；

计算结果如图12所示，x轴为采样时刻，从2009年1月1日12：59：59为起始点累计到当前时间，单位为秒，y轴为光轴夹角。

(1.6)计算星敏感器的光轴夹角变化差值的频率和振幅的初值。所述步骤1.6具体包括：

将光轴夹角变化差值的主波形设为y＝f(t)＝L cos(w×(t-t₀))，其中L为振幅，w为频率，t为时间，t₀为初始时刻偏移；

计算星敏感器光轴夹角变化差值的主波形的振幅L的初值为：

其中，ξ＝max(Δθ_i)，为光轴夹角变化差值的最大值，为光轴夹角变化差值的最小值；

查找两个光轴夹角变化差值的最大值ξ₁和ξ₂的对应的时间间隔T₁和T₂，则主波的周期的初值为ΔT＝T₁-T₂，频率的初值为w＝2π/ΔT。

如图12所示，可以看出波形变化，将光轴夹角变化差值的主波形设为y＝f(t)＝Lcos(w×(t-t₀))，如图13所示，其中L为振幅，w为频率，t为时间，t₀为初始时刻偏移。计算主波形的主要参数包括振幅和频率等，为步骤3提供初值。

首先计算星敏感器光轴夹角变化差值的主波形的振幅初值，选取略大于一个信号周期的数据，计算变化量的最大值ξ＝max(Δθ_i)＝10角秒，计算变化差值的最小值则根据波形特征，星敏感器光轴夹角变化差值的振幅的初值为：然后计算光轴夹角变化差值的频率初值，查找两个变化最大值ξ₁和ξ₂的对应的时间95783638.5和95777966.69(单位：秒)，则主波的周期为ΔT＝T₁-T₂＝5671.81秒，w＝2π/ΔT＝0.0011Hz。

由于星敏感器姿态可以表示成光轴和旋转角，所以光轴的夹角变化可以反映出三轴姿态也应呈现出相似的变化规律，比如二者的波形相似，频率相近，但振幅不等。步骤(1.6)计算的长周期误差的频率和振幅初值即为长周期误差的频率和振幅初值，利用光轴夹角差值变化规律的探测解决了姿态无静态参考量无法探测的问题，为步骤3的计算提供初值。

所述的步骤(1)主要采用下传全弧段原始星敏感器、陀螺等姿态测量数据。利用两个或多个星敏感器安装关系几乎不变的特性，采用两个星敏感器光轴夹角来探测星敏感器姿态测量中存在的长周期变化探测，解决了由于准确参考导致星上姿态系统误差变化无法探测的问题。

步骤2，利用地面控制点和外方位观测值计算姿态改正量。

所述步骤(2)针对长周期姿态误差随时间变化规律特性，从精确重构误差模型的角度，设计地面控制数据分布，确定启算景间距、数量以及启算景内地面控制点分布等，利用加入姿态变化项的严密成型模型，反算获取近乎等纬度间距离散姿态误差项。

从查阅到文献情况看，现有研究工作一般假定卫星在同一时刻可获得2个以上已知地面控制点信息，利用双矢量(多矢量)定姿反推得出卫星姿态，与卫星确定姿态对比获得系统误差信息。从原理看该方法只适用于相机为框幅式成像原理的卫星，对目前在航天遥感领域通常使用的线阵推扫式遥感相机并不适用，因为线阵式相机为行扫描成像，很难保证在卫星同一时刻获得2个以上地面点信息。针对这一问题，该步骤采用定向片平差方法，本发明从线阵相机的成像原理出发，以星敏感器系统误差为未知量，以卫星到已知地面点视线矢量为观测量直接建立观测方程，多个已知地面点可以建立多个观测方程，由最小二乘解算得到系统误差信息。该方法可以有效避开利用双矢量(多矢量)方法求解卫星姿态，因此并不要求卫星在同一时刻同时获得多个地面点信息。

步骤2主要包括以下子步骤：

(2.1)将相机坐标系转换到卫星本体坐标系；

当传感器线阵列以其卫星轨道方向为轴向倾斜Ψ_y角、传感器阵列在其卫星轨道面内以垂直轨道方向为轴倾斜Ψ_x角时

传感器在卫星本体上的安装矩阵为线阵相机CCD每一个探元在相机坐标系中有一个三维指向，如图2所示，它可以通过实验室给定的相机主距、主点，CCD尺寸大小和个数等进行计算，也可经过在轨几何检校重新标定，提供给每个CCD的指向角文件，可供用户几何定位使用。指向角的表达式如下：

${\stackrel{\→}{u}}_i^{\′}=\left(\begin{array}{c}tg\left({\mathrm{\ψ}}_X\right)\\ tg\left({\mathrm{\ψ}}_Y\right)\\ -1\end{array}\right)---\left(11\right)$

在严密成像模型中进行归一化得：

${\stackrel{\→}{u}}_i=\frac{{\stackrel{\→}{u}}_i^{\′}}{\left|\right|{\stackrel{\→}{u}}_i^{\′}\left|\right|}---\left(12\right)$

线阵上每个CCD均给出对应的指向角，第i个CCD的指向角为

传感器坐标系的定义：传感器坐标系的原点在线阵投影中心，沿扫描线方向为X轴，沿飞行方向为Y轴，Z轴按照右手法则确定。

卫星本体坐标系的定义：本体坐标系的原点位于卫星的质心，X轴、Y轴和Z轴分别为卫星的三个主惯量轴。X轴沿着卫星横轴，Y轴沿着纵轴指向卫星飞行方向，Z轴按照右手法则确定。本体坐标系是卫星上仪器设备安装的基准参考坐标系。

(2.2)将卫星本体坐标系转换到J2000惯性坐标系；

在成像时刻t，星敏感器本体坐标系到J2000惯性坐标系的旋转矩阵为：

$R_{star}^{J2000}\left(t\right)=\left[\begin{array}{ccc}{q_1}^2-{q_2}^2-{q_3}^2+{q_4}^2& 2(q_1{q}_2+q_3{q}_4)& 2(q_1{q}_3-q_2{q}_4)\\ 2(q_1{q}_2-q_3{q}_4)& -{q_1}^2+{q_2}^2-{q_3}^2+{q_4}^2& 2(q_2{q}_3+q_1{q}_4)\\ 2(q_1{q}_3+q_2{q}_4)& 2(q_2{q}_3-q_1{q}_4)& -{q_1}^2-{q_2}^2+{q_3}^2+{q_4}^2\end{array}\right]---\left(13\right)$

其中，q1，q2，q3，q4分别表示t时刻星敏感器在J2000惯性坐标系下的姿态四元数，具体如公式(1)(2)所描述。

本体到J2000惯性坐标系的转换：

$R_{body}^{J2000}=R_{star}^{J2000}{R}_{body}^{star}---\left(14\right)$

其中，星敏感器在卫星本体系下的安装矩阵为卫星发射前实验室测定并提供给用户，为已知。

空间固定惯性J2000参考系定义：如图3所示，空间固定惯性参考系简称CIS(Conventional Inertial System，CIS)，常用来描述卫星的运动，一般卫星星历的计算都在该坐标系下完成。CIS的原点为地球质心，Z轴指向天球的北极，X轴指向春分点，Y轴按照右手法则确定。由于地球围绕太阳运动，春分点和北极点经常发生变化。因此，国际组织便规定以某个时刻的春分点和北极点为基准，建立协议空间固定惯性系统。CIS一般采用国际大地测量协会和国际天文学联合会会议于1984年启用的协议天球坐标系J2000。

(2.3)将J2000惯性坐标系转换到WGS-84坐标系。

${\left[\begin{array}{c}X-X_S\\ Y-Y_S\\ Z-Z_S\end{array}\right]}_{WGS-84}=R_{J2000}^{WGS-84}{\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_{J2000}=W\left(t\right)\·R\left(t\right)\·PN\left(t\right){\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_{J2000}---\left(15\right)$

其中PN(t)为岁差和章动矩阵，R(t)为地球自转矩阵，W(t)为极移矩阵。具体形式参考IERS2000。

通过步骤2.1、步骤2.2和步骤2.3的铺垫，为步骤2.4构建附加了姿态校正的严密成像模型奠定了基础。

(2.4)基于步骤(2.1)至步骤(2.3)的坐标系转换，构建附加了姿态改正量的严密成像模型；

设待求量三轴姿态改正量为(Δω_t，Δκ_t)，(Δω_t，Δκ_t)是一个随时间变化的角度值，则基于姿态改正角计算姿态改正矩阵R_C(t)为：

$R_C\left(t\right)=\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₃由Δω_t、和Δκ_t计算，公式如下：

将姿态改正矩阵R_C(t)代入到严密成像模型中，则可以得到附加了姿态改正量的严密成像模型：

${\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_{WGS-84}={\left[\begin{array}{c}{X}_{GPS}\\ Y_{GPS}\\ Z_{GPS}\end{array}\right]}_{WGS-84}+{\mathrm{mR}}_C\left(t\right)R_{J2000}^{WGS-84}{R}_{Star}^{J2000}{R}_{body}^{star}{R}_{camera}^{body}\left[\begin{array}{c}tan\left({\mathrm{\ψ}}_X\right)\\ tan\left({\mathrm{\ψ}}_Y\right)\\ -1\end{array}\right]---\left(16\right)$

式中，为WGS-84坐标系下地面点的坐标，其中X，Y，Z分别为地面点的横坐标、纵坐标和高程值，为星上GPS装置测得的卫星质心的位置，X_CPS、Y_GPS、Z_GPS分别表示由星上载荷GPS接收机测量的卫星质心的横坐标，纵坐标和高度值，m为比例因子，其由卫星飞行高度决定，R_C(t)为t时刻的姿态改正矩阵，为J2000惯性坐标系转换到WGS-84坐标系的旋转矩阵，为卫星本体坐标系转换到J2000惯性坐标系的旋转矩阵，为相机坐标系转换到卫星本体坐标系的旋转矩阵，为探元指向角，由卫星研制方提供。

(2.5)利用步骤2.4建立的严密成像模型，基于地面控制点和外方位观测值计算姿态改正量；

其中，外方位观测值是指卫星在轨飞行所测量的线元素和角元素，即GPS和星敏陀螺测量值。

卫星的姿态改正矩阵R_C(t)可利用卫星对同一地区采用不同角度拍摄的立体影像对，通过基于定向片的光束法平差来计算得到。其流程如图6所示。步骤2.5具体包括以下子步骤：

(2.5.1)设置定向片个数，并计算控制点范围内的定向片成像时刻和定向片成像时刻之间的间隔；

针对某一控制地区，将已知地面坐标(一般表示为WGS84大地经纬度坐标表示(纬度B，经度L，高程H)的控制点，在立体影像(即两张或者三张同一地区不同视角观测的影像)上找到该控制点并标记像点坐标，已知地面经纬度坐标和影像像点坐标的点简称为像控点。判断影像中像控点在沿卫星飞行方向的最大和最小边界影像行，从卫星影像辅助元数据中获得最大和最小边界行对应行成像时刻为T_max，和T_min，设定参与计算的定向片(影像之间的连接点)的个数n，同时计算每个定向片的间隔时间T_span＝(T_max-T_min)/(n-1)，则每i个定向片的行时刻为T_i＝T₀+T_span*(i-1)。

该示例中，选取覆盖有地面控制的六个资源三号卫星影像数据，每个影像数据范围至少包括一景或以上，该遥感卫星为三线阵立体卫星。以其中一个控制点区域影像作为试验对象展开描述。该控制区域含28个像控点，分别确定像控点坐标在前，正，后三个线阵影像中沿轨方向范围为(76841，129342)，(129558，218503)，(77490，129955)(单位：像素)。根据行时文件，获取控制区域在三个影像上对应的行时范围(97931508.321636394，97931534.082430676)，(97931479.770396858，97931505.481362402)，(97931451.377631396，97931477.038832352)。以三视最大时间间隔(97931451.377631396，97931534.082430676)，设定4个定向片，则时间间隔：(97931534.082430676-97931451.377631396)/(4-1)＝27.5682秒。则第一个定向片行时为97931451.3776313秒，第二个定向片行时为97931451.3776313960+27.5682＝97931478.9458978180秒，第三个定向片行时为：97931451.3776313960+27.5682*2＝97931506.514164239秒，第四个定向片行时为：97931534.0824306610秒。

(2.5.2)把控制点的经纬度坐标(B，L，H)转换为地心直角坐标(X，Y，Z)；

由于公式(16)中为WGS84地心直角坐标系，为了满足严格模型计算的需要，在利用公式(16)前，把控制点大地经纬度坐标(B，L，H)转换为地心直角坐标(X，Y，Z)。

地心直角坐标系原点位于地心；Z轴为极轴，向北为正；X轴穿过本初子午线与赤道的交点；Y轴穿过赤道与东经90°的交点。

X＝(v+H)cosBcosL

Y＝(v+H)cosBsinL

Z＝((1-e²)v+H)sinB (18)

式中：v为WGS84椭球面纬度B处的卯酉曲率半径，v＝a/(1-e²sin²B)^0.5。B和L分别为坐标点的纬度和经度，H为相对WGS84椭球面的高度。E为椭球第一偏心率，a为WGS84椭球长半轴，a＝6378137.0米，b为WGS84椭球短半轴，b＝6356752.3142米。

该示例中，将第一个WGS84坐标系下大地坐标(114.3659，38.1715，82.133)(其中前两项单位为度，后一项单位为米)转化为地心坐标系为(-2071371.47923，4573549.931328，3920478.894641)(单位：米)。

(2.5.3)获得定向片时刻处姿态值，并设定定向片姿态改正数初值；

利用下传的卫星姿态数据，该姿态数据中包含采样时刻和对应的J2000坐标系下姿态四元数(q值表示)，利用第i个定向片处的行时刻Ti，采用球面内插如公式(19)内插出对应的姿态值由于本发明中仅计算定向片时刻的姿态改变量(即改正数)，因此，直接利用每个定向片内插出的姿态四元数值作为理论不变值。由于姿态改正量一般很小，迭代解算前的初值设为0即可。

$SLERP(t,{\stackrel{\·}{q}}_1,{\stackrel{\·}{q}}_2)=\stackrel{\·}{q}\left(t\right)={\stackrel{\·}{q}}_1{\left({\stackrel{\·}{q}}_1^{-1}{\stackrel{\·}{q}}_2\right)}^t---\left(19\right)$

其中，t为内插时刻，为内插结果，为内插时刻相邻四元数值。

根据步骤(2.5.1)得到的四个定向片的行时，其中第一个定向片行时为：97931451.3776313秒，第二个定向片行时为：97931478.9458978秒，第三个定向片行时为：97931506.5141642秒，第四个定向片行时为：97931534.0824306秒。

分别查找下传的姿态文件中距离待内插时间(四个定向片的时刻)最近的两组四元数，利用这两组最近的四元数进行四元数球面内插，可以获得四个定向片时刻的姿态值：

第一个定向片：q₄＝0.16294326

第二个定向片：q₄＝0.17031543

第三个定向片：q₄＝0.17775321

第四个定向片：q₄＝0.18494362。

(2.5.4)计算各控制点在立体像对(即左中右三张影像)的姿态和轨道；

同一控制点在立体像对上的两个像点对称为同名点，分别在立体影像上的像点坐标为(x1，y1)和(x2，y2)，x为垂轨行方向，y为沿轨列方向。卫星研制单位提供给用户沿轨方向每行的成像时间，利用立体像点沿轨方向坐标y1，y2，分别内插出各自对应的时刻t1，t2。

利用卫星下传的轨道位置文件，利用拉格朗日公式(如公式20)分别内插出当前立体像点对应的卫星位置和

利用设定的离散间隔T_span的采样时刻的定向片的姿态，利用多项式内插出立体影像上连接点的姿态如公式(21)所示。

区间[a，b]上任意一点x的n阶拉格朗日多项式的代数表达式为：

$f\left(x\right)=\underset{i=0}{\overset{n}{\Σ}}{y}_i\underset{j\≠i}{\overset{n}{\underset{j=0}{\Π}}}\left(\frac{x-x_i}{x_i-x_j}\right)---\left(20\right)$

其中，y＝f(x)是区间[a，b]上的一个实函数，x_i(i＝0，1，…，n)是[a，b]上n+1个互异实数，且y＝f(x)在x_i的值为y_i＝f(x_i)，则点x_i(i＝0，1，…，n)称为插值节点，包含插值节点的区间[a，b]称为差值区间。

n阶多项式拟合如下式所示，

$f\left(t\right)=\underset{i=0}{\overset{n}{\Σ}}{d}_i{t}^i-x_0+d_1{t}^1+...+d_n{t}^n---\left(21\right)$

式中，t为时间变量，n为多项式阶数，d₀～d_n为多项式系数。一般多项式阶数根据实际情况确定。

该示例中，步骤(2.5.2)中第一个WGS84坐标系下地面大地坐标(114.3659，38.1715，82.133)对应在资源三号卫星前正后三视上的像点坐标分别如下：

表2 像点坐标 单位：像素

在卫星辅助行时文件，可以查到后视相机该控制点像点在沿轨y方向(行方向)前后的累积秒，示例具体如下：

表3 后视像点累积秒 单位：秒

基于以上行方向的累积秒参考辅助文件，可以内插出该点的时刻t：97931248.1588635423。

利用轨道文件和公式(20)，计算出t时刻的轨道

[PX PY PZ]＝[-2028198.3 3501873.4 5549723.3]

利用定向片已知的姿态和公式(21)，计算出t时刻的姿态

q₄＝0.291054376

在卫星辅助行时文件，可以查到下视相机该控制点像点在沿轨y方向(行向)前后的累积秒，示例具体如下：

表4 下视像点累积秒 单位：秒

行方向	累积秒
		212162	97931248.0897352695
212163	97931248.0900262743

基于以上行方向的累积秒参考辅助文件，可以内插出该点的时刻t：97931248.089897534。

利用轨道文件和公式(20)，计算出t时刻的轨道

[PX PY PZ]＝[-2028912.3 3501297.5 5550265.8]

利用定向片已知的姿态和公式(21)，计算出t时刻的姿态

q₄＝0.291032164

在卫星辅助行时文件，可以查到前视相机该控制点像点在沿轨y方向(行方向)前后的累积秒，示例具体如下：

表5 前视像点累积秒 单位：秒

行方向	累积秒
		126199	97931248.0853606015
126200	97931248.0858536065

基于以上行方向的累积秒参考辅助文件，可以内插出该点的时刻t：97931248.0855964523。

利用轨道文件和公式(20)，计算出t时刻的轨道

[PX PY PZ]＝[-2028903.8 3501289.5 5550261.3]

利用定向片已知的姿态和公式(21)，计算出t时刻的姿态

$\stackrel{\‾}{q}=\[\begin{array}{ccc}{q}_1& q_2& q_3\end{array}\]={\left[\begin{array}{ccc}-0.77635714& 0.54870657& 0.29101164\end{array}\right]}^T,$ q₄＝0.2910012164

(2.5.5)获取控制点所在位置的CCD像元对应的指向角；

如步骤2.5.4所述，控制点垂直轨道方向(即沿卫星CCD方向)像点坐标为x1，x2，根据提供给用户的每片CCD的探测指向角值，进行相邻两个CCD的指向值进行线性内插。

基于相邻的CCD指向角线性内插如下所示：

$y=y_i\frac{x-x_{i+1}}{x_i-x_{i+1}}+y_{i+1}\frac{x-x_i}{x_{i+1}-x_i}---\left(22\right)$

其中，x_i(i＝0，…，n)为线阵CCD第i个像元的排列位置，x为待内插点位置，y为内插结果，x_i、x_i+1分别为i、i+1位置值，y_i、y_i+1分别为x_i、x_i+1处的CCD的指向角。n为线阵探元个数，y_i为x_i处像元获取影像的指向角信息。

该示例中，基于后视行列方向的像方坐标(3200.434 125615.633)，结合辅助文件可得相邻像方坐标位置处的整像素指向角值：

表6 后视整像素指向角 单位：秒

CCD行号	tan(ψX)	tan(ψY)
			3200	-0.031270036689982	-0.383183680349274
3201	-0.031263698193461	-0.383183677797243

根据以上行号对应的指向角可以内插出当前CCD的指向角分别为：tg fy＝-0.031267285782491886，tg fx＝-0.383183678733838377。

基于正视行列方向的像方坐标(5045.813 212162.438)，结合辅助文件可得相邻像方坐标位置处的整像素指向角值：

表7 正视整像素指向角 单位：秒

CCD行号	tan(ψX)	tan(ψY)
			5045	-0.029821872593494	-0.000802497768906
5046	-0.029817756788962	-0.000802499993538

根据以上行号对应的指向角可以内插出当前CCD的指向角分别为：tg fy＝-0.029818526444409484，tg fx＝-0.000802498743294816。

基于前视行列方向的像方坐标(3621.661 126199.547)，结合辅助文件可得相邻像方坐标位置处的整像素指向角值：

表8 前视整像素指向角 单位：秒

CCD行号	tan(ψX)	tan(ψY)
			3621	-0.029007144301723	0.383723636718796
3622	-0.029000804318194	0.383723639574706

根据以上行号对应的指向角可以内插出当前CCD的指向角分别为：tan(ψ_X)＝-0.029002953572610331，tan(ψ_Y)＝0.38372363828097877。

(2.5.6)将由卫星方提供的传感器在卫星本体上的安装矩阵代入到公式(16)中，并在计算中作为已知且不变，设定定向片姿态改正量初值为0；

(2.5.7)利用姿态解算模型，如公式(16)，对定向片姿态改正参数进行线性化，得到像控点对应的姿态改正数的误差方程；由于像控点的像点和对应的地面坐标已知，未知参数仅包括姿态改正数。

V＝AX-L (23)

其中，

L＝f-f₀， (25)

其中，f0为将未知数迭代前一次初值代入公式(24)计算得到的。

$f={\left[\begin{array}{c}{X}_{GPS}-X\\ Y_{GPS}-X\\ Z_{GPS}-X\end{array}\right]}_{WGS-84}+{\mathrm{mR}}_C\left(t\right)R_{J2000}^{WGS-84}{R}_{Star}^{J2000}{R}_{body}^{star}{R}_{camera}^{body}\left[\begin{array}{c}tan\left({\mathrm{\ψ}}_X\right)\\ \tan \left({\mathrm{\ψ}}_Y\right)\\ -1\end{array}\right],---\left(26\right)$

式中，为WGS-84坐标系下的与空间坐标的差值，为星上GPS装置测得的卫星质心的位置，为测得的地面控制点在WGS84坐标系下的坐标，m为比例因子，R_C(t)为t时刻的姿态改正矩阵，为J2000惯性坐标系转换到WGS-84坐标系的旋转矩阵，为卫星本体坐标系转换到J2000惯性坐标系的旋转矩阵，和共称为卫星测量的外方位元素，为相机坐标系转换到卫星本体坐标系的旋转矩阵，为探元指向角；

(2.5.8)利用最小二乘解算未知数，根据下视和前后视分辨率不同以及对平差结果的影响，给定不同的权，对误差方程进行法化，解算定向片的姿态改正量。对于遥感卫星，经过事后精密定轨，卫星轨道位置精度在5cm左右，定位误差主要来自于姿态角的误差，在进行外方位元素答解时，可视GPS坐标观测值为已知，直接计算姿态角误差改正数。对于n个采样点，每个采样点有3个待求未知参数，即姿态改正量(Δω_t，Δκ_t)，共3×n个未知参数。地面控制至少需要(3×n)/2个地面控制点，即可解求这n个采样点对应的(Δω_t，Δκ_t)

如果平差的未知数改正数小于阈值1e-10度，则认为迭代收敛，否则更新未知数的值继续迭代。整个迭代求解的过程如图9所示，即首先根据姿态误差初步规律确定定向片的间隔，把控制点坐标转换为地心坐标，设定各定向片改正数的初值为0，利用前方交会计算各连接点的地心坐标，列出各控制点、连接点在各影像的误差方程，利用最小二乘解算定向片时刻的姿态改正数和各加密点的坐标，代入迭代，当满足阈值则输出定向片的姿态改正值，否则更新改正值和加密点坐标值重新列出各控制点、连接点在各影像的误差方程进行循环计算知道满足要求。

输出平差结果。把定向片的姿态改正量和成像时间作为平差结果如下表所示。

表9 某地区4个定向片计算的姿态误差

同样方法，选取另外两个控制区域及覆盖影像数据，同样取4个定向片，基于上述步骤计算对应的三轴姿态误差如下表所示。

表10 某地区4个定向片计算的姿态误差

表11 某地区4个定向片计算的姿态误差

(3)利用步骤1计算的初值和步骤2计算的姿态改正量拟合出长周期误差，获得长周期误差修正模型，其中，步骤(3)具体包括以下步骤：

(3.1)针对步骤2解求的姿态改正量(Δω_t，Δκ_t)，构建卫星姿态长周期误差模型：

为了描述方便，这里将模型进行简化。给定阶数m＝1，且略去b_{Δω_i}sin(ixω_Δω)，b_{Δκ_i}sin(ixω_Δκ)部分，以上公式简化为：

其中，振幅a_{Δω_1}、a_{Δκ_1}，截距c_Δω、c_Δκ，频率ω_Δω、ω_Δκ为待求参数。

(3.2)基于最小二乘原理对振幅a_{Δω_1}、a_{Δκ_1}，截距c_Δω、c_Δκ，频率ω_Δω、ω_Δκ参数进行解算，将解算得到的最终的参数值a_{Δω_1}，a_{Δκ_1}，ω_Δω，ω_Δκ，c_Δω，c_Δκ代入到所述卫星姿态长周期误差模型，即得到了长周期误差修正模型。

其中，步骤(3.2)具体为：

根据步骤(3.1)中的卫星姿态长周期误差模型，列出误差方程形式：

其中，

f_{ω_i}，f_{κ_i}表示ω，k分别对应的第i个测量值，Δa_{Δω_1}，Δa_{Δθ_1}，Δω_{Δω_1}，Δω_{Δκ_1}，Δc_{Δω_1}，Δc_{Δκ_1}表示待求的未知参数，为振幅、频率和截距的初值，其中，振幅、频率的初值为在步骤1中计算获得的振幅和频率的初值，截距c_Δω、c_Δκ的初值设为0。

由步骤(1.5)获得的周期初值为5671.81秒，x为时间Δt＝t-t₀，这里t₀＝95777428.69秒，则角频率ω_Δω、和ω_Δκ初值设为2*phi/5671.81，，振幅a_{Δω_1}、和a_{Δκ_1}均设为

根据最小二乘原理，最小二乘的解为：

将各参数初值代入上式，根据最小二乘计算求得未知参数改正量Δa_{Δω_1}，Δa_{Δκ_1}，Δω_Δω，Δω_Δκ，Δc_Δω，Δc_Δκ，并计算出对应的残差V_Δω，V_Δκ，若结果残差中误差小于给定阈值K＝0.00001度，则计算结束；否则累加未知参数值Δa_{Δω_1}，Δa_{Δκ_1}，Δω_Δω，Δω_Δκ，Δc_Δω，Δc_Δκ得到新的参数a_{Δω_1}，a_{Δκ_1}，ω_Δω，ω_Δκ，c_Δω，c_Δκ，重新代入上式中解算。计算流程如图6所示。

迭代结束后，得到最终的参数值a_{Δω_1}，a_{Δκ_1}，ω_Δω，ω_Δκ，c_Δω，c_Δκ，代入到卫星姿态长周期误差模型(下式)，即得到了本发明最终待求的长周期误差修正模型。

整个步骤3可以利用图10所示来描述具体流程。即，首先给定阶数，然后给定频率，进行平差，在平差结果满足精度要求的情况下结束，如果不满足改变阶数、频率，重新进行平差，直到满足结果，通过步骤3构建了精确的卫星姿态长周期误差修正模型，利用该模型即可完成姿态长周期误差的抑制。

由于姿态长周期误差的周期较差，传统的方法无法探测，更无法消除其对卫星测绘精度的影响，如表9至表11表示，若使用未经本发明处理的姿态，长周期绝对姿态误差达到3角秒(均值)，对于500km轨道高度，其定位误差则约为7.5m；采用本方法构建模型后，如表12所示，拟合残差达到0.25角秒，地面定位精度达到约1m。采用本方法进行卫星姿态长周期误差的消除有助于卫星测绘几何定位精度的提高，为卫星测绘向大比例尺测图的发展提供有力支撑。

表12 解算的模型参数

将表12的参数，代入到长周期误差模型中，则求得的长周期误差模型为：

所述步骤(3)选取符合误差规律变化趋势的高精度拟合模型，采用最小二乘估计方法对参数进行最优估计。在解决模型选取、高精度采样数据以及参数估计方法的基础上，实现星敏感器长周期误差项拟合与修正，从而提高与陀螺联合定姿的绝对精度。

以上内容仅为本发明的较佳实施例，对于本领域的普通技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，本说明书内容不应理解为对本发明的限制。

Claims (8)

1.一种基于地面控制的卫星长周期姿态误差修正方法，其特征在于，所述方法包括以下步骤：

步骤1，基于星敏感器光轴夹角计算长周期姿态误差的频率和振幅初值，其中，步骤1具体包括以下子步骤：

步骤1.1，随机抽取一段卫星下传的数据，得到一组离散的以归一化四元数形式表示的姿态数据，所述姿态数据为第一星敏感器和第二星敏感器测量的姿态数据，其中，所述第一星敏感器测量的姿态数据优选的用四元数表示为：

$q=\left[\begin{array}{c}\stackrel{\‾}{q}\\ q_4\end{array}\right]=\left[\begin{array}{c}{q}_1\\ q_2\\ q_3\\ q_4\end{array}\right],$

其中，为矢量，q₄为标量，且满足约束方程：

q₁ ²+q₂ ²+q₃ ²+q₄ ²＝1；

所述第二星敏感器测量的姿态数据用四元数表示为：

$q=\left[\begin{array}{c}{\stackrel{\‾}{q}}^{\′}\\ q_4^{\′}\end{array}\right]=\left[\begin{array}{c}{q}_1^{\′}\\ q_2^{\′}\\ q_3^{\′}\\ q_4^{\′}\end{array}\right],$

其中，为矢量，q′₄为标量，且满足约束方程：

q′₁ ²+q′₂ ²+q′₃ ²+q′₄ ²＝1。

步骤1.2，计算所述第一星敏感器和第二星敏感器的光轴矢量，其中：

所述第一星敏感器在某一个时刻的光轴矢量A为：

A＝(2(q₁q₃+q₂q₄) 2(q₂q₃-q₁q₄) -q₁ ²-q₂ ²+q₃ ²+q₄ ²)；

所述第二星敏感器在同一时刻的光轴矢量B为：

B＝(2(q′₁q′₃+q′₂q′₄) 2(q′₂q′₃-q′₁q′₄) -q′₁ ²-q′₂ ²+q′₃ ²+q′₄ ²)；

步骤1.3，基于所述第一星敏感器和第二星敏感器的光轴矢量，计算每一个时刻的所述第一星敏感器和第二星敏感器的光轴夹角实测值；

步骤1.4，基于所述第一星敏感器和第二星敏感器在卫星本体上的安装矢量，计算所述第一星敏感器和第二星敏感器的光轴夹角安装值；

步骤1.5，计算星敏感器的光轴夹角变化差值；

步骤1.6，计算星敏感器的光轴夹角变化差值的频率和振幅的初值，并将所述星敏感器的光轴夹角变化差值的频率和振幅的初值作为所述长周期误差的频率和振幅初值。

步骤2，利用地面控制点和外方位观测值计算姿态改正量；其中，步骤2具体包括以下子步骤：

步骤2.1，将相机坐标系转换到卫星本体坐标系；

步骤2.2，将卫星本体坐标系转换到J2000惯性坐标系；

步骤2.3，将J2000惯性坐标系转换到WGS-84坐标系；

步骤2.4，基于步骤2.1至步骤2.3的坐标系转换，构建附加了姿态改正量的严密成像模型；

步骤2.5，利用步骤2.4建立的严密成像模型，基于地面控制点和外方位观测值计算姿态改正量。

步骤3，利用步骤1计算的初值和步骤2计算的姿态改正量拟合出长周期误差，获得长周期误差修正模型；其中，步骤3具体包括以下步骤：

步骤3.1，针对步骤2解求的姿态改正量构建卫星姿态长周期误差模型：

其中，振幅a_{Δω_1}、a_{Δκ_1}，截距c_Δω、c_Δκ，频率ω_Δω、ω_Δκ为待求参数。

步骤3.2，基于最小二乘原理对振幅a_{Δω_1}、a_{Δκ_1}，截距c_Δω、c_Δκ，频率ω_Δω、ω_Δκ进行解算，将解算得到的最终的参数值a_{Δω_1}，a_{Δκ_1}，ω_Δω，ω_Δκ，c_Δω，c_Δκ代入到所述卫星姿态长周期误差模型，即得到了长周期误差修正模型。

2.根据权利要求1所述的方法，其特征在于，所述步骤1.3具体包括：

根据公式(2.1)计算光轴夹角实测值θ：

$\mathrm{\θ}=arctan\frac{A\×B}{A\·B}---\left(2.1\right)$

其中，A和B分别为所述第一星敏感器和第二星敏感器的光轴矢量；

根据公式(2.1)，依次计算每一个时刻的第一星敏感器和第二星敏感器的光轴夹角θ_t1，θ_t2，θ_t3，…，θ_tm，…，θ_tN，其中，tm表示第m个采样时刻，θ_tm表示计算的tm时刻的第一星敏感器和第二星敏感器的光轴夹角，tN表示第N个采样时刻，N表示光轴夹角的总数。

3.根据权利要求2所述的方法，其特征在于，所述步骤1.4具体包括：

根据公式(3.1)计算所述第一星敏感器和第二星敏感器的光轴夹角安装值所述光轴夹角安装值为光轴夹角理论值；

$\stackrel{\‾}{\mathrm{\θ}}=arctan\frac{ST1\×ST2}{ST1\·ST2}---\left(1.3\right)$

其中，ST1和ST2分别为第一星敏感器和第二星敏感器在卫星本体上的安装矢量。

4.根据权利要求3所述的方法，其特征在于，所述步骤1.5具体包括：

依次计算每一时刻的星敏感器光轴夹角实测值θ_t1，θ_t2，θ_t3，…，θ_tm，…，θ_tN与光轴安装值的差Δθ_t1，Δθ_t2，Δθ_t3，…，Δθ_tm，…，Δθ_tN：

${\mathrm{\Δ\θ}}_{tm}={\mathrm{\θ}}_{tm}-\stackrel{\‾}{\mathrm{\θ}}$

其中，tm表示第m个采样时刻，Δθ_tm表示计算的tm时刻的第一星敏感器和第二星敏感器的光轴夹角差值，N表示计算的光轴夹角差值的总数。

5.根据权利要求1所述的方法，其特征在于，所述步骤1.6具体包括：

将光轴夹角变化差值的主波形设为y＝f(t)＝L cos(w×(t-t₀))，其中L为振幅，w为频率，t为时间，t₀为初始时刻偏移；

计算星敏感器光轴夹角变化差值的主波形的振幅L的初值为：

其中，ξ＝max(Δθ_t)，为光轴夹角变化差值的最大值，为光轴夹角变化差值的最小值；

查找两个光轴夹角变化差值的最大值ξ₁和ξ₂的对应的时间间隔T₁和T₂，则主波的周期的初值为ΔT＝T₁-T₂，频率的初值为w＝2π/ΔT。

6.根据权利要求1所述的方法，其特征在于，步骤2.4具体包括：

设待求量三轴姿态改正量为则基于姿态改正角计算姿态改正矩阵R_C(t)为：

$R_C\left(t\right)=\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₃由Δω_t、和Δκ_t计算，公式如下：

将姿态改正矩阵R_C(t)代入到严密成像模型中，则可以得到附加了姿态改正量的严密成像模型：

${\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_{WGS-84}={\left[\begin{array}{c}{X}_{GPS}\\ Y_{GPS}\\ Z_{GPS}\end{array}\right]}_{WGS-84}+{\mathrm{mR}}_C\left(t\right)R_{J2000}^{WGS-84}{R}_{Star}^{J2000}{R}_{body}^{star}{R}_{camera}^{body}\left[\begin{array}{c}\tan \left({\mathrm{\ψ}}_X\right)\\ \tan \left({\mathrm{\ψ}}_Y\right)\\ -1\end{array}\right]$

式中，为WGS-84坐标系下地面点的坐标，其中X，Y，Z分别为地面点的横坐标、纵坐标和高程值，为星上GPS装置测得的卫星质心的位置，X_GPS、Y_GPS、Z_GPS分别表示由星上载荷GPS接收机测量的卫星质心的横坐标，纵坐标和高度值，m为比例因子，其由卫星飞行高度决定，R_C(t)为t时刻的姿态改正矩阵，为J2000惯性坐标系转换到WGS-84坐标系的旋转矩阵，为卫星本体坐标系转换到J2000惯性坐标系的旋转矩阵，为相机坐标系转换到卫星本体坐标系的旋转矩阵，为探元指向角。

7.根据权利要求1所述的方法，其特征在于，步骤2.5具体包括以下步骤：

步骤2.5.1，设置定向片个数，并计算控制点范围内的定向片成像时刻和定向片成像时刻之间的间隔；

步骤2.5.2，将控制点的经纬度坐标(B，L，H)转换为地心直角坐标(X，Y，Z)；

步骤2.5.3，获得定向片时刻处姿态值，并设定定向片姿态改正数初值；

步骤2.5.4，计算各控制点在立体像对的姿态和轨道；

步骤2.5.5，获取控制点所在位置的CCD像元对应的指向角；

步骤2.5.6，将由卫星方提供的传感器在卫星本体上的安装矩阵代入到下述公式(16)中，并在计算中作为已知且不变，设定定向片姿态改正量初值为0：

${\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_{WGS-84}={\left[\begin{array}{c}{X}_{GPS}\\ Y_{GPS}\\ Z_{GPS}\end{array}\right]}_{WGS-84}+{\mathrm{mR}}_C\left(t\right)R_{J2000}^{WGS-84}{R}_{Star}^{J2000}{R}_{body}^{star}{R}_{camera}^{body}\left[\begin{array}{c}\tan \left({\mathrm{\ψ}}_X\right)\\ \tan \left({\mathrm{\ψ}}_Y\right)\\ -1\end{array}\right]$

步骤2.5.7，利用公式(16)，对定向片姿态改正参数进行线性化，得到像控点对应的姿态改正数的误差方程；

步骤2.5.8，利用最小二乘解算未知数，根据下视和前后视分辨率不同以及对平差结果的影响，给定不同的权，对误差方程进行法化，解算定向片的姿态改正量。

8.根据权利要求1所述的方法，其特征在于，所述步骤3.2具体包括：

步骤(3.2)具体为：

根据步骤(3.1)中的卫星姿态长周期误差模型，列出误差方程形式：

其中，

f_{ω_i}，f_{κ_i}表示ω，k分别对应的第i个测量值，Δa_{Δω_1}，Δa_{Δθ_1}，Δω_{Δω_1}，Δω_{Δκ_1}，Δc_{Δω_1}，Δc_{Δκ_1}表示待求的未知参数，为振幅、频率和截距的初值，其中，振幅、频率的初值为在步骤1中计算获得的振幅和频率的初值，截距c_Δω、c_Δκ的初值设为0。

根据最小二乘原理，最小二乘的解为：

将各参数初值代入上式，根据最小二乘计算求得未知参数改正量Δa_{Δω_1}，Δa_{Δκ_1}，Δω_Δω，Δω_Δκ，Δc_Δω，Δc_Δκ，并计算出对应的残差V_Δω，V_Δκ，若结果残差中误差小于给定阈值K＝0.00001度，则计算结束；否则累加未知参数值Δa_{Δω_1}，Δa_{Δκ_1}，Δω_Δω，Δω_Δκ，Δc_Δω，Δc_Δκ得到新的参数a_{Δω_1}，a_{Δκ_1}，ω_Δω，ω_Δκ，c_Δω，c_Δκ，重新代入上式中解算。迭代结束后，得到最终的参数值a_{Δω_1}，a_{Δκ_1}，ω_Δω，ω_Δκ，c_Δω，c_Δκ，代入到卫星姿态长周期误差模型，即得到了长周期误差修正模型。