1 引言

空间目标状态异常主要包括下述三种情况:1)目标轨道异常,即目标从其稳定运行轨道变换至临时运行轨道;2)目标失效,即目标出现翻滚、旋转等不受控运动;3)目标工作状态发生变化,即目标姿态显著不同于以往观测时的姿态。其中,目标失效以及目标工作状态发生变化这两种情况可以利用光学观测数据进行检测。检测原则可总结为从当前的观测数据曲线中找出与历史观测数据曲线不匹配的,将其判定为异常。为了确保历史观测数据能够较完整地呈现出空间目标的运动规律,需要庞大的历史数据库作为支撑。

近年来,随着我国广角望远镜阵观测数据的指数级增长以及仿真模拟技术的逐渐完善,空间目标光学观测系统已经积累了大量不同目标的历史观测数据^[1]。传统的异常检测方法主要是基于特定的物理模型,根据各模型对当前状态下卫星光度特性的描述效果进行分类。但这种方法成本高,对目标的先验信息要求苛刻,导致处理能力有限。随着监视目标数量的不断增加,传统的检测方法已难以适用于整个卫星群。目前,国内外大多数检测方法以数据驱动的方式,通过分析观测数据的特征对空间目标的异常状态进行判断与分类。这些方法不需要精确地反演出目标的具体工作状态,无需建立具体的物理模型来描述目标的运动特征,因而得到了广泛关注。

Linares等^[2]提出了基于深度卷积神经网络(CNN)的空间目标异常分类方法。该方法基于数据驱动,从光度曲线中提取特征,并结合姿态动力学和光学散射特性的约束,在曲线库中按设定标签对空间目标进行异常分类。仿真实验表明该方法对空间目标的异常分类准确率高达99.6%。随后,Linares等^[3]基于实测数据对该方法进行了验证,结果显示,异常分类准确率达到了75%。波尔航天公司的Regan^[4]将模块化的神经网络用于空间目标监视中,着重解决目标发现、目标分类及决策三方面的问题,这一方法在面对不同的任务类型和参数变化时体现出了良好的自适应性。

Xu等^[5]讨论了地球同步轨道(GEO)卫星表面材料镜面反射的振幅变化特性和光束宽度,并对正常工作状态下由太阳运动引起的卫星镜面反射持续时间进行了定量分析,提出了将卫星镜面反射持续时间长短作为GEO卫星工作状态是否异常的判断依据。然而,由于低轨道卫星的镜面反射持续时间很短,通常无法观测到,因此该方法难以应用于低轨卫星的异常检测。

现有基于光度数据的空间目标异常检测大都面向GEO轨道目标。由于GEO轨道目标相对地面静止,所以目标的位置变化对观测数据的影响较小。而对于相对地面运动速度较快的低轨目标来说,若不考虑光学几何关系而直接对历史观测数据进行特征提取,将造成特征无法提取或提取错误。因为实际上目标的观测数据是由目标的照射面和观测面以及入射角和观测角决定的,这一几何位置关系不仅与空间目标自身的工作状态有关,还与太阳、目标、测站的相对位置关系有关^[6]。所以直接采用数据驱动的方法进行大数据学习,不适用于低轨目标。

因此,本文提出了一种基于动态时间规整(DTW)距离的异常检测方法,用该方法来诊断空间目标的工作状态,解决低轨目标的异常检测问题。首先,利用与观测距离无关的光学散射截面积(OCS)来描述空间目标的光度特性;其次,基于OpenGL拾取技术^[7]建立不同观测场景、不同观测时间下的空间目标光度观测数据集;最后,根据太阳、目标、测站三者之间的相对位置关系,从历史观测数据集中找到与当前观测几何条件相似的数据后,将当前观测数据与历史数据进行比较。本文的异常判断依据为计算当前OCS曲线与历史OCS曲线的DTW距离,若其大于所设阈值,则判定为异常工作状态。

2 观测几何条件

2.1 观测几何条件的定义

观测几何条件可以概括为太阳、目标、测站之间的相对位置关系,通常情况下考虑为相位角,即空间目标指向太阳的向量与空间目标指向测站的向量之间的夹角β,β∈[0,π]。但经过探究发现存在相位角相同但观测几何条件不一致的情况,因此仅仅靠相位角一个参量来表示观测几何条件是不够的。本文以目标质心为原点,建立质心轨道坐标系,X轴指向目标速度方向,Z轴指向地心,Y轴满足右手定则。将相位角细分为太阳高低角θ_S、太阳方位角φ_S、测站高低角θ_F和测站方位角φ_F,其中:θ_S,θ_F∈[0,π];φ_S,φ_F∈[0,2π],如图2所示。

图 1. 观测几何条件示意图

Fig. 1. Schematic of observation geometry

下载图片 查看所有图片

2.2 观测几何条件的变化对目标OCS的影响

OCS为目标自身固有的属性,与目标形状、尺寸、表面材料以及运行状态有关^[8]。此外,在实际观测中,随着空间目标的运行,太阳、目标、测站之间的相对几何关系在不断改变,几何关系的不同会导致空间目标的受射面和测站的观测面不同,因而得到的OCS也是不同的。

对于低轨空间目标而言,空间目标相对测站的运动速度较快,观测几何关系变化较为剧烈,对观测结果的影响不可忽视。另外,低轨目标受到的地球非球形摄动力和大气摄动力较强,导致目标轨道在惯性空间发生了偏移和形变。再加上地球本身的自转,星下点轨迹会发生明显的偏移,导致观测窗口的长度发生明显变化,如图2所示,因而导致光度曲线在时域上发生偏移和形变。图3给出了观测几何条件不同时,目标OCS的变化曲线。从图3中可以看出,对于不同观测几何条件下的弧段,OCS曲线的差异是比较明显的。

图 2. 低轨目标观测。(a) 卫星星下点轨迹;(b) 不同窗口的观测几何条件随时间的变化

Fig. 2. Low-orbit objects observation. (a) Track of subsatellite point; (b) observation geometry of different windows varies with time

下载图片 查看所有图片

图 3. 不同观测几何条件下的空间目标OCS曲线。(a)三轴稳定;(b)自旋

Fig. 3. OCS curves of space objects at different observation geometry. (a) Three-axis stabilization; (b) spinning

下载图片 查看所有图片

因此,对于低轨目标来说,考虑观测几何条件对观测数据的影响尤为重要。对于目标是否异常的判断,需要用当前观测弧段的OCS数据与历史数据库中观测几何条件相似弧段内的OCS数据进行对比。

2.3 观测几何条件相似的弧段搜索

2.3.1 动态规整算法

DTW由日本学者Itakura在20世纪60年代提出,其整体思想就是通过动态规划算法(DP)对两曲线进行点对点的最优映射,以描述两条曲线的相似度。相对于欧氏距离,它的优势在于可以对时域上不等长的曲线之间的距离进行表示。在对未知曲线进行伸长或缩短变换,直到与参考曲线长度一致时,未知曲线的时间轴会发生扭曲或弯折,这样便于使其特征量与参考曲线对应^[9]。

对两条不等长的曲线进行等间隔采样,得到两个不等长的时间序列:A={a₁,a₂,…,a_i,…,a_m},长度为m;B={b₁,b₂,…,b_j,…,b_n},长度为n。首先,计算A和B上每两点之间的欧氏距离,公式为

d(Ai,Bj)=|ai-bj|,i∈[1,m],j∈[1,n]。(1)

d(A_i,B_j)也被称为局部距离,将A和B中所有点之间的局部距离对应填入一个m×n的矩阵D中。矩阵D被称为路径规整距离矩阵,其中D(i,j)=d(A_i,B_j)。定义规整路径函数W(k)=(w_a(k),w_b(k)),k表示第k次匹配,k∈[1,l],l为规整路径的长度,满足l∈[max(m,n),m+n-1],w_a(k)、w_b(k)分别为序列A、B中元素的下标,w_a(k)∈[1,m],w_b(k)∈[1,n]。W(k)表示序列A中第w_a(k)个元素和序列B中第w_b(k)个元素进行映射。通过DP算法在规整距离矩阵D中找到一条规整路径W来表示序列A和B之间的映射,使得序列A和B之间的局部距离之和最小。规整距离越小,A和B的相似度就越高。从D(1,1)出发沿着某条路径到达D(m,n),规整路径应满足的约束条件为边界性、单调性和连续性^[10]。

结合单调性和连续性的约束: wak≤wa(k+1)≤wa(k)+1wbk≤wb(k+1)≤wb(k)+1,可知每一个节点下一步有三个前进方向:(w_a(k+1),w_b(k+1))=(wa(k)+1,wb(k))(wa(k),wb(k)+1)(wa(k)+1,wb(k)+1),如图4中箭头所示。

图 4. 规整路径节点前进方向示意图

Fig. 4. Schematic of the forward direction of nodes in a warping path

下载图片 查看所有图片

要找出序列A和B之间的最优映射,应使规整路径最短,对累积距离C[w_a(k),w_b(k)]进行定义:在规整距离矩阵D中从(1,1)点开始对序列A和序列B进行局部距离的计算,到达点(w_a(k),w_b(k))时,共进行了k次计算,这k个局部距离之和就是累积距离C[w_a(k),w_b(k)];而矩阵D中到达某个点的路径不止一种,对应着不等的累积距离,当到达终点(m,n)后,最短累积距离 C-(m,n)所对应的路径就是规整路径。最短累积距离以下称为“DTW距离”。

由路径的可逆性可知,受到单调性和连续性的约束,D中任意一点(w_a(k),w_b(k))最多可由(w_a(k)-1,w_b(k))、(w_a(k),w_b(k)-1)和(w_a(k)-1,w_b(k)-1)三个点得到,则最短距离 C-[w_a(k),w_b(k)]可以表示为

C-[wa(k),wb(k)]=d[wa(k),wb(k)]+min{C[wa(k)-1,wb(k)],C[wa(k),wb(k)-1],C[wa(k)-1,wb(k)-1]},(2)

最后可通过DP算法求得最短总距离 C-(m,n),此即为序列A和B相似度的评判依据。

2.3.2 基于DTW的搜索算法

前文指出,空间目标观测几何条件变化规律由太阳高低角θ_S、太阳方位角φ_S、测站高低角θ_F和测站方位角φ_F的变化规律决定。因此搜索历史观测数据库中与当前观测数据具有相似观测几何条件变化规律的算法如下。

输入:①当前观测access(β_t、θ_St、θ_Ft、φ_St和φ_Ft),②历史观测数据集Dom 1(β、θ_S、θ_F、φ_S和φ_F)。

Step 1:For access in Dom 1

相位角DTW距离为β_DTW(access)=DTW[β_t,β(access)]。

Step 2:依据β_DTW,选取o₁个DTW距离最小的access作为Dom 2。

Step 3:For access in Dom 2

太阳高低角DTW距离为θ_SDTW(access)=DTW[θ_St,θ_S(access)],测站高低角DTW距离为θ_FDTW(access)=DTW[θ_Ft,θ_F(access)],太阳方位角DTW距离为φ_SDTW(access)=DTW[φ_St,φ_S(access)],测站方位角DTW距离为φ_FDTW(access)=DTW[φ_Ft,φ_F(access)]。

Step 4:将 Dom 2中的access分别依据θ_SDTW、θ_FDTW、φ_SDTW、φ_FDTW执行4种排序(DTW距离由小到大排列),分别作为矩阵A_mn的各列(m=o₁,n=4)。

Step 5: For i in range (1, o₁)

取前i行首次出现4次的access作为与当前access观测几何条件最相似的历史弧段bestaccess。

输出:bestaccess。

3 空间目标工作状态异常检测仿真研究

将测站设定在兴隆(25.49°E,101.17°N,2046.50 m),空间目标设定在半长轴为6785 km、倾角为20°的圆轨道上。定义三轴稳定为正常工作状态,且本体轴Z轴指向地心,Y轴(卫星帆板方向)沿轨道方向,X轴满足右手定则,如图5所示,一般情况下θ=0°。定义7种异常工作状态(type 1,X轴指向偏离30°;type 2,X轴指向偏离90°;type 3,X轴指向偏离130°;type 4,X轴对日定向;type 5,以0.12 r/min的角速度绕Z轴自旋;type 6,以0.3 r/min的角速度绕Z轴自旋;type 7,以1.8 r/min的角速度绕Z轴自旋),如图6所示,图中标注了各异常姿态的偏转角。

由于形状不同的空间目标的OCS也不同,因此本文分别探究风云、GPS和天宫三种空间目标的异常检测情况。基于3D Max建立三种空间目标的3ds模型,如图7所示。本体包覆聚酰亚胺薄膜,帆板采用GaSa材料。本文采用改进Phong模型对材料的光学特性进行描述,相应参数见表1^[11],其中:ρ_d和ρ_s分别为表面材料的漫反射率和镜面反射率;α为镜向指数;δ为观测方向与镜面反射方向的夹角(单位为rad);系数a>0,用以调节菲涅耳反射的强度;系数b>0,用以调节镜面反射分量的增降速度。

图 5. 空间目标正常工作状态示意图

Fig. 5. Schematic of space object with normal working status

下载图片 查看所有图片

图 6. 空间目标异常工作状态示意图。(a)偏转30°;(b)偏转90°;(c)偏转130°;(d)绕Z轴自旋;(e)对日定向

Fig. 6. Schematic of space object with abnormal working status. (a) Offset of 30°; (b) offset of 90°; (c) offset of 130°; (d) spinning about body Z; (e) sun alignment with nadir constraint

下载图片 查看所有图片

图 7. 三种空间目标模型。(a)风云;(b) GPS;(c)天宫

Fig. 7. Three space object models. (a) Fengyun; (b) GPS; (c) Tiangong

下载图片 查看所有图片

历史数据为2018和2019年两年内空间目标正常运行的OCS观测数据;异常数据为2020年空间目标分别以以上7种异常状态运行的OCS观测数据(本文以采用系统工具箱STK和OpenGL拾取技术生成的相应的OCS仿真数据为研究对象)。判定依据为,计算历史数据库中2018年与2019年各相似弧段上OCS曲线间的DTW距离,将其作为判定是否异常的阈值(其中2019年观测数据作为当前观测(2020年)的参考数据库,2018年仅用于计算阈值)。具体步骤如下:

1)历史数据库的建立

设置STK场景时间为1 Jan. 2019 00:00 to 31 Dec. 2019 00:00,步长为1 s,状态为正常工作状态;输出空间目标在每一时刻的观测几何条件信息,并通过可观测条件(即空间目标被太阳直射且测站处于全影区或半影区)筛选得到248条可观测弧段;剔除时长小于100 s的弧段后,将剩余的186条弧段的观测几何信息与空间目标3ds模型一并导入到基于OpenGL拾取技术的空间目标光学横截面积计算方法程序^[7]中,得到186条OCS观测曲线,并将其存入历史数据库中。

2)异常判定阈值的建立

设置STK场景时间为1 Jan. 2018 00:00 to 31

Dec. 2018 00:00,步长为1 s,状态为正常工作状态;筛选后得到185条时长大于100 s的可观测弧段,并基于这些弧段生成相应的OCS曲线。

相似弧段的搜索:针对2019年场景中的每一条可观测弧段,按前文介绍的搜索算法,在2018年中找到与其最相似的弧段bestaccess。通过搜索算法共确定185对相似弧段,对应相位角曲线的DTW距离如图8所示,典型相似弧段如图9所示。

图 8. 相似弧段相位角曲线间的DTW距离

Fig. 8. DTW distance between phase angle curves of similar access

下载图片 查看所有图片

图 9. 典型相似弧段的相位角变化情况

Fig. 9. Typical examples of phase angle variation in similar access

下载图片 查看所有图片

基于DTW距离的阈值确定:针对以上搜索到的185对相似弧段,分别计算它们对应的OCS曲线之间的DTW距离,计算得到的风云、GPS、天宫的DTW距离如图10所示,典型相似弧段对应的OCS曲线的变化如图11所示。

由图10可知,在正常的工作状态下,风云、GPS、天宫与历史数据的相似弧段对应的OCS曲线的DTW距离最大分别为0.896,4.193,4.317 m²。因此,异常判定的DTW阈值分别为η_Fengyun=0.896 m²、η_GPS=4.193 m²、η_Tiangong=4.317 m²,即当前弧段与参考弧段的OCS曲线的DTW距离超过阈值时,判定当前工作状态为异常。

3)对7种异常工作状态的检测

分别对1 Jan. 2020 00:00 to 31 Dec. 2020 00:00时间内三种模型的7种异常工作状态进行仿真。基于以上介绍的搜索算法找出历史数据中的相似弧

图 10. 相似弧段对应的OCS曲线间的DTW距离

Fig. 10. DTW distance between OCS curves of similar access

下载图片 查看所有图片

图 11. 典型相似弧段的OCS变化情况

Fig. 11. Typical OCS variation of similar access

下载图片 查看所有图片

图 12. 风云的异常工作状态检测。(a)指向异常;(b)旋转异常

Fig. 12. Abnormal working status detection of Fengyun. (a) Abnormal pointing; (b) abnormal rotation

下载图片 查看所有图片

段,并计算每对相似弧段对应的OCS曲线的DTW距离,如图12~14所示。

定义异常识别率公式为

ξ=nN×100%(3)

式中:ξ为异常识别率;n为识别出工作状态为异常的弧段数目;N为待测弧段总数目。计算得到的异常识别率见表2。

由表1可知:当工作状态仅为姿态指向改变时,异常识别率为84.68%;当工作状态由三轴稳定变为自旋时,能够较准确地识别出异常,识别率为97.96%,并且随着转速变快,识别率增大;空间目标

图 13. GPS的异常工作状态检测。(a)指向异常;(b)旋转异常

Fig. 13. Abnormal working status detection of GPS. (a) Abnormal pointing; (b) abnormal rotation

下载图片 查看所有图片

图 14. 天宫的异常工作状态检测。(a)指向异常; (b)旋转异常

Fig. 14. Abnormal working status detection of Tiangong. (a) Abnormal pointing; (b) abnormal rotation

下载图片 查看所有图片

模型为风云、GPS和天宫时的识别率分别为91.27%、86.80%、93.05%;总体识别率可以达到90.37%。这一结果体现了本文检测方法的有效性。

因此,可得出结论:异常检测的识别率与空间目标具体的工作状态有关,工作状态变化越剧烈,异常识别率就越高;同时,异常检测的识别率还与空间目标的形状结构有关,其表面结构越复杂,形状特征越突出,工作状态的改变就越容易反映到其光学特性的变化上,就越容易识别出异常。

4 结论

针对低轨空间目标工作状态异常检测的问题,本文首先通过DTW算法从历史数据中找出与当前观测几何条件相似的弧段,进而计算两相似弧段上正常工作空间目标OCS曲线间的DTW距离,并将其作为判断空间目标是否异常的阈值。最后利用DTW阈值判断风云、GPS、天宫三个空间目标的7种工作状态,平均准确率可达到90.37%,说明了本文方法的有效性。