El problema lineal inverso en PDS y la descomposición en valores singulares

5/9/2018 El problema lineal inverso en PDS y la descomposici n en valores singulares - s...

http://slidepdf.com/reader/full/el-problema-lineal-inverso-en-pds-y-la-descomposicion-en-valo

El problema lineal inverso en PDS y la descomposicion

en valores singulares

20 de abril de 2009

MSc. Miguel Santana Justiz, Dr. Roberto Henry Herrera Marreroemail: (msantana, henry)@ucf.edu.cu

Resumen

En este artıculo se presenta una alternativa para la solucion del problema inverso inher-ente a la supresion de ruido (denoising) en senales directamente en el dominio temporal.El metodo consiste, esencialmente en la transformacion del kernel en una matriz Toeplitz

y su descomposicion en valores singulares SVD . Los resultados experimentales muestranbuenos indicadores de calidad de la senal en presencia de altos niveles de ruido blancogaussiano en sistemas LTI.Se realizan experimentos sobre senales simuladas con un kernel definido positivo.Se realiza una comparacion con una variante del filtro de Tikhonov en terminos de laSVD tambien.

Introduccion

El proceso de degradacion de una senal se puede representar como un operador (o un sistema)h(t) que junto a un termino aditivo de ruido η(t) actua sobre una senal de entrada x(t) paraproducir una senal degradada y(t) como se muestra en la Fig.1. La deconvolucion puedeentenderse como el proceso de obtener de forma aproximada x(t), a partir de la observaciony(t) y del conocimiento de la degradacion en la forma del operador h(t). Se supone que elconocimiento del ruido η(t) se limita a una informacion de naturaleza estadıstica.

Figura 1: Esquema del proceso de degradacion de una senal.

1



El modelo matematico que representa este fenomeno fısico es,

y(t) = x(t) ∗ h(t) + η(t). (1)

El interes de un problema inverso es conocer o estimar la senal real x(t), de forma que elresiduo η(t), sea mınimo, donde x(t) ∗ h(t), representa el producto de convolucion discretode cada punto xt, con el pulso h(t).Dada la senal observada y(t) y una estimacion para el pulso o PSF h(t), se desea extraer lafuncion de reflexividad x(t) . La solucion del problema inverso que se enuncia previamenteen (1) conduce a que la relacion entre estas senales se encuentra minimizando la funcionconvexa, (Vogel, 2002).

η2 =1

2mın y(t) − x(t) ∗ h(t)2. (2)

Existen diversos enfoques para la solucion del problema inverso que encierra el fenomenode la deconvolucion de senales unidimensionales o de la restauracion en el caso de senales

bidimensionales. Estos se refieren a distintos criterios de optimizacion que intentan solucionarel problema directamente sobre la senal en el dominio temporal y otros que lo hacen sobrela senal en un dominio transformado.

1. Desarrollo

En (Vogel, 2002) se propone transformar el kernel h(t) 1D en una matriz de convolucion, H ,a partir de una matriz Toeplitz.Una matriz Toeplitz es una matriz cuadrada, simetrica respecto a su diagonal principal deforma que se cumpla que:

∀hi,j ∈ H ; hi,j = hi+1,j+1 (3)

de donde resulta una matriz fuertemente simetrica.Los sistemas formados por matrices Toeplitz aparecen en una gran variedad de aplicacionesen matematica e ingenierıa. En el procesamiento de senales, la solucion de un sistema Toeplitzse requiere con vista a obtener los coeficientes de los filtros en el dise no de los filtros digitalesrecursivos y en esta aplicacion.

1.1. Descomposicion SVD

El problema matematico que encierra la descomposicion en valores singulares se ha tratado

por varios autores (Highan, 1996; Nash, 1990).De forma general, cualquier matriz A∈ R

m×n tiene descomposicion en valores singulares(SVD)1 (Nash, 1990), tomando cualesquiera valores positivos arbitrarios tales que: B = U Σdonde U T U = I , entonces para la matriz Am×n existe una descomposicion SVD, tal que,A = U ΣV T , donde Σ es una matriz diagonal que contiene los valores singulares tales que:σ1 ≥ σ2 ≥ . . . ≥ σk ≥ 0 , para k = 1, n y U y V son los vectores singulares izquierdos y

1SVD: del ingles, singular value decomposition

2



derechos respectivamente. Expresando la solucion del Sistema Lineal mediante la SVD,

xest =n

j=1

uT j y

σ jv j ∀σ j > 0. (4)

Este sistema es bien-definido (well-posed ) si se satisface la condicion de Picard (Hansen,2002):

||xest||2 =n

j=1

|uT j y|2σ2 j

< ∞.

Esta condicion fallara si σ j → 0 porque no se puede garantizar una dependencia continuaentre la estimacion y la observacion, esto se conoce como problema mal-definido (ill-posed ).

De lo anterior se deduce que sera mal-condicionado el sistema cuyo operador H , tengavalores singulares que decaen gradualmente a cero y la razon de su valores extremos seagrande.

cond(H ) = max{σ j}min{σ j} ; σ j > 0. (5)

Este representa el numero de condicion de la matriz, valor por el cual se multiplicara cualquieralteracion en la entrada del sistema. La solucion en (4) para un ruido aditivo η presente enel sistema quedara (Lavarello, Kamalabadi, & O’Brien, 2006):

xest =n

j=1

uT j y

σ jv j +

n j=1

uT j η

σ jv j ∀σ j > 0. (6)

Si H = I , se estara en presencia de un problema de supresion de ruido (denoising ) y si η = 0

sera un problema de deconvolucion (deblurring ) 2.

1.2. Metodo de Regularizacion

La regularizacion consiste basicamente en introducir alguna clase de informacion a priori

para estabilizar el problema en el sentido de Hadamard. En (Acar & Vogel, 1994) se muestranresultados en senales electricas a partir de la regularizacion de Tikhonov y de la VariacionTotal. (Stefan, Garnero, & Renaut, 2005) expone el uso de las mismas en la reconstrucci on desenales en estudios sısmicos, otros autores como (Herrera, Orozco, Moreno, & Calas, 2005b)utilizan la regularizacion a traves de una transformada wavelet con un filtro de Wiener y

(Hansen, 2002) utiliza la regularizacion con matrices Toeplitz.

1.2.1. Truncamiento de los valores singulares, TSVD

La manera mas natural de regularizar un problema inverso, serıa eliminar los valores singu-lares menos estables, cuando el nivel de ruido en los datos se incrementa. Este metodo esconocido como Truncamiento de la Descomposicion en Valores Singulares (TSVD, Truncated

SVD ) (Lavarello et al., 2006). Para ello se necesita un parametro de regularizacion α, que

2El termino deblurring es equivalente a desenfoque

3



actue como nivel de umbral para los valores singulares. De esta forma la expresi on siguiente(Vogel, 2002):

wα(σ2) =

1, σ2

j > α,

0, σ2 j ≤ α,

(7)

define un filtro wα, que actua sobre los valores singulares, dejando pasar solo aquellos mayoresque α, como se muestra en la Fig. 2, donde desplazamiento a la izquierda del filtro por ladisminucion del parametro de regularizacion α, implica que en la solucion del problema seincluyen los menores valores singulares. Entonces la estabilidad del sistema dependerıa engran medida de la seleccion apropiada del parametro de regularizacion y aun en ese caso setiene la indeterminacion entre supresion del ruido y perdida de informacion util. Si incluimos

10−5

10−4

10−3

10−2

10−1

100

101

0

0.2

0.4

0.6

0.8

1

1.2Funciones del Filtro para TSVD

s2

w α

( s 2 )

α=10−1

α=10−2

α=10−3

Figura 2: Filtro wα(σ2), para la regularizacion TSVD.

la funcion del filtro wα(σ2) en (4) se tendrıa:

xest =n

j=1

wα(σ2 j )

uT j y

σ jv j ∀σ j > 0. (8)

1.3. Situacion experimental, Materiales y Metodos

En todos los experimentos, se usa la misma arquitectura. No se emplea paralelizacion de losmetodos, por cuanto los resultados obtenidos pudieran ser mejorados en trabajos futuros. Seusan cuatro terminos que dan una medida comparativa del comportamiento de los algoritmos:la relacion Senal a Ruido Desenfocada, BSNR (Raimondo & Stewart, 2007), la relacionSenal a Ruido de la senal estimada, SNR (Ma, Wang, & Du, 2006), la mejora en la relacionsenal-ruido, ISNR(Neelamani, Choi, & Baraniuk, 2004) y el error cuadratico medio, MSE

(Donoho, 1995; Ma et al., 2006).

4



1.3.1. Deconvolucion de senales simuladas, caso 1D

En esta seccion se basa el analisis en los resultados previos de (Vogel, 2002). El codigo enMatLab esta disponible en su sitio web3, que es referencia en investigaciones actuales sobre eltema (Cui, Lamm, & Scofield, 2007). Para los experimentos, se considera la ecuacion integralde Fredholm de primer tipo (Vogel, 2002):

g(x) =

1

0

h(x − x)f (x)dx def = Hf (x); 0 ≤ x ≤ 1, (9)

donde f (x) es la senal original, representando en el ejemplo de (Vogel, 2002), la intensidad deuna fuente luminosa con respecto a la posicion x; h, es el kernel de la convolucion, que modelalos efectos de la turbulencia atmosferica en la propagacion de la luz, definido Gaussiano; y g

que es la senal observada, representa la version desenfocada y ruidosa.El modelo matematico que describe al kernel Gaussiano tiene la forma,

h(x − x) =1

σ√

2πexp

− 1

2

(x − x)2

σ2

y f (x) esta dado por la funcion

f (x) =

0,75 ; 0,1 ≤ x ≤ 0,250,25 ; 0,3 ≤ x ≤ 0,32

sin(2πx)4 ; 0,5 ≤ x ≤ 10 para otros

y se le anade un ruido aleatorio con media µ = 0 y desviacion tıpica ση = 10−3, ver (Kerky-

acharian, Picard, & Raimondo, 2005)), BSNR = 40dB. La matriz de convolucion, formadaa partir del kernel, tiene un numero de condicion de 2,5450e + 016, lo que es indicador delmal-condicionamiento del problema inverso.La reconstruccion de la senal usando en Metodo TSVD descrito en la seccion 1.2.1, se mues-tra en Fig.3. El parametro α hace la funcion de filtro, dejando pasar a la solucion de ladeconvolucion, componentes mas contaminados en la medida en que este mas cercano a cero.Sin embargo en la medida que este crece hacia 1,las oscilaciones se parecen mas a funcionessuaves en forma de senos y cosenos. El MSE =0,0046 para el mejor caso en la seleccion delparametro α optimo; con este se obtuvo una SNR= 16,81 dB, y un ISNR = 3,9 dB.El codigo del TSVD fue modificado con la funcion del filtro, para obtener la reconstruccion

de Tikhonov en las mismas condiciones del experimento anterior. El error en la estimaciones menor que en el TVSD (MSE= 0,0044), observe el seguimiento de la senal estimada conrespecto a la original, sobre todo en el componente sinusoidal.En ambos casos (TSVD y Tikhonov) se ha recurrido a un algoritmo iterativo, variando α,para encontrar el mejor caso, de forma similar a como se hizo en (Cui et al., 2007). La sen-sibilidad de ambos metodos al parametro de regularizacion, hace que sean poco robustos acambios de SNR, como se muestra en la Tabla 1.

3http://www.math.montana.edu/∼vogel/Book/Codes/Ch1.

5



10−2

100

0

5

10

α

| | e

α | |

Norma del Error para TSVD

0 0.5 1−0.5

0

0.5

1

1.5

x normalizada

fα

( x )

α = 0.034551

0 0.5 1−2

0

2

4

x normalizada

fα

( x )

α = 0.001

0 0.5 1−0.5

0

0.5

1

x normalizada

fα

( x )

α = 0.83768

Figura 3: Deconvolucion usando TSVD

Cuadro 1: Indicadores SNR, ISNR y MSE para diferentes BSRN = 40, 30, 15 dB.BSNR40 BSNR30 BSNR15

Metodo SNR ISNR MSE SNR ISNR MSE SNR ISNR MSETSVD 16,84 3,92 0,0046 17,08 2,78 0,0059 12,29 1,72 0,0131Tikhonov 17,08 4,16 0,0044 16,24 3,23 0,0053 12,08 1,48 0,0138

Conclusiones

1. El estudio referativo realizado permite tratar el proceso de deconvolucion a partirdel prob lema inverso que este representa como un fenomeno con un marcado mal-condicionamiento en el sentido de Hadamard.

2. El mal-condicionamiento tiene como consecuencia principal la amplificacion del ruidoen la extraccion de la senal, en la misma magnitud que indica el numero de condicionde la matriz del kernel .

3. La descomposicion SVD permite determinar el numero de condicion, en primer lugar,y a partir de la seleccion de los vectores principales , o sea, aquellos que tienen aso-ciado valores singulares grandes, la reconstruccion de la senal dada una observacion

contaminada.

4. La implementacion de la descomposicion truncada de los valores singulares TSVD, sirvecomo un metodo de regularizacion factible para la solucion del mal-condicionamientodel problema inverso.

5. Los resultados experimentales nos muestran la validez de este enfoque para la solucionde problemas inversos en el procesamiento digital, dado la comparacion con el filtro deTikhonov y sus resultados similares.

6



Referencias

Acar, R., & Vogel, C. (1994). Analysis of bounded variation penalty method for ill-posed

problems.

Cui, C., Lamm, P., & Scofield, T. (2007). Local regularization for n-dimensional integralequations with applications to image processing. Inverse Problems, 23 , 1611–1633.(doi:10.1088/0266-5611/23/4/014)

Donoho, D. (1995). De-noising by soft-thresholding. IEEE Trans. Inform., Theory , 613-627.

Hansen, C. (2002). Deconvolution and regularization with toeplitz matrices. Numerical

Algorithms, 29 , 323-378.

Herrera, R., Orozco, R., Moreno, E., & Calas, H. (2005b). Deconvolucion de senales ul-trasonicas por regularizacion wavelet del filtro de Wiener. En Simposio Ing. Elect.

Santa Clara, Cuba. (ISBN: 959250201-3)

Highan, N. J. (1996). Accuracy and stability of numerical algorithms. Dr. Dobb’s Journal.

Kerkyacharian, G., Picard, D., & Raimondo, M. (2005). Adaptive boxcar deconvolution onfull lebesgue measure sets. Journal of the Royal Statistical Society B , 57 (3), 301-369.

Lavarello, R., Kamalabadi, F., & O’Brien, W. (2006). A regularized inverse approach toultrasonic pulse-echo imaging. IEEE Transactions on Medical Imaging , 25 (6), 712–722.

Ma, Q., Wang, X., & Du, S. (2006). Method and application of wavelet shrinkage denoising

based on genetic algorithm. J Zhejiang Univ SCIENCE A, 7 (3), 361–367.

Nash, J. (1990). Compact numerical methods for computers. Adam Hilger.

Neelamani, R., Choi, H., & Baraniuk, R. (2004). ForWaRD: Fourier-wavelet regularizeddeconvolution for ill-conditioned systems. IEEE Trans. Ultrason., Ferroelect., Freq.

Contr., 52 (2), 418-432.

Raimondo, M., & Stewart, M. (2007). The WaveD Transform in R: Performs Fast Translation-Invariant Wavelet Deconvolution. Journal of Statistical Software, 21 (2), 1–28.

Stefan, W., Garnero, E., & Renaut, R. (2005). Signal restoration through deconvolutionapplied to deep mantle seismic probes. Geophysical Journal International , en revision.(http://math.asu.edu/˜rosie)

Vogel, C. (2002). Computational methods for inverse problems. SIAM, Frontier in Applied

Mathematics.

7

El problema lineal inverso en PDS y la descomposición en valores singulares

Documents

Transcript of El problema lineal inverso en PDS y la descomposición en valores singulares