Inversion des données géophysiques de champ potentiel pour la modélisation 3D et l'estimation des réserves (Exemple de la mine de Hajjar, massif de Guemassa, Maroc) : cas des données magnétiques et gravimétriques

Comptes Rendus. Géoscience, Tome 352 (2020) no. 2, pp. 139-155.

Saâd Soulaimani ; Saïd Chakiri ; Ahmed Manar ; Ayoub Soulaimani ; Abdelhalim Miftah ; Mustapha Bouiflane

(traduction de DeepL post-éditée par N. Bacaër, suggestions d'amélioration : nicolas.bacaer@ird.fr)

Résumé

L'inversion des données géophysiques est un outil qui peut être utilisé pour récupérer la distribution des propriétés physiques sous la surface à partir de données de terrain. Chaque type de données géophysiques peut être inversé en utilisant un ou plusieurs algorithmes d'inversion. Dans cet article, un ensemble de données géophysiques magnétiques et gravimétriques de la région de Hajjar au Maroc, couvrant une étendue de 3,2×1,6km2, ont été utilisés pour réaliser un modèle 3D d'un gisement et pour estimer la réserve minérale par inversion des données géophysiques du champ potentiel et estimation de la masse excédentaire. Nous encourageons ainsi le développement et l'application de l'inversion des données géophysiques du champ potentiel à l'aide des logiciels Geosoft Oasis Montaj et Voxi Earth Modelling™ et l'évaluation de sa puissance par rapport à la méthode d'estimation de la masse excédentaire. Le processus d'inversion commence par le traitement des données, puis passe à l'analyse et à l'interprétation, et se termine par l'inversion de cellules de coupe cartésiennes sans contrainte. Les résultats montrent une variation de -0,22 mGal à 1,59 mGal pour la carte d'anomalie résiduelle de gravité, conduisant à des variations de densité de 2,45g/cm3 à 4,22g/cm3, et une variation de -232 nT à 1018 nT pour la carte d'anomalie magnétique réduite.

De plus, l'inversion des données a permis de créer un modèle 3D du gisement et de la formation géologique adjacente, et d'estimer les différents paramètres caractérisant le gisement à partir des résultats de l'inversion, qui ont été confirmés par les données de l'enquête : (profondeur ≈160 m ; densité maximale ≈4.22g/cm3 ; densité minimale ≈3g / cm3 ; densité moyenne ≈3.61g/cm3 ; épaisseur de la couche de couverture ≈120 m ; immersion ≈45° ; morphologie ≈ objectif ; volume ≈4.8×106m3).

Il a donc été possible d'évaluer la réserve, et de valider la fiabilité de l'inversion en ayant une erreur quadratique moyenne entre la réserve exploitée et la réserve calculée de 13,5%, c'est-à-dire une différence insignifiante entre les réponses réelles et calculées des corps minéralisés magnétiques et gravimétriques, ce qui confirme la validité des résultats.

1. Introduction

L'inversion des données magnétiques et gravitationnelles est utilisée pour estimer la susceptibilité et la distribution de la densité dans le sous-sol. Plusieurs techniques d'inversion existent pour les champs potentiels [Boulanger et Chouteau 2001, Camacho et al. 2000, Li 2001, Li et Oldenburg 1996, 1995, Moraes et Hansen 2001]. L'inversion des données magnétiques et gravitationnelles peut être considérée comme linéaire ou non linéaire. Dans certains cas, il s'agit seulement de rechercher les susceptibilités et densités des blocs, car les géométries sont fixes. Dans d'autres cas, la question est de déterminer les géométries [Camacho et al. 2000, Montesinos et al. 2003], toutes ces méthodes d'inversion peuvent conduire à des approches linéaires et non linéaires selon le choix de la méthode ou du problème à résoudre. L'autre difficulté rencontrée dans l'inversion est la non-unicité de la solution, c'est un problème inhérent à toutes les inversions.

L'algorithme de minimisation utilisé dans la majorité des cas provient de l'algorithme développé par Beiner [1970]. Le principe est de construire une fonction d'erreur. A chaque itération de l'inversion, le modèle géologique est calculé. La fonction d'erreur évalue la différence entre les données mesurées et la réponse du modèle. Tant que cette différence reste importante, les paramètres du modèle sont modifiés l'un après l'autre jusqu'à ce qu'un ensemble de paramètres minimisant la fonction d'erreur soit trouvé. De cette façon, cette approche est répétée pour un certain nombre d'itérations jusqu'à ce que l'on obtienne le modèle dont la réponse satisfait le mieux les données. Ce principe peut être comparé à l'analyse d'une surface topographique [Beiner 1970, Fischer et Le Quang 1981] à la recherche du point le plus bas. Chaque itération permet de descendre le long de cette topographie et de converger vers l'objectif qui représente le modèle final. Un terme de lissage est ajouté à l'approche précédente pour éviter que le modèle ne présente des contrastes irréalistes, cependant, l'approche utilisée dans notre travail était basée sur une inversion linéaire. Les techniques d'inversion linéaire sont formées par les méthodes dites de gradient. Le problème à résoudre peut être réduit de manière significative par la linéarisation du problème en avant. Par la linéarisation, les fonctions de réponse peuvent être grandement simplifiées. Les méthodes d'inversion linéaire sont basées sur la solution d'ensembles d'équations linéaires, qui peuvent être résolues relativement rapidement. Ces méthodes prédominantes sont utilisées pour plusieurs problèmes géophysiques, le point principal de ces méthodes est la linéarisation du problème à terme.

L'inversion à structure minimale, basée sur les voxels, consiste à inverser directement les données gravitationnelles pour récupérer un modèle à structure minimale, la Terre est modélisée en utilisant un grand nombre de cellules rectangulaires de densité constante, et la distribution de densité finale est obtenue en minimisant une fonction objectif du modèle sous réserve d'ajustement des données observées. La fonction objectif du modèle a la flexibilité d'incorporer des informations préalables et ainsi le modèle construit non seulement s'adapte aux données mais aussi est en accord avec des contraintes géophysiques et géologiques supplémentaires. Une pondération de la profondeur est appliquée dans la fonction objectif pour contrecarrer la décomposition naturelle des grains de sorte que l'inversion donne des informations sur la profondeur. Les applications des algorithmes aux données synthétiques et de terrain produisent des modèles de densité représentatifs des structures réelles [Li et Oldenburg 1998]. La pondération des paramètres de profondeur dépend de la discrétisation du modèle, mais elle est facile à calculer. Lorsque l'on sait que les données gravitationnelles sont produites par des contrastes de densité positifs, on peut les incorporer dans l'inversion et il a été démontré que cela améliore la solution. Les applications de cette méthode à des ensembles de données synthétiques ont produit des modèles de densité représentatifs des structures réelles, et l'inversion des données de terrain a produit un modèle de densité cohérent avec la géologie [Li et Oldenburg 1998].

Afin de traiter les données et d'effectuer une inversion sans contrainte, de nombreuses étapes de traitement ont été réalisées en utilisant Geosoft Oasis montaj et Voxi Earth Modelling™. Elles sont basées sur la cellule de coupe cartésienne pour l'inversion des données du nuage. Le point de départ est l'inversion non contrainte avec traitement et analyse des données, jusqu'à l'inversion complète.

2. Contexte géologique

2.1. Cadre géologique

La province métallogénique de Jebilets-Guemassa est constituée de deux grands complexes miniers. Ces districts sont situés dans la meseta occidentale marocaine, respectivement au sud-ouest et au nord de Marrakech (figure 1A). La zone centrale des Jebilets est caractérisée par un continuum du cycle hercynien [Bordonaro 1983, Huvelin 1977]. Au cours de ce cycle, cette zone a connu une intense altération hydrothermale [Ben Aissi 2008, Essaifi 1995] et une forte activité magmatique [Aarab 1984, 1995, Bordonaro 1983, Huvelin 1977, Jadid 1989]. Ainsi, les intrusions acides et basiques qui forment le plutonisme bimodal subcontemporain [Essaifi 1995], sont composées de corps de plusieurs mètres d'épaisseur, organisés en linéaments NNE "submersibles", parallèles aux structures hercyniennes [Ben Aissi 2008, Essaifi 1995]. Ces intrusions contiennent des métallotectes visibles à la surface, à l'intérieur des calottes de fer, répartis selon trois axes majeurs. Autour des corps magmatiques, des déformations régionales et un métamorphisme de contact se développent [Bordonaro 1983, Huvelin 1977, Zaim 1990], ce qui a profondément modifié les structures primaires et les dépôts associés. Le massif de Guemassa est composé d'un bloc tectonique inclus dans le mort-terrain plioquaternaire. Il est composé de terrains sédimentaires et volcano-sédimentaires, qui appartiennent à l'âge du Haut-Visean-Namurien [Hathouti 1990, Hibti 1993]. Le socle hercynien est recouvert de sédiments tabulaires plio-quaternaires, rattachés au groupe de demi-horses de l'Atlas nord [Zouhry 1999]. Les séries d'acides volcaniques sont bien développées avec des calottes de fer [Hathouti 1990]. Elles contiennent des tufs distaux et proximaux ainsi que des coulées de lave [Hibti 2001, Hibti 1993].
Figure 1. A : Carte structurelle du Maroc montrant les principaux domaines de délimitation. Les flèches indiquent le sens des cisaillements pour les structures variscaines tardives selon [Admou et al. 2018]. B : Carte géologique et structurelle du domaine central de la ceinture hercynienne selon [Admou et al. 2018].

2.2. Minéralisation

Cette partie est tirée de [Jarni et al. 2015], qui décrivent en détail la minéralisation du Hajjar. L'analyse des données géologiques montre qu'il s'agit d'un gisement polymétallique de sulfures massifs volcanogènes, avec des ressources en minerai de 18 millions de tonnes contenant 4 à 10 % de zinc, 2 à 4 % de plomb, 0,4 à 0,8 % de cuivre et 60 ppm d'argent. L'amas de sulfure a une texture dominante à grains de pyrrhotite (90 à 95 %) minéralisée avec de la sphalérite, de la chalcopyrite, du cuivre natif et de la galène argentifère. Ce gisement est allongé NNW-SSE, avec une extension longitudinale de 250 à 500 m, une extension latérale de 50 à 250 m et une épaisseur de 40 à 70 m (figures 2, 3). La minéralisation est logée dans des formations volcaniques et volcano-sédimentaires du Carbonifère. Elle est située immédiatement sur le toit du dôme rhyolitique, avec des zones d'altération claires. Le gisement est couvert par un chapeau de fer, qui présente une zone de cimentation très riche en cuivre natif (40%) et en plomb avec une veine de galène argentée de 10 à 20 cm [Jarni et al. 2015].
Figure 2. Disposition des trous de forage à partir de la surface et géométrie du gisement au niveau -235 m modifiée par rapport à [Hathouti 1990].
Figure 3. Section de la minéralisation de sulfure de Hajjar à partir des données de forage (tableaux 1 et 2).

3. Matériaux, méthodes et calcul

3.1. Recherches géophysiques et collecte de données

En complétant l'étude de la gravité géophysique présentée par Soulaimani et al. [2017], nous avons ajouté des mesures magnétiques géophysiques. D'où l'importance d'utiliser ensemble les données gravimétriques et magnétiques pour localiser les amas de sulfures cachés en profondeur. Après la découverte par le forage HS1 de la minéralisation de sulfure de Hajjar (figures 2, 3), un levé gravimétrique et magnétique détaillé a été effectué par le Service géophysique de la Direction de la géologie du Maroc en mai 1984. L'extension de la zone est de 3,2 km de long et 1,6 km de large. Au total, 1089 stations gravimétriques et 1210 stations magnétiques ont été utilisées (figure 4a). L'espacement entre les stations est de 25 m. Cependant, l'espacement entre les profils est de 200 m, sauf pour les 5 profils centraux, pour lesquels l'espacement est réduit à 100 m (figure 4a). L'acquisition des données a été réalisée à l'aide d'un gravimètre Lacoste et Romberg et d'un magnétomètre à césium, avec une résolution de 0,01 mGal et 1 nT, respectivement.
Figure 4. a : Profils des lignes de levés gravimétriques et magnétiques, b : carte des anomalies magnétiques résiduelles, c : carte des anomalies magnétiques résiduelles réduites au pôle.

3.2. Méthodologie

3.2.1. Données magnétiques

Pour étudier l'anomalie magnétique, le champ IGRF a été supprimé, le champ magnétique résiduel est calculé comme la différence entre le champ total et l'IGRF calculé. Enfin, les données magnétiques résiduelles ont été réduites au pôle en utilisant une inclinaison de 46° N et une déclinaison de -7,40° W. La transformation est effectuée dans le domaine des fréquences en utilisant l'opérateur de filtre [Grant 1972], l'interpolation a été effectuée par la méthode de la grille à courbure minimale.

3.2.2. Données gravimétriques

Les données gravimétriques de Hajjar que nous avons pu traiter comprenaient tous les calculs détaillés effectués par le Service Géologique de la Direction de la Géologie de Rabat. De plus, nous avons créé une base de données qui contient l'anomalie de Bouguer pour une densité régionale de 2,67 g/cm3. Les 2,67 g/cm3 et 2,3 g/cm3 sont les densités couramment utilisées pour les levés gravimétriques au Maroc.

La densité de 2,67 g/cm3 a été utilisée pour la région de Hajjar car elle n'avait pas été utilisée dans les études précédentes, et la correction de Bouguer dépend de la géologie/composition dominante du substratum rocheux de la région (une formation volcano-sédimentaire). De plus, les changements maximums d'altitude dans la région couverte par les données gravimétriques sont de 117 m, ce qui est faible comparé à l'épaisseur de la couverture sédimentaire (124,4 m) et à la profondeur du gisement (158 m) [Bellott et al. 1991]. La densité moyenne de la couverture sédimentaire (environ 2,67 g/cm3) a été utilisée.

Par conséquent, une des premières étapes de calcul consiste à individualiser l'anomalie intéressante, afin de proposer des modèles de densité qui puissent en tenir compte. Pour cela, nous avons d'abord défini le champ régional. Plusieurs techniques peuvent révéler l'apparition d'une anomalie régionale. Pour ce faire, nous avons décomposé le signal par séparation linéaire, ou par une transformée de Fourier fonctionnelle qui représente l'anomalie de Bouguer. Dans le cadre de notre étude, nous nous intéressons à la localisation d'anomalies de quelques dixièmes de mGal. La technique que nous avons utilisée pour le traitement des données a été initialement choisie après avoir testé plusieurs méthodes. Elle consiste en une séparation linéaire. Cependant, le champ résiduel (RES) est obtenu par la soustraction du champ régional obtenu (REG) de l'anomalie de Bouguer (BA) ; similaire à la méthode utilisée par El Azzab et al. [2019], l'interpolation a été réalisée par la méthode de la grille à courbure minimale.

3.2.3. Masse excédentaire (loi de Gauss)/filtre de densité apparente

La loi de Gauss stipule que la masse totale d'une région est proportionnelle à la composante normale de l'attraction gravitationnelle intégrée sur la limite fermée de la région [Jackson 1996]. Parmi les applications géophysiques les plus courantes de ce concept, on trouve l'approximation de l'excès de masse sous une surface sur laquelle la composante normale de la gravité est connue [Hammer 1945, Jackson 1996]. Étant donné une masse volumétriquement contrainte, aucune autre hypothèse concernant sa profondeur, sa forme ou sa distribution de densité n'est nécessaire. Cependant, la masse étudiée doit être petite par rapport aux dimensions globales de la zone d'étude [Jackson 1996]. L'application de la loi de Gauss de cette manière entraîne certaines complications. Les données de gravité sur un plan infini n'existent pas, et par conséquent, la meilleure approximation est d'utiliser un ensemble de données qui s'étend nettement au-delà des limites de toute source d'intérêt majeure. L'isolement complet des sources n'est malheureusement pas possible dans le monde naturel. Pour cette raison, la séparation des anomalies gravitationnelles d'intérêt des autres sources locales et régionales est souvent complexe [Jackson 1996].

La loi de Gauss permet de calculer la masse totale en excès qui produit une anomalie de gravité, quelle que soit la distribution de cette masse. Green a trouvé cette solution à l'excès de masse à l'intérieur de la Terre en utilisant un hémisphère pour la surface hypothétique, S [Siqueira et al. 2017]. Supposons que le sommet de l'hémisphère soit la surface de la Terre (z = 0). L'intégrale de la surface peut alors être divisée entre le sommet et l'hémisphère [Connor et Connor 2009]. En supposant que le rayon de l'hémisphère est extrêmement grand, certains ajustements sont faits avec les limites de l'intégration pour simplifier le problème, en pratique, la masse anormale est trouvée en intégrant numériquement les données de gravité maillées sur une zone [Connor et Connor 2009] : \begin{equation*} {M}=\frac {{1}}{{2}{\pi }{G}}\sum _{{i}={1}}^{N}\sum _{{j}= {1}}^{M}{\Delta }{g}({x},{y})\cdot {\Delta }{x}\cdot {\Delta }{y} \end{equation*} où \(\Delta g(x,y)\) est l'anomalie de gravité, N et M sont le nombre de points de grille dans les directions X et Y, respectivement, et Δx et Δy est l'espacement de grille dans les directions X et Y, et G est la constante gravitationnelle 6,67384 × 10-11 m3.kg-1.s-2.

Figure 5. a : Carte de l'anomalie de Bouguer pour la densité de 2,67 g/cm3, b : carte de l'anomalie de gravité résiduelle, c : carte de densité.
Néanmoins, un filtre de densité apparente calcule la densité apparente du sol qui donnerait lieu au profil de gravité observé. La densité suppose que le profil de gravité est dû à un ensemble de prismes rectangulaires avec un sommet au niveau de l'observation du profil de gravité, un fond à la profondeur t, et une longueur de frappe infinie ; nous devons fournir l'épaisseur du modèle de la terre et la densité du fond [Gupta et Grant 1985] : \begin {equation*} {L}({\omega })=\frac {{\omega }}{{2}{\pi }{G}\cdot ({1}- {e}^{-{t}\cdot {\omega }})} \end {equation*} où t est l'épaisseur en unités terrestres du modèle de la terre, d est la densité de fond dans g/cm3 et G est la constante gravitationnelle.

3.2.4. Inversion par VOXI earth modelling™

La modélisation et l'inversion à base de voxels sont devenues un outil familier au cours des deux dernières décennies, en raison de la réduction spectaculaire de ses coûts de calcul et du fait que l'industrie de l'exploration tente d'interpréter des données géophysiques de plus en plus complexes associées à des cibles plus profondes. La simplicité des représentations de la Terre basées sur des voxels les rend attrayantes. Elles sont faciles à comprendre et à programmer, faciles à afficher et conformes aux architectures informatiques actuelles. Cependant, elles présentent un défaut important : elles ne sont pas bien adaptées pour représenter les caractéristiques géologiques courantes dans les projets d'exploration, car les caractéristiques géologiques sont souvent décrites par des surfaces [Ellis et MacLeod 2013].

VOXI Earth Modelling™ est un module de calcul en nuage et en grappe de Geosoft Oasis Montaj qui permet l'inversion des données géophysiques (par exemple, la gravité, le magnétisme) en 3D. Il utilise un algorithme d'inversion de la cellule de coupe cartésienne (CCC) développé par Ingram et al. [2003]. L'algorithme a été simplifié par Ellis et MacLeod [2013] afin de représenter les surfaces géologiques avec une plus grande précision. Les inversions peuvent être réalisées sans contrainte par Roy et al. [2017]. Les données gravimétriques et magnétiques (données au sol ou aéroportées, mesures du champ total ou du gradient) peuvent être inversées en utilisant différentes méthodes : Barbosa et Pereira [2013], Ellis et MacLeod [2013], Farhi et al. [2016], Ingram et al. [2003].

Figure 6. a : section du modèle magnétique inverse, b : section du modèle de gravité inverse, c : modèle magnétique inverse relatif à la minéralisation, d : modèle de gravité inverse relatif à la minéralisation.
La source des données peut être fournie en format base de données ou grille. Dans notre cas, nous avons utilisé les grilles résiduelles pour les cas ci-dessus, la gravité et les données magnétiques comme grille de capteurs. Nous avons donc été obligés de définir l'élévation des capteurs comme des fichiers de grille pour les capteurs d'élévation et de choisir le nombre d'échantillons par cellule. Au départ, tous les effets tels que le mouvement, le terrain, la marée, le champ régional, etc. doivent être supprimés avant la modélisation pour les deux cas "gravité" et "données magnétiques". Ensuite, nous avons fixé le critère d'ajustement (sensibilité du modèle) pour ajuster les données jusqu'à ce que la différence entre la prédiction du modèle (l'ajustement) et les données mesurées soit en moyenne inférieure au critère d'ajustement. L'augmentation du critère d'ajustement augmente le lissage des résultats de l'inversion. Pour les méthodes d'ajustement, nous pouvons utiliser plusieurs options comme les erreurs absolues, les erreurs relatives ou les fractions d'écart type [Bhattacharyya 1965, Dampney 1969, Dransfield 2007, Pedersen et Rasmussen 1990]. Nous avons choisi la valeur par défaut familière aux géophysiciens, qui est égale à 5% de l'écart type des données. Pour notre étude de cas, l'écart-type des données gravimétriques est égal à 0,008504 mGal et celui des données magnétiques est égal à 8,94 nT.

4. Synthèse et résultats

4.1. Cartes magnétiques

L'anomalie magnétique résiduelle de Hajjar (figure 4b) a une forme bipolaire avec un minimum au nord et un maximum au sud. Son axe longitudinal est-ouest est long de 1,4 km. L'amplitude de l'anomalie magnétique résiduelle est de 654 nT. On peut voir sur la carte de l'anomalie magnétique résiduelle que cette anomalie existe dans un contexte magnétique calme, ce qui a été un facteur favorable à son identification. Afin de pouvoir définir le corps minéralisé et le comparer à l'anomalie de gravité, nous avons calculé l'anomalie magnétique réduite au pôle (RTP).

La carte RTP (figure 4c) révèle une anomalie magnétique unipolaire très importante, qui coïncide avec le gisement présentant une amplitude maximale de 1018 nT. Ce maximum peut être interprété comme étant situé au-dessus du centre du gisement (d'où la raison du positionnement à cet endroit du forage HS1 pendant la phase d'exploration dans les années 80 et 90).

4.2. Cartes gravimétriques

Figure 7. Les cartes d'anomalies résiduelles des gisements. a : gravité, c : magnétique. Accompagnées des réponses des modèles inverses. b : réponse de la gravité, d : réponse magnétique.
Figure 8. Carte des différences entre les données observées et les données calculées, pour les deux. a. Données magnétiques, b. Données gravimétriques.
Nous avons affiché l'anomalie de Bouguer sous forme de carte en couleurs. Cette carte de l'anomalie de Bouguer (figure 5a) montre une anomalie positive nord-sud située au nord-ouest. Cette carte est affichée avec des contours espacés de 0,2 mGal. La forme générale de la carte montre des anomalies négatives orientées approximativement SE-NO avec une augmentation du gradient de gravité du sud-est au nord-ouest. Plusieurs anomalies sont mises en évidence avec des longueurs d'onde et des amplitudes différentes liées aux structures géologiques dont les dimensions et les densités sont variables.

La carte des anomalies de Bouguer ne permet pas de dessiner une interprétation des structures peu profondes car un effet de champ régional est évident. Si cette carte nécessite une interprétation, la grande amplitude de champ positif dans la partie nord-ouest peut s'expliquer par l'effet régional dû aux structures géologiques régionales. Dans notre cas, il peut être interprété comme un épaississement de la couverture sédimentaire dans la partie sud-est et un amincissement dans la partie nord-ouest.

Sur la carte des anomalies de gravité résiduelle (figure 5b), nous avons trouvé une anomalie positive très importante au centre gauche avec une amplitude de 1,59 mGal. Elle a une forme circulaire avec deux allongements : une extension directionnelle E-O et une direction NW-SE bordée par des bandes d'anomalies de moyenne amplitude et moins étroites. Cette anomalie occupe une surface de 333 000 m2. En analysant la carte des anomalies de gravité résiduelle (figure 5b), nous avons pu identifier les anomalies positives.

4.3. Carte de l'excès de masse/densité

Après la partie résiduelle-régionale de l'ensemble de données, en utilisant une séparation linéaire, une zone entourant immédiatement l'anomalie de gravité de Hajjar a été sélectionnée dans le but d'effectuer un calcul de masse excédentaire. Pour effectuer efficacement ce calcul, les données gravimétriques doivent d'abord être interpolées sur une grille régulière, étape déjà effectuée.
Figure 9. Disposition des trois gisements. a : Anomalie magnétique RTP. b : Anomalie de gravité résiduelle.
Une valeur seuil de gravité ("fond") est alors déterminée visuellement à partir de la carte de gravité. Cette valeur de fond, dans le cas de l'anomalie de Hajjar, a été déterminée comme étant d'environ 0,52 mGal. En utilisant un logiciel personnel pour effectuer une intégration numérique des données gravitationnelles maillées à travers la zone d'intérêt sélectionnée, la masse anormale est déterminée en soustrayant la valeur de fond assignée à chaque point de grille qui coïncide avec le corps minéralisé, puis en additionnant toutes les valeurs ensemble. La somme est ensuite multipliée par Δx Δy / 2πG.

Sur la base de la valeur de fond sélectionnée de 0,52 mGal, la masse excédentaire calculée pour l'anomalie de Hajjar est de 7,32 × 109 kg : \begin {equation*} {M}_{1}=\frac {307173.56\times 10^{-5}}{2\pi {G}}~\mathrm {kg} =7.32\times 10^{9}~\mathrm {kg}=7.32~\mathrm {Mt}. \end {equation*} Néanmoins, nous avons réalisé une carte de densité en appliquant le filtre de densité apparente pour convertir les cartes de gravité en cartes de densité (figure 5c) afin de déterminer la distribution de densité dans la région en adoptant une densité de fond de 2,67 g/cm3 et une épaisseur de couche modèle de 25 m (la moyenne de l'espacement de la grille). On observe une densité allant de 2,45 g/cm3 à 4,22 g/cm3 (figure 5c).

4.4. Inversion

La figure 6 montre les modèles obtenus à partir des inversions des données magnétiques et gravimétriques pour l'étendue du gisement dans le secteur de Hajjar en utilisant une maille de modèle de 25 m × 25 m × 12,5 m. La figure 6a montre une coupe transversale du modèle final ainsi que le champ magnétique total (sans RTP) où la source de l'anomalie magnétique est colorée en rouge. En se basant sur la carte du champ magnétique RTP, il a été possible de fixer la susceptibilité de la source perturbatrice indiquée sur la figure 6c.

De même, la figure 6b montre une coupe transversale du modèle de gravité obtenu (densité variant de 2,45 g/cm3 à 4,22 g/cm3) accompagnée de l'anomalie gravitationnelle résiduelle. La figure 6d montre le modèle de gravité pour les densités supérieures à 3 g/cm3.

Afin de valider les résultats, la procédure d'inversion établie par Voxi earth modeling™ comprenait auparavant un algorithme qui permet de calculer la réponse du modèle résultant. Les réponses gravitationnelles et magnétiques ont toutes deux été calculées, elles sont présentées sur la figure 7.

Tableau 1. Résumé des données d'enquête "trous de sonde" [Hathouti 1990]


Tableau 2. Mesure de la densité et de la susceptibilité magnétique des formations géologiques du Hajjar [Hathouti 1990] (note SGN/GPH no 628, J. Corpel, 07.12.1984)
D'après les résultats obtenus pour les deux cas, on observe une concordance presque totale entre eux, soit pour la forme, soit pour les valeurs obtenues. Pour le cas de la gravité (figures 7a, 7b et 8b), on note une petite différence qui ne dépasse pas 0,04 mGal. On note également une petite différence qui ne dépasse pas 12,16 nT pour les données magnétiques (figures 7c, 7d et 8a). On peut en déduire que les modèles obtenus sont valables ; en effet, la réserve calculée à partir du modèle inversé était de 17,32 Mt : \begin {equation*} {M}_{2}=\rho \times {V}=3.61~\mathrm {g}/\mathrm {cm}^{3} \times 4.8\times 10^{6}~\mathrm {m}^{3}\approx 17.32~\mathrm {Mt}. \end {equation*}

5. Discussion

Après le traitement et l'inversion des données, nous avons pu cartographier la densité souterraine et la distribution de la susceptibilité le long des 3 composantes (x,y,z), permettant ainsi d'avoir le modèle 3D du gisement et de la formation géologique adjacente, et les différents paramètres de décélération qui caractérisent le gisement à partir des modèles inversés non contraints. Ceci a été validé après une comparaison avec les valeurs mesurées sur les échantillons d'étude (tableaux 1, 2) et une comparaison avec l'étude de Bellott et al. [1991], (profondeur ≈ 160 m, densité maximale ≈ 4.22 g/cm3, densité minimale ≈ 3 g/cm3, densité moyenne ≈ 3.61 g/cm3, épaisseur du mort-terrain ≈ 120 m, pendage ≈ 45°, morphologie ≈ lentille, volume ≈ 4.8 × 106 m3).

Nous avons donc pu calculer la réserve du gisement avec la méthode de la masse en excès (7,32 Mt) et à partir du modèle inversé (17,32 Mt), ce qui a une différence de 12,68 Mt pour la première méthode et de 3 Mt pour la seconde méthode par rapport à la "réserve réelle" exploitée déclarée par les rapports du ministère de l'énergie et des mines du Maroc sur les notes et mémoires de Michard et al. [2011] et Jarni et al. [2015]. Notre résultat reste valable tant qu'il ne dépasse pas la valeur RMS (Root mean square error) de 20% [Soulaimani et al. 2017], une condition qui a été vérifiée pour le second résultat : \begin {eqnarray*} \mathrm {RMS}_{1}({\% })&=&\frac {|\text {Real Reserve}- \text {Estimated Reserve}|}{\text {Real Reserve}}\\ &=& \frac {|20~\mathrm {Mt}-7.32~\mathrm {Mt}|}{20~\mathrm {Mt}} =63.4\% \\ \mathrm {RMS}_{2}({\% })&=&\frac {|\text {Real Reserve}- \text {Estimated Reserve}|}{\text {Real Reserve}}\\ &=&\frac {|20~\mathrm {Mt} -17.32~\mathrm {Mt}|}{20~\mathrm {Mt}}=13.4\%. \end {eqnarray*} Il ressort de cette discussion que notre approche est valable par rapport à la première méthode, qui ne reste valable que pour certains cas, et qui donne une idée préliminaire de la masse du gisement.

Notez que le gisement réel est formé de 3 corps minéralisés. Dans notre cas, le gisement obtenu comprend les trois gisements, car nous ne pouvons pas les différencier en raison de la faible distance qui les sépare et des effets de la gravité et du champ magnétique potentiel (figure 9). Pour éliminer cet effet, il est recommandé d'utiliser des opérateurs tels qu'un filtre puissant, basé sur des calculs de tenseurs, comme les systèmes invariants et huit horizontaux [Pedersen et Rasmussen 1990].

6. Conclusion

Cette étude nous a permis de mettre en évidence la puissance du traitement des données pour l'inversion cartésienne sans contrainte des cellules coupées en géophysique, pour la modélisation 3D et l'estimation des réserves, en utilisant des données magnétiques et gravimétriques. Cela nous a permis de déterminer les différents paramètres du gisement étudié (profondeur ≈ 160 m, densité maximale ≈ 4,22 g/cm3, densité minimale ≈ 3 g/cm3, densité moyenne ≈ 3,61 g/cm3, épaisseur des morts-terrains ≈ 120 m, pendage ≈ 45°, morphologie ≈ lentille, volume ≈ 4,8 × 106 m3), ainsi que de définir la réserve avec un niveau de confiance acceptable 17,32 Mt avec une erreur d'estimation de la moyenne quadratique de 13,4%, validant ainsi notre approche. L'objectif de cette étude était d'obtenir une estimation détaillée des réserves à partir des données géophysiques et de construire un modèle 3D du gisement, non seulement pour disposer d'un modèle géophysique et géologique tridimensionnel et évaluer les réserves du gisement, mais aussi pour illustrer l'utilité, la puissance, la fiabilité et l'importance de l'utilisation conjointe des données gravimétriques et magnétiques par une inversion cartésienne sans contrainte des cellules de coupe. En fait, cette inversion conjointe est très puissante pour la caractérisation et la modélisation en 2D et 3D de gisements minéraux similaires. Elle peut permettre de réduire les coûts d'exploration.

Remerciements

Les auteurs remercient les deux réviseurs anonymes qui ont contribué par leurs commentaires constructifs à améliorer la qualité de ce travail. Un remerciement particulier pour son aide va également au professeur Jaffal Mohammed de la FST de Marrakech, de l'université Cadi Ayyad.

Annexe

Voir la figure 10.
Figure 10. Élévation régionale.

Bibliographie