Full Waveform Spectral IP: implementazione e valutazione di algoritmi di processing ed inversione, locale e globale.

Il lavoro svolto ha riguardato la messa a punto e la valutazione di algoritmi per l’elaborazione, e la loro inversione, di dati acquisiti con un sistema di Polarizzazione Indotta Full-Waveform. La Polarizzazione Indotta (IP) è una metodologia che trova frequente, se non sistematica applicazione nell’ambito dell’esplorazione mineraria; recentemente questa metodica ha trovato applicazione anche in contesti diversi, quali quello idrogeologico e quello ambientale. La IP Spettrale utilizza la risposta dispersiva del semispazio all’applicazione di una corrente variabile nel tempo, al fine di recuperare le caratteristiche tessiturali e mineralogiche del corpo o dei corpi causativi. Nella Polarizzazione Indotta classica, sia time che frequency-domain, la curva di corrente non è registrata ma viene assunta a priori; la curva di tensione viene campionata ad intervalli relativamente ampi ed un parametro di sintesi, la polarizzabilità, viene calcolato in maniera diretta. La Polarizzazione Indotta Spettrale permette di meglio descrivere il comportamento dispersivo del semispazio attraverso, ad esempio, la relazione di Pelton (Pelton et al., 1978) che è una variante dell’equazione di Cole-Cole. Si tratta di una funzione che descrive la resistività complessa e i cui parametri possono essere messi in relazione con le proprietà del sottosuolo: \rho(\omega)=\rho_{0}\prod_{i=1}^2\left[1-m_{i}\left(1-\frac{1}{1+(j\omega\tau_{i})^{c_{i}}}\right)\right] Dove \rho_0 è la resistività in continua, m_i è la caricabilità, c_i è la frequency dependence, \tau_i è la costante di tempo e ω è la frequenza angolare. Il valore massimo assunto dall’indice i (generalmente non maggiore di 2) determina il tipo di dispersione: singola o doppia. I valori di m, τ e c vengono ottenuti tramite inversione degli spettri di ampiezza e fase, e sono i parametri che possono consentire di caratterizzare ulteriormente il semispazio, almeno rispetto alla polarizzabilità l.s. In generale però, specialmente lavorando nel dominio dei tempi, i dati misurati consentono analisi limitate in termini di ampiezza di banda e soffrono delle assunzioni aprioristiche riguardanti la curva di corrente. È stato recentemente mostrato (Ollson et al., 2015) che registrando e prendendo in considerazione sia la curva di corrente che la curva di tensione, utilizzando perciò l’approccio full-waveform IP, si può estendere la band-width utile del segnale e quindi ricavare, dall’inversione, una stima più robusta dei cosiddetti parametri Cole-Cole. Con l’approccio FW-IP si considera il semispazio come un sistema lineare in cui la funzione di trasferimento sarebbe l’equazione (1), l’input del sistema è la corrente e l’output è la tensione. Tramite una semplice divisione spettrale è possibile quindi ricavare la funzione di trasferimento. Su dati sintetici senza rumore l’operazione non presenta particolari problemi; in presenza di rumore sia il calcolo della funzione di trasferimento che la successiva inversione risulteranno affetti da effetti spuri. Considerato che sui dati reali può essere presente un livello elevato di rumore, sono stati implementati, in Matlab, alcuni step di elaborazione preliminari volti a migliorare il rapporto segnale/rumore. Questi sono stati testati su serie sintetiche a cui è stato sovrapposto rumore reale, registrato cioè sul campo con configurazioni elettrodiche appropriate. Gli algoritmi implementati sono il detrending, che permette di sottrarre alla curva di tensione un trend dovuto per esempio ad effetti spuri quali polarizzazione di elettrodo oppure ai potenziali spontanei, e lo stacking, effettuato nel dominio spettrale mediando, dopo correzione di fase per lo shift temporale, le trasformate di Fourier dei singoli cicli di registrazione. Una volta ottenuti gli spettri, relativi ad un singolo ciclo dell’onda quadra utilizzata, si passa alla divisione spettrale tra lo spettro della tensione e lo spettro della corrente, il tutto moltiplicato per un coefficiente geometrico determinato dalla geometria di acquisizione. Vengono così ottenuti degli spettri di resistività che, nonostante il pre-processing, mostrano ancora, specialmente alle alte frequenze, distorsioni dovute al rumore. Per migliorare ulteriormente il rapporto segnale/rumore si effettua quindi una decimazione passando dallo spettro DFT, le cui frequenze sono spaziate in maniera lineare, ad uno spettro con spaziatura logaritmica. La decimazione è stata effettuata tramite convoluzione dei valori spettrali con una finestra di Parzen. Oltre a questa finestra, sono state utilizzate ulteriori funzioni di pesatura, calcolate tramite il reciproco dello spettro di potenza del rumore (normalizzato) e la funzione di coerenza quadratica calcolata a partire dagli spettri di tensione e corrente. Il risultato finale consiste di spettri in cui l’effetto del rumore alle alte frequenze (es. maggiori di 100 Hz) è fortemente attenuato, e quindi adeguati come input per gli algoritmi di inversione. Per il problema dell’inversione sono stati utilizzati due approcci diversi; un approccio locale tramite l’algoritmo Levenberg-Marquardt (LM), implementato ex-novo in Matlab, ed un approccio globale, Particle Swarm Optimization (PSO), utilizzando un algoritmo del Matlab Optimization Toolbox. Come detto la relazione di Pelton può essere utilizzata con modello di dispersione singolo (i=1) o doppio (i=2). Per quanto riguarda l’inversione con un singolo modello, sia l’algoritmo LM che il PSO risultano efficaci, e la scelta fra i due dipenderà quindi da considerazioni di efficienza computazionale. L’inversione con un doppio modello, dunque la stima di 7 parametri, risulta complicata dalla effettiva band-width a disposizione, e dipende, significativamente, da quanto sono diverse le τ reali, cioè da quanto i due modelli sono separati negli spettri. Dato che l’algoritmo locale può rimanere “intrappolato” in minimi locali, si è cercato di ridurre tale rischio cercando di stimare dagli spettri, tramite algoritmi, sia i bounds che lo starting model. Tuttavia, nel caso di doppio modello l’inversione locale sembra faticare a convergere, soprattutto nel caso in cui i due modelli non siano abbastanza separati, cioè qualora si abbia uno spettro di fase costituito da un solo picco. Data la difficoltà intrinseca di definire la complessità del problema tramite l’analisi della funzione oggetto in 7 dimensioni si è deciso di testare l’algoritmo globale PSO. Confrontando i risultati dei due algoritmi si è quindi dedotto che il problema, nel caso di doppio modello, non risiede tanto nella presenza di minimi locali quanto piuttosto in un problema di equivalenza. Per accertare meglio ciò si è utilizzata una decomposizione ai valori singolari che ha mostrato come, all’aumentare del numero dei parametri e con due modelli molto ravvicinati, siamo in presenza di un problema di forte mal condizionamento. In particolare, l’analisi degli autovettori dello spazio dei modelli ha mostrato che alcuni parametri influenzano molto poco il dato osservato, il che si traduce nel fatto che diverse combinazioni degli stessi parametri fitteranno in maniera equivalente il dato. Le considerazioni di cui sopra sono state effettuate lavorando su dati sintetici contaminati da rumore reale; per testare gli algoritmi su dati reali è stata effettuata un’acquisizione in un test-site presso Fenice Capanne nel comune di Massa Marittima (Gr). Questo sito è stato scelto per la presenza di un giacimento a pirite e solfuri misti sfruttato sino alla fine degli anni ’70 e per l’elevata presenza di rumore (lo stesso utilizzato per contaminare i dati sintetici). La geometria di acquisizione scelta è del tipo polo-dipolo con distanza inter-elettrodica di 25 metri. È stato utilizzato un transmitter IP generante un’onda quadra tipo time-domain, con periodo di 8 s, la quale è stata registrata tramite uno shunt posizionato lungo il circuito di corrente. Il primo canale (canale zero) del receiver (16 canali a 24-bit Simultaneous-Sampling) ha registrato la curva di corrente, mentre gli altri canali hanno registrato le curve di tensione agli elettrodi di tensione. La qualità del dato, a causa della strumentazione in parte sperimentale (shunt e receiver, che hanno richiesto alcuni aggiustamenti duranti le fasi di acquisizione) non ha permesso di effettuare un’intera linea di acquisizione, così si sono avuti a disposizione solo i dati registrati dal canale 0, che misura la curva di corrente, e dal canale 1 che registra la curva di tensione ad una distanza di 25 metri dall’elettrodo di corrente. I dati acquisiti permettono comunque di arrivare ad alcune conclusioni. Gli step di processing messi a punto risultano efficaci in termini di miglioramento del rapporto segnale/rumore, mentre l’applicazione dell’approccio full-waveform ha consentito, a partire dai dati reali, di ottenere degli spettri di ampiezza e di fase con una buona risoluzione alle alte frequenze. I due algoritmi di inversione applicati hanno confermato le osservazioni effettuate con i dati sintetici e cioè che l’inversione Cole-Cole per un singolo modello non ponga alcun problema, né in termini di convergenza, né in termini di scelta del modello di partenza, mentre per l’inversione con un doppio modello si osservano problemi di equivalenza causati dal mal condizionamento del problema. Quindi, considerata la probabile assenza di minimi locali, o l’influenza relativamente ridotta che questi possono avere sulla convergenza, risulta più efficace ed efficiente, con questo tipo di dati, l’utilizzo di un algoritmo di inversione locale piuttosto che globale