CN107449416A

CN107449416A - 基于矢量累积的恒星拖尾星点提取方法

Info

Publication number: CN107449416A
Application number: CN201710467087.7A
Authority: CN
Inventors: 于起峰; 候旺; 雷志辉; 朱宪伟; 包云娜; 尚洋
Original assignee: National University of Defense Technology
Current assignee: National University of Defense Technology
Priority date: 2017-06-20
Filing date: 2017-06-20
Publication date: 2017-12-08
Anticipated expiration: 2037-06-20
Also published as: CN107449416B

Abstract

本发明涉及基于矢量累积的恒星拖尾星点提取方法。首先，从原始图像中提取特征图像；然后分别计算拖尾星点的旋转角和长度；之后，将整幅星图沿着拖尾建造矢量以拖尾步长进行累积，提高拖尾星点的信噪比，得到能量增强的目标，获得结果图像；最后对恒星目标通过传统质心法对目标进行检测并进行位置估计。本发明能够有效提高对拖尾星点定位精度；可以增强拖尾星点的信噪比和极限星等；在飞行器快速转动条件下可有效工作；在星敏感器非线性转动的条件下可稳定工作。

Description

基于矢量累积的恒星拖尾星点提取方法

技术领域

本发明属于天文导航系统技术领域，针对拖尾现象提高星点定位精度解算飞行器姿态，获得高精度的姿态信息。

背景技术

天文导航系统是一种重要的自主导航系统。对比其它导航方法，这种方法有四个优点：第一，天文导航系统被动接受天体(例如，太阳，月亮，卫星，或其它恒星)辐射的信号。于是，它就成为一种自主的和被动的导航方法。第二，天文导航系统接受全波段信号。因此天文导航有着良好的抗干扰能力。这是由于自然天体的运动规律不变，所以不存在累积误差。第三，由于天文导航系统不受空间和时间的限制，因此有着广泛的应用。第四，天文导航系统设备比较简单，价格便宜。在天文导航系统中，星敏感器被用于获取精确的姿态信息。星敏感器的内部处理算法包括图像预处理、星点位置估计、星图识别、星图匹配和姿态估计。星点位置估计是星图识别和匹配的基础。只有在星点位置估计准确且可靠的条件下，星敏感器才能输出高精度的姿态信息。在动态条件下，星点会产生图像拖尾现象。图像拖尾现象对飞行器姿态估计有非常重要的影响。

在星点位置估计中图像拖尾主要产生两个问题：首先，快速的飞行器转动导致信噪比明显降低，导致星点定位精度和极限星等也会下降。第二，飞行器非线性转动导致姿态估计偏差。因此，抑制图像拖尾对高精度的星点位置估计非常关键。

按照传统方法解算，飞行器姿态是通过拖尾星点质心进行解算的。传统方法认为此姿态对应的时刻为曝光中心时刻。但是，当飞行器进行非匀速转动时，一般条件下，曝光中心时刻的星点位置不为拖尾星点质心位置。于是，估计的姿态与时刻没有一一对应关系，会导致姿态估计的偏差。

发明内容

为此，本发明特提出一种基于矢量累积的拖尾星点提取方法，针对星敏感器转动条件下恒星目标信噪比下降导致的定姿精度下降问题，以及拖尾恒星加速转动条件下提点精度误差增加的问题，能够有效提高对深空背景下的星点定位精度。

为达到上述目的，本发明采用如下技术方案：首先，从原始图像中提取特征图像；然后分别计算拖尾星点的旋转角和长度；之后，通过矢量累积算法，将整幅星图沿着拖尾建造矢量以拖尾步长进行累积，提高拖尾星点的信噪比，得到能量增强的目标，获得结果图像；最后对恒星目标通过传统质心法对目标进行检测并进行位置估计。

其中：为计算拖尾星点的旋转角和长度，建立基于线性拖尾模型的星点数学模型：

本发明从姿态方程W＝AV出发分析拖尾星点的形状，提出一种分析拖尾星点的线性模型，建立拖尾星点的数学模型。其中A是星敏感器输出的飞行器姿态矩阵，W和V分别表示体坐标系和惯性坐标系下的单位矢量。W和V的值可由以下方程算得：

其中(α,δ)表示恒星在天球上的赤经和赤纬，(x,y)表示星点在焦平面上位置，f表示相机的焦距；

设(u,v)表示星点在焦平面上的速度，t为时间，而相机的焦距与时间无关，因此，在大多数情况下，(x²+y²)/f小于0.01。于是，以时间为自变量对姿态方程求导可得近似方程：

在视场较小和飞行器姿态变化较小条件下，星点在焦平面上每个位置运动速度可以认为近似相等，即拖尾星点形状在焦平面上每个位置都是近似相同的。这些发现证明了单帧星图上的所有星点拖尾方向和长度都是相同的。因此，提出一种分析拖尾星点数学模型。

建立在二维高斯函数模型以及拖尾星点成因的基础上，提出一种基于线性拖尾模型的星点数学模型如下：

其中I₀表示星敏感器上的星点总能量，(-x₀,0)为曝光开始时焦平面上星点位置，(x₀,0)为曝光结束时焦平面上星点位置，σ为星点的高斯半径，拖尾星点的中心位于坐标原点，在曝光时间沿着x轴方向移动。这个模型为标准星点数学模型；

上面公式可以写成两个独立函数的乘积形式：I(x,y)＝I₁(x)I₂(x)，其中：

星点数学模型的一个例子如附图2所示。在附图2(a)和(b)中分别显示了线性拖尾星点模型的三维曲面和二维投影图像。附图2(c)为I₁(x)函数形状，附图2(d)为I₂(y)函数形状。由于实际拖尾星点存在平移和旋转，因此，我们可通过标准拖尾模型的几何变换(旋转和平移)来描述实际拖尾星点。我们的模型是一个在小视场和飞行器姿态变化较小条件下的线性模型，并且此线性模型有明确的表达式。

矢量累积算法

本发明提供矢量累积算法估计星点位置，包括五个子步骤：

1)提取星图I_s中高对比度的局部图像I_f；

2)计算拖尾星点旋转角；

本发明用Random变换计算拖尾星点的旋转角。Radon变换将二维图像转化到直线参数空间。星图中每条直线在直线参数空间上都有一个峰值。这个特性可以用于图像处理中的直线检测。在二维欧几里德空间中的基本Radon变换为：

其中I为二维图像灰度分布函数，θ为投影角，δ为冲击函数，ρ为(x,y)平面内直线到原点的距离。若进行基本Radon变换，当特征图像投影在拖尾方向时，能量最集中，如附图3(a)所示。沿着别的方向能量较分散，如附图3(b)所示。如果投影方向垂直于拖尾方向，投影能量最分散，如附图3(c)所示。

但是，由于图像的形状为矩形，在相同θ条件下积分路径随着ρ的不同而不同，即，投影值不仅仅与灰度值相关，而且与积分路径相关。因此，基本Radon变换是在不同积分路径下比较投影值。结果不能准确反映图像上的直线信息。本发明提出一种归一化的Radon变换算法解决这个问题。

如公式(3.2)所示，积分路径被归一化，于是，在Radon参数空间中的投影值可以在相同积分路径下进行比较；

在归一化Radon变换的基础上，定义目标函数f(θ)：

其中：

目标函数通过以下方式建立：(I)归一化投影曲线；(II)对该曲线进行积分。θ越接近拖尾旋转角时，函数值越小。

如附图4所示，阴影部分的面积为f(θ)函数的值。图中的三条曲线分别为r_N0(ρ,θ_a)，r_N0(ρ,θ₁)，r_N0(ρ,θ₂)归一化函数。很容易得知，当θ＝θ_a时，r(ρ,θ)曲线的能量最集中。因此，阴影部分的面积最小。数学证明如下：

如附图5所示，模糊星点的图像在θ＝θ_c和θ≠θ_c的条件下被分别绘制在(a)和(b)图上。然后，在正下方绘制出归一化投影曲线r(ρ,θ_c)和r(ρ,θ)。

当θ＝θ_c时：

当θ≠θ_c时，

最后，通过下式计算旋转角：

θ_c＝arg{min[f(θ)]} (3.7)

3)计算拖尾星点长度；

待解算拖尾长度的图像为拖尾方向为水平的特征图像，可以用标准星点模型进行分析。积分长度用L表示，在水平方向上对公式(2.4)进行积分并对其归一化：

将公式(2.5)、(2.6)代入公式(3.8)，得到：

分析沿x轴方向的积分，x轴方向的子函数为：

不同积分步长L条件下，g₁(x)有不同性质，因此定义函数：

提出一个以L为自变量的基于残差平方和的目标函数e(L)来估计f_v(L)以得到高精度的拖尾星点长度。当该函数值达到最小时，L值最接近实际拖尾星点长度：

其中，参数和通过下式计算：

其中，L∈[1,L₀]。推导过程如下：

设计分段函数f_p(L)以估计函数f_v(L)：

模型的参数使用最小二乘法来估计。

求解上式得到参数和

附图6给出一个示例。拖尾的实际长度为20.0，使用总残差函数计算出的长度为21.0。该结果在误差范围之内。

4)矢量累积算法

整幅星图通过拖尾星点长度沿着旋转角方向累积，得到结果图像I_r(x,y)。在解算的旋转角θ_c和长度L_c基础上，将原始图像上沿着图像中心旋转θ_c，并以长度L_c进行累积，表达式为：

I_r(x,y)＝∫_ζI_s(ξ,η)ds (3.17)

其中，ζ为一条直线，起点为(x,y)，终点为(x+cosθ_cL_c,y+sinθ_cL_c)。经过矢量累积后，获得结果图像I_r(x,y)。此图像的特点为：该图像上的拖尾星点经过能量累积，成为一个单峰函数，峰值为拖尾恒星的起始点。附图7为矢量累积的示意图。左图圆圈中心为拖尾星点的端点，右图圆圈中心为对原始图像进行累积后的峰值位置。这两个圆圈在图像上的位置相同。

5)对原始图像进行方向累积后的结果图像上，使用质心法计算星点坐标；

其中(x,y)为图像坐标点，u_xy为此点的强度值，Ω为星点所在的区域。由矢量累积方法可知，理论上，质心位置为拖尾星点的起始点位置。

最后，提出一个数学模型估计拖尾星点位置的精度，模型示意图如附图8，最优旋转角和解算的旋转角方向存在偏差，如图中左下方所示。圆圈A的位置为实际质心位置。圆圈B的位置通过存在偏差的旋转角和存在偏差的拖尾星点长度解算。A和B之间的偏差为星点位置估计偏差，它由旋转角和拖尾长度的解算偏差造成。δ_x表示x轴方向的偏差，δ_y为y轴方向的偏差。分析解算的旋转角偏差和拖尾星点长度偏差带来的拖尾星点位置偏差误差δ_x和δ_y：

其中(x_A,y_A)为实际星点质心坐标；

下面从理论上推导此偏差。在解算得到拖尾旋转角和星点长度的条件下，星点的质心可以通过下式获得，可得公式：

其中，x₂＝x₁+cos△θ_c△L_c,y₂＝y₁+sin△θ_c△L_c,△θ_c＝θ_c-θ_r和△L_c＝L_c-L_r。θ_c和θ_r分别为拖尾星点的解算旋转角和实际旋转角。L_c和L_r分别为拖尾星点的长度解算值和实际值；

I_r(x,y)的解析表达式为：

当x₁<-x₀且x₂<-x₀,

I_r(x,y)＝G₁(x₁,x₂) (3.21)

当x₁<-x₀且–x₀≤x₂<x₀,

I_r(x,y)＝G₁(x₁,-x₀)+G₂(-x₀,x₂) (3.22)

当x₁<-x₀且x₂≥x₀,

I_r(x,y)＝G₁(x₁,-x₀)+G₂(-x₀,x₀)+G₃(x₀,x₂) (3.23)

当–x₀≤x₁<x₀且–x₀≤x₂<x₀,

I_r(x,y)＝G₂(x₁,x₂) (3.24)

当–x₀≤x₁<x₀且x₂≥x₀,

I_r(x,y)＝G₂(x₁,x₀)+G₃(x₀,x₂) (3.25)

当x₁≥x₀且x₂≥x₀,

I_r(x,y)＝G₃(x₁,x₂) (3.26)

为简化计算，分别在三个区间[-∞,-x₀]，[-x₀,x₀]和[x₀,∞]上定义三个函数G₁(x₁,x₂)、G₂(x₁,x₂)、G₃(x₁,x₂)。下面推导这三个函数：

当△θ_c为零时,公式(3.20)可以被简化为：

当x₁<-x₀且x₂<-x₀时，得：

当-x₀≤x₁<x₀且–x₀≤x₂<x₀时,

当x₁≥x₀且x₂≥x₀,

当△θ_c不为零时，可得：

当x₁<-x₀且x₂<-x₀,

当–x₀≤x₁<x₀且–x₀≤x₂<x₀,

当x₁≥x₀且x₂≥x₀,

在前文中证明了在星图中所有的条痕都可以使用相同的模型建模。然后分析了图像上所有星点的精度。在不同的星点旋转角和长度条件下进行仿真实验以验证I_r(x,y)的解析表达式。由于篇幅所限，只列出其中一颗星的精度结果作为示例，如附图9和附图10所示。需要注意的是，所有恒星的精度分析结果都是相似的；

仿真拖尾星点的基本设计参数如下：I₀为50.0，x₀为10.0，σ为5.0。附图9和附图10显示了δ_x和δ_y随着△L_c和△θ_c的变化。│δ_y│随着△θ_c增加而增加。同时│δ_y│随着△L_c增加而增加。当△L_c小于等于零时，│δ_x│随着△θ_c值的增加而增加。但是当△L_c大于零时，│δ_x│随着△θ_c的增加先减少后增加。根据公式(3.20)至(3.24)解算出理论值用于对比。从附图9和附图10这两幅图像中可知，理论值与实验值符合的很好。

本发明可以达到以下的技术效果

1、适用于天文导航系统深空恒星星点定位，能够有效提高对拖尾星点定位精度；

2、可以增强拖尾星点的信噪比和极限星等；

3、在飞行器快速转动条件下可有效工作；

4、在星敏感器非线性转动的条件下可稳定工作。

附图说明

图1基于矢量累积方法的星点位置估计流程图；

图2线性拖尾星点模型；

图3特征图像在不同方向的投影示意图；

图4 f(θ)函数示意图；

图5 f(θ)函数性质的数学证明

图6总残差值e(L)曲线示例图；

图7矢量累积示意图；

图8拖尾星点位置估计精度的模型示意图；

图9以△Lc和△θc为变量的函数δx；

图10以△Lc和△θc为变量的函数δy。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细说明：

本发明提供一种深空背景下的恒星星点拖尾目标检测方法，如附图1所示，该方法通过以下步骤实现：

1)从姿态方程出发分析拖尾星点的形状，提出一种分析拖尾星点的线性模型，建立拖尾星点数学模型；

2)实行矢量累积算法估计星点位置，包括五个子步骤：提取星图中高对比度的局部图像；通过归一化Radon变换计算实际拖尾模型和标准拖尾模型之间的旋转角；通过这个旋转角旋转局部图像，计算拖尾星点的长度；整幅星图通过拖尾星点长度沿着旋转角方向进行累积得到结果图像；使用质心法处理结果图像求解星点坐标；

3)提出一个数学模型估计拖尾星点位置的精度，分析解算的旋转角偏差和拖尾星点长度偏差带来的拖尾星点位置偏差误差。

Claims

1.基于矢量累积的恒星拖尾星点提取方法，针对星敏感器转动条件下恒星目标信噪比下降导致的定姿精度下降问题，以及拖尾恒星加速转动条件下提点精度误差增加的问题，其特征在于，首先，从原始图像中提取特征图像；然后分别计算拖尾星点的旋转角和长度；之后，通过矢量累积算法，将整幅星图沿着拖尾建造矢量以拖尾步长进行累积，提高拖尾星点的信噪比，得到能量增强的目标，获得结果图像；最后对恒星目标通过传统质心法对目标进行检测并进行位置估计。

2.根据权利要求1所述的基于矢量累积的恒星拖尾星点提取方法，其特征在于，所述为计算拖尾星点的旋转角和长度，建立基于线性拖尾模型的星点数学模型：

从姿态方程W＝AV出发分析拖尾星点的形状，提出一种分析拖尾星点的线性模型，建立拖尾星点的数学模型，其中A是星敏感器输出的飞行器姿态矩阵，W和V分别表示体坐标系和惯性坐标系下的单位矢量，W和V的值可由以下方程算得：

<mrow> <mi>V</mi> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mi>c</mi> <mi>o</mi> <mi>s</mi> <mi>&alpha;</mi> <mi>c</mi> <mi>o</mi> <mi>s</mi> <mi>&delta;</mi> </mtd> </mtr> <mtr> <mtd> <mi>s</mi> <mi>i</mi> <mi>n</mi> <mi>&alpha;</mi> <mi>c</mi> <mi>o</mi> <mi>s</mi> <mi>&delta;</mi> </mtd> </mtr> <mtr> <mtd> <mi>sin</mi> <mi>&delta;</mi> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2.2</mn> <mo>)</mo> </mrow> </mrow>

设(u,v)表示星点在焦平面上的速度，t为时间，而相机的焦距与时间无关，因此，在大多数情况下，(x²+y²)/f小于0.01，于是，以时间为自变量对姿态方程求导可得近似方程：

在视场较小和飞行器姿态变化较小条件下，星点在焦平面上每个位置运动速度认为近似相等，即拖尾星点形状在焦平面上每个位置都是近似相同的，因此，提出一种基于线性拖尾模型的星点数学模型如下：

<mrow> <mi>I</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mfrac> <msub> <mi>I</mi> <mn>0</mn> </msub> <mrow> <mn>2</mn> <msup> <mi>&pi;&sigma;</mi> <mn>2</mn> </msup> </mrow> </mfrac> <mi>exp</mi> <mrow> <mo>(</mo> <mfrac> <mrow> <mo>-</mo> <msup> <mrow> <mo>(</mo> <mi>x</mi> <mo>-</mo> <msub> <mi>x</mi> <mn>0</mn> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> <mrow> <mn>2</mn> <msup> <mi>&sigma;</mi> <mn>2</mn> </msup> </mrow> </mfrac> <mo>)</mo> </mrow> <mi>exp</mi> <mrow> <mo>(</mo> <mfrac> <mrow> <mo>-</mo> <msup> <mi>y</mi> <mn>2</mn> </msup> </mrow> <mrow> <mn>2</mn> <msup> <mi>&sigma;</mi> <mn>2</mn> </msup> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>,</mo> <mi>x</mi> <mo>&GreaterEqual;</mo> <msub> <mi>x</mi> <mn>0</mn> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mfrac> <msub> <mi>I</mi> <mn>0</mn> </msub> <mrow> <mn>2</mn> <msup> <mi>&pi;&sigma;</mi> <mn>2</mn> </msup> </mrow> </mfrac> <mi>exp</mi> <mrow> <mo>(</mo> <mfrac> <mrow> <mo>-</mo> <msup> <mi>y</mi> <mn>2</mn> </msup> </mrow> <mrow> <mn>2</mn> <msup> <mi>&sigma;</mi> <mn>2</mn> </msup> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>,</mo> <mo>-</mo> <msub> <mi>x</mi> <mn>0</mn> </msub> <mo>&le;</mo> <mi>x</mi> <mo>&le;</mo> <msub> <mi>x</mi> <mn>0</mn> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mfrac> <msub> <mi>I</mi> <mn>0</mn> </msub> <mrow> <mn>2</mn> <msup> <mi>&pi;&sigma;</mi> <mn>2</mn> </msup> </mrow> </mfrac> <mi>exp</mi> <mrow> <mo>(</mo> <mfrac> <mrow> <mo>-</mo> <msup> <mrow> <mo>(</mo> <mi>x</mi> <mo>+</mo> <msub> <mi>x</mi> <mn>0</mn> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> <mrow> <mn>2</mn> <msup> <mi>&sigma;</mi> <mn>2</mn> </msup> </mrow> </mfrac> <mo>)</mo> </mrow> <mi>exp</mi> <mrow> <mo>(</mo> <mfrac> <mrow> <mo>-</mo> <msup> <mi>y</mi> <mn>2</mn> </msup> </mrow> <mrow> <mn>2</mn> <msup> <mi>&sigma;</mi> <mn>2</mn> </msup> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>,</mo> <mi>x</mi> <mo><</mo> <mo>-</mo> <msub> <mi>x</mi> <mn>0</mn> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2.4</mn> <mo>)</mo> </mrow> </mrow>

其中I₀表示星敏感器上的星点总能量，(-x₀,0)为曝光开始时焦平面上星点位置，(x₀,0)为曝光结束时焦平面上星点位置，σ为星点的高斯半径，拖尾星点的中心位于坐标原点，在曝光时间沿着x轴方向移动，这个模型为标准星点数学模型。

3.根据权利要求1所述的基于矢量累积的恒星拖尾星点提取方法，其特征在于，所述矢量累积算法，包括五个步骤：

1)提取星图I_s中高对比度的局部图像I_f；

2)计算拖尾星点旋转角；

用Random变换计算拖尾星点的旋转角，Radon变换将二维图像转化到直线参数空间，星图中每条直线在直线参数空间上都有一个峰值，这个特性用于图像处理中的直线检测，在二维欧几里德空间中的基本Radon变换为：

<mrow> <mi>r</mi> <mrow> <mo>(</mo> <mi>&rho;</mi> <mo>,</mo> <mi>&theta;</mi> <mo>)</mo> </mrow> <mo>=</mo> <msubsup> <mo>&Integral;</mo> <mrow> <mo>-</mo> <mi>&infin;</mi> </mrow> <mi>&infin;</mi> </msubsup> <msubsup> <mo>&Integral;</mo> <mrow> <mo>-</mo> <mi>&infin;</mi> </mrow> <mi>&infin;</mi> </msubsup> <mi>I</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mi>&delta;</mi> <mrow> <mo>(</mo> <mi>x</mi> <mi> </mi> <mi>c</mi> <mi>o</mi> <mi>s</mi> <mi>&theta;</mi> <mo>+</mo> <mi>y</mi> <mi> </mi> <mi>s</mi> <mi>i</mi> <mi>n</mi> <mi>&theta;</mi> <mo>-</mo> <mi>&rho;</mi> <mo>)</mo> </mrow> <mi>d</mi> <mi>x</mi> <mi>d</mi> <mi>y</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3.1</mn> <mo>)</mo> </mrow> </mrow>

其中I为二维图像灰度分布函数，θ为投影角，δ为冲击函数，若进行基本Radon变换，当特征图像投影在拖尾方向时，能量最集中，沿着别的方向能量较分散，如果投影方向垂直于拖尾方向，投影能量最分散，

但是，由于图像的形状为矩形，在相同θ条件下积分路径随着ρ的不同而不同，即，投影值不仅仅与灰度值相关，而且与积分路径相关，基本Radon变换结果不能准确反映图像上的直线信息，提出一种归一化的Radon变换算法解决这个问题，如公式(3.2)所示，积分路径被归一化，于是，在Radon参数空间中的投影值在相同积分路径下进行比较；

<mrow> <msub> <mi>r</mi> <mrow> <mi>N</mi> <mn>0</mn> </mrow> </msub> <mrow> <mo>(</mo> <mi>&rho;</mi> <mo>,</mo> <mi>&theta;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <msubsup> <mo>&Integral;</mo> <mn>0</mn> <mi>H</mi> </msubsup> <msubsup> <mo>&Integral;</mo> <mn>0</mn> <mi>W</mi> </msubsup> <mi>I</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mi>&delta;</mi> <mrow> <mo>(</mo> <mi>x</mi> <mi> </mi> <mi>c</mi> <mi>o</mi> <mi>s</mi> <mi>&theta;</mi> <mo>+</mo> <mi>y</mi> <mi> </mi> <mi>s</mi> <mi>i</mi> <mi>n</mi> <mi>&theta;</mi> <mo>-</mo> <mi>&rho;</mi> <mo>)</mo> </mrow> <mi>d</mi> <mi>x</mi> <mi>d</mi> <mi>y</mi> </mrow> <mrow> <msubsup> <mo>&Integral;</mo> <mn>0</mn> <mi>H</mi> </msubsup> <msubsup> <mo>&Integral;</mo> <mn>0</mn> <mi>W</mi> </msubsup> <mi>&delta;</mi> <mrow> <mo>(</mo> <mi>x</mi> <mi> </mi> <mi>cos</mi> <mi>&theta;</mi> <mo>+</mo> <mi>y</mi> <mi> </mi> <mi>s</mi> <mi>i</mi> <mi>n</mi> <mi>&theta;</mi> <mo>-</mo> <mi>&rho;</mi> <mo>)</mo> </mrow> <mi>d</mi> <mi>x</mi> <mi>d</mi> <mi>y</mi> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3.2</mn> <mo>)</mo> </mrow> </mrow>

在归一化Radon变换的基础上，定义目标函数f(θ)：

<mrow> <mi>f</mi> <mrow> <mo>(</mo> <mi>&theta;</mi> <mo>)</mo> </mrow> <mo>=</mo> <msubsup> <mo>&Integral;</mo> <mrow> <mo>-</mo> <mi>&infin;</mi> </mrow> <mi>&infin;</mi> </msubsup> <msub> <mi>r</mi> <mrow> <mi>N</mi> <mn>0</mn> </mrow> </msub> <mrow> <mo>(</mo> <mi>&rho;</mi> <mo>,</mo> <mi>&theta;</mi> <mo>)</mo> </mrow> <mi>d</mi> <mi>&rho;</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3.3</mn> <mo>)</mo> </mrow> </mrow>

目标函数通过以下方式建立：(I)归一化投影曲线；(II)对该曲线进行积分，θ越接近拖尾旋转角时，函数值越小，通过归一化Radon变换计算实际拖尾模型和标准拖尾模型之间的旋转角，定义解算出的旋转角为θ_c，实际旋转角为θ_a，其中：

<mrow> <msub> <mi>r</mi> <mrow> <mi>N</mi> <mn>0</mn> </mrow> </msub> <mrow> <mo>(</mo> <mi>&rho;</mi> <mo>,</mo> <mi>&theta;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <msub> <mi>r</mi> <mi>N</mi> </msub> <mrow> <mo>(</mo> <mi>&rho;</mi> <mo>,</mo> <mi>&theta;</mi> <mo>)</mo> </mrow> <mo>-</mo> <munder> <mi>min</mi> <mi>&theta;</mi> </munder> <mo>&lsqb;</mo> <msub> <mi>r</mi> <mi>N</mi> </msub> <mrow> <mo>(</mo> <mi>&rho;</mi> <mo>,</mo> <mi>&theta;</mi> <mo>)</mo> </mrow> <mo>&rsqb;</mo> </mrow> <mrow> <munder> <mrow> <mi>m</mi> <mi>a</mi> <mi>x</mi> </mrow> <mi>&theta;</mi> </munder> <mo>&lsqb;</mo> <msub> <mi>r</mi> <mi>N</mi> </msub> <mrow> <mo>(</mo> <mi>&rho;</mi> <mo>,</mo> <mi>&theta;</mi> <mo>)</mo> </mrow> <mo>&rsqb;</mo> <mo>-</mo> <munder> <mi>min</mi> <mi>&theta;</mi> </munder> <mo>&lsqb;</mo> <msub> <mi>r</mi> <mi>N</mi> </msub> <mrow> <mo>(</mo> <mi>&rho;</mi> <mo>,</mo> <mi>&theta;</mi> <mo>)</mo> </mrow> <mo>&rsqb;</mo> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3.4</mn> <mo>)</mo> </mrow> </mrow>

当θ＝θ_c时：

<mrow> <mtable> <mtr> <mtd> <mrow> <mi>f</mi> <mrow> <mo>(</mo> <msub> <mi>&theta;</mi> <mi>c</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <msubsup> <mo>&Integral;</mo> <mrow> <mo>-</mo> <mi>&infin;</mi> </mrow> <mi>&infin;</mi> </msubsup> <mfrac> <mrow> <msub> <mi>r</mi> <mi>N</mi> </msub> <mrow> <mo>(</mo> <mi>&rho;</mi> <mo>,</mo> <msub> <mi>&theta;</mi> <mi>c</mi> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <msubsup> <mo>&Integral;</mo> <mrow> <mo>-</mo> <mi>&infin;</mi> </mrow> <mrow> <mo>-</mo> <msub> <mi>x</mi> <mn>0</mn> </msub> </mrow> </msubsup> <mfrac> <msub> <mi>I</mi> <mn>0</mn> </msub> <mrow> <mn>2</mn> <msup> <mi>&pi;&sigma;</mi> <mn>2</mn> </msup> </mrow> </mfrac> <mi>exp</mi> <mrow> <mo>(</mo> <mfrac> <mrow> <mo>-</mo> <msup> <mrow> <mo>(</mo> <mi>x</mi> <mo>-</mo> <msub> <mi>x</mi> <mn>0</mn> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> <mrow> <mn>2</mn> <msup> <mi>&sigma;</mi> <mn>2</mn> </msup> </mrow> </mfrac> <mo>)</mo> </mrow> <mi>d</mi> <mi>x</mi> <mo>+</mo> <mfrac> <msub> <mi>I</mi> <mn>0</mn> </msub> <mrow> <mn>2</mn> <msup> <mi>&pi;&sigma;</mi> <mn>2</mn> </msup> </mrow> </mfrac> <mn>2</mn> <msub> <mi>x</mi> <mn>0</mn> </msub> <mo>+</mo> <msubsup> <mo>&Integral;</mo> <msub> <mi>x</mi> <mn>0</mn> </msub> <mi>&infin;</mi> </msubsup> <mfrac> <msub> <mi>I</mi> <mn>0</mn> </msub> <mrow> <mn>2</mn> <msup> <mi>&pi;&sigma;</mi> <mn>2</mn> </msup> </mrow> </mfrac> <mi>exp</mi> <mrow> <mo>(</mo> <mfrac> <mrow> <mo>-</mo> <msup> <mrow> <mo>(</mo> <mi>x</mi> <mo>-</mo> <msub> <mi>x</mi> <mn>0</mn> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> <mrow> <mn>2</mn> <msup> <mi>&sigma;</mi> <mn>2</mn> </msup> </mrow> </mfrac> <mo>)</mo> </mrow> <mi>d</mi> <mi>x</mi> </mrow> </mfrac> <mi>d</mi> <mi>&rho;</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>=</mo> <msubsup> <mo>&Integral;</mo> <mrow> <mo>-</mo> <mi>&infin;</mi> </mrow> <mi>&infin;</mi> </msubsup> <mi>exp</mi> <mrow> <mo>(</mo> <mfrac> <mrow> <mo>-</mo> <msup> <mi>y</mi> <mn>2</mn> </msup> </mrow> <mrow> <mn>2</mn> <msup> <mi>&sigma;</mi> <mn>2</mn> </msup> </mrow> </mfrac> <mo>)</mo> </mrow> <mi>d</mi> <mi>x</mi> <mo>=</mo> <msqrt> <mrow> <mn>2</mn> <mi>&pi;</mi> </mrow> </msqrt> <mi>&sigma;</mi> </mrow> </mtd> </mtr> </mtable> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3.5</mn> <mo>)</mo> </mrow> </mrow>

当θ≠θ_c时，

<mrow> <mtable> <mtr> <mtd> <mrow> <mi>f</mi> <mrow> <mo>(</mo> <mi>&theta;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mn>2</mn> <msub> <mi>x</mi> <mn>0</mn> </msub> <mi>sin</mi> <mrow> <mo>(</mo> <mo>|</mo> <mrow> <mi>&theta;</mi> <mo>-</mo> <msub> <mi>&theta;</mi> <mi>c</mi> </msub> </mrow> <mo>|</mo> <mo>)</mo> </mrow> <mo>+</mo> <mfrac> <mrow> <msubsup> <mo>&Integral;</mo> <mrow> <mo>-</mo> <mi>&infin;</mi> </mrow> <mn>0</mn> </msubsup> <msubsup> <mo>&Integral;</mo> <mrow> <mo>-</mo> <mi>&infin;</mi> </mrow> <mi>&infin;</mi> </msubsup> <msub> <mi>I</mi> <mn>0</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mi>d</mi> <mi>y</mi> <mi>d</mi> <mi>x</mi> <mo>+</mo> <msubsup> <mo>&Integral;</mo> <mrow> <mn>2</mn> <msub> <mi>x</mi> <mn>0</mn> </msub> </mrow> <mi>&infin;</mi> </msubsup> <msubsup> <mo>&Integral;</mo> <mrow> <mo>-</mo> <mi>&infin;</mi> </mrow> <mi>&infin;</mi> </msubsup> <msub> <mi>I</mi> <mrow> <mi>r</mi> <mn>0</mn> </mrow> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mi>d</mi> <mi>y</mi> <mi>d</mi> <mi>x</mi> </mrow> <mrow> <msubsup> <mo>&Integral;</mo> <mrow> <mo>-</mo> <mi>&infin;</mi> </mrow> <mi>&infin;</mi> </msubsup> <msub> <mi>I</mi> <mrow> <mi>r</mi> <mn>0</mn> </mrow> </msub> <mrow> <mo>(</mo> <mn>0</mn> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mi>d</mi> <mi>y</mi> </mrow> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>&ap;</mo> <mn>2</mn> <msub> <mi>x</mi> <mn>0</mn> </msub> <mi>sin</mi> <mrow> <mo>(</mo> <mo>|</mo> <mrow> <mi>&theta;</mi> <mo>-</mo> <msub> <mi>&theta;</mi> <mi>c</mi> </msub> </mrow> <mo>|</mo> <mo>)</mo> </mrow> <mo>+</mo> <msqrt> <mrow> <mn>2</mn> <mi>&pi;</mi> </mrow> </msqrt> <mi>&sigma;</mi> <mo>></mo> <msqrt> <mrow> <mn>2</mn> <mi>&pi;</mi> </mrow> </msqrt> <mi>&sigma;</mi> <mo>=</mo> <mi>f</mi> <mrow> <mo>(</mo> <msub> <mi>&theta;</mi> <mi>c</mi> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3.6</mn> <mo>)</mo> </mrow> </mrow>

通过下式计算旋转角：

θ_c＝arg{min[f(θ)]} (3.7)，

3)计算拖尾星点长度；

待解算拖尾长度的图像为拖尾方向为水平的特征图像，用标准星点模型进行分析，积分长度用L表示，在水平方向上对星点数学模型公式(2.4)进行积分并对其归一化：

<mrow> <mi>g</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>L</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <mn>2</mn> <msup> <mi>&pi;&sigma;</mi> <mn>2</mn> </msup> </mrow> <mrow> <msub> <mi>I</mi> <mn>0</mn> </msub> <mi>L</mi> </mrow> </mfrac> <msubsup> <mo>&Integral;</mo> <mi>x</mi> <mrow> <mi>x</mi> <mo>+</mo> <mi>L</mi> </mrow> </msubsup> <mi>I</mi> <mrow> <mo>(</mo> <mi>&xi;</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mi>d</mi> <mi>&xi;</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3.8</mn> <mo>)</mo> </mrow> </mrow>

将公式(2.5)、(2.6)代入公式(3.8)，得到：

<mrow> <mi>g</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>L</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <mn>2</mn> <msup> <mi>&pi;&sigma;</mi> <mn>2</mn> </msup> </mrow> <mrow> <msub> <mi>I</mi> <mn>0</mn> </msub> <mi>L</mi> </mrow> </mfrac> <msub> <mi>I</mi> <mn>2</mn> </msub> <mrow> <mo>(</mo> <mi>y</mi> <mo>)</mo> </mrow> <msubsup> <mo>&Integral;</mo> <mi>x</mi> <mrow> <mi>x</mi> <mo>+</mo> <mi>L</mi> </mrow> </msubsup> <msub> <mi>I</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>&xi;</mi> <mo>)</mo> </mrow> <mi>d</mi> <mi>&xi;</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3.9</mn> <mo>)</mo> </mrow> </mrow>

分析沿x轴方向的积分，x轴方向的子函数为：

<mrow> <msub> <mi>g</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>L</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <mn>2</mn> <msup> <mi>&pi;&sigma;</mi> <mn>2</mn> </msup> </mrow> <mrow> <msub> <mi>I</mi> <mn>0</mn> </msub> <mi>L</mi> </mrow> </mfrac> <msubsup> <mo>&Integral;</mo> <mi>x</mi> <mrow> <mi>x</mi> <mo>+</mo> <mi>L</mi> </mrow> </msubsup> <msub> <mi>I</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>&xi;</mi> <mo>)</mo> </mrow> <mi>d</mi> <mi>&xi;</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3.10</mn> <mo>)</mo> </mrow> </mrow>

不同L下条件下，g₁(x)有不同性质，因此定义函数：

提出一个以L为自变量的基于残差平方和的目标函数e(L)来估计f_v(L)以得到高精度的拖尾星点长度，当该函数值达到最小时，L值最接近实际拖尾星点长度：

<mrow> <mi>e</mi> <mrow> <mo>(</mo> <mi>L</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>L</mi> </munderover> <msup> <mrow> <mo>&lsqb;</mo> <mover> <msub> <mi>a</mi> <mn>0</mn> </msub> <mo>^</mo> </mover> <mo>-</mo> <msub> <mi>f</mi> <mi>v</mi> </msub> <mrow> <mo>(</mo> <mi>l</mi> <mo>)</mo> </mrow> <mo>&rsqb;</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>l</mi> <mo>=</mo> <mi>L</mi> </mrow> <msub> <mi>L</mi> <mn>0</mn> </msub> </munderover> <msup> <mrow> <mo>&lsqb;</mo> <mover> <msub> <mi>a</mi> <mn>1</mn> </msub> <mo>^</mo> </mover> <msup> <mrow> <mo>(</mo> <mi>l</mi> <mo>-</mo> <mi>L</mi> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <mover> <msub> <mi>a</mi> <mn>0</mn> </msub> <mo>^</mo> </mover> <mo>-</mo> <msub> <mi>f</mi> <mi>v</mi> </msub> <mrow> <mo>(</mo> <mi>l</mi> <mo>)</mo> </mrow> <mo>&rsqb;</mo> </mrow> <mn>2</mn> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3.12</mn> <mo>)</mo> </mrow> </mrow>

其中，参数和通过下式计算：

<mrow> <mover> <msub> <mi>a</mi> <mn>0</mn> </msub> <mo>^</mo> </mover> <mo>=</mo> <mfrac> <mn>1</mn> <mi>L</mi> </mfrac> <munderover> <mi>&Sigma;</mi> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>L</mi> </munderover> <msub> <mi>f</mi> <mi>v</mi> </msub> <mrow> <mo>(</mo> <mi>l</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3.13</mn> <mo>)</mo> </mrow> </mrow>

<mrow> <mover> <msub> <mi>a</mi> <mn>1</mn> </msub> <mo>^</mo> </mover> <mo>=</mo> <mfrac> <mrow> <munderover> <mi>&Sigma;</mi> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>L</mi> <mn>0</mn> </msub> </munderover> <msub> <mi>f</mi> <mi>v</mi> </msub> <mrow> <mo>(</mo> <mi>l</mi> <mo>)</mo> </mrow> <mo>-</mo> <mover> <msub> <mi>a</mi> <mn>0</mn> </msub> <mo>^</mo> </mover> <mrow> <mo>(</mo> <msub> <mi>L</mi> <mn>0</mn> </msub> <mo>-</mo> <mi>L</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> <mrow> <munderover> <mi>&Sigma;</mi> <mrow> <mi>l</mi> <mo>=</mo> <mi>L</mi> </mrow> <msub> <mi>L</mi> <mn>0</mn> </msub> </munderover> <msup> <mrow> <mo>(</mo> <mi>l</mi> <mo>-</mo> <mi>L</mi> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3.14</mn> <mo>)</mo> </mrow> </mrow>

其中，L∈[1,L₀]，推导过程：设计分段函数f_p(L)以估计函数f_v(L)：

模型的参数使用最小二乘法来估计，

<mrow> <mtable> <mtr> <mtd> <mrow> <munderover> <mi>&Sigma;</mi> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>L</mi> </munderover> <msup> <mrow> <mo>&lsqb;</mo> <msub> <mi>a</mi> <mn>0</mn> </msub> <mo>-</mo> <msub> <mi>f</mi> <mi>v</mi> </msub> <mrow> <mo>(</mo> <mi>l</mi> <mo>)</mo> </mrow> <mo>&rsqb;</mo> </mrow> <mn>2</mn> </msup> <mo>&RightArrow;</mo> <mi>min</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <munderover> <mi>&Sigma;</mi> <mrow> <mi>l</mi> <mo>=</mo> <mi>L</mi> </mrow> <msub> <mi>L</mi> <mn>0</mn> </msub> </munderover> <msup> <mrow> <mo>&lsqb;</mo> <msub> <mi>a</mi> <mn>1</mn> </msub> <msup> <mrow> <mo>(</mo> <mi>l</mi> <mo>-</mo> <mi>L</mi> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>-</mo> <msub> <mi>f</mi> <mi>v</mi> </msub> <mrow> <mo>(</mo> <mi>l</mi> <mo>)</mo> </mrow> <mo>&rsqb;</mo> </mrow> <mn>2</mn> </msup> <mo>&RightArrow;</mo> <mi>min</mi> </mrow> </mtd> </mtr> </mtable> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3.16</mn> <mo>)</mo> </mrow> </mrow>

4)矢量累积；

整幅星图通过拖尾星点长度沿着旋转角方向累积，得到结果图像I_r，在解算的旋转角θ_c和长度L_c基础上，将原始图像上沿着图像中心旋转θ_c，并以长度L_c进行累积，表达式为：

I_r(x,y)＝∫_ζI_s(ξ,η)ds (3.17)

其中，ξ为一条直线，起点为(x,y)，终点为(x+cosθ_cL_c,y+sinθ_cL_c)，经过矢量累积后，获得结果图像I_r(x,y)，

5)使用质心法处理结果图像I_r求解星点坐标；

<mrow> <mover> <mi>x</mi> <mo>^</mo> </mover> <mo>=</mo> <mfrac> <mrow> <munder> <mi>&Sigma;</mi> <mi>&Omega;</mi> </munder> <msub> <mi>xu</mi> <mrow> <mi>x</mi> <mi>y</mi> </mrow> </msub> </mrow> <mrow> <munder> <mi>&Sigma;</mi> <mi>&Omega;</mi> </munder> <msub> <mi>u</mi> <mrow> <mi>x</mi> <mi>y</mi> </mrow> </msub> </mrow> </mfrac> <mo>,</mo> <mover> <mi>y</mi> <mo>^</mo> </mover> <mo>=</mo> <mfrac> <mrow> <munder> <mi>&Sigma;</mi> <mi>&Omega;</mi> </munder> <msub> <mi>yu</mi> <mrow> <mi>x</mi> <mi>y</mi> </mrow> </msub> </mrow> <mrow> <munder> <mi>&Sigma;</mi> <mi>&Omega;</mi> </munder> <msub> <mi>u</mi> <mrow> <mi>x</mi> <mi>y</mi> </mrow> </msub> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3.18</mn> <mo>)</mo> </mrow> </mrow>

其中(x,y)为图像坐标点，u_xy为此点的强度值，Ω为星点所在的区域，由矢量累积方法可知，理论上，质心位置为拖尾星点的起始点位置。