Благодарим автора на блога за неговата презентация
Въведение
Тази забележка е написана в следствие на конференцията и съдържа резюме на доклада, с който аз се изказах. Събитието се проведе през юни в град Санкт Петербург. Като се има предвид, че работя само на квартал от мястото на провеждане на конференцията, не можех да не посетя това събитие. През 2016 и 2017 година слушах докладите на конференцията, а тази година изнесох доклад. На първо място, стана интересна (както ми се струва) тема, която развиваме с , а на второ място, след дълго изучаване на законодателството на Руската Федерация относно санкционната политика, в предприятието, в което работя, се появиха две лицензии .
Преди да премина към темата на моето изказване, бих искал да подчертая добрата организация на събитието. На визитната страница на конференцията се използва изображение на Казанския събор. Съборът е една от основните атракции на Петербург и е много добре видим от залата, в която се проведе конференцията.
На входа на СПбГЭУ участниците бяха посрещнати от помощници между студентите – не позволяваха да се изгубим. По време на регистрацията бяха раздавани малки сувенирчета (играчка – мигащ спайки, химикалка, стикери с символиката на Wolfram). Обядът и кафето също бяха включени в програмата на конференцията. За вкусното кафе и сладкиши вече споменах на стената на групата – готвачите се справиха отлично. С тази уводна част бих искал да подчертая, че самото събитие, неговият формат и мястото на провеждане вече предизвикаха положителни емоции.
Докладът, който беше подготвен от мен и Кирил Белов, се нарича «Използване на Wolfram Mathematica за решаване на задачи от приложна геофизика. Спектрален анализ на сеизмологични данни или «къде тичаха древните реки». Съдържанието на доклада обхваща две части: на първо място, това е използване на алгоритми, налични в за анализ на геофизични данни, а на второ място, как геофизичните данни да бъдат поставени в Wolfram Mathematica.
Сейсморазведка
За начало е необходимо да се направи малък преглед на геофизиката. Геофизиката е наука, която изучава физическите свойства на скалните породи. И тъй като породите имат различни свойства: електрически, магнитни, еластични, съществуват съответните методи на геофизиката: електрическа разузнаване, магнитно разузнаване, сейсмично разузнаване… В контекста на тази статия ще се фокусираме само върху сейсмичното разузнаване. Сейсмичното разузнаване е основният метод за търсене на нефт и газ. Методът се основава на възбуждане на еластични колебания и последваща регистрация на отговорите от скалите, които образуват изследваната територия. Възбуждането на колебанията се извършва на сушата (с динамит или неексплозивни вибрационни източници на еластични колебания) или в морето (с пневмопушки). Еластичните колебания се предават през дебелината на скалите, пречупвайки се и отразявайки се на границите на слоевете с различни свойства. Отразените вълни се връщат на повърхността и се регистрират от геофони на сушата (обикновено това са електродинамични уреди, основани на движението на магнит, висящ в бобина) или хидрофони в морето (основани на пьезоефекта). По времето на пристигане на вълните може да се съди за дълбочините на геологичните слоеве.
Сейсмичният кораб буксува оборудване
Пневмопушката възбужда еластични колебания
Вълните преминават през дебелината на скалните породи и се регистрират от хидрофони
Научно-изследователският кораб за геофизично разузнаване „Иван Губкин“ на кея при Благовещенския мост в Санкт Петербург
Модел на сейсмичния сигнал
Скалните породи имат различни физически свойства. За сейсмичното разузнаване преди всичко са важни еластичните свойства — скоростта на разпространение на еластичните колебания и плътността. Ако два слоя имат идентични или близки свойства, вълната „не ще усети“ границата между тях. Но ако скоростите на вълните в слоевете се различават, на границата на слоевете ще настъпи отражение. Колкото по-голяма е разликата в свойствата, толкова по-интензивно е отражението. Неговата интензивност се определя от коефициента на отражение (rc):
където ρ — плътност на породите, ν — скорост на вълните, 1 и 2 обозначават горния и долния слой.
Една от най-простите и често използвани модели на сеизмичен сигнал е свързващата модел, при която регистрираната сеизмична траса се представя като резултат от свързването на последователност от коефициенти на отражение с проучвателен импулс:
където s(t) — сеизмична траса, т.е. всичко, което е записал хидрофон или геофон в рамките на фиксирано време на регистрация, w(t) — сигнал, който генерира пневматичен пистолет, n(t) — случайен шум.
Нека да изчислим синтетична сеизмична траса за пример. Като изходен сигнал ще използваме широко използвания в сеизмичната разузнаване импулс на Рикер.
length=0.050; (*Дължина на сигнала*)
dt=0.001;(*Честота на проби на сигнала*)
t=Range[-length/2,(length)/2,dt];(*Време на сигнала*)
f=35;(*Централна честота*)
wavelength=(1.0-2.0*(Pi^2)*(f^2)*(t^2))*Exp[-(Pi^2)*(f^2)*(t^2)];
ListLinePlot[wavelength, Frame->True, PlotRange->Full, Filling->Axis, PlotStyle->Black,
PlotLabel->Style["Изначален вейвлет", Black, 20],
LabelStyle->Directive[Black, Italic],
FillingStyle->{White, Black}, ImageSize->Large, InterpolationOrder->2]Изходен сеизмичен импулс
Ще зададем две граници на дълбочини от 300 ms и 600 ms, а коефициентите на отражение ще бъдат случайни числа
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["Коефициенти на отражение", Black, 20],
LabelStyle->Directive[Black, Italic]]Последователност на коефициентите на отражение
Ще изчислим и визуализираме сеизмичната траса. Понеже коефициентите на отражение имат различни знаци, получаваме два знакоподобни отражения на сеизмичната траса.
traceExamle=ListConvolve[wavelength[[1;;;;1]],rcExample];
ListPlot[traceExamle,
PlotStyle->Black, Filling->0, Frame->True, Axes->False,
PlotLabel->Style["Сеизмична траса", 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 was generated at the Institute …nimum velocity of 1500 m/s and a maximum of 5500 m/s)»
За да визуализираме самата геоложка модел, можем да се обърнем към сеизмичните трасета чрез ключа «traces» (една от функциите на пакета е независимостта на ключовете от регистъра):
ArrayPlot[Transpose[marmousi["traces"]], PlotTheme -> "Detailed"]Моделът Marmousi
Понастоящем пакетът също така позволява да се зареждат данни на части от големи файлове, което прави възможна обработката на файлове с размери, достигащи десетки гигабайти. Освен това, функциите на пакета включват функции за експортиране на данни в .segy и частично дописване в края на файла.
Отделно заслужава да се отбележи функционалността на пакета при работа с комплексната структура на .segy файловете. Тъй като той позволява не само да се получават данни по ключове и индекси за отделни трасета, заглавия, но и да се променят с последващо записване в файла. Много технически детайли за реализацията на GeologyIO надхвърлят обхвата на тази статия и вероятно заслужават отделно описание.
Актуалността на спектралния анализ в сеизмичната разузнавателна дейност
Възможността за импорт на сеизмични материали в Wolfram Mathematica позволява използването на вградения функционал за обработка на сигнали за експериментални данни. Тъй като всяка сеизмична линия представлява времеви ред, то един от основните инструменти за тяхното изучаване е спектралният анализ. Сред предпоставките за анализ на честотния състав на сеизмичните данни можем да споменем например следните:
- Различните типове вълни се характеризират с различен честотен състав. Това позволява да се отделят полезните вълни и да се заглушат вълни-намеси.
- Такива свойства на скалите, като порьозност и насищане, могат да влияят на честотния състав. Това позволява да се отделят скали с по-добри свойства.
- Слоевете с различна дебелина предизвикват аномалии в различни честотни диапазони.
Третата точка е основна в контекста на тази статия. По-долу е приведен фрагмент от код за изчисление на сеизмични линии в случай на слой с променяща се дебелина — модел на клина. Този модел традиционно се изследва в сеизмологичните проучвания за анализ на интерференционните ефекти, когато вълни, отразени от много слоеве, се наклаждат една върху друга.
nx=200;(* Брой точки в решетката в X направление*)
ny=200;(* Брой точки в решетката в Y направление*)
T=2;(* Общо време на разпространение*)
(* Скорост и плътност*)
modellv=Table[4000,{i,1,ny},{j,1,nx}];(* Скорост на P-вълната в м/с*)
rho=Table[2200,{i,1,ny},{j,1,nx}];(* Плътност в г/см^3, използва се постоянна плътност*)
Table[modellv[[150-Round[i*0.5];;,i]]=4500;,{i,1,200}];
Table[modellv[[;;70,i]]=4500;,{i,1,200}];
(* Изчертаване на модела *)
MatrixPlot[modellv,PlotLabel->Style["Модел на слой",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 Hz, 55 Hz и 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]Набор от исходни сигнали с честоти 35 Hz, 55 Hz, 85 Hz
След извършване на изчисления на сеизмичните трасета и построяване на графики на амплитудите на отразената вълна, можем да видим, че за различни честоти аномалията се наблюдава при различни дебелини на пласта.
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["Амплитуди на отражения",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["Сейсмично трасе", 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; (*координата X на началото и края на снимането и стъпката на трасето*)
startXL=1165;EndXL=1615;stepXL=2; (*координата Y на началото и края на снимането и стъпката на трасето*)
numIL=(EndIL-startIL)/stepIL+1; (*броят трасета по оста X*)
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["Разлагане на вейвлети",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->{"Време","Честота(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}]; След прилагане на функцията появяват се нови масиви от данни, съответстващи на избраните честоти. В посочения по-горе пример това са честотите: 38Hz, 33Hz, 27Hz. Изборът на честоти обикновено се прави на база тестове - получават се резултатни карти за различни честотни комбинации и се избира най-информативната с оглед на геолога.
Ако е необходимо да се споделят резултатите с колеги или да се предоставят на клиента, то може да се използва функцията SEGYExport[] от пакета GeologyIO.
outputdata=seismic3DSEGY;
outputdata["traces",1;;-1]=tracesCWD[[All,3]];
outputdata["textheader"]="Резултат от вълновата декомпозиция";
outputdata["binaryheader","NumberDataTraces"]=Length[tracesCWD[[All,3]]];
SEGYExport["D:result.segy",outputdata];Разполагаща с три такива куба (нисък, среден и високочестотен компонент), обикновено се използва RGB смесване за съвместна визуализация на данните. На всеки от компонентите се присвоява свой цвят — червен, зелен, син. В Wolfram Mathematica това може да се направи, използвайки функцията .
В резултат на това получаваме изображения, по които може да се извършва геоложка интерпретация. Меандрите, които се визуализират на среза, позволяват да се очертаят палеорислата, които с по-голяма вероятност могат да бъдат резервоари и да съдържат находища на нефт. Търсенето и анализът на съвременни аналози на такава речна система позволяват да се определи най-перспективните части на меандра. Собствените русла се характеризират с мощни слоеве от добре сортиран пясъчник и са добри резервоари за нефт. Зоните извън „въжените“ аномалии са аналогични на съвременните поймени отложения. Пойменните отложения основно са представени от глинести породи и сондажите в тези зони ще бъдат неефективни.
RGB срез на куба от данни. В центъра (малко наляво от центъра) може да се проследи меандрящата река.
RGB срез на куба от данни. В лявата част може да се проследи меандрящата река.
В някои случаи качеството на сеизмичните данни позволява получаването на значително по-ясни изображения. Това зависи от методологията на полевите работи, оборудването, прилаганите алгоритми за потискане на шума. В такива случаи се виждат не само фрагменти от речните системи, но и цели дълги палеоречки.
RGB смесване на трите компонента на куба от сеизмични данни (хоризонтален срез). Дълбочината е приблизително 2km.
Изображение от сателит на река Волга в района на Саратов.
Заключение
В Wolfram Mathematica можете да анализирате сеизмични данни и да решавате практическите задачи, свързани с търсенето на полезни изкопаеми, а пакетът GeologyIO позволява да направите този процес по-удобен. Структурата на сеизмичните данни е такава, че използването на вградени методи за ускоряване на изчисленията (, ,…) е много ефективно и позволява обработка на големи обеми данни. В немалка степен това се дължи на особеностите на съхранение на данни в пакета GeologyIO. Струва си да се спомене, че пакетът може да бъде използван не само в областта на приложната сеизморазведка. Практически същите типове данни се използват и в георадиолокацията и сеизмологията. Ако имате предложения как да подобрите резултата, кои алгоритми за анализ на сигнал от арсенала на Wolfram Mathematica са приложими за тези данни или имате критични забележки – оставяйте коментари.
Източник: habr.com
