深層学習を使用した時系列異常検出

この例では、シーケンスデータまたは時系列データで異常を検出する方法を示します。

シーケンスデータまたは時系列データの集合で異常または異常な領域を検出するには、自己符号化器を使用できます。自己符号化器は、入力を低次元空間に変換 (符号化ステップ) し、低次元表現から入力を再構成 (復号化ステップ) することによって、入力を複製するように学習させるモデルの一種です。自己符号化器の学習に、ラベル付けされたデータは必要ありません。

自己符号化器自体は異常を検出しません。代表的なデータのみを使用して自己符号化器に学習させると、代表的なデータから学習した特徴のみを使用して入力データを再構成できるモデルが得られます。自己符号化器を使用して観測値が異常かどうかをチェックするには、観測値をそのネットワークに入力し、元の観測値と再構成された観測値の間の誤差を測定します。元の観測値と再構成された観測値の間の誤差が大きいことは、自己符号化器の学習に使用されたデータを代表するものではない特徴が元の観測値に含まれており、異常があることを示しています。元のシーケンスと再構成されたシーケンスの間の要素単位で誤差を観測することで、異常の局所領域を特定できます。

以下のイメージは、異常な領域を強調表示したシーケンスの例を示しています。

この例では、波形データセットを使用します。このデータセットには、3 つのチャネルの異なる波長の合成生成波形が 2,000 個含まれています。

学習データの読み込み

WaveformData.mat から波形データセットを読み込みます。観測値は numChannels 行 numTimeSteps 列の配列です。ここで、numChannels と numTimeSteps はそれぞれシーケンスのチャネル数とタイムステップ数です。

load WaveformData

最初のいくつかのシーケンスのサイズを表示します。

data(1:5)

ans=5×1 cell array
    {3×103 double}
    {3×136 double}
    {3×140 double}
    {3×124 double}
    {3×127 double}

チャネル数を表示します。ネットワークに学習させるには、各シーケンスのチャネル数が同じでなければなりません。

numChannels = size(data{1},1)

numChannels = 3

最初のいくつかのシーケンスをプロットに可視化します。

figure
tiledlayout(2,2)
for i = 1:4
    nexttile
    stackedplot(data{i}',DisplayLabels="Channel " + (1:numChannels));
    title("Observation " + i)
    xlabel("Time Step")
end

データを学習区画と検証区画に分割します。データの 90% を使用してネットワークに学習させ、10% を検証用に残しておきます。

numObservations = numel(data);
XTrain = data(1:floor(0.9*numObservations));
XValidation = data(floor(0.9*numObservations)+1:end);

学習用データの準備

この例で作成されるネットワークは、データの時間の次元を係数 2 で繰り返しダウンサンプリングし、次に出力を係数 2 で同じ回数だけアップサンプリングします。ネットワークがシーケンスを明確に再構成して必ず入力と同じ長さにできるようにするには、 $2^{K}$ の最も近い倍数の長さになるようにシーケンスを切り捨てます。ここで $K$ はダウンサンプリング演算の回数です。

入力データを 2 回ダウンサンプリングします。

numDownsamples = 2;

シーケンスを 2^numDownsamples の最も近い倍数に切り捨てます。ネットワーク入力層の最小シーケンス長を計算できるように、シーケンス長を含むベクトルも作成します。

sequenceLengths = zeros(1,numel(XTrain));

for n = 1:numel(XTrain)
    X = XTrain{n};
    cropping = mod(size(X,2), 2^numDownsamples);
    X(:,end-cropping+1:end) = [];
    XTrain{n} = X;
    sequenceLengths(n) = size(X,2);
end

同じ手順で、検証データを切り捨てます。

for n = 1:numel(XValidation)
    X = XValidation{n};
    cropping = mod(size(X,2),2^numDownsamples);
    X(:,end-cropping+1:end) = [];
    XValidation{n} = X;
end

ネットワークアーキテクチャの定義

次のネットワークアーキテクチャを定義します。これは、データのダウンサンプリングとアップサンプリングによって入力を再構成します。

シーケンス入力に対して、入力チャネル数と一致する入力サイズでシーケンス入力層を指定します。z スコア正規化を使用してデータを正規化します。ネットワークが学習データを確実にサポートするように、MinLength オプションを学習データ内で最も短いシーケンスの長さに設定します。
入力をダウンサンプリングするために、1 次元畳み込み層、ReLU 層、およびドロップアウト層から成る反復するブロックを指定します。符号化された入力をアップサンプリングするために、1 次元転置畳み込み層、ReLU 層、およびドロップアウト層から成るブロックを同数含めます。
畳み込み層に対して、フィルターの減少数をサイズ 7 に指定します。出力が係数 2 で必ず均等にダウンサンプリングされるようにするには、ストライドを 2 に指定し、Padding オプションを "same" に設定します。
転置畳み込み層に対して、フィルターの増加数をサイズ 7 に指定します。出力が係数 2 で必ず均等にアップサンプリングされるようにするには、ストライドを 2 に指定し、Cropping オプションを "same" に設定します。
ドロップアウト層に対して、ドロップアウトの確率を 0.2 に指定します。
入力と同じ数のチャネルをもつシーケンスを出力するには、入力のチャネル数と一致する数のフィルターをもつ 1 次元転置畳み込み層を指定します。出力シーケンスが層入力と必ず同じ長さになるようにするには、Cropping オプションを "same" に設定します。
最後に、回帰層を含めます。

ダウンサンプリング層やアップサンプリング層の数を増減するには、学習用データの準備セクションで定義されている変数 numDownsamples の値を調整します。

minLength = min(sequenceLengths);
filterSize = 7;
numFilters = 16;
dropoutProb = 0.2;

layers = sequenceInputLayer(numChannels,Normalization="zscore",MinLength=minLength);

for i = 1:numDownsamples
    layers = [
        layers
        convolution1dLayer(filterSize,(numDownsamples+1-i)*numFilters,Padding="same",Stride=2)
        reluLayer
        dropoutLayer(dropoutProb)];
end

for i = 1:numDownsamples
    layers = [
        layers
        transposedConv1dLayer(filterSize,i*numFilters,Cropping="same",Stride=2)
        reluLayer
        dropoutLayer(dropoutProb)];
end

layers = [
    layers
    transposedConv1dLayer(filterSize,numChannels,Cropping="same")
    regressionLayer];

ネットワークを対話的に表示または編集するには、ディープネットワークデザイナーを使用できます。

deepNetworkDesigner(layers)

学習オプションの指定

学習オプションを指定します。

Adam ソルバーを使用して学習させます。
学習を 120 エポック行います。
すべてのエポックでデータをシャッフルします。
検証データを使用してネットワークを検証します。シーケンスを入力とターゲットの両方として指定します。
学習の進行状況をプロットに表示します。
詳細出力を非表示にします。

options = trainingOptions("adam", ...
    MaxEpochs=120, ...
    Shuffle="every-epoch", ...
    ValidationData={XValidation,XValidation}, ...
    Verbose=0, ...
    Plots="training-progress");

ネットワークの学習

関数 trainNetwork を使用してネットワークに学習させます。自己符号化器に学習させる場合、入力とターゲットは同じです。学習データを入力とターゲットの両方として指定します。

net = trainNetwork(XTrain,XTrain,layers,options);

Figure Training Progress (17-Jan-2023 16:43:54) contains 2 axes objects and another object of type uigridlayout. Axes object 1 contains 8 objects of type patch, text, line. Axes object 2 contains 8 objects of type patch, text, line.

ネットワークのテスト

検証データを使用してネットワークをテストします。検証シーケンスごとに、シーケンスと再構成シーケンスの間の平均絶対誤差 (MAE) を計算します。

YValidation = predict(net,XValidation);

MAEValidation = zeros(numel(XValidation),1);
for n = 1:numel(XValidation)
    X = XValidation{n};
    Y = YValidation{n};
    MAEValidation(n) = mean(abs(Y - X),"all");
end

ヒストグラムで MAE 値を可視化します。

figure
histogram(MAEValidation)
xlabel("Mean Absolute Error (MAE)")
ylabel("Frequency")
title("Representative Samples")

最大 MAE を異常検出のベースラインとして使用できます。検証データから最大 MAE を決定します。

MAEbaseline = max(MAEValidation)

MAEbaseline = 0.5003

異常シーケンスの特定

異常な領域を含むように検証シーケンスの一部を手動で編集して、一連の新しいデータを作成します。

検証データのコピーを作成します。

XNew = XValidation;

変更するシーケンスを 20 個ランダムに選択します。

numAnomalousSequences = 20;
idx = randperm(numel(XValidation),numAnomalousSequences);

選択したシーケンスごとに、データのパッチ XPatch を 4*abs(Xpatch) に設定します。

for i = 1:numAnomalousSequences
    X = XNew{idx(i)};

    idxPatch = 50:60;
    XPatch = X(:,idxPatch);
    X(:,idxPatch) = 4*abs(XPatch);

    XNew{idx(i)} = X;
end

新しいデータで予測を行います。

YNew = predict(net,XNew);

予測ごとに、入力シーケンスと再構成シーケンスの間の MAE を計算します。

MAENew = zeros(numel(XNew),1);
for n = 1:numel(XNew)
    X = XNew{n};
    Y = YNew{n};
    MAENew(n) = mean(abs(Y - X),"all");
end

プロットで MAE 値を可視化します。

figure
histogram(MAENew)
xlabel("Mean Absolute Error (MAE)")
ylabel("Frequency")
title("New Samples")
hold on
xline(MAEbaseline,"r--")
legend(["Data" "Baseline MAE"])

最大 MAE 値をもつ上位 10 個のシーケンスを特定します。

[~,idxTop] = sort(MAENew,"descend");
idxTop(1:10)

最大 MAE 値をもつシーケンスとその再構成をプロットで可視化します。

X = XNew{idxTop(1)};
Y = YNew{idxTop(1)};

figure
t = tiledlayout(numChannels,1);
title(t,"Sequence " + idxTop(1))

for i = 1:numChannels
    nexttile

    plot(X(i,:))
    box off
    ylabel("Channel " + i)

    hold on
    plot(Y(i,:),"--")
end

nexttile(1)
legend(["Original" "Reconstructed"])

異常な領域の特定

シーケンス内で異常な領域を検出するには、入力シーケンスと再構成シーケンスの間の MAE を求め、しきい値を超える誤差がある領域を強調表示します。

入力シーケンスと再構成シーケンスの間の誤差を計算します。

MAE = mean(abs(Y - X),1);

タイムステップのウィンドウサイズを 7 に設定します。検証データを使用して特定された最大誤差値を MAE の値が少なくとも 10% 上回るタイムステップが含まれるウィンドウを特定します。

windowSize = 7;
thr = 1.1*MAEbaseline;

idxAnomaly = false(1,size(X,2));
for t = 1:(size(X,2) - windowSize + 1)
    idxWindow = t:(t + windowSize - 1);

    if all(MAE(idxWindow) > thr)
        idxAnomaly(idxWindow) = true;
    end
end

シーケンスをプロットに表示し、異常な領域を強調表示します。

figure
t = tiledlayout(numChannels,1);
title(t,"Anomaly Detection ")

for i = 1:numChannels
    nexttile
    plot(X(i,:));
    ylabel("Channel " + i)
    box off
    hold on

    XAnomalous = nan(1,size(X,2));
    XAnomalous(idxAnomaly) = X(i,idxAnomaly);
    plot(XAnomalous,"r",LineWidth=3)
    hold off
end

xlabel("Time Step")

nexttile(1)
legend(["Input" "Anomalous"])