Аникиев Д.В. Методы обращения сейсмических волновых полей. Удк методы обращения сейсмических волновых полей
 Скачать 1.61 Mb.
|
|
Технологии сейсморазведки, № 1, 2014, с. 38–58 http://ts.sbras.ru УДК МЕТОДЫ ОБРАЩЕНИЯ СЕЙСМИЧЕСКИХ ВОЛНОВЫХ ПОЛЕЙ Д.В. Аникиев, В.В. Казей, Б.М. Каштан, А.В. Пономаренко, В.Н. Троян, Р.А. Шигапов Санкт-Петербургский государственный университет, Санкт-Петербург, Петергоф, ул. Ульяновская, 1, Россия, e-mail: d.anikiev@spbu.ru, v.kazei@spbu.ru, b.kashtan@spbu.ru, a.v.ponomarenko@spbu.ru, v.troyan@spbu.ru, В статье описаны методы получения оценок пространственного распределения упругих параметров геологической среды поданным сейсморазведки. Обзор включает широкий спектр технологий решения обратной задачи сейсмики – от лучевой сейсмической томографии и обращения дисперсионных кривых, позволяющих восстановить трендовую составляющую модели среды, до метода полного обращения волновых полей и миграции в обратном времени, позволяющих значительно повысить разрешение изображения. Метод полного обращения волновых полей является на данный момент одним из наиболее актуальных и интенсивно развивающихся направлений сейсморазведки. Рассматриваются различные современные модификации метода и обсуждаются теоретические и практические проблемы. Метод полного обращения волновых полей, миграция с обращением времени, линеаризованное обращение, продолжение волновых полей, обратная кинематическая и динамическая задачи сейсмики, оптимизация, поверхностные волны, затухание, анизотропия, методы изображения в высоком разрешении OF SEISMIC WAVEFORM INVERSION D.V. Anikiev, V.V. Kazei, B.M. Kashtan, A.V. Ponomarenko, V.N. Troyan, R.A. Shigapov Saint Petersburg State University, 198504, 1, Ulyanovskaya str., Peterhof, Saint Petersburg, Russia, e-mail: d.anikiev@spbu.ru, v.kazei@spbu.ru, b.kashtan@spbu.ru, a.v.ponomarenko@spbu.ru, v.troyan@spbu.ru, r.shigapov@spbu.ru The paper concerns methods for elastic medium parameter estimation using seismic data. The overview covers a wide range of seismic inversion technologies, from ray-based seismic tomography and inversion of dispersion curves, which allows recovering general trends of the medium model, to full waveform inversion and reverse-time migration providing high-resolution imaging. The full waveform inversion method is cur- rently one of the most important and rapidly developing lines in seismic exploration. The discussion also ad- dresses various up-to-date modifications of the method and relevant theoretical and practical issues. Full-waveform inversion, reverse-time migration, linearized inversion, wavefield continuation, kinematic and dynamic inversion, optimization, surface waves, attenuation, anisotropy, high-resolution imaging ВВЕДЕНИЕ Основной проблемой сейсморазведки является решение обратной задачи, те. определение внутреннего строения Земли на основании сейсмических волновых полей, регистрация которых осуществляется на дневной поверхности, в приповерхностном слое или в скважинах. Первые методы обращения наблюдаемых сейсмических волновых полей относятся к началу ХХ в, когда технические характеристики сейсмоприемников ограничивали возможность изучения динамических атрибутов сейсмических волновых полей, связанных с формой волнового пакета, частотой и энергией. В связи с этим первые практические методы решения обратной задачи сейсмики использовали лишь кинематические атрибуты, связанные с формой волновых фронтов, траекториями лучей, с распределением времен распространения сейсмических волн. Стечением времени регистрирующая аппаратура совершенствовалась, стали появляться хорошо откалиброванные сейсмограммы. Эти обстоятельства сделали возможным применение методов обращения волновых полей, использующих не только кинематические, но иди нами ческие сейсмические атрибуты. Во многом благодаря совершенствованию сейсмометрического оборудования и развитию ЭВМ большую популярность обрело научное направление, связанное с решением обратной задачи сейсмики. Значительный прогресс в развитии сейсморазведки произошел с переходом к цифровой регистрации сейсмических данных и их компьютерной обработке. Если ранее было необходимо непосредственное участие сейсморазведчика во всех этапах обработки и интерпретации, то сейчас предпринимаются шаги к автоматизации процесса обработки данных вплоть до геологической интерпретации. Представленный в настоящей статье обзор позволяет проследить все этапы развития современных способов решения обратной задачи сейсмики. Одним из наиболее прогрессивных методов как сточки зрения удобства автоматизации обработки данных, таки сточки зрения точности оценки физических параметров исследуемой среды является метод полного обращения сейсмических волновых полей. Научное направление, связанное с решением динамической об Д.В. Аникиев, В.В. Казей, Б.М. Каштан, А.В. Пономаренко, В.Н. Троян, Р.А. Шигапов, 2014 ратной задачи сейсмики, интенсивно развивается, и хотя подробный обзор метода полного обращения волновых полей предложен в работах [Чеверда, 2009; Virieux, Operto, 2009], появилось большое количество исследований, способствовавших дальнейшему развитию этой методики. В статье рассматриваются различные современные модификации метода полного обращения волновых полей, обсуждаются существующие проблемы и последние достижения в теоретической и практической областях. КИНЕМАТИЧЕСКИЕ МЕТОДЫ Впервые кинематическая обратная задача для дифференциального уравнения эйконала | ∇τ| = V –1 (x) была сформулирована в 1907 г. в работах [Herglotz, 1907; Wiechert, 1907], где для сферически симметричной модели Земли, те. в предположении, что скорость монотонно зависит только от глубины, были доказаны существование и однозначность определения искомой скорости V по известной функции τ – времени прихода сейсмического сигнала из одной точки поверхности в другую. Позже подход был адаптирован для целей сейсморазведки в работе Чибисов, 1934]. В дальнейшем теория одномерной кинематической задачи как для рефрагированных, таки для отраженных волн получила развитие в работах [Бессонова и др, 1973; Гейко, 1970; Гервер, Маркушевич, 1967; Гольдин, 1977; Молчан, 1977; Павленкова, 1973; Риз- ниченко, Первые исследования обратной кинематической задачи для горизонтально-неоднородной модели среды выполнены [Лаврентьев, Романов, 1966]. На их основе был создан численный алгоритм решения трехмерной обратной кинематической задачи Алексеев и др, 1969]. Многомерная обратная кинематическая задача рассматривается также в работах Ани- конов, 1975, 1978; Бейлькин, 1979; Бернштейн, Гер- вер, 1978]. В работах [Глоговский, 2011; Glogovsky et al., 2009] рассмотрены вопросы корректности и структурной устойчивости обратной кинематической задачи сейсмики. ЛУЧЕВАЯ СЕЙСМИЧЕСКАЯ ТОМОГРАФИЯ Развитие в 1970-x гг. вычислительной томографии и ее применение в различных областях науки привело к появлению аналога такого метода в геофизике – лучевой сейсмической томографии (ЛСТ). Основы лучевой томографии описаны в работах [Андерсон, Дзевонский, 1984; Aki et al., 1977; Aki, Lee, 1976; Dziewonski et al., 1977; Dziewonski, Anderson, 1981]. Подход получил широкое применение в сейсмологии для изучения крупных континентальных областей. В области прикладной сейсмики известны работы Рыжиков, Троян, 1985; Burmakov et al., 1984; Dorbath et al., 1986; Ellsworth, Koyagani, 1977; Grasso et al., 1983; Maguire et al., 1985; Mitchell et al., 1977; Nerces- sian et al., 1984; Reasenberg et al., 1980]. Подробный обзор литературы по лучевой сейсмической томографии за е гг. представлен [Nolet, 1987]. Все редине х гг. СВ. Гольдин опубликовал серию работ [1996а,б,в, 1997, 1999а,б], в которых была детально исследована проблема единственности и устойчивости решения томографической задачи. В ЛСТ при межскважинных наблюдениях в качестве исходных данных для восстановления скорости обычно используются времена прихода первых вступлений прямых или преломленных волна в случае наблюдений по методу отраженных волн – времена прихода выделенных отраженных волн [Lo, Inderwiesen, 1994]. Типичная модель лучевой томографии записывается в виде где в сейсмическом случае τ i – момент прихода сейсмической волны, распространяющейся вдоль луча L i , f(x) = V –1 (x) – обратная скорость или функция медленности случайная компонента модели. Решение для однородной опорной модели среды можно найти при помощи обратного преобразования Радона [Radon, 1917], однако этот способ не получил широкого распространения в сейсморазведке. На практике чаще всего пользуются алгебраическим методом реконструкции, который сводится к следующему неизвестную функцию f(x) представляют в виде разложения по базисным функциям Ψ j (x): f j j j M ( ) ( ) x x = = ∑ Θ где Θ j – набор неизвестных параметров. Данные наблюдений записываются в этом случае следующим образом: τ ε i ij j i j M C = + = ∑ Θ 1 , где C dl ij j L i = ∫ Ψ ( Выбор базисных функций зависит от априорных представлений о свойствах среды (симметрия, слоистость и т. д. На практике чаще всего используются базисные функции ) , , x x x = ∈ ∉ ⎧ ⎨ ⎪ ⎩⎪ 1 В этом случае все исследуемое пространство (в случае 3D) или плоскость (в случае 2D) разбивается на элементарные элементы объема или плоскости Элементы матрицы ∅ ∩ = либо равны константе C, если луч L i пересекает элемент объема или плоскости Ω j , либо равны нулю, если не пересекает ( ∩ означает пересечение, ∅ – нулевое множество). Окончательно модель измерений запишется в виде = C Θ + Предположим, что случайная компонента ε и вектор параметров Θ подчиняются нормальным законам распределения ( ε∈N(0, R ε ) с нулевым математическим ожиданием и матрицей ковариации R ε , ΘN(<Θ>,R Θ ) с математическим ожиданием < Θ> и матрицей ковариации R Θ . Задача состоит в получении оценки для вектора параметров Θ. Для этого применяется метод максимума апостериорной вероятности [Гольцман, 1971; Троян, 1982; Tarantola, 1987, 2005]: Θ Θ Θ Θ ττ Θ Θ Θ Θ = ( ) ( ) ⎡⎣ ⎤⎦ max где p( Θ) – априорная плотность распределения, а p( τ/Θ) – условная плотность распределения. В случае нормального распределения векторов ε и такой критерий сводится к == ττ Θ Θ ττ Θ Θ Θ Θ Θ Θ Θ Θ Θ Θ Θ Θ Θ Θ −−11 Явная запись искомой оценки и ее ковариационная матрица могут быть представлены в двух видах = + ( ) + ( ) − − − C R C R C R R T T ε ε 1 1 1 , R Θ Θ Θ Θ −−11 ( ) = + ( ) − − C R C R T ε 1 В частном случае некоррелированных векторов R I ε ε σ = 2 и σ 2 оценка может быть представлена в виде = + ( ) + ( ) − C а соответствующая ей ковариационная матрица – R Θ Θ Θ Θ 22 ( ) = + ( ) − σ α C где α σ σ ε = 2 2 Θ – регуляризующий параметр Тихонов, Арсенин, Ковариационная матрица R Θ Θ ( ) может быть использована для оценки качества интерпретации. Дисперсии ошибок восстановления расположены на главной диагонали матрицы ) , при этом внедиаго- нальные элементы описывают корреляционные связи между компонентами оцениваемых параметров. Лучевая томография подразумевает использование лишь кинематических атрибутов волнового поля, что заведомо накладывает на искомое решение обратной задачи ограничения, связанные как с предельной разрешающей способностью, обусловленной локальными размерами первой зоны Френеля (см, напр [Yilmaz, 1987]), таки с проблемой многозначности времен пробега сейсмических волн, обусловленной распространением сейсмической энергии по различным путям в сложнопостроенных средах. Тем не менее результаты лучевой томографии по восстановлению функции скорости в среде часто используются для построения опорной модели среды, необходимой для последующего, более детального ее изображения. Одним из необходимых шагов в процессе автоматизации обработки данных является отказ от визуального определения первых вступлений (годографов) сейсмических волн, обыкновенно предваряющего томографические методы скоростного анализа. Автоматизация пикирования позволяет существенно увеличить скорость обработки данных [Billette et al., Полнота использования информации, содержащейся в сейсмических измерениях, однозначно связана с областью применимости прямой задачи. Лучевая сейсмическая томография предполагает в качестве информативных параметров времена прихода волн определенного типа. Подобное сворачивание информации в один параметр существенно обедняет содержание эксперимента, поэтому важным шагом к более полному и точному описанию внутреннего строения Земли и к дальнейшему развитию сейсморазведки является использование динамических атрибутов регистрируемых сейсмических волн. ЯВНЫЕ ДИНАМИЧЕСКИЕ МЕТОДЫ * Впервые динамическая обратная задача сейсмики для горизонтально-однородной среды (так называемая обратная задача Лэмба) детально исследована в работе Алексеев, 1962], где автор строит решение в спектральной области на базе математического аппарата теории обратного рассеяния, развитого в работах [Гельфанд, Левитан, 1951; Крейн, 1954; Марченко, В статье Благовещенский, 1966] впервые рассмотрено решение динамической обратной задачи во временной области и предложен метод нелинейных вольтерровских систем интегральных уравнений. Независимо предложили альтернативное интегральное уравнение (также во временной области) в задаче восстановления формы речевого тракта человека по акустическим измерениям. Нелинейные уравнения во временной области также использованы в [Symes, 1979]. В статье [Burridge, 1980] исследовано применение уравнений Гельфан- да–Левитана–Крейна–Марченко для теории упругости во временной области и установлена связь между ними и уравнением, представленным [Gopinath, Son- dhi, В спектральной области, помимо работ Алексеев, уравнение [Марченко, 1955] используется в [Ware, Aki, 1969]. В статье [Berryman, Greene, 1980] представлен численный метод решения обратной задачи на основе работы [Ware, Aki, 1969] и показано, что он в точности совпадает с матричным подходом [Goupillaud, 1961] для слоистой среды. Альтернативные подходы в спектральной области изложены, Достаточно обширный обзор зарубежной литературы гг. по явным методам решения динамических обратных задач представлен [Newton, 1981]. Однако ряд важных результатов, в основном связанных с решением динамической обратной задачи для гиперболических уравнений го порядка, принадлежит отечественным математикам [Аниконов, Ерохин, 1981; Благовещенский, 1971; Белишев, 1978, а, б б Бухгейм, Кардаков, 1978; Запреев, Чеверда; Кабанихин, 1984, 1988; Романов, 1972, 1983; Яхно, Численная реализация спектрального метода, изложенного в статье [Carroll, Santosa, 1981], была предложена, который использовал дискретный аналог уравнения Гельфанда–Левитана для восстановления акустического импеданса и оценил погрешность решения обратной задачи при наличии * Методы, основанные на явном обращении нелинейного оператора обратной задачи и решении возникающих при этом интегральных уравнений, получили название явные шума. Стабильность такого решения рассмотрена [Carroll, Santosa, В работе [Carrion, 1985] показано, что трехмерная обратная задача сейсмики для горизонтально-одно- родной среды может быть сведена к одномерной задаче путем преобразования Радона и решена для определенной пары параметров ( τ, p). Здесь же представлена схема численного решения обратной задачи сейсмики на основе метода характеристик во временной области и показано, что использование метода характеристик позволяет восстанавливать структуру среды с кусочно-непрерывным распределением пара- метров. В монографии [Кабанихин, 1988] подробно рассмотрены основные итерационные методы решения систем интегральных уравнений и изучена их сходимость. В статье [Кабанихин, Баканов, 1999] представлен дискретный аналог метода Гельфанда–Левитана в двумерной обратной задаче для скалярного гиперболического уравнения. В работах Романов, 2007, 2008] исследуется проблема устойчивости обратной задачи для уравнения упругости в случае импульсного источ- ника. Несмотря на строгий математический подход, существенным недостатком перечисленных явных методов решения динамической обратной задачи сейсмики является их ориентация на вертикально-неодно- родные среды. Теоретически возможно распространить полученные результаты на среды с горизонтальной неоднородностью. Например, в работе [Blagovestchen- skii et al., 2006] предлагается численная схема решения обратной динамической задачи для среды со слабой горизонтальной неоднородностью, однако практическое применение таких подходов сопряжено со значительными требованиями к входным сейсмическим данным. Одним из перспективных методов решения обратной динамической задачи для произвольной неоднородной среды является метод граничного управления (так называемый Boundary Control method, или метод. Основы ВС-метода для решения многомерной динамической обратной задачи представлены в работе [Белишев, б, далее метод получил развитие в работах [Белишев, 1990, 2008; Belishev, 1997, 2001]. Однако практической реализации метода для решения обратной задачи сейсморазведки на данный момент не осуществлено. ДИНАМИЧЕСКИЕ МЕТОДЫ РЕШЕНИЯ ЛИНЕАРИЗОВАННОЙ ОБРАТНОЙ ЗАДАЧИ Скоростную модель среды условно можно разделить на две составляющие плавно меняющуюся трен- довую составляющую, которая характеризует время распространения волн между удаленными точками среды, и резко меняющуюся составляющую, которая включает все локальные возмущения параметров среды, не оказывающие значительного влияния на времена пробега волн между удаленными точками среды. Зачастую обращение волновых полей рассматривается в рамках теории обратного рассеяния, то есть в предположении, что трендовая составляющая модели известна, а неизвестны лишь локальные неоднородности резкие возмущения физических параметров среды. Именно наличие таких возмущений приводит к возникновению регистрируемых отраженных, дифрагированных и рассеянных волн. В предположении заданной опорной модели можно линеаризовать задачу, воспользовавшись приближением Борна, и тем самым записать систему линейных уравнений относительно модели резких локальных возмущений. Линеаризованная обратная задача для волнового уравнения в случае вертикально-неоднородной модели среды впервые исследована в работе Романов, 1969]. Здесь, а также в [Лаврентьев и др, 1969; Романов показано, что задача восстановления локальных неоднородностей скорости распространения волн сводится к одной из множества задач дифференциальной геометрии. Решение одномерной обратной задачи в линеаризованной постановке, по-видимому, было впервые получено [Cohen, Bleistein, 1977] для однородной акустической вмещающей среды, после чего развито в работах [Bleistein, Cohen, 1979; Cohen, Bleistein, 1979]. Впоследствии метод получил название борновское обращение (Born inversion) и был обобщенна случай вертикально-неоднородной вмещающей среды [Bleistein, Gray, 1985] и произвольной трехмерной вмещающей среды [Bleistein et al., 1987]. Однако впервые неоднородность опорной модели среды была учтена в [Clayton, Stolt, 1981], где за основу взято лучевое представление для функции Грина в акустической среде. Структура возникающего при этом интегрального оператора и особенности численной реализации подхода исследованы [Cohen et al., 1986]. Аспекты, связанные счисленной реализацией линеаризованного обращения на основе лучевого представления для функции Грина, исследованы [Bleistein et al., 1985]. Использование параксиального лучевого метода для приближенного вычисления тензора Грина в случае трехмерной упругой среды обсуждается в статье [Beydoun, Mendes, 1989]. Линеаризованное обращение в совокупности с лучевой теорией, используемой для построения аппроксимации волновых полей, привело к появлению термина лучевое борновское обращение (Ray + Born inversion) [Beydoun, Men des, 1989; Forgues et al., 1994; Lambare et al., В рамках динамической теории упругости линеаризованное обращение данных многократных перекрытий для вертикально-неоднородной опорной среды рассмотрел В.Г. Романов [1983]. Он получили проанализировал систему интегральных уравнений и доказал единственность определения с ее помощью упругих параметров Ламе и плотности. Впоследствии похожий подход к обратной задаче динамической теории упругости применил Т. Ikelle [1995] в спектральной области, позже обобщив его на случай анизотропной среды [Ikelle, 1997]. G. Beylkin [1985] предложил метод решения линеаризованной обратной задачи для акустического уравнения, основанный на обобщенном обратном преобразовании Радона и впоследствии расширенный на случай упругой среды (см [Beylkin, Burridge, 1987, Совершенствование регистрирующей аппаратуры и вычислительной мощности позволило перейти от метода ЛСТ к более полной и потому более сложной методике дифракционной томографии (см Рыжиков, Троян, 1985, а, б, 1989, 1994; Троян, Ки- селев, 2000; Devaney, 1984, 1986; Troyan, Kiselev, 2010]). Метод, как и все ранее перечисленные подходы, основан на теории обратного рассеяния в предположении, что известное распределение скорости распространения сейсмических волн содержит неизвестные малые локальные скоростные аномалии. Термин дифракционная томография применительно к сейсморазведке введен [Devaney, 1984] с целью обозначить принципиальное отличие метода от ЛСТ. В общем случае искомыми параметрами могут являться все три параметра упругой среды – плотность, скорости продольных и поперечных волн или же комбинации этих параметров (см Рыжиков, Троян, 1994; Троян, Киселев, 2000]). Исследование разрешающей способности метода дифракционной томографии проведено [Троян, Крауклис, 1996]. [Devaney, 1984] установлена связь между про- странственно-временным спектром рассеянного поля в приближении Борна и пространственным спектром для малой локальной неоднородности в двумерной однородной опорной среде. Автор показал, что наблюдений по системе ВСП недостаточно для полноценного восстановления локальных неоднородностей. [Wu, ц, 1987] заметили, что если использовать данные ВСП в совокупности с наземными наблюдениями или данными межскважинного просвечивания, то результаты могут быть существенно улучшены. В статье [Mora, 1989] спектральный анализ применен к задаче двух полупространств и показано, что поверхностных наблюдений может быть вполне достаточно для восстановления полного спектра неоднородности, если имеются данные для достаточно больших выно- сов источник–приемник. Данная модель более подробно изучена в [Kazei et al., 2012, 2013a], где указано на необходимость включения в рассматриваемое поле головных волн и волн шепчущей галереи в случае наличия градиента в нижнем полупространстве. Следует отметить, что связь, установленная [Devaney, 1984], существует лишь в случае, когда среда, вмещающая исследуемую локальную аномалию, однородна. Использование локального спектра неоднородности и локального представления поля в виде суперпозиции плоских волн позволяет численно обобщить подход [De vaney, 1984] на неоднородную вмещающую среду [Xie et al., 2006]. Локальное разложение также может быть выполнено аналитически в горизонтально-одно- родной модели с постоянным градиентом [Kazei et al., 2013б]. Поскольку термин томография исторически связан с определением именно скорости распространения волн, то и термин дифракционная томография в настоящее время используется в основном для задач восстановления скоростных аномалий в опорной среде (из недавних работ см, напр [Wu, 2007]). Если же ставится задача нахождения локальных аномалий упругих параметров среды, то используется термин линеаризованное обращение (linearized in- version; см, напр [Hak, Mulder, В работе [Tarantola, 1984] линеаризованная обратная динамическая задача сейсмики трактуется сточки зрения теории оптимизации, что позволяет находить решение путем последовательных итераций, связанных с минимизацией среднеквадратического функционала невязки между наблюденными модельным волновыми полями. Метод получил название итеративное линеаризованное обращение (iterative linearized inversion). Подобные статистические модели интерпретации сейсмических данных предложены ранее [Гольцман, 1971; Троян, 1982]. Итеративный метод асимптотического линеаризованного обращения на базе метода наименьших квадратов обсуждается [Jin et al., 1992]. Метод был применен также к реальным сейсмическим данным [Lambar é et al., 1992]. Два эффективных алгоритма линеаризованного обращения волновых полей (во временной и частотной областях) в случае упругой среды представлены в [Geller, Hara, Рассмотрим основные теоретические аспекты линеаризованного обращения волновых полей. Пусть уравнения движения в среде заданы системой = |