A propos de la température (1) : analyse brute (de pomme)

1. Analyse brute (de pomme)
2. Découpage en 4 bandes
3. Variation quotidienne de la température : modèle en modulation d'amplitude
4. Variation quotidienne de la température : ajustement style modulation de fréquence
5. Variation quotidienne de la température : et pourquoi pas un modèle en modulation du rapport cyclique jour/nuit?
6. Synthèse en modèle simplifié de la température (pour des usages bateaux)

---

Non, non, ce blog n'est pas mort! Il dormait profondément, et il est grand temps de le réveiller un peu.

Honnêtement, mener en parallèle le "faire" et le "dire" sur plusieurs fronts à la fois, et sans que ça tourne au superficiel ou à la répétition facile, c'pas évident! D'aucuns diront que c'est normal car un homme / une femme (rayez la mention que vous préférez) ne sait pas faire deux choses à la fois. Je vous laisse jauger ça à l'aune de la phrase qui précède ;-)

Un dernier petit mot de circonstance.

En cette époque étrange au pays de Descartes, réactiver ce blog est aussi "Charlie". Chacun à son niveau, défendre la liberté de "faire" et "dire" dans le respect de l'autre, faire briller "Les Lumières" dans l'obscurantisme, combattre l'entropie encore et toujours...

 Cet article est le premier d'une série à propos de la température extérieure.
Les motivations sont multiples :

• se donner un éclairage qualitatif et quantitatif vis à vis des informations de température qui sortent presque toujours dans la littérature technique tel le lapin du chapeau du magicien. Ou alors c'est copié betement sur le voisin.
• se donner les moyens de manipuler le plus proprement possible l'info qu'on a à disposition. Que ce soit pour extraire des caractéristiques en vue de les injecter dans un autre calcul, ou pour alimenter une simulation.
• (re)voir des fondamentaux de traitement du signal, en mettant la théorie au service de la pratique
• (re)voir des fondamentaux du logiciel scientifique scilab, qui est gratuit et magnifiquement utile/puissant en parallèle du calcul crayon-papier

En fait, je suis déjà tombé sur cette problématique lors du boum d'activité sur la future feuille de calcul de bilan hygrothermique mensualisé.

Je suis retombé cet été sur ce problème de "lapin magique" en buchant sur le puits canadien. En particulier après avoir constaté un peu de "olé-olé" dans les données d'entrée et hypothèses de calcul dans plusieurs papiers. Juste avant de tomber sur un article qui met carrément les pieds dans le plat!

Dans les 2 cas, il faut alimenter le biniou avec des données qui tiennent la route. Les variations de température à la station-météo de Honfleur ne sont pas les memes qu'au village machin-chose dans les hauteurs de Nice... Corollaire, il faut pouvoir jauger/analyser -un peu- si des données météo "sorties du chapeau" tiennent la route.

Exemple : Météo-France, ZE autorité nationale en la matière. Leurs relevés et cie sont payants, typés selon l'utilisation ciblée, et la tambouille n'est pas bien claire.
Autre exemple : récupération d'un relevé brut (amateur ou payant) sous forme synthétique. Genre température moyenne/min/max du jour ou du mois.
Ou encore, mon super-relevé météo maison à base de sonde USB ou de mini-station météo automatique.

Allez, c'parti 

1- Récupération et mise en forme des données

Les données récupérées sont un relevé des températures de 2015 toutes les 10 minutes d'un site dans le département de la Loire.

Le relevé présentait une quinzaine de trous de 6 à 30h. Les trous ont été comblés « au mieux » avec les données des jours précédent et suivant.

D'une année sur l'autre, bien évidemment, la température réelle n'a aucune raison de recoller parfaitement  ! Dans le relevé, on observe sans surprise une discontinuité entre 31 décembre 2015 et 1er janvier 2015 pris comme 1er janvier de « l'année suivante » . On redresse cette tendance supposée linéaire.

N=length(temp);
x=(0:N-1)';
// Nombre de samples par jour
x0=N/365;
// Segment de droite de pente équivalente à la droite passant par les 2
// extrémités du relevé, et passant par 0 au milieu
temp_tendance = (temp($)-temp(1))*x/(N-1) - (temp($)-temp(1))/2;
// on redresse autour du point milieu
temp_redressee = temp – temp_tendance;

// Dans la suite, on va travailler sur signal centré
temp_centree = temp-mean(temp);

La température annuelle moyenne est de 11,54₃°C, et n'est pas affectée par le redressement de la tendance « climatique ».

Comme on peut le voir, la température annuelle suit à peu près une sinusoïde, avec quelques accidents saisonniers, une sorte d'oscillation plus ou moins hebdomadaire, et enfin une oscillation journalière un peu chaotique.

2- Analyse spectrale de la température annuelle

On commence par mettre de côté la température annuelle moyenne, pour travailler sur un signal centré.

On applique au signal une transformée de Fourier (FFT), ainsi que l’estimateur de densité spectrale de puissance par la méthode de Welch.

Pour mémoire, le pas entre 2 points du spectre vaut l’inverse de la durée d’échantillonage. Autrement dit, pour N échantillons représentant 1 an, la résolution de la FFT vaut 1 « fois par an », soit (365*24*3600)^-1 Hz.

Au fait, la méthode de Welch, késako ? Appliquée correctement, c’est une façon de booster le « contraste » et la lisibilité de la FFT, sans tomber dans de la bidouille infâme. A la base, c’est une méthode d’analyse d’un signal bruité. Cette méthode permet typiquement de faire remonter des raies noyéees dans du bruit. Elle peut aussi faire ressortir des fondamentales cachées, c’est à dire absentes de la FFT mais exprimées implicitement par superposition constructive d’harmoniques. Elle rend aussi la résolution plus grossière, c’est à dire qu’elle rend les petits lobes plus pointus. On l’utilise ici comme une sorte d’autocorrélation, en ayant à l’esprit qu’elle donnera par construction d’autant moins d’information supplémentaire vis à vis de la FFT que la fréquence sera basse. Pourquoi? Comment? Pour avoir un résultat intéressant, on applique la méthode de Welch avec les paramètres classiques suivants : Fenêtre mobile de longueur K pris comme la puissance de 2 immédiatement inférieure au signal → K ≤ N < 2K ≤ 2N Recouvrement 50 % entre 2 fenêtres mobiles Fenêtre de Hamming, qui augment la résolution de x1,3 contre x1 pour la fonction porte. Pour mémoire, le rôle principale de ce fenêtrage est de solutionner les phénomènes de Gibbs en lissant les discontinuités du signal périodique sous-jacent qui apparaîtraient inévitablement avec une fonction porte.1 On allonge le signal à 50 répliques. Avec pour effet : Une exploitation des données plus complète/homogène, et moins biaisée grâce à la périodisée. Travailler sur une unique période aurait mécaniquement éjecté du calcul les N-K échantillons finaux et renforcé la contribution des échantillons entre K/2→K ! De plus, même avec le filtrage par Hamming, la fenêtre mobile travaille sur un signal tronqué/déformé, d’où un biais d’autant plus important que la fréquence est basse. La périodisée augmente le nombre de contributions décorrélées et atténue ce problème – au détriment du temps de calcul et des ressources machine. Entre la fenêtre de Hamming et la taille de la fenêtre mobile, la résolution passe de 1 à 1.3*N/K « fois par an. Vu les bornes sur K, elle est donc comprise entre 1,3 et 2,6 « fois par an ». Suite à cette perte de résolution : On se gardera d’observer ou conclure quoi que ce soit sur le Welch pour la variation annuelle, c’est à dire de la fréquence annuelle jusqu’à la fréquence trimestrielle (1+2,6 = ~ 4 x par an). La résolution de la « variation saisonnière » à 45 jours (soit 365/45=~8 fois par an), est en pratique de 365 / (365/45 ± 2,6) = 34 à 66 jours De la même manière, la « variation hebdomadaire » à 5 jours verra une résolution sur le Welch de 5 jours ± 4h ⇒ effet négligeable A fortiori, pour la fréquence quotidienne, résolution de 1 jour ± 10 minutes ⇒ effet insignifiant. Quelques remarques supplémentaires : ajouter des 0 à la fin du signal centré augmente -un peu- le confort visuel de façon artificielle. Plus précisemment, cela permet d’augmenter le nombre de points, sans affecter la courbe qui supporte les points. En particulier, cela n’améliore pas la résolution de la courbe. Passer la fft sur la périodisée au lieu d’une unique période facilite -un peu- l'identification de l'extrémité de certains lobes. Mais hormis en basse fréquence, cette info est déjà dans le Welch (1) Sur le principe, on aurait pu « corriger le climat » par la même méthode. Mais cela revient en pratique à forcer la température « redressée » vers 11,5°C en janvier et décembre. Pertinent pour étudier des variations quotidiennes, limite pour du mensuel. Le reste, on oublie ! Si on disposait d’un relevé sur une dizaine d’années, là, ça deviendrait pertinent -jusqu’aux variations annuelles, pas plus !

Revenons aux courbes obtenues.

// analyse spectrale 1 : FFT toute simple
F=fft(temp_centree);
mag=abs(F)/N;
phi=phasemag(F);

// analyse spectrale 2 : Welch sur périodisée
periodisee = repmat(temp_centree,50,1);
K=2^int(log2(N));
welch=pspect(K/2,K,'hm',temp_centree)/N;

A titre indicatif, voici les lignes de code scilab utilisées pour tracer et mettre en forme les courbes.
 

xlim = round(5.25*365);
xklim = round(xlim/N*K);
xk = (0:xklim-1)';
figure(2);

subplot(3,1,1); plot(x(1:xlim),mag(1:xlim));
xgrid;
xtitle('Spectre de la température `corrigée du climat`','','Magnitude de la FFT (°C)');
a1=get("current_axes");
a1.data_bounds=[0;xlim;0.01; 10;0;0];
a1.tight_limits="on";
a1.log_flags="nln"; 

subplot(3,1,2); plot(xk/K*N, welch(1:xklim));
xgrid; xtitle('Estimateur de la densité spectrale de puissance','Fréquence (x fois par an)','Welch (°C²)');
a2=get("current_axes");
a2.data_bounds=[0;xlim;0.0001;1;0;0];
a2.tight_limits="on";
a2.log_flags="nln";

subplot(3,1,3); plot(x(1:xlim),phi(1:xlim),'r');
xgrid; xtitle('','','Phase de la FFT (°)');
a3=get("current_axes");
a3.data_bounds=[0;xlim;-180;180;0;0];
a3.tight_limits="on";

La courbe de phase apparaît aléatoire. On n’en tirera aucun renseignement, ni sur la phase φ, ni sur le retard de groupe t_g= - dφ/dω.

Côté magnitude, hormis les raies annuelles et journalières attendues, on constate qu’il y a des composantes non négligeables entre les deux, sans qu’il ressorte de façon évidente une raie ou même un lobe. Après la partie journalière, seules les harmoniques de la composante journalière sont significatives. Au-dessus de l’harmonique journalière 3, le signal est intégralement sous le plancher du « bruit » annuel à journalier.

On fait un petit zoom jusqu’à la fréquence journalière.

Rien de flagrant. Le Welch rend juste le spectre un peu plus lisible.
Ah si, un truc à noter : la frontière entre variations "quotidienne" et "hebdomaire" est apparemment en fc=266 (= tous les ~1,37 jours)

Zoom supplémentaire à 3 jours soit F=365/3=~122 fois par an. Là, on attaque la variation "hebdomadaire"...

Dans le Welch, on voit apparaître deux petits lobes intéressants autour de 90 et 110.

Certains creux de la FFT (vers F= 69, 48, 41) ne sont pas clairement confirmés par le Welch. Dans cette zone, le Welch garde encore de la pertinence. Cela incite à la prudence sur ce qu’on peut faire dire à la FFT dans ces zones.

Pour finir, zoom sur les variations annuelle à semi-mensuelles. On abandonne le Welch au passage, ça évitera qu’il « ne bouche le port de Marseille » telle la sardine bien connue.

Ah, si, finalement on peut tirer quelque chose de la phase. Elle suggère que :

• soit la variation annuelle s’exprime avec la fondamentale F=1 jusqu’aux harmoniques F=4, suivie d’un lobe saisonnier entre F=5 et 8

• soit la variation annuelle se construit de la fondamentale F=1 jusqu’à l’harmonique F=7 ou 8

Prochain épisode : découpage du spectre en 4 bandes