I. INTRODUCCIÓN

En el 2005 se diseñó, midió y procesó una red geodésica nacional de 34 vértices; esto mediante la metodología GPS. Este conjunto de vértices constituyó la Red Geodésica Nacional de Primer Orden y materializó el Marco Geodésico Nacional de Referencia de Costa Rica denominado CR05 Figure 1. A su vez, el diseño de esta red consideró observaciones de vínculo a un conjunto de estaciones internacionales GNSS pertenecientes al Sistema de Referencia Geocéntrico para las Américas (SIRGAS) y al Servicio Internacional de GNSS (IGS). Lo anterior permitió definir las coordenadas de CR05 vinculado con el Marco Internacional Terrestre de Referencia (ITRF2000), reducido a la época de referencia 2005,83. La exactitud general lograda en las coordenadas geodésicas ajustadas de estos vértices estuvo entre ± 1 cm y ± 3 cm. El CR05 fue oficializado como sistema de referencia geodésico nacional mediante el Decreto Ejecutivo número 33797-MJ-MOPT del 30 de marzo de 2007 y vigente a partir del 06 de junio del 2007 [6]. De los 34 vértices que definieron el CR05, 33 fueron pasivos y solamente el vértice ETCG, localizado en la Escuela de Topografía, Catastro y Geodesia de la Universidad Nacional, se contempló como única estación GPS de operación continua Figure 1.

Asimismo, a inicios del 2018 el Instituto Geográfico Nacional de Costa Rica (IGNCR), mediante el Decreto Ejecutivo número D-40962-MJP del 17 de abril de 2018, actualizó el Marco Geodésico Nacional de Costa Rica a través el denominado CR-SIRGAS; este corresponde a un marco de referencia moderno vinculado directamente con el al ITRF2008 por medio de ocho estaciones GNSS de operación continua administradas por el IGNCR y, que a su vez, están integradas a la Red GNSS de Operación Continua de SIRGAS (SIRGAS-CON) [7]. CR-SIRGAS se redujo a la semana 1803, específicamente a la época de referencia 2014,59 [8]. Con la implementación de CR-SIRGAS se da un avance importante en el país, debido a que, por primera vez, se cuenta con una referencia geodésica cinemática, ya que las coordenadas geocéntricas de las estaciones de referencia son calculadas semanalmente por los diferentes Centros de Procesamiento y Combinación de SIRGAS y vinculadas con la solución más reciente del marco internacional. La nueva Red Geodésica Nacional de Primer Orden consta de 42 vértices, ocho de ellos son estaciones activas y el resto corresponde a vértices pasivos Figure 1. Adicionalmente, dentro de las ventajas en la implementación de CR-SIRGAS se encuentra su actualización constante, la cual permite al IGNCR adecuar e incorporar a esta referencia geodésica nacional los cambios y actualizaciones que sucedan a nivel internacional por medio de la publicación de resoluciones y directrices técnicas.

2. MARCO TEÓRICO

A. Relación entre marcos geodésicos de referencia

La relación entre dos marcos geodésicos de referencia se establece por medio de un conjunto de siete parámetros que pueden derivarse empíricamente a partir de las coordenadas cartesianas geocéntricas de un grupo de puntos comunes a ambos marcos, llamados generalmente puntos homólogos o idénticos. En la Figure 2 se representa esta situación en la que se tienen dos marcos geodésicos de referencia denominados como marco A (color azul) y marco B (color negro) y un punto P (color gris), el cual posee coordenadas tridimensionales en los dos marcos.

Respecto a la transformación estándar entre dos marcos, esta es una similitud Euclidiana de siete parámetros [10] y [11]; los cuales relacionan a los dos marcos mediante tres valores lineales que representan las diferencias del geocentro de un marco origen respecto al marco destino denominadas traslaciones (DX, DY, DZ); tres valores angulares dados por los giros de cada uno de los tres ejes coordenados del marco origen denominados rotaciones (R_X, R_Y, R_Z), cuyo resultado genera que los dos marcos geodésicos sean paralelos. Y finalmente una relación de proporcionalidad denominada factor de escala (m), aplicado al marco origen, de manera que ambos tengan el mismo tamaño. Este parámetro de escala se expresa generalmente en unidades de milímetros por kilómetro [ppm]. Conocidos los siete parámetros de transformación, el marco de referencia origen A se traslada, se rota y se escala para llevarlo al marco de referencia destino B.

B. Modelo general de transformación

La transformación de sistemas coordenados tridimensionales ha cobrado mayor relevancia y aplicabilidad desde que las diferentes técnicas satelitales hicieron posible relacionar sistemas locales y globales [12]. Cualquier sistema cartesiano puede ser transformado a otro sistema cartesiano a través de las rotaciones respectivas. Lo anterior solo si sus orígenes son iguales y si ambos sistemas están definidos de la misma manera; por ejemplo: sistemas de mano derecha o de mano izquierda [13] (ver ecuación (1)).

Dicho modelo de transformación asume que el factor de escala es el mismo en todas las direcciones, lo cual implica semejanza o similitud [14] entre los dos marcos a relacionar, esto quiere decir que, tras haber aplicado el proceso de transformación, no hay distorsión angular. El modelo es conocido como transformación de similitud o semejanza de Helmert; a su vez, de acuerdo con algunos trabajos y adecuaciones posteriores, también se le llama con el nombre de transformación Bursa-Wolf [15] y [16].

A partir de los siete parámetros de transformación (ver Figure 2), el modelo matemático que relaciona los dos marcos de referencia se muestra en la ecuación 1 [13], [17], [18] y [19]. Ahí, el vector T cuenta con los valores de las tres traslaciones, m es el factor de escala general y R es una matriz de rotación ortogonal.

Por su parte, en la ecuación (2) se exponen las componentes tanto del vector de traslaciones T como las de la matriz de rotación R, cuya suma de cuadrados de cualquier fila o columna es igual a la unidad y su inversa es igual a la transpuesta R^-1 = R^T.

En cuando a convención, se asume que los ángulos de giro son pequeños y positivos en sentido antihorario cuando se mira a lo largo del eje en dirección del origen del sistema [12] y [13]. En el conjunto de ecuaciones (3) se presentan los componentes de las matrices de rotación aplicadas a cada uno de los tres ejes coordenados X, Y, Z indicados con los números 1, 2 y 3 respectivamente; y cuyo producto matricial R₃, R₂, R₁ dará como resultado la matriz de rotación R [20] y [21].

Al expresar los ángulos de giro en radianes, el valor del coseno será aproximadamente la unidad, el valor del seno será aproximadamente igual al ángulo y el producto de los senos de dos ángulos es prácticamente nulo. Debido a estas propiedades, los elementos de la matriz de rotación R se pueden simplificar en la forma presentada en la ecuación (2)[22] y [13]. En este documento el cálculo de los parámetros de rotación se hizo originalmente en radianes; sin embargo, los resultados se presentan en unidades de miles de arco segundos [mas] (1 más = 0,001”).

Este modelo de transformación entre sistemas de referencia es empleado para relacionar las diferentes soluciones globales del Marco Internacional Terrestre de Referencia ITRF desde su primera solución [23]. Desde entonces, se ha establecido como una convención en la relación de las diferentes soluciones ITRF [10], [11] y [24] hasta la más reciente solución ITRF2014 [25]. En el caso específico de los conjuntos de coordenadas provistos por las técnicas de la geodesia espacial, se asume la forma linealizada de la matriz de rotaciónR [24], tal como se muestra en la ecuación (4).

Así, al sustituir los componentes del vector de traslaciones y de la matriz de rotación (ecuación (2)) y, posteriormente desarrollar las operaciones respectivas, el cálculo de las coordenadas del marco A (origen) en el marco B (destino) se ofrece en la ecuación (5).

La ampliación de este modelo de transformación integrando las variaciones temporales de sus parámetros es el empleado en las diferentes soluciones del ITRF (ver ecuación (1)). De esta manera, se pasa de un modelo estático a uno dinámico de 14 parámetros [26].

En la práctica este modelo de transformación desprecia los términos de segundo orden (diferenciales del parámetro de escala y de la matriz de rotación respecto al vector de coordenadas del sistema origen) ya que sus valores representan unos 0,1 mm en 100 años, el factor de escala y las rotaciones son del orden de 10^-5 y las componentes del vector de velocidades es de unos 10 cm por año, como se muestra en la ecuación (7)[24] y [25].

C. Modelo de transformación baricéntrico

Si se aplica el modelo de transformación descrito anteriormente en regiones limitadas y no de manera global se obtiene una alta correlación entre los parámetros de traslación y rotación; para eliminar este efecto se acostumbra a reducir las coordenadas del marco origen a su centro de gravedad o baricentro. Esta variante es conocida como el modelo de transformación Molodensky-Badekas [3] y [4]. Asimismo, el modelo baricéntrico es recomendado para regiones. Por tanto, se considera adecuado para Costa Rica, pues el país se ubica a una latitud promedio de 10° norte y presenta una extensión territorial continental aproximada de 51100 Km².

Respecto a las coordenadas cartesianas geocéntricas del baricentro (X₀, Y₀, Z₀), estas se calculan como el promedio aritmético simple. Por ejemplo, en una zona con n puntos considerados las coordenadas del baricentro se calculan tal como se muestra en la ecuación (8), donde el superíndice hace referencia al marco origen A.

Por su parte, la relación matemática del modelo de transformación baricéntrico se presenta en la ecuación (9) y la forma desarrollada para cada una de las coordenadas se expone en la ecuación (10). Se aprecia en ambas ecuaciones la reducción de las coordenadas del sistema a transformar u origen a su respectivo baricentro.

Respecto a la ecuación (10), se tienen como datos las coordenadas cartesianas geocéntricas de los puntos idénticos, mientras que las incógnitas son precisamente los siete parámetros de transformación, los cuales se agrupan en el denominado vector de incógnitas que se presenta de forma transpuesta en la ecuación (11). La solución de este vector requiere conocer como mínimo tres puntos idénticos. No obstante, cuando el sistema de ecuaciones (10) está sobre determinado, la solución se obtiene a través de un ajuste por mínimos cuadrados.

El modelo matemático del ajuste por mínimos cuadrados implica plantear las observaciones en función de las incógnitas siguiendo. Por ejemplo, para una cantidad p de puntos idénticos, sus ternas coordenadas cartesianas geocéntricas serán asumidas como las observaciones, siendo 3p ecuaciones en total. Las diferencias en las coordenadas de los puntos idénticos (ver ecuación (12)) se asumen como errores casuales o aleatorios, los cuales son compensados en el proceso de ajuste.

Las 3p ecuaciones de observación deben linealizarse y los respectivos coeficientes se agrupan en la matriz de configuración A, tal como se muestra en la ecuación (13)[17] y [18].

La solución del vector de incógnitas se exhibe en la ecuación (14) el cual está en función de la matriz de configuración A, de la matriz de pesos P y del vector de observaciones L que contiene las diferencias entre las ternas de coordenadas del marco destino menos el marco origen. El modelo estocástico asume generalmente una matriz de pesos identidad P = I.

La desviación estándar de la unidad de pesos a posteriori se determina considerando un total de f = 3p – 7 de grados de libertad (ver ecuación (15)).

En aplicaciones que relacionen marcos de referencia determinados por medio de mediciones satelitales y convencionales, es necesario emplear el modelo de transformación según el caso requerido y determinando los siete parámetros de transformación. Sin embargo, cuando se relacionan dos marcos geodésicos derivados de mediciones satelitales como GNSS, la definición de la escala es consistente con unidad de longitud; es decir, el metro que está definido con base en la velocidad de la luz. Lo anterior implica que no hay variación temporal para la relación de escalas y, por lo tanto, una transformación de similitud no debe incluir el parámetro de factor de escala [5]. En este trabajo se aplicó el modelo de transformación baricéntrico; se excluyó el factor de escala del conjunto de parámetros.

3. METODOLOGÍA

A. La referencia geodésica de Costa Rica y las diferentes soluciones ITRF

La implementación de CR05 se hizo con vínculo al ITRF2000, específicamente en la semana 1342, mientras que la oficialización de CR-SIRGAS, como el actual marco geodésico nacional de referencia, se redujo a la semana 1803. Entre estas dos referencias nacionales se tiene un total de 8,76 años, periodo en el cual se publicaron y adoptaron globalmente diferentes soluciones del ITRF, específicamente la solución ITRF2005/IGS2005 en la semana 1400 [27], la solución ITRF2008/IGS2008 en la semana 1632 [28], la versión IGb2008 en la semana 1709 hasta la actual solución ITRF/IGS2014 a partir de la semana 1934 [29].

En la Figure 3 se muestra la situación anterior mediante una línea de tiempo, en la cual, adicionalmente, en la parte superior se presenta el establecimiento del CR05 a finales del año 2005, la adopción de CR-SIRGAS a mediados del año 2014 y la reciente actualización del CR-SIRGAS al marco ITRF/IGb2014. Además, se incluye el inicio de la más reciente solución ITRF/IGS2014 a inicio del año 2017 y la solución IGb14.

B. De CR05 a CR-SIRGAS

Uno de los objetivos fundamentales de este trabajo consistió en determinar un conjunto de parámetros de transformación que permitan llevar la información geoespacial vinculada al anterior marco CR05, al nuevo y oficial CR-SIRGAS, procurando mantener campos puntuales con la menor deformación posible. En cada una de las determinaciones se hizo un análisis de los residuales en sus correspondientes coordenadas topocéntricas [N, E, U]. Por lo tanto, se determinaron rangos de tolerancia para el análisis de los residuales de los procesos de transformación, con la finalidad de identificar los eventuales puntos inconsistentes y, además, mantener mayoritariamente la semejanza entre los marcos a relacionar. Dichos rangos permitieron encontrar incongruencias que hubieran incidido en la determinación de los diferentes conjuntos de parámetros de transformación. En la componente horizontal se establecieron como límites máximos de tolerancia ± 100 mm y ± 200 mm, usando para el límite inferior el criterio de tres veces el error promedio que se obtuvo en las coordenadas geocéntricas de CR05, el cual fue procesado con un software comercial; mientras que el límite superior se determinó basado en el error máximo de la actual cartografía urbana nacional a escala 1:1000. En la componente vertical se definió un rango de tolerancia entre ± 200 mm y ± 300 mm considerando que esta componente, en general, tiene una incertidumbre mayor a la horizontal en un factor de dos y tres veces.

4. RESULTADOS

A. Transformación directa

La determinación directa del conjunto de seis parámetros de transformación entre CR05 y CR-SIRGAS se hizo considerando las coordenadas cartesianas geocéntricas de un total de 34 puntos idénticos entre estaciones activas y puntos pasivos. En la Figure 4 se muestra la ubicación de los puntos idénticos dentro del territorial continental de Costa Rica, junto con un gráfico que muestra una primera cuantificación de las diferencias de coordenadas topocéntricas entre las dos épocas, donde se nota que la mayoría de los puntos se movieron en dirección noreste.

El procedimiento tomó las épocas de referencia de los dos marcos a relacionar, 2005,83 y 2014,59 respectivamente. El proceso general de cálculo de los parámetros de transformación no contempla los cambios o eventos que se hayan dado entre esos dos instantes de tiempo, ya que las coordenadas del conjunto de puntos idénticos simplemente representan dos realidades consideradas en distintos momentos, lo cual no permite incorporar dentro en los resultados las deformaciones que hayan sufrido los puntos, implicando por lo tanto una eventual falta de semejanza entre los marcos a relacionar. Sin embargo, en aras de ofrecer una solución práctica los residuales de las transformaciones fueron analizados con los criterios de tolerancia establecidos. Es así, como el conjunto de parámetros permite relacionar dichas épocas de referencia asumiendo adicionalmente que los dos marcos son semejantes desde el punto de vista matemático y no desde la realidad. Las dos variantes de la transformación directa entre el anterior marco CR05 y el actual CR-SIRGAS se describen a continuación:

• Determinación del conjunto de parámetros de transformación con base en las coordenadas geocéntricas de los 34 puntos idénticos activos y pasivos y contemplando una tolerancia máxima en el análisis de los residuales de ± 200 mm en las dos componentes horizontales y de ± 300 mm en la componente vertical. El resultado final excluyó 3 puntos; la estación activa NICY ubicada en la ciudad de Nicoya y dos puntos pasivos, LUCÍA en el noroeste del país y LAUREL en el sector suroeste (ver Figure 4). Los resultados de este conjunto de parámetros se presentan en la Table 1.

• Determinación del conjunto de parámetros de transformación con base en las coordenadas geocéntricas de los 34 puntos idénticos activos y pasivos y contemplando una tolerancia máxima en el análisis de los residuales de ± 100 mm en las dos componentes horizontales y de ± 200 mm en la componente vertical. El resultado final excluyó un total de 12 puntos idénticos, en su mayoría ubicados en la zona de Península de Nicoya y el Pacífico Central y Sur de Costa Rica (ver Figure 5 B). Los resultados de este conjunto de parámetros se presentan en la Table 2.

Los valores de los residuales horizontales y verticales derivados de los dos procesos de determinación de los parámetros de transformación se presentan en el conjunto de gráficos de barras de la Figure 6. En ellos, los dos superiores corresponden con los resultados contemplando 31 puntos idénticos, mientras que los dos inferiores son para los resultados usando 22 puntos idénticos.

Para conocer el comportamiento de los residuales horizontales y verticales en el conjunto de mapas de la Figure 7 se muestran los resultados de una interpolación de estos residuales, separados igualmente en los procesos de transformación directa usando 31 y 22 puntos idénticos respectivamente.

B. Transformación en escalera

Anteriormente, se expuso que desde la puesta en marcha del marco CR05 hasta el actual CR-SIRGAS se publicaron mundialmente diferentes soluciones del ITRF, las cuales no se implementaron oficialmente dentro de la referencia geodésica nacional. Sin embargo, los usuarios en las etapas de procesamientos de las observaciones GNSS, integraron estas diferentes soluciones ITRF mediante el uso de los archivos de órbitas precisas y de coordenadas finales SIRGAS, obteniendo sus resultados en un marco diferente al oficial. Además de este aspecto geométrico, desde el punto de vista físico el sismo de Nicoya en agosto de 2012 provocó importantes afectaciones y deformaciones en el marco CR05 las cuales deben ser consideradas para un proceso de transformación válido para las diferentes épocas de referencia.

A este proceso se le denominó transformación en escalera ya que, por medio de los diferentes conjuntos de parámetros de transformación, se puede llevar la información desde el anterior marco CR05 a diferentes referencias según sea la necesidad. La transformación en escalera se aplicó en dos variantes y se usaron como insumos exclusivamente las coordenadas geocéntricas de 9 estaciones GNSS de operación continua, ocho estaciones administradas por el IGNCR denominadas CIQU, LIBE, LIMN, NICO, NEIL, PUNT, RIDC y SAGE más la estación ETCG. En el análisis de los residuales se establecieron valores de tolerancia menores a los empleados en el método directo, definiéndose un rango entre los ± 30 mm y ± 50 mm en horizontal y entre los ± 50 mm y ± 100 mm en vertical, debido a que exclusivamente se usaron estaciones activas en la determinación de los parámetros de transformación y sus procesamientos se hicieron de forma rigurosa con programas científicos. Para la transformación en escalera se realizaron dos variantes, las cuales se describen a continuación:

a. Parámetros de transformación considerando como origen el anterior marco geodésico nacional de referencia CR05, época 2005,83. La situación se presenta en la Figure 8, donde se esquematiza por medio de líneas las épocas destino para un total de 5 conjuntos de parámetros de transformación. La solución internacional IGS05 inició en mayo de 2006, específicamente en la semana 1400. Sin embargo, para este estudio se hace la indicación a la semana 1591, la cual corresponde con la primera determinación de las coordenadas geocéntricas de las ocho estaciones GNSS de medición continua administradas por el IGNCR. El proceso a CR-SIRGAS ya fue contemplado en las dos variantes del apartado anterior y finalmente se indica la semana de inicio de la actual solución ITRF/IGS2014.

En este proceso de transformación en escalera se estableció una tolerancia para el análisis de los residuales un valor de ± 50 mm en las dos componentes horizontales y de ± 100 mm en la componente vertical. En la Table 3 se presentan los valores y exactitudes obtenidos de este proceso.

En el cálculo de los parámetros de transformación relativos a la semana 1705 se usaron los mismos valores de tolerancia horizontal y vertical; sin embargo, luego del análisis de los respectivos residuales se tuvo que descartar un total de seis estaciones, lo que imposibilitó seguir con el procedimiento. Para continuar con este cálculo, se duplicaron los valores iniciales de la tolerancia pasando a ± 100 mm en las dos componentes horizontales y a ± 200 mm en la componente vertical para las semanas 1705 y 1709. Este aumento de tolerancia permitió obtener los respectivos conjuntos de parámetros de transformación, cuyos valores se presentan en la Table 4.

De las 9 estaciones GNSS consideradas inicialmente en el proceso de transformación en escalera, el análisis de los residuales individuales de cada uno de los cálculos implicó que la estación CIQU fuera excluida para las semanas 1591, 1632 y 1703 y la estación NICY excluida para la semana 1709. Adicionalmente, no se contó con datos de la estación LIMN para la semana 1703 y 1705, generando que el proceso de transformación se efectuara contemplando un total de 8 puntos idénticos para las semanas 1591, 1632 y 1709, mientras que para las semanas 1703 y 1705 se usaron solamente 7 estaciones comunes. Los valores de los residuales derivados de esta primera variante de la transformación en escalera se presentan en la Table 5.

b. Parámetros de transformación considerando como origen y destino épocas de referencia continuas en las cuales hubo un cambio o actualización del ITRF. La situación se presenta en la Figure 9, donde se esquematiza por medio de líneas las épocas origen y destino para un total de 5 conjuntos de parámetros de transformación. Este proceso inició en la semana 1591 y se incluye adicionalmente un conjunto de parámetros de transformación entre las semanas 1703 y 1705, una semana antes y una semana después del terremoto de agosto de 2012 cuyas coordenadas geocéntricas ajustadas se obtuvieron de [30].

En el proceso de cálculo de los parámetros de transformación en escalera con épocas de referencia comunes, los residuales de cada determinación se analizaron con base en una tolerancia de ± 30 mm en las dos componentes horizontales y de ± 50 mm en la componente vertical. Aunque en estos cálculos se inició con un grupo de 9 estaciones GNSS, para la semana 1703 se excluyeron las estaciones CIQU y LIBE, mientras que para la semana 1709 adicionalmente se excluyó la estación NICY. Además, para este conjunto de parámetros de transformación, se tuvo que aumentar la tolerancia horizontal y vertical pasando a ± 75 mm y ± 100 mm respectivamente.

Entre las semanas 1709 y 1803 solo se pudo contar con un total de 7 estaciones comunes. Los resultados y exactitudes de estos cuatro conjuntos de parámetros de transformación se presentan en la Table 6, mientras que en la Table 7 se tiene el conjunto de parámetros entre las semanas 1703 y 1705 en la cual se trabajó con una tolerancia de ± 100 mm en las dos componentes horizontales y de ± 200 mm en la componente vertical. Estos valores fueron aumentados debido a la afectación del terremoto de agosto de 2012.

En la Table 8 se presentan los valores de los residuales obtenidos del proceso de la segunda variante de la transformación en escalera.

C. Parámetros de transformación en la zona de Nicoya

Los resultados derivados de la transformación directa en su versión de 22 puntos idénticos condujeron al cálculo de un conjunto de parámetros de transformación con base en los 12 puntos rechazados y cuya ubicación se ve en la Figure 5B. Se definieron dos zonas, denominadas Nicoya con un total de 8 puntos (LIBE, TAMARINDO 2, NICY, LUCIA, MONTEVERDE, PUNT y COLES) y la Pacífico Sur con un total 4 puntos (DOMINICAL, MAIZ, SIRENA y LAUREL). Es evidente que para la zona Pacífico Sur se tiene una configuración muy pobre en cuanto a cantidad de puntos, por lo que no se determinaron los parámetros respectivos. En el caso de la zona Nicoya se inició el cálculo con tolerancias de ± 100 mm en la parte horizontal y ± 200 mm en la vertical; como resultado, se descartaron 4 puntos, dejando el cálculo solo con 4 puntos idénticos (BLANCA, COLES, PUNT y MONTEVERDE). Cuando se usó el segundo criterio de tolerancia, es decir ± 200 mm en horizontal y ± 300 mm en vertical, se eliminó solamente la estación NICY. En la Table 9 se expone el conjunto de parámetros resultantes.

D. Consistencia de los parámetros de transformación

Los dos conjuntos de seis parámetros de transformación directa entre el anterior marco CR05 y el actual marco CR-SIRGAS presentados en las Table 1 y Table 2 fueron usados para calcular las coordenadas transformadas de un conjunto de 8 estaciones GNSS de operación continua: ETCG, LIBE, LIMN, NEIL, NICY, PUNT, RIDC y SAGE. Estas estaciones fueron calculadas diariamente por los Centros de Procesamiento SIRGAS e incluidas dentro de las soluciones semanales finales precisamente en la semana 1803; de hecho, excluyendo a la estación ETCG, las restantes son las estaciones GNSS que definen el actual marco nacional geodésico de referencia CR-SIRGAS [8].

La evaluación de la consistencia se hizo comparando directamente de las ternas de coordenadas geodésicas transformadas con los dos conjuntos de parámetros y los resultados de las coordenadas semanales finales de SIRGAS entre las semanas 1804 y 1933 [7], es decir, una semana después de la implementación de CR-SIRGAS hasta una semana previa a la implementación mundial del ITRF2014, para un total de 129 semanas. En cada una de las estaciones se hizo la comparación de las coordenadas geocéntricas y las diferencias respectivas, que se transformaron a coordenadas topocéntricas y cuyo comportamiento para las estaciones ETCG, LIBE, LIMN y NEIL se presenta en el conjunto de gráficos de la Figure 10; mientras que los resultados para las estaciones NICY, PUNT, RIDC y SAGE en muestran en la Figure 11. En ambas figuras, la línea vertical de color negro separa los dos resultados; a la izquierda la comparación con los parámetros obtenidos con el conjunto de parámetros determinados usando 31 puntos idénticos; y a la derecha, la comparación con el conjunto de parámetros de transformación obtenidos contemplando 22 puntos idénticos.

Los conjuntos de parámetros obtenidos de la variante de transformación en escalera entre épocas continuas fueron usados para calcular las coordenadas de los puntos 26 puntos pasivos. El proceso consistió en aplicar los diferentes conjuntos de parámetros de transformación desde la época 2005,83 hasta la época 2014,59 comparando con las coordenadas geocéntricas oficiales de CR-SIRGAS. Una interpolación con las diferencias en coordenadas topocéntricas se presentan en la Figure 12.

5. DISCUSIÓN Y CONCLUSIONES

• En el país se tiene una gran cantidad de información geoespacial vinculada al anterior marco geodésico de referencia CR05, la cual, dependiendo de las aplicaciones, requiere de un proceso correcto para su utilización dentro del actual y oficial marco geodésico de referencia CR-SIRGAS. Todos los conjuntos de 6 parámetros de transformación presentados abarcan exclusivamente el período entre 2005,83 y 2014,59. Un campo puntual cuyas coordenadas se obtengan como producto de la aplicación de un conjunto de parámetros de transformación, nunca podrá sustituir a una medición y determinación directa.

• La determinación de un conjunto de parámetros de transformación depende exclusivamente de la calidad de las coordenadas del conjunto de puntos idénticos y de su distribución espacial. El modelo de transformación Helmert ya sea la variante geocéntrica o baricéntrica contempla adicionalmente mantener la similitud entre los dos campos relacionados. En el caso específico de los conjuntos de 6 parámetros de transformación determinados en este trabajo, se partió de las coordenadas geocéntricas del marco CR05 cuyas exactitudes fueron aproximadamente del orden centimétrico en comparación con exactitudes de orden subcentimétricas en las determinaciones para las semanas 1400, 1632, 1709 y 1803, implicando una falta de homogeneidad en la calidad de los puntos idénticos relacionados. Sin embargo, los resultados presentados deben tomarse como conjuntos de 6 parámetros de transformación que permitirán realizar actualizaciones de coordenadas con un orden de error menor a 20 cm (ver conjuntos de mapas de la Figure 7), dependiendo del instante en el que se capturó la información y la referencia a la cual se desea llevar esa información, generalmente al actual CR-SIRGAS.

• Los valores obtenidos en todos los conjuntos de parámetros de transformación están totalmente correlacionados con las posiciones de los vértices en el momento de realizar las mediciones. La cuantificación inicial de las diferencias de coordenadas entre los 34 puntos idénticos presentó rangos de variación de 62 cm y 38 cm en las componentes norte y este respectivamente, y 25 cm en la componente vertical. En la Figure 4 se muestra que prácticamente todos los puntos presentaron un desplazamiento en sentido noreste, inclusive se nota un mayor cambio en la componente norte de todos los puntos; exceptuando algunos vértices ubicados en el sector occidental del país específicamente los puntos pasivos BLANCA, GRANDE, LUCÍA, TAMARINDO 2 y las estaciones activas LIBE y NICY, los cuales están dentro de la zona de afectación del terremoto de agosto de 2012, el cual se localizó apenas a 24 Km al sur-suroeste de playa Sámara, frente a la Península de Nicoya y a una profundidad de 15,4 Km [1]. Lo anterior implica que todos los conjuntos de 6 parámetros de transformación calculados en esta investigación solo ofrecen una aproximación a las distintas realidades geodésicas y físicas específicas que se tuvo entre mediados del año 2005 y mediados del año 2014.

• Se brindan fundamentalmente dos grupos de 6 parámetros de transformación estáticos entre el anterior marco CR05 y el actual CR-SIRGAS. El primer conjunto calculado directamente entre las épocas de referencia 2005,83 y 2014,59 partiendo de 34 puntos idénticos entre estaciones pasivas y activas, el cual a su vez fue desarrollado en dos variantes obteniendo como resultado dos campos puntuales finales de 31 y 22 vértices respectivamente. Y un segundo conjunto calculado considerando fundamentalmente los cambios de ITRF y afectación por el terremoto de agosto de 2012, denominada internamente como transformación en escalera. Este último se determinó igualmente en dos variantes con la diferencia que exclusivamente se consideraron estaciones activas GNSS, debido a que sobre los puntos pasivos no se realizaron mediciones en las épocas respectivas. En la primera variante calculó parámetros desde la época 2005,83 a las diferentes referencias; y la segunda variante, calculó los conjuntos de parámetros entre las épocas de referencia comunes.

• Los parámetros de transformación directa presentados en las Table 1 y Table 2 muestran traslaciones significativas con valores promedio de +11,0 cm ± 1 cm en X, de +15 cm ± 1 cm en Y, +12 cm ±1 cm en Z, mientras que las rotaciones son significativas para el eje Y, el eje Z y no significativa para la rotación en el eje X.

• El comportamiento de los residuales del conjunto de parámetros determinados con 31 puntos idénticos, presentaron valores en la componente horizontal entre 5 mm y 155 mm con la mayor diferencia localizada en el sector sur de la Península de Nicoya, en la parte del Pacífico central y sur y en la zona del Caribe sur (ver Figure 6A). En la componente vertical las diferencias están entre -50 mm y 100 mm con un comportamiento mayoritariamente entre -25 mm y 25 mm y con algunas zonas en el occidente, centro oriente y sur del país con valores superiores a los 25 mm (ver Figure 6B).

• Los residuales del conjunto de parámetros determinados con 22 puntos idénticos presentaron valores en la componente horizontal entre 25 mm y 115 mm con un comportamiento casi que simétrico en la zona norte del país con valores entre 40 mm y 70 mm aproximadamente, mientras que el sector sur los resultados están entre 70 mm y 100 mm (ver Figure 6C). La componente vertical presentó valores entre -80 mm y 80 mm con una distribución heterogénea con resultados, entre -40 mm y 40 mm (ver Figure 6D).

• En el conjunto de gráficos presentados en la Figure 10 y Figure 11 se muestran la consistencia en la aplicación los dos conjuntos de parámetros derivados de la transformación directa. En estos gráficos se evidencia que para el período comprendido entre la semana 1804 (mediados de del año 2014) y la semana 1933 (inicio del año 20147) las diferencias respecto a las soluciones semanales final de SIRGAS llegan a valores absolutos de 10 cm para la estación ETCG, de 12 cm para las estaciones PUNT y RIDC, de 14 cm para LIMN y NEIL, de 20 cm para la estación LIBE, de 25 cm para la estación SAGE y de 50 cm para la estación NICY. Dichos valores demuestran que la aplicación perpetua de parámetros de transformación estáticos no es correcta, ya que, por ejemplo, en un período de 129 semanas las diferencias llegan a valores de medio metro. Por lo tanto se evidencian dos alternativas: una, el cálculo de un nuevo conjunto de parámetros como los publicados en [31] y dos, la potencial determinación de un conjunto dinámico de parámetros. Sin embargo, lo correcto es usar en los procedimientos, las observaciones reales reducidas a la época de observación y las respectivas coordenadas finales de SIRGAS.

• Los conjuntos de parámetros de la transformación en escalera (variante a) brindaron resultados prácticamente significativos (ver Table 3 y Table 4). Por su parte, los residuales presentaron un comportamiento en la componente norte entre -15 mm y +18 mm, entre -38 mm y +20 mm en el componente este y entre los -16 mm y +7 mm en la componente vertical hasta la semana 1703, antes del terremoto de agosto de 2012. Luego de esta afectación, los residuales aumentaron significativamente llegando a extremos entre -96 mm y +75 mm en el norte, -61 mm y 76 mm en el este y -44 mm y +40 mm en altura (ver Table 5). Por su parte, los conjuntos de parámetros de la transformación (variante b) brindaron igualmente resultados significativos (ver Table 6 y Table 7), mientras que los residuales presentaron un comportamiento relativamente menor respecto a la variante a; presentando específicamente en la componente norte valores entre -82 mm y 35 mm, mientras que para la componente este se tuvo resultados entre -61 mm y 48 mm y en la componente vertical entre -20 mm y 11 mm (ver Table 8). Se tiene igualmente una consistencia en los resultados al comparar con los cambios sufridos por los puntos entre las dos épocas de referencia, especialmente en la componente norte.

• Los conjuntos de parámetros de la transformación en escalera (variante b) permitieron detectar diferencias en valor absoluto sobre los puntos pasivos de la red nacional de primer orden de 30 cm en la componente horizontal y vertical. En la Figure 12A y Figure 12B, se muestra como las mayores diferencias se presentan occidental del país, específicamente en el vértice LUCÍA y en el sector sur en el punto LAUREL; mientras que en vertical las mayores discrepancias están en sector occidental de Costa Rica comprendido entre los 10° N y 11° N y -85° W y -86° W principalmente.

• Los resultados mostraron como los diferentes criterios de tolerancia para el cálculo de los parámetros de transformación permitieron detectar que prácticamente todos los vértices idénticos del sector occidental del país sufrieron cambios de más de ±100 mm. Si específicamente se considera la exactitud relativa de ± 20 cm dada por la cartografía nacional urbana a escala 1:1000, como el criterio para la aplicación de un determinado grupo de parámetros de transformación, quedará nuevamente en los usuarios la aplicación de un análisis previo del comportamiento de los residuales de cada conjunto.

• Los puntos excluidos de la transformación directa entre CR05 y CR-SIRGAS ubicados en el sector occidental de Costa Rica, se usaron para el cálculo de un conjunto de parámetros de transformación en la zona denominada Nicoya. Los resultados presentados en la Table 9 muestran una traslación significativa en dirección Y de +10 cm ±3 cm, así como rotaciones significativas en los tres ejes. Ante la poca cantidad de puntos idénticos en la zona denominada Pacífico Sur no se calculó el conjunto de parámetros respectivos y eventualmente integrar como refuerzo de esta zona por ejemplo el punto pasivo ALEGRE y las estaciones activas SAGE y NEIL (ver Figure 5B) implicaría una combinación de exactitudes relativas que provocaría su exclusión del proceso de transformación. De la misma manera una determinación de parámetros de transformación considerando estos 12 puntos presenta una configuración inadecuada, recomendándose en su lugar usar la variante directa con 31 puntos idénticos. Se recalca el hecho que este conjunto de parámetros de transformación se brinda para que el usuario disponga de información que pueda eventualmente utilizar y sobre todo analizar los resultados de su aplicación como un insumo para una toma adecuada de decisiones.

• Es por lo anteriormente expuesto, que el usuario final debe conocer los metadatos de su información origen (fundamentalmente fecha de captura, datum asociado, elipsoide de referencia y la referencia de las estaciones usadas como vínculo), con el objetivo de que se aplique el conjunto de parámetros correspondientes que le permita transformar la información origen al destino requerido. La intención de los autores es, por lo tanto, ofrecer a los distintos usuarios grupos de parámetros de transformación que, de acuerdo con la referencia de la información origen y sus necesidades específicas, le ayuden a una adecuada y correcta georreferenciación; siguiendo con lo especificado en [8]. Además, se le da también información sobre las potenciales diferencias en los resultados que se puedan obtener dependiendo de cuál conjunto de parámetros sea aplicado y en qué zona geográfica del país.

• La utilización de CR-SIRGAS como marco nacional geodésico de referencia en Costa Rica es de carácter obligatorio desde el año 2018, para todas las actividades técnicas que involucren tanto la captura como el procesamiento de información geoespacial nacional. Durante la redacción de este trabajo, el IGNCR publicó la primera actualización de CR-SIRGAS en cuanto a época de referencia y marco internacional, siendo respectivamente época 2019,24 y marco ITRF/IGS2014. Además, se ofrecen un conjunto de 7 parámetros de transformación entre las épocas 2014,59 y 2019,24 [31]. De esta manera, se tiene actualmente la posibilidad de transformar información vinculada a CR05 hasta la primera solución de CR-SIRGAS, semana de referencia 1803 en las dos variantes directas y las dos variantes en escalera. Y de esta semana, hasta la semana 2046 correspondiente con la época 2019,24.

• Finalmente, es importante que los usuarios nacionales estén relacionados con los cambios en el marco internacional como la próxima solución ITRF2020 [32] prevista para 2021, la cual implicará una nueva y necesaria actualización en el marco geodésico nacional.

REFERENCIAS

L. Linkimer et al., “El terremoto de Sámara (Costa Rica) del 5 de setiembre del 2012 (Mw 7,6)”, Rev. Geológica América Cent., no. 49, Jan. 2011, doi: 10.15517/rgac.v0i49.13104.

L. Sánchez y H. Drewes, “Geodetic Monitoring of the Variable Surface Deformation in Latin America”, in International Association of Geodesy Symposia, Springer, Berlin, Heidelberg, 2020.

M. S. Molodensky, Methods for study of the external gravitational field and figure of the earth. 1962.

J. Badekas, Investigations related to the establishment of a world geodetic system, Ohio State University, 1962.

H. Drewes, “How to Fix the Geodetic Datum for Reference Frames in Geosciences Applications?”, En Geodesy for Planet Earth, International Association of Geodesy Symposia, Vol 136., S. Kenyon,M. Pacino, y U. Marti, Eds. 2012, pp. 67–76.

Decreto Ejecutivo Número 33979-MJ-MOPT, “Declara como datum horizontal oficial para Costa Rica el CR05”, La Gaceta No 108, Jun. 06, 2007

SIRGAS, “Coordenadas semanales de las estaciones SIRGAS-CON”, 2020.http://www.sirgas.org/es/sirgas-con-network/coordinates/weekly-positions/ (accesado en Ago. 10, 2020)

Decreto Ejecutivo Número 40962-MJP, “Actualizacion del Sistema Geodésico de Referencia Horizontal Oficial para Costa Rica”, La Gaceta No 66, Apr. 17, 2018.

P. Wessel et al., “The Generic Mapping Tools Version 6”, Geochemistry, Geophys. Geosystems, vol. 20, no. 11, pp. 5556–5564, Nov. 2019, doi: 10.1029/2019GC008515.

D. McCarthy, “IERS Technical Note 21. IERS Conventions (1992)”, Paris, 1996. [en línea]. Disponible:https://www.iers.org/IERS/EN/Publications/TechnicalNotes/tn21.html?nn=94912

D. McCarthy and G. Petit, “IERS Conventions (2003). IERS Technical Note 32. Verlag des Bundesamts für Kartographie und Geodäsie, Frankfurt am Main 2004”, Paris, 2004. [en línea]. Disponible:https://www.iers.org/IERS/EN/Publications/TechnicalNotes/tn32.html?nn=94912

A. Leick, GPS Satellite Surveying, 4th ed. New Jersey: Jhon Wiley and Sons, 2015.

G. Xu and Y. Xu, GPS Theory, Algorithms and Applications, 3rd ed. Berlin, Heidelberg: Springer Berlin Heidelberg, 2016.

W. Torge, Geodesy, 3rd ed. Berlin: de Gruyter, 2001.

M. Bursa, “The theory for the determination of the non-parallelism of the minor axis of the referenceellipsoid and the inertial polar axis of the earth, and the planes of the initial astronomic and geodeticmeridians from observations of artificial earth satellites”, Stud. Geophys. Geod., vol. 6, pp. 209–214,1962.

H. Wolf, “Geometric connection and re-orientation of three-dimensional triangulation nets”, Bull.géodésique, vol. 68, no. 1, pp. 165–169, Jun. 1963, doi: 10.1007/BF02526150.

B. Hofmann-Wellenhof, H. Lichtenegger, y J. Collins, GPS Theory and Practice, 5th ed. Wien: Springer-Verlag Wein New York, 2001.

B. Hofmann-Wellenhof, H. Lichtenegger, y E. Wasle, GNSS - Global Navigation Satellite System GPS, GLONASS, Galileo and more. Wein: Springer Wein New York, 2008.

M. Hooijiberg, Geometrical Geodesy. Berlin, Heidelberg: Springer Berlin Heidelberg, 2008.

J. Subirana, J. M. J. Zornoza, y M. Hernández Pajares, GNSS Data Processing Volume I: Fundamentals and Algorithms (ESA TM-23/1; May 2013). ESA Communications, 2013.

C. Jekeli y O. Montenbruck, “Time and Reference Systems”, En Springer Handbook of Global Navigation Satellite Systems, Cham: Springer International Publishing, 2017, pp. 25–58.

G. Xu y J. Xu, “Coordinate and Time Systems”, in Orbits, Berlin, Heidelberg: Springer Berlin Heidelberg, 2013, pp. 5–26.

C. Boucher y Altamimi, “IERS Technical Note 6. ITRF 89 and other realizations of the IERS Terrestrial Referencia System for 1989”, Paris, 1989. [en línea]. Disponible: https://www.iers.org/IERS/EN/Publications/TechnicalNotes/tn06.html?nn=94912

G. Petit y B. Luzum, “IERS Conventions (2010). IERS Technical Note 36. Verlag des Bundesamts für Kartographie und Geodäsie, Frankfurt am Main 2010”, Paris, 2010. [en línea]. Disponible:https:// www.iers.org/IERS/EN/Publications/TechnicalNotes/tn36.html?nn=94912

Z. Altamimi, P. Rebischung, L. Métivier, y X. Collilieux, “ITRF2014: A new release of the International Terrestrial Reference Frame modeling nonlinear station motions”, J. Geophys. Res. Solid Earth, vol. 121, no. 8, pp. 6109–6131, Aug. 2016, doi: 10.1002/2016JB013098.

Z. Altamimi, P. Rebischung, L. Métivier, y X. Collilieux, “ITRF2014: A new release of the International Terrestrial Reference Frame modeling nonlinear station motions”, J. Geophys. Res. Solid Earth, vol. 121, no. 8, pp. 6109–6131, Aug. 2016, doi: 10.1002/2016JB013098.

T. Soler y J. Marshall, “A note on frame transformations with applications to geodetic datums”, GPS Solut., vol. 7, no. 1, pp. 23–32, May 2003, doi: 10.1007/s10291-003-0044-8.

Z. Altamimi, X. Collilieux, J. Legrand, B. Garayt, y C. Boucher, “ITRF2005: A new release of the International Terrestrial Reference Frame based on time series of station positions and Earth Orientation Parameters”, J. Geophys. Res., vol. 112, no. B9, p. B09401, Sep. 2007, doi: 10.1029/2007JB004949.

Z. Altamimi, X. Collilieux, y L. Métivier, “ITRF2008: an improved solution of the international terrestrial reference frame”, J. Geod., vol. 85, no. 8, pp. 457–473, Aug. 2011, doi: 10.1007/s00190-011-0444-4.

Y. Céspedes, “Análisis de la influencia del marco de referencia SIRGAS en las soluciones semanales semilibres de la red de estaciones GNSS del Registro Nacional”, Universidad Nacional, Heredia, Costa Rica, 2016.

Instituto Geográfico Nacional, “Directriz Nº DIG-001-2020 Parámetros de transformación para pasar de las épocas 2014,59 a la 2019,24 en el ITRF14 correspondiente con CR-SIRGAS”, La Gaceta No 223, Sep. 04, 2020

Z. Altamimi, X. Collilieux, P. Rebischung, L. Metivier, y K. Chanard, “The ITRF Three decades of research and development, usages and applications”, presentado en Implementation of the UN-GA Resolution on the GGRF in Latin America, Argenitna, 16-20 setiembre, 2019.