О РЕГУЛЯРИЗАЦИИ ЗАДАЧИ ТОМОГРАФИИ ИОНОСФЕРЫ.

Д.В.Кирьянов

МГУ им. М.В.Ломоносова



Обсуждается метод томографической реконструкции поля ионосферной электронной плотности, основанный на классической тихоновской регуляризации. Приводится соответствующий вычислительный алгоритм и результаты численного эксперимента.
Введение. Методы спутниковой радиотомографии получили в последнее время значительное распространение для ионосферных исследований. Основная идея лучевой томографии состоит в численном решении системы большого числа интегральных уравнений, дискретизация который приводит к в виде системе линейных уравнений с достаточно разреженной матрицей А [Куницын, Терещенко]:
AN=f (1)
где N(J) – искомый вектор значений электронной концентрации в пространственных узлах сетки, содержащий MJ элементов; K – номер луча (всего лучей MK), правая часть уравнений fK представляет собой результаты измерений (чаще всего, доплеровского смещения частоты). По повторяющимся индексам здесь и далее предполагается суммирование.
Система (1) является переопределенной, т.е. количество уравнений превышает (как правило, в несколько раз) число неизвестных. Поэтому искомый вектор можно интерпретировать только как «псевдорешение» [Воеводин] или «обобщенное решение» [Рябенький] системы (1), т.е. вектор из пространства RM, минимизирующий функционал невязки
. (2)
Здесь двойными вертикальными чертами обозначена евклидова норма вектора невязки системы (1).
Для решения системы (1), ввиду большой размерности задачи, используются, в основном, приближенные итерационные методы. Одним из наиболее часто применяемых алгоритмов являются итерационные алгоритмы семейства ART [Ценсор], которые своим развитием обязаны технологии интерпретации изображений. В частности, в ионосферной томографии широко применяется алгоритм ART3. Поскольку матрица А является плохо обусловленной, то точное решение (1) оказывается крайне неустойчивым и, как правило, весьма далеким от реального поля электронной концентрации в ионосфере. Применение алгебраических алгоритмов реконструкции, в которых учитывается, помимо прочего, положительная определенность вектора , как показало многочисленное моделирование [Kunitsyn et.al.’1995], позволяет в значительной степени обойти указанные трудности. Однако задача (1) остается некорректной, поэтому вопросы о единственности и устойчивости наеденного решения остаются открытыми, что дает недоброжелателям ионосферной томографии повод усомниться в достоверности результатов реконструкции [Данилкин, Тольский].
Настоящая работа связана с предложением нового алгоритма реконструкции, основанного на методе классической тихоновской регуляризации. Такой подход позволяет корректно сформулировать задачу ионосферной томографии и гарантирует существование, единственность и устойчивость получаемой реконструкции.
Теория. Как известно, задачу поиска обобщенного решения исходной вещественной системы (1) можно свести к задаче поиска классического решения системы уравнений [Рябенький]:
. (3)
(Символ «Т» обозначает операцию транспонирования.)
Однако задача (3), так же, как и (1), остается некорректной вследствие плохой обусловленности их матриц. Плохая обусловленность системы означает, что погрешности коэффициентов ее матрицы (имеющие место ввиду ошибок аппроксимации) и погрешности правых частей (шум измерений) сильно искажают решение. Поскольку задача неустойчива, то (точное) решение задачи с неточно заданными матрицей и правой частью может не иметь ничего общего с искомым вектором электронной концентрации.
Общеупотребительный подход к решению некорректных задач заключается в их регуляризации [Тихонов, Арсенин]. При этом используется дополнительная априорная информация о решении, которая может быть как качественной, так и количественной. Например можно искать решение, максимально близкое к некоторому профилю N0(z), т.е. к некоторому вектору N0 Концепция регуляризации применительно к (1) сводится к замене задачи (2) на задачу о минимизации функционала
, (4)
где a – малый положительный параметр регуляризации.
Как показано в [Калиткин], задача о минимизации функционала (4) эквивалентна системе, содержащей, как и система (3), MK линейных уравнений:
, (5)
где символ E обозначает единичную матрицу размера MKхMK. Решая систему (5), можно получить регуляризованное решение, зависящее от a. Из (5) хорошо ясен смысл параметра a: при малых a~0 обусловленность системы (5) близка к плохой обусловленности (1), а при больших a система (5) обусловлена хорошо, но ее решение далеко от решения искомой задачи томографии. А именно, чем больше a, тем ближе решение к N0. Очевидно, что на практике необходимо выбирать промежуточные a. Принцип выбора конкретного значения параметра a будет обсужден ниже при анализе результатов численного эксперимента.
Описанный регуляризационный подход сводит некорректную задачу томографии к условно-корректной (по Тихонову) [Тихонов, Арсенин] задаче отыскания решения системы (5), которое, в силу линейности задачи, является единственным и устойчивым.
Компьютерное моделирование. Проиллюстрируем сказанное конкретным вычислительным примером. Зададим некоторую модель ионосферы ( рис.1), которая образована наложением на регулярный высотный профиль неоднородности, расположенной в центре, и промоделируем процесс доплеровских измерений. Высота ИСЗ в компьютерном эксперименте составляла 900 км, два приемника располагались в точках с координатами 0 км и 500 км по поверхности Земли в плоскости пролета ИСЗ. Дискретизация интегральных уравнений проводилась [Куницын, Терещенко] методом кусочно-планарной аппроксимации на сетке 21х21 узел. Такая крупная (по сравнению с реальными задачами томографии ионосферы) сетка была взята нами для того, чтобы проверить выдвигаемую методику с помощью точного метода решения (5), не прибегая к приближенным итерационным алгоритмам. Поскольку число арифметических операций, требующихся компьютеру для завершения точного метода, порядка (MK) 6 , то для расчета на сетке 100х100 узлов необходимо время в 15,000 раз большее. В скобках заметим, что на решение системы (5) методом Гаусса с выбором главного элемента [Рябенький] Pentium–200 затрачивает около минуты, поэтому на практике следует, конечно, использовать итерационные методы. На наш взгляд, наилучшим может оказаться, например, метод сопряженных градиентов, поскольку в его рамках можно организовать эффективный спуск по параметру a.
Задав реалистичный профиль N0(z) в качестве априорной информации, проиллюстрируем предлагаемый подход точным(!) решением регуляризованной задачи. Для того, чтобы определить оптимальный параметр регуляризации, проведем серию расчетов с различными a, контролируя невязку. На рис.2 изображена зависимость нормы невязки
(6)
решения регуляризованной задачи от параметра регуляризации a. На том же рисунке изображена суммарная ошибка (т.е. норма отклонения полученного решения от модельного распределения).
Для реконструкции можно использовать значение, соответствующее глобальному минимуму зависимости e(a) ( рис.3 ), либо применить т.н. принцип невязки [Гласко], который требует выбора a, с которым невязка приблизительно равна сумме погрешностей измерений (т.е. заданий правой части) и аппроксимации (результат реконструкции изображен на рис.4 ). Из сравнения рисунков видно, что наилучшее совпадение с моделью демонстрирует последний рисунок, соответствующий более оптимальному a.
В заключение заметим, что применяемые на практике методы реконструкции неявно также используют элементы регуляризации (для метода ART3 параметрами регуляризации являются число итераций и индексы релаксации, графики, подобные рис.3 , приведены в [Kunitsyn et.al.]). Однако такую регуляризацию нельзя считать доказательством устойчивости решения.
Конечно, разработанные к настоящему времени итерационные алгоритмы реконструкции отлично справляются с приведённой моделью ионосферы. Дальнейшие компьютерные эксперименты и сопоставление с результатами других алгоритмов, видимо, позволит выделить круг задач, в которых метод регуляризации может быть наиболее эффективным по сравнению с остальными.

Литература:
  1. Воеводин В.В. Вычислительные основы линейной алгебры. М., Наука, 1977.
  2. Гласко В.Б. Обратные задачи математической физики. Изд. МГУ.
  3. Данилкин Н.П., Тольский К.Л. Достоверность ионосферной томографии. Геомагнетизм и Аэрономия (в печати).
  4. Калиткин Н.Н. Численные методы. М., Наука, 1978.
  5. Куницын В.Е., Терещенко Е.Д. Томография ионосферы. М.,1991.
  6. Рябенький В.С. Введение в вычислительную математику. М., Наука, 1994.
  7. Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. М., Наука, 1974.
  8. Ценсор Я. ТИИЭР, т.71 (1983) №3, с.148.
  9. Kunitsyn V.E., et.al. Methods and Algorithms of ray tomography for ionospheric research. Ann. Geophys. 13(1995).