CN103617613A - 一种微小卫星非合作目标图像处理方法 - Google Patents

Abstract

本发明公开了一种微小卫星非合作目标图像处理方法，通过计算目标星在一帧初始图像中所占的像素个数N，判断像素个数N是否大于一预设阈值，若否，则基于反馈直线分离的方法识别出所述帧的初始图像中目标星的特征；若是，则基于反馈直线聚合的方法识别出所述帧的初始图像中目标星的特征，能够直接提取目标星本身的特征、无需对目标星进行加装改造、对目标星本身无特别要求；相对于合作目标图像处理方法，执行非合作目标图像处理方法的相关硬件结构具有低功耗、小重量；针对目标星相对距离不同自适应的采用不同的特征提取算法、实现大范围特征提取，只需要目标星有适合的光照就可以提取所需的目标星的特征。

Description

一种微小卫星非合作目标图像处理方法

技术领域

本发明属于航天器测量与图像处理技术领域，具体地说,涉及一种微小卫星非合作目标图像处理方法。

背景技术

在空间交会对接过程中，需要精确测量追踪航天器与目标星之间的相对位置、相对姿态角、相对速度和相对姿态角速度等相对运动信息。为获得准确的被测目标姿态信息，现有的合作目标算法需要对目标进行一定的改装或者加装处理，增加了目标星的功耗、重量，因此，目前亟需一种能够有效的克服这些缺点的非合作目标图像处理方法。

发明内容

针对目标星相对距离不同自适应的采用不同的特征提取算法、实现大范围特征提取，本发明提出一种微小卫星非合作目标图像处理方法。

本发明解决其技术问题所采用的技术方案是:通过计算目标星在一帧初始图像中所占的像素个数N，判断像素个数N是否大于一预设阈值，若否，则基于反馈直线分离的方法识别出所述帧的初始图像中目标星的特征；若是，则基于反馈直线聚合的方法识别出所述帧的初始图像中目标星的特征，能够直接提取目标星本身的特征、无需对目标星进行加装改造、对目标星本身无特别要求；相对于合作目标图像处理方法，执行非合作目标图像处理方法的相关硬件结构具有低功耗、小重量；针对目标星相对距离不同自适应的采用不同的特征提取算法、实现大范围特征提取，只需要目标星有适合的光照就可以提取所需的目标星的特征。

本发明微小卫星非合作目标图像处理方法，其特点是包括以下步骤:

步骤1.计算目标星在一帧初始图像中所占的像素个数N，判断像素个数N是否大于一预设阈值，若否，则基于反馈直线分离的方法识别出初始图像中目标星的特征；若是，则基于反馈直线聚合的方法识别出初始图像中目标星的特征；

步骤2.计算目标星在相面中所占的像素个数N，采用大律法计算目标星在相面中所占的像素个数N；

步骤3.基于反馈直线分离的方法识别出初始图像中目标星的特征:对初始图像的每个像素点的像素值f(x,y)进行中值滤波得到滤波后的图像的每个像素点的像素值f1(x,y)；对滤波后的图像进行双线性插值得到插值后的图像；对插值后的图像的每个像素点进行sobel算子滤波、非极大抑制、二值化得到二值化后的边缘图像的每个像素点的像素值G(x,y)；对二值化后的边缘图像的每个像素点的像素值G(x,y)进行Hough直线检测以获取目标星的每个边缘i的边缘直线（ρi，θi）作为目标星的直线特征，i为正整数；

步骤4.对初始图像的每个像素点的像素值f(x,y)进行中值滤波，得到滤波后的图像的每个像素点的像素值f1(x,y):获取待进行中值滤波的像素点的第一3×3阶邻域窗口；对所述3×3阶邻域窗口中的9个像素点分别分配一个初值为0的计数器，并将所述第一3×3阶邻域窗口中的每个待比较像素点的像素值分别与其它8个像素点的像素值进行比较，每次若大于其它8个像素点中的一个像素点的像素值，则每次该待比较像素点的计数值+1，否则，每次该待比较像素点的计数值不变；对所述第一3×3阶邻域窗口中的每个像素点的计数器的计数值进行排序，获取排序为中间的计数值作为待进行中值滤波的像素点的滤波后的像素值f1(x,y)；

步骤5.对滤波后的图像进行双线性插值得到插值后的图像:将滤波后的图像分割为若干个4×4阶邻域窗口；

根据公式 $\begin{array}{c}{f}_2(x+u,y+v)=(1-u)(1-v)f_1(x,y)+(1-u){\mathrm{vf}}_1(x,y+1)+\\ u(1-v)f_1(x+1,y)+{\mathrm{uvf}}_1(x+1,y+1)\end{array}$ 获取每个4×4阶邻域窗口内插值后的点的像素值f₂(x,y)，其中u，v是[0，1]区间的浮点数，（u，v）取值分别为(0，0.5)、(0.5，0)、(0.5，0.5)、(0.5，1)、(1，0.5)，将得到的插值后的图像中的每个滤波后的像素点和插值后的点的像素值统一计为f₂’(x,y)；

步骤6.对插值后的图像的每个像素点进行sobel算子滤波:在第二3×3阶邻域窗口内根据公式g(x,y)=|g_x(x,y)+g_y(x,y)|获取每个待进行sobel算子滤波的像素点的梯度幅值g(x,y)，其中，

g_x(x,y)=-f′₂(x-1,y-1)-2·f′₂(x-1,y)-f′₂(x-1,y+1)+

f′₂(x+1,y-1)+2·f′₂(x+1,y)+f′₂(x+1,y+1)，

g_y(x,y)=f′₂(x-1,y-1)+2·f′₂(x,y-1)+f′₂(x+1,y-1)-

f′₂(x-1,y+1)-2·f′₂(x,y+1)-f′₂(x+1,y+1)，

g_x(x,y)，g_y(x,y)分别为每个待进行sobel算子滤波的像素点的X方向和Y方向的梯度幅值；

步骤7.对插值后的图像的每个像素点进行非极大抑制:根据sobel算子滤波的结果判断插值后的图像的每个像素点C的梯度幅值是否大于梯度方向上相邻的两点dTmp1和dTmp2的梯度幅值，若是，则像素点C的梯度幅值不变；若否，将该像素点C的梯度幅值置0；

步骤8.对插值后的图像的每个像素点进行二值化得到二值化后的边缘图像的所有像素点的像素值G(x,y):判断非极大抑制后的每一个像素点是否大于预设的梯度幅值阈值，若是，则将该像素点的梯度幅值置为1，若否，则将该像素点的梯度幅值置为0；提取每个梯度幅值为1像素点的像素值作为二值化后的边缘图像的像素点的像素值G(x,y)；

步骤9.对二值化后的边缘图像的每个像素点的像素值G(x,y)进行Hough直线检测以获取目标星的每个边缘i的边缘直线（ρi，θi）作为目标星的直线特征:设Hough平面累加器（ρi，θi），其中θi从0到180°，θi的步长为1°，ρi为

其中，w为图像宽度、h为图像高度，ρi的步长为一个像素；对二值化后的边缘图像的每个像素点进行Hough变换：依次将θi从0到180°带入直线极坐标方程ρi=xcosθi+ysinθi中计算，每次计算得（ρi，θi）在累加器的相应单元的计数值+1；Hough变换处理完成后，设定阈值T_hough，累加器中大于T_hough的单元对应的（ρi，θi）作为检测出的边缘直线的直线特征，其中，T_hough=0.5·(Amax-Amin)，Amin、Amax分别为累加器的所有单元的计算值的最小值和最大值；

步骤10.对二值化后的边缘图像的每个像素点的像素值G(x,y)进行Hough直线检测以获取目标星的每个边缘i的边缘直线（ρi，θi）作为目标星的直线特征，从第二帧的初始图像开始，减少Hough空间的搜索范围；

步骤11.对二值化后的边缘图像的每个像素点的像素值G(x,y)进行Hough直线检测以获取目标星的每个边缘i的边缘直线（ρi，θi）作为目标星的直线特征，并对G(x,y)进行harris角点提取以获取目标星的角点特征；

步骤12.对G(x,y)进行harris角点提取以获取目标星的角点特征:对二值化后的边缘图像的每个像素点的像素值G(x,y)进行x、y和xy方向离散求导得Ix、Iy和Ixy:

Ix=-2·G(x-2,y)-G(x-1,y)+G(x+1,y)+2·G(x+2,y)，

Iy=-2·G(x,y-2)-G(x,y-1)+G(x,y+1)+2·G(x,y+2)，

Ixy=Ix·Iy；

对每个Ix和Iy做平方得到Ix²、Iy²；采用5阶高斯滤波器对Ix²、Iy²和Ixy分别进行5阶高斯滤波，得到滤波后的结果:Ix²*、Iy²*、Ixy*，其中，5阶高斯滤波器如下:h=[15851；

       52134215；

       83456348；

       52134215；

       15851]；

对二值化后的边缘图像的每个像素点(x,y)求对应的R=（AB-CD）²-k（A+B）²，其中，A=Ix²*(i,j)、B=Iy²*(i,j)、C=D=Ixy*(i,j)，k为响应系数，取0.04；选取角点阈值，如果该像素点对应的R大于所述角点阈值，并且对应的R在该像素点的3×3阶邻域窗口内最大，则将该像素点作为目标星的角点特征作为角点，其中角点阈值为该像素点的3×3阶邻域窗口内最大R值的0.01倍；

步骤13.基于反馈直线聚合的方法识别出所述帧的初始图像中目标星的特征:对所述帧的初始图像的每个像素点的像素值f(x,y)进行中值滤波得到滤波后的图像的每个像素点的像素值f1(x,y)；对滤波后的图像的每个像素点进行sobel算子滤波、非极大抑制、二值化得到二值化后的边缘图像的每个像素点的像素值G(x,y)；对二值化后的边缘图像的每个像素点的像素值G(x,y)进行Hough直线检测以获取目标星的每个边缘i的边缘直线（ρi，θi），i为正整数；对目标星的表示同一边缘i的多条边缘直线（ρi，θi）进行聚合以获取一条聚合后的边缘直线（ρi’，θi’）作为目标星的直线特征；

步骤14.对所述帧的初始图像的每个像素点的像素值f(x,y)进行中值滤波得到滤波后的图像的每个像素点的像素值f1(x,y):获取待进行中值滤波的像素点的第一3×3阶邻域窗口；对所述3×3阶邻域窗口中的9个像素点分别分配一个初值为0的计数器，并将所述第一3×3阶邻域窗口中的每个待比较像素点的像素值分别与其它8个像素点的像素值进行比较，每次若大于其它8个像素点中的一个像素点的像素值，则每次该待比较像素点的计数值+1，否则，每次该待比较像素点的计数值不变；对所述第一3×3阶邻域窗口中的每个像素点的计数器的计数值进行排序，获取排序为中间的计数值作为待进行中值滤波的像素点的滤波后的像素值f1(x,y)；

步骤15.对滤波后的图像的每个像素点进行sobel算子滤波:将每个待进行sobel算子滤波的像素点的第二3×3阶邻域窗口；在第二3×3阶邻域窗口内根据公式g(x,y)=|g_x(x,y)+g_y(x,y)|获取每个待进行sobel算子滤波的像素点的梯度幅值g(x,y)，其中，

g_x(x,y)=-f₁(x-1,y-1)-2·f₁(x-1,y)-f₁(x-1,y+1)+

f₁(x+1,y-1)+2·f₁(x+1,y)+f₁(x+1,y+1)，

g_y(x,y)=f₁(x-1,y-1)+2·f₁(x,y-1)+f₁(x+1,y-1)-

f₁(x-1,y+1)-2·f₁(x,y+1)-f₁(x+1,y+1)，

g_x(x,y)，g_y(x,y)分别为每个待进行sobel算子滤波的像素点的X方向和Y方向的梯度幅值；

步骤16.对滤波后的图像的每个像素点进行非极大抑制:根据sobel算子滤波的结果判断插值后的图像的每个像素点C的梯度幅值是否大于梯度方向上相邻的两点dTmp1和dTmp2的梯度幅值，若是，则像素点C的梯度幅值不变；若否，将该像素点C的梯度幅值置0；

步骤17.对滤波后的图像的每个像素点进行二值化得到二值化后的边缘图像的每个像素点的像素值G(x,y):判断非极大抑制后的每一个像素点是否大于一预设的梯度幅值阈值，若是，则将该像素点的梯度幅值置为1，若否，则将该像素点的梯度幅值置为0；提取每个梯度幅值为1像素点的像素值作为二值化后的边缘图像的像素点的像素值G(x,y)；

步骤18.对二值化后的边缘图像的每个像素点的像素值G(x,y)进行Hough直线检测以获取目标星的每个边缘i的边缘直线（ρi，θi）:设定Hough平面累加器（ρi，θi），其中θi从0到180°，θi的步长为1°，ρi为

其中，w为图像宽度、h为图像高度，ρi的步长为一个像素；对二值化后的边缘图像的每个像素点进行Hough变换：依次将θi从0到180°带入直线极坐标方程ρi=xcosθi+ysinθi中计算，每次计算的（ρi，θi）在累加器的相应单元的计数值+1；Hough变换处理完成后，设定阈值T_hough，累加器中大于T_hough的单元对应的（ρi，θi）作为检测出的边缘直线，其中，T_hough=0.5·(Amax-Amin)，Amin、Amax分别为累加器的所有单元的计算值的最小值和最大值；

步骤19.对目标星的表示同一边缘i的多条边缘直线（ρi，θi）进行聚合以获取一条聚合后的边缘直线（ρi’，θi’）作为目标星的直线特征:

（1）令从所述Hough直线检测出的所有边缘直线（ρi，θi）作为直线集L；

（2）从直线集L中取出第一条直线L₁，该直线L₁表示为（ρ_L1，θ_L1），Ln为直线集L中除第一条直线L₁外的所有其它直线，直线Ln表示为（ρ_L1，θ_L1），计算第一条直线L1与直线集L中的每一根直线Ln的d_θ=|θ_Ln-θ_L1|，d_ρ=|ρ_Ln-ρ_L1|，将第一条直线L1与直线集L中的每一根直线Ln进行比较，若有d_θ≤30且d_ρ≤5，则每次将直线Ln放入直线L₁的待聚合集中；否则，每次将直线Ln放入非当前待聚合集中；

（3）令所述非当前待聚合集中为直线集L，判断当前的直线集L中的直线数据是否小于等于1，若是，则转到下步骤（4），若否，则转到上述步骤（2）；

（4）分别将每一个直线L₁的待聚合集中的多条直线进行聚合，以获取每一边缘的一条聚合后的边缘直线（ρi’，θi’）作为目标星的直线特征；

步骤20.对二值化后的边缘图像的每个像素点的像素值G(x,y)进行Hough直线检测以获取目标星的每个边缘i的边缘直线（ρi，θi）的步骤中，从第二帧的初始图像开始，减少Hough空间的搜索范围；

步骤21.对二值化后的边缘图像的每个像素点的像素值G(x,y)进行Hough直线检测以获取目标星的每个边缘i的边缘直线（ρi，θi）的步骤之后，对G(x,y)进行harris角点提取以获取目标星的角点特征；

步骤22.对G(x,y)进行harris角点提取以获取目标星的角点特征:对二值化后的边缘图像的每个像素点的像素值G(x,y)进行x、y和xy方向离散求导得Ix、Iy和Ixy:

Ix=-2·G(x-2,y)-G(x-1,y)+G(x+1,y)+2·G(x+2,y)，

Iy=-2·G(x,y-2)-G(x,y-1)+G(x,y+1)+2·G(x,y+2)，

Ixy=Ix·Iy；

对每个Ix和Iy做平方得到Ix²、Iy²；采用5阶高斯滤波器对Ix²、Iy²和Ixy分别进行5阶高斯滤波，得到滤波后的结果Ix2*、Iy2*、Ixy*，其中，5阶高斯滤波器如下:h=[15851；

       52134215；

       83456348；

       52134215；

       15851]；

对二值化后的边缘图像的每个像素点(x,y)求对应的R=（AB-CD）²-k（A+B）²，其中，A=Ix2*(i,j)、B=Iy2*(i,j)、C=D=Ixy*(i,j)，k为响应系数，取0.04；选取角点阈值，如果该像素点对应的R大于所述角点阈值，并且对应的R在该像素点的3×3阶邻域窗口内最大，则将该像素点作为目标星的角点特征作为角点，其中角点阈值为该像素点的3×3阶邻域窗口内最大R值的0.01倍。

有益效果

本发明提出的一种微小卫星非合作目标图像处理方法，通过计算目标星在一帧初始图像中所占的像素个数N，判断像素个数N是否大于一预设阈值，若否，则基于反馈直线分离的方法识别出所述帧的初始图像中目标星的特征；若是，则基于反馈直线聚合的方法识别出所述帧的初始图像中目标星的特征，能够直接提取目标星本身的特征、无需对目标星进行加装改造、对目标星本身无特别要求；相对于合作目标图像处理方法，执行非合作目标图像处理方法的相关硬件结构具有低功耗、小重量；针对目标星相对距离不同自适应的采用不同的特征提取算法、实现大范围特征提取，只需要目标星有适合的光照就可以提取所需的目标星的特征。

附图说明

下面结合附图和实施方式对本发明一种微小卫星非合作目标图像处理方法作进一步详细说明。

图1为本发明微小卫星非合作目标图像处理方法的原理图。

图2为本发明微小卫星非合作目标图像处理方法的流程图。

图3为图2的步骤S2的详细流程图。

图4为4×4阶邻域窗口内插值后的像素点的示意图。

图5为梯度幅值判断的示意图。

图6为图2的步骤S3的详细流程图。

图7为边缘直线进行聚合的原理图。

具体实施方式

本实施例是一种微小卫星非合作目标图像处理方法。

参阅图1-图7，本发明提出的微小卫星非合作目标图像处理方法，包括步骤S1～步骤S3；

步骤S1，计算目标星在一帧初始图像中所占的像素个数N，判断像素个数N是否大于一预设阈值，若否，则转到步骤S2，若是，则转到步骤S3；

计算目标星在相面中所占的像素个数N的步骤中，使用大律法计算目标星在相面中所占的像素个数N。具体的是，使用大津法计算目标所占的像素个数N，利用相机模型计算切换算法的所需的预设阈值N0，当N≤N0时采用的是基于反馈直线分离的特征识别方法，当N>N0时采用的是基于反馈直线聚合的特征识别方法。

步骤S2，基于反馈直线分离的方法识别出所述帧的初始图像中目标星的特征；若像素个数N是否小于等于所述预设阈值，则说明目标星相对追踪时距离较远，则需要执行步骤S2。

优选的，如图3所示，步骤S2包括步骤S21～步骤S24。

步骤S21，对所述帧的初始图像的每个像素点的像素值f(x,y)进行中值滤波得到滤波后的图像的每个像素点的像素值f1(x,y)；本步骤是对获得的灰度图像进行平滑滤波去除噪声。采用的是中值滤波，中值滤波是对图像上每一个像素点的灰度值设置为该点某邻域窗口内的所有像素点灰度值的中值。

优选的，步骤S21包括:获取待进行中值滤波的像素点的第一3×3阶邻域窗口；对所述3×3阶邻域窗口中的9个像素点分别分配一个初值为0的计数器，并将所述第一3×3阶邻域窗口中的每个待比较像素点的像素值分别与其它8个像素点的像素值进行比较，每次若大于其它8个像素点中的一个像素点的像素值，则每次该待比较像素点的计数值+1，否则，每次该待比较像素点的计数值不变；对所述第一3×3阶邻域窗口中的每个像素点的计数器的计数值进行排序，获取排序为中间的计数值所述作为待进行中值滤波的像素点的滤波后的像素值f1(x,y)。具体的是，设f(x,y)、f1(x,y)分别为滤波前后的图像的像素点的像素值，W为3×3阶滤波窗口，则有:f₁(x,y)=med{f(x-k,y-l),(k,l∈W)}，例如，k,l∈0,1,2，其中med为取中值，f(x-k,y-l)为滤波窗口内所有点的灰度值。步骤S21可在FPGA中实现，其实现过程是:在FPGA中，开始第一3×3阶邻域窗口中图像传输时，对第一3×3阶邻域窗口中的原始图的第一行像素数据和第二行像素数据以移位寄存器的方式进行缓存，而对第三行进行2个最新像素的缓存。当第三行的第三个像素到达时，构成3×3结构的缓存数据。对该结构中的9个像素进行排序。其排序过程是:在第一个时钟周期到来时，对9个像素，每个像素分配一个初值为0的计数器，并将该像素分别与其它8个像素进行比较，若大于，则该像素的计数值+1，否则计数值不变。在第二个时钟周期到来时，对每个像素计数器进行统计。统计后，取排序后中间的计数值对应的像素点的像素值作为输出像素值。上述过程中，两个时钟的计算过程通过流水线实现，通过这种方式，保证每当有一个图像像素数据达到时，都会有一个计算后得到的中值输出。

步骤S22，对滤波后的图像进行双线性插值得到插值后的图像；当目标星距离较远时，目标星在图像上占的像素较少，提取的特征直线较少，当检测精度较低时，一些特征直线可能检测不出，这里对灰度图像进行插值，提高图像特征直线的检测精度，从而可分辨出更多的特征直线。

优选的，步骤S22包括:将滤波后的图像分割为若干个4×4阶邻域窗口；

根据公式 $\begin{array}{c}{f}_2(x+u,y+v)=(1-u)(1-v)f_1(x,y)+(1-u){\mathrm{vf}}_1(x,y+1)+\\ u(1-v)f_1(x+1,y)+{\mathrm{uvf}}_1(x+1,y+1)\end{array}$

获取每个4×4阶邻域窗口内插值后的像素点的像素值f2(x,y)，其中u，v是[0，1]区间的浮点数，本步骤中在每个4×4邻域增加5个像素，（u，v）取值分别为(0，0.5)、(0.5，0)、(0.5，0.5)、(0.5，1)和(1，0.5)，

如图4所示，其中点1、3、7、9为滤波后的点f₁(i,j)、f₁(i,j+1)、f₁(i+1,j)、f₁(i+1,j+1)，其余点为插值后的点:

f₂(x+0.5,y)=0.5f₁(i,j)+0.5f₁(i,j+1)、

f₂(x,y+0.5)=0.5f₁(i,j)+0.5f₁(i+1,j)、

f₂(x+0.5,y+0.5)=0.5f₁(i,j)+0.5f₁(i,j+1)+0.5f₁(i+1,j)+0.5f₁(i+1,j+1)、

f₂(x+1,y)=0.5f₁(i,j+1)+0.5f₁(i+1,j+1)、

f₂(x+0.5,y+1)=0.5f₁(i+1,j)+0.5f₁(i+1,j+1)，

最后将得到的插值后的图像中的每个滤波后的像素点和插值后的点的像素值统一计为f₂’(x,y)。

步骤S23，对插值后的图像的每个像素点进行sobel算子滤波、非极大抑制、二值化得到二值化后的边缘图像的每个像素点的像素值G(x,y)；具体的是，本步骤即是对插值后的图像f₂’(x,y)上每一个点进行灰度梯度计算，对其结果进行二值化提取图像的边缘信息，本实施例中可采用的是改进的canny算法。在步骤S23中，对插值后的图像的每个像素点进行sobel算子滤波的步骤包括:

将每个待进行sobel算子滤波的像素点的第二3×3阶邻域窗口；

在第二3×3阶邻域窗口内根据公式g(x,y)=|g_x(x,y)+g_y(x,y)|获取每个待进行sobel算子滤波的像素点的梯度幅值g(x,y)，其中，

g_x(x,y)=-f′₂(x-1,y-1)-2·f′₂(x-1,y)-f′₂(x-1,y+1)+

f′₂(x+1,y-1)+2·f′₂(x+1,y)+f′₂(x+1,y+1)，

g_y(x,y)=f′₂(x-1,y-1)+2·f′₂(x,y-1)+f′₂(x+1,y-1)-

f′₂(x-1,y+1)-2·f′₂(x,y+1)-f′₂(x+1,y+1)，

g_x(x,y)，g_y(x,y)分别为每个待进行sobel算子滤波的像素点的X方向和Y方向的梯度幅值。本步骤可在FPGA中实现，其实现过程是:在FPGA中，对第二3×3阶邻域窗口内的第一行像素数据和第二行像素数据以移位寄存器的方式进行缓存，而对第三行进行2个最新像素的缓存。当第三行的第三个像素到达时，构成3×3结构的缓存数据，代入式中求得g(x,y)。使用这种方式，每个像素周期得到一个计算后的输出结果。

优选的步骤S23中，对插值后的图像的每个像素点进行非极大抑制的步骤包括:根据sobel算子滤波的结果判断插值后的图像的每个像素点C的梯度幅值是否大于梯度方向上相邻的两点dTmp1和dTmp2的梯度幅值，若是，则像素点C的梯度幅值不变；若否，将该像素点C的梯度幅值置0。本步骤为边缘细化过程，可以使边缘更加清晰细致。每一个点的梯度方向进行梯度幅值比较，如果此点的梯度比它相邻的两个点梯度幅值都大，则不变；否则将其置0并排除出边缘点。C点在梯度方向上相邻的两点为dTmp1、dTmp2，则极大值抑制结果为:如果C点梯度幅值比dTmp1、dTmp2点都大，则C点梯度幅值不变；否则，将C点梯度幅值置0。其中，dTmp1、dTmp2点的梯度幅值用插值由g1g2g3g4的灰度值求出:

dTmp1=w·g1+(1-w)·g2

dTmp2=w·g3+(1-w)·g4

其中w为权重，其计算公式为:

当g_x(x,y)＞g_y(x,y)时，

当g_x(x,y)＜g_y(x,y)时，

在FPGA中用g(x,y)和dTmp1,dTmp2进行比较大小时，为避免进行除法运算，将两边同时乘以分母，用三个乘法器在一个周期内得到三个相乘结果，第二个周期完成比较，得到实时的非极大值抑制后的图像数据。

步骤S23中，对插值后的图像的每个像素点进行二值化得到二值化后的边缘图像的所有像素点的像素值G(x,y)的步骤包括:判断非极大抑制后的每一个像素点是否大于一预设的梯度幅值阈值T₁，若是，则将该像素点的梯度幅值置为1，若否，则将该像素点的梯度幅值置为0；

提取每个梯度幅值为1像素点的像素值作为二值化后的边缘图像的像素点的像素值G(x,y)。此步骤的目的是提高边缘和其它位置的对比度，突出边缘。步骤为设定一个梯度幅值阈值，对每一个像素点如果其梯度幅值大于阈值则此像素点梯度幅值设置为1，否则设置为0，最终得到边缘提取结果G(x,y)。梯度幅值阈值T1的选择可采用计梯度幅值直方图的方法计算:首先对整个图像梯度幅值进行统计，按照从小到大排列，设梯度幅值范围为Gmin～Gmax，其中Gmin为最小梯度幅值，Gmax为最大梯度幅值；然后计算梯度幅值阈值为:

T₁=（Gmax-Gmin）×0.79

此步骤可在FPGA中实现，实现过程是:在每一帧的数据传输过程中，设置两个变量Gmin和Gmax。每个像素到达时，对其判断是否大于Gmax，是否小于Gmin。若是，更新相应的Gmax或Gmin的值，否则Gmax和Gmin保持不变。完成图像传输后，计算T₁值，作为下一帧的梯度幅值阈值。

步骤S24，对二值化后的边缘图像的每个像素点的像素值G(x,y)进行Hough直线检测以获取目标星的每个边缘i的边缘直线（ρi，θi）作为目标星的直线特征，i为正整数。

步骤S24包括:设定Hough平面累加器（ρi，θi），其中θi从0到180°，θi的步长为1°，ρi为

其中，w为图像宽度、h为图像高度，ρi的步长为一个像素；对二值化后的边缘图像的每个像素点即所有非0点进行Hough变换:依次将θi从0到180°带入直线极坐标方程ρi=xcosθi+ysinθi中计算，每次计算得到的（ρi，θi）在累加器的相应单元的计数值+1；

Hough变换处理完成后，设定阈值T_hough，累加器中大于T_hough的单元对应的（ρi，θi）作为检测出的边缘直线的直线特征，其中T_hough的选取采用的是直方图形式，取所有累加器单元的值，其范围为Amin～Amax，则有T_hough=0.5·(Amax-Amin)，Amin、Amax分别为所述累加器的所有单元的计算值的最小值和最大值。

在步骤S24中，从第二帧的初始图像开始，减少Hough空间的搜索范围。例如，从第二帧开始，Hough空间可减少为（ρi±3）及（θi±10）。具体的是，对图像进行插值在第一帧时，Hough参数空间为全图全角度空间，检测到的特征直线根据直线参数不同可将其分离、识别。从第二帧开始，对于每个边缘点其Hough空间可由上一帧反馈的特征直线参数对Hough空间进行压缩、减少，设上一帧反馈的特征直线参数为ρiθi，此时Hough空间为（ρi±3）及（θi±10）。

步骤S24之后还包括对G(x,y)进行harris角点提取以获取目标星的角点特征。G(x,y)进行harris角点提取以获取目标星的角点特征的步骤包括:对二值化后的边缘图像的每个像素点的像素值G(x,y)进行x、y和xy方向离散求导得Ix、Iy和Ixy:

Ix=-2·G(x-2,y)-G(x-1,y)+G(x+1,y)+2·G(x+2,y)，

Iy=-2·G(x,y-2)-G(x,y-1)+G(x,y+1)+2·G(x,y+2)，

Ixy=Ix·Iy；

对每个Ix和Iy做平方得到Ix²、Iy²；采用5阶高斯滤波器对Ix²、Iy²和Ixy分别进行5阶高斯滤波，得到滤波后的结果:Ix²*、Iy²*、Ixy*，其中，5阶高斯滤波器如下:

         h=[15851；

            52134215；

            83456348；

            52134215；

            15851]；

对二值化后的边缘图像的每个像素点(x,y)求对应的R=（AB-CD）²-k（A+B）²，其中，A=Ix²*(i,j)、B=Iy²*(i,j)、C=D=Ixy*(i,j)，k为响应系数，取0.04；

选取角点阈值T₂，如果该像素点对应的R大于所述角点阈值T₂，并且对应的R在该像素点的3×3阶邻域窗口内最大，则将该像素点作为目标星的角点特征作为角点，其中，角点阈值T₂为该像素点的3×3阶邻域窗口内最大R值的0.01倍。即满足下式即为角点:

R(i,j)>0.01*Rmax&&R(i,j)>R(i-1,j-1)&&

R(i,j)>R(i-1,j)&&R(i,j)>R(i-1,j+1)&&

R(i,j)>R(i,j-1)&&R(i,j)>R(i,j+1)&&

R(i,j)>R(i+1,j-1)&&R(i,j)>R(i+1,j)&&R(i,j)>R(i+1,j+1)

步骤S3，基于反馈直线聚合的方法识别出所述帧的初始图像中目标星的特征。若像素个数N是否大于所述预设阈值，则说明目标星相对追踪时距离较近，则需要执行步骤S3。例如，当目标星在相面所占像素大于50×50的预设阈值时，由分离算法切换为直线聚合的特征识别算法。

如图6所示，步骤S3包括步骤S31～步骤S34。

步骤S31，对所述帧的初始图像的每个像素点的像素值f(x,y)进行中值滤波得到滤波后的图像的每个像素点的像素值f1(x,y)；

步骤S31包括:获取待进行中值滤波的像素点的第一3×3阶邻域窗口；对所述3×3阶邻域窗口中的9个像素点分别分配一个初值为0的计数器，并将所述第一3×3阶邻域窗口中的每个待比较像素点的像素值分别与其它8个像素点的像素值进行比较，每次若大于其它8个像素点中的一个像素点的像素值，则每次该待比较像素点的计数值+1，否则，每次该待比较像素点的计数值不变；对所述第一3×3阶邻域窗口中的每个像素点的计数器的计数值进行排序，获取排序为中间的计数值作为待进行中值滤波的像素点的滤波后的像素值f1(x,y)。

步骤S32，对滤波后的图像的每个像素点进行sobel算子滤波、非极大抑制、二值化得到二值化后的边缘图像的每个像素点的像素值G(x,y)；

优选的，步骤S32中，对滤波后的图像的每个像素点进行sobel算子滤波的步骤包括:将每个待进行sobel算子滤波的像素点的第二3×3阶邻域窗口；在第二3×3阶邻域窗口内根据公式g(x,y)=|g_x(x,y)+g_y(x,y)|获取每个待进行sobel算子滤波的像素点的梯度幅值g(x,y)，其中，

g_x(x,y)=-f₁(x-1,y-1)-2·f₁(x-1,y)-f₁(x-1,y+1)+

f₁(x+1,y-1)+2·f₁(x+1,y)+f₁(x+1,y+1)，

g_y(x,y)=f₁(x-1,y-1)+2·f₁(x,y-1)+f₁(x+1,y-1)-

f₁(x-1,y+1)-2·f₁(x,y+1)-f₁(x+1,y+1)，

g_x(x,y)，g_y(x,y)分别为每个待进行sobel算子滤波的像素点的X方向和Y方向的梯度幅值。

优选的，步骤S32中，对滤波后的图像的每个像素点进行非极大抑制的步骤包括:根据sobel算子滤波的结果判断插值后的图像的每个像素点C的梯度幅值是否大于梯度方向上相邻的两点dTmp1和dTmp2的梯度幅值，若是，则像素点C的梯度幅值不变；若否，将该像素点C的梯度幅值置0。

对滤波后的图像的每个像素点进行二值化得到二值化后的边缘图像的每个像素点的像素值G(x,y)的步骤包括:判断非极大抑制后的每一个像素点是否大于一预设的梯度幅值阈值，若是，则将该像素点的梯度幅值置为1，若否，则将该则此像素点的梯度幅值置为0；提取每个梯度幅值为1像素点的像素值作为二值化后的边缘图像的像素点的像素值G(x,y)。

步骤S33，对二值化后的边缘图像的每个像素点的像素值G(x,y)进行Hough直线检测以获取目标星的每个边缘i的边缘直线（ρi，θi），i为正整数；

优选的，步骤S33包括:设定Hough平面累加器（ρi，θi），其中，θi从0到180°，θi的步长为1°，ρi为

其中，w为图像宽度、h为图像高度，ρi的步长为一个像素；对二值化后的边缘图像的每个像素点进行Hough变换:依次将θi从0到180°带入直线极坐标方程ρi=xcosθi+ysinθi中计算，每次计算得到的（ρi，θi）在累加器的相应单元的计数值+1；

Hough变换处理完成后，设定阈值T_hough，累加器中大于T_hough的单元对应的（ρi，θi）作为检测出的边缘直线，其中，T_hough=0.5·(Amax-Amin)，Amin、Amax分别为所述累加器的所有单元的计算值的最小值和最大值。

优选的，步骤S33，从第二帧的初始图像开始，减少Hough空间的搜索范围。例如，从第二帧开始，Hough空间可减少为（ρi±3）及（θi±10）；

步骤S33之后还包括:对G(x,y)进行harris角点提取以获取目标星的角点特征。对G(x,y)进行harris角点提取以获取目标星的角点特征的步骤包括:对二值化后的边缘图像的每个像素点的像素值G(x,y)进行x、y和xy方向离散求导得Ix、Iy和Ixy:

Ix=-2·G(x-2,y)-G(x-1,y)+G(x+1,y)+2·G(x+2,y)，

Iy=-2·G(x,y-2)-G(x,y-1)+G(x,y+1)+2·G(x,y+2)，

Ixy=Ix·Iy；

对每个Ix和Iy做平方得到Ix²、Iy²；采用5阶高斯滤波器对Ix²、Iy²和Ixy分别进行5阶高斯滤波，得到滤波后的结果:Ix2*、Iy2*、Ixy*，其中，5阶高斯滤波器如下:

       h=[15851；

          52134215；

          83456348；

          52134215；

          15851]；

对二值化后的边缘图像的每个像素点(x,y)求对应的R=（AB-CD）²-k（A+B）²，其中，A=Ix2*(i,j)、B=Iy2*(i,j)、C=D=Ixy*(i,j)，k为响应系数，取0.04；

选取角点阈值，如果该像素点对应的R大于所述角点阈值，并且对应的R在该像素点的3×3阶邻域窗口内最大，则将该像素点作为目标星的角点特征作为角点，其中角点阈值为该像素点的3×3阶邻域窗口内最大R值的0.01倍。

步骤S34，对目标星的表示同一边缘i的多条边缘直线（ρi，θi）进行聚合以获取一条聚合后的边缘直线（ρi’，θi’）作为目标星的直线特征。具体的是，当目标星距离较近时，目标星在图像上占的像素较多，提取的特征直线较多，一个边缘可能检测出多条特征直线，需要采用特征直线的聚合算法。聚合算法是将同一边缘的直线归为一组、聚合放入一个集合中，从而能够识别所检测的直线，并且输出所需要的特征直线。这里将同一边缘上的直线当成一类从所有直线中选出，然后对这些直线拟合输出一条直线作为边缘的直线。这种方法既能够识别直线，也可以提高直线的提取精度。

优选的，步骤S34包括:

步骤一，令从所述Hough直线检测出的所有边缘直线（ρi，θi）作为直线集L；

步骤二，从直线集L中取出第一条直线L₁，该直线L₁表示为（ρ_L1，θ_L1），Ln为直线集L中除第一条直线L₁外的所有其它直线，n为正整数，直线Ln表示为（ρ_L1，θ_L1），计算第一条直线L₁与直线集L中的每一根直线Ln的d_θ=|θ_Ln-θ_L1|， d_ρ=|ρ_Ln-ρ_L1|，将第一条直线L₁与直线集L中的每一根直线Ln进行比较，若有d_θ≤30且d_ρ≤5，则每次将直线Ln放入直线L₁的待聚合集中；否则，每次将直线Ln放入非当前待聚合集中；

步骤三，令所述非当前待聚合集L₀中为直线集L，判断当前的直线集L中的直线数据是否小于等于1，若是，则转到步骤四，若否，则转到所述步骤二；

步骤四，分别将每一个直线L1的待聚合集中的多条直线进行聚合，以获取每一边缘的一条聚合后的边缘直线（ρi’，θi’）作为目标星的直线特征；

其中直线集L为检测到的所有直线，每输出一类边缘直线时，此直线集L的直线数减少；L₁是一类边缘的直线集，这些直线拟合可以得到属于同一边缘的特征直线；L₀是不属于当前边缘的直线。步骤S34的原理如图7所示。

本发明通过计算目标星在一帧初始图像中所占的像素个数N，判断像素个数N是否大于一预设阈值，若否，则基于反馈直线分离的方法识别出所述帧的初始图像中目标星的特征；若是，则基于反馈直线聚合的方法识别出所述帧的初始图像中目标星的特征，能够直接提取目标星本身的特征、无需对目标星进行加装改造、对目标星本身无特别要求；相对于合作目标图像处理方法，执行非合作目标图像处理方法的相关硬件结构具有功耗低、重量小；针对目标星相对距离不同自适应的采用不同的特征提取算法、实现大范围特征提取，只需要目标星有适合的光照就可以提取所需的目标星的特征。

Claims (1)

1.一种微小卫星非合作目标图像处理方法，其特征在于包括以下步骤:

步骤1.计算目标星在一帧初始图像中所占的像素个数N，判断像素个数N是否大于一预设阈值，若否，则基于反馈直线分离的方法识别出初始图像中目标星的特征；若是，则基于反馈直线聚合的方法识别出初始图像中目标星的特征；

步骤2.计算目标星在相面中所占的像素个数N，采用大律法计算目标星在相面中所占的像素个数N；

步骤3.基于反馈直线分离的方法识别出初始图像中目标星的特征:对初始图像的每个像素点的像素值f(x,y)进行中值滤波得到滤波后的图像的每个像素点的像素值f1(x,y)；对滤波后的图像进行双线性插值得到插值后的图像；对插值后的图像的每个像素点进行sobel算子滤波、非极大抑制、二值化得到二值化后的边缘图像的每个像素点的像素值G(x,y)；对二值化后的边缘图像的每个像素点的像素值G(x,y)进行Hough直线检测以获取目标星的每个边缘i的边缘直线（ρi，θi）作为目标星的直线特征，i为正整数；

步骤4.对初始图像的每个像素点的像素值f(x,y)进行中值滤波，得到滤波后的图像的每个像素点的像素值f1(x,y):获取待进行中值滤波的像素点的第一3×3阶邻域窗口；对所述3×3阶邻域窗口中的9个像素点分别分配一个初值为0的计数器，并将所述第一3×3阶邻域窗口中的每个待比较像素点的像素值分别与其它8个像素点的像素值进行比较，每次若大于其它8个像素点中的一个像素点的像素值，则每次该待比较像素点的计数值+1，否则，每次该待比较像素点的计数值不变；对所述第一3×3阶邻域窗口中的每个像素点的计数器的计数值进行排序，获取排序为中间的计数值作为待进行中值滤波的像素点的滤波后的像素值f1(x,y)；

步骤5.对滤波后的图像进行双线性插值得到插值后的图像:将滤波后的图像分割为若干个4×4阶邻域窗口；

根据公式 $\begin{array}{c}{f}_2(x+u,y+v)=(1-u)(1-v)f_1(x,y)+(1-u){\mathrm{vf}}_1(x,y+1)+\\ u(1-v)f_1(x+1,y)+{\mathrm{uvf}}_1(x+1,y+1)\end{array}$ 获取每个4×4阶邻域窗口内插值后的点的像素值f₂(x,y)，其中u，v是[0，1]区间的浮点数，（u，v）取值分别为(0，0.5)、(0.5，0)、(0.5，0.5)、(0.5，1)、(1，0.5)，将得到的插值后的图像中的每个滤波后的像素点和插值后的点的像素值统一计为f₂’(x,y)；

步骤6.对插值后的图像的每个像素点进行sobel算子滤波:在第二3×3阶邻域窗口内根据公式g(x,y)=|g_x(x,y)+g_y(x,y)|获取每个待进行sobel算子滤波的像素点的梯度幅值g(x,y)，其中，

g_x(x,y)=-f′₂(x-1,y-1)-2·f′₂(x-1,y)-f′₂(x-1,y+1)+

f′₂(x+1,y-1)+2·f′₂(x+1,y)+f′₂(x+1,y+1)，

g_y(x,y)=f′₂(x-1,y-1)+2·f′₂(x,y-1)+f′₂(x+1,y-1)-

f′₂(x-1,y+1)-2·f′₂(x,y+1)-f′₂(x+1,y+1)，

g_x(x,y)，g_y(x,y)分别为每个待进行sobel算子滤波的像素点的X方向和Y方向的梯度幅值；

步骤7.对插值后的图像的每个像素点进行非极大抑制:根据sobel算子滤波的结果判断插值后的图像的每个像素点C的梯度幅值是否大于梯度方向上相邻的两点dTmp1和dTmp2的梯度幅值，若是，则像素点C的梯度幅值不变；若否，将该像素点C的梯度幅值置0；

步骤8.对插值后的图像的每个像素点进行二值化得到二值化后的边缘图像的所有像素点的像素值G(x,y):判断非极大抑制后的每一个像素点是否大于预设的梯度幅值阈值，若是，则将该像素点的梯度幅值置为1，若否，则将该像素点的梯度幅值置为0；提取每个梯度幅值为1像素点的像素值作为二值化后的边缘图像的像素点的像素值G(x,y)；

步骤9.对二值化后的边缘图像的每个像素点的像素值G(x,y)进行Hough直线检测以获取目标星的每个边缘i的边缘直线（ρi，θi）作为目标星的直线特征:设Hough平面累加器（ρi，θi），其中θi从0到180°，θi的步长为1°，ρi为

其中，w为图像宽度、h为图像高度，ρi的步长为一个像素；对二值化后的边缘图像的每个像素点进行Hough变换：依次将θi从0到180°带入直线极坐标方程ρi=xcosθi+ysinθi中计算，每次计算得（ρi，θi）在累加器的相应单元的计数值+1；Hough变换处理完成后，设定阈值T_hough，累加器中大于T_hough的单元对应的（ρi，θi）作为检测出的边缘直线的直线特征，其中，T_hough=0.5·(Amax-Amin)，Amin、Amax分别为累加器的所有单元的计算值的最小值和最大值；

步骤10.对二值化后的边缘图像的每个像素点的像素值G(x,y)进行Hough直线检测以获取目标星的每个边缘i的边缘直线（ρi，θi）作为目标星的直线特征，从第二帧的初始图像开始，减少Hough空间的搜索范围；

步骤11.对二值化后的边缘图像的每个像素点的像素值G(x,y)进行Hough直线检测以获取目标星的每个边缘i的边缘直线（ρi，θi）作为目标星的直线特征，并对G(x,y)进行harris角点提取以获取目标星的角点特征；

步骤12.对G(x,y)进行harris角点提取以获取目标星的角点特征:对二值化后的边缘图像的每个像素点的像素值G(x,y)进行x、y和xy方向离散求导得Ix、Iy和Ixy:

Ix=-2·G(x-2,y)-G(x-1,y)+G(x+1,y)+2·G(x+2,y)，

Iy=-2·G(x,y-2)-G(x,y-1)+G(x,y+1)+2·G(x,y+2)，

Ixy=Ix·Iy；

对每个Ix和Iy做平方得到Ix²、Iy²；采用5阶高斯滤波器对Ix²、Iy²和Ixy分别进行5阶高斯滤波，得到滤波后的结果:Ix²*、Iy²*、Ixy*，其中，5阶高斯滤波器如下:h=[15851；

          52134215；

          83456348；

          52134215；

          15851]；

对二值化后的边缘图像的每个像素点(x,y)求对应的R=（AB-CD）²-k（A+B）²，其中，A=Ix²*(i,j)、B=Iy²*(i,j)、C=D=Ixy*(i,j)，k为响应系数，取0.04；选取角点阈值，如果该像素点对应的R大于所述角点阈值，并且对应的R在该像素点的3×3阶邻域窗口内最大，则将该像素点作为目标星的角点特征作为角点，其中角点阈值为该像素点的3×3阶邻域窗口内最大R值的0.01倍；

步骤13.基于反馈直线聚合的方法识别出所述帧的初始图像中目标星的特征:对所述帧的初始图像的每个像素点的像素值f(x,y)进行中值滤波得到滤波后的图像的每个像素点的像素值f1(x,y)；对滤波后的图像的每个像素点进行sobel算子滤波、非极大抑制、二值化得到二值化后的边缘图像的每个像素点的像素值G(x,y)；对二值化后的边缘图像的每个像素点的像素值G(x,y)进行Hough直线检测以获取目标星的每个边缘i的边缘直线（ρi，θi），i为正整数；对目标星的表示同一边缘i的多条边缘直线（ρi，θi）进行聚合以获取一条聚合后的边缘直线（ρi’，θi’）作为目标星的直线特征；

步骤14.对所述帧的初始图像的每个像素点的像素值f(x,y)进行中值滤波得到滤波后的图像的每个像素点的像素值f1(x,y):获取待进行中值滤波的像素点的第一3×3阶邻域窗口；对所述3×3阶邻域窗口中的9个像素点分别分配一个初值为0的计数器，并将所述第一3×3阶邻域窗口中的每个待比较像素点的像素值分别与其它8个像素点的像素值进行比较，每次若大于其它8个像素点中的一个像素点的像素值，则每次该待比较像素点的计数值+1，否则，每次该待比较像素点的计数值不变；对所述第一3×3阶邻域窗口中的每个像素点的计数器的计数值进行排序，获取排序为中间的计数值作为待进行中值滤波的像素点的滤波后的像素值f1(x,y)；

步骤15.对滤波后的图像的每个像素点进行sobel算子滤波:将每个待进行sobel算子滤波的像素点的第二3×3阶邻域窗口；在第二3×3阶邻域窗口内根据公式g(x,y)=|g_x(x,y)+g_y(x,y)|获取每个待进行sobel算子滤波的像素点的梯度幅值g(x,y)，其中，

g_x(x,y)=-f₁(x-1,y-1)-2·f₁(x-1,y)-f₁(x-1,y+1)+

f₁(x+1,y-1)+2·f₁(x+1,y)+f₁(x+1,y+1)，

g_y(x,y)=f₁(x-1,y-1)+2·f₁(x,y-1)+f₁(x+1,y-1)-

f₁(x-1,y+1)-2·f₁(x,y+1)-f₁(x+1,y+1)，

g_x(x,y)，g_y(x,y)分别为每个待进行sobel算子滤波的像素点的X方向和Y方向的梯度幅值；

步骤16.对滤波后的图像的每个像素点进行非极大抑制:根据sobel算子滤波的结果判断插值后的图像的每个像素点C的梯度幅值是否大于梯度方向上相邻的两点dTmp1和dTmp2的梯度幅值，若是，则像素点C的梯度幅值不变；若否，将该像素点C的梯度幅值置0；

步骤17.对滤波后的图像的每个像素点进行二值化得到二值化后的边缘图像的每个像素点的像素值G(x,y):判断非极大抑制后的每一个像素点是否大于一预设的梯度幅值阈值，若是，则将该像素点的梯度幅值置为1，若否，则将该像素点的梯度幅值置为0；提取每个梯度幅值为1像素点的像素值作为二值化后的边缘图像的像素点的像素值G(x,y)；

步骤18.对二值化后的边缘图像的每个像素点的像素值G(x,y)进行Hough直线检测以获取目标星的每个边缘i的边缘直线（ρi，θi）:设定Hough平面累加器（ρi，θi），其中θi从0到180°，θi的步长为1°，ρi为

其中，w为图像宽度、h为图像高度，ρi的步长为一个像素；对二值化后的边缘图像的每个像素点进行Hough变换：依次将θi从0到180°带入直线极坐标方程ρi=xcosθi+ysinθi中计算，每次计算的（ρi，θi）在累加器的相应单元的计数值+1；Hough变换处理完成后，设定阈值T_hough，累加器中大于T_hough的单元对应的（ρi，θi）作为检测出的边缘直线，其中，T_hough=0.5·(Amax-Amin)，Amin、Amax分别为累加器的所有单元的计算值的最小值和最大值；

步骤19.对目标星的表示同一边缘i的多条边缘直线（ρi，θi）进行聚合以获取一条聚合后的边缘直线（ρi’，θi’）作为目标星的直线特征:

（1）令从所述Hough直线检测出的所有边缘直线（ρi，θi）作为直线集L；

（2）从直线集L中取出第一条直线L₁，该直线L₁表示为（ρ_L1，θ_L1），Ln为直线集L中除第一条直线L₁外的所有其它直线，直线Ln表示为（ρ_L1，θ_L1），计算第一条直线L1与直线集L中的每一根直线Ln的d_θ=|θ_Ln-θ_L1|，d_ρ=|ρ_Ln-ρ_L1|，将第一条直线L1与直线集L中的每一根直线Ln进行比较，若有d_θ≤30且d_ρ≤5，则每次将直线Ln放入直线L₁的待聚合集中；否则，每次将直线Ln放入非当前待聚合集中；

（3）令所述非当前待聚合集中为直线集L，判断当前的直线集L中的直线数据是否小于等于1，若是，则转到下步骤（4），若否，则转到上述步骤（2）；

（4）分别将每一个直线L₁的待聚合集中的多条直线进行聚合，以获取每一边缘的一条聚合后的边缘直线（ρi’，θi’）作为目标星的直线特征；

步骤20.对二值化后的边缘图像的每个像素点的像素值G(x,y)进行Hough直线检测以获取目标星的每个边缘i的边缘直线（ρi，θi）的步骤中，从第二帧的初始图像开始，减少Hough空间的搜索范围；

步骤21.对二值化后的边缘图像的每个像素点的像素值G(x,y)进行Hough直线检测以获取目标星的每个边缘i的边缘直线（ρi，θi）的步骤之后，对G(x,y)进行harris角点提取以获取目标星的角点特征；

步骤22.对G(x,y)进行harris角点提取以获取目标星的角点特征:对二值化后的边缘图像的每个像素点的像素值G(x,y)进行x、y和xy方向离散求导得Ix、Iy和Ixy:

Ix=-2·G(x-2,y)-G(x-1,y)+G(x+1,y)+2·G(x+2,y)，

Iy=-2·G(x,y-2)-G(x,y-1)+G(x,y+1)+2·G(x,y+2)，

Ixy=Ix·Iy；

对每个Ix和Iy做平方得到Ix²、Iy²；采用5阶高斯滤波器对Ix²、Iy²和Ixy分别进行5阶高斯滤波，得到滤波后的结果Ix2*、Iy2*、Ixy*，其中，5阶高斯滤波器如下:h=[15851；

              52134215；

              83456348；

              52134215；

              15851]；

对二值化后的边缘图像的每个像素点(x,y)求对应的R=（AB-CD）²-k（A+B）²，其中，A=Ix2*(i,j)、B=Iy2*(i,j)、C=D=Ixy*(i,j)，k为响应系数，取0.04；选取角点阈值，如果该像素点对应的R大于所述角点阈值，并且对应的R在该像素点的3×3阶邻域窗口内最大，则将该像素点作为目标星的角点特征作为角点，其中角点阈值为该像素点的3×3阶邻域窗口内最大R值的0.01倍。

Images