人工衛星テレメトリのデータセットtelemanomを利用する

NASA JPLのHundman et al 氏が 2018年に発表した自身の論文「Detecting Spacecraft Anomalies Using LSTMs and Nonparametric Dynamic Thresholding」で、Github上で公開しているデータセットtelemanomの使い方について紹介します。

この論文自体は、2022年現在500件以上引用されており、センサーや機器データの異常検知の論文の検証用データセットとして広く利用されているようです。

本記事の内容

Telemanomのデータセットの入手方法
Telemanomデータセットの解説
データセット中の異常個所を簡単に可視化する方法

データセットの概要

NASAが打ち上げている2つの人工衛星（宇宙機）SMAPとMSLのデータセットです。

telemanomには、これらの2つの宇宙機の正常なテレメトリと異常のテレメトリが含まれており、異常のテレメトリに関しても異常区間が教師ラベル付きで与えられており、非常に扱いやすいデータセットとなっています。

ちなみにテレメトリーとは宇宙機の内部状態や内部ソフトウェアのパラメータ値などのデータのことをいいます。

SMAP(Soil Moisture Active/Passive Mission)は、2015年に打ち上げられた地球観測衛星で、Lバンドの合成開口レーダーによって地球の土壌水分を計測する衛星です。一方、MSL(Mars Science Laboratory)は、2011年に打ち上げられた火星探査機です。

2つとも宇宙機ですが、SMAPは地球周回しているため比較的地球近傍にいる人工衛星であり、MSLは火星に向かう探査機であり打ち上げ以降は地球を離れるといった違うはありますが、それ以外の人工衛星のバス部分の基本的な機器の仕組みは、データ解析をする上でそこまで意識する必要はないと思われます。

データセットの入手方法

まず、正解ラベルが入ったcsvシートと一連のデータをダウンロードします。

まず正解ラベルのcsvは下記で入手できます。

curl -O https://raw.githubusercontent.com/khundman/telemanom/26831a05d47857e194a7725fd982d5dea5402dd4/labeled_anomalies.csv

続いて、一連のデータセットです。Github上にあるように、クラウドにデータが挙げられているのでそちらのURLからダウンロードします。

curl -O https://s3-us-west-2.amazonaws.com/telemanom/data.zip && unzip data.zip && rm data.zip

データセットの解説

データセットの中身を見ていきます。まず、必要なライブラリをインポートします。

import pandas as pd
import numpy as np

import matplotlib.pyplot as plt

まず、最初にダウンロードしたlabeled_anomalies.csvの中身を見てみます。

このファイルは時系列データのどの時点(index)が異常であったかを示しており、このファイルの内容をもとに、時系列の異常箇所を見ていくことになります。

df = pd.read_csv('labeled_anomalies.csv')
df

各系列とその概要はこのようになっています。

chan_id	データセットの系列	それぞれ秘匿化されているが、先頭のアルファベットが衛星のコンポーネントを示しているアルファベットとコンポーネントの対応は下の表を参照
spacecraft	宇宙機の種類	SMAP, MSL の2種類
anomaly_sequences	異常値のインデックス	[[5300, 5747]] のように、範囲で異常値のindexを指定している
class		contextual, point
num_values	testデータセットの数	test 用の時系列データの点数

chan_id は、A-1やP-1のように各アルファベット1文字と数字から構成されています。

最初のアルファベットは宇宙機の構成要素（電源系なり温度なり）を示しているようです。

Github上のIssue上や論文での内容をもとに、簡単な対応表をまとめました。一部不明なところがあるので、分かり次第追記していこうと思います。

A	不明
P	Power (電源系)
T	Temperature(温度)
R	Radiation (放射線)
M	Memory
F	不明
D	不明
I	Instrumentation (計装)

telemanomに含まれる異常データ

正解ラベルファイルの内容がわかったので、実際のデータセットを見て、異常データの確認をしています。

ダウンロードしたzipファイルを展開すると、下記にようになっています。trainとtestにデータセットが分かれており、channelごとにデータを扱えるようになっています。非常に整理されており、使いやすいデータセットです。

data
├ 2018-05-19_15.00.10
├ test
│ ├ A-1.npy
│ ├ A-2.npy
│ ├ A-3.npy
│ └ ...
└ train
　 ├ A-1.npy
　 ├ A-2.npy
　 ├ A-3.npy
　 └ ...

SMAP電源系の異常データE-2を見る

telemanomデータセットには、82種類の異常のデータセットが含まれています。

試しにSMAPのE-2チャネルの異常を見ていきます。

まず、E-2のtrainデータとtestデータを読み込み、データ点数を確認します。

train_df = np.load('./data/train/E-2.npy')
test_df = np.load('./data/test/E-2.npy')

print(train_df.shape)
print(test_df.shape)

(2880, 25)
(8532, 25)

訓練データが、2880点、テストデータが8532点あることがわかりました。

これらを可視化して確認していきます。

まず、正常な訓練データです。train, test共に25系列ありますが、最初の0列目が実際のテメトリーの値のようです。そのほかの系列は、おそらくコマンドがいつ送られたかを示しています。(エンコードされており、どのようなコマンドがSMAPに送信されたかは分からないようになっています。

telemetries = train_df[:, 0]

x = np.arange(0, telemetries.shape[0], 1)
fig, ax = plt.subplots(figsize=(20, 4), dpi=120)
ax.plot(x, telemetries)

心電図のような、周期的な波形から構成される時系列が確認できます。

続いて異常が含まれる、テストデータを可視化します。

telemetries = test_df[:, 0]
x = np.arange(0, telemetries.shape[0], 1)
fig, ax = plt.subplots(figsize=(20, 4), dpi=120)

ax.plot(x,telemetries)

5000~6000付近で、波形の形状が1つだけおかしい点がありますね。ここが異常箇所のように思われます。

labeled_anomalies.csvには、testデータに含まれる異常値の正解ラベルとなるindexが含まれるので、実際に異常の場所を確認してみましょう。

df[df['chan_id']=='E-2']["anomaly_sequences"]

[[5598, 6995]]

5598~6995のあたりが異常値になっているようで、予想が当たりました。

この辺りをもっと拡大してみていきます。

異常値を示すラベルが文字列形式になっており少し扱いにくいのでtokenizerとparser関数を（超適当に）実装し、Pythonから扱いやすいようにします。

def tokenizer(labels):
    tokens = []
    i = 0
    while (i < len(labels)):
        val = labels[i]
        if val == " ":
            i += 1
            continue
            
        if (val == "[" or val == "]" or val == ","):
            tokens.append(val)
            i += 1
            continue
        
        if val.isdigit():
            j = 1
            while (labels[i+j].isdigit()):
                j += 1
            
            digit = labels[i:i+j]
            tokens.append(digit)
            
            i += j
            continue
            
        print("un known format: {}".format(val))
        break
            
    return tokens

def parser(token):
    i = 0
    stack = []
    result = []
    while (i < len(token)):
        val = token[i]
        if val.isdigit():
            stack.append(val)
            
        elif (val == "]"):
            try:
                end = int(stack.pop())
                start = int(stack.pop())
                result.append([start, end])
            except IndexError:
                return result
                
        i += 1
    return result

ここまでで準備は整いました。

あとは、これらを利用することで異常コードを読み込んでグラフ化することができます。

target = "E-2" # 任意に指定
anomaly_sequences = df[df['chan_id']==target]["anomaly_sequences"].values[0]
token = tokenizer(anomaly_sequences)
labels = parser(token)

test_df = np.load("./data/test/{}.npy".format(target))
telemetries = test_df[:, 0]
x = np.arange(0, telemetries.shape[0], 1)

fig, ax = plt.subplots(figsize=(20, 4), dpi=120)

ax.plot(x,telemetries)

for _, label in enumerate(labels):
    ax.axvspan(label[0], label[1], color="red", alpha=0.2)

異常個所をハイライトして、時系列を可視化することができました。

他の異常個所も諸々見てみる

異常データが82種類もあるので、アルゴリズムや機械学習モデルの検証にかなり使えそうだなと思いました。

他の異常の例もいくつか掲載します。

F-2の異常パターン

P-14の異常パターン

さまざまな異常パターンがありますね。SMAPでは合計55個の異常データ、MSLでは27ヶ所の異常データがあるので、一気にまとめてみていきましょう。

tests = df[df["spacecraft"] == "SMAP"]
num_col = 2

width, height = 20, 3 * (len(tests) // num_col + 1)
fig, ax = plt.subplots((len(tests) + 1) //num_col, num_col, figsize=(width, height), facecolor='w')

for idx, row  in tests.iterrows():
    target = row["chan_id"]
    c_row = idx // num_col
    c_col = idx % num_col
    
    anomaly_sequences = df[df['chan_id']==target]["anomaly_sequences"].values[0]
    token = tokenizer(anomaly_sequences)
    labels = parser(token)

    test_df = np.load("./data/test/{}.npy".format(target))
    telemetries = test_df[:, 0]
    x = np.arange(0, telemetries.shape[0], 1)
    ax[c_row, c_col].set_title(target)
    ax[c_row, c_col].plot(x,telemetries)
        
    for _, label in enumerate(labels):
        ax[c_row, c_col].axvspan(label[0], label[1], color="red", alpha=0.2)

plt.show()

これで一気に全部のテレメトリを可視化することができました。

画像サイズの都合で画質を荒くしているので、自分の手元でコードを動かして、確認してみることをお勧めします。

telemanomには、テレメトリのデータセットだけでなく、エンコードされたコマンドの入力も含まれており、いつどのような指令が衛星や探査機に送出されたかも分かっています。

エンコードされているので、あまり有用に使えないかもしれませんが、気になる人は是非確認してみてください。

telemanomのライセンスについて

ライセンスについては、配布元に記載されているので、利用する際にはしっかり確認するようにしましょう。