WO2005019779A1

WO2005019779A1 - Método de visualización de flujos multifásicos usando tomografía de capacitancia eléctrica

Info

Publication number: WO2005019779A1
Application number: PCT/MX2003/000067
Authority: WO
Inventors: José Carlos ORTIZ ALEMÁN; Roland Martin; José Carlos GAMIO ROFFÉ
Original assignee: Instituto Mexicano Del Petróleo
Priority date: 2003-08-22
Filing date: 2003-08-22
Publication date: 2005-03-03
Also published as: US7496450B2; US20070133746A1; AU2003256141A1

Abstract

La presente invención se refiere a una nueva técnica de reconstrucción de imágenes para la visualización de flujos multifásicos usando tomografía de capacitancia eléctrica (TEC), basada en métodos heurísticos nolineales de optimización global usando recocido simulado y algoritmos genéticos. Se obtienen datos de capacitancia eléctrica medidos entre electrodos colocados sobre la superficie exterior de un tubo(eléctricamente aislante) que contiene los fluidos. Dichos datos dependen de la distribución de los fluidos dentro del tubo. Los datos se procesan para reconstruir una imagen de la distribución espacial de la permitividad eléctrica relativa (también llamada constante dieléctrica) dentro del tubo, pozo o tanque, la cual refleja la distribución de las distintas fases presentes en el flujo.

Description

MÉTODO DE VISUALIZACIÓN DE FLUJOS MULTIFÁSICOS USANDO TOMOGRAFÍA DE CAPACITANCIA ELÉCTRICA

DESCRIPCIÓN

CAMPO TÉCNICO DE LA INVENCIÓN La presente invención se relaciona con un nuevo método de reconstrucción de imágenes para la visualización de flujos multifásicos usando tomografía de capacitancia eléctrica (TCE), para emplearse en medidores tomográficos de flujo multifásico en la cuantificación de los diferentes fluidos (como gas, aceite y agua) producidos por un pozo petrolífero. También se usa para la visualización y monitoreo de otros flujos y procesos multifásicos que ocurren en la industria tales como: optimización del diseño y operación de separadores, optimización del diseño de unidades de desintegración de fluidización catalítica (FCC) y sistemas de lecho fluidizado, y optimización del diseño de reactores de lecho fijo.

ANTECEDENTES DE LA INVENCIÓN En términos generales, la tomografía sirve para obtener una imagen del corte transversal de un objeto en un plano determinado. La tomografía de rayos X fue la primera en ser desarrollada (en la década de 1970) y actualmente se usa de manera rutinaria en el área médica, así como en algunas aplicaciones industriales (por ejemplo, inspección interna de piezas y detección de fallas en materiales). Posteriormente, se han venido desarrollando nuevos métodos de tomografía orientados más bien a procesos industriales, conocidos con el nombre genérico de tomografía de procesos (Williams y Beck, 1995). La finalidad de estos métodos, que comenzaron a evolucionar a mediados de la década de 1980, es obtener una imagen de la distribución espacial de las fases o componentes en un proceso industrial, usando únicamente sensores externos y sin causar ninguna perturbación en el proceso, tal como se indica en la figura 1 , que muestra un sistema de tomografía de procesos compuesto de un tanque o tubería (A), un sistema de adquisición de datos (B) y una computadora de reconstrucción de imágenes (C). En otras palabras, la tomografía de procesos proporciona una manera de 'mirar' desde afuera al interior de los mismos, sin necesidad de alterarlos físicamente, y esto representa un modo radicalmente distinto de obtener información estructural del proceso a nivel global, a diferencia de los métodos tradicionales, que se basan en el muestreo local en un determinado número de puntos. El proceso puede ocurrir en un reactor o mezcladora, en un lecho fluidizado, en el interior de un separador, o dentro de una tubería transportando flujo multifásico, por poner algunos ejemplos. Existe toda una gama de principios y técnicas que pueden ser explotados en la tomografía de procesos, incluyendo métodos eléctricos, basados en la medición de impedancia, ultrasónicos, de resonancia magnética, ópticos, y a base de radiación ionizante (rayos X o gama). En general, los métodos de radiación ionizante producen las imágenes de mejor definición, pero son relativamente lentos. Por otro lado, los métodos eléctricos producen imágenes de baja resolución pero son mucho más rápidos, de construcción robusta y de costo relativamente bajo. De ahí la necesidad de contar con mejores métodos de reconstrucción de imágenes que, tales como los de la presente invención, permitan obtener imágenes con una mayor definición y fidelidad. Particularmente en el caso de la tomografía de impedancia eléctrica o, simplemente, tomografía eléctrica, se ha dado un progreso muy notable en los últimos años. Este tipo de tomografía tiene dos modalidades principales: la tomografía de capacitancia y la tomografía de resistencia. En un sistema de tomografía de capacitancia (Beck et al., 1997; Gamio, 1997; Plaskowski et al., 1995), (tal como se indica en la figura 2, formado por un sensor (A), un sistema de adquisición de datos (B) y una computadora de reconstrucción de imágenes (C)), normalmente usado con mezclas donde la fase continua es aislante, se utiliza un sensor, tal como se indican en las figuras 3 y 4, formado por un arreglo circular de electrodos (B) distribuidos alrededor de la sección transversal a examinar, y se mide la capacitancia eléctrica de todos los posibles pares de electrodos. Por medio de una computadora y un algoritmo de reconstrucción de imágenes apropiado, los datos medidos son usados para producir un mapa (o sea, una imagen) de las variaciones de la constante dieléctrica (o permitividad relativa) dentro del sensor, proporcionando así una indicación gráfica de la distribución interna de los diferentes componentes de la mezcla. Con relación a las figuras 3 y 4, los electrodos (B) pueden estar colocados en el exterior de un tubo aislante (A), a fin de simplificar la fabricación y evitar el contacto directo con los fluidos (E) del proceso. Un segundo tubo metálico externo eléctricamente aterrizado (C) sirve como blindaje o pantalla eléctrica y para proporcionar resistencia mecánica. El sensor también cuenta con dos electrodos de guarda cilindricos (D) que se mantienen asimismo aterrizados. La patente británica GB 2 2 4640, del 6 de septiembre de 989, describe un sistema de tomografía de capacitancia eléctrica en el cual el método empleado para realizar la reconstrucción de imágenes es un algoritmo de proyección inversa lineal, conocido en inglés como 'linear back-projection' (o 'LBP'). Sin embargo, dicho método de reconstrucción produce imágenes de calidad relativamente baja. Los métodos de reconstrucción de imágenes propuestos en la presente invención permiten superar esta seria limitación. La tomografía de resistencia, por otro lado, está orientada a mezclas donde la fase continua es conductora (Plaskowski et al., 1995; Williams y Beck, 1995). En este caso, los electrodos se colocan al ras de la superficie interior del recipiente que contiene la mezcla, y en contacto directo con los fluidos de la misma. Se aplican una serie de corrientes de excitación y se miden los voltajes producidos, y a partir de ellos se encuentra la distribución de conductividad dentro del sensor, la cual indica la distribución de los diferentes componentes. En principio, la tomografía de capacitancia eléctrica (TCE) tiene importantes aplicaciones en la medición de flujo multifásico, y en particular de flujo bifásico gas-aceite, tal como se presenta comúnmente a la salida de muchos pozos petrolíferos. La manera tradicional de cuantificar los diferentes fluidos producidos por un pozo consiste en separar la mezcla de los mismos por gravedad usando grandes tanques, para posteriormente medir cada componente en forma individual usando medidores de flujo homogéneo convencionales. En las últimas dos décadas han surgido medidores de flujo multifásico que permiten cuantificar la producción de un pozo sin necesidad de separar la mezcla (Thorn et al., 1997). Sin embargo, los medidores multifásicos disponibles actualmente sufren de una sensibilidad a cambios en el patrón de flujo, a menos que se utilicen con dispositivos mezcladores o acondicionadores de flujo, lo cual introduce pérdidas irrecuperables de presión (que finalmente se traducen en pérdidas de energía). La limitación anterior se podría superar mediante la TCE, ya que ésta permite conocer el patrón de flujo y ese conocimiento se podría utilizar para compensar la respuesta de los medidores multifásicos convencionales, o, alternativamente, también es posible en principio diseñar un nuevo tipo de medidor multifásico tomográfico, basado en el análisis de las imágenes de TCE obtenidas en dos planos transversales de la tubería próximos entre sí (Hammer eí al., 1997; Plaskowski et al., 1995). Adicionalmente, la TCE tiene aplicaciones potenciales en la visualización, monitoreo y control de numerosos procesos multifásicos industriales. Sin embargo, el factor principal que hasta ahora ha limitado la aplicación práctica de la TCE es la poca fidelidad o exactitud de las imágenes obtenidas usando los métodos de reconstrucción de imágenes disponibles, de ahí la importancia de contar con métodos mejorados como el de la presente invención (Yang y Peng, 2003). Los métodos directos simples como el de proyección inversa lineal (LBP por sus siglas en inglés) producen imágenes relativamente deficientes que sólo proporcionan una indicación cualitativa de la distribución de componentes dentro del sensor. Tal es el caso de los métodos que se describen en la patente británica GB 2214640 del 6 de septiembre de 1989. Por otro lado, los métodos más sofisticados, basados en técnicas iterativas de optimización local, por lo general necesitan de uno o varios parámetros de regularización cuyo valor óptimo depende precisamente de la imagen que se pretende reconstruir (y que no se conoce), además de que la regularización introducida tiene el efecto de suavizar los contornos de la imagen, haciéndola más difusa. Tal es el caso de los métodos que se describen en la patente británica GB 2 329 476 del 24 de marzo de 1999. Existe pues la necesidad urgente de contar con métodos de reconstrucción mejores y con un menor error. Como un ejemplo de ello, recientemente fue patentado un método de reconstrucción de imágenes basado en la aplicación de redes neuronales artificiales (patente norteamericana US 6,577,700 B1 del 10 de junio de 2003). En la presente invención, se describen nuevos métodos de reconstrucción de imágenes basados en recocido simulado y algoritmos genéticos. De acuerdo con las leyes de la física, en especial de la electrostática, el sensor usado en TCE (tal como se observa en la figura 4) puede ser considerado como un caso particular de un sistema de conductores en un medio dieléctrico, cuya teoría fue desarrollada originalmente por J. C. Maxwell (1873). En nuestro caso, los electrodos detectores (B) actúan como los conductores mientras que el tubo aislante (A) y el contenido del sensor actúan como el medio dieléctrico. Para un sensor de n electrodos, la carga eléctrica q¡ inducida en los electrodos está relacionadas con el potencial eléctrico v¡ de los mismos por medio del siguiente conjunto de ecuaciones lineales simultáneas:

#2 = c₂ v_λ + ^C22^V2 + + ^C2n^Vn (1 ) <ln = ^C»lVl + ^C„₂V₂ + + c mi V n

donde los cu son los coeficientes de capacitancia propia (o simplemente la capacitancia propia) del electrodo i, mientras los demás, cy , con / ≠ j, son los coeficientes de capacitancia mutua (o simplemente las capacitancias mutuas) entre los electrodos i y j. Puesta en forma de matrices, la ecuación (1 ) queda como:

q = C v (3)

Dado que las capacitancias mutuas tienen la propiedad de reciprocidad, esto es, c = C_β , sólo hay m = |«(n-l) capacitancias mutuas independientes, correspondientes a la matriz triangular inferior o superior de C, y correspondientes también a cada uno de los m pares distintos de electrodos que es posible formar en el sensor. El valor de las capacitancias mutuas es una función no lineal muy complicada de la geometría del sistema de conductores y de la distribución espacial de la constante dieléctrica o permitividad relativa ε (que llamaremos también simplemente 'permitividad') de los materiales dieléctricos que los separan. En el caso del sensor de TCE, la geometría de los electrodos, la del tubo aislante y el valor de su constante dieléctrica, son todos fijos. Por lo tanto, podemos decir que las capacitancias mutuas son una función no lineal sólo de la distribución espacial de la constante dieléctrica en el interior del sensor, ε(x,y). Al problema de calcular las capacitancias mutuas correspondientes a una distribución específica de permitividad dentro del sensor se le llama problema directo. El empleo de los electrodos de guarda cilindricos (D) mostrados en la figura

3, junto con la suposición de que la distribución espacial de la permitividad cambia poco en el sentido axial, permite representar el sensor por medio de un modelo bidimensional (Xie er a/., 1989). A menos que se indique explícitamente otra cosa, en lo sucesivo utilizaremos el modelo bidimensional del sensor. Por lo tanto, las cargas eléctricas g¡ y las capacitancias c_ü y c de las ecuaciones (1) y (2) serán ahora cantidades por unidad de longitud de los electrodos en sentido axial. Utilizaremos la tilde (o sea, q_¡ , c_¡¡ y c_& ) para denotar las cantidades totales que resultan al considerar la longitud real de los electrodos en la práctica. Las variables anteriores se relacionan entre sí por medio de la longitud de los electrodos L, de acuerdo con

*' ⁼ I ' ^c" ⁼ T ^{y c}* ⁼ T ⁽⁴⁾ Si dividimos el interior del sensor bidimensional en p regiones de igual área (que llamaremos pixeles) donde la permitividad se considera constante, entonces la versión discreta del problema directo es

donde c es el vector de capacitancias mutuas (por unidad de longitud), f son funciones no lineales que no se conocen explícitamente y ε es el vector de permitividades correspondiente a las p regiones o pixeles dentro de la zona de exploración. Aplicando la Ley de Gauss, las capacitancias mutuas por unidad de longitud axial de los electrodos se pueden calcular como

donde ε₀ es una constante física llamada permitividad del vacío, igual a 8.854 10 ^"12 faradios por metro, F¡ es una curva cerrada que rodea al electrodo i, di es un vector (normal) que representa un elemento infinitesimal de la curva r¿ , di es un elemento de longitud de dicha curva, el símbolo ' ^• ' representa el producto escalar de dos vectores, y φ¹ es el potencial electrostático que se produce en el sensor al aplicar un voltaje de V voltios en el electrodo j (al que llamaremos electrodo fuente o emisor) y 0 voltios en todos los demás (a los que llamaremos electrodos receptores). El potencial φi está determinado por la solución de la siguiente ecuación diferencial parcial

V - ε(x,y)Vφ^J = 0 (7)

sujeta a las condiciones de frontera (a) φⁱ= V voltios en el electrodo fuente y (b) φ^J= 0 en los electrodos receptores y en la pantalla externa. En general, la ecuación (7) no tiene solución analítica y es necesario resolverla numéricamente. Al problema de estimar cuál es la distribución espacial de permitividad dentro del sensor que corresponde a un conjunto determinado de valores de capacitancia mutua, se le conoce como problema inverso, y es el problema que deben resolver los métodos de reconstrucción de imágenes. Normalmente, la estimación de la permitividad se hace en forma discreta, representándola por un vector ε como el de la ecuación (5), que se debe calcular a partir de un vector de capacitancias mutuas observadas c, obtenido por medio de un aparato de medición apropiado. Para resolver el problema inverso, la mayoría de los sistemas de TCE emplean el método de proyección inversa lineal (LBP) (Plaskowski et al., 1995; Yang y Peng, 2003; Xie et al., 1989, 1992), que se describe a continuación. Como primer paso, para cada uno de los m = in(n-l) pares de electrodos que es posible formar, se calcula un mapa de sensibilidad definido por _{s¡(k) =} ^c*(*) -^c/(«*» _para ¿ = i₉ . . , _Λ (8) ^C¡ (fitll) ~ ^Ci (emp)

donde k es el número de pixel (de 1 ap), cι(k) es la capacitancia medida con el par de electrodos i cuando el área correspondiente al pixel k está llena de un material de alta permitividad y el resto del sensor contiene material de baja permitividad, mientras que Cif_uU) y c^_mp) son las capacitancias del par de electrodos i cuando el interior del sensor está lleno de material de alta y baja permitividad, respectivamente. Generalmente estos mapas de sensibilidad se calculan resolviendo numéricamente la ecuación (7) y aplicando la ecuación (6). Ya teniendo los mapas de sensibilidad, éstos se pueden usar para obtener una imagen de permitividad a partir de cualquier vector de m = n(n-Y) capacitancias mutuas medidas, c, correspondientes a alguna distribución desconocida de material dentro del sensor. Para ello se deben primero normalizar las lecturas de capacitancia medidas de acuerdo con _{λ¡ (}k) = °i - ^CHe_mp _(Q) ^Ci (jull) ~ ^Ci {emp)

donde λ¡ es la capacitancia mutua normalizada correspondiente al par de electrodos i y c¡ es el valor de capacitancia medido en la práctica para ese mismo par de electrodos. La fórmula básica del método LBP calcula un 'nivel de gris¹ g(k) para cada pixel como

ι=l En principio, se supone que ese nivel de gris está relacionado linealmente con la permitividad, con g = l y g = 0 correspondiendo a la permitividad de los materiales de alta y baja permitividad, respectivamente. El método LBP se basa en realizar una aproximación lineal de un problema que, como ya se ha mencionado, en realidad es esencialmente no lineal (Gamio y Ortiz-Aleman, 2003). Por lo tanto, este método de reconstrucción de imágenes provoca errores considerables, los cuales son particularmente graves si las diferencias de permitividad son grandes. Hasta ahora, la principal alternativa al método LBP ha sido el empleo de métodos iterativos que buscan minimizar una función objetivo, usando técnicas de optimización local como el método de Newton-Raphson regularizado o variantes del mismo (Yang y Peng, 2003). Como un ejemplo de estos métodos está el que se usa en el programa de cómputo EIT2D (Vauhkonen et al., 2001 ), desarrollado por investigadores de Finlandia y el Reino Unido. Su método de basa en minimizar, con respecto a ε, el funcional || c_meαs - c_cα/c || ² + ² || L ε || ² (11 )

donde α es un parámetro de regularización, L es una matriz que contiene algún tipo de información a príorí sobre la suavidad de c, y c_cα¡_c= f(ε) es el vector de las capacitancias mutuas calculadas para una distribución de permitividad específica dentro del sensor. En el programa EIT2D, c_cafc se calcula resolviendo el problema directo por medio del método de elementos finitos. Partiendo de una estimación inicial ε₀, la minimización se realiza por medio del siguiente procedimiento iterativo (básicamente un método de tipo Newton con regularización de Tikhonov) ε_A+ι = H + [ J/ J_t + «^{2 T}L ]^_1 { j,^τ [ c_meαs - f(ε,) ] - Λ^τL ε_t } (12) donde J_k se conoce como la matriz jacobiana, que contiene las derivadas parciales de f(ε), evaluadas en ε_k

J* = J(⁸ _t) = ( i = \ ...m, j = \ ...p ) (13) dε_j ^εt, Sin embargo, estos métodos de reconstrucción de imágenes tienen el problema de que requieren, para su buena operación, de uno o más parámetros de regularización cuyo valor correcto tiene una fuerte dependencia, precisamente, de las características de la distribución de permitividad que se pretende reconstruir, lo cual implicaría conocer de antemano la solución del problema a resolver. Asimismo, generalmente estos métodos producen imágenes distorsionadas, debido a que la regularización tiene un efecto de suavizamiento excesivo sobre la permitividad obtenida. Si la regularización es muy fuerte ocurrirá el suavizamiento mencionado, y si es muy débil el método puede volverse inestables y/o no converger a la solución buscada. Estos algoritmos de optimización local, durante su búsqueda, exploran sólo un pequeño sector del dominio de soluciones, restringido a la vecindad que circunda la solución inicial. Si la solución óptima del problema, es decir, el mínimo absoluto de la función objetivo, se encuentra alejada de la solución inicial, difícilmente será alcanzada debido a la presencia de mínimos relativos interpuestos en su camino, lugares donde suelen quedar entrampados estos métodos. Los métodos más utilizados correspondientes a esta categoría son la inversión lineal por mínimos cuadrados y las técnicas que emplean el gradiente de la función objetivo, como el de máxima pendiente (en inglés 'steepest descent') y el gradiente conjugado. En general, los métodos de búsqueda local explotan la escasa información derivada de la comparación de una pequeña cantidad de modelos, evitando así una búsqueda extensiva en el espacio de modelos (Sambridge y Drijkoningen, 1992). Los métodos de búsqueda global, como su nombre lo indica, exploran todo el dominio de soluciones a lo largo del proceso de inversión. Hacen un rastreo exhaustivo en el espacio de modelos. De esta manera, a pesar de existir soluciones parciales del problema, es mayor la probabilidad de que la solución final corresponda al mejor ajuste entre los datos observados y los datos sintéticos. Esta clase de métodos, en contraste con las técnicas locales, no requiere de la información proporcionada por las derivadas de la función objetivo debido a que en ellos el problema no requiere ser linealizado. Los algoritmos de optimización global emplean criterios estocásticos para explorar simultáneamente todo el espacio de soluciones en la búsqueda del modelo óptimo. El más conocido de los métodos globales es el de Monte Cario, que realiza una búsqueda puramente aleatoria y no sesgada. Es decir, al generar cada nuevo modelo, no aprovecha la información obtenida de los modelos previamente evaluados (Gallagher et al., 1991 ). La aleatoriedad no encauzada es el rasgo más característico de este método, que lo distingue del resto de los métodos globales. Entre las técnicas de optimización global, se encuentran también los métodos de algoritmos genéticos y de recocido simulado. Ambos fueron concebidos como analogías con sistemas de optimización que ocurren en la naturaleza. Los algoritmos genéticos emulan los mecanismos de la evolución biológica mientras que la técnica de recocido simulado tiene su base en la termodinámica. Ambos métodos son intrínsecamente no lineales y, por lo tanto, se prestan naturalmente a su aplicación a la tomografía de capacitancia, un problema no lineal. ESPECIFICACIÓN DE LA INVENCIÓN La presente invención se refiere a un método de reconstrucción de imágenes para la visualización de flujos multifásicos usando tomografía de capacitancia eléctrica, el cual está basado en técnicas heurísticas no lineales de optimización global, en particular al método de recocido simulado y los algoritmos genéticos. Asimismo, considera el arreglo circular de electrodos metálicos rectangulares contiguos alrededor de la pared externa de un tubo hecho de un material eléctricamente aislante, formando un sensor. A través este sensor fluye una mezcla de fluidos en forma de flujo multifásico, cuya distribución espacial en el interior del sensor se desea conocer. Se obtienen datos a base de efectuar mediciones de capacitancia mutua entre todos los electrodos de pared que es posible formar. Dichos datos dependen de la distribución de los fluidos dentro de la tubería o tanque. Por lo tanto, un objeto de la presente invención es proporcionar un método para procesar datos a base de efectuar mediciones de capacitancia mutua entre todos los electrodos de pared para reconstruir una imagen de la distribución espacial de las fases o componentes de la mezcla multifásica, que fluye a través del sensor, utilizando los métodos de recocido simulado y/o de algoritmos genéticos. Otro objeto de la presente invención es que se utiliza en medidores tomográficos de flujo multifásico para cuantificar los diferentes fluidos (como gas, aceite y agua) producidos en un pozo petrolífero. Aún otro objeto más de la presente invención es que se usa para la visualización y monitoreo de otros flujos y procesos multifásicos que ocurren en la industria, tales como, optimización del diseño y operación de separadores, optimización del diseño de unidades de desintegración de fluidización catalítica (FCC) y sistemas de lecho fluidizado, y optimización del diseño de reactores de lecho fijo.

Referencias Bibliográficas Beck M. S., Byars M., Dyakowski T., Waterfall R., He R., Wang S. M. y Yang W. Q. 1997, Principies and industrial applications of electrical capacitance tomography, Measurement + Control, Vol. 30, pp. 197-200. Cruz-Atienza V. M. 1999, Inversión global con algoritmos genéticos y cristalización simulada aplicada a funciones de receptor: modelos estructurales de velocidades para la corteza en la República Mexicana. Tesis, Facultad de Ingeniería, UNAM. Gallagher K., Sambridge M. y Drijkoningen G. 1991 , Genetic algorithms: an evolution from Monte Cario Methods for strongly non-linear geophysical optimization problems, Geophys. Res. Lett, 18, pp. 2177-2180. Goldberg D. E. 1989, Genetic Algorithms in: Search, Optimization, and Machine Learnlng, Addison-Wesley, Reading, MA. Gamio J. C. 1997, A High-sensitivity Flexible-excitation Electrical

Capacitance Tomography System, PhD Thesis, University of Manchester Institute of Science and Technology, Inglaterra. Gamio J. C. y Ortiz-Aleman J. C. 2003, An ¡nterpretation of the linear back- projection algorithm used in electrical capacitance tomography, 3rd World Congress on Industrial Process Tomography, Banff, Canadá (en prensa). Hammer E. A. y Johansen G. A. 1997, Process tomography in the oil industry: state of the art and future possibilities, Measurement + Control, Vol. 30, pp.212-216. Holland J. H. 1975, Adaptation in Natural and Artificial Systems, University of Michigan Press. Maxwell J. C. 1873, A Treatise on Electricity and Magnetism (Vol. I), Clarendon Press, pp. 88-97. Metrópolis N., Rosenblueth A., Rosenblueth M., Teller A. y Teller E. 1953, Equation of state calculations by fast computing machines, J. Chem. Phys., 21 , pp. 1087-1092. Ortiz-Alemán C, Iglesias-Mendoza A., Cruz-Atienza V. M. 1999, Inversión of site response at México City by using genetic algorithms and simulated annealing, EOS, Transactlons of the American Geophysical Union, 80, 46, F708. Ortiz-Alemán C, Urrutia-Fucugauchi J. y Pilkington M. 2001 , Three- dimensional modeling of aeromagnetic anomalies over the Chicxulub cráter, Lunar and Planetary Science Conference, Proceedings CD Volume, 32, Houston, Texas. Ortiz-Alemán C, Urrutia-Fucugauchi J. e Iglesias-Mendoza A. 2002,

Inversión de la estructura del cráter de chicxulub empleando métodos de inversión global, Revista Geofísica, 57, pp. 59-79. Ortiz-Alemán C. y Urrutia-Fucugauchi J., 2003, Central zone structure and magnetic sources in the Chicxulub cráter as derived from three-dimensional modeling of aeromagnetic anomalies, Earth, Planets and Space (en prensa). Pilkington M. 1997, 3-D magnetic imaging using conjúgate gradients, Geophysics, 62, pp. 1132-1142. Plaskowski A., Beck M. S., Thorn R. y Dyakowski T. 1995, Imaging Industrial Flows: Applications of Electrical Process Tomography, Institute of Physics Publishing, UK. Rodríguez-Zúñiga J. L., Ortiz-Alemán C, Padilla G. y Gaulon R. 1996, Application of genetic algorithms to constrain shallow elastic parameters using in situ ground ¡nclination measurements, Soil Dyn and Earth Eng, 16 (3), pp. 223- 234. Sambridge M. y Drijkoningen G. 1992, Genetic algorithms in seismic waveform inversión, Geophys J. Int., 109, pp. 323-342. Sen M. K. y Stoffa P. L. 1995, Global Optimization Methods in Geophysical

Inversión, Elsevier Science Publishers, Amsterdam, The Netherlands. Stoffa P. L. y Sen M. K. 1991 , Nonlinear multiparameter optimization using genetic algorithms: inversión of plane-wave seismograms, Geophysics, 56, pp. 1794-1810. Thom R., Johansen G. A. y Hammer E. A. 1997, Recent developments in three phase flow measurement, Measurement Science and Technology, Vol. 8, pp. 691-701. Vasudevan K., Wilson W. G. y Ladilaw W. 1991 , Simulated annealing statics computation using an order-based energy function, Geophysics, vol. 56, pp. 1831- 1839. Vauhkonen M., Lionheart W. R. B., Heikkinen L. M., Vauhkonen P. J. y Kaipio J. P. 2001 , A MATLAB package for the EIDORS project to reconstruct two- dimensional EIT images, Phys. Measurement, 22, pp. 107-111. Williams R. A. y Beck M. S. (eds) 1995, Process Tomography - Principies, Techniques and Applications, Butterworth Heinemann. Xie C. G., Plaskowski A. y Beck M. S. 1989, 8-electrode capacitance system for two-component flow identification. Part 1 : Tomographic flow imaging, IEE Proceedings A, 136 (4), pp. 173-183. Xie C. G., Huang S. M., Hoyle B. S., Thom R., Lenn C, Snowden D. y Beck M. S. 1992, Electrical capacitance tomography for flow imaging: System model for development of image reconstruction algorithms and design of primary sensors, lEE Proc.-G, 139 (1), pp. 89-98. Yamanaka H. e Ishida H. 1996, Application of genetic algortihms to an inversión of surface-wave dispersión data, Bulletin of the Seismological Society of America, 36, pp. 436-444. Yang W. Q. y Peng L. 2003, Image reconstruction algorithms for electrical capacitance tomography, Measurement Science and Technology, 14(1 ), pp. R1- R13.

BREVE DESCRIPCIÓN DE LOS DIBUJOS DE LA INVENCIÓN

La figura 1 representa esquemáticamente un sistema de tomografía de procesos utilizado en la presente invención. La figura 2 representa esquemáticamente un sistema de tomografía de capacitancia eléctrica empleado en la presente invención. La figura 3 muestra un dibujo del sensor empleado en tomografía de capacitancia usado en la presente invención. La figura 4 muestra un corte transversal del sensor a la altura de la zona de medición utilizado en la presente invención. La figura 5 presenta un diagrama esquemático del método de recocido simulado empleado en la presente invención. La figura 6 muestra el diagrama de flujo de la realización del método de recocido simulado usado en la presente invención. La figura 7 presenta un diagrama esquemático del método de algoritmos genéticos utilizado en la presente invención. La figura 8 muestra esquemáticamente el proceso de cruza de modelos usado en el método de algoritmos genéticos empleado en la presente invención. La figura 9 muestra esquemáticamente el proceso de mutación de modelos usado en el método de algoritmos genéticos en la presente invención. La figura 10 muestra el diagrama de flujo de la realización del método de algoritmos genéticos utilizado en la presente invención. La figura 11 muestra los resultados obtenidos al reconstruir imágenes de un flujo anular simulado de gas y aceite, usando el método de recocido simulado de la presente invención. La figura 12 muestra los resultados obtenidos al reconstruir imágenes de un flujo estratificado simulado de gas y aceite, usando el método de recocido simulado de la presente invención. La figura 13 muestra los resultados obtenidos al reconstruir imágenes de un flujo simulado de aceite con burbujas de gas, usando el método de recocido simulado de la presente invención. DESCRIPCIÓN DETALLADA DE LA INVENCIÓN Los nuevos procedimientos de reconstrucción de imágenes de la presente invención están basados en métodos heurísticos de optimización global, específicamente el método de recocido simulado (llamado en inglés 'simulated annealing') y los algoritmos genéticos. Una pluralidad de n electrodos metálicos se coloca alrededor de la periferia de una región a visualizar y se efectúan mediciones eléctricas de resistencia o de capacitancia entre ellos. O sea, se obtienen m = ^n(n-l) mediciones o datos. Las mediciones eléctricas son preferiblemente valores de capacitancia, aunque también podrían ser valores de resistencia. Como se muestra en la figura 3, los electrodos (B) pueden ser de forma rectangular y estar colocados sobre la pared externa de un tubo (A) hecho de material eléctricamente aislante, formando así un sensor, el cual contiene un flujo multifásico o multicomponente (E). El objetivo es utilizar dichas mediciones para reconstruir una imagen de la distribución espacial del valor de la constante dieléctrica (o permitividad) en la región a visualizar. Dicha distribución de la permitividad refleja la distribución de las sustancias que ocupan el interior del sensor, donde se encuentra la región a visualizar. Consideraremos que dicha región está dividida enp partes iguales.

Método de Recocido Simulado El método de recocido simulado está basado en una analogía con el proceso termodinámico de la cristalización. Un fluido mineral que se enfría lentamente hasta alcanzar un estado de energía bajo, da lugar a la formación de cristales bien definidos. Si, por el contrario, la sustancia abandona su estado de equilibrio térmico con un enfriamiento súbito o parcial, entonces el cristal resultante tendrá muchos defectos, en caso de que la sustancia no forme un vidrio, caracterizado por el desorden metaestable de sus moléculas. Este concepto es utilizado en el contexto de los métodos de optimización para reconocer configuraciones o modelos potencialmente útiles. Los átomos de cada configuración molecular equivalen a los parámetros del modelo en el problema inverso (o sea, la permitividad en los distintos pixeles de una imagen). La energía del sistema para dicha configuración se relaciona con la función de costo (o de desajuste) asociada al conjunto de parámetros que intervienen en el modelo. En nuestro caso para la presente invención, la energía del sistema se asocia con la siguiente norma L₂: m Σi ) meas ~ ) calcY _{jE = Í2} = _-___ -^ (₁ = 1 , . . , I») (14) í=l donde c(í)_meas son las m capacitancias mutuas medidas y c(i)_caι_c son las calculadas al resolver el problema directo para una distribución de permitividad determinada ε. Partiendo de una distribución inicial de permitividad, el método genera una gama de configuraciones o combinaciones de parámetros considerando una cierta temperatura T para el proceso. Para este propósito se emplea el criterio de Metrópolis et al. (1953), que consiste en desplazar un parámetro, en cada iteración, una distancia aleatoria y pequeña. Este desplazamiento provoca un cambio ΔE en la energía total del sistema. Si ΔE es menor o igual a cero, el desplazamiento del parámetro es aceptado y la configuración resultante se considera como la nueva configuración actualizada. Cuando existe un incremento en la energía del sistema (ΔE es mayor que cero), la probabilidad de aceptación o rechazo para el desplazamiento se determina como: P(AE) = e- ^E/τ (15)

Para saber si es o no admitido el cambio de posición que implica el aumento de la energía del sistema, se elige en forma aleatoria un número entre cero y uno, que se compara con el valor de la probabilidad correspondiente a ΔE. Si es menor dicho número, se admite el desplazamiento y se considera a la nueva configuración como la actualizada; si es mayor, no se admite y se conserva la configuración que se tenía antes del movimiento. Repitiendo sucesivamente este procedimiento se está simulando el movimiento térmico de los átomos de un sistema (que se encuentra en equilibrio térmico), a una temperatura fija T. Si lo que se desea es alcanzar el estado base del sistema, es decir, el estado de menor energía y mayor ordenamiento, entonces se debe disminuir muy lentamente la temperatura para simular un proceso cuasiestático. Esto significa que, durante el enfriamiento, el sistema debe experimentar una sucesión de estados infinitesimalmente alejados del estado de equilibrio térmico. El método de recocido simulado tiene tres componentes básicos (Vasudevan, 1991): una función de energía o función de costo, una función de orden (criterio de Metrópolis) y un parámetro que controla la temperatura del sistema. El proceso consta de tres ciclos anidados. En la figura 5 se muestra un diagrama de funcionamiento para el método en la presente invención. El ciclo externo (3) regula la temperatura del sistema. Cada vez que se cumple un ciclo, la temperatura del proceso disminuye al ser multiplicada por un factor R_τ que normalmente es cercano a uno (0 <R_T< 1). De esta manera se lleva a cabo el enfriamiento lento y gradual que se propone. El ciclo intermedio (2) se encarga de actualizar los valores, independientes entre sí, de una serie de constantes K¡ asociadas a cada parámetro. Dichas constantes determinan el máximo incremento que podrá tener cada parámetro al momento de ser perturbado en el ciclo más interno (1) del proceso. Los valores que adoptan estas constantes dependen de la cantidad de veces que haya sido aceptado el modelo actual al término de cada secuencia de ciclos internos (1) según el criterio de Metrópolis. En el ciclo interno (1 ) se perturban los valores de los parámetros empleando los factores K_t , definidos en el ciclo intermedio (2). La perturbación se realiza multiplicando a cada parámetro por el resultado del producto de su correspondiente K¡ por un valor aleatorio entre menos uno y uno. Posteriormente se calcula la respuesta sintética del modelo actual y se evalúa el cambio de energía en el sistema asociado a la nueva configuración de parámetros. Dicha variación de energía corresponde al desajuste entre la curva de datos sintética y la observada o medida. Si el desajuste decrece, entonces la nueva configuración será aceptada como la actual y perturbada de la misma manera. Si, por el contrario, la perturbación aleatoria produce un crecimiento en el desajuste, asociado a un incremento en la energía E del sistema, a dicha configuración se le asigna una probabilidad de aceptación de acuerdo con el criterio de Metrópolis. Los ciclos (1 ), (2), y (3) de la figura 5, se repiten, mientras la temperatura del proceso disminuye progresivamente. Conforme disminuye la temperatura, las variaciones de los parámetros son cada vez más pequeñas. De esta forma, la búsqueda en el dominio de soluciones tiende a confinarse hacia los modelos asociados con el mínimo absoluto de la función de desajuste E. Εl resultado final es un conjunto de valores para los parámetros (o sea, la permitividad en los distintos pixeles que forman la imagen) cuya respuesta sintética reproduce los datos observados, salvo por un error suficientemente pequeño. Como un ejemplo, una realización específica, pero no la única, del método de recocido simulado se presenta en la figura 6. Εn el bloque (1 ), partimos de una serie de valores iniciales de la temperatura, las constantes de perturbación, las permitividades y la función de costo (T₀ , iζ_(o) y £¡- _o) (con i = 1, ..., p) y E₀). Al mismo tiempo, también en el bloque (1) se inicia el contador del ciclo externo, denotado por k. Posteriormente, en el bloque (2) se inician los contadores de los ciclos interno e intermedio, / y j, respectivamente, y entramos al ciclo interno. Como parte de dicho ciclo, en el bloque (3) se perturba aleatoriamente cada uno de los parámetros (o permitividades), uno tras otro. También en el bloque (3), cada vez que se perturba un parámetro, se resuelve el problema directo y se calcula la función de costo o desajuste, E, aplicando la ecuación (14). Si hubo un decremento de E con respecto a la evaluación anterior, se acepta el valor perturbado del parámetro como el nuevo valor del mismo, se incrementa el contador de parámetros i del ciclo interno, y se procede a perturbar el siguiente parámetro (si lo hay). Si, por el contrario, hubo un incremento de E, entonces se aplica el criterio de Metrópolis, en el bloque (4), para decidir si se acepta o no el valor perturbado del parámetro como su nuevo valor actualizado. Si, de acuerdo con dicho criterio, se acepta el nuevo valor, entonces se incrementa el contador del ciclo interno i y se procede a perturbar el siguiente parámetro (si lo hay). Si, de acuerdo con el criterio de Metrópolis, no se acepta el valor perturbado, entonces no se incrementa el contador / y se procede a perturbar nuevamente el mismo parámetro. Una vez que se han actualizado todos los parámetros ε¡ , en el bloque (5) se ajusta el valor de las constantes K_t (que determinan la forma en la cual se perturba a los parámetros en el ciclo interno) y se incrementa el contador del ciclo intermedio, . N_κ determina el número de veces que se repetirá el ciclo interno sin que haya una disminución de la temperatura. O sea, el ciclo intermedio consiste en la repetición del ciclo interno N_κ veces, pero con diferentes valores de K¡ . Esto se hace con el fin de evitar que la temperatura descienda demasiado rápidamente, lo cual puede ser contraproducente en algunas aplicaciones del recocido simulado. Sin embargo, en la presente invención en particular no ocurre así, y el valor de N_κ se toma como uno. Al terminar el ciclo intermedio, en el bloque (6) se reduce la temperatura como se indicó en párrafos anteriores y se incrementa el contador del ciclo externo k. El procedimiento anterior completo se repite hasta que k alcanza el límite de iteraciones N_T, O antes si se alcanza un valor suficientemente bajo de la función de costo E.

La presente invención también se refiere a un método de reconstrucción de imágenes basado en algoritmos genéticos, el cual se describe a continuación.

Método de Algoritmos Genéticos Los algoritmos genéticos, originalmente propuestos por Holland (1975), representan una evolución del método de Monte Cario para problemas de optimización fuertemente no lineales. La búsqueda del modelo óptimo se lleva a cabo explorando simultáneamente la totalidad del espacio de soluciones, empleando una regla de transición probabilística para guiar dicha búsqueda. El proceso se inicia a partir de un conjunto de modelos seleccionados aleatoriamente. Los parámetros de cada modelo se transforman en código binario para formar cadenas denominadas cromosomas, sobre las cuales se aplican criterios de selección natural y genética. Los procesos de selección, cruza y mutación actualizan la población de modelos, dando lugar a una nueva generación de cromosomas, emulando la forma en que los sistemas biológicos evolucionan para producir organismos mejor adaptados al entorno. El proceso completo se repite hasta que la media de la función de ajuste se acerca al máximo ajuste para toda la población. El diagrama de flujo de la figura 7 resume el proceso utilizado para aplicar un esquema de inversión basado en algoritmos genéticos, similar al descrito por Rodríguez-Zúñiga et al. (1996) y Ortiz-Alemán et al. (2002, 2003) para la inversión de parámetros elásticos del subsuelo a partir de datos de inclinación del terreno y para la inversión de datos aero magnéticos, respectivamente. Los pasos básicos del proceso se describen brevemente a continuación, con referencia a la figura 7.

Discretización En el bloque (1 ) de la figura 7, los parámetros se representan por medio de un vector de incógnitas ε, que representa la distribución de permitividad desconocida. La función de costo, que determina el ajuste entre los datos observados y la respuesta sintética de un modelo dado, se denota como E(ε). La codificación de los parámetros se realiza tomando en cuenta la extensión necesaria de la búsqueda en el espacio de modelos y la resolución deseada (Stoffa y Sen, 1991). De esta manera, la extensión se define para cada parámetro estableciendo un par de cotas α_¿ y b¡ , es decir, ai < ε¡ < b¡ . La resolución se controla con el intervalo de discretización d¡ , definido como

donde N¡ es la cantidad de posibles valores para el parámetro durante el proceso (Sambridge y Drijkoningen, 1992). Los modelos permitidos, ε, definidos por el conjunto de parámetros ε¡ , están restringidos al dominio de valores ε_i = a, + jd_¡ para / = 0, ... , N_; (17)

Población inicial También en el bloque (1) de la figura 7, a partir de combinaciones aleatorias de los parámetros, se construye una población inicial de modelos, cuya dimensión depende del problema particular a resolver. Cada combinación se traduce en un conjunto de índices enteros, definidos por la ecuación (17). Estos valores enteros, que establecen el valor particular de cada parámetro del modelo, posteriormente se codifican como cadenas binarias llamadas cromosomas. Estas cadenas están formadas por grupos consecutivos de bits, llamados genes, que representan el valor de los diferentes parámetros ε¡ . Los modelos que forman las poblaciones en generaciones posteriores se crean a partir de los tres mecanismos evolutivos esenciales: la selección, la cruza y la mutación.

Problema directo y evaluación de la función de costo El problema directo, en el bloque (2) de la figura 7, consiste en calcular la respuesta teórica o sintética para cada modelo en el proceso iterativo.

Posteriormente, en el bloque (3) de la figura 7, dicha respuesta sintética se compara con las observaciones o datos mediante alguna medida de semejanza conocida como función de costo. El criterio utilizado en este trabajo es la norma L₂ definida anteriormente en la ecuación (14), aunque puede haber otros. La norma seleccionada para evaluar el ajuste debe ser sensible a la forma y complejidad de las curvas observada y calculada.

Selección En el bloque (4) de la figura 7, a partir de una población de Q individuos y de sus respectivas funciones de costo E(ε_k) (k= l, ...,Q), se determina para cada modelo una probabilidad ('acumulada') de selección P(ε_k) que dependerá de su nivel de ajuste. Una fórmula que puede ser usada para determinar la probabilidad acumulada de selección es:

P(ε_k) = P(ε_H) + ^{E ~}_^{E (8¿)} (18) zí V . max prom

donde E_max, y E_pιvm son las funciones de costo máxima y promedio de la generación, respectivamente, y Q es el número de individuos de la población. Como siguiente paso, se utiliza un procedimiento de ruleta sesgada (Goldberg, 1989) para seleccionar una nueva población de modelos. Se generan Q números aleatorios r_k entre cero y uno. Si P(ε_k.χ) < r_k<P(ε_k), entonces ε_k es seleccionado para formar parte de la nueva población. Esto implica que en la nueva población puede haber modelos 'gemelos'. Adicionalmente, se pueden añadir a la población de Q modelos, 'clones' de los mejores modelos, los cuales no se someterán a la cruza ni a la mutación, con el objeto de asegurar que dichos modelos deseables no se pierdan al pasar por esos procesos que son de naturaleza aleatoria. Alternativamente, la probabilidad de selección puede determinarse como (Sambridge y Drijkoningen, 1992)

P(ε_k) = a - b E(ε_k) (19)

que describe una distribución de probabilidades lineal, y P(ε_k) = Ae^{~BE z} (20)

que corresponde a una distribución exponencial. Los valores que suelen tomar las constantes a, b, A y B son los siguientes b = Q^~l (E - E prom r / ¹ , ? a > ~ bE

donde E_max, E_prom y E_σ son las funciones de costo máxima, promedio, y la desviación estándar de todos los desajustes de la población inicial, respectivamente. Stoffa y Sen (1991) proponen un criterio de selección basado en una probabilidad de actualización. Εl criterio consiste en comparar el desajuste de cada modelo perteneciente a la generación actual con el de un modelo de la generación anterior, elegido al azar. Si el desajuste del modelo actual es menor, entonces éste se conserva. De lo contrario, se considera un valor P_u que establece la probabilidad de sustitución del modelo anterior por el actual. Este procedimiento controla la influencia de los ajustes de generaciones previas sobre la población actual. El valor sugerido por Stoffa y Sen (1991) para la probabilidad de sustitución P_u es del 90%.

Cruza Los nuevos modelos se engendran en el bloque (5) de la figura 7 a partir de una generación progenitora de Q modelos. En forma aleatoria se integran Q/2 parejas de individuos. Cada pareja es potencialmente capaz de cruzarse. Para determinar cuales de ellas llevarán a cabo la cruza, se le asocia a cada una un número aleatorio entre cero y uno. Si dicho número es menor que la probabilidad de cruza P_c entonces la pareja correspondiente efectúa el apareamiento. De lo contrario, la pareja se preserva intacta en la siguiente generación. El mecanismo de cruza se basa en elegir aleatoriamente la posición de un gen para ambas cadenas. Las cadenas se parten en ese punto para intercambiar información entre ellas, tal como se observa en la figura 8, donde se muestra en (A) la partición de los cromosomas padres en cuatro gametos, y en (B) el intercambio y unión de dos gametos formando dos cigotos hijos. El propósito de cruzar dos cadenas diferentes es explorar nuevas regiones del dominio de soluciones, donde pudiera ubicarse el mínimo absoluto. En el proceso normal de cruza, las parejas que se aparean tienen dos hijos, y la población de modelos se mantiene automáticamente en Q individuos. Alternativamente, se puede considerar que las parejas tienen sólo un hijo, en cuyo caso la población se debe completar con modelos generados aleatoriamente hasta volver a tener Q individuos.

Mutación La mutación, al igual que la reproducción sexual (o cruza), propicia la diversidad genética en una población. La mutación hace posible que la búsqueda prospere cuando se encuentra confinada en los alrededores de un mínimo local. La mutación se realiza en el bloque (6) de la figura 7, mediante el cambio de paridad de un bit seleccionado al azar dentro de la cadena binaria (cromosoma). El porcentaje de modelos a ios cuales se aplica el proceso de mutación, al igual que en la cruza, depende de un parámetro denominado probabilidad de mutación, simbolizado por P_m . Este mecanismo previene la convergencia prematura del método, cuando la población es excesivamente homogénea e incapaz de continuar el proceso evolutivo. En la figura 9 se representa una cadena arbitraria y se ejemplifica la mutación del valor de un bit aleatorio, en (A) se muestra el cromosoma original y en (B) el cromosoma mutado. Una alternativa para definir la probabilidad de mutación P_m fue propuesta por Yamanaka e Ishida (1996). Consiste en determinar el nivel de homogeneidad de los individuos en cada generación mediante el cálculo de un coeficiente de variación promedio , para cada parámetro, mediante la expresión

donde ? es la cantidad de parámetros, ε_¡ es el promedio del z^'-ésimo parámetro, y σ¿ es la desviación estándar. A continuación se define P_m como función de /, es decir

P_M para y > 0.1 P_m = < 0Λ para 0.02 < 7 < 0.1 (22) 0.2 para γ < 0.02 donde P_ini es la probabilidad de mutación inicial. Con la mutación concluye la secuencia de operaciones que define a un algoritmo genético. Dicha secuencia se repite hasta satisfacer alguna tolerancia preestablecida. Como ejemplo, una realización específica, pero no la única, del método de algoritmos genéticos se presenta en la figura 10. En el bloque (1), se parte de una población de Q modelos generados aleatoriamente, cada uno representando una distribución de permitividad. Al mismo tiempo, también en el bloque (1) se inicia el contador de iteraciones z. A seguir, se inicia el contador de modelos k y comienza un ciclo por medio del cual, para cada uno de los Q modelos, en el bloque (2) se resuelve el problema directo y se calculan la función de desajuste E(ε_k) y la probabilidad acumulada de selección P(ε_k). Posteriormente, en el bloque (3) se realiza la selección de modelos con menor E, de acuerdo con el valor de P(ε_k) y aplicando el procedimiento de ruleta sesgada ya descrito con anterioridad. Luego, en el bloque (4) se pasa a la cruza de modelos. En dicho bloque (4), se forman aleatoriamente Qfl parejas, que se cruzan selectivamente de acuerdo con el valor de P_c , produciendo sólo un hijo por pareja. En ese mismo bloque (4), las parejas no cruzadas se conservan (ambos integrantes) iguales y la población se completa a Q individuos con nuevos modelos aleatorios. Por último, en el bloque (5) se calculan el coeficiente de variación promedio γ y la probabilidad de mutación P_m , y se realiza la mutación de modelos de acuerdo con la descripción dada en párrafos anteriores. Todo el procedimiento anterior se repite continuamente, incrementando i en cada iteración. El algoritmo se termina cuando se alcanza un número preestablecido de iteraciones N o cuando la función de desajuste es lo suficientemente pequeña.

Formulación del problema directo. El problema directo consiste en el cálculo de las capacitancias mutuas cy , i ≠ j, que resultan de la presencia de una distribución de permitividad ε en el interior del sensor. Tanto el método de recocido simulado como el de algoritmos genéticos requieren de la solución repetida del problema directo. Por lo anterior, es importante contar con un método adecuado para resolver dicho problema, que logre un equilibrio razonable entre exactitud (o precisión) y rapidez. En el contexto de esta invención, el problema directo se puede resolver utilizando una rutina optimizada desarrollada por los autores basada en el método de volúmenes finitos, que es descrita muy brevemente a continuación. Esta rutina es más eficiente que las reportadas hasta hoy en la literatura (Yang y Peng, 2003) pues es comparable en precisión con implantaciones basadas en el método de elementos finitos con mallas de 9,000 elementos triangulares. La velocidad de ejecución es superior a las de los métodos de elementos finitos y diferencias finitas. La rutina está escrita en Fortran 90 y es totalmente transportable (se ha probado en computadoras tipo PC y 'clusters' de PCs basados en Windows y Linux, en estaciones de trabajo tipo SUN y Alpha y en supercomputadoras Cray). La rutina puede extenderse al caso tridimensional sin mayores modificaciones y es fácilmente paralelizable. El problema directo se resuelve mediante el método de volúmenes finitos en una configuración cilindrica. De esta manera se eliminan las soluciones indeterminadas en el centro del disco (como sucede en el método de diferencias finitas) y el refinamiento de la malla se vuelve más flexible en comparación con los métodos de elementos finitos. Se resuelven la ecuación (7), que se repite a continuación

V - ε(x,y)Vφ^k = 0

donde ε es la permitividad y φ^k es el potencial electrostático que se genera cuando el electrodo k es el emisor (o fuente). La ecuación está sujeta a las condiciones de frontera (a) φ^k=V voltios en el electrodo fuente y (b) φ^k = 0 en los electrodos receptores y en la pantalla externa. Definiendo las coordenadas radial y angular como r y θ, y utilizando el método de volúmenes finitos, la ecuación discreta es formulada en forma conservadora en cada celda Ω,₇ como f V - (εVφ^k) dΩ_ϋ = 0 para i = \, . . , N_r y j = l, . . , N_θ (23)

donde i y / se refieren a la discretización en r y θ, respectivamente, y N_r y N_θ son el número de secciones en las que se subdivide el radio y la circunferencia del sensor, respectivamente. Aplicando el teorema de Gauss en coordenadas polares, las ecuaciones discretas pueden escribirse como l e^φ dT_v = 0 (24)

donde Ty es la frontera de la celda de volumen finito Ω . La frontera Ty se define por ry y T_E a lo largo de las coordenadas radiales, y I> y r^ a lo largo de las coordenadas angulares. La ecuación (24) puede expresarse como la suma de los flujos a través de las caras T_N, r_s , T_E y IV

De (25), el término correspondiente a los flujos en el radio cero desaparece y el problema es equivalente a resolver las ecuaciones en la cercanía del centro sobre triángulos que tienen un vértice en el centro. Entonces, el sistema discreto de ecuaciones para el problema directo resulta bien planteado. El sistema completo es similar a un sistema de ecuaciones Laplaciano y debe resolverse un sistema diagonal bandeado que incluye las condiciones de frontera periódicas impuestas por la geometría del problema. La matriz correspondiente es positiva definida y no simétrica, características que se pueden aprovechar para seleccionar específicamente los métodos óptimos de solución. Finalmente, las capacitancias mutuas se calculan integrando los gradientes del potencial a lo largo de una curva cerrada que rodea los electrodos, de acuerdo con la ecuación (6), que se repite a continuación

donde ε₀ es la permitividad del vacío (8.854 x 10^"12 faradios por metro), r¡ es una curva cerrada que rodea al electrodo i, di es un vector (normal) que representa un elemento infinitesimal de la curva r¡, di es un elemento de longitud de dicha curva y φ^k es el potencial electrostático que se produce en el sensor al aplicar un voltaje de V voltios en el electrodo k (fuente) y cero voltios en todos los demás (receptores). La integración se realiza utilizando una regla trapezoidal y los gradientes de potencial se calculan al cuarto orden. Durante el procedimiento de reconstrucción de una imagen de permitividad por el método de recocido simulado, es necesario resolver el problema directo y encontrar el potencial eléctrico repetidamente para sucesivas distribuciones de permitividad relativamente similares, mientras el método converge a la solución final. En virtud de que el potencial correspondiente a dichas distribuciones sucesivas cambia relativamente poco, es posible acelerar todo el proceso si se emplea un método iterativo para resolver el problema directo, tomando como estimación inicial del potencial a la solución correspondiente a la configuración de permitividad anterior. Como la estimación inicial del potencial estará muy cerca de la solución, dicho método iterativo convergerá en pocas iteraciones, alcanzando rápidamente una exactitud aceptable.

EJEMPLO Con objeto de evaluar el desempeño de los métodos de reconstrucción de imágenes descritos anteriormente, calculamos un conjunto de datos TCE sintéticos con la subrutina del problema directo. Para emular las fuentes principales de incertidumbre en los datos producidos por uri instrumento de medición operando en condiciones de trabajo (esto es, los errores aleatorios inherentes al proceso de medición de cualquier variable física y los errores debidos a la precisión propia del instrumento de medición), los cálculos de las capacitancias sintéticas para el modelo ideal fueron evaluadas con una precisión numérica del orden de 10^"11 en el método iterativo utilizado en el problema directo para el cálculo de los potenciales. Con objeto de simular la imprecisión (sistemática) asociada al sensor TCE, durante el proceso de inversión (estimación de la distribución de permitividad eléctrica en el interior del tubo) se utilizó una precisión sensiblemente menor (10^"5) en el cálculo de los potenciales en el problema directo. Las interpretaciones sin restricciones de datos de campos potenciales pueden ser de muy poco interés práctico debido a la gran ambigüedad presente entre las observaciones y las soluciones estimadas. La no unicidad en los problemas con campos potenciales surge principalmente de dos fuentes: la primera es la ambigüedad inherente provocada por la física del problema que permite la existencia de una infinidad de soluciones que reproducen la anomalía de campo potencial; la segunda resulta de la utilización de un número finito de datos que están contaminados por errores y que pueden contener información insuficiente para construir una solución única del problema. Las estrategias que permiten superar esta no unicidad consisten en la incorporación de suficiente información a priori para restringir las soluciones resultantes a una región del espacio de parámetros que es considerada físicamente razonable (Pilkington, 1997). En el caso particular de la tomografía de capacitancia eléctrica se cuenta con información acerca de los valores típicos de la permitividad eléctrica para el gas, aceite y agua. El conocimiento de estos valores, así como el contraste significativo entre las propiedades del gas y el aceite, permiten aplicar las técnicas de inversión que aquí se discuten a datos de TCE para flujos bifásicos contaminados por errores relativos menores al 2% (siendo este el nivel máximo de error que producen los sistemas de adquisición de datos empleados normalmente en TCE). Por otra parte, si se cuenta con una estimación estadística precisa de las ¡ncertidumbres en los datos, entonces es posible construir un modelo que considere las ¡ncertidumbres en los datos y en las estimaciones de los parámetros, utilizando el esquema que proponen diversos autores para la inversión de datos de campos potenciales (por ejemplo, Sen y Stoffa, 1995). Usando la subrutina del problema directo a base de volúmenes finitos, se calcularon las capacitancias sintéticas considerando tres distribuciones específicas de permitividad, correspondientes a tres patrones de flujo típicos (anular, estratificado y de burbujas). Simulamos un sensor TCE de 12 electrodos y los valores de capacitancia para todas las combinaciones posibles de electrodos fueron calculados. Esos fueron los datos simulados. Consideramos una distribución de dos componentes con un material de permitividad baja de 1 (gas) y otro de permitividad alta de 2.5 (aceite). Nuestras pruebas numéricas se han hecho con datos libres de ruido y con datos contaminados con errores aleatorios de hasta 1%. El algoritmo está codificado en Fortran 90 y corre en una computadora Pentium 4 de 1.7 gigahertz con 512 megabytes de memoria RAM. Probamos con una malla de 120x60 elementos para reducir los tiempos de la inversión (~ 30 minutos para 60,000 iteraciones), pero los resultados son válidos para dimensiones más grandes. La validez de los resultados de estas simulaciones fue corroborada con resultados experimentales de laboratorio obtenidos empleando modelos físicos y un sistema real de tomografía de capacitancia. Después de una parametrización adecuada, tanto el método de algoritmos genéticos como el de recocido simulado produjeron resultados muy similares, para los tres patrones de flujo considerados. La calidad de las imágenes de permitividad reconstruidas depende principalmente del número de iteraciones del método, tal y como sucede para muchas otras aplicaciones (por ejemplo, Ortiz- Alemán et al., 1999, 2001, 2002, 2003; Cruz-Atienza, 1999). En la figura 11 se presenta la imagen ideal o 'verdadera' de un flujo anular simplificado (A), y las reconstrucciones obtenidas después de 30,000 (B) y 60,000 (C) cálculos del problema directo (o iteraciones del método). En las figuras 12 y 13 se presentan figuras similares para un flujo estratificado y un flujo de burbujas, respectivamente, mostrando en (A) la imagen ideal, en (B) la imagen reconstruida con 30,000 iteraciones y en (C) la imagen reconstruida con 60,000 iteraciones. En todos los casos mostrados en las figuras se usó el método de recocido simulado (con los algoritmos genéticos no hubo diferencias significativas en los resultados) y se partió de una distribución de permitividad inicial homogénea. Estos resultados muestran claramente que las imágenes reconstruidas son muy parecidas a las imágenes ideales o 'verdaderas'. La exactitud de estas reconstrucciones es considerablemente mejor que las reportadas hasta ahora en la literatura (Yang y Peng, 2003). Si bien los métodos de esta invención no requieren, para converger, que la estimación inicial de la distribución de permitividad esté próxima a la solución, es posible acelerar un poco el proceso si se emplea como estimación inicial la imagen resultante de un método directo simple como el LBP.

Claims

NOVEDAD DE LA INVENCIÓNHabiendo descrito la presente invención, ésta se considera como novedad y por lo tanto, se reclama como propiedad el contenido en las siguientes cláusulas:

1.- Un método de reconstrucción de imágenes que comprende: (a) la obtención de datos medidos correspondientes a la imagen; (b) el procesamiento de dichos datos medidos por medio de métodos heurísticos no lineales de optimización global, específicamente el método de recocido simulado o los algoritmos genéticos, para obtener una imagen; y (c) el desplegado de dicha imagen en un dispositivo de visualización.

2.- Un método de conformidad con la cláusula 1 , caracterizado porque los datos medidos son parámetros eléctricos en el perímetro de una región tal como el interior de una tubería, el de un pozo o el de un tanque, para obtener una imagen de la distribución espacial de la permitividad eléctrica (o constante dieléctrica) o de la conductividad eléctrica en la que se refleja la distribución espacial de los materiales o sustancias, tales como gases y/o líquidos.

3.- Un método de conformidad con la cláusula 2, caracterizado porque los parámetros eléctricos son valores de capacitancia eléctrica medidos entre los distintos electrodos un sensor formado por una pluralidad de electrodos colocados a lo largo del perímetro de la región (tubería, pozo, tanque).

4.- Un método de conformidad con las cláusulas 2 y 3, caracterizado porque la imagen obtenida es una distribución espacial de la permitividad eléctrica (o constante dieléctrica) en la región de interés (tubería, pozo, tanque), la cual refleja la distribución espacial de los materiales o sustancias, tales como gases y/o líquidos, que ocupan dicha región (tubería, pozo, tanque).

5.- Un método de conformidad con las cláusulas 2 a 4, caracterizado porque la imagen está formada por una cantidad finita de subregiones o pixeles en la visualización, cuyo número depende de la resolución deseada.

6.- Un método de conformidad con la cláusula 3, caracterizado porque el sensor está formado por un tubo fabricado de un material eléctricamente aislante, sobre cuya pared externa se coloca un anillo de electrodos metálicos rectangulares.

7.- Un método de conformidad con las cláusulas 3 y 6, caracterizado porque el sensor contiene un flujo multifásico o multicomponente, y obtiene imágenes 5 que muestran la distribución de las fases o componentes, tales como gases y/o líquidos de acuerdo a lo dicho en antecedentes.

8.- Un método de conformidad con la cláusula 1 , caracterizado porque emplea el método de recocido simulado como el método de optimización.

9.- Un método de conformidad con la cláusula 8, caracterizado porque se0 minimiza iterativamente una función de costo asociada con la energía del sistema, con respecto a ε (el vector de los valores de la permitividad en cada pixel de la imagen), y siendo de la forma:

5 donde c ^βαr son las m capacitancias mutuas medidas y c ^c(ε_k) son las calculadas al resolver el problema directo para una distribución de permitividad determinada ε_k.

10.- Un método de conformidad con la cláusula 9, caracterizado porque para0 efectuar la minimización se emplea el criterio de Metrópolis.

11.- Un método de conformidad con la cláusula 10, caracterizado porque en el proceso iterativo de minimización se emplea como estimación inicial de la distribución de permitividad, ε, tanto una distribución homogénea como la que resulta de aplicar a los datos medidos el método de proyección inversa5 lineal (LBP).

12.- Un método de conformidad con la cláusula 9, caracterizado porque el cálculo de las capacitancias calculadas c ^c(ε_k) , conocido como el problema directo, se realiza por medio del método de volúmenes finitos.

13.- Un método de conformidad con la cláusula 12, caracterizado porque para lao solución del sistema de ecuaciones que resulta al resolver el problema directo (cálculo de c ^c(ε_k) ), se utilizan métodos iterativos que convergen rápidamente al emplear como estimación inicial del potencial electrostático, el potencial que resultó de la solución del problema directo en la iteración anterior del problema inverso.

14.- Un método, de conformidad con la cláusula 1 , caracterizado porque se emplea el método de algoritmos genéticos como el método de optimización.

15.- Un método de conformidad con la cláusula 14, caracterizado porque, a partir de una población inicial de Q modelos de permitividad ε_kip) (k=\,...,Q), aplicando mecanismos evolutivos como selección, cruza y mutación, se obtienen iterativamente nuevas poblaciones.

16.- Un método de conformidad con la cláusula 15, caracterizado porque los individuos ε_/< poseen una función de costo o de desajuste pequeña, donde la función de desajuste está dada por:

donde c'"^eas son las m capacitancias mutuas medidas y c ^c(ε_k) son las calculadas al resolver el problema directo para el modelo (o distribución de permitividad) ε¿.

17.- Un método de acuerdo con las cláusulas 15 y 16, caracterizado porque la probabilidad acumulada de selección para un modelo determinado ε_k está dada por

P(ε_k) = P(ε_k__l) + ^E^ ^{~ E}W Q i.^Emaχ - ^E rom ) donde E_max, y E_prom son las funciones de costo máxima y promedio de la generación, respectivamente, y Q es el número de individuos de la población.

18.- Un método de conformidad con la cláusula 17, caracterizado porque se puede emplear el procedimiento de la ruleta sesgada para decidir cuales modelos son seleccionados en cada iteración.

19.- Un método de conformidad con la cláusula 18, caracterizado porque la cruza y la mutación de modelos se realizan aleatoriamente de acuerdo las probabilidades de cruza y de mutación, P_c yP_m -

20.- Un método de conformidad con la cláusula 19, caracterizado porque la probabilidad de mutación P_m está determinada mediante el coeficiente de variación promedio γ, dado por

donde p es la cantidad de parámetros, ε_t es el promedio del /-ésimo parámetro, y σ, es la desviación estándar.

21.- Un método de conformidad con las cláusulas 19 y 20, caracterizado porque P_m se define como función de γ, es decir

P_m para γ > 0.\ P_m = 0.1 para 0.02 < ;κ < 0.1 0.2 para ; < 0.02 donde P_ιm es la probabilidad de mutación inicial.

22.- Un método de conformidad con las cláusulas 15 y 16, caracterizado porque el cálculo de c™^!c(ε_k) , conocido como el problema directo, se realiza por medio del método de volúmenes finitos.