Dziękujemy autorowi bloga za jego wystąpienie
Wprowadzenie
Niniejsza notatka jest napisana po konferencji i zawiera streszczenie wystąpienia, które wygłaszałem. Wydarzenie odbyło się w czerwcu w Petersburgu. Biorąc pod uwagę, że pracuję kwartał od miejsca konferencji, nie mogłem tego wydarzenia pominąć. W 2016 i 2017 roku słuchałem wykładów konferencji, a w tym roku wystąpiłem jako prelegent. Po pierwsze, pojawił się interesujący (wg mnie) temat, który rozwijamy z , a po drugie, po długim studiowaniu przepisów prawa RF w zakresie polityki sankcyjnej, w firmie, w której pracuję, pojawiły się aż dwie licencje .
Zanim przejdę do tematu mojego wystąpienia, chciałbym podkreślić dobrą organizację wydarzenia. Na stronie rejestracyjnej konferencji wykorzystano obraz Kazan Cathedral. Katedra jest jednym z głównych zabytków Petersburga i jest doskonale widoczna z sali, w której odbywała się konferencja.
Na wejściu do SPbGEU uczestników witały asystentki wśród studentów – nie pozwalały się zgubić. Podczas rejestracji rozdawano małe upominki (zabawka – migający spajk, długopisy, naklejki z logo Wolfram). Obiad i przerwa na kawę były również uwzględnione w harmonogramie konferencji. O smacznym kawie i ciastkach już wspominałem na stronie grupy – kucharze spisali się na medal. Tą wstępną częścią chciałbym podkreślić, że samo wydarzenie, jego format i lokalizacja przynoszą już pozytywne emocje.
Wykład, który przygotowałem razem z Kiryłem Biełowem, nazywa się «Wykorzystanie Wolfram Mathematica do rozwiązywania problemów zastosowanej geofizyki. Analiza spektrogramu danych sejsmicznych lub „gdzie płynęły starożytne rzeki”. Treść wykładu obejmuje dwie części: po pierwsze, to wykorzystanie algorytmów dostępnych w do analizy danych geofizycznych, a po drugie, jak umieścić dane geofizyczne w Wolfram Mathematica.
Sejsmika
Na początek warto zrobić mały wstęp do geofizyki. Geofizyka to nauka badająca fizyczne właściwości skał. Ponieważ skały mają różne właściwości: elektryczne, magnetyczne, elastyczne, istnieją odpowiednie metody geofizyki: elektrobadania, magnetobadania, sejsmiki... W kontekście tego artykułu szczegółowo omówimy tylko badania sejsmiczne. Badania sejsmiczne są główną metodą poszukiwania ropy naftowej i gazu. Metoda ta opiera się na wzbudzaniu fal elastycznych i następnej rejestracji odpowiedzi od skał stanowiących badaną powierzchnię. Wzbudzenie fal odbywa się na lądzie (dynamitem lub nieekspansyjnymi źródłami wibracyjnymi fal elastycznych) lub na morzu (pneumatycznymi działkami). Fale elastyczne rozprzestrzeniają się przez warstwy skał, załamując się i odbijając na granicach warstw o różnych właściwościach. Odbite fale wracają na powierzchnię i są rejestrowane przez geofony na lądzie (zazwyczaj są to urządzenia elektrodynamiczne oparte na ruchu magnesu zawieszonego w cewce) lub hydrofony w morzu (oparte na efekcie piezoelektrycznym). Na podstawie czasu przybycia fal można ocenić głębokość warstw geologicznych.
Statek sejsmiczny ciągnie sprzęt
Działka pneumatyczna wzbudza fale elastyczne
Fale przechodzą przez warstwy skał i są rejestrowane przez hydrophony
Statek badawczo-naukowy do geofizyki 'Iwan Gubkin' przy molo u Mostu Błagowieszczeńskiego w Petersburgu
Model sygnału sejsmicznego
Skały mają różne właściwości fizyczne. Do badań sejsmicznych szczególnie ważne są właściwości elastyczne — prędkość rozprzestrzeniania się fal elastycznych oraz gęstość. Jeśli dwie warstwy będą miały identyczne lub podobne właściwości, fala 'nie zauważy' granicy między nimi. Jeżeli natomiast prędkości fal w warstwach będą się różnić, na granicy warstw nastąpi odbicie. Im większa różnica w właściwościach, tym intensywniejsze będzie odbicie. Jego intensywność będzie określana współczynnikiem odbicia (rc):
gdzie ρ — gęstość skał, ν — prędkość fal, 1 i 2 oznaczają górną i dolną warstwę.
Jednym z najprostszych i najczęściej stosowanych modeli sygnału sejsmicznego jest model konwolucyjny, w którym zarejestrowana ścieżka sejsmiczna przedstawiana jest jako wynik konwolucji sekwencji współczynników odbicia z impulsem sondującym:
gdzie s(t) — ścieżka sejsmiczna, tzn. wszystko, co zarejestrował hydrofon lub geofon w czasie trwania rejestracji, w(t) — sygnał generowany przez pistolet pneumatyczny, n(t) — szum losowy.
Obliczmy przykładową sztuczną ścieżkę sejsmiczną. Jako sygnał początkowy wykorzystamy szeroko stosowany w sejsmice impuls Ricker'a.
length=0.050; (*Długość sygnału*)
dt=0.001; (*Częstotliwość próbkowania sygnału*)
t=Range[-length/2,(length)/2,dt]; (*Czas sygnału*)
f=35; (*Częstotliwość centralna*)
wavelet=(1.0-2.0*(Pi^2)*(f^2)*(t^2))*Exp[-(Pi^2)*(f^2)*(t^2)];
ListLinePlot[wavelet, Frame->True, PlotRange->Full, Filling->Axis, PlotStyle->Black,
PlotLabel->Style["Początkowy falowód", Black, 20],
LabelStyle->Directive[Black, Italic],
FillingStyle->{White, Black}, ImageSize->Large, InterpolationOrder->2]Początkowy impuls sejsmiczny
Ustalimy dwie granice na głębokościach 300 ms i 600 ms, a współczynniki odbicia będą losowymi liczbami.
rcExample=ConstantArray[0,1000];
rcExample[[300]]=RandomReal[{-1,0}];
rcExample[[600]]=RandomReal[{0,1}];
ListPlot[rcExample, Filling->0, Frame->True, Axes->False, PlotStyle->Black,
PlotLabel->Style["Współczynniki odbicia", Black, 20],
LabelStyle->Directive[Black, Italic]]Sekwencja współczynników odbicia
Obliczmy i przedstawmy ścieżkę sejsmiczną. Ponieważ współczynniki odbicia mają różne znaki, otrzymujemy na ścieżce sejsmicznej dwa odbicia zmiennego znaku.
traceExamle=ListConvolve[wavelet[[1;;;;1]], rcExample];
ListPlot[traceExamle,
PlotStyle->Black, Filling->0, Frame->True, Axes->False,
PlotLabel->Style["Ścieżka sejsmiczna", Black, 20],
LabelStyle->Directive[Black, Italic]]Modelowana ścieżka
W przypadku tego przykładu należy zaznaczyć, że w rzeczywistości głębokości warstw określane są w metrach, a obliczenia ścieżki sejsmicznej prowadzi się w obszarze czasowym. Lepiej byłoby ustalić głębokości w metrach i obliczyć czasy przybycia, znając prędkości w warstwach. W tym przypadku od razu ustaliłem warstwy na osi czasowej.
Jeśli mówimy o badaniach w terenie, w wyniku takich obserwacji rejestrowana jest ogromna ilość podobnych szeregów czasowych (tras sejsmicznych). Na przykład, podczas badania obszaru o długości 25 km i szerokości 15 km, gdzie w wyniku prac każda trasa charakteryzuje komórkę o wymiarach 25x25 metrów (taka komórka nazywana jest binem), ostateczny zbiór danych będzie zawierał 600000 tras. Przy kroku dyskretyzacji w czasie równym 1 ms, czasie zapisu wynoszącym 5 sekund, ostateczny plik danych przekroczy 11 GB, a objętość materiału „surowego” może wynieść setki gigabajtów.
Jak pracować z nimi w ?
Pakiet
Początek prac nad pakietem nastąpił na ścianie grupy VK rosyjskojęzycznego wsparcia. Dzięki odpowiedziom społeczności rozwiązanie zostało znalezione bardzo szybko. W rezultacie przekształciło się w poważny projekt. Odpowiedni nawet został wyróżniony przez moderatorów. W chwili obecnej pakiet obsługuje pracę z następującymi typami danych, które są aktywnie wykorzystywane w branży geologicznej:
- import danych kartograficznych w formacie ZMAP i IRAP
- import pomiarów w otworach w formacie LAS
- wejście i wyjście plików sejsmicznych w formacie
Aby zainstalować pakiet, należy postępować zgodnie z instrukcjami na stronie pobierania zbudowanego pakietu, tzn. wykonać następujący kod w dowolnym :
If[PacletInformation["GeologyIO"] === {}, PacletInstall[URLDownload[
"https://wolfr.am/FiQ5oFih",
FileNameJoin[{CreateDirectory[], "GeologyIO-0.2.2.paclet"}]
]]]Po tym pakiet zainstaluje się w domyślnym folderze, którego ścieżkę można uzyskać w następujący sposób:
FileNameJoin[{$UserBasePacletsDirectory, "Repository"}]Dla przykładu zaprezentujemy podstawowe możliwości pakietu. Wywołanie odbywa się tradycyjnie dla pakietów w języku Wolfram:
Get["GeologyIO`"]Pakiet jest rozwijany z użyciem . Umożliwia to dołączenie dokumentacji do głównych funkcji pakietu, która pod względem formatu przedstawienia nie różni się od dokumentacji samej Wolfram Mathematica, oraz dostarczenie pakietowi plików testowych do pierwszego zapoznania się.
Plik "Marmousi.segy" jest syntetycznym modelem geologicznego profilu, opracowanym przez francuski instytut naftowy. Używając tego modelu, programiści testują własne algorytmy symulacji pola falowego, przetwarzania danych, inwersji śladów sejsmicznych itd. Sam model Marmousi jest przechowywany w repozytorium, z którego został pobrany pakiet. Aby uzyskać plik, wykonujemy następujący kod:
If[Not[FileExistsQ["Marmousi.segy"]],
URLDownload["https://wolfr.am/FiQGh7rk", "Marmousi.segy"];]
marmousi = SEGYImport["Marmousi.segy"]Wynik importu – obiekt SEGYData
Format SEGY zakłada przechowywanie różnorodnych informacji o obserwacjach. Po pierwsze, to komentarze tekstowe. Zawierają one informacje o miejscu pracy, nazwy firm, które przeprowadziły pomiary itp. W naszym przypadku ten nagłówek jest wywoływany zapytaniem z kluczem TextHeader. Tutaj przedstawiono skrócony tekstowy nagłówek:
Short[marmousi["TextHeader"]]„Zbiór danych Marmousi został wygenerowany w Instytucie … minimalna prędkość 1500 m/s, a maksymalna 5500 m/s)”
Właściwą model geologiczny można wyświetlić odwołując się do sejsmicznych śladów za pomocą klucza „traces” (jedną z cech pakietu jest niezależność kluczy od wielkości liter):
ArrayPlot[Transpose[marmousi["traces"]], PlotTheme -> "Detailed"]Model Marmousi
Obecnie pakiet umożliwia również ładowanie danych częściami z dużych plików, co pozwala na przetwarzanie plików, których rozmiar może osiągać dziesiątki gigabajtów. W skład funkcji pakietu wchodzą również funkcje do eksportu danych do .segy i częściowego dopisywania na końcu pliku.
Warto również zwrócić uwagę na funkcjonalność pakietu przy pracy ze złożoną strukturą plików .segy. Umożliwia on nie tylko dostęp do poszczególnych śladów, nagłówków po kluczach i indeksach, ale także ich modyfikację z późniejszym zapisem do pliku. Wiele szczegółów technicznych realizacji GeologyIO wykracza poza zakres tego artykułu i prawdopodobnie zasługuje na osobne omówienie.
Znaczenie analizy spektralnej w sejsmoeksploracji
Możliwość importu danych sejsmicznych do Wolfram Mathematica pozwala na korzystanie z wbudowanych funkcji przetwarzania sygnałów dla danych eksperymentalnych. Ponieważ każda ścieżka sejsmiczna jest szeregami czasowymi, jednym z głównych narzędzi ich analizy jest analiza spektrogramu. Do przesłanek analizy składu częstotliwości danych sejsmicznych można zaliczyć, na przykład, następujące:
- Różne typy fal charakteryzują się odmiennym składem częstotliwości. Umożliwia to wyróżnienie użytecznych fal i stłumienie fal zakłócających.
- Cechy skał, takie jak porowatość i nasycenie, mogą wpływać na skład częstotliwości. Umożliwia to wyróżnianie skał o lepszych właściwościach.
- Warstwy o różnej grubości powodują anomalie w różnych zakresach częstotliwości.
Trzeci punkt jest kluczowy w kontekście tego artykułu. Poniżej znajduje się fragment kodu do obliczenia ścieżek sejsmicznych w przypadku warstwy o zmiennej grubości — model klina. Ten model jest tradycyjnie badany w sejsmice w celu analizy efektów interferencyjnych, gdy fale odbite od wielu warstw nakładają się na siebie.
nx=200;(* Liczba punktów siatki w kierunku X *)
ny=200;(* Liczba punktów siatki w kierunku Y *)
T=2;(* Czas całkowity propagacji *)
(* Prędkość i gęstość *)
modellv=Table[4000,{i,1,ny},{j,1,nx}];(* Prędkość fali P w m/s *)
rho=Table[2200,{i,1,ny},{j,1,nx}];(* Gęstość w g/cm³, używana jako stała gęstość *)
Table[modellv[[150-Round[i*0.5];;,i]]=4500;,{i,1,200}];
Table[modellv[[;;70,i]]=4500;,{i,1,200}];
(* Rysowanie modelu *)
MatrixPlot[modellv,PlotLabel->Style["Model warstwy",Black,20],
LabelStyle->Directive[Black,Italic]]Model klina
Prędkość fal wewnątrz klina wynosi 4500 m/s, poza klinem 4000 m/s, a gęstość przyjęto jako stałą 2200 g/cm³. Dla tego modelu obliczymy współczynniki odbicia i ścieżki sejsmiczne.
rc=Table[N[(modellv[[All,i]]-PadLeft[modellv[[All,i]],201,4000][[1;;200]])/(modellv[[All,i]]+PadLeft[modellv[[All,i]],201,4500][[1;;200]])],{i,1,200}];
traces=Table[ListConvolve[wavelet[[1;;;;1]],rc[[i]]],{i,1,200}];
starttrace=10;
endtrace=200;
steptrace=10;
trasenum=Range[starttrace,endtrace,steptrace];
traserenum=Range[Length@trasenum];
tracedist=0.5;
Rotate[Show[
Reverse[Table[
ListLinePlot[traces[[trasenum[[i]]]]*50+trasenum[[i]]*tracedist,Filling->{1->{trasenum[[i]]*tracedist,{RGBColor[0.97,0.93,0.68],Black}}},PlotStyle->Directive[Gray,Thin],PlotRange->Full,InterpolationOrder->2,Axes->False,Background->RGBColor[0.97,0.93,0.68]],
{i,1,Length@trasenum}]],ListLinePlot[Transpose[{ConstantArray[45,80],Range[80]}],PlotStyle->Red],PlotRange->All,Frame->True],270Degree]Ścieżki sejsmiczne dla modelu klina
Sekwencja tras sejsmicznych przedstawiona na tym rysunku nazywana jest przekrojem sejsmicznym. Jak można zauważyć, jej interpretacja może być przeprowadzona na poziomie intuicyjnym, ponieważ geometria odbitych fal jednoznacznie odpowiada wcześniej ustalonej modelowi. Analizując bardziej szczegółowo trasy, można zauważyć, że trasy od 1 do około 30 nie różnią się – odbicie od górnej granicy warstwy i od dolnej nie nakładają się na siebie. Od 31. trasy odbicia zaczynają interferować. I choć w modelu współczynniki odbicia nie zmieniają się w poziomie – trasy sejsmiczne zmieniają swoją intensywność przy zmianie grubości warstwy.
Rozważmy amplitudę odbicia od górnej granicy warstwy. Od 60. trasy intensywność odbicia zaczyna wzrastać, a na 70. trasie osiąga maksimum. Tak ujawnia się interferencja fal od górnej i dolnej granicy warstwy, prowadząc w niektórych przypadkach do znaczących anomalii w zapisie sejsmicznym.
ListLinePlot[GaussianFilter[Abs[traces[[All,46]]],3][[;;;;2]],
InterpolationOrder->2,Frame->True,PlotStyle->Black,
PlotLabel->Style["Amplitude of reflection",Black,20],
LabelStyle->Directive[Black,Italic],
PlotRange->All]Wykres amplitudy odbitej fali od górnej krawędzi klina
Logiczne jest, że gdy sygnał ma niższą częstotliwość, interferencja zaczyna ujawniać się przy większych grubościach warstwy, a w przypadku sygnału wysokoczęstotliwościowego interferencja występuje przy mniejszych grubościach. Następny fragment kodu tworzy sygnał z częstotliwościami 35 Hz, 55 Hz i 85 Hz.
waveletSet=Table[(1.0-2.0*(Pi^2)*(f^2)*(t^2))*Exp[-(Pi^2)*(f^2)*(t^2)],
{f,{35,55,85}}];
ListLinePlot[waveletSet,PlotRange->Full,PlotStyle->Black,Frame->True,
PlotLabel->Style["Set of wavelets",Black,20],
LabelStyle->Directive[Black,Italic],
ImageSize->Large,InterpolationOrder->2]Zestaw sygnałów źródłowych z częstotliwościami 35 Hz, 55 Hz, 85 Hz
Przeprowadzając obliczenia tras sejsmicznych i budując wykresy amplitud odbitej fali, możemy zobaczyć, że dla różnych częstotliwości anomalia jest obserwowana przy różnych grubościach warstwy.
tracesSet=Table[ListConvolve[waveletSet[[j]][[1;;;;1]],rc[[i]]],{j,1,3},{i,1,200}];
lowFreq=ListLinePlot[GaussianFilter[Abs[tracesSet[[1]][[All,46]]],3][[;;;;2]],InterpolationOrder->2,PlotStyle->Black,PlotRange->All];
medFreq=ListLinePlot[GaussianFilter[Abs[tracesSet[[2]][[All,46]]],3][[;;;;2]],InterpolationOrder->2,PlotStyle->Black,PlotRange->All];
highFreq=ListLinePlot[GaussianFilter[Abs[tracesSet[[3]][[All,46]]],3][[;;;;2]],InterpolationOrder->2,PlotStyle->Black,PlotRange->All];
Show[lowFreq,medFreq,highFreq,PlotRange->{{0,100},All},
PlotLabel->Style["Amplitudes of reflection",Black,20],
LabelStyle->Directive[Black,Italic],
Frame->True]Wykresy amplitud fali odbitej od górnej krawędzi klina dla różnych częstotliwości
Możliwość wnioskowania o grubości warstwy na podstawie wyników obserwacji sejsmicznych jest niezwykle przydatna, gdyż jednym z głównych zadań przy poszukiwaniu złóż ropy naftowej jest ocena najbardziej obiecujących miejsc do wiercenia (tj. tych obszarów, gdzie warstwa ma większą grubość). Oprócz tego w przekroju geologicznym mogą występować obiekty, które swoim genezą powodują nagłą zmianę grubości warstwy. To sprawia, że analiza spektralna jest skutecznym narzędziem do ich badania. W następnej części artykułu przyjrzymy się tym obiektom geologicznym bliżej.
Dane eksperymentalne. Gdzie zostały uzyskane i czego w nich szukać?
Materiał, który jest analizowany w artykule, został pozyskany na terenie Zachodniej Syberii. Region ten, jak zapewne wszyscy wiedzą, jest głównym regionem wydobycia ropy w naszym kraju. Aktywna eksploatacja złóż rozpoczęła się w tym regionie w latach 60. XX wieku. Główną metodą poszukiwania złóż ropy jest sejsmika. Interesująco jest przyjrzeć się zdjęciom satelitarnym tego obszaru. Przy małej skali można zauważyć ogromną liczbę bagien i jezior, powiększając mapę można zobaczyć wiertnicze pola wierceń, a przy maksymalnym powiększeniu rozróżnić też przecinki profili, po których przeprowadzono obserwacje sejsmiczne.
Zdjêcie satelitarne z map Yandex — okolice miasta Noyabrsk
Sieć wiertniczych pól na jednym z złóż
Rudy naftowe Zachodniej Syberii leżą w szerokim zakresie głębokości — od 1 km do 5 km. Główny wolumen skał zawierających ropę uformował się w okresie jury i kredy. Okres jury zapewne jest znany wielu ze względu na film o tej samej nazwie. znacząco różnił się od współczesności. W encyklopedii Britannica znajduje się seria paleokart, które charakteryzują każdą epokę geologiczną.
Obecny czas
Okres jurajski
Zauważ, że w czasie jury terytorium Zachodniej Syberii stanowiło nadmorską linię brzegową (ląd przecinany rzekami i płytkim morzem). Ponieważ klimat był komfortowy, można przypuszczać, że typowy krajobraz tego okresu wyglądał następująco:
Syberia okresu jurajskiego
Na tym obrazku istotne są nie tyle zwierzęta i ptaki, ile przedstawienie rzeki na tle. Rzeka to ten sam obiekt geologiczny, na którym wcześniej się skupiliśmy. Chodzi o to, że działalność rzek pozwala na gromadzenie dobrze sortowanych piaskowców, które następnie staną się zbiornikiem dla ropy naftowej. Te zbiorniki mogą mieć osobliwe, skomplikowane kształty (jak koryto rzeki) i zmienną grubość — przy brzegach grubość jest niewielka, a bliżej środka koryta lub w miejscach meandrów rośnie. Zatem, rzeki uformowane w okresie jurajskim obecnie znajdują się na głębokości około trzech kilometrów i są przedmiotem poszukiwań zbiorników ropy naftowej.
Dane eksperymentalne. Obrazowanie i wizualizacja
Zróbmy od razu zastrzeżenie dotyczące materiałów sejsmicznych przedstawionych w artykule — z uwagi na to, że objętość danych wykorzystanych do analizy jest znaczna — w treści artykułu zamieszczono jedynie fragment oryginalnego zbioru tras sejsmicznych. Pozwoli to wszystkim chętnym odtworzyć podane obliczenia.
Pracując z danymi sejsmicznymi, geofizyk zazwyczaj korzysta ze specjalistycznego oprogramowania (istnieje kilku liderów branży, których rozwiązania są szeroko stosowane, na przykład Petrel lub Paradigm), które umożliwia analizę różnych typów danych i posiada wygodny interfejs graficzny. Mimo że takie oprogramowanie jest bardzo wygodne, ma również swoje wady — na przykład wdrożenie nowoczesnych algorytmów w stabilnych wersjach zajmuje dużo czasu, a możliwości automatyzacji obliczeń są zazwyczaj ograniczone. W takiej sytuacji bardzo przydatne staje się wykorzystanie systemów komputerowej matematyki i języków programowania wysokiego poziomu, które pozwalają korzystać z szerokiej bazy algorytmów i jednocześnie przejmują wiele rutynowych zadań. Takie podejście zbudowane jest na bazie pracy z danymi sejsmicznymi w Wolfram Mathematica. Nie ma sensu pisać bogatej funkcjonalności interaktywnej pracy z danymi - ważniejsze jest zapewnienie możliwości ładowania z powszechnie akceptowanego formatu, zastosowanie do nich pożądanych algorytmów i ponowne wyeksportowanie w zewnętrznym formacie.
Zgodnie z przedstawionym schematem, załadujemy oryginalne dane sejsmiczne i wyświetlimy je w :
Get["GeologyIO`"]
seismic3DZipPath = "seismic3D.zip";
seismic3DSEGYPath = "seismic3D.sgy";
If[FileExistsQ[seismic3DZipPath], DeleteFile[seismic3DZipPath]];
If[FileExistsQ[seismic3DSEGYPath], DeleteFile[seismic3DSEGYPath]];
URLDownload["https://wolfr.am/FiQIuZuH", seismic3DZipPath];
ExtractArchive[seismic3DZipPath];
seismic3DSEGY = SEGYImport[seismic3DSEGYPath]Załadowane i zaimportowane w ten sposób dane to trasy zarejestrowane na obszarze o rozmiarze 10 na 5 kilometrów. W przypadku, gdy dane zostały pozyskane w metodzie trójwymiarowego sejsmicznego badania (rejestracja fal odbywa się nie wzdłuż poszczególnych profili geofizycznych, a na całym obszarze jednocześnie), możliwe jest uzyskanie sześcianów danych sejsmicznych. Są to obiekty trójwymiarowe, których przekroje pionowe i poziome pozwalają szczegółowo badać środowisko geologiczne. W omawianym przykładzie mamy do czynienia z danymi trójwymiarowymi. Niektóre informacje możemy uzyskać z nagłówka tekstowego, na przykład tak
StringPartition[seismic3DSEGY["textheader"], 80] // TableFormC 1 TO JEST PLIKI DEMO DLA TESTU PAKIETU GEOLOGYIO
C 2
C 3
C 4
C 5 DATA UŻYTKOWNIKA: UŻYTKOWNIK WOLFRAM
C 6 NAZWA BADANIA: GDZIEŚ W SYBERII
C 7 TYP PLIKU 3D SEISMIC VOLUME
C 8
C 9
C10 Z ZAKRES: PIERWSZE 2200M OSTATNIE 2400M
Ten zbiór danych wystarczy, aby zaprezentować podstawowe etapy analizy danych. Ślady w pliku zapisane są sekwencyjnie, a każdy z nich wygląda mniej więcej jak na następnym rysunku - to rozkład amplitud odbitych fal wzdłuż osi pionowej (osi głębokości).
ListLinePlot[seismic3DSEGY["traces"][[100]], InterpolationOrder -> 2, PlotStyle -> Black, PlotLabel -> Style["Seismic trace", Black, 20], LabelStyle -> Directive[Black, Italic], PlotRange -> All, Frame -> True, ImageSize -> 1200, AspectRatio -> 1/5]Jeden z śladów sejsmicznego przekroju
Znać liczbę śladów w każdym kierunku badanej sekcji, można utworzyć trójwymiarową macierz danych i wyświetlić ją za pomocą funkcji Image3D[]
traces=seismic3DSEGY["traces"];
startIL=1050;EndIL=2000;stepIL=2; (* współrzędna X początku i końca pomiaru oraz krok śladów *)
startXL=1165;EndXL=1615;stepXL=2; (* współrzędna Y początku i końca pomiaru oraz krok śladów *)
numIL=(EndIL-startIL)/stepIL+1; (* liczba śladów na osi X *)
numXL=(EndXL-startXL)/stepIL+1; (* liczba śladów na osi Y *)
Image3D[ArrayReshape[Abs[traces/Max[Abs[traces[[All,1;;;;4]]]]],{numIL,numXL,101}],ViewPoint->{-1, 0, 0},Background->RGBColor[0,0,0]]Trójwymiarowy obraz sześcianu danych sejsmicznych. (Oś pionowa — głębokość)
W przypadku gdy obiekty geologiczne, które są interesujące, tworzą intensywne anomalie sejsmiczne, można użyć narzędzi wizualizacyjnych z przezroczystością. "Niezważane" sekcje rejestracji można uczynić niewidocznymi, pozostawiając widoczne tylko anomalie. W Wolfram Mathematica można to zrobić za pomocą i .
data = ArrayReshape[Abs[traces/Max[Abs[traces[[All,1;;;;4]]]]],{numIL,numXL,101}];
Graphics3D[{Opacity[0.1], Raster3D[data, ColorFunction->"RainbowOpacity"]}, Boxed->False, SphericalRegion->True, ImageSize->840, Background->None]Obraz sześcianu danych sejsmicznych z użyciem funkcji Opacity[] i Raster3D[]
Jak w przypadku syntetycznego przykładu, na przekrojach oryginalnego sześcianu można wyróżnić pewne granice geologiczne (warstwy) z zmiennym ukształtowaniem.
Głównym narzędziem analizy spektralnej jest transformata Fouriera. Dzięki niej można ocenić amplitudowo-czestotliwościowy spektrum każdej ścieżki lub grupy ścieżek. Jednak po przetworzeniu danych na obszar częstotliwości traci się informacje o tym, w jakich momentach (czytaj, na jakich głębokościach) zmienia się częstotliwość. Aby móc lokalizować zmiany sygnału na osi czasowej (głębokości) stosuje się okienkową transformatę Fouriera oraz rozkład falkowy. W tym artykule wykorzystano rozkład falkowy. Technologia analizy falkowej zaczęła być powszechnie stosowana w sejsmologii w latach 90. XX wieku. Jej przewagą nad okienkową transformacją Fouriera jest lepsza rozdzielczość czasowa.
Za pomocą następującego fragmentu kodu można przeprowadzić rozkład na poszczególne składniki jednej z sejsmicznych ścieżek:
cwd=ContinuousWaveletTransform[seismicSection["traces"][[100]]]
Show[
ListLinePlot[Re[cwd[[1]]],PlotRange->All],
ListLinePlot[seismicSection["traces"][[100]],
PlotStyle->Black,PlotRange->All],ImageSize->{1500,500},AspectRatio->Full,
PlotLabel->Style["Rozkład falkowy",Black,32],
LabelStyle->Directive[Black,Italic],
PlotRange->All,
Frame->True]Rozkład ścieżki na składniki
Aby ocenić, jak energia odbicia jest rozłożona w różnych momentach przybycia fal, stosuje się skalogramy (analog spektrogramu). Zazwyczaj w praktyce nie ma potrzeby analizować wszystkich składników. Zwykle wybiera się składnik niskoczęstotliwościowy, średnioczęstotliwościowy i wysokoczęstotliwościowy.
freq=(500/(#*contWD["Wavelet"]["FourierFactor"]))&/@(Thread[{Range[contWD["Octaves"]],1}]/.contWD["Scales"])//Round;
ticks=Transpose[{Range[Length[freq]],freq}];
WaveletScalogram[contWD,Frame->True,FrameTicks->{{ticks,Automatic},Automatic},FrameTicksStyle->Directive[Orange,12],
FrameLabel->{"Czas","Częstotliwość(Hz)"},LabelStyle->Directive[Black,Bold,14],
ColorFunction->"RustTones",ImageSize->Large]Skalogram. Wynik funkcji
W Wolfram Language do transformacji falkowej używana jest funkcja . Zastosowanie tej funkcji do całego zbioru ścieżek realizuje się za pomocą funkcji . Należy zauważyć, że jedną z mocnych stron Wolfram Mathematica jest możliwość wykorzystania równoległego przetwarzania. . W podanym przykładzie nie ma potrzeby stosowania równoległego przetwarzania — rozmiar danych nie jest duży, ale w przypadku pracy z eksperymentalnymi zbiorami danych zawierającymi setki tysięcy ścieżek jest to koniecznością.
tracesCWD=Table[Map[Hilbert[#,0]&,Re[ContinuousWaveletTransform[traces[[i]]][[1]]][[{13,15,18}]]],{i,1,Length@traces}]; Po zastosowaniu funkcji Pojawiają się nowe zbiory danych odpowiadające wybranym częstotliwościom. W podanym powyżej przykładzie są to częstotliwości: 38 Hz, 33 Hz, 27 Hz. Wybór częstotliwości najczęściej dokonuje się na podstawie testów — uzyskuje się efektywne mapy dla różnych kombinacji częstotliwości i wybiera się najbardziej informatywną z perspektywy geologa.
Jeśli wyniki trzeba udostępnić kolegom lub przedstawić klientowi, można skorzystać z funkcji SEGYExport[] pakietu GeologyIO.
outputdata=seismic3DSEGY;
outputdata["traces",1;;-1]=tracesCWD[[All,3]];
outputdata["textheader"]="Wynik dekompozycji fal";
outputdata["binaryheader","NumberDataTraces"]=Length[tracesCWD[[All,3]]];
SEGYExport["D:result.segy",outputdata];Mając do dyspozycji trzy takie sześciany (niskoczęstotliwościową, średnioczęstotliwościową i wysokoczęstotliwościową komponentę), zazwyczaj stosuje się mieszanie RGB do wspólnej wizualizacji danych. Każdej z komponent przypisuje się swój kolor - czerwony, zielony, niebieski. W Wolfram Mathematica można to zrobić za pomocą funkcji .
W rezultacie uzyskuje się obrazy, na podstawie których można przeprowadzać interpretację geologiczną. Meandry, które są widoczne na przekroju, pozwalają określić paleorusła, które z większym prawdopodobieństwem mogą być zbiornikami i zawierać złoża ropy naftowej. Poszukiwanie i analiza współczesnych analogów takiego systemu rzekomego pozwala określić najbardziej obiecujące części meandrów. Właściwe rusła charakteryzują się grubymi warstwami dobrze sortowanego piaskowca i są dobrym zbiornikiem dla ropy. Obszary poza „sznurkowymi” anomaliami są podobne do współczesnych osadów zalewowych. Osady zalewowe głównie składają się z glin i wiercenie w tych strefach będzie nieefektywne.
RGB przekrój sześcianu danych. W centrum (nieco na lewo od centrum) można dostrzec meandrującą rzekę.
RGB przekrój sześcianu danych. Po lewej stronie można dostrzec meandrującą rzekę.
W niektórych przypadkach jakość danych sejsmicznych pozwala uzyskać znacznie wyraźniejsze obrazy. Zależy to od metodologii prac polowych, sprzętu oraz stosowanych algorytmów tłumienia hałasu. W takich przypadkach widoczne są nie tylko fragmenty systemów rzek, ale i całe rozległe paleorzeki.
Mieszanie RGB trzech komponent sześcianu danych sejsmicznych (przekrój poziomy). Głębokość wynosi około 2 km.
Obraz satelitarny rzeki Wołgi w rejonie Saratowa.
Podsumowanie
W Wolfram Mathematica można analizować dane sejsmiczne i rozwiązywać zastosowania związane z poszukiwaniem surowców naturalnych, a pakiet GeologyIO sprawia, że ten proces jest wygodniejszy. Struktura danych sejsmicznych jest taka, że wykorzystanie wbudowanych metod przyspieszania obliczeń (, ,…) jest bardzo efektywne i pozwala na przetwarzanie dużych zbiorów danych. W znacznym stopniu sprzyjają temu cechy przechowywania danych pakietu GeologyIO. Co ciekawe, pakiet może być stosowany nie tylko w dziedzinie zastosowanej sejsmologii. Praktycznie te same typy danych są wykorzystywane w georadarach i sejsmologii. Jeśli masz pomysły, jak poprawić wyniki, jakie algorytmy analizy sygnału z arsenalu Wolfram Mathematica można zastosować do takich danych lub masz krytyczne uwagi - zostaw komentarze.
Źródło: habr.com
