Analyse van de FAIRMODE CT6 methode voor PM2.5 fijnstofmetingen
We analyseren in een aantal sessies de voorbeeld fijnstof kalibratie source code, die we ontvingen van het RIVM: sensorFusion_21Jan2021.R en interpolate.R
In deze eerste sessie kijken we naar het overall verhaal;in volgende sessies kijken we naar details
Allereerst wat definitiekwesties:
kalibratie: de nauwkeurigheid (van meetinstrumenten) bepalen aan de hand van een standaard en evt het meetinstrument ijken aan de hand van die standaard
Het huidige algoritme voorziet van een set van de de goedkope fijnstofsensoren ruwe meetresulaten
van een betrouwbaarder waarde, in de wandelgangen ook wel de “gekalibreerde waarde” genoemd.
Strict genomen is dit niet correct, want we laten statistiek los op een hele set en niet op op individuele
sensoren.
Om de verwarring niet groter te maken blijven we we de term “kalibreren” gebruiken in het vervolg, waar
we eigenlijk “groeps correctie” bedoelen.
We kijken naar het huidige algoritme as is, met twee doelen:
- transparatie beorderen door begrip
- ideeen opdoen voor alternatieve correctie methodes
Gebruikte afkortingen:
CSS = Citizen Science Station
RRS = RIVM Reference Station
NA = Niet Available , er is voor dat moment geen meting.
De CSS meetwaardes
-wijken af van de RRS metingen, (omdat het goedkope sensoren ipv professionele aparatuur betreft)
-wijken af van nabij gelegen CSS van dezelfde fabrikant / andere fabrikant
-hebben afhankelijkheden van externe factoren (bv relatieve vochtigheid)
-hebben mogelijk lange termijn drift
Hoe plak je nu op al die onbetrouwbaarheid iets betrouwbaarders ?
De kern van het algoritme is een niet triviale statische resample methode, genaamd bootstrapping
Het voordeel van deze methode is dat je een schatting kunt maken hoe groot de fout is die je maakt
door met onbetrouwbare meetdata te werken.
Een kleine terzijde is dat in de huidige RIVM interface de berekende kal waardes kunnen worden
opgevraagd, maar niet de bijbehorende statistische onzekerheid.
Elk uur wordt onderstaand algorithme uitgevoerd, waarbij deze interpretatie geheel voor onze
eigen rekening komt, en wellicht nog niet 100% juist is.
Main flow
-
lees voor alle Nederlandse CSS-sensors de PM2.5 metingen voor dit specifieke uur
In de rivm operationele praktijk worden alleen sensors van het type SDS011, gekalibreerd.
(99% van alle in Nedeland aanwezige CSS zijn van het type SDS011)In de voorbeeldcode wordt dat type onderscheid niet gemaakt en worden bv ook PMS sensoren meegenomen
NB hieruit volgt dus dat de voorbeeldcode zelf nu niet draait op de operationele rivm server-
Verwijder PM2.5 waardes die groter zijn dan 500 ug/m^3
-
Als er meerdere metingen op een enkele locatie zijn, neem de gemiddelde waarde ,waarbij
vooraf NA’s verwijderd zijn
-
-
Lees de RRS sensor data van alle RRS in Nederland voor hetzelfde uur
-
Verwijder PM2.5 waardes kleiner dan 0 (duidt op een RSS error)
-
Als er meerdere RRS op een enkele locatie staan, neem het gemiddelde en gooi de NA’s weg
-
-
Maak voor elke RRS een groep CSS als voldaan is aan de eis dat:
“er tenminste 2 CSS zijn,binnen een straal van 7.5 km van de RRS.”
Als er niet voldaan wordt aan die eis, dan krijgt de betreffende RRS geen “CSS-buurgroep”
en doet niet mee met het kalibreren van de CSS.Als er meer dan 10 CSS voldoen, neem dan de 10 CSS die het meest nabijgelegen zijn bij betreffende RRS
Een enkele CSS kan evt in meerdere RSS groepen zitten.We gaan nu verder met die RRS die wel een CSS buurgroep hebben
RRS die geen CSS buurgroep hebben, komen in de Bootstrap procedure in 7.4 aan bod -
Gooi alle rijen met NAs weg van de betreffende uur data, zowel voor CSS als RRS
-
Outlier analysis : Gooi alle rijen weg uit CSS meetwaardes die buiten de landelijke the 5% and 95%
percentiles liggen. Hiermee worden voor de kalibratie dus de extreem lage en de extreem
hoge CSS waardes buiten beschouwing gelaten. -
Initialiseer de calibration data frames
-
Kalibratie factor per RRS sensor met 0
-
Kalibrated CSS sensor waardes met de ruwe meetwaardes van de sensoren binnen de groep
-
-
Doe de Bootstrap procedure een aantal keer (bv 100).
De Bootstraprocedure heeft als doel om de onzekerheid in de berekende kalibratie te schatten.
De Bootstraprocedure wordt hieronder apart uitgewerkt -
Vanuit de bootstrap procesdure kunneen we door van de herhaalde de sampeling de 5% en 95% grenzen te
bepalen, een confidence inteval bepalen , zowel voor de RRS calibratie factor , als voor de CSS
gekalibreerde waardes. -
Save RRS calibration factors en calibrated CSS PM2.5 meetdata in csv files
-
Creeer een plaatje van Nederland met in kleur gecodeerd de calibration factor voor elke RRS
Bootstrap procedure (in een run bv 100 keer uitgevoerd)
7.1. Sample alle CSS data van het betreffende uur (uit stap 5)
7.2. Sample een gedeelte (80%) van de RRS data met teruglegging voor de volgende loop
7.3. For each RRS (met buur-groep) loop door de sensor group:
1. Neem de fijnstof uurwaarde waarde van he RRS als reference value
2. Kalibreer dusdanig dat het gemiddelde van de metingen binnen de groep gelijk is aan de bijbehorden RRS.
Deze loop_ bootcalibratiefactor is gelijk aan: RRSwaarde / (average of the CSS data for this group)
De aanname is dat er geen offset is tussen het gemiddelde van de CSS groep en het bijbehorende RRS,
maar alleen een factor !
Bjvoorbeeld : als de RSS 10 geeft en de gemiddelde bijbehorende CSSgroup 15, dan is de
vermenigvuldigingsfactor voor de CSS: 10/15
NB Er wordt in dit voorbeeld dus niet 5 afgetrokken van de CSS groep
3. Sla de berekende kalibratiefactor per RRS op
7.4. …
7.5. Voor elke CSS in the sample vermenigvuldig de gemeten waarde met een factor die verkregens wordt door gewogen interpolatie van alle nederlandse RRS uur waardes met de kalibratiefactor uit 7.3.2.
De weging is zo dat de correctie afneemt met de afstand, dwz hoe verder de CSS af ligt van de RRS hoe minder de relatieve invloed van de RRS op de CSS kalibratiecorrectie is.
In grote lijnen de formule om de uurcorrectiefactor voor een CSS te berekenen
Een CSSi heeft vaste afstand Dij tot een RRSj
Voor elk uur
heeft RRSj een korrectie factor Rj
heeft een CSSi een PM25 meetwaarde PM25ruw
dan wordt de Correctiefactor voor CSSi een op afstand Dij gewogen som van een functie op de RCj
en is de gekalibreerde fijnstof PM25 waarde CSSi = PM25ruw * CNi
De weging met afstand Dij vindt plaats met de IDW formule, zie bv:
Al met al een ingewikkelde procedure, statistisch verantwoord, maar het is mij niet helemaal
duidelijk hoe de vertaling te maken naar de fysica.
De functie formule in de code bevat een aantal parameters. Aan elke van die parameters kun je draaien.
Een volgende keer kijken we steeksproefgewijs wat het effect daarvan is op het eind resultaat.
Intussen denken we of we een eenvoudiger model kunnen bedenken.
. Lasso techniek zit dacht ik standaard in de Python math lib. Mi is de groepscorrectie methode veelbelovend. Maar helpt alleen bij grotere aantallen in de buurt bivakkerende meetkits :-(.