Содержание к диссертации
Введение
2. Многочастотное обобщение функционального метода решения обратной двумерной монохроматической задачи рассеяния 24
2.1. Формализм комплексных волновых векторов. Основные понятия, уравнения и обозначения функциональных методов решения обратных задач 24
2.2. Модифицированный двумерный алгоритм Новикова в монохроматическом и полихроматическом режимах 31
2.3. Алгебраизация уравнений при численной реализации алгоритма... 51
2.4. Решение прямой задачи рассеяния для цилиндрических рассеивателей с произвольным показателем преломления 55
2.5. Численное моделирование двумерного модифицированного алгоритма
2.6. Выводы 72
3. Трехмерная обратная акустическая задача рассеяния - алгоритм новикова-хенкина 74
3.1. Основные обозначения и уравнения алгоритма 74
3.2. Алгебраизация уравнений при сборе экспериментальных данных и численной реализации алгоритма 89
3.3. Результаты численного моделирования алгоритма 97
3.4. Выводы 105
4. Трехмерная обратная акустическая задача рассеяния - модифицированный алгоритм новикова 107
4.1. Основные обозначения и уравнения алгоритма 107
4.2. Численное моделирование алгоритма 119
4.3. Выводы 142
5. Основные результаты и выводы 144
Литература 146
- Модифицированный двумерный алгоритм Новикова в монохроматическом и полихроматическом режимах
- Численное моделирование двумерного модифицированного алгоритма
- Алгебраизация уравнений при сборе экспериментальных данных и численной реализации алгоритма
- Численное моделирование алгоритма
Введение к работе
Работа посвящена моделированию функционально-аналитических методов решения- двумерных и трехмерных обратных задач акустического рассеяния. В настоящее время такие методы представляются более мощными и перспективными, чем итерационные; изначально они были ориентированы на решение обратных задач квантового рассеяния, но в последние годы появились общие функциональные методы, пригодные для решения обратных задач рассеяния различной физической природы (в том числе акустической). Анализ возможности применения» подобных методов на практике очень, важен. В представляемой диссертационной работе исследуется возможность, расширения уже хорошо зарекомендовавшего себя двумерного алгоритма Новикова-Рриневича на немонохроматический случай, а также ранее не исследовавшиеся, с точки зрения акустических приложений, трехмерные алгоритм Новикова-Хенкина и новый алгоритм Новикова.
Актуальность темы и цели исследования. Структура диссертационной работы Актуальность темы
Решение обратных волновых задач представляет собой активно развивающееся- направление- как математической, так и прикладной- физики (акустики и оптики). В акустике под обратными волновыми задачами понимается восстановление источников звука (обратная задача излучения) или характеристик неоднородностей, рассеивающих первичное поле (обратная? задача рассеяния); по измерениям акустического поля. Наиболее актуальными направлениями применения на практике результатов разработок в теории обратных задач рассеяния являются медицинская диагностика и акустика океана. Проблемам разработки медицинских акустических томографов, решающих задачу ранней диагностики, более безопасных, чем рентгеновские, и менее дорогостоящих, чем ЯМР-томографы, посвящено в последнее время большое количество как теоретических, так и экспериментальных исследований. Помимо медицинских приложений, результаты исследований различных обратных задач рассеяния имеют широкое применение в задачах создания систем промышленной диагностики (дефектоскопии), решения прикладных проблем геоакустики и акустики океана.
Существует два крупных класса подходов к решению обратной задачи-рассеяния один из которых основан на итерационных методах, а другой - на методах функционального анализа. Преимуществом итерационных методов является то, что они могут базироваться на фрагментарных данных, полученных при различной геометрии эксперимента, различных частотах. В этих методах не накладывается жестких требований на полноту данных рассеяния отдельно для, каждой из частот или каждой конфигурации падающего поля в эксперименте, любая априорная и апостериорная информация может быть также использована. Основные же преимущества функционального подхода заключаются как в получении строгого (или почти строгого) решения, так и (в ряде случаев) в меньшем количестве вычислительных операций, по сравнению с итерационными методами.
Функционально-аналитические методы имеют свои корни в математических и физико-теоретических работах, исследовавших процессы квантовой" теории рассеяния и, соответственно, решения прямых и обратных задач потенциального рассеяния. В силу ряда принципиальных отличий природы потенциального рассеяния в квантовой механике и рассеяния волн на неоднородностях фазовой скорости, такие методы не всегда» применимы в полной мере к обратным задачам акустики. Однако в последние годы наблюдается стремление разработать общие подходы, пригодные для решения обратных задач рассеяния различной физической природы (как квантово-механических, так и акустических, электродинамических, эластодинамических), в том числе изоэнергетических (монохроматических) обратных задач.-Проведенные исследования опираются на методы функционального и многомерного комплексного анализа1, и отличаются высокой математической строгостью.
Анализ возможности применения на практике разработанных в последние десятилетия функционально-аналитических методов решения двумерных обратных задач рассеяния был подробно проведен в последнее время в работах сотрудников кафедры акустики физического факультета МГУ. Созданы работоспособные компьютерные программы, позволяющие достаточно быстро восстанавливать двумерные рассеиватели любой формы по данным рассеяния на основе функционально-аналитического двумерного монохроматического алгоритма Новикова-Гриневича [1,2,3], и показано, что область работоспособности этого алгоритма реально намного больше, чем первоначально предполагалось [4]. Тем не менее, при увеличении силы, рассеяния, на модельных численных решениях двумерных задач наблюдалось возникновение неустойчивости и повышенной чувствительности к ошибкам в данных рассеяния [4]. Основная и принципиальная причина этого явления, кроется в размерностной безызбыточности данных рассеяния для двумерной монохроматической задачи. А именно, массив таких данных соответствует двумерной параметрической области углов падения и рассеяния (0,27г)х (0,2тг). Той- же размерностью характеризуется искомая двумерная неоднородность -конечных размеров: Отсутствие- размерностной избыточности данных для монохроматической задачи приводит к тому, что условия многочисленных теорем и утверждений относительно единственности решения таких задач (одна из ранних работ в этом направлении - [5]) носят ограничительный характер [6, 7].
Импульсный, или, во многом-эквивалентный ему, многочастотный режим снимает это ограничение (поскольку в параметризации данных рассеяния появляется дополнительное измерение - частота). Однако простое аддитивное объединение множества решений монохроматических задач не является эффективным методом и не эквивалентно решению единой задачи, в котором, используется факт общности рассеивателя для каждой монохроматической задачи. Поэтому необходим более полный метод объединения решений, который не лежит «на поверхности» и нуждается в дополнительном исследовании. Обратные-задачи рассеяния в трех измерениях обладают размерностной избыточностью даже в монохроматическом случае. Поэтому проблемы обеспечения единственности, свойственные двумерной задаче, здесь не возникают [6, 7, 8]. В настоящее время основным методом синтеза трехмерных акустических томограмм является простое объединение послойных двумерных изображений. В этом методе пренебрегается «обменом» рассеянными полями между слоями. Между тем, ошибки, вызываемые этим пренебрежением, имеют порядок второго борновского члена (при разложении рассеянного поля в ряд Борна-Неймана), т.е. достаточно существенны.
Строгое решение трехмерной обратной задачи рассеяния функциональными методами не удается получить простой модификацией двумерных алгоритмов, что требует разработки новых методов решения. Результаты этих исследований появились в самое последнее время, но их практическая пригодность в акустических системах отнюдь не очевидна, а требуемый объем вычислений в настоящее время представляется чрезвычайно большим. Тем не менее, быстрый прогресс в данной области позволяет поставить вопрос о первых попытках реализации этих новых подходов при решении трехмерных обратных задач акустического рассеяния (хотя бы самых простых). Любой способ решения обратных задач рассеяния на практике требует детального исследования алгоритмов восстановления рассеивателеи с целью их наилучшего функционального и технического согласования с измерительной установкой (томографом), методом получения первичных данных рассеяния и способами отображения итоговых результатов. При этом с одной стороны, функциональные методы являются наиболее продвинутыми и мощными в теоретическом и принципиальном плане, а с другой стороны, их применение для решения обратных задач рассеяния классических полей находится только в начальной стадии исследования, что делает несомненно важным дальнейшее продвижение в этом направлении. Поэтому актуальность представляемой работы заключается в создании метода решения единой акустической многочастотной двумерной задачи томографического типа, основанной на функциональном подходе. Это открыло возможность применения данного метода для импульсного режима медицинских томографов, повысив, тем самым, информативность и расширив область их работоспособности. Вторая часть работы носит исключительно пионерский характер и открывает цикл исследований по применению в решении трехмерных обратных задач акустического рассеяния строгих функциональных методов их решения. Эти исследования направлены, в конечном счете, на создание трехмерных систем акустоскопии принципиально нового типа.
В связи с этим, можно выделить следующие основные цели диссертационной работы.
• Найти и апробировать на численных моделях метод органичного объединения многочастотных данных и методов решения множества монохроматических двумерных обратных задач акустического рассеяния в виде единого функционального алгоритма, базирующегося1 на монохроматических вариантах решения, развитых в работах С.П.Новикова, Л.Д.Фаддеева, П.Г.Гриневича, С.В.Манакова и включающих последние результаты Р.Г.Новикова.
• Реализовать и апробировать на простейших моделях функциональный метод решения обратной трехмерной задачи акустического рассеяния, основываясь на результатах исследований и алгоритмах Г.М.Хенкина и Р.Г.Новикова.
• Провести- сравнительную оценку практической области работоспособности указанных подходов.
В связи с двумя основнымиг направлениями работы основные задачи диссертационной работы можно разбить на две группы.
По первой части работы (немонохроматическому двумерному функциональному алгоритму):
1. Найти метод объединения последовательности операций при решении монохроматических частных задач в единый взаимосвязанный процесс.
2. Продемонстрировать расширение области работоспособности немонохроматического двумерного алгоритма решения обратной задачи акустического рассеяния по сравнению с простой суммой монохроматических решений.
3. Оценить вычислительную сложность алгоритмической реализации этого метода и физические ограничения на область его применимости.
По второй части работы (трехмерным функциональным алгоритмам) задачи соответствуют начальному этапу исследований:
4. Реализовать трехмерные функциональные алгоритмы в виде конкретных работающих программ и исследовать с их помощью восстановление характерных рассеивателей простейшей формы.
5. Оценить возможности и перспективы реализации этих алгоритмов для практических целей.
Научная новизна работы:
Г. Поставлена не рассматривавшаяся ранее задача поиска методов объединения невзаимосвязанных решений множества монохроматических задач акустического рассеяния в единый процесс согласованного использования всей совокупности данных рассеяния. Найдено, обобщение двумерного метода Новикова-Гриневича» и модифицированного метода Новикова на немонохроматический случай.
2. Проведен цикл модельных исследований возможностей предложенного метода решения полихроматической обратной задачи томографического типа. Полученные результаты позволяют прийти к выводу о перспективности использования модифицированного полихроматического двумерного алгоритма Новикова для прикладных задач акустического томографирования медицинской направленности. Полихроматический алгоритм более устойчив и информативен, чем результат аддитивного синтеза одночастотных решений при тех же исходных данных.
3. Впервые выполнена конкретная реализация математических алгоритмов решения трехмерной обратной задачи акустического рассеяния и проведено исследование возможностей этих методов на примерах различных модельных задач на основе алгоритмов Новиоква-Хенкина и нового алгоритма Новикова.
4. Проведен цикл численных модельных экспериментов, продемонстрировавший широкие возможности нового функционально -10-аналитического трехмерного алгоритма Новикова и выявивший его работоспособность для рассеивателеи произвольной силы.
Обоснование и достоверность полученных результатов: Достоверность результатов,- представленных в диссертации, подтверждается решением обратных модельных задач, давших оценки, близкие к исходным характеристикам двумерных и трехмерных рассеивателеи, использованным при синтезе тестовых данных рассеяния.
Практическая ценность работы
1. Показана практическая реализуемость и. широкие прикладные, возможности предложенного обобщения двумерного функционально-аналитического метода на немонохроматический случай, что открывает возможность его применения В реальных ультразвуковых томографах, работающих в импульсном или многочастотномрежиме.
2. Показана практическая реализуемость алгоритма- решения трехмерной монохроматической обратной задачи, а также возросшая устойчивость решения? по сравнению с двумерной задачей, что позволяет, говорить о практической перспективности данного направления, требующего, однако, высокой производительности используемых вычислительных средств.
3. На численных примерах исследована помехоустойчивость перечисленных алгоритмов, которая оказалась достаточно высокой для практических целей медицинской диагностики.
Личный вклад автора заключается» в проведении- физического анализа основных методов решения обратной задачи акустического рассеяния в двумерном и трехмерном пространстве, позволившего найти органичное объединение монохроматических методов в единый процесс нахождения полихроматического решения, в анализе физического смысла операций многошагового процесса решения трехмерной обратной задачи и в разработке конкретных численных схем всех обсуждаемых в работе подходов и их конкретной реализации. Большинство работ по моделированию и анализу полученных результатов проведены им лично.
Основные положения, выносимые на защиту:
1. Совокупность отдельных (независимых) процедур решения монохроматических обратных задач акустического рассеяния может быть объединена предложенным в диссертации образом во взаимосвязанный единый процесс поиска решения, приспособленного для применения в акустических томографических системах, работающих в импульсном режиме. Результат такого объединения - возросшая информативность и помехоустойчивость томограмм.
2. Трехмерная,обратная задача рассеяния имеет практически реализуемый и промоделированный в диссертации путь решения, учитывающего всю сложность многократного рассеяния сильными неоднородностями, пригодный к практическому применению в системах акустоскопии, снабженных многоэлементной приемно-излучающей системой и высокопроизводительными вычислительными устройствами.
Апробация работы
Материалы диссертации докладывались на конференции «Ломоносов-2004» (Москва, апрель 2004), XV сессии Российского Акустического Общества (Нижний Новгород, ноябрь 2004), конференции «Ломоносов-2007» (Москва, апрель 2007), XIX сессии Российского Акустического Общества (Нижний Новгород, сентябрь 2007) и семинарах кафедры акустики физического факультета МРУ.
Публикации: основные результаты диссертации изложены в семи работах [102-108] (две из них - в рецензируемых журналах), приведенных в списке литературы.
Структура и объем диссертации
Диссертация состоит из пяти глав объемом 155 страниц, включающих 130 страниц текста и 23 рисунка, а также списка цитируемой литературы.
Модифицированный двумерный алгоритм Новикова в монохроматическом и полихроматическом режимах
В этом и последующих разделах главы 2 осуществляется и анализируется расширение модифицированного двумерного алгоритма Новикова на полихроматический случай. В настоящее время; единственным- алгоритмом решения; обратных двумерных; волновых акустических- монохроматических задач рассеяния средствами функционального- анализа, доведенным до стадии численного моделирования, является алгоритм Новикова-Гриневича. Он обеспечивает строгое решение задачи; с учетом эффектов многократного рассеяния, однако накладывает ограничения (которые будут обсуждены далее) на класс восстанавливаемых рассеивателей. Впоследствии алгоритм Новикова-Гриневича был модифицирован Р.Г.Новиковым [79, 88]. Этот новый вариант двумерного монохроматического алгоритма, обсуждению преимуществ которого и обобщению на полихроматический случай посвящена данная глава, будет далее называться модифицированным двумерным алгоритмом Новикова. Надо отметить, что в главе 4 обсуждается модифицированный трехмерный алгоритм Новикова, имеющий достаточно мало общих черт с модифицированным двумерным алгоритмом Новикова. Рассматриваемая модификация двумерного алгоритма приведена в. упоминавшихся работах его автора в терминах комплексных- переменных. В представляемой же работе все этапы алгоритма описываются.в терминах классических углов и координат. Еще раз надо подчеркнуть, что соотношение (2.2.6) справедливо при условии отсутствия рассеяния назад обобщенных волн. Это то же самое ограничение, какое имело место и для алгоритма Новикова-Гриневича.
Соотношения (2.2.19), (2.2.20) как и (2.2.6), справедливы при условии отсутствия рассеяния назад обобщенных волн. -39-Итак, последовательность реализации модифицированного алгоритма в монохроматическом варианте, т.е. при фиксированном со , такова. Из данных /((р,ф ;сз7) находятся /г±(ф,ф ; соу) решением линейных систем (2.2.3) в пространстве углов при каждом фиксированном ср. Из / (ф 1; со.) строятся функции / (!,ф,ф ; со7) и (г ф ; со )» согласно (2.2.5). Нахождение классического поля из линейной системы (2.2.6) равносильно нахождению угловых гармоник этого поля р!с1(г,#; со7) из системы (2.2.14) с учетом выражения (2.2.15) для В, в котором Q (r,q, qf; со ) рассчитываются применением к О (г,ф,ф ; со ) Фурье-преобразования по обоим углам. Далее предельные значения обобщенного поля непосредственно рассчитываются на основе значений классических полей из (2.2.16) или, что удобнее, ,(— 1)-ая угловая гармоника обобщенного поля \T{r,q = — 1; со,) рассчитывается из (2.2.18). Функция рассеивателя находится из (2.2.19) или (2.2.20), где операция дифференцирования может быть выполнена с помощью пространственного Фурье-преобразования, однако лучший результат дает метод конечных разностей. Модифицированный алгоритм Новикова сохраняет все достоинства алгоритма Новикова-Гриневича: локальность по пространственной координате (функция рассеивателя восстанавливается в фиксированной точке г независимо от других точек), возможность получения решения безытерационным способом (поскольку решаемые системы линейны относительно неизвестных), учет эффектов перерассеяния. В монохроматическом режиме качество восстановления рассеивателей данным алгоритмом и алгоритмом Новикова-Гриневича эквивалентны. Оба монохроматических алгоритма имеют одни- и те же пределы работоспособности, с точкич зрения ограничения класса восстанавливаемых рассеивателей.
Численное моделирование двумерного модифицированного алгоритма
Во всех рассматриваемых далее примерах с цилиндрическими рассеивателями радиуса а, показатель преломления рассеивателя полагался равным w = с0/с = 1.1 (с - скорость звука внутри цилиндра), т.е. рассеиватель задавался фокусирующим. Эффекты перерассеяния поля внутри фокусирующего рассеивателя (и 1) выражены значительно сильнее, по сравнению с дефокусирующим рассеивателем (п 1), при одной и той же абсолютной величине относительного контраста Ас /с0. Поэтому именно фокусирующие рассеиватели были выбраны для иллюстрации возможного эффекта возникновения неустойчивости решения обратной задачи в монохроматическом режиме и стабилизации решения в полихроматическом режиме за счет объединения одночастотных решений через условия связи. Радиус всех цилиндров составлял а = ЗА-о , где -60 A,3 = 8 е.д.д. — длина волны, соответствующая самой низкой частоте. Тогда A\\f=l.2n при ,(,.= ,0 , и А\[/ увеличивается пропорционально росту частоты. Такой достаточно большой набег фазы в сочетании с фокусирующим характером рассеивателя приводит к сильно выраженным эффектам перерассеяния.
Для рассеивателя с поглощением показатель преломления полагался равным п + iqa. Мнимая добавка qa задавалась независящей от частоты со , что соответствовало линейной зависимости от со, коэффициента поглощения внутри рассеивателя a(r,coy) =k0jqa. Это нисколько не ограничивает общности получаемых далее результатов, но удобно при численном моделировании, поскольку Imv(r,co ) приобретает квадратичную частотную зависимость Шу, так же, как и Rev(r,co ) (в пренебрежении малой поправкой а2(г,со )- [103]). Тогда в (2.2.1) v(r,co ,)/оОу не зависит от частоты, что позволяет использовать уравнения связи (2.2.28), приводящие к (2.2.30). В то же время, такой вид частотной зависимости не редок в биологических тканях [87, 91].
В последующих этапах алгоритма везде участвовала функция / (фэф ; со,)JP( ф — ф ) вместо (ф 1; со ). Численное моделирование показало, что фильтрация /г±(ф,ф ; со ) значительно эффективнее фильтрации /(ф,ф ; Шу). Это можно объяснить тем, что /г±(ф,ф ; со ) находится на основе произвольного /(ф,ф ; со ) из строгих уравнений Фаддеева (212.3), а далее классическое волновое поле находится на основе /г±(ф,ф ; со ) из соотношений, предполагающих отсутствие рассеяния назад обобщенных полей. При фильтрации именно /2±(ф,ф ; соу) влияние этого обратного рассеяния устраняется в большей степени. Оценки v(r,co )/&о , полученные модифицированным алгоритмом в монохроматическом режиме при соответствующей длине волны X0j, приведены на рис.За. Помехи в данных рассеяния /(ф,ф ; соу) отсутствовали, однако качество восстановления в выбранном диапазоне длин волн оказалось неудовлетворительным: оценки искажены сильными флуктуациями и даже резкими выбросами при А,0-= 7.5е.д.д. (Из-за выбросов, оценка при Х0 = 7.5 е.д.д. изображена в диапазоне вертикальной шкалы, отличном от диапазона других рисунков.) Уменьшение шага дискретизации по углам ф и ф не изменяет характера оценок, и, следовательно, искажения вызваны именно неустойчивостью монохроматического решения из-за большой силы рассеивателя. Как показало моделирование различных рассеивателей, возникновение неустойчивости монохроматической оценки в конкретной точке г сопровождается уменьшением обратной обусловленности матрицы + 27і-5(г,#,—g ; ce j) (Е - единичная матрица) при неизвестных угловых гармониках классического поля (Iе1 в системе (2.2.14) до = (3-т-5)х10-3.
Под обратной обусловленностью понимается отношение минимального собственного значения матрицы к максимальному.) Для рефракционно-поглощающего рассеивателя неустойчивость в рассматриваемом диапазоне длин волн возникает не только вблизи границы рассеивателя, но и в глубине него (рис.За).
Весовой коэффициент Р для уравнений связи (2.2.30), добавляемых к системе (2.2.14) в полихроматическом режиме, должен быть настолько большим, чтобы итоговая оценка wo]y(r) (2.2.31) была одинаковой при всех XQ-.
Дальнейшее увеличение 3 не изменяет этой оценки. Далее везде полагалось Р = 5; регуляризация МНК-решения системы не требовалась. Тем самым, устойчивость, монохроматического решения не всегда монотонно возрастает при снижении рабочей частоты. Это объясняется тем, что за счет сложного характера эффектов перерассеяния, эффективная сила рассеяния может как уменьшаться, так и увеличиваться (хотя амплитуда функции рассеивателя (2.2.1) уменьшается). Качество полихроматического решения w 1у, использующего рассматриваемые девять длин волн, достаточно хорошее (рис.4б, слева). Среднеарифметическое vaver (г) набора соответствующих монохроматических решений (рис.4б, справа) дает более сильные осцилляционные погрешности в области границы, хотяв данном случае различие между wpoly и vaver выражено не так сильно, как для предыдущего рассеивателя, поскольку четыре монохроматических решения из девяти близки к истине. В еще более длинноволновом диапазоне XQj- от 7 до 8 е.д.д. монохроматические оценки становятся устойчивыми и в отсутствие помех воспроизводят рассеиватель практически идеально (рис.5а). Интересно отметить, что для рассеивателя с тем же контрастом скорости, но обладающим существенным поглощением, восстановление в том же длинноволновом диапазоне уже становится неустойчивым, приводя к сильным искажениям воспроизведенных значений характеристик рассеивателя (рис.За).
Алгебраизация уравнений при сборе экспериментальных данных и численной реализации алгоритма
В данном параграфе приводится алгебраизация (дискретизация) используемых в алгоритме уравнений, необходимая для их компьютерного моделирования. Функции Лежандра POT(cos6) в (3.1.7) вычисляются так: специальной операцией метематической программной среды вычисляются полиномы Гегенбауэра (они же ассоциированные функции Лежандра) r,„(cos0n), после чего у получившегося столбца берется первое значение - это и есть функция Лежандра: Рт (z) = Т% (z).
Собрав «экспериментальные» данные, можно приступать к решению обратной задачи - численной реализации алгоритма Новикова-Хенкина: Первым (а в приближении сферической симметрии - единственным) этапом алгоритма является получение обобщенной амплитуды рассеяния; из классической путем решения уравнения Фаддеева. Рассмотрим дискретизацию уравнения Фаддеева (3:1:6) сначала для; общего случая, а затем в приближении сферической-симметрии.
Далее, с учетом анализа, проведенного в предыдущем параграфе, рассмотрим вид уравнения Фадцеева в дискретном виде с учетом сферической симметрии рассеивателя. В наших обозначениях угол а - это угол между к и 1 (он же угол 0;). Угол между m и Г обозначим как а .
Непосредственно в программной реализации моделирования решения системы уравнений Фаддеева сначала рассчитывается четырехмерный массив Ґ{пп- «ф,,п «,п ) (в процессе чего вычисляются углы а1 и значения функции -96-Хевисайда). Далее вычисляется матрица А$(п ,,п ;пр«,п ), а затем из (3.2.1) находятся искомые значения И$(п ,,иф). Наконец, с учетом, что Г=1 при 9/.=9/=0, получается нужная для дальнейших этапов нахождения характеристикрассеивателя функция h+ (п ) = h$ {пр =np,nVf = 0). Следует отметить, что геометрия задачи позволяет существенно уменьшить число неизвестных в решаемой системе уравнений. Из-за присутствия функции Хевисайда под интегралом в интегральном уравнении Фаддеева (2.1.9) присутствуют не все возможные т, а только их половина. Тогда при формировании системы относительно функции h + (к, V) следует рассматривать не все Iі, а только их соответствующую половину (поскольку при решении системы перебираются все векторы V, принимающие те же ориентации, что и векторы m под интегралом). Таким образом, количество неизвестных уменьшается в 2 раза. Это общее положение вне зависимости от вида рассеивателя (т.е. симметрия рассеивателя не играет роли).
Следовательно, в рассматриваемом случае фиксированных векторов к и у+ (ориентированных, как оговорено выше), векторы m и V пробегают только значения с координатами (pw є (0, ти) и ф/, є (0, тг) соответственно. При этом Пцт r=0,N(?/2-l; четырехмерные массивы /+ {пр,, пщ , пр«, п ) и AsinP "Ф.Jпр« "Фш) имеютразмерность NpxN ?/2xNpxN /2. Для центрально симметричных рассеивателей существует дополнительная возможность уменьшения общего количества неизвестных еще в 2 раза, поскольку для таких рассеивателей функция h + (к, Iі) одинакова для тех двух различных ориентации вектора V, которые по отношению друг к другу являются зеркально отраженными относительно плоскости векторов к и у+ (это видно непосредственно из уравнения Фаддеева, т.к. /(т,Г) зависит только от угла а между m и Г, и под интегралом осуществляется перебор т): / (ир.,иф,) = /г(ир.,-Л ф/2-1-иф/1).
На этапе компьютерного моделирования член Zx (поверхностный интеграл) не рассматривался, что существенно упрощало численную реализацию алгоритма и, одновременно, обеспечивало устойчивость решения. Однако, за счет данного приближения, процессы многократного рассеяния учитываются при восстановлении не в полной мере, и решение - оценка пространственного спектра v(-p) и соответствующая ему функция рассеивателя v(r) - утрачивает строгость. Полагалось Imk -»0, и, следовательно, в реализованном приближенном варианте алгоритма спектр v(-p) оценивался на основе предельных значений Н обобщенной амплитуды рассеяния только для пространственных компонент р 2к0: у(-р)«Я (к.,р). (3.3.1) Поэтому результат восстановления v(r) сравнивался с функцией vcut(r), пространственный спектр которой совпадает со спектром истинного рассеивателя v(r) внутри сферы радиуса 2к0 и равен нулю вне этой сферы.
Это четко видно на примере чисто рефракционных рассеивателей, описываемых действительной функцией v(r). В случае слабых (борновских) рассеивателей, для которых процессы перерассеяния пренебрежимо малы, обобщенная амплитуда рассеяния совпадает с классической, и обе оценки совпадают между собой v(r) = vbom (г), отличаясь от истинной функции v(r) только за счет обнуленных при восстановлении высокочастотных компонент v(-p) при р 2к0. Однако, как только рассеиватель перестает быть слабым, в борновской оценке возникает ложная мнимая часть Imvborn(r), возрастающая по мере увеличения силы рассеивателя.
С другой стороны, наблюдаемая погрешность восстановления рассеивателя определяется именно значением Ац/, а не размером рассеивателя и его контрастом в отдельности. Надо учесть, что в этих численных экспериментах помеха, вносимая во все данные, соответствующие одному и тому же значению угла рассеяния (угол между к и I), полагалась одинаковой. Следовательно, избыточность трехмерной задачи при восстановлении не использовалась, и такая коррелированная помеха увеличивала значение погрешности \х, по сравнению со случаем некоррелированной помехи. В то же время, при некоррелированных ошибках измерения, размерностная избыточность данных /(к,1) (к,ІЕЖ3) в трехмерной задаче приводит к повышению помехоустойчивости решения. Действительно, одному и тому же значению вектора р = к -1 соответствует континуум пар векторов k,l ER , различающихся по направлению. Так, конец вектора к может вращаться по окружности в плоскости, перпендикулярной вектору р.
Численное моделирование алгоритма
В случае рассеивателей произвольной формы или конфигурации решение прямой задачи может быть только приближенным [99, 100]. Поэтому работоспособность трехмерного алгоритма Новикова иллюстрируется на примере рассеивателей в виде шара радиуса а со скоростью звука с и коэффициентом поглощения или усиления а внутри шара. Шарообразный рассеиватель не является идеальным модельным примером для иллюстрации качества решения обратной задачи рассеяния, так как он имеет резкие границы. Выбор такого вида объекта связан только с тем, что для него существует строгое аналитическое решение прямой задачи, т.е. возможность вычисления /(к, Г) с помощью сходящегося ряда [96, тії, с.452-455]. Необходимо подчеркнуть, что алгоритм Новикова может быть применен для решения обратной задачи рассеяния от практически любых рассеивателей, и использование при численном моделировании сферически симметричных объектов не ограничивает общности получаемых результатов. При рассмотрении рассеивателей, обладающих как рефракционной компонентой, так и поглощающей или усиливающей, данные рассеяния моделировались на основе тех же аналитических выражений, что и для чисто рефракционного рассеивателя, но показатель преломления п = с0/с становился комплексной величиной n + iqa; тогда ct = &0 7a - амплитудный коэффициент поглощения в случае qa 0 или усиления в случае qa 0. Рефракционная компонента такого рассеивателя создает дополнительный набег фазы волны вдоль диаметра рассеивателя А\/= la \ kQ — k\=2ak0\\ — п\, где к = аз/с.
В случае сильных рассеивателей, создающих большой дополнительный набег фазы Avj/ тс или обладающих сильным поглощением (усилением), после первой итерации наблюдается сильное отклонение оценки v =1(г) от начальной пространственно-профильтрованной оценки v0(r). Однако при дальнейшем итерировании оценка сходится к итоговой практически без осцилляции. Как уже упоминалось, присутствие дополнительной фильтрующей функции F(\ р () заметно уменьшает количество необходимых итераций и ускоряет сходимость. Сама итерационная процедура оказалась устойчивой к выбору начального приближения для итерирования Iя й 1 нелинейного члена Wj-\ Нj-\ ] в (4-1.14). Рассмотрим алгебраизацию уравнений при численной реализации алгоритма. Экспериментальными данными рассеяния, получаемыми путем решения прямой задачи, является классическая амплитуда рассеяния /(к,1) (к,ІєІЕІ ), из которой на первом этапе алгоритма путем решения уравнения Фаддеева (4.1.1) находится обобщенная амплитуда рассеяния /zy±(k,l) (k,IeIR3). В разделах 3.1, 3.2 уже приводилось решение прямой задачи для случая рассеивателя в виде однородного шара, а также дискретизация уравнения Фаддеева - в общем случае и случае сферической симметрии рассеивателя. Поскольку к=1=[Г=А:0, 17 1=1 (смысл вектора I , использумого при решении уравнения Фаддеева, объяснен в разделе 3.1), то результат решения системы уравнений Фаддеева - предельные значения обобщенной амплитуды рассеяния - зависят, согласно (3.1.6), только от углов, характеризующих векторы к, 1 и у+: h(Q +,ф +;Qk,qk;0/.,ф,.); для вектора у" рассуждения аналогичны.
В обсуждаемом в данной главе модифицированном алгоритме Новикова обобщенная амплитуда рассеяния обозначается как Н± (k, р) = hy± (к, 1 = к — р). В сферически симметричном случае Н± (к, р) = Я ( р ) = h ( р ), и, кроме того, /г+(р) = /Г(р); поэтому речь идет о функции, которую можно обозначить просто как Hs(\ р ) = Я ( р ). В дискретном виде это по-прежнему одномерный массив Hs (пр ).
Хотя результат на каждой итерации j и итоговый результат зависят только от р в случае сферической симметрии, для реализации алгоритма необходимо определить конкретные координаты вектора р в пространстве. Используемый в алгоритме вектор v задается сонаправленным оси Oz, т.е. v= 1; 9V=0; (pv значения не имеет. Частично результаты численного моделирования приведены в [106]. Обсуждаемые ниже некоторые характерные примеры свидетельствуют, что итоговая итерационная оценка v(r), получаемая из найденного пространственного спектра v(% = -р) = Н(\ = 0, р), оказывается наиболее близкой к пространственно-профильтрованной функции истинного рассеивателя vcut (г). Для борновских рассеивателей v(r) совпадает с vcut(r). Для более сильных рассеивателей могут возникать отличия в виде осцилляции, которые объясняются тем, что рассматриваемые шарообразные рассеиватели создают рассеяние назад за счет резкого изменения акустических -126-характеристик на границе шара. Поэтому пространственные спектры рассеивателя и его вторичных источников содержат компоненты при Ь, 2kQ, являющиеся своего рода помехой для алгоритма Новикова и нарушающие, в случае неборновских рассеивателей, полное совпадение восстановленной оценки v(r) и ожидаемой в идеале функции vcut(r). Задавались скорость звука с0= 1500м/сек и длина волны А,0=3х10 3м, соответствующие реальной ситуации в ультразвуковой томографии мягких биологических тканей. На рисунках 15-23 изображены центральные сечения истинного рассеивателя v и его пространственно-профильтрованной версии vcut, а также итоговой итерационной оценки v, сравниваемой с фаддеевской vfad и борновской vbom. Параметры т, х{ фильтрующей функции F(\ р ) оговорены в тексте для каждого конкретного рассеивателя.
