CN108827279A - 一种基于地球反照的模拟式太阳敏感器抗干扰方法 - Google Patents

Abstract

本发明提供了一种基于地球反照的模拟式太阳敏感器抗干扰方法，涉及属于敏感器抗干扰领域，本发明建立地球反照率模型，计算太阳敏感器接收的辐照度，并由此得到辐照度统计特性，通过获取太阳光与地球反照光判别函数，即可作为太阳敏感器对地球反照的抗干扰依据，进行太阳光和地球反照光的区分。本发明针对基于地球反照的太阳敏感器抗干扰问题，不仅分析了地球反照情况下辐照度与入射角之间的关系，而且设计线性分类器，确定太阳光与地球反照光之间的线性判别函数，能够充分利用测量信息提高抗干扰的正确率，本发明的分析结果表明该抗干扰方法的合理性，能为太阳敏感器测量量的正确使用提供参考依据。

Description

一种基于地球反照的模拟式太阳敏感器抗干扰方法

技术领域

本发明涉及属于敏感器抗干扰领域，尤其是一种太阳敏感器的抗干扰方法，防止将地球反照光误判为太阳光。

背景技术

在微小卫星姿态确定与控制系统(ADCS)中，太阳敏感器作为主要姿态测量部件，在使用过程中，难以避免地受到地球反照光的干扰，可能将其误认为太阳光进行计算，造成姿态确定系统的错误输入，导致星上太阳能板、通讯天线等器件指向存在偏差，从而影响卫星的任务执行甚至对其生存构成威胁。因此，要使用太阳敏感器进行姿态确定，必须增强其对地球反照光的抗干扰能力，消除或者减少地球反照光对太阳敏感器姿态测量的影响以保证姿态确定系统的正确输入，对于卫星控制系统非常重要。

针对地球反照光的抗干扰问题，目前已经对地球反照光的影响机理进行详细的分析，(杨阳，宏力，陆敬辉，姜伟，地气光对星敏感器星提取精度影响分析[J].光电工程，2016.4,43(4):8-14)，建立了地球反照光对敏感器的辐照度影响模型，可以用来进行太阳敏感器的抗地球反照干扰分析。

现有方法主要集中在硬件方面，在太阳敏感器的探测单元上增加光谱滤光器或者遮光罩(张春明，解永春，王立，赵春晖，钟红军。地球反照对星敏感器的影响分析[J].激光与红外,2012.9,42(9):1011-1015)，减小接受段的光谱范围或者降低进入敏感器的辐照度，消除杂光光谱能量，这种方法可以有效减小杂光影响，但会增加太阳敏感器的机械复杂度，提高器件成本。

而对于软件层面的研究内容较少(韩柯,金仲和,王昊.基于太阳能电池板的皮卫星最优姿态确定算法[J].太阳能学报,2010,44(9):1719-1723；王恩宏,胡以华,刘伟.数字式太阳姿态敏感器抗干扰特性研究[J].红外技术,2007.4,29(4):218-221)，主要通过仿真分析，根据统计结果得到电流固定阈值，尽量减小入射的地球反照光与太阳光亮度比值，但是由于入射光比值的不确定性，固定阈值没有足够的适应性，若电流阈值选择过小，当地球反照光入射角减小，接收到的地球反照光会增强，可能导致地球反照光产生的电流大小超过阈值，同时，若电流阈值选择过大，过设计会导致敏感器的使用次数减少，降低姿态确定的速度和精度，因此固定阈值法不能完全有效的进行地球反照光的抗干扰。

现有方法虽然进行了基本的地球反照抗干扰研究，但仅从单独的电流角度或者视场角度考虑，而且阈值限于固定，不能适应不同入射角下的地球反照光情况。

发明内容

为了克服现有技术的不足，本发明针对地球反照光的抗干扰方法设计，也就是根据太阳敏感器接收辐照度和入射角之间关系确定太阳光与地球反照光之间的线性判别函数关系，避免将地球误认为太阳的情况，为姿态确定系统提供正确的太阳矢量输入。相对于当前已有技术，本发明所设计的线性判别函数不同于固定阈值，有效增加太阳敏感器的正确使用情况次数，并且能够根据轨道情况调整线性判别函数，保证在不同轨道情况下的适应性。

本发明解决其技术问题所采用的技术方案的详细步骤如下：

步骤1：建立地球反照率模型

将地球按照纬度间隔Δφ，经度间隔Δθ进行网格划分，根据全年全球反照率的实测数据，计算各数据点的全年反照率平均值，作为地球反照率模型；

详细步骤如下：

将地球按照纬度间隔Δφ，经度间隔Δθ进行网格划分，获得l×p个数据点，对应2维的数据空间Φ×Θ：

式中,Φ为数据点的纬度集合，Θ为数据点的经度集合，φ_a为数据点的纬度，且φ_a∈Φ，a∈{1,2,...l}，θ_b为数据点经度，且θ_b∈Θ，b∈{1,2,...p}；

对每一数据点统计n天的反照率数据，进行平均值的计算：

式中，ρ_j(φ_a,θ_b)为数据点(φ_a,θ_b)在第j天的反照率，j∈{1,2,…n}，ρ(φ_a,θ_b)为数据点(φ_a,θ_b)的平均反照率；

步骤2：计算太阳敏感器接收的辐照度

根据黑体法则和朗伯体假设，每个数据点按照φ±Δφ/2和θ±Δθ/2扩展为微元，得到从单一微元到太阳敏感器的地球反照辐射通量密度，将地球表面光照区和太阳敏感器视场的交集区域中所有微元相加，计算卫星接收的总体地球反照光辐照度和太阳光辐照度；

详细步骤如下：

步骤2.1：地球反照表面区域为光照区与敏感器视场的交集，利用以下计算公式判断微元是否同时处于光照区和敏感器视场内：

式中，α_sun为太阳光在地表某点上的入射角，α_sensor为反照光在太阳敏感器上的入射角，r_gsat为地表某点指向卫星的距离矢量，L_s为太阳光线的单位矢量，n_g为地表某点法线单位矢量，n_sensor为太阳敏感器法线单位矢量；

当α_sun＜90°，微元处于光照区；当α_sensor＜FOV且α_sensor＜90°，微元处于太阳敏感器视场内，其中FOV为太阳敏感器半视场角；

步骤2.2：计算太阳敏感器接收的反照光辐照度

每个数据点按照φ±Δφ/2和θ±Δθ/2扩展为微元，微元面积A_c为：

式中：R为地球平均半径；

从某一微元到太阳敏感器的地球反照辐射通量密度E_valid为：

其中，太阳平均辐照度E_AM0＝1366.5W/m²，α_sat为太阳光在地表某点上的漫反射角；

根据是否处于光照区与敏感器视场区域，可得到太阳敏感器接收的反照光辐照度E_rft：

式中，V_sat为太阳敏感器视场区域，V_sun为地球表面光照区域；

步骤2.3：计算接收的太阳光辐照度

根据太阳位置、敏感器位置和太阳光入射角大小，计算接收的太阳光辐照度E_sun：

其中，β_sensor为太阳光直射在太阳敏感器上的入射角；

步骤2.4：计算太阳敏感器接收的总辐照度E、入射角γ和太阳指示标志σ，计算公式如下：

步骤3：辐照度统计特性

利用蒙特卡洛法，随机生成太阳矢量、卫星位置以及太阳敏感器方向，根据太阳敏感器给定数值设定太阳敏感器视场FOV，根据卫星运行轨道设定轨道高度，根据步骤2计算的球反照光辐照度和太阳光辐照度，继而获得每一次蒙特卡洛试验的总辐照度、入射角和太阳指示标志；将总辐照度和入射角画出x-y散点图中，即得到总辐照度与入射角之间的统计关系图；

步骤4：获取太阳光与地球反照光判别函数

根据步骤3的总辐照度与入射角之间的统计关系图，依据太阳指示标志σ_k的取值判断样本点是否为太阳光，按照样本点是否为太阳光将样本点划分为两类，将包含总辐照度和入射角所有样本点设为训练样本，引入最小均方误差方法，获得线性判别函数的权向量，将权向量代入判别函数，获得线性分类器，即可作为太阳敏感器对地球反照的抗干扰依据，进行太阳光和地球反照光的区分。

所述步骤3中利用蒙特卡洛试验计算总辐照度、入射角和太阳指示标志的详细步骤如下：

步骤3.1：在第k次蒙特卡洛试验时，随机生成太阳矢量、太阳敏感器位置、太阳敏感器方向，并设定太阳敏感器FOV、轨道高度h、轨道倾角i；

在[0,1]之间按照均匀分布生成随机数m_azi,k和m_ele,k，根据黄赤交角θ_ooe，生成太阳矢量的方位角和高低角

得到太阳光线单位矢量L_s,k为：

在[0,1]之间按照均匀分布生成二次随机数m_a′_zi,k和m_e′_le,k，并根据轨道倾角i，生成太阳敏感器位置的方位角和高低角：

得到太阳敏感器位置r_esat,k为：

其中，h是太阳敏感器到地面的高度，继而得到r_gsat,k＝r_esat,k-R·n_g，其中r_gsat,k为第k次蒙特卡洛试验时地表某点指向卫星的距离矢量；

在[0,1]之间按照均匀分布生成三次随机数m″_azi,k和m″_ele,k，生成太阳敏感器指向的方位角和高低角

得到太阳敏感器指向矢量n_sensor,k为：

步骤3.2：根据式(6)计算太阳敏感器接收的反照光辐照度E_rft,k，根据式(7)计算太阳敏感器接收的太阳光辐照度E_sun,k，并利用式(8)获得第k次总辐照度E_k、入射角γ_k和太阳指示标志σ_k；

重复步骤3.1和3.2，直到试验次数达到设定的蒙特卡洛试验次数L为止，即试验次数k＝1,2，…L，得到相应的第k次蒙特卡洛试验的总辐照度为E_k、入射角为γ_k和太阳指示标志为σ_k，其中太阳指示标志σ_k表征此样本点是否为太阳光，总辐照度E_k和入射角γ_k组成一个样本点，进而获得所有总辐照度与所有入射角的统计关系图；

所述步骤4的实现步骤如下：

利用LMSE(least mean-square error，最小均方误差)算法，选取判别函数的形式，根据LMSE算法获得判别函数的权向量，设定LMSE算法的迭代次数为t，设太阳光照与地球反照光的判别函数g为：

g＝α₁γ+α₂E+α₃ (11)

式中，α₁，α₂和α₃为待确定系数；

设权向量a及训练样本x为：

将所有样本点设为训练样本，其中，太阳指示标志σ_k＝1的训练样本为x_a，

x_b，……，训练样本相应的判别函数g＞0；太阳指示标志σ_k＝0的训练样本为x_a′，

x_b′，……，该训练样本相应的判别函数g＜0，且有：

x_a＝[γ_a,E_a,1]x_b＝[γ_b,E_b,1]…

x_a′＝[γ_a′,E_a′,1]x_b′＝[γ_b′,E_b′,1]…

其中，_γa,E_a分别为样本x_a对应的入射角和总辐照度，_γb,E_b分别为样本x_b对应的入射角和总辐照度；γ_a′,E_a′分别为样本x_a′对应的入射角和总辐照度；γ_b′,E_b′分别为样本x_b′对应的入射角和总辐照度；

步骤4.1：设正值矢量b的初值b₁＝[1 1…1 1…]^T，训练样本集X阵增广形式为：

步骤4.2：计算权向量a的权集初值a₁为：

a₁＝(X^TX)^-1X^Tb₁ (14)

步骤4.3：循环迭代计算，迭代公式如下：

其中，e_k是一个差量，a_k+1为第k+1次迭代的权向量，b_k+1为第k+1次迭代的正值矢量；

步骤4.4：当步骤4.3中公式(16)的迭代次数达到t次或e_k＝0，则获得最终权向量a和判别函数g；将待判定数据的总辐照度E_k、入射角γ_k代入判别函数g，如g＜0，则判定太阳敏感器将地球反照光误判为太阳光，否则敏感数据无效。

本发明的有益效果是针对基于地球反照的太阳敏感器抗干扰问题，不仅分析了地球反照情况下辐照度与入射角之间的关系，而且设计线性分类器，确定太阳光与地球反照光之间的线性判别函数，能够充分利用测量信息提高抗干扰的正确率，本发明的分析结果表明该抗干扰方法的合理性，能为太阳敏感器测量量的正确使用提供参考依据。

附图说明

图1为地球反照示意图。

图2为模拟式太阳敏感器示意图。

图3为全球反照率全年平均值。

图4为辐照度与入射角关系图。

图5为太阳光与地球反照光的判别函数示意图。

其中，α_sun为太阳光在地表某点上的入射角，α_sat为太阳光在地表某点上的漫反射角；α_sensor为反照光在太阳敏感器上的入射角；β_sensor为太阳光直射在太阳敏感器上的入射角；r_gsat为地表某点指向卫星的距离矢量；L_s为太阳光线的单位矢量；n_g为地表某点法线单位矢量；n_sensor为太阳敏感器法线单位矢量；θ为光线入射角；为光线方位角；O,X,Y,Z分别为太阳敏感器坐标系原点以及X，Y，Z轴；I₁,I₂,I₃,I₄分别为四块太阳能板的电流信息；H为太阳敏感器的上下表面的高度差。

具体实施方式

下面结合附图和实施例对本发明进一步说明。

本发明进行的太阳敏感器的抗干扰工作是整个姿态确定与控制系统中的重要环节，解决的问题是基于地球反照光模型，根据一定输入条件，包括轨道高度、倾角等，利用线性可分几何分类法确定合适的辐照度与入射角之间的线性判别函数关系，进而规避地球反照光的干扰。

本发明在现有技术的基础上，建立了以年份为单位的地球反照模型，结合太阳敏感器的视场情况，分析了辐照度与入射角之间的关系，运用线性分类器，确定太阳光与地球反照光之间的线性判别函数关系，能比较合理地规避太阳敏感器的地球反照光干扰。

为实现上述目的，在本发明中采用的具体实施如下：

步骤1：将地球按照纬度间隔Δφ＝1°，经度间隔Δθ＝1.25°进行网格划分，利用美国戈达德空间飞行中心的TOMS全球反照率数据(ozoneaq.gsfc.nasa.gov/data/reflectivity)，根据全年全球反照率的实测数据，计算各数据点的全年反照率平均值，作为简化的地球反照率模型；

步骤1.1：网格划分

将地球按照纬度间隔Δφ＝1°，经度间隔Δθ＝1.25°进行网格划分，获得180×288个数据点，对应2维的数据空间Φ×Θ。

式中：

Φ——数据点的纬度集合；

Θ——数据点的经度集合；

φ_a——数据点纬度，φ_a∈Φ，a∈{1 2…180}；

θ_b——数据点经度，θ_b∈Θ，b∈{1 2…288}。

步骤1.2：计算平均值

利用美国戈达德空间飞行中心的TOMS全球反照率数据(ozoneaq.gsfc.nasa.gov/data/reflectivity)，对每一数据点统计n天的反照率数据，进行平均值的计算：

式中：ρ_j(φ_a,θ_b)为数据点(φ_a,θ_b)在第j天的反照率，j∈{1 2…n}；ρ(φ_a,θ_b)为数据点(φ_a,θ_b)的平均反照率；

步骤2：计算太阳敏感器接收的辐照度

根据黑体法则和朗伯体假设，参考(Dan D.V.Bhanderi,Thomas Bak.Modelingearth albedo for satellites in earth[M].AIAA Guidance,Navigation and ControlConference and Exhibit,2005.8)，图1为地球反照示意图，考虑太阳光在太阳敏感器上的入射角，每个数据点按照φ±Δφ/2和θ±Δθ/2扩展为微元，得到从单一微元到太阳敏感器的地球反照辐射通量密度，将地球表面光照区和太阳敏感器视场的交集区域中所有微元相加，计算卫星接收的总体地球反照光辐照度和太阳光辐照度；

步骤2.1：确定地球反照有效区域

地球反照表面区域为光照区与敏感器视场的交集，所以需要利用以下计算公式判断微元是否同时满足光照区和敏感器视场内。

其中，α_sun——太阳光在地表某点上的入射角；α_sensor——反照光在太阳敏感器上的入射角；r_gsat——地表某点指向卫星的距离矢量；L_s——太阳光线的单位矢量；n_g——地表某点法线单位矢量；n_sensor——太阳敏感器法线单位矢量；

当α_sun＜90°，微元处于光照区；当α_sensor＜FOV且α_sensor＜90°，微元处于太阳敏感器视场内，其中FOV为太阳敏感器半视场角；

步骤2.2：计算太阳敏感器接收的反照光辐照度

每个数据点按照φ±Δφ/2和θ±Δθ/2扩展为微元，微元面积A_c为

式中：R为地球平均半径。

那么，从某一微元到太阳敏感器的地球反照辐射通量密度E_valid为：

其中：太阳平均辐照度E_AM0＝1366.5W/m²，α_sat为太阳光在地表某点上的漫反射角。

根据是否处于光照区与敏感器视场的交集分别讨论，可以得到太阳敏感器接收的反照光辐照度E_rft：

式中：V_sat为太阳敏感器视场区域；V_sun为地球表面光照区域。

步骤2.3：计算接收的太阳光辐照度

根据太阳位置、敏感器位置和太阳光入射角大小，计算接收的太阳光辐照度E_sun：

其中，β_sensor为太阳光直射在太阳敏感器上的入射角；

步骤2.4：计算太阳敏感器接收的总辐照度E、入射角γ和太阳指示标志σ：

步骤3：辐照度统计特性

利用蒙特卡洛法，随机生成太阳矢量、卫星位置以及太阳敏感器方向，根据太阳敏感器说明书给定的数值设定太阳敏感器视场FOV，根据卫星运行轨道设定轨道高度，根据步骤2计算地球反照光辐照度和太阳光辐照度，继而获得每一次的总辐照度、入射角和太阳指示标志；将每一次得到的总辐照度和入射角画入x-y散点图中，即得到总辐照度与入射角之间的统计关系图；

利用蒙特卡洛法，计算不同情况下太阳敏感器接收到的总辐照度和入射角，获得地球反照光统计特性。图2为模拟式太阳敏感器示意图。

具体步骤为：

步骤3.1：在第k次蒙特卡洛试验时，随机生成太阳矢量、太阳敏感器位置、太阳敏感器方向，并设定太阳敏感器FOV、轨道高度h、轨道倾角i；

在[0,1]之间按照均匀分布生成随机数m_azi,k和m_ele,k，并考虑黄赤交角θ_ooe，生成太阳矢量的方位角和高低角

可以得到太阳光线单位矢量L_s,k为：

在[0,1]之间按照均匀分布生成二次随机数m_a′_zi,k和m_e′_le,k，并考虑轨道倾角i，生成太阳敏感器位置的方位角和高低角

可以得到太阳敏感器位置r_esat,k为：

其中，h是太阳敏感器到地面的高度，继而得到r_gsat,k＝r_esat,k-R·n_g，r_gsat,k为第k次地表某点指向卫星的距离矢量；

在[0,1]之间按照均匀分布生成三次随机数m″_azi,k和m″_ele,k，生成太阳敏感器指向的方位角和高低角

可以得到太阳敏感器指向矢量n_sensor,k为：

步骤3.2：根据式(6)计算太阳敏感器接收的反照光辐照度E_rft,k，根据式(7)计算太阳敏感器接收的太阳光辐照度E_sun,k，并利用式(8)获得第k次总辐照度E_k、入射角γ_k和太阳指示标志σ_k。

重复步骤3.1和3.2，直到试验次数达到设定的蒙特卡洛试验次数L为止，即试验次数k＝1,2，…L，得到相应的第k次蒙特卡洛试验的总辐照度为E_k、入射角为γ_k和太阳指示标志为σ_k，其中太阳指示标志σ_k表征此样本点是否为太阳光，总辐照度E_k和入射角γ_k组成一个样本点，进而获得所有总辐照度与所有入射角的统计关系图；

步骤4：获取太阳光与地球反照光判别函数

根据步骤3的总辐照度与入射角之间的统计关系图，依据太阳指示标志σ_k取值判断样本点是否为太阳光，按照样本点是否为太阳光将样本点划分为两类，将包含总辐照度和入射角所有样本点设为训练样本，引入最小均方误差方法，获得线性判别函数的权向量，将权向量代入判别函数，获得线性分类器，即可作为太阳敏感器对地球反照的抗干扰依据，进行太阳光和地球反照光的区分。

利用LMSE(least mean-square error，最小均方误差)算法设定判别函数的形式，确定线性判别函数的权向量，设定LMSE迭代次数为t，而对于判别函数的训练过程，就是确定该函数的权集，设太阳光照与地球反照光的判别函数g为：

g＝α₁γ+α₂E+α₃ (27)

式中，α₁，α₂和α₃为待确定系数。

设权向量a及训练样本x为：

将所有样本点设为训练样本，其中，太阳指示标志σ_k＝1的训练样本为x_a，

x_b，……，训练样本相应的判别函数g＞0；太阳指示标志σ_k＝0的训练样本为x_a′，

x_b′，……，该训练样本相应的判别函数g＜0，且有：

x_a＝[γ_a E_a 1]x_b＝[γ_b E_b 1]…

x_a′＝[γ_a′ E_a′ 1]x_b′＝[γ_b′ E_b′ 1]…

步骤4.1：设正值矢量b的初值b₁＝[1 1…1 1…]^T，训练样本集X阵增广形式为：

步骤4.2：计算权向量a权集初值a₁

a₁＝(X^TX)^-1X^Tb₁ (30)

步骤4.3：循环迭代计算，迭代公式如下：

其中，e_k是一个差量，a_k+1为第k+1次迭代的权向量，b_k+1为第k+1次迭代的正值矢量；

步骤4.4：当步骤4.3中公式(16)的迭代次数达到t次或e_k＝0，则获得最终权向量a和判别函数g；将待判定数据的总辐照度E_k、入射角γ_k代入判别函数g，当g＜0，则判定太阳敏感器将地球反照光误判为太阳光，否则敏感数据无效。

本实例中，以某低轨道卫星为例，在任意姿态情况下进行太阳光与地球反照光的判别，轨道信息如表1所示：

表1轨道信息

轨道偏心率	0
		轨道高度	300km
轨道倾角	40.901°
		太阳矢量	随机生成
卫星姿态	随机生成

a.建立地球反照率模型

①网格划分

将地球按照纬度Δφ_g＝1°，经度Δθ_g＝1.25°进行网格划分，得到的180×288个数据点，对应2维的数据空间D为

②计算全年平均值

按照同一微元，统计365天的反照率数据，进行平均值的计算

其中，在两极附近由于极夜存在，部分时间段数缺失，仅采用有效数据进行分析，其中n＜365。如：对于地球表面微元(经度-59.375°，纬度6.5°)，其中ρ(φ_g)＝0.2827。图3为全球反照率全年平均值。

b.计算太阳敏感器接收的辐照度

获得太阳敏感器相关数据，确定地球反照有效区域后，计算太阳敏感器接收的反照光辐照度和太阳光辐照度，最后计算太阳敏感器接收的总辐照度E_sum和入射角γ_sum。

以某次仿真情况为例：地心指向太阳敏感器单位矢量为[0.5558 -0.80290.2155]^T，地心指向太阳单位矢量为[0.9889 -0.1138 -0.0959]^T，当轨道高度较低时，太阳敏感器指向太阳单位矢量与之相同为[0.9889 -0.1138 -0.0959]^T，太阳敏感器法线指向[-0.4175 -0.7873 -0.4537]^T，太阳敏感器视场大小FOV＝60°、高度h＝300km、轨道倾角i＝40.901°。

①确定地球反照有效区域

计算可得：地球表面微元(经度-59.375°，纬度6.5°)，对应α_sun＝54.0570°＜90°，α_sensor＝46.3905°＜60°，说明微元处于反照有效区域内，并计算其他符合条件的微元。

②计算太阳敏感器接收的反照光辐照度

以此微元为例进行计算示例，地球表面微元(经度-59.375°，纬度6.5°)，每个数据点按照φ_g±0.5°和θ_g±0.625°扩展为微元，扩展后对应面积为

A_c＝15496.1559km²

对应反照率为0.2827，α_sat＝7.1746°，太阳敏感器指向微元的单位矢量为[-0.5555 -0.1063 -0.8247]^T。那么，从单一微元到太阳敏感器的地球反照辐射通量密度E_valid为：

此微元指向太阳敏感器的通量密度矢量E_valid为：

E_valid＝[0.5555 0.1063 0.8247]^T×0.0736E_AM0

将所有同时处于光照区与敏感器视场的微元进行矢量加和，得到：

E_rft＝[0.1496 0.0335 0.1474]^TE_AM0

③计算太阳敏感器接收的太阳光辐照度为0；

④计算太阳敏感器接收的总辐照度E_sum和入射角γ_sum和太阳指示标志σ：

c.计算辐照度统计特性

①按照上述条件随机生成太阳矢量、太阳敏感器位置、太阳敏感器方向，输入太阳敏感器视场大小FOV＝60°、轨道高度h＝300km、轨道倾角i＝40.901°。

例如：按照均匀分布生成[0,1]间的随机数m_azi,1＝0.4和m_ele,1＝0.1，并考虑黄赤交角θ_ooe＝23.4397°，利用其生成太阳矢量的方位角和高低角。

可以得到太阳光线单位矢量L_s为：

L_s＝[0.7661 -0.2601 -0.5878]

随机生成固定高度和倾角的太阳敏感器位置、利用本体系与北东地坐标系随机生成太阳敏感器方向的生成方法相似。

②根据式(6)计算太阳敏感器接收的反照光辐照度E_rft,k，根据式(7)计算太阳敏感器接收的太阳光辐照度E_sun,k，并利用式(8)获得总辐照度E_k、入射角γ_k和太阳指示标志σ_k。

③设定的蒙特卡洛试验次数为10000，获得E-γ关系图如图4所示。按照数据记录，辐照度较大部分为太阳光，辐照度较小部分为地球反照光，具备线性可分性。

d.获取太阳光与地球反照光判别函数

从图4中，选择训练样本，太阳光照与地球反照光的判别函数为：g＝α₁γ+α₂E+α₃，选择训练样本为x_a，x_b，…，x_a′，x_b′，…

σ_k＝1:x_a＝[1.7504 0.9995 1]x_b＝[20.1980 0.9385 1]…

σ_k＝0:x_a′＝[1.3041 0.4411 1]x_b′＝[10.1216 0.4909 1]…

①设正值矢量b的初值b₁＝[1 1…1 1…]^T，训练样本集X阵增广形式为：

②计算权向量a权集初值a₁：

a₁＝(X^TX)^-1X^Tb₁

＝[0.0303 3.7436 -2.9583]^T

③循环计算E_k和A_k+1、b_k+1。

④循环100次结果为：

a₁₀₀＝[0.0042 2.0246 -1.2500]^T

获得的判别函数为：

g＝0.0042γ+2.0246E-1.2500

太阳光与地球反照光的判别函数的判别效果如图5所示，可以看出，可以有效区分太阳光和地球反照光。若需电流与入射角关系判别函数，可根据敏感器特性对辐照度转换，即可得到。

Claims (3)

1.一种基于地球反照的模拟式太阳敏感器抗干扰方法，其特征在于包括下述步骤：

步骤1：建立地球反照率模型

将地球按照纬度间隔Δφ，经度间隔Δθ进行网格划分，根据全年全球反照率的实测数据，计算各数据点的全年反照率平均值，作为地球反照率模型；

详细步骤如下：

将地球按照纬度间隔Δφ，经度间隔Δθ进行网格划分，获得l×p个数据点，对应2维的数据空间Φ×Θ：

式中,Φ为数据点的纬度集合，Θ为数据点的经度集合，φ_a为数据点的纬度，且φ_a∈Φ，a∈{1,2,...l}，θ_b为数据点经度，且θ_b∈Θ，b∈{1,2,...p}；

对每一数据点统计n天的反照率数据，进行平均值的计算：

式中，ρ_j(φ_a,θ_b)为数据点(φ_a,θ_b)在第j天的反照率，j∈{1,2,…n}，ρ(φ_a,θ_b)为数据点(φ_a,θ_b)的平均反照率；

步骤2：计算太阳敏感器接收的辐照度

根据黑体法则和朗伯体假设，每个数据点按照φ±Δφ/2和θ±Δθ/2扩展为微元，得到从单一微元到太阳敏感器的地球反照辐射通量密度，将地球表面光照区和太阳敏感器视场的交集区域中所有微元相加，计算卫星接收的总体地球反照光辐照度和太阳光辐照度；

详细步骤如下：

步骤2.1：地球反照表面区域为光照区与敏感器视场的交集，利用以下计算公式判断微元是否同时处于光照区和敏感器视场内：

式中，α_sun为太阳光在地表某点上的入射角，α_sensor为反照光在太阳敏感器上的入射角，r_gsat为地表某点指向卫星的距离矢量，L_s为太阳光线的单位矢量，n_g为地表某点法线单位矢量，n_sensor为太阳敏感器法线单位矢量；

当α_sun＜90°，微元处于光照区；当α_sensor＜FOV且α_sensor＜90°，微元处于太阳敏感器视场内，其中FOV为太阳敏感器半视场角；

步骤2.2：计算太阳敏感器接收的反照光辐照度

每个数据点按照φ±Δφ/2和θ±Δθ/2扩展为微元，微元面积A_c为：

式中：R为地球平均半径；

从某一微元到太阳敏感器的地球反照辐射通量密度E_valid为：

其中，太阳平均辐照度E_AM0＝1366.5W/m²，α_sat为太阳光在地表某点上的漫反射角；

根据是否处于光照区与敏感器视场区域，可得到太阳敏感器接收的反照光辐照度E_rft：

式中，V_sat为太阳敏感器视场区域，V_sun为地球表面光照区域；

步骤2.3：计算接收的太阳光辐照度

根据太阳位置、敏感器位置和太阳光入射角大小，计算接收的太阳光辐照度E_sun：

其中，β_sensor为太阳光直射在太阳敏感器上的入射角；

步骤2.4：计算太阳敏感器接收的总辐照度E、入射角γ和太阳指示标志σ，计算公式如下：

步骤3：辐照度统计特性

利用蒙特卡洛法，随机生成太阳矢量、卫星位置以及太阳敏感器方向，根据太阳敏感器给定数值设定太阳敏感器视场FOV，根据卫星运行轨道设定轨道高度，根据步骤2计算的球反照光辐照度和太阳光辐照度，继而获得每一次蒙特卡洛试验的总辐照度、入射角和太阳指示标志；将总辐照度和入射角画出x-y散点图中，即得到总辐照度与入射角之间的统计关系图；

步骤4：获取太阳光与地球反照光判别函数

根据步骤3的总辐照度与入射角之间的统计关系图，依据太阳指示标志σ_k的取值判断样本点是否为太阳光，按照样本点是否为太阳光将样本点划分为两类，将包含总辐照度和入射角所有样本点设为训练样本，引入最小均方误差方法，获得线性判别函数的权向量，将权向量代入判别函数，获得线性分类器，即可作为太阳敏感器对地球反照的抗干扰依据，进行太阳光和地球反照光的区分。

2.根据权利要求1所述的一种基于地球反照的模拟式太阳敏感器抗干扰方法，其特征在于：

所述步骤3中利用蒙特卡洛试验计算总辐照度、入射角和太阳指示标志；

步骤3.1：在第k次蒙特卡洛试验时，随机生成太阳矢量、太阳敏感器位置、太阳敏感器方向，并设定太阳敏感器FOV、轨道高度h、轨道倾角i；

在[0,1]之间按照均匀分布生成随机数m_azi,k和m_ele,k，根据黄赤交角θ_ooe，生成太阳矢量的方位角和高低角

得到太阳光线单位矢量L_s,k为：

在[0,1]之间按照均匀分布生成二次随机数m′_azi,k和m′_ele,k，并根据轨道倾角i，生成太阳敏感器位置的方位角和高低角

得到太阳敏感器位置r_esat,k为：

其中，h是太阳敏感器到地面的高度，继而得到r_gsat,k＝r_esat,k-R·n_g，其中r_gsat,k为第k次蒙特卡洛试验时地表某点指向卫星的距离矢量；

在[0,1]之间按照均匀分布生成三次随机数m″_azi,k和m″_ele,k，生成太阳敏感器指向的方位角和高低角

得到太阳敏感器指向矢量n_sensor,k为：

步骤3.2：根据式(6)计算太阳敏感器接收的反照光辐照度E_rft,k，根据式(7)计算太阳敏感器接收的太阳光辐照度E_sun,k，并利用式(8)获得第k次总辐照度E_k、入射角γ_k和太阳指示标志σ_k；

重复步骤3.1和3.2，直到试验次数达到设定的蒙特卡洛试验次数L为止，即试验次数k＝1,2，…L，得到相应的第k次蒙特卡洛试验的总辐照度为E_k、入射角为γ_k和太阳指示标志为σ_k，其中太阳指示标志σ_k表征此样本点是否为太阳光，总辐照度E_k和入射角γ_k组成一个样本点，进而获得所有总辐照度与所有入射角的统计关系图。

3.根据权利要求1所述的一种基于地球反照的模拟式太阳敏感器抗干扰方法，其特征在于：

所述步骤4中，利用LMSE算法，选取判别函数的形式，根据LMSE算法获得判别函数的权向量，设定LMSE算法的迭代次数为t，设太阳光照与地球反照光的判别函数g为：

g＝α₁γ+α₂E+α₃ (11)

式中，α₁，α₂和α₃为待确定系数；

设权向量a及训练样本x为：

将所有样本点设为训练样本，其中，太阳指示标志σ_k＝1的训练样本为x_a，x_b，……，训练样本相应的判别函数g＞0；太阳指示标志σ_k＝0的训练样本为x_a′，x_b′，……，该训练样本相应的判别函数g＜0，且有：

x_a＝[γ_a,E_a,1]x_b＝[γ_b,E_b,1]…

x_a′＝[γ_a′,E_a′,1]x_b′＝[γ_b′,E_b′,1]…

其中，γ_a,E_a分别为样本x_a对应的入射角和总辐照度，γ_b,E_b分别为样本x_b对应的入射角和总辐照度；γ_a′,E_a′分别为样本x_a′对应的入射角和总辐照度；γ_b′,E_b′分别为样本x_b′对应的入射角和总辐照度；

步骤4.1：设正值矢量b的初值b₁＝[1 1 … 1 1 …]^T，训练样本集X阵增广形式为：

步骤4.2：计算权向量a的权集初值a₁为：

a₁＝(X^TX)^-1X^Tb₁ (14)

步骤4.3：循环迭代计算，迭代公式如下：

其中，e_k是一个差量，a_k+1为第k+1次迭代的权向量，b_k+1为第k+1次迭代的正值矢量；

步骤4.4：当步骤4.3中公式(16)的迭代次数达到t次或e_k＝0，则获得最终权向量a和判别函数g；将待判定数据的总辐照度E_k、入射角γ_k代入判别函数g，如g＜0，则判定太阳敏感器将地球反照光误判为太阳光，否则敏感数据无效。