ВведениеСовременный этап интерпретации данных сейсмической разведки требует применения комплексных алгоритмов, обеспечивающих согласованные и геологически достоверные результаты. Ключевой вызов возникает на этапе прогноза геологических свойств: сегодня недостаточно только восстановления упругих параметров (акустического импеданса, плотности и др.), необходим прямой прогноз литофаций для построения детальных моделей. Несмотря на наличие ряда высокоэффективных алгоритмов синхронной инверсии (СИ) [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):Ai=12(1+tan2θi),Bi=-4(Vs2Vp2)sin2θi,Ci=12(1-4(Vs2Vp2)sin2θi),A_i=\frac{1}{2}(1+\tan^2\theta_i), B_i=-4(\frac{V_s^2}{V_p^2})\sin^2\theta_i, C_i=\frac{1}{2}(1-4(\frac{V_s^2}{V_p^2})\sin^2\theta_i)x=log(Vp),y=log(Vs),z=log(p)x=\log(V_p), y=\log(V_s), z=\log(p).       (3)Тогда матрица оператора инверсии представляет собой матрицу (4):.       (4)В итоге получаем СЛАУ Sobs = Gm, которая решается методом наименьших квадратов с добавлением регуляризующих компонентов. Обновление модели проводится через поправки с помощью СЛАУ (5):(GTG+R)Δm=GT(Sobs-S)(G^TG+R)Δm=G^T(S_{obs}-S).       (5)1.3. Градиентный метод оптимизации с линеаризацией целевой функцииОсновные предположения для реализации метода взяты из источников [2] и [4].Каждое значение коэффициента отражения вычисляется исходя из значения физических свойств в текущей точке и в следующей:ri = f(Vpi, Vpi+1, Vsi, Vsi+1, ρi, ρi+1).Минимизируется функционал вида (6):||ds-d||2=||w*r-d||2→min||d_s-d||^2=||w*r-d||^2→min.       (6)Для трассы коэффициентов отражения выполним разложение в ряд Тейлора до первого члена: r(m) = r(m0) + FΔm, где r(m) — трасса коэффициентов отражения, r(m0) — начальная трасса коэффициентов отражения, построенная по низкочастотной модели, F — якобиан, Δm — вектор поправок физических свойств.В качестве уравнения для коэффициента отражения используется уравнение Нотта — Цеппритца.Подставив разложение в (6), получим выражение (7):||w*FΔm-Δd0||2→min,Δd0=d-w*r(m0)||w*FΔ \mathbb{m}- Δd_0 ||^2→min,Δ\mathbb{d_0}=d-w*r(m_0 ).       (7)Приращение m вычисляется методом Nesterov ADAM [7]: ∆m=-α(g+r)∆m= -α(g+r),       (8)g=(WF)T(Wr(m0)-d)g=(WF)^T (Wr(m_0 )-d),где αα — шаг градиентного спуска, 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.860.960.39Аппроксимация Аки — Ричардса и решение СЛАУ (1.2)0.840.977.09Градиентный метод оптимизации с линеаризацией целевой функции через уравнение Нотта — Цеппритса (1.3)0.970.9913.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​, m​Rho] — вектор текущих значений упругих свойств в точке (акустический импеданс, сдвиговой импеданс, плотность), полученных на текущей итерации инверсии;μ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. Метрики расчета синтетического разрезаМетрикаЗначениеCorrAi0.8811CorrSI0.7147CorrRho0.6126F1-scoreнеколлектор0.9361F1-scoreколлектор0.6599Относит. СКО между синтетикой и сейсмикой по итерациям литологической инверсии0 : 1.00711 : 0.06042 : 0.02903 : 0.017584 : 0.01135 : 0.00786 : 0.00567 : 0.00428 : 0.00349 : 0.002710 : 0.0023 Из таблицы 2 можно сделать вывод, что алгоритм сходится к сейсмическим данным. Чем больше число итераций литологической инверсии, тем точнее литологическая картина. Коэффициенты корреляции в таблице 2 значительно ниже, чем в пункте 1.3, что объясняется упрощенной регуляризацией синхронной инверсии рассматриваемого алгоритма. Коэффициенты регуляризации фиксированы на протяжении всех итераций и могут быть неоптимальными, что открывает направление для дальнейших исследований по их адаптивному подбору.ВыводПо результатам сравнительного анализа трех реализаций синхронной инверсии выбран градиентный метод оптимизации с линеаризацией целевой функции. Он обеспечивает наилучший баланс точности и скорости. Его использование в методе литологической инверсии полностью оправдано.Разработанный на этой основе алгоритм продемонстрировал хорошее соответствие синтетической модели. Дальнейшие исследования будут направлены на апробацию алгоритмов на реальных геолого-геофизических данных.