Zona de estudio

En este trabajo se propone analizar la zona geológicamente relacionada con el Granito de Durania en Norte de Santander (Figura 1). Esta área tiene un interés geológico debido a su potencial minero-económico asociado (^{Mantilla et al., 2009}, 2011) y por la presencia conocida de zonas de difícil acceso, boscosas y de vegetación espesa, topografía irregular y cuerpos de agua, en la cual es clave una correcta definición de los puntos para la medición de datos gravimétricos y magnéticos.

Figura 1 Mapa geológico de la zona de estudio con una extensión de 269,8 km². Geográficamente cercano al municipio de Durania, Norte de Santander. Adaptado de ^{Fúquen et al. (2010)}. 

Geológicamente este cuerpo corresponde a un granito blanco moscovítico (^{Ward et al., 1973}) originado por la inyección de magmas mantélicos en la corteza. La edad de este cuerpo se ha determinado por el método U-Pb en circones, resultando 442,6 Ma (+7,4/-6,0 Ma) (^{Botello-Becerra, 2014}) correspondiente al Triásico-Jurásico. Adicionalmente, se considera que el Granito de Durania intruyó el basamento metamórfico del Macizo de Santander (^{Alarcón et al., 2020}) y se emplazó durante el Cretácico (^{Rojas-Barbosa, 2013}) en las formaciones sedimentarias. Estas consideraciones geológicas sugieren, además del potencial extractivo mencionado, altos contrastes en las propiedades físicas (densidad y susceptibilidad magnética) entre la litología intrusiva y la sedimentaria que requieren de una geometría de exploración adecuada que permitan una correcta interpretación.

Clasificación de coberturas y predicción de variables del terreno

Clasificación de coberturas: la clasificación de coberturas se realizó usando una imagen Sentinel-2 y EUROSAT aplicando aprendizaje profundo a través de la potencialidad que tienen estos métodos en la tarea de clasificación, y especialmente con el uso de redes convolucionales sobre imágenes (^{Guo et al., 2017}). Para Sentinel-2, se consideró una adquisición con cobertura de nubes del 3% con fecha 4 de febrero de 2021; disponible en el portal Earth Explorer del Servicio Geológico de Estados Unidos (Figura 2) (Copernicus Sentinel Data, 2021) y el conjunto de imágenes Sentinel-2 clasificadas EUROSAT (22000 imágenes) (^{Helber et al., 2018}, 2019) donde hay imágenes dentro de las siguientes etiquetas: cultivo anual, bosque, vegetación herbácea, pasto, cultivo permanente, residencial, río y mar-lago.

Figura 2 División de la imagen satelital en ventanas de clasificación. Las ventanas de clasificación en bordes blancos se delimitaron según las dimensiones del conjunto de datos de entrenamiento EUROSAT, cuyas dimensiones eran de 64 x 64 píxeles y 3 bandas (R, G, B) para un total de 667 ventanas. 

En esta investigación se usó el aprendizaje profundo para analizar la favorabilidad de los siguientes factores condicionantes: 1) las coberturas vegetales debido al efecto multitrayecto que generan zonas espesas sobre el posicionamiento GPS y a las zonas de propiedad privada donde el acceso sería restringido, 2) las características del terreno como la inestabilidad, la irregularidad y la inclinación de zonas que dificultarían el acceso y afectarían la calidad del dato, y 3) los cuerpos de agua debido a la carencia de equipos adecuados para adquisiciones en medio acuático y sus efectos sobre la calidad del dato.

Para la clasificación de coberturas se implementó una Red Neuronal Convolucional (CNN), donde se incluyeron las imágenes de EUROSAT como datos de entrenamiento. Así, y a partir de pruebas de verificación y evaluación hechas durante la arquitectura, se determinó un 75% (16500) de las imágenes para el grupo de entrenamiento (subdividido a 12375 de entrenamiento y 4125 de validación) y, el 25% (5500) restante, se consideró como grupo de prueba. El diseño de la red se hizo de tipo secuencial (un solo valor de salida para un valor de entrada) donde se incluyó un optimizador de tipo RMSprop, una función de pérdida tipo categorical Crossentropy y una métrica de precisión accuracy. Estos parámetros permiten a la CNN desarrollar una clasificación multiclase como se describe en ^{Guo et al. (2017)} y ^{Sultana et al. (2018)}; y se configuró de acuerdo con el mejor desempeño obtenido siguiendo los detalles de implementación mostrados en ^{Developers (2022)}. En la Figura 3 se muestran los tamaños y cantidad de filtros, además de cada operación convolucional dentro de la red, finalizando con una capa de salida completamente conectada (fully-connected).

Variables del terreno: para llevar a cabo la predicción de las variables asociadas al terreno se implementó una Red Neuronal Artificial Perceptrón Multicapa (RNAPM) como se muestra en la Figura 4, la cual contiene una capa de entrada con 5 neuronas, una capa oculta con 12 neuronas que usan la función de activación ReLu, y una capa de salida con la función de salida Sigmoid, para los cálculos del gradiente y su propagación hacia atrás se usó el optimizador Adam (^{Kingma y Ba, 2015}), y la función de pérdida Binary Crossentropy con una métrica de precisión accuracy. Esta configuración permite a la red clasificar a la zona de estudio entre 0 y 1, donde los valores superiores a 0 corresponden a las zonas donde se consideraría favorable la adquisición. Esta arquitectura de red completamente conectada se diseñó siguiendo las definiciones y parámetros trabajados en ^{Goyes-Peñafiel y Hernandez-Rojas (2021a}, 2021b).

Para el entrenamiento se usaron 500 coordenadas con datos extraídos del modelo de elevación digital DEM (^{ASF DAAC, 2015}). Sobre esos datos se consideró una división de la siguiente forma: 60% datos de entrenamiento y 40% de validación. Esta división es producto de la verificación y evaluación de rendimiento computacional de la red durante las pruebas. Para los datos extraídos por las coordenadas, se consideraron 5 variables asociadas al terreno: pendiente, rugosidad, orientación, curvatura tangencial y curvatura de perfil. La incorporación de estas variables permite eliminar el sesgo univariable asociado a la predicción. Asimismo, previo al entrenamiento de la red, se realizó una clasificación binaria sobre la variable pendiente (Figura 5) para la solución del diseño de adquisición. La clasificación binaria se realizó a partir de pruebas de verificación y evaluación donde se estimó que las pendientes favorables estaban en un rango inferior a 10° y que están relacionadas con la ejecución real de proyectos de adquisición en zonas montañosas, donde la logística se ve afectada principalmente por la pendiente del terreno.

Cuerpos de agua: el análisis de factores condicionantes presentes en la zona de estudio se complementó con la ubicación e identificación de cuerpos de agua presentes en la zona. Así, se consideró la información georreferenciada disponible desde el Instituto Geográfico Agustín Codazzi (IGAC) (planchas 98IB, 98ID, 98IIA, 98IIC, 98IIIB y 98IVA) donde se delimitaron los ríos Pamplonita y Zulia, al igual que dos lagunas no nombradas (^{Instituto Geográfico Agustín Codazzi, 2018}). Con esta información georreferenciada se complementa y fortalece la estimación de la favorabilidad de puntos de adquisición geofísica implementada en la clasificación resultante. Así, se garantizará que las coordenadas de adquisición no coincidan con la ubicación de los cuerpos de agua.

Figura 3 Arquitectura de la Red Neuronal Convolucional. La función de activación empleada en las capas ocultas es de tipo ReLu, mientras que la de salida es la softmax. 

Figura 4 Modelo de Red Neuronal Artificial Perceptrón Multicapa para predecir la favorabilidad de la medición en puntos de adquisición a partir del análisis de terreno con 5 variables de entrada. 

Figura 5 Puntos de entrenamiento sobre la variable pendiente del terreno. 

Evaluación del diseño de adquisición Para evaluar las geometrías obtenidas, se usaron los puntos de muestreo sobre una adquisición simulada en los levantamientos de gravimetría y magnetometría con la librería SimPEG (^{Cockett et al., 2015}). La construcción del modelo 3D se realizó con una malla refinada sobre el cuerpo intrusivo, esta estrategia permite optimizar la discretización de modelos geologicalmente complejos (^{Torres y Goyes-Peñafiel, 2013}). Posteriormente, y siguiendo un flujo típico en la interpretación de datos geofísicos, se aplicaron técnicas de interpolación para generar imágenes geofísicas que cubran el área de estudio y determinar anomalías asociadas al levantamiento. Las técnicas de interpolación usadas fueron Kriging Ordinario (OK), Ponderación de Distancia Inversa (IDW), y un método basado en machine learning llamado Bosque Aleatorio (RF). Así, mediante los algoritmos de la librería PyKrige en Python (OK), el módulo de IDW en QGIS y la librería Sklearn de Python (RF), se ejecutaron las interpolaciones de datos simulados. Cabe destacar que para IDW y RF se establecieron los valores por defecto para la ejecución de los códigos. Los resultados finales se evaluaron mediante tres métricas de error que relacionan al dato simulado con el resultado de la interpolación. Las tres métricas de error implementadas fueron Error Absoluto Medio (MAE), Error Cuadrático Medio (MSE) y coeficiente de determinación (r²).

Clasificación de coberturas y predicción de variables del terreno

Una vez entrenadas la CNN y la RNAPM, se obtuvieron las clasificaciones correspondientes para la cobertura del suelo y la predicción de la favorabilidad en función de las variables del terreno. De esta forma, el etiquetado de coberturas obtenido arrojó zonas de bosque, pasto, vegetación herbácea, residencia, cultivos anuales y cultivos permanentes (Figura 6).

Figura 6 Clasificación de la cobertura del suelo con la red neuronal convolucional e imágenes satelitales. 

La predicción realizada en función de las variables del terreno (pendiente, rugosidad, orientación, curvatura tangencial y curvatura de perfil) se ejecutó en 200 épocas con una exactitud de prueba 0,87 y una pérdida de prueba de 0,3. La predicción arrojó un ráster multibanda de 1791305 píxeles donde se encontraron tres rangos de favorabilidad: 1) píxeles con un valor mayor a 0,1507 favorables, 2) entre 0,027 y 0,1507, medianamente favorables y 3) menores a 0,027 zonas de baja favorabilidad para la adquisición (Figura 7).

Para el diseño final de la geometría de adquisición, se trazaron 31 transectas sobre zonas donde surgiera una continuidad entre probabilidades medias y altas (amarillo y verde, respectivamente) de tal manera que se lograra una extensión de terreno aceptable para realizar las adquisiciones de los métodos de gravimetría y magnetometría; además de evitar los cuerpos de agua georreferenciados. Estas transectas suponen recorridos para ambos métodos con el fin de optimizar la logística de adquisición. Sobre las transectas trazadas, se puntualizaron zonas de adquisición cada 100 m generando un total de 234 puntos sugeridos para las adquisiciones (Figura 8). Durante el trazado de puntos se ubicaron 5 transectas sobre el afloramiento del cuerpo intrusivo (transectas 10, 15, 16, 20 y 26) de acuerdo con los puntos de acceso y la favorabilidad. En las zonas donde no aflora el Granito de Durania, se trazaron 26 transectas en función de los puntos de acceso, cercanía a los municipios y la favorabilidad de la zona. Por ejemplo, al oeste de Durania se trazaron 7 transectas, mientras que al este otras 5. Al norte de Chinácota se ubicaron 2, mientras que en cercanías de Bochalema se trazaron 12. Cabe mencionar que las vías de acceso están trazadas, en su mayoría, siguiendo la trayectoria del Río Pamplonita y conectando a los municipios. Por esto, la mayoría de las transectas buscan las zonas de favorabilidad en cercanía a los municipios.

Figura 7 Análisis binario del terreno. Las zonas verdes corresponden con una favorabilidad alta en la medición de datos puntuales, y las zonas rojas indican una mayor restricción. 

Evaluación del diseño: simulación e interpolación espacial de datos

La evaluación del diseño se realizó mediante la simulación de una adquisición de gravimetría y magnetometría sobre la zona de estudio. Para esto, se utilizó la geometría obtenida (Figura 8) y se elaboró un modelo tridimensional del Granito de Durania para evaluar el impacto de la limitada cantidad de puntos de adquisición diseñados respecto al área de estudio y realizar una interpolación de mediciones geofísicas usada con frecuencia para la elaboración de mapas asociados a la interpretación de datos geofísicos.

El modelo 3D de la Figura 9 consiste en la construcción computacional del Granito de Durania donde se incluyen las propiedades físicas de interés como la densidad y susceptibilidad magnética, el cual está basado en el corte geológico de la plancha 98 (^{Fúquen et al., 2010}) que se muestra en la Figura 1. La reconstrucción computacional y la simulación de datos se realizaron con la librería SimPEG (^{Cockett et al., 2015}). Con el fin de establecer un contraste sobre el cuerpo intrusivo, se establecieron dos litologías en la zona: una sedimentaria (densidad: 2,42 y susceptibilidad magnética de 0,0001 SI; ^{Lowrie, 2007}) y otra intrusiva (2,61 y 0,01 SI; Lowrie, 2007), con un contraste positivo de densidad de 0,19 g/cm³. Adicionalmente, para la simulación de la respuesta del campo magnético se asumió la dirección de un campo inductivo con inclinación (31,06°), declinación (-9,25°) e intensidad horizontal (27117 nT), correspondientes a la ubicación geográfica del área de estudio (^{NCEI Geomagnetic Modelling Team, 2009}).

Figura 8 Transectas de adquisición. Las zonas poco favorables para la adquisición se muestran en color rojo con probabilidad < 0,0276 y las cuales se excluyen completamente del análisis. En la parte inferior se observa un acercamiento a una de las zonas y se muestra un ejemplo del resultado de dividir una transecta en puntos separados cada 100 m. 

Figura 9 Modelo de densidades y susceptibilidades en los cortes z=-500, y=0, x=2500. La discretización usada corresponde a una malla refinada sobre el Granito de Durania. 

De la simulación se obtuvo un conjunto de datos de gravimetría y magnetometría (Figuras 10A y 10B) en las posiciones de la geometría obtenida por el método propuesto sobre los 234 puntos de muestreo. Además, la medición fue contaminada con ruido gaussiano con desviación estándar de 0,001 mGal más 1% de la medición máxima para simular condiciones de campo típicas de este tipo de adquisiciones. Posteriormente, se llevó a cabo la interpolación de las mediciones aplicando los métodos de Kriging Ordinario (OK) (Figuras 11A, 11B), Ponderación de Distancia Inversa (IDW) (Figuras 11C, 11D) y Bosque Aleatorio (RF) (Figuras 11E, 11F). En el caso de OK se implementó un modelo de semi-variograma esférico y un número de intervalo promedio de 30, mientras que para IDW se utilizó el coeficiente de distancia de 2 y un número de píxeles de 250. Para RF se usaron los valores predeterminados en la librería de Python.

Figura 10 Datos geofísicos simulados con la geometría de adquisición obtenida a partir del análisis de restricciones ambientales y topográficas. A. Anomalía gravimétrica y B. Anomalía magnética de campo total. 

Figura 11 Interpolación de datos simulados. A y B corresponden a la interpolación con OK; C y D a los resultados con IDW; E y F representan la interpolación con RF. 

El método de interpolación que obtuvo mejor métrica fue el IDW donde los valores promedio fueron 0,047 (MAE), 0,01 (MSE) y 0,99 (r ²). Por su parte, el método OK obtuvo valores de 0,0474 (MAE), 1,7 (MSE) y 0,989 (r ²); mientras que el método RF obtuvo 3,95 (MAE), 28,95 (MSE) y 0,76 (r ²).

En la Figura 11 se destaca una zona con valores de 2,5 a 5,0 mGal en la zona central del afloramiento del Granito de Durania para la anomalía gravimétrica (Figuras 11A, 11C y 11E). Se observa que la zona central, donde aflora el cuerpo intrusivo, tiene una anomalía más alta que el entorno sedimentario que lo rodea, coincidiendo con el contraste positivo (Δρ = 0,19 g/cm ³ ) resultante de la asignación de los valores de densidad para las dos litologías (Δρ = ρ ₂ - ρ₁; Δρ = 2,61 g/cm³ -2,42 g/cm³). Por su parte, la anomalía de campo total en el caso de magnetometría (Figuras 11B, 11D y 11F) muestra un comportamiento de tipo dipolo entre la zona central y sur del afloramiento, sugiriendo una unidad magnética amplia.