Дзякуем аўтару блога за яго даклад
Увядзенне
Гэта нататка напісана па слядах канферэнцыі і змяшчае канспект даклада, з якім я выступаў. Мерапрыемства адбылося ў чэрвені ў горадзе Санкт-Пецярбургу. Улічваючы тое, што працую я ў квартале ад месца правядзення канферэнцыі, я не мог не наведаць гэтую падзею. У 2016 і 2017 гадах я слухаў даклады канферэнцыі, а сёлета выступіў з дакладам. Во–первых, з'явілася цікавая (як мне здаецца) тэма, якую мы развіваем з , а па-другое, пасля працяглага вывучэння заканадаўства РФ у частцы санкцыйнай палітыкі, на прадпрыемстве, дзе я працую з'явілася аж цэлых дзве ліцэнзіі .
Перш, чым перайсці да тэмы майго выступу, я б хацеў адзначыць добрую арганізацыю мерапрыемства. На візітнай старонцы канферэнцыі выкарыстоўваецца выява Казанскага сабора. Сабор з'яўляецца адной з галоўных славутасцяў Пецярбурга і яго вельмі добра відаць з залы, у якой праходзіла канферэнцыя.
На ўваходзе ў СПбГЭУ удзельнікаў сустракалі памочнікі з ліку студэнтаў - не дазвалялі заблудзіцца. Падчас рэгістрацыі выдаваліся невялікія сувенірчыкі (цацка – які мігціць знітоўкі, ручка, налепкі з сімволікай Wolfram). Абед і кава-брэйк таксама былі ўключаны ў расклад канферэнцыі. Пра смачную каву і піражкі я ўжо адзначаў на сцяне групы – кухары малайцы. Гэтай уступнай часткай я б хацеў падкрэсліць, што само мерапрыемства, яго фармат і месца правядзення ўжо прыносяць станоўчыя эмоцыі.
Даклад, які быў падрыхтаваны мною і Кірылам Бяловым называецца «Выкарыстанне Wolfram Mathematica для рашэння задач прыкладной геафізікі. Спектральны аналіз сейсмічных дадзеных ці "куды беглі старажытныя рэкі". Змест даклада пакрывае дзве часткі: па-першае, гэта выкарыстанне алгарытмаў, якія ёсць у для аналізу геафізічных дадзеных, а па-другое, гэта тое як геафізічныя дадзеныя змясціць у Wolfram Mathematica.
Сейсмаразведка
Для пачатку трэба зрабіць невялікі экскурс у геафізіку. Геафізіка - гэта навука, якая вывучае фізічныя ўласцівасці горных парод. Ну а паколькі пароды маюць розныя ўласцівасці: электрычныя, магнітныя, пругкія, то існуюць адпаведныя метады геафізікі: электраразведка, магнітаразведка, сейсмаразведка… У кантэксце дадзенага артыкула падрабязней закранем толькі сейсмічную разведку. Сейсмаразведка з'яўляецца асноўным метадам пошуку нафты і газу. Метад заснаваны на ўзрушанасці пругкіх ваганняў і наступнай рэгістрацыі водгуку ад горных парод складалых доследную тэрыторыю. Узрушанасць ваганняў ажыццяўляецца на сушы (дынамітам або невзрывными вібрацыйнымі крыніцамі пругкіх ваганняў) або на моры (пнеўмапушкамі). Пругкія ваганні распаўсюджваюцца праз тоўшчу горных парод пераламляючыся і адлюстроўваючыся на межах пластоў з рознымі ўласцівасцямі. Адлюстраваныя хвалі вяртаюцца на паверхню і рэгіструюцца геафонамі на сушы (як правіла гэта электрадынамічныя прыборы, заснаваныя на руху магніта, падвешанага ў шпульцы) або гідрафонамі ў моры (заснаваны на п'езаэфекце). Па часе прыходу хваль можна меркаваць аб глыбінях геалагічных пластоў.
Сейсмічнае судна буксіруе абсталяванне
Пнеўмапушка ўзбуджае пругкія ваганні
Хвалі праходзяць праз тоўшчу горных парод і рэгіструюцца гідрафонамі
Навукова-даследчае судна геафізічнай разведкі "Іван Губкін" на прычале каля Благавешчанскага моста ў Пецярбургу
Мадэль сейсмічнага сігналу
Горныя пароды маюць розныя фізічныя ўласцівасці. Для сейсмавыведкі перш за ўсё важныя пругкія ўласцівасці - хуткасць распаўсюджвання пругкіх ваганняў і шчыльнасць. Калі два пласта будуць мець аднолькавыя ці блізкія ўласцівасці то хваля "не заўважыць" мяжу паміж імі. Калі ж хуткасці хваль у пластах будуць адрознівацца, то на мяжы пластоў паўстане адлюстраванне. Чым больш розніца ва ўласцівасцях тым больш інтэнсіўнае адлюстраванне. Яго інтэнсіўнасць будзе вызначацца каэфіцыентам адлюстравання (rc):
дзе ρ - шчыльнасць парод, ν - хуткасць хваль, 1 і 2 абазначаюць верхні і ніжні пласты.
Адной з найболей простых і часта выкарыстоўваных мадэляў сейсмічнага сігналу з'яўляецца скруткавая мадэль, калі зарэгістраваная сейсмічная траса ўяўляецца, як вынік скруткі паслядоўнасці каэфіцыентаў адлюстравання з зандаваным імпульсам:
дзе s(t) - сейсмічная траса, г.зн. усё што запісаў гідрафон ці геафон на працягу фіксаванага часу рэгістрацыі, w(t) - Сігнал, які генеруе пнеўмапушка, n(t) - выпадковы шум.
Разлічым для прыкладу сінтэтычную сейглядасу. У якасці зыходнага сігналу будзем выкарыстоўваць шырока выкарыстоўваны ў сейсмавыведцы імпульс Рыкера.
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]Зыходны сейсмічны імпульс
Дзве межы зададзім на глыбінях 300 мс і 600 мс, а каэфіцыенты адлюстравання будуць выпадковымі лікамі.
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]]Паслядоўнасць каэфіцыентаў адлюстравання
Разлічым і адлюструем сейсмічную трасу. Паколькі каэфіцыенты адлюстравання маюць розныя знакі, то і на сейсмічнай трасе атрымліваем два знакапераменных адлюстравання.
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]]Змадэляваная траса
Для дадзенага прыкладу трэба зрабіць агаворку - у рэчаіснасці глыбіня пластоў вызначаецца вядома ў метрах, а разлік сейсмічнай трасы адбываецца для часовай вобласці. Правільней было б задаць глыбіні ў метрах і разлічыць часы прыходу ведаючы хуткасці ў пластах. У дадзеным выпадку я адразу задаў пласты на часовай восі.
Калі казаць аб палявых даследаваннях, то ў выніку такіх назіранняў рэгіструецца велізарная колькасць падобных часавых шэрагаў (сейсмічных трас). Напрыклад, пры даследаванні ўчастку даўжынёй 25 км і шырынёй 15 км, дзе ў выніку прац кожная траса характарызуе вочка памерам 25х25 метраў (такое вочка завецца бін), фінальны масіў дадзеных будзе ўтрымоўваць 600000 трас. Пры кроку дыскрэтызацыі па часе роўным 1 мс, часу запісу 5 секунд канчатковы файл дадзеных складзе больш за 11 Гб, а аб'ём зыходнага "волкага" матэрыялу можа скласці сотні гігабайт.
Як працаваць з імі ў ?
Пакет
Пачаткам распрацоўкі пакета стаў на сцяне VK гурта рускамоўнай падтрымкі. Дзякуючы адказам суполкі рашэнне было знойдзена вельмі хутка. І ў выніку перарасло ў сур'ёзную распрацоўку. Адпаведны нават быў адзначаны мадэратарамі. На дадзены момант пакет падтрымлівае працу з наступнымі тыпамі дадзеных, якія актыўна выкарыстоўваюцца ў геалагічнай галіне:
- імпарт картаграфічных даных фармату ZMAP і IRAP
- імпарт вымярэнняў у свідравінах фармату LAS
- увод і вывад сейсмічных файлаў фармату
Каб усталяваць пакет неабходна прытрымлівацца інструкцый на старонцы загрузкі сабранага пакета, г.зн. выканаць наступны код у любым :
If[PacletInformation["GeologyIO"] === {}, PacletInstall[URLDownload[
"https://wolfr.am/FiQ5oFih",
FileNameJoin[{CreateDirectory[], "GeologyIO-0.2.2.paclet"}]
]]]Пасля чаго пакет усталюецца ў тэчку па змаўчанні, шлях да якой можна атрымаць наступным чынам:
FileNameJoin[{$UserBasePacletsDirectory, "Repository"}]Для прыкладу прадэманструем асноўныя магчымасці пакета. Выклік ажыццяўляецца традыцыйна для пакетаў у Wolfram Language:
Get["GeologyIO`"]Пакет распрацоўваецца з выкарыстаннем . Гэта дазваляе суправадзіць асноўны функцыянал пакета дакументацыяй, якая па фармаце падання не адрозніваецца ад дакументацыі самой Wolfram Mathematica і забяспечыць пакет тэставымі файламі для першага знаёмства.
Такім файлам, у прыватнасці з'яўляецца файл "Marmousi.segy" - гэта сінтэтычная мадэль геалагічнага разрэзу, якую распрацаваў французскі інстытут нафты. Выкарыстоўваючы гэтую мадэль распрацоўшчыкі тэстую ўласныя алгарытмы мадэлявання хвалевага поля, апрацоўкі дадзеных, інверсіі сейсмічных трас і г.д. Сама мадэль Marmousi захоўваецца ў рэпазітары, адкуль быў запампаваны сам пакет. Для таго, каб атрымаць файл - выканаем наступны код:
If[Not[FileExistsQ["Marmousi.segy"]],
URLDownload["https://wolfr.am/FiQGh7rk", "Marmousi.segy"];]
marmousi = SEGYImport["Marmousi.segy"]Вынік імпарту - аб'ект SEGYData
Фармат SEGY мяркуе захоўванне рознай інфармацыі аб назіраннях. Па-першае, гэта тэкставыя каментарыі. Сюды заносяць інфармацыю аб месцы правядзення работ, назвы кампаній, якія выканалі вымярэнні і г.д. У нашым выпадку гэты загаловак выклікаецца запытам з ключом TextHeader. Тут прыведзены скарочаны тэкставы загаловак:
Short[marmousi["TextHeader"]]"The Marmousi data set была generated ў Institute ...nimum velocity of 1500 m/s and maximum of 5500 m/s)"
Адлюстраваць уласна гелагічную мадэль можна звярнуўшыся да сейсмічных трас па ключы «traces» (адной з фіч пакета з'яўляецца незалежнасць ключоў ад рэгістра):
ArrayPlot[Transpose[marmousi["traces"]], PlotTheme -> "Detailed"]Мадэль Marmousi
На бягучы момант пакет таксама дазваляе загружаць дадзеныя часткамі з вялікіх файлаў, дзякуючы чаму становіцца магчыма апрацоўка файлаў, памер якіх можа дасягаць дзясяткаў гігабайт. Таксама ў лік функцый пакета ўваходзяць функцыі для экспарту дадзеных у .segy і частковай дазапісы ў канец файла.
Асобна варта адзначыць функцыянальнасць пакета пры працы са складанай структурай. Segy-файлаў. Бо ён дазваляе не толькі звяртацца па ключах і азначнікам да асобных трас, загалоўкам, але і змяняць іх з наступным запісам у файл. Многія тэхнічныя дэталі рэалізацыі GeologyIO выходзяць за рамкі гэтага артыкула і, верагодна, заслугоўваюць асобнага апісання.
Актуальнасць спектральнага аналізу ў сейсмаразведцы
Магчымасць імпарту сейсмічных матэрыялаў у Wolfram Mathematica дазваляе выкарыстоўваць убудаваны функцыянал апрацоўкі сігналаў для эксперыментальных дадзеных. Паколькі кожная сейсмічная траса ўяўляе сабою часавы шэраг, то адным з асноўных прылад іх вывучэння з'яўляецца спектральны аналіз. Сярод перадумоў аналізу частотнага складу сейсмічных дадзеных можна назваць, напрыклад, наступныя:
- Розныя тыпы хваль характарызуюцца розным частотным складам. Гэта дазваляе вылучыць карысныя хвалі і здушыць хвалі-перашкоды.
- Такія ўласцівасці горных парод, як сітаватасць і насычанасць могуць уплываць на частотны склад. Гэта дазваляе выдзяляць горныя пароды з лепшымі ўласцівасцямі.
- Пласты з рознай таўшчынёй абумоўліваюць анамаліі ў розных частотных дыяпазонах.
Трэці пункт з'яўляецца асноўным у кантэксце гэтага артыкула. Ніжэй прыведзены фрагмент кода для разліку сейсмічных трас у выпадку пласта з якая змяняецца таўшчынёй - мадэль кліну. Гэтая мадэль традыцыйна вывучаецца ў сейсмавыведцы для анлізу інтэрферэнцыйных эфектаў, калі хвалі адлюстраваныя ад шматлікіх пластоў накладваюцца сябар на сябра.
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]]Мадэль выкліноўваецца пласта
Хуткасць хваль усярэдзіне кліну складае 4500 м/з, па-за клінам 4000м/з, а шчыльнасць прынятая сталай 2200 г/гл³. Для такой мадэлі разлічым каэфіцыенты адлюстравання і сейсмічныя трасы.
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]Сейсмічныя трасы для мадэлі кліну
Паслядоўнасць сейсмічных трас намаляваная на гэтым малюнку завецца сейсмічным разрэзам. Як можна заўважыць, яго інтэрпрэтацыя можа выконвацца і на інтуітыўным узроўні, паколькі геаметрыя адлюстраваных хваль адназначна адпавядае мадэлі, якая была зададзена раней. Калі больш дэталёва аналізаваць трасы, то можна заўважыць, што трасы з 1-й па, прыкладна, 30-ю не адрозніваюцца - адлюстраванне ад даху пласта і ад падэшвы не накладваюцца сябар на сябра. Пачынальна з 31-й трасы адлюстравання пачынаюць інтэрфераваць. І, хоць, у мадэлі каэфіцыенты адлюстравання не змяняюцца па гарызанталі - сейсмічныя трасы змяняюць сваю інтэнсіўнасць пры змене таўшчыні пласта.
Разгледзім амплітуду адлюстравання ад верхняй мяжы пласта. Пачынальна з 60-й трасы інтэнсіўнасць адлюстравання пачынае ўзрастаць і на 70-й трасе становіцца максімальнай. Так выяўляецца інтэрферэнцыя хваль ад даху і падэшвы для пластоў, прыводзячы ў некаторых выпадках да істотных анамаліяў сейсмічнага запісу.
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]Графік амплітуды адлюстраванай хвалі ад верхняй абзы кліну
Лагічна, што калі сігнал больш нізкачашчынны, то інтэрферэнцыя пачынае выяўляцца пры вялікіх таўшчынях пласта, а ў выпадку высокачашчыннага сігналу інтэрфрэнцыя ўзнікае пры малодшых таўшчынях. Наступны фрагмент кода стварае сігнал з частотамі 35 Гц, 55 Гц і 85 Гц.
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]Набор зыходных сігналаў з частотамі 35 Гц, 55Гц, 85Гц
Выканаўшы разлік сейсмічных трас і пабудаваўшы графікі амплітуд адлюстраванай хвалі, мы можам убачыць, што для розных частот анамалія назіраецца пры розных таўшчынях пласта.
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]Графікі амплітуд адлюстраванай хвалі ад верхняй абзы кліну для розных частот
Магчымасць рабіць высновы аб таўшчыні пласта па выніках сейсмічных назіранняў з'яўляецца вельмі карыснай, бо адной з галоўных задач пры разведцы радовішча нафты з'яўляецца ацэнка найбольш перспектыўных кропак для залажэння свідравіны (г.зн. тых участкаў, дзе пласт мае вялікую таўшчыню). Акрамя гэтага ў геалагічным разрэзе могуць сустракацца такія аб'екты, якія сваім генезісам абумоўліваюць рэзкую змену таўшчынь пласта. Гэта робіць спектральны аналізу эфектыўным інструментам іх вывучэння. У наступнай частцы артыкула разгледзім такія гелагічныя аб'екты падрабязней.
Эксперыментальныя даныя. Дзе атрыманы і што ў іх шукаць?
Матэрыялы, якія аналізуюцца ў артыкуле атрыманы на тэрыторыі Заходняй Сібіры. Рэгіён, як напэўна ведаюць усе без выключэння, з'яўляецца асноўным нафтаздабыўным рэгіёнам нашай краіны. Актыўная распрацоўка метараджэнняў пачалася ў рэгіёне ў 60-я гады мінулага стагоддзя. Асноўным метадам пошуку радовішчаў нафты з'яўляецца сейсмаразведка. Цікава разглядаць спадарожнікавыя здымкі гэтай тэрыторыі. Пры малым маштабе можна адзначыць велізарную колькасць балот і азёр, павялічваючы карту можна ўбачыць куставыя пляцоўкі свідравання свідравін, а павялічыўшы карту да мяжы можна адрозніць і прасекі профіляў, па якіх выконваліся сейсмічныя назіранні.
Спадарожнікавы здымак Яндэкс карт — раён горада Лістапад.
Сетка кустовых пляцовак на адным з метараджэнняў
Нафтаносныя пароды Запаной Сібіры залягаюць у шырокім дыяпазоне глыбінь - ад 1км да 5км. Асноўны аб'ём парод змяшчаюць нафту сфарміраваны ў юрскі і крэйдавы час. Юрскі перыяд напэўна шмат каму вядомы па аднайменным фільме. істотна адрозніваўся ад сучаснага. У энцыклапедыі Брытаніка ёсць серыя палеакарт, якія характарызуюць кожную гелагічную эпоху.
Цяперашні час
Юрскі перыяд
Звярніце ўвагу, што ў юрскі час тэрыторыя Заходняй Сібіры ўяўляла сабой марское ўзбярэжжа (суша перасякальная рэкамі і плыткае мора). Паколькі клімат быў камфортным, можна меркаваць, што тыповы пейзаж таго часу выглядаў наступным чынам:
Сібір юрскага перыяду
На гэтым малюнку для нас важныя не гэтулькі звяры і птушкі, колькі выява ракі на заднім плане. Рака - гэта той самы гелагічны аб'ект, на якім мы спыняліся раней. Справа ў тым, што дзейнасць рэк дазваляе назапашвацца добра сартаваным пяшчанікам, якія затым стануць рэзервуарам для нафты. Гэтыя рэзервуары могуць мець мудрагелістую, складаную форму (як і рэчышча ракі) і яны маюць зменлівую таўшчыню - у берагоў таўшчыня малая, а бліжэй да цэнтра рэчышча або на ўчастках меандр узрастае. Такім чынам, сфармаваныя ў юрскі час рэкі зараз знаходзяцца на глыбіні каля трох кіламетраў і з'яўляюцца аб'ектам пошуку рэзервуараў нафты.
Эксперыментальныя даныя. Апрацоўка і візуалізацыя
Зробім адразу агаворку, адносна сейсмічных матэрыялаў паказаных у артыкуле - у сілу таго, што аб'ём дадзеных выкарыстаных для аналізу значны - у тэкст артыкула змешчаны толькі фрагмент арыгінальнага набору сейсмічных трас. Гэта дазволіць усім жадаючым прайграць прыведзеныя вылічэнні.
Працуючы з сейсмічнымі дадзенымі геафізік, як правіла выкарыстоўвае спецыялізаванае праграмнае забеспячэнне (ёсць некалькі перадавікоў галіны, чыімі распрацоўкамі актыўна карыстаюцца, напрыклад Petrel або Paradigm), якое дазваляе аналізаваць розныя тыпы дадзеных і мае зручны графічны інтэрўфейс. Нягледзячы на ўсю зручнасць такія віды ПА маюць і свае недахопы - напрыклад укараненне сучасных алгарытмаў у стабільныя версіі займае шмат часу, а магчымасці аўтаматызацыі разлікаў, як правіла, абмежаваныя. У такой сітуацыі вельмі зручным становіцца выкарыстанне сістэм кампутарнай матэматыкі і моў праграмавання высокага ўзроўня, якія дазваляюць выкарыстоўваць шырокую алгарытмічную базу і, разам з тым, бяруць на сябе шмат руціны. Па такім прынцыпе і пабудавана праца з сейсмічнымі дадзенымі ў 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]Загружаныя і імпартаваныя такім чынам даныя гэта трасы зарэгістраваныя на ўчастку памерам 10 на 5 кіламетраў. У тым выпадку, калі дадзеныя атрыманы па методыцы трохмернай сейсмавыведкі (рэгістрацыя хваль вядзецца не ўздоўж асобных геафізічных профіляў, а на ўсім пляцы адначасова) становіцца магчымым атрымаць кубы сейсмічных дадзеных. Гэта трохмерныя аб'екты, вертыкальныя і гарызантальныя зрэзы якіх дазваляюць дэталёва вывучыць геалагічнае асяроддзе. У разгледжаным прыкладе мы маем справу як раз трохмернымі дадзенымі. Некаторыя звесткі мы можам атрымаць з тэкставага загалоўка, напрыклад так
StringPartition[seismic3DSEGY["textheader"], 80] // TableFormC 1 THIS IS DEMO FILE FOR GEOLOGYIO PACKAGE TEST
C 2
C 3
C 4
C 5 DATE USER NAME: WOLFRAM USER
C 6 SURVEY NAME: SOMEWHERE IN SIBERIA
C 7 FILE TYPE 3D SEISMIC VOLUME
C 8
C 9
C10 Z RANGE: FIRST 2200M LAST 2400M
Гэтага набору дадзеных нам будзе дастаткова, каб прадэманстраваць асноўныя этапы аналізу даных. Трасы ў файле запісаны паслядоўна і кожная з іх выглядае прыкладна як на наступным малюнку -гэтае размеркаванне амплітуд адлюстраваных хваль уздоўж вертыкальнай восі (восі глыбінь).
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]Адна з трас сейсмічнага разрэзу
Ведаючы тое, якая колькасць трас размешчана ў кожным напрамку вывучанага ўчастка можна сфарміраваць трохмерны масіў дадзеных і адлюстраваць яго з дапамогай функцыі 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]]Трохмернае выява куба сейсмічных дадзеных. (Вертыкальная вось - глыбіня)
У тым выпадку калі геалагічныя аб'екты, уяўлялыя цікавасць, ствараюць інтэнсіўныя сейсмічныя анамаліі, то можна выкарыстоўваць прылады візуалізацыі з празрыстасцю. «Няважныя» участкі запісу можна зрабіць нябачнымі, пакідаючы бачнымі толькі анамаліі. У Wolfram Mathematica гэта можна зрабіць з дапамогай и .
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]Выява куба сейсмічных дадзеных з выкарыстаннем функцый Opacity[] і Raster3D[]
Як і на сінтэтычным прыкладзе, на зрэзах арыгінальнага куба можна вылучыць некаторыя геалагічныя межы (пласты) з зменлівым рэльефам.
Асноўнай прыладай спектральнага аналізу з'яўляецца пераўтварэнне Фур'е. З яго дапамогай можна ацаніць амплітудна-частотны спектр кожнай трасы ці групы трас. Аднак, пасля перакладу дадзеных у вобласць частот губляецца інфармацыя аб тым, на якіх часах (чытай на якіх глыбінях) змяняецца частата. Для таго, каб мець магчымасць лакалізаваць змены сігналу на часовай (глыбіннай) восі выкарыстоўваюць аконнае пераўтварэнне Фур'е і вэйвлет разлажэнне. У дадзеным артыкуле выкарыстоўваецца вэйвлет разлажэнне. Тэхналогія вейвлет аналізу стала актыўна прымяняцца ў сейсмаразведцы ў 90-х гадах. Перавагай перад аконным пераўтварэннем Фур'е лічыцца лепшая часовая дазволенасць.
З дапамогай наступнага фрагмента кода можна выканаць раскладанне на асобныя кампаненты адной з сейсмічнах трас:
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]Раскладанне трасы на кампаненты
Для ацэнкі таго, як размеркавана энергія адлюстравання на розных часах прыходу хваляў выкарыстоўваюцца скалаграмы (аналаг спектраграмы). Як правіла на практыцы няма неабходнасці аналізаваць усе кампаненты. Звычайна выбіраюць нізка-, сярэдне- і высокачашчынны складнік.
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]Скалаграма. Вынік функцыі
У Wolfram Language для вэйвлет пераўтварэнні выкарыстоўваецца функцыя . А прымяненне гэтай функцыі да ўсяго набору трас адбудзецца з выкарыстаннем функцыі . Тут варта адзначыць адну з моцных бакоў Wolfram Mathematica магчымасць выкарыстоўваць разпраралельванне . У прыведзеным прыкладзе ў распаралельванні няма неабходнасці - аб'ём дадзеных не вялікі, але пры працы з эксперыментальнымі наборамі дадзеных якія змяшчаюць сотні тысяч трас гэта з'яўляецца неабходнасцю.
tracesCWD=Table[Map[Hilbert[#,0]&,Re[ContinuousWaveletTransform[traces[[i]]][[1]]][[{13,15,18}]]],{i,1,Length@traces}]; Пасля прымянення функцыі з'яўляюцца новыя масівы дадзеных якія адпавядаюць абраным частотам. У прыведзеным вышэй прыкладзе гэта частоты: 38Гц, 33Гц, 27Гц. Выбар частот ажыццяўляецца часцей за ўсё на аснове тэставання -атрымліваюць дзейсныя карты для розных частотных камбінацый і выбіраюць найболей інфарматыўны з пункта гледжання геолага.
Калі вынікамі неабходна падзяліцца з калегамі ці падаць іх замоўцу, то можна выкарыстоўваць функцыю SEGYExport[] пакета 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];Маючы ў распараджэнні тры такіх куба (нізкачашчынную, сярэднечашчынную і высокачашчынную кампаненты), як правіла, выкарыстаюць RGB змешванне для сумеснай візуалізацыі дадзеных. Кожнай з кампанент прысвойваецца свой колер - red, green, blue. У Wolfram Mathematica гэта можна зрабіць з выкарыстаннем функцыі .
У выніку атрымліваюцца выявы, па якіх можна выконваць гелагічную інтэрпрэтацыю. Меандры, якія фісуюцца на зрэзе дазваляюць сконтурить палеаруслы, якія з большай верагоднасцю могуць быць рэзервуарамі і ўтрымоўваць запасы нафты. Пошук і аналіз сучасных аналагаў такой рачной сістэмы дазваляе вызначыць найбольш перспектыўныя часткі меандр. Уласна рэчышча характарызуюцца магутнымі пластамі добра сартаванага пяшчаніку і з'яўляюцца добрым рэзервуарам для нафты. Участкі за межамі «шнуркавых» анамалій аналагічныя сучасным пойменным адкладам. Пойменныя адклады ў асноўным прадстаўлены гліністымі пародамі і бурэнне ў гэтыя зоны будзе неэфектыўным.
RGB зрэз куба дадзеных. У цэнтры (некалькі лявей за цэнтр) можна прасачыць раку, якая меандруе.
RGB зрэз куба дадзеных. У левай частцы можна прасачыць меандруючую раку.
У некаторых выпадках якасць сейсмічных дадзеных дазваляе атрымліваць істотна больш выразныя выявы. Гэта залежыць ад методыкі палявых работ, абсталявання, прымяюцца алгарытам у шумапрыглушэнні. У такіх выпадках бачныя не толькі фрагменты рачных сістэмы, але і цэлыя працяглыя палеарэкі.
RGB змешванне трох кампанент куба сейсмічных дадзеных (гарызантальны зрэз). Глыбіня прыкладна 2км.
Выява са спадарожніка ракі Волгі ў раёне Саратава
Заключэнне
У Wolfram Mathematica можна аналізаваць сейсмічныя дадзеныя і вырашаць прыкладныя задачы злучаныя з пошукам карысных выкапняў, а пакет GeologyIO дазваляе зрабіць гэты працэс зручнейшым. Структура сейсмічных дадзеных такая, што выкарыстанне убудаваных спосабаў паскарэння разлікаў (, ,…) з'яўляецца вельмі эфектыўным і дазваляе апрацоўваць вялікія аб'ёмы дадзеных. У немалой ступені гэтаму спрыяюць асаблівасці захоўвання дадзеных пакета GeologyIO. Дарэчы сказаць, пакет можа выкарыстоўвацца не толькі ў вобласці прыкладной сейсмавыведкі. Практычна такія ж тыпы дадзеных выкарыстоўваюцца ў геарадыёлакацыі і сейсмалогіі. Калі ў вас ёсць прапановы як палепшыць вынік, якія алгарытмы аналізу сігналу з арсенала Wolfram Mathematica дастасавальныя да такіх дадзеных ці ў вас маюцца крытычныя заўвагі - пакідайце каментары.
Крыніца: habr.com