автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.16, диссертация на тему:Математическое моделирование кинематики и динамики сейсмических полей в двумерно-неоднородных средах

РОССИЙСКАЯ АКАДЕМИЯ НАУК Сибирское отделение

Институт Вычислительной Математики и Математической Геофизики

На правах рукописи

Белоносов Андрей Сергеевич

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ КИНЕМАТИКИ И ДИНАМИКИ СЕЙСМИЧЕСКИХ ПОЛЕЙ В ДВУМЕРНО-НЕОДНОРОДНЫХ СРЕДАХ

05.13.16 - применение вычислительной техники математического моделирования и математических методов в научных исследованиях

Диссертация на соискание ученой степени кандидата физико-математических наук

Научные руководители: академик РАН А.С. Алексеев д.ф.-м.н. В.А. Цецохо

Новосибирск, 1999

Оглавление

Введение ................................................... 4

Глава 1. Обратная кинематическая задача сейсмики

в 2Э постановке .................................. 15

1.1. Постановка задачи ............................... 15

1.2. Начальные данные ............................... 17

1.3. Алгоритм решения обратной кинематической задачи 19

1.3.1. Построение локально-линейной скоростной модели . . 22

1.3.2. Нахождение оптимальной линейной скорости ....... 22

1.3.2.1. Вывод явного выражения для времени между двумя точками среды с линейной скоростью .............. 23

1.3.2.2. Постановка и решение задачи минимизации ........ 24

1.3.3. Реконструкция скорости в вертикальных сечениях . . 27 1.3.3.1. Решение интегрального уравнения ................. 30

1.3.4. Построение двумерной скорости ................... 34

1.3.5. Усиленная гипотеза компенсируемости ............. 35

1.3.6. Локально одномерный метод ...................... 38

1.3.7. Сравнительные эксперименты на модельных данных 39

Глава 2. Об одном приложении метода обратной фильтрации

к линейным задачам вычислительной математики ... 44

Глава 3. Комплекс алгоритмов и программ

гладкого восполнения сеточных функций,

заданных на нерегулярных сетках в Яп ............ 60

3.1. Постановка задачи.

Предварительное описание метода ................. 60

3.1.1. Основы применения (разбиения единицы) ........... 61

3.1.2. Применение к аппроксимации ..................... 62

3.1.3. Применение к восполнению........................ 63

3.2. Алгоритм г-гладкого приближения

функций многих переменных ...................... 66

3.3. Алгоритм г-гладкого восполнения функций, заданных приближенно в узлах хаотической

сетки на плоскости ............................... 75

3.4. Применение к построению карты изолиний эффективных скоростей ........................... 84

Заключение .................................................. 89

Литература .................................................. 90

Приложения .................................................. 96

Введение

Методы математического моделирования и вычислительного эксперимента эффективно применяются в разных научных и прикладных дисциплинах.

Важную роль эти методы играют в геофизике, где объектом исследований часто оказываются недоступные визуальному наблюдению физические свойства глубинных слоев Земли и внутренние источники геофизических полей.

Надежность результатов исследований в этих случаях существенно зависит от трех факторов: реалистичности и общности физико-математической модели изучаемого объекта, от количества наблюдаемой информации и от корректности математических методов обработки этой информации.

Каждый из этих факторов проявляется в любых задачах дистанционных исследований, когда об изучаемом объекте имеется лишь косвенная информация. В геофизических задачах эти факторы имеют свою специфику.

1. Фактор реалистичности и общности модели в геофизике, в частности, в сейсмических исследованиях, приводит к необходимости рассматривать многомерные (двумерные и трехмерные) модели пространственного строения объектов. Введение слоистых моделей с усредненными характеристиками слоев эффективны лишь в сейсморазведке, где имеется большая априорная информация о среде за

у "П ^

счет использования измерении в скважинах. В сейсмологии и при глубинных сейсмических зондированиях (ГСЗ) необоснованное упро-

щение модели путем ее усреднения может приводить к ошибкам интерпретации наблюдений. В частности, это иногда приводило на начальном этапе развития метода ГСЗ (1950 - 1970 гг.) к ошибкам распознавания природы рефрагированных волн в Земной коре за счет пренебрежения малым вертикальным градиентом скоростей сейсмических волн в слоях на больших глубинах. Поэтому для повышения качества математических моделей приходится повышать их сложность сообразуясь, однако, с объемом и точностью исходных данных и с точностью применяемых алгоритмов обработки.

Переход к многомерным моделям и к соответствующим численным методам в геофизике начался в 1960-1970 годах в работах сибирских математиков и геофизиков М.М.Лаврентьева, А.С.Алексеева, В.Г.Романова [1]-[4]. Этот переход был связан с существенным развитием математики в СОАН применительно к многомерным задачам и относился к широкому кругу прикладных задач сейсмологии и сейсморазведки: к прямым кинематическим и динамическим задачам, а также к обратным кинематическим и динамическим задачам сейсмики.

В предлагаемой работе в главе 1 рассмотрена обратная кинематическая задача сейсмики в двумерной непрерывно неоднородной среде.

Если в прямой задаче сложность среды ограничена лишь возможностью ее численного представления и аппроксимации, то в обратной задаче имеются принципиальные трудности. Они связаны с вопросами единственности решения обратной задачи в пробном классе сред и, даже если единственность имеется, - с проблемой устойчивости численного решения при выполнении большого числа шагов по глубине. Целью исследования в этой главе является получение достаточно устойчивого численного метода решения обратной задачи до определенной глубины.

2. Фактор полноты исходной (практической) информации в геофизике наталкивается на существенные ограничения. Для протяженных исследуемых областей Земной коры с размерами иногда в сотни километров с неоднородностями среды размерами в десятки-сотни метров требуется иметь густую сеть регистрирующих геофизические поля приборов.

Чрезвычайно высокая цена каждой точки наблюдения заставляет искать экономичные методы аппроксимации моделей и обработки данных с целью выбора рациональных систем наблюдений. В сейсморазведке, в ГСЗ и в сейсмологии обнаружение полезных волн обычно производят путем корреляции в пространстве распространяющихся вдоль профилей "сгустков энергий" волн - волновых фаз. При этом каждая новая точка пространственной системы наблюдений требует установки (хотя бы временной) специального прибора и удорожает работы при малом шаге между приборами. Между тем, каждый уже установленный сейсмограф записывает с весьма высокой детальностью (с малым шагом дискретизации) функцию времени, описывающую волны в фиксированной точке наблюдений. Детальность записи по времени достигается без больших экономических потерь. Возникает тенденция использовать более эффективно динамические характеристики волнового поля. Было показано ранее (А.С.Алексеев [16], [2]), что использование динамической информации путем решения обратных динамических задач существенно повышает разрешающую способность и надежность сейсмических исследований.

В данной работе, в главе 2, предлагается метод выделения высокочастотных компонент волн из волнового поля, зарегистрированного на более низких доминирующих частотах, но с весьма большой точностью. Здесь наблюдается повышение разрешающей способности обработки за счет использования высокой точности регистрации

сигналов во времени.

3. Фактор корректности обработки наблюдений определяется не только реалистичностью модели, точностью алгоритмов, но и точностью представления больших массивов наблюдаемых данных. Особую сложность процедура представления данных имеет в случае больших хаотических (нерегулярных) сетей точек наблюдения. Аппроксимация многомерных полей (функций) по таким сеткам должна производиться без потери аналитических свойств этих функций.

В главе 3 диссертации излагается алгоритм аппроксимации и восполнения функций многих переменных, заданных на хаотических сетках с большого порядка.

Отмеченные выше три фактора, присутствующие в методах математического моделирования в геофизике, подробно поясняются в соответствующих разделах диссертации. На актуальных прикладных задачах иллюстрируется их содержание и численные подходы.

Перейдем к описанию содержания диссертации по главам.

В главе 1 рассматривается обратная кинематическая задача сей-смики (ОКЗ) в двумерной постановке. Физическая постановка ОКЗ, следуя [3], заключается в следующем: в области пространства х = (а^, а'з), ограниченой некоторой поверхностью 5, рассматривается волновой процесс, порожденный сосредоточенными источниками возмущений, приложенными в точках х3 поверхности 5. Волны от источников возмущений распространяются в области Б с конечной скоростью зависящей от координат точки пространства. В точках границы приборы фиксируют время пробега т{х,х3) волн от источника х3 € £ до приемника х £ Б. Требуется по времени г(ж,ж-?), где х1 хз—произвольные точки поверхности найти скорость у(х) распространения волн внутри области В.

Известно, что функция т(х,худовлетворяет внутри области И

уравнению эйконала

(1)

и условию

Г (я, xi) = 0(\х — хЦ): х —> xi .

Поэтому математическая постановка ОКЗ может быть сформулирована как задача об отыскании функции п(х), входящей в уравнение эйконала, по известному для различных точек поверхности 5 (области И) решению задачи (1), (2)

Функция t(x,xi) при фиксированной точке xi называется поверхностным годографом волны (продольной или поперечной). Таким образом, данными ОКЗ является семейство поверхностных годографов.

ОКчЗ принадлежит к числу классически некорректных задач математической физики и является обьектом исследования многих ученых, как математиков, так и геофизиков.

Первая ОКЗ была поставлена и решена в 1905-1907 г.г. немецкими учеными: математиком Г. Герглотцем [5] и геофизиком Е. Вихер-том [6]. Полученный ими результат позволял по годографу рефраги-рованной волны определять одномерный закон изменения скорости в Земле в предположении сферически-симметричного ее строения и монотонного возрастания скорости с глубиной. Для сред неоднородных только по вертикали годограф рефрагированых волн т{х, £) допускает прямое обращение в скоростной разрез v(z), если в среде отсутствуют волноводы.

C.B. Чибисовым [7], [8] для целей сейсморазведки был обобщен метод Герглотца-Вихерта. Им был впервые разработан метод интерпретации годографов с разрывами.

ос^ — i>(,^oc^ X^^^ ос^ ^ ос s

Интенсивное развитие теории и методов решения обратных задач связано с начавшимся применением методов сейсмической разведки месторождений нефти и газа. В работах Г.А. Гамбурцева, Ю.В. Ри-зниченко, И.С. Берзон, H.H. Пузырева были решены обратные кинематические задачи [9]—[13] для отраженных и головных волн в слоисто-неоднородных средах. Наиболее общим результатом можно считать метод полей времен, созданный Ю.В. Ризниченко [9].

Начиная с работ А.Н. Тихонова [14], [15], в которых был сформулирован общий подход к решению некорректных задач, различные постановки обратных задач (а так же обратных кинематических задач) были сформулированы и исследованы в работах A.C. Алексеева, М.М. Лаврентьева, В.Г. Романова, Ю.Е. Аниконова, С.П. Шишат-ского, А.Л. Бухгейма, C.B. Гольдина, Т.И. Облогиной, Т.Б. Яновской и их учеников [16]—[30].

Методы решения обратных кинематических задач можно разделить на три группы.

1. Оптимизационные методы: машинный перебор, градиентные методы минимизации функционала невязок и т.д. В основе этих методов лежит возможность эффективного решения прямой задачи. Первый из названных методов позволяет получить первоначальное представление о строении среды, хотя достижение наилучшего совпадения между рассчитанными и наблюденными данными не гарантирует истинности полученного результата. Эффективность второго существенно зависит от выбора нулевого приближения. Эти методы впервые использовались в работе [31].

2. Метод линеаризации (проектирование на конечномерное подпространство). В сейсмических исследованиях приближенная математическая модель может быть зарание подобрана на основании априорных сведений или получена приближенными методами по ин-

формации о наблюденном волновом процессе. Если предположить, что функция п(х) = ^у представлена в виде

п(х) = щ(х) + щ(х) , (3)

где п0(х) - описанная выше (подобранная) функция, а щ(х) мала по сравнению с щ{х). Тогда, обозначив через г0(ж,Жо) время побега волны вдоль луча Г(х,ж0), отвечающего скорости v0(x) = получим, что

т(х, х0) и tq(x, х0) + тг(х, х0) ,

где

Ti(x,Xq) = J rii(x)ds. (4)

г

Определение добавки rii(x) к известной функции щ(х) сводится к решению уравнения (4). Это уравнение является линейным интегральным уравнением и (4) является задачей линейной сейсмической томографии. Такая постановка, предложенная А.С.Алексеевым и М.М.Лаврентьевым, рассмотрена, например, в [4]. Таким образом, задача сводится к отысканию подинтегральной функции п1(х) через известные интегралы от неё вдоль заданного семейства кривых Г(^,ж0). Преимущество описанного метода состоит в том, что для решения линейных задач разработаны вычислительные методы и эффективные алгоритмы. Один из них - проектирование бесконечного пространства моделей на специально выбранное конечномерное подпространство. Далее алгоритм расчета зависит от конкретного сочетания щ(х) и щ(х), [20]—[21], [32]—[35].

Замечание. Если в (3) щ = const, mo лучи P(rr, Xq) - прямые линии. Такие постановки используются в медицинской томографии. В 80-х годах обратные кинематические задачи стали называть задачами "сейсмической томографии" [32]-[35]. Сходство

терминологии с медицинской связано с тем, что постановки задач и системы сбора сейсмических данных близки к используемым в медицине. Однако имеется существенное различие между этими задачами. Задачи сейсмической томографии являются нелинейными; луч зависит от реконструируемой скорости; лучи криволинейны и образуют нерегулярное семейство; практически невозможно создать системы наблюдений, обеспечивающие полное покрытие лучами области реконструкции скорости. Кроме того, обратные кинематические задачи принадлежат к классу некорректных задач и в этом их особенность.

3. Численный алгоритм решения двумерной ОКЗ в точной постановке впервые был предложен в работе [30]. Он основан на разностной аппроксимации задачи Коши для нелинейного эволюционного уравнения в частных производных первого порядка относительно функции г(х1,х2, г) ~ минимального времени пробега возмущения между точками (х1,г) и (х2,г), расположенными на одной глубине

г, = ±фг\хъ г) - (г^)2 ± г) - (т^ (5)

где

п(хъ г) = -Тх1(х1, хъ п(х2} г) = т'Х2(х2, хъ г) .

Позднее для решения этого уравнения М.Е. Романов [36] использовал метод характеристик. В работах [37]-[38] алгоритм модифицирован для случая линеаризованной постановки.

Алгоритмы третьей группы основаны на использовании эволюционного уравнения (5) для послойного пересчета времени т{х^х2, г) с глубины г = ¿к на глубину 2 = гк+1 = + к на основе конечно-разностной аппроксимации, либо методом характеристик. После каждого пересчета скорость V на новой глубине вычисляется по г для близких точек, когда лучи практически прямые. В алгоритмах этого

типа приходится бороться с существенной вычислительной неустойчивостью, суть которой состоит в том, что скорость на глубине в точке (ж, х) вычисляется как отношение ПРИ х\ и х2 близ-

ких к х (обычно х — (х1 + /2), при этом величина получается из некоторого т(х1,х2,г) - времени вдоль большого луча, проходящего через точки (х1,г) и (х2,2), путем многократного вычитания малых величин (в процессе послойного пересчета). В результате этой процедуры величина т(х[, х2, г) при больших 2 может оказаться состоящей только из ошибок округления.

Наш метод относится к классу методов, занимающих промежуточное положение между указанными (первым и третьим). В нем присутствуют и пересчет функции т(жь х2: и оптимизация. Пересчет, однако, осуществляется на глубину (Нг - шаг дискретизации уравнения (5) по переменной обеспечивающий аппроксимацию этого уравнения). Скорость же в слое находится с привлечением процедуры оптимизации по подобластям с горизонтальным размером Д <С Ь, где Ь - горизонтальный размер области, в которой расположены источники и приемники. Полученные локальные приближения по какому-либо методу склеиваются в функцию, принимаемую за искомую скорость в слое.

Если в качестве локальных приближений в г-ж подобласти (г = 1,2,..., А^) ограничиться нулевыми приближениями в виде линейных скоростей

г)г-(ж, х) = а{х + Ь¡2 + сг- ,

где аг-, 6?:, сг- определяются оптимизационным методом, то мы получим предварительное грубое представление о предлагаемом алгоритме ([39], [40]).

С учетом устойчивости алгоритмов, установленной на модельных примерах, решена практическая задача интерпретации серии годографов рефрагированных волн, описанная в приложении А.

В главе 2 для некоторого класса линейных задач математической физики предложен метод численного формирования "высокочастотных" решений на основе изв