SVILUPPO DI ALGORITMI GENETICI IN LINGUAGGIO C CON APPLICAZIONI A PROBLEMI DI FULL WAVEFORM INVERSION

Un problema geofisico inverso consiste nell’ottenere il modello per il quale i dati predetti meglio si adattino a quelli osservati. Il problema sismico di full waveform inversion mira a calcolare modelli di velocità ad alta risoluzione minimizzando la differenza tra le forme d’onda sismiche sintetiche e quelle osservate. Attualmente esistono diversi strumenti che possono essere utilizzati allo scopo, la maggior parte dei quali sono però distribuiti con licenza di tipo commerciale. Per risolvere il problema di FWI possono essere utilizzati diversi metodi quali metodi di ottimizzazione locale, come ad esempio il metodo Gauss-Newton, che però producono spesso soluzioni intrappolate in minimi locali o metodi di ottimizzazione globale, quali i metodi simulated annealing, neighborhood algorithm o gli algoritmi genetici. Gli algoritmi genetici in informatica (Holland 1992, Holland 1995) sono metodi di ricerca stocastici basati sull’imitazione del processo naturale di evoluzione biologica introdotto da Darwin in “L’origine delle Specie” (Darwin 1859). Facendo riferimento ad alcuni recenti lavori sviluppati sull’argomento (Sajeva et al. 2014) e a librerie Matlab esistenti (Pohlheim 2006) si è voluto sviluppare un software in linguaggio C allo scopo di risolvere il problema sismico di FWI elastica monodimensionale tramite l’utilizzo di algoritmi genetici. Gli obiettivi principali di questo lavoro sono la possibilità di svincolarsi da software proprietario, ottenere prestazioni computazionali migliori, assicurarsi maggiore portabilità e condivisione, ottenere l’opportunità di adeguare gli algoritmi esistenti ad uno specifico uso personale implementando nuove soluzioni non già previste. Il metodo di lavoro affrontato per lo sviluppo di questo codice ha previsto una prima parte di studio teorico degli algoritmi seguita da uno studio approfondito più tecnico delle librerie Matlab esistenti, allo scopo di selezionare quali funzioni fossero necessarie e come integrare le stesse con nuovi metodi per rispettare le nostre specifiche. La parte di sviluppo ha previsto lo studio e la codifica delle funzioni fondamentali degli algoritmi genetici (quali ad esempio la selezione, la ricombinazione, la mutazione, la migrazione) adeguando ogni funzione alle nostre necessità; successivamente a questa fase la codifica è stata estesa ad uno studio delle funzioni oggetto di tipo analitico da utilizzare per una preliminare fase di test. Il software è stato testato su problemi di minimo utilizzando funzioni oggetto analitiche quali le funzioni di De Jong e di Rastrigin e quelle più “complesse” di Schwefel e di Michalewicz e successiva mente è stata effettuata una comparazione del software elaborato con l’algoritmo Matlab esistente, eseguendo una statistica per confrontare l’efficienza dell’algoritmo nei due diversi ambienti al variare dei parametri fondamentali quali la dimensione dello spazio dei modelli, il numero di iterazioni, l’accuratezza. Si è verificato che il software codificato in C permette di ottenere identici risultati con un tempo di elaborazione di quasi 2 ordini di grandezza inferiore. In questa fase di test preliminare si è anche verificato la congruenza dei risultati nei due diversi ambienti, proponendo un metodo di verifica valido caso per caso: impostando per entrambi gli algoritmi la stessa generazione di numeri pseudo-casuali i risultati ottenuti indipendentemente dai due algoritmi risultano identici (alla precisione richiesta di 64-bit, la codifica standard utilizzata durante tutta l’elaborazione). Alla fase di test è seguita una fase di ulteriore sviluppo riguardante l’algoritmo da utilizzare nel caso geofisico di FWI. Si è considerato un modello geologico mono-dimensionale; supponendo di conoscere lo spessore degli strati e le proprietà della colonna d’acqua si è invertito per le velocità P e S e per la densità relative ad 8 strati, quindi con un totale di 21 parametri da invertire. In questa fase si è introdotto l’uso del software OASES (Schmidt 1987) per la soluzione del problema diretto, unendo l’utilizzo dei due programmi tramite linguaggio di Bash Scripting. Durante l’elaborazione ci si è imbattuti in un problema fondamentale di tipo computazionale, problema che si verifica anche in ambiente Matlab; la soluzione del problema diretto richiede la più grossa percentuale di risorse computazionali, rendendo quasi trascurabile il miglioramento dell’efficienza temporale ottenuta con la conversione da Matlab a C. Si è allora introdotta nell’algoritmo una procedura di parallelizzazione. Anziché generare un singolo modello a iterazione vengono generati più modelli contemporanemante indirizzando il calcolo di ogni modello su un processore diverso. La procedura permette una riduzione effettiva del tempo di calcolo che aumenta all’aumentare del numero dei processori. Si è infine utilizzato il codice elaborato per indagare lo spazio dei parametri nel caso geofisico specifico. Sono stati prodotti diversi test, introducendo del rumore per rendere il dato più simile al caso reale e mettendo l’algoritmo in condizioni di ricerca difficili, ad esempio allargando i range di ricerca o modificando la centratura della soluzione rispetto ai limiti di ricerca. I test hanno dimostrato che alcune condizioni rendono l’indagine meno efficiente ma intervenendo con la modifica opportuna di alcuni parametri, quale ad esempio la pressione di selezione, si può affinare la ricerca. Il programma sviluppato presenta i seguenti vantaggi rispetto a software analogo esistente: - è possibile risolvere il problema geofisico di FWI per dati sintetici con prestazioni computazionali migliori - il codice non utilizza nessun tipo di software esterno (a meno del software OASES, che però è distribuito liberamente) e ciò permette di non dover sostenere elevati costi di acquisto di programmi proprietari - la codifica del programma, indipendente dallo sviluppo di altri programmatori, permette di adeguare le funzioni alle nostre ricerche specifiche senza dover richiedere modifiche e integrazioni a terze parti

High-frequency seismic interferometry: broadband measurements of surface-wave phase velocities using “large-N” arrays

High-frequency seismic surface waves sample the top few tens of meters to the top few kilometres of the subsurface. They can be used to determine three-dimensional distributions of shear-wave velocities and to map the depths of discontinuities (interfaces) within the crust. Passive seismic imaging, using ambient noise as the source of signal, can thus be an effective tool of exploration for mineral, geothermal and other resources, provided that sufficient high-frequency signal is available in the ambient noise wavefield and that accurate, high-frequency measurements can be performed on this signal. Ambient noise imaging using the ocean-generated noise at 5-30 s periods is now a standard method, but less signal is available at frequencies high enough for deposit-scale imaging (0.2-30 Hz), and few studies have reported successful measurements in broad frequency bands. Here, we develop a workflow for the measurement of high-frequency, surface-wave phase velocities in very broad frequency ranges. Our workflow comprises (1) a new noise cross-correlation procedure that accounts for the non-stationary properties of the high frequency noise sources, removes bandpass filtering, replaces temporal normalization with short time window stacking, and drops the explicit spectral normalization by adopting cross-coherence; (2) a new phase-velocity measurement method that extends the bandwidth of reliable measurements by exploiting the (resolved) 2π ambiguity of phase-velocity measurements; (3) interstation-distance-dependent quality control that uses the similarity of subgroups of dispersion curves to reject outliers and identify the frequency ranges with accurate measurements. The workflow is highly automated and applicable to large arrays. Applying our method to data from a large-N array that operated for one month near Marathon, Ontario, Canada, we use rectangular subarrays with 150-m station spacing and, typically, 1 hour of data and obtain Rayleigh-wave phase-velocity measurements in a 0.55-23.8 Hz frequency range, spanning over 5.4 octaves, nearly twice the typical frequency range of 1.5-3 octaves in previous studies. Phase-velocity maps and the subregion-average 1D velocity models they constrain show a high-velocity anomaly consistent with the known, west-dipping gabbro intrusions beneath the area. The new structural information can improve our understanding of the geometry of the gabbro intrusions, hosting the Cu-PGE Marathon deposit

Seismic imaging at a mineral-deposit scale using highfrequency surface waves (0.5–24 Hz) in ambient noise wavefield

Ambient noise, surface-wave tomography (ANSWT) is now a routine technique for imaging crustal and upper mantle structure at a regional scale. Its cost efficiency and environmental friendliness also make ANSWT an attractive method for mineral exploration. However, the application of the technique in mineral exploration requires the retrieval of wide-band, high-frequency surface waves from seismic noise, so as to obtain highresolution images of shallow structures. We present a new workflow optimized to extend the bandwidth of highfrequency surface waves retrieved. It comprises short time-window stacking, cross-coherence and an improved phase velocity measurement method. We tested the workflow on data from a large-N array over a Cu-PGE deposit in Ontario, Canada, and successfully measured phase velocities for numerous inter-station pairs in a broad frequency ranges from 0.55 Hz to 23.8 Hz. Analysis of the phase velocity maps reveals a west-dipping highvelocity anomaly that matches the west-dipping, multi-staged gabbro intrusions associated with the deposit

Optimal resolution tomography with error tracking and the structure of the crust and upper mantle beneath Ireland and Britain

The classical Backus–Gilbert method seeks localized Earth-structure averages at the shortest length scales possible, given a data set, data errors, and a threshold for acceptable model errors. The resolving length at a point is the width of the local averaging kernel, and the optimal averaging kernel is the narrowest one such that the model error is below a specified level. This approach is well suited for seismic tomography, which maps 3-D Earth structure using large sets of seismic measurements. The continual measurement-error decreases and data-redundancy increases have reduced the impact of random errors on tomographic models. Systematic errors, however, are resistant to data redundancy and their effect on the model is difficult to predict. Here, we develop a method for finding the optimal resolving length at every point, implementing it for surface-wave tomography. As in the Backus–Gilbert method, every solution at a point results from an entire-system inversion, and the model error is reduced by increasing the model-parameter averaging. The key advantage of our method stems from its direct, empirical evaluation of the posterior model error at a point. We first measure inter- station phase velocities at simultaneously recording station pairs and compute phase-velocity maps at densely, logarithmically spaced periods. Numerous versions of the maps with varying smoothness are then computed, ranging from very rough to very smooth. Phase-velocity curves extracted from the maps at every point can be inverted for shear-velocity (V S ) profiles. As we show, errors in these phase-velocity curves increase nearly monotonically with the map roughness. We evaluate the error by isolating the roughness of the phase-velocity curve that cannot be explained by any Earth structure and determine the optimal resolving length at a point such that the error of the local phase-velocity curve is below a threshold. A 3-D V S model is then computed by the inversion of the composite phase-velocity maps with an optimal resolution at every point. The estimated optimal resolution shows smooth lateral variations, confirming the robustness of the procedure. Importantly, the optimal resolving length does not scale with the density of the data coverage: some of the best-sampled locations display relatively low lateral resolution, probably due to systematic errors in the data. We apply the method to image the lithosphere and underlying mantle beneath Ireland and Britain. Our very large data set was created using new data from Ireland Array, the Irish National Seismic Network, the UK Seismograph Network and other deployments. A total of 11 238 inter-station dispersion curves, spanning a very broad total period range (4–500 s), yield unprecedented data coverage of the area and provide fine regional resolution from the crust to the deep asthenosphere. The lateral resolution of the 3-D model is computed explicitly and varies from 39 km in central Ireland to over 800 km at the edges of the area, where the data coverage declines. Our tomography reveals pronounced, previously unknown variations in the lithospheric thickness beneath Ireland and Britain, with implications for their Caledonian assembly and for the mechanisms of the British Tertiary Igneous Province magmatism

SISMIKO:emergency network deployment and data sharing for the 2016 central Italy seismic sequence

At 01:36 UTC (03:36 local time) on August 24th 2016, an earthquake Mw 6.0 struck an extensive sector of the central Apennines (coordinates: latitude 42.70° N, longitude 13.23° E, 8.0 km depth). The earthquake caused about 300 casualties and severe damage to the historical buildings and economic activity in an area located near the borders of the Umbria, Lazio, Abruzzo and Marche regions. The Istituto Nazionale di Geofisica e Vulcanologia (INGV) located in few minutes the hypocenter near Accumoli, a small town in the province of Rieti. In the hours after the quake, dozens of events were recorded by the National Seismic Network (Rete Sismica Nazionale, RSN) of the INGV, many of which had a ML > 3.0. The density and coverage of the RSN in the epicentral area meant the epicenter and magnitude of the main event and subsequent shocks that followed it in the early hours of the seismic sequence were well constrained. However, in order to better constrain the localizations of the aftershock hypocenters, especially the depths, a denser seismic monitoring network was needed. Just after the mainshock, SISMIKO, the coordinating body of the emergency seismic network at INGV, was activated in order to install a temporary seismic network integrated with the existing permanent network in the epicentral area. From August the 24th to the 30th, SISMIKO deployed eighteen seismic stations, generally six components (equipped with both velocimeter and accelerometer), with thirteen of the seismic station transmitting in real-time to the INGV seismic monitoring room in Rome. The design and geometry of the temporary network was decided in consolation with other groups who were deploying seismic stations in the region, namely EMERSITO (a group studying site-effects), and the emergency Italian strong motion network (RAN) managed by the National Civil Protection Department (DPC). Further 25 BB temporary seismic stations were deployed by colleagues of the British Geological Survey (BGS) and the School of Geosciences, University of Edinburgh in collaboration with INGV. All data acquired from SISMIKO stations, are quickly available at the European Integrated Data Archive (EIDA). The data acquired by the SISMIKO stations were included in the preliminary analysis that was performed by the Bollettino Sismico Italiano (BSI), the Centro Nazionale Terremoti (CNT) staff working in Ancona, and the INGV-MI, described below

Le attività del gruppo operativo INGV "SISMIKO" durante la sequenza sismica "Amatrice 2016",

SISMIKO è un gruppo operativo dell’Istituto Nazionale di Geofisica e Vulcanologia (INGV) che coordina tutte le Reti Sismiche Mobili INGVPublishedLecce3T. Sorgente sismica4T. Sismicità dell'Italia8T. Sismologia in tempo reale1SR TERREMOTI - Sorveglianza Sismica e Allerta Tsunami2SR TERREMOTI - Gestione delle emergenze sismiche e da maremoto3SR TERREMOTI - Attività dei Centr

Hot Upper Mantle Beneath the Tristan da Cunha Hotspot From Probabilistic Rayleigh-Wave Inversion and Petrological Modeling

Understanding the enigmatic intraplate volcanism in the Tristan da Cunha region requires knowledge of the temperature of the lithosphere and asthenosphere beneath it. We measured phase-velocity curves of Rayleigh waves using cross-correlation of teleseismic seismograms from an array of ocean-bottom seismometers around Tristan, constrained a region-average, shear-velocity structure, and inferred the temperature of the lithosphere and asthenosphere beneath the hotspot. The ocean-bottom data set presented some challenges, which required data-processing and measurement approaches different from those tuned for land-based arrays of stations. Having derived a robust, phase-velocity curve for the Tristan area, we inverted it for a shear wave velocity profile using a probabilistic (Markov chain Monte Carlo) approach. The model shows a pronounced low-velocity anomaly from 70 to at least 120 km depth. VS in the low velocity zone is 4.1-4.2 km/s, not as low as reported for Hawaii (∼4.0 km/s), which probably indicates a less pronounced thermal anomaly and, possibly, less partial melting. Petrological modeling shows that the seismic and bathymetry data are consistent with a moderately hot mantle (mantle potential temperature of 1,410-1,430°C, an excess of about 50-120°C compared to the global average) and a melt fraction smaller than 1%. Both purely seismic inversions and petrological modeling indicate a lithospheric thickness of 65-70 km, consistent with recent estimates from receiver functions. The presence of warmer-than-average asthenosphere beneath Tristan is consistent with a hot upwelling (plume) from the deep mantle. However, the excess temperature we determine is smaller than that reported for some other major hotspots, in particular Hawaii

Broadband Surface Wave Tomography of Ireland, Britain and Other Regions

Over the last decades, seismic surface-wave studies have produced increasingly detailed images of the Earth s structure at a regional scale. In this study, we have tuned well-established techniques and when required implemented new ones in order to investigate regions in which important debates are still ongoing, regarding the structure and the evolution of the Earth beneath them. Several studies suggested that the Paleogene uplift of parts of Britain and Ireland was caused by a lateral branch of the Iceland mantle plume, which played a fundamental role in the evolution of the North Atlantic Ocean over the past 60 M.y. Alternatively, among competing hypothesis, it was suggested that the uplift could be due to the far-field stress associated with the Alpine and Pyrenees Orogenies, with reactivation of old Variscan and Caledonian faults across Ireland and Britain. A major part of this study is aimed at gaining new insights into the seismic structure of the British Isles, which can help us answer these open questions. Teleseismic earthquakes and ambient noise, recorded at densely spaced seismic stations in the region, were used to determine the surface-wave dispersion across the British Isles and construct detailed images of the seismic structure beneath the area. The measurements, obtained using independent surface-wave analysis techniques (cross-correlation of teleseismic surface waves, multimode waveform fitting, and ambient noise interferometry), were applied to produce the first 3D shear-velocity model of the lithosphere and the asthenosphere of the entire region including Ireland, Britain, and the Irish Sea. The application of different methodologies yielded complementary frequency bands of the measurements, sensitive to different depths, from the shallow crust to the deep upper mantle. Abundant, newly available data was used to image the region with higher resolution than previously. The highly uneven station coverage resulted in a considerably irregular distribution of the measurements in the area; this, and the effects of errors on the measurements, required the development of a new, multi-resolution tomographic scheme. This scheme allows us to maximize the information extracted from the data and reach an optimal target resolution of the model at each knot, minimizing the effects of uneven data sampling and of the propagation of systematic errors.The multi-resolution phase-velocity maps, obtained at densely spaced periods, were inverted, point by point, for shear-velocity structure in order to produce a 3D, shear-velocity model of the lithosphere and asthenosphere. The optimal resolution tomography offers important new insights into the structure and evolution of the British Isles. A robust, low-velocity anomaly beneath the Irish Sea and its surroundings persists in the models from ~60 to at least 140 km depth, indicating an anomalously thin lithosphere. The area that exhibits the low velocity anomaly corresponds to where uplift and volcanism are evidenced by geological data. Our results also show a striking correlation with proposed underplating thickness and denudation, gravity, and thermochronological measurements, and rule out the once common assumption of a constant lithospheric thickness across Britain and Ireland. At lithospheric depths, a clear contrast in seismic velocities between Ireland and Britain could possibly explain why the seismicity is nearly absent in Ireland, while modest but clearly higher in Britain. The higher velocities beneath most of Ireland indicate thicker lithosphere and colder geotherms, likely resulting in a higher-strength lithosphere, resisting deformation. In the lithospheric mantle, the model displays an elongated high-velocity anomaly stretching W-E approximately along the Iapetus Suture Zone in Ireland, which may be the expression of a remnant of the Caledonian Iapetus slab beneath the suture or, alternatively, fragments of thick continental lithosphere incorporated into the Irish landmass in the course of the Caledonian Orogeny. Another part of this study was on using surface-wave analysis to investigate the lithosphere-asthenosphere system beneath the Tristan da Cunha Hotspot, with the goal of understanding the enigmatic intraplate volcanism in the region. Surface-wave analysis was applied in a challenging setting, as this work involved the use of data recorded by ocean-bottom seismometers, which required data-processing and measurement approaches substantially different from those tuned for land-based arrays of stations. We constrained a region-average, shear-velocity structure, using two-station, cross-correlation measurements across the area, and inferred the temperature of the lithosphere and asthenosphere beneath the area by means of petrological modeling. Seismic inversion and petrological modeling show a lithospheric thickness of only 65 70 km, confirming the previous estimates obtained from receiver functions. Our observations are consistent with a hot plume from the deep mantle, but the excess temperature estimated is much smaller than that reported for some other major hotspots, in particular, Hawaii

Seismic Thermography

ABSTRACT Variations in temperature within the Earth are of great interest because they indicate the thickness and, consequently, mechanical strength of the lithosphere and density variations and convection patterns in the sublithospheric mantle. Seismic tomography maps seismic velocity variations in the mantle, which strongly depend on temperature. Temperatures are, thus, often inferred from tomography. Tomographic models, however, are nonunique solutions of inverse problems, regularized to ensure model smoothness or small model norm, not plausible temperature distributions. For example, lithospheric geotherms computed from seismic velocity models typically display unrealistic oscillations, with improbable temperature decreases with depth within shallow mantle lithosphere. The errors due to the intermediate-model nonuniqueness are avoided if seismic data are inverted directly for temperature. The recently developed thermodynamic inversion methods use computational petrology and thermodynamic databases to jointly invert seismic and other data for temperature and composition. Because seismic velocity sensitivity to composition is much weaker than to temperature, we can invert seismic data primarily for temperature, with reasonable assumptions on composition and other relevant properties and with additional inversion parameters such as anisotropy. Here, we illustrate thus-defined seismic thermography with thermal imaging of the lithosphere and asthenosphere using surface waves. We show that the accuracy of the models depends critically on the accuracy of the extraction of structural information from the seismic data. Random errors have little effect but correlated errors of even a small portion of 1% can affect the models strongly. We invert data with different noise characteristics and test a simple method to estimate phase velocity errors. Seismic thermography builds on the techniques of seismic tomography and relies on computational petrology, but it is emerging as a field with its scope of goals, technical challenges, and methods. It produces increasingly accurate models of the Earth, with important inferences on its dynamics and evolution.</jats:p