Leer · Uso de técnicas de optimización global, para resolver problemas de inversión sísmica, segundo grupo de exploración

Página 1 de 169Uso de técnicas de optimización global, para resolver problemas de inv…

p. 1

USO DE T ´ECNICAS DE OPTIMIZACI ´ON GLOBAL PARA RESOLVER

PROBLEMAS DE INVERSI ´ON S´ISMICA - SEGUNDO GRUPO DE

EXPLORACI ´ON

JAIR ABRAHAM BELE ˜NO DAVILA

LIBARDO ANDR ´ES ESCALANTE ROJAS

YUDY ANDREA SARMIENTO L ´OPEZ

UNIVERSIDAD INDUSTRIAL DE SANTANDER

FACULTAD DE INGENIER´IAS F´ISICO - MEC ´ANICAS

ESCUELA DE INGENIER´IAS EL ´ECTRICA, ELECTR ´ONICA Y

TELECOMUNICACIONES

BUCARAMANGA

2015

p. 2

USO DE T ´ECNICAS DE OPTIMIZACI ´ON GLOBAL PARA RESOLVER

PROBLEMAS DE INVERSI ´ON S´ISMICA - SEGUNDO GRUPO DE

EXPLORACI ´ON

JAIR ABRAHAM BELE ˜NO DAVILA

LIBARDO ANDR ´ES ESCALANTE ROJAS

YUDY ANDREA SARMIENTO L ´OPEZ

TRABAJO DE GRADO PARA OPTAR AL T´ITULO DE INGENIERO

ELECTR ´ONICO

DIRECTORES:

Dr.Ing. OSCAR MAURICIO REYES TORRES Ph.D(c) SERGIO ALBERTO ABREO CARRILLO

UNIVERSIDAD INDUSTRIAL DE SANTANDER

FACULTAD DE INGENIER´IAS F´ISICO - MEC ´ANICAS

ESCUELA DE INGENIER´IAS EL ´ECTRICA, ELECTR ´ONICA Y

TELECOMUNICACIONES

BUCARAMANGA

2015

p. 3

DEDICATORIA

Dedico este trabajo de grado a Dios y a mis padres. A Dios porque ha estado conmigo a cada paso que doy, cuid´andome y d´andome fortaleza para continuar, a mis padres, quienes a lo largo de mi vida han velado por mi bienestar y educaci´on siendo mi apoyo en todo momento. Es por ellos que soy lo que soy ahora. Los amo con mi vida.

Andrea

p. 4

DEDICATORIA

Dedico este trabajo de grado con todo mi amor y mi cari˜no para las personas que hicieron todo en la vida para que yo pudiera lograr y alcanzar este sue˜no, por motivarme, por su paciencia y comprensi´on, sacriﬁcaron su tiempo para que yo pudiera cumplir con el m´ıo, por entregarme su apoyo cuando sent´ıa que el camino se terminaba.

Mi eterna gratitud a esas personas importantes en mi vida, que siempre estuvieron listas para brindarme toda su ayuda, ahora me toca regresar un poco de todo lo inmenso que me han otorgado. A mi querida madre Lucila, quien me ha entregado el regalo de la vida y su incansable amor que me hace sentir el hijo m´as afortunado. A mis bellos hermanos, Diana y Camilo, quienes representan mis pilares en la vida, son mi mayor alegr´ıa y siempre estaremos juntos. A mi pr´ıncipe de ojos azules, Santiago, quien desde el d´ıa de su nacimiento me entrego la luz de su sonrisa. A el amor de mi vida, Jainy, quien llor´o y ri´o en cada momento junto a m´ı y fue capaz de contenerme cuando todo iba mal. Gracias por amarme como solo t´u lo puedes hacer, por entregarme esa hermosa sonrisa con la que iluminas todos y cada uno de mis d´ıas, por hacer de mi el hombre m´as feliz y por vivir este amor tan hermoso. Son mi familia, son mi raz´on de vivir, son los que hacen que mi vida tenga alg´un signiﬁcado y me brindan su amor y eso es lo m´as valioso para mi. Por y para siempre vivir´an en mi coraz´on. Los amo demasiado. Cinco sonrisas. Cinco corazones.

“Antes de que el mundo fuera creado, ya exist´ıa la palabra. La palabra era la fuente de vida, y esta vida trajo la luz a la humanidad. La luz brilla en la oscuridad, y la oscuridad nunca ha sido apagada”. - Hiroya Oku.

Jair

p. 5

DEDICATORIA

Al buen Dios de amor que cada d´ıa no encontr´o qu´e excusas inventarse para llamar mi atenci´on, para conquistarme y mostrar su respaldo y amor. Para ´El, que nos habilit´o con salud, sabidur´ıa y paciencia para afrontar cada reto. Por su sostenimiento y cuidado cuando se decidi´o partir de nuestros hogares para afrontar esta nueva etapa de la vida. Porque nunca ´El fue escaso y dej´o que siempre disfrut´aramos de su compa˜n´ıa. Por dejarme cada d´ıa respirar y parpadear... ¡Inﬁnitas gracias Viejo!

A mi Pap´a Libardo Escalante Quintero que me ense˜naba a nunca bajar la guardia y me recordaba que ten´ıa muchas capacidades para afrontar cada situaci´on. A mi mam´a Maria Delia Rojas, que cada d´ıa me recordaba .Estar de rodillas delante de Dios para estar de pie delante de los hombres”. Para ellos quienes nunca dejaron de estar al pie del ca˜n´on, no teniendo en cuenta el tiempo ni la distancia, ni los recursos, ni aun cuando pensamos en desfallecer, ellos estuvieron ah´ı esperanzados en que la recompensa ser´ıa muy grande. Por amarme, apoyarme y orar cada d´ıa por mi vida, gracias.

A mis abuelitas, Elodia Quintero y Etelvina D´ıaz, quienes a pesar de sus limitaciones y la distancia, nunca dejaron de preguntar con ternura “¿Ya casi termina?”, buscando siempre as´ı una respuesta alentadora y esperando volver a disfrutar tiempo juntos.

A la familia Esparza Escalante, Mart´ın, mi hermana Mayra, sus hijas, Salito y Luci, quienes me acogieron en su hogar no midiendo nada, sino por el contrario, gracias por aportar esa cuota de compa˜n´ıa y amor. De igual forma a las familias Escalante, Rojas, Rangel, Borja, Giraldo, gracias por su apoyo y por sus recursos.

p. 6

Al movimiento estudiantil y profesional Alfa y Omega, y la familia del CENTI, representados en l´ıderes, disc´ıpulos y hermanos en la fe, que me acogieron como un integrante m´as permiti´endome compartir un v´ınculo Eterno. Por sus oraciones y compa˜n´ıa, gracias.

A mis amigos y hermanos, gracias por su apoyo, por hacer que su compa˜n´ıa y cada encuentro, fuera un factor importante para recobrar fuerzas. A mis profesores, grupo de investigaci´on CPS, director y co-director de proyecto, gracias por sus aportes y por dar lo mejor en esa hermosa labor. “Ser´a como ´arbol plantado junto a corrientes de aguas, Que da su fruto en su tiempo, Y su hoja no cae; Y todo lo que hace, prosperar´a”. Salmo 1:3. Para ti, Gracias.

“Una sola vida tenemos y esta hay que invertirla”. - N´estor Chamorro Pesantes. Libardo

p. 7

AGRADECIMIENTOS

Este trabajo fue realizado con la dedicaci´on y la ilusi´on de lograr un objetivo, de llevar a cabo un sue˜no y hacerlo realidad.

Agradecemos a Dios por habernos guiado por el camino de la felicidad hasta ahora. Tambi´en, a cada uno de los que son parte de nuestras familias por su fuerza y apoyo incondicional que nos han ayudado y llevado hasta donde estamos ahora. Por ´ultimo a nuestros compa˜neros de tesis porque en esta armon´ıa grupal lo hemos logrado.

Y agradecemos a los profesores Oscar Reyes Y Sergio Abreo por apoyarnos y entregarnos su conﬁanza en el transcurso del trabajo.

Andrea, Jair y Libardo

p. 8

CONTENIDO

INTRODUCCI ´ON . . . . . . . . . . . . . . . . . . . . . . . . . .

DESCRIPCI ´ON Y PLANTEAMIENTO DEL SEMINARIO

. . . . . .

1.1

INTRODUCCI ´ON . . . . . . . . . . . . . . . . . . . . . . .

1.2

DESCRIPCI ´ON DEL SEMINARIO. . . . . . . . . . . . . . . .

1.3

MODELO MATEM ´ATICO EMPLEADO. . . . . . . . . . . . . .

1.4

PLANTEAMIENTO DEL PROBLEMA . . . . . . . . . . . . . .

MODELADO Y ADQUISICI ´ON

. . . . . . . . . . . . . . . . .

2.1

INTRODUCCI ´ON . . . . . . . . . . . . . . . . . . . . . . .

2.2

CONSTRUCCI ´ON DE MODELOS GEOF´ISICOS EN SU . . . . . .

2.3

ADQUISICI ´ON S´ISMICA EN SU . . . . . . . . . . . . . . . .

2.3.1

Trazado de rayos . . . . . . . . . . . . . . . . . . . . . . .

2.3.2

Adquisici´on y trazas s´ısmicas. . . . . . . . . . . . . . . . . .

M ´ETODOS DE SOLUCI ´ON METAHEUR´ISTICOS . . . . . . . . .

3.1

INTRODUCCI ´ON . . . . . . . . . . . . . . . . . . . . . . .

3.2

M ´ETODOS DE SOLUCI ´ON

. . . . . . . . . . . . . . . . . .

p. 9

3.3

SIMULATED ANNEALING . . . . . . . . . . . . . . . . . . .

3.3.1

Descripci´on de par´ametros . . . . . . . . . . . . . . . . . . .

3.3.2

Diagrama de ﬂujo de SA . . . . . . . . . . . . . . . . . . . .

3.3.3

Descripci´on e implementaci´on de las funciones de prueba en SA . .

3.4

ARTIFICIAL BEE COLONY

. . . . . . . . . . . . . . . . . .

3.4.1

Descripci´on de par´ametros . . . . . . . . . . . . . . . . . . .

3.4.2

Diagrama de ﬂujo de ABC . . . . . . . . . . . . . . . . . . .

3.4.3

Descripci´on e implementaci´on de las funciones de prueba en ABC .

3.5

INTERFAZ ENTRE MATLAB Y SEISMIC UNIX . . . . . . . . . .

3.5.1

SA y ABC implementados en el problema de inversi´on s´ısmica . . .

PRUEBAS Y RESULTADOS. . . . . . . . . . . . . . . . . . .

4.1

METODOLOG´ıA . . . . . . . . . . . . . . . . . . . . . . .

4.2

RESULTADOS . . . . . . . . . . . . . . . . . . . . . . . .

CONCLUSIONES Y RECOMENDACIONES. . . . . . . . . . . .

5.1

CONCLUSIONES. . . . . . . . . . . . . . . . . . . . . . .

5.2

RECOMENDACIONES

. . . . . . . . . . . . . . . . . . . .

REFERENCIAS BIBLIOGR ´AFICAS . . . . . . . . . . . . . . . . . .

p. 10

BIBLIOGRAF´IA . . . . . . . . . . . . . . . . . . . . . . . . . . .

ANEXOS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

p. 11

LISTA DE FIGURAS

Descripci´on de los temas abordados a nivel geof´ısico. . . . . . . . . .

Modelo del subsuelo . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Conjunto de trazas s´ısmicas

. . . . . . . . . . . . . . . . . . . . . . .

Descripci´on del seminario de investigaci´on

. . . . . . . . . . . . . . .

Primer Modelo

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Modelos de dos capas paralelas . . . . . . . . . . . . . . . . . . . . .

Modelos de cuatro capas paralelas . . . . . . . . . . . . . . . . . . . .

Modelos de cuatro capas no uniformes

. . . . . . . . . . . . . . . . .

Modelos de cuatro capas no uniformes con una variante . . . . . . . .

Trazado de rayos sobre el modelo de capas paralelas

. . . . . . . . .

Trazado de rayos sobre el modelo 3, modelo de 4 capas del subsuelo

Trazado de rayos sobre el modelo 4, modelo de 4 capas con intrusi´on

Representaciones de las trazas s´ısmicas

. . . . . . . . . . . . . . . .

Descripci´on de los m´etodos de soluci´on . . . . . . . . . . . . . . . . .

Diagrama de ﬂujo de SA . . . . . . . . . . . . . . . . . . . . . . . . . .

Funciones prueba en SA

. . . . . . . . . . . . . . . . . . . . . . . . .

p. 12

Diagrama de ﬂujo de ABC . . . . . . . . . . . . . . . . . . . . . . . . .

Diagrama de comunicaci´on entre MATLAB y SU

. . . . . . . . . . . .

Promedio de los coeﬁcientes de la correlaci´on cruzada con respecto al n´umero de iteraciones en la primera prueba con SA.

. . . . . . . .

Promedio de los coeﬁcientes de la correlaci´on cruzada con respecto al n´umero de ge´ofonos en la primera prueba con SA.

. . . . . . . . .

Promedio de sloths en la segunda prueba con SA. . . . . . . . . . . .

Promedio de los coeﬁcientes de la correlaci´on cruzada con respecto al n´umero de ge´ofonos en la segunda prueba con SA. . . . . . . . . .

p. 13

LISTA DE TABLAS

Velocidades s´ısmicas

. . . . . . . . . . . . . . . . . . . . . . . . . . .

Valores de sloth para cada una de las capas del tercer modelo . . . .

Tabla de funciones prueba . . . . . . . . . . . . . . . . . . . . . . . . .

Especiﬁcaciones del equipo de c´omputo usado en los experimentos sobre las funciones de prueba

. . . . . . . . . . . . . . . . . . . . . .

Par´ametros bajo los cuales se obtuvieron los mejores resultados en la implementaci´on de SA en el conjunto de funciones de prueba

. . .

Promedios y desviaci´on de las funciones prueba de SA . . . . . . . .

Par´ametros bajo los cuales se obtuvieron los mejores resultados en la implementaci´on de ABC en el conjunto de funciones de prueba

. .

Promedios y desviaci´on de las funciones prueba de ABC

. . . . . . .

Especiﬁcaciones del equipo de c´omputo usado en los experimentos sobre el proceso de inversi´on s´ısmica . . . . . . . . . . . . . . . . . .

Par´ametros de la primera prueba con SA

. . . . . . . . . . . . . . . .

Alimento obtenido por la colonia en la prueba realizada con ABC . . .

Soluciones exploradas por la colonia en la prueba realizada con ABC

p. 14

Comparaci´on de sloths ﬁnales entre SA y ABC con respecto al modelo de referencia . . . . . . . . . . . . . . . . . . . . . . . . . . . .

p. 15

RESUMEN

TITULO:

USO DE T ´ECNICAS DE OPTIMIZACI ´ON GLOBAL PARA RESOLVER

PROBLEMAS DE INVERSI ´ON S´ISMICA - SEGUNDO GRUPO DE

EXPLORACI ´ON *

AUTORES:

JAIR ABRAHAM BELE ˜NO DAVILA

LIBARDO ANDR ´ES ESCALANTE ROJAS

YUDI ANDREA SARMIENTO L ´OPEZ**

PALABRAS CLAVES:

Implementaci´on, Metaheur´ıstica, Modelo, Optimizaci´on, S´ısmica, T´ecnica.

DESCRIPCI ´ON:

Este libro contiene el desarrollo del trabajo de grado bajo la metodolog´ıa de seminario de investigaci´on llamado Uso de t´ecnicas de optimizaci´on global para resolver problemas de inversi´on s´ısmica, donde el objetivo principal consiste en el estudio, implementaci´on y an´alisis de t´ecnicas metaheur´ısticas, aplicadas a problemas de inversi´on s´ısmica.

Como inicio hacia la comprensi´on del problema propuesto, se investiga la deﬁnici´on de conceptos b´asicos en geof´ısica, en inversi´on s´ısmica, en optimizaci´on global y sus t´ecnicas metaheur´ısticas. Adicionalmente se presenta el modelo matem´atico y sus par´ametros. Asimismo, se presenta la caracterizaci´on e implementaci´on de dos t´ecnicas metaheur´ısticas denominadas como Optimizaci´on mediante Recocido Simulado y Colonia Artiﬁcial de Abejas, las cuales fueron probadas sobre funciones matem´aticas de prueba. Una vez el funcionamiento de las t´ecnicas son ajustadas, se deﬁne una m´etrica de comparaci´on para aplicar las t´ecnicas metaheur´ısticas adecuando los datos modelados con respecto a los datos observados. Los resultados de mayor importancia en la implementaci´on de las t´ecnicas metaheur´ısticas sobre el problema propuesto, se registran mediante tablas, gr´aﬁcas y bases de datos digitales. Se presenta un an´alisis, discusiones, recomendaciones y conclusiones. Se generan unos documentos adicionales conformados por actas y memorias que describen los temas tratados en cada sesi´on del seminario y anexos que contienen el desarrollo de los c´odigos implementados. Estos documentos y archivos constituyen un conjunto de referencias para la realizaci´on de trabajos futuros en esta ´area.

*Trabajo de Grado en la Modalidad de Seminario de Investigaci´on. **Facultad de Ingenier´ıas F´ısico-Mec´anicas. Escuela de Ingenier´ıas El´ectrica, Electr´onica y Telecomunicaciones. Director Dr.Ing. Oscar Mauricio Reyes Torres. Codirector PhD(c). Sergio Alberto Abreo Carrillo.

p. 16

ABSTRACT

TITTLE:

USING GLOBAL OPTIMIZATION TECHNIQUES TO SOLVE SEISMIC

INVERSION PROBLEMS - SECOND GROUP OF EXPLORATION*

AUTHORS:

JAIR ABRAHAM BELE ˜NO DAVILA

LIBARDO ANDR ´ES ESCALANTE ROJAS

YUDI ANDREA SARMIENTO L ´OPEZ**

KEYWORDS:

Implementation, Metaheuristic, Model, Optimization, Seismic, Technical.

This book contains the development of the project work under the methodology of research seminar called Using global optimization techniques for solving seismic inversion, where the main objective is the study, implementation and analysis of metaheuristic techniques, applied to seismic inversion problems. To understand the proposed problem, the deﬁnition of basic concepts in geophysics, seismic inversion, global optimization and metaheuristic techniques are investigated. Additionally mathematical models and their parameters are presented. Also, the characterization and implementation of two metaheuristic optimization techniques called Simulated Annealing and Artiﬁcial Bee Colony were tested on math test functions. Once operating techniques are adjust, a metric of comparison is deﬁned to implement the metaheuristic techniques modeled to ﬁt the data with respect to the observed data. The results of major importance in the implementation of metaheuristic techniques on the proposed problem were register using tables, graphs and digital databases. Analysis, discussions, recommendations and conclusions are presented. An additional document is provided by made up of records and reports describing the topics covered in each seminar session and annexes containing the development of codes generated implemented. These documents and ﬁles are a set of references for conducting future work in this area.

*Project Work in Seminar Research Mode.

**Faculty of Physicomechanical Engineers. School of Electrical Engineers, Electronics and Telecommunications. Director Dr.Ing. Oscar Mauricio Reyes Torres. Codirector PhD(c). Sergio Alberto Abreo Carrillo.

p. 17

INTRODUCCI ´ON

El contenido de este documento presenta el desarrollo del trabajo de grado bajo la modalidad de seminario de investigaci´on llamado Uso de t´ecnicas de optimizaci´on global para resolver problemas de inversi´on s´ısmica, propuesto con el objetivo de estudiar, analizar e implementar algoritmos de optimizaci´on global basados en estas t´ecnicas metaheur´ısticas para ser aplicados a problemas de inversi´on s´ısmica. Este documento est´a constituido por cinco cap´ıtulos en donde se presenta el trabajo realizado para solucionar la problem´atica identiﬁcada; la cual corresponde a la incorporaci´on de dos algoritmos de optimizaci´on global para resolver el problema de inversi´on s´ısmica. Durante el desarrollo del seminario de investigaci´on se deﬁnen conceptos necesarios para el entendimiento de la geof´ısica para identiﬁcar las caracter´ısticas del subsuelo. Adem´as, de poder revisar t´ecnicas para la comparaci´on de datos s´ısmicos, implementar t´ecnicas metaheur´ısticas de optimizaci´on global y aplicar el procesamiento de se˜nales para la comparaci´on de las se˜nales digitales para obtener la estructura litol´ogica de una porci´on de corteza terrestre. Con base en esto se realizan conclusiones y recomendaciones que se complementan en los anexos para presentar en detalle los temas expuestos por medio de c´odigos implementados y bases de datos conformadas por actas y memorias.

El primer cap´ıtulo expone los conceptos b´asicos en geof´ısica, inversi´on s´ısmica, optimizaci´on y t´ecnicas metaheur´ısticas, para poder comprender el planteamiento del problema que se discute posteriormente, esto con el ﬁn de aclarar el contexto bajo el cual se desarrolla este seminario de investigaci´on.

p. 18

En el segundo cap´ıtulo se presenta a Seismic Unix como una de las herramientas computacionales necesarias para llevar a cabo el seminario de investigaci´on, sabiendo que ´esta contribuye de manera signiﬁcativa en las simulaciones enfocadas al ´area de la geof´ısica; teniendo como objetivo principal el desarrollo de modelos de velocidades de una secci´on del subsuelo y la adquisici´on s´ısmica a implementar sobre el mismo.

En el tercer cap´ıtulo se revisan generalidades y conceptos b´asicos de las t´ecnicas metaheur´ısticas seleccionadas, las cuales son Simulated Annealing (SA) y Artiﬁcial Bee Colony (ABC). Cada una de ellas se centra en aspectos como: la analog´ıa que tiene cada una con la naturaleza, la parametrizaci´on, su respectiva evaluaci´on sobre funciones de prueba y su empalme con el proceso de inversi´on s´ısmica. El cuarto capitulo esta conformado por las simulaciones y pruebas realizadas en el transcurso del seminario de investigaci´on con SA y ABC en el proceso de inversi´on s´ısmica. El quinto y ´ultimo cap´ıtulo comprende las conclusiones y recomendaciones relacionadas con los resultados obtenidos en la implemantaci´on de las t´ecnicas metaheur´ısticas en el problema de inversi´on s´ısmica. Finalmente se anexan de forma digital los documentos que evidencian el proceso del seminario de investigaci´on que hacen referencia a actas semanales en donde se explican brevemente los temas tratados, las discusiones generadas, los compromisos correspondientes a las pr´oximas sesiones, los asistentes, expositores y relatores en cada sesi´on. Por otra parte, se encuentran las memorias que presentan los temas tratados durante cada sesi´on, junto con el an´alisis y desarrollo de los avances presentados. La estructura de cada memoria consta de una introducci´on, temas tratados, implementaci´on, an´alisis de resultados, discusi´on y conclusiones.

p. 19

Como documentos auxiliares se encuentran los ap´endices, que contienen conceptos b´asicos, la secuencia de los c´odigos implementados, y pruebas realizadas con las implementaciones para cada algoritmo. Dichos ap´endices se presentan como una extensi´on de las memorias y a su vez como informaci´on complementaria del documento principal.

La metodolog´ıa llevada a cabo a lo largo del seminario de investigaci´on consisti´o en una sesi´on semanal durante el semestre con duraci´on de dos horas cada una. En cada una de ellas se presentaron las bases te´oricas e implementaciones que se iban requiriendo a medida que se avanzaba en la soluci´on del problema conforme al cronograma estipulado para el desarrollo del seminario de investigaci´on. Cada sesi´on cont´o con espacios de presentaci´on de la tem´atica a tratar acompa˜nada de discusi´on, an´alisis y ajustes en relaci´on con el contenido presentado. De igual manera, cont´o con los aportes propuestos por los directores y los asistentes a dichas sesiones.

p. 20

DESCRIPCI ´ON Y PLANTEAMIENTO DEL SEMINARIO

1.1

INTRODUCCI ´ON

Este cap´ıtulo contiene los principios necesarios (con base en una revisi´on bibliogr´aﬁca de conceptos y deﬁniciones b´asicas) que permiten entender el problema discutido en el seminario de investigaci´on. La descripci´on se hace de una manera gr´aﬁca para facilitarle al lector la comprensi´on de los conceptos en geof´ısica y la forma en que se trabaj´o durante el seminario. La primera parte de la revisi´on se hace de manera secuencial, donde se presentan algunas deﬁniciones de las ciencias correspondientes al tema de estudio (Figura 1) y a continuaci´on se expone el modelo matem´atico que se utiliza durante todo el seminario de investigaci´on. Finalmente se presenta el planteamiento del problema para comprender de manera general el trabajo que se realiz´o.

p. 21

1.2

DESCRIPCI ´ON DEL SEMINARIO

Figura 1. Descripci´on de los temas abordados a nivel geof´ısico. La geof´ısica estudia los fen´omenos f´ısicos de la Tierra, donde se comprenden dos ramas principales, las cuales son la geof´ısica interna y la geof´ısica externa. La segunda estudia las propiedades f´ısicas del entorno terrestre. La primera es necesaria para el an´alisis del interior de la corteza terrestre. Una de sus ramas, la sismolog´ıa, es aquella que comprende el estudio de la propagaci´on de las ondas s´ısmicas naturales que viajan a trav´es del subsuelo. Con la sismolog´ıa viene anidado un proceso de generaci´on de ondas artiﬁciales la cual es denominada s´ısmica, donde su teor´ıa abarca gran cantidad de conceptos; entre ellos se encuentra la ley de la elasticidad, impedancia el´astica e impedancia ac´ustica, velocidades s´ısmicas (ondas P, ondas S, ondas de Rayleight y ondas Love) y ´angulos de incidencia, entre

p. 22

otros. El objetivo principal de la s´ısmica es la localizaci´on de yacimientos de hidrocarburos para la industria del petr´oleo. Sus principales m´etodos empleados son la s´ısmica de reﬂexi´on y s´ısmica de refracci´on; m´etodos de gran importancia para la exploraci´on del subsuelo en el descubrimiento de estos yacimientos (Anexo A). La exploraci´on s´ısmica es el m´etodo geof´ısico que permite realizar descripciones cuantitativas del subsuelo a trav´es de la generaci´on de ondas s´ısmicas artiﬁciales usando alg´un tipo de agente externo (fuente). La exploraci´on s´ısmica est´a dividida en cuatro etapas que son: la etapa topogr´aﬁca que es la encargada del estudio para la ubicaci´on de los sensores (ge´ofonos) y fuentes (martillo, vibro, etc.); la etapa de perforaci´on, se encarga de crear peque˜nos pozos para la ubicaci´on de los sensores y las fuentes; la etapa de registro en donde se lleva a cabo la recolecci´on de la informaci´on por medio de los sensores y ﬁnalmente se efect´ua la etapa de restauraci´on de la zona en donde el proceso de exploraci´on fue desarrollada. El prop´osito de la exploraci´on s´ısmica es el registro de datos e informaci´on del subsuelo. A partir de ellos se puede conseguir un modelo de impedancias el cual es el objetivo principal de la inversi´on s´ısmica, donde un problema puede ser planteado en t´erminos de hallar una funci´on no lineal multidimensional o funci´on objetivo. El prop´osito de la inversi´on s´ısmica es la transformaci´on de observaciones s´ısmicas en propiedades cuantitativas de rocas que describan un reservorio y/o yacimiento. Por lo tanto, el principio de la inversi´on s´ısmica se basa en la resoluci´on de un problema de estimaci´on de par´ametros, sabiendo que ´estos dan una serie de caracter´ısticas (por ejemplo, la elasticidad de los materiales que conforman cierta porci´on de tierra) con las cuales se modelan las capas del subsuelo (Figura 2).

p. 23

Figura 2. Modelo de capas del subsuelo.

Fuente: Editora Digital on Mirror. ¿Qu´e es el fracking y c´omo afecta al medio ambiente? Disponble en: http://www.ecoosfera.com/2013/12/que-es-el-fracking-ycomo-afecta-al-medio-ambiente Consultado el 15 de Febrero de 2015. Una traza s´ısmica es la representaci´on gr´aﬁca de la transici´on de diferentes frentes de ondas que viajan desde un emisor (fuente) hasta un receptor (ge´ofono) a trav´es del subsuelo. Por esto, para cada receptor ubicado a lo largo de una distancia horizontal existe un tiempo recepci´on de datos y cuando se agrupan se tiene un conjunto de trazas s´ısmicas (Figura 3). Este conjunto puede interpretarse como, en qu´e instante de tiempo y cu´an fuerte (amplitud) los receptores perciben perturbaciones por parte del frente de onda que regresa a la superﬁcie.

p. 24

Figura 3. Conjunto de trazas s´ısmicas.

La inversi´on s´ısmica se divide en dos procesos fundamentales, la inversi´on basada en la traza y la inversi´on basada en el modelo matem´atico; ambas tienen un objetivo ﬁnal com´un que es la estimaci´on de un modelo de capas del subsuelo que pueda simular los datos con el menor error posible en comparaci´on a un modelo geol´ogico.

El estudio de la geof´ısica y en especial el estudio de la sismolog´ıa, puede realizarse mediante una interpretaci´on visual por medio de herramientas computacionales como es el caso de Seismic Unix (SU), dise˜nado por el Centro de Fen´omenos Ondulatorios (CWP) de la Escuela de Minas de Colorado. SU es una herramienta que tiene un conjunto de rutinas de c´omputo cient´ıﬁco, donde se elaboran diversas tareas de ´ındole geof´ısico, entre ellas, el modelado de la propagaci´on de ondas s´ısmicas, adquisiciones s´ısmicas y el procesamiento de datos s´ısmicos. La plata-

p. 25

forma de funcionamiento usa entornos UNIX, donde su programaci´on es basada en shell-scripts los cuales brindan una reducci´on en el tiempo de c´omputo ya que se ejecutan rutinas que no requieren entornos gr´aﬁcos.

Para la generaci´on de modelos del subsuelo, SU utiliza la triangulaci´on de Delaunay, la cual permite representar de manera acertada estructuras simples o complejas del subsuelo. Este m´etodo es altamente utilizado para la generaci´on de gr´aﬁcos e im´agenes digitales. Asimismo, su mayor ventaja es la simplicidad en el c´alculo de tiempos de viaje y caminos de las ondas s´ısmicas. Para m´as informaci´on a lo relacionado con SU, consultar el anexo B.

En la Figura 4 se presenta el resumen de las etapas del seminario de investigaci´on. El cual se desarroll´o sobre dos columnas vertebrales, representadas en un componente geof´ısico y uno computacional. El primero, corresponde a las implementaciones sobre el software Seismic Unix. Sabiendo que previo a esto, se abarcan conceptos b´asicos en geof´ısica (orientado a los procesos de exploraci´on s´ısmica) e inversi´on s´ısmica; estas integraciones presentes en el componente geof´ısico, est´an descritas con mayor detalle en el anexo 5.2. El segundo componente hace referencia al uso de las t´ecnicas de optimizaci´on global, las cuales se evidencian en el cap´ıtulo 3 de este documento.

p. 26

Figura 4. Descripci´on del seminario de investigaci´on.

1.3

MODELO MATEM ´ATICO EMPLEADO

El modelo matem´atico es representado a trav´es de la ecuaci´on: (∂xΘ)2 + (∂yΘ)2 + (∂zΘ)2 = 1 C2

(1)

que corresponde a la aproximaci´on no lineal de la ecuaci´on Eikonal, Θ representa el frente de onda y C es la velocidad de propagaci´on de los frentes de onda en las capas del subsuelo. Esta ecuaci´on es sobre la cual se desarrolla el proceso de

p. 27

la adquisici´on s´ısmica [1].

Es necesario conocer el papel de la adquisici´on s´ısmica en el proceso de inversi´on (v´ease Figura 5), en este documento se resume en cuatro etapas. La primera corresponde a los datos observados, que son el punto de partida durante un proceso de inversi´on y ser´an los “datos de referencia” a comparar (recuadro verde). Debido a que no se poseen estos datos de referencia, se hace necesario implementar un modelado del subsuelo y realizar un primer proceso de adquisici´on s´ısmica (recuadro rojo), para as´ı obtener los datos de referencia. Figura 5. Descripci´on del proceso de inversi´on s´ısmica. Para la segunda etapa (recuadro fucsia), se parte de un modelo generado, que es lo aconsejado como gu´ıa para garantizar que las velocidades iniciales sean n´umeros positivos y mayores a cero; posteriormente se realizan un proceso de adquisici´on y obtenci´on de datos modelados. Sabiendo que los mencionados en la etapa anterior pertenecen a la generaci´on de los datos de referencia.

p. 28

La tercera etapa (recuadro azul) est´a marcada por la comparaci´on entre los datos observados y los datos modelados a trav´es de la m´etrica de comparaci´on seleccionada (correlaci´on cruzada); esta fue objetivo de discusi´on durante el transcurso del seminario ya que existen diferentes formas de comparar se˜nales digitales. Esta m´etrica se caracteriza porque entrega un vector o un escalar que representa la variaci´on entre las trazas s´ısmicas de los datos observados y los modelados. Por ´ultimo (recuadro anaranjado), la t´ecnica metaheur´ıstica necesita del resultado entregado por la m´etrica comparaci´on, para as´ı minimizar la variaci´on entre los datos observados y los modelados. Sabiendo que si existe variaci´on, es necesario que la t´ecnica evolucione para realizar una b´usqueda de nuevos valores de velocidades. Es en esta secci´on donde se aporta el trabajo de las t´ecnicas de optimizaci´on global en un proceso de inversi´on s´ısmica. De esta manera busca repetirse el ciclo y generar modelos hasta que se cumpla el criterio de parada, asumiendo que el modelo obtenido no presenta variaciones y es el m´as cercano al de referencia.

1.4

PLANTEAMIENTO DEL PROBLEMA

Uno de los principales objetivos de la inversi´on s´ısmica es encontrar modelos aproximados de las capas que conforman el subsuelo, para que as´ı m´as adelante puedan ser estudiados con la ayuda de interpretaciones de expertos. En este contexto, la inversi´on s´ısmica consiste en encontrar un modelo de velocidades minimizando la diferencia entre los datos modelados y observados a partir de la modiﬁcaci´on de los par´ametros del modelo matem´atico. El modelo geol´ogico se construye a trav´es de par´ametros f´ısicos que caracterizan las propiedades de las capas de ro-

p. 29

ca; principalmente, entre ellos pueden encontrarse la velocidad de onda P, velocidad de onda S, impedancia el´astica e impedancia ac´ustica [2]. Uno de los objetivos que se trabaj´o fue aplicar e implementar dos m´etodos de optimizaci´on global a problemas geof´ısicos. Aunque hacemos hincapi´e en los aspectos de la aplicaci´on de estos algoritmos, se presentan tambi´en los fundamentos te´oricos suﬁcientes para que los asistentes al seminario entiendan los principios subyacentes en que se basan estos algoritmos.

Otro de los objetivos es describir con suﬁciente detalle los fundamentos de dos m´etodos de optimizaci´on con aplicaci´on a la inversi´on s´ısmica, explorando y comparando diferentes alternativas de implementaci´on.

p. 30

MODELADO Y ADQUISICI ´ON

2.1

INTRODUCCI ´ON

Este cap´ıtulo contiene las etapas de aprendizaje y desarrollo con respecto al manejo del software Seismic Unix (SU) en los procesos de modelado y adquisici´on s´ısmica. Se presentan de manera precisa los pasos realizados en el proceso de creaci´on de cuatro modelos del subsuelo con caracter´ısticas diferentes. Adem´as, se da a conocer una serie de componentes necesarios para generar un proceso de adquisici´on s´ısmica sobre los modelos desarrollados en el transcurso del seminario de investigaci´on.

2.2

CONSTRUCCI ´ON DE MODELOS GEOF´ISICOS EN SU

El objetivo principal en esta etapa, es proveer al lector las herramientas necesarias para el manejo del software Seismic Unix (SU) en la creaci´on de modelos del subsuelo en 2D y en el proceso de adquisici´on s´ısmica. Con SU se crearon cuatro modelos diferentes del subsuelo, el cual uno de ellos es tomado como modelo de referencia. Los modelos construidos presentan el mismo procedimiento en la generaci´on del c´odigo fuente, pero todos son distintos. Los modelos desarrollados se consideran capas homog´eneas e isotr´opicas, cada una de las capas posee una velocidad de onda ac´ustica (Onda P), donde la informaci´on relevante entre cada una de las capas es el valor del sloth (velocidad), el cual est´a deﬁnido en la ecuaci´on 2.1

p. 31

y sus unidades son [s2/m2]. El valor de sloth es la caracter´ıstica de mayor relevancia en el modelado del subsuelo, constituye la parte fundamental en la analog´ıa de los modelos con respecto a una zona de la Tierra.

s = 1 V 2

(2)

V representa la velocidad con la que viaja una onda ac´ustica a trav´es de un medio, sus unidades son [m/s]. Distintos valores de velocidades se muestran en la Tabla 1.

Tabla 1. Velocidades s´ısmicas Material V[m/s] Capa Meteorizadas

300-900

Aluviones Modernos

350-1500

Arcillas

1000-2000

Areniscas

1400-4500

Conglomerados

2500-5000

Calizas

4000-6000

Dolomias

5000-6000

Sal

4500-6500

Los modelos creados fueron los siguientes: el primero de ellos est´a constituido de dos capas paralelas; el segundo con cuatro capas paralelas; el tercer modelo est´a constituido por cuatro capas no uniformes y el ´ultimo presenta una lente asociada a domos salinos, trampas de rocas, entre otros en el interior de una de sus capas. Todos los modelos fueron desarrollados en una extensi´on (distance) de 0 a

10 kil´ometros con respecto al eje x y una profundidad (depth) de 0 a 7 kil´ometros

con respecto al eje z. La explicaci´on del c´odigo fuente utilizado para la generaci´on

p. 32

de los modelos se puede observar en el anexo C.

El primer modelo creado, consta de dos capas isotr´opicas y homog´eneas en forma paralela (Figura 6). La primera capa esta limitada entre 0 y 3,5 km de profundidad, con un valor de sloth de 0,9 [s2/m2]. De igual manera, la segunda capa esta deﬁnida entre 3,5 y 7 km de profundidad y su valor de sloth fue de 0,3 [s2/m2]. La ﬁnalidad de este dise˜no es obtener el conocimiento b´asico del c´odigo fuente para la generaci´on del modelo.

Figura 6. Modelos de dos capas paralelas.

El segundo modelo consta de cuatro capas isotr´opicas y homog´eneas en forma paralela (Figura 7). La estructura de las capas y los respectivos valores de velocidad son los siguientes: la primera est´a comprendida entre 0 y 1,8 km y el valor del sloth es de 0,9 [s2/m2]; la segunda capa entre 1,8 y 3,5 km y un valor de velocidad de 0,7 [s2/m2]; la tercera capa entre 3,5 y 5,3 km con un valor de velocidad de 0,5 [s2/m2]

p. 33

de sloth y la ´ultima capa est´a deﬁnida entre 5,3 y 7 km con un valor de velocidad de 0,3 [s2/m2]. El prop´osito de este modelo era desarrollar una modiﬁcaci´on al adicionar dos capas uniformes y paralelas al primer modelo.

Figura 7. Modelos de cuatro capas paralelas.

El tercer modelo est´a constituido por cuatro capas no uniformes (Figura 8) con diferentes valores de velocidad. Fue el modelo de referencia implementado en el proceso de inversi´on s´ısmica en el seminario de investigaci´on. Las capas estas constituidas por coordenadas en un plano (x;z). En la tabla 2 se aprecian los valores de sloth utilizados en el modelo del subsuelo de cuatro capas no uniformes y las coordenadas de cada una de las capas.

p. 34

Figura 8. Modelos de cuatro capas no uniformes.

Tabla 2. Valores de sloth para cada una de las capas del tercer modelo Capas Sloth [s2/m2] Coordenadas (x;z) Capa uno

0,9

(0;3),(0;0),(10;0),(10;3), (9,5;3), (5,5;2) y (3;2,8) Capa dos

0,7

(0;3), (3;2,8), (5,5;2), (9,5;3), (10;3), (10;6), (7;5), (3;4,9) y (0;5) Capa tres

0,5

(3;4,9), (7;5), (10;6), (10;7) y (7;7) Capa cuatro

0,3

(0;5), (3;4,9), (7;7) y (0;7) El ´ultimo modelo desarrollado (Figura 9) presenta las mismas caracter´ısticas del modelo descrito anteriormente, con una diferencia detallada por una burbuja o domo salino ubicado en la segunda capa que ´esta constituido por las siguientes coordenadas: (4;4), (5,5;3,5), (7;4) y (5,5;4,5).

p. 35

Figura 9. Modelos de cuatro capas no uniformes con una variante.

2.3

ADQUISICI ´ON S´ISMICA EN SU

Durante esta secci´on se abarca la simulaci´on de procesos de adquisici´on s´ısmica desde el software SU. Entonces, el objetivo principal es brindar los conceptos y las herramientas computacionales disponibles para la generaci´on de trazados de rayos y seguidamente, para el posicionamiento de fuentes de perturbaci´on y sus receptores durante la obtenci´on de una adquisici´on s´ısmica.

2.3.1.

Trazado de rayos El trazado de rayos s´ısmicos es ampliamente utilizado como una forma de exploraci´on del subsuelo. Se fundamenta en determinar el recorrido que realiza una onda s´ısmica desde un punto donde es generada (fuente perturbadora o s´olo fuente) has-

p. 36

ta el receptor (ge´ofono). Entonces, este trazado permite ver el comportamiento de la onda que viaja por el subsuelo, siendo as´ı la primera experiencia de observaci´on en cuanto a las condiciones del terreno.

En SU, el comando para el trazado de rayos es triray. Para la generaci´on de los rayos, triray requiere de un par´ametro de entrada: siendo el modelo generado en la secci´on 2.1 del presente capitulo (v´ease Figura 8. Despu´es de ingresarlo se deﬁnen los ´angulos de disparo, y a trav´es de triangulaci´on se obtiene el trazado, teniendo entonces que:

• triray: trazado din´amico de rayos para modelos basados en triangulaci´on [3].

triray < modelﬁle > rayends [par´ametros opcionales].

Cabe aclarar que la ventaja de usar triray es que pueden encontrarse errores en la etapa de modelado, como por ejemplo, que hayan superﬁcies no cerradas o constru´ıdas inadecuadamente. Por esto se utiliza esta herramienta brindada por SU y de forma m´as detallada, se citan los par´ametros opcionales de triray [3]:

• xs: coordenada de la fuente en direcci´on del eje horizontal (superﬁcie de la

tierra).

• zs: coordenada de la fuente en direcci´on del eje vertical nangle, n´umero de

´angulos de salida.

• fangle: primer ´angulo de salida (en grados).

• langle: ´ultimo ´angulo de salida (en grados).

• nxz: n´umero de coordenadas (x,z) en el archivo de rayo.

p. 37

• nangle: n´umero de ´angulos que se desean observar.

Seguido del trazado, viene la etapa de agrupaci´on de los rayos con el modelo generado (secci´on de modelado de capas) y para esto se utilizan los siguientes comandos:

• psgraph: El cual crea un archivo PostScript(.eps) en base a un archivo binario.

• psmerge: Encargado de unir archivos PostScript.

Esto quiere decir que al ﬁnal de esta etapa se hace un empalme entre el modelo de capas y el trazado generado. Durante esta operaci´on se crea la imagen de c´omo posiblemente se propagar´an las ondas durante la adquisici´on. A continuaci´on se presenta el trazado de rayos sobre 3 modelos de capas del subsuelo (v´eanse las Figuras 10, 11 y 12), los cuales presentan ciertas similitudes pero el ´ultimo tiene un lente o intrusi´on.

Figura 10. Trazado de rayos sobre el modelo de capas paralelas. (a) Trazado con ´angulos de igual apertura.

(b) Trazado con una direcci´on espec´ıﬁca.

p. 38

Figura 11. Trazado de rayos sobre el modelo 3 (modelo de 4 capas del subsuelo). (a) Trazado sobre la primera capa del modelo 3.

(b) Trazado sobre la segunda capa del modelo 3.

p. 39

Figura 12. Trazado de rayos sobre el modelo 4 (modelo de 4 capas con intrusi´on). (a) Trazado sobre la primera capa del modelo 4.

(b) Trazado sobre la segunda capa del modelo 4.

(d) Trazado sobre la cuarta capa del modelo 4.

(e) Trazado sobre la quinta capa del modelo 4.

p. 40

2.3.2.

Adquisici´on y trazas s´ısmicas Esta secci´on se enfoca en a la obtenci´on los datos s´ısmicos, o sea que aqu´ı es cuando se tendr´a la primera simulaci´on de perturbaci´on y obtenci´on de datos s´ısmicos. Por esto se siguen los fundamentos plasmados en el el libro Seismic Data Processing with Seismic Un*x [3] donde pueden encontrarse el dise˜no de la adquisici´on, los par´ametros a tener en cuenta y la obtenci´on de trazas s´ısmicas sobre un modelo de capas. Entonces se siguen algunos de los ejemplos contenidos all´ı pero aclarando que, para nuestro seminario s´olo se efectuar´a una adquisici´on con un disparo y sin ninguna inclinaci´on espec´ıﬁca en cuanto a la direcci´on de las ondas s´ısmicas. El objetivo de la adquisici´on entonces es, brindar informaci´on de los cambios de velocidades que experiment´o la onda s´ısmica desde el momento en que sali´o de la fuente hasta llegar al receptor o ge´ofono; ese resultado se puede evidenciar a trav´es de las trazas s´ısmicas, por ejemplo, a continuaci´on se presentan dos formas en las que SU permite visualizarlas (Figura 13).

Figura 13. Representaciones de las trazas s´ısmicas.

(a) Trazas s´ısmicas en forma de imagen sencilla.

(b) Trazas s´ısmicas en forma de mapas de bits.

p. 41

En esta se relaciona la cantidad de tiempo que han grabado cada uno de los ge´ofonos, o sea que su eje horizontal se sit´uan cada uno de ellos y en su eje vertical los instantes de tiempo en que se captaron las perturbaciones. Aunque para las representaciones utilizan los mismos datos, las diferencias entre las trazas se fundamentan en la interpretaci´on bajo las cuales se est´e realizando el an´alisis geof´ısico. Por ejemplo, la representaci´on sencilla es com´unmente utilizada cuando se realizan an´alisis de propagaciones con frente de onda (permite identiﬁcar la forma de las ondas reﬂejadas y las difractadas), mientras que en forma de bits se utiliza para el trazado de rayos (permite determinar qu´e est´a captando cada ge´ofono en forma espec´ıﬁca).

En el anexo D se presentan: el c´odigo utilizado para realizar el proceso de adquisici´on s´ısmica a trav´es de SU, los par´ametros de mayor relevancia de este proceso y algunas trazas s´ısmicas obtenidas.

p. 42

M ´ETODOS DE SOLUCI ´ON METAHEUR´ISTICOS

3.1

INTRODUCCI ´ON

En este cap´ıtulo se presentan dos t´ecnicas metaheur´ısticas: Simulated Annealing (SA) y Artiﬁcial Bee Colony (ABC), utilizadas en el seminario de investigaci´on. Se realiza una descripci´on de las caracter´ısticas principales de SA y ABC, como lo son, sus analog´ıas f´ısicas, descripci´on de los par´ametros de mayor relevancia, diagramas de ﬂujo y sus respectivos pseudoc´odigos. De igual manera, se desarrolla la evaluaci´on de ambas t´ecnicas metaheur´ısticas en funciones matem´aticas multidimensionales para veriﬁcar su funcionamiento, obtener los par´ametros principales y realizar la implementaci´on en el proceso de inversi´on s´ısmica. Finalizando, se muestran los pasos en el proceso de vinculaci´on entre MATLAB y SU para resolver el problema planteado en seminario de investigaci´on.

3.2

M ´ETODOS DE SOLUCI ´ON

Los problemas de optimizaci´on pueden resolverse a trav´es de m´etodos exactos y m´etodos heur´ısticos1. En los m´etodos exactos se conoce la funci´on objetivo, la continuidad de la funci´on y su respectiva derivada. No obstante, si no se cumple con alguna de las condiciones mencionadas se procede a utilizar un m´etodo heur´ıstico. 1Son una metodolog´ıa y/o conjunto de pasos ordenado, son basados en el uso de reglas emp´ıricas para la b´usqueda de una soluci´on aproximada de un problema en particular

p. 43

Los m´etodos heur´ısticos permiten encontrar soluciones que dependen de su punto de partida y presenta un avance evolutivo con cada paso iterativo. Para la optimizaci´on de problemas y c´alculos de alta complejidad se han desarrollado m´ultiples t´ecnicas, ´estas son clasiﬁcadas en exactas y aproximadas. Las exactas se caracterizan por garantizar la b´usqueda de la soluci´on ´optima en cualquier problema, su principal inconveniente es el crecimiento exponencial del tiempo de c´omputo en su resoluci´on. Las aproximadas sacriﬁcan la garant´ıa de encontrar el resultado ´optimo de un problema a cambio de obtener una buena soluci´on en un tiempo razonable. Durante las tres ´ultimas d´ecadas se desarrollaron y utilizaron tres tipos de t´ecnicas aproximadas, se les conoce como: m´etodos constructivos, que parten de una soluci´on vac´ıa y en su evoluci´on se construye su estructura; m´etodos de b´usqueda local que usan el concepto de vecindario y recorren el espacio de b´usqueda hasta encontrar un m´ınimo local y las t´ecnicas metaheur´ısticas se fundamentan en la combinaci´on de diferentes m´etodos heur´ısticos para conseguir una exploraci´on del espacio de b´usqueda m´as eﬁciente. En la Figura 14 se ubica el esquema que relaciona los m´etodos de soluci´on con las t´ecnicas metaheur´ısticas que son el objetivo principal de este cap´ıtulo.

p. 44

Figura 14. Descripci´on de los m´etodos de soluci´on.

Existe una gran variedad de formas para la clasiﬁcaci´on de las t´ecnicas metaheur´ısticas, como lo son: basadas en la naturaleza (algoritmos bioinspirados) o no basados en ella, basadas en la memoria o sin memoria, con funci´on objetivo est´atica o din´amica, entre otras. Las t´ecnicas metaheur´ısticas pueden utilizar un ´unico punto de partida o trabajar sobre un conjunto o poblaci´on en un espacio de b´usqueda, estas t´ecnicas se dividen en dos ramas, las cuales son las t´ecnicas basadas en la trayectoria y las t´ecnicas basadas en la poblaci´on. Las t´ecnicas basadas en la trayectoria son aquellas que parten de un punto inicial y generan un proceso de actualizaci´on mediante la exploraci´on del vecindario, formando una trayectoria. Por otro lado, las t´ecnicas basadas en la poblaci´on trabajan con un conjunto de individuos que representan varias soluciones en el espacio de b´usqueda, su resultado depende fundamentalmente en la evoluci´on de la poblaci´on en cada paso o itera-

p. 45

ci´on.

A lo largo de toda la revisi´on bibliogr´aﬁca que se realiz´o antes de dar inicio al seminario de investigaci´on, se encontraron gran cantidad de t´ecnicas metaheur´ısticas, como lo son: Simulated Anneling [4], Artiﬁcial Bee Colony [5], Particle Swarm Optimization [6], Genetic Algorithms [7] y Intelligent Water Drops [8]. A cada una de ellas se le hizo su respectiva analogica f´ısica, el estudio de su estructura y su funcionamiento. Por medio de este an´alisis se pudo seleccionar las t´ecnicas adecuadas para implementarlas en el problema de inversi´on s´ısmica. Las t´ecnicas seleccionadas fueron Simulated Annealing (SA) y Artiﬁcial Bee Colony (ABC). En el caso de SA ya hab´ıa sido aplicado en un problema de inversi´on s´ısmica en una dimensi´on por Mrinal y Stoffa en el a˜no 1995 [9]. Tambi´en, fue estudiado e implementado por el grupo CEMOS2 en un trabajo de investigaci´on sobre el desarrollo de sistemas de ecuaciones. [10]. Igualmente, el algoritmo de ABC tambi´en ha sido implementado por el grupo CEMOS en algunos trabajos de investigaci´on [11], [12] y [13] pero no aplicado a la problem´atica de la inversi´on s´ısmica. La elecci´on de las t´ecnicas SA y ABC fue basada en la informaci´on previa obtenida por medio de la revisi´on bibliogr´aﬁca, por las implementaciones ya realizadas en problemas de optimizan global, por la programaci´on de script en la plataforma MATLAB y para efectuar un an´alisis comparativo entre una t´ecnica de b´usqueda basada en trayectoria con un solo punto o part´ıcula (SA) y una t´ecnica de b´usqueda basada en un conjunto de puntos o poblaci´on (ABC).

Recientemente, como consecuencia de la evoluci´on de los computadores, los m´etodos de optimizaci´on global se han aplicado a varios problemas geof´ısicos. A diferencia de los m´etodos de optimizaci´on local, estos tratan de encontrar el m´ınimo 2Grupo de Investigaci´on en control, electr´onica, modelado y simulaci´on de la Universidad Industrial de Santander

p. 46

global de la funci´on objetivo. Algunos de los algoritmos de optimizaci´on global son de naturaleza estoc´astica y utilizan informaci´on m´as general dentro del espacio de b´usqueda para actualizar su posici´on. Si bien la convergencia de estos m´etodos no est´a garantizada para cada problema de optimizaci´on, algoritmos como SA, ha reportado resultados relativamente conﬁables en la soluci´on del problema de inversi´on s´ısmica [9]. De igual manera, el algoritmo ABC presenta caracter´ısticas interesantes como la especializaci´on de agentes y la posibilidad de identiﬁcar regiones ´optimas [14].

3.3

SIMULATED ANNEALING

La t´ecnica de SA fue formulada por Scott Kirkpatrick, Daniel Gelllat y Mario Vecchi en el a˜no de 1983 [4]. Esta t´ecnica est´a inspirada en el proceso de enfriamiento de los materiales. En este proceso, a medida que baja la temperatura el material va modiﬁcando su conﬁguraci´on y cada una de las conﬁguraciones tiene asociada una energ´ıa determinada. El algoritmo de SA proviene del proceso de recocido3 del acero y cer´amicas. Se inicia con un calentamiento en los s´olidos que causa que los ´atomos aumenten su energ´ıa y puedan desplazarse de sus posiciones iniciales. Un enfriamiento lento presenta mayores probabilidades de alcanzar una conﬁguraci´on de cristalizaci´on. Mientras que, con una velocidad de enfriamiento r´apida se logran obtener s´olidos no cristalinos a causa de que partiendo de estados energ´eticos elevados la temperatura desciende r´apidamente y las mol´eculas no pueden adoptar una conﬁguraci´on m´as estable.

3Tratamiento t´ermico donde su ﬁnalidad es el ablandamiento y la recuperaci´on de estructuras en los materiales

p. 47

3.3.1.

Descripci´on de par´ametros SA resuelve gran cantidad de problemas de optimizaci´on, para ellos es necesario ser precisos a la hora de escoger las etapas del algoritmo, las ecuaciones y variables espec´ıﬁcas. El ´exito de SA depende del uso de sus par´ametros propios, como lo son: las probabilidades iniciales y ﬁnales, de ellas se generan la temperatura inicial y ﬁnal las cuales dan inicio y terminaci´on al proceso de recocido (en las ecuaciones

3 y 4 se muestran la relaciones entre la temperatura inicial y ﬁnal con respecto a

las probabilidades iniciales (pi) y ﬁnales (pf) y el factor de enfriamiento, el cual es el encargado de disminuir la temperatura en cada paso evolutivo de SA(ecuaci´on 5). Ti = −1 log(pi)

(3)

Tf = −1 log(pf)

(4)

Frac = Tf Ti

n−1

(5)

Las iteraciones se hacen de forma sucesiva, evaluando los estados vecinos hasta encontrar una aproximaci´on optima, la cual debe cumplir con los requisitos espec´ıﬁcos del problema. En cada nueva iteraci´on el par´ametro de temperatura va disminuyendo simulando el proceso de enfriamiento de un s´olido (metal) hasta llegar a una conﬁguraci´on cristalina, que en nuestro caso representa el m´ınimo global de una funci´on objetivo.

p. 48

3.3.2.

Diagrama de ﬂujo de SA Los pasos para el desarrollo del c´odigo se ilustran de manera gr´aﬁca (Figura 15) y se nombran a continuaci´on:

• Deﬁnir los valores de la probabilidad inicial y la probabilidad ﬁnal.

• Generar un nuevo vecino.

• Evaluar el vecino.

• Hallar la diferencia entre la nueva evaluaci´on y la evaluaci´on actual.

• Aceptar o rechazar la nueva soluci´on.

• En casa d rechazar la nueva soluci´on, hallar una valor de probabilidad.

• Generar un numero aleatorio entre [0,1].

• Comparar el valor de la nueva probabilidad y el numero aleatorio para aceptar

o rechazar la nueva soluci´on.

• La temperatura disminuye con cada nueva iteraci´on.

• El proceso llega a su ﬁnal cuando alcanza el criterio de parada.

• Se entrega la soluci´on obtenida.

p. 49

Figura 15. Diagrama de ﬂujo de SA.

p. 50

3.3.3.

Descripci´on e implementaci´on de las funciones de prueba en SA En la Tabla 3 se presentan las caracter´ısticas de las funciones de prueba [15] utilizadas para comprobar el funcionamiento de SA. En la tabla se incluye el nombre de cada una de las funciones, su formulaci´on, su rango y valor ´optimo, esto con el ﬁn de presentar mayor informaci´on sobre las funciones prueba y poder analizar de mejor manera los resultados obtenidos. Estas funciones presentan algunas diferencias y similitudes en ciertos aspectos. Una de las diferencias se puede evidenciar con la funci´on Ratrigin, que presenta gran variedad de m´ınimos locales y la funci´on Rosenbrock, que es la constituci´on de un valle.

Para lograr un buen funcionamiento de SA aplicado en el problema de inversi´on s´ısmica el primer paso es comprobarlo mediante un n´umero determinado de funciones de prueba, con las cuales se veriﬁca la eﬁciencia en alcanzar el m´ınimo global establecido por dichas funciones.

p. 51

Tabla 3. Tabla de funciones prueba Funci´on Formulaci´on Rango ´Optimo Booth f1(x, y) = (x + 2y −7)2 + (2x + y −5)2 −10 ≤(x, y) ≤10 f(1, 3) = 0 f2(x, y) = a(y −bx2 −cx −r)2+ Branin s(1 −t)cos(x) + s −10 ≤x ≤10 f(9.42478,2.475)=

397887

a = 1; b =

5,1

4π2 ; c = 5 π ; r = 6; s = 10; t =

8π

0 ≤y ≤15

Bukin f3(x, y) = 100p |y + 0,01x2| + 0,01|x + 10| −15 ≤x ≤−5 f(-10,1)=0 −3 ≤y ≤3 Rastrigin f5(xi) = 10n + Pn i=1[x2 i+1 −10cos(2πxi)] −5,12 ≤x ≤5,12 f(0,0)=0 Rosenbrock f6(xi) = Pn−1 i=1 [100(xi+1 −x2 i )2 + (1 −xi)2] −∞≤x ≤∞ f(1,1)=0 Los par´ametros de importancia en el algoritmo son: la probabilidad inicial, el valor de normalizaci´on, n´umero de iteraciones, n´umero de ensayos por cada iteraci´on y la probabilidad ﬁnal. El ´unico par´ametro establecido y modiﬁcado durante las pruebas realizadas fue el de la probabilidad inicial, la cual tambi´en puede ser conocida como la temperatura inicial para que el proceso de enfriamiento comience, con respecto a su analog´ıa f´ısica. Las tres probabilidades usadas fueron: 0,2 (enfriamiento r´apido), 0,5 (enfriamiento moderado) y 0,7 (enfriamiento lento), implementando SA en catorce funciones de prueba, entre las cuales se mencionan algunas de ellas como: Branin, Booth, Rosenbrock, Rastrigin y Bukin (tabla 3).Las pruebas en las cuales siempre se obtuvo el resultado correcto y esperado se muestran en la Figura 16.

p. 52

Figura 16. Funciones prueba en SA

(a) Branin (b) Booth (c) Bukin (d) Rastrigin (e) Rosenbrock Para el desarrollo del experimento con SA en la implementaci´on con las funciones prueba para corroborar su ´optimo funcionamiento en la localizaci´on del m´ınimo global, se utilizaron tres valores para la probabilidad inicial como punto de partida, se desarrollaron 10 repeticiones para cada una de las probabilidades iniciales, cada repetici´on constaba de 80 iteraciones. Cuando SA detecta que la variaci´on en las temperaturas es cero el se detiene,siendo este el criterio de parada implementado. Se tomaron como datos de mayor relevancia para el criterio de selecci´on de los va-

p. 53

lores de los par´ametros con mayor eﬁciencia, el tiempo de c´omputo y la proximidad de la part´ıcula a las coordenadas del m´ınimo global en cada una de las funciones prueba.

Para hacer la evaluaci´on del algoritmo de SA se utiliz´o MATLAB y un equipo de c´omputo cuyas caracteristicas se muestran en la Tabla 4.

Tabla 4. Especiﬁcaciones del equipo de c´omputo usado en los experimentos sobre las funciones de prueba Procesador Intel R ⃝CoreTMi5-240M CPU @ 2.30GHz

RAM

4.00 GB

Sistema operativo Microsoft R ⃝Windows R ⃝7 Home Premium de 64 bits Alimentaci´on Conectado a la red el´ectrica Plan de energ´ıa Alto rendimiento En la Tabla 5 se muestran los par´ametros bajo los cuales se obtuvieron los mejores resultados en la implementaci´on de SA en el conjunto de funciones de prueba y en la Tabla 6 se presenta los promedios y desviaciones est´andar de los experimentos realizados con cada una de las funciones de prueba. En el seminario de investigaci´on se gener´o una discusi´on con cada una de las funciones con respecto a la selecci´on de los par´ametros que ser´an utilizados en el proceso de inversi´on s´ısmica con SA. Ademas, se encuentra que existe cierta tendencia del algoritmo para encontrar el m´ınimo de la manera m´as ´optima bajo la probabilidad inicial entre 0,5 y 0,7, por esto se tiene como preferencia un valor entre estos para la siguiente fase del seminario.

p. 54

Tabla 5. Par´ametros bajo los cuales se obtuvieron los mejores resultados en la implementaci´on de SA en el conjunto de funciones de prueba Funci´on Probabilidad inicial Probabilidad ﬁnal Booth

0,5

0,001

Branin

0,2

0,001

Bukin

0,7

0,001

Rastrigin

0,7

0,001

Rosenbrock

0,7

0,001

Tabla 6. Promedios y desviaci´on de las funciones prueba de SA Funci´on Probabilidad M´ınimo promedio Desviaci´on Tiempo de c´omputo [s] Iteraciones

0,2

12,5893449

17,5

Booth

0,5

20,07959

0,25168568

27,883298

39,1

0,7

36,755046

51,7

0,2

27,7029

9,2844217

9,4

Branin

0,5

27,7029

28,2459699

35,2

0,7

27,7029

42,9912697

51,7

0,2

6,20101

2,52656112

16,6805347

21,5

Bukin

0,5

3,98197

2,68709655

32,6001855

45,2

0,7

5,08931

1,93219435

38,5398945

0,2

-1,51411

0,02145773

9,3733725

7,1

Rastrigin

0,5

2,52316

1,3645464

37,4496792

35,4

0,7

1,876054

1,51227954

51,7390166

51,6

0,2

1,775916

2,21133478

7,9211243

8,7

Rosenbrock

0,5

3,906229

3,31564537

17,7830898

22,6

0,7

2,919689

4,00180246

34,7112264

46,6

Con cada valor se determin´o el n´umero iteraciones que la part´ıcula utiliza para alcanzar la proximidad a las coordenadas del m´ınimo global y el tiempo de c´omputo que tarda dicha operaci´on. Seg´un sea la forma de la funci´on de prueba, el algoritmo puede presentar similitudes con cada una de las probabilidades iniciales, sus ´unicos cambios se presentan en el n´umero de iteraciones y el tiempo de c´omputo utilizado en la simulaci´on. En las funciones Branin y Booth (Figura 16a y 16b), donde el valor

p. 55

de las coordenadas del m´ınimo global son (1,1), se presento un comportamiento similar con respecto a cada valor de probabilidad usado con SA en las iteraciones y tiempo de computo empleado.

La probabilidad inicial de 0,7 fue la que present´o una mayor cobertura con la cercan´ıa de las coordenadas al valor del m´ınimo global, el n´umero de iteraciones es un valor cercano al l´ımite establecido de 80, demostrando el funcionamiento del criterio de parada establecido con el valor de la probabilidad ﬁnal. El tiempo de c´omputo no excede los 60 segundos para lograr la localizaci´on con mayor proximidad a las coordenadas al m´ınimo global (Tabla 6). Se advierte que en algunos casos la probabilidad inicial de 0,2 presenta un gran acercamiento al m´ınimo global, pero dado el caso de su proceso de un enfriamiento r´apido, la exploraci´on que realiza es r´apida y puede ocasionar alg´un tipo de estancamientos en un m´ınimo local (ver las Tablas E 4, E 5 y E 6 del Anexo E).

3.4

ARTIFICIAL BEE COLONY

En la naturaleza muchas especies de insectos viven en grupos de colonias o enjambres, lo que les facilita realizar sus tareas colectivas. Si bien, cada miembro del grupo presenta capacidad limitada para realizar tareas, la colonia exhibe un comportamiento emergente que surge de la interacci´on entre los individuos [5]. ABC es un algoritmo de inteligencia de enjambre desarrollado por Dervis Karaboga en 2005 [5], que busca emular el comportamiento de b´usqueda y explotaci´on de fuentes de alimento de las abejas con el ﬁn de encontrar

p. 56

buenas soluciones a problemas de optimizaci´on de alta complejidad (t´ecnicas metaheuriisticas). El algoritmo recrea una colmena artiﬁcial conformada por tres tipos de abejas employed, scouts y onlookers(obreras, exploradoras y observadoras) que establecen una zona de comunicaci´on e implementan un proceso de recolecci´on de alimento (polen). El proceso de b´usqueda de las fuentes de alimento es de manera aleatoria.

3.4.1.

Descripci´on de par´ametros ABC es una t´ecnica metaheur´ıstica dise˜nada para localizar el m´ınimo global de una funci´on objetivo en un espacio determinado, conocido como espacio de b´usqueda. El representa y corresponde a un entorno y cada punto de ese entorno representa una posible soluci´on (fuente de alimento), que el grupo de abejas artiﬁciales pueden utilizar. De igual forma, el ﬁtness (calidad) del polen de una fuente de alimento representa la cercan´ıa a una posible soluci´on. Dentro de la colonia artiﬁcial existen tres tipos de abejas:

• Employed: Explotan cada una de las fuentes especiﬁcas de alimento que han

sido exploradas con anterioridad y comparten la informaci´on sobre la calidad y cantidad de alimento al resto de la colonia.

• Scouts: Son las abejas encargadas de la b´usqueda de nuevas fuentes de

alimento en un espacio aleatorio.

• Onlookers: Son las responsables de recibir la informaci´on acerca de las

fuentes de alimento y escogen una seg´un la informaci´on de la calidad y cantidad.

La mitad de poblaci´on de la colonia artiﬁcial esta constituida por employed y la segunda mitad por scouts y onlookers. Cada uno de los ciclos de b´usqueda en ABC

p. 57

consiste en tres pasos: el envi´o de employed a las respectivas fuentes de alimento para medir el ﬁtness (en la ecuaci´on 6 se expresa su formulaci´on matem´atica) de cada una de ellas; la selecci´on de las fuentes de alimento por parte de las onlookers despu´es de compartir la informaci´on obtenida de las employed para determinar la cantidad alimento de las fuentes y la ﬁjaci´on de las scouts para ser enviadas a posibles nuevas fuentes de alimento.

fiti =         

1+fi si fi ≥0

1 + |fi|

si fi < 0

(6)

fi, corresponde al valor de la funci´on objetivo evaluada en los puntos del espacio de b´usqueda donde el ﬁtness del alimento es el mejor.

ABC genera una poblaci´on inicial con una distribuci´on aleatoria de i soluciones, donde i = 1,..., SN, j = 1,...,D, donde SN es el n´umero de fuentes de alimento e indica el tama˜no de la poblaci´on y D es el numero de par´ametros que se requieren optimizar. En cada ciclo del algoritmo, en el espacio de b´usqueda deﬁnido entre xmin ≤x ≤xmax, cada una de las employed determina una nueva fuente de alimento vecina a su fuente actual y calcula la cantidad de alimento que la nueva fuente posee mediante la ecuaci´on 7.

vij = xij + θij(xij −xki)

(7)

p. 58

θij es un n´umero aleatorio entre [-1,1], k̸ = i ∈{1,2,...,D}, xij es el par´ametro j-esimo de una soluci´on xi que fue seleccionada para ser modiﬁcada. Si la cantidad de n´ectar de esta nueva fuente de alimento es mayor a la actual, la employed se mueve a su nueva fuente de alimento. De lo contrario se mantiene. A las fuentes de alimento que han sido olvidadas se les denomina abandonadas, y si transcurren muchas iteraciones sin que la soluci´on mejore, la employed ubicada en esa fuente se convierte en una scout. En esta posici´on una nueva soluci´on es generada aleatoriamente para la scout, ecuaci´on 8:

xj i = xj min + u(xj max −xj min)

(8)

j se determina de manera aleatoria y debe ser diferente de i y u, es un n´umero aleatorio entre [-1,1].

Una onlooker elige una fuente de alimento dependiendo del valor de probabilidad (pi) asociado con esta ´ultima (ecuaci´on 9).

pi = fiti

PSN

n=1 fitn

(9)

donde SN, como ya se mencion´o, es el numero de fuentes de alimento, y fiti es el valor del ﬁtness de la soluci´on.

p. 59

Una vez que toda las employed culminan el proceso de b´usqueda, comparten la informaci´on de cada una de las fuentes de alimento y su respectiva ubicaci´on con las onlookers, quienes sean encargadas de seleccionar la fuente de alimento en funci´on del ﬁtnees. Las employed y scouts retienen la posici´on de la fuente que ofrece mayor alimento a la colonia.

Finalmente, las scouts se encargan de reemplazar todas las fuentes de alimento que ya han sido agotadas con nuevas fuentes seleccionadas al azar. El algoritmo asume que en cada ciclo de trabajo realizado, una ´unica fuente puede ser agotada y la employed pasa a ser una scouts, por lo que, si m´as de una fuente ha excedido el valor l´ımite para ser abandonada y olvidada, debe ser elegida aquella fuente que tiene el mayor valor [16].

3.4.2.

Diagrama de ﬂujo de ABC Los pasos para el desarrollo del c´odigo se ilustran de manera gr´aﬁca (Figura 17) y se nombran a continuaci´on:

• Generar un conjunto (poblaci´on) de soluciones (individuos) al problema.

• Evaluar cada soluci´on en la funci´on objetivo.

• Seleccionar las mejores soluciones de la poblaci´on seg´un el valor del ﬁtness.

• Generar nuevas soluciones a partir de las mejores seg´un su valor de ﬁtness.

• Evaluar las nuevas soluciones.

• Escoger las soluciones que entrar´an a la siguiente iteraci´on.

• Preguntar acerca de la distribuci´on de las scouts.

p. 60

• Memorizar la mejor soluci´on.

• Veriﬁcar si existen fuentes abandonadas y asignar un nuevo punto si existen.

• Preguntar si el criterio de parada se ha cumplido para concluir el proceso. De

lo contrario el proceso retorna a buscar una mejor soluci´on.

p. 61

Figura 17. Diagrama de ﬂujo de ABC.

p. 62

3.4.3.

Descripci´on e implementaci´on de las funciones de prueba en ABC La versi´on utilizada del algoritmo de ABC es la propuesta por Karaboga [17]. El proceso de veriﬁcaci´on del funcionamiento de ABC se se desarrollo en catorce funciones de prueba , entre las cuales se mencionan algunas de ellas como: Brainin, Both, Rosenbrock, Rastrigin y Bukin como se han mencionado en la secci´on de descripci´on e implementaci´on de funciones prueba en SA(Figura 16 y Tabla 3). En la Tabla 7 se muestran los par´ametros bajo los cuales se obtuvieron los mejores resultados en la implementaci´on de ABC en el conjunto de funciones de prueba y en la Tabla 8 se presenta los promedios y desviaciones est´andar de los experimentos realizados con cada una de las funciones de prueba. Los valores de la tabla generaron discusiones con respecto a la selecci´on de los par´ametros que serian utilizados en el proceso de inversi´on s´ısmica con ABC. Las discusiones se encuentran en las memorias digitales del documento. En las discusiones se resaltaron los valores de las tres pruebas realizadas con las poblaciones de 10, 20 y 30 abejas. Tabla 7. Par´ametros bajo los cuales se obtuvieron los mejores resultados en la implementaci´on de ABC en el conjunto de funciones de prueba Funci´on Tama˜no de la poblaci´on Booth

Branin

Bukin

Rastrigin

Rosenbrock

p. 63

Tabla 8. Promedios y desviaci´on de las funciones prueba de ABC Funci´on Poblaci´on M´ınimo promedio Desviaci´on Iteraciones promedio Tiempo promedio [s]

1,12919007

1,146559107

23,2

0,647

Booth

0,23016737

0.051349396

22,1

0,906

0,012744912

0,013532305

19,8

1,137

0,36216399

0,352536372

25,9

0,696

Branin

0,21329765

0,15553577

22,3

0,878

0,22446724

0,132097387

19,5

1,165

8,6784274

7,758857266

18,5

0,638

Bukin

7,5289446

6,937159233

16,3

0,795

3,5909497

2,656896372

16,4

0,997

1,48556586

1,478197364

21,9

0,643

Rastrigin

0,665324746

1,013439495

25,3

0,925

0,204889333

0,21624768

1,176

2,4230653

2,771242645

15,6

0,597

Rosenbrock

1,72072109

1,98449941

14,7

0,748

0,333023497

0,758225293

15,3

0,978

El ´unico par´ametro que exige el algoritmo para su funcionamiento es el de seleccionar un tama˜no de poblaci´on adecuada para el proceso de exploraci´on y b´usqueda dentro del algoritmo, siendo la principal caracter´ıstica de ABC, donde solo con el n´umero de abejas necesarias se logra alcanzar el m´ınimo global de la funci´on objetivo en un tiempo no superior a los 2 [s] y siempre con una aproximaci´on a las coordenadas del m´ınimo de la funci´on de prueba. Esta cercan´ıa de la poblaci´on a las coordenadas del m´ınimo global de la funci´on fue la selecci´on de criterio de eﬁciencia del par´ametro de poblaci´on inicial en cada prueba realizada, puede observarse en el anexo F.

Para el desarrollo del experimento de evaluaci´on de ABC con el grupo de funciones de pruebas previamente seleccionadas, donde su objetivo es el de encontrar su ´optimo funcionamiento en la localizaci´on del m´ınimo global en cada una de las funciones de prueba. Se utilizaron tres valores para el tama˜no de la poblaci´on, se desarrollaron 10 ensayos para cada una de las poblaciones seleccionadas, en cada

p. 64

ensayo se produc´ıa una rutina de 5 repeticiones, cada e constaba de 1000 iteraciones como l´ımite en el desarrollo del algoritmo, el tiempo de c´omputo y la proximidad a las coordenadas del m´ınimo son los factores m´as importantes en la toma de los resultados obtenidos en cada una de las simulaciones realizadas. Los valores de la poblaci´on seleccionados fueron de: 10, 20 y 30 abejas, con cada valor de poblaci´on se tomaron los datos con respecto al n´umero de iteraciones que el algoritmo utiliza para alcanzar el m´ınimo, la proximidad de la part´ıcula al m´ınimo global y el tiempo de todo el proceso de la simulaci´on. El tiempo de c´omputo es un dato relevante a causa de la gran velocidad de convergencia de ABC en comparaci´on con SA. La poblaci´on indicada para una excelente b´usqueda y exploraci´on en cualquier funci´on objetivo es entre 20 y 30 abejas, con estos n´umeros se apreciaron los resultados con mayor precisi´on como puede apreciarse en las funciones Branin, Booth, Rastrigin y Rosenbrock (Ver Tabla 8).

3.5

INTERFAZ ENTRE MATLAB Y SEISMIC UNIX

En el desarrollo del seminario de investigaci´on se cre´o y se modiﬁc´o la funci´on encargada de la comunicaci´on entre MATLAB y SU. Dicha funci´on fue denominada inversi´on autom´atica, est´a descrita en el anexo G. Cabe aclarar que este tipo de funci´on es creada en el entorno de programaci´on de MATLAB y es un complemento para los c´odigos de SA y ABC con la funci´on ﬁtness (Anexos H, I y J). El puente de comunicaci´on entre MATLAB y SU se realiza a trav´es de la funci´on inversi´on autom´atica, quien es la encargada de realizar los cambios en las l´ıneas de

p. 65

c´odigo en la terminal de UBUNTU necear´ıas para realizar los procesos de creaci´on de un modelo, su respectiva adquisici´on y su comparaci´on entre las trazas s´ısmicas. La caracter´ıstica del puente de comunicaci´on es de car´acter unidireccional desde MATLAB a SU y la forma en que se transforma en un puente bidireccional se hace en primera medida por medio de los comandos de operaciones matem´aticas que ofrece el paquete de SU. A partir de esto se procede a convertir la cadena de caracteres que son entregados por SU a un formato decimal de tal forma que MATLAB pueda leer el indicador de comparaci´on.

El diagrama de de interfaz entre MATLAB y SU se presenta en la Figura 18. Donde la m´etrica de comparaci´on y las t´ecnicas metaheur´ısticas son desarrolladas a trav´es de MATLAB. La generaci´on de modelos y el proceso de adquisici´on s´ısmica son realizados por medio de SU. El canal de comunicaci´on como ya se menciono anteriormente es la funci´on inversi´on autom´atica, esta funci´on es la encargada de realizar la actualizaci´on del modelo sint´etico, de generar el nuevo modelo con las capas actualizadas, de realizar una nueva adquisici´on s´ısmica, el procesamiento de los datos .SEGY y .SU, convertir las trazas s´ısmicas en formas de estructuras y obtener el coeﬁciente correlaci´on de las trazas s´ısmicas de referencia con las obtenidas en la nueva adquisici´on. Posteriormente, las t´ecnicas metaheur´ısticas realizan un proceso de minimizar el valor de la funci´on objetivo que es representado por el vector que la m´etrica de comparaci´on entrega en cada nueva iteraci´on.

p. 66

Figura 18. Diagrama de comunicaci´on entre MATLAB y SU.

Otra modiﬁcaci´on entre la comunicaci´on de MATLAB y SU se hace a causa de la manipulaci´on de las trazas s´ısmicas, las cuales est´an limitadas por los comandos que ofrece SU. Por esta raz´on, la manipulaci´on de las trazas s´ısmicas es implementada en MATLAB, de tal modo que se permita no solo recibir un indicador de comparaci´on de trazas por parte de SU, sino la posibilidad de implementar una m´etrica de comparaci´on que manipule las trazas s´ısmicas y arroje los datos m´as dicientes a SA y ABC. Finalmente, al llamar a la funci´on inversi´on autom´atica por el c´odigo principal, ella tiene la capacidad de hacer de manera autom´atica las siguientes funciones (Figura 18):

• Se actualiza el modelo sint´etico para la part´ıcula (SA) o una poblaci´on de

part´ıculas (ABC) de acuerdo a los sloths iniciales.

• Se genera el modelo del subsuelo con las capas actualizadas.

• Se simula la adquisici´on para obtener las trazas sint´eticas.

p. 67

• Se realiza la conversi´on del formato .SU a .SEGY y viceversa.

• Se importan las trazas s´ısmicas sint´eticas y reales a MATLAB en forma de

estructura.

• Se calculan los coeﬁcientes de correlaci´on entre las trazas s´ısmicas sint´eticas

y reales.

3.5.1.

SA y ABC implementados en el problema de inversi´on s´ısmica En el transcurso del seminario de investigaci´on surgieron cambios en la adaptaci´on con respecto al numero de dimensiones de los algoritmos de SA y ABC que fueron utilizados para la evaluaci´on de sus respectivos funcionamientos en funciones de prueba. Las modiﬁcaciones realizadas corresponden a la adaptaci´on de SA y ABC con la funci´on inversi´on autom´atica. La versi´on ﬁnal de SA y ABC utilizados en el proceso de inversi´on s´ısmica estan localizados en los anexos H y I. En primer lugar se realizaron las implementaciones de los algoritmo SA y ABC, que se encargaba de buscar un modelo sint´etico por medio de las coordenadas localizadas por la part´ıcula o el conjunto de part´ıculas, la cual se presenta como la analog´ıa a la b´usqueda del m´ınimo global de modo que presenten alguna coincidencia con los valores de sloths del modelo de referencia. Segundo, se implement´o la correlaci´on cruzada como m´etrica de comparaci´on. La m´etrica de comparaci´on entrega un vector, el cual esta constituido por los m´aximos coeﬁcientes de correlaci´on, que han sido obtenidos mediante la comparaci´on entre las trazas del modelo de referencia y las trazas del modelo sint´etico. Finalmente, se procedi´o a realizar se procedieron a realizar un n´umero de pruebas con SA y ABC, estas pruebas se muestran en el capitulo 4 y en el anexo K. En el desarrollo de las pruebas se efectuaron variaciones con relaci´on a los par´ametros de SA (probabilidades,

p. 68

temperaturas, iteraciones, ensayos, entre otros) y de ABC (poblaci´on, iteraciones y rutinas). En consecuencia, el objetivo de esta variaci´on es el de encontrar el mejor conjunto de par´ametros, que permitan hallar los sloths correspondientes a un modelo sint´etico que presente una soluci´on aproximada con relaci´on al modelo usado como referencia, el cual esta ubicado en el capitulo 2, en la secci´on 2.1 y en la Figura 8.

p. 69

PRUEBAS Y RESULTADOS

En la etapa ﬁnal del seminario se realizaron las pruebas y simulaciones con SA y ABC. Adem´as, se hizo la implementaci´on y emparejamiento de las t´ecnicas metaheur´ısticas con el proceso de inversi´on s´ısmica para encontrar el modelo de velocidades ´optimo del subsuelo. Por consiguiente, en este cap´ıtulo se habla de las simulaciones realizadas y los resultados obtenidos.

4.1

METODOLOG´ıA En el proceso de las simulaciones se abarcan aquellos par´ametros deﬁnitivos que fueron seleccionados a trav´es de las simulaciones realizadas con las funciones pruebas, aplicadas en SA y ABC. Los par´ametros son de crucial importancia como se observo en las Tablas 6 y 8 en el capitulo 3, donde se puede observar el comportamiento de los algoritmos sobre las funciones de prueba en el proceso de minimizaci´on. Teniendo en cuenta los resultados de los experimentos aplicados en las funciones de prueba (ver Anexos E y F) los par´ametros de funcionamiento de SA para dar comienzo son: una probabilidad inicial de 0,7, una probabilidad ﬁnal de 0,0001 y un n´umero de iteraciones variable. Para el caso de ABC, todo el algoritmo radica en la selecci´on del tama˜no de la poblaci´on, este n´umero es de 30 abejas. A diferencia de los trabajos realizados con ABC por el grupo CEMOS [11], [12], [13] donde el n´umero de la poblaci´on era de un valor superior en comparaci´on a nuestra selecci´on. Los valores utilizados en los trabajos de grado desarrollados por el grupo CEMOS fueron descartados a causa del funcionamiento que presento ABC sobre

p. 70

las funciones de prueba con una poblaci´on de 30 abejas (v´ease Tabla 8). El procedimiento de las simulaciones realizadas para SA y ABC fue el siguiente: de las cuatro capas del modelo principal creado a trav´es de SU los valores de sloths iniciales del modelo de referencia son 0,77; 0,50; 0,26 y 0,12, se tom´o como base a las pruebas que se realizaron, y por ´ultimo se llev´o el registro de cada uno de los respectivos casos (Anexo K).

El algoritmo que present´o las principales modiﬁcaciones en los valores de los par´ametros fue SA, donde el valor de 0,7 en su probabilidad inicial no presentaba resultados signiﬁcativos. Por consiguiente, este par´ametro fue cambiando como se pude observar en la Tabla 10 y la Tabla K 1 ubicada en el anexo K de complementos. Las variaciones de los valores iniciales de las capas al igual que sus respectivos puntos de inicio en el desarrollo de todas las simulaciones.

4.2

RESULTADOS

Los datos de SA arrojados por las pruebas realizadas, se presentan de la siguiente forma: para su interpretaci´on las pruebas relevantes que contienen los datos arrojados se encuentran ubicados en una tabla. La tabla contiene los par´ametros utilizados y que fueron variados, como lo son: la probabilidad inicial, la probabilidad ﬁnal o criterio de parada, la temperatura inicial, la temperatura ﬁnal y el paso de reducci´on en cada iteraci´on. La m´etrica utilizada fue la correlaci´on cruzada. Las especiﬁcaciones de la maquina de c´omputo utilizada para el desarrollo de las pruebas sobre el proceso de inversi´on s´ısmica se aprecian en la Tabla 9. Las pruebas que se

p. 71

presentan parten de un punto cercano para as´ı observar la convergencia que presenta el algoritmo de SA bajo el proceso de inversi´on s´ısmica. En lo que respecta a ABC, se realiza una prueba diferente, ya que aqu´ı no se tiene la b´usqueda de una part´ıcula alrededor de un campo, sino que es una colonia o grupo de puntos que son evaluados y que seguidamente buscan nuevas soluciones que dependen de un valor de rendimiento o ﬁtness.

Para la primera prueba desarrollada mediante SA se usaron los par´ametros contenidos en la Tabla 10. En esta primera prueba se tienen que cuatro de los cinco intentos realizados, donde cada intento constaba 5 iteraciones, la part´ıcula presento un movimiento cercano del valor exacto; que para el caso de la primera capa tiene el valor de 0,77. Con esto presente, se tiene que el valor de la correlaci´on no es el adecuado, ya que con cada iteraci´on presenta un alejamiento del valor de 1 el cual es el valor ideal. El alejamiento se puede evidenciar en la Figura 20, donde cada ge´ofono es utilizado sobre la adquisici´on s´ısmica y como en cada iteraci´on (Figura 19 realizada en el transcurso de la prueba para encontrar el valor de la primera capa. Tabla 9. Especiﬁcaciones del equipo de c´omputo usado en los experimentos sobre el proceso de inversi´on s´ısmica Procesador Intel R ⃝Xeon R ⃝CPU ES-2609 @ 2.40GHz Memoria RAM DDR3 de 8.00 GB Sistema operativo debian 8 de 64 bits

p. 72

Tabla 10. Par´ametros de la primera prueba con SA Par´ametros Prueba 1 Prueba 2 Intentos

Iteraciones

Probabilidad inicial

0,7

Probabilidad ﬁnal

0,001

Temperatura inicial

2,8037

Temperatura ﬁnal

0,1448

Paso de reducci´on

0,3

Tiempo de c´omputo [s]

30,11

58,85

En las Figuras 19 y 20 se observa la variaci´on del coeﬁciente de correlaci´on cruzada que se present´o en la prueba realizada con respecto al n´umero de iteraciones y g´eofonos. El valor ideal del coeﬁciente es igual a uno, pero en esta prueba se puede evidenciar que los par´ametros que se utilizaron con SA no son los indicados a causa de la lejan´ıa a los valores del coeﬁciente de correlaci´on cruzada y el valor de sloth en la primera capa.

Figura 19. Promedio de los coeﬁcientes de la correlaci´on cruzada con respecto al n´umero de iteraciones en la primera prueba con SA.

p. 73

Figura 20. Promedio de los coeﬁcientes de la correlaci´on cruzada con respecto al n´umero de ge´ofonos en la primera prueba con SA.

Para la segunda prueba se aumentaron las iteraciones con respecto a los par´ametros de funcionamiento de SA, obteniendo as´ı que tres de los cinco intentos culminan muy cercanos al valor de sloth real de la primera capa, quedando como un factor determinante a la hora de seleccionar par´ametros para futuros procesos de inversi´on. Esa cercan´ıa se ve entonces reﬂejada en algunos valores del comportamiento de la funci´on correlaci´on, pero aquellos que se encuentran lejos de la soluci´on intentan volver al valor real; no descartando entonces que ese factor exploratorio tiende a regresar despu´es de cierta cantidad de iteraciones. En las Figuras 21 y 22 se observan las variaciones tanto del valor de los sloths de la primera capa como los coeﬁcientes de correlaci´on cruzada que se presentaron en la prueba realizada con respecto al n´umero de iteraciones y g´eofonos. En el anexo K se encuentran ubicadas las pruebas complementarias con respecto a SA.

p. 74

Figura 21. Promedio de sloths en la segunda prueba con SA.

Figura 22. Promedio de los coeﬁcientes de la correlaci´on cruzada con respecto al n´umero de ge´ofonos en la segunda prueba con SA.

Con ABC se realiza una prueba diferente, ya que aqu´ı no se tiene la b´usqueda de una part´ıcula alrededor de un campo, sino que es una colonia o grupo de puntos que son evaluados y que seguidamente buscan nuevas soluciones basados en un vector de rendimiento (ﬁtness), o sea que el espacio de soluciones est´a delimitado pero en constante movimiento y selecci´on. Por algunas cuestiones de tiempo sobre el cronograma del seminario de investigaci´on, s´olo se realiza un seguimiento guiado (dando valores aleatorios a la primera capa, pero cercanos a la soluci´on de la misma) a la convergencia del algoritmo. Por ejemplo, a continuaci´on: en la Tabla 11, se presenta el cambio del alimento obtenido por la colonia para la primera capa mientras las dem´as se mantienen en un valor est´atico; en la Tabla 12, se presentan

p. 75

las soluciones exploradas por la colonia para la primera capa y en la Tabla 13, se muestra la mejor soluci´on encontrada por la colonia. La prueba fue realizada con una colonia de 30 abejas, donde aquella que tiene un valor m´as cercano al real es guardado para que aleatoriamente la colonia se mueva y as´ı se busquen soluciones, y si dado el caso no aparece un mejor valor se selecciona el m´as cercano como soluci´on. Esta informaci´on se observa en la quinta columna, donde se guarda ese primer valor de capa y uno de los dem´as es cambiado aleatoriamente, pero si el entregado como nuevo valor 0,8375, no es mejor que 0,7942, el que estaba previamente se deja como el mejor vector soluci´on encontrado. Por esto, la convergencia del algoritmo de ABC estar´ıa limitada y sujeta al tama˜no de la colonia y a la cantidad de abejas que se cambian en el proceso de b´usqueda de nuevos puntos aleatorios.

Tabla 11. Alimento obtenido por la colonia en la prueba realizada con ABC Alimento Valores Primera capa

0,7986

0,8021

0,7314

0,7377

0,7942

Segunda capa

0,5000

Tercera capa

0,2600

Cuarta capa

0,1200

Tabla 12. Soluciones exploradas por la colonia en la prueba realizada con ABC Soluciones exploradas Valores Primera Capa

0,8375

0,8021

0,7314

0,7377

0,7942

Segunda Capa

0,7000

Tercera Capa

0,7000

Cuarta Capa

0,7000

p. 76

Tabla 13. Soluciones exploradas por la colonia en la prueba realizada con ABC Sloths ﬁnales Valores Primera Capa

0,7942

Segunda Capa

0,5000

Tercera Capa

0,2600

Cuarta Capa

0,1200

En s´ıntesis, SA presenta un funcionamiento con la probabilidad inicial de 0,7 pero el algoritmo debe presentar un mayor numero de iteraciones para que el valor de los coeﬁcientes de cada ge´ofono obtengan una mayor cercan´ıa al valor ideal de 1. En lo que respecta a ABC, la poblaci´on de 30 abejas es un numero adecuado pero el algoritmo debe ser simulado con un mayor numero de iteraciones para poder apreciar con mayor detalle el trabajo realizado por toda la colonia. Finalmente es necesario contrastar los valores ﬁnales de los modelos de velocidades obtenidos mediante el uso de SA, ABC y el modelo sint´etico de referencia (tabla 14); aqu´ı se presentan los sloths ﬁnales de pruebas para cada algoritmo (valores tomados de la Figura 19 y la Tabla 13. Todo esto con el objetivo de dar indicios (as´ı se haya partido de la soluci´on y s´olo variara la primera capa) del comportamiento y evoluci´on de las t´ecnicas metaheur´ısticas durante el proceso de inversi´on s´ısmica.

Tabla 14. Comparaci´on de sloths ﬁnales entre SA y ABC con respecto al modelo de referencia Capas del modelo Sloth modelo de referencia SA

ABC

Primera Capa

0,77

0,775

0,794

Segunda Capa

0,50

Tercera Capa

0,26

Cuarta Capa

0,12

p. 77

CONCLUSIONES Y RECOMENDACIONES

En este cap´ıtulo se presentan las conclusiones y recomendaciones de lo desarrollado con las t´ecnicas metaheur´ısticas aplicadas en un problema de inversi´on s´ısmica, en el transcurso del seminario de investigaci´on.

5.1

CONCLUSIONES

Se encuentra que el comportamiento de las t´ecnicas metaheur´ısticas en la obtenci´on del modelo de velocidades del subsuelo, est´a directamente relacionado con los par´ametros de la probabilidad inicial de 0,7 para SA y una poblaci´on de 30 abejas para ABC sobre los cuales operan estas t´ecnicas, sobre todo en lo que concierne a la capacidad exploratoria de soluciones en estos algoritmos. Por ejemplo, durante la evaluaci´on de las funciones de prueba, una poblaci´on de

30 en una t´ecnica como ABC permite tener m´as informaci´on del espacio de b´usque-

da, lo cual aumenta las probabilidades de encontrar con mayor precisi´on una posible soluci´on. De igual manera, representa una reducci´on notable en el n´umero de iteraciones realizadas y en el tiempo de c´omputo.

Un valor de probabilidad inicial de 0,2 para SA y un valor de poblaci´on inferior a 10 abejas para ABC, representa un convergencia prematura como puede eviden-

p. 78

ciarse en los experimentos realizados sobre las funciones de prueba. La carencia de una estrategia en la realizaci´on de los experimentos, conlleva a una interpretaci´on err´onea de los resultados del algoritmo en donde se cree que el funcionamiento es err´oneo e incompleto como lo sucedido en las pruebas realizadas con SA y ABC en el cap´ıtulo 4.

Se entregan las herramientas te´orico-pr´acticas de todo lo desarrollado en el transcurso del seminario de investigaci´on, estas herramientas se presentan en memorias digitales que contienen los scripts de MATLAB de las dos t´ecnicas metahrur´ısticas, los modelos generados en SU, las memorias y actas de cada sesi´on.

5.2

RECOMENDACIONES

Para trabajos posteriores en la implementaci´on de los dos algoritmos de optimizaci´on global, ser´ıa de gran utilidad hacer un programa en MATLAB que realice una variaci´on de todos los par´ametros que contienen los algoritmos, con la idea de que autom´aticamente se puedan obtener resultados de todas las combinaciones posibles. Con base en esto, es posible obtener una base de datos que contenga informaci´on importante de las pruebas, en donde se permita hacer un estudio detallado del comportamiento de los algoritmos en relaci´on con los valores de los par´ametros, para determinar cu´al es la combinaci´on de par´ametros que arrojar´ıa resultados m´as precisos, en cuanto al hallazgo del m´ınimo global en funciones prueba de 2D y la

p. 79

mejor combinaci´on de sloths en lo referente a inversi´on s´ısmica. Reducir las l´ıneas de c´odigo de los algoritmos con la ﬁnalidad de optimizar recursos computacionales, por ejemplo, las l´ıneas de c´odigo que hacen la actualizaci´on de las ecuaciones de posici´on y ﬁtness del alimento en el algoritmo ABC, podr´an reducirse al asignar una matriz que represente en sus columnas y ﬁlas los valores necesarios para actualizar estas ecuaciones, quitando un ciclo de repetici´on que recorre cada una de las dimensiones para la respectiva actualizaci´on. En futuros trabajos relacionados con la implementaci´on de algoritmos en problemas de inversi´on s´ısmica, se podr´a ahondar en las tem´aticas requeridas para hacer que la combinaci´on de sloths esperada no corresponda a un modelo con geometr´ıa ﬁja como se trabaj´o durante este seminario, sino que por el contrario, generar modelos con mayor complejidad y con caracter´ısticas con mayor similitud a la superﬁcie terrestre.

La implementaci´on de m´etricas de comparaci´on de trazas s´ısmicas fueron enfocadas en la comparaci´on de se˜nales en toda su extensi´on, estas fueron, correlaci´on y correlaci´on cruzada de sismogramas. Sin embargo, no se utiliz´o una m´etrica que tuviera en cuenta la energ´ıa de la se˜nal dado que se plante´o en una de las sesiones la necesidad de considerar este par´ametro en vista de que el primer arribo de las trazas contiene una mayor energ´ıa que podr´ıa opacar el valor de las capas m´as profundas, por eso es recomendable buscar una que lo haga con el ﬁn de permitir una nueva alternativa de comparaci´on.

p. 80

Hacer pruebas en las que se combinen las m´etricas de comparaci´on, ya que es posible que en la ejecuci´on de los algoritmos se requiera la utilizaci´on de varias m´etricas como es el caso de la correlaci´on cruzada y la diferencia de m´ınimos cuadrados, que permitan una mejor comparaci´on en los diferentes momentos de la implementaci´on.

Realizar el proceso de una adquisici´on sint´etica a partir de la ecuaci´on de onda y no de trazado de rayos como se realiz´o. Asimismo, aumentar el n´umero de fuentes y ge´ofonos en el proceso de adquisici´on en SU para obtener una mejor resoluci´on en las trazas s´ısmicas.

p. 81

[1] CAICEDO, Mario y MORA Pl´acido. Temas de propagaci´on de ondas. Universidad Sim´on Bol´ıvar, 2004.

[2] P ´EREZ, Carlos, MATEO M´onica y MACI ´A Antonio. Aplicaci´on de tomograf´ıa de refracci´on s´ısmica y an´alisis de microtemores como t´ecnicas de prospecci´on geof´ısicas en estudios geot´ecnicos en ediﬁcaci´on. En: Informes de la Construcci´on [en l´ınea] vol.65, abril-junio, 2013. http://rua.ua.es/dspace/bitstream/ 10045/33408/1/2013_Perez_Mateo_Macia_InformesConstr.pdf Revisado el 20 de Marzo de 2014.

[3] FOREL, David, BENZ Thomas y PENNINGTON Waybe. Seismic data processing with Seismic Un*x. Society of Exploration Geophysicists, 2005. [4] KIRKPATRICK, Scott, GELATT Daniel y VECCHI Mario. Optimization by Simulated Annealing. American Association for the Advancement of Scienci, 1983, vol.220.

[5] KARABOGA, Dervis. An Idea Based on Honey Bee Swarm for Numerical Optimization. Technical Report, 2005, vol.6.

[6] KENNEDY, James y EBERHART Russell. Particle Swarm Optimization. Purdue School of Engineering and Technology, Indianapolis, 1995, vol.95,. p. 1942-

1948.

[7] MITCHELL, Melanie. An introduction to genetic algorithms. Massachusetts Institute of Technology, 1998.

[8] SHAA-HOSSEINI, Hamed. An Approach to Continuos Optimization by the Intelligent Water Drops Algorithm. Procedia - Social and Behavioral Sciences, 2001, vol.32,. p. 224-229.

REFERENCIAS BIBLIOGRÁFICAS

p. 82

[9] MRINAL, Sen y STOFFA Paul. Global optimization methods in geophysical inversion. Advances in exploration geophysics. Austin, Texas. Elsevier Science,

1995.

[10] ROMERO, Jorge y SUAREZ John. Soluci´on de las ecuaciones que describen el comportamiento de los modos h´ıbridos de una gu´ıa de onda Rectangular parcialmente llena mediante el m´etodo de Optimizaci´on Recocido Simulado (Simulated Annealing). Bucaramanga.: Universidad Industrial de Santander,

2013.

[11] CELIS, Julieth y RINC ´ON Francis. Evaluaci´on y comparaci´on entre los m´etodos Newton Raphson y Artiﬁcial Bee Colony (ABC) para el an´alisis del ﬂujo de carga de un sistema de potencia. Bucaramanga.: Universidad Industrial de Santander,

2013.

[12] FUENTE, Rafael y PETRO Elkin. Mantenimiento preventivo b´asico de un desﬁbrilador monof´asico mediante los m´etodos de Enjambre de Part´ıculas Mejorado y Colonia Artiﬁcial de Abeja. Bucaramanga.: Universidad Industrial de Santander, 2013.

[13] ´AVILA, Jos´e y NAVARRO Orlando. El m´etodo de Colonia Artiﬁcial de Abejas y el criterio de m´ınima entrop´ıa para el dise˜no ´optimo de un disipador de calor. Bucaramanga.: Universidad Industrial de Santander, 2014. [14] REZA, Seyyed, MALEKI Isa, HOJJATKHAH Sohrab y BAGHERINIA Ali. Evaluation the efﬁciency of artiﬁcial bee colony and the ﬁreﬂy algorithm in solving the continuos optimization problem. En: International Journal on Computational & Applications (IJCSA) [en l´ınea] vol.3, agosto, 2013. http:// arxiv.org/ftp/arxiv/papers/1310/1310.7961.pdf Revisado el 20 de Enero de 2014.

p. 83

[15] Bioinformatics Laboratory. Test functions for optimization needs [en l´ınea] http://www.bioinformaticslaboratory.nl/twikidata/pub/Education/ NBICResearchSchool/Optimization/VanKampen/BackgroundInformation/ TestFunctions-Optimization.pdf Revisado el 1 de Marzo de 2014. [16] AKAY, Bahriye and KARABOGA Dervis. A Modiﬁed Artiﬁcial Bee Colony Algorithm for Real-parameter Optimization. Information Sciences, 2012, vol.192, p.

120-142.

[17] Intelligent Systems Research Group, Departament of Computer Engineering, Erciyes University, Turkiye. Artiﬁcial Bee Colony (ABC) Algorithm Homepage [en l´ınea] http://mf.erciyes.edu.tr/abc/index.htm Revisado el 10 de Mayo de 2014.

[18] Wikipedia, La enciclopedia libre. Geof´ısica [en l´ınea] http://es.wikipedia. org/wiki/Geofisica Revisado el 2 de Agosto de 2014.

[19] Ecopetrol. Proyecto exploratorio Siriri-Catleya. Adquisici´on datos s´ısmicos [en l´ınea] http://www.ecopetrol.com.co/especiales/siriri/docs/0042.pdf Revisado el 8 de Mayo de 2014.

[20] RUIZ, Cristina. Inversi´on s´ısmica y estudio de atributos s´ısmicos post apilamiento de los niveles i3 y tu de la formaci´on oﬁcina en el campo guico guara, Estado Anzoategui. Universidad Sim´on Bol´ıvar, 2007. [21] F´ısica en L´ınea por Elba Sep´ulveda. Ley de Hooke [en l´ınea] https://sites. google.com/site/timesolar/fuerza/ley-de-hooke Revisado el 18 de Agosto de 2014.

[22] Schlumberger Oilﬁel Glosary. Seismic velocity [en l´ınea] http://www.glossary. oilfield.slb.com/es/Terms/sseismic_velocity.aspx Revisado el 18 de Agosto de 2014.

p. 84

[23] Laboratorio de Procesado de Imagen. Ondas [en l´ınea] http://www.lpi.tel. uva.es/~nacho/docencia/ing_ond_1/trabajos_06_07/io3/public_html/ Ondas/Ondas.html Revisado el 20 de Agosto de 2014.

[24] Schlumberger Oilﬁel Glosary. Impedancia el´astica [en l´ınea] http://www. glossary.oilfield.slb.com/es/Terms/e/elastic_impedance.aspx Revisado el 18 de Agosto de 2014.

[25] Schlumberger Oilﬁel Glosary. Impedancia ac´ustica [en l´ınea] http://www. glossary.oilfield.slb.com/es/Terms/a/acoustic_impedance.aspx Revisado el 18 de Agosto de 2014.

[26] Schlumberger Oilﬁel Glosary. ´Angulo de incidencia [en l´ınea] http://www. glossary.oilfield.slb.com/es/Terms/a/angle_of_incidence.aspx Revisado el 18 de Agosto de 2014.

[27] Universitat Polit`ecnica de Catalunya. Procesado de S´ısmica de Reﬂexi´on Supericial [en l´ınea] http://upcommons.upc.edu/pfc/bitstream/2099.1/3404/ 7/41205-7.pdf Revisado el 10 de Marzo de 2014.

[28] Schlumberger. Inversi´on s´ısmica: Lectura entre l´ıneas. En: Oilﬁeld Review [en l´ınea] vol.20, enero, 2008 http://www.slb.com/~/media/Files/resources/ oilfield_review/spanish08/sum08/inversion_sismica.pdf Revisado el 2 de Julio de 2014.

[29] GARC´IA, V´ıctor. Aplicaci´on de un algoritmo de inversi´on s´ısmica bayesiana pre-apilamiento para estimaci´on de propiedades el´asticas en un yacimiento gas´ıfero costa afuera, Trinidad & Tobago. Universidad Sim´on Bol´ıvar, 2006. [30] MONCAYO, Edward. Inversi´on s´ısmica mediante un algoritmo gen´etico/Seismic inversion with a genetic algorithm. Tesis Doctoral. Universidad Nacional de Colombia, 2010.

p. 85

[31] LARA, Jos´e y CAICEDO Mario. Modelado s´ısmico con seismic unix. Universidad Sim´on Bol´ıvar, 2010.

[32] HALE, Dave y COHEN Jack. Triangulated models of earth?s subsurface. School of Mines, Center of Wave Phenomena, Colorado, 1991.

p. 86

BIBLIOGRAF´IA

• ´AVILA, Jos´e y NAVARRO Orlando. El m´etodo de Colonia Artiﬁcial de Abejas y

el criterio de m´ınima entrop´ıa para el dise˜no ´optimo de un disipador de calor. Bucaramanga.: Universidad Industrial de Santander, 2014.

• CELIS, Julieth y RINC ´ON Francis. Evaluaci´on y comparaci´on entre los

m´etodos Newton Raphson y Artiﬁcial Bee Colony (ABC) para el an´alisis del ﬂujo de carga de un sistema de potencia. Bucaramanga.: Universidad Industrial de Santander, 2013.

• FOREL, David, BENZ Thomas y PENNINGTON Waybe. Seismic data

processing with Seismic Un*x. Society of Exploration Geophysicists, 2005.

• FUENTE, Rafael y PETRO Elkin. Mantenimiento preventivo b´asico de un

desﬁbrilador monof´asico mediante los m´etodos de Enjambre de Part´ıculas Mejorado y Colonia Artiﬁcial de Abeja. Bucaramanga.: Universidad Industrial de Santander, 2013.

• KARABOGA, Dervis. An Idea Based on Honey Bee Swarm for Numerical

Optimization. Technical Report, 2005, vol.6.

• KIRKPATRICK, Scott, GELATT Daniel y VECCHI Mario. Optimization by

Simulated Annealing. American Association for the Advancement of Scienci, 1983, vol.220.

• MRINAL, Sen y STOFFA Paul. Global optimization methods in geophysical

inversion. Advances in exploration geophysics. Austin, Texas. Elsevier Science,

1995.

p. 87

• ROMERO, Jorge y SUAREZ John. Soluci´on de las ecuaciones que describen

el comportamiento de los modos h´ıbridos de una gu´ıa de onda Rectangular parcialmente llena mediante el m´etodo de Optimizaci´on Recocido Simulado (Simulated Annealing). Bucaramanga.: Universidad Industrial de Santander,

2013.

p. 88

ANEXOS

ANEXO A. CONCEPTOS B ´ASICOS

A.1. GEOF´ISICA

La geof´ısica es la ciencia encargada del estudio de la estructura, condiciones f´ısicas y evoluci´on de la corteza terrestre. Para su estudio, usa m´etodos cuantitativos f´ısicos naturales o dise˜nados por el hombre, como son la reﬂexi´on y refracci´on de ondas mec´anicas, la gravedad, campos electromagn´eticos, magn´eticos o el´ectricos, fen´omenos f´ısicos y radiactivos. La geof´ısica abarca dos grandes ramas: la geof´ısica interna y externa [18].

A.1.1. Geof´ısica interna Es la encargada de analizar el interior de la Tierra y los principales temas que estudia son [18]:

• Sismolog´ıa: Su objetivo de estudiar los fen´omenos generados al interior de la

Tierra producto de la propagaci´on de ondas el´asticas.

• Geodin´amica: Se encarga de estudiar aquellas modiﬁcaciones que sufre la

corteza terrestre.

A.1.2.

Geof´ısica externa Estudia las propiedades f´ısicas del entorno terrestre, entre ellas se encuentran el geomagnetismo, paleomagnetismo, gravimetr´ıa,oceanograf´ıa, metereolog´ıa, areonom´ıa y climatolog´ıa, pero debido a que ese campo de la geof´ısica no es nuestro objeto de estudio, no se abarcar´a en este documento [18].

p. 89

A.2. S´ISMICA

Es la generaci´on artiﬁcial de ondas las cuales penetran el interior de la de la Tierra, atravesando las diferentes capas geol´ogicas, produciendo peque˜nos ecos que son detectados en la superﬁcie por instrumentos altamente sensibles.La informaci´on captada es procesada con herramientas de computaci´on e interpretada geol´ogicamente para entender la conﬁguraci´on y geometr´ıa interna de la corteza terrestre [19].

La s´ısmica dentro de la industria petrolera tiene como principal objetivo localizar aquellas trampas que contengan hidrocarburo. Existen ciertos ambientes que tienen estructuras que brindan las condiciones necesarias para encontrar excelentes trampas, sin embargo ´estas necesariamente no contienen hidrocarburo, por esta raz´on se trata de extraer la mayor cantidad de informaci´on posible de los datos s´ısmicos, para tratar de descifrar las propiedades del ´area de estudio [20]. A.2.1. Teor´ıa S´ısmica: Conceptos B´asicos

• Ley de la elasticidad de Hooke: Establece la relaci´on entre el alargamiento o

estiramiento longitudinal y la fuerza aplicada, siendo entonces, una propiedad f´ısica en la que los objetos con capaces de cambiar de forma cuando act´ua una fuerza de deformaci´on sobre ellos [21].

• Velocidades s´ısmicas: Es la velocidad con la que viaja una onda ac´ustica

a trav´es de un medio, es decir, distancia dividida por el tiempo de viaje; esta puede determinarse a partir de perﬁles s´ısmicos verticales o a partir del an´alisis de velocidad de los datos s´ısmicos. Adem´as se presentan variaciones en sentido vertical, lateral y azimutal, en los medios anisotr´opicos como las rocas; las cuales tienden a incrementarse con la profundidad en la Tierra porque la compactaci´on reduce la porosidad [22].

p. 90

En un terremoto se transmiten ondas que viajan por el interior de la tierra, siguiendo caminos curvos debido a la variada densidad y composici´on del interior de la Tierra. A este tipo de ondas se llaman ondas internas, centrales o de cuerpo, transmitidas por los temblores preliminares de un terremoto pero poseen poco poder destructivo. Este clase de ondas son divididas en dos grupos [23]:

• Ondas P: Las ondas P (Primarias) son ondas longitudinales, lo cual signiﬁca

que el suelo es alternadamente comprimido y dilatado en la direcci´on de la propagaci´on (Figura A 1).

Figura A 1. Ondas P.

Fuente: Laboratorio de Procesado de Imagen [23].

• Ondas S: Las ondas S (secundarias) son ondas transversales o de corte, es

decir, el suelo es desplazado perpendicularmente a la direcci´on de propagaci´on, alternadamente hacia un lado y hacia el otro (Figura A 2). Las ondas S pueden viajar ´unicamente a trav´es de s´olidos y usualmente la onda S tiene mayor amplitud que la onda P.

Figura A 2. Ondas S.

Fuente: Laboratorio de Procesado de Imagen [23].

p. 91

A.2.2. Impedancia el´astica e impedancia ac´ustica La impedancia el´astica es el producto entre la densidad de un medio y la velocidad de su onda de corte (onda S) [24]. En cambio la impedancia ac´ustica es el producto de la densidad por la velocidad s´ısmica, ´esta var´ıa entre las diferentes capas de rocas y se indica generalmente con el s´ımbolo Z; teniendo en cuenta que las diferencias de impedancia ac´ustica entre las capas de rocas afecta el coeﬁciente de reﬂexi´on [25]. A.2.3. ´Angulo de Incidencia Es el ´angulo agudo en el que una trayectoria s´ısmica choca con una l´ınea normal a una interface; en geof´ısica, es similar a una onda s´ısmica que choca con los estratos. La ley de Snell describe la relaci´on entre el ´angulo de incidencia y el ´angulo de refracci´on de una onda, teniendo entonces que el ´angulo de refracci´on depende de la velocidad de la onda en ese medio [26]. A.2.4. Coeﬁciente de reﬂexi´on El coeﬁciente de reﬂexi´on describe la amplitud (o la intensidad) de una onda reﬂejada respecto a la onda incidente. Se supone que la onda incidente tiene una magnitud de uno, la reﬂejada R y la transmitida 1-R. El coeﬁciente de Reﬂexi´on (R) es una funci´on de las velocidades y las densidades de dos medios adyacentes a una interfaz [20].

A.2.5. S´ısmica de refracci´on La s´ısmica de refracci´on realiz´o un gran aporte a la prospecci´on (exploraci´on del subsuelo encaminada a descubrir yacimientos minerales, petrol´ıferos, arqueol´ogicos) s´ısmica en sus comienzos. Hasta la d´ecada de los

60 fue extremadamente popular, especialmente en la exploraci´on de cuencas sedi-

mentarias4 donde condujo al descubrimiento de grandes campos de petr´oleo. 4Cuenca sedimetaria: Es una zona deprimida presente en la corteza terrestre que se caracteriza por la acumulaci´on de sedimentos.

p. 92

El m´etodo se basa principalmente en la medici´on del tiempo de viaje de las ondas refractadas cr´ıticamente en las interfaces del subsuelo con diferentes propiedades f´ısicas; fundamentalmente por el contraste entre las impedancias ac´usticas (Figura A 3). La s´ısmica de refracci´on solo considera las refracciones con ´angulo cr´ıtico ya que son las ´unicas ondas refractadas que llegan a la superﬁcie y pueden ser captadas por los ge´ofonos [27].

Figura A 3. S´ısmica de refracci´on.

Fuente: Universitat Polit`ecnica de Catalunya. Procesado de S´ısmica de Reﬂexi´on Supericial [27].

Debido a su menor costo y al tipo de informaci´on que proporciona la s´ısmica de refracci´on, es un potente m´etodo que actualmente se emplea en estudios de estructuras profundas.

A.2.6. S´ısmica de reﬂexi´on Este m´etodo se basa en las reﬂexiones del frente de onda s´ısmica sobre las distintas interfaces del subsuelo, responden al igual que en la refracci´on, a contrastes de impedancias que posteriormente se relacionaran con las distintas capas del subsuelo. Las reﬂexiones son detectadas por los receptores (ge´ofonos) que se ubican en la superﬁcie y que est´an alineados con la fuente emisora (Figura A 4).

p. 93

Figura A 4. S´ısmica de reﬂexi´on.

Fuente: Universitat Polit`ecnica de Catalunya. Procesado de S´ısmica de Reﬂexi´on Supericial [27].

Este m´etodo es una de las t´ecnicas de prospecci´on geof´ısica mas utilizada, debido a que su resultado es una imagen denominada secci´on s´ısmica en donde se aprecia la geometr´ıa de las estructuras geol´ogicas [27].

A.3. INVERSI ´ON S´ISMICA

Para la industria de Energ´ıa y Petroleo (E & P), muchas mediciones son basadas en procesos de inversi´on para su interpretaci´on. Entonces, queriendo llegar a ella se necesitan estimaciones de una serie de resultados, pero ese ´ambito es muy inestable ya que se necesitan mediciones m´ultiples (las cuales se realizan a condiciones diferentes, con diferentes operarios e incluso con recursos cient´ıﬁcos diferentes) las cuales arrojan una relaci´on matem´atica que no tendr´a una ´unica soluci´on. Por esto, la inversi´on s´ısmica es una forma matem´atica que estima una respuesta basada en informaci´on previa y da la oportunidad de ser modiﬁcada hasta que sea aceptable.

El proceso de inversi´on, como su nombre lo indica, puede ser considerado como la inversa del modelo directo, al que a veces se reﬁere simplemente como modelado. A.3.1 Modelado directo El modelo directo comienza con un modelo de las propiedades del subsuelo, luego simula matem´aticamente un experimento o un

p. 94

proceso f´ısico por ejemplo, electromagn´etico, ac´ustico, nuclear, qu´ımico u ´optico en el modelo del subsuelo y ﬁnalmente provee como salida una respuesta modelada. Si el modelo y los supuestos son precisos, la respuesta modelada se asemeja a los datos reales.

A.3.2 Modelado inverso El modelado inverso hace lo contrario, o sea que comienza con los datos medidos reales, aplica una operaci´on que retrocede el proceso a trav´es de un experimento f´ısico y produce un modelo del subsuelo. Si la inversi´on se realiza correctamente, el modelo del subsuelo es semejante al subsuelo real. El proceso de inversi´on es utilizado por muchas disciplinas de Energ´ıa & Petroleo y puede aplicarse en una amplia gama de escalas y con niveles de complejidad variables. Entre ellas se encuentran [28]:

• Extracci´on de las litolog´ıas de las capas y las saturaciones de ﬂuidos a partir

de mediciones de registros m´ultiples.

• Interpretaci´on de vol´umenes de gas, petr´oleo y agua utilizando registros de

producci´on.

• Inferencia de la permeabilidad y los l´ımites del yacimiento derivados de los

datos de presiones transitorias.

• Mapeo de los frentes de ﬂuidos a partir de mediciones electromagn´eticas entre

pozos.

A.4. GENERALIDADES

El objetivo principal en la inversi´on s´ısmica es tratar de obtener un modelo de impedancia del subsuelo a partir de la combinaci´on de datos s´ısmicos e informaci´on de pozo. De esta forma, se crea un modelo cuantitativo del

p. 95

reservorio que nos permite caracterizar el mismo y establecer criterios ´optimos de gerencia y explotaci´on con el menor riesgo posible.

A.4.1. Teor´ıa del Problema Inverso En geof´ısica, uno de los principales objetivos es hacer descripciones cuantitativas acerca del subsuelo, entonces, la ﬁnalidad principal de la inversi´on es transformar las observaciones s´ısmicas en propiedades cuantitativas de rocas que describan un reservorio.

Matem´aticamente, un modelo es un espacio que se construye por medio de parametrizaciones de diferentes propiedades del subsuelo. Este espacio es llamado espacio del modelo y se denota como “m”. De igual forma se tiene un sistema de datos obtenidos experimentalmente, los cuales, por medio del modelado directo, se obtienen modelos del subsuelo.

Los datos observados correspondientes residen en un espacio matem´atico llamado espacio de los datos y se denota como “d”(Figura A 5). Figura A 5. Diagrama convencional del problema inverso.

p. 96

En general, los modelos creados son funciones continuas del espacio con inﬁnitos grados de libertad que describen el subsuelo. Paralelo a esto, los datos son funciones discretas ya que estos representan valores obtenidos experimentalmente. Por lo tanto, la inversi´on se enfoca hacia la resoluci´on de un problema de evaluaci´on del modelo estimado en base del modelo real (Figura A 6).

Una r´apida comparaci´on de ambas variables muestra que el paso de los datos a un modelo no puede ser ´unico; es decir, que deben existir elementos del espacio del modelo que no tienen inﬂuencia alguna en el espacio de los datos. Esta carencia de unicidad representa la idea principal en la soluci´on del problema inverso. Teniendo entonces que, se puede obtener un n´umero inﬁnito de modelos que pueden reﬂejar respuestas semejantes. El proceso de inversi´on entonces, se basa principalmente en la resoluci´on de un problema de estimaci´on, partiendo de un conjunto de datos observados; el objetivo es estimar el modelo que pueda simular estos datos con el menor error posible en comparaci´on con el modelo real del subsuelo. Figura A 6. Esquema del problema inverso.

Una forma de controlar la no unicidad en las soluciones en un proceso de inversi´on puede ser restringir los modelos con distribuciones a priori, las cuales incluyen

p. 97

informaci´on adicional a los datos observados d y se basan tanto en conocimientos previos del interprete, como en datos geol´ogicos y macro modelos construidos a partir de correlaciones con registros de pozos.

Finalmente el objetivo es minimizar el error entre los datos y los datos creados a partir de un modelo m. La funci´on que describe este error es llamada funci´on objetivo. Por lo tanto, el modelo ﬁnal ser´a aquel en donde la funci´on objetivo asociada “ond´ıcula - modelo - datos s´ısmicos”sea m´ınima. Esta minimizaci´on representa el prop´osito principal del proceso de optimizaci´on en la inversi´on s´ısmica [29].

A.4.2. Inversi´on basada en la traza Con la inversi´on basada en la traza, el proceso parte esencialmente de los datos s´ısmicos con la probabilidad de usar algunos datos no - s´ısmicos como por ejemplo, la tendencia de las frecuencias bajas derivadas de las velocidades de los pozos. Entonces se requiere un conjunto de m´etodos que usen las trazas s´ısmicas para calcular las frecuencias altas de la serie de reﬂectividad, ese es el objetivo principal de este m´etodo. [30]. A.4.3. Inversi´on basada en el modelo En este tipo de inversi´on, el proceso comienza con un modelo y se da un peso adicional a los datos no s´ısmicos. Estos datos no est´an limitados a la informaci´on de pozo, por ejemplo, las distribuciones estad´ısticas tambi´en pueden ser consideradas en el modelo inicial. El t´ermino “inversi´on basada en un modelo”hace referencia a un conjunto de m´etodos que intentan obtener como resultado una resoluci´on igual o mayor que la obtenida en la s´ısmica, dando un peso importante a informaci´on a priori que la ma-

p. 98

yor´ıa de las veces consiste en registros de pozo.

El valor de la impedancia depende del valor de la reﬂectividad y del valor de la impedancia de la capa previa, por lo que peque˜nos errores crean error acumulado en la inversi´on. Adicionalmente, debido a la no unicidad de la soluci´on, diferentes modelos pueden producir una traza sint´etica que sea muy similar a la traza s´ısmica, es decir diferentes modelos de impedancias pueden obtener un error m´ınimo sin ser el modelo m´as cercano a la realidad [30].

A.5. OPTIMIZACI ´ON GLOBAL

Es un procedimiento mediante el cual se determina el mejor valor o conjunto de valores de una funci´on o procesos en particular. Por lo general, la optimizaci´on busca minimizar o maximizar funciones por medio de la obtenci´on del punto correspondiente a dicha caracter´ıstica.

p. 99

ANEXO B. INTRODUCCI ´ON A SU

B.1. SU: PROCESAMIENTO Y MODELADO S´ISMICO

Para ﬁnales de la d´ecada de los 80?s, espec´ıﬁcamente para 1987, en el Centro de Fen´omenos Ondulatorios de la Escuela de Minas de Colorado, se estaba formando un grupo de trabajo, liderado por Jack Cohen y Shuki Ronen, destinado a la generaci´on de un conjunto de herramientas computacionales para las tareas relativas al procesamiento s´ısmico, todas estas herramientas pensadas para trabajar en un ambiente tipo UNIX. Es as´ı como nace SU, el cual es un paquete de tipo software libre que trabaja como extensi´on del shell de UNIX (o cualquier derivaci´on de ´este) [31].

B.2. PRINCIPIOS EN LINUX Y SHELL-SCRIPTING

En el capitulo 2 se mencion´o que el programa que utilizaremos en el seminario para el modelado ser´a SU. Todas las aplicaciones de SU corren en el sistema operativo UNIX, lo que hace necesario manejar los fundamentos b´asicos de redirecciones y tuber´ıas que se emplean en el manejo del shell (terminal) de UNIX.

Linux Es un sistema operativo, mantenido por miles de programadores del todo el mundo que, provee al usuario de una interfaz para interactuar de alg´un modo con el software (las aplicaciones que corren a nivel de usuario) y el hardware (dispositivos f´ısicos) de la m´aquina. Linux es s´olo el kernel (n´ucleo) de la maquina, que permite la interacci´on entre el hardware y el software, y fue creado por Linus Torvalds. Formalmente, todas las distribuciones que existen actualmente del sistema operativo Linux deber´ıan ser llamadas GNU/Linux; sin embargo, se ha popularizado

p. 100

el nombre del segundo miembro de esta uni´on [31].

Linux Es un sistema operativo de open source (c´odigo abierto), lo que implica que su c´odigo fuente est´a disponible para sus usuarios, y puede ser modiﬁcado o alterado a criterio de los mismos para mejorar, crear u optimizar ciertas y determinadas tareas que el sistema requiera. El hecho que Linux sea open source ha incentivado a la formaci´on de una comunidad mundial de programadores [31]. El terminal de Linux es quien realiza la interpretacion de los comandos que permite la comunicaci´on con el kernel de la m´aquina e interactuar con los dispositivos l´ogicos que controlan el hardware.. B´asicamente lo que se puede hacer a trav´es del entorno gr´aﬁco, se puede hacer por el terminal. Adem´as, una vez aprendido los comandos b´asicos y al ganar cierta agilidad realizando tareas rutinarias en el terminal, se hace mas eﬁciente el trabajo por medio de un entorno gr´aﬁco [31]. La programaci´on de shell-scripts es fundamental a la hora de trabajar de manera efectiva en el terminal. Los llamados scripts (guiones) son simplemente archivos de texto, que contienen un n´umero determinado de comandos (enlazados o no) que son ejecutados al momento de ser interpretados por el shell [31]. El lenguaje de shell-scripting m´as ampliamente usado es Bash, el cual viene del acr´onimo Bourne-Again-Shell. En la mayor´ıa de las derivaciones de UNIX, Bash es el lenguaje por defecto del shell. El uso de shell-scripting nos facilita de una manera mas r´apida de resolver ciertos problemas computacionales [31].

p. 101

Para trabajar con SU los scripts son esenciales para aprovechar al m´aximo el alcance del paquete. Mediante los scripts son realizadas las mayorias de las tareas importantes que se pueden hacer con SU, sobre todo en aquellas donde trabajar directamente en el shell se puede convertir en algo poco funcional e inefectivo [31].

B.3. CREACI ´ON DE MODELOS DEL SUBSUELO

La necesidad de los humanos de apreciar y captar la informaci´on a trav´es de sus receptores como lo son los sentidos (vista, o´ıdo, olfato, tacto y gusto). Uno de los receptores de informaci´on m´as importantes es la vista. A trav´es de la vista se ubican, se asocian, se interpreta y se entiende gran parte de las tareas, sucesos y dem´as escenarios en los que nos desenvolvemos. La dependencia que tenemos sobre las computadoras, ha obligado a que se desarrollen maneras eﬁcientes, y computacionalmente viables, de representar gr´aﬁcamente objetos, cuya interpretaci´on visual es importante al momento de la resoluci´on de un problema dado [31].

B.3.1 Triangulaci´on de Delaunay Existen muchas herramientas de visualizaci´on que cumplen con este ´ultimo prop´osito, entre ellas tenemos se puede mencionar la Triangulaci´on de Delaunay. La triangulaci´on de Delaunay es un m´etodo ampliamente utilizado en la computaci´on para la generaci´on de gr´aﬁcos (Figura B 1), esta no es mas que la uni´on de un conjunto de puntos a trav´es de tri´angulos. Sin embargo, estos tri´angulos deben cumplir una condiciona especial: la condici´on de Delaunay [32].

p. 102

Figura B 1. Triangulaci´on de la fotograf´ıa de una iglesia. Fuente: Inversi´on s´ısmica mediante un algoritmo gen´etico/Seismic inversion with a genetic algorithm [31].

Esta condici´on dice: el interior de la circunferencia que circunscribe a cada uno de los tri´angulos debe ser vac´ıa. La triangulaci´on es ´unica y si solo 3 v´ertices entran en la circunferencia que circunscriben cualquier triangulo. Esta formalidad es insigniﬁcante en aplicaciones pr´acticas [32].

En la geof´ısica la interpretaci´on visual juega un papel de muy importante y las acertadas consideraciones que pueden o no tomarse en cuenta a la hora del estudio pertinente. En el modelado s´ısmico, espec´ıﬁcamente en el modelado de los perﬁles de velocidades, no s´olo la visualizaci´on es de gran utilidad, sino tambi´en el adecuado ajuste de las propiedades f´ısicas respectivas a cada coordenada del modelo. Que sea necesario que se represente de la manera m´as ﬁel posible, las condiciones reales del subsuelo para as´ı obtener una respuesta s´ısmica comparable con la respuesta s´ısmica real [31].

En la triangulaci´on de Delaunay simple se puede observar que los tri´angulos tienen a ser m´as equil´ateros que en la triangulaci´on ajustada (Figura B 2), adem´as que en esta ´ultima se puede cuestionar si los v´ertices m´as cercanos a los bordes del cuadrado cumplen o no efectivamente con la condici´on de Delaunay [31].

p. 103

Figura B 2. Triangulaci´on de Delaunay simple (izquierda) y ajustada (derecha). Fuente: Triangulated models of earth’s subsurface [32].

Por medio de este m´etodo se pueden representar, de manera acertada estructuras del subsuelo simples y/o complejas. Sin embargo, no es la ´unica ventaja que tiene el m´etodo. Una ventaja que presentan los moldeos triangulares es la simplicidad que ofrecen a la hora de calcular tiempos de viaje y caminos de rayos cuando se hacen por medio del trazado de rayos [32].

p. 104

ANEXO C. C ´ODIGO PARA GENERAR UN MODELO EN SU

En las siguientes lineas se presenta el c´odigo utilizado para la generaci´on de los modelos del subsuelo en SU:

1.#! /bin/sh 2.# Modelo de pr´actica: Sesion 3. Seminario de Invetigaci´on.

4.# model number 5.model=1

7.# data directory (optional, if not set data will go into current directory) 8.psﬁle=model$model.eps 9.dataﬁle=model$model.dat

10.

11.trimodel xmin=0 zmin=0 xmax=10.0 zmax=7.0

12.

1 xedge=0.0,10.0

13.

zedge=0.0,0.0

14.

sedge=0,0

15.

2 xedge=0.0,3.0,5.5,9.5,10.0

16.

zedge=3.0,2.8,2.0,3.0,3.0

17.

sedge=0,0,0,0,0

18.

3 xedge=0.00,3.0,7.0,10.00

19.

zedge=5.00,4.9,5.0,6.00

20.

sedge=0,0,0,0

21.

4 xedge=3.00,7.00

22.

zedge=4.90,7.00

p. 105

23.

sedge=0,0

24.

5 xedge=0.00,7.0,10.00

25.

zedge=7.00,7.0,7.00

26.

sedge=0,0,0 27.kedge=1,2,3,4,5 28.v1 sﬁll=0.01,1.01,0.0,0,4.02,0.0,0.0

29.

sﬁll=0.01,4.00,0.0,0,3.70,0.0,0.0

30.

sﬁll=5.00,5.50,0.0,0,2.00,0.0,0.0

31.

sﬁll=0.10,6.00,0.0,0,1.40,0.0,0.0 32.>$dataﬁle 33.spsplot < $dataﬁle > $psﬁle

34.

gedge=0.5 gtri=2.0 n=0 gmax=1

35.

title=“Mi primer Modelo”

36.

labelz=“Depth (km)”labelx=“Distance (km)”

37.

dxnum=1 dznum=1 wbox=10 hbox=7

38. exit 0

En el script se identiﬁcan las siguientes lineas: Sistema: Linea 1, llama al shell y la linea 38 sale del shell. Variables:

• En la linea 5 podemos ver un n´umero por cada vez que corre correctamente el

programa.

• En la linea 9 le asignamos un nombre a la salida del modelo binario y es usado

en la linea 32.

• En la linea 8 a la salida .eps y es usada en la linea 33.

p. 106

Programa trimodel: Lineas 11-32 crean el modelo. Trimodel llena el modelo con tri´angulos de (1 /velocidad2), mientras (1/velocidad) es llamado “slowess”, (1/velocidad2) es llamado “sloth”.

El uso de trimodel puede dividirse en 3 partes:

• La linea 11 deﬁne las dimensiones del modelo.

• Los cinco paquetes de xedge, zadge y sedge son una tripleta que deﬁnen

los limites de la capa y los gradientes de velocidad. Hay un requerimiento importante y es que cada parte de la tripleta tenga el mismo n´umero de valores.Es decir, la tripleta para la capa uno superior del modelo tiene dos valores para cada xedge, zadge y sedge (lineas 12-14), mientras que para la capa dos hay cinco valores para cada tripleta (lineas 15-17).

• La linea de sedge esta siempre de ceros porque todas las capas son

isotropicas y homogeneas.

Programa spsplot: Lineas 33-37 crean el archivo del programa .eps.

p. 107

ANEXO D. C ´ODIGO PARA GENERAR UNA ADQUISICI ´ON EN SU

El c´odigo utilizado para realizar las adquisiciones en SU, es el siguiente: #! /bin/sh # File: acq1.sh # Set messages on # #set -x # Assign values to variables num=1 nangle=201 fangle=-65 langle=65 nt=1370 dt=0.0073 # Name input model ﬁle inmodel=model num.dat # Name output seismic ﬁle outseis=seis num.su #================================================= # Create the seismic traces with “triseis” # i-loop = 40 source positions (disparos) # j-loop = 60 geophone positions (split-spread) # per shot position # k-loop = layers 2 through 10 # (do not shoot layers 1 and 11) # fs= fuentes # sx= posici´on de la fuente

p. 108

# ﬂrd= cabecera n´umero de disparos # fg= posici´on del ge´ofono # gx= posici´on de la fuente # tracl= secuencia num´erica de las trazas # tracf= secuencia num´erica de los disparos echo “ – –Begin looping over triseis.” done i=0 while [ “$i” – ne “1”] do fs=‘bc – l << – END $i * 0.1

END‘

sx=‘bc – l << – END $i * 100

END‘

ﬂdr=‘bc – l << – END

p. 109

$i + 1

END‘

j=0 while [ “$j” – ne “60”] do fg=‘bc – l << – END $i * 0.1 + $j *0.1

END‘

gx=‘bc – l << – END $i * 100 + $j * 100 – 2950

END‘

offset=‘bc – l << –END $j * 100 – 2950

END‘

tracl=‘bc – l << – END $i * 60 + $j + 1

END‘

p. 110

tracf=‘bc – l << – END $j + 1

END‘

echo “ Sx=$sx Gx=$gx ﬂdr=$ﬂdr Offset=$offset tracl=$tracl \ fs=$fs fg=$fg” k=2 while [ “$k” – ne “5”] do triseis < inmodel xs=5,7.9 xg=1.05,10.4 zs=0,0 zg=0,0 \ nangle=$nangle fangle=$fangle langle=$langle \ kreﬂect=$k krecord=1 fpeak=12 lscale=0.5 \ ns=1 fs=$fs ng=1 fg=$fg nt=$nt dt=$dt | suaddhead nt=$nt | sushw key=dt,tracl,tracr,ﬂdr,tracf,trid,offset,sx,gx \ a=4000,$tracl,$tracl,$ﬂdr,$tracf,1,$offset,$sx,$gx >> temp k

p. 111

k=‘expr $k + 1‘ done j=‘expr $j + 1‘ done i=‘expr $i + 1‘ echo “ – –End looping over triseis.” #================================================= # Sum contents of the “temp”ﬁles echo “ – –Sum ﬁles.”

susum temp2 temp3 > tempa susum tempa temp4 > outseis # Remove temp ﬁles echo “ – –Remove temp ﬁles.” rm -f temp* # Exit politely from shell script echo “ – –Finished!” exit

p. 112

Hay algunos par´ametros que se requieren para realizar este proceso de adquisici´on, entre ellos se encuentran:

• Offset m´aximo: El cual hace relaci´on a la distancia m´as lejana en la cual se

encuentra el ´ultimo receptor.

• Intervalo entre puntos de disparo ∆s: Es la separaci´on de las posiciones de

las fuentes (esto aplica cuando se realizan varias perturbaciones, para nuestro seminario s´olo se har´a una).

• Intervalo entre grupos de receptores ∆g: Es la separaci´on entre los ge´ofonos.

• Tiempo de grabaci´on tmax: Es la cantidad de tiempo en que los ge´ofonos se

encuentran grabando o tomando informaci´on. Pero en el proceso de adquisici´on tambi´en son necesarios par´ametros del subsuelo, los cuales son:

• Longitud del modelo (Lx): El cual hace referencia a la extensi´on horizontal del

modelo.

• Profundidad del modelo (Lz).

Ahora bien, se hace necesario determinar los ´ultimos par´ametros de la adquisici´on, los cuales desarrollan un trabajo computacional especial en medio del proceso iterativo, los cuales son:

• Longitud del tendido Lt = 2* offset max.

• N´umero de grupos de receptores ng =Lt/∆g+1.

p. 113

• N´umero de Disparos ns valor ajustado seg´un la longitud del modelo, la longitud

del tendido y el intervalo entre disparos. ns= (L – Lt)/∆g.

• Offset m´ax = (numero de ge´ofonos /2) * ∆g - ∆g/2

• En cada adquisici´on se simul´o el n´umero de muestras por traza fue de 1370

y el intervalo de muestreo en tiempo de 7,3 milisegundos, con el objetivo de que cuando se haga el producto de estos dos valores se tenga un tiempo de grabaci´on de 10 segundos.

• En todos los casos el tipo de levantamiento fue split-spread o sea que se

ubicaron 18 geof´onos (9 a cada lado).

Con esto presente, se realizan unas adquisiciones previas al empalme con los m´etodos de optimizaci´on global, las pruebas realizadas se presentan en la siguiente ﬁgura, all´ı se observa c´omo var´ıa el proceso de adquisici´on dependiendo de la posici´on de la fuente; aclarando que a lado y lado de la fuente se encuentran 9 ge´ofonos. En la Figura D 1. se pueden observar el conjunto de trazas s´ısmicas obtenidas por medio de la adquisici´on realizada en SU.

p. 114

Figura D 1. Adquisici´on con la fuente en 6 Km.

p. 115

ANEXO E. FUNCIONES DE PRUEBA DE SA

Tabla E 1. Funci´on Booth con Probabilidad Inicial de 0,2 y Probabilidad Final de

0,001.

Ensayo Iteraciones M´ınimo Coordenadas Tiempo[s]

(1,1)

17,031446

(1,1)

8,861365

(1,1)

4,347865

(1,1)

5,993299

(1,1)

5,468934

(1,1)

10,868697

(1,1)

19,490887

(1,1)

11,834478

(1,1)

13,53582

(1,1)

28,460658

Promedios

17,5

12,5893449

p. 116

Tabla E 2. Funci´on Booth con Probabilidad Inicial de 0,5 y Probabilidad Final de

0,001.

Ensayo Iteraciones M´ınimo Coordenadas Tiempo[s]

(1;1)

35,538999

(1;1)

27,003056

(1;1)

29,854818

(1;1)

34,428368

(1;1)

18,197358

(1;1)

23,840078

(1;1)

23,91782

(1;1)

26,344367

(1;1)

24,474879

(1;1)

35,233237

Promedios

39,1

20,7959

27,883298

p. 117

Tabla E 3. Funci´on Booth con Probabilidad Inical 0,7 y Probabilidad Final de 0,001. Ensayo Iteraciones M´ınimo Coordenadas Tiempo[s]

(1;1)

36,53786

(1;1)

37,807911

(1;1)

36,40299

(1;1)

37,06048

(1;1)

35,174376

(1;1)

36,918083

(1;1)

34,979867

(1;1)

35,127478

(1;1)

40,890218

(1;1)

36,651197

Promedios

51,7

36,755046

p. 118

Tabla E 4. Funci´on Branin con Probabilidad Inicial de 0,2 y Probabilidad Final de

0,001.

Ensayo Iteraciones M´ınimo Coordenadas Tiempo[s]

27,7029

(1;1)

7,727915

27,7029

(1;1)

10,309123

27,7029

(1;1)

10,398984

27,7029

(1;1)

13,446883

27,7029

(1;1)

13,592214

27,7029

(1;1)

4,572768

27,7029

(1;1)

13,593751

27,7029

(1;1)

5,953403

27,7029

(1;1)

8,135886

27,7029

(1;1)

5,11329

Promedios

27,7029

9,2844217

p. 119

Tabla E 5. Funci´on Branin con Probabilidad Inicial de 0,5 y Probabilidad Final de

0,001.

Ensayo Iteraciones M´ınimo Coordenadas Tiempo[s]

27,7029

(1;1)

34,966416

27,7029

(1;1)

23,02947

27,7029

(1;1)

20,822068

27,7029

(1;1)

25,23541

27,7029

(1;1)

29,565843

27,7029

(1;1)

24,825044

27,7029

(1;1)

32,377747

27,7029

(1;1)

23,61659

27,7029

(1;1)

36,037825

27,7029

(1;1)

31,983286

Promedios

35,2

27,7029

28,2459699

p. 120

Tabla E 6. Funci´on Branin con Probabilidad Inicial de 0,7 y Probabilidad Final de

0,001.

Ensayo Iteraciones M´ınimo Coordenadas Tiempo[s]

27,7029

(1;1)

40,25355

27,7029

(1;1)

36,093901

27,7029

(1;1)

39,30691

27,7029

(1;1)

41,910134

27,7029

(1;1)

41,400628

27,7029

(1;1)

67,181626

27,7029

(1;1)

42,316341

27,7029

(1;1)

41,627265

27,7029

(1;1)

38,727369

27,7029

(1;1)

41,094973

Promedios

51,7

27,7029

42,9912697

p. 121

Tabla E 7. Funci´on Bukin con Probabilidad Inicial de 0,2 y Probabilidad Final de

0,001.

Ensayo Iteraciones M´ınimo Coordenadas Tiempo[s]

5,6231

(-0,34446; 0,0042408)

12,781861

5,6469

(-0,727; 0,0022004)

18,36951

4,4757

(1; 0,011906)

18,826173

1,8116

(-0,30898; 0,00066607)

13,588849

8,9005

(0,57008; 0,010985)

13,518162

10,5613

(1; -0,00092307)

17,969944

4,6036

(-1; 0,012037)

23,224338

8,5317

(-0,82665; 0,013957)

15,97413

5,6685

(0,24785; 0,0037124)

15,984975

6.1872

(0,51001; 0,0063003)

16,567405

Promedios

21,5

6,20101

16,6805347

p. 122

Tabla E 8. Funci´on Bukin con Probabilidad Inicial de 0,5 y Probabilidad Final de

0,001.

Ensayo Iteraciones M´ınimo Coordenadas Tiempo[s]

3,126

(-0,4157; 0,0026462)

29,583638

5,6533

(1 ; 0,069272)

33,585756

2,2892

(-0,68862; 0,0042597)

31,445226

2,2338

(0,18091; -0,00012724)

32,853332

1,7582

(-0,10658; 0,00038889)

31,770657

2,3494

(0,76718; 0,0063881)

33,303359

10,8293

(-0,41356, -0,0098104)

30,179863

3,7697

(-0,64412; 0,0027975)

30,697145

3,2041

(0,96859; 0,0010567)

36,711151

4,6067

(0,012022; 0,96059)

35,871728

Promedios

45,2

3,98197

32,6001855

p. 123

Tabla E 9. Funci´on Bukin con Probabilidad Inicial de 0,7 y Probabilidad Final de

0,001.

Ensayo Iteraciones M´ınimo Coordenadas Tiempo[s]

3,3701

(-0,0174464; -0,0010144)

38,998613

3,9934

(-1; 0,0084764)

38,429666

4,1099

(-0,54346; -0,0015788)

39,667287

5,5084

(-0,40041; -0,0013262)

37,953047

6,6325

(-0,45812; 0,0063721)

36,865758

4,9635

(1; 0,0076443)

38,712096

4,2089

(0,89911; 0,0097649

37,237366

9,3724

(-1; 0,018616)

41,078588

6,0623

(0,29119; -0,0027035)

37,714282

2,6717

(0,25085; 0,0012893)

38,742242

Promedios

3,75305

38,5398945

p. 124

Tabla E 10. Funci´on Rastrigin con Probabilidad Inicial de 0,2 y Probabilidad Final de 0,001.

Ensayo Iteraciones M´ınimo Coordenadas Tiempo[s]

2,1312

(-0,98248; -0,79885)

13,815474

3,3239

(-1; -0,65024)

12,980391

3,643

(0,96166; -0,64731)

17,173177

1,9997

(0,99718; 0,1814)

5,212816

1,1731

(0,021583; -0,69467)

4,166601

2,6714

(1; -0,82165)

5,148462

2,2105

(1; -0,68854)

6,599362

1,4721

(-0,96652; 0,28468)

4,180385

1,7476

(1; -0,71537)

17,318136

3,8004

(-1; -0,63735)

7,138921

Promedios

7,1

2,41729

9,3733725

p. 125

Tabla E 11. Funci´on Rastrigin con Probabilidad Inicial de 0,5 y Probabilidad Final de 0,001.

Ensayo Iteraciones M´ınimo Coordenadas Tiempo[s]

1,5279

(-1; 0,29736)

4,021771

2,8622

(0,075733; -0,82416)

42,708105

1,524

(0,00263345; -0,81634)

48,563944

1,0283

(-0,019248, 0,31629)

49,79108

1,9937

(-1;, 0,1764)

29,465846

4,8499

(-1; 0,39146)

45,824459

3,6888

(-0,8918; -0,72224)

33,525393

1,5896

(-1; 0,3005)

37,112577

1,6931

(0,9396; 0,23515)

48,772617

4,4741

(1; 0,38372)

34,711

Promedios

35,4

2,52316

37,4496792

p. 126

Tabla E 12. Funci´on Rastrigin con Probabilidad Inicial de 0,7 y Probabilidad Final de 0,001.

Ensayo Iteraciones M´ınimo Coordenadas Tiempo[s]

1,3024

(0,012979; 0,16989)

46,332541

2,1926

(0;057069; 0,16087)

52,456074

1,5239

(-1; 0,29715)

49,984358

3,857

(0,95595; 0,1339)

52,282793

1,1846

(1; 0,2736)

47,901148

5,2316

(-0,90328; -0,64206)

55,875822

1,0627

(-1; 0,25037)

51,627547

0,58696

(-0,0051408; 0,30011)

49,924146

1,2056

(1; 0,22182)

57,238527

0,61318

(0,014077; -0,75469)

53,76721

Promedios

51,6

1,876054

51,7390166

p. 127

Tabla E 13. Funci´on Rosenbrock con Probabilidad Inicial de 0,2 y Probabilidad Final de 0,001.

Ensayo Iteraciones M´ınimo Coordenadas Tiempo[s]

2,1261

(-0,66937; 0,59)

8,192266

0,23849

(0,56256; 0,33819)

12,630277

(1; 1)

7,479052

6,9289

(-0,0055381; -0,24323)

5,662937

1,9449

(-0,22385; 0,11697)

6,135395

1,7764

(0,63764; 0,53485)

10,144735

(1;1)

2,997208

0,82797

(0,69101; 0,39191)

8,51312

(1;1)

5,449873

3,9164

(-0,60615; 0,25181)

12,00638

Promedios

8,7

1,775916

7,9211243

p. 128

Tabla E 14. Funci´on Rosenbrock con Probabilidad inicial de 0,5 y Probabilidad Final de 0,001.

Ensayo Iteraciones M´ınimo Coordenadas Tiempo[s]

1,5279

(-1; 0,29736)

4,021771

2,8622

(0,075733; -0,82416)

42,708105

1,524

(0,00263345; -0,81634)

48,563944

1,0283

(-0,019248; 0,31629)

49,79108

1,9937

(-1; 0,1764)

29,465846

4,8499

(-1; 0,39146)

45,824459

3,6888

(-0,89188; -0,72224)

33,525393

1,5896

(-1; 0,3005)

37,112577

1,6931

(0,9396; 0,23515 )

48,772617

4,4741

(1; 0,38372)

34,711

Promedios

35,4

2,52316

37,4496792

p. 129

Tabla E 15. Funci´on Rosenbrock con Probabilidad inicial de 0,7 y Probabilidad Final de 0,001.

Ensayo Iteraciones M´ınimo Coordenadas Tiempo[s]

1,3024

(0,01297; 0,16989)

46,332541

2,1926

(0,057069; 0,16087)

52,456074

1,5239

(-1; 0,29715)

49,984358

3,857

(0,95595; 0,1339)

52,282793

1,1846

(1; 0,2736)

47,901148

5,2316

(-0,90328; -0,64206)

55,875822

1,0627

(-1; 0,25037)

51,627547

0,58696

(-0,0051408; 0,30011)

49,924146

1,2056

(1; 0,22182)

57,238527

0,61318

(0,014077; -0,75469)

53,76721

Promedios

51,6

1,876054

51,7390166

p. 130

ANEXO F. FUNCIONES DE PRUEBA DE ABC

Tabla F 1. Funci´on Booth con un Tama˜no de poblaci´on de 10. Ensayo Iteraciones Tiempo[s] M´ınimo promedio Desviaci´on est´andar

0,627928

0,207973

0,393534

0,599649

0,562272

0,966198

0,576885

0,0848997

0,12107

0,670745

1,34586

2,05576

0,625248

0,293219

0,402777

0,627221

0,480007

0,540871

0,656699

1,40349

1,12775

0,674627

1,30273

2,30914

0,642411

1,64886

2,96756

0,776886

3,96259

6,00768

Promedios

23,2

0,6478299

p. 131

Tabla F 2. Funci´on Booth con un Tama˜no de poblaci´on de 20. Ensayo Iteraciones Tiempo[s] M´ınimo promedio Desviaci´on est´andar

0,808465

0,00220355

0,000314121

0,858929

0,00755903

0,0085583

0,965008

0,000995091

0,001693

0,87541

0,162044

0,120316

1,026993

0,000110482

0,000226

0,902882

0,051898

0,0688626

0,862608

0,00142141

0,0019895

0,953311

0,00342254

0,0054137

0,989224

0,000363

0,000811695

0,82366

0,000150267

0,00033

Promedios

22,1

0,906649

p. 132

Tabla F 3. Funci´on Booth con un Tama˜no de poblaci´on de 30. Ensayo Iteraciones Tiempo[s] M´ınimo promedio Desviaci´on est´andar

1,398023

0,00683423

0,006556

1,017035

0,004006267

0,0069635

1,189788

0,0441296

0,04027

0,943701

0,001355

0,00302988

1,03381

0,0272755

0,597006

1,167543

0,0154997

0,019057

1,201479

0,00822195

0,0065976

1,21231

0,0131186

0,0261644

1,20449

0,00105739

0,0017456

1,009865

0,00595088

0,0106864

Promedios

19,8

1,1378044

p. 133

Tabla F 4. Funci´on Branin con un Tama˜no de poblaci´on de 10. Ensayo Iteraciones Tiempo[s] M´ınimo promedio Desviaci´on est´andar

0,642256

0,993819

1,35432

0,639292

0,2116

0,296968

0,963216

0,0803623

0,179696

0,633387

0,651061

0,418038

0,674732

0,907386

0,737599

0,676721

0,201872

0,28258

0,693487

0,0797342

0,178291

0,720508

0,162313

0,222322

0,704135

0,253239

0,231616

0,621747

0,0802534

0,179452

Promedios

25,9

0,6969481

p. 134

Tabla F 5. Funci´on Branin con un Tama˜no de poblaci´on de 20. Ensayo Iteraciones Tiempo[s] M´ınimo promedio Desviaci´on est´andar

0,894075

0,169776

0,233075

0,964737

0,0795777

0,177941

0,781096

0,161871

0,2216784

0,843716

0,28357

0,268687

0,863056

0,161203

0,220741

1,066318

0,615778

0,211027

0,838531

0,258179

0,237825

0,881517

0,0799318

0,178733

0,80001

0,160073

0,219194

0,854382

0,163017

0,23223

Promedios

22,3

0,8787438

p. 135

Tabla F 6. Funci´on Branin con un Tama˜no de poblaci´on de 30. Ensayo Iteraciones Tiempo[s] M´ınimo promedio Desviaci´on est´andar

1,005547

0,40538

0,00813793

1,485296

0,0797275

0,178276

1,049545

0,324866

0,181702

1,365495

0,399116

0,00211089

1,250349

0,079647

0,178096

1,123669

0,0797559

0,17834

1,098408

0,22054

0,230146

1,128752

0,336165

0,190369

0,922161

0,159295

0,218123

1,226955

0,16018

0,21934

Promedios

19,5

1,1656168

p. 136

Tabla F 7. Funci´on Bukin con un Tama˜no de poblaci´on de 10. Ensayo Iteraciones Tiempo[s] M´ınimo promedio Desviaci´on est´andar

0,665159

20,2703

18,1745

0,575441

24,1911

2,58944

0,601846

2,87323

6,42474

0,602976

3,542374

5,08724

0,635585

6,14425

9,84178

0,637515

10,2145

9,54007

0,643211

1,46505

1,78629

0,670389

2,29925

3,53864

0,638906

6,48498

10,007

0,71502

9,29924

11,0564

Promedios

18,5

0,6386048

p. 137

Tabla F 8. Funci´on Bukin con un Tama˜no de poblaci´on de 20. Ensayo Iteraciones Tiempo[s] M´ınimo promedio Desviaci´on est´andar

0,756931

5,90494

4,47574

0,717696

9,17411

7,46124

0,73103

12,6827

10,3719

0,758678

0,289246

0,646774

0,821151

1,99366

2,91593

0,87101

2,6722

1,56575

0,845718

14,5974

8,72356

0,840134

21,8381

15,7532

0,81758

2,88613

2,18956

0,798047

3,25096

4,04522

Promedios

16,3

0,7957975

p. 138

Tabla F 9. Funci´on Bukin con un Tama˜no de poblaci´on de 30. Ensayo Iteraciones Tiempo[s] M´ınimo promedio Desviaci´on est´andar

1,008591

7,29274

4,60071

1,013786

7,92102

3,585996

0,869272

3,02457

4,69554

1,069895

2,64524

3,99678

0,933024

0,858328

1,91928

0,879774

3,66688

3,6021

1,101706

0,533579

1,19312

1,050762

1,43853

2,02503

0,981611

2,32733

2,49255

1,069933

6,20128

4,95076

Promedios

16,4

0,9978354

p. 139

Tabla F 10. Funci´on Rastrigin con un Tama˜no de poblaci´on de 10. Ensayo Iteraciones Tiempo[s] M´ınimo promedio Desviaci´on est´andar

0,684543

4,33021

3,05743

0,605824

1,81963

2,6258423

0,670606

0,498027

0,681531

0,633755

0,209434

0,468309

0,643395

0,0150346

0,0336183

0,651303

2,32723

2,08013

0,626677

1,17222

1,11418

0,614568

0,605687

0,899586

0,645789

3,47266

2,33551

0,660467

0,405526

0,555289

Promedios

21,9

0,6436927

p. 140

Tabla F 11. Funci´on Rastrigin con un Tama˜no de poblaci´on de 20. Ensayo Iteraciones Tiempo[s] M´ınimo promedio Desviaci´on est´andar

1,043412

0,0396587

4,47574

0,971173

0,643202

0,0819043

0,840069

0,466004

0,770038

0,8743

1,25742

0,700969

0,84731

0,410546

0,575045

0,913259

0,0219638

0,043007

0,902152

6,03388*e-7 1,34992*e-6

0,913628

0,000198715

0,000444341

1,029203

0,0156295

0,0344245

0,870698

3,1333

2,31464

Promedios

25,3

0,9205204

p. 141

Tabla F 12. Funci´on Rastrigin con un Tama˜no de poblaci´on de 30. Ensayo Iteraciones Tiempo[s] M´ınimo promedio Desviaci´on est´andar

1,196334

0,000212319

0,000474459

1,316803

0,00496159

0,0103685

1,034993

0,513015

0,7300748

1,243897

0,384496

0,483485

1,096165

0,35538

0,427005

1,275771

0,0327445

0,0,49975

1,030727

0,485644

0,998241

1,41999

0,00694052

0,0154

1,072445

0,0080354

0,0179677

1,07576

0,257464

0,438388

Promedios

1,1762885

p. 142

Tabla F 13. Funci´on Rosenbrock con un Tama˜no de poblaci´on de 10. Ensayo Iteraciones Tiempo[s] M´ınimo promedio Desviaci´on est´andar

0,636362

0,0850987

0,190287

0,592263

0,7009136

0,859342

0,562416

3,26463

4,38528

0,60259

1,10827

2,47817

0,593234

2,90247

4,83953

0,577023

2,76233

3,23293

0,628695

0,06566

0,135343

0,655905

9,40164

11,4295

0,559082

3,22297

5,02619

0,571911

0,689083

1,44936

Promedios

15,6

0,5979481

p. 143

Tabla F 14. Funci´on Rosenbrock con un Tama˜no de poblaci´on de 20. Ensayo Iteraciones Tiempo[s] M´ınimo promedio Desviaci´on est´andar

0,697612

0,978903

2,18889

0,798947

0,0569814

0,118452

0,671397

2,91043

3,30841

0,751023

0,054773

0,122476

0,763198

0,0270995

0,0605963

0,793177

4,54917

4,03817

0,773375

3,72498

3,47252

0,768195

4,58192

4,09884

0,793086

0,144857

0,211715

0,672811

0,178097

0,251738

Promedios

14,7

0,7482821

p. 144

Tabla F 15. Funci´on Rosenbrock con un Tama˜no de poblaci´on de 30. Ensayo Iteraciones Tiempo[s] M´ınimo promedio Desviaci´on est´andar

1,050322

0,0860868

0,192496

0,988424

0,0975121

0,161486

0,982141

0,251952

0,373073

1,037878

0,0303004

0,0677538

0,919991

0,00228247

0,00510376

0,935498

0,071211

0,143139

0,79611

0,151366

0,262939

1,025568

0,0630777

0,141046

1,073455

0,0943065

0,0881474

0,977257

2,48214

3,68793

Promedios

15,3

0,9786644

p. 145

ANEXO G. FUNCI ´ON AUTOM ´ATICA

function [ z ]= inversion automatica ( x )

%x =[0.77

0.5 0.26

0 . 1 2 ] ; % Valor sloths modelo de referencia

% 1. Creo el archivo que modificara al modelo

str6 = s t r c a t ( ’ /home/ User / inversion matlab / ’ , num2str ( i ) , ’ / capas . sh ’ ) ;

f i l e I D =fopen ( str6 , ’w ’ ) ;

f p r i n t f ( f i l e I D , ’ v1 s f i l l =0.01 ,1.01 ,0.0 ,0 , % f ,0.0 ,0.0 %s\n ’ , x ( 1 ) , ’ \ ’ ) ;

f p r i n t f ( f i l e I D , ’ s f i l l =0.01 ,4.00 ,0.0 ,0 , % f ,0.0 ,0.0 %s\n ’ , x ( 2 ) , ’ \ ’ ) ;

f p r i n t f ( f i l e I D , ’ s f i l l =5.00 ,5.50 ,0.0 ,0 , % f ,0.0 ,0.0 %s\n ’ , x ( 3 ) , ’ \ ’ ) ;

f p r i n t f ( f i l e I D , ’ s f i l l =0.10 ,6.00 ,0.0 ,0 , % f ,0.0 ,0.0 %s\n ’ , x ( 4 ) , ’ \ ’ ) ;

fclose ( f i l e I D ) ;

% 2.

Borro el modelo actualizado . sh para luego modificarlo

str7 = s t r c a t ( ’ export LD\ LIBRARY\ PATH=/ usr / lib64 ;

cd /home/ User / inversion matlab ; . / Clean . sh ’ ) ;

d= unix ( str7 ) ;

% 3.

Escribo la primera parte de el archivo de modelo actualizado

str8 = s t r c a t ( ’ export LD\ LIBRARY\ PATH=/ usr / lib64 ;

cd /home/ User / inversion matlab / ; cat modelado1 . sh |

%sed −n ”1 ,27p ” >> modelo\ actualizado . sh ’ ) ;

d= unix ( str8 ) ;

% 4.

Escribo las l i n e s nuevas con los nuevos valores de slowlness

str9 = s t r c a t ( ’ export LD\ LIBRARY\ PATH=/ usr / lib64 ;

cd /home/ User / inversion matlab / ; cat capas . sh |

sed −n ”1 ,4p ” >> modelo\ actualizado . sh ’ ) ;

d= unix ( str9 ) ;

% 5.

Escribo la tercera parte del archivo del modelo actualizado

str10= s t r c a t ( ’ export LD\ LIBRARY\ PATH=/ usr / lib64 ;

cd /home/ User / inversion matlab / ; cat modelado1 . sh |

sed −n ”32 ,38p ” >> modelo\ actualizado . sh ’ ) ;

d= unix ( str10 ) ;

% 6. Genero el modelo

str1 = s t r c a t ( ’ export LD\ LIBRARY\ PATH=/ usr / lib64 ;

p. 146

cd /home/ User / inversion matlab / ; sh modelo\ actualizado . sh ’ ) ;

d= unix ( str1 ) ;

% 7.

Simulo la adquisicion y obtengo las trazas

str2 = s t r c a t ( ’ export LD\ LIBRARY\ PATH=/ usr / lib64 ;

cd /home/ User / inversion matlab / ; sh acq1 . sh ’ ) ;

d= unix ( str2 ) ;

% 8. Para ver las trazas sismicas

str3 = s t r c a t ( ’ export LD\ LIBRARY\ PATH=/ usr / lib64 ;

cd /home/ User / inversion matlab / ; suxwigb < seis1 . su ’ ) ;

d= unix ( str3 ) ;

% 9. Para ver el sismograma

str4 = s t r c a t ( ’ export LD\ LIBRARY\ PATH=/ usr / lib64 ;

cd /home/ User / inversion matlab / ; suximage < seis1 . su ’ ) ;

d= unix ( str4 ) ;

% 10. Cambio de formato SU a SEGY para la traza actualizada

str5 = s t r c a t ( ’ export LD LIBRARY PATH=/ usr / lib64 ; , . . .

cd /home/ User / inversion matlab / ; sh sutosegy ’ ) ;

d= unix ( str5 ) ;

% 11. Copiar la traza actualizada . sgy en la carpeta de la p a r t i c u l a

str6 = s t r c a t ( ’ export LD LIBRARY PATH=/ usr / lib64 ; , . . .

cd /home/ User / inversion matlab / ’ , num2str ( i ) , ’ ; , . . .

cp Trazas actualizadas . sgy /home/ User / Libardo inversion matlab / ’ , num2str ( i ) , ’ / ’ ) ;

d= unix ( str6 ) ;

% 12. Importo las trazas a matlab en forma de estructuras

str7 = s t r c a t ( ’ /home/ User / inversion matlab / ’ , num2str ( i ) , ’ / Trazas actualizadas . sgy ’ ) ;

estructura1= r e a d s e g y f i l e ( str7 ) ;

str8 = s t r c a t ( ’ /home/ User / inversion matlab / ’ , num2str ( i ) , ’ / Trazas reales . sgy ’ ) ;

estructura2= r e a d s e g y f i l e ( str8 ) ;

% 13. Grafica de las trazas real y actual s wplot ( estructura1 . traces ) ;

s wplot ( estructura2 . traces ) ;

% 14.

Calculo los coeficientes de correlacion entre las trazas s i n t e t i c a s y reales .

f o r j =1:17

w=estructura1 . traces ( : , j ) ;

y=estructura2 . traces ( : , j ) ;

%calcula el c o e f i c i e n t e de correlacion cruzada

p. 147

coeficientes ( : , j )= xcorr (w, y , ’ coeff ’ ) ;

end

f o r k=1:17

max coef ( k ,1)=max( coeficientes ( : , k ) ) ;

end

z=max coef ; Listing 1. Inversion Automatica.m

p. 148

ANEXO H. ALGORITMO PARA SIMULATED ANNEALING

clc ; clear a l l ; close a l l ;

%addpath ( ’ / home/ User cps / SeisLab 10 .0301/S4M/ Geophysics 3 . 0 ’ ) ;

Algoritmo Simulated Annealing

Grupo 2.

Seminario de investigacion : Implementacion de tecnicas de

% optimizacion global en problemas de inversion sismica .

%%Fecha y hora

S = date ;

%Se i n i c i a l i z a el tiempo de computo que se va a u t i l i z a r

elapsedTime = 0;

t i c

%% Superficie

% Designando variables para los puntos de la g r a f i ca .

Xmin=0; % Limite i n f e r i o r en la capa 1

Xmax=1; % Limite superior en la capa 1

Ymin=0; % Limite i n f e r i o r en la capa 2

Ymax=1; % Limite superior en la capa 2

Zmin=0; % Limite i n f e r i o r en la capa 3

Zmax=1; % Limite superior en la capa 3

Hmin=0; % Limite i n f e r i o r en la capa 4

Hmax=1; % Limite superior en la capa 4

% Definicion fronteras

bnd=[Xmin Xmax

Ymin Ymax

Zmin Zmax

Hmin Hmax ] ;

% Posicion I n i c i a l

X i n i ( 1 , : ) = bnd ( : , 1 ) + rand ( size ( bnd , 1 ) , 1 ) . ∗( bnd ( : , 2 ) −bnd ( : , 1 ) ) ; %Modo a l e a t o r i o

% de asignar las capas

% Seccion para f i j a r capas manualmente

% X i n i (1 ,2)=0.5;

X i n i (1 ,3)=0.26;

X i n i (1 ,4)=0.12;

% Valores de las velocidades de referencia

%X i n i =[0.77

0.5 0.26

0 . 1 2 ] ;

p. 149

%% Simulated Annealing

%Numero de iteraciones

n = 5;

%Numero de ensayos por c i c l o

m = 5; %m = 50;

%Numero de soluciones aceptadas

na = 1.0;

%Normalizacion

dt =0.1;

% Probabilidad de aceptar la peor solucion al i n i c i o .

p1 = 0.7;

% Probabilidad de aceptar la peor solucion al f i n a l .

p50 = 0.001; %p50 = 0.001;

% Temperatura I n i c i a l

t1 = −1.0/ log ( p1 ) ;

% temperatura f i n a l

t50 = −1.0/ log ( p50 ) ;

%Contador para saber cuantas veces acepta valores .

contador =0;

% fraccion de reducion por cada c i c l o

frac = ( t50 / t1 ) ˆ ( 1 . 0 / ( n −1.0));

% I n i c i a l i z a c i o n de x

x ( 1 , : ) = X i n i ;

% Mejores resultados actuales por el momento .

xc = X i n i ;

fc = inversion automatica3 ( xc ) ; %Evaluacion en la inversion .

fs = zeros ( n+1 ,17); %Almacena cada uno de los valores fc aprobados

por ser la primera evaluacion , es guardado directamente como punto de partida ) .

fs ( 1 , : ) = fc ;

% Temperatura actual .

t = t1 ;

%I n i c i a l i z a c i o n de DeltaE para el avance

DeltaE = 0.1;

% Promedio de DeltaE ( Para normalizar el exponente de la exponencial

% decreciente de temperatura .

DeltaE avg = DeltaE ;

p. 150

%Matriz de unos para encontrar el error medio cuadratico .

unos=ones (17 ,1);

% Seguimiento , para enterarse de cuantas veces ingresa a las iteraciones cuando

% decrece la exponencial .

contador aceptadas igualesque cero =0;

contador aceptadas menoresque cero =0;

f o r i =1:n

i f ( t > t50 )

f o r j =1:m

% Generacion de nuevos puntos de exploracion para dos capas f i j a s y

% dos a l e a t o r i a s

x i (1) = xc (1) + ( rand ( ) −0.5)∗dt ;

x i (2) = xc (2) + ( rand ( ) −0.5)∗dt ;

x i (3) = 0.26;

x i (4) = 0.12;

100

x i (3) = xc (3) + ( rand ( ) −0.5)∗dt ;

101

x i (4) = xc (4) + ( rand ( ) −0.5)∗dt ;

102

103

% Acotamiento de fronteras para los valores maximos y minimos .

104

%Limites superiores e i n f e r i o r e s

105

x i (1) = max( min ( x i ( 1 ) , 1 . 0 ) , 0 ) ;

106

x i (2) = max( min ( x i ( 2 ) , 1 . 0 ) , 0 ) ;

107

x i (3) = max( min ( x i ( 3 ) , 1 . 0 ) , 0 ) ;

108

x i (4) = max( min ( x i ( 4 ) , 1 . 0 ) , 0 ) ;

109

110

% Interaccion con el proceso de inversion .

” inversion automatica3 ”

111

% es el archivo modificado suministrado por el profesor Abreo .

112

% Recordando que este devuelve un valor de correlacion cruzada .

113

momentaneo=inversion automatica3 ( x i ) ;

114

115

restas2=unos−abs (momentaneo ) ;

116

error medio actual= sqrt (sum( restas2 . ˆ 2 ) / 1 7 ) ;

117

118

% Evaluacion de los puntos previos para comparar .

Sabiendo que fc

119

% es el punto que c i c l o a c i c l o va quedando como el mas apropiado .

120

restas=unos−abs ( fc ) ;

121

error medio anterior = sqrt (sum( restas . ˆ 2 ) / 1 7 ) ;

122

123

% Matriz de avance para determinar la posicion del gradiente

124

%con las evaluaciones en la inversion y los nuevos puntos .

125

Matriz de avance ( i , j )=( error medio actual−error medio anterior ) ;

p. 151

% U t i l i z a d o para la normalizacion .

128

DeltaE = error medio actual−error medio anterior ;

129

130

%% Para el caso donde el error medio actual es menor que el anterior ,

131

%%quiere decir que se tiene un valor mejor , o sea que es tomado y evaluado .

132

i f ( error medio actual −error medio anterior <= 0) ,

133

134

xc2 (1) = x i ( 1 ) ;

135

xc2 (2) = x i ( 2 ) ;

136

xc2 (3) = x i ( 3 ) ;

137

xc2 (4) = x i ( 4 ) ;

138

139

% Bandera para contar los puntos aceptados .

140

contador aceptadas menoresque cero=contador aceptadas menoresque cero +1;

141

142

fc = inversion automatica3 ( xc2 ) ;

143

144

%% Para el caso donde se tiene un peor valor de error de correlacion ,

145

%% no descarta el peor punto , sino que genera una probabilidad de aceptacion

146

%%de ese punto basado en el valor a l e a t o r i o rand generado .

147

148

else

149

150

%GENERA LA PROBABILIDAD DE ACEPTACION.

151

p1 = exp(−(DeltaE / ( DeltaE avg ∗t ) ) ) ;

152

153

% Determinacion de cual es el punto que se acepta .

154

i f ( rand ( ) < p1 )

155

156

% Aceptacion de la peor solucion

157

accept = true ;

158

159

else

160

161

%No acepta la peor solucion .

162

accept = false ;

163

164

end

165

166

i f ( accept==true )

167

168

% Actualiza cada punto en la solucion aceptada .

p. 152

xc2 (1) = x i ( 1 ) ;

170

xc2 (2) = x i ( 2 ) ;

171

xc2 (3) = x i ( 3 ) ;

172

xc2 (4) = x i ( 4 ) ;

173

174

contador=contador +1; % Contador de aceptaciones

175

% Incremento del numero de soluciones aceptadas

176

na = na + 1.0;

177

% Actualizacion DeltaE avg

178

DeltaE avg = ( DeltaE avg ∗( na−1.0) + DeltaE ) / na ;

179

180

e l s e i f ( accept == false )

181

xc2 (1) = xc ( 1 ) ;

182

xc2 (2) = xc ( 2 ) ;

183

xc2 (3) = xc ( 3 ) ;

184

xc2 (4) = xc ( 4 ) ;

185

end

186

187

end %( error medio anterior >= error medio actual )

188

189

fc = inversion automatica3 ( xc2 ) ;

190

191

end %f o r j =1:m

192

193

convergencias ( i )= contador ;

194

contador de aceptaciones=convergencias ’ ;

195

contador =0;

196

xc=xc2 ; % Reemplazo del nuevo punto .

197

198

% Graba el mejor valor de X evaluado al f i n a l de cada i n te n t o .

199

x ( i +1 ,1) = xc2 ( 1 ) ;

200

x ( i +1 ,2) = xc2 ( 2 ) ;

201

x ( i +1 ,3) = xc2 ( 3 ) ;

202

x ( i +1 ,4) = xc2 ( 4 ) ;

203

204

% Almacenamiento de los puntos aceptados

205

fs ( i +1 ,:) = fc ;

206

207

% Reduccion de temperatura para el siguiente c i c l o .

208

t = frac ∗t1 ;

209

210

time . horas=toc /3600;

211

time . minutos=toc /60;

p. 153

time . segundos=toc ;

213

%−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−

214

215

SA( k ) = s t r u c t ( ’ Fecha ’ ,S , . . .

216

’ Algoritmo ’ , ’SA ’ , . . .

217

’ Arquitectura ’ , ’Osaka ’ , . . .

218

’ Valores Parametros ’ , [ k ; n ;m; p1 ; p50 ; t1 ; t50 ; dt ] , . . .

219

’ Tiempo de Computo s ’ , time , . . .

220

’ Tipo de Metrica ’ , ’ Correlacion cruzada ’ , . . .

221

’ Iteraciones realizados ’ , i , . . .

222

’ Correlacion ’ , fc , . . .

223

’ Matriz de Correlcion ’ , fs , . . .

224

’ Matriz de avance ’ , Matriz de avance , . . .

225

’ Contador de aceptaciones de iteraciones ’ , contador de aceptaciones , . . .

226

’ S l o t h s i n i c i a l e s ’ , X ini , . . .

227

’ S l o t h s f i n a l e s ’ , x , . . .

228

’ Observaciones ’ , ’ Pruebas de SA realizadas en Lenovo .

229

C r i t e r i o de parada modificado , se ha u t i l z a d o el proceso de

230

screen en la terminal de Ubuntu .

El c r i t e r i o de parada

231

no esta funcionando , se detiene solo por el numero

232

de iteraciones . ’ ) ;

233

234

235

end %f o r i =1:n

236

237

save SA

238

239

%% Resultados

240

% Imprime la solucion

241

242

disp ( [ ’ Valores de las velocidades i n i c i a l e s :

’ , num2str ( X i n i ) ] )

243

disp ( [ ’ Valores de las velocidades encontradas :

’ , num2str ( x ( end , : ) ) ] )

244

245

%% Imagen del modelo obtenido por medio de SA

246

str1xx= s t r c a t ( ’ export LD LIBRARY PATH=/ usr / lib64 ;

247

cd /home/ User / inversion matlab / ; evince model1 . eps ’ ) ;

248

d= unix ( str1xx ) ;

249

end Listing 2. SA.m

p. 154

ANEXO I. ALGORITMO PARA ARTIFICIAL BEE COLONY

clear a l l ; close a l l ; clc

Algoritmo A r t i f i c i a l Bee Colony

Seminario de investigacion : Implementacion de tecnicas de

% optimizacion global en problemas de inversion sismica .

% Fuente :

A r t i f i c i a l Bee Colony Algorithm Homepage

%http : / / mf . erciyes . edu . t r / abc /

%% Control de los parametros globales de ABC

NP=10; %Tamano de la colonia

FoodNumber=NP/ 2 ; % Fuentes de alimento

% ( Mitad de la colonia )

l i m i t =100; % Limite de intentos para que la

% obrera abandone la fuente de comida

maxCycle=10; %Numero de iteraciones

NumGeof=17; %Numero de geofonos .

%/∗ Variables ∗/

objfun= ’ inversion \ automatica3 ’ ; % Funcion o b j e t i v o

D=4; %Numero de parametros a optimizar ,

% dimensiones

ub=ones (1 ,D)∗1; % Valores maximos de f r o n t e r a

lb=ones (1 ,D) ∗( 0 ) ; % Valores minimos de f r o n t e r a

runtime =1; % Cantidad de veces que se ejecuta el

% algoritmo

contador = 0;

S=date ;

GlobalMins=zeros (1 , runtime ) ;

f o r r =1: runtime

t i c

i t e r =1;

% I n i c i a l i z a c i o n de las fuentes de comida

% Variables i n i c i a l i z a d a s en el rango [ lb , ub ] .

p. 155

Range = repmat ( ( ub−lb ) , [ FoodNumber 1 ] ) ;

Lower = repmat ( lb , [ FoodNumber 1 ] ) ;

Foods = ( rand (FoodNumber ,D) .∗Range + Lower ) ’ ; % Generacion

%de puntos a l e a t o r i o s para la busqueda de soluciones ,

%” Matriz Alimento ”

f o r s=1:FoodNumber , % Para pruebas estaticas y guiadas

% se reemplazan los valores que desee

% Foods (1 , s )=0.77;

% Foods (2 , s )=0.5;

Foods (3 , s )=0.26;

Foods (4 , s )=0.12;

end

f o r w=1:FoodNumber , % Evaluacion de las abejas en la funcion

% o bj et ivo

ObjVal ( : ,w)= feval ( objfun , Foods ( : ,w) ,w) ;

end

f o r q=1:FoodNumber , % Calculo del Fitness para cada

% componente de una abeja

Fitness vector ( : , q)= calculateFitness ( ObjVal ( : , q ) ) ;

end

f o r z=1:FoodNumber , % Organizacion de cada Fitness

% por abeja evaluada

Fitness ( : , z)= sqrt (sum ( ( Fitness vector ( : , z ) ) . ˆ 2 ) / NumGeof ) ;

end

% Reinicio del contador de intentos

t r i a l =zeros (1 ,FoodNumber ) ;

% Mejor fuente de alimento memorizado hasta el momento

[ GlobalMin , BestInd ]= min ( Fitness ) ;

p. 156

GlobalParams=Foods ( : , BestInd ) ; %Toma como parametros

%globales el minimo encontrado por ahora junto

%con los valores correspondiente de la matriz alimento

while ( ( i t e r <= maxCycle ) ) ,

%%FASE DE ABEJAS OBRERAS % %

f o r i =1:(FoodNumber)

% El parametro a cambiar es determinado

% aleatoriamente

Param2Change= f i x ( rand∗D)+1;

100

% Seleccion de una solucion a l e a t o r i a vecina para

101

% producir una posible solucion

102

neighbour= f i x ( rand ∗(FoodNumber ) ) + 1 ;

103

104

% Posible caso donde aleatoriamente se toma

105

%la solucion ” i ” que se esta evaluando , buscando

106

%cambiarlo para que otra sea la que se evalue

107

108

while ( neighbour== i )

109

neighbour= f i x ( rand ∗(FoodNumber ) ) + 1 ;

110

end ;

111

112

sol ( : , i )=Foods ( : , i ) ;

113

114

Aplicacion de la ecuacion v\ { i j }=x\ { i j }+

115

%\ phi \ { i j }∗( x\ { k j}−x\ { i j })

116

sol (Param2Change , i )=Foods (Param2Change , i )

117

+(Foods (Param2Change , i )−Foods (Param2Change , neighbour ) )

118

∗( rand −0.5)∗2;

119

120

121

Acotamiento de las fronteras

122

sol= max( min ( sol ’ , ub ( 1 ) ) ’ , lb ( 1 ) ) ;

123

124

% Evaluacion de la nueva solucion ” mutante ”

125

f o r e=1:FoodNumber ,

p. 157

ObjValSol ( : , e)= feval ( objfun , Foods ( : , e ) , e ) ;

127

end

128

129

f o r v=1:FoodNumber ,

130

Fitness \ vectorSol ( : , v)= calculateFitness

131

( ObjVal ( : , v ) ) ;

132

end

133

134

f o r h=1:FoodNumber ,

135

FitnessSol ( : , h)= sqrt (sum ( ( Fitness \ vectorSol

136

( : , h ) ) . ˆ 2 ) / NumGeof ) ;

137

end

138

139

[ GlobalMinSol , BestIndSol ]= min ( FitnessSol ) ;

140

% minimo y su posicion del vector f i t n e s s

141

GlobalParamsSol=Foods ( : , BestIndSol ) ;

142

% Seleccion de aquel vector que tiene el mejor f i t n e s s

143

144

%Una seleccion prematura es aplicada entre la solucion

145

% actual y una mutante

146

147

i f ( GlobalMinSol < GlobalMin ) % Si la solucion del

148

% mutante es mejor que la actual , reemplaza la solucion

149

% mutante y r e i n i c i a el contador de intentos

150

Foods ( : , i )= sol ( : , BestIndSol ) ;

151

Foods test ( : , BestIndSol )= sol ( : , BestIndSol ) ;

152

Fitness ( : , BestIndSol )= FitnessSol ( : , BestIndSol ) ;

153

ObjVal ( : , BestIndSol )= ObjValSol ( : , BestIndSol ) ;

154

t r i a l ( i )=0;

155

else

156

t r i a l ( i )= t r i a l ( i )+1; % Si la solucion no mejora ,

157

%incrementa en contador de intentos ” t r i a l ”

158

end ;

159

end ;

160

161

%% Calculo de probabilidades

% %

162

163

prob =(0.9.∗Fitness . / max( Fitness ) ) + 0 . 1 ;

164

165

%% % % %

FASE DE ABEJAS OBSERVADORAS % % % %

166

i =1;

167

t =0;

168

while ( t<FoodNumber)

p. 158

i f ( rand<prob )

170

t = t +1;

171

172

% El parametro a cambiar es determinado aleatoriamente

173

Param2Change= f i x ( rand∗D)+1;

174

175

% Seleccion de una solucion a l e a t o r i a vecina

176

% para producir una posible solucion mutante

177

neighbour= f i x ( rand ∗(FoodNumber ) ) + 1 ;

178

179

% Posible caso donde aleatoriamente se tome la

180

% solucion i que se esta evaluando , busca cambiarlo

181

% para que otra sea la que se evalue

182

while ( neighbour== i )

183

neighbour= f i x ( rand ∗(FoodNumber ) ) + 1 ;

184

end ;

185

186

% Los puntos de Foods son u t i l i z a d o s para generar

187

% puntos a l e a t o r i o s

188

sol ( : , i )=Foods ( : , i ) ;

189

190

% Aplicacion de la ecuacion v\ { i j }=x\ { i j }+\ phi \ { i j }

191

%∗( x\ { k j}−x\ { i j })

192

sol (Param2Change)=Foods (Param2Change , i )+

193

( Foods (Param2Change , i )−Foods (Param2Change , neighbour ) )

194

∗( rand −0.5)∗2;

195

196

% Acotamiento de las fronteras de los nuevos

197

puntos generados

198

sol= max( min ( sol ’ , ub ( 1 ) ) ’ , lb ( 1 ) ) ;

199

200

201

% Evaluacion de la nueva solucion en la

202

% funcion o bj et iv o

203

f o r o=1:FoodNumber ,

204

ObjValSol ( : , o)= feval ( objfun , Foods ( : , o ) , o ) ;

205

end

206

207

% Evaluacion de los nuevos valores de alimento

208

% en la funcion f i t n e s s

209

f o r b=1:FoodNumber ,

210

Fitness vectorSol ( : , b)= calculateFitness ( ObjVal

211

( : , b ) ) ;

p. 159

end

213

214

% Evaluaciones de los valores fitness , por cada uno

215

% de los alimentos

216

f o r l =1:FoodNumber ,

217

FitnessSol ( : , l )= sqrt (sum ( ( Fitness \

218

vectorSol ( : , l ) )

219

. ˆ 2 ) / NumGeof ) ;

220

end

221

222

% Posicion y valor del minimo en el vector FitnessSol

223

[ GlobalMinSol , BestIndSol ]= min ( FitnessSol ) ; %Llama

224

% el punto evaluado como minimo a traves de

225

% la posicion BestInd

226

227

%Seleccion de esa ” mejor ” abeja

228

GlobalParamsSol=Foods ( : , BestIndSol ) ;

229

230

%Una seleccion prematura es aplicada entre la

231

%solucion actual y una mutante

232

i f ( GlobalMinSol < GlobalMin ) % Si la solucion

233

%del mutante es mejor que la actual , reemplaza

234

%la solucion mutante y r e i n i c i a el contador de

235

%intentos

236

Foods ( : , i )= sol ( : , BestIndSol ) ;

237

Foods test ( : , BestIndSol )= sol ( : , BestIndSol ) ;

238

Fitness ( : , BestIndSol )= FitnessSol ( : , BestIndSol ) ;

239

ObjVal ( : , BestIndSol )= ObjValSol ( : , BestIndSol ) ;

240

t r i a l ( i )=0;

241

else

242

t r i a l ( i )= t r i a l ( i )+1; % Si la solucion no mejora ,

243

% incrementa en contador de intentos t r i a l

244

end ;

245

end ;

246

247

i = i +1;

248

i f ( i ==(FoodNumber)+1)

249

i =1;

250

end ;

251

end ;

252

253

% La mejor fuente de comida es guardada

p. 160

ind=BestIndSol ;

256

i f ( GlobalMinSol < GlobalMin )

257

GlobalMin=GlobalMinSol ;

258

GlobalParams=Foods ( : , ind ) ;

259

end ;

260

261

GlobalMins avg ( r , i t e r )= GlobalMin ;

262

263

264

%% %FASE DE ABEJAS BUSCADORAS % % % %

265

266

%NOTA : En el ABC basico , solo una Scout es

267

%lanzada en cada c i c l o ∗/

268

269

% Determina s i el contador de pruebas ha excedido la cantidad

270

%de intentos a traves del valor de t r i a l s .

271

ind= f i n d ( t r i a l ==max( t r i a l ) ) ; % Tuvo la mayor

272

%cantidad de busquedas .

273

ind=ind ( end ) ;

274

275

i f ( t r i a l ( ind)> l i m i t )

276

t r i a l ( ind )=0;

277

sol =(ub(1)−lb ( 1 ) ) . ∗rand (1 ,D)+ lb ;

278

279

f o r g=1:FoodNumber ,

280

ObjValSol ( : , g)= feval ( objfun , Foods ( : , g ) , g ) ;

281

end

282

283

f o r f =1:FoodNumber ,

284

Fitness vectorSol ( : , f )= calculateFitness ( ObjVal

285

( : , f ) ) ;

286

end

287

288

Foods ( : , ind )= sol ;

289

Fitness ( : , ind )= FitnessSol ( : , ind ) ;

290

ObjVal ( : , ind )= ObjValSol ( : , i n d l ) ;

291

292

end ;

293

294

% Matriz de seguimiento de las mejores

295

% ubicaciones con sus mejores individuos

296

UbicacionMejoresIndividuos=BestIndSol ;

297

MejoresIndividuos ( : , r )=Foods ( : , BestIndSol ) ;

p. 161

i t e r = i t e r +1;

300

301

end % Fi na li z a ABC

302

303

contador = 0;

304

GlobalMins ( r )= GlobalMin ;

305

SlothsMins ( : , r )= BestIndSol ;

306

307

time . horas=toc /3600;

308

time . minutos=toc /60;

309

time . segundos=toc ;

310

311

end ; %end of runs

312

313

314

315

ABC = s t r u c t ( ’ Fecha ’ ,S , . . .

316

’ Algoritmo ’ , ’ABC ’ , . . .

317

’ Arquitectura ’ , ’Osaka ’ , . . .

318

’ Poblacion ’ , NP,

319

’ Dimensiones ’ ,D,

320

’ Fuentes de comida ’ ,FoodNumber ,

321

’ S l o t h s f i n a l e s ’ ,Foods ( : , BestIndSol ) ,

322

’ Minimos por iteracion ’ , GlobalMins ,

323

’ Sloth correspondientes por iteracion ’ , SlothsMins ,

324

’ Ciclos maximos ’ ,maxCycle ,

325

’ Rutinas ’ , runtime ,

326

’ Iteraciones ’ , i t e r ,

327

’ Tiempo de Computo s ’ , time , . . .

328

’ Tipo de Metrica ’ , ’ Correlacion cruzada ’ , . . .

329

’ Alimento ’ ,Fodds , . . .

330

’ Parametros globales ’ , GlobalParams , . . .

331

’ Fitness ’ , Fitness , . . .

332

’ Avance fitness ’ , FitnessSol , . . .

333

’ Soluciones exploradas ’ , sol , . . .

334

’ Observaciones ’ , ’ Pruebas de SA realizadas en Lenovo .

335

C r i t e r i o de parada modificado , se ha u t i l z a d o el proceso de

336

screen en la terminal de Ubuntu .

Sin c r i t e r i o de parada . ’ ) ;

337

338

339

340

save ABC

p. 162

f p r i n t f ( ’ \n minimo ( Promedio de los intentos ) = \n ’ ,mean( GlobalMins , r ) ) ;

343

f p r i n t f ( ’ \n minimo ( Puntos ) = \n ’ , [ Foods ( : , BestInd ) ] ) ; Listing 3. ABC.m

p. 163

ANEXO J. FUNCI ´ON REALIMENTACI ´ON

function fFitness = calculateFitness ( fObjV )

fFitness =zeros ( size ( fObjV ) ) ;

ind= f i n d ( fObjV >=0);

fFitness ( ind ) = 1 . / ( fObjV ( ind ) + 1 ) ;

ind= f i n d ( fObjV <0);

fFitness ( ind )=1+abs ( fObjV ( ind ) ) ; Listing 4. Fitness

p. 164

ANEXO K. COMPLEMENTO DE RESULTADOS

En este anexo se presenta un complemento de las pruebas desarrolladas en el capitulo 4. Se realizan cuatro pruebas denominadas como: prueba a, prueba b, prueba c y prueba d. El par´ametro del n´umero de intentos es el ´unico que presenta un cambio en cada una de las pruebas, y el tiempo de c´omputo es tomado como un dato relevante en cada una de ellas (Tabla K 1). La probabilidad ﬁnal (0,001), la temperatura inicial (2,8037), la temperatura ﬁnal (0,1448) y el paso de reducci´on (0,3) son los par´ametros est´aticos en cada prueba. En las Figuras K 1 y K 2 se presentan los valores promedio de la correlaci´on cruzada y de sloths de la prueba a. Tabla K 1. Par´ametros de la tercera prueba con SA.

Par´ametros Prueba a Prueba b Prueba c Prueba d Intentos

Probabilidad inicial

0,7

0,5

o,5 Tiempo de c´omputo [s]

64,23

196,06

611,75

147,44

Figura K 1. Promedio de sloths en la tercera prueba con SA..

p. 165

Figura K 2. Promedio de los coeﬁcientes de la correlaci´on cruzada con respecto al n´umero de ge´ofonos en la tercera prueba con SA.

Para la prueba b se mantienen los valores de iteraci´on e intentos anteriores, pero se var´ıa el valor de probabilidad de encontrar un nuevo punto de b´usqueda aceptable. En las Figuras K 3 y K 4 se presentan los valores promedio de la correlaci´on cruzada y de los sloths.

Debido a que la probabilidad de aceptaci´on disminuye, se ve reﬂejado en la cercan´ıa de los valores ﬁnales encontrados, ya que al menos cuatro de los cinco intentos realizados convergen a valores cercanos de la soluci´on de la capa. Figura K 3. Promedio de sloths en la cuarta prueba con SA.

p. 166

Figura K 4. Promedio de los coeﬁcientes de la correlaci´on cruzada con respecto al n´umero de ge´ofonos en la cuarta prueba con SA.

Para la prueba c se mantienen los par´ametros anteriores a excepci´on de la cantidad de intentos, el cual toma ahora un valor de diez y que arroja un comportamiento similar a la prueba anterior, teniendo al menos ocho de diez intentos que convergen muy cerca al valor te´orico de la capa. En las Figuras K 5 y K 6 se presentan los valores promedio de la correlaci´on cruzada y de los sloths. Figura K 5. Promedio de sloths en la quinta prueba con SA.

p. 167

Figura K 6. Promedio de los coeﬁcientes de la correlaci´on cruzada con respecto al n´umero de ge´ofonos en la quinta prueba con SA.

Ahora bien, estas pruebas se realizaron partiendo de la soluci´on exacta de las capas del modelo generado como gu´ıa, pero para la prueba d, se dejan los primeros dos valores de las capas aleatorias y las otras dos ﬁjas. En las Figuras K 7, K 8 y K 9 se presentan los valores promedio de la correlaci´on cruzada y de los sloths. Para la primera capa (Figura K 7), se tiene un comportamiento ciertamente aleatorio pero tiene una oscilaci´on muy cercana al valor exacto de ´esta capa que es 0,77; pero puede decirse que es exploratorio pero no converge r´apidamente el algoritmo bajo esas condiciones.

Ya para la segunda capa se tiene un comportamiento distinto (Figura K 8), ´este si presenta valores m´as cercanos al real que es 0,5, pero a pesar de que parten de un valor cercano a ´el, despu´es de ciertas iteraciones llegan a un valor muy cercano, siendo as´ı que para la segunda capa se tiene un mejor comportamiento en comparaci´on con la primera capa analizada anteriormente, pero que presentaban avances muy variables.

p. 168

Figura K 7. Promedio de sloths de la primera capa en la sexta prueba con SA. Figura K 8. Promedio de sloths de la segunda capa en la sexta prueba con SA. El an´alisis de la gr´aﬁca (Figura K 9) de correlaci´on para esta prueba es muy particular, ya que, para el avance, no presenta valores cercanos a 1, pero si se puede inferir que ese valor oscilatorio de la primera capa es compensado por el valor de la segunda capa, que ciclo a ciclo, no deja que ´este valor de capa cambie signiﬁcativamente, sino que se tenga una tendencia a situarse alrededor de un valor cercano a 0,8.

Pero el comportamiento del valor de correlaci´on se comporta de manera extra˜na, ya que a pesar de que algunos valores de la capa son muy cercanos, el valor de correlaci´on se estanca en un aparente m´ınimo local y alrededor de ese valor, o sea que existe una tendencia ya que se sit´uan al menos seis de los diez. Aunque no es

p. 169

tan extra˜no este comportamiento, ya que si se interpreta lo que representa cada una de estas l´ıneas, es el error medio cuadr´atico, teniendo as´ı que, son la ra´ız cuadrada de suma de cada uno los puntos al cuadrado; debido a esa variabilidad, el valor que pudo entregar debi´o ser un valor cercano a uno, que es lo deseado te´oricamente, pero sin tener en cuenta qu´e tanto estaba lejos la respuesta del valor deseado. Figura K 9. Promedio de los coeﬁcientes de la correlaci´on cruzada con respecto al n´umero de ge´ofonos en la sexta prueba con SA.

Se observa un comportamiento muy similar a las primeras pruebas realizadas, por esto se determina que el par´ametro que no se var´ıa, que es la probabilidad, puede considerarse como un valor est´andar a la hora de realizar futuras pruebas del algoritmo.

En cuanto a lo que corresponde a la correlaci´on, se evidencia un comportamiento similar a la correlaci´on anterior pero con mayor cantidad de l´ıneas ya que se aumentaron los intentos, teniendo entonces que, la correlaci´on tiene una tendencia a alejarse del valor ideal pero debido al valor entregado por el error medio cuadr´atico, puede ser cercano a 1 pero seguir avanzando sin tener en cuenta el valor de la capa.

Cita: Beleño Dávila, Jair Abraham, Escalante Rojas, Libardo Andrés, Sarmiento López, Yudy Andrea (2015), Uso de técnicas de optimización global, para resolver problemas de inversión sísmica, segundo grupo de exploración, Universidad Industrial de Santander, p. N. https://noesis.uis.edu.co/handle/20.500.14071/32559