Diseño de una red de monitoreo óptima

Resumen

La geoestadística está conformada por un conjunto de técnicas para el análisis y predicción de valores de una propiedad distribuida en el espacio y/o en el tiempo. En contraposición con la estadística clásica o convencional, tales valores no se consideran independientes, sino que se supone de manera implícita que están correlacionados unos con otros, es decir, que existe una dependencia espacial y/o temporal.

El método propuesto para el diseño de la red de monitoreo utiliza una combinación de métodos geoestadísticos, mediante los cuales se obtiene la matriz de covarianza de los datos de carga hidráulica. Para ello se utiliza un filtro de Kalman que selecciona los puntos de muestreo con base en la reducción de la varianza, y en un método de optimización heurístico.

La red piloto inicial quedó integrada por 84 aprovechamientos distribuidos en la zona de estudio. La metodología seguida para el análisis geoestadístico de cada uno de los parámetros consiste de tres pasos que son: análisis exploratorio de los datos, el análisis estructural y las predicciones (kriging o simulaciones).

El objetivo del análisis exploratorio de los datos es caracterizar a la muestra tratando de obtener la mayor información posible a partir de los datos disponibles. El objetivo del análisis estructural es caracterizar la estructura espacial de una variable regionalizada. El modelo teórico que se utilizó en este trabajo fue el esférico con pepita (o nugget). Para la predicción de valores en puntos no muestrales, se eligió la técnica de kriging ordinario.

El diseño de la red de monitoreo piezométrica seleccionada quedó integrada por 14 pozos los cuales optimizan el número mínimo de las posiciones de los pozos que permitirán obtener una buena estimación de los niveles de la carga hidráulica en todo el acuífero de Janos, Chih.

Palabras clave: geoestadística, covarianza, red de monitoreo, optimización, filtro de Kalman

Introducción

El acuífero de Janos se localiza en la porción noroeste del estado de Chihuahua, colindando al Norte con él límite internacional de los Estados Unidos de Norteamérica, al Este con los municipios de Ascensión y Nuevo Casas Grandes, al Sur con el municipio de Casas Grandes y al Oeste con el estado de Sonora (Figura 1).

La información generada en el estudio Actualización de mediciones piezométricas de los acuíferos denominados: Janos, Buenaventura y Tabalaopa-Aldama, Chihuahua (MI, 2004), se tomó como base para el diseño de la red de monitoreo.

El método propuesto para el diseño de la red de monitoreo utiliza una combinación de métodos geoestadísticos, mediante los cuales se obtiene la matriz de covarianza de los datos de carga hidráulica. Para ello se utiliza un filtro de Kalman que selecciona los puntos de muestreo con base en la reducción de la varianza, y en un método de optimización heurístico.

Figura 1 Localización de la zona de estudio

Establecimiento de la red de monitoreo piloto

De la información histórica de mediciones piezométricas en el acuífero Janos (Guysa, 2000; MI, 2004), los dos primeros años de la Tabla 1 son reportados en el estudio realizado por la empresa Guysa, y los subsecuentes en el estudio realizado por la empresa MI. Para este acuífero se contó con muy poca información histórica, y dado que no se tiene un censo actualizado se dispuso de la información preexistente. Se optó por la piezometría del año más reciente para el análisis geoestadístico, tomando estos datos como la red piloto que sirvió como base para el diseño de la red óptima.

Tabla 1 Registro histórico de piezometría

La red piloto inicial quedó integrada por 84 aprovechamientos distribuidos en la zona de estudio. Su distribución espacial se puede observar en la Figura 2. Esta red se reporta que existe en (MI, 2004), su toma de datos no es constante como se aprecia en la Tabla 1.

Figura 2 Red de monitoreo piloto

Los pozos utilizados para el monitoreo del acuífero presentan profundidades que varían en el rango de 2.7 a 86.5 m.

Análisis geoestadístico de los datos

La metodología seguida en este trabajo para el análisis geoestadístico de cada uno de los parámetros consiste de tres pasos que son: análisis exploratorio de los datos, el análisis estructural y las predicciones (kriging o simulaciones).

Análisis exploratorio de los datos

El objetivo del análisis exploratorio de los datos es caracterizar a la muestra tratando de obtener la mayor información posible a partir de los datos disponibles. Dicho análisis se basa en técnicas estadísticas convencionales que permiten obtener información, desconocida a priori sobre la muestra bajo estudio, que es imprescindible para realizar "correctamente" cualquier análisis estadístico y en particular un análisis geoestadístico.

Para validar el análisis geoestadístico, éste debe cumplir los siguientes requisitos: que la muestra no se vea afectada por valores atípicos (outliers), tanto distribucionales como espaciales, que su distribución de probabilidad sea normal, que no exista tendencia y que la población tenga una distribución espacial homogénea (Díaz, 2002).

Detección de outliers

En este paso se deben detectar y posiblemente eliminar los outliers de cualquiera de las siguientes dos clases: los distribucionales que son los que tienen valores que se alejan significativamente del valor medio de la muestra; y los espaciales que son los que tienen un valor muy diferente al de sus vecinos más cercanos.

Prueba de normalidad

Esta prueba consiste en una revisión de la distribución de la probabilidad de los datos procurando que sea normal; de no ser así, se puede proceder eliminando los outliers detectados o bien, haciendo una transformación de los datos hasta que se tenga una distribución normal.

Análisis espacial

El análisis espacial consiste en una revisión de la distribución espacial de la muestra, es decir, de la forma como se encuentra distribuida en el área de estudio.

Estadísticos

En este paso se debe realizar un análisis de la información estadística que arrojan los estadísticos. Los valores de importancia son el valor medio, varianza, desviación estándar, coeficiente de variación, valor mínimo, valor máximo, mediana, etc.

Análisis estructural

El proceso de estimación y modelación de la función que describe la correlación espacial es conocido como "análisis estructural". El objetivo del análisis estructural es caracterizar la estructura espacial de una variable regionalizada. En este trabajo la función que se utilizó para describir la correlación espacial de la variable fue el semivariograma, por éste motivo el análisis estructural consistió en ajustar un modelo al semivariograma muestral.

Ajuste del modelo

Para llevar a cabo el ajuste del modelo se debe buscar una función analítica que represente adecuadamente los valores estimados del semivariograma muestral también llamado semivariograma experimental. Los modelos teóricos que pueden ser usados son: el esférico, el gaussiano, el exponencial, de potencia, lineal con meseta y sin meseta (Samper y Carrera, 1990).

El modelo teórico que se utilizó en este trabajo fue el esférico con pepita (o nugget). Este modelo se caracteriza por tres parámetros: el alcance o radio de correlación espacial, la meseta (sill) y la microvarianza o varianza nugget. La fórmula que define el modelo esférico es la siguiente:

donde es el alcance y la meseta.

Se realizó un ajuste a prueba y error para cada uno de los modelos y se seleccionó el que mostró un mejor ajuste según el criterio de información de Akaike, el cual es un compromiso entre la complejidad del modelo (número de parámetros) y su bondad de ajuste (valor medio cuadrático del error).

El criterio de información de Akaike, se define como:

Predicciones

Una vez realizado el análisis estructural, la predicción de valores en puntos no muestrales se puede hacer aplicando la técnica de interpolación "kriging". Para la predicción de valores en puntos no muestrales, se eligió la técnica de kriging ordinario mediante la cual se estiman los valores del parámetro y la varianza para una malla que cubre el área de estudio.

Kriging ordinario

Se estimaron los diferentes parámetros mediante kriging ordinario porque no existe una tendencia significativa. El kriging empleado en el presente trabajo es el ordinario ya que es el método usado más comúnmente por su robustez.

Validación cruzada

En este caso se siguió el procedimiento de validación cruzada mediante el método de leave one out que permite establecer si el modelo de semivariograma obtenido es representativo de la variabilidad espacial de los datos usados en la estimación del mismo. El método consiste en sacar un elemento de la muestra y estimar el valor en ese punto con kriging ordinario empleando el modelo de semivariograma obtenido.

Como resultado se tiene un mapa de las diferencias entre el valor real y el estimado.

Se realiza un análisis integral de los estadísticos de las diferencias y se aplican de manera combinada algunos de los siguientes criterios: i) el valor medio de las diferencias debe ser cercano a cero, ii) la varianza de las diferencias debe ser pequeña, iii) la varianza normalizada de las diferencias debe ser próxima a la unidad y iv) la correlación muestral entre el valor real y la diferencia dividida por la desviación típica debe ser próxima a cero.

Resultados del análisis geoestadístico de la carga hidráulica

En esta sección se presenta el análisis geoestadístico de la carga hidráulica para obtener la red de monitoreo óptima en el acuífero de Janos, Chihuahua. Los datos analizados fueron los niveles estáticos del agua subterránea del año 2004.

Al iniciar el análisis exploratorio de la carga hidráulica se tiene que para este parámetro se cuenta con 84 puntos de medición de los cuales seis se detectaron como valores atípicos (outliers distribucionales). En la Figura 3 se aprecian con el símbolo de estrella. La carga hidráulica se midió en msnm.

Figura 3 Detección de outliers distribucionales

La prueba de normalidad se llevó a cabo analizando el histograma de la muestra en el que se aprecia una asimetría positiva en los datos (Figura 4a), por lo tanto la distribución no es normal. Para resolver el problema de asimetría en la distribución se hizo una transformación logarítmica a los datos y se le quitaron seis outliers (Figura 4b). Se aprecia que el retiro de estos mejoró la simetría de la distribución. El cálculo subsecuente se realizó con estos datos.

Figura 4 Histograma de los datos de carga hidráulica, 2004

Al semivariograma experimental de carga hidráulica, se le ajustó un modelo de tipo esférico. En la Figura 5 se aprecia gráficamente que el modelo se ajusta de manera aceptable. La línea continua es el semivariograma experimental y la línea indicando el número de pares es el teórico. En la Tabla 3 se muestran los parámetros del modelo obtenido.

Figura 5 Semivariograma experimental y teórico

Tabla 3 Modelación de la función de correlación espacial de la carga hidráulica

Como resultado de la validación cruzada de la carga hidráulica sin outliers, en la Figura 6 aparecen las diferencias entre el valor real y el estimado al quitar cada uno de los datos. Los errores positivos aparecen con el símbolo "+" y los negativos con el símbolo "-". Se observa una distribución uniforme de los mismos en el área de estudio. La Tabla 4 presenta un resumen de los resultados de la validación cruzada. En dicha tabla se observa que el valor medio de los errores en la estimación en los pozos es pequeño.

Figura 6 Distribución de los errores

Tabla 4 Resumen de resultados de la validación cruzada sin seis outliers

Ya teniendo el modelo ajustado, se procedió a realizar la estimación espacial mediante kriging ordinario, el cual consiste en interpolar espacialmente el valor de la propiedad en un punto mediante una combinación lineal de sus vecinos más cercanos. Los resultados del kriging para la carga hidráulica se muestran en la Figura 7.

Figura 7 Estimación de la carga hidráulica

En la Tabla 5 se presenta un resumen de los resultados de las estimaciones. Los valores estimados tanto máximos como mínimos son muy semejantes a los datos tomados en campo. Lo anterior junto con los valores pequeños de la varianza nos indica que se tiene una buena modelación de la distribución de la carga hidráulica en el área de estudio. Los valores del coeficiente de variación y el de asimetría nos indican que la distribución de los valores estimados es simétrica.

Tabla 5 Resumen de resultados de la estimación

Diseño óptimo de la red de monitoreo

Los sitios de monitoreo tomados en cuenta para seleccionar la red de monitoreo óptima son los 78 que conforman la red de monitoreo piloto ya que se consideran aptos para muestreo del agua subterránea.

Método para elegir la red de monitoreo óptima

Para seleccionar la red de monitoreo óptima, se utilizó una combinación de métodos geoestadísticos, un filtro de Kalman y un método de optimización heurístico. El método utilizado en el diseño de la red óptima es una modificación al propuesto por Herrera en 1998 (Júnez, 2005).

El filtro de Kalman puede calcular la varianza del error de la estimación con base en la posición de los sitios de muestreo, sin necesidad de conocer el dato medido en campo. De esta manera el filtro de Kalman se puede utilizar para determinar, dados sitios de muestreo, cómo afecta añadir uno adicional a la varianza del error de la estimación resultante. Para construir la red de monitoreo óptima se minimizó la varianza del error de la estimación en todos los pozos de la red piloto, esto es, cada nuevo punto que se añadió a la red fue aquel que resultara en un valor de la varianza menor. La selección de los puntos se llevó a cabo a través de un método heurístico que trabaja secuencialmente, en cada paso se hace una evaluación del valor de la varianza resultante de añadir cada uno de los pozos seleccionados para muestreo de la red de monitoreo piloto y se elige el que da el valor mínimo de la varianza.

El filtro de Kalman utiliza como dato de entrada la matriz de covarianza del error de la estimación a priori o inicial. En este trabajo se calculó esta matriz de covarianza de la carga hidráulica a partir del semivariograma estimado en el análisis geoestadístico.

La fórmula que se utilizó para hacer dicho cálculo es

El algoritmo elige entonces aquellos pozos que al añadirse minimizan la varianza total del error de la estimación para la carga hidráulica.

En este procedimiento es necesario determinar el número de puntos que se incluirán en la red de monitoreo óptima. Para hacer esto se analizó la varianza total del error de la estimación de todos los parámetros que se obtiene conforme se van añadiendo los pozos seleccionados para formar parte de la red. En la Figura 8 se muestra la varianza total de la estimación contra el número de pozos elegidos para la red de monitoreo. Como se observa, los primeros pozos elegidos proporcionan mucha información y la varianza total se reduce significativamente, y conforme se van añadiendo pozos a la red de monitoreo, la información obtenida por cada pozo va disminuyendo. Para medir la contribución relativa de cada pozo se calculó la diferencia de la varianza total producida con pozos y la varianza total producida por pozos, dividida entre esta última, a lo que llamamos reducción relativa de la varianza total:

En la Figura 9 se muestra la gráfica de la reducción relativa de la varianza total. En esta gráfica se determinó el punto en el que se presenta el valor mínimo de la función como el número posible de pozos a incluir en la red de monitoreo. Como se aprecia en la figura pueden ser 14 y 25 los pozos que conformen la red de monitoreo debido a que allí se presentan los valles de la función analizada.

Figura 8 Varianza total vs. Núm. de pozos

Figura 9 Reducción relativa de la varianza total vs. Núm. de pozos

Para decidir cuál sería el número de pozos que conformaría la red de monitoreo se compararon los resultados obtenidos en las estimaciones con kriging ordinario con 78 datos, y las estimaciones con 14 y 25 respectivamente (Tabla 6).

Para calcular la desviación estándar de las estimaciones se empleó la siguiente fórmula:

Tabla 6 Comparación de resultados de las opciones de la red óptima

Del análisis anterior se concluye que, el utilizar una red de 14 pozos es suficiente para tener una buena estimación en el acuífero.

En la Tabla 7 se presenta la lista de los pozos que conforman la red de monitoreo óptima y en la Figura 10 se muestra la distribución espacial.

Tabla 7 Relación de pozos que conforman la red de monitoreo óptima

Los resultados del kriging para la carga hidráulica de la red se muestran en la Figura 11.

Figura 10 Red de monitoreo óptima

Figura 11 Estimación de la carga hidráulica utilizando la red de monitoreo óptima

Conclusiones

En el censo de aprovechamientos se identificaron 403, 405, 433, 822, 971, 986 aprovechamientos para los años 1971, 1972, 1979, 1988, 1999, 2000, pero en el 2002 solamente se sondearon 89 pozos de una red de monitoreo ya existente. Para el 2004 esta red quedó integrada por 84 pozos solamente, los cuales están distribuidos a lo largo de todo el acuífero.

El diseño de la red de monitoreo piezométrica seleccionada quedó integrada por 14 pozos, los cuales optimizan el número mínimo de las posiciones de los pozos, que permitirán obtener una buena estimación de los niveles de la carga hidráulica en todo el acuífero de Janos, Chih.

Referencias

Díaz, M. A., Geoestadística Aplicada, Instituto de Geofísica, UNAM, Instituto de Geofísica y Astronomía, CITMA, Cuba, 2002, 88 p.

Guysa, Simulación hidrodinámica en el acuífero de Janos, Chihuahua, Geofísica de Exploraciones Guysa, S.A. de C.V., México, 2000, 200 pp.

Herrera, G., Cost Effective Groundwater Quality Sampling Network Design, Ph. D. Dissertation, University of Vermont, Burlington, Vermont, 1998.

Júnez, H. E., Diseño de una red de monitoreo de la calidad del agua para el acuífero Irapuato-Valle, Guanajuato, Tesis de Maestría, Universidad Nacional Autónoma de México, Jiutepec, Morelos, México, 2005.

MI, Actualización piezométrica de los acuíferos denominados: Janos, Buenaventura y Tabaloapa-Aldama, Chihuahua, MORO INGENIERÍA, S. C. para la Gerencia Regional Río Bravo, México, 2004, 40 pp.

Samper, F. J. y J. Carrera, Geoestadística, aplicaciones a la hidrogeología subterránea, Centro Internacional de Métodos Numéricos en Ingeniería, Barcelona, 1990, 447 p.

Autor:

M.I. Pablo Gallardo Almanza,

M.I. Leticia Becerra Soriano