一种用于静态红外地球敏感器的高效算法

陈 璐，王 昊，王本冬，金仲和

(浙江大学航空航天学院，杭州 310027)

0 引 言

近年来，微纳卫星技术得到了迅猛发展。它有着质量轻、体积小、功耗小的特点[1-2]，在航天应用领域极具潜力。目前，微纳卫星所执行的任务种类和难度在不断增加，例如要求微纳卫星实现大角度快速机动、编队飞行、双目视觉等[3-5]。而提升微纳卫星的机动能力是实现高难度任务的基础。

微纳卫星机动能力取决于姿态敏感器和控制器的性能。姿态敏感器包括地球敏感器、太阳敏感器、星敏感器和陀螺仪等[6-7]。其中静态红外地球敏感器是微纳卫星较为常用且可靠的定姿敏感器[8-12]。

静态红外地球敏感器通过对地球凝视拍摄，计算成像矢量的方向，经过坐标转换，最终得到卫星的对地姿态[9]。它相较于动态扫描式红外地球敏感器而言，具有速度快、功耗小、寿命长和体积小等优势。文献[10]中静态红外地球敏感器的姿态确定算法是将最小二乘法与霍夫变换结合，通过计算地球圆盘圆心来确定姿态角，姿态确定精度为0.1°。文献[11] 也利用最小二乘迭代的方式求解非整圆弧的圆心。文献[12]考虑了地球扁率和红外辐射曲线随季节变化的特点，利用补偿和滤波的方式获得地球切线边缘。但以上算法都基于地球圆盘圆心的计算，这会占用较多计算资源，同时地球圆盘圆心坐标的微小偏差也会使姿态角计算产生明显偏差。

本文提出了一种基于静态全景红外地球敏感器的微纳卫星姿态角确定算法。新算法采用间隔向量积与t Location-Scale分布置信区间均值结合的算法。相较于传统的方法，新算法直接计算地球圆环的法向量，避免了计算圆环圆心位置的中间步骤，使得定姿精度提高50%，定姿速度提高87%，提升了微纳卫星的姿态确定性能。

1 算法分析

1.1 测量模型

静态红外地球敏感器的测量目标如图1所示，地心矢量与成像轴Z轴的夹角β为测定的姿态角。

图1 静态红外地球敏感器的成像示意图Fig.1 Schematic diagram of static infrared earth sensor imaging

敏感器采用全景式鱼眼红外镜头，其成像规律为[13]：

y=fθ

(1)

式中：y为入射光线所成像素点与成像中心的距离，f为焦距，θ为入射光线与镜头法线夹角。

红外敏感器的成像与θ正相关，故无论敏感器转到什么方向，地球成像面均为E面。将E面投射到光学系统成像面，其坐标为(X，Y)，系统算法的目的是通过坐标(X，Y)得到β值。

1.2 成像点预处理

1.2.1边缘腐蚀提取

由于地球红外成像是一个圆面，后续计算是基于成像边缘进行的，故首先要进行边缘提取操作。这里采用Sobel算子强边沿检测和腐蚀结合的方式。具体而言，Sobel算子图像领域为

(2)

式中：z1,z2,z3,…,z9为该位置像素点的灰度值。

计算梯度如下

[(z3+2z6+z9)-(z1+2z4+z7)]2}1/2

(3)

(4)

二值图像腐蚀即对图像A采用向量减法运算

AΘB={(X,Y)|B(X,Y)⊆A}

(5)

式中：B(X,Y)为位于(X，Y)的结构元素像素区域，运算AΘB表示:与A相应位置的像素值一致时，A相应位置的像素点保留[15]。

1.2.2解畸变

基于全景鱼眼红外镜头的成像特性[13]，将成像点(X，Y)映射至半径为R的球面上，得到的边沿点(x,y,z)构成一个截面圆环D，如图2所示。

根据镜头成像公式(1)有

(6)

图2 像素点解畸变Fig.2 Elimination of pixel distortion

根据物理成像关系，可以得到成像位置坐标与红外成像三维方向之间的关系:

(7)

(8)

z=Rcos(θ)

(9)

式中：R为映射球面半径，点(x,y,z)为光学成像点(X，Y)到球面的映射，它们构成一个空间截面圆环D。

1.3 目标角计算

当得到成像圆环后，通常利用圆心位置计算目标角度，方法包括：三点均值法，霍夫变换+最小二乘法和面积积分法[16]。本文提出间隔向量积与t Location-Scale分布置信区间均值结合的方法，与原有霍夫变换+最小二乘算法[10]对比，如图3所示。

图3 算法流程图对比Fig.3 Comparison of algorithm flow charts

间隔向量积和t Location-Scale分布拟合是新算法与传统算法的本质差别。一方面，间隔向量积采用直接计算圆环法向量的方法，避免了计算圆环圆心的中间步骤，也避免了霍夫变换对切半径的估计和全局计算。减少了关键信息的损耗，节省计算资源。另一方面，t Location-Scale分布拟合可用于总体方差未知的情况，能较好地排除异常值，模拟尾部较重的数据分布和集中趋势。由t Location-Scale分布拟合滤去干扰点，可以获得准确的姿态角。

1.3.1间隔向量积

经过解畸变处理后，光学系统成像平面上的像素点被投影到空间球面上，构成截面圆环D,如图4所示。圆环D上每个像素点有一个三维坐标(x,y,z)。然而由于边缘提取的误差问题，有些像素点可能不在截面圆上，它们被称为干扰点。

图4 像素点映射至空间位置Fig.4 Pixel mapping spatial position

任意选取截面圆环D上三点a，b和c，可以构成两个向量Vba和Vbc。由向量定理可知[15]，两向量Vba和Vbc矢积得到的向量Vβ与原两向量所在平面D垂直。则有Vβ与r平行，显然Vβ与Z轴所成夹角即为目标值β。

具体而言，若将a,b,c三点选取的距离较近，其构成的两向量将趋于共线，矢积所得垂向量趋于0，不便于后续计算，所以应间隔选取像点。将空间圆环D分为A段、B段和C段，表示为三个集合A={a1,a2,…,a3},B={b1,b2,…,b3}和C={c1,c2,…,c3}。各部分包含n个像点，如图5所示。

图5 分割截面圆环D上像点Fig.5 Partition of image points on ring D

从这三个集合中选取三个点ak,bk和ck，得到两个向量Vbak和Vbck，称为空间向量对。顺序取点，至少可以构成n个空间矢量对，它们在同一空间圆环D上。由于向量Vbak和向量Vbck共面不共线，根据向量理论[17]，得到

Vβ={Vβk|Vβk=Vabk×Vack,k=1,2,3,…,n}

(10)

式中：Vβk为向量积结果，Vβ是由n个Vβk向量构成的集合。

Vβ中向量Vβk均垂直于圆环D，故计算β角有

(11)

式中：Bt为单张红外成像得到的β角集合。

1.3.2t Location-Scale分布拟合

由于干扰点的存在，有些像点并不分布在同一平面D上，导致了Bt集合当中角度βk的差异。其典型的概率密度分布如图6所示。

图6 5°时Bt内样本概率密度分布Fig.6 Probability density distribution of samples in Bt at 5°

由图6可知，βk的分布存在偏度，即存在较重尾部。t Location-Scale分布可用于总体方差未知的情况，能够模拟出尾部较重的数据分布，拟合数据的集中趋势，同时排除异常值的干扰[18]。相比之下，受重尾的影响正态分布的均值容易偏离大多数样本所在的区域。所以选择t Location-Scale分布进行拟合，其概率密度函数如下

(12)

式中：Γ(·)表示Gamma方程，∞ <μ<+∞为位置参数，σ>0为尺度参数，ν>0为形状参数。

t Location-Scale分布的参数可以通过极大似然函数估计得到[19]则其极大似然函数有

(13)

对(μ,σ,ν) 求μ和σ的偏导数，令其等于0有

(14)

(15)

(16)

式(14)至式(16)表明，参数间存在耦合，但仍可以通过迭代收敛的方法得到估计值，其迭代式为

(17)

(18)

式中：(i)表示第i次迭代。

迭代终止条件为前后两次结果差值小于10-4。由于自由度ν未知，故进行极大似然估计的迭代

(19)

式中：φ(x)=dln(Γ(x))/dx,-φ(·)+ln(·)近似为[20]

(20)

由式(17)～式(20)可得t Location-Scale的μ和σ。

在具体工程应用中，实时通过t Location-Scale分布拟合，需占用较多计算资源，故简化计算。由于Bt∈(0,90)°，故将区间(0,90)划分为h个子区间(90(i-1)/h, 90i/h)，其中i=1,2,3,…,h，则估算μ有

(21)

地面试验可测得的t Location-Scale分布在各个角度区间的方差，由此建立方差经验表，如表1所示，实际工程中可以通过μ查表得到σ。由于t Location-Scale分布的方差结果主要受红外镜头固有性能的影响，在重复测试后，固定角度的方差估计结果没有大的波动。所以方差经验表的引入在±10°内对卫星姿态确定精度没有明显的影响，同时能减少迭代估计的运算量，从而加快运算速度。

1.3.3置信区间均值

由于Bt样本较易出现异常值，直接利用整体均值μ估计β值，会受到干扰。故选取一定置信区间内的样本均值进行估计。虽然置信区间较宽时，其置信水平较高，但是包含的干扰点较多。根据3σ准则可知，数据分布在(μ-σ,μ+σ)的概率为68.27%，是较为适中的接受域。最终估计β值为

(22)

式中：m为落入接受域中βk的个数。

2 试验测试

静态红外地球敏感器的测试系统，如图7所示。

图7 静态红外地球敏感器测试系统Fig.7 Static infrared earth sensor test system

图7中红外地球敏感器测试系统由5部分组成：红外地球模拟器、红外地球敏感器、单轴高精度转台、DSP和PC机。测试方案如下：

1)红外地球模拟器的直径为38 cm，镜头距离模拟地球的距离为3 cm，模拟950 km轨道高度。

2)PC机发出角度控制信号，控制转台转动。

3)将红外地球敏感器置于高精度转台上，转台带动红外地球敏感器转动，多角度拍摄地球模拟器。在(-10°，10°)范围内，按照 0.1°为间隔拍摄，在(-50°，-10°)和(10°，50°)范围，按照1°为间隔拍摄。得到不同角度下的地球成像。

4)DSP将成像数据返回给PC，PC机利用图像处理算法分析数据，得出地球矢量信息。

根据测试系统开展测试，测试场景如图8所示。

图8 测试场景Fig.8 Test scenario

3 测试结果分析

3.1 边缘提取

地球红外成像如图9所示，需将图像中间60°盲区去除。之后由Sobel边缘检测得到强边缘图像如图10所示。最后经过二值图像腐蚀，如图11所示。

图9 地球红外成像Fig.9 Earth infrared imaging

图10 Sobel边缘检测Fig.10 Sobel edge detection

图11 二值图像腐蚀Fig.11 Binary image corrosion

3.2 解畸变

式(6)～式(9)得到解畸变后的像点如图12所示。图中像点构成的圆环即为截面圆环D，像点疏密不均，边缘存在一些干扰点。

图12 解畸变效果Fig.12 Distortion elimination effect

3.3 间隔向量积与分布结果

对截面圆环D上的像点进行间隔向量积运算，以10°附近为例，可以得到运算结果如图13所示。

图13 10度的Bt样本Fig.13 Bt samples at 10°

之后通过t Location-Scale拟合，得到Bt样本概率密度分布和拟合结果。例如，当目标β值为60°，30°和5°时，结果如图14、图15和图6所示。

图14 60°Bt样本概率密度分布Fig.14 Probability density distribution of samples in Bt at 60°

从图14和图15可以看出，Bt样本分布呈现单峰特性，存在偏度，在60°附近，Bt样本分布右偏；在30°以内，Bt样本分布左偏。t Location-Scale分布估计的参数，如表1所示。由表1可知，随着β值的减小，其均值附近的概率密度显著增加，分布趋向于集中。之后计算在1σ区间内Bt样本均值，如图14所示。对每一张红外成像做如上处理，可以得到目标角度β的估计值。

表1 t Location-Scale分布参数估计Table 1 t Location-Scale distribution parameter estimation

值得注意的是，以下因素会对新算法造成影响。首先，t Location-Scale分布的形状参数过高，将使得分布接近于正态分布，减弱剔除干扰点的能力；反之分布曲线双侧尾部会过高，将减少算法的有效数据量。尺度参数过高会使计算包含太多干扰点；反之会减弱算法对位置参数不准确的容忍能力。其次，间隔向量积中像点距离不能随机选取，否则将可能使像点构成的两向量趋于共线，从而无法计算姿态角。另外，红外地球成像的轮廓不能过腐蚀，否则将损耗有效像点，影响算法的姿态确定精度。

3.4 算法结果比较

按照第2.2节所述的试验方法进行测试。将新算法与常用的霍夫变换+最小二乘法[10]的实测结果比较，如图16和图17所示。可知当测量角度大于40°时，新算法能够得到平滑的测试结果。

图16 原算法目标角度估计结果Fig.16 Target angle estimation results of original algorithm

图17 新算法目标角度估计结果Fig.17 Target angle estimation results of new algorithm

图18 ±10°内原算法测量误差Fig.18 Measuring error of original algorithm in ± 10°

图19 ±10°内新算法测量误差Fig.19 Measuring error of original algorithm in ± 10°

由图18和图19可知，当目标角度在±10°以内时，原算法的误差角度的标准差为0.06°，新算法为0.03°，测量精度提升了50%。进一步比较各区间内算法精度及计算速度如表2所示。

表2 在不同区间内算法的精度比较Table 2 Accuracy comparison of algorithms in different regions

由表2可知，新算法在各区间的测量精度均有所提高。另外，表2中列出了新算法工程版的性能，其单张红外成像的平均处理时间相对于原算法缩短了87%，提升了静态红外地球敏感器姿态确定效率。

4 结 论

针对静态红外地球敏感器的姿态角度计算问题，新算法采用了间隔向量积与t Location-Scale置信区间均值结合的方法。相比于原算法中霍夫变换与最小二乘法结合的方法，新算法无需给定相对投影距离，减少了计算的不确定性。试验表明，在±10°的视场内，新算法的姿态确定精度为0.03°，与原算法相比提高了50%，同时其计算时间缩短了87%。本研究提升了静态红外地球敏感器的测量精度和效率，对于微纳卫星的姿态确定有一定的应用价值。