МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ

имени М. В. Ломоносова

ФАКУЛЬТЕТ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И

КИБЕРНЕТИКИ

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

Головина Светлана Георгиевна

ПРИМЕНЕНИЕ ИНТЕГРАЛЬНЫХ УРАВНЕНИЙ В

ЧИСЛЕННЫХ МЕТОДАХ ОПРЕДЕЛЕНИЯ

ГРАНИЦ НЕОДНОРОДНЫХ СРЕД

Специальность 05.13.18 - математическое моделирование, численные методы и комплексы программ

ДИССЕРТАЦИЯ

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

Научный руководитель: доктор физико-математических наук, профессор Захаров Е.В.

Москва Содержание Введение 1 Определение границ неоднородностей по измерениям акустического поля 1.1 Численный метод решения прямой задачи для уравнения Гельмгольца........................... 1.2 Постановка обратной задачи................. 1.3 Нелинейное операторное уравнение для обратной задачи и численный метод его решения............... 1.4 Линеаризация обратной задачи. Результаты численного моделирования......................... 2 Определение контура зоны малой проницаемости в плоском слое (стационарный случай) 2.1 Прямая и обратная задача в тонком неоднородном слое.. 2.2 Интегральное уравнение для неизвестного контура.... 2.3 Вывод операторного уравнения для границы неоднородности............................... 2.4 Численный метод решения обратной задачи........ 2.5 Результаты численного моделирования........... 3 Определение неизвестных границ в неоднородном плоском слое (нестационарный случай) 3.1 Прямая и обратная задача для уравнения параболического типа............................... 3.2 Сведение краевой задачи к интегральному уравнению.. 3.3 Численный метод определения границы неоднородности. 3.4 Результаты численного моделирования........... Заключение Список литературы Введение В настоящее время быстро развивается большая область теоретических и прикладных исследований, связанная с определением неизвестных границ физических и искусственных неоднородных сред, наличия раcсеивателя (неоднородности) в среде, его формы и структуры по наблюдениям за распространением в таких средах зондирующих естественных или специально организованных полей (акустических, тепловых, электромагнитных и других) [18],[25],[32],[40].

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

Рассмотренные в диссертационной работе постановки обратных задач получили широкое распространение в сейсморазведке [16],[23],[42], инженерной геофизике [34], акустике океана, дефектоскопии, геоакустике [43],[44], физике атмосферы, медицине и многих других областях.

Одной из важных задач разведочной геофизики является поиск скопления углеводородного сырья [42] в геосреде при исследовании новых месторождений. При эксплуатации старых месторождений актуален поиск расположения зон малой проницаемости, где бурение новых скважин нецелесообразно.

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

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

Первая постановка обратных задач в акустике относится к 80-м годам девятнадцатого столетия и принадлежит Рэлею, который анализировал проблему восстановления распределения плотности в неоднородной струне на основании измерения частот ее колебаний, рассматривал низкочастотную теорию рассеяния.

Начало строгого математического исследования обратных задач акустики и сейсмики относится к первой половине 50-х годов прошлого столетия. Толчок этому дало решение обратных задач рассеяния на основе анализа уравнения Шредингера. Первые задачи этого класса исследовались в квантовой теории как задачи восстановления рассеивающего потенциала для одномерного (или сферически симметричного) случая.

Систематическое исследование этой проблемы проводилось И.М. Гельфандом, Б.М. Левитаном, М.Г.Крейном, В.А.Марченко.[15],[19],[29].

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

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

Методы решения каждой из этих обратных задач можно разбить на два направления: решение обратных задач в нелинейной постановке, когда учитывается сложная структура рассеянного поля и многократное переотражение волн [32],[52], [53] и построение линеаризованных моделей (приближение однократного рассеяния ( борновское приближение), рентгеновская томография и т.п. [34],[42]).

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

В диссертационной работе показано, что приближение однократного рассеяния (борновское приближение) хорошо моделирует лишь слабые рассеиватели и предполагает отсутствие фокусировки перерассеянного поля. Абстрактной математической моделью проблем вычислительной диагностики в данном случае служит линейное операторное уравнение I рода. Линейность модели позволяет построить численные алгоритмы, воспользовавшись хорошо разработанными методами решения линейных задач, в том числе и некорректно поставленных [21], [36],[57].

Пусть z - неизвестная характеристика исследуемой среды, u - полученные данные в результате эксперимента. Оператор A осуществляет связь где Z и U - метрические пространства. В настоящей работе показано, что решение линеаризованной обратной задачи сейсмики является первым этапом большинства итерационных алгоритмов решения обратных задач, учитывающих многократное рассеяние.

Ограниченность области измерений, погрешность экспериментальных данных приводит к некорректности соответствующих задач, для которых необходимо применять специальные методы. Принципиальная трудность в разработке методов решения таких задач состоит в том, что не выполнены стандартные в нелинейном анализе условия регулярности, что в свою очередь требует новых идей в разработке и анализе приближенных методов решения таких задач. Наибольшее распространение среди них получил метод регуляризации Тихонова [21],[36]. При этом задача определения границы неоднородности сводится к минимизации регуляризирующего функционала [37]. Для выбора параметра регуляризации обычно применяют принцип невязки [6],[37], использующий информацию о погрешности экспериментальных данных. Как показано в [57], при различных уровнях погрешности получаемые решения могут различаться между собой достаточно сильно. Повышение точности может быть достигнуто за счет учета априорной информации о форме неоднородности и последующего применения итерационного процесса, что требует значительных вычислительных затрат.

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

В первой главе диссертации рассматриваются проблемы нелинейной вычислительной диагностики н а примере задачи акустического зондирования. Строится скалярная трехмерная волновая модель для определения спектральной амплитуды акустического поля, возбуждаемого точечными гармоническими источниками в среде с локальной неоднородностью. Решены прямая задача для сферической неоднородности методом разделения переменных, прямая и обратная задача методом интегральных уравнений. Разработаны итерационные методы решения некорректной обратной задачи в борновском приближении и нелинейной постановке. Эти методы основываются на идеях итеративной регуляризации [5]-[7]. Для тестирования программы решения обратной задачи рассеяния использовались результаты решения прямой задачи [59], [60].

Вторая глава исследования посвящена изучению возможности определения зон малой проницаемости по измерениям установившегося давления в скважинах, при этом давление поперек слоя практически постоянно, за исключением некоторых зон полной непроницаемости, которые ограничены искомыми контурами. Данная задача сводится к решению обратной задачи для эллиптического уравнения. Обратная задача такого типа возникает, например, при поиске зон малой проницаемости в нефтяном пласте (т.е. зон, где бурение новых скважин нецелесообразно) по измерениям давления в имеющихся скважинах. Предполагается, что давление в нефтяном пласте не изменяется поперек пласта, а проницаемость в нем постоянна, за исключением области полной непроницаемости, размеры которой малы по сравнению с размерами пласта [61],[62]. В этой главе диссертации аналитически выписано решение прямой задачи для круга, методом интегральных уравнений решена прямая и обратная задачи. Разработан итерационный метод решения некорректной обратной задачи определения границы зоны малой проницаемости в плоском слое в линейном приближении. Для тестирования разработанного метода использовалось решение для круга выписанное в явном виде.

В третьей главе была рассмотрена задача определения границы зоны непроницаемости в плоском слое, когда в качестве входной информации для восстановления границы используются нестационарные данные о давлении в скважинах, что позволяет уменьшить количество измерений. Предложенный подход рассматривается применительно к двумерной среде, что обусловлено сложностью проведения численных экспериментов в трехмерном случае по причине больших вычислительных затрат [64],[65],[69]. Решение прямой задачи сведено к решению интегральных уравнений и в частном случае, когда неоднородность имеет форму круга, решение выписано аналитически и использовано для тестирования при численных расчетах. Предложен и реализован итерационный алгоритм решения некорректной обратной задачи, основанный на линеаризации интегральных операторов.

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

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

Основные результаты работы перечислены в заключении, опубликованы в работах [60-70], докладывались на научных семинарах и конференциях.

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

1.1 Численный метод решения прямой задачи для уравнения Гельмгольца Рассмотрим распространение акустических волн в трехмерной безграничной однородной среде, содержащей односвязную ограниченную область H R3 с достаточно гладкой границей H. Однородная среда имеет постоянное значение распределения скорости c0 - фоновое распределение скорости. Неоднородное включение описывается кусочногладкой функцией c(r) - скоростью распространения волны, r = (x, y, z), вне H c(r) = c0.

Распространение звуковых волн малой амплитуды можно описать с помощью акустического волнового уравнения:

где u(r, s, t) - поле в среде, зависящее от пространственных переменных r, s R3 и времени t 0, - оператор Лапласа, c(r) - скорость распространения волны в среде, функция источника (r s)f (t) - описывает возмущение среды точечным источником, расположенным в точке s некоторой ограниченной области S, (·) - дельта функция Дирака.

Начальные условия имеют вид:

В ряде приложений, неоднородность облучается различными точечными источниками, располагающимися в точках s R3. Источники локализованы в пределах области S R3 - области излучения. Измерение поля, рассеянного неоднородностью, может проводиться лишь в пределах области P R3 - области приема, причем замыкание области H не пересекается с замыканием областей P и S. Функция f (t), описывающая возмущение среды, известна априори (см.Рис.1.1).

Прямая задача состоит в вычислении поля в области расположения приемников P. Рассмотрим среду, в которой распространение акустических колебаний описывается источниками с временной зависимостью exp(t), тогда:

Получим уравнение Гельмгольца для поля в точке r R3 для источника в точке s R3 на частоте :

дополненное условиями излучения Зоммерфельда при | r | :

где u(r, s, ), f () - преобразования Фурье от u(r, s, t) и f (t) соответственно по переменной t, - комплексная переменная. Здесь мы ввели новый коэффициент k(r, ) предполагая, что среда однородна вне сферы достаточно большого радиуса (так как неоднородность мы считаем локальной, следовательно существует некоторый радиус R, такой, что для r(x, y, z), Приходим к граничной задаче для уравнения Гельмгольца с условиями сопряжения на границе H для функции u(r, s, ) и ее нормальных производных, где n - внешняя нормаль к границе:

[(r, s, )] = 0, [ u(r,s,) ] = 0, где [·] - обозначает разрыв значения функции на границе раздела сред.

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

Учитывая (1.6), уравнение (1.3) запишем в виде:

где коэффициент не зависит от r. В этом уравнении можно трактоc вать правую часть равенства как источник полного поля u, состоящий из источника первичного поля (r s)f () и источников рассеянного поля 2 u(r, s, )(r) - вторичных источников, наведенных в области рассеяния H полем первичных источников [18].

Полное поле в описанной выше модели можно представить как сумму первичного поля u0 (r, s, ) и вторичного (отраженного неоднородностью) поля uh (r, s, ):

которые удовлетворяют уравнению Гельмгольца с разными правыми частями:

От дифференциального уравнения (1.7) перейдем к интегральным уравнениям.

Рассмотрим функцию Грина вида:

где |r s| - расстояние между точками r, s R3, являющуюся решением уравнения (1.7) с правой частью вида (r s) и стремящуюся к 0 при Тогда, в предположении, что рассматриваемые интегралы существуют, получим интегральное уравнение:

где s S R3, r будем предполагать принадлежащим области измерений P R3, интегрирование проводится по области H, где (x) отличная от нуля функция.

Запишем уравнение (1.11) в виде:

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

Заметим, что левая часть уравнения (1.12) - поле, которое измеряется в области расположения приемников P, а интегрирование проводится по локальной области H R3, где функция (x) отлична от нуля.

Рассмотрим уравнение (1.12) отдельно для областей r H и r P, тогда вводя обозначения r = h H и r = p P запишем соответственно:

При решении прямой задачи, когда граница неоднородности H известна и известна функция (x), x H, необходимо вычислить значение функции uh (p, s, ) в области P, используя систему (1.14) при фиксированных Интегральную систему удобно решать последовательно в два этапа.

На первом этапе решается интегральное уравнение второго рода в области H - первое уравнение в системе (1.14), из которого определяется функция u(h, s, ), h H. Вычисленная функция u(h, s, ) подставляется во второе уравнение системы, откуда прямым интегрированием по H определяется искомая функция uh (p, s, ) в области P.

Для проверки построенной модели распространения акустических волн в однородной среде, содержащей локальную неоднородность, были проведены модельные расчеты для сферической неоднородности, когда прямую задачу удается решить аналитически [60]. Наличие сферической симметрии позволяет применить для решения прямой задачи метод разделения переменных и использовать представления искомых величин в виде разложения в ряд по собственным функциям волнового оператора [11], построить решения в аналитической форме и изучить основные закономерности рассматриваемых процессов, сравнить полученные результаты вычислений с функцией uh (p, s, ) из системы (1.14).

Рассмотрим неоднородность сферической формы радиуса R = a, причем внутри шара фазовая скорость c = const = c0 см.Рис. 1.2. Пусть на границе шара поставлены условия сопряжения для функции u и её норu мальной производной n :

где n - вектор нормали к границе шара H, [·] - обозначает скачок значений функций на границе раздела сред.

Таким образом мы приходим к трехмерной граничной задаче для оператора Гельмгольца вне и внутри шара:

где k0 = c0, k= c, u = u(r, s, ), s положение источника, частота излучения, A() - амплитуда.

Введем сферическую систему координат (,, ), центр которой совпадает с центром шара и ось проходит через точечный источник так, что задача становится осесимметричной по координате. Обозначим u = u (, ).

Задачу (1.15) решаем методом разделения переменных. Представим решение задачи (1.15) в виде ряда по сферическим функциям для операторов Гельмгольца:

u (, ) = где из граничных условий следуют соотношения:

и где n, n и их производные выражаются через цилиндрические функции полуцелого порядка:

n (x) = Ren (x) = Jn+1/2 (x)/ x;

n (x) = Hn+1/2 (x)/ x;

а функции n (x) связаны известными рекуррентными соотношениями:

u0 (, ) = кости = const до источника, r0 - расстояние от центра до источника, Pn (cos ) - полиномы Лежандра n-го порядка, Jn+1/2 (x) - функция Бесселя первого рода (n+1/2) порядка, Hn+1/2 (x) - функция Ханкеля первого рода (n + 1/2) порядка.

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

На рисунке (см.Рис.1.3) приведены результаты численных расчетов поля uh (p) при фиксированном источнике и одной фиксированной частоте двумя предложенными методами для неоднородности сферической формы.

Были выбраны следующие параметры и геометрия модельного эксперимента. Неоднородность радиуса 10 метров залегает на глубине метров от поверхности, внутри скорость равна 1600 м/с, а вне - м/с; источник возмущения расположен на поверхности строго над шаром и имеет частоту колебаний 400 Гц., что соответствует длине волны 4 метра внутри шара и 4.5 метра вне шара. Приемники расположены на поверхности на прямой, проходящей через точку источника на дискретной сетке из 1000 точек. Сетка покрывает отрезок [80, 80] с линейными размерами, превосходящими диаметр неоднородности в несколько раз.

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

1.2 Постановка обратной задачи Предположим, что проводилась серия возмущений однородной среды, занимающей все пространство R3 точечными источниками, расположенными в области S R3. Функция, описывающая возмущение среды, f () считается известной, будем называть ее - функцией источников.

Функцию источника можно аппроксимировать дельта-функцией. Эта аппроксимация основывается на предположении о малых линейных размерах источников по сравнению с размерами областей S, H и P. Для лись измерения излучения u(p, s, t) в области расположения приемников p P, P R3. Предполагается также, что в однородной среде с известной характеристикой c(r) = c0 = const локализована неоднородность в пределах области H R3 с характеристикой c(r) = c0, причем замыкание области H не пересекается с замыканием областей P, S. Местоположение и форма неоднородности и значение функции c(r) при r H неизвестны.(см.Рис.1.4).

В качестве изображения среды выступает скоростная функция c(r) и обратная задача сводится к нахождению коэффициента k(r, ) = c(r), характеризующего рассеивающую неоднородность, на основании измерения поля u(p, s, ) в области P. Источники считаются при этом заданными и могут варьироваться, что соответствует облучению неоднородности различными способами: с изменением частоты, ракурса облучения и т.д.

Подобная постановка обратной задачи распространена в задачах зондирования.

При решении обратной задачи имеются априорные сведения о функции скорости, например известна фоновая скорость распространения волн в среде, она может быть задана постоянной c0 = const или произвольной c0 = c0 (r). В данном случае искомая скорость c(r) = c0 + c (r), тогда решение обратной задачи сводится к определению c (r).

Разложим c2 (r) в ряд:

Таким образом, искомая характеристика неоднородности при решении обратной задачи определяется по фоновой скорости и наблюденному в области P волновому полю. Данная процедура будет устойчивой лишь в том случае, если выполняется условие c (r) c0. Ограничиваясь в разложении c2 (r) в ряд величинами первого порядка малости, получим:

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

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

Если c0 (r) - гладкая функция, а в нашем случае c0 = const, то условие (r) 1 отвечает модели слабонеоднородной среды, которая широко используется в прикладных задачах. Множество частот будем считать фиксированным и ограниченным.

1.3 Нелинейное операторное уравнение для обратной задачи и численный метод его Для решения обратной задачи нахождения неизвестной поверхности локальной неоднородности в однородной среде рассмотрим уравнение (1.12) отдельно для областей h H и p P :

где неизвестными являются функция (r) и поле вторичных источников u(r, s, ) при r H.

Разработка численных методов решения нелинейной системы (1.17) встречает ряд трудностей: данные и решение лежат в разных областях, различны размерности множества входных данных и решения, не выполнены обычные условия регулярности.

В работе [7] рассмотрен абстрактный аналог системы (1.17) в виде операторного уравнения F (x) = 0, где оператор F действует из гильбертова пространства H1 в другое гильбертово пространство H2, x = -вектор неизвестных. Для его решения использован метод Ньютона-Гаусса [2]:

xn+1 = arg min F (xn ) + Fn (x xn ) где Fn = F (xn ).

В работе [7] для решения системы (1.13) также рассмотрен итеративнорегуляризованный метод Ньютона-Гаусса, на каждом шаге которого минимизируется по x функционал:

где n -параметр регуляризации, n -некоторый элемент H1. Тогда итерационный процесс можно записать как:

В [7] приведены результаты численных расчетов с использованием предложенных методов для двумерной области, одной частоты, а искомая функция вычислялась на сетке 5x5.

При решении поставленной задачи предложенным методом возникают значительные трудности, связанные с размерностью массивов данных.

Рассмотрим оптимизированный метод решения системы (1.17) с использованием метода Ньютона-Гаусса. Учитывая специфику задачи, в предложенном методе удается при построении итерационного алгоритма уменьшить вычисления и снизить размерность используемых массивов.

Запишем систему (1.17) в операторном виде:

где интегральные операторы H и P определены следующим образом:

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

Изменение расположения источников приводит к пересчету оператора S, приемников - оператора P, изменение частот - требует вычисление трех операторов.

Выпишем схему итерационного процесса:

1 шаг: зададим начальное приближение функции 0 (r) = 0;

2 шаг: найдем функцию u0 (r, s, ) решив уравнение F1 (0 (0 ), 0 ) = используя стандартные методы;

3 шаг: решаем уравнение F2 (0 (0 ), 1 ) = 0 регуляризованным метоu дом Ньютона-Гаусса и находим 1 (r);

4 шаг: найдем функцию u1 (r, s, ) из уравнения F1 (1 (1 ), 1 ) = 0 и т.д.;

Выпишем регуляризованный метод Ньютона-Гаусса для выше указанного метода:

n+1 = n ((F2n ) ·F2n +n I)1 ((F2n ) ·F2n (n (n ), n )+n (n )), (1.20) где n - параметры регуляризации, F2n = F2 (n ). В предложенном методе на каждом шаге итерационного процесса приходится численно решать интегральное уравнение, что замедляет скорость вычисления итераций, но существенно позволяет снизить размерность обрабатываемых массивов, повысить качество вычислений.

Следует заметить, что первый шаг итерационного метода совпадает с решением обратной задачи в борновском приближении, действительно, при начальном приближении 0 = 0 из уравнения F1 (0 (0 ), 0 ) = находим u0 (0 ) = u0, тогда следует, что F2 (0, ) = P u0 U = 0.

Получили интегральное уравнение Фредгольма I рода для нахождения : P u0 (r) = U. Решение обратной задачи в борновском приближении будет подробно рассмотрен ниже.

На рисунке 1.5а приводится схема эксперимента, на верхней плоскости расположены источник и приемники, точное решение имело форму шара радиуса 10 метров. В качестве входных данных использовались результаты решения прямой задачи для сферы, выписанные в явном виде. При решении обратной задачи выше изложенным методом с нулевым приближением, исследуемая область H имела форму куба, размерностью 180 180 180. На рисунке 1.5б приводится сечение восстановленной неоднородности в плоскости OXZ. В силу симметрии задачи сечения OY Z и OXY совпадают. На рисунке 1.5в изображено добавочное поле, создаваемое неоднородностью на поверхности, на рисунке 1.5 г - график добавочного поля на поверхности вдоль оси OX. Расчеты проводились для одной частоты 400 Гц.

Далее на рисунках приводятся результаты расчетов, где в качестве входных данных использовались результаты решения прямой задачи с внесенной погрешностью. Результаты численного моделирования показали, что для одночастотного случая количество приемников должно быть велико и они должны располагаться в области большей, чем характерные размеры искомых неоднородностей.

В качестве начального приближения для всех расчетов бралось нулевое приближение 0 = 0, параметр регуляризации = 0.01. Численные расчеты показали быструю сходимость итерационного процесса. Для восстановления неизвестных неоднородностей достаточно нескольких итераций, при этом ошибка вычислений составляет менее 1%.

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

На рисунке 1.6а изображена восстановленная неоднородность предложенным итерационным методом после 10 итераций. На рисунках 1.6(б, в, г) приводятся проекции точного решения на плоскости XY, ZY и XZ соответственно. На рисунке 1.6д приведено полутоновое изображение поля u, измеренного в области приема P. Приемники Rec в количестве Nrecx Nrecy = 25 25 располагались на верхней плоскости исследуемой области Y = 18 ( рис.1.6к), при этом X и Y менялись от -1 до 19. Четыре источника имели координаты Sor1 = (6, 18, 6), Sor2 = (12, 18, 6),Sor3 = (12, 18, 6) и Sor4 = (12, 18, 12) ( рисунок 1.6и).

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

На рисунке 1.7а изображена восстановленная неоднородность после 10 итераций. На рисунках 1.7(б, в, г) приводятся проекции точного решения на плоскости XY, ZY и XZ соответственно. На рисунке 1.7д приведено полутоновое изображение добавочного поля u, измеренного в области приема P. Приемники Rec в количестве Nrecx Nrecy = располагались на верхней плоскости исследуемой области Y = 18 ( рис.1.7к), при этом X и Y менялись от -1 до 19.Четыре источника имели координаты Sor1 = (6, 18, 6), Sor2 = (12, 18, 6),Sor3 = (12, 18, 6) и Sor4 = (12, 18, 12) см. Рис 1.7и.

На рисунке 1.8 приведены результаты решения обратной задачи для случая, когда точным решением являются две неоднородности, имеющие форму куба, расположенных в исследуемой области размером см. рисунке 1.8з. Координаты одного из кубов указаны на рисунке 1.8е, второго: X1 = 10.9, X2 = 14.4 координаты Y 1, Y 2, Z1, Z3 аналогичны первой неоднородности.

На рисунке 1.8а изображены восстановленные неоднородности после 10 итераций. На рисунках 1.8(б, в, г) приводятся проекции точного решения на плоскости XY, ZY и XZ соответственно. На рисунке 1.8д приведено полутоновое изображение добавочного поля u, измеренного в области приема P. Приемники Rec в количестве Nrecx Nrecy = располагались на верхней плоскости исследуемой области (Y = 18 смотри рис.1.8к), при этом X и Y менялись от -1 до 19.Четыре источника имели координаты Sor1 = (6, 18, 6), Sor2 = (12, 18, 6),Sor3 = (12, 18, 6) и Sor4 = (12, 18, 12) см. Рис 1.8и.

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

На рисунке 1.9а изображены восстановленные неоднородности после 10 итераций. На рисунках 1.9(б, в, г) приводятся проекции точного решения на плоскости XY, ZY и XZ соответственно. На рисунке 1.9д приведено полутоновое изображение поля u, измеренного в области приема P. Приемники Rec в количестве Nrecx Nrecy = 25 25 располагались на верхней плоскости исследуемой области Y = 18 (см. рис.1.9к), при этом X и Y менялись от -1 до 19.Четыре источника имели координаты Sor1 = (6, 18, 6), Sor2 = (12, 18, 6), Sor3 = (12, 18, 6) и Sor4 = (12, 18, 12) см. Рис 1.9и. Ниже приведены результаты работы итерационного процесса, которые иллюстрируют быструю сходимость. На рисунке 1.10 изобРис. 1.10:

ражены точное решение обратной задачи (заданная неоднородность), и начальное приближение для итерационного процесса, в качестве которого бралась нулевая функция.

На рисунках 1.11 и 1.12 приводятся результаты вычислений предложенным итерационным методом. Источники располагались в точках:

(0,0,0), (200,200,0), измерения проводились на сетке приемников 25x25, покрывающих область 800x800 при z = 0. Лоцируемая область размером 200x200x200 находилась на глубине 70, в ней имелась неоднородность 90x20x90, =400Гц.

1.4 Линеаризация обратной задачи. Результаты численного моделирования Трудности исследования и решения нелинейных задач привели к тому, что в теоретических и прикладных исследованиях часто используют различные линеаризованные постановки задачи. Линеаризация задачи происходит при дополнительных требованиях, налагаемых на неизвестную функцию, которые в свою очередь, существенно сокращают класс исследуемых рассеивателей. В этом разделе проводится анализ широко используемых методов линеаризации обратной задачи: низкочастотная асимптотика для уравнения Гельмгольца и борновское приближение, анализ численных результатов.

Низкочастотная асимптотика для уравнения Гельмгольца.

Рассмотрим уравнение:

Пусть решение u(r, s, t) при t экспоненциально стремится к нулю.

От дифференциального уравнения (1.21) перейдем к интегральным уравнениям. Рассмотрим функцию Грина вида где |r s| - расстояние между точками r, s R3, являющуюся решением уравнения (1.21) с правой частью вида (r s) и стремящуюся к 0 при Так как источник поля имеет вид: f (r, s, ) = (r s)f (), то он создаст первичное поле:

Тогда, в предположении, что рассматриваемые интегралы существуют, получим интегральное уравнение:

где s S R3, r будем предполагать принадлежащим области измерений P R3, интегрирование проводится по области H, где (x) отлична от нуля. Как видно, уравнение (1.23) является нелинейным уравнением относительно неизвестных функций u(x, s, ) и (x), где x H.

Поскольку u(x, s, ) определена на полуплоскости, содержащей точку = 0, и замыкания областей P, Q не имеют общих точек с замыканием области H, то конечный предел правой части уравнения (1.23) при существует при достаточно общих условиях и равен:

Из существования конечного предела правой части следует существование конечного предела левой части уравнения (1.23), который равен:

Заметим, что приведенная задача рассматривалась для случая среды, занимающей все пространство R3. Однако, не представляет особого труда перейти к случаю среды, занимающей часть пространства R3. Для этого необходимо в приведенных выше рассуждениях рассматривать функцию Грина уравнения Гельмгольца для соответствующей части пространства. Например, в случае среды, занимающей полупространство z 0 (в декартовой системе координат Oxyz ) c граничным условием uz (r, s, t) |z=0 = 0, функция Грина имеет вид где s - точка симметричная точке s относительно плоскости z = 0.

Используя функцию Грина (1.26) и учитывая (1.25), получим следующее уравнение:

где r P, s S. Таким образом, полученное уравнение (1.27) является линейным интегральным уравнением I рода относительно неизвестной функции (x).

При проведении численных расчетов использовался регуляризированный метод простой итерации:

D, где D - максимальный диаметр области S, функция Грина (2.8) может быть записана в виде:

отсюда следует, что нормальная производная функции Грина не зависит от R0 и, следовательно, решение интегрального уравнения (2.13), т.е. плотность потенциала µ(M ; C), не зависит от радиуса внешнего контура R0. Однако расчет давления в скважине по формуле (2.12) зависит от радиуса R0, который неизвестен при решении обратной задачи.

Действительно, при больших R0 >> D0 (D0 - диаметр области, где рассчитывается давление ) функция Грина имеет вид:

В результате давление при (M, Mk ) D0, согласно (2.12), имеет вид где Присутствие множителя ln R0 существенно влияет на точность расчетов, так как при больших R0 необходимо решать интегральное уравнение (2.13) с высокой точностью. В противном случае возникает существенная погрешность.

В работе [61] предложен подход, позволяющий исключить зависимость решения задачи (2.12)-(2.13) от радиуса внешней окружности R0.

Рассмотрим предельный случай для (2.1)-(2.3):

где - первичное поле, т.е. фоновое давление в точке M с координатами (x, y) создаваемое точечным источником Mk с координатами (xk, yk ) в отсутствии зоны непроницаемости.

Будем искать решение задачи (2.19)-(2.21) в виде:

где u(M ) - давление, возникающее в результате наличия зоны непроницаемости S, а P0 - пластовое давление.

Устремим R0, тогда задачу (2.9)-(2.11) можно сформулировать, как внешнюю задачу Неймана для функции u(M ; C) :

u(M ; C) 0 на бесконечности.

Для того чтобы удовлетворить условию u(M ; C) 0 на бесконечности, решение задачи (2.22)-(2.23) будем искать в виде следующего потенциала простого слоя:

где µ(M0 ; C) - плотность потенциала, действительная функция, определенная на контуре C, точка O - произвольная фиксированная точка в области неоднородности S, не принадлежащая контуру C, расстояние (M, O) = 1. Таким образом, выбор точки O зависит лишь от области S. Выражение (2.24)является решением уравнения Лапласа и стремится Заметим, что P (M ; Mk, C)- давление при работающей k -ой скважине зависит от наличия неоднородности с границей C, а пластовое давление P0 и G(M ; Mk ) не зависят от формы локальной неоднородности.

Подставим представление (2.24 ) в граничное условие (2.23), получим интегральное уравнение для нахождения плотности потенциала µ(M ; C):

где Таким образом, при заданном контуре C решение задачи (2.9)-(2.11), а следовательно и задачи (2.1)-(2.3), сводится к решению интегрального уравнения. Подобный подход к решению задачи (2.1)-(2.3) был предложен в работе [61].

2.3 Вывод операторного уравнения для границы неоднородности Рассмотрим одну звездную область S, контур C которой задан в полярных координатах с центром в точке O:

где r - расстояние от центра неоднородности O с координатами (0, y0 ) до контура C:

Координаты точки на контуре C запишем в виде:

Выпишем вектор нормали n = { dx, 1} в полярных координатах:

Пусть точка M не лежит на контуре C и имеет координаты (, r()), точка M0 лежит на контуре C и имеет координаты (, r()), тогда уравнение (2.24) в полярных координатах примет вид:

здесь: 0 2, r() = (M, O), r() = (M0, O), l() = Аналогично, для функции µ(; r) запишем уравнение (2.25 ) в полярG(M,M0 ) где 0 2, r в обозначении µ(; r) означает, что решение получено при фиксированном контуре r(), Mk - источник с координатами (k, Rk ), а точки M и M0 лежат на контуре C и имеют координаты (, r()) и (, r()) соответственно.

Введем следующее обозначение ядра интегрального уравнения (2.30):

Будем считать, что контур C таков, что величиной r () в формуле (2.30) и (2.31) можно пренебречь, т.е. выполняется следующее условие:

Установим связь между изменением контура и изменением давления, вычисленного в точке M (, r()).

Пусть нам известно некоторое начальное приближение контура r0 ().

Линеаризуем интегральные операторы, ограничиваясь двумя членами разложения:

Тогда для функции u(; r) справедливо представление:

здесь r0 () = r() r0 () - отклонение начального приближения контура от точного решения. Для функции u(; r) ввели новое обозначение - ur (; r), показывая, что вычисление u(; r) проводится для контура r().

Уравнение (2.36) представляет собой интегральное уравнение Фредгольма первого рода относительно неизвестной функции r0 () :

где ядро:

а правая часть u (, r) = ur (; r) ur0 (; r).

Аналогично, относительно неизвестной функции µ (, r0 ) = µ(; r) µ(; r0 ) получим интегральное уравнение Фредгольма второго рода:

где 0 2, ядро T (, ; r0 ) = H(, ; r0 ) · r0 (), а правая часть:

Таким образом, уравнение (2.37) устанавливает связь между изменением контура и изменением давления в скважинах при одном источнике Mk, а уравнение (2.38) позволяет найти поправку к плотности потенциала при изменении контура [63].

Уравнения (2.37) и (2.38) можно записать в операторном виде:

Получили линейное операторное уравнение для вычисления функции r0, где Ak - оператор, устанавливающий связь между изменением давления в скважинах и изменением контура.

2.4 Численный метод решения обратной задачи Для вычисления Ak r0 () при заданном контуре r() и начальном приближении контура r0 () и k-ой работающей скважине, необходимо решить интегральное уравнение (2.38), найдя поправку к плотности потенциала простого слоя µ (; r) при изменении контура, а из уравнения (2.37) определить новую поправку r () к границе неоднородности и т.д.

Еще раз отметим, что такая постановка задачи позволяет исключить зависимость решения задачи (2.1)-(2.3) от радиуса внешней окружности R0, который не известен.

При решении обратной задачи мы учитываем, что каждая скважина становится рабочей, а остальные закрываются и в закрытых скважинах измеряется давление. Таким образом, мы имеем K(K 1) измеренных данных при различных работающих скважинах. Так как исходные данные u(; r) - измерения в скважинах, заданы приближенно, необходимо учитывать специфику уравнения (2.40), связанную с его некорректностью.

Рассмотрим операторное уравнение 1-го рода где здесь L - количество точек на контуре, K- количество пробуренных скважин, каждая из которых последовательно является источником, а в следующих измеряется давление, Ak - оператор при k-ой работающей скважине, ukj - изменение давления в j-ой скважине при k-ой работающей.

Рассмотрим сглаживающий функционал А.Н.Тихонова:

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

Используя необходимое условие экстремума функционала, получим, что если на r достигается нижняя грань M [r ], то r является решением уравнения где A - оператор, сопряженный к оператору A, тогда поправка к границе неоднородности где E-единичный оператор.

Используя свойства сопряженного оператора A и переходя в уравнении Эйлера к конечно-разностной аппроксимации, получим систему алгебраических уравнений. Для решения этой системы можно использовать стандартные численные методы.

Построим следующий итерационный процесс решения уравнения (2.41):

1. на первом этапе численного решения обратной задачи найдем начальное приближение r0 () искомого контура. Очевидно, что самым простым начальным приближением является окружность радиуса a. Задача поиска радиуса a сводится к минимизации функции одной переменной, данная задача решается стандартными методами, например, методом перебора;

2. из уравнения (2.31) для начального приближения контура r0 () вычислим плотность потенциала, а из (2.30) давление ur0 (, r), тогда u0 () = ur (; r) ur0 (; r), где ur (; r)- известная функция давления при решении обратной задачи;

3. на j - м шаге итерационного процесса вычислим поправку к контуру:

где ядро:

а правая часть uj () = ur (; r) urj (; r);

4. из следующего уравнения находим поправку к плотности потенциала µj (; rj ):

где 0 < < 2, ядро T (, ; rj ) = H(, ; rj ) · rj (), а правая часть:

Далее для вычисленной поправки µj (; rj ) найдем новую поправку к контуру rj+1 () и т.д. Если изменение давления удовлетворяет неравенству:

где - заданная точность вычислений, то искомый контур найден.

2.5 Результаты численного моделирования При исследовании результатов решения обратной задачи, в качестве входных данных использовалось вычисленное давление, являющиеся решением прямой задачи с внесенной погрешностью. Для проверки чувствительности измеряемых величин от расположения и формы неоднородности, точности решения прямой задачи, приводится анализ полученных численных результатов, который проводился для случая одной звездной области непроницаемости, контур C которой задан в полярных координатах r = r(). Тогда ядро интегрального уравнения (2.25) можно представить в виде:

H(, ; r) = где l() = После преобразования получаем ядро интегрального уравнения:

H(, ; r) = Отметим, что ядро интегрального уравнения не имеет особенности при совпадении аргументов:

Правая часть уравнения (2.25) равна M -точка с координатами (r(), ), Mk - точка с координатами (Rk, k ), где f (; r) = H(, k ; r), при r(k ) = Rk, или Учитывая,что dlM0 = l()d, можем записать интегральное уравнение (2.25) в виде:

Если ввести новую неизвестную функцию U (; r) = µ(; r)l(), получим интегральное уравнение:

Таким образом, необходимые для обратной задачи данные, вычислим по формуле:

Расчеты проводились для случая, когда контуром является эллипс:

где a, b -полуоси эллипса, эксцентриситет эллипса = b/a. При этом считалось, что точки измерения и возбуждения находятся на окружности радиуса Ru, т.е. Rm = Ru для любого m. Поэтому результаты расчета u(n) зависят от безразмерного параметра = b/a, места расположения источника k и измерения m.

На рисунке 2.2 приведены значения давления u(n) в зависимости от m при различных значениях k при = 0, 5, Ru = 3.

Для исследования чувствительности измеряемых величин к положению и размерам области непроницаемости, проводились расчеты значений u(n) в m точках окружности радиуса Ru = 3. При различных формах эллипса (меняется эксцентриситет ), вычислялась величина q() = umin /umax, где umin, umax - минимальное и максимальное значения u(n) в точках измерения соответственно.

На рисунке 2.3 приведены графики зависимости q от при различных положениях источника: 1 = 0, 2 = /2, 3 = 3/4. Источники и приемники расположены на окружности радиуса Ru = 3. Результаты расчетов при = 1, т.е. когда область непроницаемости принимает форму круга, для различных положений источников совпадают.

Из приведенных рисунков видна чувствительность измеряемых величин к положению и размерам области непроницаемости.

Для случая, когда неоднородность представляет собой круг радиуса a, решение прямой задачи в полярных координатах r = r(), [0, 2] с центром в круге, выписывается аналитически. Оно может быть представлено в виде:

где r1 = r0, 1 = 0 - точка положения сопряженного источника, {r0, 0 = 0}- точка положения источника. Функция (2.50) является решением уравнения Лапласа и стремится к нулю на бесконечности. Проверим выполнение граничного условия (2.23):

Далее:

Таким образом, доказано, что граничное условие (2.23) задачи выполняется, т.е.(2.50) является решением задачи(2.22) для круга.

Для проверки точности решения прямой задачи, проводилось сравнение численных расчетов для области непроницаемости, имеющей форму круга используя формулу (2.50) и решение, полученное с использованием интегрального уравнения (2.25) и формулы (2.24). Проведенное исследование показало, что давление, вычисленное с использованием аналитического решения задачи и полученное с использованием интегрального уравнением, совпадают с достаточно большой точностью (погрешность вычислений составляет менее 0, 0001%).

Для численной реализации итерационного метода вводится сетка i, 0 2 и сеточные аналоги функций rj (i ), i, j N и замена интегралов входящих в интегральные уравнения и формулы на квадратурные формулы.

Уравнение (2.37) заменяется конечноразностным, т.е. системой линейных алгебраических уравнений следующего вида:

где N( i, j) = N (i/n, j/m), к системе (2.51) применяется метод регуляризации.

Приведем результаты вычислительного эксперимента. Модельные расчеты проводились для случая, когда контуром C является развернутый на угол относительно центра координат, эллипс:

где a, b - полуоси эллипса, а = b/a-эксцентриситет эллипса. При решении модельной задачи источник располагался в скважине с координатами (0.0, 2.0). При этом считалось, что 7 точек измерения находятся на окружности радиуса Ru, т.е. Rm = Ru = 2 для любого m и шагом /4.

На рисунке 2.4 сплошной чертой изображен контур эллипса с полуосями a = 1.5, b = 0.8, а пунктиром изображен восстановленный контур по 25 точкам. В качестве начального приближения бралась окружность с радиусом 1.

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

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

3.1 Прямая и обратная задача для уравнения параболического типа Рассмотрим следующую задачу для уравнения теплопроводности:

где D – односвязная область c ляпуновской границей, n – единичный вектор нормали к, направленный внутрь D, (xk, yk ) R2 \D – фиксированная точка, P0 – фиксированное число.

Рассмотрим конечный набор точек (xl, yl ) R2 \D, l = 1, n и обозначим через P (t, xl, yl, xk, yk ) решение задачи (3.1) в этих точках.

Обратная задача заключается в определении контура по известным при t > 0 функциям fk (t) и P (t, xl, yl, xk, yk ), где k K 1, 2,..., n, l = 1, n, l = k.

Предположим, что давление в нефтяном пласте не изменяется поперек пласта, а проницаемость в нем постоянна, за исключением области D полной непроницаемости, размеры которой малы по сравнению с размерами пласта. В точках, имеющих координаты (x1, y1 ),...,(xn, yn ), расположены n скважин, при этом в скважине (xk, yk ), 1 k n, создается давление, изменяющееся по закону fk (t). До начала работы скважины (xk, yk ), при всех t > 0 на большом расстоянии от исследуемой области давление в пласте равно некоторому постоянному пластовому давлению P0. В этом случае давление P (t, x, y, xk, yk ) в точке (x, y) в момент времени t вне зоны непроницаемости является решением задачи (3.1).

Рассмотрим задачу (3.1) при произвольном фиксированном k K.

Пусть M и M0 - точки с координатами (x, y) и (xk, yk ) соответственно.

Введем функцию P (t, M, M0 ) = P (t, M, M0 ) P0, тогда (3.1) можем записать как:

Применим преобразование Лапласа:

Для u(p, M, M0 ) получим следующую задачу:

Отметим, что при p > 0 из принципа максимума следует единственность решения задачи (3.3). При p = s2 > 0 фундаментальным решением оператора s2 в R2, стремящимся к нулю на бесконечности, является функция где K0 (z) – функция Макдональда нулевого порядка. При p = s2 > решение (3.3) представим в виде где а функция v(s, M, M0 ) является решением задачи В случае, когда D является кругом радиуса a, задачу (3.4) можно решить аналитически. Данное решение использовано для тестирования поставленной задачи.

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

Будем искать решение задачи (3.4) в виде потенциала простого слоя с непрерывной плотностью µ(N, M0 ), удовлетворяющей уравнению при (N ):

Перейдем в уравнениях (3.5) и (3.6) к полярным координатам с центром в точке O. Будем считать, что начало координат O находится в области D с границей, a полярные координаты точек контура имеют вид (h(), ), [0; 2], h() C 1 [0; 2].

Обратная задача сводится, таким образом, к восстановлению неизвестной функции h() по известным функциям v(s, M, M0 ) и F (s2 ), заданным в одной и той же полярной системе координат с центром в точке Введем обозначения для полярных координат точек M, M0, N и Q:

(r, ), (r0, 0 ), (h(), ) и (h(), ), соответственно. Будем предполагать, что область D имеет такую форму и начало координат расположено в ней таким образом, что выполняется следующее условие:

Указанное условие позволяет переписать подынтегральные выражения в уравнениях (3.5) и (3.6) в отношении входящей неизвестной функции h(), так как в этом случае длина элемента границы может быть записана в виде:

а единичный вектор внутренней нормали к примет вид Таким образом, для вычисления неизвестной функции h() мы получили нелинейное операторное уравнение:

где а функция µ является решением уравнения где 3.3 Численный метод определения границы Вычисление неизвестной границы неоднородности заданной в полярных координатах h() состоит из двух этапов:

-определение начального приближения h0 () искомого контура;

-построение итерационной процедуры для отыскания неизвестной функции h() с заданной точностью.

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

Для решения уравнения (3.10) построим следующий итерационный процесс. Пусть hj () - функция, полученная на j-ой итерации. В качестве начального приближения h0 () неизвестного контура берем найденную ранее окружность. На каждом очередном шаге вместо исходного уравнения (3.10) будет решаться его линеаризация в окрестности функции hj1 (), полученной на предыдущем шаге.

Таким образом, получаем линейное операторное уравнение для функции j (), где j () = hj () hj1 (). Решив это уравнение и вычислив функцию j (), найдем Выпишем оператор B[hj ] и функцию j из уравнения (3.13).

где приближенно задает в точке hj1 дифференциал отображения hj µj1, определяемого уравнениями W (, ) = и удовлетворяет уравнению:

где h3 () + h3 () h3 () cos( ) 3hj1 ()h2 () cos( )+ Таким образом, для вычисления неизвестной границы h() области D построен следующий итерационный процесс:

1. задается начальное приближение h0 ();

2. из уравнения (3.18) для j = 1, 2,... вычисляем µj1 ();

3. из уравнения (3.21) для j = 1, 2,... вычисляем µj1 ();

4. из уравнения (3.17) для j = 1, 2,... вычисляем j1 () и, следовательно, hj () = hj1 () + j1 ().

3.4 Результаты численного моделирования В случае, когда D является кругом радиуса a, задачу (3.4) можно решить аналитически. Выписанное ниже решение будем использовать для тестирования поставленных задач.

Пусть начало координат совпадает с центром круга D, а (r, ) и (r0, 0 ) полярные координаты точек M и M0, соответственно. Тогда система (3.4) примет вид:

где Kl (z) и Il (z) функции Макдональда и Инфельда порядка l, соответственно.

Для проверки точности решения прямой задачи с помощью интегрального уравнения (3.5) и формулы (3.6) проводилось сравнение результатов численных расчетов по этим формулам для области D, имеющей форму круга, с точным решением для круга (3.6). Проведенное исследование показало, что значения функции v(s, M, M0 ), вычисленные с использованием аналитического решения (3.22) и с использованием интегрального уравнения, совпадают с достаточно высокой точностью:

значения отличаются менее, чем на 0,01%.

Численное исследование прямой задачи проводилось для случая, когда контур является эллипсом Точки измерения и возбуждения находились на окружностим радиуса R = 4. На рисунке 3.1 приведены графики решений для четырех вариантов расположения точки возбуждения: (r0, 0 ) равно (4, 0), (4, /6), (4, /2) и (4, ).

Для исследования зависимости u(s2, x, y, x0, y0 ) от параметра s проводилось вычисление этой функции при различных s. На рисунке 3.2 приведены результаты расчетов. Точки измерения располагались на окружности радиуса R = 4 с шагом /15, точка возбуждения имела координаты (4, 7/5).

На рисунке 3.3 приведен результат решения обратной задачи для случая, когда контур является эллипсом с полуосями 3 и 1.5, с центром в точке (1.0, 0) и углом между большей осью и осью абсцисс в /4 (изображен сплошной линией). Пунктирной линией изображено восстановленное расположение контура.

Исходными данными являлись f (t) = 1, p = s2 = 0, 25, (x0, y0 ) = (4, 0), (xj, yj ), j = 1, 7, расположенные на окружности радиуса 4 с шагом /4, начиная с точки (2 2, 2 2). В качастве значений v(s, xi, yi, x0, y0 ) были взяты результаты решения соответствующей прямой задачи с внесенной погрешностью в 2,5%.

На рисунке 3.4 приведен результат решения обратной задачи для случая, когда контур задается функцией h() = 6 + 0, 4 cos 3 + 0, 4 sin (изображен сплошной линией).

Пунктиром изображен восстановленный контур. Точки (xl, yl ), l = 1, 12 были расположены с шагом /6, начиная с точки (8.0, 0) на окружности 8 с центром в начале координат, множество K равнялось {1, 8}, f1 (t) = f8 (t) = 1, p1 = s2 = 0, 1, p2 = s2 = 0, 2 В качестве значений v(s, xl, yl, x0, y0 ) были взяты результаты численного решения соответствующей прямой задачи. Начальным приближением h0 была выбрана окружность с центром в начале координат радиуса 4.

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

2. Предложен и численно исследован метод определения структурных неоднородностей по измерениям акустического поля в ограниченной области в трехмерном пространстве.

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

Список литературы [1] Алифанов О.М. Обратные задачи теплообмена. М.: Машиностроение. 1988. C.244.

[2] Баев А.В., Солтан И.Е. Решение обратной задачи для уравнения колебаний с направленными источниками // Вест.Моск.Ун-та. Сер.

15. Вычисл. Матем.и киберн. 1996. № 1. C. 55-61.

[3] Баев А.В. Математическое моделирование рассеяния акустических волн в трещиноватых средах // Журнал выч. мат. и матем. физики. Т.52. №9. 2012. C.1676-1693.

[4] Баев А.В. Математическое моделирование рефракции акустической волны в окрестности каустики // Журнал выч. мат. и матем.

физики. Т.53. №7. 2013. C.1124-1138.

[5] Бакушинский А.Б., Гончарский А.В. Итеративные методы решения некорректных задач. М.: Наука. 1989. C. 130.

[6] Бакушинский А.Б., Гончарский А.В. Некорректные задачи. Численные методы и приложения. М.: Изд-во Моск. ун-та. 1989. С.

[7] Бакушинский А.Б., Левитан С.Ю. Некоторые модели и численные методы нелинейной вычислительной диагностики //М. сб. трудов ВНИИСИ АН СССР. 1991. 13. C. 3-25.

[8] Бакушинский А.Б., Поляк Б.Т. О решении вариационных неравенств // ДАН СССР. 1974. Т. 219. №5. C. 1038-1041.

[9] Бахвалов Н.С. Численные методы. М.: Наука. 1973. С. 631.

[10] Березин Ф.А., Шубин М.А. Уравнение Шредингера.М.: Изд-во Моск.ун-та. 1983. C. 345.

[11] Бицадзе А.В. Уравнения математической физики. М.: Наука. 1982.

[12] Воеводин В.В. Вычислительные основы линейной алгебры. М.: Наука. 1967. С. 500.

[13] Владимиров В.С. Уравнения математической физики. М.:Наука.

1981. С. 512.

[14] Гахов Ф.Д. Краевые задачи. М.: Наука. 1977. C. 640.

[15] Гельфанд И. М., Левитан Б. М. Об определении дифференциального уравнения по его спектральной функции // Изв. АН СССР.

Сер. матем. 15:4. 1951. С. 309–360.

[16] Гласко В.Б. Некоторые обратные задачи сейсмологии // Некорректные задачи естествознания. М.: Изд-во Моск. Ун-та. 1987. C.

[17] Годунов С.К. Уравнения математической физике. М.: Наука.

1971.C. 416.

[18] Горюнов А.А., Сасковец А.В. Обратные задачи рассеяния в акустике. М., Изд. МГУ. 1989. C.152.

[19] Гохберг И. Ц., Крейн М. Г. Теория вольтерровых операторов в гильбертовом пространстве и ее приложения. М.: Наука. 1967. С.508.

[20] Довгий С.А., Лифанов И.К. Методы решения интегральных уравнений. Киев.: Наукова думка. 2002.C.344.

[21] Денисов А.М. Введение в теорию решения обратных задач: Учеб.

пособие. М.: Изд-во МГУ. 1994. C. 208.

[22] Ершов А.А. Асимптотика решения уравнения Лапласа со смешанными условиями на границе// Ж. вычисл. мат. и мат.физ. 2011. 51.

№6. C.1064-1080.

[23] Жданов М.С. Теория обратных задач и регуляризации в геофизике. М.: Научный мир. 2007. C. 712.

[24] Кангужин Б.Е., Анияров А.А. Математические заметки. 2011. 89.

№6. C.856-867.

[25] Клаербоут Д.Ф. Сейсмическое изображение земных недр. М.: Недра. 1989. С.408.

[26] Колтон Д., Кресс Р. Методы интегральных уравнений в теории рассеяния. М.: МИР. 1987. С.312.

[27] Ладыженская О.А. Краевые задачи математической физики. М.:

Наука. 1973. С.408.

[28] Лебедев Н.Н. Специальные функции и их применения. М.: ФИЗМАТЛИТ. 1963. C.360.

[29] Левитан Б. М. Об асимптотическом поведении функции Грина и разложении по собственным функциям уравнения Шредингера // Матем. сб., 41(83):4. 1957. С. 439–458.

[30] Никифоров А.Ф., Уваров В.Б. Специальные функции математической физики. М.: Наука. 1984.C.344.

[31] Прямые и обратные задачи математической физики / Сборник.Под. Ред. Тихонова А.Н., Самарского А.А. М.: Изд-во Моск.

Ун-та. 1991. C. 258.

[32] Рамм А.Г. Многомерные обратные задачи рассеяния. М.: Мир.

1994.C.495.

[33] Рихтмайер Р.Принципы современной математической физики. М.:

Мир. 1982.C.488.

[34] Cтовас А.М. Борновское приближение в задачах миграции и инверсии для акустического уравнения // Физика Земли. 1992. №3.

C. 29-37.

[35] Тихонов А.Н., Самарский А.А. Уравнения математической физики. М.:Наука. 1972. С. 735.

[36] Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач.

М.: Наука. 1986. C. 288.

[37] Тихонов А.Н., Гончарский А.В., Степанов В.В., Ягола А.Г. Численные методы решения некорректно поставленных задач. М.: Наука.

1990. C. 264.

[38] Файзуллин И.С., Шапиро С.А. О затухании упругих волн в горных породах, звязанном с рассеянием на дискретных неоднородностях.

ДАН СССР. 1987. 295. 2. C. 341-344.

[39] Файзуллин И.С., Шапиро С.А. Особенности затухания сейсмических волн в случайно-неоднородных средах. ДАН СССР. 1988. 302.

№5. C. 1073-1077.

[40] Херен Г. Восстановление изображений по проекциям. Основы реконструктивной томографии. М.: Мир. 1983.

[41] Ц.На. Вычислительные методы решения прикладных граничных задач. М.: Мир. 1982.C.296.

[42] Berkhout A.J. Seismic inversion in terms pre-stack migration and multiple elimination // Proc.IEEE 1986.V.74.№3. P. 415-427.

[43] Blok H., Oristaglio M. Waveeld imaging and inversion in electromagnetics and acoustics // Delft University of Technology.

Report number. Et/EM 1995-21. 1995. Р.132.

[44] Сhadan K., Sabatier P.C. Inverse Problems in Quantum Scattering Theory. New York. Springer-Verleg. 1989. C. 499.

[45] Chandler-Wilde Simon N., Heinemeyer Eric, Potthast Roland.Acoustic Scattering by mildly raught unbounded Surfaces in three dimensions //SIAM. J.Appl.Math. 2006. 66. №3. P.1002-1026.

[46] Kress R. A factorization method for an inverse Neumann problem for harmonic vector elds//Georg.Math.J. 2003. 10. №3. P.549-560.

[47] Liu Xiaodong, Zhang B. Direct and inverse obstacle scattering problems in a piecewise homogeneous medium //SIAM. J.Appl.Math. 2010. 70.

№7-8. C.3105-3120.

[48] Masaru Ikehata and Mishio Kawashita. The enclosure method for the heat equation//Inverse Problems. 2009. V.25. P.1-10.

[49] Ramm A.G. An inversion problem for Helmgoltz’s equation // Appl.Math.Lett. 89. 2. №1. Р. 105-108.

[50] Ramm A.G. Uniqueness result for inverse problem of geophysics // Inverse problems. 1990. 6. P. 635-642.

[51] Ramm A.G., Sjostrand J.S. An inverse problem for the wave equation //Math. Zeitsch. 1991. 206. P. 119-130.

[52] Ramm A.G. Non-uniqueness of the solution to an inverse problem in geophysics // Inverse Problems. 1986. v.2. no.2. P. 123-125.

[53] Ramm A.G. Completness of the products of solutions to PDE and uniqueness theorems in inverse scattering // Inverse Problems. 1987.

v.3. no.4. P. L77-L82.

[54] Roman Chapko, Rainer Kress and Jeong-Rock Yoon. An inverse boundary value problem for the heat equation:the Neumann condition //Inverse Problems. 1999. V.15. P.1033-1046.

inversion//Inverse problems. 1988. 4. L. P. 7-10.

[56] Sergio Vessella. Quantitative estimates of unique continuation for parabolic equations,determination of unknown time-varying boundaries and optimal stability estimates// Inverse Problems. 2008.

V.24. 023001(topical review) [57] Tarantola A. Inverse Problem Theory. Amsterdam. Oxford. New York.

Tokyo: Elsevier. 1987. P. 613.

[58] Tikhonov A.N., Goncharsky A.V., Stepanov V.V., Yagola A.G.

Numerical Methods for the Solution of ill-posed Problems. Kluwer.

Dordrecht. 1995. P.324.

[59] Головина С.Г., Романов С.Ю., Степанов В.В. Об одной обратной задаче сейсмики// Вест.Моск.Ун-та. Сер. 15. Вычисл. Матем.и киберн. 1994. № 4. C. 16-21.

[60] Головина С.Г., Никитина Е.В. Численный анализ методов определения спектральной амплитуды акустического поля. Вест.Моск.Унта. Сер. 15. Вычисл.Матем.и киберн. 1997. № 3. C. 20-23.

[61] Головина С.Г., Денисов А.М., Дмитриев В.И. Об обратной задаче определения зон малой проницаемости в нефтяном пласте // Прикладная математика и информатика. М.: МАКС Пресс. 2005. № 21.

[62] Головина С.Г., Денисов А.М., Дмитриев В.И. Определение границ зон малой проницаемости// Материалы международной конференции Тихонов-100. Москва. Изд-во МГУ. № 4. C.31- [63] Головина С.Г. Метод линеаризации в обратной задаче определения зон малой проницаемости в нефтяном пласте// Вест.Моск.Ун-та.

Сер. 15. Вычисл. Матем.и киберн. 2008. № 1. C. 5-9.

[64] Головина С.Г., Разборов А.Г. О восстановлении неизвестной границы в смешанной задаче для уравнения теплопроводности // Прикладная математика и информатика. М.: МАКС Пресс. 2008. № 30.

C. 18-24.

[65] Головина С.Г., Разборов А.Г. Об определении границы двумерной области по решению внешней начально-краевой задачи для уравнения теплопроводности // Прикладная математика и информатика.

М.: МАКС Пресс. 2009. № 33. C. 69-74.

[66] Головина С.Г., Разборов А.Г. Об определении зон малой проницаемости в нефтяном пласте по давлению в скважинах // Проблемы динамического управления. Серия:Сборник научных трудов факультета ВМиК МГУ им.Ломоносова. М.: МАКС Пресс. 2010. Том 5. C. 42-48.

[67] Головина С.Г. Определение границы локальной неоднородности по измерениям акустического поля//Ломоносовские чтения - 2013. Москва. МГУ. ВМК. URL:http://www.cs.msu.su/ sites/cmc/les/docs/programma _i_tezisy_2013_3.pdf. 2013. С.32Golovina S.G., Denisov А.М., Dmitriev V.I. // Book of abstracts of international conference.Moscow State Lomonosov University Moscow.

Rassia. P.56-57.

[69] Golovina S.G., Razborov A.G. Finding the unknown boundary in the initial boundary-value problem for the heat equation// Computation Mathematics and Modeling (Springer New York) 2009. V.20. № 3. с.

231-236.

[70] Golovina S.G., Razborov A.G. Determining the boundary of a twodimensional region from the solution of the external initial boundaryvalue problem for the heat equation // Computation Mathematics and Modeling (Springer New York) 2010. V.21. № 2. C. 178-183.