Хабрахабр

Wolfram Mathematica в Геофизике

Благодарим автора блога Антона Екименко за его доклад

Введение

Мероприятие состоялось в июне в городе Санкт-Петербурге. Эта заметка написана по следам конференции Wolfram Russian Technology Conference и содержит конспект доклада, с которым я выступал. В 2016 и 2017 годах я слушал доклады конференции, а в этом году выступил с докладом. Учитывая то, что работаю я в квартале от места проведения конференции, я не мог не посетить это событие. Во–первых, появилась интересная (как мне кажется) тема, которую мы развиваем с Кириллом Беловым, а во-вторых, после длительного изучения законодательства РФ в части санкционной политики, на предприятии, где я тружусь появилось аж целых две лицензии Wolfram Mathematica.

На визитной страничке конференции используется изображение Казанского собора. Прежде, чем перейти к теме моего выступления я бы хотел отметить хорошую организацию мероприятия. Собор является одной из главных достопримечательностей Петербурга и его очень хорошо видно из зала, в котором проходила конференция.

image

Во время регистрации выдавались небольшие сувенирчики (игрушка – мигающий спайки, ручка, наклейки с символикой Wolfram). На входе в СПбГЭУ участников встречали помощники из числа студентов – не позволяли заблудиться. Про вкусный кофе и пирожочки я уже отмечал на стене группы – повара молодцы. Обед и кофе-брейк также были включены в расписание конференции. Этой вводной частью я бы хотел подчеркнуть, что само мероприятие, его формат и место проведения уже приносят положительные эмоции.

Спектральный анализ сейсмических данных или «куда бежали древние реки». Доклад, который был подготовлен мною и Кириллом Беловым называется «Использование Wolfram Mathematica для решения задач прикладной геофизики. Содержание доклада покрывает две части: во-первых, это использование алгоритмов, имеющихся в Wolfram Mathematica для анализа геофизических данных, а во-вторых, это то как геофизические данные поместить в Wolfram Mathematica.

Сейсморазведка

Геофизика — это наука, изучающая физические свойства горных пород. Для начала нужно сделать небольшой экскурс в геофизику. Сейсморазведка является основным методом поиска нефти и газа. Ну а поскольку породы имеют разные свойства: электрические, магнитные, упругие, то существуют соответствующие методы геофизики: электроразведка, магниторазведка, сейсморазведка… В контексте данной статьи подробнее затронем только сейсмическую разведку. Возбуждение колебаний осуществляется на суше (динамитом или невзрывными вибрационными источниками упругих колебаний) или на море (пневмопушками). Метод основан на возбуждении упругих колебаний и последующей регистрации отклика от горных пород слагающих исследуемую территорию. Отражённые волны возвращаются на поверхность и регистрируются геофонами на суше (как правило это электродинамические приборы, основанные на движении магнита, подвешенного в катушке) или гидрофонами в море (основаны на пьезоэффекте). Упругие колебания распространяются через толщу горных пород преломляясь и отражаясь на границах слоёв с разными свойствами. По времени прихода волн можно судить о глубинах геологических слоёв.

Сейсмическое судно буксирует оборудование

Пневмопушка возбуждает упругие колебания

Волны проходят через толщу горных пород и регистрируются гидрофонами

Научно-исследовательское судно геофизической разведки «Иван Губкин» на причале у Благовещенского моста в Петербурге

Модель сейсмического сигнала

Для сейсморазведки прежде всего важны упругие свойства — скорость распространения упругих колебаний и плотность. Горные породы имеют разные физические свойства. Если же скорости волн в слоях будут отличаться, то на границе слоёв возникнет отражение. Если два слоя будут иметь одинаковые или близкие свойства то волна «не заметит» границу между ними. Его интенсивность будет определяться коэффициентом отражения (rc): Чем больше разница в свойствах тем интенсивнее отражение.

где ρ — плотность пород, ν — скорость волн, 1 и 2 обозначают верхний и нижний слои.

Одной из наиболее простых и часто используемых моделей сейсмического сигнала является свёрточная модель, когда зарегистрированная сейсмическая трасса представляется, как результат свёртки последовательности коэффициентов отражения с зондирующим импульсом:

все что записал гидрофон или геофон в течение фиксированного времени регистрации, w(t) — сигнал, который генерирует пневмопушка, n(t) — случайный шум. где s(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->,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 Гб, а объём исходного «сырого» материала может составить сотни гигабайт.

Как работать с ними в Wolfram Mathematica?

Благодаря ответам сообщества решение было найдено очень быстро. Началом разработки пакета стал вопрос на стене VK группы русскоязычной поддержки. Соответствующий пост на стене Wolfram Comunity даже был отмечен модераторами. И в результате переросло в серьёзную разработку. На текущий момент пакет поддерживает работу со следующими типами данных, которые активно используются в геологической отрасли:

  1. импорт картографических данных формата ZMAP и IRAP
  2. импорт измерений в скважинах формата LAS
  3. ввод и вывод сейсмических файлов формата SEGY

выполнить следующий код в любом блокноте Mathematica: Чтобы установить пакет необходимо следовать инструкциям на странице загрузки собранного пакета, т.е.

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 и снабдить пакет тестовыми файлами для первого знакомства. Пакет разрабатывается с использованием Wolfram Workbench.

Используя эту модель разработчики тестирую собственные алгоритмы моделирования волнового поля, обработки данных, инверсии сейсмических трасс и т.д. Таким файлом, в частности является файл «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 позволяет использовать встроенный функционал обработки сигналов для экспериментальных данных. Среди предпосылок анализу частотного состава сейсмических данных можно назвать, например, следующие:

  1. Разные типы волн характеризуются разным частотным составом. Это позволяет выделить полезные волны и подавить волны-помехи.
  2. Такие свойства горных пород, как пористость и насыщенность могут влиять на частотный состав. Это позволяет выделять горные породы с лучшими свойствами.
  3. Слои с разной толщиной обуславливают аномалии в разных частотных диапазонах.

Ниже приведён фрагмент кода для расчета сейсмических трасс в случае слоя с меняющейся толщиной — модель клина. Третий пункт является основным в контексте данной статьи. Эта модель традиционно изучается в сейсморазведке для анлиза интерференционных эффектов, когда волны отраженные от многих слоёв накладываются друг на друга.

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]

Сейсмические трассы для модели клина

Как можно заметить, его интерпретация может выполняться и на интуитивном уровне, поскольку геометрия отражённых волн однозначно соответствует модели, которая была задана ранее. Последовательность сейсмических трасс изображенная на этом рисунке называется сейсмическим разрезом. Начиная с 31-й трассы отражения начинают интерферировать. Если более детально анализировать трассы, то можно заметить, что трассы с 1-й по, примерно, 30-ю не отличаются — отражение от кровли пласта и от подошвы не накладываются друг на друга. И, хотя, в модели коэффициенты отражения не меняются по горизонтали — сейсмические трассы меняют свою интенсивность при изменении толщины пласта.

Начиная с 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. В такой ситуации очень удобным становится использование систем компьютерной математики и языков программирования высокого уровня, которые позволяют использовать широкую алгоритмическую базу и, вместе с тем, берут на себя много рутины. Нецелесообразно писать богатый функционал интерактивной работы с данными -важнее обеспечить загрузку из общепринятого формата, применить к ним желаемые алгоритмы и выгрузить обратно во внеший формат.

Следуя предложенной схеме, загрузим оригинальные сейсмические данные и отобразим их в 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] // TableForm

C 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 это можно сделать с помощью Opacity[] и 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]

Изображение куба сейсмических данных с использованием функций 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]

Результат функции WaveletScalogram[]
Скалограмма.

А примение этой функции ко всему набору трасс осуществется с испольшованием функции Table[]. В Wolfram Language для вейвлет преобразования используется функция ContinuousWaveletTransform[]. В приведённом примере в распаралеливании нет необходимости — объём данных не велик, но при работе с экпериментальными наборами данных содержащими сотни тысяч трасс это является необходимостью. Здесь стоит отметить одну из сильных сторон Wolfram Mathematica возможноть использовать распрараллеливание ParallelTable[].

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

В приведенном выше примере это частоты: 38Гц, 33Гц, 27Гц. После применения функции ContinuousWaveletTransform[] появляются новые массивы данных соответствующие выбранным частотам. Выбор частот осущетвляется чаще всего на основе тестирования -получают результативные карты для разных частотных комбинаций и выбирают наиболее информативный с точки зрения геолога.

Если результатами необходимо поделиться с коллегами или предоставить их заказчику, то можно использовать функцию 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];

Каждой из компонент присваивается свой цвет — red, green, blue. Имея в распоряжении три таких куба (низкочастотную, среднечастотную и высокочастотную компоненты), как правило, используют RGB смешивание для совместной визуализации данных. В Wolfram Mathematica это можно сделать с использованием функции ColorCombine[].

Меандры, которые фисируются на срезе позволяют оконтурить палеорусла, которые с большей вероятностью могут быть резервуарами и содержать запасы нефти. В результате получаются изображения, по которым можно выполнять гелогическую интерпретацию. Собственно русла характеризуются мощными слоями хорошо сортированного песчаника и являются хорошим резервуаром для нефти. Поиск и анализ современных аналогов такой речной системы позволяет определить наиболее перспективные части меандр. Пойменные отложения в основном представлены глинистыми породами и бурение в эти зоны будет неэффекивным. Участки за пределами «шнурковых» аномалий аналогичны современным пойменным отложениям.

В центре (несколько левее центра) можно проследить меандрирующую реку.

RGB срез куба данных. RGB срез куба данных. В левой части можно проследить меандрирующую реку.

Это зависит от методики полевых работ, оборудования, примяемых алгоритомв шумоподавления. В некоторых случаях качество сейсмических данных позволяет получать существенно более чёткие изображения. В таких случаях видны не только фрагменты речных системы, но и целые протяженные палеореки.

Глубина примерно 2км.

Изображение со спутника реки Волги в районе Саратова
RGB смешивание трёх компонент куба сейсмических данных (горизонтальный срез).

Заключение

Структура сейсмических данных такова, что использование встроенных способов ускорения расчетов (ParallelTable[], ParallelDo[],...) является очень эффективным и позволяет обрабатывать большие объёмы данных. В Wolfram Mathematica можно анализировать сейсмические данные и решать прикладные задачи связанные с поиском полезных ископаемых, а пакет GeologyIO позволяет сделать этот процесс более удобным. К слову сказать, пакет может использоваться не только в области прикладной сейсморазведки. В немалой степени этому способствуют особенности хранения данных пакета GeologyIO. Практически такие же типы данных используются в георадиолокации и сейсмологии.Если у вас есть предложения как улушить результат, какие алгоритмы анализа сигнала из арсенала Wolfram Mathematica применимы к таким данным или у вас имеются критические замечания -оставляйте комментарии.

Теги
Показать больше

Похожие статьи

Добавить комментарий

Ваш e-mail не будет опубликован. Обязательные поля помечены *

Кнопка «Наверх»
Закрыть