Численный метод и алгоритм решения обратных коэффициентных задач акустического зондирования функционально-градиентных материалов тема диссертации и автореферата по ВАК РФ 05.13.18, кандидат наук Темьянов Булат Каримович

• Темьянов Булат Каримович
• кандидат наук
• 2017

Оглавление диссертации кандидат наук Темьянов Булат Каримович

отражателей

1.1.1 Лучевая и дифракционная томографии

1.2 Алгоритмы ультразвуковой томографии, для средних и сильных отражателей

1.2.1 Итерационные алгоритмы восстановления неоднородностей

1.2.2 Функционально-аналитические алгоритмы восстановления неоднородностей

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

1.4 Метод граничного управления

1.5 Постановка цели и задач исследования

ГЛАВА 2. Численный метод определения функции профиля модуля упругости акустической среды

2.1 Постановка обратной коэффициентной задачи

2.2 Основной измерительный алгоритм

2.2.1 Основные соотношения и этапы измерительного алгоритма

2.2.2Математическая модель одномерной неоднородной акустической среды. Вычисление акустического адмиттанса и ядра интегрального уравнения

2.2.3 Модификации ядра и правой части интегрального уравнения

2.2.4 Метод оценивания начального приближения основного измерительного алгоритма

2.2.5 Регуляризация и численный алгоритм решения интегрального уравнения

2.3 Численная реализация алгоритма определения распределения модуля упругости неоднородной акустической среды

2.3.1 Программная реализация измерительного алгоритма

2.3.2 Результаты численного моделирования измерительного алгоритма

2.4 Выводы

ГЛАВА 3. Экспериментальное исследование численного метода определения функции профиля модуля упругости неоднородной акустической среды

3.1 Цели и задачи эксперимента

3.2 Экспериментальное оборудование

3.2.1 Ультразвуковой дефектоскоп TomoScan FOCUS LT

3.2.2 Датчик на основе ультразвуковой фазированной решетки

3.3 Измерение акустического адмиттанса неоднородной акустической среды

3.4 Экспериментальное исследование измерительного алгоритма

3.4.1 Схема и методика эксперимента

3.4.2 Результаты экспериментального исследования измерительного алгоритма на первом образце

3.4.3 Результаты экспериментального исследования измерительного алгоритма на втором образце

3.5 Выводы

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

акустической среды

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

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

4.3 Многомерный алгоритм определения распределения модуля упругости неоднородной акустической среды

4.4 Оценка пространственной разрешающей способности измерительного алгоритма

4.5 Выводы

Заключение

Список литературы

Приложение

Приложение

Приложение

з

Введение

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

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

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

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

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

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

Цель и задачи исследования

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

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

Научная задача исследования - разработка алгоритма решения обратной задачи акустики для сред с непрерывными переменными по одной пространственной координате (глубине) механическими свойствами во всем объекте зондирования. Для достижения поставленной цели и решения научной задачи необходимо решение следующих задач:

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

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

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

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

Научная новизна работы

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

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

Практическая ценность

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

Реализация результатов работы

Результаты диссертационной работы использовались:

1. В научно-исследовательской работе, проведенной на кафедре Технологического оборудования медицинской и легкой промышленности ФГБОУ

ВО «Казанский национальный исследовательский технологический университет» в рамках федеральной целевой программы «Научные и научно-педагогические кадры инновационной России» по теме: «Современная высокотехнологическая ультразвуковая аппаратура для медицинских исследований» (шифр заявки «20101.1-235-073-015) по Государственному контракту от «11» июня 2010 г. №02

2. В учебном процессе кафедры Радиоэлектроники и информационно-измерительной техники (РИИТ) ФГБОУ ВО «Казанский национальный исследовательский технический университет им. А.Н. Туполева - КАИ» в рамках дисциплины «Программные комплексы» при подготовке магистров по направлению 11.04.01 «Радиотехника»

Научные положения, выносимые на защиту

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

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

2. Измерительный алгоритм, реализующий предложенный численный метод.

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

Введение диссертации (часть автореферата) на тему «Численный метод и алгоритм решения обратных коэффициентных задач акустического зондирования функционально-градиентных материалов»

Апробация работы

Результаты работы докладывались и обсуждались на международных, всероссийских и республиканских научно-технических (НТК) и научно-практических конференциях (НПК): на двух международных НПК «Инженерные, научные и образовательные приложения на базе технологий National Instruments» (Москва, 2011, 2012гг.); на двух международных НТК «Проблемы техники и технологий телекоммуникаций» (Казань, 2011, 2014 гг.); на региональной НПК «Современные методы электрофизической диагностики» (Казань, 2013г.); на

международной НТК «Нигматуллинские чтения» (Казань, 2013г.); на Всероссийской НТК «Информационные технологии в электротехнике и электроэнергетике» (Чебоксары, 2012 г.); на трех Всероссийских НТК «Динамика нелинейных дискретных динамических систем» (Чебоксары, 2013, 2015 2017 гг.); на международной НТК «Прикладная электродинамика, фотоника и живые системы» (Казань, 2015 г.).

Публикации

Основные научные и практические результаты диссертационной работы опубликованы в 20 работах, в том числе в 5 статьях (в рецензируемых журналах, рекомендованных ВАК), 15 - в сборниках материалов международных и всероссийских конференций.

Соответствие содержания диссертации паспорту научной специальности

Диссертация соответствует паспорту специальности 05.13.18 «Математическое моделирование, численные методы и комплексы программ»:

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

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

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

Структура и объем работы

Диссертационная работа содержит введение, четыре главы, заключение, список литературы из 118 наименований и три приложения. Материал изложен на 260 страницах, содержит 81 рисунок и 4 таблицы в тексте.

ГЛАВА 1. Численные методы и алгоритмы ультразвуковой диагностики материалов с непрерывными переменными механическими свойствами

1.1 Алгоритмы ультразвуковой реконструктивной томографии для слабых

отражателей

Традиционные методы ультразвуковой дефектоскопии [6] занимаются определением положения несплошностей в пространстве объекта контроля по отраженным от них ультразвуковым импульсам. Такое направление в дефектоскопии получило название отражательной акустической томографии [7]. Акустическая томография в данном случае понимается как направление акустики, занимающееся формированием (получением, восстановлением) акустических томограмм, т.е. масштабных изображений акустических характеристик в определенных сечениях трехмерных объектов [7,8]. К акустическим характеристикам относят параметры объектов, влияющие на распространение, отражение, рассеяние акустических волн (плотность, скорость звука, отражающая способность и т. п.).

Следует отметить, что в отражательной томографии, как и в ряде других методов (методы синтезированной апертуры), вообще нельзя установить вид отображаемой неоднородности. Считается, что фиксируется эффективная отражательная способность, хотя это очень условное и нестрогое понятие [7]. Нас же интересует определение распределений акустических параметров, характеризующих механические свойства (плотность, модуль упругости) рассматриваемой зондируемой среды. Подобные параметры определяются (восстанавливаются) в реконструктивной акустической томографии.

1.1.1 Лучевая и дифракционная томографии

Реконструктивная акустическая томография имеет свои истоки в рентгеновской томографии, отсюда же возникла трактовка томографии как

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

1. Учет эффектов дифракции.

2. Облучения объекта с разных пространственных положений (многоракурсность).

3. Использование первых приближений Борна или Рытова.

4. Сведение обратной задачи определения неоднородности к ряду преобразований Фурье или смежных преобразований (вычисление свертки или преобразование Радона)

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

Выделим основные семейства алгоритмов лучевой и дифракционной реконструктивных томографий. В лучевой реконструктивной томографии различают три основных семейства алгоритмов восстановления функции неоднородности: алгоритмы, основанные на разложении в конечные ряды (series

expansion methods) [9-21], алгоритмы, основанные на прямом преобразовании Фурье (direct Fourier transform algorithms) [22-29] и алгоритм обратного проецирования с фильтрацией (filtered backprojection algorithm) [30-37].

В дифракционной томографии также существуют три семейства алгоритмов [7]: алгоритмы, основанные на разложении в конечные ряды [38] и алгоритмы, основанные на интегральных преобразованиях, которые в свою очередь делятся на два семейства в зависимости от того в пространственной [39-43] или в частотной [39,40] области проводится реконструкция изображений.

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

За последние 30 лет появились новые алгоритмы решения обратных задач акустики, преодолевающие ограничения первых приближений Борна и Рытова, позволяющие восстанавливать среднеконтрастные и сильноконтрастные неоднородности в изучаемых объектах. Учитывая это, мы начнем рассмотрение численных методов и алгоритмов ультразвуковой диагностики с этих методов.

1.2 Алгоритмы ультразвуковой томографии, для средних и сильных

отражателей

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

внешнего источника /0(к,г)и опустив при этом зависимость всех функций от волнового числа к:

V2 £(г)+ к2 8 (г) = е(г)8 (г) + /0(г), (1.2.1)

где £ (г) - акустическое давление, б(т) - функция неоднородности. Применяя к уравнению (1.2.1) функцию Грина 00(г,г')получим следующие уравнения

8 (г) = (г) + \ 00 (г, г ' )е(г') 8 (г ' )йг', (1.2.2)

К

8т (г) = \ О0(г, г') У0(г), (1.2.3)

я

где К - область, содержащая неоднородности, Я - поверхность на которой находятся источники излучения. Уравнение (1.3.2) в теории рассеяния известно как уравнение Липпмана - Швингера [44]. Как отмечено в [44] полное поле £ (г) известно лишь на приемной апертуре г еУ, в противном случае, если бы поле было известно всюду, то уравнение (1.2.2) представляло бы собой уравнение Фредгольма первого рода, что упростило бы его решение. Учитывая, что 8 (г) зависит от е(г), задача их нахождения из уравнения (1.2.2) является нелинейной.

Уравнение (1.2.2) можно более наглядно представить системой следующего вида [44]

8(У) = 8 т(у) +1О0(У,г')е(г')8(г ')йг', у е У, (1.2.4)

К

8(г) = 8т(г) + \О0(г,г')Б(г')8(г')йг', г е К. (1.2.5)

К

Система (1.2.4) - (1.2.5) является наиболее общей формой записи обратной задачи рассеяния [44].

Для упрощения дальнейших рассуждений, введем три интегральных оператора Р0, Q0, К0 как было сделано в [44]. Оператор Р0 переводит источники поля f0(г) в первичное излучаемое поле 8 0(г)

80(г) = Р0/0(г') = /О0(г, г')/0(г')йг', г е К. (1.2.6)

я

Оператор К0 переводит вторичные источники поля ^(г) в поле 8(г)

#(г) = Я0Яг') = |в0(г, г')Г(г')йг\ г е Я. (1.2.7)

л

Оператор Q0 переводит источники, возникающие в области Я, в рассеянное поле 8(у), у е У на приемной апертуре

8(у) = ^^(г ') = \в0(у,г')^(г ')йг', у еУ. (1.2.8)

Я

Выражения (1.2.4), (1.2.5) с использованием введенных операторов могут быть записаны в компактной форме

8 - 8 п = , (1.2.9)

8 - 8 ы = Я088, (1.2.10)

при этом 8 и 8 п в левой части (1.2.9) зависят от у е У, а в левой части уравнения

(1.2.10) зависят от г е Я. В правой части обоих уравнений £ и 8 зависят от г е Я.

Используя введенный способ обозначений систему (1.2.9) - (1.2.10) можно

записать в альтернативной форме [44]

8 = - PJ0, (1.2.11)

где Е - единичный оператор. Как видно из (1.2.11), выражение в квадратных скобках при £Я01| < 1 можно разложить в ряд по итерированным ядрам

ад

8 = 8¡п + +... = £ ^(^Г18 п. (1.2.12)

п=1

Такой ряд также называется рядом Борна - Неймана [8,44,45]. В случае слабых рассеивателей, то есть если |£Я0| << 1, достаточно удержать только первый член.

Алгоритмы восстановления в таком случае сводятся к уже рассмотренным алгоритмам дифракционной томографии. В случае рассеивателей средней силы одного члена в разложении уже недостаточно, задача становится нелинейной и решать ее нужно с помощью других подходов, если же рассеиватель сильный, то есть ||£Я0|| > 1, то ряд Борна - Неймана расходится. В таком случае приходится

расширять область применимости алгоритмов для рассеивателей средней силы [44] или создавать алгоритмы на новых принципах.

1.2.1 Итерационные алгоритмы восстановления неоднородностей

При рассмотрении рассеивателей средней силы необходимо учитывать многократные рассеяния. Одним из возможных подходов их учета является применение итерационных методов. В работах [44,46-48] приводится итерационный метод нахождения функции неоднородности на основе попеременных решений уравнений (1.2.9), (1.2.10). Перепишем их в более подходящей форме с учетом изменения частоты с и вариации апертуры датчика

g (а, р, с) = S0(P,o)sg (а, с), (1.2.13)

g (а, с) = R (с)е g (а, с) + g0(a, с), (1.2.14)

где а и р - конфигурационные параметры, ответственные за пространственные изменения апертур в процессе зондирования [44], a gs (r) = g(r) - g in(r) при r e Y. Алгоритм в данном случае заключается в следующем: на каждой итерации выражение (1.2.14) разрешается относительно g(а, с) при фиксированном s, после этого решается (1.2.13), то есть находится s при найденном g(а,с). В работе [44,47] предлагаются также и модификации этого алгоритма.

Уравнение (1.2.2) можно также переписать в терминах операторов T(k, k0s0") переводящих плоскую падающую волну с волновым вектором k0s^ в рассеянную плоскую волну с волновым вектором k1 = (kx, ky):

+ад

T(k1,k0s0) = s(k1 -k0s0) + ¡G0(k[)T{k[,k0s"0 )s(k1 -k[)dk[, (1.2.15)

—ад

где s(k1) и G0(kj)- пространственные спектры функций s(r) и G0(r) соответственно, s0 = s0(а) = (cosа, sin а) - единичный вектор, характеризующий направление распространения падающей волны, k0- волновое число падающей волны. В таком случае итерационный процесс решения обратной задачи сводится к оценке s(kl)из экспериментальных данных и данных о T(k1,k0s0"), а затем к

оценке T(k1, k0 s0°) из уже найденных значений s(k1), в соответствии со следующими выражениями [49]

+ад

~т) (ку0 — к0 з" )= / к к0 з") ~0(к[)Г. (к[ , к0 З^к — к'Ж, (1.2.16)

—ад

+ад

Т]+1 к, Кз" ) = у5 (к1 — к0< ) +| 0~0 (к; )Т (к[, к0з" )8(Г> к — к[)бк' , (1.2.17)

—ад

где з0а - единичный вектор, характеризующий направление распространения рассеянной волны, / (к0 з0а, к0 з") - амплитуда рассеяния, полученная в эксперименте, I и у - номера итераций для функции неоднородности ¿~(к1) и оператора Т(к1, к0з") соответственно. Достоинством этого варианта является возможность широкого применения быстрого преобразования Фурье (БПФ) [44].

В работах[44,47] предложена следующая итерационная процедура решения обратной задачи рассеяния, основанная на (1.2.11):

8з (", а, ю) = Яе(Д ю)Е — (ю)]1 8т (", а)£{. (1.2.18)

Сходимость данного алгоритма определяется модулем отношения рассеянного поля к падающему в каждой точке области рассеяния и имеет место при 18з (г)/ 8(г)| < 1. При этом характер сходимости экспоненциальный.

I N

Относительная ошибка оценки е(т) после N шагов итераций ~ 8Я (г)/ 8(г)

[44,47]. Также рассмотрен способ расширения области сходимости данного алгоритма, позволяющий восстанавливать сильные рассеиватели [47,50].

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

Необходимо отметить, что итерирование системы двух уравнений (1.2.13) -(1.2.14), либо (1.2.16) - (1.2.17) имеет то преимущество, что не требует трудоемкой и неустойчивой операции решения СЛАУ [44]. Вычислительные операции также облегчаются сверточным характером операторов Я0, Р0, К0 [44].

Сложность данного вида алгоритмов заключается в нахождении большого количества вспомогательных переменных - рассеянного поля внутри рассеивателя [44]. Итерирование же уравнения (1.3.18) приводит к необходимости обращения матрицы СЛАУ на каждом итерационном шаге, хотя

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

В работе [46] кроме рассмотренного выше алгоритма попеременных итераций предлагается алгоритм решения обратной задачи рассеяния на основе задачи минимизации функционала

2 и 2

Пб) = 1/2£|*, - V*! = 1/2Е к, , (1.2.19)

-/0t

о,ф о, ф

где ф- угол излучения плоской волны, при условии удовлетворения g(r) уравнению (1.2.5) либо (1.2.10), выступающих в качестве ограничений для минимизации. Поиск минимума (1.2.19) производится методом Ньютона - Гаусса, согласно которому на каждой его n - ой итерации решается задача минимизации

min УI\j(n) 8s(n] (r) + О 2, (1.2.20)

О,ф

где Ss n (r) - поправка к функции неоднородности s(r) на n -ой итерации, J матрица Якоби функции r

„( n ) ' о,ф

8г (")

•С, = = -^0 [е - я0 (с)б(п) (г)]-1 *(п) (г, с, р), (1.2.21)

ОБ (г)

где Бп)(г) и *(п)(г,с,р) - диагональные операторы, полученные из вектор-функций Бп)(г) и *(п)(г,со,р) при фиксированных со и р.

Таким образом, алгоритм, приведенный в [46], состоит из следующих шагов:

1. Решение прямой задачи для уравнения (1.2.5) либо (1.2.10), нахождение полного поля *(п)(г,с,р) по известным начальным или полученным из предыдущей итерации значениям Бп-1)(г);

( ) II 2 2

2. Вычисление значений невязок Ср. Если ^ С , то алгоритм

останавливается и Б п-1)(г)- искомая величина. В противном случае производится переход к третьему шагу;

3. Решение задачи (1.2.20);

4. Вычисление следующего приближения s(n\r) = s(n_1)(r) + Ss(n)(r). Данный алгоритм сходится даже в том случае, если алгоритм, основанный

на попеременных итерациях, уже не может быть применен из-за сильной контрастности или больших размеров неоднородности [46]. Проведена также оценка сложности алгоритма, которая зависит от размерности массива функции неоднородности s(r), количества углов р, количества итераций метода Гаусса -Ньютона, количества итераций метода решения прямой задачи и вычисления обратного оператора в (1.2.21) и т. д. [46].

В работах Колтона, Монка и Кресса [51-53] представлен алгоритм решения обратной задачи рассеяния, в котором неоднородность показателя преломления n(r), определяется по данным рассеяния в дальней зоне. При этом

неоднородность находится в пределах сферы с радиусом a, где a - некоторое действительное число, то есть n(r) = 1 при |r| > a. Для работы алгоритма

необходимо измерить диаграммы направленности F (г, к, v) рассеянного излучения вдоль направлений задаваемых вектором r = r /| r| при облучении

объекта измерения плоской волной направленной вдоль единичного вектора v. В соответствии с алгоритмом требуется также вычислить диаграмму направленностиFA (Г,к,v) решения h(r) следующей вспомогательной краевой задачи

Ah + к2h = 0, при Q0, (1.2.22)

h(r) = exp [/kr • v]+ hs (r), (1.2.23)

dh

— + iAh = 0, при аедП0, (1.2.24)

да

(dhs Л

lim а--ikh = 0, (1.2.25)

^ да )

где а = |r|, Q0 = (r: |r| < b, b > a), A- параметр, выбираемый из условия существования для его значения решения задачи (1.2.22) - (1.2.25), -

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

;р-1

1 (г, к, V) - ^ (г, к, V )]е рд (г)бу(г) = —Ур (г),

02 к

(1.2.26)

где ^ (г) - так называемое ядро Герглотца [80-82], Ур (г) - сферические функции

[51,53] с параметрами р = 1,Р и р = - р, р, 5 - единичная сфера, б&(г) - элемент площади. На втором шаге определяется функция ш(т) = 1 - п(т), характеризующая неоднородность, решением нелинейной задачи минимизации функционала

Р р N

I = Т Е Е

р=1 р=- р г=1

да

+ ¡л^ррр - к? 5 ш^рр + К )

(1.2.27)

Ь2 (дПо)

где г = 1, N - номера учитываемых волновых чисел к, функции ирр, ирр задаются следующим образом

ирр(т) = Iехр(- гкт■ б(б)Ж(б), б е^3, (1.2.28)

5 2

ир (т) = И{1)(ка)Ур (г), (1.2.29)

где Ярх)(ка) - сферическая функция Ханкеля 1 рода, р -го порядка, ирр (т) также называется волновой функцией Герглотца [51-53]. Функция wpq удовлетворяет уравнению Липпмана - Швингера в следующем операторном виде

w =и - к2 Я0 w , т еО0.

рр рр 0 рр 0

(1.3.30)

Задача минимизация выше упомянутого функционала осуществляется методом Левенберга - Марквардта [51].

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

2

сферической функции в последнем численно значимом члене разложения данных рассеяния в дальней зоне в ряд по сферическим функциям [51-53].

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

В работе Гутмана и Клибанова [54] представлен квазиньютоновский метод решения обратной задачи рассеяния отличающийся систематическим использованием уравнения Гельмгольца в отличие от ранее рассмотренных алгоритмов, основанных на уравнении Липпмана - Швингера. Данный метод рассмотрен в применении к слабоконтрастным неоднородностям, хотя только его первая итерация совпадает с приближением Борна, а дальнейшие итерации лишь улучшают качество получаемого приближения. В качестве основного исходного выражения используется следующая краевая задача:

где V(г)- искомая функция, б(г)- неизвестная функция неоднородности, Ф(г)-функция, специальным образом построенная из данных рассеяния, Г - граница трехмерного куба , окружающего неоднородность [54]. Искомые функции б(г ) и V (г) определяются как предел последовательности следующего итерационного процесса

(1.2.31)

(1.2.32)

(1.2.33)

Л(б, V) = AV + к V + к 2б(У + Ф) + (АФ + к 2ф). (1.2.35)

В качестве начального приближения используются следующие выражения

й0 =(0, и0 - Ф), (1.2.36)

и0 = ехр(/к(г, у)), (1.2.37)

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

направлена плоская волна (1.2.37). Учитывая выражение (1.2.35), соотношение (1.2.34) можно переписать в следующем виде

АУ + к У + к 2епи0 = к 2Бп_хи0 - к Ч-! (У- + Ф)- (АФ + к Ф). (1.2.38) Для дальнейших вычислений выбиралось положительное целое число N такое, что ^ < 2к (ЛN < 2к для двумерного случая). Затем формировалось множество Р точек X е^3 таких что, координаты г1, г2, г3 принадлежащие Z и

Р, были целочисленными и удовлетворяли неравенству |г;| < N, [ = 1,2,3. Также фиксировалось положительное число д такое, что д2 + 2к2 было нецелым числом. Аналогично координатам г использовались также и координаты У = (у1, у2, у3)е Я3, такие что |у;| < М , [ = 1,2,3, где м - некоторое положительное целое число.

Для каждого X е Р выбирался единичный вектор уг е S2 удовлетворяющий соотношению:

2к(уг, ^ = 2. (1.2.39)

При умножении (1.2.38) на ехр(/(х, X - куг))и интегрировании получающегося выражение по объему П1 с учетом граничных условий (1.2.32), (1.2.33), и (1.2.39) было получено следующее соотношение

Б = Я (г)ехр((г,)йг = -11/(г)ехр((г,X - куг))йг, (1.2.40)

к

/(г) = к2Бп_х(г)и0(г) - к2Бп-1 (г)УП-1 (г) + Ф(г))- (АФ(г) + к2Ф(г)), (1.2.41)

то есть n -ое приближение функции s(r) можно вычислить по его Фурье-коэффициентам. Функция Vn (r) определялась в форме

Vn (r, vz) = exp(fc + ik)r, vz)W(r, vz). (1.2.42)

При умножении (1.2.38) на exp(-(^ + ik)r,vz) + i(r,Y)) и интегрировании данного соотношения по объему Qp было получено выражение для Фурье-

коэффициентов функции Wn (r, vz)

WY

P(k, Y, vz)

k2 jen exp(- g(r,vz) + i{r, Y))dr +

где

j f (r) exp(- + ik X r, vz) + ^r, Y) )dr

p(k ,Y, vz )= 2k(vz , Y) - 2i^ vz , Y) - |Y|2 + + 2igk

(1.2.43)

(1.2.44)

Тогда искомая функция Vn (r, vz) определяется следующим выражением

Vn(r,vz) = ^expfe + ik)r,vt))EWn1, exp(- i(r,Y>).

алгоритма проводились для

(1.2.45) функции

£(r) =

(1.2.46)

Численные испытания неоднородности вида

а, при Ь < г < а Р, при 0 < г < Ь , 0, при г > а

Относительная погрешность решения при 12 итерациях и различных значениях параметров (1.2.46) лежала в диапазоне от 1,22% до 3,86%. Значительная часть работы [54] посвящена математическому обоснованию алгоритма.

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

1

осуществлялось в N конфигурациях. Функция неоднородности определялась как линейная комбинация

м

б(г ) = ) , (1.2.47)

т=1

где \^т(г), г е^2,т = 1,М] - семейство линейно независимых функций выбранное либо исходя из априорной информации о функции неоднородности б(г), либо, если таковой информации не имеется, в виде полиномов или функций формы метода конечных элементов. Аналогичным образом вводится второе семейство линейно независимых функций {^2(у), у е У, ] = 1, J} имеющих область

определения на граничной поверхности У, на которой происходит измерение. Таким образом, обратная задача сводится к определению М коэффициентов разложения б = (б1,б1,..ем )т . Рассеянное поле разлагалось по базису функций

Ф2( у)

gsi = {{ gsiЖ) Y ,( Y '•••'( Y J, d-2-48)

где (y) - скалярное произведение по области определения Y [55]. Также вводились следующие интегральные операторы

m _

R(s) gt = 1 eM gt, i = 1, N, (1.2.49)

m=1

L(e)g > = g > - RB)g >, i = 1, N, (1.2.50)

s(ft f = {{So gi£-4>ï)Y ,(S g^2)Y ,-(So ge-ФЖ), J. (1.2.51)

Ha n -ой итерации рассматриваемого алгоритма приближения gin и en вычислялись следующим образом

~ ~ гпс Л г\

8I^o = 8г , Б0 = 0

8{,п = 8i,n-1 + ^ г 1n-1, Бп = Бп-1 + А^ (1 .2.52)

г> = - Цеп ), г2п = г - ^)Бп ,

где 8ilnc - первичная излученная волна в г -той конфигурации, а коэффициенты а\ и р1п выбираются из условия минимума выражения

^г.1 2 Т г.'

|| ¡.п ^2 Ц ■,

У

^ = -¡=-+ -. (1.2.53)

п Пс 2 ^п „2 У ^

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

г

б 1 =

5 5 5

дБ1 д£2 5бм у

N ,,

Т\\>;2

¡■=1

. (1.2.54)

81=81 ,п-1,Е=Еп-1

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

Численное испытание алгоритма проводилось для разных типов среднеконтрастных неоднородностей как непрерывных, так и разрывных кусочно-постоянных. При этом погрешность для непрерывных функций неоднородностей составляла менее 1%, для разрывных - около 8%. Погрешность решения уменьшалась при повышении отношения размеров зондируемой среды к волновому числу излучаемой волны. Представлен анализ поведения решений итерационного метода при добавлении к модельным данным шума со среднеквадратическим значением в 10%, при этом погрешность нахождения функции неоднородности относительно слабо зависит от добавленного шума, изменение нормы погрешности составило около 1% при общей погрешности равной примерно 8%. Изучался вопрос о частотном ограничении спектра функции неоднородности разработанным алгоритмом, для этого на первом этапе восстанавливался разрывный кусочно-постоянный профиль, полученное приближение затем использовалось в качестве искомого для работы алгоритма на втором этапе. Погрешность восстановления для разрывного профиля составила около 6%, а для частотно-ограниченного - менее 3%.

В работе Наттерера и Вюббелинга [56] представлен итерационный алгоритм отличающийся простотой и достаточно малой сложностью. Данный

I

алгоритм есть обобщение алгебраического алгоритма реконструкции (ART), который уже был нами рассмотрен, на нелинейный случай. Согласно рассматриваемому алгоритму формируется множество единичных векторов

Список литературы диссертационного исследования кандидат наук Темьянов Булат Каримович, 2017 год

(или источники)

Рассеиватель окружен двумя

Источники замкнутыми, всюду гладкими,

(или приемники) „

4 ^ у выпуклыми контурами О и 12 г

уеОу у

произвольной формы, на одном из

Область рассеяния Я

--которых размешены точечные

Рис. 1 Геометрия задачи решаемой источники, а на Другом точечные

методом Марченко - Ньютона - Роуза „ т/- I ¿ ¿ I

у приемники. Как указано в [66], по

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

Грина [66,67] неоднородной среды подчиняются уравнениям Гельмгольца с правой частью равной 8 -функции Дирака:

А0+(а,г,+ 0+(а,г,у) = 8(г- у), у еПу, (1.2.68)

с (г) у

2

,2

АО (а, г, х) + -0- О (а, г, х)=8(г - х), х е Я. (1.2.69)

с (г)

Следствием этих уравнений является уравнение МНР для полного поля в частотном представлении

О (а, у,х)-0+(а,х, у)= {йа |д0 (а, г у 1 о (а, г,х)-

м дп

(1.2.70)

/ \дО (а, г, х)

0 (а, г, у)-^-

дп

где О±(а,х,у)и О (а,г,х) - искомые поля, О+ (а,г,у)- рассеянное поле,

д0(а, г,■) , 0 . - - производная по внешней нормали (по отношению к контуру 12г)

дп

г

нормали пг в точке г, йа - элемент площади. Учитывая, что искомые поля

О± (а, х, у) и О- (а, г, х) имеют различные аргументы, требуется применить так

называемую процедуру продления полей в пределах однородного пространства [66]. Для этой процедуры осуществляется переход к угловым гармоникам полей уравнения (1.2.70) [66]. Во временном представлении для рассеянных полей g ±= G±- Go; (G0; - падающее поле), при S- образном внешнем воздействии уравнение (1.3.70) принимает вид

g(t, у,х) - g+ (t,х, у)-Jdrjda g*(t-T, z,x) -

- G± (r, z, y) = daz j^^^ Gj (t-r, z, x) - (1.2.71)

dnz J Qj l dnz

; 5Go+ (t -r, z, x)) - g - o, z, y) —^-

dnz J

Дополнительно используется причинная связь

I . I r rr\

\r - r \ \r - r

g + (t, rr") = 0 при Vt <r+ -L, g ~ (t, r', r") = 0 при Vt >r = J-L. (1.2.72)

max c(r) max c(r)

r r

Также в работе [66] предложена модификация алгоритма МНР, которая допускает любой спектр облучающего поля. Конечные выражения уравнений МНР и условий причинности после дискретизации приводится к плохо обусловленной системе линейных алгебраических уравнений, решать которые можно только лишь с применением регуляризации. Как показало численное моделирование, алгоритм МНР обладает неединственностью решений. Полученные оценки решений отличные от точных также являются решениями. Следует правда отметить, что в работе [66] не указано каким принципом выбора параметра регуляризации руководствуются при решении уравнений МНР, указано лишь то, что он определяется экспериментально, при этом, если этот параметр выбран неверно, полученные оценки могут существенно отличаться от точных решений. Также в [66] предпринимается попытка устранения неединственности путем использования дополнительных уравнений связи типа Липпмана -Швингера при сохранении линейности решаемой задачи, но неоднозначность решения подобным путем также не устраняется.

В последующей части работы [66] представляются модифицированный алгоритм МНР, а также второй основной функционально-аналитический алгоритм данной работы - алгоритм Новикова - Гриневича [66-69]. Для их рассмотрения представляется формализм комплексных волновых векторов, являющийся неотъемлемой частью функционально-аналитических методов. При этом обратная задача рассеяния предполагает восстановление скалярной функции неоднородности вида

г л а2 а2 аа (г, а)

б(г) = — --^г - 12а а , (1.2.73)

С с (г) с(г)

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

к = к Я + к,, 1 = 1Я + й1, к Я 1 к ,, к2 = к2, кЯ ,к, е^п , (1.2.74)

где Шп - множество пар или троек действительных чисел при п = 2 или 3 соответственно. С физической точки зрения этот переход означает переход от плоских однородных волн к неоднородным. В двумерной задаче имеется только две ортогональные ориентации ненулевой мнимой части волнового вектора относительно его действительной составляющей. Переход к комплексным волновым векторам приводит к обобщению функций от них зависящих. Обобщенное волновое поле £ (г, к) при этом удовлетворяет следующему уравнению:

а£(г,к) + к2£(г,к) = \(г)£(г,к), ке Сп, к2 = к2, (1.2.75)

где Сп -множество пар или троек комплексных чисел при размерности п равной 2 или 3 соответственно. Используется также прием снятия несущей частоты [66-68] с помощью введения функции /л(г,к) = ехр(-/кг) £ (г,к), подчиняющейся уравнению Липпмана - Швингера следующего вида

((г,к) = 1 + | Оц(г - г' ,к)((г ' ,к) £ (г ,к)йг', (1.2.76)

г 'еЕ"

где Ом(г,к)- функция Грина, введенная Л. Д. Фаддеевым [66-69]. Классические

значения О0± (г,к = кЕ) функций Грина получаются предельным переходом при

|к7| ^ 0, если кЕ и кг сонаправлены, либо направлены противоположно.

Процедура восстановления характеристик двумерного рассеивателя в соответствии с модифицированным алгоритмом МНР состоит из следующих этапов [67]:

1. Расчет обобщенной амплитуды рассеяния Н± (к, 1) = Н (ф, ф') (при к 1 = к± ^ 0) из экспериментальных данных

ч i

2л

h+ (ф,Ф') -Л jh+ (ф,ф') ■ вн[+ sin{f - (фп,<p')dp" = f {ф,ф'), (1.2.77)

0

где вн - функция Хевисайда. Уравнение (1.2.77) линейно относительно неизвестной функции h±, учитывает перерассеяния волн на неоднородности среды и нелинейно относительно данных рассеяния f (p,p).

2. Из h ±(p, p) вычисляется функция p(k,l) = p(p,p) - ядро в интегральном уравнении связи предельных значений g± (r,p) волновой функции g(r,k). Данное уравнение также называется р - соотношением и является модифицированным уравнением МНР [66,67]

2л

g+(r, ф) - g (r, p) = jp(p,f)g - (r,f)df, (1.2.78)

0

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

2л

р{ф,ф') + Л |р(ф,ф" )fi(ффвн [sin(ф - ф W = Л(фф), (1.2.79)

0

fi(ф, ф) = вн [- sin(ф - ф)У (ф, ф)- вн [sin(ф - ф)]ь (ф, ф). (1.2.80)

Важнейшую роль в алгоритме МНР играет соотношение Сохоцкого [66,67] ( k = [kx, ky }е C2). Оно характеризует связь между полем ¿u(r,k) = ¡л(г,ф) и его

предельными значениями ¡л±

1 2гл{и+ (р -м (г,р)}ехр(р")р

Мг,к) = 1 + -^Т{м+к^^УУ^^у (1.2.81)

2л 0 (. „) к + к * ехр(р )--

к

Тогда линейная система, состоящая из модифицированных уравнений МНР (1.2.78) и одного из уравнений Сохоцкого (1.2.81),рассматриваемого для ¿и^ М± (при к 1 = к± ^ 0), обеспечивает единственность восстановления внутренних полей [66,67]. Данный алгоритм имеет следующие преимущества: а) сохраняется свойство линейности относительно неизвестных полей и или их угловых гармоник, б) эта система не требует знания «нефизических» данных рассеяния в отличие от (1.2.70) - (1.2.72), в) система однозначно разрешима, что обеспечивается неоднородностью уравнений Сохоцкого и тем обстоятельством, что оно выделяет из возможных решений только то, которое удовлетворяет требуемым свойствам аналитичности.

На конечном этапе проводится восстановление функции неоднородности рассеивателя из уравнения Липпмана - Швингера или уравнения Гельмгольца

*(г) = к2 + А±±(гу). (1.2.82)

(г, р)

Из-за неточности в вычислении £1 (г,р) оценка е(г) может зависеть от ракурса р, простейшим выходом из этой ситуации является усреднение по углу р, однако есть возможность более органичного объединения информации полученной с различных направлений облучения. Такую возможность предоставляет алгоритм Новикова - Гриневича, первые два этапа которого совпадают с этапами модифицированного алгоритма МНР. Поэтому рассмотрение алгоритма Новикова - Гриневича начнем с его третьего этапа:

3. Оцененные ^(к,1) позволяют найти значения разностного поля

К (г, р) = м+ (г, р) - и (г, р)

2л

К(г,р) = | Др,р')ехр(/к{х(со8р' - со8р)}+ у{&тр - 8тр))х

0 - / ч П (1.2.83)

1+±ч к(г,1 ,

2л О ехр(р ) - (1 + 0)ехр(р')

Выражение (1+0) в (1.2.83) означает присутствие бесконечно малой положительной добавки к единице.

4. Поле К (г, ф) связано линейно с искомой функцией неоднородности соотношением

. ч к е(г) = -—

2л

г.д дл I— + —

V дх ду у

2л

IК (г,ф)ехр(ф)ф (1.2.84)

В [66,67] показывается, что « линеаризация» задачи достигается благодаря свойствам симметрии предельных значений обобщенных функций Грина -Фаддеева (г, к) относительно направления вектора к еЖ2 и интегрированию по всем углам падения плоской волны ф. При этом пространственный спектр Т к) = | Е + (г,к) ехр(- к классических внутренних источников

Е + (г, к) = а(г)£+ (г,к) должен быть локализован внутри круга радиуса 2к

Т(£к) = 0 при |£-к| > 2к. (1.2.85)

Особенности алгоритма Новикова - Гриневича заключаются в следующем:

1. Учитываются эффекты многократного рассеяния волн, однако алгоритм остается линейным относительно искомой функции неоднородности;

2. Алгоритм имеет существенное специфическое ограничение - требование отсутствия рассеяния назад. Данное условие есть требование на устойчивость решения двумерной монохроматической задачи рассеяния;

3. Для медицинских приложений диапазон приемлемых частот лежит в пределах от десятков кГц до нескольких МГц;

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

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

0

В то же время алгоритм не лишен и недостатков:

1. Алгоритм не допускает простого обобщения на трехмерное пространство;

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

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

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

2л

h±(р,р';с)-л ÍИг(р,р';с)-вн\+sin(р"-р)]х

7 о 7 . (1.2.86)

х f (, р;со7 )d( = f (р, р;со7), 7 = 1, J Отличием данного уравнения от аналогичного уравнения (1.2.77) является введение частотной зависимости, которое понадобится при его многочастотном (полихроматическом) обобщении.

2. Нахождение классического запаздывающего волнового поля [70,71] gcl (r, k, с) , а точнее модулирующей функции ¡ucl (r, k; с) = exp(-ikr) gcl (r, k; с)

. Для этого последовательно находят вспомогательные функции

И± (r, р, р; с) = И (р, р; с) exp\ik7 {x(cosp - cosp) + _y(sin р - sin р)}], (1.2.87)

Q+ (r, р, (; с) = h+ (r, р, (; с )вн \+ sin(( - р)], (1.2.88)

- 2л

B(r, р, (; с) = - j Q (r, (, (; с)Z+ (р - ()d( -

л 0 , Ц.2.89)

- - jQ+ (r)z (р-()d(

2 о

где x± (() = 1/\1 - (1 + 0)exp(()]. Далее решается система линейных уравнений

2л

jucl (r, р; с7) + л j B(r, р, р; с )^cl (r, р; с )dp = 1. (1.2.90)

о

3. Расчет предельного значения обобщенного поля / (г,ф;с.) по выражению

2л

л (г, ф; с.) = л (г, ф; с.) + л' 10- (г, ф, ф; с. )/ (г, ф; с. )ф. (1.2.91)

4. Вычисление искомой функции неоднородности

д д

£(г,Ю.) = ( ' .) 2л

к С я Я \2л

'- + ■

V

(т,ф,с. )ехр('ф)^ф. (1.2.92)

дх ду

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

рассчитать пространственный спектр вторичных источников Ес (г,к;с..) = а(т,с.)(г,к;с.)и оценить размеры области его локализации.

Численное испытание алгоритма показало, что эффекты перерассеяния приводят к увеличению области расположения пространственного спектра вторичных источников по сравнению с Борновской оценкой. Еще одним преимуществом данного алгоритма является возможность обобщения на многочастотный случай. В соответствии с этим обобщением в расчет алгоритма вводились уравнения связи, вытекающие из общности рассеивателя для различных частот. В простейшем случае одинаковой частотной зависимости Яе б(т,с] ) и \шб(т, с.)

уравнение связи имеет вид:

е(г,с.)/ к. -е(г,с.+1)/. = 0, ] = 1, ] -1. (1.2.93)

В итоге результатом расчета является частотно-независимая оценка функции неоднородности:

йро1у (г) = 1ро1у (г, с)/ к2, (1.2.94)

0

где Vро1у(г,а.)- оценка частотно зависящей функции неоднородности получаемой

из уравнения (1.2.92) с учетом уравнения связи (1.2.93). Для численной проверки данного алгоритма специально решалась прямая задача для цилиндрического рассеивателя, обладающего как рефракцией, так и поглощением. Выбор рассеивателя цилиндрической формы был обусловлен имеющимся точным аналитическим решением прямой задачи для такого случая. Как и в случае с алгоритмом Новикова - Гриневича тестирование показало, что рассеивание назад искажает восстанавливаемую картину функции неоднородности обратной задачи. Для уменьшения его влияния использовалась фильтрация обобщенной амплитуды рассеивателя. В целом, при достаточно большом количестве частот полихроматический алгоритм обеспечивает значительно более устойчивое восстановление, чем монохроматические решения или их среднеарифметическое и предоставляет возможность более органичного объединения многочастотных данных [70,71].

Второй основной алгоритм, рассмотренный в [70-73],называемый трехмерным алгоритмом Новикова - Хенкина, связывает пространственный спектр рассеивателя £(%) = £(-р), где р е^3 - фиксированный вектор,

удовлетворяющий условиям к2 = к2, 2кр = р2, с обобщенными данными рассеяния Н (к,р) = й(к,1 = к - р) следующим соотношением

е(-р) = Н (к. ,р) + ^ [Н (к, р)|^ ^, к. |+ ^ 2 |Н (к, р)|^ ^ ,к ,| (1.2.94)

где

[н (к, р)|^ ,к.|= \ дн (к ' , р) л Кр (к ' , к.), (1.2.95)

к' еМ

Z2[Н(к,р)|к ,к= |{Н+ (к ' ,р) - Н (к ' ,р)}Кр(к ' ,к.). (1.2.96)

к' едМ

В выражениях (1.2.94) - (1.2.96) к„ - произвольное фиксированное значение вектора из множества всех векторов к е С3, интегральное ядро Кр(к ' ,к„) ид -производная дН (к ' ,р) функции Н (к ' ,р) определяются из выражений (черта над вектором означает комплексное сопряжение)

к„(к' ,к.)

ёе1 (р,к' + к,к' -к.) 2т\к' -к .|2

¿к

¿к

¿к'

р к' - р к' р к' - р к' р к' - р к'

Г У 7 у Г 7 X г% 7 ГХ У -ТУ X

(1.2.97)

Бх (р, к ' , к. )йкх + в у (р, к ' , к. )йку + вг (р, к ' , к. )^у = 2 вп (р, к, к. )йкп,

п=1

ъи (к',р) = ЪИ^ к + анск^р) ¿р + жскул) ¿р = анск^р)

лт / X Г У ^ т ^ 7 ^^ ■» ' т ' 4 7

5к'

ъку

ак'

аку

где

а

1

ак 2

а

+ г

а

а(Яе к') а(1тк')

^ т ^ ^ т ^

Верно также следующее соотношение

называемое а - уравнением [70-74]

аи (к ', р) = -2л | (с ¿ку+с/ку+^ )н (к ' и (к '+с,р+о-.

СеИ3

.(1.2.99)

+ 2к '

Индексы п,т = 1,2,3 соответствуют у-, у -,; -компонентам векторов к ' и к' соответственно. Знак л в выражении (1.2.95) обозначает внешнее произведение [75]:

3 3

аи (к у , р) л к р (к у , к.) = аи(ку ,р) вп (р, к у , к. ж л ¿кп,

т=1 п=1

ак

¿к' л ¿к' = 21й(Иек' )<!(1тк'),

т т V т / V т Р

¿к' л ¿к' = ц У(Ие к' У(Ие к') + й(1тк' ^(1тк')] +

т п ' 1 Ь V т / V п / \ т / \ п /Л

(1.2.100)

(1.2.101)

+1

& (Ие ктУ(1ткп)- а(1ткт^(Ие к'п)]

где

1 для пар {т = у, п = у}, {т = у, п = ;}, {т = п = у}, 1 -1 для пар {т = у, п = у}, {т = п = у}, {т = у, п = ;}.

(1.2.102)

В (1.2.96) граница дЫ соответствует действительным векторам к ' е Я3, удовлетворяющим условиям к ' 2 = к2, 2к 'р = р2. Действительные векторы р, кК = к ' и бесконечно малая мнимая добавка к 1 = ±0([р х к ']), определяющие функции И ± (к ' , р) = И (к ' ± Ю([р х к ' ],р), образуют правостороннюю или левостороннюю систему координат соответственно [70].

Как видно из выражений (1.2.94) - (1.2.102), алгоритм очень сложен для практической реализации. Самым сложным в вычислительном плане и к тому же неустойчивым при нахождении является член Z1 выражения (1.2.94). Именно поэтому при численном моделировании алгоритма данный член не рассматривался, что приводит к определенным искажениям в результате восстановления. Более того, моделировался случай сферически симметричного рассеивателя, что означает H+ (k,p) = H- (k,p) и следовательно третий член Z2

равен нулю, поэтому моделирование сводилось в основной своей части к решению уравнений Фадеева (1.2.86) [70-72]. В силу того, что в реализованном варианте алгоритма пространственный спектр s(-p) оценивался на основе

предельных значений H± обобщенной амплитуды рассеяния только для пространственных компонент |p| < 2к, результат восстановления s(r) сравнивался с функцией scut (r), пространственный спектр которой совпадает со спектром истинного рассеивателя s(r) в пределах |p| < 2к и равен нулю вне их [70-72]. За

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

Третий основной алгоритм работы [70] - новый трехмерный алгоритм Новикова основывается на предельном соотношении

s(4 = -p) = lim H (k, p), (1.2.103)

|k|

причем чем сильнее рассеиватель, тем больше должно быть значение |k|, для

которого асимптотика (1.2.103) оказывается верной [70].

Для построения функционально - аналитических алгоритмов использовалась комплексная X- плоскость, между каждой точкой которой Яе Cи

двумерным комплексным вектором к е С , удовлетворяющим условиям к = к , 2кр = р2, существует взаимно-однозначное соответствие. Подобно двумерному случаю и в трехмерии вводится своя Х- плоскость [70]. Для этого формируется множество векторов Ьч сонаправленных оси задаваемой некоторым вектором V. Для любого фиксированного вектора р, непараллельного V, определяется тройка взаимно ортогональных векторов р, 0х(р),ах(р), где \в\ = 1, \ = 1:

V х р

Щ(р)

^р) = £Щ>, p \ Lv,

(1.2.104)

^х р|' |р|

где знак х обозначает векторное произведение. Тогда переход от векторов к е С2 к х - плоскости осуществляется по правилу:

к -[\(р) + ¿юДр)]

Л(к,р) =

(1.2.105)

4k2 - |р2 /4 '

Обратный переход (А,р)^(к,р) осуществляется следующим образом:

кЦ, р) = кх (Л, рЩ (р) + к2 (Л, p)Wi (р) + р/2, (1.2.106)

где

/

(Л, р) = кЩ(р) =

^1(Л,р) = к®1 (р) =

1

л

Л + v Ц

4k2 -I p\2/4

2

/

1

л

-Л vЛ У

i4k2 -I p|2 / 4

2

(1.2.107)

(1.2.108)

В (1.2.105) - (1.2.108) Л^0, р е R3 \ Lv. Далее согласно [70] функция H (к,р) = H (Л,р) рассматривается в переменных (Л,р).

В соответствии с алгоритмом Новикова вначале рассчитываются предельные значения функции H± (к, р) = h± (к, l = к - р), то есть H± (Л,р) = lim H(Л,р) из трехмерного аналога уравнения Фаддеева (1.2.86) ,

Л—>1+0

после чего используются в нахождении функции H (Л,р) уже при произвольном Л е C и |р| < 2k исходя из соотношения:

H (Л,р) = H 0(Л,р) + P({H, H };Л,р). (1.2.109)

Член Н0 (Л, р) вычисляется с использованием следующих выражений [70]:

(1.2.110)

-Л

Л

И0 (Л, р) = | И при 0 <|Л< 1,

2—1^=1-0 д-Л

И,

1 Л

(Л,р) = -— | ---И (др^д, при |Л|> 1. (1.2.111)

2— -1=1+0 д(д - Л)

Член Р({и, И},Л,р) в (1.2.109) получается интегрированием функции {И, И}—р) = аИ(д,р) п0 двумерной области комплексных значений:

ад

Р({И,и};Л,р)=-- |Г{И,И}—p)d(Иед^(1тд), при 0<|Л!< 1, (1.2.112)

— д|<1 д - Л

1 Л

Р({и,И};Л,р) =— |Г---{и,И¡—р)^(Reg)d(1тд), при |Л> 1.(1.2.113)

—>1—(д - Л)

Функция {И, И} находится на основе аналога а - уравнения (1.2.99), записанного в терминах переменных на Л - плоскости:

{и, и Л) = дид-p) = --1

4 - —

ЪЛ

V

к2 - М-(Л - (соф-1)-

(1.2.114)

Н • ф

81П ф

Л

х И (;1,-£) И (;2,р + £)д^фх, Ле С, Ц ф 0,1, \р\ < 2к, р е Ц.

В (1.2.114) положение вектора £ на окружности, задаваемой уравнением £2 + 2к£ = 0, характеризуется углом ф1 е [-——) [101]:

£ = Rek(cosф1 -1) + Иек1 Бтф, где к1 = 1тк х Иек/^тЦ. (1.2.115)

Значения ; и г2 определяются исходя из выражений [70]:

гДЛрф)

к \в1(-£) + 1о1(-£)]

к - £ /4

72 (Л,р,ф)

(к + £)-[^(р + £) + шх(р + £)]

(1.2.116)

(1.2.117)

4к2 - |р + £2 / 4 '

Вектора к,£,6^(-£), ¿у1(-£), 01(р + £), ¿уДр + £) находятся из соотношений (1.2.106), (1.2.115), (1.2.104) соответственно.

Следует отметить, что нахождение {Н, Н)(Х,р) остается строгим до тех пор,

пока Н(к, р) = 0 при |р| > 2к . Это условие является ограничением на область

локализации пространственного спектра рассеивателя и его вторичных источников. Однако такое ограничение не сужает существенным образом область практического применения алгоритма. Заключительным этапом является оценка пространственного спектра рассеивателя из соотношения (1.2.103) принимающего в терминах переменных Х -плоскости вид

£(-р) = Цш Н (Х,р). (1.2.118)

Затем определяется само распределение функции неоднородности е(т) с помощью преобразования Фурье.

По сравнению с алгоритмом Новикова - Хенкина данный алгоритм существенно более удобен в вычислительном плане и не содержит неустойчивых операций [70]. Главная трудность алгоритма Новикова заключается в использовании нелинейного интегрального соотношения (1.2.114), что говорит о том, что необходимо применить итерационные методы нахождения Н(Х,р). При численной реализации вводится также дополнительная фильтрующая функция ^(|р|), при использовании которой схема итерационного решения выглядит следующим образом [70,71]

Н (Х,р) = ^ (|р|) Н 0 (Х,р) + ^ (|р|) Р ^Н,^ Н 7--Х|;Х,р ^, (1.2.119)

гдеН;=0(Х,р) = ^(|р|)Н0(Х,р), Н = (1 -ех)Н— + £1Н], 0<е1 < 1- весовой множитель, Н] - оценка функции Н (Х,р) получаемая на итерации с номером

] > 1 [70,71]. Весовой множитель использовался для избежания возможной «раскачки» итерационного процесса в случае достаточно сильных рассеивателей. Введение дополнительной фильтрующей функции ^(|р|) улучшает сходимость

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

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

полученной пространственной фильтрацией точного спектра функцией ^ (|р|). Для

более сильных рассеивателей могут возникать осцилляции, которые объясняются тем, что шарообразные рассеиватели создают рассеяние назад и имеют компоненты в пространственном спектре при |р| > 2к [70]. Отмеченная

особенность является своего рода помехой для алгоритма Новикова но, тем не менее, не выявлено жестких ограничений на силу рассеивателя, необходимых для обеспечения работоспособности алгоритма, все объекты были восстановлены с хорошим качеством. Помехоустойчивость алгоритма по результатам испытаний с присутствием шумов достаточно высока и пригодна для практических целей диагностики [70]. Относительная среднеквадратическая погрешность восстановления составляла 8«0.036 и 8«0.11 при среднеквадратических ошибках в данных 8В = 0.03 и 8В = 0.09 соответственно.

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

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

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

(Г)) = 1 £ J J(gm (r, t) - gbm (r, t)) drdt, (1.3.1)

2 m=0 0 Sm

где s2- множество искомых распределений неоднородностей, S - поверхность измерения, T - отрезок времени, в течение которого проводится измерение граничных полей акустического давления gbm (r, t) в m -ой конфигурации, gm (r, t)

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

^(r) = *2Л (r) Ф'£(г2А (r)), (1.3.2)

где sln (r)- приближений распределений неоднородности на n -ой итерации, Ф^ (^2 n (r)) - градиент функционала (1.3.1), тп -шаг спуска вычисляемый из условия

минимума (1.3.2) с помощью одномерной оптимизации.

В различных работах градиент функционала (1.3.1) находится по-разному. Так в работах А.В.Гончарского, С.Ю.Романова, С.Ю. Сережникова, С.А. Харченко [78-82] составляющие градиента функционала (1.3.1) для скорости звука c(r) и коэффициента поглощения a(r) для m-ой конфигурации одного точечного источника и приемников имеют вид

, , чч гdw (r,t) dg (r,t) 7 „ „ чч Tdw (r,t) Л ч 7 „ „ Ф'c.m(c(r)) = J ^ ; dt, Фa,m (a(r)) = J ^' ; Agm (r, t)dt , (1.3.3)

где gm (r, t) - решение прямой задачи

С2 (r) + a(r) d (Agm ) Ag m (r, t) = S(r - r„) f (t), (1.3.4)

ar dt

(r, t)

gm t),

dt

0, d ng|sT = Pm (r, t). (1.3.5)

t=0

0

где gm (r, t)-поле акустического давления в области q4 с , N = 2,3, окружающей неоднородность. Область Q4 окружена поверхностью Sm, r1 -радиус-вектор точки расположения источника, dng\ - производная функции g

П SmT

вдоль нормали к поверхности Sm в области Sm х (0, T), p(r, t) - известная функция. Предполагалось также, что вне области неоднородности скорость звука равна c(r) = c0 = const.

Функция wm (r, t) есть решение сопряженной начально-краевой задачи

С 2 (r) + 4*(r)dw^ 1 - AWm (r, t) = 0, (1.3.6)

dt [ dt J

dw (r, t)

m ^ " '

w (r, t) =

m 4 ' ylt=T

= 0, д w I = g(r,t)| - gh (r,t). (1.3.7)

" n m\S t^ IS T ° h,m 4 " / v ^

t=T

dt

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

Решение прямой и сопряженной задач осуществлялось методом конечных разностей с явными разностными схемами. Данный численный метод относится к SIMD - алгоритмам (Single Instruction - Multiple Data) и имеет высокую степень параллелизма. Такие задачи могут быть эффективно решены на кластерных системах, как с процессорами общего назначения, так и с графическими ускорителями (GPU) [80] . Расчеты для работ [80,82] в которых приведены данные вычислительных экспериментов были проведены на вычислительном комплексе «Ломоносов» суперкомпьютерного комплекса МГУ имени М.В. Ломоносова.

В работах [76,77] выражения для градиента функционала и начально-краевых задач другие. Так в [76] компонента градиента для скорости звука c(r) имеет вид

1 м IT d2 g (r t)

ФС(c(r)) = — wm(r,T-1) gdm(2 , )dt, (1.3.8)

С (r) m=0 0 dt

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

ду(г, г)

дг

dp(r, t)

= -Vg (r, t), (1.3.9)

дг

= -V • v(r,t) + 4л f s(r,t')dt', (1.3.10)

g (r, t) = с2 (r)

p(r, t) + z(r) |-(- Ap(r, t ))y/2-1 + r(r)(- Ap(r, t)) dt

( y+1) / 2-1

(1.3.11)

где g (r, t), v(r, t), p(r, t) - акустические поля давления колебательной скорости и плотности соответственно, s(r, t) - функция источника. Функции г(г) и r(r) соответствуют акустическому поглощению и дисперсии при распространении звука.

Во втором случае решалась начально-краевая задача для неоднородного волнового уравнения

1 д2

Ag(r,t) - -—■— g(r,t) = -4лs(r,t), (1.3.12)

с (r) дt

начальные и граничные условия приведены не были. Аналогично решалась сопряженная начально-краевая задача для уравнения

1 д2

Aw(r, t) —w(r, t) = -^(r, t), (1.3.13)

с (r) дt

где rx(r,t) = g(r,T -1) - gb (r,T -1). Все задачи решались с помощью псевдоспектрального метода к -пространств [84,85]. Стоит отметить, что в выражениях (1.3.9)-(1.3.13) не проставлен индекс m, вследствие того что в работе [76] широко применялся метод кодирования источника (source encoding method). Данный метод подробно освещается в работе [83] и заключается в следующем. Функционал (1.3.1) заменяется на следующий

ф(^(г)) = 1 £ff(em (t) *(g (s2(r), sm (r, t))-gb,m (r, t)))2 drdt, (1.3.14)

2 m=0 0 S

0

где ет - кодирующая последовательность, а знак * означает свертку. Так как решение 8т (х, г) линейно относительно функции источника ^ (х, г) при т-ой конфигурации, то функционал (1.3.14) можно переписать

1 Т (( ( м Л м Л Л2

Ф(*2 (г)) = 211(I 8 ^ (г), Е ет(г) * ^(г, г) ] - ; ет(г)8Ьт(г, г) йЫг. (1.3.15)

\ т=0 У т=0 ] ]

Фактически это означает, что вместо решения прямой и сопряженной задач для каждой конфигурации источников и приемников в отдельности можно решить

м

только одну задачу для сложного источника Е(г, г) = ; ет (г) * (г, г) и граничных

т=0

М

данных 8ы.=; ет (г) * 8Ьт (г, г), что позволяет значительно сократить

т=0

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

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

Ях(с(т))= |

^дс(г) + дс(т)Л

дх

ду

йхйу, Я2 (с(г)) = |

У

8 +