INTRODUCCIÓN
El piedemonte llanero colombiano viene siendo sometido a una agriculturización intensiva (Romero, Flantua, Tansey, y Berrio, 2012); consecuentemente, la sostenibilidad de los servicios ecosistémicos está amenazada (Lavelle et al., 2014). Para mitigar este impacto, es imprescindible contar con técnicas que permitan determinar de manera rápida, precisa y económica, la variabilidad espacial de propiedades del suelo (Hempel et al., 2008).
Técnicas geoestadísticas y locales han sido reconocidas como técnicas convencionales de interpolación espacial para propiedades del suelo (Webster y Oliver, 2007); sin embargo, avances en teledetección han generado información auxiliar de factores formadores del suelo, que podrían aportar a la predicción espacial de propiedades del mismo. A partir de estos avances, se han desarrollado alternativas de técnicas de interpolación “mixtas”, en las cuales combinan kriging e información auxiliar; numerosos trabajos han reportado que las técnicas de interpolación mixtas pueden tener alta precisión de predicción (Bishop y McBratney, 2001; Goovaerts, 1999; Oliver, 2010). Sin embargo, no se tiene evidencia de trabajos que hayan evaluado la eficiencia de predicción entre técnicas geoestadísticas, locales y mixtas, para las condiciones del piedemonte llanero colombiano.
El objetivo de este trabajo fue comparar técnicas de interpolación espacial geoestadísticas, locales y mixtas de propiedades del suelo en un área experimental localizada en Villavicencio, correspondiente al piedemonte llanero colombiano, como se aprecia en la Figura 1.
METODOLOGÍA
Área experimental
Este trabajo se llevó a cabo en un área agrícola de 9.38 has, localizada en Villavicencio, Meta (lat: 4.072; long: -73.580; Figura 1). Los suelos del área de estudio provienen de aluviones y sedimentos aluviales heterométricos recientes, desde finos hasta gruesos y se encuentran en las distintas formas del relieve, la capacidad de retención de agua se considera de media a alta (Rippstein, Escobar y Motta, 2001); el área experimental se divide territorialmente en zonas de acuerdo a los siguientes sistemas productivos: pasturas, plátano, cultivos semestrales (maíz, soya y/o arroz), frutales y sistemas agroforestales.
Información auxiliar del suelo
Un modelo digital de elevación (MDE) fue utilizado como información auxiliar del suelo, a partir de este, índices de terreno fueron calculados utilizando el procedimiento de análisis de terreno de SAGA-GIS v2.0 (Olaya y Conrad, 2009); la Tabla 1 muestra una breve descripción de cada índice (Hengl y Reuter, 2008). El MDE se obtuvo mediante la interpolación de curvas de nivel de acuerdo a la metodología ANUDEM descrita por Hutchinson (1989). La medición de curvas de nivel se realizó con una estación total SOKKIA CX-105 (SOKKIA, Overland Park, KS). La interpolación de las curvas de nivel se realizó con la función “Topo to Raster” del módulo de Spatial Analyst de ArcGis v10.2, el MDE elaborado tuvo una resolución espacial de 10 m.
Índices † | Descripción | |
---|---|---|
Elevación (m) | Elevación en metros. | |
Flujo Acumulado (Flujo_Acum) | Calcula el valor de flujo acumulado en cada pixel, a partir de los valores de flujo acumulado aledaños. Muestra cómo se genera una red de flujos acumulados, determinando los patrones de variabilidad espacial de la acumulación del agua en un área específica. | |
Índice Topográfico de Humedad (TWI) | Los valores de sus celdas describen que tan propenso es un sitio a estar saturado. Tiene en cuenta tanto la geometría de la pendiente local como la pendiente y el área de captación. | |
Índice de Convergencia (Ind_Converg) | Calcula un índice de la convergencia o divergencia de un flujo dentro de un paisaje. Sus resultados se dan en porcentajes. Valores negativos corresponden a condiciones de flujo convergentes y viceversa. | |
Factor de Perdida de Suelo (LS_Factor) | Representa el efecto de la geomorfometría de una zona sobre la perdida de suelos. Combina los efectos de la longitud de la pendiente con su gradiente. | |
Profundidad de valle (Prof_valle) | Calcula la distancia vertical de un punto a un nivel base de una red de canales o flujos de agua. |
†. Índices geomorfométricos e hidrológicos calculados a partir del modelo digital de elevación
Fuente: elaboración propia
Esquema de muestreo de suelos y análisis de laboratorio
El algoritmo hipercubo latino condicionado (HCLc) fue utilizado para determinar la localización óptima de los puntos de muestreo a partir de la elevación; HCLc ha sido utilizado en la generación de cartografía digital de suelos a escala de lote (Castro, Costa, Peralta, y Aparicio, 2015; Schmidt et al., 2014). Los propósitos de utilizar HCLc como esquema de muestreo fueron: (a) evitar el uso de decisiones subjetivas en la selección de localización de muestras, (b) darle un enfoque estadístico al esquema de muestreo de suelos. Para este trabajo se decidió trabajar con 50 muestras de suelo, como cantidad suficiente para evaluar el desempeño predictivo de las técnicas de interpolación espacial (Kerry y Oliver, 2003).
Cada muestra de suelo se completó a partir de tres submuestras de 0-30 cm de profundidad, todas las muestras fueron analizadas en el laboratorio de suelos de la Universidad de los Llanos para: contenido de materia orgánica (MO; Walkley y Black, 1934), fósforo disponible (; Bray y Kurtz, 1945), pH (1:1) con potenciómetro, aluminio intercambiable (Al) por el método KCL, calcio (Ca), magnesio (Mg), potasio (K), sodio (Na; extracción con acetato de amonio 1N y neutro), cobre (Cu), hierro (Fe) y manganeso (Mn; extracción con DTPA).
Selección de índices de terreno
Ckg y REML-EBLUP son métodos de interpolación espacial mixtos que requieren de variables auxiliares para optimizar la predicción espacial; debido a esto, se aplicó un método de selección de índices con mayor relación a cada propiedad del suelo (Guyon y Elisseeff, 2003).
Para la selección de índices, se implementó el algoritmo de bosques aleatorios o random forest (RF; Breiman, 2001), utilizando la librería “RandomForest” del software R v3.2.3. (R Development Core Team, 2015). RF calibra los parámetros “ntree” y “mtry” utilizando como método de estimación del error el “bootstrap” (Breiman, 2001). Este método consiste en tomar una cantidad establecida de muestras aleatorias (bootstrap samples) de un conjunto de datos de calibración; posteriormente, RF ajusta un árbol para cada bootstrap sample. Por cada bootstrap sample Xi, un tercio de los datos de calibración son dejados fuera de la muestra (out-of-bag; OOB), la tasa de error OOB (ERROOB) es estimada a medida que se agregan árboles al bosque. El “mtry” y “ntree” que generen el menor ERROOB deben ser elegidos para el modelo (Breiman, 2001).
La clasificación de importancia de índices se realizó mediante el cálculo de un índice de importancia, el cual se basa en el incremento del cuadrado medio del error (%incMSE), para cada árbol de regresión dentro de cada “bosque”, cuando cada índice fue permutado aleatoriamente en la OOB. Mayor detalle del proceso de selección de variables con RF se encuentra en Behrens, Zhu, Schmidt, y Scholten, (2010).
Técnicas de interpolación espacial
Las técnicas de interpolación espacial seleccionadas fueron: distancia inversa ponderada (1) y Spline (2) como técnicas locales, kriging ordinario (3) y kriging universal (4) como técnicas geoestadísticas, y cokriging (5) y mejor predicción lineal insesgada empírica con máxima verosimilitud restringida (6) como técnicas mixtas. Esta selección se basó en el tipo de dato, la complejidad computacional de cada técnica y los resultados de trabajos previos (Minasny y McBratney, 2007; Peña, Rubiano, Peña, y Chaves, 2009). Las bibliotecas “automap”, “gstat”, “sp” y “fields” del software R v3.2.3 fueron utilizadas para ejecutar todas las técnicas de interpolación. En la Tabla 2 se muestra una breve descripción de las técnicas de interpolación seleccionadas.
Evaluación de predicción
En orden de comparar la efectividad de las técnicas de interpolación, se utilizaron como métricas de error de predicción: el error medio absoluto (MAE) según la ecuación (7), la mediana del error predictivo (MdPE), mediante la ecuación (8) y la raíz del cuadrado medio del error (RMSE) por medio de la ecuación (9).
MAE determina la precisión general de la interpolación espacial, mientras RMSE es una métrica de error muy común para medir precisión de predicción, que da mayor peso a los valores extremos de los resultados; por su parte, MdPE fue utilizado en conjunto con MAE, como una métrica de error alternativa de precisión general de interpolación. En general, MdPE puede reducir el impacto que tienen los valores atípicos sobre la MAE; adicionalmente, se evaluó la diferencia entre valores interpolados y medidos en cada punto de muestreo, para esto, se calculó el error estadístico (SE) usando la ecuación (10).
RESULTADOS
Los resultados sugieren que HCLc capturó significativamente la variabilidad espacial de la elevación utilizando 50 muestras distribuidas como se muestra en la Figura 2; por lo tanto, en este trabajo se corrobora el supuesto de que HCLc puede capturar eficientemente la estructura de distribución espacial de covariables relacionadas con factores formadores del suelo y, por ende, es una alternativa para mejorar los esquemas de muestreo de suelos en condiciones del piedemonte llanero colombiano.
La Tabla 3 muestra los parámetros de estadística descriptiva para las propiedades del suelo, los valores mínimos y máximos de pH y MO eran esperables debido a la génesis de los suelos de la zona y a las condiciones de intemperancia a las cuales se encuentra sometidos.
†. Coeficientes de asimetría y curtosis son adimensionales; ††. N: Distribución normal; L: Distribución Lognormal; ~N: Distribución con Tendencia Normal
Fuente: elaboración propia
Excepto pH y MO, el coeficiente de variación (CV; Giraldü, 2014) fue mayor al 23% para todas las propiedades del suelo; estos resultados fueron similares a los reportados por Orozco, Flores, y Sanabria. (2015), quienes plantearon que las propiedades químicas de los suelos del piedemonte llanero colombiano suelen tener alta heterogeneidad. 𝑃 𝑑𝑖𝑠𝑝 , Mg y Ca, en ese orden, presentaron los mayores CV, estos resultados se esperaban debido a las diferencias de uso y manejo del suelo, la aplicación de fertilizantes fosfatados, la aplicación de cal dolomita y la diferencia de tasa de absorción de P, Ca y Mg que tienen los cultivos implantados en el área experimental. Excepto MO, todos los CV de las propiedades de suelo analizadas fueron similares a los reportados por (Peña et al., 2009), ratificando los efectos de las prácticas de uso y manejo de suelo de los sistemas agroproductivos convencionales sobre las propiedades del suelo en el piedemonte llanero colombiano.
Las propiedades pH, Ca y K presentaron distribución normal, mientras que Al y Mg presentaron distribución con tendencia normal; las demás presentaron distribución lognormal (Tabla 3). Los coeficientes de curtosis calculados para todas las propiedades del suelo fueron diferentes a los reportados por Peña et al., (2009); a partir de este resultado, se plantea que aunque los CV puedan ser similares, la distribución de cada propiedad del suelo puede llegar a tener patrones espaciales característicos para cada área experimental.
Importancia de los índices de terreno
La Figura 3 muestra la clasificación de importancia de los índices de terreno por propiedad del suelo, a partir de la aplicación de RF. Los índices más importantes para cada propiedad fueron utilizados para Ckg y REML-EBLUP.
En la mayoría de las propiedades de suelo, el índice profundidad de valle (Prof_valle) e índice de convergencia (Ind_converg) fueron los de mayor y menor importancia respectivamente; el índice factor de pérdida de suelo (LS Factor) fue el segundo más importante para todas las variables, por otro lado, el índice topográfico de humedad (TWI) solo fue importante para Na. Estos índices, podrían ser los que más contribuyen a la explicación de la varianza y patrones espaciales de cada propiedad del suelo y, por ende, podrían ser utilizados dentro de funciones de pedotransferencia o para determinar patrones espaciales de indicadores de calidad del suelo (Castro et al., 2015).
Precisión de interpolación espacial
La Tabla 4 muestra los MAE, MdPE y RMSE de las predicciones espaciales por propiedad del suelo y técnica de interpolación.
En general, los valores de MAE fueron moderados en todas las técnicas de interpolación, para todas las propiedades del suelo; por su parte, los valores de MdPE tuvieron mayor variabilidad, para la mayoría de las bases intercambiables. Para las interpolaciones de 𝑃 𝑑𝑖𝑠𝑝 y Fe, los valores de RMSE fueron considerablemente más altos y consistentes en magnitud. Las diferencias de valores de RMSE, MAE y MdPE entre propiedades, se deben al efecto del uso y manejo del suelo sobre cada propiedad del suelo, dentro del área experimental (Orozco et al., 2015).
Los valores mínimos de MAE y RMSE de la mayoría de las propiedades del suelo, se presentaron con Ckg y REML-EBLUP, al respecto, Minasny y McBratney (2007) reportaron que REML-EBLUP fue eficiente para interpolar propiedades del suelo, cuando la tendencia espacial es fuerte y el número de muestras es bajo (<200); asimismo, concluyeron que para mejorar la predicción espacial con REML-EBLUP, es prioritario hacer énfasis en las metodologías de muestreo y análisis de laboratorio de las muestras. Por su parte, Ckg, utilizando elevación como información auxiliar, ha sido escasamente reportado como una técnica de interpolación espacial eficiente de propiedades del suelo (Song, Zhang, y Wang., 2014); los resultados sugieren que técnicas de interpolación mixtas que involucren índices de terreno como información auxiliar, podrían tener mayor capacidad predictiva que las técnicas de interpolación geoestadísticas y locales, en condiciones de suelo del piedemonte llanero colombiano.
La Figura 4, la Figura 5 y la Figura 6 muestran las diferencias entre valores predichos y observados por propiedad del suelo y técnica de interpolación; en las técnicas geoestadísticas y mixtas, la mayoría de las figuras de valores predichos observados tienen tendencia de homocedasticidad. Este resultado era esperable, ya que ha sido ampliamente discutido que IDW y Spline suelen marcar tendencia de patrones espaciales, mas no ser muy precisos en la interpolación de valores puntuales (Li y Heap, 2011). En la predicción espacial de Mg, K y Al, únicamente Ckg y EBLUP-RMSE tuvieron mayor capacidad predictiva en comparación con las demás técnicas de interpolación. En las demás propiedades, no se presentó mucha diferencia en capacidad predictiva entre las técnicas geoestadísticas y las mixtas.
La Figura 7 muestra el error estadístico (SE) entre valores interpolados y medidos en cada punto de muestreo para todas las técnicas de interpolación.
Excepto para MO, Fe y Pdisp, la SE fue alta para todas las técnicas de interpolación; por otro lado, excepto para Al, la mayoría de técnicas de interpolación presentaron tendencia a sobreestimación de valores en cada punto. Como se esperaba, las técnicas de interpolación local presentaron tendencia a predecir los valores más extremos para la mayoría de propiedades del suelo, en comparación con las demás técnicas de interpolación. En MO, entre los puntos cinco y veinte se presentó tendencia leve a la sobreestimación, mientras que en los puntos 30 al 40 a la subestimación.
Como se observa en la Figura 2, de los puntos cinco al veinte tienden a estar en zonas con cultivos transitorios, mientras que los puntos 30 al 40 en zonas con cultivos perennes o sistemas agroforestales. Estos resultados sugieren que en zonas con uso y manejo más intensivo del suelo, puede haber una tendencia a tener procesos degradativos más intensos de la MO y, por ende, esto puede afectar los patrones espaciales de la misma, (Peña et al., 2009).
En Pdisp, entre los puntos 40 y 50 se presentó tendencia leve a la sobreestimación, en comparación con los demás puntos de muestreo; de los puntos 40 a 50 se concentran en zonas con cultivos transitorios, los cuales vienen siendo fertilizados de manera continua y espacialmente desordenada. Esto, sin duda, es un cofactor a tener en cuenta para evaluar el desempeño predictivo de las técnicas de interpolación. Por último, Fe no tuvo un patrón claro de sobre o subestimación en la predicción de valores, para todas las técnicas de interpolación, lo cual era esperable debido a las fuertes condiciones de intemperancia en las que se encuentran sometidos los suelos del área experimental.
Comparación de patrones espaciales
La Figura 8 y la Figura 9 muestran las predicciones espaciales de cada propiedad; en general, la mayoría de estas presentaron diferencias entre técnicas de interpolación. Un análisis visual determina que la predicción espacial de MO y de la mayoría de bases intercambiables presentaron mayor heterogeneidad espacial con Ckg y REML-EBLU, coincidiendo parcialmente con lo reportado por Minasny y McBratney (2007). Por su parte Al, Mg y Cu, presentó mayor heterogeneidad espacial con KU, coincidiendo parcialmente con lo reportado por Bishop y McBratney (2001). Cu, presentaron similar patrón espacial con todas las técnicas de interpolación, mientras que Ca presentó mayor heterogeneidad espacial con KO y REML-EBLUP.
Las interpolaciones con IDW y Spline tuvieron mayor homogeneidad espacial para la mayoría de propiedades del suelo, mientras que las interpolaciones con KU presentaron patrones espaciales más heterogéneos que con KO. En general, los resultados sugieren que la inclusión de índices de terreno como variables auxiliares en técnicas de interpolación mixta, puede aportar a la caracterización de los patrones espaciales de todas las propiedades del suelo en condiciones de la altillanura colombiana, excepto pH y Cu.
Los resultados permiten establecer un punto de referencia del uso de técnicas de interpolación espacial en condiciones del piedemonte llanero colombiano; no obstante, estos resultados deben ser validados para otras áreas experimentales aledañas, ya que responden exclusivamente a condiciones climatológicas, de uso y manejo del suelo y edafológicas del área experimental.
CONCLUSIONES
En este trabajo, seis técnicas de interpolación de tipo local (IDW y Spline), geoestadística (KO y KU) y mixta (Ckg y REML-EBLUP) fueron comparadas para predecir propiedades del suelo en condiciones del piedemonte llanero colombiano; los resultados validan la utilidad de HCLc como esquema de muestreo y soportan el supuesto que HCLc captura adecuadamente la distribución de información auxiliar del suelo. Por su parte, RF demostró ser adecuado para determinar los índices de terreno con mayor capacidad predictiva para cada una de las propiedades del suelo.
En general, los resultados sugieren que las técnicas de interpolación mixtas teniendo como información auxiliar del suelo e índices de terreno, proporcionaron una mejora significativa de la predicción de propiedades del suelo, en comparación con las demás técnicas; este resultado es importante para implementar estrategias de mejoramiento de la sostenibilidad de los servicios ecosistémicos del suelo, en condiciones del piedemonte llanero colombiano. Los resultados obtenidos requieren ser validados para otras zonas y con otras fuentes de información auxiliar.