多视场星/地球敏感器通过在透镜系统前安装分束装置，将视场偏转为多个方向，通过共用成像系及信号处理系统，实现对恒星、地球的同时成像、检测，降低了设备体积和功耗，已在当前航天领域得到广泛应用.近年来，国内众多科研机构在该领域进行了深入的研究和实践，先后研制了涵盖不同轨道高度、不同视场形式、不同成像波长等多种地球^[1-3]、月球敏感器^[4].“对视场成像进行精确仿真，并进行有针对性的算法设计”已成为各类成像式敏感器设计过程中的共性、关键性问题.文献[5-6]首次系统讨论了大气散射的计算模型，根据大气密度变化来仿真大气散射效果，能够有效仿真敏感器观测结果.文献[7-9]则进一步改进了散射模型并应用到不同星球大气的仿真中.在地球敏感器成像检测方法中，文献[10]利用低分辨可见光成像，通过检验地球中心到边缘的距离搜索最优圆心，方法稳定可靠.当前，随着CCD传感器性能的不断提高，新型星敏感器的分辨率已经达到1 024×1 024以上，以地球同步轨道为例，地球扁率造成的边界偏移可达3像素左右，采用已有圆拟合模型无法克服偏移误差，只有构建有效的椭圆拟合算法才能进一步提升地心矢量的检测精度.文献[11]为进一步提高地心检测精度，提出根据已知轨道参数，估计地球长轴方向，利用长轴信息将椭圆边界点映射到圆周上，一定程度上提高了拟合精度.文献[12-13]在物理仿真的基础上，研究了针对微小卫星的可见光地球敏感器，相对于物理仿真，本文研究通过软件仿真，以进一步研究敏感器对于地球表面、轮廓特征、星月背景干扰的适应性.本文在系统分析敏感器成像检测的关键环节的基础上，借鉴航天员模拟器等^[7]已有仿真方法，着重研究高精度椭球成像模型、地球大气轮廓成像仿真、晨昏线仿真方法，构建了具有通用性的地球敏感器成像仿真平台，并应用到北京控制工程研究所某型敏感器研发中，依托该平台对不同地心检测算法进行比较分析，根据系统在轨运行特点，提出了一种分步地球图像检测算法，一方面克服单纯圆拟合方法带来的地球扁率误差，另一方面利用轨道先验知识克服椭圆拟合过程中的长轴不稳定性.有效缩短地球敏感器的研制周期、显著提升了检测性能及系统可靠性.

1 地球敏感器仿真计算模型

如图 1所示，地球敏感器仿真计算可分为“地球光学建模”、“相机成像建模”、“图像检测算法”3部分.其中相机成像建模是一个通用问题，本文着重讨论与敏感器问题具体相关的地球光学建模及相应的图像检测算法.

1.1 地球敏感器成像仿真

成像式地球敏感器首先检测地球轮廓边缘，然后进行拟合，从而得到地球地心矢量.对于大气轮廓成像的准确仿真，是各型敏感器设计实验的难点和关键.文献[5]系统地讨论了地表反射、大气散射作用下的成像仿真，较好模拟了外太空对地观测形态，文献[7]则在此基础上进一步探讨了通过GPU提高算法实时性，并用于航天员训练仿真.本文在借鉴文献[7]的基础上，根据地球敏感器的检测特点，一方面忽略水面、云层建模等对于精度影响较小的地物目标仿真，另一方面加入了对于地球轮廓的精确描述，构建适用于敏感器算法测试的图像仿真系统，主要考虑3方面因素：

1) 地球形状.地球是一个低扁率的椭球体，不考虑地球扁率的情况下，将导致敏感器精度下降，因此在仿真过程中，必须呈现不同视角下椭球轮廓投影.

2) 地、日、月时空位置及光照模型.地、日、月的空间位置及太阳光照模型是地球轮廓、晨昏线仿真的基础.通过月球及恒星仿真，模拟背景星空对敏感器的干扰.

3) 地表、大气反射模型.由于地球大气的存在，形成了一个对太阳光线反射、散射、吸收的复杂过程，造成地球轮廓边界的模糊过渡，这一干扰直接影响地球轮廓的检测精度.

1.2 地心检测算法研究

地心检测算法主要分为边界点检测、圆/椭圆拟合两部分.在边界点检测过程中，关键在于克服地球大气轮廓造成的光晕干扰.圆拟合、椭圆拟合等一般性问题在图像处理、应用数学等领域已得到深入研究，具体到地球边缘拟合上，具有以下特点：地球扁率低，加之受到大气轮廓的干扰，采用圆拟合鲁棒性高而精度低，采用椭圆拟合精度高但鲁棒性低.如何利用轨道特点及先验知识，从系统的角度提升地心矢量的计算精度是敏感器图像处理算法的关键.

实际应用中，轨道参数误差将对检测算法精度产生影响.本文基于图像仿真平台，模拟不同参数误差下的成像状态，进而比较分析不同拟合算法精度及鲁棒性，为算法的优化选择提供数据支持.

2 地球成像仿真数学模型

如图 2所示，构造地球椭球模型，然后基于光路可逆原理，将像素v逐点映射到椭球模型上，获取像素亮度.映射过程分为3种状态：1)映射到地球表面u，则根据u点的光照状态确定亮度；2)映射到大气边缘，则根据大气边缘光照部分u₁-u₂散射状态确定亮度；3)当v点未能映射到地球上，则为图像背景(月球、恒星)部分.

本系统中，儒略日、恒星时计算采用International Astronomical Union提供的开源代码SOFA.太阳轨道计算采用Paul Schlyter的中等精度计算方法，月球轨道计算采用Astronomical Algorithms所述算法实现，地球纹理及星表可根据需要进行更换.图像仿真的关键就在于从v到u的轮廓投影及u点的亮度估计.

2.1 轮廓投影

对于目标(地球)本体坐标系下的齐次椭圆方程，记

则椭球面方程为, 其中u为变量，A为系数矩阵，令L=A^TA写成二次型形式为u^TLu=R², 不同的系数代表不同的二次曲面.

设敏感器坐标到目标本体坐标系的变换为u_o=T_os(u_s-d_os)，令m_o=T_osd_os，则有

式中：T_os为敏感器-目标本体坐标系变换矩阵，m_os为敏感器坐标原点到本体坐标系原点的平移向量(地心向量).由u_o^TLu_o=R²得到观测坐标系中的曲面方程为

(1)

记像平面上一点为v=(x, y, -f), 过原点及v的向量记为αv，设αv与椭球有交点，代入式(1)有

整理得到关于α的二次方程

令：

要使α有解，解得到成像区域方程为

(2)

令，当v=(x, y, -f)^T满足

(3)

时，v属于成像区域，对应的球坐标为

2.2 亮度估计

影响成像亮度的因素主要有:1)太阳距离及入射角度；2)u处的光线吸收性质.

曲面法向量为

u₀处的切平面方程为n₀^tu-u₀=0，当阳光照到该点时有

式中s为太阳中心向量，当cos β=0时为晨昏线.当无遮挡时，该处光强为Icos β，亮度I的计算采用文献[5-6]中的方法

地球大气轮廓是一个与地球扁率相同，半径增加的椭球，适用以上方程.采用以下步骤计算大气成像：

1) 大气边缘判断.根据式(2)计算像点是否属于大气轮廓且不属于地球表面;

2) 亮度估计.如果属于大气边缘，则根据式(3)中的D值，估计大气厚度(方程两根间的距离).

2.3 成像仿真试验

如图 3所示，成像仿真流程分为背景星空仿真、地、月、日仿真、视场叠加3部分.

背景星空仿真中，利用星表计算个星点坐标是否能投射到像平面上，如果能，则根据星等进行亮度仿真.地、日、月仿真中，根据轮廓投影描述的光路可逆过程，逐像素计算点v的投影坐标，顺序判断是否映射到地球、大气轮廓、月球、太阳，并按映射类型进行亮度仿真.最后进行镜头畸变及多视场叠加.

通过修改矩阵L，可以仿真地、日、月等不同任意形式的二次曲面.采用从像点到地球的逐点映射模型，其优点在于可以以任意精度生成仿真图像，敏感器图像仿真的关键就在于地球轮廓的准确描述.因此v点映射到地球轮廓附近时，对其进行亚像素插值，获取高分辨的图像，然后滑动滤波获取逼真的平滑边界.

仿真效果如图 4所示，从图 4(a)、(b)对比可以看出，相对于google Earth，该算法能够生成更加真实、平滑的大气轮廓.从图 4(c)、(d)对比中，可以看出，算法能有效仿真晨昏线，但云层、气旋等大气变化则相对单调，进一步改进中，可以构建单独的云层贴图与地形纹理贴图配合使用.

该方法的缺点在于需要逐点计算、运算量大，利用Intel i7 2.4 G CPU不做深度优化的情况下，单帧1 024×1 024图像仿真时间为0.4 s，测试前，通常需要根据预设参数，事先生成样本集，不适用于一些实时性高的场合.

3 地心检测

地球敏感器的测量原理及过程详见文献[14-15]，本文重点研究地球图像中的地心检测问题，流程如图 5所示，主要分为“初步圆拟合”和“二次椭圆拟合”两个过程.

初步圆拟合中，快速搜索到若干地球轮廓边界点，进行圆拟合，获取初步圆心半径；二次椭圆拟合中，搜索圆周轮廓边界，剔除野点后进行椭圆拟合，以获得高精度的地心坐标及地球半径.

3.1 地心初定位

网格搜索如图 6(a)所示，在图像上以一定间隔利用Sobel算子做水平/垂直搜索，找到候选边界点.然后对各点邻域特征进行进一步判别.文献[7]对判别规则进行了详细讨论，主要包括：前景点数、背景点数、前景背景灰度均值差、前景背景中心距离，利用这些规则可以剔除大部分的恒星干扰、地球纹理干扰以及晨昏线干扰.

判别出有效边界点后，进行圆拟合，以初步确定圆心、半径.该环节的关键在于保证各种干扰环境下的算法的鲁棒性.因此采用Hough变换的思想进行计算：任意3点拟合一个圆，检测其他点到该圆心的距离，小于误差阈值则投票，投票最多的圆为拟合结果，其优点在于鲁棒性强、野点超过70%时也能有效拟合，因而能有效识别、剔除晨昏线.

该方法的缺点在于运算量大，算法实现过程中，根据轨道先验，仅选取相互距离较大的点进行拟合测试，并采用部分投票等方法，能减少运算量.

在敏感器连续运行过程中，利用前后帧图像的相关性，即可获得圆心大致位置信息.初定位仅用于初始帧的检测，实时性要求不高，因此采用该运算量大、但鲁棒性高的算法是恰当的.

3.2 二次拟合

圆周边界搜索如图 6(b)所示，获取地球轮廓大致位置信息后，在轮廓附近，顺序获取地球边界点.相比网格搜索得到的边界点，顺序边界点具有以下优点：仅在轮廓附近搜索、干扰少；边界点彼此相邻、相关有序.

利用这些信息，在圆周边界点中进一步甄选出最稳定、可靠的点进行椭圆拟合：1)边界强度统计.对边界强度进行直方图统计，剔除异常点; 2)边界位置判别.计算边界点到圆心之间距离，剔除异常点；3)相关性判别.相邻边界点的边界强度及到圆心距离不应产生跳变，否则剔除.

利用椭圆模型进行拟合，可以克服圆拟合带来的扁率误差.但由于地球扁率低，直接椭圆拟合难以准确估计椭圆长轴方向，导致精度严重下降，尤其是在地球可见弧段较短时，无法保证鲁棒性，因此需要引入先验知识，给出长轴方向后，再进行拟合.

4 地心检测试验

针对某型高轨道卫星为例进行试验，对比实验已知长轴方向和未知长轴方向下的椭圆最小二乘拟合，参数设定如下：地球敏感器设定为正交双视场传感器，图像分辨率1 024×1 024.网格搜索设定为36×36，圆周边界搜索设定为360个点.

利用图像仿真生成样本：为测试在卫星偏航角未知，即图像中地球椭圆模型长短轴方向未知的情况下，算法的精度.设定卫星偏航角设定为0.5°、1.0°、2.0°、3.0°产生4组样本，以模拟卫星轨道的误差.

敏感器大部分时间不能得到完整的地球边缘，且不同时刻对应的弧长不同、不同的季节对应的弧段位置不同，考虑不同季节太阳照射位置不同，将弧段位置分为春、夏、秋、冬4个季节，每个季节样本1 000张，合计4 000张；考虑不同轨道高度，设定为3~4×10⁴ km，每间隔2 000 km为一级，每组合计6级共计24 000张样本.

根据样本可见轮廓的长短进行分类，可见轮廓大于2/3的为整圆，1/3~2/3的为半圆，统计圆心坐标误差及半径误差.分别对3种方法进行实验：1)最小二乘圆拟合实验，结果见表 1；2)未知轨道先验知识的椭圆拟合：假设椭圆长轴方向未知、扁率未知，直接进行五参数椭圆拟合测试，结果见表 2；3)已知轨道先验知识的椭圆拟合：假设椭圆长轴方向已知、偏航角为0、扁率已知.进行三参数椭圆拟合，结果见表 3.

通过以上算法比较看出，圆拟合方法对于整圆圆心估计精度尚可，但对于半径估计误差较大，尤其是偏航角增大时.与之相比，对于算法2当长短轴未知时，直接进行椭圆拟合精度较差，尤其是半圆情况下.这与已有文献[10-11]中提出的椭圆鲁棒性差的观点是一致的.分析其原因，由于椭圆拟合未能正确找到长轴方向，错误椭圆模型非但不能带来拟合精度的提高，相反将导致误差放大，这也是已有文献中，不建议直接采用椭圆拟合的原因.对于算法3，利用了轨道先验知识，并模拟了实际的误差状态，在小角度偏航角的情况下，利用先验知识确定轨道长轴及扁率.在长轴方向误差较小时，椭圆模型依然能提高拟合精度.因此在实际应用中，利用先验知识可以有效保证椭圆算法鲁棒性.

实验过程中，根据轮廓长度(有效边界点数)与圆周长度的比例，统计不同比例下的地球中心拟合误差及地球半径拟合误差，如图 7所示.从图 7中可以看出，当可见轮廓超过40%时，现有算法能较准确计算中心，当可见轮廓超过80%，误差基本控制在0.1像素以下.当可见轮廓少于30%时误差迅速增大.半径拟合误差在可见轮廓大于35%时基本保持稳定.可见轮廓小于30%时，误差增大，但相对中心误差增速较缓慢，可见在可见轮廓较低时，中心误差是算法的主要问题.

5 结论

1) 针对地球敏感器检测算法的特点，研究地球椭球成像、大气轮廓、晨昏线及背景星空成像，构建了高精度地球成像仿真平台.

2) 在地球成像仿真平台上，验证测试了针对不同应用场合的敏感器算法，并根据卫星在轨运行特点，提出了分步地球图像检测算法.该方法一方面克服单纯圆拟合方法带来的地球扁率误差，另一方面利用轨道先验知识克服椭圆拟合过程中的长轴不稳定性.实验证明，在卫星存在小角度偏航角误差的情况下，采用椭圆拟合依然能够有效提高了敏感器检测精度及系统鲁棒性.通过构建该图像仿真及算法测试平台，将有效缩短各型敏感器的研发周期、保障产品质量.