深層学習を使用した時系列異常検出
この例では、シーケンス データまたは時系列データで異常を検出する方法を示します。
シーケンス データまたは時系列データの集合で異常または異常な領域を検出するには、自己符号化器を使用できます。自己符号化器は、入力を低次元空間に変換 (符号化ステップ) し、低次元表現から入力を再構成 (復号化ステップ) することによって、入力を複製するように学習させるモデルの一種です。自己符号化器の学習に、ラベル付けされたデータは必要ありません。
自己符号化器自体は異常を検出しません。代表的なデータのみを使用して自己符号化器に学習させると、代表的なデータから学習した特徴のみを使用して入力データを再構成できるモデルが得られます。自己符号化器を使用して観測値が異常かどうかをチェックするには、観測値をそのネットワークに入力し、元の観測値と再構成された観測値の間の誤差を測定します。元の観測値と再構成された観測値の間の誤差が大きいことは、自己符号化器の学習に使用されたデータを代表するものではない特徴が元の観測値に含まれており、異常があることを示しています。元のシーケンスと再構成されたシーケンスの間の要素単位で誤差を観測することで、異常の局所領域を特定できます。
以下のイメージは、異常な領域を強調表示したシーケンスの例を示しています。
この例では、波形データ セットを使用します。このデータ セットには、3 つのチャネルの異なる波長の合成生成波形が 2,000 個含まれています。
学習データの読み込み
WaveformData.mat
から波形データ セットを読み込みます。観測値は numTimeSteps
行 numChannels
列の配列です。ここで、numtimeSteps
と numChannels
はそれぞれシーケンスのタイム ステップ数とチャネル数です。
load WaveformData
最初のいくつかのシーケンスのサイズを表示します。
data(1:5)
ans=5×1 cell array
{103×3 double}
{136×3 double}
{140×3 double}
{124×3 double}
{127×3 double}
チャネル数を表示します。ネットワークに学習させるには、各シーケンスのチャネル数が同じでなければなりません。
numChannels = size(data{1},2)
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 回ダウンサンプリングします。
numDownsamples = 2;
シーケンスを 2^numDownsamples
の最も近い倍数に切り捨てます。ネットワーク入力層の最小シーケンス長を計算できるように、シーケンス長を含むベクトルも作成します。
sequenceLengths = zeros(1,numel(XTrain)); for n = 1:numel(XTrain) X = XTrain{n}; cropping = mod(size(X,1), 2^numDownsamples); X(end-cropping+1:end,:) = []; XTrain{n} = X; sequenceLengths(n) = size(X,1); end
同じ手順で、検証データを切り捨てます。
for n = 1:numel(XValidation) X = XValidation{n}; cropping = mod(size(X,1),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")];
ネットワークを対話的に表示または編集するには、ディープ ネットワーク デザイナーを使用できます。
deepNetworkDesigner(layers)
学習オプションの指定
学習オプションを指定します。
Adam ソルバーを使用して学習させます。
学習を 120 エポック行います。
すべてのエポックでデータをシャッフルします。
検証データを使用してネットワークを検証します。シーケンスを入力とターゲットの両方として指定します。
学習の進行状況をプロットに表示します。
詳細出力を非表示にします。
options = trainingOptions("adam", ... MaxEpochs=120, ... Shuffle="every-epoch", ... ValidationData={XValidation,XValidation}, ... Verbose=0, ... Plots="training-progress");
ネットワークの学習
関数trainnet
を使用してニューラル ネットワークに学習させます。自己符号化器に学習させる場合、入力とターゲットは同じです。回帰の場合は、平均二乗誤差損失を使用します。既定では、関数 trainnet
は利用可能な GPU がある場合にそれを使用します。GPU を使用するには、Parallel Computing Toolbox™ ライセンスとサポートされている GPU デバイスが必要です。サポートされているデバイスについては、GPU 計算の要件 (Parallel Computing Toolbox)を参照してください。そうでない場合、関数は CPU を使用します。実行環境を指定するには、ExecutionEnvironment
学習オプションを使用します。
net = trainnet(XTrain,XTrain,layers,"mse",options);
ネットワークのテスト
検証データを使用してネットワークをテストします。関数minibatchpredict
を使用して予測を行います。既定では、関数 minibatchpredict
は利用可能な GPU がある場合にそれを使用します。
YValidation = minibatchpredict(net,XValidation);
各テスト シーケンスについて、予測とターゲットの間の平方根平均二乗誤差 (RMSE) を計算します。ターゲット シーケンスの長さを参照に使用し、予測シーケンス内のパディング値をすべて無視します。
for n = 1:numel(XValidation) T = XValidation{n}; sequenceLength = size(T,1); Y = YValidation(1:sequenceLength,:,n); err(n) = rmse(Y,T,"all"); end
ヒストグラムで RMSE 値を可視化します。
figure histogram(err) xlabel("Root Mean Square Error (RMSE)") ylabel("Frequency") title("Representative Samples")
最大 RMSE を異常検出のベースラインとして使用できます。検証データから最大 RMSE を決定します。
RMSEbaseline = max(err)
RMSEbaseline = single
0.6140
異常シーケンスの特定
異常な領域を含むように検証シーケンスの一部を手動で編集して、一連の新しいデータを作成します。
検証データのコピーを作成します。
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 = minibatchpredict(net,XNew);
予測ごとに、入力シーケンスと再構築シーケンスの間の RMSE を計算します。
errNew = zeros(numel(XNew),1); for n = 1:numel(XNew) T = XNew{n}; sequenceLength = size(T,1); Y = YNew(1:sequenceLength,:,n); errNew(n) = rmse(Y,T,"all"); end
プロットで RMSE 値を可視化します。
figure histogram(errNew) xlabel("Root Mean Square Error (RMSE)") ylabel("Frequency") title("New Samples") hold on xline(RMSEbaseline,"r--") legend(["Data" "Baseline RMSE"])
最大 RMSE 値をもつ上位 10 個のシーケンスを特定します。
[~,idxTop] = sort(errNew,"descend");
idxTop(1:10)
ans = 10×1
1
88
62
39
14
58
31
56
75
83
最大 RMSE 値をもつシーケンスとその再構成をプロットで可視化します。
X = XNew{idxTop(1)}; sequenceLength = size(X,1); Y = YNew(1:sequenceLength,:,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"])
異常な領域の特定
シーケンス内で異常な領域を検出するには、入力シーケンスと再構成シーケンスの間の RMSE を求め、しきい値を超える誤差がある領域を強調表示します。
入力シーケンスと再構成シーケンスの間の誤差を計算します。
RMSE = rmse(Y,X,2);
タイム ステップのウィンドウ サイズを 7 に設定します。検証データを使用して特定された最大誤差値を RMSE の値が少なくとも 10% 上回るタイム ステップが含まれるウィンドウを特定します。
windowSize = 7; thr = 1.1*RMSEbaseline; idxAnomaly = false(1,size(X,1)); for t = 1:(size(X,1) - windowSize + 1) idxWindow = t:(t + windowSize - 1); if all(RMSE(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,1)); XAnomalous(idxAnomaly) = X(idxAnomaly,i); plot(XAnomalous,"r",LineWidth=3) hold off end xlabel("Time Step") nexttile(1) legend(["Input" "Anomalous"])
強調表示された領域は、誤差値が最大誤差値より少なくとも 10% 高いタイム ステップのウィンドウを示しています。
参考
trainnet
| trainingOptions
| dlnetwork
| sequenceInputLayer
| convolution1dLayer
| transposedConv1dLayer