gazya.ru страница 1
скачать файл

А.Л. Карчевский, К.Т. Искаков, Б.Б. Шолпанбаев


Метод послойного пересчёта для решения обратных задач геофизики


(Евразийский национальный университет им. Л.Н.Гумилёва, г.Астана)
О
братные задачи играют огромную роль в математическом моделировании и интерпретации наблюдаемых данных. Что такое обратная задача? В общем виде это можно пояснить схемой:

Т.е., если известна причина и надо установить следствия, то это есть прямая задача, если известны следствия и надо установить причину этих следствий, то это есть обратная задача.

Математически можно обратную задачу сформулировать как

Aq=g.


На практике очень часто обратная задача решается при помощи минимизации функционала невязки:

У этого способа решения обратной задачи есть один недостаток: скорость решения обратной задачи зависит от скорости решения прямой задачи, потому что во время минимизации функционала невязки приходится многократно (часто несколько тысяч раз) решать прямую задачу. Какие есть пути решения данной проблемы? Это



  • создание суперкомпьютеров,

  • распараллеливание алгоритмов,

  • создание быстрых алгоритмов решения прямых задач.

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

Одним из первых технологичных алгоритмов послойного пересчета для решения дифференциального уравения второго порядка для горизонтально-слоистой однородной среды был алгоритм Тихонова А.Н. и Шахсуварова Д.Н. [1]. Однако он имел ограничения: при его численной реализации существовали выражения, записанные с участием экспонент, имеющих показатели с положительными действительными частями, что приводило к накоплению ошибок округления при послойном пересчете.

Далее, идея послойного пересчета была реализована в следующем виде. Хорошо известно, что дифференциальное уравнение или система дифференциальных уравнений второго порядка могут быть сведены с помощью специальной замены функций к дифференциальному уравнению Риккати или дифференциальному матричному уравнению Риккати. Замечательным является тот факт, что когда коэффициенты дифференциального уравнения Риккати или дифференциального матричного уравнения Риккати являются постоянными, тогда они имеют решения, которые могут быть записаны в аналитическом виде. По всей видимости, впервые для построения численных алгоритмов, активно применяемых в геофизической практике, этот прием был использован в работе Дмитриева В.И. [2] для дифференциального уравнения второго порядка для решения прямой задачи электроразведки. Теперь уже стало общепризнано [3], что метод послойного пересчета является наиболее подходящим при численном решении обратной задачи при помощи метода минимизации функционала невязки в случае, когда среда является горизонтально-слоистой и однородной.

Процесс послойного пересчета имеет также название метода «прогонки», который описан в работе Гельфанда И.М. и Локуциевского О.В. [4].

Идеи работы Дмитриева В.И. [2] развивались далее им и его соавторами для других задач электродинамики (см., например, [5]). Для задач теории упругости эта методика была применена в работах Аккуратова Г.В. и Дмитриева В.И. [6,7] для получения следа решения системы дифференциальных уравнений теории упругости на поверхности z=0. Авторы рассмотрели систему дифференциальный уравнений теории упругости для продольных и поперечных смещений для горизонтально-слоистой изотропной среды. В данном случае уже системы дифференциальных уравнений сводилась к дифференциальному матричному уравнению Риккати. Замечательно то, что дифференциальное матричное уравнение Риккати, у которого матрицы-коэффициенты постоянны, также имеет решения, которые могут быть выписаны в аналитическом виде.

Далее свое применение для системы дифференциальных уравнений теории упругости алгоритм послойного перечета получил в работах Фатьянова А.Г. и Михайленко Б.Г. [8-10]. В работе [8,9] рассуждения ведутся для получения следа решения системы дифференциальных уравнений теории упругости для изотропной среды с поглощением, в работе [10] - для трансверсально-изотропной среды, когда ось бесконечной симметрии направлена по оси Oz.

В работах [11,12] представлен алгоритм для решения системы дифференциальных уравнений теории упругости, а в работе [13] данная методика применена для уравнений Максвелла для горизонтально-слоистых сред любого вида анизотропии.
Метод послойного пересчёта для решения дифференциального уравнения
второго порядка

Рассмотрим среду - n-слойную структуру с границами раздела ; m-ый слой находится в интервале , последний N+1 (подстилающий) слой есть . Физические свойства каждого слоя характеризуются величинами rk , то есть r - кусочно-постоянная функция переменной z , 0 < z < ∞.

Опишем метод послойного пересчета, использующий переход к дифференциальному уравнению Риккати, на следующем примере. Имеем следующую прямую задачу:

(1)

(2)

(3)

(4)

Здесь точка z* - это точка, в которой находится источник, . Будем считать, что , это никак не ограничивает нас, это предположение сделано для простоты и определённости.

Функция r(z) является кусочно-постоянной функцией. Будем обозначать rk - значения функции r(z) в интервале , точки zk - точки разрыва среды, , значение rN+1 будет соответствовать значению функции r(z) в полупространстве . Здесь использовано обозначение для значения скачка функции u(z) в точке zk.

Введем функции x(z) и s(z) следующими равенствами:





(5)

Подстановка (5) в (1) приведёт нас к дифференциальным уравнениям Риккати





(6)

Действительно,



Известно [13], что дифференциальное уравнений Риккати имеет три решения. Дифференциальное уравнений Риккати с постоянными коэффициентами замечательно тем, что имеет решения, которые могут быть выписаны в аналитическом виде. В частности, два решения уравнения Риккати в каждом интервале имеют вид:



где . Найдём третье решение. Введём новую функцию l(z) соотношением:

x(z)=l(z)+rk,

следовательно, l(z) удовлетворяет уравнению



Введём функцию



,

которая удовлетворяет уравнению

w'-2rkw=1.

Это уравнение решается методом вариации произвольной постоянной, т.е.



,

откуда


Находим c(z) , потом w(z) , l(z) и



Постоянную c0 находим, положив z=zk , и учитывая x(zk)=xk , следовательно,



(7)

Аналогичным способом получаем выражение для решения уравнения Риккати на интервале , если начальное условие находится в точкеzn-1:



(8)

Из условий склейки (3) получаем:



(10)

Для полупространства можно положить

x(z)=-rN+1.

Действительно, решение u(z) на полупрямой имеет вид , поскольку должно быть удовлетворено второе краевое условие (2), то .

Из (10) условий склейки следует, что известно xN, следовательно, используя (7), можем найти решение уравнения Риккати на интервале , и так далее, т.е. имеем рекуррентную формулу:


(11)


Из первого краевого условия (2) и (5) следует , следовательно, из (8) получаем



(12)

Из условий склейки (5) в точке z* следует:



откуда


(13)

т.е. получены начальные условия, чтобы решить дифференциальные уравнения (5). Рассмотрим интервал [0,z*] , из (5) имеем:



Следовательно,



(14)

Найдём u(z) на интервале [zk-1,zk] . Из (5) следует:



откуда получаем



(15)

где Для uk, положив z=zk} , напишем рекуррентную формулу






Порядок действий при нахождении u(0)

Чтобы найти u(0) действуем по следующей схеме:



  • на полупрямой положим , в силу условий склейки (9) получаем , следовательно, можем найти

  • на интервале [zk-1,zk] , поскольку знаем xk , можем найти , в силу условий склейки (9) имеем xk-1 , т.е. ведём пересчёт по рекуррентной формуле (11), следовательно, знаем x*;

  • вычисляем s* по формуле (12);

  • вычисляем по формуле (13);

  • вычисляем u(0) по формуле (14).

Если необходимо знать u(z) на интервале [zn-1,zn] продолжаем действия:

  • вычисляем по формуле (13);

  • вычисляем u1 и uk () по рекуррентным формулам (16)

  • вычисляем u(z) в точке z () по формуле (15).

Описанный выше алгоритм есть метод послойного пересчёта.
Градиент функционала невязки

Откуда взялась постановка прямой задачи (1)-(4)? Например, из уравнения акустики



Считаем, что плотность ρ=const . Сделаем преобразование Лапласа по временной переменной и преобразование Фурье по горизонтальным переменным



Здесь p=-α+iw - параметр преобразования Лапласа, ν1 и ν2 - параметры преобразования Фурье по горизонтальным переменным, , тогда и f=f(p). Можно привести и другие примеры практических постоновок задач, следующих из сейсморазведки, электроразведки, акустики и т.д..

Cчитаем, что относительно решения прямой задачи (1)-(4) известна следующая дополнительная информация:

(17)

Будем считать, что скорость v2 в первом слое и на полупрямой известна.

Обратная задача (1)-(4), (17) может решаться численно при помощи минимизации функционала невязки:

(18)

Придадим приращение функции (поскольку является кусочно-постоянной функцией, то естественно считать приращение тоже кусочно-постоянной функцией), следовательно, функция u(z) получит приращение . Если написать постановку прямой задачи для u+ и вычесть из неё постановку прямой задачи для функции u , то получим



(19)

(20)

(20)

Здесь было использовано



Рассмотрим приращение функционала невязки:



(22)

Получим постановку сопряжённой задачи. С этой целью выпишем функционал Лагранжа

L[]=J[]+Re

и рассмотрим его приращение



(23)

Рассмотрим интеграл




С учетом (24) можно переписать (23):



Приравнивая к нулю слагаемые с независимыми приращениями, получим постановку сопряжённой задачи



Рассмотрим приращение функционала невязки



Учтём


(последних два слагаемых равны нулю в силу условий склейки (21) и (28), тогда



Используя условия склейки (21) и (28), получим



тогда


Заменим и , используя (19) и (26), тогда



Получим аналитическое выражение для градиента функционала невязки



Обратим внимание, что решение сопряжённой задачи (26)-(28) может быть найдено при помощи метода послойного пересчёта, следовательно,



(29)

где функция y вводится соотношением . Нетрудно видеть, что x(z)=y(z) для всех z . Учитывая это равенство, получаем

J

где


P=

Т.о., имеем



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

Т.о. получили градиент функционала невязки для поиска неизвестной скорости в слоях.
Доказательство дифференцируемости функционала невязки
по координате точки разрыва среды и её аналитическое выражение

Зафиксируем все точки zk кроме одной zs .

Будем рассматривать кусочно-постоянную функцию ()с индексом + , если её значение в слоях совпадают со значениями функции , все точки разрывов ( k=1,...,s-1,s+1,...,N ) совпадают. Функция () в точке непрерывна, а разрыв находится в точке . Также рассмотрим кусочно-постоянную функцию () индексом - , если их значения в слоях совпадают со значениями функции , все точки разрывов ( k=1,...,s-1,s+1,...,N ) совпадают. Функция ()в точке непрерывна, а разрыв находится в точке (см. Рис.1). Соответствующие новым коэффициентам решения прямой задачи (1)-(4) снабдим теми же индексами.



Рис.2 Пример кусочно-постоянных функций с индексами + и - .

Для функционала невязки (18) нетрудно получить



Если существует



, (30)

тогда получим



(30)

Докажем существование предела (30). Это будет предел справа.

Функция есть решение следующей задачи:

Здесь введено обозначение , которое, например, для кусочно-постоянной функции принимает следующие значения:



Нетрудно видеть, что



Введём новую функцию



в этом случае можем записать следующую систему:



(36)

где

Пусть



тогда в силу условий склейки (33) получим

Решение на интервале



следовательно,



В силу условий склейки (34) имеем



Запишем


В этом равенстве можем перейти к пределу справа при :



отсюда и из (32)-(35) следует постановка прямой задачи:



Если теперь рассмотреть



и по схеме, приведённой выше, доказать существование предела



то окажется, что значения пределов (36) и (41) совпадают. Совпадение в некоторой точке левой и правой производной для непрерывной функции означает существование производной в этой точке. Т.о. доказали существование производной функционала невязки по координате точки разрыва среды и получили для неё выражение



где функция w есть решение задачи (37)-(40).

Осталось получить аналитическую формулу для вычисления этой производной. Введём новые функции следующими соотношениями:

Значения пересчитываются по рекуррентным соотношениям:





В точке , благодаря условиям склейки (39), получаем



что позволяет найти решение второго дифференциального уравнения:



Т.о. получено выражение для производной по координате точки разрыва среды. Можно объединить формулы градиента по скорости и производной по координате точки разрыва среды



и решать общую обратную задачу по определению скорости в слоях и местоположения точек разрыва среды (см., например, [15,16]).


ЛИТЕРАТУРА


  1. Тихонов А.Н., Шахсуваров Д.Н., Метод расчета электромагнитных полей, возбуждаемых переменным током в слоистых средах // Известия АН СССР, сер. Геофизическая, 1956, № 3, с. 251-254

  2. Дмитриев В.И., Общий метод расчета электромагнитного поля в слоистой среде // Вычислительные методы и программирование, 1968, т. 10, с. 55-65.

  3. Somersalo Е., Cheney M., Isaacson D., Isaacson E., Layer stripping: a direct numerical method for impedance imaging // Inverse Problems, 1991, v. 7, p. 899-926.

  4. Гельфанд И.М., Локуциевский О.В., Метод «прогонки» для решения разностных уравнений. В кн.: Годунов С.К., Рябенький В.С., Введение в теорию разностных схем. М.: Наука, 1962. С. 283-309.

  5. Дмитриев В.И., Федорова Э.А., Численные исследования электромагнитных полей в слоистых средах // Вычислительные методы и программирование, 1980, т. 32, с. 150-183.

  6. Аккуратов Г.В., Дмитриев В.И., Метод расчета поля установившихся упругих колебаний в слоистой среде // В кн.: Численные методы в геофизике. М.: МГУ, 1979, с. 3-12.

  7. Аккуратов Г.В., Дмитриев В.И., Метод расчета поля установившихся упругих колебаний в слоистой среде // Журнал вычислительной математики и математической физики, 1984, т. 24, с. 272-286.

  8. Фатьянов А.Г., Михайленко Б.Г., Метод расчета нестационарных волновых полей в неупругих слоисто-неоднородных средах // Доклады АН СССР, 1988, т. 301, с. 834-839.

  9. Фатьянов А.Г., Нестационарные сейсмические волновые поля в неоднородных анизотропнных средах с поглощением энергии. Новосибирск: Препринт ВЦ СО АН, 1989, № 857.

  10. Фатьянов А.Г., Полуаналитический метод решения прямых динамических задач в слоистых средах // Доклады АН СССР, 1990, т. 310, с. 323-327.

  11. Карчевский А.Л., Метод численного решения системы упругости для горизонтально слоистой анизотропной среды // Геология и Геофизика, 2005, т. 46, № 3, с. 339-351. (Перевод: Karchevsky A.L., A numerical solution to a system of elasticity equations for layered anisotropic media // Russian Geology and Geophysics, 2005, v. 46, n. 3, p. 339-351.)

  12. Карчевский А.Л., Прямая динамическая задача сейсмики для горизонтально-слоистых сред // Сибирские Электронные Математические Известия, 2005, т. 2, с. 23-61. (pdf-файл: http://semr.math.nsc.ru/v2/p23-61.pdf)

  13. Каpчевcкий А.Л., Аналитичеcкое pешение уpавнений Макcвелла в чаcтотной облаcти для гоpизонтально-cлоиcтыx анизотpопныx cpед // Геология и Геофизика, 2007, т. 48, № 8, с. 889-898. (Перевод: Karchevsky A.L., A frequency-domain analytical solution of Maxwell's equations for layered anisotropic media // Russian Geology and Geophysics, 2007, v. 48, n. 8, p. 689-695.)

  14. Kamke E., Diffrentiatialgleichungen Lösungsmenthoden und Lösungen. I Gewöhnliche Diffrentiatialgleichungen, 6, Verbesserte Auflage, Leipzig, 1959. (Русский перевод: Камке Э., Справочник по обыкновенным дифференциальным уравнениям. М.: Наука, 1976.)

  15. Карчевский А.Л., Горизонтально-слоистая среда: дифференцирование по координате точки разрыва среды // Технологии сейсморазведки, 2011, № 3, c. 17-22.

  16. Karchevsky A.L., Reconstruction of pressure velocities and boundaries of thin layers in thinly-stratified layers // Journal of Inverse and Ill-Posed Problems, 2010, v. 18, p. 371-388.


Карчевский А.Л., Искаков К.Т., Шолпанбаев Б.Б.

Геофизиканың кері есептерін шешуге арналған қабатты есептеу әдісі

Кері есептердің, математикалық модельдеуде және зерттелетін мәліметтер интерпретациясында маңызы өте зор. Кері есеп деген не? Оны келесі сызбамен түсіндіруге болады:




Яғни, себеп белгілі болса, онда нәтижені орнату керек. Бұл тура есеп. Егер натиже белгілі болып, ондай натиженің себептерін анықтау керек болса, онда кері есеп болады.

Математикалық тұрғыда кері есепті келесідей құруға болады:

Aq=g.


Тәжірибеде кері есеп үйлеспеушілік функционалының минимизациясы көмегімен шешіледі:

Кері есепті шешуде қолданылатын бұл тәсілдің бір кемшілігі бар: кері есепті шешу жылдамдығы тура есепті шешу жылдамдығына тәуелді, себебі, үйлеспеушілік функционалының минимизациясы кезінде тура есепті бірнеше рет (бірнеше мың рет) шешу керек болады. Бұл мәселені шешудің қандай жолдары бар? Олар:



  • суперкомпьтерлерді құру;

  • алгоритмдерді параллелдеу;

  • тура есептерді шешудің жылдам алгоритмдерін құру.


Karchevskiy A.L., Iskakov K.T., Sholpanbaev B.B.

Method layer recalculation for decision of the inverse problems of the geophysics

T
he Inverse problems play the enormous role in mathematical modeling and interpretation observed data. What is an inverse problem? In general type this possible explain scheme:

That is to say if known reason and it is necessary to install effects, that this there is straight line problem if known effects and it is necessary to install reason these effect, that this there is inverse problem.

Mathematically possible inverse problem to formulate as

Aq=g.

In practice much often inverse problem dares at minimization function:



There is one defect Beside this way of the decision of the inverse problem: velocity of the decision of the inverse problem depends on velocities of the decision of the direct problem since during minimization function happens to repeatedly (often several thousand once) solve the direct problem. What there is way of the decision given problems? This:



  • creation super computer;

  • multisequencing algorithm;

  • creation quick algorithm decisions of the direct problems.

One way of the making the quick decision of the direct problems and is dedicated to this work. The Givenned method is identified the method an layer recalculation. On our glance, this best method for decision of the differential equations and systems of the differential equations of the second order. The Method layer recalculation allows to get the expressions in such type that at recalculation with layer on layer of the mistake of the calculations are not accumulated.

One of the first technology algorithm layer of the recalculation for decision differential equation second order for horizontally-flaky uniform ambience was an algorithm Tihonov A.N. and Shahsuvarov D.N. [1]. However he had a restrictions: under his(its) numerical realization existed the expressions recorded with participation exhibitor, having factors with positive real parts that brought about accumulation rounding error at layer recalculation.
скачать файл



Смотрите также:
А. Л. Карчевский, К. Т. Искаков, Б. Б. Шолпанбаев
153.33kb.