太阳敏感器与地球敏感器结合[2-3]，可以确定卫星三轴姿态。然而基于光电二极管的太阳敏感器，其输出的电压是太阳光和地球反照光等杂光辐照的叠加，如果忽略地球反照光的影响，直接将测量电压转换为太阳矢量，会导致转换结果存在较大误差，进而影响姿态估计的精度。为此需要考虑地球反照光，建立精确的光电二极管量测模型，才能实现高精度的卫星姿态估计。
文献[4]提出了地球反照光的数学模型及其计算方法，该方法需要实时的球面积分运算，计算较为复杂。文献[5-8]对其进行了简化。文献[5]假设地球表面的反照率仅与地球纬度有关，与地球经度无关，建立地球反照光的查找表并将其拟合成曲面，在已知卫星姿态和卫星轨道位置的情况下可以快速计算地球反照光影响。这种方法前期计算较为复杂且仅适用于单一轨道。文献[6]采用随机实验的方法得到地球反照光统计特性，然后根据统计的均值和方差补偿太阳矢量模型，这种方法对太阳矢量的动态范围有很大限制。文献[7-8]假设地球反照光为与卫星天底方向反向平行的单一矢量，简化地球反照光的计算，但仍存在计算复杂和精度较低的问题。
本文结合地球反照光模型，提出光电二极管动态偏置混合高斯噪声的量测模型，并应用滑窗估计和随机权重策略动态更新地球反照光偏置和噪声方差。在权重更新过程中，引入Huber影响函数处理量测残差的异常值，从而提高了算法鲁棒性和参数估计精度。同时，针对每个光电二极管的地球反照光分布不一，采用多比例因子分别估计每个光电二极管的地球反照光影响。基于本文建立的地球反照光校正的光电二极管太阳敏感器量测模型，结合地球敏感器和陀螺等传感器，采用无迹卡尔曼滤波(Unscented Kalman Filter, UKF)算法[9]实现了卫星的高精度姿态估计。
1 地球反照光校正的光电二极管量测模型 1.1 光电二极管工作原理 光电二极管的输出电压和太阳矢量与光电二极管敏感轴夹角的余弦成正比，根据电压输出，可以感测太阳方位。光电二极管常安装于卫星表面，易受周围地球反照光等杂光的干扰。图 1为地球反照光几何示意图，照射在积分区域(太阳照射与卫星可视方位区域交集)的太阳光，经地球表面漫反射后，共同作用于卫星方向。地球反照光的计算不仅与地球表面的反照率有关，还和太阳、卫星、地球相对位置有关，此外，地球反照光还与卫星本身的姿态有关。因此地球反照光的分布函数较为复杂，其计算较为繁琐。
 |
| 图 1 地球反照光几何示意图 Fig. 1 Geometric sketch of the earth's albedo |
| 图选项 |
地球反照光的影响较大，可达到直射太阳光总量的20%~30%，不能简单地忽略其影响。考虑地球反照光的影响，单个光电二极管的输出电压V可以表示为[10]
 | (1) |
 | (2) |
 | (3) |
式中：Vd为太阳照射分量；Va为地球反照光分量；vV为模型误差的零均值高斯白噪声；n =[cos ? cos θ cos ?sin θ sin ?]T为光电二极管敏感轴的单位矢量方向，可由安装高度角?和方位角θ表示；sb是太阳矢量在本体坐标系下的表示，可由在惯性坐标下太阳矢量经姿态矩阵变换得到；ψ为光电二极管半视场角；A为地球表面卫星可视区域和太阳照射的交集区域；α为地球表面dA的反照率；nA为地球表面面元dA的法向单位矢量；s⊕为地球到太阳的矢量方向；rAB为地球表面面元dA到卫星的矢量方向；B表示卫星当前轨道位置；S表示在地球阴影区的轨道。
1.2 地球反照光校正 式(3)中地球反照光分量Va计算需要球面积分，运算复杂，难以在实际中得到应用。本文将地球反照光设成一动态偏置项，补偿在光电二极管的量测模型中，同时将动态偏置估计的误差和光电二极管本身量测的误差统一为混合高斯模型，建立光电二极管的动态偏置混合高斯量测模型：
 | (4) |
式中：下标k表示采样时间序列；rk为地球反照光动态偏置；vk为非高斯噪声，可用式(5)表示:
 | (5) |
其中：pN(vk)为已知的均值为0、方差为R1, k的高斯噪声概率密度函数；qN(vk)为未知的均值为0、方差为R2, k的污染噪声概率密度函数；参数ε∈(0, 1)为污染系数，用来控制污染噪声的强弱，本文取ε=0.05；p(vk)为整体量测噪声vk的概率密度函数，其方差为Rk，满足Rk=(1－ε)R1, k+εR2, k。
单个光电二极管仅能测量一个太阳矢量分量，至少需要有3个有效的光电二极管量测值才能求取完整的太阳矢量。对于多个光电二极管的量测模型，可以表示为
 | (6) |
矢量表示形式为
 | (7) |
式中: H =[n 1T n 2T … nMT]T，ni(i=1, 2, …, M)表示光电二极管i的安装方位, M为光电二极管的数量。
2 模型参数的在线估计和更新方法 式(7)中地球反照光校正的光电二极管量测模型的地球反照光偏置项 rk和噪声方差 Rk是未知的，需要在卫星姿态估计过程中在线估计并动态更新。本文采用滑窗估计和随机权重策略动态更新模型参数，如图 2所示，假设模型参数在窗口采样时间序列k－1, k－2, …, k－N范围内不变，对于历史数据计算的模型参数赋予不同的权重，获取当前时刻模型估计的参数。由式(7)，结合滑窗估计和随机权重算法，此时可以得到动态偏置和噪声方差的估计公式分别为
 | (8) |
 | (9) |
式中：
 |
| 图 2 滑窗估计和随机权重算法示意图 Fig. 2 Schematic diagram of sliding window estimation and random weighting algorithm |
| 图选项 |
随着量测噪声的统计变化，量测残差向量 ek－j将存在偏置，它的幅值将会增加[11]，可以选用量测残差的模值作为权重，如式(10)和式(11)所示：
 | (10) |
 | (11) |
然而光电二极管噪声分布较为复杂，当出现残差较大的异常值时，直接采用残差向量作为模值会引起权重分配不合理，进而导致偏差和方差估计产生较大误差。基于Huber影响函数的鲁棒技术可以有效地处理非高斯噪声的情况，其更改了后验噪声方差矩阵，降低了对异常值的灵敏度。因此，本文引入Huber影响函数处理量测残差的异常值，其表达式如下：
 | (12) |
 | (13) |
式中：引入 Tk使 ηk满足关于概率密度对称和边缘概率密度条件[12]。ψ(·)是Huber函数，其表达式为
 | (14) |
其中:sgn(·)是符号函数；kε的取值和污染系数ε有关。
另外，考虑到每个光电二极管的地球反照光的分布情况不同，当所有的光电二极管都采用统一的权重系数会降低对地球反照光的跟踪特性，因此本文采用多重比例因子分别估计地球反照光对每个光电二极管的影响，此时可得
 | (15) |
 | (16) |
式中：diag(·)表示将一个向量转换为对角矩阵。
3 基于UKF的卫星姿态估计 3.1 姿态估计状态方程 卫星姿态运动学[13]可以用四元数表示为
 | (17) |
 | (18) |
式中：ω =[ωx ωy ωz]T为卫星三轴角速度；q (t)为四元数[14]，以绕某一固定旋转轴 e 旋转一个角度?描述姿态旋转变换，其定义为q =[qvTq4]T，其中qv=[q1 q2 q3]T= e ·sin ?、q4=cos ?。可以用四元数表示从惯性坐标系到本体坐标系的旋转变换矩阵A (q):
 | (19) |
陀螺常用的数学模型为[15]
 | (20) |
式中：ω 为真实的相对惯性的角速度；
 | (21) |
其中：δ(t－τ)为Dirac delta函数。
此时可以建立卫星姿态估计的连续状态方程为
 | (22) |
其离散状态方程可通过龙格库塔方法实现。然而，如果直接应用姿态动力学直接应用于UKF, 预测的四元数不能保证其模值依然保持为1。常用的解决方法是使用3个无约束的四元数误差向量表示4个元素的四元数[16]。定义误差四元数δ q =[δ qv T δq4]T，一般使用罗德里格斯(GRP)表示误差四元数为
 | (23) |
式中：a为0到1区间的参数；f为放大因子。a=0和f=1表示的是Gibbs向量，a=1和f=1表示的修正罗德里格斯向量。从δ p 到δ q 的逆变换为
 | (24) |
 | (25) |
3.2 姿态估计量测方程 本文采用光电二极管和地球敏感器两种姿态传感器，需要将其测量值融入到量测方程中。光电二极管的量测方程另一种表达形式为
 | (26) |
式中：sref为太阳矢量在惯性下的表示，可以查找星历表获得。
地球敏感器测量天底方向，其量测模型可以表示为
 | (27) |
式中：rearth为地球矢量在惯性系下的表示，可以由卫星轨道参数求解；εk的噪声方差为σST2 I3×3，I3×3为单位阵。
结合2个传感器测量模型可得量测方程：
 | (28) |
3.3 UKF算法实现 当已知状态更新的状态方程模型，且建立了状态和量测方程的噪声和误差统计模型，卡尔曼滤波方法采用递推的方式，从量测信息中实时提取出被估计量信息并存储在估计值中[17]。UKF是对线性卡尔曼滤波的改进，其不需要对状态方程和量测方程线性化，常用于姿态估计等非线性滤波算法中。UKF通过sigma点捕获系统真实的均值和方差，其精度可以达到泰勒展开式三阶近似。
定义离散系统的非线性状态方程和量测方程为
 | (29) |
式中：Xk=[δ qkT βk T]T为k时刻的状态量；Zk=[bk T Vk T]T为k时刻的量测量；qk =[rk T 01×3]T为系统动态偏置，rk为光电二极管模型的偏置参数，地球敏感器的偏置为0；Wk－1和 Vk分别为过程和量测方程的加性噪声, 其方差分别为 Qk和 Gk，表达式分别为
 | (30) |
 | (31) |
式中：Δt为采样时间间隔；Rk为光电二极管的噪声方差；σST2 I3×3为地球敏感器的噪声方差。UKF算法的具体实现可以参考文献[9]。
4 仿真校验 4.1 仿真条件 本文数值仿真选取卫星轨道高度400 km、倾角90°的圆轨道，轨道周期约为90 min, 选取卫星刚出地球阴影区到进入地球阴影区前的2 500 s时间段进行仿真。陀螺的输出采样频率为10 Hz, 常值漂移σv=1×10-4(°)/
 |
在出地球阴影时，结合光电二极管和地球敏感器测量的2个矢量，使用QUEST算法[18]给定卫星初始姿态。陀螺漂移初始值 β0=[0 0 0]T。
图 3为3个光电二极管的理想电压和量测电压的对比图，地球反照光影响为量测电压与理想电压的差值。由图中可以看出，对于同一个光电二极管，在轨道的不同时段, 地球反照光与太阳直射光的比值不同, 比值大约为0~25%。对于不同光电二极管，每个光电二极管的地球反照光分布情况不同。此外，在卫星刚出背光面和刚入背光面时，理想光电二极管和实际光电二极管的电压相差不大，地球反照光较弱，可以选择卫星刚出背光面时刻作为太阳矢量估计或者姿态估计的初始时刻。
 |
| 图 3 光电二极管1、2和3的理想电压与量测电压对比 Fig. 3 Comparison of ideal voltage and measured voltage for photodiode 1, 2 and 3 |
| 图选项 |
4.2 仿真结果 为了量化仿真结果，选取三轴欧拉角的模值作为评判指标：
 | (32) |
式中：θk、φk和ψk分别为横滚角、俯仰角和偏航角。为了保证结果的可靠性，参数使用蒙特卡罗仿真50次。基于本文建立的地球反照光校正的光电二极管太阳敏感器量测模型，结合地球敏感器和陀螺等传感器，采用UKF算法进行卫星姿态估计。
图 4对比了固定权重(如均值)、量测残差模值、量测残差经Huber影响函数处理后的模值、本文方法等4种权重选取策略的效果。从整体来看，在初始三轴卫星姿态1.5°左右时，采用本文建立的地球反照光校正模型和UKF算法，三轴姿态精度可以很快的收敛到0.5°甚至更高精度，验证了本文简化地球反照光模型具有一定的可行性。比较不同权重选择策略，可以看出采用本文方法精度较高，三轴姿态精度可以达到0.2°~0.3°。
 |
| 图 4 卫星三轴姿态估计误差对比 Fig. 4 Comparison of satellite three-axis attitude estimation error |
| 图选项 |
此外，还可以以地球反照光建模的动态偏置电压估计精度作为评价标准，地球反照光估计精度越高，光电二极管测量的太阳矢量精度就越高，进而姿态的估计精度越高。单个光电二极管的偏置估计误差可以通过多次蒙特卡罗方法求均值获得。表 1为光电二极管偏置估计误差均方根(RMS)结果，可以看出，本文权重选取策略可以有效提高偏置估计精度。图 5为3个光电二极管的偏置估计误差。从图 5(b)可以明显看出, 偏置估计的误差随时间推移而明显减小，偏置估计的精度越来越高。
表 1 光电二极管偏置估计误差RMS Table 1 RMS of photodiode bias estimation error
| mV | |||
| 权重选取策略 | 光电二极管1 | 光电二极管2 | 光电二极管3 |
| 固定权重 | 3.18 | 3.46 | 1.53 |
| 残差模值 | 2.97 | 2.95 | 1.46 |
| Huber模值 | 2.90 | 2.88 | 1.48 |
| 本文方法 | 2.24 | 2.79 | 1.13 |
表选项
 |
| 图 5 光电二极管1、2和3的偏置估计误差 Fig. 5 Bias estimation errors of photodiode 1, 2 and 3 |
| 图选项 |
5 结论 针对地球反照光的影响，本文建立了光电二极管动态偏置混合高斯量测模型，简单有效，通过对地球反照光的准确估计，提高了光电二极管对太阳矢量的测量精度。实验表明：
1) 应用地球反照光校正的光电二极管和地球敏感器组合定姿，可以消弱地球地球反照光的干扰，快速提高三轴姿态精度，精度可以达到0.2°~0.3°。
2) 在应用滑窗估计和随机权重估计量测模型参数过程中，采用多比例因子和Huber影响函数的权重处理方法，可以有效提高地球反照光动态偏置电压的估计精度。
本文提出的地球反照光校正方法，可以推广应用于光电二极管与其他姿态传感器(如地磁)的组合定姿。
参考文献
| [1] | POST M A, LI J Q, LEE R. A low-cost photodiode sun sensor for CubeSat and planetary microrover[J]. International Journal of Aerospace Engineering, 2013, 2013: 549080. |
| [2] | UNHELKAR V V, HABLANI H B. Spacecraft attitude determination with sun sensors, horizon sensors and gyros:Comparison of steady-state Kalman filter and extended Kalman filter[M]. Advances in Estimation, Navigation, and Spacecraft Control.Berlin: Springer, 2015: 413-437. |
| [3] | GARCIA R V, KUGA H K, ZANARDI M C F P S. Unscented Kalman filter applied to the spacecraft attitude estimation with euler angles[J]. Mathematical Problems in Engineering, 2012, 2012: 985429. |
| [4] | FLATLEY T W, MOORE W A.An earth albedo model: A mathematical model for the radiant energy input to an orbiting spacecraft due to the diffuse reflectance of solar radiation from the Earth below: NASA-TM-104596[R].Washington, D.C., NASA, 1994. |
| [5] | APPEL P. Attitude estimation from magnetometer and earth-albedo-corrected coarse sun sensor measurements[J]. Acta Astronautica, 2005, 56(1-2): 115-126. DOI:10.1016/j.actaastro.2004.09.001 |
| [6] | 韩柯, 金仲和, 王昊. 基于太阳能电池板的皮卫星最优姿态确定算法[J]. 浙江大学学报(工学版), 2010, 44(9): 1719-1723. HAN K, JIN Z H, WANG H. Optimal attitude determination method forpico-atellite using solar panels[J]. Journal of Zhejiang University(Engineering Science), 2010, 44(9): 1719-1723. DOI:10.3785/j.issn.1008-973X.2010.09.015 (in Chinese) |
| [7] | BHANDERI D D V. Spacecraft attitude determination with earth albedo corretted sun sensor measurements[M]. Aalborg: Aalborg University, 2005. |
| [8] | HAAVE H R.Simulating sun vector estimation and finding gyroscopes for the NUTS project[D].Trondheim: Norwegian University of Science and Technology, 2016. |
| [9] | JULIER S J, UHLMANN J K. Unscented filtering and nonlinear estimation[J]. Proceedings of the IEEE, 2004, 92(3): 401-422. DOI:10.1109/JPROC.2003.823141 |
| [10] | O'KEEFE S A, SCHAUB H. Consider-filter-based on-orbit coarse sun sensor calibration sensitivity[J]. Journal of Guidance, Control, and Dynamics, 2016, 40(5): 1300-1303. |
| [11] | GAO S, HU G, ZHONG Y. Windowing and random weighting-based adaptive unscented Kalman filter[J]. International Journal of Adaptive Control and Signal Processing, 2015, 29(2): 201-223. DOI:10.1002/acs.v29.2 |
| [12] | MASRELIEZ C, MARTIN R. Robust Bayesian estimation for the linear model and robustifying the Kalman filter[J]. IEEE transactions on Automatic Control, 1977, 22(3): 361-371. DOI:10.1109/TAC.1977.1101538 |
| [13] | 章仁为. 卫星轨道姿态动力学与控制[M]. 北京: 北京航空航天大学出版社, 1998: 157-176. ZHANG R W. Satellite dynamics and control of spacecraft[M]. Beijing: Beihang University Press, 1998: 157-176. (in Chinese) |
| [14] | KIM S G, CRASSIDIS J L, CHENG Y, et al. Kalman filtering for relative spacecraft attitude and position estimation[J]. Journal of Guidance, Control, and Dynamics, 2007, 30(1): 133-143. DOI:10.2514/1.22377 |
| [15] | FARRENKOPF R L. Analytic steady-state accuracy solutions for two common spacecraft attitude estimators[J]. Journal of Guidance, Control, and Dynamics, 1978, 1(4): 282-284. DOI:10.2514/3.55779 |
| [16] | ZANETTI R, DEMARS K J. Fully multiplicative unscented kalman filter for attitude estimation[J]. Journal of Guidance, Control, and Dynamics, 2017, 41(5): 1183-1189. |
| [17] | 秦永元, 张洪钺, 汪叔华. 卡尔曼滤波与组合导航原理[M]. 西安: 西北工业大学出版社, 2015: 33-57. QIN Q Y, ZHANG H Y, WANG S H. Kalman filtering and integrated navigation principle[M]. Xi'an: Northwestern Polytechnical University Press, 2015: 33-57. (in Chinese) |
| [18] | MARKLEY F L, MORTARI D. Quaternion attitude estimation using vector observations[J]. Journal of the Astronautical Sciences, 2000, 48(2): 359-380. |
