Saçıcıların elastik dalga saçılımı ile görüntülenmesi
Imaging of scatterers through elastic wave scattering
Yazar: UTKU HARMANKAYA
Danışman: PROF. DR. AYŞE KAŞLILAR ŞİŞMAN
İstanbul Teknik Üniversitesi / Lisansüstü Eğitim Enstitüsü / Jeofizik Mühendisliği Ana Bilim Dalı / Jeofizik Mühendisliği Bilim Dalı
Derece: Doktora, 2023
Konu: Jeofizik Mühendisliği = Geophysics Engineering

Özet:
Sismik yöntemler ile yeraltı yapılarının incelenmesi, yer içinde yayılan sismik dalgaların yansıma, kırılma, saçılma ve dispersiyon özellikleri üzerinden değerlendirilmesi ile yapılmaktadır. Yeraltı yapılarının geometrik ve fiziksel özelliklerinin kestiriminde, kayıt edilen sismik dalga alanları üzerinde uygulanan veri işlem, görüntüleme, düz ve ters çözüm tekniklerinden faydalanılmaktadır. Çalışmanın amacına ve ölçeğine bağlı olarak yeraltının yüzeye yakın veya derin kısımları incelenebilir. Yeraltında süreksizlik teşkil eden yapıların sismik yöntemler ile araştırılmasında sıklıkla saçılmış dalgalar kullanılmaktadır. Sismik saçılma, yer içerisinde ilerleyen sismik dalganın baskın dalga boyu ile yeraltında bulunan saçıcının boyutlarının birbirine yakın olduğu koşullarda gerçekleşir. Yeraltı saçıcılarına örnek olarak, tüneller, doğal veya insan yapımı gömülü objeler, boşluklar, yüzeye yakın derinliklerdeki magma sokulumları, sediman havzalar gibi yapılar verilebilir. Bu tür saçıcılardan kaynaklanan saçılmış cisim ve/veya yüzey dalgaları ile de bu yapıların özellikleri incelenebilir. Sediman havzalar petrol, kömür gibi fosil enerji kaynakları, magma sokulumları ise jeotermal kaynaklar açısından önem taşımakta olup, boşluk ve tünel benzeri doğal veya insan yapımı yeraltı yapıları ise bina, yol, tünel, demiryolu ve benzeri inşaat çalışmaları açısından tehlike oluşturabileceği için önem taşımaktadırlar. Bu nedenle, yapılar için tehlike oluşturan yüzeye yakın saçıcıların tespiti ve karakterizasyonu çevre risklerinin azaltılması ve can, zaman ve ekonomik kayıpların önlenmesi açısından önemlidir.
Bu tez çalışması ile amaçlanan, yer içerisinde bulunan yüzeye yakın saçıcıların Lamé ve yoğunluk kontrastlarının birlikte kullanıldığı bir sismik modelleme (düz çözüm) ve görüntüleme (ters çözüm) tekniği geliştirmektir. Burada kullanılacak olan yöntem daha önce sadece yoğunluk kontrastları için geliştirilmiş olup, bu çalışmada Lamé parametre kontrastlarının da eklenmesiyle genişletilmiştir. Sismik dalga yayılımı ortamın katı, sıvı veya gaz gibi fiziksel durumundan etkilendiği gibi, kayacın boşluk oranı, akışkan içeriği, gömülü olduğu derinlik gibi özelliklerinden de etkilenmektedir. Lamé ve yoğunluk kontrastlarının birlikte modellenmesi ve sonrasında ters çözüm ile kestirimi, saçıcının petrofiziksel özellikleri ile ilişkili olan Poisson oranı, Young modülü, Bulk modülü gibi bilgilerin ve sismik hız bilgilerinin elde edilmesine olanak sağlayabilmektedir. Bu bilgiler, saçıcının içinde gömülü bulunduğu artalan ortama göre olan fiziksel yapısı (zayıf/güçlü, boşluk/boşluk olmayan vb.) hakkında detaylı bilgi verebilir. Bu nedenle saçılmış dalga alanının, yoğunluk ve Lamé parametrelerinin birlikte dikkate alınarak modellenmesi ve ters çözümü saçıcının hem petrofiziksel özelliklerinin hem de konumunun daha doğru şekilde kestirilmesine katkı sağlayabilir. 
Tez çalışmasının ilk aşamasında, mevcut durumda saçıcıları sadece yoğunluk kontrastı ile tanımlayan modelleme tekniğinin Lamé parametre kontrastları eklenerek genişletilmesi ele alınmıştır. Kullanılan yöntem, yatay tabakalı ortamda elastodinamik dalga yayılımının Born yaklaşımı ile modellenmesine dayanmaktadır. Dalga alanının integral ifadesi ile hesaplanmasına dayanan bu modelleme tekniğinde, saçıcının gömülü olduğu yatay tabakalı ortamda yayılan dalga alanı iki ayrı dalga alanının toplamı olarak hesaplanmaktadır. Bunlar, saçıcının olmadığı durum için doğrudan ve saçıcının olduğu durum için saçılmış dalga alanlarıdır. Born yaklaşımı nedeniyle, saçılmış dalga alanının hesaplanmasında saçılma her bir saçıcı için tekil saçılma olarak hesaplanmakta, saçıcılar arası çoklu saçılmalar ise ihmal edilmektedir. Green tensörlerinin hesaplandığı modellemelerde, gereken tensörlerin her bir modelleme süreci için baştan hesaplatılabildiği alışılagelmiş modelleme akışı kullanılabildiği gibi, söz konusu tensörlerin önceden bir tablo halinde hesaplatılıp daha sonraki modellerde tablodan çağrılmaları da mümkündür. Bu tabloların kullanımı modellemeleri önemli ölçüde hızlandırsa da artalan parametrelerini ve modelin maksimum uzunluğu ve derinliğini de sabitlemektedir. Bu tez çalışmasında gerçekleştirilen bütün modelleme ve görüntüleme çalışmalarında önceden hesaplanmış Green tabloları kullanılmaktadır. Tez çalışması kapsamında Lamé parametrelerinin türetimlere eklenmesinden önce Fortran 77 kodlarından oluşan orijinal program Fortran 90/95 standardına geçirilmiştir. Programın güncel hesaplama sistemlerinden daha iyi faydalanabilmesi için paralel programlama tekniklerinden de faydalanılmıştır. Lamé kontrastlarının integral ifadelere eklenmesinin ardından, programın doğru sonuç verdiğini teyit etmek için mevcut bir sonlu farklar modelleme programının ürettiği sismogramlar ile kıyaslamalar yapılmıştır. Her iki programda da yarı sonsuz ve tabakalı ortam içine gömülü ve her üç parametrede de kontrast veren saçıcılı modeller için elde edilen sismogramlar iz bazında karşılaştırılmış ve her iki programın da sonuçları arasındaki benzerliğin yüksek olduğu gözlenmiştir. Bu çalışmada geliştirilen modelleme kısmına Lamé parametrelerinin eklenmesi ile birlikte, yüzeye yakın saçıcılar için her bir parametre kontrastının saçılmış dalga alanı üzerindeki etkisini incelemek için ayrıca modelleme çalışmaları da gerçekleştirilmiştir. Elde edilen sonuçlara göre, üç parametre içinde yüzey dalgalarının baskın olduğu sismogramlarda en az etkinin λ parametre kontrastına ait olduğu görülmüştür. Bu gözlem yüzeye yakın saçıcılar için yapılacak görüntüleme çalışmalarında izlenecek yol için önem taşımaktadır. Ayrıca, yüzeyde bulunan bir kaynak-alıcı çifti için her bir parametre kontrastının duyarlılık çekirdekleri (sensitivity kernel) hesaplanmış ve kesitler olarak grafiklenmişlerdir. Bu çekirdeklerden her bir parametre kontrastı için saçıcının konumu ve derinliğinin ters çözüme olan etkisi irdelenmiştir.
Tezin ikinci aşamasında, programın ters çözüm kısmı da öncelikle programlama tekniği açısından benzer bir şekilde geliştirilmiştir. Ters çözüm tekniği olarak iteratif eşlenik gradyan yöntemi kullanılmaktadır. Burada kullanılan iteratif ters çözüm şeması, farklı parametre cinslerini (yoğunluk ve Lamé parametreleri) aynı anda güncelleyecek şekilde genişletilmiştir. Programda birden fazla parametre kontrast cinsi için ters çözüm yapmak amacıyla yüzeye yakın bir ve iki saçıcılı yarı sonsuz modeller için ters çözüm denemeleri gerçekleştirilmiştir. Bir saçıcılı modelde, yoğunluk ve Lamé parametrelerinde %10 kontrast veren ancak bu parametrelerdeki aynı oranda kontrasttan dolayı sismik hızlarda kontrast göstermeyen bir saçıcı kullanılmıştır. Bu, yeraltı yapıları ile ilgili araştırmalarda karşılaşılabilen ve sismik hızlarda zayıf kontrast veren ancak yoğunluk ve Lamé parametrelerinde yüksek kontrastlı yapıları temsil eden bir model olup, yöntemin Lamé parametrelerini kestirmesinin olası bir avantajını da göstermektedir. İki saçıcılı modelde ise saçıcılar içi hava ve su dolu olmak üzere yüksek parametre kontrastı verecek şekilde ayarlanmışlardır. Her iki modelde de ters çözümler her bir parametre için ve parametre çiftlerinden oluşan birleşimler için (ρ-λ, ρ-μ ve λ-μ olacak şekilde) yapılmışlardır. Elde edilen sonuçlarda, tek saçıcılı model için üç parametrede de kontrast oluşturan saçıcının sadece tek parametre ile ters çözümlerinde yoğunluk (ρ) ve μ parametreleri ile saçıcı konumu kestirilebilmiş ancak gerçek parametre kontrastları elde edilememiştir. İki saçıcılı modelde ise saçıcı konumu kestirimlerinde ancak sınırlı bir başarı elde edilmiştir. λ parametresi için yapılan ters çözümlerde ise saçıcı konumun da kestirilemediği görülmüştür. Bunlara karşılık, çift parametreli ters çözümlerden elde edilen sonuçların genel olarak daha iyi olduğu görülmüştür. İki saçıcılı model de her ikisinin de konumlarının başarıyla kestirilebildiği bu örnekler arasında en iyi sonuçların ise ρ-μ parametre çiftinden elde edildiği görülmektedir. Genel olarak, kontrast olarak Born yaklaşımına uymayan saçıcılar (özellikle de iki saçıcılı modeldekiler) için başarılı konum kestirimleri gerçekleştirilmiş ancak gerçek parametre kontrast değerleri yaklaşık olarak kestirilmiştir. Buna karşılık, özellikle iki saçıcılı modelde saçıcıların kontrast parametreleri arasındaki fark ilişkisinin sonuçlara yansıdığı görülmektedir. λ parametresi ile konum ve kontrast kestirimlerinden genelde iyi sonuç alınamamasının nedeni de modelleme kısmında incelenmiş olan izlerden ve duyarlılık çekirdeklerinden görüldüğü üzere parametrenin yüzey dalgalarına olan göreceli sınırlı katkısıdır. 
Bu tez çalışması ile saçıcıların konumlarını, hacimsel koordinatlarını (şekli) ve Lamé ve yoğunluk kontrastlarının belirlenmesini sağlayacak bir görüntüleme tekniği elde edilmiştir. Sismik hız kestirimleri için de geliştirilebilecek bu yöntem için tez çalışmasında verilen örnekler mühendislik ölçeğinde yüzeye yakın saçılma örnekleri olup yöntem daha büyük ölçekli sismolojik çalışmalara da genişletilerek uygulanabilir.

Abstract:
The investigation of the subsurface structures with seismic methods is conducted by recording the propagating P, S, Rayleigh, Love, and similar waves and using their reflection, refraction, scattering and dispersion characteristics. The physical and geometrical properties of the subsurface structures are estimated with the help of signal processing, imaging and inversion techniques which are applied to the wavefield recorded and sampled both in time and spatial domains. For both engineering and global scale studies, shallower or deeper parts of the subsurface can be investigated depending on objective of the study. Scale determines the investigation depth, which in turn determines frequency content and receiver sampling. One of the tools for investigation of subsurface structures are scattered waves. Seismically, scattering occurs when the scatterer and the seismic wavelength have similar dimensions. Using scattered body and/or surface waves, one can obtain information on the physical properties of the scatterers. Examples of near surface scatterers, in descending order of scale, include sedimentary basins, magmatic intrusions, karstic and mining related cavities and buried objects (waste areas, archaeological ruins etc.). Sedimentary basins can be important in terms of energy sources (oil, coal etc.), mantle plumes in terms of geothermal sources, while cavities, archaeological ruins etc. are potential risks to construction of buildings, roads, highways and railways. Detection and monitoring of scatterers in engineering scale is considered to be important in reduction of environmental risks, economic losses and loss of life.
The aim of this thesis is to develop a method for modeling and imaging of subsurface scatterers through their density and Lamé parameter contrasts. The basis of this method was originally established for a previous study concerning only density contrasts of the scatterers. Here, the method is improved with the addition of Lamé parameter contrasts in its formulations for both modelling and inversion parts. As the seismic wave propagation is affected by factors such as the physical state of the media, fluid content, porosity, and depth, modelling and inversion for both density and Lamé contrast parameters may enable us to extract information on scatterer's seismic velocity and its other properties such as Poisson ratio, Young modulus and Bulk modulus. These in turn can provide insight on the physical structure of the scatterer in relation to its surrounding medium (weak/strong zone, cavity/non-cavity etc.). While density characterizes the wave amplitude of the scattered wavefield, Lamé and density parameters together provide phase information on the scattered wavefield. Phase information can contribute to the accuracy of the scatterer location estimation. Therefore, modelling the scattered wavefield by considering both Lamé and density parameter contrasts and using it for the inverse problem would lead to better estimation of the location and petrophysical properties of scatterers. 
This thesis is divided into two main parts. First part deals with the modeling method and related tests, whereas second part deals with the inversion method and imaging examples with newly added parameter contrasts. For the modeling, the method considers three dimensional elastodynamic wave propagation and scattering in a layered medium considering domain-type integral representation with Born approximation. By using the Born approximation, embedded scatterer is assumed have low parameter contrasts to background medium. The wavefield is modelled in separate parts, where one is the incident wavefield calculated for the medium in the absence of scatterers, and the other is the scattered wavefield due to an embedded subsurface scatterer. While the method can handle varying geometries for the aforementioned scatterer, multiple scattering due to interactions of two or more scatterers is omitted due to Born approximation. As the modelling involves the calculations of Green tensors, this can either be modelled 'on the fly', where every requiring Green tensor is calculated in each modeling instance, or they can be calculated separately and saved in a table for any subsequent modelling effort. Using such Green tables significantly speeds up wavefield calculations, although it requires fixed background model parameters (halfspace or layered) and maximum horizontal and vertical distances. All of the modeling and inversion examples in this thesis made use of such pre-calculated Green tables. As the previous version of the method utilized codes written in Fortran 77 standard, here they are updated into Fortran 90/95 standard for better compatibility with modern computing systems. To improve the general performance of the code, parallel computing techniques are also implemented into the program. As the method calculates the wavefields for discrete frequency values, the problem is highly suitable for distributing across multiple processors. Therefore, we were able to obtain good performance scaling with number of utilized processors. Parallel processing was applied to both 'on the fly' wavefield modeling and Green table calculations to accelerate the calculations in both modeling alternatives for our method.
After revamping the program code, the Lamé parameter contrasts are added into the formulations for scattered wavefield modeling. Correct inclusion of Lamé parameters is then tested by wavefield comparisons with the three dimensional finite-difference wavefield modeling program SOFI3D. Two subsurface models, one with higher background velocities and one with lower, are chosen for these tests. For both models, the background medium is a halfspace, while scatterer is set to have contrast in all three parameters. As the density and Lamé parameter contrast have the same ratio, this corresponds to no contrasts in velocities. Comparable wavefields are obtained from using the same subsurface model, source/receiver geometry and source wavelet in both programs. Trace comparisons between the seismograms show agreement with high semblance value for forward and backward scattered wave arrivals, thus affirming the correct inclusion of Lamé parameter contrasts. Any difference between the seismograms is explained by differences in modeling techniques, especially with regards to implementation of attenuation. 
Further testing of our method involves the examination of the effects of different parameter contrast on scattered arrivals. Here, selected traces from scattered wavefields due to contrasts from each parameter is examined in detail. Comparisons are made for both single and pairs of parameters, where it is observed that among the three, the λ has the least, almost negligible, contribution in surface wave scattering due to a near surface scatterer. This is due to dominance of Rayleigh waves in the near surface environment, which are primarily affected by ρ and μ parameters. Another discussion about the parameter effects on wavefields is made over the sensitivity kernels. By calculating and plotting the kernel slices of each parameter for a source-receiver pair, we examine the effect of scatterer positioning on the wavefield. Here, along with the expected decay of Rayleigh wave sensitivities with depth, we observed the insensitivity of surface wave dominated kernels to λ parameter. This is in line with the previous findings from trace analyses, suggesting that for inversions of data gathered for this kind of source/receiver/scatterer geometries, one will need to mostly rely on the contributions of ρ and μ parameters.
For the inversion part of our method, the programming aspects of the code are similarly updated. Inversion process uses an iterative conjugate gradient scheme that updates the model parameter to minimize residual between data and model. Both the modeling and parameter updating subroutines of the inversion process is expanded to update density and Lamé parameter contrasts simultaneously. To test the improved inversion scheme, we use two separate near surface scatterer models. First one has only one near surfacer scatterer with relatively lower (10%) contrast in all three parameters (but no contrast in velocities), while second model has two scatterers (one filled with air, other with water) with strong parameter contrasts. For both subsurface models, the model parameters to be updated during the inversion encompasses a volume that is larger than the scatterers. As is established in other studies concerning full waveform inversions, proper initial models can be required for successful imaging. Therefore, we use initial models that prelocates the scatterers with contrasts that are 25% percent of the actual contrast values. For the initial models of two scatterer case, initial parameters for the scatterers are extended into the volume between the two, thereby starting the inversion with slightly incorrect models. The inversions are done for both single parameter cases and parameter pairs in the forms of ρ-λ, ρ-μ and λ-μ. The unused parameter in each inversion is not updated from its initial model values. Thus, the accuracy of the inversions is expected to be negatively affected by omitted parameters. This is reflected in the results, where single parameter inversions generally have less successful (compared to parameter pair ones) estimates for scatterer locations and contrasts. As for parameter pair inversions, while they are found to be generally better than the single parameter cases, the best results are obtained for the parameter pair ρ-μ. As it was established in the modeling section of this thesis, for this kind of subsurface geometry, λ is the least contributor among the three parameters. Therefore, we observed little improvement from using λ as one of the updated parameters, whereas ρ-μ was sufficient to represent the subsurface structure. From these examples, we conclude that for inversion of surface wave scattering due to near surface heterogeneities, it may be possible to utilize only ρ and μ parameters instead of all three to save computational cost. It should also be noted that even as the initial models have scatterer contrast values in the volume separating the actual scatterers, the method is able to correctly estimate the scatterers as separate entities. While scatterer locations are successfully estimated, future improvements to the inversion could be made for better estimations of the actual contrast values themselves.
With this thesis, a modeling and inversion method is obtained to help determine the location and density and Lamé parameters of a subsurface scatterer. The method, which can also be further developed by implementation of seismic velocity estimations, has the potential to become a cost-effective tool in mitigating hazards in engineering projects and even possibly be expanded to investigations of larger scale scatterers such as imaging sedimentary basins or mantle plumes.