[Paper Review] Detecting Spacecraft Anomalies Using LSTMs and Nonparametric Dynamic Thresholding

본 논문은 2018년 NASA에서 발표한 논문으로 이후 많은 시계열 이상치 탐지 논문들에서 인용한 유명한 논문입니다. 저도 업무 중 이상치 탐지 모델을 계발하는 데 이 논문의 방법을 사용했기 때문에 이 기회에 논문의 리뷰를 해보도록 하겠습니다.

1. Introduction

우주선은 온도, 방사선, 전력 등과 같은 feature들을 측정하는 수천 개의 채널을 가진 매우 복잡하고 비싼 기계입니다. 이러한 채널을 모니터링하는 것은 복잡성과 비용을 고려했을 때 중요하고 필수적인 요소입니다. 즉, 잠재적 위험을 미리 감지하는 것은 중요합니다.
기존의 이상치 탐지 방법은 rule 기반으로 사람이 정한 특정 기준을 벗어나면 anomaly라고 판단하는데, 이는 전문인력의 도메인이 필요하며 특정 기준이 계속해서 update되어야 하는데 이는 너무 많은 비용이 듭니다.

1.1 Contributions

본 논문에서는 spacefraft data 뿐만 아니라 일반적인 multivariate time series data에서도 적용할 수 있는 방법론을 제안합니다.
특히 높은 예측성능과 해석력을 모두 얻기 위해 LSTM을 사용한 예측모델을 제안하고, nonparametric, dynamic한 비지도 thresholding 방법론을 제안합니다.

2. Methods

다음 방법은 LSTM을 사용하는 비지도 이상 탐지 접근법의 핵심 구성 요소를 설명합니다. 새로운 비지도 임계값 방법을 사용하여 측정 데이터의 수천 채널을 자동으로 평가하고 예측 오류가 이상을 나타내는지 결정합니다.

2.1 Telemetry Value Prediction with LSTMs

2.1.1 Single-Channel Models

각 채널에 대한 단일 모델이 생성되고 각 모델은 해당 채널의 값을 예측하는 데 사용됩니다. 즉, multivariate model을 사용하는 것이 아니라 univariate model을 사용합니다. 논문에서는 m-dim output을 한꺼번에 예측하는 모델을 만들면, m이 클 때 LSTM의 성능이 떨어지고 이는 anomaly를 놓칠 수 있기 때문이라고 설명합니다. 하지만 제 생각에 이는 각 채널들의 영향력을 무시할 수 있다고 생각합니다.

2.1.2 Predicting Values for a Channel

각각의 time step마다 바로 다음 시점만 예측하고 error를 비교합니다. 즉, 한 시점만 예측하는 모델인 셈인데 이는 processing time이 적은 모델이 필요하기 때문이라고 합니다.

2.2 Errors and Smoothing

우선 예측값 $\hat{y}^{(t)}$에 대해, error는 $e^{(t)}=|y^{(t)}-\hat{y}^{(t)}|$로 정의됩니다. 즉, error vector $\mathbf{e}$는 다음과 같이 정의됩니다.
$\mathbf{e}=[e^{(t-h)},…,e^{(t-l_s)},…,e^{(t-1)}, e^{(t)}]$, $h$ : 과거 error값의 개수
error vector $\mathbf{e}$는 급증되는 error를 억제하기 위해 EWMA로 smoothing처리를 해줍니다.
$\mathbf{e}_s=[e_s^{(t-h)},…,e_s^{(t-l_s)},…,e_s^{(t-1)}, e_s^{(t)}]$

2.3 Threshold Calculation and Anomaly Scoring

본 논문에서는 레이블이나 통계적 가정이 필요없고 낮은 오버헤드로 높은 성능을 달성하는 비지도 방법을 제안합니다.
threshold $\epsilon$은 다음과 같이 $\boldsymbol{\epsilon}$ 집합에서 선택되어 집니다.
$\boldsymbol{\epsilon}=\mu(\boldsymbol{e}_s) + \boldsymbol{z}\sigma(\boldsymbol{e}_s)$

$\begin{align} \epsilon=\argmax(\boldsymbol{\epsilon})={ {\vartriangle \mu(\boldsymbol{e}s)/\mu(\boldsymbol{e}_s)+\vartriangle \sigma(\boldsymbol{e}_s)/\sigma(\boldsymbol{e}_s)} \over {|\boldsymbol{e}_a|+|\boldsymbol{E}{seq}|^2}} \end{align}$ such that:
$\begin{align} \& \vartriangle \mu(\boldsymbol{e}s)=\mu(\boldsymbol{e}_s)-\mu({e_s \in \boldsymbol{e}_s | e_s < \epsilon}) \& \vartriangle \sigma(\boldsymbol{e}_s)=\sigma(\boldsymbol{e}_s)-\sigma({e_s \in \boldsymbol{e}_s | e_s < \epsilon}) \& \boldsymbol{e}_a={e_s \in \boldsymbol{e}_s|e_s> \epsilon} \& E{seq}=\textnormal{continuous sequences of} \; e_a \in \boldsymbol{e}_a \end{align}$

위 수식을 코드로 이해해 봅시다. 다음의 데이터를 smoothing한 error라고 가정해보면,
먼저 $\mu(\boldsymbol{e}s)$와 $\sigma(\boldsymbol{e}_s)$를 정의하고, $\vartriangle \mu(\boldsymbol{e}_s)$와 $\vartriangle \sigma(\boldsymbol{e}_s)$, $\boldsymbol{e}_a$, $E{seq}$를 정의해줍니다.
그 전에 적절한 z값을 설정해줘야 하는데 원래는 for문을 돌며 찾습니다. 공식 코드에서는 2.5~10 사이를 0.5간격으로 돕니다. 일단 2.5로 두고 시작해봅니다.

z = 2.5

m = np.mean(s_error)
std = np.std(s_error)

th = m + z * std

delta_m = m - np.mean(s_error[s_error < th])
delta_std = std - np.std(s_error[s_error < th])

ea = np.sum(s_error > th)
idx = np.where(s_error > th)[0]

$E_{seq}$는 mit의 consecutimve_groups를 이용해 구하는데 이 함수는 연속된 인덱스를 묶어주는 함수입니다. 실행시키면 다음과 같은 리스트가 생성됩니다.

import more_itertools as mit

eseq = [list(group) for group in mit.consecutive_groups(idx)]
for i in range(len(eseq)):
    print(eseq[i])
    print('')

[493, 494]
[496, 497, 498, 499, 500, 501, 502]
[569, 570, 571, 572, 573, 574, 575, 576, 577, 578, 579, 580, 581, 582, 583, 584, 585, 586, 587, 588, 589, 590, 591, 592, 593, 594, 595, 596, 597, 598, 599, 600]
[604, 605]
[702, 703, 704, 705, 706, 707, 708, 709, 710, 711, 712, 713, 714, 715, 716, 717, 718, 719, 720, 721, 722]
[932, 933, 934, 935, 936, 937, 938, 939, 940, 941]
[997, 998, 999, 1000, 1001, 1002, 1003, 1004, 1005, 1006, 1007, 1008, 1009, 1010, 1011, 1012, 1013, 1014, 1015, 1016]
[2134, 2135, 2136, 2137, 2138, 2139, 2140, 2141, 2142, 2143, 2144, 2145, 2146, 2147, 2148, 2149, 2150, 2151, 2152, 2153, 2154, 2155, 2156]
[2381]
...

그러면 시각화를 통해 th를 넘늠 continuous sequences인 $E_{seq}$를 시각화해봅니다.

위와 같은 과정들로 각각의 값들을 정의하게 됩니다. 그 후, 분모, 분자를 정의하게 되고 numerator/denominator가 가장 크게 되는 값을 threshold로 결정하게 됩니다.

numerator = ((m - delta_m) / m) + ((std - delta_std) / std)
denominator = ea + eseq ** 2

그러면 numerator/denominator의 최대값은 무엇을 의미할까요?
분모를 작게 하는 것은 이상치의 값과 th를 넘는 연속된 이상치들의 sequence 갯수를 최소화하는 것이고, 분자를 최대화하는 것은 이상치를 제거했을 때 평균과 분산을 크게 하는 것이라고 볼 수 있습니다. 하지만 본 논문에서는 이 방벙의 근거를 충분히 제시하지 않아 아쉬운 점으로 남아있습니다.