Introducción
En la actualidad existen numerosas poblaciones que se encuentran expuestas a la ocurrencia de movimientos de ladera repentinos, que pueda incluir tanto masas de suelos como de rocas, a lo largo de la pendiente de una ladera. A estos movimientos se los conoce como deslizamientos. Los cuales son eventos naturales, que por lo general no se pueden evitar, pero si se puede reducir el daño ocasionado realizando modelos de predicción.
Cuando existen grandes incertidumbres sobre la probabilidad de ocurrencia de un fenómeno, es difícil tomar decisiones sobre las medidas, la ubicación, etc. El coste de un proyecto puede resultar muy alto o se requieren asumir riesgos de características y magnitudes no determinadas. (Moral, 2014)
A nivel regional son diversos los trabajos que se han realizado sobre deslizamiento o también conocidos como movimientos de laderas cuya finalidad es contar con un insumo que sirva de línea base para estudios relacionados con la gestión de riesgos a desastres y el ordenamiento territorial de ciudades, cuencas y municipios a lo largo de América Latina tal y como lo detallan (Quesada-Román, y otros, 2017). Para el caso de México existen documentos como el trabajado por (Paz-Tenorio, 2012) utilizando el método heurístico con base en mapas cualitativos y mejorados en el 2015 aplicando el análisis multicriterio. Así mismo en Ecuador hay evidencias del uso de modelos probabilísticos para modelar la ocurrencia de deslizamientos utilizando sistemas de información geográficos y geoestadística de acuerdo con (Hermosa, Avilés, Almeida, & D´Howitt, 2010).
Para la región centroamericana se cuenta con diversos estudios como para el que detalla (Quesada-Román, 2016) para Costa Rica donde presenta un análisis morfométrico y morfogenético para definir los distintos niveles de susceptibilidad a deslizamientos e inundaciones. En Nicaragua según (Úbeda, 2016) se han desarrollado mapas de susceptibilidad a deslizamientos usando el método heurístico geomorfológico por medio de sistemas de información geográfico. En los últimos años (Reyes & Rivera, 2019) han trabajado en la construcción del mapa de susceptibilidad a deslizamientos para todo El Salvador utilizando la técnica Multivariate Adaptative Regression Splines (MARS).
Para el caso de Honduras (Harp, Castañeda, & Held, 2002) desarrollaron un mapa del inventario de deslizamientos ocasionados en Tegucigalpa producto del paso del Huracán Mitch y también se han identificado zonas susceptibles a deslizamientos haciendo uso del modelo de talud infinito (Suárez & Domínguez-Cuesta, 2020)
El presente informe se presenta la aplicación de un modelo de predicción de susceptibilidad a deslizamientos, el cual fue aplicado para un sector de la zona occidental de Honduras, para el cual se hizo uso de la metodología MARS.
Materiales y métodos
El estudio se desarrolló en el departamento de Ocotepeque tomando como unidad de análisis la Cordillera del Merendón (Mapa 1) que en su parte más alta alcanza casi a los 2,500 m.s.n.m. que se extiende desde el municipio de Mercedes al sur del departamento cruzando diversos municipios hasta llegar a Dolores Merendón.
Dicha cordillera es importante porque dentro de ella se encuentran la Reserva Biológica Güisayote y la Reserva Biológica el Pital con un área de 116.97 km² y 33.24 km² respectivamente, ambas forman parte del sistema de áreas protegidas de la zona occidental de Honduras.
Metodología MARS
Para poder desarrollar el análisis de la regresión se utilizó la metodología MARS (Multivariate Adaptive Regression Splines) (Friedman, 1991). Dado un conjunto de variables predictoras, MARS lo que hace es ajustar un modelo a la forma de una expansión de productos de funciones bases de predictores elegidos durante una estrategia recursiva particionada forward y backwars. MARS es una herramienta flexible que automatiza la construcción de modelos de predicción, permitiendo la selección de variables relevantes, la transformación de las variables predictoras, establecer las interacciones de las variables predictoras, el tratamiento de los valores perdidos y un autotest para protegerse del sobreajuste. Finalmente, puede revelar patrones y relaciones que es difícil, si no imposible, que otros métodos puedan revelar. (Vanegas & Vásquez, 2016)
Inventario de deslizamientos
El inventario de deslizamientos para el área seleccionada consta de 645 puntos, (Mapa 2) los cuales han sido ubicados en las partes altas de los deslizamientos, es decir, en la corona. Dicho inventario fue creado a partir de imágenes satelitales utilizando Google Earth. Es importante mencionar que a nivel de nuestro país no existe un inventario que reúna todos los deslizamientos que han ocurrido, por lo cual de los puntos seleccionados solo se cuentan con la ubicación de estos, no se conoce la fecha en que sucedieron, ni la tipología, es decir, si es flujo o deslizamiento.
Los puntos se encuentran ubicados entre los municipios de Ocotepeque, Sinuapa, La Labor, San Francisco del Valle, San Marcos y Mercedes correspondientes al departamento de Ocotepeque, ubicado en la parte occidental de Honduras. En el siguiente mapa se puede observar la ubicación de los puntos.
Variables utilizadas
La regresión realizada por la metodología MARS, permite el uso de variables continuas y categóricas. Las variables continuas no cuentan con un rango establecido por lo cual se utilizaron para representar la elevación del terreno y los productos derivados del modelo de elevación digital. En cambio, las variables categóricas cuentan con un rango definido, utilizado para representar la geología y la cobertura forestal o uso del suelo ver tabla 1.
Variables continuas
Código | Nombre | Unidad |
---|---|---|
ele | Elevación | metros |
slo | Pendiente | grados |
plc | Curvatura en planta | Radianes/metros |
prc | Curvatura de perfil | Radianes/metros |
twi | Índice topográfico de humedad | metros |
Fuente: Elaboración propia
Todas las variables continuas se generaron a partir de un modelo de elevación digital (DEM) de 12.5 metros, el cual fue obtenido del Centro de Archivo Activo Distribuido de Instalaciones Satelitales de Alaska (ASF DAAC) del cual se descargó el DEM utilizando las imágenes proporcionadas por ALOS PALSAR que es uno de los múltiples recursos cartográficos disponibles dentro de los productos del satélite ALOS de Agencia Japonesa de Exploración Aeroespacial (JAXA), que adquirió imágenes radar entre 2006 y 2011. La elevación se obtuvo a través del DEM utilizado que cuenta con una resolución de 12.5 metros.
Las variables de pendiente, aspecto, curvatura en planta y curvatura de perfil se obtuvieron utilizando la herramienta “Slope, Aspect, Curvature” del programa SAGA GIS (2.3.2), haciendo uso. (Mapa 3, 4, 5 y 6)
del método “9 parameter 2nd order polynom (Zerenbergen & Thorne 1987)” ya que este realiza un Análisis cuantitativo de la topografía de la superficie terrestre. Las unidades del aspecto y de la pendiente fueron calculados en grados.
Para el cálculo del índice topográfico de humedad (Topographic Wetness Index) se utilizó la herramienta “Topographic Wetness Index” de SAGA GIS (2.3.2). Para poder utilizar esta herramienta se utilizan el raster de elevación, el raster de pendiente (slope), raster de Catchment. (Mapa 7)
Area que corresponde a Flow Accumulation, el área de conversión se utiliza “no conversión (area already given as specific catchment area)” y el método TWI se selecciona “Standard”.
Para realizar el cálculo del Catchment Area se generó un DEM que no contenga depresiones para el cual se utilizó la herramienta de “Fill Sinks (Wang & Liu)”. “Este método permitir la creación de modelos hidrológicos de elevación, es decir, no solo para llenar las depresiones, sino también para preservar una pendiente descendente a lo largo de la ruta del flujo” (Wichmann, 2007) se aprecia en la tabla 2.
Variables categóricas
Código | Nombre | Clases de la variable |
---|---|---|
ASP | Aspecto | ASP1 (N: 0 - 22.5 / 337.5 - 360) |
ASP2 (NE: 22.5 - 67.5) | ||
ASP3 (E: 67.5 - 112.5) | ||
ASP4 (SE: 112.5 - 157.5) | ||
ASP5 (S: 157.5 - 202.5) | ||
ASP6 (SW: 202.5 - 247.5) | ||
ASP7 (W: 247.5 - 292.5) | ||
ASP8 (NW: 292.5 - 337.5) | ||
LCL | Geoformas | LCL1 (Strems) |
LCL2 (Midslope drainages) | ||
LCL3 (Upland drainages) | ||
LCL4 (Valleys) | ||
LCL5 (Plains) | ||
LCL6 (Open slopes) | ||
LCL7 (Upper slopes) | ||
LCL8 (Local ridges) | ||
LCL9 (Midslope ridges) | ||
LCL10 (High ridges) | ||
USO | Cobertura del suelo | Reclasificado en 8 categorías |
GEO | Geología | Reclasificado en 7 categorías |
Fuente: Elaboración propia
Las geoformas (Mapa 8) fueron obtenidas a partir del DEM de 12.5 metros, para esto se utilizó la herramienta “TPI Based Landform Classification” de SAGA GIS (2.3.2). Como resultado se obtiene un raster que se muestra a continuación con las geoformas clasificadas en 10 categorías tabla 3.
Código | Categoría |
---|---|
1 | Strems |
2 | Midslope drainages |
3 | Upland drainages |
4 | Valleys |
5 | Plains |
6 | Open slopes |
7 | Upper slopes |
8 | Local ridges |
9 | Midslope ridges |
10 | High ridges |
Fuente: Elaboración propia
Para la generación del raster de uso de suelo se utilizó el mapa de cobertura forestal de Honduras del año 2018 elaborado por el Instituto Nacional de Conservación y Desarrollo Forestal, Áreas Protegidas y Vida Silvestre (ICF), el mapa fue reclasificado en 8 categorías tal y como se muestra en la Tabla 4 en base al tipo de cobertura y el resultado de la reclasificación se puede observar en el Mapa 9.
Reclasificación | Descripción |
---|---|
1 | Bosque Latifoliado Húmedo |
1 | Bosque Latifoliado Deciduo |
1 | Bosque Mixto |
2 | Bosque de Conífera Denso |
2 | Bosque de Conífera Ralo |
3 | Vegetación Secundaria Húmeda |
3 | Vegetación Secundaria Decidua |
4 | Arboles Dispersos |
Fuente: Elaboración propia
Mapa geológico de Honduras.
Para la elaboración del raster de geología se hizo uso del Mapa Geológico de Honduras que tiene una escala de 1:500,000 que fue creado en 1991 por el Instituto Geográfico Nacional de Honduras, la reclasificación se realizó en 7 categorías cuyo detalle se muestra en el Mapa 10 y el detalle de la litología de muestra en la tabla 5, donde se observa que para el area de estudio predominan las rocas piroclásticas andesíticas y riolíticas, volcánicas que pertenecen a la Formación Tmp del Grupo Padre Miguel.
Reclasificación | Formación | Litología |
---|---|---|
1 | Tpm | Rocas piroclásticas andesíticas y riolíticas, volcánicas |
2 | Qal | Cantos rodados, adoquines, grava, arena, barro |
3 | Pzm | Esquisto, filita, gneis, cuarcita, mármol, vetas de cuarzo |
4 | Ky | Lutitas calcáreas, calizas, margas, dolomitas |
5 | Kva | Camas rojas heterogéneas, caliza Jaitique |
6 | Kti | Granito, granodiorita, diorita, tonalita |
7 | JKhg | Lutitas, areniscas, carbones |
Fuente: Elaboración propia
Análisis de multicolinealidad
La multicolinealidad se produce cuando hay relación entre las variables independientes del modelo. Aunque difícilmente habrá una correlación perfecta (no podrían estimarse los parámetros del modelo), la presencia de multicolinealidad inexacta o imperfecta indica que las variables independientes están compartiendo información y a la hora de utilizarlas para predecir una variable dependiente se produce un fenómeno de redundancia: estamos usando varias veces lo mismo para pronosticar algo. Esto se traduce en una mayor imprecisión en las estimaciones. (Universidad Complutense de Madrid, 2014)
Para realizar este análisis de utilizó el paquete “usdm” en R con el objetivo de establecer la dependencia entre las variables. Se consideró el factor de influencia de varianza (Variance Influencer Factot, VIF) como indicador de dependencia, por lo que se tomó un valor aceptable de VIF<10, Si el valor para una variable no cumple con este criterio se elimina de las variables utilizadas para la creación del modelo.
Construcción del modelo
Para realizar la construcción de un mapa que muestre la susceptibilidad a deslizamientos se debe considerar el análisis de cada tipo de deslizamiento por separado, pero debido a que no se cuenta con un inventario de deslizamiento que cuente con dicha separación se utilizó el inventario completo. Se han generado un total de 10 modelos diferentes en los cuales se varían los datos de deslizamiento y no deslizamiento para tomar en cuenta la sensibilidad que tiene el modelo a los datos de entrada.
Para cada uno de los modelos se ha tomado el 70% del grupo de datos de deslizamiento de forma aleatoria y se ha tomado a misma cantidad de puntos con la diferencia que corresponderán a datos de no deslizamiento. El 30% de los datos restantes se utilizaron para verificar el desempeño con el que cuenta el modelo para poder predecir deslizamientos que no fueron usados.
Como base se utilizó el script en R que ha desarrollado la Cooperación Italiana en el marco del Proyecto Regional de Formación Aplicada a los Escenarios de Riesgos con Vigilancia y Monitoreo de los Fenómenos Volcánicos, Sísmicos Hidrogeológicos en Centro América - RIESCA. El script utiliza como entrada las capas en formato raster con un tamaño de celda de 12.5 metros. En el caso de los puntos de deslizamientos, la geología y la capa de cobertura forestal fueron convertidos al mismo tamaño de celda y los que son de tipo vectorial fueron convertidos a raster.
Verificación del ajuste del modelo
Con el objetivo de verificar de qué manera se ajusta el modelo a los datos utilizados se emplean los siguientes parámetros:
El análisis del área bajo la curva haciendo uso de las curvas ROC que fueron obtenidas al usar el paquete pRoc de R. Este paquete es una herramienta para visualizar, suavizar y comparar las características de funcionamiento del receptor (curvas ROC). El área (parcial) debajo de la curva (AUC) se puede comparar con pruebas estadísticas basadas en estadísticas U o bootstrap. (Robin, Turck, & Hainard, 2020)
La matriz de confusión que incluye la precisión (accuracy), sensibilidad (sensivity), especificidad (specificity), capacidad de predicción de posyovos (positive prediction value (PPV)) y la capacidad de predicción de negativos (negative prediction value (NPV)). La matriz de confusión se define de la siguiente manera ver tabla 6:
Total= () | Referencia | ||
---|---|---|---|
No deslizamiento (0) | Deslizamiento (1) | ||
Predicción | No deslizamiento (0) | TN | FN |
Deslizamiento (1) | FP | TP |
Fuente: Elaboración propia
Dónde:
TN |
significa True Negatives (Verdaderos Negativos) |
FN |
significa False Negatives (Falsos Negativos) |
FP |
significa False Positives (Falsos Positivos) |
TP |
significa True Positives (Verdaderos Positivos) |
Los indicadores se calculan con las siguientes ecuaciones:
Accuracy=(TP+TN/(TP+FP+TN+FN)
Sensivity=TP/(TP+FN)
Specificity=TP(FP+TN)
PPV=TP/(TP+FP)
NPV=TN/(TN+FN)
El área bajo la curva es un indicador de la capacidad del modelo para poder separar la ocurrencia de deslizamientos (positivos) y de no deslizamientos (negativos), el valor del área bajo la curva debe estar entre 0 y 1. Si se obtiene un valor de 0 significa que el modelo cambia losa datos, es decir que convierte los valores positivos en negativos y viceversa, un valor de 1 indica que el modelo separa los casos completamente y si se obtiene un valor de 0.5 significa que el modelo no tiene la capacidad de separar los casos.
Verificación del desempeño del modelo
Par poder validad el desempeño que tiene el modelo para predecir los deslizamientos, se utilizó el 30% de los datos del inventario que se generó, calculando los mismos indicadores que se utilizaron para la verificación del ajuste del modelo.
Resultados y discusión
Para el caso del análisis de multicolinealidad como se muestra en la tabla 7, realizado mediante el Variance Influence Factor (VIF) muestra que ningún valor es mayor a 10, por lo que se utilizaron todas las variables.
No. | Variable | VIF |
---|---|---|
1 | ele | 4.820756 |
2 | asp | 3.709486 |
3 | plc | 4.837379 |
4 | prc | 4.850403 |
5 | slo | 8.316688 |
6 | lcl | 6.273585 |
7 | twi | 4.697236 |
8 | uso | 1.735107 |
9 | geo | 3.985416 |
Fuente: Elaboración propia
Para el caso de las curvas ROC que se obtuvieron haciendo uso del paquete pROC de R, se crearon dos curvas, una que corresponde a los datos de calibración y otra correspondiente a los datos de validación. Para el caso de los datos de calibración se obtuvo un área bajo la curva de AUC=0.948 y para los datos de validación se obtuvo un UAC= 0.934. Las curvas que se muestran a continuación muestran el promedio de las 10 curvas ROC generadas para cada modelo (ver figura 1 y 2).
La matriz de confusión obtenida para los datos de calibración se muestra en las siguientes tablas 8 y 9:
Predicción | Total= 4,510 | Referencia | |
No deslizamiento (0) | Deslizamiento (1) | ||
No deslizamiento (0) | 3767 | 283 | |
Deslizamiento (1) | 741 | 4277 |
Fuente: Elaboración propia
Los parámetros de la matriz de confusión son:
Accuracy | 0.8863 |
Sensitivity | 0.9373 |
Specificity | 0.8353 |
Pos Pred Value | 0.8505 |
Neg Pred Value | 0.9301 |
Fuente: Elaboración propia
La matriz de confusión para los datos de validación se muestra en la siguiente tabla 10 y 11:
Predicción | Total= 1,930 | Referencia | |
No deslizamiento (0) | Deslizamiento (1) | ||
No deslizamiento (0) | 1672 | 217 | |
Deslizamiento (1) | 258 | 1713 |
Fuente: Elaboración propia
Los parámetros de la matriz de confusión son:
Accuracy | 0.8769 |
Sensitivity | 0.8876 |
Specificity | 0.8663 |
Pos Pred Value | 0.8691 |
Neg Pred Value | 0.8851 |
Fuente: Elaboración propia
Mapa de susceptibilidad a deslizamientos
Para cada pixel del raster de susceptibilidad y para cada set de datos se ha generado un valor de probabilidad entre 0 y 1. El Mapa 11 muestra la susceptibilidad a deslizamientos inducidos por lluvia, cuyo resultado final representa el promedio de los valores de probabilidad de cada uno de los 10 modelos creados. El mapa se ha reclasificado en cinco categorías usando el método de “Natural Braeks” implementado en ArcGIS.
Conclusiones
La metodología MARS permite realizar de manera automática los aspectos de modelación de la regresión clásica, seleccionando las variables predictoras, estimando los valores perdidos, transformando variables, detectando interacciones, contrastando y asegurando la correcta construcción del modelo, y permitiendo resultados más exactos y completos.
El mapa de susceptibilidad inducidos por lluvia generado presenta una buena base como un modelo de predicción de deslizamiento, el cual puede ser mejorar se genera un inventario de deslizamientos más detallado.
Los resultados pueden ser todavía más finos si se incluyen diversas variables, como ser fallas, lugares poblados y demás, con la finalidad que estos mapas sirvan como punto de partido para la correcta gestión de riesgos.