CN105928526B - 一种基于可见光地球敏感器确定卫星姿态的方法 - Google Patents

Abstract

一种基于可见光地球敏感器确定卫星姿态的方法，利用可见光地球敏感器确定卫星对地的两轴姿态。首先通过对CMOS成像器件所成的地球图像处理，提取到地球图像的有效边缘，然后应用最小二乘法拟合提取到的有效边缘的中心。最后在卫星对地姿态角为小角度的假设下，近似计算卫星的对地滚动角和俯仰角，从而最终确定卫星对地的两轴姿态角。相对于传统的利用红外地球敏感器确定卫星对地的滚动和俯仰两轴姿态的方法，本发明方法提高了卫星对地两轴姿态角测量精度，并且该方法适用于小型化的小卫星姿态测量，顺应当前小卫星的发展潮流。

Description

一种基于可见光地球敏感器确定卫星姿态的方法

技术领域

本发明属于卫星控制领域，涉及一种确定卫星姿态的方法。

背景技术

目前，在以地球作为参考矢量确定卫星姿态角的敏感器中，红外地球敏感器得到了广泛的应用，并且日趋成熟，因此利用红外地球敏感器确定卫星对地的两轴姿态角成为了常见方法。

现有的利用红外地球敏感器确定卫星对地两轴姿态角的方法是：一方面，利用红外地球敏感器扫描机构进出地球二氧化碳红外辐射带以及红外扫描机构到达基准点的时间，计算出红外地球敏感器扫描地球的弦宽和地心矢量E在红外地球敏感器的测量坐标系x_SEoy_SE平面中的投影与x_SE轴的夹角，由上述计算所得的夹角以及球面三角的知识可以计算出地心矢量E在红外地球敏感器测量坐标系中的坐标表示，又由红外地球敏感器的安装矩阵可得地心矢量E在星体坐标系中的坐标表示；另一方面，地心矢量E在卫星轨道坐标系的坐标已知，由未知卫星姿态角计算的旋转矩阵可计算得地心矢量E在星体坐标系的坐标表示，由上述两种途径转换成的地心矢量在星体坐标系下的不同表达式，即可求出卫星对地的滚动和俯仰角。

由上述红外地球敏感器确定卫星对地三轴姿态的原理可知，卫星姿态角的测量精度受限于时间测量精度和扫描机构扫描角速度的测量精度，以及卫星轨道参数测量精度(视不同的安装矩阵)。由于上述测量物理量均具有不确定性，因此造成卫星姿态的测量精度波动较大。

另外，当前卫星正朝着小型化、智能化、高精度、高稳定性、低功耗和长寿命方向发展，而主要的红外地球敏感器有圆锥扫描和摆动扫描等，这类红外地球敏感器由于存在扫描机构，其质量大、功耗多，也不适合当前卫星的发展潮流。

发明内容

本发明解决的技术问题是：克服现有技术的不足，提供了一种基于可见光地球敏感器确定卫星姿态的方法，解决了现有红外地球敏感器确定卫星姿态中时间测量精度、扫描机构角速度测量精度、卫星轨道参数测量精度导致卫星姿态角测量精度下降的问题，并且该方法适用小型化的小卫星姿态测量，顺应当前小卫星的发展潮流。

本发明的技术解决方案是：一种基于可见光地球敏感器确定卫星姿态的方法，包括如下步骤：

(1)利用可见光地球敏感器中的CMOS成像器件对地球进行成像，获取输入图像；

(2)对所述的输入图像进行预处理，去除噪声，平滑图像；

(3)对步骤(2)获取的图像进行边缘检测，获取有效边缘；

(4)对所述有效边缘进行拟合，获取有效边缘中心和半径；

(5)计算获得卫星对地的滚动角φ和俯仰角θ，其中

(C_x,C_y)为有效边缘中心坐标，R为有效边缘的半径，f为CMOS成像器件的焦距。

所述步骤(1)中的输入图像应满足：CMOS成像器件所成地球图像地平轮廓有效弧长占到总弧长的50％。所述步骤(2)中平滑图像的方法为3×3卷积模板的中值滤波。所述步骤(3)中边缘检测的方法为：首先采用Sobel算子检测所有的边缘，然后采用Hough变换检测出有效边缘。所述步骤(4)中对有效边缘进行拟合的方法为最小二乘法。

本发明与现有技术相比的优点在于：

(1)本发明方法通过可见光地球敏感器确定卫星对地的两轴姿态角，相比于传统的红外地球敏感器确定卫星的对地姿态角，该方法消除了时间测量精度、扫描机构角速度测量精度、轨道参数测量精度对红外地球敏感器测量卫星姿态角的影响，从而提高卫星姿态测量的精度，而且该方法算法简单，实现容易，对卫星平台的控制要求不高；

(2)本发明方法的重点在于图像有效边缘的提取，边缘提取的好坏直接决定方法的精度。本发明采取的提取边缘方法不仅能精确提取有效边缘，满足计算量的要求，而且算法的鲁棒性很好；

(3)本发明方法可集成到控制器中，与其他部组件共同形成轻质量、精度满足要求的可见光敏感器，为确定卫星对地两轴姿态角提供另外一种解决方案。

附图说明

图1为本发明方法的流程框图；

图2为本发明CMOS相机成像示意图；

图3为本发明卫星姿态坐标系示意图；

具体实施方式

随着微小型化CMOS成像器件的快速发展，不仅为地球可见光敏感器的真正使用提供了极大可能，而且顺应了卫星的发展趋势。本发明就基于可见光地球敏感器如何确定卫星对地两轴姿态，提出一套完整的方法。

如图1流程所示，本发明方法主要包括以下几个步骤：第一，CMOS成像器件对地球成像后输入图像；第二，图像预处理；第三，提取图像有效边缘；第四，应用最小二乘拟合图像有效边缘中心；第五，在卫星对地姿态角为小角度假设前提下，根据有效边缘中心近似计算卫星滚动和俯仰姿态角。以下分别详细介绍上述五个步骤。

第一，可见光地球敏感器中的CMOS成像器件对地球成像后输入图像。

为了保证最终确定的卫星俯仰角和滚动角精度，要求CMOS成像器件所成地球图像地平轮廓有足够长的弧度。一般要求成像平面有效弧长占到总弧长的50％。

CMOS成像器件所拍摄的图像满足弧长要求后输入图像，进行下一步的处理。

第二，图像预处理。

对图像进行预处理的目的在于去除噪声，平滑图像。图像预处理的方法有很多，如中值滤波，均值滤波等，本发明中采用中值滤波。

中值滤波采用3×3卷积模板，具体做法如下：将中心像素及其3×3邻域内8个像素的灰度值按从小到大排序，取第5个元素灰度值为中心像素的灰度值。中值滤波不仅能去除噪声点，为以后边缘提取减少计算量，还能保护图像的边缘，获得较满意的图像预处理效果。

第三，检测有效边缘。

这一步是本发明方法为解决上述所提出问题的关键步骤，也是核心部分。检测有效边缘方法的好坏直接决定了本发明方法所需计算量的大小和所确定姿态角的精度。

目前检测边缘的算法有很多种，常见的有：微分算子、拉普拉斯高斯算子、Canny算子、拟合法等。本发明采用Sobel微分算子和Hough变换检测有效边缘，Sobel算子是检测所有的边缘，包括有效边缘和无效边缘，而Hough变换则是剔除Sobel算子所检测到的无效边缘，保留其有效边缘。下面分别介绍这两种方法。

Sobel算子首先计算某一像素点处的灰度梯度幅值，然后将此值与所设的阀值进行比较，如果该点处灰度梯度幅值大于阀值，则保留该点，否则剔除该点。该方法的中心思想是图像边缘像素点邻域内灰度值变化比较剧烈，应用此方法确定出邻域内灰度值变化比较大的像素点，从而检测出图像边缘。Sobel算子中所设阀值的大小不变，阀值大小需要根据实际图像确定，一般阀值的大小选取为不大于整个图像灰度值平均值的最大正整数。

下面介绍Sobel算子计算像素点处灰度梯度幅值的方法。

任一点处梯度幅值计算如下：

其中s_x，s_y分别表示在该点处x、y方向导数，而对于数字图像，可以用一阶差分代替一阶微分。通常计算某一像素点处灰度梯度幅值有3×3模板、5×5模板等，本发明采用3×3模板，故s_x，s_y可写为：

s_x＝(a₂+2a₃+a₄)-(a₀+a₇+a₆) (2)

s_y＝(a₀+2a₁+a₂)-(a₆+2a₅+a₄) (3)

其中a₀…a₇代表该像素点处3×3邻域内元素的灰度值，a₀…a₇各像素点灰度值在3×3邻域内位置如下：

a<sub>0</sub>	a<sub>1</sub>	a<sub>2</sub>
			a<sub>7</sub>	[i,j]	a<sub>3</sub>
a<sub>6</sub>	a<sub>5</sub>	a<sub>4</sub>

注：[i,j]代表所求灰度梯度幅值像素点坐标。

如前所述，Sobel算子检测到的边缘包括有效边缘和无效边缘，其中有效边缘是地球的轮廓边缘，无效边缘包括地球晨昏线边缘、图像背景明亮星体边缘、月球边缘等。为了提高所确定姿态角的精度，需要对Sobel算子检测到的边缘进一步检测，剔除无效边缘，保留有效边缘。下面介绍检测有效边缘方法之Hough变换。

Hough变换有已知图像圆半径Hough变换和未知图像圆半径Hough变换两种。对于本问题，当然属于后者，因为地球成像后所形成图像圆的半径是无法提前精确获取的，一种快速的未知图像圆半径的Hough变换方法如下：

首先，对于圆的检测，利用梯度信息，Ballard用下式检测圆

x₀＝x±rcos(h) (4)

y₀＝y±rsin(h) (5)

式中，(x₀,y₀)为圆心坐标，h为圆边缘点(x,y)处梯度方向角，r为图像圆半径。

基于上述公式，检测步骤如下：

1.首先根据先验知识或者已知数据估计参数半径r的变化范围r∈[r_min,r_max]以及步长Δr大小。

对于本问题，地球成像后成像面上近似圆形的半径r范围确定如下：

首先，相机成像示意图如图2所示，图中，R为相机可分辨到的大气层高度，变化范围R∈[6371,6371+120]；H为卫星高度，已知值；为半张角，大小R_E为地球半径，则由此可确定半张角的范围；f为相机焦距，已知值；r为地球成像后成像面上近似圆形的半径，由可确定r的范围r∈[r_min,r_max]；步长Δr由敏感器控制器的计算能力以及精度要求等决定，如果控制器的计算能力强，那步长就可以选的小一点，不仅满足姿态精度要求，同时满足姿态实时性的要求。

从半径r₀＝r_min开始，每一次循环半径增加Δr，直至r＝r_max，对于任一循环，半径r_i＝r_min+Δr×i,(i＝0,1…m),其中m为循环的总次数，表达式将每一边缘坐标点(x_j,y_j)(j＝1,…n)及其梯度方向角h_j依次代入上述(4)和(5)Ballard公式中，计算出圆心坐标(x_0ij,y_0ij)，其中j表示有效和无效边缘点总个数，然后将计算所得的圆心坐标四舍五入取整，上述各个边缘点坐标是在以图像左下角为坐标原点，向右为x轴，向上为y轴的坐标系下表示。为防止计算得到的整数圆心坐标为负数，将原有坐标系分别沿x轴，y轴负方向平移1024个像素点，平移后得到新的正的整数圆心坐标值，以新坐标值为二维数组的坐标值，二维数组的值统计累加该二维坐标值出现的次数，遍历所有的边缘点，最后由累加次数峰值确定此次循环r＝r_i时的圆参数c_i＝{x_0i,y_0i,r_i}，每一次循环确定一组参数c_i＝{x_0i,y_0i,r_i}和累加次数K_i。通过此次循环累加次数K_i和上一次循环累加次数K_i-1的比较，记录比较高的累加次数K_max以及相应的圆参数c_max＝{X₀,Y₀,R}，当r_i＝r_max时，最终确定的圆参数c_max＝{X₀,Y₀,R}为有效边缘像素级圆心坐标{X₀,Y₀}以及半径R。

(4)和(5)Ballard公式均取正号即可，灰度梯度方向角h的正弦和余弦计算如下：

其中s_x,s_y由(2)和(3)计算可得。

2.未知半径Hough变换得到像素级圆心坐标{X₀,Y₀}以及半径R，如果某边缘点为有效边缘点，即圆周上的边界点，则该点与像素级圆心坐标{X₀,Y₀}之间距离满足：

不符合(8)关系式的干扰边缘点将被删除，从而获得有效边缘。

第四，应用最小二乘拟合有效边缘中心；

得到有效边缘之后，接下来应用最小二乘拟合有效边缘中心。

对给定的一组数据点(x_k,y_k)(k＝1,…l)，l为数据点个数，至少为3，设与该组数据点的距离的平方和为最小的圆的方程为：

(x-x₀)²+(y-y₀)²＝r² (9)

式中(x,y)表示上述圆周上的坐标值，(x₀,y₀)表示上述圆的圆心坐标，r表示上述圆的半径。

又设a₁＝∑x_k，b₁＝∑y_k，a₂＝∑(x_k)²，b₂＝∑(y_k)²，a₃＝∑(x_k)³，b₃＝∑(y_k)³，c₁₁＝∑x_k*y_k，c₁₂＝∑x_k*y_k ²，c₂₁＝∑(x_k)²*y_k。

设f(x₀,y₀,R)＝∑((x_k-x₀)²+(y_k-y₀)²-r²)²，分别令 即

-4∑((x_k-x₀)²+(y_k-y₀)²)(x_k-x₀)＝0

-4∑((x_k-x₀)²+(y_k-y₀)²-R²)(y_k-y₀)＝0

-4∑((x_k-x₀)²+(y_k-y₀)²)r＝0

将a₁，b₁，a₂，b₂，a₃，b₃，c₁₁，c₁₂，c₂₁代入上面三个方程整理得：

a₁((x₀)²+(y₀)²)-2a₂x₀-2c₁₁y₀+a₃+c₁₂-r²a₁＝0

b₁((x₀)²+(y₀)²)-2c₁₁x₀-2b₂y₀+c₂₁+b₃-r²b₁＝0

n((x₀)²+(y₀)²)-2a₁x₀-2a₁x₀+a₂+b₂-lr²＝0

解得圆心坐标及半径为：

(10)、(11)、(12)计算所得结果是坐标数值结果，没有实际物理意义，第五步中需要有物理意义的结果，设单位像素长度为d，给上述(10)、(11)、(12)计算结果都乘以单位像素长度d，所得结果为有长度物理意义的圆心坐标和半径：

C_x＝x₀×d (13)

C_y＝y₀×d (14)

R＝r×d (15)

第五，计算卫星对地的滚动角和俯仰角。

在计算卫星姿态角之前，首先建立卫星姿态坐标系O-XYZ(如附图3)和敏感器坐标器O-X_TY_TZ_T。卫星姿态坐标系O-XYZ定义如下：原点O为卫星质心，OX轴指向卫星飞行方向，OZ轴指向地心，OY轴由右手定则确定。敏感器坐标器O_T-X_TY_TZ_T各个坐标轴均与卫星姿态坐标系坐标轴平行，唯一不同原点O_T位于敏感器的质心。

根据上述坐标系的定义，姿态角有如下说明：卫星绕X轴旋转为滚动角，绕Y轴旋转为俯仰角，绕Z轴旋转为偏航角，均以绕X轴，Y轴，Z轴正向旋转为正的姿态角。

在此基础上，确定卫星的滚动角φ和俯仰角θ：附图3中，为半张角，由附图2可得式中，R由(15)求得，f已知，则半张角有：

设地心在成像面中的坐标为(M_x,M_y)，由图像几何中心即最小二乘所得的有效边缘中心坐标为(C_x,C_y)，其中地心在成像面中的坐标(M_x,M_y)有：

M_x＝f*tanθ (17)

M_y＝f*tanφ (18)

推导有效边缘中心坐标(C_x,C_y)的计算表达式，为说明问题简单起见，假设卫星滚动角φ＝0，俯仰角θ≠0，则成像面上在X轴方向距离有效边缘中心最远点和最近点的坐标分别为则有效边缘中心X坐标为：C_x＝X_max-X_min，即：

同理可得：

在滚动角φ和俯仰角θ为小角度下，由(17)和(19)、(18)和(20)可得：

故由(17)和(21)、(18)和(22)以及半张角计算公式(16)可得卫星滚动角φ和俯仰角θ有：

首先卫星携带的可见光地球敏感器的CMOS成像器件对空间中地球成像，然后输入图像，其次检测所得图像边缘以及剔除干扰边缘，从而得到地球成像图像的最终有效边缘，再者应用最小二乘拟合出有效边缘的中心位置，接着，根据图像中心位置坐标及卫星姿态角为小角度假设下近似计算卫星滚动和俯仰两轴姿态角。至此，卫星对地的两轴姿态角均已确定，为计算卫星的对地偏航角提供数据支撑。

本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。

Claims (1)

1.一种基于可见光地球敏感器确定卫星姿态的方法，其特征在于包括如下步骤：

(1)利用可见光地球敏感器中的CMOS成像器件对地球进行成像，获取输入图像；输入图像应满足：CMOS成像器件所成地球图像地平轮廓有效弧长占到总弧长的50％；

(2)对所述的输入图像进行预处理，去除噪声，平滑图像；平滑图像的方法为3×3卷积模板的中值滤波，具体做法如下：将中心像素及其3×3邻域内8个像素的灰度值按从小到大排序，取第5个元素灰度值为中心像素的灰度值；

(3)对步骤(2)获取的图像进行边缘检测，获取有效边缘；边缘检测的方法为：首先采用Sobel算子检测所有的边缘，采用Sobel算子计算像素点处灰度梯度幅值的方法为其中s_x，s_y分别表示在像素点处x、y方向导数，s_x＝(a₂+2a₃+a₄)-(a₀+a₇+a₆)，s_y＝(a₀+2a₁+a₂)-(a₆+2a₅+a₄)，其中a₀…a₇代表该像素点处3×3邻域内元素的灰度值，a₀…a₇各像素点灰度值在3×3邻域内位置如下：

a₀ a₁ a₂ a₇ [i,j] a₃ a₆ a₅ a₄

[i,j]代表所求灰度梯度幅值像素点坐标；

然后采用Hough变换检测出有效边缘，Hough变换检测采用未知图像圆半径Hough变换，具体如下：

A、首先根据先验知识或者已知数据估计参数半径r的变化范围r∈[r_min,r_max]以及步长Δr大小；

H为卫星高度，为半张角，R_E为地球半径，f为相机焦距，r为地球成像后成像面上近似圆形的半径，由此确定r的范围r∈[r_min,r_max]；

步长Δr由敏感器控制器的计算能力以及精度要求决定；

从半径r₀＝r_min开始，每一次循环半径增加Δr，直至r＝r_max，对于任一循环，半径r_i＝r_min+Δr×i,(i＝0,1…m),其中m为循环的总次数，表达式将每一边缘坐标点(x_j,y_j)(j＝1,…n)及其梯度方向角h_j依次代入x₀＝x±rcos(h)和y₀＝y±rsin(h)，(x₀,y₀)为圆心坐标，h为圆边缘点(x,y)处梯度方向角，r为图像圆半径，计算出圆心坐标(x_0ij,y_0ij)，其中j表示有效和无效边缘点总个数，然后将计算所得的圆心坐标四舍五入取整，上述各个边缘点坐标是在以图像左下角为坐标原点，向右为x轴，向上为y轴的坐标系下表示；为防止计算得到的整数圆心坐标为负数，将原有坐标系分别沿x轴，y轴负方向平移1024个像素点，平移后得到新的正的整数圆心坐标值，以新坐标值为二维数组的坐标值，二维数组的值统计累加该二维坐标值出现的次数，遍历所有的边缘点，最后由累加次数峰值确定此次循环r＝r_i时的圆参数c_i＝{x_0i,y_0i,r_i}，每一次循环确定一组参数c_i＝{x_0i,y_0i,r_i}和累加次数K_i；通过此次循环累加次数K_i和上一次循环累加次数K_i-1的比较，记录比较高的累加次数K_max以及相应的圆参数c_max＝{X₀,Y₀,R}，当r_i＝r_max时，最终确定的圆参数c_max＝{X₀,Y₀,R}为有效边缘像素级圆心坐标{X₀,Y₀}以及半径R；

B、未知半径Hough变换得到像素级圆心坐标{X₀,Y₀}以及半径R，如果某边缘点为有效边缘点，即圆周上的边界点，则该点与像素级圆心坐标{X₀,Y₀}之间距离满足不符合该关系式的干扰边缘点将被删除，从而获得有效边缘；

(4)对所述有效边缘进行拟合，获取有效边缘中心和半径；对有效边缘进行拟合的方法为最小二乘法，具体圆心坐标为：C_x＝x₀×d，C_y＝y₀×d，半径为：R＝r×d，d为单位像素长度，

(x_k,y_k)，k＝1,…l为给定的一组数据点，l为数据点个数，至少为3，a₁＝∑x_k，b₁＝∑y_k，a₂＝∑(x_k)²，b₂＝∑(y_k)²，a₃＝∑(x_k)³，b₃＝∑(y_k)³，c₁₁＝∑x_k*y_k，c₁₂＝∑x_k*y_k ²，c₂₁＝∑(x_k)²*y_k；

(5)计算获得卫星对地的滚动角φ和俯仰角θ，其中

(C_x,C_y)为有效边缘中心坐标，R为有效边缘的半径，f为CMOS成像器件的焦距。