Takk til forfatteren av bloggen for rapporten hans
Innledning
Dette notatet ble skrevet i kjølvannet av konferansen og inneholder et sammendrag av rapporten jeg ga. Arrangementet fant sted i juni i St. Petersburg. Med tanke på at jeg jobber et stykke unna konferansesiden, kunne jeg ikke la være å delta på dette arrangementet. I 2016 og 2017 hørte jeg på konferanserapporter, og i år holdt jeg en presentasjon. For det første har det dukket opp et interessant (synes det for meg) tema, som vi utvikler oss med , og for det andre, etter en lang studie av lovgivningen i Den russiske føderasjonen angående sanksjonspolitikk, ved bedriften der jeg jobber, dukket det opp så mange som to lisenser .
Før jeg går over til temaet for talen min, vil jeg bemerke den gode organiseringen av arrangementet. Konferansens besøksside bruker et bilde av Kazan-katedralen. Katedralen er en av hovedattraksjonene i St. Petersburg og er veldig godt synlig fra salen der konferansen fant sted.
Ved inngangen til St. Petersburg State Economic University ble deltakerne møtt av assistenter blant studentene - de lot dem ikke gå seg vill. Under registreringen ble det gitt ut små suvenirer (et leketøy - en blinkende pigg, en penn, klistremerker med Wolfram-symboler). Lunsj og kaffepauser var også inkludert i konferanseplanen. Jeg har allerede notert om deilig kaffe og paier på gruppens vegg - kokkene er flotte. Med denne innledende delen vil jeg understreke at selve arrangementet, dets format og plassering allerede bringer positive følelser.
Rapporten som ble utarbeidet av meg og Kirill Belov heter «Using Wolfram Mathematica for å løse problemer i anvendt geofysikk. Spektralanalyse av seismiske data eller "hvor eldgamle elver rant." Innholdet i rapporten dekker to deler: For det første bruk av algoritmer tilgjengelig i for å analysere geofysiske data, og for det andre er dette hvordan man legger geofysiske data inn i Wolfram Mathematica.
Seismisk leting
Først må du gjøre en kort ekskursjon i geofysikk. Geofysikk er vitenskapen som studerer de fysiske egenskapene til bergarter. Vel, siden bergarter har forskjellige egenskaper: elektriske, magnetiske, elastiske, er det tilsvarende metoder for geofysikk: elektrisk prospektering, magnetisk prospektering, seismisk prospektering... I sammenheng med denne artikkelen vil vi bare diskutere seismisk prospektering mer detaljert. Seismisk leting er hovedmetoden for å lete etter olje og gass. Metoden er basert på eksitering av elastiske vibrasjoner og påfølgende registrering av responsen fra bergartene som utgjør studieområdet. Vibrasjoner utløses på land (med dynamitt eller ikke-eksplosive vibrasjonskilder for elastiske vibrasjoner) eller til sjøs (med luftkanoner). Elastiske vibrasjoner forplanter seg gjennom bergmassen, brytes og reflekteres ved grensene til lag med forskjellige egenskaper. Reflekterte bølger går tilbake til overflaten og registreres av geofoner på land (vanligvis elektrodynamiske enheter basert på bevegelsen til en magnet suspendert i en spole) eller hydrofoner i havet (basert på den piezoelektriske effekten). Ved tidspunktet for ankomst av bølger kan man bedømme dybden av geologiske lag.
Slepeutstyr for seismisk fartøy
Luftpistolen stimulerer elastiske vibrasjoner
Bølgene går gjennom bergmassen og registreres av hydrofoner
Geofysisk undersøkelsesfartøy "Ivan Gubkin" ved brygga nær Blagoveshchensky-broen i St. Petersburg
Seismisk signalmodell
Bergarter har forskjellige fysiske egenskaper. For seismisk leting er elastiske egenskaper først og fremst viktige - hastigheten på forplantning av elastiske vibrasjoner og tetthet. Hvis to lag har samme eller lignende egenskaper, "vil ikke bølgen legge merke til" grensen mellom dem. Hvis bølgehastighetene i lagene er forskjellige, vil refleksjon oppstå ved grensen til lagene. Jo større forskjell i egenskaper, jo mer intens blir refleksjonen. Dens intensitet vil bli bestemt av refleksjonskoeffisienten (rc):
hvor ρ er bergartens tetthet, ν er bølgehastigheten, 1 og 2 indikerer øvre og nedre lag.
En av de enkleste og mest brukte seismiske signalmodellene er konvolusjonsmodellen, når den registrerte seismiske trasen er representert som et resultat av konvolusjon av en sekvens av refleksjonskoeffisienter med en sonderende puls:
hvor s(t) — seismisk spor, dvs. alt som ble tatt opp av en hydrofon eller geofon i løpet av en fast opptakstid, w(t) - signalet som genereres av luftpistolen, n(t) - tilfeldig støy.
La oss beregne et syntetisk seismisk spor som et eksempel. Vi vil bruke Ricker-pulsen, mye brukt i seismisk leting, som startsignal.
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]Innledende seismisk impuls
Vi vil sette to grenser på dybder på 300 ms og 600 ms, og refleksjonskoeffisientene vil være tilfeldige tall
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]]Sekvens av refleksjonskoeffisienter
La oss beregne og vise det seismiske sporet. Siden refleksjonskoeffisientene har forskjellige fortegn, får vi to vekslende refleksjoner på det seismiske sporet.
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]]Simulert spor
For dette eksemplet er det nødvendig å reservere seg - i virkeligheten bestemmes dybden av lagene, selvfølgelig, i meter, og beregningen av det seismiske sporet skjer for tidsdomenet. Det ville være mer riktig å sette dybdene i meter og beregne ankomsttidene med kjennskap til hastighetene i lagene. I dette tilfellet setter jeg umiddelbart lagene på tidsaksen.
Hvis vi snakker om feltforskning, blir et stort antall lignende tidsserier (seismiske spor) registrert som et resultat av slike observasjoner. For eksempel, når man studerer et sted som er 25 km langt og 15 km bredt, hvor hvert spor som et resultat av arbeid karakteriserer en celle som måler 25x25 meter (en slik celle kalles en bin), vil den endelige datamatrisen inneholde 600000 1 spor. Med en samplingstid på 5 ms og en opptakstid på 11 sekunder, vil den endelige datafilen være mer enn XNUMX GB, og volumet av det originale "råmaterialet" kan være hundrevis av gigabyte.
Hvordan jobbe med dem ?
pakke
Utviklingen av pakken startet på VK-veggen til den russisktalende støttegruppen. Takket være fellesskapets svar ble en løsning funnet veldig raskt. Og som et resultat vokste det til en alvorlig utvikling. Tilsvarende Det ble til og med markert av moderatorer. For øyeblikket støtter pakken arbeid med følgende datatyper som brukes aktivt i den geologiske industrien:
- import av kartdata i ZMAP- og IRAP-formater
- import av målinger i LAS-format brønner
- input og output av seismiske filer format
For å installere pakken må du følge instruksjonene på nedlastingssiden til den sammensatte pakken, dvs. kjør følgende kode i en hvilken som helst :
If[PacletInformation["GeologyIO"] === {}, PacletInstall[URLDownload[
"https://wolfr.am/FiQ5oFih",
FileNameJoin[{CreateDirectory[], "GeologyIO-0.2.2.paclet"}]
]]]Deretter vil pakken bli installert i standardmappen, banen som kan hentes som følger:
FileNameJoin[{$UserBasePacletsDirectory, "Repository"}]Som et eksempel vil vi demonstrere hovedfunksjonene til pakken. Samtalen gjøres tradisjonelt for pakker på Wolfram-språket:
Get["GeologyIO`"]Pakken er utviklet vha . Dette lar deg følge hovedfunksjonaliteten til pakken med dokumentasjon, som når det gjelder presentasjonsformat ikke skiller seg fra dokumentasjonen til Wolfram Mathematica selv, og gi pakken testfiler for første bekjentskap.
Spesielt en slik fil er filen "Marmousi.segy" - dette er en syntetisk modell av en geologisk seksjon, som ble utviklet av det franske petroleumsinstituttet. Ved å bruke denne modellen tester utviklere sine egne algoritmer for bølgefeltmodellering, databehandling, seismisk sporinversjon, etc. Selve Marmousi-modellen er lagret i depotet der selve pakken ble lastet ned. For å få filen, kjør følgende kode:
If[Not[FileExistsQ["Marmousi.segy"]],
URLDownload["https://wolfr.am/FiQGh7rk", "Marmousi.segy"];]
marmousi = SEGYImport["Marmousi.segy"]Importer resultat - SEGYData-objekt
SEGY-formatet innebærer å lagre ulike opplysninger om observasjoner. For det første er dette tekstkommentarer. Dette inkluderer informasjon om plasseringen av arbeidet, navnene på firmaene som utførte målingene mv. I vårt tilfelle kalles denne overskriften opp av en forespørsel med TextHeader-tasten. Her er en forkortet tekstoverskrift:
Short[marmousi["TextHeader"]]"Marmousi-datasettet ble generert ved instituttet ... minimumshastighet på 1500 m/s og maksimalt 5500 m/s)"
Du kan vise den faktiske geologiske modellen ved å få tilgang til de seismiske sporene ved å bruke "spor"-tasten (en av funksjonene til pakken er at nøklene ikke skiller mellom store og små bokstaver):
ArrayPlot[Transpose[marmousi["traces"]], PlotTheme -> "Detailed"]Modell Marmousi
Foreløpig lar pakken deg også laste inn data i deler fra store filer, noe som gjør det mulig å behandle filer hvis størrelse kan nå titalls gigabyte. Pakkens funksjoner inkluderer også funksjoner for å eksportere data til .segy og delvis legge til på slutten av filen.
Separat er det verdt å merke seg funksjonaliteten til pakken når du arbeider med den komplekse strukturen til .segy-filer. Siden det lar deg ikke bare få tilgang til individuelle spor og overskrifter ved hjelp av nøkler og indekser, men også å endre dem og deretter skrive dem til en fil. Mange av de tekniske detaljene i GeologyIOs implementering er utenfor rammen av denne artikkelen og fortjener sannsynligvis en egen beskrivelse.
Relevans av spektralanalyse i seismisk leting
Muligheten til å importere seismiske data til Wolfram Mathematica lar deg bruke innebygd signalbehandlingsfunksjonalitet for eksperimentelle data. Siden hvert seismisk spor representerer en tidsserie, er et av hovedverktøyene for å studere dem spektralanalyse. Blant forutsetningene for å analysere frekvenssammensetningen til seismiske data, kan vi for eksempel nevne følgende:
- Ulike typer bølger er preget av ulik frekvenssammensetning. Dette lar deg fremheve nyttige bølger og undertrykke interferensbølger.
- Bergartsegenskaper som porøsitet og metning kan påvirke frekvenssammensetningen. Dette gjør det mulig å identifisere bergarter med de beste egenskapene.
- Lag med forskjellig tykkelse forårsaker anomalier i forskjellige frekvensområder.
Det tredje punktet er det viktigste i sammenheng med denne artikkelen. Nedenfor er et kodefragment for beregning av seismiske spor ved lag med varierende tykkelse - en kilemodell. Denne modellen er tradisjonelt studert i seismisk utforskning for å analysere interferenseffekter når bølger som reflekteres fra mange lag legges over hverandre.
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]]Modell av en pinch-out formasjon
Bølgehastigheten inne i kilen er 4500 m/s, utenfor kilen 4000 m/s, og tettheten antas å være konstant 2200 g/cm³. For en slik modell beregner vi refleksjonskoeffisienter og seismiske spor.
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]Seismiske spor for kilemodellen
Rekkefølgen av seismiske spor vist i denne figuren kalles et seismisk snitt. Som du kan se, kan tolkningen også utføres på et intuitivt nivå, siden geometrien til de reflekterte bølgene helt klart tilsvarer modellen som ble spesifisert tidligere. Hvis du analyserer sporene mer detaljert, vil du legge merke til at spor fra 1 til omtrent 30 ikke er forskjellige - refleksjonen fra taket av formasjonen og fra bunnen overlapper ikke hverandre. Fra det 31. sporet begynner refleksjonene å forstyrre. Og selv om refleksjonskoeffisientene i modellen ikke endres horisontalt - de seismiske sporene endrer intensiteten når tykkelsen på formasjonen endres.
La oss vurdere amplituden av refleksjon fra den øvre grensen til formasjonen. Fra den 60. ruten begynner intensiteten av refleksjonen å øke og ved den 70. ruten blir den maksimal. Dette er hvordan interferensen av bølger fra taket og bunnen av lagene manifesterer seg, noe som i noen tilfeller fører til betydelige anomalier i den seismiske registreringen.
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 over amplituden til den reflekterte bølgen fra den øvre kanten av kilen
Det er logisk at når signalet er lavere frekvens, begynner interferens å vises ved store formasjonstykkelser, og i tilfellet med et høyfrekvent signal oppstår interferens ved mindre tykkelser. Følgende kodebit lager et signal med frekvenser på 35 Hz, 55 Hz og 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]Et sett med kildesignaler med frekvenser på 35 Hz, 55Hz, 85Hz
Ved å beregne seismiske spor og plotte grafer av reflekterte bølgeamplituder, kan vi se at det for forskjellige frekvenser observeres en anomali ved forskjellige formasjonstykkelser.
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]Grafer over amplitudene til den reflekterte bølgen fra den øvre kanten av kilen for forskjellige frekvenser
Evnen til å trekke konklusjoner om tykkelsen av formasjonen fra resultatene av seismiske observasjoner er ekstremt nyttig, fordi en av hovedoppgavene i oljeleting er å vurdere de mest lovende punktene for å legge en brønn (dvs. de områdene der formasjonen er tykkere). I tillegg kan det i den geologiske delen være gjenstander hvis tilblivelse forårsaker en skarp endring i tykkelsen på formasjonen. Dette gjør spektralanalyse til et effektivt verktøy for å studere dem. I neste del av artikkelen vil vi vurdere slike geologiske objekter mer detaljert.
Eksperimentelle data. Hvor fikk du tak i dem og hva skal du se etter i dem?
Materialene som ble analysert i artikkelen ble hentet i Vest-Sibir. Regionen, som alle uten unntak sikkert vet, er den viktigste oljeproduserende regionen i landet vårt. Aktiv utvikling av forekomster begynte i regionen på 60-tallet av forrige århundre. Hovedmetoden for å lete etter oljeforekomster er seismisk leting. Det er interessant å se på satellittbilder av dette territoriet. I liten skala kan du se et stort antall sumper og innsjøer; ved å forstørre kartet kan du se klyngebrønnboringssteder, og ved å forstørre kartet til det ytterste kan du også skille mellom lysningene av profilene langs hvilke seismikk observasjoner ble utført.
Satellittbilde av Yandex-kart - Noyabrsk byområde
Et nettverk av brønnputer ved et av feltene
Oljeførende bergarter i Vest-Sibir forekommer i et bredt spekter av dybder - fra 1 km til 5 km. Hovedvolumet av bergarter som inneholder olje ble dannet i jura- og kritttiden. Juratiden er nok kjent for mange fra filmen med samme navn. var vesentlig forskjellig fra den moderne. Encyclopedia Britannica har en serie paleomap som karakteriserer hver helogisk epoke.
Nåværende tid
Jura perioden
Vær oppmerksom på at i juratiden var territoriet til Vest-Sibir en havkyst (land krysset av elver og et grunt hav). Siden klimaet var behagelig, kan vi anta at et typisk landskap på den tiden så slik ut:
Jurassic Sibir
I dette bildet er det viktige for oss ikke så mye dyrene og fuglene, men bildet av elven i bakgrunnen. Elva er det samme geologiske objektet som vi stoppet ved tidligere. Faktum er at aktiviteten til elver gjør at velsorterte sandsteiner kan samle seg, som da vil bli et reservoar for olje. Disse reservoarene kan ha en bisarr, kompleks form (som et elveleie) og de har variabel tykkelse - nær bredden er tykkelsen liten, men nærmere midten av kanalen eller i meanderområder øker den. Så elvene dannet i juraen er nå på en dybde på omtrent tre kilometer og er gjenstand for et søk etter oljereservoarer.
Eksperimentelle data. Bearbeiding og visualisering
La oss umiddelbart ta forbehold om de seismiske materialene vist i artikkelen - på grunn av at mengden data som brukes til analysen er betydelig - er bare et fragment av det opprinnelige settet med seismiske spor inkludert i artikkelens tekst. Dette vil tillate hvem som helst å reprodusere beregningene ovenfor.
Når du arbeider med seismiske data, bruker en geofysiker vanligvis spesialisert programvare (det er flere industriledere hvis utviklinger brukes aktivt, for eksempel Petrel eller Paradigm), som lar deg analysere forskjellige typer data og har et praktisk grafisk grensesnitt. Til tross for all bekvemmeligheten har denne typen programvare også sine ulemper - for eksempel tar implementeringen av moderne algoritmer i stabile versjoner mye tid, og mulighetene for å automatisere beregninger er vanligvis begrenset. I en slik situasjon blir det veldig praktisk å bruke datamatematikksystemer og programmeringsspråk på høyt nivå, som tillater bruk av en bred algoritmisk base og samtidig tar på seg mye rutine. Dette er prinsippet som brukes for å arbeide med seismiske data i Wolfram Mathematica. Det er upassende å skrive rik funksjonalitet for interaktivt arbeid med data - det er viktigere å sikre lasting fra et generelt akseptert format, bruke de ønskede algoritmene til dem og laste dem opp tilbake til et eksternt format.
Etter det foreslåtte opplegget vil vi laste de originale seismiske dataene og vise dem inn :
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]Dataene som lastes ned og importeres på denne måten er rutene registrert på et område som måler 10 x 5 kilometer. Hvis dataene innhentes ved hjelp av en tredimensjonal seismisk undersøkelsesmetode (bølger registreres ikke langs individuelle geofysiske profiler, men over hele området samtidig), blir det mulig å få seismiske datakuber. Dette er tredimensjonale objekter, hvorav vertikale og horisontale deler tillater en detaljert studie av det geologiske miljøet. I det betraktede eksemplet har vi å gjøre med tredimensjonale data. Vi kan få litt informasjon fra tekstoverskriften, som dette
StringPartition[seismic3DSEGY["textheader"], 80] // TableFormC 1 DETTE ER DEMO-FIL FOR GEOLOGYIO-PAKKETEST
C 2
C 3
C 4
C 5 DATO BRUKERNAVN: WOLFRAM BRUKER
C 6 UNDERSØKELSENS NAVN: ET STED I SIBERIA
C 7 FILTYPE 3D SEISMISK VOLUM
C 8
C 9
C10 Z REKKEVIDDE: FØRSTE 2200M SISTE 2400M
Dette datasettet vil være nok for oss til å demonstrere hovedstadiene i dataanalyse. Sporene i filen registreres sekvensielt, og hver av dem ser omtrent ut som følgende figur - dette er fordelingen av amplitudene til reflekterte bølger langs den vertikale aksen (dybdeaksen).
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]Et av de seismiske snittsporene
Når du vet hvor mange spor som er plassert i hver retning av det studerte området, kan du generere en tredimensjonal datamatrise og vise den ved hjelp av Image3D[]-funksjonen
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-bilde av en seismisk datakube. (Vertikal akse - dybde)
Hvis de geologiske trekkene av interesse skaper intense seismiske anomalier, kan visualiseringsverktøy med gjennomsiktighet brukes. "Uviktige" områder av opptaket kan gjøres usynlige, slik at bare uregelmessigheter er synlige. I Wolfram Mathematica kan dette gjøres ved hjelp av и .
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]Seismisk datakubebilde ved hjelp av Opacity[] og Raster3D[] funksjoner
Som i det syntetiske eksemplet kan man på deler av den opprinnelige kuben identifisere noen geologiske grenser (lag) med variabelt relieff.
Hovedverktøyet for spektralanalyse er Fourier-transformasjonen. Med dens hjelp kan du evaluere amplitude-frekvensspekteret til hvert spor eller gruppe av spor. Etter overføring av dataene til frekvensdomenet går det imidlertid tapt informasjon om på hvilke tidspunkter (les i hvilke dybder) frekvensen endres. For å kunne lokalisere signalendringer på tids (dybde) aksen, brukes vindusert Fourier transformasjon og wavelet dekomponering. Denne artikkelen bruker wavelet-dekomponering. Wavelet-analyseteknologi begynte å bli aktivt brukt i seismisk leting på 90-tallet. Fordelen i forhold til den vindusbaserte Fourier-transformasjonen anses å være bedre tidsoppløsning.
Ved å bruke følgende kodefragment kan du dekomponere et av de seismiske sporene til individuelle komponenter:
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]Dekomponering av et spor til komponenter
For å vurdere hvordan refleksjonsenergien er fordelt ved ulike bølgeankomsttider, brukes skalogrammer (analogt med et spektrogram). Som regel er det i praksis ikke nødvendig å analysere alle komponentene. Vanligvis velges lav-, mellom- og høyfrekvente komponenter.
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. Funksjonsresultat
Wolfram-språket bruker funksjonen for wavelet-transformasjon . Og bruken av denne funksjonen på hele settet med spor vil bli utført ved hjelp av funksjonen . Her er det verdt å merke seg en av styrkene til Wolfram Mathematica – evnen til å bruke parallellisering . I eksemplet ovenfor er det ikke behov for parallellisering - datavolumet er ikke stort, men når man jobber med eksperimentelle datasett som inneholder hundretusenvis av spor, er dette en nødvendighet.
tracesCWD=Table[Map[Hilbert[#,0]&,Re[ContinuousWaveletTransform[traces[[i]]][[1]]][[{13,15,18}]]],{i,1,Length@traces}]; Etter å ha brukt funksjonen Nye datasett vises tilsvarende de valgte frekvensene. I eksemplet ovenfor er disse frekvensene: 38Hz, 33Hz, 27Hz. Valget av frekvenser utføres oftest på grunnlag av testing - de får effektive kart for forskjellige frekvenskombinasjoner og velger den mest informative fra en geologs synspunkt.
Hvis du trenger å dele resultatene med kolleger eller gi dem til kunden, kan du bruke funksjonen SEGYExport[] til GeologyIO-pakken
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];Med tre av disse kubene (lavfrekvente, mellomfrekvente og høyfrekvente komponenter), brukes RGB-blanding vanligvis til å visualisere dataene sammen. Hver komponent er tildelt sin egen farge - rød, grønn, blå. I Wolfram Mathematica kan dette gjøres ved hjelp av funksjonen .
Resultatet er bilder som kan foretas geologisk tolkning. Slyngningene som er registrert på seksjonen gjør det mulig å avgrense paleokanaler, som mer sannsynlig er reservoarer og inneholder oljereserver. Søket og analysen av moderne analoger av et slikt elvesystem lar oss bestemme de mest lovende delene av buktningene. Selve kanalene er preget av tykke lag med velsortert sandstein og er et godt reservoar for olje. Områder utenfor "snøre"-anomaliene ligner på moderne flommarkavsetninger. Flomflateavsetninger er hovedsakelig representert av leirholdige bergarter, og boring i disse sonene vil være ineffektivt.
RGB-stykke av datakuben. I midten (litt til venstre for midten) kan du spore den buktende elven.
RGB-stykke av datakuben. På venstre side kan du spore den slyngende elven.
I noen tilfeller gir kvaliteten på seismiske data betydelig klarere bilder. Dette avhenger av feltarbeidsmetodikken, utstyret som brukes av støyreduksjonsalgoritmen. I slike tilfeller er ikke bare fragmenter av elvesystemer synlige, men også hele utvidede paleo-elver.
RGB-blanding av tre komponenter i en seismisk datakube (horisontal skive). Dybde ca 2 km.
Satellittbilde av Volga-elven nær Saratov
Konklusjon
Wolfram Mathematica lar deg analysere seismiske data og løse anvendte problemer knyttet til mineralutforskning, og GeologyIO-pakken gjør denne prosessen mer praktisk. Strukturen til seismiske data er slik at bruk av innebygde metoder for å fremskynde beregninger (, ,...) er svært effektiv og lar deg behandle store mengder data. I stor grad tilrettelegges dette av datalagringsfunksjonene til GeologyIO-pakken. Forresten kan pakken brukes ikke bare innen anvendt seismisk leting. Nesten de samme typene data brukes i jordpenetrerende radar og seismologi. Hvis du har forslag til hvordan du kan forbedre resultatet, hvilke signalanalysealgoritmer fra Wolfram Mathematica-arsenalet som kan brukes på slike data, eller hvis du har noen kritiske kommentarer, vennligst Legg igjen en kommentar.
Kilde: www.habr.com