Dziękuję autorowi bloga za jego raport
Wprowadzenie
Niniejsza notatka została napisana w następstwie konferencji i zawiera streszczenie raportu, który przekazałem. Wydarzenie odbyło się w czerwcu w Petersburgu. Biorąc pod uwagę, że pracuję niedaleko miejsca konferencji, nie mogłem nie wziąć udziału w tym wydarzeniu. W 2016 i 2017 roku wysłuchałem relacji z konferencji, w tym roku wygłosiłem prezentację. Po pierwsze pojawił się ciekawy (wydaje mi się) temat, nad którym się rozwijamy , a po drugie, po długim badaniu ustawodawstwa Federacji Rosyjskiej dotyczącego polityki sankcji, w przedsiębiorstwie, w którym pracuję, pojawiły się aż dwie licencje .
Zanim przejdę do tematu mojego wystąpienia, chciałbym zwrócić uwagę na dobrą organizację wydarzenia. Na stronie odwiedzającej konferencję znajduje się wizerunek katedry kazańskiej. Katedra jest jedną z głównych atrakcji Petersburga i jest bardzo dobrze widoczna z sali, w której odbywała się konferencja.
Przy wejściu do Państwowego Uniwersytetu Ekonomicznego w Petersburgu na uczestników czekały asystentki spośród studentów, które nie pozwoliły im się zgubić. Podczas rejestracji rozdawane były drobne upominki (zabawka – świecący kolec, długopis, naklejki z symbolami Wolframa). W programie konferencji przewidziano także lunch i przerwy kawowe. O pysznej kawie i ciastach wspominałam już na ścianie grupy – szefowie kuchni są świetni. Już na wstępie pragnę podkreślić, że samo wydarzenie, jego format i lokalizacja budzą już pozytywne emocje.
Raport, który przygotowaliśmy wspólnie z Kirillem Biełowem, nosi tytuł „Wykorzystanie programu Wolfram Mathematica do rozwiązywania problemów w geofizyce stosowanej. Analiza widmowa danych sejsmicznych lub „tam, gdzie płynęły starożytne rzeki”. Treść raportu obejmuje dwie części: po pierwsze, wykorzystanie algorytmów dostępnych w do analizy danych geofizycznych, a po drugie, jak umieścić dane geofizyczne w Wolfram Mathematica.
Eksploracja sejsmiczna
Najpierw musisz odbyć krótką wycieczkę do geofizyki. Geofizyka to nauka zajmująca się badaniem właściwości fizycznych skał. Cóż, ponieważ skały mają różne właściwości: elektryczne, magnetyczne, sprężyste, istnieją odpowiednie metody geofizyki: poszukiwania elektryczne, poszukiwania magnetyczne, poszukiwania sejsmiczne... W kontekście tego artykułu omówimy bardziej szczegółowo jedynie poszukiwania sejsmiczne. Główną metodą poszukiwania ropy i gazu są badania sejsmiczne. Metoda polega na wzbudzeniu drgań sprężystych i późniejszej rejestracji reakcji ze skał tworzących obszar badań. Wibracje wzbudzane są na lądzie (za pomocą dynamitu lub niewybuchowych źródeł drgań sprężystych) lub na morzu (za pomocą wiatrówek). Drgania sprężyste rozchodzą się w górotworze, ulegają załamaniu i odbiciu na granicach warstw o różnych właściwościach. Odbite fale wracają na powierzchnię i są rejestrowane przez geofony na lądzie (najczęściej urządzenia elektrodynamiczne oparte na ruchu magnesu zawieszonego w cewce) lub hydrofony w morzu (bazujące na efekcie piezoelektrycznym). Do czasu nadejścia fal można ocenić głębokość warstw geologicznych.
Sprzęt do holowania statków sejsmicznych
Wiatrówka wzbudza drgania sprężyste
Fale przechodzą przez górotwór i są rejestrowane przez hydrofony
Statek badawczy do badań geofizycznych „Iwan Gubkin” na molo w pobliżu mostu Błagowieszczeńskiego w Petersburgu
Model sygnału sejsmicznego
Skały mają różne właściwości fizyczne. Przy badaniach sejsmicznych istotne są przede wszystkim właściwości sprężyste – prędkość propagacji drgań sprężystych oraz gęstość. Jeśli dwie warstwy mają takie same lub podobne właściwości, to fala „nie zauważy” granicy między nimi. Jeśli prędkości fal w warstwach będą się różnić, wówczas na granicy warstw nastąpi odbicie. Im większa różnica właściwości, tym intensywniejsze odbicie. Jego intensywność będzie określona przez współczynnik odbicia (rc):
gdzie ρ to gęstość skał, ν to prędkość fali, 1 i 2 oznaczają warstwę górną i dolną.
Jednym z najprostszych i najczęściej stosowanych modeli sygnału sejsmicznego jest model splotowy, w którym zarejestrowany ślad sejsmiczny jest reprezentowany w wyniku splotu ciągu współczynników odbicia z impulsem sondującym:
gdzie s(t) — ślad sejsmiczny, tj. wszystko co zostało zarejestrowane przez hydrofon lub geofon w ustalonym czasie nagrywania, w(t) - sygnał generowany przez wiatrówkę, n(t) - losowy hałas.
Jako przykład obliczmy syntetyczny ślad sejsmiczny. Jako sygnał początkowy wykorzystamy impuls Rickera, szeroko stosowany w badaniach sejsmicznych.
length=0.050; (*Signal lenght*)
dt=0.001;(*Sample rate of signal*)
t=Range[-length/2,(length)/2,dt];(*Signal time*)
f=35;(*Central frequency*)
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["Initial wavelet",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ą liczbami losowymi
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["Reflection Coefficients",Black,20],
LabelStyle->Directive[Black,Italic]]Sekwencja współczynników odbicia
Obliczmy i wyświetlmy ślad sejsmiczny. Ponieważ współczynniki odbicia mają różne znaki, otrzymujemy dwa naprzemienne odbicia na śladzie sejsmicznym.
traceExamle=ListConvolve[wavelet[[1;;;;1]],rcExample];
ListPlot[traceExamle,
PlotStyle->Black,Filling->0,Frame->True,Axes->False,
PlotLabel->Style["Seismic trace",Black,20],
LabelStyle->Directive[Black,Italic]]Symulowany tor
Dla tego przykładu konieczne jest dokonanie rezerwacji – w rzeczywistości głębokość warstw określa się oczywiście w metrach, a obliczenie śladu sejsmicznego następuje dla dziedziny czasu. Bardziej poprawne byłoby podawanie głębokości w metrach i obliczanie czasów przybycia, znając prędkości w warstwach. W tym przypadku od razu ustawiam warstwy na osi czasu.
Jeśli mówimy o badaniach terenowych, to w wyniku takich obserwacji rejestrowana jest ogromna liczba podobnych szeregów czasowych (śladów sejsmicznych). Przykładowo badając teren o długości 25 km i szerokości 15 km, gdzie w wyniku pracy każdy ślad charakteryzuje komórkę o wymiarach 25x25 metrów (taka komórka nazywa się koszem), ostateczna tablica danych będzie zawierać 600000 1 śladów. Przy czasie próbkowania 5 ms i czasie nagrywania 11 sekund ostateczny plik danych będzie miał ponad XNUMX GB, a objętość oryginalnego „surowego” materiału może sięgać setek gigabajtów.
Jak z nimi pracować ?
Pakiet
Rozpoczęto prace nad pakietem na ścianie VK rosyjskojęzycznej grupy wsparcia. Dzięki odpowiedziom społeczności bardzo szybko znaleziono rozwiązanie. W rezultacie przerodziło się to w poważny rozwój. Odpowiedni Zostało to nawet oznaczone przez moderatorów. Aktualnie pakiet wspiera pracę z następującymi typami danych, które są aktywnie wykorzystywane w branży geologicznej:
- import danych mapowych w formatach ZMAP i IRAP
- import pomiarów w studniach w formacie LAS
- wejście i wyjście w formacie plików sejsmicznych
Aby zainstalować pakiet, należy postępować zgodnie z instrukcjami na stronie pobierania zmontowanego pakietu, tj. wykonaj następujący kod w any :
If[PacletInformation["GeologyIO"] === {}, PacletInstall[URLDownload[
"https://wolfr.am/FiQ5oFih",
FileNameJoin[{CreateDirectory[], "GeologyIO-0.2.2.paclet"}]
]]]Po czym pakiet zostanie zainstalowany w folderze domyślnym, do którego ścieżkę można uzyskać w następujący sposób:
FileNameJoin[{$UserBasePacletsDirectory, "Repository"}]Jako przykład zademonstrujemy główne możliwości pakietu. Wywołanie odbywa się tradycyjnie dla pakietów w języku Wolfram:
Get["GeologyIO`"]Pakiet jest rozwijany przy użyciu . Pozwala to na dołączenie do głównej funkcjonalności pakietu dokumentacji, która pod względem formatu prezentacji nie odbiega od dokumentacji samego Wolfram Mathematica oraz dostarczenie pakietu plików testowych przy pierwszym zapoznaniu się.
Takim plikiem jest w szczególności plik „Marmousi.segy” - jest to syntetyczny model przekroju geologicznego opracowany przez Francuski Instytut Naftowy. Korzystając z tego modelu, programiści testują własne algorytmy modelowania pola falowego, przetwarzania danych, inwersji śladów sejsmicznych itp. Sam model Marmousi przechowywany jest w repozytorium, z którego został pobrany sam pakiet. Aby pobrać plik, uruchom 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 polega na przechowywaniu różnych informacji o obserwacjach. Po pierwsze, są to komentarze tekstowe. Obejmuje to informacje o lokalizacji prac, nazwach firm, które wykonywały pomiary itp. W naszym przypadku nagłówek ten wywoływany jest poprzez żądanie z kluczem TextHeader. Oto skrócony nagłówek tekstu:
Short[marmousi["TextHeader"]]„Zbiór danych Marmousi został wygenerowany w Instytucie... prędkość minimalna 1500 m/s i maksymalna 5500 m/s)”
Możesz wyświetlić rzeczywisty model geologiczny, uzyskując dostęp do śladów sejsmicznych za pomocą klucza „ślady” (jedną z cech pakietu jest to, że w kluczach nie jest rozróżniana wielkość liter):
ArrayPlot[Transpose[marmousi["traces"]], PlotTheme -> "Detailed"]Modelka Marmusi
Obecnie pakiet umożliwia także ładowanie danych w częściach z dużych plików, dzięki czemu możliwa jest obróbka plików, których rozmiar może sięgać kilkudziesięciu gigabajtów. Funkcje pakietu obejmują również funkcje eksportu danych do .segy i częściowego dopisywania na końcu pliku.
Osobno warto zwrócić uwagę na funkcjonalność pakietu podczas pracy ze złożoną strukturą plików .segy. Ponieważ umożliwia nie tylko dostęp do poszczególnych śladów i nagłówków za pomocą kluczy i indeksów, ale także ich zmianę, a następnie zapisanie ich do pliku. Wiele szczegółów technicznych implementacji GeologyIO wykracza poza zakres tego artykułu i prawdopodobnie zasługuje na osobny opis.
Znaczenie analizy spektralnej w badaniach sejsmicznych
Możliwość importowania danych sejsmicznych do Wolfram Mathematica pozwala na wykorzystanie wbudowanej funkcji przetwarzania sygnałów dla danych eksperymentalnych. Ponieważ każdy ślad sejsmiczny reprezentuje szereg czasowy, jednym z głównych narzędzi do ich badania jest analiza spektralna. Wśród warunków wstępnych analizy składu częstotliwości danych sejsmicznych możemy wymienić na przykład:
- Różne rodzaje fal charakteryzują się różnym składem częstotliwości. Pozwala to wyróżnić fale użyteczne i stłumić fale zakłócające.
- Właściwości skał, takie jak porowatość i nasycenie, mogą wpływać na skład częstotliwości. Dzięki temu możliwe jest zidentyfikowanie skał o najlepszych właściwościach.
- Warstwy o różnej grubości powodują anomalie w różnych zakresach częstotliwości.
Trzeci punkt jest najważniejszy w kontekście tego artykułu. Poniżej fragment kodu do obliczania śladów sejsmicznych w przypadku warstwy o różnej grubości – model klinowy. Model ten jest tradycyjnie badany podczas badań sejsmicznych w celu analizy efektów interferencji, gdy fale odbite od wielu warstw nakładają się na siebie.
nx=200;(* Number of grid points in X direction*)
ny=200;(* Number of grid points in Y direction*)
T=2;(*Total propagation time*)
(*Velocity and density*)
modellv=Table[4000,{i,1,ny},{j,1,nx}];(* P-wave velocity in m/s*)
rho=Table[2200,{i,1,ny},{j,1,nx}];(* Density in g/cm^3, used constant density*)
Table[modellv[[150-Round[i*0.5];;,i]]=4500;,{i,1,200}];
Table[modellv[[;;70,i]]=4500;,{i,1,200}];
(*Plotting model*)
MatrixPlot[modellv,PlotLabel->Style["Model of layer",Black,20],
LabelStyle->Directive[Black,Italic]]Model formacji pinch-out
Prędkość fali wewnątrz klina wynosi 4500 m/s, na zewnątrz klina 4000 m/s, a gęstość przyjmuje się jako stałą 2200 g/cmXNUMX. Dla takiego modelu obliczamy współczynniki odbicia i ślady 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]Ślady sejsmiczne dla modelu klinowego
Sekwencja śladów sejsmicznych pokazana na tym rysunku nazywana jest sekcją sejsmiczną. Jak widać, jego interpretacji można dokonać także na poziomie intuicyjnym, gdyż geometria odbitych fal wyraźnie odpowiada modelowi określonemu wcześniej. Jeśli dokładniej przeanalizujesz ślady, zauważysz, że ślady od 1 do około 30 nie różnią się od siebie – odbicie od stropu formacji i od dna nie nakładają się na siebie. Począwszy od 31. śladu odbicia zaczynają się zakłócać. I choć w modelu współczynniki odbicia nie zmieniają się w poziomie, to ślady sejsmiczne zmieniają swoją intensywność wraz ze zmianą miąższości formacji.
Rozważmy amplitudę odbicia od górnej granicy formacji. Począwszy od 60. trasy intensywność odbicia zaczyna rosnąć, a na 70. trasie osiąga maksimum. W ten sposób objawia się interferencja fal ze stropu i dna warstw, prowadząc w niektórych przypadkach do znacznych 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 fali odbitej od górnej krawędzi klina
Logiczne jest, że gdy sygnał ma niższą częstotliwość, zakłócenia zaczynają pojawiać się przy dużych grubościach formacji, a w przypadku sygnału o wysokiej częstotliwości zakłócenia pojawiają się przy mniejszych grubościach. Poniższy fragment kodu tworzy sygnał o częstotliwościach 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 o częstotliwościach 35 Hz, 55 Hz, 85 Hz
Obliczając ślady sejsmiczne i sporządzając wykresy amplitud fal odbitych, możemy zobaczyć, że dla różnych częstotliwości obserwuje się anomalię przy różnych grubościach formacji.
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
Umiejętność wyciągania wniosków na temat miąższości formacji na podstawie wyników obserwacji sejsmicznych jest niezwykle przydatna, ponieważ jednym z głównych zadań poszukiwań ropy naftowej jest ocena najbardziej perspektywicznych miejsc do ułożenia odwiertu (czyli obszarów, w których formacja jest grubszy). Ponadto w przekroju geologicznym mogą występować obiekty, których geneza powoduje gwałtowną zmianę miąższości formacji. To sprawia, że analiza spektralna jest skutecznym narzędziem do ich badania. W dalszej części artykułu przyjrzymy się bliżej takim obiektom geologicznym.
Dane eksperymentalne. Skąd je kupiłeś i czego w nich szukać?
Materiały analizowane w artykule pozyskano z zachodniej Syberii. Region, jak zapewne wiedzą wszyscy bez wyjątku, jest głównym regionem wydobywczym ropy naftowej w naszym kraju. Aktywny rozwój złóż rozpoczął się w regionie w latach 60-tych ubiegłego wieku. Główną metodą poszukiwania złóż ropy naftowej są badania sejsmiczne. Interesujące jest spojrzenie na zdjęcia satelitarne tego terytorium. W małej skali można zauważyć ogromną liczbę bagien i jezior, powiększając mapę można zobaczyć miejsca wierceń studni klastrowych, a powiększając mapę do granic możliwości można także rozróżnić prześwity profili, wzdłuż których przebiegają ruchy sejsmiczne. przeprowadzono obserwacje.
Zdjęcie satelitarne map Yandex - obszar miasta Noyabrsk
Sieć studni na jednym z pól
Skały roponośne zachodniej Syberii występują na różnych głębokościach – od 1 km do 5 km. Główna część skał zawierających ropę powstała w okresie jurajskim i kredowym. Okres jurajski jest zapewne znany wielu z filmu o tym samym tytule. różnił się znacząco od współczesnego. Encyklopedia Britannica zawiera serię paleomap charakteryzujących każdą epokę helogiczną.
Czas teraźniejszy
Okres jurajski
Należy pamiętać, że w czasach jurajskich terytorium zachodniej Syberii było wybrzeżem morskim (kraj przecinany rzekami i płytkim morzem). Ponieważ klimat był komfortowy, możemy założyć, że typowy krajobraz tamtych czasów wyglądał następująco:
Jurajska Syberia
Na tym zdjęciu ważne są dla nas nie tyle zwierzęta i ptaki, ile obraz rzeki w tle. Rzeka to ten sam obiekt geologiczny, przy którym zatrzymaliśmy się wcześniej. Faktem jest, że działalność rzek pozwala na gromadzenie się dobrze wysortowanych piaskowców, które następnie staną się zbiornikiem ropy. Zbiorniki te mogą mieć dziwaczny, skomplikowany kształt (jak koryto rzeki) i mieć zmienną miąższość - w pobliżu brzegów miąższość jest niewielka, ale bliżej środka koryta lub w obszarach meandrowych wzrasta. Tak więc rzeki utworzone w okresie jurajskim znajdują się obecnie na głębokości około trzech kilometrów i są przedmiotem poszukiwań złóż ropy.
Dane eksperymentalne. Przetwarzanie i wizualizacja
Należy od razu zgłosić zastrzeżenie dotyczące przedstawionych w artykule materiałów sejsmicznych – z uwagi na fakt, że ilość danych wykorzystanych do analizy jest znacząca – w tekście artykułu uwzględniono jedynie fragment oryginalnego zestawu śladów sejsmicznych. Dzięki temu każdy może odtworzyć powyższe obliczenia.
Podczas pracy z danymi sejsmicznymi geofizyk zwykle korzysta ze specjalistycznego oprogramowania (jest kilku liderów branży, z których opracowań aktywnie korzysta, na przykład Petrel czy Paradigm), które pozwala analizować różne rodzaje danych i posiada wygodny interfejs graficzny. Mimo całej wygody tego typu oprogramowanie ma też swoje wady – np. implementacja nowoczesnych algorytmów w stabilnych wersjach zajmuje dużo czasu, a możliwości automatyzacji obliczeń są zwykle ograniczone. W takiej sytuacji bardzo wygodne staje się korzystanie z systemów matematyki komputerowej i języków programowania wysokiego poziomu, które pozwalają na wykorzystanie szerokiej bazy algorytmicznej, a jednocześnie wymagają dużej rutyny. Jest to zasada stosowana podczas pracy z danymi sejsmicznymi w programie Wolfram Mathematica. Niewłaściwe jest pisanie bogatej funkcjonalności do interaktywnej pracy z danymi - ważniejsze jest zapewnienie ładowania z ogólnie przyjętego formatu, nałożenie na nie pożądanych algorytmów i wgranie z powrotem do formatu zewnętrznego.
Zgodnie z zaproponowanym schematem załadujemy oryginalne dane sejsmiczne i wyświetlimy je :
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]Pobrane i zaimportowane w ten sposób dane to trasy zarejestrowane na obszarze o wymiarach 10 na 5 kilometrów. Jeżeli dane pozyskiwane są metodą trójwymiarowych badań sejsmicznych (rejestracja fal nie odbywa się wzdłuż poszczególnych profili geofizycznych, ale jednocześnie na całym obszarze), możliwe staje się uzyskanie kostek danych sejsmicznych. Są to trójwymiarowe obiekty, których przekroje pionowe i poziome pozwalają na szczegółowe badanie środowiska geologicznego. W rozważanym przykładzie mamy do czynienia z danymi trójwymiarowymi. Możemy uzyskać pewne informacje z nagłówka tekstu, na przykład poniżej
StringPartition[seismic3DSEGY["textheader"], 80] // TableFormC 1 TO JEST PLIK DEMONICZNY DO TESTU PAKIETU GEOLOGYIO
C 2
C 3
C 4
C 5 DATA NAZWA UŻYTKOWNIKA: UŻYTKOWNIK WOLFRAM
C 6 NAZWA BADANIA: GDZIEŚ NA SYBERII
C 7 TYP PLIKU OBJĘTOŚĆ SEJSMICZNA 3D
C 8
C 9
ZASIĘG C10 Z: PIERWSZE 2200M OSTATNIE 2400M
Ten zestaw danych wystarczy nam do zademonstrowania głównych etapów analizy danych. Ślady w pliku zapisywane są sekwencyjnie i każdy z nich wygląda mniej więcej tak, jak na poniższym rysunku - jest to rozkład amplitud fal odbitych 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 ze śladów sekcji sejsmicznej
Wiedząc, ile śladów znajduje się w każdym kierunku badanego obszaru, możesz wygenerować trójwymiarową tablicę danych i wyświetlić ją za pomocą funkcji Image3D[]
traces=seismic3DSEGY["traces"];
startIL=1050;EndIL=2000;stepIL=2; (*координата Х начала и конца съёмки и шаг трасс*)
startXL=1165;EndXL=1615;stepXL=2; (*координата Y начала и конца съёмки и шаг трасс*)
numIL=(EndIL-startIL)/stepIL+1; (*количество трасс по оис Х*)
numXL=(EndXL-startXL)/stepIL+1; (*количество трасс по оис Y*)
Image3D[ArrayReshape[Abs[traces/Max[Abs[traces[[All,1;;;;4]]]]],{numIL,numXL,101}],ViewPoint->{-1, 0, 0},Background->RGBColor[0,0,0]]Obraz XNUMXD kostki danych sejsmicznych (oś pionowa – głębokość)
Jeśli interesujące cechy geologiczne powodują intensywne anomalie sejsmiczne, można zastosować narzędzia do wizualizacji z przezroczystością. „Nieistotne” obszary nagrania można uczynić niewidocznymi, pozostawiając widoczne jedynie anomalie. W Wolfram Mathematica można to zrobić za pomocą и .
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 kostki danych sejsmicznych przy użyciu funkcji Opacity[] i Raster3D[].
Podobnie jak w przykładzie syntetycznym, na przekrojach pierwotnej sześcianu można zidentyfikować pewne granice (warstwy) geologiczne o zmiennej rzeźbie.
Głównym narzędziem analizy widmowej jest transformata Fouriera. Za jego pomocą można ocenić widmo amplitudowo-częstotliwościowe każdego śladu lub grupy śladów. Jednak po przesłaniu danych do dziedziny częstotliwości tracona jest informacja o tym, w jakich momentach (czytaj na jakich głębokościach) zmienia się częstotliwość. Aby móc zlokalizować zmiany sygnału na osi czasu (głębokości), wykorzystuje się okienkową transformatę Fouriera i dekompozycję falkową. W tym artykule wykorzystano rozkład falkowy. Technologia analizy falkowej zaczęła być aktywnie wykorzystywana w badaniach sejsmicznych w latach 90-tych. Za przewagę nad okienkową transformatą Fouriera uważa się lepszą rozdzielczość czasową.
Korzystając z poniższego fragmentu kodu, można rozłożyć jeden ze śladów sejsmicznych na poszczególne składowe:
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["Wavelet decomposition",Black,32],
LabelStyle->Directive[Black,Italic],
PlotRange->All,
Frame->True]Rozkład śladu na składowe
Aby ocenić rozkład energii odbicia przy różnych czasach nadejścia fali, stosuje się skalogramy (analogicznie do spektrogramu). Z reguły w praktyce nie ma potrzeby analizowania wszystkich komponentów. Zazwyczaj wybierane są komponenty o niskiej, średniej i wysokiej częstotliwości.
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->{"Time","Frequency(Hz)"},LabelStyle->Directive[Black,Bold,14],
ColorFunction->"RustTones",ImageSize->Large]Skalogram. Wynik funkcji
Język Wolfram wykorzystuje funkcję transformacji falkowej . A zastosowanie tej funkcji do całego zestawu śladów zostanie wykonane za pomocą tej funkcji . W tym miejscu warto zwrócić uwagę na jedną z mocnych stron Wolfram Mathematica - możliwość wykorzystania równoległości . W powyższym przykładzie nie ma potrzeby równoległości – wolumen danych nie jest duży, jednak przy pracy ze zbiorami danych eksperymentalnych zawierającymi setki tysięcy śladów jest to konieczność.
tracesCWD=Table[Map[Hilbert[#,0]&,Re[ContinuousWaveletTransform[traces[[i]]][[1]]][[{13,15,18}]]],{i,1,Length@traces}]; Po zastosowaniu funkcji Pojawią się nowe zestawy danych odpowiadające wybranym częstotliwościom. W powyższym przykładzie są to częstotliwości: 38 Hz, 33 Hz, 27 Hz. Doboru częstotliwości najczęściej dokonuje się na podstawie badań – uzyskują one efektywne mapy dla różnych kombinacji częstotliwości i wybierają tę najbardziej pouczającą z punktu widzenia geologa.
Jeśli chcesz podzielić się wynikami ze współpracownikami lub przekazać je klientowi, możesz skorzystać z funkcji SEGYExport[] pakietu GeologyIO
outputdata=seismic3DSEGY;
outputdata["traces",1;;-1]=tracesCWD[[All,3]];
outputdata["textheader"]="Wavelet Decomposition Result";
outputdata["binaryheader","NumberDataTraces"]=Length[tracesCWD[[All,3]]];
SEGYExport["D:result.segy",outputdata];W przypadku trzech takich kostek (składników o niskiej częstotliwości, średniej częstotliwości i wysokiej częstotliwości) do łącznej wizualizacji danych zwykle stosuje się mieszanie RGB. Każdemu komponentowi przypisany jest własny kolor - czerwony, zielony, niebieski. W Wolfram Mathematica można to zrobić za pomocą funkcji .
W rezultacie powstają obrazy, na podstawie których można dokonać interpretacji geologicznej. Zarejestrowane na przekroju meandry pozwalają na wytyczenie paleokanałów, które najprawdopodobniej są zbiornikami i zawierają zasoby ropy. Poszukiwanie i analiza współczesnych analogów takiego systemu rzecznego pozwala określić najbardziej obiecujące fragmenty meandrów. Same kanały charakteryzują się grubymi warstwami dobrze wysortowanego piaskowca i są dobrym zbiornikiem ropy. Obszary poza anomaliami „koronkowymi” są podobne do współczesnych osadów zalewowych. Osady zalewowe reprezentowane są głównie przez skały ilaste i wiercenia w tych strefach będą nieefektywne.
Wycinek RGB kostki danych. W centrum (nieco na lewo od centrum) można prześledzić meandrującą rzekę.
Wycinek RGB kostki danych. Po lewej stronie można prześledzić meandrującą rzekę.
W niektórych przypadkach jakość danych sejsmicznych pozwala uzyskać znacznie wyraźniejsze obrazy. Zależy to od metodyki prac terenowych, sprzętu wykorzystywanego przez algorytm redukcji szumów. Widoczne są wówczas nie tylko fragmenty systemów rzecznych, ale także całe rozbudowane rzeki paleo.
Mieszanie RGB trzech składników kostki danych sejsmicznych (przekrój poziomy). Głębokość około 2 km.
Zdjęcie satelitarne rzeki Wołgi w pobliżu Saratowa
wniosek
Wolfram Mathematica pozwala analizować dane sejsmiczne i rozwiązywać stosowane problemy związane z poszukiwaniem minerałów, a pakiet GeologyIO ułatwia ten proces. Struktura danych sejsmicznych jest taka, że zastosowanie wbudowanych metod przyspieszających obliczenia (, ,…) jest bardzo wydajny i pozwala na przetwarzanie dużej ilości danych. W dużej mierze ułatwiają to funkcje przechowywania danych pakietu GeologyIO. Nawiasem mówiąc, pakiet można wykorzystać nie tylko w zakresie stosowanych badań sejsmicznych. Prawie te same typy danych wykorzystuje się w radarach penetracyjnych i sejsmologii.Jeśli masz sugestie jak poprawić wynik, jakie algorytmy analizy sygnałów z arsenału Wolfram Mathematica można zastosować do takich danych lub masz jakieś krytyczne uwagi, proszę zostaw komentarz.
Źródło: www.habr.com