УДК 550.8.053
В. И. Голубев
Московский физико-технический институт (государственный университет)
Методика отображения и интерпретации результатов полноволновых сейсмических расчётов
Исследуется проблема наглядного представления и интерпретации результатов компьютерного моделирования динамических процессов, инициируемых в гетерогенных геологических средах в процессе сейсмической разведки на нефтегазовых месторождениях. Автором предложена методика построения синтетических сейсмограмм по заданным профилям на основе расчётных данных как в 2D-, так и в ЗБ-случаях. Она реализована в виде компьютерной программы «Seismograph», производящей полуавтоматическое построение наглядных графиков и их сохранение в виде графических файлов для дальнейшего анализа. В работе приведены примеры практически значимых геофизических результатов, полученных с использованием «Seismograph» на этапе интерпретации. Описан формат хранения данных численных расчётов, эффективный как на этапе сохранения, так и на этапе постпроцессинга. Отдельно рассмотрены эффективные алгоритмы сохранения данных на диск в случае однопроцессорного и многопроцессорного вариантов расчёта.
Ключевые слова: математическое моделирование, сеточно-характеристический численный метод, сейсморазведка, интерпретация экспериментальных данных.
1. Введение
Согласно сегодняшним прогнозам потребление энергии в мире будет продолжать стремительно расти в связи с развитием науки и техники. Несмотря на развитие альтернативных источников энергии (ветряные, геотермальные, солнечные, гидроэнергетические и биотопливные), а также рост значимости ядерной энергетики, традиционные источники энергии, такие как нефть и природный газ, остаются важнейшим ресурсом. В то время как запасы традиционной (легко извлекаемой) нефти существенно истощены, во всём мире активно начинается разработка нетрадиционной (тяжелой) нефти. Основной сложностью данного процесса является отсутствие эффективных подходов по извлечению значительного процента от объёма всего месторождения. Кроме того, остаётся открытым вопрос о методах эффективного поиска трещиноватых кластеров, залегающих на большой глубине и содержащих нефте- и газопродукты.
Одним из подходов к решению данных проблем является разработка новых компьютерных программ и комплексов по расчёту процесса сейсмического исследования недр Земли, которые позволят с высокой степенью точности провести моделирование полевых экспериментов. При этом появляется возможность задания определённой геометрии среды (количество и ориентация флюидосодержащих трещин, слоёв различных пород и т.д.) с последующим получением синтетического отклика. Для регистрации сигнала-отклика на практике широко применяются специальные приборы — сейсмоприёмники. Они позволяют регистрировать смещение земной поверхности в точке их установки как функцию времени. Одним из видов их классификации является деление по типу регистрации: однокомпо-нентные и многокомпонентные датчики. При однокомпонентной регистрации записывается лишь вертикальная составляющая смещения. При многокомпонентной регистрации также фиксируются проекции смещения на север и восток. Кроме того, в зависимости от физического принципа, на котором основан датчик, возможна регистрация смещения, скорости или ускорения дневной поверхности.
Таким образом, с применением компьютерного моделирования могут быть решены следующие научно-технические задачи:
2.
• оценка возможности регистрации полезного отклика от трегциноватои структуры с заданными параметрами на фоне помех;
• разработка новых методов и методик обработки данных полевых измерений с целью выделения сигнала, несущего информацию о структуре месторождения;
• поиск оптимальной расстановки источников возмущений и регистрирующих сейсмо-приёмников для изучения геологической структуры с заданными свойствами.
Компьютерное моделирование сейсмических процессов в геологических средах
На кафедре информатики Московского физико-технического института (государственного университета) ведётся разработка суперкомпьютерных программных комплексов для моделирования динамических процессов в геологических средах. Для описания их поведения при динамических нагрузках используется явное выделение всех разномасштабных неоднородностей на этапе построения расчётной сетки и определяющие уравнения линейно-упругого тела, решамемые сеточно-характеристическим численным методом [1]. Определяющая система уравнений в случае линейно-упругого тела имеет вид [2]:
Р • Щ = агз,
= \(е 11 + е2 2 + £з з + 2ре\3,
(1) (2) (3)
Здесь р - плотность среды, ^ - компоненты скорости смещения, а^ - компоненты тензора напряжений, - компоненты тензора деформацпй, А и р - упругие постоянные среды. К сожалению, в общем случае система уравнений (1)-(3) не может быть решена аналитически. Для её решения используется, во-первых, дискретизация на расчётную сетку (неструктурную треугольную или тетраэдральную, структурную параллелепипедную или гексаэдральную). Во-вторых, расщепление по координатам, которое позволяет преобразовать исходную систему в набор одномерных систем уравнений. Поскольку данные системы гиперболические, они обладают полным набором собственных векторов с действительными собственными значениями. При переходе в систему координат, связанную с собственными векторами, каждая из них распадается на набор независимых скалярных уравнений переноса, каждое из которых решается с помощью метода характеристик. После этого общее решение восстанавливается в каждом узле расчётной сетки на новом временном слое. Схематически данный подход изображён на рис. 1.
Рис. 1. Переход на следующий временной слой методом характеристик
Для обеспечения высокой точности решения при сохранении приемлемой скорости расчётов широко распространены два основных подхода: расширение сеточного шаблона и использование продолженных систем уравнений. Рассмотрим второй из них на примере двумерного уравнения переноса:
ди ди ди ^
Ы х дх у ду
Введём следующие обозначения:
du du d2u .g.
dx' dy7 dxdy
Алгоритм решения состоит из двух аналогичных шагов расщепления по координатам. При этом, например, на шаге расщепления по X решается следующая система уравнений:
' du I _ du =0 dt + Cxdx =
I r 9v =0 dt + Cxdx = 0' dw | _ dw =o dt + Cx dx = 0'
§P + с §P =0 dt + Cxdx =
Комбинированием уравнений (1-2) и (3-4) записанной системы получаем две независимые продолженные системы уравнений, которые решаются с третьим порядком точности на двухточечном шаблоне. Итоговый порядок полученной схемы не понижается из-за расщепления и равен 3.
Необходимо отметить, что результатом численного расчёта являются значения всех компонент тензора напряжений и всех компонент скорости смещения в каждом узле сетки в последовательные моменты времени. Таким образом, появляется возможность проведения симуляции процесса сейсмической разведки в геологической среде с заданной структурой и сравнения результатов расчета с реальными данными наземных измерений. Целью работы являлась разработка эффективного подхода и инструмента для построения синтетических сейсмограмм, адаптированных под ведущиеся на кафедре разработки.
3. Сохранение результатов численного расчёта
Одной из особенностей компьютерного моделирования на расчётной сетке (особенно задач в ЗО-постановке) явлется большой объём обрабатываемой информации. Например, расчётная сетка из 1 миллиона узлов с хранящейся в них информацией при использовании типа float для упругой задачи в 2D занимает порядка 40 Mb и 80 Mb в 3D соответственно. При переходе к типу double эти цифры удваиваются. При среднем числе шагов по времени порядка 1000 получается огромная величина. Поскольку операция обращения к жёсткому диску является наиболее затратной, необходимость сохранения результатов может наложить существенные ограничения также и на время расчёта. Кроме проблем с быстродействием и хранением возникает также и проблема постобработки результатов. Отметим, что в полевом эксперименте единственной наблюдаемой величиной является сейсмограмма. В связи с этим представляется перспективным построение синтетических сейсмограмм по результатам численного расчёта, например, для проведения сравнения с экспериментальными данными. Для этого достаточно сохранения всех компонент скорости смещения в заранее фиксированных узлах сетки - местах расстановки сейсмодатчиков. Для сохранения данных с сейсмометров в разрабатываемых программных комплексах было решено использовать формат CSV (Comma-Separated Values) с сохранением числовых значений в ASCII-формате. Таким образом, стандартная структура файла (ЗО-случай) имеет следующий вид:
Time;Vx;Vy;Vz;
0;0;0;0
0.01;0;0;0
0.02;1.45е-2;1е-2;1.2е-5
Необходимо отметить, что задачи полноволнового моделирования сейсмических процессов отличаются повышенной вычислительной сложностью. Для обеспечения возможности расчёта на сетках достаточной детальности за приемлемое время используются параллельные алгоритмы, ориентированные на работу как на системах с общей, так и с распределённой памятью (кластерные системы). Их программная реализация основывается на
v
использовании технологий ОрепМР и MPI. Остановимся подробнее на процессе сохранения данных с сейсмоприёмников. Такой процесс при расчёте в один поток не представляет сложностей. На начальном этапе производится поиск номеров узлов, в которых устанавливаются приёмники. При этом производится полный перебор всех узлов расчётной сетки и находится узел, ближайший к месту установки приёмника. Далее на каждом временном шаге производится последовательное сохранение в файл на диске всех компонент вектора скорости в каждом из отобранных ранее узлов. В случае, когда работает более одного потока, необходима синхронизация между ними для устранения race-condition при попытке записи на диск. Кроме того, так как в сейсмограмме важна последовательность объединения отдельных сейсмотрасс, каждый из процессов должен знать, слева или справа от группы приёмников, обсчитываемых каждым из остальных процессов, находится его группа датчиков. Для преодоления данных сложностей был реализован следующий алгоритм. Изначально принимается решение о сохранении каждой из сейсмотрасс в отдельный файл. Далее, при задании расстановки сейсмодатчиков они указываются в той последовательности, которую они должны сохранять при объединении в косу. При этом порядковый номер приёмника записывается в название файла, отвечающего за хранение данных с него на диске. После окончания расчёта образуется набор файлов, который затем объединяется в финальную сейсмограмму согласно установленной нумерации. Описанный подход позволяет эффективно сохранять данные для построения синтетических сейсмограмм как при последовательном, так и параллельном выполнении расчётных программ (в 2D- и 3D-случаях).
4. Структура программы «Seismograph» и её возможности
Автором разработана программа «Seismograph», позволяющая в полуавтоматическом режиме проводить построение синтетических сейсмограмм по результатам компьютерного моделирования. Для обеспечения кросс-платформенности совместно с высокой скоростью работы и простотой поддержки кода в качестве языка разработки выбран Python [3]. Для реализации GUI (Graphic User Interface) выбрана библиотека Tkinter [3] в силу кросс-платформенности, наличия версий для Linux, Microsoft Windows, Mac OS, а также вхождения в стандартную библиотеку Python. Использование командной Open Source утилиты Gnuplot позволило с минимальными усилиями реализовать отрисовку сейсмограммы по подготовленным входным данным и сохранение её в формате PNG (Portable Network Graphics) для дальнейшего использования в презентациях, публикациях и т.д. Разработанная программа была протестирована и успешно работает как на системах семейства Windows, так и Unix.
Одной из задач, стоящих перед автором на момент начала работы над проектом, являлось создание простого и одновременно интуитивно понятного интерфейса. Одновременно с этим хотелось предоставить пользователю возможность настройки параметров отображения сейсмограммы под конкретные особенности расчётной задачи (разномасштабные сейсмические отклики, 20/30-область интегрирования, изменяемая плотность расстановки сейсмометров). При реализации пользовательского интерфейса был избран минималисти-ческий путь. Программа содержит всего одно окно и запуск построения сейсмограммы контролируется одной кнопкой «Draw». Настройка всех параметров вынесена в один столбец в левой части экрана, а остальная часть рабочего пространства отведена под отображение построенного графика (см. рис. 2).
Опишем подробнее влияние значения каждого параметра на вид результирующей сейсмограммы. Во-первых, доступен выбор «Draw Quantity» (отображаемой величины) = Vx, Vv, Vz (ЗО-случай), Vabs - любой компоненты вектора скорости смещения и его модуля. Для наглядного изображения знакопеременных значений (отдельных компонент вектора скорости) реализованы настройки «Paint Level», «Down Limit», «Up Limit», «Paint below/above». Последнее поле может принимать значения «below» и «above», определяющие, как закрашивать отрицательную или положительную часть каждой сейсмотрассы соответственно.
При этом нулевой уровень закраски может быть изменен заданием значения из диапазона [ 1, 11 в поле «Paint Level». Значение «попе» означает выбор значения но умолчанию, равного 0 для данной настройки. Возможность закраски нижней части графика бывает удобна при анализе отдельных компонент отклика при интенсивных поперечных волнах, тогда как для анализа модуля скорости, очевидно, может быть использовано лишь значение «above» (при нулевом «Paint Level»). Поддерживается также задание верхних) («Up Level») и нижних) («Down Level») порогов отображения еейемотраее, что может быть использовано для одновременного отображения на сейсмограмме волн с большой разницей в амплитуде. По умолчанию пределы не устанавливаются, и в случае, если показания приемника превышают расстояние между ближайшими прямыми, происходит наложение данных соседних еейемотраее, которое обычно затрудняет интерпретацию результатов.
Draw
"Distance" le-4 Amplification
\l
Paint Level none
Down Limit none
Up Limit [none Start Receiver
\q
Receivers Numbei
M I
Every К Receiver
Draw Quantity
Vx
From Time Numb
start Until Time Numt"
end
Paint below/abov above
Time Tics
fo.i
Рис. 2. Главное окно программы «Seismograph»
Часто бывает удобным проводить детальный анализ определенной части (ограниченной но номерам сейсмонриемников и диапазону времени записи) всех сохраненных данных. Возможность ограничения отображаемого временного диапазона реализована параметрами «From Time Number» и «Until Time Number». В них задаются номера отсчетов времени датчиков, между которыми будут отображены все моменты времени. Кроме целых значений ноле также может принимать строковые значения «start» и «end», означающие начало записи и конец записи (расчета) соответственно. Умелое использование данных параметров позволяет с высокой точностью выделить интересующие интерпретатора сигналы отклики. Предусмотрена также возможность пространственного прореживания информации. Кроме построения всей сейсмотраммы но расположенным через равные промежутки вдоль прямой линии сейсмодатчикам возможен выбор: первого датчика в серии («Start Receiver»), количества датчиков в серии («Receivers Number») и промежутка пропускаемых датчиков («Every К Receiver»).
Поля «Distance» и «Amplification» позволяют изменять расстояние между сейсмотрас-еами с ближайших датчиков и усиление каждой сейсмотрассы. Они могут быть использованы для получения более наглядной волновой картины путём надлежащих) масштабирования. Поле «Time Tics» задаёт шах' шкалы но оси ординат, что может быть использовано, например, для повышения точности ручного пикирования первых вступлений. Поле «Gap (m)» позволяет задать в имени сохраняемого файла информацию о расстоянии между ближайшими сейсмотрассами в полученной сейсмограмме, что может быть полезно при обработке множества данных в широком диапазоне параметров для систематизации результатов обработки.
Процедура отрисовки сейсмограммы по заданным параметрам состоит из следующих шагов: подготовка входных данных для Gnuplot (создание временных файлов), вызов Gnuplot, заканчивающийся сохранением изображения в формате PNG в каталог Results с соответствующим именем файла, и отображение построенной сейсмограммы в главном окне программы. На первом шаге в соответствии с заданными пользователем параметрами происходит выборка анализируемых значений (или их генерация для случая «Vabs») из входного файла «receivers.csv» и их сохранение в промежуточный файл с данными. Часть настроек для Gnuplot (например, размер сохраняемой картинки в пикселях), скрытая от обычного пользователя, хранится статически в файле с настройками «configuration.ini». Остальные настройки (например, имя файла, шаг шкалы времени и т.д.) генерируются интерактивно согласно введённым пользователем значениям в GUI и совместно с фиксированной частью сохраняются в промежуточный файл настроек. На втором шаге происходит создание отдельного процесса в системе, в котором выполняется Gnuplot с параметрами, переданными в промежуточных файлах. Таким образом, после завершения его работы на диске образуется графический файл с результатами, который на третьем шаге стандартными средствами Tkinter подгружается в главное окно программы. Описанный выше «раздельный» (модульный) подход, во-первых, позволяет в существенной мере использовать наработки Open Source сообщества, а также даёт возможность опытному пользователю, знакомому со скриптовыми языками, производить свободную модификацию «Seismograph».
5. Полученные с использованием «Seismograph» важные научно-прикладные результаты
Первая версия программы (вместо GUI использовался консольный интерфейс, поддерживалась только Microsoft Windows) была разработана летом 2010 года. Примерно через полгода была осуществлена доработка GUI, а ещё через полгода - портирование «Seismograph» на Linux. Программа с самого начала получила высокую востребованность в научной группе на кафедре информатики, в которой автор работает над кандидатской диссертацией. Она активно используется научным коллективом для визуализации результатов численного моделирования задач сейсморазведки и глобальной сейсмики в 20-постановке на неструктурных треугольных сетках, а также в ЗО-постановке на структурных криволинейных и неструктурных тетраэдральных сетках.
Так, в работе [4] интерпретацией синтетических сейсмограмм, построенных по данным численного 2D-моделирования, установлено, что
• кластер однородно ориентированных субвертикальных макротрещин формирует при падении на него плоской волны продольных колебаний многофазные фронты рассеянных продольных и обменных волн;
• в условиях субгоризонтальных границ, возбуждения плоского фронта продольных колебаний и флюидонасыщенного продуктивного трещиноватого пласта имеется возможность записи Х-компоненты регистрации фронта рассеянных поперечных колебаний как основного отклика от кластера субвертикальных макротрещин;
• структура волновых откликов от кластера характеризуется наличием интенсивных многофазных горизонтальных фронтов (в пределах ширины кластера) на времени прихода продольной дифрагированной волны при регистрации вертикальной компоненты и на времени прихода обменной дифрагированноый волны при регистрации горизонтальной компоненты. Вне зон этих фронтов регистрируются на порядок более слабые фрагменты гипербол дифрагированных волн от краевых макротрещин кластера;
• отмечается резкое различие в характере отклика в зависимости от насыщения системы (кластера) макротрещин жидкостью (водой или нефтью) или газом. В первом
случае они наилучшим образом проявляются при регистрации горизонтальной (X) компоненты записи в виде фронта рассеянных обменных волн. Во втором - при регистрации вертикальной (Z) компоненты в виде многофазного фронта рассеянных продольных волн.
В работе [5] проведено численное моделирование сейсмических откликов от кластера субвертикальных газонасыщенных и флюидонасыщенных макротрещин в 20-постановке. Анализ синтетических сейсмограмм, построенных в программе «Seismograph», позволил получить следующие практически важные результаты:
• протяжённость плоской части сейсмического отклика приблизительно равна горизонтальной протяжённости кластера как для газонасыщенных, так и для флюидонасыщенных макротрещин;
• зависимость периода сферической волны, распространяющейся от краев кластера, от его вертикальной протяжённости близка к линейной. Отличие между различными типами заполнителя выражается в изменении коэффициента наклона;
• зависимость периода плоской части сейсмического отклика от расстояния между трещинами близка к линейной. Отличие между различными типами заполнителя выражается в изменении коэффициента наклона.
В работе [6] исследовалось влияние неидеальности макротрещин (смещения отдельных звеньев), изменчивости наклона отдельных трещин, а также смены типа начального возмущения на «точечный источник» на сейсмический отклик от кластера субвертикальных газонасыщенных и флюидонасыщенных макротрещин в 20-постановке. Интерпретация синтетических сейсмограмм позволила выявить следующие эмпирические закономерности:
• при боковых смещениях отдельных звеньев от оси трещины порядка 20-40 см изменение энергии отклика незначительно и составляет порядка 5-10 %;
• неоднородность в ориентации трещин в пределах кластера и небольшие (единицы градусов) отклонения от единого направления слабо влияют на формирование фронта обменных рассеянных волн;
• преимущество регистрации Х-компоненты фронта обменных рассеянных волн как показателя зоны развития трещин сохраняется при «точечном» возбуждении в системе МОГТ при вводе стандартных кинематических поправок и раздельном суммировании правых и левых относительно пункта взрыва частей записей.
В работе [7] исследовался процесс инициации сейсмической активности в геологическом вмещающем массиве в ЗО-постановке. Источник возмущения выбирался в виде стандартной модели очага землетрясения («подвижка по разлому») с заданной плоскостью разрыва и направлением проскальзывания. Анализ синтетических сейсмограмм подтвердил зависимость времени прихода продольных и поперечных волн от скоростной модели среды.
В работе [8] проведено исследование структуры упругих возмущений, распространяющихся из очага землетрясения в 20-постановке. В результате анализа синтетических сейсмограмм выявлено дополнительное воздействие переотражённых и дифрагированных волн на дневную поверхность при наличии приповерхностного карстового отложения. Экспериментально показана возможность регистрации всех кратных волн, возникающих на границах между геологическими слоями с различными упругими характеристиками.
6. Заключение
В статье предложена методика наглядного представления результатов численных расчётов динамических процессов в неоднородных геологических средах на основе построения синтетических сейсмограмм в полуавтоматическом режиме. Описана разработанная
на языке Python программа, названная «Seismograph», позволяющая производить визуализацию как в случае 2D-, так и ЗО-расчётов. Основными отличительными особенностями «Seismograph» являются: кросс-платформенность (Microsoft Windows, Linux), широкий спектр настроек, существенно упрощающий процесс постпроцессинга и интерпретации результатов расчётов, простота расширения и адаптации под конкретные особенности регистрируемых и моделируемых процессов (например, разномасштабная амплитудность и временная продолжительность откликов). Приведены примеры использования «Seismograph» на этапе интерпретации результатов моделирования для получения практически значимых результатов. Продемонстрирована высокая востребованность описанного подхода и важность полученных с его использованием эмпирических зависимостей.
Исследование выполнено при частичной финансовой поддержке РФФИ в рамках научного проекта № 14-07-31181 мол_а и стипендии Президента РФ СП-2548.2013.5.
Литература
1. Магомедов K.M., Холодов A.C. Сеточно-характеристические численные методы. — М. : Наука, 1988. - С. 288.
2. Седов Л. И. Механика сполшной среды. — М. : Наука, 1970. — С. 492.
3. John Е. Grayson Python and Tkinter Programming. — Manning Publications, 1999. — P. 658.
4. Левянт В. В., Петров И. В., Муратов М.В. Численное моделирование волновых откликов от системы (кластера) субвертикальных макротрещин // Технологии сейсморазведки. — 2012. — № 1. — С. 5-21.
5. Муратов М. В., Петров И. Б. Расчет волновых откликов от систем субвертикальных макротрещин с использованием сеточно-характеристического метода // Математическое моделирование. — 2013. — Т. 25, № 3. — С. 89-104.
6. Левянт В. В., Петров И. В., Муратов М.В., Выко С. А. Исследование учтойчивости образования фронта рассеянных обменных волн от зоны макротрещин // Технологии сейсморазведки. — 2013. — № 1. — С. 32-45.
7. Голубев В. П., Петров И. В., Хохлов Н. И. Численное моделирование сейсмической активности сеточно-характеристическим методом // Журнал вычислительной математики и математической физики. — 2013. — Т. 53, № 10. — С. 1709-1720.
8. Голубев В. П., Квасов И. Е., Пет,ров И. Б. Воздействие природных катастроф на наземные сооружения // Математическое моделирование. — 2011. — Т. 23, № 8. — С. 46-54.
Поступим в редакцию 12.01.2014■