Reconstrucci
´
on de Im
´
agenes Optoac
´
usticas: Efecto
de los Sensores Piezoel
´
ectricos de Banda Ancha
P. Massaro Rocca
∗
, L. Ciocci Brazzano
∗†
, E. Acosta
∗
, M. G. Gonz
´
alez
∗†1
∗
Universidad de Buenos Aires, Facultad de Ingenier
´
ıa,
Grupo de L
´
aser,
´
Optica de Materiales y Aplicaciones Electromagn
´
eticas (GLOMAE)
Paseo Col
´
on 850, C1063ACV, Buenos Aires, Argentina
†
Consejo Nacional de Investigaciones Cient
´
ıficas y T
´
ecnicas, (CONICET)
Godoy Cruz 2290, C1425FQB, Buenos Aires, Argentina
1
mggonza@fi.uba.ar
Abstract—Optoacoustic imaging is a hybrid technique with
a rapid evolution within the area of biomedical imaging. It
combines the rich and versatile optical contrast with the high
spatial resolution associated with the low-dispersion ultrasonic
wave propagation nature. Its application raises a series of
difficult problems that require further developments. One
of them, the objective of this work, is to analyze how the
optoacoustic image reconstruction algorithms are affected by
the use of non-ideal (real) broadband sensors. In particular,
as an original contribution, we study the case of integrating
ultrasonic sensors based on piezoelectric polymer thin films. In
order to achieve this, a temporary back-projection algorithm
was implemented for line detectors. Then, using a parametric
model for this type of sensors, we studied the influence of
the use of a real device on the efficiency of the implemented
algorithm.
Resumen— La obtenci
´
on de im
´
agenes optoac
´
usticas es
una t
´
ecnica h
´
ıbrida con una r
´
apida evoluci
´
on dentro del
´
area
de im
´
agenes biom
´
edicas.
´
Esta se beneficia tanto del rico y
vers
´
atil contraste
´
optico como de la alta resoluci
´
on espacial
asociada con la naturaleza de baja dispersi
´
on de propagaci
´
on
de las ondas ultras
´
onicas. Su aplicaci
´
on plantea una serie de
dif
´
ıciles problemas que exigen ulteriores desarrollos. Uno de
ellos, el objetivo de este trabajo, es analizar c
´
omo afecta a
los algoritmos de reconstrucci
´
on de im
´
agenes optoac
´
usticas
el uso de sensores no ideales (reales) de banda ancha. En
particular, como aporte original, se estudia el caso de los
sensores ultras
´
onicos extensos basados en pel
´
ıculas delgadas
de pol
´
ımeros piezoel
´
ectricos. Para lograr esto, se implement
´
o
un algoritmo de retroproyecci
´
on temporal para detectores
lineales. Luego, usando un modelo param
´
etrico para este tipo
de sensores, se estudi
´
o la influencia del uso de un dispositivo
real sobre la eficiencia del susodicho algoritmo.
I. INTRODUCCI
´
ON
La tecnolog
´
ıa para la obtenci
´
on de im
´
agenes optoac
´
usticas
(IOA) de origen biol
´
ogico est
´
a actualmente en pleno desa-
rrollo, generando nuevos enfoques t
´
ecnicos y aplicaciones.
En comparaci
´
on con la formaci
´
on de im
´
agenes ecogr
´
aficas
por ultrasonido (IEU), las amplitudes de las se
˜
nales OA son
relativamente bajas pero su contenido espectral es amplio,
abarcando frecuencias desde varias decenas de kHz hasta
un centenar de MHz para estructuras a escala microm
´
etrica.
De hecho, se ha demostrado que es posible generar IOA a
partir de tejido vivo con una resoluci
´
on espacial del orden
de los micrones, lo que permite la obtenci
´
on de informaci
´
on
bioqu
´
ımica relevante [1].
Una configuraci
´
on t
´
ıpica OA consta de tres elementos
esenciales: una fuente de luz, un sistema de detecci
´
on de
las ondas ac
´
usticas y uno para el procesamiento de las
se
˜
nales. Por ejemplo, en el modo de excitaci
´
on pulsada,
el tejido es iluminado por un l
´
aser que emite pulsos de
luz monocrom
´
atica con una duraci
´
on t
´
ıpica de algunos
nanosegundos. Las tasas de repetici
´
on de pulsos son del
orden de unas pocas decenas de Hertz, con energ
´
ıas en
el rango de microJoules por pulso. Cuando esta t
´
ecnica
es utilizada para realizar tomograf
´
ıa (TOA), los perfiles de
presi
´
on generados por la excitaci
´
on
´
optica son captados con
sensores que rodean la zona de inter
´
es. Por lo tanto, se puede
decir que la tecnolog
´
ıa de formaci
´
on de IOA es h
´
ıbrida, lo
que le permite beneficiarse tanto del rico y vers
´
atil contraste
´
optico como de la alta resoluci
´
on espacial (limitada por
difracci
´
on) asociada con la naturaleza de baja dispersi
´
on de
propagaci
´
on de las ondas ultras
´
onicas, en comparaci
´
on con
las electromagn
´
aticas [2].
Es posible clasificar los sensores ultras
´
onicos usados
en la obtenci
´
on de IOA en dos categor
´
ıas: transductores
piezoel
´
ectricos, en los cuales la presi
´
on medida es directa-
mente proporcional a la se
˜
nal el
´
ectrica, y detectores
´
opticos,
los cuales son sensibles a cambios en el largo del camino
´
optico inducido por ondas de presi
´
on [3]. Los transductores
piezoel
´
ectricos son los m
´
as com
´
unmente utilizados y est
´
an
basados en materiales polim
´
ericos (sensores de banda ancha)
o cer
´
amicos (sensores resonantes). Desde el punto de vista
geom
´
etrico, estos sensores son de dos tipos: de peque
˜
na
apertura (puntuales) o de gran apertura (extensos). En los
de peque
˜
na apertura los efectos de difracci
´
on deben ser
compensados por los algoritmos de reconstrucci
´
on de la
imagen y, por otro lado, tienen la ventaja de que es posible
hacer un arreglo con m
´
ultiples elementos. De esta forma
se lleva a cabo una r
´
apida adquisici
´
on de im
´
agenes con
resoluci
´
on adecuada. En el caso de un detector de gran
apertura se puede conseguir una imagen de alta resoluci
´
on
angular y, si el detector tiene un gran ancho de banda, se
tiene una buena discriminaci
´
on en distancias [4].
Para conseguir im
´
agenes a partir de se
˜
nales OA se deben
resolver dos problemas inversos: uno ac
´
ustico y otro
´
optico
[5]. En ambos casos se parte de las se
˜
nales ac
´
usticas
medidas. En el problema inverso ac
´
ustico se intenta mapear
la energ
´
ıa depositada en la muestra, mientras que el objetivo
del problema inverso
´
optico es conseguir la imagen del
coeficiente de absorci
´
on. La TOA aplicada a la obtenci
´
on
Revista elektron, Vol. 1, No. 2, pp. 58-65 (2017)
ISSN 2525-0159
58
Recibido: 25/09/17; Aceptado: 13/11/17
de im
´
agenes de objetos vivos es la aplicaci
´
on que mayores
desaf
´
ıos presenta para resolver ambos problemas inversos.
´
Opticamente, las grandes variaciones en los coeficientes
de dispersi
´
on y absorci
´
on de los tejidos vivos lleva a
problemas inversos no lineales muy complejos. Por otro
lado, ac
´
usticamente, la geometr
´
ıa del sistema de detecci
´
on,
as
´
ı como tambi
´
en la heterogeneidad y p
´
erdidas usualmente
presentes en la muestra, conducen a distorsiones y artefactos
en las im
´
agenes obtenidas [5].
Existe una gran cantidad de t
´
ecnicas de inversi
´
on para
obtenci
´
on de im
´
agenes en TOA que pueden ser agrupadas
dentro de 4 categor
´
ıas: I) algoritmos en el dominio del
tiempo (retroproyecci
´
on); II) algoritmos en el dominio de
la frecuencia; III) algoritmos de inversi
´
on del tiempo y IV)
aquellos basados en modelos matriciales [5]. Las primeras
dos categor
´
ıas, en gran medida, permiten obtener soluciones
cerradas al problema de inversi
´
on. El enfoque que mejores
resultados ha tenido en los trabajos experimentales es el
algoritmo de retroproyecci
´
on, el cual es muy simple de
implementar y es aceptable para la mayor
´
ıa de las configu-
raciones experimentales para obtenci
´
on de IOA, a
´
un cuando
no es exacta. Por este motivo es que en este trabajo nos
abocaremos a este enfoque.
Los algoritmos en el dominio del tiempo se basan en
proyectar cada una de la se
˜
nales temporales OA unidimen-
sionales en el espacio tridimensional de una forma que
es consistente con el principio de tiempo de vuelo [5].
El proceso de retroproyecci
´
on generalmente involucra tres
pasos. El primer paso es el pre-procesamiento sobre cada
una de la se
˜
nales ac
´
usticas medidas, cuya operaci
´
on depende
de la implementaci
´
on particular del algoritmo. Luego se
realiza la retroproyecci
´
on, la cual consiste en tomar cada una
de las se
˜
nales pre-procesadas y proyectarlas sobre esferas
conc
´
entricas en el espacio 3D. Finalmente, para obtener la
imagen, se realiza la suma sobre todos los datos calculados
en el paso anterior.
En la Fig.1 se puede observar un esquema del proceso
de retroproyecci
´
on 2D. Una se
˜
nal ac
´
ustica originada en la
posici
´
on r es captada por 4 sensores (Fig. 1.a). Las se
˜
nales
medidas poseen un retardo proporcional al m
´
odulo de la
distancia entre la ubicaci
´
on de la fuente y el sensor (Fig.
1.b). Cada una de las se
˜
nales se proyecta en arcos de
circunferencia con un radio igual a la distancia mencionada.
La
´
unica posici
´
on en donde todos los arcos se intersectan y
donde la contribuci
´
on de todas las se
˜
nales retroproyectadas
se suman coherentemente es en la posici
´
on r. Si la f
´
ormula
de reconstrucci
´
on es exacta y el n
´
umero de sensores es
suficientemente grande, se puede recuperar el punto fuente.
Existen muchas implementaciones de este enfoque, siendo
la m
´
as completa el algoritmo de retroproyecci
´
on universal
desarrollado por M. Xu y coautores [6] que entrega una
soluci
´
on exacta para las superficies de detecci
´
on m
´
as co-
munes: esf
´
erica, cil
´
ındrica y plana.
Los algoritmos de reconstrucci
´
on de imagen, general-
mente, suponen sensores ideales. En este trabajo se define
a un sensor ideal como el caso l
´
ımite donde su se
˜
nal de
salida sigue instant
´
aneamente a la presi
´
on ac
´
ustica. Entre los
pol
´
ımeros que se usan para estos sensores, se puede encon-
trar el fluoruro de polivinilideno (PVDF) y sus copol
´
ımeros,
que presentan gran inter
´
es debido a sus propiedades f
´
ısicas
r
r
s1
r
s2
r
s3
r
s4
r-r
s3
r-r
s4
r-r
s1
r-r
s2
(a)
r
r
s1
r
s2
r
s3
r
s4
r
s1
r
s2
r
s3
r
s4
r-r
s1
r-r
s2
(c)
(b)
r-r
s3
r-r
s4
ct
ct
ct
ct
Fig. 1. Ilustraci
´
on del proceso de retroproyecci
´
on 2D. (a) Un punto fuente
es captado por 4 sensores. (b) Se
˜
nales temporales medidas por los sensores.
(c) Cada una de las se
˜
nales se proyecta en arcos de circunferencia con un
radio igual a la distancia entre el punto fuente y el sensor.
que permiten su utilizaci
´
on en muchas aplicaciones. Estos
materiales son flexibles, est
´
an disponibles como pel
´
ıculas
delgadas, tienen un gran ancho de banda ac
´
ustica, y sus
valores de impedancia ac
´
ustica est
´
an pr
´
oximos a los del
agua y los tejidos biol
´
ogicos (a las frecuencias de inter
´
es)
[7]. Las caracter
´
ısticas del PVDF y sus copol
´
ımeros han
sido extensamente estudiadas, destac
´
andose que sus ventajas
son relevantes para transductores de recepci
´
on de ultrasonido
[8].
Sensores del tipo descripto distan de ser ideales. Den-
tro de nuestro conocimiento no hay trabajos que estu-
dien detalladamente las consecuencias del uso de sensores
pol
´
ımericos piezoel
´
ectricos no ideales en los algoritmos de
reconstrucci
´
on de imagen.
En este trabajo nos enfocaremos en los aspectos t
´
ecnicos
del problema inverso ac
´
ustico en TOA, estudiando c
´
omo
las distorsiones introducidas por el sensor afectan a
los algoritmos para reconstrucci
´
on de las im
´
agenes. En
particular, se analiza la influencia del uso de sensores
polim
´
ericos piezoel
´
ectricos de gran apertura sobre la
eficiencia del algoritmo de retroproyecci
´
on universal que
se describe en la siguiente secci
´
on.
II. ALGORITMO DE RETROPROYECCI
´
ON UNIVERSAL
Cuando la duraci
´
on del pulso l
´
aser es mucho m
´
as corta
que los tiempos de relajaci
´
on t
´
ermica y mec
´
anica involu-
crados en el proceso OA, la presi
´
on en una posici
´
on r =
(x, y, z) responde a la ecuaci
´
on de onda [9],
∂
2
∂t
2
− c
2
∇
2
p(r, t) = 0 (1)
con las siguientes condiciones iniciales,
Revista elektron, Vol. 1, No. 2, pp. 58-65 (2017)
ISSN 2525-0159
59
http://elektron.fi.uba.ar
Fig. 2. Esquema para obtenci
´
on de IOA 2-D con sensores lineales.
p(r, 0) = p
0
(r)
∂p
dt
(r, 0) = 0 r R
3
donde c es la velocidad del sonido en el medio y p
0
(r) es
el perfil espacial de presi
´
on inicial generado por la muestra.
La clave del problema inverso es reconstruir la presi
´
on
inicial p
0
(r) partiendo de los datos medidos. La se
˜
nal OA
captada por un sensor puntual en la posici
´
on r
s
se puede
expresar como [4]:
p
d
(r
s
, t) =
∂
∂t
"
t
4π
Z Z
|r
s
−r|=ct
p
0
(r)dΩ
#
(2)
donde dΩ es el elemento de
´
angulo s
´
olido del vector r
respecto de r
s
.
Si se ubican sensores puntuales ideales sobre una super-
ficie S plana o cil
´
ındrica alrededor del objeto irradiado, se
puede utilizar el algoritmo de retroproyecci
´
on universal [6]
para obtener el perfil inicial de presi
´
on:
p
o
(r) = −
2
c
2
Ω
S
Z
S
f
ubp
(r
s
, |r − r
s
|) cos(α)dS (3)
f
ubp
= (t
−1
∂p
∂t
)
donde Ω
S
es 2π o 4π ya sea una superficie plana o
cil
´
ındrica, respectivamente, y α el
´
angulo entre el vector
normal de la superficie de detecci
´
on n
s
y el vector r
s
− r.
f
ubp
es una funci
´
on de pre-procesado utilizada para resaltar
caracter
´
ısticas de la se
˜
nal. Si bien existen diversas formas
funcionales de f
ubp
, todas ellas cumplen con la finalidad
de filtrar las componentes de baja frecuencia y dejar pasar
las altas frecuencias, con un comportamiento lineal en el
medio [6]. Esto es muy necesario en retroproyecci
´
on ya
que el m
´
etodo en s
´
ı mismo realza las componentes de
baja frecuencia. Por lo tanto, con f
ubp
, se acent
´
uan las
caracter
´
ısticas de contraste de la imagen (altas frecuencias),
mientras que se minimiza el efecto de im
´
agenes borrosas
(bajas frecuencias).
A. Algoritmo de retroproyecci
´
on 2-D
Si en lugar de detectores puntuales sobre una superficie
S, se dispusieran sensores lineales a lo largo de una circun-
ferencia C y solidarios con el eje z (ver Fig. 2), la se
˜
nal
de presi
´
on captada por uno de estos sensores ubicado en la
posici
´
on (x
s
, y
s
, z) es:
p
s
(x
s
, y
s
, t) =
+L/2
Z
−L/2
p(x
s
, y
s
, z, t)dz (4)
Bajo la suposici
´
on de que la distribuci
´
on inicial de
presi
´
on tiene soporte acotado, es decir, ∂p(r)/∂z = 0 y
∂
2
p(r)/∂z
2
= 0 para |z| ≥ L/2, se puede establecer que:
Z
+∞
−∞
∂
2
p(r)
∂
2
z
dz =
Z
L/2
−L/2
∂
2
p(r)
∂
2
z
dz =
∂p(r)
∂z
dz
z=L/2
z=−L/2
= 0
(5)
Utilizando la igualdad anterior, se puede probar que
p
s
(x
s
, y
s
, t) satisface la ec. (1) bidimensional:
∂
2
∂t
2
− c
2
∂
2
∂x
2
− c
2
∂
2
∂y
2
p
s
(x
s
, y
s
, t) = 0 (6)
p
s
(x
s
, y
s
, 0) = p
s
0
(x
s
, y
s
) =
Z
+∞
−∞
p
0
(x
s
, y
s
, z)dz
∂p
s
∂t
(x
s
, y
s
, 0) = 0
Por lo tanto, a partir de las se
˜
nales captadas por detectores
lineales dispuestos sobre la curva C, es posible reconstruir
la proyecci
´
on 2-D de la distribuci
´
on de presi
´
on inicial en
el plano xy. Este m
´
etodo est
´
a descripto en el trabajo de P.
Burgholzer, et al. [9]. A continuaci
´
on se detallan los pasos
para adaptar la ec. (3) al esquema planteado en la Fig. 2.
Integrando a lo largo del eje z a ambos lados de la ec.
(3):
p
0
(x, y) = −
2
c
2
Ω
S
Z
∞
−∞
Z
S
f
ubp
(r
s
, |r − r
s
|)cos(α)dS
dz
(7)
siendo ahora cos(α) definido como:
cos(α) = n
s
·
r
s
− r
|r
s
− r|
= cos(γ)
ρ
p
ρ
2
+ (z − z
s
)
2
(8)
donde γ es el
´
angulo entre el vector normal a la curva C,
n
c
, y (x
s
− x, y
s
− y), y ρ = [(x − x
s
)
2
+ (y − y
s
)
2
]
1/2
es
la distancia entre el punto (x, y) y la ubicaci
´
on del sensor
(x
s
, y
s
).
Cambiando el orden de integraci
´
on de la ec. (7) y rees-
cribiendo el diferencial de superficie, p
o
(x, y) es igual a:
−2
Ω
S
c
2
Z
C
Z
∞
−∞
Z
∞
−∞
f
ubp
(r
s
, β)
β
dz dz
s
ρ cos(γ)dC
(9)
β =
p
ρ
2
+ (z − z
s
)
2
Teniendo en cuenta la simetr
´
ıa de la disposici
´
on de los
sensores respecto del eje z y sustituyendo z − z
s
= ¯z, la
Revista elektron, Vol. 1, No. 2, pp. 58-65 (2017)
ISSN 2525-0159
60
http://elektron.fi.uba.ar
integral interna de la ec. (9) puede evaluarse a trav
´
es de la
siguiente expresi
´
on:
2
Z
∞
0
f
ubp
(x
s
, y
s
,
p
ρ
2
+ ¯z
2
)d¯z (10)
Si la velocidad del sonido c en el medio es constante y
homog
´
enea, se puede establecer una relaci
´
on sencilla entre
el m
´
odulo de la distancia y el tiempo. De esta forma, se
realiza un cambio de variables integrando en t en lugar de
¯z obteniendo:
p
0
(x, y) =
−4
c
2
Ω
S
Z
C
Z
∞
ρ
f
ubp
(x
s
, y
s
, t)
p
t
2
− ρ
2
/c
2
dt
!
w(ρ, γ)dC
(11)
donde w = ρ cos(γ) y se denomina factor de peso.
Cabe destacar que el m
´
etodo detallado entrega una
soluci
´
on exacta si la superficie de detecci
´
on S = CxR es
un cilindro o un plano infinito. Por lo tanto, la contraparte
bidimensional ser
´
a exacta si C es una circunferencia o una
l
´
ınea infinita.
B. Discretizaci
´
on del algoritmo
La ec. (11) requiere un conjunto continuo de datos. Dado
que en una aplicaci
´
on pr
´
actica esto no es posible, resulta
necesario discretizarla.
En primer lugar se parametriza la curva C como
(x
s
, y
s
) = (x
s
(ξ), y
s
(ξ)), utilizando un par
´
ametro ξ que
pertenece a un intervalo I. Entonces la ec. (11) puede
reescribirse de la siguiente forma:
4
c
2
Ω
S
Z
I
Z
∞
ρ/c
f
ubp
(x
s
, y
s
, t)
p
t
2
− ρ
2
/c
2
dt
!
w(ρ, γ) h
C
(12)
h
C
=
s
dx
s
dξ
2
+
dy
s
dξ
2
dξ
Luego, se realiza un muestreo uniforme temporal y espa-
cial del eje de tiempo y de la curva C, respectivamente,
t
m
= m∆t m ∈ [0, . . . , N
t
− 1]
ξ
k
= k∆ξ k ∈ [0, . . . , N
ξ
− 1]
donde ∆t y ∆ξ son los pasos de tiempo y de longitud (curva
C), respectivamente. N
t
representa la cantidad de muestras
del eje de tiempo y N
ξ
la cantidad de sensores empleados.
Por lo tanto, la presi
´
on del sensor k para el tiempo discreto
t
m
es: p
m,k
= p(x
s
(ξ
k
), y
s
(ξ
k
), t
m
).
El siguiente paso es discretizar la integral interna de la
ec. (12). La funci
´
on f
ubp
es aproximada usando el m
´
etodo
de diferencias finitas de segundo orden:
q
m,k
=
p
m+1,k
t
m+1
−
p
m−1,k
t
m−1
/(2c
2
∆t) (13)
con q
0,k
= 0.
Para evaluar la funci
´
on en puntos no muestreados, se
realiza una interpolaci
´
on lineal en el intervalo,
L
k
[q](t) = q
m,k
+
t − t
m
∆t
(q
m+1,k
− q
m,k
) (14)
donde t ∈ [t
m
, t
m+1
]
Luego, la integral interna de la ec. (12) es evaluada para
ρ/c = t
m
y la funci
´
on de pre-procesado es reemplazada por
su versi
´
on discreta e interpolada:
ˆq
m,k
:=
Z
∞
t
m
L
k
[q]
p
t
2
− (t
m
)
2
dt =
N
t
−1
X
m
0
=m
q
m
0
,k
Z
t
m
0
+1
t
m
0
dt
p
t
2
− (t
m
)
2
!
+
+
N
t
−1
X
m
0
=m
q
m
0
+1,k
− q
m
0
,k
∆t
Z
t
m
0
+1
t
m
0
(t − t
m
0
)dt
p
t
2
− (t
m
)
2
!
(15)
Para que la implementaci
´
on num
´
erica sea eficiente y
precisa, es crucial que las integrales de la ec. (15) sean
evaluadas de forma anal
´
ıtica:
ˆq
m,k
=
N
t
−1
X
m
0
=m
a
m
m
0
q
m
0
,k
+
1
h
r
N
t
−1
X
m
0
=m
b
m
m
0
(q
m
0
+1,k
− q
m
0
,k
)
(16)
a
m
m
0
:= [log(t +
p
t
2
− (t
m
)
2
)]
t
m
0
+1
t=t
m
0
b
m
m
0
:= t
m
0
a
m
m
0
+ [
p
t
2
− (t
m
)
2
]
t
m
0
+1
t=t
m
0
Para evaluar la funci
´
on ˆq en puntos no muestreados, se
realiza tambi
´
en una interpolaci
´
on lineal como el de la ec.
(14).
El
´
ultimo paso consiste en evaluar la integral de l
´
ınea
correspondiente a la curva C de forma num
´
erica:
¯p
0
(x
i
, y
i
)
∼
=
p
i
0
= −
4
c
2
Ω
S
N
ξ
X
k=0
L
k
[ˆq](ρ
k,i
)w
k,i
h
k
C
(17)
donde (x
i
, y
i
), i ∈ [0, . . . , N ]
2
, representan los puntos
pertenecientes a una red uniforme del
´
area que encierra la
curva C y ρ
k,i
= [(x
i
−x
k
s
)
2
+(y
i
−y
k
s
)
2
]
1/2
es la distancia
entre el punto (x
i
, y
i
) y la posici
´
on del sensor (x
k
s
, y
k
s
).
Debido a que las diferencias finitas, la interpolaci
´
on lineal
y el m
´
etodo trapezoidal son aproximaciones de segundo
orden, una cota error del algoritmo se puede expresar como
[9]:
|p
i
0
− ¯p
0
(x
i
, y
i
)| ≤ κ · max{∆t
2
, ∆ξ
2
} + (N
t
∆t) (18)
donde κ es un constante independiente de i, ∆t y ∆ξ
(N
t
∆t) representa el error de truncamiento de la integral
interna de la eq. (11)
Revista elektron, Vol. 1, No. 2, pp. 58-65 (2017)
ISSN 2525-0159
61
http://elektron.fi.uba.ar
PVDF
Substrato
Ondas
Acústicas
Ondas
Acústicas
Fig. 3. Esquema (no a escala) del sensor. Vista lateral (izq.) y superior
(der.)
III. MODELO DEL SENSOR POLIM
´
ERICO
PIEZOEL
´
ECTRICO
Para modelar este tipo de transductores se usa el modelo
par
´
ametrico descrito en [10]. Este modelo se caracteriza por:
i) suponer que cada punto de la superficie del sensor se com-
porta como un detector puntual que puede ser aproximado
como un sistema lineal e invariante en el tiempo (LTI); ii)
los sistemas LTI que componen la superficie son id
´
enticos;
iii) las principales caracter
´
ısticas del sensor se modelan a
trav
´
es de funciones transferencia o filtros aplicadas sobre
cada uno de los sistemas LTI.
Nuestro estudio se enfoca en sensores de pel
´
ıcula delgada
de PVDF de gran apertura y cuya
´
area de detecci
´
on puede
aproximarse a una l
´
ınea (cuyo largo L es mucho mayor que
su ancho A). Se denomina cara anterior del film a aquella
que se encuentra del lado de donde provienen las ondas
ac
´
usticas generadas en la muestra. Consecuentemente, la
cara opuesta, es la cara posterior (ver Fig. 3).
Es necesario agregar que, dentro de nuestro conocimiento,
no hay modelos comerciales de sensores piezoel
´
ectricos
con geometr
´
ıa lineal para obtenci
´
on de IOA. Sin embargo,
existen numerosos trabajos de investigaci
´
on donde han sido
aplicados, incluso en sistemas de detecci
´
on para TOA [11].
En este trabajo consideramos dos aspectos f
´
ısicos del
sensor: las propiedades del material polim
´
erico que lo com-
ponen y aquellos relativos con su forma y la relaci
´
on con
su entorno.
Una presi
´
on p
s
(r
s
, t) incidente en la cara anterior del
sensor induce un desplazamiento el
´
ectrico D(r
s
, t) en el
mismo puede ser calculado a trav
´
es de la siguiente ecuaci
´
on
[10]:
D(r
s
, ω) = P(r
s
, ω) · d
33
(ω) (19)
donde d
33
es el coeficiente piezoel
´
ectrico y D(r
s
, ω) y
P (r
s
, ω) las transformadas de Fourier del desplazamiento
el
´
ectrico y la presi
´
on, respectivamente. En el caso del PVDF,
y suponiendo que la onda ac
´
ustica solo viaja en la direcci
´
on
perpendicular a su superficie (material de reacci
´
on normal)
[12], d
33
puede ser descripto por una funci
´
on Havriliak-
Negami [10]:
d
33
(ω) = d
∞
33
+
∆d
33
(1 + (iτ
0
ω)
α
)
β
(20)
donde d
∞
33
es el valor l
´
ımite del coeficiente piezoel
´
ectrico
a altas frecuencias, ∆d
33
la intensidad del proceso de
relajaci
´
on, τ
0
el tiempo caracter
´
ıstico de la relajaci
´
on y α
y β son los par
´
ametros de forma que describen la asimetr
´
ıa
y ensanchamiento de la relajaci
´
on, respectivamente. Estos
par
´
ametros de forma son ambos n
´
umeros positivos y su
producto es siempre menor que 1 [13].
Suponiendo que la velocidad de propagaci
´
on del sonido
en el PVDF, v
p
, no depende de la frecuencia, la densidad
de carga el
´
ectrica, σ
p
, puede ser determinada a partir de la
siguiente expresi
´
on [10]:
σ
p
(r
s
, ω) = D(r
s
, ω) · W (ω, Γ
p
, Γ
a
, δ, v
p
, a
p
) (21)
La funci
´
on W incorpora las reflecciones de las ondas
ac
´
usticas en las caras posterior y anterior del sensor as
´
ı
como tambi
´
en el espesor δ y la atenuaci
´
on ac
´
ustica del film
piezoel
´
ectrico:
W (ω, Γ
p
, Γ
a
, δ, v
p
, a
p
) =
F
v
p
δ
Y
v
p
t
δ
+ Γ
p
Y
v
p
t − δ
δ
+ (22)
+Γ
p
Γ
a
Y
v
p
t − 2δ
δ
e
−a
p
v
p
t
donde F es el operador transformada de Fourier, a
p
el factor
de atenuaci
´
on de la onda ac
´
ustica al propagarse a trav
´
es del
film, Π(
z
δ/v
p
) una ventana rectangular entre 0 y δ/v
p
, y
Γ
a
y Γ
p
los coeficientes de reflexi
´
on de la cara anterior
y posterior, respectivamente. Esta ecuaci
´
on supone que Γ
a
,
Γ
p
y a
p
son independientes de la frecuencia. Cabe destacar
que esta hip
´
otesis fue corroborada experimentalmente en el
trabajo [10]. Las sucesivas reflexiones de la onda ac
´
ustica
entre las caras anterior y posterior del film son reemplazadas
por el equivalente de propagaci
´
on de la onda a trav
´
es de
varias capas. En esta ecuaci
´
on son consideradas solo dos
reflexiones ya que se supone que el material polim
´
erico
presenta grandes p
´
erdidas.
Finalmente, la carga el
´
ectrica q(t) en los electrodos del
film de PVDF se determina a trav
´
es de la siguiente ecuaci
´
on:
q(t) = A
Z
L/2
−L/2
F
−1
{σ
p
(r
s
, w)} dz (23)
donde F
−1
es la transformada de Fourier inversa. La se
˜
nal
que entrega el sensor al sistema de detecci
´
on es proporcional
a q(t). Es necesario agregar que, este tipo de sensores suele
emplearse en conjunto con un amplificador de transimpedan-
cia. En el caso que el ancho de banda del mismo sea mucho
mayor que el del sensor, q(t) se puede obtener simplemente
integrando la se
˜
nal medida.
IV. SIMULACI
´
ON Y RESULTADOS
En este trabajo se considera que el tiempo caracter
´
ıstico
del pulso l
´
aser que irradia a la muestra es mucho menor que
cualquiera de los otros tiempos involucrados en el proceso
OA. En este caso el objeto dentro de la muestra que absorbe
la luz irradiada por el l
´
aser es una esfera sumergida en agua
a una temperatura T
a
= 298K. La esfera tiene un radio a
y su centro se ubica en la posici
´
on r
a
. La presi
´
on ac
´
ustica
Revista elektron, Vol. 1, No. 2, pp. 58-65 (2017)
ISSN 2525-0159
62
http://elektron.fi.uba.ar
Par
´
ametros Fijos Par
´
ametros Variables
valor de referencia valor inicial
δ 25 µm Γ
a
-0.44
τ
0
140 ns Γ
p
-0.99,-0.10
∆d
33
34.5 pC/N a
p
2·10
4
m
−1
d
∞
33
1.1 pC/N a 10-100 µm
v
p
2200 m/s R
C
0.1-1 mm
α 0.5 ∆z 1-10 um
β 0.68 ∆t 1ns
c 1500 m/s N
ξ
128
L 20 mm N
t
7000
TABLA I
PAR
´
AMETROS FIJOS Y VARIABLES USADOS EN LA SIMULACI
´
ON. TODOS
LOS VALORES SON PARA LA TEMPERATURA DE TRABAJO T
a
= 298 K.
generada por esta esfera en la posici
´
on r en funci
´
on del
tiempo se calcula a trav
´
es de la siguiente expresi
´
on:
p(r, t) = Υ
a
2
− (|r − r
a
| − ct)
2
2 |r − r
a
|
H (a − ||r − r
a
| − ct|)
(24)
donde H es la funci
´
on de Heaviside y Υ una constante de
proporcionalidad cuyo valor depende de las propiedades del
l
´
ıquido que rodea a la esfera y de la energ
´
ıa l
´
aser absorbida
por la misma en el proceso OA.
En todas las simulaciones, con el objetivo de minimizar
el efecto de discretizaci
´
on de la curva C, se us
´
o un valor de
N
ξ
≥ 128 . Para calcular la presi
´
on captada por el sensor,
se discretiza su largo de forma uniforme con un muestreo
espacial ∆z = 0.1 · a. Por otra parte, la curva C, donde se
ubican los detectores, es una circunferencia de radio R
C
=
10·a. Estos dos par
´
ametros se eligen en funci
´
on del tama
˜
no
de la fuente esf
´
erica para facilitar las comparaciones entre
simulaciones. El valor del par
´
ametro N
t
depende del ∆t
elegido y del tiempo final de integraci
´
on.
La respuesta de un sensor ideal es:
q(t) = A d
33
(0)
L/∆z
X
k=1
p(r
k
, t) (25)
Para obtener la respuesta de un sensor no ideal se utiliz
´
o
el modelo param
´
etrico detallado en la secci
´
on anterior. Se
consideraron los siguientes aspectos f
´
ısicos del sensor : i)
su espesor; ii) reflexiones en la cara posterior y anterior (de-
sadaptaci
´
on de impedancias ac
´
usticas); iii) propiedades del
material polim
´
erico (p
´
erdidas electromec
´
anicas y relajaci
´
on
piezoel
´
ectrica). Para mostrar el efecto de la desadaptaci
´
on
de impedancias en la cara posterior sobre la reconstrucci
´
on
de la imagen, se hicieron un conjunto de simulaciones para
dos casos de inter
´
es experimental: Γ
p
= −0.99 (aire) y
Γ
p
= −0.10 (acr
´
ılico).
Los valores de los distintos par
´
ametros utilizados en las
simulaciones se muestran en la tabla I. Los valores de
referencia relacionados con el PVDF fueron obtenidos de
[14].
La integral de la ec. 22 se calcula num
´
ericamente usando
la funci
´
on quadv disponible en GNU/octave.
10
6
10
7
10
8
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
Amplitud normalizada
Frecuencia [Hz]
Espesor
Espesor + reflexión
Espesor + reflexión
+ material
Fig. 4. M
´
odulo normalizado de las transferencias, con Γ
p
= −0.99 .
10
6
10
7
10
8
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Frecuencia [Hz]
Amplitud normalizada
Espesor
Espesor + reflexión
Espesor + reflexión
+ material
Fig. 5. M
´
odulo normalizado de las transferencias, con Γ
p
= −0.1 .
En las Figs. 4 y 5 se observan las funciones transferencias
de los aspectos f
´
ısicos del sensor para los dos valores de
Γ
p
. En l
´
ınea continua se muestra el efecto de un espesor
finito. La curva con trazo a rayas representa tanto el espesor
como las reflexiones en las caras posterior y anterior del
sensor. Finalmente, la curva de puntos es la transferencia que
incluye tambi
´
en las propiedades del material piezoel
´
ectrico.
El espesor del sensor funciona como un filtro pasabajos.
Por otro lado, las reflexiones producen grandes cambios
en el espectro de la se
˜
nal respecto del de un sensor ideal.
Por consiguiente, es esperable tener im
´
agenes reconstruidas
con artefactos, en especial cuando Γ
p
es elevado (ver Fig.
4). En la Fig 5 se puede apreciar el efecto de adaptar las
impedancias ac
´
usticas entre la cara posterior de la pel
´
ıcula
delgada de PVDF y el substrato donde est
´
a apoyada (menor
valor de Γ
p
). Se observa una reducci
´
on importante del
efecto distorsivo de las reflexiones con respecto al caso
desadaptado (Fig. 4). Finalmente, la funci
´
on transferencia
que introduce las propiedades del material, limita la re-
spuesta a frecuencias elevadas (≥ 1MHz), imponiendo una
frecuencia de corte a
´
un menor que en el caso de la funci
´
on
que representa el espesor del film.
Revista elektron, Vol. 1, No. 2, pp. 58-65 (2017)
ISSN 2525-0159
63
http://elektron.fi.uba.ar
a = 25μm
Γ
p
= -0.1
Espesor-reflexión-material
Γ
p
= -0.99
Ideal
x [m]
y [m]
-1 0 1
x 10
-4
-1
0
1
x 10
-4
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
x [m]
y [m]
-1 0 1
x 10
-4
-1
0
1
x 10
-4
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
x [m]
y [m]
-1 0 1
x 10
-4
-1
0
1
x 10
-4
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
Espesor-reflexión-material
Fig. 6. Reconstrucci
´
on 2-D de la distribuci
´
on inicial de presi
´
on ini-
cial ¯p
o
(r), considerando las propiedades f
´
ısica de sensores polim
´
ericos
piezoel
´
ectricos de banda ancha con geometr
´
ıa lineal.
En la Fig. 6 se muestra un ejemplo de las reconstrucciones
de la distribuci
´
on inicial de presi
´
on ¯p
o
(r), en valores nor-
malizados, para a = 25µm. La primera imagen de la figura
muestra el caso ideal. Las otras reconstrucciones presentan
los efectos producidos al considerar el espesor, el material
piezoel
´
ectrico y las reflexiones con Γ
p
= −0.99 (imagen
central) y Γ
p
= −0.1 (imagen inferior). De la figura se
puede apreciar que la introducci
´
on de los aspectos f
´
ısicos del
sensor afectan notoriamente la reconstrucci
´
on de la imagen.
Para facilitar el an
´
alisis, en las Figs. 7, 8 y 9 se muestran
las reconstrucciones 2-D de ¯p
o
(x, y = 0) producidas por
esferas con radios 100 µm, 25 µm y 10 µm, respectivamente.
-3 -2 -1 0 1 2 3
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
x[μm]
Presión normalizada
Ideal
Espesor
Espesor + reflexión
p
=-0.99
Espesor + reflexión
p
=-0.1
Espesor + reflexión
+ material
Γ
p
=-0.99
Espesor + reflexión
+ material
Γ
p
=-0.1
Γ
Γ
0
0
0
0
0
0
0
0
0
0
0
0
Fig. 7. Reconstrucci
´
on 2-D de ¯p
o
(x, y = 0) para una fuente de presi
´
on
esf
´
erica de radio a = 100µm.
-8 -6 -4 -2 0 2 4 6 8
0
0.2
0.4
0.6
0.8
1
000
0
0 0
0
0
x[μm]
Presion normalizada
Ideal
Espesor
Espesor + reflexión
Γ
p
=-0.99
Espesor + reflexión
Γ
p
=-0.1
Espesor + reflexión
+ material Γ
p
=-0.99
Espesor + reflexión
+ material Γ
p
=-0.1
Fig. 8. Reconstrucci
´
on 2-D de ¯p
o
(x, y = 0) para una fuente de presi
´
on
esf
´
erica de radio a = 25µm.
-6 -4 -2 0 2 4 6
0
0.2
0.4
0.6
0.8
1
00
0
0 0
0
Presiónnormalizada
Ideal
Espesor
Espesor + reflexión
p
=-0.99
Espesor + reflexión
p
=-0.1
Espesor + reflexión
+ material
p
=-0.99
Espesor + reflexión
+ material
p
=-0.1
Γ
Γ
Γ
Γ
x[μm]
Fig. 9. Reconstrucci
´
on 2-D de ¯p
o
(x, y = 0) para una fuente de presi
´
on
esf
´
erica de radio a = 10µm.
Revista elektron, Vol. 1, No. 2, pp. 58-65 (2017)
ISSN 2525-0159
64
http://elektron.fi.uba.ar
Para el caso donde se tiene en cuenta un espesor finito
(curvas en color rojo), se observa un efecto de filtrado de las
componentes de alta frecuencia (pasabajos). En concordan-
cia con esto, a medida que se reduce el tama
˜
no de la fuente
de presi
´
on, la atenuaci
´
on y el radio de la esfera reconstruida
aumentan en relaci
´
on al caso ideal.
Cuando se consideran tambi
´
en las reflexiones con Γ
p
=
−0.99 (curvas con trazo a rayas color negro) se pro-
ducen artefactos en la reconstrucci
´
on como por ejemplo la
aparici
´
on de presiones negativas. Esto es debido a que la
transferencia de este aspecto f
´
ısico (ver Fig. 4) posee una
importante atenuaci
´
on a bajas frecuencias (<1MHz) y una
fuerte resonancia ac
´
ustica en el rango entre 10 MHz y 100
MHz. Es importante destacar que para el caso particular de
a = 25µm (Fig. 8), la mayor parte del contenido espectral
del pulso de presi
´
on coincide con la susodicha resonancia,
generando una amplificaci
´
on en la presi
´
on.
Las curvas con trazo a rayas color verde son las mismas
simulaciones que para el caso anterior, pero con una mejor
adaptaci
´
on de impedancias, es decir Γ
p
= −0.1. Se puede
ver una reducci
´
on de los efectos adversos de la reflexiones
en las distintas interfases entre los materiales que componen
al detector.
Las curvas de puntos de colores violeta y marr
´
on
representan las reconstrucciones obtenidas al agregar las
propiedades del material. La introducci
´
on de un factor de
atenuaci
´
on distinto de cero y de la relajaci
´
on piezoel
´
ectrica
reducen las distorsiones provocadas por las reflexiones as
´
ı
como tambi
´
en la amplitud de la se
˜
nal entregada por el
sensor. Adem
´
as se puede observar que, a pesar de reducir las
reflexiones, el radio de las esferas reconstruidas no mejora,
siendo mayores que para el caso considerando solo el
espesor. Esto denota que la resoluci
´
on espacial est
´
a limitada
principalmente por las propiedades del material.
Otra cosa que puede apreciarse observando las Figs. 7, 8
y 9 es que la forma de las curvas tambi
´
en es dependiente del
radio de la esfera. Este fen
´
omeno se debe a que a menores
valores de a, mayor es el ancho de banda de la se
˜
nal de
presi
´
on y, por lo tanto, para radios peque
˜
nos, se hace m
´
as
notorio el efecto de atenuaci
´
on de las componentes de alta
frecuencia de la imagen a ser reconstruida.
V. CONCLUSIONES
En este trabajo se estudi
´
o c
´
omo los aspectos f
´
ısicos de
un sensor no ideal afectan la reconstrucci
´
on de im
´
agenes
OA. En particular se analiz
´
o el caso de sensores polim
´
ericos
piezoel
´
ectricos de banda ancha con geometr
´
ıa lineal y su
efecto sobre el algoritmo de retroproyecci
´
on universal.
El estudio se centr
´
o sobre tres aspectos f
´
ısicos del sensor:
espesor finito, reflexiones por desadaptaci
´
on de impedancia
ac
´
usticas y las propiedades del material polim
´
erico usado
para captar las ondas de presi
´
on (pel
´
ıcula delgada de PVDF).
Al tener en cuenta el espesor se observ
´
o una disminuci
´
on
de la resoluci
´
on espacial del sistema. El an
´
alisis sobre las
reflexiones mostr
´
o que la elecci
´
on del sustrato es esencial
para lograr im
´
agenes con baja distorsi
´
on. Es importante
hacer notar que cuando se agregan las propiedades del
material, se observa una reducci
´
on de los efectos de las
reflexiones ac
´
usticas en la cara posterior. Esto pone en
evidencia que no es necesario lograr una perfecta adaptaci
´
on
de impedancias ac
´
usticas para obtener una perfomance sat-
isfactoria, verificando lo mencionado en [10].
En resumen, los resultados obtenidos en este trabajo
indican que es necesario caracterizar a priori el sensor que
ser
´
a usado para captar las ondas ultras
´
onicas en un sistema
OA. La informaci
´
on es relevante para determinar si el mismo
es o no apto para la aplicaci
´
on particular. Actualmente
nos encontramos trabajando en el desarrollo de un modelo
inverso que permita minimizar los efectos adversos intro-
ducidos por los aspectos intr
´
ınsecos del detector.
AGRADECIMIENTOS
Este trabajo fue apoyado por los subsidios de la Universi-
dad de Buenos Aires (UBACyT 20020160100052BA) y de
la ANPCyT (PICT 2016-2204).
REFERENCIAS
[1] V. Ntziachristos and D. Razansky. “Molecular imaging by means
of multispectral optoacoustic tomography (MSOT)”. Chem. Rev.,
110:2783–2794, 2010.
[2] V. Ntziachristos, J. Ripoll, L. Wang, and R. Weissleder. “Looking
and listening to light: the evolution of whole-body photonic imaging”.
Nat. Biotechnol., 23:313–320, 2005.
[3] C. Lutzweiler and D. Razansky. “Optoacoustic imaging and tomogra-
phy: reconstruction approaches and outstanding challenges in image
performance and quantification”. Sensors, 13:7345–7384, 2013.
[4] M. Xu and L. V. Wang. “Photoacoustic imaging in biomedicine”.
Rev. Sci. Instrum., 77:041101–1–22, 2006.
[5] A. Rosenthal, V. Ntziachristos, and D. Razansky. “Acoustic inversion
in optoacoustic tomograph”. Nat. Biotechnol., 9:318–336, 2013.
[6] M. Xu and L. Wang. “Universal back-projection algorithm for
photoacoustic computed tomography”. Phys. Rev. E, 71:016706–1–
13, 2005.
[7] L. F. Brown. “Design considerations for piezoelectric polymer
ultrasound transducers”. IEEE Trans. Ultrason. Ferroelect. Freq.
Contr., 47:1377, 2000.
[8] Q. X. Chen and P. A. Payue. “Industrial applications of piezoelectric
polymer transducers”. Meas. Sci. Technol., 6:249, 1995.
[9] P. Burgholzer, J. Bauer-Marchallinger, H Gruen, M. Haltmeier, and
G. Paltauf. “Temporal back-projection algorithms for photoacoustic
tomography with integrating line detectors”. Inverse Problems,
23:S65–S80, 2007.
[10] A. Fern
´
andez Vidal, L. Ciocci Brazzano, C. Matteo, P. Sorichetti, and
M. G. Gonz
´
alez. “Parametric modeling of wideband piezoelectric
polymer sensors: design for optoacoustic applications”. Rev. Sci.
Instrum., 88, 2017.
[11] G. Paltauf, P Hartmair, G Kovachev, and R. Nuster. “Piezoelectric
line detector array for photoacoustic tomography”. Photoacoustics,
8:28–36, 2017.
[12] M. Haltmeier and G. Zangerl. “Spatial resolution in photoacoustic
tomography: effects of detector size and detector bandwidth”. Inverse
Problems, 26:125002, 2010.
[13] S. Havriliak and S. Negami. “A complex plane analysis of dispersions
in some polymer systems”. J.Polym.Sci., Polym. Sym., 14:99, 1966.
[14] M. G. Gonz
´
alez, P. A. Sorichetti, L. Ciocci Brazzano, and G. D. San-
tiago. “Electromechanical characterization of piezoelectric polymer
thin films in a broad frequency range”. Polym. Test., 37:163, 2014.
Revista elektron, Vol. 1, No. 2, pp. 58-65 (2017)
ISSN 2525-0159
65
http://elektron.fi.uba.ar
