Moscow, Russian Federation
Russian Federation
VAK Russia 1.2.1
VAK Russia 1.2.2
This paper addresses the current challenge of improving the reliability of seismic data interpretation by transition from the restoration of only elastic parameters to direct lithofacies prediction. Special attention is paid to the promising method of joint impedance and facies inversion (lithologic inversion), which integrates a priori facies data directly into the inversion process, ensuring consistent geological modeling. Since the synchronous inversion stage is the most computationally intensive part of the lithologic inversion algorithm, the aim of the study is to conduct a comparative analysis of modern inversion algorithms and select the most efficient one for subsequent implementation within this method.
imultaneous inversion, lithofacies prediction, seismic interpretation, geological modeling, inversion algorithms.
Введение
Современный этап интерпретации данных сейсмической разведки требует применения комплексных алгоритмов, обеспечивающих согласованные и геологически достоверные результаты. Ключевой вызов возникает на этапе прогноза геологических свойств: сегодня недостаточно только восстановления упругих параметров (акустического импеданса, плотности и др.), необходим прямой прогноз литофаций для построения детальных моделей. Несмотря на наличие ряда высокоэффективных алгоритмов синхронной инверсии (СИ) [1–4], они, как правило, решают лишь первую часть этой задачи. Полученные кубы упругих параметров требуют дополнительного этапа классификации для перехода к фациям, что вносит субъективность и увеличивает неопределенность. Таким образом, актуальной проблемой является разработка методов, интегрирующих инверсию и классификацию в единый процесс для повышения достоверности интерпретации.
В этой связи особый интерес представляет метод совместной инверсии на импеданс и фации [5, 6]. Его ключевое преимущество — интеграция априорных данных о фациях из скважин (РИГИС) непосредственно в инверсионный процесс. Это позволяет одновременно уточнять упругие свойства и прогнозировать литофациальную модель, что существенно повышает ее геологическую точность и детальность. Рассматриваемый алгоритм назван авторами литологической инверсией.
Поскольку этап синхронной инверсии является вычислительно наиболее затратным звеном в алгоритме литологической инверсии, от его эффективности напрямую зависит практическая применимость всего метода.
Цель данного исследования — провести сравнительный анализ современных алгоритмов синхронной инверсии и обосновать выбор наиболее эффективного из них для последующей реализации в рамках перспективного метода литологической инверсии.
1. Математическое описание алгоритмов синхронной инверсии
В рамках данного исследования были программно реализованы и протестированы три различных подхода к решению задачи синхронной сейсмической инверсии. Их краткое описание представлено ниже.
1.1. Стохастический адаптивный градиентный спуск (NADAM) с гибридной целевой функцией
Метод использует целевую функцию F в виде комбинации среднеквадратической ошибки (MSE) между реальными и синтетическими трассами, антикросс-корреляции χ между реальными и синтетическими трассами и регуляризирующего функционала θ (1):
F = w · MSE + (1 - w) · χ + θ, (1)
где w — управляющий вес от 0 до 1.
Каждая синтетическая трасса вычисляется как свертка трассы коэффициентов отражения r(t) с трассой импульсов w(t): sсинт.(t) = r(t) × w(t).
В качестве уравнения для коэффициента отражения r(t) используется уравнение Нотта — Цеппритца. Целевая функция минимизируется с помощью метода Nesterov ADAM [7].
1.2. Линеаризация задачи инверсии через аппроксимацию коэффициента отражения и решение системы линейных алгебраических уравнений (СЛАУ)
При построении математической модели алгоритма использовались положения из диссертации [1] и книги [4].
Для каждого угла θi записывается аппроксимация коэффициента отражения через линейную комбинацию (2):
ri = D · (Ai · x + Bi · y + Ci · z), (2)
где D — матрица Теплица, соответствующая оператору первой разностной производной, Аi, Bi, Ci — диагональные квадратные матрицы с коэффициентами, зависящие от угла θi и отношения Vs/Vp.
Выбрана аппроксимация по формуле Аки — Ричардса [3] вида (4):
,
. (3)
Тогда матрица оператора инверсии представляет собой матрицу (4):
. (4)
В итоге получаем СЛАУ Sobs = Gm, которая решается методом наименьших квадратов с добавлением регуляризующих компонентов. Обновление модели проводится через поправки с помощью СЛАУ (5):
. (5)
1.3. Градиентный метод оптимизации с линеаризацией целевой функции
Основные предположения для реализации метода взяты из источников [2] и [4].
Каждое значение коэффициента отражения вычисляется исходя из значения физических свойств в текущей точке и в следующей:
ri = f(Vpi, Vpi+1, Vsi, Vsi+1, ρi, ρi+1).
Минимизируется функционал вида (6):
. (6)
Для трассы коэффициентов отражения выполним разложение в ряд Тейлора до первого члена: r(m) = r(m0) + FΔm, где r(m) — трасса коэффициентов отражения, r(m0) — начальная трасса коэффициентов отражения, построенная по низкочастотной модели, F — якобиан, Δm — вектор поправок физических свойств.
В качестве уравнения для коэффициента отражения используется уравнение Нотта — Цеппритца.
Подставив разложение в (6), получим выражение (7):
. (7)
Приращение m вычисляется методом Nesterov ADAM [7]:
, (8)
,
где — шаг градиентного спуска, r — вектор регуляризационных поправок.
Предложенный алгоритм допускает гибкую адаптацию под различные параметризации целевой функции. В частности, при переходе от параметров (Vp, Vs, Rho) к (AI, SI, Rho) с использованием аппроксимации Фатти для коэффициентов отражения, точность решения мало изменяется, но это дает небольшую прибавку к скорости вычислений.
2. Тестирование алгоритмов синхронной инверсии
Для тестирования алгоритмов была создана синтетическая модель на основе скважинных данных с известными упругими свойствами (Vp, Vs, Rho). Путем свертки рассчитанных коэффициентов отражения (с использованием полного уравнения Нотта — Цеппритца) с импульсом Риккера был получен угловой сейсмический объем, который в дальнейшем использовался как входные данные для инверсии. Тестирование проводилось на углах до 35 градусов, где ошибки, связанные с линеаризацией, минимизированы. Для низкочастотной модели (НЧМ) «истинные» параметры были зашумлены и сглажены фильтром Баттерворта. Сравнительные результаты по точности восстановления параметров и времени вычислений отражены на рисунке 1 и в таблице 1. На рисунке 2 представлены кроссплоты для пар параметров (vp, vs) и (vp, rho) для скважинных данных и рассчитанных по каждому методу.
Таблица 1. Усредненные значения коэффициентов корреляции Пирсона, рассчитанные для каждого из трех упругих параметров (Vp, Vs, Rho), корреляции между синтетическими и истинными сейсмическими трассами, скорость работы алгоритма (трасс/сек.)
|
Метод решения |
Средняя корреляция по параметрам среды |
Средняя корреляция с сейсмическим угловым стеком |
Скорость алгоритма (трасс/сек.) |
|
NADAM с гибридной целевой функцией (1.1) |
0.86 |
0.96 |
0.39 |
|
Аппроксимация Аки — Ричардса и решение СЛАУ (1.2) |
0.84 |
0.97 |
7.09 |
|
Градиентный метод оптимизации с линеаризацией целевой функции через уравнение Нотта — Цеппритса (1.3) |
0.97 |
0.99 |
13.56 |
Рисунок 1. Результаты расчета синхронной инверсии по методам 1.1–1.3 (синий цвет) по сравнению с истинными значениями в «скважине» (оранжевый цвет) и с графиком НЧМ (черный цвет). Коэффициенты корреляции Пирсона, корень из среднего квадратичного отклонения (СКО) и относительное СКО между модельной скважиной и рассчитанной трассой по каждому из трех упругих свойств
Рисунок 2. Кроссплоты (vp, vs), (vp, rho) для модельной скважины и результатов инверсии по методам 1.1–1.3
Метод 1.3, основанный на линеаризации с использованием уравнения Нотта — Цеппритца, продемонстрировал наилучший баланс между точностью и временем расчета. Как видно на профилях упругих свойств в зависимости от времени (рис. 1) и кроссплотах (рис. 2), именно этот подход обеспечил одновременное восстановление упругих параметров.
3. Сейсмическая инверсия с привлечением данных РИГИС
Алгоритм литологической инверсии построен на фундаментальных принципах, изложенных в [5, 6], и представляет собой реализацию итерационной процедуры EM (Expectation — Maximization). Итерационная схема включает в себя следующие этапы:
-
Предварительное построение НЧМ каждого тренда фации на основе данных ГИС/РИГИС. В работе [6] предлагается линейно связывать каждое упругое свойство по времени через полином 1 степени (9):
. (9)
-
Начальная инициализации полей упругих свойств и поля фаций.
-
Построение литологии через нахождение вероятностей каждой фации (Expectation шаг) через формулу наивного Байесовского классификатора (10):
, (10)
где m = [mAI, mSI, mRho] — вектор текущих значений упругих свойств в точке (акустический импеданс, сдвиговой импеданс, плотность), полученных на текущей итерации инверсии;
μp,f(t) — трендовое (ожидаемое) значение свойства для фации в момент времени t, которое вычисляется по линейному тренду из пункта I;
σp,f — стандартное отклонение свойства p внутри фации, оцененное по скважинным данным. Оно характеризует разброс значений свойства вокруг тренда для данной фации.
-
Пересчет полей упругих свойств по данным классификации и их обработка. Вероятность каждой фации в точке умножается на значение из тренда по каждому свойству по формуле (11). Таким образом связываются поля вероятностей и поля упругих свойств:
. (11)
-
Извлечение информации из сейсмических данных посредством синхронной инверсии (Maximization шаг). В работе [6] выбрана аппроксимация Фатти для уравнения коэффициента отражения, поэтому для использовании в методе литологической инверсии адаптируем алгоритм из пункта 1.3 под параметризацию (AI, SI, Rho).
Этапы I–V — это одна полноценная итерация алгоритма литологической инверсии. Каждая итерация улучшает результат предыдущей.
4. Тестирование алгоритма литологической инверсии
Метод был протестирован на синтетической модели. В качестве исходных данных использована модельная скважина со слабым уровнем шума. На ее основе для каждой фации построены полиномиальные тренды. Начальные поля упругих свойств получаются из сглаженных фильтром Баттерворта данных ГИС. На рисунке 3 показаны тренды, построенные для каждой фации.
Рисунок 3. Тренды (AI, SI, Rho) для каждой фации
С помощью библиотеки Synthoseis [8] сгенерирован куб с литологией, приближенной к реальной геологии. С помощью свертки с импульсом Риккера получен сейсмический угловой стек.
Расчет разреза размером 324 × 834 занял 38 секунд и 10 итераций. На рисунке 4 показан истинный литологический разрез и рассчитанный.
Рисунок 4. Истинный разрез вероятностей коллектора (слева), полученный из synthoseis, и рассчитанный (справа) через метод литологической инверсии
Как видно из рисунка 4, алгоритм восстанавливает общую структуру истинного разреза, но слабо выделяет мощные слои коллектора. Этот артефакт исправляется повышением числа итераций для алгоритма синхронной инверсии.
Коэффициенты корреляции полей упругих свойств, мера f1-score для коллектора/неколлектора, среднее относительное СКО между сейсмическими и синтетическими данными для всего разреза представлены в таблице 2.
Таблица 2. Метрики расчета синтетического разреза
|
Метрика |
Значение |
|
CorrAi |
0.8811 |
|
CorrSI |
0.7147 |
|
CorrRho |
0.6126 |
|
F1-scoreнеколлектор |
0.9361 |
|
F1-scoreколлектор |
0.6599 |
|
Относит. СКО между синтетикой и сейсмикой по итерациям литологической инверсии |
0 : 1.0071 1 : 0.0604 2 : 0.0290 3 : 0.01758 4 : 0.0113 5 : 0.0078 6 : 0.0056 7 : 0.0042 8 : 0.0034 9 : 0.0027 10 : 0.0023 |
Из таблицы 2 можно сделать вывод, что алгоритм сходится к сейсмическим данным. Чем больше число итераций литологической инверсии, тем точнее литологическая картина. Коэффициенты корреляции в таблице 2 значительно ниже, чем в пункте 1.3, что объясняется упрощенной регуляризацией синхронной инверсии рассматриваемого алгоритма. Коэффициенты регуляризации фиксированы на протяжении всех итераций и могут быть неоптимальными, что открывает направление для дальнейших исследований по их адаптивному подбору.
Вывод
По результатам сравнительного анализа трех реализаций синхронной инверсии выбран градиентный метод оптимизации с линеаризацией целевой функции. Он обеспечивает наилучший баланс точности и скорости. Его использование в методе литологической инверсии полностью оправдано.
Разработанный на этой основе алгоритм продемонстрировал хорошее соответствие синтетической модели. Дальнейшие исследования будут направлены на апробацию алгоритмов на реальных геолого-геофизических данных.
1. Li Q. Development of noise-immune algorithms for dynamic seismic inversion [Doctoral dissertation] / Li Q. — Moscow, 2017.
2. Liu J. Seismic simultaneous inversion using a multi-damped subspace method / Liu J., Wang Y. // Geophysics. — 2020. — No. 85(1). — Pp. R1–R10.
3. Niu L. Linearized AVO inversion with modified Aki-Richards’ approximation / Niu L., Geng J., Wu X. // 80th EAGE Conference & Exhibition, Copenhagen, Denmark, 2018, June. — Th P1 12.
4. Wang Y. Seismic inversion: Theory and applications / Wang Y. — John Wiley & Sons, 2017.
5. Barajas-Olalde C. Joint impedance and facies inversion of time-lapse seismic data for improving monitoring of CO2 incidentally stored from CO2 EOR / Barajas-Olalde C., Mur A., Adams D. C. [et. al] // 15th International Conference on Greenhouse Gas Control Technologies (GHGT-15), 2021, March.
6. Gunning J. Joint impedance and facies inversion — seismic inversion redefined / Gunning J., Kemper M. // First Break. — 2014. — P. 89–95.
7. Zhao H. Nesterov-accelerated adaptive momentum estimation-based wavefront distortion correction algorithm / Zhao H., An J., Yu M. [et. al] // Optics Express. — 2021. — No. 8. — Pp. 7177–7185.
8. https://sede-open.github.io/synthoseis/datagenerator.html
