CC BY
37
УДК 550.834 Е.А. Хогоев
ИНГГ СО РАН, Новосибирск
АЛГОРИТМ ПОВЫШЕНИЯ ТОЧНОСТИ РЕШЕНИЯ ЗАДАЧ ДВУМЕРНОЙ СЕЙСМОТОМОГРАФИИ ПРИ ДИАГНОСТИКЕ ВЕРХНЕЙ ЧАСТИ РАЗРЕЗА
E.A. Hogoev
Trofimuk Institute of Petroleum Geology and Geophysics SB RAS Koptyug, 3 , Novosibirsk, 630090, Russian Federation
AN ALGORITHM OF ENHANCEMENT OF ACCURACY OF THE SOLUTION OF 2-D SEISMIC TOMOGRAPHY PROBLEMS AT DIAGNOSTICS OF THE UPPER PART OF A SECTION
The algorithm of solution of 2-D inverse of seismic kinematics problem is considered. For the solution of an inverse problem the Newton's method is used. The solution of an inverse (local) problem consists in the reduction of the integral equation to system of the linear algebraic equations and solution of a task of the least squares for the received system by an iterative method LSQR. The numerical experiment shows, that the method has resolvability.
Определение скоростного строения среды по временам прихода рефрагированных волн (обратная кинематическая задача) является одной из ключевых задач сейсмики. Так в глубинных сейсмических зондированиях (ГСЗ) скорость распространения волн в земной коре является предметом исследования, а в сейсморазведке МОВ - ОГТ от учета скорости в верхней части разреза зависит корректность вводимых статических поправок, которые влияют на прослеживаемость границ и точность структурных построений.
Обратная кинематическая задача сейсмики в своей математической части является классической задачей математической физики. Задача относится к нелинейным, поскольку время прихода волны определяется как интеграл от скорости по кривой (траектории луча), которая сама зависит от распределения скорости. Поэтому решение обратной кинематической задачи сейсмики имеет как практическое значение, связанное с определением скоростного строения среды при сейсморазведке и соответственно повышением точности структурных построений, так и теоретическое, поскольку служит развитию методов решения нелинейных задач, имеющих вид интегрального уравнения первого рода [1].
Данная работа продолжает численные исследования алгоритма последовательной линеаризации задачи с включением на каждом шагу полученного решения (метод Ньютона), опубликованные в [2-4].
Рассмотрим полупространство z > 0 пространства R2(x,z), заполненное средой со скоростью распространения сейсмического импульса V(x, z). Функцию V(x, z) предполагаем непрерывно дифференцируемой. Этому распределению скорости при каждом положении источника £°(х{,0) и
приемника ^(х.,0) отвечает геодезическая Г(о°,о1) = Г. Введем обозначения:
п(х,г) = XIV - функция медленности, и т - время пробега по лучевым траекториям Г в среде с распределением медленности п. Для каждой пары источник - приемник
фо°,о1 )= .
г
У нас есть последовательность соотношений для M пар источник -приемник:
| п с1г - ф = О
г1
| пск-ф1 = 0 . (1)
Г«
| п^-фА = О
Здесь п(х, 2) искомая функция медленности, Гт - зависящая от п(х, 2)
геодезическая, соединяющая пару точек (о0,о^т, ф - известное время пробега сигнала, т = 1,...М.
Переведем (1) в конечномерный вид. Разобьем область, в которой распространяются лучи, прямоугольной сеткой размером (КХ,К2), и для
каждой пары (о0 ,о1)т рассчитаем лучи в среде с задаваемой нами начальной моделью медленности по( х, 2). Лучи рассчитываются методом пристрелки. В
результате мы определяем Г” и невязку времени Аф. В сквозной нумерации
ячеек сетки / = 1,.. N х N 2 длина т-го луча в /-й ячейке а1П1. При соблюдении условий линеаризации система уравнений (1) представляется в виде системы линейных алгебраических уравнений, в которой коэффициенты есть длина отрезков каждого луча в ячейках сетки:
Аф=ф-ф = | пс1г » ^ (пд -п)с1г= | Ап с!г ''у'а^Ащ,
^гп ^ГП ^ГП у
ИЛИ
ААп = Аф . (2)
Для полученной системы (2) проекционно-градиентным методом LSQR (Пейджа - Саундерса) решается задача наименьших квадратов. Критерием останова итераций служит достижение минимума (шт = 10-4)
относительного изменения решения. Решив эту систему уравнений, получаем уточненную модель среды
п = щ+ Ап.
Полученное на очередном шаге решение п заменяет опорную модель п0 для следующей линеаризации. Если на первом шаге линеаризации начальную функцию скорости мы можем выбрать виде формулы, то после решения системы уравнений (2) получаем значения скорости на сетке. Главную сложность итеративного решения представляет собой расчет лучей по сеточно заданной функции медленности. Полученное решение Ап сглаживается посредством медианной фильтрации и последующего усреднения в скользящем окне. При этом не столь важно сохранить в точности найденное решение, сколько обеспечить устойчивость пристрелки лучей. Расчет новых траекторий лучей в уточненной на предыдущем шаге модели среды является ключевым моментом для следующей линеаризации.
Так как расчет лучей производится численным интегрированием системы нелинейных дифференциальных уравнений, нам необходимо знать в каждой точке не только значение скорости, но и ее первых производных по координатам. Для того чтобы обеспечить это требование, дискретно заданная скорость интерполируется двумерным полиномом Лагранжа. Нами использовался следующий шаблон: по 16 точкам значений двумерной функции % = п^, zk), j = 0, ... 3, к = 0, ... 3 строится полином 3-й степени
В большинстве рассмотренных случаях точность вычисления функции и ее производных по координатам оказывается удовлетворительной.
В численном эксперименте [2] было установлено как повышение точности решения обратной задачи, так и увеличение отношения восстанавливаемой части скоростного распределения к его основной известной составляющей. Доказано, что при амплитуде локальной аномалии до 30 % от основной известной составляющей скорости задача решается успешно за три итерации. В дальнейшем удалось показать сходимость итерационного процесса даже при величине восстанавливаемой части (положительного знака) скоростного распределения соизмеримой с его основной известной составляющей [4].
Нами моделировалась система наблюдений с геометрией, близкой к реальной ситуации сейсморазведки МОВ-ОГТ. Расстояние между соседними источниками и приемниками 50 м, расстояние между первой парой источник-приемник Ь0 = 500 м. Всего пар источник - приемник 30, при 9 расстановках - 270 точек наблюдения. Длина базы наблюдения АЬ = 3.4 км, смещение при смене расстановки 5Ь = 0.1 км, общая длина профиля 4.2 км.
В численном эксперименте сетка восстановления охватывает весь профиль по оси X(от -1700 до 2500 м), по оси глубины Z размер оценивается
7 2 2 ~
\/ / 4 “Ь1 / ос 1/ , А Т/ — максимальное удаление, ос — вертикальный градиент в начальной модели, ZM = 800 м. Число ячеек Ы2 = 19 по Z, N = 29
/> у 1^к
точек по X, Nv = Nz Nx = 551, что соответствует длине ячейки сетки по X ~ 150 м, по Z ~ 40 м. Число лучей (размерность данных) 9 расстановок по 30 измерений, всего 270, Nr = 270.
Модель среды имитирует две скоростные аномалии верхней части разреза, задаваемые формулой (координаты в м, скорость в м/c).
Va(x,z) = 100-(exp(-3-10"6-(x - 1 000)2 --2-10"5-(z - 100)2) + exp(-3-10-6-(x -1 000)2 - 2-10"5-(z - 100)2)),
скорость в опорной среде линейно зависит от глубины:
Vq(x,z) = 1000 + a -z, а = 5-10'4,
общая скорость в среде выражается как сумма скорости в опорной среде и аномальной
V (x, z) = Vo( x, z) +
Уа (x, 2).
В результате
проведенного вычислительного эксперимента получены
результаты, представленные на рис. 1.
б
в
г
Характерные черты выделяются сразу после обращения исходного поля времен. Проведено две итерации с включением получаемого на предыдущем шаге решения, которые ведут как к уточнению величины аномалий скорости, так и к детализации скоростного распределения. В ходе итераций происходит
изменение траекторий лучей, что выражается в некотором сокращении зоны
«освещенности».
В предыдущих работах [2-4] отмечалась, что метод
последовательной линеаризации (метод итераций Ньютона), при регуляризации промежуточного решения является устойчивым. В настоящей статье для численного эксперимента моделируется более сложная среда,
Рис. 1. Итеративное восстановление зон скоростной аномалии:
а - первая итерации, с начальным приближением в виде I о (х, г), б - вторая итерации, в - третья
итерация, г - восстанавливаемая среда. По вертикали глубина, по горизонтали координата профиля, в м. Скорость задана полутонами, шкала справа, в м^
400
200
600
400
включающая две разделенные зоны с повышенной скоростью. Результаты численного эксперимента показывают, что разработанный алгоритм обладает разрешающей способностью, и можно достичь удовлетворительного решения, восстановив как координаты разделенных экстремумов, так и значения истинной скорости.
Как отмечается в [4], предел достижимой точности при сходимости итерационного процесса определяется тем, что для решения прямой задачи требуется регуляризация поля медленности, полученного на предыдущей итерации. Сглаживание искажает поле медленности, вследствие чего относительное среднеквадратическое отклонение от точного решения, по численным экспериментам, как правило (в зависимости от сложности строения среды), составляет от 5 до 10 %.
Работа поддержана грантами РФФИ 07-07-00251-а, 08-07-00240-а и интеграционным проектом № 10 СО РАН.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Лаврентьев М.М., Романов В.Г. О трех линеаризированных обратных задачах для уравнений гиперболического типа. // Докл. АН СССР. - 1966. - Т. 171. - № 6. - С. 1279-1281.
2. Хогоев Е.А. Численное исследование итеративного подхода к обращению поля времен рефрагированных волн // Сборник трудов международной научной конференции «Сейсмические исследования земной коры». 23-25 ноября 2004 г. - Новосибирск, Академгородок, С. 190-195.
3. Зеркаль С.М., Хогоев Е.А. Итерационная технология сейсмотомографической диагностики на основе кинематики рефрагированных волн // Доклады РАН. - № 4 (401). -2005. - С. 526-528.
4. Хогоев Е.А. Диагностика скоростных неоднородностей по кинематике рефрагированных волн// Материалы 2-го междунар. симп. «Активный геофизический мониторинг литосферы Земли», Новосибирск, 12-16 сент. 2005 г. - Новосибирск, С. 272 -276.
© Е.А. Хогоев, 2008
