Aplicación del método Persistent Scatterer Interferometry (PSI) en la ciudad de Limón, Costa Rica

Application of Persistent Scatterer Interferometry method (PSI) in Limón, Costa Rica

Diana Paniagua-Jiménez1

Universidad Nacional, Costa Rica, Costa Rica

José Valverde-Calderón2

Universidad Nacional, Costa Rica, Costa Rica

Paula Molina-Calderón3

Universidad Nacional, Costa Rica, Costa Rica

Gustavo Barrantes-Castillo4

Universidad Nacional, Costa Rica, Costa Rica

Introducción

La teledetección es el proceso de detectar las características físicas de un área midiendo y analizando el comportamiento de la radiación emitida desde una plataforma como un satélite o un avión y reflejada por los elementos en la superficie terrestre. Como resultado, se obtiene imágenes que ayudan a los investigadores a monitorear objetos de interés como incendios forestales, actividad volcánica, tormentas de polvo, crecimiento de una ciudad y cambios en tierras de cultivo o bosques, durante varios años o décadas, entre otras aplicaciones .

El Radar de Apertura Sintética (Synthetic Aperture Radar, SAR por sus siglas en Inglés) es una herramienta de observación de la Tierra que utiliza frecuencias de microondas en el rango de longitudes de onda de 2 cm hasta 1 m. Los sistemas SAR se operan desde aviones o plataformas satelitales en órbita (Cho et al., 2006) y proporcionan imágenes de terreno. Una de las ventajas de las imágenes SAR es que almacenan la información de la fase, dato que puede ser luego usado para aplicaciones como cuantificar la deformación de la superficie (Homer et al., 1996).

Algunas de las aplicaciones de imágenes SAR con gran auge en la actualidad son:

•Mapeo de la cobertura forestal: Se ha usado imágenes para el monitoreo de la cobertura boscosa de grandes extensiones de bosques de difícil acceso en el terreno. Los sensores ópticos solo pueden detectar la cubierta superior, donde se produce la absorción y reflexión de partes del espectro de luz visible e infrarroja mientras que los sensores de radar se han utilizado para el mapeo de bosques. Debido a la longitud de onda más larga de un sensor de radar, la radiación penetra en la capa superior de la vegetación hasta cierto punto y se dispersa por tallos, ramas, hojas o agujas (Balzter, 2001).

•Mapeo de cuerpos de agua: SAR se puede utilizar para obtener información sobre las tasas de arrastre del suelo o para general mapas de inundaciones, aun en condiciones meteorológicas desfavorables. Los patrones espaciales de erosión por inundación y sedimentación se pueden mapear a través de la interferometría y, en algunos casos, se pueden determinar sus volúmenes. SAR es utilizado por los hidrólogos como una herramienta para monitorear el agua superficial debido a su capacidad para obtener imágenes de alta resolución en condiciones climáticas adversas y de poca iluminación natural.

Esta característica es útil, por ejemplo, cuando se observan inundaciones durante los períodos de nubosidad constante o se acorta la luz del día. Este enfoque se ha utilizado con éxito para obtener mapas precisos de las inundaciones de los ríos en Siberia y el sur de Francia (Smith, 2002).

En el caso de Costa Rica, tras la revisión bibliográfica, se determinó que hay pocas publicaciones o información acerca del uso de imágenes de radar en el país. Una aplicación es el caso del mapeo de cobertura forestal en la Península de Osa, donde la NASA realizó un inventario forestal en el Parque Nacional Corcovado y en el manglar de Bahía Drake, utilizando el sensor EcoSAR (NASA, 2014).

Otra aplicación fue la realizada también por la NASA, en colaboración con el Programa de Investigaciones Aerotransportadas y Sensores Remotos (PRIAS) para el mapeo de zonas volcánicas en Costa Rica usando el sensor UAVSAR, donde se mapeo los volcanes Turrialba, Irazú, Poas y Arenal. Más información sobre este proyecto se puede consultar en el sitio https://uavsar.jpl.nasa.gov/

La técnica Radar de Apertura Sintética Interferométrica (InSAR) se basa en la comparación de las diferencias de fase entre dos o más imágenes SAR separadas ya sea por tiempo o espacio. InSAR tiene aplicaciones prácticas que van desde medir la deformación de la superficie hasta monitorear procesos geomorfológicos. Usando InSAR, por ejemplo, los científicos han podido medir con precisión los efectos posteriores de los terremotos, mapear el flujo de glaciales y corrientes oceánicas, monitoreo de sitios mineros y deslizamientos de tierra, y tiene un valor potencial en cualquier área donde existan riesgos de desplazamiento natural (Uys, 2016).

En los últimos 20 años, InSAR han demostrado una gran capacidad para mapear la deformación de la superficie de la Tierra inducida por varios mecanismos geofísicos (Samieie-Esfahany et al., 2009). Ejemplo de esto es el trabajo de Delgado-Blasco et al., (2019) quienes utilizaron imágenes de radar para la determinación de la subsidencia en el área urbana de Roma, Italia o el trabajo de (Yhokha et al., 2018) quienes utilizaron InSAR para el monitoreo de deslizamientos en la zona de Nainital, India.

La Interferometría de Dispersión Persistente (PSI, por sus siglas en Ingles) es una técnica basada en radar que representa una clase específica de las técnicas de Radar de Apertura Sintética Interferométrica Diferencial (DInSAR, por sus siglas en Ingles). Dichas técnicas usan la información contenida en la fase de la señal de al menos dos imágenes SAR, adquiridas en diferentes momentos sobre la misma área, que se utilizan para formar un interferograma (Crosetto et al, 2016).

Otro de los resultados que se pueden obtener con InSAR es la generación de series de tiempo. Estas técnicas se pueden dividir en aquellas que se basan en el análisis de interferogramas, todo con respecto a la misma imagen maestra, comúnmente conocida como métodos de dispersión persistente (Hooper, 2010).

Las series temporales de deformación representan el producto PSI más avanzado. Proporcionan el historial de deformación durante el período observado. Las medidas de deformación con PSI se refieren a la Línea de Vista (LOS, por sus siglas en inglés) del sensor SAR, es decir, la línea imaginaria que conecta el sensor con el objeto en la superficie (Crosetto et al, 2016).

Con esta técnica es posible observar detalladamente la morfodinámica costera a largo plazo (Sajinkumar et al., 2020), determinación de la subsidencia de la superficie (Osmanoğlu et al, 2011), estudios de la atmósfera como en (Haghighi, 2017).

Área de estudio

El presente caso de estudio se enmarca en el Caribe Sur (provincia de Limón), entre las coordenadas geográficas 9°44’18.28” - 10° 0’6.57” latitud norte y 82°47’27.54” - 83°11’32.24” de longitud oeste. En la Figura 1 se puede observar cuál es el área de estudio.

Figura 1. Área de estudio

Fuente: Elaboración propia.

Características Generales

El Caribe Sur de Costa Rica está geográficamente ubicado en la provincia de Limón y de acuerdo con la división que realiza el Ministerio de Planificación y Política Económica (MIDEPLAN), corresponde con la región Huetar Caribe. En esta región, el Índice de Desarrollo Social (IDS) de 2017 indica que de los 29 distritos que conforman la región Huetar Caribe, el 86% ocupa las posiciones de mayor desventaja social, destacándose los altos niveles de desempleo en comparación con la media del país (MIDEPLAN, 2018).

El clima de la región se caracteriza por no presentar una estación seca definida, mostrando dos periodos en que disminuyen las precipitaciones, el primero de febrero a marzo y el segundo entre setiembre y octubre (IMN, 2017).

El sistema fluvial del cantón de Limón corresponde a la subvertiente Caribe de la vertiente Caribe, el cual pertenece a las cuencas de los ríos Estrella, Matina, Banano, Moín y Bananito

En el tema económico, el turismo es una de las principales actividades del Caribe Sur destacando entre los lugares más atractivos algunas playas como Cahuita, Manzanillo y Puerto Viejo, lo que permite un turismo de naturaleza, el ecoturismo y el turismo rural (Morera-Beita & Miranda-Alvarez, 2016; Trejos & Chiang, 2009).

Desde un punto de vista socio-natural, la región Caribe Sur se encuentra expuesta a riesgos relacionados con amenazas naturales como lo son la sismicidad, las inundaciones, la erosión costera y los oleajes severos, entre otros aspectos que repercuten en su dinámica local y en su condición socioeconómica (Barrantes-Castillo et al., 2019).

En cuanto a la geología y la geomorfología, la costa Caribe puede dividirse en dos sectores: el norte y el sur. La parte norte presenta características de un margen continental pasivo y un régimen tectónico extensivo (fallas normales) mientras que el sur presenta un régimen compresivo con ocurrencia de terremotos importantes (Alvarado & Cárdenas, 2016), como el ocurrido en 1991 con epicentro en el Valle de la Estrella, el cual tuvo una magnitud de 7.7 Mw (Montero et al., 1994) el cual provocó un levantamiento máximo 1.85 m según se indica en (Denyer et al., 1994).

El origen de la ciudad de Limón comienza después de dos siglos de la llegada de Cristóbal Colón a la isla Uvita en 1502, cuando el 6 de junio de 1870 se decretó la creación de la comarca de Limón; un año después inician los trabajos de construcción del ferrocarril al Atlántico que unió a San José con Puerto Limón, con el fin de facilitar las exportaciones a Europa del café costarricense. Con la construcción de las obras del ferrocarril llegaron los primeros inmigrantes jamaiquinos, quienes aportaron su lengua, cultura, religión y gastronomía.

En la actualidad, la ciudad de Limón tiene una población aproximada de 105000 habitantes, posee el puerto de mayor importancia del país y se encuentra ubicada en el litoral Caribe de Costa Rica.

Marco teórico-conceptual

Para comprender cómo trabaja los radares de imágenes por satélite, primero se debe entender el concepto de un sensor remoto activo. Un sensor de microondas activo es un sensor que emite un patrón directo de energía que ilumina una porción de la superficie terrestre y luego recibe esta porción de vuelta al instrumento; estos sensores generan su propia energía, por lo que su uso es limitado por pocas restricciones y puede ser utilizado bajo una gran variedad de condiciones operacionales. Por otro lado, los sensores pasivos son sensibles a las variaciones de la iluminación solar, su uso es determinado por la hora del día y el clima (Campbell & Wynne, 2011).

Debido a la capacidad de los radares para teledetección, se pueden adquirir imágenes aun cuando no hay luz solar, a través de nubes y durante eventos climatológicas, lo que provee la capacidad de observar áreas de interés, que con otros tipos de sensores serian inobservables debido a su lejanía y condiciones atmosféricas (Campbell & Wynne, 2011).

El Radar Interferométrico (InSAR) utiliza dos imágenes SAR de la misma región adquirida en dos posiciones diferente. La figura 2 muestra el fundamento conceptual: se tiene una imagen obtenida en un tiempo t1 y una imagen en el tiempo t2. Se conoce para t1 y t2 la distancia entre el satélite y el objeto en la superficie de la Tierra (R1 y R2 respectivamente). A partir de la diferencia ∆R (siendo esta R1 – R2) se puede obtener alturas del terreno o determinar deformaciones verticales (Campbell & Wynne, 2011; Kerle, Janssen, & Huurneman, 2004; Soergel, 2010).

La técnica InSAR tiene grandes capacidades sin embargo, una de sus mayores limitaciones es la falta de coherencia de la imagen debido a la decorrelación espacio - temporal. La decorrelación espacial es el resultado de aumentar la línea base espacial (elemento B en la Figura 2), que es la línea que une las posiciones del sensor en dos tomas distintas mientras que la decorrelación temporal ocurre cuando las propiedades de dispersión de la superficie de la Tierra dentro de una celda cambian de una imagen a otra esto como consecuencia de que las tomas se hacen en momentos distintos (Zebker & Villasenor, 1992).

Como solución a estas limitaciones, se desarrollaron una variedad de propuestas de diferentes métodos basados en InSAR, en los cuales se utiliza objetos coherentes (objetos que no presentan cambios en el tiempo siendo estos objetos llamados dispersores persistentes (PS, por sus siglas en inglés) (Osmanoğlu et al., 2011). Esto ha llevado al uso de las técnicas Persistent Scatterer Interferometry (PSI), en particular el conocido algoritmo StaMPS, desarrollado en la Universidad de Stanford (Hooper et al., 2010).

A diferencia de InSAR, donde comúnmente solo se tiene un interferograma, PSI proporciona una serie temporal que representa el desplazamiento promedio de los objetos contenidos en el pixel, permitiendo así investigar la variación temporal en el patrón de subsidencia o levantamiento de la superficie.

Esta técnica se ha aplicado con éxito en diferentes estudios de deformación de la superficie, incluida la deformación de volcanes, sobreexplotación de aguas subterráneas, deformación en áreas de deslizamientos de tierra, análisis de presas y movimiento tectónico (Osmanoğlu et al., 2011; Crosetto et al., 2016).

Para el análisis PS-InSAR, las imágenes se corregistran conjuntamente con una imagen llamada imagen maestra. La imagen maestra se elige de modo que minimice los efectos de decorrelación espacial y temporal, minimizando entre el conjunto de imágenes la magnitud de las líneas base perpendiculares (Bper en la figura 2) y teóricamente hablando, maximizando la coherencia entre los interferogramas (Kampes, 2005). Luego, se generan los interferogramas y se eliminan la influencia de la topografía (fase topográfica) mediante el uso de un modelo de elevación.

La coherencia es una medición del grado de correlación entre un pixel en dos imágenes distintas. Este valor fluctúa entre 0 y 1. Si el valor de la coherencia es menor a 0.3, se considera que el interferograma está afectado por ruido, lo que impedirá aislar la señal de deformación (Kampes, 2005).

Figura 2. Geometría básica de una adquisición SAR

Fuente: Elaboración propia.

Dentro de un interferograma, la fase de cada píxel depende de los retornos de fase de todos los scatterers en el elemento correspondiente en el suelo. Algunos píxeles contienen objetos estables cuyo retorno en fase permanece constante a lo largo del tiempo a pesar de los cambios el área circundante (Hooper et al., 2007).

Marco metodológico

La metodología para estimar la velocidad en LOS de los PS se puede agrupar en tres conjuntos de procesos: descargar las imágenes, el preprocesamiento de estas y procesamiento final para la estimación de la velocidad. Para los fines del presente trabajo, el procesamiento de las imágenes se realizó en el programa SNAP, de la ESA, el cual es un programa de acceso libre y el proceso final se realizó usando StaMPS, consiste en un conjunto de scripts ejecutados desde el software Matlab. En la Figura 3 se muestra el esquema de procesos que se llevó a cabo.

Figura 3. Procedimiento seguido para la estimación de la velocidad en LOS

Fuente: Elaboración propia.

Previo a aplicar la metodología indicada anteriormente y debido a la inexistencia de información bibliografía sobre el uso de la técnica PSI en Costa Rica y a la curva de aprendizaje por parte de los autores, se realizó una prueba de sensibilidad del algoritmo, en específico se analizó la influencia de dos variables de cálculo que son:

a) La amplitud de dispersión (DA), parámetro que condiciona la cantidad de PS preliminares que serán considerados en el procesamiento. El DA está asociado con la calidad de la fase de cada PS

b) El “density rand”, el cual regula la cantidad de PS por kilómetro cuadrado. Entre mayor sea este valor, mayor cantidad de PS se seleccionarán.

Con el fin de realizar la prueba, se hicieron cuatro procesamientos, cada uno con los mismos 29 interferogramas. Los resultados se muestran en el apartado de resultados en la Tabla 1.

Se describe a continuación el procedimiento mostrado en la Figura 3: el primer proceso consiste en la descarga de las imágenes SLC (estas imágenes son las que contienen la información de la fase y por ende las requeridas para interferometría). Actualmente hay varios servicios Web desde los cuales se puede descargar las imágenes.

El oficial es el Data Hub de Copernicus, sin embargo, para este trabajo las imágenes se obtuvieron desde el servidor espejo gestionado por el ASF Alaska, por cuanto, después de varias pruebas, se determinó que el tiempo de descarga es menor en comparación con el servidor oficial.

Se descargaron las imágenes disponibles para los años 2016, 2017 y 2018; esto dio como resultado que se obtuvieron 62 imágenes con orbita descendente. El periodo que cubren las imágenes es desde el 06 de octubre de 2016 hasta el 25 de diciembre de 2018. Se inició con las imágenes desde octubre de 2016 por cuanto imágenes anteriores a esta fecha fueron adquiridas con geometría distintas (tipo de órbita y área de cobertura).

Tras analizar el procedimiento que se debe seguir, se decidió trabajar solamente con 30 imágenes, ya que el procesamiento de imágenes de Radar es un procesamiento intensivo desde un punto de vista computacional, se consideró que no se dispone de una computadora con la capacidad requerida para trabajar con las 62 imágenes.

Una vez descargadas las imágenes y seleccionadas cuales se procesarían, se procedió a preprocesar las en el programa SNAP. Este preprocesamiento se describe a continuación:

a) Se procede a seleccionar la imagen máster, usando para este fin la herramienta “Stack Overview”. Como resultado, se obtiene que la imagen máster es la imagen del 11 de enero de 2018

b) El siguiente proceso es aplicar el split a las imágenes y corregirlas por el efecto de desviaciones entre la órbita teórica y la órbita real del satélite. Con el split se reduce la imagen al área de interés, por cuanto una imagen SLC completa cubre 250 km por 250 km, así con el split se selecciona para trabajar únicamente el área de interés En cuanto a las órbitas, este proceso se realiza para corregir el efecto de las perturbaciones de la órbita del satélite, las cuales son causadas principalmente por el campo de gravedad de la Tierra. Para aplicar esta corrección se requiere de conexión a Internet, por cuanto hay que descargar los archivos con la información de la órbita real seguida por el satélite

c) Una vez seleccionada el área de trabajo y corregido el efecto de las orbitas, se genera un stack. Esto consiste en generar un solo archivo que contiene a la imagen máster y el resto de las imágenes (hasta este paso cada imagen es un archivo individual). Debido a las limitaciones computacionales ya comentadas, se decidió dividir este proceso en tres, es decir, generar tres archivos: dos contienen las imágenes master y diez imágenes más, mientras que uno contiene la imagen master y nueve imágenes

d) Al resultado del proceso anterior se le aplicó un deburst. Este proceso se debe llevar a cabo por la forma como el sensor en el satélite realiza el barrido de la superficie, mismo que depende del modo de escaneo. Este efecto se observa en la Figura 4 y consiste en las franjas horizontales. Una vez aplicada esta corrección, dichas franjas no se observan más en las imágenes.

Figura 4. Ejemplo previa aplicación del deburst

Fuente: Elaboración propia a partir del procesamiento de una imagen SLC en el programa SNAP.

e) El siguiente proceso es generar los interferogramas. En la figura 5 se muestra un ejemplo de un interferograma.

Figura 5. Ejemplo de interferograma

Fuente: Elaboración propia a partir del procesamiento del programa SNAP.

Como se desprende de la Figura 5, el interferograma no muestra ningún patrón evidente de deformación. Esto se debe a que típicamente en presencia de patrones de deformación lentos y por la longitud de onda del satélite, no es posible separar la señal de deformación del ruido. Es por ello que se recurre a métodos de DInSAR para generar series de tiempo, ya que de esta forma si es posible obtener información sobre la deformación de la superficie. El interferograma de la Figura 5 fue realizado con las imágenes del 23 de noviembre de 2016 y el 11 de enero de 2018.

f) El último paso del preprocesamiento es remover la llamada fase topográfica. Esta fase se denomina así, debido a que está relacionada con el efecto de la topografía en el barrido que hace el sensor, por lo que para removerla se requiere de un modelo de elevación. En este trabajo y para este fin, se usó el modelo de la misión Shuttle Radar Topography Mission (SRTM) con una resolución de 3 segundos de arco. Una vez removida la fase topográfica, se considera que la fase interferogramétrica contiene la señal de deformación.

Una vez realizados los pasos anteriores, se procede a preparar la información para procesar la con StaMPS. El programa SNAP tiene incorporada la opción de exportar los resultados al formato requerido por StaMPS, por lo que se utilizó esta opción. Una vez se dispone de la información, se realiza el procesamiento, mismo que consiste en siete procesos, los cuales se describen a continuación:

a) Estimar el ruido de la fase: Este es un paso iterativo que estima el valor del ruido de la fase para cada píxel candidato en cada interferograma.

b) Selección de los PS: Los pixeles son seleccionados en base a sus características de ruido. En este paso también se estima el porcentaje de pixeles aleatorios (no-PS) en una imagen del cual la densidad por kilómetro cuadrado pudo ser obtenida. El fundamento para este paso es que los pixeles ruidosos pueden causar dificultades en el procesamiento y dadas sus características, la contribución no puede ser suavizada usando pixeles de bajo ruido (con fase estable). En este paso se tiene la opción de descartar los pixeles ruidosos si estos exceden un umbral definido por el usuario.

c) Remover los PS: Este paso toma como información de entrada el resultado del paso anterior. Los pixeles que se han identificado como ruidosos son eliminados, de forma que tras este paso quedan los pixeles que cumplen con la calidad definida por el usuario.

d) Corregir la fase: La fase de los pixeles resultantes del paso anterior son corregidos por el error asociado al modelo de elevación, denominado SULA. Este error está asociado a que la fase esta correlacionada con la topografía, por lo que se estima este efecto y se elimina de la fase.

e) Unwrapping de la fase: Una de las limitaciones de InSAR es que la fase se mide en ciclos desde -π hasta +π o desde 0 hasta +2π según la convención. Sin embargo, para estimar la deformación real, se debe sumar la cantidad de ciclos en el interferograma. Este proceso se denomina unwrapping y para este caso se utilizó el programa Snaphu.

f) Estimación del AOE: Finalmente, el último paso es estimar el error de la atmósfera y órbita (AOE) de la imagen máster.

g) Generación de las series de tiempo: una vez concluido el procesamiento StaMPS, se exportan los resultados a un archivo formato csv separado por comas. Para la visualización de las series de tiempo, se usó la aplicación Stamps-Visualizer, la cual es una aplicación que se ejecuta en el programa R.

Resultados

Para facilitar la lectura, los resultados están divididos en bloques: el primero hace referencia a una prueba de sensibilidad del algoritmo, donde se modificaron dos parámetros con el fin de determinar la influencia de estos en los resultados. El siguiente punto, muestra la forma como se estimó la velocidad en LOS de la estación LIMN, única estación GNSS dentro del área de estudio, la cual aporta un valor para identificar un punto de referencia. Finalmente, se muestran los principales resultados del procesamiento.

Prueba de sensibilidad del algoritmo

En la tabla 1 se muestran los resultados de la prueba de sensibilidad del algoritmo.

Tabla 1. Sensibilidad del algoritmo

Parametro	Calculo_1	Calculo_2	Calculo_3	Calculo_4
PS_candidates	929315	929315	636484	636484
PS_selected_initially	523793	495093	512403	486303
PS_after_reestimation_coherence	154447	124982	147187	121375
PS_kept_after_dropping_noisy_pixels	44174	41421	37855	35997
ref_PS_selected	44174	41421	37855	35997

Fuente: Elaboración propia.

De la Tabla 1:

a) PS_candidates es el número de PS candidatos

b) PS_selected_initially: es la cantidad de PS seleccionados luego del paso 2 en StaMPS

c) PS_after_reestimation_coherence: es la cantidad de PS seleccionados luego del paso 3 en StaMPS

d) PS_kept_after_dropping_noisy_pixels: es la cantidad de PS seleccionados luego del paso 4 en StaMPS

e) ref_PS_selected: es la cantidad de PS seleccionados luego del paso 5 en StaMPS

El Gráfico 1 muestra los resultados de los valores dados en la Tabla 1. Se observa cómo tras cada paso en el procesamiento, la cantidad de PS disminuye.

Gráfico 1. Resultados de la prueba de sensibilidad del algoritmo

Fuente: Elaboración propia.

Estimación de la velocidad en LOS de la estación GNSS LIMN

PSI es un método relativo, por lo que, para la estimación de las velocidades, se le debe indicar al algoritmo un punto el cual se considera estable o al menos del cual se conoce su variación en el tiempo. Para el presente estudio se usó la estación GNSS LIMN, la cual está ubicada en las instalaciones del Banco de Costa Rica en Limón. De esta estación se conoce los componentes de la velocidad (Valverde-Calderón, 2020). La velocidad geocéntrica fue proyectada a LOS aplicando la fórmula dada en (Hanssen, 2001):

Donde:

dr: velocidad en LOS

du: velocidad en up

de: velocidad en este

dn: velocidad en norte

θinc: ángulo de incidencia

αh: azimut del satélite.

Para el cálculo se usaron los siguientes valores:

• θinc = 33.94°, obtenido del archivo ps1.mat, generado por StaMPS

• α = -167.98°, obtenido de los metadatos de la imagen master

• Latitud: 9.9930920° N

• Longitud: -83.0263684° W

Como resultado se obtuvo una velocidad en LOS de +3.0290 mm/yr.

Aplicación del método PSI en Limón y estimación de las velocidades en LOS

Una vez se analizaron los resultados de la prueba de sensibilidad del algoritmo y se estimó la velocidad de la estación LIMN, se procedió a aplicar el método PSI usando el algoritmo StaMPS.

En el cálculo inicial, se trabajó con 30 imágenes, lo que permite la formación de 29 interferogramas. En la Tabla 2 se muestra la información de cada interferograma (identificado con la fecha) en relación con la imagen master: Bperp, que es la magnitud de la línea base perpendicular (según se explicó en el apartado del marco conceptual), Btemp, que es la cantidad de días entre la adquisición de la imagen i y la imagen master y ρ que es un estimado de la coherencia entre la imagen i y la imagen master.

Tabla 2. Información de cada interferograma

Fecha	Bperp (m)	Btemp (días)	ρ	Fecha	Bperp (m)	Btemp (días)	ρ
2016-11-23	-35.07	414.0	0.60	2018-02-04	54.37	-24.0	0.93
2017-01-10	47.17	366.0	0.64	2018-02-28	43.28	-48.0	0.92
2017-02-09	67.22	336.0	0.65	2018-03-24	-95.01	-72.0	0.86
2017-03-05	-84.28	312.0	0.66	2018-04-17	-30.24	-96.0	0.89
2017-04-10	72.71	276.0	0.70	2018-05-11	59.16	-120.0	0.85
2017-05-04	78.46	252.0	0.72	2018-06-04	-65.59	-144.0	0.82
2017-05-28	-42.72	228.0	0.76	2018-06-28	20.44	-168.0	0.83
2017-06-21	35.68	204.0	0.79	2018-07-22	11.65	-192.0	0.82
2017-07-15	-24.40	180.0	0.82	2018-08-15	-18.10	-216.0	0.79
2017-08-08	33.20	156.0	0.83	2018-09-08	33.31	-240.0	0.76
2017-09-01	22.50	132.0	0.86	2018-10-02	-33.12	-264.0	0.74
2017-09-25	21.01	108.0	0.88	2018-10-26	14.68	-288.0	0.73
2017-10-19	-11.82	84.0	0.91	2018-11-19	4.37	-312.0	0.71
2017-11-12	-4.82	60.0	0.94	2018-12-13	-62.47	-336.0	0.66
2017-12-06	-54.83	36.0	0.92

Fuente: Elaboración propia.

La coherencia cambia de píxel a píxel y de interferograma a interferograma, por lo que se generó la banda de coherencia, aplicando un filtro que deje los pixeles con coherencia mayores a 0.6, esto con el fin de identificar interferogramas ruidosos. En la Figura 6 se muestra la banda de coherencia para el interferograma formado con las imágenes del 23 de noviembre de 2016 y el 11 de enero de 2018 (esta es la imagen master).

Figura 6. Ejemplo de coherencia en el casco urbano de Limón

Fuente: Elaboración propia a partir del procesamiento en el programa SNAP.

Otro de los resultados son los interferogramas. Estos son una representación del cambio de fase entre la imagen usada y la imagen máster. Se llama interferograma enrollado cuando la fase está en el rango de -π a +π. En la figura 7 se muestra un ejemplo de un interferograma enrollado elaborado con las imágenes del 04 de febrero de 2018 y la imagen máster, como se observa en la simbología la escala va de -3.1 a 3.1 siendo estos valores la fase de cada PS.

Figura 7. Ejemplo de interferograma enrollado

Fuente: Elaboración propia a partir del procesamiento en StaMPS.

Se llama interferograma desenrollado cuando la fase está en el rango de -n*π a +n*π, donde n es un número entero. Este número hace referencia a la cantidad de ciclos completos acumulados. En la Figura 8 se muestra el interferograma desenrollado para las imágenes del 04 de febrero de 2018 y la imagen máster, donde los valores en la escala de colores representan la cantidad de ciclos acumulados para ese PS.

Figura 8. Ejemplo de interferograma desenrollado

Fuente: Elaboración propia a partir del procesamiento en StaMPS.

Tras aplicar el procesamiento comentado en el marco metodológico, el resultado principal es el mapa que muestra la velocidad en LOS para los PS en el área de trabajo, en milímetros por año (ver Figura 9). De la anterior figura las áreas de interés en tonos amarillos representan velocidades negativas en LOS, mientras que los valores que tienden al azul son velocidades positivas en LOS.

Figura 9. Velocidad en LOS

Fuente: Elaboración propia a partir del procesamiento en StaMPS.

En la Figura 10 se muestra el error en la velocidad en LOS, también en milímetros por año. En este caso los PS en color rojo son los que tienen menor error y conforme cambia el color hacia los tonos azules el error crece.

Figura 10. Error de la velocidad en LOS

Fuente: Elaboración propia a partir del procesamiento en StaMPS.

Finalmente, junto con las velocidades para el área de trabajo, se pueden generar la serie de tiempo para cada PS de interés. En la Figura 11 se muestra un ejemplo de la serie de un PS, generada usando la aplicación StaMPS-Visualizer (Höser, 2018).

Figura 11. Ejemplo de serie de tiempo de un PS

Fuente: Elaboración propia a partir del uso de la aplicación StaMPS-Visualizer.

Discusión de resultados

El presente artículo busca compartir la primera experiencia en el país en el uso de imágenes de radar para determinar deformaciones de la superficie. La determinación y el análisis de las deformaciones de la superficie han sido realizados en varios proyectos, como por ejemplo en Valverde-Calderón (2020) pero, se ha limitado al uso de estaciones GNSS, con las limitaciones de estos métodos relacionados con la densidad espacial de las estaciones y la disponibilidad de datos.

Se logró estimar la velocidad en LOS para el área de interés. El rango de velocidades va desde -11 mm/yr hasta +20 mm/yr (recordar que es la velocidad en LOS, no la velocidad vertical), con errores que van desde los 0.7 mm/yr hasta los 13 mm/yr.

Al inicio se trabajó con 30 imágenes y 29 interferogramas, sin embargo, conforme se fue avanzando en el flujo de trabajo, se fueron descartando los interferogramas con más ruido, por lo que los resultados mostrados son los generados a partir de 22 interferogramas, que cubren el periodo desde el 10 de enero de 2017 hasta el 18 de noviembre de 2018.

Luego del procesamiento, se identificaron las bondades y limitaciones de la técnica. Teóricamente la técnica permite determinar deformaciones de áreas de interés con un gran detalle sin la necesidad de colocar instrumental en la superficie terrestre y con una precisión comparable a los métodos geodésicos disponibles, sin embargo, se hace necesario disponer de información de otras técnicas como los métodos GNSS, a efectos de tener información para la validación de los resultados; en este proyecto se utilizó como estación de referencia la estación LIMN, sin embargo, no se dispone de otras estaciones GNSS para la validación de los resultados, lo que supone un desafío para futuras aplicaciones.

La meta original del proyecto es obtener a partir de la técnica PSI desplazamientos verticales; sin embargo, debido a que solo se procesaron imágenes con órbitas descendentes, solo se pudo estimar la velocidad en LOS. Para futuros trabajos, se debe procesar imágenes de órbita ascendente y descendente y combinar los resultados de estos procesamientos para obtener la velocidad de la deformación de la superficie en el componente altura y en la dirección este – oeste.

En el presente proyecto no se corrigió el efecto la Pantalla de Fase Atmosférica (Atmospheric Phase Screen o APS, por sus siglas en ingles). Este efecto consiste en el error que introduce la atmósfera sobre la señal emitida por el satélite y depende del contenido de vapor de agua en la atmósfera. En una segunda etapa del proyecto que dio origen a esta investigación se corregirá dicho efecto.

De acuerdo con Barrantes-Castillo et al., (2019) en el Caribe sur de Costa Rica se presenta un proceso extendido de erosión costera. De acuerdo con estos autores la geodinámica local puede ser una de las causas de este fenómeno al sumarse al ascenso de nivel del mar, similar al que juega en la costa de la Península de Nicoya. En el caso del Caribe Sur, aunque la geodinámica y tectónica de este margen no se corresponde con una zona de subducción propiamente, no se ha estudiado con el detalle el mecanismo de acumulación de esfuerzos que podrían producir la subsidencia de la costera en los periodos entre grandes sismos que levantan la costa. En este sentido la técnica utilizada puede contribuir en la compresión de la geodinámica local o al menos cuantificar su posible contribución al proceso de erosión costera.

Conclusiones

Se logró estimar la velocidad en LOS para el área de interés, lo que permite generar una primera aproximación sobre el comportamiento de la superficie terrestre en el área de interés. Sin embargo, debido a que el mayor movimiento de la corteza es en el componente horizontal, no es posible con estos resultados estimar la velocidad en el componente vertical. Para estimar la velocidad vertical, es necesario el combinar resultados de procesamientos donde se usen imágenes con órbita ascendente y descendente (Hu et al., 2014).

A pesar de que la bibliografía recomienda que se debe usar mínimo 20 interferogramas, tras el procesamiento y siendo una posibilidad el tener que descartar información durante el proceso de los datos, se concluye que se debe procesar la mayor cantidad posible de imágenes.

Sin embargo, esto representa un desafío desde el punto de vista computacional, ya que entre más imágenes se usen en el proceso, se requiere de mayor cantidad de recursos, como memoria y espacio en el disco duro. Por ejemplo, para este proyecto se usó una computadora con 16 GB de memoria RAM, cantidad que en algunos procesos resulto ser el valor mínimo de operación, de forma que los cálculos se pudieran realizar.

El uso de estos métodos tiene un potencial para ser usados en Costa Rica, particularmente por la dinámica de la superficie. En este aspecto, la Escuela de Topografía, Catastro y Geodesia, junto con la Escuela de Ciencias Geográficas, a través del proyecto FIDA 0433-18 “Erosión costera, geodinámica regional y gobernanza para la gestión del riesgo socioambiental en el Caribe Sur de Costa Rica” aplicarán este método con el fin de generar insumos para el entendimiento del proceso erosivo en el Caribe Sur de Costa Rica.

Referencias bibliográficas

Alvarado, G. & Cárdenas, G. (2016). Chapter 3 Geology, Tectonics, and Geomorphology of Costa Rica: A Natural History Approach. En M. Kappelle, Costa Rican Ecosystems (p. 744). London.

Balzter, H. (2001). Forest mapping and monitoring with interferometric synthetic aperture radar (InSAR). Progress in physical geography, 25(2), pp. 159-177.

Barrantes-Castillo, G., Arozarena-Llopis, I., Sandoval-Murillo, L. F., & Valverde-Calderón, J. F. (2019). Playas críticas por erosión costera en el caribe sur de Costa Rica, durante el periodo 2005-2016. Revista Geográfica de América Central, 1(64), pp. 95-122. https://doi.org/10.15359/rgac.64-1.4

Campbell, J. B. & Wynne, R. H. (2011). Introduction to remote sensing. Guilford Press. Chicago Press.

Cho, B. L., Kong, Y. K., Park, H. G., & Kim, Y. S. (2006). Automobile-based SAR/InSAR system for ground experiments. IEEE Geoscience and Remote Sensing Letters, 3(3), pp. 401-405.

Crosetto, M., Monserrat, O., Cuevas-González, M., Devanthéry, N., & Crippa, B. (2016). Persistent scatterer interferometry: A review. ISPRS Journal of Photogrammetry and Remote Sensing, 115, pp. 78-89.

Delgado-Blasco, J. M., Foumelis, M., Stewart, C., & Hooper, A. (2019). Measuring urban subsidence in the rome metropolitan area (italy) with sentinel-1 snap-stamps persistent scatterer interferometry. Remote Sensing, 11(2), p. 129. https://doi.org/10.3390/rs11020129

Denyer, P., Arias, O., & Personius., S. (1994). Efectos tectonicos del terremoto de Limon, Costa Rica. Rev. Geol. Amer. Central, Volumen especial, pp. 39-52.

Haghighi, M., & Motagh, M. (2017). Sentinel-1 InSAR over Germany: Large-scale interferometry, atmospheric effects, and ground deformation mapping. ZfV: Zeitschrift für Geodäsie, Geoinformation und Landmanagement, 2017(4), 245-256.

Hanssen, R.F. (2001). Radar interferometry: data interpretation and error analysis (Vol.2). Springer Science & Business Media.

Homer, J., Longstaff, I. D., & Callaghan, G. (1996). High resolution 3-D SAR via multi-baseline interferometry. In IGARSS’96. 1996 International Geoscience and Remote Sensing Symposium (Vol. 1, pp. 796-798). IEEE.

Hooper, A. (2010). A statistical-cost approach to unwrapping the phase of InSAR time series. In Proceedings of the International Workshop on ERS SAR Interferometry, Frascati, Italy (Vol. 30).

Hooper, A., Bekaert, D., & Spaans, K. (2010). StaMPS/MTI manual. Delft Institute of Earth Observation and Space Systems Delft University of Technology, Kluyverweg, 1, p. 2629.

Hooper, A., Segall, P., & Zebker, H. (2007). Persistent scatterer InSAR for crustal deformation analysis, with application to Volcán Alcedo, Galápagos. Journal of Geophysical Research, 112(B07407), p. 19.

Höser, Thorsten. (2018). Analysing the Capabilities and Limitations of InSAR using Sentinel-1 Data for Landslide Detection and Monitoring. https://doi.org/10.13140/rg.2.2.35085.59362

Hu, J., Li, Z. W., Ding, X. L., Zhu, J. J., Zhang, L., & Sun, Q. (2014). Resolving three-dimensional surface displacements from InSAR measurements: A review. Earth-Science Reviews, 133, pp. 1-17.

Instituto Meteorológico Nacional [IMN]. (2017). Clima de Costa Rica: el clima y las regiones climáticas de Costa Rica. Recuperado el 07 de noviembre de 2017, de Instituto Meteorológico Nacional: https://www.imn.ac.cr

Kampes, B. M. (2005). Displacement Parameter Estimation using Permanent Scatterer Interferometry, Ph.D. thesis, Delft University of Technology.

Kerle, N., Janssen, L. L., & Huurneman, G. C. (2004). Principles of remote sensing. ITC, Educational textbook series, 2, p. 250.

Ministerio de Planificación Nacional y Política Económica. [MIDEPLAN]. (2018). Área de Análisis del Desarrollo. Índice de desarrollo social 2017 / Ministerio de Planificación Nacional y Política Económica. -- San José, CR: MIDEPLAN, 2018.

Montero, W., Pardo, M., Ponce, L., Rojas, W., & Fernández, M. (1994). Evento principal y replicas importantes del terremoto de Limón. Revista Geológica de América Central. Vol. espec. Terremoto de Limón, pp. 93-102.

Morera-Beita, C. & Miranda-Álvarez, P. (2016). De la geografía del turismo al análisis territorial del turismo: el rastro en Costa Rica. Revista Geográfica de América Central, 1(54), 15-43. https://doi.org/10.15359/rgac.1-54.1

Municipalidad de Limón. (2019). Municipalidad del Cantón Central de Limón. Obtenido de Historia: http://www.municlimon.go.cr/index.php/mn-conozcanos/mn-micanton/mn-historiacanton

NASA. (2014). EcoSAR. Recuperado el 29 de Enero de 2019, de https://ecosar.gsfc.nasa.gov/campaigns/costa-rica

Osmanoğlu, B., Dixon, T. H., Wdowinski, S., Cabral-Cano, E., & Jiang, Y. (2011). Mexico City subsidence observed with persistent scatterer InSAR. International Journal of Applied Earth Observation and Geoinformation, 13(1), pp. 1-12.

Sajinkumar, K. S., Bincy, H. S., Bouali, E. H., Oommen, T., Vishnu, C. L., Anilkumar, Y. & Keerthy, S. (2020). Picturing beach erosion and deposition trends using PSInSAR: an example from the non-barred southern west coast of India. Wetlands Ecology and Management, pp. 1-14.

Samieie-Esfahany, S., Hanssen, R., van Thienen-Visser, K., & Muntendam-Bos, A. (2009). On the effect of horizontal deformation on InSAR subsidence estimates. In Proceedings of the Fringe 2009 Workshop, Frascati, Italy (Vol. 30).

Smith, L. C. (2002). Emerging applications of interferometric synthetic aperture radar (InSAR) in geomorphology and hydrology. Annals of the Association of American Geographers, 92(3), pp. 385-398.

Soergel, U. (2010). Review of radar remote sensing on urban areas. In Radar remote sensing of urban areas. Springer, Dordrecht, pp. 1-47.

Trejos, B. & Chiang, L. H. N. (2009). Local economic linkages to community based tourism in rural Costa Rica. Singapore journal of tropical geography, 30(3), pp. 373-387.

U.S. Geological Survey. (2019). USGS science for a changing world. Obtenido de Mapping, Remote Sensing, and Geospatial Data: https://www.usgs.gov/faqs/what-remote-sensing-and-what-it-used?qt-news_science_products=0#qt-news_science_products

Uys, Duan. (2016). InSAR: an introduction. Preview. 182, 43-48 https://doi.org/10.1071/PVv2016n182p43

Valverde-Calderón, J. (2020). Estudio del efecto de un terremoto sobre un marco geodésico de referencia. Uniciencia, 34 (1), 1-2019, pp. 1-19. http://dx.doi.org/10.15359/ru.34-1.1

Yhokha, A., Goswami, P.K., Chang, CP. et al. (2018). Application of Persistent Scatterer Interferometry (PSI) in monitoring slope movements in Nainital, Uttarakhand Lesser Himalaya, India. J Earth Syst Sci 127, 6 https://doi.org/10.1007/s12040-017-0907-y

Zebker, H., & Villasenor, J. (1992). Decorrelation in interferometric radar echoes. IEEE Transactions on Geoscience and Remote Sensing 30 (5), pp. 950-959.

---