Type your text Compendio de Trabajos de Investigación CNDG – Biblioteca Instituto Geofísico del Perú. V. 6 (2005) p. 9 - 18 RELACIÓN ESPACIO-TEMPORAL ENTRE ÍNDICES HIDROLÓGICOS EN EL PERÚ Y VARIABLES CLIMÁTICAS GLOBALES JOE TIPIANI MONTES E.A.P. Ingeniería Mecánica de Fluidos Universidad Nacional Mayor de San Marcos joet@huascaran.igp.gob.pe, jtipianim@yahoo.com Prácticas dirigidas por: Ing. Grace Trasmonte Soto Centro de Predicción Numérica del Tiempo y Clima RESUMEN En el presente trabajo se analizó si existe una posible relación espacio-temporal, mediante el uso de mapas de correlación lineal (r), entre dos series hidrológicas ubicadas en el norte del Perú (descargas del río Piura) y el Altiplano (nivel del lago Titicaca), con algunas variables climáticas a escala global (Temperatura Superficial del Mar -TSM, Presión superficial a Nivel del Mar -PNM y la Radiación de Onda Larga –ROL, que es un índice de lluvias). Los principales resultados sugieren que las descargas (por ende las lluvias) en Piura, durante el periodo de lluvias (enero a abril) están altamente correlacionadas (r =+/- 0.7 a 0.8) con sectores del Pacífico ecuatorial en general, en forma inversa con el sector occidental en comparación a los sectores oriental y central; en el caso de los sectores central y oriental, particularmente alrededor de la longitud 140° oeste y entre 80° y 100° de longitud oeste, se obtuvo una relación directa con las anomalías de TSM, inversa con las anomalías de PNM y directa/inversa con las anomalías de lluvias/ROL; el sector occidental (principalmente entre 140 y 160 de longitud este), por el contrario presentó una relación inversa con las anomalías de TSM, directa con las de PNM e inversa/directa con las de las lluvias/ROL. Por otro lado, si bien las correlaciones que se han encontrado entre los niveles del lago Titicaca (por ende también las lluvias en el Altiplano peruano) y las variables climáticas globales, en el mejor de los casos son regulares (r =+/- 0.4 a 0.5), podría haber una relación con el sector del Atlántico norte (especialmente el sector entre 0° y 20° de latitud norte y 40° y 60° de longitud oeste) durante el verano y otoño, inversa con la TSM, directa con la PNM e inversa/directa con las lluvias/ROL, así también se encontró una relación inversa (probablemente por efecto de teleconexión), con las lluvias de sectores subtropicales en el Pacífico central sur, entre invierno y primavera (promedio de julio a setiembre y de agosto a setiembre). INTRODUCCIÓN La disponibilidad del recurso agua es de interés general, por su escasa disponibilidad espacial, alta variabilidad de su ciclo e importancia en el quehacer humano. Su alta variabilidad puede ser explicado por cambios o mecanismos de la atmósfera o del sistema climático océano- atmósfera-suelo, a diferentes escalas espaciales y temporales. En el presente estudio se analiza la posible influencia de factores climáticos, océano-atmosféricos (Temperatura Superficial del Mar -TSM, Presión superficial a Nivel del Mar -PNM y la Radiación de Onda Larga –ROL, que es un índice de lluvias), de diferentes áreas geográficas localizadas entre el Pacífico y el Atlántico, a la variabilidad de dos cuencas hidrográficas del Perú, como son la del río Piura, ubicada en la costa norte (departamento de Piura) y la del Lago Titicaca, ubicada en la sierra sur (departamento de Puno). El Río Piura nace en el Ecuador y desemboca en el Océano Pacífico, es un río de corta longitud, muy torrentoso, cuyo caudal se ve incrementado en el verano, J. Tipiani 10 principalmente entre los meses de diciembre a marzo. El lago Titicaca es el centro de una gran altiplanicie de unos 200,000 km2, conocida como Altiplano o meseta del Collao, situada a 3809 msnm, tiene una superficie de 8562 km2, el volumen total de agua se calcula en 903 km3, alimentado por la descarga de más de 25 ríos, el mayor de ellos es el Ramis, que drena cerca del 20% de la cuenca del Titicaca, desde el extremo noroeste, el Desaguadero, en el sur del lago, es el único río que drena al Titicaca y lo conecta con el Lago Poopó, por este río sale sólo el 5% del agua que entra al lago, por otro lado, se pierde una cantidad importante de humedad por evaporación (calculado en hasta 600m3/s), debido al intenso sol y los fuertes vientos altiplánicos; el nivel del lago varía estacionalmente, y en ciclos de varios años; durante la estación lluviosa - entre diciembre y marzo- el nivel del lago se incrementa, bajando normalmente durante los meses de invierno. ANTECEDENTES La relación entre variables hidrológicas y factores atmosféricos a escala regional global, ha sido estudiado por diferentes investigadores. Montesinos y Garreaud (1994), hallaron una asociación significativa entre las anomalías de la TSM en el Pacífico central y oriental, con el régimen hidrológico de la zona central de Chile, que permitió el desarrollo de modelos de pronostico de las condiciones pluviométricas en dicho país. Trasmonte (1998), analizando información de lluvias y descargas de la cuenca del río Chancay-Lambayeque en el departamento de Chiclayo, encontró que las lluvias en la cuenca baja y media eran altamente correlacionadas con la variabilidad océano-atmosférica en la costa norte peruana (campos de TSM en Puerto Chicama y de PNM, entre Piura, Talara y Chiclayo); dicha correlación se intensificaba durante los eventos El Niño/La Niña. Woodman (1999) encontró una correlación alta entre las lluvias de Piura y las temperaturas superficiales marinas en varios sectores del Pacífico oriental, frente a las costas del Perú, diseñando un modelo de regresión no lineal para el pronóstico de las lluvias en dicha región. Cárdenas (2001), correlacionó datos de lluvias de la cuenca del río Chancay- Huaral en la sierra central peruana, con algunos índices climáticos de El Niño (Indice de Oscilación Sur, TSM en las áreas Niño 3 y 1+2 y el Indice ROL); encontró que la TSM en las áreas Niño 1+2 y 3, es la variable mejor correlacionada con las precipitaciones en la cuenca (con valores de r del orden de 0.5). Gálvez (2002), indicó que los sistemas que estimulan la ocurrencia de lluvias intensas generalizadas en toda la cuenca del Mantaro y en la sierra central peruana, podrían estar asociadas a periodos Relación espacio-temporal entre índices hidrológicos en el Perú y variables climáticas globales 11 ligeramente fríos en el Pacífico central y occidental (Areas Niño 4 y 3.4). También encontró otras posibles áreas de influencia como sectores del Atlántico norte y sur. DATOS Y METODOLOGÍA En el presente trabajo se utilizaron datos mensuales provenientes del reanalysis del NCEP/NCAR (Kalnay, 1996) para las variables globales, estos son datos océano- atmosféricos de todo el globo terrestre, son una combinación de datos observados (de estaciones puntuales y de satélites) y otros generados mediante modelos, los cuales vienen en formato binario y en grilla. La información disponible del reanálisis puede ser para algunas horas específicas -como las horas sinópticas (00, 06, 12, 18 UTC)-, diaria o mensual. En la mayoría de variables océano-atmosféricas se cuenta con el periodo de información de enero de 1948 a la fecha, aunque hay casos como las del ROL, cuyo inicio de su periodo de información es de 1958. Para el presente estudio se eligió como periodo de análisis entre 1958 y 1998, espacialmente se va a trabajar con el área comprendida entre los meridianos 140° este y 20° este, y latitudes 60° norte y 70° sur, abarcando gran parte del continente americano y los océanos adyacentes- Pacífico y Atlántico-; el límite temporal superior fue dado por la disponibilidad de los datos del índice de radiación de onda larga ROL y el límite temporal inferior fue dado por la accesibilidad hasta dicho periodo en el banco de datos del IGP. En el caso de las variables hidrológicas locales se trabajó con el mismo periodo temporal, es decir entre enero de 1958 y diciembre de 1998. Se contó con datos del nivel del lago Titicaca (Ntiticaca), de la estación Muelle Enafer (15°50’ latitud sur, 70°01’ longitud oeste y altitud 3809.92 msnm), promedios mensuales dados en metros, información de la Dirección de Hidrografía de la Marina- DHN; así mismo, con datos de descargas del río Piura (Qpiura), de la estación Puente Sánchez Cerro (05°11’55” latitud sur, 80°37’20” longitud oeste y 23 msnm), promedios mensuales dados en m3/s, cuya fuente es el Ministerio de Agricultura. Inicialmente los datos fueron analizados en su calidad y consistencia, además se tuvo que realizar algunas transformaciones de formato, como de excel a formato ascii y de ascii a formato binario, en este último caso mediante el uso de un programa en fortran. Para el procesamiento estadístico y gráfico de la información, principalmente en la obtención de los mapas de correlación entre las variables globales y locales, se contó con el software conocido como GrADS (Grid Analysis and Display Sistem, 2004), con programas que fueron diseñados y elaborados por personal del CPNTC/IGP. J. Tipiani 12 La primera parte del trabajo esta dedicada a la obtención de los mapas de correlación, mediante el uso de coeficientes de correlación lineal, que son índices estadísticos, encargados de medir y representar la interacción lineal entre dos variables cuantitativas, con el objetivo de determinar el grado de asociación que existe entre estas (Calzada, 1981). Para el caso específico del presente trabajo, se obtiene la correlación entre la variable de una estación puntual (hidrológica) y el campo de datos de una variable global (variable océano-atmosférica). Los datos hidrológicos fueron correlacionados con las anomalías de los campos espaciales de TSM, PNM y ROL, en forma mensual (enero, febrero, etc.), trimestral (promedios enero a marzo, febrero a abril y así sucesivamente), y promedios de mayor número de meses (enero a abril, para el periodo lluvioso y junio a setiembre para el periodo seco), la idea de promediar varios meses permite en algunos casos mejorar la correlación, a la vez de incorporar el factor de desfase en el tiempo, de una posible influencia de los campos océano-atmosféricos. La climatología utilizada para obtener las anomalías de las variables globales fue de todo el periodo de registro de información, es decir de 1948 a 1998 para la TSM, PNM y de 1958 a 1998 para el ROL. Finalmente, se analizaron los resultados teniendo en cuenta como van a verse afectadas nuestras variables de estudio al presentarse índices de correlación significativos y según el tipo de dependencia resultante, directa o indirecta (positiva o negativa). RESULTADOS Y DISCUSIÓN En general se encontraron correlaciones altas entre las variables océano- atmosféricas y las descargas del río Piura (valor absoluto de r entre 0.7 y 0.8), entre primavera y verano y por el contrario, correlaciones moderadas, en el mejor de los casos, con el nivel del lago Titicaca (valor absoluto de r entre 0.4 a 0.5), en algunos periodos del año, no con tanta persistencia como en el caso anterior. Analizando los valores obtenidos con las anomalías de TSM (Figura 1), se pudo observar que en promedio para el período húmedo (lluvias) en el Perú, entre enero y abril, la correlación es alta y directa (positiva) entre QPiura (y por ende las lluvias en Piura) y un buen sector del Pacífico central y oriental (140° longitud oeste -80° longitud oeste, 0°- 10° latitud sur), con valores de r entre 0.6 y 0.7, ésta relación se debilita entre invierno y primavera (julio a octubre), con valores de r inferiores a +0.4. En el caso del Ntiticaca, la correlación más evidente, aunque con valores entre -0.3 y –0.4, se dio para este mismo periodo con un sector Relación espacio-temporal entre índices hidrológicos en el Perú y variables climáticas globales 13 del Atlántico norte (2° -20° latitud norte y 40°- 60° longitud oeste, al norte del continente sudamericano), dicha correlación se mantiene entre marzo y mayo. Figura 1. Campos de Correlación entre Variables Hidrológicas y Temperaturas superficiales -Promedio de enero a marzo. Los resultados con el índice de PNM, presentan por lo general una correlación inversa a la indicada por la TSM. Ratifican la importante influencia del Pacífico ecuatorial central y oriental sobre Qpiura, entre verano y otoño y su debilitamiento entre invierno y primavera, pero además descubre una relación de tipo directo (es decir a mayor presión, mayor descarga y viceversa) con sectores del Pacífico occidental (aproximadamente entre las longitudes 150° este y 160° oeste y latitudes 10°sur y 15°norte), con valores de r hasta +0.8 para el promedio de enero a abril (Figura 2, panel superior). En el caso de Ntiticaca, también se verifica la influencia del Atlántico norte en el período húmedo (entre diciembre y abril), es decir J. Tipiani 14 a mayores presiones en el Atlántico norte se posibilitaría mayores valores del nivel del lago Titicaca, aun cuando las correlaciones oscilan entre +0.3 y +0.4 (Figura 2, panel inferior). Figura 2. Campos de Correlación entre Variables Hidrológicas y Presiones al nivel de superficie - Promedio enero a abril. En el caso del ROL, las más altas correlaciones (valor absoluto de r entre 0.7 y 0.8) con Qpiura, se encuentran similarmente a los casos anteriores entre diciembre y abril, de tipo inverso con el sector central y oriental del Pacífico ecuatorial, particularmente alrededor de la longitud 140° oeste y entre 80° y 100° de longitud oeste (véase Figura 3) -es decir, a menor valor del ROL (mayor precipitación) en dichas zonas, se esperaría mayor descarga en el río Piura-, así como también se ha encontrado la relación contraria a la anterior con el Relación espacio-temporal entre índices hidrológicos en el Perú y variables climáticas globales 15 Pacífico occidental (r hasta +0.7)- sector comprendido entre 1as latitudes 0° y 10° norte y longitudes 140° y 160° este. Con el Ntiticaca, se esboza también una correlación directa con el Atlántico norte (Figura3, panel inferior), con valor entre +0.3 y +0.4 entre verano y otoño -los mayores valores se dan en los meses de enero y mayo-, esto quiere decir que lluvias mayores de lo normal en algunos o varios sectores del Atlántico norte, principalmente el sector entre 0° y 10° latitud norte y 40° a 60° longitud oeste, podrían estar relacionadas con menores lluvias y descargas en el Altiplano y viceversa. Figura 3. Campos de Correlación entre Variables Hidrológicas y Radiaciones de Onda Larga- Promedio febrero a abril. J. Tipiani 16 Adicionalmente se encontró otras correlaciones interesantes y que vale la pena mencionarlas, como el caso de una posible influencia extratropical, por teleconexión, del Pacífico central y oriental (aproximadamente entre las latitudes 50° y 65° sur y longitudes 180° - 160° oeste y 120° - 80° oeste), en las descargas del Lago Titicaca (Figura 4, panel inferior), entre invierno e inicios de primavera (promedio de junio a agosto, julio a setiembre) observado en el ROL, con valores de r hasta +0.5 (los mayores valores obtenidos con la serie del lago Titicaca), así como también correlaciones importantes entre los caudales de Piura y la Zona de Convergencia del Pacífico sur (ZCPS) acentuado durante el otoño (marzo a mayo), observados en el campo de presiones, con valores de r entre +0.5 y +0.6 (Figura 4, panel superior). El primer caso indicaría que a mayor/menor ROL/lluvias en la zona subtropical del Pacífico en las áreas indicadas, podría incrementar los niveles del lago Titicaca y viceversa; en el segundo caso, a mayores presiones en la zona de la ZCPS, podría aumentar también las descargas en el río Piura, y viceversa. CONCLUSIONES Se ha encontrado una alta correlación (valor absoluto de r entre 0.6 y 0.8) entre el Pacífico ecuatorial y las descargas del río Piura en el periodo de lluvias (diciembre y abril), disminuyendo sustantivamente dicha relación entre julio y octubre. Se distingue una influencia opuesta entre sectores del Pacífico central y oriental y el occidental, así, las descargas del río Piura estarían relacionadas directamente con las anomalías de TSM y lluvias (representados por el ROL) del Pacífico ecuatorial oriental y central, especialmente entre los sectores 80° y 100° de longitud oeste y alrededor de 140° longitud oeste, e inversamente con las anomalías de TSM y lluvias del Pacífico occidental, sobre todo el sector comprendido entre 1as latitudes 0° y 10° norte y longitudes 140° y 160° este. Las relaciones con las anomalías de PNM son contrarias a las indicadas. Las correlaciones encontradas con el nivel del lago Titicaca son de bajas a moderadas(valor absoluto de r hasta 0.5), sin embargo, se pudo definir algunas posibles áreas de influencia al período de lluvias, en el Atlántico norte, especialmente el sector entre 0° y 20° de latitud norte y 40° y 60° de longitud oeste, así como el sector subtropical del Pacífico central y oriental (aproximadamente entre las latitudes 50° y 65° sur y longitudes 180° y 160° oeste y 120° y 80° oeste) entre invierno e inicios de primavera. Los relativamente bajos valores de correlación indicarían la existencia de un mayor número de factores que intervienen en la variabilidad del sector del Altiplano, y que no están siendo tomados en cuenta en el estudio. Relación espacio-temporal entre índices hidrológicos en el Perú y variables climáticas globales 17 Figura 4. Campos de Correlación entre Variables Hidrológicas y Variables Climáticas. Promedios marzo a mayo y junio a setiembre. RECOMENDACIONES Se podrían realizar estudios similares en las áreas de estudio, utilizando variables hidrológicas complementarias, como lluvias por ejemplo, que podrían mejorar la definición de las mismas u otras posibles áreas de influencia en las características hidrológicas. Se recomienda realizar estudios de correlación espacial con desfasajes temporales, entre 2 y 12 meses (antes y después), esto puede ser utilizado para mejorar o aumentar el número de posibles variables predictoras que podrían trabajarse en un modelo de pronóstico estadístico, de regresiones múltiples. Esto sería más factible para las descargas del J. Tipiani 18 río Piura, por las mejores correlaciones que se han obtenido. Además también se podrían realizar estudios similares para otras zonas de interés del país. AGRADECIMIENTOS Agradezco al Centro de Predicción Numérica del Tiempo y del Clima del Instituto Geofísico Nacional por la oportunidad de realizar mis prácticas pre- profesionales y de manera muy especial a la Ing. Grace Trasmonte S. por su asesoramiento, guía y múltiples contribuciones y recomendaciones para la elaboración del presente estudio. También hacer participes de mis agradecimientos al Met. José Manuel Gálvez y al Lic. Raúl Chávez, por compartir sus programas los que me fueron de gran utilidad y a todas aquellas personas quienes de alguna manera hicieron posible la realización de este estudio. BIBLIOGRAFÍA Calzada, B. (1981). Métodos Estadísticos para la Investigación. Editorial Milagros S.A; Lima, Perú. 463 p. Cárdenas, J. (2001). Predicción de las Precipitaciones sobre la cuenca Chancay- Huaral. En: Tavera, H. (Ed.), Compendio de Trabajos de Investigación realizados por estudiantes, Volumen 2, Instituto Geofísico del Perú, pp. 31 – 40. Gálvez, J. (2002). Estudio de las Precipitaciones en la Cuenca del Mantaro y sus Causas. http://huascaran.igp.gob.pe/users/jose/T1- mantaro/index1.html. GrADS (Grid Analysis and Display Sistem) (2004). http://grads.iges.org/grads/grads.html Kalnay, R. (1996). The NCEP/NCAR 40- year Reanalysis Project. Bull. Amer. Meteor. Soc. 77, pp. 437-471. Montecinos, A y R. Garreaud (1994). Pronostico Estacional Del Régimen Fluviométrico en Chile Central. XVI Congreso Latinoamericano de Hidráulica. Santiago, Chile. Trasmonte, G. (1998). Régimen de Lluvias y Descargas de la Cuenca del río Chancay-Lambayeque y su relación con el Fenómeno El Niño/ La Niña. Tesis para optar el título profesional de Ingeniero Mecánico de Fluidos. UNMSM. Lima. 105 p. Woodman, R. (1999). El Fenómeno El Niño y el Clima en el Perú. En los Albores del Siglo XXI, Ciclo de conferencias 1997-1998. Ediciones del Congreso del Perú. Lima-Perú. pp. 201-242. Compendio de Trabajos de Investigación CNDG – Biblioteca Instituto Geofísico del Perú. V. 6 (2005) p. 9 - 18 Compendio de Trabajos de Investigación CNDG – Biblioteca Instituto Geofísico del Perú. V. 6 (2005) p. 19 - 28 GENERACION DE GRAFICAS TRIDIMENSIONALES PARA EL ANALISIS Y PRONOSTICO DEL TIEMPO Y CLIMA– USO DEL VIS5D JOSÉ ANTONIO INGA MONTAÑEZ Especialidad Ing. Física Facultad de Ciencias Universidad Nacional de Ingeniería jaim_111@yahoo.es, jaim103@hotmail.com Prácticas dirigidas por: Ing. Grace Trasmonte Soto Centro de Predicción Numérica del Tiempo y Clima RESUMEN En este trabajo se explica principalmente el manejo y uso del programa VIS5D (Visualización de Datos en cinco Dimensiones), el cual nos permite generar campos en tres dimensiones, siendo esto de gran utilidad en las tareas de evaluación de la circulación atmosférica. Se muestra como se manejan los datos de entrada al programa, así como también el uso de las diversas herramientas que se pueden utilizar para la visualización y tener una mejor perspectiva de la investigación que se esta realizando. INTRODUCCION Realizar el análisis del tiempo o clima consiste en observar, medir, colectar, transmitir, procesar e interpretar gran cantidad de datos que pueden abarcar incluso a todo el globo terrestre. Estos datos deben ser analizados cuidadosamente para tener una visión más clara de las condiciones de la atmósfera. Como la atmósfera cambia continuamente, este análisis por lo general, debe ser realizado en el menor tiempo posible. Mientras que la tecnología en hardware y software avanza, nuevas herramientas se están desarrollando para manipular tal cantidad de información, simular las condiciones atmosféricas actuales y las futuras, dadas por los pronósticos numéricos del tiempo o clima. En este contexto, el proceso de manipulación y visualización de toda esa información ha ido desarrollándose, contándose con herramientas cada vez más eficientes, que permiten un manejo tridimensional (3D) y hasta multidimensional de los sistemas y pronósticos atmosféricos. Una herramienta útil, para visualizar la información atmosférica, fue desarrollada por la Universidad de Wisconsin (University of Wisconsin's Space Science and Engineering Center -UWSSEC), como software libre, conocido con el nombre de VIS5D, el cual es tema del presente trabajo. Se describirán las características más importantes de este software, así como su manejo, aplicado a información para análisis y pronostico del tiempo y clima, que podría ser utilizado en las actividades del Centro de Predicción numérica del Tiempo y Clima (CPNTC) del IGP. J. Inga 20 EL VIS5D El VIS5D fue desarrollado por Bill Hibbard en el Centro de Ciencia del Espacio e Ingeniería de la Universidad de Wisconsin (University of Wisconsin's Space Science and Engineering Center UWSSEC), para visualizar ordenes dadas en cinco dimensiones, las cuales son latitud, longitud, altura, tiempo, y otro u otros campos físicos (por ejemplo: viento, presión, temperatura, humedad, etc.), es decir, en nuestro caso, información meteorológica sobre rejillas regulares y geográficamente referenciadas. Tiene herramientas simples y de fácil manejo para crear visualizaciones en 3D, lo cual se realiza usando conceptos matemáticos de cortes transversales, representación de volumen, vectores multiespaciales, isosuperficies (es decir, superficies tridimensionales que demuestran todos los puntos de referencias para una variable que tiene el mismo valor), así como otras aplicaciones. También incluye las instalaciones para las visualizaciones de las líneas proyectadas sobre mapas (limites geográficos) de países, departamentos, etc., así como de información topográfica superficial e imágenes provenientes de satélites, como por ejemplo el de cobertura de las nubes. El programa VIS5D en estos momentos es un software libre, el cual cada versión fue evolucionando de acuerdo a las sugerencias y las necesidades de los usuarios; en la actualidad hay versiones disponibles desde la 4.2 (con algunas dificultades por que no viene con algunas librerías o subprogramas que pueden ser de interés para el operador) hasta la 5.0. La versión mejorada que se puede obtener de la red, es la 5.1, que incluye más herramientas en el panel de control, lo cual hace más fácil su operación y la visualización de los datos que se están procesando. Hardware y software básico Para poder trabajar con el programa VIS5D es necesario los siguientes requisitos de equipo de computo y sistema operativo para su ejecución y operación: • Estación de trabajo, que puede ser: Silicon Graphics con IRIX 5.x o superior, IBM RS/6000 con AIX 3 o superior (apoyado por hardware OpenGL-based 3-D, Sun con SunOS 5.x o superior), HP, HP-UX A.09.01 o superior (apoyado por el hardware PEX-based 3-D) ó DEC Alpha, con OSF/1 V1.3 o superior. • Computadoras personales, IBM PC compatibles con Linux v1.2 o superior, Pentium 90MHz o con CPU superior. • Por lo menos se recomienda 32 MB de RAM en todos los casos. • Por lo menos se recomienda para el color 8-bit y si fuera posible de 24- bit. • Sistema operativo: Windows NT ó OS/2, que funcione con Intel. Generación de gráficas tridimensionales para el análisis y pronóstico del Tiempo y Clima– uso del VIS5A 21 Datos de entrada Para poder usar el VIS5D, los datos a graficarse deben ingresar en un formato especial o similar. No todos los datos atmosféricos o salidas de los modelos numéricos están en el formato VIS5D, para ello, existen diversos programas o librerías que realizan la transformación al formato VIS5D, lo cual se verá en detalle en el siguiente subcapítulo. El VIS5D puede trabajar con los datos organizados como rectángulo en cinco dimensiones (5-D). Las primeras tres dimensiones son espaciales: filas, columnas, y niveles (o latitud, longitud, y altura). La cuarta dimensión es el tiempo. La quinta dimensión es la enumeración de variables físicas múltiples, tales como temperatura, presión, humedad, etc. Además de los datos en el programa mismo, hay un número de parámetros necesarios para describir un banco de datos del VIS5D: los tamaños de las cinco dimensiones (número de filas, de columnas, de niveles, de pasos de tiempo, y de variables), de la posición y de la orientación geográfica de los datos (proyección del mapa), de los nombres de las variables, de los tiempos reales y de las fechas que se asociaron a cada paso de tiempo, etc. Convertir datos al formato de VIS5D Una de las dificultades que se puede tener para usar este programa es convertir los datos que se desea visualizar al formato VIS5D, en el caso de pronósticos, hay modelos numéricos que tienen programas incluidos para dicho fin, es decir convertir a dicho formato. A continuación veremos tres opciones para poder convertir datos atmosféricos al formato VIS5D. Los programas que son utilizados para este propósito están hechos en lenguaje Fortran y C, además existen algunas subrutinas ya realizadas o bien algunas librerías ya definidas, como ejemplo veremos solo lo utilizado en Fortran, de manera similar se trabaja con el lenguaje C. Se observa que al momento de utilizar los programas mencionados, se crean los archivos en el formato de VIS5D (esto es en una carpeta del mismo programa con el nombre de convert). a.- Usando Subrutinas: En este caso se utilizan los siguientes programas: - foo_to_v5d.f.- Asume una proyección rectangular del mapa de lat/lon y un sistema lineal igualmente espaciado de coordenada vertical. - foo2_to_v5d.f.- Permite cualquier proyección vertical del mapa, así como un diverso número de los niveles verticales para cada variable. En cualquier caso, cada programa de la conversión utiliza tres funciones para escribir el archivo en formato VIS5D: - V5dCreateSimple: Se utiliza para crear los archivos de v5d que especifican solamente los parámetros más básicos. J. Inga 22 - V5dCreate: Se utiliza para crear los archivos de v5d que permite utilizar parámetros más complicados. - V5dWrite: Se utiliza para escribir una sola rejilla tridimensional de datos a un archivo de v5d. Por rejilla es identificado un paso de tiempo y un número de variable física. - V5dClose: Cierra el archivo de v5d después de que se haya escrito la rejilla previa. Es el final del programa. El usuario puede leer, ejecutar y escribir los programas mencionados de acuerdo a sus intereses de investigación. Hay diversas versiones de estas funciones para los programas de Fortran así como también para C. A continuación, un ejemplo de los valores típicos que se pueden asignar a cada uno de las variables, si uno utiliza el programa de foo_to_v5d.f; de manera similar con algunas modificaciones, se trabaja con el foo2_to_v5d.f.: Asignación Contenido numtimes = 5 5 pasos de tiempo numvars = 4 4 variables físicas nr = 30 30 filas en cada rejilla 3-D ( por ejemplo latitud) nc = 40 40 columna en cada rejilla 3-D (por ejemplo longitud) nl = 20 20 niveles 3-D (por ejemplo nivel vertical) varname(1) = "U" U (este / oeste) componente de viento varname(2) = "V" V (norte / sur) componente de viento varname(3) = "T" Temperatura varname(4) = "P" Presión timestamp(1) = 140000 2:00:00 pm timestamp(2) = 141500 2:15:00 pm timestamp(3) = 143000 2:30:00 pm timestamp(4) = 144500 2:45:00 pm timestamp(5) = 150000 3:00:00 pm datestamp(1) = 94036 36th day of 1994 (Febrero 5) datestamp(2) = 94036 " datestamp(3) = 94036 " datestamp(4) = 94036 " datestamp(5) = 94036 " northlat = 60.0 límite norte de la caja está en 30 grados de latitud latinc = 1.0 hay 1 grado de la latitud entre cada uno de las 30 filas westlon = 100.0 límite oeste de la caja 3- D está en 100 grados de longitud loninc = 0.5 0.5 grado de la longitud entre cada uno de las 40 columnas bottomhgt = 0.0 el fondo de la caja está en los 0km (nivel del mar) hgtinc = 1.0 1 Km. entre cada uno de los 20 niveles (tope en 19.0km) b.- Uso de librerías Grib2v5d: El programa grib2v5d convierte el contenido de un archivo de GRIB (Un formato ideado por la Organización Mundial de Meteorología, usado para intercambiar información meteorológica en forma comprimida y binaria), en un archivo que se pueda ver con el VIS5D. El programa se escribe principalmente en FORTRAN 90 y se deriva de una versión anterior hecha por Paolo Patruno (mayor información visitar a la pagina Web de grib2v5d.sourceforge.net); hace uso de una librería en Fortran77 para el acceso de manejo a los archivos del grib, una biblioteca de C para la lectura y un subprograma de C (modificado proporcionado por Vis5d) para los archivos de la topografía y mapa. Generación de gráficas tridimensionales para el análisis y pronóstico del Tiempo y Clima– uso del VIS5A 23 c.- Uso de librerías de h5_utils: Este paquete de h5_utils es un sistema de las utilidades que se utiliza para la visualización y la conversión de datos científicos en el formato HDF5( Un formato estándar, libre, binario para los datos científicos multidimensionales). Además proporciona una herramienta simple para la visualización de imágenes, h_5utils también incluye programas para convertir el banco de datos HDF5 en los formatos requeridos por otro software de visualización libre (como el VIS5D). Ejecución del VIS5D Cuando se ejecuta el VIS5D, el panel de control y la ventana de exhibición aparecen (Figura 1). Dentro del panel de control el usuario puede cambiar las diversas opiniones para visualizar los datos, elegir qué datos desea observar, animarlo, agregando una información del mapa o de la topografía, así como la incorporación de nuevos datos, en si, hay diversas herramientas que se utilizan para tener una buena información de lo que se requiere analizar. Entre las principales herramientas que se encuentran en el panel de control son: a.- Cantidades escalares: Aquí están las variables físicas como T, P, V, U, etc. b.- Herramientas escalares del campo: Se utiliza para ver las variables físicas en contorno tanto verticales como horizontales y volúmenes (Isosurf, Contour Horiz, Sílice Vert, etc). c.- Herramientas vectoriales del campo: Se utiliza para ver la trayectoria vectorial de vientos tanto vertical como horizontal (Hwind1, Vwind1, Hstrean, etc.). d.- Modo de interacción: Se utiliza básicamente para mover los cortes transversales, acercar y alejar la imagen, entre otros. (Normal, Sílice, Label, etc.) Las otras herramientas de la parte superior son utilizadas para ver la topografía, los mapas, el reloj, vistas diversas (este, oeste, norte, sur), guardar imágenes, animaciones entre otros (topo, map, anime, etc) VIS5D permite que el usuario modifique la imagen para requisitos particulares, mientras que se pueda observar los datos. En la línea de comando hay diversas opciones que pueden ser manipuladas. El programa permite que el usuario cambie tamaños de la ventana, exhiba la fecha, la cantidad de memoria, etc. La mayoría de opciones pueden también ser cambiadas mientras que el programa está funcionando. Seleccionando el botón de la exhibición del panel de control, permite que el usuario cambie los valores prefijados. El usuario puede cambiar un nombre o unidad de variables, los tiempos, las fechas, la proyección, sistema de coordenadas vertical, o los niveles bajos. Además, el paquete de VIS5D incluye v5dimport, (Figura 2) que sirve para convertir y leer archivos de VIS5D. Con el J. Inga 24 v5dimport podemos seleccionar que tipo de variable física deseamos visualizar (temperatura, velocidad de viento, presión, etc.) y dándole otro nombre para que no se confunda con el archivo principal, podemos visualizar los datos escogidos. Figura 1. Panel de Control y Visualización del VIS5D Figura 2. V5dimport para convertir y leer archivos de VIS5D Generación de gráficas tridimensionales para el análisis y pronóstico del Tiempo y Clima– uso del VIS5A 25 APLICACIÓN A continuación mostraremos algunos ejemplos del manejo del programa VIS5D. Se ha utilizado la información de una simulación atmosférica, del sistema LAMPS (Sistema de Predicción a escala limitada de mesoescala), sectorizando el área del Perú, y océano adyacente, entre las longitudes 90° oeste y 52° oeste y latitudes 20° norte y 20° sur. Estos son datos que corresponden al mes de abril de 1989, el primer día de pronóstico (viernes) consta de 22 horas y del segundo día (sábado) se tiene 12 horas. Se puede ver en el grafico (Figura 3) las temperaturas del aire a 8 kilómetros del nivel del mar, con valores entre –24°C y – 52°C. Esta grafica se obtiene haciendo clic al campo de temperatura en el panel de control (Contour Horiz.). Allí le indica a que altura se desea mostrar la temperatura y la cantidad de líneas (isoterma) que se desea poner y dándole con el clic secundario se cambia el color de las líneas. Figura 3. Temperaturas a 8 Km. J. Inga 26 De la diversas herramientas que ofrece el VIS5D se encuentra la visualización de datos de los diversos campos que se esta trabajando, para puntos específicos en nuestra área de estudio. En la Figura 4, se muestra el campo de nubes en la zona de observación para un momento dado, así también, escogiendo en la pantalla un punto de información que deseamos ver (en el ejemplo un lugar en el departamento de Loreto de coordenadas: 7.824° latitud sur y 22.356° longitud oeste), nos muestra en la parte inferior izquierda toda su información. En este caso se esta trabajando con 10 campos los que son: • U: vientos zonales (este/oeste) • V: vientos meridionales (norte/sur) • W: vientos verticales • T: temperatura • P: presión • S: humedad específica • CWAT: concentración de nubes • RWAT: concentración de lluvias • SPD: velocidad del viento horizontal Figura 4. Concentración de Nubes y datos de campos Generación de gráficas tridimensionales para el análisis y pronóstico del Tiempo y Clima– uso del VIS5A 27 Usando VIS5D, también se puede hacer cortes transversales. Un corte transversal demuestra una sección de datos a una altura o un nivel particular. Los cortes transversales pueden ser horizontales o verticales. Usando esto, en el ejemplo (Figura 5), se muestra la concentración de lluvias con sus respectivas intensidades en un corte realizado a 6.9° latitud sur, observándose una mayor concentración de lluvias en la amazonía brasilera, también se puede observar la componente vertical de los vientos asociados a dicha concentración. Figura 5. Lluvias y Vientos CONCLUSIONES • VIS5D permite que los usuarios visualicen y analicen información tridimensional de manera sencilla y fácil. • Es muy útil en las aplicaciones en meteorología, por la cantidad de información que se maneja, y porque esta puede ser visualizada a diferentes niveles verticales de J. Inga 28 la atmósfera. Es útil para el análisis de información en tiempo real, así como también información proveniente de modelos numéricos. • El programa permite a los usuarios modificar la forma de visualización según requisitos particulares, tal que cada aspecto de la simulación puede ser explorado, también nos permite añadir la topografía del terreno y un mapa de la región, lo que le ayuda a localizar áreas del interés. • Entre las dificultades que se presentan al momento de trabajar con el VIS5D es el espacio de memoria requerido por la maquina al momento de procesar la información, así como el tamaño de los archivos resultantes, ello depende de cuantos campos se desea visualizar, y se hace más evidente cuando se realizan animaciones. Por otro lado, al ser un software libre, no siempre están disponibles todas las librerías para poder trabajar con este programa. AGRADECIMIENTOS El autor agradece a la Ing. Grace Trasmonte por el apoyo durante la elaboración de este informe. Así mismo, al Ing. Javier Viglanzoni por su colaboración en la logística. BIBLIOGRAFÍA http://www.ssec.wisc.edu/~billh/vis5d.htm http://www.ssec.wisc.edu/~billh/view5d.ht ml www.ssec.wisc.edu/mug/vis5d_guide/Curr ent/Vis5D-60.HTML http://grib2v5d.sourceforge.net Compendio de Trabajos de Investigación CNDG – Biblioteca Instituto Geofísico del Perú. V. 6 (2005) p. 29 - 40 ANÁLISIS DE LA VARIACIÓN ESPACIAL Y TEMPORAL DE LA TEMPERATURA Y PRECIPITACIÓN EN LA CUENCA DEL MANTARO CON LA BASE DE DATOS DE LAS ESTACIONES AUTOMÁTICAS Y SATÉLITE TANIA GIOVANNA CAYCHO BUSTAMANTE. Facultad de Ciencias Naturales y Matemática Escuela Profesional de Física Universidad Nacional del Callao tany_cb@yahoo.com Prácticas dirigidas por: Dr. Pablo Lagos E. Dra. Yamina Silva V. Centro de Predicción Numérica del Tiempo y Clima RESUMEN El presente trabajo comprende la evaluación y análisis de las variables meteorológica como la precipitación y temperatura del aire en la cuenca del Mantaro. Se evalúa la variación espacial utilizando datos de satélite y temporal con datos históricos obtenidos de la base de datos del IRI disponible en Internet y datos de la estación automática de Lircay. En general, los valores máximos de precipitación se observan durante los meses de diciembre a marzo, y mínimos en los meses de junio y julio. La temperatura más baja se observa en el mes de junio mientras que el valor máximo en el mes de noviembre. El procesamiento y análisis de los datos observados en las estaciones meteorológicas e imágenes del satélite GOES, se realizaron usando los lenguajes de programación FORTRAN y GRADS. INTRODUCCIÓN En los últimos años los daños causados por eventos hidrometeorológicos se han incrementado provocando inundaciones, principalmente debido a la ocurrencia de eventos extremos como precipitaciones muy intensas. Estos eventos han producido mayores daños debido al crecimiento poblacional y la urbanización en sitios de potencial peligro. El problema se puede enfrentar mediante medidas de prevención y dar una respuesta más realista a la previsión del riesgo hidrológico utilizando modelos numéricos para el pronostico de lluvias. El estudio del clima involucra la construcción o reconstrucción de series temporales de datos climáticos. La variación de los fenómenos climáticos a lo largo del tiempo provee una medida ya sea cualitativa o cuantitativa de la variabilidad del clima. La variación espacial y temporal de los datos climáticos que incluyen la temperatura, la precipitación (lluvias) para la cuenca del Mantaro varía con las estaciones y con los años Durante la historia más reciente, los científicos han sido capaces de construir series temporales del clima a partir de datos instrumentales observados por estaciones meteorológicas o usando imágenes de satélites meteorológicos para la evaluación del clima. T. Caycho 30 OBJETIVOS Objetivos Generales Análisis de la variación espacial y temporal de la temperatura y precipitación en la cuenca del Mantaro con la base de datos de las estaciones automáticas y de satélite. Objetivos Específicos * Elaborar programas para el procesamiento de datos meteorológicos. * Analizar la variación espacial de la precipitación utilizando datos obtenidos del IRI1, así como su distribución mensual climatológica utilizando histogramas. * Analizar los datos diarios de temperatura y precipitación de la estación automática de Lircay de los años de 2001, 2002 y 2003. * Observar y analizar la precipitación estimada del satélite GOES12 para la cuenca del Mantaro. MARCO TEÓRICO Instrumentos de medición meteorológica La meteorología utiliza instrumentos esenciales, como el barómetro, el pluviómetro, el termómetro y el higrómetro para las mediciones de presión atmosférica, precipitaciones, temperatura y el grado de humedad del aire [1]. Las variables meteorológicas son medidas por la estación meteorológica convencional, la estación meteorológica 1 International Research Institute for Climate Predicction automática y el satélites meteorológicos [2]. Descripción teórica de las variables climáticas. Tiempo y clima: Las actas de la Conferencia Mundial del Clima (1979) adoptaron las siguientes definiciones para el tiempo como el estado de la atmósfera en un lugar y momento determinado; y el clima como la sucesión periódica de tipos de tiempo. Elementos del tiempo y clima: La atmósfera en sus niveles inferiores, se caracteriza por un gran dinamismo que influye en el surgimiento de los llamados fenómenos meteorológicos [3]. Las circunstancias atmosféricas de un lugar y en un momento determinado, son el resultado de la combinación de varios elementos comunes al tiempo y al clima, como son la temperatura que refleja el estado energético del aire, el cual se traduce en un determinado grado de calentamiento; la presión atmosférica del aire determina el peso de las masas de aire que ejerce la atmósfera sobre la superficie de la tierra, expresada en milibares; el aire que se desplaza paralelamente a la superficie terrestre se llama viento y los movimientos verticales del aire se llaman corrientes; la humedad atmosférica es el estado que presenta la atmósfera en relación con el vapor de agua que contiene. Mientras que las precipitaciones Análisis de la variación espacial y temporal de la Temperatura y Precipitación en la Cuenca del Mantaro 31 es el producto de la condensación atmosférica, que puede ser sólida o líquida. La cantidad de lluvia de un día se mide en el pluviómetro y se expresa en mm. Procedimiento para la adquisición de datos climatológicos. Selección y adquisición de datos: Existen diferentes instituciones que se encargan de recopilar y almacenar datos climatológicos, estas instituciones pueden ser nacionales como internacionales. En el presente trabajo se seleccionaron datos de tres fuentes diferentes: a) Las series históricas de precipitación mensual fueron obtenidas de la base de datos del IRI. Los datos de precipitación mensual son analizados para algunas estaciones ubicadas dentro de la cuenca del Mantaro, cuyas coordenadas geográficas están limitadas por 13.5°S - 10°S y 76.8°W - 74.1°W. La lista de estaciones obtenidas de la base de datos del IRI se muestra en el Tabla1. b) Datos diarios de temperatura del aire y precipitación registrados por la estación automática ó llamada DCPs (Data Collection Platform) de Lircay ubicado a los 12.98 de latitud Sur y 74.72 de longitud Oeste, perteneciente al proyecto multinstitucional titulado "Mejoramiento de la capacidad de pronostico y evaluación del fenómeno El Niño para la prevención y mitigación de desastres en el Perú" con el apoyo y cooperación del Banco Mundial (BM), para los años de 2001, 2002 y 2003. c) Datos diarios de precipitación estimada en base a la imagen infrarroja del satélite GOES 12 operados por la NOAA de diciembre 2003, enero y febrero 2004. Estos datos son experimentales y aun no han sido validados para el Perú, sin embargo es una información importante porque muestra la variación espacial de las lluvias estimadas [6]. La técnica de estimación de lluvia por satélite se denomina Auto estimado, donde, la taza de precipitación es obtenida a partir de la temperatura de brillo del tope de la nube derivada de la banda infrarroja (10.7 µm) del satélite GOES (Geostationary Operational Environmental Satellite) de resolución espacial 4 km x 4 km [5]. Las estimaciones están en base a intervalos de temperatura (ºK) de la atmósfera, y la lluvia precipitada (mm) obteniéndose una curva exponencial como se muestra en la Figura 1. T. Caycho 32 Tabla 1. Lista de estaciones meteorológicas ubicadas en la cuenca del Mantaro y alrededores que se encuentran en la base de datos del IRI. Nº Id de la estación Nombre de la estación Longitud (º) Latitud (º) Elevació n (m) Fecha (inicio) Fecha (final) 1 84570000 Cerro_de_Pasco 76.25W 10.68S 4334 Jul 1949 Feb1980 2 84570002 Fundición 76.27W 10.75S 4268 Jan1957 May1972 3 84570007 Laguna_Huarón 76.42W 11.02S 4649 Jan1957 Dec1974 4 84570015 Quiulacocha 76.28W 10.7S 4299 Nov1952 Dec1974 5 84570016 S.J.Pallanga 76.45W 11.15S 4649 Jan1957 Dec1974 6 84570018 Shelby 76.23W 10.82S 4131 Jan1957 Dec1974 7 84570019 Upamayo 76.28W 10.92S 4079 Oct1963 Dec1974 8 84600001 Atocsaico 76.07W 11.3S 4199 Jan1957 Dec1974 9 84600003 Corpacancha 76.22W 11.37S 4249 Jan1957 Dec1974 10 84600005 Huallacocha 76.1W 11.77S 4399 Jan1957 Dec1974 11 84600009 Huascacocha 76.08W 11.58S 4499 Jan1957 Dec1974 12 84600012 Malpaso 76.03W 11.4S 3749 Jan1943 Dec1974 13 84600013 Marcapomacocha 76.33W 11.4S 4412 Oct1964 May1981 14 84600015 Morococha 76.1W 11.63S 4539 Jan1943 Dec1974 15 84600016 Pachachaca 76.0W 11.62S 3999 Jan1949 Nov1977 16 84600019 Pinascocha_HDA 75.83W 11.82S 4299 Apr1957 Apr1974 17 84600021 Pomacocha 76.13W 11.73S 4265 Mar1938 Dec1977 18 84600023 Punabamba 76.08W 11.48S 4099 Apr1957 Apr1974 19 84600025 San_Cristobal 76.05W 11.73S 4699 Jan1957 Dec1974 20 84630000 Huayao/Huancay 75.3W 12.0S 3350 Jan1922 Dec1980 21 84630003 Comas 75.08W 11.72S 3299 Dec1963 Sep1978 22 84630004 Consav_HDA 75.63W 11.98S 3882 Apr1957 Apr1974 23 84630005 Huancalpi 75.25W 12.57S 3879 Jan1965 Aug1981 24 84630006 Huantan 75.82W 12.45S 3271 Oct1963 Dec1976 25 84630007 Jauja 75.5W 11.8S 3410 Jan1935 Dec1979 26 84630008 Laive_HDA 75.4W 12.3S 4095 Sep1963 Jun1981 27 84630009 Matibamba 74.82W 12.08S 2199 Aug1963 Jul1977 28 84630010 Pachacayo_HDA 75.72W 11.77S 3599 Apr1957 Dec1974 29 84630011 Pampas_Colonia 75.88W 12.63S 3378 Sep1963 Dec1976 30 84630012 Pesqueria_ING 75.27W 11.88S 3399 Aug1963 Dec1981 31 84630014 Ricran 75.52W 11.53S 3729 Apr1965 Dec1981 32 84630015 Runatullo 75.02W 11.62S 3149 Sep1966 Apr1981 33 84630017 San_Lorenzo 74.78W 12.3S 2599 Sep1963 Dec1979 34 84630018 San_Pedro_Chucl 75.5W 11.75S 3379 Sep1963 May1978 35 84630020 Tarma 75.7W 11.42S 3050 Nov1963 Dec1981 36 84630021 Vilca 75.83W 12.12S 3815 Sep1963 Jan1977 37 84630022 Viquez 75.23W 12.17S 3184 Jan1963 Jun1981 38 84630023 Yauricocha 75.72W 12.32S 4521 Jan1943 Jun1972 39 84630024 Yauyos 75.92W 12.37S 2870 Oct1963 Dec1977 40 84673000 Ayacucho 74.2W 13.2S 2761 Apr1962 Nov2000 41 84673005 Lircay 74.72W 12.98S 3270 Dec1964 Dec1979 42 84673006 Luricocha 74.27W 12.9S 2579 Nov1963 Dec1977 43 84673008 Paucarbamba 74.53W 12.55S 3360 Dec1964 Mar1982 44 84673010 San_Miguel 73.98W 13.02S 2660 Sep1964 Sep1977 45 84673009 La_Quinua 74.13W 13.05S 3099 Jan1964 Mar1981 46 84673011 San_Pedro_Cachi 74.4W 13.08S 3187 Apr1965 Sep1981 47 84680003 Choclococha 75.07W 13.15S 4549 Jul1958 Dec1980 48 84680006 Huancavelica 74.98W 12.78S 3669 Dec1963 Sep1974 49 84680012 San_Genaro 75.1W 13.2S 4569 May1958 May1975 50 84680017 Telepaccha 75.3W 12.75S 4399 Sep1963 Oct1981 Fuente:IRI. [4] Análisis de la variación espacial y temporal de la Temperatura y Precipitación en la Cuenca del Mantaro 33 Figura 1. Relación de la estimación de la lluvia por el radar en función a la temperatura [5] Los puntos representan el promedio de la estimación de la lluvia por radar para cada intervalo de temperatura, y la curva representa la regresión dada por: [ ])2.1*(**)2*(*10*6382.3exp*)11(*10*1183.1 TR −−= (1) Donde R es la proporción de la lluvia y T la temperatura en grados Kelvin. * Cuando el punto indica el máximo valor, representa la cima de la nube relativamente alta y el punto se localiza en ambientes más fríos, por lo tanto la proporción de la lluvia dada por la curva de regresión permanece inalterable. * Si el punto muestra el mínimo valor, indica la presencia de la cima de la nube relativamente baja y en ambientes cálidos, entonces la proporción de lluvia tiende a cero. * Si el punto se localiza en promedio entre el máximo y el mínimo, y a temperatura ambiente, la proporción de la lluvia tiende a cero [5]. d) Control de calidad de datos: El control de calidad de los datos tiene como objetivo eliminar o detectar los datos dudosos o defectuosos. Para este estudio el proceso de detectar y luego eliminar los errores o defectos se hace mediante el método subjetivo visual, verificando los datos provenientes del IRI y DCPs de formato ASCII y de satélite en formato binario. Para ello se verifica el formato de la fuente, el sitio y ubicación del dato con respecto a las fechas, la coherencia en la magnitud y las unidades en las que se expresan las variables medidas, se verifica la transcripción correcta a la base de datos digital y posteriormente se omite en caso que tenga un valor dudoso. Análisis estadístico de los datos. Una vez que se han recolectado información se realiza el análisis de los datos climáticos. El objetivo del análisis estadístico es identificar el comportamiento sistemático en un conjunto de datos, que permita obtener una señal en los datos que pueda distinguirse del ruido. Para ello el análisis descriptivo está encaminado únicamente a documentar aspectos particulares de las variaciones presentadas en la serie de datos (señal). Los índices calculados incluirán la media y la varianza (o la desviación estándar). También se notará la ocurrencia de eventos extremos, ciclos y tendencias. Mientras que el análisis investigativo se dirige a verificar hipótesis predefinidas. T. Caycho 34 Las hipótesis deben tener a priori una sólida base científica Cualquier conjunto de datos usado para análisis estadístico debe ser representativo de los procesos físicos relevantes; lo suficiente en cantidad para soportar el método estadístico usado y ser preciso y confiable (homogéneo). Visualización y control de los datos. Para el sistema de selección de datos se elaboran programas en el lenguaje de programación FORTRAN (FORmula TRANslator) 90 para la lectura y el procesamiento de datos de varios años, con el fin del cambio de formato ASCII a formato Binario. El control y visualización de los datos se desarrollo usando el programa GRADS (GRid Analysis and Display System) el cuál lee los datos en formato binario y la ubicación geográfica de las estaciones para los tres tipos fuentes de datos. Según la variable meteorológica se realizó el análisis estadístico y la visualización de los datos. RESULTADOS Análisis de los datos de precipitación obtenidos de la base de datos del IRI En la cuenca del Mantaro se encuentran instalados varias estaciones meteorológicas tanto en la zona norte, centro y sur. El promedio multianual de la precipitación mensual, indica que estas se inician en el mes de agosto. Las máximas precipitaciones ocurren en los meses de enero a marzo y las mínimas en junio- julio. La estación de Marcapomacocha comparada con las otras estaciones alcanza un valor máximo de 240 mm en el mes de febrero y en junio se presentan las mínimas precipitaciones llegando hasta los 10 mm como se muestra en la Figura 2. La Figura 3, muestra la variabilidad mensual de precipitación. Los máximos valores se presentan a comienzos y final de cada año (estación de verano) y un mínimo cerca a cero coincidente con la estación de invierno, es decir a mediados de cada año. La máxima precipitación con valor de 480mm se observa en la estación de Marcapomacocha en el mes de enero de 1978. Compendio de Trabajos de Investigación CNDG – Biblioteca Instituto Geofísico del Perú. V. 6 (2005) p. 29 - 40 Figura 2. Climatología mensual de precipitación (mm) para 6 estaciones de la cuenca del Mantaro T. Caycho 36 Figura 3. Precipitación mensual en 6 estaciones representativas de la cuenca del Mantaro Figura 3. Continua.../. La anomalía mensual de la precipitación en la estación de Atocsaico con respecto a las otras estaciones muestra valores por debajo de su normal (-100 mm) en el verano austral del año 1970 y valores sobre su mensual (+255 mm) en el verano del año 1974 como muestra la Figura 4. Análisis de la variación espacial y temporal de la Temperatura y Precipitación en la Cuenca del Mantaro 37 Figura 4. Anomalía mensual de la precipitación en 6 estaciones de la Cuenca del Mantaro Análisis de la base de datos de precipitación y temperatura de la estación automática de Lircay. La Figura 5, muestra la variabilidad en la suma anual de la precipitación durante los 2001 y 2003. En el año 2002 se registró el mayor acumulado en comparación a los otros dos años, con un valor máximo de 765 mm. Figura 5. Precipitación acumulada anual para los años 2001, 2002 y 2003 Análisis de la precipitación mensual: La estación de Lircay en enero del año 2001 del mes registró una precipitación de 180 mm, mas de 100mm superior a los registros de los años 2002 y 2003. Los mínimos valores, cercanos a cero mm, s registraron en el mes de junio para los años 2002 y 2003. Figura 6. Precipitación acumulada mensual en la estación automática de Lircay para el año 2001,2002 y 2003 Análisis de temperatura mensual del aire. La estación de Lircay en el mes de junio del 2001, 2002 y 2003, muestra valores mínimos cercanos a los 2 ºC y registra un valor máximo de 24 ºC a fines del mes de setiembre del año 2003 debido a la variación de la temperatura y el cambio estacional de la tierra (Figura 7). T. Caycho 38 Figura 7. Temperaturas máxima y mínima mensuales en la estación automática de Lircay para el año 2001, 2002 y 2003. Análisis de la precipitación estimada del satélite georeferenciado Goes12 Los datos de precipitación estimada por satélite muestran para el mes de diciembre del año 2003, máximos valores (140mm) en la parte centro-occidental de la cuenca del Mantaro y un valor mínimo de 26 mm al sur este de la cuenca como se muestra en la Figura 8. Figura 8. Precipitación estimada acumulado para el mes de diciembre del 2003 En la Figura 9, se observa un valor máximo de 220 mm de precipitación estimada para enero del 2004 en la región norte y centro occidental y valores mínimos de 25 mm en la parte centro- oriental y sur. Figura 9. Precipitación estimada, acumulado mensual para el mes de enero del 2004 El análisis de la imagen satelital (Figura 10), para la cuenca del Mantaro registra en el mes de febrero del 2004, mínimos valores (20 mm) en la parte norte y sur- este mientras que los valores máximos (160mm) en la región centro y sur- occidental. Figura 10. Precipitación estimada acumulado para el mes febrero del 2004 CONCLUSIONES El procesamiento de datos provenientes de diversas fuentes tanto de las estaciones automáticas o de satélites tiene el mismo procedimiento ya que para su evaluación y análisis se usan los lenguajes de programación FORTRAN y GRADS. Análisis de la variación espacial y temporal de la Temperatura y Precipitación en la Cuenca del Mantaro 39 La climatología calculada como el promedio multianual de los datos de precipitación para las estaciones obtenidas de la base de datos del IRI, muestra valores máximos en el mes de febrero y mínimos en el mes de junio. Los datos mensuales presentan máximos valores a comienzos y finales de cada año, observándose el máximo valor de 480 mm en el mes de enero de 1978 y mínimos cercano a cero en el mes de junio de los años 1967, 1978 y 1979 para la estación de Marcopomacocha. La máxima anomalía negativa presenta valores de -100 mm en el verano austral del año 1970 y la máxima anomalía positiva (+255) en el verano del año 1974 en la estación de Atocsaico. Los datos de precipitación registrados en la estación automática de Lircay en los últimos 3 años, indican que en suma anual, el año 2002 fue mas lluvioso que el 2001 y 2003, sin embargo en enero del año 2001 se registró el valor máximo (180 mm) de precipitación. Los valores mínimos de la temperatura media en la estación de Lircay, se registraron durante el mes de junio con valores cercanos a 2ºC y los máximos a inicios del mes de noviembre, siendo el año 2003 que presentó el valor máximo (24 ºC) con respecto a los otros años. La precipitación mensual estimada con datos de satélite registró para el mes enero del 2004 máximas de 220 mm en la parte norte y centro occidental y mínimas en febrero del 2004 de 20 mm en la región norte y sur este con respecto a la cuenca del Mantaro. Las variables climatológicas de temperatura y precipitación varían de acuerdo a las estaciones del año y su ubicación geográfica, con respecto a las diversas altitudes tanto para los datos registrados por estaciones o monitoreados en tiempo real por el satélite. AGRADECIMIENTOS Mi agradecimiento al INSTITUTO GEOFÍSICO DEL PERÚ (IGP), por haberme dado la oportunidad de realizar mis practicas Pre_Profesionales en el Centro de Predicción Numérica del Tiempo y Clima. Al Dr. Pablo Lagos, Dra. Yamina y al Fis. Berlín Segura, en el asesoramiento del presente informe a su vez al equipo CPNTC por su constante apoyo. A mi familia y amigos por la confianza en mi persona. BIBLIOGRAFIA [1]http://www.imn.ac.cr/educa/instrument os/tmaxmin.htm [2]http://www.geocities.com/joanballester/ Espanyol/Meteo/Satelits_esp.htm T. Caycho 40 [3]http:// www.fca.unl.edu.ar/Clima/01- TiempoyClima.pdf [4]http://ingrid.ldeo.columbia.edu [5]http://www.rsd.gsfc.nasa.gov/goeseast/ [6]http://orbit35i.nesdis.noaa.gov/arad/ht/f f/gilberto.html Análisis de la variación espacial y temporal de la Temperatura y Precipitación en la Cuenca del Mantaro 41 Compendio de Trabajos de Investigación CNDG – Biblioteca Instituto Geofísico del Perú. V. 6 (2005) p. 41 - 48 DISEÑO Y CONSTRUCCIÓN DE LA ANTENA PERSEUS VHF CÉSAR DE LA JARA SÁNCHEZ Especialidad de Ingeniería Electrónica Facultad de Ciencias e Ingeniería Pontificia Universidad Católica del Perú cdelajara@jro.igp.gob.pe Prácticas dirigidas por. Dr. Ronald Woodman Dr. Martín Sarango Radio Observatorio de Jicamarca - IGP RESUMEN Durante los meses de Agosto y Setiembre del año 2004 NASA realizó una campaña de cohetería en la que se lanzaron 14 cohetes suborbitales en el Pacífico del Sur y en la que se realizaron distintos experimentos que tenían por objetivo estudiar las alteraciones en la ionosfera causadas por la interacción entre el sol y el campo magnético de la tierra. En el presente trabajo se describen el diseño y fabricación de la antena VHF utilizada en el experimento peruano (experimento PERSEUS) que formó parte de la campaña de cohetería mencionada. Se realiza un análisis del tipo de antena empleado así como sus fundamentos físicos y las ventajas e inconvenientes que se presentan al utilizarlas. También se dan algunos detalles de la construcción y ensamblaje de la antena. INTRODUCCIÓN El objetivo era diseñar una antena que pudiera ser montada en la superficie de un cohete y cuyo propósito sería la emisión de una onda de radio de alta frecuencia durante el vuelo del mismo. El cohete llegaría hasta la ionósfera, realizando un vuelo de varios minutos de duración, soportando esfuerzos mecánicos muy grandes, por lo que era necesario que la antena fuera diseñada para resistir intensas aceleraciones, altas temperaturas y fuertes vibraciones y sacudidas. Dos antenas irían montadas en cada cohete, y emitiría una señal de 37.866 Mhz. A esta frecuencia la longitud de onda (λ) de la señal transmitida es de 7.2 metros, lo que significa que la construcción de una antena de un cuarto de longitud de onda requeriría de 1.98 metros, algo irrealizable debido a las limitaciones de espacio. Se disponía de aproximadamente 1 metro sobre el cohete para el montaje de cada antena, lo que obligaba a que estas tuvieran una longitud menor a ¼ de λ convirtiéndose así en antenas eléctricamente cortas. ANTENAS CORTAS Cuanto más baja sea la frecuencia que se utilice en un sistema de transmisión, se encontrará una mayor dificultad en construir una antena de ¼ de longitud de onda. No es absolutamente necesario construir una antena entera, es decir, se puede utilizar una antena cuya longitud sea una fracción de ¼ λ, pero al hacer esto se deben buscar ciertos compromisos que C. De la Jara 42 permitan mantener la eficiencia, impedancia y ancho de banda en valores aceptables(ARRL, 2003). Una antena corta presenta una resistencia de radiación baja. A medida que la relación Largo de la Antena / Longitud de onda se reduce, el transmisor tiene que entregar corrientes cada vez más altas para radiar una cantidad significativa de potencia. Una antena corta presenta una impedancia con componente reactiva grande. La expresión Z = R + jX indica que la impedancia puede ser modelada como una resistencia R, en serie con una reactancia X, con la misma corriente fluyendo a través de ambas. En un sistema como este, la corriente no está en fase con el voltaje. Los máximos de corriente y voltaje no ocurren al mismo tiempo por lo que no se hace un uso eficiente de la potencia, es decir, el producto de la corriente por el voltaje no alcanza el valor que podría si estas estuviesen en fase. Figura 1. Corriente y voltaje fuera de fase La energía transferida a la antena por la componente activa tiene dos efectos: producir calentamiento en la estructura de la antena, es decir que exista la presencia e pérdidas por efecto Joule; y generar la emisión de ondas RF (Booker, 59). Otra desventaja de estas antenas es que cuanto más corta sea, más sensible se hace a los ajustes, además el valor final de la impedancia dependerá también de la posición de objetos que la rodean. CONSTRUCCION DE LA ANTENA PERSEUS VHF En el diseño de la antena existían dos parámetros sobre los que no existía mucho control: 1) La longitud de la antena. Esta debía ajustarse a las dimensiones del cohete 2) La separación entre la antena y la superficie del cohete Debido a las limitaciones de longitud, fue necesario diseñar una antena cuyo largo sea menor a ¼ de longitud de onda. A medida que se reduce la longitud, la magnitud de la reactancia capacitiva se incrementa y la antena se comporta como un condensador cuya capacitancia es proporcional al largo de la varilla. Esta es una de las razones por las que la antena no es un radiador eficiente; no solo tiene una baja resistencia de radiación, lo que incrementa el flujo de corriente, sino que también tiende a presentar una impedancia distinta a la característica resistiva del transmisor. Además, las altas corrientes Diseño y construcción de la antena Perseus VHF 43 necesarias para radiar potencia, disipan calor por efecto Joule(ARRL, 2003). Aún si la reactancia capacitiva es compensada adecuadamente, la antena es ineficiente, posee un ancho de banda estrecho y tiene propensión a perder sintonía por la más pequeña alteración en la distribución de objetos en su vecindad. No es posible tener a la vez una antena corta, eficiente y de banda ancha, por lo que el diseño exigía la búsqueda de compromisos. Siendo la limitación de espacio el factor condicionante, se buscó lograr un equilibrio adecuado entre eficiencia y ancho de banda. La escasa separación entre la varilla de la antena y la superficie del cohete (menor a 10 cm), ocasiona que la antena tenga un comportamiento semejante al de una línea de transmisión en que las corrientes opuestas cancelan la energía radiada. A pesar de que esta separación es una pequeña fracción de la longitud de onda, es lo suficientemente grande como para permitirle a la antena radiar energía. CARGA CAPACITIVA Para compensar la longitud faltante de la antena corta, se le cargó con un condensador en extremo opuesto al punto de alimentación, este condensador permitió llevar el punto de resonancia de la antena al valor requerido. Se diseñó un capacitor con dos placas de acero, material cerámico y arandelas de acero (Figura 2), se utilizaron tornillos para sujetarlo al cohete. El material utilizado como dieléctrico debía ser capaz de resistir las altas temperaturas que produce la fricción del aire durante el vuelo del cohete. Se tuvieron dos alternativas, teflón, que no es muy costoso y es fácil de conseguir; y cerámica, que es costosa y había que comprarla fuera del país. Finalmente se utilizó la cerámica porque esta puede soportar temperaturas más altas que el teflón. Figura 2. Capacitor de placas de acero Para inmovilizar las arandelas de cerámica se perforaron cavidades en la placa superior de acero (Figura 3). C. De la Jara 44 Figura 3. placa superior de acero. Vista inferior (arriba) y vista superior(abajo) Haciendo un análisis detallado observamos que el condensador está formado por la suma de varias capacitancias en paralelo (Figura 4) Figura 4. Secciones del capacitor C1 es una capacitancia cilíndrica formada por el tornillo y la placa de acero superior, tiene dieléctrico de cerámica. C2 es una capacitancia de placas en paralelo, formada por las arandela de acero y la placa de acero superior, tiene dieléctrico de cerámica. C3 es una capacitancia cilíndrica formada por el tornillo y la placa de acero superior, tiene dieléctrico de aire C4 es una capacitancia de placas en paralelo, formada por las 2 placas de acero, tiene dieléctrico de aire. C5 es una capacitancia cilíndrica, formada por el tornillo y la placa de acero inferior, tiene dieléctrico de cerámica. C6 es una capacitancia de placas en paralelo, formada por las placas de acero superior e inferior, tiene dieléctrico de cerámica. La capacitancia total viene dada por: Ctot = 4(C1+C2+C3+C5+C6) + C4 Debido a la curvatura del cohete, la placa inferior de acero debía fabricarse con una curvatura en la cara inferior, cuyo radio era igual al radio del cohete (Figura 5) La otra cara (la superior) era plana y es donde se apoyan las arandelas de cerámica. Figura 5. Placa inferior de acero El capacitor tiene además un trimmer ubicado en la placa de acero superior (Figura 6). Con este se puede variar la capacitancia para acomodar el punto de Diseño y construcción de la antena Perseus VHF 45 resonancia exactamente en el valor deseado. Para inmovilizar el trimmer a la placa se utilizó un prisionero. que antes del vuelo del cohete se fijaba con sellador de tornillos. Figura 6. Trimmer DISEÑO MECÁNICO El material seleccionado para la fabricación de la antena fue acero inoxidable 316, este tipo de acero es usado en una gran variedad de aplicaciones industriales, además de ser altamente resistente a la corrosión y capaz de soportar temperaturas de hasta 1200oC. La antena completa está compuesta por (Figura 7): - Una varilla de acero de 80 cm de longitud y 8 mm de diámetro. - Un capacitor de carga (de placas paralelas, de acero) en un extremo, que también sirve para sujetar la antena al cohete. - Una base de antena Spike1 que servía como punto de alimentación y como elemento de sujeción (Figura 7). Figura 7. Antena completa Como puede verse en la Figura 7, la varilla de acero no es recta, se realizó un diseño curvo para mantener a la antena lo más alejada posible de la superficie del cohete. La varilla se soldó al capacitor de carga para formar una sola pieza. Para sujetar la varilla a la base spike se fabricó un vástago de doble hilo (Figura 10). El doble hilo es necesario para que ambas piezas queden alineadas luego de ser ajustadas (Figura 9) 1 Las Spike son un tipo de antena que ha sido utilizada durante décadas para la transmisión de ondas de radio en aplicaciones aeroespaciales. Su longitud, en correspondencia con las frecuencias que emplean, son de unos pocos centímetros: 5.6 cm @ 550 Mhz, hasta 25.4 cm @ 216 Mhz, lo que las hace inútiles para nuestra aplicación. La antena completa no podía ser utilizada para transmitir 37.866 Mhz, pero su base resultaba adecuada para fijar nuestra varilla al cohete. Estos dispositivos nos fueron proporcionados por NASA. C. De la Jara 46 Figura 8. Base spike Figura 9. Vástago de doble hilo Figura 10. Varilla de acero unida a la base spike por medio de un vástago de acero de doble hilo RENDIMIENTO Y RESULTADOS Para poder realizar las pruebas de rendimiento a la antena y con el objetivo de obtener un buen grado de aproximación a las condiciones reales, se construyó una estructura de acero y aluminio de dimensiones similares a la de los cohetes que serían empleados en los lanzamientos (Figura 11). Figura 11. Estructura de Aluminio utilizada para probar las antenas Las pruebas se hicieron montando dos antenas sobre la estructura, opuestas una a la otra y alejándolas de cualquier otra estructura metálica que pudiera producir acoplamiento. Diseño y construcción de la antena Perseus VHF 47 Figura 12. Antenas montadas sobre estructura metálica Las antenas mostraron un buen comportamiento desde los primeros ensayos. La impedancia se acercaba a los 50 ohmios a la frecuencia de resonancia, y el ancho de banda era de 110Khz, que se encontraba dentro de los límites permisibles para nuestra aplicación (Figura 13). Figura 13. Curva de frecuencia vs Pérdidas de retorno de las antenas. Se observa un ancho de banda de 110 Khz y una impedancia de casi 50 ohm (RL = -20dBm) También se realizaron pruebas con la antena en el laboratorio Goddard de NASA (Figura 14), las que nos ayudaron a implementar mejoras al diseño. Figura 14. Antena durante las pruebas de vibración en NASA CONCLUSIONES Hemos visto que una antena corta no posee un alto grado de eficiencia, no tienen un buen ancho de banda y requiere el uso de corrientes altas para radiar energía, pero también sabemos que si hacemos uso de compromisos podemos mantener estos parámetros en niveles aceptables para una determinada aplicación. Las antenas desarrolladas en el Radio Observatorio de Jicamarca, junto con un sistema de transmisión, fueron montadas en 6 cohetes de NASA (Figura 15) durante la campaña de cohetería EQUIS II del 2004, mostrando un buen rendimiento C. De la Jara 48 durante el vuelo y logrando obtener valiosa información que luego serviría para completar el primer proyecto aeroespacial peruano (Figura 16). Figura 15. Antenas montadas en uno de los cohetes, listo para ser lanzado Figura 16. Lanzamiento de uno de los cohetes AGRADECIMIENTOS Mi gratitud a los doctores Martín Sarango y Jorge Chau, por haber depositado su confianza en mí y al doctor Ronald Woodman por compartir su experiencia y conocimientos. BIBLIOGRAFÍA ARRL (2003) The national association for Amateur Radio, 20th edición, 2003. The ARRL Antenna Book. Booker, H. (1959) An Approach to Electrical Science. Cornell University. Compendio de Trabajos de Investigación CNDG – Biblioteca Instituto Geofísico del Perú. V. 5 (2004) p. 49 - 60 DESARROLLO DE CONTENEDORES PARA LOS TRANSMISORES DEL PROYECTO ESPACIAL PERSEUS I JORGE ANTONIO CHOCOS NÚÑEZ & KARIM MILAGROS KUYENG RUIZ Sección Electricidad y Electrónica Facultad de Ciencias e Ingeniería Pontificia Universidad Católica del Perú jchocos@jro.igp.gob.pe, kkuyeng@jro.igp.gob.pe Prácticas dirigidas por: Dr. Martín F. Sarango Radio Observatorio de Jicamarca -IGP RESUMEN Entre Agosto y Setiembre de 2004 se llevó a cabo la campaña de cohetería científica EQUIS-2 en la Base de Lanzamientos. Experimentales Ronald Reagan, en el Pacífico Ecuatorial. El Radio Observatorio de Jicamarca del Instituto Geofísico del Perú, participó en esta campaña con el primer experimento espacial peruano, destinado a medir el Contenido Total de Electrones (TEC: Total Electron Content) en la Ionosfera. El proyecto denominado PERSEUS (Peruvian Rocket Sounding Experiment form Upper-Atmosphere Studes – Experimento peruano de cohetería científica para estudios de la alta atmósfera) comprendió la concepción, diseño y desarrollo de un experimento basado en la técnica de fase diferencial, utilizando radio-transmisores transportados por cohetes. El presente trabajo contiene la descripción del sistema de contenedores (o bandejas) aplicables para la electrónica de los transmisores del proyecto PERSEUS I. En total se desarrollaron ocho transmisores de frecuencia dual: 37.866 MHz (VHF) y 567.99 MHz (UHF); cada uno de los cuales constaba de tres bandejas contenedoras para los osciladores, los amplificadores y los acondicionadores de impedancias respectivamente. A continuación se describe brevemente el experimento, para luego exponer los detalles de diseño, fabricación y montaje de las mencionadas bandejas. INTRODUCCIÓN Durante agosto y setiembre del año 2004 se realizó en el atolón Kwajalein, República de las Islas Marshall, la campaña de cohetería científica EQUIS-2. Esta fue organizada por la NASA, junto a otras universidades e instituciones de los Estados Unidos, con la participación del Instituto Geofísico del Perú. El objetivo de esta campaña fue realizar estudios de la alta atmósfera por medio de instrumentos científicos a bordo de cohetes suborbitales. La participación del IGP fue a través del proyecto PERSEUS: PEruvian Rocket Sounding Experiment for Upper- Atmosphere Studies, que constituye el primer experimento espacial peruano. Este proyecto consistió en la concepción, diseño y desarrollo de un instrumento para medir el Contenido Total de Electrones (TEC) en la Ionosfera por medio de un experimento de Fase Diferencial. El experimento fue programado para ser lanzado a bordo de ocho cohetes de NASA. Por tal motivo se construyeron ocho radiotransmisores con dos frecuencias de salida: una en VHF -37.866 MHz- y la otra en UHF –567.99 MHz- La J. Chocos & K. Kuyeng frecuencia mayor se obtiene al multiplicar la frecuencia menor por un factor de 15, así ambas frecuencias tienen una diferencia de fase constante al momento de ser transmitidas. Cuando el cohete que lleva el transmisor atraviesa la ionosfera, la fase de la frecuencia VHF se altera en forma proporcional al TEC, pero la frecuencia UHF no sufre esta alteración. Es así como, utilizando receptores digitales en tierra para determinar la diferencia de fase entre las dos frecuencias, se obtiene el TEC presente en la atmósfera durante el vuelo. En el presente trabajo se incluye una breve descripción del sistema electrónico, el proceso de selección del material utilizado en las bandejas, la distribución de los componentes en cada bandeja, la topología del conjunto diseñada con ayuda de programas de diseño asistido por computadora (CAD), la fabricación propiamente dicha y los pasos seguidos en el montaje de las diferentes bandejas de un transmisor. Es importante resaltar que este trabajo forma parte del desarrollo del primer instrumento peruano lanzado al espacio, y que su diseño y construcción ha estado en su totalidad en manos de ingenieros y científicos peruanos. DISEÑO DE LAS BANDEJAS Descripción del Transmisor PERSEUS Los elementos de un transmisor del proyecto PERSEUS son los siguientes: • Bandeja de Osciladores • Bandeja de Amplificadores • Bandeja Adaptadora de Impedancias • Dos antenas de VHF • Dos antenas de UHF • Cables de interconexión. La bandeja de osciladores genera una frecuencia de 37.866MHz con un nivel de 10 dBm, y otra frecuencia de 567.99MHz con un nivel de 15 dBm. Ambas son acondicionadas en la Bandeja de Amplificadores obteniéndose 40 dBm de potencia de salida para ambas frecuencias. La señal de salida VHF es entregada a la Bandeja Adaptadora de Impedancias, para luego ser conectada a las dos antenas de VHF. La señal de UHF es llevada a las antenas de UHF a través de una red de adaptación de impedancias. La descripción funcional y de diseño de los transmisores está contenida en el informe técnico del proyecto PERSEUS (Sarango et al., 2004). Desarrollo de contenedores para los transmisores del Proyecto Espacial Perseus I Figura 1. Diagrama de conexiones del transmisor PERSEUS Elección de Materiales El aluminio fue el material seleccionado para la fabricación de las bandejas de los transmisores, debido a que posee las siguientes características: • Es liviano • Resistente • Excelente conductor de calor y electricidad • Dúctil y maleable. Distribución de Componentes Bandeja de Osciladores Está compuesta por los siguientes elementos (Figura 2):  1 OCXO (Oscilador controlado por temperatura) de 37.866 MHZ  1 Divisor de potencia  1 Amplificador de bajo ruido  1 Multiplicador de frecuencia x 5  1 Multiplicador de frecuencia x 3  1 Regulador de Voltaje J. Chocos & K. Kuyeng  2 Conectores SMA  1 Conector DEMA 9 Pines  8 Adaptadores SMA  Cable coaxial semirígido RG-402 Para la distribución de componentes se tomó en cuenta criterios de disipación de calor y las limitaciones de espacio que se tenían. Figura 2. Distribución de Componentes en la Bandeja de Osciladores Se utilizaron dos tipos de osciladores: tipo ULNO y tipo Sprinter, siendo la diferencia que el primero tiene un nivel de ruido de fase muy pequeño, del orden de los -165 dBc. Bandeja de Amplificadores Está compuesta por los siguientes elementos (Figura 3):  1 Amplificador de Potencia UHF  1 Amplificador de Potencia VHF  1 Conversor DC-DC  1 Circuito Impreso para el Conversor  2 Circuitos Impresos para los amplificadores  1 Conector DEMA 9 Pines  4 Conectores SMA Hembra  Cable coaxial semirígido RG-402 Los amplificadores se ubicaron centrados en la bandeja para que el área de disipación sea mayor. Figura 3. Distribución de Componentes en la Bandeja de Amplificadores de Potencia Bandeja Adaptadora de Impedancias Compuesta por los siguientes elementos (Figura 4):  3 Conectores SMA Hembra de Panel  1 Grillete de cables eléctricos  Cable coaxial RG-179 (75 ohms)  Cable coaxial RG-400 (50 ohms) Figura 4. Bandeja de adaptación de impedancias Con la finalidad que los cables de la red de adaptación de impedancias no se enreden o estorben en el espacio asignado, se los arrolló en forma de espiral dentro de ésta bandeja. Desarrollo de contenedores para los transmisores del Proyecto Espacial Perseus I Selección de Dimensiones El espacio disponible para el transmisor fue proporcionado por NASA. El transmisor debería de montarse dentro de un cilindro de 5.56” de alto y 12” de diámetro (Figura 5). Figura 5. Espacio disponible en el cohete para el Experimento PERSEUS Tomando en cuenta el espacio disponible, el tamaño de los componentes del transmisor y las distribuciones escogidas, las bandejas se diseñaron con dimensiones similares. En el caso de la bandeja adaptadora de impedancias, la altura no es de 1.5” como las otras sino de 1.05” (Figura 6). Figura 6. Tamaño de bandejas (oscilador y amplificador de potencia): 7.7” x 6.8” x 1.5 “ Las dimensiones del transmisor completo pueden ser observado en la Figura 7. Figura 7. Transmisor completo: 7.7” x 6.8” x 4.05“ Diseño de Bandejas en Autocad Se siguieron los siguientes pasos para diseñar las bandejas en AUTOCAD: • Obtener del fabricante y/o a través de mediciones las dimensiones exactas de cada componente • Construir un modelo en AUTOCAD de cada componente • Distribuir dichos modelos en las bandejas de acuerdo a los criterios vistos anteriormente • Diseñar las rutas para los canales de cables de alimentación. • Para evitar que la vibración excesiva de los componentes se ubicaron al ras del nivel superior de las bandejas, de esta manera su movimiento estará limitado por los tornillos en las bases y por la tapa superior de la bandeja. J. Chocos & K. Kuyeng Bandeja de Osciladores ULNO Figura 8. Medidas de la bandeja de Oscilador ULNO (vista superior) Figura 9. Sólido de la bandeja de Oscilador ULNO (vista superior) Figura 10. Sólido de la bandeja de Oscilador ULNO (vista inferior) Bandeja de Osciladores Sprinter Figura 11. Medidas de la bandeja de Oscilador SPRINTER (vista superior) Figura 12. Sólido de la bandeja de Oscilador SPRINTER (vista superior) Figura 13. Sólido de la bandeja de Oscilador SPRINTER (vista inferior) Desarrollo de contenedores para los transmisores del Proyecto Espacial Perseus I Bandeja de Amplificadores de Potencia Figura 14. Medidas de bandeja de Amplificadores de potencia (vista superior) Figura 15. Sólido de bandeja de Amplificadores de Potencia (vista superior) Figura 16. Sólido de bandeja de Amplificadores de Potencia (vista inferior) Bandeja para adaptación de impedancias Figura 17. Medidas de bandeja para adaptación de impedancias (vista superior) Figura 18. Sólido de la bandeja para adaptación de impedancias Fabricación de Bandejas de Aluminio Fabricación de prototipo Cuando el diseño en AUTOCAD de las bandejas estuvo terminado, se mandó a fabricar un prototipo para: • Determinar la existencia de posibles errores de diseño • Familiarizarse con el montaje de los componentes • Realizar pruebas de transmisión de potencia en campo • Determinar si la disipación de calor era la adecuada J. Chocos & K. Kuyeng El prototipo se fabricó en Lima usando fresadoras y tornos de control manual. Fabricación de las bandejas del Transmisor Las bandejas de osciladores y de amplificadores fueron fabricadas en un taller especializado utilizando maquinas de control numérico. Esto, debido a la complejidad del instrumento y a que era necesario un buen ajuste entre las bandejas y sus componentes para evitar vibraciones excesivas. Las bandejas para adaptación de impedancias se fabricaron en una fundición debido a la simplicidad de su diseño (Figura 19). Figura 19. Bandeja de Osciladores fabricada con máquinas de control numérico Montaje de los transmisores Montaje de Bandeja de Osciladores Los puntos a resaltar en el montaje de la bandeja de osciladores son los siguientes (Figura 20): • Se verifica que cada uno de los módulos del transmisor encaje dentro de su respectiva ubicación en la bandeja. En los casos en que existan imperfecciones del maquinado se utilizan limas y esmeriles para eliminarlas. • Se perforan los huecos para los tornillos que fijan los módulos a la bandeja de aluminio. Se perforan también los huecos para los tornillos de los conectores externos y de las tapas. • Se hacen las roscas, a los agujeros taladrados en la bandeja, de acuerdo al tipo tornillo utilizado. • Se interconectan el oscilador, el divisor de potencia, el amplificador y los multiplicadores con sus respectivos conectores. • Se agrega líquido fijador de tornillos a los conectores para evitar que se desajusten por la vibración. • Se sueldan cables de suficiente longitud en los terminales de alimentación de los módulos, de manera que alcancen hasta el regulador de voltaje. • Se comprueba que no existan corto circuitos o soldaduras frías y se verifica el funcionamiento de cada módulo. Figura 20. Interconexión de módulos en la Bandeja de Osciladores Desarrollo de contenedores para los transmisores del Proyecto Espacial Perseus I • Se montan los módulos en la bandeja • Los cables de alimentación se llevan, para ser conectados, al regulador de voltaje a través de las canaletas. • Se colocan los conectores SMA de panel en las salidas VHF y UHF. • Se sueldan los cables de alimentación de voltaje no regulado en el conector DEMA • Se monta el conector DEMA en la bandeja. • Se conectan, al regulador de voltaje, los cables de alimentación de los módulos y los cables provenientes del conector DEMA. • Se rellenan todos los espacios vacíos de la bandeja con silicona RTV para reducir el efecto de las vibraciones • Por último, se colocan las tapas a las bandejas. Figura 21. Bandeja de Osciladores ya montada y sin RTV Montaje de Bandeja de Amplificadores Los puntos a resaltar en el montaje de la bandeja de amplificadores son los siguientes: • Se suelda la tarjeta de circuito impreso sobre el conversor DC-DC y luego se montan los componentes correspondientes, de forma superficial.(Figura 22) Figura 22. Montaje de componentes del circuito impreso del conversor DC-DC • Se perforan los huecos para los tornillos que fijan el conversor DC-DC, los amplificadores, los circuitos impresos, los conectores externos y las tapas. • Se hacen las roscas, a los agujeros de acuerdo al tipo tornillo utilizado. • Se sueldan los conectores SMA de 90º y los amplificadores a sus respectivas tarjetas (Figura 23). • Se sueldan cables en los terminales de alimentación de los amplificadores, de manera que alcancen hasta el conversor DC-DC • Se comprueba que no existan corto circuitos o soldaduras frías. Figura 23. Detalle de conexión de cables y conectores a los amplificadores J. Chocos & K. Kuyeng • La parte inferior del amplificador se embadurna con pasta térmica de silicona para permitir una mejor conducción de calor entre el amplificador y la bandeja de aluminio. • Se coloca líquido fijador en los conectores y tornillos para evitar que se desajusten por la vibración. • Se montan los amplificadores en su respectiva ubicación (Figura 24). Figura 24. Detalle del montaje de un amplificador • Se colocan los conectores SMA de panel en las salidas VHF y UHF (Figura 25). Figura 25. Detalle del montaje de los conectores de panel • Se conectan, al conversor DC-DC, los cables de alimentación de los amplificadores y los cables provenientes del conector DEMA. • Se rellenan todos los espacios vacíos de la bandeja con silicona RTV para reducir el efecto de las vibraciones (Figura 26) Figura 26. Bandeja de Amplificadores con RTV de protección • Se colocan las tapas rotuladas (Figura 27). Figura 27. Bandeja de Amplificadores terminada y con tapas Montaje de la bandeja para adaptación de impedancias Para diseñar las longitudes de los cables que componen esta bandeja se empleó la Carta de Smith. A cada antena se le agregó una longitud “L” de cable coaxial RG-179 Desarrollo de contenedores para los transmisores del Proyecto Espacial Perseus I de 75 Ω de tal manera que la impedancia de 50Ω es transformada en 100Ω. Luego, al estar las dos antenas en paralelo, la impedancia resultante son los 50Ω que queríamos. A una de las antenas se le agregó media longitud de onda (λ /2) para balancear el flujo de corrientes. El diagrama de conexión de la red de adaptación de impedancias es la siguiente (Figura 28): Figura 28. Conexiones de la Red de Adaptación de Iimpedancias Los puntos a resaltar en el montaje de la bandeja de amplificadores son los siguientes: • Los cables se acomodaron en el interior de la bandeja arrollados en forma de espiral. • Se usaron pernos especiales, llamados grilletes, para facilitar la conexión de las mallas de tierra de los cables. • Se agregó silicona RTV a los espacios vacíos para minimizar las vibraciones RENDIMIENTO Y RESULTADOS • Una vez construidos los transmisores, fueron sometidos a diversas y exigentes pruebas de pre-vuelo Los resultados obtenidos fueron satisfactorios. • Como ya se ha mencionado, se construyeron ocho transmisores, pero por motivos de orden logístico, sólo seis fueron lanzados al espacio. • Los seis instrumentos funcionaron perfectamente y sin inconvenientes. • Durante el vuelo de uno de los cohetes se presentaron problemas en el sistema de seguimiento de NASA. Sin embargo, el transmisor PERSEUS funcionó correctamente, permitiendo determinar la trayectoria del cohete. Esta información ha sido entregada a NASA como parte de la colaboración IGP-NASA. CONCLUSIONES Varios miembros del equipo peruano participaron en las pruebas de pre-vuelo y en la integración de los transmisores en los cohetes. Estas actividades se realizaron en los laboratorios de NASA Wallops Flight Facility en Virginia, Estados Unidos. Y se constituyeron en una oportunidad inmejorable de entrenamiento y capacitación. La experiencia y los conocimientos adquiridos en este proyecto son invaluables: Se tuvo la oportunidad de incursionar en el desarrollo de tecnología aerospacial, afianzar los conocimientos de electrónica y radiofrecuencia, interactuar con científicos y técnicos de la NASA, visitar laboratorios y universidades en los Estados Unidos entre otras actividades (Figura 29). J. Chocos & K. Kuyeng Figura 29. Ingenieros del ROJ en los laboratorios de NASA mostrando un transmisor en el interior de uno de los cohetes Los resultados científicos del experimento serán comunicados en publicaciones posteriores. Sin embargo, los análisis preliminares indican que el experimento PERSEUS fue exitoso, tanto en su componente tecnológica como en la parte científica. El éxito del experimento ha abierto las puertas para la realización de otros proyectos futuros en el campo de la tecnología aeroespacial en el Perú. AGRADECIMIENTOS Nuestra gratitud al Dr. Martín Sarango, Jefe del Proyecto PERSEUS, por confiar en nosotros y por su dirección y apoyo durante la realización del presente trabajo. Asimismo, al personal del ROJ por su ayuda y comprensión. BIBLIOGRAFÍA ARRL (2002). The ARRL Handbook for Radio Amateurs, 69 th edition. Sarango, M.F., R.F. Woodman, J.L. Chau, J. Chocos, C. De La Jara, K. Kuyeng, G. Michhue y C. Minchola (2004): Informe Técnico del Proyecto PERSEUS I, Reporte Técnico Interno, Radio Observatorio de Jicamarca - IGP, Diciembre. Compendio de Trabajos de Investigación CNDG – Biblioteca Instituto Geofísico del Perú. V. 5 (2004) p. 61 - 68 AUTOMATIZACIÓN DE LOS TRANSMISORES DE ALTA POTENCIA DEL RADIO OBSERVATORIO DE JICAMARCA – SEGUNDA ETAPA JUAN CARLOS ESPINOZA GUERRA Programa de Ingeniería Mecánica Eléctrica Facultad de Ingeniería Universidad de Piura jespinoza@jro.igp.gob.pe Prácticas dirigidas por: Ing. Otto Castillo Gonzáles Dr. Justo Oquelis Cabredo Radio Observatorio de Jicamarca - IGP RESUMEN El presente proyecto tiene por objetivo principal la ejecución de la segunda etapa de automatización y control de los transmisores de alta potencia del Radio Observatorio de Jicamarca. La primera etapa esta ya concluida y ahora el sistema de transmisores del ROJ cuenta con un PLC (Program Logic Controller) y un sistema SCADA (Supervisory, Control and Data Adquisition) como elementos de control y supervisión respectivamente. En la segunda etapa del proyecto se tiene previsto automatizar algunos sistemas que aún se operan de forma local y manual como son el encendido de las etapas PA y DRIVER de transmisión, el control de la tensión de filamento de la etapa PA y el control de la tensión de pantalla de las etapas PA y DRIVER. La segunda etapa también involucra la incorporación al sistema de unos módulos inteligentes de la marca Acromag los cuales poseen entradas analógicas y capacidad de comunicación con el PLC vía protocolo MODBUS estándar RS485. INTRODUCCIÓN El sistema de control de los transmisores de alta potencia del Radio Observatorio de Jicamarca (ROJ) fue concebido para operar con tecnología electromecánica, su funcionamiento y lógica de control se hacía por medio de contactores y relés. Principalmente, el control se realizaba sobre los parámetros del sistema de enfriamiento de agua, tanto de alta como de baja presión; así como sobre algunos parámetros eléctricos del mismo. El ROJ, siempre está innovando y actualizando su tecnología de funcionamiento, por tanto con el fin de lograr este propósito se puso en marcha la automatización con PLC del sistema de transmisores. Actualmente el sistema de transmisores de alta potencia del ROJ cuenta con un sistema de control automático y supervisión para los siguientes procesos: - Activación de la alimentación primaria de baja tensión (120V) - Sistema de enfriamiento el cual incluye dos circuitos de agua uno de baja y otro de alta presión, cuatro intercambiadores de calor agua – aire y dos ventiladores auxiliares. - Selección de los transmisores a utilizar. - Activación de las fuentes de alta tensión (20kV) para las etapas PA de cada transmisor. - Activación de la señal de RF. - Sistema de alarmas. J. Espinoza 62 Con la ejecución de la segunda etapa el sistema de transmisores del ROJ quedará totalmente automatizado con lo que se espera mejorar su operación. FUNDAMENTOS Cuando el ROJ opera con su máxima potencia (4x1.5MW) la señal que es enviada a la antena pasa por diferentes etapas de amplificación como se muestra en la Figura 1, durante la ejecución de la segunda etapa del proyecto de automatización se trabajará básicamente en las etapas PA y DRIVER de transmisión. 20KV Placa Pantalla Grilla Filamento PA ANTENAEXCITADORES Filamento Grilla Pantalla Placa DRIVER 10KV 8W 10KW 100KW 1.5MW 1 - 5W Filamento Pantalla Placa Grilla MST Modificado 16.5V14.8V 900V 480V 900V 300V Figura 1. Etapas de Amplificación. Etapa PA La etapa de salida de los transmisores del ROJ, tienen como elemento principal un tubo de vacío de alta potencia tetrodo modelo 8973 de la marca Eimac, el cual es usado como amplificador de RF, este tubo cuenta con cuatro electrodos: el filamento (cátodo), la placa (ánodo), control de grilla y la grilla de pantalla. Los parámetros controlados en esta etapa corresponden a la tensión de filamento, tensión de grilla y tensión de pantalla. De los cuales la tensión de filamento y la tensión de pantalla son variables. En cuanto la tensión de grilla es fija. El panel de control de la etapa PA cuenta con dos juegos de pulsadores ON/OFF, el primero para el encendido y apagado del sistema y el segundo para el encendido y apagado de la tensión de pantalla del tubo de radiofrecuencia; también cuenta con dos selectores para ajustar la tensión de filamento; un variac para ajustar la tensión de grilla y finalmente otro variac para ajustar la tensión de pantalla. Tanto los controles de encendido y apagado como las tensiones de filamento y de pantalla deberán ser controlados de manera remota, utilizando el PLC para su control, además de funcionar de manera local como se ha venido haciendo hasta ahora. En el caso de la tensión de grilla, esta no requiere de control automático, ya que trabaja en un valor fijo. Etapa DRIVER La etapa DRIVER también cuenta con un tubo de vacío el 4CX40,000G. El cual posee cuatro electrodos: el cátodo o filamento, la grilla de control, la grilla de pantalla y el ánodo o placa. La operación de la etapa DRIVER está definida por la tensión de pantalla, la cual dependiendo del nivel de potencia que se quiera transmitir se ajusta entre los valores a elegir 0 y 500 VDC y variable hasta 750 VDC como máximo. Principalmente la automatización se realizará sobre el panel de control de esta etapa, el cual cuenta con dos juegos de pulsadores ON/OFF para el encendido y apagado del sistema y de la tensión de Automatización de los transmisores de alta potencia del ROJ 63 pantalla respectivamente además de un variac para regular la tensión de pantalla. Estos controles deberán ser controlados de manera remota, utilizando para su control el PLC. Además de funcionar de manera local como se hace actualmente. DESCRIPCION DEL PROYECTO La segunda etapa del proyecto de automatización de los transmisores principales, en un principio comprende los siguientes puntos: - Control de los pulsadores de encendido y apagado de las etapas PA y DRIVER. - Control de la tensión de filamento de la etapa PA. - Control de la tensión de pantalla de las etapas PA y DRIVER. - Implementación de un sistema de medición de los parámetros eléctricos de la red principal como son la tensión, corriente y potencia primaria. Además de la potencia de salida de los transmisores. - Incorporación al sistema de unos módulos I/O (input/output) inteligentes de la serie 900MB de la marca Acromag con interfase de comunicación vía protocolo MODBUS estándar RS-485. Control de los pulsadores ON/OFF de la etapa PA Para el control de los pulsadores de encendido y apagado del sistema y de la tensión de pantalla es necesario utilizar relés auxiliares los cuales serán activados y desactivados utilizando un módulo de salidas tipo relé del PLC. Es necesario utilizar relés auxiliares para conmutar los nuevos circuitos con los circuitos existentes, esto para que sea posible la operación de manera local y sin utilizar el PLC. Estos relés auxiliares serán activados por un selector LOCAL/REMOTE ubicado en el panel de control. Por seguridad los pulsadores OFF del panel de control podrán ser activados siempre sin importar la posición del selector LOCAL/REMOTE. Control de los pulsadores ON/OFF de la etapa DRIVER Al igual que la etapa PA se utilizaran relés auxiliares los cuales serán activados por un módulo de salidas tipo relé del PLC, también se utilizaran relés para conmutar los circuitos nuevos de los existentes, esto por la misma razón que en la etapa PA. Por seguridad los pulsadores OFF del panel de control podrán ser activados siempre sin importar la posición del selector LOCAL/REMOTE. Control de la tensión de filamento de la etapa PA La tensión de filamento es proporcionada por una fuente, la cual cuenta con un transformador variable que tiene como tensión de entrada 440 V trifásicos y como salida una tensión que puede ir de 0 V a 20 V, también trifásicos. Esta tensión de salida es rectificada y conectada al cátodo (filamento) del tubo de vacío del transmisor. J. Espinoza 64 El transformador variable de la fuente de la tensión de filamento esta acoplado a un motor síncrono de control bidireccional, con el cual es posible ajustar la tensión de filamento haciéndolo girar en un sentido ó en otro. El control del giro del motor síncrono se realiza utilizando dos selectores del panel de control. Un primer selector de dos posiciones Auto/Manual es utilizado para incrementar de forma automática la tensión de filamento hasta un límite prefijado. El segundo selector de tres posiciones Raise/Neuter/Lower, es utilizado para incrementar ó disminuir la tensión de filamento, cuando el selector Auto/Manual se encuentra en la posición Manual. El control de la tensión de filamento actualmente se basa en la activación de unos contactores que impiden que la tensión de filamento supere los límites preestablecidos, quitando la alimentación del motor síncrono. En el nuevo sistema, el control de la tensión de filamento se efectuará utilizando el PLC, el cual por medio de un módulo de salidas tipo relé activará o desactivará unos relés auxiliares que alimentan el motor síncrono. Es necesario utilizar relés para conmutar los circuitos nuevos y los existentes, ya que el sistema debe permitir la operación local (como se viene realizando actualmente) La retroalimentación del sistema se realizará utilizando los circuitos acondicionadores existentes, que adaptan el nivel de la tensión de filamento a un nivel que puede ser leído por un módulo de entradas analógicas del PLC (Figura 2). RELES 440 V TENSIÓN DE FILAMENTO 0 - 20V CIRCUITO ACONDICIONADOR0 - 5V MOTOR SÍNCRONO TRANSFORMADOR VARIABLE DL305 DL35 0 A /D IN PU T D IG ITA L O U TPU T Figura 2. Diagrama de bloques del sistema de control de la tensión de filamento de la etapa PA Control de la tensión de pantalla de las etapas PA y DRIVER La tensión de pantalla de la etapa PA es proporcionada por una fuente DC, la cual entrega una tensión variable de 0 VDC a 1350 VDC, que se ajusta utilizando un transformador variable en su entrada AC. El valor nominal de operación de la tensión de pantalla es de 900 V. La tensión de pantalla en la etapa DRIVER, también es proporcionada por una fuente DC, que entrega una tensión variable entre 0 VDC y 1200 VDC. Para ajustar esta tensión también se utiliza un transformador variable. Los valores de operación de la tensión de pantalla en la etapa DRIVER pueden estar en el rango de 500 VDC a 750 VDC, dependiendo de la potencia de salida que se requiera. El elemento a controlar es el transformador variable con el que se ajusta la tensión de pantalla, tanto de la etapa PA como de la etapa DRIVER. Para su control se utilizará un motor paso a paso el cual a Automatización de los transmisores de alta potencia del ROJ 65 su vez será controlado por una estación remota (Figura 3). El PLC se comunicará con la estación remota enviándole los datos necesarios para hacer girar los motores a pasos. La comunicación entre PLC y estaciones remotas se basará en el estándar RS485 protocolo MODBUS en donde el PLC actuará como maestro y las estaciones remotas como esclavos. DL305 DL350 PC Convertidor RS232/RS485 RS232 R S2 32 PLC ADAPTADOR DE COMUNICACIONES PIC16F873 DRIVER1 ESTACION REMOTA TX1 BUS RS485 DRIVER2 POWER SUPPLY MOTOR PA MOTOR DRV ADAPTADOR DE COMUNICACIONES PIC16F873 DRIVER1 ESTACION REMOTA TX2 DRIVER2 POWER SUPPLY MOTOR PA MOTOR DRV ADAPTADOR DE COMUNICACIONES PIC16F873 DRIVER1 ESTACION REMOTA TX3 DRIVER2 POWER SUPPLY MOTOR PA MOTOR DRV ADAPTADOR DE COMUNICACIONES PIC16F873 DRIVER1 ESTACION REMOTA TX4 DRIVER2 POWER SUPPLY MOTOR PA MOTOR DRV Figura 3. Diagrama de bloques del sistema de control de la tensión de pantalla La estación remota será la encargada de controlar el giro de los motores a pasos que ajustarán la tensión de pantalla de las etapas PA y DRIVER al valor adecuado en cada caso (Figura 4). La estación remota tendrá como elemento principal el circuito integrado PIC16F873A, además contará con un adaptador de comunicaciones, las respectivas etapas de potencia para alimentar los motores a pasos y su respectiva fuente. DL305 D L350 A/D PIC 16F873MAX485 BUS RESET POWER SUPPLY ADAPTADOR DE COMUNICACIONES 5VDC5VDC V+ 110V TENSIÓN DE PANTALLA PA 0 - 1300 VDC MOTOR PAP VARIAC V+ A B C D a b c d CIRCUITO ACONDICIONADOR DRIVER V+ 110V MOTOR PAP VARIAC V+ A B C D a b c d DRIVER ESTACION REMOTA PA DRIVER TENSIÓN DE PANTALLA DRIVER 0 - 1500 VDC CIRCUITO ACONDICIONADOR 0 - 5VDC 0 - 5VDC PLC 110V Figura 4. Diagrama de bloques del sistema de control de la tensión de filamento de la etapa PA EL microcontrolador PIC 16F873 cuenta con 4096 palabras de memoria de programa, 192 bytes de RAM, 128 bytes de flash EEPROM y 22 I/O de propósito general, además de una interfase de comunicaciones USART (Universal Addressable Synchronous Asynchronous Receiver Transmitter) implementada en hardware [Microchip, 2003]. El programa principal del PIC constara de un lazo cerrado que compara la posición actual del motor con la enviada por el PLC en caso de ser diferente, la posición del motor es actualizada (Figura 5). INICIO CONFIGURACION DE PUERTOS Y USART RESET MOTOR_PA Y MOTOR_DRV MOTOR_PA = MOTOR_PA_RX MOTOR_DRV = MOTOR_DRV_RX SI SI MAIN ACTUALIZA MOTOR_PA ACTUALIZA MOTOR_DRV NO NO Figura 5. Diagrama de flujo del programa principal del PIC J. Espinoza 66 INTERRUPCCION RECEPCION DE DATOS CODIGO MODBUS OK EJECUTA ACCION QUE INDICA EL CODIGO MODBUS ERROR EN COMUNICACION SI NO SI DATOS MODBUS OK NO NO SI ENVIA RESPUESTA SALE DE INTERRUPCION ENVIA ERROR Figura 6. Diagrama de flujo de la rutina de interrupción del PIC Cuando el PIC recibe datos del PLC, el programa principal se interrumpe y se atiende la rutina interrupción que recibe los datos, ejecuta la acción que indica el código y envía la respuesta al PLC. Las rutinas de actualización de los motores consisten en la generación de los pulsos necesarios para hacerlo girar en un sentido ó en el otro hasta alcanzar la posición indicada por el PLC. Como se menciono antes la comunicación a nivel físico se basará en el estándar RS485 y el PIC cuenta con una interfase para comunicaciones seriales USART que opera con niveles TTL por lo que para poder conectar el PIC al bus de comunicaciones RS485 es necesario adaptar las señales TTL a las utilizadas en el estándar RS485, para ello se utiliza un transceptor 485. Este integrado contiene un receptor diferencial con histéresis y un transmisor, también diferencial, que puede ser puesto en estado de alta impedancia. Esto permite que varios transceptores estén conectados a una sola línea de bus formada por 2 conductores (A y B) además de un blindaje. La señal de un uno lógico se traduce a que la línea A sea más positiva que la línea B en un valor mayor a 200mV y en el caso contrario B>A tendremos un cero lógico [EIA, 1997]. El protocolo MODBUS utilizado por el PLC es un protocolo maestro esclavo, donde solo el maestro puede iniciar la comunicación con cualquiera de los esclavos. El protocolo MODBUS define un simple PDU (Protocol Data Unit) que consta de dos campos la función y los datos, cuando el protocolo trabaja en un bus ó una red se introducen campos adicionales sobre el ADU (Application Data Unit) (Figura 7). DIRECCION FUNCION DATOS CRC ADU PDU Figura 7. General MODBUS frame En donde la dirección (1 byte) identifica al esclavo con el que se va a iniciar la comunicación, puede tomar valores desde 1 hasta 247 decimal. La función (1 byte) indica la acción a ejecutar, por ejemplo leer un registro o escribir en un registro del esclavo. Los datos (2xN bytes) son los Automatización de los transmisores de alta potencia del ROJ 67 bytes necesarios para ejecutar la acción que indica el código y por último el CRC (Cyclic Redundancy Check) constan de dos bytes enviados para comprobar si hubo errores en la comunicación [Modicon, 2002]. RESULTADOS Y DISCUSIÓN El proyecto se encuentra en su etapa de desarrollo e implementación y hasta ahora se han realizado los siguientes avances. Control de los pulsadores ON/OFF de las etapas PA y DRIVER Para determinar el control de estos pulsadores se estudió la secuencia detallada de encendido y apagado tanto de la etapa PA como de la etapa DRIVER. Se modificaron los planos de los circuitos de control de ambas etapas en donde se muestran los nuevos circuitos a implementar. También se identificaron las nuevas entradas y salidas del sistema, en total se determino que son necesarias 44 salidas y 32 entradas, lo que corresponde a tres módulos de 16 salidas tipo relé (D3-16TR) y dos módulos de 16 entradas digitales de 110V (D3-16NA) del PLC. Control de la tensión de pantalla de las etapas PA y DRIVER En esta parte del proyecto se trabajó básicamente con el microcontrolador PIC, el cual es el elemento principal de la estación remota que controla la tensión de pantalla. Lo primero que se realizo fue un circuito grabador con el cual es posible programar el PIC. Luego se trabajo en la implementación del protocolo MODBUS en el PIC, hasta ahora se tiene una primera versión del programa el cual soporta las funciones de lectura y escritura de registros en la RAM del PIC, pudiendo así cambiar el valor de las variables utilizadas en el programa de control de los motores. El programa también detecta errores en la comunicación ya que calcula los bytes del CRC y genera respuestas de error para las funciones no soportadas así como también por registros y datos erróneos. Para probar el programa del PIC se monto una pequeña red utilizando una PC como estación maestra y dos microcontroladores PIC como estaciones remotas, dado que la PC cuenta con un puerto de comunicaciones RS232 fue necesario utilizar un convertidor auto soportado RS232/RS485. En la PC se utilizo el software LookOut (software SCADA) como Terminal MODBUS que genera los códigos y los envía al PIC Para simular la posición de los motores a pasos se utilizaron juegos de led’s. El sistema utilizado en las pruebas se muestra en la Figura 8: J. Espinoza 68 PC Convertidor RS232/RS485 R S2 32 R S4 85 PIC 16F873 MAX 485 LED’S PIC 16F873 MAX 485 LED’S Figura 8. Sistema utilizado en las pruebas con el PIC CONCLUSIONES Las pruebas realizadas con el microcontrolador indican que es posible montar una red de estaciones remotas capaces de comunicarse con una estación maestra utilizando el protocolo MODBUS. El siguiente paso es montar la red pero utilizando el PLC como estación maestra luego se procederá con el diseño de la estación remota. Se está estudiando la posibilidad de utilizar los módulos Acromag como interfase entre el PLC y la estación remota en lo que se refiere al control de la tensión de pantalla. Esto debido a que los módulos ya cuentan con el protocolo de comunicación implementado, además de cuatro salidas digitales con las cuales se podría enviar las ordenes al microcontrolador. AGRADECIMIENTOS En primer lugar un agradecimiento muy especial al Dr. Jorge L. Chau por brindarme la oportunidad de realizar mis practicas en el Radio Observatorio de Jicamarca. Asi mismo, agradezco al Dr. Justo Oquelis y al Ing. Otto Castillo por su paciencia y constante apoyo. BIBLIOGRAFÍA EIA, (1997) “EIA-485 Standard for electrical characteristics of generators and receivers for use in balanced digital multipoint systems”. Web: http://global.his.com. Microchip, (2003). “PIC16F87XA Data Sheet”. Web: http://www.microchip.com. Modicon, (2004). “MODBUS Application Protocol Specification”, version 1.1a. Web: http://www.modbus.org. Modicon, (2002) “MODBUS over a Serial Line Specification & Implementation guide”, version 1.0. Web: http://www.modbus.org. Compendio de Trabajos de Investigación CNDG – Biblioteca Instituto Geofísico del Perú. V. 5 (2004) p. 69 - 76 ESTUDIO DE LOS ECOS DE METEOROS DISCRIMINANDO LA SEÑAL DEL ELECTROCHORRO ECUATORIAL FREDDY RONALD GALINDO PALOMINO Escuela Profesional de Ingeniería Física Facultad de Ciencias Universidad Nacional de Ingeniería fgalindo@jro.igp.gob.pe Prácticas dirigidas por: Dr. Jorge Chau Radio Observatorio de Jicamarca - IGP RESUMEN Las observaciones de los ecos de meteoros datan del año 1940 aproximadamente, mas estas observaciones empiezan a ser un tópico de investigación recién en los últimos años. La búsqueda de una caracterización del meteoro en parámetros físicos tales como la velocidad, rango de altura o dirección, se desarrolló haciendo uso del radar como herramienta de adquisición de información de estos eventos; el cual conjuntamente con un procesamiento de la señal permite extraer la información requerida en el estudio. El presente trabajo tendrá como objetivo determinar la velocidad radial de los meteoros, cuya señal adquirida será afectada por la influencia del Electrochorro Ecuatorial (Ecuatorial Electrojet - EEJ) INTRODUCCIÓN El estudio de ecos de meteoros en la atmósfera limeña es una investigación que se viene realizano en el Radio Observatorio de Jicamarca desde hace algunos años [Chau y Woodman, 2004]. Los experimentos permiten almacenar información de una gran cantidad eventos para su posterior procesamiento y extracción de los parámetros físicos a investigar (Figura 1). La dificultad que presenta el Radio Observatorio en la adquisición de ecos de meteoro, se debe a una fuerte presencia de ecos debido al Electrochorro Ecuatorial [Yerson et al., 2004] presente sobre Jicamarca. Debido a esa desventaja las observaciones se limitan a aquellos momentos en el cual la influencia del Electrochorro Ecuatorial es muy débil o se presenta esporádicamente. Este trabajo pretende generar un programa de procesamiento que permitirá discriminar la señal generada por el Electrochorro Ecuatorial en la señal recibida, para poder estudiar solamente la señal generada por el meteoro, lo cual permitiría un mayor tiempo de investigación. Para fines de estudio y análisis, se hará uso del modelo generado por el Dr. J. Chau, el cual simula ecos de meteoro sin la presencia del Electrochorro Ecuatorial. Esto servirá para comprender el funcionamiento del algoritmo de procesamiento de la señal. Posteriormente, se desarrollará un modelo que tenga en cuenta la influencia del Electrochorro Ecuatorial, lo cual nos permitirá validar el procesamiento desarrollado para trabajar la señal. Así también se tomaran como base los programas desarrollados por el Dr. J. Chau para el procesamiento de señales de meteoros. F. Galindo 70 Figura 1. Gráfico que ilustra la ubicación del meteoro en función de la altura, Notemos que el evento empezó a ser estudiado a una altura de 110Km. FUNDAMENTOS El funcionamiento de un radar, en términos generales, consiste de un proceso de transmisión de un pulso de energía, y la recepción o “escucha” de un eco [Sarango, 1996]. Esta recepción o “escucha” se debe a la reflexión que sufre parte de la energía del pulso producto de un obstáculo (o fuente del eco) presente en su camino (Figura 2). La distancia de la fuente que produce el eco, se puede determinar calculo y el tiempo transcurrido entre la transmisión y recepción. Aunque los radares tuvieron sus primeros usos en la búsqueda de objetos sólidos, su rango de utilidad se amplio, siendo utilizados hoy, en el estudio de “hidrometeoros”. Figura 2. Esquema ilustrativo de la operación de un radar. Notemos que el pulso de energía transmitido será reflejado por el obstáculo y su la energía reflejada será recibida como una señal de referencia. La información adquirida (Figuras 3 y 4) para este trabajo fue obtenida mediante el uso del radar de Apertura-50Mhz del Radio Observatorio de Jicamarca, usando en la transmisión un pulso de energía codificado. El código binario utilizado es el código Barker-13 [Brookner, 1977], el cual permite lograr buena resolución sin perder sensibilidad, además la fase del código (0-180°) es implementado fácilmente. En general el código Barker (Figura 5) se caracteriza por tener, un pico de respuesta igual a N, donde N es la longitud del código. El signo menos indica una fase de cero grados y el signo positivo una de 180 grados. Estudio de los Ecos de meteoros discriminando la señal del electrochorro ecuatorial 71 Figura 3. Señal recibida en ausencia de un objeto que refleje el pulso de energía emitido. Figura 4. Señal recibida en presencia de un objeto que refleje el pulso de energía emitido. Figura 5. Gráfico de la auto-correlación del código Barker-13, notemos la diferencia del pico central respecto a los demás. El voltaje complejo (eco) recepcionado mediante el radar, es almacenado como datos crudos (Raw Data) para su posterior procesamiento. El proceso de decodificación de la información almacenada, se realiza mediante la convolución de la señal con el código, siendo una forma práctica y rápida de implementar este proceso, mediante el uso la transformada de Fourier. A continuación se presenta el algoritmo matemático de la decodificación: }}{}{{ }{}{}{ *1 * CodeFSFS CodeFSFSF CodeSS inout inout inout F −= = ⊗= donde, F{ } simboliza la transforma de Fourier, F-1{ } la transformada inversa de Fourier, ⊗ representa el proceso de convolución, * representa la conjugada F. Galindo 72 compleja, Sin la señal de entrada, Sout la señal obtenida luego de la decodificación y Code, es el código usado en el proceso de decodificación, es nuestro caso es el Código Barker. Debido a un cambio Doppler en la frecuencia de la señal recibida, el cual se genera por la reflexión del pulso de energía sobre un objeto en movimiento, sea acercando o alejándose al radar (Doppler positivo y negativo respectivamente) el proceso de decodificación no puede realizarse directamente entre la señal y el código debido al desfase. Una forma óptima de solucionar es hacer un barrido de todas las posibles frecuencias doppler. La que genere la mejor decodificación es la frecuencia Doppler mas aproximada. Una vez establecida la frecuencia Doppler, se procede a determinar el rango de altura del evento que provoco el eco. Para esto se analiza el espectro de la señal en el rango de las alturas (Figura 6), tal como se mostramos a continuación. Figura 6. Gráfico del Espectros de Potencia de la Señal decodificada. El eje “x” representa el rango de alturas y el eje “y” el valor del espectro de potencia de la señal. Aquella región que presente el mayor pico de espectro, es la altura a la cual se dio el evento. Notemos en la Figura 6 que el pico se da en la posición 27 la cual representa una altura de 109.5 Km. RESULTADOS Y DISCUSION El procesamiento estudiado permite ser una etapa de inicio, en la introducción al procesamiento de la señal recibida. Se observaron distintos eventos de meteoros simulados por computadora, lo cual nos permite tener un conocimiento de las características que dominan el procedimiento. Por ejemplo tenemos, el espectro de la señal para obtener la frecuencia Doppler o el rango de altura, además podemos establecer un valor de voltaje en presencia y ausencia de la fuente del eco, con fines de un posterior ajuste. Determinado el valor de la frecuencia Doppler, se desarrollo un programa, que busca refinar los resultados obtenidos, por medio de un ajuste de la señal codificada, tomando como valores de inicio los determinados el procesamiento de data. La Figura 7 nos ilustra gráficamente el ajuste de la señal [Mathews et al., 2003], notemos la presencia de dos señales, debido a la representación compleja de la señal. El ajuste se ha graficado dentro de mismo gráfico. El problema principal encontrado en este proceso se debe a la estabilidad del algoritmo, ya que un mal estimado en los valores de inicio nos da resultados Estudio de los Ecos de meteoros discriminando la señal del electrochorro ecuatorial 73 erróneos. Así el factor mas importante en el ajuste es la frecuencia Doppler, la amplitud se puede estimar del valor máximo de la señal y la fase de cada señal no afecta al ajuste. Otro punto dentro de estos errores se da cuyo ajustamos una señal que es propiamente ruido, mas este problema se resuelve de manera rápida, filtryo el proceso de ajuste a señales que contengan ecos. Las Figuras 7 y 8 ilustran estas ideas. Citaremos a continuación algunos valores numéricos de la frecuencia luego del proceso de ajuste de la señal. Frequencia Inicial Frecuencia ajustada 10937.5 hz 11232.4 hz 9375.00 hz 10400.4 hz 10937.5 hz 10261.0 hz 9375.00 hz 10109.7 hz Figura 7. Gráfico del Espectro de Potencia de la Señal decodificada. El eje “x” representa el rango de alturas y el eje “y” el valor del espectro de potencia de la señal. Figura 8. El ajuste de la señal, puede generar errores, tal como se muestra a continuación. Notemos que la señal ajustada es propiamente ruido. Determinada la frecuencia Doppler. Se puede obtener la velocidad radial del meteoro mediante la siguiente relación: 2 * λ doppleroptima fV = donde, Voptima representa la velocidad radial del meteoro, Fdoppler la frecuencia F. Galindo 74 Doppler y λ la longitud de onda de la señal transmitida. Las Figuras 9 y 10 son esquemas que muestran gráficamente las velocidades obtenidas luego del ajuste realizado para un conjunto eventos. La línea continua, representan los valores obtenidos antes del ajuste, es decir son los valores calculados del espectro de la señal. Las cruces representan los valores de ajustados. Notemos que existen algunas partes donde no existen cruces o la línea cae a cero, en este caso el algoritmo filtra aquellas señales que tiene una relación señal-ruido (SNR) demasiada baja, por eso que estas regiones no se realiza ajuste. Figura 9. Gráfico de las velocidades radiales de los meteoros estudiados, la línea continua simboliza la velocidad obtenida del análisis del espectro y las cruces la velocidad después del procesote ajuste. Una segunda forma de afinar los resultados, es realizó un ajuste de la señal en el dominio de frecuencia (Figura 11). Figura 10. Gráfico similar al de la Figura 9, la diferencia es el nivel del filtro señal-ruido establecido. Figura 11. El gráfico de la izquierda es el espectro de la señal. La gráfica de la derecha es el ajuste realizado alrededor del máximo. Una comparación cuantitativa entre el ajuste realizado en el dominio de tiempo y del dominio de frecuencia, parecen no mostrar mayor diferencia. En otras palabras es posible afinar los resultados en cualquiera de los dominios. Una etapa que también se ha desarrollado es un programa de filtrado (Filtro Pasa Alta), la Figura 12 muestra la utilización del filtro, para eliminar componentes de baja frecuencia. Estudio de los Ecos de meteoros discriminando la señal del electrochorro ecuatorial 75 Figura 12. Filtrado de las componentes de baja frecuencia de la señal. CONCLUSIONES Los resultados muestran que el ajuste influir, en la búsqueda de resultados más exactos generyo valores que se aproximan a caso ideal, así también la mala estimación de la frecuencia como valor de entrada, altera el resultado obtenido, pudiendo generar valores que no están de acuerdo con el valor buscado. Por eso el ajuste debe estar limitado solamente a señales que contenga realmente ecos de meteoros, para tal fin el valor de la relación señal-ruido debe de oscilar entre 4.2 a 4.7 como valor. Dado que el trabajo todavía se encuentra en desarrollo no podemos emitir algún resultado final, más es importante mencionar que la búsqueda final del trabajo es lograr un programa óptimo, que me permita procesar información de Meteoros de manera rápida y confiable, bajo la influencia o no influencia del Electrochorro Ecuatorial. Para tal fin se están desarrollo los programas en el lenguaje de programación científica IDL 6.0, por su fácil manejo, amplia librería y libertad de trabajo. AGRADECIMIENTOS Un agradecimiento al Dr. Jorge L. Chau por brindarme la oportunidad, confianza, libertad y motivación para el desarrollo de este trabajo, a su vez por todo su apoyo y conocimientos brindados a mi persona. A todo el personal del Radio Observatorio de Jicamarca por su amabilidad y apoyo. BIBLIOGRAFÍA Brookner E. (1977): Radar Technology, 1st Edition, Artech House, pp. 135, 144-145. Chau J. y Woodman R. (2004): Observations of meteor-head echoes using the Jicamarca 50MHz radar in interferometer mode, Atmosphere, Chemistry y Physics, 4, pp. 511-521 Mathews J.D., Doherty J., Wen C. H., Briczinski S.J., Janches D., Meisel D.D. (2003): An update on UHF radar meteor observations y associated signal processing techniques at Arecibo Observatory, Journal of Atmospheric y Solar-Terrestrial Physics 65 . pp 1141- 1143 Sarango M. (1996): Sistema de controlador/procesador multi-DSP para el radar MST perfilador de vientos de la Estación Científica Antartica Machu Picchu, Tesis Doctoral, Departamento de Teoría de la Señal y Comunicaciones, Universidad Politécnica de Cataluña, pp. 0.5-0.6 F. Galindo 76 Yerson D., Anghel A., Chau J.L. y Veliz O. (2004): Daytime vertical ExB drift velocities inferred from ground-based magnetometer observations at low latitudes, Space Weather, vol. 2, pp 1-3 Compendio de Trabajos de Investigación CNDG – Biblioteca Instituto Geofísico del Perú. V. 5 (2004) p. 77 - 86 ANALISIS DINAMICO EN IMÁGENES DE LA CAPA IONOSFERICA FREDDY RONALD GALINDO PALOMINO Escuela Profesional de Ingeniería Física Facultad de Ciencias Universidad Nacional de Ingeniería fgalindo@jro.igp.gob.pe Prácticas dirigidas por: Dr. Jorge Chau Radio Observatorio de Jicamarca - IGP RESUMEN El estudio de la Capa Ionosférica como ciencia, se inicia alrededor de los años 1924-1925, desarrollándose un conjunto de investigaciones con la finalidad de describirla físicamente. Su gran utilidad en las radiocomunicaciones, y el conocimiento obtenido sobre su comportamiento, buscan alcanzar un objetivo fundamental, su máximo aprovechamiento en las radiocomunicaciones. Un claro ejemplo del esfuerzo puesto por el hombre para cuantificar las características de la capa ionosférica, son los estudios realizados por Appleton y Barnet(1925), la teoría desarrollada por Chapman(1931) [Kohl et al., 1996], y el Modelo de Ionosfera Internacional de Referencia (IRI). Dentro de este marco de investigación, nuestro trabajo, estudiará las propiedades dinámicas de la ionosfera, centrándonos en la determinación del Campo de Velocidad de la Capa Ionosférica, para lograr este objetivo utilizaremos el algoritmo desarrollado para el seguimiento de partículas o estudios de flujos, denominado PIV (Particle Image Velocimetry) o DPIV (Digital Particle Image Velocimetry). INTRODUCCIÓN El desarrollo de técnicas para el estudio del movimiento de fluidos, se inicia con la técnica denominada PIV, el cual con la aparición y avance de la tecnología digital empezó a denominarse DPIV. Esta técnica, permite determinar en tiempo real el campo de velocidad de un fluido. Su aplicación hoy en día tiene distintos fines, en aerodinámica, se utiliza para conocer como actúan la presión del fluido sobre las diferentes partes de un objeto, con movimiento relativo al fluido; de manera similar, en la industria automotriz, permite conocer aquellas regiones donde el aire actúa generyo mayor fricción, otro estudio importante se da en las fluidos que se desplazan por tubos, aquí la técnica nos permite conocer aquellas regiones que serán sometidas a mayor o menor presión. Desde un punto de vista genérico, esta técnica se divide en etapas bien marcadas: Adquisición de las imágenes (digitalización de la información), procesamiento de las imágenes, determinación del desplazamiento y generación del campo de velocidades. Haciendo una analogía con nuestro trabajo de investigación, la diferencia fundamental y única, radica en la adquisición de las imágenes, el método DPIV, usa cámaras CCD para adquirir la imagen, mientras que las imágenes de la capa ionosférica (Figura 1) fueron adquiridas mediante el uso del radar, en el Radio Observatorio de Jicamarca. El procesamiento de la imagen, se da para mejorar algunas características de la imagen, ya sea para remarcar rasgos F. Galindo 78 o eliminar aquellas regiones generadas por ruido. La parte mas importante del desarrollo de este trabajo, esta centrada en la determinación del desplazamiento realizado por la Capa Ionosférica, es así que la utilización del método denominado Correlación Cruzada (Cross Correlation) nos permitirá determinar de manera rápida el desplazamiento, para finalmente generar el Campo de Velocidad en la región de la Capa Ionosférica estudiada. Figura 1. Imagen digital de una región de la Capa Ionosférica estudiada, adquirida mediante radar el 27 de Noviembre del 2003. FUNDAMENTOS El cálculo del campo de velocidades se basa en la definición de la magnitud física denominada velocidad media: t r v ∆ ∆ = donde, ν representa la velocidad media, ∆r el desplazamiento del objeto entre las imágenes, y ∆t el intervalo de tiempo entre las imágenes, este último parámetro es un valor constante y definido por el proceso de adquisición de la imagen. El desplazamiento se determina mediante el uso de la correlación, es así que la confiabilidad de nuestros resultados estarán determinados por el valor obtenido del desplazamiento. Correlación Cruzada La correlación cruzada es una herramienta que me permite determinar el grado de similitud entre dos imágenes, obteniéndose gran correlación si dichas imágenes corresponden a objetos con gran similitud, así los objetos de las imágenes se encuentren desplazados dentro de las imágenes, caso contrario la correlación es baja si los objetos no son iguales o presentan marcada diferencia. La correlación de una imagen con ella misma se denomina auto-correlación. Utilizando este concepto, de esta manera podemos discriminar y determinar el movimiento del objeto, ya que la correlación no solo me determina con certeza la similitud sino además me permite calcular el desplazamiento de los objetos entre las imágenes. Matemáticamente la correlación cruzada [Oppenheim y Schafer, 1975] se define en su forma discreta como: Análisis Dinámico en Imágenes de la Capa Ionosférica 79 ∑ − = += 1 0 )()()( N tgftC τ ττ donde, C(t) representa la correlación de las funciones f y g, N es el número datos a trabajar y t es el llamado “lag” o retraso. La forma continua de la correlación cruzada es: τττ dtgftC ∫ ∞ ∞− += )()()( Una breve representación del proceso es ilustrado en las Figuras 2, 3 y 4. Figura 2. Las gráficas muestran objetos que se desplazan con el transcurrir del tiempo. | Figura 3. Las gráficas mostradas anteriormente, pero superpuestas. Notamos más claramente el pequeño desplazamiento realizado. Figura 4. Función correlación obtenida para esta secuencia de imágenes. La ubicación del pico en el plano, me permite afirmar que dicha imagen se desplazo 1 píxel a la derecha horizontalmente y 4 píxeles hacia arriba. Una forma practica y útil de utilizar el algoritmo de la correlación cruzada es implementarla haciendo uso de la transformada rápida de Fourier (Fast Fourier Transform, FFT) [Brookner, 1977], la cual genera un ahorro importante, en la carga de procesamiento, reduciendo esta, de N2 a N*log2N operaciones, donde N representa el número de datos que contiene la imagen (o la señal) a analizar y el logaritmo es de base 2. El algoritmo utilizado para generar le correlación es: )}}({)}({{ *1)( tgFtfFFtC −= donde, F{ f } simboliza la transforma de Fourier de la función f, * representa la conjugada compleja y F-1{ } la transformada inversa de Fourier. El pico de la correlación (Figura 5) me informa que existe una gran similitud entre los datos (imágenes) que estamos trabajyo, de ahí podemos extraer, que el objeto analizado es el mismo. La ubicación del F. Galindo 80 pico dentro de la gráfica (Figura 6) determinará el desplazamiento realizado por el objeto entre imágenes, para esto se establece un sistema de referencia que me permite relacionar la posición del pico con el desplazamiento. Figura 5 .Gráfica de la función correlación, notemos la diferencia entre el pico máximo y los demás pico. Esto permite definir la similitud entre los objetos de las imágenes, así también el desplazamiento del objeto. Figura 6. Gráfica de contorno de la función correlación, el máximo determinado en la parte inferior derecha, (Punto color negro) determina la dirección y valor del desplazamiento. Dentro de las condiciones necesarias para el uso de esta técnica, debemos mencionar la no deformación del objeto como tal, entre las imágenes o una alteración mínima, así como también debe existir una constancia de la densidad, y el movimiento del objeto debe ser traslación pura, por ello que estos factores conllevan a establecer preguntas importantes, tal como cual debe ser el tiempo optimo entre imágenes, o el tamaño necesario de la ventana de análisis, a fin de obtener buenos resultados. Esquema del proceso para determinar el Campo de Velocidad en imágenes de la Capa Ionosférica La determinación del campo de velocidades, en las imágenes ionosféricas se realizara mediante las siguientes etapas. Primeramente la imagen se descompondrá en sus colores básicos (Figura 7) y se procederá a trabajar con componente. Figura 7. Descomposición de la imagen en sus componentes fundamentales (RGB), el gráfico muestra una componente, la intensidad en la escala de grises es un indicador del peso que tiene la componente sobre la imagen. Análisis Dinámico en Imágenes de la Capa Ionosférica 81 Como segunda etapa la imagen será procesada con el fin de resaltar bien su forma (Figura 8). Figura 8. La gráfica muestra el proceso de realce realizado a una imagen, a fin de mejorar su notoriedad. En el tercer paso cada imagen obtenida será dividida, en pequeñas ventanas cuadradas (Figura 9), denominadas ventanas de interrogación, cuyas dimensiones deberán ser potencia de dos, debido al uso de la FFT en los algoritmos. Figura 9. División de la imagen, en las denominadas ventanas de interrogación. La cuarta etapa es utilizar el algoritmo de la correlación cruzada en todas las ventanas, a fin de determinar el desplazamiento. Una vez determinado el desplazamiento en cada componente de una misma ventana de interrogación, se determina la velocidad de esa ventana, como promedio ponderado sus valores y la intensidad total en cada componente. El proceso se aplica a toda la imagen generándose el campo de velocidades (Figura 10). Figura 10. Campo de Velocidades generado por el algoritmo descrito en este trabajo. RESULTADOS Y DISCUSIÓN En líneas generales la estimación del movimiento del plasma ionosférico, trajo consigo resultados alentadores, con miras a perfeccionar la técnica, en oposición a esto se hallo cierta dificultad a la hora de trabajar con la información a ionosférica, debido a las restricciones impuestas por la técnica. A continuación mencionamos los resultados obtenidos en el trabajo realizado. El sistema de referencia que me permite calcular el valor del desplazamiento, lo puedo establecer en el centro del gráfico F. Galindo 82 de contorno de la correlación, este sistema me permite relacionar la ubicación del máximo y el desplazamiento realizado, es importante mencionar que dicha relación fue derivada a partir de datos que simulaban el movimiento de un conjunto de puntos, así conociendo a priori el desplazamiento se llego a concluir la siguiente relación: MydyMxdx MyLMxLSi =∧=→ >∧> 2/2/: MydyLMxdx MyLMxLSi =∧−=→ >∧< 2/2/: LMydyLdx MyLMxLSi −=∧=→ <∧> 2/2/: LMydyLMxdx MyLMxLSi −=∧−=→ <∧< 2/2/: donde: Mx = Coordenada x del pico máximo My = coordenada y del pico máximo L = Ancho de la ventana dx = Desplazamiento del objeto en x dy = Desplazamiento del objeto en x De estas relaciones, conociendo la ubicación del máximo, determinamos el valor del desplazamiento. La discretización de la información nos genera un problema, en la determinación de la posición del pico de máxima correlación, para eliminar este problema se pretende realizar un ajuste alrededor del máximo de correlación, a fin de obtener resultados más precisos. El tamaño de la ventana de interrogación esta relacionada con la densidad de información contenida en la ventana y con el módulo del desplazamiento realizado en la imagen (Figura 11). Figura 11. Efecto de la densidad sobre el valor del desplazamiento. El tamaño de la ventana de interrogación para los diferentes gráficos observados, debido al desplazamiento de la imagen de la Capa Ionosférica sugiere trabajar con un rango de valores de 16 a 32 píxeles. Un problema observado, se presenta cuyo la región de la Capa Ionosférica tiene una dimensión mayor que la ventana de interrogación (Figura 12), debido a esto el algoritmo calcula un desplazamiento cero, ya que el objeto seguirá conteniendo toda la ventana de interrogación, ya que las regiones dentro de esta capa parecerán no sufrir cambios como para que el algoritmo identifique eso. Análisis Dinámico en Imágenes de la Capa Ionosférica 83 Figura 12. Notemos la marcada diferencia entre la dimensión de la ventana y la Capa Ionosférica, al desplazarse unos cuantos píxeles el algoritmo determinara desplazamiento cero. Otro factor importante que debemos mencionar es la falta de continuidad de una imagen con su consecutiva. En el análisis notamos que ciertas regiones de la imagen desaparecen súbitamente, como producto del proceso de adquisición que se desarrolla temporalmente, lo cual equivaldría a una falta de constancia de la densidad o distribución de la información en la imagen. A continuación se muestra algunos Campos de Velocidad obtenido en las imágenes trabajadas, debemos resaltar que aquellas regiones que contienen un punto representan un valor nulo de la velocidad y que las flechas direccionales son proporcionales al desplazamiento determinado en píxeles (Figuras 13 al 18). Figura 13. Campo de velocidad para una misma región de la imagen, analizada con distintos parámetros. Figura 14. Campo de velocidad generada con una ventana de interrogación de 16 píxeles. F. Galindo 84 Figura 15. Campo de velocidad generada con una ventana de interrogación de 16 píxeles. Figura 17. Campo de velocidad con una ventana de interrogación de 32 píxeles, comparémosla con la Figura 14. Figura 16. Campo de velocidad generada con una ventana de interrogación de 16 píxeles. Figura 18. Campo de velocidad con una ventana de interrogación de 32 píxeles, comparémosla con la Figura 15. Los programas desarrollados para los fines ya mencionados fueron generados mediante el software para programación científica IDL 6.0, de la empresa Research Análisis Dinámico en Imágenes de la Capa Ionosférica 85 Systems Inc. (RSI) que por su fácil manejo permitío desarrollar una librería para el presente trabajo. CONCLUSIONES El denominado método PIV o DPIV es una técnica que a evolucionado mucho con el transcurrir de los años y los resultados que con él se pueden obtener son bastantes óptimos, usyo una lógica similar; esperaríamos que el desarrollo de este trabajo, el cual utiliza los mismos criterios establecidos por PIV o DPIV, genere resultados adecuados, mas existe una marcada diferencia, mas allá de ser simplemente fluidos; que nos hace reflexionar y estimula a seguir estudiyo de manera mas profunda y detallada los parámetros que me permitan caracterizar la velocidad de la ionosfera. Es así que este trabajo sirve como inicio, para una investigación más exhaustiva del problema. Como mencionamos líneas mas arriba los resultados son alentadores, en muchos casos notamos que desplazamiento determinado coincidía con el movimiento observado en las imágenes ionosféricas, mas la diferencia entre los criterios aplicados para el desarrollo del PIV o DPIV y su uso en imágenes ionosféricas trae consigo errores, que se manifiestan en los resultados del Campo de Velocidad. Se espera desarrollar criterios que permitan perfeccionar o refinar aquellos parámetros que afectan de manera directa el cálculo de la velocidad. Es por eso que este trabajo todavía esta en una etapa de desarrollo y se espera en un futuro lograr resultados óptimos. AGRADECIMIENTOS Un agradecimiento al Dr. Jorge L. Chau por brindarme la oportunidad, confianza, libertad y motivación para el desarrollo de este trabajo, a su vez por todo su apoyo y conocimientos brindados a mi persona. A todo el personal del Radio Observatorio de Jicamarca por su amabilidad y apoyo. BIBLIOGRAFÍA Brookner E. (1977): Radar Technology, 1st Edition, Artech House, pp. 149-154 Kohl, H. Ruster, R. y Schlegel K.(1996): Modern Ionospheric Sciencie, 1st Edition, Europian Geophysical Society, pp. 7-8. Oppenhiem A. y Schafer R. (1975): Digital Signal Processing, 1st Edition, Prentice-Hall., pp. 554-562. F. Galindo 86 Compendio de Trabajos de Investigación CNDG – Biblioteca Instituto Geofísico del Perú. V. 5 (2004) p. 87 - 100 MEDICIÓN DE PARÁMETROS SCATTERING DE UN TRANSISTOR PARA EL DISEÑO DE UN AMPLIFICADOR DE BAJO RUIDO WILBERT JESUS VILLENA GONZALES Especialidad de Ingeniería Electrónica Facultad de Ciencias e Ingeniería Pontificia Universidad Católica del Perú wvillena@jro.igp.gob.pe Prácticas dirigidas por: Ing. Fernando Villanueva Ruiz Ing. Otto Castillo Gonzáles Radio Observatorio de Jicamarca - IGP RESUMEN Un amplificador de bajo ruido (LNA) juega un papel de suma importancia en el diseño de cualquier tipo de receptor, ya que su función principal es la de amplificar señales extremadamente pequeñas tratando de añadir la menor cantidad de ruido posible, esto es, preservando el nivel de relación señal a ruido (SNR) del sistema, por lo que el diseño de un buen LNA dependerá de cuan bien se conoce al dispositivo central o corazón del amplificador y que mejor manera sino definiendo sus parámetros S los cuales describen de manera completa el comportamiento del sistema a la frecuencia de trabajo requerida. Por tanto, el presente trabajo trata de explicar de una manera sencilla la teoría detrás de los parámetros S y la forma mas adecuada para su medición usando elementos relativamente fáciles de encontrar en cualquier centro dedicado a las telecomunicaciones o ramas afines. INTRODUCCIÓN Uno de los puntos mas difíciles en el diseño de un LNA es que debido a la característica no-lineal del dispositivo éste puede exhibir comportamientos muy distintos en cada rango de frecuencias, especialmente en frecuencias altas. Esto debido principalmente a la presencia de capacitancias parásitas propias del amplificador las cuales van influyendo cada vez mas y mas a medida que se incrementa la frecuencia de trabajo. Es por eso necesario tener una caracterización del elemento (transistor, tubo, nuvistor, etc.) en la frecuencia a la que se utilizará, ya que de esa forma se obtendrán los mejores resultados en el diseño. La forma de caracterizar un dispositivo es mediante una red de dos puertos, es decir, mediante un conjunto de relaciones de entrada-salida, lo cual tiene muchas ventajas entre las cuales está su uso práctico sin la necesidad de conocer la estructura interna del sistema por lo que esta metodología a modo de “caja negra” ha tenido un enorme apego por parte de muchos ingenieros cuya mayor preocupación está en el funcionamiento global del sistema y no en el análisis de cada uno de sus componentes individuales por lo que este enfoque es especialmente importante en los circuitos de RF y microondas, donde la solución completa de las ecuaciones de Maxwell son muy difíciles de obtener o dan mas información de la que generalmente se necesita en diseños prácticos de por ejemplo filtros, resonadores ó amplificadores [Ludwig y Bretchko., 2000]. W. Villena 88 Para definir una red dos puertos existen varios parámetros que se pueden usar de los cuales los mas reconocidos para el diseño de circuitos son los parámetros Scattering o “S”, no solo por el tipo de información que proveen sino, por la facilidad con la que pueden obtenerse a comparación de otros parámetros cuyo proceso de medición implica que el dispositivo debe estar perfectamente abierto o cortocircuitado lo que puede ser demasiado difícil de lograr, especialmente a altas frecuencias donde las inductancias y capacitancias propias del dispositivo hacen muy difícil de obtener el valor correcto. Por otro lado, la ventaja de los parámetros S radica en que al implicar ondas viajeras éstas no varían en magnitud a lo largo de las líneas de transmisión, a diferencia de voltajes y/o corrientes, por lo que los parámetros S de un dispositivo pueden ser medidos estando éste a una determinada distancia de los equipos de medición. FUNDAMENTOS Los parámetros “S” son simplemente descriptores de potencia de una onda que nos permiten definir relaciones de entrada- salida de una red en términos de ondas viajeras incidente y reflejada. De acuerdo a la Figura 1 se puede ver que una onda viajera es aquella generada por una determinada fuente de voltaje (V1+) y que viaja al puerto 1 de una red a través de una línea de transmisión cuya impedancia característica es Zo. Cuando la onda alcanza la red dos nuevas ondas viajeras se generan. Una aparece en el puerto 2 alejándose de la red (V2+) y la otra aparece en el puerto 1 viajando de regreso a la fuente (V1-), por lo que los parámetros “S” caracterizan la red indicando la cantidad de potencia reflejada en ambos puertos ( oZV /1− , oZV /2− ) en relación a la cantidad de potencia incidente en cada uno de ellos ( oZV /1+ , oZV /2+ ) [Swartz, 1995]. Figura 1. Esquema de la transmisión y-división de una onda viajera a su paso por una red de dos puertos. Como los parámetros “S” están basados en las características de reflexión y por tanto, en las relaciones de potencia en una red, uno de sus enfoques más útiles es la representación de una red en términos de potencias. Así, se define un nuevo juego de parámetros (a1, b1) y (a2, b2) donde a y b son potencias normalizadas al valor de impedancia caraterística Zo de la línea de conexión usada. )( 2 1 )( 2 1 non o n non o n IZV Z b IZV Z a −= += Entonces, oZ V a + = 1 1 , oZ Vb − = 1 1 Medición de parámetros scattering de un transistor para el diseño de un amplificador de bajo ruido 89 oZ V a + = 2 2 , oZ Vb − = 2 2 Por tanto, y concluyendo lo anteriormente dicho, cualquier onda viajera está constituida por dos componentes: incidente y reflejada. −+ += nnn VVV Figura 2. Convención usada para definir los parámetros S de una red de dos puertos. Por ejemplo, la potencia total que fluye a una carga consiste de la porción de (a2) que es reflejada a la salida de la red y de la porción de (a1) que es transmitida a través de la red. Asimismo, la potencia total que fluye de la entrada de la red de dos puertos hacia la fuente consiste de la porción de (a1) que es reflejada a la entrada de la red y de la fracción de (a2) que es transmitida a través de la red. Por lo que para la red de la Figura 2 se obtienen las siguientes relaciones: 2221212 2121111 aSaSb aSaSb += += Significado de los parámetros S En base a la convención de direcciones de la Figura 2 y a las relaciones mostradas anteriormente, la definición de cada uno de los parámetros S es: 1 1 021 1 11 puertoelenincidentepotencia puertoelenreflejadapotencia a bS a == = 1 2 021 2 21 puertoelenincidentepotencia puertoalatransmitidpotencia a bS a == = 2 2 012 2 22 puertoelenincidentepotencia puertoelenreflejadapotencia a bS a == = 2 1 012 1 12 puertoelenincidentepotencia puertoalatransmitidpotencia a bS a == = Esto quiere decir que S11 es una medida del coeficiente reflexión a la entrada de la red, es decir, indica la porción de la onda proveniente de la fuente que retorna a ella por lo que se debe garantizar que a2=0 y que la carga sea igual a la impedancia característica Zo en caso se quiera hallar su valor. S21, llamado coeficiente de transmisión directa, indica el grado de amplificación de la red para un determinado nivel de onda a su paso por ella, por lo que también se debe garantizar a2=0 y Zcarga=Zo para hallar su valor. S22 es una medida del coeficiente de reflexión a la salida de la red, es decir, indica la porción de la onda proveniente de la carga que retorna a ella, por lo que se debe garantizar que a1=0 y que la impedancia de la fuente sea igual a Zo para poder calcularlo. Finalmente, S12, llamado coeficiente de transmisión inversa, indica el nivel de amplificación de la red para una onda que ingresa por su puerto de salida (puerto 2). En el caso de transistores utilizados como amplificadores, cuanto mas pequeño sea el valor de S12, mejor será la estabilidad y rendimiento del W. Villena 90 amplificador, ya que lo que hace este valor es reducir la ganancia total del dispositivo, es decir, es una “ganancia negativa” (Anderson et al. 1997). Cálculo de parámetros S Para poder calcular los parámetros S es necesario asegurar una condición de adaptación perfecta de impedancias, esto es, asegurar que las impedancias de la carga y fuente sean iguales a la impedancia característica de la línea de transmisión de tal modo que no exista onda reflejada en el lado opuesto del grupo de parámetros a medir ya que de esa forma se asegura que el cálculo del parámetro requerido sea independiente de lo que hay al otro lado de la red. Por ejemplo, tal como se dijo anteriormente, para hallar S11 y S21 se necesita que a2=0, esto quiere decir que la onda que sale por el puerto 2 (b2) debe ser totalmente absorbida por la carga. Igualmente, para el caso de S22 y S12, se debe asegurar que a1=0, el cual viene a ser la onda reflejada por la carga de la fuente colocada a la entrada de la red. Como se estableció anteriormente, S11 y S22 representan los coeficientes de reflexión ( Γ ) a la entrada y salida de la red respectivamente, por lo que hay dos formas de poderlos calcular. La primera es utilizando directamente la relación para el cálculo de coeficientes de reflexión, la cual relaciona las impedancias de entrada y salida de la red (Zin, Zout ) con Zo pero esto resulta muy difícil cuando no se sabe exactamente lo que hay dentro de la red o cuando lo que está dentro es un circuito o sistema demasiado difíciles de analizar (Figura 3). oin oin in ZZ ZZ + − =Γ , oout oout out ZZ ZZ + − =Γ Por eso, en este caso se recurre directamente a las relaciones de potencia de cada parámetro que se establecieron anteriormente. 012 2 012 2 22 021 1 021 1 11 == == == == Eii r a Eii r a E E a bS E E a bS Por lo que la dificultad de hallar S11 y S22 se reduce únicamente a medir los voltajes reflejado e incidente en los puertos 1 y 2 de la red respectivamente. Para el cálculo de S21 y S12 también se usan las relaciones mostradas en la sección anterior pero para este caso se tiene que hacer un pequeño desarrollo matemático. 1 2 0 11 22 021 2 21 2 )( 2 1 )( 2 1 22 G VI o o o o a V V IZV Z IZV Z a bS = + − == == − = ++ 2 1 0 22 11 012 1 12 2 )( 2 1 )( 2 1 11 G VI o o o o a V V IZV Z IZV Z a bS = + − == == − = −− Medición de parámetros scattering de un transistor para el diseño de un amplificador de bajo ruido 91 Por lo tanto, el cálculo de S21 y S12 se reduce nuevamente a mediciones de voltaje en entrada y salida Figura 3. (a) Medición de S11 y S21 adaptando la impedancia de línea Zo en el puerto 2 a su correspondiente impedancia de carga ZL=Zo. (b) Medición de S22 y S12 adaptando la impedancia de línea Zo en el puerto 1 a su correspondiente impedancia de fuente ZG=Zo. Medición de los parámetros S La forma mas básica y fácil de medir los parámetros “S” es mediante tres elementos: un generador de señales, un acoplador direccional y un voltímetro vectorial (VVM). El generador de señales es usado para proporcionar la señal con la frecuencia necesaria para las mediciones. El acoplador direccional es aquel instrumento necesario para poder aislar las ondas incidente y reflejada y poder así medirlas por separado mediante el voltímetro vectorial. Tiene esencialmente tres puntos de conexión: entrada, salida y un punto de muestreo que proporciona un nivel de señal equivalente (pero no igual) a la señal de entrada, es decir, con la misma frecuencia pero con un nivel de atenuación determinado. El Voltímetro Vectorial es un equipo de medición que ante dos señales de entrada, éste otorga principalmente dos datos importantes: Desfase entre dichas señales y relación de magnitud entre ellas. Medición S11 y S22 La conexión básica es la mostrada en la Figura 4 en la que la lectura del voltaje en el canal A del VVM (AD) es proporcional a la amplitud de la onda entrante al dispositivo (a1D). Similarmente, el voltaje en el canal B (BD) es proporcional a la amplitud del voltaje reflejado del dispositivo (b1D), por lo que se puede escribir: W. Villena 92 D BD D AD bKB aKA 1 1 = = donde KA y KB son constantes que dependen de los cables de conexión. Ya que a2D es cero por la presencia de la carga ZL=Zo en el puerto 2, S11 está dado por: A D B D D D K A K B a bS == 1 1 11 Para encontrar los valores de KA y KB es necesario hacer una segunda medición con un dispositivo de prueba (DUT- Dispositive Under Test) conocido. A esto se llama una medida de calibración. Si el DUT es removido y reemplazado por un cortocircuito lo que teóricamente debe resultar es que el voltaje en el canal A y en el canal B sean iguales y que estén desfasados 180º pero debido a la presencia del comportamiento no-ideal de los cables y de los mismos conectores los valores obtenidos son: S BS S AS bKB aKA 1 1 = = donde a1S es la amplitud del voltaje entrante al cortocircuito y b1S es la amplitud del voltaje reflejado del mismo en el puerto 1. Sin embargo, ya que están desfasados 180º el cociente de estas amplitudes es -1. 1 1 1 −== A S B S S S K A K B a b por tanto: S S A B A B K K −=              −= S S D D A B A B S11 El mismo procedimiento se emplea para hallar S22 con la salvedad que ahora la señal de la fuente ingresa por el puerto 2. Figura 4. Interconexión de equipos para la medición de S11 y S22 Medición de parámetros scattering de un transistor para el diseño de un amplificador de bajo ruido 93 Medición S21 y S12 El procedimiento necesario para hallar estos dos parámetros es prácticamente el mismo que el descrito para el caso de S11 y S22, sólo que ahora se deben muestrear los puertos de entrada y salida al mismo tiempo ya que se debe ver cuan grande es la salida respecto del nivel de entrada, por lo que el esquema de interconexión es el mostrado en la Figura 5 en el que se puede ver que el DUT es conectado directamente entre dos acopladores direccionales, el voltaje en el canal A es proporcional al voltaje que entra al dispositivo bajo medición (DUT) y el voltaje leído en el canal B es proporcional al voltaje transmitido a través del DUT. Para S21 se obtiene: D B D D A D bKB aKA 2 1 = =  A D B D D D K A K B a bS == 1 2 21 Para encontrar ahora las constantes se debe hacer otra medida de calibración. En lugar del DUT se coloca una línea que conecte los dos acopladores. En este caso, si los dos acopladores tienen el mismo factor de acoplamiento, las magnitudes de los voltajes en los canales A y B deben ser las mismas. Por lo tanto, E B E E A E bKB aKA 2 1 = = donde 1 1 1 == A E B E E E K A K B a b E E A B A B K K −=∴              −= E E D D A B A B S21 Para el caso de S12 el procedimiento es exactamente el mismo que S21, sólo que ahora la señal entra por el puerto 2 y sale por el puerto 1. Figura 5. Interconexión de equipos para la medición de S21 y S1 W. Villena 94 Mediciones de calibración Para el caso de S11 y S22, el DUT debe ser reemplazado por un cortocircuito o por un circuito abierto. En el primero de los casos el VVM debe dar por resultado o1801∠ en un caso ideal ya que en este caso la señal al llegar al cortocircuito retornará con su mismo nivel de amplitud pero desfasada 180º. En el segundo caso, la amplitud continúa siendo la misma así como la fase, por lo que el resultado del VVM es o01∠ en el caso ideal. Para el caso de S21 y S12 el DUT debe ser reemplazado por un cable, de tal modo que el nivel de señal y la fase en los canales A y B serán los mismos, por lo que el resultado dado por el VVM debe ser o01∠ en un caso ideal. Polarización de dispositivo En el caso de transistores y tubos es necesario seleccionar el punto de reposo en el que trabajará el dispositivo y diseñar el circuito para lograr dicho punto ya que los parámetros dependen directamente de los niveles de voltaje y corriente presentes. El grado de dependencia dependerá de las características de cada dispositivo. RESULTADOS Y DISCUSIÓN En líneas generales, la medición de parámetros S viene a ser un método completamente rutinario ya que simplemente se interconectan los equipos y el VVM muestra el valor requerido. Como forma de mostrar los métodos de medición descritos anteriormente se medirán los parámetros S de un transistor, el MGF1302, el cual es un transistor GaAsFET de bajo ruido cuyos parámetros sólo están dados por el fabricante de 500MHZ a 12GHz por lo que se harán mediciones a frecuencias menores, comprendidas entre 50MHz y 250MHz para ver el comportamiento del dispositivo a frecuencias menores. Transistor MGF1302 Lo primero que se hace es seleccionar el punto de polarización en el que trabajará el dispositivo y diseñar el circuito para ello. En el caso específico del MGF1302 el punto escogido es: VVDE 3= , mAI D 10= ya que en este punto están los parámetros dados por el fabricante, siendo el circuito obtenido el mostrado en la Figura 6 en donde se puede observar la presencia de bobinas de choke (RFC) de tal modo que éstas actúen como impedancias altas en AC y así no permitan que la fuente de DC interfiera en los resultados. Asimismo, se usan capacitores de acoplo y desacoplo (CB) cuya función es actuar como cortocircuitos en AC y permitir únicamente el flujo de la señal a lo largo del dispositivo mas no así del nivel DC. Medición de parámetros scattering de un transistor para el diseño de un amplificador de bajo ruido 95 Los resultados de las mediciones de parámetros S son las mostradas en la Tabla 1 así como en la Figura 7 en donde además están los datos dados por el fabricante hasta 12GHz para así poder ver la tendencia en el comportamiento del dispositivo conforme aumenta la frecuencia. Figura 7. Circuito de polarización del FET (Vds=3V, Ids=10mA, Rs=75Ω, Rd=100Ω, RG=10kΩ, CB=0.1uF, VDD=5V) Con respecto a S11 se puede ver que conforme aumenta la frecuencia su valor decrece lo cual quiere decir que a mayores frecuencias el valor en magnitud de la impedancia de entrada del dispositivo se va acercando al valor de la impedancia de la fuente (para este caso es de 50Ω. Con respecto a S22 se puede ver que sucede lo contrario de S11, ya que la impedancia de salida va disminuyendo conforme disminuye la frecuencia. Tabla 1. Resultados de las mediciones de parámetros S del FET MGF1302 para el rango de frecuencias entre 50MHz y 500MHz expresados en magnitud (absoluta) y fase (en grados Sexagesimales). Las mediciones se hicieron de 50MHz a 250MHz, mientras que el valor a 500MHz es el valor dado por el fabricante, ya que se empalmarán los resultados medidos con los dados por el fabricante. Frec. (MHz) MAG ANG MAG ANG 500 0.997 -13.3 0.664 -10.3 250 0.991 -11.8 0.64 -9.9 200 0.987 -11.5 0.64 -9.7 150 0.976 -8 0.61 -8.3 100 0.975 -6 0.575 -8 50 0.97 -5.3 0.57 -3 S11 S22 Frec. (MHz) MAG ANG MAG ANG 500 3.809 167.6 0.019 80.1 250 3.83 174 0.003 157 200 3.93 179 0.003 170 150 3.92 -177 0.008 -123 100 4.01 -173 0.01 -118 50 5.88 -135 0.013 13 S21 S12 Con respecto a S21 se puede ver que la ganancia del dispositivo aumenta conforme disminuye la frecuencia por lo que conviene usar a este transistor como un amplificador de alta ganancia a frecuencias relativamente bajas en comparación con frecuencias del orden de los Gigahertz. Con respecto a S12, se puede ver que conforme aumenta la frecuencia aumenta la magnitud de S12, esto quiere decir que aumenta la realimentación entre la entrada y la salida pudiéndose llegar incluso a la oscilación en caso no se tomen las medidas pertinentes al momento de diseñar el amplificador. W. Villena 96 S11 S22 S21 S12 Figura 7. Parámetros “S” del FET MGF1302 para el rango de 50MHz a 12GHz Asimismo, resulta interesante poder comparar los valores medidos con los obtenidos teóricamente mediante extrapolación lineal de los valores dados por el fabricante, tal como se observa en la Figura 8 en donde se puede comprobar la no-linealidad del transistor. Por lo que si en un diseño específico se usan valores teóricos puede llegarse a resultados erróneos tal como es el caso de S11 cuyo valor para bajas frecuencias es mayor a uno lo que significaría que la onda reflejada a su entrada es mayor que la incidente lo que en pocas palabras quiere decir que el componente es altamente inestable y con muchas probabilidades de oscilar a esas frecuencias lo cual en realidad no sería así ya que como se ve, su verdadero valor (en magnitud) es menor a la unidad. Asimismo, para el caso de S22, Medición de parámetros scattering de un transistor para el diseño de un amplificador de bajo ruido 97 aunque pareciera que hay mucha diferencia entre los valores medido y calculado (especialmente en la zona de bajas frecuencias), la escala de la gráfica nos dice que no es así, ya que por ejemplo para la frecuencia de 50MHz la diferencia entre valores es de 0.1. Esto no hace nada mas que ratificar la veracidad de las mediciones hechas y al mismo tiempo probar la hipótesis inicial acerca de la No- linealidad del dispositivo en distintos rangos de frecuencia. Figura 8. Comparación de los parámetros S medidos con los obtenidos mediante extrapolación a partir de los valores dados por el fabricante a partir de los 500MHz. W. Villena 98 Por otro lado, la diferencia entre los valores medido y calculado para el caso de S12 se va incrementando a medida que aumenta la frecuencia (el valor medido tiende a decrecer mientras el calculado tiende a aumentar) lo que quiere decir que el transistor ha sido diseñado para un rango de frecuencias intermedio entre el rango de RF y de las microondas, ya que como se dijo anteriormente, un transistor ideal es aquel en el que su ganancia inversa (S12) es cero. CONCLUSIONES En general, la medición de parámetros S es un procedimiento rutinario y simple una vez que se tiene un entendimiento básico de lo que son las ondas reflejada e incidente por lo que lo tal vez lo único difícil sea tener todos los elementos necesarios para la medición. Asimismo, habría que tener dos puntos importantes respecto a las mediciones y son que la medición del parámetro S12 resulta algo difícil de realizar ya que los niveles de señal son tan pequeños que el Voltímetro Vectorial no da un valor fijo o estable, es decir, hay mucha variación por lo que es necesario hacer varias mediciones de este parámetro y tomar (de ser necesario). Esto se puede solucionar teniendo a la mano un VVM de mayor sensibilidad o lo mas factible sería tener un acoplador direccional cuyo factor de acoplamiento sea mas pequeño para que de esa forma llegue mas señal al VVM y así obtener una medición mas exacta. El otro punto que hay que tener en cuenta es que conforme se hagan mediciones a frecuencias mayores, éstas pueden resultar mas inexactas ya que empiezan a influir efectos capacitivos e inductivos propios de los conectores y de las mismas pistas del circuito de polarización del transistor por lo que se debe procurar el menor número de conexiones posibles y que el tamaño de las pistas sea lo mas pequeño posible. Finalmente, habría que agregar que el método descrito en el presente trabajo es uno de varios métodos existentes entre los que destacan el uso de un Analizador de Redes (Network Analyzer) el cual puede es capaz de medir magnitud y fase de una red simple o dual, por lo que su ventaja radica en que cada una de las unidades funcionales asociadas al Voltímetro Vectorial están incorporadas en un único instrumento para una medición completamente automatizada del dispositivo de RF o MW. Otro método es mediante la prueba de tonos la cual usa un tono sinusoidal de prueba combinado con otras señales, obteniéndose los parámetros mediante un posterior procesamiento de señales. Pero este método es generalmente usado para dispositivos de potencia los cuales son afectados muchas veces por el ruido de la red o del medio ambiente. Medición de parámetros scattering de un transistor para el diseño de un amplificador de bajo ruido 99 AGRADECIMIENTOS Mi agradecimiento al Dr. Jorge L. Chau por otorgarme desde un inicio la oportunidad y confianza necesarias para realizar esta clase de proyectos. Asimismo, a Fernando Villanueva, Darwin Córdova y Otto Castillo por su paciencia y apoyo constantes en todas las dudas que se presentaban durante el proyecto. BIBLIOGRAFÍA Anderson, D., Smith, L. y Gruzynski, J. (1997). S-Parameter Techniques for Faster, more Accurate Network Design. Hewlett Packard -Application Note 95-1. Ludwig, R. y Bretchko, P. (2000). RF Circuit Design: Theory y Applications. Worcester Polytechnic Institute. Pags.168- 189 Swartz, T. (1995). The Development of an Automated Measurement System. EEAP 399. W. Villena 100 Compendio de Trabajos de Investigación CNDG – Biblioteca Instituto Geofísico del Perú. V. 5 (2004) p. 101 - 112 APLICACIÓN DE LA DISTRIBUCIÓN DE POISSON PARA EL CALCULO DEL PERIODO DE RETORNO DE LOS SISMOS JUAN CARLOS VILLEGAS Escuela Profesional de Ingeniería Geofísica Facultad de Geología, Geofísica y Minas Universidad Nacional de San Agustín - Arequipa jcvillegas@axil.igp.gob.pe Prácticas dirigidas por: Dr. Hernando Tavera H. Centro Nacional de Datos Geofísicos RESUMEN La sismicidad en el Perú es producto de la convergencia entre las placas de Nazca y la Sudamericana. A fin de conocer las características de la sismicidad en la Región Central del Perú (-9.5º, -14.5º y -73º, - 81º) y la posible ocurrencia de sismos de magnitud elevada, en el presente estudio, se realiza un análisis sismo-estadístico para estimar los siguientes parámetros: el periodo medio de retorno, probabilidad de ocurrencia, frecuencia sísmica e índice de sismicidad de grandes sismos que puedan afectar a esta región. Para tal objetivo se hace uso de la relación Gutemberg - Richter (Log N = a – bM), con la cual se calculan los valores de las constantes a y b. Los datos utilizados en el presente estudio fueron extraídos del catalogo del IGP (1980-2000). Los resultados indican la ocurrencia de importante actividad sísmica en esta región caracterizada con un valor de b de 1.54, típico de áreas orogénicas donde el potencial sísmico es muy alto. Además, en relación a la probabilidad de ocurrencia, sismos de mb = 5.0 ocurrirían en periodos de 10 años con una probabilidad de 98%, sismos con mb = 5.5 cada 20 años con una probabilidad de 80 % y eventos de mb = 6.0 cada 55 años con una probabilidad del 50%. INTRODUCCIÓN El Perú se encuentra ubicado en la parte central del borde occidental de Sudamérica y se caracteriza por ser una de las regiones sísmicas mas activas del mundo. Esta sismicidad se debe al proceso de subducción de la placa de Nazca bajo la placa Sudamericana con una velocidad de 8cm/año en dirección N80º (Minster y Jordan 1978), en el interior del continente la sismicidad, en mayor porcentaje, se debe a las deformaciones corticales o sistemas de fallas En la parte Central del Perú, la sismicidad se ha caracterizado por la ocurrencia continua de sismos, como los de 1586 (8.1), 1655 (7.4 mb), 1664 (7.8 mb), 1687 (8.2 mb), 1746 (8.6), 1806 (7.5 mb), 1940 (8.2 mb), 1947 (7.5 mb), 17 de octubre de 1966 (7.8 Ms), 31 de Mayo de 1970 (7.8 Ms) y 1974 (7.5 Ms). Todos estos sismos se originaron en el proceso de fricción de placas. En el presente estudio se realiza un análisis de las características espaciales de la sismicidad de la Región Central del Perú (Figura 1) a fin de estimar los siguientes parámetros: Frecuencia Sísmica, índice de Sismicidad, Periodo medio de Retorno y Probabilidad de Ocurrencia (riesgo sísmico) de sismos de J. Villegas 102 una determinada magnitud. El análisis de recurrencia se realiza utilizando la distribución de Poisson ya que dicho procedimiento se constituye como uno de los modelos mas sencillos para estudiar la recurrencia de los terremotos. La información sísmica utilizada fue extraída del Catálogo Sísmico del Instituto Geofísico del Perú, para el periodo comprendido entre los años comprendida entre 1980 y 2000 (Mb ≥ 3.5). Figura 1. Ubicación geográfica de la región de estudio. CARACTERÍSTICAS TECTONICAS Y SISMICIDAD En la Región Central del Perú, la mayor parte de sismos se generan por la energía liberada del contacto entre placas y su cantidad dependerá de la velocidad con que se desplacen y la fricción que exista entre ellas. Este proceso, conocido como subducción ha dado origen a la formación de la Cordillera de los Andes que se extiende de Norte a Sur a lo largo del continente Sudamericano, la cual tiene un ancho de 250 Km, en general esta se conforma de una ancha franja de cadenas montañosas con un vulcanismo importante que se distribuye en forma paralela al margen de subducción; por su actual topografía se distingue cinco unidades estructurales bien definidas, Figura 2 (Tavera y Buforn 1998). Estas unidades son: La Zona Costera (ZC).- Zona estrecha de aproximadamente 40 km de ancho que se extiende de norte a sur y esta constituida en su mayoría por suaves plegamientos volcánicos y rocas sedimentarias del Mesozoico. La zona sur está formada por basamentos de rocas cristalinas plegadas y sujetas a deformación desde el precámbrico. La Cordillera Occidental.- En el Perú se distribuye paralelo a la costa de norte a sur. La parte más elevada de esta Cordillera (4200-4500 m), esta formada por series del Mesozoico, plegadas y cubiertas de manera heterogénea por capas volcánicas del Cenozoico. Esta Cordillera aumenta su anchura en el Sur del Perú. El Altiplano.- Situada entre las Cordilleras Occidental y Oriental. En la región Sur su anchura es de 200 km. y se extiende hacia el norte hasta los 9 °S, en donde alcanza un ancho de unos 50 Km y después desaparece. Esta unidad esta formada por una serie de cuencas Aplicación de la distribución de poisson para el calculo del periodo de retorno de los sismos 103 intramotañosas del Cenozoico que se prolongan hacia el Altiplano. La Cordillera Oriental.- Menos elevada que la Cordillera Occidental (3700-4000 m.) y corresponde a un extenso anticlinal formado por depósitos intrusivos del Precambrico. En la región sur, la Cordillera se curva en dirección E-W para luego continuar paralela a las otras unidades. La Zona Subandina.- Zona de anchura variable donde se amortiguan las estructuras andinas. Esta zona se localiza entre la Cordillera Andina y la Llanura Amazónica y está formada por una cobertura de sedimentos del Mesozoico y Cenozoico. La Dorsal de Nazca, es una cordillera submarina que forma parte de la placa oceánica y esta ubicada a la altura de 15º S con una altura de 2000 a 4000 m orientado en dirección NE – SW perpendicular a la fosa Perú – Chile. La Fractura de Mendaña se cree que comprendería una porción de una antigua zona de divergencia de placas, está localizada a la altura de 10º S con una altura aproximada de 2000 m, se acepta que ambas estructuras están compuestas de basaltos. Figura 2. Principales unidades estructurales y estructuras Submarinas presentes en el borde occidental del Perú. (Tavera y Buforn 1998). El área claro indica alturas menores a 1500 m.; áreas grises, alturas entre 1500 y 4000 m; áreas oscuras, alturas mayores a 4000 m. En los últimos 60 años, la región Central del Perú fue afectada por 5 sismos de gran magnitud (Mw ≥ 7.0), los cuales presentaron procesos complejos de ruptura. De norte a sur se puede citar los siguientes. sismos: 24 de Mayo de 1940 (Mw = 7.8), 24 de Agosto de 1942 (Mw = 8.2), 17 de Octubre de 1966 (Mw = 8.0), 31 de Mayo de 1970 (Mw = 7.8), 3 de Octubre de 1974 (Mw = 8.0) y el de Nazca de 1996 (Mw = 7.7). En la Figura 3, se muestra la distribución espacial de los sismos ocurridos en la región Central del Perú entre los años 1980-2000. Los círculos rojos corresponden a sismos con foco superficial (h = 60 km) y los verdes a sismos con foco intermedio (61 < h ≤ 300). También se observa que la mayoría de sismos se J. Villegas 104 distribuyen entre la línea de costa y la fosa Perú - Chile, siendo estos de foco superficial e intermedio. En el interior del continente, se nota sismos de foco superficial e intermedio. Los primeros asociados al proceso de deformación cortical y a sistemas de fallas tales como Huaytapallana, Satipo, Amauta y Ayacucho. Los segundos están asociados al proceso de subducción a niveles de profundidad mayores y a la deformación interna de la placa de Nazca. Figura 3. Distribución espacial de sismicidad para la región Central del Perú Los círculos oscuros corresponden a sismos con foco superficial y los círculos claros a sismos con foco intermedio. RIESGO SISMICO En general, para estudios de riesgo sísmico se utilizan dos enfoques: El Enfoque Determinístico, supone que la sismicidad futura será igual a la pasada, siendo el máximo terremoto ocurrido el máximo previsible. El Enfoque Probabilístico, en el cual las distancias a las fuentes sísmicas potenciales y las magnitudes generadas por estas, se tratan como variables aleatorias. El resultado será una curva de peligro que representa la excedencia de un valor pre-especificado del movimiento en un lugar dado.(Escalante 2002). En el presente estudio se toma el enfoque probabilístico debido a la necesidad de evaluar la PO de sismos que puedan afectar la Región Central Perú. Relación Frecuencia Magnitud Esta relación forma parte del método probabilístico y fue propuesta por Gutenberg y Richter (1944) a fin de establecer una relación entre el número de sismos que ocurren en una región y sus respectivas magnitudes, ósea que, de Aplicación de la distribución de poisson para el calculo del periodo de retorno de los sismos 105 existir una disminución en la frecuencia de los sismos, la magnitud del sismo que pudiera ocurrir iría en aumento ya que se produce una mayor acumulación de energía sísmica, y de existir un aumento en la frecuencia de los sismos la magnitud del sismo a ocurrir disminuiría ya que la energía sísmica es liberada continuamente. Por ello la frecuencia de los sismos frente a su tamaño (N), tiende a tener una forma lineal fija (Figura 4): MeMF β−=1)( ; M > 0....(1) Figura 4. Relación lineal en función del logaritmo N y la Magnitud M de un sismo. )()( MbaMNLog −= ... (2) Donde N es el número de sismos que ocurre en cierto periodo de tiempo y M la magnitud de los sismos y puede ser dada en mb o Ms, los parámetros “a “y “b” son dos constantes que representan el número de sismos de magnitud mayor que cero, y la proporción de sismos con magnitudes pequeñas y grandes, respectivamente (Udias 1999). Los valores de las constantes dependen del periodo de datos considerado, del área,. de las propiedades físicas del medio y son indicativos del nivel de sismicidad. Índice de Sismicidad El índice de sismicidad define el número anual medio espectado de sismos con magnitud (M) mayor que una M determinada. Para ello es importante conocer la frecuencia acumulativa, para lo cual se integra la ecuación 2. )(')( MbaMNLog −= ... (3) donde 'a = a - Log (b Ln 10) a ! = 'a - Log (T) T = Periodo de Observación MbaN ⋅−⋅= 1010 ! ...(4) Periodo de Retorno En base a los valores de las constantes “a ” y “b” se puede calcular el periodo de retorno de futuros terremotos para un rango determinado de magnitud. El valor inverso de la ecuación 4 permite encontrar el periodo de retorno. N TR 1= ...(5) Probabilidad de Ocurrencia La probabilidad de ocurrencia de uno o más sismos de magnitud mayor que una J. Villegas 106 determinada, durante un periodo de tiempo (T) dado, se puede deducir de la distribución de probabilidad discreta que se ajusta a la forma exponencial e-ht con la sgte. ecuación: )exp(1)( TNMP ir ⋅−−= ...(6) Donde: N es el numero de sismos esperado por año y T el periodo de tiempo en años. Distribución de Poisson La distribución de Poisson es una distribución de probabilidad discreta perteneciente a ciertas variables aleatorias N que cuentan, con un número de ocurrencias discretas que toman lugar durante un intervalo de tiempo largo. La distribución de Poisson adquiere valores de X = 0, 1, 2, 3 ... y se utiliza a menudo como modelo para el número de eventos en un período de tiempo específico. ! )( k ekNP kλλ−== Donde: - e es la base del logaritmo natural (e = 2.71828), - k! es el factorial de k - λ es un número real positivo, igual al número espectado de casos que ocurren durante un intervalo dado. La distribución de probabilidad de la variable aleatoria de Poisson X que representa el número de resultados que ocurren en un intervalo de tiempo dado esta indicado por t: [ ] )exp( ! )(, t n ttnNP n ⋅−⋅== λλ ...(7) n = 0, 1, 2, ....; λ es el numero medio de ocurrencias por intervalo de tiempo, t puede ser 1 año. Para una mejor comparación de la distribución de Poisson con los números observados se recurre a la siguiente relación que permite graficar el número esperado K de años dentro de los cuales un numero de eventos n puede ocurrir. ],[],[ tnNPktnNK =⋅== ...(8) Donde k = numero de intervalos de tiempo usados. Estos conceptos y relaciones son usados en el presente estudio para estimar el riesgo sísmico para la región central del Perú. CATALOGO SISMICO Y TRATAMIENTO DE DATOS Un Catálogo Sísmico, se constituye, como una base de datos válida para realizar cualquier estudio de sismología. Esta base Aplicación de la distribución de poisson para el calculo del periodo de retorno de los sismos 107 tiene mayor validez cuando los parámetros que caracterizan a un sismo se calculan en las mismas condiciones y así constituir un catálogo sísmico homogéneo. A la vez deberá contener los siguientes parámetros sísmicos: fecha, hora, origen, latitud, longitud, profundidad, magnitud e intensidad. El catalogo sísmico utilizado en el presente estudio corresponde al Instituto Geofísico del Perú y comprende el periodo 1980 – 2000 para las latitudes -9.5º, -14.5º y longitudes -73º W, -81º W, todos con una magnitud mb > 3.5. En la Figura 5 se observa la distribución temporal de la sismicidad para el periodo de tiempo considerado, en la cual se nota que a partir del año 1994 se incrementa el registro de sismos esto debido a que en 1993 el IGP inicio el mejoramiento y ampliación de la Red Sísmica Nacional. Distribución Temporal de la Sismicidad 2 3 4 5 6 7 197 8 198 0 198 2 198 4 198 6 198 8 199 0 199 2 199 4 199 6 199 8 200 0 Tiempo en Años M ag ni tu d (m b) Figura 5. Distribución temporal de la sismicidad en la región de estudio (1980-2000). A partir de la distribución temporal se determina el “Umbral Mínimo” de magnitud para el catalogo elaborando gráficos de frecuencia sísmica para diferentes intervalos de magnitud (0.1, 0.25, 0.5, 0.75). Así, se consideró que para un intervalo de magnitud de 0.25 (Tabla 1), los datos muestran buena tendencia del número de sismos y magnitud coherente con el patrón sismotectónico de la región Central del Perú. Tabla 1. Valores utilizados en la curva frecuencia magnitud con un intervalo de mag. 0.25. MAGNITUD N Log(N) 3.5 218 2.34 3.8 134 2.13 4.0 156 2.19 4.3 97 1.99 4.5 174 2.24 4.7 119 2.08 5.0 38 1.58 5.3 12 1.08 5.5 9 0.95 5.7 2 0.30 6.0 1 0.00 J. Villegas 108 0 0.5 1 1.5 2 2.5 2 3 4 5 6 7 Magnitud (mb) Lo g (N ) M C = 4.5 a Figura 6. Distribución Frecuencia-Magnitud Catálogo del IGP (1980 - 2000). Constituida la curva de frecuencia - magnitud se observa una magnitud de completeza de 4.5, logrando seleccionar un total de 697 eventos. Para la determinación de los parámetros “a ” y ”b” se emplea el método gráfico según la Figura 6. El valor de “b” es representado por la pendiente de la curva y el valor de “a ” por la siguiente relación: bmNa 22log += ...(9) Los resultados de los parámetros son:: a = 9.3 y b = 1.55. Para el cálculo de estos parámetros también se emplea el Método de Mínimos Cuadrados, el cual realiza un ajuste a la pendiente de la distribución frecuencia – magnitud de sismos a partir del punto de máxima curvatura representado por la magnitud de completeza (MC) hasta la máxima magnitud de los datos sísmicos. Para este efecto se hace uso de las siguientes relaciones (Tabla 2): ( )( ) ( )∑ ∑ = = − −− = n i i n i ii LSE mm LogNLogNmm b 1 2 1 ...(10) mbLogNa LSELSE ⋅+= ...(11) donde: ∑ = = n i iLogNNN LogN 1 1 1 y ∑ = = n i ii NmN m 1 1 Para realizar el calculo por el método de mínimos cuadrados se resuelve la tabla 2, en la cual se presentan los resultados de los parámetros para calcular a y b: Tabla 2. Parámetros a utilizar en el calculo de las constantes a y b por el método de mínimos cuadrados. mi Ni Ni mi Log Ni Ni Log Ni mi – m (mi – m)2 Log Ni - Log N (Log Ni - Log N)* (mi - m) 4.5 174 783.00 2.24 389.86 -0.18 0.03 0.22 -0.04 4.7 119 559.30 2.08 246.99 0.02 0.0003 0.05 0.00 5.0 38 190.00 1.58 60.03 0.32 0.10 -0.45 -0.14 5.3 12 63.60 1.08 12.95 0.62 0.38 -0.95 -0.58 5.5 9 49.50 0.95 8.59 0.82 0.67 -1.07 -0.87 5.7 2 11.40 0.30 0.60 1.02 1.03 -1.72 -1.75 6.0 1 6.00 0.00 0.00 1.32 1.73 -2.03 -2.67 37 355 1662.8 8.23 719.02 3.91 3.94 -5.95 -6.05 Aplicación de la distribución de poisson para el calculo del periodo de retorno de los sismos 109 Luego de remplazar los valores obtenidos en la tabla anterior, en las ecuaciones 10 y 11 se obtiene los siguientes resultados: a = 9.19 y b = 1.54 Los resultados obtenidos por ambos métodos son bastante similares, si embargo el método por mínimos cuadrados presenta mejores limites de confianza (Wiemer y Wyss 1994), ya que permite obtener valores de la desviación estándar o sea los errores asociados a los parámetros de la ecuación de regresión lineal. Empleando las ecuaciones 3 y 4 se obtiene como resultado la siguiente expresión 12, con la cual se calcula el número anual de sismos para distintas magnitudes.: )(54.134.7 10.10 mN −= ...(12) El valor inverso de la ecuación 9 (ec. 5) permite obtener el periodo medio de retorno para distintas magnitudes. En la Tabla 3 y Figura 7 se muestra los periodos de retorno sismos de distinta magnitud que pueden ocurrir en la región de estudio. Tabla 3. Periodo de retorno. Mag (mb) N T(R) 4.5 2.57 0.4 5.0 0.44 2 5.5 0.07 13 6.0 0.01 79 0.1 10.0 1000.0 2 3 4 5 6 7 MAGNITUD (mb) Pe rio do d e R et or no Figura 7. Periodo de retorno para la Región Central del Perú. Luego haciendo uso de la ecuación 6 se calcula la probabilidad de ocurrencia. En la Figura 8 se presentan las curvas de probabilidad de ocurrencia obtenidas para sismos con magnitudes (mb) iguales a: 5.0, 5.5, 6.0, 6.5. 0 0.2 0.4 0.6 0.8 1 10 30 50 70 90 11 0 13 0 15 0 17 0 19 0 21 0 Años P( t) mb = 5.0 mb = 5.5 mb = 6.0 mb = 6.5 Figura 8. Probabilidad de ocurrencia de un evento sísmico (mb>5.0) para un periodo de 210 años. J. Villegas 110 Distribución de Poisson El análisis de la distribución de Poisson sugieren que para una MC = 4.5 es difícil obtener una buena distribución de Poisson, debido al alto numero de eventos de baja magnitud. A fin de evaluar la distribución de Poisson se procede a construir una base de datos con sismos de mb ≥ 5.0 y cuyos valores han sido mostrados en la Tabla 4. Ploteando estos datos, se obtiene una buena distribución normal; es decir, en forma de campana lo cual es una característica importante en una distribución de gauss (Figura 9). Tabla 4. Número de eventos por año de terremotos con mb ≥ 5.0 (1980 - 2000). Año N° Año N° Año N° 1980 3 1987 4 1994 3 1981 5 1988 2 1995 3 1982 4 1989 2 1996 1 1983 3 1990 3 1997 0 1984 3 1991 3 1998 0 1985 4 1992 5 1999 2 1986 4 1993 3 2000 1 0 5 10 0 2 4 6 n Fr ec ue nc ia Figura 9. Grafica de distribución normal de eventos. Con los datos de la tabla 4 se procede a calcular la distribución aleatoria P de Poisson aplicando las ecuaciones 7 y 8. Tabla 5. Valor P de la distribución de Poisson y Número esperado K. n Eventos P K 0 2 0.06 1.33 1 2 0.17 3.66 2 3 0.24 5.06 3 8 0.22 4.66 4 4 0.15 3.22 5 2 0.08 1.78 Con la relación de Poisson (ecuación 7), se calcula la probabilidad de ocurrencia de un evento para un periodo de tiempo determinado, esto expresado en porcentaje. 0 0.05 0.1 0.15 0.2 0.25 0.3 0 1 2 3 4 5 6 n P Figura 10. Distribución de Poisson. Probabilidad de que un evento tenga lugar en un intervalo de tiempo dado. En la Figura 10 se muestra la distribución de Poisson en forma de campana, similar a la campana Gausiana, y es válida para poder calcular la probabilidad de ocurrencia de un próximo sismo. Con los resultados obtenidos en la Tabla 5, se plotea los valores de K versus la clase (n) y se obtiene la curva de la Figura. 11, la misma que es similar a la obtenida en la distribución de Poisson, esto indica que los datos utilizados son confiables y por lo tanto, es posible hacer el calculo de probabilidad de ocurrencia. Aplicación de la distribución de poisson para el calculo del periodo de retorno de los sismos 111 0 1 2 3 4 5 6 0 1 2 3 4 5 6 n K Figura 11. Numero esperado K de años en los que n eventos ocurren. Estos resultados sugieren que la probabilidad de ocurrencia de un sismo de una determinada magnitud es constante en el tiempo. También es necesario considerar que los intervalos pequeños entre sismos son más probables que los largos; la probabilidad de que se den en forma simultanea es muy pequeña. Finalmente, usando la selección de datos y las relaciones que se muestran en el capitulo anterior se tiene los siguientes resultados: Los límites críticos de ocurrencia de un próximo sismo con magnitud de 5.0 y 5.5 mb son: para un sismo con mb = 5.0, la probabilidad de ocurrencia dentro de 10 años es del 98%, para un sismo con mb = 5.5, la probabilidad de ocurrencia dentro de 20 años es del 80% y para un sismo con mb = 6.0, la probabilidad de ocurrencia dentro de 55 años es de 50%. La distribución de Poisson muestra que para eventos con magnitud mb mayor a 5.0, estos obedecen a una buena distribución y por lo tanto, se pueden aplicar para la predicción de sismos. DISCUSIÓN E INTERPRETACIÓN Los sismos que se presentan en el Perú están directamente asociados al proceso de subducción de la placa de Nazca bajo la Sudamericana, a la vez se encuentran influenciados por la Fractura de Mendaña, la Dorsal de Nazca y el proceso de deformación cortical que ocurre en el continente. Estos procesos ocurren a diferentes intervalos de profundidad y han dado lugar a la formación de la Cordillera de los Andes. Según los datos empleados en el presente estudio, los valores hallados para las constantes a = 9.19 y b = 1.54, sirvieron para definir el número anual medio de sismos espectados de una magnitud determinada. Luego con este criterio se analizo el periodo de retorno teniendo así que los limites críticos de ocurrencia de un próximo evento son para magnitudes de 5.0 y 5.5 mb. La probabilidad de ocurrencia de sismos para la Región Central del Perú es: para una magnitud mb = 5.0, la probabilidad de ocurrencia dentro de 10 años es del 98%, para un evento mb = 5.5, la probabilidad de ocurrencia dentro de 20 años es del 80%, para un evento mb = 6.0, la probabilidad de ocurrencia dentro de 55 años es de 50%. Los resultados obtenidos para a y b son altos e indican que la región de estudio es propensa a la ocurrencia de sismos de baja Compendio de Trabajos de Investigación CNDG – Biblioteca Instituto Geofísico del Perú. V. 5 (2004) p. 113 - 122 EL CALCULO DE HIPOCENTROS: MÉTODO DE GEIGER Y ELABORACIÓN DE DROMOCRONAS IVAN LUIS CCUNO CHALLCO Escuela Profesional de Ingeniería Geofísica Universidad Nacional de San Agustín de Arequipa iccuno@axil.igp.gob.pe Practicas dirigidas por: Dr. Hernando Tavera Centro Nacional de Datos Geofísicos - Sismología RESUMEN En este estudio se realiza un análisis del método propuesto por Geiger en 1911 para la localización de hipocentros y de los parámetros que definen la ubicación de un sismo en el espacio. Este método es parte esencial de los diversos programas de calculo hipocentral utilizados frecuentemente: Hypo71, Hypoinverse, Fasthypo, entre otros. Asimismo, se construye las dromocronas para dos sismos de magnitud moderada ocurridos en Perú en Febrero del 2004. Estos sismos tienen curvas dromocronicas representadas por rectas, característica típica de sismos ocurridos a distancias regionales. INTRODUCCIÓN El interés por determinar el punto en el cual se produce un sismo (foco o hipocentro) y su proyección sobre la superficie (epicentro), se remonta a los primeros trabajos de investigación realizados por sismologos como Mallet y Omori al final del siglo XIX e inicios del siglo XX. Con el desarrollo de la instrumentación y la creación de observatorios sismológicos a nivel global, la determinación de los focos sísmicos se realizó con mayor precisión. El foco sísmico representa el punto inicial, en el espacio y tiempo, en donde se inicia la ruptura y se generan las ondas sísmicas. Los cuatro parámetros que definen la localización del foco en el espacio y en el tiempo son: la hora origen (t), las coordenadas geográficas (x, y) y la profundidad del foco (z). La localización de sismos, utilizando información de una red de estaciones, es usualmente formulada como un problema no lineal de mínimos cuadrados. La suma de los cuadrados de los residuales de los tiempos de llegada del valor teórico y el calculado para una serie de estaciones debe ser minimizada, lo cual conduce directamente al problema de la optimización no lineal. Los métodos de optimizaciones no lineales tienen muchas aplicaciones científicas y en el caso de la sismología el método propuesto por Geiger en 1911, resulta ser la herramienta más efectiva a la fecha. Este método permite reducir las diferencias existentes entre los tiempos de llegada teóricos y calculados para ondas sísmicas P y S, de manera que los residuales sean minimizados en un cierto sentido mediante mínimos cuadrados. Así, la localización de los sismos, como un problema de optimización no lineal, es resuelto. I. Ccuno 114 CALCULO DE HIPOCENTROS El primer problema en la sismología es determinar los parámetros de los sismos; es decir, el tiempo de origen, la localización del hipocentro, y su tamaño. Para determinar estos parámetros se dispone de N observaciones de tiempos de llegada (ti) de las ondas P y S a N estaciones de coordenadas (x, y), todas expresadas como (t0, x0, y0, z0); además se debe considerar un modelo de velocidad para dichas ondas. Si un sismo ocurre en un tiempo de origen t0 y el hipocentro se localiza en (x0 ,y0 ,z0), una serie de tiempos de llegada de ondas puede ser obtenido de una red de estaciones. Usando estos datos, mediante un procedimiento inverso, se calcula el tiempo de origen y el hipocentro del sismo. Este problema ha sido muy estudiado en sismología y a continuación, se describe como el método de Geiger permite dar solución al problema de la localización de sismos. La localización de sismos se realiza en un espacio de cuatro dimensiones: tiempo (t), y coordenadas espaciales x, y, z (latitud, longitud y profundidad), y en este espacio, un vector puede ser definido como (Lee, 1975): χ = (t , x , y , z )T (1) donde, el exponente T denota la transpuesta. En un espacio euclidiano n- dimensional se usa x para denotar un vector en que las coordenadas son x1, x2 , . . . .,xn, en reemplazo de las coordenadas espaciales x, y, y z. Para localizar un sismo se usa una serie de tiempos de llegada de las ondas τk a las estaciones en posiciones (xk , yk , zk) , k = 1, 2, . . . ,m. Asimismo, se debe asumir un modelo de tierra en la que los tiempos de viaje de las ondas teóricas TK de un hipocentro de ensayo (x*, y*, z*) puedan ser calculados. Dado el tiempo de origen y el hipocentro como un vector de ensayo χ* en un espacio euclidiano de cuatro dimensiones, se tiene: χ* = (t*, x*, y*, z*)T (2) Los tiempos de llegada de las ondas teóricas tk de χ* a la k-esima estación, representan el tiempo de llegada Tk mas el de ensayo para el tiempo de origen t* , tk (χ*) = TK (χ*) + t* para k = 1, 2, . . . ,m. (3) En este caso, tk no depende de t* y Tk se expresa como Tk(χ*) por conveniencia en la notación. Ahora, los residuales en la k- esima estación (rk), se define como la diferencia entre los tiempos de llegada observados y teóricos, rk = (χ*) = ( )*χτ kk t− = ( ) ** tTKk −− χτ para = 1, 2, . . . ,m. (4) MÉTODO DE GEIGER Geiger (1911) fue el primero en aplicar el método de Gauss-Newton para resolver el El calculo de hipocentros: método de Geiger y elaboración de dromocronas 115 problema de la localización de sismos, usando como alternativa la reducción por mínimos cuadrados, F (χ*) = [ ]∑ = m k kr 1 2*)(χ (5) donde, la residual rk(χ*) fue definido en la ecuación (6), siendo m el numero total de observaciones. Se puede considerar a los residuales rk(χ*), k = 1, 2, . . .,m como las componentes de un vector en un espacio euclidiano m-dimensional y puede ser expresado como: r = ( r1 (χ*), r2 (χ*),..., rm (χ*))T (6) siendo el ajuste del vector definido como, δχ =( δt , δx , δy , δz )T En el método de Gauss-Newton, la serie de ecuaciones lineales puede ser fue resuelto paso a paso para el ajuste del vector. En este caso, la serie de ecuaciones lineales puede ser escrito como: AT Aδχ = -AT r (7) donde, la matriz Jacobiana A esta definido por: A =               ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ z r y r x r t r z r y r x r t r z r y r x r t r mmmm MMMM 2222 1111 y las derivadas parciales están evaluadas en el vector ensayo (χ*) (Zill, 1988). Usando los tiempos de llegada residuales la matriz Jacobiana A se convierte en: A = -               ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ z m y m x m zyx zyx TTT TTT TTT 1 1 1 222 111 MMMM (8) El signo menos en el extremo derecho de la ecuación anterior surge de la definición de los tiempos de llegada de las ondas. Sustituyendo la ecuación (7) en la ecuación (8) y realizando la operación de la matriz, se tiene una serie de cuatro ecuaciones lineales simultaneas con cuatro incógnitas: G δχ = ρ (9) donde: G =           ∑ ∑ ∑∑ ∑∑ ∑ ∑ ∑ ∑ ∑∑ ∑ ∑ ∑ 2 2 2 iiiiii iiiiii iiiiii iii ccbcac cbbbab cabaaa cbam ρ = ( )Τ∑ ∑ ∑ ∑ kkkkkkk rcrbrar ,,, y la sumatoria ∑ es para k = 1, 2, . . .,m. Aquí se ha introducido, las siguientes expresiones: (10) I. Ccuno 116 Estas ecuaciones se refieren al sistema normal de ecuaciones para el problema de la localización de los sismos: dado una serie de tiempos de llegada de las ondas P y S, a partir de los cuales se calcula las derivadas del vector ensayo χ*, es posible resolver el vector δχ. Entonces se reemplaza χ* por χ* + δχ y se soluciona el problema no lineal por un procedimiento iterativo que implica solamente soluciones de una serie de cuatro ecuaciones lineales. En vez de solucionar un sistema uniforme determinado de cuatro ecuaciones lineales, se puede derivar un sistema determinado excedente equivalente para m ecuaciones lineales como: A δχ = - r (11) donde, la matriz Jacobiana A es una serie de m ecuaciones en la matriz. Entonces, las relaciones finales pueden ser expresadas como: (12) = rk (χ*) para k = 1 , 2 , . . . , m. Usando la inversión generalizada, según lo discutido en la sección anterior, el sistema de m ecuaciones puede ser resuelto directamente mediante los ajustes δt, δx, δy, y δz. CALCULO DE DROMOCRONAS Las dromocronas resultan de plotear el tiempo de llegada de la onda sísmica (minutos) versus la distancia recorrida por las ondas del epicentro a la estación (grados). Estas curvas permiten identificar las diferentes fases sísmicas registradas en un sismograma, siendo posible apreciar el diferente recorrido de las mismas cuando atraviesan el interior de la Tierra (Payo, 1986). Para la construcción de la dromocrona para un determinado sismo se tiene que identificar los tiempos de llegada de las diferentes fases sísmicas a cada una de las estaciones (Vilca, 2002), por ejemplo la Red Sísmica Nacional del Perú (CAM, QUI, PAR, GUA, ZAM, NNA, CVE, MIS, SGR, LYAR, CUS). Seguidamente, con los tiempos de llegada y el tiempo de origen del sismo se calcula ∆t, ∆t = ti – t0 (13) donde, ti : tiempo de arribo de la fase t0 : tiempo de origen del sismo Por ejemplo en las Tablas 1 y 2, se presentan los valores calculados para el tiempo t en minutos y la distancia d en grados para dos sismos ocurridos en Perú: el 24 y 25 de febrero del 2004. Los epicentros de estos sismos estuvieron ubicados en las localidades de Ocoña *** ,, χχχ z Tc y Tb x Ta KkKkKk ∂ ∂=∂ ∂=∂ ∂= **2 χχχ δδδ z Ty y Tx x Tt KKK ∂ ∂+∂ ∂+   ∂ ∂+ El calculo de hipocentros: método de Geiger y elaboración de dromocronas 117 (Arequipa) y Chilca (Lima). En este caso, t se ha convertido de segundos a minutos y los grados fueron normalizados a una distancia de 111.1 Km =1º. Tabla 1. Valores de t y ∆t de las ondas P para el Sismo Chilca-Lima del 24 de Febrero del 2004 Estaciones Sísmicas Tiempo Origen (hh:mm:ss.ss) Tiempo Arribo (seg.) t (minutos) Distancia "d"(grados) CAM 14:46:30.70 41.68 0.18 0.44 QUI 14:46:30.70 45.26 0.24 0.67 PAR 14:46:30.70 56.28 0.43 1.45 GUA 14:46:30.70 62.04 0.52 1.87 ZAM 14:46:30.70 70.95 0.67 2.52 NNA 14:46:30.70 43.80 0.22 0.55 Tabla 2. Valores de t y ∆t para las ondas P del sismo de Ocona-Caraveli-Arequipa del 25 de febrero del 2004 Estaciones Sísmicas Tiempo Origen (hh:mm:ss.ss) Tiempo Arribo (seg.) t (minutos) Distancia "d"(grados) CAM 20:27:5.22 92.97 1.46 5.73 QUI 20:27:5.22 76.75 1.19 4.72 PAR 20:27:5.22 67.32 1.04 3.98 GUA 20:27:5.22 60.01 0.91 3.51 ZAM 20:27:5.22 52.86 0.79 2.90 CVE 20:27:5.22 37.80 0.54 1.81 MIS 20:27:5.22 41.40 0.60 2.03 SGR 20:27:5.22 20.70 0.26 0.77 LYAR 20:27:5.22 56.60 0.86 3.14 CUS 20:27:5.22 63.10 0.96 3.57 Dromocronas para el sismo del 24 de Febrero del 2004 En las Figuras 1 y 2 se representan las dromocronas para las ondas P y S del sismo del 24 de Febrero del 2004 y en ellas se observa que las estaciones cercanas al epicentro del sismo se encuentran a una distancia de 0.4º (44 km) y la más lejana a 2.5º (277 km). A pesar del poco numero de datos, se observa que los valores se distribuyen sobre una recta, característica de sismos cercanos cuyos rayos se propagan en la corteza. 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 0.0 2.0 4.0 6.0 8.0 distancia (grados) tie m po (m in ut os ) Figura 1. Dromocronas de las ondas P para el Sismo de Chilca-Lima del 24 de Febrero del 2004. I. Ccuno 118 0.00 0.50 1.00 1.50 2.00 2.50 3.00 0.00 2.00 4.00 6.00 8.00 distancia (grados) tie m po (m in ut os ) Figura 2. Dromocronas de las ondas S para el Sismo de Chilca-Lima del 24 de Febrero del 2004. Dromocronas para el sismo del 25 de Febrero del 2004 Las Figuras 3 y 4 representan las dromocronas para el grupo de las ondas P y S correspondientes para el sismo del 25 de Febrero del 2004 y en ella se puede observar que la estación más cercana se encuentra aproximadamente a los 0.7º (77 km) y la más lejana a 5.7º (633 km). Al igual que el sismo anterior, las ondas viajan por la corteza. 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0 2 4 6 8 distancia (grados) tie m po (s eg un do s) Figura 3. Dromocronas de las ondas P para el Sismo de Ocoña-Arequipa del 25 de Febrero del 2004. Figura 4. Dromocronas de las ondas S para el Sismo de Ocoña-Arequipa del 25 de Febrero del 2004. Utilizando la información de las figuras anteriores, es posible estimar la variación de la velocidad en función de la distancia para cada fase sísmica, esto mediante el trazado de líneas tangentes que permitan construir una recta cuya inversa definirá la velocidad relativa de cada onda. Así, las velocidades relativas para la onda P del sismo de Chilca-Lima es de 4.2 km/s; mientras que, para la onda S es de 2.4 km/s. Para el sismo de Ocoña-Caraveli- Arequipa, la velocidad para la onda P es de 4.2 km/s y para la onda S, de 2.1 km/s. Estas velocidades para las ondas son menores a las definidas teóricamente, en razón que al ser los sismos cercanos a las estaciones sísmicas, las ondas viajan completamente por la corteza superior. RELACION VP/VS En la Figura 5, se considera una varilla de longitud original l y sección A que se extiende un incremento de longitud (l) debido a la aplicación de una fuerza tensional (F) en ambas caras (Figura 5a). 0.00 0.50 1.00 1.50 2.00 2.50 3.00 0.00 2.00 4.00 6.00 8.00 distancia (grados) tie m po (m in ut os ) El calculo de hipocentros: método de Geiger y elaboración de dromocronas 119 El módulo de elasticidad relevante es el denominado Módulo de Young (E), y esta definido por (Udias y Mescua, 1997): Modulo de young E = l l A F ∆ (14) Nótese que la extensión de esta varilla es acompañada de una reducción en su diámetro debido a que la varilla sufre deformación tanto longitudinal como lateral. La razón entre la deformación lateral y longitudinal es conocida como el coeficiente de Poisson (s), este valor oscila entre 0.05 y 0.40, siendo su valor medio para la mayoría de sólidos elásticos 0.25, aproximadamente. El módulo de compresibilidad (K) expresa la razón entre el esfuerzo y la deformación y para el caso de una presión hidrostática simple (P) aplicada a un elemento cúbico (Figura 5b), se define como: K νν∆ = P (15) De un modo similar el módulo de cizalla (µ) se define como la razón entre el esfuerzo de cizalla () y la deformación de cizalla resultante (tan ): modulo de cizalla = θ τµ tan = (16) El módulo axial () define la razón del esfuerzo longitudinal y la deformación longitudinal para el caso de que no hay deformación lateral (el material está constreñido a deformarse de un modo uniaxial): modulo axial = l l A F ∆= (17) Finalmente, la relación entre la deformación longitudinal ε1 y la deformación transversal εa , es conocida como el coeficiente de Poisson, l l a a a ∆ ∆ == 1ε εν (18) Figura 5. Esquemas que muestran la deformación de una barra metálica y un cubo. I. Ccuno 120 En un medio sólido, dos tipos de ondas fundamentales pueden propagarse: ondas compresionales (P) y ondas de cizalla (S). La velocidad de propagación de una onda compresional en cualquier material viene dada por el modulo elástico apropiado del material sobre la densidad del material: ρ µ=V (19) de este modo, la velocidad de las ondas P (VP) involucra una deformación compresiva que viene dada por el modulo axial sobre la densidad del material: ρ ψ=PV (20) con µψ .34+Κ= (21) La velocidad de la onda P es definida como: ρ µ.34+Κ=PV (22) Así, en la ecuación que define la velocidad de la onda P, E es la variable más importante ya que controla la velocidad de las ondas sísmicas en un medio, ( ) ( )( )ννρ ν 211 1 −+ −Ε=PV (23) y la velocidad de la onda S, que involucra deformación de cizalla pura viene dada por: ρ µ=sV (24) entonces,Vs es proporcional a la razón del esfuerzo longitudinal sobre un cuerpo a la extensión longitudinal producida (modulo de young) por una unidad incrementándose el coeficiente de poisson sobre la densidad del material. ( )νρρ µ +Ε== 121.SV (25) Ahora, si se relaciona las ecuaciones anteriores, se obtiene que las ondas P y ondas S son directamente proporcionales a los parámetros elásticos de los materiales e inversamente proporcional a sus densidades. ( ) ( )ν ν .21 1.2 − −= S P V V (26) y como el coeficiente de Poisson para las rocas consolidadas es típicamente del orden de 0.25, se obtiene un valor para la relación de velocidades de, 73.1= S P V V . (27) Este valor, es típico para una medio elástico, homogéneo e indefinido. El calculo de hipocentros: método de Geiger y elaboración de dromocronas 121 CONCLUSIONES El objetivo principal del método de Geiger es reducir las diferencias existentes entre los valores teóricos y los observados para los tiempos de llegada de las ondas sísmicas a fin de encontrar una solución satisfactoria para las expresiones de las series de Taylor alrededor del punto tor0, x0, y0, z0, definido como hipocentro inicial de ensayo. Después mediante un ajuste de mínimos cuadrados y minimizando los residuales de TP y TS teórico-observado. La solución puede ser evaluada mediante el RMS del ajuste. Todas las soluciones desarrolladas para el calculo de hipocentros sísmicos, consideran el método de Geiger (Hypoellipse, Hypo71, Hypoinverse, Fasthypo, etc.). Para los sismos de Chilca–Lima y Caraveli-Arequipa se ha construido las curvas dromocronicas para las ondas P y S. Considerando que los sismos son superficiales y ocurridos a distancias cortas, las domocronas son líneas rectas. La velocidad para la onda P es de 4.2 km/seg y para la onda S de 2.4 km/seg. La relación Vp/Vs =1.73, es equivalente a un coeficiente de Poisson de 0.25; es decir, para un medio homogéneo e isotópico. AGRADECIMIENTOS Mi agradecimiento al Dr. Hernando Tavera por su asesoramiento durante la realización del presente estudio. A Rocio Parillo por su ayuda en la redacción del texto. A todo el personal del CNDG- Sismología por su amistad y compañerismo que hicieron fructífera mi estancia en Lima. Al Instituto Geofísico del Perú por recibirme en su seno y permitirme realizar mis practicas pre- profesionales. BIBLIOGRAFÍA Kulhanek, O. (1990). Anatomy of Seismograms. Seismological Section. University of Uppsala. 125 p. Lee, W. (197 5). Generalized Inversion and No-Linear Optimization. Determination of Origin Time and Hypocenter, 105-135. Payo, G. (1986). Introducción al Análisis de Sismogramas. Instituto Geográfico Nacional, Madrid, 90 p. Udias, A. y Mescua, J. (1997). Fundamentos de Geofísica, Primera edición, UCM. Madrid, 420 p. Vilca, R. (2002). Estimación de Velocidades Relativas para Sismos de Gran Magnitud a partir de Curvas Dromocronicas. Informe de Practicas. Instituto Geofísico del Perú, Lima-Peru. Zill, D. (1988). Ecuaciones Diferenciales con Aplicaciones, México, 2-10. I. Ccuno 122 ANEXO PRINCIPALES FASES SÍSMICAS P, S Ondas longitudinales y transversales en el manto. Pg, Sg Ondas compresionales que viajan en la capa granítica de la corteza terrestre. PmP, SmS Ondas compresionales que se reflejan en la discontinuidad de moho. Pn, Sn Ondas compresionales que viajan a lo largo de la discontinuidad de moho. P*, S* (o Pb, Sb) Ondas compresionales que viajan a lo largo de la discontinuidad de conrad. PP, PPP, SS, SSS, etc, Ondas P y S reflejadas una, dos o más veces en la superficie de la tierra. PcP, ScS, PcP2, ScS2 , etc. Ondas reflejadas una dos o más veces en el núcleo externo. P’ o PKP Ondas reflejadas en el núcleo externo. SKS Ondas transversales transformadas en ondas longitudinales por refracción dentro del núcleo y transformadas nuevamente en onda transversal es por refracción desde el núcleo a la superficie terrestre de la tierra. PKKP, SKKS, PKKKP, Ondas longitudinales o transversales reflejadas una, dos, o PmKP (m-1) veces en la superficie interior del núcleo externo. P’P’ o PKPPKP Onda PKP reflejada en la superficie de la tierra del lado opuesto del foco. Suele ser visible cuando la P’ corresponde a la zona de la cáustica P’P’P’ Onda PKP reflejada dos veces en la superficie de la tierra en dirección del arco mayor del circulo máximo. L, Ondas superficiales sin especificar el tipo. G Ondas tipo love de gran periodo y amplitud que se propagan en el manto. LO o LQ Ondas superficiales Love. LR Ondas superficiales Rayleigh. Figura. Trayectoria de las ondas P y S por el interior de la tierra (Kulhanek, 1990 Compendio de Trabajos de Investigación CNDG – Biblioteca Instituto Geofísico del Perú. V. 6 (2004) p. 123 - 136 MODELAMIENTO MATEMÁTICO PARA LA PREDICCION Y POSICIÓN DE CUERPOS CELESTES DESDE CUALQUIER PUNTO DEL PLANETA JASON MÉNDEZ CÓRDOVA Escuela de Física Facultad de Ciencias Naturales y Matemática Universidad Nacional del Callao jmendez@axil.igp.gob.pe Prácticas dirigidas por: Dr. Mutsumi Ishitsuka Komaki Observatorio de Ancón RESUMEN En este trabajo de investigación se desarrollo los modelos y ecuaciones matemáticas que rigen el movimiento de los cuerpos celestes como: Planetas, asteroides, cometas, sol, lunas. Teniendo como centro de coordenadas nuestro planeta. En la astronomía siempre es necesario saber la posición de los cuerpos celestes en cualquier instante (hh:mm:ss). Con el modelo matemático hallado o desarrollado se genera un seudocodigo de programa para generar su software correspondiente llamado “Ephemeris” INTRODUCCIÓN Las fórmulas para estos cálculos pueden ser complicadas, pero si se toman en cuenta algunas simplificaciones es posible obtener unos cálculos sencillos con una exactitud de un minuto de arco = 1/60 grados. En la astronomía todos los cuerpos celestes están regidos por sus elementos orbitales que son: Consisten en 6 cantidades las cuales definen completamente una orbita circular, elíptica, parabólica o hiperbólica. Tres de estas cantidades describen la forma, tamaño y posición del planeta en la orbita: [Chauvenet, 1891] a: Distancia media e: excentricidad T: Tiempo en el perihelio Los tres restantes elementos definen la orientación de la orbita en el espacio: (Figura 1) i: Inclinación de la orbita con relación a la eclíptica. N: Longitud del nodo ascendente. W: Es el ángulo desde el nodo Ascendente al perihelio a lo largo de la orbita.(Figura 1) Figura 1. Elementos Orbitales J. Mendez 124 Existen otros elementos que se obtienen a partir de los valores anteriores que son usadas en los cálculos que son los siguientes: AU: Unidad Astronómica es la distancia media de la Tierra al Sol. Para describir la posición en la orbita, se usa tres ángulos Anomalía Media, Anomalía Verdadera y Anomalía Excéntrica. Ellas son cero cuando el planeta esta en el perihelio. M: Anomalía Media, este ángulo se incrementa uniformemente hasta 360º por período orbital. v: Anomalía Verdadera , es el actual ángulo entre el planeta y el perihelio como es visto desde el Sol. Este ángulo no se incrementa uniformemente, sino que cambia más rápidamente en el perihelio, Figura 2 Figura 2 . Anomalía verdadera E: Anomalía Excéntrica, este es un ángulo auxiliar usado por la ecuación de Kepler. Figura 3 Figura 3. Anomalía excéntrica w1 = N + w = Longitud del perihelio L = M + w1 = Longitud media q = a * (1 - e) = Distancia del perihelio Q = a * (1 + e) = Distancia del afelio P = a^1,5 = período orbital (en años) T = época de M - (M / 360) / P = Tiempo del perihelio Con estos datos se pueden aplicar las ecuaciones fundamentales de kepler para los movimientos de los planetas (Figura 1). GEOMETRÍA, SISTEMAS DE COORDENADAS Y EFEMÉRIDES Secciones y Ecuación Cónica La ecuación cónica puede considerarse de la Figura 4 de la siguiente manera: wxvre r e p ==− cos Entonces la ecuación será: )cos1/()cos1( veproverp +=−= También según la Figura q puede definirse en la coordenada positiva Xw donde Yw=0 (radio de la intersección de la cónica con la axisa Xw), r=q, tendremos (Figura 4, 5, 6): )1( eqp += Modelamiento matemático para la predicción y posición de cuerpos celestes desde cualquier punto 125 N(t-T)=M M=E-esinE )cos1( sin1 )(cos 2 Ear Eeay eEax w w −= −= −= Figura 4. Definición de directriz en una cónica Figura 5. Secciones cónicas de focos constantes Figura 6. Excentricidades (e) Ecuación de la Elipse, Parábola e Hipérbola Las coordenadas polares r y v no necesariamente están introducidos e las relaciones de tiempo y posición (en la ecuación de Kepler). Para esto será necesario introducir ciertas cantidades auxiliares D, E y F. En las Figuras 7a y 7b se muestran la generación de la anomalía media [Baker, 1967]:M que esta evaluado por )( Ttk − para la parábola y la anomalía media Mh. Figura 7a. Posición para un tiempo dado en una cónica arbitraria T t M 2/3/ akn µ= E Xw Yw r Figura 7b. Posición para un tiempo dado en una elipse La Elipse: Las orbitas elípticas, la cantidad auxiliar es E llamada anomalía excéntrica (Figura 8) y las ecuaciones de coordenadas referidas a la orbita esta definido por [Baker, 1967]: J. Mendez 126 )cos1( sin1 )(cos 2 Eear EeaY eEaX −= −= −= ω ω De la Figura tenemos: )(coscoscos eEaaeEaXvr −=−== ω Reemplazando con las ecuaciones (Figura 8): )cos1( )(cos)1( 2 Eea eEaeeaeXpr −= −−−=−= ω Finalmente EeaY EeaXrY w w sin)1( sin)1( 2 222222 −= −=−= ω Figura 8. Anomalía excéntrica De la ecuación de Kepler, orbita elíptica tenemos: 2 3 sin)()( a k n donde EeEttnMTtnM oo µ= −=−+=−= La Parábola: aquí la cantidad auxiliares D y esta definido por µ/rr & 222 2222 )( rsvr zzyyxxrrrr zyxrrs &&& &&&&& &&&&&& −= ++=•= ++=•= (Figura 9) Figura 9. Velocidad de componentes Relacionando en este caso y asiendo para el caso de e=1 y a=infinito, entonces D=(2q)0.5tan(v/2) entonces tendremos: 2 2 2 2 2 Dqr qDY DqX += = −= ω ω si análogamente lo relacionamos en la ecuación de Kepler estaría formulado por: 6 )()( 3DqD ttkMTtkM opop += −+=−∆ µµ La Hipérbola: para este caso asumimos como a negativo, las formulas elípticas pueden ser transformadas sustituyendo a E por iF y observando la siguiente convención 1−+=i y 11 22 −+=− eie . En este caso f esta definido con la definición de E Modelamiento matemático para la predicción y posición de cuerpos celestes desde cualquier punto 127 )cosh1( sinh1 )(cosh 2 Fear FeaY eFaX −= −−= −= ω ω (Ver Figura 10) Desarrollando en la ecuación de Kepler en el caso hiperbólico seria: 2 3 )( )sinh( )()( a k n donde FFe ttnMTtnM h ohhohh − = −= −+=−= µ Figura 10. Comparación de E y F SISTEMAS DE COORDENADAS Coordenadas esféricas horarias - Ecuador celeste: es el plano que pasa por el observador y es perpendicular al eje del mundo. - Círculos horarios: Son los círculos máximos de la esfera celeste que pasan por los polos celestes. - Paralelos celestes: Son los círculos menores paralelos al Ecuador . Las coordenadas de un punto A (x' , y', z') se llaman rectilíneas horarias (Figura 11) Figura 11. Sistemas de Coordenadas -Ángulo Horario H: es el arco de ecuador medido desde la recta de intersección del plano meridiano y el plano del ecuador, 0h. a 24h. -Declinación D: Es el ángulo medido desde el círculo del astro al Ecuador. Positiva hacia el Norte y negativa hacia el Sur Varia pues entre 90º y -90º - Distancia Polar P: Es el complementario de la Declinación. Por tratarse de coordenadas esféricas se cumplirá:         =         senDr senHDr HDr Z Y X ' .cos' cos.cos' ' ' ' J. Mendez 128 - Ascensión (δ ) y Declinación recta (α ): Esta relacionado con las coordenadas rectangulares. (Figura 12) x = r cosδ cosα y = r cosδsen α z = r senδ Figura 12. Ascensión recta y declinación Recta Clasificación de los sistemas de referencia Los sistemas de referencia se clasifican según la elección del origen en: Coordenadas topocéntricas: Centradas en el observador. Coordenadas geocéntricas: Centradas en el centro de la Tierra. Coordenadas heliocéntrica: Centradas en el centro del sistema Solar. Coordenadas galácticas: Centradas en el centro de la Galaxia. Atendiendo a que sus valores dependan o no de la posición del observador las coordenadas se clasifican en : Locales : Coordenadas Horizontales y Horarias No Locales: Coordenadas Ecuatoriales , Eclípticas, Galácticas ESCALA DE TIEMPO La escala del tiempo está formulada para contar días. Las horas, minutos y segundos son expresados en fracción de días. El día cero comienza el 31 de diciembre de 1999 a las 12:00 am UT (Universal Time, es decir, Tiempo Universal). d= 367*y - int( (7 * (y +(int( (m + 9) / 12) ) ) ) / 4) + int( 275*m / 9 ) + D - 730530 +UT / 24 [Serafino Zani – Milano 1996] Donde y = año (cuatro dígitos); m = mes; D = día y UT = en horas + decimales. int () es una función que sólo toma la parte entera de la división Algunos datos de Elementos Orbítales Elementos Orbítales del Sol (Tabla 1): N=0i=0 w=282,9404+4,70935E-5*d a=1(AU) e=0,016709-1,151E-9*d M = 356,0470 + 0,9856002585 * d Elementos Orbitales de la Luna: N=125,1228-0,0529538083*d i=5,1454 w=318,0634+0,1643573223*d a=60,2666(radioTerrestre) e=0,0549 M = 115,3654 + 13,0649929509 * d Elementos Orbitales de Marte: N=49,5574+2,11081E-5*d i=1,8497-1,78E-8*d w=286,5016+2,92961E-5*d a=1,523688(AU) e=0,093405+2,516E-9*d M = 18,6021 + 0,5240207766 * d Modelamiento matemático para la predicción y posición de cuerpos celestes desde cualquier punto 129 TIEMPO SIDERAL Y ÁNGULO HORARIO. AZIMUT Y ALTITUD TIEMPO SIDERAL (SIDTIME) Tiempo Sideral en el Meridiano de Greenwich (GMST) y el Tiempo Universal (UT), todos en horas + decimales, son necesarios calcularlos para obtener la altitud y azimut en nuestro sitio de observación. La Altitud y Azimut, cambian a medida que transcurre el tiempo y depende del sitio geográfico del observador. SIDTIME = GMST + UT + Longitud / 15 GMST = ang(L + 180) / 15 = ang( ang(427,9895482) +180) / 15 = ang(67,98954819 + 180) / 15 GMST = ang(247,9895482) / 15 = 247,9895482 / 15 = 16,53263655 UT = Hora local - Diferencia Horaria = 11,75 - (-4) = 15,75 SIDTIME = 16,53263655 + 15,75 + (- 66,9166666667) / 15 = 27,8215254389 Si SIDTIME es negativo, sumamos 24h. Si es mayor a 24h, entonces restamos 24h SIDTIME = 27,8459698833 - 24 = 3,8215254389 Ángulo Horario (HA) incrementa con el tiempo (a menos que se mueva más rápido que la rotación de la Tierra, como por ejemplo, satélites). SIDTIME y RA tienen que estar en la misma unidad, por tanto hay que pasar SIDTIME a grados, multiplicándolo por 15. HA = SIDTIME - RA Entonces, se traslada estos datos a coordenadas rectangulares, donde el eje X apunta al ecuador celeste, el eje +Y apunta hacia el horizonte oeste y el eje Z apunta al polo norte celeste. x = cos(HA) * cos(Decl) y = sin(HA) * cos(Decl) z = sin(Decl) Se retorna este sistema de coordenadas xhor = x * sin(lat) - z * cos(lat) yhor = y zhor = x * cos(lat) + z * sin(lat) Finalmente, se pasa del sistema de coordenadas rectangulares al sistema esférico donde el radio puede ser la unidad. Azimut = atan2(xhor; yhor) + 180º Altitud = asin(zhor) ALGUNOS SEUDOCODIGOS Y CÁLCULOS Cálculo de la anomalía media del planeta La anomalía media del planeta se calcula mediante la siguiente fórmula: M = n * d + L - v n es el movimiento diario d es el número de días desde la fecha de los elementos L es la longitud media v es la longitud del perihelio Cálculo de la anomalía verdadera del planeta [Lacruz-Mayo, 1993] v = M + 180/pi * [(2 * e - (e^3) /4) * sin (M)+ 5/4 * e^2 * sin (2*M)+ 13/12 * e^3 * sin (3*M)] v es la anomalía verdadera M es la anomalía media J. Mendez 130 e es la excentricidad pi es 3.14159... Calculo del radio Vector del Planeta [Lacruz-Mayo 1993] r = a * (1 - e^2) / [1 + e * cos (v)] a es el semieje mayor e es la excentricidad v es la anomalía verdadera Tabla 1. Datos de algunos elementos Orbitales ( a, e y T) Planeta a: semieje mayor e: excentricidad T: periodo Mercurio 0,387 0,206 87,97d Venus 0,723 0,007 224,8d Tierra 1,000 0,017 365,3d Marte 1,524 0,093 1,881a Asteroides (Ceres) 2,768 0,078 1681,6d Júpiter 5,204 0,048 11,881a Saturno 9,575 0,052 29,458a Urano 19,31 0,050 0,773 Neptuno 30,2 0,004 164,79a Plutón 39,91 0,257 248,4a APLICACIÓN Y SEUDOCODIGO DE PROGRAMA “EPHEMERIS DE LA LUNA” CÁLCULO DE LA POSICIÓN DE LA LUNA Para calcular la Posición de la Luna con exactitud hay que considerar cientos de términos periódicos [Paz Soldán , 1898], por esta razón es importante saber que términos hay que considerar para cometer errores de unos 10" en la longitud, 3" en la latitud y 0",2 en el paralaje. Conocido el paralaje la distancia a la Luna D se puede obtener en radios ecuatoriales mediante: Sin P=1/D. El método aquí descrito está tomado de 30 posiciones de la luna descrita por el método numérico de Jean Meuus. Los siglos julianos de 36525 días, transcurridos desde 0,5 de Enero de 1900 hasta la fecha son: T = (JD- 2415020) / 36525 Sus cuadrado y cubo: T2 = T * T; T3 = T2 * T Como T se expresa en siglos hay que tomar en los cálculos un número elevado de decimales. La longitud media de la Luna: L1= 270.434164 + 481267.8831 * T - .001133 * T2 + .0000019 * T3 La anomalía media del Sol: M = 358.475833 + 35999.0498 * T - .00015 * T2 - .0000033 * T3 La anomalía media de la Luna: M1 = 296.104608+477198.8491 *T + .009192 * T2 + 0000144 *T3 La elongación media de la Luna: D=350.737486 + 445267.1142 * T - .001436 * T2 + .0000019*T3 Modelamiento matemático para la predicción y posición de cuerpos celestes desde cualquier punto 131 La distancia media de la Luna desde el Nodo ascendente: F = 11.250889 + 483202.0251 * T - .003211 * T2 - .0000003 * T3 La longitud del nodo ascendente de la Luna: OM=259.183275- 1934.142*T+.002078*T2+0000022*T3 A estos valores medios deben añadirse algunas variaciones periódicas, llamados Términos Aditivos L1 = L1 + .000233 * sin((51.2 + 20.2 * T) * PI / 180) M=M-.001778*sin((51.2 + 20.2 * T)*PI/180) M1=M1+.000817*sin((51.2+20.2*T)*PI/1 80); D = D + .002011 * sin((51.2 + 20.2 * T) * PI / 180) Estos cuatro términos tienen un periodo de 1782 años. El término: S = .003964 * sin((346.56 + 132.87 * T - .0091731 * T2) * PI / 180) El día Juliano en un método numérico seria representado por: Fecha juliana Hora = Hora + (Minutos / 60); GGG = 1; if (YY <= 1585) GGG = 0; JD = -1 * entero(7 * (entero((MM + 9) / 12) + Año) / 4); S = 1;if ((Mes 9)<0) S=-1; A = abs(Mes 9); J1 = entero(año + S * entero(A / 7)); J1 =-1* entero((entero(J1 / 100) + 1)*3/4); JD = JD + entero(275 * Mes / 9) + DD + (GGG * J1); JD = JD + 1721027 + 2 * GGG + 367 * YY - 0.5; JD = JD + (Hora / 24); J2=JD; COMPARACIÓN DE DATOS CON ALMANAQUES NÁUTICOS Los datos mostrados en la tabla muestra la comparación entre los datos de los almanaques Náuticos, o datos de efemérides de la Nasa con los datos obtenidos usando el software desarrollado “Ephemeris”, siendo el % de error de 0.41225502%[Tabala 2][Young, 1989] J. Mendez 132 Tabla 2. Datos de Ephemeris y su comparación (error%) 27-12-2004 (Planeta Mercurio) 12:00 m Minutos Nasa Ephemeris Porcentaje (%) 2 2.707 2.7823 2.706394 3 4.917 4.9924 1.5102957 4 7.127 7.2024 1.0468733 5 9.337 9.4124 0.8010709 6 11.547 11.6223 0.6478924 7 13.757 13.8322 0.543659 8 15.967 16.0421 0.4681432 9 18.177 18.2519 0.4103682 10 20.386 20.4617 0.3699595 11 22.596 22.6715 0.3330172 12 24.806 24.8812 0.3022362 13 27.015 27.0909 0.2801679 14 29.225 29.3006 0.2580152 15 31.435 31.5102 0.2386529 16 33.644 33.7198 0.2247937 17 35.854 35.9293 0.2095783 18 38.063 38.1388 0.1987477 19 40.273 40.3483 0.186625 20 42.482 42.5578 0.1781107 21 44.691 44.7672 0.1702139 22 46.901 46.9765 0.1607187 23 49.11 49.1859 0.1543125 24 51.319 51.3952 0.1482629 25 53.529 53.6045 0.1408464 26 55.738 55.8137 0.1356298 27 57.947 58.0229 0.1308104 % Error 0.412255 Nota: Utilizando mínimos cuadrados se puede hallar la ecuación de la grafica; para ambos casos las ecuaciones tienes una ligera variación [Figuras 13a,13b] Almanaque “Nasa” y = 2.1009x + 0.287 Programa “Ephemeris” y = 2.103x + 0.324 Datos Ephemeris Nasa y = 2.1009x + 0.287 0 20 40 60 80 0 5 10 15 20 25 30 Minutos se gu nd os Figura 13a. Grafica comparativa de datos de Ephemeris Nasa y Software Ephemeris Datos "Ephemeris-Programa" y = 2.103x + 0.324 0 20 40 60 80 0 5 10 15 20 25 30 M inutos Figura 13b. Grafica comparativa de datos de Ephemeris Nasa y Software Ephemeris RESULTADOS Y DISEÑOS DEL SOFTWARE El conjunto de figuras e imágenes que se presenta a continuación, se muestra el diseño del Software “Ephemeris” para el calculo y predicción de cuerpos celestes, el mismo que está en desarrollo, siendo el presente estudio parte de los resultados obtenidos a la fecha. Modelamiento matemático para la predicción y posición de cuerpos celestes desde cualquier punto 133 Diseño Detalle El software calcula la predicción de Posición de cuerpos celestes, así como también predicción de eclipses solares y Lunares, posición de satélites galileanos, sol, Luna, Planetas, Asteroides , cometas, etc. Esta imagen muestra el calculo de la posición del planeta Marte en coordenadas Geocéntricas (dia/mes/año– hh:mm) J. Mendez 134 Este es el diseño para el cálculo de la predicción de Eclipses Solares, mostrando los datos necesarios. Esta imagen muestra el lugar por donde se podrá observar el eclipse solar 8 abril 2005 (imagen q muestra el software “Ephemeris”) Esta imagen muestra el diseño del calculo y predicción de Eclipses Lunares, mostrando los datos necesarios Modelamiento matemático para la predicción y posición de cuerpos celestes desde cualquier punto 135 Esta imagen muestra la trayectoria por la umbra y penumbra del Eclipse Lunar del 24 de abril 2005 Esta imagen muestra el calculo y el tiempo de duración de las estaciones del año Esta imagen muestra el calculo y predicción de las fases lunares para un mes dado El programa desarrollado se le ha llamado Ephemeris y en él se detalla la posición exacta de cada planeta en coordenadas geocéntricas, ecuatoriales, heliocéntricas y J. Mendez 136 topo céntricas a una fecha (día, mes, año), hora y minutos Nota: El trabajo de investigación presentado aun esta en proceso y solo se muestra un avance del mismo. AGRADECIMIENTOS Primero a Dios que me da fuerzas cada día de mi vida, a mis padres José Morales y Zenobia Córdova que siempre me apoyan y me dan aliento, a mis hermanos Jose Eduardo y Luis Alfredo que siempre me alientan, al Dr. Mutsumi Ishitsuka Director del Observatorio de Ancón por darme la oportunidad de hacer mis investigaciones, darme su apoyo y comprensión, a mi gran amor Edin Milagros Chang que siempre esta ahí para apoyarme, y a todo el personal que trabaja e investiga en el Observatorio de Ancon BIBLIOGRAFIA Baker, R. (1967) Astrodynamics, Second Edition, Computer Sciences Corporation and University of California , Los Angeles. Chauvenet, W. (1891).Manual of Spherical and Practical Astronomy-Vol I. Paz Soldan, M. (1848).Tratado Elemental de Astronomía – Teórica y Práctica, Volumen I. Paris. Young, Ch. (1989). Nautical Almanac 1990-2003, Issued by the Nautical Almanac Office United States, Naval Observatory Paginas Web utilizadas y consultadas -http://asteroidi.uai.it/manuale.pdf Observatorio Astronómico Serafino Zani – Guía de Observaciones -http://www.terra.es/personal6/achernar- /efemerides/efemerides.htm Asociación Astronómica Sirio de Pontevedra http://www.astrogea.org/foed/efemerides/e femerides.htm Efemérides Astronómicas 2004 Softwares Utilizados para la comparación de datos: - Software “Total Eclipse” versión 2.0 , Zephyr Services 1900 Murray Avenue Pittsburg. PA 15217 - Software “Mica” Multiyear Interactive Computer Almanac – Version 1.50. U. S. Naval Observatory Compendio de Trabajos de Investigación CNDG – Biblioteca Instituto Geofísico del Perú. V. 6 (2004) p. 137 - 146 DESCRIPCIÓN DEL MÉTODO DE PROYECCIÓN PARA EL CONTEO DE MANCHAS SOLARES EN LA ESTACIÓN SOLAR UNICA-IGP ADITA NEHEMIAS QUISPE QUISPE Facultad de Ciencias Escuela Académica Profesional de Física Universidad Nacional San Luis Gonzaga de Ica aditaquispe@yahoo.es Prácticas dirigidas por: Mutsumi Ishitsuka Komaki Observatorio de Ancón RESUMEN El objetivo de este trabajo es describir el Método de Proyección para el conteo de manchas solares realizados por los estudiantes de la Universidad Nacional San Luis Gonzaga de Ica en convenio con el Instituto Geofísico del Perú, para ello se utilizó un telescopio refractor instalado en las instalaciones de la universidad de Ica. Para los registros de manchas solares realizados en la Estación Solar UNICA-IGP, se tomo las consideraciones de ubicación y orientación del telescopio de acuerdo a la hora de registro de los datos, este método de proyección para el conteo de manchas solares tiene por objetivo realizar el dibujo y conteo de manchas solares, que por medio del uso del telescopio con un ocular apropiado nos permite proyectar la imagen ampliada del Sol a 15 centímetros de diámetro y a través de las transparencias determinaremos sus coordenadas heliográficas, longitud relativa (E-W), longitud verdadera (L) y latitud(φ). INTRODUCCION La superficie del sol conocida como fotosfera es una capa delgada de solo 550 Km. de espesor, su temperatura efectiva es de 6000K, (Priest, 1982). La fotosfera tiene una serie de características propias, nos referimos a: *La granulación, que es un efecto de la parte más externa de las celdas convectivas, tienen un tamaño que varía entre los 700Km y 1500Km, Estas celdas traen la energía solar del interior del sol, (Priest, 1982). *Las Fáculas fotosferica son una estructura en forma de fibras más brillantes que la fotosfera, rodeando los bordes de las granulaciones *En fotografías de alta resolución podemos observar puntos oscuros que se caracterizan por contener campos magnéticos muy intensos (Priest, 1982). La actividad solar se hace evidente en la Fotosfera con la presencia de manchas solares que tiene un ciclo de 11 años. Las manchas solares es un fenómeno conocido desde mucho tiempo atrás (inclusive antes de la invención del telescopio), y fue redescubierto por Galileo en sus observaciones en 1610 (Kaufmann, 1994) Las manchas solares tienen una considerable influencia en las condiciones de nuestro planeta. Los periodos de gran actividad solar tienen efectos en la perturbación del campo magnético terrestre y provocan tormentas geomagnéticas. Si las tormentas solares están asociadas a una eyección de masa coronal (CME), pueden que ocasionen apagones en plantas eléctricas y daños en los satélites que orbitan la Tierra. El estudio y análisis de esta actividad es muy importante, dado que su influencia sobre la atmósfera terrestre se considera cada vez más significativa. A. Quispe 138 Para el estudio de las manchas solares existe diversos métodos de los cuales describiremos el método de proyección que consiste en dibujar cada estructura que presenta el sol en una imagen del disco solar proyectada por el telescopio de 15 cm. de diámetro En la actualidad este método es considerado estándar y es usado por muchos observatorios solares en el mundo, con el objeto de contabilizar las manchas solares observados en el disco solar y obtener el índice de la actividad solar que se denomina Número Relativo de Manchas Solares MANEJO DEL TELESCOPIO. Figura 1. Telescopio Takahashi, proyectando la imagen del disco solar A continuación describiremos los detalles mas importantes que debe de tenerse en cuenta para un manejo adecuado del telescopio. El telescopio con el que se ha trabajado es un refractor marca Takahashi, modelo FCT150 con las siguientes características: • Abertura del lente objetivo principal es de 15cm.de diámetro. • Su distancia focal es de 1050mm. • Instalado sobre una montura ecuatorial modelo EM-500. Durante su manipulación debemos tener bastante cuidado al momento de observar con el telescopio. Nunca se debe de observar directamente al sol por que quedaríamos totalmente ciegos. Para el manejo del telescopio debemos considerar los siguientes puntos: 1. Retirar la tapa exterior, sin tocar el lente 2. Observaremos al lente objetivo, e inmediatamente colocamos el anillo de reducción para reducir la abertura del telescopio de 15cm a 10 cm. Por medio de su rosca que nos permite colocarlo por encima del objetivo. 3. Colocar el soporte y luego la hoja para proceder a graficar lo que se observa en la proyección (ver Figura 1 y 2). 4. Con la ayuda del ocular del tipo ortoscópico de 25mm. de distancia focal podemos proyectar la imagen del sol del tamaño adecuado. Figura 2. En detalle la imagen del disco solar sobre la pantalla CALIBRACIÓN Debemos estar seguros que el telescopio se encuentre perfectamente operativo y adecuadamente calibrado antes de realizar las observaciones y registrar los datos. Para esto el eje principal de la montura ecuatorial debe orientarse paralelo al eje polar o eje de rotación de la Tierra Por otro lado debemos asegurarnos que el telescopio con todas las partes a usarse debe de estar en equilibrio, para ello usaremos el sistema de contrapesos y Descripción del método de proyección para el conteo de manchas solares en la Estación Solar única-IGP 139 debemos mover la posición del telescopio con mucho cuidado para este propósito. DESCRIPCIÓN DEL MÉTODO DE PROYECCIÓN El método consiste en el uso de un telescopio, que por medio de un ocular apropiado, nos permite proyectar la imagen del disco solar. La imagen que observamos corresponde a la fotosfera solar que tiene la forma de un disco blanco brillante, cabe resaltar que la imagen esta invertida (Figura 2). Teniendo en cuenta estas consideraciones se podrá observar las machas solares mediante el método de proyección. Considerando los siguientes pasos: Primer paso.- Orientar el telescopio al sol, y con ayuda del buscador o del ocular nos guiaremos por la sombra que proyecta hay que dirigir el telescopio al sol y moverlo hasta que su sombra sea circular, luego realizar pequeños ajustes en el foco con la finalidad de obtener la mejor calidad de imagen. Dibujar cada estructura que se observa con la mayor fidelidad posible, no debemos olvidar de graficar las fáculas si se observa alguna (Figura 2). Segundo paso. A continuación, se escoge alguna pequeña mancha y se le marca por medio de un punto. Paramos el motor de seguimiento y esperamos un tiempo para que la imagen del Sol se desplace a causa del movimiento diurno, Entonces marcaremos otra vez con un punto la posición final de la mancha considerada. Si nosotros trazamos una recta entre los puntos que define el momento inicial y final, estaremos describiendo la dirección E-W producida por la rotación de la tierra en el punto de observación, como muestras las Figuras 3 y 4 Trazar una recta paralela a la línea que define la dirección E-W. por el centro del disco solar, para luego trazar su perpendicular que cruce el centro del disco, esta nueva recta representa el eje de rotación de la tierra. Figura 3. Seguimiento de una mancha para orientar el N-W terrestre Figura 4. Nos muestra el N-S y el E-W terrestre Tener en cuenta que si por algún motivo se mueve la hoja de recolección de datos, esta ya no se debe considerar, es recomendable volver a realizarlo. DETERMINACIÓN DE LAS COORDENADAS HELIOGRÁFICAS USANDO LA PLANTILLA El primer paso es conocer los valores de P, Bo, L0 que serán explicados más adelante. Ahora, rotamos el valor del ángulo P que nos permita establecer la dirección del eje de rotación del Sol, téngase en cuenta que si P es positivo rotamos en el sentido de A. Quispe 140 las agujas del reloj, en caso contrario será en dirección opuesta. La forma más sencilla de medir la posición de las manchas es utilizar un conjunto de plantillas sobre el dibujo de manchas solares. El valor de B0 define la plantilla a usarse. El meridiano de la plantilla debe coincidir con el eje de rotación solar y, por ello, debemos orientar la imagen de la manera más precisa posible. Usamos el dibujo, de las manchas solares por el método de proyección y en él marcaremos la posición heliográfica de las manchas solares (García La Rosa, 1989). Por lo general, se considera el centro de las manchas para determinar la ubicación de ellas. Figura 5. Plantilla que nos permite ubicar mejor las coordenadas heliográficas de cada grupo de manchas observada. Posteriormente buscaremos en el MICA1 el valor de P (ángulo de posición entre el eje terrestre y el eje solar), y llevándolo sobre el dibujo lo tendremos orientado. 1MICA (Multiyear Interactive Computer Almanac). es un sistema de software que provee de datos astronómicos de alta precisión calcula en tiempo real, en un intervalo de 15 años. Para fijar la posición de un punto sobre la superficie terrestre utilizamos dos coordenadas: la latitud y la longitud. La latitud se mide desde el ecuador hacia los polos y desde 0º hasta 90º. Para la longitud, Sin embargo, no hay que perder de vista que en el Sol no hay estructuras estáticas y que, por lo tanto, lo normal es que las coordenadas y las formas van variando con el tiempo. Se consideraron las variables heliográficas: g (grupo de manchas), f (número de manchas), T (tipo de mancha), Φ (latitud relativa), λ (longitud relativa), L (longitud heliocéntrica), la visibilidad, la fecha y la hora de la toma de datos. Si los ejes de rotación del Sol y la Tierra fuesen paralelos entre sí y perpendiculares a la eclíptica, el eje de rotación solar siempre se vería en dirección Norte-Sur y el ecuador solar sería una línea recta que pasaría por el centro del disco aparente. ANÁLISIS DE LOS DATOS OBSERVADOS Ahora describiremos el trabajo que debe realizarse con cada una de las imágenes del sol obtenidas por los medios ya descritos. Hay que tomar consideración a los siguientes datos: - Fecha de observación. - Hora de observación. - Visibilidad: ™ Pobre 1 ™ Regular 2 ™ Excelente 3 - Nombre del observador. - Los valores de P. B0 y L0 1. Normalizar el valor de las coordenadas para poder simplificar por un lado nuestro sistema de unidades y poder referirnos en términos de radio solar. 2. Usando el valor del ángulo P, que define el ángulo que existe entre el eje de rotación de la tierra con respecto a la del sol. Puesto que nuestro interés es ver el Descripción del método de proyección para el conteo de manchas solares en la Estación Solar única-IGP 141 movimiento de rotación del sol, debemos de rotar el valor de las coordenadas (x ; y) a las nuevas coordenadas solares (x’ , y’), usando para ello la siguiente matriz de rotación:     −= psenp senpp A cos cos Entonces, ( ) ( )    −= psenp senpp yxyx cos cos '' O lo que es lo mismo que (x’ , y’)= (x, y)A; donde A es la matriz de rotación. El valor de p lo buscaremos en las efemérides. El origen de este ángulo se debe a la inclinación del eje de rotación terrestre con respecto a la eclíptica varia entre +26.4° y -26.4°, siendo el valor positivo cuando esta en sentido N-E. Figura 6. Valores de P durante un año, desviaciones del eje solar respecto a la dirección Norte-Sur. 3. Los valores del centro del disco solar representa al ángulo Bo, en el eje vertical y el ángulo L0 en el eje horizontal, se define el ángulo Bo como el valor de la latitud heliográfica en el centro del disco solar, este valor varia debido a la inclinación del eje de rotación solar con respecto de la eclíptica, el valor cambia entre +7.2° y -7.2°. 4. El valor L0 nos define el valor de longitud heliocéntrica en el centro del disco solar. La coordenada L0 esta relacionada directamente con la rotación del sol y definido entre 360° y 0°. Aproximadamente cada 27 días debido a la rotación del sol, se inicia un nuevo ciclo es decir una nueva rotación solar, de acuerdo a la rotación de Carrington. Debemos subrayar que la rotación de Carrington corresponde a la rotación siderea que es un poco menor que la rotación sinódica. Se define la rotación siderea como el periodo de rotación del sol, suponiendo que nos encontramos en el centro de masa del sol. Ejemplo de una hoja de recolección de datos para el registro de mancha solares en la Estación Solar UNICA-IGP. Compendio de Trabajos de Investigación CNDG – Biblioteca Instituto Geofísico del Perú. V. 6 (2004) p. 137 - 146 ESTACIÓN SOLAR UNICA-IGP FECHA:…………………………. HORA:………(UT) VISIBILIDAD:..………………....OBSERVADOR: ……………………..…………… λ longitud relativa (E-W) L longitud verdadera Φ Declinación Figura 7. Formato de Recolección de datos Esta hoja es colocada en el tablero para dibujar cada estructura, como se muestra en la Figura 2. Cuando nos referimos al Número relativo de manchas Solares, nos referimos a la famosa ecuación definida por Rudolf Wolf en 1868: R = k (10f + g) Donde: • R, representa el número relativo de manchas solares, • f, es el conteo total de manchas observadas, • g, es el número de grupo de manchas observadas, • k, es la constante de cada observador. La ecuación Rudolf Wolf, permite comparar los valores generados por la misma formula con otros observatorios mundiales en este caso con el International Relative Sunspot Number de (Brussela). Datos considerados (Figura 7) 1. Los datos fueron tomados en la Estación Solar UNICA-IGP Se ha considerado las observaciones desde el 24 de noviembre del 2003 al 24 de marzo del 2004. 2. El objetivo seria poder contar con el número relativo de manchas solares. Considerando la distribución latitudinal de las manchas solares, y el conteo total de manchas. Ya realizados los dibujos mediante el método de proyección, el observador va contando el número de grupos de manchas solares y de este modo se obtienen los siguientes elementos: • La apariencia que tuvo el disco solar en el momento de su Ica Nº f tipo Φ λ L P=………… B0=………. L0=………. Descripción del método de proyección para el conteo de manchas solares en la Estación Solar única-IGP 143 observación que se sintetiza por el número relativo de manchas solares. • La posición aparente de las manchas en el disco solar. • Número de manchas de cada grupo. ™ Norte: azul. ™ Sur: rojo. • Latitud heliográfica del grupo de mancha φ, donde: ™ Φ > 0 Norte ™ ™ Φ < 0 Sur • Tipo de mancha basado en la clasificación de Zurich. • Número contado de manchas en cada grupo. • El resultado de 10g+f Ejemplo de las observaciones realizadas en la Estación de Ica. Figura 8. Ejemplo de uno de las tomas de datos en ICA COMO LLENAR EL CUADRO DE LA IZQUIERDA (Figura 8) Primero.- Introducir los datos como fecha, hora en tiempo universal al programa MICA para obtener el valor de P, B0, L0. Segundo.- Introducir los datos como fecha, hora en tiempo universal al programa MICA para obtener el valor de P, B0, L0, como lo describimos anteriormente. Tercero.-Trazar una recta formando un ángulo P con el eje terrestre, para trazar el eje solar. Cuarto.-Empezar a enumerar cada grupo de mancha, los del sur de color rojo y los A. Quispe 144 del norte de color azul y llenar la columna donde señala (Ica Nº). Quinto.-Contabilizar cuantas manchas hay en cada grupo y llenar la columna donde señala (f). Sexto.-Determinar el tipo de mancha según la clasificación de Zurich. Sétimo.-Se coloca la plantilla (Figura 5) sobre el dibujo y se determina su ubicación latitudinal S-N y se va llenando la columna donde señala (φ). Octavo.-Sobre la plantilla determinaremos su longitud relativa E-W y luego llenar la columna donde señala (λ). Noveno.-Para llenar la columna L realizamos dos pequeñas ecuaciones donde: L = L0 – λ, se resta si λ esta en el E. L = L0 + λ, se suma si λ esta en el W. CONCLUSIÓN Podemos concluir que las observaciones realizadas en la estación solar UNICA-IGP por el método de proyección solar nos ha permitido realizar un buen registro de manchas solares, es un método muy práctico de realizarlo. La secuencia diaria de estos registros resulta fundamental para determinar la trayectoria de las manchas solares. Sin embargo estos datos pueden contener errores que pueden ser por muchos factores, entre estos tenemos: • Debido a la calidad del seguimiento. • A la transparencia que el cielo puede presentar al momento de realizar las observaciones. • Debido a las turbulencias que presentan como consecuencia del movimiento de la atmósfera. • Debido a la estabilidad del tablero, donde se dibuja las manchas solares. • A la experiencia del observador que le permite describir correctamente los grupos de manchas solares. AGRADECIMIENTOS Agradecimiento al Director de la Estación Telemétrica y Satelital de Ancón perteneciente al Instituto Geofísico del Perú Dr. Mutsumi Ishitsuka Komaki por permitirme realizar mis prácticas Pre- Profesionales en el Área de Astronomía y por sus enseñanzas impartidas durante el desarrollo del presente informe. Al profesor Hugo Trigoso Aviles por su apoyo y enseñanzas. BIBLIOGRAFIA García de la Rosa, (1989). Apuntes de física solar, Instituto Astrofísico de Canarias, España. Descripción del método de proyección para el conteo de manchas solares en la Estación Solar única-IGP 145 Kaufmann; W. (1994). El Universo; Cuarta Edición. Priest, E. (1982). Solar Magneto Hydrodynamics. D. Reídle Publishing Company. Paginas Web consultadas: Página de Heliofísica de la Agrupación Astronómica de Cáceres (estadísticas mensuales) Parhelio (Registro de Actividad e Información Heliofísica de AstroRed). Home Page de SOHO (Solar and Heliospheric Observatory) Home Page de TRACE (Transition Region and Coronal Explorer) http://www.spaceweather.com/ http://ciencia.nasa.gov A. Quispe 146 Compendio de Trabajos de Investigación CNDG – Biblioteca Instituto Geofísico del Perú. V. 6 (2004) p. 147 - 156 PROCESAMIENTO Y COMPARACION DE LOS DATOS DE MANCHAS SOLARES REGISTRADOS EN EL OBSERVATORIO DE LA UNIVERSIDAD NACIONAL DE ICA CECILIA LÓPEZ CÓRDOVA Facultad de Ciencias Naturales y Matemáticas Universidad Nacional del Callao cecifisica@yahoo.es Prácticas dirigidas por: Dr. Mutsumi Ishitsuka Komaki Observatorio de Ancón RESUMEN En el siguiente trabajo se describe el método empleado en analizar los datos de manchas solares registrados por los estudiantes de la Facultad de Física de la Universidad Nacional San Luis Gonzaga de Ica (UNICA). Se basa en el procesamiento y comparación de los mencionados datos. El procesamiento consiste en el cálculo de los parámetros observacionales como son la apariencia del disco solar al momento de su observación, que se expresa mediante el cálculo del Número Relativo de Manchas Solares (R) y el cálculo de la constante observacional (K). La comparación consiste en evaluar el grado de error de los datos de manchas solares registrados en la UNICA en relación con los datos registrados en el observatorio de Huancayo y los datos registrados por el International Relative Sunspot Numbers de Bruselas. Como resultado del presente trabajo la base de datos que obtuvimos presentan las tendencias iniciales del grupo que observa al Sol usando esta técnica que por primera vez se esta instaurando en la UNICA. INTRODUCCION La superficie visible del Sol se denomina Fotosfera cuyo aspecto es cambiante y activa, aquí se observan una serie de perturbaciones que se conocen como actividad solar, que se define como el conjunto de fenómenos que podemos observar en la fotosfera del Sol (Priest, 1982). Cuando hablamos de actividad solar nos referimos a una gran actividad magnética en el Sol, estas perturbaciones magnéticas se pueden evidenciar mediante una diversidad de procesos de naturaleza magnética una de ellas y la mas conocida es la mancha solar1. 1 Quien tuviera mayor interés en estudiar la diversidad de los procesos magnéticos que se Las Manchas Solares son el resultado de un intenso campo magnético; cuando tales campos llegan a la superficie solar, frenan el movimiento del material propiciando un descenso de la temperatura, tal diferencia de temperatura es suficiente para ver por contraste esta región como una mancha oscura (García de la Rosa, 1989; Kazunari, y Masamitsu, 2004). La región central oscura es conocida como umbra rodeada por la penumbra mas clara, cuyo diámetro es en promedio 2.5 veces la umbra. El campo magnético en el centro de la umbra es de 2000 a 3000 Gauss y puede alcanzar alrededor de 4000 Gauss. observan en la fotosfera puede leer el capítulo 1 pp 10-23 y pp 46-55.Priest, E (1982). C. López 148 La temperatura efectiva de la mancha solar en la umbra es 3700K en comparación con la temperatura de la penumbra en donde la temperatura es aproximadamente 5600K. Para medir la actividad solar se utiliza el número de Wolf o Número Relativo de Manchas Solares (R), expresada del siguiente modo: R = (10g + f) K Al construir un gráfico en el que se representen los números de Wolf en función del tiempo, se puede constatar que la actividad solar fluctúa en ciclos de 11,2 años en promedio, de manera que en cada ciclo hay un máximo en la actividad solar (si el número de manchas, fulguraciones y protuberancias es mayor) y un mínimo en la actividad solar (con casi una ausencia total de manchas solares). Las manchas solares tienen una considerable influencia en las condiciones de nuestro planeta. Por lo que el estudio y análisis de esta actividad es muy importante, y se considera cada vez más significativa. ANTECEDENTES El Sunspot Index Data Center (SIDC) es el centro mundial de datos para establecer el índice de manchas solares, comenzó en Enero de 1981 la confección de un índice de manchas solares, llamado Número Internacional de Manchas Solares (Internacional Sunspot Number) (Ri). Su director es Pierre Cugnon situado en el Observatorio Real de Bélgica en Bruselas. El Centro está también asociado estrechamente con el departamento de radioastronomía y física solar del Observatorio Real de Bélgica, Observatorios astronómicos como los de Zúrich y Greenwich (Trigoso, 2004). NÚMERO RELATIVO DE MANCHAS SOLARES Cuando nos referimos al Número Relativo de Manchas Solares, nos referimos a la famosa ecuación definida por Rudolf Wolf en 1848. ),10( gfkR += Donde: k Es el factor de corrección del observador. R Representa al número relativo de manchas solares. f Es el número total de manchas observadas. g Es el número total de grupos de manchas observados. Esta formula nos permite comparar los valores generados con la misma fórmula por otros observatorios y/o centros mundiales de datos. Procesamiento y comparación de los datos de manchas solares registrados en el Observatorio de la Universidad Nacional de Ica 149 En donde cada parámetro observacional como son R, k, f, g se calculan en base al método de proyección. Que es un método estándar que tienen por objeto realizar el dibujo y conteo de estructuras que son posibles de ver en la imagen en luz blanca del sol, esto es, la imagen de la fotosfera solar. Sin embargo estos datos pueden contener errores, que pueden ser debido a muchos factores como: • Percepción personal (k) • A la transparencia que el cielo pueda presentar al momento de realizar la observación. • Debido a las turbulencias, que se presentan como consecuencia del movimiento convectivo de la atmósfera terrestre. • A la resolución del telescopio, el mismo que podemos calcularlo mediante la siguiente fórmula: )( 206265)(22.1 cmD cmR ∗= λ DATOS Y METODOLOGIA Características de los datos Los datos empleados en el procesamiento corresponden al periodo de Noviembre del 2003 a Marzo del 2004 y .se realizaron con ayuda de un telescopio refractor de 10cm de abertura y 1050 mm. de distancia focal marca Takahashi, y un ocular, para proyectar una imagen del disco solar de 15 cm. de diámetro. Instalado en la facultad de Física de la Universidad Nacional de Ica (UNICA), por el Instituto Geofísico del Perú (IGP) bajo un convenio realizado en abril del 2001. Para el proceso de comparación los datos empleados corresponden a: • Los registrados por los estudiantes de la universidad Nacional de Ica. • Los registrados en el observatorio de Huancayo. • Los datos del Número Relativo Internacional de Manchas Solares de Bruselas. PROCESAMIENTO DE LOS DATOS DE MANCHAS SOLARES La metodología empleada se basa en el método de los mínimos cuadrados para el cálculo de la constante observacional K. Las ecuaciones básicas empleadas son: 1) yi = K xi Vi = K xi – yi 2) K[x] – [y] = 0 [ ] [ ]x yK = Se cumple que el error o residuo: mínimoVV iV ∑ ==∑ ≈∑ 2 0 C. López 150 Además deberemos calcular el grado de dispersión que presentan los datos para esto debemos calcular el factor de correlación R2 para cada uno de los observadores. Donde: [ ] = Σ Vi = Es el error (residuo) x = Es el número de Wolf del observador (RObs) y = Es el número de Wolf del Centro Mundial de Manchas Solares de Bruselas (R Brus) De la ecuación de Wolf, tenemos: ( ) f 10gK BrusR Obs R 43421 += [ ]   = ObsR BrusRK ' BrusRObsRKiV −= ' El valor de k es un indicador de la aptitud que tiene cada observador al momento de realizar las observaciones, por lo que este valor debe calcularse por separado para cada observador, y lo designamos como Kobs1, Kobs2, Kobs3. Inicialmente se considero la unidad para la constante observacional (k),y posteriormente se calculo las constantes individualmente Kobs1, Kobs2, Kobs3. con la finalidad de obtener el número de Wolf de cada observador. CÁLCULO DEL ÍNDICE DE ACTIVIDAD SOLAR (NÚMERO DE WOLF) OBSERVADOR 1 (Figura 1): Kobs1 = 015.0443.0 ± Robs1 = Kobs1 (10g + f), para 47 datos 0=∑ iV R2=0.800 0 50 100 150 200 250 0 20 40 60 80 100 120 140 Correlación de los datos del Observador1 Y=0.443X K=0.443 ±0.015 R2=0.800 N=47 Número Robs1 R b ru se la s Figura 1. Correlación entre el Número Internacional de Manchas Solares de Bruselas con el número Internacional de Ica que corresponde al observador1. Observamos que la mayoría de los puntos están cercanos a la línea de tendencia. Procesamiento y comparación de los datos de manchas solares registrados en el Observatorio de la Universidad Nacional de Ica 151 OBSERVADOR 2 (Figura2): Kobs2 = 0.019 0.428 ± , Robs2 = Kobs2 (10g + f), para 40 datos 0=∑ iV R2 = 0.782 0 50 100 150 200 250 300 0 20 40 60 80 100 120 140 Y=0.428X K=0.428 ±0.019 R2=0.782 N=40 Correlación de los datos del Observador2 R b ru se las Número Robs2 Figura 2. Correlación entre el Número Internacional de Manchas Solares de Bruselas con el número Internacional de Ica que corresponde al observador2. Observamos que la dispersión es bastante alta y se mantiene una simetría aceptable. . OBSERVADOR 3 (Figura 3): Kobs3 = 0.0165790. ± , Robs3 = Kobs3 (10g + f), para 38 datos 0=∑ iV R2 =0.905 0 50 100 150 200 250 300 0 20 40 60 80 100 120 140 Correlación de los datos del Observador3 Y=0.579X K=0.579 ±0.016 R2=0.905 N=38 R br us el as Número Robs3 Figura 3. Correlación entre el Número Internacional de Manchas Solares de Bruselas con el número Internacional de Ica para el observador 3. A pesar de observar una dispersión baja, apreciamos una asimetría muy fuerte en relación a los resultados anteriores. C. López 152 0 20 40 60 80 100 120 140 160 180 0 20 40 60 80 100 120 140 Número R Huancayo R b ru se las Correlación de los datos Huancayo Y = 0.213856+0.76 K = 0.76±0.04 R2=0.858 N = 61 Figura 4. En la siguiente figura muestra la correlación de los datos de Huancayo, podemos observar que en la dispersión de los datos el valor de la constante k para señor Armengol de la Cruz es de 0.76±0.04 y con una correlación de R2=0.86. de donde podemos apreciar una alta simetría. COMPARACION DE LOS DATOS DE MANCHAS SOLARES En la fase de comparación se evalúa el grado de error de los datos registrados, analizando la tendencia que presenta cada observador en relación con los datos del centro mundial de datos de manchas solares de Bruselas. OBSERVADOR 1 (Figura 5): 17/11/03 07/12/03 27/12/03 16/01/04 05/02/04 25/02/04 0 20 40 60 80 100 120 140 Kobs1=0.443 FECHA NU M ER O RE LA TI VO IN TE RN AC IO NA L Bruselas Obs1 Figura 5. Comparación gráfica del Número Internacional Definitivo de Manchas Solares de Bruselas, con el número de Ica del observador1 para los meses de noviembre a diciembre del 2003 y enero a febrero del 2004. Observaremos que los datos tienen tendencias aproximadas, observándose en algunos puntos una disminución en el número de Wolf de Ica. Procesamiento y comparación de los datos de manchas solares registrados en el Observatorio de la Universidad Nacional de Ica 153 OBSERVADOR 2 (Figura 6): 17/11/03 07/12/03 27/12/03 16/01/04 05/02/04 25/02/04 0 20 40 60 80 100 120 140 Kobs2=0.428 N U M ER O R EL AT IV O IN TE R N AC IO N AL FECHA Bruselas Obs2 Figura 6. Comparación gráfica del Número Internacional Definitivo de Manchas Solares de Bruselas, con el número de Ica del observador2 para los meses de noviembre a diciembre del 2003 y enero, febrero y marzo del 2004. Se observa que la tendencia es aproximada aunque presenta mayor dispersión en relación al Obs1. OBSERVADOR 3 (Figura 7): 17/11/03 07/12/03 27/12/03 16/01/04 05/02/04 25/02/04 0 20 40 60 80 100 120 140 160 Kobs3=0.579 Bruselas Obs3 N U M ER O R EL AT IV O IN TE R N AC IO N AL FECHA Figura 7. Comparación gráfica del Número Internacional Definitivo de Manchas Solares de Bruselas, con el número de Ica observador3 para los meses de noviembre a diciembre del 2003 y enero, febrero del 2004. Notaremos que el valor de k es mayor sin embargo presenta una mayor asimetría. C. López 154 01/03/03 16/03/03 31/03/03 15/04/03 30/04/03 15/05/03 30/05/03 14/06/03 29/06/03 0 20 40 60 80 100 120 140 160 KHuancayo=0.76 FECHA N U M ER O R EL AT IV O IN TE R N AC IO N AL Bruselas Huancayo Figura 8. Comparación gráfica del Número Internacional Definitivo de Manchas Solares de Bruselas, con el número de Wolf de Huancayo que corresponde a los meses de Marzo a Junio 2003 Se puede observar que la tendencia es muy aproximada, y mientras los datos cuestionados del observatorio de Huancayo parecen interpretar una disminución del número de Wolf, en realidad hay un incremento de esta variable en los demás puntos se observa una buena correlación de los datos .por lo que podemos concluir que el observador 3 presenta un mejor índice de aproximación. Trigoso. H., (2004). CONCLUSIONES * Las observaciones registradas diariamente presentan una secuencia en la grafica, sin embargo la tendencia será discontinua al no registrarse datos en forma diaria, y los datos obtenidos no serán confiables debido a que el observador pierde la continuidad. * Los valores calculados de la constante observacional (K) que corresponde a cada observador son aceptables dentro del rango establecido. * Al realizar la comparación grafica del Número Relativo Internacional o número de Wolf, del Centro mundial de datos de mancha solares ubicado en Bruselas y el número de relativo de manchas solares de Ica, notamos que existe bastante aproximación con un mínimo de error para este periodo de tiempo. Así mismo al realizar la comparación con el Observatorio de Huancayo, este ultimo presenta una mejor aproximación debido a la continuidad de sus registros (Figuras 4 y 8). * Los resultados obtenidos por los obs1,obs2,obs3 son datos preliminares no son determinantes para realizar un análisis mas exhaustivo. Observación: La calidad de los datos registrados en el Observatorio de Huancayo supera a los registrados en Ica, esto es debido a la secuencia de sus registros. Procesamiento y comparación de los datos de manchas solares registrados en el Observatorio de la Universidad Nacional de Ica 155 AGRADECIMIENTOS Al Instituto Geofìsico del Perù, por contribuir al desarrollo de la investigaciòn cientìfica. Mi especial reconocimiento al director del Observatorio de Ancòn perteneciente al Instituto Geofìsico del Perù Dr. Mutsumi Ishitsuka Komaki por su apoyo, confianza y por brindarme la oportunidad de dessarrollarme profesionalmente. Al profesor Hugo Trigoso Aviles por su constante apoyo. A mis padres porque son la razon de mi esfuerzo. BIBLIOGRAFIA Garcia de la Rosa, I. (1989): Apuntes de Fisica Solar; Instituto Astrofísico de Canarias. Kazunari, S. y Masamitsu, O. (2004): The Sun. SHOKABO Pub.Co.Ltd., Japan. Priest, E. (1982): Solar Magnetophydrodynamics, D. Reídle Publishing Company, England. Trigoso, H. (2004): Informe Técnico Análisis de los datos de Observaciones de Manchas Solares, en el observatorio de Huancayo, IGP Páginas WEB consultadas SIDC (Sunspot Index Data Center) http://www.oma.be/KSB-ORB/SIDC/ Solar Data Analysis Center http://umbra.nascom.nasa.gov/sdac.html Página de Heliofísica de la Agrupación Astronómica de Cáceres (estadísticas mensuales) http://personal.redestb.es/intercom/pages/s ol_est.htm Página del Sol del departamento de Astrofísica de la Universidad Complutense de Madrid http://www.ucm.es/OTROS/Astrof/sol.htm l Parhelio (Registro de Actividad e Información Heliofísica de AstroRed) http://www.astrored.net/helio. Home Page de SOHO (Solar and Heliospheric Observatory) http://sohowww.nascom.nasa.gov/ Home Page de TRACE (Transition Region and Coronal Explorer) http://vestige.lmsal.com/TRACE/ Noticias e información sobre el estado del Sol http://www.spaceweather.com/ C. López 156