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

3 МДР 1ЯЯ7

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

ЗЕРКАЛЬ Сергей Михайлович

РЕШЕНИЕ СЕЙСМОТОМОГРАФИЧЕСКИХ ЗАДАЧ МЕТОДАМИ ИТЕРАТИВНОЙ И ДЕСКРИПТИВНОЙ РЕГУЛЯРИЗАЦИИ

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

АВТОРЕФЕРАТ диссертации на соискание ученой степени доктора технических наук

Новосибирск - 1997

Работа выполнена в Институте математики им. С. Л. Соболева СО РАН

Научный консультант

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

профессор Воск'обойников Ю.Е.

Официальные оппоненты: •

Защита состоится "19" марта 1997 года в 10.00 часов на заседании диссертационного совета Д 063.34.03. в Новосибирском государственном техническом университете по адресу: 630087 Новосибирск, пр. Маркса, 20.

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

доктор технических наук профессор Спекгор А.А.

доктор технических наук профессор Кузнецов В.В.

доктор физико-математических наук профессор Трофимов О.Е.

Ведущая организация:

Объединенный институт физики Земли им. О.Ю. Шмидта РАН

Автореферат разослан

Ученый секретарь диссертационного совета к.т.н., доцент

Общая характеристика работы

Актуальность работы.

В последние два десятилетия возникло и интенсивно развивается, проникая в различные области человеческого знания, новое направление в науке - вычислительная (компьютерная) томография. Метод компьютерной томографии позволяет по характеристикам пропущенного (трансмиссионная томография) через исследуемый объект или его собственного (эмиссионная томография) излучения 'судить о внутреннем устройстве, физическом состоянии или химическом составе этого объекта. К настоящему времени в геофизическую терминологию прочно вошел термин "сейсмическая томография". Это обусловлено тем, что постановки ряда основных задач в сейсмике являются, по существу, классическими для томографии.

Типичным примером такой задачи является обратная кинематическая задача (ОКЗ), результат решения которой определяет скорость распространения зондирующего сигнала в исследуемой среде, что весьма важно не только для геофизики, но и, в частности, для оптики, гидроакустики. Возникновение обратной кинематической задачи в различных областях науки и человеческой * деятельности характеризует фундаментальность и большое прикладное значение данной задачи. Впервые обратная кинематическая задача была решена в предположении сферической симметрии Земли Г.Герглощем и Е.Вихертом (1907 г.), и только в 60-х годах, начиная с работ М.М.Лаврентъева и В.Г.Романова, началось систематическое исследование неодномерной ОКЗ. Результат М.М. Лаврентьева и В.Г. Романова был сформулирован для случая двух переменных и относится к линеаризованной постановке, характер некорректности полученной задачи сильный, как у задачи Коши для уравнения Лапласа. Конструктивное решение ОКЗ в общей "классической" постановке довольно затруднительно и до сих пор не получено.

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

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

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

Исследуемые в диссертации обратные задачи - задачи некорректные, и для их решения необходимо использовать метод регуляризации, основы которого заложены в работах А. Н. Тихонова, М. М. Лаврентьева и других известных специалистов в этой области. Если при решении слабо некорректной ОКЗ с полными данными применяется простая дескриптивная регуляризация, использующая априорную информацию о гладкости решения, векторе измерений, и заключающаяся в сглаживании массива данных и результата решения обратной задачи кубическими сплайнами, то решение ОКЗ с неполными данными и задачи определения правой части волнового уравнения в исследуемой дискретной постановке требует более сложной комбинации итеративной и дескриптивной регуляризации.

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

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

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

в) разработка устойчивых эффективных численных алгоритмов, реализующих решение задач а) и б);

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

Научная новизна диссертационной работы определяется следующими результатами:

- Рассмотрена новая (сейсмотомографическая) постановка обратной кинематической задачи, позволяющая снять ее переопределенность и обеспечивающая слабый характер некорректности задачи.

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

- Разработан и численно исследован алгоритм решения обратной кинематической задачи на основе метода вычислительной томографии.

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

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

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

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

Основными положениями, выносимыми на защиту, являются:

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

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

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

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

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

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

7. Развитие на основе проведенных исследований нового направления геофизике - сейсмической томографии.

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

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

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

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

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

Внедрение результатов работы. Исследования по теме диссертации проводились в основном в соответствии с координационными планами НИР ВЦ СО АН СССР, Новосибирского госуниверситета с 1983 по 1986 гг. и Института математики СО РАН с 1987 по 1995 гг.

Созданное программное обеспечение внедрено в МНПО "Спектр", п/я В-2962, в/ч 69104, что подтверждается соответствующими актами о внедрении в приложении к диссертации.

Апробация работы. Основные результаты докладывались и обсуждались на:

- Всесоюзной школе-семинаре "Теория и методы решения некорректно поставленных задач и их приложения" (Самарканд, 1983);

Всесоюзных симпозиумах по вычислительной томографии (Новосибирск, 1983; Куйбышев, 1985; Киев, 1986; Звенигород (Моск. обл.), 1991);

- Всесоюзном семинаре по оптической томографии (Таллин, 1988);

- Республиканской научно-технической конференции "Сейсмические методы поиска и разведки месторождений полезных ископаем ых"(Киен, 1991);

Республиканской научно-технической конференции "Дифференциальные уравнения и их приложения" (Таджикистан, Куляб, 1991);

- Международной школе-семинаре "Методы оптимизации и их приложения" (Иркутск, 1995).

На семинарах ВЦ СО РАН, ИТПМ СО РАН, ИГФ СО РАН (руководитель чл. - корр., проф. C.B. Гольдин), ИМ СО РАН (руководитель академик РАН, проф. М.М. Лаврентьев) и других научных семинарах.

Публикации. Общее количество печатных работ более 30. Основные результаты диссертации опубликованы в 15 научных статьях, в том числе в центральных научных журналах "Доклады АН СССР", "Геология и геофизика", Известия РАН (серия "Физика Земли").

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

Содержание работы

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

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

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

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

Точная математическая постановка обратной кинематической задачи формулируется следующим образом [1].

Пусть в области и пространства И > 1, ограниченной

поверхностью происходит волновой процесс, порожденный точечным источником колебаний, приложенным в точке Х°, и пусть известна функция Т(х°,х) для всевозможных х° еБ, X еВ . Требуется найти скорость передачи сигналов У(х) для Х° £ £>.

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

Впервые обратная кинематическая задача сейсмики была рассмотрена Г. Герглотцем и Е. Вихертом (в 1905-1907 гг.) в предположении сферической

симметрии Земли. Решение задачи оказалось непрерывным из С1 —> С . При решении практических задач функция Т задана дискретно. Одномерная обратная кинематическая задача является некорректной задачей и характер ее некорректности такой же, как в задаче численного дифференцирования таблично заданной функции. Неоднозначность задачи была изучена М.Л. Гервером и В.М. Маркушевичем (1965-1967 гг.).

До середины 60-х годов все усилия были направлены на исследования одномерной задачи и только начиная с работы [2] М.М.Лаврентьева и В.Г.Романова началось систематическое исследование многомерных постановок обратной кинематической задачи.

Первый результат для неодномерной обратной кинематической задачи был получен М.М.Лаврентьевым, В.Г.Романовым и сформулирован в двумерном пространстве .X; и Хг для области, представляющей. собой

полуплоскость Х2 ^ 0. Этот результат относится к линеаризованной постановке. Предполагается, что скорость исследуемого сигнала раскладывается в следующую сумму:

У(х) = УоШ+Г](х), х б Я2, 8

где Уо(х%) известна, а У](х) мала по сравнению Уо(х). Тогда с точностью до малой более высокого порядка, чем У*(х), справедливо

Го(х0,х')

здесь То - время пробега сигнала вдоль луча Го в среде со скоростью

V- Уо(х).

Было показано, что для случая

У0(х)=АХ2 + В, А > О, В > 0, (3)

У](х) однозначно определяется из (2) по данным Т(х°¡X1), для пар точек X°,Х1 расположенных на прямой Х2 = 0. Установлено, что задача определения V) из уравнения (2) сильно некорректна. Характер ее некорректности такой же, как у задачи Коши для уравнения Лапласа.

Коллективом в составе А.С.Алексеева, М.М.Лаврентьева, Р.Г.Мухометова, Н.Л.Нерсесова и В.Г.Романова в 1969-1971 гг. на основе этого результата успешно была выполнена численная обработка сейсмологических данных и получено распределение скорости продольных волн на интервале большого круга от Памира до Байкала на глубину в среднем до 350 км. Под большим кругом пошшается сечение Земли плоскостью, проходящей через ее центр, на границе этого круга 0 I

располагаются точки X и X .

Дальнейшие теоретические исследования многомерной обратной кинематической задачи проводились Ю.Е.Аниконовым, Р.Г.Мухометовым и В.Г.Романовым, Г.Я.Бейлькиным, И.Н.Бернштейном и М.Л.Гервером.'

Постановка обратной кинематической задачи в диссертации формулируется для трехмерно-неоднородной среды, допускающей линеаризацию аналогично формулам (1) -(3), только теперь X =(х ¡,Х2,2). В этом случае задача определения У](х) по

является

переопределенной: по функции шести переменных ^определяется функция трех переменных У/(х), если источники и приемники сигнала расположить на плоскости 2 — 0, то все равно переопределенность сохранится: по функции четырех переменных определяется функция трех переменных.

В диссертации в качестве системы наблюдений предлагается окружность на плоскости г = О

хх2+Х22 = Г2. (4)

В случае такой системы наблюдений переопределенность снимается, так как теперь ^-функция двух полярных углов ф\ и фг на источник и приемник соответственно и радиуса Г, то есть по функции трех переменных Т(<Р\, фг, 1") восстанавливается функция трех переменных У(х) .

Заметим, что множество лучей Го в случае (Уу—А +Вг), натянутых на окружность (4), образует поверхность шарового сегмента (см. рис. 1).

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

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

Итак, находясь в условиях линеаризованной постановки, воспользуемся следующим трехмерным аналогом уравнения (2):

Тх(х°,х1) = ГО0,*1)- ГоСЛ*1) и (5)

Го(Л*')

здесь Л7| (х) =---——, а Го(х°,Х1) определяется из системы

У{х) у0(х)

уравнений

н2 =Р2 ,/?=oonst,p>0,

(x,v)-p~Q,v= (sin0-cos в),р- const

Геометрический смысл величин р, р и 6 показан на рис. 2.

Перейдем в выражении (5) от интегрирования по Г0(х0,х') к интегрированию вдоль Обозначим:

Будем рассматривать случай 2 £ О, тогда перед корнем следует выбрать знак "+". Для ск имеем (см.рис.2.)

ds = dly¡l + (z¡?.

Уравнение прямой l(x°,xl): JCj sin 9— Х^ COS 9 — р = О, тогда

di _ di . _ z¡ = —cos —sin 0,

di _ ___д:2

ds = di.

1 +

COS в + Х2 SU10

Г~2 2 2 VР ~*1 ~*2 J

_ - sin 6>- eos fl)* _ У/72 - р1

* А

z +в

* А

2 +s

di.

Подставляем полученное для fife выражение в (5) и получаем

фг»*.) = Г д

\1Р2~Р2 'М

в

или

сй.

г +-

В

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

в

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

'Фр^-'Ь!/ +^яА}{р,е), (7)

где Я ' - оператор обратного преобразования Радона.

Определяя искомую функцию по формуле обращения (7) в круге

~> 2^2

.V] + Л'2 ь г , мы тем самым определяем на поверхности шарового сегмента, образованного лучами Гр, опирающимися на окружность (4).

Меняя г, получаем систему вложенных шаровых сегментов, на поверхности которых известна функция /1], и тем самым получаем решение трехмерной ■задач» в объеме, ограниченном поверхностью внешнего сегмента (с наибольшим Г) и поверхностью внутреннего (с наименьшим г).

Обозначим х = и рассмотрим следующее уравнение:

\с^х\р2)5{(х,у)-р)ф,(1^к\г))ск = /{р,9) , х еЯ2,

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

(параметр/" фиксирован) допускал факторизацию вида

о(|х|,р) = а1(|л|),о2(/72)- (8)

Теорема. Вес допускает факторизацию вида (8) тогда и только

тогда, когда

п0(г) = (А + Вк)~1.

Отметим, что при выводе формулы (7)

„ а2{р2) = 4р2-р2.

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

угол в] =0-1)Ав, У = 1..... М , Л9= 180°/М, а М * число

проекций. Матрица размером N X М , столбцами которой являются проекции (Л^ - число пар источник - приемник), используется в качестве входных данных для блока решения обратной задачи и называется проекционной матрицей.

В нашем случае реконструкция искомой функции производилась по алгоритму Ерохина - Шнейдерова. На первом шаге регуляризация задачи осуществлялась применением статистического сглаживания данных прямой задачи кубическими сплайнами, учитывающими уровень погрешности исходных данных [3]. Исследования проводились на средах двух типов:

квадратичная добавка - П{(х) = СОоХ2 , добавка в виде гауссиана -

-(со,(х0-х))2 _3

П](х) = ,(Од =С0П81, СО - постоянный вектор, X &К .

Во всех вариантах щ(хз) = (4+ 0.5 Хз)" .

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

Рассмотрим первую среду. Пусть для нее выполняются следующие условия:

пшЬ:,] = тн^л^] = 1, 0 <х, <1, сц =ссгб1

Очевидно, что при этих условиях размер добавки зависит от й)о.

Качество восстановления с увеличением П\ ухудшается. Начиная с й)о = 0.15 , график восстанавливаемой функции заметно отклоняется от исходного, сохраняя при этом достаточно информации о поведении искомой

функции. При СОо = 0.2 ближе к краям, где П] достигает своего наибольшего значения в исследуемой области, на графике появляются изгибы, переходящие при СОо ~ 0.255 в осцилляции. В данном случае относительный размер добавки следует ограничить 20 % от По.

При наложении случайной помехи на значения 7/ = Т - То уровень шума задавался до 10%. Несмотря на то, что уровень помехи выбирался высоким для данной некорректной задачи, функция П] согласно разработанному алгоритму восстанавливается достаточно устойчиво. Следует отметить, что кроме случайного шума присутствуют погрешности систематического характера, вносимые в Т при решении прямой задачи:

ошибка метода Рунге — Кутта при решении системы дифференциальных уравнений;

-ошибка пристрелки;

- ошибка интерполяции;

- приближенный характер формулы.

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

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

диктуется как экономико-материальными требованиями, так и трудностями чисто технического характера. Следует отметить, что вклад источников и проекций (направлений просвечивания) в качество получаемого решения обратной задачи различен и экономия одного источника не эквивалентна экономии одной проекции даже в случае N = М. Более того, для конкретных типов сред существуют для чисел N и М нижние и верхние границы У, л/ и n\ м" •, то есть такие значения, что для N < Л^ и М < л/ происходит полная потеря информации, а в случае N > N и м > м* существенного улучшения качества решения обратной задачи не происходит.

Рассмотрим вторую среду. Зафиксируем N = 15 и будем изменять М. Эксперименты проводились для ¿Уо=0.1, £У = (-JÍ5,-ЙЗ, 1), Хо — (0,0,-1). Значительное отличие установлено только в случае М = 3 и М = 5 и обусловливается присутствием артефактов. Артефакты - это волнистая структура у подножия восстановленного гауссиана. На рис. 3, А артефакты показаны стрелками, количество их равно 7 (А/ = 7). Количество артефактов находится в прямой зависимости от количества проекций. Артефакты следует отнести к собственным шумам алгоритма. С увеличением числа направлений просвечивания количество артефактов возрастает, а амплитуда их уменьшается.

Если проанализировать восстановление при одном числе проекций М = 9, но при существенно разном количестве источников N - 5 ч N = 15, то можно отметить, что в случае пяти источников высота гауссиана составляет.*) 33 , что это « 60 % от истинного значения -.0 55, кроме того, восстановленный гауссиан значительно шире у подножия, чем исходный. Иная картина в случае пятнадцати источников, здесь высота гауссиана составляет.0 58, что ~ 105 % от истинного значения, и восстановленный гауссиан достаточно хорошо аппроксимирует исходный.

При исследовании на ЭВМ алгоритма решения обратной кинематической задачи проведен следующий анализ: а) исследовалась зависимость качества восстановления от величины добавки относительно линейного сигнала; б) с помощью генератора псевдослучайных чисел вводился шум в исходные данные, то есть в значения разности T¡ — Т - То, пропорцианально их локальному значению; в) варьировалась система наблюдений - изменялось количество источников и число проекций.

По этим пунктам можно сделать следующие выводы. Верхняя граница размера добавки n¡ — п - щ, когда формула (7) дает удовлетворительные

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

При исследовании алгоритма решения обратной задачи рассматривались функции п¡, симметричные относительно координатной оси Хз. Это объясняется большим объемом требовавшихся вычислений на ЭВМ, а в исследованных случаях достаточно было насчитать только половину первой проекции, из которой потом генерировались остальные. Ниже определяется функция П/ в случае, когда ее ось симметрии не проходит через центр системы координат. В такой ситуации требуется отдельное вычисление каждой проекции, что естественно увеличивает время решения прямой задачи. Рассматриваемый случай моделирует определение локального возмущения линейного показателя преломления в произвольной части исследуемой области трехмерного пространства.

При проведении численного эксперимента в качестве П/ выбиралась следующая функция:

1*00 = оле-15*?"15^-0-3)2-^'-1)2.

На рис. 3, А-В, приведены результаты определения функции П] в точках поверхности, образованной лучами, соединяющими источники и приемники для системы наблюдений из 15 источников и соответственно 15 приемников, число проекций в данном случае равнялось 7. Полная погрешность измерения Т/ — Т - То оценивается в 3 % от локального значения.

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

изображение исходной и восстановленной функции П\ соответственно. На рис. 3, В, показано совместное изображение сечений исходной и восстановленной функции nj плоскостью Xj — 0 .

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

Обе рассматриваемые в диссертации обратные задачи - задачи некорректные, и для их решения необходимо использовать ' метод регуляризации, основы которого заложены в работах А.Н.Тихонова, М.М.Лаврентьева и ряда других известных специалистов в этой области. Если при решении слабо некорректной задачи главы II применяется простая дескриптивная регуляризация, использующая априорную информацию о гладкости решения, векторе измерений и заключающаяся в сглаживании массива данных и результата решения обратной задачи кубическими сплайнами, то решение задач глав III и IV требует более сложной комбинации итеративной и дескриптивной регуляризации.

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

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

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

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

17

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

В §3.1 излагаются постановка обратной кинематической задачи с неполными данными, обусловленными локальной непрозрачностью в исследуемой области, и вывод основной системы линейных алгебраических уравнений (СЛАУ), получаемой в результате алгебраизации интегрального преобразования Радона.

В § 3.2 - 3.3 приводится используемый в вычислениях регуляризующий итерационный алгоритм последовательного учета строк матрицы основной СЛАУ, дается описание нелинейного комбинированного алгоритма, позволяющего преодолевать трудности, возникающие в результате итерационного решения данной задачи, а именно подавлять как собственные шумы алгоритма, так и случайные искажения в канале приема сигнала. Также формулируются рекомендации по планированию эксперимента и определению оптимальных параметров наблюдений, гарантирующих достоверность получаемого решения обратной задачи. Приводится численный эксперимент, моделирующий решение поставленной задачи.

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

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

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

считаются обратная задача определения коэффициента при операторе Лапласа и обратная задача определения правой части.

Задачи, решению которых посвящена данная диссертация, объединяются как обратные задачи определения коэффициента V(x) или правой части f(x, t) волнового уравнения

4т- V2{x)bu = f{x,t) , (9)

дГ

Л - оператор Лапласа, по амплитудно-фазовой информации о его решении, заданной в конечном числе точек.

Задача определения коэффициента V(x), - имеющего физический смысл скорости распространения волнового сигнала в рассматриваемой среде, решается в кинематической постановке. Рассмотрим следующее уравнение:

Au+k2n2(x,z)u= S(x -x°)S(z)t

х eR2, Z eRl , n = v~\ {Щ

являющееся стационарным следствием из (9). При достаточно большой частоте к решение уравнения Щ имеет, как известно, следующий вид:

и = А(х, z, х°) ехр z, j, где А удовлетворяет уравнению переноса, а фаза X -уравнению эйконала:

[grad т]2 =n2(x,z).

Упрощающее предположение относительно частоты к носит асимптотический характер. Уравнение (5) может быть получено в результате линеаризации уравнения эйконала. В случае определения правой части, в качестве упрощающего предположения, имеющего также асимптотический характер, мы считаем, что расстояние между источниками и приемниками сигнала г » 1, и фактически исследуем задачу при г —> оо, что позволяет свести ее к решению системы трансцендентных уравнений специального вида. Подлежащая определению функция f(x,t) имеет вид линейной комбинации точечных источников:

т .,

/<Х0= Z AjS(x-Xj)e (Ц)

7=1

здесь т - число источников, Aj - мощность, а .Xj - местоположение у-го

источника, к = 2я7 Я, Л - длина волны.

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

Задачи локации подразделяются на два класса - пассивную, и активную. Пассивная локация состоит в определении характеристик излучателя (дискретного либо непрерывного) по его полю, замеренному на удаленной от излучателя антенне - непрерывной или дискретной. Активная локация заключается в определении характеристик неоднородной среды по зондирующему сигналу, замеренному на некоторой антенне. Источники, создающие первичное (облучающее) поле, считаются известными.

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

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

криволинейный вид . Повысить устойчивость ее численного решения можно следующими путями:

1. Учетом всей априорной информации о решении.

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

3. Увеличением числа наблюдений.

Развитый в диссертации алгоритм определения /(.х,1') вида 01) в сложной, неустойчивой помехо-сигнальной ситуации применим как к прямолинейным, так и к криволинейным антеннам и не требует эквидистантности.

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

В §4.1 приводятся математическая постановка задачи и вывод основной системы нелинейных трансцендентных уравнений. Задача ставится как задача определения правой части специального вида волнового уравнения по известным функционалам от его решения в конечном числе точек дуги окружности (антенне). Скорость распространения сигнала в среде считается постоянной.

В §4.2 излагается алгоритм для представляющей практический интерес криволинейной антенны. В основу алгоритма положен метод наименьших квадратов. Для решения полученной системы уравнений применяются различные регуляризованные методы (метод Гаусса-Ньютона, метод проекции градиента). Алгоритм использует информацию после предварительной статистической обработки.

В §4.3 показывается, что алгоритм допускает редукцию к прямолинейной антенне и не требует эквидистантности.

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

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

Основные результаты и выводы

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

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

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

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

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

6. Проведено численное исследование решения задачи диагностики (разрешения) источников волнового поля различными методами, использующими итеративную и дескриптивную регуляризацию.

7. Разработанный алгоритм решения обратной кинематической задачи работает устойчиво к случайным искажениям входных данных, граница применимости линеаризации, определяемая относительным размером нелинейной трехмерной части показателя преломления к его основной линейной составляющей, может доходить до 20 %. Параметры системы наблюдений оказывают активное воздействие на качество решения обратной задачи, однако можно сделать вывод, что в исследованных моделях для получения качественной информации о решении обратной задачи может быть достаточно 3^5 проекций с 5 -г 7 элементами.

8. При решении задачи диагностики источников волнового поля установлено, что наилучший эффект получен в результате применения метода проекции градиента. Регуляризованные модификации метода Гаусса -Ньютона оказались малопригодными. Ранее при решении подобной задачи угловое расстояние между источниками составляло не менее 9° [4], применение же разработанного алгоритма, имеющего за счет использования регуляризации повышенную разрешающую способность, позволило уверенно исследовать ситуации с угловыми расстояниями до 1,5° при отличии амплитуд источников на порядок.

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

10. Рассмотренная в диссертации постановка обратной кинематической задачи относится к нелинейной трансмиссионной топографии, задача диагностики источников излучения - к эмиссионной томографии, а в целом эти исследования развивают новое направление в геофизике — сейсмическую томографию.

Литератора

1. Лаврентьев М.М., Романов В.Г., Шишатский С.П. Некорректные задачи математической физики и анализа. М.: Наука, 1980. 286 с.

2. Лаврентьев М.М., Романов В.Г. О трех линеаризованных обратных задачах для гиперболических уравнений // Докл. АН СССР. 1966. Т 171, № 6. С. 1279-1281.

3. Пикапов В. В. Комплекс программ томографической реконструкции локальных характеристик плазмы // Комплексы программ математической физики (СКПФМ - VII). Новосибирск, 1982. С. 132 - 135.

4. Алексеев В.И., Гительсон B.C., Глебов Г.М., Каленов E.H., Тихонравов В.Н. Точность определения параметров источников случайных акустических сигналов методом прямого разрешения // Акустический журнал. 1981. Т. XXVII, вып.1. С. 30-35.

Основные результаты диссертации опубликованы в следующих работах :

5. Зеркаль С.М. Трехмерная обратная кинематическая задача сейсмики в линеаризованной постановке и ее численное исследование // Теория и методы решения некорректно поставленных задач и их приложение. Новосибирск, 1983. С. 104-105.

6. Зеркаль С.М.. Хогоев Е.А. Результаты численного решения обратной трехмерной линеаризованной задачи восстановления коэффициента

преломления. Новосибирск, 1984. 26 с.(Препринт/АНСССР. Сиб. Огд-ин. Вычислительный центр; №524).

7.3еркаль С.М. Численное решение системы нелинейных уравнений обратной задачи распространения волн // Обратные задачи математической физики. Новосибирск, 1985. С. 76-89.

8. Зеркаль С.М. Регуляризованный алгоритм решения обратной задачи распространения волн // II Всесоюзный симпозиум по вычислительной томографии, Куйбышев, 1985. С. 60-61.

9. Зеркаль С.М.. Хогоев Е.А. Обратная задача определения малых локальных неоднородностей квазилинейного показателя преломления // Линейные и нелинейные задачи вычислительной томографии. Новосибирск, 1985. С. 61-66.

10. Зеркаль С.М. Малоракурсная нелинейная томография // Реконструктивная томография. Куйбышев, 1987. С. 22-27.

11. Зеркаль С.М.. Зеркаль Л.Б.. Брагин М.Ю. Вычисление проекций в нелинейной томографии // Реконструктивная томография. Куйбышев, 1987. С. 28-33.

12. Зеркаль С.М. Численное решение обратной трехмерной кинематической задачи сейсмики в линеаризованной постановке // Геология и геофизика. 1988. № 11. С.126-133.

13. Зеркаль С.М.. Хахлютнн В.П, Об определении функции по ее интегралам вдоль некоторого семейства лучей // Оптическая томография. Таллин, 1988. С. 102-104.

14. Зеркаль С.М. Кинематическая томография // V Всесоюзный симпозиум по вычислительной томографии. М., 1991. С. 101.

15. Зеркаль С.М. Определение непрозрачных зон в Земле методом компьютерной томографии в кинематической постановке // Докл. АН СССР. 1991. Т. 317, №2. С. 330-333.

16. Гарин В.П.. Зеркаль С.М.. Хогоев Е.А. К решению обратной задачи сейсморазведки методом компьютерной томографии // Сейсмические методы поиска и разведки месторождений полезных ископаемых. Киев, 1991. С. 24.

17. Бронников A.B.. Воскобойников Ю.Е.. Зеркаль С.М. Решение обратной кинематической задачи методом компьютерной томографии для сред, содержащих непрозрачные включения и допускающих линеаризацию Н Тезисы докладов Республиканской научной конференции "Дифференциальные уравнения и их приложения". Куляб, 1991. С. 32. .

18. Лаврентьев М.М.. Бронников A.B.. Воскобойников Ю.Е..Зеркаль С.М.. Хогоев Е.А. Сейсмическая томография сред с квазилинейным изменением скорости, содержащих поглощающие включения // Изв. РАН. Сер. Физика Земли. 1995. №6. С. 26-31.

19. Лаврентьев М.М.. Воскобойников Ю.Е.. Зеркаль С.М.. Мухометов Р.Г. Сейемотомография в кинематических исследованиях глубинного строения Земли // Методы оптимизации и их приложения. Тезисы докладов. Иркутск, 1995. С. 268.

Приоритетными являются работы, опубликованные без соавторов.

Рис. 1.

Система наблюдений, формирование лучами поверхности шарового сегмента

Рис. 2.

Геометрическая иллюстрация к выводу формулы обращения

Рис. 3.

Решение обратной задачи в случае смещенной локальной неоднородности