CN115659196A - 基于非线性偏差演化的天基光学观测短弧关联与聚类方法 - Google Patents

基于非线性偏差演化的天基光学观测短弧关联与聚类方法 Download PDF

Info

Publication number
CN115659196A
CN115659196A CN202211594600.6A CN202211594600A CN115659196A CN 115659196 A CN115659196 A CN 115659196A CN 202211594600 A CN202211594600 A CN 202211594600A CN 115659196 A CN115659196 A CN 115659196A
Authority
CN
China
Prior art keywords
observation
arc
space
matrix
clustering
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202211594600.6A
Other languages
English (en)
Other versions
CN115659196B (zh
Inventor
杨震
李嘉胜
罗亚中
王�华
朱阅訸
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
National University of Defense Technology
Original Assignee
National University of Defense Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by National University of Defense Technology filed Critical National University of Defense Technology
Priority to CN202211594600.6A priority Critical patent/CN115659196B/zh
Publication of CN115659196A publication Critical patent/CN115659196A/zh
Application granted granted Critical
Publication of CN115659196B publication Critical patent/CN115659196B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Radar Systems Or Details Thereof (AREA)
  • Image Analysis (AREA)

Abstract

本发明公开了一种基于非线性偏差演化的天基光学观测短弧关联与聚类方法，包括：首先获取天基短弧光学观测数据，进行数据预处理；根据观测弧段特征信息与先验信息，划定各观测弧段对应容许域；在容许域内对两观测弧段间最小马氏距离进行优化；根据两观测弧段间最小马氏距离判别两者是否关联；根据观测弧段两两关联结果，构建观测弧段关联矩阵；利用BEA算法将观测弧段关联矩阵变换为观测弧段聚类矩阵；根据观测弧段聚类矩阵元素排列特征进行分割，得到最终关联聚类结果。本发明应用于太空态势感知领域，解决现有技术中难以属于同一空间目标的观测弧段进行关联并聚类的问题，同时兼顾算法的计算正确率与计算效率。

Description

基于非线性偏差演化的天基光学观测短弧关联与聚类方法
技术领域
本发明涉及太空态势感知技术领域，具体是一种基于非线性偏差演化的天基光学观测短弧关联与聚类方法。
背景技术
随着航天事业的不断发展，在轨空间目标的数量正急剧增加，例如正在开展的“星链”计划，预计完成后将部署超过4万颗卫星。截止至2022年10月13日，可被空间目标监视网(U.S. Space Surveillance Network, SSN)追踪的直径大于10cm的在轨空间目标总数已达26174个，其中有效载荷数量为9719个，仅占总数的37.13%。据估计，有超过30万个直径大于1cm的空间目标在轨运行，直径小于1cm的更是数以百万计。空间目标的观测与编目是空间态势监视和碰撞预警的重要基础，对于维护在轨资产安全与空间安全具有重要意义，数量如此庞大的空间目标为空间目标观测编目工作的准确性提出了更高的要求。
在空间目标的观测工作中，天基光学观测由于其观测精度高，抗干扰性强等独特优势正受到越来越多的青睐，但由于天基观测卫星与被观测目标相对速度通常较大，导致单一观测弧段时长很短，通常不超过两分钟，被称为短弧观测片段。由于单个弧段时间长度较短，轨道确定的精度难以保证，很难直接进行目标编目，一般需要累积多个观测弧段。此时就存在观测弧段的关联和聚类问题，需要对属于同一空间目标的观测弧段进行识别。现有针对天基光学观测弧段的关联方法较少，且多数集中在观测弧段两两之间的关联匹配，关联准确性仍有进一步提升的空间，并且针对多观测弧段之间如何进行聚类的问题，当前领域研究仍有所欠缺。
发明内容
针对上述现有技术中的不足，本发明提供一种基于非线性偏差演化的天基光学观测短弧关联与聚类方法，解决现有技术中难以属于同一空间目标的观测弧段进行关联并聚类的问题，同时兼顾算法的计算正确率与计算效率。
为实现上述目的，本发明提供一种基于非线性偏差演化的天基光学观测短弧关联与聚类方法，包括如下步骤：
步骤1，利用天基光学观测卫星进行天基光学观测，获取多组分别属于不同空间目标的原始观测短弧片段，也称观测弧段，每个观测弧段数据包括若干个观测数据点，每个观测数据点由被观测目标相对于低轨光学观测卫星的赤经、赤纬、观测时刻以及观测平台的位置速度信息组成；
步骤2，分别对每个观测弧段中赤经、赤纬关于时间的函数式进行拟合，得到赤经、赤纬随时间的变化率信息，并对存在明显异常的观测数据点进行剔除；
步骤3，根据步骤2处理后得到的各观测弧段中赤经、赤纬随时间的变化率信息，结合天基观测卫星运行轨道与被观测目标大致运行轨道区间等先验信息，对各观测弧段对应目标轨道在斜距与斜距变化率平面上的容许域范围进行划定；
步骤4，对待关联观测弧段两两在步骤3中划定的容许域范围内对斜距与斜距变化率组合进行优化，结合运用航天器轨道预报与偏差演化算法，找到使赤经、赤纬预报值与实际观测值马氏距离最小的斜距与斜距变化率组合，并对两观测弧段间最小马氏距离进行记录；
步骤5，以步骤4中所记录的两观测弧段间最小马氏距离为关联判别依据，得到观测弧段两两关联匹配结果；
步骤6，根据步骤5中所得观测弧段两两关联匹配结果，构建观测弧段关联矩阵，并利用BEA算法对观测弧段关联矩阵进行行列变换，将观测弧段关联矩阵转换成观测弧段聚类矩阵；
步骤7，根据观测弧段聚类矩阵行列元素排列特征对观测弧段聚类矩阵进行分割，得到观测弧段关联聚类结果，实现对属于同一空间目标的观测弧段的关联聚类。
在其中一个实施例，步骤1的实现过程为：
已知通过光学观测卫星上安装的天基观测设备，进行天基光学观测后，得到的多组分别属于不同空间目标的天基测角数据
Figure 100002_DEST_PATH_IMAGE001
，即所述多组分别属于不同空间目标的观测弧段为
Figure 100002_DEST_PATH_IMAGE002
，
Figure 100002_DEST_PATH_IMAGE003
，
Figure 100002_DEST_PATH_IMAGE004
，其中，
Figure 100002_DEST_PATH_IMAGE005
为空间目标数量，
Figure 100002_DEST_PATH_IMAGE006
为观测弧段的数量；
Figure 100002_DEST_PATH_IMAGE007
为第
Figure 100002_DEST_PATH_IMAGE008
个空间目标的第
Figure 100002_DEST_PATH_IMAGE009
个观测弧段，
Figure 100002_DEST_PATH_IMAGE010
，其中，
Figure 100002_DEST_PATH_IMAGE011
为第
Figure 100002_DEST_PATH_IMAGE012
个空间目标的第
Figure 100002_DEST_PATH_IMAGE013
个观测弧段的数据行数，下标
Figure 100002_DEST_PATH_IMAGE014
表示观测弧段中第
Figure 100002_DEST_PATH_IMAGE015
行数据，
Figure 100002_DEST_PATH_IMAGE016
为观测历元时刻，
Figure 100002_DEST_PATH_IMAGE017
为赤经，
Figure 100002_DEST_PATH_IMAGE018
为赤纬，
Figure 100002_DEST_PATH_IMAGE019
与
Figure 100002_DEST_PATH_IMAGE020
为分别为每行数据观测历元时刻对应的观测卫星的位置和速度矢量。
在其中一个实施例，步骤2的实现过程为：
步骤2.1，采用二次多项式分别对每个观测弧段中赤经、赤纬关于时间的函数式进行拟合，设赤经
Figure 100002_DEST_PATH_IMAGE021
和赤纬
Figure 100002_DEST_PATH_IMAGE022
对时间的函数
Figure 100002_DEST_PATH_IMAGE023
、
Figure 100002_DEST_PATH_IMAGE024
分别表示为：
Figure 100002_DEST_PATH_IMAGE025
（1）
其中，
Figure 100002_DEST_PATH_IMAGE026
、
Figure 100002_DEST_PATH_IMAGE027
、
Figure 100002_DEST_PATH_IMAGE028
、
Figure 100002_DEST_PATH_IMAGE029
、
Figure 100002_DEST_PATH_IMAGE030
、
Figure 100002_DEST_PATH_IMAGE031
为多项式待定系数，各待定系数的初值取为：
Figure 100002_DEST_PATH_IMAGE032
（2）
Figure 100002_DEST_PATH_IMAGE033
对
Figure 100002_DEST_PATH_IMAGE034
、
Figure 100002_DEST_PATH_IMAGE035
、
Figure 100002_DEST_PATH_IMAGE036
的偏导数以及
Figure 100002_DEST_PATH_IMAGE037
对
Figure 100002_DEST_PATH_IMAGE038
、
Figure 100002_DEST_PATH_IMAGE039
、
Figure 100002_DEST_PATH_IMAGE040
的偏导数分别为：
Figure 100002_DEST_PATH_IMAGE041
（3）
因此可以使用最小二乘法得到对
Figure 100002_DEST_PATH_IMAGE042
、
Figure 100002_DEST_PATH_IMAGE043
、
Figure 269802DEST_PATH_IMAGE036
初值的改进量
Figure 100002_DEST_PATH_IMAGE044
、
Figure 100002_DEST_PATH_IMAGE045
、
Figure 100002_DEST_PATH_IMAGE046
，以及为对
Figure 100002_DEST_PATH_IMAGE047
、
Figure 100002_DEST_PATH_IMAGE048
、
Figure 100002_DEST_PATH_IMAGE049
初值的改进量
Figure 100002_DEST_PATH_IMAGE050
、
Figure 100002_DEST_PATH_IMAGE051
、
Figure 100002_DEST_PATH_IMAGE052
：
Figure 100002_DEST_PATH_IMAGE053
（4）
其中，
Figure 100002_DEST_PATH_IMAGE054
是
Figure 100002_DEST_PATH_IMAGE055
的矩阵，
Figure 100002_DEST_PATH_IMAGE056
为
Figure 100002_DEST_PATH_IMAGE057
的转置矩阵，上标-1表示矩阵的求逆运算，
Figure 100002_DEST_PATH_IMAGE058
是
Figure 100002_DEST_PATH_IMAGE059
维的向量，
Figure 100002_DEST_PATH_IMAGE060
为赤经的多项式预测值，
Figure 100002_DEST_PATH_IMAGE061
是
Figure 100002_DEST_PATH_IMAGE062
维的向量，
Figure 100002_DEST_PATH_IMAGE063
为赤纬的多项式预测值；
则将
Figure 100002_DEST_PATH_IMAGE064
、
Figure 100002_DEST_PATH_IMAGE065
、
Figure 100002_DEST_PATH_IMAGE066
以及
Figure 100002_DEST_PATH_IMAGE067
、
Figure 100002_DEST_PATH_IMAGE068
、
Figure 100002_DEST_PATH_IMAGE069
更新为：
Figure 100002_DEST_PATH_IMAGE070
（5）
重复式(4)-(5)的最小二乘法计算以及
Figure 100002_DEST_PATH_IMAGE071
、
Figure 100002_DEST_PATH_IMAGE072
、
Figure 100002_DEST_PATH_IMAGE073
、
Figure 100002_DEST_PATH_IMAGE074
、
Figure 100002_DEST_PATH_IMAGE075
、
Figure 100002_DEST_PATH_IMAGE076
更新过程直至
Figure 100002_DEST_PATH_IMAGE077
、
Figure 100002_DEST_PATH_IMAGE078
小于设定的阈值，阈值取为
Figure 100002_DEST_PATH_IMAGE079
，最终即得到拟合出的
Figure 100002_DEST_PATH_IMAGE080
、
Figure 100002_DEST_PATH_IMAGE081
、
Figure 100002_DEST_PATH_IMAGE082
以及
Figure 100002_DEST_PATH_IMAGE083
、
Figure 100002_DEST_PATH_IMAGE084
、
Figure 100002_DEST_PATH_IMAGE085
；
步骤2.2，定义一个观测弧段的中间时刻为
Figure 100002_DEST_PATH_IMAGE086
，其中，
Figure 100002_DEST_PATH_IMAGE087
表示对应观测弧段的中间行序号，由此对于每一个观测弧段
Figure 100002_DEST_PATH_IMAGE088
都有一个对应的中间时刻数据点
Figure 100002_DEST_PATH_IMAGE089
，为：
Figure 100002_DEST_PATH_IMAGE090
（6）
其中，
Figure 100002_DEST_PATH_IMAGE091
为中间时刻赤经，
Figure 100002_DEST_PATH_IMAGE092
为中间时刻赤纬，
Figure 100002_DEST_PATH_IMAGE093
为中间时刻赤经变化率，
Figure 100002_DEST_PATH_IMAGE094
为中间时刻赤纬变化率，
Figure 100002_DEST_PATH_IMAGE095
，
Figure 100002_DEST_PATH_IMAGE096
分别为中间时刻对应光学观测卫星的位置矢量与速度矢量；其计算式如下：
Figure 100002_DEST_PATH_IMAGE097
（7）
步骤2.3，对于每个观测时刻的观测数据点，可以通过式(1)得到对应时刻的赤经、赤纬拟合值，将对应时刻的赤经、赤纬拟合值与真实观测值作差，可以得到赤经、赤纬的残差，根据总体标准差计算公式：
Figure 100002_DEST_PATH_IMAGE098
（8）
由此可以计算得到一个弧段的拟合值与实际观测值残差的标准差
Figure 100002_DEST_PATH_IMAGE099
，其中，
Figure 100002_DEST_PATH_IMAGE100
表示第
Figure 100002_DEST_PATH_IMAGE101
个观测数据的残差，
Figure 100002_DEST_PATH_IMAGE102
为残差均值，
Figure 100002_DEST_PATH_IMAGE103
为观测数据点个数。若某观测数据点的残差大于
Figure 100002_DEST_PATH_IMAGE104
则认定该点为坏点，将该观测数据点从对应观测弧段中剔除，否则可能会影响后续轨道关联与聚类的效果。
在其中一个实施例，步骤3的实现过程为：
步骤3.1，估算被观测目标半长轴的取值区间
Figure 100002_DEST_PATH_IMAGE105
、偏心率的取值区间
Figure 100002_DEST_PATH_IMAGE106
、相对于观测卫星的斜距
Figure 100002_DEST_PATH_IMAGE107
的取值区间
Figure 100002_DEST_PATH_IMAGE108
以及斜距变化率
Figure 100002_DEST_PATH_IMAGE109
的取值区间
Figure 100002_DEST_PATH_IMAGE110
，具体包括如下步骤：
步骤3.1.1，根据被观测目标大致运行轨道区间先验信息，可估算被观测目标半长轴的取值区间
Figure 100002_DEST_PATH_IMAGE111
与偏心率的取值区间
Figure 100002_DEST_PATH_IMAGE112
。如被观测目标为近GEO目标，则其半长轴与偏心率的取值区间可以分别取为
Figure 100002_DEST_PATH_IMAGE113
与
Figure 100002_DEST_PATH_IMAGE114
。
步骤3.1.2，根据天基观测卫星运行轨道与被观测目标大致运行轨道区间等先验信息，可估算被观测目标相对于观测卫星的斜距
Figure 100002_DEST_PATH_IMAGE115
的取值区间
Figure 100002_DEST_PATH_IMAGE116
以及斜距变化率
Figure 100002_DEST_PATH_IMAGE117
的取值区间
Figure 100002_DEST_PATH_IMAGE118
。斜距与斜距变化率的取值区间
Figure 100002_DEST_PATH_IMAGE119
与
Figure 100002_DEST_PATH_IMAGE120
可按下式进行估算：
Figure 100002_DEST_PATH_IMAGE121
（9）
其中，
Figure 100002_DEST_PATH_IMAGE122
与
Figure 100002_DEST_PATH_IMAGE123
分别代表位置与速度的大小，上标s与t分别代表天基观测卫星与被观测目标，下标peri与apo分别代表近地点与远地点，如
Figure 100002_DEST_PATH_IMAGE124
表示天基观测卫星在远地点处速度的大小。由于被观测目标在近地点与远地点处的准确位置速度无从得知，此处采用大致的估计值即可。
需要注意的是，若被观测目标无任何可用先验信息，则按绕地运行卫星应满足的基本条件对上述取值区间进行估算即可。
步骤3.2，根据观测弧段的中间时刻赤经
Figure 100002_DEST_PATH_IMAGE125
、中间时刻赤纬
Figure 100002_DEST_PATH_IMAGE126
及其变化率信息
Figure 100002_DEST_PATH_IMAGE127
、
Figure 100002_DEST_PATH_IMAGE128
以及被观测目标半长轴的取值区间
Figure 100002_DEST_PATH_IMAGE129
划定观测弧段对应被观测目标在斜距
Figure 100002_DEST_PATH_IMAGE130
与斜距变化率
Figure 100002_DEST_PATH_IMAGE131
平面（简称为
Figure 100002_DEST_PATH_IMAGE132
平面）上的容许域范围，具体包括：
为了更好的理解如何对容许域进行划定，需要对一些将会用到的变量符号进行介绍：
设观测弧段对应被观测目标位置与速度矢量分别为与，则其与天基观测卫星的位置与速度矢量
Figure 100002_DEST_PATH_IMAGE133
与
Figure 100002_DEST_PATH_IMAGE134
存在如下关系：
Figure 100002_DEST_PATH_IMAGE135
（10）
其中，
Figure 100002_DEST_PATH_IMAGE136
与
Figure 100002_DEST_PATH_IMAGE137
代表被测目标相对于天基观测卫星的位置与速度矢量；
相对位置速度
Figure 100002_DEST_PATH_IMAGE138
与
Figure 100002_DEST_PATH_IMAGE139
可以用斜距
Figure 100002_DEST_PATH_IMAGE140
、赤经
Figure 100002_DEST_PATH_IMAGE141
、赤纬
Figure 100002_DEST_PATH_IMAGE142
及其变化率
Figure 100002_DEST_PATH_IMAGE143
、
Figure 100002_DEST_PATH_IMAGE144
来表示，为：
Figure 100002_DEST_PATH_IMAGE145
（11）
其中，中间参数
Figure 100002_DEST_PATH_IMAGE146
、
Figure 100002_DEST_PATH_IMAGE147
与
Figure 100002_DEST_PATH_IMAGE148
的定义如下：
Figure 100002_DEST_PATH_IMAGE149
（12）
此外还需要定义一系列辅助标量，为：
Figure 100002_DEST_PATH_IMAGE150
（13）
经推导，观测弧段对应被观测目标的斜距与斜距变化率应当满足式(14)，有关容许域的具体推导过程可阅读以下文献：Milani A, Gronchi G F, De’ Michieli Vitturi,M, Knežević Z. Orbit Determination with Very Short Arcs. I Admissible Regions[J]. Celestial Mechanics and Dynamical Astronomy, 2004, 90(1-2):59-87.
Figure 100002_DEST_PATH_IMAGE151
（14）
其中，
Figure 100002_DEST_PATH_IMAGE152
为地心引力常数，函数关系
Figure 100002_DEST_PATH_IMAGE153
定义如下：
Figure 100002_DEST_PATH_IMAGE154
（15）
当式(14)中的半长轴
Figure 100002_DEST_PATH_IMAGE155
分别取区间
Figure 100002_DEST_PATH_IMAGE156
的上界和下界时，会在
Figure 100002_DEST_PATH_IMAGE157
平面上分别得到两条曲线，设在
Figure 100002_DEST_PATH_IMAGE158
平面上这两条曲线之间所圈定的区域为
Figure 100002_DEST_PATH_IMAGE159
，被测目标的
Figure 100002_DEST_PATH_IMAGE160
与
Figure 100002_DEST_PATH_IMAGE161
只能够在区域
Figure 100002_DEST_PATH_IMAGE162
内选取。
步骤3.3，根据观测弧段的中间时刻赤经
Figure 100002_DEST_PATH_IMAGE163
、中间时刻赤纬
Figure 100002_DEST_PATH_IMAGE164
及其变化率信息
Figure 100002_DEST_PATH_IMAGE165
、
Figure 100002_DEST_PATH_IMAGE166
以及被观测目标偏心率的取值区间
Figure 100002_DEST_PATH_IMAGE167
，划定观测弧段对应被观测目标在斜距
Figure 100002_DEST_PATH_IMAGE168
与斜距变化率
Figure 100002_DEST_PATH_IMAGE169
平面上的容许域范围，具体包括：
首先对一些将会用到的辅助矢量进行定义：
Figure 100002_DEST_PATH_IMAGE170
（16）
此外还需要定义一系列辅助标量：
Figure 100002_DEST_PATH_IMAGE171
（17）
经推导，观测弧段对应被观测目标的斜距与斜距变化率应当满足式(16)，有关容许域的具体推导过程可阅读以下文献：Milani A, Gronchi G F, De’ Michieli Vitturi,M, Knežević Z. Orbit Determination with Very Short Arcs. I Admissible Regions[J]. Celestial Mechanics and Dynamical Astronomy, 2004, 90(1-2):59-87.
Figure 100002_DEST_PATH_IMAGE172
（18）
其中，
Figure 100002_DEST_PATH_IMAGE173
为地心引力常数，函数关系
Figure 100002_DEST_PATH_IMAGE174
与
Figure 100002_DEST_PATH_IMAGE175
定义如下：
Figure 100002_DEST_PATH_IMAGE176
（19）
当式(16)中的偏心率
Figure 100002_DEST_PATH_IMAGE177
分别取区间
Figure 100002_DEST_PATH_IMAGE178
的上界和下界时，会在
Figure 100002_DEST_PATH_IMAGE179
平面上分别得到两条曲线，设在
Figure 100002_DEST_PATH_IMAGE180
平面上这两条曲线之间所圈定的区域为
Figure 100002_DEST_PATH_IMAGE181
，被测目标的
Figure 100002_DEST_PATH_IMAGE182
与
Figure 100002_DEST_PATH_IMAGE183
只能够在区域
Figure 100002_DEST_PATH_IMAGE184
内选取。
步骤3.4，设斜距
Figure 100002_DEST_PATH_IMAGE185
的取值区间
Figure 100002_DEST_PATH_IMAGE186
以及斜距变化率
Figure 100002_DEST_PATH_IMAGE187
的取值区间
Figure 100002_DEST_PATH_IMAGE188
在斜距
Figure 100002_DEST_PATH_IMAGE189
与斜距变化率
Figure 100002_DEST_PATH_IMAGE190
平面上所圈定的区域为
Figure 100002_DEST_PATH_IMAGE191
，则各观测弧段对应目标轨道在斜距与斜距变化率平面上的容许域范围为区域
Figure DEST_PATH_IMAGE192
、区域
Figure DEST_PATH_IMAGE193
与区域
Figure DEST_PATH_IMAGE194
的交集，即：
Figure DEST_PATH_IMAGE195
（20）
其中，
Figure DEST_PATH_IMAGE196
为观测弧段对应目标轨道在斜距与斜距变化率平面上的容许域范围。
在其中一个实施例，步骤4的实现过程为：
步骤4.1，对待关联观测弧段两两在步骤3中划定的容许域范围内对斜距与斜距变化率组合
Figure DEST_PATH_IMAGE197
进行优化，找到使得赤经赤纬预报值与实际观测值马氏距离
Figure DEST_PATH_IMAGE198
最小的斜距与斜距变化率组合
Figure DEST_PATH_IMAGE199
。各种优化方法已经是航天领域乃至整个科学界常用的工具，常见的优化方法有梯度下降法、牛顿法和拟牛顿法、共轭梯度法等等。在Matlab自带的优化工具箱（Optimization Tool）中有多个可供直接使用的优化函数，还可进一步自行选择不同的优化算法。比如，调用其中的fmincon函数就可实现对该问题的优化，此处不再对优化算法的实现过程进行详述；
步骤4.2，由斜距与斜距变化率组合
Figure DEST_PATH_IMAGE200
计算得到优化指标马氏距离
Figure DEST_PATH_IMAGE201
，其计算步骤为：
步骤4.2.1，根据斜距
Figure DEST_PATH_IMAGE202
、斜距变化率
Figure DEST_PATH_IMAGE203
、中间时刻赤经
Figure DEST_PATH_IMAGE204
、中间时刻赤纬
Figure DEST_PATH_IMAGE205
及其变化率
Figure DEST_PATH_IMAGE206
、
Figure DEST_PATH_IMAGE207
计算观测弧段对应轨道状态。设选取的两个待关联观测弧段分别为E与F，两弧段对应容许域分别为
Figure DEST_PATH_IMAGE208
与
Figure DEST_PATH_IMAGE209
，对于容许域
Figure DEST_PATH_IMAGE210
中所挑选出的一组斜距与斜距变化率
Figure DEST_PATH_IMAGE211
，可以计算出观测弧段E在弧段中间时刻
Figure DEST_PATH_IMAGE212
所对应的一组轨道状态
Figure DEST_PATH_IMAGE213
，计算式如下：
Figure DEST_PATH_IMAGE214
（21）
步骤4.2.2，构建观测弧段E在
Figure DEST_PATH_IMAGE215
时刻对应轨道状态在当地轨道坐标系下的轨道状态协方差矩阵
Figure DEST_PATH_IMAGE216
。通过对整个观测弧段的数据点进行多项式拟合所得到的中间时刻赤经
Figure DEST_PATH_IMAGE217
、中间时刻赤纬
Figure DEST_PATH_IMAGE218
及其变化率
Figure DEST_PATH_IMAGE219
、
Figure DEST_PATH_IMAGE220
的标准差可根据原始数据单点观测标准差进行估算，估算公式为：
Figure DEST_PATH_IMAGE221
（22）
其中，
Figure DEST_PATH_IMAGE222
与
Figure DEST_PATH_IMAGE223
分别为原始数据赤经赤纬单点观测标准差，
Figure DEST_PATH_IMAGE224
是观测弧段的数据点数量，
Figure DEST_PATH_IMAGE225
是观测弧段首尾数据点横跨的时间长度。则观测弧段E对应轨道状态在观测空间内的协方差矩阵
Figure DEST_PATH_IMAGE226
可表示为：
Figure DEST_PATH_IMAGE227
（23）
则
Figure DEST_PATH_IMAGE228
可通过
Figure DEST_PATH_IMAGE229
由下式计算得到，为：
Figure DEST_PATH_IMAGE230
（24）
其中，
Figure DEST_PATH_IMAGE231
与
Figure DEST_PATH_IMAGE232
分别为观测空间到地心惯性系与地心惯性系到当地轨道坐标系的转换矩阵。
Figure DEST_PATH_IMAGE233
的计算式为：
Figure DEST_PATH_IMAGE234
（25）
式中
Figure DEST_PATH_IMAGE235
、
Figure DEST_PATH_IMAGE236
、
Figure DEST_PATH_IMAGE237
和
Figure DEST_PATH_IMAGE238
的定义为：
Figure DEST_PATH_IMAGE239
（26）
Figure DEST_PATH_IMAGE240
的计算式为：
Figure DEST_PATH_IMAGE241
（27）
其中，变量
Figure DEST_PATH_IMAGE242
到
Figure DEST_PATH_IMAGE243
的定义如下：
Figure DEST_PATH_IMAGE244
（28）
步骤4.2.3，利用航天器轨道预报与偏差演化算法将观测弧段E在
Figure DEST_PATH_IMAGE245
时刻对应轨道状态
Figure DEST_PATH_IMAGE246
与轨道状态协方差矩阵
Figure DEST_PATH_IMAGE247
预报至观测弧段F对应弧段中间时刻
Figure DEST_PATH_IMAGE248
得到预报轨道状态
Figure DEST_PATH_IMAGE249
与预报轨道状态协方差矩阵
Figure DEST_PATH_IMAGE250
。航天器轨道预报与偏差演化算法是航天领域的成熟算法，有多重基于不同模型的算法，此处采用更为贴合实际的非线性轨道预报与偏差演化算法可提升最终关联聚类精度，有关非线性轨道预报与偏差演化算法的详细内容可以参考以下文献：杨震.非线性轨道机动瞄准与偏差演化分析方法[D]. 长沙：国防科技大学研究生院博士学位论文，2018,04.
步骤4.2.4，将经预报后得到的
Figure DEST_PATH_IMAGE251
与
Figure DEST_PATH_IMAGE252
重新转换至观测空间，得到
Figure DEST_PATH_IMAGE253
时刻赤经赤纬的预测值
Figure DEST_PATH_IMAGE254
、
Figure DEST_PATH_IMAGE255
以及观测空间内的预报协方差矩阵
Figure DEST_PATH_IMAGE256
。赤经赤纬预测值的计算式如下：
Figure DEST_PATH_IMAGE257
（29）
预报协方差矩阵
Figure DEST_PATH_IMAGE258
可通过
Figure DEST_PATH_IMAGE259
由下式计算得到：
Figure DEST_PATH_IMAGE260
（30）
其中，
Figure DEST_PATH_IMAGE261
与
Figure DEST_PATH_IMAGE262
分别为当地轨道坐标系到地心惯性系与地心惯性系到观测空间的转换矩阵。
Figure DEST_PATH_IMAGE263
的计算式为：
Figure DEST_PATH_IMAGE264
（31）
其中，变量
Figure DEST_PATH_IMAGE265
到
Figure DEST_PATH_IMAGE266
的定义如下：
Figure DEST_PATH_IMAGE267
（32）
Figure DEST_PATH_IMAGE268
的计算式为：
Figure DEST_PATH_IMAGE269
（33）
式中变量
Figure DEST_PATH_IMAGE270
、
Figure DEST_PATH_IMAGE271
、
Figure DEST_PATH_IMAGE272
的定义如下：
Figure DEST_PATH_IMAGE273
（34）
需要注意的是，由于计算马氏距离时仅用到了赤经赤纬而没有涉及赤经赤纬变化率，步骤4.2.4中计算得到的预报协方差矩阵
Figure DEST_PATH_IMAGE274
为
Figure DEST_PATH_IMAGE275
矩阵。
步骤4.2.5，计算由观测弧段E预报得到的
Figure DEST_PATH_IMAGE276
时刻赤经赤纬预测值
Figure DEST_PATH_IMAGE277
、
Figure DEST_PATH_IMAGE278
与观测弧段F在
Figure DEST_PATH_IMAGE279
时刻赤经赤纬拟合值
Figure DEST_PATH_IMAGE280
、
Figure DEST_PATH_IMAGE281
的马氏距离
Figure DEST_PATH_IMAGE282
。马氏距离
Figure DEST_PATH_IMAGE283
计算式如下：
Figure DEST_PATH_IMAGE284
（35）
其中，上标T表示矩阵转置。计算得到的马氏距离
Figure 233909DEST_PATH_IMAGE283
是工程中一种被广泛用作评定数据之间的相似度的无量纲指标，因此此处不对马氏距离进行详细介绍。
步骤4.3，将所有待关联观测弧段两两之间按步骤4.2所述计算马氏距离
Figure 587530DEST_PATH_IMAGE283
，并通过优化到两待关联观测弧段间最小马氏距离
Figure DEST_PATH_IMAGE285
，将每组观测弧段间的最小马氏距离进行记录并保存。
在其中一个实施例，步骤5的实现过程为：
以步骤4中所记录的两观测弧段间最小马氏距离为关联判别依据，一一进行判别，得到观测弧段两两关联匹配结果。可采用工程上常用的马氏距离判别依据进行判别，即
Figure DEST_PATH_IMAGE286
（36）
若最小马氏距离小于等于5，则认为这两个观测弧段之间关联成功，可能为对同一空间目标进行观测所产生的观测弧段。
其中一个实施例，步骤6的实现过程为：
步骤6.1，构建观测弧段关联矩阵
Figure DEST_PATH_IMAGE287
，其中，为待关联弧段数量。关联矩阵
Figure 933192DEST_PATH_IMAGE287
中第
Figure DEST_PATH_IMAGE288
行第
Figure DEST_PATH_IMAGE289
列元素
Figure DEST_PATH_IMAGE290
按照如下规则进行取值：
Figure DEST_PATH_IMAGE291
（37）
其中，
Figure DEST_PATH_IMAGE292
表示第
Figure DEST_PATH_IMAGE293
个待关联弧段与第
Figure DEST_PATH_IMAGE294
个待关联弧段间的最小马氏距离；
步骤6.2，利用BEA（Bond Energy Algorithm）算法对观测弧段关联矩阵
Figure DEST_PATH_IMAGE295
进行行列变换，将观测弧段关联矩阵
Figure 882824DEST_PATH_IMAGE295
变换成观测弧段聚类矩阵
Figure DEST_PATH_IMAGE296
。BEA算法是一种广泛应用在分布式数据库系统中大型表的纵向划分的算法，还可以实现矩阵元素的聚类。有关BEA算法的原理与具体实现步骤，可以参考以下文献：Ozsu M T, Valduriez P.Principles of distributed database systems[M]. [S.l.]：Prentice-Hall，1999.
其中一个实施例，步骤7的实现过程为：
步骤7.1，构建聚类分割辅助序列
Figure DEST_PATH_IMAGE297
与
Figure DEST_PATH_IMAGE298
。为了实现对聚类矩阵
Figure 704367DEST_PATH_IMAGE296
的分割，首先定义两个各具有
Figure DEST_PATH_IMAGE299
个元素的序列
Figure 963310DEST_PATH_IMAGE297
与
Figure 248798DEST_PATH_IMAGE298
，序列
Figure 516968DEST_PATH_IMAGE297
与
Figure 322245DEST_PATH_IMAGE298
中元素的取值与聚类矩阵
Figure 548827DEST_PATH_IMAGE296
中的元素有关，存在如下关系：
Figure DEST_PATH_IMAGE300
（38）
其中，
Figure DEST_PATH_IMAGE301
表示聚类矩阵
Figure DEST_PATH_IMAGE302
中的第
Figure DEST_PATH_IMAGE303
行第
Figure DEST_PATH_IMAGE304
列元素，
Figure DEST_PATH_IMAGE305
表示序列
Figure DEST_PATH_IMAGE306
中的第
Figure DEST_PATH_IMAGE307
个元素，
Figure DEST_PATH_IMAGE308
表示序列
Figure DEST_PATH_IMAGE309
中的第
Figure 806764DEST_PATH_IMAGE307
个元素；
步骤7.2，根据序列
Figure DEST_PATH_IMAGE310
与
Figure DEST_PATH_IMAGE311
中元素变化规律对聚类矩阵
Figure DEST_PATH_IMAGE312
进行分割。当
Figure 222833DEST_PATH_IMAGE310
与
Figure 397462DEST_PATH_IMAGE311
中元素变化规律满足以下条件时对聚类矩阵
Figure 467049DEST_PATH_IMAGE312
进行分割：
Figure DEST_PATH_IMAGE313
（39）
满足以上条件的
Figure DEST_PATH_IMAGE314
值为分割点。若仅存在一个分割点
Figure 540179DEST_PATH_IMAGE314
，则聚类矩阵
Figure DEST_PATH_IMAGE315
以第
Figure DEST_PATH_IMAGE316
行第
Figure 228780DEST_PATH_IMAGE316
列元素为界，被分割为由第
Figure DEST_PATH_IMAGE317
行第
Figure DEST_PATH_IMAGE318
列元素构成的
Figure DEST_PATH_IMAGE319
与由第
Figure DEST_PATH_IMAGE320
行第
Figure DEST_PATH_IMAGE321
列元素构成的
Figure DEST_PATH_IMAGE322
两个聚类子矩阵。若存在多个分割点，分割方式以此类推。
步骤7.3，对于位于同一聚类子矩阵内的观测弧段，视为聚类成功，认定这些观测弧段是对同一空间目标进行观测所产生的观测弧段。由此最终得到观测弧段关联聚类结果，实现对属于同一空间目标观测弧段的关联聚类。
本发明适用于短弧天基光学观测条件下的观测弧段关联与聚类，通过构建观测弧段的容许域，利用优化方法、BEA算法、轨道及其偏差预报技术实现了对观测弧段的关联与聚类。由于容许域包含广泛的特性，降低了方法的漏警率；由于选取马氏距离这一无量纲量为判别依据，巧妙规避了关联检测阈值设计问题；由于本发明并未对采用何种优化方法与轨道及其偏差预报方法作出限制，因此使用者可以根据实际需求选取不同的计算方法，从而实现计算效率与计算准确率的兼顾。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图示出的结构获得其他的附图。
图1为本发明实施例中基于非线性偏差演化的天基光学观测短弧关联与聚类方法的流程示意图。
本发明目的的实现、功能特点及优点将结合实施例，参照附图做进一步说明。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明各个实施例之间的技术方案可以相互结合，但是必须是以本领域普通技术人员能够实现为基础，当技术方案的结合出现相互矛盾或无法实现时应当认为这种技术方案的结合不存在，也不在本发明要求的保护范围之内。
本实施例公开了一种基于非线性偏差演化的天基光学观测短弧关联与聚类方法，主要用于解决现有技术中难以属于同一空间目标的观测弧段进行关联并聚类的问题，同时兼顾算法的计算正确率与计算效率。
参考图1，本实施例中的基于非线性偏差演化的天基光学观测短弧关联与聚类方法具体包括如下步骤1-步骤7。
步骤1，获取天基短弧光学观测数据。具体地，利用天基光学观测卫星进行天基光学观测，获取多组分别属于不同空间目标的原始观测短弧片段，也称观测弧段，每个观测弧段数据包括若干个观测数据点，每个观测数据点由被观测目标相对于低轨光学观测卫星的赤经、赤纬、观测时刻以及观测平台的位置速度信息组成。
假设利用运行在轨道高度为800km太阳同步轨道上的某低轨光学观测卫星对4颗运行在近GEO轨道上的空间目标进行为期7天的光学观测，角度观测误差为3个角秒，观测起止时间分别为2019.12.21.12:00:00至2019.12.28.12:00:00，观测得到58组原始观测短弧片段，也称观测弧段。其中，第1~15号观测弧段由对同一空间目标进行观测得到，第16~29号观测弧段由对同一空间目标进行观测得到，第30~44号观测弧段由对同一空间目标进行观测得到，第45~58号观测弧段由对同一空间目标进行观测得到。该低轨光学观测卫星在
Figure DEST_PATH_IMAGE323
初始时刻的轨道根数为：
Figure DEST_PATH_IMAGE324
。每个观测弧段数据包括若干个观测数据点，每个观测数据点由被观测目标相对于低轨光学观测卫星的赤经、赤纬、观测时刻以及观测平台的位置速度信息组成。
因此本实施例中步骤1获取天基短弧光学观测数据的过程为：
已知通过光学观测卫星上安装的天基观测设备，并进行弧段关联匹配后，得到的多组分别属于不同空间目标的天基测角数据
Figure DEST_PATH_IMAGE325
，即多组分别属于不同空间目标的观测弧段为
Figure DEST_PATH_IMAGE326
，
Figure DEST_PATH_IMAGE327
，
Figure DEST_PATH_IMAGE328
，其中，
Figure DEST_PATH_IMAGE329
为空间目标数量，
Figure DEST_PATH_IMAGE330
为观测弧段的数量；
Figure DEST_PATH_IMAGE331
为第
Figure DEST_PATH_IMAGE332
个空间目标的第
Figure DEST_PATH_IMAGE333
个观测弧段，
Figure DEST_PATH_IMAGE334
，其中，
Figure DEST_PATH_IMAGE335
为第
Figure 290539DEST_PATH_IMAGE332
个空间目标的第
Figure 796607DEST_PATH_IMAGE333
个观测弧段的数据行数，下标
Figure DEST_PATH_IMAGE336
表示观测弧段中第
Figure 560295DEST_PATH_IMAGE336
行数据，
Figure DEST_PATH_IMAGE337
为观测历元时刻，
Figure DEST_PATH_IMAGE338
为赤经，
Figure DEST_PATH_IMAGE339
为赤纬，
Figure DEST_PATH_IMAGE340
与
Figure DEST_PATH_IMAGE341
为分别为每行数据观测历元时刻对应的观测卫星的位置和速度矢量。
步骤2，进行数据预处理。具体地，采用二次多项式分别对每个观测弧段中赤经、赤纬关于时间的函数式进行拟合，得到赤经、赤纬随时间的变化率信息，并对存在明显异常的观测数据点进行剔除。其具体实施过程为：
步骤2.1，采用二次多项式分别对每个观测弧段中赤经、赤纬关于时间的函数式进行拟合，设赤经
Figure DEST_PATH_IMAGE342
和赤纬
Figure DEST_PATH_IMAGE343
对时间的函数
Figure DEST_PATH_IMAGE344
、
Figure DEST_PATH_IMAGE345
分别表示为：
Figure DEST_PATH_IMAGE346
（1）
其中，
Figure DEST_PATH_IMAGE347
、
Figure DEST_PATH_IMAGE348
、
Figure DEST_PATH_IMAGE349
、
Figure DEST_PATH_IMAGE350
、
Figure DEST_PATH_IMAGE351
、
Figure DEST_PATH_IMAGE352
为多项式待定系数，各待定系数的初值取为：
Figure DEST_PATH_IMAGE353
（2）
Figure DEST_PATH_IMAGE354
对
Figure DEST_PATH_IMAGE355
、
Figure DEST_PATH_IMAGE356
、
Figure DEST_PATH_IMAGE357
的偏导数以及
Figure DEST_PATH_IMAGE358
对
Figure DEST_PATH_IMAGE359
、
Figure DEST_PATH_IMAGE360
、
Figure DEST_PATH_IMAGE361
的偏导数分别为：
Figure DEST_PATH_IMAGE362
（3）
因此可以使用最小二乘法得到对
Figure DEST_PATH_IMAGE363
、
Figure DEST_PATH_IMAGE364
、
Figure DEST_PATH_IMAGE365
初值的改进量
Figure DEST_PATH_IMAGE366
、
Figure DEST_PATH_IMAGE367
、
Figure DEST_PATH_IMAGE368
，以及为对
Figure DEST_PATH_IMAGE369
、
Figure 168432DEST_PATH_IMAGE360
、
Figure 52074DEST_PATH_IMAGE361
初值的改进量
Figure DEST_PATH_IMAGE370
、
Figure DEST_PATH_IMAGE371
、
Figure DEST_PATH_IMAGE372
：
Figure DEST_PATH_IMAGE373
（4）
其中，
Figure DEST_PATH_IMAGE374
是
Figure DEST_PATH_IMAGE375
的矩阵，
Figure DEST_PATH_IMAGE376
为
Figure DEST_PATH_IMAGE377
的转置矩阵，上标-1表示矩阵的求逆运算，
Figure DEST_PATH_IMAGE378
是
Figure DEST_PATH_IMAGE379
维的向量，
Figure DEST_PATH_IMAGE380
为赤经的多项式预测值，
Figure DEST_PATH_IMAGE381
是
Figure DEST_PATH_IMAGE382
维的向量，
Figure DEST_PATH_IMAGE383
为赤纬的多项式预测值；
则将
Figure DEST_PATH_IMAGE384
、
Figure DEST_PATH_IMAGE385
、
Figure DEST_PATH_IMAGE386
以及
Figure DEST_PATH_IMAGE387
、
Figure DEST_PATH_IMAGE388
、
Figure DEST_PATH_IMAGE389
更新为：
Figure DEST_PATH_IMAGE390
（5）
重复式(4)-(5)的最小二乘法计算以及
Figure DEST_PATH_IMAGE391
、
Figure DEST_PATH_IMAGE392
、
Figure DEST_PATH_IMAGE393
、
Figure DEST_PATH_IMAGE394
、
Figure DEST_PATH_IMAGE395
、
Figure DEST_PATH_IMAGE396
更新过程直至
Figure DEST_PATH_IMAGE397
、
Figure DEST_PATH_IMAGE398
小于设定的阈值，阈值取为
Figure DEST_PATH_IMAGE399
，最终即得到拟合出的
Figure 591027DEST_PATH_IMAGE391
、
Figure 560120DEST_PATH_IMAGE392
、
Figure 715158DEST_PATH_IMAGE393
以及
Figure 718886DEST_PATH_IMAGE394
、
Figure 363494DEST_PATH_IMAGE395
、
Figure 570616DEST_PATH_IMAGE396
；
步骤2.2，定义一个观测弧段的中间时刻为
Figure DEST_PATH_IMAGE400
，其中，
Figure DEST_PATH_IMAGE401
表示对应观测弧段的中间行序号，由此对于每一个观测弧段
Figure DEST_PATH_IMAGE402
都有一个对应的中间时刻数据点
Figure DEST_PATH_IMAGE403
，为：
Figure DEST_PATH_IMAGE404
（6）
其中，
Figure DEST_PATH_IMAGE405
为中间时刻赤经，
Figure DEST_PATH_IMAGE406
为中间时刻赤纬，
Figure DEST_PATH_IMAGE407
为中间时刻赤经变化率，
Figure DEST_PATH_IMAGE408
为中间时刻赤纬变化率，
Figure DEST_PATH_IMAGE409
，
Figure DEST_PATH_IMAGE410
分别为中间时刻对应光学观测卫星的位置矢量与速度矢量；其计算式如下：
Figure DEST_PATH_IMAGE411
（7）
步骤2.3，对于每个观测时刻的观测数据点，可以通过式(1)得到对应时刻的赤经、赤纬拟合值，将对应时刻的赤经、赤纬拟合值与真实观测值作差，可以得到赤经、赤纬的残差，根据总体标准差计算公式：
Figure DEST_PATH_IMAGE412
（8）
由此可以计算得到一个弧段的拟合值与实际观测值残差的标准差
Figure DEST_PATH_IMAGE413
，其中，
Figure DEST_PATH_IMAGE414
表示第
Figure DEST_PATH_IMAGE415
个残差，
Figure DEST_PATH_IMAGE416
为残差均值，
Figure DEST_PATH_IMAGE417
为观测数据点个数。若某观测数据点的残差大于
Figure DEST_PATH_IMAGE418
则认定该点为坏点，将该观测数据点从对应观测弧段中剔除，否则可能会影响后续轨道关联与聚类的效果。
步骤3，根据观测弧段特征信息与先验信息，划定各观测弧段对应容许域，具体地，根据步骤2处理后得到的各观测弧段中赤经、赤纬随时间的变化率信息，结合天基观测卫星运行轨道与被观测目标大致运行轨道区间等先验信息，对各观测弧段对应目标轨道在斜距与斜距变化率平面上的容许域范围进行划定。其具体实施过程为：
步骤3.1，估算被观测目标半长轴的取值区间
Figure DEST_PATH_IMAGE419
、偏心率的取值区间
Figure DEST_PATH_IMAGE420
、相对于观测卫星的斜距
Figure DEST_PATH_IMAGE421
的取值区间
Figure DEST_PATH_IMAGE422
以及斜距变化率
Figure DEST_PATH_IMAGE423
的取值区间
Figure DEST_PATH_IMAGE424
，具体包括如下步骤：
步骤3.1.1，根据被观测目标大致运行轨道区间先验信息，可估算被观测目标半长轴的取值区间
Figure DEST_PATH_IMAGE425
与偏心率的取值区间
Figure DEST_PATH_IMAGE426
。
本实施例中，被观测目标均运行在近GEO轨道，因此设定半长轴与偏心率的取值区间分别取为
Figure DEST_PATH_IMAGE427
与
Figure DEST_PATH_IMAGE428
。
步骤3.1.2，根据天基观测卫星运行轨道与被观测目标大致运行轨道区间等先验信息，估算被观测目标相对于观测卫星的斜距
Figure DEST_PATH_IMAGE429
的取值区间
Figure DEST_PATH_IMAGE430
以及斜距变化率的取值区间
Figure DEST_PATH_IMAGE431
。斜距与斜距变化率的取值区间
Figure DEST_PATH_IMAGE432
与
Figure DEST_PATH_IMAGE433
可按下式进行估算：
Figure DEST_PATH_IMAGE434
（9）
其中，
Figure DEST_PATH_IMAGE435
与
Figure DEST_PATH_IMAGE436
分别代表位置与速度的大小，上标s与t分别代表天基观测卫星与被观测目标，下标peri与apo分别代表近地点与远地点，如
Figure DEST_PATH_IMAGE437
表示天基观测卫星在远地点处速度的大小。由于被观测目标在近地点与远地点处的准确位置速度无从得知，此处采用大致的估计值即可。
本实施例中，估算得到斜距
Figure DEST_PATH_IMAGE438
与斜距变化率
Figure DEST_PATH_IMAGE439
的大致取值区间分别为：
Figure DEST_PATH_IMAGE440
与
Figure DEST_PATH_IMAGE441
。
步骤3.2，根据观测弧段的中间时刻赤经
Figure DEST_PATH_IMAGE442
、中间时刻赤纬
Figure DEST_PATH_IMAGE443
及其变化率信息
Figure DEST_PATH_IMAGE444
、
Figure DEST_PATH_IMAGE445
以及被观测目标半长轴的取值区间
Figure DEST_PATH_IMAGE446
划定观测弧段对应被观测目标在斜距
Figure DEST_PATH_IMAGE447
与斜距变化率
Figure DEST_PATH_IMAGE448
平面（简称为
Figure DEST_PATH_IMAGE449
平面）上的容许域范围，具体为：
为了更好的理解如何对容许域进行划定，需要对一些将会用到的变量符号进行介绍：
设观测弧段对应被观测目标位置与速度矢量分别为
Figure DEST_PATH_IMAGE450
与
Figure DEST_PATH_IMAGE451
，则其与天基观测卫星的位置与速度矢量
Figure DEST_PATH_IMAGE452
与
Figure DEST_PATH_IMAGE453
存在如下关系：
Figure DEST_PATH_IMAGE454
（10）
其中，
Figure DEST_PATH_IMAGE455
与
Figure DEST_PATH_IMAGE456
代表被测目标相对于天基观测卫星的位置与速度矢量；
相对位置速度
Figure DEST_PATH_IMAGE457
与
Figure DEST_PATH_IMAGE458
可以用斜距
Figure DEST_PATH_IMAGE459
、赤经
Figure DEST_PATH_IMAGE460
、赤纬
Figure DEST_PATH_IMAGE461
及其变化率
Figure DEST_PATH_IMAGE462
、
Figure DEST_PATH_IMAGE463
来表示，为：
Figure DEST_PATH_IMAGE464
（11）
其中，中间参数
Figure DEST_PATH_IMAGE465
、
Figure DEST_PATH_IMAGE466
与
Figure DEST_PATH_IMAGE467
的定义如下：
Figure DEST_PATH_IMAGE468
（12）
此外还需要定义一系列辅助标量，为：
Figure DEST_PATH_IMAGE469
（13）
经推导，观测弧段对应被观测目标的斜距与斜距变化率应当满足式(14)，有关容许域的具体推导过程可阅读以下文献：Milani A, Gronchi G F, De’ Michieli Vitturi,M, Knežević Z. Orbit Determination with Very Short Arcs. I Admissible Regions[J]. Celestial Mechanics and Dynamical Astronomy, 2004, 90(1-2):59-87.
Figure DEST_PATH_IMAGE470
（14）
其中，
Figure DEST_PATH_IMAGE471
为地心引力常数，函数关系
Figure DEST_PATH_IMAGE472
定义如下：
Figure DEST_PATH_IMAGE473
（15）
当式（14）中的半长轴
Figure DEST_PATH_IMAGE474
分别取区间
Figure DEST_PATH_IMAGE475
的上界和下界时，会在
Figure DEST_PATH_IMAGE476
平面上分别得到两条曲线，设在
Figure DEST_PATH_IMAGE477
平面上这两条曲线之间所圈定的区域为
Figure DEST_PATH_IMAGE478
，被测目标的
Figure DEST_PATH_IMAGE479
与
Figure DEST_PATH_IMAGE480
只能够在区域
Figure DEST_PATH_IMAGE481
内选取。
步骤3.3，根据观测弧段的中间时刻赤经
Figure DEST_PATH_IMAGE482
、中间时刻赤纬
Figure DEST_PATH_IMAGE483
及其变化率信息
Figure DEST_PATH_IMAGE484
、
Figure DEST_PATH_IMAGE485
以及被观测目标偏心率的取值区间
Figure DEST_PATH_IMAGE486
，划定观测弧段对应被观测目标在斜距
Figure DEST_PATH_IMAGE487
与斜距变化率
Figure DEST_PATH_IMAGE488
平面上的容许域范围，具体包括：
首先对一些将会用到的辅助矢量进行定义：
Figure DEST_PATH_IMAGE489
（16）
此外还需要定义一系列辅助标量：
Figure DEST_PATH_IMAGE490
（17）
经推导，观测弧段对应被观测目标的斜距与斜距变化率应当满足式(16)，有关容许域的具体推导过程可阅读以下文献：Milani A, Gronchi G F, De’ Michieli Vitturi,M, Knežević Z. Orbit Determination with Very Short Arcs. I Admissible Regions[J]. Celestial Mechanics and Dynamical Astronomy, 2004, 90(1-2):59-87.
Figure DEST_PATH_IMAGE491
（18）
其中，
Figure DEST_PATH_IMAGE492
为地心引力常数，函数关系
Figure DEST_PATH_IMAGE493
与
Figure DEST_PATH_IMAGE494
定义如下：
Figure DEST_PATH_IMAGE495
（19）
当式中的偏心率
Figure DEST_PATH_IMAGE496
分别取区间
Figure DEST_PATH_IMAGE497
的上界和下界时，会在
Figure DEST_PATH_IMAGE498
平面上分别得到两条曲线，设在
Figure 52980DEST_PATH_IMAGE498
平面上这两条曲线之间所圈定的区域为
Figure DEST_PATH_IMAGE499
，被测目标的
Figure DEST_PATH_IMAGE500
与
Figure DEST_PATH_IMAGE501
只能够在区域
Figure DEST_PATH_IMAGE502
内选取。
步骤3.4，设斜距
Figure DEST_PATH_IMAGE503
的取值区间
Figure DEST_PATH_IMAGE504
以及斜距变化率
Figure DEST_PATH_IMAGE505
的取值区间
Figure DEST_PATH_IMAGE506
在斜距
Figure DEST_PATH_IMAGE507
与斜距变化率
Figure DEST_PATH_IMAGE508
平面上所圈定的区域为
Figure DEST_PATH_IMAGE509
，则各观测弧段对应目标轨道在斜距与斜距变化率平面上的容许域范围为区域
Figure DEST_PATH_IMAGE510
、区域
Figure DEST_PATH_IMAGE511
与区域
Figure DEST_PATH_IMAGE512
的交集，即：
Figure DEST_PATH_IMAGE513
（20）
其中，
Figure DEST_PATH_IMAGE514
为观测弧段对应目标轨道在斜距与斜距变化率平面上的容许域范围。
步骤4，在容许域内对两观测弧段间最小马氏距离进行优化。具体地，对待关联观测弧段两两在步骤3中划定的容许域范围内对斜距与斜距变化率组合进行优化，结合运用航天器轨道预报与偏差演化算法，找到使赤经、赤纬预报值与实际观测值马氏距离最小的斜距与斜距变化率组合，并对两观测弧段间最小马氏距离进行记录。其具体实施过程为：
步骤4.1，对待关联观测弧段两两在步骤3中划定的容许域范围
Figure DEST_PATH_IMAGE515
内对斜距与斜距变化率组合
Figure DEST_PATH_IMAGE516
进行优化，找到使得赤经赤纬预报值与实际观测值马氏距离
Figure DEST_PATH_IMAGE517
最小的斜距与斜距变化率组合
Figure 678259DEST_PATH_IMAGE516
。此处直接调用Matlab自带优化工具箱（Optimization Tool）中的fmincon函数以实现对该问题的优化，此处不再对优化算法的实现过程进行详述。
步骤4.2，由斜距与斜距变化率组合
Figure 493768DEST_PATH_IMAGE516
计算得到优化指标马氏距离
Figure 188186DEST_PATH_IMAGE517
，其计算步骤为：
步骤4.2.1，根据斜距
Figure DEST_PATH_IMAGE518
、斜距变化率
Figure DEST_PATH_IMAGE519
、中间时刻赤经
Figure DEST_PATH_IMAGE520
、中间时刻赤纬
Figure DEST_PATH_IMAGE521
及其变化率
Figure DEST_PATH_IMAGE522
、
Figure DEST_PATH_IMAGE523
计算观测弧段对应轨道状态。设选取的两个待关联观测弧段分别为E与F，两弧段对应容许域分别为
Figure DEST_PATH_IMAGE524
与
Figure DEST_PATH_IMAGE525
，对于容许域
Figure DEST_PATH_IMAGE526
中所挑选出的一组斜距与斜距变化率
Figure 763655DEST_PATH_IMAGE516
，可以计算出观测弧段E在弧段中间时刻
Figure DEST_PATH_IMAGE527
所对应的一组轨道状态
Figure DEST_PATH_IMAGE528
，计算式如下：
Figure DEST_PATH_IMAGE529
（21）
步骤4.2.2，构建观测弧段E在
Figure DEST_PATH_IMAGE530
时刻对应轨道状态在当地轨道坐标系下的轨道状态协方差矩阵
Figure DEST_PATH_IMAGE531
。通过对整个观测弧段的数据点进行多项式拟合所得到的中间时刻赤经
Figure DEST_PATH_IMAGE532
、中间时刻赤纬
Figure DEST_PATH_IMAGE533
及其变化率
Figure DEST_PATH_IMAGE534
、
Figure DEST_PATH_IMAGE535
的标准差可根据原始数据单点观测标准差进行估算，估算公式为：
Figure DEST_PATH_IMAGE536
（22）
其中，
Figure DEST_PATH_IMAGE537
与
Figure DEST_PATH_IMAGE538
分别为原始数据赤经赤纬单点观测标准差，
Figure DEST_PATH_IMAGE539
是观测弧段的数据点数量，
Figure DEST_PATH_IMAGE540
是观测弧段首尾数据点横跨的时间长度。则观测弧段E对应轨道状态在观测空间内的协方差矩阵
Figure DEST_PATH_IMAGE541
可表示为：
Figure DEST_PATH_IMAGE542
（23）
则
Figure DEST_PATH_IMAGE543
可通过
Figure DEST_PATH_IMAGE544
由下式计算得到，为：
Figure DEST_PATH_IMAGE545
（24）
其中，
Figure DEST_PATH_IMAGE546
与
Figure DEST_PATH_IMAGE547
分别为观测空间到地心惯性系与地心惯性系到当地轨道坐标系的转换矩阵。
Figure DEST_PATH_IMAGE548
的计算式为：
Figure DEST_PATH_IMAGE549
（25）
式中
Figure DEST_PATH_IMAGE550
、
Figure DEST_PATH_IMAGE551
、
Figure DEST_PATH_IMAGE552
和
Figure DEST_PATH_IMAGE553
的定义为：
Figure DEST_PATH_IMAGE554
（26）
Figure DEST_PATH_IMAGE555
的计算式为：
Figure DEST_PATH_IMAGE556
（27）
其中，变量
Figure DEST_PATH_IMAGE557
到
Figure DEST_PATH_IMAGE558
的定义如下：
Figure 72801DEST_PATH_IMAGE244
（28）
步骤4.2.3，利用航天器轨道预报与偏差演化算法将观测弧段E在
Figure DEST_PATH_IMAGE559
时刻对应轨道状态
Figure DEST_PATH_IMAGE560
与轨道状态协方差矩阵
Figure DEST_PATH_IMAGE561
预报至观测弧段F对应弧段中间时刻
Figure DEST_PATH_IMAGE562
得到得到预报轨道状态
Figure DEST_PATH_IMAGE563
与预报轨道状态协方差矩阵
Figure DEST_PATH_IMAGE564
。航天器轨道预报与偏差演化算法是航天领域的成熟算法，有多重基于不同模型的算法，此处采用更为贴合实际的非线性轨道预报与偏差演化算法可提升最终关联聚类精度，有关非线性轨道预报与偏差演化算法的详细内容可以参考以下文献：杨震.非线性轨道机动瞄准与偏差演化分析方法[D]. 长沙：国防科技大学研究生院博士学位论文，2018,04.
步骤4.2.4，将经预报后得到的
Figure DEST_PATH_IMAGE565
与
Figure DEST_PATH_IMAGE566
重新转换至观测空间，得到时刻赤经赤纬的预测值
Figure DEST_PATH_IMAGE567
、
Figure DEST_PATH_IMAGE568
以及观测空间内的预报协方差矩阵
Figure DEST_PATH_IMAGE569
。赤经赤纬预测值的计算式如下：
Figure DEST_PATH_IMAGE570
（29）
预报协方差矩阵
Figure DEST_PATH_IMAGE571
可通过
Figure DEST_PATH_IMAGE572
由下式计算得到：
Figure DEST_PATH_IMAGE573
（30）
其中，
Figure DEST_PATH_IMAGE574
与
Figure DEST_PATH_IMAGE575
分别为当地轨道坐标系到地心惯性系与地心惯性系到观测空间的转换矩阵。
Figure DEST_PATH_IMAGE576
的计算式为：
Figure DEST_PATH_IMAGE577
（31）
其中，变量
Figure DEST_PATH_IMAGE578
到
Figure DEST_PATH_IMAGE579
的定义如下：
Figure DEST_PATH_IMAGE580
（32）
Figure DEST_PATH_IMAGE581
的计算式为：
Figure DEST_PATH_IMAGE582
（33）
式中变量
Figure DEST_PATH_IMAGE583
、
Figure DEST_PATH_IMAGE584
、
Figure DEST_PATH_IMAGE585
的定义如下：
Figure 526390DEST_PATH_IMAGE273
（34）
需要注意的是，由于计算马氏距离时仅用到了赤经赤纬而没有涉及赤经赤纬变化率，步骤4.2.4中计算得到的预报协方差矩阵
Figure DEST_PATH_IMAGE586
为
Figure DEST_PATH_IMAGE587
矩阵。
步骤4.2.5，计算由观测弧段E预报得到的时刻赤经赤纬预测值
Figure DEST_PATH_IMAGE588
、
Figure DEST_PATH_IMAGE589
与观测弧段F在
Figure DEST_PATH_IMAGE590
时刻赤经赤纬拟合值
Figure DEST_PATH_IMAGE591
、
Figure DEST_PATH_IMAGE592
的马氏距离
Figure DEST_PATH_IMAGE593
。马氏距离
Figure DEST_PATH_IMAGE594
计算式如下：
Figure DEST_PATH_IMAGE595
（35）
其中，上标T表示矩阵转置。计算得到的马氏距离
Figure DEST_PATH_IMAGE596
是工程中一种被广泛用作评定数据之间的相似度的无量纲指标，因此此处不对马氏距离进行详细介绍。
步骤4.3，将所有待关联观测弧段两两之间按步骤4.2所述计算马氏距离
Figure 973683DEST_PATH_IMAGE596
，并通过优化到两待关联观测弧段间最小马氏距离
Figure DEST_PATH_IMAGE597
，将每组观测弧段间的最小马氏距离进行记录并保存。
步骤5，根据两观测弧段间最小马氏距离判别两者是否关联。具体地，以步骤4中所记录的两观测弧段间最小马氏距离为关联判别依据，得到观测弧段两两关联匹配结果。其具体实施过程为：
以步骤4中所记录的两观测弧段间最小马氏距离为关联判别依据，一一进行判别，得到观测弧段两两关联匹配结果。可采用工程上常用的马氏距离判别依据进行判别，即
Figure DEST_PATH_IMAGE598
（36）
若最小马氏距离小于等于5，则认为这两个观测弧段之间关联成功，可能为对同一空间目标进行观测所产生的观测弧段。
步骤6，根据步骤5中所得观测弧段两两关联匹配结果，构建观测弧段关联矩阵，并利用BEA算法对观测弧段关联矩阵进行行列变换，将观测弧段关联矩阵转换成观测弧段聚类矩阵。其具体实施过程为：
步骤6.1，构建观测弧段关联矩阵
Figure DEST_PATH_IMAGE599
。关联矩阵
Figure DEST_PATH_IMAGE600
中第
Figure DEST_PATH_IMAGE601
行第
Figure DEST_PATH_IMAGE602
列元素
Figure DEST_PATH_IMAGE603
按照如下规则进行取值：
Figure DEST_PATH_IMAGE604
（37）
其中，
Figure DEST_PATH_IMAGE605
表示第
Figure DEST_PATH_IMAGE606
个待关联弧段与第
Figure DEST_PATH_IMAGE607
个待关联弧段间的最小马氏距离；
步骤6.2，利用BEA（Bond Energy Algorithm）算法对观测弧段关联矩阵
Figure DEST_PATH_IMAGE608
进行行列变换，将观测弧段关联矩阵
Figure 431471DEST_PATH_IMAGE608
变换成观测弧段聚类矩阵
Figure DEST_PATH_IMAGE609
。BEA算法是一种广泛应用在分布式数据库系统中大型表的纵向划分的算法，还可以实现矩阵元素的聚类。有关BEA算法的原理与具体实现步骤，可以参考以下文献：Ozsu M T, Valduriez P.Principles of distributed database systems[M]. [S.l.]：Prentice-Hall，1999.
经变换后，聚类矩阵
Figure 749451DEST_PATH_IMAGE609
中行列序号与原观测弧段序号不在一一对应，而是对应变换为：58 56 50 46 51 49 57 55 53 45 48 47 52 54 36 34 32 44 42 35 33 43 4039 37 31 41 38 30 28 26 22 27 25 29 20 21 19 18 17 16 23 24 11 13 7 5 4 3 121 14 8 6 10 2 15 9。
步骤7，根据观测弧段聚类矩阵元素排列特征进行分割，得到最终关联聚类结果。其具体实施过程为：
步骤7.1，构建聚类分割辅助序列
Figure DEST_PATH_IMAGE610
与
Figure DEST_PATH_IMAGE611
。为了实现对聚类矩阵
Figure 47709DEST_PATH_IMAGE609
的分割，首先定义两个各具有58个元素的序列
Figure 965986DEST_PATH_IMAGE610
与
Figure 883258DEST_PATH_IMAGE611
，序列
Figure DEST_PATH_IMAGE612
与
Figure 570591DEST_PATH_IMAGE611
中元素的取值与聚类矩阵
Figure 898804DEST_PATH_IMAGE609
中的元素有关，存在如下关系：
Figure DEST_PATH_IMAGE613
（38）
其中，
Figure DEST_PATH_IMAGE614
表示聚类矩阵
Figure DEST_PATH_IMAGE615
中的第
Figure DEST_PATH_IMAGE616
行第
Figure DEST_PATH_IMAGE617
列元素，
Figure DEST_PATH_IMAGE618
表示序列
Figure DEST_PATH_IMAGE619
中的第
Figure DEST_PATH_IMAGE620
个元素，
Figure DEST_PATH_IMAGE621
表示序列
Figure DEST_PATH_IMAGE622
中的第
Figure DEST_PATH_IMAGE623
个元素；
步骤7.2，根据序列
Figure DEST_PATH_IMAGE624
与
Figure DEST_PATH_IMAGE625
中元素变化规律对聚类矩阵
Figure DEST_PATH_IMAGE626
进行分割。当
Figure DEST_PATH_IMAGE627
与
Figure DEST_PATH_IMAGE628
中元素变化规律满足以下条件时对聚类矩阵
Figure DEST_PATH_IMAGE629
进行分割：
Figure DEST_PATH_IMAGE630
（39）
满足以上条件的
Figure DEST_PATH_IMAGE631
值为分割点。经计算，得到3个分割点，分别为
Figure DEST_PATH_IMAGE632
、
Figure DEST_PATH_IMAGE633
与
Figure DEST_PATH_IMAGE634
。
步骤7.3，根据3个分割点的取值
Figure DEST_PATH_IMAGE635
、
Figure DEST_PATH_IMAGE636
与
Figure DEST_PATH_IMAGE637
可知，
Figure DEST_PATH_IMAGE638
被分割为4个聚类子矩阵，分别为第
Figure DEST_PATH_IMAGE639
行第
Figure DEST_PATH_IMAGE640
列元素构成的聚类子矩阵、第
Figure DEST_PATH_IMAGE641
行第
Figure DEST_PATH_IMAGE642
列元素构成的聚类子矩阵、第
Figure DEST_PATH_IMAGE643
行第
Figure DEST_PATH_IMAGE644
列元素构成的聚类子矩阵与第
Figure DEST_PATH_IMAGE645
行第
Figure DEST_PATH_IMAGE646
列元素构成的聚类子矩阵。对于位于同一聚类子矩阵内的观测弧段，视为聚类成功，认为这些观测弧段是对同一空间目标进行观测所产生的观测弧段。
最终得到观测弧段关联聚类结果如下表所示：
表1实施例测试结果展示表
Figure DEST_PATH_IMAGE647
通过与测试答案进行对比可知，总共58个观测弧段，其中关联聚类正确个数为58，正确率为
Figure DEST_PATH_IMAGE648
，本方法可以实现对属于同一空间目标观测弧段的关联聚类。
以上所述仅为本发明的优选实施例，并非因此限制本发明的专利范围，凡是在本发明的发明构思下，利用本发明说明书及附图内容所作的等效结构变换，或直接/间接运用在其他相关的技术领域均包括在本发明的专利保护范围内。

Claims (10)

1.一种基于非线性偏差演化的天基光学观测短弧关联与聚类方法，其特征在于，包括如下步骤：
步骤1，获取多组分别属于不同空间目标的观测弧段，每个观测弧段包括若干个观测数据点，每个观测数据点由被观测目标相对于观测卫星的赤经、赤纬、观测时刻以及观测平台的位置速度信息组成；
步骤2，分别对每个观测弧段中赤经、赤纬关于时间的函数式进行拟合，得到赤经、赤纬随时间的变化率信息，并对存在明显异常的观测数据点进行剔除；
步骤3，根据各观测弧段中赤经、赤纬随时间的变化率信息与先验信息，对各观测弧段对应目标轨道在斜距与斜距变化率平面上的容许域范围进行划定；
步骤4，对待关联观测弧段两两在容许域范围内对斜距与斜距变化率组合进行优化，基于非线性偏差演化找到使赤经、赤纬预报值与实际观测值马氏距离最小的斜距与斜距变化率组合，并对两观测弧段间最小马氏距离进行记录；
步骤5，以两观测弧段间最小马氏距离为关联判别依据，得到观测弧段两两关联匹配结果；
步骤6，根据观测弧段两两关联匹配结果，构建观测弧段关联矩阵，并对观测弧段关联矩阵进行行列变换，将观测弧段关联矩阵转换成观测弧段聚类矩阵；
步骤7，根据观测弧段聚类矩阵行列元素排列特征对观测弧段聚类矩阵进行分割，得到观测弧段关联聚类结果，实现对属于同一空间目标的观测弧段的关联聚类。
2.根据权利要求1所述的基于非线性偏差演化的天基光学观测短弧关联与聚类方法，其特征在于，步骤1中，所述多组分别属于不同空间目标的观测弧段为
Figure DEST_PATH_IMAGE001
，
Figure DEST_PATH_IMAGE002
，
Figure DEST_PATH_IMAGE003
，其中，
Figure DEST_PATH_IMAGE004
为空间目标数量，
Figure DEST_PATH_IMAGE005
为观测弧段的数量；
Figure DEST_PATH_IMAGE006
为第
Figure DEST_PATH_IMAGE007
个空间目标的第
Figure DEST_PATH_IMAGE008
个观测弧段，
Figure DEST_PATH_IMAGE009
，其中，
Figure DEST_PATH_IMAGE010
为第
Figure 824124DEST_PATH_IMAGE007
个空间目标的第
Figure 210106DEST_PATH_IMAGE008
个观测弧段的数据行数，下标
Figure DEST_PATH_IMAGE011
表示观测弧段中第
Figure 927526DEST_PATH_IMAGE011
行数据，
Figure DEST_PATH_IMAGE012
为观测历元时刻，
Figure DEST_PATH_IMAGE013
为赤经，
Figure DEST_PATH_IMAGE014
为赤纬，
Figure DEST_PATH_IMAGE015
与
Figure DEST_PATH_IMAGE016
为分别为每行数据观测历元时刻对应的观测卫星的位置和速度矢量。
3.根据权利要求2所述的基于非线性偏差演化的天基光学观测短弧关联与聚类方法，其特征在于，步骤2中，采用二次多项式分别对每个观测弧段中赤经、赤纬关于时间的函数式进行拟合，具体为：
设赤经
Figure DEST_PATH_IMAGE017
和赤纬
Figure DEST_PATH_IMAGE018
对时间的函数
Figure DEST_PATH_IMAGE019
、
Figure DEST_PATH_IMAGE020
分别表示为：
Figure DEST_PATH_IMAGE021
其中，
Figure DEST_PATH_IMAGE022
、
Figure DEST_PATH_IMAGE023
、
Figure DEST_PATH_IMAGE024
、
Figure DEST_PATH_IMAGE025
、
Figure DEST_PATH_IMAGE026
、
Figure DEST_PATH_IMAGE027
为多项式待定系数，各待定系数的初值取为：
Figure DEST_PATH_IMAGE028
Figure DEST_PATH_IMAGE029
对
Figure DEST_PATH_IMAGE030
、
Figure DEST_PATH_IMAGE031
、
Figure DEST_PATH_IMAGE032
的偏导数以及
Figure DEST_PATH_IMAGE033
对
Figure DEST_PATH_IMAGE034
、
Figure DEST_PATH_IMAGE035
、
Figure DEST_PATH_IMAGE036
的偏导数分别为：
Figure DEST_PATH_IMAGE037
使用最小二乘法得到对
Figure DEST_PATH_IMAGE038
、
Figure DEST_PATH_IMAGE039
、
Figure DEST_PATH_IMAGE040
初值的改进量
Figure DEST_PATH_IMAGE041
、
Figure DEST_PATH_IMAGE042
、
Figure DEST_PATH_IMAGE043
，以及为对
Figure DEST_PATH_IMAGE044
、
Figure DEST_PATH_IMAGE045
、
Figure DEST_PATH_IMAGE046
初值的改进量
Figure DEST_PATH_IMAGE047
、
Figure DEST_PATH_IMAGE048
、
Figure DEST_PATH_IMAGE049
：
Figure DEST_PATH_IMAGE050
其中，
Figure DEST_PATH_IMAGE051
是
Figure DEST_PATH_IMAGE052
的矩阵，
Figure DEST_PATH_IMAGE053
为
Figure DEST_PATH_IMAGE054
的转置矩阵，上标-1表示矩阵的求逆运算，
Figure DEST_PATH_IMAGE055
是
Figure DEST_PATH_IMAGE056
维的向量，
Figure DEST_PATH_IMAGE057
为赤经的多项式预测值，
Figure DEST_PATH_IMAGE058
是
Figure DEST_PATH_IMAGE059
维的向量，
Figure DEST_PATH_IMAGE060
为赤纬的多项式预测值；
则将
Figure DEST_PATH_IMAGE061
、
Figure DEST_PATH_IMAGE062
、
Figure DEST_PATH_IMAGE063
以及
Figure DEST_PATH_IMAGE064
、
Figure DEST_PATH_IMAGE065
、
Figure DEST_PATH_IMAGE066
更新为：
Figure DEST_PATH_IMAGE067
重复最小二乘法计算以及
Figure DEST_PATH_IMAGE068
、
Figure DEST_PATH_IMAGE069
、
Figure DEST_PATH_IMAGE070
、
Figure DEST_PATH_IMAGE071
、
Figure DEST_PATH_IMAGE072
、
Figure DEST_PATH_IMAGE073
的更新过程直至
Figure DEST_PATH_IMAGE074
、
Figure DEST_PATH_IMAGE075
小于设定的阈值，即得到拟合出的
Figure DEST_PATH_IMAGE076
、
Figure DEST_PATH_IMAGE077
、
Figure DEST_PATH_IMAGE078
以及
Figure DEST_PATH_IMAGE079
、
Figure DEST_PATH_IMAGE080
、
Figure DEST_PATH_IMAGE081
。
4.根据权利要求2所述的基于非线性偏差演化的天基光学观测短弧关联与聚类方法，其特征在于，步骤2中，所述对存在明显异常的观测数据点进行剔除，具体为：
对于每个观测时刻的观测数据点，得到其对应时刻的赤经、赤纬拟合值，将对应时刻的赤经、赤纬拟合值与真实观测值作差得到赤经、赤纬的残差，计算一个观测弧段的拟合值与实际观测值残差的标准差
Figure DEST_PATH_IMAGE082
，为：
Figure DEST_PATH_IMAGE083
其中，
Figure DEST_PATH_IMAGE084
表示第
Figure DEST_PATH_IMAGE085
个观测数据的残差，
Figure DEST_PATH_IMAGE086
为残差均值，
Figure DEST_PATH_IMAGE087
为观测数据点个数；
将观测弧段中残差大于
Figure DEST_PATH_IMAGE088
的观测数据点剔除。
5.根据权利要求1至4任一项所述的基于非线性偏差演化的天基光学观测短弧关联与聚类方法，其特征在于，步骤3具体包括：
估算被观测目标半长轴的取值区间
Figure DEST_PATH_IMAGE089
、偏心率的取值区间
Figure DEST_PATH_IMAGE090
、相对于观测卫星的斜距
Figure DEST_PATH_IMAGE091
的取值区间
Figure DEST_PATH_IMAGE092
以及斜距变化率
Figure DEST_PATH_IMAGE093
的取值区间
Figure DEST_PATH_IMAGE094
；
根据观测弧段的中间时刻赤经
Figure DEST_PATH_IMAGE095
、中间时刻赤纬
Figure DEST_PATH_IMAGE096
及其变化率信息
Figure DEST_PATH_IMAGE097
、
Figure DEST_PATH_IMAGE098
以及被观测目标半长轴的取值区间
Figure DEST_PATH_IMAGE099
，划定观测弧段对应被观测目标在斜距
Figure DEST_PATH_IMAGE100
与斜距变化率
Figure DEST_PATH_IMAGE101
平面上的区域
Figure DEST_PATH_IMAGE102
；
根据观测弧段的中间时刻赤经
Figure DEST_PATH_IMAGE103
、中间时刻赤纬
Figure DEST_PATH_IMAGE104
及其变化率信息
Figure DEST_PATH_IMAGE105
、
Figure DEST_PATH_IMAGE106
以及被观测目标偏心率的取值区间
Figure DEST_PATH_IMAGE107
，划定观测弧段对应被观测目标在斜距
Figure DEST_PATH_IMAGE108
与斜距变化率
Figure DEST_PATH_IMAGE109
平面上的区域
Figure DEST_PATH_IMAGE110
；
设斜距
Figure DEST_PATH_IMAGE111
的取值区间
Figure DEST_PATH_IMAGE112
以及斜距变化率
Figure DEST_PATH_IMAGE113
的取值区间
Figure DEST_PATH_IMAGE114
在斜距
Figure DEST_PATH_IMAGE115
与斜距变化率
Figure DEST_PATH_IMAGE116
平面上所圈定的区域为
Figure DEST_PATH_IMAGE117
，则各观测弧段对应目标轨道在斜距与斜距变化率平面上的容许域范围为区域
Figure DEST_PATH_IMAGE118
、区域
Figure DEST_PATH_IMAGE119
与区域
Figure DEST_PATH_IMAGE120
的交集，即：
Figure DEST_PATH_IMAGE121
其中，
Figure DEST_PATH_IMAGE122
为观测弧段对应目标轨道在斜距与斜距变化率平面上的容许域范围。
6.根据权利要求1至4任一项所述的基于非线性偏差演化的天基光学观测短弧关联与聚类方法，其特征在于，步骤4中，根据两观测弧段间最小马氏距离的确定过程为：
设选取的两个待关联观测弧段分别为E与F，两观测弧段对应的容许域范围分别为
Figure DEST_PATH_IMAGE123
与
Figure DEST_PATH_IMAGE124
；
对于容许域范围
Figure DEST_PATH_IMAGE125
中选出的一组斜距与斜距变化率
Figure DEST_PATH_IMAGE126
，根据观测弧段E的中间时刻赤经
Figure DEST_PATH_IMAGE127
、中间时刻赤纬
Figure DEST_PATH_IMAGE128
及其变化率
Figure DEST_PATH_IMAGE129
、
Figure DEST_PATH_IMAGE130
，计算观测弧段E在弧段中间时刻
Figure DEST_PATH_IMAGE131
所对应的轨道状态
Figure DEST_PATH_IMAGE132
；
构建观测弧段E在
Figure DEST_PATH_IMAGE133
时刻对应轨道状态在当地轨道坐标系下的轨道状态协方差矩阵
Figure DEST_PATH_IMAGE134
；
将观测弧段E在
Figure DEST_PATH_IMAGE135
时刻对应的轨道状态
Figure DEST_PATH_IMAGE136
与轨道状态协方差矩阵
Figure DEST_PATH_IMAGE137
预报至观测弧段F对应弧段中间时刻
Figure DEST_PATH_IMAGE138
，得到预报轨道状态
Figure DEST_PATH_IMAGE139
与预报轨道状态协方差矩阵
Figure DEST_PATH_IMAGE140
；
将经预报后得到的
Figure DEST_PATH_IMAGE141
与
Figure DEST_PATH_IMAGE142
转换至观测空间，得到
Figure DEST_PATH_IMAGE143
时刻赤经赤纬的预测值
Figure DEST_PATH_IMAGE144
、
Figure DEST_PATH_IMAGE145
以及观测空间内的预报协方差矩阵
Figure DEST_PATH_IMAGE146
；
计算由观测弧段E预报得到的
Figure DEST_PATH_IMAGE147
时刻赤经赤纬预测值
Figure DEST_PATH_IMAGE148
、
Figure DEST_PATH_IMAGE149
与观测弧段F在
Figure DEST_PATH_IMAGE150
时刻赤经赤纬拟合值
Figure DEST_PATH_IMAGE151
、
Figure DEST_PATH_IMAGE152
的马氏距离
Figure DEST_PATH_IMAGE153
，为：
Figure DEST_PATH_IMAGE154
通过对斜距与斜距变化率组合进行优化，即能得到两待关联观测弧段E与F间的最小马氏距离
Figure DEST_PATH_IMAGE155
。
7.根据权利要求1至4任一项所述的基于非线性偏差演化的天基光学观测短弧关联与聚类方法，其特征在于，步骤5具体为：
若两待关联观测弧段之间的最小马氏距离小于或等于5，则判定两观测弧段之间关联成功。
8.根据权利要求7所述的基于非线性偏差演化的天基光学观测短弧关联与聚类方法，其特征在于，步骤6中，所述观测弧段关联矩阵为
Figure DEST_PATH_IMAGE156
，
Figure DEST_PATH_IMAGE157
为待关联弧段数量；
关联矩阵
Figure DEST_PATH_IMAGE158
中第
Figure DEST_PATH_IMAGE159
行第
Figure DEST_PATH_IMAGE160
列元素
Figure DEST_PATH_IMAGE161
的取值为：
Figure DEST_PATH_IMAGE162
其中，
Figure DEST_PATH_IMAGE163
表示第
Figure DEST_PATH_IMAGE164
个待关联弧段与第
Figure DEST_PATH_IMAGE165
个待关联弧段间的最小马氏距离；
利用BEA算法对观测弧段关联矩阵
Figure DEST_PATH_IMAGE166
进行行列变换，将观测弧段关联矩阵
Figure DEST_PATH_IMAGE167
变换成观测弧段聚类矩阵
Figure DEST_PATH_IMAGE168
。
9.根据权利要求8所述的基于非线性偏差演化的天基光学观测短弧关联与聚类方法，其特征在于，步骤6中，利用BEA算法对观测弧段关联矩阵
Figure DEST_PATH_IMAGE169
进行行列变换，将观测弧段关联矩阵
Figure DEST_PATH_IMAGE170
变换成观测弧段聚类矩阵
Figure DEST_PATH_IMAGE171
。
10.根据权利要求9所述的基于非线性偏差演化的天基光学观测短弧关联与聚类方法，其特征在于，步骤7具体包括：
构建聚类分割辅助序列
Figure DEST_PATH_IMAGE172
与
Figure DEST_PATH_IMAGE173
，序列
Figure DEST_PATH_IMAGE174
与
Figure DEST_PATH_IMAGE175
中元素的取值为：
Figure DEST_PATH_IMAGE176
其中，
Figure DEST_PATH_IMAGE177
表示观测弧段聚类矩阵
Figure DEST_PATH_IMAGE178
中的第
Figure DEST_PATH_IMAGE179
行第
Figure DEST_PATH_IMAGE180
列元素，
Figure DEST_PATH_IMAGE181
表示序列
Figure DEST_PATH_IMAGE182
中的第
Figure DEST_PATH_IMAGE183
个元素，
Figure DEST_PATH_IMAGE184
表示序列
Figure DEST_PATH_IMAGE185
中的第
Figure DEST_PATH_IMAGE186
个元素；
以
Figure DEST_PATH_IMAGE187
与
Figure DEST_PATH_IMAGE188
中满足
Figure DEST_PATH_IMAGE189
的
Figure DEST_PATH_IMAGE190
为分割点，将观测弧段聚类矩阵
Figure DEST_PATH_IMAGE191
分割为若干聚类子矩阵；
判定同一聚类子矩阵内的观测弧段是对同一空间目标进行观测所产生的观测弧段，由此得到观测弧段关联聚类结果。
CN202211594600.6A 2022-12-13 2022-12-13 基于非线性偏差演化的天基光学观测短弧关联与聚类方法 Active CN115659196B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211594600.6A CN115659196B (zh) 2022-12-13 2022-12-13 基于非线性偏差演化的天基光学观测短弧关联与聚类方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211594600.6A CN115659196B (zh) 2022-12-13 2022-12-13 基于非线性偏差演化的天基光学观测短弧关联与聚类方法

Publications (2)

Publication Number Publication Date
CN115659196A true CN115659196A (zh) 2023-01-31
CN115659196B CN115659196B (zh) 2023-06-23

Family

ID=85019537

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211594600.6A Active CN115659196B (zh) 2022-12-13 2022-12-13 基于非线性偏差演化的天基光学观测短弧关联与聚类方法

Country Status (1)

Country Link
CN (1) CN115659196B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117249835A (zh) * 2023-11-09 2023-12-19 南京航空航天大学 一种天基无源协同多目标观测数据关联定位方法

Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105224737A (zh) * 2015-09-22 2016-01-06 中国人民解放军63921部队 一种空间目标轨道改进初值修正方法
CN107153209A (zh) * 2017-07-06 2017-09-12 武汉大学 一种短弧段低轨导航卫星实时精密定轨方法
CN108761507A (zh) * 2018-05-21 2018-11-06 中国人民解放军战略支援部队信息工程大学 基于短弧定轨和预报的导航卫星轨道快速恢复方法
CN109506630A (zh) * 2018-11-02 2019-03-22 北京空间飞行器总体设计部 一种甚短弧高频仅角度观测值的初轨确定方法
CN110017832A (zh) * 2019-03-19 2019-07-16 华中科技大学 一种基于Gauss解群优选的短弧初轨确定方法
CN110986962A (zh) * 2019-12-09 2020-04-10 中国科学院国家授时中心 基于高轨通信卫星的低轨卫星全弧段测定轨方法
CN111578950A (zh) * 2020-06-09 2020-08-25 中国人民解放军63921部队 一种面向天基光学监视的geo目标自主弧段关联与定轨方法
CN112945182A (zh) * 2021-01-26 2021-06-11 深圳市微视星辰数据网络科技有限公司 一种观测数据-编目目标关联匹配方法
CN113204917A (zh) * 2021-04-25 2021-08-03 中国科学院国家空间科学中心 针对geo目标的天基光学测角弧段初定轨及关联方法
CN114396953A (zh) * 2021-11-30 2022-04-26 中国西安卫星测控中心 一种天基短弧光学测轨数据的关联方法

Patent Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105224737A (zh) * 2015-09-22 2016-01-06 中国人民解放军63921部队 一种空间目标轨道改进初值修正方法
CN107153209A (zh) * 2017-07-06 2017-09-12 武汉大学 一种短弧段低轨导航卫星实时精密定轨方法
CN108761507A (zh) * 2018-05-21 2018-11-06 中国人民解放军战略支援部队信息工程大学 基于短弧定轨和预报的导航卫星轨道快速恢复方法
CN109506630A (zh) * 2018-11-02 2019-03-22 北京空间飞行器总体设计部 一种甚短弧高频仅角度观测值的初轨确定方法
CN110017832A (zh) * 2019-03-19 2019-07-16 华中科技大学 一种基于Gauss解群优选的短弧初轨确定方法
US20210199439A1 (en) * 2019-03-19 2021-07-01 Huazhong University Of Science And Technology Short arc initial orbit determining method based on gauss solution cluster
CN110986962A (zh) * 2019-12-09 2020-04-10 中国科学院国家授时中心 基于高轨通信卫星的低轨卫星全弧段测定轨方法
CN111578950A (zh) * 2020-06-09 2020-08-25 中国人民解放军63921部队 一种面向天基光学监视的geo目标自主弧段关联与定轨方法
CN112945182A (zh) * 2021-01-26 2021-06-11 深圳市微视星辰数据网络科技有限公司 一种观测数据-编目目标关联匹配方法
CN113204917A (zh) * 2021-04-25 2021-08-03 中国科学院国家空间科学中心 针对geo目标的天基光学测角弧段初定轨及关联方法
CN114396953A (zh) * 2021-11-30 2022-04-26 中国西安卫星测控中心 一种天基短弧光学测轨数据的关联方法

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117249835A (zh) * 2023-11-09 2023-12-19 南京航空航天大学 一种天基无源协同多目标观测数据关联定位方法
CN117249835B (zh) * 2023-11-09 2024-03-29 南京航空航天大学 一种天基无源协同多目标观测数据关联定位方法

Also Published As

Publication number Publication date
CN115659196B (zh) 2023-06-23

Similar Documents

Publication Publication Date Title
Schmidt Application of state-space methods to navigation problems
Viegas et al. Discrete-time distributed Kalman filter design for formations of autonomous vehicles
US8260732B2 (en) Method for identifying Hammerstein models
CN103439731A (zh) 基于无迹卡尔曼滤波的gps/ins组合导航方法
US10768312B1 (en) Integrity analysis method based on kinematic-to-kinematic relative positioning scenario
CN104317990A (zh) 一种基于风险的多阶段任务航天器可靠性改进方法
CN115659196A (zh) 基于非线性偏差演化的天基光学观测短弧关联与聚类方法
CN106767780A (zh) 基于Chebyshev多项式插值逼近的扩展椭球集员滤波方法
Quine A derivative-free implementation of the extended Kalman filter
CN105629272A (zh) 短弧段批处理卫星自主定轨方法及装置
CN115523927B (zh) 基于光学传感器观测的geo航天器机动检测方法
EP3176608B1 (en) Efficient covariance matrix update
Sarychev et al. Optimal regressors search subjected to vector autoregression of unevenly spaced TLE series
Lei et al. Identification of uncatalogued LEO space objects by a ground-based EO array
Girija et al. Tracking filter and multi-sensor data fusion
CN114396953A (zh) 一种天基短弧光学测轨数据的关联方法
Han et al. Rapid algorithm for covariance ellipsoid model based collision warning of space objects
CN110703284B (zh) 一种基于稀疏核学习的单站gnss瞬时速度和加速度构建方法
Jauberthie et al. An optimal input design procedure
Hall Implementation recommendations and usage boundaries for the two-dimensional probability of collision calculation
Goh et al. Real-time estimation of satellite's two-line elements via positioning data
Schroeder et al. Removal of ambiguous wind directions for a Ku-band wind scatterometer using three different azimuth angles
Li et al. Sequential unbiased converted measurement non‐linear filter with range rate in ECEF coordinates
Feldhacker et al. Multi-element trajectory models for satellite tour missions
Gehly et al. Lp-norm batch estimation as applied to orbit determination

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant