Дякуємо автору блогу за його доповідь
Запровадження
Ця нотатка написана слідами конференції і містить конспект доповіді, з якою я виступав. Захід відбувся у червні у місті Санкт-Петербурзі. Зважаючи на те, що працюю я в кварталі від місця проведення конференції, я не міг не відвідати цієї події. У 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 був створений в Institute …nimum velocity of 1500 m/s і 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