Exemple interactif

Polluez!

// Code de la classe Track to track Kalman Filter

/** Class representing a track to track Kalman filter and accepting an arbitrary number of tracks*/
class T2TKFilter {
    tracks;
    n_tracks;
    w;

    /**
     * Initialize class
     * All parameters arrays must be of equal length
     * @param {Array<number>} init_estimates - Initial estimates, array of initial estimates, one value per track
     * @param {Array<number>} init_variances - Initial uncertainties, array of initial estimates, one value per track
     * @param {Array<number>} measure_variances - Initial measure variances, array of initial estimates, one value per track
     * @param {number} process_noise - Process noise of modeled phenomenon, 1 value for the whole filter
     * */
    constructor(init_estimates, init_variances, measure_variances, process_noise) {
        this.n_tracks = init_estimates.length;
        this.tracks = [];
        for (let i = 0; i < this.n_tracks; i++) {
            this.tracks.push({
                xp: init_estimates[i], // previous x measure estimate
                pp: init_variances[i], // previous p variance estimate
                r: measure_variances[i],
                xn: NaN,
                pn: NaN
            });
        }

        this.w = process_noise;
    }


    /**
     * Prediction step of the Kalman filter
     * @param {Array<number>} values - Measured values, one value per track
     * @param {bool} treat_outliers - Remove outliers outside of tolerance
     * @param {number} tolerance - % (from n standard deviation ) of mahalanobis distance tolerated, only used if treat_outliers=true
     * */
    predict(values, treat_outliers = true, tolerance = 0.9) {       
        for (let i = 0; i < this.n_tracks; i++) {
            let z = values[i]
            let xp = this.tracks[i].xp;
            let pp = this.tracks[i].pp;
            let r = this.tracks[i].r;

            /**
             * If mahalanobis distance > tolerance
             * set variance to infinity, 
             * thus making xn = xp and pn = pp
            */
            if (treat_outliers && (Mahalanobis1D(xp, z, pp) > 1 / Math.sqrt(1 - tolerance))) {
                r = Infinity;
            }

            let Kn = pp / (pp + r); // Calculate kalman gain

            this.tracks[i].xn = xp + Kn * (z - xp);// calculate estimate
            this.tracks[i].pn = (1 - Kn) * pp;
        }
    }

    /**
     * Update step of the Kalman filter
     * @return {number} the predicted value
     * */
    update() {
        /* 
            This method return a fused estimates of each track.
            It also updates the system
        */

        for (var track of this.tracks) { // Update
            track.xp = track.xn;
            track.pp = track.pn + this.w;
        }

        //fuse all track estimates
        return this.tracks.reduce((acc, track) => 
          acc + track.xn / track.pn, 0) / this.tracks.reduce((acc, track) => acc + 1 / track.pn, 0);

    }

    /**
     * Shorthand for update then predict steps of the Kalman filter
     * @return {number} the predicted value
     * */
    predict_update(values) {
        this.predict(values);
        return this.update();
    }
}


/**
* Calculates a single dimension mahalanobis 
* distance between an an estimate, a measure and its variance
* @return {number} the distance in number of standard deviations
* */
function Mahalanobis1D(x, z, p) {
    return Math.sqrt((x - z) ** 2 / p)
}

Problème

Comment faire pour suivre la qualité de l'air réelle d'une localisation alors que nous ne disposons que de capteurs très peu fiables, installés dans des conditions non idéales ?

Filtre de Kalman

Présentation

L'intuition du filtre de Kalman est élémentaire :

1. On émet une supposition sur une donnée et son incertitude
2. On regarde la donnée mesurée
3. On émet les supposition pour la prochaine mesure que l'on ajuste en fonction de la différence entre l'incertitude supposée et le résultat mesuré.
4. On répète .2.

Ainsi, un filtre rudimentaire ne suivant qu'une seule sortie capteur peut être schématisé ainsi :

  %%{
  init: {
    'theme': 'base',
    'themeVariables': {
      'primaryColor': '#ececed',
      'primaryTextColor': '#404040',
      'primaryBorderColor': '#000',
      'lineColor': '#F8B229',
      'secondaryColor': '#006100',
      'tertiaryColor': '#fff'
    }
  }
}%%
flowchart TD;
subgraph Données;
Z;
I;
end;

subgraph Filtre;
Z("Mesure capteur") -->|"$$\LARGE{z_n}$$"| B("Actualiser<br>
    $$\begin{gather} \hat{x}_{n,n}  = {x}_{n,n-1} + K_n (z_n - \hat{x}_{n,n}) \\ K_n = \frac{p_{n,n-1}}{p_{n,n-1}+r_n} \\ p_n = (1-K_n)p_{n, n-1}\end{gather}$$");

I("Initialiser<br>$$\begin{gather}\hat{x}\text{ l'estimation} \\ p \text{ l'incertitude} \\  q \text{ le bruit du phénomène} \\ r \text{ la variance du capteur}\end{gather}$$") --> |"$$\large\hat{x}_0 \text{ ,  } p_0 \text{ et  } q_0$$"| C;
I --> |"$$\LARGE{r}$$"|B;
B --> C("Prédire <br> $$\begin{gather} \hat{x}_{n+1,n} = \hat{x}_{n,n} \\ p_{n+1,n} = p_{n,n}+q_n \end{gather}$$");
C --> |"$$\large\hat{x}_{n,n-1} \text{ et } p_{n,n-1}$$"| B;
end;
C --o D(("$$\LARGE\hat{x}_{n,n}$$"));

En pratique

Par défaut, le modèle suppose que la variable qu'on cherche à mesurer est constante dans le temps. C'est rarement le cas dans la réalité. Un bruit "process noise" a donc été ajouté : il permet au filtre d'accepter une certaine variance autour de l'estimation. C'est le choix de ce paramètre qui va décider le plus (avec la variance du capteur) de l'efficacité du filtre.

En faisant varier le bruit du process, on modifie la vitesse avec laquelle le filtre considère qu'un changement de valeur constitue une innovation réelle. Bien sûr, pour filtrer le bruit d'un seul capteur, cet outil n'est pas spécialement efficient.

Ce type de filtre se prête davantage au suivi de systèmes complexes comme un modèle CMAQ de qualité de l'air rectifié avec des capteurs (Djalalova et al., 2015).

Fusion de senseurs

Il est possible de fusionner l'échantillonage de plusieurs senseurs avec les filtres de Kalman. La fusion peut avoir lieu au moment de la mesure ou de la prédiction. La seconde méthode est utilisée ici.

  %%{
  init: {
    'theme': 'base',
    'themeVariables': {
      'primaryColor': '#ececed',
      'primaryTextColor': '#404040',
      'primaryBorderColor': '#000',
      'lineColor': '#F8B229',
      'secondaryColor': '#006100',
      'tertiaryColor': '#fff'
    }
  }
}%%
flowchart LR;

subgraph Données
Z1
A1

Z2
A2
end

Z1("Capteur A") --> B1 ;
A1("Initialisation") --> C1;
A2("Initialisation") --> C2;
Z2("Capteur B") --> B2;

subgraph "filtre A"
C1("Prédiction") --> B1("Actualisation")
B1 --> C1
end

subgraph "filtre B"
C2("Prédiction") --> B2("Actualisation")
B2 --> C2

end

C1 --> |"$$\LARGE \begin{gather}\hat{x}_n^{(a)} \\ p_n^{(a)}\end{gather}$$ ____"|F("Fusion de données<br>$$\hat{x}_n = {{\hat{x}_n^{(a)}p_n^{(a) -1}+\hat{x}_n^{(b)}p_n^{(b) -1}} \over {p_n^{(a) -1}+p_n^{(b) -1}}} $$")
C2 --> |"$$\LARGE \begin{gather}\hat{x}_n^{(b)} \\ p_n^{(b)}\end{gather}$$ ____"|F
F --o X(("$$\LARGE{\hat{x}_n}$$"))

Voici un exemple de fusion issu des mesures PM10 sur 24h de 2 microcapteurs installés à Thorigné-Fouillard tirés de la plateforme de science collaborative sensor community.