Аникиев Д.В. Методы обращения сейсмических волновых полей. Удк методы обращения сейсмических волновых полей
 Скачать 1.61 Mb.
|
|
f, (где L – волновой оператор, содержащий информацию о параметрах среды, u – вектор смещения, а f соответствует функции источника. Пусть также модель среды m есть сумма гладкой опорной модели m 0 (трендовой составляющей) и модели резких локальных возмущений δm: m = m 0 + Линейный оператор L представим в виде суммы операторов L = L 0 + δL, где L 0 связан только сопор- ной моделью m 0 , а δL связан только с моделью возмущений. Аналогично, поле u представляется в виде суммы полей u = u 0 + δu, где поле u 0 удовлетворяет уравнению L 0 u 0 = f, а рассеянное волновое поле δu удовлетворяет уравнению δLδu = f. Таким образом, уравнение запишется в виде L 0 u 0 + δLu 0 + L 0 δu + δL δu = f. В приближении Борна пренебрегаем слагаемым го порядка малости, что в итоге дает систему уравнений 0 0 0 0 u f u u = = − ⎧ ⎨ ⎩ Тогда для рассеянного поля δu будет справедливо уравнение ω ω ξ δ ξ ω ξ ξ u x x G u u s X , ; ( , ; ) ( ) ( , ) ( ) = − ⊂ ∑ 0 0 L , (где суммирование выполняется по пространственной области X, в которой задана модель возмущений δm, а G 0 – тензор Грина для опорной модели среды Векторное уравнение можно записать в виде системы линейных уравнений = Aδm. (Таким образом, изначально нелинейная зависимость сейсмических данных от параметров среды ли- неаризуется. Перепишем уравнение в виде = A T δu, (что более удобно, так как матрица A T A квадратная и можно формально записать решение для δm в виде = (Если рассмотреть задачу сточки зрения теории оптимизации, согласно которой минимизируется среднеквадратический функционал невязки вида J = || δu –Aδm|| 2 , то матрица A T A в уравнении будет представлять собой матрицу Гессе (Гессиан, составленную из вторых производных функционала J по модельным параметрам (с учетом линеаризации. Выражение по-разному называется и интерпретируется в зависимости от теории, в рамках которой изучается обратная задача. В дифракционной томографии вводится понятие томографического функционала Рыжиков, Троян, 1994], который имеет смысл чувствительности отдельных измерений к изменению параметров в различных пространственных областях среды. [Dahlen et al., 2000; Liu, Tromp, 2008; Tromp et al., 2005] для описания подобного функционала используют термин ядро чувствительности (sensitivity kernel). Наконец, в работах, использующих для построения модели среды оптимизационные методы, фигурирует понятие градиента функционала невязки (gradient of a misfit functional) [Plessix, 2006; Virieux, Op- erto, В ЛСТ томографический функционал является сингулярными локализуется вдоль луча, соединяющего источники приемник. В дифракционной томографии носитель функционала в случае однородной опорной среды и точечного источника локализуется в параболическом слое, когда падающее поле является плоской волной, ив эллиптическом слое, когда падающее поле является сферической волной, если источник является точечным. При этом толщина слоя зависит от частоты зондирующего импульса. Поскольку матрица Гессе A T A близка к вырожденной, то вместо явного обращения используется псевдообращение Мура–Пенроуза [Golub, Kahan, 1965] либо регуляризация по Тихонову Тихонов, Ар- сенин, 1974] или в статистической постановке по Трояну [Троян, Киселев, 2000]. Однако размер матрицы Гессе для реалистичных моделей среды настолько велик, что непосредственное ее вычисление, хранение или обращение является чрезвычайно ресурсоемкой задачей. В качестве альтернативы непосредственному вычислению элементов матрицы Гессе можно использовать итеративный метод сопряженных градиентов [Hestenes, Stiefel, 1952], который позволяет избежать явного вычисления матрицы Гессе. На данный момент в научном сообществе активно обсуждается проблема совместного восстановления упругих параметров среды (см [Аникиев и др, 2013; Anikiev et al., 2013, 2014; Gholami et al., 2013; Jeong et al., 2012; Prieux et al., 2013]). [Tura, Johnson, 1993] исследовали стабильность линеаризованного обращения для случая упругой модели среды и предложили подход, согласно которому различные параметры среды могут быть восстановлены из соответствующих типов волн рассеянного поля. В работе [De Nicolao et al., 1993] рассмотрена линеаризованная обратная задача для упругости сточки зрения сингулярного разложения (SVD) оператора, связывающего модельные параметры и входные данные. Их анализ показал, что на основе данных отраженных продольных волн можно дать оценку только для одного из трех параметров упругой среды. Позже [Lebrun et al., 2001] также проанализировали сингулярное разложение и показали, что параметры упругой среды немо- гут быть восстановлены с равной степенью достоверности, причем наибольшие затруднения вызывает восстановление плотности. Линеаризация предполагает, что опорная модель среды близка к истинному распределению параметров. Важным аспектом при этом является скорость сходимости итеративной схемы. Число итераций зависит от обусловленности исходной обратной задачи. Чтобы ускорить сходимость метода сопряженных градиентов, можно использовать предобусловливание [Капорин, 2011; Axelsson, 1996], которое заключается в домножении обеих частей уравнения слева на симметричную положительно определенную матрицу- предобусловливатель P: PA T A δm = Матрица P должна быть близка к обращенной матрице Гессе, но если в качестве предобусловливате- ля взять саму матрицу (A T A) –1 , то метод сопряженных градиентов сойдется за одну итерацию, то есть в этом случае попросту отсутствует необходимость в итерационном методе. Таким образом, скорость сходимости не единственный важный аспект предобусловли- вания. На совокупное качество предобусловливания влияет ряд факторов затраты на предварительное вычисление предобусловливателя и на его применение на каждой итерации сокращение числа итераций, необходимых для получения достаточного качества решения [Капорин, В последнее время появилось множество подходов к проблеме предобусловливания для задач сейсморазведки (см [Gelis et al., 2007; Hak, Mulder, 2008; Hu et al., 2011; Kormann et al., 2013; Korta et al., 2013; Liu et al., 2013; Plessix, Mulder, 2004; Tang, 2009; Tang, Lee, 2010]). Среди них наиболее часто используемым является метод Якоби, заключающийся в построении предобусловливателя в виде обращенной диагональной части матрицы Гессе. При подобной аппроксимации учитывается лишь взаимнокорреляционная функция полей, рассеянных от одной пространственной точки, ноне берется в расчет взаимное влияние полей, рассеянных от соседних пространственных точек. Тем не менее подобное предобусловливание позволяет существенно улучшить изображение среды для целевых геологических объектов. Например, на рис. 1 представлено сравнение результата стандартной миграции, результата итеративного линеаризованного обращения с предобусловли- ванием по методу Якоби и истинного распределения отражающей способности для локальной области модели. Рисунки взяты из работы [Tang, В случае, когда метод Якоби неэффективен, например, при обращении волновых полей в случае упругой модели среды [Аникиев и др, 2013; Anikiev et al., 2013], могут быть использованы более сложные подходы для построения предобусловливателя. В частности, метод неполной факторизации, метод неполного разложения Холецкого и метод приближенного обратного треугольного разложения. Эффективный вариант предобусловливателя для линеаризованного обращения в случае упругой изотропной среды, основанный на неполном обратном треугольном разложении разреженной аппроксимации матрицы Гессе [Капорин, 2011, 2012], представлен в [Anikiev et al., 2014]. Стоит отметить, что более сложные методы предобу словливания опираются на априорные предположения о структуре разреженности матрицы Гессе. Развитие методов предобусловливания является важным научным направлением в задаче совместной оценки параметров упругой среды в областях с неравномерной освещенностью СВЯЗЬ ЛИНЕАРИЗОВАННОГО ОБРАЩЕНИЯ С МИГРАЦИОННЫМИ МЕТОДАМИ Резкие возмущения параметров, восстанавливаемые при линеаризованном подходе, формируют локальные отражающие и рассеивающие объекты, обеспечивающие резкое изменение направления распространения сейсмической энергии, вплоть до ее возвращения на свободную поверхность. На получение изображений таких геологических объектов также направлены все существующие методы сейсмических миграционных преобразований, основы которых заложены в работах [Тимошин, 1962, 1972, 1978; Claer- bout, 1971; Hagedoorn, 1954; Mursgrave, 1961; Rockwell, При удачном выборе опорной модели среды миграционное преобразование корректно определяет расположение и геометрию соответствующих отражающих границ [Claerbout, 1985]. Тем не менее структурное изображение, получаемое с помощью миграции, не позволяет судить об изменчивости параметров среды вдоль этих границ [Zhu et al., Как показано [Lailly, 1983; Tarantola, 1984], первый шаг итеративного поиска решения обратной задачи является эквивалентным процедуре миграции до суммирования [Тимошин, 1972; Claerbout, 1971; Schneider, 1978]. Связь между решением обратной задачи в линеаризованной постановке и миграционными преобразованиями также исследована [Beylkin, 1985; Bleistein, 1987]. Было показано, что на основе лучевого представления для функции Грина можно построить весовые коэффициенты, использование которых позволяет получить изображение в истинных амплитудах. Под термином истинная амплитуда, как правило, понимается величина коэффициента отражения. Сходство и одновременное независимое развитие миграционных методов и методов линеаризованного обращения волновых полей привело к появлению термина “миграция/инверсия” (в англоязычной литературе – migration/inversion) [Beydoun, Mendes, 1989; Lambare et al., 2003; ∅stmo et al., 2002; Tang, 2009; Xu et al., 1999]. Достаточно полный и подробный анализ современных проблем, связанных с получением волновых изображений среды в истинных амплитудах, выполнен [Hertweck, Рассмотрим эквивалентность миграции и линеаризованного обращения подробнее на примере миграции с обращением времени (reverse-time migration) и линеаризованной обратной задачи для случая упругой изотропной среды. Систему уравнений движения в спектральной области запишем следующим образом [Slaughter, 2002]: ω ρ ω λ ω μ ω ω 2 ( ) ( , ) ( , ) ( x u x x u x I x u x u + + ⎡ ⎣⎢ + + div ( )div ( , ) ( ) grad T grad ,, ) , , x x ( ) ⎤ ⎦ где u( ω, x) – волновое поле, f(ω, x s ) – функция источника единичный тензор, ω – частота, ρ(x) – плотность) и μ(x) – параметры Ламе. Тогда томографический функционал (или ядро чувствительности, характеризующий плотность ρ, можно записать как g ρ ω ω ω ω δ ω ( ) Re , ; , ; , ; , , x G x x u x x u x x x x = ( ) ( ) ( ) ∑ 2 T g s g s s g 0 где G 0 ( ω, x; x g ) – тензор Грина для опорной среды, u 0 ( ω, x; x s ) – поле от источника в опорной среде, а δu(ω, x g ; x s ) есть разница наблюденного поля u obs ( ω, x g ; x s ) и смоделированного поля u 0 ( ω, x g ; x s ) в точке приемника x g . Черта над δu в формуле обозначает комплексное сопряжение, которое эквивалентно обращению знака времени во временной области. Таким образом, для построения градиента в точке среды x необходимо просуммировать произведения спектров Рис. 1. равнение результата стандартной миграции (а, результата итеративного линеаризованного обращения с предобусловливанием по методу Якоби (20 итераций) (б, истинного распределения отражающей способности (в) для локальной области (красный прямоугольник) модели Sigsbee2A (г двух волновых полей (прямого и обратного) для каждой частоты ω, источника x s и приемника x g . Прямое поле u 0 ( ω, x; x s ) распространяется в опорной модели от источника с функцией f( ω, x s ), расположенного в точке x s . Обращенное во времени поле G x x u x x g x g s r 0 T ω δ ω , ; , ; ( ) ( ) ∑ распространяется также в опорной модели, но сразу от нескольких приемников, в каждом из которых в качестве сигнала используется разность зарегистрированного и смоделированного значений поля, записанная в обратном времени. Томографические функционалы, связанные с чувствительностью обращения по параметрам Ламе, выглядят иначе (см, напр Рыжиков, Троян, 1994; Fichtner, 2010; Liu, Tromp, 2008]), но по-прежнему в их основе лежит произведение прямого и обращенного во времени волновых полей. Как можно заметить, формула вида (5) напоминает изображающий принцип в миграции с обращением времени [Тимошин, 1972; Baysal et al., 1983; Claerbout, 1971]: I ( ) , ; , ; , , x u x x u x x f s r g x где u f – прямое поле, распространяющееся от источников в опорной среде, а u r – обратное поле, получающееся продолжением в обратном времени колебаний, зарегистрированных приемниками. При этом полагается, что данные для каждого приемника записаны для достаточного интервала времени, а обращение во времени получается путем перестановки временных отсчетов сейсмической трассы в обратном порядке. Существенным отличием миграционного преобразования от томографического функционала является то, что в миграции строится продолжение в обратном времени самого зарегистрированного волнового поля, тогда как в алгоритме обращения строится продолжение разницы зарегистрированного и смоделированного волновых полей. Как видно из представленного анализа, ив миграции, ив обращении сейсмических данных ключевой является задача продолжения волновых полей в опорную среду. Этой проблеме посвящена монография [Петрашень, Нахамкин, 1973], где авторы предложили оригинальный подход к анализу и обработке данных сейсморазведки. Выполненный авторами анализ значительно опередил аналогичные зарубежные исследования в данном направлении [Berkhout, 1981; Berryhill, 1979, 1984]. Таким образом, работа Петра- шень, Нахамкин, 1973] предвосхитила появившиеся впоследствии практические алгоритмы продолжения волновых полей. В настоящем номере журнала Технологии сейсморазведки публикуется ряд статей, отражающих основные теоретические и практические аспекты, представленные в указанной монографии. МЕТОД ПОЛНОГО ОБРАЩЕНИЯ ВОЛНОВЫХ ПОЛЕЙ Развитие в середине – второй половине ХХ в. теории оптимального управления Алексеев и др, 1979; Понтрягин, 1959] оказало значительное влияние на методы решения обратной динамической задачи сей- смики. Теория оптимального управления подразумевает сведение любой задачи к минимизации или максими- зации некоторого функционала, который и определяет условие оптимальности решения задачи. В обратной динамической задаче сейсмики основным критерием определения набора физических параметров m = { m 1 , m 2 , …, m N } является степень соответствия наблюденных и моделированных данных. Наиболее естественным функционалом, определяющим соответствие зарегистрированных и смоделированных в опорной среде сейсмических данных как с математической, таки с вычислительной точки зрения, является норма в функциональном пространстве, поскольку сейсмограмма есть нечто иное, как функция координат источника x s , приемника x g и времени t, определенная для дискретного набора значений переменных. Применение теории управления к одномерной обратной задаче для волнового уравнения подробно изложено в работах [Bamberger et al., 1977, 1979], авторы которых предложили использовать в качестве функционала невязки полей распространенную в математической физике L 2 норму u m u m x x - u x x s g s g x x s r ( ) , , , , , , , ( ) = ( ) ( ) ∑ 1 2 2 t t t obs , (где u есть моделированное, а u obs – зарегистрированное волновые поля. Начиная поиск оптимального набора параметров с некоторой стартовой модели, [Bamberger et al., 1977, 1979] предложили использовать серию линейных спусков для продвижения к глобальному минимуму функционала. В настоящее время предложенный [Bamberger et al., 1977, 1979] метод в англоязычной литературе известен как полное обращение волновых полей (full waveform inversion; FWI). Как происходит с любым развивающимся направлением, с названием метода связана некоторая неопределенность. Термин обращение полных волновых полей (full-waveform inver- sion), вероятно, впервые появился в работе [Pan et al., 1988], где методика, введенная [Bamberger et al., 1977, 1979], применялась для восстановления структуры среды на основе исходных сейсмограмм. Однако довольно часто при обращении сейсмических волновых полей таили иная часть информации, содержащейся в данных, не используется. Например, данные обращаются только в определенном частотном диапазоне или же игнорируются поверхностные волны. Тем не менее термин FWI используется ив этих случаях, поскольку по-прежнему используется подход к обратной задаче сейсмики как к нелинейной проблеме оптимизации. Таким образом, в настоящее время слово полное относится скорее не к полноте используемых сейсмических данных, а к самому методу, позволяющему учесть всю информацию, содержащуюся в волновых полях. Существенной особенностью метода при этом являются составление функционала невязки на основе полных волновых полей (пусть и после процедур фильтрации, а не их отдельных кинематических или динамических атрибутов, и последовательное приближение к модели, наиболее правдоподобно описывающей зарегистрированные волновые поля. Исходя из этих соображений, наиболее логично назвать FWI методом полного обращения волновых полей. В работах [Bamberger et al., 1977, 1979] предложены способ итеративной минимизации функционала невязки и процедура вычисления градиента функционала невязки на каждом шаге минимизации ВЫЧИСЛЕНИЕ ГРАДИЕНТА |