• CORRELAZIONE (AUTO e CROSS)

• DECONVOLUZIONE

• MODELLO CONVOLUZIONALE DEGLI STRATI

• ATTENUAZIONE MULTIPLE PRIMARIE

Il grado di similitudine tra due funzioni f1(t) ed f2(t) può essere misurato moltiplicando le due funzioni tra loro e integrando il risultato. Se le due funzioni hanno fasi molto diverse risulterà un numero negativo, ma se vogliamo una informazione connessa con la forma dei segnali converrà ripetere la moltiplicazione e l’integrazione dopo avere shiftato un segnale rispetto all’altro di un particolare shift temporale.

Questa operazione si chiama Cross – Correlazione:

AUTOCORRELAZIONE E CROSS - CORRELAZIONE

dtfftx )()()( 21

Se f1 = f2 = f la misura di similarità di f rispetto a se stessa si chiama A(t) = Autocorrelazione

Ovviamente, la funzione di A avrà valore massimo a tempo di shift uguale a 0, ma se Il segnale ha spiccate periodicità, la funzione di A mostrerà dei massimi a shift Temporali uguali a questi periodi. La trasformata di Fourier della funzione di Aautocorrelazione è la potenza spettrale del segnale stesso

• La convoluzione:

dtfftx )()()( 21

La convoluzione è strutturalmente simile alla correlazione sono diverse per il segno – di

Sistema(f2)f1 x

Esempi convoluzione/cross - correlazione

Date le sequenze a e b In Matlab:

AUTOCORRELAZIONE E CROSS - CORRELAZIONE

a=[2 2 2 2 2];b=[3 3 3 ]; c=conv(a,b)c = 6 12 18 18 18 12 6 figure(1) plot(c) stem(c) holdCurrent plot held plot(c)

c=xcorr(a,b)c = Columns 1 through 6 4.4409e-016 4.4409e-016 6.0000e+000 1.2000e+001 1.8000e+001 1.8000e+001 Columns 7 through 9 1.8000e+001 1.2000e+001 6.0000e+000 figure(2) plot(c) holdCurrent plot held stem(c)

DECONVOLUZIONE

Se la sorgente fosse un pulso ideale, una serie di pulsi riflessi molto stretti sarebbero ricevuti dalle varie discontinuità di impedenza del sottosuolo.

Si è precedentemente visto che le sorgenti reali danno segnali complessi, che saranno riflessi a ciascuna interfaccia.

Questa struttura dei segnali tenderà ad oscurare la registrazione dei treni d’onda provenienti dalle varie riflessioni.

Possiamo considerare la forma d’onda riflessa ricevuta come la convoluzione della forma d’onda di ingresso con la risposta all’impulso del sottosuolo. Per riottenere la risposta alla forma d’onda ideale potremo “de-convolvere”

Tale deconvoluzione equivale a trovare un filtro che opera sulla forma d’onda della sorgente e la trasforma in un singolo pulso (filtro inverso).

Un filtro, come abbiamo visto, è un sistema nel quale avviene la convoluzione tra il segnale d’ingresso e la risposta all’impulso del filtro.

Nel caso dei segnali sismici reali più che funzioni analogiche, si opera su sequenze numeriche, che sono la rappresentazione dei segnali elaborati dai computer. Qui useremo una rappresentazione simbolica compatta.Obiettivo: Progettare un operatore D (filtro numerico) che, quando viene convoluto, con il segnale S produce un singolo impulso :

D*S=operatore di convoluzione

Il segnale riflesso osservato R sarà la convoluzione della risposta del sottosuolo (earth) all’impulso E (insieme dei segnali riflessi) con il segnale sorgente:

R=E*S

E

S R

Quindi se si applica l’operatore D al segnale osservato (ricevuto) si ha:

D * R = D * E * S = E * D * S = E * = E (Risposta del suolo)

Dove abbiamo sfruttato la proprietà commutativa della convoluzione ed anche il fatto che la convoluzione di un segnale con l’impulso dà il segnale stesso.

Qundi:

Se progettiamo un operatore D che riesca a collassare il segnale sorgente in un singolo pulso, possiamo applicarlo al segnale riflesso per rimuovere gli effetti della forma del segnale estraendo così la vera risposta del sottosuolo.

Cerchiamo quindi il migliore operatore, in senso statistico, di lunghezza finita, che minimizzi lo scarto quadratico medio tra il segnale ideale e quello ottenuto. Tale operatore è dato dalla soluzione del set di equazioni:

)()(1

jfJ xx

n

Jxz

J = 0, 1………….n

xx = funzione di autocorrelazione di ingrassof = coefficienti dell’operatore convoluzionexz = cross-correlazione tra ingresso e uscita idelale

Nel seguito della presentazione troviamo il procedimento analitico per la determinazione dell’ operatore (filtro inverso) che applicato al segnale emesso (wavelet source reale o firma) riesce a collassarlo in un pulso ideale ( diverso da zero solo in un punto).

Ora daremo solo la parte applicativa:Supponiamo di avere un wavelet source WS definito dalla sequenza di valori WS= [0 7 -3 1] ; vedi Nota

Supponiamo di aver ottenuto, con il procedimento statistico (minimi quadrati), riportato nelle pagine seguenti

l’ operatore D=[0.145 0.062 0.008]; v.Nota

Eseguendo ora la convoluzione di D con WS otteniamo il WI ideale: WI=D*WS ovvero: WI=conv(D,WS); v Nota,

WI= 0 1.0150e+000 -1.0000e-003 1.5000e-002 3.8000e-002 8.0000e-003

Nota:le espressioni in colore possono essere incollate in Matlab ed eseguite in linea.I comandi per eseguire i grafici

Per ottenere i grafici; figure(1);plot(WS); hold;stem(WS);figure(2);plot(WI);hold;stem(WI)

Supponiamo di voler calcolare un operatore di lunghezza 2 (n=1), il cui effetto sia per quanto possibile una spike (impulso) nell’origine.

Tenendo conto che la funzione di autocorrelazione è di tipo pari )1()1( xxxx

0)1()0()1(

1,0

0

7)0()1()0(

1,0

0

10

10

xzxxxx

xzxxxx

ff

J

ff

J

Avremo quindi:

Calcoliamo i valori dell’autocorrelazionexx (0) = 7 x 7 + (-3) x (-3)+1 x 1 = 59 xx (1) = 7x (-3) +(-3 x 1) = -24 xx (2) = 7 x 1 = 7

Sostituendo nel sistema: +59 f0 -24 f1 = 7 f0 = 0.142-24 f0 +59 f1 = 0 f1 = 0.058

Esempio numerico: una forma d’onda (7, -3, 1) che può essere assimilata ad un segnale reale:

7

3

1 1 0 0

0 t

Cosicché l’operatore richiesto è (0.142, 0.058), che produce una forma d’onda in uscita:

D = (0.142, 0.058)

S = (7, -3, 1)

SD

Eseguendo la convoluzione abbiamo:t = 0 7 x 0,142 = 0,99t = 1 (-3) x 0.142 +7 x 0.058 = 0.02t = 2 (1) x 0.142 +(-3) x 0.058 = 0,03t = 3 (1) x 0.058 = 0.06

Che è una buona approssimazione dell’impulso. Teniamo conto che:

...0,, 10 aaA ...0,,, 210 bbbB

La convoluzione è y = A*B

0

)(k

kn knbay n = 0, 1, 2, 3, 4…

....

....

....

....

1221303

0211202

1201101

2211000

bababay

bababay

bababay

bababay

t=0t=1

t=2

t=3

b0 b1 b2

a0a1

=a0b0

a0a1

=a1b0 + a0b1

a0a1

=a1b1 + a0b2

a0a1

= a1b2

0

)(k

kn knbay n = 0, 1, 2, 3, 4…

Se scegliessimo un operatore di lunghezza 3 invece che 2, riusciremmo a concentrare meglio l’energia in una singola spike, avremmo le seguenti equazioni:

+59 f0 -24 f1 +7 f2 = 7

-24 f0 +59 f1 -24 f2 = 7

+7 f0 -24 f1 +59 f2 = 7

D= (0.145, 0.062, 0.008)L’effetto di questo operatore sulla forma d’onda è:D*S = (1.015, 0.015, 0.038, 0.008)

Quindi conoscendo la forma d’onda della sorgente possiamo trovare un opportuno operatore che, applicato alle tracce sismiche, ci fornisce la migliore risoluzione, cioè la risposta dell’impulso. Risulta quindi necessario conoscere la funzione di autocorrelazione del segnale sorgente. Ma considerato che la sequenza dei coefficienti di riflessione di un suolo stratificato può essere considerata random e che il rumore può essere trascurato (rapporto segnale – rumore molto grande), la funzione di autocorrelazione della generica traccia (segnale di uscita) è uguale alla funzione di autocorrelazione del segnale sorgente. In realtà, a causa delle modificazioni che il segnale ha nell’attraversamento degli strati, quali l’assorbimento delle alte frequenze, sarà necessario calcolare la funzione di autocorrelazione a tratti lungo la traccia sismica.

TRASFORMATA DI FOURIER

dtetfc tj

)(

2

1)(

TRASFORMATA DI FOURIER, LAPLACE E TRASFORMATA Z

Esiste se If(t)I è integrabile(vincolo importante)

TRASFORMATA DI LAPLACE

dtetfsc st 0

)()(

s j

dteetfsc tjt 0

))(()(

La e-t va abbastanza velocemente a 0 per cui: dtetf t

0)(

Quasi tutti i segnali fisici soddisfano la precedente relazione.

Nel discreto, cioè nel caso di segnali campionati uniformemente, che sono quindi costituiti da sequenze di numeri è applicata la trasformata Z:

0

)(n

nnzxzX

TRASFORMATA Z

Valida per IzI > R e stabile per R<1 e con js eez

Una sequenza:

[xn] = [x0, x1, …………………. XM]

Ha come trasformata Z:

MMn zxzxzxxxz ......)( 2

21

10

I problemi che possono risolvere con la deconvoluzione sono anche:

1) GHOSTS

MULTIPLE PRINCIPALI, SECONDARIE E GOST su TRACCE SISMICHE

superficie

sorgenteghost

diretto

ghost

diretto

ricevitore

Il segnale ghost (fantasma) è uguale al segnale utile, nella forma, può avere quasila stessa intensità e fase invertita rwair= 1 rwbottomed un ritardoche può essere:- al minimo 2 per distanza ricevitore - superficie, - al massimo quest’ultima più 2 volte distanza sorgente superficie. Per distanze intendiamo i percorsi ottici dei raggi.I segnali ghost si possono eliminare sia con una regolazione geometrica, cioè portando la sorgente ed il ricevitore quasi in superficie ovvero con un filtraggio, sia silenziando il ricevitore per un tempo opportuno.

Multipla secondaria

Multipla primaria

2) MULTIPLE PRIMARIE

La multipla si genera da una seconda riflessione

superficie

mare

fondo

S A B

r

-r

Ricordiamo che A e B coincidono e che i raggi sono inclinati solo per comodità di visualizzazione, poiché in realtà sono verticali.Nelle registrazioni le multiple s riconoscono perché sono posizionate ad una profondità,che è esattamente il doppio di quella dell’eco diretta ed inoltre, le pendenze sonoraddoppiate.Per fondali profondi le multiple sono stampate al di sotto del segnale utile, ma nel caso di bassi fondali, sono sovrapposte al segnale e tendono a mascherarlo. Anche le multiple posso essere filtrate, ovviamente non con un filtro tradizionale (di banda), perché esso eliminerebbe anche il segnale utile, ma con uno numerico.

3)MULTIPLE SECONDARIEPiu difficili da eliminare ma piu’ deboli.

-1 -1

-r2r2

r

Eliminazione delle multiple tramite operatore DT

Anche in questo caso se vogliamo tralasciare teoria e formalismi possiamo affrontare in modo speditivo il problema (Jones ,pag85).Calcoliamo la risposta all’ impulso dello strato d’ acqua compreso tra aria e fondale . Consideriamo che la sequenza WT dei valori ottenuti dall’ idrofono e quindi dal registratore (dopo conversione A/D) supponendo che il coefficiente di riflessione acqua-aria e’ -1, quello acqua-fondale è r e la prima riflessione ha ampiezza A (arbitraria).I coefficienti sono 1 impulso primario, 2 prima multipla di riflessione, 3 seconda multipla etc.,l’ intervallo di tempo tra i due segnali è il tempo di andata ritorno attraverso lo strato d’ acqua.WT=[ (A) (-2*A*r) (3*A*r^2) (-4*A*r^3) (5*A*r^4)];

La sequenza inversa di WT ha solo elementi diversi da zero:

DT=[ (1) (2*r) (r^2)];

convolvendo le due sequenze:T=DT*WT=[A 0 0 0 0];

ovvero in Matlab :T=conv(DT,WT); con A=1:T=[1 0 0 0 0 0 0];

L’approccio nel dominio del tempo é utile sia quando il sistema e descritto in termini di equazioni differenziali sia quando ne è definita la risposta impulsiva. In particolare é utile la descrizione temporale nei casi in cui i dati sperimentali sono rappresentati da serie temporali, cioè dalle tracce sismiche in cui l’ordinata é il segnale trasdotto dal rivelatore e l'ascissa il tempo (discretizzato nei sistemi di acquisizione digitali).Il modello convoluzionale del suolo consiste nel considerare che il sismogramma registrato sia la convoluzione dei segnale sorgente sn, con la risposta all’impulso unitario del terreno hn, più un termine di rumore o disturbo d(t). Si ammettono per solito le ipotesi di cui in appresso:

a) il processo sismico (attraversamento di uno o più strati di suolocon riflessioni, soddisfa la teoria di propagazione delle onde elastiche neimezzi isotropi ed omogenei, senza perdite e non dispersivi e quindi conconservazione della forma durante Ia propagazione stessa.b) alla rappresentazione continua del segnale

che tiene conto delle successive riflessioni di energia alle varie discontinuità, si sostituisce quella convertita in digitale (sequenza numerica} nel rispetto del teorema del campionamento. Sara allora t = k t essendo t il tempo di Campionamento o tempo tra due campioni successivi e k = 1, 2, 3, …..

0

)()()(n

nnn tdtshty (1)

c) Il ritardo in tempo n può essere rappresentato da un multiplo intero dit. Cosicchè n =nt; n = 0, 1, 2, … per cui la 1 si sostituisce con la seguente:

0

)()()(n

nn tkdtnkshtky

kytky )( nkstnks )(

(2)

definendo

e la kt = k, la 15 diviene:

0nknknk dshy

Il sismogramma (yk) inteso come sequenza numerica è la convolzione del segnale sorgente (sk) con la risposta all’impulso del terreno (hk) più il rumore (dk)

+sorgente

sk

hk = risposta

all’impulso del sist. Acqua-Fondale

xk

dk = rumore

Yk = sismogramma

Fig 3 – schematizzazione corrispondente alla la convolzione del segnale sorgente (sk) con la risposta all’impulso del terreno (hk) più il rumore (dk)

Tralasciamo la presenza del rumore e consideriamo il modello matematico di suolo costituito da un solo strato omogeneo ed isotropo (acqua) tra l’aria e il bedrock yk = xk , fig. 3

H1

Sorgente / rilevatoreAria z00

r0=1

r

1

Acqua z01 = 1v1

Suolo z02 = 2v2

i = densità dello strato (i = 1, 2, ….);Vi = velocità di propagazione (onde e s) nello strato(i = 1, 2, ….);ri = coefficiente di riflessione all’interfaccia strato acqua – fondale (sedimento – roccia);r0 = coefficiente di riflessione all’interfaccia tra aria e lo strato di acqua;ti = coefficiente di trasmissione 1+r1;H = spessore dello strato;ri, r0, ti, sono numeri reali e Ir0I e IriI <=1 mentre ItiI <=2.

Fig 3

Fig 3

Supponendo di applicare come sorgente sk un impulso k, avremo al rilevatore dellafig 5, nel quale i percorsi sono rappresentati inclinati solo per chiarezza, la situazione in appresso indicata.

Sorgente S k=k

ritardo 2n 4n 6n

r0

r1

r0 r0

r1 r1 r1

riv. = xk

Fig. 5

La situazione è quindi la medesima della fig. 4 cioè incidenza normale.. Considerandoche n = nt è il tempo di attraversamento dello strato H, ed si ha essendo r0=-1

....;63140

2121 nknknkk rrrrx (4)

Per cui il comportamento del suolo può essere descritto da una struttura di ritardi e somme.

Eseguendo la trasformata di Laplace z che consente una agevole rappresentazione simbolica dei ritardi convolutivi, si ottiene:

....;)( 620

31

40

21

21

nnn zrrzrrzrzx (5)

Se r0 = +1 come si verifica per le vibrazioni se il mezzo superiore è l’aria, posto 2n=si ha:

....;)( 31

2211

zrzrzrzx (6)

La (18) è operativamente equivalente alla struttura riportata in fig 6, corrispondente ad un filtro generalizzato di tipo FIR (Finite Impulse Response) non ricorsivo (la trasf. z dik è 1)

z- z- z-

r1 r12 r1

3

+ +

1

X(z)

La funzione X(z) rappresenta la risposta all’impulso del sistema suolo. Ponendo in evidenza r1z si ottiene:

.......])(1[)( 2111

zrzrzrzX (7)

serie di potenze in r1 ed in z.Se ora invece dell’impulso unitario si considera come sorgente un generico segnale sk

si avrà:

......6314

2121 nknknkk srsrsrx (8)

Trasformando nel dominio delle Z abbiamo

........)(1)( 2111 zrzrzrzx

)()()()( 331

2211 zszrzszrzszrzx

ovvero

X(z) = S(z) H(z) ; H(z) = “risposta Z”

........)(1)(

)()( 2

111 zrzrzrzs

zxzH

Poiché 1 è negativo I1I<1, posto II abbiamo:

Poiché la somma della serie:

1........1 432 xperxxxx

è

x11

abbiamo

)(1

)(

1)(

zSzr

zrzX

zr

zrzH

quindi

Noi vorremmo eliminare tutti i contributi dovuti alla riverberazione che corrispondonoa lasciare solo , infatti la risposta ideale sarebbe:zr

)(1

)(

)()('

zSz

zzX

zSzzX

(A) ; la reale è

Cerchiamo un operatore F(z) tale che

F(z) X(z) = X’(z)

Dividendo membro a membro (A) e (B)

(B)

)()1()(

1)(

)(

'

'

zXzrzX

zrzX

zX

Quindi l’operatore cercato è

zrzF 1)(