CC BY
38
УДК 518.4+517.421.1
О. Н. Белоусова
Институт автоматики и электрометрии СО РАН пр. Акад. Коптюга, 1, Новосибирск, 630090, Россия
E-mail: [email protected]
ИНТЕРАКТИВНОЕ МОДЕЛИРОВАНИЕ НА СЕРВЕРЕ «МАТЕМАТИЧЕСКИЕ ПРОБЛЕМЫ ГЕОФИЗИКИ»
В данной работе представлены результаты исследований автора в области использования Интернет-ресурсов в образовательных и аналитических информационных системах, в том числе и разработки, связанные с созданием демонстрационных версий более крупных программных продуктов, претворяющих математическое моделирование некоторых важных геофизических задач.
Ключевые слова: математическое моделирование, алгоритмы, интерактивное моделирование, сервер, вычислительный эксперимент, имитационные модели, прямая и обратная кинематические задачи геофизики, динамический анализ.
Создание, разработка, исследование и развитие моделей процессов и явлений окружающего нас мира составляют основу научного поиска и прогресса в науке. Моделирование (область знаний, занимающаяся теорией и использованием моделей) тесно связано в своем развитии с информатикой и является одним из ее основных инструментов. Очень важно, чтобы описание объекта исследований было выполнено максимально точно и не допускало принципиальных разночтений. Этим требованиям в наибольшей степени отвечают математические модели. В случае, когда математическое моделирование затруднено или вообще невозможно, прибегают к имитационным (моделирование динамических объектов и явлений) и текстовым моделям, использующим по возможности однозначные понятия.
Интерактивное математическое моделирование тесно связано с вычислительным экспериментом. Одно из определений описывает интерактивное (диалоговое) взаимодействие «человек - компьютерная программа» как взаимодействие, при котором сообщение связано с множеством предыдущих сообщений и с отношениями между ними. Одним их самых сложных уровней интерактивности традиционно считается уровень физической модели объекта, при котором система позволяет пользователю изменять те или иные параметры, и при этом меняется поведение системы. Тем не менее современный уровень развития программного обеспечения для распределенных систем позволяет создавать программные комплексы, реализующие численное компьютерное моделирование в режиме реального времени, обладающие развитым и гибким пользовательским интерфейсом, причем результат моделирования, как правило, допускает визуализацию и служит важным элементом в наглядном представлении научных знаний.
В статье приводятся важные для геофизики интерактивные численные модели кинематической сейсмотомографии и динамического анализа [1; 2] в теории распространения волн. Они интегрированы в созданный в 2006 г. при поддержке Сибирского отделения Российской академии наук (междисциплинарные проекты № 10, 35) и РФФИ (проект № 07-07-00251) и развивающийся по настоящее время информационно-аналитический сервер «Математические проблемы геофизики» (http://emf.ru/mpg) и являются существенным методическим аспектом в области получения знаний и в образовательном процессе по соответствующим специальностям.
Численное моделирование прямой кинематической задачи
Постановка задачи. Данная программа численного моделирования реализует решение прямой кинематической задачи с использованием системы дифференциальных уравнений
ISSN 1818-7900. Вестник НГУ. Серия: Информационные технологии. 2008. Том 6, выпуск 2 © О. Н. Белоусова, 2008
луча в случае двух пространственных переменных, выписанной в 1967 г. в работе [3] как следствие из уравнения эйконала. Решение ищется понижением порядка дифференциальных уравнений (со второго до первого) введением дополнительных переменных, что приводит к увеличению числа уравнений в системе, в которую добавляется еще одно уравнение для расчета времени, после чего применяется метод Рунге-Кутта. Вычисления выполняются согласно методике, подробно описанной в монографии [4]. В приведенной ниже численной иллюстрации решения прямой кинематической задачи она решается в общей постановке, которая при решении обратной кинематической задачи в томографической постановке может быть адаптирована к используемой круговой системе наблюдения.
Программа реализует решение прямой задачи геометрической сейсмики для среды, в которой скорость определяется как
V = 1 + Ахх + Ауу + А
Траектория луча определяется из решения системы нелинейных дифференциальных уравнений:
„ дlnV , , дln V , , . 2 л.дlnV х =-х у +-х z + (х -1)-,
ду dz дх
, , д ln V , , д ln V , , , ,2 1Л д ln V
у = ^^ху zу + (у -1)^—,
дх дz ду
, , д ln V , д ln V , ( ,2 1) д ln V
z =-х z +--у z + (z -1)-.
дх ду дz
Здесь х = х(^), у = у(^), z = z(s), х' = ёх / ds , у' = ёу / ds, z' = dz / ds, s - расстояние вдоль луча. При этом начальные данные:
х^) = х10, у^о) = у10, zi(so) = zi0,
х'(s0) = cos ф, у' (s0) = cos Р, z '(s0) = cos Y, s0 = {х, у, z},
где ф, p, у - направляющие углы выхода луча из точки s0.
Для вычисления в программе задаются следующие параметры в соответствующих окошках:
х, у, z - точки выхода луча из источника (0 < х < 4; - 0,5 <у < 0,5; 0 < z < 1); Ах, Ау, Az - коэффициенты глобального градиента скорости (0 < Ах < 0,2; 0 < Ау < 0,5; 0,5 < Аг < 2);
a, b, c - координаты вектора выхода луча из точки s0, после операции «Normalize» получаем величины
a' = a /V a2 + b2 + c2,
b' = b/4a2 + b2 + c2,
c' = c/4a2 + b2 + c2; H - шаг в методе Рунге-Кутта (0,001 < H < 0,1).
Описание демоверсии программы. Демонстрационная версия программы [5] написана на языке Java и запускается на компьютере пользователя в любом браузере, например Internet Explorer (рис. 1). Программный продукт позволяет задавать в окошках интерфейса численные значения следующих параметров: координаты источника (X, Y, Z), значения амплитуд (Ax, Ay, Az), направляющих косинусов (a, b, c) и H. В зависимости от заданных параметров на экране происходит отрисовка траектории луча в двух графических окнах по различным плоскостям: XZ и XY. При этом каждый параметр может задаваться в определенных численных ограничениях, что также поясняется рядом с окошком задания значений. Работа с демо-версией программы происходит с помощью следующих кнопок:
• Normalize - нормализация величин a, b, c;
• Clear - очистка экрана;
• Start - вычисление и отрисовка луча.
3 Prnhlij? Miunnotf If>1fllfl<l1 ~xplnrc и
fci* td« view ^ff/ontes 1м1я öelp О
г Hack -> . j Stop 1 a & >f.li Kfcirri» Sewnh Fn ■М 3 i.iirjpi Kislrjry Uflil Sue Рпгя s Trlt . iäl ■ QtKUU
Afidre« jej JiUp /fa'seivermcih n KniflPP/CT/gira/seiSMOj№OflL,'l HTM 3 Г1»
«Г~ -1 Storch (\wbü äÖtietore ince " CSShppf ling * фТюг*« фьпгпев
*
X[IHW4J 1' F 5 у
Y (-0 S<YiO 5} F
i:i )
AifOiArtO?) t
Ay(0<Ay<0.5f jo 4
At (ОЯАиД 1' 7
в |0.ei00Ce£32«lKl?5 /
b ||l иЧ'ЛУНЬЛЧУ ЭП.У I /
с jD5H67SS5iee^70g
H(GQ6l*Hcl) JoCDl ___ г-
NaimoSie EI I
амг I.....'
Normalize нормализация величин a.b.с
С tear очистка экрана
Start вычисление и отрисовка луча 'I
* «i«>i
Рис. 1. Пример работы демоверсии программы «DirectProblems» в браузере Internet Explorer
Численное моделирование обратной кинематической задачи в томографической постановке
Постановка задачи. Рассмотрим трехмерно-неоднородную среду с показателем преломления п(х, у, г) = у-\х, у, z), где У(х, у, г) - скорость распространения колебаний в среде. В точке «0 = (хо, Уо, 0) генерируется сигнал, а в точке « = (хь у1, 0) регистрируется время прихода рефрагированной волны - Г(«0, «1). Обратная кинематическая задача заключается в определении функции У(х, у, г) по заданной функции Т(«0, «1). В общем случае обратная кинематическая задача является переопределенной: по функции четырех переменных Т(х0, у0, х1, у1) определяется функция трех переменных У(х, у, г). В переопределенной постановке обратная кинематическая задача исследовалась В. Г. Романовым [6] и Р. Г. Мухометовым [7].
Организация системы наблюдений в форме окружности с радиусом г (с центром в точке (0, 0)) позволяет снять переопределенность задачи, так как теперь Т - функция радиуса г и двух углов ф1, фг (в полярных координатах) на источник и приемник соответственно. Исследуемая среда предполагается регулярной. Это означает, что изменения скорости в среде таковы, что паре точек источник-приемник («0 и «1) соответствует одна геодезическая линия (луч) Г(«0, «1). Следующим важным предположением является то, что скорость У(х, у, г) представима в виде
V (х ) = V,( z) + х ), х = (х, у, z),
V0 >> |V1|, V)(z) = А + Bz, А = const, B = const, A > 0, B > 0.
Соотношение величин А и B таково, что обеспечивается достаточное заглубление луча при заданной базе наблюдения (расстояние между источником и приемником). Функция V0(z) считается известной, т. е. числа А и B заданы, определению подлежит функция ^(х, у, z).
Для дальнейших рассуждений, используя метод линеаризации обратной кинематической задачи для многомерных сред, систематически применяемый начиная с работ [8-10], приходим к формуле
Tl(So,S1) = { n1ds,
Г>( S>A)
где n =---, T1 = T - T0. Значения Т считаются известными (результат решения прямой
V V0
задачи, а на практике это вектор измерений), значения T0 в случае V0(z) = А + Bz вычисляются в явном виде, а Г0 - дуги окружностей, лежащих на сфере R2 = (z--)2 + |х| ,
B
х = (х, у) . На данном этапе задача свелась к определению по функции T1 функции n1 из интегрального уравнения.
Использование системы наблюдений в виде окружности не только снимает переопределенность, но и существенно формирует томографическую постановку исследуемой задачи. Лучи Г0, «натянутые» на окружность системы наблюдений, образуют поверхность шарового сегмента (рис. 2). Изменение радиуса r позволяет получить систему вложенных шаровых сегментов, заполняющих объем исследуемой области в R3. Определяя n1 на поверхности таких шаровых сегментов, получаем решение трехмерной задачи. Заметим, что используемое послойное изучение объекта исследования и методика снятия проекционных данных ставят решаемую задачу в один ряд с известными задачами классической томографии.
Методика вычислительного эксперимента. Решению обратной задачи предшествует решение прямой, в данном случае краевой (двухточечной) кинематической задачи. Двухточечная прямая кинематическая задача решается методом пристрелки в комбинации с интерполяцией времен, рассчитанных для трех точек в окрестности приемника в точку приема. Используется многократное решение задачи Коши для системы дифференциальных уравнений луча с выбором начальных данных, основанном на методе секущих. Алгоритм решения такой задачи подробно описан в монографии [9].
В результате в соответствии с использованием метода вычислительной томографии формируется проекционная матрица, элементы которой являются временами прихода рефрагированных волн с соответствующими весами, определяемыми переходом от интегрирования вдоль луча Г0 дуги окружности к стягивающей эту дугу хорде, что осуществляет переход к двумерной задаче [9]. Далее полученный на поверхности круг разбивается на пиксели, элементы проекционной матрицы представляются как лучевые суммы, и, применяя метод алгебраической ре-конструкции, получаем искомое решение томографической задачи (в квадрате, описанном вокруг окружности системы наблюдения) и тем самым исходной обратной кинематической задачи. Используется фиксированная система наблюдений с 9-ю парами источник-приемник с 8-ю направлениями, угол проекций 45 °. Сетка восстановления - фиксированная.
Рис. 2. Система сбора данных, используемых при решении обратной кинематической задачи в томографической постановке
Демонстрационная версия программы, реализующая решение обратной кинематической задачи, написана на Java и запускается на компьютере пользователя в любой программе-браузере (например, Internet Explorer) (рис. 3).
Рис. 3. Иллюстрация работы демоверсии программы «InverseProblems», реализующей численное решение обратной кинематической задачи в томографической постановке (в Internet Explorer)
Описание демоверсии программы. Демонстрационная версия данной программы [11], как и предыдущая, написана на языке Java и запускается на компьютере пользователя в любом браузере, например Internet Explorer. Программный продукт позволяет задавать в окошках интерфейса численные значения параметров для референтной среды (объединены словом «Global»), и для локальной неоднородности (объединены словом «Local»).
При обработке данных необходимо обязательно произвести следующие операции (соответствующие кнопки интерфейса):
• Solution - решение системы линейных алгебраических уравнений (СЛАУ), полу-
ченных из предыдущих операций;
• Median Filter - медианная фильтрация;
• Average - фильтрация с помощью скользящего среднего (сглаживание);
• Interpolation - интерполяция для визуализации.
Также в интерфейсе программы предусмотрены следующие операции (соответствующие кнопки):
BVP
Cycle Projection
насчитывает одну проекцию, при этом в Projection задается любое значение угла проекции; просвечивание с разных сторон;
угол наклона всей системы (все насчитывается в среде с аномальным включением - скоростная неоднородность, затем насчитывается в среде без включения; считывание всех управляющих параметров, заданных пользователем в окошках программы);
• Parametrs - раздел задания параметров включает как Global (параметры
референтной среды), так и Local (параметры локальной неоднородности);
• Correct decision - вычисление точного решения;
• In left (right) panel - выбор левой (правой) графической панели, на которой будет
происходить визуализация точного и полученного численного решений;
• Map image - визуализация точного и численного решений.
Численное моделирование расчета синтетических сейсмограмм
в условиях горизонтально-слоистой среды
Постановка задачи. Для изучения особенностей глубинного строения земной коры и свойств глубоко погруженных геологических объектов, влияющих на формирование залежей полезных ископаемых, в настоящее время активно проводятся сверхглубинные сейсмические исследования методом общей средней точки (ОСТ) по сети опорных профилей, охватывающих всю территорию России. В волновой картине глубинного разреза ОСТ явление интерференции элементарных волн становится определяющим в сейсмической трассе. В качестве объекта для изучения эффекта затухания амплитуд сейсмических волн в земной коре предполагается использовать локальный волновой пакет, полученный в результате интерференции отраженных, дифрагированных и рассеянных волн на локальной неоднородности геологической среды и обладающий устойчивой (повторяющейся) формой в пределах этой неоднородности. Традиционно при динамическом анализе сейсмических разрезов ОСТ используются фрагменты сейсмической трассы или набора трасс, которые, не являясь формой полного сейсмического отклика от локального участка среды, не позволяют извлекать достаточную полезную информацию, особенно в условиях сложной геологии земной коры. Локальная форма волнового пакета представляет собой существенно иную по смыслу характеристику отражающей способности среды и является обобщенной формой сейсмического колебания, полученной в результате регуляризации форм отдельных волновых пакетов (синфазного суммирования с последующим осреднением), что повышает устойчивость и достоверность установления необходимых зависимостей динамических характеристик сейсмического отклика от локальных свойств неоднородного объекта гетерогенной среды.
Метод и вычислительный алгоритм. В качестве одного из возможных путей решения обратной динамической задачи сейсмики для сверхглубинных данных ОСТ предлагается метод получения количественной оценки пространственно-зависимого затухания энергии волн в земной коре на основе исследования амплитудных спектров локальных волновых пакетов. Метод основан на использовании следующего алгоритма вычисления локальной формы волнового пакета, действующего в пределах локального участка D сейсмического разреза ОСТ с характерным размером L. Алгоритм осуществляет вычисление формы локального волнового пакета Sloc(t) путем синфазного суммирования отдельных волновых пакетов Sn(t), выделенных на соответствующих фрагментах трасс. При этом признаком отдельного волнового пакета на трассе ОСТ является неразделимая интерференция отраженных, дифрагированных и рассеянных волн, имеющая начало и конец на сейсмической трассе. Форма отдельных волновых пакетов Sn(t) непрерывно изменяется от точки к точке на сверхглубинном сейсмическом разрезе ОСТ под влиянием случайных микро- и мезомасштабных неоднородностей гетерогенной среды. В этих условиях устойчивостью и постоянством формы обладает локальный волновой пакет Sloc(t), вычисленный в режиме «бегущего окна» на разрезе ОСТ с размером L, равным наблюдаемым среднемасштабным неоднородностям сейсмического волнового поля. Метод реализован в оригинальном специализированном вычислительном комплексе.
Описание демоверсии программы. Демонстрационная версия программы «DynamicSeis» (рис. 4) предназначена для расчета синтетических сейсмограмм в условиях горизонтально-слоистой среды.
Для изучения конкретных геологических объектов применяют специальные системы наблюдения, ориентированные на эффективное выделение определенных типов волн. Про-
грамма «Вупашю8е18» позволяет моделировать центральную систему наблюдений, используемую в методе общей глубинной точки (ОГТ), при глубинном сейсмическом зондировании (ГСЗ), а также в методе КМПВ.
Рис. 4. Скрин-шот работающей в Интернет демоверсии программы «Бупаш1с8е18»
Панель «Параметры модели» рабочего окна предназначена для задания конкретной расстановки пунктов приема наблюдений. Она позволяет задать количество слоев модели на полупространстве (число слоев +1), количество пунктов приема 1111 (нечетное число), интервал между пунктами приема в метрах, шаг дискретизации сейсмограммы в миллисекундах. Принято, что в симметричной расстановке взрыв располагается в центре. Для дальнейших расчетов необходимо также определить форму падающей волны путем выбора одного из трех предлагаемых в списке затухающих «сигналов», различающихся длительностью: один период, полтора периода, два периода. В результирующей сейсмограмме присутствуют следующие типы волн: прямая поверхностная, продольные - головная и отраженная.
Панель задания модели среды содержит таблицу «Задайте скорости, мощности и частоты». Заполнение таблицы позволяет конструировать различные физико-геологические ситуации, возникающие в условиях плоскопараллельных сред, с помощью задания мощности слоев, их пластовой скорости и основной частоты отклика на вызванные сейсмические колебания. По умолчанию таблица содержит один слой с пластовой скоростью 2 000 м/с и основной частотой 30 Гц. Вторая строка таблицы при этом описывает свойства нижележащего полупространства: V = 6 000 м/с, ^ = 20 Гц.
Панель «Параметры визуализации» дает возможность определить название расчетной модели для экрана и печати, а также задать, какие типы волн включать в сейсмограмму (прямая поверхностная, продольные - головная, отраженная). Для достижения выразительного вида синтетической сейсмограммы можно выбрать автоматическую регулировку амплитуд и заливку положительной части сигнала черным цветом. Функция «Хсоп81 = 0» программы «Бупаш1с8е18» позволяет создавать синтетическую сейсмограмму нулевых удалений.
Таким образом, в программе «Бупаш1с8е18» динамические параметры волн регулируются скоростными параметрами слоев, их толщиной, глубиной удаления от поверхности, формой «возбуждаемого» импульса, его амплитудой, частотой и длительностью. После задания всех необходимых параметров модели нажатие на кнопку «ОК» позволяет увидеть результат моделирования - синтетическую сейсмограмму в окне визуализации программы.
Заключение
Результаты данной работы - это информационно-вычислительные модели, развивающие подход к решению геофизических проблем, основанный на теории обратных задач сейсмической диагностики геологических сред. Проведенные исследования открывают новые возможности применения информационно-вычислительных технологий для более глубокого анализа данных сейсмического эксперимента. Системное представление полученных знаний на информационном пространстве разрабатываемых Интернет-проектов значительно повышает эффективность их использования и расширяет круг внедрения этих знаний.
Список литературы
1. Гошко Е. Ю, Зеркаль С. М., Хогоев Е. А. Вычислительная томография и динамический анализ в сейсмике. Новосибирск, 2007. 171 с.
2. Гошко Е. Ю., Зеркаль С. М. Алгоритмическое обеспечение исследования гетерогенных сред земной коры по динамическим характеристикам локальных волновых пакетов // Вестн. Новосиб. гос. ун-та. Серия: Информационные технологии. 2006. Т. 3, вып. 1. С. 33-44.
3. Белоносова А. В., Алексеев А. С. Об одной постановке обратной кинематической задачи сейсмики для двумерной неоднородной среды // Некоторые методы и алгоритмы интерпретации геофизических данных. М., 1967. С. 137-154.
4. Lavrentiev M. M., Zerkal S. M., Trofmov O. E. Computer Modelling in Tomography and Ill-Posed Problems. Inverse and Ill-Posed Problems Series, The Netherlands, VSP. 2001.128 p.
5. Белоусова О. Н, Кисленко Н. П., Хогоев Е. А. Развитие методической составляющей раздела «Геотомография» сервера «Методы решения условно-корректных задач» с использованием демоверсий численного решения томографических задач // Тез. докл. 61-й научно-технической конференции. Новосибирск: НГАСУ, 2004. С. 126-127.
6. Романов В. Г. Обратные задачи для дифференциальных уравнений. Новосибирск, 1978. 88 с.
7. Мухометов Р. Г. Обратная кинематическая задача сейсмики на плоскости // Математические проблемы геофизики / Под ред. М. М. Лаврентьева, А. С. Алексеева. Новосибирск, 1975. С.243-254.
8. Белоносова А. В., Алексеев А. С. Об одной постановке обратной кинематической задачи сейсмики для двумерной неоднородной среды // Некоторые методы и алгоритмы интерпретации геофизических данных. М., 1967. С. 137-154.
9. Lavrentiev M. M., Zerkal S. M., Trofimov O. E.. Computer Modelling in Tomography and Ill-Posed Problems. Inverse and Ill-Posed Problems Series, The Netherlands: VSP, 2001. 128 p.
10. Алексеев А. С., Лаврентьев М. М., Мухометов Р. Г., Романов В. Г. Численный метод решения трехмерной обратной кинематической задачи сейсмики // Математические проблемы геофизики / ВЦ СО АН СССР. Новосибирск, 1969. Вып. 1. С. 179-201.
11. Белоусова О. Н, КисленкоН. П., Хогоев Е. А. Имитационная модель кинематической сейсмотомографии в Интернет-проекте «Геотомография» // Тр. НГАСУ. 2005. Т. 8, № 3 (33). С.16-21.
Материал поступил в редколлегию 09.07.2008
O. N. Belousova
Interactive Modelling on a Server «Mathematical Problems of Geophysics»
In the given work results of researches of the author in the field of use of Internet resources in educational and analytical information systems, including the workings out connected with creation of demonstration versions of larger software products, some important geophysical problems realising mathematical modelling are represented.
Keywords: mathematical modelling, algorithms, interactive modelling, server, computing experiment, imitating models, direct and inverse kinematic problems of geophysics, dynamic analysis.
