Img

IPOGEO - Inverse Problems and Optimization in Geophysics

ipogeo


Indice

Modelli matematici per la prospezione geofisica

La prospezione geofisica utilizza tecniche non invasive e non distruttive per indagare le caratteristiche del sottosuolo ed eventuali disomogeneità. Tali tecniche sono di grande importanza in molti ambiti: ingegneria civile, identificazione di inquinanti nell'ambiente, agricoltura di precisione, localizzazione di mine inesplose in zone di guerra, esplorazione di siti archeologici sepolti, etc.

Figura 1: Stima della conducibilità elettrica di una zona del Parco di Molentargius (Cagliari) sulla base di misure elettromagnetiche.
Image fetta

Il terreno viene attraversato da un campo primario di onde sismiche o elettromagnetiche, a bassa o alta frequenza (georadar), provocate artificialmente. Gli strumenti misurano il campo secondario distorto dal terreno e consentono di dedurre le proprietà del sottosuolo (densità, conducibilità elettrica, etc.) da un confronto con il campo primario.

I modelli matematici che descrivono questi fenomeni sono spesso classici (equazioni delle onde o di Maxwell), ma richiedono approcci diversi da quelli usuali perché non interessa tanto l'andamento della soluzione, quanto la struttura del mezzo di propagazione. Questo conduce a modelli integrali, sistemi di equazioni nonlineari e, più in generale, problemi inversi.

Risoluzione numerica di equazioni integrali

Un'equazione integrale è un'equazione in cui la soluzione f è una funzione che compare sotto il segno di integrale. Si parla di equazione integrale lineare di prima specie se essa si presenta solo sotto il segno di integrale

$\displaystyle \int_{{D}}^{}$k(x, y)f (x)dx = g(y),    y∈D.

Si fa riferimento, invece, a un'equazione integrale lineare di seconda specie se la soluzione è presente anche fuori l'integrale

f (y) + $\displaystyle \int_{{D}}^{}$k(x, y)f (x)dx = g(y),    y∈D.

In entrambi i casi le funzioni k e g sono dette nucleo e termine noto, rispettivamente, e D rappresenta il dominio della funzione f (un intervallo finito, un intervallo illimitato, una curva, ...).

Diversi problemi applicativi, come ad esempio il semplice problema di Dirichlet o Neumann per l'equazione di Laplace, si possono riformulare in termini di equazioni integrali. Nelle applicazioni, la soluzione f è il segnale da ricostruire; la funzione g, di solito nota solo in un set di punti, rappresenta i dati sperimentali; il nucleo k invece, di solito nota analiticamente, rappresenta la risposta all'impulso dell'apparecchiatura sperimentale.

Figura 2: Funzioni rappresentanti di Riesz per un particolare sistema di equazioni integrali risolte in uno spazio di Hilbert con nucleo di riproduzione (RKHS).
Image etaHcap

La ricerca mira a sviluppare metodi numerici per approssimare la soluzione di equazioni integrali di prima e seconda specie. Tale sviluppo si articola nei seguenti step:

  • Sulla base della regolarità delle funzioni coinvolte, si sceglie lo spazio funzionale in cui approssimare la soluzione.
  • Si approssima l'operatore integrale mediante una opportuna formula di quadratura o cubatura e si valuta l'equazione risultate in determinati punti. Si ottiene così un sistema di equazioni. La scelta dello schema di quadratura dipende dalle eventuali patologie delle funzioni coinvolte.
  • Si studia teoricamente il metodo in termini di stabilità e convergenza, fornendo, se possibile, anche delle stime teoriche dell'errore.
  • Si effettua una vasta sperimentazione numerica al fine di verificare se le stime teoriche sono confermate dalle prove numeriche.

Problemi inversi e metodi di regolarizzazione

Ci troviamo di fronte ad un problema inverso quando vogliamo ricostruire un segnale a partire da misurazioni indirette, legate al primo tramite una qualche relazione funzionale. Un esempio è fornito dalle equazioni integrali di prima specie. In generale, questo tipo di problema ammette soluzioni multiple ed è molto sensibile alla presenza delle inevitabili perturbazioni presenti nei dati. Inoltre, si tratta spesso di problemi di grandi dimensioni e dall'esoso costo computazionale.

Figura 3: Sinogramma prodotto dalla trasformata di Radon, modello utilizzato nella TAC, applicata ad un oggetto modello.
Image IP

Lo sviluppo di metodi numerici per la loro risoluzione che siano al tempo stesso, accurati, stabili rispetto alla presenza di errori nei dati, e computazionalmente efficienti, è un compito impegnativo, ma interessante e per il quale esiste una forte richiesta dalle scienze applicate. Gli algoritmi richiedono particolari tecniche di regolarizzazione, destinate a rendere i problemi ben posti e ridurre l'influenza della propagazione degli errori. Tali tecniche dipendono da parametri che devono essere stimati in modo opportuno.

Tecniche di ottimizzazione

Figura 4: Funzioni Griewank e Michalewics.
Image FigLera2 Image FigLera3

I problemi di ottimizzazione abbondano nella maggior parte dei campi della scienza, dell'ingegneria e della tecnologia. In molti di questi problemi è necessario calcolare l'ottimo globale (o una buona approssimazione) di una funzione multivariabile, spesso multiestremale, non differenziabile e difficile da valutare. Le variabili che definiscono la funzione da ottimizzare possono essere continue e/o discrete e, inoltre, spesso soddisfano determinati vincoli. Il problema generale si definisce come

Determinare x*∈S, tale che F(x*)≤F(x), ∀x∈S,

dove F : S→R è una funzione assegnata in S⊆RN.

I problemi di ottimizzazione globale appartengono alla classe di complessità dei problemi NP-hard. Tali problemi sono molto difficili da risolvere. I tradizionali algoritmi di ottimizzazione di discesa basati su informazioni locali non sono adeguati per risolvere questi problemi. Nella maggior parte dei casi di interesse pratico il numero di ottimi locali aumenta, in media, in modo esponenziale con la dimensione del problema (numero di variabili).

Figura 5: Funzione GKLS e curva di Peano.
Image FigLera1 Image FigLera4

Esistono diverse possibili tecniche di risoluzione.

  • Approccio deterministico. Garanzia di successo. Teoremi di convergenza. Condizioni, a volte restrittive, per la funzione obiettivo.
  • Approccio stocastico. Generalmente, ipotesi molto deboli sulla funzione obiettivo. Ampliamento della classe delle funzioni obiettivo. Non assoluta garanzia di successo del metodo. Convergenza in probabilità.
  • Approccio euristico. Spesso riescono a trovare la soluzione globale. Non riescono a dimostrare teoricamente l'ottimalità della soluzione trovata!

Un possibile approccio deterministico introdotto da Strongin nel 1972 prevede l'utilizzo delle curve “riempi-spazio” di Peano per ridurre la dimensione N del problema originale.

In particolare Strongin ha dimostrato che il minimo globale di un problema N-dimensionale Lipschitziano è uguale al minimo globale della funzione monodimensionale

f (x) = F(p(x)),   x∈[0, 1],

dove p(x) è la curva di Peano. Inoltre la funzione f diventa Hölderiana

| F(x) - F(y)|≤H|| x - y||1/N,

con H = 2L$\sqrt{{N+3}}$, costante di Hölder.

In Fig. 4 possiamo vedere una approssimazione di livello 5 della curva di Peano-Hilbert in 2 dimensioni. Per quanto riguarda l'utilizzo delle curve riempi-spazio di Peano, nelle recenti ricerche sono state utilizzate due differenti approssimazioni della curva: la prima, la piecewise-linear approximation, più comunemente considerata nell'ottimizzazione globale, e un secondo tipo, detta non-univalent approximation, meno conosciuta, che ci ha permesso di accelerare il processo per la determinazione dell'ottimo globale.

Figura 6: 3-RPR robot, a sinistra, e Delta-robot, a destra.
Image 3-RPR_robot Image delta_robot

Recentemente le suddette tecniche sono state applicate per la risoluzione di sistemi di disequazioni per la determinazione delle aree di lavoro di robot industriali paralleli. La Fig. 4 mostra due esempi di robot paralleli di cui abbiamo individuato l'area di lavoro (definita fissando eventuali parametri). La Fig. 4 riporta una rappresentazione grafica della area di lavoro del robot 3-RPR, in particolare a sinistra abbiamo uno zoom che evidenzia la differenza tra le due approssimazioni della curva utilizzate e a destra la variazione dell'area di lavoro al variare dei parametri.

Figura 7: Aree di lavoro del robot 3-RPR.
Image robot1 Image robot2

Estrapolazione e stima di parametri

Un problema che si ritrova spesso nelle applicazioni è quello di determinare il limite di una successione (scalare o vettoriale) {sn}n∈$\scriptstyle \mathbb {N}$. Nella realtà, i dati a disposizione sono chiaramente finiti, pertanto è necessario stimare il limite conoscendo solo N elementi s1,..., sN.

Le tecniche di estrapolazione si possono interpretare come una generalizzazione della più nota interpolazione. Tuttavia, mentre nell'interpolazione i dati sono utilizzati per stimare il comportamento di una funzione all'interno dell'intervallo in cui sono contenuti, con l'estrapolazione si studia il comportamento della funzione all'esterno di tale intervallo.

Figura 8: Risoluzione regolarizzata dell'equazione integrale di Phillips: la soluzione di Tikhonov basata su un metodo di estrapolazione risulta più accurata di quella ottenuta mediante la curva-L.
Image phill

Lo studio delle tecniche di estrapolazione ha lo scopo di una risoluzione più efficiente e veloce e, data la sua formulazione estremamente generale, ha un campo di applicazione molto vasto. Nel caso dei problemi inversi, essa può essere utilizzata per stimare il parametro di regolarizzazione e per migliorare la qualità di una soluzione approssimata.

Metodi numerici per la prospezione FDEM

Uno dei campi di applicazione in cui CaNA ha lavorato più a lungo riguarda la prospezione elettromagnetica nel dominio della frequenza (FDEM). Si tratta di una tecnica non-invasiva basata sull'uso di un Ground Conductivity Meter, strumento relativamente poco costoso e semplice da usare. Esso crea, tramite una bobina trasmittente, un campo elettromagnetico primario che attraversando il terreno genera, per induzione elettromagnetica, un campo secondario in eventuali conduttori presenti nel sottosuolo. Il rapporto tra i campi secondario e primario viene misurato da una bobina ricevente.

Esistono modelli integrali sia lineari che nonlineari che legano le caratteristiche elettromagnetiche del suolo alle misure. Entrambi sono malcondizionti e propagano in modo distruttivo gli errori di misura.

Figura 9: Raccolta dati di induzione elettromagnetica presso l'aeroporto di Elmas.
Image gppe_elmas

Il modello lineare è stato studiato con vari approcci matematici. Quello più efficace fa ricorso a spazi di Hilbert con nucleo di riproduzione (RKHS). Per il modello nonlineare, particolarmente complicato, è stata dapprima trovata l'espressione dello Jacobiano rispetto alla conducibilità e alla permeabilità magnetica. È stato poi sviluppato un algoritmo risolutivo regolarizzato basato sul metodo di Gauss-Newton e la singular value decomposition (SVD). Il software Matlab sviluppato è stato pubblicato in forma open source.