Abstract

Les indices colorimétriques en télédétection ont un large domaine d’application pour caractériser différents types de surface. Bien souvent, les formes d’indices de télédétection sont définies empiriquement, que ce soit en sélection des longueurs d’onde, de la forme de l’équation, ou des coefficients dans l’équation permettant de calculer le dit indice. Et un très petit nombre d’études optimisent les paramètres utilisés pour ces indices. Ces indices sont alors utilisés tel-quel pour de la classification (souvent linaire comme les RandomForest). Mais aucune étude ne semble chercher la meilleure forme de façon automatique afin d’optimiser la classification et/ou la segmentation.

Notons qu’il est important d’optimiser l’indice en amont, car la transformation engendre une perte d’information et de caractéristique essentielle à la classification.
Un “mauvais indice” ne permet donc pas une bonne classification, pour pallier à cet effet les études privilégie presque toujours l’utilisation de plusieurs indices et de méthode de classification avancée

Le site www.indexdatabase.de référence ces différents indices (choix des longueurs d’ondes et des coefficients) en fonction du matériel d’acquisition. Souvent, le dit matériel est satellitaire (Landsat-5) donc moins influencé par les différences de luminosité, ce qui pose des soucis dans d’autres applications telles qu’en proxy-détection ce qui est notre cas.

L’approche standard pour sélectionner l’indice qui correspond le mieux a ce que l’on cherche a caractériser consiste a tester tous les indices disponibles en fonction de nos longueurs d’ondes avec une simple corrélation entre les indices et une vérité terrain Svitlana Kokh. Ce n’est pas une mauvaise approche, mais elle reste sous-optimal et fastidieuse, car elle nécessite de coder l’ensemble des indices et leurs différentes versions et le niveau de corrélation n’est pas le meilleur estimateur.

Ainsi, l’objectif de cette étude est de développer une méthode de recherche de l’indice optimal grâce à une approche statistique par descente de gradient, sur différentes formes d’équation générique. Nous allons également tester de nouvelles approches par traitement du signal et analyse d’image. L’apprentissage est effectué à travers des mécanismes de deep-learning Tensorflow 2.0.

À titre d’illustration, l’indice le plus utilisé pour la végétation est NDVI qui offre une mean_iou de $72.7\%$, tandis que les “meilleurs” indices empiriques offre une mean_iou au mieux de $82.3\%$, la ou nos des formes d’indice généré automatiquement propose une mean_iou de $87\%$. Cette différence est suffisamment significative pour améliorer la segmentation et la robustesse de l’indice a différents facteurs extérieurs, ainsi que la forme des éléments détectée.

Materiel

Les images ont été acquises grâce à la caméra multispectrale à six bandes Airphen. C’est une caméra scientifique multispectrale développée par des agronomes pour des applications agricoles. Elle peut être intégrée dans différents types de plates-formes telles que des drones, des robots de phénotypage, etc.

La caméra a été configurée en utilisant les bandes 450/570/675/710/730/850 nm avec un FWHM de 10 nm. La focal de chaque objectif est de 8 mm. Leurs résolutions brutes, pour chaque bande spectrale sont de 1280x960 px avec une précision de 12 bits. Enfin, l’appareil est équipé d’une antenne GPS interne, qui peut être utilisée pour obtenir la distance par rapport au sol.

L’alignement est affiné en deux étapes, avec (i) une estimation approximative de l’enregistrement affine et (ii) un enregistrement en perspective pour le raffinement et la précision grâce à la détection et la mise en correspondance de points clé. La méthode utilisée pour l’enregistrement est basée sur des travaux antérieurs Two-step Multi-spectral Registration, où la méthode montre une précision de l’enregistrement jusqu’au sous-pixel.

Données

Les données ont été acquises sur le site de l’IRSTEA de Montoldre, dans l’Allier (03), en France dans le cadre de l’ANR Challenge Rose.

12 images présentant des caractéristiques très distinctes en terme d’illumination (ombre, matin, soir, plein soleil, nuageux, …) ont été acquises.

Ainsi 3 vérités terrain ont été définies :

Attention : cet article correspond à de premiers travaux sur le sujet en utilisant des images prises en bonnes conditions.
Des résultats complémentaires sont disponibles dans la publication disponible ici.

Exemple pour la végétation :

Image originale	Masque

@Article{rs13122261,
AUTHOR = {Vayssade, Jehan-Antoine and Paoli, Jean-Noël and Gée, Christelle and Jones, Gawain},
TITLE = {DeepIndices: Remote Sensing Indices Based on Approximation of Functions through Deep-Learning, Application to Uncalibrated Vegetation Images},
JOURNAL = {Remote Sensing},
VOLUME = {13},
YEAR = {2021},
NUMBER = {12},
ARTICLE-NUMBER = {2261},
URL = {https://www.mdpi.com/2072-4292/13/12/2261},
ISSN = {2072-4292},
DOI = {10.3390/rs13122261}
}

Pré-traitement de l’information

Les bandes spectrales présentent par nature un bruit important associé au capteur CCD, ce qui pose des soucis lors de la normalisation.
Pour pallier à cet effet, 1% du signal minimum et maximum est supprimé par le calcul des quantiles, puis chaque bande est normalisée dans l’intervalle $[0,1]$.

Nous ajoutons également quelques transformations des bandes spectrales afin d’enrichir le pool d’information et de prendre en compte les gradients dans l’image.
Le choix s’est orienté vers 5 informations importantes :

1. L’image de la déviation standard entre toutes les bandes $\lambda_{dev}$
2. Le gradient via un filtre de Sharr sur $\lambda_{dev}$ noté $\lambda_{grad}$
3. Les valeurs propres maximales de la matrice Hessienne de $\lambda_{dev}$ notées $\lambda_{ridge_{max}}$
4. Les valeurs propres minimales de la matrice Hessienne de $\lambda_{dev}$ notées $\lambda_{ridge_{min}}$
5. Le laplacien de $\lambda_{dev}$ noté $\lambda_{laplace}$

Déviation standard	Filtre de Sharr	Valeurs propres max	Valeurs propres min	Laplacien

En effet, ces 5 transformations offrent d’une part une information importante sur le mélange spectral grâce à la déviation standard. Le filtre de Sharr et les valeurs propres maximum donne une information spatiale importante sur la rupture des gradients, donc sur la limite extérieure des objets, ce qui aidera probablement a la convergence. Les valeurs propres minimum, aussi appeler ridge, quant à elle semble détecté facilement les éléments fins tels que les monocotylédones pour des images de végétations.

Indices et forme d’équation

Les indices colorimétriques ont prouvé leur efficacité dans la description des surfaces, ainsi tous les ans différents articles de revue concernant les indices colorimétriques sortent, souvent dans le cadre d’un premier article de thèse appliqué à un domaine d’étude spécifique (Xue Jinru, Çağatay Tanrıverdi, Jiri Mezera). Nous allons repartir de la base grâce au site www.indexdatabase.de, à partir de cette base de données, 89 indices de végétation ont été identifiés comme compatibles avec nos longueurs d’ondes, ils seront ainsi testés. D’autre part, ils seront comparés à nos indices conçus automatiquement. Six formes d’équations simples ont été extraites à partir de l’ensemble des 519 indices de la base de données :

En analysant ces différentes équations, nous pouvons définir de nouvelles formes d’équation générique qui les synthétisent, et qui prennent en compte l’ensemble des bandes spectrales ainsi qu’un certain voisinage. Deux autres formes d’équations sont également intéressantes à optimiser, qui sont respectivement les approximations de fonction continue par développement de Taylor, ainsi que les approximations de fonction continue par morceaux grâce aux opérateurs morphologiques. Ces équations seront alors optimisées permettant de définir de nouveaux indices automatiquement, dans la seconde partie. L’ensemble de ces méthodes est développé via TensorFlow et les paragraphes suivants présentent les différents modèles :

Linéaire : Pour synthétiser les équations 1 et 2, nous pouvons définir une simple équation linéaire telle que $I = \sum_{i=0}^{N}{\alpha_i \lambda_i}$. Cette équation peut se généraliser au domaine 2D par une Convolution avec dans sa forme simple une taille de noyaux $k=1,1$.

Différence Linéaire : Pour synthétiser les équations 3, 4 et 5, un modèle simple reposant sur une division de deux modèles précédents est possible. De la même façon, cette forme est généralisable au domaine 2D et correspond alors à deux Convolution2D, l’une pour le numérateur, l’autre pour le dénominateur. Nous fixons dans un premier temps la taille des noyaux à $k=1,1$ pour ne pas prendre en considération le voisinage. Nous avons utilisé l’opérateur “nan div” pour remplacer les résultats “not a number” par zéro.

Polynomial : D’après le théorème de Stone–Weierstrass, toute fonction continue définie sur un segment peut être approchée uniformément par des fonctions polynomiales. Puisque l’on restreint le segment au domaine [0-1], les polynômes de Bernstein en sont une démonstration commune. Ainsi, toutes formes d’indices colorimétriques peut être approximée par un polynôme $I = \sum_{i=0}^{N}{\alpha_i \lambda_i}^{\delta_i}$. Pour des raisons d’implémentations, deux opérateurs sont définis. D’une part l’exposant $\lambda_i^\delta$, puis l’équation linéaire définie par une Convolution2D pour garder la possibilité de l’étendre au domaine 2D.

Différence Polynomial : En suivant l’idée de la forme précédente et pour respecter la forme d’équations 6, nous introduisons simplement une différence de polynôme, l’un au numérateur permettant d’optimiser la classe recherchée et l’autre au dénominateur jouant le rôle de normalisateur $\frac{\sum_{i=0}^{N} \alpha_i \lambda_i ^ {\delta_i} + A} {\sum_{j=0}^{N} \alpha_j \lambda_j ^ {\delta_j} + B}$. (voir les réultats si pertinant -> puisque thérème de StoneWeierstrass)

Approximation Universelle de Fonction : En utilisant les développements de Taylor, nous pouvons décomposer n’importe quelle fonction $f(x)$ en $f(x) = f(0) + f’(x) x + \frac{1}{2} f’‘(x) x^2 + \frac{1}{6} f’‘’(x) x^3 + o(x^3)$. Une approche pour apprendre cette forme de développement est proposée par Huang et al., 2017 que l’on appelle communément DenseNet et correspond alors à la somme de la concaténation du signal et de ses dérivées $\mathbf{x} \to \left[\mathbf{x}, f_1(\mathbf{x}), f_2(\mathbf{x}, f_1(\mathbf{x})), f_3(\mathbf{x}, f_1(\mathbf{x}), f_2(\mathbf{x}, f_1(\mathbf{x})), \ldots\right]$. Les dérivées sont disponibles par l’utilisation d’une Convolution2D.

https://d2l.ai/chapter_convolutional-modern/densenet.html

Filtre d’entrée : Pour supprimer une partie du signal qui ne serait pas indispensable, nous étudierons l’ajout d’un filtre passe-bande, en amont du réseau. Effectivement, pourquoi prendre en considération l’ensemble des valeurs d’un signal ? Il est très probable que seulement une partie de ce dernier caractérise ce que l’on cherche. Un bon exemple concerne les indices de végétation, seules les valeurs fortes dans le vert et le proche infra-rouge, ainsi que les valeurs faibles dans le rouge et le bleu sont caractéristiques de la végétation :

Note : C’est d’ailleurs le principe de l’indice NDVI, en raison des couches spongieuses qui se trouvent sur leur face arrière, les feuilles réfléchissent beaucoup de lumière dans le proche infrarouge, ce qui contraste fortement avec la plupart des objets non-végétaux. Lorsque la plante est déshydratée ou stressée, la couche spongieuse se compresse et les feuilles réfléchissent moins de lumière dans le proche infrarouge, rejoignant les valeurs du rouge, dans le domaine visible. Ainsi, la combinaison mathématique de ces deux signaux peut aider à différencier les plantes des objets non-végétaux et les plantes saines des plantes malades. Cependant, cet indice est alors moins intéressant lorsqu’il s’agit de détecter seulement la végétation et est fortement influencé par l’ombrage ou la chaleur.

Nous allons donc ajouter un filtre dans les équations précédentes pour supprimer les valeurs “indésirables” par l’utilisation de deux seuils a et b, qui seront également appris. S’il s’avère que l’ensemble du signal est intéressant, ces deux paramètres ne changeront pas et les valeurs seront a=0 et b=1, dans le cas contraire, a et b changeront. Pour appliquer le premier seuil, on utilise l’équation $y = \max(x-a,0) \div (1-a)$ et permet donc de supprimer les valeurs basses. À partir de ce nouveau signal, on applique l’équation $z = \max(b-y,0) \div b$ pour supprimer les valeurs hautes.

Ces deux paramètres permettent aussi de gérer des paramètres extérieurs, en effet lorsque a devient négatif, le paramètre a pour effet de rehausser le signal, ce dernier est alors compris dans l’intervalle $[a,1]$. Le second paramètre b permet par exemple de supprimer les effets de spéculaires ou de bruit trop important.

Approximation de fonction par érosion-dilatation : En étendant l’idée de sélection de bandes, nous pouvons non pas en sélectionner une, mais plusieurs, grâce à des opérateurs morphologiques notés $x \oplus s = \max_k(x_k + s_k)$ et $x \ominus s = \max_k(s_k - x_k)$ avec $x_k$ les différentes bandes spectrales et les $s_k$ le coefficient d’érosion de la bande spectrale $\ominus$ ou de dilatation $\oplus$. En augmentant le nombre d’érosion-dilatation, c’est-à-dire, en définissant plusieurs $s_k$ noté $s_{k,i}$, nous pouvons approximer n’importe quelle fonction continue par morceaux (Ranjan Mondal). On pose alors $z_i^{+} = x \oplus s_i = \max_k(x_k + s_{k,i})$ qui est le neurone structurant de dilatation $i$ et $z_i^{-} = x \ominus s_i = \max_k(s_{k,i} - x_k)$ qui est le neurone structurant d’érosion. Pour obtenir la sortie $I = \sum_{i=0}^{N}{z_i^{+}w_i^{+}} + \sum_{i=0}^{N}{z_i^{-}w_i^{-}}$ dont les $w_i^{+}$ et les $w_i^{-}$ sont les coefficients de combinaison linéaire obtenue par Convolution 2D pour rester dans le domaine 2D.

Raffinement spatial de l’indice : Pour prendre en compte différentes échelles dans l’image, nous étudions l’ajout d’une partie en aval du réseau. Appelé “Spatial Pyramide Refinement Block” et consiste en la somme de différentes Convolution2D dont les tailles de noyaux ont été fixées à 1,3,5,7. Le résultat de chaque convolution est concaténé, et le résultat final est donné par une combinaison linéaire de ces différentes échelles modélisées encore une fois par une Convolution2D avec k=1.

Fonction de perte

Dans cette partie, nous présenterons les différentes fonctions de pertes possibles en relation à la classification binaire. Cette partie est inspirée de https://lars76.github.io/neural-networks/object-detection/losses-for-segmentation/. Pour obtenir un indice et faciliter la convergence, on s’intéressera uniquement aux valeurs comprises entre 0 et 1 en sortie de la dernière couche, via une fonction d’activation de type ClippedReLU ($0<=x<=1$), ce qui est négatif ou nul sera donc indésirable et supérieur ou égal à 1 la classe recherchée et entre les deux la frontière d’indécision à optimiser. Les valeurs correspondent alors à la probabilité que le pixel soit la surface recherchée $P(Y=1)=1=p$ ou non $P(Y=0) = 1-p$. Dans ce cas, 4 fonctions sont communément utilisées : Cross-Entropy, Weigthed-Cross-Entropy, Balanced-Cross-Entropy et Focal-Loss (dans la suite nous noterons $p$ la vérité terrain et $\hat{p}$ la prédiction).

Cross-Entropy : La cross-entropy peut être définie par $CE(p,\hat{p}) = -p\log(\hat{p}) -(1-p)\log(1-\hat{p})$. C’est une fonction de perte standard qui ne prend pas en compte le ratio des classes (unbalanced). Elle permet de mesurer la séparabilité des classes. Lorsque $CE \rightarrow 0$ il n’y a plus d’entropie et $\hat{p} \rightarrow p$, dans le cas contraire, si $CE \rightarrow 1$ alors $\hat{p}$ et $p$ sont indissociables. La fonction de perte $CE$ ne semble pas être adaptée à notre cas, puisque l’on utilise une vérité terrain dont le ratio entre surface recherchée et le reste, des différences pour chaque image.

Weigthed-Cross-Entropy : C’est une variante de $CE$ où tous les membres de classes positive ($p=1$) sont pondérés par un coefficient $\beta$, lequel permet de réguler le ratio des classes. On écrit alors $WEC(p,\hat{p}) = -\beta p\log(\hat{p}) -(1-p)\log(1-\hat{p})$. Pour décrémenter l’effet de la classe $Y=0$ on fixe $\beta>1$ et $\beta<1$ dans le cas contraire. Le coefficient $\beta$ peut être fixé pour maximiser une classe spécifique, ou contrebalancer le ratio des classes. Le coefficient $\beta$ pourrait néanmoins être calculé indépendamment pour chaque image.

Balanced-Cross-Entropy : Cette variante est similaire à WCE, on introduit cependant $1-\beta$ pour la classe négative. On a donc l’équation $BCE(p,\hat{p}) = -\beta p\log(\hat{p}) -(1-\beta)(1-p)\log(1-\hat{p})$. Cette forme est intéressante, puisque $\beta$ peut soit être fixé comme pour $WCE$, soit être calculé pour chaque image, dans ce cas $\beta = \sum(\frac{P(Y=1)}{P(Y=0)})$.

Focal Loss : Cette fonction, proche de BCE essaie de décrémenter l’influence des éléments “faciles” pour améliorer en priorité les éléments difficilement séparables. On utilise alors deux coefficients $\alpha$ et $\gamma$ sur la distribution de $\hat{p}$. La fonction s’écrit alors $FL(p,\hat{p}) = -(\alpha (1-\hat{p})^\gamma) p\log(\hat{p}) -(1-\alpha)\hat{p}^\gamma(1-p)\log(1-\hat{p})$. La fonction de perte $FL$ ne permet pas de prendre en considération le ratio des classes différents pour chaque image. La fonction ne permet pas d’optimiser l’indice dans l’intervalle [0-1] avec efficacité et donc l’optimisation des métriques définis plus loin. Bien que le coefficient $\alpha$ puisse être calculé automatiquement, le coefficient $\gamma$ est difficile à fixer.

Trois autres fonctions de perte sont intéressantes à étudier quand il s’agit de segmentation : Dice-Loss, TverskyIndex-Loss et IntersectionOverUnion-Loss souvent utilisés dans leurs versions métrique.

Dice-Loss : La fonction de perte est définie par $\text{DL}(p, \hat{p}) = 1 - \frac{2p\hat{p}}{p + \hat{p}}$

IntersectionOverUnion-Loss : Récemment, Y.Wang et al ont proposé une solution pour optimiser une approximation de l’intersection sur l’union dans le cas de segmentation binaire. La fonction de perte est définie par $IoU = 1 - \frac{I(p, \hat{p})}{U(p, \hat{p})}$ avec $I(p, \hat{p}) = p\hat{p}$ et $U(p, \hat{p}) = p+\hat{p} - p\hat{p}$. Les performances de cette fonction de perte semblent plus efficaces que les méthodes “simples” précédemment citées (Gellért Máttyus, Dingfu Zhou, Zhaohui Zheng).

Le choix se fera donc sur la fonction de perte $BCE$ avec un calcul automatique de $\beta$.

Métrique d’évaluation

Lorsque le nombre d’éléments entre chaque classe est fortement déséquilibré, beaucoup de métriques standard sont inefficaces et donc inadaptées. Par exemple, notre vérité terrain sol/végétation comporte ~83% de sol et donc ~17% de végétation, une mauvaise métrique comme l’accuracy montrera de “bonnes” performances. Dans le cas de l’accuracy $\frac{tp+tn}{tp+tn+fp+fn}$ si l’évaluation ne détecte que du sol alors les performances seront de ~83% ce qui n’est donc pas représentatif. Ce que l’on cherche alors est de trouver des métriques permettant de prendre en compte le ratio de chaque classe. Communément, les performances des indices colorimétriques est calculée par une cross-entropy $-\frac{1}{N} \sum_{i=0}^{N}{y_{true}\log(y_{pred}) + (1-y_{true})\log(1-y_{pred})}$ entre l’indice et une vérité terrain. Comme nous l’avons vu précédemment, cette métrique n’est pas non plus adaptée car elle ne prend pas en compte le ratio des classes. Il existe beaucoup d’autres métriques (Abdel Aziz, David Martin Ward Powers, Takaya Saito, László A. Jeni) :

Balanced Accuracy : Aussi appelé Youden Statistics pour prendre en compte le ratio des classes, la forme de l’accuracy est transformée afin de calculer séparément les performances de chaque classe et devient donc $\frac{1}{2}(\frac{tp}{tp+fp} + \frac{tn}{tn+fn})$

Dice : Le coefficient de Dice, également appelé indice de recouvrement, est communément utilisé pour vérifier les performances des algorithmes de segmentation $\frac{2p\hat{p}}{p+\hat{p}}$

Mean Intersection Over Union : C’est une autre métrique souvent utilisée pour mesurer les performances de la segmentation, notée $\frac{p\hat{p}}{p+\hat{p} - p\hat{p}}$

Precision : Correspond à une partie de l’équation de balanced accuracy $\frac{tp}{tp+fp}$ et permet de répondre à la question : Quelle proportion des identifications positives a été effectivement correcte ?

Recall : Quelle proportion de vrais positifs a été identifiée correctement ? $\frac{tp}{tp+fn}$

Matthews correlation : $\frac{tptn - fpfn}{\sqrt((tp+fp)(tp+fn)(tn+fp)*(tn+fn))}$

Comparaison avec les indices standards

Pour effectuer une comparaison juste/équitable, il est nécessaire d’optimiser chaque indice standard. Pour ce faire, un réseau de neurones minimal est utilisé afin d’apprendre une droite de régression. Le réseau est donc composé de l’indice spectral, suivi d’une normalisation $x=(x-min)/(min-max)$, puis d’une convolution 2D avec une taille de noyau de $k=1,1$. Puisque l’indice génère une seule dimension, l’équation de sortie est alors $I = \alpha * NormalizedIndex + \beta$. Pour effectuer la classification de la même façon que notre méthode, une activation de type ClipedRelu est utilisée. Évidemment, les mêmes métriques et fonction de perte sont utilisées.

Réseau de neurones simple

Resultats sur les indices de végétations

Performance des indices généré (balanced-accuracy)

Performance des indices généré (dice)

Performance des indices généré (mean-iou) :

Performances des indices standars :

Conclusion

A partir de ces résultats nous pouvons conclure que n'importe quelle combinaison linéaire simple permet d’obtenir des indices de végétation tout aussi performants que les indices définis empiriquement, de plus des modèles plus avancés permettent d'obtenir des résultats plus performant et adaptés aux données récoltées.
Notons également que les indices NDVI = $\frac{n-r}{n+r}$, EnhancedVegetationIndex = $2.5*\frac{n-r}{n+6r-7.5b+1}$, EnhancedVegetationIndex2 = $2.4\frac{n-r}{n+r+1}$ et EnhancedVegetationIndex3 = $2.5\frac{n-r}{n+2.4r+1}$ sont très proche dans leurs formes et propose des performances très différentes après optimisation. On conviendra donc qu'il est plus important d'optimiser la forme de l'équation et les coefficients des bandes spectrales qu'une simple régression, c'est pourquoi nos indices générés automatiquement sont bien plus performants !