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

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

    Введение


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

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

    image

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

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

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

    Пакет GeologyIO


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

    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]

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



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

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



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

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

    После применения функции ContinuousWaveletTransform[] появляются новые массивы данных соответствующие выбранным частотам. В приведенном выше примере это частоты: 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 это можно сделать с использованием функции ColorCombine[].

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

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



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



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

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



    Изображение со спутника реки Волги в районе Саратова



    Заключение


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

    Комментарии 35

      0

      Как говорится, к вопросу, так ли все просто в нефтяные.

        0
        Вопроса три.
        1. Используются ли фазовые методы для улучшения различимости?
        2. Насколько хорошо видны рудные образования и их состав?
        3. Кроме сейсморазведки, используется нейтринная и другие радиационные виды, в том числе с использованием космического излучения?
          +1

          По п. 1 не очень понял этого вопроса — я не знаком с термином "фазовые методы". Если вы перефразируете вопрос, то попробую ответить.
          По п.2 что касается рудных объектов. Поиск и разведка рудных объектов существенно отличаются от таковых для нефти и газа. Для рудных полезных ископаемых на первое место выходят такие методы как электроразведка и магниторазведка. Как правило рудные тела не так контрастны по значениям упругих свойств.
          По. п3. В общем можно утверждать, что сейчас разведка месторождений нефти и газа — это комплексирование сейсморазведки и бурения. Под бурением я имею в виду и запись вдоль ствола скважины физических полей (электрокаротаж, измерения естественная радиоактивности...). В контексте поиска полезных ископаемых я не слышал о тех методах, что вы упомянули.

            +1
            1. В фазовом методе учитывается и фаза волны, пришедшей на приемник, при последующей математической обработке данных.
            2. Как раз учет фаз позволит четко разглядеть и рудные месторождения, в том числе.
            3. Примером «радиационного метода» является просвечивание египетской пирамиды космическим излучением. Приёмники нейтрино (мюонов или других видов излучения) ставятся вокруг объекта и через несколько дней (месяцев) накопления данных, постепенно проявляется внутреннее строение объекта. То же, что и рентген. Только в рентгеновском исследовании используется «свой» источник излучения, а в описанном случае используется космическое излучение.
              0

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

                0
                3. Примером «радиационного метода» является просвечивание египетской пирамиды космическим излучением. Приёмники нейтрино (мюонов или других видов излучения) ставятся вокруг объекта и через несколько дней (месяцев) накопления данных, постепенно проявляется внутреннее строение объекта. То же, что и рентген. Только в рентгеновском исследовании используется «свой» источник излучения, а в описанном случае используется космическое излучение.
                Слабо себе представляю как можно установить приёмники нейтрино вокруг подземного месторождения )) Научная фантастика.
                  0
                  Ну, зачем же «вокруг месторождения». Можно «невдалеке», за пару-тройку тыщ километров. Мы же не доставляем телескоп на Луну, чтобы рассмотреть её получше. И того что видно с Земли уже достаточно.
                  Насчет нейтрино, да, не знал что наука сильно отстала. Но мюонное просвечивание вроде бы освоено — "Просвечивание пирамиды продолжается".
                    –2
                    (флуд, с тараканами в чужих головах не разговариваю)
                      –1
                      Вы с фазовыми антенными решетками знакомы?
                      В данном случае, я сейчас должен для вас бесплатно изобрести методику (не знаю, существует ли она и не велосипед ли я буду изобретать). Вот прям сейчас мне что-то не хочется. То ли я устал, то ли разочаровался в жизни…
                        –2
                        Корона не жмёт?
                          –1
                          Да, слегка маловата. Надо больше, больше! :)))
                          Но, вообще, вы меня знаете? Или так, «с первого взгляда»?
                      0
                      Размер и цена нейтринного датчика обычно такие, что ну его нахрен такое исследование, см. супер-Камиоканде, например.
                        0
                        Ну, во первых, задачи разные. А во вторых, это не тот инструмент, который можно держать дома, для того чтобы померить, «чего это седни голова болит, не многовато ли нейтрино на улице».
                        Ситуация как и с первыми ЭВМ — Их мало, но можно арендовать время. Только вовремя заготовить и принести с собой перфокарты.
                    0
                    3. Вокруг объекта на получится установить приемники. Так как изучается толща земли. А если просто поставить на поверхности, то получим данные строения всего земного шара.
                      0
                      Вот и хорошо! Не надо никуда ездить. С одного места и так всё увидим. А если таких мест будет два, то получим стереокартинку. Вообще замечательно!
                      (Эх, какой я гениальный ...!).
                0
                Вертикальная ось — глубина

                Глубина, или время прихода отраженной волны?

                Дело хорошее, но остается открытым вопрос зачем, если есть пакеты типа OpendTect, который емнип даже в бесплатной версии умеет загрузить кубик, показать его, посчитать спекдекомп и сделать RGB Blend.
                  +2

                  Это конечно вы правильно заметили! Измеряем времена прихода волн. Поэтому сделал в статье оговорочку на этот счёт.


                  Про OpenedTect — OpenedTect, на мой взгляд, очень хорошее ПО и именно спектральная декомпозиция в нём сделана удобно. ПО действительно бесплатное, однако его установка и использование на рабочей станции возбраняются (по крайней и мере в некоторых случаях) политикой безопасности. Не готов ответить почему так, но такие сложности с OpenedTect (OT)имеются! Также, насколько я знаю, OT не распараллеливает расчёт этой задачи. Но ПО хорошее, особенно возможность дополнять функционал плагинами.


                  Wolfram Mathematica мною используется, как более универсальный инструмент. Помимо декомпозии использую для статистического анализа сейсмических данных и параметров пласта, для сейсмичекого моделирования.

                    +1
                    Согласен, самостоятельно написанное решение всегда гибче, тем более если стоит задача посчитать много комбинаций и выгрузить много карт.

                    Я сам с WM дела не имел никогда, в ней можно посчитать классификацию сейсмического сигнала?
                      0

                      Да, классификацию пробовал делать. Точнее один разок попробовал функцию FindClusters[] на синтетических данных.
                      https://my.pcloud.com/publink/show?code=XZYFHv7ZV6Ac1zi2IYHMfyr4oG31CpbU3GBk
                      На оригинальных данных ничего не делал пока что.

                        0
                        Интересно, надо будет попробовать. Я для этой задачи смотрю в сторону python и tf, но судя по всему есть смысл попробовать в wm
                  0
                  куб сейсмических данных

                  А почему вы называете параллелепипед кубом? Вопрос без подвоха, мы и сами в институте так делаем. Вот и интересно, это совпадение или геофизический сленг?
                    0

                    Я, признаться, даже и не задумывался об этом никогда! Куб и куб:). Думаю, к сленгу можно отнести!
                    Конечно все эти наборы данных имеют разный размер по осям. Если по X и Y это десятки километров, то по Z первые километры.

                    0
                    Наверное стоило бы сказать есть ли аналогичное в других математических пакетах, например в Matlab? Без этого непонятно: есть ли достижения или это только воспроизведение чего-то похожего.
                    Я, например, с ходу нашел OpenSeismoMatlab www.mathworks.com/matlabcentral/fileexchange/67069-openseismomatlab
                      +1

                      Согласен с тем, что аналоги есть. Matlab один из близких и популярен среди геофизиков. Поскольку я готовил доклад с целью демонстрации возможностей Wolfram на конференцию Wolfram Technologies, поэтому и счёл возможным не писать про аналогичные инструменты. Думаю, при других условиях такой обзор был бы необходим.
                      Есть достижение или нет. Как сказать!:) Считаю, что пакет GeologyIO, действительно полезная разработка. Хотя в публикациях я и ранее видел сейсмические данные в Wolfram.
                      Применение декомпозии к сейсмотрассам технология не новая (но эффективная), наиболее часто цитируют статью
                      Partyka, G., J. Gridley, and J. A. Lopez, 1999, Interpretational applications of spectral decomposition in reservoir characterization: The Leading Edge
                      Мне хотелось показать, что Wolfram может использоваться, как рабочий и удобный инструмент геофизика.

                      +1
                      здорово, что есть примеры геологических кейсов, особенно на хабре.
                      нефтянка — такое непаханное поле для таких примеров
                        –1
                        Фантазии про то, как двигались континенты и как из деревьев и динозавров образовывалась нефть IMHO лучше бы отделить от технического описания математических пакетов и поместить в отдельную статью ну хотя бы из уважения к памяти Николая Кудрявцева, который, кстати, окончил Ленинградский горный институт в 1922 году.

                        Кудрявцев заложил основы абиогенной теории происхождения нефти, ставшей популярной в советской, а затем в российской и украинской науке (из Википедии ).
                          +1

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


                          Что же касается движения литосферных плит, то как мне кажется, слово "фантазия" не очень уместно. Всё таки модель движущихся плит объясняет многие наблюдаемые явления, в частности зоны спрединга и субдукции, возраст и магнитные аномалии океанической коры. Были ранее и есть сейчас альтернативные модели — океанизация Земли, например.


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

                            0
                            А по факту, вам часто встречались резервуары реки? Или все же они больше похожи на озера или отдельные заливы древнего моря. Именно продуктивные резервуары.
                              +1

                              Все перечисленные типы имеют место! Получается так: в среднеюрское время была суша, то резервуары связаны с руслами палеореками; в поздней юре условия сменились на мелкое море и стали накапливаться песчаники вдольбереговых баров (как и сейчас на морских побережьях) и месторождения приурочены к этим этим барам. Озерные отложения также фиксируются в скважинах, но они как правило более глинистого состава, имеют меньшие толщины и поэтому менее интересны. Это про Западную Сибирь.
                              В других регионах (Тимано-Печорская провинция, Волго-Урал, Прикаспий) разрез сложен карбонатными породами и мы "видим" соответствующие объекты — барьерные и одиночные рифы.

                                0
                                Причем рифы зачастую видно просто на амплитудах, для них даже спекдекомп считать не обязательно
                                  +1

                                  Да, рифы часто хорошо видны "невооруженным глазом" — и в обычных амплитудах и рельефе поверхности. Мне кажется декомпозиция наиболее информативна для палеоречек, каньонов.

                                    0
                                    Пожалуй стоит пояснить для не-экспертов, в рифах существенно выше скорость волны, чем в других карбонатах, из-за этого в кровле рифа получается сильный скачок импеданса (акустический импеданс — произведение скорости на плотность, первая формула в статье), что приводит к высокой амплитуде отражения.

                                    Кстати из-за повышенной скорости в рифах зачастую отражения ниже рифа «подтягиваются» вверх, и весьма заметно. Это один из хороших диагностических признаков рифов на разрезах амплитуд, но важно учесть это при преобразовании из временного домена в глубинный, а то можно получить артефакты и «фиктивные» структуры, которых по факту нет.
                                      +1

                                      Согласен по всем пунктам.
                                      Скорости продольных волн в песчаниках ~3000м/с, а известняках ~5000м/с.

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


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

                          А глупости про нефть из мантии я даже комментировать не буду.
                            –1
                            Я больше скажу, эти фантазии основаны на палеомагнитных данных, позволяющих восстановить координаты данной точки в момент формирования горной породы.

                            Можно отследить на какое-то время назад, но на достаточно большой промежуток времени это уже скорее фантазия. Картинка из энциклопедии Британика на 152 млн лет назад — фантазия, базирующаяся на фантазиях Вегенера. Сам Вегенер заявлял, что эта фантазия — всего лишь возможный вариант, нужны дополнительные исследования. Бедняга тогда, сидя на какой-то станции во льдах не имел глобуса, поэтому фантазировал исходя из плоской модели Земли. Может проблемы пространственного восприятия у него, я не знаю. Сами можете посмотреть его манипуляции с континентами — ни один пазл не уходит за пределы картинки, чтобы появиться с другого края.

                            А глупости про нефть из мантии я даже комментировать не буду.

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

                            Так что описание программных продуктов — это одна песня, а споры про истинность разных картинок о том, как оно было сотни миллионов лет назад — это другая песня. Зачем их смешивать в исполнении одного номера — непонятно.

                          Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

                          Самое читаемое