CN105910607B - 基于地面控制的卫星长周期姿态误差修正方法 - Google Patents
基于地面控制的卫星长周期姿态误差修正方法 Download PDFInfo
- Publication number
- CN105910607B CN105910607B CN201610208310.1A CN201610208310A CN105910607B CN 105910607 B CN105910607 B CN 105910607B CN 201610208310 A CN201610208310 A CN 201610208310A CN 105910607 B CN105910607 B CN 105910607B
- Authority
- CN
- China
- Prior art keywords
- star sensor
- attitude
- satellite
- value
- optical axis
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000002715 modification method Methods 0.000 title abstract 2
- 230000003287 optical effect Effects 0.000 claims abstract description 106
- 238000012937 correction Methods 0.000 claims abstract description 87
- 238000000034 method Methods 0.000 claims abstract description 44
- 238000004364 calculation method Methods 0.000 claims abstract description 24
- 239000013598 vector Substances 0.000 claims description 45
- 230000008859 change Effects 0.000 claims description 41
- 239000011159 matrix material Substances 0.000 claims description 34
- 238000003384 imaging method Methods 0.000 claims description 26
- 238000009434 installation Methods 0.000 claims description 24
- 238000005070 sampling Methods 0.000 claims description 20
- 238000005259 measurement Methods 0.000 claims description 16
- 238000006243 chemical reaction Methods 0.000 claims description 5
- 239000000523 sample Substances 0.000 claims description 5
- 230000009466 transformation Effects 0.000 claims description 5
- 241000287196 Asthenes Species 0.000 claims 1
- 238000013507 mapping Methods 0.000 abstract description 9
- 238000001514 detection method Methods 0.000 abstract description 2
- 230000001737 promoting effect Effects 0.000 abstract 1
- 238000001444 catalytic combustion detection Methods 0.000 description 20
- 230000036544 posture Effects 0.000 description 15
- 238000001914 filtration Methods 0.000 description 10
- 230000001186 cumulative effect Effects 0.000 description 9
- 238000010586 diagram Methods 0.000 description 5
- 230000008569 process Effects 0.000 description 5
- 238000013461 design Methods 0.000 description 4
- 238000011161 development Methods 0.000 description 3
- 238000005457 optimization Methods 0.000 description 3
- 230000003068 static effect Effects 0.000 description 3
- 206010034719 Personality change Diseases 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 230000002441 reversible effect Effects 0.000 description 2
- 241000243251 Hydra Species 0.000 description 1
- 230000002567 autonomic effect Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000033228 biological regulation Effects 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000002401 inhibitory effect Effects 0.000 description 1
- 230000005764 inhibitory process Effects 0.000 description 1
- 238000007689 inspection Methods 0.000 description 1
- QRXWMOHMRWLFEY-UHFFFAOYSA-N isoniazide Chemical compound NNC(=O)C1=CC=NC=C1 QRXWMOHMRWLFEY-UHFFFAOYSA-N 0.000 description 1
- 238000011545 laboratory measurement Methods 0.000 description 1
- 238000012821 model calculation Methods 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 230000008520 organization Effects 0.000 description 1
- 230000000737 periodic effect Effects 0.000 description 1
- 238000012805 post-processing Methods 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 230000002829 reductive effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000012827 research and development Methods 0.000 description 1
- 230000001629 suppression Effects 0.000 description 1
- 230000001360 synchronised effect Effects 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
- G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
- G01C21/24—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 specially adapted for cosmonautical navigation
Landscapes
- Engineering & Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Astronomy & Astrophysics (AREA)
- Automation & Control Theory (AREA)
- General Physics & Mathematics (AREA)
- Navigation (AREA)
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，随机抽取一段卫星下传的数据，得到一组离散的以归一化四元数形式表示的姿态数据，所述姿态数据为第一星敏感器和第二星敏感器测量的姿态数据，其中，所述第一星敏感器测量的姿态数据优选的用四元数表示为：
其中，为矢量，q4为标量，且满足约束方程：
q1 2+q2 2+q3 2+q4 2＝1；
所述第二星敏感器测量的姿态数据用四元数表示为：
其中，为矢量，q′4为标量，且满足约束方程：
步骤1.2，计算所述第一星敏感器和第二星敏感器的光轴矢量，其中：
所述第一星敏感器在某一个时刻的光轴矢量A为：
A＝(2(q1q3+q2q4) 2(q2q3-q1q4) -q1 2-q2 2+q3 2+q4 2)；
所述第二星敏感器在同一时刻的光轴矢量B为：
步骤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)计算光轴夹角实测值θ：
其中，A和B分别为所述第一星敏感器和第二星敏感器的光轴矢量；
根据公式(2.1)，依次计算每一个时刻的第一星敏感器和第二星敏感器的光轴夹角θt1，θt2，θt3，…，θtm，…，θtN，其中，tm表示第m个采样时刻，θtm表示计算的tm时刻的第一星敏感器和第二星敏感器的光轴夹角，tN表示第N个采样时刻，N表示光轴夹角的总数。
优选的，所述步骤1.4具体包括：
根据公式(3.1)计算所述第一星敏感器和第二星敏感器的光轴夹角安装值所述光轴夹角安装值为光轴夹角理论值；
其中，ST1和ST2分别为第一星敏感器和第二星敏感器在卫星本体上的安装矢量。
优选的，所述步骤1.5具体包括：
依次计算每一时刻的星敏感器光轴夹角实测值θt1，θt2，θt3，…，θtm，…，θtN与光轴安装值的差Δθt1，Δθt2，Δθt3，…，Δθtm，…，ΔθtN：
其中，tm表示第m个采样时刻，Δθtm表示计算的tm时刻的第一星敏感器和第二星敏感器的光轴夹角差值，N表示计算的光轴夹角差值的总数。优选的，所述步骤1.6具体包括：
将光轴夹角变化差值的主波形设为y＝f(t)＝L cos(w×(t-t0))，其中L为振幅，w为频率，t为时间，t0为初始时刻偏移；
计算星敏感器光轴夹角变化差值的主波形的振幅L的初值为：
其中，ξ＝max(Δθi)，为光轴夹角变化差值的最大值，为光轴夹角变化差值的最小值；
查找两个光轴夹角变化差值的最大值ξ1和ξ2的对应的时间间隔T1和T2，则主波的周期的初值为ΔT＝T1-T2，频率的初值为w＝2π/ΔT。
优选的，步骤2.4具体包括：
设待求量三轴姿态改正量为(Δωt，Δκt)，则基于姿态改正角计算姿态改正矩阵RC(t)为：
其中，a1、a2、a3、b1、b2、b3、c1、c2和c3由Δωt、和Δκt计算，公式如下：
将姿态改正矩阵RC(t)代入到严密成像模型中，则可以得到附加了姿态改正量的严密成像模型：
式中，为WGS-84坐标系下地面点的坐标，其中X，Y，Z分别为地面点的横坐标、纵坐标和高程值，为星上GPS装置测得的卫星质心的位置，XGPS、YGPS、ZGPS分别表示由星上载荷GPS接收机测量的卫星质心的横坐标，纵坐标和高度值，m为比例因子，其由卫星飞行高度决定，RC(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：
步骤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)的原始测量数据，双星敏感器在卫星平台上固联安装，两个星敏感器的光轴呈一定夹角，相对关系理论上认为不变。
第一星敏感器测量的姿态数据用四元数表示为：
其中，为矢量，一组q1、q2和q3值描述了卫星的一种姿态，q4为标量，且满足约束方程：
q1 2+q2 2+q3 2+q4 2＝1 (2)
同理，第二星敏感器测量的姿态数据用四元数表示为：
其中，为矢量，一组q′1、q′2、和q′3、值描述了卫星的一种姿态，q′4为标量，且满足约束方程：
双星敏感器采样时间同步，以第一帧下传数据为例：(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个矢量。
从而，第一星敏感器的测量的姿态数据可以表示为：
q4＝0.66743721。
第二星敏感器的测量的姿态数据可以表示为：
q4′＝0.63957798。
步骤1.2，计算所述第一星敏感器和第二星敏感器的光轴矢量；
将观测的姿态数据以四元数表示，在直角坐标系下，惯性系与星敏感器本体坐标系下姿态转换关系如下式：
其中，为惯性坐标系下姿态，为星敏感器本体坐标系下姿态。R(q)为旋转矩阵，且R(q)为：
星敏感器本体坐标系的Zb轴即为星敏感器光轴，Zb轴在惯性坐标系中的坐标为：
(2(q1q3+q2q4) 2(q2q3-q1q4) -q1 2-q2 2+q3 2+q4 2)。
下面对第一星敏感器和第二星敏感器的光轴矢量分别说明：
第一星敏感器的旋转矩阵为：
第二星敏感器的旋转矩阵为：
第一星敏感器在当前时刻的光轴指向(即光轴的方向)用矢量A表示：
A＝(2(q1q3+q2q4) 2(q2q3-q1q4) -q1 2-q2 2+q3 2+q4 2)＝(-0.296532528，0.917861147，0.263816932)；
第二星敏感器在同一时刻的光轴指向用矢量B表示：
(1.3)基于所述第一星敏感器和第二星敏感器的光轴矢量，计算所述第一星敏感器和第二星敏感器的光轴夹角实测值；
两星敏感器的安装关系如图5所示。
根据两个矢量的点积公式如下：
A·B＝||A|| ||B||cosθ (4)
其中θ为光轴矢量A，B的夹角。
根据两个矢量的叉积公式如下：
A×B＝||A|| ||B||sinθ (5)
对叉积取模然后除以点积：
所以：
其中，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°)]；
则第一星敏感器和第二星敏感器的光轴安装值(即光轴夹角理论值)为：
其中，ST1和ST2分别为第一星敏感器和第二星敏感器在卫星本体上的安装矢量。
(1.5)计算星敏感器的光轴夹角变化差值；
具体的，依次计算每一时刻的星敏感器光轴夹角实测值θt1，θt2，θt3，…，θtm，…，θtN与光轴安装值(光轴夹角理论值)的差Δθt1，Δθt2，Δθt3，…，Δθtm，…，ΔθtN：
其中，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-t0))，其中L为振幅，w为频率，t为时间，t0为初始时刻偏移；
计算星敏感器光轴夹角变化差值的主波形的振幅L的初值为：
其中，ξ＝max(Δθi)，为光轴夹角变化差值的最大值，为光轴夹角变化差值的最小值；
查找两个光轴夹角变化差值的最大值ξ1和ξ2的对应的时间间隔T1和T2，则主波的周期的初值为ΔT＝T1-T2，频率的初值为w＝2π/ΔT。
如图12所示，可以看出波形变化，将光轴夹角变化差值的主波形设为y＝f(t)＝Lcos(w×(t-t0))，如图13所示，其中L为振幅，w为频率，t为时间，t0为初始时刻偏移。计算主波形的主要参数包括振幅和频率等，为步骤3提供初值。
首先计算星敏感器光轴夹角变化差值的主波形的振幅初值，选取略大于一个信号周期的数据，计算变化量的最大值ξ＝max(Δθi)＝10角秒，计算变化差值的最小值则根据波形特征，星敏感器光轴夹角变化差值的振幅的初值为：然后计算光轴夹角变化差值的频率初值，查找两个变化最大值ξ1和ξ2的对应的时间95783638.5和95777966.69(单位：秒)，则主波的周期为ΔT＝T1-T2＝5671.81秒，w＝2π/ΔT＝0.0011Hz。
由于星敏感器姿态可以表示成光轴和旋转角，所以光轴的夹角变化可以反映出三轴姿态也应呈现出相似的变化规律，比如二者的波形相似，频率相近，但振幅不等。步骤(1.6)计算的长周期误差的频率和振幅初值即为长周期误差的频率和振幅初值，利用光轴夹角差值变化规律的探测解决了姿态无静态参考量无法探测的问题，为步骤3的计算提供初值。
所述的步骤(1)主要采用下传全弧段原始星敏感器、陀螺等姿态测量数据。利用两个或多个星敏感器安装关系几乎不变的特性，采用两个星敏感器光轴夹角来探测星敏感器姿态测量中存在的长周期变化探测，解决了由于准确参考导致星上姿态系统误差变化无法探测的问题。
步骤2，利用地面控制点和外方位观测值计算姿态改正量。
所述步骤(2)针对长周期姿态误差随时间变化规律特性，从精确重构误差模型的角度，设计地面控制数据分布，确定启算景间距、数量以及启算景内地面控制点分布等，利用加入姿态变化项的严密成型模型，反算获取近乎等纬度间距离散姿态误差项。
从查阅到文献情况看，现有研究工作一般假定卫星在同一时刻可获得2个以上已知地面控制点信息，利用双矢量(多矢量)定姿反推得出卫星姿态，与卫星确定姿态对比获得系统误差信息。从原理看该方法只适用于相机为框幅式成像原理的卫星，对目前在航天遥感领域通常使用的线阵推扫式遥感相机并不适用，因为线阵式相机为行扫描成像，很难保证在卫星同一时刻获得2个以上地面点信息。针对这一问题，该步骤采用定向片平差方法，本发明从线阵相机的成像原理出发，以星敏感器系统误差为未知量，以卫星到已知地面点视线矢量为观测量直接建立观测方程，多个已知地面点可以建立多个观测方程，由最小二乘解算得到系统误差信息。该方法可以有效避开利用双矢量(多矢量)方法求解卫星姿态，因此并不要求卫星在同一时刻同时获得多个地面点信息。
步骤2主要包括以下子步骤：
(2.1)将相机坐标系转换到卫星本体坐标系；
当传感器线阵列以其卫星轨道方向为轴向倾斜Ψy角、传感器阵列在其卫星轨道面内以垂直轨道方向为轴倾斜Ψx角时
传感器在卫星本体上的安装矩阵为线阵相机CCD每一个探元在相机坐标系中有一个三维指向，如图2所示，它可以通过实验室给定的相机主距、主点，CCD尺寸大小和个数等进行计算，也可经过在轨几何检校重新标定，提供给每个CCD的指向角文件，可供用户几何定位使用。指向角的表达式如下：
在严密成像模型中进行归一化得：
线阵上每个CCD均给出对应的指向角，第i个CCD的指向角为
传感器坐标系的定义：传感器坐标系的原点在线阵投影中心，沿扫描线方向为X轴，沿飞行方向为Y轴，Z轴按照右手法则确定。
卫星本体坐标系的定义：本体坐标系的原点位于卫星的质心，X轴、Y轴和Z轴分别为卫星的三个主惯量轴。X轴沿着卫星横轴，Y轴沿着纵轴指向卫星飞行方向，Z轴按照右手法则确定。本体坐标系是卫星上仪器设备安装的基准参考坐标系。
(2.2)将卫星本体坐标系转换到J2000惯性坐标系；
在成像时刻t，星敏感器本体坐标系到J2000惯性坐标系的旋转矩阵为：
其中，q1，q2，q3，q4分别表示t时刻星敏感器在J2000惯性坐标系下的姿态四元数，具体如公式(1)(2)所描述。
本体到J2000惯性坐标系的转换：
其中，星敏感器在卫星本体系下的安装矩阵为卫星发射前实验室测定并提供给用户，为已知。
空间固定惯性J2000参考系定义：如图3所示，空间固定惯性参考系简称CIS(Conventional Inertial System，CIS)，常用来描述卫星的运动，一般卫星星历的计算都在该坐标系下完成。CIS的原点为地球质心，Z轴指向天球的北极，X轴指向春分点，Y轴按照右手法则确定。由于地球围绕太阳运动，春分点和北极点经常发生变化。因此，国际组织便规定以某个时刻的春分点和北极点为基准，建立协议空间固定惯性系统。CIS一般采用国际大地测量协会和国际天文学联合会会议于1984年启用的协议天球坐标系J2000。
(2.3)将J2000惯性坐标系转换到WGS-84坐标系。
其中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)是一个随时间变化的角度值，则基于姿态改正角计算姿态改正矩阵RC(t)为：
其中，a1、a2、a3、b1、b2、b3、c1、c2和c3由Δωt、和Δκt计算，公式如下：
将姿态改正矩阵RC(t)代入到严密成像模型中，则可以得到附加了姿态改正量的严密成像模型：
式中，为WGS-84坐标系下地面点的坐标，其中X，Y，Z分别为地面点的横坐标、纵坐标和高程值，为星上GPS装置测得的卫星质心的位置，XCPS、YGPS、ZGPS分别表示由星上载荷GPS接收机测量的卫星质心的横坐标，纵坐标和高度值，m为比例因子，其由卫星飞行高度决定，RC(t)为t时刻的姿态改正矩阵，为J2000惯性坐标系转换到WGS-84坐标系的旋转矩阵，为卫星本体坐标系转换到J2000惯性坐标系的旋转矩阵，为相机坐标系转换到卫星本体坐标系的旋转矩阵，为探元指向角，由卫星研制方提供。
(2.5)利用步骤2.4建立的严密成像模型，基于地面控制点和外方位观测值计算姿态改正量；
其中，外方位观测值是指卫星在轨飞行所测量的线元素和角元素，即GPS和星敏陀螺测量值。
卫星的姿态改正矩阵RC(t)可利用卫星对同一地区采用不同角度拍摄的立体影像对，通过基于定向片的光束法平差来计算得到。其流程如图6所示。步骤2.5具体包括以下子步骤：
(2.5.1)设置定向片个数，并计算控制点范围内的定向片成像时刻和定向片成像时刻之间的间隔；
针对某一控制地区，将已知地面坐标(一般表示为WGS84大地经纬度坐标表示(纬度B，经度L，高程H)的控制点，在立体影像(即两张或者三张同一地区不同视角观测的影像)上找到该控制点并标记像点坐标，已知地面经纬度坐标和影像像点坐标的点简称为像控点。判断影像中像控点在沿卫星飞行方向的最大和最小边界影像行，从卫星影像辅助元数据中获得最大和最小边界行对应行成像时刻为Tmax，和Tmin，设定参与计算的定向片(影像之间的连接点)的个数n，同时计算每个定向片的间隔时间Tspan＝(Tmax-Tmin)/(n-1)，则每i个定向片的行时刻为Ti＝T0+Tspan*(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-e2)v+H)sinB (18)
式中：v为WGS84椭球面纬度B处的卯酉曲率半径，v＝a/(1-e2sin2B)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即可。
其中，t为内插时刻，为内插结果，为内插时刻相邻四元数值。
根据步骤(2.5.1)得到的四个定向片的行时，其中第一个定向片行时为：97931451.3776313秒，第二个定向片行时为：97931478.9458978秒，第三个定向片行时为：97931506.5141642秒，第四个定向片行时为：97931534.0824306秒。
分别查找下传的姿态文件中距离待内插时间(四个定向片的时刻)最近的两组四元数，利用这两组最近的四元数进行四元数球面内插，可以获得四个定向片时刻的姿态值：
第一个定向片：q4＝0.16294326
第二个定向片：q4＝0.17031543
第三个定向片：q4＝0.17775321
第四个定向片：q4＝0.18494362。
(2.5.4)计算各控制点在立体像对(即左中右三张影像)的姿态和轨道；
同一控制点在立体像对上的两个像点对称为同名点，分别在立体影像上的像点坐标为(x1，y1)和(x2，y2)，x为垂轨行方向，y为沿轨列方向。卫星研制单位提供给用户沿轨方向每行的成像时间，利用立体像点沿轨方向坐标y1，y2，分别内插出各自对应的时刻t1，t2。
利用卫星下传的轨道位置文件，利用拉格朗日公式(如公式20)分别内插出当前立体像点对应的卫星位置和
利用设定的离散间隔Tspan的采样时刻的定向片的姿态，利用多项式内插出立体影像上连接点的姿态如公式(21)所示。
区间[a，b]上任意一点x的n阶拉格朗日多项式的代数表达式为：
其中，y＝f(x)是区间[a，b]上的一个实函数，xi(i＝0，1，…，n)是[a，b]上n+1个互异实数，且y＝f(x)在xi的值为yi＝f(xi)，则点xi(i＝0，1，…，n)称为插值节点，包含插值节点的区间[a，b]称为差值区间。
n阶多项式拟合如下式所示，
式中，t为时间变量，n为多项式阶数，d0～dn为多项式系数。一般多项式阶数根据实际情况确定。
该示例中，步骤(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时刻的姿态
q4＝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时刻的姿态
q4＝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时刻的姿态
q4＝0.2910012164
(2.5.5)获取控制点所在位置的CCD像元对应的指向角；
如步骤2.5.4所述，控制点垂直轨道方向(即沿卫星CCD方向)像点坐标为x1，x2，根据提供给用户的每片CCD的探测指向角值，进行相邻两个CCD的指向值进行线性内插。
基于相邻的CCD指向角线性内插如下所示：
其中，xi(i＝0，…，n)为线阵CCD第i个像元的排列位置，x为待内插点位置，y为内插结果，xi、xi+1分别为i、i+1位置值，yi、yi+1分别为xi、xi+1处的CCD的指向角。n为线阵探元个数，yi为xi处像元获取影像的指向角信息。
该示例中，基于后视行列方向的像方坐标(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-f0， (25)
其中，f0为将未知数迭代前一次初值代入公式(24)计算得到的。
式中，为WGS-84坐标系下的与空间坐标的差值，为星上GPS装置测得的卫星质心的位置，为测得的地面控制点在WGS84坐标系下的坐标，m为比例因子，RC(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Δω_isin(ixωΔω)，bΔκ_isin(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-t0，这里t0＝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，随机抽取一段卫星下传的数据，得到一组离散的以归一化四元数形式表示的姿态数据，所述姿态数据为第一星敏感器和第二星敏感器测量的姿态数据，其中，所述第一星敏感器测量的姿态数据用四元数表示为：
其中，为矢量，q4为标量，且满足约束方程：
q1 2+q2 2+q3 2+q4 2＝1；
所述第二星敏感器测量的姿态数据用四元数表示为：
其中，为矢量，q′4为标量，且满足约束方程：
步骤1.2，计算所述第一星敏感器和第二星敏感器的光轴矢量，其中：
所述第一星敏感器在某一个时刻的光轴矢量A为：
A＝(2(q1q3+q2q4) 2(q2q3-q1q4) -q1 2-q2 2+q3 2+q4 2)；
所述第二星敏感器在同一时刻的光轴矢量B为：
步骤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)计算光轴夹角实测值θ：
其中，A和B分别为所述第一星敏感器和第二星敏感器的光轴矢量；
根据公式(2.1)，依次计算每一个时刻的第一星敏感器和第二星敏感器的光轴夹角θt1,θt2,θt3,L,θtm,L,θtN，其中，tm表示第m个采样时刻，θtm表示计算的tm时刻的第一星敏感器和第二星敏感器的光轴夹角，tN表示第N个采样时刻，N表示光轴夹角的总数。
3.根据权利要求2所述的方法，其特征在于，所述步骤1.4具体包括：
根据公式(3.1)计算所述第一星敏感器和第二星敏感器的光轴夹角安装值所述光轴夹角安装值为光轴夹角理论值；
其中，ST1和ST2分别为第一星敏感器和第二星敏感器在卫星本体上的安装矢量。
4.根据权利要求3所述的方法，其特征在于，所述步骤1.5具体包括：
依次计算每一时刻的星敏感器光轴夹角实测值
θt1,θt2,θt3,L,θtm,L,θtN与光轴安装值的差
Δθt1,Δθt2,Δθt3,L,Δθtm,L,ΔθtN：
其中，tm表示第m个采样时刻，Δθtm表示计算的tm时刻的第一星敏感器和第二星敏感器的光轴夹角差值，N表示计算的光轴夹角差值的总数。
5.根据权利要求1所述的方法，其特征在于，所述步骤1.6具体包括：
将光轴夹角变化差值的主波形设为y＝f(t)＝Lcos(w×(t-t0))，其中L为振幅，w为频率，t为时间，t0为初始时刻偏移；
计算星敏感器光轴夹角变化差值的主波形的振幅L的初值为：
其中，ξ＝max(Δθi)，为光轴夹角变化差值的最大值，为光轴夹角变化差值的最小值；
查找两个光轴夹角变化差值的最大值ξ1和ξ2的对应的时间间隔T1和T2，则主波的周期的初值为ΔT＝T1-T2，频率的初值为w＝2π/ΔT。
6.根据权利要求1所述的方法，其特征在于，步骤2.4具体包括：
设待求量三轴姿态改正量为则基于姿态改正角计算姿态改正矩阵RC(t)为：
其中，a1、a2、a3、b1、b2、b3、c1、c2和c3由Δωt、和Δκt计算，
公式如下：
将姿态改正矩阵RC(t)代入到严密成像模型中，则可以得到附加了姿态改正量的严密成像模型：
式中，为WGS-84坐标系下地面点的坐标，其中X,Y,Z分别为地面点的横坐标、纵坐标和高程值，为星上GPS装置测得的卫星质心的位置，XGPS、YGPS、ZGPS分别表示由星上载荷GPS接收机测量的卫星质心的横坐标，纵坐标和高度值，m为比例因子，其由卫星飞行高度决定，RC(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，将由卫星方提供的传感器在卫星本体上的安装矩阵代入到下述公式中，并在计算中作为已知且不变，设定定向片姿态改正量初值为0：
式中，为WGS-84坐标系下地面点的坐标，其中X,Y,Z分别为地面点的横坐标、纵坐标和高程值，为星上GPS装置测得的卫星质心的位置，XGPS、YGPS、ZGPS分别表示由星上载荷GPS接收机测量的卫星质心的横坐标，纵坐标和高度值，m为比例因子，其由卫星飞行高度决定，RC(t)为t时刻的姿态改正矩阵，为J2000惯性坐标系转换到WGS-84坐标系的旋转矩阵，为卫星本体坐标系转换到J2000惯性坐标系的旋转矩阵，为相机坐标系转换到卫星本体坐标系的旋转矩阵，为探元指向角；
步骤2.5.7，利用步骤2.5.6的上述公式，对定向片姿态改正参数进行线性化，得到像控点对应的姿态改正数的误差方程；
步骤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Δκ，代入到卫星姿态长周期误差模型，即得到了长周期误差修正模型。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610208310.1A CN105910607B (zh) | 2016-04-07 | 2016-04-07 | 基于地面控制的卫星长周期姿态误差修正方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610208310.1A CN105910607B (zh) | 2016-04-07 | 2016-04-07 | 基于地面控制的卫星长周期姿态误差修正方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN105910607A CN105910607A (zh) | 2016-08-31 |
| CN105910607B true CN105910607B (zh) | 2018-11-13 |
Family
ID=56744799
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201610208310.1A Active CN105910607B (zh) | 2016-04-07 | 2016-04-07 | 基于地面控制的卫星长周期姿态误差修正方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN105910607B (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115147313A (zh) * | 2022-09-01 | 2022-10-04 | 中国科学院空天信息创新研究院 | 椭圆轨道遥感图像的几何校正方法、装置、设备及介质 |
Families Citing this family (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108571981B (zh) * | 2018-03-28 | 2021-08-13 | 西安爱生技术集团公司 | 一种侦察无人机目标定位集成检校方法 |
| CN108820253B (zh) * | 2018-04-16 | 2020-10-02 | 上海微小卫星工程中心 | 一种轨道短时失效情况下对地定向姿态的计算方法 |
| CN109033604B (zh) * | 2018-07-18 | 2022-10-11 | 哈尔滨工业大学 | 含旋转载荷的卫星动力学建模及轴承处受力的确定方法 |
| CN109696179B (zh) * | 2018-11-15 | 2022-10-18 | 上海航天控制技术研究所 | 一种遥感卫星星敏感器热弹性误差估计方法 |
| CN109581428B (zh) * | 2018-12-05 | 2022-02-18 | 上海航天计算机技术研究所 | 一种基于光学影像的在轨自修正的定位方法 |
| CN111552751B (zh) * | 2020-04-22 | 2021-04-27 | 中国人民解放军战略支援部队信息工程大学 | 三维地标控制点生成及应用方法、生成及应用装置 |
| CN114061382B (zh) * | 2021-09-17 | 2023-05-12 | 中国人民解放军63875部队 | 一种基于中远距离下中轴矢量交会姿态测量的精度预估仿真方法 |
| CN114396934B (zh) * | 2022-01-24 | 2022-12-09 | 自然资源部国土卫星遥感应用中心 | 顾及卫星周期误差的姿态优化方法 |
| CN114858133B (zh) * | 2022-04-21 | 2023-01-17 | 武汉大学 | 一种恒星观测模式下姿态低频误差修正方法 |
| CN114972078B (zh) * | 2022-05-09 | 2023-05-05 | 安徽大学 | 应用sar影像提升国产光学卫星影像无控几何质量方法及系统 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102506893A (zh) * | 2011-09-29 | 2012-06-20 | 北京控制工程研究所 | 一种基于地标信息的星敏感器低频误差补偿方法 |
| CN103323026A (zh) * | 2013-05-30 | 2013-09-25 | 北京控制工程研究所 | 星敏感器和有效载荷的姿态基准偏差估计与修正方法 |
| CN104729537A (zh) * | 2015-03-19 | 2015-06-24 | 北京控制工程研究所 | 一种星敏感器低频误差在轨实时补偿方法 |
-
2016
- 2016-04-07 CN CN201610208310.1A patent/CN105910607B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102506893A (zh) * | 2011-09-29 | 2012-06-20 | 北京控制工程研究所 | 一种基于地标信息的星敏感器低频误差补偿方法 |
| CN103323026A (zh) * | 2013-05-30 | 2013-09-25 | 北京控制工程研究所 | 星敏感器和有效载荷的姿态基准偏差估计与修正方法 |
| CN104729537A (zh) * | 2015-03-19 | 2015-06-24 | 北京控制工程研究所 | 一种星敏感器低频误差在轨实时补偿方法 |
Non-Patent Citations (3)
| Title |
|---|
| 卫星摄影姿态测定系统低频误差补偿;王任享等;《测绘学报》;20160229;第45卷(第2期);第127-130页 * |
| 基于严格成像模型的卫星姿态角系统误差补偿;雷玉飞等;《航天器工程》;20120229;第21卷(第1期);第25-30页 * |
| 基于地标信息的星敏感器低频误差标定方法;熊凯等;《空间控制技术与应用》;20120630;第38卷(第3期);第11-15页 * |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115147313A (zh) * | 2022-09-01 | 2022-10-04 | 中国科学院空天信息创新研究院 | 椭圆轨道遥感图像的几何校正方法、装置、设备及介质 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN105910607A (zh) | 2016-08-31 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN105910607B (zh) | 基于地面控制的卫星长周期姿态误差修正方法 | |
| CN102741706B (zh) | 地理参照图像区域的方法 | |
| CN104764443B (zh) | 一种光学遥感卫星严密成像几何模型构建方法 | |
| KR102005620B1 (ko) | 정지궤도 관측위성 영상 기하보정 장치 및 방법 | |
| Di et al. | A Self-Calibration Bundle Adjustment Method for Photogrammetric Processing of Chang $^{\prime} $ E-2 Stereo Lunar Imagery | |
| CN107144293A (zh) | 一种视频卫星面阵相机的几何定标方法 | |
| CN105698764A (zh) | 一种光学遥感卫星影像时变系统误差建模补偿方法及系统 | |
| Poli | A rigorous model for spaceborne linear array sensors | |
| Oberst et al. | The Phobos geodetic control point network and rotation model | |
| US9383210B2 (en) | Image navigation and registration (INR) transfer from exquisite systems to hosted space payloads | |
| CN111473802A (zh) | 一种基于线阵推扫的光学传感器内方位元素定标方法 | |
| CN111912430B (zh) | 高轨光学卫星的在轨几何定标方法、装置、设备及介质 | |
| CN106885585A (zh) | 一种基于光束法平差的星载摄影测量系统一体化检校方法 | |
| Tang et al. | Geometric accuracy analysis model of the Ziyuan-3 satellite without GCPs | |
| Poli | General model for airborne and spaceborne linear array sensors | |
| CN111156969A (zh) | 一种宽幅遥感影像立体测绘方法及系统 | |
| Pan et al. | Systematic geolocation errors of FengYun-3D MERSI-II | |
| Zhan et al. | Adaptive celestial positioning for the stationary mars rover based on a self-calibration model for the star sensor | |
| RU2723199C1 (ru) | Способ и система определения ориентации космического аппарата в пространстве с автономной коррекцией эффекта аберрации света | |
| CN115200573B (zh) | 空间目标的测量装备定位方法、系统和存储介质 | |
| Jiang et al. | Correction of distortions in YG-12 high-resolution panchromatic images | |
| Paluszek et al. | Optical navigation system | |
| Radhadevi et al. | An algorithm for geometric correction of full pass TMC imagery of Chandrayaan-1 | |
| Zhong et al. | Construction and verification of rigorous imaging model of spaceborne linear array vertical orbit ring scanning sensor | |
| Monay et al. | Diwata-2 targeting assessment and attitude error determination using a quaternion-based transformation system |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CP01 | Change in the name or title of a patent holder | ||
| CP01 | Change in the name or title of a patent holder |
Address after: 101300 No. 28 Lianhuachi West Road, Haidian District, Beijing Patentee after: Ministry of Natural Resources Land Satellite Remote Sensing Application Center Address before: 101300 No. 28 Lianhuachi West Road, Haidian District, Beijing Patentee before: Satellite Surveying and Mapping Application Center, NASG |