Detekce anomálií v datech o znečištění ovzduší

Chytré lampy na Karlínském náměstí poskytují data o znečištění ovzduší polétavým prachem PM₁₀. V této práci v nich detekuji anomálie pomocí algoritmů strojového učení pro predikci časových řad.

Otevřená data přináší velký prostor pro aplikaci metod strojového učení a pro zlepšení kvality života nebo přírodního prostředí. Tato práce se zabývá detekcí anomálií v datech o polétavém prachu z chytrých lamp, které jsou poskytovány hlavním městem Praha jako otevřená data.

Pevné částice (nebo polétavý prach, anglicky particulate matter, PM) jsou tuhé nebo kapalné částice v ovzduší o velikosti v rozsahu 1 nm až 100 µm. Data použitá v tomto projektu obsahují PM₁₀ (částice menší než 10 µm) a PM_2,5 (částice menší než 2,5 µm). PM_2,5 je tedy podmnožina měření PM₁₀, a proto budu hledat anomálie pouze v PM₁₀.

V této práci chci prozkoumat přístup k měření pomocí predikce časové řady následovanou prahovnáním. Chci, aby se zvolený model naučil pravidelnosti v datech a pokud se stane něco neočekávaného, tak to pomocí prahování dokáži odhalit. Jako algoritmy strojového učení jsem vyzkoušel lineární regresi a LSTM rekurentní neuronovou síť. Ukázalo se, že lineární regrese pro tento problém postačuje, protože je schopna časovou řadu koncentrace PM₁₀ nejlépe předpovídat.

Tato zpráva je rozdělena do několika částí, které se postupně zabývají explorací dat, předzpracovním, metrikami, selekcí příznaků, porovnáním modelů a testováním výsledného algoritmu.

Explorace data

Pražská datová platforma Golemio poskytuje data z pilotního provozu Senzorické sítě veřejného osvětlení, v rámci kterého bylo nainstalováno 92 chytrých pouličních LED lamp v blízkosti Karlínského náměstí na Praze 8. Některé z těchto lamp mají senzory pro měření a sběr dat o hluku, prašnosti a množství dalších polutantů.

Datové zdroje

V katalogu datové platformy Golemio jsou k dispozici dva datové CSV soubory: jeden za 2. pololetí roku 2018 a druhý za zatím uplynulou část 1. pololetí roku 2019. Data z roku 2018 obsahují celkem 473529 záznamů. Z toho jsou dostupné odečty PM₁₀ pouze v 129807 případech. V souboru za roku 2019 je celkem 226672 záznamů a z toho 62518 měření PM₁₀.

Každý záznam se skládá z následujících atributů:

V tabulce ppb znamená anglicky parts per bilion.

Datová sada by měla obsahovat odečty z 43 senzorů, které jsou rozmístěny na vybraných lampách, ale ve skutečnosti jsou k dispozici data pouze z 22 senzorů.

Z těchto 22 senzorů měří podle dostupných dat polétavý prach PM₁₀ pouze 6 senzorů a z nich nejvíce měření má stanice s identifikátorem: y7e4onsytkb3ydonflz5kcbcigkh5ulo.

Trénovací a testovací data

Vezmeme-li pouze data z této lampy a ponecháme si atributy o3, pm2p5, no2, so2 a pm10, tak dostaneme pět různých dimenzí příznaků, které obsahují 32872 záznamů v roce 2018 a 12709 záznamů v roce 2019. To je dohromady 45581 záznamů za oba roky. V této práci pracuji s časovou řadou, a proto zachovám jako index datum a čas měření (atribut starttimestamp).

Pro potřebu verifikace natrénovaného modelu se nabízí rozdělit data na trénovací a testovací množinu podle roku měření. Takové rozdělení má 72,1% záznamů v trénovací množině a 27,9% záznamů v množině testovací.

Vizualizace

V následující části jsem se zaměřil na porozumění datům o znečištění ovzduší skrz vizualizace. Vizualizace nezpracovaných časových řad jednotlivých atributů ukazuje, že se v datech zřejmě anomálie vyskytují. Stačí pouze data vhodně předzpracovat a vybrat algoritmy strojového učení, které dokážou anomálie nejlépe detekovat. Dále je vidět, že data na první pohled neobsahují žádné sezóního chování.

Rozložení absolutních velikostí PM₁₀ odhaluje boxplot. Mnoho bodů se nachází nad 3. kvartilem. Průměr příznaku je 14,093698 s relativně malou standardní odchylkou hodnoty 18,931824. Samozřejmě všechny tyto body nemusí být anomálie, protože mnou navržená detekce záleží na předchozích měřeních a ne pouze na absolutní velikosti.

Lineární tvar lag grafu ukazuje, že vhodnou metodou pro predikci časové řady bude nějaký autoregresní model (hodnotu předpovídá na základě několika předchozích pozorování). Graf odhaluje zhruba 6 anomálních měření, které leží mimo diagonálu.

Pokud se v datech nacházejí lineární závislosti je důležité je odhalit, protože následně stačí použít jednoduchých lineárních modelů a je možné dosáhnout vysoké efektivnosti i přesnosti.

Lineární závislosti mezi atributy odhaluje korelační diagram (anglicky scatter plot). Protože hledám anomálie v PM₁₀, tak mě zajímá korelace příznaku pm10 s ostatními. Za povšimnutí stojí, že anomálie v PM₁₀ se vyskytují pouze, pokud je nízká hladina SO₂ a O₃.

Každopádně tento vztah zřejmě není lineární, protože korelační koeficient je v obou případech blízký nule (0,03 resp. -0,07) v porovnání s největším korelační koeficient, který má s příznakem pm2p5, a to 0,53:

Přesto je z korelačního grafu PM₁₀ a PM_2,5 vidět, že některé anomálie nastávají, jak při zvýšené hladině PM_2,5, tak při jejím normálním stavu.

Intervaly měření

Podle specifikace je vyčítací frekvence 15 minut z brány do platformy. Data ze souborů pravidelné intervaly neobsahují. Průměrně je dostupné nové měření asi každých 7 minut a 52 sekund se standardní odchylkou 18 minut a 20 sekund. Výpadky ale mohou systém odstavit i na delší dobu. V našich datech je nejdelší prodleva přes 1 den a 15 hodin. Následují graf ukazuje histogram intervalů měření.

Log–log plot

Poslední vizualizační pomůckou, která je vhodná pro hledání anomálií z hlediska jejich absolutní velikosti je log–log graf. Především v jeho pravé dolní části je vidět, že data obsahují některé abnormálně velké hodnoty a bude úkolem námi naučených klasifikátorů tyto hodnoty najít i vzhledem k časové závislosti na bodech předchozích.

Předzpracování dat

Původní data tedy obsahují problémy, které by mohly zamezit správnému natrénování algoritmů strojového učení. Za prvé nejsou data samplovány v rovnoměrných intervalech, a proto je v rámci předzpracování data přesamplujeme každých 15 minut (vyčítací frekvence z platformy). K doplnění chybějících hodnot jsem použil metodu fill forward, která doplní vždy předchozí hodnotu.

Za druhé data musím data přetransformovat do podoby vhodné pro předpovídání časových řad, abych potom mohl detekovat anomálie. Takže ze všech příznaků měření v čase \(t\) se budu snažit předpovědět hodnotu pm10 v čase \(t + 1\).

Poté data rozdělím na trénovací a validační množinu. Validační množina slouží k výběru modelu a jeho hyperparametrů. Konzervativně jsem se rozhodnul, že data rozdělím na půlku, aby odhad byl co nejpřesnější.

Nakonec použiji standardní škálování, aby různé příznaky byly zhruba stejných rozsahů. Každý příznak trénovací množiny bude mít nulový průměr a jednotkovou standardní odchylku. Takto předzpracovaná data umožní správnou selekci příznaků pomocí linearní regrese s L1 regularizací, kterou se zabývám v další části.

Metrika úspěšnosti

Pro měření úspěšnosti predikce časové řady je vhodná odmocnina střední kvadratické chyby (anglicky root mean squared error neboli RMSE), kterou zvolený model minimalizuje, ale není žádoucí dosáhnout nulové chyby, protože potom by nebylo možné detekovat anomálie. RMSE je definována takto: \[\mathrm{RMSE} = \sqrt{\frac{1}{N} \sum^{N}_{i = 1}(\hat{y_i} - y_i)^2},\] kde \(N\) je počet měření, \(\hat{y_i}\) je predikce modelu a \(y_i\) je skutečně naměřená hodnota. Chci aby se model naučil pouze pravidelnosti v datech a nikoliv se přeučit, aby si dokázal zapamatovat anomálie.

V měření úspěšnosti detekce anomálií je problém absence označení anomálií v datech. Detekce tedy musí být kontrolovány člověkem (nejlépe doménovým expertem).

Selekce příznaků

Selekce podmnožiny příznaků může být přínosná pro snížení komunikační zátěže při online monitorování, kdy není potřeba skrz API dotazovat všechny příznaky. Vhodné příznaky pro lineární modely můžeme vybrat lineární regrese s L1 regularizací (tzv. lasso). L1 regularizace totiž koeficienty u nevhodných příznaků drží blízko nule.

Z výše zmíněných důvodů jsem zkusil natrénovat lasso pro 23 různých hodnot regularizačního parametru alpha na logaritmické škále mezi 10^-10 a 10¹. Výsledné grafy níže ukazují, že důležitý je pouze příznak pm10. Koeficienty ostatních příznaků jsou blízké nule, pokud je RMSE nízká. Pro testované lineární modely stačí tedy uvažovat pouze tento příznak.

Porovnání modelů

V této práci detekuji anomálie pomocí předpovědi časové řady, a to stejně jako v článku Time Series Anomaly Detection. Při detekování anomálií v časových řadách se používá metoda, kdy je daná řada modelem předpověděna dopředu (například pomocí rekurentní neuronové sítě) a následně je porovnána se skutečným měření. Z tohoto porovnání jsou určeny anomálie (například pomocí prahování a dalších pravidel).

Jako vhodné modely jsem zvolil následující tři:

1. základní model (anglicky baseline),
2. lineární regresní model a
3. rekurentní neuronovou síť (konkrétně dnes nejpopulárnější Long Short-Term Memory).

Základní model

Základní model poskytuje odrazový můstek, ke kterému budeme moci vztahovat výsledky ostatních modelů. Navíc je jednoduchý na implementaci a dosahuje relativně dobrého RMSE: 0,6492.

Lineární regrese

Druhým modelem je lineární regrese bez regularizace, protože selekce příznaků ukázala, že stačí příznak pm10. Pokud tento model předpovídá pouze z bezprostředně předcházející hodnoty, jeho úspěšnost v RMSE je 0,6402.

Lepší je nechat lineární regresi předpovídat z více předcházejících měření. Graf níže ukazuje výsledek pro velikosti okna předcházejících měření velikosti 1 až 100. Jako nejlepší se ukazuje velikost okna velikosti 2, kde je RMSE rovna 0,6380. Následně se RMSE zhorší a konverguje k hodnotě zhruba 0,71.

Long Short-Term Memory (LSTM)

Nejpokročilejším prediktivním modelem, který na data použiji, je rekurentní neuronová síť, a to konkrétně architektura LSTM. Rekurentní sítě se ukázaly být velice efektivní i v aproximování nelineárních vztahů v datech, proto nyní zachovám všech pět příznaků jako vstup rekurentní sítě. To ukazuje blog The Unreasonable Effectiveness of Recurrent Neural Networks.

Konkrétní architekturu jsem zvolil jednovrstvou LSTM s dimenzí skrytého vektoru 64. Výstupem LSTM musí být pouze jedno číslo, které predikuje měření PM₁₀, a proto je na vlastní buňce LSTM ještě lineární vrstva, která produkuje predikci PM₁₀ v následujícím čase.

Tuto LSTM jsem trénoval pomocí MSE ztrátové funkce: \[\mathrm{MSE} = \frac{1}{N} \sum^{N}_{i = 1}(\hat{y_i} - y_i)^2,\] kde \(N\) je počet měření, \(\hat{y_i}\) je predikce modelu a \(y_i\) je skutečně naměřená hodnota, na vstupních sekvencích délky 100 po celkem 50 epoch. Hodnoty ztátové funkce na trénovacím a validačním setu ukazuje graf níže. Velký rozdíl mezi trénovací a validační ztrátovou funkcí nemusí v tomto regresním problému představovat přeučení, protože ve validačním setu je jedna velká anomálie, která může hodnotu ztrátové funkce výrazně zhoršit. Výsledné RMSE na validační setu je 1,0588.

Výsledky

Následující tabulka srovnává zvolené predikční modely z hlediska metriky RMSE měřené na validačním setu:

Nejlepších výsledků dosahuje lineární regresní model, který má na vstupu dvě měření z historie, a proto tento model použijeme pro detekci anomálií.

Detekce anomálií

Nyní po výběru vhodného prediktivního modelu mohu detekovat anomálie postupem stejným jako v článku Time Series Anomaly Detection. Postup je následující:

1. lineární regresní model s oknem 2 predikuje další měření PM₁₀,
2. predikovanou hodnotu porovnám se skutečně naměřenou hodnotou (vypočítám rozdíl v absolutní hodnotě, \(\mid\hat{y}- y\mid\)),
3. pomocí prahování určím, zda se jedná o anomálii.

Histogram absolutní hodnoty rozdílu predikce \(\hat{y}\) a skutečné hodnoty \(y\) je na obrázku níže. Většina hodnot je malé velikosti, protože pochází z předvídatelné distribuce, kterou dokáže lineární regrese zachytit.

Teď bych potřeboval doménového experta, který by dokázal určit, které rozdíly jsou abnormální a které přirozené. Takový expert ale není k dispozici, a tak zvolím jako hodnotu pro pahování hodnotu okraje prvního binu z histogramu výše: 3,9617.

Takto zvolený práh odhalil ve validačních datech celkem 9 anomálií, které jsou zobrazeny v obrázku níže.

Testování

Postup pro detekci anomálii je vyvinutý a stačí jej otestovat na testovacích datech. Jak ukazuje graf níže, jsou tyto data daleko složitější a obsahují zřejmě daleko více anomálií. Pro trénování jsem použil všechna data z roku 2018, která předzpracuji stejně jako v předchozích pokusech a následně aplikuji postup popsaný v předchozí části.

Podle histogramu níže můj algoritmus opravdu správně detekoval daleko více anomálií než ve validačním setu (pro testovací množinu je hodnota prahu: 6,8293).

Pomocí stanoveného prahu jsem detekoval celkem 113 anomálií, které jsou zobrazené na grafu níže. Veškeré abnormálně vysoké špičky byly detekovány správně. Také je vidět, že některé špičky detekovány nebyly, to může být tím, že nárůst je tak pozvolný a za tak dlouhý časový okamžik, že se vlastně o anomálii nejedná, ale je to předpovídatelné chování.

Obrázek níže zobrazuje 9 náhodně vybraných anomálií z celkem 113 detekovaných.

Závěr

V této práci jsem se zabýval detekcí anomálií v otevřených datech z chytrých lamp o polétavém prachu PM₁₀. Jako nejlepší se ukázalo předpovídat časovou řadu pomocí lineární regrese (LSTM bylo horší), která na vstup dostane dvě poslední měření. Následně jsem předpověděnou hodnotu porovnal se skutečnou. Posledním krokem je prahování, které určí zda hodnoto je z hlediska modelu normální nebo ne.

Na testovacích datech se ukázalo, že tento postup je schopen detekovat většinu anomálií a také bere v úvahu momentální trendy v datech pokud trvají delší dobu a nestanou se náhlé změny.