Regresión-kriging

Técnica de predicción espacial

En estadística aplicada y geoestadística , la regresión-kriging ( RK ) es una técnica de predicción espacial que combina una regresión de la variable dependiente sobre variables auxiliares (como parámetros derivados de modelado digital de elevación, teledetección/imágenes y mapas temáticos) con interpolación ( kriging ) de los residuos de regresión. Es matemáticamente equivalente al método de interpolación denominado kriging universal y kriging con deriva externa , donde los predictores auxiliares se utilizan directamente para resolver los pesos de kriging. ^[1]

BLUP para datos espaciales

El modelo universal del esquema de variación espacial.

La regresión kriging es una implementación del mejor predictor lineal insesgado (BLUP) para datos espaciales, es decir, el mejor interpolador lineal que asume el modelo universal de variación espacial. Matheron (1969) propuso que un valor de una variable objetivo en alguna ubicación puede modelarse como una suma de los componentes deterministas y estocásticos: ^[2]

${\displaystyle Z(\mathbf {s} )=m(\mathbf {s} )+\varepsilon '(\mathbf {s} )+\varepsilon ''}$

que denominó modelo universal de variación espacial . Los componentes deterministas y estocásticos de la variación espacial se pueden modelar por separado. Al combinar los dos enfoques, obtenemos:

${\displaystyle {\hat {z}}(\mathbf {s} _{0})={\hat {m}}(\mathbf {s} _{0})+{\hat {e}}(\mathbf {s} _{0})=\sum \limits _{k=0}^{p}{{\hat {\beta }}_{k}\cdot q_{k}(\mathbf {s} _{0})}+\sum \limits _{i=1}^{n}\lambda _{i}\cdot e(\mathbf {s} _{i})}$

donde es la parte determinista ajustada, es el residuo interpolado, son los coeficientes deterministas estimados del modelo ( es la intersección estimada), son los pesos de kriging determinados por la estructura de dependencia espacial del residuo y donde es el residuo en la ubicación . Los coeficientes de regresión se pueden estimar a partir de la muestra mediante algún método de ajuste, por ejemplo, mínimos cuadrados ordinarios (MCO) o, de manera óptima, utilizando mínimos cuadrados generalizados (GLS): ^[3] ${\displaystyle {\hat {m}}(\mathbf {s} _{0})}$ ${\displaystyle {\hat {e}}(\mathbf {s} _{0})}$ ${\displaystyle {\hat {\beta }}_{k}}$ ${\displaystyle {\hat {\beta }}_{0}}$ ${\displaystyle \lambda _{i}}$ ${\displaystyle e(\mathbf {s} _{i})}$ ${\displaystyle {\mathbf {s} }_{i}}$ ${\displaystyle {\hat {\beta }}_{k}}$

${\displaystyle \mathbf {\hat {\beta }} _{\mathtt {GLS}}=\left(\mathbf {q} ^{\mathbf {T} }\cdot \mathbf {C} ^{-\mathbf {1} }\cdot \mathbf {q} \right)^{-\mathbf {1} }\cdot \mathbf {q} ^{\mathbf {T} }\cdot \mathbf {C} ^{-\mathbf {1} }\cdot \mathbf {z} }$

donde es el vector de coeficientes de regresión estimados, es la matriz de covarianza de los residuos, es una matriz de predictores en las ubicaciones de muestreo y es el vector de valores medidos de la variable objetivo. La estimación GLS de coeficientes de regresión es, de hecho, un caso especial de la regresión ponderada geográficamente. En este caso, los pesos se determinan objetivamente para tener en cuenta la autocorrelación espacial entre los residuos. ${\displaystyle \mathbf {\hat {\beta }} _{\mathtt {GLS}}}$ ${\displaystyle \mathbf {C} }$ ${\displaystyle {\mathbf {q} }}$ ${\displaystyle \mathbf {z} }$

Una vez que se ha estimado la parte determinista de la variación (parte de regresión), el residuo se puede interpolar con kriging y agregar a la tendencia estimada. La estimación de los residuos es un proceso iterativo: primero se estima la parte determinista de la variación utilizando MCO, luego se utiliza la función de covarianza de los residuos para obtener los coeficientes GLS. A continuación, estos se utilizan para volver a calcular los residuos, a partir de los cuales se calcula una función de covarianza actualizada, y así sucesivamente. Aunque muchos geoestadísticos recomiendan este procedimiento como el adecuado, Kitanidis (1994) demostró que el uso de la función de covarianza derivada de los residuos de MCO (es decir, una única iteración) suele ser satisfactorio, porque no es lo suficientemente diferente de la función derivada después de varias iteraciones; es decir, no afecta mucho a las predicciones finales. Minasny y McBratney (2007) informan resultados similares: parece que utilizar más datos de mayor calidad es más importante que utilizar métodos estadísticos más sofisticados. ^[4]

En notación matricial, la regresión-kriging se escribe comúnmente como: ^[5]

${\displaystyle {\hat {z}}_{\mathtt {RK}}(\mathbf {s} _{0})=\mathbf {q} _{\mathbf {0} }^{\mathbf {T} }\cdot \mathbf {\hat {\beta }} _{\mathtt {GLS}}+\mathbf {\lambda } _{\mathbf {0} }^{\mathbf {T} }\cdot (\mathbf {z} -\mathbf {q} \cdot \mathbf {\hat {\beta }} _{\mathtt {GLS}})}$

donde es el valor predicho en la ubicación , es el vector de predictores y es el vector de pesos de kriging utilizados para interpolar los residuos. El modelo RK se considera el mejor predictor lineal de datos espaciales . ^{[5] [6]} Tiene una varianza de predicción que refleja la posición de nuevas ubicaciones (extrapolación) tanto en el espacio geográfico como en el de características: ${\displaystyle {\hat {z}}({\mathbf {s} }_{0})}$ ${\displaystyle {\mathbf {s} }_{0}}$ ${\displaystyle {\mathbf {q} }_{\mathbf {0} }}$ ${\displaystyle p+1}$ ${\displaystyle \mathbf {\lambda } _{\mathbf {0} }}$ ${\displaystyle n}$

${\displaystyle {\hat {\sigma }}_{\mathtt {RK}}^{2}(\mathbf {s} _{0})=(C_{0}+C_{1})-\mathbf {c} _{\mathbf {0} }^{\mathbf {T} }\cdot \mathbf {C} ^{-\mathbf {1} }\cdot \mathbf {c} _{\mathbf {0} }+\left(\mathbf {q} _{\mathbf {0} }-\mathbf {q} ^{\mathbf {T} }\cdot \mathbf {C} ^{-\mathbf {1} }\cdot \mathbf {c} _{\mathbf {0} }\right)^{\mathbf {T} }\cdot \left(\mathbf {q} ^{\mathbf {T} }\cdot \mathbf {C} ^{-\mathbf {1} }\cdot \mathbf {q} \right)^{\mathbf {-1} }\cdot \left(\mathbf {q} _{\mathbf {0} }-\mathbf {q} ^{\mathbf {T} }\cdot \mathbf {C} ^{-\mathbf {1} }\cdot \mathbf {c} _{\mathbf {0} }\right)}$

donde es la variación del umbral y es el vector de covarianzas de los residuos en la ubicación no visitada. ${\displaystyle C_{0}+C_{1}}$ ${\displaystyle {\mathbf {c} }_{0}}$

Árbol de decisiones para seleccionar un modelo de predicción espacial adecuado.

Muchos (geo)estadísticos creen que solo hay un modelo de predicción lineal imparcial óptimo para datos espaciales (por ejemplo, regresión-kriging), todas las demás técnicas como el kriging ordinario, la correlación ambiental, el promedio de valores por polígonos o la interpolación de distancia inversa pueden verse como sus casos especiales. Si los residuos no muestran autocorrelación espacial (efecto nugget puro), la regresión-kriging converge a regresión lineal múltiple pura, porque la matriz de covarianza ( ) se convierte en una matriz identidad. Del mismo modo, si la variable objetivo no muestra correlación con los predictores auxiliares, el modelo de regresión-kriging se reduce al modelo de kriging ordinario porque la parte determinista es igual al valor medio (global). Por lo tanto, el kriging puro y la regresión pura deben considerarse solo como casos especiales de regresión-kriging (ver figura). ${\displaystyle \mathbf {C} }$

RK y Reino Unido/KED

La literatura geoestadística utiliza muchos términos diferentes para lo que son esencialmente las mismas técnicas o al menos muy similares. Esto confunde a los usuarios y los distrae de usar la técnica correcta para sus proyectos cartográficos. De hecho, tanto el kriging universal, el kriging con deriva externa y el kriging de regresión son básicamente la misma técnica.

Matheron (1969) denominó originalmente la técnica Le krigeage universel , sin embargo, la técnica fue pensada como un caso generalizado de kriging donde la tendencia se modela como una función de coordenadas. Por lo tanto, muchos autores reservan el término kriging universal (UK) para el caso en el que solo se utilizan las coordenadas como predictores. Si la parte determinista de la variación ( deriva ) se define externamente como una función lineal de algunas variables auxiliares, en lugar de las coordenadas, se prefiere el término kriging con deriva externa (KED) (según Hengl 2007, "Acerca de la regresión-kriging: de las ecuaciones a los estudios de casos"). En el caso de UK o KED, las predicciones se realizan como con kriging, con la diferencia de que la matriz de covarianza de los residuos se extiende con los predictores auxiliares. Sin embargo, la deriva y los residuos también se pueden estimar por separado y luego sumarse. Este procedimiento fue sugerido por Ahmed et al. (1987) y Odeh et al. (1995) más tarde lo denominó regresión-kriging , mientras que Goovaerts (1997) utiliza el término kriging con un modelo de tendencia para referirse a una familia de interpoladores, y se refiere a RK como kriging simple con medias locales variables . Minasny y McBratney (2007) simplemente llaman a esta técnica Predictor empírico lineal imparcial mejor, es decir, E-BLUP . ^{[7] [8] [9] [4]}

En el caso de KED, las predicciones en nuevas ubicaciones se realizan mediante:

${\displaystyle {\hat {z}}_{\mathtt {KED}}(\mathbf {s} _{0})=\sum \limits _{i=1}^{n}w_{i}^{\mathtt {KED}}(\mathbf {s} _{0})\cdot z(\mathbf {s} _{i})}$

para

${\displaystyle \sum \limits _{i=1}^{n}w_{i}^{\mathtt {KED}}(\mathbf {s} _{0})\cdot q_{k}(\mathbf {s} _{i})=q_{k}(\mathbf {s} _{0})}$

para o en notación matricial: ${\displaystyle k=1,\ldots ,p}$

${\displaystyle {\hat {z}}_{\mathtt {KED}}(\mathbf {s} _{0})=\mathbf {\delta } _{\mathbf {0} }^{\mathbf {T} }\cdot \mathbf {z} }$

donde es la variable objetivo, 's son las variables predictoras, es decir, los valores en una nueva ubicación , es el vector de pesos KED ( ), es el número de predictores y es el vector de observaciones en las ubicaciones primarias. Los pesos KED se resuelven utilizando las matrices extendidas: ${\displaystyle z}$ ${\displaystyle q_{k}}$ ${\displaystyle ({\mathbf {s} }_{0})}$ ${\displaystyle {\mathbf {\delta } }_{\mathbf {0} }}$ ${\displaystyle w_{i}^{\mathtt {KED}}}$ ${\displaystyle p}$ ${\displaystyle \mathbf {z} }$ ${\displaystyle n}$

${\displaystyle \mathbf {\lambda } _{\mathbf {0} }^{\mathtt {KED}}=\left\{w_{1}^{\mathtt {KED}}(\mathbf {s} _{0}),\ldots ,w_{n}^{\mathtt {KED}}(\mathbf {s} _{0}),\varphi _{0}(\mathbf {s} _{0}),\ldots ,\varphi _{p}(\mathbf {s} _{0})\right\}^{\mathbf {T} }=\mathbf {C} ^{{\mathtt {KED}}-1}\cdot \mathbf {c} _{\mathbf {0} }^{\mathtt {KED}}}$

donde es el vector de pesos resueltos, son los multiplicadores de Lagrange, es la matriz de covarianza extendida de residuos y es el vector extendido de covarianzas en la nueva ubicación. ${\displaystyle {\mathbf {\lambda } }_{\mathbf {0} }^{\mathtt {KED}}}$ ${\displaystyle \varphi _{p}}$ ${\displaystyle {\mathbf {C} }^{\mathtt {KED}}}$ ${\displaystyle {\mathbf {c} }_{\mathbf {0} }^{\mathtt {KED}}}$

En el caso de KED, la matriz de covarianza extendida de residuos se ve así (Webster y Oliver, 2007; p. 183): ^[10]

${\displaystyle \mathbf {C} ^{\mathtt {KED}}=\left[{\begin{array}{ccccccc}C(\mathbf {s} _{1},\mathbf {s} _{1})&\cdots &C(\mathbf {s} _{1},\mathbf {s} _{n})&1&q_{1}(\mathbf {s} _{1})&\cdots &q_{p}(\mathbf {s} _{1})\\\vdots &&\vdots &\vdots &\vdots &&\vdots \\C(\mathbf {s} _{n},\mathbf {s} _{1})&\cdots &C(\mathbf {s} _{n},\mathbf {s} _{n})&1&q_{1}(\mathbf {s} _{n})&\cdots &q_{p}(\mathbf {s} _{n})\\1&\cdots &1&0&0&\cdots &0\\q_{1}(\mathbf {s} _{1})&\cdots &q_{1}(\mathbf {s} _{n})&0&0&\cdots &0\\\vdots &&\vdots &\vdots &\vdots &&\vdots \\q_{p}(\mathbf {s} _{1})&\cdots &q_{p}(\mathbf {s} _{n})&0&0&\cdots &0\end{array}}\right]}$

y así: ${\displaystyle \mathbf {c} _{\mathbf {0} }^{\mathtt {KED}}}$

${\displaystyle \mathbf {c} _{\mathbf {0} }^{\mathtt {KED}}=\left\{C(\mathbf {s} _{0},\mathbf {s} _{1}),\ldots ,C(\mathbf {s} _{0},\mathbf {s} _{n}),q_{0}(\mathbf {s} _{0}),q_{1}(\mathbf {s} _{0}),\ldots ,q_{p}(\mathbf {s} _{0})\right\}^{\mathbf {T} };q_{0}(\mathbf {s} _{0})=1}$

Por lo tanto, KED se ve exactamente como un kriging ordinario, excepto que la matriz/vector de covarianza se extiende con valores de predictores auxiliares.

Aunque a primera vista la KED parece ser computacionalmente más sencilla que la RK, los parámetros del variograma para la KED también deben estimarse a partir de los residuos de regresión, lo que requiere un paso de modelado de regresión independiente. Esta regresión debería ser GLS debido a la probable correlación espacial entre los residuos. Nótese que muchos analistas utilizan en cambio los residuos MCO, que pueden no ser demasiado diferentes de los residuos GLS. Sin embargo, no son óptimos si hay alguna correlación espacial y, de hecho, pueden ser bastante diferentes para puntos de muestra agrupados o si el número de muestras es relativamente pequeño ( ). ${\displaystyle \ll 200}$

Una limitación de KED es la inestabilidad de la matriz extendida en el caso de que la covariable no varíe suavemente en el espacio. RK tiene la ventaja de que separa explícitamente la estimación de tendencia de la predicción espacial de residuos, lo que permite el uso de formas de regresión arbitrariamente complejas, en lugar de las técnicas lineales simples que se pueden utilizar con KED. Además, permite la interpretación separada de los dos componentes interpolados. El énfasis en la regresión es importante también porque el ajuste de la parte determinista de la variación (regresión) es a menudo más beneficioso para la calidad de los mapas finales que el ajuste de la parte estocástica (residuos).

Software para ejecutar regresión-kriging

Ejemplo de un marco genérico para la predicción espacial de variables del suelo basado en regresión-kriging. ^[9]

La regresión kriging se puede automatizar, por ejemplo, en el entorno de cálculo estadístico R, mediante el uso del paquete gstat o geoR. Las entradas y salidas típicas incluyen:

ENTRADAS:

• Conjunto de interpolación (mapa de puntos) — en ubicaciones primarias; ${\displaystyle z(\mathbf {s} _{i})}$ ${\displaystyle i=1,\ldots ,n}$
• Valores mínimos y máximos esperados y precisión de la medición ( ); ${\displaystyle \Delta z}$
• Predictores continuos (mapa rasterizado) — ; en nuevas ubicaciones no visitadas ${\displaystyle q(\mathbf {s} )}$
• Predictores discretos (mapa de polígonos);
• Conjunto de validación (mapa de puntos) — (opcional); ${\displaystyle z*(\mathbf {s} _{j})}$ ${\displaystyle j=1,\ldots ,l}$
• Espaciado de retardo y distancia límite (necesaria para ajustar el variograma);

SALIDAS:

• Mapa de predicciones y error relativo de predicción;
• Mejor subconjunto de predictores y significancia de correlación (R-cuadrado ajustado);
• Parámetros del modelo de variograma (por ejemplo , , , ) ${\displaystyle C_{0}}$ ${\displaystyle C_{1}}$ ${\displaystyle R}$
• Coeficientes de deriva del modelo GLS;
• Precisión de la predicción en los puntos de validación: error medio de predicción (MPE) y error cuadrático medio de predicción (RMSPE);

Aplicación de la regresión-kriging

El kriging de regresión se utiliza en diversos campos aplicados, desde meteorología, climatología, cartografía de suelos, cartografía geológica, modelado de distribución de especies y similares. El único requisito para utilizar el kriging de regresión en lugar del kriging ordinario, por ejemplo, es que existan una o más capas de covariables y que estén significativamente correlacionadas con la característica de interés. Algunas aplicaciones generales del kriging de regresión son:

• Mapeo geoestadístico: la regresión-kriging permite el uso de técnicas geoestadísticas híbridas para modelar, por ejemplo, la distribución espacial de las propiedades del suelo.
• Reducción de escala de mapas: la regresión kriging se puede utilizar como marco para reducir la escala de varios mapas en cuadrícula existentes. En este caso, las capas de covariables deben estar disponibles con una resolución mejor (que corresponde a la intensidad de muestreo) que los datos de puntos originales. ^[11]
• Propagación de errores : los mapas simulados generados mediante un modelo de regresión kriging se pueden utilizar para probar escenarios y estimar la incertidumbre propagada.

Simulaciones de concentraciones de zinc derivadas mediante un modelo de regresión Kriging. Este modelo utiliza una covariable continua (distancia al río) y una categórica (frecuencia de inundaciones). El código utilizado para generar estos mapas está disponible aquí.

Los algoritmos basados ​​en kriging de regresión desempeñan un papel cada vez más importante en la geoestadística porque el número de posibles covariables aumenta cada día. ^[1] Por ejemplo, los DEM están ahora disponibles en varias fuentes. Ahora se pueden solicitar imágenes detalladas y precisas de la topografía a sistemas de teledetección como SPOT y ASTER ; SPOT5 ofrece el escáner estereoscópico de alta resolución (HRS), que se puede utilizar para producir DEM con resoluciones de hasta 5 m. ^[12] También se pueden obtener diferencias más finas en la elevación con escáneres láser aerotransportados. El costo de los datos es gratuito o disminuye a medida que avanza la tecnología. La NASA registró la mayor parte de la topografía del mundo en la Misión Topográfica Radar del Transbordador en 2000. ^[13] Desde el verano de 2004, estos datos han estado disponibles (por ejemplo, a través del FTP del USGS) para casi todo el globo con una resolución de aproximadamente 90 m (para el continente norteamericano con una resolución de aproximadamente 30 m). Asimismo, las imágenes multiespectrales MODIS están disponibles para su descarga gratuita en resoluciones de 250 m. También se puede descargar un gran repositorio gratuito de imágenes Landsat a través del Global Land Cover Facility (GLCF).

Referencias

1. Pebesma, Edzer J (1 de julio de 2006). "El papel de las variables externas y las bases de datos SIG en el análisis geoestadístico" (PDF) . Transactions in GIS . 10 (4): 615–632. doi :10.1111/j.1467-9671.2006.01015.x. S2CID  22146107.
2. Matheron, Georges (1969). "Parte 1 de Cahiers du Centre de morphologie mathématique de Fontainebleau". El krigeage universal . Escuela Nacional Superior de Minas de París.
3. Cressie, Noel (2012). Estadísticas para datos espacio-temporales . Hoboken, NJ: Wiley. ISBN 978-0-471-69274-4.
4. Minasny, Budiman; McBratney, Alex B. (31 de julio de 2007). "Predicción espacial de las propiedades del suelo utilizando EBLUP con la función de covarianza Matérn". Geoderma . 140 (4): 324–336. Bibcode :2007Geode.140..324M. doi :10.1016/j.geoderma.2007.04.028.
5. Christensen, Ronald (2001). Modelado lineal avanzado: datos multivariados, de series temporales y espaciales; regresión no paramétrica y maximización de la superficie de respuesta (2. ed.). Nueva York, NY [ua]: Springer. ISBN 978-0-387-95296-3.
6. Goldberger, AS (1962). "Mejor predicción lineal imparcial en el modelo de regresión lineal generalizada". Revista de la Asociación Estadounidense de Estadística . 57 (298): 369–375. doi :10.1080/01621459.1962.10480665. JSTOR  2281645.
7. Ahmed, Shakeel; De Marsily, Ghislain (1 de enero de 1987). "Comparación de métodos geoestadísticos para estimar la transmisividad utilizando datos sobre transmisividad y capacidad específica". Investigación de recursos hídricos . 23 (9): 1717. Bibcode :1987WRR....23.1717A. doi :10.1029/WR023i009p01717.
8. Odeh, IOA; McBratney, AB; Chittleborough, DJ (31 de julio de 1995). "Resultados adicionales sobre la predicción de propiedades del suelo a partir de atributos del terreno: cokriging heterotópico y kriging de regresión". Geoderma . 67 (3–4): 215–226. Código Bibliográfico :1995Geode..67..215O. doi :10.1016/0016-7061(95)00007-B.
9. Hengl, Tomislav; Heuvelink, Gerard BM; Stein, Alfred (30 de abril de 2004). "Un marco genérico para la predicción espacial de variables del suelo basado en la regresión-kriging" (PDF) . Geoderma . 120 (1–2): 75–93. Bibcode :2004Geode.120...75H. doi :10.1016/j.geoderma.2003.08.018.
10. Webster, Richard; Oliver, Margaret A. (2007). Geoestadística para científicos ambientales (2.ª ed.). Chichester: Wiley. ISBN 978-0-470-02858-2.
11. Hengl, Tomislav; Bajat, Branislav; Blagojević, Dragan; Reuter, Hannes I. (1 de diciembre de 2008). "Modelado geoestadístico de la topografía mediante mapas auxiliares" (PDF) . Computers & Geosciences . 34 (12): 1886–1899. Bibcode :2008CG.....34.1886H. doi :10.1016/j.cageo.2008.01.005.
12. Toutin, Thierry (30 de abril de 2006). "Generación de DSM a partir de datos estéreo HRS en trayectoria y HRG transversales a trayectoria SPOT-5 mediante espaciotriangulación y autocalibración". ISPRS Journal of Photogrammetry and Remote Sensing . 60 (3): 170–181. Bibcode :2006JPRS...60..170T. doi :10.1016/j.isprsjprs.2006.02.003.
13. Rabus, Bernhard; Eineder, Michael; Roth, Achim; Bamler, Richard (31 de enero de 2003). "La misión de topografía por radar del transbordador: una nueva clase de modelos digitales de elevación adquiridos por radar espacial". ISPRS Journal of Photogrammetry and Remote Sensing . 57 (4): 241–262. Bibcode :2003JPRS...57..241R. doi :10.1016/S0924-2716(02)00124-7.

Lectura adicional

• Capítulo 2, Regresión-kriging , en Tomislav Hengl (2009), A Practical Guide to Geostatistical Mapping , 291 pág., ISBN 978-90-9024981-0 . [1] 
• Hengl T.; Heuvelink GBM; Rossiter DG (2007). "Acerca de la regresión-kriging: de ecuaciones a estudios de casos". Computers & Geosciences . 33 (10): 1301–1315. Bibcode :2007CG.....33.1301H. doi :10.1016/j.cageo.2007.05.001.

Enlaces externos

• Paquete Gstat (implementa KED)
• Paquete GeoR (implementa KED)