СПОСОБ ОПРЕДЕЛЕНИЯ УПРУГИХ СВОЙСТВ ГОРНЫХ ПОРОД НА ОСНОВЕ ПЛАСТОВОЙ АДАПТИВНОЙ ИНВЕРСИИ СЕЙСМИЧЕСКИХ ДАННЫХ Российский патент 2014 года по МПК G01V1/28 

Изобретение относится к способам расчета упругих свойств горных пород: продольного импеданса (Ip), поперечного импеданса (Is) и плотности (ρ) по сейсмическим данным.

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

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

Прогнозирование свойств коллекторов по данным сейсморазведки выполняется способами динамической интерпретации (ДИ), основным элементом которой является сейсмическая инверсия, осуществляющая оптимизационный подбор модели среды. На сегодняшний день практически во всех имеющихся отечественных и зарубежных системах используются способы инверсии, основанные на тонкослоистых моделях среды. В этих моделях временная мощность пластов постоянна и равна шагу дискретизации, что соответствует толщинам слоев порядка 3-4 м. Поскольку вследствие регуляризации получаемые такими способами кривые импедансов и плотностей являются сглаженными, уместно назвать эти способы термином «непрерывная инверсия». Теоретические исследования и практические примеры опробования показали, что непрерывная инверсия дает надежные оценки акустических свойств достаточно мощных (15-20 м) однородных коллекторов, тогда как тонкие пласты-коллекторы мощностью до 10 м либо не отображаются на результатах инверсии, либо их акустические свойства оцениваются с недопустимо большими погрешностями (до 30-40%).

Для обработки 2D сейсмических данных нами разработана система динамической интерпретации OTRI (в русской аббревиатуре ОТДИ), в состав которой входит оригинальная процедура акустической инверсии, основанная на применении пластовой модели среды («пластовая инверсия»). Эта технология позволяет оценивать акустические свойства пластов с высокой точностью (единицы процентов). Как показали наши исследования в случае различия акустических свойств коллекторов и вмещающих пород, применение пластовой инверсии позволяет оценивать эффективные мощности коллекторов, начиная с Н_эф=2-3 м. Данная точность достаточна для построения детальных геологических моделей, используемых не только на этапе разведки, но и на этапе эксплуатации месторождений, например при дизайне гидроразрыва пласта.

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

Выявление коллекторов и изучение распределения их емкостных свойств по площади, что необходимо для разведки и рациональной эксплуатации месторождений, осуществляется в результате применения динамической интерпретации (ДИ) сейсморазведки в комплексе с геофизическим исследованием скважин (ГИС). Основным элементом динамической интерпретации сейсморазведки является сейсмическая инверсия, осуществляющая трансформацию сейсмических волновых полей в детальные модели среды, описывающие распределение в пространстве акустических параметров: скоростей V_p и V_s продольных и поперечных волн, а также плотностей ρ. Совокупная интерпретация полученных по данным сейсморазведки прогнозных акустических параметров и имеющихся данных ГИС позволяет прогнозировать литологию пород, емкостные свойства коллекторов и характер их насыщения (газ или жидкость).

В практике интерпретации данных сейсморазведки используются разнообразные способы одномерной (акустической) и двумерной (упругой) сейсмических инверсий. Одномерная сейсмическая инверсия обрабатывает трассы суммарных сейсмических разрезов с целью получения акустических импедансов I_p=ρ*V_p и используется для оценки емкостных свойств коллекторов. Двумерная инверсия обрабатывает исходные сейсмограммы (до суммирования) продольных волн, получает все три акустических параметра и используется для прогнозирования литологического состава пород и характера насыщения коллекторов.

Все программы сейсмической инверсии, содержащиеся в широко применяемых зарубежных интерпретационных системах, основаны на тонкослоистых моделях среды, в которых минимальные мощности пластов составляют 3-4 м. Казалось бы, что именно такие программы являются наиболее эффективными для прогнозирования свойств маломощных коллекторов, имеющих подобные толщины. Однако вследствие неминуемого ограничения частотного диапазона сейсморазведки, ее реальная разрешающая способность (обычно 10-15 м) не позволяет надежно оценивать параметры отдельных пластов, слагающих используемые тонкослоистые модели [Veeken P., Da Silva. M., 2004]. Теоретические исследования и многочисленные практические эксперименты показали, что эти программы инверсии позволяют с высокой точностью (единицы процентов) определять параметры однородной группы тонких слоев, объединяющихся в более мощные пласты (от 15-20 м и более) [Buland A., Omre H., 2003, Bosch et al., 2010]. То есть такие способы инверсии дают хорошие результаты при прогнозировании свойств достаточно мощных коллекторов, но не пригодны для исследования коллекторов, мощность которых не превышает 10-15 м. На полученных практических результатах видно, что ошибки оценки акустических свойств отдельных тонких пластов с резко отличными акустическими свойствами, мощность которых не превышает 5-10 м, могут достигать 30% и более.

Выход из казалось бы неразрешимой ситуации с исследованием маломощных коллекторов состоит в использовании при выполнении сейсмической инверсии пластовых моделей среды, в которых минимальная мощность пластов соответствует реальной разрешающей способности сейсморазведки в конкретных сейсмогеологических условиях. При этом имеется возможность оценить с высокой точностью акустические параметры пластов, мощность которых составляет 15-20 м. Это оказывается достаточным для того, чтобы выявить наличие в таких пластах маломощных коллекторов (от 2-3 м), отличающихся по акустическим свойствам от вмещающих пород, а также количественно оценить их эффективные толщины Н_эф, а в некоторых ситуациях - их линейные емкости Е=Н_эф*К_п (К_п - коэффициент пористости).

Нами разработана оригинальная отечественная система одномерной инверсии ОТДИ, основанная на использовании пластовой модели среды [Кондратьев И.К. и др., 2005, Рыжков В.И. и др., 2008]. Эта система при практическом применении с использованием 2D сейсморазведки на десятках месторождений, расположенных в России и за рубежом (Казахстан, Узбекистан, Китай, Алжир), показала высокую эффективность прогнозирования маломощных коллекторов, в том числе и в сложных сейсмогеологических условиях Восточной Сибири, а также в подсолевых отложениях Прикаспия [Кондратьев И.К. и др., 2010].

Однако на сегодня существует только одномерная пластовая инверсия, применяемая для обработки данных 2D сейсморазведки. Именно из-за этого она не находит широкого применения, и нефтяные компании вынуждены использовать импортные программные продукты. Необходимо развить подход, реализующий пластовую инверсию, на вариант двумерной инверсии и сделать пластовую инверсию пригодной для обработки данных 3D сейсморазведки. Необходимо также опробовать новые варианты пластовой инверсии в основных нефтегазоносных провинциях России и разработать методические рекомендации по ее использованию для эффективного прогнозирования маломощных коллекторов.

Известен «СПОСОБ ПОИСКА, РАЗВЕДКИ, ИССЛЕДОВАНИЯ, ОЦЕНКИ И ПРОГНОЗИРОВАНИЯ РАЗРАБОТКИ ЗАЛЕЖИ И МЕСТОРОЖДЕНИЯ ПОЛЕЗНЫХ ИСКОПАЕМЫХ» (RU 2001120889, опубл. 2004, Миколаевский Э.Ю.).

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

Предлагаемый нами метод предназначен для изучения пластов малой мощности и основывается на модели физики пласта (rock physics).

Известна «АВТОМАТИЗИРОВАННАЯ СИСТЕМА ПОИСКА, РАЗВЕДКИ, ИССЛЕДОВАНИЯ, ОЦЕНКИ И ПРОЕКТИРОВАНИЯ РАЗРАБОТКИ ЗАЛЕЖИ И МЕСТОРОЖДЕНИЙ ПОЛЕЗНЫХ ИСКОПАЕМЫХ» (RU 22831U, Миколаевский Э.Ю.).

Данная система предназначена для реализации способа поиска и разведки, описанного в заявке RU 2001120889. Автоматическая обработка большого количества разнородных параметров требует мощного логического блока, который является центральным звеном данной системы. Устойчивые статистические связи разнородных параметров характерны только для массивных однородных пластов. Сложноустроенные коллекторы малой мощности не находят отклик в большинстве геофизических полей, регистрируемых на поверхности, и имеют сложные корреляционные связи геофизического параметр-петрофизического свойства.

Предлагаемый нами метод направлен на повышение разрешающей способности сейсморазведки при измерении конкретного петрофизического свойства - акустического импеданса в тонких пластах.

Известен патент CN 2577296 «UNDERGROUND MULTI-FREQUENCY ACOUSTIC COMBINED DETECTOR», HOU CHUNHUI [CN]; CHU ZEHAN [CN]; LI JIANHAO [CN].

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

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

Известен патент EP 0889331 «METHODE POUR MODELISER EN 3D L'IMPEDANCE D'UN MILIEU HETEROGENE», Grizon, Laurent, Leger Michel, Richard Vincent и др.

В решении описан метод оптимизации априорно заданной трехмерной модели импедансов, полученной в результате сейсмической интерпретации и интерполяции данных ГИС имеющихся скважин.

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

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

Известен «СПОСОБ ДЛЯ ОБНАРУЖЕНИЯ ПОТЕНЦИАЛЬНОГО ВОДНОГО ПОТОКА С МАЛОЙ ГЛУБИНОЙ ЗАЛЕГАНИЯ» (RU 2319983, за авт. Надир Дуттаб, Субхашис Маллик).

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

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

Известна международная заявка WO 2010/092084 «FINITE ELEMENT MODELLING OF BOREHOLE SEISMIC RESPONSES IN LAYERED ANISOTROPIC FORMATIONS AND ITS USE FOR ELASTIC INVERSION», JORGENSEN Ole.

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

Наша работа посвящена способу решения обратной задачи (прогноз неизвестных свойств среды на основе зарегистрированного волнового поля).

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

Для достижения поставленной цели необходимо решить следующие задачи:

a) разработка алгоритмов и методики пластовой акустической инверсии данных 3D сейсморазведки;

b) разработка алгоритмов и методики пластовой синхронной инверсии данных 3D сейсморазведки;

c) разработка методических рекомендаций оптимального комплексирования акустической и синхронной пластовой инверсии.

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

Указанный технический результат достигается за счет того, что способ определения упругих свойств горных пород на основе пластовой адаптивной инверсии сейсмических данных, характеризующийся применением пластовых моделей среды, в которых минимальные временные мощности τ_min пластов соответствуют реальной разрешающей способности сейсморазведки и геологии осадконакопления и вычисляются согласно формуле:

τ_min(мс)= $$\frac{1}{4}\ast \frac{1000}{\Delta f}$$ , где Δf - рабочая полоса частот, в которой, после всех примененных процедур обработки, уровень сигналов превышает уровень; причем в качестве исходной модели используют откорректированную геоакустическую модель импедансов соответствующей скважины, либо трассу импедансов полученного ранее сейсмоакустического разреза по секущему профилю; а при построении рабочей модели преобразуют сейсмический временной разрез в детальную пластовую модель акустических импедансов (сейсмоакустический разрезом), по которой свидетельствуют о диапазонах изменения акустических параметров слоев.

Благодаря тому что способ основан на применении пластовых моделей среды, в которых минимальные временные мощности τ_min пластов соответствуют реальной разрешающей способности сейсморазведки и геологии осадконакопления, эмпирическая формула, связывающая τ_min с качеством исходных сейсмических материалов (отношение сигнал/помеха, реальная полоса частот сейсмических сигналов), применение пластовой сейсмической инверсии обеспечивают максимальную точность оценки упругих свойств горных пород. Другие способы сейсмической инверсии основаны на применении эквидистантных моделей среды с постоянной мощностью пластов, равной интервалу дискретности Δt сейсмической записи. Реализуемое в пластовой инверсии согласование детальности получаемых сейсмоакустических разрезов с качеством исходных материалов и границ геологических пластов позволяет эффективно использовать метод для прогнозирования коллекторов в сложных сейсмогеологических условиях.

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

Осуществление изобретения

Перечень исходных данных, требуемых для реализации способа:

τ_min - минимальная временная мощность пластов получаемых моделей, соответствующая реальной разрешающей способности сейсморазведки;

$$\overline{q}$$ - исходная откорректированная пластовая модель среды, соответствующая сейсмической трассе временного разреза для пластовой акустической инверсии и сейсмограмме для пластовой синхронной инверсии;

$$t_0^k(x)$$ - линия t₀(х) для k-го реперного пласта, полученная в результате корреляции отражений от кровли соответствующего пласта на временном разрезе;

$$I{p}_k^{\min },I{p}_k^{\max }$$ - граничные значения продольных импедансов для k-го реперного пласта;

$$I{s}_k^{\min },I{s}_k^{\max }$$ - граничные значения поперечных импедансов для k-го реперного пласта;

ΔIp_max,ΔIs_max - максимальное значение изменения импедансов на одну итерацию.

Для задания τ_min используют эмпирическую формулу, характеризующую оптимальное значение этой величины по критерию точности аппроксимации «истинной» детальной модели соответствующей трассой полученного сейсмоакустического разреза:

τ_min(мс)= $$\frac{1}{4}\ast \frac{1000}{\Delta f}$$

В этой формуле Δf - рабочая полоса частот, в которой, после всех примененных процедур обработки, уровень сигналов превышает уровень.

Исходной моделью $$\overline{q}$$ является либо откорректированная геоакустическая модель импедансов соответствующей скважины, либо трасса импедансов полученного ранее сейсмоакустического разреза по секущему профилю. Отметим, что в случае невысокого качества исходных сейсмических разрезов, последняя модель также нуждается в корректировке, выполняемой в интерактивном режиме.

Главным вопросом в задании линий $$t_0^k(x)$$ является выбор реперных горизонтов, осуществляемых в соответствии со следующими условиями. Во-первых, реперные пласты должны быть распространены на всей площади и характеризоваться выдержанными упругими свойствами (обычно плащевидные глины или маркирующие карбонаты). Во-вторых, желательно, чтобы эти пласты располагались на границах (или вблизи от них) контрастирующих по акустическим свойствам мощных толщ (сотни метров). В случае выполнения первого условия, диапазоны изменения импеданса ΔIp= $$I{p}_p^{\max }-I{p}_p^{\min }$$ в реперных пластах составляют обычно 800 - 1500 м/с·г/см³(ед). Диапазоны изменения Ip в межреперных пространствах ( $$I{p}_{k,k+1}^{\min },I{p}_{k,k+1}^{\max }$$ ) определяются в результате анализа кривых Ip(H) по всем имеющимся скважинам. Величина ΔIp_max - максимально допускаемая корректировка Ip в пластах модели при одном шаге интерактивной оптимизации, задается обычно в диапазоне (0.35-0.5) от минимальной величины ΔIp (диапазонов изменения Ip в реперных пластах). Аналогично ограничения задаются для поперечных импедансов Is.

Алгоритм пластовой акустической инверсии

Постановка задачи пластовой акустической инверсии

Введем следующие обозначения:

$$\overline{s}$$ =[s₁, s₂, …, s_N]^Т - сейсмическая трасса, представленная вектором - столбцом отсчетов s(t) с шагом дискретизации Δt, в интервале t_нач, t_кон;

$$\overline{p}$$ =[Ip₁, Ip₂, …, Ip_m, t₁, t₂, …, t_m]^Т= $$\left[\overline{Ip},\overline{t}\right]$$ - произвольный вектор параметров модели, имеющий ту же структуру, что и $$\overline{q}$$ ;

$$\overline{m}(\overline{p})$$ =[m₁, m₂, …, m_N]^Т - модельная трасса для вектора $$\overline{p}$$ ;

$$\overline{y}=\overline{s}-\overline{m}(\overline{q})$$ - разностная трасса между $$\overline{s}$$ и модельной трассой для начального приближения $$\overline{q}$$ ;

$$\overline{x}=\overline{p}-\overline{q}$$ - корректирующий вектор, приводящий $$\overline{q}$$ к вектору $$\overline{p}=\overline{q}+\overline{x}$$ .

Целью акустической инверсии является преобразование сейсмического временного разреза s(t, х) в детальную пластовую модель акустических импедансов Iр(t, х), называемую сейсмоакустическим разрезом. Основой нашего алгоритма является обобщенная линейная инверсия, которая модифицирована нами для введения априорной информации о диапазонах изменения акустических параметров слоев. Ниже приводится описание математической основы модифицированного алгоритма обобщенной линейной инверсии (МОЛИ).

Базовые уравнения обобщенной линейной инверсии

Реальную трассу $$\overline{s}$$ представим в виде:

$$\overline{s}=\overline{m}(\overline{p})+\overline{n}(\overline{p})$$ , (1)

где $$\overline{n}(\overline{p})$$ - помеха, под которой понимается все то, что отличает $$\overline{s}$$ от $$\overline{m}(\overline{p})$$ при выбранной модели $$\overline{p}$$ . Если полагать, что $$\overline{p}$$ близка к $$\overline{q}$$ , можно считать, что $$\overline{m}(\overline{p})-\overline{m}(\overline{q})$$ линейно зависит от $$(\overline{p}-\overline{q})$$ . Для того чтобы выразить это математически, используют разложение $$\overline{m}(\overline{p})-\overline{m}(\overline{q})$$ в ряд Тэйлора с отбрасыванием всех членов, кроме первого, что дает:

$$\overline{m}(\overline{p})\cong \overline{m}(\overline{q})+{\displaystyle \sum _i\frac{\partial \overline{m}(\overline{q})}{\partial q_i}}(p_i-q_i)$$ . (2)

Уравнение (2) можно представить в матричной форме:

$$$$ $$\overline{m}(\overline{p})=\overline{m}(\overline{q})+A\overline{x}$$ , (3)

или:

$$\left[\begin{array}{c}{m}_1(\overline{p})\\ \\ m_2(\overline{p})\\ \\ .\\ \\ .\\ \\ m_N(\overline{p})\end{array}\right]=\left[\begin{array}{c}{m}_1(\overline{q})\\ \\ m_2(\overline{q})\\ \\ .\\ \\ .\\ \\ m_N(\overline{q})\end{array}\right]\begin{array}{c}\end{array}+\left[\begin{array}{ccccc}\frac{\partial m_1}{\partial q_1}& \frac{\partial m_1}{\partial q_2}& .& .& \frac{\partial m_1}{\partial q_M}\\ \frac{\partial m_2}{\partial q_1}& \frac{\partial m_2}{\partial q_2}& .& .& \frac{\partial m_2}{\partial q_M}\\ .& .& .& .& \\ .& .& .& .& \\ \frac{\partial m_N}{\partial q_1}& \frac{\partial m_N}{\partial q_2}& .& .& \frac{\partial m_N}{\partial q_M}\end{array}\right]\left[\begin{array}{c}{p}_1-q_1\\ p_2-q_2\\ .\\ .\\ p_M-q_M\end{array}\right]$$ (4)

Далее под $$\overline{p}$$ будем понимать оптимальное решение, которое находится в результате применения нашего алгоритма. Одно из условий оптимальности:

$${\Vert \overline{m}(\overline{p})-\overline{s}\Vert }^2=\min$$ (5)

( $$\Vert \overline{u}\Vert$$ означает норму пространства L₂, т.е. $$\Vert \overline{u}\Vert$$ ²= $${\displaystyle \sum _i{u}_i^2}$$ - сумма квадратов компонент $$\overline{u}$$ ).

Второе условие оптимальности - определенные ограничения на отклонение компонент $$\overline{p}$$ , т.е. параметров p_i, от соответствующих компонент q_i исходного приближения $$\overline{q}$$ . В результате в алгоритме МОЛИ решается задача минимизации следующего функционала:

$$F(\overline{p})={\Vert \overline{m}(\overline{p})-\overline{s}\Vert }^2+{\displaystyle \sum _i{\alpha }_i{(p_i-q_i)}^2}={Ф}_1^{}+{Ф}_2^{}$$ (6)

Подставляя $$\overline{m}(\overline{p})$$ из (3) и учитывая выражение для $$\overline{у}$$ в приведенных вначале раздела обозначениях, получим тот же функционал в другой форме:

$$Ф(\overline{х})={\Vert А(\overline{х})-\overline{у}\Vert }^2+{\displaystyle \sum _i{\alpha }_i{х}_i{}^2}={Ф}_1^{\text{'}}+{Ф}_2^{\text{'}}$$ (7)

Минимизируя Ф( $$\overline{x}$$ ), мы находим корректирующий вектор $$\overline{x}$$ , дающий оптимальное решение $$\overline{p}=\overline{q}+\overline{x}$$ . Для нахождения минимума функционала Ф( $$\overline{х}$$ ) составляем систему уравнений:

$$\frac{\partial Ф(\overline{х})}{\partial {х}_1}=0$$ , $$\frac{\partial Ф(\overline{х})}{\partial {х}_2}=0$$ , …, $$\frac{\partial Ф(\overline{х})}{\partial {х}_M}=0$$ ,

которая в матричной форме имеет вид:

$$\left[A^{Т}A+E\overline{\alpha }\right]\cdot \overline{х}={А}^{Т}\overline{у}$$ , (8)

где $$Е\overline{\alpha }$$ - диагональная матрица, в которой отличны от нуля лишь диагональные элементы ( $${\alpha }_1$$ , $${\alpha }_2$$ , …, $${\alpha }_M$$ ).

Расчет демпфирующих множителей $${\alpha }_i$$

Основной особенностью модифицированного алгоритма инверсии является способ задания демпфирующих множителей $${\alpha }_i$$ , задачей которого является не позволить параметрам $${р}_i$$ оптимальных моделей выйти далеко за границы установленных диапазонов. Как следует из (6), (7), величины $${\alpha }_i$$ регулируют степень свободы изменения $$q_i$$ , т.к. чем больше $${\alpha }_i$$ , тем сильнее вклад корректировки $$q_i$$ , т.е. величины $$x_i^2={(p_i-x_i)}^2$$ , в функционалы, которые необходимо минимизировать. Величины $${\alpha }_i$$ выбираются в МОЛИ, исходя из условий:

А. Оба слагаемых в функционалах (6), (7) должны быть примерно одинаковы по величине.

В. Величины $${\alpha }_i$$ должны быть тем больше, чем ближе $$q_i$$ к заданным ограничениям $$q_{гр}$$ .

Предварительно заметим, что в (6) и (7) Ф $${}_1^{\text{'}}$$ =Ф₁ и Ф $${}_2^{\text{'}}$$ =Ф₂, т.к. (7) получено из (6) путем простой замены буквенных обозначений.

Для того чтобы определить $${\alpha }_i$$ , постараемся вначале определить величину $${Ф}_1^{}={\Vert \overline{m}(\overline{p})-\overline{s}\Vert }^2$$ . Она заранее не известна, т.к. $$\overline{p}$$ - это модель, которая получится после выполнения соответствующей итерации в процессе оптимизации. Однако можно предположить, что энергия ошибки Ф₁ после оптимизации будут представлять некоторую часть Cn(Cn<1) энергии начальной ошибки $${\Vert \overline{m}(\overline{q})-\overline{s}\Vert }^2$$ , где n - номер итерации. Предполагая, что амплитуды ошибки (точнее, помехи $$\overline{n}$$ в формуле (1)) в результате оптимизации в среднем уменьшатся в 2 раза на первом шаге итерации, задаем Сn=0,25. Поскольку на последующих шагах итерации, по мере приближения к истинной модели, уменьшение ошибки будет замедляться, задаем следующую последовательность величин Сn:

С₁=0,25; С₂=0,5; С₃=0,7. (9)

Тогда величину Ф₁ можно оценить формулой:

$${Ф}_1=Сn\cdot {\Vert \overline{m}(\overline{q})-\overline{s}\Vert }^2=Cn\cdot {\Vert \overline{y}\Vert }^2$$ . (10)

В предыдущем разделе говорилось о том, что оптимизация выполняется на каждом шаге итерации последовательно: вначале по вектору $$\overline{t},$$ затем по вектору $$\overline{Ip}$$ . Соответственно, по-разному оцениваются величины $${\alpha }_i$$ . Описание способов получения этих оценок начнем с процесса оптимизации акустических параметров $$\overline{Ip}$$ (далее для упрощения обозначений - величина $$q_i$$ ).

Вначале рассмотрим частный случай, когда значения $${х}_i$$ корректирующего вектора $$\overline{x}$$ одинаковы по модулю, а именно:

$$\left|x_i\right|=x_{ср}=0.5\Delta q_{\max }$$ . (11)

В данном случае Δq_max=ΔIp_max - максимальное изменение $$q_i$$ на 1 итерацию

Тогда, согласно условию А, требуя равенство Ф₁ и Ф₂ в формуле (7) и учитывая (10), получим одинаковые значения $${\alpha }_i$$ = $${\alpha }_{ср}$$ в виде:

$${\alpha }_{ср}=\frac{Сn\cdot {\Vert \overline{y}\Vert }^2}{M\cdot x_{ср}^2}$$ , (12)

где М - число слоев в модели.

Теперь перейдем к общему случаю, задавая $${\alpha }_i$$ в соответствии с условием В, ориентируясь на пример конкретного положения оптимизируемого параметра $$q_i$$ на оси q (см. Фиг.1).

На Фиг.1:

$$q_i$$ - начальное приближение параметра Iр в i-ом слое,

$$q_{гр}^{\min ,\max }$$ - $$I_p^{\min }$$ , $$I_p^{\max }$$ для любого, в т.ч. нереперного слоя,

$$\Delta q_i=\min \left\{(q_i-q_{гр}^{\min }),\begin{array}{c}\end{array}(q_{гр}^{\max }-q_i)\right\}$$ , (13)

δ - минимальное допускаемое расстояние между $$q_i$$ и любым из $$q_{гр}$$ , обеспечивается путем редактирования $$\overline{q}$$ (см. ниже), определяется формулой

δ= $$\frac{1}{6}\Delta q_{\max }$$ , (14)

где $$\Delta q_{\max }$$ - максимальное изменение Ip на 1 итерацию.

Теперь, для выполнения условия В, зададим $${\alpha }_i$$ со структурой, аналогичной (12), но с учетом реального положения $$q_i$$ :

$${\alpha }_i=\frac{\beta \cdot Cn\cdot {\Vert \overline{y}\Vert }^2}{M\cdot \Delta q_i^2}$$ . (15)

В выражении (15) β находится путем тестирования при обработке конкретных реальных материалов по отдельным профилям, исходя из условия обеспечения устойчивости результата при внесении вариаций в исходную модель $$\overline{q}$$ (β=1 ÷ 4).

Посмотрим на отдельных примерах, как меняется $${\alpha }_i$$ в зависимости от положения $$q_i$$ . При $$q_i$$ = $$q_{ср}^{}$$ , $$\Delta q_i=0.5(q_{гр}^{\max }-q_{гр}^{\min })$$ - мы имеем минимальное $${\alpha }_i$$ . В ситуации, когда $$\Delta q_i=x_{ср}=0.5\Delta q_{\max }$$ (полагаем, что это наиболее вероятная ситуация), $${\alpha }_i$$ = $${\alpha }_{ср}$$ (для выбранного фиксированного β). При минимальном $$\Delta q_i$$ =δ (см. (14)), $$\Delta q_i=\frac{1}{3}{х}_{ср}$$ , следовательно, $${\alpha }_1$$ будет примерно на порядок больше, чем $${\alpha }_{ср}$$ .

Рассмотрим теперь способ определения демпфирующих множителей $${\alpha }_i$$ при расчете корректирующих значений Δt_i вектора $$\overline{t},$$ входящего в начальное приближение $$\overline{q}=\left[\overline{Ip},\overline{t}\right]$$ . Учитывая, что Δt_i не должен быть велик изначально (Δt_i ≅ Δt, Δt - интервал дискретности), полагаем, что все $${\alpha }_i$$ одинаковы и, по аналогии с (12), выражаются формулой:

$${\alpha }_i={\alpha }_{ср}=\frac{Cn\cdot {\Vert \overline{y}\Vert }^2}{M\cdot {\theta }^2}.$$ (16)

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

$$\theta =\left[\frac{{\tau }_{\min }-1}{3}\right]+1$$ , (17)

где [] означает целую часть.

Рассчитанные в результате решения уравнения (8) величины Δt_i (компоненты вектора $$\overline{x}$$ являются рациональными числами, не кратными Δt. Величины Δt_i ограничиваются согласно условию

$$\left|\Delta t_i\right|\le 0.4{\tau }_{\min }$$ , (18)

после чего они добавляются в значения t_i исходной модели $$\overline{q}$$ для получения оптимизированной модели $$\overline{р}$$ . Таким же образом при оптимизации Iр, рассчитанные в результате решения (8) значения ΔIp_i, прибавляются к исходным значениям Ip_i вектора $$\overline{q}$$ .

Редактирование вектора $$\overline{q}$$ , поступающего на оптимизацию

Учитывая итерационность оптимизированного процесса, полученные после некоторого шага оптимизации векторы $$\overline{р}$$ подаются вновь на вход процесса оптимизации, т.е. становятся вновь начальным приближением $$\overline{q}$$ . Все начальные приближения, в том числе и исходные модели, подвергаются редактированию перед процессом оптимизации.

Редактирование вектора $$\overline{t}$$ начального приближения $$\overline{q}$$

Компоненты вектора $$\overline{t}$$ - величины t_i редактируются либо путем сдвига по оси t, либо в результате объединения пластов, при этом:

а) для реперных пластов с заданными временами t_k (k=1, 2, …, K) допускается объединение только с подстилающими пластами;

б) при объединении пластов скорость нового пласта рассчитывается как средневзвешенная скоростей компонент с весами, равными временной мощности τ_i объединяемых пластов.

Редактирование вектора $$\overline{t}$$ проводится следующим образом:

- проверяется выполнение условия τ_i>Δt, если нет, то два (или более) соседних пластов объединяются до выполнения этого условия;

- границы всех пластов выставляются на ближайшие точки шкалы дискретности;

- пласты с временной мощностью τ_i<τ_min объединяются с одним из двух соседних пластов, у которого Ip ближе к Ip_i (в отношении реперных пластов смотри а)).

Редактирование вектора $$\overline{Ip}$$ начального приближения $$\overline{q}$$

Выполняется проверка: если $$I{p}_i>I{p}_i^{\max }-\delta$$ , то задаем

$$I{p}_i=I{p}_i^{\max }-\delta$$ ,

а где $$I{p}_i^{\max }$$ - граничное значение Ip для i-го пласта, δ задано формулой (14). Аналогично, если $$I{p}_i<I{p}_i^{\min }+\delta$$ , то задаем

$$I{p}_i=I{p}_i^{\min }+\delta$$ .

Редактирование вектора $$\overline{Ip}$$ в окончательной модели $$\overline{p}$$

Окончательную модель $$\overline{р}$$ получаем после выполнения последнего шага итерации. Редактируются значения Ip_i согласно условиям:

если $$I{p}_i>I{p}_i^{\max }$$ , то задаем

$$I{p}_i=I{p}_i^{\max }$$ ;

если $$I{p}_i<I{p}_i^{\min }$$ , то задаем

$$I{p}_i=I{p}_i^{\min }$$ .

Алгоритм пластовой синхронной инверсии

Алгоритм синхронной инверсии является некоторым обобщением алгоритма акустической инверсии с увеличением размерности исходных данных (сейсмические трассы) и параметров модели. В частности теперь вектор $$\overline{s}$$ - сейсмограмма - описывается выражением

$$\overline{s}={\left[\overline{s_1},\overline{s_2}\mathrm{,...,}\overline{s_L}\right]}^{\begin{array}{c}\end{array}T}$$ , (19)

где компоненты $$\overline{s}$$ представляют собой трассы s(t, ξ₁), s(t, ξ₂), …, s(t, ξ_L). Соответственно, развернутое выражение для начального приближения $$\overline{q}$$ имеет вид (М-слойная модель)

$$\overline{q}$$ =[Vp₁, Vp₂, …,Vp_M, Vs₁, Vs₂, …, Vs_M, ρ₁, ρ₂, …, ρ_M]=

=[q₁, q₂, …, q_M, q_M+1, q_M+2, …, q_2M, q_2M+1, q_2M+2, …, q_3M]. (20)

Вопрос о структуре вектора $$\overline{q}$$ требует специального объяснения. Синхронная пластовая инверсия выполняется как второй этап после проведения акустической инверсии. При этом временные мощности пластов Δt_i в точке х профиля отсчитываются с сейсмоакустических разрезов как фиксированные параметры, не требующие оптимизации. Это объясняет отсутствие вектора $$\overline{t}$$ в структуре $$\overline{q}$$ . Теперь надо объяснить, откуда берутся значения Vp, Vs, ρ пластов в начальном приближении. Для получения приближенных оценок этих параметров исследуемая пачка пластов разбивается на несколько интервалов с близкими литофизическими свойствами (карбонатные, терригенные - преимущественно глинистые или преимущественно песчанистые). Для каждого из этих интервалов, используя данные ГИС по региону, где находится площадь исследования, определяют средние величины плотностей ρ, которые приписываются всем пластам интервала. С полученного ранее сейсмоакустического разреза отсчитывают значения Ip пластов и вычисляют Vp=Ip/ρ. Также по анализу имеющейся геолого-геофизической информации для каждого из выделенных интервалов задают величины γ=Vs/Vp, что позволяет рассчитать по значениям Vp значения Vs для каждого пласта модели.

По аналогии с акустической инверсией, далее нам требуется рассчитать модельную сейсмограмму $$\overline{m}(\overline{q})$$ ,

$$\overline{m}(\overline{q})={\left[\overline{m_1}(\overline{q}),\overline{m_2}(\overline{q})\mathrm{,...,}\overline{m_L}(\overline{q}),\right]}^{\begin{array}{c}\end{array}T}$$ , (21)

трассы которой m(t, ξ₁), m(t, ξ₂), …, m(t, ξ_L), рассчитываются для тех же удалений ξ, что и в соответствующих трассах реальной сейсмограммы $$\overline{s}$$ .

Для расчета $$\overline{m}(\overline{q})$$ потребуется форма сейсмического сигнала h(t, ξ). Для ее нахождения вначале определяет форму h(t, 0) на нулевом удалении взрыв-прибор. Для этого рассчитывают трассу коэффициентов отражений r(t, 0) по модели Ip(t), выбираемой с сейсмоакустического разреза Ip(t, x) в нужной точке х профиля. Сигнал h(t, 0) определяется по известным формулам как оптимальный винеровский фильтр, преобразующий r(t, 0) в трассу s(t, ξ₁), где ξ₁ - минимальное удаление в сейсмограмме s(t, ξ). Затем, по известным из обработки кинематическим поправкам для кровли Δt_кр(ξ) и подошвы Δt_под(ξ) исследуемой пачки слоев, растягиваем сигналы h(t, 0) согласно формуле:

$$h(t,\xi )=h\left[\frac{\Delta t_{под}(\xi )}{\Delta t_{кр}(\xi )}t\mathrm{,0}\right]$$ . (22)

Учитывая сложность и времяемкость вычислительного процесса синхронной инверсии, мощность исследуемой пачки слоев ΔНобщ ограничивают сравнительно небольшим интервалом, включающим основные продуктивные пласты. При этом ΔНобщ обычно на порядок меньше глубины Н до кровли пачки. В этой ситуации можно пренебречь преломлением лучей внутри пачки и считать, что угол падения θ одинаков для всех плоскопараллельных границ пластов, слагающих исследуемую пачку. При этом реальное двойное время пробега волн в i-ом пласте выражается формулой:

$$\Delta t_i(\theta )=\frac{\Delta t_i(0)}{\cos \theta }$$ , (23)

где Δt_i(0) - временная мощность пласта в модели Ip(t), выбираемой из сейсмоакустического разреза. Как сказано в разделе 1.1, в исходную сейсмограмму $$\overline{s}$$ при обработке вводились статические и кинематические поправки, приводящие годографы отражений от каждой i-ой границе к константе: t_i(ξ)=const=t_i. Это условие должно выполняться и для модельной сейсмограммы $$\overline{m}(\overline{q})$$ . При этом Δt_i(θ) заменяются на Δt_i(0), тогда используемую для расчета $$\overline{m}(\overline{q})$$ импульсную трассу на удалении взрыв-прибор ξ можно определить выражением:

$$r(t,\theta )={\displaystyle \sum _i{r}_i(\theta )\cdot \delta (t-t_i)}$$ , (24)

где r_i(θ) - коэффициент отражения от кровли i-го пласта при падении луча под углом θ,

t_i - время прихода отражения от кровли i-го пласта в модели Ip(t),

δ(t) - единичный импульс.

При обычных для практики небольших углах наклона β отражающих границ (до 10°) угол θ определяется по формуле

$$\theta (\xi )=arctg\frac{\xi }{t_1\cdot V_{ср}}$$ , (25)

где V_ср - средняя скорость до кровли пачки, определяется в процессе предшествующих структурных построений.

Коэффициент отражения R_i(θ) описывается формулой

$$Ri(\theta )=\frac{1}{2}\left(\frac{\Delta V_p}{{\overline{V}}_p}+\frac{\Delta \rho }{\overline{\rho }}\right)+\left(\frac{\Delta V_p}{2{\overline{V}}_p}-4\frac{V_s^2}{V_p^2}\frac{\Delta V_s}{{\overline{V}}_s}-2\frac{V_s^2}{V_p^2}\frac{\Delta \rho }{\overline{\rho }}\right){\sin }^2\theta +\frac{1}{2}\frac{\Delta V_p}{{\overline{V}}_p}{\sin }^2\theta {\tan }^2\theta$$ (26)

где $$\Delta V=V(t_i)-V(t_{i-1}),\overline{V}=\frac{1}{2}\left(V(t_i)+V(t_{i-1})\right)$$

После этих промежуточных вычислений, трассы m(t, ξ) - компоненты вектора $$\overline{m}(\overline{q})$$ (21), рассчитываются по формуле:

m(t, ξ)=h(t, ξ) * r[t, θ(ξ)]. (27)

Введем обозначения:

$$\overline{y}=\overline{s}-\overline{m}(\overline{q})$$ - разность между исходной сейсмограммой и модельной сейсмограммой для параметров исходной модели $$\overline{q}$$ ,

$$\overline{x}=\overline{p}-\overline{q}$$ - корректирующий вектор, приводящий $$\overline{q}$$ к оптимальному решению $$\overline{p}=\overline{q}+\overline{x}$$ .

Применяя в случае синхронной инверсии тот же метод обобщенной линейной инверсии, как и для акустической инверсии, корректирующий вектор $$\overline{x}$$ находится в результате решения системы уравнений (8), имеющий вид:

$$\left[A^{Т}A+E\overline{\alpha }\right]\cdot \overline{х}={А}^{Т}\overline{у}$$ (28)

В данном случае матрица А представляется в виде:

$$А\begin{array}{c}\end{array}=\begin{array}{c}\end{array}\left[\begin{array}{cccccc}\frac{\partial \overline{m_1}(\overline{q})}{\partial q_1}& \frac{\partial \overline{m_1}(\overline{q})}{\partial q_2}& .& .& .& \frac{\partial \overline{m_1}(\overline{q})}{\partial q_{3M}}\\ \frac{\partial \overline{m_2}(\overline{q})}{\partial q_1}& \frac{\partial \overline{m_2}(\overline{q})}{\partial q_2}& .& .& .& \frac{\partial \overline{m_2}(\overline{q})}{\partial q_{3M}}\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ \frac{\partial \overline{m_L}(\overline{q})}{\partial q_1}& \frac{\partial \overline{m_L}(\overline{q})}{\partial q_2}& .& .& .& \frac{\partial \overline{m_L}(\overline{q})}{\partial q_{3M}}\end{array}\right]$$ (29)

В этом выражении $$\frac{\partial \overline{m_{\ell }}(\overline{q})}{\partial q_i}$$ - частная производная $$\ell$$ -ой трассы сейсмограммы $$\overline{m}(\overline{q})$$ по i-ому параметру вектора $$\overline{q}$$ (20). По аналогии со случаем акустической инверсии, она представляется в виде:

$$\frac{\partial \overline{m_{\ell }}(\overline{q})}{\partial q_i}=\frac{\Delta r_{i\text{'}}[\theta ({\xi }_{\ell })]\cdot h(t-t_{i\text{'}},{\xi }_{\ell })+\Delta r_{i\text{'}+1}[\theta ({\xi }_{\ell })]\cdot h(t-t_{i\text{'}+1},{\xi }_{\ell })}{\Delta q_i}$$ , (30)

где i' - номер слоя (i'=1, 2, …, M), соответствующего параметру q_i - см. (20),

Δr_i', Δr_i'+1 - изменения коэффициентов отражения от кровли и подошвы слоя i' при замене q_i на q_i+Δq_i,

t_i', t_i'+1 - положение на оси t кровли и подошвы слоя i' в модели Ip(t), полученной в результате акустической инверсии.

Для составления уравнения (28), так же, как и в случае акустической инверсии, нужно задать демпфирующие множители $${\alpha }_i$$ . Они определяются точно таким же методом, как и для акустической инверсии. По аналогии с (15), в данном случае выражение для $${\alpha }_i$$ можно представить в виде:

$${\alpha }_i=\frac{\beta \cdot Cn\cdot {\Vert \overline{y}\Vert }^2}{3M\cdot \Delta q_i^{}}$$ , (31)

где $$\Delta q_i=\min \left(\left|q_i-q_i^{\min }\right|,\left|q_i-q_i^{\max }\right|\right)$$ , а граничные значения для q_i заданы в исходной информации.

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

τ_min(мс)= $$\frac{1}{4}\ast \frac{1000}{\Delta f}$$ , где Δf - рабочая полоса частот. В качестве исходной модели используют откорректированную геоакустическую модель импедансов соответствующей скважины, либо трассу импедансов полученного ранее сейсмоакустического разреза по секущему профилю. При построении рабочей модели преобразуют сейсмический временной разрез в детальную пластовую модель акустических импедансов (сейсмоакустический разрезом), по которой свидетельствуют о диапазонах изменения акустических параметров слоев. Технический результат - повышение точности оценки упругих свойств горных пород. 1 ил.

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