# Sur la modélisation du bruit en détection d'anomalie magnétique

Jean-Pierre Tanguy, Abdel Boudraa

<figure><img src="image_1.png" /></figure>

## To cite this version:

Jean-Pierre Tanguy, Abdel Boudraa. Sur la modélisation du bruit en détection d'anomalie magnétique. Colloque Francophone de Traitement du Signal et des Images, GRETSI, Sep 2017, Juan Les Pins, France. hal-04691099

HAL Id: hal-04691099

https://hal.science/hal-04691099v1

Submitted on 11 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Sur la modélisation du bruit en détection d'anomalie magnétique

JEAN-PIERRE TANGUY¹, ABDEL BOUDRAA²

¹DGA TN Brest

²IRENav (EA 3634) Ecole Navale/Arts-Métiers ParisTech, BCRM Brest, CC 600, 29240 Brest

(jean-pierre.tanguy@intradef.gouv.fr, boudra@ecole-navale.fr)

**Résumé:** En détection magnétique aéroportée, la signature de la cible recherchée est en général entachée d'un bruit, dit magnétique, résultant d'origines physiques diverses. Il est utile de disposer d'un modèle réaliste de ce bruit, afin d'estimer les performances des systèmes de détection. Un modèle de bruit obtenu par filtrage linéaire d'un bruit blanc a déjà été étudié. Ce modèle n'est pas entièrement satisfaisant, car l'écart entre la densité spectrale (dsp) du bruit de sortie et la DSP de référence basée sur des mesures physiques réelles augmente avec la fréquence. En outre, nous n'avons pas la garantie que les échantillons temporels obtenus soient représentatifs d'un véritable bruit magnétique. Le but ici est de passer en revue les moyens susceptibles d'aboutir à une modélisation plus satisfaisante et de mettre en évidence les difficultés rencontrées.

**Abstract:** In airborne magnetic detection, the signature of the searched target is usually corrupted by a noise referred to as magnetic noise, and which has various physical origins. It is useful to have a realistic model of this noise, in order to assess the performances of detection systems. A noise model obtained by linear filtering of a white noise has already been investigated. This model is not completely satisfactory, as the discrepancy between the power spectral density of the output noise and the reference standard noise power spectral density based on real measurements increases with frequency. Besides, there is no guarantee that the time samples obtained are representative of a true magnetic noise. It is intended here to make an overview of the means likely to result in a more satisfactory modeling and to point out the difficulties encountered.

## 1 Introduction

La détection d'objets ferromagnétiques cachés (sous-marins,...) présente un intérêt, tant civil que militaire. La présence de tels objets produit une perturbation du champ magnétique terrestre, désignée par anomalie magnétique. La détection de ces sources magnétiques cachées, a conduit à différentes méthodes, telle que la détection d'anomalie magnétique (MAD). Cette méthode est passive et présente l'avantage de la discrétion. Le détecteur MAD est en général un magnétomètre scalaire RMN. La signature reçue par le capteur est entachée d'un bruit, dit magnétique, qui résulte d'origines physiques diverses. Pour l'analyse des systèmes MAD, la cible est souvent modélisée par un dipôle magnétique et le bruit par un processus aléatoire (PA) stationnaire gaussien. L'hypothèse d'un dipôle magnétique est valable si la distance cible capteur est grande devant la dimension de la cible. Plus récemment, des modèles plus précis faisant intervenir des multi-pôle ont récemment été étudiés [1],[2]. La détection MAD repose sur le filtrage adapté. Des mesures magnétiques ont montré que le bruit, principalement d'origine géomagnétique, présente souvent une densité spectrale (dsp) en $1/ν^α$ avec $0 < α < 2$ [3]. Une autre DSP du bruit, peut-être plus réaliste a été introduite dans [4]. La recherche d'un modèle de bruit, s'appuyant sur cette dernière DSP, est présentée dans [5]. En considérant le bruit comme la sortie d'un filtre linéaire excité par un bruit blanc, les méthodes de synthèse des filtres numériques peuvent être utilisées pour générer un PA stationnaire de DSP donnée. Cependant, le modèle

obtenu n'étant pas entièrement satisfaisant, on se propose ici de l'améliorer, l'accent étant mis sur la prise en compte de la loi de probabilité du bruit. Une connaissance fine du bruit est nécessaire pour une modélisation réaliste qui contribuera à une évaluation théorique correcte des performances (en portée) d'un système MAD. Celles-ci dépendent bien sûr du type de traitement du signal utilisé, de la qualité de la modélisation du signal utile (anomalie) et du bruit (bruit magnétique). L'objectif est d'aboutir à un modèle réaliste permettant l'évaluation de systèmes de détection magnétique et leur optimisation en situation opérationnelle.

## 2 Données du problème

La donnée disponible est une DSP, $γ(ν)$, dit gabarit standard [4], obtenu à partir de mesures effectuées en un lieu donné. Il est défini pour deux cas "fort bruit" et "faible bruit" correspondant respectivement à des petits fonds et des grands fonds [5] (Figure 1).

<figure><img src="image_3.png" /><figcaption>**Figure 1:** Gabarit standard pour les cas bruit fort et faible bruit.</figcaption></figure>
La dsp est donnée par $\gamma(\nu) = A_1 + \gamma_2(\nu)$, avec

$$ \gamma_2(\nu) = \frac{B_1\nu^2 + C_1}{\nu^4} = \frac{B_1}{\nu^2} + \frac{C_1}{\nu^4}. $$

Un terme correctif $d^4$ a été ajouté au dénominateur pour que la dsp soit définie même à $\nu = 0$. La dsp modifiée s'écrit $\gamma_m(\nu) = \frac{B_1\nu^2 + C_1}{\nu^4 + d^4} \quad \forall \nu$, ne s'écarte pas significativement de la vraie l'expression. On peut vérifier que la TF inverse de $\gamma_m$, contrairement à celle de $\gamma_2$ (qui, à $\infty$, tend vers $\infty$, au lieu de tendre vers 0) représente une fonction d'autocorrélation (FAC) d'expression

$$ \Gamma_m(\tau) = \left(\frac{A}{\alpha}\nu_0 + \frac{B}{\alpha}\right)e^{-\alpha|\tau|}\cos(2\pi\nu_0\tau) - \frac{A}{2\pi}e^{-\alpha|\tau|}\sin(2\pi\nu_0|\tau|) \\ \text{ave } \alpha = \frac{B_1 - C_1}{2d\sqrt{2}} \quad \beta = \frac{C_1}{2d^2} \quad \nu_0 = \frac{d}{\sqrt{2}} \quad A = 4\pi^2\alpha \quad B = 4\pi^2\beta $$

On peut vérifier que $\Gamma_m(\tau)$ possède les trois propriétés d'une FAC.

### 3 Hypothèses

Pour la recherche d'un modèle de bruit, on suit l'approche exposée dans [6]. Bien que le bruit soit non stationnaire sur une longue période, il raisonnable de considérer que le bruit est localement stationnaire, sur la durée du signal utile, car ce dernier est présent seulement sur une très courte durée.

#### 3.1 Bruit blanc ou coloré

Le bruit n'est manifestement pas blanc, d'après son allure graphique (Fig. 1). L'hypothèse d'un bruit blanc (faite dans [6] pour le bruit géomagnétique) serait valable si la bande du signal utile était très petite devant celle du bruit et si la dsp du bruit avait des variations lentes sur la bande du signal utile. Ce n'est pas le cas ici, car la bande passante du signal utile (estimée entre 0.3 et 0.4 Hz) est très proche de celle du bruit (0.5 Hz).

#### 3.2 Densité de probabilité

Si la probabilité d'avoir une amplitude inférieure ou égale à $3\sigma$ est d'environ 99% [6], il y a des chances pour que le bruit soit gaussien. L'hypothèse d'un bruit laplacien est parfois faite [7]. Pour un bruit Laplacien centré d'écart-type $\sigma$, la densité de probabilité est donnée par $f(x) = \frac{1}{\sigma\sqrt{2}}e^{-\frac{\sqrt{2}}{\sigma}|x|}$. Ainsi, on peut aisément générer des échantillons $x_i$ distribués selon une loi laplacienne par les relations [6]

$$ x_i = \frac{\sigma}{\sqrt{2}} \ln(2u_i) \quad \text{si} \quad 0 < u_i < \frac{1}{2} $$

$$ x_i = -\frac{\sigma}{\sqrt{2}} \ln[2(1-u_i)] \quad \text{si} \quad \frac{1}{2} < u_i < 1 $$

où $u_i$ est une suite de nombres pseudo-aléatoires uniformément distribués entre 0 and 1. Une non linéarité judicieusement choisie permet donc d'obtenir une loi de probabilité donnée à partir d'une loi gaussienne. Cette non linéarité entraîne malheureusement une distorsion de la dsp du PA que l'on cherche à générer. Cette distorsion est cependant peu importante si la caractéristique de la non-linéarité est relativement lisse. Il est possible d'exprimer sous forme d'une somme de série [8] la dsp du PA en sortie de la non-linéarité, c.à.d. la dsp du PA ayant la loi de probabilité souhaitée. Un calcul numérique sera forcément approché, puisqu'il faut limiter le nombre de termes de la série. Cependant, une expression analytique de la dsp peut être obtenue si la caractéristique de non linéarité est une fonction erreur. On trouve dans la littérature des méthodes de génération de PAs ayant une loi de probabilité donnée [9]. Le bruit est modélisé par une superposition de sinusoïdes aléatoires pour lesquelles les amplitudes $A_i$, les phases $\varphi_i$ et éventuellement les fréquences $f_i$ sont des variables

aléatoires. $$ x(n) = \sum_{i=1}^{p} A_i \cos(2\pi f_i n + \varphi_i) $$

Dans une autre méthode sur la synthèse de bruits gaussiens et non gaussiens, le spectre d'amplitude est supposé connu (c'est le cas de notre application) [10]. Des manipulations du spectre de phase et les moments d'ordre supérieur (Skewness, Kurtosis) permettent de générer des bruits avec différentes lois de probabilité.

### 4 Méthodes de génération du bruit

Nous retenons les hypothèses suivantes pour obtenir un modèle approprié de bruit magnétique : le bruit est considéré stationnaire, au moins sur la durée du signal utile, il est centré et coloré dont les échantillons ont la même loi de probabilité. Deux approches sont considérées pour générer un bruit remplissant ces hypothèses [5].

#### 4.1 Approche temporelle

Dans cette approche, le bruit cherché résulte du filtrage linéaire d'un bruit blanc par un filtre de gain complexe

$$ H(z), \text{ donné par } H(z) = \frac{b_0 + b_1 z^{-1}}{1 + a_1 z^{-1} + a_2 z^{-2}} $$

Les coefficients $b_k$ et $a_k$ sont obtenus par des techniques classiques de synthèse des filtres numériques. La généralisation de cette expression est une fonction de transfert rationnelle de la forme

$$ H(z) = \frac{b_0 + b_1 z^{-1} + \dots + b_q z^{-q}}{1 + a_1 z^{-1} + \dots + a_p z^{-p}} $$

modèle pôle zéro appelé aussi processus autorégressif à moyenne adaptée (ARMA(p,q)). Le cas particulier où tous les coefficients $b_k$ sont nuls, sauf $b_0 = 1$,
correspond au modèle tout-pôle, ou modèle autorégressif d'ordre p (AR(p)), pour lequel la fonction de transfert n'a que des pôles, et s'exprime par

$$H(z) = \frac{1}{1 + a_1 z^{-1} + \dots + a_p z^{-p}} \quad . \text{p est l'ordre du}$$

modèle et les paramètres du modèle sont les coefficients $a_1, a_2, \dots, a_p$. L'avantage du modèle AR est son aptitude à modéliser n'importe quelle DSP non uniforme, si le nombre p de paramètres est suffisamment élevé. En outre, les processus AR sont plus faciles à manipuler que les processus ARMA. Des algorithmes bien connus (Levinson, Burg,...) permettent d'estimer les coefficients d'un modèle AR [11]. Une des difficultés de la modélisation autorégressive réside dans la détermination de l'ordre p. Il existe différents critères (mentionnons le critère bien connu de l'erreur finale de prédiction Final Prediction Error) et le critère d'information de Akaike), mais qui ne donnent pas nécessairement de résultats satisfaisants. Une règle pratique souvent employée consiste à choisir l'ordre p

tel que $\frac{N}{3} \prec p \prec \frac{N}{2}$ si N est le nombre d'échantillons temporels disponibles. On peut ainsi générer des échantillons de bruit et effectuer un filtrage linéaire de ces échantillons selon l'équation récursive

$$y(n) = x(n) - a_1 y(n-1) - a_2 y(n-2) - \dots - a_p y(n-p)$$

Si $x(n)$ est gaussien, il en est de même de $y(n)$, qui résulte d'une transformation linéaire de $x(n)$. Si $x(n)$ n'est pas gaussien, des transformations non linéaires doivent être effectuées pour obtenir la loi de probabilité souhaitée [12]. Ceci est illustré par la figure 2 avec l'exemple d'un filtre AR pour obtenir un bruit coloré. Bien sûr, n'importe quel filtre de mise en forme peut être employé. $g_1$ et $g_2$ désignent respectivement les transformations utilisées pour obtenir une variable aléatoire uniformément distribuée à partir d'une variable aléatoire gaussienne et pour obtenir une variable aléatoire laplacienne à partir d'une variable aléatoire uniformément distribuée. Une expression explicite de $g_2$ est donnée en III-C.

$y(n)$ est un bruit gaussien centré de fonction de répartition $F_Y(y)$. La variable aléatoire $u = F_Y(y)$ suit une loi uniforme sur l'intervalle $[0,1]$. (c'est vrai pour n'importe quelle fonction de répartition $F_Y(y)$) [12]. D'après un résultat bien connu, la fonction de répartition d'une variable aléatoire gaussienne Y est

$$F_Y(y) = \frac{1}{2} + \frac{1}{2} \operatorname{erf} \left( \frac{y}{\sigma_Y \sqrt{2}} \right) \quad \text{if } \sigma_Y \text{ est l'écart-type de Y.}$$

La transformation $g_1$ est $g_1(y) = F_Y(y)$

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Génération d'un bruit Laplacien – approche temporelle</figcaption></figure>

## 4.2 Approche fréquentielle

Cette approche ne fait pas intervenir le filtrage linéaire d'un bruit blanc. La DSP de référence (ici la DSP modifiée cf II-A) est considérée comme la DSP estimée par la méthode du périodogramme simple. On peut ainsi l'écrire $\hat{\gamma}(n) = \frac{1}{N} |X(n)|^2$ où $X(n)$ désigne la transformée de Fourier discrete (TFD) des échantillons

$$x(k) \quad X(n) = \sum_{k=0}^{N-1} x(k) e^{-j \frac{2\pi ik}{N}} \quad n = 0, 1, \dots, N-1$$

On a seulement accès au module de la TFD, et non pas à la TFD elle-même, qui est définie avec une phase inconnue $\Psi(n)$.

$$X(n) = |X(n)| e^{j\Psi(n)} = |X(n)| \beta(n) \quad \text{avec } \beta(n) = e^{j\Psi(n)}$$

On peut générer une suite de nombres complexes aléatoires $\beta(n)$ tels que $|\beta(n)|=1$, puis définir la suite $X(n)$ dans le domaine fréquentiel par $X(n) = \sqrt{N\hat{\gamma}(n)}\beta(n)$. Les échantillons temporels cherchés peuvent théoriquement être obtenus par transformée de Fourier discrète inverse (TFDI) de $X(n)$

$$x(k) = \frac{1}{N} \sum_{n=0}^{N-1} X(n) e^{j \frac{2\pi ik}{N}} \quad k = 0, 1, \dots, N-1$$

On note que les échantillons $x(k)$ sont en général à valeur complexe. Pour surmonter cette difficulté, on peut réarranger les échantillons fréquentiels $X(n)$ de façon à faire apparaître la symétrie hermitienne, conduisant ainsi à des échantillons temporels à valeur réelle. Une manière alternative est de considérer la DSP comme le spectre d'un signal analytique. La partie réelle de la TFDI de $X(n)$ est ainsi le processus temporel cherché. Jusqu'ici, nous ne sommes pas sûrs que les échantillons temporels obtenus de cette façon soient représentatifs du bruit magnétique. Il peut être approprié d'estimer leur loi de probabilité, et, si nécessaire, d'effectuer une transformation non linéaire pour obtenir des échantillons laplaciens.

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Génération d'un bruit Laplacien – approche fréquentielle</figcaption></figure>
# 5 Application à un bruit réel

On dispose de deux enregistrements réels de bruit, effectués à une fréquence $F_e=5$Hz lors de vols MAD. Pour ces deux enregistrements, on donne les moments d'ordres 1 à 4 ainsi que les DSP (périodogrammes simple et moyenné). Le Kurtosis est voisin de 3, ce qui est plutôt en faveur d'une loi gaussienne. Une hypothèse gaussienne simplifie considérablement la génération d'échantillons temporels. Pour l'approche temporelle, un modèle AR convient, la difficulté résidant alors dans le choix de son ordre. Pour l'approche fréquentielle, la procédure présentée en 4.2 s'applique sans difficulté.

<figure><img src="image_6.png" /><figcaption>Enregistrement de bruit de 375 échantillons</figcaption></figure>

Moyenne= 0.0107 ; Variance= 0.0230
Skewness= 0.2073 ; Kurtosis= 3.0909

<figure><img src="image_7.png" /><figcaption>Enregistrement de bruit de 3394 échantillons</figcaption></figure>

Moyenne= 6.09 e-4 ; Variance= 0.0074
Skewness= -0.0510 ; Kurtosis= 3.3612

# 6 Conclusion

Un travail antérieur sur la génération d'échantillons temporels d'un bruit magnétique a mis en évidence quelques difficultés sur un sujet, qui, à première vue, paraît classique. Pour surmonter ces difficultés, deux approches sont suggérées. Dans la première, le filtrage autorégressif d'un bruit blanc permet de générer un processus temporel de DSP donnée. Dans la deuxième, la

DSP voulue est considérée comme le périodogramme simple des données cherchées. Le processus temporel est généré par transformée de Fourier inverse. Pour ces deux approches, des hypothèses à la fois simples et raisonnables (réalistes ?) sont faites pour permettre une mise en œuvre aisée du processus de génération des échantillons temporels. Ces hypothèses sont liées essentiellement à la stationnarité, à la nature du bruit (blanc ou coloré) et à la loi de probabilité. La prochaine étape consistera à illustrer ces deux approches par des exemples de résultats.

# REFERENCES

[1] P. Pepe, S. Zozor, L.L. Rouve, J.L. Coulomb, Ch. Serviere et J. Muley, "Etude des performances de la détection d'anomalie magnétique généralisée au cas des cibles multipolaires," *Colloque GRETSI*, pp. 1-4, 2015.

[2] P. Pepe, S. Zozor, L.L. Rouve, J.L. Coulomb, Ch. Serviere and J. Muley, "Generalization of GLRT-based magnetic anomaly detection, EUSIPCO, pp. 1975-1979, 2015.

[3] A. Sheinker, A. Shkalim, N. Salomonski, B. Ginzburg, L. Frumkin and B.Z. Kaplan, "Processing of a scalar magnetometer signal contaminated by 1/f noise", *Sensors and Actuators A*, vol. 138, pp. 105-111, 2007.

[4] R. Blanpain, Traitement en temps réel du signal issu d'une sonde magnétométrique pour la détection d'anomalies magnétiques. Thèse INPG, 1979.

[5] J.P. Tanguy, "Modélisation du bruit magnétique," *Colloque GRETSI*, pp. 1-4, 2015.

[6] S.M. Kay, Fundamentals of statistical signal processing. Practical algorithm development. Prentice Hall, 2013.

[7] J. Schweiger, Evaluation of geomagnetic activity in the MAD frequency band (0.04 to 0.06 Hz), Masters thesis, Naval postgraduate school, 1982.

[8] U. Gujar and R. Kavanagh, "Generation of random signals with specified probability density functions and power density spectra," *IEEE Trans. Automatic Control*, vol. 13, no. 6, pp. 716-719, 1968.

[9] W.E. Thomson, "The response of a non-linear system to random noise," *Proc. IEE*, vol. 102, no. 1, pp. 46-48, 1955.

[9] S.M. Kay, "Representation and generation of non-gaussian wide-sense stationary random processes with arbitrary PSDs and a class of PDFs," *IEEE Trans. Sig. Proc.*, vol. 58, no. 7, pp. 3448-3458, 2010.

[10] K.D. Hsueh, R.P. Hamernik, "A generalized approach to random noise synthesis: Theory and computer simulation," *JASA*, vol. 87, no. 3, pp. 1207-1217, 1990.

[11] S.M. Kay, Modern spectral estimation. Theory and application, Prentice Hall 1988.