Děkuji autorovi blogu za jeho zprávu
úvod
Tato poznámka byla napsána v návaznosti na konferenci a obsahuje shrnutí zprávy, kterou jsem podal. Akce se konala v červnu v Petrohradě. Vzhledem k tomu, že pracuji blok od místa konference, nemohl jsem se této akce nezúčastnit. V roce 2016 a 2017 jsem poslouchal zprávy z konferencí a letos jsem měl prezentaci. Za prvé se objevilo zajímavé (zdá se mi) téma, které rozvíjíme a za druhé, po dlouhém studiu legislativy Ruské federace ohledně sankční politiky se v podniku, kde pracuji, objevily až dvě licence .
Než přejdu k tématu mého vystoupení, rád bych poznamenal dobrou organizaci akce. Návštěvní stránka konference používá obrázek kazaňské katedrály. Katedrála je jednou z hlavních atrakcí Petrohradu a je velmi dobře viditelná ze sálu, ve kterém se konference konala.
U vchodu do Petrohradské státní ekonomické univerzity na účastníky čekali asistenti z řad studentů - nedovolili jim zabloudit. Při registraci byly rozdány drobné upomínkové předměty (hračka - blikající bodec, propiska, samolepky se symboly Wolfram). Do programu konference byly zahrnuty i obědy a přestávky na kávu. Už jsem se zmínil o lahodné kávě a koláčích na stěně skupiny - kuchaři jsou skvělí. Touto úvodní částí bych rád zdůraznil, že již samotná akce, její formát a umístění přináší pozitivní emoce.
Zpráva, kterou jsme připravili já a Kirill Belov, se jmenuje „Použití Wolfram Mathematica k řešení problémů v aplikované geofyzice. Spektrální analýza seismických dat nebo „kde tečou staré řeky“. Obsah zprávy zahrnuje dvě části: za prvé, použití algoritmů dostupných v pro analýzu geofyzikálních dat a za druhé, takto lze vložit geofyzikální data do Wolfram Mathematica.
Seismický průzkum
Nejprve si musíte udělat krátkou exkurzi do geofyziky. Geofyzika je věda, která studuje fyzikální vlastnosti hornin. No a jelikož horniny mají různé vlastnosti: elektrické, magnetické, elastické, existují odpovídající metody geofyziky: elektrická prospekce, magnetická prospekce, seismická prospekce... V rámci tohoto článku se budeme podrobněji zabývat pouze seismickou prospekcí. Seismický průzkum je hlavní metodou hledání ropy a plynu. Metoda je založena na buzení elastických vibrací a následném záznamu odezvy od hornin tvořících studovanou oblast. Vibrace jsou buzeny na souši (s dynamitem nebo nevýbušnými zdroji vibrací elastických vibrací) nebo na moři (se vzduchovými zbraněmi). Elastické vibrace se šíří horninovým masivem, lámou se a odrážejí na hranicích vrstev s různými vlastnostmi. Odražené vlny se vracejí na povrch a jsou zaznamenávány geofony na souši (obvykle elektrodynamická zařízení založená na pohybu magnetu zavěšeného v cívce) nebo hydrofony v moři (na základě piezoelektrického jevu). Podle doby příchodu vln lze posuzovat hloubky geologických vrstev.
Tažné zařízení pro seismická plavidla
Vzduchová pistole vyvolává elastické vibrace
Vlny procházejí skalním masivem a jsou zaznamenávány hydrofony
Geofyzikální průzkumná výzkumná loď "Ivan Gubkin" na molu poblíž Blagoveščenského mostu v Petrohradu
Model seismického signálu
Horniny mají různé fyzikální vlastnosti. Pro seismický průzkum jsou důležité především elastické vlastnosti – rychlost šíření elastických vibrací a hustota. Pokud mají dvě vrstvy stejné nebo podobné vlastnosti, pak si vlna „nevšimne“ hranice mezi nimi. Pokud se rychlosti vlnění ve vrstvách liší, dojde k odrazu na hranici vrstev. Čím větší je rozdíl ve vlastnostech, tím intenzivnější je odraz. Jeho intenzita bude určena koeficientem odrazivosti (rc):
kde ρ je hustota horniny, ν je rychlost vlny, 1 a 2 označují horní a spodní vrstvy.
Jedním z nejjednodušších a nejčastěji používaných modelů seismických signálů je konvoluční model, kdy zaznamenaná seismická stopa je reprezentována jako výsledek konvoluce sekvence koeficientů odrazu se snímacím impulsem:
kde s(t) — seismická stopa, tzn. vše, co bylo zaznamenáno hydrofonem nebo geofonem během pevně stanovené doby nahrávání, w(t) - signál generovaný vzduchovou pistolí, n(t) - náhodný hluk.
Vypočítejme syntetickou seismickou stopu jako příklad. Jako počáteční signál použijeme Rickerův puls, široce používaný při seismickém průzkumu.
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]Počáteční seismický impuls
Nastavíme dvě hranice v hloubkách 300 ms a 600 ms a koeficienty odrazu budou náhodná čísla
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]]Posloupnost koeficientů odrazu
Pojďme vypočítat a zobrazit seismickou stopu. Protože koeficienty odrazu mají různá znaménka, dostáváme dva střídající se odrazy na seismické stopě.
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]]Simulovaná trať
U tohoto příkladu je nutné provést rezervaci - v reálu se hloubka vrstev určuje samozřejmě v metrech a výpočet seismické stopy nastává pro časovou oblast. Správnější by bylo nastavit hloubky v metrech a vypočítat doby příchodu se znalostí rychlostí ve vrstvách. V tomto případě rovnou nastavím vrstvy na časovou osu.
Pokud mluvíme o terénním výzkumu, pak je v důsledku takových pozorování zaznamenáno obrovské množství podobných časových řad (seismických stop). Například při studiu místa 25 km dlouhého a 15 km širokého, kde v důsledku práce každá stopa charakterizuje buňku o rozměrech 25x25 metrů (takové buňce se říká bin), bude konečné datové pole obsahovat 600000 1 stop. S dobou vzorkování 5 ms a dobou záznamu 11 sekund bude mít výsledný datový soubor více než XNUMX GB a objem původního „surového“ materiálu může být stovky gigabajtů.
Jak s nimi pracovat ?
Balíček
Vývoj balíčku začal na VK zdi ruskojazyčné podpůrné skupiny. Díky ohlasům komunity se velmi rychle našlo řešení. A v důsledku toho to přerostlo ve vážný vývoj. Odpovídající Dokonce to bylo označeno moderátory. V současné době balíček podporuje práci s následujícími datovými typy, které se aktivně používají v geologickém průmyslu:
- import mapových podkladů ve formátech ZMAP a IRAP
- import měření v jamkách formátu LAS
- vstup a výstup formátu seismických souborů
Pro instalaci balíčku je třeba postupovat podle pokynů na stránce stahování sestaveného balíčku, tzn. v libovolném spusťte následující kód :
If[PacletInformation["GeologyIO"] === {}, PacletInstall[URLDownload[
"https://wolfr.am/FiQ5oFih",
FileNameJoin[{CreateDirectory[], "GeologyIO-0.2.2.paclet"}]
]]]Poté bude balíček nainstalován do výchozí složky, jejíž cestu lze získat následovně:
FileNameJoin[{$UserBasePacletsDirectory, "Repository"}]Jako příklad si ukážeme hlavní schopnosti balíčku. Hovor se provádí tradičně pro balíčky v jazyce Wolfram:
Get["GeologyIO`"]Balíček je vyvinut pomocí . To umožňuje doprovázet hlavní funkcionalitu balíčku dokumentací, která se z hlediska formátu prezentace neliší od dokumentace samotného Wolfram Mathematica, a poskytnout balíčku testovací soubory pro první seznámení.
Takovým souborem je zejména soubor „Marmousi.segy“ - jedná se o syntetický model geologického řezu, který byl vyvinut Francouzským ropným institutem. Pomocí tohoto modelu vývojáři testují své vlastní algoritmy pro modelování vlnového pole, zpracování dat, inverzi seismické stopy atd. Samotný model Marmousi je uložen v úložišti, odkud byl stažen samotný balíček. Chcete-li získat soubor, spusťte následující kód:
If[Not[FileExistsQ["Marmousi.segy"]],
URLDownload["https://wolfr.am/FiQGh7rk", "Marmousi.segy"];]
marmousi = SEGYImport["Marmousi.segy"]Výsledek importu - objekt SEGYData
Formát SEGY zahrnuje ukládání různých informací o pozorováních. Za prvé jsou to textové komentáře. Patří sem informace o místě díla, názvy firem, které měření prováděly atp. V našem případě je tato hlavička volána požadavkem s klíčem TextHeader. Zde je záhlaví zkráceného textu:
Short[marmousi["TextHeader"]]"Datový soubor Marmousi byl generován v institutu ...minimální rychlost 1500 m/s a maximální 5500 m/s)"
Skutečný geologický model můžete zobrazit přístupem k seismickým stopám pomocí klávesy „trasy“ (jednou z funkcí balíčku je, že klíče nerozlišují malá a velká písmena):
ArrayPlot[Transpose[marmousi["traces"]], PlotTheme -> "Detailed"]Modelka Marmousi
Aktuálně balíček také umožňuje načítat data po částech z velkých souborů, čímž je možné zpracovávat soubory, jejichž velikost může dosahovat desítek gigabajtů. Mezi funkce balíčku patří také funkce pro export dat do .segy a částečné připojení na konec souboru.
Samostatně stojí za zmínku funkčnost balíčku při práci se složitou strukturou souborů .segy. Protože umožňuje nejen přistupovat k jednotlivým stopám a hlavičkám pomocí klíčů a indexů, ale také je měnit a následně zapisovat do souboru. Mnoho technických detailů implementace GeologyIO přesahuje rámec tohoto článku a pravděpodobně si zaslouží samostatný popis.
Význam spektrální analýzy při seismickém průzkumu
Schopnost importovat seismická data do Wolfram Mathematica vám umožňuje používat vestavěnou funkci zpracování signálu pro experimentální data. Protože každá seismická stopa představuje časovou řadu, jedním z hlavních nástrojů pro jejich studium je spektrální analýza. Mezi předpoklady pro analýzu frekvenčního složení seismických dat můžeme jmenovat například následující:
- Různé typy vln se vyznačují různým frekvenčním složením. To vám umožní zvýraznit užitečné vlny a potlačit interferenční vlny.
- Vlastnosti hornin, jako je pórovitost a nasycení, mohou ovlivnit frekvenční složení. To umožňuje identifikovat horniny s nejlepšími vlastnostmi.
- Vrstvy s různou tloušťkou způsobují anomálie v různých frekvenčních rozsazích.
Třetí bod je hlavní v kontextu tohoto článku. Níže je uveden fragment kódu pro výpočet seismických stop v případě vrstvy s různou tloušťkou - klínový model. Tento model je tradičně studován při seismickém průzkumu za účelem analýzy interferenčních efektů, když se vlny odražené od mnoha vrstev na sebe navrství.
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 pinch-out formace
Rychlost vlny uvnitř klínu je 4500 m/s, vně klínu 4000 m/s a hustota se předpokládá konstantní 2200 g/cm³. Pro takový model vypočítáme koeficienty odrazu a seismické stopy.
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]Seismické stopy pro klínový model
Posloupnost seismických stop zobrazená na tomto obrázku se nazývá seismická sekce. Jak vidíte, jeho interpretace může být také provedena na intuitivní úrovni, protože geometrie odražených vln jasně odpovídá modelu, který byl specifikován dříve. Pokud budete analyzovat stopy podrobněji, všimnete si, že stopy od 1 do přibližně 30 se neliší - odraz od střechy útvaru a ze dna se navzájem nepřekrývají. Od 31. stopy začnou odrazy rušit. A i když v modelu se koeficienty odrazu horizontálně nemění - seismické stopy mění svou intenzitu se změnou tloušťky útvaru.
Uvažujme amplitudu odrazu od horní hranice útvaru. Počínaje 60. cestou se intenzita odrazu začíná zvyšovat a u 70. cesty se stává maximální. Takto se projevuje interference vlnění ze střechy a dna vrstev, vedoucí v některých případech k výrazným anomáliím v seismickém záznamu.
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]Graf amplitudy odražené vlny od horní hrany klínu
Je logické, že při nižším kmitočtu se rušení začíná objevovat při velkých tloušťkách útvaru a u vysokofrekvenčního signálu při menších tloušťkách. Následující fragment kódu vytváří signál s frekvencemi 35 Hz, 55 Hz a 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]Sada zdrojových signálů s frekvencemi 35 Hz, 55 Hz, 85 Hz
Výpočtem seismických stop a vynesením grafů amplitud odražených vln můžeme vidět, že pro různé frekvence je pozorována anomálie při různých tloušťkách formace.
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]Grafy amplitud odražené vlny od horního okraje klínu pro různé frekvence
Schopnost vyvozovat závěry o tloušťce útvaru z výsledků seismických pozorování je mimořádně užitečná, protože jedním z hlavních úkolů při průzkumu ropy je posoudit nejslibnější body pro položení vrtu (tj. oblasti, kde se formace nachází tlustší). V geologickém řezu se navíc mohou nacházet objekty, jejichž geneze způsobí prudkou změnu mocnosti souvrství. Díky tomu je spektrální analýza účinným nástrojem pro jejich studium. V další části článku se budeme těmito geologickými objekty zabývat podrobněji.
Experimentální data. Kde jsi je sehnal a co v nich hledat?
Materiály analyzované v článku byly získány v západní Sibiři. Region, jak asi každý bez výjimky ví, je hlavním těžařským regionem naší země. Aktivní rozvoj ložisek začal v regionu v 60. letech minulého století. Hlavní metodou hledání ložisek ropy je seismický průzkum. Je zajímavé podívat se na satelitní snímky tohoto území. V malém měřítku si můžete všimnout obrovského množství bažin a jezer; zvětšením mapy můžete vidět místa vrtání klastrových vrtů a zvětšením mapy na maximum můžete také rozlišit mýtiny profilů, podél kterých seismické byla provedena pozorování.
Satelitní snímek map Yandex - oblast města Noyabrsk
Síť podložek studní na jednom z polí
Roponosné horniny západní Sibiře se vyskytují v širokém rozmezí hloubek - od 1 km do 5 km. Hlavní objem hornin obsahujících ropu vznikl v době jury a křídy. Období jury znají pravděpodobně mnozí ze stejnojmenného filmu. se výrazně lišila od té moderní. Encyclopedia Britannica má řadu paleomap, které charakterizují každou helogickou éru.
V dnešní době
jura
Vezměte prosím na vědomí, že v dobách jury bylo území západní Sibiře mořské pobřeží (země protkané řekami a mělkým mořem). Vzhledem k tomu, že klima bylo příjemné, můžeme předpokládat, že typická krajina té doby vypadala takto:
Jurská Sibiř
Na tomto obrázku pro nás nejsou tak důležitá zvířata a ptáci, ale obraz řeky v pozadí. Řeka je stejný geologický objekt, u kterého jsme se zastavili dříve. Činnost řek totiž umožňuje hromadění dobře vytříděných pískovců, které se pak stanou zásobárnou ropy. Tyto nádrže mohou mít bizarní, složitý tvar (jako koryto řeky) a mají proměnlivou tloušťku - u břehů je tloušťka malá, ale blíže ke středu kanálu nebo v meandrech se zvětšuje. Řeky vzniklé v juře jsou tedy nyní v hloubce asi tří kilometrů a jsou předmětem hledání ropných rezervoárů.
Experimentální data. Zpracování a vizualizace
Okamžitě udělejme výhradu k seismickým materiálům uvedeným v článku - vzhledem k tomu, že množství dat použitých pro analýzu je významné - je v textu článku zahrnut pouze fragment původního souboru seismických stop. To umožní komukoli reprodukovat výše uvedené výpočty.
Při práci se seismickými daty geofyzik obvykle používá specializovaný software (existuje několik lídrů v oboru, jejichž vývoj aktivně využívá, například Petrel nebo Paradigm), který vám umožňuje analyzovat různé typy dat a má pohodlné grafické rozhraní. Přes veškerou pohodlnost mají tyto typy softwaru také své nevýhody - například implementace moderních algoritmů ve stabilních verzích zabere spoustu času a možnosti automatizace výpočtů jsou obvykle omezené. V takové situaci je velmi výhodné používat počítačové matematické systémy a programovací jazyky na vysoké úrovni, které umožňují použití širokého algoritmického základu a zároveň zabírají spoustu rutiny. Toto je princip používaný pro práci se seismickými daty ve Wolfram Mathematica. Psát bohatou funkcionalitu pro interaktivní práci s daty je nevhodné – důležitější je zajistit načtení z obecně uznávaného formátu, aplikovat na ně požadované algoritmy a nahrát je zpět do externího formátu.
Podle navrženého schématu načteme původní seismická data a zobrazíme 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]Takto stažená a importovaná data jsou trasy zaznamenané na ploše o rozměrech 10 krát 5 kilometrů. Pokud jsou data získávána metodou trojrozměrného seismického průzkumu (vlny nejsou zaznamenávány podél jednotlivých geofyzikálních profilů, ale po celém území současně), je možné získat seismické datové kostky. Jedná se o trojrozměrné objekty, jejichž vertikální a horizontální řezy umožňují podrobné studium geologického prostředí. V uvažovaném příkladu máme co do činění s trojrozměrnými daty. Můžeme získat nějaké informace z textové hlavičky, jako je tato
StringPartition[seismic3DSEGY["textheader"], 80] // TableFormC 1 TOTO JE DEMO SOUBOR PRO TEST GEOLOGICKÉHO BALÍČKU
C 2
C 3
C 4
C 5 DATUM UŽIVATELSKÉ JMÉNO: WOLFRAM USER
C 6 NÁZEV PRŮZKUMU: NĚKDE NA SIBIŘI
C 7 TYP SOUBORU 3D SEISMICKÝ OBJEM
C 8
C 9
ŘADA C10 Z: PRVNÍ 2200M POSLEDNÍ 2400M
Tento soubor dat nám postačí k demonstraci hlavních fází analýzy dat. Stopy v souboru jsou zaznamenávány sekvenčně a každá z nich vypadá asi jako na následujícím obrázku - to je rozložení amplitud odražených vln podél svislé osy (osa hloubky).
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]Jedna ze stop seismického řezu
Když víte, kolik stop se nachází v každém směru studované oblasti, můžete vygenerovat trojrozměrné datové pole a zobrazit je pomocí funkce 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]]XNUMXD obrázek seismické datové krychle (svislá osa - hloubka)
Pokud geologické prvky zájmu vytvářejí intenzivní seismické anomálie, lze použít transparentní vizualizační nástroje. „Nedůležité“ oblasti záznamu lze zneviditelnit, takže budou viditelné pouze anomálie. Ve Wolfram Mathematica to lze provést 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 seismické datové krychle pomocí funkcí Opacity[] a Raster3D[]
Stejně jako v syntetickém příkladu lze na řezech původní krychle identifikovat některé geologické hranice (vrstvy) s proměnlivým reliéfem.
Hlavním nástrojem pro spektrální analýzu je Fourierova transformace. S jeho pomocí můžete vyhodnotit amplitudově-frekvenční spektrum každé stopy nebo skupiny stop. Po přenosu dat do frekvenční domény se však ztratí informace o tom, v jakých časech (čti v jakých hloubkách) se frekvence mění. Aby bylo možné lokalizovat změny signálu na časové (hloubkové) ose, používá se okénková Fourierova transformace a vlnkový rozklad. Tento článek používá vlnkový rozklad. Technologie vlnové analýzy se začala aktivně používat při seismickém průzkumu v 90. letech. Za výhodu oproti okénkové Fourierově transformaci se považuje lepší časové rozlišení.
Pomocí následujícího fragmentu kódu můžete rozložit jednu ze seismických stop na jednotlivé komponenty:
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]Rozklad stopy na složky
K posouzení, jak je energie odrazu distribuována při různých dobách příchodu vlny, se používají skalogramy (analogické ke spektrogramu). V praxi zpravidla není potřeba analyzovat všechny komponenty. Obvykle se volí nízkofrekvenční, střední a vysokofrekvenční složky.
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. Výsledek funkce
Wolfram Language používá funkci pro vlnkovou transformaci . A aplikace této funkce na celou sadu tras bude provedena pomocí funkce . Zde stojí za zmínku jedna ze silných stránek Wolfram Mathematica – schopnost používat paralelizaci . Ve výše uvedeném příkladu není potřeba paralelizace – objem dat není velký, ale při práci s experimentálními datovými sadami obsahujícími statisíce stop je to nutnost.
tracesCWD=Table[Map[Hilbert[#,0]&,Re[ContinuousWaveletTransform[traces[[i]]][[1]]][[{13,15,18}]]],{i,1,Length@traces}]; Po aplikaci funkce Objeví se nové datové sady odpovídající vybraným frekvencím. Ve výše uvedeném příkladu jsou tyto frekvence: 38Hz, 33Hz, 27Hz. Volba frekvencí probíhá nejčastěji na základě testování – získávají efektivní mapy pro různé kombinace frekvencí a vybírají tu nejinformativnější z pohledu geologa.
Pokud potřebujete výsledky sdílet s kolegy nebo je poskytnout zákazníkovi, můžete využít funkci SEGYExport[] balíčku 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];Se třemi z těchto krychlí (nízkofrekvenční, středofrekvenční a vysokofrekvenční složky) se k vizualizaci dat obvykle používá prolnutí RGB. Každá součástka má přiřazenou svou barvu – červená, zelená, modrá. Ve Wolfram Mathematica to lze provést pomocí funkce .
Výsledkem jsou snímky, ze kterých lze provést geologickou interpretaci. Meandry, které jsou na úseku zaznamenány, umožňují vymezit paleochanály, které jsou spíše rezervoáry a obsahují zásoby ropy. Hledání a analýza moderních analogů takového říčního systému nám umožňuje určit nejslibnější části meandrů. Samotné kanály se vyznačují silnými vrstvami dobře tříděného pískovce a jsou dobrým rezervoárem ropy. Oblasti mimo „krajkové“ anomálie jsou podobné moderním náplavovým ložiskům. Lužní ložiska jsou zastoupena především jílovitými horninami a vrty do těchto zón budou neúčinné.
RGB řez datové krychle. Uprostřed (trochu vlevo od středu) můžete sledovat meandrující řeku.
RGB řez datové krychle. Po levé straně můžete sledovat meandrující řeku.
V některých případech umožňuje kvalita seismických dat výrazně jasnější snímky. To závisí na metodice práce v terénu, zařízení používaném algoritmem redukce šumu. V takových případech jsou vidět nejen fragmenty říčních systémů, ale i celé rozšířené paleoříky.
RGB míchání tří složek seismické datové krychle (horizontální řez). Hloubka cca 2 km.
Satelitní snímek řeky Volhy poblíž Saratova
Závěr
Wolfram Mathematica vám umožňuje analyzovat seismická data a řešit aplikované problémy související s průzkumem nerostů a balíček GeologyIO tento proces usnadňuje. Struktura seismických dat je taková, že použití vestavěných metod urychluje výpočty (, ,…) je velmi efektivní a umožňuje zpracovávat velké množství dat. To je do značné míry usnadněno funkcemi ukládání dat balíčku GeologyIO. Mimochodem, balíček lze využít nejen v oblasti aplikovaného seismického průzkumu. Téměř stejné typy dat se používají v pozemních radarech a seismologii. Pokud máte návrhy, jak zlepšit výsledek, jaké algoritmy analýzy signálů z arzenálu Wolfram Mathematica lze na taková data použít, nebo pokud máte nějaké kritické připomínky, prosím Zanechat komentář.
Zdroj: www.habr.com