DiallelX: Sveliamo i Segreti dei Terremoti con il Fortran Moderno!

Amici appassionati di scienza e tecnologia, preparatevi perché oggi vi porto nel cuore pulsante della ricerca sismica, un campo dove i dati sono sovrani e la capacità di analizzarli rapidamente può fare la differenza. Parleremo di microterremoti, quei piccoli sussulti della Terra che, sebbene spesso impercettibili per noi, sono miniere d’oro di informazioni per capire la geometria delle faglie attive, la loro attività e persino i processi che precedono grandi terremoti. E come facciamo a scovare questi flebili segnali in un mare di rumore? Con una tecnica chiamata Matched Filter Analysis (MFA), ma soprattutto, con strumenti software all’avanguardia. Ed è qui che entra in gioco il protagonista della nostra storia: DiallelX.

La Sfida: Cacciare Microterremoti in Oceani di Dati

Immaginate di dover setacciare anni e anni di registrazioni sismiche continue, provenienti da decine di sensori, alla ricerca di segnali che assomigliano a migliaia di “eventi modello” (i nostri template). Un lavoro da far tremare i polsi, vero? La chiave è calcolare il cosiddetto coefficiente di cross-correlazione di rete (NCC), una misura di quanto le forme d’onda continue assomiglino ai nostri template. Se l’NCC supera una certa soglia, bingo! Abbiamo un potenziale evento.

Il problema è che questo processo, per tutte le combinazioni possibili, è incredibilmente dispendioso in termini di tempo e risorse computazionali. Pensate che calcolare l’NCC per ogni canale richiede operazioni dell’ordine di (O(M log {W})), dove M è il numero di campioni nei dati continui e W la lunghezza delle finestre di analisi. Moltiplicate questo per K canali e capirete perché analizzare, ad esempio, 10 anni di dati con migliaia di template può richiedere supercomputer con centinaia di GPU e costi esorbitanti, come nel caso dello studio di Ross et al. (2019).

Certo, si potrebbe pensare di ridurre la frequenza di campionamento dei dati per alleggerire il carico, ma questo va bene per terremoti a bassa frequenza. Per i piccoli terremoti “regolari”, quelli con frequenze dominanti sopra i 10 Hz, un downsampling troppo aggressivo significherebbe perdere preziose informazioni, rendendoli invisibili. C’era bisogno di una soluzione più agile, accessibile anche a chi non dispone di infrastrutture da capogiro.

DiallelX: Una Ventata d’Aria Fresca nel Mondo Fortran

Ed ecco che entro in scena io, o meglio, il team di cui faccio parte, con DiallelX. Abbiamo sviluppato questo codice in Fortran moderno (sì, il Fortran è ancora vivo e vegeto, e più performante che mai in certi ambiti!) proprio per affrontare questa sfida. L’idea di base? Invece di puntare al calcolo esatto dell’NCC, che è costoso, DiallelX ne calcola un’approssimazione. “Un’approssimazione?”, direte voi. “E l’accuratezza?”. Tranquilli, le nostre indagini statistiche hanno dimostrato che questa approssimazione è più che sufficiente per i nostri scopi, e i vantaggi sono enormi!

DiallelX è progettato per essere:

• Efficiente nell’uso della memoria: Niente più workstation che implorano pietà.
• Con un output di dimensioni ridotte: Dimenticatevi file giganteschi che intasano i dischi.
• Facilmente parallelizzabile: Sfrutta al massimo i moderni processori multi-core.

Il codice è open source e disponibile su GitHub, perché la scienza è condivisione! Il nome “DiallelX” deriva da un termine della genetica che indica tutte le possibili ricombinazioni tra diversi genotipi; un po’ come noi che incrociamo tutti i dati con tutti i template.

Come Funziona la Magia di DiallelX (Senza Diventare Matti)

Per capire DiallelX, dobbiamo fare un piccolo passo indietro e parlare di MASS (Mueen’s Algorithm for Similarity Search), un algoritmo veloce per calcolare le funzioni di cross-correlazione (CF) usando la Trasformata Discreta di Fourier (DFT). MASS divide la forma d’onda continua in finestre parzialmente sovrapposte. DiallelX prende spunto da questo, ma con una differenza cruciale: mentre MASS punta al calcolo esatto, noi, come detto, approssimiamo.

In DiallelX, la lunghezza della finestra di analisi (W) è uguale alla lunghezza del template (L). Calcoliamo la CF nel dominio della frequenza. Certo, trattare la finestra come una funzione periodica introduce un’approssimazione, ma se c’è un segnale sismico nella finestra, i valori di CF saranno comunque significativamente più alti rispetto a quelli tra il template e il semplice rumore.

E qui, amici, arriva una delle magie di DiallelX, specialmente per i dati multi-canale (e in sismologia, abbiamo spesso almeno 3 canali per stazione, e più stazioni per localizzare un evento, quindi facilmente 12 o più canali!):
La formula per l’NCF approssimato è:
$$ text{NCF}_{text{approx}} = text{IDFT} left[ sum_{k=1}^{K} frac{widehat{u}^{(k)} cdot overline{widehat{v}^{(k)}}}{left| u^{(k)} right|_{L^2} cdot left| v^{(k)} right|_{L^2}} right] $$
Notate qualcosa? La sommatoria su K (i canali) è dentro la Trasformata Inversa di Fourier (IDFT). Questo significa che possiamo ridurre il numero di IDFT da K a 1, normalizzando prima nel dominio della frequenza e scambiando l’ordine delle sommatorie. Poiché le IDFT sono computazionalmente non trascurabili, questo trucchetto ci fa guadagnare un sacco di tempo, specialmente con molti canali. Immaginate: con 12 canali, riduciamo i calcoli di IDFT a un dodicesimo!

Non Affoghiamo nei Dati: Ottimizzazione dell’Output

Un altro problema con l’analisi di grandi moli di dati è… la mole di dati in output! Se salvassimo tutti i massimi di NCF per ogni combinazione finestra-template, ci ritroveremmo con file enormi. Per esempio, 24 ore di dati a 100 Hz e 10.000 template potrebbero generare oltre 1 GB di dati di massimi! DiallelX adotta due strategie furbe:

1. Per ogni finestra, salviamo solo il risultato del template che le assomiglia di più. L’output diventa indipendente dal numero di template!
2. Confrontando i risultati di finestre adiacenti (che spesso si sovrappongono e contengono lo stesso pacchetto d’onda), salviamo solo quello con l’NCC più alto.

Il risultato? Un file di output finale (i “candidati”) che, anche per anni di dati e migliaia di template, raramente supera i 200 MB. Geniale, no?

C’è un piccolo dettaglio tecnico da gestire: la periodicità introdotta dalla DFT può creare ambiguità sulla posizione esatta del picco di correlazione. DiallelX risolve questo problema ricalcolando i candidati con una tecnica chiamata “zero padding” (aggiungendo zeri alla fine del segnale), ma solo per i candidati già selezionati, così da non appesantire troppo i calcoli iniziali.

Il Parametro “a”: Trovare il Giusto Equilibrio

Quando si usa DiallelX, si può specificare un parametro opzionale ‘a’, che è il rapporto tra la lunghezza della finestra W e lo “stride” (delta) (cioè di quanto si sposta la finestra ad ogni passo). Un valore di ‘a’ più grande significa più sovrapposizione tra le finestre, il che probabilisticamente migliora l’accuratezza nel trovare il punto esatto di inizio del segnale. Però, aumenta anche il numero di finestre e quindi il costo computazionale. Il valore ottimale di ‘a’ dipende dai dati; di solito un valore come 4 offre un buon compromesso, garantendo un’autocorrelazione superiore a 0.8 con alta probabilità in casi pratici.

Numeri alla Mano: Prestazioni e Confronti

Abbiamo messo alla prova DiallelX su una workstation con una CPU AMD Ryzen Threadripper 3990X. I risultati sono stati entusiasmanti! Per esempio, analizzare un giorno di dati (8.64 milioni di campioni) con 1000 template da 1024 campioni, su 32 canali, ha richiesto circa 32 minuti usando 120 thread. La cosa notevole è che il tempo di esecuzione non cresce linearmente con il numero di canali, grazie alla furbata della riduzione delle IDFT.
Nello studio di Naoi et al. (2025), DiallelX è stato usato per analizzare dati da esperimenti di idrofratturazione: nel caso più imponente, l’analisi di dati equivalenti a 5.7 anni di registrazioni sismiche a 100 Hz con oltre 11.000 template è stata completata in 18 ore! Questo apre la porta all’analisi di interi dataset ventennali, come quelli della rete Hi-net in Giappone.

E rispetto ad altri codici come SEC-C (una versione multi-canale di MASS basata su MATLAB)? Beh, un confronto diretto è difficile perché SEC-C è in MATLAB/Python e DiallelX in Fortran. Ma, per darvi un’idea, su dati identici, DiallelX in single-thread è risultato oltre 30 volte più veloce di SEC-C in MATLAB. E il consumo di RAM? SEC-C ha richiesto circa 68 GB, DiallelX solo 2.6 GB! Questa differenza è abissale e rende DiallelX utilizzabile su workstation comuni, anche per dataset impegnativi.

Orizzonti Futuri: DAS e GPU

Le potenzialità di DiallelX non si fermano qui. Stiamo pensando a come potrebbe rivoluzionare l’analisi dei dati da sensori DAS (Distributed Acoustic Sensing), quelle fibre ottiche che possono trasformarsi in migliaia di sensori sismici. I segnali DAS sono spesso molto rumorosi, ma trattandoli come un sistema multi-canale, l’NCC può emergere dal rumore, e l’efficienza di DiallelX sarebbe cruciale per gestire la vastità di questi dati.
E il futuro? Sicuramente il supporto GPU. La metodologia di DiallelX si presta bene al calcolo parallelo su GPU, anche se dovremo affrontare la sfida della velocità di comunicazione con la memoria di sistema, magari con tecnologie come GPUDirect Storage.

Insomma, DiallelX si è dimostrato un vero e proprio game-changer per l’analisi di grandi dataset sismici, permettendo di identificare un numero significativamente maggiore di eventi rispetto ad altre tecniche, come dimostrato nello studio di Naoi et al. (2025) dove ha superato anche approcci basati sul deep learning. È la dimostrazione che con algoritmi intelligenti e un buon vecchio Fortran ottimizzato, si possono ancora fare grandi cose nella scienza!

Fonte: Springer