Modélisation GPR

Georadar (dispositif radio-technique de sondage souterrain, GPR, Ground Penetrating Radar), qui est actuellement très utilisé - de la cartographie des trous de lapin et de l' étude des lézards à la recherche de mines , reste un plaisir assez cher.


Affichage GPR (cadre de l'émission de télévision britannique "Command of the Time")

Mais il est possible d'évaluer ses capacités et d'étudier divers aspects de l'interaction du champ électromagnétique du géoradar avec l'environnement sans acquérir / louer un appareil «en fer». Le package gprMax ( gpr - de l'abréviation GPR, Max - les premières lettres du nom de James Clerk Maxwell, qui a jeté les bases de l'électrodynamique), distribué sous la licence GNU GPL v3, nous aidera.
Les auteurs de ce projet, commencé en 1996, sont Craig Warren de l'Université de Northumbria et Antonis Giannopoulos de l'Université d'Edimbourg. Le package a été initialement développé en C, puis complètement réécrit en combinaisons Python 3 / Cython.


L'installation du package nécessite des compilateurs installés qui prennent en charge OpenMP (Microsoft Visual C ++ 2015 Build Tools (cette version est recommandée!) Pour Windows / gcc pour Linux), la bibliothèque NumPy et le compilateur Cython. Après avoir téléchargé à partir du référentiel sur GitHub et décompressé le code source du projet, accédez au dossier racine et exécutez les commandes:

python setup.py build python setup.py install 

En guise de «démarrage rapide», nous envisageons de travailler avec le boîtier à l'aide d'un simple exemple bidimensionnel - l'antenne d'émission T d'un radar pulsé (impulsion GPR) émet une impulsion électromagnétique, dont une partie de l'énergie atteint directement l'antenne de réception R sous la forme d'une onde directe (DW - onde directe), et certaines pénètrent à travers le sable, il est réfléchi depuis la surface du cylindre conducteur et atteint l'antenne réceptrice sous la forme d'une onde réfléchie (RW - onde réfléchie):


Format de fichier d'entrée
Créez le dossier des modèles dans le dossier racine du projet, dans lequel nous plaçons le fichier texte hello.in contenant les commandes pour effectuer la simulation (les commandes suivantes correspondent à la (troisième) version actuelle du projet).

Chaque équipe a la forme:

 #: _1 _2 _3 ... 

Une seule commande peut être écrite sur une seule ligne et le premier caractère de la ligne contenant la commande doit être #.

Les commandes peuvent être accompagnées de commentaires:

 ## 

L'ordre des commandes est important pour les commandes de construction d'objets - ces commandes sont exécutées dans l'ordre où elles apparaissent dans le fichier d'entrée.

Forme d'impulsion

Une impulsion électromagnétique émise par un géoradar dure quelques fractions de nanoseconde, de plus, trois formes d'impulsions sont le plus souvent utilisées:

1. une période d'onde sinusoïdale (sinusoïdale)
2. Moment gaussien (gaussien)
3. Le chapeau mexicain (ricker) est proportionnel à la dérivée seconde de la fonction gaussienne, la forme de la courbe d'une telle impulsion ressemble à un sombrero (cette forme d'impulsion a été utilisée par le géophysicien américain Norman Ricker en 1953 pour étudier les signaux sismiques)
   

Pour notre exemple, nous sélectionnons une impulsion gaussienne (type d'impulsion - gaussienne) avec une fréquence centrale $$$f_c = 1 \\ GHz = 1 \ cdot 10 ^ 9 Hz$ défini par la commande:

 #waveform: gaussian 1 1e9 pulse 

(1 - amplitude d'impulsion conditionnelle, impulsion - étiquette d'impulsion)

Dans ce cas, l'élan utilisé dans la simulation est décrit par l'expression:

$$

$$W (t) = e ^ {- 2 \ cdot {\ pi} ^ 2 \ cdot {f_c} ^ 2 {(t- {1 \ over {f_c}}) ^ 2}$$

Modèle d'environnement et système de coordonnées

Dans la modélisation 2D, la zone étudiée est divisée en cellules d'une taille donnée, et le système de coordonnées du modèle ressemble à ceci - les axes X et Y forment le plan calculé (avec une largeur $$$w$ et grand $$$h$ ), la longueur du modèle le long de l'axe Z a une valeur égale à l'étape d'échantillonnage $$$\ Delta$ .


Lorsque vous choisissez une étape d'échantillonnage, vous pouvez suivre la règle de base - la taille de l'étape ne doit pas dépasser un dixième de la plus petite longueur d'onde étudiée dans le modèle («10 cellules par longueur d'onde»):

$$

$$\ Delta \ le 0,1 \ cdot {\ lambda} _ {min}$$

Pour déterminer la longueur d'onde, il faut connaître la fréquence maximale prise en compte dans le spectre du signal émis et la vitesse des ondes dans le milieu considéré.

Vitesse de propagation d'une onde électromagnétique dans un milieu à constante diélectrique relative $$${\ epsilon} _r$ , exprimé en centimètres par nanoseconde - l'unité de vitesse adoptée dans le radar est déterminée par l'expression:

$$

$$v = {30 \ over {\ sqrt {{\ epsilon} _r}}}$$

La longueur d'onde en centimètres est déterminée par l'expression:

$$

$$\ lambda = {v \ over f}$$

( $$$f$ - fréquence en GHz).
Pour afficher la forme de l'impulsion gaussienne et son spectre, vous pouvez utiliser la commande:

 python -m tools.plot_source_wave gaussian 1 1e9 5e-9 1e-12 -fft 

(gaussien - type d'impulsion, 1 - amplitude d'impulsion, 1e9 - fréquence centrale (1 GHz), 5e-9 - durée d'affichage des impulsions (5 ns), 1e-12 - pas de temps (1 ps))


Selon le spectre d'une impulsion gaussienne avec $$$f_s = 1 \ GHz$ déterminer qu'à -40 dB la fréquence $$$f \ environ 3 \ GHz$ .
Dans cet exemple, le milieu dans lequel se trouve l'objet sondé, on sélectionne du sable sec avec une permittivité relative $$${\ epsilon} _r = 3$ .
Trouvez la vitesse de propagation de l'onde électromagnétique dans le sable:

$$

$$v = {30 \ over {\ sqrt {3}}} = 17,3 \ cm / ns$$

Définissez la longueur d'onde dans le sable:

$$

$$\ lambda = {v \ over f} = {17,3 \ over {3}} = 5,8 cm = 58 mm$$

Sur cette base, nous sélectionnons le pas qui est le même pour tous les axes ( $$$\ Delta = {\ Delta} _x = {\ Delta} _y = {\ Delta} _z$ ) et égal à 2 mm = 0,002 m pour des raisons de commodité (un nombre entier de marches s'inscrit dans 1 cm):

 #dx_dy_dz: 0.002 0.002 0.002 

Limitez la zone simulée à un rectangle de largeur $$$w$ égal à 80 cm = 0,8 m et hauteur $$$h$ égal à 60 cm = 0,6 m:

 #domain: 0.60 0.60 0.002 

(pour un modèle à deux dimensions, il est nécessaire d'indiquer une épaisseur égale à un pas (0,002))

La taille de la zone de simulation et la taille du pas spatial déterminent le nombre de cellules du modèle et, par conséquent, les besoins en mémoire de l'ordinateur.

Nous décrivons le sable avec une conductivité spécifique $$$\ sigma = 0,01 \ \ cm / m$ et constante diélectrique relative $$${\ epsilon} _r = 3$ commande:

 #material: 3 0.01 1 0 sand 

(1 correspond à la perméabilité magnétique relative $$${\ mu} _r$ égal à l'unité (pas de propriétés magnétiques), 0 - pas de perte magnétique, et étiquette de sable de ce matériau).

Remplissez le sable avec la majeure partie de la zone simulée (de y = 0 à y = 38 cm = 0,38 m):

 #box: 0 0 0 0.80 0.38 0.002 sand 

(0 0 0 - coordonnées du coin inférieur gauche, 0,80 0,38 0,002 - coordonnées du coin supérieur droit (0,002 - pas d'échantillonnage))
Le reste est de l'espace libre (label free_space), presque équivalent en propriétés de l'air ( $$${\ epsilon} _r = 1$ , $$${\ mu} _r = 1$ , $$${\ sigma} = 0$ )
Les limites de la zone de simulation sont présentées sous la forme de conditions aux limites absorbantes (ABC).

Comme cible, nous choisissons un cylindre parmi un conducteur idéal (réfléchissant complètement le rayonnement électromagnétique) avec un rayon de 6 cm = 0,06 m avec un centre situé en un point avec les coordonnées x = 25 cm = 0,25 m et y = 10 cm = 0,1 m:

 #cylinder: 0.250 0.1 0 0.250 0.1 0.002 0.06 pec 

(le pec est un matériau parfaitement conducteur)

Antennes

Le GPR simulé est équipé de deux antennes - émission et réception.
Dans notre étude de cas, nous représentons une antenne émettrice avec un dipôle Hertz dont la longueur est égale à l'étape d'échantillonnage (dans le cas tridimensionnel, vous pouvez sélectionner une antenne à partir d'une vaste bibliothèque) gisant dans le sable (travaillant en contact avec le milieu sondé) à une distance de 5 cm à gauche du milieu de la région (x = 35 cm = 0,35 m, y = 38 cm = 0,38 m):

 #hertzian_dipole: z 0.35 0.380 0.0 pulse 

(z est l'axe de polarisation dipolaire (pour le cas bidimensionnel (mode 2D TMz) seul z est valide), pulse est l'étiquette de la forme de l'impulsion rayonnée par l'antenne)

L'antenne de réception est généralement située à une petite distance constante de la réception, qui est appelée la base de l'unité d'antenne (cette option pour la position relative des antennes est appelée "décalage commun"). Comme base, sélectionnez une distance de 10 cm, donc la coordonnée horizontale est 35 + 10 = 45 cm = 0,45 m (5 cm à droite du milieu de la région):

 #rx: 0.45 0.38 0.0 

Intervalle de simulation

Le choix de la fenêtre temporelle pour la modélisation est déterminé de sorte que le signal réfléchi par la cible ait le temps d'atteindre l'antenne de réception.

Nous déterminons le temps approximatif requis par le signal dans le cas considéré, en prenant la distance des antennes à la cible $$$h = 18 \ \ cm$ :
$$$t \ approx {2 \ cdot h \ over {v}} = {2 \ cdot 18 \ over 17.3} = 2.1 \ \ ns$
Étant donné que le sommet de l'impulsion gaussienne du radar avec une fréquence centrale de 1 GHz est décalé par rapport au début de l'axe du temps de 1 ns, nous sélectionnons une fenêtre temporelle de 5 nanosecondes:

 #time_window: 5e-9 

Modélisation

Ainsi, le contenu du fichier d'entrée est le suivant:


Nous commençons le processus de modélisation:

 python -m gprMax models\hello.in 

Pour effectuer la simulation, la méthode du domaine temporel à différence finie (FDTD, domaine temporel à différence finie) est utilisée (l'algorithme de base de la méthode a été proposé par Kane Yee), après quoi les cellules dans lesquelles le modèle est divisé sont nommées cellules Yee . Une solution numérique est obtenue dans le domaine temporel en résolvant les équations de Maxwell pour chaque cellule.

Dans le cas bidimensionnel (mode 2D TMz), seule la composante est calculée $$$E_z$ champ électrique et composants $$$H_x$ et $$$H_y$ champ magnétique.

Si la quantité de mémoire disponible est dépassée, un message sur le manque de mémoire est émis pour effectuer la simulation:


Si l'un des objets sort du cadre de la simulation, un message d'erreur s'affiche:


Pour effectuer la simulation avec le modèle décrit, il a fallu seulement environ 56 Mo de RAM (si vous réduisez le pas de moitié - jusqu'à 1 mm - les besoins en mémoire augmentent à 99 Mo).

Une fois la simulation terminée, le fichier hello.out apparaît dans le dossier des modèles , contenant les résultats de la simulation au format HDF5 , spécialement conçu pour stocker des données numériques.

Construction de voies

Pour visualiser les résultats, nous construisons les traces:

 python -m tools.plot_Ascan models\hello.out 

Chaque piste (A-scan) présentée dans la fenêtre qui s'ouvre affiche un graphique temporel d'une des composantes du champ électromagnétique à l'emplacement de l'antenne de réception:


Une onde directe provenant directement de l'antenne d'émission (DW) et une onde réfléchie par la cible (RW) sont visibles sur les trajets.

L'axe horizontal du temps est lié à la profondeur de la cible reflétant le signal à travers la vitesse de l'onde électromagnétique dans la substance.

Mais que se passe-t-il si vous placez la commande cylindre devant la commande box dans le fichier d'entrée?


Le signal réfléchi par le cylindre a disparu - le sable a absorbé le cylindre (c'est un exemple de l'importance de l'ordre dans lequel les objets sont construits).

Construction de profil

Mais plus informatif est le radarogramme (radargram) - profil (B-scan), qui est une combinaison de nombreuses pistes construites lors du déplacement du géoradar dans une direction donnée - c'est la procédure même pour déplacer un chariot avec un radar le long de la zone étudiée.

Modifiez la description des antennes en les déplaçant au début de l'axe horizontal:

 #hertzian_dipole: z 0.10 0.38 0.0 pulse #rx: 0.20 0.38 0.0 

Nous fixons le pas de déplacement des antennes égal à 1 cm = 0,01 m:

 #src_steps: 0.01 0 0 #rx_steps: 0.01 0 0 

Ainsi, le contenu du fichier d'entrée est le suivant:


Exécutez la simulation en mode batch:

 python -m gprMax models\hello.in -n 50 

(50 est le nombre de pas que le radar déplace).

Après le démarrage, la simulation est effectuée séquentiellement pour 50 positions GPR:


Après la fin de la simulation, il y a 50 fichiers hello1.out ... hello50.out dans le dossier modèle.
Combinez ces fichiers dans le fichier hello_merged.out avec la commande:

 python -m tools.outputfiles_merge models/hello 

Créez un profil:

 python -m tools.plot_Bscan models\hello_merged.out Ez 

(Ez est la composante du champ électromagnétique par laquelle nous construisons le profil - la composante directement convertie en tension)


Comme vous pouvez le voir, la barre horizontale causée par l'onde directe est affichée sur le profil ci-dessus, et l'hyperbole caractéristique causée par l'onde réfléchie et montrant le changement de la distance de la cible au GPR lors du déplacement est affichée ci-dessous.

La légende de droite montre les couleurs et les intensités de champ correspondantes. $$$E_z$ (montré sur la piste):
rouge - $$$E_z> 0$
blanc - $$$E_z = 0$
bleu - $$$E_z <0$
En analysant de tels profils, nous pouvons tirer des conclusions sur la profondeur, la taille et même la forme des cibles, et des réseaux de neurones sont également utilisés pour cela.


Itinéraires et profil sur l'écran GPR (cadre de l'émission de télévision britannique «Command of the Time»)

Mais la réflexion se produit non seulement à partir d'objets conducteurs, mais aussi à la frontière de deux couches avec des constantes diélectriques différentes.

Créer une deuxième couche de sable avec une constante diélectrique dans la partie inférieure du modèle $$${\ epsilon} _r = 9$ :



Comme on peut le voir, sous la «trace» de l'onde directe (DW), un segment linéaire de l'onde (RW) réfléchi par l'interface des deux couches de sable est apparu.

En dehors du cadre de cet exemple, les fonctionnalités de gprMax telles que la modélisation tridimensionnelle, les formes de surface complexes, les modèles d'antennes détaillés, expliquant la dispersion des ondes électromagnétiques ... De plus, le package peut être utilisé non seulement pour simuler le GPR, mais aussi pour étudier la propagation des ondes électromagnétiques dans divers environnements, et accélérer les calculs à l'aide de la technologie NVIDIA CUDA , qui offre une multiplication par dix de la vitesse par rapport au calcul parallèle sur un processeur utilisant OpenMP. Pour augmenter la flexibilité de la modélisation, vous pouvez placer des blocs de code en Python dans le fichier d'entrée.

Quelques exemples d'utilisation du package gprMax:

• évaluation de l'effet des propriétés du sol sur l'efficacité du GPR dans la recherche de mines antipersonnel
• vérifier la possibilité d'utiliser un radar pénétrant dans le sol pour une détermination non destructive de la profondeur de plantation des semences de maïs
• l'amélioration de la qualité d'interprétation des données de terrain GPR dans l'étude de la résistance des structures des passages souterrains à travers l'autoroute
• étude de la relation entre l'amplitude du signal réfléchi et le diamètre des armatures en béton
• évaluation de l'aptitude du GPR à détecter les vides dans les troncs d'arbres

Liens utiles:

Site officiel de gprMax
Guide d'utilisation de GprMax
gprMax sur YouTube