Программный комплекс решения обратных кинематических задач микросейсмического мониторинга Текст научной статьи по специальности «Математика»

УДК 004.021, 550.34.06

Ф. Д. Шмаков

Югорский научно-исследовательский институт информационных технологий ул. Мира, 151, Ханты-Мансийск, 628011, Россия

E-mail: sfd@uriit.ru

ПРОГРАММНЫЙ КОМПЛЕКС РЕШЕНИЯ ОБРАТНЫХ КИНЕМАТИЧЕСКИХ ЗАДАЧ МИКРОСЕЙСМИЧЕСКОГО МОНИТОРИНГА

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

Ключевые слова: обратная кинематическая задача, локация, источник, месторождение углеводородов, сейсмическая эмиссия, численные методы, алгоритм.

Введение

Локация - определение местоположения объекта посредством анализа волнового поля, отраженного от объекта или излучаемого им.

Принципы локации источников хорошо описаны в классических работах [ 1; 2]. При этом рассматриваются задачи локации в однородных средах с известными спектрами источника и помех. Примером могут служить задачи локации в гидроакустике в однородной среде с одним типом волн и известной частотой. Наиболее сложными в решении являются задачи локации источников по регистрационным записям упругих волн, распространяющихся в неоднородных средах, что обусловлено наличием различных типов волн: продольной, поперечной, Рэлея, Лява и т. д.

Определение положения источников сейсмического излучения относится к обратным кинематическим задачам. Один из способов их решения основывается на использовании методов направленного приема. К ним можно отнести сейсмоэмиссионную томографию, когда известны источник и время начала излучения волн, сейсмическую локацию бокового обзора (СЛБО). Это наиболее разработанные методы [3], использующие решение обратной кинематической задачи для определения положения источников сейсмического излучения. Данные методы позволяют производить исследование сейсмической эмиссии на разрабатываемых месторождениях УВ.

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

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

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

ISSN 1818-7900. Вестник НГУ. Серия: Информационные технологии. 2010. Том 8, выпуск 2 © Ф. Д. Шмаков, 2010

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

Технология микросейсмического мониторинга месторождений УВ основана на пассивной схеме наблюдений. Регистрация сейсмической эмиссии осуществляется на дневной поверхности над залежью УВ посредством сейсмической антенны (группы сейсмоприемников). Применяемые в технологии микросейсмического мониторинга месторождений УВ методы при обработке и интерпретации регистрационных записей сейсмической антенны позволяют:

• выделять пространственные зоны микросейсмической активности;

• анализировать изменение интенсивности излучения энергии в процессе разработки месторождения;

• анализировать связь микросейсмической активности с интенсивностью добычи / закачки флюида;

• оценивать изменение конфигурации каналов фильтрации;

• наблюдать за продвижением фронта воды от нагнетательных скважин;

• выявлять области активных разломов, зоны трещиноватости и т.п.

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

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

Рассматриваемая обратная кинематическая задача является некорректной ввиду отсутствия единственности решения в достаточно широком классе функций. Именно поэтому решение ищется в классе параметрических моделей среды. Для решения задачи с приближенными данными используется идеология поиска квазирешения методом наименьших квадратов [4].

Метод локации

В методе локации применяется обработка регистрационных записей сейсмической антенной, расположенной на дневной поверхности, с сейсмоприемниками, установленными друг от друга на удалениях не более 100 метров. Количество датчиков ограничено. Обычно сейсмическая антенна содержит от 18 до 36 сейсмоприемников с апертурой приемной группы 400-800 метров. Антенна располагается в области над исследуемой зоной предполагаемых источников. Местоположение каждого датчика сейсмической антенны определяется с помощью дифференциального С^-приемника. Регистрирующая аппаратура синхронизируется по спутниковой привязке [5].

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

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

Решение обратной кинематической задачи

Рассмотрим модель изотропной горизонтально-слоистой среды с известной зависимостью скорости волн \>(г) = \>, от глубины. Такая постановка предполагает, что на исследуемом участке месторождения УВ проведены работы вертикального сейсмического профилирова-

ния (ВСП), результатом которых является зависимость интегральных и интервальных скоростей от глубины.

По смыслу решаемых задач в технологии микросейсмического мониторинга месторождений УВ эндогенные источники можно рассматривать как точечные. Для них выполняется неравенство

2

где Ь - линейный размер источника; \ - длина волны; г0 - расстояние от источника до пункта наблюдения.

Исследование неравенства показало, что оно выполняется для частот свыше 0,1 Гц, диаметров источника менее 50 метров и расстояний от пункта наблюдений до источника более 0,5 км. Таким образом, модель точечного источника подходит для наблюдений за тектоническими процессами на месторождении УВ с помощью локальной сейсмической антенны [1].

Для изотропной горизонтально-слоистой среды время пробега можно определять как интеграл вдоль прямолинейного пути. Таким образом, определению подлежит всего четыре параметра: три пространственных координаты источника (две по латерали - хв, уй одна по глубине - и момента начала излучения волн - ¿0 или скоростной характеристики среды V;,.

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

По регистрационным записям определяем время задержки сигнала хтп между всеми парами точек наблюдения сейсмической антенны, установленной на дневной поверхности. Оценку хтп вычисляем в каждый дискретный момент времени. Время между парами точек наблюдения определяем как разницу времен пробега от источника колебаний с координатами х, У, 2 к двум точкам наблюдения с координатами хт, ут, и хп, уп, со скоростью распространения сейсмической волны :

_р(я,у =Т

/у / у тпУ 5'

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

Для всех пар точек наблюдения вычисляем взаимно корреляционные функции Стп. Используем рекуррентную формулу:

С (7 + 1 ) = С (0-- X а-к)-х а-к)-х (7 + 1)-х (7 + 1) ,

тп 4 тп 4 7 т 4 п 4 т 4 п 4 "

к

где С - вычисляемая взаимная корреляция; х , х - данные точек наблюдения, вычисления производятся в окне длиной к.

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

Р(х , у , V ) = У Т(х,у,г,у)-х 2 «С2 -»тт. (1)

^ ^ ? ✓ ^ ? з' г ^ ^^^ тп 4 ^ ? ✓ ^ ? э' г ^ тп тп х ,у ,\>

Рассмотрим метод минимизации суммы взвешенных квадратов невязок между наблюдаемыми и теоретическими значениями при определении координат источника.

Пусть а = (х, у, г, \>_) - набор параметров. Тогда минимизируемый функционал (1) пере-

Е2 2

Тт„ (а) ~%гт ' Стп ■ Если в качестве пробного значения а выбрать

т,п

о — дТ 0

а , достаточно близкое к а, так, чтобы —— почти не изменялось в интервале между а и

да.

]

а, то

О = ■

дР

- дс{

с1т1(а°) + (ак-а°к)-^(а0)

V),

(2)

лк J

где <Лтп(а) = 1тп-Ттп{а) -Стп. Отсюда следует, что поправку р = а -а, которую нужно внести в модель, можно получить из невязок ё, если (2) записать в виде

Сс1 = ССр,

(3)

м

где С = ——(а ) и О =0 .

^ тп^^ V ' тпJ .тп

В теории оценивания параметров методом наименьших квадратов уравнение (3) известно

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

Для минимизации (1) предлагаем два итерационных метода: метод наименьших квадратов, метод градиента. Пусть требуется минимизировать величину

(4)

где Тш и А - исходные данные и параметры модели соответственно. Имея начальную

пробную модель Л0, проведем линеаризацию связи между исходными данными и парамет-

гдТ \

рами модели, полагая Т(А) ~ Т(А ) + (¡а. где О =

, дА. . ,

V з )А=А

, а = А-А . Подставляя

<1 — т — Т(А '), получаем /■' ~ |с/ — Оа\ . Минимизируя /'. находим решение методом наименьших квадратов: аь—(СС)1Сс{. Процедура отыскания решения итерационным методом в случае линейности строится путем перехода ко второй модели и использования уравнения

А1 =А°+аь и проверки Т(А) ~Т(А1) + Са, где Отп =

/дТ л

тп

ч дА. , ,

V з

. Далее процесс повторяет-

ся. Если такой процесс не сходится, вторую модель можно задать в виде А1 =А° + Ка , где 0<К <1.

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

ЕдР о

-(А — А ) = 0, направ-

,■ дА. 3 3

ление быстрейшего изменения, которое нормально к этой плоскости, дается вектором а° с

/ Л А°

компонентами <-->, взятыми в точке А = А .

дА,

2

а

дР" дТ дТ

Из (4) имеем -= 2 V т -Т (А) = 2V й = А и находим, что

у 3 т-п 3 т'п

аО параллельно Ой.

Коррекция модели по методу градиента всегда будет сходящимся процессом. Однако процесс может сходиться весьма медленно. С другой стороны, в методе наименьших квадратов сходимость, если она имеется, всегда оказывается быстрой. Векторы аО и аь часто оказываются почти перпендикулярными друг другу. В этом случае выбирают среднее направление т — (СС + Е21)1Сс{. При малых е2 процесс может расходиться, а при больших

сходиться, но медленно. Оптимальное значение е2, выбранное методом подбора, обеспечит быструю сходимость.

Естественное обобщение этих двух методов - метод, согласно которому на первом этапе отыскивается зона минимума методом Монте-Карло, а затем минимизируется ¥ в окрестности зоны минимума [1].

Минимизацию функционала ^(х ух, г. , ) осуществляем методом сопряженных градиентов с учетом значения функции взаимной корреляции. Задача минимизации (1) не решается, если значения максимумов функции взаимной корреляции меньше некоторой наперед заданной величины, например 0,7 [6]. В случае, когда поле полученных решений (поле координат источников) фокусируется при уменьшении невязки, решение считаем удовлетворительным.

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

Реализация алгоритма

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

Система архивации и визуализации данных. Огромные объемы первичной сейсмической информации, собираемой при регистрации сейсмической эмиссии в полевых условиях, предполагают наличие современных средств ее сбора, передачи и обработки. Например, объем информации, собранной при пассивном наблюдении за источниками сейсмической эмиссии на месторождении «Лебяжье», составил 0,4 терабайта. Нами разработаны соответствующие средства автоматизации сбора и архивирования первичной микросейсмической информации. Система состоит из следующих программных модулей.

• Платформенно-независимый сервер реального времени для сбора и хранения полевой информации.

• Архивный клиент-визуализатор.

Программа локации источников сейсмической эмиссии.

Система реализует следующие функции:

• прием информации со станций Reftek das130-1 в режиме реального времени посредством телеметрии по стандартному протоколу RTP;

• передачу данных клиентам реального времени по своему собственному протоколу для дальнейшей обработки;

• работу с архивами данных в формате Reftek (индексирование, поиск данных в базе);

• мониторинг состояния станций в реальном времени;

• управление процессом оцифровки станций RefTек;

• визуализацию данных и их предварительную обработку;

• спектральный анализ архивных данных.

Сервер обеспечивает прием данных по телеметрической системе, реализуя стандартный протокол RTP (Ref Tek Protocol), архивирование информации в стандартном представлении, совместимом с пакетом PASSCAL, выдачу данных по запросу клиента, используя протокол TCP/IP. Сервер содержит средства управления и контроля станциями Reftek [7].

База данных первичной информации микросейсмического мониторинга месторождений УВ. Данные базы разделены на архивы данных. Архив - это данные в формате Reftek и файл собственного PASSCAL-формата, хранящий индексную информацию, создаваемую СУБД при процедуре импорта данных из архива Reftek и использующуюся для быстрого поиска нужных сейсмических данных в архиве Reftek. Файл с индексной информацией мал по объему (сотни килобайт при сотнях гигабайт «первичных» данных в архиве Reftek), вся информация из него загружается в память ЭВМ при открытии архива. Оригинальный алгоритм индексирования позволяет обеспечить быстрый доступ клиентам базы данных к архивам с временным позиционированием на трассы (с точностью до одного отчета АЦП). Алгоритм индексирования реализован в составе «Базы данных первичной информации микросейсмического мониторинга месторождений углеводородов ХМАО» [8].

В состав системы входит архивный клиент [9]. Архивный клиент-визуализатор предназначен для работы с базой данных первичной информации микросейсмического мониторинга месторождений УВ. Он позволяет:

• выбирать источник данных: архив, станцию, канал;

• выбирать временной диапазон;

• визуализировать выбранный интервал данных в окне программы и менять параметры отрисовки: уровень сигнала, удаление математического ожидания, масштаб отрисовки;

• просматривать спектр сигнала на некотором участке: размер шага и окно задается опционально;

• рассчитывать взаимную корреляцию между каналами;

• просматривать полученные графики корреляции и ковариации;

• преобразовывать данные в популярные форматы для дальнейшей обработки внешними программами.

Программа расчета положения источников сейсмической эмиссии. Алгоритм локации источников сейсмической эмиссии реализован в составе «Программы выделения зон повышенной сейсмической эмиссии по трехкомпонентным записям» [10]. Как было отмечено выше, алгоритм определения координат источников сейсмической эмиссии основан на квазиньютоновской схеме и методе сопряженных градиентов. Алгоритмы реализованы в виде динамической библиотеки для Win32 Borland Pascal Library (BPL) и оформлены согласно программному интерфейсу компоненты TQuasiNewton для языков C++ и Borland Pascal. В настоящее время существуют версии интерфейсов для Delphi 6, Delphi 7 и Delphi 2007.

Причины объединения двух алгоритмов в одну компоненту следующие: большой объем общего кода, связанного с процедурами одномерного поиска с квадратичной (SP2) и кубической (SP3) интерполяцией;

• однотипные процедуры обработки исключительных ситуаций (try ... except)',

• однотипные процедуры реализации приоритетов и возможностей фонового решения задачи;

возможность быстрого переключения методов в зависимости от специфики задачи.

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

Входными данными для программы локации источников сейсмической эмиссии являются сейсмотрассы, полученные из системы сбора. Фактически это синхронизированные данные сейсмических станций Reftek, специально подготовленные программой - архивным клиен-

том, с паспортом каждой трассы. Программа локации осуществляет расчет пространственного положения источников сейсмической эмиссии. Для каждого источника вычисляются его латеральные координаты, глубина, средняя скорость до приемника, время регистрации. Критерием выбора участков записей является невязка решения обратной задачи. На рис. 1 представлено распределение источников сейсмической эмиссии при производстве ГРП, а именно множество решений (координат источников) в проекции на дневную поверхность, полученное при обработке регистрационных записей длительностью 20 минут с количеством используемых точек наблюдения - 36: а - без учета невязки; Ъ - с невязкой меньше, чем 10~5;

с - с невязкой при минимизации функционала в пределах Ю-4 —10~5.

300т -

250 -

О | | I I I | > I I I | I " I I | ' I ^« » | > I I | | I I ■ ■ | ■ I I I |

0 50 100 150 200 250 300 350 400 450 500 т 300т -

О 50 100 150 200 250 300 350 400 450 500 т

Рис. 1. Проекция зоны ГРП на дневную поверхность.

Фокусировка источников сейсмической эмиссии при уменьшении невязки

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

Результатом обработки данных микросейсмического мониторинга являются 2В- и 30-модели распределения источников сейсмического излучения. Нами разработаны специа-

лизированные программы визуализации, которые дают возможность производить анализ распределения источников сейсмической эмиссии в пространстве месторождения УВ. Результаты работы программы визуализации изображены на рис. 2, где показано распределение источников сейсмической эмиссии при закачке жидкости в пласт в проекции на дневную поверхность, полученное при обработке регистрационных записей длительностью 4 часа. В процессе закачки жидкости в пласт отбирались события с невязкой решения обратной кинематической задачи не более 10~б.

Рис. 2. Распределение источников сейсмической эмиссии при закачке жидкости в пласт в проекции на дневную поверхность

Заключение

Разработаны быстродействующий алгоритм локации источников сейсмической эмиссии и программное обеспечение для решения обратных кинематических задач в технологии микросейсмического мониторинга месторождений УВ.

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

Данный комплекс алгоритмов и программ прошел апробацию при обработке микросейсмических данных, зарегистрированных на месторождениях Ханты-Мансийского автономного округа: «Лебяжье», «Приразломное», «Ханты-Мансийское», «Западно-Малобалыкское», «Вать-Еганское», «Потанайское».

В результате обработки данных микросейсмического мониторинга построены 2D- и 3D-модели распределения источников сейсмического излучения на указанных месторождениях.

Список литературы

1. Аки К., РичардсП. Количественная сейсмология: В 2 т. М.: Мир, 1983.

2. Бендат Дж., Пирсол А. Применения корреляционного и спектрального анализа. М.: Мир, 1983. 312 с.

3. Кузнецов О. Л., Чиркин И. А., Курьянов Ю. А., Рогоцкий Г. В., Дыбленко В. П. Сейсмо-акустика пористых и трещиноватых геологических сред. Экспериментальные исследования. М.: Гос-й науч. центр Российской Федерации - ВНИИгеосистем, 2004. Т. 2. 362 с.

4. Тихонов А. Н, Арсенин В. Я. Методы решения некорректных задач. М.: Наука, 1979. 288 с.

5. Ерохин Г. Н, Майнагашев С. М., Бортников П. Б., Кузьменко А. П., Родин С. В. Способ контроля разработки залежей углеводородов по микросейсмической эмиссии. Патент РФ №2309434.

6. Бортников П. Б., Майнагашев С. М. Обратные задачи микросейсмического мониторинга // Информационные технологии и обратные задачи рационального природопользования. Ханты-Мансийск, 2005. С. 79-83.

7. Симаков Д. Е. RTPD-J Server (КЭСР 226). Свидетельство о государственной регистрации программы для ЭВМ № 2006613623.

8. Симаков Д. Е., Шмаков Ф. Д. База данных первичной информации микросейсмического мониторинга месторождений углеводородов Ханты-Мансийского автономного округа. Свидетельство о государственной регистрации базы данных № 2009620134.

9. Симаков Д. Е. RTPD-J Archive Client (КЭСР 226). Свидетельство о государственной регистрации программы для ЭВМ № 2006613621.

10. Шмаков Ф. Д. Выделение зон повышенной сейсмической эмиссии по трехкомпонент-ным записям (ВЗПОСЭ). Свидетельство о государственной регистрации программы для ЭВМ № 2009611712.

Материал поступил в редколлегию 15.03.2010

F. D. Shmakov

PROGRAM COMPLEX OF THE SOLUTION OF AN INVERSE KINEMATICS OF MICROSEISMIC MONITORING

The solution of an inverse kinematics of a seismic wave transmission for the purpose of location of seismic emission sources directed on research of tectonic activity of hydrocarbon deposits is offered. The algorithm and the software for the solution of an inverse kinematics in methods of microseismic monitoring of hydrocarbon deposits are created. The software includes the system of data archiving and visualization, the program of computation of unknown parameters of seismic emission sources.

Keywords: inverse kinematics, location, source, hydrocarbon deposits, seismic emission, numerical methods, algorithm.