Wolfram Mathematica i geofysik

Tack till bloggens författare Anton Ekimenko för hans rapport

Inledning

Denna anteckning skrevs i kölvattnet av konferensen Wolfram rysk teknikkonferens och innehåller en sammanfattning av rapporten jag gav. Evenemanget ägde rum i juni i St. Petersburg. Med tanke på att jag jobbar ett kvarter från konferensplatsen kunde jag inte låta bli att delta i detta evenemang. Under 2016 och 2017 lyssnade jag på konferensrapporter och i år höll jag en presentation. För det första har det dykt upp ett intressant (förefaller det mig) ämne som vi utvecklar med Kirill Belov, och för det andra, efter en lång studie av Ryska federationens lagstiftning angående sanktionspolitik, på företaget där jag arbetar, dök upp så många som två licenser upp Wolfram Mathematica.

Innan jag går vidare till ämnet för mitt anförande vill jag notera den goda organisationen av evenemanget. Konferensens besökssida använder en bild av Kazan-katedralen. Katedralen är en av S:t Petersburgs främsta attraktioner och är mycket tydligt synlig från hallen där konferensen ägde rum.

Wolfram Mathematica i geofysik

Vid ingången till St. Petersburg State Economic University möttes deltagarna av assistenter bland studenterna - de tillät dem inte att gå vilse. Under registreringen delades små souvenirer ut (en leksak - en blinkande spik, en penna, klistermärken med Wolfram-symboler). Även lunch och fika ingick i konferensschemat. Jag har redan noterat om utsökt kaffe och pajer på gruppens vägg - kockarna är fantastiska. Med den här inledande delen vill jag betona att själva evenemanget, dess format och plats redan väcker positiva känslor.

Rapporten som utarbetades av mig och Kirill Belov heter ”Att använda Wolfram Mathematica för att lösa problem inom tillämpad geofysik. Spektralanalys av seismiska data eller "där gamla floder rann." Innehållet i rapporten omfattar två delar: för det första användningen av algoritmer tillgängliga i Wolfram Mathematica för att analysera geofysiska data, och för det andra är detta hur man lägger in geofysiska data i Wolfram Mathematica.

Seismisk utforskning

Först måste du göra en kort utflykt till geofysik. Geofysik är vetenskapen som studerar de fysiska egenskaperna hos bergarter. Jo, eftersom bergarter har olika egenskaper: elektriska, magnetiska, elastiska, finns det motsvarande metoder för geofysik: elektrisk prospektering, magnetisk prospektering, seismisk prospektering... I samband med denna artikel kommer vi bara att diskutera seismisk prospektering mer i detalj. Seismisk prospektering är den huvudsakliga metoden för att leta efter olja och gas. Metoden är baserad på excitation av elastiska vibrationer och efterföljande registrering av responsen från de stenar som utgör studieområdet. Vibrationer exciteras på land (med dynamit eller icke-explosiva vibrationskällor för elastiska vibrationer) eller till sjöss (med luftpistoler). Elastiska vibrationer fortplantar sig genom bergmassan, bryts och reflekteras vid gränserna för lager med olika egenskaper. Reflekterade vågor återvänder till ytan och registreras av geofoner på land (vanligtvis elektrodynamiska enheter baserade på rörelsen av en magnet upphängd i en spole) eller hydrofoner i havet (baserat på den piezoelektriska effekten). Vid tidpunkten för vågornas ankomst kan man bedöma djupen av geologiska skikt.

Bogserutrustning för seismisk fartyg
Wolfram Mathematica i geofysik

Luftpistolen exciterar elastiska vibrationer
Wolfram Mathematica i geofysik

Vågorna passerar genom bergmassan och registreras av hydrofoner
Wolfram Mathematica i geofysik

Geofysisk undersökningsforskningsfartyg "Ivan Gubkin" vid piren nära Blagoveshchensky-bron i St. Petersburg
Wolfram Mathematica i geofysik

Seismisk signalmodell

Bergarter har olika fysiska egenskaper. För seismisk utforskning är elastiska egenskaper i första hand viktiga - hastigheten för utbredning av elastiska vibrationer och densitet. Om två lager har samma eller liknande egenskaper, kommer vågen inte att "märka" gränsen mellan dem. Om våghastigheterna i skikten skiljer sig, kommer reflektion att ske vid gränsen av skikten. Ju större skillnaden är i egenskaper, desto intensivare blir reflektionen. Dess intensitet kommer att bestämmas av reflektanskoefficienten (rc):

Wolfram Mathematica i geofysik

där ρ är stendensiteten, ν är våghastigheten, 1 och 2 anger de övre och nedre skikten.

En av de enklaste och mest använda seismiska signalmodellerna är faltningsmodellen, när det registrerade seismiska spåret representeras som resultatet av faltning av en sekvens av reflektionskoefficienter med en sonderingspuls:

Wolfram Mathematica i geofysik

där s(t) — seismiskt spår, dvs. allt som spelades in med en hydrofon eller geofon under en fast inspelningstid, w(t) - signalen som genereras av luftpistolen, n(t) - slumpmässigt brus.

Låt oss beräkna ett syntetiskt seismiskt spår som ett exempel. Vi kommer att använda Ricker-pulsen, som ofta används i seismisk utforskning, som den initiala signalen.

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]

Initial seismisk impuls
Wolfram Mathematica i geofysik

Vi kommer att sätta två gränser på djupen 300 ms och 600 ms, och reflektionskoefficienterna kommer att vara slumptal

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 reflektionskoefficienter
Wolfram Mathematica i geofysik

Låt oss beräkna och visa det seismiska spåret. Eftersom reflektionskoefficienterna har olika tecken får vi två alternerande reflektioner på det seismiska spåret.

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]]

Simulerad bana
Wolfram Mathematica i geofysik

För det här exemplet är det nödvändigt att göra en reservation - i verkligheten bestäms djupet av lagren, naturligtvis, i meter, och beräkningen av det seismiska spåret sker för tidsdomänen. Det skulle vara mer korrekt att ställa in djupen i meter och beräkna ankomsttiderna med kännedom om hastigheterna i lagren. I det här fallet ställer jag genast in lagren på tidsaxeln.

Om vi ​​talar om fältforskning, som ett resultat av sådana observationer registreras ett stort antal liknande tidsserier (seismiska spår). Till exempel, när man studerar en plats 25 km lång och 15 km bred, där varje spår, som ett resultat av arbete, karakteriserar en cell som mäter 25x25 meter (en sådan cell kallas en bin), kommer den slutliga datamatrisen att innehålla 600000 1 spår. Med en samplingstid på 5 ms och en inspelningstid på 11 sekunder kommer den slutliga datafilen att vara mer än XNUMX ​​GB, och volymen av det ursprungliga "råmaterialet" kan vara hundratals gigabyte.

Hur man arbetar med dem Wolfram Mathematica?

paket GeologiIO

Utvecklingen av paketet började fråga på den rysktalande stödgruppens VK-vägg. Tack vare communityns svar hittade man en lösning mycket snabbt. Och som ett resultat växte det till en allvarlig utveckling. Motsvarande Wolfram Community väggstolpe Det markerades till och med av moderatorer. För närvarande stöder paketet arbete med följande datatyper som aktivt används inom den geologiska industrin:

  1. import av kartdata i ZMAP- och IRAP-format
  2. import av mätningar i brunnar i LAS-format
  3. in- och utmatning av seismiska filformat SEGY

För att installera paketet måste du följa instruktionerna på nedladdningssidan för det sammansatta paketet, d.v.s. exekvera följande kod i någon Mathematica anteckningsbok:

If[PacletInformation["GeologyIO"] === {}, PacletInstall[URLDownload[
    "https://wolfr.am/FiQ5oFih", 
    FileNameJoin[{CreateDirectory[], "GeologyIO-0.2.2.paclet"}]
]]]

Därefter kommer paketet att installeras i standardmappen, sökvägen till vilken kan erhållas enligt följande:

FileNameJoin[{$UserBasePacletsDirectory, "Repository"}]

Som ett exempel kommer vi att visa paketets huvudfunktioner. Anropet görs traditionellt för paket på Wolfram-språket:

Get["GeologyIO`"]

Paketet är utvecklat med hjälp av Wolfram arbetsbänk. Detta gör att du kan komplettera paketets huvudfunktionalitet med dokumentation, som när det gäller presentationsformat inte skiljer sig från dokumentationen för Wolfram Mathematica själv, och att förse paketet med testfiler för den första bekantskapen.

Wolfram Mathematica i geofysik

Wolfram Mathematica i geofysik

En sådan fil, i synnerhet, är filen "Marmousi.segy" - detta är en syntetisk modell av en geologisk sektion, som utvecklades av det franska petroleuminstitutet. Med hjälp av denna modell testar utvecklare sina egna algoritmer för vågfältsmodellering, databehandling, seismisk spårinversion, etc. Själva Marmousi-modellen lagras i arkivet där själva paketet laddades ner. Kör följande kod för att hämta filen:

If[Not[FileExistsQ["Marmousi.segy"]], 
URLDownload["https://wolfr.am/FiQGh7rk", "Marmousi.segy"];]
marmousi = SEGYImport["Marmousi.segy"]

Importera resultat - SEGYData-objekt
Wolfram Mathematica i geofysik

SEGY-formatet innebär att olika information om observationer lagras. För det första är det textkommentarer. Det gäller bland annat information om platsen för arbetet, namnen på de företag som utfört mätningarna m.m. I vårt fall anropas denna rubrik av en begäran med TextHeader-nyckeln. Här är en förkortad textrubrik:

Short[marmousi["TextHeader"]]

"Marmousi-datauppsättningen genererades vid institutet ... högsta hastighet på 1500 m/s och maximalt 5500 m/s)"

Du kan visa den faktiska geologiska modellen genom att komma åt de seismiska spåren med hjälp av "spår"-tangenten (en av funktionerna i paketet är att nycklarna är skiftlägesokänsliga):

ArrayPlot[Transpose[marmousi["traces"]], PlotTheme -> "Detailed"]

Modell Marmousi
Wolfram Mathematica i geofysik

För närvarande låter paketet dig även ladda data i delar från stora filer, vilket gör det möjligt att bearbeta filer vars storlek kan nå tiotals gigabyte. Paketets funktioner inkluderar även funktioner för att exportera data till .segy och delvis lägga till i slutet av filen.

Separat är det värt att notera paketets funktionalitet när man arbetar med den komplexa strukturen av .segy-filer. Eftersom det låter dig inte bara komma åt enskilda spår och rubriker med hjälp av nycklar och index, utan också att ändra dem och sedan skriva dem till en fil. Många av de tekniska detaljerna i GeologyIO:s implementering ligger utanför ramen för denna artikel och förtjänar förmodligen en separat beskrivning.

Relevans av spektralanalys vid seismisk utforskning

Möjligheten att importera seismisk data till Wolfram Mathematica gör att du kan använda inbyggd signalbehandlingsfunktion för experimentella data. Eftersom varje seismiskt spår representerar en tidsserie, är ett av huvudverktygen för att studera dem spektralanalys. Bland förutsättningarna för att analysera frekvenssammansättningen av seismiska data kan vi till exempel nämna följande:

  1. Olika typer av vågor kännetecknas av olika frekvenssammansättning. Detta låter dig markera användbara vågor och undertrycka interferensvågor.
  2. Bergegenskaper som porositet och mättnad kan påverka frekvenssammansättningen. Detta gör det möjligt att identifiera bergarter med de bästa egenskaperna.
  3. Skikt med olika tjocklek orsakar anomalier i olika frekvensområden.

Den tredje punkten är den viktigaste i denna artikels sammanhang. Nedan visas ett kodfragment för beräkning av seismiska spår i fallet med ett lager med varierande tjocklek - en kilmodell. Denna modell studeras traditionellt inom seismisk utforskning för att analysera interferenseffekter när vågor som reflekteras från många lager överlagras på varandra.

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-formation
Wolfram Mathematica i geofysik

Våghastigheten inuti kilen är 4500 m/s, utanför kilen 4000 m/s, och densiteten antas vara konstant 2200 g/cm³. För en sådan modell beräknar vi reflektionskoefficienter och seismiska spår.

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]

Seismiska spår för kilmodellen
Wolfram Mathematica i geofysik

Sekvensen av seismiska spår som visas i denna figur kallas en seismisk sektion. Som du kan se kan dess tolkning också utföras på en intuitiv nivå, eftersom geometrin hos de reflekterade vågorna tydligt motsvarar modellen som specificerades tidigare. Om du analyserar spåren mer detaljerat kommer du att märka att spår från 1 till cirka 30 inte skiljer sig - reflektionen från formationens tak och från botten överlappar inte varandra. Från och med det 31:a spåret börjar reflektionerna störa. Och även om reflektionskoefficienterna i modellen inte ändras horisontellt - de seismiska spåren ändrar sin intensitet när tjockleken på lagret ändras.

Låt oss överväga reflektionsamplituden från formationens övre gräns. Från och med den 60:e rutten börjar reflektionens intensitet öka och vid den 70:e rutten blir den maximal. Detta är hur störningen av vågor från taket och botten av lagren manifesterar sig, vilket i vissa fall leder till betydande anomalier i den seismiska posten.

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 över amplituden för den reflekterade vågen från kilens övre kant
Wolfram Mathematica i geofysik

Det är logiskt att när signalen är lägre frekvens börjar störningar uppträda vid stora formationstjocklekar, och i fallet med en högfrekvent signal uppstår störningar vid mindre tjocklekar. Följande kodsnutt skapar en signal med frekvenser på 35 Hz, 55 Hz och 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]

En uppsättning källsignaler med frekvenser på 35 Hz, 55Hz, 85Hz
Wolfram Mathematica i geofysik

Genom att beräkna seismiska spår och plotta grafer av reflekterade vågamplituder kan vi se att för olika frekvenser observeras en anomali vid olika formationstjocklekar.

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 över amplituderna för den reflekterade vågen från kilens övre kant för olika frekvenser
Wolfram Mathematica i geofysik

Förmågan att dra slutsatser om formationens tjocklek från resultaten av seismiska observationer är extremt användbar, eftersom en av huvuduppgifterna vid oljeprospektering är att bedöma de mest lovande punkterna för att lägga en brunn (dvs de områden där formationen är tjockare). Dessutom kan det i den geologiska sektionen finnas föremål vars tillkomst orsakar en kraftig förändring av formationens tjocklek. Detta gör spektralanalys till ett effektivt verktyg för att studera dem. I nästa del av artikeln kommer vi att överväga sådana geologiska objekt mer i detalj.

Experimentella data. Var fick du tag i dem och vad ska du leta efter i dem?

Materialet som analyserades i artikeln erhölls i västra Sibirien. Regionen är, som alla utan undantag säkert vet, den viktigaste oljeproducerande regionen i vårt land. Aktiv utveckling av fyndigheter började i regionen på 60-talet av förra seklet. Den huvudsakliga metoden för att söka efter oljefyndigheter är seismisk prospektering. Det är intressant att titta på satellitbilder av detta territorium. I liten skala kan du notera ett stort antal träsk och sjöar; genom att förstora kartan kan du se klusterbrunnars borrplatser, och genom att förstora kartan till det yttersta kan du också urskilja röjningarna av profilerna längs vilka seismiken observationer utfördes.

Satellitbild över Yandex-kartor - Noyabrsk stadsområde
Wolfram Mathematica i geofysik

Ett nätverk av brunnsplattor vid ett av fälten
Wolfram Mathematica i geofysik

Oljebärande stenar i västra Sibirien förekommer på ett brett spektrum av djup - från 1 km till 5 km. Huvudvolymen av stenar som innehåller olja bildades under jura och krita. Juratiden är nog känd för många från filmen med samma namn. Jurassiskt klimat skilde sig väsentligt från den moderna. Encyclopedia Britannica har en serie paleomap som kännetecknar varje helogisk era.

föreliggande
Wolfram Mathematica i geofysik
juraperioden
Wolfram Mathematica i geofysik

Observera att under juratiden var västra Sibiriens territorium en havskust (land som korsades av floder och ett grunt hav). Eftersom klimatet var behagligt kan vi anta att ett typiskt landskap på den tiden såg ut så här:

Jurassic Sibirien
Wolfram Mathematica i geofysik

På den här bilden är det viktiga för oss inte så mycket djuren och fåglarna, utan bilden av floden i bakgrunden. Floden är samma geologiska objekt som vi stannade vid tidigare. Faktum är att flodernas aktivitet tillåter att välsorterade sandstenar samlas, som sedan kommer att bli en reservoar för olja. Dessa reservoarer kan ha en bisarr, komplex form (som en flodbädd) och de har varierande tjocklek - nära bankerna är tjockleken liten, men närmare mitten av kanalen eller i meanderområden ökar den. Så floderna som bildas i Jurassic är nu på ett djup av cirka tre kilometer och är föremål för ett sökande efter oljereservoarer.

Experimentella data. Bearbetning och visualisering

Låt oss omedelbart göra en reservation beträffande de seismiska materialen som visas i artikeln - på grund av att mängden data som används för analysen är betydande - ingår bara ett fragment av den ursprungliga uppsättningen av seismiska spår i artikelns text. Detta kommer att tillåta vem som helst att reproducera ovanstående beräkningar.

När man arbetar med seismiska data använder en geofysiker vanligtvis specialiserad programvara (det finns flera branschledare vars utvecklingar används aktivt, till exempel Petrel eller Paradigm), som låter dig analysera olika typer av data och har ett bekvämt grafiskt gränssnitt. Trots all bekvämlighet har dessa typer av programvara också sina nackdelar - till exempel tar implementeringen av moderna algoritmer i stabila versioner mycket tid, och möjligheterna att automatisera beräkningar är vanligtvis begränsade. I en sådan situation blir det väldigt bekvämt att använda datormatematiksystem och högnivåprogrammeringsspråk, som tillåter användningen av en bred algoritmisk bas och som samtidigt tar på sig mycket rutin. Detta är principen som används för att arbeta med seismiska data i Wolfram Mathematica. Det är olämpligt att skriva rik funktionalitet för interaktivt arbete med data - det är viktigare att säkerställa laddning från ett allmänt accepterat format, tillämpa de önskade algoritmerna på dem och ladda upp dem tillbaka till ett externt format.

Efter det föreslagna schemat kommer vi att ladda de ursprungliga seismiska data och visa dem Wolfram Mathematica:

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]

Data som laddas ner och importeras på detta sätt är de rutter som registrerats på ett område som mäter 10 gånger 5 kilometer. Om data erhålls med hjälp av en tredimensionell seismisk undersökningsmetod (vågor registreras inte längs enskilda geofysiska profiler, utan över hela området samtidigt), blir det möjligt att erhålla seismiska datakuber. Dessa är tredimensionella objekt, vars vertikala och horisontella sektioner möjliggör en detaljerad studie av den geologiska miljön. I det övervägda exemplet har vi att göra med tredimensionell data. Vi kan få lite information från textrubriken, så här

StringPartition[seismic3DSEGY["textheader"], 80] // TableForm

C 1 DETTA ÄR DEMO-FIL FÖR GEOLOGYIO-PAKETEST
C 2
C 3
C 4
C 5 DATUM ANVÄNDARNAMN: WOLFRAM ANVÄNDARE
C 6 UNDERSÖKNINGENS NAMN: NÅGONSTANS I SIBERIEN
C 7 FILTYP 3D SEISMISK VOLYM
C 8
C 9
C10 Z Räckvidd: FÖRSTA 2200M SISTA 2400M

Denna datamängd kommer att räcka för att vi ska demonstrera huvudstadierna av dataanalys. Spåren i filen registreras sekventiellt och var och en av dem ser ut ungefär som följande figur - detta är fördelningen av amplituderna för reflekterade vågor längs den vertikala axeln (djupaxeln).

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]

Ett av de seismiska snittspåren
Wolfram Mathematica i geofysik

Genom att veta hur många spår som finns i varje riktning av det studerade området kan du generera en tredimensionell datamatris och visa den med hjälp av funktionen 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-bild av en seismisk datakub. (Vertikal axel - djup)
Wolfram Mathematica i geofysik

Om de geologiska egenskaperna av intresse skapar intensiva seismiska anomalier, kan visualiseringsverktyg med transparens användas. "Oviktiga" områden av inspelningen kan göras osynliga och lämnar endast avvikelser synliga. I Wolfram Mathematica kan detta göras med hjälp av Opacitet[] и Raster3D[].

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 datakubbild med funktionerna Opacity[] och Raster3D[] Wolfram Mathematica i geofysik

Som i det syntetiska exemplet kan man på delar av den ursprungliga kuben identifiera några geologiska gränser (lager) med variabel relief.

Det huvudsakliga verktyget för spektralanalys är Fouriertransformen. Med dess hjälp kan du utvärdera amplitud-frekvensspektrumet för varje spår eller grupp av spår. Men efter att ha överfört data till frekvensdomänen går information förlorad om vid vilka tidpunkter (läs på vilka djup) frekvensen ändras. För att kunna lokalisera signalförändringar på tids-(djup)axeln används den fönsterförsedda Fouriertransformen och wavelet-sönderdelningen. Den här artikeln använder wavelet-nedbrytning. Wavelet-analysteknik började användas aktivt i seismisk utforskning på 90-talet. Fördelen jämfört med den fönsterförsedda Fouriertransformen anses vara bättre tidsupplösning.

Med hjälp av följande kodfragment kan du dekomponera ett av de seismiska spåren i enskilda 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]

Nedbrytning av ett spår till komponenter
Wolfram Mathematica i geofysik

För att bedöma hur reflektionsenergin är fördelad vid olika vågankomsttider används skalogram (analogt med ett spektrogram). Som regel finns det i praktiken inget behov av att analysera alla komponenter. Normalt väljs låg-, mellan- och högfrekvenskomponenter.

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. Funktionsresultat WaveletScalogram[]
Wolfram Mathematica i geofysik

Wolfram-språket använder funktionen för wavelet-transformation ContinuousWaveletTransform[]. Och tillämpningen av denna funktion på hela uppsättningen spår kommer att utföras med funktionen Tabell[]. Här är det värt att notera en av styrkorna med Wolfram Mathematica - förmågan att använda parallellisering Parallelltabell[]. I exemplet ovan finns det inget behov av parallellisering - datamängden är inte stor, men när man arbetar med experimentella datamängder som innehåller hundratusentals spår är detta en nödvändighet.

tracesCWD=Table[Map[Hilbert[#,0]&,Re[ContinuousWaveletTransform[traces[[i]]][[1]]][[{13,15,18}]]],{i,1,Length@traces}]; 

Efter att ha tillämpat funktionen ContinuousWaveletTransform[] Nya datamängder visas som motsvarar de valda frekvenserna. I exemplet ovan är dessa frekvenser: 38Hz, 33Hz, 27Hz. Valet av frekvenser utförs oftast på basis av tester - de får effektiva kartor för olika frekvenskombinationer och väljer den mest informativa från en geologs synvinkel.

Om du behöver dela resultaten med kollegor eller tillhandahålla dem till kunden kan du använda funktionen SEGYExport[] i GeologyIO-paketet

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 dessa kuber (lågfrekventa, mellanfrekventa och högfrekventa komponenter) används vanligtvis RGB-blandning för att visualisera data tillsammans. Varje komponent tilldelas sin egen färg - röd, grön, blå. I Wolfram Mathematica kan detta göras med hjälp av funktionen ColorCombine[].

Resultatet är bilder från vilka geologisk tolkning kan göras. De meandrar som registreras på sektionen gör det möjligt att avgränsa paleokanaler, som mer sannolikt är reservoarer och innehåller oljereserver. Sökningen och analysen av moderna analoger av ett sådant flodsystem gör att vi kan bestämma de mest lovande delarna av slingrarna. Kanalerna i sig kännetecknas av tjocka lager av välsorterad sandsten och är en bra reservoar för olja. Områden utanför "spets"-avvikelserna liknar moderna översvämningsavlagringar. Översvämningsavlagringar representeras huvudsakligen av leriga bergarter och borrning i dessa zoner kommer att vara ineffektiv.

RGB-del av datakuben. I mitten (något till vänster om mitten) kan du spåra den slingrande floden.
Wolfram Mathematica i geofysik
RGB-del av datakuben. På vänster sida kan du spåra den slingrande floden.
Wolfram Mathematica i geofysik

I vissa fall möjliggör kvaliteten på seismiska data betydligt tydligare bilder. Detta beror på fältarbetsmetodik, utrustningen som används av brusreduceringsalgoritmen. I sådana fall är inte bara fragment av flodsystem synliga, utan även hela utsträckta paleo-floder.

RGB-blandning av tre komponenter i en seismisk datakub (horisontell skiva). Djup ca 2 km.
Wolfram Mathematica i geofysik
Satellitbild av floden Volga nära Saratov
Wolfram Mathematica i geofysik

Slutsats

Wolfram Mathematica låter dig analysera seismiska data och lösa tillämpade problem relaterade till mineralprospektering, och GeologyIO-paketet gör denna process mer bekväm. Strukturen av seismiska data är sådan att man använder inbyggda metoder för att påskynda beräkningar (Parallelltabell[], ParallelDo[],...) är mycket effektiv och låter dig bearbeta stora mängder data. Till stor del underlättas detta av datalagringsfunktionerna i GeologyIO-paketet. Förresten, paketet kan användas inte bara inom området för tillämpad seismisk utforskning. Nästan samma typer av data används i markpenetrerande radar och seismologi. Om du har förslag på hur du kan förbättra resultatet, vilka signalanalysalgoritmer från Wolfram Mathematica-arsenalen som är tillämpliga på sådan data, eller om du har några kritiska kommentarer, vänligen Lämna en kommentar.

Källa: will.com

Lägg en kommentar