AMBL株式会社 Advent Calendar 2022の11日目の記事です．本稿ではAutoEncoderを用いた画像の異常検知を解説します．
不良品の発見を自動化したい
画像の異常検知
工業製品の製造ラインにおいて，疵などの欠陥を検知するために人工知能が使われています．正常な製品の中から不良品を自動で検出する技術のことを異常検知と呼びます．製造ライン上に設置したカメラで撮影した画像に対して異常検知を行うため，画像に対して人工知能を適用することになります．
従来は人手でやっていた不良品の検出を人工知能で代替することにより，人件費の削減や検出の精度向上などが期待できます．異常検知が人工知能に完全に代替されれば，人間の労働が不要になるので素晴らしいですネ！そのような素敵な可能性を秘めている人工知能を利用した異常検知について解説します．
不良品のデータ収集は難しい
画像に対してよく使われる人工知能の手法として深層学習があります．深層学習を行うためには一般的に大量のデータを必要とします．
ここで正常と異常の分類を深層学習で行うことを考えます．学習をうまく行うためには正常なデータの数と異常なデータの数を同じ程度用意して学習する必要があります．ところが，正常なデータと比較すると異常なデータは少ないという問題があります．工場などではそれほど多くの不良品が製造されることは考えづらいからです．また不良品のバリエーションは多岐にわたるため，全ての不良品のバリエーションについてデータを網羅的に集めるのが難しいという問題もあります．
そこで深層学習の手法の一つであるAutoEncoderを応用することが提案されています．AutoEncoderでは「自分が学習した内容のみを出力できる」ということが重要です．次の章で詳しく解説します．
AutoEncoderの異常検知への応用
次元削減と再構成を行う
AutoEncoderは2006年にジェフリー・ヒントンらによって提案された次元削減のアルゴリズムですが，異常検知にも応用されています．AutoEncoderで使われているEncoderとDecoderという構造はノイズ除去や自動着色，超解像などにも応用されているため覚えておくと様々なタスクへの応用が利きます．
下図に示すようにAutoEncoderはEncoderとDecoderという構造に大きく分けて考えることができます．
Encoderの役割
Encoderは入力画像の次元を削減する役割があります．上図の例のように，サイズが28×28［px］の画像を入力したとします．次元の削減とは，28×28=784次元のデータをより少ない次元数，例えば64次元で表現できるように変換することを意味します．ただ闇雲に次元を減らすのではなく元のデータが持つ特徴をできるだけ損なわないように小さな次元に変換する必要があります．
Decoderの役割
DecoderはEncoderで次元を削減した状態から，元の入力画像に近づけるように画像を復元する役割があります．Decoderにより出力された画像のことを再構成画像と呼びます．
全体としてみると，入力画像をEncoderで次元削減し，Decoderにより次元削減したデータから入力画像に近づけるように画像を復元するという流れになっています．学習する際には，入力画像と同じ画像を正解ラベルとして使用します．
覚えたことだけを再現できる
では，なぜこの仕組みが異常検知に適用できるのでしょうか．ポイントはAutoEncoderは学習した画像のみを再構成できるということです．すなわち学習した画像と似たタイプの画像であればAutoEncoderは入力画像とソックリな画像を再構成できます．言い換えれば，学習した画像と似ていない画像は正しく再構成できません．次のような例を考えます．
AutoEncoderを正常な製品の画像（正常画像）で学習したとします．このAutoEncoderに正常画像を入力して推論したら入力画像とソックリな再構成画像が得られます．他方で欠陥のある製品の画像（異常画像）を入力して推論したらどうなるでしょうか？AutoEncoderは学習データに存在しなかった異常画像の再構成画像をうまく作れないはずです．この復元できない部分が異常であると予測できます．
上述のことからAutoEncoderの学習には正常画像と比較して収集が難しい異常画像は必要ないと言えます．異常画像はテストに用いる少量だけを用意すればよいです．これこそがAutoEncoderを使うご利益です．
さて，異常検知をするためには，「うまく再構成できない」ということを定量的に表現する必要があります．これを表す指標を異常度と呼びます．次の節で簡単に説明します．
どれほど異常であるかを数値で表す
異常度としてはL1距離やL2距離を使うことができます．また近年では画像の類似度を測る指標であるSSIMが異常度として有用であるという報告1もあります．ここではL1距離を例にとって計算の手順（ステップ1，2）を解説します．
ステップ1：異常度のヒートマップを算出
まず，入力画像と再構成画像の各画素値の差の絶対値を計算します．これが各画素値に対する異常度です．異常度を各画素ごとに計算すると，入力画像や再構成画像と同じサイズを持ち，異常度が格納されているテンソルが得られます．各画素値の差が大きければ大きいほど再構成がうまくできていない（＝異常度が高い）ことを意味します．したがって，この各画素値の差の絶対値が格納されたテンソルは異常度のヒートマップとして考えることができます．差の絶対値が大きいほど異常度が高いので，グレースケールの場合は白い部分ほど異常度が高いというヒートマップになります．
ステップ2：画像全体としての異常度を算出
次に，入力画像と再構成画像の各画素値の差の絶対値が格納されたテンソルの成分について和を計算します．つまり，画像全体で考えると異常度はどれくらいかということ計算します．本稿では，このように画像毎に異常度を計算することで，正常と異常の判定を行います．
次で解説する実装では損失関数にL2距離，異常度にL1距離を用います．
Pythonで異常検知を実装してみよう！
ここまで説明したことを体感していただくために，AutoEncoderを用いた簡単な異常検知にチャレンジします．手書き文字の認識用データセットとして有名なMNISTを用います．正解ラベルが「1」の画像を正常画像，正解ラベルが「9」の画像を異常画像と定義して実験を行います．Google Colaboratoryで実行することを想定します．
Googleドライブ内のデータの読み書きを可能にする
from google.colab import drive
drive.mount('/content/drive')
必要なライブラリのインポート
import matplotlib.pyplot as plt
import numpy as np
import tensorflow as tf
import os
import csv
from sklearn.model_selection import train_test_split
from tensorflow.keras import layers, losses
from tensorflow.keras.datasets import mnist
from tensorflow.keras.models import Model, load_model
from tensorflow.keras.callbacks import ModelCheckpoint, EarlyStopping, TensorBoard
学習・テスト用のデータを準備
学習・テスト用のデータとしてMNISTを用います．MNISTの中で正解ラベルが「1」のデータを正常，正解ラベルが「9」のデータを異常と定義します．学習には正解ラベルが「1」のデータのみ，テストには正解ラベルが「1」と「9」のデータを用います．以前，手書き数字の「1」と「9」を読み間違えて大変なことになったことがあるので，過去への復讐のためにこの設定にしました．
# MNISTのデータをダウンロード
(x_train, y_train), (x_test, y_test) = mnist.load_data()
# 画素値を0～255から0～1に変換
x_train = x_train.astype('float32') / 255.0
x_test = x_test.astype('float32') / 255.0
# 学習データとして「1（正常）」のデータを抽出
train_input = []
for i in range(len(x_train)):
if y_train[i] == 1:
train_input.append(x_train[i])
x_train = np.array(train_input)
print('学習データ数: ', x_train.shape[0])
# テストデータとして「1（正常）」と「9（異常）」の画像を抽出
test_input = []
test_label = []
for i in range(len(x_test)):
if y_test[i] == 1 or y_test[i] == 9:
test_input.append(x_test[i])
test_label.append(y_test[i])
x_test = np.array(test_input)
y_test = np.array(test_label)
print('テストデータ数: ', x_test.shape[0])
AutoEncoderのアーキテクチャを定義
AutoEncoderのアーキテクチャはTensorFlowのチュートリアル2をもとに構成しました．このAutoEncoderには28×28=784次元のデータを64次元までエンコードし，元のサイズ（784次元）にデコードする働きがあります．全結合層のみで構成されています．
latent_dim = 64
class Autoencoder(Model):
def __init__(self, latent_dim):
super(Autoencoder, self).__init__()
self.latent_dim = latent_dim
self.encoder = tf.keras.Sequential([
layers.Flatten(),
layers.Dense(392, activation='relu'),
layers.Dense(latent_dim, activation='relu'),
])
self.decoder = tf.keras.Sequential([
layers.Dense(392, activation='relu'),
layers.Dense(784, activation='sigmoid'),
layers.Reshape((28, 28))
])
def call(self, x):
encoded = self.encoder(x)
decoded = self.decoder(encoded)
return decoded
autoencoder = Autoencoder(latent_dim)
モデルの保存先やコールバックの設定などを行う
# 学習済みモデルとTensorBoard用のログの保存先を作成
save_model_dir = '/content/drive/MyDrive/AutoEncoder/AdvenCalendar/Model/'
log_dir = '/content/drive/MyDrive/AutoEncoder/AdvenCalendar/LOG/'
os.makedirs(save_model_dir, exist_ok=True)
os.makedirs(log_dir, exist_ok=True)
# コールバックを定義
my_callbacks = [
ModelCheckpoint(filepath=save_model_dir, monitor='val_loss', save_best_only=False),
EarlyStopping(monitor='val_loss', min_delta=0, patience=10),
TensorBoard(log_dir=log_dir)
]
モデルをコンパイル
オプティマイザはAdam，損失関数はL2距離をミニバッチ数で割った値（平均二乗誤差，Mean Squared Error）を用います．
autoencoder.compile(optimizer='adam', loss=losses.MeanSquaredError())
学習を実行
autoencoder.fit(x_train, x_train, epochs=300, shuffle=True,
validation_data=(x_test, x_test), callbacks=my_callbacks)
テストデータに対する推論を実行
# モデルの読み込み
saved_model = load_model(save_model_dir)
# テストデータに対する推論
encoded_imgs = saved_model.encoder(x_test).numpy()
decoded_imgs = saved_model.decoder(encoded_imgs).numpy()
推論結果の分析1
下図のように，入力画像（original，上段），再構成画像（reconstructed，中段），入力画像と再構成画像の差分画像（下段）を表示します．
差分画像は白い部分ほど異常度が高いことを表すヒートマップとして考えることができます．また，差分画像の上に表示される数値は画像毎の異常度の総和（L1距離）です．以降では，この数値を用いて画像毎に異常と正常を判定します．
結果を見ると，「9（異常）」は「1（正常）」と比較して再構成がうまくいっておらず，異常度も高くなっています．期待通りの結果が得られました．この分析を実行するコードは次の通りです．
n = 5
plt.figure(figsize=(10, 7))
for i in range(n):
# 入力画像を表示（上段）
ax = plt.subplot(3, n, i + 1)
plt.title("original")
plt.imshow(tf.squeeze(x_test[i]))
plt.gray()
ax.get_xaxis().set_visible(False)
ax.get_yaxis().set_visible(False)
# 再構成画像を表示（中段）
bx = plt.subplot(3, n, i + n + 1)
plt.title("reconstructed")
plt.imshow(tf.squeeze(decoded_imgs[i]))
plt.gray()
bx.get_xaxis().set_visible(False)
bx.get_yaxis().set_visible(False)
# 入力画像と再構成画像の差分から異常部位を示すヒートマップを作成（下段）
# ヒートマップでは色が白い部分ほど異常度が高い
# 画像の上部の数字は画像の異常度の総和を計算した値
cx = plt.subplot(3, n, i + 2 * n + 1)
# 入力画像と再構成画像の画素値の差の絶対値を計算 → 異常度のヒートマップに対応
diff_img = np.abs(tf.squeeze(x_test[i]) - tf.squeeze(decoded_imgs[i]))
# 画像の異常度の総和を計算．
anomalous = np.sum(diff_img)
plt.title(anomalous)
plt.imshow(tf.squeeze(diff_img))
plt.gray()
cx.get_xaxis().set_visible(False)
cx.get_yaxis().set_visible(False)
plt.show()
推論結果の分析2
画像毎の異常度のヒストグラムを作成することで正常画像と異常画像における異常度の分布を確認します．つまりデータ全体で見たときに正常画像と異常画像が分離しているかどうかを確かめるということです．その結果が下図です．横軸が異常度，縦軸が頻度（画像枚数）を表しています．
ヒストグラムより，正常画像「1」（赤色で表現）と異常画像「9」（青色で表現）が画像毎の異常度という尺度を用いることである程度分離できることが分かりました．目分量では異常度のしきい値を20くらいに設定すればうまく分類できそうです．次の分析3では，正常と異常のしきい値をどこに決めたらよいか，計算により決定します．ここまでの分析を実行するコードは次の通りです．
# 画像ごとの異常度をヒストグラムとして出力
output_path = '/content/drive/MyDrive/AutoEncoder/AdvenCalendar/anomaly_score/'
txt_name_normal = 'normal.txt'
txt_name_anomaly = 'anomaly.txt'
input_normal_path = output_path + txt_name_normal
input_anomaly_path = output_path + txt_name_anomaly
output_histgram = '/content/drive/MyDrive/AutoEncoder/AdvenCalendar/anomaly_score/'
normal_score =[]
with open(input_normal_path, 'r') as f:
reader = csv.reader(f)
for row in reader:
row = int(float(row[0]))
normal_score.append(row)
anomaly_score =[]
with open(input_anomaly_path, 'r') as f:
reader = csv.reader(f)
for row in reader:
row = int(float(row[0]))
anomaly_score.append(row)
plt.title("Anomaly Score Histgram")
plt.xlabel("Anomaly Score")
plt.ylabel("Frequency")
plt.hist(normal_score, bins=10, alpha=0.3, histtype='stepfilled', color='r', label="normal: 1")
plt.hist(anomaly_score, bins=10, alpha=0.3, histtype='stepfilled', color='b', label='anomaly: 9')
plt.legend(loc=1)
plt.savefig(output_histgram + 'histgram.png')
plt.show()
plt.close()
推論結果の分析3
正常と異常を最もよく分離できる異常度のしきい値がいくつであるかを計算により決定します．ここではF1値を最大にするようなしきい値が最適であるとします．次のような結果が得られました．
しきい値: 19
Precision: 0.9920212765957447
Recall: 0.9859030837004406
F1値: 0.9889527176314628
「1（正常）」と「9（異常）」をうまく分類できるようなしきい値を決定することができました．この計算は次のコードで実行できます．
f1_best = 0
threshold_best = 0
precision_best = 0
recall_best = 0
for threshold in range(0, 40):
# Positiveを正常，Negativeを異常と設定
tp_count = 0
fn_count = 0
for score in normal_score:
# 正常であるのに異常と間違えて判定
if score >= threshold:
fn_count += 1
# 正常であると正しく判定
else:
tp_count += 1
tn_count = 0
fp_count = 0
for score in anomaly_score:
# 異常であると正しく判定
if score >= threshold:
tn_count += 1
# 正常であると間違えて判定
else:
fp_count += 1
try:
# precision
p = tp_count / (tp_count + fp_count)
# Recall
r = tp_count / (tp_count + fn_count)
except ZeroDivisionError:
continue
# F1
f1 = 2.0 * r * p / (r + p)
if f1 > f1_best:
threshold_best = threshold
f1_best = f1
precision_best = p
recall_best = r
print('しきい値: ', threshold_best)
print('Precision: ', precision_best)
print('Recall: ', recall_best)
print('F1値: ', f1_best)
おわりに
AutoEncoderを利用した異常検知について解説しました．MNISTではうまく異常検知ができましたが，実データに対してうまくいくかどうかは覚束ないところがあります．MNISTはあまりにも簡単なデータであるからです．今後は，実データに近い異常検知用データセットであるMVTecのデータセット3を用いた場合はどれくらいうまくいくのか実験してみたいです．
-
Improving Unsupervised Defect Segmentation by Applying Structural Similarity to Autoencoders，https://arxiv.org/abs/1807.02011 ↩
-
Intro to Autoencoders，https://www.tensorflow.org/tutorials/generative/autoencoder ↩
-
THE MVTEC ANOMALY DETECTION DATASET (MVTEC AD)，https://www.mvtec.com/company/research/datasets/mvtec-ad ↩