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

ОБЪЕДИНЕННЫЙ ИНСТИТУТ ЯДЕРНЫХ ИССЛЕДОВАНИЙ

КОСАРЕВ Евгений Леонидович

ОБРАТНЫЕ ЗАДАЧИ ЭКСПЕРИМЕНТАЛЬНОЙ ФИЗИКИ И МЕТОДЫ ИХ РЕШЕНИЯ

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

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

10-94-411

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

УДК 621.391

Дубна 1994

Работа выполнена в Институте Физических Проблем им.П.Л.Капицы, Российской академии наук.

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

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

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

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

Ю.Е.Нестерихин Г.А.Ососкок

И. М. Соболь

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

Физический институт им.П.Н.Лебедева, Российской академии наук, г.Москва

Ур С.О

Защита состоится — » 1994 года в * часов

на заседании Специализированного ученого совета Д 047.01.04 при Лаборатории вычислительной техники и автоматизации Объединенного института ядерных исследований но адресу: 141980, г.Дубна, Московской обл.

С диссертацией можно ознакомиться в библиотеке Объединенного института ядерных исследований

Автореферат разослан ^_'' 1994 года.

Ученый секретарь Специализированного Совета кандидат физико-математических наук

З.М.Иванченко

Введение

Актуальность темы

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

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

Из сказанного ясно, что обе указанные задачи взаимосвязаны, и фактически именно шумы определяют величину возможного увеличения разрешения (или сверхразрешения) измерительного прибора. В отсутствие шумов разрешение возможно увеличивать неограниченно. Несмотря на то, что указанная точка зрения известна в физике ухе более 40 лет со времени работ Г.С.Горелика и И.Л.Бернштейна [2] и [3], практические алгоритмы восстановления были весьма несовершенны, что препятствовало их широкому использованию в экспериментальных исследованиях.

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

Af(y) = F(x), (1)

где А обозначает интегральный оператор (возможно и нелинейный)

ь

Af(y) = J K[x,y,f(y)}dy, (2)

а

определяемый конкретной задачей, f(y) - искомая функция a F(x) есть измеряемая на опыте зависимость. Целью является определение неизвестной функции f(y) по измеряемой F(x).

Однако, уже в первом издании монографии Куранта и Гильберта [4] было отмечено по отношению к интегральным уравнениям первого рода вида

/W = / K(s,t)<p(t)dt, (3)

что: " ... Затруднения для теории интегральных уравнении первого рода лежат в том обстоятельстве, что при непрерывном ядре K(s, t) многообразие всех кусочно-непрерывных функций if(s) преобразовывается в часть того же многообразия, так как

1Кавычками отмечена цитата из доклада В.А.Амбарцумяна [J] на международном симпозиуме по фундаментальным проблемам теоретической и математической физики в Дубне в 1979 г.

все получающиеся таким образом функции /(а) во всяком случае непрерывны. Если ядро Я(л,г) дифференцируемо, то всякая кусочно-непрерывная функция, даже всякая только интегрируемая функция преобразуется в дифференцируемую. Стало

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

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

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

Эта неустойчивость в свою очередь объясняется тем, что в бесконечномерном пространстве вполне непрерывные операторы, не могут иметь непрерывных обратных. К классу вполне непрерывных операторов относится оператор (3) при, например, непрерывном по обеим переменным ядре А"(з,(). Однако, как было показано в 1943 г. А.Н.Тихоновым [5] при ограничении области возможных решений компактным множеством М из взаимной однозначности отображения следует непрерывность обратного оператора.

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

Ах = у, х е Х,у 6 V,

а постулируется существование некоторого множества М С X в пространстве всех возможных решений X, такого, что решение единственно на этом множестве, т.е. на множестве N = А(М) определен обратный оператор А'1, и решение устойчиво в области Ы, т.е. обратный оператор А'1 непрерывен в соответствующей метрике, определенной в области N. Если М - компактное множество, то условие корректности по Тихонову следует из его теоремы 1943 г.

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

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

Многообразие и важность практических применений обусловило появление многочисленных обзоров и монографий по обратным задачам [6-22), включая также обзор автора [12]. Однако, несмотря на такое обилие литературы по данному предмету, безусловно свидетельствующему об актуальности темы, проблема адекватного

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

Цель работы

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

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

Основные направления исследований

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

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

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

В диссертации помучен ряд новых результатов, которые выносятся на защиту:

1. Предложен новый метод решения интегральных уравнений Абеля и Шле-мильха, основанный на разложении искомой функции по полиномам. Для получения коэффициентов разложения искомого решения входные данные разлагаются по соответствующим ортогональным базисным функциям. При этом при вычислении очередных коэффициентов Фурье разлагается только остаток разложения (метод "снятия шкурок с луковицы"). Этот метод имеет наименьший рост вычислительной погрешности с ростом степени полинома по сравнению с классическим способом, основанным на вычислении коэффициентов Фурье всей разлагаемой функции. Метод решения интегрального уравнения Абеля непосредственно применим в геометро-оптических задачах физики плазмы с цилиндрической симметрией как без учета рефракции, так и с ее учетом, и как составная часть полной диагностической задачи при учете дифракции и поглощения диагностического пучка в плазме.

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

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

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

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

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

С использованием предложенного нового метода решения интегрального уравнения Шлемильха измерены эмиссионные характеристики борид-лантанового (ЬаВь) катода в полях до 500 кУ/см и температурах до 1720"С в ускоряющем резонаторе микротрона, что позволило с большей точностью получить оценки предельных токов микротрона, определяемых электродинамическим взаимодействием ускоряемого пучка и ускоряющего резонатора микротрона.

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

астрономы согласились с выводом автора.

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

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

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

С использованием этих данных и разработанного нового метода создан газовый детектор полного поглощения для ультрамягкого рентгена (УМР) на основе многопроволочной пропорциональной камеры, сочетающей 100% эффективность регистрации с рекордным энергетическим разрешением < 6% в указанном диапазоне. Этим детектором впервые измерен спектр УМР из катодного пятна вакуумной дуги.

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

Указанный комплекс программ находит широкое применение в экспериментальных исследованиях как в нашей стране (ИФП им.П.Л.Капицы, ИРЭ РАН), так и за рубежом (Daresbury Laboratory SERC, England, Photon Factory, Tsukuba, Japan).

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

Диссертация содержит результаты 18 оригинальных работ [23] - [40], опубликованных в России и за рубежом в период с 1972 по 1993 годы. Работы докладывались на семинарах в ИФП имени П.Л.Капицы РАН, ФИАН им.П.Н.Лебедева, МГУ им.М.В.Ломоносова, ИАЭ им.И.В.Курчатова, МИФИ, ЛВТА ОИЯИ, ИППИ РАН, Института математики Сибирского отделения РАН в Новосибирске, ИЯФ им.Г.И.Будкера СО РАН, ГЕОХИ им.В.И.Вернадского РАН, у чл.-корр. РАН проф. С.М.Рытова (РИАН им.А.Л.Минца) в России, а также в статистической лаборатории Кембриджского университета, Англия у проф. Д.Кендалла, в Дарсберийской лаборатории в Англии, в ЦЕРНе, Швейцария у проф. Дж.Шарпака, в Photon Factory, КЕК, Tsukuba, Japan, а также на международной школе по вычислительной физике в Чехословакии в 1979 году, на международной конференции по приборам для сян-

хротронного излучения SRI-82 в Гамбурге, Германия, на Всесоюзной конференции по томографии в Новосибирске в 1984 г., на международном совещании по методу максимума энтропии MAXENT Workshop, Cambridge, UK, 1988 и на международной конференции по обратным задачам в Осаке, Япония в 1990 г.

Основные результаты получены в Институте физических проблем имени П.Л.Капицы РАН частично в соавторстве с ЛА.Вайнштейном, E.Pantos, В.Д.Песковым, В.И.Гельфгатом и Е.Р.Подоляком.

Публикации

По теме диссертации опубликована одна обзорная работа [12) и 18 статей [2340] в отечественных и иностранных изданиях. Список научных работ по теме диссертации приведен в конце автореферата.

Личный вклад автора

Работы по новым методам решения интегральных уравнений Абеля и Шлемиль-ха, а также их применениям выполнены автором самостоятельно. В работе по применению разработанного метода решения интегрального уравнения Абеля для задачи по диагностике плазмы в условиях наличия рефракции и дифракции волнового пучка основная роль принадлежит покойному чл.-корр. АН СССР профессору Л.А.Вайнштейну. Работа по оптимальной фильтрации выполнена автором совместно с E.Pantos. Работы по разработке новых рентгеновских детекторов и новых методов обработки данных, получаемых из этих детекторов, выполнены автором совместно с В.Д.Песковым и Е.Р.Подоляком, а по многодырочной камере также и с Г.Ф.Карабаджаком, А.В.Кашлюком и С.В.Переверзевым. Программный комплекс RECOVERY создан автором совместно с В.И.Гельфгатом и Е.Р.Подоляком. Работы о шенноновском пределе сверхразрешения выполнены автором самостоятельно. Автор выражает искреннюю благодарность и признательность всем упомянутым лицам.

Структура и объем диссертации

Диссертация состоит из 9 глав. Она изложена на 232 стр., содержит 61 рисунок, 11 таблиц, тексты 12 программ для ЭВМ на 35 стр. и список литературы из 273 наименований.

Содержание диссертации

Оригинальные результаты, представленные в диссертации, опубликованы в работах [23-40] и могут быть быть кратко сформулированы следующим образом:

Глава 1. Введение

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

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

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

Глава 2. Задача о токе эмиссии и ортогональные полиномы

Задача сводится к решению интегрального уравнения Шлемильха

«/2 -г/2

./(£) = — / 1(Есоъ<р)4<р= - I /(Евт^сЬр. (1)

-1г/2 О

в котором 1(Е) искомая функция, а J(E) экспериментальные данные. Точное решение уравнения (1) в виде

г/г

/(£) = 2.7(0)+ I У{Е*т>р)<1<р. (2)

о

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

Основой для предлагаемого нами нового способа численного решения интегрального уравнения (1) является следующий факт. Функции Ет при т = 1,2,... являются собственными функциями линейного интегрального оператора в правой части уравнения (1). В самом деле, пусть

1(Е) = Ет,

тогда

г/г

У(Е) = - [ ЕГыпт<р<1ч> = ХпЕ",

7Г J

О

где собственные числа Дт, равные

. 1 / . т . .11 т

1 - / ет рам = -—^-V,

2 г(т + 0

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

Ло=1/2. Л| = 1/тг, >т = Лт_2(1 - 1/т), то = 2,3..

Нетрудно проверить, что функции Ет являются собственными функциями и для преобразования от 1 к / согласно уравнению (2).

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

м к-0

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

Сначала мы действуем стандартно: для нахождения разложения экспериментальной функции у(х,) по обычным полиномам находим разложение этой функции по полиномам Чебышева {^(х,)}, к = 0,1,2,...,IV — 1, ортогональным на заданной дискретной совокупности входных точек {г;}, г = 1,2, ...,ЛГ, после чего линейным преобразованием коэффициентов разложения по полиномам Чебышева находим разложение по обычным полиномам. Однако при получении каждого очередного М 4-1-го коэффициента разложения по ортогональным полиномам мы не пользуемся стандартной формулой

ом+1 = (у^м-н) (3)

(здесь и далее в формулах (5) и (6) круглые скобки означают скалярное произведение), а вместо этого сначала находим остаток аппроксимации

м

Ям(х) = у(х)-^ск<рк(х) (4)

к-О

и только затем следующий коэффициент разложения

см-и = (Ям,рм+0- (5)

(см. [25]). Аппроксимация продолжается вплоть до такой степени М, при которой впервые выполняется условие

(Ям, Ям) < (¿2/, ¿2/)- (6)

По предложению нобелевского лауреата профессора А.Клуга из лаборатории молекулярной биологии Кембриджского университета Англии совокупность формул (4) и (5) может быть названа процессом "снятия шкурок, с луковицы".

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

Задача о токе эмиссии в микротроне была решена автором в 1972 г. [23], что позволило с большим пониманием подойти к оценке предельных токов в микротроне, связанных с взаимодействием пучка с ускоряющим резонатором. Фактически это была первая обратная задача экспериментальной физики, рассмотренная автором. На рис.1 показаны входные данные и выходные результаты задачи о токе эмиссии. Из рисунка видно, что величина погрешности решения примерно равна погрешности входных данных, что указывает на хорошую обусловленность разработанного нового способа решения интегрального уравнения Шлемильха.

л А- ,3

> -п*

Г у

у

У

й /00 200 ХР А00 ¿00

Напряженность поля {кУ/см)

Рис. 1 Входные данные (слева) и выходные результаты (справа) задачи о токе эмиссии. Цифры на графиках соответствуют различной температуре борид-лантанового катода: 1 = 1600°С, 2 = 1650°С, 3 = 1680°С, 4 = 1720"С.

Глава 3. Интегральное уравнение Абеля и звездная задача

Интегральное уравнение Абеля

Р{х)= 1 ^у====У О < V < х, (7)

в котором Р(х) заданная, а /(у) искомая функция, впервые было введено и решено Н.Абелем в 1826 г. в связи с обобщением задачи механики о таутохроне, и фактически это была первая задача физики, непосредственно приводящая к интегральному уравнению первого рода. Это уравнение имеет столь многообразные применения в различных областях науки, что их полное перечисление заняло бы здесь слишком много места. Отметим только два, рассматриваемых в диссертации: в физике плазмы при диагностике цилиндрических разрядов волновыми пучками [26] и астрономии при определении пространственной плотности звезд в шаровых звездных скоплениях [27], [28].

Оба интегральных уравнения: Шлемильха (1) и Абеля (7) тесно связаны и заменой переменных могут быть преобразованы одно в другое. Точное решение

интегрального уравнения Абеля

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

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

измерений, так и аналитические свойства преобразования Абеля, выражаемые формулой (7). Такой способ был предложен автором в 1973 году [24]. Он очень близок к описанному выше способу решения интегрального уравнения Шлемильха и основан на разложении искомой функции /(ж) по собственным функциям интегрального оператора

)/(* - у)

которыми являются степенные функции }„(х) = с" для п =0, 1, 2,... Соответствующие собственные значения вычисляются по рекуррентной формуле

А0= 1, АП = А„_,/(Ц-1/2п), п= 1,2,....

Таким образом, если /(¡с) = хп, то Р^) = 2Х„хпу/х. Отличие F(x) от последнего выражения может быть ооусловлено только ошибками измерений.

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

Я»)-Ел*\ 0 < х < а. (8)

1.-0

В этом случае экспериментальная функция Р(х) необходимо представляется в виде

т

^(аО-М'ОХЛЛ«*, (9)

А-0

где <р(х) = 2у/х.

Разложение (9) находится методом среднеквадратичной аппроксимации функциями

Фц(®) = ср{х)Рк(х), ортогональными на заданном отрезке 0 < х < а с весом

Цг) ~ 1/[№(х)]2,

где 6Р(х) > 0 - среднеквадратичная погрешность измерения функции ^ в точке х, а Р*(х) - ортогональные полиномы на отрезке [0,а] с весовой функцией д(х) = ■ш(х)^2(х).

В зависимости от того, как задана функция F(x), строятся разные системы полиномов. При непрерывном задании F(x) и v>(x) = const этими полиномами являются классические полиномы Якоби Р^0,|'((2е - а)/а). При задании .F(z) в дискретной системе точек zi,z2,...,xn строится указанная выше система полиномов Чебышева, ортонормированных на указанном множестве точек. Степень то аппроксимирующего полинома выбирается методом регрессионного анализа с использованием статистики

к-О

¡5F(x.

4'

которая распределена как Хп-т-1-

Одновременно с вычислением ортогональных полиномов определяются их коэффициенты и тем самым искомое разложение

F{x) = ¥>(z) £ Fkxk

к-О

(10)

После этого мы вычисляем

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

Sf(x) =2

£ {DFUzax»/{\aА,)

«,/9-0

1/2

(И)

В этом выражении ИР - ковариационная матрица (матрица ошибок) коэффициентов определяемая при получении разложения (10).

Этот коридор ошибок, или неустранимая погрешность, приближенно соответствует 95%-ной доверительной зоне для решения /(г). В практических задачах основным источником погрешности восстановления является наличие случайных ошибок в исходных данных. Поэтому оценка погрешности (11) неизбежно имеет вероятностный характер.

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

В задаче о шаровых звездных скоплениях в предположении сферической симметрии пространственная плотность звезд 1(т) связана с наблюдаемой поверхностной плотностью J(R) (которая в этом предположении будет иметь круговую симметрию) следующим интегральным уравнением первого рода:

Ro

J(R)-2f:

I(r)rdr

, 0 < R < Ro,

в котором Яо есть радиус скопления такой, что при R > Ro

I(R)=J{R) = 0.

Это уравнение заменой переменных сводится к интегральному уравнению Абеля (7).

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

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

статистикой.

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

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

Рис.2 Пространственная плотность вспыхивающих звезд в Плеядах.

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

3

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

Глава 4. Оптимальная фильтрация данных с шумами

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

Пусть мы имеем измеренные значения /,• функции }(х) в дискретных точках ¡ti, 12,..., ijv с соответствующими ошибками измерений (шумом) e¡

Л =/(*;)i - 1, 2,..., N и рассмотрим случай, когда нет систематических ошибок, т.е.

e¡ = 0,

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

Т^] = a2S<j.

Здесь Sij - есть символ Кронекера и <г2 - есть дисперсия шума, одинаковая во всех точках.

Функция /(г) может быть разложена в ряд по некоторым ортогональным базисным функциям {уД®)}

« = 1,2,... (12)

а

со следующим соотношением ортогональности

1

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

Хорошо известно, что минимальная ошибка аппроксимации

пипЯд^, а = 1, 2,..,, М,...

(13)

где

rm = ÜL

, М =1,2,...

получается если для коэффициентов разложения са используются коэффициенты Фурье

Сс, - (/,¥>«) = £

Поскольку из наблюдений нам известны только значения /¡, включающие случайный шум то минимизация (13) не может быть выполнена в явном виде. Вместо этого мы можем искать минимум среднего значения тт В?м, где

Щл ~

с* = (/,*>«) = 53 /.^(АО-

Минимум величины Я2М соответствует минимуму разности (в среднеквадратичном смысле) между неизвестной функцией f(x) и ее аппроксимацией.

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

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

/|(г) = 53*«г<^°(г:)

а

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

й2 = 53 Я^) - ^.ёаУеЛ®;

«I. о

а

£>(<:„) = (Яс)«,.

Здесь обозначено

Действуя стандартным методом,

мы находим эти коэффициенты

с1 + р(сау

(14)

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

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

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

Минимальное значение К2 равно

— с2 д 2 . = у >_

и для погрешности аппроксимации функции fi(x) мы имеем выражение

*/,(*)- 2а ' .

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

}{х), которую мы хотим узнать. Множитель 2 в этой формуле приблизительно соответствует обычно принимаемому 95% доверительному интервалу.

Важно отметить, что непосредственное использование формулы фильтра Винера (14) для целей оптимальной фильтрации невозможно, поскольку в эту формулу входят неизвестные нам амплитуды с». Если бы мы знали эти амплитуды, мы могли бы тогда сконструировать оптимальный фильтр для их выделения из шума, но в этом случае их выделять вовсе не нужно, поскольку мы считаем их известными!

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

Мы предлагаем следующий алгоритм для использования оптимальной фильтрации. Вместо неявной формулы (14) для фильтра Винера мы используем приближенную явную формулу

i.__

а % + mv

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

¿I > D(ca) для a < ао, и простую прямолинейную экстраполяцию

ln(c^) & Aia 4- Si дня a<¡ < a < a¡ для той части спектра, где сигнал слабее, чем шум

с I < D(ca)

и полностью тонет в нем.

Мы убедились, чго такой простой подход вполне достаточен для большинства спектроскопических данных. Однако имеются случаи, например, для данных ЕХ-AFS экспериментов (представляющих собой затухающие синусоиды), где скорость уменьшения амплитуд спектра Фурье лучше аппроксимируется выражением

1п(г2)^Л21па+В2.

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

Приведенный здесь вывод винеровского фильтра для неслучайного сигнала и случайного шума вместе с полным текстом соответствующей программы оптимальной фильтрации для ЭВМ были впервые опубликованы автором и Е.Пантосом в препринте Дарсберийской лаборатории в 1982 г. Впоследствии примерно то же самое, ю без указания конкретного алгоритма экстраполяции и без ссылки на работы [29] i [30], было опубликовано (по-видимому, совершенно независимо) в полезной книге 'Numerical Recipes" (The Art of Scientific Computing) by W.H.Press et all, Cambridge Jniversity Press, Cambridge, U.K., 1 изд. 1989, 2 изд. 1992. Программа оптимальной фильтрации с использованием быстрого преобразования Фурье на основе описанного 1Лгоритма широко используется как в нашей стране, так и за рубежом. Текст ее гриведен в гл.9 диссертации.

Глава 5. Многодырочная камера для мягкого рентгена

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

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

Описанные в этой главе результаты получены автором совместно с Г.Ф.Карабад-гаком, А.В.Кашлюком, С.В.Переверзевым и В.Д.Песковым [34].

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

Работа, описываемая в этой главе, выполненная совместно с В.Д.Песковым и i.Р.Подол яком [31], также бьша напрямую связана с исследованиям по физике лазмы газового разряда в мощных СВЧ полях, проводившихся П.Л.Капицей. Уже

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

Это разрешение представлялось недостаточным для решения поставленных за-

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

Для решения первой из указанных задач были собраны и критически проанализированы все имевшиеся экспериментальные данные по фотопоглощению рентгеновского излучения в молекулярном водороде, а Е.Р.Подоляком и В.И.Алавердовым были выполнены дополнительные измерения. Результаты анализа опубликованы е работах (32) и [33], где приведена аппроксимация полученных данных, имеющая погрешность около 5% с доверительной вероятностью 95% в диапазоне 17-400 эВ.

Задача восстановления спектра рентгеновского излучения по измерению поглощения сводится к решению следующего интегрального уравнения первого рода

Щх) = J ДЕ)[1 - е-а(в)4]е-*'в'*(гя, (15)

в котором х - обозначает глубину слоя поглощения рентгеновского излучения в газе N(x) - интенсивность излучения на глубине х, ¡(Е) - искомый спектр рентгенов ского излучения на входе слоя поглощения газа толщиной х, а{Е) - коэффициент поглощения рентгеновского излучения как функция энергии Е, Д - размер одно{ секции пропорционального счетчика для регистрации излучения, так что множител!

w(E) = 1 — ехр[—а(В)Д]

представляет собой эффективность регистрации, пропорциональную вероятности по глощения рентгеновского фотона с энергией Е в одной секции счетчика толщиной Д. Интегрирование в (15) производится по диапазону энергий фотонов ДВ, гд( w(E) > 0.

Интегральное уравнение (15) есть частный случай общего уравнения (1) и: Введения к автореферату, где входные данные являются случайными, поскольк; в большинстве случаев N(x) имеет пуассоновское распределение при каждом х После дискретизации интегрального уравнения (15) получаем систему линейны: алгебраических уравнений для определения неизвестных ]k = J(Ek), fc=l,2,...,m и. входных данных Ni=N(xi), ¡=1,2,...,п

Ni = jrhPik, ¿=l,2,...,n, (16

в которой Pik = SE -[1 - е-°(£»)д] • е""^»)1', а SE = Ек+\ - Ек, k= 1,2,...,m - 1 ша дискретизации по энергии фотонов.

Подчеркнем, что число уравнений п в системе (16) может быть как больше, Tai и меньше количества неизвестных т. Для случайных входных данных JV; = Nix, в подавляющем большинстве случаев система (16) не имеет решения. Именш случайность значений и приводит к несовместности системы (16) в отличие а уравнения (I) при точной правой части.

Для придания смысла решению системы (16) примем во внимание, что случайны! величины {А^} в совокупности независимы и могут рассматриваться имеющим! полиномиальное распределение, при котором вероятность V получить из опыт;

точно указанные выше значения равна

ЛЛ "

п щ .-I i-l

Здесь введены следующие обозначения

.. А .r ^ Pik _ ехр[-а(£*)х,-]

i-l к-1 2 2 ехр[-а(Д*)^1

i-l i-1

i-l

Решением системы (16) согласно методу максимума правдоподобия (МП), назовем вектор {дь}, при котором вероятность V или logP достигают максимальных значений. Величину L = logP будем называть далее функцией правдоподобия.

Максимум функции правдоподобия всегда существует, так как она непрерывна, а выражение

п п

L = const 4- N¡log Р< = const N 53 /»'°S Р» > .-i i-i

-де fi = N,/N, ограничено сверху, поскольку

53 fM pí - 23 /.- (?.//.) < о ¿ i »

s силу неравенства теории информации.

Поиск максимума функции правдоподобия производится с помощью итерационной гроцедуры, представляющей собой модификацию итерационной процедуры, ранее 1редлохенной М.З.Тараско. В этом методе не делается никаких предположений о шде восстанавливаемой функции 1(E), кроме того, что 1(E) > 0. Не требуется также шформации о начальном приближении, в качестве которого мы обычно выбираем

Г(Е) = const > 0.

В явном виде итерационная формула может быть записана в следующей форме

/ \

Зк -Як

л -i

(17)

!десь з =1,2,... - номер итерации, к - длина шага в пространстве искомых векторов д}. При к = 1 итерационная формула (17) совпадает с итерационной формулой -1.3.Тараско.

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

Итерации продолжаются пока величина

* ~ к Ъ '

где ¿А^- = N(pi — /;), больше, чем уровень гДе Я - уровень значимости х2

критерия, обычно принимаемый равному 5%, и должны останавливаться, как только эта величина уменьшится до уровня

где Рг = (80 — 95)%.

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

Оценка погрешности восстановленного спектра 61(Е) производится методом статистического моделирования повторением процедуры восстановления 10-20 раз для различных случайных наборов исходных данных с одной и той же функцией распределения. На этом закончено краткое описание процедуры восстановления методом МП для входных данных с биномиальным или с пуассоновским распределениями. Модификация итерационной формулы (17) для входных данных с гауссовским распределением описана в гл.7 диссертации.

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

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

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

На рисунке 3 ниже приведены рентгеновские спектры, полученные в опьгге с рентгеновской трубкой при напряжениях на трубке 410, 448 и 470 вольт. Спектры восстановлены описанным методом из данных по поглощению. На рисунке видно, что верхняя граница спектра равна напряжению на рентгеновской трубке в каждом из трех независимых опытов. Кроме того, в каждом из опытов видна одна и та же линия углерода СУ„, Е = 275 эВ. Полуширина линии, равная разрешению счетчика, примерно равна 15 эВ (<6%). Лучшего значения трудно было ожидать, так как данные по поглощению, определяющие ядро интегрального уравнения (15), известны только с точностью около 5%. Указанное энергетическое разрешение для газовых

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

W

ВЫ)

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

Глава 7. Пакет программ восстановления сигналов RECOVERY и его применения

Материал данной главы основан на совместной работы автора диссертации с В.И.Гельфгатом и Е.Р.Подоляком [38], депонированной в ВИНИТИ в 1991 г., а также опубликованной на английском языке в 1993 г. в журнале Computer Physics Communications [39].

Восстановление сигналов сводится к решению линейного интегрального уравнения первого рода

? К{*,у)ШЛу = Рй(х), c<x<d, (18)

Ja

i котором fo(y) искомая функция (сигнал), Р0(е) результат преобразования (иска-кения) сигнала измерительным прибором, характеризуемым аппаратной функцией К{х,у) - ядром интегрального уравнения, причем результатом измерения является ie функция Fo(i), а случайная функция F(x), лишь среднее значение которой по мзным реализациям шума равно правой части Fo(x) этого интегрального уравнения

F(z) = Fo(z). (19)

Это условие фактически означает отсутствие систематических ошибок при измере-щях. Кроме того, принимается, что шум некоррелирован в соседних точках

6F(xi) ■ SFix,-} = <r? • 6ц, где SF(x) = F (s) - Fo(z),

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

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

измерительным прибором в присутствии шума. Шум в отдельных точках независим и может иметь гауссовское, биномиальное или пуассоновское распределения. Программы комплекса реализуют решение интегрального уравнения свертки с ядром К{х,у) = К(х — у) дая коррекции на аппаратную функцию произвольного вида, зависящую только ог разности аргументов (эта задача весьма распространена в спектроскопии любого рода), решение интегрального уравнения Лапласа в действительной области

оо

f(t) = I j(A)exp(—А()(Ц, о

в котором /(t) заданная функция, a 5(A) искомая функция для разложения заданной экспериментальной функции на экспоненты с непрерывным спектром декрементов (эта задача имеет многочисленные применения, например, в ядерной физике или оптике при исследовании временной структуры процессов флуоресценции), и решение интегрального уравнения (15) (см.выше) для восстановления спектра ультрамягкого рентгеновского излучения из измерений его поглощения в газе. Программы основаны на методе максимума правдоподобия, при использовании которого достигается максимальное улучшение разрешения сигналов.

Комплекс программ состоит из набора трех главных подпрограмм MLU8, MLGS и MLP8, реализующих итерационный способ нахождения максимума функции правдоподобия для неотрицательных решений соответственно для случаев гауссовского распределения ошибок во входных данных с постоянной величиной дисперсии ctq во всех экспериментальных точках, гауссовского распределения для случая, когда в разных точках дисперсия может быть разной <т?, г =1,2, ...,п, и для пуассоновско-го распределения входных данных. Подпрограмма MLPS требует положительности аппаратной функции, для подпрограмм MLU& и MLGS аппаратная функция может быть знакопеременной.

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

В работах [38] и [39] приведены полные тексты на Фортране всех программ и подпрограмм комплекса RECOVERY, реализующие метод скорейшего подъема при поиске максимума функции правдоподобия. Это дает возможность применять эти программы на любых ЭВМ, имеющих фортрановский компилятор. Огедует упомянуть единственное ограничение, состоящее в том, что графическая часть пакета рассчитана на использование IBM-совместимых персональных ЭВМ и компилятора IBM Fortran/2.

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

э

В

о

г

1

Н5Ь 38^ 512 Ь<*0 768

Рис.4 Восстановление сигнала, программой йсопу 2 для аппаратной функции с лорснтцев-схим профилем шириной В » 60, отношением сигнал/шум 20 дБ при числе данных N = 329. 1 - входные данные, 2 - аппаратная функция, 3 - результат восстановления, 4 ■ точный ответ. В целях графической ясности кривая 1 сдвинута по вертикали вверх. Достигнутое значение Х2-крптерпя равно х= 0.903

Глава 8. Шенноновский предел сверхразрешения при восстановлении сигналов

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

Впервые результаты этой главы бьии доложены на международном совещании по методу максимума энтропии в Англии в 1988 г. и опубликованы в работах [35-37] в 1989-1990 гг. Ввиду сравнительной новизны и важности этого результата приведем здесь его простой эвристический вывод, следуя работе [36]. Подробное формальное доказательство опубликовано в [37] и воспроизведено в главе 8 диссертации.

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

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

Для случая, когда аппаратная функция К(х - у) имеет ограниченную ширину тространственного спектра Фурье Ж, существование предела разрешения близких :игналов /о(у) следует из известной теореме К.Шеннона 1949 г. о предельной

К(х,у) = К(х-у)

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

(В1Х)1аят/си) = Wc„- log2( 1 + Е./Еп). (20)

Здесь обозначено через X = d—с, В/X есть максимальное число бит информации на единицу измерения переменной х, различимых при данных уровнях энергии сигнала

Е. = Г Ft(x)dx

J-со

и шума

£„ = п • <х2,

п. равно числу точек измерения экспериментальной функции

F(x)~Fd(x) + N(X),

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

Если мы определим ширину аппаратной функции как эквивалентную ширину ее квадрата из соотношения

Д= Г K2(x)dx , (21)

J-оо

при условии, что if(0) = 1, и обозначим через S минимальное расстояние, на котором могут быть различимы близкие сигналы /о(у), то коэффициентом сверхразрешения (Super Resolution) мы будем называть отношение

SR^A/S. (22)

Доказательство теоремы Шеннона, данное самим Шенноном, основано на разложении правой части F(x) в ряд Когельникова по базисным функциям с ограниченным спектром

•i=i-2'3.....

7г(2Wx - г)

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

X = d - с

имеется лишь конечное число отсчетов

п = 2WX.

Для аппаратных функций с неограниченным спектром (как, например, для ло-рентцевского или гауссовского профилей) можно вместо разложения по ряду Котельникова использовать разложение по любому другому подходящему базису. Этот метод решения интегральных уравнений первого рода описан, в работе [12], где он назван методом ортогональных разложений (МОР), причем, как показано в работах [29] и [30] для получения минимальной среднеквадратичной погрешности решения необходимо использовать во всех рядах сравнительно небольшое число m. < п членов разложения.

А.Н.Колмогоровым указано, что для справедливости теоремы Шеннона фактически достаточно конечномерности пространства сигналов, что имеет место практически

для всех типов экспериментальных данных F(x), поскольку функция F(x) всегда измеряется в конечном числе экспериментальных точек.

Естественно что формула Шеннона должна быть изменена для базисных функций с неограниченным спектром. Из соображений размерности вместо (20) мы можем

записать , ,, „ .

Д/6 = Const х log2(l + Е./Еп).

Численное значение константы, которое вообще говоря, может зависеть как от формы аппаратной функции, так и от формы различаемых сигналов не может быть определено только из анализа размерности. Для этой цели в работах [35-37J используются два способа:

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

- численные эксперименты.

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

= |log2(l

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

6

0 130 260 390 520

'ис.5 Восстановление двух ланий на расстояниях больше и меньше порога разрешения

На рис.5 показан пример восстановления двух линий, свернутых с лорентцевской аппаратной функцией с параметром D = 60 и наличии шума 20 дБ при расстоянии между линиями: a) S2 — 35 (больше порога разрешения) и б) 52 = 28 (меньше порога разрешения).

Глава 9. Тексты программ

В этой главе приводятся тексты программ, описанных в диссертации:

- Программы аппроксимации экспериментальных данных полиномом методом ортогональных разложений с новым способом нахождения коэффициентов Фурье ("снятие шкурок с луковицы") CODAI, CODA2 на Алголе для ЭВМ НР-1000 и CODAF на Фортране;

- Программы доя решения интегрального уравнения Шлемильха САТЗ и САТ6 на Алголе для ЭВМ НР-1000;

- Программа для решения задачи диагностики цилиндрического плазменного разряда с учетом рефракции на основе решения интегрального уравнения Абеля PHASE на Алголе для ЭВМ НР-1000;

- Программа оптимальной фильтрации FILT и основная подпрограмма FILTER, реализующая алгоритм описанный в гл.4, на Фортране;

- Главные программы пакета RECOVERY: Dconv, Dexpa и Datten для решения интегрального уравнения свертки, интегрального уравнения Лапласа в действительной области и для восстановления спектра рентгеновского излучения по измерениям его поглощения в газе (IBM Fortran/2).

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

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

2. Разработанный математический аппарат и алгоритмы решения обратных задач, изложенные в диссертации, дали возможность создать комплексы программ обработки и восстановления сигналов, применяемые как в нашей стране (ИФП им.П.Л.Капицы, ИРЭ РАН, НИИ Автоматических Систем, МКБ "Факел"), так и за рубежом (Daresbury Laboratory SERC, England, Photon Factory, Tsukuba, Japan).

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

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

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

Литература

[1] Амбарцумян В.А. Обратные задачи астрофизики. В сборнике "Фундаментальные проблемы теоретической и математической физики", стр,9-20, D 12831, Дубна, 1979

[2] Горелик Г.С. О применении модуляционного метода в оптической интерферометрии. ДАН СССР, том 83, стр.549-552, 1952

[3] Бернштейн И.Л., Горелик Г.С. К теории звездного интерферометра Майкельсона. ДАН СССР, т.86, стр.47-50, 1952

[4] Курант Р. и Гильберт Д. Методы математической физики. Пер. с нем., том 1, стр.147, ГТТИ, М., 1932

[5] Тихонов А.Н. Об устойчивости обратных задач. ДАН СССР, том 39, N 5, стр.195-198, 1943

[6] Раутиан С.Г. Реальные спектральные приборы. УФН, т.66, вып.З, стр.475-517, 1958

[7] Турчин В.Ф., Козлов В.П., Малкевич М.С. Использование методов математической статистики для решения некорректных задач. УФН, т. 102, вып.З, стр.345-386, 1970

[8] Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. Наука, ФМ, М., 1 изд., 1974, 3 изд., 1986

[9] Иванов В.К., Васин В.В., Танана В.П. Теория линейных некорректных задач и ее приложения. Наука, М., 1978

10] Frieden B.R. Image Enhancement and Restoration. In "Picture Processing and Digital Filtering", pp.179-248, ed.T.S.Huang, Springer-Verlag, Derlin - N.Y., 1979

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

12] Kosarev E.L. Application of the 1-st Kind Integral Equations in Experimental Physics. Comput. Phys, Comm. Vol.20., No.l, pp. 69-75, 1980

13] Blass W.E., Halsey B.W. Deconvolution of Absorption Spectra. - N.Y. : Acad. Press, 1981. - 160 P.

14] Mendel J.M. Optimal Seismic Deconvolution. Academic Press, N.Y., 1983

[15] Deconvolution with application in spectroscopy. Ed. P.A.Jansson. Academic Press, N.Y., 1984

[16] Алгоритмы и программы восстановления зависимостей / В.Н. Вапник, Т.Г. Глазкова, В.А. Кощеев и др.; Под ред. В.Н. Вапника. - М. : Наука, 1984. -815 с.

[17] Верлань А.'Ф., Сизиков B.C. Интегральные уравнения. Методы. Алгоритмы. Программы : Справочное пособие. - Киев : Наукова думка, 1986. - 543 с.

[18] Василенко Г.И., Тараторкин A.M. Восстановление изображений. - М. : Радио и связь, 1986. - 302 с.

[19] Реконструкция изображений. Под ред. Г.Старка. Пер. с англ. изд. Academic Press, 1987, Мир, М., 1992

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

[21] Пытьев Ю.П. Методы анализа и интерпретации эксперимента. Изд. МГУ, М., 1990

[22] Федотов A.M. Некорректные задачи со случайными ошибками в данных. Наука, Сиб.отд., Новосибирск, 1990

Основные работы автора по теме диссертации

[23] Е.Л.Косарев. Измерение катодных потерь в микротроне. ЖТФ, т.42, N4, 841 -843, 1972.

[24] Косарев Е.Л. О численном решении интегрального уравнения Абеля. ЖВМ и МФ, том 13, N 6, стр. 1591-1596, 1973

[25] Е.Л.Косарев. Аппроксимации экспериментальных данных полиномом методом ортогональных разложений. Алгоритмы N П 000257 в сб.: Алгоритм и программы, N1, 1973, изд. ВНТИЦЕНТР, М., 1973.

[26] Л.А.Вайнштейн, Е.С.Биргер, Н.Б.Конюхова, Е.Л.Косарев, Г.П.Прудковский. Коротковолновая диагностика плазменного шнура. Физика плазмы, т. 2, вып. 4, стр. 658 - 671, 1976.

[27] Е.Л.Косарев. Новый метод восстановления пространственной плотности звезд в шаровых скоплениях и его применение к вспыхивающим звездам в Плеядах. Письма <з астрономический журнал, т.6, N7, стр. 408-413. 1980

[28] Е.Л.Косарев. Пространственное распределение вспыхивающих звезд в Плеядах, вычисленное методом ортогональных разложений по данным наблюдений. Переменные звезды, т. 21, N3, стр. 321-363, 1980

[29] Kosarev E.L., Pantos Е. Optimal Smoothing of Noisy Data by Fast Fourier Transform. J. Phys. E. : Sei. Instrum. Vol.16, pp. 537-543, 1983

[30] Косарев Е.Л., Пантос E. Оптимальное сглаживание данных с шумом, использующее быстрое преобразование Фурье. Приборы и техника эксперимента -1985. - Вып. 3. - С. 92-95.

¡31] Косарев Е. Л., Песков В. Д., Подоляк Е. Р. Восстановление спектра ультрамягкого рентгеновского излучения из измерений его поглощения в газе. ЖТФ. -19S3. - Т. 53, вып. 6. - С. 1101-1114.

E.L.Kosarev, V.D.Peskov and E.R.Podolyak. High resolution soft X-ray spectrum reconstruction by MVVPC attenuation measurements. NIM, Vol.208, pp.637-645, 1983

[32] Косарев Е.Л., Подоляк E.P. Фотопоглощение рентгеновского излучения в молекулярном водороде в диапазоне энергий фотонов 17-400 эВ. Оптика и спектроскопия, т.56, вып.4, стр.643-646, 1984

[33] E.L.Kosarev, E.R.Podolyak. Comments on the data on photoabsorption in molecular hydrogen for low-energy X-rays. Nuclear Instruments and Methods in Physics Research, A261, 161-162, 1987

[34] Г.Ф.Карабаджак, А.В.Кашлюк, Е.Л.Косарев, С.В.Переверзев, В.Д.Песков. Многодырочная камера-обскура для наблюдения слабого импульсного рентеновского излучения. ПТЭ, вып.2, стр.181-186, 1988

[35] Косарев Е. Л. Шенноновский предел сверхразрешения и его достижение при восстановлении сигналов. ПТЭ. - 1989. - Вып.4. - С. 84-87.

[36] E.L.Kosarev. Superresolution limit for signal recovery. In J.Skilling (ed.), Maximum Entropy and Bayesian Methods, Kluwer Academic Publ., Dordrecht, 1989, pp.475-480

[37] Косарев E. Л. О пределе сверхразрешения при восстановлении сигналов. Радиотехника и электроника, т.35, вып.1, стр. 68-87, 1990

E.L.Kosarev. Shannon's superresolution limit for signal recovery. Inverse Problems, Vol.6, No.l, pp. 55-76, 1990

38] Гельфгат В.И., Косарев Е.Л., Подоляк E.P. Комплекс программ восстановления сигналов из зашумленных данных методом максимума правдоподобия. Приборы и техника эксперимента, вып.5, стр.93-98, 1991

39] Гельфгат В.И., Косарев Е.Л., Подоляк Е.Р. Комплекс программ восстановления сигналов из зашумленных данных методом максимума правдоподобия. Депонированная статья в ВИНИТИ, номер 2635-В91, 112 стр., 1991

40] V.I.Gelfgat, E.L.Kosarev, E.R.Podolyak. Programs for signal recovery from noisy data using the maximum likelihood principle.

Part 1: General description. CPC, Vol.74, No.3, pp.335-348, 1993 Part 2: Program implementation, ibid, pp.349-357, 1993

Рукопись поступила в издательский отдел 21 октября 1994 года.